快速傅立葉變換(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ū)間長(zhǎng)度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)算方法,這就是快速傅立葉變換簡(jiǎ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)有限長(zhǎng)序列,其DFT正變換為
其反變換(IDFT)
x(n)=
二者的差別只在于
考慮x(n)為復(fù)數(shù)序列的一般情況,每計(jì)算一個(gè)X(k),需要N次復(fù)數(shù)乘法以及(N-1)次復(fù)數(shù)加法。因此,對(duì)所有N個(gè)k值,共需N2次復(fù)數(shù)乘法及
N(N-1)次復(fù)數(shù)加法運(yùn)算。所以直接計(jì)算DFT,乘法次數(shù)和加法次數(shù)都是和N2成正比的,當(dāng)N很大時(shí),運(yùn)算量是很可觀的,因而需要改進(jìn)對(duì)DFT的計(jì)算方法,以減少運(yùn)算次數(shù)。
下面討論減少運(yùn)算工作量的途徑。仔細(xì)觀察DFT的運(yùn)算就可看出,利用系數(shù)
(1)
(2)
(3)
由此可得:
二、時(shí)域抽取法基-2 FFT原理
先設(shè)序列點(diǎn)數(shù)為N=2M,M為整數(shù)。如果不滿足這個(gè)條件,可以人為地加上若干零值點(diǎn),使之達(dá)到這一要求。這種N為2的整數(shù)冪的FFT稱基-2 FFT。
設(shè)輸入序列長(zhǎng)度為N=2M (M為正整數(shù)) ,將該序列按時(shí)間順序的奇偶分解為越來越短的子序列,稱為按時(shí)間抽取(DIT )的FFT算法。也稱Cooley - Tukey算法。
將N=2M的序列x(n)先按n的奇偶分成以下兩組:
x(2r+1)=x2(r)
則x(n)的DFT為:
=
=
=
=
式中
由(4.2.4)式可看出,一個(gè)N點(diǎn)DFT已分解成兩個(gè)N/2點(diǎn)的DFT,它們按(4.2.4)式又組合成一個(gè)N點(diǎn)DFT。
現(xiàn)討論
式中用了
同理可得:
再考慮
所以前半部分:
后半部分:
=
這樣,只要求出0到(N/2-1)區(qū)間的所有
采用蝶形運(yùn)算符號(hào)表示的圖示法,可將上面討論的分解過程表示于下圖中。此圖表示N=23=8的情況,其中輸出值X(0)到X(3)是由(4.2.7)式給出的,而輸出值X(4)到X(7)是由(4.2.8)給出。
既然如此,由于N=2M,因而N/2仍是偶數(shù),可以進(jìn)一步把每個(gè)N/2點(diǎn)子序列再按其奇偶部分分解為兩個(gè)N/4點(diǎn)的子序列。
x1(2r+1)=x4(r)
=
=
且
同理,
且
其中,
進(jìn)一步具體化:N=8=23,
=x(0)+x(4)
注意上式中
同理,
三、DIT-FFT算法與直接計(jì)算DFT運(yùn)算量的比較
由上圖得,當(dāng)N=2M時(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ù)乘法為
當(dāng)N=1024時(shí),復(fù)數(shù)乘之比:
顯然,N越大,優(yōu)越性就越明顯。
四、按時(shí)間抽選的FFT算法的特點(diǎn):
1、 原位運(yùn)算
由圖4.2.4可以看出,DIT-FFT的運(yùn)算過程很有規(guī)律。N=2M點(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ǔ)單元中便依次存放
2、 倒序規(guī)律
由圖4.2.4看出,按原位計(jì)算時(shí),FFT的輸出
造成倒序的原因是輸入x(n)按標(biāo)號(hào)n的偶奇的不斷分組而造成。由于N=2M,所以倒序數(shù)可用M位二進(jìn)制數(shù)(nM-1nM-2…n0)2(當(dāng)N=8=23時(shí),二進(jìn)制為三位)表示。第一次分組,標(biāo)為n0。 n為偶數(shù)在上半部分,用n0=0表示,n為奇數(shù)在下半部分,用n0=1表示;第二次分組,標(biāo)為n1。偶數(shù)部分再分為偶(0)奇(1),奇數(shù)部分再分為偶(0)奇(1)…。依次類推,直到M次分組,最后所得二進(jìn)制倒序數(shù)如圖示。
下表列出了N=8時(shí)以二進(jìn)制數(shù)表示的順序數(shù)和倒序數(shù),由表顯而易見,只要將順序數(shù)(n2n1n0)的二進(jìn)制位倒置,則得對(duì)應(yīng)的二進(jìn)制倒序值(n0n1n2)。
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=2M,共有第一級(jí)蝶形運(yùn)算,第二級(jí)蝶形運(yùn)算,…,第M級(jí)蝶形運(yùn)算。用L表示第某級(jí)運(yùn)算,每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)的“距離”為B=2L-1。
5、 旋轉(zhuǎn)因子的變化規(guī)律
仍觀察圖4.2.4,每級(jí)都有N/2個(gè)蝶形,每個(gè)蝶形都要乘以因子
L=1,第一級(jí):
L=2,第二級(jí):
L=3,第三級(jí):
所以對(duì)于N=2M的一般情況,第L級(jí)的
由于2L=2M2L-M=N2L-M ,所以
所以,第L級(jí)的
例如:
L=1,第一級(jí):2L-1=1,J=0,
L=2,第二級(jí):2L-1=2,J=0,1,
L=3,第三級(jí):2L-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)算可表示成如下形式:
式中
L=1,…,M。
DIT-FFT運(yùn)算和程序框圖如下:
同一旋轉(zhuǎn)因子對(duì)應(yīng)著間隔為2L點(diǎn)的2M-L個(gè)蝶形。
五、按頻率抽選(DIF)的基2 FFT算法
設(shè)序列點(diǎn)數(shù)為N=2M ,M為整數(shù)。
先把輸入按n 的順序分成前后兩半:
=
=
=
-1,k為奇數(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),
=
(4.2.16)用下圖表示,N=8。一次分解流圖。
由于N=2M,N/2仍是偶數(shù),繼續(xù)將N/2點(diǎn)DFT分成偶數(shù)組和奇數(shù)組。圖4.2.12表示N=8時(shí)二次分解運(yùn)算流圖。
最后完整的分解流圖為下圖:
這種算法是對(duì)
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)=
與離散傅立葉正變換(
如果希望直接調(diào)用FFT子程序計(jì)算IFFT,則可用下面的方法:
由于 x(n)=
∴ x*(n)=
x(n)=
=
這樣可以先將
4.3 進(jìn)一步減少運(yùn)算量的措施
研究進(jìn)一步減少運(yùn)算量的途徑,以程序的復(fù)雜度換取計(jì)算量的進(jìn)一步提高
一、 多類蝶形單元運(yùn)算
由圖4.2.4已得出結(jié)論,N=2M點(diǎn)FFT共需要MN/2次復(fù)數(shù)乘法。
由
綜上所述,先除去第一、第二兩級(jí)后,所需復(fù)數(shù)乘法次數(shù)
進(jìn)一步考慮各級(jí)中的無關(guān)緊要旋轉(zhuǎn)因子。當(dāng)L=3時(shí),有兩個(gè)無關(guān)緊要的旋轉(zhuǎn)因子
因此
在基2 FFT程序中,若包含了所有旋轉(zhuǎn)因子,則稱該算法為一類蝶形單元運(yùn)算;若去掉
二、 旋轉(zhuǎn)因子的生成
在FFT運(yùn)算中,旋轉(zhuǎn)因子
三、 實(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。
聯(lián)系客服