基于小波分析的信號去噪

Matlab與通信系統(tǒng)分析基于小波變換的信號去噪仿真院別

計算機與通信工程學院通信工程2012 年06月9日目錄概述 1小波變換 1小波變換與信號去噪 1小波分析基本理論 2小波去噪原理 2小波去噪步驟 3小波去噪的閾值 4閾值函數(shù) 4閾值的選取 5三種閾值去噪形式 6小波去噪的MATLAB實現(xiàn) 7MATLAB中小波去噪的函數(shù)集合 7小波去噪的MATLAB仿真 95.結(jié)語 13參考文獻 14人員分工 14概述小波變換J.Morlet1974A.CalderonHardy空間的原子分解和無條件基的深入研究J.O.Stromberg構(gòu)造了歷史上非常類似于現(xiàn)在的小1986S.Mallat立了構(gòu)造小波基的統(tǒng)一方法——多尺度分析之后，小波分析才開始蓬勃發(fā)展起來，之I.Daubechies的《小波十講》對小波的普及起了重要的推動作用。現(xiàn)在，小波變換已經(jīng)在科技信息產(chǎn)業(yè)領域取得了令人矚目的成就。小波變換具有Fourier變換不能解決的許多困難問題，因而小波變化被譽為“數(shù)學顯微鏡B超、CT、核磁共振成像的時間，提高分辨率等。小波變換與信號去噪效果了。1995年Donoho和Johnstone提出了小波收縮去噪的技術，他們研究的是在疊加性高1斯白噪聲環(huán)境下檢測出真實信號的情況，利用正交小波變換和高斯隨機變量的性質(zhì)對信號的小波分解系數(shù)做閾值量化，無失真的還原出真實信號。由小波變換的特性可知，高斯噪聲的小波變換仍然是高斯分布的，它均勻分布在頻率尺度空間的各部分，而信號由于其帶限性，它的小波系數(shù)僅僅集中在頻率尺度空間上的有限部分。相比于傳統(tǒng)的傅里葉變換，小波分析方法是一種窗口大?。创翱诿娣e）固定但其形狀可改變，時間窗和頻率窗都可改變的時頻局部化分析方法，它在低頻部分具有較寬的時間窗和較高的頻率分辨率，在高頻部分具有較寬的頻率窗和較高的時間分辨率。正是這種特性，使小波變換具有對信號的自適應性。本文對小波去噪方法進行了研究，驗證了小波去噪的可靠性并對比了傅里葉去噪和小波去噪的效果。小波分析基本理論小波去噪原理疊加性高斯白噪聲是最常見的噪聲模型，受到疊加性高斯白噪聲“污染”的觀測信號可以表示為：yfi i

zi

i1,., (2.1)yz~N，i i i i為噪聲水平，信號長度為n.為了從含噪信號yi

中還原出真實信號fi

，可以利用信號和噪聲在小波變換下的不同的特性，通過對小波分解系數(shù)進行處理來達到信號和噪聲分離的目的。在小波域，有效信號所產(chǎn)生的小波系數(shù)其模值往往較大；而噪聲經(jīng)過小波變換有白化趨勢（高斯白噪聲經(jīng)過小波變換，在小波域仍然表現(xiàn)為很強的隨機性，通常仍認為是高斯分布），其小波域?qū)禂?shù)模值很小。而且有用信號通常表現(xiàn)為低頻信號或是一些比較平穩(wěn)的信號，而噪聲信號則通常表現(xiàn)為高頻信號，所以我們可以先對含噪信號進行小波分解（如進行三層分解）：SCA1

CD1CA CD CD2 2 1

(2.2)CA3

CD1

CD2

CD32圖2.1三層小波分解示意圖其中CAi為分解的近似部分，CDii1,2,則噪聲部分通常包含在CD1CD2CD3中，用門限閾值對小波系數(shù)進行處理（零），最后重構(gòu)信號即可達到去噪的目的。小波去噪步驟總結(jié)去噪過程，可以分成以下三個步驟：對觀測數(shù)據(jù)作小波分解變化：WyW f Wz (2.3)0 0 0yyy1 2

，?yn

fff1 2

，?fn

z是zz1 2

，?zn

，其中用到了小波分解變換是線性變換的性質(zhì)。對小波系數(shù)W0

作門限閾值處理（根據(jù)具體情況可以使用軟閾值處理或硬閾值22lognt n

(2.4)門限閾值處理可以表示為 可以證明當n趨于無窮大時使用閾值公式(2.4)對小波tn系數(shù)作軟閾值處理可以幾乎完全去除觀測數(shù)據(jù)中的噪聲。3對處理過的小波系數(shù)作逆變換

1重構(gòu)信號：0fW

0

Wtn

d (2.5)即可得到受污染采樣信號去噪后的信號。小波去噪的閾值閾值函數(shù)在對小波系數(shù)作門限閾值處理操作時，可以使用軟閾值處理方法或硬閾值處理方法，硬閾值處理只保留較大的小波系數(shù)并將較小的小波系數(shù)置零，表達式見式(3.1)。 (,t) ,

wwt

(3.1)軟閾值處理將較小的小波系數(shù)置零但對較大的小波系數(shù)向零作了收縮，表達式見式(3.2)。wt,0, (w,t)0,Swt,

wtwt (3.2)wt直觀形式見圖3.1（圖中取t=0.4）從圖上我們可以看出軟閾值處理是一種更為平滑小波系數(shù)連續(xù)性差,在某些點上會產(chǎn)生間斷，可能會引起重構(gòu)信號的振蕩。軟閾值處理實質(zhì)上是在硬閾值的基礎上將邊界出現(xiàn)的不連續(xù)點收斂到零，對小波分解系數(shù)作了收縮，從而Donoho-Johnstone將這種去噪技術稱之為小波收縮。這樣可以有效避免間斷，獲得比較光滑的重建新號。針對硬閾值函數(shù)的不連續(xù)性和軟閾值函數(shù)存在恒定偏差,以及閾值函數(shù)中參數(shù)選擇的問題,提出了漸近半軟閾值函數(shù)。半軟閾值形式較為復雜，它將軟閾值和硬閾值的優(yōu)點結(jié)合起來，更有效地抑制了噪聲，提高了信噪比，最大限度地提高了去噪的效果。由于篇幅所限，這里不多加介紹。4硬閾值信號10.50-0.5-10 10 20 30 40 50 60 70 80 90 100軟閾值信號10.50-0.5-10 10 20 30 40 50 60 70 80 90 100閾值的選取

圖3.1硬閥值和軟閥值目前常見的閾值選取規(guī)則主要有固定閾值估計，極值閾值估計，無偏似然估計和啟發(fā)式閾值估計幾種。sqtwolog:2log(Length(2log(Length(N))其中，N為信號長度。minimaxisqtwolog量是實現(xiàn)在最壞條件下最大均方誤差最小的任選量。0, n32t[0.39360.1829log2n],n32

(3.5)rigrsure:采用史坦的無偏似然估計原理進行閾值選擇，首先得到一個給定閾值5t的思然估計，再將非似然的t最小化，就得到了所選的閾值。heursure：選擇啟發(fā)式閾值，它是sqtwolog和rigrsure的綜合，是最優(yōu)預測變量heursure采用固定閾值sqtwolog一般來講，極值閾值估計和無偏似然估計方法比較保守，當噪聲在信號的高頻段分布較少時，這兩種閾值估計方法效果較好可以將微弱的信號提取出來。而固定閾值估計和啟發(fā)式閾值估計法去噪比較徹底，在去噪時顯得更為有效，但是也容易把有用的信號誤認為噪聲而去掉。三種閾值去噪形式在實際應用中，處理信號和圖像常借助計算機工具。MATLAB相應的為我們提供了小波分析工具箱，對于利用小波變換方法對數(shù)據(jù)的處理提供了方便。利用小波分析工具箱中的函數(shù)對信號和圖像進行閾值去噪大致可分為三種：比較平滑，但容易丟失信號的有用成分。默認閾值去噪處理。該方法利用''ddencmp''用''wdencmp''函數(shù)進行去噪處理。給定軟(或硬)閾值去噪處理。在實際的去噪處理中，閾值往往可以通過經(jīng)驗公''wthresh''函數(shù)進行。在三種方法中，強制去噪處理方法較簡單，去噪后的信號較平滑，但容易丟失有用信號。對于第三種方法，閾值需要通過經(jīng)驗公式獲得。因此此方法本文也不考慮。對于初學者來說，自己設定閾值是一件較困難的事情，為此，對于閾值需要根據(jù)經(jīng)驗設定的選項不予以考慮。所以本文中采取默認閾值方法消噪。另外在利用閾值進行小波系數(shù)的處理時，可以選擇全局閾值處理和部分閾值處理。全局閾值即在每層分解的小波系數(shù)進行處理時采用相同的閾值。局部閾值即在每層分解的小波系數(shù)進行處理時根據(jù)不同層的噪聲估計值不斷調(diào)整閾值，采取不同的閾值。全局閾值對于不同特征的信號區(qū)分度差，不論信號的低頻域分量小，還是高頻分量小，都是由默認閾值來處理，有時效果會較差。比如對于含高頻分量較大，而低頻分量較小的信號，采用全局閾值算法消噪往往會丟失大量的高頻部分信息，從而引起偏差。6但是全局閾值較為簡單，且可由''ddencmp''直接獲得，若采用局部閾值處理，仍然需要按照經(jīng)驗獲取閾值向量，本文不予考慮。綜上，本文采用全局默認閾值的處理方法。小波去噪的Matlab實現(xiàn)MATLAB噪方面實現(xiàn)Donoho-Johnstone等的去噪算法，而且可以選擇使用圖形界面操作工具或者去噪函數(shù)集合兩種形式，這里利用小波去噪函數(shù)集合編寫MATLAB程序，完成仿真。4.1Matlab中小波去噪的函數(shù)集合下面是幾個最為常用的與小波去噪有關的函數(shù)：wnoise函數(shù)：調(diào)用方式：x=wnoise(fun,n,snr)；作用：產(chǎn)生Donoho-Johnstone設計的6種用于測試小波去噪效果的典型測試數(shù)據(jù)，函數(shù)根據(jù)輸入?yún)?shù)fun的值輸出名為''blocks''，''bumps''，''heavy''，''doppler''，''quadchirp''或''mishmash''的6種函數(shù)數(shù)據(jù)，數(shù)據(jù)長度為2n。信噪比為''snr''。這6種測試數(shù)據(jù)在驗證和仿真實驗時非常有用。wden函數(shù)：調(diào)用方式：[XD,CXD,LXD]=wden(X,TPTR,SORH,SCAL,N,''wname'')；[XD,CXD,LXD]=wden(C,L,TPTR,SORH,SCAL,N,''wname'')；作用：''wden''是最主要的一維小波去噪函數(shù)。其中輸入?yún)?shù)''x''為輸入需要的信號''TPTR''為3.2節(jié)中4種閥值形式，''SORH''設定為''s''表示用軟門限閾值處理表示用硬門限閾值處理3.2節(jié)中提到的4種閾值形式是基于信號水平為1的高斯白噪聲模型推導得到的，當噪聲不是白噪聲時，必須在小波分解的不同層次估計噪聲水平，SCAL=''one''不進行重新估計，SCAL=''sln''只根據(jù)第一層小波分解系數(shù)估計噪聲水平，SCAL=''mln''在每個不同的小波分解層次估計噪聲水平根據(jù)''SCAL''參數(shù)的設定''wden''函數(shù)決定最終應用于每一個小波分解層次的閾值函數(shù)。最后兩個參數(shù) ''N''和''wname''表示利用名為''wname''的小波對信號進行N層分解除''wden''外還有功能更強大的用于一維或二維小波去噪或壓縮的函數(shù)''wdencmp''。ddencmp函數(shù)：7調(diào)用方式：[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,IN2,X)；[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,''wp'',X)；[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,''wv'',X)；作用：函數(shù)''ddencmp''用于獲取信號在消噪或壓縮過程中的默認閾值。輸入?yún)?shù)''X''為一維或二維信號；''IN1''取值為''den''或''cmp''，''den''表示進行去噪，''cmp''表示進行壓縮；''IN2''取值為''wv''或''wp''，''wv''表示選擇小波，''wp''表示選擇小波包。返回值''THR''是返回的閾值；''SORH''是軟閾值或硬閾值選擇參數(shù)；''KEEPAPP''表示保存低頻信號；''CRIT''是熵名（只在選擇小波包時使用）。wdencmp函數(shù)：調(diào)用方式：[XC,CXC,LXC]=wdencmp(''gbl'',X,''wname'',N,THTR,SORH,KEEPAPP)；[XC,CXC,LXC]=wdencmp(''lvd'',X,''wname'',N,THTR,SORH)；[XC,CXC,LXC]=wdencmp(''lvd'',C,L,''wname'',N,THTR,SORH)；''wdencmp''''wname''是所用的小波函數(shù)，''gbl''(''global''的縮寫)表示每層采用不同的閾值進行處理，''N''表示小波分解的層數(shù)，''THR''為閾值向量，對于后面兩種格式，''THR''的長度為N，''SORH''表示選擇軟閾值或硬（分別取值為''s''和參數(shù)''KEEPAPP''取值為1''XC''是''XC''的小波分解結(jié)構(gòu)。MATLAB中與小波去噪有關的函數(shù)總結(jié)見表4.1和表4.2。表4.1小波分解和重構(gòu)函數(shù)分類分類函數(shù)名appcoefdetcoefdwtwavedecupcoefupwlevidwtwaverecwrcoef說明多尺度一維小波分解一維小波系數(shù)的單尺度重構(gòu)單尺度一維離散小波逆變換多尺度一維小波重構(gòu)8表4.2壓縮和消噪函數(shù)分類名稱說明閾值獲取函數(shù)ddencmp獲取在消噪和壓縮中的默認閾值thselect去噪的閾值選擇wbmpen獲取一維小波去噪閾值wdcbm用Birge-Massart算法獲取小波變換閾值去噪函數(shù)wden用小波變換對一維信號自動消噪wdencmp用小波進行消噪或壓縮閾值處理函數(shù)wthcoef一維信號小波系數(shù)的閾值處理wthresh軟閾值或硬閾值處理4.2小波去噪的MATLAB仿真以下為仿真源代碼：clc;clearsnr=4; %設置信噪比%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MATLAB中用''wnoise''產(chǎn)生測試信號%原始信號為xref，含高斯白噪聲的信號為x%信號類型為blocks（由函數(shù)中參數(shù)1決定）%長度均為2^11（由函數(shù)中的參數(shù)11決定）%信噪比snr=4（由函數(shù)中的參數(shù)snr決定）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[xref,x]=wnoise(1,11,snr);xref=xref(1:2000); %2000點x=x(1:2000); %2000點%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用全局默認閾值進行去噪處理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[thr,sorh,keepapp]=ddencmp('den','wv',x); %獲取全局默認閾值xd=wdencmp('gbl',x,'sym8',3,thr,sorh,keepapp); %利用全局默認閾值對信號去噪%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9%下面是作圖函數(shù)，作出原始信號和含噪聲信號的圖%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figuresubplot(211);plot(xref); %title('原始信號');subplot(212);plot(x);title('含噪聲信號'); %畫出含噪聲信號的圖%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面用傅里葉變換進行原信號和噪聲信號的頻譜分析%通過在頻域上的分析選擇合適的濾波器對信號進行去噪%最后將這一經(jīng)典的去噪結(jié)果與小波去噪結(jié)果進行比較%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%dt=1/(2^11); %時域分辨率Fs=1/dt; %計算頻域分辨率df=Fs/2000;xxref=fft(xref); %對原始信號做快速傅里葉變換xxref=fftshift(xxref); %將頻譜圖平移xxref=abs(xxref); %取傅里葉變換的幅值xx=fft(x); %對含噪聲信號做快速傅里葉變換xx=fftshift(xx); %將頻譜搬移absxx=abs(xx); %取傅里葉變換的幅值ff=-1000*df:df:1000*df-df; %figure(2)subplot(211);plot(ff,xxref); %title('原始信號的頻譜圖');subplot(212);plot(ff,absxx);title('含信號噪聲的頻譜圖'); %畫出含噪聲信號的頻譜圖%進行低通濾波，濾波頻率為0~200的相對頻率indd2=1:800; %0頻左邊高頻率系數(shù)置零10xx(indd2)=zeros(size(indd2));indd2=1201:2000;xx(indd2)=zeros(size(indd2)); %0頻右邊高頻系數(shù)置零xden=ifft(xx); %濾波后的信號作傅里葉逆變換xden=abs(xden); %figure(3)subplot(211);plot(xd); %title('小波去除噪后的信號');subplot(212);plot(xden); %title('傅里葉分析去噪的信號');以下是仿真結(jié)果：原始信號和含高斯白噪聲的信號波形見圖4.1。原始信號和含高斯白噪聲的信號經(jīng)4.24.3。原始信號151050-50 200 400 600 800 1000 1200 1400 1600 1800 2000含噪聲信號20100-100 200 400 600 800 1000 1200 1400 1600 1800 2000圖4.1原始信號和含高斯白噪聲的信號11100008000600040002000

原始信號的頻譜圖0-1500100005000