用STM32实现FFT

前言

电子设计大赛已经过去很久了，一直想写一篇关于FFT的文章也没有来得及，现在写一下来记录分享一下。
本篇文章不讲复杂的FFT原理，只讲如何在stm32里面怎么实现FFT

---

一、FFT是什么，能干啥？

FFT(fast Fourier transform)快速傅里叶变换，把时域信号变换到频域，至于换到频域后能做啥，那就多了,你能搜到这篇文章，说明你已经知道他是干啥的了，，，，

二、先用matlab来认识一下FFT

1.上代码

代码如下：

fs = 1024;           %采样频率
NPT = 1024;          %样点个数
t=0:1/fs:1;          %采样步长,模拟单片机采样频率
y= 2+200*cos(2*pi*15*t)+40*cos(2*pi*100*t+pi*90/180)+15*cos(2*pi*250*t+pi*90/180);
plot(t,y);		     %画出未变换前的时域信号波形
fbl=fs/NPT;          %分辨率
f=(0:NPT)*fbl;       %每个点所代表的频率
Y=fft(y);		     %进行FFT变换，matlab中fft函数，输入多少采样数据，就输出多少幅频值(是个复数)
				     %所以length(y)=length(Y); 
figure(2)
plot(f(1:NPT/2),abs(Y(1:NPT/2)));
				     %以每个点所代表的频率画出其对应的赋值

2.代码关键字讲解

思想都是一样的，matab知道了，在32里面也差不多

采样频率： 用单片机ADC采样，结果是离散的数字信号，参考采样定理，采样频率至少大于信号频率的两倍，在这里我们采样频率1024,远大于基波频率2，选取1024也是为了方便自己，可以更直观的观察输出（2的整数倍都可）。
分辨率： 输出时，最小频率间隔(Fs/N)，第n个点所代表的频率即为Fn=(n-1)*Fs/N,在这个程序中，采样频率是1024，采样点数是1024，所以，分辨率为1Hz
输出赋值 fft输出就是复数，取模值即可

3.程序结果

未变换前的时域信号波形图：

转换后输出的幅频函数波形图：

由图可看出输入信号大概有三个频率，15，100，250.与输入信号相符

二、淦! 上32

STM32官方给了三个函数库文件，分别是：
cr4_fft_64_stm32.s
cr4_fft_256_stm32.s
cr4_fft_1024_stm32.s
（这些都是汇编语言写的，好在我们不需要看懂，只需要会调用就好。）
分别对应这采样点数为64，256，1024。库函数的获取，可以在官网上获取
我们的目标是计算周期为1ms波形的THD值！，使用采样点数为64位的函数

1、代码移植

新建文件夹，把
cr4_fft_64_stm32.s
cr4_fft_256_stm32.s
cr4_fft_1024_stm32.s
和stm32_dsp.h
文件添加到工程里（就像添加.C文件一样），注：stm32_dsp.h文件里第27行需要根据自己单片机型号来修改

2，ADC采集

ADC采集没有什么好说的，我想让大家注意一点就是采样周期，在这里，为了方便使用，我们设ADC采样频率为64KHz，注意，这里直接调用正点原子的Get_Adc();函数的话，他的执行周期很长，需要改短（自己跳进Get_Adc(）函数里，看函数定义，再改），要不然，在定时器中，ADC采集时间会超过中断事件，导致采样周期不是64KHz,最好在中断中设置一个LED灯反转的程序，检测ADC采样周期，确保ADC采集频率稳定

3、上函数

在上函数之前，我们把函数的每一个输入输出变量的含义说一下
在我的程序里是这样调用的：cr4_fft_64_stm32(output,input,NPT);
她的输入变量是input：是一个int32_t类型的数组，高十六位是实部，第十六位是虚部，显然，我们ADC采集只能采集实数。所以我对采集到的数据进行这样的处理：

	 ADC_Value[i]=Get_Adc(ADC_Channel_5);
	 input[i] = ((s16)ADC_Value[i])<<16;

至于虚部，就不用管啦，移位后默认是0
他的另一个输入是NPT ，采样点数。最好是定义成宏定义，因为要用到他定义的数组太多，这里是64
他的输出是output,也是个复数，高十六位代表实部，低十六位代表虚部，他们的模值就是我们要的幅值
呐，官方给了专门求赋值的函数

void dsp_asm_powerMag(int32_t *IBUFOUT)
{
	s16 lx,ly;
	u32 i;
	for(i=0;i<NPT/2;i++)//由于FFT的频域结果是关于奈奎斯特频率对称的，所以只计算一半
	{
		lx = (IBUFOUT[i] << 16)>>16;      //取低十六位，虚部
		ly = (IBUFOUT[i] >> 16);          //取高十六位，实部
		float X = NPT*((float)lx)/32768;
		float Y = NPT*((float)ly)/32768;
		float Mag = sqrt(X*X+Y*Y)/NPT;
		OUTPUT_MAG[i] = (u32)(Mag*65536);
	}//这些就是计算振幅IBUFOUT[i] = sqrt(lx*lx+ly*ly)*2/NPT
	 //之所以先除以32768后乘以65536是为了符合浮点数的计算规律
	OUTPUT_MAG[0] = OUTPUT_MAG[0]/2;//这个是直流分量，不需要乘以2
}

直接用就很OK了

	dsp_asm_powerMag(output);		

那么我们调用了函数，也转换成幅值了，所以，我该怎么用这些输出的数呢。他输出的是什么***，接下来就需要对应起来，
我们需要知道的是 每个频率点对应的幅值，现在有幅值，怎么把他对应的点找到，这样我们就大功告成了
还记得前面讲的分辨率嘛，对，我们采样点是64个，采样频率也是64KHz,那么分辨率就是fs/NPT即为1KHz，也就是说输出的output每隔一个对应着一个相差位1Hz的频率，那么output[1]就对应着1KHz，output[2]就对应着2KHz，output[3]就对应着3KHz，，，，我们只取前NPT/2个，不要问为什么，问奈奎斯特吧。
还有就是output[0]对应着直流分量

我只用了前几个意思意思。

	cr4_fft_64_stm32(output,input,NPT);
	dsp_asm_powerMag(output);		

	LCD_ShowxNum(100,0,    OUTPUT_MAG[0]  ,7,16,0);	
	LCD_ShowxNum(100,16,   OUTPUT_MAG[1]  ,7,16,0);		
	LCD_ShowxNum(100,32,   OUTPUT_MAG[2]  ,7,16,0);		
	LCD_ShowxNum(100,48,   OUTPUT_MAG[3]  ,7,16,0);		
	LCD_ShowxNum(100,64,   OUTPUT_MAG[4]  ,7,16,0);		
	LCD_ShowxNum(100,80,   OUTPUT_MAG[5]  ,7,16,0);	
    THD=100*(float)sqrt((Voltage[4]*Voltage[4]+Voltage[2]*Voltage[2]+Voltage[3]*Voltage[3]+Voltage[5]*Voltage[5]))/(float)(Voltage[1]);

总结

YSJX