实验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阶微分方程𝑎3𝑦(3)(𝑡)+𝑎2𝑦(2)(𝑡)+𝑎1𝑦′(𝑡)+𝑎0y(𝑡)=𝑏3𝑦(3)(𝑡)+𝑏2𝑦(2)(𝑡)+𝑏1𝑦′(𝑡)+𝑏0y(𝑡)可用a=[a3,a2,a1,a0]；b=[b3,b2,b1,b0];sys=tf(b,a)获得连续时间LTI模型。注意微分方程中为零的系数一定要写入向量a和b中。【例2-1】描述某力学系统中物体位移y(t)与外力f(t)的关系为m𝑑2𝑦(𝑡)𝑑𝑡2+𝑓𝑑dy(t)dt+𝑘𝑠y(t)=x(t)物体质量m=lkg，弹簧的弹性系数ks=100N/m，物体与地面的摩擦系数fd=2N·s/m，系统的初始储能为零，若外力x(t)是振幅为10、周期为1的正弦信号，求物体的位移y(t)。解：由已知条件，系统的输入信号为x(t)=10sin（2πt），系统的微分方程为𝑑2𝑦(𝑡)𝑑𝑡2+2dy(t)dt+100y(t)=x(t)计算物体位移y(t)的MATLAB程序如下：%program2_1微分方程求解ts=0;te=5;dt=0.01;sys=tf([1],[12100]);t=ts:dt:te;x=10*sin(2*pi*t);y=lsim(sys,x,t);plot(t,y);xlabel('Time(sec)')ylabel('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),系统的微分方程可写成:𝑑2ℎ(𝑡)𝑑𝑡2+2dh(t)dt+100h(t)=10δ(t)物体位移y(t)即系统的冲激响应，计算其的MATLAB程序如下：%program3_2连续时间系统的冲激响应clearclcts=0;te=5;dt=0.01;sys=tf([10],[12100]);t=ts:dt:te;y=impulse(sys,t);plot(t,y);xlabel('Time(sec)')ylabel('h(t)')00.511.522.533.544.55-0.25-0.2-0.15-0.1-0.0500.050.10.150.2Time(sec)y(t)图2-2连续时间系统的冲激响应3.离散的时间系统零状态相应的求解大量的离散时间LTI系统都可以用如下的线性常系数差分方程描述：∑𝑎𝑖𝑦[𝑘−𝑖]=𝑛𝑖=0∑𝑏𝑗𝑥[𝑘−𝑗]𝑚𝑗=0其中a0=1,x[k]、y[k]分别表示系统的输入和输出，n是差分方程的阶数。已知差分方程的n个初始状态和输入x[k]，就可以编程由下式迭代计算出系统的输出：y[k]=-∑aiy[k-i]ni=1+∑𝑏𝑗𝑥[𝑘−𝑗]𝑚𝑗=0在零初始状态下，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点滑动平均(MovingAverage)系统的输入与输出关系为y[k]=1M∑𝑥[𝑘−𝑛]𝑀−1𝑛=0试编程实现M点滑动平均系统对受噪声干扰的信号去噪。解：系统的输入信号x[kl含有有用信号s[k]和噪声信号d[k]。噪声信号d[k]可以用rand函教产生，将其叠加在有用信号s[k]上，即得到受噪声干扰的输入信号x[k]。下面的程序实现了对信号x[k]去噪，取M=5。%program2_3SignalSmoothingbyMovingAverageFilterclearclcR=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('Timeindexk');legend('d[k]','s[k]','x[k]');00.511.522.533.544.55-0.8-0.6-0.4-0.200.20.40.60.81Time(sec)h(t)M=5;b=ones(M,1)/M;a=1;y=filter(b,a,x);figure(2);plot(k,s,'b--',k,y,'g-');xlabel('Timeindexk');legend('s[k]','y[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]分别是差分方程左、右端的系数向量，k05101520253035404550-1012345678Timeindexkd[k]s[k]x[k]05101520253035404550012345678Timeindexks[k]y[k]表示输出序列的取值范围，h就是系统的单位脉冲响应。【例2-4】用impz函数求离散时间LTI系统6y[k]+5y[k-1]+y[k-2]=10x[k]的单位脉冲响应h[k]。解：MATLAB程序如下：%program3_4离散系统地单位脉冲响应clearclck=0:10;a=[651];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_5clearclcx=[1,2,3,4];y=[1,1,1,1,1];z=conv(x,y);012345678910-1.5-1-0.500.511.52N=length(z);stem(0:N-1,z);图2-5离散序列的卷积Conv函数也可以用来计算两个多项式的积。例如多项式）（323ss和)23(2ss的乘积可通过下面的MATLAB语句求出：a=[l，0，2，3]；b=[l，3，2]；c=conv(a，b)语句a=[1，0，2，3]和b=[l，3，2]分别是多项式）（323ss和)23(2ss的向量表示。注意，在用向量表示多项式时，应将多项式各项包括零系数项的系数均写入向量的对应元素中。如多项式）（323ss中2次方的系数为零，故向量a的第2个元素也为零。如果表示成a=[l，2，3]，则计算机将认为表示的多项式为s2+2s+3。上面语句运行的结果为c=1349136即61393)23)(32(24523ssssssss三．实验内容：1．使用matlab计算如下序列和的卷积和，绘出他们的时域波形。01234567012345678910nf1nf2nf其它,01,10,21,11nnnnf其它,022,12nnf2．已知某LTI离散系统，其单位响应4nununh，求该系统在激励为3nununf时的零状态响应ny，并绘出其时域波形。四．思考题：1．分析实验内容1中序列nf1和nf2的时域宽度与nf的时域宽度的关系。