More FFT algorithms

An FFT splits a DFT into small ones, one stage per factor of $N$. Guess first: how much does it save at $N=31$, a prime?

The factors of N set the stages

The FFT (14.1) split a DFT by 2: even-indexed samples in one half, odd-indexed in the other. That trick is not special to 2. If $N=N_1N_2$, write $n=N_2n_1+n_2$. The DFT becomes $N_2$ DFTs of length $N_1$, a twiddle on each output, then $N_1$ DFTs of length $N_2$. With $N_2=2$ this is 14.1’s split again, because $n_2=0$ picks the even samples and $n_2=1$ the odd ones.

Repeat the split on the small DFTs and the result is a chain of stages, one for each factor in $N=N_1N_2\cdots$. The size $N_i$ of a stage is its radix, and a length whose factors differ is called mixed radix. A stage of radix $N_i$ is $N/N_i$ small DFTs of $N_i^2$ multiplies each, so it costs $N\cdot N_i$, and the whole FFT costs

multiplies. This is Cooley and Tukey’s count. It also counts the multiplies by $\pm1$ that 14.1 skipped, so at $N=32$ it says 320 where 14.1’s $\tfrac N2\log_2N$ said 80.

The picture at the top of the page draws these stages. For $N=12=2\times2\times3$ there are three, and the work is $12\times(2+2+3)=84$, against 144 directly. For $N=32=2\times2\times2\times2\times2$ there are five stages of radix 2, and the work is 320 against 1024. $N=30=2\times3\times5$ is mixed radix, three stages, and 300 against 900: $N$ need not be a power of two.

And the guess: $N=31$ is prime, so there is no factor to split by. It stays one 31-point DFT, 961 multiplies, the same as direct. The figure further down rescues it.

When the clip ends, drag up or down on the plot, or use the arrow keys, to choose $N$ from 2 to 32; the caption names the factors and the work.

Try $N=16$, which gives four stages of radix 2: work 128 against 256. Then try $N=17$, one more sample and a prime: 289 against 289, no saving at all. A length with one large factor, like $N=62=2\times31$, saves little: 2046 against 3844. The saving comes from small factors.

A prime length, by three FFTs

A prime length has no factors to split by, but it can be rewritten as a convolution, and a convolution can be done with a longer FFT. The rewriting starts from an identity for the product in the exponent of $W_N^{kn}$:

Check it with $k=3$ and $n=5$: $\tfrac12(9+25-4)=15$, which is $3\cdot5$. Put it in the exponent of $W_N^{kn}$:

The DFT now has three steps. Multiply $x[n]$ by the chirp $e^{-j\pi n^2/N}$, a tone whose frequency rises steadily. Convolve the result with the opposite chirp, $e^{+j\pi n^2/N}$. Multiply each output by the chirp again. The convolution in the middle is the one from Circular vs linear convolution (13.4): do it with FFTs, zero-padded to any power of two that is at least $2N-1$ so that it stays linear.

For $N=1009$, a prime, the first power of two at least $2\cdot1009-1=2017$ is 2048.

The count follows from 14.1. One 2048-point FFT is $\tfrac{2048}{2}\cdot11=11\,264$ complex multiplies, and the convolution needs three of them: 33 792. Add 2048 products of the two spectra and 2018 chirp multiplies, 1009 before and 1009 after, and the total is 37 858. Computing the DFT directly takes $1009^2=1\,018\,081$, about 27 times as many.

Two other ideas get one sentence each. Rader’s method also rewrites a prime-length DFT as a convolution, using a different re-indexing. Prime-factor (Good and Thomas) indexing splits $N=N_1N_2$ with no twiddles at all, when $N_1$ and $N_2$ share no divisor.

Two real FFTs for the price of one

Most signals are real, and the FFT is built for complex input. That wastes half of every transform: the imaginary part of the input is zero. So put one real signal in the real part and another in the imaginary part:

and take one FFT, $G[k]$. Because $x$ is real, its DFT has $X[N-k]=X[k]^*$ (The DFT, 13.2), and so does $y$‘s. The $*$ is the conjugate, the mirror image of Complex numbers for signals (3.3). Take the conjugate of $G$ at the mirrored bin:

Adding this to $G[k]=X[k]+jY[k]$ leaves $2X[k]$, and subtracting it leaves $2jY[k]$. So

with $G[N]$ meaning $G[0]$. One complex FFT carries two real DFTs, and the mirrored conjugate pulls them apart.

Before you watch, think of two voices on one phone line, one in each ear of a stereo pair. Which part of $G$ would tell you which voice is which?

One complex FFT of $g=x+jy$ gives $G=8,4,4,0,0,0,-4,4$.

Mirrored and conjugated, $G[N-k]^*=8,4,-4,0,0,0,4,4$. At $k=2$ and 6 the two disagree in sign. Read bin 2 yourself: $G[2]=4$ and $G[6]^*=-4$. Half their sum is 0, and half their difference divided by $2j$ is $4/j=-4j$.

Half the sum is $X$: 8 at $k=0$ (the 1), and 4 at $k=1$ and 7 (the cosine), while bins 2 and 6 cancel. Half the difference, divided by $2j$, is $Y$: $-4j$ at $k=2$ and $+4j$ at $k=6$ (the sine). Two real DFTs come from one FFT, and they match two separate ones.

Notice the bars at $k=2$ and $k=6$. They cancel in $X$ and double in $Y$, because the mirrored conjugate of $G$ keeps the sign of the real signal’s bins and flips the sign of the imaginary signal’s.

The same graph, run backwards

Decimation in time, from 14.1, splits the inputs into even and odd samples. Decimation in frequency splits the outputs. First form two sequences of length $N/2$ from the two halves of $x$:

The even-indexed outputs $X[0],X[2],\dots$ are the $N/2$-point DFT of $a$, and the odd-indexed ones are that of $b$. The reason is one line each. For even outputs $W_N^{2m(n+N/2)}=W_N^{2mn}$, so the two halves of the sum add. For odd outputs $W_N^{(2m+1)(n+N/2)}=-W_N^{(2m+1)n}$, so the second half is subtracted, and what is left of the exponent is $W_N^n$ times $W_{N/2}^{mn}$.

The graph is 14.1’s mirrored: inputs in order, the twiddle after the subtraction, outputs in bit-reversed order. In the figure, solid wires add and dashed wires are the ones subtracted.

The numbers beside the rows are the sums and the differences before their twiddles. The sums are 1, 3, 2, 0. The differences are 1, 1, −2, −2, and the twiddles $W_8^0,\dots,W_8^3$ then turn them into $1$, $0.707-0.707j$, $2j$ and $1.414+1.414j$.

Computing in place

Each butterfly reads two values and writes two values to the same two places. So the whole FFT needs no second array: it works in place. libsig’s fft.c does this. It reorders the input into bit-reversed order first and then overwrites it stage by stage.

What libraries choose

Which radix is best is a counting question. Radix-4 handles two stages of radix 2 at once, and its multiplies by $\pm j$ are free. Split-radix uses radix 2 for the even-indexed outputs and radix 4 for the odd-indexed ones. To compare them, count real operations. A multiply by $\pm1$ or $\pm j$ is free. A multiply by $e^{\pm j\pi/4}$ or $e^{\pm j3\pi/4}$ costs two multiplies and two adds. Any other twiddle costs four multiplies and two adds. The table shows the exact counts for $N=1024$.

Radix-4 saves 11.7 % of the operations of radix-2 and split-radix saves 15.0 %. These are small gains in total, but nearly all of the saving comes from multiplies. Radix-4 needs a length that is a power of 4. Split-radix works for every power of 2.

FFTW, and NumPy’s pocketfft, go further. They factor any $N$, use radix 2, 3, 4, 5 and larger blocks of hand-tuned code, and fall back on Rader’s method or Bluestein’s trick for large primes. FFTW also times several candidate plans on your machine and keeps the fastest. So $N$ need not be a power of two, but large primes are slow, and it is worth avoiding them. libsig’s FFT is radix-2 only, which is why the instruments on this site use powers of two.

Worked example

Every number below was recomputed in NumPy.

1. Factor counts, $N(N_1+\cdots)$ against $N^2$. For $N=12$: 84 against 144. For $N=30$: 300 against 900. For $N=32$: 320 against 1024. For $N=1000=2^3\cdot5^3$: 21 000 against 1 000 000. For $N=1024$: 20 480 against 1 048 576. For the prime $N=1009$: 1 018 081 either way.
2. Bluestein at $N=1009$. FFT size 2048. Three FFTs cost 33 792 complex multiplies, plus 2048 products and 2018 chirp multiplies: 37 858 against 1 018 081.
3. Two for one. The packed spectrum is $G=8,4,4,0,0,0,-4,4$. Its mirrored conjugate is $8,4,-4,0,0,0,4,4$. The results are $X=8,4,0,0,0,0,0,4$ and $Y=0,0,-4j,0,0,0,4j,0$.
4. DIF first stage on $x=1,2,0,-1,0,1,2,1$: sums $1,3,2,0$, differences $1,1,-2,-2$ before the twiddles.
5. Exact counts at $N=4096$, total real operations: radix-2 204 816, radix-4 179 552, split-radix 172 040. At $N=1024$ the totals are 40 976, 36 192 and 34 824.

Where you’ll meet this

Fast convolution (14.3) uses power-of-two FFTs to run long filters, and Goertzel and the chirp z-transform (14.4) comes back to the chirp for zooming into a band. Chapter 15 starts with Windowing (15.1), and spectrum analysers there take real-input FFTs.

Reference card