FFT 与窗函数：蝶形计算、复杂度与频谱泄漏

上一篇推导了 DFT 公式并手算了一个 $N=4$ 的例子。但直接按公式计算，$N$ 个输出每个都要做 $N$ 次乘加，总共 $N^2$ 次运算。$N = 1024$ 时是 100 万次，$N = 48000$ 时是 23 亿次——实时播放根本扛不住。

这一篇讲两件事：FFT 怎么把 DFT 加速到 $O(N\log N)$（用 N=2、N=4 完整演示），以及窗函数为什么必不可少（频谱泄漏从哪来、汉宁窗怎么推导）。

---

1. FFT 的核心思想：分而治之

1.1 旋转因子的周期性

记：

$$W_N = e^{-i2\pi/N}$$

它是”每 $N$ 步转一圈”的旋转因子。关键性质：

$$W_N^{k+N} = W_N^k, \qquad W_N^{2k} = W_{N/2}^{k}$$

第一个性质：指数以 $N$ 为周期；第二个性质：半长因子可以复用。FFT 就靠这两个性质，把大问题拆成小问题。

1.2 奇偶拆分

把 DFT 求和按 $n$ 的奇偶拆开：

$$X_k = \sum_{n \text{ even}} x[n] W_N^{kn} + \sum_{n \text{ odd}} x[n] W_N^{kn}$$

偶数项令 $n = 2m$，奇数项令 $n = 2m+1$：

$$X_k = \sum_{m=0}^{N/2-1} x[2m] W_N^{2mk} + W_N^k \sum_{m=0}^{N/2-1} x[2m+1] W_N^{2mk}$$

利用 $W_N^{2mk} = W_{N/2}^{mk}$：

$$X_k = \underbrace{\sum_{m=0}^{N/2-1} x[2m] W_{N/2}^{mk}}_{A_k} + W_N^k \underbrace{\sum_{m=0}^{N/2-1} x[2m+1] W_{N/2}^{mk}}_{B_k}$$

$A_k$ 和 $B_k$ 分别是偶数样本和奇数样本的 N/2 点 DFT。也就是说：一个大 DFT = 两个小 DFT 再合并。

---

2. N = 2 的蝶形

2.1 直接算

$N=2$ 时 $W_2^0 = 1$、$W_2^1 = -1$：

$$X_0 = x[0] + x[1]$$ $$X_1 = x[0] - x[1]$$

就两行，一加一减。这已经是一个完整的蝶形。

2.2 蝶形的图形含义

把输入 $x[0]$、$x[1]$ 画在左边，输出 $X_0$、$X_1$ 画在右边：

• 实线：$x[0] \to X_0$（加）；
• 实线：$x[1] \to X_0$（加）；
• 虚线：$x[1] \to X_1$（减）。

这就是”蝴蝶”形状的由来。

---

3. N = 4 的完整计算

3.1 拆成两个 N=2

输入 $x = [1, 2, 0, -1]$：

• 偶数样本：$x[0]=1,\ x[2]=0$；
• 奇数样本：$x[1]=2,\ x[3]=-1$。

偶数部分的 N=2 DFT：

$$A_0 = 1 + 0 = 1, \qquad A_1 = 1 - 0 = 1$$

奇数部分的 N=2 DFT：

$$B_0 = 2 + (-1) = 1, \qquad B_1 = 2 - (-1) = 3$$

3.2 合并

旋转因子（$N=4$）：$W_4^0 = 1$、$W_4^1 = -i$。

$$X_0 = A_0 + W_4^0 B_0 = 1 + 1 = 2$$ $$X_1 = A_1 + W_4^1 B_1 = 1 + (-i)\cdot 3 = 1 - 3i$$ $$X_2 = A_0 - W_4^0 B_0 = 1 - 1 = 0$$ $$X_3 = A_1 - W_4^1 B_1 = 1 + 3i$$

结果和上一篇直接按 DFT 公式手算完全一致。注意：蝶形法只做了 4 次复数乘加，直接法是 16 次。

---

4. 复杂度：为什么是 O(N log N)

设 $T(N)$ 是 N 点 FFT 的运算量。由拆分公式：

$$T(N) = 2T(N/2) + O(N)$$

含义：两个 N/2 点子问题，加上 N 次合并。展开：

$$T(N) = 2T(N/2) + N$$ $$= 2(2T(N/4) + N/2) + N = 4T(N/4) + 2N$$ $$= \dots = N \cdot T(1) + N\log_2 N = O(N\log N)$$

对比：

N	DFT 直接法 $N^2$	FFT $N\log_2 N$	加速比
16	256	64	4×
256	65536	2048	32×
1024	1048576	10240	102×
48000	23 亿	75 万	3072×

这就是为什么实时音频分析能用 FFT。

---

5. 频谱泄漏：直接截断的问题

5.1 问题

FFT 处理的是一段有限长的数据，等价于把无限长信号乘以一个”矩形窗”：

$$w_{\text{rect}}[n] = \begin{cases} 1, & 0 \le n < N \\ 0, & \text{其它} \end{cases}$$

如果信号的频率不是 bin 频率的整数倍，截断处波形”断”了。在频域里，这个断点会扩散成一大片旁瓣，能量从真正的频率泄漏到相邻 bin——这就是频谱泄漏（spectral leakage）。

5.2 数学表现

矩形窗的频谱是 $\text{sinc}$ 函数：

$$W_{\text{rect}}(f) = \frac{\sin(\pi N f/f_s)}{\pi f/f_s}$$

主瓣很窄但旁瓣很高（只衰减约 13 dB），所以泄漏明显。

---

6. 窗函数：让两端平滑

6.1 思路

与其让数据”哐当”一下截断，不如先把两端乘以一个平滑衰减到 0 的函数，再接缝就没了。加窗后的信号：

$$y[n] = x[n] \cdot w[n]$$

6.2 汉宁窗（Hann window）

最常用的窗：

$$w[n] = 0.5 - 0.5\cos\left(\frac{2\pi n}{N-1}\right), \qquad n = 0, 1, \dots, N-1$$

它两端为 0、中间为 1，正好是一个”半个余弦波抬起来”的形状。

图里：

• 绿色曲线是窗函数；
• 橙色是原始正弦；
• 白色是加窗后的信号：两端平滑衰减到 0，中间保持原样。

6.3 窗的代价

加窗后主瓣变宽（频率分辨率略降），但旁瓣大幅压低（约 31 dB），泄漏大大减少。这是”分辨率 vs 泄漏”的权衡。

常见窗对比：

窗	主瓣宽度	旁瓣衰减	用途
矩形	最窄	~13 dB	对瞬态敏感时
汉宁	中	~31 dB	一般音乐分析
汉明	中	~43 dB	语音
布莱克曼	宽	~58 dB	需要极低旁瓣

播放器用 AnalyserNode 时，浏览器内部已经应用了平滑窗（smoothing），但对自定义 FFT 就要自己加窗。

---

7. 与播放器的关系

播放器每一帧：

1. 从音频流取 1024 个采样点；
2. 乘以汉宁窗（或依赖浏览器内部处理）；
3. 做 FFT，得到 512 个 bin；
4. 按频段求和、取对数、归一化；
5. 指数平滑后驱动地形。

窗函数决定了频谱的”清晰度”：窗选得好，一个 440 Hz 的音只点亮一个 bin；窗选得差，会点亮一排相邻 bin，地形看起来”糊”。

---

8. 小结

1. FFT 利用旋转因子周期性，把大 DFT 拆成两个小 DFT 再合并；
2. N=2 蝶形：一加一减；N=4 是两个 N=2 加一次合并；
3. 复杂度从 $O(N^2)$ 降到 $O(N\log N)$；
4. 直接截断 = 乘矩形窗，产生频谱泄漏；
5. 汉宁窗 $w[n] = 0.5 - 0.5\cos(2\pi n/(N-1))$ 让两端平滑，压低旁瓣；
6. 窗是”分辨率 vs 泄漏”的权衡。

到这里，从声波到频谱的数学基础就完整了。后续系列文章会继续讲：频段特征与节拍检测、把频谱映射成 3D 地形、封面取色算法，以及可交互的正弦波混合实验。

8. 深入专题：递归树、位反转、窗的频响与常见误区

FFT 和窗函数还能再往下挖一层：递归树长什么样、位反转和就地计算是怎么回事、不同窗的频率响应差多少、以及 N=8 的完整步骤。

8.1 FFT 递归树

把 $T(N) = 2T(N/2) + N$ 画成树：

• 根节点是 $N$ 点问题；
• 每层拆成两个子问题；
• 叶子是 1 点问题（$X_0 = x[0]$，不需要计算）；
• 每层合并的总代价是 $N$；
• 层数 $\log_2 N$；
• 总代价 $N\log_2 N$。

8.2 位反转与就地计算

递归 FFT 在实现时通常改成迭代：先把输入按”位反转”重新排列，然后从最底层开始逐层合并。

$N=8$ 的位反转：3 位二进制反转。

索引	二进制	反转	结果
0	000	000	0
1	001	100	4
2	010	010	2
3	011	110	6
4	100	001	1
5	101	101	5
6	110	011	3
7	111	111	7

重排后，相邻两个点组成 N=2 蝶形，再相邻四个点组成 N=4，依此类推。整个过程在原数组上就地完成，不需要额外内存。

8.3 窗函数的频率响应对比

时域上看，三种窗的区别是两端”圆滑程度”：

• 矩形窗：两端突变；
• 汉宁窗：两端平滑到 0；
• 布莱克曼窗：更圆滑。

频域上的代价是主瓣宽度：

窗	主瓣宽度（bin）	旁瓣峰值
矩形	2	−13 dB
汉宁	4	−31 dB
汉明	4	−43 dB
布莱克曼	6	−58 dB

主瓣越宽，两个相近频率越难分开；旁瓣越低，远处泄漏越少。应用时按需选择：乐器分离用窄主瓣，一般频谱显示用汉宁。

8.4 N=8 蝶形步骤表

输入 $x = [x_0, x_1, \dots, x_7]$，位反转后顺序 $[x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7]$。

第 1 层（N=2）：

$$\begin{aligned} a_0 &= x_0+x_4, & a_1 &= x_0-x_4 \\ a_2 &= x_2+x_6, & a_3 &= x_2-x_6 \\ a_4 &= x_1+x_5, & a_5 &= x_1-x_5 \\ a_6 &= x_3+x_7, & a_7 &= x_3-x_7 \end{aligned}$$

第 2 层（N=4）：每 4 个一组，用 $W_4$ 合并。

第 3 层（N=8）：两组 N=4 用 $W_8$ 合并。

每层恰好 4 个蝶形 = $N/2$ 个，共 $\log_2 N = 3$ 层，总蝶形数 $N/2\log_2 N = 12$，每个蝶形一次复数乘法。这就是 $O(N\log N)$ 的具体来源。

8.5 常见误区

误区一：“FFT 会丢失信息”

FFT 是精确可逆变换，逆变换能无损还原原始时域信号（浮点误差除外）。

误区二：“窗函数越多越好”

加窗会损失两端数据（乘 0），主瓣变宽、分辨率下降。没有免费午餐：泄漏和分辨率只能取舍。

误区三：“fftSize 越大频谱越好”

fftSize 大 → 频率分辨率细，但时间分辨率差（需要更多采样时间），而且低频 bin 更密、高频 bin 没变。实时分析要在时间/频率分辨率之间平衡。

误区四：“频谱图上峰值对应精确频率”

不加窗或窗不好时，峰值可能落在两个 bin 之间，显示偏低的 bin。精细频率估计需要抛物线插值或更多采样。

误区五：“所有 FFT 实现都自动加窗”

AnalyserNode 内部有平滑处理，但自定义 FFT 必须自己加窗，否则泄漏明显。

8.6 扩展练习

1. 写出 N=8 位反转表并验证第 2 层合并的索引规律。
2. 计算 N=16 时蝶形总数，验证 $N/2\log_2 N$。
3. 为什么布莱克曼窗旁瓣最低但主瓣最宽？用窗函数公式解释。
4. 在播放器里把 FFT 长度从 1024 改成 4096，频率分辨率变化多少？时间窗口变多长？

练习

1. 用 $T(N) = 2T(N/2) + N$ 写出 $T(8)$ 的展开过程。
2. 为什么 $W_N^{2k} = W_{N/2}^k$？用定义验证。
3. 汉宁窗在 $n=0$ 和 $n=N-1$ 的值是多少？为什么这样能避免泄漏？
4. 播放器频谱显示”糊”或”抖”，分别可能是窗和哪一步的问题？