It is well known that the output of a filter can be obtained by convolving the filters input with its impulse response. (Here we are assuming that the filter is a linear time invariant system.) It is also well known that the the Fourier Transform (discrete of otherwise) of the output is obtained by multiplying the Fourier Transform of the input by the Fourier Transform of the filters impulse response (the filters 'frequency response'.) There may also be a scaling factor to consider, depending on how the Fourier Transform was defined.
In this Annex we show how use the DFT (FFT) as an efficient convolution algorithm.
Given two discrete signals f and h which are defined for all integral n, the convolution of these signals is defined as:
In effect, this expression is simply the sum of shifted copies of h (centred on n'), which are weighted by the value f at n'. So we can think of h as the impulse response and f as the signal. However, this distinction is artificial and unnecessary. A simple substitution shows that convolution commutes:
If we substitute:
and note that the summation covers all n', and therefore all n'' (independent of n):
If both f and h are periodic, with period N, then their convolution (as defined above) will yield an infinite result (assuming the signals are non-zero). However, we can usefully re-define the convolution of such signals:
It is easy to see that g is also periodic, with period N. (Strictly, this only requires that h is periodic, but we also require that f is periodic to ensure circular convolution commutes.)
If we take the DFT of g:
The last step above is the first trick. See what we've done?
A little re-arrangement reveals:
The second trick is to notice that the inner summation is independent of n', despite the fact that n' appears in both factors. Both factors are periodic, and the sum includes every possible product, so it doesn't matter what order we use to sum the products:
So, the final result is:
It is worth noting that had we not used normal DFT scaling, but instead used the 'raw' (un-scaled) FFT, we would have obtained the result:
So time domain convolution is accomplished by a simple multiplication in the frequency domain. This property is one of the most useful features of the FFT. An N point circular convolution has algorithmic complexity O(N2), but using FFT's we can reduce this to O(N.logN). So large convolutions are often performed with FFT's. It is also simple to show (and not in the least surprising) that the converse is also true. Time domain multiplication is accomplished by frequency domain convolution:
If we return to the original (aperiodic) definition of convolution, we can observe that if both f and h are of finite length then so is their convolution g. Here, we assume that the sequences are zero other than in some defined range. Specifically, if f has length Nf and h has length Nh then g has length Ng =Nh +Nf -1. To see this, recall the definition of convolution:
The minimum and maximum values of n' which yield non-zero f(n') are 0 and Nf -1 respectively.
For n'=0, the minimum and maximum values of n which yield non-zero h(n-n') are 0 and Nh -1 respectively. Similarly, for n'=Nf -1, the minimum and maximum values of n which yield non-zero h(n-n') are Nf -1 and (Nf -1)+(Nh -1)=Nf +Nh -2 respectively. So..
Suppose that instead of being finite, the sequences f and h were both periodic, with period N. (This could be accomplished by padding the sequences with the required number of zeros between periods.) Circular convolution would result in a summed periodic repetition of g. Provided N>=Ng there would be no overlap between periods, so an Ng point 'snap shot' of the resulting summed periodic sequence will be the same as for the aperiodic convolution. This means that we can use an N point FFT to perform this convolution, where N>=Nf +Nh -1.
For convienience, we would probably want to use a Radix 2 FFT, so a good choice for performing convolutions by FFT (in a Finite Impule Response (FIR) filter for example) is to design a filter of Nh=2p+1 taps, and operate on block sizes of Nf=2p , the resulting minimum size of FFT is then N=Nf +Nh -1=2p+2p+1-1=2p+1. (The output from such a scheme would be overlapping blocks of 2p+1 points at 2p point spacing which need to be added.)
©1999 - Engineering Productivity Tools Ltd.