STM32F103 如何實現(xiàn) FFT?

筆者能力有限，如果文中出現(xiàn)錯誤的地方，還希望各位朋友能夠給我指出來，我將不勝感激，謝謝~

引言

數(shù)字信號在我們生活中隨處可見，自然而然地就會涉及到對于數(shù)字信號的處理，最為典型的一個應(yīng)用就是示波器，在使用示波器的過程當(dāng)中，我們會通過示波器測量到信號的頻率以及幅值，同時我們也可以通過示波器對測量到的信號進(jìn)行 FFT ，從而能夠觀察到待測信號的頻譜，方便直觀的看出信號的高頻分量和低頻分量，從而幫助我們?nèi)コ盘栔袛y帶的噪聲。而在嵌入式方面的應(yīng)用，我們可以直接使用 DSP 芯片對信號進(jìn)行處理，同時， ARM 公司推出的 Cortex-M4F 內(nèi)核是帶有 FPU ,DSP 和 SIMD 單元的，針對于這些單元也增加了專用的指令，指令如下圖所示：

ARM 官方也對此專門做了一個 DSP 方面的庫，方便用戶調(diào)用。關(guān)于 Cortex M4 的信號處理本文暫不進(jìn)行闡述，相反本文的對象是 Cortex M3 ,基于 STM32F103 的 FFT，而在上述圖中，我們看到針對于 Cortex M3 來說，是不帶 FPU 以及 DSP 的，那有如何來進(jìn)行 FFT 呢？

FFT 的提出

在數(shù)字信號處理中常常需要使用到離散傅里葉變換(DFT)，從而能夠獲取到信號的頻域特征。盡管傳統(tǒng)的 DFT 算法能夠獲取到信號的頻域特征，但是算法計算量大，耗時長，不利于進(jìn)行計算機(jī)實時對信號進(jìn)行處理。因此才有了 FFT 的出現(xiàn)。需要強(qiáng)調(diào)的是，F(xiàn)FT 并不是一種新的頻域特征獲取方式，而是 DFT 的一種快速實現(xiàn)算法。FFT 之所以能夠改善運算量，是因為其充分運用了 DFT 運算中的對稱性和周期性，從而能夠?qū)⑦\算量從 N^2 減少到 N*log2(N),其中 N 為待計算的序列的長度。當(dāng) N 非常大的時候，這種優(yōu)化在時間維度上的提升是非常顯著的。下面是關(guān)于 DFT 和 FFT 所需乘法次數(shù)的比較曲線。

FFT 變換之后和原始信號的對應(yīng)關(guān)系

假設(shè)我們對一個波形進(jìn)行了采樣，采樣了 N 個點，經(jīng)過 FFT 之后，就可以得到 N 個點的 FFT 結(jié)果，每一個點就對應(yīng)著一個頻率點。這個點的模值，就是該頻率下的幅度特性。具體的關(guān)系就是如果原始信號的峰值為 A ，那么 FFT 的結(jié)果的每個點的模值就是 A 的 N / 2 倍。而第一個點就是直流分量，它的模值是直流分量的 N 倍。而每個點的相位就是在該頻率下的信號的相位，第一個點表示的是直流分量，也就是 0 HZ的點，而最后一個點 N 的再下一個點（實際這個點是不存在的），也就是 N+1 個點則表示的是采樣頻率 Fs，這中間被 N - 1 個點平均分成 N 等份，每一個點的頻率依次增加。也就是如果要計算某個點的頻率，那么就只需要這樣計算即可：Fn = (n - 1) * Fs / N。
從上述所展示的公式，我們可以知道 Fn 所能夠分辨的頻率為 Fs / N，如果采樣頻率 Fs 為 1024Hz,采樣點數(shù)為 1024 點，則可以分辨到 1 HZ。也就是說采樣 1s 時間的信號并做 FFT ，則結(jié)果可以分析精確到 1 Hz，如果采樣 2 s 時間的信號并做 FFT，則結(jié)果可以分析精確到 0.5 Hz，所以也就說明了一個道理，如果要提高頻率分辨率，則必須增加采樣點數(shù)，也就是采樣時間，下面這張圖更能夠清晰地表示這種關(guān)系：

將原信號變換之后的頻譜的寬度與原始信號也存在一定的關(guān)系。根據(jù) Nyquist采樣定律，F(xiàn)FT 之后的頻譜寬度最大只能是原始信號采樣率的 1/2，如果原始信號的采樣頻率為 4GS/s，那么 FFT 之后的頻寬最多只能是 2GHz，這還只是理想情況。所以也能夠得出一個結(jié)論：時域信號的采樣率乘上一個固定系數(shù)即是變換之后頻譜的寬度，可以用如下所示的一張圖清晰說明：

經(jīng)過上述的分析，我們有了如下的結(jié)論：更高的頻譜分辨率需要有更長的采樣時間，更寬的頻譜分布需要提高對于原始信號的采樣率，那我們在實際的使用過程中，當(dāng)然是希望頻譜更寬，分辨率更加精確，那么示波器的長存儲就是必要的。

F103 如何進(jìn)行 FFT

FFT 匯編庫介紹

在本文的開頭敘述了 ARM Cortex M4 具有 FPU 以及 DSP 指令，同時 ARM 官方也出了 DSP 方面的庫來進(jìn)行數(shù)字信號處理方面的工作，那么針對于 ARM Cortex M3 的 STM32F103 又是如何進(jìn)行 FFT 的呢，顯然，如果我們用 C 語言直接編寫 FFT 算法，那樣子的效率是極其低下的，因此，本文采用的方法是 ST 官方匯編 FFT 庫的應(yīng)用，由于官網(wǎng)現(xiàn)在找不到這個軟件包，可以在公眾號后臺回復(fù) FFT 獲取軟件包。
簡單介紹一下，這個庫是由匯編實現(xiàn)的，而且是基 4 算法，所以實現(xiàn) FFT 在速度上較快。如果 X[N]是采樣信號的話，使用 FFT 時必須滿足如下兩條：

• N 得滿足 4^n (n = 1,2,3…),也就是以 4 為基數(shù)。

• 采樣信號必須是 32 位數(shù)據(jù)，高 16 位存實部，低 16 位存虛部(這個是針對大端模式)，小端模式是高位存虛部，低位存實部，一般常用的是小端模式。

匯編 FFT 的實現(xiàn)主要包括以下三個函數(shù)：

1. cr4_fft_64_stm32 : 實現(xiàn) 64 點 FFT

2. cr4_fft_256_stm32: 實現(xiàn) 256 點 FFT

3. cr4_fft_1024_stm32: 實現(xiàn) 1024 點 FFT

FFT 匯編庫移植

將我們下載到的文件進(jìn)行解壓得到如下所示的文件：

進(jìn)一步的我們需要將文件加入到我們的 keil 工程，加入工程之后的圖如下所示：

因為本文是針對于 256 點的 FFT ，因此只需要將cr4_fft_256_stm32 添加進(jìn)來即可，加進(jìn)來之后，再使用到 FFT 的文件里添加相關(guān)路徑就可以。下面講述具體的代碼實例。

代碼實例

FFT 計算幅值

首先我們定義采樣的點數(shù)，以及 FFT 的輸入數(shù)組，輸出數(shù)組，以及各個諧波的幅值：

#define  NPT  256             /* 采樣點數(shù) */
uint32_t lBufInArray[NPT];    /* FFT 運算的輸入數(shù)組 */
uint32_t lBufOutArray[NPT/2]; /* FFT 運算的輸出數(shù)組 */
uint32_t lBufMagArray[NPT/2]; /* 各次諧波的幅值 */

在上述中，F(xiàn)FT 的輸出數(shù)組和各次諧波的幅值的數(shù)組中只有采樣點數(shù)的一半，是因為 FFT 計算出來的數(shù)據(jù)是對稱的，因此通常而言只取一半的數(shù)據(jù)。
下面是波形采樣并進(jìn)行 FFT 的代碼：

void adc_sample(void)
{
    for (i = 0; i < NPT; i++)
    {
        lBufInArray[i] = ADC_ConvertedValue[i];
    }
    cr4_fft_256_stm32(lBufOutArray,lBufInArray, NPT);
    GetPowerMag();
}

for循環(huán)里將 ADC 采集的數(shù)據(jù)存儲到 FFT 運算的輸入數(shù)組中去，在這里需要注意的是 STM32 是小端模式，因此采樣數(shù)據(jù)是高位存虛部，低位存實部。緊接著就是調(diào)用匯編函數(shù)進(jìn)行 FFT，F(xiàn)FT 運算之后，就進(jìn)行幅值的計算，幅值的計算函數(shù)如下所示：

void GetPowerMag(void)
{
    signed short lX,lY;
    float X,Y,Mag;
    unsigned short i;
    for(i=0; i<NPT/2; i++)
    {
        lX  = (miniscope.freqency.lBufOutArray[i] << 16) >> 16;
        lY  = (miniscope.freqency.lBufOutArray[i] >> 16);

        //除以32768再乘65536是為了符合浮點數(shù)計算規(guī)律
        X = NPT * ((float)lX) / 32768;
        Y = NPT * ((float)lY) / 32768;
        Mag = sqrt(X * X + Y * Y)*1.0/ NPT;
        if(i == 0)    
            miniscope.freqency.lBufMagArray[i] = (unsigned long)(Mag * 32768);
        else
            miniscope.freqency.lBufMagArray[i] = (unsigned long)(Mag * 65536);
    }
}

上述代碼中，lx 和 ly 的計算中，分別取的是FFT 的輸出數(shù)組的高位和低位。進(jìn)一步的，在計算 x 和 y 的時，除以 32768 是為了符合浮點數(shù)計算規(guī)律，至于為什么要進(jìn)行浮點化，是因為浮點化就好像 10 進(jìn)制里面的科學(xué)計數(shù)法。32768 = 2 的 15 次。除以 32768 也就是去除了浮點數(shù)后面的那個基數(shù)，只剩下前面的。比如 1991 改寫成 1.991 * 10 的三次冪，除以 10 的三次方，只剩下 1.991，方便于下面的運算。而在后面又乘以 32768 和 65536 是因為要恢復(fù)到原先數(shù)據(jù)的大小，為什么下標(biāo)為 0 的乘以 32768，而大于 0 的乘以 65535，是因為下標(biāo)為 0 的代表的是直流分量，而剩余的是求出的乘以 2 才是實際模值。

FFT 計算頻率

在本文的前面，筆者給出了這樣一個公式用來計算 FFT 變換之后每個點對應(yīng)的頻率：
Fn = (n - 1) * Fs / N
N 是采樣的點數(shù)，F(xiàn)s 是采樣頻率，采樣點數(shù)已經(jīng)知道，還剩下采樣頻率未知，采樣頻率說白了也就是采樣一個點的時間，也是 1 s 鐘采樣的點數(shù)，而這個該怎么確定呢?，F(xiàn)有兩種方法，第一種方法是在單片機(jī)進(jìn)行 ADC 采集時，通過延時的方法每隔一段時間進(jìn)行讀取轉(zhuǎn)換得到的數(shù)據(jù)，而這個延時的時間就是采樣頻率，這樣聽起來略顯粗糙。另一種比較精確的方法，是通過 DMA + TIM 的方法，也就是通過 TIM 產(chǎn)生 PWM ，通過 PWM 觸發(fā) ADC 進(jìn)行采集，這個時候，PWM 的頻率也就是 ADC 的采樣率，只需要控制 PWM 的頻率就可以控制 ADC 的采樣率，采集的數(shù)據(jù)通過 DMA 搬運至內(nèi)存，當(dāng)采樣的點數(shù)達(dá)到規(guī)定的采樣點數(shù)時，觸發(fā) DMA 中斷，在中斷里給出數(shù)據(jù)處理的信號，進(jìn)一步進(jìn)行 FFT，具體的原理及代碼參考筆者的這篇文章：STM32 定時器觸發(fā) ADC 多通道采集，DMA搬運至內(nèi)存。下面是一個簡單的代碼示例：

void wave_frequency_calculate(void)
{
    sample_frequency = 72000000.0 / (float)((sample_arr + 1) * (sample_psc + 1));
    wave_frequency = frequency_max_position * sample_frequency / NPT;
}

上述第一行代碼是根據(jù)公式計算 pwm 的頻率的公式，而 pwm 的頻率也就是我們所需要的的采樣頻率。第二條代碼中的 frequency_max_position 是除了直流分量幅值最大的點在數(shù)組中的位置，而這個點所對應(yīng)的頻率也就是我們采樣波形的頻率，至此，我們就計算出了采樣波形的頻率。

結(jié)論

上述就是關(guān)于 STM32F103 中實現(xiàn) FFT 的一個基本方法，通過 FFT 計算出了波形的頻譜，能夠在不借助 DSP 芯片的前提下比較快的實現(xiàn)了 FFT ，對我們在 F103 平臺上進(jìn)行信號處理提供了很大的幫助，這就是本次分享的內(nèi)容啦。

如果覺得我的文章對你所有幫助，歡迎轉(zhuǎn)發(fā)，點贊，再看三連呀~