MATLAB基础及实例进阶课件.ppt

《MATLAB基础及实例进阶课件.ppt》由会员分享，可在线阅读，更多相关《MATLAB基础及实例进阶课件.ppt（91页珍藏版）》请在三一办公上搜索。

1、2023/1/23,机械工业出版社,Page 1,第12章 MATLAB信号处理,【学习目标】了解信号处理的概念领会与掌握用信号处理工具箱提供的各种 方法设计滤波器理解统计信号处理的概念及应用会用信号处理GUI对信号进行分析与处理全面了解工具箱中的其他信号处理方法,2023/1/23,机械工业出版社,Page 2,第12章 MATLAB信号处理,2023/1/23,机械工业出版社,Page 3,12.1 信号处理工具箱基础,信号处理工具箱是一个基于MATLAB环境的工具集合，可以解决诸如波形产生、滤波器设计及实现、参数建模和谱分析等一大类信号处理问题。,2023/1/23,机械工业出版社,Pa

2、ge 4,1.工具箱简介,1）命令行函数MATLAB信号处理工具箱提供了命令行函数用于解决诸如数字滤波器的设计、分析及实现、模拟滤波器的设计、分析及实现、线性系统变换、窗函数、谱分析和倒谱分析、变换、统计信号处理、参数建模、线性预测、多速率信号处理和波形产生等信号处理问题。,2023/1/23,机械工业出版社,Page 5,2）图形用户界面MATLAB信号处理工具箱提供的交互式图形用户界面用于解决以下3种问题：,滤波器设计和分析窗的设计和分析信号的作图及分析、谱分析和滤波,2023/1/23,机械工业出版社,Page 6,3）支持的数据类型 信号处理工具箱仅支持双精度类型的输入数据。若输入数据

3、为单精度浮点型或单精度整型，则大多数情况下会产生错误的结果。滤波器设计工具箱与定点工具箱结合在一起，可用于单精度浮点型和定点型的滤波问题和滤波器设计问题。,2023/1/23,机械工业出版社,Page 7,2.交互式工具,由于以下直观易用的交互式工具的使用，信号处理工具箱的功能得以极大扩充。,1）滤波器设计和分析工具（FDATool）滤波器设计和分析工具（FDATool）和Filterbuilder为滤波器设计提供了一个功能完善的平台。FDATool和Filterbuilder为其他的滤波器设计方法、量化特点、C代码生成和其它增强过滤的滤波器设计工具箱产品提供了无缝连接。若配备有Filter

4、Design HDL Coder软件，则可以由FDATool和Filterbuilder 生成HDL代码。,2023/1/23,机械工业出版社,Page 8,2）滤波器可视化工具（fvtool）提供用于查看、注释和打印滤波器响应曲线的图形环境。3）信号处理工具（sptool）提供用于信号观察、滤波器设计和谱分析的丰富的图形环境。4）窗函数涉及和分析工具（wintool）提供用于设计和对比窗函数的环境。5）常函数可视化工具（wvtool）提供用于查看、注释和打印窗函数曲线的图形环境。,2023/1/23,机械工业出版社,Page 9,3.基本的信号处理概念,1）信号的表示数字阵列：MATLAB环

5、境中的中央数据结构，一个二维或多维的有序实数集合或复数集合。所涉及的基本数据类型（一维信号或序列、多通道信号和二维信号）通常都适于用阵列描述。,向量表示：MATLAB通常用1n维或n1维的向量表示一维的采样信号或序列，n是采样点数。在MATLAB中产生一个序列的一种方法是在命令窗口直接将序列元素罗列出来。,2023/1/23,机械工业出版社,Page 10,2）波形产生时间向量：假设产生信号所用的采样频率为1000 Hz，则可用下面的代码产生一个时长为1秒（间隔1毫秒）的时间信号：,t=(0:0.001:1);,产生一个由两个正弦信号（频率为50 Hz和120 Hz，幅度为1和2）构成的采样信

6、号y：,y=sin(2*pi*50*t)+2*sin(2*pi*120*t);,2023/1/23,机械工业出版社,Page 11,对信号y加离散的正态分布的白噪声，并画出其前50个点,randn(state,0);yn=y+0.5*randn(size(t);plot(t(1:50),yn(1:50),2023/1/23,机械工业出版社,Page 12,下面的代码用来产生单位脉冲信号、单位阶跃信号、单位斜坡函数和方波：,t=(0:0.001:1);imp=1;zeros(99,1);%单位脉冲信号 unit_step=ones(100,1);%单位阶跃信号 ramp_sig=t;%单位斜坡信

7、号 quad_sig=t.2;%时间信号的二次波 sq_wave=square(4*pi*t);%占空比为50%的方波,2023/1/23,机械工业出版社,Page 13,多通道信号：可由MATLAB中的矩阵描述。例如，下面代码用以产生上面的代码中最后三个信号构成的3通道信号：,z=ramp_sig quad_sig sq_wave;,多通道信号也可由下面的方法产生：先产生一个列向量，再将列向量进行复制以构成一个矩阵，该矩阵就代表一个多通道信号。例如，下面的代码先产生一个6元素的列向量（首元素为1，其他5个元素为0），再将该向量复制构成一个3通道信号：,2023/1/23,机械工业出版社,Pa

8、ge 14,a=1 zeros(1,5);c=a(:,ones(1,3)c=1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,常用的周期信号 信号处理工具箱提供了用于产生广泛应用的周期信号的函数，如：,sawtooth产生锯齿波，幅度为1，周期为2；使用该函数时，可以用参数“width”指定信号的最大值出现的位置。,square产生周期为2的方波；使用时，可以指定方波的占空比（方波信号的正值所占周期的比例）。,2023/1/23,机械工业出版社,Page 15,4.信号产生实例,1.序列的产生：分别用行向量和列向量表示序列,x=4 3 7-9 1 x=4 3 7-9 1

9、x=x%序列的转置 x=4 3 7-9 1,2023/1/23,机械工业出版社,Page 16,对于单通道信号而言，通常用列向量表示，因为列向量容易扩展到多通道状态。对于多通道信号，矩阵的每一列代表一个通道的信号，而矩阵的每一行则对应一个采样点。例如，在例12-1的基础上产生下面的3通道信号，代码如下：,y=x 2*x x/pi y=4.0000 8.0000 1.2732 3.0000 6.0000 0.9549 7.0000 14.0000 2.2282-9.0000-18.0000-2.8648 1.0000 2.0000 0.3183,2023/1/23,机械工业出版社,Page 17

10、,若信号值中有复数元素，则转置操作将对复元素取共轭。因此，若想在对一个由复数元素构成的行向量进行转置时，复元素不取共轭，可以用“.”实现。例如：,x=1-i 3+i 2+3*i 4-2*i;x1=x x1=1.0000+1.0000i 3.0000-1.0000i 2.0000-3.0000i 4.0000+2.0000i,x2=x.x2=1.0000-1.0000i 3.0000+1.0000i 2.0000+3.0000i 4.0000-2.0000i,2023/1/23,机械工业出版社,Page 18,2.周期信号的产生：产生一个时长为1.5秒、频率为50 Hz、采样频率为10 kHz的

11、锯齿波信号，并画出前0.2秒的波形。,fs=10000;t=0:1/fs:1.5;x=sawtooth(2*pi*50*t);plot(t,x);hold on;axis(0 0.2-1 1);,2023/1/23,机械工业出版社,Page 19,12.2 滤波器设计与实现,滤波器设计是根据指定的滤波性能要求确定滤波器系数的过程，滤波器的实现则是根据这些系数选择合适的滤波器结构。滤波器的设计和实现都完成以后，即可对输入数据进行滤波。本节介绍如何利用信号处理工具箱设计和实现滤波器。,2023/1/23,机械工业出版社,Page 20,1.滤波器的实现和分析,1）卷积和滤波滤波的数学基础就是卷积运

12、算。MATLAB中的函数conv可以用来实现两个一维向量的卷积运算，下面的代码可以实现序列 和其自身的卷积：,conv(1 1 1,1 1 1)ans=1 2 3 2 1,2023/1/23,机械工业出版社,Page 21,假设一个滑动平均滤波器的单位脉冲响应为输入序列x(k)为一个长度为5的随机向量，若在MATLAB中分别用变量h和x描述，则用下面的代码即可计算出滤波器的输出y：,x=randn(5,1);h=1 1 1 1/4;y=conv(h,x);,2023/1/23,机械工业出版社,Page 22,2）滤波器和传递函数 若用Y(z)和X(z)分别表示数字滤波器的输出序列y(n)和输入

13、序列x(n)的Z变换，H(z)表示滤波器的脉冲响应的Z变换，则式子,即是滤波器的传递函数。b(i)和a(i)是滤波器的系数，滤波器的阶数则由n和m中的较大者描述。,2023/1/23,机械工业出版社,Page 23,滤波器的系数和名字,n=0：b为一个标量，滤波器称为无限脉冲响应（IIR）滤波器、全极点型滤波器、递归系统或 自回归（AR）滤波器。m=0：a为一个标量，滤波器称为有限脉冲响应（FIR）滤波器、全零点型滤波器、非递归系 统或滑动平均（MA）滤波器。n和m都大于0：滤波器是一个IIR、零极点型的递 归系统或自回归滑动平均（ARMA）滤波器。,2023/1/23,机械工业出版社,Pag

14、e 24,用函数filter进行滤波,若某单极点数字滤波器（低通）的传递函数为,则该滤波器在MATLAB中可用向量b和a描述,b=1;a=1-0.9;,2023/1/23,机械工业出版社,Page 25,2.频率响应,数字滤波器的频率响应可由MATLAB中的函数freqz求解，若一个数字滤波器的频率响应用下式描述：,下面的代码求出滤波器的p点复频率响应：,h,w=freqz(b,a,p);,2023/1/23,机械工业出版社,Page 26,函数freqz还可以指定采样频率或任意的频率点数。例如：若要求一个12阶切比雪夫I型滤波器的256点频率响应，且指定采样频率为1000 Hz，实现代码如下

15、：,b,a=cheby1(12,0.5,200/500);h,f=freqz(b,a,256,1000);,若采样频率用Fs描述，则上述代码中的返回值f是一个256点的频率向量，其元素在0到Fs/2之间。,2023/1/23,机械工业出版社,Page 27,若要画出滤波器的幅度相应和相位响应，则在调用函数freqz时省去输出即可。例如，一个9阶的巴特沃斯低通滤波器的截止频率为400 Hz，采样频率为2000 Hz，求其256点的复频率响应，并画出其幅度相应和相位响应，实现代码如下：,b,a=butter(9,400/1000);freqz(b,a,256,2000),2023/1/23,机械工

16、业出版社,Page 28,2023/1/23,机械工业出版社,Page 29,观察滤波器的幅度响应和相位响应也可以经由fvtool完成。在MATLAB的命令窗口输入下面的代码：,fvtool(b,a),在打开的窗口中单击工具栏上的，即可在同一窗口中观察滤波器的幅度相应和相位响应，如下图所示。,2023/1/23,机械工业出版社,Page 30,2023/1/23,机械工业出版社,Page 31,3.滤波器的设计,1）IIR滤波器设计,2023/1/23,机械工业出版社,Page 32,下面的代码分别设计了不同类型的低通、带通、高通和带阻滤波器：,b,a=butter(5,0.4);%巴特沃斯低

17、通滤波器b,a=cheby1(4,1,0.4 0.7);%切比雪夫I型带通滤波器b,a=cheby2(6,60,0.8,high);%切比雪夫II型高通滤波器b,a=ellip(3,1,60,0.4 0.7,stop);%椭圆型带阻滤波器,2023/1/23,机械工业出版社,Page 33,传统的IIR数字滤波器设计包括下面3个步骤：1）设计截止频率为1的模拟低通滤波器（归一化模 拟低通滤波器），并将该原型滤波器转换为实 际的模拟滤波器；2）将模拟滤波器变换到数字域；3）将滤波器离散化。,在上述设计过程中所用到的函数及功能如下：,2023/1/23,机械工业出版社,Page 34,模拟低通原型

18、滤波器设计buttap：设计巴特沃斯型低通模拟原型滤波器；cheb1ap：设计切比雪夫I型低通模拟原型滤波器besselap：设计Bessel型低通模拟原型滤波器；ellipap：设计椭圆型低通模拟原型滤波器；cheb2ap：设计切比雪夫II型低通模拟原型滤波器,频率变换lp2lp：将低通模拟原型滤波器转换为模拟低通滤波器lp2hp：将低通模拟原型滤波器转换为模拟高通滤波器lp2bp：将低通模拟原型滤波器转换为模拟带通滤波器lp2bs：将低通模拟原型滤波器转换为模拟带阻滤波器,离散化bilinear：双线性变换法设计IIR数字滤波器impinvar：脉冲响应不变法设计IIR数字滤波器,2023

19、/1/23,机械工业出版社,Page 35,2）FIR滤波器设计,信号处理工具箱提供了许多用于FIR滤波器设计的函数，如：fir1、fir2、kaiserord、firls、firpm、firpmord、fircls、fircls1、cfirpm和firrcos等。在这些函数中，cfirpm可以设计任意一种线性相位滤波器和非线性相位滤波器，而其他的函数都只能设计线性相位滤波器。,2023/1/23,机械工业出版社,Page 36,以函数fir1为例，fir1是用窗化法设计FIR滤波器，其语法形式为：,b=fir1(n,Wn)b=fir1(n,Wn,win),Wn：截止频率（取值在0,1之间，1

20、对应奈奎斯 特频率）；b：长为n+1的行向量；win：长为n+1的向量，win省略则默认使用海明窗,2023/1/23,机械工业出版社,Page 37,4.滤波器设计实例,用双线性变换法设计一个IIR数字带通滤波器。,clc;clear all;close all;z,p,k=cheb1ap(5,3);b,a=zp2tf(z,p,k);fs=2;u1=2*fs*tan(0.1*(2*pi/fs)/2);u2=2*fs*tan(0.5*(2*pi/fs)/2);Bw=u2-u1;Wo=sqrt(u1*u2);bt,at=lp2bp(b,a,Wo,Bw);bz,az=bilinear(bt,at,

21、fs);fvtool(bz,az);,2023/1/23,机械工业出版社,Page 38,2023/1/23,机械工业出版社,Page 39,设计一个48阶带通滤波器，通带为,b=fir1(48,0.35 0.65);fvtool(b,1),2023/1/23,机械工业出版社,Page 40,5.用Filterbuilder GUI设计滤波器,用filterbuilder设计一个带通IIR滤波器的步骤,1）在MATLAB的命令窗口输入filterbuilder，打开如图所示的滤波器响应类型选择对话框；,2023/1/23,机械工业出版社,Page 41,2）在列表框中选择【Bandpass】，

22、然后单击，打开如图所示的带通滤波器设计界面；,2023/1/23,机械工业出版社,Page 42,3）将【Impulse response】项选为“IIR”，【Order mode】项选为“Minimum”，【Design method】项选为“Butterworth”，单击窗口右上角的，即可观察设计出的IIR滤波器的幅频响应。,2023/1/23,机械工业出版社,Page 43,6.滤波器设计和分析GUI,1）用FDATool设计滤波器 以设计FIR带通滤波器为例，步骤如下：,在MATLAB的命令窗口输入fdatool，打开滤波 器设计窗口；,指定滤波器的类型、设计方法和各种设计指标；,单击

23、窗口下方的，即在窗口中显示 滤波器的幅频响应曲线。,2023/1/23,机械工业出版社,Page 44,2023/1/23,机械工业出版社,Page 45,2）用FDATool分析滤波器,改变图像显示区域所显示的曲线种类,菜单项【Analysis】的下拉菜单用来显示滤波器的各种响应类型。其中，下拉菜单项【Overlay Analysis】可用来同时显示多种响应类型。,工具栏上的按钮组 也可用来改变图形显示区域所显示的响应类型，从左到右依次是：幅度相应、相位响应、幅度和相位响应、群时延、相时延、脉冲响应、阶跃响应、零极点图、滤波器系数、滤波器信息、幅度相应估计和舍入噪声功率谱。,2023/1/2

24、3,机械工业出版社,Page 46,显示数据点信息,在曲线的某一点上单击鼠标右键，即可显示该点的数据信息，已显示信息的点用黑色方块标注；在黑色方块上再次单击鼠标右键，即可通过右键单击菜单改变数据信息的显示位置、字体大小、是否可移动选项、插值方法、是否将信息输出到变量空间选项、删除该数据点信息选项和删除所有数据点信息选项。,2023/1/23,机械工业出版社,Page 47,改变分析参数和采样频率,用鼠标右键单击图形显示区域的空白处，显示下拉菜单【Analysis Parameters】、【Sampling Frequency】和【Whats This?】，分别用来设置分析参数、改变采样频率和显

25、示帮助信息。,在独立窗口显示响应曲线,单击菜单【View】【Filter Visualization Tool】即可打开fvtool，用独立的窗口显示滤波器的响应曲线。,2023/1/23,机械工业出版社,Page 48,滤波器的其他分析方法,在滤波器设计窗口的左边界下部有一个如下图所示的工具栏（设计窗口上是由上至下显示的），从上到下依次是：创建多速率滤波器、滤波器类型转换、设置量化参数、滤波器实现模型、零极点编辑器和从变量空间导入滤波器。,2023/1/23,机械工业出版社,Page 49,12.3 统计信号处理,统计信号处理是指，利用信号与干扰的不同统计特性，从统计的观点出发，研究如何对这

26、种随机性观测结果进行处理，以便尽可能地抑制干扰而提取有用的信号。统计信号处理是一个重要的科学分支，其主要内容是信号的检测理论、估计和滤波理论及应用。,2023/1/23,机械工业出版社,Page 50,1.相关和协方差,两个随机过程的互相关和互协方差可由函数xcorr和xcov计算，特殊地，函数xcorr和xcov也可以计算自相关和自协方差。与卷积的计算类似，相关的计算函数xcorr采用的算法思想也是FFT。不同的是，在进行相关计算前，应先将两个子序列中的一个序列进行翻转。,2023/1/23,机械工业出版社,Page 51,2.谱估计,2023/1/23,机械工业出版社,Page 52,用非

27、参数化方法中的周期图法，对叠加噪声的两个正弦信号的和xn做功率谱估计。,clc;clear all;close all;fs=1000;t=(0:fs)/fs;A=1 2;f=150;140;xn=A*sin(2*pi*f*t)+0.1*randn(size(t);Hs=spectrum.periodogram(Hamming);psd(Hs,xn,Fs,fs,NFFT,1024,SpectrumType,onesided);,2023/1/23,机械工业出版社,Page 53,2023/1/23,机械工业出版社,Page 54,12.4 信号处理GUI,在MATLAB的命令窗口输入sptoo

28、l，得右图。,2023/1/23,机械工业出版社,Page 55,1.信号观察器,信号观察器完成的工作如下：,分析和对比向量信号或阵列（矩阵）信号放大信号的一部分测量信号数据的多种特征对比多通道信号在音频硬件上播放部分信号打印信号波形,2023/1/23,机械工业出版社,Page 56,观察信号的步骤,1）选中【Signals】列表框中的一个信号（例如 chirp）；,2）单击下方的显示的波形如右图所示。,2023/1/23,机械工业出版社,Page 57,2.滤波器可视化工具,FVTool的启动是从SPTool开始的，方法如下：,单击SPTool窗口上的【Filters】列表框下方的：此时F

29、VTool显示已选择滤波器的幅度相应，且和SPTool相联系，SPTool中的任何改变将引起FVTool中相应的改变。,单击FVTool中的菜单【File】【New Filter Analysis】：此时FVTool是一个独立的窗口，与SPTool没有任何联系。,2023/1/23,机械工业出版社,Page 58,3.频谱观察器,频谱观测器可以完成下列工作：,分析和对比谱密度曲线用不同的谱估计方法（Burg、协方差、FFT、修正协方差、MTM、MUSIC、Welch和Yule-Walker AR）创建功率谱修改功率谱密度参数（如FFT长度、窗的类型和 采样频率）打印谱曲线,2023/1/23,

30、机械工业出版社,Page 59,若要查看SPTool窗口的【Spectra】列表框中已有信号的功率谱，则将该信号选中，然后单击下方的 即可。,若【Spectra】列表框中没有要查看的信号，则创建步骤如下：,1）在【Signals】列表框中选择一个信号（如 mtlb）；,2）单击【Spectra】列表框下方的,3）在打开的频谱观测器窗口中单击左下方的,2023/1/23,机械工业出版社,Page 60,2023/1/23,机械工业出版社,Page 61,4.噪声的滤波和分析,1）从SPTool中导入一个信号,在命令窗口输入下面的代码；,randn(state,0);x=randn(5000,1)

31、;,输入sptool打开SPTool；,单击菜单【File】【Import】；,在打开的Import to SPTool窗口中选择信号x，单击，将x导入到右侧的【Data】区域，在右下方的【Sampling Frequency】处输入5000，在【Name】处输入noise，单击。,2023/1/23,机械工业出版社,Page 62,2）设计滤波器,2023/1/23,机械工业出版社,Page 63,3）将设计的滤波器应用于信号noise,选择信号noise vector和滤波器fir1design；,单击【Filters】列表框下方的,打开Apply Filter窗口，如下图所示；,在【Ou

32、tput Signal】处输入blnoise；,单击,2023/1/23,机械工业出版社,Page 64,4）分析信号,对比信号noise和blnoise。,按下Shift键的同时依次用鼠标单击信号noise 和blnoise，将它们同时选中；,单击信号观察器工具栏上的“”，在打开的窗 口中选择信号blnoise，然后单击，信 号观察器中即同时显示出noise和blnoise，如 下图所示。,2023/1/23,机械工业出版社,Page 65,2023/1/23,机械工业出版社,Page 66,5）观察信号noise和blnoise的功率谱,选择信号noise，单击【Spectra】列表框下方

33、 的；,在打开的Spectrum Viewer窗口中单击，即显示出信号noise的功率谱密度；,同样的方法显示信号blnoise的功率谱密度，在 SPTool的【Spectra】列表框中显示出 spect1auto和spect2auto。,重新激活SPTool窗口，同时选中spect1auto 和spect2auto，单击【Spectra】列表框下方 的，即在图形显示区同时显示了 noise和blnoise的功率谱密度，如下图所示。,2023/1/23,机械工业出版社,Page 67,2023/1/23,机械工业出版社,Page 68,12.5 工具箱中的其他常用处理方法,除了滤波器的设计和实

34、现以及统计信号处理外，信号处理工具箱还提供了多种多样的函数，用于其它常用的信号处理理论，如：窗函数、参数建模、重采样、倒谱分析、通信中的调制与解调以及几种特殊的信号变换等。,2023/1/23,机械工业出版社,Page 69,1.窗,窗的选择在滤波器的设计和谱估计中起着非常重要的作用，主要用来抑制对无限长序列进行截短所产生的吉布斯效应的影响。信号处理工具箱提供了用于设计和分析窗的图形用户界面工具wintool和窗的实现工具wvtool。以wintool为例，在MATLAB的命令窗口输入wintool，即打开窗的设计和分析窗口，默认显示海明窗，如下图左所示。wintool可用来同时显示多个窗函数

35、的波形及幅度响应，下图右将Flat Top窗、汉宁窗、布莱克曼窗和海明窗同时显示出来以作对比。,2023/1/23,机械工业出版社,Page 70,2023/1/23,机械工业出版社,Page 71,2.参数建模,参数建模用于寻找一个数学模型的参数来描述一个信号、系统或过程，用系统的已知信息创建模型。参数建模的应用包括语音和音乐合成、数据压缩、高分辨率谱估计、通信、工业和仿真。1）可用的参数建模函数,时域建模：arburg、arcov、armcov、aryule、lpc、levinson、prony和stmcb。频域建模：invfreqz和invfreqs。,2023/1/23,机械工业出版社

36、,Page 72,2）时域建模,线性预测建模假设信号x(k)的每一个输出样本均是过的去n个输出的线性组合，即,则，信号x的一个n阶全极点模型可以用下面的代码产生：,a=lpc(x,n);,用下面的代码产生一个全极点滤波器的脉冲响应序列，并添加噪声：,randn(state,0);x=impz(1,1 0.1 0.1 0.1 0.1,10)+randn(10,1)/10;,2023/1/23,机械工业出版社,Page 73,则，系统的一个4阶全极点模型的系数可由下面的代码生成：,a=lpc(x,4)a=1.0000 0.2574 0.1666 0.1203 0.2598,3）频域建模,设计一个低

37、通巴特沃斯滤波器；,b,a=butter(4,0.4)b=0.0466 0.1863 0.2795 0.1863 0.0466a=1.0000-0.7821 0.6800-0.1827 0.0301,2023/1/23,机械工业出版社,Page 74,建模,b4,a4=invfreqz(h,w,4,4)b4=0.0466 0.1863 0.2795 0.1863 0.0466a4=1.0000-0.7821 0.6800-0.1827 0.0301,频率向量w的单位是rad/s，不一定要均匀划分。对于计算出来的w，函数invfreqz可以设计任意阶次的滤波器与之相对应。例如，下面的代码建立一个

38、3阶模型：,b4,a4=invfreqz(h,w,3,3)b4=0.0464 0.1785 0.2446 0.1276a4=1.0000-0.9502 0.7382-0.2006,2023/1/23,机械工业出版社,Page 75,函数invfreqz用下式来确定所建立的模型是最好的,其中，A(w(k)和B(w(k)是多项式系数向量a和b的傅里叶变换，n是频率向量w和脉冲响应h(k)的长度，wt(k)是不同频率之间的相对误差的加权。此时，模型的建立可由下面的代码实现：,invfreqz(h,w,n,m,wt);,但是，此模型的稳定性不能确保。,2023/1/23,机械工业出版社,Page 76

39、,函数invfreqz提供了一种使实频率响应点之间的均方误差加权和最小的算法，即,继续在MATLAB的命令窗口输入下面的代码以建立新的模型：,wt=ones(size(w);b30,a30=invfreqz(h,w,3,3,wt,30)b30=0.0464 0.1829 0.2572 0.1549a30=1.0000-0.8664 0.6630-0.1614,此时，所建立的模型总是稳定的。,2023/1/23,机械工业出版社,Page 77,将由上述两种方法所得到的模型的频率响应和原始巴特沃斯滤波器的频率响应，同时用fvtool显示出来并加以比较，如下图所示。,2023/1/23,机械工业出版

40、社,Page 78,3.重采样,信号处理工具箱提供了许多函数，用于以较高的采样率或较低的采样率对信号进行重采样，有：,upfirdn：将重采样应用于FIR滤波器spline：三次样条插值decimate：抽取interp：插值interp1：一维插值resample：以新的采样率重采样,2023/1/23,机械工业出版社,Page 79,以函数resample为例进行说明，其语法形式为：,y=resample(x,p,q);,即，函数以原始采样频率的p/q倍速率对信号x进行重采样，返回结果y的长度是x长度的p/q倍。,MATLAB中有一长度为4001的音频信号向量，采样率为7418 Hz，可在

41、命令窗口通过下面的命令查看,clear load mtlb whos Name Size Bytes Class Attributes Fs 1x1 8 double mtlb 4001x1 32008 double FsFs=7418,2023/1/23,机械工业出版社,Page 80,若要以8192 Hz的采样率将此信号播放出来，需对其进行重采样。首先，用函数rat获得重采样所需的因子p和q，重采样后进行放音。代码如下：,p,q=rat(8192/Fs,0.0001)p=127q=115 y=resample(mtlb,p,q);sound(y,8192),2023/1/23,机械工业出版

42、社,Page 81,4.倒谱分析及实例,倒谱分析可由工具箱中的函数cceps完成，用以估计输入序列的复倒谱，其语法形式为：,xhat=cceps(x),返回值xhat是一个和输入序列x同长度的实序列。,以100 Hz的采样率产生一个45 Hz的正弦波,t=0:0.01:1.27;s1=sin(2*pi*45*t);,在信号开始后的0.2秒给信号加一个回音，幅度为信号幅度的一半,s2=s1+0.5*zeros(1,20)s1(1:108);,2023/1/23,机械工业出版社,Page 82,c=cceps(s2);plot(t,c),计算新信号的复倒谱：,2023/1/23,机械工业出版社,P

43、age 83,5.通信应用及实例,1）调制,信号处理工具箱中的函数modulate可实现对信号的调制，常用的语法形式为：,y=modulate(x,fc,fs,method,opt)y,t=modulate(x,fc,fs,method,opt),其中，x是原信号，fc是载波频率，fs是采样频率，method是调制方法选择参数，opt其他的参数选项。,2023/1/23,机械工业出版社,Page 84,2）解调,解调由函数demod实现，常用的语法形式为：,x=demod(y,fc,fs,method,opt),各参数的含义同modulate。若y是一个矩阵，则demod对y按列进行解调。,调

44、制与解调应用示例。,clc;clear all;close all;t=(0:1/1000:2);x=sin(2*pi*50*t);y=modulate(x,200,1000,am);z=demod(y,200,1000,am);subplot(221),plot(t(1:150),x(1:150);title(Original Signal);subplot(212),plot(t(1:150),y(1:150);title(Modulated Signal);subplot(222),plot(t(1:150),z(1:150);title(Demodulated Signal);,202

45、3/1/23,机械工业出版社,Page 85,2023/1/23,机械工业出版社,Page 86,6.特殊变换及实例,信号处理中用到的变换主要是傅里叶变换。此外，还有其他一些变换，如：调频Z变换、离散余弦变换、希尔伯特变换（也称正交移相器）和Walsh-Hadamard变换等。信号处理工具箱提供了用于实现这些变换的函数，分别是：czt、dct（离散余弦逆变换是idct）、hilbert和fwht（反变换是ifwht）。,2023/1/23,机械工业出版社,Page 87,希尔伯特变换,clc;clear all;close all;t=(0:1/1023:1);x=sin(2*pi*60*t)

46、;y=hilbert(x);plot(t(1:50),real(y(1:50),hold onplot(t(1:50),imag(y(1:50),:);axis(0 0.05-1.1 2);legend(Real Part,Imaginary Part,location,northeast);,2023/1/23,机械工业出版社,Page 88,Walsh-Hadamard变换,clc;clear all;close all;xe=ecg(512);%产生ECG信号xr=repmat(xe,1,8);%复制信号x=xr+0.1.*randn(1,length(xr);%给信号添加噪声y=fwh

47、t(x);%对信号进行Walsh-Hadamard变换figure(color,white);subplot(2,1,1);plot(x);xlabel(Sample index);ylabel(Amplitude);title(ECG Signal);subplot(2,1,2);plot(abs(y);xlabel(Sequency index);ylabel(Magnitude);title(WHT Coefficients);,2023/1/23,机械工业出版社,Page 89,%信号的能量主要集中在1100 Hz之间的低频处%下面的代码只保存变换序列y的前1024个点，并对其作反变换%将重构信号与原信号加以比较%y(1025:length(x)=0;xHat=ifwht(y);figure(color,white);plot(x);hold on;plot(xHat,r);xlabel(Sample index);ylabel(ECG signal amplitude);legend(Original Signal,Reconstructed Signal);,2023/1/23,机械工业出版社,Page 90,2023/1/23,机械工业出版社,Page 91,