Hi, I am creating a Matlab model of a multichannel system where the input signal
is sampled below the Nyquist rate in each channel, but is still recoverable
because each channel has a uniform phase offset. In this case there are M=3
channels that sample at half the Nyquist Rate (L=2). The channels have uniform
phase offsets of [0 Tsamp/3 2*Tsamp/3]. When corrected and interleaved, the
output should be a perfect reconstruction of the input.

This is fairly simple conceptually, but I am having trouble implementing the
reconstruction filters in the frequency domain. Does anyone have any ideas?

===============================================

clear;

clc;

M = 3; %M=Sampling Channels

L = 2; %L=Undersampling Ratio for each Channel

%Define Nyquist Frequency Bandlimit for Input Signals

fn = 21;

Tn = 1/fn;

%Define Sampling Offsets

Tsamp = L*Tn;

fsamp = 1/Tsamp;

%Phase Offset Between Channels

tau_m = [0 Tsamp/3 2*Tsamp/3];

%Define Sampling Grids

tmin = 0;

tmax = 100*1/0.5;

t0 = (tmin+tau_m(1)):Tsamp:tmax;

t1 = (tmin+tau_m(2)):Tsamp:tmax;

t2 = (tmin+tau_m(3)):Tsamp:tmax;

t0 = t0(1:(length(t0)-1));

%Sample Input Signal

x0 = cos(2*pi*10*t0);

x1 = cos(2*pi*10*t1);

x2 = cos(2*pi*10*t2);

x0 = upsample(x0,L);

x1 = upsample(x1,L);

x2 = upsample(x2,L);

%FFT of Input Signal

X0 = fft(x0)/length(x0);

X1 = fft(x1)/length(x1);

X2 = fft(x2)/length(x2);

%Define Frequency Grid for Sampling Reconstruction Filters

omega_mat0 = linspace(-pi,pi-2*pi/length(X0),length(X0));

omega_mat1 = linspace(-pi,pi-2*pi/length(X1),length(X1));

omega_mat2 = linspace(-pi,pi-2*pi/length(X2),length(X2));

%Reconstruction Filtering (Frequency Domain)

OUT0 = X0.*fftshift(2/3); %No Shift

OUT1 = X1.*fftshift(2/3*exp(j.*omega_mat1*2/3)); %Phase Correction for
Interleaving

OUT2 = X2.*fftshift(2/3*exp(j.*omega_mat2*4/3)); %Phase Correction for
Interleaving

OUT = OUT0+OUT1+OUT2;

figure

subplot(2,1,1)

stem(omega_mat2,fftshift(real(OUT)))

title('FFT OUT Total After Filters');

xlabel('Frequency (rad/s)');

ylabel('Real');

subplot(2,1,2)

stem(omega_mat2,fftshift(imag(OUT)),'r')

title('FFT OUT Total After Filters');

xlabel('Frequency (rad/s)');

ylabel('Imag');