# Fast Fourier Transform
An asymptotically-improved \(O(N \log N)\) algorithm for implementing the Discrete Fourier Transform (naively \(O(N^2)\)).
# 1 Libraries
- FFTW, fast, but CPU only; also some possible issues with concurrency in shared plugin hosts.
# 2 Decimation In Time
Define \(DFT(N/T, t, x)\) to be the \(T\) size \(N/T\) DFTs with stride offset \(t \in [0, T)\).
Then the \(DFT(N/1, 0, x)\) can be decomposed into \(T\) smaller DFTs thus:
\[ \begin{aligned} X &= DFT(N/1, 0, x) \\ X_k &= \sum_{n=0}^{N-1} e^{- 2 \pi i n k / N} x_n \\ &= \sum_{t=0}^{T-1} \sum_{m=0}^{N/T-1} e^{- 2 \pi i (T m + t) k / N} x_{T m + t} \\ &= \sum_{t=0}^{T-1} e^{-2 \pi i k t / N} \sum_{m=0}^{N/T-1} e^{- 2 \pi i k m / (N / T)} x_{T m + t} \\ &= \sum_{t=0}^{T-1} e^{-2 \pi i k t / N} DFT(N/T, t, x) \\ Y_{t,k} &= DFT(N/T, t, x) \\ X_k &= \sum_{t=0}^{T-1} e^{-2 \pi i k t / N} Y_{t, k} \\ \end{aligned} \]
Continuing, it turns out the recombination is in fact a set of size \(T\) DFTs:
\[ \begin{aligned} X_k &= \sum_{t=0}^{T-1} e^{-2 \pi i k t / T} \frac{e^{-2 \pi i k t / N}}{e^{-2 \pi i k t / T}} Y_{t, k} \\ Z_{t,k} &= \frac{e^{-2 \pi i k t / N}}{e^{-2 \pi i k t / T}} Y_{t, k} \\ X_k &= \sum_{t=0}^{T-1} e^{-2 \pi i k t / T} Z_{t. k} \\ &= DFT(T, 0, Z_{t, k}) \end{aligned} \]
That is, X, the length N DFT of x, can be calculated by first doing T-many length N/T DFTs interleaved of x, multiplying by twiddle factors, then recombining with N/T-many length T DFTs. Inverse DFT follows the same process in reverse.
Applying this recursively gives an implementation of the Fast Fourier Transform algorithm, reducing the cost of DFT from \(O(N^2)\) to \(O(N \log N)\).
# 3 Parallel FFTW
Once you have “planned” your FFT using the FFTW library, you can execute it from another thread. One strategy for large FFTs is to split it into a set of smaller FFTs as above, and execute each in parallel (perhaps using OpenMP).
I implemented it in a phase vocoder time stretcher:
16 threads is fastest in terms of waiting, but is about 35% less efficient in terms of output duration vs CPU seconds than when using 2 threads (the minimum my current code supports). 8 threads is probably the sweet spot: 15% less efficient than 2 threads with waiting times 3x shorter (fast enough for realtime streaming).
# threads output/CPU seconds realtime speed
2 0.21 0.420
4 0.20 0.746
8 0.18 1.232
16 0.13 2.015
Compared with using FFTW3’s built-in OpenMP support:
# threads output/CPU seconds realtime speed
1 0.285 0.285
2 0.251 0.482
4 0.208 0.762
8 0.156 1.149
16 0.125 1.259
Throughput is higher at higher thread counts for manual splitting/threading, though it’s much harder to implement.