FFT代码说明
FFT为Fast Fourier Transformation，即快速傅里叶变换，本项目中，FFT的目标是识别频率为形如式的一个正弦信号：
       
其中，；因为单片机通过ADC接口读取该正弦信号的电压值，而12位精度的ADC的值范围在0-4096之间，如信号经过放大器后映射到0-3.3V之间，则振幅A的取值0-4096之间。
假设，信号经过放大器后，其电压值最大为3.3V，最小为0V，则此信号的振幅为1.65V，对应A=2048，即该信号为：
       
本文中给出的例程即通过FFT识别式这种正弦信号。
假设采样频率为Fs，信号频率Fn，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。
第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：
       
频率分辨率（）等于采样时间的倒数。例如要分辨0.1Hz，则需要采集10s。
1基于STM32官方DSP库的FFT算法
工程文件中包含三个函数库，分别为：
cr4_fft_64_stm32.s
cr4_fft_256_stm32.s
cr4_fft_1024_stm32.s
分别对应数据点数为64,256和1024时的FFT算法。下面将以数据点数是1024为例，说明FFT的实现过程。
1.1通过函数生成原始数据
for(i=0;i<NPT;i++)
  {
Fx=2048+2048*sin(PI2*i*21/Fs+20)+1000*sin(PI2*i*15/Fs)+100*sin(PI2*i*25/Fs)
+rand()%100;    //振幅2048，频率21HZ，为主要的谐波分量，21HZ也是所需信号的频率，rand()%100为0-100之间的随机噪声
    lBUFIN[i] = ((s16)fx)<<16;    //高位为实部，低位为虚部
  }
Fs为采样频率，此处设为102.4Hz，NPT为数据长度（即前文中提到的N），为了能够分辨0.1Hz，数据长度NPT设为1024，则频率分辨率=102.4/1024=0.1Hz。
1.2对原始数据进行FFT
STM32 库中给出的FFT函数格式如下：
Void cr4_fft_1024_stm32(lBUFOUT, lBUFIN, NPT);
其中，NPT为数据点个数，通过修改宏定义来修改NPT代表的值，此处为1024；LBUFOUT[NPT] 和LBUFIN[NPT]均是长度为NPT的32位长整型的数组，高16位为实部，低16位为虚部，LBUFOUT[NPT]用于保存FFT之后的输出值；LBUFIN[NPT]用于保存输入的数值。
对原始数据进行FFT之后，需要对输出数据进行处理。
1.2.1频率计算
每个频率点处代表的真实频率为
for(i=0;i<NPT/2;i++)
{
Fn=i*Fs/NPT        //由于此处i是从0开始的，所以不需要再减1
}
1.2.2幅值计算
经过FFT后，每个频率点处的真实幅值的计算公式如下。
       
void dsp_asm_powerMag(void)
{
  s16 lX,lY;
  u32 i;
  for(i=0;i<NPT/2;i++)      //由于FFT的频谱结果是关于奈奎斯特频率对称的，所以只计算一半的点即可
  {
    lX  = (lBUFOUT[i] << 16) >> 16;        //取低16位，虚部
    lY  = (lBUFOUT[i] >> 16);            //取高16位，实部
    {
        float X    = NPT * ((float)lX) /32768;
        float Y    = NPT * ((float)lY) /32768;
        float Mag = sqrt(X*X + Y*Y)/NPT;
        lBUFMAG[i]    = (u32)(Mag * 65536);
    } //真实振幅lBUFMAG[i]=    sqrt(lX*lX + lY*lY)*2/NPT，先除以32768，又乘以65536是为了符合浮点数的计算规律   
  }
lBUFMAG[0]= lBUFMAG[0]/2;    //直流分量不需要乘以2，所以按照乘以2计算后，要除以2
}
得到每个频率点所代表的真实频率Fn，以及该真实频率所代表的原始信号的振幅A，将其通过串口打印到电脑上。
for(i=0;i<NPT/2;i++)
  {
    printf("%4d,    %.2f    %10d\n",i,((float)i*Fs/NPT),lBUFMAG[i]);  }
结果如下图所示。
图 1 f=0hz 处幅值
单片机printf函数图 2 f=15Hz处幅值
图 3 f=21hz处幅值
图 4 f=25hz处幅值
2非官方FFT算法
另外程序中给出了非官方FFT计算方法，文件名为FFT2.c，相关的函数如下：
dsp_g2_init();  //用于生成计算所需原始数据，
dsp_g2_test();  //非STM32库中的FFT算法
printf_dsp_g2_Mag();//将各个频率点上的振幅打印到电脑上
timetest_g2FFT(); //通过定时器测试该FFT算法消耗的时间，包括计算各个频率点的振幅计算时间
非官方的FFT算法的计算结果详见”非官方FFT结果.txt”附件
3耗时测试
通过配置TIM6定时器进行计时，单位为ms。对官方FFT算法以及非官方FFT算法进行了耗时测试，测试10次FFT总时间，求取平均值，其中包括计算各个频率点幅值的时间，将测试结果打印到电脑上。
如图 5所示，官方FFT算法耗时9ms，非官方耗时411ms。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。