小波分析及其MATLAB实现

第8章 小波分析及其MATLAB 实现

本文档节选自：

Matlab 在数学建模中的应用，卓金武等编著， 北航出版社，2011年4月出版

8.1小波分析基本理论

8.1.1 Fourier 变换的局限性

8.1.1.1变换的含义

我们把那些定义域和因变域都不是数值或常量的函数称为变换或算子，它们是定义域和值域本身为函数集的函数，如傅里叶变换（Fourier Transform ）和拉普拉斯变换（Laplace Transform ），其定义域是时间的函数，而因变域是频率的函数。简单地说，变换的基本思想仍然是映射，变换是函数的函数。

8.1.1.2 Fourier 变换的局限性

信号分析的主要目的是寻找一种简单有效的信号变换方法，使信号包含的重要特征能显示出来。在小波变换兴起之前，Fourier 级数展开和Fourier 分析是刻画函数空间、求解微分方程、进行数学计算的有效方法和有效的数学工具。从物理直观上看，一个周期性振动的量可以看成是具有简单频率的简谐振动的叠加。Fourier 级数展开则是这一物理过程的数学描述，Fourier 变化和Fourier 逆变换公式如下：

函数)()(1

R L t f ∈的连续Fourier 变换定义为t t f f t d e )()(?i -?+∞

∞

-=ωω

)(?ωf

的Fourier 逆变换定义为 ωωωd e )(?π21)(i ?+∞∞

-=t f t f 从公式上看，Fourier 变换是域变换，它把时间域和频率域联系起来，在时间域上难以

观察到的现象和规律，在频域往往能十分清楚地显示出来。频谱分析的本质就是对)(?ωf

的加工、分析和滤波等处理。Fourier 变换是平稳信号分析的最重要的工具。然而在实际运用

中，所遇到的信号大多数并不平稳，而是时变频率信号，这时人们需要知道信号在突变时刻所对应的频率成分，显然Fourier 变换的积分作用平滑了非平稳过程的突变部分，作为积分核t

ω-i e

的幅值在任何情况下均为1，因此，在频谱)(?ωf

的任一频点值是由时间过程)(t f 在整个时间域上)(∞+-∞，

上的贡献决定的；反之，时间过程)(t f 在某一时刻的状态也是由)(?ωf 在整个频率域)(∞+-∞，上贡献决定的。也就是说，Fourier 变换是整体变换，只能反

映信号的整体特征，而对信号的局部（奇异性）反映不敏感，对描述信号激烈震荡的细节无

能为力[1]。

1946年诺贝尔奖获得者D.Gabor 引入了短时傅里叶变换(Short-time Fourier Transform)，虽然短时傅里叶变换随着窗口在时间轴t 上滑动可以抠取频谱上的所有信息，但从本质上讲，短时傅里叶变换是一种单一分辨率的信号分析方法，因为它使用一个固定的短时窗函数（常用Gauss 函数），窗口大小缺乏自适应功能，窗口位置可以移动，但是形状不能改变。

下面以声波为例来说明傅里叶变换的局限性。声音的特征与振幅、时间、频率以及波形有关联。音强对应于振幅大小；音长对应于声波持续时间；音高对应于频率高低；音质不同时波形有别。Fourier 变换无法反映信号在哪一时刻有高音，在哪一时刻有低音，因此所有的音符都挤在了一起。小波变换有效地克服了Fourier 变换的这一缺点，信号变换到小波域后，小波不仅能检测到高音与低音，而且还能将高音与低音发生的位置与原始信号相对应。可以非常直观地知道什么时候发出了什么频率的音符，如图8-1所示。

图8-1 声音信号经过Fourier 变换和Wavelet 变换

8.1.2 伸缩平移和小波变换

从事石油信号处理的法国工程师J.Morlet 在研究地下岩石油层分布时，于1974年首先提出了小波变换的概念，并且J.Morlet 构建了光滑柔性的、时频域局部性能都较好的Morlet 小波函数。数学家Meyer ，Mallat ，Daubechies ，K.Chui 等人的工作为小波分析的诞生和发展奠定了基础[2]。同时计算机的发展，也为小波分析的成长提供了有利土壤，有些不能用解析式表示的小波，借助于计算机得到了广泛的应用。小波分析是纯数学、应用数学和工程技术的完美结合。小波分析的核心是小波函数的构建和多尺度分析。小波函数的主要特质是快速衰减性和震荡性，其子函数都是相互正交。这里正交的概念不是狭义上的垂直，而是指不能用任意1-N 个子小波来表征第N 个子小波（假如母小波函数经过伸缩平移得到N 个子小波）。就像平面向量永远无法表征立体向量一样。

8.1.2.1 小波的定义

小波(Wavelet)是由英文单词“Wave ”和“小”的后缀“-let ”构成，表示的是一种长度有限，平均值为0的波形。小波函数的确切定义为：设)(t ψ为一平方可积函数，即

)()(2R L t ∈ψ，若其Fourier 变换)(?ωψ

满足条件： ?∞<=R

d C ωωωψ

ψ||)(?2

（8-1）

则称)(t ψ为一个基本小波、母小波或者容许小波，我们称式（8-1）为小波函数的可容许条件。)(2

R L 表示满足

?

+∞

dt t f 2

)(的函数空间。更一般地，)(R L p 表示满足

?

+∞

p

dt t f )(的函数空间。小波)(t ψ有两个基本特点：

一是“小”，在)(2

R L 空间内，我们常选取紧支集或近似紧支集（具有时域的局限性）且具有正则性（具有频域的局限性）的实数或复数做为小波母函数，这样的小波母函数在时频域都会具有较好的局部特性，即快速衰减性。

二是正负交替的“波动性”，也即直流分量为零，即震荡性。 如下图8-2中图（a ）正弦函数不是小波函数因为它不具有快速衰减性，是无限延伸的；而图（b ）也不是小波函数因为它既不具有正负交替的震荡性也不具有衰减性；图（c ）和图（d ）为小波函数，分别是4阶复Gaussian 小波的实部和虚部部分的函数图像。

图8-2 波形示意图

8.1.2.2母小波伸缩平移

将小波母函数)(t ψ进行伸缩和平移，就可以得到子函数)(,t b a ψ： )(

1)(,a

b

t a

t b a -=

ψψ, R b a ∈,,0>a （8-2） 我们称)(,t b a ψ为依赖于参数a ,b 的小波家族子函数。上式中，a 为伸缩因子，用来控制家族子小波)(,t b a ψ图像的“体形”（“胖瘦”和“高矮”），b 为平移因子，用来控制小波子函数的中心位置。由于尺度因子a 和平移因子b 是连续变化的值，因此称)(,t b a ψ为连续小波函数基。它们是由同一母函数)(t ψ经伸缩和平移得到的一组函数序列。下图8-3为Morlet 小波经伸缩平移后得到的几何图像。

-10

010

x

y

图（a ）

-10

010

x

y

图（b ）

-5

5

x

y

图（c ）

-5

5

x

y

图（d ）

图8-3 Morlet 小波母函数经过伸缩平移在时频域的图像

关于式（8-2）两点说明：

（1）由同一母小波经过伸缩和平移得到的小波家族子函数序列均是相互正交的。 （2）)(

1)(,a

b

t a

t b a -=ψψ而不是)()(,a b t t b a -=ψψ是为了使所有子小波以及其母函数具有相等的能量。

8.1.2.3连续小波变换

定义符号>??<，为内积运算符。假设在区间),(21t t 定义两个连续时间实能量信号)(t x 和)(t y ，则此两个信号的内积表达式为：?

>=

<2

1

)()()(),(t t dt t y t x t x t y 。连续信号)(t f 在

小波基下展开称为连续小波变换（Continue Wavelet Transform ，CWT ），其表达式为：

?

->=

=

b a f dt a

b

t t f a

t t f b a WT )(

)(1)(),(),(*,ψψ 其中，

)(*

a

b t -ψ是)(a b

t -ψ的共轭函数。所谓“共轭”就是实部相同虚部相反。显然，当)(,t b a ψ为实小波函数时，有：)()(*

a

b

t a b t -=-ψψ。 8.1.3小波变换入门和多尺度分析

8.1.3.1小波变换的入门----平均和细节

小波知识艰深难懂,一般初学者难以掌握小波分析的核心思想, 孙延奎教授编著的文献

[3]

将小波本质解说的简明透彻，通俗易懂，读者可以系统参阅一下。

设{}21,x x 是一个由两个元素组成的信号，定义这两个元素的平均和细节为:

2

/)(2/)(2121x x d x x a -=+=

则可将{}d a ,作为原信号另一种表示，且原信号{}21,x x 可由{}d a ,恢复如下：

d

a x d a x -=+=21

信号{}d a ,是{}21,x x 的线性变换结果，可以看出，当1x 与2x 非常接近时，d 会很小，这时{}21,x x 可近似地用{}a 来表示，由此实现信号的压缩。重构信号为{}a a ,，误差信号为

{}{}d d a x a x

,,21

=--。特别地，当21x x =时，0=d ，

这时可以实现信号的无损重构。这表明，平均值a 可以看成是原信号的整体信息，而d 可以看成是用a 表示原始信号所丢失的细节信息。用{}a 近似地表示{}21,x x 可以实现信号压缩，而且丢失的很小细节信息对最终信号的重构不会造成影响。

对于较多元素的信号{}4321,,,x x x x ，通过如下平均和细节运算， 平均：

2

/)(2/)(431,1210,1x x a x x a +=+=

细节：

2

/)(2/)(431,1210,1x x d x x d -=-=

可得到原信号的另一种表示{}1,10,1,1,10,1,,d d a a 。若细节信号1,10,1,d d 都很小，则丢失细节信号，可得压缩后的信号为{}

1,10,1,a a 。用类似的方法对压缩后的信号{}

1,10,1,a a 进行压缩：

平均： 2/)(1,10,10,0a a a += 细节： 2/)(1,10,10,0d d d +=

如果0,01,10,1,,d d d 都非常小，则我们可以用{}

0,0a 代替原信号{}4321,,,x x x x 。事实上，

4/)(2/)2/)(2/)((2/)(432143211,10,10,0x x x x x x x x a a a +++=+++=+=

即0,0a 为整个信号所有元素的平均值，它保留了信号最基本的信息。

易知，由{}

0,00,0,d a 可以重构{}1,10,1,a a ，由{}1,10,1,a a 和{}

1,10,1,d d 可以重构

{}4321,,,x x x x ，于是，由{}1,10,10,00,0,,,d d d a 可恢复{}4321,,,x x x x ，故我们找到了原信号

序列的另一种表示：{}

1,10,10,00,0,,,d d d a 。该序列称为原始序列{}4321,,,x x x x 的小波变换。为体现小波分析中多分辨率表示的思想，我们称原信号{}4321,,,x x x x 为最高分辨率信息，而根据求平均与细节的不同层次将平均信息（也称为低频信息）与细节信息（也称为高频信息）与分辨率联系起来，具体如下：

{}4321,,,x x x x ---- 最高分辨率信息；{}1,10,1,a a ---- 次高分辨率低频信息；

{}1,10

,1,d d

---- 次高分辨率细节信息；{}0,0a ---- 最低分辨率低频信息；

{}0

,0d ---- 最低分辨率细节信息；

那么，{}4321,,,x x x x 的小波变换{}

1,10,10,00,0,,,d d d a 由信号整体平均以及两个不同分辨率级的细节信息组成，因而是{}4321,,,x x x x 的一个多分辨率表示。信号的小波变换过程可以用金字塔(Mallat)算法表示，Mallat 算法是小波变换的快速算法。以{}4213，，，-为例，它的金字塔算法如图8-3所示。

???

???23---- 最低分辨率低频信息； ?

??

???21---- 最低分辨率细节信息； {}1,2---- 次高分辨率低频信息； {}3,1----- 次高分辨率细节信息； {}4213，，，----- 最高分辨率信息；

图8-4 简易塔式算法

上述方法精准而形象地表现出了小波多尺度分析的基本思想,实际工程应用的小波变换过程与上述算法是不同的，但是上面算法很忠实地表现了小波分析的基本理念，理解它再去理解小波分析，思路就会豁然开朗。

8.1.3.2多尺度分析

首先看一个例子。用haar 小波两层子小波函数和尺度函数逼近（重构）信号

)2sin(1)(t t f π+=。haar 小波是最简单的小波，它的尺度函数在[0 1]闭区间等于1，小波

母函数在[0 0.5）半闭区间等于1，在[0.5 1]闭区间等于-1，如下图8-5所示：

图8-5 Haar 小波尺度函数和小波母函数

符号约定[4]：

ij ? 第i 层且平移j 个单位的尺度函数；ij ψ 第i 层且平移j 个单位的小波函数；

ij c 第i 层且平移j 个单位的尺度系数；ij d 第i 层且平移j 个单位的小波系数；

显然，当0=i 时，ij ?和ij ψ表示的是拉伸平移之前的尺度函数和母小波函数。为了得到各层小波系数、尺度系数以及对信号进行重构，根据连续小波变换的定义，将ij ?和ij ψ与连续信号)(t f 作内积：

1)]2sin(1[)(),(1

0000=+>==

6366.0)]2sin(1[)]2sin(1[)(),(1

2

1210

0000=+-+>==

1573.1)]2sin(1[2)(),(210

1010=+>==

dt t t t f c π?

2570.0)]2sin(1[2)(),(1

2

11111=+>==

0)]2sin(1[2)]2sin(1[2)(),(214

1410

1010=+-+>==

?dt t dt t t t f d ππψ

0)]2sin(1[2)]2sin(1[2)(),(1

4

3432

11111=+-+>==

dt t dt t t t f d ππψ

类似地有：

?????????

???=????????????=1817.01817.08183.08183.023*******c c c c c ?

?????

?????

?--=????????????=1319.01319.01319.01319.023*******d d d d d 则原信号)(t f 可以表示为：

0.2

0.4

0.6

0.8

1

x

y

Haar 小波尺度函数

x

y

)()()()()()2sin(1)(20201111101000000000t d t d t d t d t c t t f ψψψψ?π++++≈+=

)()()(232322222121t d t d t d ψψψ+++

即：)(1319.0)(6366.0)()(200000t t t t f ψψ?-+≈

)(1319.0)(1319.0)(1319.0232221t t t ψψψ

-++ 那么，用haar 小波逼近)(t f 的效果如下图8-6所示：

图8-6 haar 小波逼近效果

更一般地，假设连续信号)(t f 位于0V 空间，则)(t f 可以分解到下列空间[5]：

122110W W V W V V ⊕⊕=⊕=

1231233W W W W V W W W V n n ⊕⊕??⊕=??=⊕⊕⊕=

“⊕”表示的是“直和”之意，它有两重意思。以110W V V ⊕=为例，表示1V 和1W 都是0V 的子集，且1V 与1W 是相互正交。这里的)1(≥i V i 和)1(≥j W j 分别表示低频概貌和高频细节所在的子空间。

对连续信号)(t f 进行采样后，对这个信号进行小波多尺度分解，其实质就是把采集到的信号分成两个部分，即高频部分和低频部分，而低频部分通常包含了信号的主要信息，高频部分则是信号的细节信息，常常与噪声和扰动信息联系在一起。根据实际需要，可以继续对所得的低频部分进行再分解。信号分解的层数不是任意的，对于长度为N 的信号最多分解成N 2log 层。

8.1.4小波窗函数自适应分析

把小波)(t ψ看成一个窗函数，利用时间—频率窗来理解小波变换的时频局部化能力[6~8]。其中，||?表示取模运算，||||?表示空间内的范数，在)(2

R L 中，范数

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

00.511.5

2x

y

haar 小波函数逼近y=1+sin(2*pi*t)

2

12

2

1])([,||||?=>=

dx x f f f f

设小波母函数)(t ψ具有有限支撑即小波函数能量集中的区间，则称

dt t t t R

2

2

*

|)(|||

||1?=

ψψ (8-3)

为时窗中心。称

2

/12

2

*]|)(|)([||

||1dt t t t t R

?-=

?ψψ (8-4)

为时窗半径，称

ω

ωψωωψωd R

2

2*

|)(?|||)(?||1

?=

(8-5)

为频窗中心，称

2/122

*]|)(?|)([||)(?||1ωωψ

ωωωψωd R ?-=? (8-6)

为频窗半径。

下面计算)(,t b a ψ的时窗和频窗中心及半径。利用小波函数的基本原理，容易知道

1||)(||||)(||22,==t t b a ψψ。由式（8-3）~式（8-6），得

=-==??dt a b t a t

dt t t t R

R

b a 2

2

,*

|)(|1|)(|ψψ

=

+=+???

du u b du u au du u b au R

R

R

22

2

|)(||)(||)(|)(ψψψb at +*

ψ

(8-7)

式中，*

ψt 是)(t ψ的窗中心。

=-=??2

/12

,2

*]

|)(|)([dt t t t t R

b a ψ2/12

2

*

]|)(|1)

([dt a

b

t a b at t R ?---ψψ ψψψψψt a du u t u a

du u b at b au R

R

?=-=--+=?

?2/12

2

*2

2

/12

2

*

]|)(|)([]

|)(|)([

(8-8)

式中，ψt ?是)(t ψ的窗半径。类似地，有

22

2

2,||)(?|||)(?||)(?|||)(?||ωψωωψ

ωωψωψ

ω===??-d a a d e a a R

R

ib b a

于是有

*

2

22

,2,*

1

|)(?|||)(?||1

|)(?|||)(?||1

ψωωωψ

ωωψ

ωωψωωψ

ωa

d a a d R

R b

a b a ===

?

?

(8-9)

式中，*

ψω是)(?ωψ

的窗中心。 =-=??2/12,*

,]|)(?|)([||)(?||1ωωψ

ωωωψ

ωd R b a b a

ψψωωωψ

ωωωψ?=-?a

d a a a R 1

]|)(?|)1([||)(?||12/12

*

(8-10)

用紧支集或快速衰减为零的)(,t b a ψ乘以信号)(t f ，可以形象地称为将信号)(t f 开了一个窗。显然，窗口的性能优劣关系着对信号分析能力的强弱。函数)(,t b a ψ时间—频域窗口中心是),

(0

a

b ω，其中0ω为小波母函数)(t ψ的频率。在时—频相平面上，对于任意固定的

平移尺度b ，都可以作出)(,t b a ψ相应的时—频窗（2t ?为底，2ω?为高矩行窗口），如图8-7，其面积为：

244414 22,≥=??=????=??=??=???=S t a a

t

t a t a t S b a ωωωωωψψψψ

(8-11)

由上式可知，小波母函数)(t ψ与)(,t b a ψ函数窗口面积是相等的，是一个定值。当时窗变宽则频窗变小，这种现象称为海森堡（Heisenberg ）测不准原理。当且仅当

a

b t iat

e

a

ce t 4)(2

21)(--=πψ时等号成立。

图8-7 小波窗函数示意图

尺度倒数1/a 在一定意义上对应于)(?,ωψ

b a 的频率ω，因为主频是a

ω，而

a

ω是时—

频窗中心的纵坐标。尺度越小，对应的频率越高，时间窗越窄，小波函数越“瘦”，衰减得

越快。尺度越大，对应的频率越低，时间窗越宽，小波函数越“胖”，衰减得越慢。所以处理高频信号时，ω?变大；处理低频信号时，t ?变小。这种自适应功能很方便对信号进行各种处理。

基于小波变换的小波分析利用一个可以伸缩平移、面积不变形状可变的视窗能够聚焦到信号的任意细节进行时频域处理，既可看到信号的全貌，又可分析信号的细节，并且可以保留数据的瞬时特性，非常适合探测正常信号中夹带的瞬态反常现象并展示其成分。 8.1.5梦的凝缩与小波多尺度分解

在第7章遗传算法中，我们已经创造性将遗传算法与弗洛伊德梦的解析法联想在一起，体现了不同领域的知识之间融会贯通的发散思维，在这里我们将简单讲解“梦的凝缩”与“小波多尺度分解”之间的耦合关系。在介绍梦的凝缩作用之前，需要预先介绍一下自然界中的“逆反定律”。所谓“逆反定律”可以自定义如下：当原来的状态或平衡被打破时，物体或系统要反抗现在的力或状态，欲保持原来的力或状态的一种现象。自然界中符合这条法则的物体和系统很多，可逆反应中的勒夏特列原理：当化学平衡的条件（温度、压强、浓度等）被外界打破时，化学平衡一定会朝着减弱这种反应的方向进行移动。可逆反应的化学平衡移动只能减弱但不能抵消。可逆反应本身的属性决定平衡会发生移动但是不会恢复到平衡移动之前的状态。物理学中的楞次定律，当条形磁铁穿过螺旋形的导电线圈时，线圈产生的磁场力方向一定是阻碍（但不是阻止）条形磁铁运动的方向。光的本质是电磁破，当光在照射过程中如果遇到障碍物，就会发生衍射，衍射的方向一定是朝着阻碍光线传播的方向进行扩展。从众多例子当中，我们可以看见一个共性，当系统原有的状态或平衡被打破时，系统虽然一定会朝着外界驱动合力的方向进行运动，但是一定会受到一个阻力的作用，一个反抗当前外界驱动合力的阻力的作用。“逆反定律”在社会科学领域和自然科学领域都会有很大的应用。实际生活中，我们经常会遇见企业中下属不支持上司的做法或决定，但是下属又不得不依照上司的吩咐去做，这样下属虽然表面依照上司吩咐行事，但并不是心悦诚服，内心深处会自觉地排斥上司的这个决定。

弗洛伊德认为，一个人永远无法确定地说他已经把整个梦完完全全地解析出来。尽管所作的解释都没有什么问题，令人满意，但他仍可能再由同一个梦找出另外一个意义出来，对已经解析令人满意的梦仍可以继续分析下去，就会挖掘出更多藏在梦里的意义，也就是说梦可以一层一层往下解析。可见梦不是表面看起来那样浅显简单，而是经过“凝缩”的，如图8-8所示。

图8-8 梦的逐层分解

因此梦的内容 Dream=显意1∪显意2∪显意3∪显意4∪…… 所有看似秘密性梦的隐意部分经过逐层解析都会变成显意。但是实际当中，梦的完整解

析是十分困难的事情，原因有两点：一方面是因为随着梦的解析逐渐深入，梦的显意和梦的隐意之间相互转化和相互嵌套，涉及梦的情节越来越复杂；另一方面是梦在解析过程中，符合“逆反定律”，会受到梦的情节看似不符合逻辑、梦醒后遗忘了梦里某个重要细节等各方面的阻抗作用以阻止梦继续深入解析下去。所以要想彻头彻尾解析一个梦几乎是一个不可能事件。从理论上讲，如果梦一直分解下去，最终会挖掘到一个人道貌岸然下最原始的劣根性。从这个意义上说，梦不可能也没有必要一直解析下去。

小波多尺度分析就是把信号分解成不同频率的高频信号和低频信号，然后把低频信号再次分解成次高频信号和次低频信号，如此循环下去直到分解到可分解的层数，与用解析法解梦如出一辙。

8.2小波分析MATLAB程序设计

8.2.1 小波分析工具箱函数指令

这里只介绍常用的几种函数指令及其最常用的语法格式[9]。

（1）waveinfo函数

【语法格式】waveinfo('wname')

【实现功能】查询小波函数的基本信息。

例如查询Haar小波函数，在命令窗口输入：>> waveinfo('haar')

程序输出的结果如下：

H AARINFO Information on Haar wavelet.

Haar Wavelet

General characteristics: Compactly supported

wavelet, the oldest and the simplest wavelet.

scaling function phi = 1 on [0 1] and 0 otherwise.

wavelet function psi = 1 on [0 0.5[, = -1 on [0.5 1] and 0 otherwise.

……

其它小波函数信息查询同样使用此命令，比如查询Morlet小波函数和Mexihat函数，可以键入：>> waveinfo('morl') 和>> waveinfo('mexh')。通过此命令，可以访问到小波函数的表达式、有效支撑等基本性质。

（2）wfilters函数

【语法格式】[Lo_D,Hi_D,Lo_R,Hi_R]= wfilters (’ wname’)

[F1,F2]= wfilters (’ wname’,’type’)

【实现功能】小波滤波器。第一种格式用来计算正交小波或双正交小波’ wname’相关的4个滤波器。这4个滤波器分别为：

Lo_D——分解低通滤波器；Hi_D——分解高通滤波器；

Lo_R——重构低通滤波器；Hi_R——重构高通滤波器；

第二种格式返回以下滤波器：

如果‘type’=’d’，则返回分解滤波器Lo_D和Hi_D；

如果‘type’=’r’，则返回重构滤波器Lo_R和Hi_R；

如果‘type’=’l’，则返回低通滤波器Lo_D和Lo_R；

如果‘type’=’h’，则返回高通滤波器Hi_D和Hi_R；

wfilters函数应用程序清单：

%本例需要重点掌握wfilters函数、stem函数的用法和底层绘图技法的属性设置

clc;clear all;clf;%清屏、清内存以及清除当前图形 %’db1’, ’db2’,……’db45’,’haar ’小波都是正交小波 [Lo_D,Hi_D,Lo_R,Hi_R]= wfilters ('db45'); %stem 实现画出杆状图

subplot(221);stem(Lo_D,'color','r');xlim([0 95]); title('

分

解

低

通

滤

波

器

','fontsize',10);axis

tight;xlabel('x');ylabel('y');

subplot(222);stem(Hi_D,'color','r');xlim([0 95]); title('

分

解

高

通

滤

波

器

','fontsize',10);axis

tight;xlabel('x');ylabel('y');

subplot(223);stem(Lo_R,'color','r');xlim([0 95]); title('

重

构

低

通

滤

波

器

','fontsize',10);axis

tight;xlabel('x');ylabel('y');

subplot(224);stem(Hi_R,'color','r');xlim([0 95]); title('

重

构

高

通

滤

波

器

','fontsize',10);axis

tight;xlabel('x');ylabel('y');

程序输出的结果如图8-9所示：

图8-9 小波滤波器

（3）dwt 函数

【语法格式】[cA,cD]=dwt(X,’wname ’)或[cA,cD]=dwt(X, Lo_D,Hi_D)

【实现功能】dwt 命令使用特定的小波’wname ’或者特定的小波分解滤波器Lo_D 和Hi_D 执行单层一维小波分解。

下面来看dwt 函数的一个应用程序实例：

%本程序需要重点掌握dwt 函数用法以及图像分区画法的窍门 clc;clear all;close all;clf; a=randn(1,256); b=1.5*sin(1:256); s=a+b;

[cal,cd1]=dwt(s,'haar');

分解低通滤波器

x

y

分解高通滤波器

x

y

重构低通滤波器

x y

重构高通滤波器

x y

subplot(311);plot(s,'k-');title('原始信号','fontsize',10); axis tight;xlabel('x');ylabel('y');

subplot(323);plot(cal,'k-');title('haar 低频系数','fontsize',10); axis tight;xlabel('x');ylabel('y');

subplot(324);plot(cd1,'k-');title('haar 高频系数','fontsize',10); axis tight;xlabel('x');ylabel('y');

%计算两个相关的分解滤波器，并直接使用该滤波器计算低频和高频系数 [Lo_D,Hi_D]=wfilters('haar','d'); [cal,cd1]=dwt(s,Lo_D,Hi_D);

%进行单尺度db2离散小波变换并观察最后系数的边缘效果 [ca2,cd2]=dwt(s,'db2');%db2也是一种小波函数

subplot(325);plot(ca2,'k-');title('db2低频系数','fontsize',10); axis tight;xlabel('x');ylabel('y');

subplot(326);plot(cd2,'k-');title('db2高频系数','fontsize',10); axis tight;xlabel('x');ylabel('y');

图8-10 单层一维小波分解示意图

程序输出的结果上图8-10所示。

dwt2函数实现单尺度二维离散小波分解的功能，语法格式为：

[cA, cH , cV cD]=dwt2(X,’wname ’)或[cA, cH , cV cD]=dwt2(X, Lo_D,Hi_D) 【相关链接】idwt 函数实现单尺度一维小波重构，语法格式为： X=idwt(cA,cD, ’wname ’)或X=idwt(cA,cD, Lo_R,Hi_R) idwt2函数实现单尺度二维小波重构，语法格式为：

X=idwt2(cA, cH , cV cD, ’wname ’)或X=idwt2(cA, cH , cV cD, Lo_R,Hi_R) （4）wavedec 函数

【语法格式】[C,L]= wavedec (X,N,’wname ’);

[C,L]= wavedec (X,N, Lo_D,Hi_D)

【实现功能】使用给定的小波’wname ’或者滤波器Lo_D 和Hi_D 进行多尺度一维小波

分解。第一种格式返回信号X 在N 层的小波分解。N 必须是正整数。输出的结果包含分解

50100

150

200

250

原始信

号

x

y

haar

低频系数

x

y haar 高频系

数

x

y

db2

低频系数

x

y

db2高频系数

x

y

向量C 和相应的记录向量L 。第二种格式使用指定的低通和高通分解滤波器，返回分解结构。

下面来看wavedec 函数的应用实例：

%本程序需要重点掌握wavedec 函数和wavread 函数的用法 clear all;clc;close all;

%函数wavread 读入声音文件，扩展名.wav 不能省略。6144为采样点数量。

s=wavread('bhk_med.wav',6144);%bhk_med.wav 放在当前默认work 文件目录下 %函数subplot(2,1,1)与subplot(211)实现效果相同

figure;subplot(2,1,1);plot(4000:length(s),s(4000:length(s)),'k-'); xlim([4000 6200]);ylim([-0.08 0.08]); xlabel('信号序列');ylabel('信号值'); title('原始声音图像','fontsize',11);

[c,l]=wavedec(s(4000:length(s)),2,'db2');%使用db2小波进行两层分解 subplot(212);plot(c,'k-');xlim([0 2200]);ylim([-0.17 0.17]); title('小波分解结构','fontsize',11);

xlabel('低频系数&第二层及第一层高频的信号序列','fontsize',11); ylabel('信号值');

程序输出的结果如图8-11所示。

图8-11 多尺度一维小波分解示意图

wavedec2函数实现二维多尺度分解，语法格式为：

[C,S]= wavedec 2(X,N,’wname ’)或[C,S]= wavedec 2(X,N, Lo_D,Hi_D)

格式一使用小波’wname ’返回矩阵X 尺度为N 时的小波分解。输出是分解向量C 和相应的记录矩阵S 。N 必须是正整数。

【相关链接】waverec 函数实现一维小波重构，语法格式为： X= waverec (C,L,’wname ’)或X= waverec (C,L, Lo_R,Hi_R)

waverec 使用指定小波’wname ’或者重构滤波器Lo_R 和Hi_R 进行一维多尺度小波重构，它返回的是原信号X 。

waverec2函数实现二维小波重构，语法格式为：

X= waverec2 (C,S,’wname ’)或X= waverec2 (C,S, Lo_R,Hi_R) 用法与waverec 类似。 （

5）wrcoef 函数

信号序列

信号值

原始声音图像

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

小波分解结构

低频系数&第二层及第一层高频的信号序列

信号值

【语法格式】X =wrcoef (‘type ’,C,L,’wname ’,N)

X =wrcoef (‘type ’,C,L, Lo_R,Hi_R,N)

【实现功能】wrcoef 基于小波分解结构[C,L]，以及指定的小波’wname ’或者重构滤波器Lo_R 和Hi_R 进行一维小波系数的单支重构。第一种格式基于小波分解结构[C,L]在N 层计算重构系数向量。N 为正整数。‘type ’决定重构的系数是低频（‘type ’=’a ’）还是高频（‘type ’=’d ’）。第二种格式必须根据指定的重构滤波器进行系数重构。

wrcoef 函数应用实例：

%本程序需要重点掌握wrcoef 函数的用法和分段函数的图形画法

clc;close all;clear all;N=1000;t=1:N;sig1=sin(0.3*t);%生成正弦信号 sig2(1:500)=((1:500)-1)/500;sig2(501:N)=(1000-(501:1000))/500;%生成三

角波信号

x=sig1+sig2;[c,l]=wavedec(x,2,'db6');%进行两层小波分解

a2=wrcoef('a',c,l,'db6',2);a1=wrcoef('a',c,l,'db6',1);%重构第1~2层逼近

系数

d2=wrcoef('d',c,l,'db6',2);d1=wrcoef('d',c,l,'db6',1);%重构第1~2层细节

系数

subplot(511);plot(x,'linewidth',2);ylabel('原始信号');xlabel('信号序列

');

subplot(512);plot(a2,'linewidth',2);ylabel('a2');xlabel('信号序列'); subplot(513);plot(a1,'linewidth',2);ylabel('a1');xlabel('信号序列'); subplot(514);plot(d2,'linewidth',2);ylabel('d2');xlabel('信号序列'); subplot(515);plot(d1,'linewidth',2);ylabel('d1');xlabel('信号序列');

程序输出的结果如图8-12所示。

图8-12 一维小波系数单支重构示意图

【相关链接】wrcoef2函数实现二维小波系数的单支重构，用法与wrcoef 类似，语法格式：X =wrcoef 2(‘type ’,C,S,’wname ’,N)或X =wrcoef2 (‘type ’,C,S, Lo_R,Hi_R,N)。

8.2.2小波分析程序设计综合案例

本案例是对RGB 图像进行多尺度分解然后重构。内容涉及数字图像的程序载入、图像

原始

信号

信号序列

a 2

信号序列

a 1

信号序列

d 2

信号序列

d 1

信号序列

显示、格式转换、wavedec2函数以及wrcoef2函数等用法。小波分析程序实现的核心在于将原始数据或图片以及视频文件进行分解，对分解后的分量进行复杂的预处理，然后反变换合成。MATLAB 小波分析工具箱已经将大量实用的命令集成，调用很方便。

下面来看一段程序：

clc;close all;clear all; %从D 盘读入原始RGB 图像

X0=imread('d:\MathLogo.jpg');

X=rgb2gray(X0);%将真彩色图像转换成灰度图像因为真彩色图像是三维数据 [c,s]=wavedec2(X,4,'sym5'); %重构尺度1~4的低频部分

a1=wrcoef2('a',c,s,'sym5',1);a2=wrcoef2('a',c,s,'sym5',2); a3=wrcoef2('a',c,s,'sym5',3);a4=wrcoef2('a',c,s,'sym5',4); %绘制尺度1~4的图形,并隐藏边框和坐标轴

subplot(4,2,1); image(a1); title('1尺度低频'); box off;axis off; subplot(4,2,2); image(a2); title('2尺度低频'); box off;axis off; subplot(4,2,3); image(a3); title('3尺度低频'); box off;axis off; subplot(4,2,4); image(a4); title('4尺度低频'); box off;axis off; %重构尺度为1时的高频各分量

hd1=wrcoef2('h',c,s,'sym5',1);%重构尺度为1时的水平方向上的高频分量 vd1=wrcoef2('v',c,s,'sym5',1);%重构尺度为1时的垂直方向上的高频分量 dd1=wrcoef2('d',c,s,'sym5',1);%重构尺度为1时的对角方向上的高频分量 %重构尺度为2时的高频各分量

hd2=wrcoef2('h',c,s,'sym5',2);%重构尺度为2时的水平方向上的高频分量 vd2=wrcoef2('v',c,s,'sym5',2);%重构尺度为2时的垂直方向上的高频分量 dd2=wrcoef2('d',c,s,'sym5',2);%重构尺度为2时的对角方向上的高频分量 %重构尺度为3时的高频各分量

hd3=wrcoef2('h',c,s,'sym5',3);%重构尺度为3时的水平方向上的高频分量 vd3=wrcoef2('v',c,s,'sym5',3);%重构尺度为3时的垂直方向上的高频分量 dd3=wrcoef2('d',c,s,'sym5',3);%重构尺度为3时的对角方向上的高频分量 %重构尺度为4时的高频各分量

hd4=wrcoef2('h',c,s,'sym5',4);%重构尺度为4时的水平方向上的高频分量 vd4=wrcoef2('v',c,s,'sym5',4);%重构尺度为4时的垂直方向上的高频分量 dd4=wrcoef2('d',c,s,'sym5',4);%重构尺度为4时的对角方向上的高频分量 %画出尺度为4时高频各分量的图像

subplot(4,2,5);imshow(X0); title('原始RGB 图像');box off;axis off; subplot(4,2,6);imshow(hd4);title('4尺度水平高频');box off;axis off; subplot(4,2,7);imshow(vd4);title('4尺度垂直高频');box off;axis off; subplot(4,2,8);imshow(dd4);title('4尺度对角高频');box off;axis off;

本程序在相同的环境下运行了多次，总共使用了景物格式完全相同但是像素大小不同的两幅图像，像素（Pixel ）大小分别是3711112?和12063664?，程序输出的结果分别对应图8-13和图8-14，由此两组对比图像可以得出以下结论：

（1）图像被分解的尺度越高，清晰度越差。

（2）低频反映的是图像的基本轮廓，而高频只是细节。

（3）图像的像素越大，对多尺度分解具有相当大的阻抗作用。高分辨率图像的清晰度

随着分解尺度的加大，图像失真效果不明显。

1112?的图像中高频部分基本看不出轮廓，（4）高频与低频是相对的。在像素大小371

3664?的图像中高频部分却可以依稀看见轮廓分割线。

但在1206