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).
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:
MATLAB's Signal Processing toolbox implements this algorithm as the function czt(). See the code in the corresponding M-file for implementation details.
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]
arr = [1.0, 2.0, 3.0], and not arr = [1.0, 2.0, 3.0, 2.5, 3.5].
$\endgroup$
If you want to code a DFT of size $N=5$ you have several options:
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.
Are you aware of FFTW? It's been available for free for years; I have used it and can attest to its ease of use and speed.
