0
$\begingroup$

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:

Frequency Domain vs Time Domain Comparison

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.

Improve this question
$\endgroup$
1
$\begingroup$

Attempt 2

I am not quite sure what is the interpretation of your filtering, in special, the IFFT part, you are taking a slice of the spectrum shifted by 3*Nfft/4, I think it should be (circularly) shifted by 7*Nfft/8.

  % your method
  xa = ifft(X((3*Nfft/4)+1:4*Nfft/4),Nfft/4) / 4;
  % what I would do
  xb = ifft(fftshift(fftshift(X)(3*Nfft/8+1:5*Nfft/8))) / 4;

Notice that you factor is (Nfft/2)^2 because you multiplied the spectra by Nfft.

To make sure I am not missing any constant I plotted the histogram of the correlations between the fft(trimmed(signal)) and decimated(ifft(signal)) for your method and the modified formula (correlation of 1 is exact up to an affine transformation).

corelation analysis

The complete code I used to check is this

cora = []
corb = [];
for ii = 1:10000
  x = zeros(Nfft, 1);
  x(1:4:end) = randn(Nfft/4, 2)*[1;1j];
  x = cumsum(x);
  X = fft(x);
  % your method
  xa = ifft(X((3*Nfft/4)+1:4*Nfft/4),Nfft/4) / 4;
  % what I would do
  xb = ifft(fftshift(fftshift(X)(3*Nfft/8+1:5*Nfft/8))) / 4;

  % make sure this is not for a constant factor
  cora(ii) = abs(corr(xa, x(1:4:end))); % 0.66
  corb(ii) = abs(corr(xb, x(1:4:end))); % 0.98
end
figure
subplot(211),hist(cora, 25),xlim([0, 1]);
title('correlation histogram for your method'), xlabel('correlation')
subplot(212),hist(corb),xlim([0,1])
title('correlation histogram for the proposed method'), xlabel('correlation')

Working your example.

A few modifications that makes your frequency and time domain filtering consistent

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;
sig=randn(size(sig)); # white noise for testing
P=41;%filter length
h1 = exp(-((1:P)/10).^2); % (Octave) without designfilt

%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);% low pass prototype
H1c=fft(h1_c,Nfft); %complex bandpass filter
%implement overlap save 
wLen = Nfft - (P - 1);
BLOCKS = 10;
y_osa = zeros(1, wLen * BLOCKS);
for i=1:BLOCKS % This code snippet compute the OSA filter with decimation by 4
    sig_win = sig(wLen*(i-1)+(1:Nfft));
    X_sig=fft(sig_win,Nfft);
    Y_sig=H1c.*X_sig;%multiply with complex bandpass filter
    y_sigd = ifft(Y_sig);
    y_osa(wLen*(i-1)+1:wLen*i)=y_sigd(P:Nfft);
end

y1 = conv(h1_c, sig(1:(BLOCKS+1)*wLen));
y1 = filter(h1_c, 1, sig);
y1 = y1(P:wLen*BLOCKS+P-1);

figure
hold on
plot(real(y1))
plot(real(y_osa))
```
Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you for your response. My understanding of the overlap save method is that firstly, no zero padding of the signal is required, therefore passing sig_win with Nfft points. My idea was to pass the first chunk as is and then use the index pointer to overlap the P-1 points in the next chunk and so on. Also for the output y_sig the result is to be decimated hence using IFFT of N/4 points. $\endgroup$
    – malik12
    Mar 4 '21 at 9:57
  • $\begingroup$ Sorry I mixed overlap-save and overlap-add here. Now I tested before give you an answer :) check if that solves. $\endgroup$
    – Bob
    Mar 4 '21 at 13:50
  • $\begingroup$ No issue and thanks again for your elaborated response. As far as the interpretation of the filtering is concerned, the input is basically a complex I/Q input hence valid frequency range is 0 to Fs and currently I am interested in 750Khz to 1000Khz frequencies therefore taking only the last 64 bins (3*Nfft/4+1: Nfft) for IFFT to perform decimation as well. As this approach or interpretation incorrect? $\endgroup$
    – malik12
    Mar 4 '21 at 14:04
  • $\begingroup$ So you have a problem, your sampling frequency is 1000KHz, you can only recover frequencies under the Nyquist frequency (fs/2 = 500KHz). $\endgroup$
    – Bob
    Mar 4 '21 at 14:13
  • $\begingroup$ But I am using IQ sampling therefore the non-aliased range is 1000Khz that is the same as sampling frequency is it not? $\endgroup$
    – malik12
    Mar 4 '21 at 14:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.