3、快速傅里叶变换：Cooley-Tukey算法原理与蝶形运算
好，咱们今天来啃一块硬骨头——快速傅里叶变换。说实话，我刚入行那会儿，看到FFT的推导公式就头大。但后来在项目中做实时频谱分析，发现不用FFT根本跑不动，硬着头皮啃下来才发现，这东西其实很巧妙。
FFT的核心思想，说白了就是四个字：分而治之。Cooley-Tukey算法正是这个思想的经典实现。你想想看，一个N点的DFT，直接算要O(N²)的复杂度，N=1024时就是一百多万次运算。但用FFT，复杂度降到O(N log₂N)，同样1024点，只需要一万次左右。差了整整两个数量级。
3.1 从DFT到FFT：问题在哪？
先回顾一下DFT的公式：
X[k] = Σ_{n=0}^{N-1} x[n] · W_N^{kn}, k = 0,1,...,N-1其中W_N = e^{-j2π/N}，我们叫它旋转因子。每个X[k]都要做N次复数乘法和N-1次复数加法。N个点就是N²次乘法。N=1024时，就是一百万次。嗯，这个量级在嵌入式系统里基本是灾难。
我在做一款便携式振动分析仪时就遇到过这个问题。用STM32F4直接算1024点DFT，采样率10kHz，结果CPU占用率飙到90%以上，其他任务全被卡死。后来换成FFT，CPU占用率直接降到15%。这就是差距。
3.2 Cooley-Tukey算法的核心思想
Cooley-Tukey算法的精髓，就是利用旋转因子的对称性和周期性，把大DFT拆成小DFT。
旋转因子有两个重要性质：
- 对称性：W_N^{k+N/2} = -W_N^k
- 周期性：W_N^{k+N} = W_N^k
基于这两个性质，我们可以把N点DFT拆成两个N/2点DFT。怎么拆？按奇偶位置拆分输入序列。
假设N是2的整数次幂（这是最常用的场景），我们把x[n]分成：
- 偶数索引：x[0], x[2], x[4], ... → 记作x_even[r]
- 奇数索引：x[1], x[3], x[5], ... → 记作x_odd[r]
然后DFT公式可以写成：
X[k] = Σ_{r=0}^{N/2-1} x[2r] · W_N^{2rk} + Σ_{r=0}^{N/2-1} x[2r+1] · W_N^{(2r+1)k}
= Σ_{r=0}^{N/2-1} x_even[r] · W_{N/2}^{rk} + W_N^k · Σ_{r=0}^{N/2-1} x_odd[r] · W_{N/2}^{rk}
= X_even[k] + W_N^k · X_odd[k]这里X_even[k]和X_odd[k]分别是N/2点的DFT。你看，一个N点DFT变成了两个N/2点DFT，外加一次复数乘法。
关键点：这个分解可以递归进行。N/2点DFT又可以拆成两个N/4点DFT，一直拆到2点DFT为止。2点DFT的计算非常简单，不需要乘法。
3.3 蝶形运算：FFT的基本单元
刚才的公式里，X_even[k]和X_odd[k]都是N/2点的DFT，所以k的范围是0到N/2-1。那k从N/2到N-1的部分怎么办？
利用旋转因子的周期性：
X_even[k + N/2] = X_even[k]
X_odd[k + N/2] = X_odd[k]再利用对称性：W_N^{k+N/2} = -W_N^k
所以：
X[k] = X_even[k] + W_N^k · X_odd[k], k = 0,1,...,N/2-1
X[k + N/2] = X_even[k] - W_N^k · X_odd[k], k = 0,1,...,N/2-1这个结构就是蝶形运算。你画出来看看，像不像一只蝴蝶的两只翅膀？
我的经验：刚开始学蝶形运算时，我建议你亲手画一遍8点FFT的蝶形图。从输入到输出，每一级都标清楚旋转因子。画完你就彻底理解了。我在培训新人时，这是必做的练习。
3.4 8点FFT的完整流程
咱们以8点FFT为例，走一遍完整流程。N=8，需要log₂8=3级蝶形运算。
第一级：把8点DFT拆成两个4点DFT
- 偶数索引：x[0], x[2], x[4], x[6]
- 奇数索引：x[1], x[3], x[5], x[7]
第二级：每个4点DFT再拆成两个2点DFT
- 偶数索引的偶数：x[0], x[4]
- 偶数索引的奇数：x[2], x[6]
- 奇数索引的偶数：x[1], x[5]
- 奇数索引的奇数：x[3], x[7]
第三级：2点DFT直接计算
每一级的计算量：N/2次复数乘法和N次复数加法。3级总共：3 × (4次乘法 + 8次加法) = 12次乘法 + 24次加法。
对比直接DFT：8²=64次乘法，56次加法。差距已经很明显了。
| N | 直接DFT乘法次数 | FFT乘法次数 | 加速比 |
|---|---|---|---|
| 8 | 64 | 12 | 5.3x |
| 64 | 4096 | 192 | 21.3x |
| 1024 | 1,048,576 | 5120 | 204.8x |
| 4096 | 16,777,216 | 24576 | 682.7x |
看到没？N越大，加速效果越恐怖。N=4096时，加速比接近700倍。这就是为什么所有实时信号处理系统都在用FFT。
3.5 位反转排序
刚才的分解过程里，输入序列被反复按奇偶拆分。最终你会发现，输入的顺序被打乱了。比如8点FFT，输入顺序变成了：
原始索引: 0 1 2 3 4 5 6 7
FFT输入: 0 4 2 6 1 5 3 7这个顺序怎么来的？把原始索引写成二进制，然后反转比特位：
- 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
这就是位反转排序。在实现FFT时，我们通常先对输入做位反转重排，然后从第一级开始逐级做蝶形运算。
注意：我曾经在实现时犯过一个低级错误——忘了做位反转，结果输出全乱套了。排查了半天才发现是输入顺序没处理好。所以写代码时，位反转这一步一定要单独测试验证。
3.6 代码实现：C语言版基2FFT
下面给一个最简化的基2FFT实现。这个版本只做原位计算，节省内存。
#include <math.h>
#include <complex.h>
// 位反转函数
unsigned int bit_reverse(unsigned int x, int log2n) {
unsigned int y = 0;
for (int i = 0; i < log2n; i++) {
y = (y << 1) | (x & 1);
x >>= 1;
}
return y;
}
// 基2 FFT，原位计算
void fft(complex double *x, int n) {
int log2n = (int)(log(n) / log(2) + 0.5);
// 位反转重排
for (int i = 0; i < n; i++) {
int j = bit_reverse(i, log2n);
if (i < j) {
complex double tmp = x[i];
x[i] = x[j];
x[j] = tmp;
}
}
// 蝶形运算
for (int len = 2; len <= n; len <<= 1) {
double angle = -2 * M_PI / len;
complex double wlen = cos(angle) + I * sin(angle);
for (int i = 0; i < n; i += len) {
complex double w = 1 + 0 * I;
for (int j = 0; j < len / 2; j++) {
complex double u = x[i + j];
complex double v = x[i + j + len / 2] * w;
x[i + j] = u + v;
x[i + j + len / 2] = u - v;
w *= wlen;
}
}
}
}这段代码的核心就是三层循环：最外层控制蝶形级数，中间层控制每组蝶形的起始位置，最内层做具体的蝶形运算。每次迭代，len翻倍，代表当前处理的DFT长度。
优化建议：在实际项目中，我一般会预计算所有旋转因子，存成一个查找表。这样能省掉每次计算cos/sin的开销。对于实时系统，这个优化能带来20%-30%的性能提升。
3.7 避坑指南
最后分享几个我在项目中踩过的坑：
- 输入长度必须是2的幂：基2FFT要求N=2^m。如果数据长度不是2的幂，需要补零到下一个2的幂。我曾经在音频处理中忘了补零，结果频谱里出现了奇怪的伪影。
- 旋转因子的精度：在定点DSP上实现时，旋转因子的量化误差会累积。我建议用Q15或Q31格式，并做适当的舍入处理。
- 原位计算的陷阱：蝶形运算中，u和v的计算顺序很重要。必须先保存u的原始值，再更新。否则会覆盖掉后面要用到的数据。
- 复数乘法的实现：一个复数乘法需要4次实数乘法和2次实数加法。在硬件加速时，可以用3次乘法+5次加法的优化方案，减少乘法器资源。
嗯，Cooley-Tukey算法和蝶形运算就讲到这里。下一章咱们聊聊FFT在卷积神经网络里的具体应用——怎么用FFT加速卷积计算。那个话题更有意思，到时候见。