MATLAB入门及其信号处理应用 1 概述 2 基本数值运算 3 基本语句 4 MATLAB函数 5 MATLAB在信号处理中的应用举例.

Presentation on theme: "MATLAB入门及其信号处理应用 1 概述 2 基本数值运算 3 基本语句 4 MATLAB函数 5 MATLAB在信号处理中的应用举例."— Presentation transcript:

1 MATLAB入门及其信号处理应用 1 概述 2 基本数值运算 3 基本语句 4 MATLAB函数 5 MATLAB在信号处理中的应用举例

2 1 概 述 1.1 MATLAB程序设计语言简介 matlab语言是由美国的Clever Moler博士于1980年开发的设计者的初衷是为解决“线性代数”课程的矩阵运算问题取名MATLAB即Matrix Laboratory 矩阵实验室的意思。 MATLAB已经不仅仅是一个“矩阵实验室”了，它集科学计算、信号处理、图象处理；声音处理于一身，并提供了丰富的Windows图形界面设计方法

3 对于线性系统有Ax=b 例、用一个简单命令求解线性系统 3x1+ x2 - x3 = 3.6 x1+2x2+4x3 = 2.1
A=[3 1 -1;1 2 4;-1 4 5]; b=[3.6;2.1;-1.4]; x=A\b x =1.4818 0.3848 对于线性系统有Ax=b

4 1.2 MATLAB应用入门 1. MATLAB的安装与卸载 2. MATLAB的启动与退出（ Ctrl+Q ） 1.3. MATLAB界面简介

5 Matlab 的工作界面 当前工作目录 命令 提示符 命令窗口 当前工作空间 输入命令的历史记录

6 1) 菜单栏 菜单栏中包括File、Edit、View、Web、Window和Help六个菜单项。这里着重介绍File项。  File项是数据输入/输出的接口， 包括10个子项， 这里重点介绍其中的5个子项：  New: 新建文件项。 有四个选择: MFile(.M，文本格式的MATLAB程序文件， 可以直接通过文件名的方式在MATLAB环境下解释运行); Figure(图形); Model(仿真模型文件)和GUI(可视化界面文件)。

7 Open: 打开所有MATLAB支持的文件格式，系统将自动识别并采用相应的程序对文件进行处理。例如, 打开一个
Open: 打开所有MATLAB支持的文件格式，系统将自动识别并采用相应的程序对文件进行处理。例如, 打开一个.m文件， 系统将自动打开M文件编辑器对它进行编辑。  Import Data...: 导入用于MATLAB处理的数据函数，包括各种图像文件、声音文件和.mat文件。  Save Workspace As...: 将工作空间的变量以.mat(二进制)或ASCII文本的形式存入文件。  Set Path...: 设置工作路径。可以打开路径设置(Set Path)对话框(图2)，将用户自己建立的目录加入MATLAB的目录系统中， 以便所编制的文件能够在MATLAB环境中直接调用。

8 图 2 路径设置对话框

9 单击Add Folder. 按钮可以将你的一个文件夹加入到系统路径中; Add with Subfolders
单击Add Folder... 按钮可以将你的一个文件夹加入到系统路径中; Add with Subfolders... 允许把一个文件夹包括其所有的子文件夹加入到系统路径中。这两种操作均可以直观地在右侧的路径栏内看到结果。 选中一个加入的文件夹， 你可以利用Move to Top(移至所有路径的最前面), Move Up(上移一个), Move Down(下移一个), Move to Bottom(移至所有路径的最后面)等四个按钮将改变文件在系统路径中的排列位置以利于对文件的搜索使用， 也可以利用Remove按钮将其删除。对路径操作完毕后，按Save按钮予以保存; 按Close按钮关闭本对话框; 按Revert按钮取消所有未保存的改动; 按Default按钮将还原到MATLAB安装时的路径设置； 按Help按钮则启动帮助系统解答疑难。

10 2) 命令行区 对输入命令的解释MATLAB按以下顺序进行：  ① 检查它是否是工作空间中的变量， 是则显示变量内容。 ② 检查它是否是嵌入函数， 是则运行之。 ③ 检查它是否是子函数。 ④ 检查它是否是私有函数。  ⑤ 检查它是否是位于MATLAB搜索路径范围内的函数文件或脚本文件。  请注意，如果有两个以上的方案与输入的命令相匹配， MATLAB将只执行第一个匹配。

11 4. MATLAB常用命令 表1 MATLAB常用命令

12

13 MATLAB帮助 MATLAB Help 完善的HTML格式联机帮助系统，非常全面。使用方法： 1.从help菜单中选取； 2.在命令窗口中执行helpdesk或doc。 PDF文档 用Adobe acrobat reader阅读。

14 Help系列 helpwin

15 Help ?help HELP topics: matlab\general - General purpose commands.
matlab\ops Operators and special characters. matlab\lang Programming language constructs. matlab\elmat Elementary matrices and matrix manipulation. matlab\elfun Elementary math functions. matlab\specfun Specialized math functions. matlab\matfun Matrix functions - numerical linear algebra. matlab\datafun Data analysis and Fourier transforms. … … For more help on directory/topic, type "help topic".

16 Help + 函数（类）名 ?help general General purpose commands.
MATLAB Toolbox Version Dec-1997 General information help On-line help, display text at command line. helpwin On-line help, separate window for navigation. helpdesk - Comprehensive hypertext documentation and troubleshooting. demo Run demonstrations. ver MATLAB, SIMULINK, and toolbox version information. … … See also PUNCT.

17 Lookfor ?lookfor diff SETDIFF Set difference.
DIFF Difference and approximate derivative. POLYDER Differentiate polynomial. ODE113 Solve non-stiff differential equations, variable order method. ODE15S Solve stiff differential equations, variable order method. ODE23 Solve non-stiff differential equations, low order method. ODE23S Solve stiff differential equations, low order method. … …

18 Lookfor命令搜索路径中每个M文件的第一行，扫描是否包含所要查询的字符串。
帮助机理 Help命令显示相应M文件的注释区 Lookfor命令搜索路径中每个M文件的第一行，扫描是否包含所要查询的字符串。 ?help abs ABS Absolute value. ABS(X) is the absolute value of the elements of X. When X is complex, ABS(X) is the complex modulus (magnitude) of the elements of X.

19 常用輔助指令 load, save 資料存取 (.mat檔) clc 清屏幕 clear 清除所有变量 close all 关闭所有图形
Exist 变量检验函数 What 目录中文件列表 Who 内存变量列表 Whos 内存变量详细信息 Which 确定文件位置

20 可以先输入命令的前几个字符，再按上下键缩小搜索范围
几个小技巧 Matlab 的命令记忆功能：上下箭头键 可以先输入命令的前几个字符，再按上下键缩小搜索范围 命令补全功能： Tab 键 用 Esc 键 删除命令行 低级语言包括机器语言和汇编语言。机器语言就是计算机指令的集合，它与计算机同时诞生，是第一代的计算机语言；汇编语言是用符号来表示计算机指令，被称为第二代语言 Ctrl + C 可中斷程式執行

21 联机演示系统 基本介绍 Intro

22 演示 demo

23 命令行键盘技巧  调用上一行 home 光标置于当前行开头  调用下一行 end 光标置于当前行末尾
 光标左移一个字符 esc 清除当前输入行  光标右移一个字符 del 删除光标处的字符 Ctrl＋  光标左移一个单词 backspace 删除光标前的字符 Ctrl＋  光标右移一个单词 alt＋backspace 恢复上一次的删除

24 标点 ： 具有多种应用功能 . 小数点及域访问符 ；区分行，取消运行显示等 … 续行符 ，区分列，函数参数分隔符 % 注释标记
： 具有多种应用功能 . 小数点及域访问符 ；区分行，取消运行显示等 … 续行符 ，区分列，函数参数分隔符 % 注释标记 () 指定运算先后次序 ! 调用操作系统运算 [] 矩阵定义标志 ＝ 赋值标记 {} 用于构成单元数组 ‘ 字符串标示符

25 2 基本数值运算 2.1 MATLAB内部特殊变量和常数 MATLAB内部有很多变量和常数， 用以表达特殊含义。常用的有：
(1) 变量ans： 指示当前未定义变量名的答案。  (2) 常数eps：表示浮点相对精度， 其值是从1.0到下一个最大浮点数之间的差值。该变量值作为一些MATLAB函数计算的相对浮点精度，按IEEE标准， eps=2-52,　近似为2.2204e-016。

26 (3) 常数Inf： 表示无穷大。 当输入或计算中有除以0时产生Inf。 
(4) 虚数单位i,j： 表示复数虚部单位， 相当于 。 (5) NaN： 表示不定型值， 是由 0/0 运算产生的。  (6) 常数pi： 表示圆周率π， 其值为 …。

27 2.2 变量类型 1． 变量命名规则 MATLAB中对变量的命名应遵循以下规则：  （1） 变量名可以由字母、　数字和下划线混合组成， 但必须以字母开头。 (用isvarname检验，例isvarname var1，是为1) （2）长度不超过 63 个字符（6.5 版本以前为 19 个） 。 (可用namelengthmax查看)  （3） 变量命名区分大小写。

28 2． 局部变量和全局变量 局部变量是指那些每个函数体内自己定义的，不能从其他函数和MATLAB工作空间访问的变量。  全局变量是指用关键字“global”声明的变量。 全局变量名应尽量大写，并能反映它本身的含义。如果需要在工作空间和几个函数中都能访问一个全局变量，必须在工作空间和这几个函数中都声明该变量是全局的。

29 2.3 矩阵及其运算 MATLAB具有强大的矩阵运算和数据处理功能， 对矩阵的处理必须遵从代数规则。  1． 矩阵生成 1） 一般矩阵的生成 对于一般的矩阵MATLAB的生成方法有多种。 最简单的方法是从键盘直接输入矩阵元素。直接输入矩阵元素时应注意： 各元素之间用空格或逗号隔开，用分号或回车结束矩阵行，用中括号把矩阵所有元素括起来。

30 4 5 6 7 8 9] 例1 在工作空间产生一个3×3矩阵A可用MATLAB语言描述如下： 
] 运行结果： A=  

31 2） 特殊矩阵的生成 对于特殊的矩阵可直接调用MATLAB的函数生成。  用函数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。

32 2． 矩阵的运算 矩阵的运算有基本运算和函数运算两种类型。基本运算包括矩阵的加、减、乘、除、乘方、求转置、求逆等，其主要特点是通过MATLAB提供的基本运算符＋、－、*、/（\）、^等即可完成。函数运算主要是通过调用MATLAB系统内置的运算函数来求取矩阵的行列式（det(A)）, 求秩（rank(A)）, 求特征值和特征向量（［V， D］＝eig(A)）, 求Jordan标准形（jordan(A)）和矩阵分解等。需要用时可以参阅联机帮助和相关参考书。

33 例 2 矩阵的基本运算。 A=[1, 2, 3; 4, 5, 6]; B =[6, 5, 4; 3, 2, 1];
例 2 矩阵的基本运算。 A=[1, 2, 3; 4, 5, 6]; B =[6, 5, 4; 3, 2, 1]; C =A+B %计算两个矩阵的和 D =B’ %计算矩阵B的转置 E=A*D %做矩阵乘法，必须要满足 % 矩阵乘法的基本要求 %E应该是2阶方阵 F=det(E) %求E的行列式值 G=E^(-1) %求E的逆 输出结果：  C=  D= 6 3 5 2 4 1 E= 28 10 73 28 F=54 G= 

34 3.2 绘图语句 常用的MATLAB绘图语句有figure、plot、subplot、stem等， 图形修饰语句有title、axis、text等。  1． figure figure有两种用法，只用一句figure命令，会创建一个新的图形窗口，并返回一个整数型的窗口编号。figure(n)表示将第n号图形窗口作为当前的图形窗口， 并将其显示在所有窗口的最前面; 如果该图形窗口不存在， 则新建一个窗口，并赋以编号n。

35 figure(2); plot(x,sin(x),'bo:'); 2． plot
线型绘图函数。用法plot(x,y,’s’)。参数x为横轴变量，y为纵轴变量，s控制图形的基本特征如颜色、线型等，可省略，常用方法如下所示。 b blue point solid g green o circle : dotted r red x x-mark dashdot c cyan plus dashed m magenta * star (none) no line y yellow s square k black d diamond w white v triangle (down) ^ triangle (up) < triangle (left) > triangle (right) p pentagram h hexagram x=0:0.1*pi:2*pi; %定义x向量 figure(1); %创建图形窗口，编号1 plot(x,sin(x),‘r*-'); %画图 figure(2); plot(x,sin(x),'bo:');

36 subplot(2,2,1); %在第2个窗口中作图 plot(x,sin(x),'rp'); %画一正弦波， 红色
3． Stem 绘制离散序列图，常用格式stem(y)和stem(x,y)分别和相应的plot函数的绘图规则相同，只是用stem命令绘制的是离散序列图。  figure(3); stem(x,sin(x),'b.:'); 4． subplot subplot(m,n,i) 图形显示时分割窗口命令，把一个图形窗口分为m行，n列，m×n个小窗口，并指定第i个小窗口为当前窗口。 subplot(2,2,1); %在第2个窗口中作图 plot(x,sin(x),'rp'); %画一正弦波， 红色 subplot(2,2,2); stem(x,sin(x),'b.');

37 figure plot(x,sin(x)); %画图 title('正弦线'); %给图形加标题 xlabel('X');
5． 绘图修饰命令 在绘制图形时，我们通常需要为图形添加各种注记以增加可读性。 在plot语句后使用title(′标题′)可以在图形上方添加标题， 使用xlabel(′标记′)或ylabel(′标记′)为X轴或Y轴添加说明，使用text(X值、Y值、′想加的标示′)可以在图形中任意位置添加标示。 figure plot(x,sin(x)); %画图 title('正弦线'); %给图形加标题 xlabel('X'); ylabel('SIN(X)'); text(4,0.5,'注记');

38 MATLAB 语句：  x=0:0.1*pi:2*pi; %定义x向量 figure(1); %创建一个新的图形窗口， 编号为1
subplot(2,2,1); %将窗口划分为2行, 2列, 在第1个窗口中作图 plot(x,sin(x)); %画图 title('正弦线'); %给图形加标题 subplot(2,2,2); %在第2个窗口中作图 plot(x,sin(x),'r'); %画一正弦波， 红色 xlabel('X');

39 ylabel('SIN(X)'); %给y轴加说明
subplot(2,2,3); %在第2个窗口中作图 % plot(x,sin(x),'--'); %画一正弦波, 由%变成注释 h=stem(x,sin(x), 'b.'); %绘制蓝色离散序列图 subplot(2,2,4); %在第2个窗口中作图 plot(x,sin(x),'r+'); %画一正弦波， 红色破折线 % k=stem(x,sin(x)); text(4,0.5,'注记');

40 MATLAB 语句执行结果

41 4 MATLAB函数 4.1 函数及其调用方法 在MATLAB语言中，M文件有两种形式：脚本和函数。 
脚本没有输入/输出参数，只是一些函数和命令的组合。它可在MATLAB环境下直接执行，也可以访问存在于整个工作空间内的数据。由脚本建立的变量在脚本执行完后仍将保留在工作空间中可以继续对其进行操作，直到使用clear命令对其清除为止。 函数是MATLAB语言的重要组成部分。MATLAB提供的各种工具箱中的M文件几乎都是以函数的形式给出的。函数接收输入参数，返回输出参数，且只能访问该函数本身工作空间中的变量，从命令窗或其他函数中不能对其工作空间的变量进行访问。

42 1． 函数结构 MATLAB语言中提供的函数通常由以下五个部分组成： （1） 函数定义行; 例：function ［Y,I］=max(x) （2） H1行；例： MAX Largest component. （3） 函数帮助文件； （4） 函数体； （5） 注释。

43 这五个部分中最重要的是函数定义行和函数体。 
函数定义行：MATLAB语言在M文件的第一行用关键字“function”把M文件定义为一个函数，并指定它的名字（必须和文件名相同），同时也定义了函数的输入和输出参数。 函数定义行是一个MATLAB函数所必需的，其他各部分的内容可以没有， 这种函数称为空函数。  例如： 求最大值函数“max”的定义行可描述为 function ［Y,I］=max(x)

44 其中, “max”为函数名， 输入参数为“x”， 输出参数为“Y”和“I”。 
函数体：函数体是函数的主体部分，它包括进行运算和赋值的所有MATLAB程序代码。函数体中可以包括流程控制、输入/输出、计算、赋值、注释以及函数调用和脚本文件调用等。 在函数体中完成对输出参数的计算。

45 2． 函数调用 函数调用的过程实际上就是参数传递的过程。例如，在一个脚本文件里调用函数“max”可采用如下方式：  n=1:20; a=sin(2*pi*n/20); ［Y,I］=max(a); 该调用过程把变量“a”传给了函数中的输入参数“x”,然后把函数运算的返回值传给输出参数“Y”和“I”。其中，Y是a序列的最大值，I是最大值Y对应的坐标值。

46 axis([0,1.1,-1.2,1.2]) 4.2 常用数字信号处理函数 1． 信号产生函数
1） 三角波或锯齿波发生函数： sawtooth() 语法格式：sawtooth(t,width)。产生以2π 为周期幅值范围在［-1，+1］之间的三角波或锯齿波。参数t为时间向量； width是［0，1］之间的数， 它决定函数在一个周期内上升部分和下降部分的比例。width=0.5产生三角波，width=1产生锯齿波，此时函数可简写为：sawtooth(t)。 t=0:0.001:1; x0=sawtooth(2*pi*5*t); %在［0,1］之间产生5个周期的锯齿波 figure plot(t,x0) axis([0,1.1,-1.2,1.2])

47 2）方波发生函数： square() 语法格式：square(t)。产生以2π为周期幅值范围在［-1， +1］之间的方波，参数t为时间向量。  3）sinc发生函数： sinc() 语法格式： sinc(t) t≠0 t=0

48 例 6 信号产生举例 clear all t=0:0.001:1;
例 6 信号产生举例 clear all t=0:0.001:1; x0=sawtooth(2*pi*t); %在［0， 0.1］之间产生1个周期的锯齿波 figure plot(t,x1) x1=sawtooth(2*pi*5*t); %在［0， 1］之间产生5个周期的锯齿波 subplot(221) x2=sawtooth(2*pi*5*t,0.5); %在［0，1］之间产生5个周期的三角波 subplot(222) plot(t,x2) x3=square(2*pi*5*t); %在［0， 1］之间产生5个周期的方波 subplot(223) plot(t,x3) axis(［0, 1,-1.2,1.2］)

49 (a) 锯齿波 (b) 三角波 (c) 方波 (d) 抽样函数 t=-4:0.1:4; x4=sinc(t); %产生抽样函数
subplot(224) plot(t,x4) (a) 锯齿波 (b) 三角波 (c) 方波 (d) 抽样函数

50 2．常用窗的MATLAB函数表示 表 3 常用窗的MATLAB函数表示 窗名称 MATLAB函数 矩形窗 boxcar(N) 哈明窗
hamming(N) 三角窗 triang(N) 布莱克曼窗 blackman(N) 汉宁窗 hanning(N) 凯塞-贝尔窗 kaiser(N， BETA)

51 说明：除凯塞-贝尔窗外其他窗函数的使用方法相同。 函数的参数N是窗长度，调用结果为一个列向量。 
例 产生50点的哈明窗可用MATLAB语言表示为： y=hamming(50);plot(y)凯塞-贝尔窗函数是一组可调窗函数。其语法格式为： Kaiser(N， BETA), 返回一个N点的Kaiser窗，参数BETA是窗函数表达式中的参数β，其含义参照理论部分介绍。

52 Hphase=unwrap(Hphase); %解卷绕
3． 滤波器分析与实现函数 1） 取绝对值： abs() 语法格式： abs(x)。 当x为实数时计算x的绝对值；x为复数时得到的是复数的模值；x为字符串时得到各字符的ASCII码。 2） 取相角：angle() 语法格式：angle(z)。 求复矢量或复矩阵的相角，结果为一个以弧度为单位介于-π和+π之间的值。 Hphase=angle(H); Hphase=unwrap(Hphase); %解卷绕

53 语法格式：conv(x,y)。 求矢量x和y的卷积，若x(n)和y(n)的长度分别为M和N， 则返回值是长度为M+N-1的矢量。
例 7 x(n)=［3 4 5］; h(n)=［ ］，求其线性卷积。 MATLAB语句如下：  x=［3 4 5］; h=［ ］; y=conv(x,h) 运行结果: y= 两序列的相关运算 。MATLAB实现：y=xcorr(x1,x2)。

54 4） 利用指定的数字滤波器对数据进行滤波：filter()
常用语法格式：y=filter(b,a,x)。函数filter利用数字滤波器对数据进行滤波时，采用直接Ⅱ型结构实现，因而适用于IIR和FIR两种滤波器。参数：a=［a0 a1 a2 … aM］, b=［b0 b1 b2 … bN］是滤波器系数，x为输入序列矢量，y为滤波后的输出。即： 滤波器的系统函数为：  标准形式中取a0=1，若输入滤波系数a中a0≠1时， MATLAB会自动归一化系数；若a０=0，系统给出出错信息。

55 例 8 在语音信号处理中，常利用周期脉冲信号通过AR(10)模型来近似合成浊音信号。若信号周期T=46,AR模型的系数a=［1, -1
用MATLAB语句实现如下： T=46; a= [1, ,1.2594, , , ,0.3651, ,　0.3339,　 , ］;  for i=0:5 x(i*T+1)=1;  end y=filter(1,a,x); plot(y)

56 图 5 合成浊音波形

57 5） 计算数字滤波器H(z)的频率响应H(ejω):freqz()
语法格式：［H,W］=freqz(B,A,N) 得到数字滤波器的N点的频率向量W和与之相对应的N点的频率响应向量H，计算所得的N个频率点均匀的分布在［0, π］上。参数Ａ=［a0 a1 a2 … aM］, B=［b0 b1 b2 … bN］是滤波器系数, 即滤波器H(z)形式如下：

58 参数N(与上式中阶次N的含义不同)最好选用2的整数次幂， 以便使用FFT进行快速运算，N的缺省值为512。freqz(B,A,N)直接绘制频率响应图，而不返回任何值。 
H=freqz(B,A,W) 返回W向量中指定的频率范围内的频率响应。其中，W以弧度为单位在［0, π］范围内。  ［H,F］=freqz(B,A,N,Fs) 对H(ejω)在［0,Fs/2］上等间隔采样N点，采样点频率及相应的频率响应值分别记录在F和H中。

59 [H,w]=freqz(b,a,N,'whole',Fs)
频率响应函数：freqz.m 已知A(z)、B(z), 求系统的频率响应。基本的调用格式是： [H,w]=freqz(b,a,N,'whole',Fs) N是频率轴的分点数，建议N为2的整次幂；w是返回频率轴座标向量，绘图用；Fs是抽样频率，若Fs＝1，频率轴给出归一化频率；’whole’指定计算的频率范围是从0～FS,缺省时是从0～FS/2. 幅频率响应：Hr=abs(H); Hphase=angle(H); Hphase=unwrap(Hphase);

60 h = impz(b, a, N) 或 [h,t]=impz(b,a,N) N是所需的的长度。前者绘图时n从1开始， 而后者从0开始。
impz.m 6） 计算数字滤波器H(z)的单位脉冲响应h(n):impz() 语法格式：［H,T］=impz(B,A)。 滤波器用传递函数模型限定，参数B, A分别为H(z)分子分母多项式的系数，函数返回滤波器的冲击响应列向量H和时间即采样间隔列向量T。 h = impz(b, a, N) 或 [h,t]=impz(b,a,N) N是所需的的长度。前者绘图时n从1开始， 而后者从0开始。

61 若系统的h(n)已知，由y(n)=x(n)*h(n)，用conv.m文件可求出y(n) 。 y=conv (x, h)
与逆Z变换 相关的matlab 函数 1．filter.m 本文件用来求离散系统的输出y(n) 。 若系统的h(n)已知，由y(n)=x(n)*h(n)，用conv.m文件可求出y(n) 。 y=conv (x, h) filter文件是在A(z)、B(z)已知，但不知道h(n)的情况下求y(n)的。 调用格式是: y=filter(b, a, x) x, y, a 和 b都是向量。

62 与逆Z变换 相关的matlab 函数 2.impz.m 在 A（z）、B（z）已知情况下, 求系统的单 位抽样响应 h(n)。调用格式是： h = impz(b, a, N) 或 [h,t]=impz(b,a,N) N是所需的的长度。前者绘图时n从1开始， 而后者从0开始。

63 3. residuez.m 将X(z) 的有理分式分解成简单有理分式的和，因此可用来求逆变换。调用格式： [r,p,k]= residuez(b,a) 假如知道了向量r, p和k，利用residuez.m还可反过来求出多项式A(z)、B(z)。格式是 [b,a]= residuez(r,p,k)。

64 Z变换 F=ztrans( f ) 对f(n)进行Z变换，其结果为F(z)
. 在MATLAB语言中有专门对信号进行正反Z变换的函数ztrans( ) 和itrans( )。其调用格式分别如下：   F=ztrans( f )      对f(n)进行Z变换，其结果为F(z)   F=ztrans(f,v)   对f(n)进行Z变换，其结果为F(v)   F=ztrans(f,u,v)    对f(u)进行Z变换，其结果为F(v)   f=itrans ( F )      对F(z)进行Z反变换，其结果为f(n)   f=itrans(F,u)   对F(z)进行Z反变换，其结果为f(u)   f=itrans(F,v,u )    对F(v)进行Z反变换，其结果为f(u) 注意： 在调用函数ztran( )及iztran( )之前，要用syms命令对所有需要用到的变量（如t,u,v,w）等进行说明，即要将这些变量说明成符号变量

65 Z变换 例 求数列 fn=e-n的Z变换及其逆变换。命令如下： syms n z fn=exp(-n);
Fz=ztrans(fn,n,z) %求fn的Z变换 f=iztrans(Fz,z,n) %求Fz的逆Z变换 例①．用MATLAB求出离散序列 的Z变换MATLAB程序如下： syms k z f=0.5^k;          %定义离散信号 Fz=ztrans(f)       %对离散信号进行Z变换 运行结果如下： Fz = 2*z/(2*z-1)

66 Z变换 例②．已知一离散信号的Z变换式为 , 求出它所对应的离散信号f(k). MATLAB程序如下： syms k z
Fz=2* z/(2*z-1);       %定义Z变换表达式 fk=iztrans(Fz,k)        %求反Z变换 运行结果如下：fk = (1/2)^k 例③：求序列的Z变换. syms n hn=sym('kroneckerDelta(n, 1) + kroneckerDelta(n, 2)+ kroneckerDelta(n, 3)') Hz=ztrans(hn) Hz=simplify(Hz) 运行结果如下：Fz= (z^2 + z + 1)/z^3

67 4． 变换函数 1） 一维快速离散Fourier变换： fft() 语法格式：y=fft(x)。 y是计算信号x的快速离散傅里叶变换。 当x为矩阵时，计算x中每一列信号的离散傅里叶变换。当x的长度为2的幂时，用基2算法；否则，采用较慢的分裂基算法。  y=fft(x,n)。计算n点的FFT。当x的长度大于n时，截断x; 否则补零。 2） 一维快速离散Fourier逆变换： ifft() 语法格式：y=ifft(x)。y是计算信号x的快速离散傅里叶变换的逆变换。  y=ifft(x,n)。计算n点的快速离散傅里叶变换的逆变换

68  3） 离散余弦变换（DCT）： dct() 语法格式： y=dct(x)。 计算信号x的离散余弦变换。  y=dct(x,n)。 计算n点的离散余弦变换。  当x的长度大于n时， 截断x； 否则补零。  离散余弦逆变换可由函数idct实现。

69 例 9 计算信号 , n=1, 2, …, 100 的DCT。 用MATLAB语言可实现如下：  N=100; n=1:N; x=n+20*sin(2*pi*n/20); y=dct(x); z=idct(y); subplot(311); stem(x,′.′); ylabel(′原始信号′); subplot(312); stem(y,′.′); ylabel(′DCT信号′); subplot(313); stem(z,′.′); ylabel(′IDCT信号′)

70 图 6 信号的DCT和IDCT变换

71 5 MATLAB在信号处理中的应用举例 5.1 IIR滤波器的设计与实现
基于模拟滤波器变换原理IIR滤波器的经典设计: 首先, 根据模拟滤波器的指标设计出相应的模拟滤波器; 然后, 将设计好的模拟滤波器转换成满足给定技术指标的数字滤波器。常用算法有脉冲响应不变法和双线性变换法。在MATLAB的数字信号处理工具箱中, 提供了相应的设计函数。常用的有：buttord，buttap，lp2lp，impinvar，bilinear

72 滤波器的技术指标 幅度响应指标、相位响应指标 通带要求： 阻带要求： 通带最大衰减： 阻带最小衰减： 数字低通滤波器的幅度特性

73 1． Butterworth滤波器阶数选择函数
［N,Wn］=buttord(Wp,Ws,Rp,Rs) 输入参数： Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减；  输出参数：N符合要求的滤波器最小阶数，Wn为Butterworth滤波器固有频率（3 dB）。 2．创建butterworth模拟原型滤波器  ［Z,P,K］=buttap(N);

74 3．零极点增益模型到传递函数模型的转换 ［num,den］=zp2tf(Z,P,K); 输入参数：Z,P,K分别表示零极点增益模型的零点、　极点和增益；  输出参数：num,den分别为传递函数分子和分母的多项式系数。  4．从低通原型向低通截至频率为Wn的滤波器转换 ［b,a］=lp2lp(Bap,Aap,Wn); 功能：把模拟滤波器原型转换成截至频率为Wn的低通滤波器。

75 5． 双线性变换函数 ［bz,az］=bilinear(b,a,Fs); 功能：把模拟滤波器的零极点模型转换为数字滤波器的零极点模型。其中，Fs是采样频率。 例 用双线性变换法设计一个Butterworth低通滤波器， 要求其通带截止频率100 Hz，阻带截止频率200 Hz，通带衰减Rp小于2 dB，阻带衰减大于15 dB，采样频率Fs=500 Hz。  wp=100*2*pi/Fs; ws=200*2*pi/Fs;Rp=2;Rs=15; Fs=500;Ts=1/Fs;

76 MATLAB实现程序： %把模拟滤波器原型转换成截至频率为Wn的低通滤波器 [b,a]=lp2lp(Bap,Aap,Wn); %用双线性变换法实现模拟滤波器到数字滤波器的转换 [bz,az]=bilinear(b,a,Fs); %绘制频率响应曲线 [H,W]=freqz(bz,az); plot(W*Fs/(2*pi),abs(H)); grid xlabel('频率/Hz') ylabel('频率响应幅度') %求模拟滤波器参数 Rp=2; Rs=15; Fs=500; Ts=1/Fs; wp=100*2*pi/Fs; ws=200*2*pi/Fs; Fs=1; wap=2*Fs*tan(wp/2); was=2*Fs*tan(ws/2); %选择滤波器的最小阶数 [N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %创建butterworth模拟滤波器 [Z,P,K]=buttap(N); %把滤波器零极点模型转化为传递函数模型 [Bap,Aap]=zp2tf(Z,P,K);

77 图11 Butterworth低通滤波器的频率响应

78 5.冲激响应不变法 （1） 冲激响应不变法设计IIR数字滤波器的基本原理： （2）MATLAB信号处理工箱中的专用函数impinvar( )： 格式：[BZ,AZ] =impinvar(B,A,Fs) 功能：把具有[B,A]模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的 传递函数模型[BZ,AZ]。采样频率Fs的默认值为Fs=1。

79 上例脉冲激响应不变法MATLAB实现程序：
%把模拟滤波器原型转换成截至频率为Wn的低通滤波器 [b,a]=lp2lp(Bap,Aap,Wn); %用脉冲响应不变法实现模拟滤波器到数字滤波器的转换 [bz,az]=impinvar(b,a,Fs); %绘制频率响应曲线 [H,W]=freqz(bz,az); plot(W*Fs/(2*pi),abs(H)); grid xlabel('频率/Hz') ylabel('频率响应幅度') %求模拟滤波器参数 Rp=2; Rs=15; Fs=500; Ts=1/Fs; wp=100*2*pi/Fs; ws=200*2*pi/Fs; % Fs=1; wap=wp*Fs;was=ws*Fs; %选择滤波器的最小阶数 [N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %创建butterworth模拟滤波器 [Z,P,K]=buttap(N); %把滤波器零极点模型转化为传递函数模型 [Bap,Aap]=zp2tf(Z,P,K);

80 5.冲激响应不变法 Butterworth低通滤波器的频率响应

81 5.2 FIR滤波器的设计与实现 FIR滤波器的设计方法有窗函数法和频率采样法两种，在MATLAB的数字信号处理工具箱中提供了函数fir1。fir1是采用经典窗函数法设计线性相位FIR数字滤波器，且具有标准低通、带通、高通和带阻等类型。  语法格式：  B=fir1(n,Wn) B=fir1(n,Wn, 'ftype') B=fir1(n,Wn,window) B=fir1(n,Wn, 'ftype',window)

82 其中，n为FIR滤波器的阶数，对于高通、带阻滤波器n取偶数。Wn为滤波器截止频率，取值范围0～1。对于带通、带阻滤波器，Wn=［W1,W2］,且W1<W2。 ′ftype′为滤波器类型。 缺省时为低通或带通滤波器, 为′high′时是高通滤波器, 为′stop′时是带阻滤波器。 window为窗函数，列向量，其长度为n+1; 缺省时，自动取hamming窗。输出参数B为FIR滤波器系数向量， 长度为n+1。

83 例 用窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率Wp=0. 2π，阻带截止频率Ws=0
例 用窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率Wp=0.2π，阻带截止频率Ws=0.3π，阻带衰减不小于40 dB，通带衰减不大于3 dB。 实现程序： wp=0.2*pi; ws=0.3*pi; wdelta=ws-wp; N=ceil(8*pi/wdelta); Wn=( )*pi/2; b=fir1(N,Wn/pi,hanning(N+1)); freqz(b,1,512)

84 图 10 FIR滤波器的幅频和相频特性图

85 5.3 快速离散Fourier变换函数 1） 一维快速离散Fourier变换： fft()
语法格式：y=fft(x)。 y是计算信号x的快速离散傅里叶变换。 当x为矩阵时，计算x中每一列信号的离散傅里叶变换。当x的长度为2的幂时，用基2算法；否则，采用较慢的分裂基算法。  y=fft(x,n)。计算n点的FFT。当x的长度大于n时，截断x; 否则补零。 2） 一维快速离散Fourier逆变换： ifft() 语法格式：y=ifft(x)。y是计算信号x的快速离散傅里叶变换的逆变换。  y=ifft(x,n)。计算n点的快速离散傅里叶变换的逆变换

86 5.3 利用离散傅里叶变换（DFT）分析信号的频谱
例 已知序列x(n)=2 sin(0.48πn)+cos(0.52πn) 0≤n<100, 试绘制x(n)及它的离散傅里叶变换|X(k)|图。  MATLAB实现程序： xlabel('n');ylabel('x(n)'); title('x(n) N=100'); subplot(1,2,2) k=0:length(magXK)-1; stem(k,magXK,'.'); xlabel('k');ylabel('|X(k)|'); title('X(k) N=100'); clear all N=100; n=0:N-1; xn=2*sin(0.48*pi*n)+cos(0.52*pi*n); XK=fft(xn,N); magXK=abs(XK); phaXK=angle(XK); subplot(1,2,1) plot(n,xn)

87 图8 信号及其离散傅里叶变换

88 5.4 利用FFT实现线性卷积 若序列x1(n)、x2(n)为长度分别为N1、N2的有限长序列，yc(n)=x1(n)○x2(n)，yl(n)=x1(n) * x2(n)。由DFT的性质可知：当N≥N1+N2-1时有yl(n)=yc(n)=IDFT{DFT[x1(n)]·DFT[x2(n)]}。序列较长时DFT运算通常用快速算法FFT实现。在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅里叶变换和逆变换。

89 N1=length(n); N2=length(m); xn=0.8.^n; %生成x(n)
例 用FFT实现如下两序列的线性卷积。  实现程序： 0≤n≤11 其他 0≤n≤5 clear all n=[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); %实序列的循环卷积仍然为实序列 end x=0:N-1; stem(x,yn,'.')

90 图 9 利用FFT实现线性卷积

91 5.5 fftfilt.m 用叠接相加法实现长序列卷积。格式:
y=fftfilt(h,x) 或 y=fftfilt(h, x，N) 记 的长度为 ， 的长度为 。 若采用第一个调用方式，程序自动地确定对 分段的长度 及做FFT的长度 ,显然， 最接近 ,且是的2的整次幂。分的段数为 。 采用第二个调用方式，使用者可自己指定做FFT的长度。建议使用第一个调用方式。

92 例:令x(n)为一正弦加白噪声信号，长度为500， h(n)是用fir1. m文件设计出的一个低通FIR滤波器，长度为11
例:令x(n)为一正弦加白噪声信号，长度为500， h(n)是用fir1.m文件设计出的一个低通FIR滤波器，长度为11.试用fftfilt实现长序列的卷积 clear; % 用叠接相加法，计算滤波器系数h和输入信号x的卷积 % 其中h为10阶hanning窗，x是带有高斯白噪的正弦信号 h=fir1(10,0.3,hanning(11));% h: is the impulse response N=500;p=0.05;f=1/16; % of a low-pass filter. u=randn(1,N)*sqrt(p); % u:white noise s=sin(2*pi*f*[0:N-1]); % s:sine signal x=u(1:N)+s; % x: a long sequence; y=fftfilt(h,x); % y=x*h subplot(211) plot(x); subplot(212) plot(y);

93

94 5 MATLAB在信号处理中的应用举例 5.6 线性卷积与圆周卷积的计算 例 已知两序列： 0≤n≤11 其他 0≤n≤5
例 已知两序列： 0≤n≤11 其他 0≤n≤5 求它们的线性卷积yl(n)=h(n)*x(n)和N点的圆周卷积y(n)=h(n)○N x(n), 并研究两者之间的关系。 N

95 实现程序：  计算圆周卷积的函数。 function yc=circonv(x1,x2,N) %直接计算圆周卷积 y=circonv(x1,x2,N) %输出参数： 圆周卷积结果y %输入参数： 需要计算圆周卷积的序列x1,x2和圆周卷积的点数N if length(x1)>N error('N must not be less than length of x1'); end if length(x2)>N error('N must not be less than length of x2'); %以上语句判断两个序列的长度是否小于N

96 x1=[x1,zeros(1,N-length(x1))]; %填充序列x1(n)使其长度为N1+N2-1 % (序列x1(n)的长度为N1， 序列x2(n)的长度为N2)
x2=x2(mod(-n,N)+1); %生成序列x2((-n))N H=zeros(N,N); for n=1:1:N H(n,:)=cirshiftd(x2,n-1,N); %该矩阵的k行为x2((k-1-n))N end yc=x1*H'; %计算循环卷积

97 function y=cirshiftd(x,m,N)
%directly realize circular shift for sequence x %y=cirshiftd(x,m,N); %x:input sequence whose length is less than N %m:how much to shift %N:circular length %y:output shifted sequence if length(x)>N error('the length of x must be less than N'); end x=[x,zeros(1,N-length(x))]; n=[0:1:N-1]; y=x(mod(n-m,N)+1);

98 (2)　研究两者之间的关系。 clear all; n=[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) yln=conv(xn,hn); %直接用函数conv计算线性卷积 ycn=circonv(xn,hn,N1); %用函数circonv计算N1点圆周卷积 ny1=(0:1:length(yln)-1); ny2=[0:1:length(ycn)-1]; subplot(2,1,1); stem(ny1,yln); %画图 ylabel('线性卷积') subplot(2,1,2);stem(ny2,ycn);ylabel('圆周卷积'); axis([ ])

99 图 7 线性卷积和圆周卷积的比较