小波域噪聲抑制報告

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

2、使時間-頻率變窄或變寬，即在低頻部分具有較高的頻率分辨率和較低的時間分辨率， 在高頻部分具有較高的時間分辨率和較低的頻率分布率， 很適合 于探測正常信號中加帶的瞬態(tài)反?，F(xiàn)象并展示其成分， 所以被譽為分析信號的顯微鏡。 原則 上講， 傳統(tǒng)上使用傅立葉分析的地方， 都可以用小波分析取代。 小波分析優(yōu)于傅立葉變換的 地方是，它在時域和頻域同時具有良好的局部化性質(zhì)。1. 連續(xù)小波變換設(shè) x L2 R ，其傅立葉變換為 ? w ，當 ? 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. 離散小波變換在實際應用中，尤其是在計算機上實現(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)相應的離散小波變換為Cf j,k R f x j,k x dx取a0 2,b0 1，則 j,k x 2j 2 2jx k ，就形成離散小波變換中最常用的二進小波。3. 二維小波變換 由于圖像是二維信號，考慮二維尺度函數(shù)是可分離的情況，也就是：(x,y) (x) (y)(4-7)其中 (x) 是一維尺度函數(shù)。 (x) 是與之相對應的小波函數(shù)， 那么就可以得到二維小波 變換的三個基本小波：1(x,y) (x) (y)2(x,y) (x) (y)(4-8)3 (x,y) (x) (y)需要注意的是這里使用的上

5、標只是索引而不是指數(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 方向分別進行一維濾波得到。于是圖像的正交小波分解可以理解為一組獨立的空間有向的頻率通道上的信號分解。 每 一尺度分解成四個子帶 LL 、HL、LH 和 HH ，分別表征圖像的低頻信息及水平、垂直和斜 方向上的細節(jié)。圖 4.1 所示為圖像的三層小波分解的結(jié)果。圖 4.1 圖像的三層小波分解示意圖從小波變換的能量觀點來看， 高斯噪聲的小波變換仍然是高斯分布， 它是均勻分

6、布在相 空間的各個部分， 而信號由于其帶限性， 它的小波變換系數(shù)僅僅集中在相空間的上的一小部 分。在小波域上， 所有的小波系數(shù)都對噪聲有貢獻， 也就是噪聲的能量分布在所有的小波系 數(shù)上， 而只有一小部分小波系數(shù)對信號能量有貢獻， 所以可以把小波系數(shù)分為兩類： 第一類 小波系數(shù)僅僅由噪聲變換得到， 這類小波系數(shù)幅值小， 數(shù)目較多； 第二類小波系數(shù)由信號變 化得到，并包含噪聲的變換結(jié)果，這類小波系數(shù)幅值大，數(shù)目較少。從這點出發(fā)，可以通過 系數(shù)幅值上的差異構(gòu)造來達到去噪的效果。小波閾值去噪分三個步驟：利用小波變換對圖像進行分解； 設(shè)置一定的閾值對小波系數(shù)進行處理； 小波逆變換得到空域結(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為圖像的長度， 為噪聲標準方差。 Donoho 提出的統(tǒng)一閾值計算簡單， 已證明在高斯白噪聲下對光滑信號是漸進最

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

9、值處理子程序流程圖：輸入圖像，將圖像像素 點數(shù)記為圖像大小 N計算噪聲標準方差 o計算去噪閾值 t輸入圖像，將圖像像 素點數(shù)記為圖像大小 N令噪聲標準 差 o=20計算去噪閾值 t利用公式 4-10進行閾值處理輸出硬閾值結(jié)果利用公式 4-11 進行閾值處理輸出軟閾值結(jié)果四、實驗程序主程序：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;五、實驗結(jié)果原圖硬