第8章 小波与小波变换小波分析是近十几年才发展起来并迅速应用到图像处理和语音分析等众多领域的一种数学工具。它是继110多年前的傅立叶(Joseph Fourier)分析之后的一个重大突破,无论是对古老的自然学科还是对新兴的高新技术应用学科都产生了强烈冲击。
小波理论是应用数学的一个新领域。要深入理解小波理论需要用到比较多的数学知识。本章企图从工程应用角度出发,用比较直观的方法来介绍小波变换和它的应用,为读者深入研究小波理论和应用提供一些背景材料。
8.1 小波介绍
8.1.1 小波简史傅立叶理论指出,一个信号可表示成一系列正弦和余弦函数之和,叫做傅立叶展开式。用傅立叶表示一个信号时,只有频率分辨率而没有时间分辨率,这就意味我们可以确定信号中包含的所有频率,但不能确定具有这些频率的信号出现在什么时候。为了继承傅立叶分析的优点,同时又克服它的缺点,人们一直在寻找新的方法。
20世纪初,哈尔(Alfred Haar)对在函数空间中寻找一个与傅立叶类似的基非常感兴趣。1909年他发现了小波,并被命名为哈尔小波(Haar wavelets),他最早发现和使用了小波。
20世纪70年代,当时在法国石油公司工作的年轻的地球物理学家Jean Morlet提出了小波变换WT(wavelet transform)的概念。
进入20世纪80年代,法国的科学家Y.Meyer和他的同事开始为此开发系统的小波分析方法。Meyer于1986年创造性地构造出具有一定衰减性的光滑函数,他用缩放(dilations)与平移(translations)均为(的整数)的倍数构造了空间的规范正交基,使小波得到真正的发展。
小波变换的主要算法则是由法国的科学家Stephane Mallat在1988年提出[1]。他在构造正交小波基时提出了多分辨率的概念,从空间上形象地说明了小波的多分辨率的特性,提出了正交小波的构造方法和快速算法,叫做Mallat算法[1]。该算法统一了在此之前构造正交小波基的所有方法,它的地位相当于快速傅立叶变换在经典傅立叶分析中的地位。
Inrid Daubechies,Ronald Coifman和Victor Wickerhauser等著名科学家把这个小波理论引入到工程应用方面做出了极其重要的贡献。例如,Inrid Daubechies于1988年最先揭示了小波变换和滤波器组(filter banks)之间的内在关系[2],使离散小波分析变成为现实。在信号处理中,自从S.Mallat和Inrid Daubechies发现滤波器组与小波基函数有密切关系之后,小波在信号(如声音信号,图像信号等)处理中得到极其广泛的应用。
经过十几年的努力,这门学科的理论基础已经基本建立,并成为应用数学的一个新领域。这门新兴学科的出现引起了许多数学家和工程技术人员的极大关注,是国际科技界和众多学术团体高度关注的前沿领域。
8.1.2 小波概念小波是定义在有限间隔而且其平均值为零的一种函数,它的波形如图8-01(b)所示。图8-01(a)是大家所熟悉的正弦波,图8-01(b)是从许多使用比较广泛的小波中挑选出的几种一维小波。
在图8-01(b)所示的小波中,缩放函数和小波函数的名称大多数是以开发者的名字命名的.。例如,Moret小波函数是Grossmann和Morlet在1984年开发的,db6缩放函数和db6小波函数是Daubechies开发的开发几种小波之一;Meyer缩放函数和Meyer小波函数是Meyer开发的。但也有不少例外,例如Sym6缩放函数和sym6小波函数则是symlets的简写,是Daubechies提议开发的几种对称小波之一;coif2缩放函数和coif2小波函数是Daubechies应R,Coifman的请求而开发的几种小波之一。
与图8-01(a)相比,图8-01(b)所示的小波具有有限的持续时间和突变的频率和振幅,波形可以是不规则的,也可以是不对称的,在整个时间范围里的幅度平均值为零。而正弦波和余弦波具有无限的持续时间,它可从负无穷扩展到正无穷,波形是平滑的,它的振幅和频率也是恒定的。


(b) 部分小波图8-01 正弦波与小波在众多的小波中,选择什么样的小波对信号进行分析是一个至关重要的问题。使用的小波不同,分析得到数据也不同,这是关系到能否达到使用小波分析的目的问题。如果没有现成的小波可用,那么还需要自己开发合用的小波。
顺便要提及的是,小波函数在时域和频域中都应该具有某种程度的平滑度(smoothness)和集中性(concentration),这个复杂的概念在数学上使用消失矩(vanishing moments)来描述,用N表示小波的消失矩的数目。例如,Daubechies小波简写成dbN,如db1,db2,……,db9,从Daubechies小波波形来看,N数目的大小反映了Daubechies小波的平滑度和集中性。
8.1.3 小波分析信号分析一般是为了获得时间和频率域之间的相互关系。傅立叶变换提供了有关频率域的信息,但时间方面的局部化信息却基本丢失。与傅立叶变换不同,小波变换通过平移母小波(mother wavelet)可获得信号的时间信息,而通过缩放小波的宽度(或者叫做尺度)可获得信号的频率特性。对母小波的缩放和平移操作是为了计算小波的系数,这些系数代表小波和局部信号之间的相互关系。本节将介绍小波分析中常用的三个基本概念:连续小波变换、离散小波变换和小波重构。
1,连续小波变换傅立叶分析是把一个信号分解成各种不同频率的正弦波,因此正弦波是傅立叶变换的基函数。同样,小波分析是把一个信号分解成将原始小波经过移位和缩放之后的一系列小波,因此小波同样可以用作表示一些函数的基函数。可以说,凡是能够用傅立叶分析的函数都可以用小波分析,因此小波变换也可以理解为用经过缩放和平移的一系列函数代替傅立叶变换的正弦波。
仔细观察图8-02所示的正弦波和小波可以发现,用不规则的小波分析变化激烈的信号也许比用平滑的正弦波更有效,或者说对信号的基本特性描述得更好。

图8-02 傅立叶分析与小波分析使用的基函数数学上傅立叶分析的过程实际上是用傅立叶变换表示,

这个式子的含义就是,傅立叶变换是信号与复数指数()之积在信号存在的整个期间里求和。傅立叶变换的结果是傅立叶系数,它是频率的函数。
同样,连续小波变换(continuous wavelet transform,CWT )用下式表示,

这个式子的含义就是,小波变换是信号与被缩放和平移的小波函数之积在信号存在的整个期间里求和。CWT变换的结果是许多小波系数,这些系数是缩放因子(scale)和位置(position)的函数。
对缩放因子可这样来理解。如果用字母表示缩放因子,例如,对于正弦函数,
; 它的缩放因子
它的缩放因子
它的缩放因子
CWT的变换过程可分成如下5个步骤:
步骤1,把小波和原始信号的开始部分进行比较。
步骤2,计算系数。该系数表示该部分信号与小波的近似程度。系数的值越高表示信号与小波越相似,因此系数可以反映这种波形的相关程度。
步骤3,把小波向右移,距离为,得到的小波函数为,然后重复步骤1和2。再把小波向右移,得到小波,重复步骤1和2。按上述步骤一直进行下去,直到信号结束。
步骤4,扩展小波,例如扩展一倍,得到的小波函数为。
步骤5,重复步骤1~4。
CWT的整个变换过程如图8-03所示。

图8-03 连续小波变换的过程小波变换完成之后得到的系数是在不同的缩放因子下由信号的不同部分产生的。这些小波系数、缩放因子和时间之间的关系和它们的含义可以用图8-04(a)表示,该图是用MATLAB软件绘制的。图8-04(a)是用二维图像表示的小波变换分析图,轴表示沿信号的时间方向上的位置,轴表示缩放因子,每个点的颜色表示小波系数的幅度大小。图8-04(b)是用三维图像表示的小波变换分析图,轴表示小波变换之后的系数。

(a) 二维图

(b) 三维图图8-04 连续小波变换分析图[4]
小波的缩放因子与信号频率之间的关系可以这样来理解。缩放因子小,表示小波比较窄,度量的是信号细节,表示频率比较高;相反,缩放因子大,表示小波比较宽,度量的是信号的粗糙程度,表示频率比较低。
2,离散小波变换在计算连续小波变换时,实际上也是用离散的数据进行计算的,只是所用的缩放因子和平移参数比较小而已。不难想象,连续小波变换的计算量是惊人的。为了解决计算量的问题,缩放因子和平移参数都选择(.>0的整数)的倍数。使用这样的缩放因子和平移参数的小波变换叫做双尺度小波变换(dyadic wavelet transform),它是离散小波变换(discrete wavelet transform,DWT)的一种形式。从文献看,离散小波变换通常指的就是双尺度小波变换。
使用离散小波分析得到的小波系数、缩放因子和时间关系如图8-05所示。图(a)是20世纪40年代使用Gabor开发的短时傅立叶变换(short time Fourier transform,STFT)得到的时间-频率关系图,图(b)是20世纪80年代使用Morlet开发的小波变换得到的时间-缩放因子(反映频率)关系图。

图8-05 离散小波变换分析图执行离散小波变换的有效方法是使用滤波器。该方法是Mallat在1988年开发的,叫做Mallat算法[1]。这种方法实际上是一种信号的分解方法,在数字信号处理中称为双通道子带编码。
用滤波器执行离散小波变换的概念如图8-06所示。图中,S表示原始的输入信号,通过两个互补的滤波器产生A和D两个信号,A表示信号的近似值(approximations),D表示信号的细节值(detail)。在许多应用中,信号的低频部分是最重要的,而高频部分起一个“添加剂”的作用。犹如声音那样,把高频分量去掉之后,听起来声音确实是变了,但还能够听清楚说的是什么内容。相反,如果把低频部分去掉,听起来就莫名其妙。在小波分析中,近似值是大的缩放因子产生的系数,表示信号的低频分量。而细节值是小的缩放因子产生的系数,表示信号的高频分量。

图8-06 双通道滤波过程由此可见,离散小波变换可以被表示成由低通滤波器和高通滤波器组成的一棵树。原始信号通过这样的一对滤波器进行的分解叫做一级分解。信号的分解过程可以叠代,也就是说可进行多级分解。如果对信号的高频分量不再分解,而对低频分量连续进行分解,就得到许多分辨率较低的低频分量,形成如图8-07所示的一棵比较大的树。这种树叫做小波分解树(wavelet decomposition tree)。分解级数的多少取决于要被分析的数据和用户的需要。

(a)信号分解 (b)系数结构 (c)小波分解树图8-07 小波分解树小波分解树表示只对信号的低频分量进行连续分解。如果不仅对信号的低频分量连续进行分解,而且对高频分量也进行连续分解,这样不仅可得到许多分辨率较低的低频分量,而且也可得到许多分辨率较低的高频分量。这样分解得到的树叫做小波包分解树(wavelet packet decomposition tree),这种树是一个完整的二进制树。图8-08表示的是一棵三级小波包分解树。小波包分解方法是小波分解的一般化,可为信号分析提供更丰富和更详细的信息。例如,小波包分解树允许信号S表示为


图8-08 三级小波包分解树随便要提及的是,在使用滤波器对真实的数字信号进行变换时,得到的数据将是原始数据的两倍。例如,如果原始信号的数据样本为1000个,通过滤波之后每一个通道的数据均为1000个,总共为2000个。于是,根据尼奎斯特(Nyquist)采样定理就提出了降采样(downsampling)的方法,即在每个通道中每两个样本数据取一个,得到的离散小波变换的系数(coefficient)分别用cD和cA表示,如图8-09所示。图中的符号表示降采样。

图8-09 降采样过程
3,小波重构离散小波变换可以用来分析或者叫做分解信号,这个过程叫做分解或者叫做分析。把分解的系数还原成原始信号的过程叫做小波重构(wavelet reconstruction)或者叫做合成(synthesis),数学上叫做逆离散小波变换(inverse discrete wavelet transform,IDWT)。
在使用滤波器做小波变换时包含滤波和降采样两个过程,在小波重构时要包含升采样(upsampling)和滤波过程。小波重构的方法如图8-10所示,图中的符号表示升采样。

图8-10 小波重构方法升采样是在两个样本数据之间插入“0”,目的是把信号的分量加长。升采样的过程如图8-11所示。

图8-11 升采样的方法重构过程中滤波器的选择也是一个重要的研究问题,这是关系到能否重构出满意的原始信号的问题。在信号的分解期间,降采样会引进畸变,这种畸变叫做混叠(aliasing)。这就需要在分解和重构阶段精心选择关系紧密但不一定一致的滤波器才有可能取消这种混叠。低通分解滤波器(L)和高通分解滤波器(H)以及重构滤波器(L'和H')构成一个系统,这个系统叫做正交镜像滤波器(quadrature mirror filters,QMF)系统,如图8-12所示。

图8-12 正交镜像滤波器系统
8.1.4 小波定义在数学上,小波定义为对给定函数局部化的函数。小波可由一个定义在有限区间的函数来构造,称为母小波(mother wavelet)或者叫做基本小波。一组小波基函数,,可通过缩放和平移基本小波来生成,

其中,为进行缩放的缩放参数,反映特定基函数的宽度(或者叫做尺度);为进行平移的平移参数,指定沿轴平移的位置。
当和的情况下,一维小波基函数序列定义为
,或者
本教材将采用下面的表示法,

其中,为平移参数,为缩放因子。
函数以小波为基的连续小波变换定义为函数和的内积,

在1984年,A.Grossman和J.Morlet指出,连续小波的逆变换为,

其中,为母小波的允许条件(admissible condition),

其中,为的傅立叶变换,而是在平方可积的实数空间。
8.2 哈尔函数哈尔小波是小波系列中最简单的小波,因此本节将从哈尔小波入手,首先介绍用来构造任意给定信号的哈尔基函数,然后介绍用来表示任意给定信号的哈尔小波函数,最后介绍函数的规范化和哈尔基的构造。
8.2.1 哈尔基函数基函数是一组线性无关的函数,可以用来构造任意给定的信号。例如,用基函数的加权和表示。最简单的基函数是哈尔基函数(Haar basis function)。哈尔基函数在1909年提出,它是由一组分段常值函数(piecewise-constant function)组成的函数集。这个函数集定义在半开区间上,每一个分段常值函数的数值在一个小范围里是“1”,其他地方为“0”,现以图像为例并使用线性代数中的矢量空间来说明哈尔基函数。
如果一幅图像仅由=1个像素组成,这幅图像在整个区间中就是一个常值函数。用表示这个常值函数,用表示由这个常值函数生成的矢量空间,即
:
它的波形如图8-13所示。

图8-13 的波形这个常值函数也叫做框函数(box function),它是构成矢量空间的基。
如果一幅图像由=2个像素组成,这幅图像在区间中有两个等间隔的子区间:和,每一个区间中各有1个常值函数,分别用表示。用表示由2个子区间中的常值函数生成的矢量空间,即
:,
它们的波形如图8-14所示。

图8-14 的波形这2个常值函数就是构成矢量空间的基。
如果一幅图像由个像素组成,这幅图像在区间中被分成个等间隔的子区间:,它们的常值函数分别用,,和表示,用表示由4个子区间中的常值函数生成的矢量空间,即

它们的波形如图8-15所示。

图8-15 ,,和的波形这4个常值函数就是构成矢量空间的基。
我们可以按照这种方法继续定义基函数和由它生成的矢量空间。总而言之,为了表示矢量空间中的矢量,每一个矢量空间都需要定义一个基(basis)。为生成矢量空间而定义的基函数也叫做尺度函数(scaling function),这种函数通常用符号表示。哈尔基函数定义为
 (8-1)
哈尔基尺度函数定义为
 (8-2)
其中,为尺度因子,改变使函数图形缩小或者放大;为平移参数,改变使函数沿轴方向平移。
空间矢量定义为
 (8-3)
其中,表示线性生成(linear span)。
需要注意的是,有些文章使用负整数来定义尺度函数并且使用不同的符号规则。因此在阅读有关小波方面的文章时要注意作者使用的符号规则。
由于定义了基和矢量空间,就可以把由个像素组成的一维图像看成为矢量空间中的一个矢量。由于这些矢量都是在单位区间上定义的函数,所以在矢量空间中的每一个矢量也被包含在矢量空间中。这说明矢量空间是嵌套的,即

矢量空间的这个性质可写成

8.2.2 哈尔小波函数小波函数通常用表示。与框函数相对应的小波称为基本哈尔小波函数(Haar wavelet functions),并由下式定义,
 (8-4)
哈尔小波尺度函数定义为
 (8-5)
用小波函数构成的矢量空间用表示,
 (8-6)
其中,表示线性生成;为尺度因子,改变使函数图形缩小或者放大;为平移参数,改变使函数沿轴方向平移。
根据哈尔小波函数的定义,可以写出生成,和等矢量空间的小波函数。
生成矢量空间的哈尔小波:

它的波形如图8-16所示。

图8-16 的波形生成矢量空间的哈尔小波:
 
它们的波形如图8-17所示。

图8-17 的波形生成矢量空间的哈尔小波:

它们的波形如图8-18所示。

图8-18 、、和的波形用哈尔小波生成的矢量空间包含在矢量空间,这个性质用下式表示,

8.2.3 函数的规范化在小波变换中,有时要对基函数和小波函数进行规范化(normalization)。在半开区间中,如果函数的内积
 (8-7)
 (8-8)
则称是规范化的函数。哈尔基和哈尔小波分别使用下面两个式子进行规范化,
 (8-9)
 (8-10)
其中,常数因子用来满足标准内积(inner product)等于1的条件。如果小波函数不是在区间中定义的函数,常数因子将改变。
8.2.4 哈尔基的结构使用哈尔基函数和哈尔小波函数生成的矢量空间和具有下面的性质,
 (8-11)
其中,符号“”表示直和。这就是说,在矢量空间中,根据所选择的内积,生成矢量空间的所有函数与生成矢量空间的所有函数都是正交的,即子空间是子空间的正交补(orthogonal complement)。式(8-11)表明,在矢量空间中,矢量空间中的小波可用来表示一个函数在矢量空间中不能表示的部分。因此,小波可定义为生成矢量空间的一组线性无关的函数的集合。这些基函数具有两个重要性质:
生成矢量空间的基函数与生成矢量空间的基函数构成矢量空间的一个基。
生成矢量空间的每一个基函数与生成矢量空间的每一个基函数正交。
8.3 哈尔小波变换小波变换的基本思想是用一组小波函数或者基函数表示一个函数或者信号,例如图像信号。为了理解什么是小波变换,下面用一个具体的例子来说明小波变换的过程。
1,求有限信号的均值和差值
[例8,1] 假设有一幅分辨率只有4个像素的一维图像,对应的像素值或者叫做图像位置的系数分别为:
[9 7 3 5]
计算它的哈尔小波变换系数。
计算步骤如下:
步骤1:求均值(averaging)。计算相邻像素对的平均值,得到一幅分辨率比较低的新图像,它的像素数目变成了2个,即新的图像的分辨率是原来的1/2,相应的像素值为:
[8 4]
步骤2:求差值(differencing)。很明显,用2个像素表示这幅图像时,图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像素组成的原始图像,就需要存储一些图像的细节系数(detail coefficient),以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像素对的平均值,或者使用这个像素对的差值除以2。在这个例子中,第一个细节系数是(9-8)=1,因为计算得到的平均值是8,它比9小1而比7大1,存储这个细节系数就可以恢复原始图像的前两个像素值。使用同样的方法,第二个细节系数是(3-4) = -1,存储这个细节系数就可以恢复后2个像素值。因此,原始图像就可以用下面的两个平均值和两个细节系数表示,
[8 4 1 -1]
步骤3:重复步骤1和2,把由第一步分解得到的图像进一步分解成分辨率更低的图像和细节系数。在这个例子中,分解到最后,就用一个像素的平均值6和三个细节系数2,1和-1表示整幅图像:
[6 2 1 -1]
这个分解过程如表8-1所示。
表8-1 哈尔变换过程分辨率
平均值
细节系数
4
[9 7 3 5]
2
[8 4]
[1 -1]
1
[6]
[2]
由此可见,通过上述分解就把由4像素组成的一幅图像用一个平均像素值和三个细节系数表示,这个过程就叫做哈尔小波变换(Haar wavelet transform),也称哈尔小波分解(Haar wavelet decomposition)。这个概念可以推广到使用其他小波基的变换。
从这个例子中我们可以看到:
变换过程中没有丢失信息,因为能够从所记录的数据中重构出原始图像。
对这个给定的变换,我们可以从所记录的数据中重构出各种分辨率的图像。例如,在分辨率为1的图像基础上重构出分辨率为2的图像,在分辨率为2的图像基础上重构出分辨率为4的图像。
通过变换之后产生的细节系数的幅度值比较小,这就为图像压缩提供了一种途径。例如,去掉一些微不足道的细节系数并不影响对重构图像的理解。
2,哈尔小波变换在例8.1中,求均值和差值的过程实际上就是一维小波变换的过程,现在用数学方法重新描述小波变换的过程。
(1) 用中的哈尔基表示图像=[9 7 3 5]有像素,因此可以用生成矢量空间中的框基函数的线性组合表示,

其中的系数是4个正交的像素值[9 7 3 5],因此,

可用图8-19表示。

图8-19 用中的哈尔基表示
(2) 用和中的函数表示生成矢量空间的基函数为和,生成矢量空间的小波函数为和。根据式(8-11),

因此,可表示成

其中,系数和是分辨率为2时的像素平均值,和为细节系数;和是中的哈尔基函数,和是中的哈尔小波函数。可用图8-20表示。

图8-20 用和中的函数表示
(3) 用和中的函数表示生成矢量空间的基函数为,生成矢量空间的小波函数为,生成矢量空间的小波函数为和。根据(8-11),

可表示成

并可用图8-21表示。

图8-21 用和中的函数表示其中,4个系数,,和就是原始图像通过哈尔小波变换所得到的系数,用来表示整幅图像的平均值和不同分辨率下的细节系数。4个函数,,和就是构成空间的基。
8.4 规范化算法规范化的小波变换与非规范化的小波变换相比,唯一的差别是在规范化的小波变换中需要乘一个规范化的系数。下面用一个例子说明。
[例8.2] 对函数作哈尔小波变换。
哈尔小波变换实际上是使用求均值和差值的方法进行分解。我们把看成是矢量空间中的一个向量,尺度因子,因此最多可分解为3层,如图8-22所示。

图8-22 小波分解的层次结构分解过程如下。
步骤1:

步骤2:

步骤3:

根据这个例子,我们可以归纳出规范化的哈尔小波变换的一般算法。假设一维阵列有个元素,等于2的幂,执行一维哈尔小波变换的伪代码如下:
******************************************************************************
proc DecomposeArray(C, array of color):
while h > 1 do,

for  to  do:


end for

end while
end proc
******************************************************************************
8.5 二维哈尔小波变换前面已经介绍了一维小波变换的基本原理和变换方法。这节将结合具体的图像数据系统地介绍如何使用小波对图像进行变换。
我们已经知道,一幅图像可看成是由许多像素组成的一个大矩阵,在进行图像压缩时,为降低对存储器的要求,人们通常把它分成许多小块。例如以8×8个像素为一块,并用矩阵表示,然后分别对每一个图像块进行处理。在小波变换中,由于小波变换中使用的基函数的长度是可变的,因此虽然无须像以离散余弦变换(DCT)为基础的JPEG标准算法那样,把输入图像进行分块,以避免产生JPEG图像那样的“块效应”,但为便于理解小波变换的奥妙,还是从一个小的图像块入手,并且继续使用哈尔小波对图像进行变换。
8.5.1 二维小波变换举例假设有一幅灰度图像,其中的一个图像块用矩阵表示为,

使用灰度表示的图像如图8-23所示。

图8-23 图像矩阵A的灰度图一个图像块是一个二维的数据阵列,进行小波变换时可以对阵列的每一行进行变换,然后对行变换之后的阵列的每一列进行变换,最后对经过变换之后的图像数据阵列进行编码。
1,求均值与求差值为读者对用小波变换压缩图像有一个完整的概念,还是从求均值(averaging)与求差值(differencing)开始。在图像块矩阵中,第一行的像素值为,
R0,[64 2 3 61 60 6 7 57]
步骤1:在R0行上取每一对像素的平均值,并将结果放到新一行N0的前4个位置,其余的4个数是R0行每一对像素的第一个数与其相应的平均值之差。这个变换过程如下所示,
R0,[64 2 3 61 60 6 7 57]
   
求差值,
N0,[33 32 33 32 31 -29 27 -25]
64-33=31 3-32=-29 60-33=27 7-32=-25
步骤2:对行N0的前4个数使用与第一步相同的方法,得到两个平均值和两个系数,并放在新一行N1的前4个位置,其余的4个细节系数直接从行N0复制到N1的相应位置上。整个过程如下所示,
N0,[33 32 33 32 31 -29 27 -25]
 
求差值,
N1,[32.5 32.5 0.5 0.5 31 -29 27 -25]
33-32.5=0.5 33-32.5=0.5
步骤3:用与步骤1和2相同的方法,对剩余的一对平均值求平均值和差值,
N1,[32.5 32.5 0.5 0.5 31 -29 27 -25]

求差值,
N1,[32.5 0 0.5 0.5 31 -29 27 -25]
32.5-32.5=0
2,图像矩阵的计算使用求均值和求差值的方法,对矩阵的每一行进行计算,得到,
=
其中,每一行的第一个元素是该行像素值的平均值,其余的是这行的细节系数。使用同样的方法,对的每一列进行计算,得到,
=
其中,左上角的元素表示整个图像块的像素值的平均值,其余是该图像块的细节系数。根据这个事实,如果从矩阵中去掉表示图像的某些细节系数,事实证明重构的图像质量仍然可以接受。具体做法是设置一个阈值,例如的细节系数就把它当作“0”看待,这样经过变换之后的矩阵就变成,
=
与相比,“0”的数目增加了18个,也就是去掉了18个细节系数。这样做的好处是可提高小波图像编码的效率。对矩阵进行逆变换,得到了重构的近似矩阵,

矩阵的数据用图8-24(a)表示,矩阵的数据用图8-24(b)表示。对比图8-24(b)和图8-24(a),如果不事先告诉你,图8-24(a)是原图,而图8-24(b)是经过变换并且去掉某些细节系数之后重构的图,也许你很难断定那一幅是原图,那一幅是重构图。这说明图像质量的损失还是能够接受的。

图8-24 原图与重构图像的比较
3,使用线性代数由于图像可用矩阵表示,使用三个矩阵和同样可以对矩阵进行求平均值和求差值。这三个矩阵分别是第一、第二和第三次分解图像时所构成的矩阵。现以矩阵的第一行的数据为例,
[64 2 3 61 60 6 7 57]
步骤1,使用计算
[64 2 3 61 60 6 7 57]
=[64 2 3 61 60 6 7 57]
=[33 32 33 32 31 -29 27 -25]
步骤2,使用计算
[33 32 33 32 31 -29 27 -25]
=[33 32 33 32 31 -29 27 -25]
=[32.5 32.5 0.5 0.5 31 -29 27 -25]
步骤3,使用计算
[32.5 32.5 0.5 0.5 31 -29 27 -25]
=[32.5 32.5 0.5 0.5 31 -29 27 -25]
=[32.5 0 0.5 0.5 31 -29 27 -25]
为简化计算,可把三个矩阵相乘之后再计算。三个矩阵相乘的结果为,

如果使用规格化的小波进行计算,则

以上计算哈尔小波正变换的步骤就简化成,

计算哈尔小波的逆变换简化为,

4,变换实例为进一步理解小波变换的基本原理和在图像处理中的应用,我们可使用MATLAB软件中的小波变换工具箱(Wavelet Toolbox)编写小波变换程序,对原始图像进行分解和重构。
图8-25表示图像分解和重构过程。利用小波变换,用户可以按照应用要求获得不同分辨率的图像。如图8-26所示,图8-26(a)表示原始的Lena图像,图8-26(b)表示通过一级(level)小波变换可得到1/4分辨率的图像,图8-26(c)表示通过二级小波变换可得到1/8分辨率的图像,图8-26(d)表示通过三级小波变换可得到1/16分辨率的图像。

图8-25小波图像分解与重构

图8-26使用小波分解产生多种分辨率图像
阈值处理可用于去除图像中的噪声。在取不同阈值的情况下重构图像,可观察到图像质量发生的变化,如图8-27所示,图8-27(a)表示原始的Lena图像,图8-27(b)表示阈值5的重构图像,图8-27(c)表示阈值10的重构图像,而图8-27(d)表示阈值20的重构图像。从图中可以看到,随着阈值的增大,图像质量也随之降低。

(a) 原始图像 (b)5

(c) 10 (d)20
图8-27 不同阈值下的Lena图像
8.5.2 二维小波变换方法用小波对图像进行变换有两种方法,一种叫做标准分解(standard decomposition),另一种叫做非标准分解(nonstandard decomposition)。标准分解方法是指首先使用一维小波对图像每一行的像素值进行变换,产生每一行像素的平均值和细节系数,然后使用一维小波对这个经过行变换的图像的列进行变换,产生这个图像的平均值和细节系数。标准分解的过程如下,
******************************************************************************
procedure StandardDecomposition(C,array [1.,, h,1.,, w] of reals)
for row 1 to h do
Decomposition(C [row,1,,, w])
end for
for col 1 to w do
Decomposition(C [1,,, h,col])
end for
end procedure
******************************************************************************
图8-28表示标准分解方法,图像是使用MATLAB编写的程序分解得到的。

图8-28 图像的标准分解方法非标准分解是指使用一维小波交替地对每一行和每一列像素值进行变换。首先对图像的每一行计算像素对的均值和差值,然后对每一列计算像素对的均值和差值。这样得到的变换结果只有1/4的像素包含均值,再对这1/4的均值重复计算行和列的均值和差值,依此类推。非标准分解的过程如下:
*************************************************************************
procedure NonstandardDecomposition(C,arrayof reals)
 (normalize input coefficients)
while h > 1 do
for row 1 to h do
DecompositionStep(C [row,1,,, h])
end for
for col 1 to h do
DecompositionStep(C [1,,, h,col])
end for

end while
end procedure
**************************************************************************
图8-29表示非标准分解方法,图像是使用MATLAB编写的程序分解得到的。

图8-29 图像的非标准分解方法标准分解方法和非标准分解相比,它们得到的变换结果是完全相同的,只是非标准算法的计算量少可以一些。
练习与思考题写出矢量空间的哈尔小波并画出它的波形。
写出4×4哈尔小波变换矩阵。
使用MATLAB中的多级一维小波分解函数例程(function)wavedec,对例8.2所示的函数作小波变换。
使用规范化的小波变换算法,用MATLAB编写一个M文件,重新计算的哈尔小波变换。
参考文献和站点
Mallat,S,G.,A Theory for Multiresolution Signal Decomposition,The Wavelet Representation,IEEE Trans,PAMI,vol,11,no,7,July 1989,pp,674-693.
Daubechies,I,Orthonormal Bases of Compactly Supported Wavelets,Comm,Pure and Applied Math.,vol,41,Nov,1988,pp,909-996.
Eric J,Stollnitz,Tony D,DeRose,and David H,Salesin,Wavelets for computer graphics,A primer,part 1,IEEE Computer Graphics and Applications,15(3),76–84,May 1995.
http://www.mathworks.com/access/helpdesk/help/toolbox/wavelet/wavelet_tocframe.shtml (浏览日期:2001年6月16日)
Eric?J,Stollnitz,Tony?D,DeRose,and David?H,Salesin.,Wavelets for Computer Graphics,A Primer,IEEE?Computer Graphics and Applications,15(3),76-84,May?1995 (part?1),and 15(4),75-85,July?1995 (part?2).
Peggy Morton and Arne Petersen,Image Compression Using the Haar Wavelet Transform,Math 45,College of the Redwoods,Dev.19,1997.
http://amath.colorado.edu/courses/4720/2000Spr/Labs/openingpage.html (浏览日期:2001年7月1日)
Sarkar,T.K.,Su,C.,Adve,R.,Salazar-Palma,M.,Garcia-Castillo,L.,and BOIX,R.R.,A Tutorial on Wavelets from an Electrical Engineering Perspective,Part 1,Discrete Wavelet Techniques,IEEE Antennas & Propagation Magazine,Vol,40,No,5,pp,49-70,Oct,1998,
SARKAR,T.K.,SU,C.,A Tutorial on Wavelets from an Electrical Engineering Perspective,Part 2,The Continuous Case,IEEE Antennas & Propagation Magazine,Vol,40,No,6,pp,36-48,Dec,1998,