# Convolution to do FIR with initial and final conditions!

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!

• The difference between the outputs of a FIR filter with and without initial states would only occur at the start of the output (same length as b, which I will call N). So you could add the initial states to the beginning of x and crop of the first N samples of y, or do conv normally and replace the first N values of y with the results of filter. – fibonatic Aug 13 '16 at 14:54
• @fibonatic, thanks I got it working an hour after I posted this question, exactly the way you described. – Ali Aug 13 '16 at 19:02

## 2 Answers

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)));

• I'm running Matlab 2012 on a RedHat server. It makes sense that conv beats filter in execution time for an FIR case but I didn't expect that big of a difference. – Ali Aug 13 '16 at 19:06
• After posting my question I tried writing my own C code and it barely beat conv function. Then I searched a bit more and came across fftfilt and it's amazingly fast. But then again I had the problem to preserve state from past runs, I did similar tricks as yours and it's working great. Thanks a lot. – Ali Aug 13 '16 at 19:13

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).

• as I said in my other comments, I did expect conv to perform faster but not an order of magnitude faster. Thanks for describing overlap-add method, most likely fftfilt function is implemented very much what you described. Regarding your comment about Matlab, for this particular question I have to use Matlab. I have been using Matlab and C since 1996, and Python for the past few years, for simulations and numerical analysis and I totally understand if someone prefers one over others. – Ali Aug 13 '16 at 19:23