This method uses a convolution to evaluate an FFT. Take a look at Annex C if your not familiar with convolution. As FFT's go, this method is fairly slow, but has the advantage that it will work with any length of sequence, including primes. The method works by expressing the DFT as a convolution, and using 1 FFT and 1 IFFT to perform this convolution (assuming certain constants are available as look up tables). These forward and inverse transforms can be performed with any FFT/IFFT routines you happen to have available.
The method described here is a simplified version of a more general transform called the 'Chirp Z Transform' (CZT). Since this document is concerned with the DFT, which is a special case of the CZT, the full detail of the CZT is not presented.
Recall the definition of the DFT of a sequence of N samples:
The complex exponential factors can by rearranged:
Using this form in the DFT yields:
Now, all our other FFT definitions had the scaling factor dropped. Since what we're going to develop here is another FFT, we shall also drop the scaling factor in the above expression, for the sake of consistency:
The summation looks like a convolution (it is), but we need to be careful about how this expression is interpreted. We can't regard this as a circular convolution with period N, because the complex exponentials do not in general have period N (but they do if N is even). In any case, (for even N) a circular convolution would need an N point FFT. Since the whole purpose of this exercise is to derive an N point FFT, that's not much help. Fortunately, the fact that we're only interested in values for k in the range k=0..N-1 enables us to evaluate the summation as the convolution of two finite length (aperiodic) sequences. Annex C describes how to do this with any convenient FFT (such as a Radix 2) which is sufficiently long.
The FFT becomes:
where C(k) is the convolution..
This is the aperiodic convolution of 2 finite length sequences of lengths N (x) and 2N-1 (y). The length of the convolved result will be N+(2N-1)-1=3N-2. To be precise, C(k) will be non-zero for: k=-(N-1)..2(N-1). In Annex C it is shown that this convolution can be 'circularised' provided the repetition period was at least as great as the length of the convolved result (so that no overlap occurred). However, this was based on the assumption that we want the whole convolved result. In this case we don't, we only want the C(k) for k=0..N-1. So some overlap is acceptable, provided the overlapping images (aliases) don't intrude on this region. This is important because we don't want to use an FFT which is any larger than necessary.
Suppose we chose an NF point FFT. The convolved result will have periodicity NF. We require that the boundaries of overlapped images (aliases) don't intrude on the region k=0..N-1.
In other words, for the lower bound of the upper image, we require:
Similarly, for the upper bound of the lower image, we require:
Both these constraints are effectively the same constraint on NF. But we also have an additional constraint, that NF is sufficient to accomodate the sequence y (which has length 2N-1). Clearly, this is automatically satisfied assuming the above constraint is also. So the length of FFT we must use is constrained by:
Well that's the hard part done. We know what size FFT we need, all that remains to be done is to outline an algorithm based on this.
To evaluate an N point DFT by convolution, we chose a suitable size NF of FFT, such that:
This will usually be an integral power of 2, so that a Radix 2 algorithm can be used.
Generate an array y[n], n=0..NF-1:
Evaluate the FFT of y[n]:
(If more than one N point FFT is required, the array Y is a constant which can be re-used.)
Generate an array x[n], n=0..NF-1:
Evaluate the FFT of x[n]:
Multiply X and Y element by element (NF elements total) and take the IFFT:
Finally, extract the FFT we want from the first N elements of C (which has NF elements):
The scaling factor (1/NF) has been included explicitly in the above expression, on the assumption that the (so called) IFFT is unscaled. The net result is an unscaled (just like all our other FFT's) N point FFT algorithm.
If we assume that we want to use Radix 2 FFT's in this algorithm, then a generalisation of the best and worst cases are N=2p-1 (best case) and N=2p+1 (worst case). In the best case we can use NF=2p+1, so we need an FFT and an IFFT, both about twice the length of our sequence. So the very best performance of this algorithm is at least 4 times slower than a straight 2p point FFT. In the worst case we need to use NF=2p+2, so were looking at a factor of at least 8 difference here.
The worst case situation could be mitigated by chosing a suitable mixed radix algorithm (3 point transforms aren't that hard). However, there's no avoiding the fact that an FFT and IFFT, each of about twice the length of the required FFT, will be needed in the end.
©1999 - Engineering Productivity Tools Ltd.