小波分析在图像处理中的应用

第二章　小波分析应用示例

本章将探讨小波分析在实际中的一些应用[7,8,9]。

第一节　理解小波分析

小波变换是克服其他信号处理技术缺陷的一种分析信号的方法。小波由一族小波基函数构成，它可以描述信号时间（空间）和频率（尺度）域的局部特性。采用小波分析最大优点是可对信号进行实施局部分析，可在任意的时间或空间域中分析信号。小波分析具有发现其他信号分析方法所不能识别的、隐藏于数据之中的表现结构特性的信息，而这些特性对机械故障和材料的损伤等识别是尤为重要的。如何选择小波基函数目前还没有一个理论标准，常用的小波函数有Haar、Daubechies(dbN)、Morlet、Meryer、Symlet、Coiflet、Biorthogonal小波等15种。但是小波变换的小波系数为如何选择小波基函数提供了依据。小波变换后的系数比较大，就表明了小波和信号的波形相似程度较大；反之则比较小。另外还要根据信号处理的目的来决定尺度的大小。如果小波变换仅仅反映信号整体的近似特征，往往选用较大的尺度；反映信号细节的变换则选用尺度不大的小波。由于小波函数家族成员较多，进行小波变换目的各异，目前没有一个通用的标准。

根据实际运用的经验，Morlet小波应用领域较广，可以用于信号表示和分类、图像识别、特征提取；墨西哥草帽小波用于系统识别；样条小波用于材料探伤；Shannon正交基用于差分方程求解。

现在对小波分解层数与尺度的关系作如下解释：

是不是小波以一个尺度分解一次就是对小波进行一层的分解？

比如：[C,L]=wavedec(X,N,'wname')中，N为尺度，若为1，就是进行单尺度分解，也就是分解一层。但是W=CWT(X,[2:2:128],'wname','plot')的分解尺度又是从2～128以2为步进的，这里的“分解尺度”跟上面那个“尺度”的意思一样吗？

[C,L]=wavedec(X,N,'wname')中的 N 为分解层数，不是尺度，以 wname是DB小波为例，如DB4，4为消失矩，则一般滤波器长度为8，阶数为7。

wavedec针对于离散，CWT是连续的。

多尺度又是怎么理解的呢？

多尺度的理解：如将0～pi定义为空间V0，经过一级分解之后V0被分成0～pi/2的低频子空间V1和pi/2～pi的高频子空间W1，然后一直分下去……得到VJ+WJ+…+W2+W1。因为VJ和WJ是正交的空间，且各 W 子空间也是相互正交的，所以分解得到了是相互不包含的多个频域区间，这就是多分辨率分析，即多尺度分析。

当然多分辨率分析是有严格数学定义的，但完全可以从数字滤波器角度理解它，也可以从函数空间角度理解。

是不是说分解到W3、W2、W1、V3就是三尺度分解？

简单的说尺度就是频率，不过是反比的关系。确定尺度关键还要考虑你要分析信号的采样频率大小，因为根据采样频率大小才能确定你的分析频率是多少。（采样定理）。然后再确定你到底分多少层。

假如有一个10Hz和50Hz的正弦混合信号，采样频率是500Hz，是不是就可以推断出10Hz和50Hz各自对应的尺度了呢？意思是，是不是有一个频率和尺度的换算公式？实际频率＝小波中心频率×采样频率/尺度。

在小波分解中，若将信号中的最高频率成分看作是1，则各层小波小波分解便是带通或低通滤波器，且各层所占的具体频带为（三层分解）a1:0～0.5，d1: 0.5～1; a2:0～0.25，d2: 0.25～0.5; a3: 0～0.125，d3:0.125～0.25。

如果要得到频率为0.125～0.25的信号信息，是不是直接对d3的分解系数直接重构之后就是时域信息了？这样感觉把多层分解纯粹当作滤波器来用了，又怎么是多分辨分析？怎样把时频信息同时表达出来？

确实是把它当成了滤波器来用了，也就是说只看重了小波分析的频域局部化的特性。但是很多人都忽略其时域局部化特性，因为小波是变时频分析的方法，根据测不准原理如果带宽大，则时窗宽度就要小。那么也就意味着如果要利用其时域局部化特性就得在时宽小的分解层数下研究，也就是低尺度下。这样就可以更容易看出信号在该段时间内的细微变化，但是就产生一个问题，这一段的频率带很宽，频率局部化就体现不出来了。

对d3进行单支重构就可以得到0.125～0.25的信号了，当然频域信息可能保存的比较好，但如果小波基不是对称的话，其相位信息会失真。小波变换主要也是用在高频特征提取上。

层数不是尺度，小波包分解中，N应该是层数，小波分解的尺度为a，分解层次为j。如果是连续小波分解尺度即为a。离散小波分解尺度严格意义上来说为a＝2j，在很多书上就直接将j称为尺度，因为一个j就对应者一个尺度a。其实两者是统一的。

小波基：一般从线性相位、消失矩、相似性、紧支撑等来选择。

第二节　小波分析在信号处理中的应用

一、监测信号的自相似性

直观上讲，小波分解系数表示了信号与小波之间的“相似指数”，如果相似程度越高，则相似指数越大。因此如果一个信号的不同的尺度之间相似，则小波系数在不同的尺度上也应该相似。因此可以通过小波分解检测信号的自相似性，即检测信号的分形特征。

下面通过一个简单的例子来说明小波分析在检测信号自相似性中的应用（见图2.1），待检测的信号是经过反复迭代生成的信号，因此具有自相似性。程序代码如下：

图2.1　小波分析在检测信号自相似性中的应用

从图2.1中可以看出分解后的小波系数在许多尺度上看上去都非常相似。

二、检测信号的奇异性

信号的突变点和奇异点等不规则部分通常包含重要信息。

一般信号的奇异性分为两种情况：（1）信号在某一时刻其幅值发生突变，引起信号的非连续，这种类型的突变称为第一类型的间断点；（2）信号在外观上很光滑，幅值没有发生突变，但是信号的一阶微分有突变发生且一阶微分不连续，这种类型的突变称为第二类型的间断点。

应用小波分析可以检测出信号中的突变点的位置、类型以及变化的幅度。下面介绍小波分析在信号奇异性检测中的应用。

1. 第一类型间断点的检测

在本例中，信号的不连续是由于低频特征的正弦信号在后半部分突然有高频特征的正弦信号加入。首先利用傅里叶变换分析对信号在频域进行分析，没发现检测突变点；接着利用小波分析进行分析，结果证明它能够准确地检测出信号幅值突变的位置，即高频信号加入的时间点。

程序代码如下：

源信号及分析结果如图2.2所示。由图2.2中可以看出，由于傅里叶变换不具有时间分辨力，因此无法检测信号的间断点。而在小波分析的图中，在信号的小波分解的第一层高频系数d1和第二层高频系数d2中，可以非常清楚地观察到信号的不连续点，用db1小波比用db6小波要好。

图2.2　检测信号幅值突变的位置

这个例子也表明小波分析在检测信号的奇异点时具有傅里叶变换无法比拟的优越性，利用小波分析可以精确地检测出信号的突变点。

在信号处理中，信号中通常都包含噪声，而噪声的存在增加了辨别信号不连续点的难度。一般来说，如果信号小波分解的第一层能够估计出噪声的大体位置，则信号的间断点就能够在小波分解的更深层次上表现出来。

2. 应用小波分析识别某一频率区间上的信号

在本例中，使用小波分析一个由三个不同频率的正弦信号叠加的信号，看是否能将这三个正弦信号区分开来，结果证明小波分析可以很好地识别某一频率区间的信号（见图2.3）。

图2.3　小波分析识别某一频率区间上的信号

程序代码如下：

在本例中，该信号是由周期分别为200、20、2的信号组成的，它们的采样周期均为1，为方便起见，在此分别称为低频、中频和高频的正弦信号。从图2.3中可以看出，低频、中频和高频信号分别对应于分解的近似信号a4、细节信号d4以及细节信号d1。

第三节　小波分析在图像处理中的应用

一、小波分析用于图像压缩

1. 局部压缩

基于离散余弦变换的图像压缩算法，其基本思想是在频域对信号进行分解，驱除信号点之间的相关性，并找出重要系数，滤掉次要系数，以达到压缩的效果，但该方法在处理过程中并不能提供时域的信息，在比较关心时域特性的时候显得无能为力。

但是这种应用的需求是很广泛的，比如遥感测控图像，要求在整幅图像有很高压缩比的同时，对热点部分的图像要有较高的分辨率；例如医疗图像，需要对某个局部的细节部分有很高的分辨率，单纯的频域分析的方法显然不能达到这个要求，虽然可以通过对图像进行分块分解，然后对每块作用不同阈值或掩码来达到这个要求，但分块大小相对固定，有失灵活。

在这个方面，小波分析就优越得多，由于小波分析固有的时频特性，可以在时频两个方向对系数进行处理，这样就可以对感兴趣的部分提供不同的压缩精度。

下面利用小波变化的时频局部化特性，举一个局部压缩的例子，程序代码如下：

从运行结果可以看出，小波域的系数表示的是原图像各频率段的细节信息，并且提供了一种位移相关的信息表述方式，通过对局部细节系数处理来达到局部压缩的效果。

在本例中，把图像中部的细节系数都置零，从压缩图像中可以很明显地看出只有中间部分变得模糊（比如在原图中很清晰的围巾的条纹不能分辨），而其他部分的细节信息仍然可以分辨得很清楚。

最后需要说明的是本例只是为了演示小波分析应用在图像局部压缩的方法，在实际的应用中，可能不会只做一层变换，而且作用阈值的方式可能也不会是将局部细节系数全部清除，更一般的情况是在 N 层变换中通过选择零系数比例或能量保留成分作用不同的阈值，实现分片的局部压缩。而且，作用的阈值可以是方向相关的，即在三个不同方向的细节系数上作用不同的阈值。

二维小波分析用于图像压缩是小波分析应用的一个重要方面。它的特点是压缩比高，压缩速度快，压缩后能保持图像的特征基本不变，且在传递过程中可以抗干扰。小波分析用于图像压缩具有明显的优点。

基于小波分析的图像压缩方法很多，比较成功的有小波包、小波变换零树压缩、小波变换矢量量化压缩等。

下面再给出用wdenemp函数对一个图像（文件名tire.mat）进行压缩的程序。

输出结果如下：

利用二维小波变换进行图像压缩时，小波变换将图像从空间域变换到时间域，它的作用与以前在图像压缩中所用到的离散余弦（DCT）、傅里叶变换（FFT）等的作用类似。但是要很好的进行图像的压缩，需要综合的利用多种其他技术，特别是数据的编码与解码算法等，所以利用小波分析进行图像压缩通常需要利用小波分析和许多其他相关技术共同完成。

2. 二维信号压缩中的阈值的确定

由于阈值处理只关心系数的绝对值，并不关心系数的位置，所以二维小波变换系数的阈值化方法同一维情况大同小异，为了方便用户使用小波工具箱，对某些阈值化方法提供了专门的二维处理命令，下面举例说明，在这个例子中同时采用求默认阈值的ddencmp命令和基于经验公式的wdcbm2命令对图像进行压缩，并对压缩效果进行比较。

从显示结果可见分层阈值化压缩方法同全局阈值化方法相比，在能量损失不是很大的情况下可以获得最高的压缩化，这主要是因为层数和方向相关的阈值化方法能利用更精细的细节信息进行阈值化处理。

3. 基于小波包变换的图像压缩

小波分析之所以在信号处理中有着强大的功能，是基于其分离信息的思想，分离到各个小波域的信息除了与其他小波域的关联，使得处理的时候更为灵活。全局阈值化方法作用的信息粒度太大，不够精细，所以很难同时获得高的压缩比和能量保留成分，在作用的分层阈值以后，性能明显提高，因为分层阈值更能体现信号固有的时频局部特性。

但是小波分解仍然不够灵活，分解出来的小波树只有一种模式，不能完全地体现时频局部化信息。而压缩的核心思想即是尽可能去除各小波域系数之间的信息关联，最大限度体现时频局部化的信息，因此，实际的压缩算法多采用小波包算法，而小波树的确定则是根据不同的信息论准则，以达到分解系数表达的信息密度最高。

下面通过一个例子来说明小波包分析在图像压缩中的应用，并给出性能参数以便于同基于小波分析的压缩进行比较。

在将小波包用于信号压缩的过程中，ddencmp命令返回的最优小波树标准都是阈值化标准。据这个标准确定的最优小波树可以使得压缩过程的零系数成分最高，并自动降低计算量。

最后需要说明的一点，对高频成分很多的图像，小波包的分解细节信息的特点尤其能发挥其优势。正因为这点，FBI的指纹库就是采用的基于小波包的压缩算法WSQ。

图像压缩是应用非常广泛的一类问题，所以其机器实现效率是至关重要的，在实际的应用中，如JPEG2000，一般不采用通常的mallat算法做小波分解，而是应用特定的双正交小波，利用其滤波器分布规则的特性，用移位操作来实现滤波操作。

二、小波分析用于图像去噪

噪声可以理解为妨碍人的视觉器官或系统传感器对所接收图像源进行理解或分析的各种因素。一般噪声是不可预测的随机信号，它只能用概率统计的方法去认识。噪声对图像处理十分重要，它影响图像处理的输入、采集、处理的各个环节以及输出结果的全过程。特别是图像的输入、采集的噪声是个十分关键的问题，若输入伴有较大噪声，必然影响处理全过程及输出结果。因此一个良好的图像处理系统，不论是模拟处理还是计算机处理无不把减少最前一级的噪声作为主攻目标。去噪已成为图像处理中极其重要的步骤。

对二维图像信号的去噪方法同样适用于一维信号，尤其是对于几何图像更适合。二维模型可以表述为

s(i,j)=f(i,j)+δ·e(i,j) i,j=0,1,…，m−1

其中，e是标准偏差不变的高斯白噪声。二维信号用二维小波分析的去噪步骤有3步：

（1）二维信号的小波分解。选择一个小波和小波分解的层次N，然后计算信号s到第N层的分解。

（2）对高频系数进行阈值量化。对于从1到N的每一层，选择一个阈值，并对这一层的高频系数进行软阈值量化处理。

（3）二维小波的重构。根据小波分解的第 N 层的低频系数和经过修改的从第一层到第N层的各层高频系数计算二维信号的小波重构。

在这3个步骤中，重点是如何选取阈值和阈值的量化。

下面给出一个二维信号（文件名为detfinger.mat），并利用小波分析对信号进行去噪处理。Matlab的去噪函数有ddencmp，wdencmp等，其去噪过程可以按照如下程序进行。

从输出结果可以看出，Matlab中的ddencmp和wdencmp函数可以有效地进行去噪处理。

二维信号在应用中一般表现为图像信号，二维信号在小波域中的降噪方法的基本思想与一维情况一样，在阈值选择上，可以使用统一的全局阈值，又可以分作三个方向，分别是水平方向、竖直方向和对角方向，这样就可以把在所有方向的噪声分离出来，通过作用阈值抑制其成分。

三、小波分析用于图像增强

小波分析在二维信号（图像）处理方面的优点主要体现在其时频分析特性，前面介绍了一些基于这种特性的一些应用的实例，但对二维信号小波系数的处理方法只介绍了阈值化方法一种，下面将介绍一下以前在一维信号中用到的抑制系数的方法，这种方法在图像处理领域主要应用于图像增强。

图像增强问题的基本目标是对图像进行一定的处理，使其结果比原图更适用于特定的应用领域，这里“特定”这个词非常重要，因为几乎所有的图像增强问题都是与问题背景密切相关的，脱离了问题本身的知识，图像的处理结果可能并不一定适用，比如某种方法可能非常适用于处理X射线图像，但同样的方法可能不一定也适用于处理火星探测图像。

在图像处理领域，图像增强问题主要通过时域（沿用信号处理的说法，空域可能对图像更适合）和频域处理两种方法来解决。时域方法通过直接在图像点上作用算子或掩码来解决，频域方法通过修改傅里叶变换系数来解决。这两种方法的优劣很明显，时域方法方便快速但会丢失很多点之间的相关信息，频域方法可以很详细地分离出点之间的相关，但需要做两次数量级为nlogn的傅里叶变换和逆变换的操作，计算量大得多。

小波分析是以上两种方法的权衡结果，建立在如下的认识基础上：傅里叶分析在所有点的分辨率都是原始图像的尺度，但对于问题本身的要求，可能不需要这么大的分辨率，而单纯的时域分析又显得太粗糙，小波分析的多尺度分析特性为用户提供了更灵活的处理方法。可以选择任意的分解层数，用尽可能少的计算量得到满意的结果。

小波变换将一幅图像分解为大小、位置和方向都不同的分量。在做逆变换之前可以改变小波变换域中某些系数的大小，这样就能够有选择地放大所感兴趣的分量而减小不需要的分量。下面给出一个图像增强的实例。

给定一个wmandril.mat图像信号。由于图像经二维小波分解后，图像的轮廓主要体现在低频部分，细节部分体现在高频部分，因此可以通过对低频分解系数进行增强处理，对高频分解系数进行衰减处理，从而达到图像增强的效果。具体程序清单如下：

下面将主要讨论图像增强中的钝化和锐化两种方法，钝化操作主要是提出图像中的低频成分，抑制尖锐的快速变化成分；锐化操作正好相反，将图像中尖锐的部分尽可能地提取出来，用于检测和识别等领域。将以例子说明这两种方法在Matlab中的实现，并对于基于傅里叶变换的传统频域方法同小波方法做一下比较。

1. 图像钝化

图像钝化在时域中的处理相对简单，只需要对图像作用一个平滑滤波器，使得图像中的每个点与其相邻点做平滑即可，现在来介绍一下基于傅里叶变换的频域处理方法。

下面以chess信号为例，通过两种方法对图像钝化的结果做一下比较。

从运行结果中可以看出，采用DCT在频域做滤波的方法得到钝化结果更为平滑，这是因为其分辨率最高，而小波方法得到的结果在很多地方有不连续的现象，因为对系数做放大或抑制在阈值两侧有间断，而且分解层数很低，没有完全分离出频域的信息。而且在做系数放大或抑制的时候，采用的标准是根据系数绝对值的大小，没有完全体现出其位置信息，但是在小波系数中，很容易在处理系数的过程中加入位置信息。

2. 图像锐化

与图像钝化所做的工作相反，图像锐化的任务是突出高频信息，抑制低频信息，从快速变化的成分中分离出标识系统特性或区分子系统边界的成分，以便于进一步的识别、分割等操作。

在时域（空域）中，锐化的方法不外乎是进行掩码或做差分，同钝化的道理一样，无论是掩码和差分都很难识别点之间的关联信息，下面的例子同样是在频域完成的，用传统的傅里叶分析方法（这里采用的是DCT变换）得到的频域系数。

从运行结果中可以看出，使用DCT方法进行高通滤波得到的高频结果比较纯粹，完全是原图像上的边缘信息，而在小波方法得到的结果中，不只有高频成分，还有变换非常缓慢的低频成分，这是因为两者同样在小波系数上体现为绝对值较低的部分，但这些成分的存在对进行进一步分析并无多大影响。

最后来比较一次这两个例子的时间复杂度，对DCT方法，需要做两次复杂度为O（nlogn）的DCT变换，中间系数处理部分复杂度为O（n），而对小波变换，无论是分解和重构还有系数处理的复杂度都是O（n），所以时间复杂度的优势非常明显。

四、小波分析用于图像融合

图像融合是将同一对象的两个或更多的图像合成在一幅图像中，以便它比原来的任何一幅图像更容易地为人们所理解。这一技术可应用于多频谱图像理解以及医学图像处理等领域。在这些场合，同一物体部件的图像往往是采用不同的成像机理得到的。

下面用二维小波分析将woman.mat和wbarb.mat两幅图像融合在一起。程序清单：

五、小波分析用于图像分解

回顾从一维离散小波变换到二维的扩展，二维静态小波变换采用相似的方式。对行和列分别采用高通和低通滤波器。这样分解的结果仍然是四组图像、近似图像、水平细节图像、竖直细节图像和对角图像，与离散小波变换不同的只是静态小波分解得到的四幅图像与原图像尺寸一致，道理与一维情况相同。

二维离散小波变换同样只提供了一个函数 swt2，因为它不对分解系数进行下采样，所以单层分解和多层分解的结果是一样的。不需要另外提供多层分解的功能。

下面举一个用命令行使用swt命令的例子，大家可以对比它和dwt处理结果的区别，在命令行下键入：

由于分解过程中没有改变信号的长度，所以在显示近似和细节系数时不需要重建。

同idwt2类似，Matlab对二维静态小波重建提供了iswt2命令，同idwt一维情况类似，对经过重建滤波后的信号不做上采样（因为近似系数和细节系数大小都与原信号一致）。

同一维的静态小波重建一样，将用例子说明如何将iswt2单纯用做滤波器来实现各层系数的重建，与一维的情况不同的只是为了重建第j层近似系数，需要4次用到iswt2作为重建滤波器对第j+1层的系数进行滤波，在对某一个近似系数滤波的过程中，同样需要把其他的三个系数指定为零。

为了便于比较，本例接上面的二维静态分解的例子，直接利用对noiswom的分解结果，从中重建各级系数。