STM32-FFT代码说明

FFT代码说明FFT为FastFourierTransformation，即快速傅里叶变换，本项目中，FFT的目标是识别频率为形如式(1.1)的一个正弦信号：()sin(2)fxAft(1.1)其中，21Hz5%f；因为单片机通过ADC接口读取该正弦信号的电压值，而12位精度的ADC的值范围在0-4096之间，如信号经过放大器后映射到0-3.3V之间，则振幅A的取值0-4096之间。假设，信号经过放大器后，其电压值最大为3.3V，最小为0V，则此信号的振幅为1.65V，对应A=2048，即该信号为：()20482048sin(221)fxt(1.2)本文中给出的例程即通过FFT识别式(1.2)这种正弦信号。假设采样频率为Fs，信号频率Fn，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：1*/FnnFsN(1.3)频率分辨率（/fFsN）等于采样时间的倒数。例如要分辨0.1Hz，则需要采集10s。1基于STM32官方DSP库的FFT算法工程文件中包含三个函数库，分别为：cr4_fft_64_stm32.scr4_fft_256_stm32.scr4_fft_1024_stm32.s分别对应数据点数为64,256和1024时的FFT算法。下面将以数据点数是1024为例，说明FFT的实现过程。1.1通过函数生成原始数据for(i=0;iNPT;i++){Fx=2048+2048*sin(PI2*i*21/Fs+20)+1000*sin(PI2*i*15/Fs)+100*sin(PI2*i*25/Fs)+rand()%100;//振幅2048，频率21HZ，为主要的谐波分量，21HZ也是所需信号的频率，rand()%100为0-100之间的随机噪声lBUFIN[i]=((s16)fx)16;//高位为实部，低位为虚部}Fs为采样频率，此处设为102.4Hz，NPT为数据长度（即前文中提到的N），为了能够分辨0.1Hz，数据长度NPT设为1024，则频率分辨率=102.4/1024=0.1Hz。1.2对原始数据进行FFTSTM32库中给出的FFT函数格式如下：Voidcr4_fft_1024_stm32(lBUFOUT,lBUFIN,NPT);其中，NPT为数据点个数，通过修改宏定义来修改NPT代表的值，此处为1024；LBUFOUT[NPT]和LBUFIN[NPT]均是长度为NPT的32位长整型的数组，高16位为实部，低16位为虚部，LBUFOUT[NPT]用于保存FFT之后的输出值；LBUFIN[NPT]用于保存输入的数值。对原始数据进行FFT之后，需要对输出数据进行处理。1.2.1频率计算每个频率点处代表的真实频率为for(i=0;iNPT/2;i++){Fn=i*Fs/NPT//由于此处i是从0开始的，所以不需要再减1}1.2.2幅值计算经过FFT后，每个频率点处的真实幅值的计算公式如下。0(LBUFOUT[0])/N(LBUFOUT[i])2/N0iAabsAabsi(1.4)voiddsp_asm_powerMag(void){s16lX,lY;u32i;for(i=0;iNPT/2;i++)//由于FFT的频谱结果是关于奈奎斯特频率对称的，所以只计算一半的点即可{lX=(lBUFOUT[i]16)16;//取低16位，虚部lY=(lBUFOUT[i]16);//取高16位，实部{floatX=NPT*((float)lX)/32768;floatY=NPT*((float)lY)/32768;floatMag=sqrt(X*X+Y*Y)/NPT;lBUFMAG[i]=(u32)(Mag*65536);}//真实振幅lBUFMAG[i]=sqrt(lX*lX+lY*lY)*2/NPT，先除以32768，又乘以65536是为了符合浮点数的计算规律}lBUFMAG[0]=lBUFMAG[0]/2;//直流分量不需要乘以2，所以按照乘以2计算后，要除以2}得到每个频率点所代表的真实频率Fn，以及该真实频率所代表的原始信号的振幅A，将其通过串口打印到电脑上。for(i=0;iNPT/2;i++){printf(%4d,%.2f%10d\n,i,((float)i*Fs/NPT),lBUFMAG[i]);}结果如下图所示。图1f=0hz处幅值图2f=15Hz处幅值图3f=21hz处幅值图4f=25hz处幅值2非官方FFT算法另外程序中给出了非官方FFT计算方法，文件名为FFT2.c，相关的函数如下：dsp_g2_init();//用于生成计算所需原始数据，dsp_g2_test();//非STM32库中的FFT算法printf_dsp_g2_Mag();//将各个频率点上的振幅打印到电脑上timetest_g2FFT();//通过定时器测试该FFT算法消耗的时间，包括计算各个频率点的振幅计算时间非官方的FFT算法的计算结果详见”非官方FFT结果.txt”附件3耗时测试通过配置TIM6定时器进行计时，单位为ms。对官方FFT算法以及非官方FFT算法进行了耗时测试，测试10次FFT总时间，求取平均值，其中包括计算各个频率点幅值的时间，将测试结果打印到电脑上。如图5所示，官方FFT算法耗时9ms，非官方耗时411ms。图5两种算法所消耗时间4结论官方FFT算法与非官方FFT算法数据点为1024时，幅值精度均很高。但是官方FFT算法很快，仅耗时9ms，非官方库耗时较长，为411ms。