基于小波分析的信號(hào)去噪

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

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

zi

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

中還原出真實(shí)信號(hào)fi

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

CD1CA CD CD2 2 1

(2.2)CA3

CD1

CD2

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

，?yn

fff1 2

，?fn

z是zz1 2

，?zn

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

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

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

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

0

Wtn

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

wwt

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

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

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

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

原始信號(hào)的頻譜圖0-1500100005000