DIT基-2 FFT原理与实现：蝶形运算、位序倒置及工程细节

做信号处理这些年，我碰到的第一个"看起来不难、但没多少人敢说自己完全搞懂"的算法，大概就是快速傅里叶变换（FFT）。标题里的"按时间抽取（DIT）基-2 FFT"，是所有FFT版本里用得最多、也最基础的一种：它把一个N点DFT按时间下标n的奇偶关系拆成两个N/2点DFT，再继续拆下去，直到剩下最基本的2点蝶形运算。别看公式就一页纸，里面的分治思路、位序倒置、蝶形流图，每一处都是后面理解其他FFT变种、甚至各种开源频谱分析库的地基。

一次调试音频频谱可视化工具，我用最朴素的DFT双重循环去算1024点频谱，循环里全是复数乘法，CPU占用直接拉满，帧率掉到没法看。后来才发现直接算DFT要N²次复数乘法，而基-2 FFT只需要(N/2)log₂N次，1024点就差了两百多倍。这篇文章就从"FFT到底优化掉了什么"讲起，再手推一遍DIT基-2的蝶形原理，给一份可以直接跑的标准实现，最后聊工程里没人提醒你的采样率、窗函数和归一化问题。适合正在学数字信号处理的学生，也适合想把底层搞明白的工程师。

1. 从N²到NlogN：FFT到底省掉了什么

1.1 直接DFT的计算瓶颈

先回到定义。N点离散傅里叶变换（DFT）公式是：

[ X(k)=\sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N},\quad k=0,1,\dots,N-1 ]

把 (e^{-j2\pi nk/N}) 记作旋转因子 (W_N^{nk})，DFT就是 (X(k)=\sum x(n)W_N^{nk})。写两层循环去算，每个输出点都要做N次复数乘法，N个输出点总共是N²次复数乘法，加上N(N-1)次复数加法。

数据量小的时候无所谓，一到工程场景就露馅。1024点DFT要一百多万次复数乘法，4096点就要一千六百多万次，如果做的是流式频谱分析，比如每帧4096点、每秒还要更新几十帧，普通单线程几乎跑不动。为什么音频可视化、振动分析、雷达信号处理这些场景几乎离不开FFT，原因就在这个计算量差距。

基-2 FFT借助分治思想把复杂度压到 (O(N\log_2 N)) 的量级。复数乘法次数大约是 ((N/2)\log_2 N) 次，加法和旋转因子运算也在同一量级。直接感和直觉为什么差距这么大，看这张表就清楚了：

点数越大，优势越离谱。65536点在普通电脑上用FFT只需要几毫秒到几十毫秒，直接DFT则是几十亿次乘法，等得人心态崩。

1.2 旋转因子里的三个性质，才是分治能成立的支点

FFT不是靠"魔法"把计算变快，而是靠旋转因子的三个代数性质：

• 周期性：(W_N^{k+N}=W_N^k)，所以同一旋转因子在多个位置能复用
• 对称性：(W_N^{k+N/2}=-W_N^k)，这一条直接诞生了蝶形运算里的加/减结构
• 可约性：(W_N^2=W_{N/2})，意思是旋转因子角度翻倍，等价于把变换点数N减半

可约性尤其关键。DIT能不断把N点DFT拆成N/2点DFT，靠的就是这个性质。很多教材直接抛结论"按奇偶抽取"，但没解释为什么能用 (W_{N/2}^{rk}) 代替 (W_N^{2rk})。其实道理很简单：(e^{-j2\pi(2r)k/N}=e^{-j2\pi rk/(N/2)})，分母从N变成N/2，所以原来的N点旋转因子就变成了N/2点旋转因子，推导中自然而然地降维了。

这三个性质不是数学游戏，它们决定了蝶形运算的硬件/软件实现方式：对称性让减法变成加法，可约性让上一层算过的 (W_{N/2}^{k}) 能直接复用，周期性让旋转因子的查找表只需要存N/2个值。我后来写优化版本时把旋转因子预先算好放在表里，省掉了大量实时三角函数计算，靠的就是周期性这个性质。

2. 蝶形运算与位序倒置：DIT基-2的原理是怎么来的

2.1 按时间抽取到底"抽取"了什么

DIT的全称是Decimation In Time，按时间抽取。这里的抽取就是把原始序列x(n)按下标n的奇偶分成两组：

[ \text{偶数下标：}x(0),x(2),x(4),\dots=x(2r) ] [ \text{奇数下标：}x(1),x(3),x(5),\dots=x(2r+1) ]

其中 (r=0,1,\dots,N/2-1)。把DFT公式按这个分组展开：

[ X(k)=\sum_{r=0}^{N/2-1}x(2r)W_N^{2rk}+\sum_{r=0}^{N/2-1}x(2r+1)W_N^{(2r+1)k} ]

利用可约性 (W_N^{2rk}=W_{N/2}^{rk})，第二项再把公因子 (W_N^k) 提出来：

[ X(k)=\underbrace{\sum_{r}x(2r)W_{N/2}^{rk}}{G(k)}+W_N^k\underbrace{\sum{r}x(2r+1)W_{N/2}^{rk}}_{H(k)} ]

这里G(k)是原始序列偶数下标那N/2个点构成的N/2点DFT，H(k)是奇数下标那N/2个点构成的N/2点DFT。因为G和H都是以N/2为周期的周期函数，所以可以用前半段的G(k)、H(k)算出后半段的X(k)：

[ X(k)=G(k)+W_N^kH(k),\quad k=0,1,\dots,N/2-1 ] [ X(k+N/2)=G(k)-W_N^kH(k),\quad k=0,1,\dots,N/2-1]

后半段用到的关键就是对称性：(W_N^{k+N/2}=-W_N^k)。于是"两个N/2点DFT再加加减减"一次能算两个输出点，数据流画出来就像一只蝴蝶，所以叫蝶形运算。

更直白一点：原来N点DFT直接算N²次乘法，现在拆成两个N/2点DFT需要约2×(N/2)²次，外加N/2个旋转因子乘法，总量大约 (N^2/2+N/2)。看起来只省了约一半，但迭代拆分下去就不得了，每拆一层规模减半，一共拆 (\log_2 N) 层，最后成本变成Nlog₂N量级。

2.2 手推一个N=4的例子，胜过背十遍公式

拿N=4、输入 (x=[1,0,-1,0]) 验证整套推导。这个信号是 (x(n)=\cos(2\pi n/4)) 一个周期的采样，频率正好落在k=1上，理论频谱应该是k=1和k=3处各有幅度2，其余为零。

第一步，按奇偶拆。偶数下标取出 (x(0)=1,x(2)=-1)，做2点DFT得G：

[ G(0)=1+(-1)=0,\quad G(1)=1-(-1)=2 ]

奇数下标取出 (x(1)=0,x(3)=0)，做2点DFT得H：

[ H(0)=0,\quad H(1)=0 ]

第二步，做两级蝶形 (X(k)=G(k)+W_4^kH(k))，(X(k+2)=G(k)-W_4^kH(k))。旋转因子 (W_4^0=1)，(W_4^1=e^{-j\pi/2}=-j)。H全为0，所以：

[ X(0)=G(0)=0,\quad X(1)=G(1)=2 ] [ X(2)=G(0)=0,\quad X(3)=G(1)=2 ]

结果完全符合直觉：cos分量在正频率k=1和"负频率折返"的k=3处各贡献2。这个例子同时解释了为什么工程里看实信号频谱要关注N-k位置的镜像，也解释了为什么DFT幅度和N有关而不是直接在峰值显示信号的原始振幅，后面工程部分会再展开。

2.3 位序倒置是怎么冒出来的

N=4时没问题，N=8时你会发现一个现象：输入不是按0,1,2,3,4,5,6,7的顺序进入蝶形网络，而是按0,4,2,6,1,5,3,7。这就是位序倒置（bit-reversal），也是DIT初学最容易踩的坑。

为什么会出现倒序？因为每次递归都要把当前区间按奇偶下标重新排列。N=8先分成偶数区[0,2,4,6]和奇数区[1,3,5,7]；每一半继续按奇偶分，偶数区变成[0,4]和[2,6]，奇数区变成[1,5]和[3,7]。三次拆分后得到的顺序正好是0,4,2,6,1,5,3,7。把这个顺序写成二进制再看：原索引1（二进制001）变成了位置4（二进制100），原索引3（011）变成了位置6（110），本质上是把索引二进制位左右颠倒。

位序倒置不是bug，是"按时间抽取"的自然结果。你要么在进入蝶形前先把数组重排，要么把倒置计算放到蝶形里面，总归逃不掉。这个倒序规律对任意N=2^m都成立，所以实现FFT时第一步就是按位倒序重排输入数组。

2.4 信号流图级别看8点蝶形

用N=8把级数理一遍。总共 (\log_2 8=3) 级：

• 第1级：数组按位倒序后，相邻两个点一组形成4个2点蝶形，旋转因子全是 (W_8^0=1)，这一级只做加减
• 第2级：4个点一组形成2个块，每块里做2个蝶形，旋转因子分别是 (W_8^0) 和 (W_8^2)（等价于 (W_4^0) 和 (W_4^1)）
• 第3级：8个点一个大块，做4个蝶形，旋转因子是 (W_8^0,W_8^1,W_8^2,W_8^3)

第二级的旋转因子为什么可以写成 (W_8^2)，本质上还是可约性。不同教材画出的流图可能长得不一样，但点数、级数、旋转因子位置这三个要素是不变的。看懂这个结构后写代码就不会糊涂：每一级的块大小就是当前蝶形跨度len，块与块之间完全独立，可以顺序扫描。

3. 照着写就能跑：标准DIT基-2 FFT实现与逐段拆解

3.1 先说代码骨架，再说为什么这么搭

我用C语言写一个最小但完整的DIT基-2 FFT。代码只依赖标准库math.h，复数用结构体保存实部和虚部，便于初学者看清每一步。

#include <math.h> #define PI acos(-1.0) typedef struct { double re; double im; } Complex; /* 把低bits位的二进制数按位反转 */ static unsigned int reverse_bits(unsigned int x, int bits) { unsigned int y = 0; for (int i = 0; i < bits; i++) { y = (y << 1) | (x & 1); x >>= 1; } return y; } /* 输入数组重排为位倒序 */ void bit_reverse(Complex *x, int n) { int bits = 0; for (int t = n; t > 1; t >>= 1) bits++; for (int i = 0; i < n; i++) { unsigned int j = reverse_bits((unsigned int)i, bits); if (j > i) { Complex tmp = x[i]; x[i] = x[j]; x[j] = tmp; } } } /* 蝶形运算，n必须是2的整数次幂 */ void fft_dit(Complex *x, int n) { bit_reverse(x, n); for (int len = 2; len <= n; len <<= 1) { double theta = -2.0 * PI / len; Complex wn = {cos(theta), sin(theta)}; int half = len / 2; for (int i = 0; i < n; i += len) { Complex w = {1.0, 0.0}; for (int j = 0; j < half; j++) { Complex u = x[i + j]; /* v = x[i+half+j] * w */ Complex v = { x[i + j + half].re * w.re - x[i + j + half].im * w.im, x[i + j + half].re * w.im + x[i + j + half].im * w.re }; x[i + j].re = u.re + v.re; x[i + j].im = u.im + v.im; x[i + j + half].re = u.re - v.re; x[i + j + half].im = u.im - v.im; /* 下一个蝶形的旋转因子：w = w * wn */ Complex next = { w.re * wn.re - w.im * wn.im, w.re * wn.im + w.im * wn.re }; w = next; } } } }

3.2 三层循环的作用域，千万别搞反

很多初学FFT的人死在三层循环的顺序上。最外层的len是当前蝶形跨度，从2开始每轮翻倍，直到len=n。中间层是块起点i，每次跳len个点。最内层j遍历块内前一半的点，用旋转因子和对应后一半的点做蝶形。

为什么块步进要用len而不是别的值？因为第level级恰好有 (n/(2\cdot len)) 个块，每个块内部做half个蝶形，所以中间层i += len，内层j < half。很多人先算好所有旋转因子再用二维索引，没转过弯来，其实最稳的写法就是三层循环，每层负责一个维度：

• len：控制当前在分治的第几层
• i：定位块的起始位置
• j：块内第几个蝶形，同时决定用旋转因子表里的哪个值

旋转因子的更新逻辑也值得说明。内层j从0递增，w从1开始每轮乘一次wn，wn=W_len^1，于是第j个蝶形自然用的是W_len^j。这个写法虽然比"查表"稍慢，但逻辑最贴近蝶形流图，适合第一版调试。

3.3 两个测试样例，一分钟验证正确性

写完代码先别急着处理真实信号。我用两个边界测试验证：

第一个测试，输入冲激函数 (x(0)=1)，其余为0。冲激的DFT理论值应该全部是1。执行代码后，如果输出全为1（实部1、虚部0），说明整体结构正确。这是DFT最基础的"一致性测试"。

第二个测试，输入 (x[n]=\cos(2\pi\cdot2\cdot n/8))，也就是8点序列里恰好包含2个完整周期的余弦信号。执行N=8的FFT，理论结果应该是 X[2]=4、X[6]=4，其余接近0。为什么是4而不是2，前面N=4的推导已经解释过：单频余弦在正负频率各贡献一个幅度为N/2的峰，所以这里每个峰是8/2=4。

如果这两个测试都通过，FFT核心基本没问题。如果第一个测试通过、第二个不通过，问题一般出在旋转因子的符号或者输入序列的频率换算上；如果第一个测试就不通过，八成是位倒序没做对，或者蝶形里的加减写反了。这个排错顺序我用了很多次，比瞎看代码高效得多。

4. 工程落地之前：采样率、窗函数、归一化这几个坑不能跳过

4.1 频率分辨率和补零的真相

FFT的横轴不是"点数k"，而是频率 (f=k\cdot F_s/N)，其中Fs是采样率，N是FFT点数。相邻两条谱线之间的间隔叫频率分辨率：

[ \Delta f = \frac{F_s}{N} ]

这个公式意味着分辨率只取决于"采样率"和"FFT点数"，而FFT点数又受限于"实际信号时长"T=N/Fs。结论是：分辨率由你能拿到多长时间的信号决定。采样率44.1kHz、N=1024时，(\Delta f\approx43)Hz，10kHz的信号和10.043kHz的信号在频谱上根本分不开。

网上常有人问"我补零到8192为什么分辨率没变？"——补零只是在两个真实谱线之间插值，让曲线更光滑，并没有增加信号时长，所以不会提高真实分辨率。想要真正分辨两个频率相近的信号，必须延长采样时长。我做过一个振动监测项目，开始以为补零能"看得更清楚"，补到很大点数才发现峰位纹丝不动，后来把采样时间翻了一倍才真正把两个峰分开。

4.2 频谱泄漏与窗函数选择

FFT对有限长信号做周期延拓。如果你的N个采样点内部不是整数个周期，首尾接缝处就会产生跳变，频谱能量从信号真实频率"漏"到旁边一堆谱线上，这叫频谱泄漏。

矩形窗对应老老实实截断N个点，它的旁瓣峰值只比主瓣低13dB左右，泄漏效应相当明显。工程上通常先对时域数据乘一个窗函数，把首尾压小，再接缝处的不连续趋势平滑掉。常用窗函数的性能对比如下：

选窗的本质是权衡：主瓣越窄，频率分辨能力越强；旁瓣越低，泄漏抑制越好。没有任何窗同时做到两者最优。我日常做音频频谱可视化，默认先用Hann窗；如果做的是单频信号的幅度校准，会用Blackman-Harris，因为它旁瓣极其干净，读幅度误差更小。

4.3 归一化：为什么峰值不直接等于信号幅度

很多第一次用FFT的人会问："我输入幅度1的余弦信号，频谱图里峰值为什么不是1？"前面已经提到，DFT的幅度是叠加性的。实信号N点FFT，单频cos在正负频率各产生N/2的峰。要还原真实幅度，工程做法分两步：

第一步，除以N，把DFT的"叠加和"变成"平均幅度"。 第二步，改成单边谱时，把正频率部分的幅度乘2，直流分量（k=0）和奈奎斯特点（k=N/2）不乘2。

组合起来就是：对实信号的正频率部分做 (2|X(k)|/N)，直流部分做 (|X(0)|/N)。这样输入幅度1的余弦信号，在对应频率谱线上读出来的就是1。

加窗之后还要再乘一个能量恢复系数。因为窗函数会压低时域信号的总体能量，导致频谱幅度整体变小。经验做法是：

[ \text{恢复系数}=\frac{N}{\sum_{n=0}^{N-1}w(n)} ]

Hann窗的系数大约是2，Blackman-Harris窗的系数约1.7左右。具体数值随N略有浮动，我一般在代码里先算一次窗总和，再乘恢复系数，严谨且普适。

5. 从基-2往两边看：DIF、分裂基与实序列优化的取舍思路

5.1 和DIT对应的DIF，到底差在哪

DIT按时间抽取，是按输入下标n的奇偶拆分；另一大类DIF（Decimation In Frequency）按频率抽取，是按输出下标k的奇偶拆分。DIF的蝶形结构是"先加减、后乘旋转因子"，DIT是"先乘旋转因子、后加减"。两者的计算量是同一个量级，但数据流方向刚好相反：

工程中很多现成库用DIF，是因为它可以自然序输入，位倒置放在最后，某些流水线硬件实现里更省事。但我建议初学者先学DIT，理由只有一个：DIT的分解过程更直观，从信号本身出发，一步拆成奇偶两组，理解成本最低。DIF看懂了DIT之后，半小时就能改出来。

5.2 基-2不是终点：实序列优化与分裂基

如果你只处理实信号，直接用N点复FFT有点浪费。实序列的DFT有共轭对称性，可以利用对称关系，把N点实序列拼成一个N/2点复序列，跑一次N/2点复FFT，再把结果拆成N点实序列频谱，计算量大约减半。这个方法在频谱仪这类高频测量工具里很常见。

另一个方向是分裂基（split-radix）。它不是纯基-2，而是把偶数下标用基-2、奇数下标用基-4混合拆分，在基-2基础上再省约三分之一的复数乘法。如果你看过一些FFT库的源码，会发现它们内部的蝶形种类比教科书里复杂，常见原因就是用了混合基。理解标准DIT基-2之后，读混合基代码也能看出脉络，无非是多算几种子蝶形的旋转因子组合。

5.3 什么时候该自己写，什么时候直接调库

自己手写FFT的价值在理解，不在性能。工程上，我建议绝大多数场景直接用成熟的开源FFT库，它们做了Cache优化、SIMD向量化，甚至对不同N自动切换拆分策略，性能比自己第一版代码高好几倍。但前提是你要知道库用的是正负号约定还是反约定、归一化在内部做没做、输入是自然序还是倒序。这些信息不手写一遍FFT很难真正理解，而一旦理解，换个库、换种配置就只是查文档的事。

我自己保留手写版本的场景是教学演示和嵌入式小平台。在某些单片机环境里，外部库过于臃肿，或不能引入浮点SIMD，那么一份两三页的DIT基-2实现稳定可靠、方便剪裁，反而是最省事的方案。

附：我踩过几次坑之后，最想留下的排查清单

最后按经验排序留一批检查项。FFT结果不对，90%的情况逃不出这几条：

• 位倒序忘做或做错。先用冲激输入x(0)=1验证输出是否全1，这个测试能立刻暴露倒序问题
• 旋转因子符号搞反。DFT指数是负号，对应theta=-2π/len。如果用了+2π/len，频率轴会反转
• 层的步进写错。中间层i必须i+=len，不是i+=2或i+=half
• 归一化时机不对。FFT内部很少除以N，归一化要放在后续幅度计算里，否则会误判库的输出
• 窗函数和幅度恢复系数不匹配。换了窗却没换系数，测出来的幅度就会系统性偏移

这套代码现在我还在一个音频分析工具里用。STFT按帧处理，256点一帧配合Hann窗，在普通电脑上能做到实时刷新。写完FFT后的那种感觉是：算法本身并不复杂，复杂的是那些藏在公式背后的工程细节。希望这篇从原理到实现的拆解，能让你少走一些弯路。