基于matlab的自动识别谱峰的程序设计 毕业论文
目录摘要1一绪论211几种常用寻峰方法的简单说明212小波变换413MATLAB小波分析工具箱6二小波分析基本原理721一维连续小波分析722一维离散小波分析823信号的初步去噪12三基于MATLAB的程序设计1431设计流程1432程序设计要点16附完整程序23结论31参考文献32致谢33摘要本文通过对染噪信号特征分析,设计了一种自动识别谱峰的程序。对比傅立叶分析的缺陷,小波方法在抑制噪音和局部分析中有着优异的性能。通过研究一维小波变换的基本原理,及其在信号去噪中的应用,基于MATLAB设计出的程序很好地完成了预期目标。关键词小波变换MATLAB去噪谱峰识别ABSTRACTBASEDONANALYSISOFTHECHARACTERISTICSOFNOISEDSIGNAL,THISPAPERISTODESIGNPROCEDURETOIDENTIFYPEAKSCONTRASTTOFOURIERANALYSIS,WAVELETSPROVIDEEXCELLENTCOMPRESSIONANDLOCALIZATIONPROPERTIESBYSTUDYINGTHETHEBASICPRINCIPLESOFONEDIMENSIONALWAVELETTRANSFORMANDTHEAPPLICATIONINSIGNALDENOISING,THEPROGRAMDESIGNEDBASEDONMATLABCOMPLETESTHETARGETSWELLKEYWORDSWAVELETTRANSFORMMATLABSIGNALDENOISINGIDENTIFYINGPEAKS第一章绪论在我们的周围,每天都有大量的信号需要我们进行分析,例如我们说话的声音,机器的振动,金融变化数据,地震信号,音乐信号,医疗图像等。相当多的信号需要进行有效的编码,压缩,消噪,重建,建模和特征提取。对于实际应用中得到的光谱信号,我们有时需要对其进行去噪、识别、定位以及寻峰等操作。本章简要说明了几种常用谱峰识别方法,然后对本文用到的小波变换和MATLAB小波工具箱作基本的说明。11几种常用寻峰方法的简单说明1)蒙特卡罗算法蒙特卡罗MONTECARLO算法,也称为统计实验方法,应用在寻峰算法中又称为质心探测法。原理为利用蒙特卡罗算法,把数据采集卡DAQ采集的波形曲线数据进行分峰截幅后,作为质量非均匀的曲线段处理。波形数据中每点的横坐标值相当于质点系中各质点的位矢,纵坐标值相当于质点系中各质点的质量大小,质点系的质心横坐标可由质心定义式的蒙特卡罗算法求出。在波形轴对称或陡峭时,质心位置与波形峰值位置一致。2)直接比较法直接比较法即线形插值微分法。原理为对数据采集卡采集的波形曲线数据进行分峰截幅后,应用一阶数值微分,微分值0点的位置即为原波形峰值位置。直接比较法应用前差公式或后差公式进行线形插值微分。两种公式的效果相同。3)二次插值数值微分法二次插值数值微分法为非线性插值数值微分法,其峰值位置判定原理和直接比较法一致。二次插值数值微分法应用中点公式,即三点公式进行非线性插值数值微分。4)一般多项式拟合法一般多项式拟合法原理为对数据采集卡采集的波形曲线数据进行分峰截幅后,采用一般多项式作拟合函数,最小二乘法作为判定,得到拟合式1拟合多项式的一阶微分解析式为2对应的一阶微分方程式为3方程的解即对应拟合函数的峰值位置。5多项式高斯公式拟合法多项式高斯公式拟合法原理为对数据采集卡采集的波形曲线数据进行高斯函数多项式变换,采用一般多项式拟合法的原理得到峰值位置。高斯函数为4对高斯函数进行对数变换,则有5通过一般多项式拟合法中的二次多项式拟合,可得到高斯函数变换后的系数6其中。把数据采集卡采集的波形数据分峰截幅后,取,用二次多项式拟合得到峰值位置。6高斯公式非线性曲线拟合法高斯公式拟合法原理为把数据采集卡采集的波形曲线数据直接作为高斯函数进行拟合处理,不经过多项式变换。运用莱文伯马克特LEVENBERGMARQUARDTLM算法和最小二乘法判定,拟合得出高斯函数4式的一组参数,满足输入数据点(X,Y)。LM算法提供了一个求解函数最小值的数值方法,是在高斯牛GAUSSNEWTONGN算法和非线性梯度下降算法之间插值,LM算法相比较GN算法,更能抵制噪声的影响,即使在初始值远离最终最小值的时候也可以精确得到解。在最小二乘法中,设给定实验数据X和FX,以及拟合曲线PX,要求拟合最佳,则要求满足最小二乘准则7利用LM算法可以得到关于函数的的最小值的解P,从而获得高斯公式非线性曲线拟合的各个系数,最终得到峰值位置。7对以上6种方法的总结输入信号的信噪比对于寻峰算法中算法误差的影响最大。信噪比大小的改变对算法精度的影响超过相同信噪比时各种算法间精度的差异,且信噪比对所有算法的影响基本一致。6种算法的误差随噪声幅度的增大而增大。仿真结果表明,在噪声幅度/信号幅度为00010080的范围内,算法误差与信噪比呈线性关系。由此可以得到,排除噪声对分析过程影响可以有效提高寻峰的精确度。本文采用小波分析方法,在MATLAB环境下,通过对小波系数处理,可以有效地抑制噪声对结果的影响,在谱峰识别精度上得到提高。12小波变换1从傅里叶变换到小波变换总所周知,自从1822年傅里叶发表“热传导解析理论”以来,傅里叶变换一直是信号处理领域中应用最广泛的一种分析手段。傅里叶变换的基本思想是将信号分解成一系列不同频率的连续正弦波的叠加,或者从另外一个角度来说是将信号从时间域转换到频率域。对于许多情况,傅里叶分析能很好地满足分析要求的。但是傅里叶变换有一个严重的不足,那就是在做变换时丢掉了时间信息,无法根据傅里叶变换的结果判断一个特定的信号是在什么时候发生的。也就是说,傅里叶变换只是一种纯频域的分析方法,它在频域里的定位是完全准确的,而在时域无任何定位性。如果要分析的信号是一种平稳信号,这一点也许不是很重要。然而实际中,大多数信号均含有大量的非稳态成分,例如偏移、突变等情况,而这些情况往往是相当重要的,反映了信号的重要特征。对这一类时变信号进行分析,通常需要提取某一时间段的频域信息或某一频率段所对应的时间信息。因此,需要寻求一种具有一定的时间和频率分辨率的基函数来分析时变信号。图11对平稳信号的频谱分析图12一个非频谱信号及其傅里叶谱为了研究信号在局部时间范围的频域特征,1946年GABOR提出了著名的GABOR变换,之后进一步发展成为短时傅里叶变换(STFT)。其基本思路是给信号加一个小窗,信号的傅里叶变换主要集中在对小窗里的信号进行变换,因此可以反映出信号局部特征。但由于STFT在许多领域获得了广泛的应用。但由于STFT的定义决定了其窗函数的大小和形状均与时间和频率无关而保持固定不变,这对分析时变信号来说是不利的。高频信号一般持续时间很短,而低频信号持续时间较长。因此,我们期望对于高频信号采用小时间窗,对于低频信号则采用大时间窗进行分析。在进行信号分析时,这种变时间窗的要求同STFT的固定时间窗的特性是相矛盾的。这表明STFT在处理这一类问题时已无能为力了。以上GABOR变换的不足之处,恰恰是小波变换的特长所在。小波变换不但继承和发展了STFT的局部化思想,而且克服了窗口大小不随频率变化、缺乏离散正交基的缺点,是一种比较理想的信号处理方法。2)小波变换小波(WAVELET),即小区域的波,是一种特殊的长度有限、平均值为0的波形。它有两个特点一是“小”;二是正负交替的“波动性”,也即直流分量为零。傅里叶分析是将信号分解成一系列不同频率的正弦波的叠加,同样小波分析是将信号分解成一系列小波函数的叠加,而这些小波函数都是由一个母小波函数经过平移和尺度伸缩得来的。由直觉,我们想到用不规则的小波函数来逼近尖锐变化的信号显然要比光滑的正弦曲线要好,同样,信号局部的特性用小波函数来逼近显然要比光滑的正弦函数来逼近更好。小波变换的定义是把某一被称为基本小波的函数做位移后,再在不同尺度下与待分析的信号FT做内积如上所述,小波分析的一个主要优点就是能够分析信号的局部特征,这为我们用小波变换的方法精确进行谱峰识别提供了思路。13MATLAB小波分析工具箱MATLAB,意即“矩阵实验室”。它是一种以矩阵作为基本数据单元的程序设计语言,提供了数据分析、算法实现与应用开发的交互式开发环境。一方面,MATLAB可以方便实现数值分析、优化分析、数据处理、自动控制、信号处理等领域的数学计算,另一方面,也可以快捷实现计算可视化、图形绘制、场景创建和渲染、图像处理、虚拟现实和地图制作等分析处理工作。MATLAB分为总包和若干个工具箱,小波工具箱使用小波分析技术进行信号分析、压缩和去噪处理。使用MATLAB语言进行小波分析,可以有两种基本方式命令方式和图形窗口方式。本文使用MATLAB辅助小波分析,对光谱信号进行小波变换后,通过对系数的选取处理,达到精确谱峰识别的目的。第二章小波分析基本原理21一维连续小波变换从数学上讲,小波是一个随时间振荡并很快衰减的函数,“小”是指它的衰减性,“波”是指它的波动性。211连续小波基函数如果函数,并满足以下“容许性”条件那么称是一个“母小波”(MOTHERWAVELET),其中是的FOURIER变换。我们很容易推导出,或者等价地在时域里有也就是说必须具有带通性质,且必是有正负交替的振荡波形,使得其平均值为零,这就是称为小波的原因。也就是说小波母函数满足两个条件1小它们在时域都具有紧支集或近似紧支集。2波动性由上面的容许性条件可知,小波的直流分量为0,因此可以断定小波具有正负交替的波动性。图21常见小波函数将小波母函数进行伸缩和平移,设伸缩因子(又称为尺度因子)为A,平移因子为,伸缩、平移后的函数为,则小波基函数定义为依赖于参数,。由于尺度因子T和平移因子是连续变化的,因此为连续小波基函数。它们为由同一个母函数经过拉伸和平移后获得的函数系。212连续小波变换连续小波变换的定义如下设FX是平方可积函数,是小波函数,则称为FX的小波变换(WT)。因为参数和是连续变化的,所以得到的就是连续小波变换(CWT)。与傅里叶变换类似,CWT也是一种积分变换,我们称为小波变换系数。由于小波基不同于傅里叶基,因此小波变换与傅里叶变换有许多不同之处,其中最重要的是,小波基具有尺度,平移两个参数,因此,将函数在小波基下展开,就意味着将一个时间函数投影到二维的时间尺度相平面上。213使用小波工具箱实现连续小波分析这里仅应用命令行方式。小波工具箱只需要一个函数来实现连续小波分析CWT。其格式为COEFSCWTS,SCALE,WNAME,PLOTS为待分析的信号。SCALES为连续小波变换的尺度向量,WNAME为小波函数名,PLOT用来画出小波变换后的系数的图形,系数的大小是以灰度的深浅来表示。22一维离散小波分析虽然从提取特征的角度看,常常还需要采用连续小波变换,但是在每个可能的尺度离散点都去计算小波系数,那将是个巨大工程,并且产生一大堆令人讨厌的数据。如果只取这些尺度的一小部分,以及部分时间点,将会大大减轻我们的工作量,并且不失准确性。在连续小波中,将尺度参数和平移参数分别离散化,取离散化的小波系数为逆变换为221单尺度一维离散小波分解1)分解X为被分析的离散信号,WNAME为分解所用到的小波函数,CA和CD分别为返回的低频系数和高频系数向量。CA,CDDWTX,WNAME2)从系数中重构低频和高频部分从第一步中产生的系数CA和CD构造第一层的低频和高频(A和D)系数YUPCOEFO,X,WNAME,N,L该函数用于计算向量X向上N步的重构小波系数。如果OA,则是对低频系数进行重构,如果OD,则是对高频系数进行重构。3显示低频和高频部分SUBPLOT1,2,1PLOTATITLEAPPROXIMATIONASUBPLOT1,2,2PLOTDTITLEDETAILD4由小波逆变换恢复信号XIDWTCA,CD,WNAME222多尺度一维分解1)WAVEDEC为多尺度一维小波分解函数,其格式为C,LWAVEDECX,N,WNAME它用小波完成对信号X的一维多尺度分解。N为尺度,且为严格的正函数。输出参数C由CAJ,CDJ,CDJ1,CD1组成,L由CAJ的长度,CDJ的长度,X的长度组成。2)提取系数的低频和高频部分APPCOEF用于从小波分解结构C,L中提取一维信号的低频系数,格式为AAPPCOEFC,L,WNAME,N其中,C,L为小波分解结构,WNAME为小波分解函数,N为尺度。DETCOEF用于从小波分解结构C,L中提取一维信号的高频系数DDETCOEFC,L,N3重构信号WRCOEF是对一维信号的分解结构C,L用指定的小波函数或重构滤波器进行重构的函数XWRCOEFTYPE,C,L,WNAME,N当TYPEA时,对信号的低频部分进行重构;当TYPED时,对信号的高频部分进行重构。4)显示多尺度一维分解结果其方法是SUBPLOT2,2,1PLOTA3TITLEAPPROXIMATIONA3SUBPLOT2,2,2PLOTD1TITLEDETAILD1SUBPLOT2,2,3PLOTD2TILTEDETAILD2SUBPLOT2,2,4PLOTD3TITLEDETAILD3图23多尺度一维分解结果5)重构原始信号XWAVERECC,L,WNAME这里WAVEREC为多尺度一维小波重构函数,用指定的小波函数对小波分解结构C,L进行重构。23信号的初步去噪一般来说,现实中的信号都是带噪信号,在对信号做进一步分析之前,需要将有效信号提取出来。目前,人们已根据噪声的统计特征和频谱分布规律,开发出了多种多样的信号去噪方法。其中最为直观的一种方法是,根据噪声能量一般集中于高频,而信号频谱分布于一个有限区间的特点,用傅里叶变换将含噪信号变换到频域,然后采用低通滤波器进行滤波。当信号和噪声的频带相互分离时这种方法比较有效,但当信号和噪声的频带相互重叠(比如含有白噪声)时,则效果较差,因为低通滤波器在抑制噪声的同时,也将信号的边缘部分变得模糊。利用小波对信号去噪,首先要识别出信号的哪一部分或哪些部分包含噪声,然后舍弃这些部分进行信号重构。对信号进行多尺度分解,在舍弃高频分量后,信号的低频部分变得越来越“纯洁”,但同时也丢失了初始信号的瞬变特征。因此,更优化的去噪是阈值去噪,即只去除那些超过某一设定值的细节部分。231小波分析用于消噪的基本原理一个含噪的一维信号模型可表示为如下形式SKFKPEK,K0,1,N1其中,SK为含噪信号,FK为有用信号,EK为噪声信号。这里我们认为EK是一个1级高斯白噪声,通常表现为高频信号,而工程实际中FK通常为低频信号,或者是一些比较平稳的信号。因此,我们可以按如下的方法进行消噪处理;首先对信号进行小波分解,一般地,噪声信号多包含在具有较高频率的细节中,从而可利用门限阈值等形式对所分解的小波系数进行处理,然后对信号进行小波重构即可达到对信号的消噪的目的。对信号消噪实质上是抑制信号中的无用部分,恢复信号中有用部分的过程。一般地,一维信号消噪的过程可分为如下三个步骤(1)一维信号的小波分解。选择一个小波并确定分解的层次,然后进行分解计算;(2)小波分解高频系数的阈值量化。对各个分解尺度下的高频系数选择一个阈值进行软阈值量化处理;(3)一维小波的重构。根据小波分解的最底层低频系数和各层高频系数进行一维小波重构。这三个步骤中,最关键的是如何选择阈值及如何进行阈值量化,在某种程度上,它关系到信号消噪的质量。232应用函数进行信号消噪处理消噪处理中两个非常有用的函数WDEN和WDENCMP。1WDEN的最简单用法如下SDWDENS,TPTR,SORH,SCAL,N,WAVENAME它所返回的是经过对原始信号S进行消噪处理后的信号SD。其中,TPTR指定阈值选取规则,SORH指定选取软阈值(SORHS)或硬阈值(SORHH),N为小波分解层数,WAVENAME指定分解时所用的小波。SCAL是阈值尺度改变的比例,它有如下三种选择(1)SCALONE,表示基本模式。(2)SCALSLN,表示未知尺度的基本模式,且仅根据第一层的小波分解系数来估计噪声的层次,并只进行一次估计,以此来变换阈值的尺度。(3)SCALMLN,表示非白噪声的基本模式,且在每个不同的小波分解的层次上都估计噪声的层次,以此来变换阈值的尺度。2WDENCMP是一种用的更为普遍的函数,它可以直接对一维或二维信号进行消噪或压缩,处理方法也是通过对小波分解系数进行阈值量化来实现。其使用方法如下XDWDENCMPOPT,X,WAVENAME,N,THR,SORH,KEEPAPP其中(1)OPTGBL,THR0,则阈值为全局阈值;OPTLVD,THR是向量,则阈值是在各层上大小不同的数值。(2)KEEPAPP1,不对小波分解后的系数作任何处理;KEEPAPP0,对小波分解后的低频系数也进行阈值量化处理。(3)X是待处理信号。(4)XD是处理后的信号。其余参数同函数WDEN中的参数。3)小波分析进行消噪处理选取阈值一般有下述三种方法(1)默认阈值消噪处理。该方法利用函数DDENCMP生成信号的默认阈值,然后利用函数WDENCMP进行消噪处理。(2)给定阈值消噪处理。在实际的消噪处理过程中,阈值往往可以通过经验公式获得,且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可用函数WTHRESH(3)强制消噪处理。该方法是将小波分解结构中的高频系数全部置0,即滤掉所有高频部分,然后对信号进行小波重构。这种方法比较简单,且消噪后的信号比较光滑,但是容易丢失信号中的有用成分。第三章基于MATLAB的程序设计本文选用光谱仪实际测得数据作为仿真对象,信号变量名为SIG图31。图3131设计流程对SIG运用软阈值法去噪多尺度变换;对高频系数作软阈值处理;重构信号SIG_DENOISE载入信号SIG建立结构数组CHAIN,存放处理链条(代表谱峰位置)所需的各种不同类型变量(链条总数COUNT_CHAIN,链条中点数目CHAINNCOUNT,极值点属性CHAINMTX等)。初始化CHAIN。在尺度2上,找出所有极大值系数所在点,设定每个点为一根链条的起点,建立链条的解析关系,将链条属性存入CHAIN中。遍历余下尺度,同上方法寻找各个尺度上极值点位置。比较极值点同每个链条的距离,若最小距离小于13,则该点归入相应链条;否则,建立新的链条,方法同上。在极值点画O,显示所有标志谱峰位置的链条。建立数组PEAKINFO记录谱峰信息,类似同尺度上寻极值方法,在每个链条上找到最佳变换尺度,将坐标、尺度值和对应系数存入PEAKINFO。对SIG_DENOISE一维连续小波变换,得到二维系数矩阵CW_SIG,每个元素代表对应点在对应尺度下变换系数,即小波函数与信号变化趋势的相关程度。32程序设计要点321载入并画出原始信号SIG。LOADSPECTRASIGMATFIGURESUBPLOT2,1,1PLOTSIGTITLE原始信号322使用软阈值法对信号去噪。具体思路是1)选用SYM8小波、尺度N10对SIG进行一维多尺度变换,得到参数矩阵C和L(见(1)。N10W_NAMESYM8C,LWAVEDECSIG,N,W_NAME(1)其中C由CA10,CD10,CD9,CD1组成,CA10代表低频系数,CD10,CD9,CD1代表各个尺度上的高频系数。L由CAJ的长度,CDJ的长度,X的长度组成。2)保持参数L不变,对系数矩阵C进行软阈值处理,得到CC矩阵。阈值处理过程如下对于C中前5个元素复制到CC中相应位置;确定阈值THRES见函数WAVE_THRES,遍历C中剩余元素当C(I)小于等于THRES时,CCI置零;当CI大于THRES时,CCICITHRES;当CILENGTHSIGWDLENGTHSIGX0ELSEIFX0WDTHRESCCICITHRESELSEIFCIN0WIDTHD_TMP1WIDTH1AKWIDTHKD_TMPWIDTH2N0KWIDTH1AK1N0D_TMPN0KWIDTH1N2AN0ELSED_TMPAKWIDTHKWIDTHENDENDIFABSD_TMPWIDTH1MAXD_TMP1E6LOCMK0IIDX0CONTINUEENDENDKFINDI0IDIK结论本文首先比较研究了一些常用谱峰识别的方法,进而从傅立叶分析引出更适合局部分析的小波方法。通过研究一维小波变换及其在去噪中的应用,基于MATLAB编写的程序实现了光谱图的小波降噪滤波,并在此基础上实现了谱峰的准确定位和谱峰宽度、强度等信息的提取。参考文献1徐金明,张孟喜,丁涛MATLAB实用教程清华大学出版社20052飞思科技产品研发中心MATLAB65辅助小波分析与应用电子工业出版社20033朱来东,廉小亲,江远志小波变换在信号降噪中的应用及MATLAB实现北京工商大学学报自然科学版20094茅耀斌小波信号处理及应用5王晓荣,程明霄谱峰识别的计算机设计与实现南京化工大学学报20016朱浩瀚,秦海琨,张敏,赖淑蓉,廖延彪光纤布拉格光栅传感解调中的寻峰算法中国激光20087TTONYCAI,DONGMAOZHANG,ANDDORBENAMOTZENHANCEDCHEMICALCLASSICATIONOFRAMANIMAGESUSINGMULTIRESOLUTIONWAVELETTRANSFORMATION致谢本文始终是在导师蔡志坚老师的悉心指导和热情帮助下完成的。在整个研究过程中,蔡老师都耐心为我排解疑问、提供建议和思路,在此论文完成之际,谨向蔡老师表达我最真挚的敬意和感谢在愉快的大学时光中,我认识了很多老师和同学,和他们一起学习、成长是我一生中永不褪色的回忆。谨向这些朝夕相处的老师和好友致意最后,我还要感谢我的家人,是他们的支持和关怀,使我一路成长,让我能够走到今天