基于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变换，获取指定频率的幅值信息呢？下面举个例子：

使用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计算是查询预先计算好的三角函数值，以节省算力。