STM32实现FFT，求取幅度频谱

FFT不太对劲的理解

FFT的原理比较复杂，因为32使用FFT不用去管算法是如何运作的，我在这里就进行简单的介绍了。

因为是简单介绍，就只介绍下幅度频谱图，不考虑相位频谱图。

​ FFT可以将一个信号从时域变换到频域，比如一个1VPP的1k的正弦信号，它的时域和频域的示意图如下：

​ 频域为我们观察信号提供了一个新的视角。比如下面是1k和2k信号的叠加。

​ 从时域上看，1k+2k的波形不容易进行处理，也不好猜出来这个波形到底有什么特性（当然这个例子其实还是比较好猜测的，复杂情况就不好看了）。可是变换到频域后，特性非常的明显，处理起来就方便了。

STM32实现FFT

添加DSP库

STM32 DSP库的快速添加 基于cubemx 调用，使用DSP库_四臂西瓜的博客-CSDN博客_cubemx dsp

要采用第二种方法添加DSP，第一种方式添加的DSP库比较老，支持的FFT函数用起来不方便，这篇文章介绍的FFT函数老版本不支持。

设置ADC采集交流+串口重定向

STM32HAL ADC+TIM+DMA采集交流信号 基于cubemx_四臂西瓜的博客-CSDN博客

代码编写

DSP库里面的FFT支持32-4096，同时是$$2^n$$个整数的傅里叶变换。首先定义下面这些变量。其中只有fft_inputbuf，fft_outputbuf是用来进行FFT的。

#define FFT_LENGTH 1024

__IO uint8_t AdcConvEnd = 0;
uint16_t adcBuff[FFT_LENGTH];

float fft_inputbuf[FFT_LENGTH * 2];  
float fft_outputbuf[FFT_LENGTH];  

​ 定义完成变量后进行服务内容的书写。首先是ADC进行数据的采集，这一部分再赘述：

HAL_ADCEx_Calibration_Start(&hadc1);
HAL_ADC_Start_DMA(&hadc1, (uint32_t *)adcBuff, FFT_LENGTH);
HAL_TIM_Base_Start(&htim3);

while (!AdcConvEnd)       //等待转换完毕
	;

​ 傅里叶变换分实数和虚数两种，使用最为频繁的是虚数。我们需要把ADC采集到的数据以虚数的形式存放到傅里叶变换的输入数组fft_inputbuf。虚数的存放方式如下

​ 比如ADC采集到的第二个数据adcBuff[1]，它的数据存入到fft_inputbuf[2],它的虚部是fft_inputbuf[3]，补零。

​ 我们调用arm_cfft_f32(&arm_cfft_sR_f32_len1024, fft_inputbuf, 0, 1);来对输入数据进行傅里叶变换。它的输入参数含义如下

• arm_cfft_sR_f32_len1024为傅里叶变换结构体，1024是要运算的点数。我们在进行其他点数的计算时，比如说32个点的FFT，就可以用arm_cfft_sR_f32_len32。
• fft_inputbuf为傅里叶变换需要处理的数据的首地址
• 第三个参数0是正变换，1是反变换
• 第四个参数一般是1

​ 经过傅里叶变换后的结果仍然为复数，虚部和实部的比可以计算出频率点的相位，这个在这里不进行考虑。我们直接对复数取模。

​ 取模是实部和虚部的平方和取平均来算，不过我们没必要自己去这样写，因为DSP库为我们提供了取模函数：arm_cmplx_mag_f32(fft_inputbuf, fft_outputbuf, FFT_LENGTH);参数含义如下：

• fft_inputbuf源数据，形式为复数
• fft_outputbuf取完模后的数据，形式为实数
• FFT_LENGTH是取模的点数

​ 他们背后的运算是：$$fftoutputbuf[i]=\sqrt{fftinputbuf[i\ast 2{]}^2+fftinputbuf[i\ast 2+1{]}^2}$$

for (int i = 0; i < FFT_LENGTH; i++)
{
    fft_inputbuf[i * 2] = adcBuff[i] * 3.3 / 4096;//实部赋值，* 3.3 / 4096是为了将ADC采集到的值转换成实际电压
    fft_inputbuf[i * 2 + 1] = 0;//虚部赋值，固定为0.
}

arm_cfft_f32(&arm_cfft_sR_f32_len1024, fft_inputbuf, 0, 1);
arm_cmplx_mag_f32(fft_inputbuf, fft_outputbuf, FFT_LENGTH); 

​ 运算完成的结果还需要进行下面的处理，背后的原理就跟FFT算法有关了，这里不做解释。

​ 我们进行的是1024个点的FFT变换，那么fft_outputbuf下标0需要除以1024，其余的数除以512。

fft_outputbuf[0] /= 1024;

for (int i = 1; i < FFT_LENGTH; i++)//输出各次谐波幅值
{
    fft_outputbuf[i] /= 512;
}

最后再把运算结果打印出来即可。

printf("FFT Result:\r\n");

for (int i = 0; i < FFT_LENGTH; i++)//输出各次谐波幅值
{
	printf("%d:\t%.2f\r\n", i, fft_outputbuf[i]);
}

总的代码如下图

运行结果

​ ADC以100k的频率去采集信号发生器产生的976hz的正弦信号，信号VPP=2v，直流偏置为2V。（注意别超过内部ADC的测量范围0-3.3V）

FFT Result:
0:	1.97
1:	0.00
2:	0.00
3:	0.00
4:	0.00
5:	0.00
6:	0.00
7:	0.00
8:	0.00
9:	0.01
10:	1.05
11:	0.01
12:	0.00
13:	0.00
14:	0.00
...（后面的数据都为0）

结果分析

​ FFT计算出来的数据是对称的，我们只取前一半的数据，此时的前一半数据是512个。

​ FFT输出数组的下标表示的频率，计算关系为：
$$频率=数组下标\ast \frac{采样率}{fft计算的点数}$$
​ 利用这个公式分析下上方的运行结果。数组下标0对应的是0hz，也就是直流偏置，幅度为1.97，和输入信号的2V符合。数组下标10对应的频率为$$\frac{100k}{1024}\ast 10\approx 976.5hz$$，对应幅度为1.05V，和输入信号的2VPP相符。

精度问题

不知道读者有没有注意到待测频率为976hz，而不是我们平时常见的1k，2k？这是为了这个频率正好落在数组下标10点上。10对应的是$$\frac{100k}{1024}\ast 10\approx 976.5hz$$，11对应的$$\frac{100k}{1024}\ast 11\approx 1075hz$$，1k落在了两点中间，这样就引起了栅栏效应，给人的直观感受就是能量分散了。下面是把待测信号改成1k后的运行结果。

FFT Result:
0:	1.98
1:	0.02
2:	0.03
3:	0.03
4:	0.03
5:	0.04
6:	0.05
7:	0.07
8:	0.10
9:	0.18
10:	0.95
11:	0.30
12:	0.13
13:	0.09
...

可以看到10，11，9等下标都分到了电压（能量）。实际应用中应尽可能避免栅栏效应。

栅栏效应下的补偿

如果不可避免得碰到了栅栏效应，是可以通过数据处理尽可能的还原求取待测信号的幅度值。方法是把频率点附近的能量聚集起来，将附近频率点的幅度平方求和，再取平均。

比如上面1k的情况，就可以把8，9，10，11，12的能量聚集起来。
$$\sqrt{0.1^2+0.18^2+0.95^2+0.3^2+0.13^2}=1.026$$
经过补偿后的幅度值就跟VPP2V真实值更加接近了。

我这里只取了5个点，如果不同主要频率点下标相差比较大，可以取更多；反之更少。

后记

本文章收录于：

唐承乾的电赛小站

本文为系列文章中的冰山一角，欢迎进入小站查看。

配套程序：

STM32进行FFT傅里叶变换——0积分下载

然后，打开工程，新建 DSP_LIB 分组，并将 arm_cortexM4lf_math.lib 添加到工程里面，如图：

2， 添加头文件包含路径
添加好.lib 文件后，我们要添加头文件包含路径，将第一步拷贝的 Include 文件夹和 DSP_LIB文件夹，加入头文件包含路径。
3， 添加全局宏定义
最后，为了使用 DSP 库的所有功能，我们还需要添加几个全局宏定义：
1，__FPU_USED
2，__FPU_PRESENT
3，ARM_MATH_CM4
4，__CC_ARM
5，ARM_MATH_MATRIX_CHECK
6，ARM_MATH_ROUNDING
添加方法：点击C/C++选项卡，然后在 Define 里面进行设置，如图

这里，两个宏之间用“,”隔开。并且，上面的全局宏里面，我们没有添加__FPU_USED，因为这个宏定义在 Target 选项卡设置 Code Generation 的时候，选择了：Use FPU（如果没有设置 Use FPU，则必须设置！！），故 MDK 会自动添加这个全局宏，因此不需要我们手动添加了。
至此，STM32F4 的 DSP 库运行环境就搭建完成了。

照一定大小采样频率 F 去采集信号，采集 N 个点，那么通过对这 N 个点进行 FFT 运算，就可以得到这个信号的频谱特性。 这里还涉及到一个采样定理的概念：在进行模拟/数字信号的转换过程中，当采样频率 F 大于信号中最高频率 fmax的 2 倍时(F>2*fmax)，采样之后的数字信号完整地保留了原始信号中的信息，采样定理又称奈奎斯特定理。举个简单的例子：比如我们正常人发声，频率范围一般在 8KHz以内，那么我们要通过采样之后的数据来恢复声音，我们的采样频率必须为 8KHz 的 2 倍以上，也就是必须大于 16KHz 才行。

模拟信号经过 ADC 采样之后，就变成了数字信号，采样得到的数字信号，就可以做 FFT变换了。N 个采样点数据，在经过 FFT 之后，就可以得到 N 个点的 FFT 结果。为了方便进行FFT 运算，通常 N 取 2 的整数次方。

假设采样频率为 F，对一个信号采样，采样点数为 N，那么 FFT 之后结果就是一个 N 点的复数，每一个点就对应着一个频率点（以基波频率为单位递增），这个点的模值（sqrt(实部虚部2)）就是该频点频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为 A，那么 FFT 的结果的每个点（除了第一个点直流分量之外）的模值就是 A 的 N/2倍，而第一个点就是直流分量，它的模值就是直流分量的 N 倍。

这里还有个基波频率，也叫频率分辨率的概念，就是如果我们按照 F 的采样频率去采集一个信号，一共采集 N 个点，那么基波频率（频率分辨率）就是 fk=F/N。 这样,第 n 个点对应信号频率为：F*(n-1)/N；其中 n≥1，当 n=1 时为直流分量。
关于 FFT 就介绍到这。

arm_status arm_cfft_radix4_init_f32( 
arm_cfft_radix4_instance_f32 * S, 
uint16_t fftLen,uint8_t ifftFlag,uint8_t bitReverseFlag)  
void arm_cfft_radix4_f32(const arm_cfft_radix4_instance_f32 * S,float32_t * pSrc) 
void arm_cmplx_mag_f32(float32_t * pSrc,float32_t * pDst,uint32_t numSamples) 

第一个函数 arm_cfft_radix4_init_f32，用于初始化 FFT 运算相关参数，其中：fftLen 用于指定 FFT 长度（16/64/256/1024/4096），本章设置为1024；ifftFlag 用于指定是傅里叶变换(0)还是反傅里叶变换(1)，本次设置为0；bitReverseFlag 用于设置是否按位取反，本章设置为 1；最后，所有这些参数存储在一个 arm_cfft_radix4_instance_f32 结构体指针 S 里面。

第二个函数 arm_cfft_radix4_f32 就是执行基 4 浮点 FFT 运算的，pSrc 传入采集到的输入信号数据（实部+虚部形式），同时 FFT 变换后的数据，也按顺序存放在 pSrc 里面，pSrc 必须大于等于 2 倍 fftLen 长度。另外，S 结构体指针参数是先由 arm_cfft_radix4_init_f32 函数设置好，然后传入该函数的。

第三个函数 arm_cmplx_mag_f32 用于计算复数模值，可以对 FFT 变换后的结果数据，执行取模操作。pSrc 为复数输入数组（大小为 2*numSamples）指针，指向 FFT 变换后的结果；pDst为输出数组（大小为 numSamples）指针，存储取模后的值；numSamples 就是总共有多少个数据需要取模。
通过这三个函数，我们便可以完成 FFT 计算，求出模值后便可以求出幅值、相位，知道采样率和FFT长度就可求出频率。

fft.c 代码

#include "arm_math.h"
#include "common.h"
#include "fft.h"


/*********************************************************************************
************************STM32F407核心开发板******************************
**********************************************************************************
* 文件名称: fft.c                                                             *
* 文件简述：DSP FFT使用                                                           *
* 创建日期：2020.08.21                                                           *
* 版    本：V1.0                                                                 *
* 作    者：近视未看清人心                                                              *
* 说    明：完成FFT信号分析需要用到的函数集，包括求模，幅值、频率、相位
注意：FFT的长度，也就是采样点的个数在fft.h 宏定义中完成* 
**********************************************************************************
*********************************************************************************/	

arm_cfft_radix4_instance_f32 scfft;	//设置FFT参数的结构体变量（基4）

/****************************************************************************
* 名    称: void FFTx4_init(u8 ifftFlag,u8 bitReverseFlag)
* 功    能：FFT初始化函数（基4）
* 入口参数：ifftFlag：于指定是傅里叶变换(0)还是反傅里叶变换(1)
						bitReverseFlag：是否按位取反
                 
* 返回参数：无
* 说    明：
****************************************************************************/

void FFTx4_init(u8 ifftFlag,u8 bitReverseFlag)
{
	arm_cfft_radix4_init_f32(&scfft,length,ifftFlag,bitReverseFlag);//初始化scfft结构体，设定FFT相关参数
}




/****************************************************************************
* 名    称: void cFFTx4(float * pSrc)
* 功    能：FFT计算函数（基4）
* 入口参数：pSrc：pSrc 传入采集到的输入信
									号数据（实部+虚部形式），同时 FFT 变换后的数据，也按顺序存放在 pSrc 里面                
* 返回参数：无
* 说    明：pSrc的长度必须>=2*length
****************************************************************************/

void cFFTx4(float * pSrc)
{
	arm_cfft_radix4_f32(&scfft,pSrc);	//FFT计算（基4）
}





/****************************************************************************
* 名    称: void cmplxFFTx4(float * pSrc,float * pDst)
* 功    能：对FFT结果求模值（基4）
* 入口参数：pSrc：需要求模的复数指针（实部+虚部形式）
					 pDst：存放模值
					

* 返回参数：无
* 说    明：pSrc的长度必须>=2*length
****************************************************************************/

void cmplxFFTx4(float * pSrc,float * pDst)
{
	arm_cmplx_mag_f32(pSrc,pDst,length);	//计算复数模值

}





/****************************************************************************
* 名    称: void all_result_x4(float * pSrc,float * Amp,float * rate,float * Phase,u32 Fs)
* 功    能：计算各点幅值、频率、相位
* 入口参数：pSrc：原始信号的复数指针（实部+虚部形式）
					 Amp：存放幅值
					 rate：存放各点的频率
					 Phase：存放各点的相位
					 size：采样点的个数

* 返回参数：无
* 说    明：pSrc的长度必须>=2*size
****************************************************************************/
void all_result_x4(float * pSrc,float * Amp,float * rate,float * Phase,u32 Fs)
{
	float pDst[length];
	u32 i;
	arm_cfft_radix4_f32(&scfft,pSrc);	//FFT计算（基4）
	arm_cmplx_mag_f32(pSrc,pDst,length);	//计算复数模值
	//计算各点幅值、频率、相位
	for(i=0;i<length;i++)
	{
		if(i==0)	Amp[0]=pDst[0]/length;	//第一个点直接除以size(z直流分量)
		else Amp[i]=pDst[i]*2/length;	//其他点除以size/2
		rate[i]=Fs/length*i;//	各点频率
		Phase[i]=atan2(pSrc[2*i+1], pSrc[2*i]); /* atan2求解的结果范围是(-pi, pi], 弧度制 */
	}
	
}

fft.h 代码

#ifndef _FFT_H
#define _FFT_H
#include "arm_math.h"
#include "common.h"

#define length 1024 //FFT长度（16/64/256/1024/4096）,即采样点的个数

void FFTx4_init(u8 ifftFlag,u8 bitReverseFlag);//FFT初始化函数（基4）
void cFFTx4(float * pSrc);//FFT计算函数（基4）
void cmplxFFTx4(float * pSrc,float * pDst);//对FFT结果求模值（基4）
void all_result_x4(float * pSrc,float * Amp,float * rate,float * Phase,u32 Fs);//计算各点幅值、频率、相位


#endif

main.c 代码

#include "common.h"
#include "key.h"
#include "fft.h"
#include "usart3.h"

/*********************************************************************************
********************** STM32F407应用开发板(高配版)*************************
**********************************************************************************
* 文件名称: DSP FFT使用                                                         *
* 文件简述：                                                                     *
* 创建日期：2020.8.21                                                          *
* 版    本：V1.0                                                                 *
* 作    者：近视未看清                                                               *
* 说    明：将时域序列信号FFT后幅值、频率、相位通过串口发送到PC                         * 
*                                     *
**********************************************************************************
*********************************************************************************/

float input[length*2];//前一个数表示实部，后一个表示虚部
float Amp[length];//存放幅值
float rate[length];//存放频率
float Phase[length];
int main(void)
{ 
	u16 i;
  NVIC_PriorityGroupConfig(NVIC_PriorityGroup_2); //设置系统中断优先级分组2
	delay_init();		  //初始化延时函数
	KEY_Init();	      //IO初始化
	uart3_init(9600);
	FFTx4_init(0,1);//FFT初始化函数（基4）
	//生成信号序列
	for(i=0;i<length;i++)
	{
		input[i*2]=100+10*arm_sin_f32(2*PI*i/length)+30*arm_sin_f32(2*PI*i*4/length)+50*arm_cos_f32(2*PI*i*8/length);	//生成输入信号实部
		input[2*i+1]=0;//时域，虚部为零
	}
	
  while(1)
	{
		key_scan(0);
		if(keydown_data==KEY0_DATA)
		{
			all_result_x4(input, Amp, rate, Phase,1024);//计算各点幅值、频率、相位
			printf("A:       F;    P\n");
			for(i=0;i<length;i++)
			printf("%.2f: %.2f: %.2f\n",Amp[i],rate[i],Phase[i]);//向串口发送幅值、频率、相位
		}
    delay_ms(5);
	}
}

由于 FFT 变换后的结果具有对称性，所以，实际上有用的数据，只有前半部分，后半部分和前半部分是对称关系，比如 第2个点和 最后一个点，5 和 1021,9和 1017 等，就是对称关系，因此我们只需要分析前半部分数据即可。这样，就只有第 1、2、5、9 这四个点幅值比较大，其他点的幅值几乎都是0，重点分析。
采样频率为 1024Hz，那么总共采集 1024 个点，频率分辨率就是 1Hz，对应到频谱上面，两个点之间的间隔就是 1Hz（跟图中第二列完全吻合）。第 1 点，即直流分量，幅值为100，第2、5、9点对应幅值分别为10，30，50，频率分别为1Hz，4Hz，8Hz。因此，可以看出跟上面生成的信号：
100+
10*sin(2*PI*i/1024)+
30*sin(2*PI*i*4/1024)+
50*cos(2*PI*i*8/1024)
几乎完全吻合。

基于STM32F4的音乐频谱分析

1.制作材料：

stm32f407核心板

16*16的可见灯点阵

蓝牙音频模块及其配套运放模块

若干杜邦线等等、

2.制作思路

​ 先将整体电路搭建好，也就是如下图类似的电路整体图。（由于博主自己搭的卖相不是很好，所以就从网上找了一个类似的光立方图）

​ 列控制采用SS8550进行频谱幅度的控制，而行控制通过SS8050进行循环控制，由于人眼能够识别的频率分界线为25HZ，因此我们可以将行循环时间定义为956us，时间循环太短了无法进行控制，io口可能会发生错误显示。（此原因可能是由于stm32f4cpu运行速度太高，io口工作频率达不到。）

3.软件设计思路

​ 博主直接移植的stm32f1的DSP库里的程序，没有用F4的dsp库（后来仔细看了下，F4与f1的dsp库差别在于，一个为整型，另外一个为浮点型，F4运行速度会快一点，但由于之前先做的频谱分析仪是用F1做的，所以就懒得改了。）

​ 其中FFT的核心代码就放在这里了，有需求的可以摘抄看看。

4.总结

​ 其实音乐频谱总体思路都挺简单的，但由于刚开始做的时候，不知道蓝牙还有音频模块，还以为需要自己写下协议啥的，导致看了几天蓝牙协议部分，当然最耗费时间的还是LED显示电路部分，采用的是杜邦线连接的，导致各种问题都出现了，接触不良，显示效果不好等等。