Matching outputs of FIR filters based on time domain convolution method and overlap-save method

Question

I'm trying to match(to the best extent possible) the outputs of FIR filter based on Frequency domain overlap-save method and time domain convolution method resp.

I was able to get most of the output to match except the tail as shown in the figure below. The overlap-save method always seems to lose the last half cycle for this particular signal.

It seems easier to deal with issues like filter delay and edge artifacts in time domain approach by pre-loading the filter with a part of the input stream and extracting a portion of the output signal, but is there any way to do the same in the frequency-domain/overlap-save approach for the full duration of filtered signal and also generalize for any signal length, cycles and sampling rate?

The Octave code for the two filters and the test script:

%% test filters

clc; clear all; close all;
pkg load signal

h = fir1(28, 6/24);

figure; stem(0:length(h)-1, h); title('LP Filter Impulse Response');

F = 48e3;
t = (0:1:2000)*(1/F);
f1 = 1e3;
f2 = 15e3;

signalIn = sin(2*pi*f1*t) + sin(2*pi*f2*t);
figure; plot(t, signalIn, 'b');

signalOutTDM =  FirFilterTimeD(signalIn, h);
signalOutFDM =  FirFilterFreqD(signalIn, h);

figure;
plot(1:length(signalOutTDM), signalOutTDM, 'r'); hold on;
plot(1:length(signalOutFDM), signalOutFDM, 'g');
title('Signal Out (index vs amp)'); legend('TimeDomainFIR','FreqDomainFIR');

%% Overlap-save method
function out =  FirFilterFreqD(timeDomain, filterTaps)
  
  idxOffset = 1;
  tdSize = length(timeDomain);
  numTaps = length(filterTaps);
  
  FFTlen = 2^(ceil(log2(numTaps)) + 5); % FFT length is typically set to a value much greater than numTaps 

  if FFTlen > tdSize
     FFTlen = 2^(floor(log2(tdSize))-2);  
  endif
    
  overlap = numTaps - 1;
  
  numZpad = overlap + floor(numTaps/2);
  if size(timeDomain,1) > size(timeDomain,2) 
     timeDomain = [zeros(numZpad,1); timeDomain; zeros(floor(numTaps/2),1)];
  else
     timeDomain = [zeros(1,numZpad) timeDomain zeros(1, floor(numTaps/2))];
     
  tdSizePadded = length(timeDomain);
  out = zeros(tdSizePadded,1);
  
  H = fft(filterTaps, FFTlen);
  
  pos = 0;
  stepLength = FFTlen - overlap;
  while (pos + FFTlen) <= tdSizePadded
        out_ifft   = ifft(fft(timeDomain(pos + (0:FFTlen-1) + idxOffset), FFTlen) .* H, FFTlen);
        out_select = out_ifft((numTaps-1:FFTlen-1) + idxOffset); % Take last stepLength values or discard first numTaps-1 values
        out(pos + (0:stepLength-1) + idxOffset) = out_select;       
        pos += stepLength;        
  endwhile

  out = out((overlap-1:tdSizePadded-numTaps-1) + idxOffset); 
  
end

%% Time Domain FIR filter
function out =  FirFilterTimeD(timeDomain, filterTaps)

idxOffset = 1;
tdSize = length(timeDomain);
out = zeros(tdSize,1);
idx = 0;
idxout = 0;
numTaps = length(filterTaps);
if mod(numTaps,2)==0
  delay = numTaps/2;
else
  delay = (numTaps-1)/2;
end
  
tdi = flip(timeDomain);

for i=0:tdSize-1
  idxout = tdSize - delay - i;
  if (idxout<0)
    idxout = tdSize + idxout;
  end
  
  out(idxout+idxOffset) = 0.0; 
  for x = 0:numTaps-1
    idx = mod((i+x),length(tdi));
    out(idxout + idxOffset) += tdi(idx + idxOffset) * filterTaps(x + idxOffset);
  end
  
end

end

Answer 1

The code below should solve your frequency domain problem but there is a problem with your time domain implementation too. For the frequency domain problem I have appended more zeros at the end.

%% test filters

clc; clear all; close all;
pkg load signal

h = fir1(28, 6/24);
Lgd=floor(length(h)/2);

figure; stem(0:length(h)-1, h); title('LP Filter Impulse Response');

F = 48e3;
t = (0:1:2000)*(1/F);
f1 = 1e3;
f2 = 15e3;

signalIn = sin(2*pi*f1*t) + sin(2*pi*f2*t);
figure; plot(t, signalIn, 'b');

signalOutTDM =  FirFilterTimeD(signalIn, h);
signalOutFDM =  FirFilterFreqD(signalIn, h);
signalOut =  conv(h,signalIn);
signalOut = signalOut(Lgd:end-(Lgd+1));

figure;
plot(1:length(signalOutTDM), signalOutTDM, 'r'); hold on;
plot(1:length(signalOutFDM), signalOutFDM, 'g-o');
plot(1:length(signalOut), signalOut, 'bo');
title('Signal Out (index vs amp)'); 
legend('TimeDomainFIR','FreqDomainFIR','Filter');

figure
plot(abs(signalOutFDM-signalOut.'))
title('Error sigmal: FDM implementation - Build-in filter funciton')
figure
plot(abs(signalOutTDM-signalOut.'))
title('Error signal: TDM implementation  - Build-in filter funciton')

%% Overlap-save method
function out =  FirFilterFreqD(timeDomain, filterTaps)
  
  idxOffset = 1;
  tdSize = length(timeDomain);
  numTaps = length(filterTaps);
  
  % FFTlen = 2^(ceil(log2(numTaps)) + 5); % FFT length is typically set to a value much greater than numTaps 
  FFTlen = 2^nextpow2(8*numTaps); % FFT length is typically set to a value much greater than numTaps 

  if FFTlen > tdSize
     FFTlen = 2^(floor(log2(tdSize))-2);  
  endif
    
  overlap = numTaps - 1;
  numZpad = overlap + floor(numTaps/2);
    if size(timeDomain,1) > size(timeDomain,2) 
     timeDomain = [zeros(numZpad,1); timeDomain; zeros(overlap+FFTlen,1)];
  else
     timeDomain = [zeros(1,numZpad) timeDomain zeros(1, overlap+FFTlen)];
  end
     
  tdSizePadded = length(timeDomain);
  out = zeros(tdSizePadded,1);
  
  H = fft(filterTaps, FFTlen);
  
  pos = 0;
  stepLength = FFTlen - overlap;
  i=1;
  while (pos + FFTlen) <= tdSizePadded
         out_ifft   = ifft(fft(timeDomain(pos + (0:FFTlen-1) + idxOffset), FFTlen) .* H, FFTlen);
        out_select = out_ifft((numTaps-1:FFTlen-1) + idxOffset); % Take last stepLength values or discard first numTaps-1 values
        out(pos + (0:stepLength-1) + idxOffset) = out_select;       
        pos += stepLength;        
  endwhile
 
  out = out(overlap:tdSizePadded-FFTlen-overlap-floor(numTaps/2) - 1); 
  
 end