畢業(yè)論文---自動(dòng)識(shí)別譜峰的程序設(shè)計(jì)

1、<p><b>　　目 錄</b></p><p>　　摘要…………………………………………………………………………… 1</p><p>　　一 緒論……………………………………………………………………… 2</p><p>　　1.1幾種常用尋峰方法的簡單說明……………………………………… 2</p><p&

2、gt;　　1.2小波變換……………………………………………………………… 4</p><p>　　1.3 MATLAB小波分析工具箱…………………………………………… 6</p><p>　　二 小波分析基本原理……………………………………………………… 7</p><p>　　2.1一維連續(xù)小波分析…………………………………………………… 7</p>

3、<p>　　2.2一維離散小波分析…………………………………………………… 8</p><p>　　2.3信號(hào)的初步去噪………………………………………………………12</p><p>　　三 基于MATLAB的程序設(shè)計(jì)…………………………………………… 14</p><p>　　3.1設(shè)計(jì)流程………………………………………………………………14</

4、p><p>　　3.2程序設(shè)計(jì)要點(diǎn)…………………………………………………………16</p><p>　　附：完整程序…………………………………………………………………23</p><p>　　結(jié)論 ………………………………………………………………………… 31</p><p>　　參考文獻(xiàn)………………………………………………………………………32

5、</p><p>　　致謝……………………………………………………………………………33</p><p><b>　　摘 要</b></p><p>　　本文通過對(duì)染噪信號(hào)特征分析，設(shè)計(jì)了一種自動(dòng)識(shí)別譜峰的程序。對(duì)比傅立葉分析的缺陷，小波方法在抑制噪音和局部分析中有著優(yōu)異的性能。通過研究一維小波變換的基本原理，及其在信號(hào)去噪中的應(yīng)用，基于MAT

6、LAB設(shè)計(jì)出的程序很好地完成了預(yù)期目標(biāo)。</p><p>　　關(guān)鍵詞：小波變換 MATLAB 去噪 譜峰識(shí)別</p><p><b>　　Abstract</b></p><p>　　Based on analysis of the characteristics of noised-signal,this paper is to des

7、ign 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 applicat

8、ion in signal denoising, the program designed based on MATLAB completes the targets well.</p><p>　　Keywords wavelet transform; MATLAB; signal denoising; </p><p>　　identifying peaks</p>

9、<p>　　第一章 緒 論</p><p>　　在我們的周圍，每天都有大量的信號(hào)需要我們進(jìn)行分析，例如我們說話的聲音，機(jī)器的振動(dòng)，金融變化數(shù)據(jù)，地震信號(hào)，音樂信號(hào)，醫(yī)療圖像等。相當(dāng)多的信號(hào)需要進(jìn)行有效的編碼，壓縮，消噪，重建，建模和特征提取。對(duì)于實(shí)際應(yīng)用中得到的光譜信號(hào)，我們有時(shí)需要對(duì)其進(jìn)行去噪、識(shí)別、定位以及尋峰等操作。</p><p>　　本章簡要說明了幾種常用譜峰

10、識(shí)別方法，然后對(duì)本文用到的小波變換和MATLAB小波工具箱作基本的說明。</p><p>　　1.1幾種常用尋峰方法的簡單說明</p><p><b>　　1）蒙特卡羅算法</b></p><p>　　蒙特卡羅(Monte Carlo) 算法 ,也稱為統(tǒng)計(jì)實(shí)驗(yàn)方法,應(yīng)用在尋峰算法中又稱為質(zhì)心探測(cè)法。原理為利用蒙特卡羅算法,把數(shù)據(jù)采集卡(DAQ)

11、 采集的波形曲線數(shù)據(jù)進(jìn)行分峰截幅后,作為質(zhì)量非均勻的曲線段處理。波形數(shù)據(jù)中每點(diǎn)的橫坐標(biāo)值相當(dāng)于質(zhì)點(diǎn)系中各質(zhì)點(diǎn)的位矢,縱坐標(biāo)值相當(dāng)于質(zhì)點(diǎn)系中各質(zhì)點(diǎn)的質(zhì)量大小,質(zhì)點(diǎn)系的質(zhì)心橫坐標(biāo)可由質(zhì)心定義式的蒙特卡羅算法求出。在波形軸對(duì)稱或陡峭時(shí),質(zhì)心位置與波形峰值位置一致。</p><p><b>　　2）直接比較法</b></p><p>　　直接比較法即線形插值微分法。原理為對(duì)數(shù)

12、據(jù)采集卡采集的波形曲線數(shù)據(jù)進(jìn)行分峰截幅后,應(yīng)用一階數(shù)值微分,微分值0 點(diǎn)的位置即為原波形峰值位置。直接比較法應(yīng)用前差公式或后差公式進(jìn)行線形插值微分。兩種公式的效果相同。</p><p>　　3）二次插值數(shù)值微分法</p><p>　　二次插值數(shù)值微分法為非線性插值數(shù)值微分法,其峰值位置判定原理和直接比較法一致。二次插值數(shù)值微分法應(yīng)用中點(diǎn)公式,即三點(diǎn)公式進(jìn)行非線性插值數(shù)值微分。</p&

13、gt;<p>　　4）一般多項(xiàng)式擬合法</p><p>　　一般多項(xiàng)式擬合法原理為對(duì)數(shù)據(jù)采集卡采集的波形曲線數(shù)據(jù)進(jìn)行分峰截幅后,采用一般多項(xiàng)式作擬合函數(shù),最小二乘法作為判定,得到擬合式</p><p><b>　　(1)</b></p><p>　　擬合多項(xiàng)式的一階微分解析式為</p><p><b&

14、gt;　　(2)</b></p><p>　　對(duì)應(yīng)的一階微分方程式為</p><p><b>　　(3)</b></p><p>　　方程的解即對(duì)應(yīng)擬合函數(shù)的峰值位置。</p><p>　　5)多項(xiàng)式-高斯公式擬合法</p><p>　　多項(xiàng)式-高斯公式擬合法原理為對(duì)數(shù)據(jù)采集卡采集的波

15、形曲線數(shù)據(jù)進(jìn)行高斯函數(shù)-多項(xiàng)式變換,采用一般多項(xiàng)式擬合法的原理得到峰值位置。</p><p><b>　　高斯函數(shù)為</b></p><p><b>　　(4)</b></p><p>　　對(duì)高斯函數(shù)進(jìn)行對(duì)數(shù)變換,則有</p><p><b>　　(5)</b></p&g

16、t;<p>　　通過一般多項(xiàng)式擬合法中的二次多項(xiàng)式擬合,可得到高斯函數(shù)變換后的系數(shù)</p><p><b>　　(6)</b></p><p>　　其中。把數(shù)據(jù)采集卡采集的波形數(shù)據(jù)分峰截幅后,取，用二次多項(xiàng)式擬合得到峰值位置。</p><p>　　6)高斯公式非線性曲線擬合法</p><p>　　高斯公式擬

17、合法原理為把數(shù)據(jù)采集卡采集的波形曲線數(shù)據(jù)直接作為高斯函數(shù)進(jìn)行擬合處理, 不經(jīng)過多項(xiàng)式變換。運(yùn)用萊文伯-馬克特(Levenberg-Marquardt (L-M) ) 算法和最小二乘法判定,擬合得出高斯函數(shù)( (4) 式) 的一組參數(shù),滿足輸入數(shù)據(jù)點(diǎn)（x,y）。</p><p>　　L-M 算法提供了一個(gè)求解函數(shù)最小值的數(shù)值方法,是在高斯-牛( Gauss-Newton ( G-N) ) 算法和非線性梯度下降算法之

18、間插值,L-M 算法相比較G-N 算法,更能抵制噪聲的影響,即使在初始值遠(yuǎn)離最終最小值的時(shí)候也可以精確得到解。</p><p>　　在最小二乘法中,設(shè)給定實(shí)驗(yàn)數(shù)據(jù)x 和f ( x) ,以及擬合曲線p ( x) , 要求擬合最佳, 則要求滿足最小二乘準(zhǔn)則</p><p><b>　　(7)</b></p><p>　　利用L-M 算法可以得到關(guān)于

19、函數(shù)的</p><p>　　的最小值的解P , 從而獲得高斯公式非線性曲線擬合的各個(gè)系數(shù),最終得到峰值位置。</p><p>　　7)對(duì)以上6種方法的總結(jié)</p><p>　　輸入信號(hào)的信噪比對(duì)于尋峰算法中算法誤差的影響最大。信噪比大小的改變對(duì)算法精度的影響超過相同信噪比時(shí)各種算法間精度的差異,且信噪比對(duì)所有算法的影響基本一致。6 種算法的誤差隨噪聲幅度的增大而增大

20、。仿真結(jié)果表明,在噪聲幅度/ 信號(hào)幅度為0. 001～0. 080的范圍內(nèi),算法誤差與信噪比呈線性關(guān)系。由此可以得到，排除噪聲對(duì)分析過程影響可以有效提高尋峰的精確度。本文采用小波分析方法，在MATLAB環(huán)境下，通過對(duì)小波系數(shù)處理，可以有效地抑制噪聲對(duì)結(jié)果的影響，在譜峰識(shí)別精度上得到提高。</p><p><b>　　1.2小波變換</b></p><p>　　1).從

21、傅里葉變換到小波變換</p><p>　　總所周知，自從1822年傅里葉發(fā)表“熱傳導(dǎo)解析理論”以來，傅里葉變換一直是信號(hào)處理領(lǐng)域中應(yīng)用最廣泛的一種分析手段。傅里葉變換的基本思想是將信號(hào)分解成一系列不同頻率的連續(xù)正弦波的疊加，或者從另外一個(gè)角度來說是將信號(hào)從時(shí)間域轉(zhuǎn)換到頻率域。對(duì)于許多情況，傅里葉分析能很好地滿足分析要求的。</p><p>　　但是傅里葉變換有一個(gè)嚴(yán)重的不足，那就是在做變換

22、時(shí)丟掉了時(shí)間信息，無法根據(jù)傅里葉變換的結(jié)果判斷一個(gè)特定的信號(hào)是在什么時(shí)候發(fā)生的。也就是說，傅里葉變換只是一種純頻域的分析方法，它在頻域里的定位是完全準(zhǔn)確的，而在時(shí)域無任何定位性。</p><p>　　如果要分析的信號(hào)是一種平穩(wěn)信號(hào)，這一點(diǎn)也許不是很重要。然而實(shí)際中，大多數(shù)信號(hào)均含有大量的非穩(wěn)態(tài)成分，例如偏移、突變等情況，而這些情況往往是相當(dāng)重要的，反映了信號(hào)的重要特征。對(duì)這一類時(shí)變信號(hào)進(jìn)行分析，通常需要提取某一時(shí)

23、間段的頻域信息或某一頻率段所對(duì)應(yīng)的時(shí)間信息。因此，需要尋求一種具有一定的時(shí)間和頻率分辨率的基函數(shù)來分析時(shí)變信號(hào)。</p><p>　　圖1.1 對(duì)平穩(wěn)信號(hào)的頻譜分析</p><p>　　圖1.2 一個(gè)非頻譜信號(hào)及其傅里葉譜</p><p>　　為了研究信號(hào)在局部時(shí)間范圍的頻域特征，1946年Gabor提出了著名的Gabor變換，之后進(jìn)一步發(fā)展成為短時(shí)傅里葉變換（

24、STFT）。其基本思路是給信號(hào)加一個(gè)小窗，信號(hào)的傅里葉變換主要集中在對(duì)小窗里的信號(hào)進(jìn)行變換，因此可以反映出信號(hào)局部特征。但由于STFT在許多領(lǐng)域獲得了廣泛的應(yīng)用。但由于STFT的定義決定了其窗函數(shù)的大小和形狀均與時(shí)間和頻率無關(guān)而保持固定不變，這對(duì)分析時(shí)變信號(hào)來說是不利的。高頻信號(hào)一般持續(xù)時(shí)間很短，而低頻信號(hào)持續(xù)時(shí)間較長。因此，我們期望對(duì)于高頻信號(hào)采用小時(shí)間窗，對(duì)于低頻信號(hào)則采用大時(shí)間窗進(jìn)行分析。在進(jìn)行信號(hào)分析時(shí)，這種變時(shí)間窗的要求同ST

25、FT的固定時(shí)間窗的特性是相矛盾的。這表明STFT在處理這一類問題時(shí)已無能為力了。</p><p>　　以上Gabor變換的不足之處，恰恰是小波變換的特長所在。小波變換不但繼承和發(fā)展了STFT的局部化思想，而且克服了窗口大小不隨頻率變化、缺乏離散正交基的缺點(diǎn)，是一種比較理想的信號(hào)處理方法。</p><p><b>　　2）小波變換</b></p><

26、p>　　小波（wavelet）,即小區(qū)域的波，是一種特殊的長度有限、平均值為0的波形。它有兩個(gè)特點(diǎn)：一是“小” ；二是正負(fù)交替的“波動(dòng)性”，也即直流分量為零。傅里葉分析是將信號(hào)分解成一系列不同頻率的正弦波的疊加，同樣小波分析是將信號(hào)分解成一系列小波函數(shù)的疊加，而這些小波函數(shù)都是由一個(gè)母小波函數(shù)經(jīng)過平移和尺度伸縮得來的。由直覺，我們想到用不規(guī)則的小波函數(shù)來逼近尖銳變化的信號(hào)顯然要比光滑的正弦曲線要好，同樣，信號(hào)局部的特性用小波函數(shù)來

27、逼近顯然要比光滑的正弦函數(shù)來逼近更好。</p><p>　　小波變換的定義是把某一被稱為基本小波的函數(shù)做位移τ后，再在不同尺度α下與待分析的信號(hào)f(t)做內(nèi)積：</p><p>　　如上所述，小波分析的一個(gè)主要優(yōu)點(diǎn)就是能夠分析信號(hào)的局部特征，這為我們用小波變換的方法精確進(jìn)行譜峰識(shí)別提供了思路。</p><p>　　1.3 MATLAB小波分析工具箱</p>

28、;<p>　　MATLAB,意即“矩陣實(shí)驗(yàn)室” 。它是一種以矩陣作為基本數(shù)據(jù)單元的程序設(shè)計(jì)語言，提供了數(shù)據(jù)分析、算法實(shí)現(xiàn)與應(yīng)用開發(fā)的交互式開發(fā)環(huán)境。一方面，MATLAB可以方便實(shí)現(xiàn)數(shù)值分析、優(yōu)化分析、數(shù)據(jù)處理、自動(dòng)控制、信號(hào)處理等領(lǐng)域的數(shù)學(xué)計(jì)算，另一方面，也可以快捷實(shí)現(xiàn)計(jì)算可視化、圖形繪制、場(chǎng)景創(chuàng)建和渲染、圖像處理、虛擬現(xiàn)實(shí)和地圖制作等分析處理工作。</p><p>　　MATLAB分為總包和若干個(gè)

29、工具箱，小波工具箱使用小波分析技術(shù)進(jìn)行信號(hào)分析、壓縮和去噪處理。使用MATLAB語言進(jìn)行小波分析，可以有兩種基本方式：命令方式和圖形窗口方式。</p><p>　　本文使用MATLAB輔助小波分析，對(duì)光譜信號(hào)進(jìn)行小波變換后，通過對(duì)系數(shù)的選取處理，達(dá)到精確譜峰識(shí)別的目的。</p><p>　　第二章 小波分析基本原理</p><p>　　2.1一維連續(xù)小波變換&l

30、t;/p><p>　　從數(shù)學(xué)上講，小波是一個(gè)隨時(shí)間振蕩并很快衰減的函數(shù)，“小”是指它的衰減性，“波”是指它的波動(dòng)性。</p><p>　　2.1.1連續(xù)小波基函數(shù)</p><p>　　如果函數(shù)，并滿足以下“容許性”條件：</p><p>　　那么稱是一個(gè)“母小波”（mother wavelet） ，其中是的Fourier變換。</p>

31、<p>　　我們很容易推導(dǎo)出，或者等價(jià)地在時(shí)域里有：</p><p>　　也就是說必須具有帶通性質(zhì)，且必是有正負(fù)交替的振蕩波形，使得其平均值為零，這就是稱為小波的原因。</p><p>　　也就是說小波母函數(shù)滿足兩個(gè)條件：</p><p>　　1. ?。核鼈?cè)跁r(shí)域都具有緊支集或近似緊支集。</p><p>　　2. 波動(dòng)性：由上

32、面的容許性條件可知，小波的直流分量為0，因此可以斷定小波具有正負(fù)交替的波動(dòng)性。</p><p>　　圖2.1 常見小波函數(shù)</p><p>　　將小波母函數(shù)進(jìn)行伸縮和平移，設(shè)伸縮因子（又稱為尺度因子） 為a，平移因子為τ，伸縮、平移后的函數(shù)為，則小波基函數(shù)定義為：</p><p>　　依賴于參數(shù)α，τ 。由于尺度因子t和平移因子τ是連續(xù)變化的，因此為連續(xù)小波基函

33、數(shù)。它們?yōu)橛赏粋€(gè)母函數(shù)經(jīng)過拉伸和平移后獲得的函數(shù)系。</p><p>　　2.1.2 連續(xù)小波變換</p><p>　　連續(xù)小波變換的定義如下：設(shè)f(x)是平方可積函數(shù)，是小波函數(shù)，則</p><p>　　稱為f(x)的小波變換（WT）。因?yàn)閰?shù)α和τ是連續(xù)變化的，所以得到的就是連續(xù)小波變換（CWT）。與傅里葉變換類似，CWT也是一種積分變換，我們稱為小波變換系數(shù)

34、。由于小波基不同于傅里葉基，因此小波變換與傅里葉變換有許多不同之處，其中最重要的是，小波基具有尺度α，平移τ兩個(gè)參數(shù)，因此，將函數(shù)在小波基下展開，就意味著將一個(gè)時(shí)間函數(shù)投影到二維的時(shí)間-尺度相平面上。</p><p>　　2.1.3 使用小波工具箱實(shí)現(xiàn)連續(xù)小波分析</p><p>　　這里僅應(yīng)用命令行方式。小波工具箱只需要一個(gè)函數(shù)來實(shí)現(xiàn)連續(xù)小波分析： cwt。其格式為：</p>

35、<p>　　coefs = cwt(s,scale,’wname’,’plot’)</p><p>　　s為待分析的信號(hào)。Scales為連續(xù)小波變換的尺度向量，wname為小波函數(shù)名，plot用來畫出小波變換后的系數(shù)的圖形，系數(shù)的大小是以灰度的深淺來表示。</p><p>　　2.2 一維離散小波分析</p><p>　　雖然從提取特征的角度看，常常還

36、需要采用連續(xù)小波變換，但是在每個(gè)可能的尺度離散點(diǎn)都去計(jì)算小波系數(shù)，那將是個(gè)巨大工程，并且產(chǎn)生一大堆令人討厭的數(shù)據(jù)。如果只取這些尺度的一小部分，以及部分時(shí)間點(diǎn)，將會(huì)大大減輕我們的工作量，并且不失準(zhǔn)確性。</p><p>　　在連續(xù)小波中，將尺度參數(shù)α和平移參數(shù)τ分別離散化，取</p><p><b>　　離散化的小波系數(shù)為</b></p><p>

37、;<b>　　逆變換為</b></p><p>　　2.2.1單尺度一維離散小波分解</p><p><b>　　1）分解</b></p><p>　　X為被分析的離散信號(hào)，wname為分解所用到的小波函數(shù)，cA和cD分別為返回的低頻系數(shù)和高頻系數(shù)向量。</p><p>　　[cA,cD]=dwt(

38、X,’wname’)</p><p>　　2）從系數(shù)中重構(gòu)低頻和高頻部分</p><p>　　從第一步中產(chǎn)生的系數(shù)cA和cD構(gòu)造第一層的低頻和高頻（A和D）系數(shù)：</p><p>　　Y=upcoef(O,X,’wname’,N,L);</p><p>　　該函數(shù)用于計(jì)算向量X向上N步的重構(gòu)小波系數(shù)。如果O=a,則是對(duì)低頻系數(shù)進(jìn)行重構(gòu)，如果O

39、=d,則是對(duì)高頻系數(shù)進(jìn)行重構(gòu)。</p><p>　　3)顯示低頻和高頻部分</p><p>　　subplot(1,2,1);plot(A);title(‘Approximation A’)</p><p>　　subplot(1,2,2);plot(D);title(‘Detail D‘)</p><p>　　4)由小波逆變換恢復(fù)信號(hào)<

40、;/p><p>　　X=idwt(cA,cD,’wname’);</p><p>　　2.2.2多尺度一維分解</p><p>　　1）wavedec為多尺度一維小波分解函數(shù)，其格式為：</p><p>　　[C,L]=wavedec(X,N,’wname’)</p><p>　　它用小波完成對(duì)信號(hào)X的一維多尺度分解。N為

41、尺度，且為嚴(yán)格的正函數(shù)。輸出參數(shù)C由[cAj,cDj,cDj-1,…,cD1]組成，L由[cAj的長度，cDj的長度，…，X的長度]組成。</p><p>　　2）提取系數(shù)的低頻和高頻部分</p><p>　　appcoef用于從小波分解結(jié)構(gòu)[C,L]中提取一維信號(hào)的低頻系數(shù)，格式為：</p><p>　　A=appcoef(C,L,’wname’,N)</p

42、><p>　　其中，[C,L]為小波分解結(jié)構(gòu)，wname為小波分解函數(shù)，N為尺度。</p><p>　　detcoef用于從小波分解結(jié)構(gòu)[C,L]中提取一維信號(hào)的高頻系數(shù)：</p><p>　　D=detcoef(C,L,N)</p><p><b>　　3)重構(gòu)信號(hào)</b></p><p>　　wr

43、coef是對(duì)一維信號(hào)的分解結(jié)構(gòu)[C,L]用指定的小波函數(shù)或重構(gòu)濾波器進(jìn)行重構(gòu)的函數(shù)：</p><p>　　X=wrcoef(‘type’,C,L,’wname’,N)</p><p>　　當(dāng)type=a時(shí)，對(duì)信號(hào)的低頻部分進(jìn)行重構(gòu)；當(dāng)type=d時(shí)，對(duì)信號(hào)的高頻部分進(jìn)行重構(gòu)。</p><p>　　4）顯示多尺度一維分解結(jié)果</p><p>&

44、lt;b>　　其方法是：</b></p><p>　　subplot(2,2,1);plot(A3);title(‘Approximation A3’)</p><p>　　subplot(2,2,2);plot(D1);title(‘Detail D1’)</p><p>　　subplot(2,2,3);plot(D2);tilte(‘Deta

45、il D2’)</p><p>　　subplot(2,2,4);plot(D3);title(‘Detail D3’)</p><p>　　圖2.3 多尺度一維分解結(jié)果</p><p><b>　　5）重構(gòu)原始信號(hào)</b></p><p>　　X=waverec(C,L,’wname’)</p><

46、p>　　這里waverec為多尺度一維小波重構(gòu)函數(shù)，用指定的小波函數(shù)對(duì)小波分解結(jié)構(gòu)[C,L]進(jìn)行重構(gòu)。</p><p>　　2.3 信號(hào)的初步去噪</p><p>　　一般來說，現(xiàn)實(shí)中的信號(hào)都是帶噪信號(hào)，在對(duì)信號(hào)做進(jìn)一步分析之前，需要將有效信號(hào)提取出來。目前，人們已根據(jù)噪聲的統(tǒng)計(jì)特征和頻譜分布規(guī)律，開發(fā)出了多種多樣的信號(hào)去噪方法。其中最為直觀的一種方法是，根據(jù)噪聲能量一般集中于高頻，

47、而信號(hào)頻譜分布于一個(gè)有限區(qū)間的特點(diǎn)，用傅里葉變換將含噪信號(hào)變換到頻域，然后采用低通濾波器進(jìn)行濾波。當(dāng)信號(hào)和噪聲的頻帶相互分離時(shí)這種方法比較有效，但當(dāng)信號(hào)和噪聲的頻帶相互重疊（比如含有白噪聲）時(shí)，則效果較差，因?yàn)榈屯V波器在抑制噪聲的同時(shí)，也將信號(hào)的邊緣部分變得模糊。</p><p>　　利用小波對(duì)信號(hào)去噪，首先要識(shí)別出信號(hào)的哪一部分或哪些部分包含噪聲，然后舍棄這些部分進(jìn)行信號(hào)重構(gòu)。對(duì)信號(hào)進(jìn)行多尺度分解，在舍棄高頻

48、分量后，信號(hào)的低頻部分變得越來越“純潔”，但同時(shí)也丟失了初始信號(hào)的瞬變特征。因此，更優(yōu)化的去噪是閾值去噪，即只去除那些超過某一設(shè)定值的細(xì)節(jié)部分。</p><p>　　2.3.1 小波分析用于消噪的基本原理</p><p>　　一個(gè)含噪的一維信號(hào)模型可表示為如下形式：</p><p>　　S(k) = f(k) + p * e(k), k = 0,1,….,n-1&

49、lt;/p><p>　　其中，s(k)為含噪信號(hào)，f(k)為有用信號(hào)，e(k)為噪聲信號(hào)。這里我們認(rèn)為e(k)是一個(gè)1級(jí)高斯白噪聲，通常表現(xiàn)為高頻信號(hào)，而工程實(shí)際中f(k)通常為低頻信號(hào)，或者是一些比較平穩(wěn)的信號(hào)。因此，我們可以按如下的方法進(jìn)行消噪處理；首先對(duì)信號(hào)進(jìn)行小波分解，一般地，噪聲信號(hào)多包含在具有較高頻率的細(xì)節(jié)中，從而可利用門限閾值等形式對(duì)所分解的小波系數(shù)進(jìn)行處理，然后對(duì)信號(hào)進(jìn)行小波重構(gòu)即可達(dá)到對(duì)信號(hào)的消噪的

50、目的。對(duì)信號(hào)消噪實(shí)質(zhì)上是抑制信號(hào)中的無用部分，恢復(fù)信號(hào)中有用部分的過程。一般地，一維信號(hào)消噪的過程可分為如下三個(gè)步驟：</p><p>　　一維信號(hào)的小波分解。選擇一個(gè)小波并確定分解的層次，然后進(jìn)行分解計(jì)算；</p><p>　　小波分解高頻系數(shù)的閾值量化。對(duì)各個(gè)分解尺度下的高頻系數(shù)選擇一個(gè)閾值進(jìn)行軟閾值量化處理；</p><p>　　一維小波的重構(gòu)。根據(jù)小波分解的

51、最底層低頻系數(shù)和各層高頻系數(shù)進(jìn)行一維小波重構(gòu)。</p><p>　　這三個(gè)步驟中，最關(guān)鍵的是如何選擇閾值及如何進(jìn)行閾值量化，在某種程度上，它關(guān)系到信號(hào)消噪的質(zhì)量。</p><p>　　2.3.2 應(yīng)用函數(shù)進(jìn)行信號(hào)消噪處理</p><p>　　消噪處理中兩個(gè)非常有用的函數(shù)wden和wdencmp。</p><p>　　1)wden的最簡單用法如

52、下：</p><p>　　sd=wden(s,tptr,sorh,scal,n,wavename)</p><p>　　它所返回的是經(jīng)過對(duì)原始信號(hào)s進(jìn)行消噪處理后的信號(hào)sd。其中，tptr指定閾值選取規(guī)則，sorh指定選取軟閾值（sorh=’s’）或硬閾值（sorh=’h’）,n為小波分解層數(shù)，wavename指定分解時(shí)所用的小波。scal是閾值尺度改變的比例，它有如下三種選擇：</

53、p><p>　?。?）scal=’one’,表示基本模式。</p><p>　　（2）scal=’sln’,表示未知尺度的基本模式，且僅根據(jù)第一層的小波分解系數(shù)來估計(jì)噪聲的層次，并只進(jìn)行一次估計(jì)，以此來變換閾值的尺度。</p><p>　?。?）scal=’mln’,表示非白噪聲的基本模式，且在每個(gè)不同的小波分解的層次上都估計(jì)噪聲的層次，以此來變換閾值的尺度。</

54、p><p>　　2)wdencmp是一種用的更為普遍的函數(shù)，它可以直接對(duì)一維或二維信號(hào)進(jìn)行消噪或壓縮，處理方法也是通過對(duì)小波分解系數(shù)進(jìn)行閾值量化來實(shí)現(xiàn)。其使用方法如下：</p><p>　　xd=wdencmp(opt,x,wavename,n,thr,sorh,keepapp)</p><p><b>　　其中：</b></p>&

55、lt;p>　　opt=’gbl’,thr>0,則閾值為全局閾值；</p><p>　　opt=’lvd’,thr是向量，則閾值是在各層上大小不同的數(shù)值。</p><p>　　keepapp=1,不對(duì)小波分解后的系數(shù)作任何處理；</p><p>　　keepapp=0,對(duì)小波分解后的低頻系數(shù)也進(jìn)行閾值量化處理。</p><p>&l

56、t;b>　　x是待處理信號(hào)。</b></p><p>　　xd是處理后的信號(hào)。其余參數(shù)同函數(shù)wden中的參數(shù)。</p><p>　　3）小波分析進(jìn)行消噪處理選取閾值一般有下述三種方法：</p><p>　　默認(rèn)閾值消噪處理。該方法利用函數(shù)ddencmp生成信號(hào)的默認(rèn)閾值，然后利用函數(shù)wdencmp進(jìn)行消噪處理。</p><p&g

57、t;　　給定閾值消噪處理。在實(shí)際的消噪處理過程中，閾值往往可以通過經(jīng)驗(yàn)公式獲得，且這種閾值比默認(rèn)閾值的可信度高。在進(jìn)行閾值量化處理時(shí)可用函數(shù)wthresh.</p><p>　　強(qiáng)制消噪處理。該方法是將小波分解結(jié)構(gòu)中的高頻系數(shù)全部置0，即濾掉所有高頻部分，然后對(duì)信號(hào)進(jìn)行小波重構(gòu)。這種方法比較簡單，且消噪后的信號(hào)比較光滑，但是容易丟失信號(hào)中的有用成分。</p><p>　　第三章 基于MAT

58、LAB的程序設(shè)計(jì)</p><p>　　本文選用光譜儀實(shí)際測(cè)得數(shù)據(jù)作為仿真對(duì)象，信號(hào)變量名為sig(圖3.1)。</p><p><b>　　圖3.1</b></p><p><b>　　3.1 設(shè)計(jì)流程</b></p><p><b>　　3.2程序設(shè)計(jì)要點(diǎn)</b></p

59、><p>　　3.2.1載入并畫出原始信號(hào)sig。</p><p>　　load spectrasig.mat; figure; </p><p>　　subplot(2,1,1); plot(sig); title('原始信號(hào)');</p><p>　　3.2.2 使用軟閾值法對(duì)信號(hào)去噪。具體思路是：</p>&l

60、t;p>　　1）選用‘sym8’小波、尺度N=10對(duì)sig進(jìn)行一維多尺度變換，得到參數(shù)矩陣C和L（見（1））。</p><p><b>　　N = 10;</b></p><p>　　w_name = 'sym8'; </p><p>　　

61、[C, L] = wavedec(sig, N, w_name); （1）</p><p>　　其中C由[cA10,cD10,cD9,…,cD1]組成，cA10代表低頻系數(shù),cD10,cD9,…,cD1代表各個(gè)尺度上的高頻系數(shù)。L由[cAj的長度，cDj的長度，…，X的長度]組成。</p><p>　　2）保持參數(shù)L不變，對(duì)系數(shù)矩陣C進(jìn)行軟閾值處理，得到cc矩

62、陣。閾值處理過程如下：對(duì)于C中前5個(gè)元素復(fù)制到cc中相應(yīng)位置；確定閾值thres(見函數(shù)wave_thres)，遍歷C中剩余元素：</p><p>　　當(dāng)C（i）小于等于thres時(shí)，cc(i)置零；</p><p>　　當(dāng)C(i)大于thres時(shí)，cc(i) = C(i) – thres；</p><p>　　當(dāng)C(i)<-thres時(shí)，cc(i) = C(

63、i) + thres。</p><p>　　3）使用waverec對(duì)系數(shù)cc和L在‘sym8’下重構(gòu)信號(hào)，命名為sig_denoise（見（2））并顯示去噪后的信號(hào)（見圖3.2）。</p><p>　　sig_denoise = waverec(cc, L, w_name);</p><p>　　subplot(2,1,2);plot(sig_denoise);&l

64、t;/p><p>　　title('去噪聲信號(hào)'); （2）</p><p>　　3.2.3 使用小波‘mexh’，在尺度為1～60下，對(duì)信號(hào)sig_denoise進(jìn)行一維連續(xù)小波變換，得到二維的系數(shù)矩陣cw_sig（見（3））。</p><p>　　cw_sig=cwt(sig_denois

65、e,1:LV,'mexh','plot') （3）</p><p>　　由小波變換的原理，一維連續(xù)小波變換是一種積分變換，小波基具有尺度α,平移τ兩個(gè)參數(shù)，因此，將函數(shù)在小波基下展開，就意味著將一個(gè)時(shí)間函數(shù)投影到二維的時(shí)間-尺度相平面上。即可得到二維的系數(shù)矩陣，其中橫坐標(biāo)上代表時(shí)域，縱坐標(biāo)代表各個(gè)尺度（見圖3.3）。矩陣中每個(gè)系數(shù)代表在相應(yīng)尺度下相應(yīng)時(shí)間點(diǎn)處變換系數(shù)

66、，這個(gè)系數(shù)表示了小波函數(shù)同該點(diǎn)處附近信號(hào)趨勢(shì)的相關(guān)性，系數(shù)值越大代表相關(guān)性越高。</p><p><b>　　圖 3.2</b></p><p>　　3.2.4 如3中所述，二維系數(shù)矩陣cw_sig中每個(gè)系數(shù)代表在相應(yīng)尺度下相應(yīng)時(shí)間點(diǎn)處的變換系數(shù)，這個(gè)系數(shù)表示了小波函數(shù)同該點(diǎn)處附近信號(hào)趨勢(shì)的相關(guān)性，系數(shù)值越大代表相關(guān)性越高。整個(gè)空間信號(hào)有若干個(gè)譜峰，對(duì)應(yīng)了單個(gè)尺度（縱

67、坐標(biāo)）上不同點(diǎn)處的若干個(gè)系數(shù)極大值；在單個(gè)點(diǎn)處，不同尺度下的變換均得到極大值系數(shù)。也就是說，在可能的譜峰位置，所有尺度上均可找到一個(gè)模極大值系數(shù)；不同尺度上系數(shù)有波動(dòng)，其中的最大值處的尺度對(duì)應(yīng)最佳變換尺度，即在此尺度上的變換可得到信號(hào)與小波函數(shù)有最高相關(guān)性。以下工作就是在各個(gè)尺度上遍歷所有點(diǎn)找到并標(biāo)記極大值系數(shù)，之后將各個(gè)極值系數(shù)對(duì)應(yīng)的標(biāo)記在各個(gè)尺度中順次按一定規(guī)則連成鏈條，鏈條即對(duì)應(yīng)可能的譜峰位置。</p><p&

68、gt;　　3.2.5建立結(jié)構(gòu)數(shù)組chain用來存放處理鏈條過程中不同類型屬性的變量empty,Mtx,ncount,a,b,c等，初始化各種變量值（見（4））。</p><p>　　for idx = 1:100</p><p>　　chain(idx).empty = 1;chain(idx).Mtx = zeros(100, 3);</p><p>　　chai

69、n(idx).ncount = 0;</p><p>　　chain(idx).a = 0;chain(idx).b = 0;chain(idx).c = 0;</p><p><b>　　end</b></p><p>　　count_chain = 0; （4）&l

70、t;/p><p><b>　　圖 3.3</b></p><p>　　3.2.6 尋找系數(shù)極大值點(diǎn)并標(biāo)志出代表峰值的模極大值系數(shù)鏈條的位置。過程如下:</p><p>　　1）在尺度2上，遍歷所有點(diǎn)上系數(shù)值，可以得到若干模極大值及它們所在位置。實(shí)現(xiàn)途徑是運(yùn)行函數(shù)[LM, IM] = my_localmax(A, th, w)，其中A = cw_si

71、g(2, :)；IM是一維數(shù)組，用來存放各個(gè)模極大值所在的空間坐標(biāo)。建立數(shù)組XM和YM，將IM中元素復(fù)制到XM中，YM中存入當(dāng)前處理的尺度即2，它們長度均等于IM長度。設(shè)定各個(gè)極值點(diǎn)是一個(gè)代表譜峰的鏈條的一部分，這樣表示鏈條數(shù)量的count_chain=Nm=length(IM)。遍歷所有極值點(diǎn)，將各個(gè)對(duì)應(yīng)坐標(biāo)，尺度和系數(shù)值存入數(shù)組chain_Mtx中（見（5））。</p><p>　　for idx = 1:N

72、m</p><p>　　chain(idx).empty = 0;</p><p>　　K = [IM(idx), 2, A(IM(idx))]; chain(idx).Mtx(1, 1:3) = K;</p><p>　　chain(idx).ncount = 1;</p><p>　　chain(idx).a = 1;chain(idx)

73、.b = 0;chain(idx).c = -IM(idx);</p><p>　　AJ(idx) = 1;BJ(idx) = 0;CJ(idx) = -IM(idx);</p><p><b>　　end</b></p><p>　　count_chain = Nm; （

74、5）</p><p>　　2）如4中所述，在可能的譜峰位置，所有尺度上均可找到一個(gè)模極大值系數(shù)；不同尺度上系數(shù)有波動(dòng)，其中的最大值處的尺度對(duì)應(yīng)最佳變換尺度。同上，利用[LM, IM] = my_localmax(A, th, w)分別找到其他所有尺度上的模極大點(diǎn)，將IM和尺度值分別存入XM和YM中。這些極大值點(diǎn)可以按規(guī)則連成若干鏈條，對(duì)應(yīng)信號(hào)中各個(gè)譜峰。</p><p>　　3）將各尺度上

75、系數(shù)極值點(diǎn)形成鏈條思路：結(jié)構(gòu)數(shù)組chain中保存了尺度2上的各個(gè)極值點(diǎn)坐標(biāo)，即Nm條直線的表達(dá)式x * AJ + y * BJ + CJ=0，其中AJ=1,BJ=0,CJ=-IM，表示chain(idx),idx=1:Nm。在尺度3上，比較各個(gè)極值點(diǎn)與Nm條直線距離，將最小值和對(duì)應(yīng)鏈條序號(hào)存入mD,iD中。設(shè)定閾值距離th_distance=13,若mD<th_distance，將該點(diǎn)歸入對(duì)應(yīng)鏈條；否則，鏈條總數(shù)count_cha

76、in加1，并將該點(diǎn)坐標(biāo)，尺度和系數(shù)值存入chain.Mtx中（見（6））。對(duì)于mD<th_distance情況下，鏈條上每增加一個(gè)點(diǎn)，chain.ncount加1；建立數(shù)組XLine和YLine，XLine存放鏈條上各個(gè)點(diǎn)橫坐標(biāo)，YLine存放尺度值；若var(XLine)<2，chain中CJ=-mean(XLine)，否則對(duì)XLine和YLine運(yùn)用最小二乘法進(jìn)行曲線擬合得到的系數(shù)p(1)和p(2)賦給AJ和CJ，BJ=

77、-1。</p><p>　　for idx = 1:Nm</p><p>　　x = IM(idx);</p><p><b>　　y = lev;</b></p><p>　　D = abs(x * AJ + y * BJ + CJ)./(sqrt(AJ.^2+BJ.^2)); </p><p>

78、;　　[mD, iD] = min(D);</p><p>　　K = [x, y, A(x)];</p><p>　　if mD < th_distance </p><p>　　chain(iD).ncount = chain(iD).ncount + 1;</p><p>　　chain(iD).Mtx(ch

79、ain(iD).ncount, 1:3) = K;</p><p>　　XLine = chain(iD).Mtx(1:chain(iD).ncount, 1);</p><p>　　YLine = chain(iD).Mtx(1:chain(iD).ncount, 2);</p><p><b>　　else</b></p>&

80、lt;p>　　count_chain = count_chain + 1;</p><p>　　chain(count_chain).empty = 0;</p><p>　　chain(count_chain).Mtx(1, 1:3) = K;</p><p>　　chain(count_chain).ncount = 1; （

81、6）</p><p>　　4）在（XM，YM）表示的所有點(diǎn)處畫出‘o’(見（7）)，可以清楚地看到代表譜峰的各個(gè)鏈條（見圖3.4）。</p><p>　　figure; plot(XM, YM, 'o'); （7）</p><p><b>　　圖 3.4</b></p><

82、;p>　　3.2.7建立數(shù)組PeakInfo記錄譜峰信息：類似同尺度上尋極值方法，在每個(gè)鏈條上找到最佳變換尺度，將坐標(biāo)、尺度值和對(duì)應(yīng)系數(shù)存入PeakInfo。</p><p>　　PeakInfo = zeros(count_chain, 3);</p><p>　　for idx = 1:count_chain</p><p>　　NJ = chain(id

83、x).ncount;</p><p>　　I_chain = chain(idx).Mtx(1:NJ, 3);</p><p>　　[LM_c, IM_c] = my_localmax(I_chain, th, w);</p><p>　　if isempty(IM_c) == 1 [LM_c, k] = max(I_chain);</p><

84、p>　　else k = IM_c(1); end</p><p>　　K = chain(idx).Mtx(k, 1:3);</p><p>　　PeakInfo(idx, 1:3) = K(1:3);</p><p><b>　　end</b></p><p>　　3.2.8 驗(yàn)證峰附近信號(hào)信噪比<

85、/p><p>　　不規(guī)則噪音疊加到信號(hào)上，如果噪音瞬時(shí)擾動(dòng)過大則產(chǎn)生新的譜峰，顯然這些譜峰是必須要排除的，思路是：對(duì)于已經(jīng)檢測(cè)到的譜峰，考察其附近信號(hào)的信噪比，若低于閾值，則認(rèn)為這個(gè)峰是由噪音擾動(dòng)形成的。由上述PeakInfo記錄了所有檢測(cè)到的譜峰信息（坐標(biāo)，尺度，系數(shù)值），對(duì)于每個(gè)峰，生成附近信號(hào)信噪比SNR_d，設(shè)定閾值3,若SNR_d<3，則放棄這個(gè)位置的譜峰。</p><p>　

86、　SNR_d = (p.a)/(r.rmse); </p><p>　　if SNR_d < 3</p><p>　　PeakFlag(idx) = 0;</p><p><b>　　else</b></p><p>　　PeakFlag(idx) = 1;</p><p>　　3.2.

87、9 畫出檢測(cè)到的譜峰（見圖3.5）。</p><p>　　for idx = 1:count_chain</p><p>　　if PeakFlag(idx) == 1</p><p>　　X = [PeakInfo(idx, 1), PeakInfo(idx, 1)];</p><p>　　Y = [0, Max_sig];</p&g

88、t;<p>　　plot(X, Y, 'blue');</p><p><b>　　end</b></p><p><b>　　end</b></p><p><b>　　圖 3.5</b></p><p><b>　　附：完整程序<

89、;/b></p><p><b>　　1）main.m</b></p><p>　　clear all; clc; close all;</p><p>　　% load original signal and plot;</p><p>　　load spectrasig.mat;</p>&l

90、t;p>　　%load simuB.mat;</p><p>　　figure; subplot(2,1,1); plot(sig); title('原始信號(hào)');</p><p><b>　　% denoise</b></p><p><b>　　N = 10;</b></p><

91、;p>　　w_name = 'sym8';</p><p>　　[C, L] = wavedec(sig, N, w_name);</p><p>　　K = floor(0.5 * N);</p><p>　　cc = wave_thres(C, L, K);</p><p>　　sig_denoise = waver

92、ec(cc, L, w_name);</p><p>　　subplot(2,1,2); plot(sig_denoise); title('去噪聲信號(hào)');</p><p>　　% cwt by mexican hat wavelet</p><p><b>　　LV = 60;</b></p><p>

93、;　　figure; cw_sig = cwt(sig_denoise, 1:LV, 'mexh', 'plot');</p><p>　　% initializee chain</p><p>　　for idx = 1:100</p><p>　　chain(idx).empty = 1;</p><p>

94、　　chain(idx).Mtx = zeros(100, 3);</p><p>　　chain(idx).ncount = 0;</p><p>　　chain(idx).a = 0;chain(idx).b = 0; chain(idx).c = 0;</p><p><b>　　end</b></p><p>　

95、　count_chain = 0;</p><p>　　% get local max,draw max lines</p><p>　　L = length(sig_denoise);</p><p>　　A = cw_sig(2, :);th = 0.01;w = 1;</p><p>　　[LM, IM] = my_localmax(A

96、, th, w);</p><p>　　Nm = length(IM);</p><p>　　XM(1:Nm) = IM(1:Nm);YM(1:Nm) = 2 * ones(1, Nm);</p><p><b>　　tm = Nm;</b></p><p>　　for idx = 1:Nm</p><

97、;p>　　chain(idx).empty = 0;</p><p>　　K = [IM(idx), 2, A(IM(idx))];</p><p>　　chain(idx).Mtx(1, 1:3) = K;chain(idx).ncount = 1;</p><p>　　chain(idx).a = 1;chain(idx).b = 0;chain(idx)

98、.c = -IM(idx);</p><p>　　AJ(idx) = 1;BJ(idx) = 0;CJ(idx) = -IM(idx);</p><p><b>　　end</b></p><p>　　count_chain = Nm;</p><p>　　th_distance = 13;</p><

99、;p>　　for lev = 3:LV</p><p>　　A = cw_sig(lev, :); w = floor(lev * 0.7);</p><p>　　[LM, IM] = my_localmax(A, th, w);</p><p>　　Nm = length(IM);</p><p>　　XM((tm+1):(tm+Nm

100、)) = IM(1:Nm);</p><p>　　YM((tm+1):(tm+Nm)) = lev;</p><p>　　tm = tm + Nm;</p><p>　　for idx = 1:Nm</p><p>　　x = IM(idx);y = lev;</p><p>　　D = abs(x * AJ + y

101、* BJ + CJ)./(sqrt(AJ.^2+BJ.^2)); </p><p>　　[mD, iD] = min(D);</p><p>　　K = [x, y, A(x)];</p><p>　　if mD < th_distance </p><p>　　chain(iD).ncount = chain(i

102、D).ncount + 1;</p><p>　　chain(iD).Mtx(chain(iD).ncount, 1:3) = K;</p><p>　　XLine = chain(iD).Mtx(1:chain(iD).ncount, 1);</p><p>　　YLine = chain(iD).Mtx(1:chain(iD).ncount, 2);</p

103、><p>　　if var(XLine) < 2</p><p>　　chain(iD).a = 1;AJ(iD) = 1;</p><p>　　chain(iD).b = 0; BJ(iD) = 0;</p><p>　　chain(iD).c = - mean(XLine);</p><p>　　CJ(iD) =

104、 -mean(XLine);</p><p><b>　　else</b></p><p>　　p = polyfit(XLine, YLine, 1);</p><p>　　chain(iD).a = p(1);AJ(iD) = p(1);</p><p>　　chain(iD).b = -1；BJ(iD) = -1;

105、</p><p>　　chain(iD).c = p(2)；CJ(iD) = p(2); </p><p><b>　　end</b></p><p><b>　　else</b></p><p>　　count_chain = count_chain + 1;<

106、/p><p>　　chain(count_chain).empty = 0;</p><p>　　chain(count_chain).Mtx(1, 1:3) = K;</p><p>　　chain(count_chain).ncount = 1;</p><p>　　chain(count_chain).a = 1;</p>&

107、lt;p>　　chain(count_chain).b = 0;</p><p>　　chain(count_chain).c = -IM(idx);</p><p>　　AJ(count_chain) = 1;</p><p>　　BJ(count_chain) = 0;</p><p>　　CJ(count_chain) = -IM

108、(idx);</p><p><b>　　end</b></p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　figure; plot(XM, YM, 'o');</p><p>　

109、　w = 2;th = 0.001;</p><p>　　PeakInfo = zeros(count_chain, 3);</p><p>　　for idx = 1:count_chain</p><p>　　NJ = chain(idx).ncount;</p><p>　　I_chain = chain(idx).Mtx(1:NJ,

110、3);</p><p>　　[LM_c, IM_c] = my_localmax(I_chain, th, w);</p><p>　　if isempty(IM_c) == 1</p><p>　　[LM_c, k] = max(I_chain);</p><p><b>　　else</b></p>&

111、lt;p>　　k = IM_c(1);</p><p><b>　　end</b></p><p>　　K = chain(idx).Mtx(k, 1:3);</p><p>　　PeakInfo(idx, 1:3) = K(1:3);</p><p><b>　　end</b></p&

112、gt;<p>　　% verify the SNR of the peak</p><p>　　for idx = 1:count_chain</p><p>　　X0 = PeakInfo(idx, 1);scale = PeakInfo(idx, 2);I0 = PeakInfo(idx, 3);</p><p>　　sigma = floor(s

113、cale/1.66) + 1;wd = 3 * sigma;</p><p>　　if X0 + wd > length(sig)</p><p>　　wd = length(sig) - X0;</p><p>　　elseif X0 - wd < 1</p><p>　　wd = X0 - 1;</p><

114、p><b>　　end</b></p><p>　　OrigData = sig((X0-wd):(X0+wd));</p><p>　　OrigData = OrigData - min(OrigData); </p><p>　　OrigData = reshape(OrigData, length(OrigData), 1);

115、 </p><p>　　XLL = 1:(2*wd + 1);XLL = XLL'; </p><p>　　Xp = wd + 1; C1 = scale/1.66; </p><p>　　f = fittype('k * x + b + a * exp(-((x-Xp)/C1)^2)', 'problem', {&

116、#39;Xp', 'C1'});</p><p>　　[p, r] = fit(XLL, OrigData, f, 'problem', {Xp, C1});</p><p>　　SNR_d = (p.a)/(r.rmse); </p><p>　　if SNR_d < 3</p><p>　

117、　PeakFlag(idx) = 0;</p><p><b>　　else</b></p><p>　　PeakFlag(idx) = 1;</p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　%

118、 plot the detected peaks</p><p>　　figure; plot(sig, 'red'); hold on; </p><p>　　Max_sig = max(sig);</p><p>　　for idx = 1:count_chain</p><p>　　if PeakFlag(idx) ==