使用MATLAB进行信号与系统分析.ppt

利用利用MATLAB进行信号与系行信号与系统分析分析•MATLAB简介 •信号的信号的MATLABMATLAB表示表示 •利用利用MATLABMATLAB进行系统的时域分析进行系统的时域分析 •利用利用MATLABMATLAB进行信号的频域分析进行信号的频域分析 •利用利用MATLABMATLAB分析系统的频率特性分析系统的频率特性 •利用利用MATLABMATLAB进行连续系统的进行连续系统的S S域分析域分析 •利用利用MATLABMATLAB进行离散系统的进行离散系统的Z Z域分析域分析 •利用利用MATLABMATLAB进行系统的状态变量分析进行系统的状态变量分析MATLAB简介(Matrix Laboratory)•MATLAB的工作方式 •如何获取帮助 •表达式——变量、数值、算数运算符、算数运算符、关系运算符、关系运算符、逻辑运算符、冒号运算符运算符、冒号运算符 •数数组及其运算及其运算 •函数文件函数文件 •循循环（（FOR、、 WHILE 循循环）） •基本基本绘图语句句MATLAB简介•1.1. MATLAB的工作方式的工作方式•（1）窗口命令方式 •（2）运行以.M为扩展名磁盘文件 •2 2．获取帮助．获取帮助 •命令窗口输入： help+函数名函数名 •例如 help plot工作方式举例•%用plot函数画一个方波•t=[-1 0 0 1 1 3]; •f=[0 0 1 1 0 0]; •plot(t,f); •xlabel('t');ylabel('f(t)'); •axis([-1 3 0 2]);直接输入以上命令建一个名字为my_file.M的文件，然后在命令窗口输入文件名及回车。

3．表达式——变量•不需要不需要变量的量的类型型说明明 •变量名的量名的第一个字符必第一个字符必须是字母是字母•变量名量名长度：度：不超不超过31个字符个字符•大写和小写的字母大写和小写的字母视为不同不同的字符的字符 • 例如：num_students = 25•特殊特殊变量：量： • pi 表示圆周率，inf 表示无穷大，NaN(Not a Number)表示不定量，如0/03．表达式——数值• MATLAB用常用常规的十的十进制表示数制表示数值 • 用用i或或j作作为后后缀来表示复数的虚部来表示复数的虚部 • 例1.235e5表示1.235105，x=2+3j •abs(x) 求复数x的模 •angle(x) 求复数x的相角(弧度) •real(x) 求复数x的实部 •imag(x) 求复数x的虚部 •conj(x) 求复数x的共轭3．表达式——运算符号算数运算符算数运算符• + +加 • - - 减 • * * 乘 • /除 • ^乘方 • ‘ 矩阵的复共轭转置3．表达式——运算符号逻辑运算符逻辑运算符• A & B逻辑与(and) • A | B逻辑或(or) • ~A 逻辑非(not)值为0时表示逻辑假(F)，其它任何非零值表示逻辑真。

3．表达式——运算符号关系运算符关系运算符• A < B小于 •A > B 大于 •A <= B 小于等于 •A >= B 大于等于 •A == B 等于 •A ~= B不等于3．表达式——运算符号冒号运算符冒号运算符•表达式 1:10 表示产生一个行向量，它的值为 • 1 2 3 4 5 6 7 8 9 10 •表达式 10:- -2:1 表示产生一个递减的行向量，它的值为 • 10 8 6 4 2 4．数组•（1）数组的构造用冒号用冒号:产生数生数组例 x=2:5 产生一个数组，它的值为 x(1)=2, x(2)=3, x(3)=4, x(4)=5x=linspace(0,2,11)将区间[0，2]均匀抽样11点作为数组x给2维数数组赋值时，用分号分号表示一行的结束，如：z=[1 2; 3 4]用用linspace产生数生数组4．数组•MATLAB 提供了一些产生基本矩阵的函数 • zeros 产生矩阵元素全为0的矩阵 • ones 产生矩阵元素全为1的矩阵 • rand 产生(0,1)均匀分布随机数矩阵 • randn 产生正态分布随机数矩阵4．数组——数组的运算•数数组和一个和一个标量相加或相乘量相加或相乘 •例 y=x-1 z=3*x •2个数个数组的的对应元素相乘除元素相乘除 .* ./ •例 z=x.*y •确定数确定数组大小的函数大小的函数 •size(A) 返回值数组A的行数和列数（二维） •length(B) 确定数组B的元素个数（一维）5.函数文件•M文件的第一行包含function•功能:建立一个函数，可以同MATLAB的库函数一样使用.例：编一个绘制图示波形的函数。

function y=tri(t) y=[ abs(t)<=1].*(1-abs(t));调用函数tri，并画出它的波形t=-2:0.05:2; plot(t,tri(t));6.For 循环•例: 编写计算s=1+2+3+...+100的MATLAB程序 •s=0;for n=1:100 s=s+n; end7.While 循环•s=0;•n=1;•eps=1e-6;while 1/(n*n) > eps s=s+1/(n*n); n=n+1; end例: 计算 的值，且误差小于10-6fprintf('s=%.5f\n',s)8.plot函数--绘图函数(continuous)•t=linspace(0,4*pi,512); •plot(t,sin(t),t,cos(t),'-.'); •title('my figure'); •xlabel('t'); •legend('sin(t)','cos(t)');9 .stem函数--绘图函数(discrete)•k=0:39; •stem(k,cos(0.9*pi*k)); •title('cos(0.9\pik)');Cos(0.9k)波形•基本信号的基本信号的MATLAB表示表示 •(指数信号Aeat ,指数序列ak , 抽样函数Sa(t),正弦型信号, 矩形脉冲信号, 三角脉冲信号 • •信号基本运算的信号基本运算的MATLAB实现 •（尺度变换、翻转、时移、 相加、相乘、差分与求和、微分与积分）信号的信号的MATLABMATLAB表示表示基本信号的MATLAB表示•指数信号Aeaty = A*exp(a*t); •指数序列ak幂运算a.^k实现•正弦型信号内部函数cos( ) 和sin( ) •抽样函数Sa(t)sinc(t) •矩形脉冲信号y = rectpuls(t,width) •三角波脉冲信号y = tripuls(t, width,skew) %decaying exponential•t=0:001:10; •A=1; •a=-0.4; •ft=A*exp(a*t); •plot(t,ft)t=0:0.1:10; A=1; a=-0.4; ft=A*exp(a*t); stem(t,ft)% rectpuls•t=0:0.001:4; •T=1; •ft=rectpuls(t-2*T,T);•plot(t,ft) •axis([0,4,-0.5,1.5])% tripuls•t=-3:0.001:3; •ft=tripuls(t,4,0.5);•plot(t,ft)ft=tripuls(t,4,1);% unit impuls sequence•k=-50:50; •delta=[zeros(1,50),1,zeros(1,50)];•stem(k,delta)function [f,k]=impseq(k0,k1,k2) %产生 f[k]=delta(k-k0)；k1<=k<=k2 k=[k1:k2];f=[(k-k0)==0];k0=0;k1=-50;k2=50; [f,k]=impseq(k0,k1,k2); stem(k,f)% unit step sequence•k=-50:50; •uk=[zeros(1,50), ones(1,51)];•stem(k,uk)function [f,k]=stepseq(k0,k1,k2) %产生 f[k]=u(k-k0);k1<=k<=k2 k=[k1:k2];f=[(k-k0)>=0];k0=0;k1=-50;k2=50; [f,k]=stepseq(k0,k1,k2); stem(k,f)信号基本运算的MATLAB实现•t=-3:0.001:3; •ft1=tripuls(2*t,4,0.5); •subplot(2,1,1) •plot(t,ft1) •title('f(2t)') •ft2=tripuls((2-2*t),4,0.5); •subplot(2,1,2) •plot(t,ft2) •title('f(2-2t)')信号的尺度变换、翻转、时移（平移）信号的尺度变换、翻转、时移（平移）已知三角波f(t)，用MATLAB画出的f(2t)和f(2-2t) 波形信号的相加与相乘信号的相加与相乘•相加用算术运算符“+”实现 •相乘用数组运算符“.*”实现 • •例画信号Aeatcos(w0t+f)的波形 •t=0:0.001:8; •A=1; a=-0.4; •w0=2*pi;phi=0; •ft1=A*exp(a*t).*sin(w0*t+phi); •plot(t,ft1) 离散序列的差分与求和离散序列的差分与求和•差分y=diff(f);•求和y=sum(f(k1:k2));连续信号的微分与积分连续信号的微分与积分•微分 y=diff(f)/h; h为数值计算所取时间间隔•定积分 quad(‘function_name’,a,b);function_name为被积函数名，a和b指定积分区间。

[例]已知三角波f(t)，画出其微分与积分的波形•%differentiation •h=0.001;t=-3:h:3; •y1=diff(f2_2(t))*1/h; •plot(t(1:length(t)-1),y1) •%integration•t=-3:0.1:3; •for x=1:length(t) • y2(x)=quad('f2_2',-3,t(x)); •end•plot(t,y2)function yt=f2_2(t) yt=tripuls(t,4,0.5);三角波f(t)微分与积分的波形利用利用MATLABMATLAB进行系统的时域分析进行系统的时域分析•连续时间系系统零状零状态响响应的求解的求解 •连续时间系系统冲激响冲激响应和和阶跃响响应的求解的求解 •离散离散时间系系统零状零状态响响应的求解的求解 •离散离散时间系系统单位脉冲响位脉冲响应的求解的求解 •离散卷离散卷积的的计算算•t 表示计算系统响应的抽样点向量a=[a3, a2, a1, a0]; b=[b3, b2, b1, b0]; sys=tf(b,a)1.连续时间系统零状态响应的求解y=lsim(sys,f,t)sys=tf(b,a)b和a分别为微分方程右端和左端各项的系数向量f 是系统输入信号向量，sys 是LTI系统模型，借助tf函数获得•连续时间系统冲激响应可用impulse函数直接求出，其调用形式为2.连续时间系统冲激响应和阶跃响应的求解y=impulse(sys,t)连续时间系统阶跃响应可用step函数直接求出，其调用形式为y=step(sys,t)t 表示计算系统响应的抽样点向量 sys 是LTI系统模型3.离散时间系统零状态响应的求解•b,a 分别是差分方程左、右端的系数向量b=[b0,b1,b2,,bM]; a=[a0,a1,a2, ,aN];可用MATLAB表示为y=filter(b,a,f)f 表示输入序列, y 表示输出序列4.离散时间系统单位脉冲响应的求解•b,a 分别是差分方程左、右端的系数向量 •k 表示输出序列的取值范围 •h 就是单位脉冲响应h=impz(b,a,k)5. 离散卷积的计算•例如：(s3+2s+3)(s2+3s+2) •可用下面MATLAB语句求出 •a =[1,0,2,3]; •b =[1,3,2];•c=conv(a,b)c=conv(a,b)式中a,b为待卷积两序列的向量表示，c是卷积结果。

conv函数也可用于计算两个多项式的积例1 求系统y”(t)+2y’(t)+100y(t)=10f(t)的零状态响应，已知f(t)=(sin2t) u(t)•%program3_1微分方程求解•ts=0;te=5;dt=0.01; •sys=tf([1],[1 2 100]); •t=ts:dt:te; •f=10*sin(2*pi*t); •y=lsim(sys,f,t); •plot(t,y); •xlabel('Time(sec)') •ylabel('y(t)')例2 求系统y”(t)+2y’(t)+100y(t)=10f(t)的零状态响应，已知f(t)=d(t)•%program3_2连续时间系统的冲激响应•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)')例3 分析噪声干扰的信号f[k]=s[k]+d[k]通过M点滑动平均系统的响应，其中s[k]=(2k)0.9k是原始信号，d[k]是噪声。

•R =51 ; d = rand(1,R)-0.5; •k=0:R-1; •s=2*k.*(0.9.^k); f=s+d; •figure(1); •plot(k,d,'r-.',k,s,'b--',k,f,'g-'); •M =5; b = ones(M,1)/M; a = 1; •y = filter(b,a,f); •figure(2); •plot(k,s,'b--',k,y,'r-');噪声干扰信号f[k]=s[k]+d[k]通过M点滑动平均系统的响应例4 求系统y[k]+3y[k-1]+2y[k-1]=10f [k]的单位脉冲响应•% program 3_4 离散系统的单位脉冲响应 •k=0:10; •a=[1 3 2]; •b=[1]; •h=impz(b,a,k); •stem(k,h)例5 计算x[k]* y[k]并画出卷积结果，已知x[k]={1,2,3,4; k=0,1,2,3}，y[k]={1,1,1,1,1; k=0,1,2,3,4} •% program 3_5•x=[1,2,3,4]; •y=[1,1,1,1,1]; •z=conv(x,y); •N=length(z); •stem(0:N-1,z);利用MATLAB进行信号的频域分析•周期信号周期信号频谱的的MATLAB实现 •用数用数值积分分析非周期信号分分析非周期信号频谱周期信号频谱的周期信号频谱的MATLABMATLAB实现实现 频谱Cn一般为复数，可分别利用abs和angle函数获得其幅频特性和相频特性。

其调用格式分别为x=abs(Cn) y=angle(Cn)周期信号的频谱Cn 为离散信号，可以用stem画出其频谱图[例1]试用MATLAB画出图示周期三角波信号的频谱解：周期信号的频谱为画三角波信号画三角波信号频谱的的MATLAB程序程序•N=8; •n1=-N:-1; %计算算n=-N到到-1的的Fourier系数系数•c1=-4*j*sin(n1*pi/2)/pi^2./n1.^2; •c0=0; %计算算n=0时的的Fourier系数系数•n2=1:N; %计算算n=1到到N的的Fourier系数系数•c2=-4*j*sin(n2*pi/2)/pi^2./n2.^2; •cn=[c1 c0 c2]; •n=-N:N; •subplot(2,1,1); •stem(n,abs(cn));ylabel('Cn的幅度的幅度'); •subplot(2,1,2); •stem(n,angle(cn)); •ylabel('Cn的相的相位位');xlabel('\omega/\omega0');程序运行程序运行结果果[例2]求周期矩形脉冲的Fourier级数表示式。

并用MATLAB 求出由前N项Fourier级数系数得出的信号近似波形取取A=1, T=2, t t=1, w w0=p p% Gibbs phenomenon•t=-2:0.001:2; •N=input('Number of harmonics= '); •c0=0.5; •fN=c0*ones(1,length(t)); %dc component•for n=0:2:N-1 % even harmonics are zero• fN=fN+cos(pi*n*t)*sinc(n/2);•end •plot(t,fN);N=5N=15N=50N=500用数值积分分析非周期信号频谱用数值积分分析非周期信号频谱数值函数积分quad8可用来计算非周期信号频谱F 是一个字符串，它表示被积函数的文件名a,b 分别表示定积分的下限和上限quad8的返回是用自适应Simpson算法得出的积分值y = quad8('F',a,b) [例]试用数值方法近似计算三角波信号的频谱F(j)= Sa2( / 2)解: 图示三角波可表示为三角波信号频谱的理论值为function y=sf1(t,w);y=(t>=-1 & t<=1).*(1-abs(t)).*exp(-j*w*t);•w=linspace(-6*pi,6*pi,512); •N=length(w);F=zeros(1,N); •for k=1:N • F(k)=quad8('sf1',-1,1,[],[],w(k)); •end•figure(1); •plot(w,real(F));title('') •xlabel('\omega');ylabel('F(j\omega)'); •figure(2); •plot(w,real(F)-sinc(w/2/pi).^2); •xlabel('\omega');title('计算误差');运行运行结果果利用利用MATLABMATLAB进行系统频域分析进行系统频域分析•连续系统频响特性的计算连续系统频响特性的计算 •周期信号通周期信号通过系系统的响的响应 •离散系统频响特性的计算离散系统频响特性的计算连续系统频响特性的计算连续系统频响特性的计算b 分子多项式系数 a 分母多项式系数 w 需计算的H()的抽样点 (数组w中少需包含两个的抽样点)。

计算频响的MATLAB函数H=freqs(b,a,w)例1 三阶归一化的Butterworth低通滤波器的系统函数为•w=linspace(0,5,200); •b=[1];a=[1 2 2 1]; •h=freqs(b,a,w); •subplot(2,1,1); •plot(w,abs(h)); •subplot(2,1,2); •plot(w,angle(h)); 试画出|H()| 和()三三阶Butterworth低通低通滤波器的幅度响波器的幅度响应和相位响和相位响应周期信号通过系统的响应周期信号通过系统的响应[例2]周期方波通过RC系统的响应%p5_2 Periodic signal pass LTI system•T=4;w0=2*pi/T;RC=0.1; •t=-6:0.01:6;N=51; •c0=0.5;xN=c0*ones(1,length(t)); %dc •for n=1:2:N % even harmonics are zero• H=abs(1/(1+j*RC*w0*n)); • phi=angle(1/(1+j*RC*w0*n)); • xN=xN+H*cos(w0*n*t+phi)*sinc(n*0.5);•end •plot(t,xN); •xlabel(['time RC=',num2str(RC)]);grid; •set(gca,'xtick',[-5 -3 -1 0 1 3 5]);周期方波通周期方波通过RC系系统的响的响应离散系统频响特性的计算离散系统频响特性的计算计算频响的MATLAB函数b 分子的系数 a 分母系数w 抽样的频率点(至少2点)， w在0~2之间幅频特性: abs, 相频特性: angleh = freqz(b,a,w)b=[1];a1=[1 -0.9]; a2=[1 0.9];w=linspace(0,2*pi,512);h1=freqz(b,a1,w);h2=freqz(b,a2,w);plot(w/pi,abs(h1),w/pi,abs(h2),':');legend('\alpha=0.9','\alpha=-0.9');利用利用MATLABMATLAB进行进行连续系统的复频域分析连续系统的复频域分析•部分分式展开的MATLAB实现 •H(s)的零极点与系的零极点与系统特性的特性的MATLAB计算一、部分分式展开的MATLAB实现[r,p,k]=residue(num,den)num,den分别为F(s)分子多项式和分母多项式的系数向量。

r为部分分式的系数，p为极点，k为多项式的系数 若为真分式，则k为零二、H(s)的零极点与系统特性的MATLAB计算计算多项式根roots的函数可用于计算H(s)的零极点r=roots(N) %计算多项式N的根H(s)零极点分布图可用pzmap函数画出，调用形式为pzmap(sys)表示画出sys所描述系统的零极点图[例1]用部分分式展开法求F(s)的反变换•%program7_1•format rat %将结果数据以分数的形式输出•num=[1 2]; •den=[1 4 3 0]; •[r,p]=residue(num,den) 运行结果为r =-1/6 ,-1/2 ,2/3 p = -3 ,-1 ,0故F(s)可展开为[例2]用部分分式展开法求F(s)的反变换•% program7_2•num=[2 3 0 5]; •den=conv([1 1],[1 1 2]); %将因子相乘的形式转换成多项式的形式•[r,p,k]=residue(num,den) •magr=abs(r) %求求r的模的模•angr=angle(r) %求求r的相角的相角运行结果为•r =-2.0000 + 1.1339i, -2.0000 - 1.1339i, 3.0000 •p =-0.5000 + 1.3229i, -0.5000 - 1.3229i, -1.0000 •k =2 •magr =2.299, 2.2991, 3.0000 •angr =2.6258, -2.6258, 0 故F(s)可展开为[例3]试画出系统 的零极点分布图，求其单位冲激响应h(t)和频率响应H(j)，并判断系统是否稳定。

•num=[1];den=[1 2 2 1]; •sys=tf(num,den); •poles=roots(den) •figure(1);pzmap(sys); •t=0:0.02:10; •h=impulse(num,den,t); •figure(2);plot(t,h) •title('Impulse Respone') •[H,w]=freqs(num,den); •figure(3);plot(w,abs(H)) •xlabel('\omega') •title('Magnitude Respone')运行运行结果果利用利用MATLABMATLAB进行离散系统的进行离散系统的Z Z域分析域分析•部分分式展开的部分分式展开的MATLABMATLAB实现 •H(z)的零极点与系的零极点与系统特性的特性的MATLABMATLAB计算算一、部分分式展开的MATLAB实现[r,p,k]=residuez(num,den)num,den分别为F(z)分子多项式和分母多项式的系数向量r为部分分式的系数，p为极点，k为多项式的系数 若为真分式，则k为零二、H(z)的零极点与系统特性的MATLAB计算利用tf2zp函数计算H(z)的零极点，调用形式为[z,p,k]=tf2zp(b,a)H(z)零极点分布图可用zplane函数画出，调用形式为zplane(b,a)b和a分别为H(z)分子多项式和分母多项式的系数向量。

返回值z为零点、p为极点、 k为增益常数[例1] 将F(z)用部分分式展开•%program8_1•num = [18]; •den = [18 3 -4 -1]; •[r,p,k] = residuez(num,den)运行结果为r =0.3600 , 0.2400 , 0.4000 p =0.5000 , -0.3333 , -0.3333 k =[]故F(z)可展开为[例2]试画出系统 的零极点分布图，求其单位冲激响应h[k]和频率响应H(ejW)•% program 8_2 •b =[1 2 1];a =[1 -0.5 -0.005 0.3]; •figure(1);zplane(b,a); •num=[0 1 2 1]; •den=[1 -0.5 -0.005 0.3]; •h=impz(num,den); •figure(2);stem(h) •xlabel('k') •title('Impulse Respone') •[H,w]=freqz(num,den); •figure(3);plot(w/pi,abs(H)) •xlabel('Frequency \omega') •title('Magnitude Respone')运行运行结果果利用利用MATLABMATLAB进行系统状态变量分析进行系统状态变量分析•微分方程到状微分方程到状态方程的方程的转换 •状状态方程系方程系统函数矩函数矩阵H(s)的的计算算 •MATLAB求解求解连续系系统的状的状态方程方程 •MATLAB求解离散系求解离散系统的状的状态方程方程微分方程到状态方程的转换微分方程到状态方程的转换•num,den 分别表示系统函数H(s)的分子和分母多项式 •A,B,C,D 分别为状态方程的矩阵。

[A,B,C,D]=tf2ss(num,den)状态方程系统函数矩阵状态方程系统函数矩阵H H( (s s) )的计算的计算•A,B,C,D 分别表示状态方程的矩阵 •K 表示函数ss2tf计算的与第k个输入相关的系统函数，即H(s)的第k列 •num 表示H(s)第k列的m个元素的分子多项式 •den 表示H(s)公共的分母多项式[num,den]= ss2tf (A,B,C,D,k)MATLABMATLAB求解连续系统的状态方程求解连续系统的状态方程•sys由函数ss构造的状态方程模型•t需计算的输出样本点, t=0:dt:Tfinal •f(:,k) 系统第k个输入在t上的抽样值•x0系统的初始状态(可缺省) •y(:,k) 系统的第k个输出 •to实际计算时所用的样本点； x系统的状态 获得连续系统状态方程的计算机表示模型sys = ss(A,B,C,D)求解状态方程[y,to,x]=lsim(sys,f,t,x0);MATLABMATLAB求解离散系统的状态方程求解离散系统的状态方程•sys 由函数ss构造的状态方程模型•f(:,k)系统第k个输入序列•x0系统的初始状态(可缺省) •y(:,k)系统的第k个输出•n序列的下标;x系统的状态获得离散状态方程的计算机表示模型sys = ss(A,B,C,D,[])求解状态方程[y,n,x]=lsim(sys,f,[],x0)[例1]写出系统 的状态方程•由[A,B,C,D]=tf2ss([1],[1 5 10])•可得所以系统的状态方程为C=0 D=0[例2]已知某连续系统的状态方程和输出方程为其初始状态和输入分别为求该系统的系统函数矩阵H(s) 和输出。

计算系统函数矩阵计算系统函数矩阵H(s)•A=[2 3;0 -1];B=[0 1; 1 0]; •C=[1 1; 0 -1];D=[1 0; 1 0]; •[num1,den1]=ss2tf(A,B,C,D,1)•[num2,den2]=ss2tf(A,B,C,D,2) •运行结果 •num1 =1 0 -1 • 1 -2 0 •den1 =1 -1 -2 •num2 =0 1 1 • 0 0 0 •den2 =1 -1 -2 •所以系统函数矩阵H(s)为计算输出计算输出•%Program 9_1 •A=[2 3;0 -1];B=[0 1; 1 0]; •C=[1 1; 0 -1];D=[1 0; 1 0]; •x0=[2 -1]; •dt=0.01;t=0:dt:2; •f(:,1)=ones(length(t),1); •f(:,2)=exp(-3*t)'; •sys=ss(A,B,C,D); •y=lsim(sys,f,t,x0); •subplot(2,1,1);plot(t,y(:,1),'r'); •ylabel('y1(t)');xlabel('t'); •subplot(2,1,2);plot(t,y(:,2)); •ylabel('y2(t)');xlabel('t');运行结果运行结果[例3]已知某离散系统的状态方程和输出方程为其初始状态和输入分别为求该系统的输出。

计算输出计算输出•%Program 9_2 •A=[0 1;-2 3];B=[0;1]; •C=[1 1; 2 -1];D=zeros(2,1); •x0=[1;-1];•N=10;f=ones(1,N); •sys=ss(A,B,C,D,[]); •y=lsim(sys,f,[],x0); •subplot(2,1,1);y1=y(:,1)'; •stem((0:N-1),y1); •xlabel('k');ylabel('y1'); •subplot(2,1,2);y2=y(:,2)'; •stem((0:N-1),y2); •xlabel('k');ylabel('y2');运行结果运行结果。