6.2 小时变成 29 毫秒
前面十六章有一个没说破的前提：你得算得动。按定义式，一次 N 点变换要 N² 次运算，2²⁰ 个点就是一万一千亿次——本机实测要跑 6.2 小时。要真是这样，前面所有的东西都是纸上谈兵。这一章讲那个改变了一切的算法，它的核心只有一句话，而且你已经知道那句话的一半了。
本机（一台笔记本，Node 24，同一份手写的、没有任何优化的代码）按定义式算一次 2048 点的变换，用了 84.6 毫秒。
问：同一台机器、同一份代码风格，算 2²⁰ = 1048576 个点的变换，走 FFT 要多久？
先看成本
N N²（定义式） N·log₂N（FFT） 倍数 ──────────────────────────────────────────────────────────── 1024 1.049 × 10⁶ 1.024 × 10⁴ 102.4 65536 4.295 × 10⁹ 1.049 × 10⁶ 4096.0 1048576 1.100 × 10¹² 2.097 × 10⁷ 52428.8 16777216 2.815 × 10¹⁴ 4.027 × 10⁸ 699050.7
这不是「快一点」，是换了一个复杂度类。N 越大差得越离谱，而现实中的 N 全都很大：一秒钟的 CD 音频是 44100 个点，一张 4K 图是八百万个像素。
核心思想：拆成两半
问题：算 N 个转速的质心，每个要 N 次乘加，一共 N² 次。
Cooley–Tukey 的观察：把 N 个采样点按下标的奇偶分成两堆，那么整个变换可以由两堆各自的（半长的）变换拼出来。
N−1
X[k] = Σ x[n] · e^(−i2πkn/N)
n=0
按奇偶拆：
N/2−1 N/2−1
= Σ x[2m]·e^(−i2πk·2m/N) + Σ x[2m+1]·e^(−i2πk(2m+1)/N)
m=0 m=0
注意关键的一步： e^(−i2πk·2m/N) = e^(−i2πkm/(N/2))
└───────┬───────┘ └────────┬────────┘
「N 点圈上走 2m 步」 = 「N/2 点圈上走 m 步」
于是：
X[k] = E[k] + e^(−i2πk/N) · O[k]
└┬┘ └┬┘
偶数点的 N/2 变换 奇数点的 N/2 变换
E 和 O 都只有 N/2 个不同的值（它们是 N/2 周期的），
所以算完 E 和 O，**每个 k 只要再花一次乘法和一次加法**。
这就是全部。用第 2 章那台机器的语言说：
把一次「N 点的绕圈」拆成两次「N/2 点的绕圈」，再花 N 次乘加把它们拼起来。
拆之前： N² 拆一次： 2 · (N/2)² + N N = 8： 8² = 64 → 2 × 4² + 8 = 40 少了三分之一 然后对那两个 N/2 再拆一次，再拆一次……拆 log₂N 层， 每层的拼接成本都是 N。 总成本： N · log₂N
请注意这里没有任何近似、没有任何取舍。FFT 算出来的和定义式算出来的是同一个数（浮点误差意义上，FFT 还更准一点——见下面那条「✗」）。它纯粹是把重复的工作认出来，只做一遍。
本机实测
同一台机器，同一份手写实现（Node 24，未做任何优化）
N = 2048
按定义式 84.6 ms
FFT 0.029 ms 快 2922 倍
N = 2²⁰ = 1048576
FFT ⟨29.0 ms⟩
按定义式（外推） 22170 秒 = ⟨6.2 小时⟩
也就是说：**FFT 算一百万个点，比定义式算两千个点还快三倍。**
实测的 7.65 × 10⁵ 倍比理论的 52429 倍还大得多，为什么？因为理论倍数只数了乘加次数，而定义式的每一项都要现算一次 cos 和 sin，那是比乘法贵得多的函数调用；FFT 则把旋转因子递推出来，几乎不调三角函数。常数因子这一层，FFT 又赢了一大截。（顺带说：这也意味着上面那个「6.2 小时」是个粗糙的外推——精确到一个数量级足够说明问题了。）
1965，和 1805
Cooley 和 Tukey 那篇论文发表于 1965 年，四页纸，是二十世纪引用最多的算法论文之一。它的背景很值得知道：
Tukey 当时在肯尼迪的科学顾问委员会里，讨论的议题是如何监测苏联的地下核试验——办法是在境外布设地震仪，用地震波的频谱来区分核爆和天然地震。而当时的计算机算不动那么长的记录。Tukey 在一次会上把这个分治的想法写在了纸上，IBM 的 Cooley 把它实现出来。
更有意思的是后来的考证：高斯在 1805 年就已经写下了同一个算法，用来从有限次观测推算小行星（谷神星和智神星）的轨道。那份手稿用拉丁文写成，没有发表，直到他死后才收进全集第三卷——而且比傅里叶 1807 年那篇奠基论文还早两年。
它被埋没了一百六十年，不是因为难，是因为当时没有值得用它去算的东西。手工计算的年代，N 顶多几十，N² 和 N log N 没有区别。这个算法必须等到「有大量数据 + 有机器」这两个条件同时成立，才成其为一件大事。
「FFT」不是一种变换，是算 DFT 的一族快速算法。这一点很多人搞混，说「用 FFT 还是用 DFT」——那就像问「用快速排序还是用排序」。
「一族」也是实话，它不止一个：
- 基 2 Cooley–Tukey：上面讲的这个，要求 N 是 2 的幂。
- 混合基：N = 2ᵃ·3ᵇ·5ᶜ 时按小质因数逐层拆。FFTW、numpy 都这么干，所以 N = 1000（= 2³·5³）跑得飞快，N = 1009（质数）就慢。
- Bluestein / Rader：N 是大质数时，用一个技巧把 DFT变成一次卷积，再用第 14 章的定理把那次卷积补零到 2 的幂来算。绕了一大圈，但仍然是 O(N log N)。
实践上的建议：如果你能自由选长度，选 2 的幂或者只含小质因数的数。numpy 的 next_fast_len 就是干这个的。同样是一万左右的长度，选 10000（= 2⁴·5⁴）和选 10007（质数），实测能差好几倍。
一个能跑的递归 FFT，十二行——写一次比读十遍有用：
import cmath
def fft(x):
N = len(x)
if N == 1:
return list(x)
E = fft(x[0::2]) # 偶数下标
O = fft(x[1::2]) # 奇数下标
out = [0] * N
for k in range(N // 2):
w = cmath.exp(-2j * cmath.pi * k / N) * O[k]
out[k] = E[k] + w # 蝶形运算
out[k + N // 2] = E[k] - w # 上下两半只差一个符号
return out
print([round(v.real, 6) + round(v.imag, 6)*1j for v in fft([1, 2, 3, 4])])
# [(10+0j), (-2+2j), (-2+0j), (-2-2j)] ← 和第 1 章那个结果一样
最后那两行 out[k] = E[k] + w 和 out[k+N/2] = E[k] − w 叫蝶形运算（因为数据流图长得像蝴蝶）。整个 FFT 就是这只蝴蝶重复 N/2 · log₂N 次。
再亲手量一次那个倍数：
import numpy as np, time
for p in (10, 13, 16, 20):
N = 2**p
x = np.random.randn(N)
t0 = time.perf_counter(); np.fft.fft(x); t1 = time.perf_counter()
print(f'N=2^{p:2d}={N:8d} FFT {1000*(t1-t0):7.2f} ms '
f'定义式要 {N*N/1e9:8.1f} × 10⁹ 次乘加')
python3 -c "import numpy as np;print(np.round(np.fft.fft([1,2,3,4]),6))"
在线跑：第一段纯标准库（cmath），python.org/shell 就行。想看蝴蝶图：Wikipedia 的 Cooley–Tukey FFT algorithm 词条里那张 8 点的数据流图，配着上面那段代码读，二十分钟能彻底通。
- 你正在听的这首歌。MP3 在 44100 Hz 立体声、窗长 1152、50% 重叠下，每秒要做 153 次变换。这还只是编码；解码端也要做同样多次。没有 FFT，随身听这个品类不会存在。
- 你手机的每一格信号。4G/5G 的 OFDM 收发两端，每个符号都要做一次 FFT/IFFT，一秒钟上万次。「OFDM 在 1960 年代就被提出，1990 年代才实用」，中间隔着的就是能不能实时算 FFT。
- 医院里的 CT 和 MRI。CT 的滤波反投影、MRI 的 k 空间重建，核心都是二维 FFT。一次 MRI 扫描要重建几百层，每层一次二维变换。
- 射电天文与 SETI。射电望远镜阵列的相关处理、脉冲星搜索、地外信号搜索，全是海量 FFT。这也是最早一批把 FFT 做进专用硬件的领域。
- 大整数乘法与密码学。第 14 章说过：大数相乘 = 数位卷积 = 频域相乘。你的 TLS 握手在算 RSA 或者椭圆曲线时，底层的大数库里跑的就是这个算法。
- Shazam、语音助手、鸟鸣 App、地震预警、示波器、频谱仪、电网监测、雷达、声呐、天气雷达……这个列表可以一直写下去。Top 10 Algorithms of the 20th Century 那份名单上，FFT 是唯一一个普通人每天在用几千次却从没听说过名字的。
「FFT 是 DFT 的快速近似。追求精度的时候应该用定义式老老实实算。」
FFT 不是近似，它算的是同一个数学量，中间没有任何截断或省略。而在浮点数意义上，FFT 通常比定义式更准，不是更差。
原因很直白：定义式对每个输出要做 N 次连加，误差按 O(√N) 甚至 O(N) 累积（这是《抹零》那本书的地盘）；FFT 只有 log₂N 层，每层一次蝶形，误差按 O(√log N) 累积。N = 2²⁰ 时，一个是二十次加法的深度，一个是一百万次。此外定义式每项现算 cos/sin，那本身也在引入误差，而 FFT 的旋转因子可以预先算好或用高精度递推。
真要挑 FFT 的毛病，只有一条：它对 N 的因数敏感。N 是大质数时，实现要绕道 Bluestein，常数因子会变大，误差也略增。判据是：担心精度的时候该做的是选一个好的长度，而不是退回定义式。
正确答案是 C：29 毫秒。比同一台机器用定义式算 2048 个点（84.6 ms）还快三倍。
A、B、D 三个错法都是在做线性或接近线性的外推：「点数是 512 倍，那时间大概也该是几百倍」。这个直觉对 O(N) 的算法成立，对这里不成立——但它错的方向很有意思：它高估了 FFT 的耗时，因为它默认「更多数据 = 更多时间」这个线性关系。而 N·log₂N 相对于 N² 的优势，恰恰在于它几乎就是线性的（log 那一项长得极慢：N 从 2048 涨到 2²⁰，log₂N 只从 11 涨到 20）。 值得把这件事记成一句话：FFT 让「做频谱分析」的成本从「和数据量的平方成正比」变成了「几乎和数据量成正比」。一个几乎线性的操作，可以放进任何实时流水线；一个平方的操作不能。这就是为什么这个算法改变的不是某个领域的效率，而是哪些领域能存在。 顺带一提，如果你选了 A，你其实猜对了另一个问题的答案：本机 FFT 算 2²⁰ 点如果不用分治而是朴素实现，需要的正是 22170 秒。你的直觉估对了数据量的增长，只是没把算法的改进算进去。这一次换域给出的不是简化，是算法从哪儿来的。分治不是随便试出来的，它是这个圈的对称性逼出来的。这也解释了为什么 N 是 2 的幂时最漂亮（每一层都能对半分），以及为什么 N 是大质数时要绕道（那个圈里没有小一号的圈可以用，只好先把它变成卷积）。
凡是有这种自相似对称性的地方，通常都藏着一个分治算法。这条经验的适用范围远比信号处理广——Strassen 的矩阵乘法、Karatsuba 的大数乘法、数论变换（NTT），都是同一个模式的产物。
这一章的一句话
FFT 不是一种变换，是「N 点的圈里藏着一个 N/2 点的圈」这个观察的直接后果；它把频谱分析的成本从平方级压到了几乎线性，而这一步决定的不是某些事情能做多快，是哪些事情能存在。
下一章处理一个具体到不能再具体的问题：你手机相册里那几万张 JPEG，为什么用的不是傅里叶变换？答案藏在一个很容易被忽略的细节里——DFT 会把一块图像的左边和右边接起来，而接缝上凭空多出的那道跳变，本机量出来是 21。