小波域噪聲抑制報告

1、實(shí)驗四 小波域噪聲抑制一、實(shí)驗?zāi)康脑谛〔ㄗ儞Q域?qū)崿F(xiàn)圖像噪聲抑制。、算法概要在分析突變信號和非平穩(wěn)信號時， Fourier 分析顯得無能為力，因此需要尋找新的分析 工具，使之適合突變信號和非平穩(wěn)信號。小波分析正是在這一背景下產(chǎn)生的。小波分析基本思想來源于可變窗口的伸縮和平移。 Mallat 提出了實(shí)現(xiàn)小波變換的金字 塔快速算法（ Mallat 算法），使得小波變換真正進(jìn)入了實(shí)際應(yīng)用。小波變換是一種信號的時 間-頻率分析方法，它具有多分辨分析的特點(diǎn)，而且在時域和頻域都具有表征信號局部特征 的能力，是一種窗口大小固定不變但其形狀可改變而且時間窗和頻率窗都可以改變的時頻局 部化分析方法。根據(jù)高頻和低頻

2、使時間-頻率變窄或變寬，即在低頻部分具有較高的頻率分辨率和較低的時間分辨率， 在高頻部分具有較高的時間分辨率和較低的頻率分布率， 很適合 于探測正常信號中加帶的瞬態(tài)反常現(xiàn)象并展示其成分， 所以被譽(yù)為分析信號的顯微鏡。 原則 上講， 傳統(tǒng)上使用傅立葉分析的地方， 都可以用小波分析取代。 小波分析優(yōu)于傅立葉變換的 地方是，它在時域和頻域同時具有良好的局部化性質(zhì)。1. 連續(xù)小波變換設(shè) x L2 R ，其傅立葉變換為 ? w ，當(dāng) ? w 滿足允許條件（完全重構(gòu)條件或恒 等分辨條件）C ? w dw（ 4-1）Rwa,b R;a 0時，稱 x 為一個基本小波。將 x 伸縮和平移后得4-2)其中 a 為

3、伸縮因子， b為平移因子。對于任意的函數(shù)f x L2 R 的連續(xù)小波變換和逆變換為1a2 R f x xab dxa,b xab dadbWf a,b f , f x C1a,b4-3)2. 離散小波變換在實(shí)際應(yīng)用中，尤其是在計算機(jī)上實(shí)現(xiàn)時，連續(xù)小波必須加以離散化。因此，有必要討論連續(xù)小波 a,b x 和連續(xù)小波變換 Wf a,b 的離散化。 這一離散化都是針對連續(xù)尺度參數(shù) a 個連續(xù)平移參數(shù) b 的，而不是針對變量 x 的。為了使離散后的函數(shù)組能覆蓋整個a,b 所表 示的平面，取 a0 1,b0 1 ，使得4-4)a a0 ,b nb 0a0 j,k Z且將 a,b x 改記為 j,k x

4、，即j,k xa0j 2x kb0a0a0 ja0j 2 a0j x kb04-5)4-6)相應(yīng)的離散小波變換為Cf j,k R f x j,k x dx取a0 2,b0 1，則 j,k x 2j 2 2jx k ，就形成離散小波變換中最常用的二進(jìn)小波。3. 二維小波變換 由于圖像是二維信號，考慮二維尺度函數(shù)是可分離的情況，也就是：(x,y) (x) (y)(4-7)其中 (x) 是一維尺度函數(shù)。 (x) 是與之相對應(yīng)的小波函數(shù)， 那么就可以得到二維小波 變換的三個基本小波：1(x,y) (x) (y)2(x,y) (x) (y)(4-8)3 (x,y) (x) (y)需要注意的是這里使用的上

5、標(biāo)只是索引而不是指數(shù)。L2(R2) 的函數(shù) f 可以寫為如下形式f(x, y)cl,k j,k(x, y)dj,k j,k(x, y)(4-9)k k 1,2,3由 Mallat 算法可知，二維圖像的分解可以通過沿 x 方向和 y 方向分別進(jìn)行一維濾波得到。于是圖像的正交小波分解可以理解為一組獨(dú)立的空間有向的頻率通道上的信號分解。 每 一尺度分解成四個子帶 LL 、HL、LH 和 HH ，分別表征圖像的低頻信息及水平、垂直和斜 方向上的細(xì)節(jié)。圖 4.1 所示為圖像的三層小波分解的結(jié)果。圖 4.1 圖像的三層小波分解示意圖從小波變換的能量觀點(diǎn)來看， 高斯噪聲的小波變換仍然是高斯分布， 它是均勻分

6、布在相 空間的各個部分， 而信號由于其帶限性， 它的小波變換系數(shù)僅僅集中在相空間的上的一小部 分。在小波域上， 所有的小波系數(shù)都對噪聲有貢獻(xiàn)， 也就是噪聲的能量分布在所有的小波系 數(shù)上， 而只有一小部分小波系數(shù)對信號能量有貢獻(xiàn)， 所以可以把小波系數(shù)分為兩類： 第一類 小波系數(shù)僅僅由噪聲變換得到， 這類小波系數(shù)幅值小， 數(shù)目較多； 第二類小波系數(shù)由信號變 化得到，并包含噪聲的變換結(jié)果，這類小波系數(shù)幅值大，數(shù)目較少。從這點(diǎn)出發(fā)，可以通過 系數(shù)幅值上的差異構(gòu)造來達(dá)到去噪的效果。小波閾值去噪分三個步驟：利用小波變換對圖像進(jìn)行分解； 設(shè)置一定的閾值對小波系數(shù)進(jìn)行處理； 小波逆變換得到空域結(jié)果圖像。 其中

7、最關(guān)鍵的就是如何選擇閾值， 因為它控制圖 像噪聲和信息的分離。 對小波系數(shù)的閾值處理有硬門限和軟門限兩種方法， 硬門限將小于閾 值的小波系數(shù)置為零，大于閾值的保留；軟門限是基于Donoho 軟門限思想的小波系數(shù)的非線性方法。硬閾值方法：4-10)W? Th W,t W0 WW tt h 0 W t軟閾值方法W? Ts W,tsgn W W t W t（4-10 W t其中 W 為小波系數(shù)， t 2log N ，這是 Donoho 基于正交離散小波變換推倒出來的通用 閾值方法，其中 N為圖像的長度， 為噪聲標(biāo)準(zhǔn)方差。 Donoho 提出的統(tǒng)一閾值計算簡單， 已證明在高斯白噪聲下對光滑信號是漸進(jìn)最

8、優(yōu)的，故得到很廣泛的應(yīng)用。其中： x, y 是像素點(diǎn)的坐標(biāo)；a 伸縮因子；b ,k 平移因子；cl,k 小波分解低頻子帶系數(shù)；dj,k 小波分解高頻子帶系數(shù)；W 小波分解子帶系數(shù)的統(tǒng)一表示；W? 閾值化處理后的各子帶小波系數(shù)； 各子帶系數(shù)標(biāo)準(zhǔn)方差；t 降噪閾值。三、算法步驟及流程圖Step1：對原圖像進(jìn)行 3 級小波分解；Step2：通過式（ 4-10）對各子帶系數(shù)進(jìn)行閾值化處理；Step3：對閾值化處理的小波系數(shù)進(jìn)行小波逆變換。流程圖如下： 主程序流程圖：調(diào)用子程序?qū)﹂撝祷幚淼男〔ㄏ禂?shù)進(jìn)行小波逆變換分別輸出原圖像、 處理后圖像、 原 圖像直方圖和去噪后圖像直方圖硬閾值處理子程序流程圖：軟閾

9、值處理子程序流程圖：輸入圖像，將圖像像素 點(diǎn)數(shù)記為圖像大小 N計算噪聲標(biāo)準(zhǔn)方差 o計算去噪閾值 t輸入圖像，將圖像像 素點(diǎn)數(shù)記為圖像大小 N令噪聲標(biāo)準(zhǔn) 差 o=20計算去噪閾值 t利用公式 4-10進(jìn)行閾值處理輸出硬閾值結(jié)果利用公式 4-11 進(jìn)行閾值處理輸出軟閾值結(jié)果四、實(shí)驗程序主程序：tic;i=imread(SAR1.tif);i = i(:,:,1);figure;imshow(uint8(i);if isa(i, double)i = im2double(i);endiSqmag=i.2;R=1000/max(max(iSqmag);iSqmag=iSqmag*R;A1,HL1,LH

10、1,HH1=wt3det_spline(iSqmag,0);A2,HL2,LH2,HH2=wt3det_spline(A1,1);A3,HL3,LH3,HH3=wt3det_spline(A2,2);A4,HL4,LH4,HH4=wt3det_spline(A3,3);%Scale 3HH3_new=Threshold(HH3);HL3_new=Threshold(HL3);LH3_new=Threshold(LH3);%Scale 2HH2_new=Threshold(HH2);HL2_new=Threshold(HL2);LH2_new=Threshold(LH2);%Scale 1HH1

11、_new=Threshold(HH1);HL1_new=Threshold(HL1);LH1_new=Threshold(LH1);% ReconstructionA2_new=iwt3det_spline(A3,HL3_new,LH3_new,HH3_new,2);A1_new=iwt3det_spline(A2_new,HL2_new,LH2_new,HH2_new,1);6A0_new=iwt3det_spline(A1_new,HL1_new,LH1_new,HH1_new,0);idespeck=abs(A0_new/R).0.5;idespeck_r=1*(idespeck-min

12、(idespeck(:)/(max(idespeck(:)-min(idespe ck(:);%Show the denoised imagedespeckledSAR=idespeck_r*255;figure(1);imshow(i*255,); title(?- ?);figure(2);imshow(despeckledSAR,); title( ￥ ? ?);figure(3);imhist(uint8(i*255);title(?-? ?);figure(4);imhist(uint8(despeckledSAR);title(￥ ? ?);toc;硬閾值處理子程序： function current= priorshrinkG(former) a,b=size(former);N=a*b;o=std2(former); t=o*sqrt(2*log10(a*b); former(find(abs(former)t)=0; current=former;return;軟閾值處理子程序：function current=Threshold(former) a,b=size(former);N=a*b;o=20;t=o*sqrt(2*log10(N); current(find(abs(former)t)=0; return;五、實(shí)驗結(jié)果原圖硬