如何使用STM32提供的DSP库进行FFT呢

　　最近，因为项目需要在STM32F103系列处理器上，对采集的音频信号进行FFT运算，然而STM32F103毕竟不是STM32F4系列的处理器，对于一般的FFT运算程序还是比较缓慢的。
　　幸亏官方提供了针对FFT的官方库，但是去官网找了半天居然找不到那个库的下载，好像官方早就把那个库下架了，估计是为了给带DSP指令集的F4系列让路。然后就只好从别人的项目中把这个官方库给扒出来了……
　　FFT的意义
　　对于很多人来说，采样频率和FFT点数之间的关系可能还是不太清楚。下面就来简单分析一下：
　　根据采样定理，采样频率必须是被采样信号最高频率的2倍。比如，需要采集音频信号，并且需要被观察到的音频频率的频率范围是20Hz到20KHz，那么使用的采样频率就必须大于40kHz。如果需要观察到的音频频率范围为0Hz到600Hz，那么使用的采样频率只需要大于1200Hz即可。
　　而FFT点数与采样频率之间有什么关系呢？本质上并没有什么关系，但是FFT点数的大小直接关系到频率分辨率。怎么来说呢？
　　假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。同时，FFT后的N个点，开始的那个点表示直流分量（即0Hz），而最后的那个点的再下一个点表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。即，某点n所表示的频率为：Fn=（n-1）*Fs/N。这就表示，Fs/N就是频率分辨率。
　　FFT官方库的使用条件
　　FFT官方库在使用上并不灵活：属于基4的FFT，即FFT点数必须是4^n。也就是说，如果要做512点或2048点的FFT，那么对不起，没法使用官方库了；
　　FFT官方库的输入输出是等长的，即256点的FFT输入也必须是256点，如果你的输入小于这个长度，是没有任何性能提升的。
　　FFT官方库的使用
　　准备工作
　　下载得到STM32的DSP库之后，就可以将其添加到自己的工程项目中了。
　　其中，stm32_dsp.h和table_fft.h两个文件是必须添加的。stm32_dsp.h是STM32的DSP库的头文件。
　　另外，对于**.s文件可以有选择的添加**（用到那个添加那个即可）。由于本文只用到了256点的FFT，所以这里只添加了cr4_fft_256_stm32.s文件。
　　
　　FFT函数说明
　　进行256点的FFT，只需要调用STM32 DSP库函数中的cr4_fft_256_stm32（）函数即可。该函数的原型为：
　　void cr4_fft_256_stm32（void *pssOUT， void *pssIN， uint16_t Nbin）;
　　其中，参数pssOUT表示FFT输出数组指针，参数pssIN表示要进行FFT运算的输入数组指针，参数Nbin表示了点数。至于该函数的具体实现，因为是用汇编语言编写的，我也不懂，这里就不妄谈了。
　　需要说明的是：按照FFT官方库的说明，pssOUT和pssIN都必须是32位的数据类型，其中高16位存储实部，低16位存储虚部。对于pssIN来说，低16位存储的虚部总是为0。
　　代码实例
　　假设ADC采样的声音数据为adc_buf［NPT］，FFT运算的输入数组为lBufInArray［NPT］。由于FFT计算出来的数据是对称的，因此通常而言输出数组取一半的数据，为lBufOutArray［NPT/2］。除此之外，还定义各次谐波幅值lBufMagArray［NPT/2］。即：
　　#define NPT 256
　　uint32_t adc_buf［NPT］={0};
　　long lBufInArray［NPT］;
　　long lBufOutArray［NPT/2］;
　　long lBufMagArray［NPT/2］;
　　调用FFT官方库的代码为：
　　//填充数组
　　for（i=0;i《NPT;i++）
　　//这里因为单片机的ADC只能测正的电压 所以需要前级加直流偏执
　　//加入直流偏执后，需要在软件上减去2048即一半，达到负半周期测量的目的（需要根据具体情况来进行配置）
　　lBufInArray［i］ = （（signed short）（adc_buf［i］-2048）） 《《 16;
　　cr4_fft_256_stm32（lBufOutArray， lBufInArray， NPT）;
　　同时，计算各次谐波幅值的函数为：
　　void GetPowerMag（）
　　{
　　signed short lX，lY;
　　float X，Y，Mag;
　　unsigned short i;
　　for（i=0; i《NPT/2; i++）
　　{
　　lX = （lBufOutArray［i］ 《《 16） 》》 16;
　　lY = （lBufOutArray［i］ 》》 16）;
　　//除以32768再乘65536是为了符合浮点数计算规律
　　X = NPT * （（float）lX） / 32768;
　　Y = NPT * （（float）lY） / 32768;
　　Mag = sqrt（X * X + Y * Y） / NPT;
　　if（i == 0）
　　lBufMagArray［i］ = （unsigned long）（Mag * 32768）;
　　else
　　lBufMagArray［i］ = （unsigned long）（Mag * 65536）;
　　}
　　}