FFT蝶形运算详解：从原理到C/Python实现与工程应用

FFT 是数字信号处理里出现频率最高的算法之一。无论是频谱分析、滤波器设计、调制解调还是设备故障诊断里的包络谱和边带分析底层几乎都落在快速傅里叶变换上。很多同学在课本上看到“蝶形运算”四个字会觉得抽象实际把它拆开看就是一组固定结构的复数加减运算再配合旋转因子做加权。理解了蝶形运算的结构和索引规律FFT 的实现和调试都会清晰很多。这篇文章不做空洞的概念铺陈直接从 DFT 的计算瓶颈讲起推导 FFT 的分解思路把蝶形运算的输入输出关系、旋转因子规律、码位倒序这些关键细节写清楚。然后给出基 2 时间抽取 FFT 的完整 C 语言和 Python 实现再补充工程上常见的 DSP 库应用、实时频谱分析、相位测量和包络谱分析思路。无论你是刚学数字信号处理的学生还是要在单片机、嵌入式平台或者 PC 上做频谱分析的工程师这篇文章都值得收藏备用。1. 核心能力速览能力项说明算法类型快速傅里叶变换FFT基 2 时间抽取为主核心结构蝶形运算由蝶形结、旋转因子、码位倒序组成计算量直接 DFT 为 O(N^2)FFT 为 O(N log2 N)输入要求序列长度通常为 2 的整数次幂如 256、512、1024、2048、4096适用信号一维实数或复数序列、时域波形、振动信号、音频信号输出内容频谱幅值、相位、功率谱、包络谱、频域滤波结果典型应用实时频谱显示、故障诊断、振动分析、电力谐波分析、通信信号处理工程落地可手写 C 实现也可调用 CMSIS-DSP、KissFFT、FFTW、SciPy 等库硬件门槛低单片机、ARM Cortex-M、PC CPU 均可运行资源占用与 N 和实现方式相关运行时内存约 O(N)堆叠递归版本需要额外栈空间对于绝大多数实时频谱分析场景先跑通一个 1024 点 FFT再根据采样率和帧长度调整参数基本就能覆盖常见的工程需求。蝶形运算理解之后不管换到哪个库原理都是相通的。2. FFT 是什么为什么需要蝶形运算离散傅里叶变换DFT的定义是X(k) Σ_{n0}^{N-1} x(n) · W_N^{kn}其中 W_N e^{-j2π/N}称为旋转因子。直接按这个公式计算 N 点 DFT每个 X(k) 都要进行 N 次复数乘法和 N-1 次复数加法全部 N 个频点累加起来复数乘法量约为 N^2。当 N1024 时N^2 约一百万次复数乘法当 N4096 时N^2 约一千六百万次。在嵌入式平台上这个开销很难做到高帧率实时处理。FFT 的核心思想不是“发明新的变换”而是利用旋转因子的三类性质把一个大 DFT 拆成多个小 DFT周期性W_N^{kN} W_N^k对称性W_N^{kN/2} -W_N^k可约性W_N^{2k} W_{N/2}^{k}蝶形运算就是这种拆解过程的最小计算单元。一次蝶形运算通常完成两个复数的线性组合组合系数是旋转因子形式固定为X1 X1 W · X2 X2 X1 - W · X2这样一次蝶形运算包含一次复数乘法、两次复数加法。整个 FFT 就是由很多层这样的蝶形结按特定规律连接起来的。理解蝶形运算重点就是理解三点每一层有几个蝶形、每个蝶形的序号间隔是多少、旋转因子取什么值。3. DFT 到 FFT 的分解推导以基 2 时间抽取DIT-FFT为例假设 N2^M把输入序列 x(n) 按 n 的奇偶分成两组x(2r) 和 x(2r1)其中 r 0, 1, ..., N/2 - 1于是 DFT 可以写成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}整理得X(k) G(k) W_N^k · H(k)其中G(k) Σ_{r0}^{N/2-1} x(2r) · W_{N/2}^{rk}H(k) Σ_{r0}^{N/2-1} x(2r1) · W_{N/2}^{rk}这里 G(k) 是偶数点序列的 N/2 点 DFTH(k) 是奇数点序列的 N/2 点 DFT。由于 G(k) 和 H(k) 都以 N/2 为周期而 W_N^{kN/2} -W_N^k所以当 k 从 0 取到 N/2-1 之后X(k) G(k) W_N^k · H(k)X(k N/2) G(k) - W_N^k · H(k)这正好是一次蝶形运算。也就是说一个 N 点 DFT 被拆成两个 N/2 点 DFT再用 N/2 个蝶形运算组合成 N 点结果。继续对 N/2 点 DFT 做同样的奇偶拆分直到拆成 2 点 DFT。2 点 DFT 的公式是X(0) x(0) x(1) X(1) x(0) - x(1)这也是最基础的蝶形结。整个分解过程一共需要 log2 N 级每一级有 N/2 次蝶形运算所以总复数乘法量约为 (N/2)·log2 N。这就是 FFT 比 DFT 快几个数量级的原因。以 1024 点为例直接 DFT 约 1048576 次复数乘法基 2 FFT 约 5120 次复数乘法相差约 200 倍。N 越大优势越明显。从频率抽取DIF也是一样的逻辑只是先把输出频点分成奇偶两组再逐级拆分蝶形结的结构和 DIT 相反。工程上 DIT 更常见本文后面都以 DIT 为主。4. 蝶形运算的三种基本结构4.1 基本蝶形结一个基本蝶形结有两个输入 A 和 B一个旋转因子 W两个输出A A W · B B A - W · B从复数运算角度如果 A、B、W 都是复数实际需要 4 次实数乘法和 6 次实数加减法。如果旋转因子的实部或虚部是 0 或 ±1还可以进一步简化。4.2 原位运算结构FFT 的工程实现大多采用“原位计算”方式即输入数据先放到数组里每一级蝶形运算的结果直接写回原数组位置。这样做的好处是只需要一个长度为 N 的复数数组不需要额外分配大块中间存储。DIT 原位运算的规律是每一级蝶形输入的“组间隔”逐级减半。第一级距离为 1第二级距离为 2第三级距离为 4依次类推。以 N8 为例第一级蝶形对是 (0,1)、(2,3)、(4,5)、(6,7)第二级是 (0,2)、(1,3)、(4,6)、(5,7)第三级是 (0,4)、(1,5)、(2,6)、(3,7)。可以看到索引跨度从 1 开始倍增这就是原位蝶形结构的核心索引规律。4.3 旋转因子查找表实际代码中如果每级实时调用 sin/cos开销较大。更常见的方法是预计算一个长度为 N/2 的旋转因子表W_N^k cos(2πk/N) - j·sin(2πk/N)然后用查表替代三角函数计算。在嵌入式平台上还可以用定点数近似旋转因子配合移位和查表来提升性能。5. 基 2 DIT-FFT 实现步骤与代码示例5.1 实现步骤一个完整的基 2 FFT 程序包含三步码位倒序bit-reversal permutation分级蝶形运算旋转因子计算或查表码位倒序的原因在于输入序列按奇偶分组后最终参加蝶形运算的序列顺序不是自然顺序而是原索引二进制位倒过来之后的顺序。例如 N8 时索引 1二进制 001会被放到位置 4二进制 100。实现时通常用一个循环交换数组元素。三级蝶形运算的伪代码结构如下for stage 1 to M len 1 stage half len 1 w_step N stage for start 0 to N-1 step len for j 0 to half-1 k start j p j * w_step W twiddle[p] temp W * X[k half] X[k half] X[k] - temp X[k] X[k] temp注意这里旋转因子的下标是按当前级步长计算而不是直接从 j 取。级数从 1 到 Mlen 从 2 开始倍增half 是当前级蝶形距离的一半。5.2 C 语言实现下面的代码实现 N1024 以内、任意 2 的幂次长度的复数 FFT。这里用 float 复数组交替存放实部和虚部输入输出共用一个数组。#include stdio.h #include math.h #define FFT_N 1024 void bit_reverse(float* re, float* im, int n) { int j 0; for (int i 0; i n - 1; i) { if (i j) { float tr re[i]; re[i] re[j]; re[j] tr; float ti im[i]; im[i] im[j]; im[j] ti; } int k n 1; while (k 0 (j k)) { j ~k; k 1; } j | k; } } void fft(float* re, float* im, int n) { bit_reverse(re, im, n); for (int len 2; len n; len 1) { int half len 1; double angle_step -2.0 * M_PI / len; for (int i 0; i n; i len) { for (int j 0; j half; j) { double angle angle_step * j; float wr (float)cos(angle); float wi (float)sin(angle); int idx i j; int idx2 idx half; float tr re[idx2] * wr - im[idx2] * wi; float ti re[idx2] * wi im[idx2] * wr; re[idx2] re[idx] - tr; im[idx2] im[idx] - ti; re[idx] re[idx] tr; im[idx] im[idx] ti; } } } } int main(void) { float re[FFT_N] {0}; float im[FFT_N] {0}; for (int i 0; i FFT_N; i) { re[i] (float)sin(2.0 * M_PI * 50 * i / FFT_N) 0.5f * (float)sin(2.0 * M_PI * 120 * i / FFT_N); } fft(re, im, FFT_N); for (int k 0; k 8; k) { float mag sqrtf(re[k] * re[k] im[k] * im[k]); printf(k%d mag%.3f\n, k, mag); } return 0; }这段代码没有做归一化输出幅值是 FFT 的原始模值。如果要做工程上的幅值谱一般需要乘以 2/N单边谱直流分量乘 1/N。使用时可以把FFT_N改成需要的点数但必须是 2 的整数次幂。5.3 Python 验证版本写一个基于 NumPy 的验证版本用于对比手写 C 实现的正确性import numpy as np def dit_fft(x): n len(x) if n 1: return x even dit_fft(x[0::2]) odd dit_fft(x[1::2]) factor np.exp(-2j * np.pi * np.arange(n // 2) / n) return np.concatenate([ even factor * odd, even - factor * odd ]) # 测试信号50Hz 120Hz 正弦叠加 fs 1024 t np.arange(fs) / fs x np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 120 * t) X dit_fft(x) mag np.abs(X[:8]) print(前8个频点幅值:, mag) print(参考结果:, np.abs(np.fft.fft(x)[:8]))这个递归版本便于理解原理但因为每次递归都创建新数组实际工程性能不如迭代版本。用迭代版本更节省内存也更容易移植到嵌入式平台。5.4 幅值谱和相位计算FFT 输出是复数通常需要进一步换算幅值mag sqrt(re^2 im^2)相位phase atan2(im, re)功率谱power re^2 im^2单边幅值谱mag_single 2 * mag / N直流分量单独取 mag / N以采样率 fs1024 HzFFT_N1024 为例频率分辨率是 fs / N 1 Hz。第 k 个频点对应的频率是 k · fs / N。如果信号包含 50 Hz 分量就会在 k50 附近出现峰值。相位测量需要注意FFT 得到的是“以窗函数起始点为参考”的相位。如果直接对一段不含整数周期采样长度的信号做 FFT相位结果会受频谱泄漏影响不能直接作为该频率成分的绝对相位。工程上常用整周期采样、加窗插值或 Goertzel 算法来测量特定频率的相位。6. DSP 库与工程应用场景6.1 嵌入式 DSP 库中的 FFT在 ARM Cortex-M 平台C MSIS-DSP 库提供了arm_cfft_f32、arm_rfft_f32等标准接口内部已经实现了码位倒序和蝶形运算用户只需要维护输入输出缓冲区和旋转因子表。典型调用流程#include arm_math.h #define FFT_SIZE 1024 float32_t input[FFT_SIZE * 2]; // 复数格式排列 float32_t magOutput[FFT_SIZE]; void process_fft(float32_t* timeDomain, float32_t* freqMag) { // 将实数时域信号填入复数数组的实部虚部清零 for (int i 0; i FFT_SIZE; i) { input[2 * i] timeDomain[i]; input[2 * i 1] 0.0f; } arm_cfft_f32(arm_cfft_sR_f32_len1024, input, 0, 1); arm_cmplx_mag_f32(input, freqMag, FFT_SIZE); }这里用的是复 FFT 接口实数信号只需要把虚部置 0。如果使用arm_rfft_f32还可以利用实信号的共轭对称性减少一半计算量输出直接就是 N/2 个频率点的幅值。6.2 实时频谱显示实时频谱分析的基本框架是以固定采样率采集时域数据。攒够 N 点数据后做 FFT。更新频谱显示或发送到上位机。比如用 F407 这类 Cortex-M4 芯片做实时频谱显示N1024 的 FFT 在开启 FPU 和 DSP 指令后通常能在毫秒级完成关键瓶颈往往在 ADC 采样率和屏幕刷新率。工程上建议把采样、FFT、显示分成三个模块用双缓冲或环形缓冲避免数据覆盖。6.3 频率测量与相位测量如果要做单频信号的频率和相位测量直接做 FFT 后找最大峰值可以粗测频率但精度受分辨率限制。提高测量精度的常见手段增加 FFT 点数分辨率变为 fs / N但会增大计算量和内存。对频谱峰值进行插值例如双谱线插值或抛物线插值可把频率精度提升到亚分辨率级别。使用 Goertzel 算法只计算关注频率点的 DFT计算量比完整 FFT 小。相位测量更稳妥的方法是先通过锁相或整周期采样使采样窗口覆盖整数个信号周期再对 FFT 峰值处的相位角取 atan2。在非整周期采样下频谱泄漏会让相位结果失真必须加窗和做相位修正。6.4 包络谱与故障诊断在旋转机械故障诊断中FFT 的一类重要应用是包络谱分析。流程通常是采集振动加速度信号。对信号做带通滤波保留与故障相关的共振频带。对滤波后的信号取包络例如使用希尔伯特变换计算解析信号幅值。对包络信号再做 FFT得到包络谱。包络谱里可以看到故障特征频率及其倍频比如轴承外圈故障特征频率、齿轮啮合频率的边带等。这个场景下的 FFT 点数一般取 1024、2048 或 4096采样率根据机器转速和故障特征频率确定。7. 资源占用与性能观察7.1 计算量对比FFT 点数 N直接 DFT 复数乘法量基 2 FFT 复数乘法量256655361024512262144230410241048576512040961677721624576这个表说明点数越大FFT 的加速比越明显。但要注意FFT 的蝶形运算需要访问数组元素随着 N 增大内存访问和缓存 miss 的影响也会增大实际耗时不会严格按理论计算量下降。7.2 内存占用迭代版 FFT 需要至少一个长度为 N 的复数数组按 float 计算占用字节数为 N × 2 × 4 8N 字节。N1024 时约 8 KBN4096 时约 32 KB。如果使用 double则翻倍。旋转因子表额外占用 N/2 个复数也就是 N×4 字节float 复数。单片机上做 FFT 之前先检查 RAM 是否足够。如果 RAM 紧张可以改用 16 位定点 FFT或者在时域窗函数计算完成后再覆盖原始数据数组减少临时变量。7.3 影响耗时的因素FFT 点数点数翻倍蝶形级数增加 1计算量近似线性翻倍再加上一层的开销。数据类型float 比 double 快定点比浮点在无 FPU 的 MCU 上更快。窗函数如果每次都实时计算窗函数系数并做乘加会额外消耗时间建议把窗函数系数预存。旋转因子实时 sin/cos 很慢用查找表可明显提速。编译器优化开启 -O2 或 -O3 后性能差异显著。缓存友好性按顺序访问数组比随机跳转快蝶形运算的索引规律已经相对规整。7.4 判断 FFT 是否耗时的思路在嵌入式平台上可以用 GPIO 翻转引脚测量 FFT 函数耗时也可以用定时器计数。PC 上可以用clock()、chrono或 Python 的time.perf_counter测量。测量时要把 FFT 单独跑多次取平均值排除其他任务干扰。8. 常见问题与排查方法问题现象可能原因排查方式解决方案FFT 输出全是 0 或幅度极小输入信号未赋值、虚部被错误清零、信号本身接近直流打印输入数组前几个值确认输入数据有效直流分量要看第 0 个频点频谱峰值位置不对采样率或频率分辨率计算错误检查 fs/N 是否和预期一致确认采样率和 FFT 点数匹配信号频率应落在频点附近峰值周围有大量泄漏采样点数不是信号周期的整数倍观察主瓣是否展宽加窗函数或调整采样长度使信号整周期截断相位结果跳动非整周期采样或噪声干扰对比加窗后的相位结果增加采样长度使用插值修正或锁相采样嵌入式平台耗时过高使用 double 计算、实时三角函数、未开启优化查看编译选项和 FT 函数内部实现改用 float、旋转因子查表、开启 FPU 和 -O2码位倒序后数据乱序索引交换逻辑写错用 N8 小序列单步调试对照 8 点码位倒序表逐项检查逆变换后波形不对忘了归一化或旋转因子符号反了先做正变换再逆变换验证检查旋转因子符号正变换用 -j逆变换用 j 并除以 N用库函数结果和自己写的结果不一致数据格式或缩放因子不同比较幅值谱形态不比较原始复数按库文档确认输入排列方式和缩放因子对于“峰值位置不对”这类问题最简单的验证方法是用一个已知频率的正弦波做测试。例如 1024 Hz 采样率下给一个 50 Hz 正弦波N1024预期峰值出现在 k50。如果出现在别的位置优先检查采样率和数据拼接逻辑。9. 最佳实践与使用建议9.1 先做小点数验证第一次写 FFT 代码不要直接上 4096 点。先用 N8 或 N16 的简单序列手算出 DFT 结果和程序输出对比。N8 的码位倒序表很短可以手动检查索引对不对。小点数跑通后再逐步放大。9.2 输入信号规范化时间序列在做 FFT 前建议先去除直流分量否则第 0 个频点的幅值会很大影响动态范围。如果只关心交流分量可以先减去均值。加窗可以减少频谱泄漏但会改变幅值需要做幅值恢复系数修正。9.3 频谱分析流程标准化一个可复用的频谱分析流程至少包含数据采集与预处理去除异常点、去直流、必要时滤波。分帧和加窗确定帧长、重叠率、窗函数类型。FFT 与幅值换算计算幅值谱、功率谱或包络谱。特征提取峰值频率、峰值幅值、边带间隔、故障特征频率匹配。结果输出保存频谱数据或绘制曲线。把这五步封装成函数后续换信号只需要改输入数据和参数。9.4 数据管理的建议原始时域数据、窗函数系数、旋转因子表、FFT 输出结果分目录或分变量存放。批量处理多个文件时输出频谱文件命名建议包含采集时间、采样率、点数等参数方便回溯。保存频谱结果时除了幅值谱尽量同时保存原始采样频率和 FFT 点数否则后续无法把频点换算成物理频率。9.5 关于库的选择纯学习或验证算法用 Python NumPy/SciPy 最方便。PC 端高性能需求可以考虑 FFTW 或 PocketFFT。ARM 嵌入式平台优先用 CMSIS-DSP 的标准库函数。需要商用发布注意检查所选库的开源许可证。CMSIS-DSP 对大部分产品友好但具体场景仍需确认授权范围。9.6 合规提醒FFT 本身只是数学工具但采集和分析的信号可能涉及隐私、版权、商业机密。例如采集设备振动数据、分析语音信号或处理通信信号时要确保已获得合法授权并且只在授权范围内使用分析结果。涉及个人音频、图像或设备内部信息的数据不要随意公开或传播。10. 总结与下一步FFT 的核心价值就一句话用 O(N log2 N) 的计算量逼近 O(N^2) 的频谱结果。蝶形运算作为 FFT 的基本单元结构固定关键在掌握每一级的索引间隔和旋转因子变化规律。基 2 时间抽取 FFT 的实现链路是码位倒序、分级蝶形、旋转因子查表跑通这三个环节对 FFT 的理解就不再停留在公式层面。建议先从 8 点或 16 点的小序列开始手写一版 C 或 Python 代码再用自己熟悉的数据比如 50 Hz 正弦波做验证。确认频谱峰值位置正确后再扩展到 1024 点或 4096 点并逐步加上窗函数、幅值修正、包络谱和相位测量。后续可以继续深入的方向有基 4 FFT 与混合基 FFT、实数序列 FFT 的优化、定点 FFT 实现、加窗插值算法、Goertzel 算法在单频测量中的应用以及把 FFT 封装成实时流式频谱分析模块。如果这篇文章帮你在 FFT 的知识点上省了时间建议收藏备用编码过程中遇到蝶形结构问题随时回来对照。