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)
        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.
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);
  % 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))
```