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

引言

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

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

FFT 的提出

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

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

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

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

經(jīng)過(guò)上述的分析，我們有了如下的結(jié)論：更高的頻譜分辨率需要有更長(zhǎng)的采樣時(shí)間，更寬的頻譜分布需要提高對(duì)于原始信號(hào)的采樣率，那我們?cè)趯?shí)際的使用過(guò)程中，當(dāng)然是希望頻譜更寬，分辨率更加精確，那么示波器的長(zhǎng)存儲(chǔ)就是必要的。

F103 如何進(jìn)行 FFT

FFT 匯編庫(kù)介紹

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

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

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

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

1. cr4_fft_64_stm32 : 實(shí)現(xiàn) 64 點(diǎn) FFT

2. cr4_fft_256_stm32: 實(shí)現(xiàn) 256 點(diǎn) FFT

3. cr4_fft_1024_stm32: 實(shí)現(xiàn) 1024 點(diǎn) FFT

FFT 匯編庫(kù)移植

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

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

因?yàn)楸疚氖轻槍?duì)于 256 點(diǎn)的 FFT ，因此只需要將cr4_fft_256_stm32 添加進(jìn)來(lái)即可，加進(jìn)來(lái)之后，再使用到 FFT 的文件里添加相關(guān)路徑就可以。下面講述具體的代碼實(shí)例。

代碼實(shí)例

FFT 計(jì)算幅值

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

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

在上述中，F(xiàn)FT 的輸出數(shù)組和各次諧波的幅值的數(shù)組中只有采樣點(diǎn)數(shù)的一半，是因?yàn)?FFT 計(jì)算出來(lái)的數(shù)據(jù)是對(duì)稱(chēng)的，因此通常而言只取一半的數(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ù)存儲(chǔ)到 FFT 運(yùn)算的輸入數(shù)組中去，在這里需要注意的是 STM32 是小端模式，因此采樣數(shù)據(jù)是高位存虛部，低位存實(shí)部。緊接著就是調(diào)用匯編函數(shù)進(jìn)行 FFT，F(xiàn)FT 運(yùn)算之后，就進(jìn)行幅值的計(jì)算，幅值的計(jì)算函數(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是為了符合浮點(diǎn)數(shù)計(jì)算規(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 的計(jì)算中，分別取的是FFT 的輸出數(shù)組的高位和低位。進(jìn)一步的，在計(jì)算 x 和 y 的時(shí)，除以 32768 是為了符合浮點(diǎn)數(shù)計(jì)算規(guī)律，至于為什么要進(jìn)行浮點(diǎn)化，是因?yàn)楦↑c(diǎn)化就好像 10 進(jìn)制里面的科學(xué)計(jì)數(shù)法。32768 = 2 的 15 次。除以 32768 也就是去除了浮點(diǎn)數(shù)后面的那個(gè)基數(shù)，只剩下前面的。比如 1991 改寫(xiě)成 1.991 * 10 的三次冪，除以 10 的三次方，只剩下 1.991，方便于下面的運(yùn)算。而在后面又乘以 32768 和 65536 是因?yàn)橐謴?fù)到原先數(shù)據(jù)的大小，為什么下標(biāo)為 0 的乘以 32768，而大于 0 的乘以 65535，是因?yàn)橄聵?biāo)為 0 的代表的是直流分量，而剩余的是求出的乘以 2 才是實(shí)際模值。

FFT 計(jì)算頻率

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

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ù)公式計(jì)算 pwm 的頻率的公式，而 pwm 的頻率也就是我們所需要的的采樣頻率。第二條代碼中的 frequency_max_position 是除了直流分量幅值最大的點(diǎn)在數(shù)組中的位置，而這個(gè)點(diǎn)所對(duì)應(yīng)的頻率也就是我們采樣波形的頻率，至此，我們就計(jì)算出了采樣波形的頻率。

結(jié)論

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

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