matlab在信号与图像处理中的应用第5章.docx

1、matlab在信号与图像处理中的应用第5章第5章 MATLAB信号处理基础知识在前面的章节中，我们对MATLAB R2007a的基本知识作了详细的介绍，目的是让读者对MATLAB的应用环境以及基本功能有个初步的了解。随着MATLAB版本不断升级，其功能日益强大，应用范围更是日益广泛，覆盖了工业、电子、医疗、建筑等众多领域。特别是MATLAB主要用于矩阵运算，在数据分析、信号处理、自动控制等方面具有非常强大的功能，已成为科研人员和工程技术人员必不可少的工具之一。从本章开始，我们将介绍MATLAB在特定领域内的应用。数字信号处理是当前的热门学科之一，相应的MATLAB R2007a所提供的信号处理

2、工具箱则对该方面的研究提供了非常全面的支持。本章将基于信号处理工具箱，介绍MATLAB信号处理的基础知识，包括信号处理工具箱的简介、信号表示方法、数据读入方式、一些典型的离散序列和常用波形如何产生等内容。学习目标 For personal use only in study and research; not for commercial use 熟悉信号处理工具箱的基本功能和特性 掌握信号的表示方法和数据的读入方式 掌握几种典型的离散序列 掌握常用波形发生器练习案例 分别生成单位抽样序列和。 分别生成单位阶跃序列和。 生成正弦型序列。 分别生成实指数序列和。 生成复指数序列。 产生均值为0.

3、5，方差为0.1的白噪声序列。 将实指数序列扩展为具有3个周期的序列。 生成一个方波信号，要求正信号所占百分比为60%。 生成0.5s的锯齿波和三角波信号，频率都为10Hz，采样频率为10kHz。 已知输入变量x，生成sinc函数波形。 生成一个线性扫频信号。 产生一个二次型扫频信号，要求二次型扫频信号频谱为凹状。 绘出狄立克莱（Dirichlet）函数图形。 绘出中心频率为50kHz的高斯正弦脉冲，要求相对带宽60%，包络下降到峰值的40dB之下。 分别产生一个对称的和一个非对称的三角脉冲。 产生一个矩形脉冲，要求脉冲宽度为2s。 产生一个锯齿形脉冲序列，要求重复频率为3Hz，锯齿宽度0.1

4、s，信号总长度1s，采样频率1kHz。 产生一个幅度递减的高斯脉冲序列，每个脉冲的中心频率10kHz，带宽40%，脉冲重复频率为1kHz，采样率50kHz，脉冲序列长度为10ms，脉冲幅度每次递减0.8。 假设一个输入信号为三角波，绘出该信号经压控振荡器的输出波形，已知采样频率10kHz。5.1 信号处理工具箱简介MATLAB工具箱中包含了许多用于解决具体问题的应用程序专用M文件，而信号处理工具箱则包含了许多执行信号处理算法的函数，如滤波器设计与实现、频谱分析、加窗、转换等等。本节将简要介绍一下信号处理工具箱的基本情况。5.1.1 什么是信号处理工具箱信号处理工具箱（Signal Proces

5、sing Toolbox）是基于MATLAB数值计算环境的一系列工具（函数）的集合。工具箱支持各种形式的信号处理操作，从波形产生到滤波器设计和实现、参数建模和谱分析等等。工具箱提供了两大类工具：命令行函数和图形用户界面（GUI），其中命令行函数主要应用于以下几个方面： 离散时间滤波器设计、分析和实现 模拟滤波器设计、分析和实现 线性系统变换 窗函数 谱分析和倒谱分析（cepstral analysis） 变换（transforms） 统计信号处理 参数建模 线性预测 多速率信号处理 波形产生而交互式的图形用户界面主要应用于： 滤波器设计和分析 窗函数设计和分析 信号作图和分析 谱分析 滤波处理

6、5.1.2 信号的表示方法在MATLAB环境中，大部分数据都是以数值阵列的形式表示，即将一组实数或虚数按一定顺序排列在两维或更多维空间内。因而采集到的基本信号（包括一维信号或序列、多通道信号、二维信号等）都要表示成阵列的形式。对于一维采样信号或序列，在MATLAB中用向量表示。所谓向量是1n或n1的阵列，这里n是序列的采样值个数。引入一个序列的方法之一是在命令提示符后输入一列元素。例如：x = 5 3 8 -6 0 7这条语句产生了一个简单的行向量，该向量由6个实数组成的序列构成。转置变换就会将该序列变成一个列向量： x = x结果为：x =538-607对于单通道信号而言，最好采用列向量进行

7、表示，这是因为列向量较易扩展到多通道。对于多通道数据而言，一个矩阵中的每一列都对应于一个通道，而矩阵中的每一行对应于一组采样点。一个包含x、2x和x/2的三通道信号可以表示为y = x 2*x x/2将上面的x值代入，得到结果：y =5.0000 10.0000 2.50003.0000 6.0000 1.50008.0000 16.0000 4.0000-6.0000 -12.0000 -3.00000.0000 0.0000 0.00007.0000 14.0000 3.50005.1.3 数据的读入方式 在MATLAB工作环境内读入在后续章节当中的一些具体范例，数据的获取一般通过以下两种

8、方式： 直接输入，即在键盘上手动输入数据； 利用MATLAB或工具箱函数，例如sin、cos、sawtooth、square等等； 在MATLAB工作环境外读入在某些应用场合，可能需要从MATLAB工作环境之外输入。根据数据的格式，可以采用以下几种方法读入： 利用MATLAB的“load”命令，从ASCII文件或MAT文件中加载数据，该函数的使用方法已在第二章中给出介绍，这里不再赘述； 利用一个低层文件I/O函数将数据读入MATLAB，例如“fopen，fread，fscanf”等； 开发一个MEX文件来读取数据； 将数据转换为MAT文件也可利用其他资源读入数据，例如高级程序语言（Fortra

9、n或者C等）可将数据写成MAT文件形式，再用MATLAB“load”命令读取。提示：信号处理工具箱仅支持双精度输入。如果输入单精度浮点型或整数型数据，将得不到正确的结果。将滤波器设计工具箱与定点工具箱联合使用，可使单精度浮点和定点类型数据能够用于滤波处理及滤波器设计。5.1.4 工具箱的核心功能信号处理工具箱函数大部分是一些以M文件表示的算法，能够实现多种信号处理功能。这些函数是在MATLAB计算与制图环境之外的专门用于信号处理领域的功能扩展。 信号与系统工具箱函数所处理的基本实体是信号和系统。这些函数专门用于处理数字（或离散）信号，而非模拟（或连续）信号。工具箱支持的主要滤波器类型是单输入单

10、输出的线性非时变滤波器。可以利用一个或几个模型（例如传递函数、状态空间、零极点增益等）来表示线性时不变系统，当然，这些不同表达形式之间可以进行相互转换。 滤波器设计分析与实现信号处理工具箱提供滤波器设计功能，即按照特定需求，使用者可以自行设计滤波器。工具箱中包含的滤波器设计功能主要包括FIR和IIR滤波器设计分析和实现、滤波阶数估计、模拟滤波器原型设计与转换等。 线性系统转换工具箱拥有大量的转换函数，包括二次剖面、状态空间、零极点、网格或梯形函数、以及传递函数之间的转换等。 窗函数工具箱提供了许多常用的窗函数，同时也提供了图形用户界面（GUI）以便于查看和比较这些窗函数，以及利用这些窗函数设计

11、滤波器。 谱分析采用参量和非参量方法，工具箱函数可用于估计功率谱密度、均方谱和伪谱。工具箱中包含的一些谱分析方法有Burg法、协方差法、特征向量法、周期图法、Welch和Yule-Walker法等等。其他的函数可用于计算功率谱密度平均能量和单边谱，以及将DC分量移动到谱中心位置等功能。 变换函数工具箱包含各种类型的变换和逆变换函数，包括傅立叶变换、Z变换、离散余弦变换、希尔伯特变换等二次剖面、状态空间、零极点、网格或梯形函数、以及传递函数之间的转换等。 统计信号处理工具箱包含计算相关、互相关、协方差和自相关等统计信息的函数。 参数建模工具箱提供了一些用于自回归参数建模的方法，有：Burg法、协

12、方差法、Yule-Walker法等，同时也提供了一些函数用于拟合模拟滤波器或离散时间滤波器的频率响应。 线性预测工具箱包含了一些函数用于计算线性预测系数，以及自相关函数与预测多项式、反射系数和线谱频率之间的转换。 多速率信号处理抽取、增降采样、重采样和三次样条插值等多速率信号处理功能，在工具箱中也有相应函数来实现。 波形产生多种周期和非周期波形也可用工具箱中的相应函数来产生，包括线性调频脉冲、高斯射频脉冲、高斯单脉冲、脉冲串、矩形波、三角波、锯齿波、方波等等。 其他功能工具箱函数除了能够实现以上各种功能之外，还包括倒谱分析、调制、解调，以及各种绘图方法。5.2 典型离散序列在实际中的信号往往是

13、连续的，经过采样和量化之后就变成离散的序列，才可送入计算机或数字滤波器进行处理，因此，我们这里关注的都是离散序列，或者说MATLAB处理的都是离散序列。首先给出一些常用的、比较典型的离散序列信号形式。5.2.1 单位抽样序列单位抽样序列的定义为单位抽样序列在离散信号处理中的作用类似于连续信号处理中的冲激函数，但不同的是是一个极限概念的信号， 在时脉宽趋于零，幅值趋于无穷大，因而并非一个现实的信号，而在时幅值为1，形式简单便于计算。 如果将在时间轴上延迟了个采样周期，则相应表达式为在MATLAB中可借助函数zeros(1,N)来实现单位抽样序列。函数zeros(1,N)产生一个长度为N的向量，向

14、量中的每个元素都为零。 下面通过具体例子看一下单位抽样序列的MATLAB生成方法。例5-1分别生成单位抽样序列和。% 实现单位抽样序列程序k=-10:10;delta=zeros(1,10),1,zeros(1,10);stem(k,delta)% 实现单位抽样序列程序k=-10:10;delta=zeros(1,15),1,zeros(1,5);stem(k,delta)程序运行结果如图5-1所示。 (a) (b) 图5-1 单位抽样序列波形 5.2.2 单位阶跃序列单位阶跃序列的定义为单位抽样序列与单位阶跃序列关系密切，可互相表示， 表达式为在MATLAB中可利用函数ones(1, N)来

15、实现单位阶跃序列。函数ones(1,N)产生长度为N的向量，向量中的每个元素都为单位1。MATLAB生成单位阶跃序列的方法通过下面的例子给出。例5-2 分别生成单位阶跃序列和。% 实现单位阶跃序列k=-10:10;delta=zeros(1,10),ones(1,11);stem(k,delta)% 实现单位阶跃序列k=-10:10;delta=zeros(1,15),ones(1,6);stem(k,delta)程序运行结果如图5-2所示。(a) (b) 图5-2 单位阶跃序列波形 提示：MATLAB产生的都是有限长序列，因此这里的单位阶跃序列为有限区间上的u(n)。5.2.3 正弦型序列正

16、弦序列的定义如下这里为幅度，为数字角频率，为初始相位。正弦序列可对连续正弦波信号采样后得到。在MATLAB中，利用函数sin或cos产生正（余）弦序列，这与连续正弦信号类似，只不过绘图采用stem函数，而非plot函数。下面是产生正弦序列的MATLAB程序范例。例5-3生成正弦型序列。% 实现正弦型序列程序n=0:50;x=sin(pi/8*n);stem(n,x)程序运行结果如图5-3所示。图5-3 正弦序列波形5.2.4 实指数序列实指数序列的定义为这里为一个实数，的值直接影响序列的波形。当时，序列收敛；当时，序列发散。在MATLAB中通过数组的幂运算来产生实指数序列。下面是产生该序列的一

17、个具体实例。例5-4分别生成实指数序列和。% 实现实指数序列程序n=0:15;x1=(3/4).n;x2=(4/3).n;figure(1);stem(n,x1)figure(2);stem(n,x2)程序运行结果如图5-4所示。 (a) (b) 图5-4 实指数序列波形 5.2.5 复指数序列复指数序列的定义式为也可以将实部和虚部分开来表示，下面通过一个例子说明利用MATLAB如何产生一个复指数序列。例5-5产生复指数序列% 实现复指数序列程序n=0:30;x=exp(.05+i*pi/6).n;xr=real(x);xi=imag(x);xm=abs(x);xa=angle(x);figu

18、re;subplot(221);stem(n,xr);title(实部);subplot(222);stem(n,xi);title(虚部);subplot(223);stem(n,xm);title(模);subplot(224);stem(n,xa);title(相角);程序运行结果如图5-5所示。从图中可以看出，复指数序列的实部和虚部都是幅度按指数增长的正弦序列。图5-5 复指数序列波形5.2.6 随机序列随机序列在MATLAB仿真分析中应用较多，常常用来模拟背景噪声。有两个函数可直接产生两类随机序列： rand(1,N)产生0，1上均匀分布的随机序列； randn(1,N)产生均值为0

19、，方差为1的高斯随机序列，即白噪声序列。当然，其他分布的随机序列也可通过以上两个函数的相应变换产生。提示：均值和方差为特定值的白噪声序列也可通过函数randn(1,N)来产生。例5-6 产生均值为0.5，方差为0.1的白噪声序列。%生成白噪声序列程序N=10;mean_v=0.5;var_v=0.1;x=mean_v+sqrt(var_v)*randn(1,N);disp(x)程序运行结果为：x =0.4410 0.7295 0.3140 1.1904 0.4569 0.5360 0.8373 0.5187 0.4698 0.23685.2.7 周期序列所谓周期序列，是满足下式的序列 式中为周

20、期序列的周期，定义为满足上式的最小正整数。在MATLAB中有很多种方法可以把一个长度为的序列扩展为有个周期的序列，这里仅给其中的一种实现方法，主要语句为：np=0:K:N-1;xp=x(mod(np,N)+1);第二条语句中的mod(np,N)为求模运算，规则是将np除以N，所得的余数为求模运算的结果。下面通过一个具体实例给出一个周期序列的产生过程。例5-7将实指数序列扩展为具有3个周期的序列。% 生成周期序列的程序n=0:15;x=(3/4).n;N=16;k=3;np=0:k*N-1;xp=x(mod(np,N)+1);stem(np,xp)程序运行结果如图5-6所示，生成的周期序列每个周

21、期内都为一个实指数序列。图5-6 实指数序列扩展为周期序列的波形 以上是一些在MATLAB中较为常用的离散序列，其他形式的序列在这里就不再一一列举，在后续章节中如有用到，再给出相应介绍。5.3 波形发生器 在MATLAB信号处理工具箱中，有多个函数用来产生各种形式的信号波形，例如方波、三角波、线性调频信号、脉冲信号等等。这里将工具箱提供的各个信号发生函数汇总于一表内，见表5-1。表5-1 信号处理工具箱中的信号发生函数函数名产生信号形式函数名产生信号形式square方波信号chirp扫频信号sawtooth锯齿波和三角波信号diricDirichlet函数波形sincsinc函数波形gausp

22、uls高斯调制正弦脉冲pulstran脉冲串gmonopuls高斯单脉冲rectpuls矩形脉冲信号tripuls三角脉冲信号vco压控振荡函数波形接下来将对这些信号发生函数逐一进行介绍，通过一些具体的实例，说明各个函数的用法。5.3.1 方波函数square函数的调用格式为 x = square(t); x = square(t,duty);第一种调用格式产生周期为2，幅值为1的方波；第二种格式产生指定形式的方波，参数duty为0到100之间的任意一个数，它限定了在一个周期内信号为正所占的百分比。例5-8生成一个方波信号，要求正信号所占百分比为60%。%方波信号生成程序t=0:0.01:8*

23、pi;x=square(t,60);plot(t,x,LineWidth,2);axis(0 8*pi -1 1.5);set(gca,FontSize,14)生成的方波信号波形如图5-7所示。图5-7 方波信号波形5.3.2 三角波函数Sawtooth函数的调用格式为： x=sawtooth(t); x=sawtooth(t,width);第一种格式产生周期为2（相对于时间向量t而言）、峰值为1的锯齿波。所谓锯齿波，是指在2的整数倍位置，函数幅值为-1，而后以1/的斜率线性递增。第二种格式可以产生三角波，参数width为0到1之间的任意一个数，决定了在0到2之间波形峰值出现的位置。函数在(0

24、, 2*width)区间内从-1线性递增到1，在(2*width, 2)区间内从1线性递减到-1。因此，当width=0.5时，就可以产生一个标准的三角波，而sawtooth(t,1)则等价于第一种格式sawtooth(t)。下面通过两个具体实例说明锯齿波和三角波的实现过程。提示：锯齿波函数sawtooth(t)和三角波函数sawtooth(t,width)中的参量t并非常规的时间向量，应将它等价于2ft0，这里t0为常规的时间向量。例5-9生成0.5s的锯齿波和三角波信号，频率都为10Hz，采样频率为10kHz。%锯齿波生成程序fs=10000;t=0:1/fs:.5;f=10;x=sawt

25、ooth(2*pi*f*t);figure(1);plot(t,x)%三角波生成程序y=sawtooth(2*pi*f*t,.5);figure(2);plot(t,y)程序运行结果如图5-8所示。图5-8 锯齿波信号（左）和三角波信号（右） 5.3.3 sinc函数函数的调用格式为 y=sinc(x);即对输入变量x计算sinc函数的值，并输出与x相同尺寸的向量y。这里sinc函数定义为：这是宽度为2、幅值为1的矩形脉冲的连续傅立叶反变换：例5-10已知输入变量x，生成sinc函数波形%产生sinc函数波形x=-10:.01:10;y=sinc(x);plot(x,y)grid on生成的波

26、形如图5-9所示。图5-9 sinc函数波形5.3.4 扫频信号函数chirpchirp函数用于产生扫频余弦信号，也称作调频余弦信号，该信号在雷达和声纳信号处理中有着广泛的应用，其调用格式为： y = chirp(t,f0,t1,f1)在时间范围t内产生扫频余弦信号，在0时刻信号的瞬时频率为f0，在t1时刻信号的瞬时频率为f1，这里频率的单位都是Hz。如果没有特别指定，参数取默认值：f0=0（除了对数扫频取e-6，其他情况都取0值）, t1=1, f1=100。 y = chirp(t,f0,t1,f1，method)method用来指定扫频方法，包括线性、二次型、对数三种，其中函数默认方法为

27、线性扫频。利用三种方法得到的瞬时频率表达式分别如下： 线性形式 fi(t)=f0+Bt，B=(f1-f0)/t1； 二次型形式 fi(t)=f0+Bt2，B=(f1-f0)/t12； 对数形式 fi(t) = f0*Bt，B = (f1/f0)1/ t1 y = chirp(t,f0,t1,f1，method, phi)参数phi用来指定初始相位，默认值为0。 y = chirp(t,f0,t1,f1，quadratic, phi, shape)参数shape用来指定二次型扫频信号频谱的形状，有concave和convex两种形式，分别表示凹状和凸状。例5-11生成一个线性扫频信号。%产生线性

28、扫频信号fs=1000;t=0:1/fs:2;f0=0;t1=2;f1=400;y=chirp(t,f0,t1,f1);figure(1);plot(t,y,LineWidth,2);axis(0 0.5 -1 1);set(gca,FontSize,14)figure(2);spectrogram(y,256,250,256,fs,yaxis) 程序中spectrogram()是采用短时傅立叶变换方法作出线性扫频信号的时频分布图。所谓时频分布是分析信号在每个时间点上的频率值，即分析信号的频率随时间的变化关系。这里采用该分析方法处理扫频信号，可以更清晰地看出扫频信号频率变化的特点。 本例中产生

29、了一个初始频率为0Hz，在2s时间点频率增为400Hz的线性扫频信号。程序运行结果如图5-10所示。左图为线性扫频信号的时域图，这里只取了前0.5s的时段，可看出随着时间的增加，信号的频率逐渐增大。右图为该信号的时频分布图，可以清晰地看到信号频率由最初的0Hz逐渐增加到400Hz的变化趋势，该图很好地反映了扫频信号的时频域特性。图5-10 线性扫频信号波形（左）及其时频分布图（右）例5-12产生一个二次型扫频信号，要求二次型扫频信号频谱为凹状。%产生二次型扫频信号t=-1:1/fs:1;f0=0;t1=1;f1=400;y = chirp(t,f0,t1,f1,q);figure(1);plot(t,y);axis(0 .5 -1 1);figure(2);spectrogram(y,128,120,128,fs,yaxis)本例中f0f1，从而产生了一个凹状的二次型扫频信号，在0时刻的频率为0，在1s时的频率为400Hz。程序运行结果如图5-11所示，左侧的时域图只截取了0.5s时间段内的波形，从右图可以清晰地看到二次型扫频信号的频率分布。图5-11 二次型扫频信号波形（左）及其时频分布图（右）5.3.5 狄立克莱（Dirichlet）函数diricdiric用于计算狄立克莱（Diric