前言：本人的課題是關(guān)于EIT采集系統(tǒng)設(shè)計(jì)，所謂的EIT，簡(jiǎn)單的說就是往人體注入特定頻率的電流信號(hào)，通過采集反饋的電壓信號(hào)，進(jìn)而使用成像算法重構(gòu)人體內(nèi)部的阻抗分布。由于采集到的電壓包含其它頻率的熱噪聲，為了只保留注入頻率的信號(hào)成分，需要對(duì)采集到的電壓信號(hào)進(jìn)行FFT處理。
在本文應(yīng)用中，F(xiàn)FT相當(dāng)于一個(gè)帶通濾波器，用于獲取指定頻率的信號(hào)信息。關(guān)于快速傅里葉變化這里不做過多的介紹，具體可參考別人寫的博客：
如何 FFT(快速傅里葉變換) 求幅度、頻率（超詳細(xì) 含推導(dǎo)過程）_Xav Zewen的博客
本文主要介紹如何在STM32F407上實(shí)現(xiàn)對(duì)特定頻率進(jìn)行FFT。重新更新并完善了一下，將代碼整理為函數(shù)，方便讀者調(diào)用。
一、使用DSP庫進(jìn)行FFT計(jì)算
1.1 DSP庫開啟
我們知道，相比于整形運(yùn)算，浮點(diǎn)運(yùn)算會(huì)大量消耗算力，若直接讓STM32強(qiáng)行計(jì)算，難以滿足實(shí)時(shí)性的需求。好在STM32F407是具有浮點(diǎn)運(yùn)算(FPU)功能，可以通過MDK配置：target->Roating Point Hardware->Use Single Precison中打開。
在進(jìn)行FFT計(jì)算時(shí)，我們還需要用到三角計(jì)算，因此我們還需要添加dsp數(shù)學(xué)庫，調(diào)用庫中的函數(shù)進(jìn)行數(shù)學(xué)運(yùn)算。DSP庫的開啟如下所示。
同時(shí)在MDK配置中添加頭文件：ARM_MATH_CM4
?通過上述操作，我們便可進(jìn)入編程環(huán)節(jié)。
1.2 調(diào)用DSP庫進(jìn)行FFT計(jì)算
那么，如何通過FFT變換，獲取指定頻率的幅值信息呢？下面舉個(gè)例子：
目的：獲取電壓數(shù)據(jù)中10kHz的幅值分量
要求：待計(jì)算的頻率、ADC的采樣頻率、采集樣本量、數(shù)據(jù)點(diǎn)N 應(yīng)該滿足以下等式：
在STM32程序中，我們可以通過中斷設(shè)定ADC的采樣頻率，進(jìn)而配平上述等式。如定義一個(gè)1us定時(shí)器中斷，在中斷任務(wù)中執(zhí)行一次ADC采集（此時(shí)采樣頻率為1MHz），一共采集1000次。
此時(shí)FFT的參數(shù)為：
????????待計(jì)算的頻率——10kHz
????????ADC的采樣頻率——1MHz
????????采集樣本量——1000
配平上述式子，求得數(shù)據(jù)點(diǎn) N = 10。
【注】 1. 若等式配不平，會(huì)導(dǎo)致頻譜泄露，造成數(shù)據(jù)失真；
? ? ? ? ? ? 2. 采樣依舊要滿足采樣定理，即：ADC的采樣頻率 >2*電壓數(shù)據(jù)中最大的頻率分量；
? ? ? ? ? ? 3. 由于FFT變換的特點(diǎn)，采集樣本量為2的指數(shù)倍能提高計(jì)算速度。?
使用DSP庫進(jìn)行FFT變換的函數(shù)如下：
// FFT計(jì)算函數(shù)
// *DATA: 導(dǎo)入待FFT計(jì)算的原始數(shù)組指針
// num：采集樣本量
// N：需要保存的第幾個(gè)數(shù)據(jù)點(diǎn)
float FFT_Calculation(float *DATA, int num, int N)
{
float array_FFT_output[num]; //儲(chǔ)存FFT變換后的512個(gè)數(shù)據(jù)
float array_arm_cmplx_mag[num]; //儲(chǔ)存FFT變換后的512個(gè)數(shù)據(jù)的幅值信息
arm_rfft_fast_instance_f32 S;
arm_rfft_fast_init_f32(&S, fftSize); //初始化結(jié)構(gòu)體S中的參數(shù)
arm_rfft_fast_f32(&S, array_f32, array_FFT_output, 0); //fft正變換
arm_cmplx_mag_f32(array_FFT_output, array_arm_cmplx_mag, num); //計(jì)算幅值
? ? return array_arm_cmplx_mag[N];
} ? 下面簡(jiǎn)單的示范一下這個(gè)函數(shù)怎么使用：
float Data_FFT[1000];? ? ? ? // 待FFT計(jì)算的原始數(shù)據(jù)(讀者自行賦值)
float result;
result = FFT_Calculation(Data_FFT,1000,10);? ? // 1000個(gè)采樣點(diǎn)數(shù)，需要保存第10個(gè)頻點(diǎn)二、FFT算法的優(yōu)化：使用DFT計(jì)算單一頻點(diǎn)信息
雖然FFT計(jì)算十分方便，但是，當(dāng)我們只需要單一頻點(diǎn)數(shù)據(jù)時(shí)，F(xiàn)FT由于計(jì)算了大量的無效數(shù)據(jù)，消耗了算力。在上面的例子中，我們只需要獲得10kHz的電壓幅值，這時(shí)代碼可以優(yōu)化為離散型傅里葉變化（DFT）。離散型傅里葉變化的計(jì)算公式如下：
離散型傅里葉變換（DFT）的計(jì)算代碼如下：
#include "arm_math.h" // 代碼中涉及了sin cos 計(jì)算, 需按1.1小節(jié)開啟DSP庫
// *Data: 導(dǎo)入待DFT計(jì)算的原始數(shù)組指針
// num：采集樣本量
// N：需要保存的第幾個(gè)數(shù)據(jù)點(diǎn)
float DFT_Calculation(float *Data, int num, int N)
{
unsigned int i;
float SUM_Re = 0; //實(shí)頻數(shù)值
float SUM_Im = 0; //虛頻數(shù)值
float result = 0; // 計(jì)算結(jié)果
//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);
}
//計(jì)算幅值
result = sqrt(SUM_Re*SUM_Re + SUM_Im*SUM_Im);
return result;
}該函數(shù)的使用方法同F(xiàn)FT一致：
float Data_DFT[1000];? ? ? ? // 待DFT計(jì)算的原始數(shù)據(jù)(讀者自行賦值)
float result;
result = DFT_Calculation(Data_DFT,1000,10);? ? // 1000個(gè)采樣點(diǎn)數(shù)，需要保存第10個(gè)頻點(diǎn)經(jīng)過測(cè)試，計(jì)算單一頻點(diǎn)信息時(shí)，使用DFT算法相比于FFT，時(shí)間節(jié)省了約51%。
進(jìn)一步提升
如果代碼對(duì)于實(shí)時(shí)性有要求，在內(nèi)存和算力的寸土寸金的單片機(jī)上，可以通過查表法，代替耗時(shí)的三角函數(shù)計(jì)算。文章來源：http://www.zghlxwxcb.cn/news/detail-521672.html
如上面的代碼， cos(2*3.1415926*N*i/num) 和 sin(2*3.1415926*N*i/num) 的計(jì)算結(jié)果是固定值，可以提前計(jì)算出N=1000，k=10，i取0~999時(shí)的cos和sin值，在DFT計(jì)算是查詢預(yù)先計(jì)算好的三角函數(shù)值，以節(jié)省算力。文章來源地址http://www.zghlxwxcb.cn/news/detail-521672.html
到了這里，關(guān)于基于STM32F407實(shí)現(xiàn)快速傅里葉變化（FFT），計(jì)算指定頻率的幅值的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
