很多做信号处理的朋友都有过这种体验numpy.fft.fft一行代码就能算出频谱STM32 里调用 DSP 库的arm_cfft_f32也能跑出结果但一旦遇到下面这些问题就很容易卡住为什么 FFT 的输出点数必须是 2 的整数次幂为什么计算出的频率轴第 0 点是直流分量第 N/2 点是奈奎斯特频率为什么加窗之后频谱会“变胖”不加窗又会“泄漏”为什么用单片机做实时频谱显示时总是担心 CPU 扛不住这些问题表面上是“FFT 使用问题”根子上却是一个东西你没有真正理解蝶形运算结构。蝶形运算是 FFT 之所以“快速”的核心也是理解频域分辨率、计算量、内存排布和实时性的钥匙。如果你只停留在调库层面遇到点数、相位、窗函数、资源占用这类问题时就只能靠试错。这篇文章围绕“数字信号处理 FFT 快速傅里叶变换蝶形运算结构方法”这条主线把 DFT 到 FFT 的演进逻辑、基 2 时间抽取DIT和频率抽取DIF两种经典结构讲透然后给出可运行的 Python 实现再结合频谱分析、相位测量、包络谱分析和 MCU 实时频谱显示这些真实场景说明蝶形结构在实际项目中到底怎么用、有哪些坑。读完你可以得到三样东西第一彻底看懂教材里那张蝴蝶交叉的信号流图第二手写一个可用的基 2 FFT而不是只会调fft函数第三在实际项目里知道自己该选多少点、加什么窗、怎么验证结果。1. 为什么 DFT 很“重”而 FFT 很“轻”先回到最原始的问题离散傅里叶变换DFT到底在算什么对一段长度为 N 的离散序列 x[n]它的 N 点 DFT 定义是X[k] sum_{n0}^{N-1} x[n] * W_N^{kn}, k 0,1,...,N-1其中旋转因子W_N^{kn} e^{-j 2π kn / N}这个公式的信息量很大。它说明DFT 的本质是把一个长度为 N 的序列分解成 N 个不同频率的复指数序列的加权和。X[0] 对应直流分量X[1] 对应基波X[k] 对应第 k 个频率分量。每一个 X[k] 的计算都需要对 N 个输入点做一次复乘和累加。所以完成全部 N 个 X[k] 的计算需要复数乘法N² 次复数加法N(N-1) 次约等于 N² 次当 N1024 时N² 1,048,576也就是约一百万次复数乘法。当 N4096 时乘法次数就变成了约 1677 万次。虽然现代 CPU 跑这个量级并不吃力但放到实时信号处理里每毫秒都要来一轮压力就完全不同了。更何况很多 MCU 根本没有硬件浮点单元一次复数乘法要拆成不知道多少条整数指令。1965 年 Cooley 和 Tukey 发表的 FFT 算法把计算量从 O(N²) 降到了 O(N log₂N)。同样是 1024 点FFT 只需要约 1024 × 10 10240 次复数运算相比直接 DFT 差了大约两个数量级。这个差距才是 FFT 能在实时系统、嵌入式设备、工业检测领域大规模落地的根本原因。2. 蝶形运算FFT 的最小组成单元要理解 FFT 为什么快必须先理解“蝶形运算”。蝶形运算是一个只有两个输入、两个输出的计算单元它的形式非常固定输入a, b 输出 X a b * W Y a - b * W其中 W 是旋转因子。之所以叫“蝶形”是因为把输入输出画成信号流图时两条支路交叉再合并形状像一只蝴蝶。蝶形运算的核心价值在于它把一个二点 DFT 的计算变成了“一次复数乘法 两次复数加法”并且在整个 FFT 的每一级里反复复用。FFT 的算法结构本质上就是把一个大 DFT 逐级拆成许多个小 DFT而每一级的小 DFT 最终都变成若干个蝶形运算。这里要特别强调一点蝶形运算里的加法和减法共享同一个乘积项 b * W。也就是说计算 a bW 和 a - bW 时复数乘法只需要做一次。这种“乘法结果复用”的思想是 FFT 能减少计算量的一个重要原因。如果你手动写过朴素 DFT会发现每算一个点都独立做 N 次乘法而 FFT 通过蝴蝶结构把中间结果反复复用避免了大重复计算。3. 基 2 时间抽取 FFTDIT-FFT的完整推导基 2 FFT 的“基 2”意思是每一级都把序列按奇偶位置分成两组直到每组只剩 2 个点。这里以最常用的按时间抽取Decimation-In-Time, DIT为例逐步拆解。3.1 第一步奇偶分离设序列长度为 N且 N 2^M。把 x[n] 分成偶数下标和奇数下标两组偶数部分x[0], x[2], x[4], ..., x[N-2] 奇数部分x[1], x[3], x[5], ..., x[N-1]代入 DFT 定义经过变形可以得到X[k] X_even[k] W_N^k * X_odd[k] X[k N/2] X_even[k] - W_N^k * X_odd[k]其中 X_even[k] 是偶数下标序列的 N/2 点 DFTX_odd[k] 是奇数下标序列的 N/2 点 DFT。这里有两个关键观察用 N/2 点 DFT 的结果可以通过加减组合出完整的 N 点 DFT。这一个“减半”就让计算量从 N² 变成了大约 2 × (N/2)² N/2是第一步优化。X[k] 和 X[kN/2] 之间恰好是一个标准的蝶形运算结构。这正是蝶形运算出现的自然来源。3.2 第二步继续对半拆分把 N/2 点的 DFT 再按奇偶拆一次变成 N/4 点的 DFT。每拆一次系统规模减半一直拆到 2 点 DFT。2 点 DFT 的公式非常简单X[0] x[0] x[1] X[1] x[0] - x[1]这正好是旋转因子为 W_2^0 1 的蝶形运算。也就是说2 点 DFT 天然就是一个蝶形。3.3 第三步码位倒序从数学推导看DIT-FFT 是“输入倒序输出顺序”。也就是说输入序列要先按照二进制位倒序重新排列然后才能逐级做蝶形运算。N8 时的输入顺序变化如下原始索引二进制位倒序重排后索引0000000010011004201001023011110641000011510110156110011371111117所以 DIT-FFT 的流程是清晰的输入序列 → 码位倒序 → 第 1 级蝶形 → 第 2 级蝶形 → ... → 第 M 级蝶形 → 输出频谱每一级蝶形的旋转因子和跨度都不一样这是初学时最容易搞混的地方。可以从 N8 的例子看规律。3.4 N8 的完整蝶形结构N8 时共有 M log₂8 3 级蝶形第 1 级蝶形跨度为 1旋转因子为 W_2^0 1第 2 级蝶形跨度为 2旋转因子为 W_4^0 1, W_4^1 -j第 3 级蝶形跨度为 4旋转因子为 W_8^0, W_8^1, W_8^2, W_8^3可以用文字描述这个结构第 3 级最关键它把前两级得到的两组 4 点结果通过 4 个蝶形组合成最终的 8 点频谱。每一级蝶形数量都是 N/24 个所以 N8 的 FFT 总共有 M × N/2 3 × 4 12 个蝶形运算。推广到一般情况N 点基 2 FFT 的蝶形总数是(N/2) * log₂N每次蝶形包含 1 次复数乘法和 2 次复数加法所以总计算量大约为复数乘法(N/2) * log₂N复数加法N * log₂N这正是 FFT 计算量 O(N log₂N) 的来源。对比 DFT 的 N²差距一目了然。4. 频率抽取 FFTDIF-FFT与 DIT-FFT 的对比除了按时间抽取的 DIT-FFT还有一种同等重要的结构按频率抽取Decimation-In-Frequency, DIF-FFT。它和 DIT-FFT 在数学上是等价的但信号流图的方向正好相反。DIF-FFT 的核心思路不是把输入按奇偶分组而是把输出按前后两半拆分。它的特点是输入是自然顺序不需要码位倒序输出是倒序的需要在使用前重新排列蝶形运算的旋转因子位置和 DIT-FFT 不同。两种结构的对比可以总结成下表对比维度DIT-FFTDIF-FFT输入顺序码位倒序自然顺序输出顺序自然顺序码位倒序旋转因子位置在蝶形输入支路在蝶形输出支路实现难度需要先做倒序需要最后做倒序使用场景通用软件实现、教材示例某些硬件流水线结构、定点实现在实际工程中DIT-FFT 更常见因为多数库函数选择“输入倒序、输出顺序”的约定。但 DSP 芯片和 FPGA 的某些 FFT IP 核可能采用 DIF 结构因为它的数据流更适合流水线。一个很关键的理解是DIT 和 DIF 不是两种不同的数学变换而是同一变换的两种计算调度方式。它们算出来的结果完全一致只是中间过程的数据流向和旋转因子的位置不同。5. 手写 Python 实现基 2 DIT-FFT 完整代码理论讲再多不如跑一段代码。下面用 Python 手写一个基 2 DIT-FFT不使用numpy.fft只借助numpy做复数数组操作。先实现码位倒序import numpy as np def bit_reverse(x): 对输入序列进行码位倒序重排。 要求 len(x) 必须是 2 的整数次幂。 n len(x) if n (n - 1) ! 0: raise ValueError(长度必须为 2 的整数次幂) 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] return x这段代码是经典的位反转算法逻辑就是用二进制位翻转的方式找到交换位置。不用记忆它直接当工具函数用即可。然后实现核心的蝶形 FFTdef fft_dit(x): 基 2 按时间抽取 FFTDIT-FFT。 x 为一维复数数组长度必须为 2 的整数次幂。 返回长度相同的复数频谱数组。 n len(x) if n (n - 1) ! 0: raise ValueError(长度必须为 2 的整数次幂) # 复制输入避免修改原数组 a np.array(x, dtypenp.complex128) a bit_reverse(a) # 开始逐级蝶形运算 size 2 while size n: half size // 2 # 当前级的旋转因子 w np.exp(-2j * np.pi * np.arange(half) / size) for start in range(0, n, size): for k in range(half): u a[start k] v a[start k half] * w[k] a[start k] u v a[start k half] u - v size 1 return a这段代码在结构上非常贴近教材里的信号流图最外层循环控制级数size 从 2 开始每轮翻倍中间层循环遍历每一组蝶形最内层循环计算同一组内的 half 个蝶形运算。运行验证一下# 生成一个包含 50Hz 和 120Hz 分量的测试信号 fs 1000 # 采样率 1000 Hz N 1024 # FFT 点数 t np.arange(N) / fs x 0.7 * np.sin(2 * np.pi * 50 * t) np.sin(2 * np.pi * 120 * t) # 计算 FFT X fft_dit(x) # 计算幅度谱 mag np.abs(X) / (N / 2) freq np.arange(N) * fs / N # 打印能量最大的 3 个频率 top_indices np.argsort(mag)[::-1][:5] for idx in top_indices[:3]: print(f频率: {freq[idx]:.2f} Hz, 幅度: {mag[idx]:.4f})预期输出类似频率: 120.00 Hz, 幅度: 1.0000 频率: 50.00 Hz, 幅度: 0.7000这说明手写 FFT 和理论完全一致。为了进一步确认还可以和numpy.fft.fft的结果对比X_np np.fft.fft(x) print(最大误差:, np.max(np.abs(X - X_np)))如果实现正确最大误差应该是 1e-12 量级来自浮点计算误差而不是一个不可忽略的大数。6. 从 Python 到嵌入式实际项目里的 FFT 怎么用理解了蝶形结构之后再回到开头提到的那些实际场景你会发现很多“工程问题”其实都是“结构问题”。6.1 参数确定采样率、点数和频率分辨率在实现一个 FFT 分析系统之前必须先明确三个参数的关系频率分辨率 Δf fs / N其中 fs 是采样率N 是 FFT 点数。如果你要区分两个频率相差 1Hz 的信号而采样率是 1000Hz那么 N 至少需要 1000。为了满足基 2 FFT工程上通常会取 1024。这就是为什么 FFT 点数往往是 1024、2048、4096 这些值——它们不是随意定的而是由“频率分辨率 → 最小 N → 向上取 2 的整数次幂”决定的。6.2 用 STM32 做实时频谱显示时的资源估算在 STM32F407 这类 MCU 上做实时频谱显示通常使用 CMSIS-DSP 库的arm_cfft_f32函数。核心调用方式大致如下#include arm_math.h #define FFT_SIZE 1024 arm_cfft_instance_f32 fft_instance; float32_t input[FFT_SIZE * 2]; // 实部虚部交错存放 float32_t output[FFT_SIZE]; void init_fft(void) { // 初始化 FFT 实例点数 1024 arm_cfft_init_f32(fft_instance, FFT_SIZE); } void compute_spectrum(void) { // input 中按 实部,虚部,实部,虚部 的顺序填充采样数据 // 实数数据虚部填 0 for (int i 0; i FFT_SIZE; i) { input[2 * i] sample_buffer[i]; input[2 * i 1] 0.0f; } // 执行 FFT第三个参数 0 表示正变换 arm_cfft_f32(fft_instance, input, 0, 1); // 计算幅度谱sqrt(real^2 imag^2) arm_cmplx_mag_f32(input, output, FFT_SIZE); // output 数组中即可得到幅度谱通常只取前 FFT_SIZE/2 个点 }这里的arm_cfft_f32定点/浮点实现内部就是按蝶形结构组织的。你在代码里只需要调用函数但理解蝶形结构能帮你回答几个实际问题为什么input要按“实部虚部交错”存放因为 SIMD 优化和内存访存的连续性要求。为什么输出只需要取前 FFT_SIZE/2 个点因为实信号的频谱关于奈奎斯特频率共轭对称后一半没有额外信息。为什么 FFT 点数越大CPU 负载越高因为蝶形总数为 (N/2) * log₂N点数翻倍蝶形数量不是翻倍而是约等于翻倍加一倍的级数运算量具体可以精确计算。从资源估算来看1024 点浮点 FFT 在 STM32F407 上通常只需要几毫秒但如果你同时还要做 LCD 刷新、ADC 采样、通信协议解析就需要仔细分配时间片。更稳妥的做法是在工程上先测量一次 FFT 的实际耗时再用定时器把采样、计算、显示切成不同的任务。6.3 包络谱分析与希尔伯特变换在预测性维护、轴承故障诊断这类场景中直接对原始振动信号做 FFT 往往看不到明显的故障特征频率。这是因为故障冲击信号是宽频带的能量分散。工程上常用的做法是先做希尔伯特变换得到解析信号再取包络最后对包络做 FFT得到包络谱。整个过程可以理解为原始振动信号 → 希尔伯特变换 → 解析信号 → 取包络 → 包络信号 → FFT → 包络谱包络谱能突出周期性冲击的特征频率配合蝶形 FFT 可以快速识别轴承故障。这也是“数字信号处理 FFT 包络谱分析”在工业预测性维护中越来越常见的原因。6.4 相位测量FFT 不仅能测幅值还能测相位。每个频点的相位可以通过复数的虚部和实部计算phase[k] arctan(imag(X[k]) / real(X[k]))在实际工程中比如用 DSP 库测量两个同频信号的相位差通常会先做 FFT找到幅度最大的频点再提取该频点的相位然后对两路信号做差。这里有一个容易踩的坑FFT 的输入起始位置会直接影响相位结果。如果采样窗口的起点不同测出来的“绝对相位”就不同但两路信号的相位差依然是稳定的。所以做相位测量时要确认两路信号使用的是同一段采样时间基准。7. 常见问题与排查方法问题现象可能原因排查方式解决方案频谱峰值位置不对采样率、FFT 点数、频率轴坐标对不上检查 fs/N 的分辨率打印频率轴确认频率轴为 arange(N) * fs / N频谱泄漏严重峰值“变胖”信号频率不是频率分辨率的整数倍用非整周期采样信号做测试加窗函数如汉宁窗、布莱克曼窗幅值偏小或偏大归一化方式不对对比已知幅度正弦信号的频谱峰值单频正弦波峰值除以 N/2MCU 上 FFT 耗时过长点数过大或使用浮点运算用定时器测量 FFT 函数耗时改用定点 FFT、降低点数、用 CMSIS-DSP 优化库输入序列长度不是 2 的整数次幂数据采集长度不合适检查 buffer 大小补零到最近的 2 的整数次幂或调整采样点数输出频率只有前一半有意义实信号频谱共轭对称检查后 N/2 点只取前 N/2 点即可相位测量结果不稳定采样窗口起点变化或没有加窗固定触发采样多次测量对比使用同步采样或多帧平均8. 最佳实践与工程建议8.1 加窗策略如果被分析信号不是整周期截断直接做 FFT 会带来频谱泄漏。通用做法是选择窗函数。做幅值精度要求高的测试推荐汉宁窗做频率分辨能力要求高的分析推荐矩形窗做动态范围要求高的场景推荐布莱克曼窗。加窗的做法很简单把采样序列和窗函数逐点相乘再送入 FFT。window np.hanning(N) x_windowed x * window X fft_dit(x_windowed)注意加窗会改变信号的总能量幅值归一化公式也要做相应修正工程上常见做法是使用窗函数的相干增益Coherent Gain做补偿。8.2 补零还是加长采样补零可以提高频谱的显示密度但不能提高真实的频率分辨率。如果你需要区分两个间隔 2Hz 的频率补零到再大也没用必须延长实际采样时间来缩短 Δf。这是一个经常被误会的点。8.3 定点 FFT 的注意事项在 MCU 上如果使用定点 FFT需要特别注意 Q 格式的选择。定点运算的动态范围有限输入信号太小会降低信噪比太大会在蝶形运算的中间级产生溢出。通常的工程做法是在每一级蝶形运算后增加移位操作或者对输入信号做归一化缩放。这部分内容如果要展开会涉及arm_cfft_q15的详细用法这里先记住一个原则定点 FFT 必须控制信号幅值防止中间级溢出。8.4 数据采集设计FFT 分析的质量高度依赖采样电路。采样率要满足奈奎斯特定理最好在被分析最高频率的 3~5 倍以上同时要加抗混叠滤波器否则高频分量折叠到低频区FFT 算得再准也没有意义。很多工程问题最后排查出来不是 FFT 代码的问题而是抗混叠滤波没有做好。8.5 多帧平均在振动信号分析、声音检测这类场景中单帧 FFT 结果方差较大容易出现虚假峰值。工程上通常对多帧幅度谱做平均常见的方法有功率谱平均和幅度谱平均。平均帧数越多方差越小但实时性会下降。这块需要根据项目实时需求做平衡。9. 总结与下一步学习建议现在回头看开头那几个问题应该都有答案了FFT 点数必须是 2 的整数次幂是因为基 2 FFT 的蝶形结构要求每级都能对半分一直分到 2 点蝶形频率轴的第 0 点是直流分量第 N/2 点是奈奎斯特频率是实信号频谱共轭对称的必然结果加窗改变频谱形态是因为非整周期截断在边界引入了不连续点MCU 上做实时频谱要不要担心 CPU可以用蝶形运算总量 (N/2) * log₂N 做预估算。蝶形运算结构是 FFT 从数学公式走向工程实现的桥梁。看懂它你就不会再把 FFT 当成一个不可拆解的黑盒。下一步的实践路径可以这样安排先把文中的 Python 代码跑通改几个不同点数观察输出和numpy.fft.fft的误差自己画一遍 N8 的 DIT-FFT 信号流图把每一级的跨度和旋转因子标出来在 MCU 或嵌入式开发板上移植一个 CMSIS-DSP 的 FFT 示例用实际采样数据验证频谱尝试做一个小项目比如“基于 FFT 的简易频谱分析仪”把采样、加窗、FFT、显示串起来有余力再研究 DIF-FFT、分裂基 FFT、实序列 FFT 优化这些都是从蝶形结构延伸出来的进阶方向。如果你正面临“会调库但不懂原理”的阶段这篇文章可以作为你补上这块拼图的起点。建议收藏备用实际做项目遇到 FFT 相关的问题时再回来对照蝶形结构和排查表会比重新翻教材更高效。 FFT快速傅里叶变换蝶形运算详解：从原理到工程实现
很多做信号处理的朋友都有过这种体验numpy.fft.fft一行代码就能算出频谱STM32 里调用 DSP 库的arm_cfft_f32也能跑出结果但一旦遇到下面这些问题就很容易卡住为什么 FFT 的输出点数必须是 2 的整数次幂为什么计算出的频率轴第 0 点是直流分量第 N/2 点是奈奎斯特频率为什么加窗之后频谱会“变胖”不加窗又会“泄漏”为什么用单片机做实时频谱显示时总是担心 CPU 扛不住这些问题表面上是“FFT 使用问题”根子上却是一个东西你没有真正理解蝶形运算结构。蝶形运算是 FFT 之所以“快速”的核心也是理解频域分辨率、计算量、内存排布和实时性的钥匙。如果你只停留在调库层面遇到点数、相位、窗函数、资源占用这类问题时就只能靠试错。这篇文章围绕“数字信号处理 FFT 快速傅里叶变换蝶形运算结构方法”这条主线把 DFT 到 FFT 的演进逻辑、基 2 时间抽取DIT和频率抽取DIF两种经典结构讲透然后给出可运行的 Python 实现再结合频谱分析、相位测量、包络谱分析和 MCU 实时频谱显示这些真实场景说明蝶形结构在实际项目中到底怎么用、有哪些坑。读完你可以得到三样东西第一彻底看懂教材里那张蝴蝶交叉的信号流图第二手写一个可用的基 2 FFT而不是只会调fft函数第三在实际项目里知道自己该选多少点、加什么窗、怎么验证结果。1. 为什么 DFT 很“重”而 FFT 很“轻”先回到最原始的问题离散傅里叶变换DFT到底在算什么对一段长度为 N 的离散序列 x[n]它的 N 点 DFT 定义是X[k] sum_{n0}^{N-1} x[n] * W_N^{kn}, k 0,1,...,N-1其中旋转因子W_N^{kn} e^{-j 2π kn / N}这个公式的信息量很大。它说明DFT 的本质是把一个长度为 N 的序列分解成 N 个不同频率的复指数序列的加权和。X[0] 对应直流分量X[1] 对应基波X[k] 对应第 k 个频率分量。每一个 X[k] 的计算都需要对 N 个输入点做一次复乘和累加。所以完成全部 N 个 X[k] 的计算需要复数乘法N² 次复数加法N(N-1) 次约等于 N² 次当 N1024 时N² 1,048,576也就是约一百万次复数乘法。当 N4096 时乘法次数就变成了约 1677 万次。虽然现代 CPU 跑这个量级并不吃力但放到实时信号处理里每毫秒都要来一轮压力就完全不同了。更何况很多 MCU 根本没有硬件浮点单元一次复数乘法要拆成不知道多少条整数指令。1965 年 Cooley 和 Tukey 发表的 FFT 算法把计算量从 O(N²) 降到了 O(N log₂N)。同样是 1024 点FFT 只需要约 1024 × 10 10240 次复数运算相比直接 DFT 差了大约两个数量级。这个差距才是 FFT 能在实时系统、嵌入式设备、工业检测领域大规模落地的根本原因。2. 蝶形运算FFT 的最小组成单元要理解 FFT 为什么快必须先理解“蝶形运算”。蝶形运算是一个只有两个输入、两个输出的计算单元它的形式非常固定输入a, b 输出 X a b * W Y a - b * W其中 W 是旋转因子。之所以叫“蝶形”是因为把输入输出画成信号流图时两条支路交叉再合并形状像一只蝴蝶。蝶形运算的核心价值在于它把一个二点 DFT 的计算变成了“一次复数乘法 两次复数加法”并且在整个 FFT 的每一级里反复复用。FFT 的算法结构本质上就是把一个大 DFT 逐级拆成许多个小 DFT而每一级的小 DFT 最终都变成若干个蝶形运算。这里要特别强调一点蝶形运算里的加法和减法共享同一个乘积项 b * W。也就是说计算 a bW 和 a - bW 时复数乘法只需要做一次。这种“乘法结果复用”的思想是 FFT 能减少计算量的一个重要原因。如果你手动写过朴素 DFT会发现每算一个点都独立做 N 次乘法而 FFT 通过蝴蝶结构把中间结果反复复用避免了大重复计算。3. 基 2 时间抽取 FFTDIT-FFT的完整推导基 2 FFT 的“基 2”意思是每一级都把序列按奇偶位置分成两组直到每组只剩 2 个点。这里以最常用的按时间抽取Decimation-In-Time, DIT为例逐步拆解。3.1 第一步奇偶分离设序列长度为 N且 N 2^M。把 x[n] 分成偶数下标和奇数下标两组偶数部分x[0], x[2], x[4], ..., x[N-2] 奇数部分x[1], x[3], x[5], ..., x[N-1]代入 DFT 定义经过变形可以得到X[k] X_even[k] W_N^k * X_odd[k] X[k N/2] X_even[k] - W_N^k * X_odd[k]其中 X_even[k] 是偶数下标序列的 N/2 点 DFTX_odd[k] 是奇数下标序列的 N/2 点 DFT。这里有两个关键观察用 N/2 点 DFT 的结果可以通过加减组合出完整的 N 点 DFT。这一个“减半”就让计算量从 N² 变成了大约 2 × (N/2)² N/2是第一步优化。X[k] 和 X[kN/2] 之间恰好是一个标准的蝶形运算结构。这正是蝶形运算出现的自然来源。3.2 第二步继续对半拆分把 N/2 点的 DFT 再按奇偶拆一次变成 N/4 点的 DFT。每拆一次系统规模减半一直拆到 2 点 DFT。2 点 DFT 的公式非常简单X[0] x[0] x[1] X[1] x[0] - x[1]这正好是旋转因子为 W_2^0 1 的蝶形运算。也就是说2 点 DFT 天然就是一个蝶形。3.3 第三步码位倒序从数学推导看DIT-FFT 是“输入倒序输出顺序”。也就是说输入序列要先按照二进制位倒序重新排列然后才能逐级做蝶形运算。N8 时的输入顺序变化如下原始索引二进制位倒序重排后索引0000000010011004201001023011110641000011510110156110011371111117所以 DIT-FFT 的流程是清晰的输入序列 → 码位倒序 → 第 1 级蝶形 → 第 2 级蝶形 → ... → 第 M 级蝶形 → 输出频谱每一级蝶形的旋转因子和跨度都不一样这是初学时最容易搞混的地方。可以从 N8 的例子看规律。3.4 N8 的完整蝶形结构N8 时共有 M log₂8 3 级蝶形第 1 级蝶形跨度为 1旋转因子为 W_2^0 1第 2 级蝶形跨度为 2旋转因子为 W_4^0 1, W_4^1 -j第 3 级蝶形跨度为 4旋转因子为 W_8^0, W_8^1, W_8^2, W_8^3可以用文字描述这个结构第 3 级最关键它把前两级得到的两组 4 点结果通过 4 个蝶形组合成最终的 8 点频谱。每一级蝶形数量都是 N/24 个所以 N8 的 FFT 总共有 M × N/2 3 × 4 12 个蝶形运算。推广到一般情况N 点基 2 FFT 的蝶形总数是(N/2) * log₂N每次蝶形包含 1 次复数乘法和 2 次复数加法所以总计算量大约为复数乘法(N/2) * log₂N复数加法N * log₂N这正是 FFT 计算量 O(N log₂N) 的来源。对比 DFT 的 N²差距一目了然。4. 频率抽取 FFTDIF-FFT与 DIT-FFT 的对比除了按时间抽取的 DIT-FFT还有一种同等重要的结构按频率抽取Decimation-In-Frequency, DIF-FFT。它和 DIT-FFT 在数学上是等价的但信号流图的方向正好相反。DIF-FFT 的核心思路不是把输入按奇偶分组而是把输出按前后两半拆分。它的特点是输入是自然顺序不需要码位倒序输出是倒序的需要在使用前重新排列蝶形运算的旋转因子位置和 DIT-FFT 不同。两种结构的对比可以总结成下表对比维度DIT-FFTDIF-FFT输入顺序码位倒序自然顺序输出顺序自然顺序码位倒序旋转因子位置在蝶形输入支路在蝶形输出支路实现难度需要先做倒序需要最后做倒序使用场景通用软件实现、教材示例某些硬件流水线结构、定点实现在实际工程中DIT-FFT 更常见因为多数库函数选择“输入倒序、输出顺序”的约定。但 DSP 芯片和 FPGA 的某些 FFT IP 核可能采用 DIF 结构因为它的数据流更适合流水线。一个很关键的理解是DIT 和 DIF 不是两种不同的数学变换而是同一变换的两种计算调度方式。它们算出来的结果完全一致只是中间过程的数据流向和旋转因子的位置不同。5. 手写 Python 实现基 2 DIT-FFT 完整代码理论讲再多不如跑一段代码。下面用 Python 手写一个基 2 DIT-FFT不使用numpy.fft只借助numpy做复数数组操作。先实现码位倒序import numpy as np def bit_reverse(x): 对输入序列进行码位倒序重排。 要求 len(x) 必须是 2 的整数次幂。 n len(x) if n (n - 1) ! 0: raise ValueError(长度必须为 2 的整数次幂) 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] return x这段代码是经典的位反转算法逻辑就是用二进制位翻转的方式找到交换位置。不用记忆它直接当工具函数用即可。然后实现核心的蝶形 FFTdef fft_dit(x): 基 2 按时间抽取 FFTDIT-FFT。 x 为一维复数数组长度必须为 2 的整数次幂。 返回长度相同的复数频谱数组。 n len(x) if n (n - 1) ! 0: raise ValueError(长度必须为 2 的整数次幂) # 复制输入避免修改原数组 a np.array(x, dtypenp.complex128) a bit_reverse(a) # 开始逐级蝶形运算 size 2 while size n: half size // 2 # 当前级的旋转因子 w np.exp(-2j * np.pi * np.arange(half) / size) for start in range(0, n, size): for k in range(half): u a[start k] v a[start k half] * w[k] a[start k] u v a[start k half] u - v size 1 return a这段代码在结构上非常贴近教材里的信号流图最外层循环控制级数size 从 2 开始每轮翻倍中间层循环遍历每一组蝶形最内层循环计算同一组内的 half 个蝶形运算。运行验证一下# 生成一个包含 50Hz 和 120Hz 分量的测试信号 fs 1000 # 采样率 1000 Hz N 1024 # FFT 点数 t np.arange(N) / fs x 0.7 * np.sin(2 * np.pi * 50 * t) np.sin(2 * np.pi * 120 * t) # 计算 FFT X fft_dit(x) # 计算幅度谱 mag np.abs(X) / (N / 2) freq np.arange(N) * fs / N # 打印能量最大的 3 个频率 top_indices np.argsort(mag)[::-1][:5] for idx in top_indices[:3]: print(f频率: {freq[idx]:.2f} Hz, 幅度: {mag[idx]:.4f})预期输出类似频率: 120.00 Hz, 幅度: 1.0000 频率: 50.00 Hz, 幅度: 0.7000这说明手写 FFT 和理论完全一致。为了进一步确认还可以和numpy.fft.fft的结果对比X_np np.fft.fft(x) print(最大误差:, np.max(np.abs(X - X_np)))如果实现正确最大误差应该是 1e-12 量级来自浮点计算误差而不是一个不可忽略的大数。6. 从 Python 到嵌入式实际项目里的 FFT 怎么用理解了蝶形结构之后再回到开头提到的那些实际场景你会发现很多“工程问题”其实都是“结构问题”。6.1 参数确定采样率、点数和频率分辨率在实现一个 FFT 分析系统之前必须先明确三个参数的关系频率分辨率 Δf fs / N其中 fs 是采样率N 是 FFT 点数。如果你要区分两个频率相差 1Hz 的信号而采样率是 1000Hz那么 N 至少需要 1000。为了满足基 2 FFT工程上通常会取 1024。这就是为什么 FFT 点数往往是 1024、2048、4096 这些值——它们不是随意定的而是由“频率分辨率 → 最小 N → 向上取 2 的整数次幂”决定的。6.2 用 STM32 做实时频谱显示时的资源估算在 STM32F407 这类 MCU 上做实时频谱显示通常使用 CMSIS-DSP 库的arm_cfft_f32函数。核心调用方式大致如下#include arm_math.h #define FFT_SIZE 1024 arm_cfft_instance_f32 fft_instance; float32_t input[FFT_SIZE * 2]; // 实部虚部交错存放 float32_t output[FFT_SIZE]; void init_fft(void) { // 初始化 FFT 实例点数 1024 arm_cfft_init_f32(fft_instance, FFT_SIZE); } void compute_spectrum(void) { // input 中按 实部,虚部,实部,虚部 的顺序填充采样数据 // 实数数据虚部填 0 for (int i 0; i FFT_SIZE; i) { input[2 * i] sample_buffer[i]; input[2 * i 1] 0.0f; } // 执行 FFT第三个参数 0 表示正变换 arm_cfft_f32(fft_instance, input, 0, 1); // 计算幅度谱sqrt(real^2 imag^2) arm_cmplx_mag_f32(input, output, FFT_SIZE); // output 数组中即可得到幅度谱通常只取前 FFT_SIZE/2 个点 }这里的arm_cfft_f32定点/浮点实现内部就是按蝶形结构组织的。你在代码里只需要调用函数但理解蝶形结构能帮你回答几个实际问题为什么input要按“实部虚部交错”存放因为 SIMD 优化和内存访存的连续性要求。为什么输出只需要取前 FFT_SIZE/2 个点因为实信号的频谱关于奈奎斯特频率共轭对称后一半没有额外信息。为什么 FFT 点数越大CPU 负载越高因为蝶形总数为 (N/2) * log₂N点数翻倍蝶形数量不是翻倍而是约等于翻倍加一倍的级数运算量具体可以精确计算。从资源估算来看1024 点浮点 FFT 在 STM32F407 上通常只需要几毫秒但如果你同时还要做 LCD 刷新、ADC 采样、通信协议解析就需要仔细分配时间片。更稳妥的做法是在工程上先测量一次 FFT 的实际耗时再用定时器把采样、计算、显示切成不同的任务。6.3 包络谱分析与希尔伯特变换在预测性维护、轴承故障诊断这类场景中直接对原始振动信号做 FFT 往往看不到明显的故障特征频率。这是因为故障冲击信号是宽频带的能量分散。工程上常用的做法是先做希尔伯特变换得到解析信号再取包络最后对包络做 FFT得到包络谱。整个过程可以理解为原始振动信号 → 希尔伯特变换 → 解析信号 → 取包络 → 包络信号 → FFT → 包络谱包络谱能突出周期性冲击的特征频率配合蝶形 FFT 可以快速识别轴承故障。这也是“数字信号处理 FFT 包络谱分析”在工业预测性维护中越来越常见的原因。6.4 相位测量FFT 不仅能测幅值还能测相位。每个频点的相位可以通过复数的虚部和实部计算phase[k] arctan(imag(X[k]) / real(X[k]))在实际工程中比如用 DSP 库测量两个同频信号的相位差通常会先做 FFT找到幅度最大的频点再提取该频点的相位然后对两路信号做差。这里有一个容易踩的坑FFT 的输入起始位置会直接影响相位结果。如果采样窗口的起点不同测出来的“绝对相位”就不同但两路信号的相位差依然是稳定的。所以做相位测量时要确认两路信号使用的是同一段采样时间基准。7. 常见问题与排查方法问题现象可能原因排查方式解决方案频谱峰值位置不对采样率、FFT 点数、频率轴坐标对不上检查 fs/N 的分辨率打印频率轴确认频率轴为 arange(N) * fs / N频谱泄漏严重峰值“变胖”信号频率不是频率分辨率的整数倍用非整周期采样信号做测试加窗函数如汉宁窗、布莱克曼窗幅值偏小或偏大归一化方式不对对比已知幅度正弦信号的频谱峰值单频正弦波峰值除以 N/2MCU 上 FFT 耗时过长点数过大或使用浮点运算用定时器测量 FFT 函数耗时改用定点 FFT、降低点数、用 CMSIS-DSP 优化库输入序列长度不是 2 的整数次幂数据采集长度不合适检查 buffer 大小补零到最近的 2 的整数次幂或调整采样点数输出频率只有前一半有意义实信号频谱共轭对称检查后 N/2 点只取前 N/2 点即可相位测量结果不稳定采样窗口起点变化或没有加窗固定触发采样多次测量对比使用同步采样或多帧平均8. 最佳实践与工程建议8.1 加窗策略如果被分析信号不是整周期截断直接做 FFT 会带来频谱泄漏。通用做法是选择窗函数。做幅值精度要求高的测试推荐汉宁窗做频率分辨能力要求高的分析推荐矩形窗做动态范围要求高的场景推荐布莱克曼窗。加窗的做法很简单把采样序列和窗函数逐点相乘再送入 FFT。window np.hanning(N) x_windowed x * window X fft_dit(x_windowed)注意加窗会改变信号的总能量幅值归一化公式也要做相应修正工程上常见做法是使用窗函数的相干增益Coherent Gain做补偿。8.2 补零还是加长采样补零可以提高频谱的显示密度但不能提高真实的频率分辨率。如果你需要区分两个间隔 2Hz 的频率补零到再大也没用必须延长实际采样时间来缩短 Δf。这是一个经常被误会的点。8.3 定点 FFT 的注意事项在 MCU 上如果使用定点 FFT需要特别注意 Q 格式的选择。定点运算的动态范围有限输入信号太小会降低信噪比太大会在蝶形运算的中间级产生溢出。通常的工程做法是在每一级蝶形运算后增加移位操作或者对输入信号做归一化缩放。这部分内容如果要展开会涉及arm_cfft_q15的详细用法这里先记住一个原则定点 FFT 必须控制信号幅值防止中间级溢出。8.4 数据采集设计FFT 分析的质量高度依赖采样电路。采样率要满足奈奎斯特定理最好在被分析最高频率的 3~5 倍以上同时要加抗混叠滤波器否则高频分量折叠到低频区FFT 算得再准也没有意义。很多工程问题最后排查出来不是 FFT 代码的问题而是抗混叠滤波没有做好。8.5 多帧平均在振动信号分析、声音检测这类场景中单帧 FFT 结果方差较大容易出现虚假峰值。工程上通常对多帧幅度谱做平均常见的方法有功率谱平均和幅度谱平均。平均帧数越多方差越小但实时性会下降。这块需要根据项目实时需求做平衡。9. 总结与下一步学习建议现在回头看开头那几个问题应该都有答案了FFT 点数必须是 2 的整数次幂是因为基 2 FFT 的蝶形结构要求每级都能对半分一直分到 2 点蝶形频率轴的第 0 点是直流分量第 N/2 点是奈奎斯特频率是实信号频谱共轭对称的必然结果加窗改变频谱形态是因为非整周期截断在边界引入了不连续点MCU 上做实时频谱要不要担心 CPU可以用蝶形运算总量 (N/2) * log₂N 做预估算。蝶形运算结构是 FFT 从数学公式走向工程实现的桥梁。看懂它你就不会再把 FFT 当成一个不可拆解的黑盒。下一步的实践路径可以这样安排先把文中的 Python 代码跑通改几个不同点数观察输出和numpy.fft.fft的误差自己画一遍 N8 的 DIT-FFT 信号流图把每一级的跨度和旋转因子标出来在 MCU 或嵌入式开发板上移植一个 CMSIS-DSP 的 FFT 示例用实际采样数据验证频谱尝试做一个小项目比如“基于 FFT 的简易频谱分析仪”把采样、加窗、FFT、显示串起来有余力再研究 DIF-FFT、分裂基 FFT、实序列 FFT 优化这些都是从蝶形结构延伸出来的进阶方向。如果你正面临“会调库但不懂原理”的阶段这篇文章可以作为你补上这块拼图的起点。建议收藏备用实际做项目遇到 FFT 相关的问题时再回来对照蝶形结构和排查表会比重新翻教材更高效。