實(shí)驗(yàn)用MATLAB計(jì)算傅里葉變換

1、實(shí)驗(yàn)二用MATLAB計(jì)算傅立葉變換(2課時(shí))、實(shí)驗(yàn)?zāi)康?、掌握用MATLAB計(jì)算DTFT及系統(tǒng)頻率響應(yīng)的方法。2、掌握用 MATLAB 計(jì)算DFT和IDFT的方法。3、掌握用DFT計(jì)算圓周卷積和線(xiàn)性卷積的方法。、實(shí)驗(yàn)設(shè)備計(jì)算機(jī)一臺(tái)，裝有 MATLAB軟件。三、實(shí)驗(yàn)原理和基本操作1. 用 MATLAB 計(jì)算 DTFTDTFT )定義為:(1)對(duì)于序列x (n),其離散時(shí)間傅立葉變換(O0X(j )二 ' x(n)e-j nn :序列的傅立葉變換(DTFT )在頻域是連續(xù)的，并且以 =2 n為周期。因此只需要知道 X(ejw ) 的一個(gè)周期，即 =0 , 2 n ，或-n , n 。就可以

2、分析序列的頻譜。用MATLAB計(jì)算DTFT ,必須在-n < n范圍內(nèi)，把，用很密的、長(zhǎng)度很長(zhǎng)的向量來(lái)近 似，該向量中各個(gè)值可用下式表示：2兀 w=k*dw=k*(2)K2其中：d =稱(chēng)為頻率分辨率。它表示把數(shù)字頻率的范圍2 n均分成K份后，每一份K的大小，k是表示頻率序數(shù)的整數(shù)向量，簡(jiǎn)稱(chēng)為頻序向量，它的取值可以有幾種方法：通常在 DTFT中，頻率取-n W< ji的范圍，當(dāng)K為偶數(shù)時(shí)，取, K K +4K彳 k ,1,T,0,1,12 2 2如果K為奇數(shù)，則取Kkk0.5/ 1,0,1/ ,0.522可以為奇偶兩種情況綜合出一個(gè)共同的確定頻序向量k的公式；K -1K -1.2 一

3、:一 2 一k=上式中 表示向下取整。在MATLAB中的向下取整函數(shù)為 floor , floor (x)的作用是把x向下(向方向)取整, 所以與(3)式等價(jià)的MATLAB語(yǔ)句為(4)k = floor (-K 0.5):( -0.5)2 2給定了輸入序列(包括序列x及其位置向量n),又設(shè)定了頻率分辨率 d及頻序向量k，則DTFT的計(jì)算式(1)可以用一個(gè)向量與矩陣相乘的運(yùn)算來(lái)實(shí)現(xiàn)。e_j 性2n1.屮 Kn1X( dX( 2), ,X( k) =X(nJ, ,x(nN)e1n2 ej 2匚eKn如果頻率向量表為 =L ! / 2,，.=k*d,而序列的位置向量為 n X= n1 ： nN,則(

4、5)式中的矩陣的指數(shù)部分可以寫(xiě)成 -j*n T* ，用MATLAB語(yǔ)句表達(dá)時(shí)，把代以w，轉(zhuǎn)置符 號(hào)nT換成MATLAB中的相應(yīng)符號(hào) n'，則求DTFT的程序可以寫(xiě)成：X =x*exp(-j *dw* n'*k)求有限序列x=2 , -1, 1, 1的DTFT，其位置向量為 nx=0 : 3。假如取64個(gè)頻點(diǎn)，畫(huà)出它在-冗_(dá) _ n范圍內(nèi)的幅頻和相頻特性。程序如下:x=2,-1,1,1; nx=0:3;K=64;dw=2*pi/K;%設(shè)定頻序向量k=floor(-K/2+0.5):(K/2-0.5);%w=li nspace(-8,8,1000);X=x*exp(-j*dw* n

5、x'*k)用(1)式計(jì)算DTFTsubplot(2,1,1); plot(k*dw,abs(X),% 畫(huà)幅頻曲線(xiàn)subplot(2,1,2) ;plot(k*dw,angle(X),% 畫(huà)相頻曲線(xiàn)2 .用MATLA B計(jì)算系統(tǒng)頻率響應(yīng) H (e")MATLA B中用于求系統(tǒng)頻率響應(yīng)的函數(shù)是freqz,其功能是：由給定的系統(tǒng)函數(shù) H(Z)的分子和分母的系數(shù)向量繪制系統(tǒng)的幅度和相位響應(yīng)。調(diào)用格式：H , w=freqz(b,a,N)，或 N 缺省H , w=freqz(b,a),此時(shí) N 取默認(rèn)值 512。H是系統(tǒng)的頻率特性，它是一個(gè)N元的復(fù)數(shù)向量。w是數(shù)字頻率向量，它把0到n均

6、分為N份，分辨率為n /N , w =0 : N-1 n /N。b和a分別為分子分母多項(xiàng)式的負(fù)幕系數(shù)向量，即多項(xiàng)式的首項(xiàng)是常數(shù)項(xiàng)，以后按e-的升幕排列，由此形成的多項(xiàng)式的系數(shù)向量。N為所選的頻率點(diǎn)數(shù)，它決定了頻率分辨率的密度。這樣求出的頻率特性是在 0W <n之間的特性，如果要讓 MATLAB計(jì)算全頻率0W <2 n 范圍的頻率特性，要增加一個(gè)輸入變?cè)獁hole ' H，w=freqz(b,a,N, 'whole '此時(shí)w的范圍是0到2n (不含2 n )。若沒(méi)有左端變量，即鍵入 freqz(b,a,N) , MATLAB將不給出數(shù)據(jù)H和w而只繪出頻率特性曲

7、 線(xiàn)。若在輸入變?cè)薪o出頻點(diǎn)向量w: freqz(b,a,w)，這種調(diào)用方法可以在自己選定的頻點(diǎn)向量w上計(jì)算頻率特性。例 2：設(shè)一個(gè) LTI 系統(tǒng)的差分方程為 y(n)-0.9y(n-1)=0.5x(n)+0.8x(n-1)求其頻率響應(yīng)。可直接寫(xiě)出HGp050.8eJ1 -0.9e 朋用下面的程序可繪出幅頻和相頻特性曲線(xiàn)。b=0.5 0.8;%分子多項(xiàng)式系數(shù)向量a=1 -0.9;%分母多項(xiàng)式系數(shù)向量H,w=freqz(b,a);%求出頻率響應(yīng)(0到pi分成500點(diǎn))subplot(2,1,1),plot(w,abs(H),grid on xlabel('頻率'),ylabel(

8、'幅度 dB')subplot(2,1,2),plot(w,angle(H),grid onxlabel('頻率'),ylabel('相角(度)')3. 用 MATLAB 計(jì)算 DFT 和 IDFT可以利用 MATLAB提供的計(jì)算快速傅立葉變換(FFT)和快速傅立葉反變換(IFFT)的函 數(shù)來(lái)計(jì)算離散傅立葉變換(DFT)和離散傅立葉反變換(IDFT )。在MATLAB信號(hào)處理工具箱中，函數(shù) fft和ifft用于快速傅立葉變換和反變換。X=fft(x)完成對(duì)序列x的L點(diǎn)DFT,其中L為序列x的長(zhǎng)度。X=fft(x , N)則指定了采用 N點(diǎn) DFT

9、。如果N>L,則程序會(huì)自動(dòng)給 x后面補(bǔ)零,使其長(zhǎng)度為N;如果N<L,則程序會(huì)自動(dòng)將x截?cái)?取 前N個(gè)數(shù)據(jù)。同樣，計(jì)算離散傅立葉反變換時(shí)，調(diào)用函數(shù) ifft方法為ifft (X )或ifft(X , N),其 輸入變?cè)囊饬x與 fft函數(shù)相同。例 3:長(zhǎng)度為 4 的有限序列：x(0)=2,x(1)= -1,x(2)=1,x(3)=1,求它的 DFT。解: x=2,-1,1,1;X=fft(x)運(yùn)行結(jié)果:4 .用MATLAB計(jì)算圓周卷積圓周卷積定理：設(shè)有限長(zhǎng)序列X,n)和X2(n)的長(zhǎng)度分別為 Ni和N? , N=maxN仆n?°x"n) 和X2(n)的N點(diǎn)DFT分

10、別為X1( k)=DFTx 1(n),X 2 (k)=DFTx 2 (n).若 X( k)= X1( k)* X 2(k)則X1 (n)與X2 (n)的N點(diǎn)圓周卷積是N Jx( n)=IDFTX(k)= ' X1(m)X2( n- m)N)RN( n)m -0可以利用此定理和函數(shù)fft、ifft計(jì)算圓周卷積例 4：設(shè) x1 (n)=1 , 2，3，x? (n)=5，4， -3，-2，計(jì)算 4 點(diǎn)圓周卷積 y(n)= x 1 (n) x 2 (n)解：寫(xiě)出MATLAB程序如下：x1=1 2 3,x2=5 4 -3 -2Xi =fft(x1 ,4); X2 =fft(x2 ,4);Y= X

11、i .* X2 ;y=ifft(Y,4)運(yùn)行結(jié)果：5利用圓周卷積計(jì)算線(xiàn)性卷積設(shè)x1 (n)為N 1點(diǎn)序列，x 2 (n)為N 2點(diǎn)序列，x1(n)和x2 (n)的線(xiàn)性卷積為Niy(n)='二 Xi(m)X2(n-m)m £y(n)%( Nj+ N 2 -1 )點(diǎn)序列。如果N 一 Ni+ N2-I，貝U N點(diǎn)圓周卷積能代表線(xiàn)性卷積。例5:設(shè)Xi(n)與X2(n)是兩個(gè)四點(diǎn)序列：xi(n)=1 , 2, 2, 1,X2(n)=1 , -1, 1, -1(1)求它們的線(xiàn)性卷積 y(n);計(jì)算圓周卷積，使得它與y(n)相等。解：線(xiàn)性卷積可以調(diào)用conv函數(shù)來(lái)求。圓周卷積可按例4的方法

12、求。程序如下：x1=1,2,2,1;x2 =1,-1,1,-1；y=c on v(x1 , x2 )X1 =fft(x1 ,7); X2 =fft(x2 ,7); Y= X1 .* X2 ;x3 =ifft(Y,7)運(yùn)行結(jié)果：6 分段卷積當(dāng)x(n)長(zhǎng)度遠(yuǎn)遠(yuǎn)大于h(n)時(shí)，要采用分段卷積的方法，以減少線(xiàn)性卷積的運(yùn)算量。分段卷積 的方法有重疊保留法和重疊相加法兩種。(1)重疊保留法設(shè)h(n)的長(zhǎng)度為M , x(n)的長(zhǎng)度為L(zhǎng)x，分段后每一段的有效數(shù)據(jù)的長(zhǎng)度為L(zhǎng)=N-M+1。編寫(xiě)出實(shí)現(xiàn)重疊保留法的函數(shù)ovrlpsav如下：fun cti on y=ovrlpsav(x,h,N)Lx=le ngth(

13、x);M=le ngth(h);M1=M-1;L=N-M1;H=fft(h,N);x=zeros(1,M1),x,zeros(1,N-1);K=floor(Lx+M1-1)/(L)+1;Y=zeros(K+1,N);for k=0:K-1xk=x(k*L+1:k*L+N);Xk=fft(xk);Y(k+1,:)=real(ifft(Xk.*H);endY=Y(:,M:N)'y=(Y(:)'(2)重疊相加法MATLAB信號(hào)處理工具箱中的函數(shù)fftfilt.m可實(shí)現(xiàn)重疊相加法。它有兩種調(diào)用格式：1) y=fftfilt(h,x)2) y=fftfilt(h,x,r)其中h代表系統(tǒng)的脈沖響應(yīng)h(n),x是輸入序列。在第一種格式中，程序自動(dòng)把輸入分成每段512個(gè)樣本。并按512點(diǎn)(如果h(n)的長(zhǎng)度比512長(zhǎng)，則按h(n)的長(zhǎng)度)FFT進(jìn)行各段的卷積運(yùn) 算。在第二種調(diào)用格式中，r是用戶(hù)指定的FFT長(zhǎng)度，而輸入x就按這個(gè)長(zhǎng)度分段。例 6：設(shè) x(n)=n+1,0 _n _9,h(n)=1,0,-1。分別用 ovrlpsav(按 N=6)和 fftfilt 函數(shù)求 x(n)與 h(n)的線(xiàn) 性卷積解：MATLAB程序?yàn)閚=0:9;x=n+1;h=1,0,-1;N=6;y1=ovrlps