3. 快速傅里叶变换（FFT）：Cooley-Tukey算法详解
聊到数字信号处理，DFT（离散傅里叶变换）是绕不开的坎。但直接算DFT，复杂度是O(N²)。N稍微大一点，比如1024点，就得算一百多万次复数乘加。这在嵌入式系统里，基本是灾难。
好在1965年，Cooley和Tukey两位大神搞出了FFT。说白了，FFT就是DFT的一种快速算法，核心思想就四个字：分而治之。我当年第一次在定点DSP上实现FFT时，看着那蝶形运算的图，愣是画了一整天才搞明白数据是怎么“飞”的。今天咱们就把这层窗户纸捅破。
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²次。
Cooley-Tukey算法怎么做的？它把N点DFT拆成两个N/2点DFT。如果N是2的整数次幂，就可以一直拆下去，直到2点DFT。这就是“基2-FFT”的由来。
拆分的数学依据是旋转因子的两个性质：
- 周期性：W_N^{k+N} = W_N^k
- 对称性：W_N^{k+N/2} = -W_N^k
利用这两个性质，我们可以把N点序列按奇偶索引分成两组：
X[k] = Σ_{r=0}^{N/2-1} x[2r] · W_{N/2}^{kr} + W_N^k · Σ_{r=0}^{N/2-1} x[2r+1] · W_{N/2}^{kr}你看，前半部分和后半部分都是N/2点的DFT。记作：
X[k] = G[k] + W_N^k · H[k], k = 0, 1, ..., N/2-1那后半部分呢？利用对称性：
X[k+N/2] = G[k] - W_N^k · H[k], k = 0, 1, ..., N/2-1这就是蝶形运算的数学基础。一个加法，一个减法，中间乘个旋转因子。嗯，这里要注意，G[k]和H[k]都是N/2点的DFT，所以k的范围只有0到N/2-1。
核心结论：一个N点DFT可以分解为两个N/2点DFT，加上N/2个蝶形运算。递归下去，复杂度就从O(N²)降到了O(N log₂N)。
3.2 蝶形运算：FFT的基本单元
蝶形运算长什么样？我画个图你就明白了。
你看，两个输入A和B，经过一个旋转因子W，产生两个输出：A+W·B 和 A-W·W。形状像不像蝴蝶的两只翅膀？这就是“蝶形运算”名字的由来。
在8点FFT里，我们需要三级蝶形运算。每一级有4个蝶形单元。我习惯把每一级叫做一个“stage”，每个stage里旋转因子的分布是有规律的。
个人经验：我在调试定点DSP上的FFT时，经常因为旋转因子的索引搞错导致结果全乱。后来我养成了一个习惯——先用手算一个4点FFT验证流程，再上板子跑。4点FFT只有2级蝶形，手算一遍就能把数据流理清楚。
3.3 8点FFT完整流程
咱们以8点FFT为例，看看数据是怎么流动的。先按奇偶索引重排输入序列，这叫“位反转排序”。
| 原始索引 | 二进制 | 位反转 | 重排后索引 |
|---|---|---|---|
| 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 |
重排后的序列是：x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]。
然后进行三级蝶形运算。每一级旋转因子的分布规律如下：
- 第1级：蝶形间距=1，旋转因子只有W_2^0 = 1
- 第2级：蝶形间距=2，旋转因子有W_4^0, W_4^1
- 第3级：蝶形间距=4，旋转因子有W_8^0, W_8^1, W_8^2, W_8^3
我画个完整的8点FFT数据流图：
你想想看，8点DFT直接算要64次复数乘加。用FFT呢？每级4个蝶形，每个蝶形1次复数乘，3级就是12次复数乘。加上加法，总共也就24次复数操作。差距一目了然。
3.4 复杂度分析
咱们来算笔账。N点基2-FFT：
- 级数：log₂N 级
- 每级蝶形数：N/2 个
- 每个蝶形：1次复数乘 + 2次复数加
- 总复数乘：(N/2) · log₂N 次
- 总复数加：N · log₂N 次
对比一下DFT和FFT的复杂度：
| N | DFT复数乘 | FFT复数乘 | 加速比 |
|---|---|---|---|
| 16 | 256 | 32 | 8x |
| 64 | 4096 | 192 | 21x |
| 256 | 65536 | 1024 | 64x |
| 1024 | 1048576 | 5120 | 205x |
| 4096 | 16777216 | 24576 | 683x |
N越大，FFT的优势越明显。1024点时，加速比已经超过200倍。这就是为什么实时频谱分析、OFDM调制解调这些应用，全都离不开FFT。
避坑指南：我曾经在一个项目中，直接用浮点FFT库处理4096点数据，结果在定点DSP上跑得奇慢无比。后来才发现，定点DSP做浮点运算要软模拟，效率极低。换成定点FFT后，速度提升了近10倍。所以选FFT实现方式时，一定要看你的硬件平台支不支持硬件浮点。
3.5 代码实现：C语言版8点FFT
光说不练假把式。我写个8点FFT的C代码，你感受一下蝶形运算在代码里长什么样。
#include <math.h>
#include <complex.h>
#define N 8
// 位反转函数
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));
// 位反转重排
for (int i = 0; i < n; i++) {
int j = bit_reverse(i, log2n);
if (j > i) {
complex double tmp = x[i];
x[i] = x[j];
x[j] = tmp;
}
}
// 蝶形运算
for (int len = 2; len <= n; len <<= 1) {
int half = len >> 1;
double angle = -2 * M_PI / len;
complex double w = cos(angle) + I * sin(angle);
for (int i = 0; i < n; i += len) {
complex double wn = 1 + 0 * I;
for (int j = 0; j < half; j++) {
complex double u = x[i + j];
complex double v = x[i + j + half] * wn;
x[i + j] = u + v;
x[i + j + half] = u - v;
wn *= w;
}
}
}
}这段代码的核心就是三层循环：最外层控制级数（len从2到N），中间层控制蝶形组，最内层控制组内的蝶形。你仔细看内层循环，就是标准的蝶形运算：u+v 和 u-v。
优化建议：在实际工程中，旋转因子可以预先计算好存成查找表，不用每次实时算cos/sin。我做过测试，查表法比实时计算快3-5倍。另外，如果N是固定的，位反转的映射表也可以提前算好存起来。
3.6 小结
Cooley-Tukey算法把DFT从O(N²)降到了O(N log₂N)，靠的就是分而治之的思想和蝶形运算这个基本单元。理解FFT的关键，我觉得就三点：
- 拆分逻辑：按奇偶索引把N点DFT拆成两个N/2点DFT
- 蝶形结构：每个蝶形做一次复数乘和两次复数加
- 位反转：输入数据要先按位反转重排，才能进行原位计算
搞懂了这些，FFT在你眼里就不再是黑盒子了。下次在示波器上看到频谱分析结果时，你可以跟同事说：“嗯，这背后就是蝶形运算在跑。”