Matlab在信号处理中的应用精要

第10章MATLAB在信号与系统分析中的应用第10章

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

引

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

MATLAB基础10.1.1 MATLAB语言的特点MATLAB语言具有以下特点：(1)编程效率高。MATLAB编程语言作为面向科学与工程计算的高级语言,允许用数学形式的语言编写程序,且比Basic、Fortran和C等语言更加接近我们书写计算公式的思维方式。用MATLAB编写程序犹如在演算纸上排列出公式与求解问题,因此,MATLAB语言也可通俗地称为演算纸式科学算法语言,它编写简单,编程效率高,易学易懂。第10章MATLAB在信号与系统分析中的应用(2)用户使用方便。MATLAB语言是一种解释执行的语言(在没被专门的工具编译之前),它灵活、方便,其调试程序手段丰富。MATLAB运行时,在命令行每输入一条MATLAB语句(命令),

包括调用M文件的语句,计算机就立即对其进行处理,完成编译、连接和运行的全过程。在运行m文件时,如果有错,计算机屏幕提示出错信息,经用户修改后再执行,直到正确为止。第10章MATLAB在信号与系统分析中的应用(3)扩充能力强。高版本的MATLAB语言有丰富的库函数,在进行复杂的

数学运算时可以直接调用。用户可以根据需要建立和扩充新的库函数,以提高MATLAB的使用效率,扩充其功能。第10章MATLAB在信号与系统分析中的应用(4)语句简单,内涵丰富。MATLAB语言中最基本、最重要的成分是函数,其一般形式为［

a,b,c

,…］=

fun

(

d,e,f

,…)即一个函数由函数名,输入变量d,e,f,…和输出变量a,b,c…组成。同一函数名F,可以有不同数目的输入变量(包括无输入变量)及不同数目的输出变量,代表着不同的含义。这不仅使MATLAB的库函数功能更丰富,而且大大减少了需要的磁

盘空间,使得MATLAB

编写的m文件简单、短小而高效。第10章MATLAB在信号与系统分析中的应用(5)高效方便的矩阵和数组运算。MATLAB语言像Basic、Fortran和C语言一样规定了矩

阵的算术运算符、关系运算符、逻辑运算符、条件运算符及赋值运算符,而且这些运算符大部分可以毫无改变地运用到数组间的运算中,有些运算符(如算术运算符)只要增加“·”就可用于数组间的运算。第10章MATLAB在信号与系统分析中的应用(6)方便的绘图功能。MATLAB有一系列绘图函数(命令),调用不同的绘图函数可方便地绘制线性坐标、对数坐标、半对数坐标及极坐标,通过相应的命令还可以在图上标出图题、XY轴标注、格(栅)等。总之,MATLAB语言的设计思想可以说代表了当前计算机高级语言的发展方向,读者在不断使用中会发现其具有巨大的潜力。第10章MATLAB在信号与系统分析中的应用MATLAB工作环境简介启动MATLAB有三种方法启动MATLAB:双击Windows桌面上的MATLAB快捷图标;通过“开始”菜单的“程序”子菜单中的MATLAB项启动;在MATLAB目录中搜索到可执行程序MATLAB.exe,双击该程序使之启动。启动后,MATLAB主界面如图10.1－1所示。第10章MATLAB在信号与系统分析中的应用图10.1－1MATLAB主界面第10章MATLAB在信号与系统分析中的应用MATLAB主界面大致包括以下几个部分:(1)菜单项;工具栏;“CommandWindow”(命令窗口),在提示符

后直接输入命令可以执行相关的命令;“LaunchPad”(分类帮助文件夹);

(5)“Workspace”(工作空间)，该程序窗口中列出了程序计算过程中产生的变量及其对应的数据的尺寸、字节和类型。选中一个变量,单击鼠标右键则可根据菜单进行相应的操作。第10章MATLAB在信号与系统分析中的应用“CommandHistory”(命令的历史记录)窗口，该窗口记录用户每次开启MATLAB

的时间,以及每次开启MATLAB后在MATLAB命令窗口中运行过的所有命令行。这些命令行记录可以被复制到命令窗口中再运行。选中该窗口中的任一命令记录,然后单击鼠标右键,则可根据弹出的菜单进行相应的操作。“CurrentDirectory”窗口，其中包含当前目录选项。第10章MATLAB在信号与系统分析中的应用2.程序编辑器1)命令文件命令文件没有输入参数,也不返回输出参数,只是一些命令行的组合。命令文件中的语句可以访问MATLAB工作空间(Workspace)中的所有数据,在运行的过程中所产生的变量均是全局变量。这些变量一旦生成,就一直保存在内存空间

中,除非用户将它们清除(用clear命令)。运行一个命令文件等价于从命令窗口中按顺序连续执行文件中的命令。由于命令文件只是一串命令的结合,因此程序不需要预先定义,而只需按命令窗口中的命令输入顺序,依次将命令编辑在命令文件

中即可。如果某个命令不需要显示结果,则在该命令后加上“;”。注意文件名一定是“.m”。命令文件的建立过程如

下:第10章MATLAB在信号与系统分析中的应用进入程序编辑器(MATLABEditor／Debug):从“File”菜单中选择“New”及“m

file”或单击“Newm

file”按钮;输入程序:在“MATLABEditor/Debug”窗口输入

MATLAB程序;保存程序:单击“Save”按钮,出现一个对话框,在文件名框中键入一个文件名,单击“保存”按钮，一个m文件便保存在磁盘上了。第10章MATLAB在信号与系统分析中的应用运行命令文件时,该m文件中的命令可以访问

MATLAB工作区中的所有变量,而且其中的所有变量也成为工作区的一部分。命令文件运行结束,所产生的变量保留在工作区,直至关闭MATLAB或用命令删除。下面是一个命令文件的例子,程序如下:%文件名example.mx=1;y=2;z=3items=x+y+zcost=x*5+y*2+z*9averagecost=cost/items第10章MATLAB在信号与系统分析中的应用当这个文件在程序编辑窗口输入并以名为example.m的m文件存盘后,只需在MATLAB命令编辑窗口键入example即

可运行,并显示同命令窗口输入命令一样的结果。在m文件中,程序的注释是以符号“％”开始直到该行结束的部分,程序执行时会自动忽略注释语句。第10章MATLAB在信号与系统分析中的应用2)函数文件如果m文件的第一行包含function,则该文件就是函数文件。每个函数文件都定义一个函数。能够像调用库函数一样方便地调用函数文件,从而可扩展MATLAB的功能。如果对于一类特定的问题,建立起许多函数m文件,就能形成工具箱。从形式上看,函数文件与命令文件的区别在于:命令文件的变量在文件执行完后保留在内存中;而函数文件内定义的变量仅在函数文件内部起作用,当函数文件执行结束后,这些内部变量将被清除。第10章MATLAB在信号与系统分析中的应用函数m文件的第一行有特殊的要求,其形式必须为function［输出变量列表］=函数名(输入变量列表)函数体语句;例:functiony=f(x)y=sin(x);注意:函数名应该和m文件名相同。第10章MATLAB在信号与系统分析中的应用学习MATLAB的基本方法1.help命令在命令窗口中使行help命令能够获得范围不同的帮助信息,例如:运行helphelp,将得到如何使用help的提示;直接运行help,会列出可以用于help显示的所有主题

(topic);运行help(topic),可获得有关该主题的帮助,比如,想对二维图形(graph2d)编程有所了解,可运行

helpgraph2d。对每个主题(topic)中的任何命令的用法,同样可以用help来查看。如对于二维图形(graph2d)命令中的plot,用help查看的方法是:

helpplot。第10章MATLAB在信号与系统分析中的应用2.lookfor命令当用户要查找具有某种功能的命令但不知道其准确名字时,help就无能为力了。而lookfor可以根据用户提供的完整或不完整的关键词,去搜索出一组与之有关的命令,用户可从列表中挑选出满足需要的命令。如利用lookfor命令查找矩阵求逆函数:>>lookforinverse第10章MATLAB在信号与系统分析中的应用doc、helpwin和helpdesk命令在命令窗口中运行doc、helpwin和helpdesk命令中的任何一个,都会打开一个名为“help”的帮助窗口。demo命令demo命令用于查看集成于MATLAB环境内的各种演示内容。在MATLAB的命令窗口键入demo命令可以得到演示界面,从而可以方便用户了解MATLAB的基本功能。帮助菜单启动MATLAB应用程序后,单击“help”主菜单,则会弹出一系列子菜单,可以根据菜单直接进行操作。第10章MATLAB在信号与系统分析中的应用10.2 信号的MATLAB表示10.2.1 连续信号的MATLAB表示严格地说,MATLAB不能处理连续信号,它是用连续信号在等间隔点的样值来近似表示连续信号的。当采样间隔足够小时,这些样值就能较好地近似表示连续信号。MATLAB提供了大量的生成基本信号的函数。最常用

的指数信号、正弦信号是MATLAB

的内部函数,即不安装任何工具箱就可调用的函数。第10章MATLAB在信号与系统分析中的应用1.单位阶跃信号单位阶跃信号的数学模型:在t=t1处跃升的阶跃信号可写为ε(t-t1),定义为第10章MATLAB在信号与系统分析中的应用单位阶跃信号m文件如下：%单位阶跃信号m文件%信号从t0到tf,在t1(t0≤t1≤tf)前为0,到t1处有一跃变,t1以后为1clear;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章MATLAB在信号与系统分析中的应用图10.2－1单位阶跃信号第10章MATLAB在信号与系统分析中的应用2.复指数信号复指数信号的数学模型:x3(t)=e(u+jω)t若ω=0,它是实指数函数;若u=0,则为虚指数函数,其实部为余弦函数,虚部为正弦函数。第10章MATLAB在信号与系统分析中的应用复指数信号[WTBZ]m文件如下：%复指数信号m文件%信号从t0到tfclear;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章MATLAB在信号与系统分析中的应用第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用图10.2－4矩形脉冲第10章MATLAB在信号与系统分析中的应用图10.2－5

三角波第10章MATLAB在信号与系统分析中的应用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无法表示无限长序列。第10章MATLAB在信号与系统分析中的应用1.单位脉冲序列单位脉冲序列的表达式:k=0k=其余延迟ks的单位脉冲序列表达式:k=k

sk=其余第10章MATLAB在信号与系统分析中的应用本例取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章MATLAB在信号与系统分析中的应用图10.2－6单位脉冲序列第10章MATLAB在信号与系统分析中的应用2.单位阶跃序列单位阶跃序列的表达式:延迟ks的单位阶跃序列表达式:第10章MATLAB在信号与系统分析中的应用本例取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章MATLAB在信号与系统分析中的应用图10.2－7单位阶跃序列第10章MATLAB在信号与系统分析中的应用3.复指数序列复指数序列的表达式:若ω=0,它是实指数序列;若α=0,则为虚指数序列,其实部为余弦序列,虚部为正弦序列。本例取α=-0.2,ω=0.5,该复指数序列的实部和虚部如图10.2－8所示，其m文件如下：第10章MATLAB在信号与系统分析中的应用％复指数序列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,100,0］)%绘图xlabel(′实部′)subplot(1,2,2),stem(k3,imag(x3),′.′);line(［0,100,0］)%绘图xlabel(′虚部′)第10章MATLAB在信号与系统分析中的应用图10.2－8复指数序列的实部、虚部第10章MATLAB在信号与系统分析中的应用10.2.3 信号基本运算的MATLAB实现1.信号的尺度变换、翻转、平移(时移)信号的尺度变换、翻转、平移运算,实际上是函数自变量的运算。在信号的尺度变换(f(at)和f［ak］)中,函数的自变量乘以一个常数,在MATLAB中可用算术运算符“*”来实现。在信号翻转运算(f(-t)和f［-k］)中,函数的自变量乘以一个负号,在MATLAB

中可以直接写出。在信号平移运算(f(t±t0)和f［k±k0］)中,函数自变量加、减一个常数,在MATLAB中可用算术运算符“+”或“-”来实现。第10章MATLAB在信号与系统分析中的应用例10.2－1

三角波f(t)如图10.2－9(a)所示,试利用MATLAB画出f(2t)和f(2-2t)的波形。解

实现f(2t)和f(2-2t)的程序如下:％program10.2－1t=-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章MATLAB在信号与系统分析中的应用图10.2－9例10.2－1图第10章MATLAB在信号与系统分析中的应用2.连续信号的微分与积分连续信号的微分可用diff近似计算,例如y=sin′(x2)=2xcos(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章MATLAB在信号与系统分析中的应用例10.2－2三角波f(t)如图10.2－10(a)所示,试利用MATLAB

画出

和

的波形。解

为便于利用quad函数计算信号的积分,将三角波f(t)写成MATLAB函数,函数名为f2_tri,程序如下:％program10.2－2functionyt=f2_tri(t)yt=tripuls(t,4,0.5);利用diff和quad函数,并调用f2_tri可实现三角波信号f(t)的微分、积分,程序如下:第10章MATLAB在信号与系统分析中的应用％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));endplot(t,y2)(a)title(′integraloff(t)′)微分、积分波形分别如图10.2－10(b)、(c)所示。第10章MATLAB在信号与系统分析中的应用图10.2-10例10.2-2图第10章MATLAB在信号与系统分析中的应用3.离散序列的差分与求和离散序列的差分MATLAB中用diff函数实现,其调用格式为,在y=diff(f)离散序列的求和

与信号相加运算不同,求和运算是把k1和k2之间的所有样本f［k］加起来,在MATLAB中用sum函数实现,其调用格式为y=sum(f(k1:k2))第10章MATLAB在信号与系统分析中的应用例10.2－3

用MATLAB计算指数信号(-1.6)kε［k］的能量。解

离散信号的能量定义为其MATLAB实现如下:％program10.2－3k=0:10;A=1;a=-1.6;fk=A*a.^k;W=sum(abs(fk).^2)运行结果为W=1.9838e+004第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用例10.3－1

求以下系统的冲激响应和阶跃响应:解

程序如下：％program10.3－1a=［746］;b=［11］;subplot(2,1,1)impulse(b,a)subplot(2,1,2)step(b,a)第10章MATLAB在信号与系统分析中的应用图10.3－1例10.3－1图第10章MATLAB在信号与系统分析中的应用方法2

应用系统函数。例10.3－2

求n阶LTI系统的冲激响应。解

设系统函数为其特性可用系统函数分子、分母系数向量b和a来表示。第10章MATLAB在信号与系统分析中的应用对于物理可实现系统，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)的分母多项式没有重根，则故有冲激响应第10章MATLAB在信号与系统分析中的应用程序如下：％program10.3－2a=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及终点tt=0:dt:tf;h=zeros(1,length(t));%h的初始化第10章MATLAB在信号与系统分析中的应用fori=1:length(a)-1

%根数为a的长度减1h=h+r(i)*exp(p(i)*t);%叠加各根分量endplot(t,h)grid运行结果(用通用程序求一个五阶系统的冲激响应,按提示输入分子、分母系数向量和时间数组):a=poly(［0,-1+2i,-1-2i,-2,-5］);b=［8,3,1］;t=0:0.2:8;第10章MATLAB在信号与系统分析中的应用因为题中给出的分母是系统的极点,而不是多项式系数,故要求出其系数可用poly函数,其格式为a=poly(p)(其中p为极点向量)即可求出h,画出的曲线如图10.3－2所示。第10章MATLAB在信号与系统分析中的应用图10.3－2例10.3－2图(5阶LTI系统的冲激响应)第10章MATLAB在信号与系统分析中的应用10.3.2 连续系统零状态响应的求解方法1

应用MATLAB工具箱提供的函数lsim。LTI连续系统以常系数微分方程描述,系统的零状态响应可通过求解初始状态为零的微分方程得到。在MATLAB中,控制系统工具箱提供了一个用于求解零初始条件微分方程数值解的函数lsim,其调用格式为y=lsim(sys,f,t)式中,t表示计算系统响应的抽样点向量;f是系统输入信号向量;sys是LTI系统模型,用来表示微分方程、差分方程、状态方程。第10章MATLAB在信号与系统分析中的应用在求解微分方程时,微分方程的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章MATLAB在信号与系统分析中的应用例10.3－3

系统的微分方程为求系统在输入为

时的零状态响应。解

MATLAB程序如下:％program10.3－3ts=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章MATLAB在信号与系统分析中的应用图10.3－3例10.3－3程序运行结果第10章MATLAB在信号与系统分析中的应用方法2

应用公式yf(t)=h(t)*f(t)。例10.3－4

设二阶连续系统的微分方程为求系统的冲激响应。若输入为f(t)=3t+cos(0.1t),求系统的零状态响应。第10章MATLAB在信号与系统分析中的应用解

求冲激响应,系统微分方程的特征方程为其特征根为p1、p2,相应的系数为r1、r2,则冲激响应为输出y(t)为输入f(t)与冲激响应h(t)的卷积。第10章MATLAB在信号与系统分析中的应用程序如下：％program10.3－4clf,cleara=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);第10章MATLAB在信号与系统分析中的应用%用极点留数法求冲激响应［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章MATLAB在信号与系统分析中的应用图10.3－4例10.3－4程序运行结果第10章MATLAB在信号与系统分析中的应用10.3.3 离散系统零状态响应的求解LTI离散系统可用线性常系数差分方程描述:其中f［k］、y［k］分别表示系统的输入和输出,n是差分方程的阶数。已知差分方程的n个初始状态和输入f［k］,就可以编程迭代计算出系统的输出:第10章MATLAB在信号与系统分析中的应用在零初始状态下,MATLAB工具箱提供了一个filter函数,计算由差分方程描述的系统响应,其调用格式为y=filter(b,a,f)式中,b=［b0,b1,b2,…,bm］,a=［a0,a1,a2,…,an］分别是差分方程左右端的系数向量,f表示输入序列,y表示输出序列。注意输出序列和输入序列的长度相同。第10章MATLAB在信号与系统分析中的应用例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程序如下:第10章MATLAB在信号与系统分析中的应用%program10.3－5R=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章MATLAB在信号与系统分析中的应用图10.3－5例10.3－5程序运行结果第10章MATLAB在信号与系统分析中的应用10.3.4 离散系统单位脉冲响应的求解在MATLAB中,求解离散系统单位脉冲响应,可用信号处理工具箱提供的函数impz,

其调用方式为h=impz(b,a,k)其中,b=［b0,b1,b2,…,bn］,a=［a0,a1,a2,…,an］分别是差分方程左、右的系数向量,k表示输出序列的取值范围,h就是系统的单位脉冲响应。第10章MATLAB在信号与系统分析中的应用例10.3－6

求离散系统的单位脉冲响应h［k］,并与理论值h［k］=-(-1)k+2(-2)k,k≥0比较。解

MATLAB程序如下:第10章MATLAB在信号与系统分析中的应用%program10.3－6k=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章MATLAB在信号与系统分析中的应用图10.3－6例10.3－6程序运行结果第10章MATLAB在信号与系统分析中的应用例10.3－7

某离散系统的差分方程为初始条件为y［0］=0,y［1］=1,激励求其单位脉冲响应、零状态响应和全响应。解：MATLAB程序如下：%program

10.3-7k=-10:20;a=[6

-5

1];b=[1];figure(1)subplot(2,1,1),第10章MATLAB在信号与系统分析中的应用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)第10章MATLAB在信号与系统分析中的应用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));endfigure(3)%递推求解stem(kj,y,′.′),title(′y(k)′),xlabel(′k′)第10章MATLAB在信号与系统分析中的应用图10.3－7例10.3－7程序运第10章MATLAB在信号与系统分析中的应用10.3.5离散卷积和的计算卷积是用来计算离散系统零状态响应的有力工具。MATLAB信号处理工具箱提供了一个计算两个离散序列卷积和的函数conv,其调用方式为c=conv(a,b)式中,a、b为待卷积和运算的两序列的向量表示,c是卷积结果。向量c的长度为向量a、b长度之和减1,即length(c)=length(a)+length(b)-1。第10章MATLAB在信号与系统分析中的应用例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－8x=［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章MATLAB在信号与系统分析中的应用图10.3－8例10.3－8程序运行结果第10章MATLAB在信号与系统分析中的应用用MATLAB实现连续系统的频域分析周期信号频域分析的MATLAB实现例10.4－1

设周期性对称三角波幅度A=1,周期T=2,试用[WTBZ]MATLAB画出其频谱。第10章MATLAB在信号与系统分析中的应用图10.4－1周期性对称三角波第10章MATLAB在信号与系统分析中的应用解

傅里叶系数该信号的傅里叶级数展开式为第10章MATLAB在信号与系统分析中的应用MATLAB程序如下:%program10.4－1N=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;第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用图10.4－2例10.4－1程序运行结果第10章MATLAB在信号与系统分析中的应用例10.4－2

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

方波第10章MATLAB在信号与系统分析中的应用解

图示周期性方波的傅里叶级数为方波f(t)的周期T=2π,由于该方波是奇对称的,在t=0～

π之间演示即可,分别计算…第10章MATLAB在信号与系统分析中的应用直到九次谐波,并绘图。MATLAB程序如下:%方波的分解N=input(′请输入谐波次数N(奇数)=′);t=0:0.01:2*pi;y=zeros(N,max(size(t)));x=zeros(size(t));fork=1:2:Nx=x+sin(k*t)/k;y((k+1)/2,:)=x;end%将各次谐波叠加绘出波形图

figure(1),subplot(3,1,1),plot(t,y(1:(N-1),:)),gridline(［0,pi+0.5］,［pi/4,pi/4］)text(pi+0.5,pi/4,′pi/4′)%分别观察各次谐波的叠加结果%加上方波幅度线及标注k=input(′请输入要观察的最高次谐波次数k=′);第10章MATLAB在信号与系统分析中的应用fori=1:(N-1)ifi==kfigure(1),subplot(3,1,2),plot(t,y(i,:)),gridendend%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:half第10章MATLAB在信号与系统分析中的应用图10.4－4例10.4－2程序运行结果第10章MATLAB在信号与系统分析中的应用10.4.2 非周期信号频域分析的MATLAB实现MATLAB提供了许多数值计算工具,可以用来进行信号

的频谱分析。quad8是计算数值积分的函数,有下面两种调用方式:y=quad8(′F′,a,b)y=quad8(′F′,a,b,［］,［］,P)其中,F是一个字符串,表示被积函数的文件名;a、b分别表示定积分的下限和上限;P表示被积函数中的一个参数。quad8的返回值是用自适应Simpson算法得出的积分值。第10章MATLAB在信号与系统分析中的应用例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程序为第10章MATLAB在信号与系统分析中的应用%program10.4－3w=linspace(-6*pi,6*pi,512);N=length(w);F=zeros(1,N);fork=1:NF(k)=quad8(′sf1′,-1,1,［］,［］,w(k));endfigure(1);plot(w,real(F));xlabel(′＼omega′);ylabel(′F(j＼omega)′);第10章MATLAB在信号与系统分析中的应用图10.4－5例10.4－3程序运行结果第10章MATLAB在信号与系统分析中的应用例10.4－4

门信号如图10.4－6所示,试计算宽度τ=1和幅度A=1的门信号p1(t)在0～fm(Hz)频谱范围内所含的信号能量。所以信号在0～fm(Hz)频谱范围内所含的信号能量第10章MATLAB在信号与系统分析中的应用图10.4－6门信号第10章MATLAB在信号与系统分析中的应用计算上式的MATLAB程序如下:

functiony=sf2(t)y=2*sinc(t).*sinc(t);%program10.4－4f=linspace(0,5,256);N=length(f);w=zeros(1,N);fork=1:Nw(k)=quad8(′sf2′,0,f(k));endplot(f,w);xlabel(′Hz′);ylabel(′E′);第10章MATLAB在信号与系统分析中的应用图10.4－7例10.4－4程序运行结果第10章MATLAB在信号与系统分析中的应用10.4.3 用MATLAB实现连续系统的频域分析当系统的频率响应H(jω)是jω的有理真分式时,有MATLAB信号处理工具箱提供的freqs函数可直接计算系统的频率响应,调用形式为H=freqs(b,a,w)其中,b是分子多项式的系数向量;a为分母多项式的系数向量;w为需计算的H(jω)的抽样点角频率矩阵(数组w中最少需包含两个w的抽样点)。第10章MATLAB在信号与系统分析中的应用例10.4－5

三阶归一化的Butterworth低通滤波器的频率响应为试画出系统的幅频响应|H(jω)|和相频响应φ(ω)。解

程序如下:%program10.4－5w=linspace(0,5,200);b=［1］;a=［1321］;H=freqs(b,a,w);%频率响应函数Subplot(2,1,1);第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用图10.4－8幅频特性、相频特性第10章MATLAB在信号与系统分析中的应用例10.4－6

RC电路如图10.4－9所示,求该电路在图10.4－10所示周期性矩形波信号作用时系统的响应。图10.4－9RC电路第10章MATLAB在信号与系统分析中的应用图10.4－10周期性矩形波信号第10章MATLAB在信号与系统分析中的应用解

该系统的频率响应函数为周期矩形波形的傅里叶系数为第10章MATLAB在信号与系统分析中的应用代入A=1,τ=2,T＝4,,得Cn=0.5Sa(0.5πn)。系统的输出响应为第10章MATLAB在信号与系统分析中的应用MATLAB计算系统响应的程序如下,程序运行结果见图10.4－11:%program10.4－6T=4;w0=2*pi/T;t=-6:0.01:6;N=51;c0=0.5;xN=c0*ones(1,length(t));%dccomponentRC=0.1;forn=1:2:N%evenharmonicsarezeroH=abs(1/(1+j*RC*w0*n));phi=angle(1/(1+j*RC*w0*n));第10章MATLAB在信号与系统分析中的应用xN=xN+H*cos(w0*n*t+phi)*sinc(n*0.5);endplot(t,xN);title(［′RC=′,num2str(RC)］);xlabel(′t(s)′);ylabel(′y(t)′);第10章MATLAB在信号与系统分析中的应用图10.4－11例10.4－6程序运行结果第10章MATLAB在信号与系统分析中的应用用MATLAB实现连续系统的S域分析MATLAB实现部分分式展开用MATLAB函数residue可以得到复杂S域表示F(s)的部分分式展开式,其调用形式为［r,p,k］=residue(num,den)其中,num、den分别为F(s)分子多项式和分母多项式的系数向量;r为部分分式的系数;p为极点;k为多项式的系数,若F(s)为真分式,则k为零。第10章MATLAB在信号与系统分析中的应用例10.5－1

用部分分式展开法求F(s)的反变换:解

程序如下:%program10.5－1formatratnum=［12］;den=［1430］;%结果数据以分数的形式输出［r,p］=residue(num,den)运行结果:r=-1/6-1/22/3p=-3-10第10章MATLAB在信号与系统分析中的应用因此F(s)可展开为所以有时F(s)表达式中分子多项式B(s)和分母多项式A(s)是以因子相乘的情况出现时,这时可用conv函数将因子相乘的形式转换成多项式的形式:C=conv(A,B)

其中,A和B是两因子多项式的系数向量,C是因子相乘所得多项式的系数向量。如果已知多项式的根,则可以利用poly函数将根式转换成多项式,其调用形式为B=poly(A)

式中,A为多项式的根,B为多项式的系数向量。第10章MATLAB在信号与系统分析中的应用例10.5－2

用部分分式展开法求F(s)的反变换解

程序如下:%program10.5－2num=［2305］;den=conv(［11］,［112］);［r,p,k］=residue(num,den)运行结果为r=-2.000+1.1339i-2.000-1.1339i3.000p=-0.5000+1.3229i-0.5000-1.3229i-1.000k=2第10章MATLAB在信号与系统分析中的应用由于上述留数r中有一对共轭复数,因此求时域表达式的计算较复杂。为了得到简洁的时域表达式,可以用cart2pol函数把共轭复数表示成模和相角形式,其调用形式为［TH,R］=cart2pol(X,Y)表示将笛卡儿坐标转换成极坐标,X、Y分别为笛卡儿坐标的横坐标和纵坐标。TH是极坐标的相角,单位为弧度;R是极坐标的模。第10章MATLAB在信号与系统分析中的应用因此在上述程序中增加下面语句,即可得到留数r的极坐标形式:［angle,mag］=cart2pol(real(r),imag(r))运行结果为-2.62582.299103.0000angle=2.6258mag=2.2991由此可得所以第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用例10.5－3

已知系统函数为求系统的零极点并画出零极点分布图。解

程序如下：%program10.5－3b=［1-1］;a=［143］;zs=roots(b);ps=roots(a);plot(real(zs),imag(zs),′o′,real(ps),imag(ps),′rx′,′axis(［-42-22］);grid;legend(′零点′,′极点′);第10章MATLAB在信号与系统分析中的应用图10.5－1例10.5－3程序运行结果第10章MATLAB在信号与系统分析中的应用MATLAB中提供了一种更简便的画出系统函数零极点

图的方法,即直接应用pzmap函数画图。pzmap函数调用形式为pzmap(sys)表示画出sys所描述系统的零极点图。LTI系统模型sys要借助tf函数获得,其调用方式为sys=tf(b,a)式中,b和a分别为系统函数H(s)分子多项式和分母多项式的系数向量。因此,上例还可用下述程序实现:第10章MATLAB在信号与系统分析中的应用%program10.5－3b=［1-1］;a=［143］;sys=tf(b,a)pzmap(sys)得到的零极点分布图如图10.5－2所示。如果已知系统函数H(s),求系统的单位冲激响应h(t)和频率响应H(jω)可以用impulse函数和freqs函数。第10章MATLAB在信号与系统分析中的应用图10.5－2零极点分布图第10章MATLAB在信号与系统分析中的应用例10.5－4

已知系统函数为试画出系统的零极点分布图,求系统的单位冲激响应h(t)和频率响应H(jω),并判断系统是否稳定。第10章MATLAB在信号与系统分析中的应用解

程序如下:%program10.5－4num=［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′)第10章MATLAB在信号与系统分析中的应用运行结果为:poles=-1.0000

-0.5000

+0.8660i

-0.5000

-0.8660i图10.5－3为例10.5－4程序的运行结果。三个极点均位于S平面左半开平面上,故该系统是稳定系统。系统函数的零极点分布图、系统的单位冲激响应和频率响应分别如图所示。第10章MATLAB在信号与系统分析中的应用图10.5－3例10.5－4程序第10章MATLAB在信号与系统分析中的应用10.5.3 用MATLAB计算拉普拉斯变换MATLAB的符号数学工具箱提供了计算拉普拉斯正反变换的函数laplace和ilaplace,其调用形式为F=laplace(f)f=ilaplace(F)式中f为信号的时域表达式的符号对象,F表示信号f的象函数表达式的符号对象。符号对象可以应用函数sym实现,其调用格式为S＝sym(A)式中,A为待分析表示式的字符串;S为符号数字或变量。第10章MATLAB在信号与系统分析中的应用例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)第10章MATLAB在信号与系统分析中的应用(2)程序如下：%program10.5－5(2)F=sym(′s^2/(s^2+1)′);ft=ilaplace(F)运行结果为ft=Dirac(t)-sin(t)第10章MATLAB在信号与系统分析中的应用10.6 用MATLAB实现离散系统的Z域分析10.6.1 部分分式展开的MATLAB实现信号的Z域表示式通常可用下面的有理分式表示为了能从信号的Z域象函数方便地得到其时域原函数,可以将F(z)展开成部分分式之和的形式,再对其取Z逆变换。第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用例10.6－1

试用MATLAB计算的部分分式展开。解

计算部分分式展开的[WTBZ]MATLAB程序如下:%program10.6－1num=［18］;den=［183-4-1］;［r,p,k］=residuez(num,den)程序运行结果为r=0.36000.24000.4000p=0.5000-0.3333-0.3333k=［］第10章MATLAB在信号与系统分析中的应用从运行结果中可以看出p(2)=p(3),表示系统有一个二阶的重极点,r(2)表示一阶极点前的系数,而r(3)就表示二阶极点前的系数。对高阶重极点,其表示方法是完全类似的,所以该F(z)的部分分式展开为第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用例10.6－2

已知一离散因果LTI系统的系统函数为求该系统的零极点。解

将系统函数改写为第10章MATLAB在信号与系统分析中的应用用tf2zp函数求系统的零极点,程序如下:%program10.6－2b=［121］;a=［1-0.5-0.0050.3］;［r,p,k］=tf2zp(b,a)

运行结果为r=-1

-1p=0.5198+0.5346ik=10.5198-0.5346i-0.5396第10章MATLAB在信号与系统分析中的应用若要获得系统函数H(z)的零极点分布图,可以直接应用zplane函数,其调用形式为zplane(b,a)

式中,b和a分别为H(z)的分子多项式和分母多项式的系数向量。它的作用是在Z平面画出单位圆、零点和极点。如果已知系统函数H(z),求系统的单位脉冲响应h［k］和频率响应H(ejΩ),则可应用impz函数和freqz函数。第10章MATLAB在信号与系统分析中的应用例10.6－3

已知一离散因果LTI系统的系统函数为试画出系统的零极点分布图,求系统的单位脉冲响应h［k］和频率响应H(ejΩ),并判断系统是否稳定。解

根据已知的H(z),用zplane函数即可画出系统的零极点分布图。利用impz函数和freqz函数求系统的单位脉冲响应和频率响应时,需要将H(z)改写成第10章MATLAB在信号与系统分析中的应用程序如下:%program10.6－3b=［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章MATLAB在信号与系统分析中的应用图10.6－1例10.6－3运行结果第10章MATLAB在信号与系统分析中的应用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章MATLAB在信号与系统分析中的应用例10.6－4

试分别用ztrans函数和iztrans函数求:(1)f［k］=cos(ak)ε(k)的Z变换;(2