《數(shù)字信號處理》實驗指導(dǎo)書--應(yīng)用試驗

1、數(shù)字信號處理實驗指導(dǎo)書電信分院電子教研室42第1章 系統(tǒng)響應(yīng)及系統(tǒng)穩(wěn)定性1.1 實驗?zāi)康膌 學(xué)會運用MATLAB求解離散時間系統(tǒng)的零狀態(tài)響應(yīng)；l 學(xué)會運用MATLAB求解離散時間系統(tǒng)的單位取樣響應(yīng)；l 學(xué)會運用MATLAB求解離散時間系統(tǒng)的卷積和。1.2 實驗原理及實例分析1.2.1 離散時間系統(tǒng)的響應(yīng)離散時間LTI系統(tǒng)可用線性常系數(shù)差分方程來描述，即 （1-1）其中，（，1，N）和（，1，M）為實常數(shù)。MATLAB中函數(shù)filter可對式（13-1）的差分方程在指定時間范圍內(nèi)的輸入序列所產(chǎn)生的響應(yīng)進(jìn)行求解。函數(shù)filter的語句格式為y=filter(b,a,x)其中，x為輸入的離散序列；y

2、為輸出的離散序列；y的長度與x的長度一樣；b與a分別為差分方程右端與左端的系數(shù)向量?！緦嵗?-1】 已知某LTI系統(tǒng)的差分方程為試用MATLAB命令繪出當(dāng)激勵信號為時，該系統(tǒng)的零狀態(tài)響應(yīng)。解：MATLAB源程序為>>a=3 -4 2;>>b=1 2;>>n=0:30;>>x=(1/2).n;>>y=filter(b,a,x);>>stem(n,y,'fill'),grid on>>xlabel('n'),title('系統(tǒng)響應(yīng)y(n)')程序運行結(jié)果如圖1-1所示

3、。圖1-1 實例1-1系統(tǒng)的零狀態(tài)響應(yīng)1.2.2 離散時間系統(tǒng)的單位取樣響應(yīng)系統(tǒng)的單位取樣響應(yīng)定義為系統(tǒng)在激勵下系統(tǒng)的零狀態(tài)響應(yīng)，用表示。MATLAB求解單位取樣響應(yīng)可利用函數(shù)filter，并將激勵設(shè)為前面所定義的impDT函數(shù)。例如，求解實例1-1中系統(tǒng)的單位取樣響應(yīng)時，MATLAB源程序為>>a=3 -4 2;>>b=1 2;>>n=0:30;>>x=impDT(n);>>h=filter(b,a,x);>>stem(n,h,'fill'),grid on>>xlabel('n

4、9;),title('系統(tǒng)單位取樣響應(yīng)h(n)')程序運行結(jié)果如圖1-2所示。圖1-2 實例1-1的系統(tǒng)單位取樣響應(yīng)MATLAB另一種求單位取樣響應(yīng)的方法是利用控制系統(tǒng)工具箱提供的函數(shù)impz來實現(xiàn)。impz函數(shù)的常用語句格式為impz(b,a,N)其中，參數(shù)N通常為正整數(shù)，代表計算單位取樣響應(yīng)的樣值個數(shù)?！緦嵗?-2】 已知某LTI系統(tǒng)的差分方程為利用MATLAB的impz函數(shù)繪出該系統(tǒng)的單位取樣響應(yīng)。解：MATLAB源程序為>>a=3 -4 2;>>b=1 2;>>n=0:30;>>impz(b,a,30),grid on&g

5、t;>title('系統(tǒng)單位取樣響應(yīng)h(n)')程序運行結(jié)果如圖1-3所示，比較圖1-2和圖1-3，不難發(fā)現(xiàn)結(jié)果相同。圖1-3 系統(tǒng)單位取樣響應(yīng)1.2.3 離散時間信號的卷積和運算由于系統(tǒng)的零狀態(tài)響應(yīng)是激勵與系統(tǒng)的單位取樣響應(yīng)的卷積，因此卷積運算在離散時間信號處理領(lǐng)域被廣泛應(yīng)用。離散時間信號的卷積定義為 （1-2）可見，離散時間信號的卷積運算是求和運算，因而常稱為“卷積和”。MATLAB求離散時間信號卷積和的命令為conv，其語句格式為y=conv(x,h)其中，x與h表示離散時間信號值的向量；y為卷積結(jié)果。用MATLAB進(jìn)行卷積和運算時，無法實現(xiàn)無限的累加，只能計算時限

6、信號的卷積。例如，利用MALAB的conv命令求兩個長為4的矩形序列的卷積和，即，其結(jié)果應(yīng)是長為7（4+4-1=7）的三角序列。用向量1 1 1 1表示矩形序列，MATLAB源程序為>>x1=1 1 1 1;>>x2=1 1 1 1;>>g=conv(x1,x2)g=1 2 3 4 3 2 1如果要繪出圖形來，則利用stem命令，即>>n=1:7;>>stem(n,g,'fill'),grid on,xlabel('n')程序運行結(jié)果如圖1-4所示。圖1-4 卷積結(jié)果圖對于給定函數(shù)的卷積和，我們應(yīng)計算卷

7、積結(jié)果的起始點及其長度。兩個時限序列的卷積和長度一般等于兩個序列長度的和減1。【實例1-3】 已知某系統(tǒng)的單位取樣響應(yīng)為，試用MATLAB求當(dāng)激勵信號為時，系統(tǒng)的零狀態(tài)響應(yīng)。解：MATLAB中可通過卷積求解零狀態(tài)響應(yīng)，即。由題意可知，描述向量的長度至少為8，描述向量的長度至少為4，因此為了圖形完整美觀，我們將向量和向量加上一些附加的零值。MATLAB源程序為>>nx=-1:5; %x(n)向量顯示范圍(添加了附加的零值)>>nh=-2:10; %h(n)向量顯示范圍(添加了附加的零值)>>x=uDT(nx)-uDT(nx-4);>>h=0.8.n

8、h.*(uDT(nh)-uDT(nh-8);>>y=conv(x,h);>>ny1=nx(1)+nh(1); %卷積結(jié)果起始點>>%卷積結(jié)果長度為兩序列長度之和減1,即0到(length(nx)+length(nh)-2)>>%因此卷積結(jié)果的時間范圍是將上述長度加上起始點的偏移值>>ny=ny1+(0:(length(nx)+length(nh)-2);>>subplot(311)>>stem(nx,x,'fill'),grid on>>xlabel('n'),tit

9、le('x(n)')>>axis(-4 16 0 3)>>subplot(312)>>stem(nh,h','fill'),grid on>>xlabel('n'),title('h(n)')>>axis(-4 16 0 3)>>subplot(313)>>stem(ny,y,'fill'),grid on>>xlabel('n'),title('y(n)=x(n)*h(n)')圖

10、1-5 利用卷積和法求解系統(tǒng)的零狀態(tài)響應(yīng)>>axis(-4 16 0 3)程序運行結(jié)果如圖1-5所示。1.3 編程練習(xí)1. 試用MATLAB命令求解以下離散時間系統(tǒng)的單位取樣響應(yīng)，并判斷系統(tǒng)的穩(wěn)定性。（1）（2）2. 已知某系統(tǒng)的單位取樣響應(yīng)為，試用MATLAB求當(dāng)激勵信號為時，系統(tǒng)的零狀態(tài)響應(yīng)。附：1. 單位取樣序列單位取樣序列，也稱為單位沖激序列，定義為 （12-1）要注意，單位沖激序列不是單位沖激函數(shù)的簡單離散抽樣，它在n=0處是取確定的值1。在MATLAB中，沖激序列可以通過編寫以下的impDT.m文件來實現(xiàn)，即function y=impDT(n)y=(n=0); %當(dāng)參

11、數(shù)為0時沖激為1,否則為0調(diào)用該函數(shù)時n必須為整數(shù)或整數(shù)向量?！緦嵗?-1】 利用MATLAB的impDT函數(shù)繪出單位沖激序列的波形圖。解：MATLAB源程序為>>n=-3:3;>>x=impDT(n);>>stem(n,x,'fill'),xlabel('n'),grid on>>title('單位沖激序列')>>axis(-3 3 -0.1 1.1)圖2-1 單位沖激序列程序運行結(jié)果如圖12-1所示。2. 單位階躍序列單位階躍序列定義為 （12-2）在MATLAB中，沖激序列可以通過

12、編寫uDT.m文件來實現(xiàn)，即function y=uDT(n)y=n>=0; %當(dāng)參數(shù)為非負(fù)時輸出1調(diào)用該函數(shù)時n也同樣必須為整數(shù)或整數(shù)向量。【實例2-2】 利用MATLAB的uDT函數(shù)繪出單位階躍序列的波形圖。解：MATLAB源程序為>>n=-3:5;>>x=uDT(n);>>stem(n,x,'fill'),xlabel('n'),grid on>>title('單位階躍序列')>>axis(-3 5 -0.1 1.1)圖2-2 單位階躍序列程序運行結(jié)果如圖12-2所示。第2章

13、 z變換及離散時間LTI系統(tǒng)的z域分析2.1 實驗?zāi)康膌 學(xué)會運用MATLAB求離散時間信號的z變換和z反變換；l 學(xué)會運用MATLAB分析離散時間系統(tǒng)的系統(tǒng)函數(shù)的零極點；l 學(xué)會運用MATLAB分析系統(tǒng)函數(shù)的零極點分布與其時域特性的關(guān)系；l 學(xué)會運用MATLAB進(jìn)行離散時間系統(tǒng)的頻率特性分析。2.2 實驗原理及實例分析2.2.1 z正反變換序列的z變換定義為 （2-1）其中，符號表示取z變換，z是復(fù)變量。相應(yīng)地，單邊z變換定義為 （2-2）MATLAB符號數(shù)學(xué)工具箱提供了計算離散時間信號單邊z變換的函數(shù)ztrans和z反變換函數(shù)iztrans，其語句格式分別為Z=ztrans(x)x=izt

14、rans(z)上式中的x和Z分別為時域表達(dá)式和z域表達(dá)式的符號表示，可通過sym函數(shù)來定義?！緦嵗?-1】 試用ztrans函數(shù)求下列函數(shù)的z變換。（1）； （2）。解：（1）z變換MATLAB源程序為>>x=sym('an*cos(pi*n)');>>Z=ztrans(x);>>simplify(Z)ans=z/(z+a)（2）z變換MATLAB源程序為>>x=sym('2(n-1)-(-2)(n-1)');>>Z=ztrans(x);>>simplify(Z)ans=z2/(z-2)/(

15、z+2)【實例2-2】 試用iztrans函數(shù)求下列函數(shù)的z反變換。（1） （2）解：（1）z反變換MATLAB源程序為>>Z=sym('(8*z-19)/(z2-5*z+6)');>>x=iztrans(Z);>>simplify(x)ans=-19/6*charfcn0(n)+5*3(n-1)+3*2(n-1)其中，charfcn0(n)是函數(shù)在MATLAB符號工具箱中的表示，反變換后的函數(shù)形式為。（2）z反變換MATLAB源程序為>>Z=sym('z*(2*z2-11*z+12)/(z-1)/(z-2)3')

16、;>>x=iztrans(Z);>>simplify(x)ans=-3+3*2n-1/4*2n*n-1/4*2n*n2其函數(shù)形式為。如果信號的z域表示式是有理函數(shù)，進(jìn)行z反變換的另一個方法是對進(jìn)行部分分式展開，然后求各簡單分式的z反變換。設(shè)的有理分式表示為 （2-3）MATLAB信號處理工具箱提供了一個對進(jìn)行部分分式展開的函數(shù)residuez，其語句格式為R,P,K=residuez(B,A)其中，B，A分別表示X(z)的分子與分母多項式的系數(shù)向量；R為部分分式的系數(shù)向量；P為極點向量；K為多項式的系數(shù)。若X(z)為有理真分式，則K為零。【實例2-3】 試用MATLAB

17、命令對函數(shù)進(jìn)行部分分式展開，并求出其z反變換。解：MATLAB源程序為>>B=18;>>A=18,3,-4,-1;>>R,P,K=residuez(B,A)R= 0.3600 0.2400 0.4000P= 0.5000 -0.3333 -0.3333K= 從運行結(jié)果可知，表示系統(tǒng)有一個二重極點。所以，X(z)的部分分式展開為因此，其z反變換為2.2.2 系統(tǒng)函數(shù)的零極點分析離散時間系統(tǒng)的系統(tǒng)函數(shù)定義為系統(tǒng)零狀態(tài)響應(yīng)的z變換與激勵的z變換之比，即 （2-4）如果系統(tǒng)函數(shù)的有理函數(shù)表示式為 （2-5）那么，在MATLAB中系統(tǒng)函數(shù)的零極點就可通過函數(shù)roots

18、得到，也可借助函數(shù)tf2zp得到，tf2zp的語句格式為Z,P,K=tf2zp(B,A)其中，B與A分別表示的分子與分母多項式的系數(shù)向量。它的作用是將的有理分式表示式轉(zhuǎn)換為零極點增益形式，即 （2-6）【實例2-4】 已知一離散因果LTI系統(tǒng)的系統(tǒng)函數(shù)為試用MATLAB命令求該系統(tǒng)的零極點。解：用tf2zp函數(shù)求系統(tǒng)的零極點，MATLAB源程序為>>B=1,0.32;>>A=1,1,0.16;>>R,P,K=tf2zp(B,A)R= -0.3200P= -0.8000 -0.2000K= 1因此，零點為，極點為與。若要獲得系統(tǒng)函數(shù)的零極點分布圖，可直接應(yīng)用z

19、plane函數(shù)，其語句格式為zplane(B,A)其中，B與A分別表示的分子和分母多項式的系數(shù)向量。它的作用是在Z平面上畫出單位圓、零點與極點?！緦嵗?-5】 已知一離散因果LTI系統(tǒng)的系統(tǒng)函數(shù)為試用MATLAB命令繪出該系統(tǒng)的零極點分布圖。解：用zplane函數(shù)求系統(tǒng)的零極點，MATLAB源程序為>>B=1,0,-0.36;>>A=1,-1.52,0.68;>>zplane(B,A),grid on>>legend('零點','極點')>>title('零極點分布圖')程序運行結(jié)果如圖

20、14-1所示?？梢姡撘蚬到y(tǒng)的極點全部在單位圓內(nèi)，故系統(tǒng)是穩(wěn)定的。圖2-1 零極點分布圖2.2.3 系統(tǒng)函數(shù)的零極點分布與其時域特性的關(guān)系與拉氏變換在連續(xù)系統(tǒng)中的作用類似，在離散系統(tǒng)中，z變換建立了時域函數(shù)與z域函數(shù)之間的對應(yīng)關(guān)系。因此，z變換的函數(shù)從形式可以反映的部分內(nèi)在性質(zhì)。我們?nèi)耘f通過討論的一階極點情況，來說明系統(tǒng)函數(shù)的零極點分布與系統(tǒng)時域特性的關(guān)系?！緦嵗?-6】 試用MATLAB命令畫出現(xiàn)下列系統(tǒng)函數(shù)的零極點分布圖、以及對應(yīng)的時域單位取樣響應(yīng)的波形，并分析系統(tǒng)函數(shù)的極點對時域波形的影響。（1） （2） （3）（4） （5） （6）（7）解：MATLAB源程序為>>b1=

21、1,0;>>a1=1,-0.8;>>subplot(121)>>zplane(b1,a1)>>title('極點在單位圓內(nèi)的正實數(shù)')>>subplot(122)>>impz(b1,a1,30);grid on;>>figure>>b2=1,0;>>a2=1,0.8;>>subplot(121)>>zplane(b2,a2)>>title('極點在單位圓內(nèi)的負(fù)實數(shù)')>>subplot(122)>>

22、;impz(b2,a2,30);grid on;>>figure>>b3=1,0;>>a3=1,-1.2,0.72;>>subplot(121)>>zplane(b3,a3)>>title('極點在單位圓內(nèi)的共軛復(fù)數(shù)')>>subplot(122)>>impz(b3,a3,30);grid on;>>figure>>b4=1,0;>>a4=1,-1;>>subplot(121)>>zplane(b4,a4)>>t

23、itle('極點在單位圓上為實數(shù)1')>>subplot(122)>>impz(b4,a4);grid on;>>figure>>b5=1,0;>>a5=1,-1.6,1;>>subplot(121)>>zplane(b5,a5)>>title('極點在單位圓上的共軛復(fù)數(shù)')>>subplot(122)>>impz(b5,a5,30);grid on;>>figure>>b6=1,0;>>a6=1,-1.2;

24、>>subplot(121)>>zplane(b6,a6)>>title('極點在單位圓外的正實數(shù)')>>subplot(122)>>impz(b6,a6,30);grid on;>>figure>>b7=1,0;>>a7=1,-2,1.36;>>subplot(121)>>zplane(b7,a7)>>title('極點在單位圓外的共軛復(fù)數(shù)')>>subplot(122)>>impz(b7,a7,30);g

25、rid on;程序運行結(jié)果分別如圖2-2的（a）、（b）、（c）、（d）、（e）、（f）、（g）所示。(a)(b)(c)(d)(e)(f)圖2-2 系統(tǒng)函數(shù)的零極點分布與其時域特性的關(guān)系(g)從圖14-2可知，當(dāng)極點位于單位圓內(nèi)時，為衰減序列；當(dāng)極點位于單位圓上時，為等幅序列；當(dāng)極點位于單位圓外時，為增幅序列。若有一階實數(shù)極點，則為指數(shù)序列；若有一階共軛極點，則為指數(shù)振蕩序列；若的極點位于虛軸左邊，則序列按一正一負(fù)的規(guī)律交替變化。2.2.4 離散時間LTI系統(tǒng)的頻率特性分析對于因果穩(wěn)定的離散時間系統(tǒng)，如果激勵序列為正弦序列，則系統(tǒng)的穩(wěn)態(tài)響應(yīng)為。其中，通常是復(fù)數(shù)。離散時間系統(tǒng)的頻率響應(yīng)定義為 （

26、2-7）其中，稱為離散時間系統(tǒng)的幅頻特性；稱為離散時間系統(tǒng)的相頻特性；是以（，若零，）為周期的周期函數(shù)。因此，只要分析在范圍內(nèi)的情況，便可分析出系統(tǒng)的整個頻率特性。MATLAB提供了求離散時間系統(tǒng)頻響特性的函數(shù)freqz，調(diào)用freqz的格式主要有兩種。一種形式為H,w=freqz(B,A,N)其中，B與A分別表示的分子和分母多項式的系數(shù)向量；N為正整數(shù)，默認(rèn)值為512；返回值w包含范圍內(nèi)的N個頻率等分點；返回值H則是離散時間系統(tǒng)頻率響應(yīng)在范圍內(nèi)N個頻率處的值。另一種形式為H,w=freqz(B,A,N,whole)與第一種方式不同之處在于角頻率的范圍由擴展到?！緦嵗?-6】 用MATLAB命

27、令繪制系統(tǒng)的頻率響應(yīng)曲線。 解：利用函數(shù)freqz計算出，然后利用函數(shù)abs和angle分別求出幅頻特性與相頻特性，最后利用plot命令繪出曲線。MATLAB源程序為>>b=1 -0.96 0.9028;>>a=1 -1.56 0.8109;>>H,w=freqz(b,a,400,'whole');>>Hm=abs(H); >>Hp=angle(H);>>subplot(211)>>plot(w,Hm),grid on>>xlabel('omega(rad/s)'),

28、ylabel('Magnitude')>>title('離散系統(tǒng)幅頻特性曲線')>>subplot(212)>>plot(w,Hp),grid on>>xlabel('omega(rad/s)'),ylabel('Phase')>>title('離散系統(tǒng)相頻特性曲線')程序運行結(jié)果如圖2-3所示。圖4-3 離散系統(tǒng)頻響特性曲線2.3 編程練習(xí)1. 試用MATLAB的residuez函數(shù)，求出的部分分式展開和。結(jié)果：MATLAB源程序為>>B=2

29、 16 44 56 32;>>A=3,3，-15,18，-12;>>R,P,K=residuez(B,A)R = -0.0177 9.4914 -3.0702 + 2.3398i -3.0702 - 2.3398iP = -3.2361 1.2361 0.5000 + 0.8660i 0.5000 - 0.8660iK = -2.66672. 試用MATLAB畫出下列因果系統(tǒng)的系統(tǒng)函數(shù)零極點分布圖，并判斷系統(tǒng)的穩(wěn)定性。（1） 結(jié)果：程序：B=2,-1.6,-0.9;A=1,-2.5,1.96,-0.48;R,P,K=tf2zp(B,A)zplane(B,A) （2）結(jié)果

30、：程序：B=1,-1;A=1,-0.9,-0.65,0.873;R,P,K=tf2zp(B,A)zplane(B,A)3. 試用MATLAB繪制系統(tǒng)的頻率響應(yīng)曲線。結(jié)果：程序： b=1; a=1,-0.75,0.125; H,w=freqz(b,a,400,'whole'); Hm=abs(H); Hp=angle(H); subplot(211) plot(w,Hm),grid on xlabel('omega(rad/s)'),ylabel('Magnitude') title('離散系統(tǒng)幅頻特性曲線') subplot(21

31、2) plot(w,Hp),grid on xlabel('omega(rad/s)'),ylabel('Phase') title('離散系統(tǒng)相頻特性曲線')實驗三 時域采樣與頻域采樣1. 實驗?zāi)康?時域采樣理論與頻域采樣理論是數(shù)字信號處理中的重要理論。要求掌握模擬信號采樣前后頻譜的變化，以及如何選擇采樣頻率才能使采樣后的信號不丟失信息；要求掌握頻率域采樣會引起時域周期化的概念，以及頻率域采樣定理及其對頻域采樣點數(shù)選擇的指導(dǎo)作用。2. 實驗原理與方法 時域采樣定理的要點是：a) 對模擬信號以間隔T進(jìn)行時域等間隔理想采樣，形成的采樣信號的頻譜是原

32、模擬信號頻譜以采樣角頻率（）為周期進(jìn)行周期延拓。公式為： b) 采樣頻率必須大于等于模擬信號最高頻率的兩倍以上，才能使采樣信號的 頻譜不產(chǎn)生頻譜混疊。 利用計算機計算上式并不方便，下面我們導(dǎo)出另外一個公式，以便用計算機上進(jìn)行實驗。 理想采樣信號和模擬信號之間的關(guān)系為： 對上式進(jìn)行傅立葉變換，得到： 在上式的積分號內(nèi)只有當(dāng)時，才有非零值，因此： 上式中，在數(shù)值上，再將代入，得到： 上式的右邊就是序列的傅立葉變換，即 上式說明理想采樣信號的傅立葉變換可用相應(yīng)的采樣序列的傅立葉變換得到，只要將自變量用代替即可。 頻域采樣定理的要點是：a) 對信號x(n)的頻譜函數(shù)X(ej)在0，2上等間隔采樣N點，

33、得到則N點IDFT得到的序列就是原序列x(n)以N為周期進(jìn)行周期延拓后的主值區(qū)序列，公式為： b) 由上式可知，頻域采樣點數(shù)N必須大于等于時域離散信號的長度M(即NM)，才能使時域不產(chǎn)生混疊，則N點IDFT得到的序列就是原序列x(n),即=x(n)。如果N>M，比原序列尾部多N-M個零點；如果N<M，z則=IDFT發(fā)生了時域混疊失真，而且的長度N也比x(n)的長度M短，因此。與x(n)不相同。 在數(shù)字信號處理的應(yīng)用中，只要涉及時域或者頻域采樣，都必須服從這兩個采樣理論的要點。 對比上面敘述的時域采樣原理和頻域采樣原理，得到一個有用的結(jié)論，這兩個采樣理論具有對偶性：“時域采樣頻譜周期

34、延拓，頻域采樣時域信號周期延拓”。因此放在一起進(jìn)行實驗。3. 實驗內(nèi)容及步驟（1）時域采樣理論的驗證。給定模擬信號， 式中A=444.128，=50，=50rad/s，它的幅頻特性曲線如圖10.2.1 圖3.1 的幅頻特性曲線 現(xiàn)用DFT(FFT)求該模擬信號的幅頻特性，以驗證時域采樣理論。 安照的幅頻特性曲線，選取三種采樣頻率，即=1kHz，300Hz，200Hz。觀測時間選。 為使用DFT，首先用下面公式產(chǎn)生時域離散信號，對三種采樣頻率，采樣序列按順序用，表示。 因為采樣頻率不同，得到的，的長度不同， 長度（點數(shù)）用公式計算。選FFT的變換點數(shù)為M=64，序列長度不夠64的尾部加零。X(k

35、)=FFTx(n) ， k=0,1,2,3,-,M-1式中k代表的頻率為 。要求： 編寫實驗程序，計算、和的幅度特性，并繪圖顯示。觀察分析頻譜混疊失真。（2）頻域采樣理論的驗證。給定信號如下： 編寫程序分別對頻譜函數(shù)在區(qū)間上等間隔采樣32和16點，得到： 再分別對進(jìn)行32點和16點IFFT，得到： 分別畫出、的幅度譜，并繪圖顯示x(n)、的波形，進(jìn)行對比和分析，驗證總結(jié)頻域采樣理論。提示：頻域采樣用以下方法容易變程序?qū)崿F(xiàn)。 直接調(diào)用MATLAB函數(shù)fft計算就得到在的32點頻率域采樣 抽取的偶數(shù)點即可得到在的16點頻率域采樣，即。 當(dāng)然也可以按照頻域采樣理論，先將信號x(n)以16為周期進(jìn)行周

36、期延拓，取其主值區(qū)（16點），再對其進(jìn)行16點DFT(FFT),得到的就是在的16點頻率域采樣。 4思考題： 如果序列x(n)的長度為M，希望得到其頻譜在上的N點等間隔采樣，當(dāng)N<M時， 如何用一次最少點數(shù)的DFT得到該頻譜采樣？5. 實驗報告及要求a) 運行程序打印要求顯示的圖形。 b) 分析比較實驗結(jié)果，簡述由實驗得到的主要結(jié)論c) 簡要回答思考題d) 附上程序清單和有關(guān)曲線。4、實驗程序清單1） 時域采樣理論的驗證程序清單% 時域采樣理論驗證程序exp2a.mTp=64/1000;%觀察時間Tp=64微秒%產(chǎn)生M長采樣序列x(n)% Fs=1000;T=1/Fs;Fs=1000;T

37、=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;alph=pi*50*20.5;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M點FFTxnt)subplot(3,2,1);stem(n,xnt,.);grid on;%調(diào)用繪圖函數(shù)stem繪制序列圖title('(a) Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk);title('(a) T*FTxa(nT),Fs=1000Hz'

38、);xlabel('f(Hz)');ylabel('幅度');axis(0,Fs,0,1.2*max(abs(Xk)%=% Fs=300Hz和 Fs=200Hz的程序與上面Fs=1000Hz完全相同。(1) Fs=1000Hz的圖片：如a所示(2) Fs=200Hz的圖片：如b所示(3) Fs=300Hz的圖片：如c所示(4) 綜合圖片如e所示2 頻域采樣理論的驗證程序清單%頻域采樣理論驗證程序exp2b.mM=27;N=32;n=0:M;%產(chǎn)生M長三角波序列x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb;

39、Xk=fft(xn,1024);%1024點FFTx(n), 用于近似序列x(n)的TFX32k=fft(xn,32);%32點FFTx(n)x32n=ifft(X32k);%32點IFFTX32(k)得到x32(n)X16k=X32k(1:2:N);%隔點抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16點IFFTX16(k)得到x16(n)subplot(3,2,2);stem(n,xn,'.');box ontitle('(b) 三角波序列x(n)');xlabel('n');ylabel('x(n)'

40、;);axis(0,32,0,20)k=0:1023;wk=k/1024;%subplot(3,2,1);plot(wk,abs(Xk);title('(a)FTx(n)');xlabel('omega/pi');ylabel('|X(ejomega)|');axis(0,1,0,200)k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box ontitle('(c) 16點頻域采樣');xlabel('k');ylabel('|X_1_6(k)

41、|');axis(0,8,0,200)n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box ontitle('(d) 16點IDFTX_1_6(k)');xlabel('n');ylabel('x_1_6(n)');axis(0,32,0,20)k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box ontitle('(e) 32點頻域采樣');xlabel('k');ylabel('

42、|X_3_2(k)|');axis(0,16,0,200)n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box ontitle('(f) 32點IDFTX_3_2(k)');xlabel('n');ylabel('x_3_2(n)');axis(0,32,0,20)結(jié)果如圖d所示實驗四 用FFT對信號作頻譜分析1實驗?zāi)康?學(xué)習(xí)用FFT對連續(xù)信號和時域離散信號進(jìn)行譜分析的方法，了解可能出現(xiàn)的分析 誤差及其原因，以便正確應(yīng)用FFT。2. 實驗原理用FFT對信號作頻譜分析是學(xué)習(xí)數(shù)字信號處理的

43、重要內(nèi)容。經(jīng)常需要進(jìn)行譜分析的信號是模擬信號和時域離散信號。對信號進(jìn)行譜分析的重要問題是頻譜分辨率D和分析誤差。頻譜分辨率直接和FFT的變換區(qū)間N有關(guān)，因為FFT能夠?qū)崿F(xiàn)的頻率分辨率是，因此要求。可以根據(jù)此式選擇FFT的變換區(qū)間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續(xù)譜，只有當(dāng)N較大時離散譜的包絡(luò)才能逼近于連續(xù)譜，因此N要適當(dāng)選擇大一些。周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。如果不知道信號周期，可以盡量選擇信號的觀察時間長一些。對模擬信號進(jìn)行譜分析時，首先要按照采樣定理將其變成時域離散信號。如果是

44、模擬周期信號，也應(yīng)該選取整數(shù)倍周期的長度，經(jīng)過采樣后形成周期序列，按照周期序列的譜分析進(jìn)行。3實驗步驟及內(nèi)容（1）對以下序列進(jìn)行譜分析。 選擇FFT的變換區(qū)間N為8和16 兩種情況進(jìn)行頻譜分析。分別打印其幅頻特性曲線。 并進(jìn)行對比、分析和討論。（2）對以下周期序列進(jìn)行譜分析。 選擇FFT的變換區(qū)間N為8和16 兩種情況分別對以上序列進(jìn)行頻譜分析。分別打印其幅頻特性曲線。并進(jìn)行對比、分析和討論。（3）對模擬周期信號進(jìn)行譜分析 選擇 采樣頻率，變換區(qū)間N=16,32,64 三種情況進(jìn)行譜分析。分別打印其幅頻特性，并進(jìn)行分析和討論。 4思考題（1）對于周期序列，如果周期不知道，如何用FFT進(jìn)行譜分析

45、？（2）如何選擇FFT的變換區(qū)間？（包括非周期信號和周期信號）（3）當(dāng)N=8時，和的幅頻特性會相同嗎？為什么？N=16 呢？5實驗報告要求（1）完成各個實驗任務(wù)和要求。附上程序清單和有關(guān)曲線。（2）簡要回答思考題。附錄： 實驗程序清單function y=mstem(x,N)k=0:N-1;wk=2*k/N;stem(wk,abs(x),'.');% 用FFT對信號作頻譜分析clear all;close all%實驗內(nèi)容(1)=x1n=ones(1,4); %產(chǎn)生序列向量x1(n)=R4(n)M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb;

46、%產(chǎn)生長度為8的三角波序列x2(n)x3n=xb,xa;X1k8=fft(x1n,8); %計算x1n的8點DFTX1k16=fft(x1n,16); %計算x1n的16點DFTX2k8=fft(x2n,8); %計算x1n的8點DFTX2k16=fft(x2n,16); %計算x1n的16點DFTX3k8=fft(x3n,8); %計算x1n的8點DFTX3k16=fft(x3n,16); %計算x1n的16點DFT%以下繪制幅頻特性曲線subplot(2,2,1);mstem(X1k8); %繪制8點DFT的幅頻特性圖title('(1a) 8點DFTx_1(n)');xl

47、abel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k8)subplot(2,2,3);mstem(X1k16); %繪制16點DFT的幅頻特性圖title('(1b)16點DFTx_1(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k16)figure(2)subplot(2,2,1);mstem(X2k8); %繪制8點DFT的幅頻特性圖title('(2a) 8點DFTx_2(n)');

48、xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k8)subplot(2,2,2);mstem(X2k16); %繪制16點DFT的幅頻特性圖title('(2b)16點DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k16)subplot(2,2,3);mstem(X3k8); %繪制8點DFT的幅頻特性圖title('(3a) 8點DFTx_3(n)');xlabel(

49、'/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k8)subplot(2,2,4);mstem(X3k16); %繪制16點DFT的幅頻特性圖title('(3b)16點DFTx_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k16)%實驗內(nèi)容(2) 周期序列譜分析=N=8;n=0:N-1; %FFT的變換區(qū)間N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8

50、=fft(x4n); %計算x4n的8點DFTX5k8=fft(x5n); %計算x5n的8點DFTN=16;n=0:N-1; %FFT的變換區(qū)間N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n); %計算x4n的16點DFTX5k16=fft(x5n); %計算x5n的16點DFTfigure(3)subplot(2,2,1);mstem(X4k8); %繪制8點DFT的幅頻特性圖title('(4a) 8點DFTx_4(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); %繪制16點DFT的幅頻特性圖title('(4b)16點DFTx_4(n)');xlabel('/