实验1 利用matlab进行系统的时域分析

实验1 利用matlab进行系统的时域分析
一．实验目的:
1．了解离散时间序列卷积与的matlab实现;
2．利用卷积与求解系统的零状态响应;
二．实验原理:
1．连续时间系统零状态响应的求解
连续时间LTI系统以常系数微分方程描述,系统的零状态响应可通过求解初始状态为零的微分方程得到。

在MATLAB中,控制系统工具箱提供了一个用于求解零初始状态微分方程数值解的函数lsim。

其调用方式为
y= lsim( sys,x,t)
式中t表示计算系统响应的抽样点向量,x就是系统输入信号向量,sys就是连续时间LTI系统模型,用来表示微分方程、差分方程、状态方程。

在求解微分方程时,微分方程的连续时间LTI系统模型sys要借助tf函数获得,其调用方式为
sys= tf(b,a)
式中b与a分别为微分方程右端与左端各项的系数向量。

例如对3阶微分方程
+++=+++
可用
a=[ a3, a2, a1, a0];b=[b3 ,b2, b1,b0]; sys=tf( b,a)
获得连续时间LTI模型。

注意微分方程中为零的系数一定要写入向量a与b中。

【例2-1】描述某力学系统中物体位移y(t)与外力f(t)的关系为
++y(t)=x(t)
物体质量m=l kg,弹簧的弹性系数ks= 100 N/m,物体与地面的摩擦系数fd=2 N·s/m,系统的初始储能为零,若外力x(t)就是振幅为10、周期为1的正弦信号,求物体的位移y(t)。

解:由已知条件,系统的输入信号为x(t)=10sin(2πt),系统的微分方程为
++100y(t)=x(t)
计算物体位移y(t)的MATLAB程序如下:
%program2_1微分方程求解
ts=0;te=5;dt=0、01;
sys=tf([1],[1 2 100]);
t=ts:dt:te;
x=10*sin(2*pi*t);
y=lsim(sys,x,t);
plot(t,y);
xlabel('Time(sec)')
ylabel('y(t)')
-0.25
-0.2-0.15-0.1-0.050
0.050.10.150.2Time(sec)
y (t )
图2-1系统的零状态响应
2、连续时间系统冲激响应与阶跃响应的求解
在MATLAB 中,求解系统冲激响应可应用控制系统工具箱提供的函数impulse,求解阶跃响应可利用函数step 。

其调用方式为
y= impulse( sys,t) y= step( sys,t)
式中t 表示计算系统响应的抽样点向量,sys 就是连续时间LTI 系统模型。

下面举例说明其应用。

【例2-2】在例2-1所述力学系统中,若外力x(t)就是强度为10的冲激信号,求物体的位移y(t)。

解:由已知条件,系统的输入信号为x(t)=10δ(t),系统的微分方程可写成:
+
+100h(t)=10δ(t)
物体位移y(t)即系统的冲激响应,计算其的MATLAB 程序如下: %program3_2连续时间系统的冲激响应 clear clc
ts=0;te=5;dt=0、01; sys=tf([10],[1 2 100]); t=ts:dt:te;
y=impulse(sys,t); plot(t,y);
xlabel('Time(sec)') ylabel('h(t)')
-0.8
-0.6
-0.4
-0.2
0.2
0.4
0.6
0.8
1
Time(sec)
h
(
t
)
图2-2 连续时间系统的冲激响应
3、离散的时间系统零状态相应的求解
大量的离散时间LTI系统都可以用如下的线性常系数差分方程描述:
其中a0=1,x[k]、y[k]分别表示系统的输入与输出,n就是差分方程的阶数。

已知差分方程的n个初始状态与输入x[k],就可以编程由下式迭代计算出系统的输出:
y[k]=-
在零初始状态下,MATLAB信号处理工具提供了一个filter函数计算由差分方程描述的系统的响应。

其调用方式为:
y= filter(b,a,x)
式中b=[b0,bl,b2,…,bM],a=[a0,a1,a2,…,aN]分别就是差分方程左、右端的系数向量,x表示输入序列,y表示输出序列。

注意输出序列的长度与输入序列长度相同。

【例2-3】受噪声干扰的信号为x[k]=s[k]+d[k],其中s[k]=(2k)0、9k就是原始信号,d[k]就是噪声。

已知M点滑动平均( Moving Average)系统的输入与输出关系为
y[k] =
试编程实现M点滑动平均系统对受噪声干扰的信号去噪。

解:系统的输入信号x[kl含有有用信号s[k]与噪声信号d[k]。

噪声信号d[k]可以用rand函教产生,将其叠加在有用信号s[k]上,即得到受噪声干扰的输入信号x[k]。

下面的程序实现了对信号x[k]去噪,取M=5。

% program2_3 Signal Smoothing by Moving Average Filter
clear
clc
R=51;
d=rand(1,R)-0、5;
k=0:R-1;
s=2*k、*(0、9、^k);
x=s+d;
figure(1);plot(k,d,'r-、',k,s,'b--',k,x,'g-');
xlabel('Time index k');legend('d[k]','s[k]','x[k]');
M=5;b=ones(M,1)/M;a=1;
y=filter(b,a,x);
figure(2);plot(k,s,'b--',k,y,'g-');
xlabel('Time index k');
legend('s[k]','y[k]');
Time index k
Time index k
图2-3M点滑动平均系统对噪声干扰信号的去噪
程序运行的结果如图3 - 25所示。

图3- 25(a)中3条曲线分别为噪声信号d[k]、有用信号s[k]与受噪声干扰的输入信号x[k]。

图3-25 (b)中。

s[k]为有用信号,y[k]就是经过5
点滑动平均系统去噪的结果。

比较这两条曲线可以瞧出,y[k]与s[k]波形除了有
21-
M
的延迟外,基本上就是相似的,这说明y[k]中的噪声信号被抑制,M点滑动平均系统实现了对受噪声干扰信号的去噪。

4、离散时间系统单位脉冲响应的求解
在MATLAB中,求解离散时间系统单位脉冲响应,可应用信号处理工具箱提供的函数impz,其调用方式为
h= impz(b,a,k)
式中b=[b0,b1,b2,…,bN,]a=[a0,a1,a2,…,aN]分别就是差分方程左、右端的系数向量,k 表示输出序列的取值范围,h就就是系统的单位脉冲响应。

【例2-4】用impz函数求离散时间LTI系统
6y[k]+5y[k-1]+y[k-2]=10x[k]
的单位脉冲响应h[k]。

解:MATLAB程序如下:
%program3_4 离散系统地单位脉冲响应
clear
clc
k=0:10;
a=[6 5 1];
b=[10];
h=impz(b,a,k);
stem(k,h)
图2-4离散系统的单位脉冲响应
5、离散卷积的计算
卷积就是用来计算系统零状态响应的有力工具。

MATLAB信号处理工具箱提供了一个计算两个离散序列卷积与的函数conv,其调用方式为
c= conv(a,b)
式中a,b为待卷积两序列的向量表示,c就是卷积结果。

向量c的长度为向量a,b长度之与减一,即length(c):length(a)+length(b) -1。

【例2-5】已知序列x[k]={1,2,3,4;k=0,1,2,3},y[k]=[1,l,1,1,1;k=0,1,2,3,4},计算x[k]*y[k]并画出卷积结果。

解:MATLAB程序如下:
% program2_5
clear
clc
x=[1,2,3,4];
y=[1,1,1,1,1];
z=conv(x,y); N=length(z); stem(0:N-1,z);
图2-5 离散序列的卷积
Conv 函数也可以用来计算两个多项式的积。

例如多项式）（323++s s 与)23(2
++s s 的
乘积可通过下面的MATLAB 语句求出:
a=[l,0,2,3]; b=[l,3,2]; c= conv(a,b)
语句a=[1,0,2,3]与b=[l,3,2]分别就是多项式）（323++s s 与)23(2
++s s 的向量表
示。

注意,在用向量表示多项式时,应将多项式各项包括零系数项的系数均写入向量的对应元
素中。

如多项式）（323
++s s 中2次方的系数为零,故向量a 的第2个元素也为零。

如果表
示成a=[l,2,3],则计算机将认为表示的多项式为s 2
+2s +3。

上面语句运行的结果为
c=
1 3 4 9 13 6 即
61393)23)(32(2
4
5
2
3
++++=++++s s s s s s s s 三． 实验内容:
1． 使用matlab 计算如下序列()n f 1与()n f 2的卷积与()n f ,绘出她们的时域波形。

()⎪⎪⎩⎪
⎪⎨
⎧==-==其它
,
01,10,21
,11n n n n f ()⎩⎨
⎧≤≤-=其它
,02
2,12n n f
2． 已知某LTI 离散系统,其单位响应()()()4--=n u n u n h ,求该系统在激励为
()()()3--=n u n u n f 时的零状态响应()n y ,并绘出其时域波形。

四． 思考题:
1． 分析实验内容1中序列()n f 1与()n f 2的时域宽度与()n f 的时域宽度的关系。