手写基2时间抽取FFT（DIT-FFT）：从DFT到蝶形运算与位反转

前面三篇笔记把频域滤波的主线推到了DFT这一步：信号变到频域，乘上滤波器传递函数，再变回时域。思路很顺，可实际操作时有个过不去的坎——DFT按定义直接算，计算量是O(N²)。我印象很深，有次做一个采集分析Demo，8K采样率、每次处理4096点，按DFT定义硬算一帧频谱，开发机都要明显卡一下；换成性能弱一些的板子，帧间隔根本压不住，做实时的念头直接泡汤。后来老老实实换成FFT，同样帧长瞬间出结果，差距完全不在一个量级。

这篇笔记就专门拆解FFT里最常用的一支：基2时间抽取FFT（DIT-FFT）。它不依赖任何库，从推导到代码都能自己写完，也是后续做频域滤波正反变换时绕不开的内核。适合刚把DFT公式背熟、想知道“计算机里到底怎么算傅里叶变换”的读者，也适合想手写FFT用于移植或教学的人。

1. DFT直接算到底贵在哪：从一次实测说起

1.1 一组数据感受一下N²的量级

DFT的定义式长这样：X(k) = Σ_{n=0}^{N-1} x(n) * W_N^{nk}，其中 W_N = e^{-j2π/N}，k 取 0 到 N-1。

直接按这个公式实现，每个 X(k) 都要做 N 次复乘，一共 N 个频点，所以总复乘数是 N²。N=8 的时候 64 次，无所谓；N=256 的时候 65536 次，也还能忍。但到了 N=1024，就是一百多万次复乘；N=8192，差不多六千七百万次；N=16384，直接奔着两亿七千万次去了。

两亿七千万次复乘听起来是个数字，放到嵌入式环境里就是秒级的处理时间。而你做频域滤波时，FFT 只是链条里的一环，前面要加窗，后面要频域相乘、再做 IFFT，任何一个环节卡住，整条链路都跑不起来。这也是为什么 Cooley 和 Tukey 在 1965 年重新整理出 FFT 算法后，数字信号处理才算真正进入工程时代——复杂度从 O(N²) 降到 O(N log₂N)，N=16384 时大约只需要 114688 次复乘，差了接近两千三百倍。

1.2 旋转因子的三个性质：周期、对称、可约

FFT 能省算力，核心不是发明了什么新运算，而是把 DFT 里大量重复的复数乘法合并掉了。重复的根源在旋转因子 W_N^{nk} 身上，它有三个教科书一定会写、但初学者不太当回事的性质：

• 周期性：W_N^{k+N} = W_N^k，意思是旋转因子按 N 循环，顺序转一圈又回到原点。
• 对称性：W_N^{k+N/2} = -W_N^k，等价于“转了半圈就是相反方向”。
• 可约性：W_{aN}^{ak} = W_N^k，上下同乘一个倍数，指数可以约分。

这三个性质看起来平淡，实际上蝴蝶效应的起点。比如 W_8^2 = W_4^1，W_16^4 = W_4^1，同一个旋转因子可以被许多频点复用。朴素的 DFT 把每个 X(k) 当独立任务，各算各的；FFT 则是把共用的旋转因子提取出来，用分治法的骨架把任务层层拆小，让一个乘法结果能被多个输出复用。

1.3 “合并同类项”这个方向是怎么来的

如果你自己动手推导过 N=4 或者 N=8 的 DFT，会注意到很多 W_N^{nk} 的值是一样的，甚至有些因为对称性可以加个负号直接复用。我最早手算 N=8 的 DFT 时，就发现某些中间乘法在多个频点里重复出现，当时只是觉得“这有点浪费”，还没想到能做成通用算法。后来才理解，FFT 的路线就是把这种“浪费”系统化：先把问题拆成子问题，子问题里的公共计算只做一次，再通过蝶形运算把结果合并起来。

2. 基2时间抽取的完整推导：一分为二、二分为四

2.1 为什么按奇偶分，而不是按前后分

“时间抽取”这个名字的来由就是：把时域序列按照序号 n 的奇偶性拆成两组。有人会问，为什么不按前半段、后半段分？答案是：按前后分也能拆出递归结构，但旋转因子没法干净地约分成 N/2 点 DFT；按奇偶分，x(2r) 和 x(2r+1) 天然让旋转因子的指数变成 2nk，正好套用可约性 W_N^{2nk} = W_{N/2}^{nk}，子问题就变成两个纯粹的 N/2 点 DFT。这一步是整个算法结构的支点。

2.2 三步展开：奇偶分列、代入约分、整理成蝶形

设 N 为 2 的幂，记 W_N = e^{-j2π/N}。把 x(n) 按 n 的奇偶拆开：

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}

第一项里 W_N^{2rk} = W_{N/2}^{rk}，第二项里 W_N^{(2r+1)k} = W_N^k * W_{N/2}^{rk}。于是：

X(k) = G(k) + W_N^k * H(k)

其中 G(k) = Σ_{r=0}^{N/2-1} x(2r) W_{N/2}^{rk}，是偶数点子的 N/2 点 DFT；H(k) = Σ_{r=0}^{N/2-1} x(2r+1) W_{N/2}^{rk}，是奇数点子的 N/2 点 DFT。

2.3 后半段频点：X(k+N/2) 的四行推导

现在关键来了。G(k) 和 H(k) 只有 N/2 个不同值（k 从 0 到 N/2-1），但 X(k) 需要 N 个输出。剩下的一半怎么来？把 k 换成 k+N/2 代入：

X(k+N/2) = G(k+N/2) + W_N^{k+N/2} H(k+N/2)

因为 G 和 H 都是 N/2 点 DFT，它们满足周期性：G(k+N/2) = G(k)，H(k+N/2) = H(k)。再配合旋转因子的对称性 W_N^{k+N/2} = -W_N^k，于是：

X(k+N/2) = G(k) - W_N^k * H(k)

这一个加号、一个减号，就是蝶形运算的原型。把 G(k) 和 W_N^k H(k) 当作两个输入，一次蝶形同时产出两个输出：X(k) 和 X(k+N/2)。整个 FFT 就是不断重复这个“两输入、两输出”的结构。

2.4 递归到底：2点DFT就是一次蝶形

子问题 N/2 点 DFT 还能继续按奇偶拆分，拆到 N/2=1 为止。但实践中拆到 2 点就够了，因为 2 点 DFT 本身就是蝶形：X(0) = x(0)+x(1)，X(1) = x(0)-x(1)，正好对应 W_2^0 = 1 的蝶形。所以基2 FFT 的“基”指的就是每层只按二进制的一位拆分，一直分到 2 为止。

从复杂度看，每分一层，子问题变成两个半尺寸 DFT，加上 N/2 次蝶形合并。递归关系式 T(N) = 2T(N/2) + N/2，解出来正好是 N/2 * log₂N 次复乘。

3. 蝶形运算与8点FFT的手工推演

3.1 蝶形单元的输入输出关系

蝶形运算的接口很固定：两个复数输入 a、b，一个旋转因子 W，输出 A = a + Wb，B = a - Wb。一次蝶形只做一次复乘、两次复加。所谓“就地计算”，就是算完后把 A 存回 a 原来的位置，B 存回 b 原来的位置，后续级数不再需要旧值。

我用生活化的比喻：DFT 是整个合唱团按乐谱各唱各的；FFT 是一群人对拍子，每一轮两两一组对一下，对完就同步一部分信息，下一轮再扩大组合。蝶形就是最小的“对拍子”动作。

3.2 旋转因子在每一级怎么摆

N=8 时一共 3 级，每级 4 个蝶形。旋转因子的位置不是乱放的，规律如下：

一般规律：第 L 级，蝶距为 2^{L-1}，组数为 N / 2^L，第 j 个蝶形（j 从 0 到 2^{L-1}-1）使用旋转因子 W_N^{j*N/2^L}。注意 W_8^2 = W_4^1，W_8^6 之类其实都被可约性消化掉了，代码里可以直接用当前级的 step 来计算。

3.3 手工走一遍N=8：用x[n]=n验证

我拿一个最简单的信号 x[n] = n（n=0..7）完整走一遍，方便你对着验证。先做位反转排序（下一节细讲），输入顺序变为 [0,4,2,6,1,5,3,7]。

第 1 级：蝶距 1，两两相加相减（W=1）。

第 2 级：蝶距 2，两组蝶形，旋转因子用 1 和 -j。中间我算到数组变成：

[12, -4+4j, -4, -4-4j, 16, -4+4j, -4, -4-4j]

第 3 级：蝶距 4，旋转因子分别是 1、W_8^1、-j、W_8^3，其中 W_8^1 ≈ 0.7071-0.7071j，W_8^3 ≈ -0.7071-0.7071j。这一步会得到最终结果。

这个结果和标准FFT库算出来完全一致。我建议你拿到代码后第一件事就跑这个用例，只要这 8 个复数对得上，说明位反转、旋转因子方向、蝶形顺序全部正确，后面再换复杂信号都稳。

3.4 计算量复核

N=8 时，直接 DFT 需要 64 次复乘；DIT-FFT 每级 4 次蝶形、共 3 级，总复乘 12 次，节省 5 倍多。N=1024 时节省接近 200 倍，这个差距会随 N 增大越来越夸张。理解这个数字，你就明白为什么频域滤波的链条里，FFT 是绕不开的底桩。

4. 位反转排序：最容易出错且藏得最深的环节

4.1 为什么会有位反转

时间抽取要求先把序列按“二进制最低位”分组：偶数放前半、奇数放后半，下一层再按次低位分组，继续递归。这个过程如果直接体现在数据位置上，最终输入数组的排列顺序就是把每个下标按二进制倒过来。比如下标 1（001）会被换到 4（100），下标 3（011）换到 6（110）。

这个现象我上第一次课的时候没反应过来，后来自己画 N=8 的递归树才看明白：每一层按一个二进制位分左右，所有层叠起来，最终位置恰好是原下标的比特反转。

4.2 N=8的位反转对照表

所以输入 [0,1,2,3,4,5,6,7] 排完位反转是 [0,4,2,6,1,5,3,7]。这个顺序我上面手工推演里直接用了，现在你明白它是怎么来的了。

4.3 用迭代算法生成位反转序号

一个很简单的方法：对每个 i 单独“求反转数再交换”，但那样效率不高。标准工程写法是用一个迭代器在循环中推导下一个反转序号，不需要查表：

int j = 0; for (int i = 0; i < N - 1; i++) { if (i < j) { swap(x[i], x[j]); } int k = N >> 1; while (j >= k) { j -= k; k >>= 1; } j += k; }

这个想法的本质是：从当前 j 出发，模拟“二进制加 1 再反转”的进位过程。建议你逐行走一遍 N=8 的情况，走完就彻底记住：它输出序列是 0,4,2,6,1,5,3,7。我第一次用这段代码时没走查证，直接套到了 N=16 上，结果频谱乱了半天，最后才发现是位反转写错了一位掩码。

4.4 如果不想显式做位反转

有一个偷懒的办法：用递归实现 FFT。递归调用天然把偶数序列和奇数序列分到两个子数组，数据位置不需要预先重排，递归返回后再合并。代价是要不断分配子数组，内存操作多，性能不如迭代版。我的建议是学习阶段用递归理解原理，工程实现还是用迭代版加位反转，速度更稳。

5. 手写DIT-FFT：C与Python双版本实现及验证

5.1 C语言版本：就地运算

#include <math.h> #include <complex.h> void dit_fft(complex float x[], int n) { // 参数n必须为2的幂 int i, j, k, L, step, half, t, m; // 位反转 j = 0; for (i = 0; i < n - 1; i++) { if (i < j) { complex float tmp = x[i]; x[i] = x[j]; x[j] = tmp; } k = n >> 1; while (j >= k) { j -= k; k >>= 1; } j += k; } // 蝶形运算 for (L = 1; (1 << L) <= n; L++) { step = 1 << L; half = step >> 1; for (t = 0; t < n; t += step) { for (m = 0; m < half; m++) { double angle = -2.0 * M_PI * m / step; complex float w = cos(angle) + I * sin(angle); int a = t + m; int b = t + m + half; complex float u = x[b] * w; x[b] = x[a] - u; x[a] = x[a] + u; } } } }

C 版本直接用 C99 的 complex float，代码很短。需要注意旋转因子角度我写成-2π*m/step，对应教科书里的 W_N^{m*N/step}，符号约定别搞反，否则频谱会沿中心翻转。

5.2 Python版本：便于验证和画图

import cmath def bit_reverse(x): n = len(x) j = 0 for i in range(1, n): bit = n >> 1 while j & bit: j ^= bit bit >>= 1 j |= bit if i < j: x[i], x[j] = x[j], x[i] def dit_fft(x): n = len(x) bit_reverse(x) levels = 0 while (1 << levels) < n: levels += 1 for L in range(1, levels + 1): step = 1 << L half = step // 2 for t in range(0, n, step): for k in range(half): w = cmath.exp(-2j * cmath.pi * k / step) a = t + k b = t + k + half u = x[b] * w x[b] = x[a] - u x[a] = x[a] + u x = [0 + 0j, 1 + 0j, 2 + 0j, 3 + 0j, 4 + 0j, 5 + 0j, 6 + 0j, 7 + 0j] dit_fft(x) print(x)

这段代码跑出来的结果，和前文表格里的 8 个复数值完全一致。我的建议是先用 Python 把每个中间数组打印出来，跟手算过程逐级对照，把蝶形顺序彻底弄清，再动手写 C 版，这样排错成本低很多。

5.3 IFFT怎么复用同一套蝶形

频域滤波最终要把结果变回时域，所以 IFFT 也得有。一个经典技巧：IFFT(X) = conj(FFT(conj(X))) / N。也就是说，把输入先取共轭，跑一遍 FFT，再取共轭、除以 N。这个方法对 DIT 结构完全适用，不需要改任何蝶形代码。

我在实际写频域滤波 Demo 时就用的这招，一份 FFT 函数包含正反变换，省事且不容易引入第二套旋转因子 bug。唯一要记住的是 IFFT 最后那个除以 N，很多刚入门的同学会漏掉，滤波输出幅度直接差 N 倍。

5.4 常见坑列表

6. 旋转因子表、归一化和其他绕不开的工程细节

6.1 预计算旋转因子表：性能翻倍的做法

上面的代码里每个蝶形都现场调用 cos/sin，这在 PC 上问题不大，但在嵌入式或数据帧率很高的场景，三角函数开销不可忽略。工程做法是提前把旋转因子存成一张表：长度为 N/2，第 k 项存 W_N^k。这样蝶形内部只需要查表，省去大量浮点函数调用。

我自己的经验：N=4096 时，预计算 VS 现场计算的耗时差距在普通开发机上大概是 20%-30%，帧率越高越明显。代价是内存多占用 N/2 个 float 复数对，N 不超过几万的话完全值得。

6.2 归一化与频谱顺序

FFT 结果不自动归一化，这也是初学者最懵的地方。比如 x[n]=1 的 8 点 FFT，X[0]=8，不是 1。做频域滤波时通常把滤波器传递函数和 FFT 结果相乘，再到 IFFT 时除以 N，所以真正严格的正逆变换要配成一对：FFT 不除，IFFT 除，整体就还原了。

频谱顺序还有个细节：FFT 输出下标 0 是直流，1 到 N/2-1 是正频率，N/2 是奈奎斯特频率，N/2+1 到 N-1 对应负频率。频域滤波如果只想保留某个频率区段，要搞清楚负频率部分的位置，否则滤完的波形会不对称。这点在前几篇笔记里提过，这里再强调一次：先 fftshift 理清频率轴，再设计滤波器响应，能省一整轮调试。

6.3 非2的幂长度：补零或换算法

基2 FFT 最吸引人，但现实信号帧长不一定恰好是 2 的幂。两个选择：一是补零到最近 2 的幂，比如 1500 点补到 2048 点，简单好用，但会引入频率分辨率变化；二是改用混合基 FFT，比如 3、5 等质因数分解，都能基于类似的分治思想实现。我建议绝大多数场景先补零，代码简单、足够可靠。如果业务对分辨率和算力都敏感，再研究混合基。

补零还有个副作用是频谱形状会稍微变“圆滑”，但那只是插值效果，不是信息变多。分析时心里有数即可。

6.4 FFT/IFFT在频域滤波链路里的完整配合

回到这篇笔记的起点——频域滤波。一条完整的链路是：时域分段加窗 → FFT → 频域乘传递函数 → IFFT → 重叠相加或去窗恢复 → 输出时域。其中 FFT 和 IFFT 一共要跑两次，所以我建议把正反变换封装成同一个接口，加一个inverse标志：

• forward：只做 FFT，不缩放；
• inverse：共轭技巧 + FFT + 除以 N。

我在做 Demo 时就是这么封的，好处是后续如果要换混合基、加窗、加帧重叠，调用层完全不用动。

6.5 最后一个我个人的小建议

调试这类算法，最忌讳上来就喂真实信号。真实信号混合了噪声、偏置、多频分量，FFT 输出一乱根本不知道哪一步错。我的固定套路是先用 x[n]=n、x[n]=sin(2π*n/N) 这类解析解已知的信号做回归测试，确认 8 点、16 点、1024 点结果全对，再上真实数据。这套流程虽然朴素，但帮我避过好几次位反转和旋转因子符号的坑。

基2时间抽取 FFT 看似公式推导只有几行，真正让它跑起来的是蝶形顺序、位反转、旋转因子方向这些细节。把这一块吃透，后面的频域滤波、加窗、重叠保留法都顺了。下一篇笔记可以接着写频率抽取 FFT 和它与时间抽取的对称性，到时候你会发现，蝶形的内核其实都是同一件事。