# Is it numerically more stable to implement filtering as multiplication or convolution?

I'm writing a program to filter a 20,000-sample signal with a fifth-order Butterworth filter offline. I've got access to an FFT implementation. There seems to be two alternatives for implementing the filtering:

• convolving the signal with the impulse response in the time domain, or
• multiplying the signal with the impulse response in the frequency domain and inverse-transforming the result

These methods would be identical in the theoretical FT case. Doing it in real life with the DFT, though, I suppose things are different. Is one of the methods numerically more stable? Are there any other problems that I should be aware of? The number of calculations is not important.

• The FFT method will be much faster to calculate for most signal lengths. Only short lengths are faster with time-domain convolution. – endolith Nov 6 '11 at 22:59

With convolution, you are not going to run into any stability problems, because there is no recursive filtering, so you are not going to accumulate any errors. In other words, the system is all zeros, no poles. I have heard anecdotally, but not verified for myself, that FFT-based convolution has lower error than time-domain convolution, simply because it has O(n log n) arithmetic operations rather than O(n^2).

Typically, as far as I am aware, Butterworth filters are implemented as recursive (IIR) filters, so that is a different topic. IIR filters have poles as well as zeros, so there can be stability issues in practice. Also, for IIR filters, FFT-based methods are not an option, but on the upside, IIR filters tend to be very low order.

As far as stability issues with IIR filters, they tend to have problems at higher orders -- I'll just throw out a number and say roughly 6th order is pushing it. Instead, they are typically implemented as cascaded biquads (2nd order filter sections). For your 5th order filter, write it as a z-domain transfer function (it'll be a 5th degree rational function), and then factor it into its 5 poles and 5 zeros. Collect the complex conjugates, and you'll have two biquads and one first-order filter. In general, stability problems tend to crop up as the poles get closer to the unit circle.

There can also be issues with noise and limit cycles in IIR filters, so there are different filter topologies (i.e. direct form I, direct form II) that have different numeric properties, but I wouldn't overthink this point -- just use double-precision and it'll almost certainly be good enough.

• What do you mean with it only working for FIR filters? I assumed that IIR filters would have to be sampled somehow anyway. Are IIR filters usually implemented in the time domain to avoid this? – Andreas Nov 7 '11 at 5:27
• As far as I am aware, IIR filters are always implemented in the time domain. An IIR filter (here, for example, a second-order filter or "biquad") is defined by a difference equation like y(n) = b0 * x(n) + b1 * x(n-1) + b2 * x(n-2) - a1 * y(n-1) - a2 * y(n-2). Note that this is a combination of the previous input samples (the x values), and the previous output samples (the y values). An FIR filter only depends on the past inputs, so it admits an efficient frequency domain implementation. An IIR filter does not, but is very efficient anyway because IIR filters tend to be much lower order. – schnarf Dec 7 '11 at 22:03
• The reason IIR filters tend to be much lower order is that the poles (the feedback of previous output samples) lets the filter get much more steep with very few coefficients as compared to an FIR filter. When I say much lower order, a typical IIR filter might be second-order (5 coefficients), while a typical FIR filter for the same task might have thousands of coefficients. – schnarf Dec 7 '11 at 22:09

In nearly all cases your best choice is neither convolution nor FFT but simply applying the IIR filter directly (using e.g. the sosfilt() function). This will be vastly more efficient in terms of CPU and memory consumption.

Whether it makes a numerical difference depends on the specific filter. The only case where some difference may creep in are if the poles are very, very close to the unit circle. Even there a few tricks that can help. DO NOT USE transfer function representation and filter() but use poles and zeroes with sosfilt(). Here is an example for the difference.

n = 2^16;  % filter length
fs = 44100; % sample rate
x = zeros(n,1); x(1) = 1;
f0 = 15; % cutoff frequency in Hz
% design with poles and zeroes
[z,p,k] = butter(5,f0*2/fs);
clf
plot(sosfilt(zp2sos(z,p,k),x));
% design with transfer function
[b,a] = butter(5,f0*2/fs);
hold on
plot(filter(b,a,x),'k');


filter() goes bad at a cutoff of about 15Hz @44.1kHz. For sosfilt() the cutoff can be well below 1/100 of Hz @44.1kHz without any problems.

IF you have stability problems the FFT doesn't help much either. Since your filter is an IIR filter, the impulse response is infinite and would have to be truncated first. At these very low frequency the impulse response gets so long that the FFT becomes impractical as well.

For example if you want a cutoff of 1/100 Hz @ 44.1 kHz and want a dynamic range in the impulse response of 100 dB, you need roughly 25 million samples !!! That's almost 10 minutes at 44.1 kHz and many, many times longer than your original signal

• This doesn't really answer the question on numerical problems, but I wasn't aware of the problems with filter -- thanks! My high-pass cutoff is 0.5 Hz @ 250 Hz. What's the reason for the problems with filter? I'm writing the implementation myself. – Andreas Nov 8 '11 at 7:52

Why do you suppose things will be different? The theoretical concepts should translate to practical applications, with the only difference being that of floating point issues, which we can't escape. You can easily verify with a simple example in MATLAB:

x=randn(5,1);
y=randn(5,1);
X=fft(x,length(x)+length(y)-1);
Y=fft(y,length(x)+length(y)-1);

z1=conv(x,y);z2=ifft(X.*Y);
z1-z2

ans =

1.0e-15 *

-0.4441
-0.6661
0
-0.2220
0.8882
-0.2220
0
-0.4441
0.8882


As you can see, the errors are of the order of machine precision. There shouldn't be any reason to not use the FFT method. As Endolith mentioned, it is much more common to use the FFT approach to filter/convolve/cross-correlate, etc. and a lot faster except for very small samples (like in this example). Not that time-domain processing is never done... it all boils down to the application, needs and constraints.

• I think the original question was drilling right into the floating-point issues inherent in FFT-based filtering versus direct implementation of the filter in the time domain. This can be a real concern for fixed-point signal processing, if you have really long filters, or if you have a bad FFT implementation, for example. You definitely won't see any effects for sequences of length 5 in double precision floating point. – Jason R Nov 7 '11 at 15:14
• @JasonR The errors are still of machine precision if you extend the length of the sequences to 1e6 in the example above. The errors that you mention crop up primarily due to poor filter design or bad FFT implementation. If those are fine, I don't see why convolution in time-domain should give a different answer from that in the frequency domain. – Lorem Ipsum Nov 7 '11 at 15:41
• It shouldn't give a different answer based on which domain you do the computations in. My only point is that the actual mathematical operations differ greatly between the two approaches, so depending upon the implementations of each that you have available, it's possible that you could see tangible differences. For double precision, assuming you have good implementations, I wouldn't expect any difference except in extreme cases. – Jason R Nov 7 '11 at 15:45