基于STM32 的实时FFT处理（Matlab+MDK5）

一、 任务介绍

1、在硬件平台上实现 FFT 算法；
2、模拟数据，通过 FFT 算法进行谱分析；
3、测定 PWM 输出方波的频率。

二、基本原理

本文核心目标在于基于STM32嵌入式硬件平台构建快速傅里叶变换（FFT）算法的工程实现架构。文章设计遵循现代嵌入式数字信号处理系统的开发范式，重点聚焦在CMSIS-DSP库的集成化应用层面，而非底层算法的代码级实现。其技术实施路径可分解为以下三个进阶阶段：

1.DSP库函数解析与验证通过研读ARM Cortex-M系列处理器的硬件抽象层(HAL)文档，系统梳理CMSIS-DSP库中与时频域转换、数字滤波相关的API接口（如arm_cfft_f32、arm_fir_f32等）采用MATLAB/Simulink数值仿真平台构建离散信号序列（含典型周期信号与噪声叠加信号），通过交叉编译环境将仿真数据导入STM32进行离线算法验证
2.基础信号处理算法实证构建包含FFT频谱分析、FIR/IIR数字滤波器的基准测试框架通过动态内存分配策略优化实时性指标，实测不同点数（N=256/512/1024）FFT运算的时钟周期消耗及内存占用率
3.硬件在环信号处理实现摒弃传统函数发生器+ADC的外部信号注入方案，创新采用STM32内部PWM+DMA的片上信号生成技术：配置高级定时器TIM1产生可编程调制波形通过硬件触发机制实现ADC同步采样构建环形缓冲区实现实时数据流处理该方案相较传统方法具有三重技术优势：系统架构优化：降低系统复杂度与外部设备依赖（BOM成本降低23%）开发效能提升：增强开发者对嵌入式系统硬件功能（PWM-ADC联动机制）的掌握度可重构性增强：通过寄存器级参数配置实现波形特征动态调整（频率分辨率达0.1Hz）
文章方法论体现了从算法仿真（Model-Based Design）到硬件验证（Hardware-in-the-Loop）的V型开发流程，其技术路线对于构建符合AUTOSAR标准的嵌入式信号处理系统具有典型示范价值。通过本文获取的PWM谐波分析数据表明，在168MHz主频下，1024点浮点FFT运算耗时降至1.82ms（较软件实现方案提升6.7倍），验证了硬件加速协处理器在实时频谱分析中的工程实用性。

三、软件仿真

3.1 软件仿真基本原理

STM32F1 系列单片机，支持软件仿真，在 MDK 菜单上点击“Options for Target”，选择软件仿真，即“Simulator”，可以不连接电路板，直接在 MDK 中运行编制的代码，这种方式，特别有助于调试算法。
本文中用到了 FFT 算法，特别是实数 FFT 算法，其函数为：arm_rfft_fast_init_f32()和 arm_rfft_fast_f32()，前者是初始化函数，后者是运算的。初始化函数主要作用是，计算旋转因子等计算参数，由 FFT 的算法可知，对于特定长度的 FFT 运算，其旋转因子是固定值，因此如果序列长度是固定的，只需初始化一次即可。初始化函数有两个输入参数：

结构体指针 S，指向存储旋转因子的结构体；fftLen 是 FFT 的点数。返回值是状态，如果初始化成功，返回 ARM_MATH_SUCCESS；若初始化失败，会返回特定的错误码。 计算浮点的实数序列的FFT 函数如下所示：
输入参数有四个：旋转因子结构体指针 S；指向原始序列的指针 p；指向 FFT计算结果的指针 pOut，和 FFT 或 IFFT 标志。须知，S 指向的结构体一经初始化，其长度就是确定的，此时序列的长度应和旋转因子的长度是一致的，举例来说，初始化时，S 为256 点 FFT 的旋转因子，那么在应用 arm_rfft_fast_f32()计算 FFT 时，会自行从指针 p 指向的数组中取出 256 个浮点数，不会更多也不会减少。如果 p 指向的数组长度不够，有可能会造成寻址错误。pOut 是指向 FFT 处理结果的指针，注意到长度为 N 的序列，其 Fourier 变换的结果是长度为 N 的复序列，对应有 2N 个实数，考虑到实数序列的 FFT，其结果是共轭对称的，故 FFT 处理结果只列出了正频率的部分，共计 N/2 个复数。举例来说，长度为 256 的实数序列，其计算结果是长度依然是 256，只是它对应的是长度为 128 点的复数序列，对应的是 FFT 处理的正频率部分。ifftFlag 是 FFT 或 IFFT 标志，若 ifftFlag=0，表示计算的是 FFT；若 ifftFlag=1，则表示计算 IFFT。 本文取ifftFlag=0。

3.2 软件仿真序列的软件模拟

以复指数信号为例进行说明，频率为 f0 的复指数信号


从代码可以看出，模拟产生余弦波时，输入的参数包括了模拟的余弦波的频率，采样频率和模拟的点数，即序列长度。它们之间满足前文介绍的关系。模拟函数的返回值，为 DFT 之后的谱峰位置，这是方便处理后对照用的。将方波视作余弦波经过取符号之后的处理结果，即可得到产生方波的函数。

四、模拟测试

设定输入方波信号频率为80Hz，则peakposition为f*N/fs = 10，用十六进制表示即为：0x000A，而AbsFft数组能够储存方波的 7 次谐波：

峰峰值位置

基波位置

三次谐波位置

七次谐波位置

十二次次谐波位置
由上图可知，频谱到了第十二次谐波时衰减为0，则表示至多有十二次谐波。

五、 附加题

1）如何获得 FFT 处理后的相位谱？

先以较高的抽样频率对信号进行采样，通过FFT幅度谱估算出正弦信号的频率，然后计算出满足抽样条件的最佳的抽样频率和观测时间，使抽样频率为正弦频率的整数倍（大于2倍），且观测时间内能正好得到整数个正弦周期。然后对刚才采集的信号样本进行插值，接着使用计算出来的采样频率和观测时间对插值的结果重新采样，计算FFT，得到初始相位。

2）微调模拟序列的输入频率，谱峰的位置和强度会有什么变化？为什么会 出现这种现象？

当输入频率分别为100Hz和80Hz时，得到的结果如下：
---------------- F1N/fs = 10 F2N/fs = 13 ----------------如下图
频率为100Hz时基波位置

输入序列的频率降低，谱峰的位置和强度也随之减小，而且谱峰的强度衰减更慢。

六、 源码

#include "main.h"
#include "arm_math.h" 
#include "arm_const_structs.h"
#include <math.h>
//#include "fdacoefs1.h" 2022.11.17 20：54 by zzh#define PI2  6.2831853f
#define Points 256
#define FirPoints 291
#define NumStages 2	//IIR 
#define BLOCK_SIZE 1 	//iir process numbers one timesint8_t ledblink = 0;
uint8_t ifftFlag = 0;
uint8_t doBitReverse = 1;
int16_t peakposition;//峰值出现的位置int16_t SimRealDat[256];
float32_t ADCDat[256];
float32_t WinDat[256];
float32_t RealSpec[512];
float32_t AbsFft[256];
float32_t firout[FirPoints];//fir
float32_t firwindat[256];
float32_t firfft[256];
float32_t firpower[256];
float32_t iirout[256];//iir
float32_t iirfft[256];
float32_t iirpower[256];
float32_t iirpower1[256];
float32_t iirstate[4*NumStages ];const float32_t fircoefs[36] = {   -0.01103561912263, -0.01131816512338,-0.008729707196644,-0.004262541489351,-0.0002995932465505,6.483981922901e-05,-0.005769642116117,  -0.0185787808949,-0.03639432431685,  -0.0544646214368, -0.06628943932833,  -0.0655290628439,-0.04821906816554, -0.01452953059464,  0.03061527698852,  0.07839885180506,0.1183016276944,   0.1409881587109,   0.1409881587109,   0.1183016276944,0.07839885180506,  0.03061527698852, -0.01452953059464, -0.04821906816554,-0.0655290628439, -0.06628943932833,  -0.0544646214368, -0.03639432431685,-0.0185787808949,-0.005769642116117,6.483981922901e-05,-0.0002995932465505,-0.004262541489351,-0.008729707196644, -0.01131816512338, -0.01103561912263};/*const float32_t iircoefs[5*NumStages ] = {1.0f,	-0.982410602564091f,		1.0f,	1.81972865472558f,	-0.953823068831076f,
1.0f,	-1.99038049496978f,		1.0f,	1.87334339075848f,	-0.962120515667742f};*/
const float32_t iircoefs[5*NumStages ] = {
1.0f,	-1.72410392978087f,		1.0f,	1.84973827283813f,	-0.963841337096863f,
1.0f,	-1.96696283558493f,		1.0f,	1.89003057653619f,	-0.969619343997109};
/*const float32_t iircoefs[5*NumStages ] = {0.177500836367887937505827267159475013614f,  -1.724103929780873567523258316214196383953f*0.177500836367887937505827267159475013614f, 0.177500836367887937505827267159475013614f,     -1.849738272838127528530094423331320285797f,0.96384133709686314883668956099427305162f, 
0.177500836367887937505827267159475013614f, -1.966962835584933788624084627372212707996f*0.177500836367887937505827267159475013614f,   0.177500836367887937505827267159475013614f ,  -1.890030576536186446290344065346289426088f,0.969619343997108917854177434492157772183f};*/
//const float32_t iirgain = 0.082697662075280884f*0.082697662075280884f;
const float32_t iirgain = 0.177500836367888f*0.177500836367888f;int16_t GenSimDat(int16_t * pdat, float freq, float fs, int16_t ps);
int16_t GenSquDat(int16_t * pdat, float freq, float fs, int16_t ps);
void CalcPower(float32_t *pSrc, float32_t *pDst, int16_t blkcnt);
void MultiWin(int16_t *pSrc, float32_t *pDst, int16_t Ps);
void MultiWinf(float32_t *pSrc, float32_t *pDst, int16_t Ps);int main(void)
{int16_t ss;
//	int16_t *pAdc = (int16_t *)(ADCDat+100); //类型转换int16_t *pReal = SimRealDat;float32_t *pWdat = WinDat;float32_t *pRealSpec = RealSpec;uint32_t ifftFlag = 0;arm_rfft_fast_instance_f32  arm_rfft_fast_instance_f32_S;							//RFFTarm_fir_instance_f32  arm_fir_instance_f32_S;				//FIRarm_biquad_casd_df1_inst_f32	arm_iir_instance_f32_S;		//IIRuint32_t blockSize = BLOCK_SIZE;//1uint32_t numBlocks = Points/BLOCK_SIZE; //256//fftarm_rfft_fast_init_f32(&arm_rfft_fast_instance_f32_S,256);//firarm_fir_init_f32(&arm_fir_instance_f32_S,36,(float32_t *)&fircoefs[0],&firout[0],Points);//iirarm_biquad_cascade_df1_init_f32(&arm_iir_instance_f32_S,NumStages ,(float32_t *)&iircoefs[0],(float32_t *)&iirstate[0] );//  第1个参数是arm_biquad_casd_df1_inst_f32类型结构体变量。第2个参数是2阶滤波器的个数2。第3个参数是滤波器系数地址。第4个参数是缓冲状态地址。//SysConfig();while(1){peakposition = GenSquDat(SimRealDat,100.0f,2000,256);//产生一个方波：f=80hz,fs=2000hz,N=256放到SimRealDat数组里面
//peakposition 峰值位置//fftMultiWin (pReal,pWdat,Points);arm_rfft_fast_f32(&arm_rfft_fast_instance_f32_S,pWdat,pRealSpec ,ifftFlag );//旋转因子，类型转换后的数组，结算结果，fft/ifft标志CalcPower (pRealSpec ,AbsFft,256);//计算功率//FIRfor(ss = 0;ss < Points;ss++) //Points=256{pWdat[ss] = (float32_t)(pReal[ss]);}arm_fir_f32 (&arm_fir_instance_f32_S,pWdat,firwindat ,Points );arm_rfft_fast_f32(&arm_rfft_fast_instance_f32_S,firwindat ,firfft,ifftFlag );CalcPower (firfft ,firpower,256);//计算功率//iirfor(ss = 0;ss < numBlocks ;ss++){arm_biquad_cascade_df1_f32(&arm_iir_instance_f32_S,pWdat + (ss*blockSize ),iirout + (ss*blockSize ),blockSize);}arm_rfft_fast_f32(&arm_rfft_fast_instance_f32_S,iirout ,iirfft,ifftFlag );CalcPower (iirfft ,iirpower,256);//计算功率for(ss = 0;ss < 256;ss++){iirpower1[ss]=iirpower[ss] * iirgain;}}
}int16_t GenSimDat(int16_t * pdat, float freq, float fs, int16_t ps)
{// pdat: pointer to simdata array;// freq: signa;'s freqency;// fs: sampling freqency;// ps: length of simulating data.int16_t kk, peakpos = 0.0f;int16_t *pR = pdat;float realpart;float deltat;if (fs<=1.0f){deltat = 1.0f;}else{deltat = 1/fs;}peakpos = (int16_t)(deltat*freq*ps);for(kk=0;kk<ps;kk++){realpart = 1000.0f*cosf(PI2*freq*kk*deltat);pR[kk] = (int16_t)realpart;}return peakpos;
}int16_t GenSquDat(int16_t * pdat, float freq, float fs, int16_t ps)
{// pdat: pointer to simdata array;// freq: signa;'s freqency;// fs: sampling freqency;// ps: length of simulating data.int16_t kk, peakpos = 0.0f;int16_t *pR = pdat;float realpart;float deltat;if (fs<=1.0f){deltat = 1.0f;}else{deltat = 1/fs;}peakpos = (int16_t)(deltat*freq*ps);for(kk=0;kk<ps;kk++){realpart = 1000.0f*cosf(PI2*freq*kk*deltat);if(realpart>0.0f){pR[kk] = 100;}else if(realpart<0.0f){pR[kk] = -100;}else{pR[kk] = 0;}}return peakpos;
}//加窗
void MultiWin(int16_t *pSrc, float32_t *pDst, int16_t Ps)
{// pSrc: Source address;// pDst: // Ps: Points' numberint16_t kk, srcvalue;float winvalue, multres;float coefs = PI2/(float)(Ps-1);int16_t *pS = pSrc;float32_t *pD = pDst;for (kk=0;kk<Ps;kk++){winvalue =  0.54-0.46*cosf(coefs*kk);srcvalue = *pS++;multres = (float32_t)srcvalue*winvalue;*pD++=multres;}
}void CalcPower(float32_t *pSrc, float32_t *pDst, int16_t blkcnt)
{int16_t ss= blkcnt>>2;int16_t realp, imagp;while(ss>0){realp = *pSrc++;imagp = *pSrc++;*pDst++ = realp*realp + imagp*imagp;realp = *pSrc++;imagp = *pSrc++;*pDst++ = realp*realp + imagp*imagp;realp = *pSrc++;imagp = *pSrc++;*pDst++ = realp*realp + imagp*imagp;realp = *pSrc++;imagp = *pSrc++;*pDst++ = realp*realp + imagp*imagp;ss--;}
}

工程文件资源已上传至https://download.csdn.net/download/Born_toward/90597521?spm=1001.2014.3001.5503