How exactly do you compute the Fast Fourier Transform?

Question

I've been reading a lot about Fast Fourier Transform and am trying to understand the low-level aspect of it. Unfortunately, Google and Wikipedia are not helping much at all.. and I have like 5 different algorithm books open that aren't helping much either.

I'm trying to find the FFT of something simple like a vector [1,0,0,0]. Sure I could just plug it into Matlab but that won't help me understand what's going on underneath. Also, when I say I want to find the FFT of a vector, is that the same as saying I want to find the DFT of a vector just with a more efficient algorithm?

Answer 1

You're right, "the" Fast Fourier transform is just a name for any algorithm that computes the discrete Fourier transform in O(n log n) time, and there are several such algorithms.

Here's the simplest explanation of the DFT and FFT as I think of them, and also examples for small N, which may help. (Note that there are alternative interpretations, and other algorithms.)

Discrete Fourier transform

Given N numbers f₀, f₁, f₂, …, f_N-1, the DFT gives a different set of N numbers.

Specifically: Let ω be a primitive Nth root of 1 (either in the complex numbers or in some finite field), which means that ω^N=1 but no smaller power is 1. You can think of the f_k's as the coefficients of a polynomial P(x) = ∑f_kx^k. The N new numbers F₀, F₁, …, F_N-1 that the DFT gives are the results of evaluating the polynomial at powers of ω. That is, for each n from 0 to N-1, the new number F_n is P(ωⁿ) = ∑_0≤k≤N-1 f_kω^nk.

[The reason for choosing ω is that the inverse DFT has a nice form, very similar to the DFT itself.]

Note that finding these F's naively takes O(N²) operations. But we can exploit the special structure that comes from the ω's we chose, and that allows us to do it in O(N log N). Any such algorithm is called the fast Fourier transform.

Fast Fourier Transform

So here's one way of doing the FFT. I'll replace N with 2N to simplify notation. We have f₀, f₁, f₂, …, f_2N-1, and we want to compute P(ω⁰), P(ω¹), … P(ω^2N-1) where we can write

P(x) = Q(x) + ω^NR(x) with

Q(x) = f₀ + f₁x + … + f_N-1x^N-1

R(x) = f_N + f_N+1x + … + f_2N-1x^2N-1

Now here's the beauty of the thing. Observe that the value at ω^k+N is very simply related to the value at ω^k:
P(ω^k+N) = ω^N(Q(ω^k) + ω^NR(ω^k)) = R(ω^k) + ω^NQ(ω^k). So the evaluations of Q and R at ω⁰ to ω^N-1 are enough.

This means that your original problem — of evaluating the 2N-term polynomial P at 2N points ω⁰ to ω^2N-1 — has been reduced to the two problems of evaluating the N-term polynomials Q and R at the N points ω⁰ to ω^N-1. So the running time T(2N) = 2T(N) + O(N) and all that, which gives T(N) = O(N log N).

Examples of DFT

Note that other definitions put factors of 1/N or 1/√N.

For N=2, ω=-1, and the Fourier transform of (a,b) is (a+b, a-b).

For N=3, ω is the complex cube root of 1, and the Fourier transform of (a,b,c) is (a+b+c, a+bω+cω², a+bω²+cω). (Since ω⁴=ω.)

For N=4 and ω=i, and the Fourier transform of (a,b,c,d) is (a+b+c+d, a+bi-c-di, a-b+c-d, a-bi-c+di). In particular, the example in your question: the DFT on (1,0,0,0) gives (1,1,1,1), not very illuminating perhaps.