前面三篇笔记把频域滤波的主线推到了DFT这一步信号变到频域乘上滤波器传递函数再变回时域。思路很顺可实际操作时有个过不去的坎——DFT按定义直接算计算量是O(N²)。我印象很深有次做一个采集分析Demo8K采样率、每次处理4096点按DFT定义硬算一帧频谱开发机都要明显卡一下换成性能弱一些的板子帧间隔根本压不住做实时的念头直接泡汤。后来老老实实换成FFT同样帧长瞬间出结果差距完全不在一个量级。这篇笔记就专门拆解FFT里最常用的一支基2时间抽取FFTDIT-FFT。它不依赖任何库从推导到代码都能自己写完也是后续做频域滤波正反变换时绕不开的内核。适合刚把DFT公式背熟、想知道“计算机里到底怎么算傅里叶变换”的读者也适合想手写FFT用于移植或教学的人。1. DFT直接算到底贵在哪从一次实测说起1.1 一组数据感受一下N²的量级DFT的定义式长这样X(k) Σ_{n0}^{N-1} x(n) * W_N^{nk}其中 W_N e^{-j2π/N}k 取 0 到 N-1。直接按这个公式实现每个 X(k) 都要做 N 次复乘一共 N 个频点所以总复乘数是 N²。N8 的时候 64 次无所谓N256 的时候 65536 次也还能忍。但到了 N1024就是一百多万次复乘N8192差不多六千七百万次N16384直接奔着两亿七千万次去了。两亿七千万次复乘听起来是个数字放到嵌入式环境里就是秒级的处理时间。而你做频域滤波时FFT 只是链条里的一环前面要加窗后面要频域相乘、再做 IFFT任何一个环节卡住整条链路都跑不起来。这也是为什么 Cooley 和 Tukey 在 1965 年重新整理出 FFT 算法后数字信号处理才算真正进入工程时代——复杂度从 O(N²) 降到 O(N log₂N)N16384 时大约只需要 114688 次复乘差了接近两千三百倍。1.2 旋转因子的三个性质周期、对称、可约FFT 能省算力核心不是发明了什么新运算而是把 DFT 里大量重复的复数乘法合并掉了。重复的根源在旋转因子 W_N^{nk} 身上它有三个教科书一定会写、但初学者不太当回事的性质周期性W_N^{kN} W_N^k意思是旋转因子按 N 循环顺序转一圈又回到原点。对称性W_N^{kN/2} -W_N^k等价于“转了半圈就是相反方向”。可约性W_{aN}^{ak} W_N^k上下同乘一个倍数指数可以约分。这三个性质看起来平淡实际上蝴蝶效应的起点。比如 W_8^2 W_4^1W_16^4 W_4^1同一个旋转因子可以被许多频点复用。朴素的 DFT 把每个 X(k) 当独立任务各算各的FFT 则是把共用的旋转因子提取出来用分治法的骨架把任务层层拆小让一个乘法结果能被多个输出复用。1.3 “合并同类项”这个方向是怎么来的如果你自己动手推导过 N4 或者 N8 的 DFT会注意到很多 W_N^{nk} 的值是一样的甚至有些因为对称性可以加个负号直接复用。我最早手算 N8 的 DFT 时就发现某些中间乘法在多个频点里重复出现当时只是觉得“这有点浪费”还没想到能做成通用算法。后来才理解FFT 的路线就是把这种“浪费”系统化先把问题拆成子问题子问题里的公共计算只做一次再通过蝶形运算把结果合并起来。2. 基2时间抽取的完整推导一分为二、二分为四2.1 为什么按奇偶分而不是按前后分“时间抽取”这个名字的来由就是把时域序列按照序号 n 的奇偶性拆成两组。有人会问为什么不按前半段、后半段分答案是按前后分也能拆出递归结构但旋转因子没法干净地约分成 N/2 点 DFT按奇偶分x(2r) 和 x(2r1) 天然让旋转因子的指数变成 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) Σ_{r0}^{N/2-1} x(2r) W_N^{2rk} Σ_{r0}^{N/2-1} x(2r1) W_N^{(2r1)k}第一项里 W_N^{2rk} W_{N/2}^{rk}第二项里 W_N^{(2r1)k} W_N^k * W_{N/2}^{rk}。于是X(k) G(k) W_N^k * H(k)其中 G(k) Σ_{r0}^{N/2-1} x(2r) W_{N/2}^{rk}是偶数点子的 N/2 点 DFTH(k) Σ_{r0}^{N/2-1} x(2r1) W_{N/2}^{rk}是奇数点子的 N/2 点 DFT。2.3 后半段频点X(kN/2) 的四行推导现在关键来了。G(k) 和 H(k) 只有 N/2 个不同值k 从 0 到 N/2-1但 X(k) 需要 N 个输出。剩下的一半怎么来把 k 换成 kN/2 代入X(kN/2) G(kN/2) W_N^{kN/2} H(kN/2)因为 G 和 H 都是 N/2 点 DFT它们满足周期性G(kN/2) G(k)H(kN/2) H(k)。再配合旋转因子的对称性 W_N^{kN/2} -W_N^k于是X(kN/2) G(k) - W_N^k * H(k)这一个加号、一个减号就是蝶形运算的原型。把 G(k) 和 W_N^k H(k) 当作两个输入一次蝶形同时产出两个输出X(k) 和 X(kN/2)。整个 FFT 就是不断重复这个“两输入、两输出”的结构。2.4 递归到底2点DFT就是一次蝶形子问题 N/2 点 DFT 还能继续按奇偶拆分拆到 N/21 为止。但实践中拆到 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 WbB a - Wb。一次蝶形只做一次复乘、两次复加。所谓“就地计算”就是算完后把 A 存回 a 原来的位置B 存回 b 原来的位置后续级数不再需要旧值。我用生活化的比喻DFT 是整个合唱团按乐谱各唱各的FFT 是一群人对拍子每一轮两两一组对一下对完就同步一部分信息下一轮再扩大组合。蝶形就是最小的“对拍子”动作。3.2 旋转因子在每一级怎么摆N8 时一共 3 级每级 4 个蝶形。旋转因子的位置不是乱放的规律如下级数 L蝶形组数每组蝶形数蝶距旋转因子1411W_8^02222W_8^0W_8^23144W_8^0W_8^1W_8^2W_8^3一般规律第 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^1W_8^6 之类其实都被可约性消化掉了代码里可以直接用当前级的 step 来计算。3.3 手工走一遍N8用x[n]n验证我拿一个最简单的信号 x[n] nn0..7完整走一遍方便你对着验证。先做位反转排序下一节细讲输入顺序变为 [0,4,2,6,1,5,3,7]。第 1 级蝶距 1两两相加相减W1。第 2 级蝶距 2两组蝶形旋转因子用 1 和 -j。中间我算到数组变成[12, -44j, -4, -4-4j, 16, -44j, -4, -4-4j]第 3 级蝶距 4旋转因子分别是 1、W_8^1、-j、W_8^3其中 W_8^1 ≈ 0.7071-0.7071jW_8^3 ≈ -0.7071-0.7071j。这一步会得到最终结果。k实部虚部02801-49.65682-443-41.65684-405-4-1.65686-4-47-4-9.6568这个结果和标准FFT库算出来完全一致。我建议你拿到代码后第一件事就跑这个用例只要这 8 个复数对得上说明位反转、旋转因子方向、蝶形顺序全部正确后面再换复杂信号都稳。3.4 计算量复核N8 时直接 DFT 需要 64 次复乘DIT-FFT 每级 4 次蝶形、共 3 级总复乘 12 次节省 5 倍多。N1024 时节省接近 200 倍这个差距会随 N 增大越来越夸张。理解这个数字你就明白为什么频域滤波的链条里FFT 是绕不开的底桩。4. 位反转排序最容易出错且藏得最深的环节4.1 为什么会有位反转时间抽取要求先把序列按“二进制最低位”分组偶数放前半、奇数放后半下一层再按次低位分组继续递归。这个过程如果直接体现在数据位置上最终输入数组的排列顺序就是把每个下标按二进制倒过来。比如下标 1001会被换到 4100下标 3011换到 6110。这个现象我上第一次课的时候没反应过来后来自己画 N8 的递归树才看明白每一层按一个二进制位分左右所有层叠起来最终位置恰好是原下标的比特反转。4.2 N8的位反转对照表原下标二进制反转后实际下标0000000010011004201001023011110641000011510110156110011371111117所以输入 [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 再反转”的进位过程。建议你逐行走一遍 N8 的情况走完就彻底记住它输出序列是 0,4,2,6,1,5,3,7。我第一次用这段代码时没走查证直接套到了 N16 上结果频谱乱了半天最后才发现是位反转写错了一位掩码。4.4 如果不想显式做位反转有一个偷懒的办法用递归实现 FFT。递归调用天然把偶数序列和奇数序列分到两个子数组数据位置不需要预先重排递归返回后再合并。代价是要不断分配子数组内存操作多性能不如迭代版。我的建议是学习阶段用递归理解原理工程实现还是用迭代版加位反转速度更稳。5. 手写DIT-FFTC与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 常见坑列表现象根因输出频谱左右颠倒旋转因子指数符号写反结果等于 0 或错乱位反转没做对输入顺序不对幅值差 N 倍IFFT 忘记除以 N长度为非 2 的幂时崩溃代码假设基2需要先补零长序列结果有微小“毛刺”浮点累计误差正常现象6. 旋转因子表、归一化和其他绕不开的工程细节6.1 预计算旋转因子表性能翻倍的做法上面的代码里每个蝶形都现场调用 cos/sin这在 PC 上问题不大但在嵌入式或数据帧率很高的场景三角函数开销不可忽略。工程做法是提前把旋转因子存成一张表长度为 N/2第 k 项存 W_N^k。这样蝶形内部只需要查表省去大量浮点函数调用。我自己的经验N4096 时预计算 VS 现场计算的耗时差距在普通开发机上大概是 20%-30%帧率越高越明显。代价是内存多占用 N/2 个 float 复数对N 不超过几万的话完全值得。6.2 归一化与频谱顺序FFT 结果不自动归一化这也是初学者最懵的地方。比如 x[n]1 的 8 点 FFTX[0]8不是 1。做频域滤波时通常把滤波器传递函数和 FFT 结果相乘再到 IFFT 时除以 N所以真正严格的正逆变换要配成一对FFT 不除IFFT 除整体就还原了。频谱顺序还有个细节FFT 输出下标 0 是直流1 到 N/2-1 是正频率N/2 是奈奎斯特频率N/21 到 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 和它与时间抽取的对称性到时候你会发现蝶形的内核其实都是同一件事。