# Implementation of Overlap Save Filtering with Decimation

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)
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.

## 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).

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);
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')


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))
$$$$
`
• 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. Mar 4 '21 at 9:57
• Sorry I mixed overlap-save and overlap-add here. Now I tested before give you an answer :) check if that solves.
– Bob
Mar 4 '21 at 13:50
• 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? Mar 4 '21 at 14:04
• So you have a problem, your sampling frequency is 1000KHz, you can only recover frequencies under the Nyquist frequency (fs/2 = 500KHz).
– Bob
Mar 4 '21 at 14:13
• But I am using IQ sampling therefore the non-aliased range is 1000Khz that is the same as sampling frequency is it not? Mar 4 '21 at 14:19