0
$\begingroup$

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.

enter image description here

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
$\endgroup$
3
  • $\begingroup$ Please post code that actually runs. Matlab doesn't have "+=" operators or "endhwile" keyword $\endgroup$
    – Hilmar
    Aug 31, 2022 at 16:05
  • $\begingroup$ @Hilmar It is Octave code as mentioned in the description and the attached plot was generated by this script. Octave does support the += operator unlike Matlab. $\endgroup$
    – Naveen
    Aug 31, 2022 at 20:32
  • 1
    $\begingroup$ Oops, my bad. Sorry about this $\endgroup$
    – Hilmar
    Sep 1, 2022 at 0:00

1 Answer 1

1
+100
$\begingroup$

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
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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