matlab在DSP中的應用（四）---離散系統的衝激響應和階躍響應
一、實驗目的
(1)加深對離散線性移不變(LSI)系統基本理論的理解，明確差分方程與系統函式之間的關係。
(2)初步瞭解用MATLAB語言進行離散時間系統研究的基本方法。
(3)掌握求解離散時間系統衝激響應和階躍響應程式的編寫方法，瞭解常用子函式。
二、實驗涉及的MATLAB子函式
1.impz
功能：求解數字系統的衝激響應。
呼叫格式：
［h，t］＝impz(b，a)；求解數字系統的衝激響應h，取樣點數為預設值。
［h，t］＝impz(b，a，n)；求解數字系統的衝激響應h，取樣點數由n確定。
impz(b，a)；在當前視窗用stem(t，h)函數出圖。
2.dstep
功能：求解數字系統的階躍響應。
呼叫格式：
［h，t］＝dstep(b，a)；求解數字系統的階躍響應h，取樣點數為預設值。
［h，t］＝dstep(b，a，n)；求解數字系統的階躍響應h，取樣點數由n確定。
dstep(b，a)；在當前視窗用stairs(t，h)函數出圖。
3.filter
功能：對數字系統的輸入訊號進行濾波處理。
呼叫格式：
y＝filter(b，a，x)；對於由向量a、b定義的數字系統，當輸入訊號為x時，對x中的資料進行濾波，結果放於y中，長度取max(na，nb)。
［y，zf］＝filter(b，a，x)；除得到結果向量y外，還得到x的最終狀態向量zf。
y＝filter(b，a，x，zi)；可在zi中指定x的初始狀態
4.filtic
功能：為filter函式選擇初始條件。
呼叫格式：
z＝filtic(b，a，y，x)；求給定輸入x和y時的初始狀態。
z＝filtic(b，a，y)；求x＝0，給定輸入y時的初始狀態
其中，向量x和y分別表示過去的輸入和輸出：
x＝［x(－1)，x(－2)，…，x(－N)
y＝［y(－1)，y(－2)，…，y(－N)］
說明：以上子函式中的b和a，分別表示系統函式H(z)中由對應的分子項和分母項係數所構成的陣列。
如下面式(5-2)所示，H(z)按z－1(或z)的降冪排列。在列寫b和a係數向量時，兩個係數的長度必須相等
在MATLAB訊號處理工具箱中，許多用於多項式處理的函式，都採用以上的方法來處理分子項和分母項係數所構成的陣列。在後面的實驗中不再說明。
三、實驗原理
1.離散LSI系統的響應與激勵
由離散時間系統的時域和頻域分析方法可知，一個線性移不變離散系統可以用線性常係數差分方程表示：
也可以用系統函式來表示：
系統函式H(z)反映了系統響應與激勵間的關係。一旦上式中的bm和ak的資料確定了，則系統的性質也就確定了。其中特別注意：a0必須進行歸一化處理，即a0＝1。
對於複雜訊號激勵下的線性系統，可以將激勵訊號在時域中分解為單位脈衝序列或單位階躍序列，把這些單元激勵訊號分別加於系統求其響應，然後把這些響應疊加，即可得到複雜訊號加於系統的零狀態響應。
因此，求解系統的衝激響應和階躍響應尤為重要。由下圖可以看出一個離散LSI系統響應與激勵的關係。
2.用impz和dstep子函式求解離散系統的單位衝激響應和階躍響應
在MATLAB語言中，求解系統單位衝激響應和階躍響應的最簡單的方法是使用MATLAB提供的impz和dstep子函式。
下面舉例說明使用impz和dstep子函式求解系統單位衝激響應和階躍響應的方法。
例5-1 已知一個因果系統的差分方程為
6y(n)＋2y(n－2)＝x(n)＋3x(n－1)＋3x(n－2)＋x(n－3)
滿足初始條件y(－1)＝0，x(－1)＝0，求系統的單位衝激響應和階躍響應。
解：首先，將y(n)項的係數a0進行歸一化，得到：
y(n)＋1/3* y(n－2)＝1/6* x(n)＋1/2* x(n－1)＋1/2* x(n－2)＋1/6* x(n－3)
則：
編寫MATLAB程式如下(取N＝32點作圖)：
a=[1 0 1/3 0];
b=[1/6 1/2 1/2 1/6];
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,n);%時域單位衝激函式
gn=dstep(b,a,n+1);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
需要注意的是：求時域單位階躍函式時，需要使用的是gn=dstep(b,a,n+1)，而不是gn=dstep(b,a,n)，在這裡前者產生32* 1的矩陣（和n一樣），後者產生31* 1矩陣。
但是，如果我就是想用gn=dstep(b,a,n)呢？也可以，只要把stem(n,gn);改成stem(n(1:length(gn)),gn);就可以了，程式碼如下：
a=[1 0 1/3 0];
b=[1/6 1/2 1/2 1/6];
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,n);%時域單位衝激函式
gn=dstep(b,a,n);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n(1:length(gn)),gn);
title('系統單位階躍函式')
axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出：
注意，如果把stem(n,gn);改成stem(gn)的話，輸出的階躍函式會向右偏移：
實際上，上面的做法都不是很標準，很容易混淆。impz和dstep這兩個函式的第三個引數是用來指定取樣點的，而上面的程式碼中傳入的是一個數組，這也是說我們自己指定了取樣點的範圍了。但是在這裡還是選擇指定取樣點N比較好理解，而且不容易出錯。
另外注意到，axis這個函式不能隨便用的，對於axis([0 N 1.1*min(hn) 1.1*max(hn)])，，如果9.5<=hn<=10，那麼1.1*min(hn)就超過了hn的最大值了，這樣就看不到影象了！！！
axis這個函式是人為用來調整的，要用之前最好先繪出影象，看了影象的數值範圍後再使用。
下面給出程式碼：
a=[1 0 1/3 0];
b=[1/6 1/2 1/2 1/6];
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,N);%時域單位衝激函式
gn=dstep(b,a,N);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 0 1.1*max(gn)])輸出結果：
例 5-2 已知一個系統函式公式
解：分析上式可知，這是一個6階系統，直接用MATLAB語言列出其bm和ak係數：
a＝［1，0，0.34319，0，0.60439，0，0.20407］;
b＝［0.1321，0，－0.3963，0，0.3963，0，－0.1321］;
程式碼：
a=[1 0 0.34319 0 0.60439 0 0.20407];
b=[0.1321 0 -0.3963 0 0.3963 0 -0.1321];
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,N);%時域單位衝激函式
gn=dstep(b,a,N);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
3.用filtic和filter子函式求解離散系統的單位衝激響應
MATLAB提供了兩個子函式filtic和filter來求解離散系統的響應。
當輸入訊號為單位衝激訊號時，求得的響應即為系統的單位衝激響應；
當輸入訊號為單位階躍訊號時，求得的響應即為系統的單位階躍響應。
例5-3 已知一個因果系統的差分方程為 6y(n)－2y(n－4)＝x(n)－3x(n－2)＋3x(n－4)－x(n－6)，滿足初始條件y(－1)＝0，x(－1)＝0，求系統的單位衝激響應和單位階躍響應。時間軸上N取32點作圖。
解：將y(n)項的係數a0進行歸一化，得到： y(n)－1/3* y(n－4)＝1/6* x(n)－1/2* x(n－2)＋1/6* x(n－4)－1/6* x(n－6)
程式碼為：
a=[1 0 0 0 -1/3 0 0];
b=[1/6 0 -1/2 0 1/2 0 -1/6];
N=32;%作圖點數
n=0:N-1;
x1= n==0;%建立輸入單位衝激訊號x1(n)
x0=0;y0=0;%初始條件
xi=filtic(b,a,y0,x0);%求等效初始條件的輸入序列
hn=filter(b,a,x1,xi);%對輸入單位衝激訊號進行濾波，求單位衝激響應
x2= n>=0;%建立輸入單位階躍訊號x2(n)
gn=filter(b,a,x2,xi);%對輸入單位衝激訊號進行濾波，求單位衝激響應
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
四、實驗任務
(1)已知離散線性時不變系統的差分方程，請分別用impz和dstep子函式、filtic和filter子函式兩種方法求解系統的衝激響應和階躍響應。
①x(n)＋x(n－6)＝y(n)
②2y(n)－3y(n－1)＋y(n－2)＝x(n－1)
解：下面均以時間軸上N取32點作圖。
①：使用impz和dstep子函式：
程式碼：
a=[1 0 0 0 0 0 0];
b=[1 0 0 0 0 0 1];
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,N);%時域單位衝激函式
gn=dstep(b,a,N);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
使用filtic和filter子函式：
程式碼：
a=[1 0 0 0 0 0 0];
b=[1 0 0 0 0 0 1];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
x1= n==0;%建立輸入單位衝激訊號x1(n)
%xi=filtic(b,a,0);%求等效初始條件的輸入序列
hn=filter(b,a,x1);%對輸入單位衝激訊號進行濾波，求單位衝激響應
x2= n>=0;%建立輸入單位階躍訊號x2(n)
gn=filter(b,a,x2);%對輸入單位衝激訊號進行濾波，求單位衝激響應
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
②：使用impz和dstep子函式：
程式碼：
a=[2 -3 1];
b=[0 1 0];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,N);%時域單位衝激函式
gn=dstep(b,a,N);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
使用filtic和filter子函式：
程式碼：
a=[2 -3 1];
b=[0 1 0];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
x1= n==0;%建立輸入單位衝激訊號x1(n)
%xi=filtic(b,a,0);%求等效初始條件的輸入序列
hn=filter(b,a,x1);%對輸入單位衝激訊號進行濾波，求單位衝激響應
x2= n>=0;%建立輸入單位階躍訊號x2(n)
gn=filter(b,a,x2);%對輸入單位衝激訊號進行濾波，求單位衝激響應
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
(2)已知離散線性時不變系統的系統函式，請分別用impz和dstep子函式、filtic和filter子函式兩種方法求解系統的衝激響應和階躍響應。
其實給出系統函式和差分方程是沒有區別的。
解：下面均以時間軸上N取32點作圖。
①：使用impz和dstep子函式：
程式碼：
a=[1 -1 1];
b=[1 -0.5 0];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,N);%時域單位衝激函式
gn=dstep(b,a,N);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
使用filtic和filter子函式：
程式碼：
a=[1 -1 1];
b=[1 -0.5 0];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
x1= n==0;%建立輸入單位衝激訊號x1(n)
%xi=filtic(b,a,0);%求等效初始條件的輸入序列
hn=filter(b,a,x1);%對輸入單位衝激訊號進行濾波，求單位衝激響應
x2= n>=0;%建立輸入單位階躍訊號x2(n)
gn=filter(b,a,x2);%對輸入單位衝激訊號進行濾波，求單位衝激響應
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
②：使用impz和dstep子函式：
程式碼：
a=[1 0 0 0 0 0];
b=[1 0.5 -0.5 -1 -0.5 1];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
hn=impz(b,a,N);%時域單位衝激函式
gn=dstep(b,a,N);%時域單位階躍函式
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
使用filtic和filter子函式：
程式碼：
a=[1 0 0 0 0 0];
b=[1 0.5 -0.5 -1 -0.5 1];
tmp=a(1);
a=a./tmp;%歸一化
b=b./tmp;
N=32;%作圖點數
n=0:N-1;
x1= n==0;%建立輸入單位衝激訊號x1(n)
%xi=filtic(b,a,0);%求等效初始條件的輸入序列
hn=filter(b,a,x1);%對輸入單位衝激訊號進行濾波，求單位衝激響應
x2= n>=0;%建立輸入單位階躍訊號x2(n)
gn=filter(b,a,x2);%對輸入單位衝激訊號進行濾波，求單位衝激響應
subplot(1,2,1),stem(n,hn);
title('系統單位衝激函式')
%axis([0 N 1.1*min(hn) 1.1*max(hn)])
subplot(1,2,2),stem(n,gn);
title('系統單位階躍函式')
%axis([0 N 1.1*min(gn) 1.1*max(gn)])輸出結果：
思考題：
①離散LSI系統的差分方程和系統函式有何聯絡？公式中的bm和ak係數在編寫程式時須注意什麼問題？
對差分方程兩邊同時進行傅立葉變換後整合化簡後的結果就是系統函式。公式中的bm和ak係數在編寫程式時須注意補上缺項以及搞清楚兩者分別對應的是輸入訊號還是輸出訊號的係數。（bm對應x[n]、 ak對應y[n]）
②簡述用子函式filter求解離散系統的單位衝激響應和單位階躍響應的基本思路。
1.根據給定的差分方程和系統函式求出bm和ak
2.用filtic函式求出初始狀態求等效初始條件的輸入序列xi ，若沒有初始條件則不用
3.用filter函式結合xi、bm、ak、單位衝激函式得到單位衝激響應（若沒有初始條件則不用xi）
4.用filter函式結合xi、bm、ak、單位階躍函式得到單位階躍響應（若沒有初始條件則不用xi）