FFT快速傅里叶变换Fast Fourier Transform是数字信号处理中最常用的算法之一。无论是振动信号频谱分析、音频处理、通信系统还是在 STM32F407 这类单片机上做实时频谱显示FFT 都是把时域波形转换为频域信息的核心步骤。很多开发者可以直接调用现成的 FFT 函数却不太理解内部蝶形运算结构为什么输入要先做二进制位倒序旋转因子按什么规律更新输出结果的顺序到底是怎样的本文围绕 FFT 快速傅里叶变换的蝶形运算结构方法展开从 DFT 的计算量谈起逐步拆解基 2 时间抽取 FFT 的蝶形图用 C 语言实现一个最小可运行程序并给出验证、排错和工程落地建议。读完这篇文章FFT 不再是一个只能接受的黑盒。1. FFT为什么要把DFT拆成蝶形结构1.1 DFT的直接计算量朴素算法难以应对实时场景离散傅里叶变换DFT的公式是X(k) Σ_{n0}^{N-1} x(n) · W_N^{kn}其中 W_N e^{-j2π/N}k 0, 1, ..., N-1按这个公式直接计算每个 X(k) 都需要 N 次复数乘法和 N-1 次复数加法计算全部 N 个频点总共需要大约 N² 次复数乘法。N1024 时直接计算约需要 1048576 次复数乘法N8192 时这个数字膨胀到约 67108864。在实时信号处理系统里这是一个几乎不可能接受的量级。FFT 的突破口在于利用旋转因子 W_N^{kn} 的对称性和周期性把一个大 N 点 DFT 逐级拆成多个小 N 点 DFT最终把复杂度降到 O(N log₂N)。同样是 1024 点FFT 大约只需要 1024 × 10 10240 次蝶形级联操作比直接计算少了两个数量级。点数 N直接 DFT 复数乘法量级FFT 复数乘法量级差距倍数1625664412816384896约 181024约 104 万约 10240约 1028192约 6711 万约 106496约 6301.2 旋转因子的两个性质对称性和周期性让 DFT 可以拆分的关键是旋转因子 W_N e^{-j2π/N} 的特殊性质。对称性W_N^{k N/2} -W_N^k周期性W_N^{k N} W_N^k这两个性质意味着在计算一个较大 N 点 DFT 时很多中间项的旋转因子是相同或相反的。通过重新组织计算顺序相同的乘法可以只算一次相反的两个结果可以合并到一个加减操作里。这就是“蝶形运算”能够压缩计算量的数学基础。基 2 时间抽取DIT是最常见的拆分方式。它把输入序列按偶数下标和奇数下标分成两组分别计算两个 N/2 点 DFT再用旋转因子 W_N^k 把两个半序列的结果合并。这样一次拆分就能让计算量减半继续递归拆分最终得到由 log₂N 级蝶形构成的运算结构。1.3 蝶形运算一次两输入两输出的基本计算单元蝶形运算是 FFT 的最小处理单元。它接收两个复数输入 A 和 B输出两个复数结果X A W · B Y A - W · B其中 W 是复数旋转因子。从信号流图上看输入两个节点输出两个节点中间交叉形状像蝴蝶所以叫蝶形运算。一次蝶形只包含一次复数乘法和两次复数加减法。如果不做拆分A W·B 和 A - W·B 会因为 W 的符号关系被重复计算两次蝶形结构把共用的 W·B 提前算一次两个输出共用从而节省了一半乘法。学习 FFT 时最容易出现的一个误解是认为蝶形结构只是某种加速技巧。实际上蝶形结构就是 DFT 数学性质在计算顺序上的直接映射理解它之后才能解释位倒序、旋转因子更新规律和输出排列顺序这些工程细节。2. 从8点例子看蝶形运算的三级结构2.1 时间抽取之后输入为什么要做二进制位倒序基 2 时间抽取 FFT 的第一步是把输入序列按下标奇偶分成两组。每次分组都会让序列的顺序发生变化。N8 时分组过程如下第 1 次分组偶数下标 x(0), x(2), x(4), x(6)奇数下标 x(1), x(3), x(5), x(7)。继续对每个组再按奇偶分组最终得到的排列顺序是x(0), x(4), x(2), x(6), x(1), x(5), x(3), x(7)这个顺序恰好等于原始下标 0 到 7 的二进制位倒序结果。把每个下标写成 3 位二进制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 时第一步通常是把输入序列按位倒序索引重新排列。这一步不是可选的它保证了后续蝶形级联时配对的两个数据正好落在同一个蝶形的两个输入位置上。2.2 第 1 级蝶形跨度 1旋转因子只有 W_2^0N8 时第一级蝶形处理块长为 2共有 4 个蝶形组。每个组内只有一个蝶形(0,1), (2,3), (4,5), (6,7)这一级对应的旋转因子是 W_2^0 1。因为 k 只取 0所以这一级实际上不产生复数乘法只做复数加减法X0 x0 x1 X1 x0 - x1这一级的意义是完成最小粒度的 DFT 合并。虽然看起来只是加减但它为下一级正确分组打好了基础。2.3 第 2 级蝶形跨度 2开始出现真正的复数乘法第二级蝶形处理块长为 4共有 2 个蝶形组。每个组内有 2 个蝶形第一组(0,2), (1,3) 第二组(4,6), (5,7)这一级的旋转因子是 W_4^0 和 W_4^1W_4^0 1 W_4^1 e^{-j2π/4} -j从这一级开始蝶形运算中会出现真正的复数乘法。每组的旋转因子数量等于当前级块长的一半也就是当前级的“蝶形跨度 m”。2.4 第 3 级蝶形跨度 4旋转因子为 W_8^0 到 W_8^3第三级蝶形处理块长为 8此时只剩下 1 个蝶形组组内有 4 个蝶形(0,4), (1,5), (2,6), (3,7)旋转因子是W_8^0, W_8^1, W_8^2, W_8^3其中W_8^1 e^{-j2π/8} cos(π/4) - j·sin(π/4) ≈ 0.7071 - 0.7071j W_8^2 -j W_8^3 -0.7071 - 0.7071j完成第三级后输出已经是按自然频率顺序排列的结果无需再重排。这一点很多初学者容易混淆输入需要位倒序输出却不需要。把三级蝶形的基本参数汇总级数块长蝶形组数每组蝶形数跨度 m旋转因子第 1 级2411W_2^0第 2 级4222W_4^0, W_4^1第 3 级8144W_8^0, W_8^1, W_8^2, W_8^32.5 为什么蝶形级数一定是 log₂N观察 8 点例子第 1 级跨度 1第 2 级跨度 2第 3 级跨度 4。每一级跨度翻倍块长也翻倍。8 点需要 3 级16 点需要 4 级N 点需要的级数是级数 log₂N每一级包含 N/2 个蝶形运算而每个蝶形运算只有 1 次复数乘法因此总计算量约为(N/2) · log₂N 次复数乘法这就是 FFT 计算复杂度 O(N log₂N) 的来源。3. 用 C 语言实现最小可运行的基2 FFT3.1 验证环境不需要开发板也能跑通很多读者误以为学 FFT 必须要有 STM32 开发板或者 DSP 开发环境。实际上理解蝶形结构只需要一台普通电脑和 C 编译器。本文用一个独立的 C 文件实现完整的基 2 时间抽取 FFT函数内部直接体现位倒序和三级蝶形循环。建议环境项目说明编译器GCC、Clang 或 MSVC 均可数学库使用 cosf/sinf链接时加 -lm运行平台Linux、macOS、Windows 都行输入数据程序内构造测试序列不依赖外部文件这个实现使用 float 类型便于后续移植到带 FPU 的 MCU如果学习环境是 PC也可以把 float 改成 double 提高精度。3.2 完整 C 代码位倒序加三级蝶形循环下面这个fft_dit函数实现了任意 2 的幂点数的基 2 时间抽取 FFT。代码故意写成教学风格避免使用复杂宏和指针技巧方便对照蝶形图阅读。#include stdio.h #include stdint.h #include math.h #ifndef M_PI #define M_PI 3.14159265358979323846 #endif typedef struct { float re; float im; } complex_t; // 将 x 的二进制位倒序后返回bits 表示有效位宽 static uint32_t bit_reverse(uint32_t x, uint32_t bits) { uint32_t y 0; for (uint32_t i 0; i bits; i) { y (y 1) | (x 1); x 1; } return y; } // 基2 时间抽取 FFT原位运算n 必须是 2 的幂 void fft_dit(complex_t *x, uint32_t n) { if (n 2 || (n (n - 1)) ! 0) { printf(error: n must be a power of 2\n); return; } uint32_t bits 0; while ((1u bits) n) { bits; } // 1. 输入按二进制位倒序重排 for (uint32_t i 0; i n; i) { uint32_t j bit_reverse(i, bits); if (j i) { complex_t tmp x[i]; x[i] x[j]; x[j] tmp; } } // 2. 蝶形运算size 是当前级的跨度 m也是每组蝶形个数 for (uint32_t size 1; size n; size 1) { uint32_t step size 1; // 当前级的块长 for (uint32_t start 0; start n; start step) { for (uint32_t k 0; k size; k) { // 旋转因子 W e^{-j2π·k / step} float angle -2.0f * (float)M_PI * (float)k / (float)step; float wr cosf(angle); float wi sinf(angle); uint32_t p start k; uint32_t q p size; float tr wr * x[q].re - wi * x[q].im; float ti wr * x[q].im wi * x[q].re; x[q].re x[p].re - tr; x[q].im x[p].im - ti; x[p].re x[p].re tr; x[p].im x[p].im ti; } } } } int main(void) { const uint32_t N 8; complex_t x[N]; // 输入为单位脉冲x(0)1其余为0 for (uint32_t i 0; i N; i) { x[i].re 0.0f; x[i].im 0.0f; } x[0].re 1.0f; fft_dit(x, N); for (uint32_t i 0; i N; i) { printf(X[%u] %8.4f %8.4fj\n, i, x[i].re, x[i].im); } return 0; }这段代码有三个关键点。第一位倒序使用bit_reverse(i, bits)计算目标位置只在j i时交换避免重复交换。第二蝶形循环是三层结构外层层级、中层块起始位置、内层蝶形序号。第三旋转因子直接在循环内用cosf/sinf计算方便理解生产环境会预计算成查找表。3.3 编译、运行与预期结果在 Linux 或 macOS 终端执行gcc -stdc11 -O2 fft_dit.c -lm -o fft_dit ./fft_dit单位脉冲的 DFT 理论上在所有频点都等于 1。运行结果应该接近X[0] 1.0000 0.0000j X[1] 1.0000 0.0000j X[2] 1.0000 0.0000j X[3] 1.0000 0.0000j X[4] 1.0000 0.0000j X[5] 1.0000 0.0000j X[6] 1.0000 0.0000j X[7] 1.0000 0.0000j由于浮点计算精度输出中可能出现极小虚部例如-0.0000j这是正常现象。单位脉冲测试能够快速发现位倒序和蝶形配对是否错误是一个非常有用的最小验证用例。注意这段 C 代码的目标是讲清楚蝶形运算原理并不是生产级 FFT 库。生产环境还需要考虑归一化方式、溢出处理、查找表优化、实数 FFT 封装等细节。4. 关键代码与参数详解4.1 位倒序函数为什么循环里要同时右移和左移bit_reverse函数是 FFT 实现中第一个容易写错的地方。static uint32_t bit_reverse(uint32_t x, uint32_t bits) { uint32_t y 0; for (uint32_t i 0; i bits; i) { y (y 1) | (x 1); x 1; } return y; }每一轮循环中x 1取出待反转数的最低有效位y 1把已经积累的位向左移动为新位腾出空间x 1把原数右移一位准备取下一位。以 3 位宽、输入 x1 为例二进制是001。三轮循环后第 1 轮取出最低位 1y 1x 右移变成 0。第 2 轮取出最低位 0y (1 1) | 0 2x 仍为 0。第 3 轮取出最低位 0y (2 1) | 0 4。结果 4正好是100也就是001的位倒序。常见错误是把右移和左移的时机写错或者忘记在交换前判断j i导致同一个元素被交换两次结果恢复原状。4.2 蝶形循环的三层结构级、块、蝶三层循环的参数直接影响结果for (uint32_t size 1; size n; size 1) { uint32_t step size 1; for (uint32_t start 0; start n; start step) { for (uint32_t k 0; k size; k) {size当前级的跨度同时也是每组内的蝶形个数。从 1 开始每级翻倍直到等于 n/2。step当前级的块长等于size 1。块与块之间不共享旋转因子。start每个块的起始位置从 0 开始以块长递增。k块内的第几个蝶形决定使用哪个旋转因子。每次蝶形操作的两个数据位置是p start k和q p size。这个配对关系对应蝶形图中的跨线连接方式。4.3 旋转因子的计算与预计算优化旋转因子的公式是W_当前级 e^{-j2π·k / step}对应代码float angle -2.0f * (float)M_PI * (float)k / (float)step; float wr cosf(angle); float wi sinf(angle);注意角度是负数对应 DFT 定义中的e^{-j2πkn/N}。如果使用正号频谱会变成镜像结果看起来类似“频率左右颠倒”这个问题在排查节会详细说明。在嵌入式设备上直接在循环内调用cosf/sinf计算耗时明显。生产项目通常把旋转因子预计算到数组float wr[N/2]; float wi[N/2]; for (uint32_t i 0; i N/2; i) { float angle -2.0f * (float)M_PI * (float)i / (float)N; wr[i] cosf(angle); wi[i] sinf(angle); }然后在蝶形内部按索引取用避免重复计算三角函数。4.4 采样率、点数和频率分辨率的换算工程中使用 FFT 时经常要把谱线序号换算成实际频率Δf Fs / N f_k k · Δf k · Fs / N其中 Fs 是采样率N 是 FFT 点数Δf 是频率分辨率k 是谱线索引。以 Fs8000Hz、N1024 为例FFT 点数采样率 Fs频率分辨率 Δf第 128 条谱线对应频率10248000 Hz7.8125 Hz1000 Hz102444100 Hz约 43.07 Hz约 5512 Hz20488000 Hz约 3.906 Hz500 Hz409644100 Hz约 10.77 Hz约 1378 Hz这里最常见的坑是直接把谱线索引当频率使用。比如在 44100Hz 采样率下索引 128 并不是 128Hz而是约 5512Hz。写代码时如果忘记除以 N 再乘以 Fs频谱图的横坐标就会完全错误。5. 验证 FFT 实现正确性的三种方法5.1 用单位脉冲验证直流分量单位脉冲 x(0)1、其余为 0 的信号其频谱在所有频点都是 1。这个特性让它可以同时验证位倒序是否正确。蝶形加法和减法是否正确。旋转因子索引是否正确。输出是否按自然频率顺序排列。运行第 3 节的 C 程序如果输出全为 1说明基本流程正确。如果输出出现乱序或不是全 1优先检查位倒序函数和蝶形配对关系。5.2 用正弦波验证峰值频率单位脉冲能发现错误但不足以证明旋转因子方向和频谱语义完全正确。更可靠的验证是输入单频正弦波检查频谱峰值是否落在预期索引上。以 Fs8000Hz、N1024、信号频率 f01000Hz 为例频率分辨率 Δf7.8125Hz峰值应落在索引 128。用 Python 可以快速验证import numpy as np fs 8000 N 1024 t np.arange(N) / fs f0 1000 x np.sin(2 * np.pi * f0 * t) X np.fft.fft(x) mag np.abs(X[:N//2]) idx np.argmax(mag) freq idx * fs / N print(fpeak index: {idx}) print(festimated frequency: {freq} Hz)输出应为peak index: 128 estimated frequency: 1000.0 Hz如果把旋转因子符号写反峰值会出现在镜像位置 N-idx 处或者虚部符号异常。这个测试比单位脉冲更严格。5.3 用 Parseval 定理验证能量守恒FFT 本质上是正交变换功率或者能量在变换前后应当守恒。Parseval 定理的离散形式是Σ|x(n)|² (1/N) · Σ|X(k)|²也就是说频域能量需要除以 N 之后才等于时域能量。用 Python 验证x np.random.randn(N) 1j * np.random.randn(N) X np.fft.fft(x) energy_time np.sum(np.abs(x)**2) energy_freq np.sum(np.abs(X)**2) / N print(energy_time, energy_freq)时域能量和频域能量应当非常接近。这个检查对自定义 C 实现的 FFT 特别有用因为它能发现蝶形运算中是否存在整体性的增益或能量泄漏问题。5.4 工程中的验证思路DSP库FFT和相位测量在 STM32F407 这类带 FPU 和 DSP 指令的 MCU 上通常会直接使用 CMSIS-DSP 库的arm_cfft_f32函数。库函数内部同样基于蝶形运算只是做了流水线优化和查找表预计算。理解蝶形结构后设置库函数参数应该关注三件事fftLen必须是 2 的幂例如 256、512、1024、2048。ifftFlag为 0 表示正变换为 1 表示反变换。bitReverseFlag为 1 表示库负责位倒序为 0 表示输入已经是倒序排列。使用 DSP 库测量相位时常见做法是计算atan2(imag, real)得到每个频点的相位。但要注意如果输入信号频率不是 Δf 的整数倍频谱泄漏会直接污染相位读数。生产环境中通常需要加窗函数例如汉宁窗再在对应主瓣附近进行插值才能得到稳定可信的相位值。6. 蝶形运算的常见错误与排查路径6.1 位倒序错误输出顺序乱单位脉冲不通过现象输入单位脉冲后输出不是全 1或者在某个位置出现异常大值。可能原因bit_reverse的位宽计算错误比如 N8 时算成 2 位。交换时没有判断j i导致交换两次顺序还原。在蝶形循环之前没有执行位倒序。检查方式打印位倒序前后的索引对应关系和 N8 的倒序表 0,4,2,6,1,5,3,7 对比。解决方式用固定 N8 单独调试倒序表确认无误后再接入蝶形循环。6.2 旋转因子方向错误频谱镜像或虚部符号异常现象输入是 1000Hz 正弦波峰值索引出现在 N/2 之后的镜像位置或者相位计算结果和预期相反。可能原因角度符号写反。FFT 使用e^{-j2πk/N}如果写成e^{j2πk/N}频谱会做共轭镜像。检查方式用单频正弦波测试打印峰值索引和atan2(imag, real)。如果峰值位置是N - 128而不是128说明符号有问题。解决方式将角度计算统一改为float angle -2.0f * (float)M_PI * (float)k / (float)step;6.3 归一化缺失幅度差 N 倍现象输入正弦波幅值为 1FFT 后在峰值频点的幅度约为 N/2而不是 1。原因这不是程序错误而是 DFT/FFT 本身没有做归一化。正向 FFT 的幅度会放大 N 倍单频正弦波的能量会分到正负频率两侧所以峰值幅度约为 N/2。处理方式如果后续只做相对比较不关心绝对幅度可以不归一化。如果要求绝对幅度单边谱需要把峰值除以 N考虑负频率对称时通常乘以 2/N。反变换时再执行 1/N 归一化。检查方式把 FFT 输出的峰值幅度除以 N看是否接近原始信号幅值。6.4 输入长度不是 2 的幂现象程序崩溃、越界、输出长度不正确。原因基 2 FFT 要求 N 必须是 2 的幂。如果直接输入 1000 点循环size n和位倒序位宽都会出错。处理方式数据补零到不小于原始长度的下一个 2 的幂。补零会改变频率分辨率但不会改变原始频率成分的位置。改用混合基 FFT支持 12、15、60 等非 2 的幂点数。直接在数据采集阶段设置 FFT 点数为 1024、2048 等 2 的幂。6.5 排错总表问题现象常见原因检查方式处理建议单位脉冲输出不是全 1位倒序错误或交换重复打印倒序索引表用 N8 固定用例单独验证峰值在镜像位置旋转因子符号写反检查角度是否带负号统一为 -2πk/step幅度比预期大 N 倍正向 FFT 未归一化除以 N 对比幅度根据单边/双边谱决定归一化程序越界或死循环输入点数不是 2 的幂检查 N 的二进制表示补零到 2 的幂或换混合基频率横坐标不对直接用索引当频率计算 idx*Fs/N换算频率再绘图7. 从“会用”到“会选”FFT变体、工程实践与扩展7.1 基2、基4、混合基与分裂基的取舍蝶形运算并不只有基 2 一种结构。基 4 FFT 每级处理 4 个输入一次蝶形包含更多复用乘法和加法次数更少混合基 FFT 支持更多非 2 的幂点数分裂基 FFT 结合基 2 和基 4是常见的实数序列高效实现方案。算法点数要求优势典型场景基 2 DIT2 的幂结构简单易理解和实现学习、教学、通用库基 44 的幂乘法次数更少硬件资源敏感信号处理混合基可分解为多个因子支持非 2 的幂点数通用 FFT 库底层分裂基2 的幂实数序列效率高高性能 DSP 库在工程选型时先确认硬件平台支持的指令集和库函数不要为了炫技手写复杂变体。多数 MCU 场景使用厂商提供的 FFT 库已经足够。7.2 复数FFT与实数FFTFFT 的数学定义基于复数输入输出但实际采集信号通常是实数序列。直接用复数 FFT 计算实数序列会浪费一半内存和一半计算量。常见优化思路有两种。第一种把一个 N 点实数序列拆成 N/2 点复数序列做一次 N/2 点 FFT再利用共轭对称性恢复结果。第二种用 N 点复数 FFT 同时处理两个 N 点实数序列一个放在实部一个放在虚部后处理拆分。CMSIS-DSP 提供的arm_rfft_fast_f32就是针对实数序列的封装底层仍然用复数蝶形但对输入输出做了重排。理解蝶形结构有助于理解这类封装函数的N/2点处理和输出排列方式。7.3 学习环境与生产环境的差异学习阶段用 PC 的 C 或 Python 跑通逻辑和生产环境之间的差异不能忽略。维度学习环境生产环境数据类型float 或 double固定点 Q15/Q31 或 float旋转因子每次现算 cos/sin预计算查找表输入数据手工构造序列ADC 采集、DMA 搬运异常处理只关注算法正确性越界、溢出、中断竞争验证手段打印结果对比频谱图、协议分析、监控告警性能要求可接受慢速实时性、流水线、缓存友好在 MCU 上做实时频谱显示时ADC 采集、DMA 传输和 FFT 计算往往需要配合双缓冲机制。一块缓冲用于采集另一块用于 FFT 计算交替使用才能避免数据被覆盖。7.4 可复用清单从理论理解到工程落地针对 FFT 蝶形运算的学习和工程落地整理一份检查清单[ ] 确认输入长度 N 是 2 的幂。[ ] 确认采样率 Fs 和频率分辨率 Δf 的关系。[ ] 检查位倒序索引表是否和理论一致。[ ] 用单位脉冲测试确认输出全 1。[ ] 用单频正弦波测试确认峰值索引和理论值一致。[ ] 用 Parseval 定理验证能量守恒。[ ] 确认归一化方式是单边谱还是双边谱。[ ] 确认旋转因子符号与 DFT 定义一致。[ ] 生产环境预计算旋转因子查找表。[ ] 生产环境将 FFT 结果用 DMA 或双缓冲保护。[ ] 加窗函数后再做幅度或相位测量。[ ] 输出频率需要重新换算不能直接使用索引。7.5 扩展方向加窗、实时频谱仪与包络谱分析理解蝶形运算后下一步通常走向工程应用。实时频谱分析需要关注窗函数对频谱泄漏的抑制常用的汉宁窗、汉明窗、布莱克曼窗各有旁瓣衰减特性。故障诊断领域常把 FFT 和希尔伯特变换结合先求包络信号再对包络做 FFT得到包络谱用于滚动轴承等旋转机械的早期故障识别。在电力系统直流采样、谐波分析、音频可视化、振动状态监测等场景中FFT 都是基础计算模块。只要把本次讲清的位倒序、旋转因子、蝶形级联和归一化四个关键点吃透后续接触任何 FFT 库或者专用 DSP 芯片都只是在此基础上做性能优化和接口封装。日常练习建议是自己写一遍基 2 FFT用单位脉冲和正弦波两种输入验证再尝试把旋转因子改成查找表对比性能变化。这个过程比单纯调用库函数更能建立对数字信号处理的直觉。