《MATLAB应用图像处理》课件第11章

第11章MATLAB小波分析在图像处理中的应用

11.1小波和小波变换

11.2离散小波变换和Mallat算法

11.3MATLAB小波分析在图像处理中的应用习题

11.1小波和小波变换

小波(Wavelet)这一术语，顾名思义，“小波”就是小区域、长度有限、均值为0的波形。所谓“小”是指它具有衰减性，而称之为“波”则是指它的波动性，其振幅正负相间的振荡形式。小波是定义在有限间隔且其平均值为零的一种函数，它的波形如图11.1所示。图11.1是从许多使用比较广泛的小波中挑选出的几种一维小波。图11.1所示的小波中，缩放函数和小波函数的名称大多数是以开发者的名字命名的，例如Moret小波函数是Grossmann和Morlet在1984年开发的，db6缩放函数和db6小波函数是Daubechies开发的几种小波之一，Meyer缩放函数和Meyer小波函数是Meyer开发的。但也有不少例外，例如sym6缩放函数和sym6小波函数则是symlets的简写，是Daubechies提议开发的几种对称小波之一，coif2缩放函数和coif2小波函数是Daubechies应R.Coifman的请求而开发的几种小波之一。图11.1小波波形图11.1所示的小波具有有限的持续时间和突变的频率和振幅，波形可以是不规则的，也可以是不对称的，在整个时间范围里的幅度平均值为零。与傅立叶分析相比，小波分析能更有效地处理不规则的激烈变化信号。傅立叶变换是把一个信号分解成各种不同频率的正弦波，因此正弦波是傅立叶变换的基函数。同样，小波分析是把一个信号分解成将原始小波经过移位和缩放之后的一系列小波，因此小波同样可以用来表示一些函数的基函数。设x(t)是平方可积函数，即x(t)∈L2(R)，则x(t)的小波变换(WaveletTransform，WT)定义为其中，ψ(t)是基本小波或母小波(MotherWavelet)函数，*为的共轭，而是基本小波ψ(t)的位移和尺度伸缩，也称为ψ(t)的生成小波。式中a＞0，称为尺度因子，b反映时间位移，其值可正可负，符号<·>表示内积。(11.1)式(11.1)中的t、a和b均为连续变量，因此式(11.1)称为连续的小波变换(ContinuousWaveletTransform，CWT)。尺度因子a的作用是对基本小波ψ(t)作伸缩，a愈大ψ(t/a)愈宽。对于一个持续时间有限的小波，ψ(t)与ψab(t)之间的关系以及不同尺度a下小波分析区间的变化可用图11.2表示。从图11.2中可以看出，小波的持续时间随a的增大而加宽，幅度则与成反比，但波形保持不变。图11.2小波的移位和伸缩

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

(11.2)它的频谱Ψ(ω)为(11.3)由于(ω)是中心频率在ω0处的高斯型函数，如图11.3(a)所示，形状近似于矩形函数，因此可以表征X(ω)在ω0附近的局部性质，X(ω)为x(t)的频谱。如果采用不同的尺度伸缩因子a，(aω)的中心频率和带宽将发生变化。例如当a = 2时，ψ(t/2)的傅立叶变换为(11.4)可见此时中心频率降到ω0/2，而带宽的比例系数由2T-1/2变为T-1/2，如图11.3(b)所示，因而|(aω)|的品质因数不变。图11.3小波中心频率和带宽随尺度的变化总之，从频域上看，用不同的尺度作小波变换，相当于用一组中心频率不同的带通滤波器对信号进行处理。带通滤波器既可以对信号进行分解，也可以用于信号的检测。图11.4表示小波变换在时频平面上的基本分析单元的特点。当a值小时，ψ(t/a)很“窄”，因此在时间轴上的观测范围小，可以“细致观察”时域波形的变化；而在频域上相当于用较高频率的小波，对信号的频谱作分辨率较低的分析。当a值较大时，ψ(t/a)变“宽”，时间轴上的观测范围大，可以“初略观察”时域波形；而在频域上相当于用低频小波对信号作分辨率较高的分析。分析频率有高有低，但在各分析频段内的品质因数Q却保持恒定。图11.4小波变换的时频特性这种分析特点是工程实际所期望的。对于高频信号，我们希望在时域有较高的分辨率，而在频域上的分辨率则允许相应地降低，因此就要求用窄的小波(即用小的尺度因子a)来“仔细观测”时域波形x(t)。反之，对于低频信号，我们希望提高频域上的分辨率，而对时域的分辨率则可以降低要求，此时，就应当用具有窄带频谱的小波(即用较大的尺度因子a)来观察X(ω)。如上所述，小波分析恰恰具有自动调整时域和频域的“视野”和分析频率高低的功能，从而保证了各分析频段内的品质因数Q的不变性。11.2离散小波变换和Mallat算法

在用计算机进行小波分析时，实际上是用离散的数据进行计算的，希望提高计算效率。不难想象，连续小波变换的计算量是惊人的。为了解决计算量的问题，缩放因子和平移参数都选择2的整数倍。使用这样的缩放因子和平移参数的小波变换叫做双尺度小波变换，它是离散小波变换(DiscreteWaveletTransform，DWT)的一种形式。总体来说，离散小波变换通常指的是双尺度小波变换。执行离散小波变换的有效方法是使用滤波器。该方法是Mallat在1988年开发的，叫做Mallat算法，这种方法实际上是一种信号的分解方法，在数字信号处理中称为双通道子带编码。双通道子带编码主要涉及滤波器族理论、尺度分析和子带编码算法。子带编码算法使用滤波器实现小波变换的算法，如图11.5所示。图11.5双通道编码滤波算法图11.5中，原始的输入信号通过两个互补的滤波器产生两个信号：大尺度低频分量和小尺度高频分量。在小波分析中，低频分量是大尺度的缩放因子产生的系数，而高频分量是小的缩放因子产生的系数。低频分量表示信号的近似值，高频分量表示信号的细节。图11.5只表示了在一个分辨率下的小波分解。当需要进行多分辨率分解时，离散小波变换可以被表示成由低通滤波器和高通滤波器组成的一棵树。信号的分解过程可以迭代，也就是说可进行多级分解。如果对信号的高频分量不再分解，而对低频分量连续进行分解，就得到许多分辨率较低的低频分量，形成如图11.6所示的一棵树，这种树叫做小波分解树。分解级数的多少取决于待分析的数据和用户的需要。图11.6多级小波分解树小波分解树表示只对信号的低频分量进行连续分解。如果不仅对信号的低频分量进行连续分解，而且对高频分量也进行连续分解，这样不仅可得到许多分辨率较低的低频分量，而且也可得到许多分辨率较低的高频分量。这样分解得到的树叫做小波包分解树(waveletpacketdecompositiontree)，这种树是一个完整的二进制树。图11.7表示的是一棵三级小波包分解树。小波包分解方法是小波分解的一般化，可为信号分析提供更丰富和更详细的信息。图11.7三级小波包分解树在使用滤波器对数字信号进行变换时，得到的数据量将是原始数据量的两倍。例如，如果原始信号的数据样本为1000个，通过滤波之后每一个通道的数据均为1000个，总共为2000个。于是，根据奈奎斯特(Nyquist)采样定理可以采用降采样(downsampling)的方法来解决问题，即在每个通道中每两个样本数据中取一个，如图11.8所示。图11.8降采样过程

MATLAB中用wavedec函数实现多级离散小波分解，wavedec的调用格式如下：

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

[C,L]=wavedec(X,N,Lo_D,Hi_D)其中X为输入信号，N为分解级数，wname为小波名称，Lo_D为低通滤波器，Lo_H为高通滤波器。返回的参数C为各级小波参数，L为各级小波参数的数据长度。图11.9为wavedec小波分解示意图。图11.9wavedec小波分解下面一段代码使用db2小波，即Daubechies小波，对一个正弦叠加信号进行3级小波分解。图11.10描述了信号的fft分解和小波分解的结果。

%调入正弦叠加信号

loadsumsin;s=sumsin;

%使用db2小波对信号进行3级分解

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

subplot(3,1,1);

plot(s);title('输入信号');subplot(3,1,2);plot(abs(fft(s)))；title('输入信号fft变换');subplot(3,1,3);plot(C);title('小波系数');axis([01000-55]);

L中的各级小波系数个数：

L=1271272525011000信号的小波重构是小波分解的逆过程，MATLAB中用idwt实现信号的小波重构。idwt的调用格式如下：

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级小波分解

idwt使用的算法如图11.11所示，其中2倍上采样是通过数据间插0实现的，低通和高通滤波器输出的数值相加后取中间的L个数值得到j - 1阶小波系数，L默认为cAj(或cDj)的长度，也可以自由设定。图11.11小波重构算法

idwt实现一级小波系数的重构，多次调用idwt可实现多级小波系数的重构。下面一段代码通过3次调用idwt并将1级细节小波系数清零，实现信号重构，计算结果如图11.12所示。由于清除了1级细节小波系数，图11.12中的重构信号比输入信号更为光滑。

%调入正弦叠加信号

loadsumsin;s=sumsin;

%使用db2小波对信号进行3级分解

[C,L]=wavedec(s,3,'db2');%构建3级小波系数cA3=C(1:1:L(1));cD3=C(L(1)+1:1:L(1)+L(2));%构建2级小波系数cA2=idwt(cA3,cD3,'db2');cD2=C(L(1)+L(2)+1:1:L(1)+L(2)+L(3));%构建1级小波系数cA1=idwt(cA2,cD2,'db2',L(5));%将1级细节小波系数清零cD1=zeros(1,L(4));%重构信号ss=idwt(cA1,cD1,'db2',L(5));subplot(2,1,1);plot(s);title('输入信号');subplot(2,1,2);plot(ss)title('重构信号');axis([01000-55]);图11.12小波信号重构11.3MATLAB小波分析在图像处理

中的应用

11.3.1小波图像分解与重构

1.二维小波分解

MATLAB中实现了单尺度二维离散小波分解、多尺度二维离散小波分解、提取二维小波分解的高频系数以及提取二维小波分解的低频系数。

1)单尺度二维离散小波分解

MATLAB中实现单尺度二维离散小波分解的函数是dwt2，其调用格式有以下两种：

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

(2)[cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D)。其中X是被分解的离散信号，wname为分解小波函数，Lo_D和Hi_D分别为分解低通滤波器和高通滤波器，两者长度必须相等。返回值cA、cH、cV和cD分别为低频系数和高频系数向量。可以按照格式[cA,cH,cV,cD]=dwt2(...,'mode',MODE)，使用扩展模式进行信号的小波分解，参数MODE可以参考函数dwtmode。如果sx=size(X)，lf为滤波器的长度，则size(cA)=size(cH)=size(cV)=size(cD)=sa，如果DWT的扩展模式为周期，此时sa=ceil(sx/2)，对于其它的模式，此时sa=floor((sx+lf-1)/2)。

2)多尺度二维离散小波分解

MATLAB中实现多尺度二维离散小波分解的函数是wavedec2，其调用格式有以下两种：

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

(2) [C,S]=wavedec2(X,N,Lo_D,Hi_D)。其中N必须为严格的正整数，wname为小波函数，Lo_D和Hi_D分别为分解低通滤波器和高通滤波器。返回分解结构[C,S]，其中C的结构为

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)为低频系数，H(N)为水平高频系数，V(N)为垂直高频系数，D(N)为斜线高频系数，所有向量均以列向量存储在矩阵C中。矩阵S的结构如下：

S(1,:)=尺度N的低频系数矩阵的大小；

S(i,:)=尺度N-i+2的高频系数矩阵的大小，i=2,…，N+1；

S(N+2,:)=size(X)。其分解过程如图11.13所示。图11.13信号的多尺度二维小波分解

3)提取二维小波分解的高频系数

MATLAB中实现提取二维小波分解的高频系数的函数是detcoef2，其调用格式如下：

D=detcoef2(O,C,S,N)其中[C,S]是信号的二维小波分解结构，尺度N必须为正整数且1≤N≤size(S,1) - 2，O取值为'h'('v'或'd')，当O为 'h' 时，表示提取水平方向系数，当O为 'v' 时，表示提取垂直方向系数，当O为 'd' 时，表示提取斜线方向系数。

4)提取二维小波分解的低频系数MATLAB中实现提取二维小波分解的低频系数的函数是appcoef2，其调用格式有以下四种：

(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对应，用于提取二维小波分解的低频系数。wname为小波函数，Lo_R和Hi_R分别为低通滤波器和高通滤波器。格式(1)和(4)提取尺度N(0≤N≤size(S,1)-2)的低频系数，格式(2)和(3)提取最后一尺度(N=size(S,1)-2)的低频系数。

2.图像的小波分解实例下面通过例11.1和11.2说明如何对图像进行单尺度分解和多尺度分解并提取多尺度分解的小波系数。例11.1在下列程序代码中，调用函数dwt2对图像进行单尺度分解，并画出图像分解的低频分量和水平、垂直及斜线方向的三个高频分量。

%加载图像

loadtire;

nbc=size(map,1);colormap(pink(nbc));%用小波函数db2对信号x进行单层分解[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('低频系数图像');subplot(3,2,4);image(wcodemat(cH1,nbc));title('水平高频系数图像');subplot(3,2,5);image(wcodemat(cV1,nbc));title('垂直高频系数图像');subplot(3,2,6);image(wcodemat(cD1,nbc));title('斜线高频系数图像');程序运行结果如图11.14所示。可以看出低频分量表现了图像的轮廓，而高频分量表现了图像的细节。图11.14单尺度图像分解例11.2下列程序代码对图像进行多尺度分解并提取图像分解的低频和高频小波系数。程序中首先利用函数wavedec2对图像进行2层分解，然后利用函数detcoef2提取小波分解的高频系数，再利用函数appcoef2提取小波分解的低频系数。

loadtire;

%用小波函数db1对信号进行2尺度分解

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

sizex=size(X)；

sizec=size(c)；val_s=s%提取尺度2的所有方向的高频系数[chd2,cvd2,cdd2]=detcoef2('all',c,s,2);sizecd2=size(chd2)；%提取尺度1的所有方向的高频系数[chd1,cvd1,cdd1]=detcoef2('all',c,s,1);sizecd1=size(chd1)；%提取尺度2的低频系数ca2=appcoef2(c,s,'db1',2);sizeca2=size(ca2)；%提取尺度1的低频系数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的低频系数');subplot(224);image(wcodemat(ca1,nbc));title('尺度1的低频系数');figure(2);subplot(231);image(wcodemat(chd2,nbc));title('尺度2水平方向的高频系数');subplot(232);image(wcodemat(cvd2,nbc));title('尺度2垂直方向的高频系数');subplot(233);image(wcodemat(cdd2,nbc));title('尺度2斜线方向的高频系数');subplot(234);image(wcodemat(chd1,nbc));title('尺度1水平方向的高频系数');subplot(235);image(wcodemat(cvd1,nbc));title('尺度1垂直方向的高频系数');subplot(236);image(wcodemat(cdd1,nbc));title('尺度1斜线方向的高频系数');程序运行结果如图11.15、图11.16所示。图11.15原始图像以及尺度1和2的低频系数图11.16尺度1和2上的高频系数

3.图像的小波重构

MATLAB中实现了单尺度二维离散小波逆变换、多尺度二维小波重构、二维小波分解的低频或高频重构、二维小波分解系数的直接重构以及二维小波分解的单尺度重构。

1)单尺度二维离散小波逆变换

MATLAB中实现单尺度二维离散小波逆变换的函数是idwt2，其调用格式有以下五种：

(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)。其中输入参数wname是小波函数，Lo_R和Hi_R分别为重构低通滤波器和高通滤波器，两者长度必须相等。[cA,cH,cV,cD]是信号的单尺度小波分解结构，返回值向量X为单尺度重构后的信号。如果设sa=size(cA)=size(cH)=size(cV)=size(cD)，滤波器的长度为lf，则size(X)=SX，此时若DWT拓展模式为周期，则SX=2*SA；若为其他模式，则SX=2*size(cA)-lf+2。格式(5)中的MODE的值参考函数dwtmode。

2)多尺度二维小波重构

MATLAB中实现多尺度二维小波重构的函数是waverec2，其调用格式有以下两种：

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

(2) X=waverec2(C,S,Lo_R,Hi_R)。其中wname是小波函数，Lo_R和Hi_R分别为重构低通滤波器和高通滤波器，两者长度必须相等。[C,S]是信号的小波分解结构，向量X为重构后的信号。

3)小波分解的低频或高频重构

MATLAB中实现函数小波分解的低频或高频重构的函数是wrcoef2，其调用格式有以下四种：

(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]是二维信号的小波分解结构，wname为所调用的小波函数，Lo_D和Hi_D分别为重构低通滤波器和高通滤波器，type是指选择进行重构信号的低频或高频部分，取值为 'a'或 'h' (或 'v'、'd' )，当type取值为 'a' 时，表示对信号的低频部分进行重构，此时N可以为0，当type取值为'h'(或'v'、'd')时，表示对信号的高频部分进行重构，此时N取值为正整数。当type='a'时，0≤N≤size(S,1)-2；当type='h'(或'v'、'd')时，1≤N≤size(S,1)-2。

4)小波分解系数的直接重构

MATLAB中实现小波分解系数的直接重构的函数是upcoef2，其调用格式有以下六种：

(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为所调用的小波函数，Lo_D和Hi_D分别为重构低通滤波器和高通滤波器。O取值为 'a' 或 'h' (或 'v'、'd' )，当O取值为 'a' 时，表示对信号的低频部分进行重构，当O取值 'h'(或 'v'、'd')时，表示对信号的高频部分进行重构。N为严格的正整数。格式(1)和(2)表示对信号X进行重构并返回中间长度为S的部分。格式(5)和(6)中，N取值为1。

5)二维小波分解单尺度重构

MATLAB中实现二维小波分解单尺度重构的函数是upwlev2，其调用格式有以下两种：

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

(2) [NC,NS,cA]=upwlev2(C,S,Lo_R,Hi_R)。其中[C,S]是信号的小波分解结构，wname是小波函数，Lo_R和Hi_R分别为低通和高通重构滤波器。函数返回分解结构[C,S]的第N步重构，即[NC,NS]是尺度N

-

1的分解结构，cA是尺度N的低频系数矩阵。

4.图像的小波重构实例下面通过几个例子说明如何进行图像的小波分解重构。例11.3下列程序代码说明如何对图像的单尺度分解进行重构。程序中首先调用函数dwt2对图像进行单尺度分解，然后调用函数idwt2对图像进行重构。

loadwoman;

sX=size(X);

%使用小波函数db4进行信号的单层分解

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

%进行小波函数的重构A0=idwt2(cA1,cH1,cV1,cD1,'db4',sX);%检查重构误差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('重构图像');重构误差为error=1.4433e-010程序运行结果如图11.17所示。图11.17单尺度二维小波重构例11.4下列程序代码说明如何对图像的多尺度小波分解重构低频或高频信号。首先使用函数wavedec2对图像进行多尺度分解，然后利用函数wrcoef2重构低频和高频信号。

loadwoman;

%对信号用小波函数sym5进行2尺度分解

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

%对小波分解结构[c,s]的低频系数进行尺度1和尺度2上的重构

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

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

%对分解结构[c,s]的高频系数进行 'h'、'v'、'd' 三个方向尺度2上的重构

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的斜线高频图像');程序运行结果如图11.18所示。图11.18二维小波系数的单支重构例11.5下列程序代码说明如何通过小波分解系数进行小波重构。首先使用函数dwt2对信号进行小波分解，然后利用函数upcoef2对高频和低频系数进行小波重构。

loadwoman;

%使用小波函数db4对信号进行2层分解

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

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

%提取尺度1的低频系数

ca1=appcoef2(c,s,'db4',1);%对尺度1的低频部分进行重构a1=upcoef2('a',ca1,'db4',1,siz);%提取尺度1的水平方向高频系数chd1=detcoef2('h',c,s,1);%对尺度1的水平方向高频部分进行重构hd1=upcoef2('h',chd1,'db4',1,siz);%提取尺度1的垂直方向高频系数cvd1=detcoef2('v',c,s,1);%对尺度1的垂直方向高频部分进行重构vd1=upcoef2('v',cvd1,'db4',1,siz);%提取尺度1的斜线方向高频系数cdd1=detcoef2('d',c,s,1);%对尺度1的斜线方向高频部分进行重构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的低频系数重构图像');subplot(324);image(wcodemat(hd1,nbc));title('尺度1的水平高频系数重构图像');subplot(325);image(wcodemat(vd1,nbc));title('尺度1的垂直高频系数重构图像');subplot(326);image(wcodemat(dd1,nbc));title('尺度1的斜线高频系数重构图像');程序运行如图11.19所示。

图11.19二维小波分解系数的直接重构例11.6下列程序代码实现对小波多尺度分解进行重构，即重构原图像(包括低频和高频)。首先使用函数wavedec2对图像进行多尺度分解，然后调用函数waverec进行图像的重构。

loadwoman;

%使用小波函数sym4对信号进行2尺度分解

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

%进行小波函数的重构a0=waverec2(c,s,'sym4');%检查重构误差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('重构图像');程序运行结果如图11.20所示。图11.20多尺度二维小波重构11.3.2小波图像压缩

MATLAB中的小波图像压缩主要包括获取压缩阈值和进行图像压缩两方面。实现获取压缩阈值的函数有ddencmp和wdcbm2，实现图像压缩的函数有wdencmp、wpdencmp和wthcoef2。

1.小波图像压缩函数

MATLAB中实现信号阈值获取的函数有ddencmp和wdcbm2，下面对它们的用法进行简单的说明。● ddencmp的调用格式有以下三种：

(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)。函数ddencmp用于获取在消噪或压缩过程中的默认阈值。输入参数X为一维或二维信号；IN1取值为 'den' 或 'cmp' ，'den' 表示去噪，'cmp' 表示压缩；IN2取值为 'wv'或 'wp'，'wv' 表示选择小波，'wp' 表示选择小波包。返回值THR是返回的阈值；SORH是软阈值或硬阈值选择参数；KEEPAPP表示保存低频信号；CRIT是熵名(只在选择小波包时用)。●函数wdcbm2的调用格式有以下两种：

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

(2) [THR,NKEEP]=wdcbm2(C,S,ALPHA,M)。函数wdcbm2用于使用Birgé-Massart算法获取二维小波变换的阈值。返回值THR是与尺度无关的阈值，NKEEP是系数的个数。[C，S]是要进行消噪或压缩的信号的在j =length(L)-2层的分解结构；ALPHA和M必须是大于1的实数；THR是关于j的向量，THR(i)是第i层在水平、垂直及斜线方向上的阈值；NKEEP也是关于j的向量，NKEEP(i)是i层的系数个数。一般压缩时ALPHA取1.5，去噪时ALPHA取3。

MATLAB中实现信号的阈值压缩的函数有wdencmp、wthcoef2以及wpdencmp。下面对它们的用法进行简单的介绍。●函数wdencmp的调用格式有下面三种：

(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)。函数wdencmp用于一维或二维信号的消噪或压缩。'gbl'(global的缩写)表示每层都采用同一个阈值进行处理，'wname'是所用的小波函数，'lvd'表示每层用不同的阈值进行处理，N表示小波分解的层数，THR为阈值向量，对于格式(2)和(3)每层都要求有一个阈值，因此阈值向量THR的长度为N，SORH表示选择软阈值或硬阈值(分别取值为's'和'h')，参数KEEPAPP取值1时，则低频系数不进行阈值量化，反之，则低频系数要进行阈值量化。

XC是消噪或压缩后的信号，[CXC,LXC]是XC的小波分解结构，PERF0和PERFL2是恢复和压缩L2范数百分比。PERFL2 = 100 × (XD的小波包系数范数/X的小波包系数)2。如果X是一维信号，小波 'wname' 是一个正交小波，则。●函数wthcoef2的调用格式有下面四种：(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)。函数wthcoef2用于二维信号小波系数的阈值处理。参数 'type' 的取值为 'h( 'v' 或 'd' )时，NC=wthcoef2( 'type',C,S,N,T,SORH)返回小波分解结构[C,S]经由向量N和T定义的阈值量化后的水平(垂直或斜线)系数。N是包含要进行阈值量化的细节尺度向量，T是相应的阈值。N和T的长度必须相等，向量N的长度必须满足1≤N(i)≤size(S,1)-2。其中SORH='s'表示软阈值处理，SORH='h'表示硬阈值处理。函数wpdencmp的调用格式有下面两种：

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

(2) [XD,TREED,PERF0,PERFL2]=wpdencmp(TREE,SORH,CRIT,PAR,KEEPAPP)。函数wpdencmp用于使用小波包变换进行信号的压缩或去噪。格式(1)返回输入信号X(一维或二维)的去噪或压缩后的信号XD。输出参数TREED是XD的最佳小波包分解树；PERFL0和PERFL2是恢复和压缩L2能量百分比。PERFL2=100*(XD的小波包系数范数/X的小波包系数)2。如果X是一维信号，'wname' 是正交小波，则。SORH的取值为 's' 或 'h'，表示软阈值或硬阈值。输入参数N是小波包分解的层数，'wname' 是小波名。函数使用由字符串CRIT定义的熵标准和阈值参数PAR实现最佳分解。如果KEEPAPP=1，则近似信号的小波系数不进行阈值量化；否则，进行量化。格式(2)与格式(1)的输出参数相同，输入选项也相同，只是它从信号的小波包分解树TREE直接进行去噪或压缩。

2.小波图像压缩实例下面将通过例11.7～11.9说明如何进行图像压缩。例11.7下列程序代码使用函数wdcbm2获取图像压缩阈值，然后采用函数wdencmp实现图像压缩：

loadjulia;

nbc=size(map,1);

%使用haar小波进行三层分解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);%对信号进行压缩[xc,cxc,lxc,PERF0,PERFL2]=wdencmp('lvd',c,s,wname,lev,thr,'h');%计算压缩比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('原始信号');subplot(1,2,2)image(wcodemat(XC,nbc))title('压缩后的信号');程序运行结果如图11.21所示。图11.21图像压缩统计小波系数中非零元素值并计算图像压缩比：cmp_rate=0.0950例11.8在本例中使用函数ddencmp获取图像压缩阈值，然后采用函数wdencmp实现图像压缩。程序代码如下：

loadjulia;

nbc=size(map,1);

%使用db2小波进行三层分解

wname='db2';lev=3;[c,s]=wavedec2(X,lev,wname);%使用ddencmp获得压缩阈值[thr,nkeep]=ddencmp('cmp','wv',X)%对信号进行压缩[xc,cxc,lxc,PERF0,PERFL2]=wdencmp('gbl',c,s,wname,lev,thr,'s',1);%计算压缩比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('原始信号');subplot(1,2,2)image(wcodemat(xc,nbc))title('压缩后的信号');程序运行结果如图11.22所示。图11.22图像压缩统计压缩计算前后，小波系数中零元素的个数比例为6.28。可见压缩处理后，小波系数中零元素的个数显著增加，这时可以利用第8章中介绍的行程编码对其进行编码压缩，从而取得较好的压缩效果。例11.9在本例中首先使用函数ddencmp获得图像压缩阈值，然后才用函数wpdencmp对图像进行小波包压缩。程序代码如下：%加载原始信号loadmask;nbc=size(map,1);%使用wpdencmp进行小波压缩%使用ddencmp获取阈值[thr,sorh,keepapp,crit]=ddencmp('cmp','wp',X);%使用阈值进行小波包压缩[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('压缩后的图像');程序运行结果如图11.23所示。图11.23图像小波包压缩11.3.3小波图像去噪

1.图像阈值去噪概述图像去噪与图像压缩有一定的相似性。不同的是图像去噪需要用二维的分析工具对图像进行分析，在对高频系数分别进行阈值量化处理时，要对图像分解的三个方向的高频系数进行阈值处理。图像去噪的一般方法如下：

(1)图像的小波分解。选择合适的小波函数以及适合的分解层次对图像进行分解。

(2)对分解后的高频系数进行阈值处理。对分解的每一层，选择合适的阈值对该层的水平、垂直和斜线三个方向的高频系数进行阈值量化处理。阈值函数和阈值选择规则同信号的阈值消噪相同。

(3)重构图像。根据小波分解的低频系数和经阈值量化处理后的高频系数进行图像重构。

MATLAB中实现了图像的阈值去噪，主要是阈值获取和阈值去噪两方面。

MATLAB中实现阈值获取的函数有ddencmp、thselect、wbmpen和wdcbm2。这些函数在前面已经介绍过。

MATLAB中实现阈值去噪的函数主要有wden、wdencmp、wpdencmp、wthresh、wpthcoef和wthcoef2。其中函数wden、wdencmp、wpdencmp、wthcoef2和wpthcoef的用法前面已介绍，下面只简单介绍一下函数wthresh的用法。函数wthresh的调用格式如下：

Y=wthresh(X,SORH,T)

Y=wthresh(X,SORH,T)返回输入向量或矩阵X经软阈值(SORH='s')或硬阈值(SORH='h')

T处理后的信号。

Y=wthresh(X,'s',T)返回的是Y = sign(X)·(|X| - T)+，即把信号的绝对值与阈值进行比较，小于或等于阈值的点变为0，大于阈值的点变为该点值与阈值的差值。

Y=wthresh(X,'h',T)返回的是Y = X·1(|X| ＞ T)，即把信号的绝对值与阈值比较，小于或等于阈值的点变为0，大于阈值的点保持不变。一般来说，用硬阈值处理后的信号比软阈值处理后的信号更粗糙。

2.图像阈值去噪实例本节通过几个例子说明如何应用MATLAB的函数进行图像的阈值去噪。例11.10使用函数wdcbm2获取图像去噪阈值，然后采用函数wdencmp实现图像去噪。程序代码如下：

loadwoman;

%产生含噪信号

init=2055615866;

randn('seed',init);x=X+18*randn(size(X));nbc=size(map,1);%使用sym4小波进行三层分解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);%对信号进行去噪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('原始信号');subplot(2,2,2);image(wcodemat(x,nbc));title('含噪信号');subplot(2,2,3);image(wcodemat(xd,nbc));title('消噪后的信号');程序运行结果如图11.24所示。

图11.24图像小波阈值去噪例11.11在本例中首先使用函数wpcoef获取图像的小波包分解系数，接着估计噪声方差，然后使用函数wpbmpen获取小波包去噪阈值，最后使用函数wpdencmp实现图像去噪。程序代码如下：

loadwoman;

%产生含噪信号

init=2055615866;

randn('seed',init);

x=X+18*randn(size(X));%对图像进行小波包分解%使用coif2分解到第三层wname='coif2';lev=3;tree=wpdec2(x,lev,wname);%估计噪声标准差det1=[wpcoef(tree,2)wpcoef(tree,3)wpcoef(tree,4)];sigma=median(abs(det1(:)))/0.6745;%选择去噪的阈值alpha=1.1;thr=wpbmpen(tree,sigma,alpha);%使用上面的阈值对图像去噪keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp);%画出原始图像和去噪后的图像nbc=size(map,1);colormap(pink(nbc));subplot(2,2,1);image(wcodemat(X,nbc));title('原始信号');subplot(2,2,2);image(wcodemat(x,nbc));title('含噪信号');subplot(2,2,3);image(wcodemat(xd,nbc));title('消噪后的信号');