信号与系统(第4版)MATLAB在信号与系统分析中的应用

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

10.0引言

一般来说,MATLAB系统包括下面五个主要部分。

(1)编程语言:是一种以矩阵和数组为基本单位的编程语言;

(2-)工作环境:包括了一系列应用工具,提供编程和调试程序的环境;

(3-)图形处理:包括绘制二维、三维图形和创建图形用户接口;

(4)数学库函数:包含了大量的数学函数,也包括复杂的功能;

(5)应用程序接口:提供接口程序,可使MATLAB与其他语言程序进行交互。

MATLAB典型的应用包括以下两方面:数值计算和符号计算;建模和动态仿真。

10.1

MATLAB基础

10.1.1

MATLAB语言的特点MATLAB语言具有以下特点:(1)编程效率高。(2-)用户使用方便。(3-)扩充能力强。

(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)编“辑Cu器rrentDirectory”窗口,其中包含当前目录选项。

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的M文件存盘后,只需在MATLAB命令编辑窗口键入example即可运行,并显示同命令窗口输入命令一样的结果。在M文件中,程序的注释是以符号“%”开始直到该行结束的部分,程序执行时会自动忽略。

2-)函数文件

如果M文件的第一行包含function,则该文件就是函数文件。每个函数文件都定义一个函数。能够像库函数一样方便地调用,从而可扩展MATLAB的功能。如果对于一类特定的问题,建立起许多函数M文件,就能形成工具箱。

从形式上看,函数文件与命令文件的区别在于:命令文件的变量在文件执行完后保留在内存中;而函数文件内定义的变量仅在函数文件内部起作用,当函数文件执行结束后,这些内部变量将被清除。

函数M文件的第一行有特殊的要求,其形式必须为

function[输出变量列表]=函数名(输入变量列表)

函数体语句;

例:functiony=f(x)

y=sin(x);

注意:函数名应该和M文件名相同。

10.1.3-学习MATLAB的基本方法

1.help命令MATLAB。

在命令窗口中使行help命令能够获得范围不同的帮助信息,例如:

(1)运行helphelp,将得到如何使用help的提示;

(2-)直接运行help,会列出可以用于help显示的所有主题(topic);

(3-)运行help(topic),可获得有关该主题的帮助,比如,想对二维图形(graph2-d)编程有所了解,可运行≫helpgraph2-d。

(4)对每个主题(topic)中的任何命令的用法,同样可以用help来查看。如对于二维图形(graph2-d)命令中的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表示1.单位阶跃信号

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

单位阶跃信号M文件如下:图10.2-1-单位阶跃信号

2-.复指数信号

复指数信号的数学模型:

若ω=0,它是实指数函数;若u=0,则为虚指数函数,其实部为余弦函数,虚部为正弦函数。图10.2-2-复指数信号实部图10.2-3-复指数信号虚部

3-.矩形脉冲

调用MATLAB的内部函数可产生一矩形脉冲信号。在MATLAB中用rectpuls函数表示矩形脉冲信号,即

y=rectpuls(t,width)

用以产生一个幅度为1,宽度为width,以t=0为对称轴的矩形波。Width的默认值为1。例如产生一个以t=2-T为对称中心的矩形脉冲信号的MATLAB程序如下,取T=1:图10.2-4矩形脉冲=

4.三角波脉冲

三角波脉冲信号在MATLAB中用tripuls函数表示,其调用形式为

y=tripuls(t,width,skew)

用以产生一个最大幅度为1,宽度为width的三角波。该函数的非零范围为(-width/2-,width/2-)。例如,用以下MATLAB语句实现三角波:图10.2-5三角波

10.2-.2-离散信号的MATLAB表示

由于MATLAB数值计算的特点,故用它来表示离散信号是非常方便的。在MATLAB中,用一个列向量来表示一个有限长序列,而这样的向量并没有包含采样位置的信息,要完全表示一个序列x(k),需用x和k两个向量,例如:

下面的箭头指示的是k=0时的采样点。该序列在MATLAB中表示为

若不需要采样位置信息或这个信息是多余的(例如该序列从k=0开始),可以只用x向量来表示。计算机的内存有限,MATLAB无法表示无限长序列。

1.单位脉冲序列

单位脉冲序列的表达式:

延迟ks

的单位脉冲序列表达式:

本例取ks=3-,此单位脉冲序列的M文件如下:图10.2-6单位脉冲序列

2-.单位阶跃序列

单位阶跃序列的表达式:

延迟ks的单位阶跃序列表达式:图10.2-7单位阶跃序列

3-.复指数序列

复指数序列的表达式:

若ω=0,它是实指数序列;若α=0,则为虚指数序列,其实部为余弦序列,虚部为正弦序列。本例取α=-0.2-,ω=0.5,该复指数序列的实部和虚部如图10.2-8所示,其M文件如下:图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(2-t)和f(2--2-t)的波形。图10.2-9例10.2-1图

2-.连续信号的微分与积分

连续信号的微分可用diff近似计算,例如y=sin'(x2-)=2-xcos(x2-)可由以下MATLAB语句近似实现:

图10.2-10例10.2-2-图

3-.离散序列的差分与求和

离散序列的差分∇f[k]=f[k]-f[k-1],在MATLAB中用diff函数实现,其调用格式为

y=diff(f)

离散序列的求和与信号相加运算不同,求和运算是把k1

和k2

之间的所有样本f[k]加起来,在MATLAB中用sum函数实现,其调用格式为

y=sum(f(k1:k2

))

例10.2-3-用MATLAB计算指数信号(-1.6)kε[k]的能量。

解离散信号的能量定义为

10.3-用MATLAB实现系统的时域分析

10.3-.1

连续系统冲激响应的求解方法1-应用微分方程。MATLAB提供了专门用于求LTI系统的冲激响应和阶跃响应的函数。系统微分方程为

例10.3-1

求以下系统的冲激响应和阶跃响应:图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)的分母多项式没有重根,则

故有冲激响应

运行结果(用通用程序求一个五阶系统的冲激响应,按提示输入分子、分母系数向量和时间数组):

因为题中给出的分母是系统的极点,而不是多项式系数,该系数可用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,其调用格式为

式中,t表示计算系统响应的抽样点向量;f是系统输入信号向量;sys是LTI系统模型,用来表示微分方程、差分方程、状态方程。在求解微分方程时,微分方程的LTI系统模型sys要借助tf函数获得,其调用方式为

sys=tf(b,a)

式中,b和a分别为微分方程右端和左端各项的系数向量。例如对三阶微分方程

可用

获得LTI模型。微分方程中系数为零也要写入向量a和b中。

例10.3-3

系统的微分方程为

求系统在输入为f(t)=10sin2-πt时的零状态响应。图10.3-3-例10.3-3-程序运行结果

例10.3-4设二阶连续系统的微分方程为

求系统的冲激响应。若输入为f(t)=3-t+cos(0.1t),求系统的零状态响应。

解求冲激响应,系统微分方程的特征方程为

其特征根为p1、p2-,相应的系数为r1、r2-,则冲激响应为

输出y(t)为输入f(t)与冲激响应h(t)的卷积。图10.3-4例10.3-4程序运行结果

10.3-.3

离散系统零状态响应的求解

LTI离散系统可用线性常系数差分方程描述:

其中f[k]、y[k]分别表示系统的输入和输出,n是差分方程的阶数。已知差分方程的n个初始状态和输入f[k],就可以编程迭代计算出系统的输出:

在零初始状态下,MATLAB工具箱提供了一个filter函数,计算由差分方程描述的系统响应,其调用格式为

式中,b=[b0,b1,b2

,…,bm],a=[a0,a1,a2

,…,an]分别是差分方程左右端的系数向量,f表示输入序列,y表示输出序列。注意输出序列和输入序列的长度相同。

例10.3-5系统的输入输出关系为

输入信号为f[k]=s[k]+d[k],其中s[k]=(2-k)0.9k,d[k]是随机信号。试用MATLAB编程求系统的零状态响应。图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比较。图10.3-6例10.3-6程序运行结

例10.3-7某离散系统的差分方程为

初始条件为y[0]=0,y[1]=1,激励求其单位脉冲响应、零状态响应和全响应。图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]*y[k],并画出卷积和结果。图10.3-8例10.3-8程序运行结果

10.4用MATLAB实现连续系统的频域分析

10.4.1周期信号频域分析的MATLAB实现例10.4-1设周期性对称三角波幅度A=1,周期T=2-,试用MATLAB画出其频谱。图10.41-周期性对称三角波图10.42-例10.41程序运行结果

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

解图示周期性方波的傅里叶级数为

方波f(t)的周期T=2-π,由于该方波是奇对称的,在t=0~π之间演示即可,分别计算图10.44例10.42-程序运行结果

10.4.2

非周期信号频域分析的MATLAB实现

MATLAB提供了许多数值计算工具,可以用来进行信号的频谱分析。quad8是计算数值积

其中,F是一个字符串,表示被积函数的文件名;a、b分别表示定积分的下限和上限;P表示被积函数中的一个参数。quad8的返回值是用自适应Simpson算法得出的积分值。

例10.4-3

试用数值积分法近似计算三角波信号f(t)=(1-|t|)g2-(t)的频谱。

解为了用quad8计算f(t)的频谱,定义如下MATLAB函数:

对不同的参数w,函数sf1将计算出傅里叶变换积分式中被积函数的值。将上述MATLAB函数用文件名sf1.m存入计算机磁盘。近似计算信号频谱的MATLAB程序为图10.45例10.43-程序运行结果

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

解该门信号的频谱为

所以信号在0~fm(Hz)频谱范围内所含的信号能量图10.47例10.44程序运行结果

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ω)|和相频响应φ(ω)。图10.48幅频特性、相频特性

例10.4-6RC电路如图10.4-9所示,激励为f(t),响应为y(t)。求该电路在图10.4-10所示周期性矩形波信号作用时系统的响应。图10.49RC电路图10.410周期性矩形波信号

解该系统的频率响应函数为

周期矩形波形的傅里叶系数为图10.411-例10.46程序运行结果

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为零。

有时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)的反变换

由于上述留数r中有一对共轭复数,因此求时域表达式的计算较复杂。为了得到简洁的时域表达式,可以用cart2-pol函数把共轭复数表示成模和相角形式,其调用形式为

[TH,R]=cart2-pol(X,Y)

表示将笛卡儿坐标转换成极坐标,X、Y分别为笛卡儿坐标的横坐标和纵坐标。TH是极坐标的相角,单位为弧度;R是极坐标的模。

10.5.2

H(s)的零极点与系统特性的MATLAB计算

系统函数H(s)通常是一个有理真分式,其分子分母均为多项式。MATLAB中提供了一个计算分子和分母多项式根的函数roots。例如多项式N(s)=s4+3-s2-+5s+7的根,可由如下语句求出:

运行结果为

调用时要注意N(s)中3次幂的系数为零。

例10.5-3已知系统函数为

求系统的零极点并画出零极点分布图。图10.51-例10.53-程序运行结果

MATLAB中提供了一种更简便的画出系统函数零极点图的方法,即直接应用pzmap函数画图。pzmap函数调用形式为

pzmap(sys)

表示画出sys所描述系统的零极点图。LTI系统模型sys要借助tf函数获得,其调用方式为

sys=tf(b,a)

式中,b和a分别为系统函数H(s)分子多项式和分母多项式的系数向量。因此,上例还可用下述程序实现:

得到的零极点分布图如图10.5-2所示。图10.0.52-零极点分布图

例10.5-4已知系统函数为

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

运行结果为:

图10.5-3为例10.5-4程序的运行结果。三个极点均位于S平面左半开平面上,故该系统是稳定系统。系统函数的零极点分布图、系统的单位冲激响应和频率响应分别如图所示。图10.53-例10.54程序运行结果

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函数求:

10.6用MATLAB实现离散系统的Z域分析

10.6.1部分分式展开的MATLAB实现信号的Z域表示式通常可用下面的有理分式表示

为了能从信号的Z域象函数方便地得到其时域原函数,可以将F(z)展开成部分分式之和的形式,然后取其Z逆变换。MATLAB的信号处理工具箱提供了一个对F(z)进行部分分式展开的函数residuez,它的调用形式如下:

其中,num、den分别表示F(z)的分子和分母多项式的系数向量;r为部分分式的系数;p为极点;k为多项式的系数。若F(z)为真分式,则k为零。借助residuze函数可以将F(z)展开成

10.6.2利用MATLAB计算H(z)的零极点与系统特性

如果系统函数H(z)的有理函数表示形式为

那么系统函数的零点和极点可以通过MATLAB函数roots得到,也可用函数tf2-zp得到,tf2-zp的调用形式为

[z,p,k]=tf2-zp(b,a)

式中,b和a分别为H(z)的分子多项式和分母多项式的系数向量,它的作用是将H(z)的有理函数表示式转换为零点、极点和增益常数的表示式,即

例10.6-2已知一离散因果LTI系统的系统函数为

求该系统的零极点。

解将系统函数改写为

若要获得系统函数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)改写成

程序运行结果如图10.6-1所示。系统函数的零极点分布图中符号o表示零点,符号o旁的数字表示零点的阶数;符号×表示极点。图中的虚线画的是单位圆。由图可知,该因果系统的极点全在单位圆内,故系统是稳定的。图10.61-例10.63-运行结果

10.6.3利用MATLAB计算Z变换

MATLAB的符号数学工具箱提供了计算Z变换的函数ztrans和Z逆变换的函数iztran