數(shù)字信號處理原理與實踐 第4章

第四章快速傅里葉變換數(shù)字信號處理原理與實踐(修訂版)@NEU離散傅里葉變換（DFT）在數(shù)字信號處理中是最常用的運算之一，由于DFT計算比較繁瑣，在相當(dāng)長的時間內(nèi)并沒有得到真正的應(yīng)用。直到1965年庫利(J.W.Cooley)和圖基(J.W.Tukey)提出了DFT的一種快速算法，又有桑德(G.Sande)和圖基的快速算法相繼出現(xiàn)，之后又出現(xiàn)了各種各樣的DFT算法，這些方法統(tǒng)稱為快速傅里葉變換（FFT）。4.1離散傅里葉變換存在的問題對于N點序列，其DFT為一般情況下，和都是復(fù)數(shù)，每計算一個值，需要作N次復(fù)數(shù)乘法和(N-1)次復(fù)數(shù)加法，求N點的值，需要N2次復(fù)數(shù)乘法，以及N(N-1)次復(fù)數(shù)加法。由于實現(xiàn)一次復(fù)數(shù)乘法需要四次實數(shù)乘法和兩次實數(shù)加法，一次復(fù)數(shù)加法需要兩次實數(shù)加法。

在實際運算過程中，運算量并沒有如此之大，因為在DFT運算中包含著大量的重復(fù)運算，同時一些系數(shù)是不需要乘法運算的。4.2按時間抽取的基2FFT算法

利用的周期性，對稱性，把長度為N點的DFT運算逐次分解為較短序列的DFT運算。這種方法叫時間抽取法(decimation-in-time,縮寫為DIT)，簡稱時選法。4.2.1算法的推導(dǎo)

設(shè)N＝2M，M為正整數(shù)，N是2的M次方，如果不滿足這個條件，可以人為的加上若干個零值點，達到這個要求。這種N為2的整數(shù)次冪的FFT稱為基2FFT。將按照奇數(shù)和偶數(shù)分解為兩個序列，將N點的DFT運算分解為兩個N/2點的DFT運算，即令n=2r，n=2r+1，r=0，1，2，….N/2-1DFT運算分為對應(yīng)的兩組，于是有

k=0，1，…N/2-1k=0，1，…N/2-1k=0，1，…N/2-1令則后半部分的這樣，N點的就可以用和來表示，分別為N/2個點的DFT。運算可以用蝶形信號流圖表示。求N點的值，需要N2次復(fù)數(shù)乘法和N(N-1)次復(fù)數(shù)加法。對于N/2點DFT有(N/2)2=N2/4次復(fù)數(shù)乘法和(N/2)(N/2-1)次復(fù)數(shù)加法，因此，兩個N/2點DFT有N2/2次復(fù)數(shù)乘法和N(N/2-1)次復(fù)數(shù)加法。把兩個N/2的DFT合成一個N點DFT時，需要N/2個蝶型運算，包括N/2次復(fù)數(shù)乘法和N次復(fù)數(shù)加法。由于N>>1，總共的復(fù)數(shù)乘法為，復(fù)數(shù)加法次數(shù)為，通過這樣的分解，運算量幾乎可以減少到一半?？梢园瓷鲜龇椒ɡ^續(xù)分解，將N/2點的DFT按奇偶部分分解為兩個N/4點的DFT，進一步減少復(fù)數(shù)乘法的運算量。以N＝8為例，將一個8點的DFT分解為兩個4點的DFT，其分解過程如圖所示同樣，再把每個4點的DFT分解為兩個點的DFT，對于兩點的DFT來說，已經(jīng)不需要再分解。因此得到以下結(jié)果：當(dāng)N＝8時，將一個N/2點的DFT分解成兩個N/4個點的DFT的分解過程如圖所示。這樣一個8點的DFT就分解成為四個2點的DFT。如果N＝16，32，或者2的更高次冪，則可以按上述方法繼續(xù)分解下去，直到只有兩點DFT為止。4.2.2算法的討論

用上面所述的方法，任何一個N＝2M點的DFT，都可以經(jīng)過M次分解，最后分成兩個點的DFT運算，如圖所示，每分一次，稱為一“級”運算，所以N點的DFT可以分為M級。由于每一級都含有N/2個蝶型單元，每一個蝶型單元又只需要一次復(fù)數(shù)乘法，兩次復(fù)數(shù)加法，則完成一次時選FFT的總計算量為：復(fù)數(shù)乘法次數(shù)：Mf=(N/2)log2N=MN/2復(fù)數(shù)加法次數(shù)：Af=Nlog2N=MN當(dāng)全部FFT完成后，變換后輸出的序列依照正序排列，即存儲單元中順序放著X(0),X(1),…,X(N-1)。但輸入?yún)s不是原來的自然順序，而是x(0),x(4),x(2),…,x(N-1)是由于作奇偶分開所產(chǎn)生的。以N＝8為例，其自然序號是0，1，2，3，4，5，6，7。第一次按奇偶分開，得到兩組N/2點的DFT，的序列號是0，2，4，6和1，3，5，7。對每一組再做奇偶分開，這時抽取后得到四組，每組的序號是0，4和2，6和1，5和3，7。這種從十進制看來很亂的輸入順序,實際上是按二進制“倒序位”排列的。仍然以N=8為例,x(0),x(1),…,x(7)對應(yīng)是x(000),x(001),x(010),x(011),x(100),x(101),x(110),x(111),將二進制碼翻轉(zhuǎn)得到x(000),x(100),x(010),x(110),x(001),x(101),x(011),x(111),它們對應(yīng)的十進制序號分別是x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)。4.3按頻率抽取基2FFT算法

與按時間抽取相對應(yīng)，常用的FFT算法還有基2頻率抽?。╠ecimation-in-frequency,縮寫為DIF）FFT算法，它不是將按時序的奇偶分解，而是將頻域的序號按奇偶分開，下面說明頻率抽取的基2FFT算法。令序列點數(shù)N＝2M，M為整數(shù)，先將N點的DFT分成上下兩個部分，即式中將X(k)分解為偶數(shù)組和奇數(shù)組,分別令k=2r,k=2r+1,r=0,1,2,…,N/2-1,則令得到兩個N/2點DFT這樣，一個N點的DFT被分解為兩個N/2點的DFT。當(dāng)N＝8時，上述過程如圖所示：將一個N點DFT分解為4個N/4點DFT的過程如圖所示。是一個完整的N＝8，按頻率抽選的FFT結(jié)構(gòu)如圖所示。4.4運算量進一步減少的方法

雖然前面討論的FFT算法，已經(jīng)大大減少了運算量，但是還可以通過一些措施，進一步減少運算量。對于N＝2M，一共需要M級運算，每級N/2個蝶型單元，每個蝶型單元進行一次復(fù)數(shù)乘法，則總共需要MN/2次復(fù)數(shù)乘法，如果＝1，－1，j，-j，則與這些因子相乘時不作實際的復(fù)數(shù)乘法運算。4.5IDFT的快速計算方法IFFT

IDFT的快速計算方法是快速傅里葉反變換，用英文縮寫IFFT來表示，無論時選法還是頻選法的FFT算法均可用于IDFT運算。對于N點序列，其DFT運算公式為其IDFT運算公式為只要將DFT運算中的改為，在將結(jié)果乘以1/N，就可以用FFT的計算方法來實現(xiàn)IFFT。頻域序列為輸入，時域序列為輸出，所以時選FFT對應(yīng)時選的IFFT，頻選FFT對應(yīng)頻選IFFT。對于系數(shù)1/N要做如下處理，把1/N分解為（1/2）M，在M的每一級運算都要乘以1/2，這樣做的目的是為了避免在運算過程中出現(xiàn)溢出。其中N為IFFT變換的點數(shù)，N＝2M，且M為整數(shù)。IFFT的流圖如圖4.6分裂基FFT算法

分裂基FFT算法又稱基2/4算法，或混合基算法，它和前面講的基2算法有關(guān)，也和基4算法有關(guān)?；?算法比基2FFT算法更快，從理論上說，用較大的基數(shù)可以進一步減少運算次數(shù)，但程序過于復(fù)雜，所以一般不取基數(shù)大于8。1984年，法國的杜梅爾和霍爾曼首先提出了“分裂基”算法，把基2算法和基4算法結(jié)合起來，其基本思路是對偶數(shù)序號輸出采用基2算法，對奇數(shù)數(shù)序號輸出采用基4算法。由于分裂基算法是目前已知的對于N＝2M算法中具有最少的加法次數(shù)和乘法次數(shù)，并且具有非常好的結(jié)構(gòu)，因此被認為是最好的FFT算法。后來的研究表明，該算法最接近理論上的所需乘法次數(shù)的最小值。首先介紹一下基4FFT算法。如果N＝16，通過上述推導(dǎo)可以把16點的DFT分成4個4點的DFT，其信號第一級流圖如圖所示?？梢杂猛瑯拥姆椒▽1(k)、X2(k)、X3(k)做第三級分解,把r為偶序號的x1(r)作N/4點的DFT,把r為奇序號的x1(r)作N8點的DFT。同樣的,對X2(k)、X3(k)做同樣的處理。以N=16為例,X(k)的第一級的分解可以得到4個分裂基,X1(k)的第二級分解可以得到2個分裂基,一個基-4的4點DFT和2個基-2的2點DFT,而X2(k)和X3(k)的第二級的分解分別是基-4的4點DFT,對N=16的分裂基FFT的示意圖如圖所示。4.7快速傅里葉變換的MATLAB程序?qū)崿F(xiàn)

MATLAB為計算快速傅里葉變換，提供了一系列豐富的數(shù)學(xué)函數(shù)，主要有fft，Ifft，fft2，ifft2等。當(dāng)所處理的數(shù)據(jù)的長度為2的冪時，采用基—2算法進行計算，計算速度會顯著增加。所以，要盡可能使所要處理的數(shù)據(jù)長度為2的整數(shù)次冪或者用添零方式來添補數(shù)據(jù)使其長度成為2的整數(shù)次冪。1．fft和ifft函數(shù)調(diào)用方式(1)Y＝fft(x)參數(shù)說明如果x是向量，則采用傅里葉變換來求解x的離散傅里葉變換；如果x是矩陣，則計算該矩陣每一列的離散傅立葉變換；(2)Y=fft(x，N)參數(shù)說明N是進行離散傅立葉變換的x的數(shù)據(jù)長度，可以通過對x進行補零或截取來獲得。(3)Y=fft(x，[]，dim)或Y＝fft(x，N，dim)參數(shù)說明在參數(shù)dim指定的維上進行離散傅立葉變換；當(dāng)x為矩陣時，dim用來指定變換的實施方向:dim=1，表明變換按列進行，dim=2表明變換按行進行。ifft的參數(shù)應(yīng)用與fft完全相同。例4-1fft(x)函數(shù)的應(yīng)用。X=[1423]Y=fft(x)運行結(jié)果：Y=10.0000-1.0000-1.0000i-4.0000-1.0000+1.0000i例4-2ifft(x)函數(shù)的應(yīng)用。Y=[10.0000-1.0000-1.0000i-4.0000-1.0000+1.0000i]X=ifft(Y)運行結(jié)果：X=14232．fft2和ifft2函數(shù)調(diào)用方式(1)Y＝fft2(x)參數(shù)說明如果x是向量，則與一維快速傅里葉變換fft相同。如果x是矩陣，則計算該矩陣的二維快速傅里葉變換，數(shù)據(jù)二維傅里葉變換fft2(x)相當(dāng)于fft(fft(x))，即先對x的列做一維傅里葉變換，然后對變換結(jié)果的行做傅里葉變換。例4-3fft2的應(yīng)用。X=[2375;2465;2456;1234];Y=fft2(X);運行結(jié)果：Y=61.0000-14.0000+7.0000i-5.0000-14.0000-7.0000i0-7.0000i-3.0000+2.0000i4.0000-1.0000i-1.0000+2.0000i7.0000-2.0000+1.0000i1.0000-2.0000-1.0000i0+7.0000i-1.0000-2.0000i4.0000+1.0000i-3.0000-2.0000i例4-4ifft2的應(yīng)用。對例4-3中的結(jié)果Y進行反變換。X=ifft2(Y)運行結(jié)果：X=23752465245612344.8基于DFT和FFT的頻譜分析

4.8.1頻譜分析的概念

頻譜分析就是計算信號各個頻率分量的幅值、相位和功率，以獲得信號的頻率結(jié)構(gòu)以及各諧波和相位的信息。在實際工程中，大部分的信號都是模擬信號。對模擬信號x(t)進行頻譜分析，就是計算信號的傅里葉變換。由于模擬信號及其傅里葉變換都是連續(xù)函數(shù)，不能用計算機實現(xiàn)，因此通常都需要通過時域采樣將模擬信號x(t)轉(zhuǎn)化為離散時間信號x(n)，然后利用DFT和FFT進行頻譜分析。各頻率分量幅值、相位以及功率的計算如下：4.8.2頻譜分析的參數(shù)選擇

在頻譜分析中，參數(shù)的選擇至關(guān)重要，這將影響頻譜分析的準確性與精度。這些重要參數(shù)包括采樣頻率fs、頻率分辨率Δf、最小的數(shù)據(jù)記錄時間Tpmin和頻域采樣點數(shù)N。為保證采樣后的信號不發(fā)生頻率混疊，采樣頻率的選擇必須遵循奈奎斯特采樣定理，即在實際頻譜分析中，采樣頻率的選擇通常為信號最高頻率的3~4倍。頻率分辨率Δf頻率分辨率（frequencyresolution）是頻率域的采樣間隔，即用DFT進行頻譜分析時，能夠分辨的兩個頻率分量之間的最小間隔。頻率域的采樣間隔為，因此頻率分辨率Δf為：3.最小的數(shù)據(jù)記錄時間Tpmin4.頻域采樣點數(shù)N在實際頻譜分析中，頻域采樣點數(shù)就是進行DFT的點數(shù)，通常選擇為2的整數(shù)次冪。例4-5對實信號進行頻譜分析，信號最高頻率為2.5kHz，要求頻率分辨率小于10Hz，試確定以下參數(shù)：最短信號記錄時間；最大采樣周期；最少采樣點數(shù)。例4-6已知信號x1(t)=sin2(30)t，x2(t)=sin2(120)t，x(t)=x1(t)+x2(t)，采樣頻率fs＝1000Hz，試利用MATLAB對信號x(t)進行頻譜分析。解：MATLAB程序如下：t=0:0.001:0.6;x1=sin(2*pi*30*t);figure;subplot(3,1,1);plot(x1);

x2=sin(2*pi*120*t);subplot(3,1,2);plot(x2);

x=x1+x2;subplot(3,1,3);plot(x);

Y=fft(x);Pyy=Y.*conj(Y)/512;f=1000*(0:256)/512;figure;plot(f,Pyy(1:257))title('Frequencycontentofy')xlabel('frequency(Hz)')結(jié)果如圖所示。例4-7假設(shè)信號含有3種頻率成分，f1=20Hz，f2=25Hz，f3=35Hz，采樣頻率為100Hz，試用頻譜分析的方法鑒別出每個頻率成分。解：MATLAB程序如下：t=0:0.01:1;x=sin(2*pi*20*t)+sin(2*pi*25*t)+sin(2*pi*35*t)

;Subplot(2，1，1)

;plot(x)

;Y=fft(x，256);f=100*(0:(128))/256;Subplot(2，1，2)

;