STM32F407上利用CMSIS-DSP实现FFT/IFFT信号还原与频域滤波

开头做电机振动监测的时候我拿着一台F407开发板和电流探头去现场采波形回来看波形心里就开始嘀咕这到底是哪个频率在抖轴承故障特征频率一般在几十赫兹到几百赫兹基频、谐波、边带混在一起光看时域波形根本分不清。以前我总习惯把数据用串口发给PC用Python或者MATLAB算FFT但现场设备要求闭环实时处理数据根本来不及往上传。这时候我才真正意识到STM32F407配上CMSIS-DSP库完全有能力在自己身上完成FFT和IFFT信号还原。我用的是CMSIS-DSP库里的arm_rfft_fast_f32函数。它厉害的地方不只是能做FFT把最后一个参数改成1它就变成了IFFT可以把频域数据还原回时域信号。这就产生了一个非常实用的玩法在单片机里做频域“手术”——采集时域信号、FFT到频域、把不想要的频率分量清掉、再IFFT还原成干净的时域波形。整条链路在F407上跑1024点FFT只需要零点几毫秒实时性完全够用。这篇文章我把完整代码、函数内部逻辑、输出数据布局以及我实际调试中踩过的坑全部讲透适合手里有F407开发板、想入门嵌入式信号处理的朋友参考。1. 为什么要在MCU上做FFT/IFFT而不是丢给上位机1.1 现场场景不允许“先上传再分析”很多人一听到FFT就觉得这是MATLAB和PC干的事但实际工程里有一大堆场景必须把频谱分析放到MCU本地。比如我做过的电机在线状态监测传感器就装在设备旁边控制器也部署在现场整台设备没有上位机只有一个串口用来输出最终的报警结果。如果每次都要把原始波形传回PC再分析不仅布线麻烦实时性也完全跟不上。设备故障的特征往往藏在频谱细节里比如轴承外圈故障的特征频率会在特定频段出现峰值齿轮磨损会带来边带成分。这些特征用FFT一看就很明显但用肉眼盯示波器波形很难发现。所以MCU本地做FFT不是炫技是真的有工程需求。另外像便携式音频降噪、电网谐波分析、振动烈度评估这类设备本质上都需要在嵌入式系统内部完成从时域到频域的处理闭环。FFT计算完马上就能判断“有没有超标”“故障特征频率在不在”然后直接输出报警或控制信号。这个优势是上位机方案给不了的。1.2 F407为什么适合干这个活STM32F407的核心配置非常适合做实数FFT主频168MHz带单精度硬件浮点单元FPU还有ARM的SIMD和DSP指令集。FFT这种大量乘加运算、大量浮点操作的任务正好打在F407的强项上。CMSIS-DSP库是ARM官方出品的数字信号处理函数库FFT只是其中一部分。它提供了多种FFT接口浮点的arm_rfft_fast_f32、arm_cfft_f32定点的arm_rfft_q15、arm_rfft_q31。用官方库最大的好处是不用自己写基2蝶形运算而且ARM针对Cortex-M系列做过多轮优化内部已经用到了FPU和指令级加速。我实测过典型的性能数据不同编译优化等级下会有出入在168MHz主频、开启硬件FPU、-O2优化的情况下256点实数FFT大概只需要几十微秒1024点大约在0.1ms到0.3ms范围内2048点也不到1ms。配上一个4.8kHz采样率、256点FFT的场景采满一帧数据要50多毫秒而FFT计算只占其中不到1%剩下的时间完全够做其他逻辑。1.3 FFT和IFFT的本质拆解与还原用生活化的方式理解时域波形就像做菜时锅里各种调料随时间变化的总味道而频域则是把这道菜拆解成“多少克盐、多少克糖、多少克醋”的配方。FFT就是把混合味道拆成配方IFFT则是按配方重新把味道还原出来。信号处理里经常需要“按配方改菜”比如原始信号里混了50Hz工频干扰和120Hz振动信号我们做FFT之后看到频谱上两个峰把120Hz对应的分量去掉再做IFFT得到的就是去掉干扰后的时域波形。这就是数字滤波器而且是频域滤波——和FIR、IIR滤波器相比思路更直观频率选择性也更容易控制。所以标题里的“信号还原”并不是简单地把FFT结果变回去而是利用“时域→频域→时域”的闭环在频域里做修改再还原出需要的时域信号。F407上跑这段逻辑关键就是把arm_rfft_fast_f32这个函数用熟。2. arm_rfft_fast_f32的底层逻辑与使用边界2.1 实数FFT为什么能省一半资源Alice通过ADC采集到的信号几乎都是实数序列比如一个数组里存了1024个float每个float都是某个时刻的采样值。对实数序列做FFT频谱有一个天然性质共轭对称。也就是说正频率部分包含了全部信息负频率部分是正频率部分的镜像没有额外信息。arm_rfft_fast_f32就是利用这个特性内部用N/2点复数FFT来实现N点实数FFT。这样计算量直接少了一半内存占用也跟着省了一半。这也是为什么它叫“fast”而且专门有个rfft前缀——它和arm_cfft_f32这种通用复数FFT是两套逻辑。具体点说如果你用arm_cfft_f32做1024点实数FFT你得先把1024个实数扩展成1024个复数虚部置0这就相当于2048个float的输入而且蝶形运算阶数也是完整的1024点。而arm_rfft_fast_f32只需要输入1024个float内部用512点复数FFT加一些拆分重组就完成了。计算量差异非常明显在实时系统里这就是能不能跑得动的区别。2.2 输出数组的内存布局一定得搞明白新手用arm_rfft_fast_f32最容易懵的地方就是输出结果的排列方式。它和文档里常见的复数FFT输出不太一样并不是“第一个bin实部、第一个bin虚部、第二个bin实部、第二个bin虚部……”这种自然交错排列而是采用了紧凑打包格式。对于N点FFT输出数组本身只需要N个float就够了排列顺序如下数组下标含义0DC分量频率0Hz实数1Nyquist频率分量采样率/2实数2、3第1个频率bin的实部、虚部4、5第2个频率bin的实部、虚部…………N-2、N-1第N/2-1个频率bin的实部、虚部这里有几个容易误解的地方。第一DC和Nyquist分量都只有实部没有虚部它们被单独放在数组最前面。第二我们实际能看的有效频率范围是0到采样率/2对应0到N/2这些bin其中bin 1到bin N/2-1是复数由实部和虚部两个float表示。第三不要傻乎乎地把整个输入数组当成普通复数数组去处理。所以我写程序时会定义一个辅助函数来读第k个bin的实部和虚部void get_bin_complex(float *fftOut, int fftSize, int k, float *re, float *im) { if (k 0) { *re fftOut[0]; *im 0.0f; } else if (k fftSize / 2) { *re fftOut[1]; *im 0.0f; } else { *re fftOut[2 * k]; *im fftOut[2 * k 1]; } }反过来设置某个bin的值时也要按这个布局写入。搞错布局是很多“FFT结果看起来不对”问题的根源。2.3 初始化函数与点数约束arm_rfft_fast_f32使用前必须先调用arm_rfft_fast_init_f32初始化实例结构体。这个初始化不是随便填个数就行它内部会计算旋转因子表twiddle factors这些表的值是后面FFT运算的基础。arm_rfft_fast_instance_f32 S; arm_rfft_fast_init_f32(S, FFT_SIZE);FFT_SIZE必须是2的幂。底层实现依赖N/2点复数FFT而arm_cfft_f32要求阶数至少是4之类的约束所以arm_rfft_fast_f32实际可用的点数也有下限通常是32点起步。上限主要受内存影响T407有192KB RAM2048点是很轻松的4096点也能跑但内存占用、耗时和实时性需要综合考虑。初始化之后实例结构体S里会存有内部FFT实例、旋转因子表指针、FFT点数等信息。这个结构体本身也占内存点数越大表越长我一般把它定义成全局变量避免在函数内反复初始化。2.4 in-place操作支持arm_rfft_fast_f32支持原地运算也就是说输入缓冲区和输出缓冲区可以是同一个指针。这样能省一块内存代码写起来也简洁arm_rfft_fast_f32(S, sampleBuf, sampleBuf, 0);但如果后续还需要原始时域数据做对比就别用in-place老老实实准备两个数组。我在后面的滤波例程里用两个数组主要是为了最后能对比原始信号和还原信号。3. 完整代码实现从时域采样到FFT/IFFT信号还原3.1 工程配置里的三件大事先说工程配置这三件事不做对后面跑出来的结果全是错的。第一必须启用FPU。IAR和Keil里需要在工程选项的Target页面勾选“Use Single Precision”或者选择带FPU的芯片型号GCC编译环境下编译选项要加上-mfpufpv4-sp-d16 -mfloat-abihard。如果忘了开FPUarm_rfft_fast_f32照样可以调用但浮点运算会走软件模拟耗时暴增几十倍而且一旦代码里有浮点类型不匹配还可能触发HardFault。第二把CMSIS-DSP库添加进工程。用STM32CubeMX建工程的话勾选CMSIS的DSP模块并启用库文件会自动添加。如果用Keil Pack方式或者直接从GitHub拉CMSIS源码需要把Source/DSP目录下的TransformFunctions和CommonTables等编译进来确保arm_rfft_fast_f32和对应的表格符号能找到。第三开启编译优化。Debug默认的-O0可以正常跑功能但性能会差很多。做FFT这种计算密集型任务建议至少用-O2实际我都是直接-O3两。优化等级对FFT耗时影响很大我实测-O2比-O0能快三倍以上。3.2 造一个测试信号为了验证FFT/IFFT链路我先生成一组合成信号模拟ADC采集到的数据。我创造的是50Hz和120Hz两个正弦波的叠加采样率设为1024Hz采样点数也取1024点。采样率除以点数得到频率分辨率正好是1Hz50Hz对应第50个bin120Hz对应第120个bin整周期采样频谱不会泄漏。#include arm_math.h #define FFT_SIZE 1024 #define SAMPLE_RATE 1024.0f float32_t input[FFT_SIZE]; float32_t fftOut[FFT_SIZE]; float32_t recon[FFT_SIZE]; void generate_test_signal(void) { for (int i 0; i FFT_SIZE; i) { float32_t t (float32_t)i / SAMPLE_RATE; input[i] 1.5f * sinf(2.0f * PI * 50.0f * t) 0.8f * sinf(2.0f * PI * 120.0f * t); } }幅度我故意设成1.5和0.8不取整数方便后面验证幅值换算是否准确。如果你有信号发生器或真实传感器也可以用ADC直接采这里先用合成信号把开发流程跑通。3.3 正向FFT把时域数据变到频域正向FFT的调用非常简单arm_rfft_fast_instance_f32 fftInst; void setup_fft(void) { arm_rfft_fast_init_f32(fftInst, FFT_SIZE); } void do_forward_fft(void) { arm_rfft_fast_f32(fftInst, input, fftOut, 0); }这里最后一个参数0表示正向变换FFT。调用完之后fftOut里存放的就是按第2章那个紧凑布局排列的频谱。我通常会先把幅值谱打出来确认一下float32_t magnitude[FFT_SIZE / 2 1]; void calc_magnitude_spectrum(void) { for (int k 0; k FFT_SIZE / 2; k) { float32_t re, im; get_bin_complex(fftOut, FFT_SIZE, k, re, im); magnitude[k] sqrtf(re * re im * im); printf(bin %3d freq %4.1f Hz mag %.3f\n, k, k * SAMPLE_RATE / FFT_SIZE, magnitude[k]); } }从打印结果能看到第50和120个bin附近幅值明显突出。注意这里直接算出来的幅值是FFT运算后的值不是信号的真实物理幅度。非DC和非Nyquist bin的物理幅值大约等于magnitude乘2再除以NDC和Nyquist则是乘1再除以N。所以50Hz分量的真实幅值约等于magnitude[50] * 2 / 1024算出来应该接近1.5。3.4 频域“手术”修改频谱后再IFFT还原现在来玩关键的一步——在频域里把120Hz分量清掉再做IFFT还原等效实现一个只保留50Hz的窄带滤波器效果。120Hz对应第120个bin我把这个bin的实部和虚部都清零为了消除泄漏带来的影响实际操作中我会把附近几根bin也一起清零这里先用单bin演示void apply_frequency_domain_filter(void) { // 清除第120个bin对应120Hz // 按紧凑布局bin不是0也不是Nyquist实部在 ffOut[2*k]虚部在 ffOut[2*k1] fftOut[2 * 120] 0.0f; fftOut[2 * 120 1] 0.0f; // 如果需要清除一个频段可以在这里加循环 }修改完成后调用同一个函数把最后一个参数改成1就变成IFFTvoid do_inverse_fft(void) { arm_rfft_fast_f32(fftInst, fftOut, recon, 1); }这里把fftOut作为输入、recon作为输出。注意传入的频域数据必须遵循正向FFT输出时相同的紧凑布局否则还原出来的时域波形的相位和幅度都是错的。库内部在ifftFlag为1时会按逆变换的方式处理这个数据并且在结果里自动做了1/N归一化。也就是说如果你对input做FFT后不做任何修改立刻IFFT还原出来的recon理论上应该和input几乎一样。3.5 结果验证怎么判断还原对不对还原结果是否准确不能只看一眼波形就完事。我会打印前10个点的原始值和还原值做对比void verify_reconstruction(void) { printf(idx original reconstructed error\n); for (int i 0; i 10; i) { float32_t err input[i] - recon[i]; printf(%3d %8.4f %8.4f %8.4f\n, i, input[i], recon[i], err); } }如果FFT和IFFT链路正确前几个点的误差应该很小通常在1e-6量级这算是浮点运算正常误差。如果误差大先查FFT点数配置、缓冲区和布局是不是对的。走完把这个滤波流程之后我会对整个FFT_SIZE个点做一次误差统计计算RMS误差来量化滤波前后信号的接近程度。120Hz被清掉后还原出来的recon应该明显不同于原始input也就是120Hz成分已经消失了剩下的波形是一个干净的50Hz正弦波。4. 实测与避坑文档里不会明说的底层问题4.1 FPU没开导致诡异结果F407的浮点运算依赖FPU但很多人做FFT之前会忽略FPU配置。症状是FFT一旦跑起来程序直接HardFault或者算出来的结果全是一堆莫名其妙的大数。Keil工程默认只有勾选了库文件和PACK选项没有主动开启FPU编译选项时生成的代码会用软浮点库性能和正确性都受影响。如果是用STM32CubeMX生成的工程在SystemInit或main开头通常会调用SCB-CPACR | ((3UL 10*2) | (3UL 11*2));来使能FPU访问控制但不同IDE这一步不完全统一。我在GCC环境下踩过坑编译命令缺了-mfloat-abihard导致FPU寄存器操作不到FFT结果稳定性极差。最稳妥的办法是在main函数初始化阶段显式使能FPU同时检查编译器的浮点ABI选项。4.2 采样率、频率分辨率和整周期采样FFT的频率分辨率等于采样率除以FFT点数。想要分辨1Hz的频率差异就得让采样率和点数的商等于1。比如用1024Hz采样率采1024点频率分辨率就是1Hz如果采样率改成2048Hz仍采1024点分辨率就变成2Hz。这个分辨率决定了你在频域能区分多近的两个峰。采样的整周期问题也很关键。如果信号频率不是频率分辨率的整数倍能量会泄漏到相邻的bin形成梳状裙边。我做合成信号时故意把50Hz和120Hz都取成整数bin就是为了先避开这个问题。真实ADC采样时很难保证整周期解决方案通常是加窗函数比如汉宁窗、布莱克曼窗。但插播一句如果后面要做IFFT还原并保持幅值准确加窗就必须在还原后除以窗函数增益否则幅值会被压低。在纯滤波场景里更简单的做法是提高频率分辨率、增长采样时间让泄漏范围变小然后滤完在时域上做必要的幅值补偿。4.3 内存分配和栈溢出FFT处理涉及大数组。1024点FFT输入缓冲区是1024个float占4KB输出缓冲区同样4KB再加上arm_rfft_fast_instance_f32内部的旋转因子表加起来也是KB级别。如果这些数组定义在函数内部局部变量里栈空间很容易被干爆然后随机HardFault。我的习惯是把FFT相关缓冲区和实例结构体全部定义成全局变量或静态变量。不要在初始化函数里循环调用arm_rfft_fast_init_f32它是重量级操作正常一次就够了。对于2048点以上FFT注意F407的RAM虽然总共192KB但有部分区域是CCM RAM能否被DMA访问需要区分。普通FFT计算用主RAM就好CCM留给中断栈这类要求低延迟的数据。4.4 幅值规格化别被FFT输出数字误导新手最容易问的一个问题是“我输入一个幅度为1.5V的正弦波为什么FFT算出来峰值是768”因为FFT结果各bin的数值和参与运算的点数N强相关。DC分量输出约等于直流平均值乘以N非DC分量的输出约等于信号幅度乘以N/2。要还原出物理幅值必须做规格化float32_t get_physical_amplitude(int k, float32_t rawMag) { if (k 0 || k FFT_SIZE / 2) { return rawMag / FFT_SIZE; } else { return 2.0f * rawMag / FFT_SIZE; } }如果不做这一步拿FFT结果做阈值判断就可能得出完全错误的结论。我调试振动阈值报警时就是靠这个换算把FFT幅值和振动加速度传感器的真实幅值对上了。4.5 IFFT输入布局错位导致还原失真做IFFT之前修改频域数据时很容易把数组中元素的下标和频率bin搞混。比如有人想清除第120个bin直接写了fftOut[120] 0但按第2章的布局第120个bin的实部在fftOut[240]虚部在fftOut[241]。这个错误不会导致程序崩溃但还原出来的波形容貌全非而且很难排查。为了避免这个问题我专门封装了设置bin值的函数void set_bin_complex(float *fftData, int fftSize, int k, float re, float im) { if (k 0) { fftData[0] re; } else if (k fftSize / 2) { fftData[1] re; } else { fftData[2 * k] re; fftData[2 * k 1] im; } }所有频域修改都走这个函数彻底避免下标算错。还是那句话宁可多写一个函数不要裸操作下标。5. 从Demo到工程接ADC和DMA的实时方案5.1 ADCDMA采集的完整链路合成信号验证完之后要处理真实ADC数据可以用定时器触发ADC采样DMA自动把结果搬进缓冲区采满一帧触发一次中断主循环里做FFT。这是嵌入式信号处理最经典的采集方案。配置要点如下ADC采样率由定时器触发频率决定。想得到第k个bin对应的频率用公式freq k * ADC_SAMPLE_RATE / FFT_SIZE反推ADC采样率。DMA传输完N个样本后触发半传输和全传输中断可以用来做双缓冲。我的做法是DMA半传输完成时处理前半段数据全传输完成时处理后半段数据这样采集和处理能流水线重叠。ADC数据左对齐或右对齐会影响数值范围。F407的ADC默认是12位右对齐读到的是0到4095之间的整数要转成float并减去2048才能去掉直流偏置。uint16_t adcRaw[FFT_SIZE]; volatile uint8_t adcDoneFlag 0; void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc) { if (hadc-Instance ADC1) { adcDoneFlag 1; } } void adc_buffer_to_float(void) { for (int i 0; i FFT_SIZE; i) { input[i] ((float32_t)adcRaw[i] - 2048.0f) / 2048.0f; } }把ADC原始码值减去2048并归一化到-1到1之间目的是去掉直流偏置也让后续FFT的数值范围保持合理。这个步骤很多人会漏掉结果频谱0Hz那里冒出一个高得离谱的直流尖峰。5.2 采样率与实时性怎么权衡做实时FFT时采样率和FFT计算耗时必须一起考虑。举个例子用8kHz采样率采1024点采满一帧信号需要128ms。而FFT计算只需要不到0.3ms占不足1%。这种情况下FFT本身不是瓶颈瓶颈反而是采样等待时间。你可以考虑提高采样率或缩短点数来提升响应速度。如果你用48kHz采样率做1024点FFT帧时长约21ms。要做实时音频频谱显示21ms一帧是能接受的但如果你用FFT/IFFT做实时降噪耳机21ms的块延迟对于对讲场景就太明显了。这种时候要缩短FFT点数或者采用重叠帧处理overlap-add/overlap-save。我在高频振动监测里用的是256点FFT帧时长远低于实时要求效果很好。5.3 在Proteus仿真和真板调试的差别我看到热搜里有“proteus没有stm32f407怎么办”。Proteus对F407的支持确实没有F103那么全很多老版本里没有对应的MCU模型。我的建议是如果主要是验证FFT/IFFT算法逻辑不关心FPU性能差异可以先在Proteus里用STM32F103模型跑同样的CMSIS-DSP库代码。F103没有FPU浮点运算靠软件模拟速度慢但算法流程是完全一样的。更推荐的做法是直接上真板。现在很多F407最小系统板几十块钱就能买到配合CMSIS-DSP库和ST-Link烧录调试比在仿真软件里猜问题高效得多。把FFT结果通过串口打印到PC用串口助手或Python脚本画频域图整个过程调试效率非常高。我在真板上调试时发现过一个问题F4的ADC采样速度很高但如果在DMA中断里直接做浮点转换和FFT可能会影响ADC采集时序导致采样间隔不均匀频谱出现杂散。解决办法是把ADC中断里只置标志位浮点转换和FFT放到主循环里执行。这个细节在仿真环境里几乎发现不了真板上一测就现原形。6. 几个让我印象深刻的扩展玩法6.1 用FFT/IFFT做简易频域滤波器我在第3章里演示了清除单个bin的滤波器实际工程中更常用的是频域掩码方式一次性把某个频段整体置零。比如做振动监测时不对电机基频下手只清除200Hz以上的高频噪声相当于低通滤波器。把1kHz到512Hz以外的bin全部清零后IFFT还原出的时域波形变得非常干净。这个方法比设计同阶数的FIR滤波器直观太多尤其适合你确实只需要处理某个固定频段的情况。同样要注意的是频域掩码做滤波会引入块效应就是帧与帧之间的信号会有拼接跳变。要缓解块效应建议采用相邻帧重叠50%再用窗函数处理后叠加这也是很多商用音频软件里所谓“频谱滤波”的标准套路。这块我还在优化但基础链路已经能用了。6.2 频谱特征提取用于故障判断除了信号还原FFT本身还有一个重要用途提取故障特征。我在电机监测里把FFT幅值谱计算出来后遍历所有bin找峰值再把峰值频率和轴承故障特征频率公式算出来的理论值做对比误差在设定范围内就认为有对应故障。整个判断过程在MCU里只需要几毫秒完全满足在线监测要求。因为arm_rfft_fast_f32的输出已经在内存里了写一个峰值搜索函数非常容易。我习惯把幅值谱保存一份到全局变量然后用滑动比较找局部最大值再做阈值判断。实际用下来比很多上位机分析实时性高得多。6.3 IFFT还原出的波形还能干什么IFFT信号还原的价值不只是滤波它还能用于波形预测、信号重构、系统辨识。比如激励信号采样后做FFT在频域里把系统幅频、相频特性修正一下再做IFFT得到修正后的时域输出这种思路在声学补偿、振动主动控制里很常见。我最近在玩一个音频相关的小项目用FFT/IFFT把采集到的发动机噪声分频段处理对特定频段做增益调整再还原回时域输出给功放实现一种简易的主动降噪模型。虽然正式产品会用更成熟的算法但用这个框架做原型验证非常快改频域增益就是改几个数组元素的事不需要重新设计时域滤波器。6.4 继续深入的学习方向如果你把arm_rfft_fast_f32摸透了接下来可以往这几个方向深入其一试试定点的arm_rfft_q31和arm_rfft_q15它们在低端MCU上同样高效数值范围需要注意其二研究arm_cfft_f32做复数FFT当你的信号是I/Q解调后的复数序列时必须用复数FFT其三把FFT和DMA、双缓冲、RTOS线程结合起来做一个自动化的实时频谱分析模块。这些都是CMSIS-DSP库自带的函数官方文档和例程都写得很清楚有前面这些基础后上手会很快。我个人在实际项目里的体会是FFT/IFFT这套东西一开始看起来像数学难题但用CMSIS-DSP之后它最难的其实是数据布局和工程配置而不是算法本身。把arm_rfft_fast_f32这个函数的输入输出格式吃到脑子里去把FPU、优化等级、缓冲区这三个坑提前填上剩下的就是业务逻辑的事。最后再分享一个小技巧在新板子上验证FFT/IFFT移植是否正确时别一上来就接传感器先用一个合成正弦波在代码里自测把FFT和IFFT往返的误差打印出来小于1e-4就说明整个链路是通的再往上接真实信号会顺畅得多。