快速傅里叶变换
概述
基于分治策略（divide and conquer）的算法，将离散傅里叶变换（DFT）的计算复杂度从O(n²)降低到O(n log n)，使大规模频谱分析从理论上可行变为实践中可行，从根本上改变了信号处理、通信、图像处理等领域的计算面貌。
关键内容
核心算法（Radix-2 Cooley-Tukey）
DFT 定义：$X_k = \sum_{j=0}^{n-1} x_j \omega_n^{jk}$，$\omega_n = e^{-2\pi i/n}$，直接计算 $O(n^2)$
分治分解（$n = 2M$，按奇偶拆分）：
$$X_k = \underbrace{\sum_{j=0}^{M-1} x_{2j}\omega_M^{jk}}{A_k} + \omega_n^k \underbrace{\sum{j=0}^{M-1} x_{2j+1}\omega_M^{jk}}_{B_k}$$
蝴蝶运算（butterfly operation）：利用周期性 $A_{k+M}=A_k$，$\omega_n^{k+M}=-\omega_n^k$：
$$X_k = A_k + \omega_n^k B_k, \quad X_{k+M} = A_k - \omega_n^k B_k$$
一次乘法（$\omega_n^k B_k$）同时产生两个输出——算法名称"蝴蝶"来自数据流图的形状。
复杂度
递推：$T(n) = 2T(n/2) + O(n)$ → $T(n) = O(n\log n)$
| $n$ | 直接 DFT（$n^2$） | FFT（$\frac{n}{2}\log_2 n$） | 加速比 |
|---|---|---|---|
| 1024 | 1,048,576 | 5,120 | 205× |
| 65536 | $4\times 10^9$ | 524,288 | 4096× |
| $10^6$ | $10^{12}$ | $10^7$ | $\approx$50000× |
算法要素
- 位反转排列（bit-reversal permutation）：就地（in-place）实现的输入重排，将 $j$ 替换为其二进制反转
- 旋转因子（twiddle factor）：$\omega_n^k$ 预计算存于查找表，避免重复计算三角函数
- 递归结构：$n = 2^m$ 时共 $m = \log_2 n$ 层，每层 $n/2$ 个蝴蝶，总 $O(n\log n)$
历史背景
冷战动因：1963年《部分禁止核试验条约》签署后，监测地下核试验需实时频谱分析全球地震数据。Richard Garwin（IBM物理学家，美国总统科学顾问委员会成员）促成了 Tukey-Cooley 合作。FFT 的诞生有深刻的地缘政治背景。
高斯先驱（鲜为人知）：卡尔·弗里德里希·高斯 早在 1805 年计算小行星轨道时就在手稿中描述了类似的快速分解方法，但直到 1866 年才以拉丁文遗著发表，未被传播——FFT 的核心思想被独立重发现了。
Danielson-Lanczos 引理（1942）：康尼利厄斯·朗佐斯与 Danielson 已将 DFT 分解表示为两个半长 DFT 的组合，但未递归到 $O(n\log n)$。
主要应用领域
| 领域 | 应用 |
|---|---|
| 通信 | OFDM（WiFi/4G/5G）、信道估计、GPS信号 |
| 图像处理 | JPEG压缩（DCT）、CT图像重建、MRI的k空间重建 |
| 音频 | MP3/AAC编码（MDCT）、音乐识别（Shazam频谱指纹）、语音识别 |
| 科学计算 | 谱方法求解PDE、分子动力学（Ewald PME）、天体物理 |
| 算法 | 快速多项式乘法 $O(n\log n)$、大整数乘法（Schönhage-Strassen） |
变体与改进
- Radix-4 / Split-Radix FFT：进一步减少乘法次数
- Bluestein 算法：处理任意长度（含素数长度）的 DFT
- 实数 FFT：利用实数信号对称性节省一半计算量
- FFTW（Fastest Fourier Transform in the West，MIT）：自适应选择最优实现的开源库
来源
- raw/books/数值分析/21_cooley_tukey_fft.md