信與系統(tǒng)實(shí)驗(yàn)報(bào)告 實(shí)驗(yàn)4非周期信的傅里葉變換實(shí)驗(yàn)

1、信號(hào)與系統(tǒng)實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)非周期信號(hào)的傅里葉變換實(shí)驗(yàn)四非周期信號(hào)的傅里葉變換一、實(shí)驗(yàn)?zāi)康母道锶~變換是通信系統(tǒng)、圖像處理、數(shù)字信號(hào)處理以及物理學(xué)等領(lǐng)域內(nèi)的一種重要的數(shù)學(xué)分析工具。通過(guò)傅里葉變換技術(shù)可以將時(shí)域上的波形分布變換為頻域上的分布，從而獲得信號(hào)的頻譜特性。MATLAB提供了專門(mén)的函數(shù)fft、ifft、fft2（即卩2維快速傅里葉變換）、ifft2以及fftshift用于實(shí)現(xiàn)對(duì)信號(hào)的傅里葉變換。本次實(shí)驗(yàn)的目的就是練習(xí)使用fft、fft以及fftshift函數(shù),對(duì)一些簡(jiǎn)單的信號(hào)處理問(wèn)題能夠獲取其頻譜特性(包括幅頻和相頻特性)。二、實(shí)驗(yàn)預(yù)備知識(shí)1.離散傅里葉變換(DFT)以及快速傅里葉變換(FFT)

2、簡(jiǎn)介設(shè)x(t)是給定的時(shí)域上的一個(gè)波形，則其傅里葉變換為X(f)=_Tx(t)e-j2nftdt(1)g顯然X(f)代表頻域上的一種分布(波形)，一般來(lái)說(shuō)X(f)是復(fù)數(shù)。而傅里葉逆變換定義為：x(t)二廣X(f)ej2Kftdf(2)g因此傅里葉變換將時(shí)域上的波形變換為頻域上的波形，反之，傅里葉逆變換則將頻域上的波形變換為時(shí)域上的波形。由于傅里葉變換的廣泛應(yīng)用，人們自然希望能夠使用計(jì)算機(jī)實(shí)現(xiàn)傅里葉變換，這就需要對(duì)傅里葉變換(即(1)式)做離散化處理，使之符合電腦計(jì)算的特征。另外，當(dāng)把傅里葉變換應(yīng)用于實(shí)驗(yàn)數(shù)據(jù)的分析和處理時(shí)，由于處理的對(duì)象具有離散性，因此也需要對(duì)傅里葉變換進(jìn)行離散化處理。而要想將

3、傅里葉變換離散化，首先要對(duì)時(shí)域上的波形x(t)進(jìn)行離散化處理。采用一個(gè)時(shí)域上的采樣脈沖序列:5(t-nT),n=0,1,2,，N_1；可以實(shí)現(xiàn)上述目的，如圖所示。其中N為采樣點(diǎn)數(shù)，T為采樣周期；f=1/T是采樣頻率。注意采樣時(shí)，采樣頻率f必須大于兩倍ss的信號(hào)頻率(實(shí)際是截止頻率)，才能避免混迭效應(yīng)。接下來(lái)對(duì)離散后的時(shí)域波形x(t)=x(t)S(tnT)=x(nT)的傅里葉變換X(f)進(jìn)行離散處理。與上述做法類似，采用頻域上的5脈沖序列：5(f-n/To),n=0,1,2,N-1；T=NT為總采樣時(shí)間x(t)5脈沖序列鳥(niǎo)X")5nT)可以實(shí)現(xiàn)傅里葉變換Xf)的離散化，如下圖示。不難看

4、出，離散后的傅里葉變換其頻率間隔(頻率軸上離散點(diǎn)的間隔，即頻域分辨率)Af=T(3)因此要增加分辨率須增加采樣點(diǎn)數(shù)目N。頻域上每個(gè)離散點(diǎn)對(duì)應(yīng)的頻率為:nnf=n-TNTN0n=0,1,2,.,N-1(4)顯然n=0的點(diǎn)對(duì)應(yīng)于直流成分。經(jīng)過(guò)以上離散化處理之后，連續(xù)積分的傅里葉變換(1)式轉(zhuǎn)變?yōu)槿缦码x散形式：7-1ifiiini8(f-n/T0)Hffllllk/F1W.fX(f):w.fX(f)混迭X(f)=x(t)e-j2兀nk/N,n=0,1,2,.,N-1(5)其中tk=kTnkk=0(k=0,1,2,N-l)代表采樣點(diǎn)時(shí)刻。X(盒)一般是復(fù)數(shù)，因此離散傅里葉變換(DFT)后變成一個(gè)N點(diǎn)(

5、采樣點(diǎn)數(shù))的復(fù)數(shù)序列。X(盒)絕對(duì)值代表振幅，其幅角代表相位，因此由(5)式可以給出DFT的振幅頻譜和相位頻譜。(5)式通常又簡(jiǎn)寫(xiě)成如下形式：X(n)=Ex(k)Wnk,n=0,1,2,.,N-1(6)Nk=0其中W=e-j2兀/n，x是采樣點(diǎn)數(shù)據(jù)，它是一個(gè)N個(gè)點(diǎn)的向量，DFTN的結(jié)果X是N個(gè)點(diǎn)的復(fù)數(shù)向量。(5)式或(6)式就是對(duì)傅里葉變換進(jìn)行數(shù)值計(jì)算的基礎(chǔ)。一般采樣點(diǎn)數(shù)N越大，DFT的結(jié)果越接近真實(shí)的情況，但是當(dāng)N較大時(shí)，(6)式的計(jì)算量很大，因?yàn)槭褂糜?jì)算機(jī)求解(6)式時(shí)，總共要執(zhí)行N2次復(fù)數(shù)乘法和Nx(N-l)次復(fù)數(shù)加法。所以直接用DFT算法(即(5)式)進(jìn)行譜分析和信號(hào)的實(shí)時(shí)處理是不切實(shí)

6、際的。為了減輕計(jì)算的壓力，人們提出了一種所謂快速傅里葉變換(FFT)的思想：取N=2m,首先將N個(gè)點(diǎn)的采樣數(shù)據(jù)x二X,xx分成兩個(gè)N/201N1點(diǎn)的序列：xx,x，,x(偶數(shù)序列)1 02N2xx,x，,x(奇數(shù)序列)2 13N1這樣處理的好處是可以把(6)式分解為兩個(gè)N/2點(diǎn)的DFT,使計(jì)算量降下來(lái)。接下來(lái)再將N/2點(diǎn)的序列仿照上述做法進(jìn)一步分裂成2個(gè)N/4點(diǎn)的序列X3和X4,另一序列X2亦做如此處理，分裂成2個(gè)N/4點(diǎn)的序列X5和工6。這3 4256樣兩個(gè)N/2點(diǎn)的序列分成了更短的4個(gè)N/4點(diǎn)的序列，依次類推，最后的結(jié)果是將一個(gè)N點(diǎn)的序列x裂成了N個(gè)點(diǎn)的單點(diǎn)序列：x0,X,x2,xN-1o

7、這樣做可以將DFT的運(yùn)算效率提高12個(gè)數(shù)量級(jí)，為數(shù)字信號(hào)處理技術(shù)應(yīng)用于各種信號(hào)的實(shí)時(shí)處理創(chuàng)造了條件，從而推動(dòng)數(shù)字處理技術(shù)的發(fā)展。由此可見(jiàn)FFT的思想實(shí)質(zhì)是不斷地把長(zhǎng)序列的DFT計(jì)算分解成若干短序列的DFT,并利用旋轉(zhuǎn)因子(即WN)的周期性和對(duì)稱性來(lái)減少DFT的運(yùn)算次數(shù)。所以FFT就是DFT的快速算法。有關(guān)FFT算法的詳細(xì)介紹和理論推導(dǎo)參見(jiàn)有關(guān)的書(shū)籍，這里不做進(jìn)一步介紹。2.FFT的MATLAB實(shí)現(xiàn)為了實(shí)現(xiàn)快速傅里葉變換，MATLAB提供了fft、ifft、fft2、ifft2以及fftshift函數(shù)，分別用于一維和二維離散傅里葉變換(DFT)及其逆變換。借助這些函數(shù)可以完成很多信號(hào)處理任務(wù)。考

8、慮到信號(hào)處理包含的領(lǐng)域很廣泛,這里只介紹一維傅里葉變換及其逆變換函數(shù)。(1) fft函數(shù)該函數(shù)使用了快速算法來(lái)實(shí)現(xiàn)時(shí)域信號(hào)的離散傅里葉變換。常用的格式：Y=fft(x)Y=fft(x,m)Y返回值(復(fù)數(shù))，返回m點(diǎn)的DFT序列，即(6)式左邊的X；m計(jì)算時(shí)使用的數(shù)據(jù)點(diǎn)數(shù)(樣本數(shù))；x時(shí)域信號(hào)x(t)在采樣點(diǎn)tk處的值，即(6)式右邊的X；若實(shí)際采樣點(diǎn)數(shù)目為N(m和N都須是2的幕次)，則x為N個(gè)元素(即長(zhǎng)度N)的向量；若向量x的長(zhǎng)度小于m,那么計(jì)算時(shí)將自動(dòng)在x序列的后面補(bǔ)0；若x的長(zhǎng)度大于m,則x自動(dòng)截?cái)?，使之長(zhǎng)度為m。對(duì)信號(hào)進(jìn)行頻譜分析時(shí)，數(shù)據(jù)樣本應(yīng)有足夠的長(zhǎng)度，一般FFT程序中所用數(shù)據(jù)點(diǎn)數(shù)(

9、m)最好與原信號(hào)含有的數(shù)據(jù)點(diǎn)數(shù)(即輸入的樣本數(shù)N)相同，這樣的頻譜圖具有較高的質(zhì)量，可減小因補(bǔ)零或截?cái)喽a(chǎn)生的影響。兩點(diǎn)說(shuō)明： 關(guān)于FFT振幅頻譜和相位頻譜的計(jì)算由于傅里葉變換的結(jié)果一般是復(fù)數(shù)，所以 對(duì)fft的結(jié)果取絕對(duì)值abs()可以得到振幅，即Amplitude=abs(Y)需要注意的是這樣得到的幅值實(shí)際并非真正的信號(hào)振幅，因其值與FFT使用的數(shù)據(jù)點(diǎn)數(shù)N有關(guān)，但不影響分析結(jié)果，在IFFT(逆變換)時(shí)已經(jīng)做了處理。要得到真實(shí)的振幅值的大小，只要將上述結(jié)果除以N/2即可。 對(duì)fft的結(jié)果使用函數(shù)angle()可以得到相位的結(jié)果。但是使用angle函數(shù)計(jì)算復(fù)數(shù)的相角時(shí)，系統(tǒng)規(guī)定一、二象限的角為0

10、冗；三、四象限的角為冗0。因此若一個(gè)角度本來(lái)應(yīng)該從0變到2冗，但計(jì)算得到的結(jié)果卻是0冗，再由-冗0,在冗處發(fā)生跳變，跳變幅度為2冗，這就叫相位的卷繞。這種相位的卷繞會(huì)使得相頻圖不連續(xù)，呈現(xiàn)鋸齒狀，為了平滑相頻圖，通常要再使用unwrap()函數(shù)進(jìn)行相位的解卷繞。因此FFT的相位頻譜圖應(yīng)該如下實(shí)現(xiàn)Phase=unwrap(angle(Y) FFT的振幅頻譜具有對(duì)稱性如下圖所示。FFT：對(duì)稱軸(Nyquist頻率)1負(fù)頻部分Tf軸N/2fsT20N-1N丄fTTs00nnf,(n0,1,2,N1)為頻率軸上的頻率點(diǎn)。nTNTNs因此用FFT對(duì)信號(hào)做譜分析，只需考察0Nyquist頻率范圍內(nèi)(共N/

11、2+1個(gè)頻率點(diǎn))的幅頻特性。(2) fftshift函數(shù)其作用是將零頻點(diǎn)移到頻譜的中間(即Nyquist頻率處)，使用格式：Y=fftshift(X)X是向量，該命令將零頻點(diǎn)移動(dòng)到頻譜X的中間，并交換頻譜X的左右兩半。將零頻點(diǎn)放到頻譜的中間對(duì)于觀察傅立葉變換是有用的。例1:對(duì)時(shí)域信號(hào)x(t)二0.5sin(2兀15t)+2sin(2兀40t)進(jìn)行頻譜分析。fs=100;%采樣頻率2倍的信號(hào)頻率N=256;%采樣點(diǎn)數(shù)目(=2的冪次)n=0:N-1;%構(gòu)造采樣點(diǎn)序列t=n/fs;%得到采樣時(shí)間序列，t=nT=n/fsx=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%產(chǎn)生

12、時(shí)域信號(hào)的樣本值，向量Y=fft(x,N);%N點(diǎn)的DFT計(jì)算mag=abs(Y);%FFT的振幅phase=unwrap(angle(Y);%FFT的相位%1.以下繪制物理頻譜圖(即正頻部分)fn=(0:N/2)*fs/N;%頻率軸上的離散頻率點(diǎn)，起始于0頻(對(duì)應(yīng)直流成分)，終%于Nyquist頻率fs/2,共N/2+1個(gè)頻率點(diǎn)八、subplot(2,2,1)%將圖形窗口分割為2X2的子窗口，并指定第1個(gè)子窗口為繪圖區(qū)plot(fn,mag(1:N/2+1)%取出前N/2+1個(gè)振幅作圖，即正頻率分量xlabel('頻率/Hz');ylabel('振幅');ti

13、tle('圖1：物理(正頻)幅頻圖');gridon%加網(wǎng)格線%2.以下繪制全頻率的幅頻圖fn1=(0:N-1)*fs/N;subplot(2,2,2)%指定第2個(gè)子窗口為繪圖區(qū)plot(fn1,mag);xlabel('頻率/Hz');ylabel('振幅');title('圖2：全頻率的幅頻圖');gridon%3.以下繪制正頻部分的相頻圖subplot(2,2,3)%指定第3個(gè)子窗口為繪圖區(qū)plot(fn,phase(1:N/2+1);xlabel('頻率/Hz');ylabel('相位');

14、title('圖3：相頻圖');grid%4.以下移動(dòng)零頻點(diǎn)Yl=fftshift(Y);%fftshift移動(dòng)頻率零點(diǎn)，并將Y的左右兩部分交換mag1=abs(Y1);%重新計(jì)算振幅fn2=fn1-fs/2;%零點(diǎn)移動(dòng)到fs/2處，故需重新標(biāo)記頻率軸subplot(2,2,4);%指定第4個(gè)子窗口為繪圖區(qū)，最終4幅圖繪制在一張圖上了plot(fn2,mag1);xlabel('頻率/Hz');ylabel('振幅');title('圖4:fftshift后的幅頻圖');grid運(yùn)行結(jié)果如下：圖說(shuō)明：1是物理譜圖（正頻部分），從中

15、看到，該信號(hào)包含兩個(gè)頻率15Hz和40Hz。由于使用的采樣頻率fs=100Hz，所以Nyquist頻率為50Hz，在圖2中明顯能看到整個(gè)頻譜圖關(guān)于Nyquist頻率對(duì)稱，不過(guò)Nyquist頻率右邊的譜圖實(shí)際上是負(fù)頻部分，沒(méi)有意義。圖4是fftshift之后的幅頻圖，由于它是圖2結(jié)果的左右交換，因此圖2右邊變成了負(fù)頻。另外，圖中的振幅不是真實(shí)的信號(hào)振幅，從信號(hào)x(t)的表達(dá)式我們知道15Hz和40Hz這兩種頻率成分的振幅分別是0.5和2。要得到真實(shí)的振幅,只需要將程序中的mag除以N/2即可。(3) ifft函數(shù)執(zhí)行離散傅里葉變換的逆變換，格式x=ifft(Y)或者x=ifft(Y,m)Y是FF

16、T的輸出結(jié)果，返回值x是時(shí)域上的結(jié)果，m仍然是計(jì)算使用的數(shù)據(jù)點(diǎn)數(shù)。在上例中若程序末尾使用：xx=ifft(YN),則得到采樣時(shí)刻點(diǎn)上，信號(hào)x(t)的樣本值。三、實(shí)驗(yàn)內(nèi)容及要求實(shí)驗(yàn)項(xiàng)目：給定采樣頻率51.2Hz及采樣點(diǎn)數(shù)N=512,計(jì)算矩形函數(shù)x(t)=|01的振幅頻譜并與理論計(jì)算結(jié)果對(duì)比。A. 顯然該信號(hào)x(t)是無(wú)限長(zhǎng)的非周期信號(hào)，因此做FFT計(jì)算時(shí)必須先將信號(hào)x(t)截?cái)酁橛邢揲L(zhǎng)度。令采樣頻率為厶，采樣點(diǎn)數(shù)目N則截?cái)嚅L(zhǎng)度是：T0=N/fs(即總的采樣時(shí)間)因此截?cái)嚅L(zhǎng)度和采樣點(diǎn)數(shù)目N成正比。對(duì)于無(wú)限長(zhǎng)的非周期信號(hào)，截?cái)嚅L(zhǎng)度應(yīng)盡可能的大，以接近實(shí)際信號(hào)，避免結(jié)果失真；如果是周期信號(hào)，則要求截?cái)嚅L(zhǎng)度為信號(hào)周期的整數(shù)倍，以免出現(xiàn)頻譜的“泄漏”。若給定采樣點(diǎn)數(shù)N=512，則時(shí)間采樣序列可用向量t表示：t=(O:Nl)/fs，矩形函數(shù)x(t)的樣本值可以使用MATLAB提供的符號(hào)函數(shù)sign(請(qǐng)使用helpsign命令查詢sign函數(shù)的定義)來(lái)表示：x=0.50.5*sign(t1)，然后使用fft命令即可獲得DFT計(jì)算結(jié)果。B. 根據(jù)(1)式不難算出上述信號(hào)x(t)的傅里葉變換的理論結(jié)果(精確值)：X(f)=sin2兀f+j(cos2兀f-1)2兀f(7)其振幅為：|X(f)*2-2cos(2兀f)/2冗f(8)然后