Non-power-of-2 FFT's?

If I have a program that can compute FFT's for sizes that are powers of 2, how can I use it to compute FFT's for other sizes?

I have read that I can supposedly zero-pad the original points, but I'm lost as to how this gives me the correct answer. If I transform a function with 5 points, I expect to get back a 5-point DFT, and I don't understand how I can extract that from the 8-point DFT (which was calculated with zero-padding).

Solution 1:

Sure, you can use a radix-2 FFT to compute FFTs for lengths not a power of 2 (but it is not as efficient as using methods specifically tailored to the factors of the sequence length). However, it's not as easy as merely padding the original array. The right way of going about it, due to Rabiner, Schafer, and Rader, is termed the "chirp-z transform". (This is a slightly condensed version of the paper previously linked to.)

To motivate the chirp-z idea, consider the DFT

$$X_k=\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{k n},\qquad k = 0,\dots,N-1$$

The key step, attributed to Bluestein, considers the algebraic identity

$$kn=\frac12(k^2+n^2-(k-n)^2)$$

Substituting into the DFT and expanding out the powers yields

$$\begin{align*}&\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{\frac12(k^2+n^2-(k-n)^2)}\\=&\left(\left(e^{-2\pi i/N}\right)^{k^2/2}\right)\sum_{n=0}^{N-1} \left(x_n \left(e^{-2\pi i/N}\right)^{n^2/2}\left(e^{-2\pi i/N}\right)^{-\frac{(k-n)^2}2}\right)\end{align*}$$

and we then recognize that the summation expression is precisely in the form of a convolution. The factor $\left(e^{-2\pi i/N}\right)^{k^2/2}$ is what is termed as a "chirp"; hence the name "chirp-z transform".

Chirp-z thus consists of three stages:

taking the Hadamard (componentwise) product of the original sequence with the chirp
convolving with the reciprocal chirp (which of course is equivalent to the inverse FFT of the product of the FFTs of the Hadamard product and the reciprocal chirp)
taking the Hadamard product of that result with a chirp. We see that three or so FFTs are needed (barring the possibility of an exploitable symmetry being present).

MATLAB's Signal Processing toolbox implements this algorithm as the function czt(). See the code in the corresponding M-file for implementation details.

The Code

I (the original poster) ended up making a Python version using SciPy, so I thought I'd edit this answer to include the code:

from scipy import *

def czt(x, m=None, w=None, a=None):
    # Translated from GNU Octave's czt.m

    n = len(x)
    if m is None: m = n
    if w is None: w = exp(-2j * pi / m)
    if a is None: a = 1

    chirp = w ** (arange(1 - n, max(m, n)) ** 2 / 2.0)
    N2 = int(2 ** ceil(log2(m + n - 1)))  # next power of 2
    xp = append(x * a ** -arange(n) * chirp[n - 1 : n + n - 1], zeros(N2 - n))
    ichirpp = append(1 / chirp[: m + n - 1], zeros(N2 - (m + n - 1)))
    r = ifft(fft(xp) * fft(ichirpp))
    return r[n - 1 : m + n - 1] * chirp[n - 1 : m + n - 1]

@vectorize  # Rounds complex numbers
def cround(z, d=None): return round(z.real, d) + 1j * round(z.imag, d)

arr = [1.0, 2.0, 3.0]

print cround(czt(arr), 4)       # [ 6.0+0.j    -1.5+0.866j -1.5-0.866j]
print cround(fft(arr), 4)       # [ 6.0+0.j    -1.5+0.866j -1.5-0.866j]

Solution 2:

If you want to code a DFT of size $N=5$ you have several options:

8-point FFT, after zero-padding
Evaluate the DFT directly
Special FFT algorithms (eg: Rader)

The zero-padding option is popular, and it's exact (in two senses: the inverse gives you back the original zero-padding sequence; and both the 8-point transform and the 5-point transform correspond to the same underlying continuous DTFT, only sampled at different frequencies). But you can't (directly, efficiently) obtain from it the 5-point DTF. So, of you really want those 5 values, this is not the alternative.

The naive alternative, to compute the 5-point DFT directly, evaluating the formula, is -of course- not a FFT (not 'fast') but for such a small size the difference might not be important.

The other alternative (a true 5-point FFT) is the most efficient but also the more complex to understand and implement. And you can't use a elementary power-of-two FFT program to implement this.