MATLAB软件的从入门到精通

1、目录前言.1第一章MATLAB基础知识1第二章MATLAB基本数值运算4第三章MATLAB的图形处理功能8第四章MATLAB的程序设计11第五章常用数字信号处理函数16第六章MATLAB在数字信号处理中的应用23实验一 常见离散信号的MATLAB产生和图形显示33实验二 离散系统的频率响应分析和零、极点分布37实验三 序列线性卷积、圆周卷积的计算及其关系的研究39实验四 利用DFT分析信号的频谱41实验五 信号时间尺度变换的研究43实验六 快速傅里叶变换及其应用47实验七 IIR滤波器的实现与应用56实验八 FIR滤波器的实现与应用6165前言数字信号处理研究数字序列信号的表示方法，并对信号进

2、行运算，以提取包含在其中的特殊信息。近年来，由于在研究及应用两方面均取得了进展，数字信号处理领域已经日趋成熟。本课程以计算机为工具，通过一定量的实验项目，以验证所学的概念和算法。由于MATLAB软件的功能十分强大，使用起来也非常方便，在工程技术中尤其是信号处理领域得到了广泛的应用，因此以MATLAB作为本计算机实验课的计算机语言工具。希望大家通过本教材的学习及上机实践，能基本掌握MATLAB程序设计知识，能利用MATLAB进行简单的数字信号处理问题，利用其提供的工具箱能进行滤波器的设计，为理论知识的实用化而奠定基础。第一章 MATLAB基础知识§1-1 MATLAB软件简介MATLA

3、B，Matrix Laboratory的缩写，是由Mathworks公司开发的一套用于科学工程计算的可视化高性能语言，具有强大的矩阵运算能力。它集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个界面友好的用户环境，在这个环境中，问题与求解都能方便地以数学的语言（主要是矩阵形式）或图形方式表示出来。与大家常用的Fortran 和C等高级语言相比，MATLAB的语法规则更简单，更贴近人的思维方式，被称为“草稿纸式的语言”。§1-2 MATLAB应用入门1 MATLAB的安装与卸载MATLAB软件在用户接口时具有较强的亲和力，其安装过程比较典型，直接运行光盘中的安装向导支撑程序SE

4、TUP.exe，按其提示一步步选择即可。MATLAB自身带有卸载程序，在其安装目录下有uninstall子目录，运行该目录下uninstall.exe的即可；也可以通过Windows系统的安装卸载程序进行卸载。 2 MATLAB的启动与退出MATLAB安装完成后，会自动在Windows桌面上生成一个MATLAB图标，它是指向安装目录下binwin32matlab.exe的链接，双击这个图标即可来到MATLAB集成环境的基本窗口；也可以在开始菜单的程序选项中选择MATLAB快捷方式；还可以在MATLAB的安装路径的bin子目录中双击可执行文件matlab.exe。MATLAB的退出与普通WIN3

5、2的程序一样，值得一提的是它有一个自身专有的快捷键Ctrl+Q。初次启动MATLAB后，将进入MATLAB默认设置下的桌面平台如图1-1所示。图1-1 MATLAB默认设置下的桌面平台3 MATLAB的桌面平台默认设置下的桌面平台包括6个窗口，分别是MATLAB主窗口、命令窗口(Command Window)、历史窗口Command History)、当前目录窗口Current Directory)、发行说明书窗口(Launch Pad)和工作间管理窗口(Workspace)。31 MATLAB主窗口MATLAB的其它几个窗口都包含在这个大的主窗口中，主窗口不能进行任何计算任务的操作，只用来进

6、行一些整体的环境参数的设置。主要包括菜单栏（File、Edit、View、Web、Window和Help共6个下拉菜单）、工具栏（10个按钮控件）等。工具栏各按钮控件及说明如下所示：32 命令窗口(Command Window)MATLAB的命令窗口如图1-2所示。其中“”为运算提示符，表示MATLAB正处在准备状态。当在提示符后输入一段运算式或命令并按Enter键后，MATLAB将给出计算结果，然后再进入准备状态。图1-2 MATLAB的命令窗口33 MATLAB常用命令MATLAB有一些嵌入函数，有时应用这些函数可以起到事半功倍的效果。MATLAB常用的控制命令见表1-1。表1-1：MAT

7、LAB常用命令命令功能cd显示或改变当前工作目录，与工具栏中同效dir列出当前目录或指定目录下的文件和子目录清单，类似于DOS命令DIRclc、home的所有显示内容，并把光标移到命令窗口的左上角clf清除MATLAB当前图形窗口中的图形clear清除内存中的变量和函数disp显示变量的内容type列出指定文件的全部内容，类似于DOS命令TYPEexit、quit退出MATLABwho列出当前工作空间中的变量whos列出当前工作空间中的变量的更多信息what列出当前目录或指定目录下的 .M文件、.MAT文件和 .MEX文件which显示指定函数或文件的路径lookfor按照指定的关键字查找所有

8、相关的 .M文件第二章 MATLAB基本数值运算§2-1 MATLAB内部特殊变量和常数MATLAB内部有很多变量和常数，用以表达特殊含义。常用的有： 变量ans：指当前未定义变量名的答案。 常数eps：表示浮点相对精度，其值是从1.0到下一个最大浮点数之间的差值。该变量值作为一些MATLAB函数计算的相对浮点相对精度，按IEEE标准，eps=2-52，近似为2.2204e-016。 常数Inf：表示无穷大。当输入或计算中有除以0时产生Inf。 虚数单位i、j：表示复数虚部单位，相当于。 NaN：表示不定型值，是由0/0运算产生的。 常数pi：表示圆周率，其值为3.141 592 6

9、53 589 7。§2-2 变量类型1变量命名规则MATLAB中对变量的命名应遵循以下规则：1） 变量名可以由字母、数字和下划线混合组成，但必须以字母开头。2） 字符长度不能大于31。3） 变量命名区分大小写。2 局部变量和全局变量局部变量是指那些每个函数体内自己定义的，不能从其它函数和MATLAB工作空间访问的变量。全局变量是指用关键字“global”声明的变量。全局变量名应尽量大写，并能反映它本身的含义。如果需要在工作空间和几个函数中都能访问一个全局变量，必须在工作空间和这几个函数中都声明该变量是全局的。§2-3 向量及其运算向量运算是矢量运算的基础，向量也是组成矩阵的

10、基本元素之一。1 向量的生成11 直接输入向量生成向量最直接的方法就是在命令窗口中直接输入。格式上的要求是，向量元素需要用“ ”括起来，元素之间可以用空格、逗号或分号分隔；需要注意的是，用空格和逗号分隔生成行向量，用分号分隔生成列向量。例 A=1,2,3 或A=1 2 3 %生成行向量 A=1;2;3 %生成列向量12 利用冒号表达式生成向量冒号表达式的基本形式为x=x0：step：xn，其中x0、step、xn分别为给定数值，x0表示向量的首元素数值，xn表示向量尾元素数值限(只有当xn- x0恰为step值的整数倍时, xn才能成为尾数)，step表示从第二个元素开始，元素数值大小与前一个

11、元素数值大小的差值。例 在命令窗口，给向量a、b、c赋值。 >> a=1:2:12a = 1 3 5 7 9 11>> b=12:-2:1b = 12 10 8 6 4 2>> c=1:2:13c = 1 3 5 7 9 11 1313 特殊向量的生成对于特殊的向量可直接调用MATLAB的函数生成。如y=linsoace(x1,x2,n)用于生成线性等分的n维行向量，使得y(1)=x1，y(n)=x2。另外，向量还可以从矩阵中提取，还可以把向量看成是1×n阶（行向量）或n×1阶（列向量）的矩阵，以矩阵形式生成。2 向量的基本运算21 加（

12、减）与数乘计算例 >> a=1,2,3,4;b=0,1,2,3;c=a-bc = 1 1 1 1>> d=a-3d = -2 -1 0 1>> 4*aans = 4 8 12 1622 对位乘、点积计算同维向量a与b的对位乘用c=a.*b实现，即c的每一个元素之值是a与b对应元素的乘积。同维向量a与b的点积，一个方法是利用函数dot来实现；另一种方法是先生成a与b的对位乘向量c， 再取c的各元素和即为a与b的点积。例 >> a.*bans = 0 2 6 12>> dot(a,b) %或sum(a.*b)ans = 20§2

13、-4 矩阵及其运算MATLAB具有强大的矩阵运算和数据处理功能，对矩阵的处理必须遵从代数规则。1 矩阵的生成（1）一般矩阵的生成对于一般的矩阵，MATLAB的生成方法有很多种。最简单的方法是从键盘直接输入矩阵元素。直接输入矩阵元素时应注意：各行的元素之间用空格或逗号隔开，行与行之间用分号或回车隔开，用中括号把矩阵所有的元素括起来。例 在工作空间产生一个3×3矩阵A可用MATLAB语言描述如下：A=1，2，3；4，5，6；7，8，9或 A= 1 2 34 5 6 7 8 9运行结果为A= 1 2 34 5 67 8 9（2）特殊矩阵的生成对于特殊的矩阵可直接调用MATLAB的函数生成用

14、函数zeros生成全0矩阵：格式为B=zeros(m,n)生成m×n的全0阵。用函数ones生成全1矩阵：格式为B=ones(m,n)生成m×n的全1阵。用函数eye生成单位阵：格式为B=eye(m,n)生成m×n矩阵，其中对角线元素全为1，其它元素为0。2矩阵的运算矩阵的运算有基本运算和函数运算两种类型。基本运算包括矩阵的加、减、乘、除、幂、转置、逆等，其主要特点是通过MATLAB提供的基本运算符+、*、/、等即可完成。函数运算主要是通过调用MATLAB系统内置的运算函数来求取矩阵的，求秩，求特征值和特征相量，等等。需要时可以参阅联机帮助和相关参考书。例 矩阵的

15、基本运算>> a=1,2,3;4,5,6;>> b=6,5,4;3,2,1;>> c=a+b %计算两个矩阵的和c =7 7 77 7 7>> d=b' %计算矩阵b的转置d = 6 35 24 1>> e=a*d %做矩阵的乘法，必须满足矩阵乘法的基本要求e =28 1073 28>> f=det(e) %求矩阵e的行列式f =54>> g=e(-1) %求矩阵e的逆g = 0.5185 -0.1852-1.3519 0.5185 第三章 MATLAB的图形处理功能 从最原始版本的MATLAB开始，图

16、形功能就已经成为基本的功能之一。随着MATLAB版本的逐步升级，MATLAB的图形工具箱从简单的点、线、面处理发展到了集二维图形、三维图形甚至四维表现图和对图形进行着色、消隐、光照处理、渲染及多视角处理等多项功能于一身的强大功能包。这里只简单讨论二维基本绘图命令及图形修饰命令。常用的绘图语句有figure、plot、subplot、stem等，图形修饰语句有title、axis、text等。1figure命令figure有两种用法，只用一句figure命令，会创建一个新的图形窗口，并返回一个整数型的窗口编号。figure(n)表示将第n号图形窗口作为当前的图形窗口，并将其显示在所有窗口的最前面

17、；如果该图形窗口不存在，则新建一个窗口，并赋以编号n。2 plot命令线性绘图函数。用法为plot(x, y,s)。参数x为横轴变量，y为纵轴变量，s用以控制图形的基本特征如颜色、粗细等的图形设置选项，通常可以省略。MATLAB语言中提供的图形控制符如表3-1所示。表3-1：MATLAB语言中的图形设置选项参数含义参数含义参数含义y黄色·点实线m紫色O圆：虚线c青色×打叉·点划线g绿色+加号-破折线b蓝色*星号向上的三角形w白色s正方形<向左的三角形k黑色d菱形>向右的三角形r红色v向下的三角形p五角星形3stem命令绘制离散序列图，常用格式为stem

18、(x)和stem(x，y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。4subplot命令subplot（m，n，i）图形显示时分割窗口命令，把一个图形窗口分为m行，n列，m×n个小窗口，并指定第i个小窗口为当前窗口。5坐标轴标注MATLAB中提供了许多关于坐标轴标注的函数，常用的函数有title、xlabel、ylabel等。其中，函数title是为图形添加标题，并将标题置于图形的顶部；而xlabel、ylabel是为x、y坐标轴添加标注，并分别将标注置于相应的坐标轴的边上。这三个函数的调用格式是大同小异的，这里仅以title为例介绍它们的调用格式

19、。title（标记）其中，标记可以文字，也可以是数学表达式。6文本标注MATLAB语言对图形进行文本注释所提供的为text。它的调用格式为text (x, y,文本)其中，(x, y)给定标注文本在图中添加的位置。7图例标注在数值计算结果的绘图中，经常会出现在同一张图形中绘制多条曲线的情况，为了能更好地区分各条曲线，MATLAB语言提供了图例标注函数legend。它能为图形中所有的曲线进行自动标注，以其输入变量作为标注文本，具体的调用格式为legend（标注1，标注2，）这里的标注1、标注2等分别对应绘图过程中按先后顺序所生成的曲线。8坐标轴的控制函数axis函数axis用来控制坐标轴的刻度范

20、围及显示形式。其最简单的调用格式为axis（xmin，xmax，ymin，ymax）即绘图时变量xxmin，xmax，变量yymin，ymax。 例 绘制信号图形。 % fun0.m 定义文件名 x=0:0.1*pi:2*pi; %定义x向量figure(1); %创建一个新的图形窗口,编号为1subplot(2,2,1); %将窗口划分为2行,2列,在第一个窗口中作图plot(x,sin(x); %画图title('正弦线'); %給图形加标题subplot(2,2,2); %在第二个窗口中作图plot(x,sin(x),'r') ; %画一正弦波红色xlab

21、el('X'); %给轴加说明ylabel('SIN(X)'); %给轴加说明subplot(2,2,3); %在第三个窗口中作图plot(x,sin(x),'-',x,cos(x),'-.m*');%画一正弦波破折线%画一余弦波紫色点划星线legend('SIN(X)','COS(X)'); %给图形加图例标注subplot(2,2,4); %在第四个窗口中作图plot(x,sin(x),'r+-'); %画一正弦波红色破折线text(4,0,'注记');以上内容以

22、文件名fun0.m存盘，在MATLAB命令窗口中执行>> fun0或将以上内容直接在MATLAB命令窗口中键入并，得到结果如下：Fig.1第四章 MATLAB的程序设计MATLAB作为一种高级计算语言，它不仅可以如前面所介绍的那样，以一种人机交互式的命令行的方式工作，还可以像BASIC、FORTRAN、C等其它高级计算机语言一样进行控制流的程序设计，即编制一种以m为扩展名的文件，简称为M文件。而且， M文件的编写具有语法简单、可读性强、调试容易、调用方便等许多优点。§4-1 M文件介绍1 M文件的特点与形式MATLAB实质上是一种解释性语言，就MATLAB本身来说，它并不

23、能做任何事情，它就像DOS操作系统的一样，本身没有实现功能而只对用户发出的指令起解释执行的作用，即命令先送到MATLAB系统内解释，再运行得到结果。因此用户可以把所要实现的指令罗列编制成文件，再统一送入MATLAB系统中解释运行，这就是以m为扩展名的M文件。用户可以使用任何字处理软件对其进行编写或修改。正是M文件的这个特点造就了MATLAB强大的可开发性和可扩展性，Mathworks公司推出的一系列工具箱就是明证。由于商用的MATLAB软件用C语言编写而成，因此，M文件的语法与C语言的十分相似。对广大的C语言爱好者来说，M文件的编写是相当容易的。M文件有两种形式，即命令式和函数式。并且要注意&

24、#183;文件扩展名一定为m。·标点符号的运用要恰到好处，建立好的书写风格，保持程序的可读性。·以符号%引导的行是注释行，不可执行，可供help命令查询；·不需要用end语句作为M文件的结束标志；·在运行此文件之前，需要把它所在目录加到MATLAB的搜索路径上去，或将文件所在目录设为当前目录。2命令式文件命令式文件就是命令行的简单叠加，MATLAB会自动按顺序执行文件中的命令，其运行相当于在命令窗口中逐行输入并运行命令，因此，用户在编制此类文件时，只需把所要执行的命令按行编辑到指定的文件中，且变量不需预先定义，也不存在文件名对应问题，也可以访问存在于整个

25、工作空间内的数据。但要注意·命令式文件在运行中所产生的所有变量均为全局变量。也就是说，这些变量一旦生成，就一直保存在内存空间中，直到用户执行clear或quit时为止。例 打开Medit窗口，编写如下程序 x=1,3,4,7;y=sum(x)/length(x);以aver.m为文件名存盘。在命令窗口键入aver并回车即可运行该程序。3函数式文件为了实现计算中的参数传递，需要用到函数式文件。函数式文件在MALTAB中应用十分广泛，MALTAB所提供的绝大多数功能函数都是由函数式文件实现的，这足以说明函数式文件的重要性。函数式文件的结构为： function 输出参数=函数名（输入参数

26、）函数体 %注释例 编制用于计算几个数的平均值的函数。打开Medit窗口，编写如下程序function y=lianxi(x)n=length(x);%确定向量x的维数if n=1 y=x;endif n>=2 y=sum(x)/n;end编写完毕后，以lianxi.m存盘。函数式文件的标志是第一行必为function语句。函数式文件可以有返回值如上例中的y，也可以只执行操作而无返回值，大多数函数式文件有返回值。函数体是函数的主体部分，它包括进行运算和赋值的所有MALTAB程序代码。函数体中可以包括流程控制、输入/输出、计算、赋值、注释以及函数调用和命令文件调用等。在函数体中完成对输出参

27、数的计算。函数式文件执行后，只保留最后结果，不保留任何中间过程，所定义的变量也仅在函数内部起作用，并随调用的结果而被清除。函数调用的过程实际上就是参数的传递过程。例 调用函数lianxi.m计算20，50，90，100，40的平均值。在命令窗口键入如下命令并回车即可。 x =20,50,90,100,40; y=lianxi(x);该调用过程把变量x传递给了输入参数，然后把函数运算的返回值传给了输出参数y。在编写函数式文件时，要特别注意函数中的·只有文件名与函数名一一对应，才能保证调用成功。·function后的语句定义函数名和输入、输出参数，在函数被调用过程中将按此输入、

28、输出格式执行。·要养成良好的注释习惯，以方便自己或其它用户的调用。§4-2 程序控制语句MALTAB的程序控制语句有循环语句、条件转移语句两种类型。1循环语句MALTAB的循环语句包括for循环和while循环两种类型。（1）for循环语法格式：for 循环变量=起始值：步长：终止值 循环体end起始值和终止值为一整型数，步长可以为整数或小数，省略步长时，默认步长为1。执行for循环时，判定循环变量的值是否大于（步长为负时则判定是否小于）终止值，不大于（步长为负时则不小于）则执行循环体，执行完毕后，再加上步长，这时循环变量的值若大于（步长为负时则小于）终止值后退出循环，否则

29、继续。for循环最大的特点是，在一般情况下，此循环语句的循环次数是预先设定好的。例 给矩阵a赋值。function a=chuzhik=5;a=zeros(k,k);for m=1:k for n=1:k a(m,n)=1/(m+n-1); endend以chuzhi.m存盘。然后在MALTAB命令窗口中执行以下命令>> a=chuzhi得到结果为：a = 1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.

30、2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111 说明：·for语句一定要有end作为结束标志，否则下面的输入都被认为是for循环之内的内容。·循环语句中的分号；可防止中间结果的输出。（2）while循环语法格式： while 表达式 循环体end其执行方式为：若表达式为真（运算值非0），则执行循环体；若表达式为假（运算结果为0），则退出循环体，执行end后的语句。while循环的特点是循环次数无法预先确定，因此在使用中要谨慎，以防止陷入死循环。例 >> a=3;b=0;while aa=a

31、-1;b=b+a;end>> bb = 3>> aa = 0同for循环比起来，while语句的判断控制可以是一个逻辑判断语句，因此它的适用范围会更广一些。例 >> a=7;b=1; c=0while (a>=1)&( b<=4) % a1与b4同时满足时执行循环体。a=a-1;b=b+2;c=3*a+2*b +c;end >> a,b,cans = 5 5 492条件转移语句MATLAB中的条件转移语句主要由if语句实现，其用法与其它高级语言相类似，基本语法格式如下：if 表达式执行语句1else执行语句2end在执行时，首

32、先计算紧跟在关键字if后的表达式，若结果为真，则执行语句；若结果为假，则执行语句2。且else子句中可以嵌套if语句而形成elseif结构，从而实现多路选择。例 欲实现一分段函数的计算。 。 function f=pdbsline(x)if x<0 f=0；elseif x<1 f=x；elseif x<2 f=2-x；else f=0；end以pdbsline.m存盘。然后在MALTAB命令窗口中执行以下命令>> pdbsline (1.36)得到结果为：ans =0. 6400第五章 常用数字信号处理函数 §5-1 波形发生器1基本信号序列的表示离散

33、时间信号是时间上不连续的“序列”，一般直接用x(n)表示。一个长度为N的序列x(n)在MATLAB中可以用一个N维行向量或列向量表示。常用的基本信号序列产生如下：1） 单位抽样序列 这一序列可用MATLAB中的zeros函数实现x=1，zeros(1,N-1)；2） 单位阶跃序列 这一序列可用MATLAB中的ones函数实现x=ones(1，N)；3） 实指数序列 MATLAB实现：n=0：N-1；x=a.n；4） 复指数序列 MATLAB实现：n=0：N-1；x=exp(c+j*w)*n)；5） 随机序列rand(1，N)产生在0,1上均匀分布的长度为N的随机序列。randn(1，N)产生长

34、度为N且均值为0方差为1的高斯随机序列，即白噪声序列。2基本波形的产生1）三角波或锯齿波发生函数：sawtooth()语法格式：sawtooth(t，width)。产生以2为周期幅值范围在-1，+1之间的三角波或锯齿波。参数t为时间向量；width是0，1之间的数，它决定于函数在一个周期内上升部分和下降部分的比例。width=0.5产生三角波，width=1产生锯齿波，此时函数可简写为sawtooth(t)。2）方波发生函数square()语法格式：square(t)。产生以2为周期幅值范围在-1，+1之间的方波，参数t为时间向量。3）sinc发生函数：sinc()语法格式：sinc(t) 例

35、 信号产生举例function fun1clear allt=0:0.0001:0.1;figure(1);x1=sawtooth(2*pi*50*t); %在0,0.1时间内产生5个周期的锯齿波subplot(2,2,1)plot(t,x1)x2=sawtooth(2*pi*50*t,0.5);%在0,0.1时间内产生5个周期的三角波subplot(2,2,2)plot(t,x2)x3=square(2*pi*50*t); %在0,0.1时间内产生5个周期的方波subplot(2,2,3)plot(t,x3)axis(0,0.1,-1.2,1.2)t=-4:0.1:4;x4=sinc(t);

36、 %产生抽样函数subplot(2,2,4)plot(t,x4)运行结果如下：§5-2 序列的操作在数字信号处理中基本的序列运算及其MATLAB实现如下表5-1所示：表5-1一些常用的序列运算及其matlab实现序列运算其MATLAB实现信号加x(n)=x1(n)+x2(n)x=x1+x2信号乘x(n)=x1(n)×x2(n)x=x1.*x2数乘y(n)=k×x(n)y=k*x折叠y(n)=x(-n)y=fliplr(x)抽样和y=sum(x(n1:n2)抽样积y=prod(x(n1:n2)信号能量Ex=sum(abs(x).2)信号功率Px=sum(abs(x)

37、.2)/n§5-3 常用数字信号处理函数对于线性离散系统,其表达方式有多种：传递函数法 零极点增益法 带余式的部分分式展开法 差分方程法 数字信号处理即研究信号通过系统后的变化，其中经常用到的函数有：1 x=roots(a)欲求多项式y(z)=a0+a1z 1+a2z 2+ 的根，可用命令x=roots(a)来实现，其中a=a0，a1，a2，是它的系数向量，x便是y(z)的根向量。例 求多项式y(z)=1-1.6z 1+0.65z 2-0.05z 3的根。编程运行如下： >> a=1,-1.6,0.65,-0.05; x=roots(a)x = 1.0000 0.5000

38、0.1000 从而原多项式可表达为y(z)=(1-z 1)(1-0.5z 1)(1-0.1z 1)。 2 b=poly(x)欲求出y(z)=(z-x1)(z-x2)的多项式表达，可用命令b=poly(x)来实现，其中x=x1,x2,是它的根向量，b是欲求多项式的系数向量。例 某滤波器的零极点增益法表达式为，求其传递函数表达式。MATLAB程序及运行结果如下：function a,b=fun2(x,y,k)b=k*poly(y); a=poly(x);>> y=0.1,0.5,1; x=0.3,2,5; k=1; b, a=fun2(y,x,k)b = 1.0000 -1.6000

39、0.6500 -0.0500a = 1.0000 -7.3000 12.1000 -3.0000从而滤波器的传递函数表达式为 。3r,p,k=residuez(b,a)由传递函数表达式求出带余式的部分分式展开式时可用命令r,p,k=residuez(b,a) 来实现，其中b、a分别是原系统传递函数表达式中的分子、分母系数向量，而r、p、k分别是该系统的带余式的部分分式展开式中的系数向量。例 对 作部分分式展开。程序和显示输出如下：>> b=-4,8; a=1,6,8; r,p,k=residuez(b,a)r = -12 8p = -4 -2k = 这就是说，系统的部分分式展开形式

40、为：。4y(n)=conv(x,h)求两个序列的线性卷积可直接采用命令y(n)=conv(x,h)来实现。其中特别要注意的是，函数conv默认序列从n=0开始。若x(n)：nx1nnx2，h(n)：nh1nnh2，则卷积结果y(n)：ny1nny2，其中ny1= nx1+ nh1，ny2= nx2+ ny2。例 已知x(n)=3，4，5，h(n)=2，6，7，8，求它们的线性卷积。>> x=3,4,5; h=2,6,7,8; n1=-1:1; n2=-1:2; y=conv(x,h) y =6 26 55 82 67 40>>nb3=n1(1)+n2(1); >&

41、gt;ne3=n1(length(x)+n2(length(h); >>n3=nb3:ne3n3 = -2 -1 0 1 2 3即序列x(n)与h(n)的线性卷积为y(n)=6，26，55，82，67，40。5X=fft(x,N) 与x=ifft(X,N)采用FFT算法计算一维序列x(n)的N点离散Fourier变换(DFT) 可以直接用命令X=fft(x,N)来实现，其逆运算用命令x =ifft(X,N)来实现。计算中若序列x（或X）的长度大于N时，截断x（或X）；否则补零。例 模拟信号x(t)=2sin(4t)+5cos(8t)，求其64点的DFT的幅值谱及相位谱。MATLAB

42、实现：function y=fun3(N)n=0:N-1;t=0.01*n;q=n*2*pi/N;x=fun30(t);y=fft(x,N);subplot(3,1,1);plot(t,x); title('原信号'); gridsubplot(3,1,2);plot(q,abs(y); title('幅值'); gridsubplot(3,1,3);plot(q,angle(y); title('相位'); gridfunction x=fun30(t)x=2*sin(4*pi*t)+5*cos(8*pi*t);>>y=fun3(6

43、4);6H,w=freqz(b,a,n,Fs)函数freqz基于FFT算法计算以传递函数形式表达的数字滤波器的频率响应。命令freqz(b,a,N,Fs)能够自动绘制以b、a为系数的数字滤波器的N点幅频和相频曲线(Fs是抽样频率，此项可缺省)，命令H,w=freqz(b,a,n)计算并返回数字滤波器的N点频率响应H（abs(H)是幅频响应，angle(H)是相频响应），相应的N个分布在(0)范围内的频率点记录在w中。例 绘制系统的幅频和相频特性曲线。>> b=1; a=1,-0.9; freqz(b,a,512) 7y=filter(b,a,x)或y,zf=filter(b,a,x

44、,xic)系统以传递函数的形式表达，b、a为其分子、分母多项式的系数向量。信号x通过该系统后的零状态响应（即该系统对信号x的滤波），可用命令y=filter(b,a,x)或y,zf= filter(b,a,x)来完成，其中zf是系统的最后状态。由于存储空间限制常常需要进行分布滤波，估算初始和最后条件对分布滤波来说是非常有用的。假设现有两部分数据，每部分有5000个点，如下：x1=randn(5000,1);x2=randn(5000,1);可能第一个序列x1对应前10分钟采集的数据，第二个序列x2对应后10分钟采集的数据。那么整个序列为x=x1;x2。如果现在没有足够的空间存储连接起来的序列x

45、。对x1、x2分别进行滤波，为了保证滤波后序列的连续性，采用滤波x1的最后条件作为滤波x2的初始条件。y1,zf=filter(b,a,x1);y2=filter(b,a,x2,zf);若计算系统在初始状态（Y,X）下的全响应，可用命令zf=filtic(b,a,Y,X) 将（Y,X）等效成输入xic，再用命令y,zf=filter(b,a,x,xic)来完成。例 求解差分方程y(n)=x(n)+x(n-1)+x(n-2)/3+0.95y(n-1)-0.9025y(n-2)，n0。其中x(n)=cos(n/3)，且y(-1)=-2，y(-2)=-3。用MATLAN语句实现如下： functio

46、n y=chafenb=1,1,1/3;a=1,-0.95,0.9025;Y=-2,-3;X=0,0;xic=filtic(b,a,Y,X); %初始状态等效成输入xicn=0:50;x=cos(pi*n/3);y=filter(b,a,x,xic); %全响应plot(n,y) %或stem(n,y)或者： function chafen %用递推法求解n=0:50; x0=cos(pi*n/3);x=0,0,x0;y=-3,-2,zeros(1,51);for i=3:53 y(i)=(x(i)+x(i-1)+x(i-2)/3+0.95*y(i-1)-0.9025*y(i-2);endz(

47、1:51)=y(3:53);plot(n,z)第六章 MATLAB在数字信号处理中的应用§6-1圆周卷积的计算对于无限长序列不能用MATLAB直接计算线性卷积，在MATLAB内部只提供了一个conv函数计算两个有限长序列的线性卷积。对于圆周卷积MATLAB内部没有提供现成的函数,我们可以按照定义式直接编程计算。例6-1 已知两序列： 求它们的线性卷积yl(n)=h(n)*x(n)和N点的圆周卷积y(n)=h(n)x(n)，并研究两者之间的关系。实现程序：（1）计算圆周卷积的函数function yc=circonv(x1,x2,N) %实现两序列x1和x2的圆周卷积if length

48、(x1)>N err0r('N must not be less than length of x1');endif length(x2)>N err0r('N must not be less than length of x2');endx1=x1,zeros(1,N-length(x1); %填充序列x1(n)使其长度为Nx2=x2,zeros(1,N-length(x2); %填充序列使x2(n)其长度为Nm=0:1:N-1;x2=x2(mod(-m,N)+1); %生成x2的圆周反转序列H=zeros(N,N);for n=1:1:N %生

49、成计算圆周卷积的矩阵H(n,:)=cirshiftd(x2,n-1,N); % x2圆周右移n位endyc=x1*H' %计算圆周卷积function y=cirshiftd(x,m,N) %序列的圆周移位if length(x)>N error('The length of x must be less than N');endx=x,zeros(1,N-length(x); 补零，长度变为Nn=0:1:N-1;y=x(mod(n-m,N)+1); %得到输出（2）研究两者之间的关系function fun5clear all;N,xn,hn=fun50;yln

50、=conv(xn,hn); %直接用函数conv计算线性卷积ycn=circonv(xn,hn,N); %用函数circonv计算N点的圆周卷积ny1=0:1:length(yln)-1;ny2=0:1:length(ycn)-1;n=0:length(xn)-1;m=0:length(hn)-1;subplot(2,2,1);stem(n,xn); xlabel('xn')subplot(2,2,2);stem(m,hn); xlabel('hn')axis(0,16,0,4);subplot(2,2,3);stem(ny1,yln); xlabel('

51、;线性卷积')subplot(2,2,4);stem(ny2,ycn); xlabel('圆周卷积')function N,xn,hn=fun50n=0:1:11;xn=0.8.n; %生成x(n)hn=ones(1,6); %生成h(n)N=12;运行结果如图所示：§6-2 利用FFT实现线性卷积若序列x1(n)、x2(n)是长度分别为N1、N2的有限长序列，N点的圆周卷积yc(n)=x1(n)x2(n)，长度为N1+N2-1的线性卷积yl(n)=x1(n)*x2(n)。由DFT的性质可知：当NN1+N2-1时有yc(n)=yl(n)=IDFTDFT x1(

52、n)·DFTx2(n)。序列较长时DFT运算通常用快速算法FFT实现。例6-2 用FFT实现上节例中的两序列的线性卷积。实现程序如下：%fun6.mn=0:1:11;m=0:1:5;N1=length(n); N2=length(m);xn=0.8.n; %生成x(n)hn=ones(1,N2); %生成h(n)N=N1+N2-1;xk=fft(xn,N);hk=fft(hn,N);yk=xk.*hk;yn=ifft(yk,N);if all(imag(xn)=0)&(all(imag(hn)=0) yn=real(yn);endx=0:N-1;stem(x,yn,'.');运行结果为：§6-3 系统的单位冲激响应线性移不变系统的特性可用它在输入信号是单位冲激序列时的输出即系统的单位抽样响应h(n)来表征。有了它，就可得到此线性移不变系统对任意输入x(n)的输出y(n)。例6-3 一个因果线性移不变系统：y(n)=0.81y(n-2)+x(n)-x(n-2)。求：（1）单位冲激响应；（2）单位阶跃响应；（3）绘制系统的幅频和相频响应。由差分方程直接得出系统函