企业官网建设流程全解析

这次我们来看数字信号处理里出场率最高的一个算法结构：FFT 的蝶形运算。

网上的 FFT 原理文章很多，但多数停留在“FFT 比 DFT 快很多”这种结论上。真正往下走就会碰到几个问题：蝶形图怎么看、旋转因子怎么算、什么是倒位序、为什么能原位计算、STM32 上调用 DSP 库做实时频谱要注意什么。这篇就把 FFT 的结构一层层拆开，重点讲蝶形运算的两种典型结构，并给出可运行的代码示例和排查思路。

文章适合这几类读者：正在学数字信号处理的大学生；需要在 MCU、实时系统里实现频谱分析的嵌入式开发者；用 MATLAB/Python 做振动分析、包络谱、故障诊断，但总搞不清 FFT 参数怎么选的人。看完之后，你应该能独立画出 8 点 FFT 蝶形图，理解旋转因子的查表逻辑，并搞清楚为什么很多 FFT 例程都要先做位序反转。

1. FFT 蝶形运算核心能力速览

先给一张速览表，把整个 FFT 蝶形结构方案的关键点列出来。

这里先说一个关键判断：无论网上代码怎么变，FFT 的蝶形运算结构都逃不出“分成上下两路，一路直接加，一路乘旋转因子后相加减”这个基本单元。把这一条刻在脑子里，后面所有内容都好理解。

2. FFT 要解决什么问题：从 DFT 复杂度说起

离散傅里叶变换的定义是：

X[k] = Σ x[n] * e^(-j2πkn/N)，n=0 到 N-1

如果直接按这个公式算，每个频点 k 都要做 N 次复乘和 N 次复加，N 个频点就是 N^2 次复数乘法。N 取 1024 时，计算量超过一百万次复乘；N 取 4096 时，直接算就非常吃力。对嵌入式场景来说，这种复杂度很难支撑实时频谱显示。

FFT 的核心思路是拆解。观察旋转因子 W_N^k = e^(-j2πk/N)，可以发现两个性质：

• 周期性：W_N^(k+N) = W_N^k
• 对称性：W_N^(k+N/2) = -W_N^k

这两个性质让很多重复计算可以合并。蝶形运算就是把一个大的 DFT 逐级拆成一半一半的小 DFT，每一级用“加/减一个旋转因子乘积”来合并结果。最终的计算量约为 (N/2) * log2(N) 次复数乘法，相比直接 DFT 是数量级的下降。

所以 FFT 并不是一种全新的变换，它只是 DFT 的高效计算方式。蝶形运算结构就是“高效”二字的落地形态。

3. 蝶形运算的基本单元

蝶形运算的名字来自数据流图的外形：两条输入线从左边进来，交叉合并后从右边出去，看起来像蝴蝶翅膀。

基本的 DIT 蝶形单元可以表示成：

用公式写就是：

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 点时：

这就是 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 控制，内层循环按同一级不同蝶形组推进。

一个容易出错的点：如果旋转因子索引算错，原位覆盖后数据会被污染，而且很难追查。排查时建议先用 N=4 或 N=8 的小点数，把每级数据打印出来和 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) // 当前蝶形在组内的位置为 j，k 为旋转因子索引 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, dtype=complex) 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 Hz，FFT 长度 1024，分辨率约 7.8125 Hz。这个数值决定了你能区分多近的两个频率峰。

如果测振动信号，采样率通常远高于轴承特征频率。预测性维护场景中常用包络谱：先让信号通过带通滤波器和希尔伯特变换求出包络波形，再对包络做 FFT 得到包络谱，用于提取轴承故障特征频率。FFT 蝶形结构本身不关心输入是原始加速度信号还是包络信号，但采样参数和分析流程差异会影响最终频率成分解释。

10. FFT 常见问题与排查方法

下面是工程中最高频的几类问题，整理成排查表：

最容易踩的坑是频谱泄漏。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 点。先用 N=8 或 N=16 的固定序列，和 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 不是终点，但蝶形运算结构是所有频域分析工具的地基，值得花时间彻底吃透。