离散数字时间系统分析

1. 掌握离散LSI系统的单位序列响应、单位阶跃响应和任意激励下响应的MATLAB求解方法。

2. 掌握离散LSI系统的频域分析方法； 3. 掌握离散LSI系统的复频域分析方法；

4. 掌握离散LSI系统的零极点分布与系统特性的关系。 二、实验原理及方法 离散LSI系统的时域分析

描述一个N阶线性时不变离散时间系统的数学模型是线性常系统差分方程，N阶LSI离散系统的差分方程一般形式为

k0aky(nk)bix(ni) （6.1）

i0NM也可用系统函数来表示

Y(z)H(z)X(z)bziMiazkk0i0Nkb(z)b0b1z1b2z2bMzM （6.2） a(z)1a1z1a2z2aNzN系统函数H(z)反映了系统响应和激励间的关系。一旦上式中ak，bi的数据确定了，系统的性质也就确定了。特别注意a0必须进行归一化处理，即a01。

对于复杂信号激励下的线性系统，可以将激励信号在时域中分解为单位序列或单位阶跃序列的线性叠加，把这些单元激励信号分别加于系统求其响应，然后把这些响应叠加，即可得到复杂信号作用于系统的

零状态响应。因此，求解系统的单位序列响应和单位阶跃响应尤为重要。由图6-1可以看出一个离散LSI系统响应与激励的关系。

x(n)h(n)y(n)x(n)*h(n)

X(z)H(z)Y(z)X(z)H(z)图6-1 离散LSI系统响应与激励的关系 (1) 单位序列响应（单位响应）

单位响应h(n)是指离散LSI系统在单位序列(n)激励下的零状态响应，因此h(n)满足线性常系数差分方程(6.1)及零初始状态，即

ah(nk)b(ni)， h(1)h(2)0 (6.3)

kik0i0NM按照定义，它也可表示为

h(n)h(n)(n) (6.4) 对于离散LSI系统，若其输入信号为x(n)，单位响应为h(n)，则其零状态响应yzs(n)为

yzs(n)x(n)*h(n) (6.5) 可见，h(n)能够刻画和表征系统的固有特性，与何种激励无关。一旦知道了系统的单位响应h(n)，就可求得系统对任何输入信号x(n)所产生的零状态响应yzs(n)。

MATLAB提供了专门用于求离散时间系统冲激响应的函数impz()，其调用格式有 [h,n]=impz(b,a)

b求解离散系统的单位脉冲响应，其中b[b0,b1,2,2

M,b，]a[1,a1,a2,,aN]，n[0,1,2,]；

[h,n]=impz(b,a,N)

求解离散系统的单位响应，采样点数由N确定，n[0,1,2,,N-1]； impz(b,a) ：在当前绘图窗口，用stem(n,h)绘出图形。 （2）单位阶跃响应

单位阶跃响应s(n)是指离散LTI系统在单位阶跃序列u(n)激励下的零状态响应，它可以表示为

s(n)u(n)h(n)mh(m) (6.6)

n上式表明，离散LSI系统的单位阶跃响应是单位脉冲响应的累加和，系统的单位阶跃响应和系统的单位响应之间有着确定的关系，因此，单位阶跃响应也能完全刻画和表征一个LSI系统。

MATLAB提供了专门用于求离散系统单位阶跃响应的函数stepz( )，其调用格式有

[s,n]=stepz(b,a) ：求解离散系统的单位阶跃响应，其中

b[bb0,b1,2,M,b][1,a1,a2,,aN]，n[0,1,2,]； ，a[s,n]=stepz(b,a,N) ：求解离散系统的单位阶跃响应，采样点数由N确定，n[0,1,2,,N-1]；

stepz(b,a) ：在当前窗口，用stem(n,s)绘出图形。 任意激励下的零状态响应

已经知道，离散LSI系统可用常系数线性差分方程(6.1)式来描述，Matlab提供的函数dlsim( )能对上述差分方程描述的离散LSI系统的响应进行仿真，该函数不仅能绘制指定时间范围内的系统响应波形

3

图，而且还能求出系统响应的数值解。其调用格式有 dlsim(b,a, x) ：求解输入序列为x的零状态响应

Matlab总是把由分子和分母多项式表示的任何系需要特别强调的是，

统都当作是因果系统。所以，利用impz (b,a)，stepz(b,a)，dlsim(b,a,x)函数求得的响应总是因果信号。

同时，卷积和也是LSI系统求解零状态响应的重要工具之一。假设系统的输入信号为x(n)，单位响应为h(n)，则系统的零状态响应yzs(n)可由（6.5）式求解。Matlab提供了专门用于求离散系统卷积和的函数conv( )，其调用格式有

y=conv(x,h) ：求解序列x，h的卷积和，若序列x的长度为n1，序列h的长度为n2，卷积和y的长度为n1+n2-1。这一点需要特别注意，否则，作图时容易造成横纵坐标长度不匹配。

例如已知输入序列x(n)0.8n(0n19)，系统的单位响应

h(n)u(n)(0n9)，求系统的零状态响应。参考程序如下：

n1=0:19; x=0.8.^n1; subplot(231) stem(n1,x,'.'); title('x(n)'); n2=0:9; h=ones(1,10); subplot(232)

4

stem(n2,h,'.'); title('h(n)');

y=conv(x,h); %求卷积和

m1=n1(1)+n2(1); %卷积和的起始位置

m2=n1(length(n1))+n2(length(n2)); %卷积和的结束位置 n3=m1:m2; % 卷积和的长度，这是需要大家特别注意的地方。 subplot(233) stem(n3,y,'.'); title('y(n)');

运行结果如图6-2所示。

10.80.60.40.2001020x(n)h(n)10.80.60.40.20051054321002040y(n)

图6-2 输入序列、单位序列及零状态响应 （4）带初始状态的任意激励下的全响应

任意激励下的离散LSI系统的全响应为零输入响应和零状态响应之和，表示为

y(n)yzi(n)yzs(n) （6.7）

在理论学习的过程中，同学们对低阶差分方程的求解已颇为头痛，高Matlab提供了用于求离散系统全响阶差分方程直接求解几乎不可能。

5

应的函数filter( )，其调用格式有 y=filter( b,a,x) ：求解零状态响应；

y=filter( b,a,x,zi) ：求解初始条件为zi的系统的全响应，zi向量的长度为max(length(a),length(b))-1，返回值为系统的全响应。

z = filtic(b,a,y,x)：将初始状态转换为初始条件，其中

x[x(1)x,(2x),(3x)，,my[y(1),y(2),y(3),,y(n)]；

z = filtic(b,a,y)：将初始状态转换为初始条件，其中x0，

y[y(1),y(2),y(3),,y(n)]。

例6-1 某LSI系统的差分方程表示式为

y(n)3y(n1)2y(n2)x(n)

满足初始状态y(1)0，y(2)0.5，用filter子函数求系统输入为

x(n)2nu(n)时的零输入、零状态及全响应。

1n2ny(n)[2(2)(1)n]u(n)，下通过解差分方程，可以得到全响应为

33面使用filter子函数对系统差分方程进行求解，同时将求解结果与理论推导出来的结果进行比较，求解结果与比较见图6-3。参考程序如下：

a=[1 3 2]; %输入系统系数a、b b=[1]; N=20; n=0:N-1;

x0=zeros(1,N); %建立零输入信号 y0=[0 0.5]; %输入初始状态

6

x=2.^n; %建立输入信号

zi=filtic(b,a,y0); %有初始状态计算初始条件 yzi=filter(b,a,x0,zi); %求零输入响应 zi1=filtic(b,a,0); %计算初始状态为0的情况 yzs=filter(b,a,x,zi1); %求零状态响应 y=filter(b,a,x,zi); %求系统的全响应

yy=2/3.*(-1).^n+1/3.*2.^n-(-2).^n; %输入系统全响应的理论值subplot(231),stem(n,x,'.'); title(‘输入信号’);

subplot(232),stem(n,yzi,'.');title(‘系统的零输入响应’); subplot(233),stem(n,yzs,'.');title('系统的零状态响应'); subplot(223),stem(n,y,'.');title('系统的全响应');

subplot(224),stem(n,yy,'.');title('理论计算的系统全响应');

7

65x 10输入信号1.516x 10系统的零输入响应5x 10系统的零状态响应42040.50-0.5-2-42005x 101020-10105x 10102001020系统的全响应理论计算的系统全响应

105500-505101520-505101520图6-3 系统的零输入响应、零状态响应及全响应

由图6-3可看出，理论计算的全响应和用filter（）函数计算出来的全响应是一致的，从而为我们求解差分方程提供了一种简便的方法。 2 离散LSI系统的复频域（Z域）分析 （1）利用Z变换解差分方程

在前面图6-1中表示了离散系统的响应与激励的关系，由图可知，系统的响应既可以用时域的方法求解，也可以用Z域的方法求解。当已知系统输入序列的Z变换X(z)，系统函数H(z)时，系统响应序列的Z变换可由Y(z)X(z)H(z)求出。Matlab提供了用于求序列Z变换和Z反变换的函数，其调用格式有

X=ztrans(x)：求无限长序列x的Z变换，返回Z变换的表达式，注意这里x，X都是符号表达式；

8

x=iztrans(X)：求X（z）的Z反变换x(n)，返回Z反变换的表达式，注意这里x，X都是符号表达式；

[r,p,c]=residuez(b,a)：把b(z)/a(z)展开成部分分式；

[b,a]=residuez(r,p,c):根据部分分式的r、p、c数组，返回有理多项式。例6-2 某离散系统，其输入为x(n)h(n)an，求当a2时系统的零状态响应。分别利用Z变换和卷积和进行求解，结果如图6-4。参考程序如下：

syms a n z %定义符号变量 x=a.^n; %输入x X=ztrans(x); %对x作Z变换 H=X;

Y=X*H; %利用Z变换求零状态响应的Z域解 y=iztrans(Y) %求Z反变换 返回值为 a^n+a^n*n 再编写画图程序： n=0:20; y=(1+n).*2.^n; stem(n,y,’.’);

利用conv()实现的过程： n=0:20; x=2.^n;

9

y=conv(x,x); stem(y(1:21),'.')

2.521.510.5002468101214161820利用Z变换得到的零状态响应x 107x 1072.521.510.500255152010利用卷积和得到的零状态响应

图6-4 利用Z变换求解零状态响应

由图6-4可知，用Z变换和卷积和的方法得到的系统零状态响应是一致的。在实际生活中，一般采用卷积和的方法，因为这种方法采用的数值运算，编程简单。

（2）系统的零极点分布与系统因果性和稳定性的关系

因果系统的单位响应h(n)一定满足当n0时，h(n)=0，那么其系统函数H(z)的收敛域一定包含点，即点不是极点，极点分布在某个圆的圆内，收敛域是圆外区域。 系统稳定要求

n|h(n)| ，对照z变换定义，系统稳定要求收敛域

包含单位圆。如果系统因果且稳定，收敛域包含点和单位圆，那么收敛域可表示为：

r|z|,0r1 (6.8)

MATLAB提供了用于求系统零极点的函数，其调用格式有 roots()：利用多项式求根函数来确定系统函数的零极点位置；

10

roots(a)：求极点位置，a为系统函数H(z)分母多项式所构成的系数向量；

roots(b)：求零点位置，b为系统函数H(z)分子多项式所构成的系数向量；

zplane(b,a)：绘制由行向量b和a构成的系统函数的零极点分布图； zplane(z,p)：绘制由列向量z确定的零点、列向量p确定的极点构成的零极点分布图。

例6-3 已知离散系统的系统函数为

41.6z11.6z24z3H(z)

10.4z10.35z20.4z3求该系统的零极点，并绘出零极点分布图，通过零极点的分布判断系统的因果性和稳定性。 Matlab参考程序如下： b=[4 -1.6 -1.6 4]; a=[1 0.4 0.35 -0.4];

z=roots(b) %求系统的零点 p=roots(a) %求系统的极点 subplot(221),zplane(b,a); title('系统的零极点分布图'); subplot(222),impz(b,a,20); title('系统的单位响应');

运行结果如下，其零极点分布图和单位响应见图6-5所示。 z =

11

-1.0000 0.7000 + 0.7141i 0.7000 - 0.7141i p =

-0.4500 + 0.7730i -0.4500 - 0.7730i 0.5000

系统的零极点分布图110系统的单位响应Imaginary Part Amplitude0.50-0.5-1-1-0.500.5Real Part15图

0-50510n (samples)156-5 离散系统的零极点分布及单位响应

由图6-5可见，该系统的所有极点均在单位圆内，系统的单位响应随着n增大而收敛，因此，该系统是一个因果稳定系统。 （3）系统的零极点分布与系统频率响应的关系 将式(6.2)因式分解，得到

1H(z)AmN (6.9)

1(1dkz)k1(1cmz1)M式中，Ab0a0，cm是H(z)的零点，dk是其极点。A参数影响频率响应的幅度大小，影响系统特性的是零点cm和极点dk的分布。下面采用

12

几何方法研究系统零极点分布对系统频率特性的影响。 将式(6.9)的分子、分母同乘以zNM，得到：

1H(z)AmNk1(1cMmz)AzNMz1)1(zcm1Nk1Mm)(1dk(zd (6.10)

k)假设系统稳定，将zej代入上式，得到频率响应

j(ecm)jMH(e)Aejj(NM)m1Nk1H(e)ejjarg[H(ej)](ed (6.11)

k)设NM，由式(6.11)得到

1H(ej)AmNk1(ecjjNm)(ed (6.12)

k)在z平面上，ejcm用一根由零点cm指向单位圆（ej）上任一点B的向量cmB表示，同样ejdk用一根由极点dk指向单位圆（ej）上任一点B的向量dkB来表示，cmB和dkB分别称为零点矢量和极点矢量，用极坐标表示为： cmBcmej，dkBdkej。

mk将cmB和dkB表示式代入式(6.12)，得到

1H(ej)AmNk1(eNjcm)k(edj1AmNck1NmBH(ej)ej() B)dk|H(e)|Ajcm1Nk1Nmd (6.13)

k 13

()mk (6.14)

m1k1NN系统或者信号的频率特性由式(6.13)和式(6.14)确定。按照式(6.13)，知道零极点的分布后，可以很容易地确定零极点位置对系统特性的影响。当B点转到极点附近时，极点矢量长度最短，因而幅度特性可能出现峰值，且极点愈靠近单位圆，极点矢量长度愈短，峰值愈高愈尖锐。如果极点在单位圆上，则幅度特性为，系统不稳定。对于零点，情况相反，当B点转到零点附近时，零点矢量长度变短，幅度特性将出现谷值，且零点愈靠近单位圆，谷值愈接近零。当零点在单位圆上时，谷值为零。综上所述，极点位置主要影响频响的峰值位置及尖锐程度，零点位置主要影响频响的谷值位置及形状。

Matlab提供了专门用于求离散系统频率响应的函数freqz(),其调用格式如下：

[H,w] = freqz(b,a,n)：返回数字系统的n点频率值（复数），这n个点均匀地分布在[0，π]上，系统默认的采样点数目为512点；

[H,f] = freqz(b,a,n,Fs)：用于对H(ej)在[0，Fs/2]上等间隔采样n点，采样点频率及相应的频响值分别存放在f和H中。

H = freqz(b,a,w)：用于对H(ej)在[0，2π]上进行采样，采样频率点由w指定。

H = freqz(b,a,f,Fs)：用于对H(ej)在[0，Fs]上进行采样，采样频率点由f指定。

freqz(b,a,n)：用于在当前图形窗口中绘制幅频和相频特性曲线。 下面介绍几个实用的函数：

14

mag=abs(H):求解系统的绝对幅频响应；

db=20*log10((mag+eps)/max(mag))：求解系统的相对幅频响应； pha=angle(H): 求解系统的相频响应； grd=grpdelay(b,a,w): 求解系统的群延迟。

下面举例说明如何利用上述函数计算并绘制系统频率响应特性曲线的编程方法。

例6-4 已知离散系统的系统函数为

0.10.2z2+0.2z40.1z6H(z)

10.6z20.4z40.04z6绘出系统的零极点分布图，系统在0~频率范围内的绝对幅频响应、相对幅频响应、相位频率响应和群延迟，并观察零极点分布与系统幅频响应的关系。 Matlab参考程序如下： b=[0.1 0 -0.2 0 0.2 0 -0.1]; a=[1 0 0.6 0 0.4 0 0.04]; 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 %绘制系统的绝对幅频响应曲线

15

axis([0 1 1.1*min(abs(H)) 1.1*max(abs(H))]); title('绝对幅频响应');

subplot(222),plot(w/pi,angle(H));grid %绘制系统的相频响应曲线 axis([0 1 1.1*min(angle(H)) 1.1*max(angle(H))]); title('相频响应');

db=20*log10(abs(H)); %求系统的对数幅频响应 subplot(223),plot(w/pi,db);grid axis([0 1 -100 5]);

title('相对幅频响应（dB）');

grd=grpdelay(b,a,w); %求系统的群延迟 subplot(224), plot(w/pi,grd);grid title('群延迟');

程序运行结果如图6-6所示。 z =

1.0000 0.8660 + 0.5000i 0.8660 - 0.5000i -0.8660 + 0.5000i -0.8660 - 0.5000i -1.0000 p =

0.4150 + 0.6435i

16

0.4150 - 0.6435i -0.4150 + 0.6435i -0.4150 - 0.6435i 0 + 0.3411i 0 - 0.3411i

零极点分布图10.80.60.4traP 0.2yrani0gam-0.2I-0.4-0.6-0.8-1-1-0.5Real Part00.51

17

绝对幅频响应0.820.60.40.20000-20.5相对幅频响应（dB）160相频响应0.5群延迟 1

4-502-10000.51000.51图6-6 离散系统的频域特性

由图6-6可看出，在0~频率范围内，零点集中在0和附近，所以在0和附近对应着幅频响应的谷值，极点集中在/2附近，所以在

/2附近对应着幅频响应的峰值，在绝对幅频响应和相对幅频响应图

中都应证了这一点。

实验内容及步骤

1. 某离散LSI系统的差分方程表示式为

y(n)1.5y(n1)0.5y(n2)x(n)

满足初始条件y(1)4，y(2)10，求系统的单位响应，单位阶跃响应，用filter子函数求系统输入为x(n)(0.25)nu(n)时的零输入、零状态及全响应。

提示：通过解差分方程，可以得到全响应为

18

1112y(n)[()n()n]u(n)，使用filter子函数对系统差分方程进行求解，

2343同时将求解结果与理论计算的结果进行比较。

2. 一简单的数字微分器

y(n)x(n)x(n1)

它计算输入序列的后向一阶差分。当数字微分器的输入为以下序 列时，计算、画出输出序列，并对该数字微分器的实用性进行评价。

（1）矩形脉冲：x1(n)5[u(n)u(n20)]

（2）三角脉冲：x2(n)n[u(n)u(n10)](20n)[u(n10)u(n20)] （3）正弦脉冲：x3(n)sin(n25)[u(n)u(n100)]

3. 已知离散系统的系统函数为

0.20.1z10.3z20.1z30.2z4H(z)

11.1z11.5z20.7z30.3z4求该系统的零极点及零极点分布图，并判断系统的因果性和稳定性。

4. 已知离散系统的系统函数为

0.10.4z2+0.4z4+0.1z6H(z) 24610.3z0.6z0.2z绘出系统的零极点分布图，系统在0~频率范围内的绝对幅频响应、相对幅频响应、相位频率响应和群延迟，并观察零极点分布与系统幅频响应的关系。

5. 观察系统零极点的位置对幅频响应的影响。已知一阶离散系统的系统函数为H(z)

zz1，（1）假设系统的零点在原点，极点分别取zp119

0.2、0.5、0.8，比较它们的幅频响应曲线，（2）假设系统的极点在原点，零点分别取0.2、0.5、0.8，比较它们的幅频响应曲线，从中总结零极点位置对幅频响应的影响。 思考题

1) 卷积和的长度与做卷积和的序列长度有何关系？如何确定卷积和的起始点和终止点？

2) 系统在原点处的零极点对系统的幅频响应有何影响？为什么？ 3) 利用原理中提到的函数分析教材中的全通滤波器和梳状滤波器。 五、实验报告要求

1) 简述实验目的和实验原理。

2) 按实验内容和步骤给出程序的源代码和运行结果，并对结果加以理论说明。

3) 总结实验中的主要结论，你的收获和体会。

20