Convolution to do FIR with initial and final conditions!

Question

I need to apply an FIR filter on chunks of very long signals. I can use MATLAB's filter function as follows to preserve the state of the filter when a new chunk of data is available:

[y, state] = filter(b, 1, x, state);

where b includes the coefficients of my FIR filter. However, b is fairly long and even each chunk of data is long and it appears to me that conv performs much faster than filter for this case. The following is a simple benchmark to compare the execution time for these two functions:

xin = rand(1, 1e6);
b = rand(1, 20e3);
tic; y = filter(b, 1, xin); toc
Elapsed time is 9.533244 seconds.

tic; y = conv(b, xin); toc
Elapsed time is 0.634058 seconds.

Unfortunately I can't preserve the state of my filter using conv function. I already tried to overlap the chunks of data and mimic the "state" for the convolution but haven't been successful yet. All suggestions and advice are greatly appreciated!

Answer 1

I tried the same benchmark and got very different results. On my machine and using Matlab R2106b (pre release), I got the following:

xin = rand(1, 1e6);
b = rand(1, 20e3);
% FILTER
tic; y1 = filter(b, 1, xin); t = toc;
fprintf('FILTER: Elapsed Time = %6.3f s\n',t);
% CONV
tic; y2 = conv(b, xin); t = toc;
fprintf('CONV: Elapsed Time = %6.3f s\n',t);
% FFTFILT
tic; y3 = fftfilt(b,xin); t = toc;
fprintf('FFTFILT: Elapsed Time = %6.3f s\n',t);

FILTER: Elapsed Time =  5.631 s
CONV: Elapsed Time =  4.921 s
FFTFILT: Elapsed Time =  0.105 s

There are a few things going on here:

1. FILTER is a generic all purpose IIR filter. Even if you don't give it any poles, the code needs to be specifically written to take advantage of this. It probably isn't since there are different functions for handling FIR filters
2. CONV is brute force and hence it's fairly slow, although it doesn't have any IIR code in it. Also note that with conv() the result is longer as it includes the "ringing" after the input has stopped (and is assumed to be zero). The excess samples can be viewed as the "state" or you filter.
3. FFTFILT is by far the fastest since it basically implements the overlap add algorithm in the frequency domain.
4. You could still speed this up by a factor of two or so by taking advantage of the fact that your input is real valued. That's beyond the scope of this question.

All of the methods allow you to "keep state" across calls, you just have to write the code for it. Below is an example that would work with FFTFILT.

%% state keeping example
xin = rand(1, 1e6)';
b = rand(1, 20e3)';
y3 = fftfilt(b,xin);

% break it up in 10 blocks
y4 = zeros(size(y3));
nBlocks = 10;
blockSize = length(xin)/nBlocks;
nTaps = length(b);
state = zeros(nTaps,1);
x2 = reshape(xin,blockSize,nBlocks);
counter = 0;
% implement overlap save
for i= 1:nBlocks
  % grab an input block, prepend it with the state and filter
  tmp = fftfilt(b,[state; x2(:,i)]);
  % save the good samples
  y4(counter + (1:blockSize)) = tmp((nTaps+1):end);
  % update the state: save the last nTaps input samples
  state = x2((end-nTaps+1):end,i);
  % advance output counter
  counter = counter + blockSize;

end
fprintf('Error = %6.2f dB\n', 10*log10(sum((y4-y3).^2)/sum(y3.^2)));

Answer 2

Your benchmark is actually quite surprising! I would've expected filter to choose an appropriate implementation of convolution to filter the input data, but thinking about this, since filter needs to pass the state in a format that isn't dependent on the convolution implementation actually used, the whole process of generating the state might take longer than the actual filtering; note that this really is just speculation. (Assumption: filter with long input data uses FFT FIR implementation internally, but the state that filter returns are taps of a direct-form FIR implementation – you basically need to reverse quite a bit of the filtering you've done to get that)

I already tried to overlap the chunks of data and mimic the "state" for the convolution but haven't been successful yet.

What you probably would like to have is a well-implemented FFT FIR! Since b doesn't seem to change, this would be a very efficient way of implementing the filtering operation.

The basic idea would be described as overlap-add filter (wikipedia), and the rough idea is this:

1. Pad your tap vector $b$ with zeros to length $N$. Get the $N$-point discrete Fourier transform of that using the FFT. Call the result $H$.
2. Take segments of length $L$ of the input signal, and pad them to $N$, transform as above.
3. point-wise multiply, since multiplication in discrete freq domain is equivalent to circular convolution in time domain.
4. inverse DFT (via IFFT)
5. By zero-padding enough, we now that a $L$ long part of the output signal can be generated by adding up the $N-L$ last samples of the last iteration to the start of the samples of this iteration.
6. save the $N-L$ last output samples for next iteration

---

personal addition: If you already know what filter you want, and you're actually after speed and don't care as much that you might need to leave Matlab to do the filtering (I personally cherish every moment that I don't have to use Matlab), look at this explanatory blog post that I've written, which actually does a lot of of band pass filters in parallel – you'd only need one of these, and you could choose a FFT FIR block to do the work; I've been doing 160 megasamples per second with these, on a "normal" PC. This is, as I'll frankly admit, based on GNU Radio and the VOLK acceleration library, because that's what I prefer for streaming samples DSP (which you seem to be doing in Matlab).