用matlab小波分析的实例

1 绪论

1.1概述

小波分析是近15年来发展起来的一种新的时频分析方法。其典型应用包括齿轮变速控制，起重机的非正常噪声，自动目标所顶，物理中的间断现象等。而频域分析的着眼点在于区分突发信号和稳定信号以及定量分析其能量，典型应用包括细胞膜的识别，金属表面的探伤，金融学中快变量的检测，INTERNET的流量控制等。

从以上的信号分析的典型应用可以看出，时频分析应用非常广泛，涵盖了物理学，工程技术，生物科学，经济学等众多领域，而且在很多情况下单单分析其时域或频域的性质是不够的，比如在电力监测系统中，即要监控稳定信号的成分，又要准确定位故障信号。这就需要引入新的时频分析方法，小波分析正是由于这类需求发展起来的。

在传统的傅立叶分析中，信号完全是在频域展开的，不包含任何时频的信息，这对于某些应用来说是很恰当的，因为信号的频率的信息对其是非常重要的。但其丢弃的时域信息可能对某些应用同样非常重要，所以人们对傅立叶分析进行了推广，提出了很多能表征时域和频域信息的信号分析方法，如短时傅立叶变换，Gabor变换，时频分析，小波变换等。其中短时傅立叶变换是在傅立叶分析基础上引入时域信息的最初尝试，其基本假定在于在一定的时间窗内信号是平稳的，那么通过分割时间窗，在每个时间窗内把信号展开到频域就可以获得局部的频域信息，但是它的时域区分度只能依赖于大小不变的时间窗，对某些瞬态信号来说还是粒度太大。换言之，短时傅立叶分析只能在一个分辨率上进行。所以对很多应用来说不够精确，存在很大的缺陷。

而小波分析则克服了短时傅立叶变换在单分辨率上的缺陷，具有多分辨率分析的特点，在时域和频域都有表征信号局部信息的能力，时间窗和频率窗都可以根据信号的具体形态动态调整，在一般情况下，在低频部分（信号较平稳）可以采用较低的时间分辨率，而提高频率的分辨率，在高频情况下（频率变化不大）可以用较低的频率分辨率来换取精确的时间定位。因为这些特定，小波分析可以探测正常信号中的瞬态，并展示其频率成分，被称为数学显微镜，广泛应用于各个时频分析领域。

全文介绍了小波变换的基本理论，并介绍了一些常用的小波函数，它们的主要性质包括紧支集长度、滤波器长度、对称性、消失矩等，都做了简要的说明。在不同的应用场合，各个小波函数各有利弊。

小波分析在图像处理中有非常重要的应用，包括图像压缩，图像去噪，图像融合，图像分解，图像增强等。文中给出了详细的程序范例，用MATLAB实现了基于小波变换的图像处理。

小波分析在图像处理中有非常重要的应用，包括图像压缩，图像去噪，图像融合，图像分解，图像增强等。文中给出了详细的程序范例，用MATLAB实现了基于小波变换的图像处理。

小波分析在图像处理中有非常重要的应用，包括图像压缩，图像去噪，图像融合，图像分解，图像增强等。文中给出了详细的程序范例，用MATLAB实现了基于小波变换的图像处理。

1.2 傅立叶变换与小波变换的比较

小波分析是傅立叶分析思想方法的发展与延拓。它自产生以来，就一直与傅立叶分析

1

武汉理工大学毕业论文

密切相关。它的存在性证明，小波基的构造以及结果分析都依赖于傅立叶分析，二者是相辅相成的。两者相比较主要有以下不同：

（1）傅立叶变换的实质是把能量有限信号f（t）分解到以{ejt}为正交基的空间上去；小波变换的实质是把能量有限信号f(t)分解到Wj（j=1，2，„，J）和Vj所构成的空间上去。

（2）傅立叶变换用到基本函数只有sin(t),cos(t),exp(it)，具有唯一性；小波分析用到的函数（即小波函数）则具有不唯一性，同一个工程问题用不同的小波函数进行分析有时结果相差甚远。小波函数的选用是小波分析应用到实际中的一个难点问题（也是小波分析研究的一个热点问题），目前往往是通过经验或不断的试验（对结果进行对照分析）来选择小波函数。

（3）在频域中，傅立叶变换具有较好的局部化能力，特别是对于那些频率成分比较简单的确定性信号，傅立叶变换很容易把信号表示成各频率成分的叠加和的形式。例如，sin(1t)0.345sin(2t)4.23cos(3t)，但在时域中，傅立叶变换没有局部化能力，即无法从信号f(t)的傅立叶变换f()中看出f(t)在任一时间点附近的性态。事实上，f()d是关于频率为的谐波分量的振幅，在傅立叶展开式中，它是由f(t)的整体性态所决定的。 （4）在小波分析中，尺度a的值越大相当于傅立叶变换中的值越小。

（5）在短时傅立叶变换中，变换系数S(,)主要依赖于信号在[,]片段中的情况，时间宽度是2（因为是由窗函数g(t)唯一确定，所以2是一个定值）。在小波变换中，变换系数Wf(a,b)主要依赖于信号在[ba,ba]片段中的情况，时间宽度是，该时间宽度是随着尺度a变化而变化的，所以小波变换具有时间局部分析能力。

（6）若用信号通过滤波器来结实，小波变换与短时傅立叶变换不同之处在于：对短时傅立叶变换来说，带通滤波器的带宽f与中心频率f无关；相反，小波变换带通滤波器的带宽f则正比于中心频率f，即

2aQffC C为常数

亦即滤波器有一个恒定的相对带宽，称之为等Q结构（Q为滤波器的品质因数，且有

Q中心频率带宽）。

1.3 小波分析与多辨分析的历史

小波理论包括连续小波和二进小波变换，在映射到计算域的时候存在很多问题 ，因

为两者都存在信息冗余，在对信号采样以后，需要计算的信息量还是相当的大，尤其是连续小波变换，因为要对精度内所有的尺度和位移都做计算，所以计算量相当的大。而二进小波变换虽然在离散的尺度上进行伸缩和平移，但是小波之间没有正交性，各个分量的信息搀杂在一起，为我们的分析带来了不便。

真正使小波在应用领域得到比较大发展的是Meyer在1986年提出的一组小波，其二进制伸缩和平移构成L2(R)的标准化正交基。在此结果基础上，1988年S.Mallat在构造正交小波时提出了多分辨分析的概念，从函数分析的角度给出了正交小波的数学解释，在空间的概念上形象的说明了小波的多分辨率特性，给出了通用的构造正交小波的方法，并将之前所有的正交小波构造方法统一起来，并类似傅立叶分析中的快速傅立叶算法，给出了小波变换的快速算法——Mallat算法。这样，在计算上变得可行以后，小波变换在各个领域才发挥它独特的优势，解决了各类问题，为人们提供了更多的关于时域分析的信息。

2

武汉理工大学毕业论文

形象一点说，多分辨分析就是要构造一组函数空间，每组空间的构成都有一个统一的形式，而所有空间的闭包则逼近L2(R)。在每个空间中，所有的函数都构成该空间的标准化正交基，而所有函数空间的闭包中的函数则构成L2(R)的标准化正交基，那么，如果对信号在这类空间上进行分解，就可以得到相互正交的时频特性。而且由于空间数目是无限可数的，可以很方便地分析我们所关心的信号的某些特性。

下面我们简要介绍一下多分辨分析的数学理论。

定义：空间L2(R)中的多分辨分析是指L2(R)满足如下性质的一个空间序列VjjZ：

（1）调一致性：VjVj1，对任意jZ

jZjZ（2）渐进完全性：IVj，closeUVjL2(R) （3）伸缩完全性：f(t)Vjf(2t)Vj1

（4）平移不变性：kZ,(2j/2t)Vjj(2j/2tk)Vj

（5）Riesz基存在性：存在(t)V0，使得j(2j/2tk)|kZ构成Vj的Risez基。关于Riesz的具体说明如下：

若(t)是V0的Risez基，则存在常数A，B，且，使得：

2Ack222ck(tk)2Bck22

对所有双无限可平方和序列ck，即

ck2kZck2

成立。

满足上述个条件的函数空间集合成为一个多分辨分析，如果(t)生成一个多分辨分析，那么称(t)为一个尺度函数。

可以用数学方法证明，若(t)是V0的Riesz基，那么存在一种方法可以把(t)转化为V0的标准化正交基。这样，我们只要能找到构成多分辨分析的尺度函数，就可以构造出一组正交小波。

多分辨分析构造了一组函数空间，这组空间是相互嵌套的，即

LV2V1V0V1V2L 那么相邻的两个函数空间的差就定义了一个由小波函数构成的空间，即

VjWjVj1

并且在数学上可以证明VjWj且ViWj，ij，为了说明这些性质，我们首先来介绍一下双尺度差分方程，由于对j,VjVj1，所以对g(x)Vj，都有g(x)Vj1，也就是说可以展开成Vj1上的标准化正交基，由于(t)V0，那么(t)就可以展开成

(t)hnZn1,n(t)

这就是著名的双尺度差分方程，双尺度差分方程奠定了正交小波变换的理论基础，从数学上可证明，对于任何尺度的j,0(t)，它在j+1尺度正交基j1,n(t)上的展开系数hn是一定的，这就为我们提供了一个很好的构造多分辨分析的方法。 在频域中，双尺度差分方程的表现形式为：

ˆ(2)H()ˆ() 如果ˆ()在=0连续的话，则有

3

武汉理工大学毕业论文

ˆ()j1H(2j)ˆ(0)

说明ˆ()的性质完全由ˆ(0)决定。

2 小波分析的基本理论

2.1 从傅立叶变换到小波变换

小波分析属于时频分析的一种，传统的信号分析是建立在傅立叶变换的基础上的，由于傅立叶分析使用的是一种全局的变换，要么完全在时域，要么完全在时域，要么完全在频域，因此无法表述信号的时频局域性质，而这种性质恰恰是非平稳信号最根本和最关键的性质。为了分析和处理非平稳信号，人们对傅立叶分析进行了推广乃至根本性的，提出并发展了一系列新的信号分析理论：短时傅立叶变换、Gabor变换、时频分析、小波变换、分数阶傅立叶变换、线调频小波变换、循环统计量理论和调幅-调频信号分析等。其中，短时傅立叶变换和小波变换也是应传统的傅立叶变换不能够满足信号处理的要求而产生的。短时傅立叶变换分析的基本思想是：假定非平稳信号在分析窗函数g（t）的一个短时间间隔内是平稳（伪平稳）的，并移动分析窗函数，使f(t)g(t)在不同的有限时间宽度内是平稳信号，从而计算出各个不同时刻的功率谱。但从本质上讲，短时傅立叶变换是一种单一分辨率的信号分析方法，因为它使用一个固定的短时窗函数。因而短时傅立叶变换在信号分析上还是存在着不可逾越的缺陷。

小波变换是一种信号的时间—尺度分析方法，它具有多分辨率分析的特点，而且在时频两域都具有表征信号局部特征的能力，是一种窗口大小固定不变但其形状可改变，时间窗和频率窗都可以改变的时频局部化分析方法。即在低频部分具有较高的频率分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于探测正常信号中夹带的瞬态反常现象并展示其成分，所以被誉为分析信号的显微镜，利用连续小波变换进行动态系统故障检测与诊断具有良好的效果。

2.1.1 傅里叶变换

在信号处理中重要方法之—是傅立叶变换(FoMierTrMsroM)，它架起了时间域和频率域之间的桥梁。

对很多信号来说，傅立叶分析非常有用。因为它能给出信号令包含的各种频率成分。但是、傅立叶变换有着严重的缺点：变换之后使信号失去了时间信息，它不能告诉人们在某段时间里发生了什么变化。而很多信号都包含有人们感兴趣的非稳态(或者瞬变)持性，如漂移、趋势项、突然变化以及信号的升始或结束。这些特性是信号的最重要部分。因此傅里叶变换不适于分析处理这类信号。

虽然傅立叶变换能够将信号的时域特征和频域特征联系起来，能分别从信号的时域和频域观察，但却不能把二者有机地结合起来。这是因为信号的时域波形中不包含任何频域信息。而其傅立叶谱是信号的统计特性，从其表达式中也可以看出，它是整个时间域内的积分，没有局部化分析信号的功能，完全不具备时域信息，也就是说，对于傅立叶谱中的

4

武汉理工大学毕业论文

某一频率，不知道这个频率是在什么时候产生的。这样在信号分析中就面临一对最基本的矛盾：时域和频域的局部化矛盾。

在实际的信号处理过程中，尤其是对非平稳信号的处理中，信号在任一时刻附近的频域特征都很重要。如柴油机缸盖表面的震动信号就是由撞击或冲击产生的，是一瞬变信号，仅从时域或频域上来分析是不够的。这就促使去寻找一种新方法，能够将时域和频域结合起来描述观察信号的时频联合特征，构成信号的时频谱。这就是所谓的时频分析法，也称为时频局部化方法。

2.1.2 短时傅里叶变换

由于标准傅立叶变换只在频域里有局部分析的能力，而在时域里不存在这种能力，Dennis Gabor于1946年引入了短时傅立叶变换。短时傅立叶变换的基本思想是：把信号划分成许多小的时间间隔，用傅立叶变换分析每一个时间间隔，以便确定该时间间隔存在的频率。其表达式为

S(,)Rf(t)g()e*jtdt (2.1)

其中*表示复共轭，g(t)是有紧支集的函数，f(t)是进入分析的信号。在这个变换中，jte起着频限的作用，g(t)起着时限的作用。随着时间的变化，g(t)所确定的“时间窗”在t轴上移动，是f（t）“逐渐”进行分析。因此，g（t）往往被称之为窗口函数， S(,)大致反映了f（t）在时刻时、频率为的“信号成分”的相对含量。这样信号在窗函数上的展开就可以表示为在[,]、[,]这一区域内的状态，并把这一区域称为窗口，和分别称为窗口的时宽和频宽，表示了时频分析中的分辨率，窗宽越小则分辨率就越高。很显然，希望和都非常小，以便有更好的时频分析效果，但还森堡测不准原理指出和是互相制约的，两者不可能同时都任意小（事实上，1，且仅当2g(t)11/4e2为高斯函数时，等号成立）

2t2由此可见，短时傅立叶变换虽然在一定程度上克服了标准傅立叶不具有局部分析能力的缺陷，但它也存在着自身不可克服的缺陷，即当窗函数g(t)确定后，矩形窗口的形状就确定了，，只能改变窗口在相平面上的位置，而不能改变窗口的形状。可以说短时傅立叶变换实质上是具有单一分辨率的分析，若要改变分辨率，则必须重新选择窗函数g(t)。因此，短时傅立叶变换用来分析平稳信号犹可，但对非平稳信号，在信号波形变化剧烈的时刻，主频是高频，要求有较高的时间分辨率（即要小），而波形变化比较平缓的时刻，主频是低频，则要求有较高的频率分辨率（即要小）。而短时傅立叶变换不能兼顾两者。

2.1.3 小波变换

小波变换提出了变化的时间窗，当需要精确的低频信息时，采用长的时间窗，当需要精确的高频信息时，采用短的时间窗。

由图1．3看出，小波变换用的不是时间-频率域，而是时间-尺度域。尺度越大，采用越大的时间窗，尺度越小，采用越短的时间窗，即尺度与频率成反比。

2.2 连续小波变换

2.2.1一维连续小波变换

5

武汉理工大学毕业论文

定义：设(t)L2(R)，其傅立叶变换为ˆ()，当ˆ()满足允许条件（完全重构条件或恒等分辨条件）

CRˆ()tba2d<  (2.2)

时，我们称(t)为一个基本小波或母小波。将母函数(t)经伸缩和平移后得

a,b(t)1a() a,bR;a0 (2.3)

称其为一个小波序列。其中a为伸缩因子，b为平移因子。对于任意的函数f(t)L2(R)的连续小波变换为

Wf(a,b)f,a,ba1/2Rf(t)(tba)dt (2.4)

其重构公式（逆变换）为

f(t)1C1a2Wf(a,b)(tba)dadb (2.5)

由于基小波(t)生成的小波a,b(t)在小波变换中对被分析的信号起着观测窗的作用，所以

(t)还应该满足一般函数的约束条件

(t)dt〈 (2.6)

故ˆ()是一个连续函数。这意味着，为了满足完全重构条件式，ˆ()在原点必须等于0，即

ˆ(0)(t)dt0 (2.7)

为了使信号重构的实现在数值上是稳定的，处理完全重构条件外，还要求小波(t)的傅立叶变化满足下面的稳定性条件：

Ajˆ(2)2B (2.8)

式中0〈AB〈

从稳定性条件可以引出一个重要的概念。

~(t)，定义（对偶小波） 若小波(t)满足稳定性条件（2.8）式，则定义一个对偶小波~ˆ()由下式给出： 其傅立叶变换~ˆ()*()j(2)（2.9）

2j注意，稳定性条件（2.8）式实际上是对（2.9）式分母的约束条件，它的作用是保证对偶

小波的傅立叶变换存在的稳定性。值得指出的是，一个小波的对偶小波一般不是唯一的，然而，在实际应用中，我们又总是希望它们是唯一对应的。因此，寻找具有唯一对偶小波的合适小波也就成为小波分析中最基本的问题。 连续小波变换具有以下重要性质：

（1）线性性：一个多分量信号的小波变换等于各个分量的小波变换之和 （2）平移不变性：若f（t）的小波变换为Wf(a,b)，则f(t)的小波变换为Wf(a,b)

（3）伸缩共变性：若f（t）的小波变换为Wf(a,b)，则f（ct）的小波变换为

6

武汉理工大学毕业论文

1cWf(ca,cb),c0，

（4）自相似性：对应不同尺度参数a和不同平移参数b的连续小波变换之间是自相似的。

（5）冗余性：连续小波变换中存在信息表述的冗余度。

小波变换的冗余性事实上也是自相似性的直接反映，它主要表现在以下两个方面： （1）由连续小波变换恢复原信号的重构分式不是唯一的。也就是说，信号f（t）的小波变换与小波重构不存在一一对应关系，而傅立叶变换与傅立叶反变换是一一对应的。

（2）小波变换的核函数即小波函数a,b(t)存在许多可能的选择（例如，它们可以是非正交小波、正交小波、双正交小波，甚至允许是彼此线性相关的）。

小波变换在不同的（a，b）之间的相关性增加了分析和解释小波变换结果的困难，因此，小波变换的冗余度应尽可能减小，它是小波分析中的主要问题之一。

2.2.2 高维连续小波变换

对f(t)L2(Rn)(n1)，公式

f(t)1C1a2Wf(a,b)(tba)dadb (2.10)

存在几种扩展的可能性，一种可能性是选择小波f(t)L2(Rn)使其为球对称，其傅立叶变换也同样球对称，

ˆ()()（2.11） 并且其相容性条件变为

C(2)20(t)2dtt（2.12）

对所有的f,gL2(gn)。

daaa,bn10Wf(a,b)Wg(a,b)dbCfa,b （2.13）

这里，Wf(a,b)=〈也可以写为

〉，(t)an/2(tba)，其中aR,a0且bRn，公式（2.6）

fC1daan10RWf(a,b)na,bdb（2.14）

如果选择的小波不是球对称的，但可以用旋转进行同样的扩展与平移。例如，在二维时，可定义

2a,b,(t)a(R(11tba))（2.15）

cos这里，a0,bR，RsinC(2)2sin，相容条件变为 cosdrr020ˆ(rcos,rsin)d（2.16）

2该等式对应的重构公式为

fC1daa30dbR220Wf(a,b,)a,b,d（2.17）

对于高于二维的情况，可以给出类似的结论。

7

武汉理工大学毕业论文

2.3 离散小波变换

在实际运用中，尤其是在计算机上实现时，连续小波必须加以离散化。因此，有必要讨论连续小波a,b(t)和连续小波变换Wf(a,b)的离散化。需要强调指出的是，这一离散化都是针对连续的尺度参数a和连续平移参数b的，而不是针对时间变量t的。这一点与我们以前习惯的时间离散化不同。在连续小波中，考虑函数：

(t)aa,b1/2(tba)

这里bR，aR，且a0，是容许的，为方便起见，在离散化中，总a只

取正值，这样相容性条件就变为

Cˆ()0d (2.18)

通常，把连续小波变换中尺度参数a和平移参数b的离散公式分别取作jjaa0,bka0b0，这里jZ，扩展步长a01是固定值，为方便起见，总是假定a01（由于m可取正也可取负，所以这个假定无关紧要）。所以对应的离散小波函数j,kj,k(t)即可写作

(t)aj/20(tka0b0aj0j)a0j/2(a0tkb0)（2.19）

j而离散化小波变换系数则可表示为

Cj,kf(t)*j,k(t)dtf,j,k（2.20）

其重构公式为

f(t)CCj,kj,k(t)（2.21）

C是一个与信号无关的常数。然而，怎样选择a0和b0，才能够保证重构信号的精度呢？显然，网格点应尽可能密（即a0和b0尽可能小），因为如果网格点越稀疏，使用的小波函数j,k(t)和离散小波系数Cj,k就越少，信号重构的精确度也就会越低。

实际计算中不可能对全部尺度因子值和位移参数值计算CWTa,b值，加之实际的观测信号都是离散的，所以信号处理中都是用离散小波变换(DwT)。大多数情况下是将尺度因子和位移参数按2的幂次进行离散。最有效的计算方法是s．Mallat于1988年发展的快小波算法(又称塔式算法)。对任一信号，离散小波变换第一步运算是将信号分为低频部分〔称为近似部分)和离散部分(称为细节部分)。近似部分代表了信号的主要特征。第二步对低频部分再进行相似运算。不过这时尺度因子已经改变。依次进行到所需要的尺度。除了连续小波(CWT)、离散小波(DWT)，还有小波包（Wavelet Packet）和小波。

2.4 小波包分析

短时傅立叶变换对信号的频带划分是线性等间隔的。多分辨分析可以对信号进行有效的时频分解，但由于其尺度是按二进制变化的，所以在高频频段其频率分辨率较差，而在低频频段其时间分辨率较差，即对信号的频带进行指数等间隔划分（具有等Q结构）。小波包分析能够为信号提供一种更精细的分析方法，它将频带进行多层次划分，对多分辨率分析没有细分的高频部分进一步分解，并能够根据被分析信号的特征，自适应地选择相应频带，使之与信号频谱相匹配，从而提高了时-频分辨率，因此小波包具有更广泛的应用价值。

8

武汉理工大学毕业论文

关于小波包分析的理解，我们这里以一个三层的分解进行说明，其小波包分解树如图 AAA3 S A1 AA2 DDA3 ADA3 DA2 DDA3 AD2 AAD3 ADA3 D1 DD2 ADD3 DDD3 图1 小波包分解树

图1中，A表示低频，D表示高频，末尾的序号数表示小波分解的层树（也即尺度数）。分解具有关系：

S=AAA3+DAA3+ADA3+DDA3+AAD3+DAD3+ADD3+DDD3。

2.4.1 小波包的定义

2在多分辨分析中，L(R)Wj ,表明多分辨分析是按照不同的尺度因子j把Hilbert

jz空间L2(R)分解为所有子空间Wj(jZ)的正交和的。其中， Wj为小波函数(t)的闭包（小波子空间）。现在，我们希望几拟议部对小波子空间Wj按照二进制分式进行频率的细分，以达到提高频率分辨率的目的。

n 一种自然的做法是将尺度空间Vj和小波子空间Wj用一个新的子空间Uj统一起来表征，若令

0UjVj1UjW jZ

jn则Hilbert空间的正交分解Vj1VjWj即可用Uj的分解统一为

Un0j1UjUj jZ （2.22）

01定义子空间Uj是函数是函数Un(t)的闭包空间，而Un(t)是函数U2n(t)的闭包空间，并令

Un(t)满足下面的双尺度方程：

u2n(t)2h(k)un(2tk)kZ (2.23) u2n1(t)2g(k)un(2tk)kZ式中，g(k)(1)kh(1k)，即两系数也具有正交关系。当n=0时，以上两式直接给出

9

武汉理工大学毕业论文

u0(t)hku0(2tk)kZ （2.24） u1(t)gku0(2tk)kZ与在多分辨分析中，(t)和(t)满足双尺度方程：

(t)hk(2tk)hkkZl2kZ （2.25） 2gkkZl(t)gk(2tk)kZ相比较，u0(t)和u1(t)分别退化为尺度函数(t)和小波基函数(t)。式（2.24）是式（2.22）

的等价表示。把这种等价表示推广到nZ（非负整数）的情况，即得到（2.23）的等价表示为

nn2n1Uj1UjUj jZ；nZ （2.26） 定义（小波包） 由式（2.23）构造的序列un(t)（其中nZ）称为由基函数u0(t)=(t)确定的正交小波包。当n=0时，即为（2.24）式的情况。

由于(t)由hk唯一确定，所以又称un(t)nZ为关于序列hk的正交小波包。

2.4.2 小波包的性质

定理1 设非负整数n的二进制表示为ni1i2i1 i=0或1

则小波包un(w)的傅立叶变换由下式给出：

un(w)mi1(w/2) (2.27) ij式中

m0(w)H(w)m1(w)G(w)1212h(k)ekjkw

g(k)ekjkw定理2 设un(t)nZ是正交尺度函数(t)的正交小波包，则un(tk),un(tl)kl，即

un(t)nZ构成L2(R)的规范正交基。

2.4.3 小波包的空间分解

令un(t)nZ是关于hk的小波包族，考虑用下列方式生成子空间族。现在令n=1，2，„；

j=1，2，„，并对（2.22）式作迭代分解，则有

WjUU2j11jU4j22j1U5j23j15j2UU,UU6j2U7j2

因此，我们很容易得到小波子空间Wj的各种分解如下：

WjUWjU2j14j2UU3j15j2U6j2U7j2

10

武汉理工大学毕业论文

„

WjU2jkkUj21jkkj„U2jkk1Uj12jkk1

„

WjUWl20U210„U021

ll2mm，m=0，1，„，2l-1；l=1，2，„。子空间序列U2j空间分解的子空间序列可写作Uj1j1m的标准正交基为2(j1)/2u2m(2jltk):kZ。容易看出，当l=0和m=0时，子空间序列U2j1l简化为U1j=Wj，相应的正交基简化为2j/2u1(2jtk)2j/2(2jtk)，它恰好是标准正交

小波族j,k(t)。

若n是一个倍频程细划的参数，即令n=2l+m，则我们有小波包的简略记号

l/2lj,k,n(t)2j/2n(2jtk)，其中，n(t)2u2m(2t)。我们把j,k,n(t)称为既有尺度指标

lj、位置指标k和频率指标n的小波包。将它与前面的小波j,k,(t)作一比较知，小波只有

离散尺度j和离散平移k两个参数，而小波包除了这两个离散参数外，还增加了一个频率参数n=2l+m。正是这个频率新参数的作用，使得小波包克服了小波时间分辨率高时频率分辨率低的缺陷，于是，参数n表示n(t)2l/2u2m(2lt)函数的零交叉数目， 也就是其波形

l的震荡次数。

定义（小波库） 由n(t)生成的函数族j,k,n(t)（其中nZ；j，kZ）称为由尺度函

数(t)构造的小波库。

推论1.1 对于每个j=0，1，2，„

232L(R)Wj=„W1W0U0U0„ （2.28）

jZ这时，族

{uj,k,un(tk)|j=„，-1，0；n=2，3，„且kZ} （2.29）

是L2(R)的一个正交基。

随着尺度j的增大，相应正交小波基函数的空间分辨率越高，而其频率分辨率越低，这正是正交小波基的一大缺陷。而小波包却具有将随j增大而变宽的频谱窗口进一步分割变细的优良性质，从而克服了正交小波变换的不足。

小波包可以对Wj进一步分解，从而提高频率分辨率，是一种比多分辨分析更加精细的分解方法，具有更好的时频特性。

2.4.4 小波包算法

n(t)Uj，则gj可表示为 下面给出小波包的分解算法和重构算法。设gnjngj(t)ndlj,nlun(2tl) (2.30)

j小波包分解算法 由dlj1,n求dlj,2n与dlj,2n1

dlj,2nakkk2ldkj1,ndkj1,ndlj,2n1b （2.31）

k2l小波包重构算法 由{dlj,2n}与dlj,2n1求dlj1,n

11

武汉理工大学毕业论文

3 几种常用的小波

1)Haar小波

A.Haar于1990年提出一种正交函数系,定义如下：

0x1/211 1/2x1 0其它H这是一种最简单的正交小波，即

(t)(xn)dx0 n1,2,„

2)Daubechies（dbN）小波系

该小波是Daubechies从两尺度方程系数hk出发设计出来的离散正交小波。一般简写为dbN，N是小波的阶数。小波和尺度函数吁中的支撑区为2N-1。的消失矩为N。除N＝1外(Haar小波)，dbN不具对称性〔即非线性相位〕；dbN没有显式表达式(除N＝1外)。但hk的传递函数的模的平方有显式表达式。假设P(y)项式的系数，则有

m0()2N1Ck0N1kkyk，其中，CkN1k为二

(cos22)P(sinN22)

其中 m0()122N1hk0ekik

3)Biorthogonal（biorNr.Nd）小波系

Biorthogonal函数系的主要特征体现在具有线性相位性，它主要应用在信号与图像的重构中。通常的用法是采用一个函数进行分解，用另外一个小波函数进行重构。Biorthogonal函数系通常表示为biorNr.Nd的形式： Nr=1 Nd=1，3，5 Nr=2 Nd=2，4，6，8 Nr=3 Nd=1，3，5，7，9 Nr=4 Nd=4 Nr=5 Nd=5 Nr=6 Nd=8

其中，r表示重构，d表示分解。 4)Coiflet（coifN）小波系

coiflet函数也是由Daubechies构造的一个小波函数，它具有coifN（N=1，2，3，4，5）这一系列，coiflet具有比dbN更好的对称性。从支撑长度的角度看，coifN具有和db3N

12

武汉理工大学毕业论文

及sym3N相同的支撑长度；从消失矩的数目来看，coifN具有和db2N及sym2N相同的消失矩数目。

5)SymletsA（symN）小波系

Symlets函数系是由Daubechies提出的近似对称的小波函数，它是对db函数的一种改进。Symlets函数系通常表示为symN（N=2，3，„，8）的形式。 6)Morlet（morl）小波

Morlet函数定义为(x)Cex7)Mexican Hat（mexh）小波

Mexican Hat函数为

2/2cos5x，它的尺度函数不存在，且不具有正交性。

(x)231/4(1x)e2x/22

它是Gauss函数的二阶导数，因为它像墨西哥帽的截面，所以有时称这个函数为墨西哥帽函数。墨西哥帽函数在时间域与频率域都有很好的局部化，并且满足

(x)dx0

由于它的尺度函数不存在，所以不具有正交性。 8)Meyer函数

Meyer小波函数和尺度函数都是在频率域中进行定义的，是具有紧支撑的正交小波。

2431/2j/2(2)esin((1))3322ˆ()(2)1/2ej/2cos((31)) 48 3322280[,]33其中，(a)为构造Meyer小波的辅助函数，且有

(2)33241/2ˆ()(2)cos((1)) 22330431/22

13

武汉理工大学毕业论文

4 小波变换在图像处理中的应用

4.1 小波分析用于图像压缩

4.1.1 基于小波变换的图像局部压缩

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

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

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

下面我们利用小波变化的时频局部化特性，举一个局部压缩的例子，大家可以通过这个例子看出小波变换在应用这类问题上的优越性。

load wbarb

％使用sym4小波对信号进行一层小波分解 [ca1,ch1,cv1,cd1]=dwt2(X,'sym4'); codca1=wcodemat(ca1,192); codch1=wcodemat(ch1,192); codcv1=wcodemat(cv1,192); codcd1=wcodemat(cd1,192);

％将四个系数图像组合为一个图像

codx=[codca1,codch1,codcv1,codcd1] ％复制原图像的小波系数 rca1=ca1; rch1=ch1; rcv1=cv1; rcd1=cd1;

％将三个细节系数的中部置零

rch1(33:97,33:97)=zeros(65,65); rcv1(33:97,33:97)=zeros(65,65); rcd1(33:97,33:97)=zeros(65,65);

14

武汉理工大学毕业论文

codrca1=wcodemat(rca1,192); codrch1=wcodemat(rch1,192); codrcv1=wcodemat(rcv1,192); codrcd1=wcodemat(rcd1,192);

％将处理后的系数图像组合为一个图像

codrx=[codrca1,codrch1,codrcv1,codrcd1] ％重建处理后的系数

rx=idwt2(rca1,rch1,rcv1,rcd1,'sym4');

subplot(221);image(wcodemat(X,192)),colormap(map);title('原始图像'); subplot(222);image(codx),colormap(map);title('一层分解后各层系数图像'); subplot(223);image(wcodemat(rx,192)),colormap(map);title('压缩图像'); subplot(224);image(codrx),colormap(map);title('处理后各层系数图像'); ％求压缩信号的能量成分 per=norm(rx)/norm(X) per =1.0000

％求压缩信号与原信号的标准差 err=norm(rx-X) err =

586.4979 运行结果如图

图2 利用小波变换的局部压缩图像

从图1可以看出，小波域的系数表示的是原图像各频率段的细节信息，并且给我们提

15

武汉理工大学毕业论文

供了一种位移相关的信息表述方式，我们可以通过对局部细节系数处理来达到局部压缩的效果。

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

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

4.1.2 小波变换用于图像压缩的一般方法

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

4.1.2.1 利用二维小波分析进行图像压缩

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

下面给出一个图像信号（即一个二维信号，文件名为wbarb.mat），利用二维小波分析对图像进行压缩。一个图像作小波分解后，可得到一系列不同分辨率的子图像，不同分辨率的子图像对应的频率是不相同的。高分辨率（即高频）子图像上大部分点的数值都接近于0，越是高频这种现象越明显。对一个图像来说，表现一个图像最主要的部分是低频部分，所以一个最简单的压缩方法是利用小波分解，去掉图像的高频部分而只保留低频部分。图像压缩可按如下程序进行处理。

%装入图像 load wbarb; %显示图像

subplot(221);image(X);colormap(map) title('原始图像'); axis square

disp('压缩前图像X的大小：'); whos('X')

%对图像用bior3.7小波进行2层小波分解 [c,s]=wavedec2(X,2,'bior3.7');

%提取小波分解结构中第一层低频系数和高频系数 ca1=appcoef2(c,s,'bior3.7',1); ch1=detcoef2('h',c,s,1); cv1=detcoef2('v',c,s,1); cd1=detcoef2('d',c,s,1); %分别对各频率成分进行重构

a1=wrcoef2('a',c,s,'bior3.7',1); h1=wrcoef2('h',c,s,'bior3.7',1); v1=wrcoef2('v',c,s,'bior3.7',1);

16

武汉理工大学毕业论文

d1=wrcoef2('d',c,s,'bior3.7',1); c1=[a1,h1;v1,d1];

%显示分解后各频率成分的信息 subplot(222);image(c1); axis square

title('分解后低频和高频信息'); %下面进行图像压缩处理

%保留小波分解第一层低频信息，进行图像的压缩 %第一层的低频信息即为ca1，显示第一层的低频信息 %首先对第一层信息进行量化编码 ca1=appcoef2(c,s,'bior3.7',1); ca1=wcodemat(ca1,440,'mat',0); %改变图像的高度 ca1=0.5*ca1;

subplot(223);image(ca1);colormap(map); axis square

title('第一次压缩');

disp('第一次压缩图像的大小为：'); whos('ca1')

%保留小波分解第二层低频信息，进行图像的压缩，此时压缩比更大 %第二层的低频信息即为ca2，显示第二层的低频信息 ca2=appcoef2(c,s,'bior3.7',2); %首先对第二层信息进行量化编码 ca2=wcodemat(ca2,440,'mat',0); %改变图像的高度 ca2=0.25*ca2;

subplot(224);image(ca2);colormap(map); axis square

title('第二次压缩');

disp('第二次压缩图像的大小为：'); whos('ca2') 输出结果如下所示： 压缩前图像X的大小：

Name Size Bytes Class

X 256x256 524288 double array

Grand total is 65536 elements using 524288 bytes 第一次压缩图像的大小为：

Name Size Bytes Class

ca1 135x135 145800 double array

Grand total is 18225 elements using 145800 bytes 第二次压缩图像的大小为：

Name Size Bytes Class

ca2 75x75 45000 double array Grand total is 5625 elements using 45000 bytes

图像对比如图所示。可以看出，第一次压缩提取的是原始图像中小波分解第一层的低

17

武汉理工大学毕业论文

频信息，此时压缩效果较好，压缩比较小（约为1/3）：第二次压缩是提取第一层分解低频部分的低频部分（即小波分解第二层的低频部分），其压缩比较大（约为1/12），压缩效果在视觉上也基本过的去。这是一种最简单的压缩方法，只保留原始图像中低频信息，不经过其他处理即可获得较好的压缩效果。在上面的例子中，我们还可以只提取小波分解第3、4、„层的低频信息。从理论上说，我们可以获得任意压缩比的压缩图像。

原始图像 分解后低频和高频信息

第一次压缩图像 第二次压缩图像

图3 利用二维小波分析进行图像压缩

下面给出一个图像信号（即一个二维信号，文件名为wbarb.mat），利用二维小波分析对图像进行压缩。一个图像作小波分解后，可得到一系列不同分辨率的子图像，不同分辨率的子图像对应的频率是不同的。高分辨率（即高频）子图像上大部分点的数值都接近于0，越是高频这种现象越明显。对一个图像来说，表现一个图像最主要的部分是低频部分，所以一个最简单的压缩方法是利用小波分解，去掉图像的高频部分而只保留低频部分。图像压缩可按如下程序进行处理。

% 调入图像

X=imread('lena.bmp’); % 归一化图像

X=double(sig)/255;

18

武汉理工大学毕业论文

% 显示图像 image(X);

colormap(map)

%对图像用bior3.7小波进行2层小波分解 [c,s]=wavedec(X,2,'bior3.7'); %设置小波系数阈制值 % thr=20;

% 提取小波分解结构中第一层的低频系数和高频系数 ca1=appcoef2(c,s,'bior3.7',1); ch1=detcoef2('h',c,s,1); cv1=detcoef2('v',c,s,1); cd1=detcoef2('d',c,s,1); %分别对各频率成分进行重构

a1=wrcoef2('a',c,s,'bior3.7',1); h1=wrcoef2('h',c,s,'bior3.7',1); v1=wrcoef2('v',c,s,'bior3.7',1); d1=wrcoef2('d',c,s,'bior3.7',1); c1=[a1,h1,v1,d1];

%下面进行图像压缩处理

%保留小波分解第一层低频信息，进行图像的压缩 %第一层的低频信息即为ca1，显示第一层的低频信息 %首先对第一层信息进行量化编码 ca1=appcoef2(c,s,'bior3.7',1); ca1=wcodemat(ca1,440,'mat',0); %改变图像的高度 ca1=0.5*ca1

%显示第一次压缩图像 image(ca1)

colormap(map);

%保留小波分解第二层低频信息，进行图像的压缩，此时压缩比更大 %第二层的低频信息即为ca2，显示第二层的低频信息 ca2=appcoef2(c,s,'bior3.7',2); %首先对第二层信息进行量化编码 ca2=wcodemat(ca2,440,'mat',0); %改变图像的高度 ca2=0.25*ca2

%显示第二次压缩图像 image(ca2); colormap(map);

下面再给出用wdenemp函数对一个图像（文件名tire.mat）进行压缩的程序。 %装入一个二维信号 load tire; %显示图像

subplot(221);image(X);colormap(map) title('原始图像');

19

武汉理工大学毕业论文

axis square

%下面进行图像压缩

%对图像用db3小波进行2层小波分解 [c,s]=wavedec2(X,2,'db3');

%使用wavedec2函数来实现图像的压缩

[thr,sorh,keepapp]=ddencmp('cmp','wv',X); %输入参数中选择了全局阈值选项‘gbl’，用来对所有高频系数进行相同的阈值量化处理

[Xcomp,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'db3',2,thr,sorh,keepapp); %将压缩后的图像与原始图像相比较，并显示出来 subplot(222);image(Xcomp);colormap(map) title('压缩图像'); axis square

disp('小波分解系数中置0的系数个数百分比：'); perf0

disp('压缩后图像剩余能量百分比：'); perfl2

输出结果如下所示：

小波分解系数中置0的系数个数百分比： perf0 =49.1935

压缩后图像剩余能量百分比： perfl2 =99.9928 图像对比如图所示：

原始图像 压缩图像

图4 利用二维小波分析对图像进行压缩

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

4.1.2.2 二维信号压缩中的阈值的确定与作用命令

由于阈值处理只关心系数的绝对值，并不关心系数的位置，所以二维小波变换系数的阈值化方法同一维情况大同小异，为了方便用户使用小波工具箱对某些阈值化方法提供了专门的二维处理命令

下面我们通过一个例子来说明二维信号的小波压缩的一般方法，在这个例子中我们同时采

20

武汉理工大学毕业论文

用求缺省阈值的ddencmp命令和基于经验公式的wdcbm2命令对图像进行压缩，并对压缩效果进行比较。

load detfingr;

％求得颜色映射表的长度，以便后面的转换 nbc=size(map,1);

％用缺省方式求出图像的全局阈值

[thr,sorh,keepapp]=ddencmp('cmp','wv',X); thr thr = 3.5000

％对图像作用全局阈值

[xd,cxd,lxd,perf0,perfl2]=wdencmp('gbl',X,'bior3.5',3,thr,sorh,kaapapp); ％用bior.3.5小波对图像进行三层分解 [c,s]=wavedec2(X,3,'bior3.5');

％指定Birge-Massart策略中的经验系数 alpha=1.5;m=2.7*prod(s(1,:)); ％根据各层小波系数确定分层阈值

[thr1,nkeep1]=wdcbm2(c,s,alpha,m); ％对原图像作用分层阈值

[xd1,cxd1,sxd1,perf01,perfl21]=wdencmp('lvd',c,s,'bior3.5',3,thr1,'s'); thr1 thr1 =

14.7026 68.4907 93.8430 14.7026 68.4907 93.8430 14.7026 68.4907 93.8430 ％将颜色映射表转换为灰度映射表 colormap(pink(nbc));

subplot(221);image(wcodemat(X,nbc));title('原始图像'); subplot(222);image(wcodemat(xd,nbc));

title('全局阈值化压缩图像');xlabel(['能量成分',num2str(perfl2),'%','零系数成分',num2str(perf0),'%']);

subplot(223);image(wcodemat(xd1,nbc));

title('分层阈值化压缩图像');xlabel(['能量成分',num2str(perfl21),'%','零系数成分',num2str(perf01),'%']); 显示结果如图所示

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

21

武汉理工大学毕业论文

图5 detfingr图像的全局阈值化压缩和分层阈值化压缩

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

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

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

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

％读入信号 load julia

％求颜色索引表长度 nbc=size(map,1);

％得到信号的阈值，保留层数，小波树优化标准 [thr,sorh,keepapp,crit]=ddencmp('cmp','wp',X) thr =

3.0000 sorh = h

keepapp = 1

22

武汉理工大学毕业论文

crit = threshold

％通过以上得到的参数对信号进行压缩

[xd,treed,perf0,perfl2]=wpdencmp(X,sorh,4,'sym4',crit,thr*2,keepapp); ％更改索引表为pink索引表 colormap(pink(nbc));

subplot(121);image(wcodemat(X,nbc));title('原始图像');

subplot(122);image(wcodemat(xd,nbc));title('全局阈值化压缩图像');

xlabel(['能量成分',num2str(perfl2),'%','零系数成分',num2str(perf0),'%']); plot(treed);

得到的压缩结果如图所示

图6 基于小波包分析的图像压缩

压缩过程中使用的最优小波数如图所示

23

武汉理工大学毕业论文

图7 最优小波树

这两个命令是Matlab小波工具箱提供的自动获取阈值和自动使用小波包压缩的命令，后者将分解阈值化和重建综合起来。在将小波包用于信号压缩的过程中，ddencmp命令返回的最优小波树标准都是阈值化标准。根据这个标准确定的最优小波树可以使得压缩过程的零系数成分最高，并且自动降低计算量。

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

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

4.2 小波分析用于图像去噪

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

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

24

武汉理工大学毕业论文

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等，其去噪过程可以按照如下程序进行。

程序清单： %装入图像 load tire

%下面进行早声的产生 init=37180252; rand('seed',init);

Xnoise=X+18*(rand(size(X))); %显示原始图像及它的含噪声的图像 colormap(map);

subplot(2,2,1);image(wcodemat(X,192)); title('原始图像X') axis square

subplot(2,2,2);image(wcodemat(X,192)); title('含噪声的图像Xnoise'); axis square

%用sym5小波对图像信号进行二层的小波分解 [c,s]=wavedec2(X,2,'sym5'); %下面进行图像的去噪处理

%使用ddencmp函数来计算去噪的默认阈值和熵标准 %使用wdencmp函数来实现图像的压缩

[thr,sorh,keepapp]=ddencmp('den','wv',Xnoise);

[Xdenoise,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'sym5',2,thr,sorh,keepapp);

%显示去噪后的图像

subplot(223);image(Xdenoise); title('去噪后的图像'); axis square

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

25

武汉理工大学毕业论文

原始图像 含噪声的图像 去噪后的图像

图8 去噪例一

再给定一个有较大白噪声的delmontl.mat图像。由于图像所含的噪声主要是白噪声，而且主要集中在图像的高频部分，所以我们可以通过全部滤掉图像中的高频部分实现图像的去噪。具体去噪过程可按照如下程序进行。程序清单：

%下面装入原始图像，X中含有被装载的图像 load wmandril; %画出原始图像

subplot(221);image(X);colormap(map); title('原始图像'); axis square %产生含噪图像

init=2055615866;randn('seed',init) x=X+38*randn(size(X)); %画出含噪图像

subplot(222);image(x);colormap(map); title('含噪声图像'); axis square;

%下面进行图像的去噪处理

%用小波函数sym4对x进行2层小波分解 [c,s]=wavedec2(x,2,'sym4');

%提取小波分解中第一层的低频图像，即实现了低通滤波去噪 a1=wrcoef2('a',c,s,'sym4'); %画出去噪后的图像

subplot(223);image(a1); title('第一次去噪图像'); axis square;

%提取小波分解中第二层的低频图像，即实现了低通滤波去噪 %相当于把第一层的低频图像经过再一次的低频滤波处理 a2=wrcoef2('a',c,s,'sym4',2); %画出去噪后的图像

subplot(224);image(a2);title('第二次去噪图像'); axis square; 输出结果如图：

26

武汉理工大学毕业论文

原始图像 含噪声图像

第一次去噪图像 第二次去噪图像

图9 去噪例二

从上面的输出结果可以看出，第一次去噪已经滤去了大部分的高频噪声，但从去噪图像与原始图像相比可以看书，第一次去噪后的图像中还是含有不少的高频噪声；第二次去噪是在第一次去噪的基础上，再次滤去其中的高频噪声。从去噪的结果可以看出，它具有较好的去噪效果。

下面再给出定一个喊有较少噪声的facets.mat图像。由于原始图像中只喊有较少的高频噪声，如果按照上一个例子把高频噪声全部滤掉的方法将损坏图像中固有的高频有用信号。因此这幅图像适合采用小波分解系数阈值量化方法进行去噪处理。程序清单：

%下面装入原始图像，X中含有被装载的图像 load facets; %画出原始图像

subplot(221);image(X);colormap(map); title('原始图像'); axis square

%产生含噪声图像

init=2055615866;randn('seed',init) x=X+10*randn(size(X)); %画出含噪声图像

subplot(222);image(X);colormap(map); title('含噪声图像'); axis square

%下面进行图像的去噪处理

%用小波画数coif3对x进行2层小波分解 [c,s]=wavedec2(x,2,'coif3');

%提取小波分解中第一层的低频图像，即实现了低通滤波去噪 %设置尺度向量n n=[1,2]

设置阈值向量p p=[10.12,23.28];

27

武汉理工大学毕业论文

%对三个方向高频系数进行阈值处理 nc=wthcoef2('h',c,s,n,p,'s'); nc=wthcoef2('v',c,s,n,p,'s'); nc=wthcoef2('d',c,s,n,p,'s');

%对新的小波分解结构[nc，s]进行重构 xx=waverec2(nc,s,'coif3'); %画出重构后图像的波形

subplot(223);image(X);colormap(map); title('去噪后的图像'); axis square 输出结果如图

原始图像 含噪声图像 去噪后图像

图10 去噪例三

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

4.3 小波分析用于图像增强

4.3.1 图像增强问题描述

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

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

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

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

28

武汉理工大学毕业论文

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

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

load wmandril %画出原始图像

subplot(221);image(X);colormap(map); title('原始图像'); axis square

%下面进行图像的增强处理

%用小波函数sym4对X进行2层小波分解 [c,s]=wavedec2(X,2,'sym4'); sizec=size(c);

%对分解系数进行处理以突出轮廓部分，弱化细节部分 for i=1:sizec(2) if(c(i)>350) c(i)=2*c(i); else

c(i)=0.5*c(i); end end

%下面对处理后的系数进行重构 xx=waverec2(c,s,'sym4'); %画出重构后的图像

subplot(222);image(xx); colormap(map);

title('增强图像'); axis square 输出结果如图所示：

原始图像 增强图像

图11 小波分析用于图像增强

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

29

武汉理工大学毕业论文

4.3.2 图像钝化

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

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

％分别保存用DCT方法和小波方法的变换系数 blur1=X; blur2=X;

％对原图像做二维离散余弦变换 ff1=dct2(X);

％对变换结果在频域做BUTTERWORTH滤波 for i=1:256 for j=1:256

ff1(i,j)=ff1(i,j)/(1+((i*j+j*j)/8192)^2); end end

％重建变换后的图像 blur1=idct2(ff1);

％对图像做2层的二维小波分解 [c,l]=wavedec2(X,2,'db3'); csize=size(c);

％对低频系数进行放大处理，并抑制高频系数 for i=1:csize(2); if(c(i)>300) c(i)=c(i)*2; else

c(i)=c(i)/2; end end

％通过处理后的小波系数重建图像 blur2=waverec2(c,l,'db3'); ％显示三幅图像

subplot(221);image(wcodemat(X,192));colormap(gray(256));title('原始图像','fontsize',18);

subplot(223);image(wcodemat(blur1,192));colormap(gray(256));title('采用DCT方法钝化图像','fontsize',18);

subplot(224);image(wcodemat(blur2,192));colormap(gray(256));title('采用小波方法钝化图像','fontsize',18);

在命令行中运行该程序后，得到如图所示的显示结果。

30

武汉理工大学毕业论文

（a）

（b）

图12 图像钝化

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

4.3.3 图像锐化

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

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

％读入chess信号 load chess;

％分别保存用DCT方法和小波方法的变换系数 blur1=X; blur2=X;

31

武汉理工大学毕业论文

％对原图像做二维离散余弦变换 ff1=dct2(X);

％对变换结果在频域做BUTTERWORTH滤波 for i=1:256 for j=1:256

ff1(i,j)=ff1(i,j)/(1+(32768/(i*i+j*j))^2); end end

％重建变换后的图像 blur1=idct2(ff1);

％对图像做2层的二维小波分解 [c,l]=wavedec2(X,2,'db3'); csize=size(c);

％对低频系数进行放大处理，并抑制高频系数 for i=1:csize(2); if(abs(c(i))<300) c(i)=c(i)*2; else

c(i)=c(i)/2; end end

％通过处理后的小波系数重建图像 blur2=waverec2(c,l,'db3'); ％显示三幅图像

subplot(221);image(wcodemat(X,192));colormap(gray(256));title('原始图像','fontsize',18);

subplot(223);image(wcodemat(blur1,192));colormap(gray(256));title('采用DCT方法锐化图像','fontsize',18);

subplot(224);image(wcodemat(blur2,192));colormap(gray(256));title('采用小波方法锐化图像','fontsize',18); 得到的结果如图所示

（a）

32

武汉理工大学毕业论文

（b）

图13 图像锐化

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

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

4.4 小波分析用于图像融合

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

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

load woman; X1=X;map1=map; %画出原始图像

subplot(221);image(X1);colormap(map1); title('woman'); axis square load wbarb; X2=X;map2=map; for i=1:256 for j=1:256

if (X2(i,j)>100)

X2(i,j)=1.2*X2(i,j); else

X2(i,j)=0.5*X2(i,j); end end end

subplot(222);image(X2);colormap(map2);

33

武汉理工大学毕业论文

title('wbarb'); axis square

%用小波函数sym4对X1进行2层小波分解 [c1,s1]=wavedec2(X1,2,'sym4');

%对分解系数进行处理以突出轮廓部分，弱化细节部分 sizec1=size(c1); for i=1:sizec1(2) c1(i)=1.2*c1(i); end

%用小波函数sym4对X2进行2层小波分解 [c2,s2]=wavedec2(X2,2,'sym4'); %下面进行小波变换域的图像融合 c=c1+c2;

%减小图像亮度 c=0.5*c;

%对融合的系数进行重构

xx=waverec2(c,s1,'sym4'); %画出融合后的图像

subplot(223);image(xx); title('融合图像'); axis square 输出结果如图：

（a）

34

武汉理工大学毕业论文

（b）

图14 小波分析用于图像融合

4.5 小波分析用于图像分解

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

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

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

load noiswom

[swa,swh,swv,swd]=swt2(X,3,'db1');

％使用db1小波对noiswom图像进行三层静态小波分解 whos

可以看出，swt2所小波分解同样不改变信号的长度，原来的96×96的图像做了三层分解以后，分解系数是12个96×96的图像。

colormap(map) kp=0;

for i=1:3

subplot(3,4,kp+1),image(wcodemat(swa(:,:,i),192)); title(['Approx,cfs,level',num2str(i)])

％显示第i层近似系数图像，以192字节为单位编码 subplot(3,4,kp+2),image(wcodemat(swh(:,:,i),192)); title(['Horiz.Det.cfs level',num2str(i)])

subplot(3,4,kp+3),image(wcodemat(swv(:,:,i),192)); title(['Vert.Det.cfs level',num2str(i)])

subplot(3,4,kp+4),image(wcodemat(swd(:,:,i),192)); title(['Diag.Det.cfs level',num2str(i)]) kp=kp+4; end

显示的结果如图所示，由于分解过程中没有改变信号的长度，所以在显示近似和细节

35

武汉理工大学毕业论文

系数时不需要重建。

图15 小波分析用于图像分解

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

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

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

load noiswom

[swa,swh,swv,swd]=swt2(X,3,'db1');

％使用db1小波对noiswom图像进行三层小波分解 mzero=zeros(size(swd)); A=mzero;

A(:,:,3)=iswt2(swa,mzero,mzero,mzero,'db1');

％使用iswt2的滤波器功能，重建第3层的近似系数，为了避免iswt的合成运算，注意在重建过程中，应保证其他各项系数为零。

H=mzero;V=mzero;D=mzero; for i=1;3

swcfs=mzero;swcfs(:,:,i)=swh(:,:,i);

H(:,:,i)=iswt2(mzero,swcfs,mzero,mzero,'db1');

36

武汉理工大学毕业论文

swcfs=mzero;swcfs(:,:,i)=swv(:,:,i);

V(:,:,i)=iswt2(mzero,mzero,swcfs,mzero,'db1'); swcfs=mzero;swcfs(:,:,i)=swh(:,:,i);

H(:,:,i)=iswt2(mzero,mzero,mzero,swcfs,'db1'); end

％分别重建1~3级的各个细节系数，同样在重建某一吸收的时候，要令其他系数为0 A(:,:,2)=A(:,:,3)+H(:,:,3)+V(:,:,3)+D(:,:,3); A(:,:,1)=A(:,:,2)+H(:,:,2)+V(:,:,2)+D(:,:,2); ％使用递推的方法建立地1层和第2层近似系数 colormap(map) kp=0;

for i=1:3

subplot(3,4,kp+1),image(wcodemat(A(:,:,i),192)); title(['第',num2str(i),'层近似系数图像'])

subplot(3,4,kp+2),image(wcodemat(H(:,:,i),192)); title(['第',num2str(i),'层水平细节系数图像']) subplot(3,4,kp+3),image(wcodemat(V(:,:,i),192)); title(['第',num2str(i),'层竖直细节系数图像']) subplot(3,4,kp+4),image(wcodemat(D(:,:,i),192)); title(['第',num2str(i),'层对角细节系数图像']) kp=kp+4; end

％画出通过手工方法重建的各级小波系数图像 err=norm(A(:,:,2)-swa(:,:,2))

％求出用这种算法重建的第2层近似系数和分解系数之间的误差 err=2.9242e+004 显示的结果如图

37

武汉理工大学毕业论文

图16 在db1小波下各级静态小波重建系数

5 全文总结

本论文主要结合小波变换的基本概念和基本原理，详细讨论小波在图像处理领域的应用，并结合MATLAB程序设计语言来说明其应用。

第一个重点是了解小波变换和小波分析的理论和方法。主要以小波的发展过程为线索，从传统的傅立叶分析入手，通过对比傅立叶分析、短时傅立叶分析和小波分析的优缺点，指出了小波分析在分析数学领域的地位及其必然性。小波的数学理论和发在科学技术界引起一场轩然大波。在数学家眼中，它被认为是调和分析，即现代傅立分析这一重要科学半个世纪以来的工作之结晶；在其它科学技术领域，特别是在信号分析、图像分析、量子物理和非线性科学等部门，它被当作近年来在工具和方法上的重大突破。

第二个重点是研究小波分析和小波变换在图像处理中的应用。小波分析的应用是与小波分析的理论研究紧密地结合在一起的。现在，它已经在科技信息领域取得了令人瞩目的成就。电子信息技术是六大高新技术中的一个重要领域，图像和信号处理又是电子信息技术领域的重要方面。现今，信号处理已经成为当代科学技术工作的重要组成部分。现在，对性质随时间稳定不变的信号，处理的理想工具仍然是傅立叶分析。但在实际应用中，绝大多数信号是非稳定的，小波分析正是适用于非稳定信号的处理工具。图像处理是针对性很强的技术，根据不同应用、不同要求需要采用不同的处理方法。采用的方法是综合各学科较先进的成果而成的，如数学、物理学、心理学、信号分析学、计算机学、和系统工程等。计算机图像处理主要采用两大类方法：一类是空域中的处理，即在图像空间中对图像进行各种处理；另一类是把空间与图像经过变换，如傅立叶变换，变到频率域，在频率域中进行各种处理，然后在变回到图像的空间域，形成处理后的图像。图像处理是“信息处理”的一个方面，这一观点现在已经为人所熟知。它可以进一步细分为多个研究方向：图

38

武汉理工大学毕业论文

片处理、图像处理、模式识别、景物分析、图像理解、光学处理等等。小波分析用在图像处理方面，主要是用来进行图像压缩、图像去噪、图像增强（包括图像钝化和图像锐化）、图像融合、图像分解。

论文中还将详细介绍用MATLAB来实现这些图像处理技术。MATLAB的语法规则比FORTRAN语言和C语言等高级语言更简单，更重要的是它贴近人思维方式的编程特点，使用MATLAB编写程序有如在便笺上列公式和求解。这一部分讨论的问题主要是小波变换从函数空间域映射到计算域的问题，也就是如何用计算机实现小波变换，因为小波的应用范围非常广泛，而每个领域都有自己特定的处理方式，所以了解这些技术对读者选用工具箱中合适的函数来解决自己的问题很有帮助。

6 参考文献

[1]胡昌华，张军波，夏军，张伟．基于MATLAB的系统分析与设计——小波分析．西安：西安电子科技大学出版社，2000

[2]董长虹，高志，余啸海．Matlab小波分析工具箱原理与应用．国防工业出版社，2004 [3]张兆礼，张春晖，梅晓丹．现代图像处理技术及Matlab实现．人民邮电出版社．2001 [4]刘贵忠，邸双亮．小波分析及其应用．西安电子科技大学出版社，1992 [5]奉前清，杨宗凯.实用小波分析.西安：西安电子科技大学出版社，2000 [6]阮秋琦.数字图像处理学.电子工业出版社，2001

[7]肖自美.图像信息理论与压缩编码技术.中山大学出版社，2000 [8]容观澳.计算机图像处理.清华大学出版社，2000

[9]姚东，王爱民，冯峰，王朝阳.MATLAB命令大全.人民邮电出版社，2000 [10]陈贺新.非线形滤波器与数字图像处理.国防工业出版社，1997

[11]Cohen A.Wavelets and multiscale signal processing. Chapman and Hall,1995 [12]Chui,C.K.An introduction to wavelets.Academic Press,1992 [13]Daubechies I. Ten lectures on wavelets.SLAM,1992

[14]Antoniadis ,Pham D T. Wavelets regression for random or irregular design.Comp.Stat.and Data Analysis,1998

39