基于STM32F407实现离散傅里叶变换（FFT、DFT），计算指定频率的幅值

前言：本人的课题是关于EIT采集系统设计，所谓的EIT，简单的说就是往人体注入特定频率的电流信号，通过采集反馈的电压信号，进而使用成像算法重构人体内部的阻抗分布。由于采集到的电压包含其它频率的热噪声，为了只保留注入频率的信号成分，需要对采集到的电压信号进行FFT处理。

在本文应用中，FFT相当于一个带通滤波器，用于获取指定频率的信号信息。关于快速傅里叶变化这里不做过多的介绍，具体可参考别人写的博客：

如何 FFT(快速傅里叶变换) 求幅度、频率（超详细 含推导过程）_Xav Zewen的博客

本文主要介绍如何在STM32F407上实现对特定频率进行FFT。重新更新并完善了一下，将代码整理为函数，方便读者调用。

一、使用DSP库进行FFT计算

1.1 DSP库开启

我们知道，相比于整形运算，浮点运算会大量消耗算力，若直接让STM32强行计算，难以满足实时性的需求。好在STM32F407是具有浮点运算(FPU)功能，可以通过MDK配置：target->Roating Point Hardware->Use Single Precison中打开。

在进行FFT计算时，我们还需要用到三角计算，因此我们还需要添加dsp数学库，调用库中的函数进行数学运算。DSP库的开启如下所示。

同时在MDK配置中添加头文件：ARM_MATH_CM4

 通过上述操作，我们便可进入编程环节。

1.2 调用DSP库进行FFT计算

那么，如何通过FFT变换，获取指定频率的幅值信息呢？下面举个例子：

目的：获取电压数据中10kHz的幅值分量

要求：待计算的频率、ADC的采样频率、采集样本量、数据点N 应该满足以下等式：

在STM32程序中，我们可以通过中断设定ADC的采样频率，进而配平上述等式。如定义一个1us定时器中断，在中断任务中执行一次ADC采集（此时采样频率为1MHz），一共采集1000次。

此时FFT的参数为：

        待计算的频率——10kHz

        ADC的采样频率——1MHz

        采集样本量——1000

配平上述式子，求得数据点 N = 10。

【注】 1. 若等式配不平，会导致频谱泄露，造成数据失真；

            2. 采样依旧要满足采样定理，即：ADC的采样频率 >2*电压数据中最大的频率分量；

            3. 由于FFT变换的特点，采集样本量为2的指数倍能提高计算速度。 

使用DSP库进行FFT变换的函数如下：

// FFT计算函数
// *DATA: 导入待FFT计算的原始数组指针
// num：采集样本量
// N：需要保存的第几个数据点
float FFT_Calculation(float *DATA, int num, int N)
{
    float array_FFT_output[num];        //储存FFT变换后的512个数据
    float array_arm_cmplx_mag[num];     //储存FFT变换后的512个数据的幅值信息

    arm_rfft_fast_instance_f32 S;
    arm_rfft_fast_init_f32(&S, fftSize);        //初始化结构体S中的参数
    arm_rfft_fast_f32(&S, array_f32, array_FFT_output, 0);          //fft正变换 
    arm_cmplx_mag_f32(array_FFT_output, array_arm_cmplx_mag, num);  //计算幅值  

    return array_arm_cmplx_mag[N];  
}    

下面简单的示范一下这个函数怎么使用：

float Data_FFT[1000];        // 待FFT计算的原始数据(读者自行赋值)
float result;

result = FFT_Calculation(Data_FFT,1000,10);    // 1000个采样点数，需要保存第10个频点

二、FFT算法的优化：使用DFT计算单一频点信息

虽然FFT计算十分方便，但是，当我们只需要单一频点数据时，FFT由于计算了大量的无效数据，消耗了算力。在上面的例子中，我们只需要获得10kHz的电压幅值，这时代码可以优化为离散型傅里叶变化（DFT）。离散型傅里叶变化的计算公式如下：

离散型傅里叶变换（DFT）的计算代码如下：

#include "arm_math.h"   // 代码中涉及了sin cos 计算, 需按1.1小节开启DSP库

// *Data: 导入待DFT计算的原始数组指针
// num：采集样本量
// N：需要保存的第几个数据点
float DFT_Calculation(float *Data, int num, int N)
{
        unsigned int i;
        float SUM_Re = 0;        //实频数值
        float SUM_Im = 0;        //虚频数值
        float result = 0;            // 计算结果
        
        //FFT展开式
        for(i=0;i<num;i++)
        {
            SUM_Re = SUM_Re + Data[i]*cos(2*3.1415926*N*i/num);
            SUM_Im = SUM_Im - Data[i]*sin(2*3.1415926*N*i/num);
        }

        //计算幅值
        result = sqrt(SUM_Re*SUM_Re + SUM_Im*SUM_Im);

        return result;
}

该函数的使用方法同FFT一致：

float Data_DFT[1000];        // 待DFT计算的原始数据(读者自行赋值)
float result;

result = DFT_Calculation(Data_DFT,1000,10);    // 1000个采样点数，需要保存第10个频点

经过测试，计算单一频点信息时，使用DFT算法相比于FFT，时间节省了约51%。

进一步提升

如果代码对于实时性有要求，在内存和算力的寸土寸金的单片机上，可以通过查表法，代替耗时的三角函数计算。

如上面的代码， cos(2*3.1415926*N*i/num) 和 sin(2*3.1415926*N*i/num) 的计算结果是固定值，可以提前计算出N=1000，k=10，i取0~999时的cos和sin值，在DFT计算是查询预先计算好的三角函数值，以节省算力。