畢業(yè)論文---自動識別譜峰的程序設(shè)計

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信號的初步去噪………………………………………………………12</p><p>　　三 基于MATLAB的程序設(shè)計…………………………………………… 14</p><p>　　3.1設(shè)計流程………………………………………………………………14</

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

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

6、LAB設(shè)計出的程序很好地完成了預(yù)期目標(biāo)。</p><p>　　關(guān)鍵詞：小波變換 MATLAB 去噪 譜峰識別</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>　　在我們的周圍，每天都有大量的信號需要我們進行分析，例如我們說話的聲音，機器的振動，金融變化數(shù)據(jù)，地震信號，音樂信號，醫(yī)療圖像等。相當(dāng)多的信號需要進行有效的編碼，壓縮，消噪，重建，建模和特征提取。對于實際應(yīng)用中得到的光譜信號，我們有時需要對其進行去噪、識別、定位以及尋峰等操作。</p><p>　　本章簡要說明了幾種常用譜峰

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

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

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

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

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

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

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

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

18、間插值,L-M 算法相比較G-N 算法,更能抵制噪聲的影響,即使在初始值遠離最終最小值的時候也可以精確得到解。</p><p>　　在最小二乘法中,設(shè)給定實驗數(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 , 從而獲得高斯公式非線性曲線擬合的各個系數(shù),最終得到峰值位置。</p><p>　　7)對以上6種方法的總結(jié)</p><p>　　輸入信號的信噪比對于尋峰算法中算法誤差的影響最大。信噪比大小的改變對算法精度的影響超過相同信噪比時各種算法間精度的差異,且信噪比對所有算法的影響基本一致。6 種算法的誤差隨噪聲幅度的增大而增大

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

36、需要采用連續(xù)小波變換，但是在每個可能的尺度離散點都去計算小波系數(shù)，那將是個巨大工程，并且產(chǎn)生一大堆令人討厭的數(shù)據(jù)。如果只取這些尺度的一小部分，以及部分時間點，將會大大減輕我們的工作量，并且不失準(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為被分析的離散信號，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ù)用于計算向量X向上N步的重構(gòu)小波系數(shù)。如果O=a,則是對低頻系數(shù)進行重構(gòu)，如果O

39、=d,則是對高頻系數(shù)進行重構(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ù)信號<

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>　　它用小波完成對信號X的一維多尺度分解。N為

41、尺度，且為嚴格的正函數(shù)。輸出參數(shù)C由[cAj,cDj,cDj-1,…,cD1]組成，L由[cAj的長度，cDj的長度，…，X的長度]組成。</p><p>　　2）提取系數(shù)的低頻和高頻部分</p><p>　　appcoef用于從小波分解結(jié)構(gòu)[C,L]中提取一維信號的低頻系數(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]中提取一維信號的高頻系數(shù)：</p><p>　　D=detcoef(C,L,N)</p><p><b>　　3)重構(gòu)信號</b></p><p>　　wr

43、coef是對一維信號的分解結(jié)構(gòu)[C,L]用指定的小波函數(shù)或重構(gòu)濾波器進行重構(gòu)的函數(shù)：</p><p>　　X=wrcoef(‘type’,C,L,’wname’,N)</p><p>　　當(dāng)type=a時，對信號的低頻部分進行重構(gòu)；當(dāng)type=d時，對信號的高頻部分進行重構(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)原始信號</b></p><p>　　X=waverec(C,L,’wname’)</p><

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

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

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

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

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

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

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

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

54、p><p>　　2)wdencmp是一種用的更為普遍的函數(shù)，它可以直接對一維或二維信號進行消噪或壓縮，處理方法也是通過對小波分解系數(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,不對小波分解后的系數(shù)作任何處理；</p><p>　　keepapp=0,對小波分解后的低頻系數(shù)也進行閾值量化處理。</p><p>&l

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

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

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

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

60、t;p>　　1）選用‘sym8’小波、尺度N=10對sig進行一維多尺度變換，得到參數(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代表各個尺度上的高頻系數(shù)。L由[cAj的長度，cDj的長度，…，X的長度]組成。</p><p>　　2）保持參數(shù)L不變，對系數(shù)矩陣C進行軟閾值處理，得到cc矩

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

63、i) + thres。</p><p>　　3）使用waverec對系數(shù)cc和L在‘sym8’下重構(gòu)信號，命名為sig_denoise（見（2））并顯示去噪后的信號（見圖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('去噪聲信號'); （2）</p><p>　　3.2.3 使用小波‘mexh’，在尺度為1～60下，對信號sig_denoise進行一維連續(xù)小波變換，得到二維的系數(shù)矩陣cw_sig（見（3））。</p><p>　　cw_sig=cwt(sig_denois

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

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

67、坐標(biāo)）上不同點處的若干個系數(shù)極大值；在單個點處，不同尺度下的變換均得到極大值系數(shù)。也就是說，在可能的譜峰位置，所有尺度上均可找到一個模極大值系數(shù)；不同尺度上系數(shù)有波動，其中的最大值處的尺度對應(yīng)最佳變換尺度，即在此尺度上的變換可得到信號與小波函數(shù)有最高相關(guān)性。以下工作就是在各個尺度上遍歷所有點找到并標(biāo)記極大值系數(shù)，之后將各個極值系數(shù)對應(yīng)的標(biāo)記在各個尺度中順次按一定規(guī)則連成鏈條，鏈條即對應(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ù)極大值點并標(biāo)志出代表峰值的模極大值系數(shù)鏈條的位置。過程如下:</p><p>　　1）在尺度2上，遍歷所有點上系數(shù)值，可以得到若干模極大值及它們所在位置。實現(xiàn)途徑是運行函數(shù)[LM, IM] = my_localmax(A, th, w)，其中A = cw_si

71、g(2, :)；IM是一維數(shù)組，用來存放各個模極大值所在的空間坐標(biāo)。建立數(shù)組XM和YM，將IM中元素復(fù)制到XM中，YM中存入當(dāng)前處理的尺度即2，它們長度均等于IM長度。設(shè)定各個極值點是一個代表譜峰的鏈條的一部分，這樣表示鏈條數(shù)量的count_chain=Nm=length(IM)。遍歷所有極值點，將各個對應(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中所述，在可能的譜峰位置，所有尺度上均可找到一個模極大值系數(shù)；不同尺度上系數(shù)有波動，其中的最大值處的尺度對應(yīng)最佳變換尺度。同上，利用[LM, IM] = my_localmax(A, th, w)分別找到其他所有尺度上的模極大點，將IM和尺度值分別存入XM和YM中。這些極大值點可以按規(guī)則連成若干鏈條，對應(yīng)信號中各個譜峰。</p><p>　　3）將各尺度上

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

76、in加1，并將該點坐標(biāo)，尺度和系數(shù)值存入chain.Mtx中（見（6））。對于mD<th_distance情況下，鏈條上每增加一個點，chain.ncount加1；建立數(shù)組XLine和YLine，XLine存放鏈條上各個點橫坐標(biāo)，YLine存放尺度值；若var(XLine)<2，chain中CJ=-mean(XLine)，否則對XLine和YLine運用最小二乘法進行曲線擬合得到的系數(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）表示的所有點處畫出‘o’(見（7）)，可以清楚地看到代表譜峰的各個鏈條（見圖3.4）。</p><p>　　figure; plot(XM, YM, 'o'); （7）</p><p><b>　　圖 3.4</b></p><

82、;p>　　3.2.7建立數(shù)組PeakInfo記錄譜峰信息：類似同尺度上尋極值方法，在每個鏈條上找到最佳變換尺度，將坐標(biāo)、尺度值和對應(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 驗證峰附近信號信噪比<

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