I am trying to implement Matlab code which performs FIR filtering with complex bandpass filter using Overlap Save method and perform decimation as well by using different length forward and inverse FFT as mentioned in this article:
https://www.iro.umontreal.ca/~mignotte/IFT3205/Documents/TipsAndTricks/MultibandFilterbank.pdf
Below is the code snip I have written:
clear all;close all;clc
fs=1e6;
sim_time=20e-3;
Nfft=256;
t=0:1/fs:sim_time-1/fs;
freq = fs/Nfft*(0:Nfft-1);
%create composite signal
fsig1= 50e3;
fsig2= 100e3;
fsig3= 300e3;
fsig4= 872e3;
sig1=2*cos(2*pi*(fsig1-5e3)*t)+4*cos(2*pi*fsig1*t)+8*cos(2*pi*(fsig1+5e3)*t) + j*(2*sin(2*pi*(fsig1-5e3)*t)+4*sin(2*pi*fsig1*t)+8*sin(2*pi*(fsig1+5e3)*t));
sig2=8*cos(2*pi*(fsig2-5e3)*t)+4*cos(2*pi*fsig2*t)+2*cos(2*pi*(fsig2+5e3)*t) + j*(8*sin(2*pi*(fsig2-5e3)*t)+4*sin(2*pi*fsig2*t)+2*sin(2*pi*(fsig2+5e3)*t));
sig3=2*cos(2*pi*(fsig3-5e3)*t)+8*cos(2*pi*fsig3*t)+4*cos(2*pi*(fsig3+5e3)*t) + j*(2*sin(2*pi*(fsig3-5e3)*t)+8*sin(2*pi*fsig3*t)+4*sin(2*pi*(fsig3+5e3)*t));
sig4=10*cos(2*pi*fsig4*t)+j*10*sin(2*pi*fsig4*t);
sig=sig1+sig2+sig3+sig4;
% Create Filter
bpFilt1 = designfilt('lowpassfir', 'FilterOrder', 40, ...
'PassbandFrequency', 100e3, 'StopbandFrequency', 150e3,...
'SampleRate', fs);
%filter impulse response
h1= bpFilt1.Coefficients;
%create complex band pass by multiplying with complex exponential
fc=800e3;
cmp_exp=exp(j*2*pi*fc*t);
h1_c=h1.*cmp_exp(1:length(h1));
%take FFT of filter impulse response
H1=fft(h1,Nfft)/Nfft;% low pass prototype
H1c=fft(h1_c,Nfft)/Nfft; %complex bandpass filter
%implement overlap save
P=41;%filter length
for i=1:10 % This code snippet compute the OSA filter with decimation by 4
if(i==1)
sig_win=sig(1:Nfft);
index_fft_in=Nfft;
else
%sig_win=sig((i-1)*256-(P-1):(i-1)*256-(P-1)+Nfft-1); %overlapped input
sig_win=sig(index_fft_in-(P-1):index_fft_in-(P-1)+Nfft-1); %overlapped input
index_fft_in=index_fft_in-(P-1)+Nfft-1;
end
X_sig=fft(sig_win,Nfft)/Nfft;
Y_sig=H1c.*X_sig;%multiply with complex bandpass filter
y_sigd(i,:)=ifft(Y_sig((3*Nfft/4)+1:4*Nfft/4),Nfft/4)*(Nfft/2)^2;%extract complete for now to test filtering
if(i==1)
y_osa=y_sigd(i,((P-1)/4)+1:end);%discard the first (P-1) points
index=length(y_osa);
else
y_osa(index:index+(Nfft/4)-(P-1)/4-1)=y_sigd(i,(P-1)/4+1:end);
index=length(y_osa);
end
end
% figure;plot(freq/1e3,abs(H1),'-*r')
% hold on;plot(freq/1e3,abs(H1c),'-ob')%plot filter responses
y1 = filter(h1_c,1,sig(1:10*Nfft)); % apply filter 1 %reference for comparison
%compare decimation in FD vs TD
y1_d = downsample(y1,4);
figure;
plot(real(y_osa),'-ob')
hold on;plot(real(y1_d),'-*g')
I am also doing the same filtering in time domain along-with downsampling to ensure that the time domain and frequency domain implementations are equivalent however there is some glitch in the output when compared as shown in image below:
Can anyone point out the mistake in the code or the method and suggest solution?
Edit: The first block of the result appears to match the reference time-domain filtered signal but some issue occurs in the subsequent blocks.
**Also, the result is fine when the same length of IFFT is performed (without any decimation) and the result matches exactly that of the time-domain reference signal. The issue is only when smaller IFFT size is used for required bins.