《MATLAB應(yīng)用圖像處理》課件第11章

第11章MATLAB小波分析在圖像處理中的應(yīng)用

11.1小波和小波變換

11.2離散小波變換和Mallat算法

11.3MATLAB小波分析在圖像處理中的應(yīng)用習(xí)題

11.1小波和小波變換

小波(Wavelet)這一術(shù)語(yǔ)，顧名思義，“小波”就是小區(qū)域、長(zhǎng)度有限、均值為0的波形。所謂“小”是指它具有衰減性，而稱之為“波”則是指它的波動(dòng)性，其振幅正負(fù)相間的振蕩形式。小波是定義在有限間隔且其平均值為零的一種函數(shù)，它的波形如圖11.1所示。圖11.1是從許多使用比較廣泛的小波中挑選出的幾種一維小波。圖11.1所示的小波中，縮放函數(shù)和小波函數(shù)的名稱大多數(shù)是以開(kāi)發(fā)者的名字命名的，例如Moret小波函數(shù)是Grossmann和Morlet在1984年開(kāi)發(fā)的，db6縮放函數(shù)和db6小波函數(shù)是Daubechies開(kāi)發(fā)的幾種小波之一，Meyer縮放函數(shù)和Meyer小波函數(shù)是Meyer開(kāi)發(fā)的。但也有不少例外，例如sym6縮放函數(shù)和sym6小波函數(shù)則是symlets的簡(jiǎn)寫(xiě)，是Daubechies提議開(kāi)發(fā)的幾種對(duì)稱小波之一，coif2縮放函數(shù)和coif2小波函數(shù)是Daubechies應(yīng)R.Coifman的請(qǐng)求而開(kāi)發(fā)的幾種小波之一。圖11.1小波波形圖11.1所示的小波具有有限的持續(xù)時(shí)間和突變的頻率和振幅，波形可以是不規(guī)則的，也可以是不對(duì)稱的，在整個(gè)時(shí)間范圍里的幅度平均值為零。與傅立葉分析相比，小波分析能更有效地處理不規(guī)則的激烈變化信號(hào)。傅立葉變換是把一個(gè)信號(hào)分解成各種不同頻率的正弦波，因此正弦波是傅立葉變換的基函數(shù)。同樣，小波分析是把一個(gè)信號(hào)分解成將原始小波經(jīng)過(guò)移位和縮放之后的一系列小波，因此小波同樣可以用來(lái)表示一些函數(shù)的基函數(shù)。設(shè)x(t)是平方可積函數(shù)，即x(t)∈L2(R)，則x(t)的小波變換(WaveletTransform，WT)定義為其中，ψ(t)是基本小波或母小波(MotherWavelet)函數(shù)，*為的共軛，而是基本小波ψ(t)的位移和尺度伸縮，也稱為ψ(t)的生成小波。式中a＞0，稱為尺度因子，b反映時(shí)間位移，其值可正可負(fù)，符號(hào)<·>表示內(nèi)積。(11.1)式(11.1)中的t、a和b均為連續(xù)變量，因此式(11.1)稱為連續(xù)的小波變換(ContinuousWaveletTransform，CWT)。尺度因子a的作用是對(duì)基本小波ψ(t)作伸縮，a愈大ψ(t/a)愈寬。對(duì)于一個(gè)持續(xù)時(shí)間有限的小波，ψ(t)與ψab(t)之間的關(guān)系以及不同尺度a下小波分析區(qū)間的變化可用圖11.2表示。從圖11.2中可以看出，小波的持續(xù)時(shí)間隨a的增大而加寬，幅度則與成反比，但波形保持不變。圖11.2小波的移位和伸縮

ψab(t)前加因子1/的目的是使在不同的a值下ψab(t)的能量保持不變。假定小波ψ(t)是高斯型的，即Morlet小波：

(11.2)它的頻譜Ψ(ω)為(11.3)由于(ω)是中心頻率在ω0處的高斯型函數(shù)，如圖11.3(a)所示，形狀近似于矩形函數(shù)，因此可以表征X(ω)在ω0附近的局部性質(zhì)，X(ω)為x(t)的頻譜。如果采用不同的尺度伸縮因子a，(aω)的中心頻率和帶寬將發(fā)生變化。例如當(dāng)a?=?2時(shí)，ψ(t/2)的傅立葉變換為(11.4)可見(jiàn)此時(shí)中心頻率降到ω0/2，而帶寬的比例系數(shù)由2T-1/2變?yōu)門-1/2，如圖11.3(b)所示，因而|(aω)|的品質(zhì)因數(shù)不變。圖11.3小波中心頻率和帶寬隨尺度的變化總之，從頻域上看，用不同的尺度作小波變換，相當(dāng)于用一組中心頻率不同的帶通濾波器對(duì)信號(hào)進(jìn)行處理。帶通濾波器既可以對(duì)信號(hào)進(jìn)行分解，也可以用于信號(hào)的檢測(cè)。圖11.4表示小波變換在時(shí)頻平面上的基本分析單元的特點(diǎn)。當(dāng)a值小時(shí)，ψ(t/a)很“窄”，因此在時(shí)間軸上的觀測(cè)范圍小，可以“細(xì)致觀察”時(shí)域波形的變化；而在頻域上相當(dāng)于用較高頻率的小波，對(duì)信號(hào)的頻譜作分辨率較低的分析。當(dāng)a值較大時(shí)，ψ(t/a)變“寬”，時(shí)間軸上的觀測(cè)范圍大，可以“初略觀察”時(shí)域波形；而在頻域上相當(dāng)于用低頻小波對(duì)信號(hào)作分辨率較高的分析。分析頻率有高有低，但在各分析頻段內(nèi)的品質(zhì)因數(shù)Q卻保持恒定。圖11.4小波變換的時(shí)頻特性這種分析特點(diǎn)是工程實(shí)際所期望的。對(duì)于高頻信號(hào)，我們希望在時(shí)域有較高的分辨率，而在頻域上的分辨率則允許相應(yīng)地降低，因此就要求用窄的小波(即用小的尺度因子a)來(lái)“仔細(xì)觀測(cè)”時(shí)域波形x(t)。反之，對(duì)于低頻信號(hào)，我們希望提高頻域上的分辨率，而對(duì)時(shí)域的分辨率則可以降低要求，此時(shí)，就應(yīng)當(dāng)用具有窄帶頻譜的小波(即用較大的尺度因子a)來(lái)觀察X(ω)。如上所述，小波分析恰恰具有自動(dòng)調(diào)整時(shí)域和頻域的“視野”和分析頻率高低的功能，從而保證了各分析頻段內(nèi)的品質(zhì)因數(shù)Q的不變性。11.2離散小波變換和Mallat算法

在用計(jì)算機(jī)進(jìn)行小波分析時(shí)，實(shí)際上是用離散的數(shù)據(jù)進(jìn)行計(jì)算的，希望提高計(jì)算效率。不難想象，連續(xù)小波變換的計(jì)算量是驚人的。為了解決計(jì)算量的問(wèn)題，縮放因子和平移參數(shù)都選擇2的整數(shù)倍。使用這樣的縮放因子和平移參數(shù)的小波變換叫做雙尺度小波變換，它是離散小波變換(DiscreteWaveletTransform，DWT)的一種形式?？傮w來(lái)說(shuō)，離散小波變換通常指的是雙尺度小波變換。執(zhí)行離散小波變換的有效方法是使用濾波器。該方法是Mallat在1988年開(kāi)發(fā)的，叫做Mallat算法，這種方法實(shí)際上是一種信號(hào)的分解方法，在數(shù)字信號(hào)處理中稱為雙通道子帶編碼。雙通道子帶編碼主要涉及濾波器族理論、尺度分析和子帶編碼算法。子帶編碼算法使用濾波器實(shí)現(xiàn)小波變換的算法，如圖11.5所示。圖11.5雙通道編碼濾波算法圖11.5中，原始的輸入信號(hào)通過(guò)兩個(gè)互補(bǔ)的濾波器產(chǎn)生兩個(gè)信號(hào)：大尺度低頻分量和小尺度高頻分量。在小波分析中，低頻分量是大尺度的縮放因子產(chǎn)生的系數(shù)，而高頻分量是小的縮放因子產(chǎn)生的系數(shù)。低頻分量表示信號(hào)的近似值，高頻分量表示信號(hào)的細(xì)節(jié)。圖11.5只表示了在一個(gè)分辨率下的小波分解。當(dāng)需要進(jìn)行多分辨率分解時(shí)，離散小波變換可以被表示成由低通濾波器和高通濾波器組成的一棵樹(shù)。信號(hào)的分解過(guò)程可以迭代，也就是說(shuō)可進(jìn)行多級(jí)分解。如果對(duì)信號(hào)的高頻分量不再分解，而對(duì)低頻分量連續(xù)進(jìn)行分解，就得到許多分辨率較低的低頻分量，形成如圖11.6所示的一棵樹(shù)，這種樹(shù)叫做小波分解樹(shù)。分解級(jí)數(shù)的多少取決于待分析的數(shù)據(jù)和用戶的需要。圖11.6多級(jí)小波分解樹(shù)小波分解樹(shù)表示只對(duì)信號(hào)的低頻分量進(jìn)行連續(xù)分解。如果不僅對(duì)信號(hào)的低頻分量進(jìn)行連續(xù)分解，而且對(duì)高頻分量也進(jìn)行連續(xù)分解，這樣不僅可得到許多分辨率較低的低頻分量，而且也可得到許多分辨率較低的高頻分量。這樣分解得到的樹(shù)叫做小波包分解樹(shù)(waveletpacketdecompositiontree)，這種樹(shù)是一個(gè)完整的二進(jìn)制樹(shù)。圖11.7表示的是一棵三級(jí)小波包分解樹(shù)。小波包分解方法是小波分解的一般化，可為信號(hào)分析提供更豐富和更詳細(xì)的信息。圖11.7三級(jí)小波包分解樹(shù)在使用濾波器對(duì)數(shù)字信號(hào)進(jìn)行變換時(shí)，得到的數(shù)據(jù)量將是原始數(shù)據(jù)量的兩倍。例如，如果原始信號(hào)的數(shù)據(jù)樣本為1000個(gè)，通過(guò)濾波之后每一個(gè)通道的數(shù)據(jù)均為1000個(gè)，總共為2000個(gè)。于是，根據(jù)奈奎斯特(Nyquist)采樣定理可以采用降采樣(downsampling)的方法來(lái)解決問(wèn)題，即在每個(gè)通道中每?jī)蓚€(gè)樣本數(shù)據(jù)中取一個(gè)，如圖11.8所示。圖11.8降采樣過(guò)程

MATLAB中用wavedec函數(shù)實(shí)現(xiàn)多級(jí)離散小波分解，wavedec的調(diào)用格式如下：

[C,L]=wavedec(X,N,'wname')

[C,L]=wavedec(X,N,Lo_D,Hi_D)其中X為輸入信號(hào)，N為分解級(jí)數(shù)，wname為小波名稱，Lo_D為低通濾波器，Lo_H為高通濾波器。返回的參數(shù)C為各級(jí)小波參數(shù)，L為各級(jí)小波參數(shù)的數(shù)據(jù)長(zhǎng)度。圖11.9為wavedec小波分解示意圖。圖11.9wavedec小波分解下面一段代碼使用db2小波，即Daubechies小波，對(duì)一個(gè)正弦疊加信號(hào)進(jìn)行3級(jí)小波分解。圖11.10描述了信號(hào)的fft分解和小波分解的結(jié)果。

%調(diào)入正弦疊加信號(hào)

loadsumsin;s=sumsin;

%使用db2小波對(duì)信號(hào)進(jìn)行3級(jí)分解

[C,L]=wavedec(s,3,'db2');

subplot(3,1,1);

plot(s);title('輸入信號(hào)');subplot(3,1,2);plot(abs(fft(s)))；title('輸入信號(hào)fft變換');subplot(3,1,3);plot(C);title('小波系數(shù)');axis([01000-55]);

L中的各級(jí)小波系數(shù)個(gè)數(shù)：

L=1271272525011000信號(hào)的小波重構(gòu)是小波分解的逆過(guò)程，MATLAB中用idwt實(shí)現(xiàn)信號(hào)的小波重構(gòu)。idwt的調(diào)用格式如下：

X=idwt(cA,cD,'?wname')

X=idwt(cA,cD,Lo_R,Hi_R)

X=idwt(cA,cD,'wname',L)

X=idwt(cA,cD,Lo_R,Hi_R,L)

圖11.10wavedec3級(jí)小波分解

idwt使用的算法如圖11.11所示，其中2倍上采樣是通過(guò)數(shù)據(jù)間插0實(shí)現(xiàn)的，低通和高通濾波器輸出的數(shù)值相加后取中間的L個(gè)數(shù)值得到j(luò)?-?1階小波系數(shù)，L默認(rèn)為cAj(或cDj)的長(zhǎng)度，也可以自由設(shè)定。圖11.11小波重構(gòu)算法

idwt實(shí)現(xiàn)一級(jí)小波系數(shù)的重構(gòu)，多次調(diào)用idwt可實(shí)現(xiàn)多級(jí)小波系數(shù)的重構(gòu)。下面一段代碼通過(guò)3次調(diào)用idwt并將1級(jí)細(xì)節(jié)小波系數(shù)清零，實(shí)現(xiàn)信號(hào)重構(gòu)，計(jì)算結(jié)果如圖11.12所示。由于清除了1級(jí)細(xì)節(jié)小波系數(shù)，圖11.12中的重構(gòu)信號(hào)比輸入信號(hào)更為光滑。

%調(diào)入正弦疊加信號(hào)

loadsumsin;s=sumsin;

%使用db2小波對(duì)信號(hào)進(jìn)行3級(jí)分解

[C,L]=wavedec(s,3,'db2');%構(gòu)建3級(jí)小波系數(shù)cA3=C(1:1:L(1));cD3=C(L(1)+1:1:L(1)+L(2));%構(gòu)建2級(jí)小波系數(shù)cA2=idwt(cA3,cD3,'db2');cD2=C(L(1)+L(2)+1:1:L(1)+L(2)+L(3));%構(gòu)建1級(jí)小波系數(shù)cA1=idwt(cA2,cD2,'db2',L(5));%將1級(jí)細(xì)節(jié)小波系數(shù)清零cD1=zeros(1,L(4));%重構(gòu)信號(hào)ss=idwt(cA1,cD1,'db2',L(5));subplot(2,1,1);plot(s);title('輸入信號(hào)');subplot(2,1,2);plot(ss)title('重構(gòu)信號(hào)');axis([01000-55]);圖11.12小波信號(hào)重構(gòu)11.3MATLAB小波分析在圖像處理

中的應(yīng)用

11.3.1小波圖像分解與重構(gòu)

1.二維小波分解

MATLAB中實(shí)現(xiàn)了單尺度二維離散小波分解、多尺度二維離散小波分解、提取二維小波分解的高頻系數(shù)以及提取二維小波分解的低頻系數(shù)。

1)單尺度二維離散小波分解

MATLAB中實(shí)現(xiàn)單尺度二維離散小波分解的函數(shù)是dwt2，其調(diào)用格式有以下兩種：

(1)?[cA,cH,cV,cD]=dwt2(X,'wname')；

(2)[cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D)。其中X是被分解的離散信號(hào)，wname為分解小波函數(shù)，Lo_D和Hi_D分別為分解低通濾波器和高通濾波器，兩者長(zhǎng)度必須相等。返回值cA、cH、cV和cD分別為低頻系數(shù)和高頻系數(shù)向量?？梢园凑崭袷絒cA,cH,cV,cD]=dwt2(...,'mode',MODE)，使用擴(kuò)展模式進(jìn)行信號(hào)的小波分解，參數(shù)MODE可以參考函數(shù)dwtmode。如果sx=size(X)，lf為濾波器的長(zhǎng)度，則size(cA)=size(cH)=size(cV)=size(cD)=sa，如果DWT的擴(kuò)展模式為周期，此時(shí)sa=ceil(sx/2)，對(duì)于其它的模式，此時(shí)sa=floor((sx+lf-1)/2)。

2)多尺度二維離散小波分解

MATLAB中實(shí)現(xiàn)多尺度二維離散小波分解的函數(shù)是wavedec2，其調(diào)用格式有以下兩種：

(1)?[C,S]=wavedec2(X,N,'wname')；

(2)?[C,S]=wavedec2(X,N,Lo_D,Hi_D)。其中N必須為嚴(yán)格的正整數(shù)，wname為小波函數(shù)，Lo_D和Hi_D分別為分解低通濾波器和高通濾波器。返回分解結(jié)構(gòu)[C,S]，其中C的結(jié)構(gòu)為

C=[A(N)|H(N)|V(N)|D(N)|...H(N-1)|V(N-1)|D(N-1)|...|H(1)|V(1)|D(1)]式中A(N)為低頻系數(shù)，H(N)為水平高頻系數(shù)，V(N)為垂直高頻系數(shù)，D(N)為斜線高頻系數(shù)，所有向量均以列向量存儲(chǔ)在矩陣C中。矩陣S的結(jié)構(gòu)如下：

S(1,:)=尺度N的低頻系數(shù)矩陣的大??；

S(i,:)=尺度N-i+2的高頻系數(shù)矩陣的大小，i=2,…，N+1；

S(N+2,:)=size(X)。其分解過(guò)程如圖11.13所示。圖11.13信號(hào)的多尺度二維小波分解

3)提取二維小波分解的高頻系數(shù)

MATLAB中實(shí)現(xiàn)提取二維小波分解的高頻系數(shù)的函數(shù)是detcoef2，其調(diào)用格式如下：

D=detcoef2(O,C,S,N)其中[C,S]是信號(hào)的二維小波分解結(jié)構(gòu)，尺度N必須為正整數(shù)且1≤N≤size(S,1)?-?2，O取值為'h'('v'或'd')，當(dāng)O為?'h'?時(shí)，表示提取水平方向系數(shù)，當(dāng)O為?'v'?時(shí)，表示提取垂直方向系數(shù)，當(dāng)O為?'d'?時(shí)，表示提取斜線方向系數(shù)。

4)提取二維小波分解的低頻系數(shù)MATLAB中實(shí)現(xiàn)提取二維小波分解的低頻系數(shù)的函數(shù)是appcoef2，其調(diào)用格式有以下四種：

(1)?A=appcoef2(C,S,'wname',N)；

(2)?A=appcoef2(C,S,'wname')；

(3)?A=appcoef2(C,S,Lo_R,Hi_R)；

(4)?A=appcoef2(C,S,Lo_R,Hi_R,N)。

appcoef2與detcoef2對(duì)應(yīng)，用于提取二維小波分解的低頻系數(shù)。wname為小波函數(shù)，Lo_R和Hi_R分別為低通濾波器和高通濾波器。格式(1)和(4)提取尺度N(0≤N≤size(S,1)-2)的低頻系數(shù)，格式(2)和(3)提取最后一尺度(N=size(S,1)-2)的低頻系數(shù)。

2.圖像的小波分解實(shí)例下面通過(guò)例11.1和11.2說(shuō)明如何對(duì)圖像進(jìn)行單尺度分解和多尺度分解并提取多尺度分解的小波系數(shù)。例11.1在下列程序代碼中，調(diào)用函數(shù)dwt2對(duì)圖像進(jìn)行單尺度分解，并畫(huà)出圖像分解的低頻分量和水平、垂直及斜線方向的三個(gè)高頻分量。

%加載圖像

loadtire;

nbc=size(map,1);colormap(pink(nbc));%用小波函數(shù)db2對(duì)信號(hào)x進(jìn)行單層分解[cA1,cH1,cV1,cD1]=dwt2(X,'db2');figure(1);subplot(3,2,1);image(wcodemat(X,nbc));title('原始圖像');subplot(3,2,3);image(wcodemat(cA1,nbc));title('低頻系數(shù)圖像');subplot(3,2,4);image(wcodemat(cH1,nbc));title('水平高頻系數(shù)圖像');subplot(3,2,5);image(wcodemat(cV1,nbc));title('垂直高頻系數(shù)圖像');subplot(3,2,6);image(wcodemat(cD1,nbc));title('斜線高頻系數(shù)圖像');程序運(yùn)行結(jié)果如圖11.14所示?？梢钥闯龅皖l分量表現(xiàn)了圖像的輪廓，而高頻分量表現(xiàn)了圖像的細(xì)節(jié)。圖11.14單尺度圖像分解例11.2下列程序代碼對(duì)圖像進(jìn)行多尺度分解并提取圖像分解的低頻和高頻小波系數(shù)。程序中首先利用函數(shù)wavedec2對(duì)圖像進(jìn)行2層分解，然后利用函數(shù)detcoef2提取小波分解的高頻系數(shù)，再利用函數(shù)appcoef2提取小波分解的低頻系數(shù)。

loadtire;

%用小波函數(shù)db1對(duì)信號(hào)進(jìn)行2尺度分解

[c,s]=wavedec2(X,2,'db1');

sizex=size(X)；

sizec=size(c)；val_s=s%提取尺度2的所有方向的高頻系數(shù)[chd2,cvd2,cdd2]=detcoef2('all',c,s,2);sizecd2=size(chd2)；%提取尺度1的所有方向的高頻系數(shù)[chd1,cvd1,cdd1]=detcoef2('all',c,s,1);sizecd1=size(chd1)；%提取尺度2的低頻系數(shù)ca2=appcoef2(c,s,'db1',2);sizeca2=size(ca2)；%提取尺度1的低頻系數(shù)ca1=appcoef2(c,s,'db1',1);sizeca1=size(ca1)；nbc=size(map,1);colormap(pink(nbc));figure(1);subplot(221);image(wcodemat(X,nbc));title('原始圖像');subplot(223);image(wcodemat(ca2,nbc));title('尺度2的低頻系數(shù)');subplot(224);image(wcodemat(ca1,nbc));title('尺度1的低頻系數(shù)');figure(2);subplot(231);image(wcodemat(chd2,nbc));title('尺度2水平方向的高頻系數(shù)');subplot(232);image(wcodemat(cvd2,nbc));title('尺度2垂直方向的高頻系數(shù)');subplot(233);image(wcodemat(cdd2,nbc));title('尺度2斜線方向的高頻系數(shù)');subplot(234);image(wcodemat(chd1,nbc));title('尺度1水平方向的高頻系數(shù)');subplot(235);image(wcodemat(cvd1,nbc));title('尺度1垂直方向的高頻系數(shù)');subplot(236);image(wcodemat(cdd1,nbc));title('尺度1斜線方向的高頻系數(shù)');程序運(yùn)行結(jié)果如圖11.15、圖11.16所示。圖11.15原始圖像以及尺度1和2的低頻系數(shù)圖11.16尺度1和2上的高頻系數(shù)

3.圖像的小波重構(gòu)

MATLAB中實(shí)現(xiàn)了單尺度二維離散小波逆變換、多尺度二維小波重構(gòu)、二維小波分解的低頻或高頻重構(gòu)、二維小波分解系數(shù)的直接重構(gòu)以及二維小波分解的單尺度重構(gòu)。

1)單尺度二維離散小波逆變換

MATLAB中實(shí)現(xiàn)單尺度二維離散小波逆變換的函數(shù)是idwt2，其調(diào)用格式有以下五種：

(1)?X=idwt2(cA,cH,cV,cD,'wname')；

(2)?X=idwt2(cA,cH,cV,cD,Lo_R,Hi_R)；

(3)?X=idwt2(cA,cH,cV,cD,wname',S)；

(4)?X=idwt2(cA,cH,cV,cD,Lo_R,Hi_R,S)；

(5)

X=idwt2(...,'mode',MODE)。其中輸入?yún)?shù)wname是小波函數(shù)，Lo_R和Hi_R分別為重構(gòu)低通濾波器和高通濾波器，兩者長(zhǎng)度必須相等。[cA,cH,cV,cD]是信號(hào)的單尺度小波分解結(jié)構(gòu)，返回值向量X為單尺度重構(gòu)后的信號(hào)。如果設(shè)sa=size(cA)=size(cH)=size(cV)=size(cD)，濾波器的長(zhǎng)度為lf，則size(X)=SX，此時(shí)若DWT拓展模式為周期，則SX=2*SA；若為其他模式，則SX=2*size(cA)-lf+2。格式(5)中的MODE的值參考函數(shù)dwtmode。

2)多尺度二維小波重構(gòu)

MATLAB中實(shí)現(xiàn)多尺度二維小波重構(gòu)的函數(shù)是waverec2，其調(diào)用格式有以下兩種：

(1)?X=waverec2(C,S,'wname')；

(2)?X=waverec2(C,S,Lo_R,Hi_R)。其中wname是小波函數(shù)，Lo_R和Hi_R分別為重構(gòu)低通濾波器和高通濾波器，兩者長(zhǎng)度必須相等。[C,S]是信號(hào)的小波分解結(jié)構(gòu)，向量X為重構(gòu)后的信號(hào)。

3)小波分解的低頻或高頻重構(gòu)

MATLAB中實(shí)現(xiàn)函數(shù)小波分解的低頻或高頻重構(gòu)的函數(shù)是wrcoef2，其調(diào)用格式有以下四種：

(1)?X=wrcoef2('type',C,S,'wname',N)；

(2)?X=wrcoef2('type',C,S,Lo_R,Hi_R,N)；

(3)?X=wrcoef2('type',C,S,'wname')；

(4)?X=wrcoef2('type',C,S,Lo_R,Hi_R)。其中[C,S]是二維信號(hào)的小波分解結(jié)構(gòu)，wname為所調(diào)用的小波函數(shù)，Lo_D和Hi_D分別為重構(gòu)低通濾波器和高通濾波器，type是指選擇進(jìn)行重構(gòu)信號(hào)的低頻或高頻部分，取值為?'a'或?'h'?(或?'v'、'd'?)，當(dāng)type取值為?'a'?時(shí)，表示對(duì)信號(hào)的低頻部分進(jìn)行重構(gòu)，此時(shí)N可以為0，當(dāng)type取值為'h'(或'v'、'd')時(shí)，表示對(duì)信號(hào)的高頻部分進(jìn)行重構(gòu)，此時(shí)N取值為正整數(shù)。當(dāng)type='a'時(shí)，0≤N≤size(S,1)-2；當(dāng)type='h'(或'v'、'd')時(shí)，1≤N≤size(S,1)-2。

4)小波分解系數(shù)的直接重構(gòu)

MATLAB中實(shí)現(xiàn)小波分解系數(shù)的直接重構(gòu)的函數(shù)是upcoef2，其調(diào)用格式有以下六種：

(1)?Y=upcoef2(O,X,'wname',N,S)；

(2)?Y=upcoef2(O,X,Lo_R,Hi_R,N,S)；

(3)?Y=upcoef2(O,X,'wname',N)；

(4)?Y=upcoef2(O,X,Lo_R,Hi_R,N)；

(5)?Y=upcoef2(O,X,'wname')；

(6)?Y=upcoef2(O,X,Lo_R,Hi_R)。其中wname為所調(diào)用的小波函數(shù)，Lo_D和Hi_D分別為重構(gòu)低通濾波器和高通濾波器。O取值為?'a'?或?'h'?(或?'v'、'd'?)，當(dāng)O取值為?'a'?時(shí)，表示對(duì)信號(hào)的低頻部分進(jìn)行重構(gòu)，當(dāng)O取值?'h'(或?'v'、'd')時(shí)，表示對(duì)信號(hào)的高頻部分進(jìn)行重構(gòu)。N為嚴(yán)格的正整數(shù)。格式(1)和(2)表示對(duì)信號(hào)X進(jìn)行重構(gòu)并返回中間長(zhǎng)度為S的部分。格式(5)和(6)中，N取值為1。

5)二維小波分解單尺度重構(gòu)

MATLAB中實(shí)現(xiàn)二維小波分解單尺度重構(gòu)的函數(shù)是upwlev2，其調(diào)用格式有以下兩種：

(1)?[NC,NS,cA]=upwlev2(C,S,'wname')；

(2)?[NC,NS,cA]=upwlev2(C,S,Lo_R,Hi_R)。其中[C,S]是信號(hào)的小波分解結(jié)構(gòu)，wname是小波函數(shù)，Lo_R和Hi_R分別為低通和高通重構(gòu)濾波器。函數(shù)返回分解結(jié)構(gòu)[C,S]的第N步重構(gòu)，即[NC,NS]是尺度N

-

1的分解結(jié)構(gòu)，cA是尺度N的低頻系數(shù)矩陣。

4.圖像的小波重構(gòu)實(shí)例下面通過(guò)幾個(gè)例子說(shuō)明如何進(jìn)行圖像的小波分解重構(gòu)。例11.3下列程序代碼說(shuō)明如何對(duì)圖像的單尺度分解進(jìn)行重構(gòu)。程序中首先調(diào)用函數(shù)dwt2對(duì)圖像進(jìn)行單尺度分解，然后調(diào)用函數(shù)idwt2對(duì)圖像進(jìn)行重構(gòu)。

loadwoman;

sX=size(X);

%使用小波函數(shù)db4進(jìn)行信號(hào)的單層分解

[cA1,cH1,cV1,cD1]=dwt2(X,'db4');

%進(jìn)行小波函數(shù)的重構(gòu)A0=idwt2(cA1,cH1,cV1,cD1,'db4',sX);%檢查重構(gòu)誤差ans=max(max(abs(X-A0)))nbc=size(map,1);colormap(pink(nbc));subplot(121);image(wcodemat(X,nbc));title('原始圖像');subplot(122);image(wcodemat(A0,nbc));title('重構(gòu)圖像');重構(gòu)誤差為error=1.4433e-010程序運(yùn)行結(jié)果如圖11.17所示。圖11.17單尺度二維小波重構(gòu)例11.4下列程序代碼說(shuō)明如何對(duì)圖像的多尺度小波分解重構(gòu)低頻或高頻信號(hào)。首先使用函數(shù)wavedec2對(duì)圖像進(jìn)行多尺度分解，然后利用函數(shù)wrcoef2重構(gòu)低頻和高頻信號(hào)。

loadwoman;

%對(duì)信號(hào)用小波函數(shù)sym5進(jìn)行2尺度分解

[c,s]=wavedec2(X,2,'sym5');

%對(duì)小波分解結(jié)構(gòu)[c,s]的低頻系數(shù)進(jìn)行尺度1和尺度2上的重構(gòu)

a1=wrcoef2('a',c,s,'sym5',1);

a2=wrcoef2('a',c,s,'sym5',2);

%對(duì)分解結(jié)構(gòu)[c,s]的高頻系數(shù)進(jìn)行?'h'、'v'、'd'?三個(gè)方向尺度2上的重構(gòu)

hd2=wrcoef2('h',c,s,'sym5',2);

vd2=wrcoef2('v',c,s,'sym5',2);

dd2=wrcoef2('d',c,s,'sym5',2);%所有的圖像都有相同的大小sX=size(X)sa1=size(a1)shd2=size(hd2)nbc=size(map,1);colormap(pink(nbc));subplot(3,2,1);image(wcodemat(X,nbc));title('原始圖像');subplot(3,2,2);image(wcodemat(a1,nbc));title('尺度1的低頻圖像');subplot(3,2,3);image(wcodemat(a2,nbc));title('尺度2的低頻圖像');subplot(3,2,4);image(wcodemat(hd2,nbc));title('尺度2的水平高頻圖像');subplot(3,2,5);image(wcodemat(vd2,nbc));title('尺度2的垂直高頻圖像');subplot(3,2,6);image(wcodemat(dd2,nbc));title('尺度2的斜線高頻圖像');程序運(yùn)行結(jié)果如圖11.18所示。圖11.18二維小波系數(shù)的單支重構(gòu)例11.5下列程序代碼說(shuō)明如何通過(guò)小波分解系數(shù)進(jìn)行小波重構(gòu)。首先使用函數(shù)dwt2對(duì)信號(hào)進(jìn)行小波分解，然后利用函數(shù)upcoef2對(duì)高頻和低頻系數(shù)進(jìn)行小波重構(gòu)。

loadwoman;

%使用小波函數(shù)db4對(duì)信號(hào)進(jìn)行2層分解

[c,s]=wavedec2(X,2,'db4');

siz=s(size(s,1),:);

%提取尺度1的低頻系數(shù)

ca1=appcoef2(c,s,'db4',1);%對(duì)尺度1的低頻部分進(jìn)行重構(gòu)a1=upcoef2('a',ca1,'db4',1,siz);%提取尺度1的水平方向高頻系數(shù)chd1=detcoef2('h',c,s,1);%對(duì)尺度1的水平方向高頻部分進(jìn)行重構(gòu)hd1=upcoef2('h',chd1,'db4',1,siz);%提取尺度1的垂直方向高頻系數(shù)cvd1=detcoef2('v',c,s,1);%對(duì)尺度1的垂直方向高頻部分進(jìn)行重構(gòu)vd1=upcoef2('v',cvd1,'db4',1,siz);%提取尺度1的斜線方向高頻系數(shù)cdd1=detcoef2('d',c,s,1);%對(duì)尺度1的斜線方向高頻部分進(jìn)行重構(gòu)dd1=upcoef2('d',cdd1,'db4',1,siz);nbc=size(map,1);colormap(pink(nbc));subplot(321);image(wcodemat(X,nbc));title('原始圖像');subplot(323);image(wcodemat(a1,nbc));title('尺度1的低頻系數(shù)重構(gòu)圖像');subplot(324);image(wcodemat(hd1,nbc));title('尺度1的水平高頻系數(shù)重構(gòu)圖像');subplot(325);image(wcodemat(vd1,nbc));title('尺度1的垂直高頻系數(shù)重構(gòu)圖像');subplot(326);image(wcodemat(dd1,nbc));title('尺度1的斜線高頻系數(shù)重構(gòu)圖像');程序運(yùn)行如圖11.19所示。

圖11.19二維小波分解系數(shù)的直接重構(gòu)例11.6下列程序代碼實(shí)現(xiàn)對(duì)小波多尺度分解進(jìn)行重構(gòu)，即重構(gòu)原圖像(包括低頻和高頻)。首先使用函數(shù)wavedec2對(duì)圖像進(jìn)行多尺度分解，然后調(diào)用函數(shù)waverec進(jìn)行圖像的重構(gòu)。

loadwoman;

%使用小波函數(shù)sym4對(duì)信號(hào)進(jìn)行2尺度分解

[c,s]=wavedec2(X,2,'sym4');

%進(jìn)行小波函數(shù)的重構(gòu)a0=waverec2(c,s,'sym4');%檢查重構(gòu)誤差ans=max(max(abs(X-a0)))nbc=size(map,1);colormap(pink(nbc));subplot(121);image(wcodemat(X,nbc));title('原始圖像');subplot(122);image(wcodemat(a0,nbc));title('重構(gòu)圖像');程序運(yùn)行結(jié)果如圖11.20所示。圖11.20多尺度二維小波重構(gòu)11.3.2小波圖像壓縮

MATLAB中的小波圖像壓縮主要包括獲取壓縮閾值和進(jìn)行圖像壓縮兩方面。實(shí)現(xiàn)獲取壓縮閾值的函數(shù)有ddencmp和wdcbm2，實(shí)現(xiàn)圖像壓縮的函數(shù)有wdencmp、wpdencmp和wthcoef2。

1.小波圖像壓縮函數(shù)

MATLAB中實(shí)現(xiàn)信號(hào)閾值獲取的函數(shù)有ddencmp和wdcbm2，下面對(duì)它們的用法進(jìn)行簡(jiǎn)單的說(shuō)明?！?ddencmp的調(diào)用格式有以下三種：

(1)?[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,IN2,X)；

(2)?[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,'wp',X)；

(3)?[THR,SORH,KEEPAPP]=ddencmp(IN1,'wv',X)。函數(shù)ddencmp用于獲取在消噪或壓縮過(guò)程中的默認(rèn)閾值。輸入?yún)?shù)X為一維或二維信號(hào)；IN1取值為?'den'?或?'cmp'?，'den'?表示去噪，'cmp'?表示壓縮；IN2取值為?'wv'或?'wp'，'wv'?表示選擇小波，'wp'?表示選擇小波包。返回值THR是返回的閾值；SORH是軟閾值或硬閾值選擇參數(shù)；KEEPAPP表示保存低頻信號(hào)；CRIT是熵名(只在選擇小波包時(shí)用)?！窈瘮?shù)wdcbm2的調(diào)用格式有以下兩種：

(1)?[THR,NKEEP]=wdcbm2(C,S,ALPHA)；

(2)?[THR,NKEEP]=wdcbm2(C,S,ALPHA,M)。函數(shù)wdcbm2用于使用Birgé-Massart算法獲取二維小波變換的閾值。返回值THR是與尺度無(wú)關(guān)的閾值，NKEEP是系數(shù)的個(gè)數(shù)。[C，S]是要進(jìn)行消噪或壓縮的信號(hào)的在j?=length(L)-2層的分解結(jié)構(gòu)；ALPHA和M必須是大于1的實(shí)數(shù)；THR是關(guān)于j的向量，THR(i)是第i層在水平、垂直及斜線方向上的閾值；NKEEP也是關(guān)于j的向量，NKEEP(i)是i層的系數(shù)個(gè)數(shù)。一般壓縮時(shí)ALPHA取1.5，去噪時(shí)ALPHA取3。

MATLAB中實(shí)現(xiàn)信號(hào)的閾值壓縮的函數(shù)有wdencmp、wthcoef2以及wpdencmp。下面對(duì)它們的用法進(jìn)行簡(jiǎn)單的介紹。●函數(shù)wdencmp的調(diào)用格式有下面三種：

(1)?[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('gbl',X,'wname',N,THR,SORH,KEEPAPP)；

(2)?[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',X,wname',N,THR,SORH)；

(3)?[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',C,L,'wname',N,THR,SORH)。函數(shù)wdencmp用于一維或二維信號(hào)的消噪或壓縮。'gbl'(global的縮寫(xiě))表示每層都采用同一個(gè)閾值進(jìn)行處理，'wname'是所用的小波函數(shù)，'lvd'表示每層用不同的閾值進(jìn)行處理，N表示小波分解的層數(shù)，THR為閾值向量，對(duì)于格式(2)和(3)每層都要求有一個(gè)閾值，因此閾值向量THR的長(zhǎng)度為N，SORH表示選擇軟閾值或硬閾值(分別取值為's'和'h')，參數(shù)KEEPAPP取值1時(shí)，則低頻系數(shù)不進(jìn)行閾值量化，反之，則低頻系數(shù)要進(jìn)行閾值量化。

XC是消噪或壓縮后的信號(hào)，[CXC,LXC]是XC的小波分解結(jié)構(gòu)，PERF0和PERFL2是恢復(fù)和壓縮L2范數(shù)百分比。PERFL2?=?100?×?(XD的小波包系數(shù)范數(shù)/X的小波包系數(shù))2。如果X是一維信號(hào)，小波?'wname'?是一個(gè)正交小波，則?！窈瘮?shù)wthcoef2的調(diào)用格式有下面四種：(1)?NC=wthcoef2('type',C,S,N,T,SORH)；(2)?NC=wthcoef2('type',C,S,N)；(3)?NC=wthcoef2('a',C,S)；(4)?NC=wthcoef2('t',C,S,N,T,SORH)。函數(shù)wthcoef2用于二維信號(hào)小波系數(shù)的閾值處理。參數(shù)?'type'?的取值為?'h(?'v'?或?'d'?)時(shí)，NC=wthcoef2(?'type',C,S,N,T,SORH)返回小波分解結(jié)構(gòu)[C,S]經(jīng)由向量N和T定義的閾值量化后的水平(垂直或斜線)系數(shù)。N是包含要進(jìn)行閾值量化的細(xì)節(jié)尺度向量，T是相應(yīng)的閾值。N和T的長(zhǎng)度必須相等，向量N的長(zhǎng)度必須滿足1≤N(i)≤size(S,1)-2。其中SORH='s'表示軟閾值處理，SORH='h'表示硬閾值處理。函數(shù)wpdencmp的調(diào)用格式有下面兩種：

(1)?[XD,TREED,PERF0,PERFL2]=wpdencmp(X,SORH,N,'wname',CRIT,PAR,KEEPAPP)；

(2)?[XD,TREED,PERF0,PERFL2]=wpdencmp(TREE,SORH,CRIT,PAR,KEEPAPP)。函數(shù)wpdencmp用于使用小波包變換進(jìn)行信號(hào)的壓縮或去噪。格式(1)返回輸入信號(hào)X(一維或二維)的去噪或壓縮后的信號(hào)XD。輸出參數(shù)TREED是XD的最佳小波包分解樹(shù)；PERFL0和PERFL2是恢復(fù)和壓縮L2能量百分比。PERFL2=100*(XD的小波包系數(shù)范數(shù)/X的小波包系數(shù))2。如果X是一維信號(hào)，'wname'?是正交小波，則。SORH的取值為?'s'?或?'h'，表示軟閾值或硬閾值。輸入?yún)?shù)N是小波包分解的層數(shù)，'wname'?是小波名。函數(shù)使用由字符串CRIT定義的熵標(biāo)準(zhǔn)和閾值參數(shù)PAR實(shí)現(xiàn)最佳分解。如果KEEPAPP=1，則近似信號(hào)的小波系數(shù)不進(jìn)行閾值量化；否則，進(jìn)行量化。格式(2)與格式(1)的輸出參數(shù)相同，輸入選項(xiàng)也相同，只是它從信號(hào)的小波包分解樹(shù)TREE直接進(jìn)行去噪或壓縮。

2.小波圖像壓縮實(shí)例下面將通過(guò)例11.7～11.9說(shuō)明如何進(jìn)行圖像壓縮。例11.7下列程序代碼使用函數(shù)wdcbm2獲取圖像壓縮閾值，然后采用函數(shù)wdencmp實(shí)現(xiàn)圖像壓縮：

loadjulia;

nbc=size(map,1);

%使用haar小波進(jìn)行三層分解wname='db2';lev=3;[c,s]=wavedec2(X,lev,wname);%使用wdcbm2獲得壓縮閾值alpha=1.5;m=3.5*prod(s(1,:));[thr,nkeep]=wdcbm2(c,s,alpha,m);%對(duì)信號(hào)進(jìn)行壓縮[xc,cxc,lxc,PERF0,PERFL2]=wdencmp('lvd',c,s,wname,lev,thr,'h');%計(jì)算壓縮比nzeros1=find(c);nzeros2=find(cxc);cmp_rate=size(nzeros2)/size(nzeros1)colormap(pink(nbc));figure(1);subplot(1,2,1);image(wcodemat(X,nbc));title('原始信號(hào)');subplot(1,2,2)image(wcodemat(XC,nbc))title('壓縮后的信號(hào)');程序運(yùn)行結(jié)果如圖11.21所示。圖11.21圖像壓縮統(tǒng)計(jì)小波系數(shù)中非零元素值并計(jì)算圖像壓縮比：cmp_rate=0.0950例11.8在本例中使用函數(shù)ddencmp獲取圖像壓縮閾值，然后采用函數(shù)wdencmp實(shí)現(xiàn)圖像壓縮。程序代碼如下：

loadjulia;

nbc=size(map,1);

%使用db2小波進(jìn)行三層分解

wname='db2';lev=3;[c,s]=wavedec2(X,lev,wname);%使用ddencmp獲得壓縮閾值[thr,nkeep]=ddencmp('cmp','wv',X)%對(duì)信號(hào)進(jìn)行壓縮[xc,cxc,lxc,PERF0,PERFL2]=wdencmp('gbl',c,s,wname,lev,thr,'s',1);%計(jì)算壓縮比zeros1=find(c==0);zeros2=find(cxc==0);zeros_rate=size(nzeros2)/size(nzeros1)colormap(pink(nbc));figure(1);subplot(1,2,1);image(wcodemat(X,nbc));title('原始信號(hào)');subplot(1,2,2)image(wcodemat(xc,nbc))title('壓縮后的信號(hào)');程序運(yùn)行結(jié)果如圖11.22所示。圖11.22圖像壓縮統(tǒng)計(jì)壓縮計(jì)算前后，小波系數(shù)中零元素的個(gè)數(shù)比例為6.28。可見(jiàn)壓縮處理后，小波系數(shù)中零元素的個(gè)數(shù)顯著增加，這時(shí)可以利用第8章中介紹的行程編碼對(duì)其進(jìn)行編碼壓縮，從而取得較好的壓縮效果。例11.9在本例中首先使用函數(shù)ddencmp獲得圖像壓縮閾值，然后才用函數(shù)wpdencmp對(duì)圖像進(jìn)行小波包壓縮。程序代碼如下：%加載原始信號(hào)loadmask;nbc=size(map,1);%使用wpdencmp進(jìn)行小波壓縮%使用ddencmp獲取閾值[thr,sorh,keepapp,crit]=ddencmp('cmp','wp',X);%使用閾值進(jìn)行小波包壓縮[xc,treed,perf0,perfl2]=wpdencmp(X,sorh,3,'db3',crit,thr,keepapp);subplot(211);colormap(pink(nbc));image(wcodemat(X,nbc))title('原始圖像');subplot(212);image(wcodemat(xc,nbc))title('壓縮后的圖像');程序運(yùn)行結(jié)果如圖11.23所示。圖11.23圖像小波包壓縮11.3.3小波圖像去噪

1.圖像閾值去噪概述圖像去噪與圖像壓縮有一定的相似性。不同的是圖像去噪需要用二維的分析工具對(duì)圖像進(jìn)行分析，在對(duì)高頻系數(shù)分別進(jìn)行閾值量化處理時(shí)，要對(duì)圖像分解的三個(gè)方向的高頻系數(shù)進(jìn)行閾值處理。圖像去噪的一般方法如下：

(1)圖像的小波分解。選擇合適的小波函數(shù)以及適合的分解層次對(duì)圖像進(jìn)行分解。

(2)對(duì)分解后的高頻系數(shù)進(jìn)行閾值處理。對(duì)分解的每一層，選擇合適的閾值對(duì)該層的水平、垂直和斜線三個(gè)方向的高頻系數(shù)進(jìn)行閾值量化處理。閾值函數(shù)和閾值選擇規(guī)則同信號(hào)的閾值消噪相同。

(3)重構(gòu)圖像。根據(jù)小波分解的低頻系數(shù)和經(jīng)閾值量化處理后的高頻系數(shù)進(jìn)行圖像重構(gòu)。

MATLAB中實(shí)現(xiàn)了圖像的閾值去噪，主要是閾值獲取和閾值去噪兩方面。

MATLAB中實(shí)現(xiàn)閾值獲取的函數(shù)有ddencmp、thselect、wbmpen和wdcbm2。這些函數(shù)在前面已經(jīng)介紹過(guò)。

MATLAB中實(shí)現(xiàn)閾值去噪的函數(shù)主要有wden、wdencmp、wpdencmp、wthresh、wpthcoef和wthcoef2。其中函數(shù)wden、wdencmp、wpdencmp、wthcoef2和wpthcoef的用法前面已介紹，下面只簡(jiǎn)單介紹一下函數(shù)wthresh的用法。函數(shù)wthresh的調(diào)用格式如下：

Y=wthresh(X,SORH,T)

Y=wthresh(X,SORH,T)返回輸入向量或矩陣X經(jīng)軟閾值(SORH='s')或硬閾值(SORH='h')

T處理后的信號(hào)。

Y=wthresh(X,'s',T)返回的是Y?=?sign(X)·(|X|?-?T)+，即把信號(hào)的絕對(duì)值與閾值進(jìn)行比較，小于或等于閾值的點(diǎn)變?yōu)?，大于閾值的點(diǎn)變?yōu)樵擖c(diǎn)值與閾值的差值。

Y=wthresh(X,'h',T)返回的是Y?=?X·1(|X|?＞?T)，即把信號(hào)的絕對(duì)值與閾值比較，小于或等于閾值的點(diǎn)變?yōu)?，大于閾值的點(diǎn)保持不變。一般來(lái)說(shuō)，用硬閾值處理后的信號(hào)比軟閾值處理后的信號(hào)更粗糙。

2.圖像閾值去噪實(shí)例本節(jié)通過(guò)幾個(gè)例子說(shuō)明如何應(yīng)用MATLAB的函數(shù)進(jìn)行圖像的閾值去噪。例11.10使用函數(shù)wdcbm2獲取圖像去噪閾值，然后采用函數(shù)wdencmp實(shí)現(xiàn)圖像去噪。程序代碼如下：

loadwoman;

%產(chǎn)生含噪信號(hào)

init=2055615866;

randn('seed',init);x=X+18*randn(size(X));nbc=size(map,1);%使用sym4小波進(jìn)行三層分解wname='sym4';lev=3;[c,s]=wavedec2(x,lev,wname);%使用wdcbm2獲得去噪閾值alpha=3;m=3.5*prod(s(1,:));[thr,nkeep]=wdcbm2(c,s,alpha,m);%對(duì)信號(hào)進(jìn)行去噪sorh='s';xd=wdencmp('lvd',c,s,wname,lev,thr,sorh);colormap(pink(nbc));figure(1);subplot(2,2,1);image(wcodemat(X,nbc));title('原始信號(hào)');subplot(2,2,2);image(wcodemat(x,nbc));title('含噪信號(hào)');subplot(2,2,3);image(wcodemat(xd,nbc));title('消噪后的信號(hào)');程序運(yùn)行結(jié)果如圖11.24所示。

圖11.24圖像小波閾值去噪例11.11在本例中首先使用函數(shù)wpcoef獲取圖像的小波包分解系數(shù)，接著估計(jì)噪聲方差，然后使用函數(shù)wpbmpen獲取小波包去噪閾值，最后使用函數(shù)wpdencmp實(shí)現(xiàn)圖像去噪。程序代碼如下：

loadwoman;

%產(chǎn)生含噪信號(hào)

init=2055615866;

randn('seed',init);

x=X+18*randn(size(X));%對(duì)圖像進(jìn)行小波包分解%使用coif2分解到第三層wname='coif2';lev=3;tree=wpdec2(x,lev,wname);%估計(jì)噪聲標(biāo)準(zhǔn)差det1=[wpcoef(tree,2)wpcoef(tree,3)wpcoef(tree,4)];sigma=median(abs(det1(:)))/0.6745;%選擇去噪的閾值alpha=1.1;thr=wpbmpen(tree,sigma,alpha);%使用上面的閾值對(duì)圖像去噪keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp);%畫(huà)出原始圖像和去噪后的圖像nbc=size(map,1);colormap(pink(nbc));subplot(2,2,1);image(wcodemat(X,nbc));title('原始信號(hào)');subplot(2,2,2);image(wcodemat(x,nbc));title('含噪信號(hào)');subplot(2,2,3);image(wcodemat(xd,nbc));title('消噪后的信號(hào)');