在 數(shù)字信號(hào)處理中常常需要用到離散傅立葉變換(DFT)，以獲取信號(hào)的頻域特征。盡管傳統(tǒng)的DFT算法能夠獲取信號(hào)頻域特征，但是算法計(jì)算量大，耗時(shí)長，不 利于計(jì)算機(jī)實(shí)時(shí)對信號(hào)進(jìn)行處理。因此至DFT被發(fā)現(xiàn)以來，在很長的一段時(shí)間內(nèi)都不能被應(yīng)用到實(shí)際的工程項(xiàng)目中，直到一種快速的離散傅立葉計(jì)算方法—— FFT，被發(fā)現(xiàn)，離散是傅立葉變換才在實(shí)際的工程中得到廣泛應(yīng)用。需要強(qiáng)調(diào)的是，F(xiàn)FT并不是一種新的頻域特征獲取方式，而是DFT的一種快速實(shí)現(xiàn)算法。 本文就FFT的原理以及具體實(shí)現(xiàn)過程進(jìn)行詳盡講解。

DFT計(jì)算公式

本文不加推導(dǎo)地直接給出DFT的計(jì)算公式：

其 中x(n)表示輸入的離散數(shù)字信號(hào)序列，WN為旋轉(zhuǎn)因子，X(k)一組N點(diǎn)組成的頻率成分的相對幅度。一般情況下，假設(shè)x(n)來自于低通采樣，采樣頻率 為fs，那么X(k)表示了從-fs/2率開始，頻率間隔為fs/N,到fs/2-fs/N截至的N個(gè)頻率點(diǎn)的相對幅度。因?yàn)镈FT計(jì)算得到的一組離散頻 率幅度只實(shí)際上是在頻率軸上從成周期變化的，即X(k+N)=X(k)。因此任意取N個(gè)點(diǎn)均可以表示DFT的計(jì)算效果，負(fù)頻率成分比較抽象，難于理解，根 據(jù)X(k)的周期特性，于是我們又可以認(rèn)為X(k)表示了從零頻率開始，頻率間隔為fs/N,到fs-fs/N截至的N個(gè)頻率點(diǎn)的相對幅度。

N點(diǎn)DFT的計(jì)算量

根 據(jù)(1)式給出的DFT計(jì)算公式，我們可以知道每計(jì)算一個(gè)頻率點(diǎn)X(k)均需要進(jìn)行N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法，計(jì)算N各點(diǎn)的X(k)共需要N^2次 復(fù)數(shù)乘法和N*(N-1)次復(fù)數(shù)加法。當(dāng)x(n)為實(shí)數(shù)的情況下，計(jì)算N點(diǎn)的DFT需要2*N^2次實(shí)數(shù)乘法，2*N*(N-1)次實(shí)數(shù)加法。

旋轉(zhuǎn)因子WN的特性

  1.WN的對稱性

        2.WN的周期性

  3.WN的可約性

根據(jù)以上這些性質(zhì)，我們可以得到式(5)的一些列有用結(jié)果

基-2 FFT算法推導(dǎo)

假設(shè)采樣序列點(diǎn)數(shù)為N=2^L，L為整數(shù)，如果不滿足這個(gè)條件可以人為地添加若干個(gè)0以使采樣序列點(diǎn)數(shù)滿足這一要求。首先我們將序列x(n)按照奇偶分為兩組如下：

于是根據(jù)DFT計(jì)算公式(1)有：

至此，我們將一個(gè)N點(diǎn)的DFT轉(zhuǎn)化為了式(7)的形式，此時(shí)k的取值為0到N-1,現(xiàn)在分為兩段來討論，當(dāng)k為0～N/2-1的時(shí)候，因?yàn)閤1(r)，x2(r)為N/2點(diǎn)的序列，因此，此時(shí)式(7)可以寫為：

而當(dāng) k取值為N/2~N-1時(shí)，k用k’+N/2取代，k’取值為0~N/2-1。對式(7)化簡可得：

綜合以上推導(dǎo)我們可以得到如下結(jié)論：一個(gè)N點(diǎn)的DFT變換過程可以用兩個(gè)N/2點(diǎn)的DFT變換過程來表示，其具體公式如式(10)所示DFT快速算法的迭代公式：

上式中X'(k’)為偶數(shù)項(xiàng)分支的離散傅立葉變換，X''(k’’)為奇數(shù)項(xiàng)分支的離散傅立葉變換。 式(10)的計(jì)算過程可以用圖1的蝶形算法流圖直觀地表示出來。

圖1 時(shí)間抽取法蝶形運(yùn)算流圖

在圖1中，輸入為兩個(gè)N/2點(diǎn)的DFT輸出為一個(gè)N點(diǎn)的DFT結(jié)果，輸入輸出點(diǎn)數(shù)一致。運(yùn)用這種表示方法，8點(diǎn)的DFT可以用圖2來表示：

圖2 8點(diǎn)DFT的4點(diǎn)分解

根據(jù)公式(10)，一個(gè)N點(diǎn)的DFT可以由兩個(gè)N/2點(diǎn)的DFT運(yùn)算構(gòu)成，再結(jié)合圖1的蝶形信號(hào)流圖可以得到圖2的8點(diǎn)DFT的第一次分解。該分解可以用以下幾個(gè)步驟來描述：

  1.將N點(diǎn)的輸入序列按奇偶分為2組分別為N/2點(diǎn)的序列

  2.分別對1中的每組序列進(jìn)行DFT變換得到兩組點(diǎn)數(shù)為N/2的DFT變換值X1和X2

  3.按照蝶形信號(hào)流圖將2的結(jié)果組合為一個(gè)N點(diǎn)的DFT變換結(jié)果

根據(jù)式(10)我們可以對圖2中的4點(diǎn)DFT進(jìn)一步分解，得到圖3的結(jié)果，分解步驟和前面一致。

圖3 8點(diǎn)DFT的2點(diǎn)分解

最后對2點(diǎn)DFT進(jìn)一步分解得到最終的8點(diǎn)FFT信號(hào)計(jì)算流圖：

圖4 8點(diǎn)DFT的全分解

從 圖2到圖4的過程中關(guān)于旋轉(zhuǎn)系數(shù)的變化規(guī)律需要說明一下。看起來似乎向前推一級(jí)，在奇數(shù)分組部分的旋轉(zhuǎn)系數(shù)因子增量似乎就要變大，其實(shí)不是這樣。事實(shí)上奇 數(shù)分組部分的旋轉(zhuǎn)因子指數(shù)每次增量固定為1，只是因?yàn)槊肯蚯巴七M(jìn)一次，該分組序列的數(shù)據(jù)個(gè)數(shù)變少了，為了統(tǒng)一使用以原數(shù)據(jù)N為基的旋轉(zhuǎn)因子就進(jìn)行了變換導(dǎo) 致的。每一次分組奇數(shù)部分的系數(shù)WN,這里的N均為本次分組前的序列點(diǎn)數(shù)。以上邊的8點(diǎn)DFT為例，第一次分組N=8,第二次分組N為4，為了統(tǒng)一根據(jù)式 (4)進(jìn)行了變換將N變?yōu)榱?，但指數(shù)相應(yīng)的需要乘以2。

N點(diǎn)基-2 FFT算法的計(jì)算量

從 圖4可以看到N點(diǎn)DFT的FFT變換可以轉(zhuǎn)為log2(N)級(jí)級(jí)聯(lián)的蝶形運(yùn)算，每一級(jí)均包含有N/2次蝶形計(jì)算。而每一個(gè)蝶形運(yùn)算包含了1次復(fù)數(shù)乘法，2 次復(fù)數(shù)加法。因此N點(diǎn)FFT計(jì)算的總計(jì)算量為：復(fù)數(shù)乘法——N/2×log2(N)    復(fù)數(shù)加法——N×log2(N)。假設(shè)被采樣的序列為實(shí)數(shù)序列，那么也只有第一級(jí)的計(jì)算為實(shí)數(shù)與復(fù)數(shù)的混合計(jì)算，經(jīng)過一次迭代后來的計(jì)算均變?yōu)閺?fù)數(shù)計(jì)算， 在這一點(diǎn)上和直接的DFT計(jì)算不一致。因此對于輸入序列是復(fù)數(shù)還是實(shí)數(shù)對FFT算法的效率影響較小。一次復(fù)數(shù)乘法包含了4次實(shí)數(shù)乘法，2次實(shí)數(shù)加法，一次 復(fù)數(shù)加法包含了2次復(fù)數(shù)加法。因此對于N點(diǎn)的FFT計(jì)算需要總共的實(shí)數(shù)乘法數(shù)量為：2×N×log2(N)；總的復(fù)數(shù)加法次數(shù) 為：2xNxlog2(N)。

N點(diǎn)基-2 FFT算法的實(shí)現(xiàn)方法

從圖4我們可以總結(jié)出對于點(diǎn)數(shù)為N=2^L的DFT快速計(jì)算方法的流程：

   1.對于輸入數(shù)據(jù)序列進(jìn)行倒位序變換。

該 變換的目的是使輸出能夠得到X(0)~X(N-1)的順序序列，同樣以8點(diǎn)DFT為例，該變換將順序輸入序列x(0)~x(7)變?yōu)槿鐖D4的 x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其實(shí)現(xiàn)方法是：假設(shè)順序輸入序列一次村在A(0)~A(N-1) 的數(shù)組元素中，首先我們將數(shù)組下標(biāo)進(jìn)行二進(jìn)制化(例：對于點(diǎn)數(shù)為8的序列只需要LOG2(8) = 3位二進(jìn)制序列表示，序號(hào)6就表示為110)。二進(jìn)制化以后就是將二進(jìn)制序列進(jìn)行倒位，倒位的過程就是將原序列從右到左書寫一次構(gòu)成新的序列，例如序號(hào)為 6的二進(jìn)制表示為110，倒位后變?yōu)榱?11，即使十進(jìn)制的3。第三步就是將倒位前和倒位后的序號(hào)對應(yīng)的數(shù)據(jù)互換。依然以序號(hào)6為例，其互換過程如下：

temp = A(6);
A(6) = A(3);
A(3) = A(6);

從上式我們可以看到不同輸入項(xiàng)所乘的旋轉(zhuǎn)因子個(gè)數(shù)(注 意這里是個(gè)數(shù)，就算是wn^0，也被考慮進(jìn)去了，因?yàn)樵跊]有放大時(shí)wn^0等于1，放大后所有旋轉(zhuǎn)因子指數(shù)模均不為1，因此需要考慮)。這就導(dǎo)致輸入不平 衡，運(yùn)算結(jié)果不正確。經(jīng)查閱相關(guān)資料，比較妥善的做法是，首先對所有旋轉(zhuǎn)因子都放大2^Q倍，Q必須要大于等于L，以保證不同旋轉(zhuǎn)因子的差異化。旋轉(zhuǎn)因子 放大，為了保證其模為1，在每一次蝶形運(yùn)算的乘積運(yùn)算中我們需要將結(jié)果右移Q位來抵消這個(gè)放大，從而得到正確的結(jié)果。之所以采用放大倍數(shù)必須是2的整數(shù)次 冪的原因也在于此，我們之后可以通過簡單的右移位運(yùn)算將之前的放大抵消，而右移位又代替了除法運(yùn)算，大大節(jié)省了時(shí)間。

假設(shè)AD位寬12，數(shù)據(jù)位寬32，符號(hào)位1位，因此有效位寬31位，采樣點(diǎn)數(shù)N，那么我們可以得到log2(N)+P+Q<=19,假設(shè)點(diǎn)數(shù)128,又Q>=L可以得到放大倍數(shù)P<=5。

#define N           128
#define POWER       6//該值代表了輸入數(shù)據(jù)首先被放大的倍數(shù)，放大倍數(shù)為2^POWER
#define P_COEF      8//該值代表了旋轉(zhuǎn)因子被放大的倍數(shù)，放大倍數(shù)為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采樣位數(shù)為12位，本可以采用s16 x[N]定義數(shù)據(jù)空間，但是為了節(jié)省存儲(chǔ)空間，fft結(jié)果也將存儲(chǔ)在該變量空間。由于受
//fft影響變量需要重新定義，定義的方式及具體原因如下：
//fft過程會(huì)乘以較大旋轉(zhuǎn)因子，因此需要32位定義
//fft過程會(huì)產(chǎn)生負(fù)數(shù)，因此需要有符號(hào)定義
//fft過程會(huì)出現(xiàn)復(fù)數(shù)，因此需要定義二維數(shù)組，[][0]存放實(shí)部，[][1]存放虛部
s32 x[N][2] = {};

//定義*p[N]是為了在第一次指針初始化以后，該數(shù)組指針按照位倒序指向數(shù)據(jù)變量x
//p[i][0]代表了指向數(shù)據(jù)的實(shí)部，p[i][1]代表了指向數(shù)據(jù)的虛部
s32 *p[N];

//定義旋轉(zhuǎn)因子矩陣
//旋轉(zhuǎn)因子矩陣存儲(chǔ)了wn^0,wn^1,wn^2...wn^(N/2-1)這N/2個(gè)復(fù)數(shù)旋轉(zhuǎn)因子 
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主函數(shù),該函數(shù)將借助指針數(shù)組p對全局變量數(shù)組x進(jìn)行fft計(jì)算
*            并將結(jié)果存儲(chǔ)在數(shù)組x中
*global var:rw->x; r->p,w。(r表示讀，w表示寫，rw表示讀寫)
*/
void fft2(void)
{
  u8 i;//i用于表示蝶形圖級(jí)聯(lián)的階數(shù)
  u8 j;//表示蝶形分組起始點(diǎn)序列，蝶形分組跨度為2^i
  u8 k;//k表示蝶形組內(nèi)偶數(shù)分支序列點(diǎn)號(hào)
  u8 gp_distance = 1;//蝶形分組跨度
  u8 temp;//temp用于記錄當(dāng)前組間距離的一半,同時(shí)也表示了同一碟形兩分支間的距離
  u8 gp_hf = 0;//記錄當(dāng)前組的中間點(diǎn)序號(hào)
  u8 delta = N;//旋轉(zhuǎn)因子下標(biāo)增量，本來下標(biāo)初始值應(yīng)該為N/2，但是。。

  s16 *pw = &(w[0][0]);

  int tp_result[2];       //用于臨時(shí)存放旋轉(zhuǎn)因子和奇數(shù)分組的乘積

  //輸入信號(hào)序列放大
  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++)//完成一組內(nèi)的所有蝶形運(yùn)算
        {  //蝶形運(yùn)算中的一組復(fù)數(shù)乘法過程
           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;
           //蝶形運(yùn)算中的2組復(fù)數(shù)加法過程
           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;
  }
}