1、實驗報告2014-2015學(xué)年第二學(xué)期開 課 單 位 適用年級、專業(yè) 13級所有專業(yè) 課 程 名 稱 數(shù)字信號處理實驗 主 講 教 師 課 程 序 號 課 程 代 碼 實 驗 名 稱 實 驗 學(xué) 時 學(xué) 號 姓 名 姓名： 學(xué)實驗1 用MATLAB產(chǎn)生時域離散信號一、.實驗?zāi)康模?、了解常用時域離散信號及其特點。2、掌握用MATLAB產(chǎn)生時域離散信號的方法。二、.實驗原理：1、時域離散信號的概念在時間軸的離散點上取值的信號，稱為離散時間信號。通常，離散時間信號用x(n)表示，其幅度可以在某一范圍內(nèi)連續(xù)取值。由于信號處理設(shè)備或裝置（如計算機、專用的信號處理芯片等）均以有限位的二進制數(shù)來表示信號的

2、幅度，因此，信號的幅度也必須離散化。我們把時間和幅度均取離散值的信號稱為時域離散信號或數(shù)字信號。在MATLAB語言中，時域離散信號可以通過編寫程序直接產(chǎn)生。三、.實驗內(nèi)容：1、閱讀并上機驗證實驗原理部分的例題程序，理解每一條語句的含義。改變例題中的有關(guān)參數(shù)（如信號的頻率、周期、幅度、顯示時間的取值范圍、采樣點數(shù)等），觀察對信號波形的影響。2、編寫程序，產(chǎn)生以下離散序列：（1）f(n)=(n) (-3n4)編寫程序：n1=-3;n2=4;n0=0;n=n1:n2;x=n=n0; stem(n,x, filled);axis(n1,n2,0,1.1*max(x);xlabel(時間(n);ylab

3、el(幅度x(n);title(單位脈沖序列);程序運行結(jié)果如圖：（2）f(n)=u(n) (-5n5)編寫程序：n1=-5;n2=5;n0=0;n=n1:n2;x=n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(時間(n);ylabel(幅度 (n);title(單位階躍序列);box程序運行結(jié)果如圖：（3）f(n)= e(0.1+j1.6)n (0n16)編寫程序：n1=16;a=0.1;w=1.6;n=0:n1;x=exp(a+j*w)*n);subplot(2,2,1);plot(n,real(x);title(復(fù)指數(shù)信號的實部

4、);subplot(2,2,3);stem(n,real(x),filled);title(復(fù)指數(shù)序列的實部);subplot(2,2,2);plot(n,imag(x);title(復(fù)指數(shù)信號的虛部);subplot(2,2,4);stem(n,imag(x),filled);title(復(fù)指數(shù)序列的虛部);box程序運行結(jié)果如圖：（4）f(n)=3sin(n/4) (0n20)編寫程序：f=1/8;Um=3;nt=2.5;T=1/f;tn=n*dt;x=Um*sin(2*f*pi*tn);subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.

5、1*max(x);ylabel(x(t);ylabel(x(n);box 程序運行結(jié)果如圖：3、一個連續(xù)的周期性方波信號頻率為200Hz，信號幅度在-1+1V之間，要求在圖形窗口上顯示其兩個周期的波形。以4kHz的頻率對連續(xù)信號進行采樣，編寫程序生成連續(xù)信號和其采樣獲得的離散信號波形。編寫程序：f=200;nt=2;N=20;T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=square(2*f*pi*tn);subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel(x(t);subplot(2,1,

6、2);stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel(x(n);box程序運行結(jié)果如圖：四、實驗總結(jié)1、通過本次實驗，我學(xué)會了MATLAB產(chǎn)生時域離散信號波形。使用MATLAB這個軟件和編寫程序的一些基本知識，同時還加深了我對各種時域離散信號公式與波形的記憶和理解。思考題：1、通過例題程序，你發(fā)現(xiàn)采樣頻率Fs、采樣點數(shù)N、采樣時間間隔dt在程序編寫中有怎樣的聯(lián)系？使用時需注意什么問題？答： dt=T/N N=Fs/f T=1/f 實驗2 離散LSI系統(tǒng)的時域分析.實驗?zāi)康模?、加深對離散系統(tǒng)的差分方程、單位脈沖響應(yīng)、單位階躍響應(yīng)和卷積分析

7、方法的理解。2、初步了解用MATLAB語言進行離散時間系統(tǒng)時域分析的基本方法。3、掌握求解離散時間系統(tǒng)的單位脈沖響應(yīng)、單位階躍響應(yīng)、線性卷積以及差分方程的程序的編寫方法，了解常用子函數(shù)的調(diào)用格式。二、實驗原理：1、離散LSI系統(tǒng)的響應(yīng)與激勵由離散時間系統(tǒng)的時域分析方法可知，一個離散LSI系統(tǒng)的響應(yīng)與激勵可以用如下框圖表示：其輸入、輸出關(guān)系可用以下差分方程描述：三、實驗內(nèi)容：1、輸入并運行例題程序，理解每一條語句的含義。2、已知描述某離散LSI系統(tǒng)的差分方程為2y(n)-3y(n-1)+y(n-2)=x(n-1)，分別用impz和dstep函數(shù)、filtic和filter函數(shù)兩種方法求解系統(tǒng)的單

8、位序列響應(yīng)和單位階躍響應(yīng)。用impz和dstep函數(shù)求解程序清單如下：a=1,-3/2,1/2;b=0,1/2,0;N=32;n=0:N-1;hn=impz(b,a,n);gn=dstep(b,a,n);subplot(1,2,1);stem(n,hn,k);title(系統(tǒng)的單位序列響應(yīng));ylabel(h(n);xlabel(n);axis(0,N,1.1*min(hn),1.1*max(hn);subplot(1,2,2);stem(n,gn,k);title(系統(tǒng)的單位階躍響應(yīng));ylabel(g(n);xlabel(n);axis(0,N,1.1*min(gn),1.1*max(gn

9、);filtic和filter函數(shù)求解程序清單如下：a=1,-3/2,1/2;b=0,1/2,0;N=32;n=0:N-1;xi=filtic(b,a,0);x1=n=0;hn=filter(b,a,x1,xi);x2=n=0;gn=filter(b,a,x2,xi);subplot(1,2,1);stem(n,hn,k);title(系統(tǒng)的單位序列響應(yīng));ylabel(h(n);xlabel(n);axis(0,N,1.1*min(hn),1.1*max(hn);subplot(1,2,2);stem(n,gn,k);title(系統(tǒng)的單位階躍響應(yīng));ylabel(g(n);xlabel(n

10、);axis(0,N,1.1*min(gn),1.1*max(gn);3、編寫程序描繪下列序列的卷積波形：（1）f1(n=u(n),f2(n)=u(n-2), (0n10) 程序清單如下：n1=0:10;N1=length(n1);f1=ones(1,N1);subplot(2,2,1);stem(n1,f1,filled);title(f1(n);n2=2:12;N2=length(n2);f2=ones(1,N2);subplot(2,2,2);stem(n2,f2,filled);title(f2(n);y=conv(f1,f2);subplot(2,1,2);stem(y,filled

11、);x(n)=sin(n/2),h(n)=(0.5)n (-3n4) 程序清單如下：n1=-3:4*pi;f1=sin(n1/2);n2=-3:4*pi;f2=0.5.n2;nys=n1(1)+n2(1);nyf=n1(end)+n2(end);y=conv(f1,f2);ny=nys:nyf;subplot(2,2,1);stem(n1,f1);subplot(2,2,2);stem(n2,f2);subplot(2,1,2);stem(ny,y);4、已知某離散LSI系統(tǒng)的單位序列響應(yīng)為h(n)=3(n-3)+0.5(n-4)+0.2(n-5)+0.7(n-6)-0.8(n-7)求輸入為x

12、(n)=e-0.5nu(n)時的系統(tǒng)響應(yīng)。程序清單如下:N=16;n=0:N-1;x=exp(-0.5*n);subplot(2,2,1);stem(n,x);a=1;b=0,0,0,3,0.5,0.2,0.7,-0.8;hn=impz(b,a,n);subplot(2,2,2);stem(n,hn);y=conv(x,hn);subplot(2,1,2);stem(y);程序運行結(jié)果如下圖：已知描述某離散LSI系統(tǒng)的差分方程為y(n)=0.7y(n-1)+2x(n)-x(n-2)，求輸入為x(n)=u(n-3)時的系統(tǒng)響應(yīng)。 程序清單如下:N=16;n1=3:N+2;f1=zeros(1,3

13、),ones(1,(N-3);subplot(2,2,1);stem(n1,f1);a=1,-0.7,0;b=2,0,-1;f2=impz(b,a,n1);subplot(2,2,2);stem(n1,f2);nys=2*n1(1);nyf=2*n1(end);y=conv(f1,f2);ny=nys:nyf;subplot(2,1,2);stem(ny,y);程序運行結(jié)果如下圖： 實驗3 離散LSI系統(tǒng)的頻域分析一、實驗?zāi)康模?、加深對離散系統(tǒng)變換域分析z變換的理解，掌握使用MATLAB進行z變換和逆z變換的常用函數(shù)的用法。2、了解離散系統(tǒng)的零極點與系統(tǒng)因果性和穩(wěn)定性的關(guān)系，熟悉使用MATL

14、AB進行離散系統(tǒng)的零極點分析的常用函數(shù)的用法。3、加深對離散系統(tǒng)的頻率響應(yīng)特性基本概念的理解，掌握使用MATLAB進行離散系統(tǒng)幅頻響應(yīng)和相頻響應(yīng)特性分析的常用方法。二、實驗原理：1、z變換和逆z變換（1）用ztrans函數(shù)求無限長序列的z變換。該函數(shù)只給出z變換的表達式，而沒有給出收斂域。另外，由于這一函數(shù)還不盡完善，有的序列的z變換還不能求出，逆z變換也存在同樣的問題。例3-1 求以下各序列的z變換 x1(n)=an x2(n)=n x3(n)=n(n-1)/2 x4(n)=ejon x5(n)=1/n(n-1)程序清單如下：syms w0 n z a;x1=an;X1=ztrans(x1)

15、x2=n;X2=ztrans(x2)x3=(n*(n-1)/2;X3=ztrans(x3)x4=exp(j*w0*n);X4=ztrans(x4)x5=1/n*(n-1);X5=ztrans(x5)程序運行結(jié)果如下：X1 =z/a/(z/a-1) X2 =z/(z-1)2X3 =1/2*z*(z+1)/(z-1)3-1/2*z/(z-1)2X4 =z/exp(i*w0)/(z/exp(i*w0)-1)X5 =z/(z-1)-ztrans(1/n,n,z)（2）用iztrans函數(shù)求無限長序列的逆z變換。例3-2 求下列函數(shù)的逆z變換。 程序清單如下：syms n z a;X1=z/(z-1);

16、x1=iztrans(X1)X2=a*z/(a-z)2;x2=iztrans(X2)X3=z/(z-1)3;x3=iztrans(X3)X4=(1-z-n)/(1-z-1);x4=iztrans(X4)程序運行結(jié)果如下：x1 =1x2 =an*nx3 =1/2*n2-1/2*nx4 =iztrans(1-z(-n)/(1-1/z),z,n)2、離散系統(tǒng)的零極點分析（系統(tǒng)極點位置對系統(tǒng)響應(yīng)的影響）例3-3 研究z右半平面的實數(shù)極點對系統(tǒng)的影響。已知系統(tǒng)的零極點增益模型分別為：求這些系統(tǒng)的零極點分布圖以及系統(tǒng)的單位序列響應(yīng)，判斷系統(tǒng)的穩(wěn)定性。程序清單如下：z1=0;p1=0.85;k=1;b1,a

17、1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title(極點在單位圓內(nèi));subplot(3,2,2);impz(b1,a1,20);z2=0;p2=1;b2,a2=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title(極點在單位圓上);subplot(3,2,4);impz(b2,a2,20);z3=0;p3=1.5;b3,a3=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title(極點在單位圓外);subplot(3,2,6);impz(b3,a3,2

18、0);程序運行結(jié)果如圖3-1所示。由圖可見，這三個系統(tǒng)的極點均為實數(shù)且處于z平面的右半平面。由圖可知，當(dāng)極點位于單位圓內(nèi)，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而收斂；當(dāng)極點位于單位圓上，系統(tǒng)的單位序列響應(yīng)為等幅振蕩；當(dāng)極點位于單位圓外，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而發(fā)散。由此可知系統(tǒng)1、2為穩(wěn)定系統(tǒng)。圖3-1例3-4 研究z左半平面的實數(shù)極點對系統(tǒng)的影響。已知系統(tǒng)的零極點增益模型分別為：求這些系統(tǒng)的零極點分布圖以及系統(tǒng)的單位序列響應(yīng)，判斷系統(tǒng)的穩(wěn)定性。程序清單如下：z1=0;p1=-0.85;k=1;b1,a1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p

19、1);title(極點在單位圓內(nèi));subplot(3,2,2);impz(b1,a1,20);z2=0;p2=-1;b2,a2=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title(極點在單位圓上);subplot(3,2,4);impz(b2,a2,20);z3=0;p3=-1.5;b3,a3=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title(極點在單位圓外);subplot(3,2,6);impz(b3,a3,20);程序運行結(jié)果如圖3-2所示。由圖可見，這三個系統(tǒng)的極點均為實數(shù)且處于z平面的

20、左半平面。由圖可知，當(dāng)極點位于單位圓內(nèi)，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而收斂；當(dāng)極點位于單位圓上，系統(tǒng)的單位序列響應(yīng)為等幅振蕩；當(dāng)極點位于單位圓外，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而發(fā)散。由此可知系統(tǒng)1、2為穩(wěn)定系統(tǒng)。圖3-2例3-5 研究z右半平面的復(fù)數(shù)極點對系統(tǒng)響應(yīng)的影響已知系統(tǒng)的零極點增益模型分別為：求這些系統(tǒng)的零極點分布圖以及系統(tǒng)的單位序列響應(yīng)，判斷系統(tǒng)的穩(wěn)定性。程序清單如下：z1=0.3,0;p1=0.5+0.7j,0.5-0.7j;k=1;b1,a1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title(極點在單位圓內(nèi));subplo

21、t(3,2,2);impz(b1,a1,20);z2=0.3,0;p2=0.6+0.8j,0.6-0.8j;b2,a2=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title(極點在單位圓上);subplot(3,2,4);impz(b2,a2,20);z3=0.3,0;p3=1+j,1-j;b3,a3=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title(極點在單位圓外);subplot(3,2,6);impz(b3,a3,20);程序運行結(jié)果如圖3-3所示。由圖可見，這三個系統(tǒng)的極點均為復(fù)數(shù)且處于z平

22、面的右半平面。由圖可知，當(dāng)極點位于單位圓內(nèi)，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而收斂；當(dāng)極點位于單位圓上，系統(tǒng)的單位序列響應(yīng)為等幅振蕩；當(dāng)極點位于單位圓外，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而發(fā)散。由此可知系統(tǒng)1、2為穩(wěn)定系統(tǒng)。圖3-3由以上三例可得結(jié)論：系統(tǒng)只有在極點處于單位圓內(nèi)才是穩(wěn)定的。例3-6 已知某離散時間系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)的零極點及零極點分布圖，并判斷系統(tǒng)的因果穩(wěn)定性。程序清單如下：b=0.2,0.1,0.3,0.1,0.2;a=1,-1.1,1.5,-0.7,0.3;rz=roots(b)rp=roots(a)subplot(2,1,1);zplane(b,a);title(系

23、統(tǒng)的零極點分布圖);subplot(2,1,2);impz(b,a,20);title(系統(tǒng)的單位序列響應(yīng));xlabel(n);ylabel(h(n);程序運行結(jié)果如下：rz = -0.5000 + 0.8660i -0.5000 - 0.8660i 0.2500 + 0.9682i 0.2500 - 0.9682irp = 0.2367 + 0.8915i 0.2367 - 0.8915i 0.3133 + 0.5045i 0.3133 - 0.5045i圖3-4由零極點分布圖可見，該系統(tǒng)的所有極點均在單位圓內(nèi)，因此該系統(tǒng)是一個因果穩(wěn)定系統(tǒng)。3、離散系統(tǒng)的頻率響應(yīng)（1）離散系統(tǒng)的頻率響應(yīng)的

24、基本概念已知穩(wěn)定系統(tǒng)傳遞函數(shù)的零極點增益模型為則系統(tǒng)的頻響函數(shù)為其中，系統(tǒng)的幅頻特性為系統(tǒng)的相頻特性為由以上各式可見，系統(tǒng)函數(shù)與頻率響應(yīng)有著密切的聯(lián)系。適當(dāng)?shù)乜刂葡到y(tǒng)函數(shù)的零極點分布，可以改變離散系統(tǒng)的頻響特性： 在原點（z=0）處的零點或極點至單位圓的距離始終保持不變，其值|ej|=1，所以，對幅度響應(yīng)不起作用； 單位圓附近的零點對系統(tǒng)幅度響應(yīng)的谷值位置及深度有明顯影響； 單位圓內(nèi)且靠近單位圓附近的極點對系統(tǒng)幅度的峰值位置及大小有明顯的影響。（2）系統(tǒng)的頻響特性分析例3-7 已知某離散時間系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)在0頻率范圍內(nèi)的相對幅頻響應(yīng)與相頻響應(yīng)。程序清單如下：b=0.1321,0,-0

25、.3963,0,0.3963,0,-0.1321;a=1,0,0.34319,0,0.60439,0,0.20407;freqz(b,a);程序運行結(jié)果如圖3-5所示。該系統(tǒng)是一個IIR數(shù)字帶通濾波器。其中幅頻特性采用歸一化的相對幅度值，以分貝（dB）為單位。圖3-5例3-8 已知某離散時間系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)在0頻率范圍內(nèi)的絕對幅頻響應(yīng)與相頻響應(yīng)。程序清單如下：b=0.2,0.1,0.3,0.1,0.2;a=1,-1.1,1.5,-0.7,0.3;n=(0:500)*pi/500;h,w=freqz(b,a,n);subplot(2,1,1);plot(n/pi,abs(h);grid;

26、axis(0,1,1.1*min(abs(h),1.1*max(abs(h);xlabel(omega/pi);ylabel(幅度);subplot(2,1,2);plot(n/pi,angle(h);grid;axis(0,1,1.1*min(angle(h),1.1*max(angle(h);xlabel(omega/pi);ylabel(相位);程序運行結(jié)果如圖3-6所示。該系統(tǒng)為一低通濾波器。圖3-6例3-9 已知某離散時間系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)在0頻率范圍內(nèi)的絕對幅頻響應(yīng)與相頻響應(yīng)、相對幅頻響應(yīng)與相頻響應(yīng)及零極點分布圖。程序清單如下：b=0.1,-0.4,0.4,-0.1;a=1,

27、0.3,0.55,0.2;n=(0:500)*pi/500;h,w=freqz(b,a,n);db=20*log10(abs(h);subplot(2,2,1);plot(w/pi,abs(h);grid;axis(0,1,1.1*min(abs(h),1.1*max(abs(h);title(幅頻特性（V）);xlabel(omega/pi);ylabel(幅度(V);subplot(2,2,2);plot(w/pi,angle(h);grid;axis(0,1,1.1*min(angle(h),1.1*max(angle(h);xlabel(omega/pi);ylabel(相位);tit

28、le(相頻特性);subplot(2,2,3);plot(w/pi,db);gridaxis(0,1,-100,5);title(幅頻特性（dB）);subplot(2,2,4);zplane(b,a);title(零極點分布);程序運行結(jié)果如圖3-7所示：圖3-7（3）一個求解頻率響應(yīng)的實用函數(shù)。在實際使用freqz進行離散系統(tǒng)頻響特性分析時。通常需要求解幅頻響應(yīng)、相頻響應(yīng)、群時延，幅頻響應(yīng)又分為絕對幅頻和相對幅頻兩種表示方法。下面定義函數(shù)freqz_m，利用該函數(shù)，可方便求出上述各項。freqz_m函數(shù)定義如下：functiondb,mag,pha,grd,w=freqz_m(b,a);H

29、,w=freqz(b,a,1000,whole);H=(H(1:501);w=(w(1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);例3-10 已知某離散時間系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)在0頻率范圍內(nèi)的絕對幅頻響應(yīng)與相頻響應(yīng)、相對幅頻響應(yīng)與相頻響應(yīng)及群時延。程序清單如下：b=0.1321,0,0.3963,0,0.3963,0,0.1321;a=1,0,-0.34319,0,0.60439,0,-0.20407;db,mag,pha,grd,w=freqz_m(b,a);subplot

30、(2,2,1);plot(w/pi,mag);gridaxis(0,1,1.1*min(mag),1.1*max(mag);title(幅頻特性（V）);xlabel(omega/pi);ylabel(幅度(V);subplot(2,2,2);plot(w/pi,pha);grid;axis(0,1,1.1*min(pha),1.1*max(pha);xlabel(omega/pi);ylabel(相位);title(相頻特性);subplot(2,2,3);plot(w/pi,db);gridaxis(0,1,-100,5);title(幅頻特性（dB）);subplot(2,2,4);pl

31、ot(w/pi,grd);gridaxis(0,1,0,10)title(群時延);程序運行結(jié)果如圖3-8所示：圖3-8三、實驗內(nèi)容：1、輸入并運行例題程序，理解每條語句的含義。2、求下列各序列的z變換：程序清單如下：symsw0na;x1=n*an;X1=ztrans(x1)x2=sin(w0*n);X2=ztrans(x2)x3=exp(-a*n)*sin(w0*n);X3=ztrans(x3)程序運行結(jié)果如下：X1=z*a/(-z+a)2X2=-sin(w0)*z/(-z2+2*z*cos(w0)-1)X3=-sin(w0)*z/exp(-a)/(-z2/exp(-a)2+2*z/exp

32、(-a)*cos(w0)-1)3、求下列函數(shù)的逆z變換程序清單如下：symsw0nza;X1=z/(z-a);x1=iztrans(X1)X2=z/(a-z)2;x2=iztrans(X2)X3=z/z-exp(j*w0);x3=iztrans(X3)X4=(1-z-3)/(1-z-1);x4=iztrans(X4)程序運行結(jié)果如下：x1=anx2=n*an/ax3=charfcn0(n)-iztrans(exp(i*w0),w0,n)x4=charfcn2(n)+charfcn1(n)+charfcn0(n)4、求下列系統(tǒng)函數(shù)所描述的離散系統(tǒng)的零極點分布圖，并判斷系統(tǒng)的穩(wěn)定性（1）z1=-0

33、.3;p1=-1+j,-1-j;k=1;b1,a1=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title(極點在單位圓內(nèi));subplot(3,2,2);impz(b1,a1,20);由圖可見當(dāng)極點位于單位圓內(nèi)，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而收斂；當(dāng)極點位于單位圓上，系統(tǒng)的單位序列響應(yīng)為等幅振蕩；當(dāng)極點位于單位圓外，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而發(fā)散。由此可知系統(tǒng)為穩(wěn)定系統(tǒng)。（2）b=4,-1.6,-1.6,4;a=1,0.4,0.35,-0.4;rz=roots(b)rp=roots(a)subplot(2,1,1);zplane(b,a

34、);title(系統(tǒng)的零極點分布圖);subplot(2,1,2);impz(b,a,20);title(系統(tǒng)的單位序列響應(yīng));xlabel(n);ylabel(h(n);rz=-1.00000.7000+0.7141i0.7000-0.7141irp=-0.4500+0.7730i-0.4500-0.7730i0.5000由零極點分布圖可見，該系統(tǒng)的所有極點均在單位圓內(nèi)，因此該系統(tǒng)是一個因果穩(wěn)定系統(tǒng)。5、已知某離散時間系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)在0頻率范圍內(nèi)的絕對幅頻響應(yīng)與相頻響應(yīng)、相對幅頻響應(yīng)與相頻響應(yīng)及群時延。b=0.187632,0,-0.241242,0,0.241242,0,-187

35、632;a=1,0,0.602012,0,0.495684,0,0.035924;db,mag,pha,grd,w=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag);gridaxis(0,1,1.1*min(mag),1.1*max(mag);title(幅頻特性（V）);xlabel(omega/pi);ylabel(幅度(V);subplot(2,2,2);plot(w/pi,pha);grid;axis(0,1,1.1*min(pha),1.1*max(pha);xlabel(omega/pi);ylabel(相位);title(相頻特性);subpl

36、ot(2,2,3);plot(w/pi,db);gridaxis(0,1,-100,5);title(幅頻特性（dB）);subplot(2,2,4);plot(w/pi,grd);gridaxis(0,1,0,10);title(群時延);四、實驗預(yù)習(xí)：1、認真閱讀實驗原理部分，明確實驗?zāi)康?，?fù)習(xí)有關(guān)離散LSI系統(tǒng)頻率響應(yīng)的理論知識。2、讀懂實驗原理部分的例題程序，熟悉與本實驗有關(guān)的MATLAB函數(shù)。3、根據(jù)實驗內(nèi)容預(yù)先編寫實驗程序，并思考本實驗提出的有關(guān)MATLAB函數(shù)在調(diào)用時應(yīng)注意哪些問題。4、預(yù)習(xí)思考題： 系統(tǒng)函數(shù)零極點的位置與系統(tǒng)單位序列響應(yīng)有何關(guān)系？ 離散系統(tǒng)的零極點對系統(tǒng)幅頻響應(yīng)

37、有何影響？答：系統(tǒng)函數(shù)零極點的位置與系統(tǒng)單位序列響應(yīng)有何關(guān)系？當(dāng)極點位于單位圓內(nèi)，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而收斂；當(dāng)極點位于單位圓上，系統(tǒng)的單位序列響應(yīng)為等幅振蕩；當(dāng)極點位于單位圓外，系統(tǒng)的單位序列響應(yīng)隨著頻率的增大而發(fā)散。離散系統(tǒng)的零極點對系統(tǒng)幅頻響應(yīng)有何影響？在原點（z=0）處的零點或極點至單位圓的距離始終保持不變，其值|ej|=1，所以，對幅度響應(yīng)不起作用；單位圓附近的零點對系統(tǒng)幅度響應(yīng)的谷值位置及深度有明顯影響；單位圓內(nèi)且靠近單位圓附近的極點對系統(tǒng)幅度的峰值位置及大小有明顯的影響。實驗4 DFS、DFT與FFT一、實驗?zāi)康模?、加深對周期序列DFS、有限長序列DFT和FFT的基

38、本概念及理論的理解。2、掌握用MATLAB語言求解DFS、DFT、FFT以及相應(yīng)反變換的方法。3、觀察周期序列的重復(fù)周期數(shù)對序列頻譜特性的影響。 二、實驗原理：1、周期序列的離散傅里葉級數(shù)（DFS）（1）DFS的基本概念：離散時間序列x(n)滿足x(n)= x(n+rN)，稱為離散周期序列，用表示。其中，N為信號的周期，x(n)稱為離散周期序列的主值。周期序列可以用離散傅里葉級數(shù)（DFS）表示：其中，是周期序列DFS第k次諧波分量的系數(shù)，也稱為周期序列的頻譜，可表示為以上兩式也是周期序列的一對傅里葉級數(shù)變換對。令，以上DFS變換對又可以寫成：與連續(xù)周期信號的傅里葉級數(shù)相比，周期序列的離散傅里葉

39、級數(shù)有以下特點： 連續(xù)周期信號的傅里葉級數(shù)由無窮多個與基波頻率成整數(shù)倍的諧波分量疊加而成，而周期為N的周期序列的傅里葉級數(shù)僅有N個獨立的諧波分量。 周期序列的頻譜也是一個以N為周期的周期序列。（2）周期序列的DFS和IDFS例4-1 已知一個周期性矩形脈沖寬度占整個周期的1/4，一個周期的采樣點數(shù)為16點，編程顯示3個周期的序列波形，并： 用傅里葉級數(shù)求信號的幅度譜和相位譜。 求傅里葉級數(shù)逆變換的圖形，并與原序列進行比較。程序清單如下：N=16;xn=ones(1,N/4),zeros(1,3*N/4);xn=xn,xn,xn;n=0:3*N-1;k=0:3*N-1;Xk=xn*exp(-j*

40、2*pi/N).(n*k);x=(Xk*exp(j*2*pi/N).(n*k)/N;subplot(2,2,1);stem(n,xn);title(x(n);axis(-1,3*N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFS|X(k)|);axis(-1,3*N,1.1*min(x),1.1*max(x);subplot(2,2,3),stem(k,abs(Xk);title(|X(k)|);axis(-1,3*N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),s

41、tem(k,angle(Xk);title(arg|X(k)|);axis(-1,3*N,1.1*min(angle(Xk),1.1*max(angle(Xk);程序運行結(jié)果如圖4-1所示：圖4-1由離散傅里葉級數(shù)逆變換圖形可見，與原序列相比，幅度擴大了32倍。這是因為周期序列為原主值序列周期的3倍，做逆變換時未做處理?？蓪⒛孀儞Q程序改為：x=(Xk*exp(j*2*pi/N).(n*k)/3*3*N;由上例可見，周期序列的DFS和IDFS是依據(jù)變換公式編程的，無論信號序列如何變化，求解的公式總是一樣的。因此，可將其編寫成通用子程序： 離散傅里葉級數(shù)變換通用子程序dfs.mfunctionXk

42、=dfs(xn,N)n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;Xk=xn*WN.nk; 離散傅里葉級數(shù)逆變換通用子程序idfs.mfunctionxn=idfs(Xk,N)n=0:N-1;k=0:N-1;WN=exp(j*2*pi/N);nk=n*k;xn=(Xk*WN.nk)/N;例4-2 利用上述兩個子程序，重做例4-1程序清單如下：N=16;xn=ones(1,N/4),zeros(1,3*N/4);n=0:N-1;k=0:N-1;Xk=dfs(xn,N);x=idfs(Xk,N);subplot(2,2,1);stem(n,xn);title(x

43、(n);axis(-1,3*N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFS|X(k)|);axis(-1,3*N,1.1*min(x),1.1*max(x);subplot(2,2,3),stem(k,abs(Xk);title(|X(k)|);axis(-1,3*N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),stem(k,angle(Xk);title(arg|X(k)|);axis(-1,3*N,1.1*min(angle(Xk),1.1*max(angl

44、e(Xk);程序運行結(jié)果如圖4-2所示。由于子程序僅適用于對主值區(qū)間進行變換，周期次數(shù)無法傳遞給子程序，因此程序執(zhí)行結(jié)果僅顯示一個周期的變換情況。圖4-2（2）周期重復(fù)次數(shù)對序列頻譜的影響理論上講，周期序列不滿足絕對可積條件，因此不能用傅里葉級數(shù)來表示。實際處理時可先取K個周期進行處理，然后令K趨于無窮大，分析其極限情況。根據(jù)這一分析思路，可以觀察序列由非周期到周期變化時，頻譜由連續(xù)譜逐漸向離散譜過渡的過程。例4-3 已知一矩形脈沖寬度占整個周期的1/2，一個周期的采樣點數(shù)為10點，用DFS求序列的重復(fù)周期數(shù)分別為1、4、7、10時的幅度譜。程序清單如下：xn=ones(1,5),zeros(

45、1,5);Nx=length(xn);Nw=1000;dw=2*pi/Nw;k=floor(-Nw/2+0.5):(Nw/2+0.5);for r=0:3 K=3*r+1; nx=0:(K*Nx-1); x=xn(mod(nx,Nx)+1); Xk=x*(exp(-j*dw*nx*k)/K;subplot(4,2,2*r+1);stem(nx,x);axis(0,K*Nx-1,0,1.1);ylabel(x(n);subplot(4,2,2*r+2);plot(k*dw,abs(Xk);axis(-4,4,0,1.1*max(abs(Xk);ylabel(X(k);end程序運行結(jié)果如圖4-3

46、所示。由圖可見，序列的重復(fù)周期數(shù)越多，頻譜越是向幾個頻點集中。當(dāng)序列的周期數(shù)趨于無窮大時，頻譜轉(zhuǎn)化為離散譜。圖4-32、離散傅里葉變換（DFT）（1）DFT與IDFT在實際中常常使用有限長序列。如果有限長序列為x(n)，則該序列的離散傅里葉變換對可表示為從離散傅里葉變換定義式可以看出，有限長序列在時域上是離散的，在頻域上也是離散的，式中，即僅在單位圓上N個等間距的點上取值，這為使用計算機進行處理帶來了方便。由有限長序列的傅里葉變換和逆變換定義可知，DFT和DFS的變換公式非常相似，因此，在程序編寫上也基本一致。例4-4 已知x(n)=0,1,2,3,4,5,6,7，求其DFT和IDFT。要求：

47、 畫出序列傅里葉變換對應(yīng)的|X(k)|和argX(k)圖形。 畫出x(n)圖形，并與IDFTX(k)圖形進行比較。程序清單如下：xn=0,1,2,3,4,5,6,7;N=length(xn);n=0:N-1;k=0:N-1;Xk=xn*exp(-j*2*pi/N).(n*k);x=(Xk*exp(j*2*pi/N).(n*k)/N;subplot(2,2,1);stem(n,xn);title(x(n);axis(-1,N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFT|X(k)|);axis(-1,N,1.1*m

48、in(x),1.1*max(x);subplot(2,2,3),stem(k,abs(Xk);title(|X(k)|);axis(-1,N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),stem(k,angle(Xk);title(arg|X(k)|);axis(-1,N,1.1*min(angle(Xk),1.1*max(angle(Xk);程序運行結(jié)果如圖4-4所示。由圖可見，與周期序列不同，有限長序列本身是僅有N點的離散序列，相當(dāng)于周期序列的主值部分。因此，其頻譜也對應(yīng)序列的主值部分，是長度為N的離散序列。圖4-4（2）DFT與DFS的聯(lián)

49、系將周期序列的傅里葉級數(shù)變換對和有限長序列的離散傅里葉變換對進行比較可見，兩者的區(qū)別僅僅是將周期序列換成了有限長序列x(n)，同時，由于式中的周期性，因而有限長序列的離散傅里葉變換實際上隱含著周期性。例4-5 已知周期序列的主值x(n)=0,1,2,3,4,5,6,7，求x(n)的周期重復(fù)次數(shù)為4次時的DFS。要求： 畫出主值序列周期序列的波形。 畫出周期序列傅里葉變換對應(yīng)的和arg的圖形。程序清單如下：xn=0,1,2,3,4,5,6,7;N=length(xn);m=0:N-1;n=0:4*N-1;k=0:4*N-1;xn1=xn(mod(n,N)+1);Xk=xn1*exp(-j*2*p

50、i/N).(n*k);subplot(2,2,1);stem(m,xn);title(x(n);subplot(2,2,2);stem(n,xn1);title(周期序列);subplot(2,2,3),stem(k,abs(Xk);title(|X(k)|);axis(-1,4*N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),stem(k,angle(Xk);title(arg|X(k)|);axis(-1,4*N,1.1*min(angle(Xk),1.1*max(angle(Xk);程序運行結(jié)果如圖4-5所示。與例4-4相比，有限長序列x

51、(n)可以看成是周期序列的一個周期；反之，周期序列可以看成是有限長序列以N為周期的周期延拓。頻域上的情況也一樣。從這個意義上說，周期序列只有有限個序列值有意義。（3）DFT與DTFT的聯(lián)系 若離散時間非周期序列為x(n)，則它的離散傅里葉變換（DTFT）對定義為其中稱為序列的頻譜?？杀硎緸?，稱為序列的幅圖4-5度譜，稱為序列的相位譜。由DTFT的定義可見，序列在時域是離散的、非周期的，在頻域是連續(xù)的、周期的。與有限長序列相比，僅在單位圓上取值，X(k)是在單位圓上N個等間距的點上取值。因此，連續(xù)譜可由離散譜X(k)經(jīng)插值后得到。例4-6 求有限長序列x(n)=0,1,2,3,4,5,6,7的D

52、TFT，將（）區(qū)間分成500份。要求： 畫出序列x(n)的圖形。 畫出由DTFT求出的幅度譜和相位譜的圖形。程序清單如下：xn=0,1,2,3,4,5,6,7;N=length(xn);n=0:N-1;w=linspace(-2*pi,2*pi,500);X=xn*exp(-j*n*w);subplot(3,1,1);stem(n,xn);box on;title(x(n);subplot(3,1,2);plot(w,abs(X);title(|X(ejomega)|);axis(-2*pi,2*pi,1.1*min(abs(X),1.1*max(abs(X);subplot(3,1,3),p

53、lot(w,angle(X);title(ephi(omega);axis(-2*pi,2*pi,1.1*min(angle(X),1.1*max(angle(X);程序運行結(jié)果如圖4-6所示：圖4-6與圖4-4相比，兩者有一定差別。主要原因在于，該例進行DTFT時，在單位圓上取250個點進行分割；而圖4-4進行DFT時，X(k)是在單位圓上以N=8的等間距取值，X(k)的序列長度與相比不夠長。若將x(n)=0,1,2,3,4,5,6,7補零拓展至長度N=100，再求其DFT，如圖4-7所示。與例4-6相比，|X(k)|和argX(k)圖形接近和的圖形。注意圖4-7對應(yīng)0,2區(qū)間。圖4-7程序

54、清單如下：N=100;xn=0,1,2,3,4,5,6,7,zeros(1,N-8);n=0:N-1;k=0:N-1;Xk=xn*exp(-j*2*pi/N).(n*k);x=(Xk*exp(j*2*pi/N).(n*k)/N;subplot(2,2,1);stem(n,xn);title(x(n);axis(-1,N,1.1*min(xn),1.1*max(xn);subplot(2,2,2);stem(n,abs(x);title(IDFT|X(k)|);axis(-1,N,1.1*min(x),1.1*max(x);subplot(2,2,3),stem(k,abs(Xk);title(

55、|X(k)|);axis(-1,N,1.1*min(abs(Xk),1.1*max(abs(Xk);subplot(2,2,4),stem(k,angle(Xk);title(arg|X(k)|);axis(-1,N,1.1*min(angle(Xk),1.1*max(angle(Xk);3、FFT（1）MATLAB提供的FFT函數(shù) 由理論學(xué)習(xí)可知，DFT是唯一在時域和頻域均離散的變換方法，它適用于有限長序列。盡管這種變換方法是可以用于數(shù)值計算的，但如果只是簡單地按照定義進行數(shù)據(jù)處理，當(dāng)序列長度很大時，將占用很大的內(nèi)存空間，且運算時間很長。 快速傅里葉變換是用于DFT運算的高效快速算法的統(tǒng)稱，

56、FFT只是其中的一種。FFT主要有時域抽取算法和頻域抽取算法，基本思想是將一個長度為N的序列分解成多個段序列，如基2算法、基4算法等，大大地縮短了DFT的時間。有關(guān)詳細理論可參考教材。MATLAB提供了進行FFT的函數(shù)fft和ifft分別用于計算DFT和IDFT。例4-7 已知一個長度為8的時域離散信號，n1=0，n2=7，在n0=4前為0，n0以后為1.對其進行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=

57、0:N-1;Xk=fft(xn,N);subplot(2,1,2);stem(k,abs(Xk);title(Xk=DFT(xn);xn1=ifft(Xk,N);subplot(2,2,2);stem(n,xn1);title(x(n)=IDFT(Xk);程序運行結(jié)果如圖4-8所示：圖4-9（2）用FFT進行頻譜分析 對有限長序列進行譜分析一個序號從n1到n2的時域有限長序列x(n)，它的頻譜定義為它的離散傅里葉變換，且在Nyquist頻率范圍內(nèi)有界并連續(xù)。序列的長度為N，則N=n2-n1+1。計算x(n)的DFT得到的是的N個樣本點。其中數(shù)字頻率為式中：為數(shù)字頻率的分辨率；k取對應(yīng)-（N-1

58、）/2到（N-1）/2區(qū)間的整數(shù)。在實際使用中，往往要求計算出信號以模擬頻率為橫坐標(biāo)的頻譜，此時對應(yīng)的模擬頻率為式中：D為模擬頻率的分辨率或頻率間隔；Ts為采樣信號的周期，Ts=1/Fs；定義信號的長度L=NTs。在使用FFT進行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,2。如果要顯示-，范圍的頻譜，則可以使用fftshift(X)進行位移。例4-8 已知有限長序列x(n)=1，2，3，2，1，其采樣頻率Fs=10Hz。請使用FFT計算其頻譜。程序清單如下Fs

59、=10;xn=1,2,3,2,1;N=length(xn);D=2*pi*Fs/N;k=floor(-(N-1)/2:(N-1)/2);X=fftshift(fft(xn,N);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é)果如圖4-10所示：圖4-10由圖4-10可知，當(dāng)有限長序列的長度N=5時，頻譜的樣本點數(shù)也為5，頻率點之間的間距非常大，即分辨率很低。即使使用了plot命令的插值

60、功能，顯示出的曲線仍是斷續(xù)的，與真實曲線有較大誤差。改變分辨率的基本方法是給輸入序列補零，即增加頻譜的密度。這種方法只是改善了圖形的視在分辨率，并不增加頻譜的細節(jié)信息。將上述有限長序列x(n)1,2,3,2,1末尾補零到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);X=fftshift(fft(xn,N);subplot(1,2,1);plot(k*D,abs(X);title(幅度頻譜);xlabe