《數(shù)字信號處理實驗》課件第14章

1.1市場與市場營銷1.2我國汽車市場的發(fā)展與現(xiàn)狀復(fù)習(xí)思考題實驗14快速傅里葉變換(FFT)一、實驗?zāi)康?/p>

(1)加深對快速傅里葉變換(FFT)基本理論的理解。

(2)了解使用快速傅里葉變換(FFT)計算有限長序列和無限長序列信號頻譜的方法。

(3)掌握用MATLAB語言進(jìn)行快速傅里葉變換時常用的子函數(shù)。二、實驗涉及的MATLAB子函數(shù)

1.fft

功能：一維快速傅里葉變換(FFT)。

調(diào)用格式：

y＝fft(x)；利用FFT算法計算矢量x的離散傅里葉變換，當(dāng)x為矩陣時，y為矩陣x每一列的FFT。當(dāng)x的長度為2的冪次方時，則fft函數(shù)采用基2的FFT算法，否則采用稍慢的混合基算法。

y＝fft(x，n)；采用n點FFT。當(dāng)x的長度小于n時，fft函數(shù)在x的尾部補零，以構(gòu)成n點數(shù)據(jù)；當(dāng)x的長度大于n時，fft函數(shù)會截斷序列x。當(dāng)x為矩陣時，fft函數(shù)按類似的方式處理列長度。

2.ifft

功能：一維快速傅里葉逆變換(IFFT)。

調(diào)用格式：

y＝ifft(x)；用于計算矢量x的IFFT。當(dāng)x為矩陣時，計算所得的y為矩陣x中每一列的IFFT。

y＝ifft(x，n)；采用n點IFFT。當(dāng)length(x)<n時，在x中補零；當(dāng)length(x)>n時，將x截斷，使length(x)＝n。

3.fftshift

功能：對fft的輸出進(jìn)行重新排列，將零頻分量移到頻譜的中心。

調(diào)用格式：

y＝fftshift(x)；對fft的輸出進(jìn)行重新排列，將零頻分量移到頻譜的中心。當(dāng)x為向量時，fftshift(x)直接將x中的左右兩半交換而產(chǎn)生y。

當(dāng)x為矩陣時，fftshift(x)同時將x的左右、上下進(jìn)行交換而產(chǎn)生y。三、實驗原理

1.用MATLAB提供的子函數(shù)進(jìn)行快速傅里葉變換

從理論學(xué)習(xí)可知，DFT是唯一在時域和頻域均為離散序列的變換方法，它適用于有限長序列。盡管這種變換方法是可以用于數(shù)值計算的，但如果只是簡單的按照定義進(jìn)行數(shù)據(jù)處理，當(dāng)序列長度很大時，則將占用很大的內(nèi)存空間，運算時間將很長?？焖俑道锶~變換是用于DFT運算的高效運算方法的統(tǒng)稱，F(xiàn)FT只是其中的一種。FFT主要有時域抽取算法和頻域抽取算法，基本思想是將一個長度為N的序列分解成多個短序列，如基2算法、基4算法等，大大縮短了運算的時間。

MATLAB中提供了進(jìn)行快速傅里葉變換(FFT)的子函數(shù)，用fft計算DFT，用ifft計算IDFT。

例14-1

已知一個長度為8點的時域離散信號，n1＝0，n2＝7，在n0＝4前為0，n0以后為1。對其進(jìn)行FFT變換，作時域信號及DFT、IDFT的圖形。

解程序如下：

n1＝0；n2＝7；n0＝4；

n＝n1：n2；N＝length(n)；

xn＝［(n－n0)>＝0］；%建立時域信號

subplot(2，2，1)；stem(n，xn)；

title(￠x(n)￠)；

k＝0：N－1；

Xk＝fft(xn，N)；%用FFT計算信號的DFT

subplot(2，1，2)；stem(k，abs(Xk))；

title(￠Xk＝DFT(x(n))￠)；

xn1＝ifft(Xk，N)；%用IFFT計算信號的IDFT

subplot(2，2，2)；stem(n，xn1)；

title(￠x(n)＝IDFT(Xk)￠)；

運行結(jié)果如圖14-1所示。

圖14-1

例14-1用FFT求有限長序列的傅里葉變換

例14-2

將例13-5已知的兩個時域周期序列分別取主值，得到x1＝［1，1，1，0，0，0］，x2＝［0，1，2，3，0，0］，求時域循環(huán)卷積y(n)并用圖形表示。

解本例將例13-5使用DFT處理的計算，改為用FFT和IFFT進(jìn)行循環(huán)卷積。

程序如下(作圖程序部分省略)：

xn1＝［0，1，2，3，0，0］；%建立x1(n)序列

xn2＝［1，1，1，0，0，0］； %建立x2(n)序列

N＝length(xn1)；

n＝0：N－1；k＝0：N－1；

Xk1＝fft(xn1，N)；%由x1(n)的FFT求X1(k)

Xk2＝fft(xn2，N)；%由x2(n)的FFT求

Yk＝Xk1.*Xk2；%Y(k)＝X1(k)X2(k)

yn＝ifft(Yk，N)；%由Y(k)的IFFT求y(n)

yn＝abs(yn)

運行結(jié)果如圖13-5所示，與例13-5用DFT計算的結(jié)果一致。

2.用FFT計算有限長序列的頻譜

1)基本概念

一個序號從n1到n2的時域有限長序列x(n)，它的頻譜X(ejw)定義為它的離散傅里葉變換，且在奈奎斯特(Nyquist)頻率范圍內(nèi)有界并連續(xù)。序列的長度為N，則N＝n2－n1＋1。計算x(n)的離散傅里葉變換(DFT)得到的是X(ejw)的N個樣本點X(ejwk)。其中數(shù)字頻率為

式中：dw為數(shù)字頻率的分辨率；k取對應(yīng)－(N－1)/2到(N－1)/2區(qū)間的整數(shù)。

在實際使用中，往往要求計算出信號以模擬頻率為橫坐標(biāo)的頻譜，此時對應(yīng)的模擬頻率為

式中：D為模擬頻率的分辨率或頻率間隔；Ts為采樣信號的周期，Ts＝1/Fs；定義信號時域長度L＝NTs。

在使用FFT進(jìn)行DFT的高效運算時，一般不直接用n從n1到n2的x(n)，而是取 的主值區(qū)間(n＝0，1，…，N－1)的數(shù)據(jù)，經(jīng)FFT將產(chǎn)生N個數(shù)據(jù)，定位在k＝0，1，…，N－1的數(shù)字頻率點上，即對應(yīng)［0，2p］。如果要顯示［－p，p］范圍的頻譜，則可以使用fftshift(X)進(jìn)行位移。

2)頻譜的顯示及分辨率問題

例14-3

已知有限長序列x(n)＝［1，2，3，2，1］，其采樣頻率Fs＝10Hz。請使用FFT計算其頻譜。

解

MATLAB程序如下：

Fs＝10；

xn＝［1，2，3，2，1］；N＝length(xn)；

D＝2*pi*Fs/N；%計算模擬頻率分辨率

k＝floor(－(N－1)/2：(N－1)/2)；%頻率顯示范圍對應(yīng) ［－p，p］

X＝fftshift(fft(xn，N))；%作FFT運算且移位p

subplot(1，2，1)；plot(k*D，abs(X)，￠o：￠)；%橫軸化 成模擬頻率作幅度譜

title(￠幅度頻譜￠)；xlabel(￠rad/s￠)；

subplot(1，2，2)；plot(k*D，angle(X)，￠o：￠)；%橫軸 化成模擬頻率作相位譜

title(￠相位頻譜￠)；xlabel(￠rad/s￠)；

程序運行結(jié)果：

absX＝

0.38202.61809.00002.61800.3820

angleX＝

－1.25662.5133[KG*4/5]0－2.51331.2566

運行結(jié)果如圖14-2所示。

圖14-2例14-3有限長序列的頻譜由圖14-2可知，當(dāng)有限長序列的長度N＝5時，頻譜的頻率樣本點數(shù)也為5，如圖上用“”表示的點位。頻率點之間的間距非常大，即分辨率很低。即使使用了plot命令的插值功能，顯示出的曲線仍是斷斷續(xù)續(xù)的，與真實曲線有較大的誤差。

改變分辨率的基本方法是給輸入序列補零，即增加頻譜的密度。注意，這種方法只是改善了圖形的視在分辨率，并不增加頻譜的細(xì)節(jié)信息。

將上述有限長序列x(n)＝［1，2，3，2，1］末尾補0到N＝1000點，將程序改為：

Fs＝10；N＝1000；

xn＝［1，2，3，2，1］；Nx＝length(xn)；

xn＝［1，2，3，2，1，zeros(1，N－Nx－1)］；

D＝2*pi*Fs/N；%計算模擬頻率分辨率

k＝floor(－(N－1)/2：(N－1)/2)；%頻率顯示范圍對應(yīng) ［－p，p］

X＝fftshift(fft(xn，N))；%作FFT運算且移位p

subplot(1，2，1)；plot(k*D，abs(X))；%橫軸化成模擬頻率作幅度譜

title(￠幅度頻譜￠)；xlabel(￠rad/s￠)；

subplot(1，2，2)；plot(k*D，angle(X))；%橫軸化成模擬頻率作相位譜

title(￠相位頻譜￠)；xlabel(￠rad/s￠)；

此時程序執(zhí)行的結(jié)果如圖14-3所示。由圖可以看出，圖形的分辨率提高，曲線幾乎是連續(xù)的頻譜了。

圖14-3將例14-2有限長序列末尾補0到N＝1000時的頻譜

3)實偶序列如何補0

例14-4

已知一個矩形窗函數(shù)序列為

采樣周期Ts＝0.5s，要求用FFT求其頻譜。

解由于該序列是一個實的偶序列，因而補0時需要仔細(xì)分析。假定按N＝32補0，則主值區(qū)域在n＝0～31，F(xiàn)FT的輸入應(yīng)為即原來n=[－5：－1]的前五個點移到n＝[27：31]中去了。

下面考慮分別用N=32，64，512，觀察不同N值代入對頻譜的影響。

程序如下，

Ts=0.5;C=[32,64,512];

%輸入不同的N值

forr＝0：2；

N＝C(r＋1)；

xn＝［ones(1，6)，zeros(1，N－11)，ones(1，5)］； %建立x(n)

D＝2*pi/(N*Ts)；

k＝floor(－(N－1)/2：(N－1)/2)；

X＝fftshift(fft(xn，N))；

subplot(3，2，2*r＋1)；plot(k*D，abs(X))；%幅度頻譜

subplot(3，2，2*r＋2)；stairs(k*D，angle(X))；%相位 頻譜

end

注意：此處相位頻譜使用了stairs，因為該相位頻譜變化率比較陡峭。

程序執(zhí)行結(jié)果如圖14-4所示。

圖14-4將例14-4有限長序列補0到N＝32、64、512時的頻譜如果將x(n)的輸入寫成

xn＝［ones(1，11)，zeros(1，N－11)］；%建立x(n－5)

相當(dāng)于起點不是取自n＝0而是n＝－5，計算的是x(n－5)的頻譜。幅度頻譜不受影響，相位頻譜引入一個線性相位－5w，如圖14-5所示。

圖14-5將有限長位移序列x(n－5)補0到N＝512時的頻譜

3.用FFT計算無限長序列的頻譜

用FFT進(jìn)行無限長序列的頻譜計算，首先要將無限長序列截斷成一個有限長序列。序列長度的取值對頻譜有較大的影響，帶來的問題是引起頻譜的泄漏和波動。

例14-5

已知一個無限長序列為

采樣頻率Fs＝20Hz，要求用FFT求其頻譜。

解

MATLAB程序如下：

Fs＝20；C＝［8，16，128］；%輸入不同的N值

forr＝0：2；

N＝C(r＋1)；

n＝0：N－1；

xn＝exp(－0.5*n)；%建立x(n)

D＝2*pi*Fs/N；

k＝floor(－(N－1)/2：(N－1)/2)；

X＝fftshift(fft(xn，N))；

subplot(3，2，2*r＋1)；plot(k*D，abs(X))；

axis(［－80，80，0，3］)；

subplot(3，2，2*r＋2)；stairs(k*D，angle(X))；

axis(［－80，80，－1，1］)；

end

運行結(jié)果如圖14-6所示。

圖14-6將無限長序列截斷為N＝8，16，128時的頻譜由圖14-6可見，N值取得越大，即序列保留得越長，曲線精度越高。

例14-6

用FFT計算下列連續(xù)時間信號的頻譜，并觀察選擇不同的Ts和N值對頻譜特性的影響。

xa(t)＝e－0.01t(sin2t＋sin2.1t＋sin2.2t)t≥0

解該題選擇了三個非常接近的正弦信號，為了將各頻率成分區(qū)分出來，在滿足奈奎斯特定理的條件下確定采樣周期，選擇三組數(shù)據(jù)，分別是Ts＝0.5s、0.25s和0.125s；再確定N值，分別選擇N＝256和2048。觀察不同Ts和N的組合對頻譜的影響。程序如下：

T0＝［0.5，0.25，0.125，0.125］；%輸入不同的Ts值

N0＝［256，256，256，2048］；%輸入不同的N值

forr＝1：4；

Ts＝T0(r)；N＝N0(r)；%賦Ts和N值

n＝0：N－1；

D＝2*pi/(Ts*N)；%計算模擬頻率分辨率

xa＝[ZK(]exp(－0.01*n*Ts).*(sin(2*n*Ts)＋sin(2.1*n*Ts)＋sin(2.2*n*Ts))；

k＝floor(－(N－1)/2：(N－1)/2)；

Xa＝Ts*fftshift(fft(xa，N))；

［r，Xa(1)］%輸出Xa(1)的數(shù)值，供誤差計算用

subplot(2，2，r)；plot(k*D，abs(Xa)，￠k￠)；

axis(［1，3，1.1*min(abs(Xa))，1.1*max(abs(Xa))］)；

end

運行結(jié)果如圖14-7所示。

圖14-7用FFT計算三個很靠近的諧波分量的頻譜圖由圖14-7可以得出以下結(jié)論：

(1)N同樣取256(如前三個圖形)，當(dāng)Ts越大時，時域信號的長度L＝NTs保留得越長，則分辨率越高，頻譜特性誤差