这次我们来看数字信号处理里出场率最高的一个算法结构FFT 的蝶形运算。网上的 FFT 原理文章很多但多数停留在“FFT 比 DFT 快很多”这种结论上。真正往下走就会碰到几个问题蝶形图怎么看、旋转因子怎么算、什么是倒位序、为什么能原位计算、STM32 上调用 DSP 库做实时频谱要注意什么。这篇就把 FFT 的结构一层层拆开重点讲蝶形运算的两种典型结构并给出可运行的代码示例和排查思路。文章适合这几类读者正在学数字信号处理的大学生需要在 MCU、实时系统里实现频谱分析的嵌入式开发者用 MATLAB/Python 做振动分析、包络谱、故障诊断但总搞不清 FFT 参数怎么选的人。看完之后你应该能独立画出 8 点 FFT 蝶形图理解旋转因子的查表逻辑并搞清楚为什么很多 FFT 例程都要先做位序反转。1. FFT 蝶形运算核心能力速览先给一张速览表把整个 FFT 蝶形结构方案的关键点列出来。能力项说明算法核心利用旋转因子的周期性和对称性把 N 点 DFT 拆解为多级蝶形运算时间复杂度直接 DFT 为 O(N^2)FFT 为 O(N log2 N)典型点数2 的幂次如 64、128、256、512、1024、4096常用结构按时间抽取 DIT-FFT、按频率抽取 DIF-FFT关键预处理输入倒位序Bit-Reversal或输出倒位序是否支持原位运算是同一组缓冲区可直接覆盖中间结果旋转因子W_N^k e^(-j2πk/N)工程上常用查表法预存典型应用实时频谱分析、振动包络谱、相位测量、卷积快速计算常见误用问题未加窗导致频谱泄漏、采样率不匹配导致镜像混叠、相位未校准这里先说一个关键判断无论网上代码怎么变FFT 的蝶形运算结构都逃不出“分成上下两路一路直接加一路乘旋转因子后相加减”这个基本单元。把这一条刻在脑子里后面所有内容都好理解。2. FFT 要解决什么问题从 DFT 复杂度说起离散傅里叶变换的定义是X[k] Σ x[n] * e^(-j2πkn/N)n0 到 N-1如果直接按这个公式算每个频点 k 都要做 N 次复乘和 N 次复加N 个频点就是 N^2 次复数乘法。N 取 1024 时计算量超过一百万次复乘N 取 4096 时直接算就非常吃力。对嵌入式场景来说这种复杂度很难支撑实时频谱显示。FFT 的核心思路是拆解。观察旋转因子 W_N^k e^(-j2πk/N)可以发现两个性质周期性W_N^(kN) W_N^k对称性W_N^(kN/2) -W_N^k这两个性质让很多重复计算可以合并。蝶形运算就是把一个大的 DFT 逐级拆成一半一半的小 DFT每一级用“加/减一个旋转因子乘积”来合并结果。最终的计算量约为 (N/2) * log2(N) 次复数乘法相比直接 DFT 是数量级的下降。所以 FFT 并不是一种全新的变换它只是 DFT 的高效计算方式。蝶形运算结构就是“高效”二字的落地形态。3. 蝶形运算的基本单元蝶形运算的名字来自数据流图的外形两条输入线从左边进来交叉合并后从右边出去看起来像蝴蝶翅膀。基本的 DIT 蝶形单元可以表示成节点输入输出上支路AA W_N^k * B下支路BA - W_N^k * B用公式写就是A A W * B B A - W * B其中 W 是旋转因子A、B 是复数数据。整个蝶形只需要一次复数乘法、两次复数加减法。从硬件或者 MCU 的角度看这个结构非常友好乘法器只用一个加减法器需要两个数据和旋转因子可以分别存储非常适合流水线实现。从软件的角度看这个结构支持原地计算A 和 B 可以直接写回 A 和 B 原来的内存位置不影响后续计算这就是“原位 FFT”的原理。蝶形单元看起来简单但多个蝶形单元组合起来时不同抽取方式会得到不同结构。下面重点看两种主流结构。4. 按时间抽取的 DIT-FFT 蝶形结构按时间抽取是把输入序列 x[n] 按照下标奇偶分成两组。偶数点组成一个 N/2 点 DFT奇数点组成另一个 N/2 点 DFT然后通过旋转因子合并。以 8 点 FFT 为例结构分三级第一级4 个 2 点蝶形第二级2 个 4 点蝶形第三级1 个 8 点蝶形每一级都有明确的旋转因子索引。8 点 FFT 的第三级旋转因子为 W_8^0、W_8^1、W_8^2、W_8^3也就是四个不同的复数旋转因子前面级数用的因子会更少。DIT-FFT 在工程中常见的处理流程是输入序列 ↓ 倒位序重排 ↓ 第 1 级蝶形间隔为 1 第 2 级蝶形间隔为 2 第 3 级蝶形间隔为 4 ↓ 输出自然顺序 X[k]倒位序的意思是把输入下标写成二进制后反转比如 8 点时原始下标二进制倒序后下标0000010014201023011641001510156110371117这就是 FFT 例程里常见bit_reverse函数做的事情。如果输入不是自然顺序而是先按倒位序排好那么每一级蝶形计算完成后输出就是自然顺序 X[0]、X[1]、X[2]……DIT-FFT 的重点理解方式把时间序列不断二分区最终每一级只处理相邻数据点之间的蝶形组合。程序实现时级数用stage表示蝶形跨度用len控制旋转因子索引用当前蝶形位置和级数联合计算。5. 按频率抽取的 DIF-FFT 蝶形结构DIF-FFT 正好和 DIT-FFT 相反。它先把序列分成前后两半先做蝶形运算再抽取偶数和奇数输出。DIF 的基本蝶形是A A B B (A - B) * W也就是说加法和减法在第一阶段完成旋转因子乘在下支路上。到了下一级整个序列被分成两个独立的子变换。DIF-FFT 的流程是输入自然顺序 x[n] ↓ 第 1 级蝶形加/减后乘旋转因子 第 2 级蝶形 第 3 级蝶形 ↓ 输出为倒位序DIF-FFT 适合硬件实现因为它不需要在输入端做倒位序可以先做流水计算最后统一倒位序输出。很多 DSP 库的内部实现会混合使用这两种结构比如某些四级拆解中前几级用 DIF后几级用 DIT用来减少旋转因子的存储量。实际工程中直接手写 FFT 时优先用 DIT-FFT因为输入倒位序之后输出顺序和频谱分析逻辑直接对应调试方便。用现成 DSP 库时不用关心内部是 DIT 还是 DIF只要知道库函数是否要求时间域数据做特殊排列。6. 原位计算与内存布局蝶形运算的一个典型优势是原位计算。因为每个蝶形只读取两个复数计算出两个新复数后可以直接覆盖原位置不需要另外开辟大型中间缓冲区。工程实现时内存布局一般是// 复数实部和虚部交错存储 float32_t fft_buffer[2 * FFT_LENGTH]; // 0: 实部0 1: 虚部0 2: 实部1 3: 虚部1如果是纯软件实现且用 float 数组可以把实部数组和虚部数组分开也可以交错存放。嵌入式常用交错存放因为 DSP 库的复数函数多数按这种布局设计。原位计算带来的好处是省内存但前提是所有中间结果不能提前覆盖。这就是为什么蝶形运算必须以“级”为单位推进不能在一个级内跨蝶形乱写。写到代码里就是外层循环按 stage 控制内层循环按同一级不同蝶形组推进。一个容易出错的点如果旋转因子索引算错原位覆盖后数据会被污染而且很难追查。排查时建议先用 N4 或 N8 的小点数把每级数据打印出来和 MATLAB/Python 的结果对比。7. 旋转因子的工程处理与查表法旋转因子 W_N^k 是一个复数W_N^k cos(2πk/N) - j * sin(2πk/N)直接实时计算三角函数效率太低尤其是嵌入式平台。工程上常用查表法预先算好长度为 N/2 的旋转因子表运行时直接取用。旋转因子查表的关键是索引映射。不同的 FFT 实现索引算法不同。以 DIT-FFT 为例第 stage 级共有2^(stage-1)个不同旋转因子蝶形跨度是2^stage。旋转因子索引通常是// 假设 N 为 2 的幂 // stage 从 1 开始M log2(N) // 当前蝶形在组内的位置为 jk 为旋转因子索引 int k j * (N stage);比如在某个级内蝶形跨度是 8 时旋转因子可能按W_N^(0 * N/8)、W_N^(1 * N/8)这样的方式递增。这个推导不复杂但很容易算错。建议先画一张完整的 8 点蝶形图把每个蝶形的旋转因子标出来再对照程序理解索引逻辑。旋转因子表的存储也有讲究。如果做的是 4096 点实数 FFT完整复数因子表大约是 2048 个复数float32 下约 16 KB。对 STM32F4 这种 RAM 在 128 KB 左右的芯片可以接受。但如果是低端 MCU就要考虑用半表只存 0 到 π/4 范围内的角度再用象限转换减少存储。8. FFT 蝶形运算编程实现与验证8.1 Python 递归 FFT 实现先从简单版本入手。Python 代码适合验证蝶形运算结构是否正确import numpy as np def fft_recursive(x): N len(x) if N 1: return np.asarray(x, dtypecomplex) even fft_recursive(x[0::2]) odd fft_recursive(x[1::2]) theta -2.0 * np.pi / N w np.exp(1j * theta * np.arange(N // 2)) T w * odd return np.concatenate([even T, even - T]) # 测试 fs 8000 N 512 t np.arange(N) / fs x np.sin(2 * np.pi * 1000 * t) 0.5 * np.sin(2 * np.pi * 2500 * t) X fft_recursive(x) freq np.arange(N) * fs / N mag np.abs(X) # 找前两个最大峰值 idx np.argsort(mag)[::-1][:4] print(峰值的 bin 索引:, idx) print(对应频率:, np.sort(freq[idx]))运行后应该能看到 1000 Hz 和 2500 Hz 对应的频率分量。如果不是这两个频率说明蝶形结构或因子索引有误。这个递归版本原理清晰但工程性能一般。它没有显式倒位序是因为递归天然完成了奇偶分解的顺序调整。理解递归版之后再看迭代版会更清楚。8.2 Python 迭代 FFT 实现迭代版更能体现蝶形运算的分级结构import numpy as np def bit_reverse_index(n, bits): return int(f{n:0{bits}b}[::-1], 2) def fft_iterative(x): N len(x) bits N.bit_length() - 1 # 输入倒位序 x [complex(x[bit_reverse_index(i, bits)]) for i in range(N)] # 多级蝶形计算 length 2 while length N: half length // 2 theta -2.0 * np.pi / length w np.exp(1j * theta * np.arange(half)) for start in range(0, N, length): for j in range(half): u x[start j] v w[j] * x[start j half] x[start j] u v x[start j half] u - v length * 2 return np.array(x) N 8 x np.array([1.0, 0.5, 0.0, 0.2, 0.1, 0.0, 0.0, 0.0]) X_ref np.fft.fft(x) X_mine fft_iterative(x) print(自实现 FFT:, np.round(np.abs(X_mine), 6)) print(NumPy FFT:, np.round(np.abs(X_ref), 6))如果两个结果一致说明 DIT-FFT 蝶形计算和倒位序都没问题。之后可以把这个流程迁移到 C 语言。8.3 C 语言 8 点 FFT 蝶形伪代码实际嵌入式中不会所有场景都用完整 DSP 库有时需要手写一个固定点数 FFT比如 64 点或 128 点。下面给出一个结构化模板#include math.h #include stdint.h #define FFT_N 8 #define FFT_M 3 // log2(FFT_N) typedef struct { float re; float im; } complex_t; // 旋转因子表长度为 FFT_N / 2 complex_t twiddle[FFT_N / 2]; void twiddle_init() { for (int k 0; k FFT_N / 2; k) { twiddle[k].re cosf(2.0f * M_PI * k / FFT_N); twiddle[k].im -sinf(2.0f * M_PI * k / FFT_N); } } uint8_t bit_reverse(uint8_t x, uint8_t bits) { uint8_t y 0; for (uint8_t i 0; i bits; i) { y (y 1) | (x 1); x 1; } return y; } void fft_radix2(complex_t *data, uint8_t N, uint8_t M) { // 1. 倒位序 for (uint8_t i 0; i N; i) { uint8_t j bit_reverse(i, M); if (i j) { complex_t tmp data[i]; data[i] data[j]; data[j] tmp; } } // 2. 多级蝶形 for (uint8_t stage 1; stage M; stage) { uint8_t length 1 stage; // 当前级蝶形跨度 uint8_t half length 1; for (uint8_t start 0; start N; start length) { for (uint8_t j 0; j half; j) { // 旋转因子索引 uint8_t k j * (N stage); complex_t w twiddle[k]; complex_t u data[start j]; complex_t v; // v w * data[start j half] v.re w.re * data[start j half].re - w.im * data[start j half].im; v.im w.re * data[start j half].im w.im * data[start j half].re; data[start j].re u.re v.re; data[start j].im u.im v.im; data[start j half].re u.re - v.re; data[start j half].im u.im - v.im; } } } }这段代码直接对应前面讲的 DIT-FFT 结构。注意旋转因子表用的是 N 点全表长度所以twiddle[k]可以直接索引如果做半表优化索引逻辑要重写。验证时可以给data输入一个已知序列比如 x [1, 0, 0, 0, 0, 0, 0, 0]输出应当是全 1 的直流频谱。再用一个正弦序列测试观察峰值频率是否正确。9. 嵌入式 DSP 库实时频谱测量很多实时频谱项目使用 M4/M7 内核单片机和 CMSIS-DSP 数学库比如常见的arm_rfft_fast_f32函数。先看接口形式#include arm_math.h #define FFT_LENGTH 1024 float32_t input[FFT_LENGTH]; float32_t fft_output[2 * FFT_LENGTH]; float32_t magnitude[FFT_LENGTH / 2]; arm_rfft_fast_instance_f32 fft_inst; void dsp_fft_init(void) { arm_rfft_fast_init_f32(fft_inst, FFT_LENGTH); } void dsp_fft_run(float32_t *time_data) { // 0 表示正变换1 表示反变换 arm_rfft_fast_f32(fft_inst, time_data, fft_output, 0); // 计算幅度谱 arm_cmplx_mag_f32(fft_output, magnitude, FFT_LENGTH / 2); }使用 DSP 库时要注意几点arm_rfft_fast_f32要求输入长度为 2 的幂。输出是复数形式实部和虚部交错排列。幅度谱只取前半部分也就是 0 到 Nyquist 频率。如果做的是 ADC 采样数据直接喂进arm_rfft_fast_f32之前不需要手动倒位序库函数内部已经处理。幅值校准和频率分辨率是嵌入式项目最容易被问到的点。频率分辨率是fs / FFT_LENGTH比如采样率 8000 HzFFT 长度 1024分辨率约 7.8125 Hz。这个数值决定了你能区分多近的两个频率峰。如果测振动信号采样率通常远高于轴承特征频率。预测性维护场景中常用包络谱先让信号通过带通滤波器和希尔伯特变换求出包络波形再对包络做 FFT 得到包络谱用于提取轴承故障特征频率。FFT 蝶形结构本身不关心输入是原始加速度信号还是包络信号但采样参数和分析流程差异会影响最终频率成分解释。10. FFT 常见问题与排查方法下面是工程中最高频的几类问题整理成排查表问题现象可能原因排查方式解决方案频谱出现多余对称峰输入了实数信号但频率轴计算错误检查频率轴是否到 fs/2实数 FFT 只取前 N/2 个频点峰值频率偏移幅度偏低非整周期采样频谱泄漏查看信号频率是否精确落在 bin 上加 Hann/Hamming 窗或改用整周期采样低频附近有明显干扰直流分量过大或加窗后旁瓣泄漏查看 0 Hz 处幅度先做去直流处理再计算 FFT相位结果不稳定未对 FFT 起点做同步触发用同一触发源采集引入同步采样或用固定相位参考点旋转因子索引越界级数 stage 与 N 关系换算错误打印各级 k 的取值重新推导 k j * (N stage)输出幅度只有理论一半单边谱未补因子 2对比频谱总能量除直流和 Nyquist 外幅度乘 24000 点 FFT 结果奇怪点数不是 2 的幂检查 N 是否是 2^M改成 4096或补零到 4096MCU 上计算时间波动大实时计算三角函数检查是否调用 sin/cos换旋转因子查表法频谱出现镜像频率采样率低于信号最高频率两倍检查前级抗混叠滤波提高采样率或增加抗混叠低通滤波最容易踩的坑是频谱泄漏。ADC 连续采样时不注意采样点数与信号周期的关系FFT 结果会出现频率分量扩散。处理方法是加窗函数工程上最常用的是 Hann 窗for (int i 0; i FFT_LENGTH; i) { float window 0.5f * (1.0f - cosf(2.0f * M_PI * i / (FFT_LENGTH - 1))); input[i] adc_value[i] * window; }加窗会降低频谱分辨率、增大主瓣宽度但能有效抑制旁瓣泄漏。做包络谱分析时窗函数的选择需要结合特征频率间隔来权衡。11. FFT 蝶形运算最佳实践11.1 先用小点数验证算法初次实现 FFT 时不要直接上 1024 点或 4096 点。先用 N8 或 N16 的固定序列和 MATLAB 或 Python 的结果逐点对比。对比时同时比较实部和虚部而不是只比较幅度。相位信息错乱往往在幅度上不明显但会在后续测量相位时暴露。11.2 旋转因子统一用查表无论用什么平台建议把旋转因子表提前算好单独放在一个数组里。这样每级蝶形只需要查表不需要调用三角函数。查表索引是 FFT 代码里最容易出错的部分一定要拿小点数逐级打印验证。11.3 数据类型按平台选择没有 FPU 的 MCU用定点 Q15 或 Q31 实现CMSIS-DSP 提供arm_cfft_q15、arm_cfft_q31。带 FPU 的 M4/M7直接用arm_rfft_fast_f32。PC 端验证优先用 Python NumPy 或 MATLAB。定点 FFT 的主要问题是溢出。蝶形计算中间值可能超过输入幅度范围通常需要对每一级或每一帧做缩放。CMSIS-DSP 的定点 FFT 函数内部已经处理了部分缩放但用户仍需关注输入信号幅值范围避免 ADC 满量程输入导致饱和。11.4 批量帧处理时管理好缓冲区实时频谱系统通常有固定采样帧长。FFT 输入缓冲区、输出复频谱、幅度谱、窗函数表建议都定义为静态数组避免在中断里动态分配内存。中断处理函数里只做数据拷贝FFT 计算放到主循环或低优先级线程。11.5 记录算法版本和采样参数FFT 的调试难点在于参数组合太多。建议在每个工程里保留一个配置文件记录采样率、FFT 点数、窗函数类型、加窗是否开启、单边谱是否补因子 2、频率轴换算公式。这样换人排查或换 MCU 移植时能快速定位问题。12. 总结与下一步FFT 蝶形运算的基本结构并不复杂核心就是一条每级计算只有“乘旋转因子 相加 相减”三步。关键是理解和控制好三件事旋转因子怎么查表、倒位序在哪里做、每一级蝶形的跨度和索引怎么算。如果是从零开始学建议按这个顺序推进先看 8 点 DIT-FFT 蝶形图再手写一个小点数迭代版 FFT然后用正弦信号验证频点和相位最后再换实时数据。直接拿现成库虽然简单但遇到“峰值对不上”“相位不对”这类问题时不理解蝶形结构会无从下手。下一步可以继续深入的方向包括实数 FFT 的优化、4 基 FFT 与混合基 FFT、加窗与频谱插值、STFT 短时傅里叶变换、以及基于 Hilbert 变换的包络谱分析。FFT 不是终点但蝶形运算结构是所有频域分析工具的地基值得花时间彻底吃透。