离散LSI系统分析

信号与系统分析实验报告

实验项目名称：离散线性时不变系统分析；

连续时间系统分析

所属课程名称：信号与系统实验教程

实验类型：验证型

指导教师：

实验日期：2013.06.04

班级：

学号：

姓名：

离散线性时不变系统分析

一、实验目的

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

2. 掌握离散线性时不变系统的频域分析方法；

3. 掌握离散线性时不变系统的复频域分析方法；

4. 掌握离散线性时不变系统的零极点分布与系统特性的关系。

二、实验原理及方法

1.

离散线性时不变系统的时域分析

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

)

()(0

i n x b k n y a M

i i N k k -=-∑∑== （2.1） 也可用系统函数来表示

12001212120

()

()()()

()1M

i

M i

i M N

N

k

N k k b z

b b z b z b z Y z b z H z X z a z a z a z a z a z

----=----=++++==

==

++++∑∑ （2.2）

系统函数()H z 反映了系统响应和激励间的关系。一旦上式中k a ，i b 的数据确定了，系统的性质也就确定了。特别注意0a 必须进行归一化处理，即01a =。

对于复杂信号激励下的线性系统，可以将激励信号在时域中分解为单位序列或单位阶跃

序列的线性叠加，把这些单元激励信号分别加于系统求其响应，然后把这些响应叠加，

即可得到复杂信号作用于系统的零状态响应。因此，求解系统的单位序列响应和单位阶跃响应尤为重要。由图2-1可以看出一个离散LSI 系统响应与激励的关系。

()()()

z X z H z =()()*()

n x n h n

图2-1 离散LSI 系统响应与激励的关系

(1) 单位序列响应（单位响应）

单位响应()h n 是指离散线性时不变系统在单位序列()n δ激励下的零状态响应，因此

()h n 满足线性常系数差分方程(2.1)及零初始状态，即

()()

N M

k

i

k i a h n k b n i δ==-=-∑∑， (1)(2)0h h -=-== (2.3)

按照定义，它也可表示为

()()()h n h n n δ=* (2.4) 对于离散线性时不变系统，若其输入信号为()x n ，单位响应为()h n ，则其零状态响应

()zs y n 为

()()*()zs y n x n h n = (2.5)

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

y n 。

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

求解离散系统的单位响应，其中012[,,,,]M b b b b b =，12[1,,,,]N

a a a a =，

[0,1,2,]n '=；

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

求解离散系统的单位响应，采样点数由N 确定，[0,1,2,

,N-1]n '=；

impz(b,a) ：在当前窗口，用stem(n,h)绘出图形。 （2）单位阶跃响应

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

()()()()n

m s n u n h n h m =-∞

=*=

∑ (2.6)

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

MATLAB 提供了专门用于求离散系统单位阶跃响应的函数stepz( )，其调用格式有 [s,n]=stepz(b,a) ：求解离散系统的单位阶跃响应，其中0

1

2

[,,,,]M

b b b b b =，

12[1,,,

,]N a a a a =，[0,1,2,]n '=；

[s,n]=stepz(b,a,N) ：求解离散系统的单位阶跃响应，采样点数由N 确定，[0,1,2,,N-1]n '=；

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

已经知道，离散线性时不变系统可用常系数线性差分方程(2.1)式来描述，Matlab 提供的函数dlsim( )能对上述差分方程描述的离散线性时不变系统的响应进行仿真，该函数不仅能绘制指定时间范围内的系统响应波形图，而且还能求出系统响应的数值解。其调用格式有

dlsim(b,a, x) ：求解输入序列为x 的零状态响应

需要特别强调的是，Matlab 总是把由分子和分母多项式表示的任何系统都当作是因果系统。所以，利用impz (b,a)，stepz(b,a)，d lsim(b,a,x)函数求得的响应总是因果信号。

同时，卷积和也是线性时不变系统求解零状态响应的重要工具之一。假设系统的输入信号为()x n ，单位响应为()h n ，则系统的零状态响应()zs y n 可由（2.5）式求解。Matlab 提

供了专门用于求离散系统卷积和的函数conv( )，其调用格式有

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

（4）带初始状态的任意激励下的全响应

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

()()()zi zs y n y n y n =+ （2.7）

Matlab 提供了用于求离散系统全响应的函数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)：将初始状态转换为初始条件，其中[(1),(2),(3),

,()]x x x x x m =----，

[(1),(2),(3),,()]y y y y y n =----；

z = filtic(b,a,):将初始状态转换为初始条件0x =,[(1),(2),(3),

,()]y y y y y n =----

2 离散线性时不变系统的复频域（Z 域）分析

（1）利用Z 变换解差分方程

在前面图2-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 都是符号表达式；

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 数组，返回有理多项式。 （2）系统的零极点分布与系统因果性和稳定性的关系

因果系统的单位响应)(n h 一定满足当0

系统稳定要求 ∞

<∑∞

-∞

=n n h |)(| ，对照z 变换定义，系统稳定要求收敛域包含单位圆。如果系统因果且稳定，收敛域包含∞点和单位圆，那么收敛域可表示为：

10,||<<∞≤

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

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

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

zplane(z,p)：绘制由列向量z 确定的零点、列向量p 确定的极点构成的零极点分布图。

（3）系统的零极点分布与系统频率响应的关系 将式(2.2)因式分解，得到

∏∏=-=---=N

k k M

m

m z d z c A z H 1

11

1)

1()

1()( (2.9) 式中，00

a b

A =，m c 是)(z H 的零点，k d 是其极点。A 参数影响频率响应的幅度大小，

影响系统特性的是零点

m c 和极点k d 的分布。

下面采用几何方法研究系统零极点分布对系统

频率特性的影响。

将式(2.9)的分子、分母同乘以M N z +，得到：

∏∏∏∏==-=-=---=--=N k k

M

m m

M

N N

k k

M

m m

d

z c

z Az z d

z c

A z H 1

11

11

1

)

()()

1()

1()( (2.10)

假设系统稳定，将ωj e z =代入上式，得到频率响应

)]

(arg[1

1)

()()

()

()(ω

ωωω

ω

ωj e

H j j N k k

j

M

m m j M N j j e e H d

e c e Ae e H =--=∏∏==- (2.11)

设M N =，由式(2.11)得到

∏∏==--=N

k k

j

N

m m j j d

e c e

A e H 1

1

)

()

()(ωω

ω (2.12)

在z 平面上，m

j c e -ω用一根由零点

m c 指向单位圆（ω

j e ）上任一点B 的向量c m 表

示，同样k

j d e -ω用一根由极点

k d 指向单位圆（ωj e ）上任一点B 的向量d k 来表示，

c m 和

B

d k 分别称为零点矢量和极点矢量，用极坐标表示为：

m

j m m e c B c α=，

k

j k k e d B d β=。

将

c m 和

d k 表示式代入式(2.12)，得到

)(1

1

1

1

)()

()

()(ω?ωωω

ωj j N

k k

N

m m

N

k k

j N

m m j j e e H d

c

A d

e c e

A e H ==--=∏∏∏∏====

∏∏===N k k

N

m m

j d

c

A

e H 1

1|)(|ω

(2.13)

∑∑==-=N

k k

N

m m 1

1

)(βαω? (2.14)

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

Matlab 提供了专门用于求离散系统频率响应的函数freqz(),其调用格式如下： [H,w] = freqz(b,a,n)：返回数字系统的n 点频率值（复数），这n 个点均匀地分布在[0，π]上，系统默认的采样点数目为512点；

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

H = freqz(b,a,w)：用于对()j

H eω在[0，2π]上进行采样，采样频率点由w指定。

H = freqz(b,a,f,Fs)：用于对()j

H eω在[0，Fs]上进行采样，采样频率点由f指定。

freqz(b,a,n)：用于在当前图形窗口中绘制幅频和相频特性曲线。

下面介绍几个实用的函数：

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

db=20*log10((mag+eps)/max(mag))：求解系统的相对幅频响应；

pha=angle(H):求解系统的相频响应；

grd=grpdelay(b,a,w): 求解系统的群延迟。

三、实验内容

1.设某LTI系统的

]5

[

]

[-

=n

n

hδ

，输入

]5

[

]

[

]

[-

-

=n

u

n

u

n

x

，求系统输出

]

[

*]

[

]

[n

h

n

x

n

y=:

主程序如下：

clear;clf;clc;

n=[-5:20]; %设定一个n的观察范围

h=delta(n-5);x=stepseq(0,-5,20)-stepseq(5,-5,20);

[y,ny]=conv_m(x,n,h,n)

subplot(2,2,1);stem(n,x);title('x[n]');

xlabel('n');axis([-5,20,0,1.2]);grid on;

subplot(2,2,2);stem(n,h);title('h[n]');

xlabel('n');axis([-5,20,0,1.2]);grid on;

subplot(2,2,3);stem(ny,y);title('y[n]');

xlabel('n');axis([-5,20,0,1.2]);grid on;

delta.m

function y=delta(x)

y=(x==0);

stepseq.m的源程序如下：

function[x,n]=stepseq(n0,n1,n2)

if nargin~=3%nargin(num ber of arguments input)；“ ~=”表示不相等 disp('Usage: Y=stepseq(n0,n1,n2)');

return;

else if((n0n2)|(n1>n2))

error('arguments must satisfy n1<=n0<=n2')%n0n2或者n1>n2 end

n=[n1:n2];

x=[(n-n0)>=0];

conv_m.m的源程序如下：

function[y,ny]=conv_m(x,nx,h,nh)

%Modified convolution routine for signal processing

%[y,ny]=conv_m(x,nx,h,nh)

%y=convolution result

%ny=support of y

%x=first signal on support nx

%nx=support of x

%h=second signal on support nh %nh=support of h if nargin~=4

disp('Usage:Y=conv_m(x,nx,h,nh)'); return ; end ;

nyb=nx(1)+nh(1); %ny ＇s begining nye=nx(length(x))+nh(length(h)); %ny ＇s end

ny=[nyb:nye]; %ny 仅仅为了计算一下结果y 对应的横坐标范围

y=conv(x,h); %Convolution and polynomial multiplication;conv 为MATLAB 的固有函数

%conv_m 函数比conv 函数多用了用于表示横坐标范围的ny 、nx 、nh,因为这里假定二维坐标范围nx 和nh 可能从负数开始，所以要重新计算y 的横坐标范围ny ；如果nx 和nh 均为[0:N]，则可以直接得出ny 为[0:2N](正如conv 函数中那样）；在conv_m 和conv 函数中，x 和h 的横坐标范围都可以实不相等的

思考题：若]

[*][][0

n n n x n y -=δ,试写出

][n y 与][n x 的关系，并对MA TLAB 的仿真结果。

解：因为[]0n n -δ函数是延时器，所以y[n]=x[n]* []0n n -δ=x[n -0n ]。

2.设某线性时不变系统的h[n]=???

?

?≤≤其它

05

n 0n

输入信号为x[n]=??

?

??≤≤其它05n 01

求输出：a ：];[*][][1

n h n x n y = b ：]5[*][][2

+=n h n x n y

（a ）主程序如下：

n=[-5:20];

u1=stepseq(0,-5,20);u2=stepseq(6,-5,20);%u1=u[n]；u2=u[n-6] %input x[n]

x=u1-u2;

%impulse response h[n] h=n.*x;

subplot(3,1,1);stem(n,x);axis([-5 20 0 2]);title('Input Sequence'); ylabel('x[n]');

subplot(3,1,2);stem(n,h);axis([-5 20 0 6]);title('Inpulse Response'); ylabel('h[n]'); %output response

[y1,ny]=conv_m(x,n,h,n);%conv_m 为自定义求卷积的函数

subplot(3,1,3);stem(ny,y1);title('Output Sequence');xlabel('n'); ylabel('y_1[n]');

（b ）主程序如下：

n=-10:20;

u1=stepseq(0,-10,20);u2=stepseq(6,-10,20);%u1=u[n]；u2=u[n-6] %input x[n] x=u1-u2;

u3=stepseq(-5,-10,20);u4=stepseq(1,-10,20);%u3=u[n+5]；u4=u[n-1] x1=u3-u4;

%impulse rsponse h[n+5] h=(n+5).*x1;