yyx第6章MATLAB在信号与系统中的应用

1第6章MATLAB在信号与系统中的应用2要求掌握MATLAB对信号与系统时频域进行分析的方法。围绕以下内容展开：1、信号表示，表达，运算（加减乘，移位，反折，尺度变换）2、系统时域求解与表达3、系统频域求解与表达4、系统频响描述以上内容又分为连续时间和离散时间3二、信号的表示和可视化，运算•MATLAB中两种方法来表示信号：•向量•符号运算的方法•连续时间信号–向量表示法（f,t两个向量）f(t)=Sa(t)4–符号运算表示法5两个问题：表示，绘图•阶跃信号u(t-t0)•冲激信号•复指数信号•正弦信号•虚指数信号(imag,real,abs,angle)•周期信号连续时间信号的时域MATLAB描述)（0tt6例6.1连续信号的MATLAB描述（1）单位冲激函数（2）单位阶跃函数：（3）复指数函数11111()()0tttxttt其余001)()(1112ttttttutxtutx)j(e)(376.1程序及执行结果8单位冲击信号functioncj(t1,t2,t0)t=t1:dt:t2;n=length(t);x=zeros(1,n);x(1,floor((t1-t0)/dt))=1/dt;stairs(t,x)9周期信号1）、函数square()功能：产生矩形波。•调用格式：x=square(t)；产生周期是2pi，幅值从-1~1的方波。•x=square(t,duty)；产生指定周期的矩形波，其中，duty用于指定脉冲宽度与整个周期的比例。2)、调用格式：x=sawtooth(t)；产生周期是，幅值从-1~1的锯齿波。x=sawtooth(t,width)；用于产生三角波，当width=0.5时，可产生对称的标准三角波；当width=1时，就产生锯齿波。10例•一连续周期矩形信号频率为10Hz，信号幅度在-1~1之间，脉冲宽度与周期的比例是1:2，用128点采样，显示三个周期的信号波形。•Matlab程序如下：f1=10;Um=1;N=128;T=1/f1;nt=3;dt=T/N;t=0:dt:nt*T;xt=Um*square(2*pi*f1*t);plot(t,xt)axis([0,nt*T,-1.1,1.1])11还有一种方法可以表示矩形波•设有一周期方波信号，幅度E=1.5V，周期T=100，脉冲宽度与周期之比为1/2，时间轴上采样点数取1000点.•例如：T=100;f1=1/T;N=1000;t=linspace(0,T,N);dt=T/(N-1);x=1.5*[ones(1,N/2),zeros(1,N/2)];12单位阶跃信号•MATLAB提供一个Heaviside函数，来表示阶跃信号u(t)。f=(t=0)•f=(t=t0);%u(t-t0)•例题：MATLAB如何实现f(t)=u(t+3)-2u(t)的实现？13连续时间信号的运算•相加，相减，相乘，对f的操作（对应时刻进行）•移位，反折，尺度变换(对t进行，并且影响到函数)14二、连续时间系统的时域求解•自由响应，强迫响应•冲激响应，阶跃响应•卷积运算•零输入响应，零状态响应15例6.2LTI系统的零输入响应•n阶线性时不变连续系统的微分方程已知y及其各阶导数的初始值为y(0),y(1)(0)，…，y(n-1)(0)，求系统的零输入响应。•解：方程的解为p1,p2…,pn是方程a1n+a2n-1+…+an+an+1=0的根,C1，…，Cn由y及其各阶导数的初始值来确定。,dddddddddd111121ubtubtubyatyatyatyammmmnnnnntpntptpnCCCtyeee)(212116例6.2LTI系统的零输入响应(续)•C1+C2+…+Cn=y0y0=y(0)•p1C1+p2C2+…+pnCn=Dy0011212111DyCpCpCpnnnnnn0100211121121DD111yyyCCCppppppnnnnnnn17例6.2LTI系统的零输入响应(续)•即V·C=Y0•其解为C=V\Y0•式中•V为范德蒙矩阵，在MATLAB的特殊矩阵库中有vander。调用方法:V=vander(p)]'D,,D,[;]',,,[0100021yyyCCCnnYC1121121111nnnnnppppppV1819例6.3n阶LTI系统的冲激响应•n阶微分方程，写成系统函数为：•冲击响应就是H(s)的拉普拉斯反变换，可以把H(s)展开为极点留数式。•其反变换为11211121)()()(nnnnmmmmasasasabsbsbsbsUsYsHnkkkpsrsH1)(nktpkkrth1e)(206.3程序•a=input('多项式分母系数向量a=(书上取poly([0,-1+2j,-1-2j,-2,-5])）');•b=input('多项式分子系数向量b=(书上取[8,3,1]）');•[r,p]=residue(b,a),%求留数•disp('解析式h(t)=Σr(i)*exp(p(i)*t)')•disp('给出时间数组t=[0:dt:tf]')•dt=input('dt=');%给定时间数组•tf=input('tf=');•t=0:dt:tf;h=zeros(1,length(t));%h的初始化•fori=1:length(a)-1%根数为a的长度减1•h=h+r(i)*exp(p(i)*t);%叠加各根分量•end•plot(t,h),grid21程序运行结果22•信号与系统中如何求解冲激响应？阶跃响应？•MATLAB提供了工具箱函数impulse和step来求解。此时要求用向量方式来表示微分方程。•如：•a=[121];b=[012];impulse(b,a)–该函数的使用方法自己利用帮助系统完善化•step(b,a)思考一下()2()()()2()rtrtrtetet23例6.4卷积的计算•根据卷积公式：•因此编程的过程为：•（1）写出h(t)的MATLAB表达式；•（2）写出u(t)的MATLAB表达式；•（3）利用MATLAB的卷积语句y=conv(u，h)求解•（4）画曲线plot(t,y)。tthuty0d)()()(246.4a显示结果25有关卷积，你还需要知道一些知识•在利用Matlab提供的卷积函数进行卷积运算时，主要是确定卷积结果的时间区间。conv函数默认两信号的时间序列从t=0开始，卷积结果对应的时间序列也从t=0开始。•任意信号的卷积运算。由于conv函数默认信号的时间从0开始，如果有信号不是从0开始，则需要编程用两个向量来确定一个信号，一个是信号波形的非零幅度样值，一个是其对应的时间向量。•通过编写函数实现上述功能。扩展函数为conv_m，其代码如下：•function[y,ny]=conv_m(x1,x2,t1,t2,dt)y=conv(x1,x2);ny1=t1(1)+t2(1);ny2=t1(length(x1))+t2(length(x2));ny=ny1:dt:ny2;0.61()(010)tfeutt26•例已知两信号分别为：求两信号的卷积和并画出卷积波形。0.61()(010)tfeutt2()(015)futt27例题的结果•程序如下：t1=0:0.1:10;f1=exp(-0.6*t1);t2=0:0.1:15;f2=ones(1,length(t2));y=conv(f1,f2);subplot(3,1,1);plot(f1)subplot(3,1,2);plot(f2)subplot(3,1,3);plot(y)28例6.5LTI系统的零状态响应•设二阶连续系统求其冲激响应。若输入为u=3t+cos(0.1t)，求其零状态响应y(t)。解：特征方程2+2+8=0•求出其特征根为p1，p2及相应的留数r1，r2，则冲激响应为：•输出y(t)可用输入u(t)与冲激响应h(t)的卷积求得。uytyty8dd2dd22tptpererth2121)(296.5程序结果30•已知系统的微分方程为：•系统的输入信号为：•求其零状态响应。介绍一下lsim函数功能：对连续系统的响应进行仿真。•调用格式：y=lsim(b,a,x,t)；将输入信号x加在由b,a定义的连续时间系统，将求出的数值解存入变量y中，•如果单独调用函数lsim(b,a,x,t)，则直接在图形区域显示零状态响应的时域仿真波形。MATLAB提供了求解零状态响应的函数lsim22()()2100()()dytdytytftdtdt()10sin(2)ftt31例：运行结果•Matlab程序如下：ts=0;te=5;dt=0.01;t=ts:dt:te;b=[1];a=[1,2,100];f=10*sin(2*pi*t);y=lsim(b,a,f,t);plot(t,y)32例6.6有重极点时的计算•n级放大器，每级的传递函数均为0/(s+0)，求阶跃响应，画出n不同时的波形和频率特性。•解：阶跃信号下系统的输出为求Y(s)的拉普拉斯反变换，即得到阶跃响应y(t)。遇到的困难是重极点，公式复杂，且结果不稳定。为了避开重极点问题，可以有意把极点拉开一些，例如设n个极点散布在-0.950到1.050之间，那样也就可当非重极点来列程序。001()()nnYsss33三、连续时间信号与系统的频域分析(傅里叶)•连续时间信号的频域分析–傅里叶级数–傅里叶变换•连续时间信号的频域分析–傅里叶变换域34例6.7方波分解为多次正弦波之和•图示的周期性方波，其傅里叶级数为•分别计算•直到9次谐波，并做图。图6.7-1输入周期性方波--xtkktttf)12sin(1213sin31sin4)()3sin31(sin4)(sin4)(31tttfttf35程序及显示366.7程序•t=0:.01:2*pi;%设定一个时间数组，有101个点•y=sin(t);plot(t,y),figure(gcf),pause%频率为w=1(f=1/2π)的正弦基波•y=sin(t)+sin(3*t)/3;plot(t,y),pause%叠加三次谐波•%用1，3，5，7，9次谐波叠加•y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;plot(t,y)•%为了绘制三维曲面，要把各次波形数据存为一个三维数组，因此必须重新定义y，重编程。•y=zeros(10,max(size(t)));x=zeros(size(t));•fork=1:2:19•x=x+sin(k*t)/k;y((k+1)/2,:)=x;•end•%将各波形迭合绘出•pause,figure(1),plot(t,y(1:9,:)),grid•line([0,pi+0.5],[pi/4,pi/4])%加上方波幅度线及标注•text(pi+0.5,pi/4,'pi/4')•%将各半波形绘成三维网格图，看出增加谐波阶次对方波逼近程度的影响•halft=ceil(length(t)/2);•pause,figure(2),mesh(t(1:halft),[1:10],y(:,1:halft)),pause%只用正半周波形•clc37例6.11方波的频谱分析•将积分上下限定为0～10s，并将t分成N等份，用求和代替积分。这样，傅立叶变换式可写为•求和可以用f(t)行向量乘以e-jtn列向量来实现。式中Δt是t的增量，在程序中，用dt表示。12j1jjj1