【STM32】使用STM32提供的DSP库进行FFT（附详细代码）

喜欢ヅ旅行 2022-02-02 14:23 3159阅读 0赞

最近，因为项目需要在`STM32F103`系列处理器上，对采集的音频信号进行`FFT`运算，然而`STM32F103`毕竟不是`STM32F4`系列的处理器，对于一般的`FFT`运算程序还是比较缓慢的。

幸亏官方提供了针对FFT的官方库，但是去官网找了半天居然找不到那个库的下载，好像官方早就把那个库下架了，估计是为了给带DSP指令集的F4系列让路。然后就只好从别人的项目中把这个官方库给扒出来了……

下载地址：[https://pan.baidu.com/s/1GiUJgEkQxDAk79iddXKsaA][https_pan.baidu.com_s_1GiUJgEkQxDAk79iddXKsaA] 提取码: dn2d

# 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]。

# 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文件。

![2019051215361584.png][]

## 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);
        }
    }

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly95bmd6bWlhby5ibG9nLmNzZG4ubmV0_size_16_color_FFFFFF_t_70][]

[https_pan.baidu.com_s_1GiUJgEkQxDAk79iddXKsaA]: https://pan.baidu.com/s/1GiUJgEkQxDAk79iddXKsaA
[FFT]: https://blog.csdn.net/iloveyoumj/article/details/53308142
[2019051215361584.png]: /images/20220202/4690cd10bc08472d9f563047fcb1ecc6.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly95bmd6bWlhby5ibG9nLmNzZG4ubmV0_size_16_color_FFFFFF_t_70]: /images/20220202/f62aaee2f27a4b73a13303a56b7fdaa7.png