Implementation of Overlap Save Filtering with Decimation
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))
```