信号与系统 分析实验报告
实验项目名称：离散线性时不变系统分析；
连续时间系统分析
所属课程名称：信号与系统实验教程 实 验 类 型 ：验证型 指 导 教 师 ：
实 验 日 期 ：2013.06.04 班 级 ： 学 号 ： 姓 名 ：
离散线性时不变系统分析
一、实验目的
1. 掌握离散线性时不变系统的单位序列响应、单位阶跃响应和任意激励下响应的MATLAB 求解方法。
2. 掌握离散线性时不变系统的频域分析方法； 3. 掌握离散线性时不变系统的复频域分析方法；
4. 掌握离散线性时不变系统的零极点分布与系统特性的关系。
二、实验原理及方法
1.
离散线性时不变系统的时域分析
描述一个N 阶线性时不变离散时间系统的数学模型是线性常系统差分方程，N 阶线性时不变离散系统的差分方程一般形式为
（2.1） N M
k =0
∑a k y (n -k ) =∑b i x (n -i )
i =0
也可用系统函数来表示
Y (z )
H (z ) ==
X (z )
∑b z
i
M
-i
（2.2）
∑a z
k k =0
i =0N
-k
b (z ) b 0+b 1z -1+b 2z -2+ +b M z -M ==
a (z ) 1+a 1z -1+a 2z -2+ +a N z -N
系统函数H (z ) 反映了系统响应和激励间的关系。一旦上式中a k ，b i 的数据确定了，系统的性质也就确定了。特别注意a 0必须进行归一化处理，即a 0=1。
对于复杂信号激励下的线性系统，可以将激励信号在时域中分解为单位序列或单位阶跃序列的线性叠加，把这些单元激励信号分别加于系统求其响应，然后把这些响应叠加，即可得到复杂信号作用于系统的零状态响应。因此，求解系统的单位序列响应和单位阶跃响应尤为重要。由图2-1可以看出一个离散LSI 系统响应与激励的关系。
(n ) x (n )*h (n )
(z ) =X (z ) H (z )
图2-1 离散LSI 系统响应与激励的关系
(1) 单位序列响应（单位响应）
单位响应h (n ) 是指离散线性时不变系统在单位序列δ(n ) 激励下的零状态响应，因此
h (n ) 满足线性常系数差分方程(2.1)及零初始状态，即
∑a h (n -k ) =∑b δ(n -i )
k
i
k =0
i =0
N M
， h (-1) =h (-2) = =0 (2.3)
按照定义，它也可表示为
h (n ) =h (n ) *δ(n ) (2.4) 对于离散线性时不变系统，若其输入信号为x (n ) ，单位响应为h (n ) ，则其零状态响应
y zs (n ) 为
y zs (n ) =x (n )*h (n ) (2.5)
可见，h (n ) 能够刻画和表征系统的固有特性，与何种激励无关。一旦知道了系统的单位响应h (n ) ，就可求得系统对任何输入信号x (n ) 所产生的零状态响应y (n ) 。
zs
MATLAB 提供了专门用于求离散系统冲激响应的函数impz()，其调用格式有 [h,n]=impz(b,a)
求解离散系统的单位响应，其中b =[b , b , b , M b ，]a =[1, a 1, a 2, , a N ]，012,
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 ) 是指离散离散时不变系统在单位阶跃序列u (n ) 激励下的零状态响应，它可以表示为
n
s (n ) =u (n ) *h (n ) =
m =-∞
∑h (m ) (2.6)
上式表明，离散线性时不变系统的单位阶跃响应是单位响应的累加和，系统的单位阶跃响应和系统的单位响应之间有着确定的关系，因此，单位阶跃响应也能完全刻画和表征一个线性时不变系统。
MATLAB 提供了专门用于求离散系统单位阶跃响应的函数stepz( )，其调用格式有 [s,n]=stepz(b,a) ：求解离散系统的单位阶跃响应，其中，b =[b , b , b , , b ]
1
2
M
a =[1, a 1, a 2, , a N ]，n =[0,1,2, ]'；
[s,n]=stepz(b,a,N) ：求解离散系统的单位阶跃响应，采样点数由N 确定，n =[0,1,2, , N-1]'；
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 ) ，则系统的零状态响应y (n ) 可由（2.5）式求解。Matlab 提
zs 供了专门用于求离散系统卷积和的函数conv( )，其调用格式有
y=conv(x,h) ：求解序列x ，h 的卷积和，若序列x 的长度为n1，序列h 的长度为n2，卷积和y 的长度为n1+n2-1。这一点需要特别注意，否则，作图时容易造成横纵坐标长度不匹配。
（4）带初始状态的任意激励下的全响应
任意激励下的离散线性时不变系统的全响应为零输入响应和零状态响应之和，表示为
y (n ) =y zi (n ) +y zs (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)：将初始状态转换为初始条件，其中x =[x (-1), x (-2), x (-3), , x (-m )]，
y =[y (-1), y (-2), y (-3), , y (-n )]；
z = filtic(b,a,):将初始状态转换为初始条件x =0, y =[y (-1), y (-2), y (-3), , 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）系统的零极点分布与系统因果性和稳定性的关系
因果系统的单位响应h (n ) 一定满足当n
，对照z 变换定义，系统稳定要求收敛域包含单位圆。系统稳定要求 ∞
n =-∞
∑|h (n ) |
如果系统因果且稳定，收敛域包含∞点和单位圆，那么收敛域可表示为：
r
MATLAB 提供了用于求系统零极点的函数，其调用格式有 roots()：利用多项式求根函数来确定系统函数的零极点位置；
roots(a)：求极点位置，a 为系统函数H (z ) 分母多项式所构成的系数向量； roots(b)：求零点位置，b 为系统函数H (z ) 分子多项式所构成的系数向量；
zplane(b,a)：绘制由行向量b 和a 构成的系统函数的零极点分布图；
zplane(z,p)：绘制由列向量z 确定的零点、列向量p 确定的极点构成的零极点分布图。
（3）系统的零极点分布与系统频率响应的关系 将式(2.2)因式分解，得到
(2.9) M
=1
H (z ) =A m N
∏(1-d k z -1)
k =1
∏(1-c m z -1)
式中，A =b a ，c 是H (z ) 的零点，d 是其极点。A 参数影响频率响应的幅度大小，
m k 00影响系统特性的是零点c 和极点d 的分布。下面采用几何方法研究系统零极点分布对系统m k 频率特性的影响。
将式(2.9)的分子、分母同乘以z N +M ，得到：
=1
H (z ) =A m N
k =1
∏(1-c
M
m
z )
=Az N -M
z -1)
-1
∏(z -c
m =1N k =1
M
(2.10)
m
) )
∏(1-d
k
∏(z -d
j ω
k
假设系统稳定，将z =e j ω代入上式，得到频率响应
H (e j ω) =Ae j (N -M )
设N =M ，由式(2.11)得到
(e ∏ω
m =1N k =1
M
(2.11)
-c m )
=H (e j ω) e j arg[H (e
k
j ω
)]
∏(e ω-d
j
)
=1
H (e j ω) =A m N
k =1
∏(e
N
j ω
(2.12)
-c m )
k
∏(e ω-d
j
)
在z 平面上，示，同样和
e j ω-c m
用一根由零点c 指向单位圆（e j ω）上任一点B 的向量表
c m m
。
e j ω-d k
和
用一根由极点d 指向单位圆（e j ω）上任一点B 的向量d k c m k
d k 分别称为零点矢量和极点矢量，用极坐标表示为：
c m =c m e j αm d k =d k e j βk
将
c m d k 表示式代入式(2.12)，得到
=1
H (e j ω) =A m N
k =1
∏(e
N
j ω
-c m )
k
∏(e ω-d
j j ω
=1
=A m N
∏c
k =1
N
m
=H (e j ω) e j ϕ(ω) B
)
∏d
k
|H (e ) |=A
∏c
m =1N k =1
N
(2.13)
m
∏d
N
N
m =1
k
ϕ(ω) =∑αm -∑βk
k =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)：用于对
H (e j ω) 在[0，Fs/2]上等间隔采样n 点，采样点频率及相应的频响值分别存放在f 和H 中。
H = freqz(b,a,w)：用于对
H (e j ω) 在[0，2π]上进行采样，采样频率点由w 指定。 H = freqz(b,a,f,Fs)：用于对
H (e j ω) 在[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 系统的
h [n ]=δ[n -5]，输入x [n ]=u [n ]-u [n -5]，求系统输出y [n ]=x [n ]*h [n ]:
主程序如下： 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(number of arguments input)；“ ~=”表示不相等
disp('Usage: Y=stepseq(n0,n1,n2)');
return ;
else if((n0n2)|(n1>n2))
error('arguments must satisfy n1n2或者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 的横坐标范围都可以实不相等的
思考题：若解：因为
y [n ]=x [n ]*δ[n -n ]
, 试写出
y [n ]与x [n ]的关系，并对MA TLAB 的仿真结果。
δ[n -n 0]函数是延时器，所以y[n]=x[n]* δ[n -n 0]=x[n -n 0]。
⎧n ⎪
2. 设某线性时不变系统的h[n]=⎨
⎪0⎩⎧10≤n ≤5⎪
输入信号为x[n]=⎨
⎪0其它⎩
0≤n ≤5
其它
求输出：a ：y [n ]=x [n ]*h [n ]; b：y [n ]=x [n ]*h [n +5]
1
2
（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;
subplot(3,1,1);stem(n,x);axis([-10 20 0 2]);title('Input sequence'); ylabel('x[n]');
subplot(3,1,2);stem(n,h);axis([-10 20 0 6]);title('Inpulse Response'); ylabel('h[n+5]'); %output response [y2,ny]=conv_m(x,n,h,n);
subplot(3,1,3);stem(ny,y2);title('Output Sequence');xlabel('n' ); ylabel('y_2[n]'); stepseq.m 的源程序如下：
function [x,n]=stepseq(n0,n1,n2); if nargin ~=3
disp('Usage:Y=stepseq(n0,n1,n2)'); elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1=0];
conv_m.m的源程序如下：
function [y,cy]=conv_m(x,nx,h,nh) if nargin~=4
disp('Usage:Y=conv_m(x,nx,h,nh)'); return ; end ;
nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(x)); ny=[nyb:nye];
y=conv(x,h);
3. 设h [n ]=(0. 9) n u [n ], 输入x [n ]=u [n ]-u [n -10]，求系统输出y [n ]=x [n ]*h [n ]:
n=-5:50;
u1=stepseq(0,-5,50);u2=stepseq(10,-5,50);%u1=u[n];u2=u[n-10] %input x[n] x=u1-u2;
%impluse response h[n] h=((0.9).^n).*u1;
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 [y,ny]=conv_m(x,n,h,n);
subplot(3,1,3);stem(ny,y);title('Output Sequence');xlabel('n' ); ylabel('y[n]');
stepseq.m 的源程序如下：
function [x,n]=stepseq(n0,n1,n2); if nargin ~=3
disp('Usage:Y=stepseq(n0,n1,n2)'); elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1=0];
conv_m.m的源程序如下：
function [y,cy]=conv_m(x,nx,h,nh) if nargin~=4
disp('Usage:Y=conv_m(x,nx,h,nh)'); return ; end ;
nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(x)); ny=[nyb:nye]; y=conv(x,h);
思考题：在此问题中h [n ]为无限长序列，在计算卷积和时应如何处理？
解：取尽可能多的包含能量分布集中的具有代表性的点
4某离散时间线性时不变系统的差分方程为：y [n ]-0. 8[n -1]=x [n ]
a. 计算n =[-10:50]时的系统冲激响应h [n ]
b. 计算n =[-10:50]时的系统阶跃响应s [n ]
a=[1,-0.8];b=1;n=[-10:50];
x=impseq(0,-10,50);
h=filter(b,a,x);%MATLAB固有函数filter 表示系统方程ay=bx；x 是冲击信号，则h 就是冲击响应 subplot(2,1,1);stem(n,h);
title('Inpulse Response');xlabel('n' );ylabel('h[n]');
x=stepseq(0,-10,50);
s=filter(b,a,x); %x是阶跃信号，则s 就是阶跃响应
subplot(2,1,2);stem(n,s);
title('Step Response');xlabel('n' );ylabel('s[n]');
impseq.m 的源程序如下：
function [x,n]=impseq(n0,n1,n2)
%Generates x[n]=delta(n-n0)；n1
%[x,n]=impseq(n0,n1,n2)＇) ；
if nargin~=3
disp('Usage:Y=impseq(n0,n1,n2)');
return ;
elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1
end
n=[n1:n2];
x=[(n-n0)==0];
stepseq.m 的源程序如下：
function [x,n]=stepseq(n0,n1,n2);
if nargin ~=3
disp('Usage:Y=stepseq(n0,n1,n2)');
elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1
end
n=[n1:n2];
x=[(n-n0)>=0];
连续时间系统分析
1. 设有两个稳定的LTI 系统，分别可由下列微分方程来描述：
dy (t ) +3y (t ) =3x (t ) ； dt
d 2y (t ) dy (t ) d 2x (t ) 3+4+y (t ) =+5x (t ) 。 22dt dt dt
请分别画出它们的系统频率响应的幅值和相位特性曲线。
lab41a.m 的源程序如下：
a=[1 3]; % y's coeff. in the system equation
b=3; % x's coeff.in the system equation
freqs(b,a); % freqs 为MATLAB 中求连续时间系统拉普拉斯变换域（s 域）的频率响应的固有函数
Lab41b.m 的源程序如下：
a=[3 4 1]; % 系数矩阵
b=[1 0 5]; % 系数矩阵
freqs(b,a); % freqs为MATLAB 中求连续时间系统拉普拉斯变换域（s 域）的频率响应的固有函数
思考题：
①求出系统a 的频率响应表达式；
H (jw ) =3 jw +3
②根据表达式，用对数坐标画出系统的幅频响应曲线以及相频响应曲线；
w=linspace(0,5,200);
b=[3];a=[1 3];
H=freqs(b,a,w);
subplot(2,1,1);plot(w,abs(H));
set(gca,'xtick' ,[0 1 2 3 4 5]);
set(gca,'ytick' ,[0 0.4 0.707 1]);
title('幅值谱|H(\omega)|');
xlabel('\omega(rad/s)');
ylabel(' 幅值' );grid on ;
subplot(2,1,2);plot(w,angle(H));
set(gca,'xtick' ,[0 1 2 3 4 5]);
title(' 相位谱' );
xlabel('\omega(rad/s)');
ylabel(' 相位' );grid on ;
运行所得到的波形图如下：
③与freqs 函数得到的结果进行对比。
与freqs 所得到的函数图像相比较，可以得出两个函数的幅值谱和相位谱的图形刚好相反对称。
0. 2s 2+0. 3s +12．有一模拟滤波器，其传递函数为H (s ) =，应用freqs 函数画出它的幅频2s +0. 4s +1
特性和相频特性曲线。
Lab42.m 的源程序如下：
a=[1 0.4 1]; % 系数矩阵
b=[0.2 0.3 1]; % 系数矩阵
freqs(b,a); % freqs 为MATLAB 中求连续时间系统拉普拉斯变换域（s 域）的频率响应的固有函数
四、实验小结
此次试验是对分析，通过实验对离散LSI 系统的单位序列响应、单位阶跃响应和任意激励下响应、系统的频域分析、复频域分析、零极点分布与系统特性的MATLAB 求解方法都有所了解并用MATLAB 提供的函数来分析。
实验中通过读和学习老师给的例题，然后仿照例题，完成了实验要求的任务。MAYLAB 功能强大，是分析系统的有力工具。
此次实验用到的函数很多，函数调用方法也很多样，一时没法全掌握，但我已全部理解，并知道如何去使用。通过实验，我对离散性时不变系统一些特性都有所掌握，对零极点分布对系统特性的影响也基本理解清楚。只有当极点在单位圆内，系统才稳定，此时如果零点在圆外，为全通滤波器，零点在圆上，为梳妆滤波器，零点在圆内，为最小相位延时系统。
信号与系统 分析实验报告
实验项目名称：离散线性时不变系统分析；
连续时间系统分析
所属课程名称：信号与系统实验教程 实 验 类 型 ：验证型 指 导 教 师 ：
实 验 日 期 ：2013.06.04 班 级 ： 学 号 ： 姓 名 ：
离散线性时不变系统分析
一、实验目的
1. 掌握离散线性时不变系统的单位序列响应、单位阶跃响应和任意激励下响应的MATLAB 求解方法。
2. 掌握离散线性时不变系统的频域分析方法； 3. 掌握离散线性时不变系统的复频域分析方法；
4. 掌握离散线性时不变系统的零极点分布与系统特性的关系。
二、实验原理及方法
1.
离散线性时不变系统的时域分析
描述一个N 阶线性时不变离散时间系统的数学模型是线性常系统差分方程，N 阶线性时不变离散系统的差分方程一般形式为
（2.1） N M
k =0
∑a k y (n -k ) =∑b i x (n -i )
i =0
也可用系统函数来表示
Y (z )
H (z ) ==
X (z )
∑b z
i
M
-i
（2.2）
∑a z
k k =0
i =0N
-k
b (z ) b 0+b 1z -1+b 2z -2+ +b M z -M ==
a (z ) 1+a 1z -1+a 2z -2+ +a N z -N
系统函数H (z ) 反映了系统响应和激励间的关系。一旦上式中a k ，b i 的数据确定了，系统的性质也就确定了。特别注意a 0必须进行归一化处理，即a 0=1。
对于复杂信号激励下的线性系统，可以将激励信号在时域中分解为单位序列或单位阶跃序列的线性叠加，把这些单元激励信号分别加于系统求其响应，然后把这些响应叠加，即可得到复杂信号作用于系统的零状态响应。因此，求解系统的单位序列响应和单位阶跃响应尤为重要。由图2-1可以看出一个离散LSI 系统响应与激励的关系。
(n ) x (n )*h (n )
(z ) =X (z ) H (z )
图2-1 离散LSI 系统响应与激励的关系
(1) 单位序列响应（单位响应）
单位响应h (n ) 是指离散线性时不变系统在单位序列δ(n ) 激励下的零状态响应，因此
h (n ) 满足线性常系数差分方程(2.1)及零初始状态，即
∑a h (n -k ) =∑b δ(n -i )
k
i
k =0
i =0
N M
， h (-1) =h (-2) = =0 (2.3)
按照定义，它也可表示为
h (n ) =h (n ) *δ(n ) (2.4) 对于离散线性时不变系统，若其输入信号为x (n ) ，单位响应为h (n ) ，则其零状态响应
y zs (n ) 为
y zs (n ) =x (n )*h (n ) (2.5)
可见，h (n ) 能够刻画和表征系统的固有特性，与何种激励无关。一旦知道了系统的单位响应h (n ) ，就可求得系统对任何输入信号x (n ) 所产生的零状态响应y (n ) 。
zs
MATLAB 提供了专门用于求离散系统冲激响应的函数impz()，其调用格式有 [h,n]=impz(b,a)
求解离散系统的单位响应，其中b =[b , b , b , M b ，]a =[1, a 1, a 2, , a N ]，012,
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 ) 是指离散离散时不变系统在单位阶跃序列u (n ) 激励下的零状态响应，它可以表示为
n
s (n ) =u (n ) *h (n ) =
m =-∞
∑h (m ) (2.6)
上式表明，离散线性时不变系统的单位阶跃响应是单位响应的累加和，系统的单位阶跃响应和系统的单位响应之间有着确定的关系，因此，单位阶跃响应也能完全刻画和表征一个线性时不变系统。
MATLAB 提供了专门用于求离散系统单位阶跃响应的函数stepz( )，其调用格式有 [s,n]=stepz(b,a) ：求解离散系统的单位阶跃响应，其中，b =[b , b , b , , b ]
1
2
M
a =[1, a 1, a 2, , a N ]，n =[0,1,2, ]'；
[s,n]=stepz(b,a,N) ：求解离散系统的单位阶跃响应，采样点数由N 确定，n =[0,1,2, , N-1]'；
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 ) ，则系统的零状态响应y (n ) 可由（2.5）式求解。Matlab 提
zs 供了专门用于求离散系统卷积和的函数conv( )，其调用格式有
y=conv(x,h) ：求解序列x ，h 的卷积和，若序列x 的长度为n1，序列h 的长度为n2，卷积和y 的长度为n1+n2-1。这一点需要特别注意，否则，作图时容易造成横纵坐标长度不匹配。
（4）带初始状态的任意激励下的全响应
任意激励下的离散线性时不变系统的全响应为零输入响应和零状态响应之和，表示为
y (n ) =y zi (n ) +y zs (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)：将初始状态转换为初始条件，其中x =[x (-1), x (-2), x (-3), , x (-m )]，
y =[y (-1), y (-2), y (-3), , y (-n )]；
z = filtic(b,a,):将初始状态转换为初始条件x =0, y =[y (-1), y (-2), y (-3), , 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）系统的零极点分布与系统因果性和稳定性的关系
因果系统的单位响应h (n ) 一定满足当n
，对照z 变换定义，系统稳定要求收敛域包含单位圆。系统稳定要求 ∞
n =-∞
∑|h (n ) |
如果系统因果且稳定，收敛域包含∞点和单位圆，那么收敛域可表示为：
r
MATLAB 提供了用于求系统零极点的函数，其调用格式有 roots()：利用多项式求根函数来确定系统函数的零极点位置；
roots(a)：求极点位置，a 为系统函数H (z ) 分母多项式所构成的系数向量； roots(b)：求零点位置，b 为系统函数H (z ) 分子多项式所构成的系数向量；
zplane(b,a)：绘制由行向量b 和a 构成的系统函数的零极点分布图；
zplane(z,p)：绘制由列向量z 确定的零点、列向量p 确定的极点构成的零极点分布图。
（3）系统的零极点分布与系统频率响应的关系 将式(2.2)因式分解，得到
(2.9) M
=1
H (z ) =A m N
∏(1-d k z -1)
k =1
∏(1-c m z -1)
式中，A =b a ，c 是H (z ) 的零点，d 是其极点。A 参数影响频率响应的幅度大小，
m k 00影响系统特性的是零点c 和极点d 的分布。下面采用几何方法研究系统零极点分布对系统m k 频率特性的影响。
将式(2.9)的分子、分母同乘以z N +M ，得到：
=1
H (z ) =A m N
k =1
∏(1-c
M
m
z )
=Az N -M
z -1)
-1
∏(z -c
m =1N k =1
M
(2.10)
m
) )
∏(1-d
k
∏(z -d
j ω
k
假设系统稳定，将z =e j ω代入上式，得到频率响应
H (e j ω) =Ae j (N -M )
设N =M ，由式(2.11)得到
(e ∏ω
m =1N k =1
M
(2.11)
-c m )
=H (e j ω) e j arg[H (e
k
j ω
)]
∏(e ω-d
j
)
=1
H (e j ω) =A m N
k =1
∏(e
N
j ω
(2.12)
-c m )
k
∏(e ω-d
j
)
在z 平面上，示，同样和
e j ω-c m
用一根由零点c 指向单位圆（e j ω）上任一点B 的向量表
c m m
。
e j ω-d k
和
用一根由极点d 指向单位圆（e j ω）上任一点B 的向量d k c m k
d k 分别称为零点矢量和极点矢量，用极坐标表示为：
c m =c m e j αm d k =d k e j βk
将
c m d k 表示式代入式(2.12)，得到
=1
H (e j ω) =A m N
k =1
∏(e
N
j ω
-c m )
k
∏(e ω-d
j j ω
=1
=A m N
∏c
k =1
N
m
=H (e j ω) e j ϕ(ω) B
)
∏d
k
|H (e ) |=A
∏c
m =1N k =1
N
(2.13)
m
∏d
N
N
m =1
k
ϕ(ω) =∑αm -∑βk
k =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)：用于对
H (e j ω) 在[0，Fs/2]上等间隔采样n 点，采样点频率及相应的频响值分别存放在f 和H 中。
H = freqz(b,a,w)：用于对
H (e j ω) 在[0，2π]上进行采样，采样频率点由w 指定。 H = freqz(b,a,f,Fs)：用于对
H (e j ω) 在[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 系统的
h [n ]=δ[n -5]，输入x [n ]=u [n ]-u [n -5]，求系统输出y [n ]=x [n ]*h [n ]:
主程序如下： 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(number of arguments input)；“ ~=”表示不相等
disp('Usage: Y=stepseq(n0,n1,n2)');
return ;
else if((n0n2)|(n1>n2))
error('arguments must satisfy n1n2或者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 的横坐标范围都可以实不相等的
思考题：若解：因为
y [n ]=x [n ]*δ[n -n ]
, 试写出
y [n ]与x [n ]的关系，并对MA TLAB 的仿真结果。
δ[n -n 0]函数是延时器，所以y[n]=x[n]* δ[n -n 0]=x[n -n 0]。
⎧n ⎪
2. 设某线性时不变系统的h[n]=⎨
⎪0⎩⎧10≤n ≤5⎪
输入信号为x[n]=⎨
⎪0其它⎩
0≤n ≤5
其它
求输出：a ：y [n ]=x [n ]*h [n ]; b：y [n ]=x [n ]*h [n +5]
1
2
（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;
subplot(3,1,1);stem(n,x);axis([-10 20 0 2]);title('Input sequence'); ylabel('x[n]');
subplot(3,1,2);stem(n,h);axis([-10 20 0 6]);title('Inpulse Response'); ylabel('h[n+5]'); %output response [y2,ny]=conv_m(x,n,h,n);
subplot(3,1,3);stem(ny,y2);title('Output Sequence');xlabel('n' ); ylabel('y_2[n]'); stepseq.m 的源程序如下：
function [x,n]=stepseq(n0,n1,n2); if nargin ~=3
disp('Usage:Y=stepseq(n0,n1,n2)'); elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1=0];
conv_m.m的源程序如下：
function [y,cy]=conv_m(x,nx,h,nh) if nargin~=4
disp('Usage:Y=conv_m(x,nx,h,nh)'); return ; end ;
nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(x)); ny=[nyb:nye];
y=conv(x,h);
3. 设h [n ]=(0. 9) n u [n ], 输入x [n ]=u [n ]-u [n -10]，求系统输出y [n ]=x [n ]*h [n ]:
n=-5:50;
u1=stepseq(0,-5,50);u2=stepseq(10,-5,50);%u1=u[n];u2=u[n-10] %input x[n] x=u1-u2;
%impluse response h[n] h=((0.9).^n).*u1;
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 [y,ny]=conv_m(x,n,h,n);
subplot(3,1,3);stem(ny,y);title('Output Sequence');xlabel('n' ); ylabel('y[n]');
stepseq.m 的源程序如下：
function [x,n]=stepseq(n0,n1,n2); if nargin ~=3
disp('Usage:Y=stepseq(n0,n1,n2)'); elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1=0];
conv_m.m的源程序如下：
function [y,cy]=conv_m(x,nx,h,nh) if nargin~=4
disp('Usage:Y=conv_m(x,nx,h,nh)'); return ; end ;
nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(x)); ny=[nyb:nye]; y=conv(x,h);
思考题：在此问题中h [n ]为无限长序列，在计算卷积和时应如何处理？
解：取尽可能多的包含能量分布集中的具有代表性的点
4某离散时间线性时不变系统的差分方程为：y [n ]-0. 8[n -1]=x [n ]
a. 计算n =[-10:50]时的系统冲激响应h [n ]
b. 计算n =[-10:50]时的系统阶跃响应s [n ]
a=[1,-0.8];b=1;n=[-10:50];
x=impseq(0,-10,50);
h=filter(b,a,x);%MATLAB固有函数filter 表示系统方程ay=bx；x 是冲击信号，则h 就是冲击响应 subplot(2,1,1);stem(n,h);
title('Inpulse Response');xlabel('n' );ylabel('h[n]');
x=stepseq(0,-10,50);
s=filter(b,a,x); %x是阶跃信号，则s 就是阶跃响应
subplot(2,1,2);stem(n,s);
title('Step Response');xlabel('n' );ylabel('s[n]');
impseq.m 的源程序如下：
function [x,n]=impseq(n0,n1,n2)
%Generates x[n]=delta(n-n0)；n1
%[x,n]=impseq(n0,n1,n2)＇) ；
if nargin~=3
disp('Usage:Y=impseq(n0,n1,n2)');
return ;
elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1
end
n=[n1:n2];
x=[(n-n0)==0];
stepseq.m 的源程序如下：
function [x,n]=stepseq(n0,n1,n2);
if nargin ~=3
disp('Usage:Y=stepseq(n0,n1,n2)');
elseif ((n0n2)|(n1>n2))
error('arguments must satisfy n1
end
n=[n1:n2];
x=[(n-n0)>=0];
连续时间系统分析
1. 设有两个稳定的LTI 系统，分别可由下列微分方程来描述：
dy (t ) +3y (t ) =3x (t ) ； dt
d 2y (t ) dy (t ) d 2x (t ) 3+4+y (t ) =+5x (t ) 。 22dt dt dt
请分别画出它们的系统频率响应的幅值和相位特性曲线。
lab41a.m 的源程序如下：
a=[1 3]; % y's coeff. in the system equation
b=3; % x's coeff.in the system equation
freqs(b,a); % freqs 为MATLAB 中求连续时间系统拉普拉斯变换域（s 域）的频率响应的固有函数
Lab41b.m 的源程序如下：
a=[3 4 1]; % 系数矩阵
b=[1 0 5]; % 系数矩阵
freqs(b,a); % freqs为MATLAB 中求连续时间系统拉普拉斯变换域（s 域）的频率响应的固有函数
思考题：
①求出系统a 的频率响应表达式；
H (jw ) =3 jw +3
②根据表达式，用对数坐标画出系统的幅频响应曲线以及相频响应曲线；
w=linspace(0,5,200);
b=[3];a=[1 3];
H=freqs(b,a,w);
subplot(2,1,1);plot(w,abs(H));
set(gca,'xtick' ,[0 1 2 3 4 5]);
set(gca,'ytick' ,[0 0.4 0.707 1]);
title('幅值谱|H(\omega)|');
xlabel('\omega(rad/s)');
ylabel(' 幅值' );grid on ;
subplot(2,1,2);plot(w,angle(H));
set(gca,'xtick' ,[0 1 2 3 4 5]);
title(' 相位谱' );
xlabel('\omega(rad/s)');
ylabel(' 相位' );grid on ;
运行所得到的波形图如下：
③与freqs 函数得到的结果进行对比。
与freqs 所得到的函数图像相比较，可以得出两个函数的幅值谱和相位谱的图形刚好相反对称。
0. 2s 2+0. 3s +12．有一模拟滤波器，其传递函数为H (s ) =，应用freqs 函数画出它的幅频2s +0. 4s +1
特性和相频特性曲线。
Lab42.m 的源程序如下：
a=[1 0.4 1]; % 系数矩阵
b=[0.2 0.3 1]; % 系数矩阵
freqs(b,a); % freqs 为MATLAB 中求连续时间系统拉普拉斯变换域（s 域）的频率响应的固有函数
四、实验小结
此次试验是对分析，通过实验对离散LSI 系统的单位序列响应、单位阶跃响应和任意激励下响应、系统的频域分析、复频域分析、零极点分布与系统特性的MATLAB 求解方法都有所了解并用MATLAB 提供的函数来分析。
实验中通过读和学习老师给的例题，然后仿照例题，完成了实验要求的任务。MAYLAB 功能强大，是分析系统的有力工具。
此次实验用到的函数很多，函数调用方法也很多样，一时没法全掌握，但我已全部理解，并知道如何去使用。通过实验，我对离散性时不变系统一些特性都有所掌握，对零极点分布对系统特性的影响也基本理解清楚。只有当极点在单位圆内，系统才稳定，此时如果零点在圆外，为全通滤波器，零点在圆上，为梳妆滤波器，零点在圆内，为最小相位延时系统。