毕业设计（论文）基于MATLAB小波变换的去噪应用

苏州大学本科生毕业设计（论文）目 录摘要 1一 绪论 21.1几种常用寻峰方法的简单说明 21.2小波变换 41.3 MATLAB小波分析工具箱 6二 小波分析基本原理 72.1一维连续小波分析 72.2一维离散小波分析 82.3信号的初步去噪12三 基于MATLAB的程序设计 14 3.1设计流程14 3.2程序设计要点16附：完整程序23结论 31参考文献32致谢33摘 要本文通过对染噪信号特征分析，设计了一种自动识别谱峰的程序。对比傅立叶分析的缺陷，小波方法在抑制噪音和局部分析中有着优异的性能。通过研究一维小波变换的基本原理，及其在信号去噪中的应用，基于MATLAB设计出的程序很好地完成了预期目标。关键词：小波变换 MATLAB 去噪 谱峰识别 AbstractBased on analysis of the characteristics of noised-signal,this paper is to design procedure to identify peaks.Contrast to Fourier analysis, Wavelets provide excellent compression and localization properties. By studying the the basic principles of one-dimensional wavelet transform and the application in signal denoising, the program designed based on MATLAB completes the targets well.Keywords wavelet transform; MATLAB; signal denoising; identifying peaks第一章 绪 论在我们的周围，每天都有大量的信号需要我们进行分析，例如我们说话的声音，机器的振动，金融变化数据，地震信号，音乐信号，医疗图像等。相当多的信号需要进行有效的编码，压缩，消噪，重建，建模和特征提取。对于实际应用中得到的光谱信号，我们有时需要对其进行去噪、识别、定位以及寻峰等操作。本章简要说明了几种常用谱峰识别方法，然后对本文用到的小波变换和MATLAB小波工具箱作基本的说明。1.1几种常用寻峰方法的简单说明1）蒙特卡罗算法蒙特卡罗(Monte Carlo) 算法 ,也称为统计实验方法,应用在寻峰算法中又称为质心探测法。原理为利用蒙特卡罗算法,把数据采集卡(DAQ) 采集的波形曲线数据进行分峰截幅后,作为质量非均匀的曲线段处理。波形数据中每点的横坐标值相当于质点系中各质点的位矢,纵坐标值相当于质点系中各质点的质量大小,质点系的质心横坐标可由质心定义式的蒙特卡罗算法求出。在波形轴对称或陡峭时,质心位置与波形峰值位置一致。2）直接比较法直接比较法即线形插值微分法。原理为对数据采集卡采集的波形曲线数据进行分峰截幅后,应用一阶数值微分,微分值0 点的位置即为原波形峰值位置。直接比较法应用前差公式或后差公式进行线形插值微分。两种公式的效果相同。3）二次插值数值微分法二次插值数值微分法为非线性插值数值微分法,其峰值位置判定原理和直接比较法一致。二次插值数值微分法应用中点公式,即三点公式进行非线性插值数值微分。4）一般多项式拟合法一般多项式拟合法原理为对数据采集卡采集的波形曲线数据进行分峰截幅后,采用一般多项式作拟合函数,最小二乘法作为判定,得到拟合式 (1)拟合多项式的一阶微分解析式为 (2)对应的一阶微分方程式为 (3)方程的解即对应拟合函数的峰值位置。5)多项式-高斯公式拟合法多项式-高斯公式拟合法原理为对数据采集卡采集的波形曲线数据进行高斯函数-多项式变换,采用一般多项式拟合法的原理得到峰值位置。高斯函数为 (4)对高斯函数进行对数变换,则有 (5)通过一般多项式拟合法中的二次多项式拟合,可得到高斯函数变换后的系数 (6)其中。把数据采集卡采集的波形数据分峰截幅后,取，用二次多项式拟合得到峰值位置。6)高斯公式非线性曲线拟合法高斯公式拟合法原理为把数据采集卡采集的波形曲线数据直接作为高斯函数进行拟合处理, 不经过多项式变换。运用莱文伯-马克特(Levenberg-Marquardt (L-M) ) 算法和最小二乘法判定,拟合得出高斯函数( (4) 式) 的一组参数,满足输入数据点（x,y）。L-M 算法提供了一个求解函数最小值的数值方法,是在高斯-牛( Gauss-Newton ( G-N) ) 算法和非线性梯度下降算法之间插值,L-M 算法相比较G-N 算法,更能抵制噪声的影响,即使在初始值远离最终最小值的时候也可以精确得到解。在最小二乘法中,设给定实验数据x 和f ( x) ,以及拟合曲线p ( x) , 要求拟合最佳, 则要求满足最小二乘准则 (7)利用L-M 算法可以得到关于函数的的最小值的解P , 从而获得高斯公式非线性曲线拟合的各个系数,最终得到峰值位置。7)对以上6种方法的总结输入信号的信噪比对于寻峰算法中算法误差的影响最大。信噪比大小的改变对算法精度的影响超过相同信噪比时各种算法间精度的差异,且信噪比对所有算法的影响基本一致。6 种算法的误差随噪声幅度的增大而增大。仿真结果表明,在噪声幅度/ 信号幅度为0. 0010. 080的范围内,算法误差与信噪比呈线性关系。由此可以得到，排除噪声对分析过程影响可以有效提高寻峰的精确度。本文采用小波分析方法，在MATLAB环境下，通过对小波系数处理，可以有效地抑制噪声对结果的影响，在谱峰识别精度上得到提高。1.2小波变换1).从傅里叶变换到小波变换总所周知，自从1822年傅里叶发表“热传导解析理论”以来，傅里叶变换一直是信号处理领域中应用最广泛的一种分析手段。傅里叶变换的基本思想是将信号分解成一系列不同频率的连续正弦波的叠加，或者从另外一个角度来说是将信号从时间域转换到频率域。对于许多情况，傅里叶分析能很好地满足分析要求的。但是傅里叶变换有一个严重的不足，那就是在做变换时丢掉了时间信息，无法根据傅里叶变换的结果判断一个特定的信号是在什么时候发生的。也就是说，傅里叶变换只是一种纯频域的分析方法，它在频域里的定位是完全准确的，而在时域无任何定位性。如果要分析的信号是一种平稳信号，这一点也许不是很重要。然而实际中，大多数信号均含有大量的非稳态成分，例如偏移、突变等情况，而这些情况往往是相当重要的，反映了信号的重要特征。对这一类时变信号进行分析，通常需要提取某一时间段的频域信息或某一频率段所对应的时间信息。因此，需要寻求一种具有一定的时间和频率分辨率的基函数来分析时变信号。 图1.1 对平稳信号的频谱分析 图1.2 一个非频谱信号及其傅里叶谱为了研究信号在局部时间范围的频域特征，1946年Gabor提出了著名的Gabor变换，之后进一步发展成为短时傅里叶变换（STFT）。其基本思路是给信号加一个小窗，信号的傅里叶变换主要集中在对小窗里的信号进行变换，因此可以反映出信号局部特征。但由于STFT在许多领域获得了广泛的应用。但由于STFT的定义决定了其窗函数的大小和形状均与时间和频率无关而保持固定不变，这对分析时变信号来说是不利的。高频信号一般持续时间很短，而低频信号持续时间较长。因此，我们期望对于高频信号采用小时间窗，对于低频信号则采用大时间窗进行分析。在进行信号分析时，这种变时间窗的要求同STFT的固定时间窗的特性是相矛盾的。这表明STFT在处理这一类问题时已无能为力了。以上Gabor变换的不足之处，恰恰是小波变换的特长所在。小波变换不但继承和发展了STFT的局部化思想，而且克服了窗口大小不随频率变化、缺乏离散正交基的缺点，是一种比较理想的信号处理方法。2）小波变换小波（wavelet）,即小区域的波，是一种特殊的长度有限、平均值为0的波形。它有两个特点：一是“小” ；二是正负交替的“波动性”，也即直流分量为零。傅里叶分析是将信号分解成一系列不同频率的正弦波的叠加，同样小波分析是将信号分解成一系列小波函数的叠加，而这些小波函数都是由一个母小波函数经过平移和尺度伸缩得来的。由直觉，我们想到用不规则的小波函数来逼近尖锐变化的信号显然要比光滑的正弦曲线要好，同样，信号局部的特性用小波函数来逼近显然要比光滑的正弦函数来逼近更好。小波变换的定义是把某一被称为基本小波的函数做位移后，再在不同尺度下与待分析的信号f(t)做内积：如上所述，小波分析的一个主要优点就是能够分析信号的局部特征，这为我们用小波变换的方法精确进行谱峰识别提供了思路。1.3 MATLAB小波分析工具箱MATLAB,意即“矩阵实验室” 。它是一种以矩阵作为基本数据单元的程序设计语言，提供了数据分析、算法实现与应用开发的交互式开发环境。一方面，MATLAB可以方便实现数值分析、优化分析、数据处理、自动控制、信号处理等领域的数学计算，另一方面，也可以快捷实现计算可视化、图形绘制、场景创建和渲染、图像处理、虚拟现实和地图制作等分析处理工作。MATLAB分为总包和若干个工具箱，小波工具箱使用小波分析技术进行信号分析、压缩和去噪处理。使用MATLAB语言进行小波分析，可以有两种基本方式：命令方式和图形窗口方式。本文使用MATLAB辅助小波分析，对光谱信号进行小波变换后，通过对系数的选取处理，达到精确谱峰识别的目的。第二章 小波分析基本原理2.1一维连续小波变换从数学上讲，小波是一个随时间振荡并很快衰减的函数，“小”是指它的衰减性，“波”是指它的波动性。2.1.1连续小波基函数如果函数，并满足以下“容许性”条件：那么称是一个“母小波”（mother wavelet） ，其中是的Fourier变换。我们很容易推导出，或者等价地在时域里有： 也就是说必须具有带通性质，且必是有正负交替的振荡波形，使得其平均值为零，这就是称为小波的原因。也就是说小波母函数满足两个条件：1. 小：它们在时域都具有紧支集或近似紧支集。2. 波动性：由上面的容许性条件可知，小波的直流分量为0，因此可以断定小波具有正负交替的波动性。图2.1 常见小波函数将小波母函数进行伸缩和平移，设伸缩因子（又称为尺度因子） 为a，平移因子为，伸缩、平移后的函数为，则小波基函数定义为：依赖于参数， 。由于尺度因子t和平移因子是连续变化的，因此为连续小波基函数。它们为由同一个母函数经过拉伸和平移后获得的函数系。2.1.2 连续小波变换连续小波变换的定义如下：设f(x)是平方可积函数，是小波函数，则称为f(x)的小波变换（WT）。因为参数和是连续变化的，所以得到的就是连续小波变换（CWT）。与傅里叶变换类似，CWT也是一种积分变换，我们称为小波变换系数。由于小波基不同于傅里叶基，因此小波变换与傅里叶变换有许多不同之处，其中最重要的是，小波基具有尺度，平移两个参数，因此，将函数在小波基下展开，就意味着将一个时间函数投影到二维的时间-尺度相平面上。2.1.3 使用小波工具箱实现连续小波分析这里仅应用命令行方式。小波工具箱只需要一个函数来实现连续小波分析： cwt。其格式为： coefs = cwt(s,scale,wname,plot)s为待分析的信号。Scales为连续小波变换的尺度向量，wname为小波函数名，plot用来画出小波变换后的系数的图形，系数的大小是以灰度的深浅来表示。2.2 一维离散小波分析虽然从提取特征的角度看，常常还需要采用连续小波变换，但是在每个可能的尺度离散点都去计算小波系数，那将是个巨大工程，并且产生一大堆令人讨厌的数据。如果只取这些尺度的一小部分，以及部分时间点，将会大大减轻我们的工作量，并且不失准确性。在连续小波中，将尺度参数和平移参数分别离散化，取离散化的小波系数为逆变换为2.2.1单尺度一维离散小波分解1）分解X为被分析的离散信号，wname为分解所用到的小波函数，cA和cD分别为返回的低频系数和高频系数向量。cA,cD=dwt(X,wname)2）从系数中重构低频和高频部分从第一步中产生的系数cA和cD构造第一层的低频和高频（A和D）系数：Y=upcoef(O,X,wname,N,L);该函数用于计算向量X向上N步的重构小波系数。如果O=a,则是对低频系数进行重构，如果O=d,则是对高频系数进行重构。3)显示低频和高频部分subplot(1,2,1);plot(A);title(Approximation A)subplot(1,2,2);plot(D);title(Detail D)4)由小波逆变换恢复信号X=idwt(cA,cD,wname);2.2.2多尺度一维分解1）wavedec为多尺度一维小波分解函数，其格式为：C,L=wavedec(X,N,wname)它用小波完成对信号X的一维多尺度分解。N为尺度，且为严格的正函数。输出参数C由cAj,cDj,cDj-1,cD1组成，L由cAj的长度，cDj的长度，X的长度组成。2）提取系数的低频和高频部分appcoef用于从小波分解结构C,L中提取一维信号的低频系数，格式为：A=appcoef(C,L,wname,N)其中，C,L为小波分解结构，wname为小波分解函数，N为尺度。detcoef用于从小波分解结构C,L中提取一维信号的高频系数：D=detcoef(C,L,N)3)重构信号wrcoef是对一维信号的分解结构C,L用指定的小波函数或重构滤波器进行重构的函数：X=wrcoef(type,C,L,wname,N)当type=a时，对信号的低频部分进行重构；当type=d时，对信号的高频部分进行重构。4）显示多尺度一维分解结果其方法是： subplot(2,2,1);plot(A3);title(Approximation A3) subplot(2,2,2);plot(D1);title(Detail D1) subplot(2,2,3);plot(D2);tilte(Detail D2) subplot(2,2,4);plot(D3);title(Detail D3) 图2.3 多尺度一维分解结果5）重构原始信号X=waverec(C,L,wname)这里waverec为多尺度一维小波重构函数，用指定的小波函数对小波分解结构C,L进行重构。2.3 信号的初步去噪一般来说，现实中的信号都是带噪信号，在对信号做进一步分析之前，需要将有效信号提取出来。目前，人们已根据噪声的统计特征和频谱分布规律，开发出了多种多样的信号去噪方法。其中最为直观的一种方法是，根据噪声能量一般集中于高频，而信号频谱分布于一个有限区间的特点，用傅里叶变换将含噪信号变换到频域，然后采用低通滤波器进行滤波。当信号和噪声的频带相互分离时这种方法比较有效，但当信号和噪声的频带相互重叠（比如含有白噪声）时，则效果较差，因为低通滤波器在抑制噪声的同时，也将信号的边缘部分变得模糊。利用小波对信号去噪，首先要识别出信号的哪一部分或哪些部分包含噪声，然后舍弃这些部分进行信号重构。对信号进行多尺度分解，在舍弃高频分量后，信号的低频部分变得越来越“纯洁”，但同时也丢失了初始信号的瞬变特征。因此，更优化的去噪是阈值去噪，即只去除那些超过某一设定值的细节部分。2.3.1 小波分析用于消噪的基本原理一个含噪的一维信号模型可表示为如下形式：S(k) = f(k) + p * e(k), k = 0,1,.,n-1其中，s(k)为含噪信号，f(k)为有用信号，e(k)为噪声信号。这里我们认为e(k)是一个1级高斯白噪声，通常表现为高频信号，而工程实际中f(k)通常为低频信号，或者是一些比较平稳的信号。因此，我们可以按如下的方法进行消噪处理；首先对信号进行小波分解，一般地，噪声信号多包含在具有较高频率的细节中，从而可利用门限阈值等形式对所分解的小波系数进行处理，然后对信号进行小波重构即可达到对信号的消噪的目的。对信号消噪实质上是抑制信号中的无用部分，恢复信号中有用部分的过程。一般地，一维信号消噪的过程可分为如下三个步骤：（1） 一维信号的小波分解。选择一个小波并确定分解的层次，然后进行分解计算；（2） 小波分解高频系数的阈值量化。对各个分解尺度下的高频系数选择一个阈值进行软阈值量化处理；（3） 一维小波的重构。根据小波分解的最底层低频系数和各层高频系数进行一维小波重构。这三个步骤中，最关键的是如何选择阈值及如何进行阈值量化，在某种程度上，它关系到信号消噪的质量。2.3.2 应用函数进行信号消噪处理消噪处理中两个非常有用的函数wden和wdencmp。1)wden的最简单用法如下：sd=wden(s,tptr,sorh,scal,n,wavename)它所返回的是经过对原始信号s进行消噪处理后的信号sd。其中，tptr指定阈值选取规则，sorh指定选取软阈值（sorh=s）或硬阈值（sorh=h）,n为小波分解层数，wavename指定分解时所用的小波。scal是阈值尺度改变的比例，它有如下三种选择：（1）scal=one,表示基本模式。（2）scal=sln,表示未知尺度的基本模式，且仅根据第一层的小波分解系数来估计噪声的层次，并只进行一次估计，以此来变换阈值的尺度。（3）scal=mln,表示非白噪声的基本模式，且在每个不同的小波分解的层次上都估计噪声的层次，以此来变换阈值的尺度。2)wdencmp是一种用的更为普遍的函数，它可以直接对一维或二维信号进行消噪或压缩，处理方法也是通过对小波分解系数进行阈值量化来实现。其使用方法如下：xd=wdencmp(opt,x,wavename,n,thr,sorh,keepapp)其中：（1） opt=gbl,thr>0,则阈值为全局阈值；opt=lvd,thr是向量，则阈值是在各层上大小不同的数值。（2） keepapp=1,不对小波分解后的系数作任何处理；keepapp=0,对小波分解后的低频系数也进行阈值量化处理。（3） x是待处理信号。（4） xd是处理后的信号。其余参数同函数wden中的参数。3）小波分析进行消噪处理选取阈值一般有下述三种方法：（1） 默认阈值消噪处理。该方法利用函数ddencmp生成信号的默认阈值，然后利用函数wdencmp进行消噪处理。（2） 给定阈值消噪处理。在实际的消噪处理过程中，阈值往往可以通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可用函数wthresh.（3） 强制消噪处理。该方法是将小波分解结构中的高频系数全部置0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且消噪后的信号比较光滑，但是容易丢失信号中的有用成分。第三章 基于MATLAB的程序设计本文选用光谱仪实际测得数据作为仿真对象，信号变量名为sig(图3.1)。 图3.13.1 设计流程对sig_denoise一维连续小波变换，得到二维系数矩阵cw_sig，每个元素代表对应点在对应尺度下变换系数，即小波函数与信号变化趋势的相关程度。对sig运用软阈值法去噪：多尺度变换；对高频系数作软阈值处理；重构信号sig_denoise载入信号sig 建立结构数组chain,存放处理链条（代表谱峰位置）所需的各种不同类型变量（链条总数count_chain,链条中点数目chain.ncount,极值点属性chain.Mtx等）。初始化chain。在尺度2上，找出所有极大值系数所在点，设定每个点为一根链条的起点，建立链条的解析关系，将链条属性存入chain中。遍历余下尺度，同上方法寻找各个尺度上极值点位置。比较极值点同每个链条的距离，若最小距离小于13，则该点归入相应链条；否则，建立新的链条，方法同上。在极值点画o，显示所有标志谱峰位置的链条。建立数组PeakInfo记录谱峰信息，类似同尺度上寻极值方法，在每个链条上找到最佳变换尺度，将坐标、尺度值和对应系数存入PeakInfo。验证峰附近信号信噪比，排除信噪比低于阈值的峰。用蓝线标志检测到的谱峰。程序结束3.2程序设计要点3.2.1载入并画出原始信号sig。load spectrasig.mat; figure; subplot(2,1,1); plot(sig); title('原始信号');3.2.2 使用软阈值法对信号去噪。具体思路是：1）选用sym8小波、尺度N=10对sig进行一维多尺度变换，得到参数矩阵C和L（见（1）。N = 10;w_name = 'sym8' C, L = wavedec(sig, 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时，cc(i)置零；当C(i)大于thres时，cc(i) = C(i) thres；当C(i)<-thres时，cc(i) = C(i) + thres。3）使用waverec对系数cc和L在sym8下重构信号，命名为sig_denoise（见（2）并显示去噪后的信号（见图3.2）。sig_denoise = waverec(cc, L, w_name);subplot(2,1,2);plot(sig_denoise);title('去噪声信号'); （2）3.2.3 使用小波mexh，在尺度为160下，对信号sig_denoise进行一维连续小波变换，得到二维的系数矩阵cw_sig（见（3）。cw_sig=cwt(sig_denoise,1:LV,'mexh','plot') （3）由小波变换的原理，一维连续小波变换是一种积分变换，小波基具有尺度,平移两个参数，因此，将函数在小波基下展开，就意味着将一个时间函数投影到二维的时间-尺度相平面上。即可得到二维的系数矩阵，其中横坐标上代表时域，纵坐标代表各个尺度（见图3.3）。矩阵中每个系数代表在相应尺度下相应时间点处变换系数，这个系数表示了小波函数同该点处附近信号趋势的相关性，系数值越大代表相关性越高。图 3.23.2.4 如3中所述，二维系数矩阵cw_sig中每个系数代表在相应尺度下相应时间点处的变换系数，这个系数表示了小波函数同该点处附近信号趋势的相关性，系数值越大代表相关性越高。整个空间信号有若干个谱峰，对应了单个尺度（纵坐标）上不同点处的若干个系数极大值；在单个点处，不同尺度下的变换均得到极大值系数。也就是说，在可能的谱峰位置，所有尺度上均可找到一个模极大值系数；不同尺度上系数有波动，其中的最大值处的尺度对应最佳变换尺度，即在此尺度上的变换可得到信号与小波函数有最高相关性。以下工作就是在各个尺度上遍历所有点找到并标记极大值系数，之后将各个极值系数对应的标记在各个尺度中顺次按一定规则连成链条，链条即对应可能的谱峰位置。3.2.5建立结构数组chain用来存放处理链条过程中不同类型属性的变量empty,Mtx,ncount,a,b,c等，初始化各种变量值（见（4）。for idx = 1:100 chain(idx).empty = 1;chain(idx).Mtx = zeros(100, 3); chain(idx).ncount = 0; chain(idx).a = 0;chain(idx).b = 0;chain(idx).c = 0;endcount_chain = 0; （4）图 3.33.2.6 寻找系数极大值点并标志出代表峰值的模极大值系数链条的位置。过程如下:1）在尺度2上，遍历所有点上系数值，可以得到若干模极大值及它们所在位置。实现途径是运行函数LM, IM = my_localmax(A, th, w)，其中A = cw_sig(2, :)；IM是一维数组，用来存放各个模极大值所在的空间坐标。建立数组XM和YM，将IM中元素复制到XM中，YM中存入当前处理的尺度即2，它们长度均等于IM长度。设定各个极值点是一个代表谱峰的链条的一部分，这样表示链条数量的count_chain=Nm=length(IM)。遍历所有极值点，将各个对应坐标，尺度和系数值存入数组chain_Mtx中（见（5）。for idx = 1:Nm chain(idx).empty = 0; K = IM(idx), 2, A(IM(idx); chain(idx).Mtx(1, 1:3) = K; chain(idx).ncount = 1; chain(idx).a = 1;chain(idx).b = 0;chain(idx).c = -IM(idx); AJ(idx) = 1;BJ(idx) = 0;CJ(idx) = -IM(idx);endcount_chain = Nm; （5）2）如4中所述，在可能的谱峰位置，所有尺度上均可找到一个模极大值系数；不同尺度上系数有波动，其中的最大值处的尺度对应最佳变换尺度。同上，利用LM, IM = my_localmax(A, th, w)分别找到其他所有尺度上的模极大点，将IM和尺度值分别存入XM和YM中。这些极大值点可以按规则连成若干链条，对应信号中各个谱峰。3）将各尺度上系数极值点形成链条思路：结构数组chain中保存了尺度2上的各个极值点坐标，即Nm条直线的表达式x * AJ + y * BJ + CJ=0，其中AJ=1,BJ=0,CJ=-IM，表示chain(idx),idx=1:Nm。在尺度3上，比较各个极值点与Nm条直线距离，将最小值和对应链条序号存入mD,iD中。设定阈值距离th_distance=13,若mD<th_distance，将该点归入对应链条；否则，链条总数count_chain加1，并将该点坐标，尺度和系数值存入chain.Mtx中（见（6）。对于mD<th_distance情况下，链条上每增加一个点，chain.ncount加1；建立数组XLine和YLine，XLine存放链条上各个点横坐标，YLine存放尺度值；若var(XLine)<2，chain中CJ=-mean(XLine)，否则对XLine和YLine运用最小二乘法进行曲线拟合得到的系数p(1)和p(2)赋给AJ和CJ，BJ=-1。for idx = 1:Nm x = IM(idx); y = lev; D = abs(x * AJ + y * BJ + CJ)./(sqrt(AJ.2+BJ.2); mD, iD = min(D); K = x, y, A(x); if mD < th_distance chain(iD).ncount = chain(iD).ncount + 1; chain(iD).Mtx(chain(iD).ncount, 1:3) = K; XLine = chain(iD).Mtx(1:chain(iD).ncount, 1); YLine = chain(iD).Mtx(1:chain(iD).ncount, 2);else count_chain = count_chain + 1; chain(count_chain).empty = 0; chain(count_chain).Mtx(1, 1:3) = K; chain(count_chain).ncount = 1; （6）4）在（XM，YM）表示的所有点处画出o(见（7）)，可以清楚地看到代表谱峰的各个链条（见图3.4）。figure; plot(XM, YM, 'o'); （7）图 3.43.2.7建立数组PeakInfo记录谱峰信息：类似同尺度上寻极值方法，在每个链条上找到最佳变换尺度，将坐标、尺度值和对应系数存入PeakInfo。PeakInfo = zeros(count_chain, 3);for idx = 1:count_chain NJ = chain(idx).ncount; I_chain = chain(idx).Mtx(1:NJ, 3); LM_c, IM_c = my_localmax(I_chain, th, w); if isempty(IM_c) = 1 LM_c, k = max(I_chain); else k = IM_c(1); end K = chain(idx).Mtx(k, 1:3); PeakInfo(idx, 1:3) = K(1:3);end3.2.8 验证峰附近信号信噪比不规则噪音叠加到信号上，如果噪音瞬时扰动过大则产生新的谱峰，显然这些谱峰是必须要排除的，思路是：对于已经检测到的谱峰，考察其附近信号的信噪比，若低于阈值，则认为这个峰是由噪音扰动形成的。由上述PeakInfo记录了所有检测到的谱峰信息（坐标，尺度，系数值），对于每个峰，生成附近信号信噪比SNR_d，设定阈值3,若SNR_d<3，则放弃这个位置的谱峰。SNR_d = (p.a)/(r.rmse); if SNR_d < 3 PeakFlag(idx) = 0; else PeakFlag(idx) = 1;3.2.9 画出检测到的谱峰（见图3.5）。for idx = 1:count_chain if PeakFlag(idx) = 1 X = PeakInfo(idx, 1), PeakInfo(idx, 1); Y = 0, Max_sig; plot(X, Y, 'blue'); endend图 3.5附：完整程序1）main.mclear all; clc; close all;% load original signal and plot;load spectrasig.mat;%load simuB.mat;figure; subplot(2,1,1); plot(sig); title('原始信号');% denoiseN = 10;w_name = 'sym8'C, L = wavedec(sig, N, w_name);K = floor(0.5 * N);cc = wave_thres(C, L, K);sig_denoise = waverec(cc, L, w_name);subplot(2,1,2); plot(sig_denoise); title('去噪声信号');% cwt by mexican hat waveletLV = 60;figure; cw_sig = cwt(sig_denoise, 1:LV, 'mexh', 'plot');% initializee chainfor idx = 1:100 chain(idx).empty = 1; chain(idx).Mtx = zeros(100, 3); chain(idx).ncount = 0; chain(idx).a = 0;chain(idx).b = 0; chain(idx).c = 0;endcount_chain = 0;% get local max,draw max linesL = length(sig_denoise);A = cw_sig(2, :);th = 0.01;w = 1;LM, IM = my_localmax(A, th, w);Nm = length(IM);XM(1:Nm) = IM(1:Nm);YM(1:Nm) = 2 * ones(1, Nm);tm = Nm;for idx = 1:Nm chain(idx).empty = 0; K = IM(idx), 2, A(IM(idx); chain(idx).Mtx(1, 1:3) = K;chain(idx).ncount = 1; chain(idx).a = 1;chain(idx).b = 0;chain(idx).c = -IM(idx); AJ(idx) = 1;BJ(idx) = 0;CJ(idx) = -IM(idx);endcount_chain = Nm;th_distance = 13;for lev = 3:LV A = cw_sig(lev, :); w = floor(lev * 0.7); LM, IM = my_localmax(A, th, w); Nm = length(IM); XM(tm+1):(tm+Nm) = IM(1:Nm); YM(tm+1):(tm+Nm) = lev; tm = tm + Nm; for idx = 1:Nm x = IM(idx);y = lev; D = abs(x * AJ + y * BJ + CJ)./(sqrt(AJ.2+BJ.2); mD, iD = min(D); K = x, y, A(x); if mD < th_distance chain(iD).ncount = chain(iD).ncount + 1; chain(iD).Mtx(chain(iD).ncount, 1:3) = K; XLine = chain(iD).Mtx(1:chain(iD).ncount, 1); YLine = chain(iD).Mtx(1:chain(iD).ncount, 2); if var(XLine) < 2 chain(iD).a = 1;AJ(iD) = 1; chain(iD).b = 0; BJ(iD) = 0; chain(iD).c = - mean(XLine); CJ(iD) = -mean(XLine); else p = polyfit(XLine, YLine, 1); chain(iD).a = p(1);AJ(iD) = p(1); chain(iD).b = -1；BJ(iD) = -1; chain(iD).c = p(2)；CJ(iD) = p(2); end else count_chain = count_chain + 1; chain(count_chain).empty = 0; chain(count_chain).Mtx(1, 1:3) = K; chain(count_chain).ncount = 1; chain(count_chain).a = 1; chain(count_chain).b = 0; chain(count_chain).c = -IM(idx); AJ(count_chain) = 1; BJ(count_chain) = 0; CJ(count_chain) = -IM(idx); end endendfigure; plot(XM, YM, 'o');w = 2;th = 0.001;PeakInfo = zeros(count_chain, 3);for idx = 1:count_chain NJ = chain(idx).ncount; I_chain = chain(idx).Mtx(1:NJ, 3); LM_c, IM_c = my_localmax(I_chain, th, w); if isempty(IM_c) = 1 LM_c, k = max(I_chain); else k = IM_c(1); end K = chain(idx).Mtx(k, 1:3); PeakInfo(idx, 1:3) = K(1:3);end% verify the SNR of the peakfor idx = 1:count_chain X0 = PeakInfo(idx, 1);scale = PeakInfo(idx, 2);I0 = PeakInfo(idx, 3); sigma = floor(scale/1.66) + 1;wd = 3 * sigma; if X0 + wd > length(sig) wd = length(sig) - X0; elseif X0 - wd < 1 wd = X0 - 1; end OrigData = sig(X0-wd):(X0+wd); OrigData = OrigData - min(OrigData); OrigData = reshape(OrigData, length(OrigData), 1); XLL = 1:(2*wd + 1);XLL = XLL' Xp = wd + 1; C1 = scale/1.66; f = fittype('k * x + b + a * exp(-(x-Xp)/C1)2)', 'problem', 'Xp', 'C1'); p, r = fit(XLL, OrigData, f, 'problem', Xp, C1); SNR_d = (p.a)/(r.rmse); if SNR_d < 3 PeakFlag(idx) = 0; else PeakFlag(idx) = 1; endend% plot the detected peaksfigure; plot(sig, 'red'); hold on; Max_sig = max(sig);for idx = 1:count_chain if PeakFlag(idx) = 1 X = PeakInfo(idx, 1), PeakInfo(idx, 1); Y = 0, Max_sig; hold on; plot(X, Y, 'blue'); endend2) function cc = wave_thres(C, L, K)cc = zeros(size(C);m = length(C);n = sum(L(1:K);cc(1:n) = C(1:n);f11 = abs(C(n+1:m);omiga = median(f11)/0.6745;thres = 0.5 * omiga * sqrt( 2*log(m) );for i = (n+1) : m if abs(C(i) < thres cc(i)=0; elseif C(i) = thres cc(i)=0; elseif C(i) > thres cc(i) = C(i) - thres; elseif C(i) < -thres cc(i) = C(i) + thres; endend3) function locM, Id = my_localmax(A, th, width)lm, Im = Lmax(A);locM = zeros(size(lm); I = zeros(size(Im);locM = lm; I = Im; N0 = length(A); N1 = length(I);N2 = 2 * width + 1; D_Tmp = zeros(1, N2);for idx = 1:N1 K = I(idx); % remove localmax < th if A(K) < th locM(K) = 0; I(idx) = 0; continue; end % get the neighboring data if N0 < 2 * width + 1 X1 = width - K + 2; X2 = N0 - K + 1 + width; D_Tmp(X1:X2) = A(1:N0); D_Tmp(1:X1) = A(1); D_Tmp(X2:N2) = A(N0); else if K < width+1 D_Tmp(width+1):N2) = A(K:(K+width); D_Tmp(width+2-K):width) = A(1:(K-1); D_Tmp(1:(width+2-K) = A(1); elseif K > N0 - width D_Tmp(1:(width+1) = A(K-width):K);