小波空域相關(guān)濾波算法分析

小波空域相關(guān)濾波算法分析摘要：小波變換已經(jīng)成為分析信號及圖像的一種重要工具，能夠描述信號在每個位置及尺度上的功率。本文首先介紹了小波空域相關(guān)濾波原理及算法框架；隨后借助Matlab工具分別利用巴特沃斯濾波、中值濾波以及空域相關(guān)濾波算法，對信噪比為19.4dB的含噪信號進(jìn)行去噪處理，結(jié)果表明空域相關(guān)濾波具有較好的去噪及保留信號突變信息的能力；最后分別對信噪比為9.1dB、18.7dB、28.8dB、40dB的含噪信號進(jìn)行分析，結(jié)果表明，空域相關(guān)濾波算法在不同信噪比下具有良好的適應(yīng)性，尤其適用于高信噪比。關(guān)鍵詞：小波變換，空域相關(guān)，濾波，信噪比1．引言小波變換從1974年提出以來，在科技信息產(chǎn)業(yè)領(lǐng)域取得了令人矚目的成就。目前小波分析的應(yīng)用領(lǐng)域越來越廣泛：數(shù)學(xué)方面的數(shù)值分析、構(gòu)造快速數(shù)值方法、曲線曲面構(gòu)造、微分方程求解、控制論等；信號分析方面的濾波、去噪、壓縮、傳遞等；圖像處理方面的圖像壓縮、分類、識別和診斷等。信號的采樣和傳輸不可避免地要受到各種噪聲和干擾的污染，因此為了獲得真實理想的信號，濾波是數(shù)據(jù)采集后的首要處理步驟。在傳統(tǒng)的基于傅里葉變換的信號去噪方法中，總是考慮使信號和噪聲的重疊部分盡可能小，在頻域通過利用不同帶寬的濾波器將信號與噪聲分開，但如果兩者重疊區(qū)域很大，或是噪聲功率過大時，就難以獲得較好的去噪效果。而小波變換以其獨特的優(yōu)勢，由于小波分析具有可變的時域、頻域分辨率，能夠?qū)π盘栠M(jìn)行精細(xì)時頻域分析的特點和優(yōu)勢，特別適用于信號的去噪。目前存在的小波濾波方法可分為貝葉斯和非貝葉斯方法，其中非貝葉斯方法又可分為三種：模極大值重構(gòu)濾波、空域相關(guān)濾波、小波域閾值濾波，三者各有所長，本文只針對空域相關(guān)濾波進(jìn)行分析。本文主要介紹了空域相關(guān)濾波的原理，對其算法步驟進(jìn)行了講解；借助Matlab工具分別利用巴特沃斯濾波、中值濾波以及空域相關(guān)濾波算法，對信噪比為19.4dB的含噪信號進(jìn)行去噪處理，結(jié)果表明空域相關(guān)濾波具有較好的去噪及保留信號突變信息的能力；最后分別對信噪比為9.1dB、18.7dB、28.8dB、40dB的含噪信號進(jìn)行分析，結(jié)果表明，空域相關(guān)濾波算法在不同信噪比下具有良好的適應(yīng)性，尤其適用于高信噪比。2．空域相關(guān)濾波原理設(shè)含噪信號的一維模型表示如下：yi=fi+zi其中y(i)為含噪信號，z(i)為噪聲，n為信號長度。為了在含噪信號中還原出真實信號f(i)，可以利用信號和噪聲在小波變換下不同的特性，通過對小波分解系數(shù)進(jìn)行處理來達(dá)到信號與噪聲分離的目的。如圖1所示，對含噪信號進(jìn)行小波分解，可知，信號的突變點具有良好的局部性質(zhì)，并且出現(xiàn)在各個尺度上，而噪聲的能量卻集中在小尺度上，其小波系數(shù)隨尺度的增加而迅速衰減，即信號和噪聲在多尺度空間上具有不同的特性，也即信號在不同尺度間具有較好的相關(guān)性，噪聲相關(guān)性較弱。根據(jù)此原理，1994年Xu設(shè)表示尺度m上位置n處含噪函數(shù)f的離散小波變換，取相鄰尺度的變換值進(jìn)行相關(guān)計算，定義Corrl(m,n)為尺度m上位置nCorrl(m,n)=i=0l-1Wm+i,n,n=1,2,?,N其中，l表示參與相關(guān)運算的尺度數(shù)，，為分解的最大尺度。由于信號突變部分的寬度隨尺度的增大而增大，當(dāng)l取值較大時會使相鄰?fù)蛔凕c相互影響，因此在實際應(yīng)用中一般取l=2。圖2給出了圖1中信號的小波系數(shù)及相關(guān)系數(shù)，可知相關(guān)系數(shù)極好的去除了噪聲的影響，并保留了初始信號。圖1模擬原信號與噪聲信號的離散二階小波變換圖2a)小波系數(shù)，b)相關(guān)系數(shù)為了使小波系數(shù)與相關(guān)系數(shù)具有可比性，將歸一化到的能量上，得到，計算公式為：（3）（4）（5）通過比較的關(guān)系可鑒別信號的重要邊緣，若，則認(rèn)為n點處的小波變換值由信號產(chǎn)生，儲存此處的的位置及大小，并將及置零，將剩余的數(shù)據(jù)再次進(jìn)行歸一化計算，大小比較，重復(fù)此過程，直到中保留下的點的能量（即）滿足某一噪聲能量閾值。隨后用剩余的點進(jìn)行小波重構(gòu)即可得到濾波后的信號。3．空域相關(guān)濾波的算法框架空域相關(guān)濾波的算法框架如圖3所示，可分為以下六個步驟（取為濾波后的值，初始值為零）：對含噪信號進(jìn)行小波變換，得到；取相鄰尺度的變換值進(jìn)行相關(guān)計算；為方便比較，將歸一化到的能量上，計算公式為：若，則認(rèn)為n點處的小波變換值由信號產(chǎn)生，將此處的賦予，并將及置零；否則認(rèn)為此處的小波變換值由噪聲產(chǎn)生，保留；重復(fù)步驟（3）與步驟（4），直到滿足一定的噪聲門限，此時中保留了去除噪聲后的小波變換；進(jìn)行小波逆變換，得到濾波后的信號。

YNN

Y圖3算法框架4．實驗分析本文通過對仿真含噪一元信號進(jìn)行去噪，說明了上述濾波算法的應(yīng)用價值，并通過與不同濾波器濾波效果進(jìn)行比較，說明了空域濾波在具有突變信號檢測中的優(yōu)勢。原信號與含噪信號如圖4所示，原信號具有較明顯的突變，噪聲為高斯白噪聲,信號含噪比SNR為19.4dB。圖4原信號與噪聲信號圖5為利用巴特沃斯2階低通濾波器濾波得到的數(shù)據(jù)，上升沿與下降沿的高頻信號被濾除一部分，在突變點處信號有嚴(yán)重失真。圖5巴特沃斯2階低通濾波結(jié)果圖6中值濾波結(jié)果圖6為中值濾波器濾波后得到的數(shù)據(jù)，可知中值濾波很好地保存了突變點處的高頻信息，相較巴特沃斯濾波器更適合檢測具有突變點的信號，從而能夠適用于圖像邊緣的檢測。圖7為利用小波空域處理前后的高頻系數(shù)，可知空域相關(guān)處理后，信號產(chǎn)生的小波系數(shù)得到保留，而噪聲產(chǎn)生的系數(shù)則極大地消除。圖8為空域相關(guān)濾波后重構(gòu)得到的信號，噪聲極大消弱，同時很大程度上保留了信號信息。a)處理前b)處理后圖7圖8空域相關(guān)濾波結(jié)果通過以上濾波結(jié)果對比可知，空域相關(guān)濾波相較巴特沃斯濾波與中值濾波，在保留信號突變點處的高頻信息的同時，還能夠更大程度的去除噪聲的影響。且此突變越大，相對信噪比越高，處理結(jié)果越好，也由此可見空域相關(guān)濾波比較適用于高信噪比的信號處理。如圖9所示給出了不同信噪比下原信號與濾波結(jié)果之間的對比情況，也說明空域相關(guān)濾波在高信噪比下更為實用。SNR=9.1dBSNR=18.7dBSNR=28.8dBSNR=40dB圖9不同信噪比下空域相關(guān)濾波結(jié)果5．總結(jié)本文詳細(xì)介紹了空域相關(guān)濾波算法，通過實驗證明了空域相關(guān)濾波算法的可行性，并通過對不同信噪比的含噪信號進(jìn)行濾波分析，表明了空域相關(guān)濾波算法在不同信噪比下具有良好的適應(yīng)性，尤其適用于高信噪比。參考文獻(xiàn):[1]YansunXu,etal.WaveletTransformDomainFilters:ASpatiallySelectiveNoiseFiltrationTechnique.IEEETrans.ImageProccessing,1994,3(6):747~758[2]黃宜軍，汪金友，小波分析在微弱信號測量中的應(yīng)用研究[J]，計量學(xué)報，2007，28（2），PP.164~166[3]王新樓，小波去噪方法分析與Matlab仿真[J]，工業(yè)控制計算機(jī)，2008，21（6），pp.55~56[4]DonohoDL．De—noisingbysoft—thresholding[J]．IEEETransOnIT，1995，41：613~627．[5]（美）DwightF.Mix，KraigJ.Olejniczak[著]，楊志華楊力華[譯]，小波基礎(chǔ)及應(yīng)用教程[M]：北京：機(jī)械工業(yè)出版社，2006.4[6]陳曉楠，索繼東.空域相關(guān)濾波法在突變檢測中的應(yīng)用.大連海事大學(xué)學(xué)報，2007(33):139~141[7]潘泉等.小波濾波方法及應(yīng)用.北京：清華大學(xué)出版社，2005附件：load('noise.mat');%載入噪聲信號I=zeros(256,1);%取信號為256點，初始化I(41:216)=1;I(70:90)=1.5;I(170:185)=1.2;noise=wgn(256,1,-20);signal=I+noise;SNR=20*log10(sqrt(sum(I.^2)/sum(noise.^2)));points=256;level=5;wf='bior1.5';%sym8,bior1.5[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);[swa,swd]=swt(signal,level,Lo_D,Hi_D);%swd是細(xì)節(jié)系數(shù)，swa是近似系數(shù)Swd_n=swd;Swd_n=swd;swd_org=swd;mask_n=zeros(size(Swd_n));%先把系數(shù)處理矩陣設(shè)置為全0。forj=1:(level-1)%在1:(level-1)分解層次上對高頻系數(shù)處理，最后一層無法求相關(guān)系數(shù)，所以不作處理。Noise_d1=swd_org(j,:);Noise_d1=Noise_d1(1:20);Noise_var=var(Noise_d1);%以信號的前20個只含有噪聲的點估計噪聲在各層的方差。Pw_var=var(swd_org(j,:));Corr=swd_org(j,:).*swd_org(j+1,:);%定義相關(guān)系數(shù)為相鄰兩層的乘積。cc=1.7;%用以設(shè)定停止迭代的噪聲能量閾值，需要根據(jù)情況調(diào)節(jié)%whilePw_var>cc*Noise_varPw=sum(abs(swd(j,:)).^2);%計算小波能量Pcorr=sum(abs(Corr).^2);%計算相關(guān)系數(shù)能量Corr_new=Corr.*((Pw/Pcorr)^0.5);%歸一化corr_mod=abs(Corr_new);w_mod=abs(swd(j,:));swd_n=swd(j,:).*(corr_mod>w_mod);%(corr_mod>w_mod)返回0或者1swd_n1=(swd_n~=0);mask_n(j,:)=mask_n(j,:)+swd_n1;%將選出的點賦給系數(shù)處理矩陣相應(yīng)位置。swd_n0=ones(size(swd_n1));swd_n0=swd_n0-swd_n1;swd(j,:)=swd(j,:).*swd_n0;%將高頻系數(shù)選出大值后的地方置0。Pw_var=var(swd(j,:));Corr_new=Corr_new.*swd_n0;%將相關(guān)系數(shù)選出大值后的地方置0。Corr=Corr_new;endendmask_max=ones(1,length(mask_n));mask_n=[mask_n((1:(level-1)),:);mask_max];%最后一層系數(shù)處理矩陣全置1。Swd_reg=swd_org.*mask_n;signal_n=iswt(swa,Swd_reg,wf);%S_mix=wden(signal_n,'sqtwolog','s','sln',5,'sym8');%rigrsure;heursure;sqtwolog;minimaxi%%%%%%%%%%%畫圖：figure;%空域法處理后的高頻系數(shù)。subplot(level+1,1,1);plot(real(signal));gridon;axistight;title('空域法處理后的高頻系數(shù)');fori=1:levelsubplot(level+1,1,i+1);plot(Swd_reg(i,:));axistight;gridon;ylabel(strcat('j=',num2str(i)));endfigure;%高頻系數(shù)處理前后的比較。fori=1:levelsubplot(level,2,2*(i)-1);plot(swd_org(i,:));axistight;gridon;ylabel(strcat('d',num2str(i)));subplot(level,2,2*(i));plot(Swd_reg(i,:));axistight;gridon;ylabel(strcat('j=