在數字信號處理中常常需要用到離散傅立葉變換(DFT)，以獲取信號的頻域特征。盡管傳統的DFT算法能夠獲取信號頻域特征，但是算法計算量大，耗時長， 不利于計算機實時對信號進行處理。因此至DFT被發現以來，在很長的一段時間內都不能被應用到實際的工程項目中，直到一種快速的離散傅立葉計算方法—— FFT，被發現，離散傅立葉變換才在實際的工程中得到廣泛應用。需要強調的是，FFT并不是一種新的頻域特征獲取方式，而是DFT的一種快速實現算法。本 文就FFT的原理以及具體實現過程進行詳盡講解。
DFT計算公式
本文不加推導地直接給出DFT的計算公式：
其中x(n)表示輸入的離散數字信號序列，WN為旋轉因子，X(k)為輸入序列x(n)對應 的N個離散頻率點的相對幅度。一般情況下，假設x(n)來自于低通采樣，采樣頻率為fs，那么X(k)表示了從-fs/2率開始，頻率間隔為fs/N,到 fs/2-fs/N截至的N個頻率點的相對幅度。因為DFT計算得到的一組離散頻率幅度值實際上是在頻率軸上從成周期變化的，即X(k+N)=X(k)。 因此任意取連續的N個點均可以表示DFT的計算效果，負頻率成分比較抽象，難于理解，根據X(k)的周期特性，于是我們又可以認為X(k)表示了從零頻率 開始，頻率間隔為fs/N,到fs-fs/N截至的N個頻率點的相對幅度。
N點DFT的計算量
根 據(1)式給出的DFT計算公式，我們可以知道每計算一個頻率點X(k)均需要進行N次復數乘法和N-1次復數加法，計算N各點的X(k)共需要N^2次 復數乘法和N*(N-1)次復數加法。當x(n)為實數的情況下，計算N點的DFT需要2*N^2次實數乘法，2*N*(N-1)次實數加法。
旋轉因子WN的特性
1.WN的對稱性
2.WN的周期性
3.WN的可約性
根據以上這些性質，我們可以得到式(5)的一系列有用結果
基-2 FFT算法推導
假設采樣序列點數為N=2^L，L為整數，如果不滿足這個條件可以人為地添加若干個0以使采樣序列點數滿足這一要求。首先我們將序列x(n)按照奇偶分為兩組如下：
于是根據DFT計算公式(1)有：
至此，我們將一個N點的DFT轉化為了式(7)的形式，此時k的取值為0到N-1,現在分為兩段來討論，當k為0～N/2-1的時候，因為x1(r)，x2(r)為N/2點的序列，因此，此時式(7)可以寫為：
而當 k取值為N/2~N-1時，k用k’+N/2取代，k’取值為0~N/2-1。對式(7)化簡可得：
綜合以上推導我們可以得到如下結論：一個N點的DFT變換過程可以用兩個N/2點的DFT變換過程來表示，其具體公式如式(10)所示DFT快速算法的迭代公式：
上式中X'(k’)為偶數項分支的離散傅立葉變換，X”(k’’)為奇數項分支的離散傅立葉變換。 式(10)的計算過程可以用圖1的蝶形算法流圖直觀地表示出來。
圖1 時間抽取法蝶形運算流圖
在圖1中，輸入為兩個N/2點的DFT輸出為一個N點的DFT結果，輸入輸出點數一致。運用這種表示方法，8點的DFT可以用圖2來表示：
圖2 8點DFT的4點分解
根據公式(10)，一個N點的DFT可以由兩個N/2點的DFT運算構成，再結合圖1的蝶形信號流圖可以得到圖2的8點DFT的第一次分解。該分解可以用以下幾個步驟來描述：
? 1.將N點的輸入序列按奇偶分為2組分別為N/2點的序列
? 2.分別對1中的每組序列進行DFT變換得到兩組點數為N/2的DFT變換值X1和X2
? 3.按照蝶形信號流圖將2的結果組合為一個N點的DFT變換結果
根據式(10)我們可以對圖2中的4點DFT進一步分解，得到圖3的結果，分解步驟和前面一致。
圖3 8點DFT的2點分解
最后對2點DFT進一步分解得到最終的8點FFT信號計算流圖：
圖4 8點DFT的全分解
從 圖2到圖4的過程中關于旋轉系數的變化規律需要說明一下。看起來似乎向前推一級，在奇數分組部分的旋轉系數因子增量似乎就要變大，其實不是這樣。事實上奇 數分組部分的旋轉因子指數每次增量固定為1，只是因為每向前推進一次，該分組序列的數據個數變少了，為了統一使用以原數據N為基的旋轉因子就進行了變換導 致的。每一次分組奇數部分的系數WN,這里的N均為本次分組前的序列點數。以上邊的8點DFT為例，第一次分組N=8,第二次分組N為4，為了統一根據式 (4)進行了變換將N變為了8，但指數相應的需要乘以2。
N點基-2 FFT算法的計算量
從 圖4可以看到N點DFT的FFT變換可以轉為log2(N)級級聯的蝶形運算，每一級均包含有N/2次蝶形計算。而每一個蝶形運算包含了1次復數乘法，2 次復數加法。因此N點FFT計算的總計算量為：復數乘法——N/2×log2(N)??? 復數加法——N×log2(N)。假設被采樣的序列為實數序列，那么也只有第一級的計算為實數與復數的混合計算，經過一次迭代后來的計算均變為復數計算， 在這一點上和直接的DFT計算不一致。因此對于輸入序列是復數還是實數對FFT算法的效率影響較小。一次復數乘法包含了4次實數乘法，2次實數加法，一次 復數加法包含了2次復數加法。因此對于N點的FFT計算需要總共的實數乘法數量為：2×N×log2(N)；總的復數加法次數 為：2xNxlog2(N)。
N點基-2 FFT算法的實現方法
從圖4我們可以總結出對于點數為N=2^L的DFT快速計算方法的流程：
?? 1.對于輸入數據序列進行倒位序變換。
該 變換的目的是使輸出能夠得到X(0)~X(N-1)的順序序列，同樣以8點DFT為例，該變換將順序輸入序列x(0)~x(7)變為如圖4的 x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其實現方法是：假設順序輸入序列一次村在A(0)~A(N-1) 的數組元素中，首先我們將數組下標進行二進制化(例：對于點數為8的序列只需要LOG2(8) = 3位二進制序列表示，序號6就表示為110)。二進制化以后就是將二進制序列進行倒位，倒位的過程就是將原序列從右到左書寫一次構成新的序列，例如序號為 6的二進制表示為110，倒位后變為了011，即使十進制的3。第三步就是將倒位前和倒位后的序號對應的數據互換。依然以序號6為例，其互換過程如下：
temp = A(6); A(6) = A(3); A(3) = A(6);
實際上考慮到執行效率，如果對于每一次輸入的數據都需要這個處理過程是非常浪費時間的。我們可以采用指向指針的指針來實現該過程，或者是采用指針數組來實現該過程。
2.蝶形運算的循環結構。
從 圖4中我們可以看到對于點數為N = 2^L的fft運算，可以分解為L階蝶形圖級聯，每一階蝶形圖內又分為M個蝶形組，每個蝶形組內包含K個蝶形。根據這一點我們就可以構造三階循環來實現蝶 形運算。編程過程需要注意旋轉因子與蝶形階數和蝶形分組內的蝶形個數存在關聯。
?3.浮點到定點轉換需要注意的關鍵問題
上 邊的分析都是基于浮點運算來得到的結論，事實上大多數嵌入式系統對浮點運算支持甚微，因此在嵌入式系統中進行離散傅里葉變換一般都應該采用定點方式。對于 簡單的DFT運算從浮點到定點顯得非常容易。根據式(1)，假設輸入x(n)是經過AD采樣的數字序列，AD位數為12位，則輸入信號范圍為 0~4096。為了進行定點運算我們將旋轉因子實部虛部同時擴大2^12倍，取整數部分代表旋轉因子。之后，我們可以按照(1)式計算，得到的結果與原結 果成比例關系，新的結果比原結果的2^12倍。但是，對于使用蝶形運算的fft我們不能采用這種簡單的放大旋轉因子轉為整數計算的方式。因為fft是一個 非對稱迭代過程，假設我們對旋轉因子進行了放大，根據蝶形流圖我們可以發現其最終的結果是，不同的輸入被放大了不同的倍數，對于第一個輸入x(0)永遠也 不會放大。舉一個更加形象的例子，還是以圖4為例。從圖中可以看出右側的X(0)可以直接用下式表示：
??? 4.計算過程中的溢出問題
最 后需要注意的一個問題就是計算過程中的溢出問題。在實際應用中，AD雖然有12位的位寬，但是采樣得到的信號可能較小，例如可能在0~8之間波動，也就是 說實際可能只有3位的情況。這種情況下為了在計算過程中不丟失信息，一般都需要先將輸入數據左移P位進行放大處理，數據放大可能會導致溢出，從而使計算錯 誤，而溢出的極限情況是這樣：假設我們數據位寬為D位(不包括符號位)，AD采樣位數B位，數字放大倍數P位，旋轉因此放大倍數Q位，FFT級聯運算帶來 的最大累加倍數L位。我們得到:
FFT代碼示例
根據以上的分析，我整理了一份128點的FFT代碼如下，該代碼在keil中仿真運行，未發現問題。
#define N 128 #define POWER 6//該值代表了輸入數據首先被放大的倍數，放大倍數為2^POWER #define P_COEF 8//該值代表了旋轉因子被放大的倍數，放大倍數為2^POWER #if (N == 4) #define L 2//L的定義滿足L = log2(N) #elif (N == 8) #define L 3//L的定義滿足L = log2(N) #elif (N == 16) #define L 4//L的定義滿足L = log2(N) #elif (N == 32) #define L 5//L的定義滿足L = log2(N) #elif (N == 64) #define L 6//L的定義滿足L = log2(N) #elif (N == 128) #define L 7//L的定義滿足L = log2(N) #endif //AD采樣位數為12位，本可以采用s16 x[N]定義數據空間，但是為了節省存儲空間，fft結果也將存儲在該變量空間。由于受 //fft影響變量需要重新定義，定義的方式及具體原因如下： //fft過程會乘以較大旋轉因子，因此需要32位定義 //fft過程會產生負數，因此需要有符號定義 //fft過程會出現復數，因此需要定義二維數組，[][0]存放實部，[][1]存放虛部 s32 x[N][2] = {}; //定義*p[N]是為了在第一次指針初始化以后，該數組指針按照位倒序指向數據變量x //p[i][0]代表了指向數據的實部，p[i][1]代表了指向數據的虛部 s32 *p[N]; //定義旋轉因子矩陣 //旋轉因子矩陣存儲了wn^0,wn^1,wn^2...wn^(N/2-1)這N/2個復數旋轉因子 s16 w[N>>1][2] = {256,0,256,-13,255,-25,253,-38,251,-50,248,-62,245,-74,241,-86,237,-98,231,-109,226,
????????????????? -121,220,-132,213,-142,206,-152,198,-162,190,-172,181,-181,172,-190,162,-198,152,
????????????????? -206,142,-213,132,-220,121,-226,109,-231,98,-237,86,-241,74,-245,62,-248,50,-251,38,
????????????????? -253,25,-255,13,-256,0,-256,-13,-256,-25,-255,-38,-253,-50,-251,-62,-248,-74,-245,-86,
????????????????? -241,-98,-237,-109,-231,-121,-226,-132,-220,-142,-213,-152,-206,-162,-198,-172,-190,-181,
????????????????? -181,-190,-172,-198,-162,-206,-152,-213,-142,-220,-132,-226,-121,-231,-109,-237,-98,-241,
????????????????? -86,-245,-74,-248,-62,-251,-50,-253,-38,-255,-25,-256,-13}; void init_pointer(void) { unsigned char i = 0; unsigned char j = 0; unsigned char k = 0; for(i = 0; i < N; i++) { j = 0; for(k = 0; k < L; k++) { j |= (((i >> k)&0x01)<<(L-k-1)); } p[i] = &x[j][0]; } } /* *description:基2fft主函數,該函數將借助指針數組p對全局變量數組x進行fft計算 * 并將結果存儲在數組x中 *global var:rw->x; r->p,w。(r表示讀，w表示寫，rw表示讀寫) */ void fft2(void) { u8 i;//i用于表示蝶形圖級聯的階數 u8 j;//表示蝶形分組起始點序列，蝶形分組跨度為2^i u8 k;//k表示蝶形組內偶數分支序列點號 u8 gp_distance = 1;//蝶形分組跨度 u8 temp;//temp用于記錄當前組間距離的一半,同時也表示了同一碟形兩分支間的距離 u8 gp_hf = 0;//記錄當前組的中間點序號 u8 delta = N;//旋轉因子下標增量，本來下標初始值應該為N/2，但是。。 s16 *pw = &(w[0][0]); int tp_result[2]; //用于臨時存放旋轉因子和奇數分組的乘積 //輸入信號序列放大 for(i = 0; i < N; i++) { x[i][0] <<= POWER; x[i][1] <<= POWER; } for(i = 0; i < L; i++) { temp = gp_distance; gp_distance <<= 1; for(j = 0; j < N; j+=gp_distance) { gp_hf = temp + j; pw = &(w[0][0]); for(k = j; k < gp_hf; k++)//完成一組內的所有蝶形運算 { //蝶形運算中的一組復數乘法過程 tp_result[0] = pw[0] * (p[k+temp][0]) - pw[1] * (p[k+temp][1]); tp_result[1] = pw[0] * (p[k+temp][1]) + pw[1] * (p[k+temp][0]); tp_result[0] >>= P_COEF; tp_result[1] >>= P_COEF; //蝶形運算中的2組復數加法過程 p[k+temp][0] = p[k][0] - tp_result[0]; p[k+temp][1] = p[k][1] - tp_result[1]; p[k][0] += tp_result[0]; p[k][1] += tp_result[1]; pw += delta; } } delta >>= 1; } }