第10章 MATLAB在信号与系统分析中的应用 10.0 引言 10.1 MATLAB基础 10.2 信号的MATLAB表示 10.3 用MATLAB实现系统的时域分析 10.4 用MATLAB实现连续系统的频域分析 10.5 用MATLAB实现连续系统的S域分析 10.6 用MATLAB实现离散系统的Z域分析 10.7 MATLAB在系统状态空间分析中的应用

10.0　引　　言 　　　一般来说,MATLAB系统包括下面五个主要部分。 　　(1)编程语言:是一种以矩阵和数组为基本单位的编程语言; 　　(2)工作环境:包括了一系列应用工具,提供编程和调试程序的环境; 　　(3)图形处理:包括绘制二维、三维图形和创建图形用户接口; 　　(4)数学库函数:包含了大量的数学函数,也包括复杂的功能; 　　(5)应用程序接口:提供接口程序,可使MATLAB与其它语言程序进行交互。

10.1　MATLAB基础 10.1.1　MATLAB语言的特点 　　MATLAB语言具有以下特点： 　　(1)编程效率高。 　　MATLAB编程语言作为面向科学与工程计算的高级语言,允许用数学形式的语言编写程序,且比Basic、Fortran和C等语言更加接近我们书写计算公式的思维方式。用MATLAB编写程序犹如在演算纸上排列出公式与求解问题,因此,MATLAB语言也可通俗地称为演算纸式科学算法语言,它编写简单,编程效率高,易学易懂。

　　(2)用户使用方便。 　　MATLAB语言是一种解释执行的语言(在没被专门的工具编译之前),它灵活、方便,其调试程序手段丰富。MATLAB运行时,在命令行每输入一条MATLAB语句(命令),包括调用M文件的语句,计算机就立即对其进行处理,完成编译、连接和运行的全过程。在运行m文件时,如果有错,计算机屏幕提示出错信息,经用户修改后再执行,直到正确为止。

　　(3)扩充能力强。 　　高版本的MATLAB语言有丰富的库函数,在进行复杂的数学运算时可以直接调用。用户可以根据需要建立和扩充新的库函数,以提高MATLAB的使用效率,扩充其功能。

　　(4)语句简单,内涵丰富。 　　MATLAB语言中最基本、最重要的成分是函数,其一般形式为 　　［a,b,c,…］=fun(d,e,f,…) 　　即一个函数由函数名,输入变量d,e,f,…和输出变量a,b,c…组成。同一函数名F,可以有不同数目的输入变量(包括无输入变量)及不同数目的输出变量,代表着不同的含义。这不仅使MATLAB的库函数功能更丰富,而且大大减少了需要的磁盘空间,使得MATLAB编写的m文件简单、短小而高效。

　　(5)高效方便的矩阵和数组运算。 　　MATLAB语言像Basic、Fortran和C语言一样规定了矩阵的算术运算符、关系运算符、逻辑运算符、条件运算符及赋值运算符,而且这些运算符大部分可以毫无改变地运用到数组间的运算中,有些运算符(如算术运算符)只要增加“·”就可用于数组间的运算。

　　(6)方便的绘图功能。 　　MATLAB有一系列绘图函数(命令),调用不同的绘图函数可方便地绘制线性坐标、对数坐标、半对数坐标及极坐标,通过相应的命令还可以在图上标出图题、XY轴标注、格(栅)等。 　　总之,MATLAB语言的设计思想可以说代表了当前计算机高级语言的发展方向,读者在不断使用中会发现其具有巨大的潜力。

10.1.2　MATLAB工作环境简介 　　1.启动MATLAB 　　有三种方法启动MATLAB: 　　(1)双击Windows桌面上的MATLAB快捷图标; 　　(2)通过“开始”菜单的“程序”子菜单中的MATLAB项启动; 　　(3)在MATLAB目录中搜索到可执行程序MATLAB.exe,双击该程序使之启动。 　　启动后,MATLAB主界面如图10.1－1所示。

图10.1－1　MATLAB主界面

　　MATLAB主界面大致包括以下几个部分: 　　(1)菜单项; 　　(2)工具栏; 　　(3)“CommandWindow”(命令窗口),在提示符后直接输入命令可以执行相关的命令; 　　(4)“LaunchPad”(分类帮助文件夹); 　　(5)“Workspace”(工作空间)，该程序窗口中列出了程序计算过程中产生的变量及其对应的数据的尺寸、字节和类型。选中一个变量,单击鼠标右键则可根据菜单进行相应的操作。

　　(6)“CommandHistory”(命令的历史记录)窗口，该窗口记录用户每次开启MATLAB的时间,以及每次开启MATLAB后在MATLAB命令窗口中运行过的所有命令行。这些命令行记录可以被复制到命令窗口中再运行。选中该窗口中的任一命令记录,然后单击鼠标右键,则可根据弹出的菜单进行相应的操作。 　　(7)“CurrentDirectory”窗口，其中包含当前目录选项。

　　2.程序编辑器 　　1)命令文件 　　命令文件没有输入参数,也不返回输出参数,只是一些命令行的组合。命令文件中的语句可以访问MATLAB工作空间(Workspace)中的所有数据,在运行的过程中所产生的变量均是全局变量。这些变量一旦生成,就一直保存在内存空间中,除非用户将它们清除(用clear命令)。运行一个命令文件等价于从命令窗口中按顺序连续执行文件中的命令。由于命令文件只是一串命令的结合,因此程序不需要预先定义,而只需按命令窗口中的命令输入顺序,依次将命令编辑在命令文件中即可。如果某个命令不需要显示结果,则在该命令后加上“;”。注意文件名一定是“.m”。命令文件的建立过程如下:

　　(1)进入程序编辑器(MATLABEditor／Debug):从“File”菜单中选择“New”及“mfile”或单击“Newmfile”按钮; 　　(2)输入程序:在“MATLABEditor/Debug”窗口输入MATLAB程序; 　　(3)保存程序:单击“Save”按钮,出现一个对话框,在文件名框中键入一个文件名,单击“保存”按钮，一个m文件便保存在磁盘上了。

　　　运行命令文件时,该m文件中的命令可以访问MATLAB工作区中的所有变量,而且其中的所有变量也成为工作区的一部分。命令文件运行结束,所产生的变量保留在工作区,直至关闭MATLAB或用命令删除。下面是一个命令文件的例子,程序如下: 　　　%文件名example.m 　　　x=1;y=2;z=3 　　　items=x+y+z 　　　cost=x*5+y*2+z*9 　　　averagecost=cost/items

当这个文件在程序编辑窗口输入并以名为example 　　当这个文件在程序编辑窗口输入并以名为example.m的m文件存盘后,只需在MATLAB命令编辑窗口键入example即可运行,并显示同命令窗口输入命令一样的结果。在m文件中,程序的注释是以符号“％”开始直到该行结束的部分,程序执行时会自动忽略注释语句。

　　2)函数文件 　　如果m文件的第一行包含function,则该文件就是函数文件。每个函数文件都定义一个函数。能够像调用库函数一样方便地调用函数文件,从而可扩展MATLAB的功能。如果对于一类特定的问题,建立起许多函数m文件,就能形成工具箱。 　　从形式上看,函数文件与命令文件的区别在于:命令文件的变量在文件执行完后保留在内存中;而函数文件内定义的变量仅在函数文件内部起作用,当函数文件执行结束后,这些内部变量将被清除。

　　函数m文件的第一行有特殊的要求,其形式必须为 　　function［输出变量列表］=函数名(输入变量列表) 　　函数体语句; 　　例:functiony=f(x) 　　 y=sin(x); 　　注意:函数名应该和m文件名相同。

10.1.3　学习MATLAB的基本方法 　　1.help命令 　　在命令窗口中使行help命令能够获得范围不同的帮助信息,例如: 　　(1)运行helphelp,将得到如何使用help的提示; 　　(2)直接运行help,会列出可以用于help显示的所有主题(topic); 　　(3)运行help(topic),可获得有关该主题的帮助,比如,想对二维图形(graph2d)编程有所了解,可运行helpgraph2d。 　　(4)对每个主题(topic)中的任何命令的用法,同样可以用help来查看。如对于二维图形(graph2d)命令中的plot,用help查看的方法是:helpplot。

　　2.lookfor命令 　　当用户要查找具有某种功能的命令但不知道其准确名字时,help就无能为力了。而lookfor可以根据用户提供的完整或不完整的关键词,去搜索出一组与之有关的命令,用户可从列表中挑选出满足需要的命令。 　　如利用lookfor命令查找矩阵求逆函数: 　　　>>lookforinverse

　　3.doc、helpwin和helpdesk命令 　　在命令窗口中运行doc、helpwin和helpdesk命令中的任何一个,都会打开一个名为“help”的帮助窗口。 　　4.demo命令 　　demo命令用于查看集成于MATLAB环境内的各种演示内容。在MATLAB的命令窗口键入demo命令可以得到演示界面,从而可以方便用户了解MATLAB的基本功能。 　　5.帮助菜单 　　启动MATLAB应用程序后,单击“help”主菜单,则会弹出一系列子菜单,可以根据菜单直接进行操作。

10.2　信号的MATLAB表示 10.2.1　连续信号的MATLAB表示 　　严格地说,MATLAB不能处理连续信号,它是用连续信号在等间隔点的样值来近似表示连续信号的。当采样间隔足够小时,这些样值就能较好地近似表示连续信号。 　　MATLAB提供了大量的生成基本信号的函数。最常用的指数信号、正弦信号是MATLAB的内部函数,即不安装任何工具箱就可调用的函数。

1.单位阶跃信号 单位阶跃信号的数学模型: 在t=t1处跃升的阶跃信号可写为ε(t-t1),定义为

　单位阶跃信号m文件如下： 　%单位阶跃信号m文件 　%信号从t0到tf,在t1(t0≤t1≤tf)前为0,到t1处有一跃变,t1以后为1 　clear;t0=0;tf=5;dt=0.1;t1=input(′t1=′); 　t=［t0:dt:tf］;　　　　　　 　%时间序列 　kt=length(t);[DW]%总的时间点数 k1=floor((t1-t0)/dt);[DW]%求t1对应的样本序号 　x2=［zeros(1,k1),ones(1,kt-k1)］; %产生阶跃信号 　subplot(2,2,3),stairs(t,x2),gridon %绘图 　axis(［0,5,0,1.1］)　%为了使方波顶部避开图框,改变图框坐标

图10.2－1　单位阶跃信号

2.复指数信号 　　复指数信号的数学模型: 　　 x3(t)=e(u+jω)t 　　若ω=0,它是实指数函数;若u=0,则为虚指数函数,其实部为余弦函数,虚部为正弦函数。

复指数信号[WTBZ]m文件如下： %复指数信号m文件 %信号从t0到tf clear; t0=0;tf=6;dt=0 　　复指数信号[WTBZ]m文件如下： 　　　%复指数信号m文件 　　　%信号从t0到tf 　　　clear; 　　　t0=0;tf=6;dt=0.05; 　　　t=［t0:dt:tf］;　　　　　　　　　 %时间序列 　　　alpha=-0.5;w=10; 　　　x3=exp((alpha+j*w)*t); %产生复指数信号 　　　subplot(2,1,1),plot(t,real(x3)),gridon %绘图 　　　subplot(2,1,2),plot(t,imag(x3)),gridon %绘图

图10.2－2　复指数信号实部

图10.2－3　复指数信号虚部

3. 矩形脉冲 调用MATLAB的内部函数可产生一矩形脉冲信号。在MATLAB中用rectpuls函数表示矩形脉冲信号,即 3.矩形脉冲 　　调用MATLAB的内部函数可产生一矩形脉冲信号。在MATLAB中用rectpuls函数表示矩形脉冲信号,即 　　　　　 y=rectpuls(t,width) 　　用以产生一个幅度为1,宽度为width,以t=0为对称轴的矩形波。Width的默认值为1。例如产生一个以t=2T为对称中心的矩形脉冲信号的MATLAB程序如下,取T=1: 　　　%矩形脉冲信号m文件 　　　t=0:0.001:4;T=1; 　　　ft=rectpuls(t-2*T,2*T); 　　　plot(t,ft)

图10.2－4　矩形脉冲

图10.2－5　三角波

10.2.2　离散信号的MATLAB表示 　　由于MATLAB数值计算的特点,故用它来表示离散信号是非常方便的。在MATLAB中,用一个列向量来表示一个有限长序列,而这样的向量并没有包含采样位置的信息,要完全表示一个序列x(k),需用x和k两个向量,例如: 　 　x(k)={2,1,-1,3,2,4,6,1} 下面的箭头指示的是k=0时的采样点。该序列在MATLAB中表示为 　　k=［-3,-2,-1,0,1,2,3,4］;　x=［2,1,-1,3,2,4,6,1］ 　　若不需要采样位置信息或这个信息是多余的(例如该序列从k=0开始),可以只用x向量来表示。计算机的内存有限,[WTBZ]MATLAB无法表示无限长序列。 ↑

　　1.单位脉冲序列 　　单位脉冲序列的表达式: k=0 k=其余 延迟ks的单位脉冲序列表达式: k=k s k=其余

本例取ks=3,此单位脉冲序列的m文件如下: 　　　％单位脉冲序列m文件 　　　clear,k0=0;kf=10;ks=3; 　　　k1=k0:kf;x1=［zeros(1,ks-k0),1,zeros(1,kf-ks)］;[WB]%单位脉冲序列的产生 　　　subplot(2,2,1),stem(k1,x1,′.′);title(′单位脉冲序列′)%绘图 所绘制的单位脉冲序列如图10.2－6所示。

图10.2－6　单位脉冲序列

2.单位阶跃序列 单位阶跃序列的表达式: 延迟ks的单位阶跃序列表达式:

本例取ks=3。 ％单位阶跃序列m文件 clear,k0=0;kf=10;ks=3; k2=k0:kf;x2=[zeros(1,ks-k0),ones(1,kf-ks+1)]; %单位阶跃序列的产生 subplot(2,2,3),stem(k2,x2,'.');title('单位阶跃序列') %绘图

图10.2－7　单位阶跃序列

3.复指数序列 　　复指数序列的表达式: 　　若ω=0,它是实指数序列;若α=0,则为虚指数序列,其实部为余弦序列,虚部为正弦序列。本例取α=-0.2,ω=0.5,该复指数序列的实部和虚部如图10.2－8所示，其m文件如下：

％复指数序列m文件 clear,k0=0;kf=20;ks=3; k3=k0:kf;x3=exp((-0. 2+0. 5j) ％复指数序列m文件 　　　clear,k0=0;kf=20;ks=3; 　　　k3=k0:kf;x3=exp((-0.2+0.5j)*k3);%复指数序列的产生 　　　subplot(1,2,1),stem(k3,real(x3),′.′);line(［0,10］,［0,0］)%绘图 　　　xlabel(′实部′) 　　　subplot(1,2,2),stem(k3,imag(x3),′.′);line(［0,10］,［0,0］)%绘图 　　　xlabel(′虚部′)

图10.2－8　复指数序列的实部、虚部

10.2.3　信号基本运算的MATLAB实现 　　 1.信号的尺度变换、翻转、平移(时移) 　　 信号的尺度变换、翻转、平移运算,实际上是函数自变量的运算。在信号的尺度变换(f(at)和f［ak］)中,函数的自变量乘以一个常数,在MATLAB中可用算术运算符“*”来实现。在信号翻转运算(f(-t)和f［-k］)中,函数的自变量乘以一个负号,在MATLAB中可以直接写出。在信号平移运算(f(t±t0)和f［k±k0］)中,函数自变量加、减一个常数,在MATLAB中可用算术运算符“+”或“-”来实现。

　　例10.2－1　三角波f(t)如图10.2－9(a)所示,试利用MATLAB画出f(2t)和f(2-2t)的波形。 　　　　解　实现f(2t)和f(2-2t)的程序如下: 　　　　　　　％program10.2－1 　　　　　　　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)′)

图10.2－9　例10.2－1图

　　2.连续信号的微分与积分 　　连续信号的微分可用diff近似计算,例如y=sin′(x2)=2xcos(x2)可由以下MATLAB语句近似实现: 　　　h=.001;x=0:h:pi; 　　　y=diff(sin(x.^2))/h; 　　连续信号的定积分可由MATLAB中的quad函数或quad8函数实现,其调用格式为 　　　quad(′function_name′,a,b) 其中,function_name为被积函数名,a和b指定积分区间。

例10. 2－2 三角波f(t)如图10. 2－10(a)所示,试利用MATLAB画出 和 　　例10.2－2　三角波f(t)如图10.2－10(a)所示,试利用MATLAB画出　　　和 的波形。 　　解　为便于利用quad函数计算信号的积分,将三角波f(t)写成MATLAB函数,函数名为f2_tri,程序如下: 　　　％program10.2－2 　　　functionyt=f2_tri(t) 　　　yt=tripuls(t,4,0.5); 　　利用diff和quad函数,并调用f2_tri可实现三角波信号f(t)的微分、积分,程序如下:

％program微分 h=0. 001;t=-3:h:3; y1=diff(f2_tri(t)) 　　　　　％program微分 　　　　　h=0.001;t=-3:h:3; 　　　　　y1=diff(f2_tri(t))*1/h; 　　　　　plot(t(1:length(t)-1),y1) 　　　　　title(′df(t)/dt′) 　　　　　%program积分 　　　　　t=-3:0.1:3; 　　　　　forx=1:length(t) 　　　　　　 y2(x)=quad(′f2_tri′,-3,t(x)); 　　　　　end 　　　　　plot(t,y2)(a) 　　　　　title(′integraloff(t)′) 　　微分、积分波形分别如图10.2－10(b)、(c)所示。

图10.2-10 例10.2-2图

3.离散序列的差分与求和 离散序列的差分 ,在MATLAB中用diff函数实现,其调用格式为 y=diff(f) 　　离散序列的求和 与信号相加运算不同,求和运算是把k1和k2之间的所有样本f［k］加起来,在MATLAB中用sum函数实现,其调用格式为 　　　y=sum(f(k1:k2))

例10.2－3 用MATLAB计算指数信号(-1.6)kε［k］的能量。 解 离散信号的能量定义为 　　％program10.2－3 　　k=0:10;A=1;a=-1.6; 　　fk=A*a.^k; 　　W=sum(abs(fk).^2) 运行结果为 　　W=1.9838e+004

　10.3　用MATLAB实现系统的时域分析 10.3.1　连续系统冲激响应的求解 　　方法1　应用微分方程。 　　MATLAB提供了专门用于求LTI系统的冲激响应和阶跃响应的函数。系统微分方程为 impulse(b,a)用于绘制向量a和b定义的LTI系统的冲激响应,step(b,a)用于绘制向量a和b定义的LTI系统的阶跃响应。其中,a和b表示系统方程中ai、bi组成的向量。

　　例10.3－1　求以下系统的冲激响应和阶跃响应: 解　程序如下： 　　　　　％program10.3－1 　　　　　a=［746］; 　　　　　b=［11］; 　　　　　subplot(2,1,1) 　　　　　impulse(b,a) 　　　　　subplot(2,1,2) 　　　　　step(b,a)

图10.3－1　例10.3－1图

方法2 应用系统函数。 例10.3－2 求n阶LTI系统的冲激响应。 解 设系统函数为 其特性可用系统函数分子、分母系数向量b和a来表示。

　　对于物理可实现系统，n≥m，即length(a)≥length(b)。length(a)－1就是系统的阶次。冲激函数的拉普拉斯变换为F(s)＝1，则系统对冲激函数的响应的拉普拉斯变换为Y(s)=H(s)F(s)=H(s)，冲激响应就是H(s)的拉普拉斯逆变换，可把H(s)展开为部分分式。如果H(s)的分母多项式没有重根，则 故有冲激响应

程序如下： ％program10.3－2 　　　a=input(′多项式分母系数向量a=′); 　　　b=input(′多项式分子系数向量b=′); 　　　［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=′); %输入dt及终点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); 　　　fori=1:length(a)-1 %根数为a的长度减1　　　　　 h=h+r(i)*exp(p(i)*t); %叠加各根分量 　　　end 　　　plot(t,h) 　　　grid 　　运行结果(用通用程序求一个五阶系统的冲激响应,按提示输入分子、分母系数向量和时间数组): 　　　a=poly(［0,-1+2i,-1-2i,-2,-5］); 　　　b=［8,3,1］; 　　　t=0:0.2:8;

　　因为题中给出的分母是系统的极点,而不是多项式系数,故要求出其系数可用poly函数,其格式为 　　　a=poly(p)(其中p为极点向量) 　　即可求出h,画出的曲线如图10.3－2所示。

图10.3－2　例10.3－2图(5阶LTI系统的冲激响应)

10.3.2　连续系统零状态响应的求解 　　方法1　应用MATLAB工具箱提供的函数lsim。 　　LTI连续系统以常系数微分方程描述,系统的零状态响应可通过求解初始状态为零的微分方程得到。 　　在MATLAB中,控制系统工具箱提供了一个用于求解零初始条件微分方程数值解的函数lsim,其调用格式为 　　y=lsim(sys,f,t) 式中,t表示计算系统响应的抽样点向量;f是系统输入信号向量;sys是LTI系统模型,用来表示微分方程、差分方程、状态方程。

在求解微分方程时,微分方程的LTI系统模型sys要借助tf函数获得,其调用方式为 sys=tf(b,a) 式中,b和a分别为微分方程右端和左端各项的系数向量。例如对三阶微分方程 可用 a=[a3,a2,a1,a0];b=[b3,b2,b1,b0];sys=tf(b,a) 获得LTI模型。微分方程中系数为零也要写入向量a和b中。

例10.3－3　系统的微分方程为 求系统在输入为 时的零状态响应。 解　MATLAB程序如下: 　　　％program10.3－3 　　　ts=0;te=5;dt=0.01; 　　　sys=tf(［1］,［1277］); 　　　t=ts:dt:te; 　　　f=10*sin(2*pi*t); 　　　y=lsim(sys,f,t); 　　　plot(t,y); 　　　xlabel(′Time(sec)′) 　　　ylabel(′y(t)′)

图10.3－3　例10.3－3程序运行结果

方法2 应用公式yf(t)=h(t)*f(t)。 例10.3－4 设二阶连续系统的微分方程为 求系统的冲激响应。若输入为f(t)=3t+cos(0.1t),求系统的零状态响应。

　　解　求冲激响应,系统微分方程的特征方程为 其特征根为p1、p2,相应的系数为r1、r2,则冲激响应为 输出y(t)为输入f(t)与冲激响应h(t)的卷积。

程序如下： 　　　％program10.3－4 　　　clf,clear 　　　a=input(′多项式分母系数向量a=　′); 　　　b=input(′多项式分子系数向量b=　′); 　　　t=input(′输入时间序列t=［0:dt:tf］′); 　　　f=input(′输入序列f=　′); 　　　%a=［1,2,8］;b=1;t=［0:0.1:5］; 　　　f=3*t+cos(0.1*t); 　　　tf=t(end); 　　　dt=tf/(length(t)-1);

%用极点留数法求冲激响应 ［r,p,k］=residue(b,a); h=r(1). exp(p(1). t)+r(2). exp(p(2) %用极点留数法求冲激响应 　　　［r,p,k］=residue(b,a); 　　　h=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t); 　　　%画出冲激响应h(t) 　　　subplot(2,1,1),plot(t,h);grid 　　　%求u和h的卷积,得输出y(t) 　　　y=conv(f,h)*dt; 　　　%画出输出y(t) 　　　subplot(2,1,2), 　　　plot(t,y(1:length(t)));grid　　%画出输出y(t) 　　运行结果(执行这个程序时,取a=［1,2,8］,b=1,t=［0:0.1:5］以及输入为f=3*t+cos(0.1*t))如图10.3－4所示。

图10.3－4　例10.3－4程序运行结果

10.3.3 离散系统零状态响应的求解 LTI离散系统可用线性常系数差分方程描述: 其中f［k］、y［k］分别表示系统的输入和输出,n是差分方程的阶数。已知差分方程的n个初始状态和输入f［k］,就可以编程迭代计算出系统的输出:

在零初始状态下,MATLAB工具箱提供了一个filter函数,计算由差分方程描述的系统响应,其调用格式为 　　　在零初始状态下,MATLAB工具箱提供了一个filter函数,计算由差分方程描述的系统响应,其调用格式为 y=filter(b,a,f) 式中,b=［b0,b1,b2,…,bm］,a=［a0,a1,a2,…,an］分别是差分方程左右端的系数向量,f表示输入序列,y表示输出序列。注意输出序列和输入序列的长度相同。

例10.3－5 系统的输入输出关系为 　　输入信号为f［k］=s［k］+d［k］,其中s［k］=(2k)0.9k,d［k］是随机信号。试用MATLAB编程求系统的零状态响应。 　　解　随机信号d［k］可由[WTBZ]MATLAB提供的rand函数产生,将其叠加在s［k］上可得到输入信号f［k］,取M=5。MATLAB程序如下:

%program10. 3－5 R=51;%输入信号的长度 %产生(-0. 5,0. 5)的离散随机数 d=rand(1,R)-0 　　　%program10.3－5 　　　R=51;%输入信号的长度 　　　%产生(-0.5,0.5)的离散随机数 　　　d=rand(1,R)-0.5; 　　　k=0:R-1; 　　　s=2*k.*(0.9.^k); 　　　f=s+d; 　　　figure(1);stem(k,d,′·′); 　　　M=5;b=ones(M,1)/M; 　　　a=1; 　　　y=filter(b,a,f); 　　　figure(2); 　　　stem(k,s,′·′); 　　　xlabel(′Timeindexk′); 　　　legend(′s［k］′,′y［k］′);

图10.3－5　例10.3－5程序运行结果

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

例10.3－6 求离散系统 的单位脉冲响应h［k］,并与理论值h［k］=-(-1)k+2(-2)k,k≥0比较。 解 MATLAB程序如下:

　　　%program10.3－6 　　　k=0:10; 　　　a=［132］; 　　　b=［1］; 　　　h=impz(b,a,k); 　　　subplot(2,1,1) 　　　stem(k,h,′.′) 　　　%title(′单位脉冲响应的近似值′); 　　　hk=-(-1).^k+2*(-2).^k; 　　　subplot(2,1,2) 　　　stem(k,hk,′.′) 　　　%title(′单位脉冲响应的理论值′);

图10.3－6　例10.3－6程序运行结果

例10.3－7　某离散系统的差分方程为 初始条件为y［0］=0,y［1］=1,激励 ,求其单位脉冲响应、零状态响应和全响应。 解：MATLAB程序如下： %program 10.3-7 k=-10:20; a=[6 -5 1]; b=[1]; figure(1) subplot(2,1,1)

　　　impz(b,a,k)　　%－10～20范围内单位脉冲响应时域波形 　　　title(′h(k)′) 　　　%单位阶跃响应 　　　kj=0:30; 　　　Uk=ones(1,length(kj)); 　　　gk=filter(b,a,Uk); 　　　subplot(2,1,2) 　　　stem(kj,gk,′.′) 　　　title(′g(k)′) 　　　%零状态响应 　　　fk=cos(kj*pi/2); 　　　figure(2)

subplot(2,1,1),stem(kj,fk,′. ′) title(′f(k)=cos(k 　　　subplot(2,1,1),stem(kj,fk,′.′) 　　　title(′f(k)=cos(k*pi/2)′) 　　　y=filter(b,a,fk); 　　　subplot(2,1,2),stem(kj,y,′.′) 　　　title(′zerostateresponse′) 　　　%全响应 　　　y(1)=0;y(2)=1;%初始值 　　　form=3:length(kj)-2; 　　 　y(m)=(1/6)*(5*y(m-1)-y(m-2)+fk(m));　　　%递推求解 　　　end 　　　figure(3) 　　　stem(kj,y,′.′),title(′y(k)′),xlabel(′k′)

图10.3－7　例10.3－7程序运行结果

10.3.5　离散卷积和的计算 　　卷积是用来计算离散系统零状态响应的有力工具。MATLAB信号处理工具箱提供了一个计算两个离散序列卷积和的函数conv,其调用方式为 　　c=conv(a,b) 式中,a、b为待卷积和运算的两序列的向量表示,c是卷积结果。向量c的长度为向量a、b长度之和减1,即length(c)=length(a)+length(b)-1。

例10.3－8　已知序列 　　x［k］={1,2,3,4;k=0,1,2,3},y［k］={1,1,1,1,1;k=0,1,2,3,4} 　　计算x［k］*y［k］,并画出卷积和结果。 　　　解　MATLAB程序如下: 　　　%program10.3－8 　　　x=［1,2,3,4］; 　　　y=［1,1,1,1,1］; 　　　z=conv(x,y); 　　　N=length(z); 　　　Stem(0:N-1,z,′.′); 　　程序运行结果为 　　　z=1　3　6　10　10　9　7　4

图10.3－8　例10.3－8程序运行结果

10. 4 用MATLAB实现连续系统的频域分析 10. 4. 1 周期信号频域分析的MATLAB实现 例10 10.4　用MATLAB实现连续系统的频域分析 10.4.1　周期信号频域分析的MATLAB实现 　　例10.4－1　设周期性对称三角波幅度A=1,周期T=2,试用[WTBZ]MATLAB画出其频谱。

图10.4－1　周期性对称三角波

解　傅里叶系数 该信号的傅里叶级数展开式为

MATLAB程序如下: %program10. 4－1 N=8; %计算n=-N到-1的傅里叶系数 n1=-N:-1; c1=-4. j MATLAB程序如下: 　　　%program10.4－1 　　　N=8; 　　　%计算n=-N到-1的傅里叶系数 　　　n1=-N:-1; 　　　c1=-4*j*sin(n1*pi/2)/pi^2./n1.^2; 　　　%计算n=0时的傅里叶系数 　　　c0=0; 　　　%计算n=1～N的傅里叶系数 　　　n2=1:N; 　　　c2=-4*j*sin(n2*pi/2)/pi^2./n2.^2;

cn=［c1c0c2］; n=-N:N; subplot(2,1,1); stem(n,abs(cn),′ 　　　cn=［c1c0c2］; 　　　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′); 　　信号频谱图如图10.4－2所示,显然这是一个离散频谱。

图10.4－2　例10.4－1程序运行结果

例10.4－2 将图10.4－3所示方波分解为多次正弦波之和。 例10.4－2　将图10.4－3所示方波分解为多次正弦波之和。 图10.4－3　方波

解　图示周期性方波的傅里叶级数为 　　方波f(t)的周期T=2π,由于该方波是奇对称的,在t=0～π之间演示即可,分别计算 …

直到九次谐波,并绘图。MATLAB程序如下: %方波的分解 N=input(′请输入谐波次数N(奇数)=′); t=0:0. 01:2 　　直到九次谐波,并绘图。MATLAB程序如下: 　　　%方波的分解 　　　N=input(′请输入谐波次数N(奇数)=′); 　　　t=0:0.01:2*pi; 　　　y=zeros(N,max(size(t)));x=zeros(size(t)); 　　　fork=1:2:N 　　　 　x=x+sin(k*t)/k;y((k+1)/2,:)=x; 　　　end 　　　%将各次谐波叠加绘出波形图 　　　figure(1),subplot(3,1,1),plot(t,y(1:(N-1),:)),grid 　　　line(［0,pi+0.5］,［pi/4,pi/4］)　%加上方波幅度线及标注 text(pi+0.5,pi/4,′pi/4′) 　　　%分别观察各次谐波的叠加结果 　　　k=input(′请输入要观察的最高次谐波次数k=′);

　　　fori=1:(N-1) 　　　　　ifi==k 　　　　　figure(1),subplot(3,1,2),plot(t,y(i,:)),grid 　　　　　end 　　　end 　　　%figure(1),subplot(2,2,3),plot(t,y(N-1,:)),grid 　　　%基波到最高次谐波的各次谐波的叠加结果 　　　%将各半波形绘成三维网格图,看出增加谐波阶次对方波逼近程度的影响 　　　halft=ceil(length(t)/2);[DW]%只用正半周波形 　　　halft=ceil(length(t))[DW]%用整周波形 　　　figure(1),subplot(3,1,3),mesh(t(1:halft),［1:N］,y(:,1:halft))

图10.4－4　例10.4－2程序运行结果

10.4.2　非周期信号频域分析的MATLAB实现 　　MATLAB提供了许多数值计算工具,可以用来进行信号的频谱分析。quad8是计算数值积分的函数,有下面两种调用方式: 　　　y=quad8(′F′,a,b) 　　　y=quad8(′F′,a,b,［］,［］,P) 其中,F是一个字符串,表示被积函数的文件名;a、b分别表示定积分的下限和上限;P表示被积函数中的一个参数。quad8的返回值是用自适应Simpson算法得出的积分值。

例10.4－3　试用数值积分法近似计算三角波信号f(t)=(1-|t|)g2(t)的频谱。 　　解　为了用quad8计算f(t)的频谱,定义如下MATLAB函数: 　　　functiony=sf1(t,w); 　　　y=(t>=-1&t<=1).*(1-abs(t)).*exp(-j*w*t); 对不同的参数w,函数sf1将计算出傅里叶变换积分式中被积函数的值。将上述MATLAB函数用文件名sf1.m存入计算机磁盘。近似计算信号频谱的MATLAB程序为

%program10. 4－3 w=linspace(-6. pi,6 　　　%program10.4－3 　　　w=linspace(-6*pi,6*pi,512); 　　　N=length(w);F=zeros(1,N); 　　　fork=1:N 　　　F(k)=quad8(′sf1′,-1,1,［］,［］,w(k)); 　　　end 　　　figure(1); 　　　plot(w,real(F)); 　　　xlabel(′＼omega′); 　　　ylabel(′F(j＼omega)′);

图10.4－5　例10.4－3程序运行结果

例10.4－4 门信号如图10.4－6所示,试计算宽度τ=1和幅度A=1的门信号p1(t)在0～fm(Hz)频谱范围内所含的信号能量。

图10.4－6 门信号

计算上式的MATLAB程序如下: functiony=sf2(t) y=2. sinc(t). sinc(t); %program10 　　计算上式的MATLAB程序如下: 　　　functiony=sf2(t) 　　　y=2*sinc(t).*sinc(t); 　　　%program10.4－4 　　　f=linspace(0,5,256); 　　　N=length(f);w=zeros(1,N); 　　　fork=1:N 　　　 　w(k)=quad8(′sf2′,0,f(k)); 　　　end 　　　plot(f,w); 　　　xlabel(′Hz′); 　　　ylabel(′E′);

图10.4－7　例10.4－4程序运行结果

10.4.3 用MATLAB实现连续系统的频域分析 当系统的频率响应H(jω)是jω的有理真分式时,有 MATLAB信号处理工具箱提供的freqs函数可直接计算系统的频率响应,调用形式为 　　　H=freqs(b,a,w) 其中,b是分子多项式的系数向量;a为分母多项式的系数向量;w为需计算的H(jω)的抽样点角频率矩阵(数组w中最少需包含两个w的抽样点)。

例10.4－5 三阶归一化的Butterworth低通滤波器的频率响应为 试画出系统的幅频响应|H(jω)|和相频响应φ(ω)。 解　程序如下: 　　　%program10.4－5 　　　w=linspace(0,5,200); 　　　b=［1］; 　　　a=［1321］; 　　　H=freqs(b,a,w);%频率响应函数 　　　Subplot(2,1,1);

plot(w,abs(H));%幅频特性曲线 set(gca,′xtick′,［012345］); set(gca,′ytick′,［00 　　　plot(w,abs(H));%幅频特性曲线 　　　set(gca,′xtick′,［012345］); 　　　set(gca,′ytick′,［00.40.7071］); 　　　grid; 　　　xlabel(′＼omega′); 　　　ylabel(′|H(j＼omega)|′); 　　　Subplot(2,1,2); 　　　plot(w,angle(H));%相频特性曲线 　　　set(gca,′xtick′,［012345］); 　　　grid; 　　　xlabel(′＼omega′); 　　　ylabel(′phi(＼omega)′);

图10.4－8　幅频特性、相频特性

例10.4－6 RC电路如图10.4－9所示,求该电路在图10.4－10所示周期性矩形波信号作用时系统的响应。

图10.4－10　周期性矩形波信号

解　该系统的频率响应函数为 周期矩形波形的傅里叶系数为

代入A=1,τ=2,T＝4, ,得Cn=0.5Sa(0.5πn)。系统的输出响应为

MATLAB计算系统响应的程序如下,程序运行结果见图10. 4－11: %program10. 4－6 T=4;w0=2 MATLAB计算系统响应的程序如下,程序运行结果见图10.4－11: 　　　%program10.4－6 　　　T=4;w0=2*pi/T; 　　　t=-6:0.01:6; 　　　N=51; 　　　c0=0.5;xN=c0*ones(1,length(t));%dccomponent 　　　RC=0.1; 　　　forn=1:2:N%evenharmonicsarezero 　　　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 　　　xN=xN+H*cos(w0*n*t+phi)*sinc(n*0.5); 　　　end 　　　plot(t,xN); 　　　title(［′RC=′,num2str(RC)］); 　　　xlabel(′t(s)′); 　　　ylabel(′y(t)′);

图10.4－11　例10.4－6程序运行结果

10.5　用MATLAB实现连续系统的S域分析 10.5.1　MATLAB实现部分分式展开 　　用MATLAB函数residue可以得到复杂S域表示F(s)的部分分式展开式,其调用形式为 　　［r,p,k］=residue(num,den) 其中,num、den分别为F(s)分子多项式和分母多项式的系数向量;r为部分分式的系数;p为极点;k为多项式的系数,若F(s)为真分式,则k为零。

例10.5－1 用部分分式展开法求F(s)的反变换: 解　程序如下: 　　　%program10.5－1 　　　formatrat　　　 　%结果数据以分数的形式输出 　　　num=［12］; 　　　den=［1430］; 　　　［r,p］=residue(num,den) 　　运行结果: 　　　r=-1/6　　-1/2　　　2/3 　　　p=-3 　　-1 0

因此F(s)可展开为 所以 　　有时F(s)表达式中分子多项式B(s)和分母多项式A(s)是以因子相乘的情况出现时,这时可用conv函数将因子相乘的形式转换成多项式的形式: 　　C=conv(A,B)其中,A和B是两因子多项式的系数向量,C是因子相乘所得多项式的系数向量。 　　如果已知多项式的根,则可以利用poly函数将根式转换成多项式,其调用形式为 　　B=poly(A)式中,A为多项式的根,B为多项式的系数向量。

例10.5－2　用部分分式展开法求F(s)的反变换 解　程序如下: 　　　%program10.5－2 　　　num=［2305］; 　　　den=conv(［11］,［112］); 　　　［r,p,k］=residue(num,den) 　　运行结果为 　　　r=-2.000+1.1339i-2.000-1.1339i3.000 　　　p=-0.5000+1.3229i-0.5000-1.3229i-1.000 　　　k=2

　　由于上述留数r中有一对共轭复数,因此求时域表达式的计算较复杂。为了得到简洁的时域表达式,可以用cart2pol函数把共轭复数表示成模和相角形式,其调用形式为 　　［TH,R］=cart2pol(X,Y) 表示将笛卡儿坐标转换成极坐标,X、Y分别为笛卡儿坐标的横坐标和纵坐标。TH是极坐标的相角,单位为弧度;R是极坐标的模。

　　因此在上述程序中增加下面语句,即可得到留数r的极坐标形式: 　　［angle,mag］=cart2pol(real(r),imag(r)) 运行结果为 　　angle=2.6258　　-2.6258　　　0 　　mag=2.2991 2.2991 3.0000 由此可得 所以

10.5.2　H(s)的零极点与系统特性的MATLAB计算 　　系统函数H(s)通常是一个有理真分式,其分子分母均为多项式。MATLAB中提供了一个计算分子和分母多项式根的函数roots。例如多项式N(s)=s4+3s2+5s+7的根,可由如下语句求出: 　　　N=［10357］ 　　　r=roots(N)运行结果为 　　　r＝0.8693+1.9219i　0.8693-1.9219i　-0.8693+0.9042i　-0.8693-0.9042i 调用时要注意N(s)中3次幂的系数为零。

例10.5－3　已知系统函数为 求系统的零极点并画出零极点分布图。 　　解　程序如下： 　　　%program10.5－3 　　　b=［1-1］; 　　　a=［143］; 　　　zs=roots(b); 　　　ps=roots(a); plot(real(zs),imag(zs),′o′,real(ps),imag(ps),′rx′,′markersize′,12); 　　　axis(［-42-22］);grid; 　　　legend(′零点′,′极点′);

图10.5－1　例10.5－3程序运行结果

　　MATLAB中提供了一种更简便的画出系统函数零极点图的方法,即直接应用pzmap函数画图。pzmap函数调用形式为 　　pzmap(sys) 表示画出sys所描述系统的零极点图。LTI系统模型sys要借助tf函数获得,其调用方式为 　　sys=tf(b,a) 式中,b和a分别为系统函数H(s)分子多项式和分母多项式的系数向量。因此,上例还可用下述程序实现:

%program10. 5－3 b=［1-1］; a=［143］; sys=tf(b,a) pzmap(sys) 得到的零极点分布图如图10 　　　%program10.5－3 　　　b=［1-1］; 　　　a=［143］; 　　　sys=tf(b,a) 　　　pzmap(sys) 得到的零极点分布图如图10.5－2所示。 　　如果已知系统函数H(s),求系统的单位冲激响应h(t)和频率响应H(jω)可以用impulse函数和freqs函数。

图10.5－2　零极点分布图

例10.5－4　已知系统函数为 　　试画出系统的零极点分布图,求系统的单位冲激响应h(t)和频率响应H(jω),并判断系统是否稳定。

　　解　程序如下: 　　　　　%program10.5－4 　　　　　num=［1］; 　　　　　den=［1231］; 　　　　　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(′ImpulseRespone′) 　　　　　［H,w］=freqs(num,den); 　　　　　figure(3);plot(w,abs(H)) 　　　　　xlabel(′＼omega′) 　　　　　title(′MagnitudeRespone′)

运行结果为: poles=-1. 0000 -0. 5000 +0. 8660i -0. 5000 -0. 8660i 图10 运行结果为: poles=-1.0000　　-0.5000　　+0.8660i　　-0.5000　　-0.8660i 　　图10.5－3为例10.5－4程序的运行结果。三个极点均位于S平面左半开平面上,故该系统是稳定系统。系统函数的零极点分布图、系统的单位冲激响应和频率响应分别如图所示。

图10.5－3　例10.5－4程序运行结果

10.5.3　用MATLAB计算拉普拉斯变换 　　MATLAB的符号数学工具箱提供了计算拉普拉斯正反变换的函数laplace和ilaplace,其调用形式为 　　　F=laplace(f) 　　　f=ilaplace(F) 式中f为信号的时域表达式的符号对象,F表示信号f的象函数表达式的符号对象。符号对象可以应用函数sym实现,其调用格式为 　　　S＝sym(A)式中,A为待分析表示式的字符串;S为符号数字或变量。

例10. 5－5 试分别用laplace和ilaplace函数求:. (1)f(t)=e-tsin(at)ε(t)的拉普拉斯变换;. (2) 例10.5－5　试分别用laplace和ilaplace函数求: (1)f(t)=e-tsin(at)ε(t)的拉普拉斯变换; (2) 的拉普拉斯反变换。 解　(1)程序如下: 　　　%program10.5－5(1) 　　　f=sym(′exp(-t)*sin(a*t)′); 　　　F=laplace(f)运行结果为 　　　F=a/((s+1)^2+a^2)

　　(2)程序如下： 　　　%program10.5－5(2) 　　　F=sym(′s^2/(s^2+1)′); 　　　ft=ilaplace(F) 运行结果为 　　　ft=Dirac(t)-sin(t)

10.6 用MATLAB实现离散系统的Z域分析 10.6.1 部分分式展开的MATLAB实现 信号的Z域表示式通常可用下面的有理分式表示 　　为了能从信号的Z域象函数方便地得到其时域原函数,可以将F(z)展开成部分分式之和的形式,再对其取Z逆变换。

　　MATLAB的信号处理工具箱提供了一个对F(z)进行部分分式展开的函数[WTBZ]residuez它的调用形式如下: 　　［r,p,k］=residuez(num,den) 其中,num、den分别表示F(z)的分子和分母多项式的系数向量;r为部分分式的系数;p为极点;k为多项式的系数。若F(z)为真分式,则k为零。借助residuze函数可以将F(z)展开成

　　例10.6－1　试用MATLAB计算 的部分分式展开。 　　解　计算部分分式展开的[WTBZ]MATLAB程序如下: 　　　%program10.6－1 　　　num=［18］; 　　　den=［183-4-1］; 　　　［r,p,k］=residuez(num,den) 程序运行结果为 　　　r=0.36000.24000.4000 　　　p=0.5000-0.3333-0.3333 　　　k=［］

　　从运行结果中可以看出p(2)=p(3),表示系统有一个二阶的重极点,r(2)表示一阶极点前的系数,而r(3)就表示二阶极点前的系数。对高阶重极点,其表示方法是完全类似的,所以该F(z)的部分分式展开为

10.6.2　利用MATLAB计算H(z)的零极点与系统特性 　　如果系统函数H(z)的有理函数表示形式为 　　那么系统函数的零点和极点可以通过MATLAB函数roots得到,也可用函数tf2zp得到,tf2zp的调用形式为 [z,p,k]=tf2zp(b,a) 式中,b和a分别为H(z)的分子多项式和分母多项式的系数向量,它的作用是将H(z)的有理函数表示式转换为零点、极点和增益常数的表示式,即

例10.6－2 已知一离散因果LTI系统的系统函数为 求该系统的零极点。 　　解　将系统函数改写为

用tf2zp函数求系统的零极点,程序如下: %program10. 6－2 b=［121］; a=［1-0. 5-0. 0050 　　用tf2zp函数求系统的零极点,程序如下: 　　　%program10.6－2 　　　b=［121］; 　　　a=［1-0.5-0.0050.3］; 　　　［r,p,k］=tf2zp(b,a)运行结果为 　　　　　r=-1　　-1 　　　　　p=0.5198+0.5346i　　0.5198-0.5346i　　-0.5396 　　　　　k=1

　　若要获得系统函数H(z)的零极点分布图,可以直接应用zplane函数,其调用形式为 　　zplane(b,a)式中,b和a分别为H(z)的分子多项式和分母多项式的系数向量。它的作用是在Z平面画出单位圆、零点和极点。 　　如果已知系统函数H(z),求系统的单位脉冲响应h［k］和频率响应H(ejΩ),则可应用impz函数和freqz函数。

例10.6－3 已知一离散因果LTI系统的系统函数为 　　试画出系统的零极点分布图,求系统的单位脉冲响应h［k］和频率响应H(ejΩ),并判断系统是否稳定。 　　解　根据已知的H(z),用zplane函数即可画出系统的零极点分布图。利用impz函数和freqz函数求系统的单位脉冲响应和频率响应时,需要将H(z)改写成

程序如下: %program10. 6－3 b=［121］; a=［1-0. 5-0. 0050 　　　程序如下: 　　　　　%program10.6－3 　　　　　b=［121］; 　　　　　a=［1-0.5-0.0050.3］; 　　　　　figure(1);zplane(b,a); 　　　　　num=［0121］; 　　　　　den=［1-0.5-0.0050.3］; 　　　　　h=impz(num,den); 　　　　　figure(2);stem(h,′.′) 　　　　　xlabel(′k′) 　　　　　title(′ImpulseRespone′) 　　　　　［H,w］=freqz(num,den); 　　　　　figure(3);plot(w/pi,abs(H)) 　　　　　xlabel(′Frequency＼omega′) 　　　　　title(′MagnitudeRespone′)

图10.6－1　例10.6－3运行结果

10.6.3　利用MATLAB计算Z变换 　　MATLAB的符号数学工具箱提供了计算Z变换的函数ztrans和Z逆变换的函数iztrans,其调用形式为 　　　F=ztrans(f) 　　　f=iztrans(F) 式中,f为信号的时域表达式的符号对象,F表示信号f的象函数表达式的符号对象。符号对象可以应用函数sym实现,其调用格式为 　　　S＝sym(A) 式中,A为待分析表示式的字符串;S为符号数字或变量。

例10. 6－4 试分别用ztrans函数和iztrans函数求: (1)f［k］=cos(ak)ε(k)的Z变换; (2) 例10.6－4　试分别用ztrans函数和iztrans函数求: 　　(1)f［k］=cos(ak)ε(k)的Z变换; 　　(2) 的Z逆变换。 　　解　(1)求f［k］的Z变换的程序如下: 　　%program10.6－4(1) 　　f=sym(′cos(a*k)′); 　　F=ztrans(f)运行结果为 　　F=(z-cos(a))*z/(z^2-2*z*cos(a)+1) 即

(2)求F(z)逆变换的程序为 %program10 　　(2)求F(z)逆变换的程序为 　　　%program10.6－4(2) 　　　F=sym(′1/(1+z)^2′); 　　　f=iztrans(F) 程序运行结果为 　　　f=Delta(n)+(-1)^n*n-(-1)^n 即

10.7　MATLAB在系统状态空间分析中的应用 10.7.1　系统微分方程到状态空间方程的转换 　　　MATLAB提供了一个tf2ss函数,它能把描述系统的微分方程转换为等价的状态空间方程,调用形式如下 [A,B,C,D]=tf2ss(num,den) 其中num，den分别表示系统函数H(s)的分子和分母多项式，A,B,C,D分别为状态空间方程的系数矩阵。

例10.7－1 一系统的微分方程为 y″(t)+5y′(t)+10y(t)=f(t) 则该系统的H(s)为 MATLAB实现状态空间方程的语句为 [A,B,C,D]=tf2ss([1],[1 5 10]) 可得 D=0 所以,系统的状态空间方程为

10.7.2　系统状态空间方程到系统函数矩阵H(s)的计算 　　利用MATLAB提供的函数ss2tf,可以计算出由状态空间方程得出的系统函数矩阵H(s),调用形式如下 [num,den]=ss2tf(A,B,C,D,k) 其中A,B,C,D分别表示状态空间方程的系数矩阵；k表示由函数ss2tf计算的与第k个输入相关的系统函数，即H(s)的第k列。Num表示H(s)第k列的m个元素的分子多项式，den表示H(s)公共的分母多项式。

例10.7－2 已知某连续时间系统的状态方程和输出方程为 　　例10.7－2　已知某连续时间系统的状态方程和输出方程为 求该系统的系统函数矩阵H(s)。

解：因 A=[2 3;0 -1]; B=[0 1;1 0];C=[1 1;0 -1]; D=[1 0;1 0]; 计算分别与H(s)两输入相关对应列的MATLAB语句为 [num1,den1]=ss2tf(A,B,C,D,1); [num2,den2]=ss2tf(A,B,C,D,2); 可得 num1 = 1 0 -1 den1 = num2 = 1 -2 0 1 -1 -2 0 1 1 0 0 0 den2 = 1 -1 -2

所以系统函数矩阵H(s)为

10.7.3 用MATLAB求解连续时间系统的状态空间方程 连续时间系统的状态空间方程的一般形式为 用sys=ss(A,B,C,D)表示状态空间方程的计算机模型，然后由lsim函数获得状态空间方程的数值解。Lsim函数的调用形式为 [y,to,x]=lsim(sys,f,t,x0) 式中:sys表示由函数ss构造的状态空间方程模型;t表示需计算的输出样本点,t=0:dt:Tfinal;f(:,k)表示系统第k个输入在t上的抽样值;x0表示系统的初始状态(可缺省);y(:,k)表示系统的第k个输出;to表示实际计算时所用的样本点;x表示系统的状态。

例10.7－3　连续时间系统的初始状态和输入分别为 试用MATLAB计算其数值解。 解　%program10.7－3 　　A[WTBZ]=［23;0-1］;B[WTBZ]=［01;10］; 　　C[WTBZ]=［11;0-1］;D[WTBZ]=［10;10］; 　　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′);

图10.7－1　例10.7－3运行结果

10.7.4 用MATLAB求解离散时间系统的状态空间方程 离散时间系统的状态空间方程的一般形式为 由sys=ss(A,B,C,D,［］)获得离散系统状态空间方程的计算机表示模型,然后由lsim函数获得其状态空间方程的数值解。lsim函数的调用形式为 　　［y,n,x］=lsim(sys,f,［］,x0) 式中:sys表示由函数ss构造的状态空间方程模型;t表示需计算的输出样本点,t=0:dt:Tfinal;f(:,k)表示系统第k个输入序列;x0表示系统的初始状态(可缺省);y(:,k)表示系统的第k个输出;n表示序列的下标;x表示系统的状态。

例10.7－4 已知离散时间系统的状态方程和输出方程为 　　例10.7－4　已知离散时间系统的状态方程和输出方程为 和 初始状态及输入为 试用MATLAB计算该系统的数值解。

　　解　程序如下： 　　　%program10.7－4 　　　A[WTBZ]=［01;-23］; 　　　B[WTBZ]=［0;1］; 　　　C[WTBZ]=［11;2-1］; 　　　D[WTBZ]=zeros(2,1); 　　　x0=［1-1］; 　　　N=10; 　　　f=ones(1,N); 　　　sys=ss(A[WTBZ],B[WTBZ],C[WTBZ],D[WTBZ],［］); 　　　y=lsim(sys,f,［］,x0); 　　　subplot(2,1,1);

y1=y(:,1)′; stem((0:N-1),y1,′ 　　　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′); 　　其数值解如图10.7－2所示。

图10.7－2　例10.7－4运行结果