第2章離散時間系統(tǒng)的時域分析2.1相關MATLAB函數(shù)

2.2用遞歸法求解差分方程

2.3卷積和系統(tǒng)零狀態(tài)響應的計算

2.4用MATLAB求離散系統(tǒng)的單位階躍響應和單位樣值響應

2.5用MATLAB求離散系統(tǒng)在任意輸入作用下的響應

2.6線性時不變系統(tǒng)的實例分析習題本章介紹離散時間系統(tǒng)的時域分析以及如何用MATLAB函數(shù)計算離散線性時不變系統(tǒng)的響應。2.1節(jié)簡要介紹離散時間系統(tǒng)時域分析用到的MATLAB函數(shù)。2.2節(jié)介紹用遞歸法求解差分方程。2.3節(jié)介紹卷積的計算和用卷積求解離散系統(tǒng)的零狀態(tài)響應。2.4節(jié)介紹用MATLAB求離散系統(tǒng)的單位階躍響應和單位樣值響應。2.5節(jié)介紹用MATLAB求離散系統(tǒng)對任意輸入的響應。2.6節(jié)為線性時不變系統(tǒng)實例分析。

離散時間系統(tǒng)的時域分析涉及到的MATLAB函數(shù)主要有：

conv

dstep

dimpulse

filter

下面簡要介紹這些函數(shù)。2.1相關MATLAB函數(shù)

1.conv

函數(shù)conv用來計算離散時間信號的卷積x(n)*h(n),用向量x和h分別定義序列x(n)和h(n)的元素,那么可以用命令

y＝conv(x,h)

計算這兩個序列的卷積。

函數(shù)conv還可以用于計算多項式的乘法。假設多項式a(s)的系數(shù)由向量a給出,多項式b(s)的系數(shù)由向量b給出,那么多項式a(s)b(s)的系數(shù)可以通過命令

ab＝conv(a,b)

賦值給向量ab。

2.dstep

函數(shù)dstep用于計算離散時間系統(tǒng)的單位階躍響應。

3.dimpulse

函數(shù)dimpulse用于計算離散時間系統(tǒng)的單位樣值響應。

4.filter

函數(shù)filter用于計算離散系統(tǒng)在任意輸入信號作用下的響應。

線性時不變離散時間系統(tǒng)可以用差分方程描述：

(2.1)

2.2用遞歸法求解差分方程式中:n是整數(shù),x(n)是輸入,y(n)是輸出,系數(shù)ai(i＝0,1,…,N)和bi(i＝0,1,…,M)均為常數(shù)。一般情況下,a0＝1,方程(2.1)可以寫做

(2.2)

差分方程可以直接利用遞歸的方法來求解。將式(2.2)改寫為

(2.3)從式(2.3)可以看出,輸出的下一個值是前N個輸出值和M＋1個輸入值的線性組合。例如,y(0)是y(－1),y(－2),…,y(－N)和x(0),x(－1),…,x(－M)的線性組合。由式(2.2)或式(2.3)確定的離散時間系統(tǒng)也被稱為遞歸離散時間系統(tǒng)。

需要注意的是,式(2.3)中至少有一個系數(shù)ai(i＝1,2,…,N)是非零的。如果所有的ai為零,那么式(2.3)可以簡化為

(2.4)用MATLAB實現(xiàn)遞歸的方法時,把式(2.3)的系數(shù)ai(i＝1,2,…,N)和bi(i＝0,1,…,M)存在向量a＝[a1a2…aN］和b＝[b0b1b2…bM］中；輸入信號x(n)存在向量x中；輸入x(n)初始值和輸出y(n)初始值存在向量x0＝[x[n0－M］,x[n0－M＋1］,…,x[n0－1］］和y0＝[y[n0－N］,y[n0－N＋1］,…,y[n0－1］］中；計算需要迭代的時間向量用n表示,這里n0表示向量n的第一個元素,如果迭代從n＝0開始,那么n0＝0；迭代輸出的值保存在向量y中,y的第一個元素對應時間向量中的n0。迭代可以從任何時間開始。式(2.3)中的和式可以寫成矩陣相乘的形式：

(2.5)

根據(jù)以上分析,可以把遞歸的計算方法寫成一個MATLAB函數(shù)：

functiony＝recur(a,b,n,x,x0,y0)；

N＝length(a)；

M＝length(b)－1；

y＝[y0zeros(1,length(n］；

x＝[x0x］；

a1＝a(length(a)∶－1∶1)；

b1＝b(length(b)∶－1∶1)；

fori＝N＋1∶N＋length(n),

y(i)＝－a1*y(i－N∶i－1)′＋b1*x(i－N∶i－N＋M)′；

end

y＝y(tǒng)(N＋1∶N＋length(n；

這個函數(shù)保存在名為recur.m的M文件中,在以后的計算中可以調用。例2-1用遞歸法求解差分方程

在0≤n≤5時的解,其中輸入信號 解：這里我們調用recur函數(shù)來求解。設方程系數(shù)向量a＝[－0.8］,b＝[2］。這是一個一階差分方程,輸出信號只有一個初始值為y(－1)＝2,向量y0＝[2］。因為輸入信號為x(n)的形式,所以不需要添加輸入信號在n<n0時的初始值,向量x0＝[］是一個空向量,但是在調用recur函數(shù)時這個向量需要定義。MATLAB程序如下：。

a＝[－0.8］；b＝[2］；

x0＝[］；y0＝[2］；

n＝0∶5；

x＝n；

y＝recur(a,b,n,x,x0,y0)

stem(n,y)

xlabel(′n′)

ylabel(′y(n)′)

得到y(tǒng)的輸出結果為：

y＝1.60003.28006.624011.299217.039423.6315

畫出輸出的波形如圖2-1所示。

圖2-1例2-1系統(tǒng)輸出的波形

離散信號x(n)和h(n)的卷積定義為

(2.6)

這個計算實際上是一個求和過程,也叫卷積和。2.3卷積和系統(tǒng)零狀態(tài)響應的計算如果信號x(n)和h(n)都是有限長的序列,設信號x(n)的定義范圍是

(2.7)

信號h(n)的定義范圍是

(2.8)

那么y(n)的定義范圍是

(2.9)

即用于保存y(n)的向量y包含Nx＋Nh－1個樣本值。在用MATLAB函數(shù)conv計算離散信號的卷積時,要注意信號x(n)和h(n)的時間起點。例如:

x＝[11111］

h＝[0123］

y＝[01366653］

如果x(n)和h(n)的時間起點都是n＝0,那么y(n)的時間起點也是n＝0。如果x(n)的時間起點是n＝－2,h(n)的時間起點是n＝－3,那么y(n)的時間起點是n＝－5。在畫出卷積結果的波形時要注意調整標號向量,使之與信號的值對應。用MATLAB計算無限長的序列的卷積時,只能截取有限長度的信號值,因此卷積得到的答案等效于信號在定義范圍外的值為零的計算結果。例如,信號x(n)的定義范圍是0≤n≤q,信號h(n)定義在0≤n≤r,則y(n)＝x(n)*h(n)的定義范圍是0≤n≤q＋r,但是只有在0≤n≤min(q,r)范圍內(nèi)的答案才是正確的。例2-2設信號x(n)和h(n)都是有限長信號:

用MATLAB計算 和

,并畫出y1(n)和y2(n)的波形。

解：根據(jù)式(2.9)的結論,信號y1(n)的范圍是0≤n≤10,信號y2(n)的范圍是－5≤n≤5,MATLAB程序如下：

n＝0∶5；

x＝ones(1,6)；

h＝n；

y＝conv(x,h)；

ny＝0∶10；

subplot(2,1,1)；

stem(ny,y)；

xlabel(′n′)；ylabel(′y_1(n)′)；title(′y_1(n)＝x(n)*h(n)′)

y2＝conv(x,h)；

ny2＝ny－5；

subplot(2,1,2)；

stem(ny2,y2)；

xlabel(′n′)；ylabel(′y_2(n)′)；title(′y_2(n)＝x[n］*h[n＋5］′)

運行的結果如圖2-2所示。

圖2-2y1(n)和y2(n)的波形圖2-2中,圖(a)為 的波形,圖(b)為y2(n)＝x(n)*h(n＋5)的波形。

例2-3設信號x(n)＝u(n),

(1)計算y(n)＝x(n)*x(n)的解析式；

(2)取信號x(n)＝u(n)(0≤n≤5),用MATLAB計算y(n),畫出y(n)在0≤n≤10的波形y1(n)和y(n)在0≤n≤5的波形y2(n),并進行比較。解：

(2)MATLAB程序如下：

x＝ones(1,6)；h＝x；

y＝conv(x,h)；

n＝0∶10；

subplot(2,1,1)；

stem(n,y(1∶length(n)

axis([－21006］)；xlabel(′n′)；ylabel(′y_1(n)′)；title(′y_1(n)＝x(n)*x(n)′)

n1＝0∶5；

subplot(2,1,2)；

stem(n1,y(1∶length(x)

axis([－21006］)；xlabel(′n′)；ylabel(′y_2(n)′)；title(′y_2(n)＝u(n)*u(n)′)

波形如圖2-3所示。

圖2-3中,圖(a)的波形在0≤n≤5范圍內(nèi)符合解析式的正確結果,從n>5開始信號的值逐漸減小,出現(xiàn)這個結果是因為向量x和y的長度不同,向量y的后5個樣本值呈下降趨勢是由計算中認為信號x(n)在n>5后的值為零造成的,這個卷積的過程可以通過圖解法進行說明。因此,一般在表示卷積波形時,像圖(b)一樣截取卷積結果的前6個樣本值來反映波形的趨勢。用卷積可以計算離散系統(tǒng)的零狀態(tài)響應。用這種方法計算響應時,首先要知道系統(tǒng)的單位樣值響應h(n),單位樣值響應與系統(tǒng)激勵的卷積就是系統(tǒng)的零狀態(tài)響應。用MATLAB求解離散系統(tǒng)單位樣值響應的問題在下一節(jié)介紹。用MATLAB求解零狀態(tài)響應時,需要注意輸出結果的取值范圍。例2-4用卷積求解差分方程

的零狀態(tài)響應,其中輸入信號x(n)＝nu(n),取x(n)的范圍是0≤n≤5。

解：用解析法求得系統(tǒng)的單位樣值響應是h(n)＝2(0.8)nu(n),用卷積求解系統(tǒng)零狀態(tài)響應的MATLAB程序如下:

n＝0∶5；

x＝n；

h＝2*(0.8).^n；

y＝conv(x,h)

stem(n,y(1∶length(n)

xlabel(′n′)；ylabel(′y_{zs}(n)′)；title(′y_{zs}(n)＝x(n)*h(n)′)

得到y(tǒng)的輸出結果為

y＝02.00005.600010.480016.384023.107218.485814.264310.36296.71743.2768

因為x(n)的范圍取0≤n≤5,所以零狀態(tài)響應取y的前6個樣本值。系統(tǒng)零狀態(tài)響應的波形如圖2-4所示。

圖2-4系統(tǒng)零狀態(tài)響應的波形

單位階躍響應是單位階躍序列作用在系統(tǒng)產(chǎn)生的零狀態(tài)響應；單位樣值響應是單位樣值序列作用在系統(tǒng)產(chǎn)生的零狀態(tài)響應。MATLAB中可以用dstep函數(shù)求系統(tǒng)的單位階躍響應,用dimpulse函數(shù)求系統(tǒng)的單位樣值響應。用這兩個函數(shù)求響應時,先把式(2.1)方程左右兩邊多項式的系數(shù)ai(i＝0,1,…,N)和bi(i＝0,1,…,M)用向量a和b表示,向量a表示與響應有關的系數(shù)：2.4用MATLAB求離散系統(tǒng)的單位階躍響應和單位樣值響應

a＝[a0a1a2…aN］

向量b表示與激勵有關的系數(shù)：

b＝[b0b1b2…bM0…0］

如果M<N,需要在向量b已有的系數(shù)后面補充N－M個0,然后用命令y＝dstep(b,a,n)或y＝dimpulse(b,a,n)求系統(tǒng)的單位階躍響應或單位樣值響應,命令中n表示需要計算的樣值點數(shù)。這里之所以這樣處理向量a和b,是因為這兩個函數(shù)中向量a和b是系統(tǒng)差分方程通過)Z)變換得到的系統(tǒng)函數(shù)H(z)的分母多項式和分子多項式的系數(shù),并且這兩個多項式是按照z的降冪排列的。例2-5求解差分方程 的單位樣值響應和單位階躍響應。

解：MATLAB程序如下：

n＝0∶10；

a＝[1－0.8］；b＝[20］；

ys＝dstep(b,a,length(n；

yi＝dimpulse(b,a,length(n；

subplot(2,1,1)；stem(n,ys)；

xlabel(′n′)；ylabel(′ys(n)′)；title(′StepResponse′)

subplot(2,1,2)；stem(n,yi)；

xlabel(′n′)；ylabel(′yi(n)′)；title(′ImpulseResponse′)

求出響應的波形如圖2-5所示。

圖2-5系統(tǒng)單位階躍響應和單位樣值響應

(a)單位階躍響應；(b)單位樣值響應

用MATLAB的filter函數(shù)可以求線性時不變離散系統(tǒng)在任意輸入作用下的響應。如果系統(tǒng)差分方程如式(2.1)所示,系數(shù)ai(i＝0,1,…,N)和bi(i＝0,1,…,M)存在向量a＝[a0a1a2…aN］和b＝[b0b1b2…bM］中,區(qū)間nx≤n≤nx＋Nx－1內(nèi)的輸入信號x(n)用向量x表示,一般nx＝0,設系統(tǒng)為零狀態(tài),那么命令y＝filter(b,a,x)就得到系統(tǒng)的輸出。用這個命令產(chǎn)生的輸出向量y包含的樣本區(qū)間的范圍與向量x定義的區(qū)間范圍是一樣的。2.5用MATLAB求離散系統(tǒng)在任意輸入作用下的響應如果系統(tǒng)為非零狀態(tài)的,filter函數(shù)在計算出向量y的第1個輸出值y(nx)時,需要輸入信號x(n)在nx－M≤n≤nx－1范圍內(nèi)的值,以及y(n)在nx－N≤n≤nx－1范圍內(nèi)的值,即系統(tǒng)的初始狀態(tài)。如果不給出這些值,filter函數(shù)假設這些樣本值全都為零。因此，在求有初始狀態(tài)的系統(tǒng)的響應時,用命令y＝filter(b,a,x,zi)得到差分方程的輸出,其中,向量zi表示系統(tǒng)的初始狀態(tài),要用狀態(tài)方程獲得?；趂ilter函數(shù)的這些特征,它不僅可以用于計算離散系統(tǒng)在任意輸入作用下的零狀態(tài)響應,還可以用于計算離散系統(tǒng)在任意輸入作用下的非零狀態(tài)的響應,即全響應。

例2-6離散系統(tǒng)的差分方程為y(n)－0.8y(n－1)＝2x(n),

(1)輸入信號x(n)＝nu(n),取x(n)的范圍是0≤n≤5,用filter函數(shù)求解系統(tǒng)的零狀態(tài)響應。

(2)輸入信號同(1),y(－1)＝2,用filter函數(shù)求解系統(tǒng)的全響應。解:

(1)求系統(tǒng)零狀態(tài)響應的MATLAB程序如下：

n＝0∶5；

a＝[1－0.8］；b＝[2］；

x＝n；

yzs＝filter(b,a,x)

得到y(tǒng)zs的輸出結果為

yzs＝02.00005.600010.480016.384023.1072

如果畫出輸出的波形,同圖2-4所示的一樣。

(2)求系統(tǒng)全響應的MATLAB程序如下：

n＝0∶5；

a＝[1－0.8］；b＝[2］；

x＝n；

zi＝[0.8*2］；

y＝filter(b,a,x,zi)

得到y(tǒng)的輸出結果為

y＝1.60003.28006.624011.299217.039423.6315

如果畫出輸出的波形,同圖2-1所示的一樣。

(2)題中給出的是一個非零初始狀態(tài)的系統(tǒng),在用filter函數(shù)時,需要給出初始狀態(tài)向量zi。這個系統(tǒng)是一階的,又因為n<0時,x(n)＝0,所以zi＝[－a1*y[－1］］,即zi＝[0.8*2］。

例2-7兔子的繁殖問題。

研究一對兔子的繁殖問題,假設每對兔子每月可以生育一對小兔,新生的小兔子要隔一個月才具有生育能力,若第一個月只有一對新生的小兔子,求第n個月兔子對的數(shù)目。2.6線性時不變系統(tǒng)的實例分析

(1)建立描述兔子繁殖問題的差分方程；

(2)用遞歸的方法求解方程,獲得系統(tǒng)的解；

(3)用經(jīng)典法求出差分方程解的表達式；

(4)比較兩種方法的結果。

解:

(1)用y(n)表示在第n個月兔子對的數(shù)目,已知y(0)＝0,y(1)＝1,y(2)＝1,y(3)＝2,…。容易得到,到第n個月時,應有y(n－2)對兔子具有生育能力,因而這批兔子會從y(n－2)對變成2y(n－2)對；此外,還有y(n－1)－y(n－2)對兔子沒有生育能力,因為它們是在第(n－1)月新生的,仍按原數(shù)目計算,可以寫出系統(tǒng)方程

y(n)＝2y(n－2)＋[y(n－1)－y(n－2)］

整理一下得

y(n)－y(n－1)－y(n－2)＝0 (2.10)

(2)將式(2.10)寫成式(2.3)的形式:

y(n)＝y(tǒng)(n－1)＋y(n－2) (2.11)

得到系統(tǒng)的遞歸關系。

(3)系統(tǒng)的差分方程是一個齊次方程,由經(jīng)典法得方程的解為

(2.12)

(4)MATLAB程序如下：

a＝[－1－1］；b＝[0］；

n＝2∶10；

x＝zeros(1,length(n；

x0＝[］；y0＝[01］；

y1＝recur(a,b,n,x,x0,y0)

y2＝(((1＋sqrt(5/2).^n－((1－sqrt(5/2).^n)/sqrt(5)

y1＝[y0y1］；y2＝[y0y2］；

n＝0∶10；

plot(n,y1,′o′,n,y2,′＋′)

title(′rabbit′)

ylabel(′y(n)′)

xlabel(′n′)

text(1,35,′o用遞歸方法求解的結果′)

text(1,30,′＋用經(jīng)典法求解的結果′)在程序中,分別用y1和y2表示用遞歸方法和解析式求得的系統(tǒng)響應向量,結果如下：

y1＝1235813213455

y2＝1.00002.00003.00005.00008.000013.000021.0000

34.000055.0000

這兩種方法的結果是相同的。畫出響應結果的圖形如圖2-6所示。

圖2-6兔子的繁殖例2-8回聲的消除。

假設有一個拾音器放在一間屋子的當中,一個人在房間的另一頭說話,錄音中將包含直接傳到拾音器的語音,以及經(jīng)過傳輸,從遠處墻壁反射回來再進入拾音器的回聲。因為回聲必須是向前傳播的,所以時間上有一定的延時；另外,因為語音有一部分被墻壁吸收,所以幅度上有一定的衰減。為了簡化分析,只考慮一次反射的回聲源。差分方程

y(n)＝x(n)＋αx(n－N) (2.13)

是一個回聲產(chǎn)生模型,其中,x(n)是未被污損的語音信號,y(n)是疊加上回聲以后的輸出信號。設延時值N＝1000,回聲衰減系數(shù)α＝0.5。

這個實例用的原始語音存儲在文件lineup.mat中,在著手分析系統(tǒng)前,需要用load命令裝入這個文件。語音波形數(shù)據(jù)存儲在向量y中。

(1)求解回聲系統(tǒng)的單位樣值響應,將它在0≤n≤4000內(nèi)的值存入向量he中,并畫出該響應的波形。

(2)計算回聲系統(tǒng)的系統(tǒng)函數(shù)H(z)。

(3)設計一個回聲消除系統(tǒng),能夠從回聲系統(tǒng)的輸出信號y(n)中恢復輸入信號x(n),寫出該系統(tǒng)的系統(tǒng)函數(shù)W(z)和差分方程。

(4)求解回聲消除系統(tǒng)的單位樣值響應,將它在0≤n≤4000樣本的近似值存入向量her中,并畫出該響應的波形。

(5)利用filter實現(xiàn)回聲消除,把消除回聲的語音數(shù)據(jù)存入向量we,利用plot畫出輸出信號。同時,利用sound聽聽輸入信號和輸出信號的聲音,并把消除回聲的語音保存在文件thegood.mat中。

圖2-7回聲系統(tǒng)的單位樣值響應

(6)將回聲系統(tǒng)和回聲消除系統(tǒng)級聯(lián)后求其單位樣值響應,將結果存入向量they中,并畫出總的單位樣值響應。

解:

(1)回聲系統(tǒng)的單位樣值響應可以用dimpulse函數(shù)求得,MATLAB程序如下：

num＝[1,zeros(1,1001－2),0.5］；

den＝[1,zeros(1,1001－1)］；

he＝dimpulse(num,den,4001)；

plot(he,′r′),title(′he′),xlabel(′n′)

axis([－500450001.2］)

回聲系統(tǒng)的單位樣值響應如圖2-7所示。

(2)把回聲系統(tǒng)方程y(n)＝x(n)＋αx(n－N)兩邊做Z變換,得系統(tǒng)函數(shù)為

(3)回聲消除系統(tǒng)的輸入信號是回聲系統(tǒng)的輸出y(n),設輸出信號z(n),消除系統(tǒng)是回聲系統(tǒng)的逆系統(tǒng),它的系統(tǒng)函數(shù)W(z)

為

回聲消除系統(tǒng)的差分方程是：

(4)這里用filter函數(shù)求回聲消除系統(tǒng)的單位樣值響應,MATLAB程序如下：

d＝[1zeros(1,4000)］；

her＝filter(1,num,d)；

plot(her,′r′),title(′her′),xlabel(′n′)

axis([－5004500－0.61.2］)

讀者可以通過(1)和(4)體會dimpulse函數(shù)和filter函數(shù)在求系統(tǒng)單位樣值響應時的不同用法。

回聲消除系統(tǒng)的單位樣值響應如圖2-8所示。

圖2-8回聲消除系統(tǒng)的單位樣值響應