第3章 离散时间LTI系统的时域分析
3.1 实验目的
? 学会运用MATLAB求解离散时间系统的零状态响应； ? 学会运用MATLAB求解离散时间系统的单位取样响应； ? 学会运用MATLAB求解离散时间系统的卷积和。
3.2 实验原理及实例分析
3.2.1 离散时间系统的响应
离散时间LTI系统可用线性常系数差分方程来描述，即
?ay(n?i)??bx(n?j) （3-1）
iji?0j?0NM其中，ai（i?0，1，?，N）和bj（j?0，1，?，M）为实常数。
MATLAB中函数filter可对式（13-1）的差分方程在指定时间范围内的输入序列所产生的响应进行求解。函数filter的语句格式为
y=filter(b,a,x)
其中，x为输入的离散序列；y为输出的离散序列；y的长度与x的长度一样；b与a分别为差分方程右端与左端的系数向量。
【实例3-1】 已知某LTI系统的差分方程为
3y(n)?4y(n?1)?2y(n?2)?x(n)?2x(n?1)
试用MATLAB命令绘出当激励信号为x(n)?(1/2)u(n)时，该系统的零状态响应。
解：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('系统响应y(n)')
n程序运行结果如图3-1所示。
图3-1 实例3-1系统的零状态响应
3.2.2 离散时间系统的单位取样响应
系统的单位取样响应定义为系统在?(n)激励下系统的零状态响应，用h(n)表示。MATLAB求解单位取样响应可利用函数filter，并将激励设为前面所定义的impDT函数。例如，求解实例13-1中系统的单位取样响应时，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'),title('系统单位取样响应h(n)')
程序运行结果如图3-2所示。
图3-2 实例13-1的系统单位取样响应
MATLAB另一种求单位取样响应的方法是利用控制系统工具箱提供的函数impz来实现。impz函数的常用语句格式为
impz(b,a,N)
其中，参数N通常为正整数，代表计算单位取样响应的样值个数。
【实例3-2】 已知某LTI系统的差分方程为
3y(n)?4y(n?1)?2y(n?2)?x(n)?2x(n?1)
利用MATLAB的impz函数绘出该系统的单位取样响应。
解：MATLAB源程序为
>>a=[3 -4 2]; >>b=[1 2]; >>n=0:30;
>>impz(b,a,30),grid on
>>title('系统单位取样响应h(n)')
程序运行结果如图3-3所示，比较图3-2和图3-3，不难发现结果相同。
图3-3 系统单位取样响应
3.2.3 离散时间信号的卷积和运算
由于系统的零状态响应是激励与系统的单位取样响应的卷积，因此卷积运算在离散时间信号处理领域被广泛应用。离散时间信号的卷积定义为
y(n)?x(n)*h(n)?m????x(m)h(n?m) （3-2）
?可见，离散时间信号的卷积运算是求和运算，因而常称为“卷积和”。
MATLAB求离散时间信号卷积和的命令为conv，其语句格式为
y=conv(x,h)
其中，x与h表示离散时间信号值的向量；y为卷积结果，它默认序列从n＝0开始。但是如果序列是从一负值开始，即
?x(n):nx1?n?nx2??h(n):nh1?n?nh2?
如果nx1<0或nh1<0就不能直接采用conv函数。其卷积结果序列为{y(n):nx1?nh1?n?nx2?nh2}，这样就可以构成一个新的卷积函数conv_m。如下所示：
function[y,ny]=conv_m(x,nx,h,nh)
ny1=nx(1)+nh(1);ny2=nx(length(x))+nh(length(h)); ny=[ny1:ny2]; y=conv(x,h)
值得注意的是用MATLAB进行卷积和运算时，无法实现无限的累加，只能计算时限信号的卷积。
【实例3-3】 已知某系统的单位取样响应为h?n??0.8n?u?n??u?n?8??，试用MATLAB求当激励信号为x(n)?u(n)?u(n?4)时，系统的零状态响应。
解：MATLAB中可通过卷积求解零状态响应，即x(n)*h(n)。由题意可知，描述h(n)向量的长度至少为8，描述x(n)向量的长度至少为4，因此为了图形完整美观，我们将h(n)向量和x(n)向量加上一些附加的零值。MATLAB源程序为
nx=-1:5; %x(n)向量显示范围(添加了附加的零值) nh=-2:10; %h(n)向量显示范围(添加了附加的零值) x=uDT(nx)-uDT(nx-4);
h=0.8.^nh.*(uDT(nh)-uDT(nh-8)); [y,ny]=conv_m(x,nx,h,nh); subplot(311)
stem(nx,x,'fill'),grid on xlabel('n'),title('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)') axis([-4 16 0 3])
程序运行结果如图3-5所示。
图3-5 利用卷积和法求解系统的零状态响应
3.3 编程练习
1.
试用MATLAB命令求解以下离散时间系统的单位取样响应。 （1）3y(n)?4y(n?1)?y(n?2)?x(n)?x(n?1) （2） 2.
已知某系统的单位取样响应为h?n??()?u?n??u?n?10??，试用MATLAB求当激励
n5y(n)?6y(n?1)?10y(n?2)?x(n) 278信号为x(n)?u(n)?u(n?5)时，系统的零状态响应。
第4章 z变换及离散时间LTI系统的z域
分析
4.1 实验目的
? ? ? ?
学会运用MATLAB求离散时间信号的z变换和z反变换； 学会运用MATLAB分析离散时间系统的系统函数的零极点；
学会运用MATLAB分析系统函数的零极点分布与其时域特性的关系； 学会运用MATLAB进行离散时间系统的频率特性分析。
4.2 实验原理及实例分析
4.2.1 z正反变换
序列x?n?的z变换定义为
X?z??Z?x?n???n????x?n?z??n （4-1）
其中，符号Z表示取z变换，z是复变量。相应地，单边z变换定义为
X?z??Z?x?n????x?n?z?n （4-2）
n?0?MATLAB符号数学工具箱提供了计算离散时间信号单边z变换的函数ztrans和z反变换函数iztrans，其语句格式分别为
Z=ztrans(x) x=iztrans(z)
上式中的x和Z分别为时域表达式和z域表达式的符号表示，可通过sym函数来定义。
【实例4-1】 试用ztrans函数求下列函数的z变换。
（1）x(n)?acos(?n)u(n)； （2）x(n)?[2解：（1）z变换MATLAB源程序为
>>x=sym('a^n*cos(pi*n)'); >>Z=ztrans(x);
>>simplify(Z) %对Z进行简化运算 ans=
z/(z+a)
nn?1?(?2)n?1]u(n)。
（2）z变换MATLAB源程序为
>>x=sym('2^(n-1)-(-2)^(n-1)');
>>Z=ztrans(x); >>simplify(Z) ans=
z^2/(z-2)/(z+2)
【实例4-2】 试用iztrans函数求下列函数的z反变换。
8z?19z(2z2?11z?12)（1）X(z)?2 （2）X(z)?
z?5z?6(z?1)(z?2)3解：（1）z反变换MATLAB源程序为
>>Z=sym('(8*z-19)/(z^2-5*z+6)'); >>x=iztrans(Z); >>simplify(x) ans=
-19/6*charfcn[0](n)+5*3^(n-1)+3*2^(n-1)
其中，charfcn[0](n)是?(n)函数在MATLAB符号工具箱中的表示，反变换后的函数形式为
x(n)??19?(n)?(5?3n?1?3?2n?1)u(n)。 6（2）z反变换MATLAB源程序为
>>Z=sym('z*(2*z^2-11*z+12)/(z-1)/(z-2)^3'); >>x=iztrans(Z); >>simplify(x) ans=
-3+3*2^n-1/4*2^n*n-1/4*2^n*n^2
其函数形式为x(n)?(?3?3?2?n1n12nn2?n2)u(n)。 44
如果信号的z域表示式X(z)是有理函数，进行z反变换的另一个方法是对X(z)进行部分分式展开，然后求各简单分式的z反变换。设X(z)的有理分式表示为
b0?b1z?1?b2z?2???bmz?mB(z) （4-3） X(z)???1?2?nA(z)1?a1z?a2z???anzMATLAB信号处理工具箱提供了一个对X(z)进行部分分式展开的函数residuez，其语句格式为
[R,P,K]=residuez(B,A)
其中，B，A分别表示X(z)的分子与分母多项式的系数向量；R为部分分式的系数向量；P为极点向量；K为多项式的系数。若X(z)为有理真分式，则K为零。
Y(z)?从
B(z)r(1)r(2)r(N)?1???...??k(1)?k(2)z?...?1?1?1 A(z)1?p(1)z1?p(2)z1?p(N)z得
出
其
时
域
信
号
为
：
而
y(n)?r(1)[p(1)]nu(n)?r(2)[p(2)]nu(n)?...?r(N)[p(N)]nu(n)?k(1)?(n)?k(2)?(n?1)?...
【实例4-3】 试用MATLAB命令对函数X(z)?展开，并求出其z反变换。
解：MATLAB源程序为
>>B=[18]; >>A=[18,3,-4,-1]; >>[R,P,K]=residuez(B,A) R= 0.3600 0.2400 0.4000 P= 0.5000 -0.3333 -0.3333 K= []
18进行部分分式
18?3z?1?4z?2?z?3从运行结果可知，p2?p3，表示系统有一个二重极点。所以，X(z)的部分分式展开为
X(z)?因此，其z反变换为
0.36?0.240.4 ?1?0.5z?11?0.3333z?1（1?0.3333z?1)2x(n)?[0.36?(0.5)n?0.24?(?0.3333)n?0.4(n?1)(?0.3333)n]u(n)
4.2.2 系统函数的零极点分析
离散时间系统的系统函数定义为系统零状态响应的z变换与激励的z变换之比，即
H(z)?如果系统函数H(z)的有理函数表示式为
Y(z) （4-4） X(z)b1zm?b2zm?1???bmz?bm?1 （4-5） H(z)?a1zn?a2zn?1???anz?an?1那么，在MATLAB中系统函数的零极点就可通过函数roots得到，也可借助函数tf2zp得到。 1， roots的格式语句为：p=roots(A),其中A为待求根的多项式的系数构成的行向量，返回向量p则包含该多项式所有的根位置列向量。 2，tf2zp的语句格式为
[Z,P,K]=tf2zp(B,A)
其中，B与A分别表示H(z)的分子与分母多项式的系数向量。它的作用是将H(z)的有理分式表示式转换为零极点增益形式，即
H(z)?k(z?z1)(z?z2)?(z?zm) （4-6）
(z?p1)(z?p2)?(z?pn)z?0.32 2z?z?0.16【实例4-4】 已知一离散因果LTI系统的系统函数为
H(z)?试用MATLAB命令求该系统的零极点。
解：1，用roots函数求系统的零极点： a=[1 0.32]; b=[1 1 0.16]; r=roots(a) r =
-0.3200 p=roots(b) p =
-0.8000 -0.2000
2，用tf2zp函数求系统的零极点，MATLAB源程序为
>>B=[1,0.32]; >>A=[1,1,0.16]; >>[R,P,K]=tf2zp(B,A) R= -0.3200 P= -0.8000 -0.2000 K= 1
因此，零点为z??0.32，极点为p1??0.8与p2??0.2。
若要获得系统函数H(z)的零极点分布图，可直接应用zplane函数，其语句格式为
zplane(B,A)
其中，B与A分别表示H(z)的分子和分母多项式的系数向量。值注意的是：求系统函数零极点时，离散系统的系统函数可能有两种形式，一种是分子分母多项式按z的降幂次序排列，另一种是按z的升幂次序排列。若是以z的降幂次序排列，则系数向量一定要由多项
?1式的最高幂次开始，一直到常数项，缺项要用0补齐；若以z的升幂次序排列，则分子分母多项式系数向量的维数一定要相同，不足的要用0补齐，否则零点或极点有可能被漏掉。
【实例4-5】 已知一离散因果LTI系统的系统函数为
?1z2?0.361?0.36z?2H(z)?2?z?1.52z?0.681?1.52z?1?0.68z?2
试用MATLAB命令绘出该系统的零极点分布图。
解：用zplane函数求系统的零极点，MATLAB源程序为
>>B=[1,0,-0.36]; >>A=[1,-1.52,0.68]; >>zplane(B,A),grid on >>legend('零点','极点') >>title('零极点分布图')
程序运行结果如图14-1所示。可见，该因果系统的极点全部在单位圆内，故系统是稳定的。
图4-1 零极点分布图
4.2.3 系统函数的零极点分布与其时域特性的关系
与拉氏变换在连续系统中的作用类似，在离散系统中，z变换建立了时域函数h(n)与z域函数H(z)之间的对应关系。因此，z变换的函数H(z)从形式可以反映h(n)的部分内在性质。我们仍旧通过讨论H(z)的一阶极点情况，来说明系统函数的零极点分布与系统时域特性的关系。
【实例14-6】 试用MATLAB命令画出现下列系统函数的零极点分布图、以及对应的时域单位取样响应h(n)的波形，并分析系统函数的极点对时域波形的影响。
（1）H1(z)?zzz （2）H2(z)? （3）H3(z)?2
z?0.8z?0.8z?1.2z?0.72
zzz （5）H5(z)?2 （6）H6(s)? z?1z?1.2z?1.6z?1z（7）H7(z)?2
z?2z?1.36（4）H4(z)?解：MATLAB源程序为
>>b1=[1,0]; >>a1=[1,-0.8]; >>subplot(121) >>zplane(b1,a1)
>>title('极点在单位圆内的正实数') >>subplot(122)
>>impz(b1,a1,30);grid on; >>figure >>b2=[1,0]; >>a2=[1,0.8]; >>subplot(121) >>zplane(b2,a2)
>>title('极点在单位圆内的负实数') >>subplot(122)
>>impz(b2,a2,30);grid on; >>figure >>b3=[0,1,0]; >>a3=[1,-1.2,0.72]; >>subplot(121) >>zplane(b3,a3)
>>title('极点在单位圆内的共轭复数') >>subplot(122)
>>impz(b3,a3,30);grid on; >>figure >>b4=[1,0]; >>a4=[1,-1]; >>subplot(121) >>zplane(b4,a4)
>>title('极点在单位圆上为实数1') >>subplot(122) >>impz(b4,a4);grid on; >>figure >>b5=[0,1,0]; >>a5=[1,-1.6,1]; >>subplot(121) >>zplane(b5,a5)
>>title('极点在单位圆上的共轭复数') >>subplot(122)
>>impz(b5,a5,30);grid on;
>>figure >>b6=[1,0]; >>a6=[1,-1.2]; >>subplot(121) >>zplane(b6,a6)
>>title('极点在单位圆外的正实数') >>subplot(122)
>>impz(b6,a6,30);grid on; >>figure >>b7=[0,1,0]; >>a7=[1,-2,1.36]; >>subplot(121) >>zplane(b7,a7)
>>title('极点在单位圆外的共轭复数') >>subplot(122)
>>impz(b7,a7,30);grid on;
程序运行结果分别如图14-2的（a）、（b）、（c）、（d）、（e）、（f）、（g）所示。
(a)
(b)
(c)
(d)
(e)
(f)
(g)
图4-2 系统函数的零极点分布与其时域特性的关系
从图14-2可知，当极点位于单位圆内时，h(n)为衰减序列；当极点位于单位圆上时，
h(n)为等幅序列；当极点位于单位圆外时，h(n)为增幅序列。若h(n)有一阶实数极点，
则h(n)为指数序列；若h(n)有一阶共轭极点，则h(n)为指数振荡序列；若h(n)的极点位于虚轴左边，则h(n)序列按一正一负的规律交替变化。
4.2.4 离散时间LTI系统的频率特性分析
对于因果稳定的离散时间系统，如果激励序列为正弦序列x(n)?Asin(n?)u(n)，则系
统的稳态响应为yss(n)?A|H(ej?)|sin[n???(?)]u(n)。其中，H(ej?)通常是复数。离散时间系统的频率响应定义为
H(ej?)?|H(ej?)|ej?(?) （14-7）
其中，|H(ej?)|称为离散时间系统的幅频特性；?(?)称为离散时间系统的相频特性；
H(ej?)是以?s（?s?2?，若零T?1，?s?2?）为周期的周期函数。因此，只要分T析H(ej?)在|?|??范围内的情况，便可分析出系统的整个频率特性。
MATLAB提供了求离散时间系统频响特性的函数freqz，调用freqz的格式主要有两种。一种形式为
[H,w]=freqz(B,A,N) 其中，B与A分别表示H(z)的分子和分母多项式的系数向量，值得注意的是，B，A采用z?1的升次幂形式；N为正整数，默认值为512；返回值w包含[0,?]范围内的N个频率等分点；返回值H则是离散时间系统频率响应H(ej?)在0~?范围内N个频率处的值。另一种形式为
[H,w]=freqz(B,A,N,’whole’)
与第一种方式不同之处在于角频率的范围由[0,?]扩展到[0,2?]。
z2?0.96z?0.9028【实例4-6】 用MATLAB命令绘制系统H(z)?2的频率响应曲线。
z?1.56z?0.8109 解：利用函数freqz计算出H(ej?)，然后利用函数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)'),ylabel('Magnitude') >>title('离散系统幅频特性曲线') >>subplot(212) >>plot(w,Hp),grid on
>>xlabel('\\omega(rad/s)'),ylabel('Phase') >>title('离散系统相频特性曲线')
程序运行结果如图4-2所示。
图4-3 离散系统频响特性曲线
4.3 编程练习
1.
2z4?16z3?44z2?56z?32试用MATLAB的residuez函数，求出X(z)?的部分分
3z4?3z3?15z2?18z?12式展开和。
试用MATLAB画出下列因果系统的系统函数零极点分布图，并判断系统的稳定性。
2.
2z2?1.6z?0.9（1）H(z)?3 2z?2.5z?1.96z?0.48（2）H(z)?z?1 432z?0.9z?0.65z?0.873z3.
z2试用MATLAB绘制系统H(z)?的频率响应曲线。
31z2?z?48
图4-3 离散系统频响特性曲线
4.3 编程练习
1.
2z4?16z3?44z2?56z?32试用MATLAB的residuez函数，求出X(z)?的部分分
3z4?3z3?15z2?18z?12式展开和。
试用MATLAB画出下列因果系统的系统函数零极点分布图，并判断系统的稳定性。
2.
2z2?1.6z?0.9（1）H(z)?3 2z?2.5z?1.96z?0.48（2）H(z)?z?1 432z?0.9z?0.65z?0.873z3.
z2试用MATLAB绘制系统H(z)?的频率响应曲线。
31z2?z?48