数字信号处理DSP领域FFT快速傅里叶变换几乎是每个人都会遇到的核心内容。很多同学在学习时都有类似的体验DFT 的公式能看懂频谱图也能画出来但一看到 8 点 FFT 的蝶形运算流图满屏交叉的连线和旋转因子符号就让人发懵更别提自己动手实现一版可运行的 FFT 程序。本文不打算只讲“怎么调用库函数”而是把 FFT 内部最关键的蝶形运算结构拆开讲透。我会从 DFT 的暴力计算出发一步步推导出蝶形流图再手写一个可以直接运行的基 2-FFT 程序最后延伸到包络谱分析、相位测量和嵌入式实时频谱这类实际工程场景。无论你是刚接触数字信号处理的学生还是在单片机、上位机上做频谱分析的开发者这篇文章都值得收藏备用。1. FFT 与蝶形运算为什么这一结构是核心1.1 FFT 解决了什么问题FFT 的全称是 Fast Fourier Transform即快速傅里叶变换。它并不是一种新的变换而是离散傅里叶变换DFTDiscrete Fourier Transform的快速计算算法。如果不做任何优化直接按照 DFT 公式计算一个 N 点的频谱需要 N 次复数乘法和 N(N-1) 次复数加法总计算量大约是 N² 的复数运算。当 N1024 时直接计算需要约 100 万次复数乘法。这个量级在普通 PC 上尚可接受但在单片机、实时信号处理系统里就很容易导致系统卡顿或无法满足实时性要求。FFT 利用旋转因子的周期性和对称性把 DFT 的计算复杂度从 O(N²) 降到了 O(N log₂N)。同样是 1024 点基 2-FFT 大约只需要 (N/2)·log₂N ≈ 5120 次复数乘法比直接计算少了大约 200 倍差距非常明显。正是这个性能提升让频谱分析在实时系统、嵌入式设备和通信系统中变得可行。1.2 蝶形运算在 FFT 中的地位FFT 之所以快核心在于它把一个大尺寸的 DFT 不断拆分成小尺寸 DFT再用一种固定形式的运算单元把结果合并起来。这个固定形式的运算单元就是“蝶形运算”。一个蝶形运算单元包含两个输入、一个旋转因子和两个输出形状像一只蝴蝶因此得名。整个 FFT 算法无论采用按时间抽取DIT还是按频率抽取DIF最终都会落到一组蝶形运算的重复执行上。理解蝶形运算结构不只是为了看懂教科书上的流图。实际工作中你会发现FFT 的优化、定点化、硬件加速器的设计、以及各种 DSP 库的源码几乎都围绕蝶形运算展开。掌握了蝶形结构你才能真正理解 FFT 的行为遇到频谱结果不对时也知道从哪里排查。2. DFT 到 FFT从暴力计算到分治思想2.1 先回顾 DFT 公式DFT 的定义式如下X[k] Σ_{n0}^{N-1} x[n] · W_N^{nk}, k 0, 1, ..., N-1其中x[n] 是时域第 n 个采样点X[k] 是频域第 k 个频率分量W_N 是旋转因子W_N e^(-j·2π/N)。如果直接按这个公式写双重循环就是最原始的 DFT 实现。时间复杂度 O(N²) 的来源就在这里。从工程角度解释这个公式的含义每个频点 k 的计算本质上是在做“输入序列与一个复指数序列的逐点相乘再累加”。N 个频点每个都要遍历全部 N 个输入点计算量自然上去了。2.2 旋转因子的两个关键性质FFT 能大幅降低计算量靠的是旋转因子 W_N 的两个数学性质对称性W_N^(k N/2) -W_N^k周期性W_N^(k N) W_N^k这两个性质意味着在计算过程中很多旋转因子是重复的或者可以取负号直接复用。蝶形运算正是把这两个性质用到了极致。2.3 分治思路奇偶分组基 2-FFT 的基本思想是分治。以按时间抽取DIT-FFT为例把 N 点的序列 x[n] 按序号 n 的奇偶分成两组偶数项x[0], x[2], x[4], ...奇数项x[1], x[3], x[5], ...分别计算这两组的 N/2 点 DFT再利用旋转因子的性质合并成 N 点 DFT。这还没完。分出来的 N/2 点 DFT可以继续按奇偶分组再拆成 N/4 点。当 N 是 2 的整数次幂时可以一直拆到 2 点 DFT 为止。2 点 DFT 是最小单元只需要一次加法、一次减法。从整体看整个过程就像一棵递归树每一层都有 N/2 个蝶形运算一共有 log₂N 层。总蝶形数就是 (N/2)·log₂N对应前面提到的计算量。3. 蝶形运算结构详解3.1 单个蝶形单元长什么样一个最基本的 DIT 蝶形单元输入是两个复数 a 和 b输出是两个复数 A 和 B中间乘一个旋转因子 W_N^r。计算公式如下A a W_N^r · b B a - W_N^r · b看起来非常简单但它揭示了一个重要的工程细节计算时必须先把 b 乘以旋转因子得到中间结果再和 a 做加减。如果直接修改 a 和 b 的存储位置后面的计算就会出错。所以蝶形运算通常需要临时变量保存中间结果。在矩阵运算或神经网络里这类操作叫“in-place 更新”FFT 同样支持原地计算只需要注意读写顺序即可。3.2 8 点 DIT-FFT 流图以 8 点 FFT 为例整个流图分三层每层 4 个蝶形运算。我用文字把结构画出来x 表示交叉连接最终数据从左往右流动输入(倒位序) 第1层 第2层 第3层 x[0] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[0] x[4] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[1] x[2] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[2] x[6] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[3] x[1] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[4] x[5] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[5] x[3] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[6] x[7] ------ 蝶形 ---- 蝶形 ---- 蝶形 ---- X[7]每一层的蝶形运算单元两个输入之间的间隔是不同的第 1 层间隔为 1 个点即相邻两两一组第 2 层间隔为 2 个点第 3 层间隔为 4 个点。如果从流图内部观察可以发现每一层的蝶形运算都遵循固定的模式。这也是写代码时最重要的依据。3.3 倒位序与旋转因子的摆放规律细心的读者会发现8 点 DIT-FFT 的输入不是 x[0], x[1], x[2], ... 的自然顺序而是 x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]这个顺序叫“倒位序”。所谓倒位序是把序号写成二进制再把二进制位反过来。例如0 000 → 000 0 1 001 → 100 4 2 010 → 010 2 3 011 → 110 6所以 DIT-FFT 的流程是先对输入做倒位序重排再进行蝶形运算输出就是自然顺序的频谱 X[k]。旋转因子的规律同样非常整齐。设当前层蝶形间隔为 m则该层的旋转因子为W_N^r, r 0, 1, ..., N/(2m) - 1以 8 点 FFT 为例第 1 层 m1旋转因子为 W_2^0第 2 层 m2旋转因子为 W_4^0, W_4^1第 3 层 m4旋转因子为 W_8^0, W_8^1, W_8^2, W_8^3。每一层的旋转因子个数正好是本层蝶形组数的一半而且每个组内共用同一个旋转因子。提前把这些旋转因子计算好并缓存成表是 FFT 性能优化的重要一步。3.4 按频率抽取DIF-FFT与 DIT 的对比除了按时间抽取还有一种常见的实现叫按频率抽取DIF-FFT。两者的关系可以这样理解DIT-FFT输入倒位序输出自然序先做倒位序重排再做蝶形运算DIF-FFT输入自然序输出倒位序先做蝶形运算再做倒位序重排。DIF 蝶形单元公式如下A a b B (a - b) · W_N^r可以看到DIF 是先加减后乘旋转因子而 DIT 是先乘后加减计算顺序正好相反。这两种结构在数学上是等价的工程上选择哪一种主要看你对输入输出排序的要求以及硬件流水线设计的方便程度。4. 手写基 2-FFT 完整实战Python4.1 环境准备本实战使用 Python 3配合 numpy 做验证但核心 FFT 代码只用标准库就能完成。建议在本地安装 numpy方便后续和 numpy.fft.fft 的结果做对比验证。文件结构很简单只需要一个 Python 脚本文件例如my_fft.py。4.2 第一步实现倒位序函数倒位序是 DIT-FFT 的入口先实现它。def bit_reverse(n, bit_width): 对一个整数 n 做 bit_width 位的二进制倒序。 例如 bit_width3 时n1(001) 返回 4(100)。 result 0 for i in range(bit_width): if (n i) 1: result | 1 (bit_width - 1 - i) return result # 验证 for i in range(8): print(i, -, bit_reverse(i, 3))运行结果0 - 0 1 - 4 2 - 2 3 - 6 4 - 1 5 - 5 6 - 3 7 - 7这里的逻辑是从最低位开始扫描如果第 i 位是 1就把它放到结果的 (bit_width-1-i) 位。这样就把二进制位顺序完全反过来了。4.3 第二步实现蝶形运算核心代码有了倒位序和蝶形单元公式完整 FFT 就很好写了import cmath def fft_dit(x): 基 2 按时间抽取 FFT。 输入 x 的长度必须是 2 的整数次幂。 N len(x) bit_width N.bit_length() - 1 # 第一步倒位序重排统一转为复数 a [complex(x[bit_reverse(i, bit_width)]) for i in range(N)] # 第二步逐层蝶形运算 length 2 # 当前蝶形运算的间隔从 2 开始每层翻倍 while length N: # 当前层使用的基础旋转因子 W_length angle -2 * cmath.pi / length w_len complex(cmath.cos(angle), cmath.sin(angle)) # 按组遍历 for i in range(0, N, length): w complex(1, 0) # 每组从 W^0 开始 half length // 2 for j in range(i, i half): u a[j] v a[j half] * w a[j] u v a[j half] u - v w * w_len # 更新旋转因子 length * 2 return a这段代码的核心结构是三层循环外层循环while length N控制层数中间循环for i in range(0, N, length)按组遍历每组长度为 length内层循环for j in range(i, i half)完成组内的 half 个蝶形运算。每个蝶形运算都遵循同一模式先保存 u再计算 v最后原地更新两个输出。4.4 第三步与 numpy 结果对比下面用一组测试信号验证计算结果import numpy as np if __name__ __main__: # 构造 8 点测试信号随便取一组数值即可 x [1.0, 2.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0] y fft_dit(x) y_np np.fft.fft(x) print(k 手写FFT numpy.fft 误差) for i in range(len(y)): err abs(y[i] - y_np[i]) print(f{i} {y[i]:.4f} {y_np[i]:.4f} {err:.2e})预期输出中手写 FFT 和 numpy 的结果应该完全一致误差在 1e-14 量级。所谓“量级”指的是浮点数运算本身带来的微小舍入误差可以忽略不计。这个对比非常重要它证明了我们手写的蝶形运算结构本身是正确的。后续如果在嵌入式上移植、优化也可以用同样的方式与标准库做交叉验证。4.5 加一行打印直观观察蝶形过程为了帮助理解蝶形运算的过程我在代码里加一个带打印的版本每一层输出当前数组的状态def fft_dit_debug(x): N len(x) bit_width N.bit_length() - 1 a [complex(x[bit_reverse(i, bit_width)]) for i in range(N)] print(倒位序重排后:, [%.3f%.3fj % (v.real, v.imag) for v in a]) length 2 layer 1 while length N: angle -2 * cmath.pi / length w_len complex(cmath.cos(angle), cmath.sin(angle)) for i in range(0, N, length): w complex(1, 0) half length // 2 for j in range(i, i half): u a[j] v a[j half] * w a[j] u v a[j half] u - v w * w_len print(f第{layer}层蝶形结果:, [%.3f%.3fj % (v.real, v.imag) for v in a]) length * 2 layer 1 return a fft_dit_debug([1.0, 2.0, 1.0, 0.5, 0.0, 0.0, 0.0, 0.0])运行后可以看到倒位序重排后: [1.0000.000j, 0.0000.000j, 2.0000.000j, 0.0000.000j, 1.0000.000j, 0.0000.000j, 0.5000.000j, 0.0000.000j] 第1层蝶形结果: [1.0000.000j, 1.0000.000j, 2.0000.000j, 2.0000.000j, 1.0000.000j, 1.0000.000j, 0.5000.000j, 0.5000.000j] ...从这些中间结果可以看出每一层蝶形运算都在逐步把时域信息“混叠”成频域信息。理解这一点后再回头看流图就不会觉得抽象了。5. FFT 的工程应用包络谱、相位测量与实时频谱5.1 包络谱分析原理与流程FFT 最常见的应用是频谱分析但在设备故障诊断领域直接对原始振动信号做频谱分析往往不够直观。比如滚动轴承出现故障时高频冲击成分会周期性出现这种周期性正是故障特征。直接看原始信号频谱时高频载波和低频调制会混在一起很难提取故障特征频率。这时候可以用包络谱分析。核心思路是先提取信号的包络波形再对包络做 FFT得到包络谱。包络提取通常借助希尔伯特变换Hilbert Transform。对原始信号 x(t) 做希尔伯特变换得到信号 ĥ(t)构造解析信号z(t) x(t) j·ĥ(t)解析信号的幅值 |z(t)| 就是原始信号的包络。把这个包络信号再做一次 FFT得到的就是包络谱。这个流程在滚动轴承故障诊断、齿轮箱监测、预测性维护中非常常见。这也是为什么在很多 DSP 资料里FFT 和希尔伯特变换总是成对出现。5.2 Python 实现包络谱用 scipy 可以非常简洁地实现包络谱分析from scipy.signal import hilbert import numpy as np import matplotlib.pyplot as plt # 构造一个调幅信号1kHz 载波100Hz 调制 fs 20000 t np.arange(0, 1, 1 / fs) carrier np.sin(2 * np.pi * 1000 * t) modulation 1 0.8 * np.sin(2 * np.pi * 100 * t) signal carrier * modulation # 希尔伯特变换提取包络 analytic_signal hilbert(signal) envelope np.abs(analytic_signal) # 对包络做 FFT env_fft np.fft.fft(envelope) freqs np.fft.fftfreq(len(t), 1 / fs) # 画图 plt.figure(figsize(10, 5)) plt.subplot(2, 1, 1) plt.plot(t[:500], signal[:500], label原始调幅信号) plt.plot(t[:500], envelope[:500], label包络波形, linewidth2) plt.legend() plt.subplot(2, 1, 2) plt.plot(freqs[:len(freqs) // 2], np.abs(env_fft[:len(freqs) // 2])) plt.xlim(0, 500) plt.xlabel(频率 (Hz)) plt.ylabel(幅值) plt.title(包络谱) plt.tight_layout() plt.show()原始信号频谱里能量集中在 1000Hz 附近而从包络谱中可以清楚看到 100Hz 的调制频率成分。这正是包络谱的价值它把高频载波“剥离”让低频调制信息暴露出来。5.3 用 FFT 测量相位FFT 不仅能看幅值还能测相位。频谱第 k 个频点的相位可以直接通过复数结果的辐角得到phase atan2(imag(X[k]), real(X[k]))在 Python 里对应np.angle(X[k])。这里有一个非常容易踩的坑对正弦信号 sin(2πft φ) 做 FFT测出的峰值 bin 相位并不是 φ而是 φ - π/2。原因是 FFT 的基函数是复指数 e^(-j2πk n/N)而正弦信号可以写成两个复指数的叠加其中正频率部分会引入 -π/2 的固定相移。示例import numpy as np fs 1024 N 1024 f0 50 phi_true 0.7 t np.arange(N) / fs x np.sin(2 * np.pi * f0 * t phi_true) X np.fft.fft(x) k int(f0 * N / fs) phi_measured np.angle(X[k]) print(真实相位:, phi_true) print(FFT测量相位(未补偿):, phi_measured) print(补偿后:, phi_measured np.pi / 2)如果 f0 正好是频率分辨率的整数倍fs/N 的整数倍补偿后可以测得很准。如果 f0 不在整数倍 bin 上频谱泄漏会引入额外的相位误差需要加窗并做相位修正或者改用 Goertzel 算法对单频点做精确测量。6. 嵌入式平台的 FFT 实现要点6.1 为什么嵌入式要用 DSP 库在 STM32 这类 MCU 上做 FFT有两条路线自己写纯 C 的蝶形运算或者直接使用官方 DSP 库。自己写的好处是理解深入、可定制性强缺点是调试周期长而且很难在性能上超越芯片厂商专门优化的库。以 ARM Cortex-M 系列为例CMSIS-DSP 库提供了现成的 FFT 函数包括复数 FFTarm_cfft_f32和实数 FFTarm_rfft_fast_f32。这些函数针对 Cortex-M 内核的浮点单元和指令流水线做了专门优化调用起来非常方便。对于 STM32F407 这类带 FPU 的 M4 内核芯片浮点 FFT 可以直接跑实时性能在多数场景下都能满足要求。如果使用不带 FPU 的 M0/M3则建议使用定点 Q15 或 Q31 版本的库函数避免浮点模拟带来的性能损失。6.2 CMSIS-DSP 库的使用思路我不打算贴一长串与特定版本绑定的代码因为 CMSIS-DSP 在不同版本之间存在差异。下面给出一个清晰的使用思路具体 API 以你所用 SDK 的官方头文件为准。以 CMSIS-DSP 的实数 FFT 为例典型流程如下// 伪代码核心流程示意具体 API 以你使用的 CMSIS-DSP 版本为准 #include arm_math.h #define FFT_SIZE 1024 float32_t input[FFT_SIZE]; // 时域输入ADC 采样结果 float32_t fft_output[FFT_SIZE]; // 频域输出 arm_rfft_fast_instance_f32 fft_inst; void fft_init(void) { arm_rfft_fast_init_f32(fft_inst, FFT_SIZE); } void fft_run(void) { // 执行实数 FFT第三个参数 0 表示正变换 arm_rfft_fast_f32(fft_inst, input, fft_output, 0); // fft_output 的前 FFT_SIZE/2 个点对应有效频谱 }要点是初始化结构体要在 FFT 之前完成输入数组长度和初始化时指定的点数一致实数 FFT 的输出布局和复数 FFT 不同取幅值前先确认存储格式如果做的是滑窗分析每次填充新数据时要注意缓存区边界的处理。6.3 实时频谱显示的工程结构在 F407 这类 MCU 上做实时频谱显示一般可以拆成这么几个模块ADC 采样 → 数据缓存 → FFT 运算 → 幅值计算 → 显示刷新实际工程里最需要注意的是“不丢数据”。ADC 采样的速率和 FFT 计算耗时不匹配时通常会采用双缓冲机制DMA 把 ADC 数据写入缓冲区 A 时CPU 对缓冲区 B 的数据做 FFT一轮完成后交换缓冲区。这样能保证采样连续不中断FFT 的计算时间只要小于一帧采样时间系统就能稳定工作。幅值计算可以进一步优化。对实数 FFT 的结果幅值通常用sqrt(re² im²)如果只是显示可以使用查找表或 fast inverse sqrt 近似减少 CPU 开销。7. 常见问题与排查思路在实际开发中FFT 结果不对是高频问题。很多情况并不是 FFT 算法本身写错而是外围参数和数据处理细节出了问题。下面整理一张常见问题排查表问题现象常见原因解决思路频谱在 0Hz 处有巨大尖峰输入信号含直流偏置FFT 前先减去信号均值峰值幅值和理论值对不上未做幅值修正单边谱非直流分量乘 2/N频率分辨率不够采样点数少或采样率过高增大 N或降低采样率频谱泄漏严重非整周期截断加 Hanning、Hamming 等窗函数输出频谱和硬件实测不一致采样率设置错误或混叠确认采样率满足奈奎斯特定理嵌入式 FFT 卡顿浮点运算量过大或中断过密改定点库优化采样与计算调度再补充几个细节0Hz 尖峰问题。这是最容易出现的现象。很多传感器输出的信号都带有直流偏置直接做 FFTX[0] 会非常大导致其他频率分量在图上几乎看不见。解决办法很简单FFT 之前对整段数据减去均值也就是去直流。幅值修正。如果只是看频谱形状不修幅值没关系但要做定量分析必须注意。对实信号取 N 点 FFT单边频谱中非直流分量的真实幅值大约是2·|X[k]|/N直流分量是|X[0]|/N。频率分辨率。频率分辨率 Δf fs/N由采样率和 FFT 点数共同决定。分辨率不够时两个接近的频率峰会被混成一个。补零只能让频谱看起来更平滑并不能真正提高物理分辨率。频谱泄漏。截断非整周期信号会产生泄漏能量从真实频率“漏”到旁边。加窗可以缓解但会牺牲幅值精度工程中需要根据需求取舍。相位测量误差。前面提到过非整周期采样时直接看峰值 bin 的相位会不准。工程上常见做法是先用窗函数抑制泄漏再做相位修正或者用 Goertzel 算法对目标频率单独计算。8. 最佳实践与工程建议结合多年的数字信号处理开发经验我建议在工程中遵守以下几点采样参数先行。开始写代码之前先确认采样率、FFT 点数和目标频率范围。频率分辨率、频谱范围、计算耗时三个指标互相制约不能只盯着一个。点数优先选 2 的幂。基 2-FFT 在嵌入式库中支持最完善性能也最优。如果你的采样率导致 N 不满足 2 的幂可以补零到下一个 2 的幂但要清楚补零只改变频域插值密度不改变物理分辨率。建立输出校准流程。使用标准正弦信号作为测试输入验证 FFT 的幅值、频率、相位三项指标是否符合预期。把校准结果记录成文档后续更换硬件或采样率时能快速发现问题。旋转因子提前缓存。如果 FFT 在循环里反复执行每次重新计算三角函数会浪费大量 CPU 时间。工程上通常预计算一张旋转因子表运行时代码只做查表。注意数据缓存与边界。滑动窗分析中每次 FFT 输入和上一帧有重叠时要确保数据搬运正确。嵌入式开发中尤其要注意缓存区的越界写问题。优先使用厂商 DSP 库。在 MCU 上做 FFT除非有明确的定制需求否则优先使用 ST、ADI 等芯片厂商提供的 DSP 库。这些库经过充分验证性能和稳定性都比自己写的初版代码好。安全边界意识。在设备故障诊断等场景中FFT 结果是后续决策的依据。如果频谱数据被噪声污染或算法参数被误改可能导致误判。因此建议对关键参数做校验对异常数据做保护避免一条异常数据把整个缓冲区的结果带偏。9. 总结与学习路线这篇文章围绕 FFT 的蝶形运算结构讲清了三个层面的内容第一FFT 为什么快。核心是利用旋转因子的对称性和周期性把 DFT 的 O(N²) 计算降为 O(N log₂N)而蝶形运算就是把这种数学优化落到具体计算步骤的载体。第二蝶形运算的规律。按时间抽取的 DIT-FFT 需要先倒位序、再逐层蝶形每一层蝶形具有相同的间隔规则和旋转因子规律。手写实现时只要三个循环结构正确结果就和标准库一致。第三工程落地方法。从包络谱分析到嵌入式实时频谱FFT 的实战价值非常广泛。理解蝶形结构后使用 DSP 库、排查频谱异常、处理幅值和相位修正都会顺手很多。接下来的学习路线建议从三个方向展开先自己手写一个 8 点和 16 点的 FFT画出对应的蝶形流图把“代码”和“流图”在脑子里对应起来再研究一下定点 FFT 的实现思路理解浮点与定点数在蝶形运算中的差异最后去找一个真实项目比如滚动轴承振动信号分析或音频实时频谱显示把 FFT 放进完整的信号链路里去用。如果你最近正在做 FFT 相关开发遇到频谱结果不对的情况建议先按顺序查三件事采样率是否满足奈奎斯特定理FFT 点数是不是 2 的幂输入数据有没有做倒位序和去直流处理。这三件事查完绝大多数基础问题都能定位。希望这篇关于 FFT 蝶形运算结构的文章能帮你在数字信号处理的路上少走一些弯路。