快速傅立葉變換(FFT)

4．1引言

快速傅立葉變換（FFT）并不是一種新的變換，而是離散傅立葉變換（DFT）的一種快速算法。

DFT的計(jì)算在數(shù)字信號(hào)處理中非常有用。例如在FIR濾波器設(shè)計(jì)中會(huì)遇到從h(n)求H(k)或由H(k)計(jì)算h(n)，這就要計(jì)算DFT；信號(hào)的譜分析對(duì)通信、圖像傳輸、雷達(dá)等都是很重要的，也要計(jì)算DFT。因直接計(jì)算DFT的計(jì)算量與變換區(qū)間長度N的平方成正比，當(dāng)N較大時(shí)，計(jì)算量太大。

自從1965年圖基（J. W. Tukey）和庫利（T. W. Coody）在《計(jì)算數(shù)學(xué)》（Math. Computation , Vol. 19, 1965）雜志上發(fā)表了著名的《機(jī)器計(jì)算傅立葉級(jí)數(shù)的一種算法》論文后，桑德(G. Sand)-圖基等快速算法相繼出現(xiàn)，又經(jīng)人們進(jìn)行改進(jìn)，很快形成一套高效運(yùn)算方法，這就是快速傅立葉變換簡稱FFT（Fast Fourier Transform）。這種算法使DFT的運(yùn)算效率提高1~2個(gè)數(shù)量級(jí)。

4．2 基2 FFT算法

一、直接計(jì)算DFT的問題及改進(jìn)的途徑

設(shè)x(n)為N點(diǎn)有限長序列，其DFT正變換為

= ， k=0，1，…，N-1

其反變換(IDFT)

x(n)= ，n=0，1，…，N-1

二者的差別只在于 的指數(shù)符號(hào)不同，以及差一個(gè)常數(shù)乘因子1/N，因而下面我們只討論DFT正變換的運(yùn)算量，反變換的運(yùn)算量是完全相同的。

考慮x(n)為復(fù)數(shù)序列的一般情況，每計(jì)算一個(gè)X(k)，需要N次復(fù)數(shù)乘法以及（N-1）次復(fù)數(shù)加法。因此，對(duì)所有N個(gè)k值，共需N²次復(fù)數(shù)乘法及
N（N-1）次復(fù)數(shù)加法運(yùn)算。所以直接計(jì)算DFT，乘法次數(shù)和加法次數(shù)都是和N²成正比的，當(dāng)N很大時(shí)，運(yùn)算量是很可觀的，因而需要改進(jìn)對(duì)DFT的計(jì)算方法，以減少運(yùn)算次數(shù)。

       下面討論減少運(yùn)算工作量的途徑。仔細(xì)觀察DFT的運(yùn)算就可看出，利用系數(shù) 以下固有特性，就可減小DFT的運(yùn)算量：

（1） 的對(duì)稱性 （ ）^*=

（2） 的周期性 = =

（3） 的可約性 = =

由此可得： = = ， =-1， =- 。這樣利用這些特性，可以將長序列的DFT分解為短序列的DFT。快速傅立葉變換算法正是基于這樣的思路而發(fā)展起來的。它的算法基本上可以分成兩大類：時(shí)域抽取法FFT（Decimation-In-Time FFT，簡稱DIT-FFT）和頻域抽取法FFT（Decimation-In-Time FFT，簡稱DIF-FFT）。

二、時(shí)域抽取法基-2 FFT原理

先設(shè)序列點(diǎn)數(shù)為N=2^M，M為整數(shù)。如果不滿足這個(gè)條件，可以人為地加上若干零值點(diǎn)，使之達(dá)到這一要求。這種N為2的整數(shù)冪的FFT稱基-2 FFT。

設(shè)輸入序列長度為N=2^M (M為正整數(shù)) ，將該序列按時(shí)間順序的奇偶分解為越來越短的子序列，稱為按時(shí)間抽取(DIT )的FFT算法。也稱Cooley - Tukey算法。

將N=2^M的序列x(n)先按n的奇偶分成以下兩組:

x(2r)=x₁(r)      r=0，1，…，N/2-1

x(2r+1)=x₂(r)

則x(n)的DFT為：

＝DFT[x(n)]= ，k=0，1，…，N-1

 = +

 = +

 = +

= +

 = + ，k=0，1，…，N-1 （4.2.4）

式中 與 分別是x₁(r)及x₂(r)的N/2點(diǎn)DFT：

      = = ，k=0，1，…，N/2-1

= = ，k=0，1，…，N/2-1

由（4.2.4）式可看出，一個(gè)N點(diǎn)DFT已分解成兩個(gè)N/2點(diǎn)的DFT，它們按（4.2.4）式又組合成一個(gè)N點(diǎn)DFT。

現(xiàn)討論 的周期性：

= = =

式中用了 =

同理可得： =

再考慮 的以下性質(zhì)：

       = =－

所以前半部分：

   = + ，k=0，1，…，N/2-1

后半部分：

      = +

= － ，（k=0，1，…，N/2-1）

這樣，只要求出0到（N/2-1）區(qū)間的所有 和 值，即可求出0到（N-1）區(qū)間內(nèi)的所有 值，這就大大節(jié)省了運(yùn)算。

= + ，k=0，1，…，     (4.2.7)

= － ，k=0，1，…， (4.2.8)

采用蝶形運(yùn)算符號(hào)表示的圖示法，可將上面討論的分解過程表示于下圖中。此圖表示N=2³=8的情況，其中輸出值X（0）到X（3）是由（4.2.7）式給出的，而輸出值X（4）到X（7）是由（4.2.8）給出。

既然如此，由于N=2^M，因而N/2仍是偶數(shù)，可以進(jìn)一步把每個(gè)N/2點(diǎn)子序列再按其奇偶部分分解為兩個(gè)N/4點(diǎn)的子序列。

     x₁(2r) = x₃(r)      r=0，1，…，N/4-1

x₁(2r+1)=x₄(r)

= +

    = +

    = + ，k=0，1，…，N/4-1

且 = - ，k=0，1，…，N/4-1

同理， 也可進(jìn)行同樣的分解：

      = + ，k=0，1，…，N/4-1

且 = - ，k=0，1，…，N/4-1

其中， = =

           = =

進(jìn)一步具體化：N=8=2³，

     =

         =x(0)+x(4) ，k=0，1

    = x(0)+ x(4)

         = x(0)+ x(4) = x(0)- x(4)

注意上式中 =

同理， = x(1)+ x(5)

           = x(1)- x(5)

三、DIT-FFT算法與直接計(jì)算DFT運(yùn)算量的比較

由上圖得，當(dāng)N=2^M時(shí)，其運(yùn)算流圖有M級(jí)蝶形，每一級(jí)都有N/2個(gè)蝶形運(yùn)算構(gòu)成。每一個(gè)蝶形運(yùn)算需要1次復(fù)數(shù)乘和2次復(fù)數(shù)加。所以每一級(jí)運(yùn)算都需要N/2次復(fù)數(shù)乘和N次復(fù)數(shù)加。

M級(jí)運(yùn)算總共需要的復(fù)數(shù)乘法次數(shù)為：

M級(jí)運(yùn)算總共需要的復(fù)數(shù)加法次數(shù)為：

如直接計(jì)算DFT，復(fù)數(shù)乘法為 次，復(fù)數(shù)加法 次。

當(dāng)N=1024時(shí)，復(fù)數(shù)乘之比： ，這樣使運(yùn)算效率提高200多倍。下圖為FFT算法和直接DFT算法所需運(yùn)算量與計(jì)算點(diǎn)數(shù)N的關(guān)系曲線：

顯然，N越大，優(yōu)越性就越明顯。

四、按時(shí)間抽選的FFT算法的特點(diǎn)：

1、        原位運(yùn)算

由圖4.2.4可以看出，DIT-FFT的運(yùn)算過程很有規(guī)律。N=2^M點(diǎn)的FFT共進(jìn)行M級(jí)運(yùn)算，每級(jí)由N/2個(gè)蝶形運(yùn)算組成。同一級(jí)中，每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)只對(duì)計(jì)算本蝶形有用，而且每個(gè)蝶形的輸入、輸出數(shù)據(jù)結(jié)點(diǎn)又同在一條水平線上，這就意味著計(jì)算完一個(gè)蝶形后，所得輸出數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的存儲(chǔ)單元。這樣，經(jīng)過M級(jí)運(yùn)算后，原來存放輸入序列數(shù)據(jù)的N個(gè)存儲(chǔ)單元中便依次存放 的N個(gè)值。這種利用同一存儲(chǔ)單元存儲(chǔ)蝶形計(jì)算輸入、輸出數(shù)據(jù)的方法稱為原位計(jì)算。原位計(jì)算可節(jié)省大量內(nèi)存，從而使設(shè)備成本降低。

2、        倒序規(guī)律

由圖4.2.4看出，按原位計(jì)算時(shí)，FFT的輸出 是按正常順序排列在存儲(chǔ)單元中，即按 ， ，…， 的順序排列，但是這時(shí)輸入x(n)卻不是按自然順序存儲(chǔ)的，而是按x(0)，x(4)，…，x(7)的順序存入存儲(chǔ)單元，看起來好象是“混亂無序”的，實(shí)際上是有規(guī)律的，我們稱之為倒序。

        造成倒序的原因是輸入x(n)按標(biāo)號(hào)n的偶奇的不斷分組而造成。由于N=2^M，所以倒序數(shù)可用M位二進(jìn)制數(shù)（n_M-1n_M-2…n₀）₂（當(dāng)N=8=2³時(shí)，二進(jìn)制為三位）表示。第一次分組，標(biāo)為n₀。 n為偶數(shù)在上半部分，用n₀=0表示，n為奇數(shù)在下半部分，用n₀=1表示；第二次分組，標(biāo)為n₁。偶數(shù)部分再分為偶（0）奇（1），奇數(shù)部分再分為偶（0）奇（1）…。依次類推，直到M次分組，最后所得二進(jìn)制倒序數(shù)如圖示。

下表列出了N=8時(shí)以二進(jìn)制數(shù)表示的順序數(shù)和倒序數(shù)，由表顯而易見，只要將順序數(shù)（n₂n₁n₀）的二進(jìn)制位倒置，則得對(duì)應(yīng)的二進(jìn)制倒序值（n₀n₁n₂）。

3、        倒序的實(shí)現(xiàn)

設(shè)原輸入序列x(n)先按自然順序存入數(shù)組A中。例如N=8，A(0)，A(1)，A(2)，A(3)，A(4)，A(5)，A(6)，A(7)中依次存放著x(0)，x(1)，x(2)，x(3)，x(4)，x(5)，x(6)，x(7)。

順序數(shù)用I表示，I=1~N-2。倒序數(shù)用J表示，與I對(duì)應(yīng)分別為4，2，6，1，5，3。當(dāng)I=J時(shí)不需要交換，I<J時(shí)調(diào)換存放內(nèi)容。

    I=1時(shí)，對(duì)應(yīng)的倒序數(shù)是4；I=2時(shí)，對(duì)應(yīng)的倒序數(shù)是2…。倒序數(shù)從4到2到…的關(guān)系可從表4.2.1得到：每次最高位加1。（注意J用十進(jìn)制數(shù)表示）。如果最高位為0，J直接加N/2，如果最高位為1，則要將最高位歸0，次高位加1。但次高位加1時(shí)也要判斷是否為1或0。程序框圖如下圖虛線框里所示：

4、        蝶形運(yùn)算兩個(gè)輸入數(shù)據(jù)的“距離”

以圖4.2.4的8點(diǎn)FFT為例，其輸入是倒位序的，輸出是自然順序的。N=2^M，共有第一級(jí)蝶形運(yùn)算，第二級(jí)蝶形運(yùn)算，…，第M級(jí)蝶形運(yùn)算。用L表示第某級(jí)運(yùn)算，每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)的“距離”為B=2^L-1。

5、        旋轉(zhuǎn)因子的變化規(guī)律

仍觀察圖4.2.4，每級(jí)都有N/2個(gè)蝶形，每個(gè)蝶形都要乘以因子 ，稱其為旋轉(zhuǎn)因子，p稱為旋轉(zhuǎn)因子的指數(shù)。第L級(jí)蝶形運(yùn)算（從左向右數(shù)），共有2^L-1個(gè)

L=1，第一級(jí)： = ，J=0

L=2，第二級(jí)： = ，J=0，1

L=3，第三級(jí)： = ，J=0，1，2，3

所以對(duì)于N=2^M的一般情況，第L級(jí)的 為：

，J=0，1，…，2^L-1-1

由于2^L=2^M2^L-M=N2^L-M ，所以 = = =

所以，第L級(jí)的 為： ，J=0，1，…，2^L-1-1

例如：

L=1，第一級(jí)：2^L-1=1，J=0，      

L=2，第二級(jí)：2^L-1=2，J=0，1，    ，

L=3，第三級(jí)：2^L-1=4，J=0，1，2，3， ， ， ，

6、蝶形運(yùn)算規(guī)律

    設(shè)序列x(n)經(jīng)時(shí)域抽選（倒序）后，存入數(shù)組X中。如果蝶形運(yùn)算的兩個(gè)輸入數(shù)據(jù)相距B個(gè)點(diǎn)，應(yīng)用原位計(jì)算，則蝶形運(yùn)算可表示成如下形式：

式中 ；J=0，1，…，2^L-1-1；

L=1，…，M。

DIT-FFT運(yùn)算和程序框圖如下：

同一旋轉(zhuǎn)因子對(duì)應(yīng)著間隔為2^L點(diǎn)的2^M-L個(gè)蝶形。

五、按頻率抽選（DIF）的基2 FFT算法

設(shè)序列點(diǎn)數(shù)為N=2^M ，M為整數(shù)。

＝ ，k=0，1，…，N-1

先把輸入按n 的順序分成前后兩半：

＝ ，k=0，1，…，N-1

= +

    = +

    = ，k=0，1，…，N-1

= ，k=0，1，…，N-1

=(-1)^k=   1，k為偶數(shù)

                  -1，k為奇數(shù)

將 分解成偶數(shù)組與奇數(shù)組，

當(dāng)k取偶數(shù)即k=2r時(shí)（r=0，1，…，N/2-1），

= =

當(dāng)k取奇數(shù)即k=2r+1時(shí)（r=0，1，…，N/2-1），

=

=

令         n=0，1，…，N/2-1

= ，r=0，1，…，N/2-1

=         (4.2.16)

、 和 之間可用下圖所示的蝶形運(yùn)算符號(hào)表示：

（4.2.16）用下圖表示，N=8。一次分解流圖。

由于N=2^M，N/2仍是偶數(shù)，繼續(xù)將N/2點(diǎn)DFT分成偶數(shù)組和奇數(shù)組。圖4.2.12表示N=8時(shí)二次分解運(yùn)算流圖。

最后完整的分解流圖為下圖：

    這種算法是對(duì) 進(jìn)行奇偶抽取分解的結(jié)果，所以稱之為頻域抽取法FFT。

    DIF-FFT算法與DIT-FFT算法類似，不同的是DIF-FFT算法輸入序列為自然順序，而輸出為倒序排列。另外，蝶形運(yùn)算略有不同，DIT-FFT蝶形先乘后加，而DIF-FFT蝶形先加后乘。

    上述兩種FFT的算法流圖形式不是唯一的。只要保證各節(jié)點(diǎn)所連支路及其傳輸系數(shù)不變，改變輸入與輸出點(diǎn)以及中間結(jié)點(diǎn)的排列順序，就可以得到其他變形的FFT運(yùn)算流圖。

六、IDFT的高效算法

   離散傅立葉反變換

x(n)= ，n=0，1，…，N-1

與離散傅立葉正變換（ = ， k=0，1，…，N-1）相比，只要將DFT中的系數(shù) 改變?yōu)?/span> ，最后乘以1/N，就是IDFT的運(yùn)算公式。流圖輸入為 ，輸出為x(n)。因此原來的DIT-FFT改為IFFT后稱為DIF-IFFT更合適；原來的DIF-FFT改為IFFT后稱為DIT-IFFT更合適。

下圖是由DIF-FFT運(yùn)算流圖改成的IFFT運(yùn)算流圖（DIT-IFFT）：

    在實(shí)際中，有時(shí)為了防止運(yùn)算過程中發(fā)生溢出，將1/N分配到每一級(jí)蝶形運(yùn)算中。由于1/N= ，所以每級(jí)的每個(gè)蝶形輸出支路均有一相乘因子1/2。如下圖示：

如果希望直接調(diào)用FFT子程序計(jì)算IFFT，則可用下面的方法：

由于    x(n)=

∴      x*(n)=

x(n)=

   = {DFT[ ]}*

這樣可以先將 取共軛，然后直接調(diào)用FFT子程序進(jìn)行DFT運(yùn)算，最后取共軛并乘以1/N得到序列x(n)。

4．3 進(jìn)一步減少運(yùn)算量的措施

研究進(jìn)一步減少運(yùn)算量的途徑，以程序的復(fù)雜度換取計(jì)算量的進(jìn)一步提高

一、   多類蝶形單元運(yùn)算

由圖4.2.4已得出結(jié)論，N=2^M點(diǎn)FFT共需要MN/2次復(fù)數(shù)乘法。

由 = ，J=0，1，…，2^L-1-1得，當(dāng)L=1時(shí)，只有一種旋轉(zhuǎn)因子 ，所以第一級(jí)不需要乘法運(yùn)算。當(dāng)L=2時(shí)，共有兩個(gè)旋轉(zhuǎn)因子： =1和 =-j，因此第二級(jí)也不需要乘法運(yùn)算。在DFT中，稱其值為±1和±j的旋轉(zhuǎn)因子為無關(guān)緊要的旋轉(zhuǎn)因子。

    綜上所述，先除去第一、第二兩級(jí)后，所需復(fù)數(shù)乘法次數(shù)

進(jìn)一步考慮各級(jí)中的無關(guān)緊要旋轉(zhuǎn)因子。當(dāng)L=3時(shí)，有兩個(gè)無關(guān)緊要的旋轉(zhuǎn)因子 和 。因?yàn)橥恍D(zhuǎn)因子對(duì)應(yīng)著2^M-L=N/2^L個(gè)蝶形運(yùn)算，所以，第三級(jí)共有2N/2³=N/4個(gè)蝶形不需要復(fù)數(shù)乘法運(yùn)算。依次類推，當(dāng)L≥3時(shí)，第L級(jí)的2個(gè)無關(guān)緊要的旋轉(zhuǎn)因子減少復(fù)數(shù)乘法的次數(shù)為2N/2^L=N/2^L-1。這樣，L=3至L=M共減少復(fù)數(shù)乘法次數(shù)為

因此

        =

在基2 FFT程序中，若包含了所有旋轉(zhuǎn)因子，則稱該算法為一類蝶形單元運(yùn)算；若去掉 =±1的旋轉(zhuǎn)因子，則稱之為二類蝶形單元運(yùn)算；若再去掉 =±j的旋轉(zhuǎn)因子，則稱為三類蝶形單元運(yùn)算；若再處理 = ，則稱之為四類蝶形運(yùn)算。我們將后三種運(yùn)算稱為多類蝶形單元運(yùn)算。顯然蝶形單元越多，編程就越復(fù)雜，但當(dāng)N較大時(shí)，乘法運(yùn)算的減少量是相當(dāng)可觀的。例如，N=4096時(shí)，三類蝶形單元運(yùn)算的乘法次數(shù)為一類蝶形單元運(yùn)算的75%。

二、   旋轉(zhuǎn)因子的生成

在FFT運(yùn)算中，旋轉(zhuǎn)因子

= ，求余弦和正弦函數(shù)值的計(jì)算量很大，所以編程時(shí)，一種方法是在每級(jí)運(yùn)算中直接產(chǎn)生，另一種方法是在FFT程序開始前預(yù)先計(jì)算好，存放在數(shù)組中，作為旋轉(zhuǎn)因子表，在程序執(zhí)行過程中直接查表得到。

三、   實(shí)序列的FFT算法

在實(shí)際工作中，數(shù)據(jù)x(n)一般都是實(shí)序列。如果直接按FFT運(yùn)算流圖計(jì)算，就是把x(n)看成一個(gè)虛部為零的復(fù)序列進(jìn)行計(jì)算，這就增加了運(yùn)算時(shí)間。處理這個(gè)問題有二種方法，一種是早期提出的用一個(gè)N點(diǎn)FFT計(jì)算N點(diǎn)實(shí)序列的FFT。第二種方法是用N/2點(diǎn)FFT計(jì)算一個(gè)N點(diǎn)實(shí)序列的DFT。