如何在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=Get_Adc（ADC_Channel_5）;
　　input = （（s16）ADC_Value）《《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 《《 16）》》16; //取低十六位，虚部
　　ly = （IBUFOUT 》》 16）; //取高十六位，实部
　　float X = NPT*（（float）lx）/32768;
　　float Y = NPT*（（float）ly）/32768;
　　float Mag = sqrt（X*X+Y*Y）/NPT;
　　OUTPUT_MAG = （u32）（Mag*65536）;
　　}//这些就是计算振幅IBUFOUT = 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］）;