離散線性時不變系統(tǒng)分析

1、實驗六 離散線性時不變系統(tǒng)分析一、 實驗目的1. 掌握離散LSI系統(tǒng)的單位序列響應(yīng)、單位階躍響應(yīng)和任意激勵下響應(yīng)的MATLAB求解方法。2.掌握離散LSI系統(tǒng)的頻域分析方法；3.掌握離散LSI系統(tǒng)的復頻域分析方法；4. 掌握離散LSI系統(tǒng)的零極點分布與系統(tǒng)特性的關(guān)系。二、實驗原理及方法1. 離散LSI系統(tǒng)的時域分析描述一個N階線性時不變離散時間系統(tǒng)的數(shù)學模型是線性常系統(tǒng)差分方程，N階LSI離散系統(tǒng)的差分方程一般形式為 （6.1）也可用系統(tǒng)函數(shù)來表示 （6.2）系統(tǒng)函數(shù)反映了系統(tǒng)響應(yīng)和激勵間的關(guān)系。一旦上式中，的數(shù)據(jù)確定了，系統(tǒng)的性質(zhì)也就確定了。特別注意必須進行歸一化處理，即。對于復雜信號激勵下

2、的線性系統(tǒng)，可以將激勵信號在時域中分解為單位序列或單位階躍序列的線性疊加，把這些單元激勵信號分別加于系統(tǒng)求其響應(yīng)，然后把這些響應(yīng)疊加，即可得到復雜信號作用于系統(tǒng)的零狀態(tài)響應(yīng)。因此，求解系統(tǒng)的單位序列響應(yīng)和單位階躍響應(yīng)尤為重要。由圖6-1可以看出一個離散LSI系統(tǒng)響應(yīng)與激勵的關(guān)系。圖6-1 離散LSI系統(tǒng)響應(yīng)與激勵的關(guān)系(1) 單位序列響應(yīng)（單位響應(yīng)）單位響應(yīng)是指離散LSI系統(tǒng)在單位序列激勵下的零狀態(tài)響應(yīng)，因此滿足線性常系數(shù)差分方程(6.1)及零初始狀態(tài)，即， (6.3)按照定義，它也可表示為 (6.4)對于離散LSI系統(tǒng)，若其輸入信號為，單位響應(yīng)為，則其零狀態(tài)響應(yīng)為 (6.5)可見，能夠刻畫和

3、表征系統(tǒng)的固有特性，與何種激勵無關(guān)。一旦知道了系統(tǒng)的單位響應(yīng)，就可求得系統(tǒng)對任何輸入信號所產(chǎn)生的零狀態(tài)響應(yīng)。MATLAB提供了專門用于求連續(xù)系統(tǒng)沖激響應(yīng)的函數(shù)impz()，其調(diào)用格式有h,n=impz(b,a)求解離散系統(tǒng)的單位響應(yīng)，其中，；h,n=impz(b,a,N) 求解離散系統(tǒng)的單位響應(yīng)，采樣點數(shù)由N確定，；impz(b,a) ：在當前窗口，用stem(n,h)繪出圖形。（2）單位階躍響應(yīng) 單位階躍響應(yīng)是指離散LTI系統(tǒng)在單位階躍序列激勵下的零狀態(tài)響應(yīng)，它可以表示為 (6.6)上式表明，離散LSI系統(tǒng)的單位階躍響應(yīng)是單位響應(yīng)的累加和，系統(tǒng)的單位階躍響應(yīng)和系統(tǒng)的單位響應(yīng)之間有著確定的關(guān)系

4、，因此，單位階躍響應(yīng)也能完全刻畫和表征一個LSI系統(tǒng)。MATLAB提供了專門用于求離散系統(tǒng)單位階躍響應(yīng)的函數(shù)stepz( )，其調(diào)用格式有s,n=stepz(b,a) ：求解離散系統(tǒng)的單位階躍響應(yīng)，其中，；s,n=stepz(b,a,N) ：求解離散系統(tǒng)的單位階躍響應(yīng)，采樣點數(shù)由N確定，；stepz(b,a) ：在當前窗口，用stem(n,s)繪出圖形。(3) 任意激勵下的零狀態(tài)響應(yīng)已經(jīng)知道，離散LSI系統(tǒng)可用常系數(shù)線性差分方程(6.1)式來描述，Matlab提供的函數(shù)dlsim( )能對上述差分方程描述的離散LSI系統(tǒng)的響應(yīng)進行仿真，該函數(shù)不僅能繪制指定時間范圍內(nèi)的系統(tǒng)響應(yīng)波形圖，而且還能求

5、出系統(tǒng)響應(yīng)的數(shù)值解。其調(diào)用格式有dlsim(b,a, x) ：求解輸入序列為x的零狀態(tài)響應(yīng)需要特別強調(diào)的是，Matlab總是把由分子和分母多項式表示的任何系統(tǒng)都當作是因果系統(tǒng)。所以，利用impz (b,a)，stepz(b,a)，dlsim(b,a,x)函數(shù)求得的響應(yīng)總是因果信號。同時，卷積和也是LSI系統(tǒng)求解零狀態(tài)響應(yīng)的重要工具之一。假設(shè)系統(tǒng)的輸入信號為，單位響應(yīng)為，則系統(tǒng)的零狀態(tài)響應(yīng)可由（6.5）式求解。Matlab提供了專門用于求離散系統(tǒng)卷積和的函數(shù)conv( )，其調(diào)用格式有y=conv(x,h) ：求解序列x，h的卷積和，若序列x的長度為n1，序列h的長度為n2，卷積和y的長度為n1

6、+n2-1。這一點需要特別注意，否則，作圖時容易造成橫縱坐標長度不匹配。（4）帶初始狀態(tài)的任意激勵下的全響應(yīng)任意激勵下的離散LSI系統(tǒng)的全響應(yīng)為零輸入響應(yīng)和零狀態(tài)響應(yīng)之和，表示為 （6.7）在理論學習的過程中，同學們對低階差分方程的求解已頗為頭痛，高階差分方程直接求解幾乎不可能。Matlab提供了用于求離散系統(tǒng)全響應(yīng)的函數(shù)filter( )，其調(diào)用格式有y=filter( b,a,x) ：求解零狀態(tài)響應(yīng)；y=filter( b,a,x,zi) ：求解初始條件為zi的系統(tǒng)的全響應(yīng)，zi向量的長度為max(length(a),length(b)-1，返回值為系統(tǒng)的全響應(yīng)。z = filtic(b,

7、a,y,x)：將初始狀態(tài)轉(zhuǎn)換為初始條件，其中，；z = filtic(b,a,y)：將初始狀態(tài)轉(zhuǎn)換為初始條件，其中，。2 離散LSI系統(tǒng)的復頻域（Z域）分析（1）利用Z變換解差分方程在前面圖6-1中表示了離散系統(tǒng)的響應(yīng)與激勵的關(guān)系，由圖可知，系統(tǒng)的響應(yīng)既可以用時域的方法求解，也可以用Z域的方法求解。當已知系統(tǒng)輸入序列的Z變換，系統(tǒng)函數(shù)時，系統(tǒng)響應(yīng)序列的Z變換可由求出。Matlab提供了用于求序列Z變換和Z反變換的函數(shù)，其調(diào)用格式有X=ztrans(x)：求無限長序列x的Z變換，返回Z變換的表達式，注意這里x，X都是符號表達式；x=iztrans(X)：求X（z）的Z反變換x(n)，返回Z反變換

8、的表達式，注意這里x，X都是符號表達式；r,p,c=residuez(b,a)：把b(z)/a(z)展開成部分分式；b,a=residuez(r,p,c):根據(jù)部分分式的r、p、c數(shù)組，返回有理多項式。（2）系統(tǒng)的零極點分布與系統(tǒng)因果性和穩(wěn)定性的關(guān)系因果系統(tǒng)的單位響應(yīng)一定滿足當時，=0，那么其系統(tǒng)函數(shù)的收斂域一定包含點，即點不是極點，極點分布在某個圓的圓內(nèi)，收斂域是圓外區(qū)域。系統(tǒng)穩(wěn)定要求，對照z變換定義，系統(tǒng)穩(wěn)定要求收斂域包含單位圓。如果系統(tǒng)因果且穩(wěn)定，收斂域包含點和單位圓，那么收斂域可表示為：(6.8)MATLAB提供了用于求系統(tǒng)零極點的函數(shù)，其調(diào)用格式有roots()：利用多項式求根函數(shù)來

9、確定系統(tǒng)函數(shù)的零極點位置；roots(a)：求極點位置，a為系統(tǒng)函數(shù)分母多項式所構(gòu)成的系數(shù)向量；roots(b)：求零點位置，b為系統(tǒng)函數(shù)分子多項式所構(gòu)成的系數(shù)向量；zplane(b,a)：繪制由行向量b和a構(gòu)成的系統(tǒng)函數(shù)的零極點分布圖；zplane(z,p)：繪制由列向量z確定的零點、列向量p確定的極點構(gòu)成的零極點分布圖。（3）系統(tǒng)的零極點分布與系統(tǒng)頻率響應(yīng)的關(guān)系將式(6.2)因式分解，得到 (6.9)式中，是的零點，是其極點。參數(shù)影響頻率響應(yīng)的幅度大小，影響系統(tǒng)特性的是零點和極點的分布。下面采用幾何方法研究系統(tǒng)零極點分布對系統(tǒng)頻率特性的影響。將式(6.9)的分子、分母同乘以，得到： (6.

10、10)假設(shè)系統(tǒng)穩(wěn)定，將代入上式，得到頻率響應(yīng) (6.11)設(shè)，由式(6.11)得到 (6.12)在平面上，用一根由零點指向單位圓（）上任一點B的向量表示，同樣用一根由極點指向單位圓（）上任一點B的向量來表示，和分別稱為零點矢量和極點矢量，用極坐標表示為： ，。將和表示式代入式(6.12)，得到 (6.13) (6.14)系統(tǒng)或者信號的頻率特性由式(6.13)和式(6.14)確定。按照式(6.13)，知道零極點的分布后，可以很容易地確定零極點位置對系統(tǒng)特性的影響。當點轉(zhuǎn)到極點附近時，極點矢量長度最短，因而幅度特性可能出現(xiàn)峰值，且極點愈靠近單位圓，極點矢量長度愈短，峰值愈高愈尖銳。如果極點在單位圓

11、上，則幅度特性為，系統(tǒng)不穩(wěn)定。對于零點，情況相反，當點轉(zhuǎn)到零點附近時，零點矢量長度變短，幅度特性將出現(xiàn)谷值，且零點愈靠近單位圓，谷值愈接近零。當零點在單位圓上時，谷值為零。綜上所述，極點位置主要影響頻響的峰值位置及尖銳程度，零點位置主要影響頻響的谷值位置及形狀。Matlab提供了專門用于求離散系統(tǒng)頻率響應(yīng)的函數(shù)freqz(),其調(diào)用格式如下：H,w = freqz(b,a,n)：返回數(shù)字系統(tǒng)的n點頻率值（復數(shù)），這n個點均勻地分布在0，上，系統(tǒng)默認的采樣點數(shù)目為512點；H,f = freqz(b,a,n,Fs)：用于對在0，F(xiàn)s/2上等間隔采樣n點，采樣點頻率及相應(yīng)的頻響值分別存放在f和H中

12、。H = freqz(b,a,w)：用于對在0，2上進行采樣，采樣頻率點由w指定。H = freqz(b,a,f,Fs)：用于對在0，F(xiàn)s上進行采樣，采樣頻率點由f指定。freqz(b,a,n)：用于在當前圖形窗口中繪制幅頻和相頻特性曲線。下面介紹幾個實用的函數(shù)：mag=abs(H):求解系統(tǒng)的絕對幅頻響應(yīng)；db=20*log10(mag+eps)/max(mag)：求解系統(tǒng)的相對幅頻響應(yīng)；pha=angle(H): 求解系統(tǒng)的相頻響應(yīng)；grd=grpdelay(b,a,w): 求解系統(tǒng)的群延遲。三、實驗內(nèi)容及步驟1. 某離散LSI系統(tǒng)的差分方程表示式為滿足初始條件，求系統(tǒng)的單位響應(yīng)，單位階躍

13、響應(yīng)，用filter子函數(shù)求系統(tǒng)輸入為時的零輸入、零狀態(tài)及全響應(yīng)。提示：通過解差分方程，可以得到全響應(yīng)為，使用filter子函數(shù)對系統(tǒng)差分方程進行求解，同時將求解結(jié)果與理論計算的結(jié)果進行比較。2. 一簡單的數(shù)字微分器它計算輸入序列的后向一階差分。當數(shù)字微分器的輸入為以下序列時，計算、畫出輸出序列，并對該數(shù)字微分器的實用性進行評價。（1）矩形脈沖：（2）三角脈沖：（3）正弦脈沖：3. 已知離散系統(tǒng)的系統(tǒng)函數(shù)為求該系統(tǒng)的零極點及零極點分布圖，并判斷系統(tǒng)的因果性和穩(wěn)定性。4. 已知離散系統(tǒng)的系統(tǒng)函數(shù)為繪出系統(tǒng)的零極點分布圖，系統(tǒng)在頻率范圍內(nèi)的絕對幅頻響應(yīng)、相對幅頻響應(yīng)、相位頻率響應(yīng)和群延遲，并觀察零

14、極點分布與系統(tǒng)幅頻響應(yīng)的關(guān)系。5. 觀察系統(tǒng)零極點的位置對幅頻響應(yīng)的影響。已知一階離散系統(tǒng)的系統(tǒng)函數(shù)為，（1）假設(shè)系統(tǒng)的零點在原點，極點分別取0.2、0.5、0.8，比較它們的幅頻響應(yīng)曲線，（2）假設(shè)系統(tǒng)的極點在原點，零點分別取0.2、0.5、0.8，比較它們的幅頻響應(yīng)曲線，從中總結(jié)零極點位置對幅頻響應(yīng)的影響。二、 實驗處理：3.1 MATLAB源代碼：系統(tǒng)的單位沖激響應(yīng)clearb=1 ,-1.5 ,0.5;a=1;N=20;n=1:N;yim=impz(b,a,N);stem(n,yim,.);單位階躍響應(yīng)clearb=1 ,-1.5 ,0.5;a=1;N=20;m=1:N;n=m;yst

15、=stepz(b,a,N);stem(n,yst,.);零輸入、零狀態(tài)及全響應(yīng)cleara=1 ,-1.5 ,0.5;b=1;N=20;n=0:N-1;x0=zeros(1,N);y0=4,10;x=0.25.n;zi=filtic(b,a,y0);yzi=filter(b,a,x0,zi);zi1=filtic(b,a,0);yzs=filter(b,a,x,zi1);y=filter(b,a,x,zi);yy=0.5.n+(1/3)*(1/4).n+2/3;subplot(2,3,1),stem(n,x,.); title(輸入信號);subplot(2,3,2),stem(n,yzi,.

16、);title(系統(tǒng)的零輸入響應(yīng));subplot(2,3,3),stem(n,yzs,.);title(系統(tǒng)的零狀態(tài)響應(yīng));subplot(2,3,4),stem(n,y,.);title(系統(tǒng)的全響應(yīng));subplot(2,3,5),stem(n,yy,.);title(理論計算的系統(tǒng)全響應(yīng));3.2 MATLAB源程序clearN1=30;N2=30;N3=120;n1=0:N1;n2=0:N2;n3=0:N3;y1=x1(n1)-x1(n1-1);y2=x2(n2)-x2(n2-1);y3=x3(n3)-x3(n3-1);subplot(3,1,1);stem(n1,y1,.);sub

17、plot(3,1,2);stem(n2,y2,.);subplot(3,1,3);stem(n3,y3,.);其中函數(shù)u(n)、x1(n)、x2(n)分別定義為：function y=u(n)y=(n0);function y1=x1(n)y1=5.*(u(n)-u(n-20);function y2=x2(n)y2=n.*(u(n)-u(n-10)+(20-n).*(u(n-10)-u(n-20);function y3=x3(n)y3=sin(pi*n/25).*(u(n)-u(n-100);3.3MATLAB源程序：b=0.2 0.1 0.3 0.1 0.2;a=1 -1.1 1.5 -

18、0.7 0.3;z=roots(b) p=roots(a) subplot(221),zplane(b,a);title(系統(tǒng)的零極點分布圖);subplot(223),impz(b,a,20);title(系統(tǒng)的單位響應(yīng));z =p = -0.5000 + 0.8660i0.2367 + 0.8915i -0.5000 - 0.8660i0.2367 - 0.8915i 0.2500 + 0.9682i0.3133 + 0.5045i 0.2500 - 0.9682i0.3133 - 0.5045i系統(tǒng)極點在單位圓內(nèi)，零點在圓上，為梳妝濾波器，系統(tǒng)穩(wěn)定。3.4MATLAB源程序：b=0.1

19、0 0.4 0 0.4 0 0.1;a=1 0 -0.3 0 0.6 0 -0.2;z=roots(b) p=roots(a) zplane(b,a);title(零極點分布圖);n=(0:500)*pi/500;figure(2)H,w=freqz(b,a,n);subplot(221),plot(w/pi,abs(H);grid axis(0 1 1.1*min(abs(H) 1.1*max(abs(H);title(絕對幅頻響應(yīng));subplot(222),plot(w/pi,angle(H);grid axis(0 1 1.1*min(angle(H) 1.1*max(angle(H)

20、;title(相頻響應(yīng));db=20*log10(abs(H); subplot(223),plot(w/pi,db);gridaxis(0 1 -100 5);title(相對幅頻響應(yīng)（dB）);grd=grpdelay(b,a,w); subplot(224), plot(w/pi,grd);gridtitle(群延遲);z = 0 + 1.6180i 0 - 1.6180i 0 + 1.0000i 0 - 1.0000i 0 + 0.6180i 0 - 0.6180ip = 0.6190 + 0.6304i 0.6190 - 0.6304i 0.5729 -0.6190 + 0.6304

21、i -0.6190 - 0.6304i -0.5729 3.5.1MATLAB源程序：b=0 1;m=0.2 0.5 0.8;for i=1:3a=1 -m(i);figure(i) z=roots(b) p=roots(a) subplot(311),zplane(b,a);title(零極點分布圖);n=(0:500)*pi/500;H,w=freqz(b,a,n);subplot(312),plot(w/pi,abs(H);grid axis(0 2 1.1*min(abs(H) 1.1*max(abs(H);title(絕對幅頻響應(yīng));db=20.*log10(abs(H); subp

22、lot(313),plot(w/pi,db);gridaxis(0 1 -10 6);title(相對幅頻響應(yīng)（dB）);endz = Empty matrix: 0-by-1p = 0.2000z = Empty matrix: 0-by-1p = 0.5000z = Empty matrix: 0-by-1p = 0.8000 MATLAB源程序：a=0 1;b=1 -0.2;/分別換做0.5、0.8z=roots(b); p=roots(a); subplot(211);zplane(a,b);title(零極點分布圖);n=(0:500)*pi/500;H,w=freqz(b,a,n);subplot(212),plot(w/pi,abs(H);grid axis(0 2 1.1*min(abs(H) 1.1*max(a