简介:这份资源面向嵌入式开发初学者与数字信号处理方向的工程师,围绕STM32F1系列单片机(ARM Cortex-M3内核)讲解并实现快速傅里叶变换,帮助读者把时域采样数据转换到频域,用于通信、音频与图像等信号分析场景。压缩包共225个文件,以60个h头文件与56个c源文件为核心,配合30个o目标文件、30个d依赖文件及27个crf等编译中间产物,另有uvprojx、uvoptx工程文件、sct分散加载脚本、hex烧录文件与map映射文件,整体约4.12MB,工程结构完整可直接编译。内容涉及Cortex-M3内核特性、Cooley-Tukey蝶形运算与位反转、SRAM内存布局、TIM定时器中断触发、CMSIS-DSP库函数调用及误差与实时性分析等知识点。目前已有10174人学习下载,适合希望掌握FFT在单片机上落地方法、对照工程源码理解算法实现与性能优化的读者参考。
1. STM32 上跑 FFT:为什么 1024 点只要几百微秒,而你的代码却卡在中断里
很多人第一次在 STM32 上做 FFT,都会经历同一个反直觉现象:明明芯片主频标称 72MHz 甚至 168MHz,算 1024 点复数 FFT 理论上不过几十万次乘加,可实际跑起来要么慢得离谱,要么结果全是噪声。问题往往不在 FFT 算法本身,而在于没搞懂单片机做频域变换时,数据是怎么从 ADC 搬进内存、定点数怎么防溢出、蝶形运算的旋转因子又从哪来。这套资源把 STM32 实现 FFT 的完整链路拆开了——从采样、加窗、定点转浮点,到基 2 蝶形、位反转、复数乘法的底层写法,再到 CMSIS-DSP 库和手写汇编的取舍。搞懂它,你就能明白为什么同样一颗芯片,有人能做到微秒级响应,有人却连实时频谱都刷不出来。
2. FFT 在单片机上到底难在哪:从 DFT 的 O(N²) 到基 2 蝶形的 O(N log N)
2.1 为什么不能直接算 DFT
离散傅里叶变换的定义很直白:对 N 个采样点,每个输出频点都要和所有输入点做一次复数乘加。N=1024 时,复数乘法次数是 N² ≈ 104 万次,每次复数乘法又要拆成 4 次实数乘法和 2 次加法。在 STM32F103 这类没有硬件浮点单元的 Cortex-M3 上,一次单精度浮点乘法靠软件模拟要几十个周期,104 万次乘加直接就是几千万周期,实时性根本无从谈起。
FFT 的核心贡献是把这 N² 次运算通过对称性和周期性降到 N/2 · log₂N 次蝶形。1024 点时,蝶形数量是 5120 个,每个蝶形包含 1 次复数乘法(4 实数乘 + 2 实数加)和 2 次复数加法。总运算量降到约 5 万次实数乘法量级,比 DFT 少了近 20 倍。这个数量级差异,才是单片机能不能做实时频谱的分水岭。
2.2 基 2 蝶形的数据流:位反转才是第一个坑
基 2 按时间抽取 FFT 的第一步不是乘法,而是位反转重排。因为蝶形运算要求输入序列按二进制位倒序排列,比如 8 点 FFT 的输入顺序是 0,4,2,6,1,5,3,7。很多手写实现结果全错,就是漏了这一步,或者位反转写成了字节反转。
常见的实现思路是用查表法:预先算好 N 个索引的位反转值存成数组,运行时直接按表搬数据。1024 点就是 1024 个 uint16_t,占 2KB Flash,换掉每次运行时逐位反转的循环开销。对于 RAM 紧张的型号,也可以原地交换:只遍历前半段,把 i 和 bitrev(i) 对应的元素互换,避免额外缓冲区。
// 位反转重排:原地交换,只处理 i < rev_i 的情况 void bit_reverse(float *real, float *imag, uint16_t n) { uint16_t j = 0; for (uint16_t i = 0; i < n - 1; i++) { if (i < j) { // 交换实部和虚部 float tr = real[i]; real[i] = real[j]; real[j] = tr; float ti = imag[i]; imag[i] = imag[j]; imag[j] = ti; } // 计算下一个位反转索引:从高位向低位进位 uint16_t k = n >> 1; while (k <= j) { j -= k; k >>= 1; } j += k; } }这段代码的关键在while (k <= j)循环:它模拟了二进制加法从最高位开始进位的过程。j始终是i的位反转值,但不需要每次重新算,而是递推得到。参数n必须是 2 的幂,否则k的移位逻辑会出错。实际项目中如果采样点数固定,更推荐查表,因为递推里的分支在 Cortex-M 上会有流水线停顿。
2.3 旋转因子的存储与精度取舍
蝶形运算里的旋转因子 W_N^k = e^(-j2πk/N),直接每次用sinf/cosf算会慢到无法接受。标准做法是预计算一张 N/2 点的表,存成 float 或 q15/q31 定点格式。1024 点浮点表占 4KB(实虚各 2KB),定点 q15 表占 2KB。
这里有个精度陷阱:如果表用 q15 定点存,但蝶形运算用 float 算,每次查表都要做 q15→float 转换,反而比直接存 float 慢。反过来,如果全用 q15 定点运算,中间结果很容易溢出——两个 q15 相乘得到 q30,累加几次就超出 32 位范围。常见的做法是蝶形内部用 q31 或 float,只在输入输出做定点转换。STM32F4 系列有 FPU,直接全 float 最省心;F0/F1 系列没有 FPU,就得在精度和速度之间做选择。
3. 从 ADC 到频谱:采样、加窗、定点化每一步都在决定成败
3.1 采样率与点数怎么定:频率分辨率不是越高越好
FFT 的频率分辨率等于采样率除以点数。比如采样率 10kHz、1024 点,分辨率约 9.77Hz。想提高分辨率,要么降采样率,要么增点数。但降采样率会压缩分析带宽,增点数会吃掉 RAM 和算力。STM32F103 只有 20KB RAM,1024 点复数 float 数组就要 8KB,再加旋转因子表 4KB,留给其他任务的空间已经很紧张。
更实际的做法是根据目标信号选参数。比如分析 50Hz 工频谐波,采样率取 3.2kHz、512 点,分辨率 6.25Hz,足够分辨到 13 次谐波,RAM 占用也降到 4KB。如果非要 4096 点,就得换 F4 系列或者用外部 SRAM。这套资源里给了不同点数下的 RAM/Flash 占用对照,选型时可以直接查。
3.2 加窗:为什么矩形窗会漏掉你的信号
直接对采样数据做 FFT,相当于加了矩形窗。矩形窗的旁瓣衰减只有 -13dB,强信号旁边的弱谐波会被旁瓣淹没。实测中如果发现频谱上某个频点附近有一片模糊的抬升,大概率就是没加窗导致的频谱泄漏。
汉宁窗是单片机里最常用的折中:旁瓣衰减 -31dB,主瓣宽度是矩形窗的两倍。代价是每个采样点要乘一个窗系数,1024 点就是 1024 次乘法。如果采样和 FFT 在同一颗芯片上做,可以在 ADC 中断里顺手乘掉,不额外增加 FFT 阶段的开销。窗系数表同样预存在 Flash 里,用 q15 格式存 1024 个点只占 2KB。
// 汉宁窗系数预计算(运行时初始化一次) // w[n] = 0.5 * (1 - cos(2*pi*n/(N-1))) void init_hanning_window(float *window, uint16_t n) { for (uint16_t i = 0; i < n; i++) { window[i] = 0.5f * (1.0f - cosf(2.0f * 3.1415926f * i / (n - 1))); } } // 采样时应用:adc_buf[i] *= window[i]注意n-1而不是n,这是汉宁窗的标准定义。如果用n,窗的两端不闭合,会引入额外的直流分量。实际项目中如果对实时性要求极高,可以改用海明窗,系数用定点近似,省掉cosf的初始化时间。
3.3 定点 vs 浮点:没有 FPU 的芯片怎么选
STM32F0/F1 没有硬件浮点单元,float 运算靠软件库,一次乘法可能 20~50 周期。这种情况下用 q15 定点做 FFT 是主流选择。q15 表示范围 [-1, 1),精度约 3e-5,对于 12 位 ADC(精度 2.4e-4)来说足够。
但定点蝶形有个经典问题:每级运算后幅度可能增长,需要右移防溢出。基 2 蝶形每级增益最大为 2,log₂N 级后总增益可达 N。1024 点就是 1024 倍,q15 根本存不下。常见做法是每级右移 1 位,总共右移 log₂N 位,最后结果缩小了 N 倍,但不会溢出。代价是信噪比下降,小信号可能被移没了。另一种做法是块浮点:每级检测最大值,动态决定移位数,复杂度高但精度好。这套资源里两种方案都有实现,可以对比着看。
4. 蝶形运算的源码级拆解:三层循环怎么写才不拖慢流水线
4.1 蝶形核心:一次复数乘法到底做了几次实数运算
基 2 蝶形的公式是:
- X[k] = A + W·B
- X[k+N/2] = A - W·B
其中 A、B 是复数,W 是旋转因子。展开后,W·B 需要 4 次实数乘法和 2 次实数加法,再加上两次复数加减,总共 4 乘 6 加。手写代码时如果把复数结构体传来传去,编译器可能生成大量冗余的 load/store。高效写法是直接用局部 float 变量展开。
// 单个蝶形:a 和 b 是输入,w 是旋转因子,输出写回 a 和 b // 所有参数都是指针,避免结构体拷贝 void butterfly(float *ar, float *ai, float *br, float *bi, float wr, float wi) { // W * B 的实部和虚部 float tr = wr * (*br) - wi * (*bi); float ti = wr * (*bi) + wi * (*br); // A + W*B 和 A - W*B float out_r = *ar + tr; float out_i = *ai + ti; *br = *ar - tr; // 注意:先算完再覆盖,顺序不能反 *bi = *ai - ti; *ar = out_r; *ai = out_i; }这里*br和*bi的赋值必须在*ar、*ai更新之前完成,否则A - W·B会用错值。这是手写蝶形最常见的翻车点,编译器不会报错,但结果全错。参数wr、wi直接传值而不是指针,是因为旋转因子在循环里是连续读取的,传值能让编译器更好地做寄存器分配。
4.2 三层循环的组织:为什么外层的步长要单独提出来
标准 FFT 有三层循环:最外层是级数(log₂N 级),中间层是每级的蝶形组数,最内层是组内蝶形。很多实现把旋转因子索引放在最内层算,每次都要做一次乘法和取模,白白浪费周期。
优化做法是把旋转因子索引的步长提到中间层。每级的步长是 N/2^s(s 是级数),组内索引递增这个步长即可。这样最内层只剩加法和数组访问,编译器更容易展开和流水线化。
// 迭代法基 2 FFT 主循环 void fft(float *real, float *imag, uint16_t n, float *wr_tab, float *wi_tab) { bit_reverse(real, imag, n); // 先做位反转 for (uint16_t s = 1; s <= log2(n); s++) { uint16_t m = 1 << s; // 当前级的蝶形跨度 uint16_t half = m >> 1; uint16_t step = n / m; // 旋转因子步长 for (uint16_t g = 0; g < n; g += m) { // 每组 for (uint16_t k = 0; k < half; k++) { // 组内 uint16_t idx = k * step; butterfly(&real[g+k], &imag[g+k], &real[g+k+half], &imag[g+k+half], wr_tab[idx], wi_tab[idx]); } } } }step = n / m是旋转因子表的索引步长,idx = k * step直接查表。如果表是按 N 点全存,这里索引不会越界;如果只存 N/2 点,需要根据idx是否超过 N/2 做符号翻转。这套资源里的表是 N/2 点存储,代码里用了一个if (idx >= n/2)分支处理,实测比全表慢约 5%,但省了一半 Flash。
4.3 用 CMSIS-DSP 还是手写:库函数的隐藏开销
CMSIS-DSP 提供了arm_cfft_f32、arm_rfft_q15等现成函数,调用简单,但有几个隐藏成本。第一,库函数为了通用性,支持任意点数(2 的幂),内部有额外的参数检查;第二,q15 版本的缩放策略是固定的,不能按信号动态调整;第三,库函数占用的 Flash 比手写大,F103 上arm_cfft_f32约 6KB,手写优化版可以压到 2KB 以内。
如果点数固定、对体积极敏感,手写更合适。如果项目里还要用 FIR、IIR 等其他 DSP 功能,CMSIS-DSP 的整体占用反而更划算。这套资源里两种实现都有,并且给了同一块板子上的实测周期数对比,可以直接看数据做决定。
5. 避坑与排查:FFT 结果不对时先查这五个地方
5.1 频谱全是直流或低频抬升
现象:FFT 输出在 0 频附近有一个巨大峰值,其他频点几乎看不见。原因通常是采样数据没有去掉直流偏置。12 位 ADC 的输出范围是 0~4095,如果直接当 float 用,相当于叠加了一个 2048 的直流分量,FFT 后 0 频点自然巨大。解决:采样后先减去均值,或者硬件上用运放做交流耦合。软件去均值可以在采样缓冲填满后做一次累加求平均,1024 点只多 1024 次加法,开销可忽略。
5.2 频谱镜像对称,但左右幅值不一致
现象:实数输入的 FFT 结果应该关于 N/2 对称,但实测左右幅值差很多。原因多半是位反转没做对,或者旋转因子表的索引算错了。排查方法:用一个已知频率的正弦波输入,看峰值是否出现在正确的 bin 上。如果峰值位置对但幅值不对,查窗函数是否应用正确;如果峰值位置都不对,查位反转和蝶形循环的边界条件。常见错误是half算成了m而不是m>>1,导致蝶形跨度过大。
5.3 定点运算结果溢出成乱码
现象:q15 定点 FFT 输出在某些频点突然跳到最大值或最小值,波形完全失真。原因是蝶形中间结果溢出。q15 乘 q15 得到 q30,两个 q30 相加可能超过 32 位有符号范围。解决:要么每级右移 1 位,要么用 q31 中间变量并在乘法后立即饱和处理。CMSIS-DSP 的 q15 FFT 内部有饱和指令,手写时如果忘了加__SSAT,溢出就会回绕。排查时可以在每级蝶形后打印最大值,看是否接近 32767。
5.4 中断里做 FFT 导致采样丢点
现象:频谱上出现不该有的杂散,或者低频段噪声异常。原因是在 ADC 中断里直接调用了 FFT 函数,FFT 执行时间超过了采样周期,导致下一次 ADC 转换完成时数据没被读走。解决:ADC 中断只负责把数据存进缓冲区,FFT 在主循环里做。如果必须实时,用 DMA 搬运 ADC 数据,中断只做缓冲区切换。STM32 的 DMA 循环模式配合半满/全满中断,可以做到零 CPU 开销采样。
5.5 旋转因子表精度不够导致底噪抬高
现象:频谱底噪比理论值高 10dB 以上,小信号被淹没。原因是旋转因子用 q15 存储,量化误差在每级蝶形中累积。1024 点有 10 级,误差累积后可能达到 -60dB 左右。解决:旋转因子表用 q31 或 float 存储,只在蝶形输入输出做定点转换。如果 Flash 够,直接存 float 表最省事。这套资源里给了 q15 表和 float 表的底噪实测对比,差距在 15dB 左右。
6. 定点与浮点的取舍:没有最优解,只有匹配场景的选择
在 STM32 上做 FFT,最终都会落到一个问题上:用定点还是浮点。这不是技术优劣的问题,而是场景匹配的问题。F4 系列有 FPU,全 float 实现 1024 点约 200~300 微秒,代码可读性好,调试方便,适合快速原型和小批量产品。F1 系列没有 FPU,float 版本要 2~3 毫秒,q15 定点可以压到 500 微秒以内,但溢出处理和精度控制需要额外精力。
我一般会这么理解:如果产品对成本敏感、出货量大,花时间优化定点实现是值得的;如果只是做验证或者出货量小,直接用带 FPU 的型号省下来的开发时间更值钱。这套资源里两种实现都有完整代码和实测数据,可以拿同一段信号分别跑,看周期数和频谱底噪的差异,再决定自己的项目走哪条路。
还有一个容易被忽略的点:FFT 只是频谱分析的一半。算出频域数据后,如果还要做峰值检测、谐波分析、总谐波失真计算,这些后处理的开销可能和 FFT 本身相当。定点实现里,后处理同样要用定点,否则来回转换的开销会吃掉 FFT 省下来的时间。资源里给了一个完整的谐波分析示例,从 ADC 采样到 THD 输出全链路都是 q15,可以参考它的数据流设计。
从那以后我每次在单片机上做频域相关的功能,都会先把采样率、点数、窗函数、数据格式这四个参数列成一张表,算清楚 RAM 占用和周期预算再动手。很多翻车不是因为算法难,而是参数没对齐。希望帮到你。
本文还有配套的精品资源,点击获取