现代信号处理教程-胡广书(清华)
第12章 双正交小波及小波包
我们在上一章给出了正交小波的构造方法。正交小波有许多好的性质,如
j,k(t),φj,k(t)=δ(k?k'),j,k(t),ψj,k(t)=δ(k?k'),j,k(t),ψj,k(t)=0 ,此
'
'
'
外,尺度函数和小波函数都是紧支撑的,有着高的消失矩等等。Daubechies给出的正交小波的构造方法可以方便的构造出所需要的小波(如DBN,SymN,CoifN)。但是,正交小波也有不足之处,即φ(t)和ψ(t)都不是对称的,尽管SymN和CoifN接近于对称,但毕竟不是真正的对称,因此,这在实际的信号处理中将不可避免地带来相位失真。φ(t)和ψ(t)的不对称性来自所使用的共轭正交滤波器组H0(z)和H1(z)的不对称性。我们已在7.8节讨论了具有线性相位的双正交滤波器组的基本概念,给出了可准确重建的双正交滤波器组的设计方法。本章,我们把这些内容引入到小波分析,给出适合小波变换的双正交滤波器组准确重建的条件,给出双正交条件下的多分辨率分析及双正交小波的构造方法,最后简要讨论小波包的基本概念
12.1 双正交滤波器组
现在,我们结合小波变换的需要来研究双正交滤波器组的内在关系及实现准确重建
的条件。所谓“小波变换的需要”是指在用H0(z)对a0(z)分解时需要将H0(z)和H1(z)的或0(n)=h0(?n),1(n)=h1(?n),系数作时间上的翻转,即用的是H0(z?1)及H1(z?1),见(10.6.1)式及图10.6.2。将图10.6.2的正变换和图10.6.3的反变换结合起来,我们可得到如图12.1.1所示的一级分解和重建的类似于两通道滤波器组的信号流图。注意,图中用于
?0(z)和H?1(z),它们分别是重建的滤波器不再是图10.6.3中的H0(z)和H1(z),而是H
H0(z)和H1(z)的对偶滤波器。有关“对偶”的概念见1.6节,在下面的讨论中将涉及对
偶滤波器的作用。
现在我们来分析该图中各信号之间的关系及实现PR的条件。由第七章关于两通道滤波器组的理论,我们有
- 352 -
图12.1.1 双正交滤波器组
)
a1(n)=a0(n)?0(2n)
=∑a0(k)h0(k?2n)=a0(k),h0(k?2n) (12.1.1a)
k
d1(n)=a0(n)?1(2n) =
∑a(k)h(k?2n)=
1
k
a0(k),h1(k?2n) (12.1.1b)
?(n)+d'(n)?h?(n) ?0(n)=a1'(n)?ha011
?(n?2l)+d(l)h=∑a1(l)h∑1?1(n?2l) 0
l
l
(12.1.2)
将(12.1.1)式代入(12.1.2)式,有
?(n?2l) ?0(n)=∑a0(k),h0(k?2l)a0
l
+
∑
l
?(n?2l) 0(k),h1(k?2l)1
(12.1.3)
(12.1.1)式是用一组向量{h0(k?2n),h1(k?2n),n,k∈Z}对a0(n)作分析,(12.1.3)式是用
?(n?2l),h?(n?2l),n,l∈Z对a(n)作综合。(12.1.3)式还可表为 一组对偶向量h010
?(n?2l)(k) ?0(n)=∑h0(k?2l),ha00
l
{}
+
∑
l
?(n?2l)(k) h1(k?2l),h10
- 353 -
(12.1.4)
显然,如果
?(n?2l)=δ(n?k) h0(k?2l),h0?(n?2l)=δ(n?k) h1(k?2l),h1
则
(12.1.5a) (12.1.5b)
?0(n)=2a0(n) a
从而实现了准确重建。(12.1.5)式的含意是,在图12.1.1中,同一条支路上的两个滤波器
?(n)或h(n),h?(n)的偶序号位移之间是正交的。但是该式没有涉及上下支路两个h0(n),h011
滤波器之间的关系。我们更关心的是这些滤波器系数的移位可否构成小波分析中的基函数。下面的两个定理清楚地回答了该问题。
定理12.1
对图12.1.1所示的两通道滤波器组,对任意的输入信号a0(n),其准确
重建的充要条件是:
*?0(ω)+H1*(ω+π)H?1(ω)=0 H0(ω+π)H
(12.1.6a) (12.1.6b)
及
?0(ω)+H1(ω)H?1(ω)=2 H0(ω)H
*
*
证明:仿照(7.1.5)式的导出,有
?(z)=1H(z?1)H?0(z)+H1(z?1)H?1(z)A0(z) A00
2
+
[]
1?(z)A(?z) ?(z)+H(?z?1)HH0(?z?1)H01102
[]
(12.1.7)
?(z)分别是a(n)和a?0(n)的z变换,A0(?z)是混迭分量。因此,为消除式中A0(z)、A00
混迭失真,应有
?(z)+H(?z?1)H?(z)=0 H0(?z?1)H011
(12.1.8a)
为保证系统的准确重建,应有
?(z)+H(z?1)H?(z)=2cz?k H0(z?1)H011
(12.1.8b)
式中c和k均为常数。令c=1,k=0,(12.1.8)式对应的频率表示是:
*?0(ω)+H1*(ω+π)H?1(ω)=0 H0(ω+π)H
- 354 -
*?0(ω)+H1*(ω)H?1(ω)=2 H0(ω)H
于是定理得证。
对比图7.1.1的两通道滤波器组,其对应的PR条件是(见(7.1.5)式):
H0(?z)G0(z)+H1(?z)G1(z)=0
H0(z)G0(z)+H1(z)G1(z)=2
(12.1.9a) (12.1.9b)
将(12.1.9)和(12.1.8)式相比较可以看出,在双正交滤波器组的情况下,我们分别用
?0(z)、H?(z)代替了G0(z)和G1(z),并在分析滤波器组中,用H0(z?1)、H1(z?1)分H1
别代替了H0(z)和H1(z)。其实,(12.1.8)式导出的原理和(12.1.9)式是完全一样的。 由(12.1.6a)式,有
??(ω)??2?H1(ω)??H?H0(ω)0
?????=?? ?
?H0(ω+π)H1(ω+π)??H1(ω)??0?
(12.1.10)
可求出
??(ω)??H?H1(ω+π)?20????=??H(ω+π)?
det()ωH0???H1(ω)?
(12.1.11)
式中
detH(ω)=H0(ω)H1(ω+π)?H1(ω)H0(ω+π)
(12.1.12)
?(z)和H?1(z)是稳定的,detH(ω)在ω=?π~π的范围内显然,为了保证对偶滤波器H0?0(z)和H?(z)是FIR的,detH(ω)应取纯延迟的形式。 应该非零。为了保证H1
仿照(7.2.16)式对G0(z)和G1(z)的定义,我们可给出在双正交条件下对偶滤波器和分析滤波器之间的关系: 或
??(ω+π) H1(ω)=e?j(2l+1)ωH0?(ω)=e?j(2l+1)ωH?(ω+π) H10
(12.1.13a) (12.1.13b)
?0(?z?1) H1(z)=z?(2l+1)H
- 355 -
(12.1.14a)
?1(z)=z?(2l+1)H0(?z?1) H
(12.1.14b)
假定l=0,它们对应的时域关系是
?(1?n) h1(n)=(?1)n+1h0?(n)=(?1)n+1h(1?n) h10
(12.1.15a) (12.1.15b)
注意,上述时域、频域关系均是在图12.1.1中的交叉方向上给出的,它正好反映了双正交滤波器组的特点。
将(12.1.13)式代入(12.1.6)式,我们可得到如下的关系:
??(ω)+H?(ω+π)H?(ω+π)=2 H0(ω)H000
(12.1.16a)
或 及
??1(ω)+H1?(ω+π)H?1(ω+π)=2 H1(ω)H
(12.1.16b)
或
??1(ω)+H0?(ω+π)H?1(ω+π)=0 H0(ω)H
(12.1.17a) (12.1.17b)
??0(ω)+H1?(ω+π)H?0(ω+π)=0 H1(ω)H
至此,我们给出了在双正交滤波器组中的若干基本关系,即
(1) 去除混迭条件:(12.1.6a)式; (2) PR条件
:(12.1.6b)式;
(3) 保证PR条件和滤波器均为FIR的情况下,四个滤波器在时域和频域的关系:
(12.1.13)式~(12.1.17)式。
回顾在共轭正交滤波器组的情况下,我们经常用到的功率互补关系,即
H0(ω)+H0(ω+π)=2,
或
H0(ω)H0(ω)+H0(ω+π)H0(ω+π)=2
?
?
22
(12.1.18)
?0(z)=H0(z),则(12.1.16a)式即变成(12.1.18)式,也即双正交滤波器组变成了显然,若H
正交滤波器组。
有了以上讨论的基础,我们可给出在小波分析中要用到的“基”的概念。
- 356 -
?0(z)和H?1(z)满足准定理12.2 如果图12.1.1中的四个滤波器H0(z),H1(z),H
确重建条件,且它们的傅里叶变换均是有界的,则
?(n?2l),h?(n?2l)},l∈Z 和 {h(n?2l),h(n?2l)},l∈Z {h0101
是L2(R)中的双正交Riesz基。
?及h?的偶序号项移位是双正交的,我们需要证明如下三个证明:为证明h0、h1、h01
关系成立: 及
?(k),h(k?2n)=δ(n) h00?(k),h(k?2n)=δ(n) h11
(12.1.19a) (12.1.19b)
?(k),h(k?2n)=h?(k),h(k?2n)=0 h0110
(12.1.19c)
由(12.1.16a)式,有
1??0(ω)+H0?(ω+π)H?0(ω+π)=1 H0(ω)H
2
该式对应的时域关系是
∞
[]
??(2n)=h00
k=?∞
?(k)h(k?2n)=δ(n)
∑h
(12.1.20)
于是(12.1.19a)式得证。同理,由(12.1.16b)式可证明(12.1.19b)式,而(12.1.17)式对应的时域关系即是(12.1.19c)式。这样,(12.1.19)式给出了三组正交关系。
?,h,h?的偶序号位移能够构成L2(R)中的双正交Riesz基,它们还需满足若h0,h011
如下的条件:
∞
11
?ω∈[?π,π],有≤∑?(ω+2kπ)≤
Bk=?∞A
2
(12.1.21)
?,h?(ω)是θ的傅里叶变换,此处θ代表h,h此即(10.2.11)式。式中A>0,B>0,θ001
?。由本定理所给的条件,即它们的傅里叶变换都是有界的,所以(12.1.21)式满足,因或h1
- 357 -
?,h及h?的偶序号移位构成L(R)中的双正交Riesz基。于是定理得证。 此h0,h011
2
我们之所以说这些序列为“双正交”基,是因为在图12.1.21中的滤波器组中,上下支
?正交,h和其对偶h?正交;同时,上下支路交叉正交,路各自是正交的,即h0和其对偶h011?。?,h正交于h即h0正交于 h在双正交滤波器中,我们并没有强调H0(z)和H1(z)0注意,11
之间的正交关系,而这一正交关系是共轭正交滤波器组中的基本关系。由此读者可搞清正交和双正交的区别。总之,在小波的多分辨率分析中,使用正交滤波器组时,分解滤波器和重建滤波器是相同的,而在双正交小波分析中,分析滤波器是H0和H1,而综合滤波器
?0和H?。 是它们的对偶,即H1
此外,(12.1.19a)和 (12.1.19b)的双正交关系与本章开头所给出的(12.1.5)式的关系是一致的,只不过(12.1.19)式更简洁。
12.2 双正交小波
上一节我们讨论了双正交滤波器的基本概念、PR条件及各滤波器时域、频域的关系。本节,我们将把双正交滤波器组的概念引入双正交小波变换,给出类似第十章的多分辨率分析。
由(9.8.18)和 (9.8.19)式,信号x(t)的离散小波变换是:
WTx(j,k)=∫x(t)ψj,k(t)dt=x(t),ψj,k(t)j,k∈Z
(12.2.1)
令dj(k)=WTx(j,k),则dj(k)称为小波系数,也即x(t)的DWT。我们可由dj(k)重建
x(t)。由(9.8.20)式,有
x(t)=∑
∞
j=0k=?∞
∑d
∞
j
?j,k(t)=∑(k)ψ
∞
j=0k=?∞
∑
∞
?j,k(t) x(t),ψj,k(t)(12.2.2)
?j,k(t)是ψj,k(t)的对偶小波。由以上两式可以看出,小波ψj,k(t)用于信号的分析,式中ψ
?j,k(t)用于信号的综合。在正交小波的情况下,ψ?j,k(t)=ψj,k(t)。 对偶小波ψ
我们在第十章详细讨论了离散小波变换的多分辨率分析,引出了尺度函数φ(t),证明
- 358 -
了在L2(R)中存在正交基φj,k(t)和ψj,k(t),给出了φj,k(t)、ψj,k(t)和正交滤波器组的关系,即二尺度差分方程和(10.4.7)和(10.4.8)式的频域关系。在双正交滤波器组的情况下,分
?)和两个小波函数?,H?1)将产生两个尺度函数(φ,φ解滤波器(H0,H1)和重建滤波器(H0
?和ψ?0,H?)。其中φ和ψ对应信号的分解,而φ?对应信号的重建。它们和H0,H(ψ,ψ1?1相应的时域和频域的关系是: 及H
及
φ(t)=2∑h0(n)φ(2t?n)
n=?∞∞
∞
(12.2.3a) (12.2.3b) (12.2.4a) (12.2.4b)
?(t)=2hφ∑?0(n)φ?(2t?n)
n=?∞∞
ψ(t)=2∑h1(n)φ(2t?n)
n=?∞∞
?(n)φ?(2t?n) ?(t)=∑hψ1
n=?∞
Φ(2ω)=
1
H0(ω)Φ(ω) 2
(12.2.5a)
?(2ω)=1H?0(ω)Φ?(ω) Φ
Ψ(2ω)=
1
H1(ω)Φ(ω) 2
(12.2.5b) (12.2.6a) (12.2.6b)
?(2ω)=1H?1(ω)Φ?(ω) Ψ
2
定理10.3给出了在正交滤波器组情况下H0(ω)和H1(ω)的关系,即(10.5.1)式。对应双正交滤波器组,这一关系变成:
*?0(ω)+H0*(ω+π)H?0(ω+π)=2 H0(ω)H
(12.2.7)
此即(12.1.6a)式。由(12.1.13)式,令l=0,则分解和重建滤波器之间有如下关系:
?0(?z?1),或H1(ω)=e?jωH?0?(ω+π) H1(z)=z?1H
- 359 -
(12.2.8a)
?1(z)=z?1H0(?z?1),或H?1(ω)=e?jωH0?(ω+π) H
(12.2.8b)
?(t)都是低通的,ψ(t)和ψ?(t)都是带通的。对同正交小波时一样,我们要求φ(t)和φ
?(z)是低通的,H(z)和H?(z)是高通的,即 应的,要求H0(z)和H011
?0(ω)ω=π=0 H0(ω)=H
?1(ω)ω=0=0 H1(ω)=H
(12.2.9a)
(12.2.9b) (12.2.10a) (12.2.10b)
∫φ(t)dt=∫φ?(t)dt=1 ∫ψ(t)dt=∫ψ?(t)dt=0
由(12.1.16)式,有
?0(ω)ω=0=2,及 H0(ω)=H?0(ω)ω=0=2 H0(ω)H
?1(ω)ω=π=2,及 H1(ω)=H?1(ω)ω=π= H1(ω)H
(12.2.11a) (12.2.11b)
类似(10.4.14)式,可由(12.2.5)式导出
H0(ω2j)
Φ(ω)=∏j=1
∞?(ω2j)H?(ω)=∏0 Φ
2j=1
∞
(12.2.12a)
(12.2.12b)
类似(10.4.15)式,可由(12.2.6)式导出
H1(ω/2)∞H0(ω2j)
Ψ(ω)=∏22j=2
?(ω2j)?(ω/2)∞HH0?(ω)=1 Ψ∏22j=2
(12.2.13a)
(12.2.13b)
由上面的讨论可知,在“双正交”的情况下,我们在第七章及第十章所讨论的滤波器
?;H,H?1;φ,φ和ψ,ψ?。组及两尺度差分方程各增加了一套“对偶”,即H0,H01
上面各节给出了它们所应满足的时域及频域关系。下面的定理将给出双正交小波基的存在
- 360 -
性。
?(ω),使得 定理12.3[42,5,8] 假定存在两个恒正的三角多项式p(ω)和p
H0()p()+H0(+π)p(+π)=2p(ω)
2222
ω
2
ωω
2
ω
(12.2.14a) (12.2.14b)
ω?(ωp(ω)+H?(ω+π)p??(ω) H(+π)=2p00
2
2
2
2
?0(ω)在?π~π内非零,则 H0(ω)、H
22
2
22
并假定
?(t)属于L(R),且满足双正交关系 1. 由(12.2.12)式定义的φ(t)和φ
?(t?n)=δ(n) (t),φ
(12.2.15)
?j,k(t)是L2(R)中的双正交Riesz基,即 2. 两个小波函数序列ψj,k(t)和ψ
?j,k(t)=δ(j?j')δ(k?k') j,k(t),ψ
'
'
(12.2.16)
该定理的证明见文献。有了L2(R)中的双正交基,我们可对x(t)作如下的分解:
x(t)=∑
2
∞
j=0k=?∞
∑
∞
∞
?j,k(t) x(t),ψj,k(t)?j,k(t)j,k(t) x(t),ψ
(12.2.17)
=∑
∞
j=0k=?∞
∑
?j,k(t)是L(R)中的Riesz基,则必然存在常数A>0,B>0,使得 既然ψj,k(t),ψ
Ax(t)≤∑x(t),ψj,k(t)
j,k
2
2
≤Bx(t)
2
2
(12.2.18a)
12
?j,k(t)x(t)≤∑x(t),ψ
Bj,k
≤
12
x(t) A
(12.2.18b)
由上面的讨论可知,在双正交的情况下,我们并不要求{ψj,k}和{ψj,k'}之间是正交的,
?}和{ψ?j',k}之间是正交的,仅要求也不要求{φj,k}和{ψj,k}之间,以及其对偶函数{φj,k
?'}之间以及{ψ}和{ψ?j',k'}之间是正交的,也即(12.2.15)和(12.2.16)式。正交{φj,k}和{φj,kj,k
- 361 -
性的放宽是使H0(z)及H1(z)具有线性相位,从而使φ(t)和ψ(t)更具有对称性,从而减小了相位失真。
在第十章的多分辨率分析中,我们假定
Vj=close{φj,k,j,k∈Z} Wj=close{ψj,k,j,k∈Z} Vj=Vj+1⊕Wj+1,Wj⊥Vj
(12.2.19a) (12.2.19b)
并有 (12.2.19c)
?将产生两个空间。除了(12.2.19a)和(12.2.19b)在双正交情况下,尺度函数φj,k及其对偶φj,k
式的关系外,还有
?,j,k∈Z} ?j=close{φ Vj,k?j=close{ψ?j,k,j,k∈Z} W
(12.2.20a) (12.2.20b)
及
?j的嵌套关系是 Vj和V
V?1?V0?V1???Vj?Vj+1?
(12.2.21a)
???V??V????V?j?V?j+? V1011
(12.2.21b)
?j,W和W?j之间有如下关系: 此时,Wj不再是Vj的正交补空间,但Vj,Vj
?j, V?j⊥Wj Vj⊥W
(12.2.22a) (12.2.22b)
?j?1=V?j⊕W?j Vj?1=Vj⊕Wj,V
由1.7节关于正交基的性质,有
k=?∞∞
?(ω+2kπ)=0
∑Φ(ω+2kπ)Φ
??
∞
(12.2.23a)
k=?∞
?(ω+2kπ)=0
∑Ψ(ω+2kπ)Ψ
- 362 -
(12.2.23b)
双正交小波下的快速算法和正交基小波下的快速算法基本相同,区别是在重建时使用的是
?(z)和H?1(z)。具体的分解方程和重建方程是: 对偶滤波器H0
aj(n)=aj?1(n)?0(2n)= dj(n)=aj?1(n)?1(2n)=
k=?∞
∑a
∞
j?1
(k)h0(k?2n) (12.2.24a)
k=?∞
∑a
∞
j?1
(k)h1(k?2n) (12.2.24b)
'?(n)+d'(n)?h?(n) aj?1(n)=aj(n)?h0j1
∞
∞
=
'
'
k=?∞
∑a
j
?(n?2k)+(k)h0
k=?∞
∑d
j
?(n?2k) (12.2.25) (k)h1
式中aj(n),dj(n)分别是aj(n),dj(n)作二插值得到的序列,见图12.1.1。
12.3 双正交小波的构造
?(t)的构造,而它们又都源于分解滤波?(t),φ(t)及φ双正交小波的构造包括ψ(t),ψ
?0(z)和H?(z)。(12.1.14)式给出了H1(z)、器H0(z)、H1(z)及用于重建的对偶滤波器H1?0(z)及H(z)的关系,因此,双正交小波构造的核心问题是H(z)和H?(z)的?(z)和HH0001
构造,这和正交小波的构造过程是一样的。如同第十一章关于正交小波的讨论,在具体给出双正交小波的构造方法之前,先讨论一下有关支撑范围、消失矩等有关的有关问题。
1. 支撑范围
?(n)都是FIR滤波器,?(t),?(t)如果h0(n)和h由(12.2.3)和(12.2.4)式,φψ(t)及ψφ(t),0?(n)的支撑范围分别是N≤n≤N,N?≤n≤N?,则将都具有有限支撑。若h0(n)和h01212
?(t)的支撑范围分别是[N,N]和N?,N?,而小波函数ψ(t)和ψ?(t)的支撑范围φ(t)和φ1212
分别是
[]
- 363 -
?+1N?N?+1??N??N+1N??N+1??N1?N221121
,,2?和??
2222????
??N?)/2 它们的长度都是(N2?N1+N21
2. 消失矩
?(ω)在ω=π处零点的数目。由定理?(t)消失矩的数目取决于H0(ω)和Hψ(t)和ψ0?(ω)在ω=π则ψ(t)有p阶消失矩。同理,若H11.1,若H0(ω)在ω=π处有p阶零点,0?(z)时,应尽量让它们?(t)有p?阶零点,则ψ?阶消失矩。因此,在构造H0(z)和H处有p0
在ω=π处有高阶的重零点。
3. 规则性
此处不再详细讨论,其一般结论是:
a) 由(12.2.4a)式,φ(t)和ψ(t)有着相同的规则性;
b) φ(t)和ψ(t)的规则性随着H0(ω)在ω=π处零点数的增加而增加; c)
?(t)和ψ?0(ω)在ω=π处零点数的增加而增加; ?(t)的规则性也是随着Hφ
?0(ω)在ω=π处有不同的零点数,?(t)的规则性则ψ(t)和ψd) 如果H0(ω)和H
也不相同。
4. 对称性
之所以使用双正交小波,其目的是使H0(z),H1(z)及其对偶滤波器具有线性相位,同时也使φ(t)和ψ(t)都具有对称性。除Haar小波外,在正交小波的情况下,上述对称性
?(n)具有奇数长且以n=0为对称,则φ(t)和φ?(t)是以是不可能实现的。如果h0(n),h0?(n)具有偶?(t)是相对位移位中心为对称的。如果h0(n),ht=0为对称的,而ψ(t)和ψ0
?(t)是以t=1/2为中心作对称,而ψ(t)和数长且以n=1/2为中心作对称,则φ(t)和φ?(t)以其位移中心作反对称。 ψ
- 364 -
显然,若h0(n),则图12.1.1中的H0(z),h1(n)是对称的,H1(z)都可改记为H0(z)和H1(z),也即在对aj(n)作分解时无需再将h0(n)和h1(n)翻转。
5.
?1?1
?(z)的构造 H0(z)及H0
?(z)具有线性相位,因此,它们的频率响应可表为: 由于要求H0(z)及H0
H0(ω)=ejkωH0(ω) ?0(ω)=ejk?ωH?0(ω) H
(12.3.1a) (12.3.1b)
?(n)为这是和Daubechies正交小波的一个主要区别。在实际工作中,我们总选取h0(n)和h0
实值序列。因此,又有
?(ω)=H?(?ω) H0(ω)=H0(?ω),H00
(12.3.2)
由(12.2.12)式,必有Φ(ω)=Φ(?ω)。同理,我们总是选择φ(t)为实函数,因此又有
?(t)。 φ(t)=φ(?t),即尺度函数φ(t)以t=0为对称。同样的结论适用于φ
?(t)以t=1/2为对称,例如,Haar小波的尺度函数即是如此。此时要求若φ(t)和φ
?0(ω)仍是偶对称,但要增加一个移位因子,即 H0(ω)、H
?0(?ω)=ejωH?0(ω) H0(?ω)=ejωH0(ω),H
(12.3.3)
?(z),使其所形成的滤波器组为双正交滤波器现在的问题是,如何找到合适的H0(z)及H0
?(t)及ψ(t),ψ?(t)的双正交条件,即满足: 组,也即保证φ(t)、φ
??0(ω)+H0?(ω+π)H?0(ω+π)=2 H0(ω)H
也即(12.1.16a)式。习惯上将该式两边取共轭,即
- 365 -
?0?(ω)+H0(ω+π)H?0?(ω+π)=2 H0(ω)H
(12.3.4)
Cohen,Daubechies给出了不同类型的双正交小波的结构方法[42, 5],其要点是:
?(ω)是(12.3.4)式的解,若H(ω)=H(?ω),则 (1). 令H0(ω)固定,假定H000
1?′?H0(ω)=H0(ω)+H0(?ω)
2
也是(12.3.4)式的解。将该式代入(12.3.4)式即可验证。
[]
?(n)是实序列,H(ω)、H?(ω)满足(12.3.2)式,所以H(ω)、(2). 因为h0(n)、h0000?(ω)均应是实系数的三角多项式,它们可分别写成 H0
ω??
H0(ω)=2?cos?P0(cosω)
2??
2l
(12.3.5a)
和
ω???0(ω)=2?H?cos?P0(cosω) 2??
2l?
(12.3.5b)
的形式。
?0(ω)按(12.3.3)式的形式对称,则它们可表为 若H0(ω)、H
H0(ω)=2e
?jω/2
ω??cos??
2??
2l+1
P0(cosω)
2l?+1
(12.3.6a)
?0(ω)=2e?jω/2?H?cos?2??
ω?
?0(cosω) P
(12.3.6b)
的形式。
(3). 将(12.3.5)和(12.3.6)式分别代入(12.3.4)式,有
2k
2k
ω?ω???0(cosω)+??0(?cosω)=2 ?cos?P0(cosω)P?sin?P0(?cosω)P
2?2???
(12.3.7)
?;对应(12.3.6)式,k=l+l?+1。 对应(12.3.5)式,k=l+l
由于?sin
??
ω?
?0(cosω)均可以表示为sin2ω的?=(1?cosω)/2,所以P0(cosω),P
22?
- 366 -
2
函数,再令
ω?ω???2ω???
P?sin2?=P0?sin2?P? 0?sin
2?2??2???
2k
2k
(12.3.8)
则(12.3.7)式可表示为:
ω?ω?ω?ω?
?cos?P(sin2)+?sin?P(cos2)=2
2?22?2??
(4). 令y=sin
2
(12.3.9)
ω
2
,则(12.3.9)式又可表为如下的Bezout方程:
k
k
(1?y)P(y)+yP(1?y)=1
(12.3.10)
该方程和(11.4.5)式是一样的,区别只是P(y)所表示的内容。只要能求出P(y),由(12.3.8)
?(y),从而可按(12.3.5)或(12.3.6)式构造出H(z)和H?(z)。 式,即可得到P0(y)和P000
(5). (12.3.10)式的解由下式给出:
?k?1+m?mk
(12.3.11) P(y)=∑?y+yR(1?2y) ???mm=0??
这和(11.4.7)式的结果是一样的,式中R(y)是一奇对称多项式,即R(y)=?R(1?y)。
k?1
?时, H(z)、H?0(z)以n=0为对称 当 k=l+l0
?+1时,H(z)、H?0(z)以n=1/2为对称 k=l+l0
?0(y)作不同的分解可得到不同类型的双正交小波。选用不同的R,对P(y)=P0(y)P
Daubechies重点给出了基于样条函数的双正交小波的构造方法,同时也给出了H0(z)、
?0(z)长度接近相等的基于样条函数的双正交小波的构造方法,现分别给以讨论。 H
12.4 双正交样条小波
样条函数是分段光滑且在连结点处具有一定光滑性的一类函数,它在数值逼近方面获得了广泛的应用。其中基数B样条(Cardinal B-Spline)函数具有最小的支撑范围且又容易在计算机上实现,因此被认为是构造小波函数的最佳候选者之一。
而m次B样条函数Nm(t)是一阶B样条函数N1(t)自身作m?1次卷积所得到的,
N1(t)正是Haar小波的尺度函数,即
- 367 -
所以
?1
N1(t)=?
?0
0≤t
其它
(12.4.1)
?t?
N2(t)=N1(t)?N1(t)=?2?t
?0?
0≤t
(12.4.2)
?t2/2?32
t??(3/2)??
N3(t)=N2(t)?N1(t)=?4
1?(t?3)2?2?0?
0≤t
1≤t
2≤t
(12.4.3)
依次类推,有
Nm(t)=Nm?1(t)?N1(t)=Nm?2(t)?N1(t)?N1(t)
=N1(t)?N1(t)???N1(t)
(12.4.4)
?(t)等于Battle和Lemarie用上述的样条函数构造了小波,其思路是令尺度函数φ?(t)往往以t=0为对称,所以令 Nm(t)。考虑到φ
m=1
?(t)=N(t) φ1
(12.4.5)
m=2
?(t)=N(t+1)=?φ2
?1?t
?0
t≤1
其它
(12.4.6)
m=3
?0.5(t+1)2?1≤t
0≤t
21≤t
?其它0?
?(t)如图12.4.1所示。由该图可以看出,N(t)是不连续的,N(t)连续但m=1,2,3时的φ12
一阶导数不连续,而N3(t)的一阶导数是连续的,曲线已比较光滑。当m增大时,Nm(t)会变得更光滑。
- 368 -
图12.4.1 由Nm(t+1)得到尺度函数
很容易证明(12.4.4)式所决定的Nm(t)的傅里叶变换是 ?jω
m
m
???1?e
????jωm/2?
jω?
=e?sin?ω/2??ω/2?? 而对移位后的φ?(t)=Nm
(t+1),其傅里叶变换为 m
Φ?(ω)=e?jεω/2
?sin?ω/2??ω/2??
如果m为偶数,式中ε=0,若m为奇数,则式中ε=1。
分析(12.4.6)式,我们发现
φ?(t)=N(t+1)=1
φ?(2t+1)+φ?(2t)+12
22
φ?(2t?1)满足我们在第十章所讨论的二尺度差分方程。同时,可求出
- 369 -
(12.4.8)
(12.4.9)
(12.4.10)
122ω? lΦ+=+ωπ(2)cos∑332l=?∞
∞
2
(12.4.11)
是有界的。当m=3时,
?(2t?1)+φ?(2t?2) ?(t)=N(t+1)=φ?(2t+1)+φ?(2t)+φφ3
同样也满足二尺度差分方程,同理可求出
1
4343414
(12.4.12)
8131?Φ(ω+2πl)=+cosω+cos2ω ∑153030l=?∞
∞
2
(12.4.13)
?(t)可构成一个多分辨率分析。由1.7节关于正交基频域的性因此,在m=1,2,3时不同的φ
?(t)的整数移位之间不构成正交基。质,由于(12.4.11) 和(12.4.13)式的右边不等于1,因此φ?(t)“正交化”由(9.8.40)式,我们可将φ,即令
?⊥(ω)=Φ
?(ω)Φ
2??∞?
?∑Φ(ω+2πl)??k=?∞?
1
2
(12.4.14)
?(t),则φ?(t?n),n∈Z可形成一族正交基。再由第十?⊥(ω)作反变换,得尺度函数φ对Φ
章的方法可得到正交归一的小波函数。
?(t)作(12.4.14)式的正交化,而直接用N(t)作适在双正交的情况下,我们可不必对φm
?(t)作为尺度函数,如(12.4.5)~(12.4.7)式所示。这样选定φ?(t)后,Daubechies当移位后的φ
?0(y)=1。P(y)=P(y),从而得到了在双正令(12.3.11)式中的R(1?2y)等于零,并令P0?(z)的系数,即 交条件下样条小波分析滤波器H0(z)和重建滤波器H0
ω??=2l? ?(ω)=2?Hcos??,N0
2??
?N
(12.4.15a)
- 370 -
ω??=2l?+1 ?0(ω)=2e?jω/2?H?cos?,N2??
ω?
Nl+l??1
?N
(12.4.15b)
m
???+?+llmω1?2
??H0(ω)=2?cos?∑?sin??,N=2l ??22m??m=0????
(12.4.16a)
H0(ω)=2e
?jω/2
ω??
cos??
2??
Nl+l?
?l+l?+m??2ω?m
???sin?,N=2l+1 ∑??2?m??m=0?
(12.4.16b)
(12.4.15a)和 (12.4.15b)分别对应(12.3.5b)式和 (12.3.6b)式,而(12.4.16)式是(12.3.5a)式、(12.3.6a)式和(12.3.11)式的结合。
?0(ω)仅和l?有关,而和l无关;由(12.4.16)式,H(ω)不由(12.4.15)式可以看出,H0
?有关,也即H(ω)取决于N和N?。给定不同的N和N?,就可求但和l有关,而且还和l0?0(ω)。将(12.4.15)式和(12.4.8)及(12.4.9)式相比较可以看出,尺度函数出一对H0(ω)和H
?0(ω)中的N?等价,也即N??1即是得到φ(t)时由φ(t)的傅里叶变换的“阶次m”和H
。 N1(t)卷积的次数,或称之为φ(t)的“阶次”
?(t)、ψ(t)和ψ?(z)、φ(t)、φ?组合情况下H0(z)、H?(t)的系数。现给出不同N和N 0
?=1,则必有l?=0,由(12.4.15b)式,有 情况1. 令N
H0(ω)=2e
?jω/2
?ejω/2+e?jω/2?2?jω
=1+e??22??
[]
所以
[注]
,
?0(z)=H
2
(1+z?1) 2
即
?(n)={0.707,0.707} h0
令N=1,则必有l=0,由(12.4.16b)式,有
- 371 -
jω/2
+e?jω/2?0?m??1?cosω??jω/2?e?H0(ω)=2e??? ∑????22??m=0?m???
m
=
2
(1+ejω) 2
=e?jω,故和本书定义的z=ejω有区别。
[注]:Daubechies在文献中令z
所以 即
H0(z)=
2
(1+z?1) 2
h0(n)={0.707,0.707}
?(t)即是Haar尺度函数,即φ?(t)=1,对0≤t≤1,其余为零。又?=1时的尺度函数φ在N
?(t)。 ?0(z)=H0(z),由(12.2.12b)式,必有φ(t)=φ?=1时的H由于在N=N
易知在该情况下的小波函数即是Haar小波,即
?1
?
?(t)=??1 ψ(t)=ψ
?0?
0≤t
我们知道, Haar小波属正交小波,即DB1,但因为它是对称的,故又属双正交小波。
?(t),?(t),?(t)分别为1φ?(t)。我们记在该情况下的φ(t),φψ(t)和ψ1.1ψ(t)及1.1ψ1.1φ(t),
?,后面的1代表N,以下均相同。 前面的1代表N
?=1的情况下,我们再令N=3,则必有l=1,由(12.4.16b)式,有 在N
jω/2
+e?jω/2??jω/2?eH0(ω)=e??2??
3
?1+m??1?cosm?
?? ∑?m???2?m=0???
1
m
2jω
=e+3+e?jω+e?j2ω
8
[]
?4?ejω?e?jω??? ?2??
所以
11111?1?
H0(z)=2??z2+z++z?1+z?2?z?3?
16221616??16
?(t)不变,1.3φ(t),1.3ψ(t)及1.3ψ?仍为1,所以, 1φ?(t)可由上一节的公式推出。因为N
但MATLAB中的wavefun.m文件可用来产生这些函数,如图12.4.2所示。
- 372 -
?=2,则l?=1,由(12.4.5a)式,有 情况2:令N
?0(ω)=2(cosω)2=2(ejω+2+e?jω) H
24
?0(z)=H
2
(z+2+z?1) 4
再令N=2,则l=1,由(12.4.16a)式,有
?1+m??2ω?
H0(ω)=cos∑?m???sin2? 2m=0?????
2ω
1
m
=
即
H0(z)=
2jω
(e+e?jω+2)(4?ejω?e?jω) 8
1131?1?
2??z2+z++z?1?z?2?
8444?8?
?0(z)及H(z),?和N的组合下的H按此方法类推,读者不难得出在不同N如表12.4.10?和N的合适组合是: 所示。N
?=1, N
?=2, N?=3, N
N=1,3,5 N=2,4,6,8 N=1,3,5,7,9
?0(z)和H0(z),特别是H0(z),若不考虑前面的2,其分母的由该表可以看出,H
?(t)都?值下,φ系数都是2的整次幂,因此有利于在计算机上快速实现。此外,在不同的N
?(t),ψ(t)及是精确已知的,这些都是样条点双正交小波的优点。在不同组合下的φ(t),φ
?=1,?(t)如图12.4.2~12.4.6所示。其中图12.4.2给出的是NN=3和5时的φ(t)及ψ(t)。ψ
?=1,N=3,图中左边标的“bior 1.3 phi-D”,即是指双正交小波的尺度函数φ(t),对应N
“D”代表分解,右边图标的“psi”指的是“ψ”。以下各图的标法均相同。另外,图中的横坐标是MATLAB按正的坐标求出的,这和12.3节所给出的.支撑范围有区别。
- 373 -
?=1,N=3,5时用于分解的φ(t)和ψ(t) 图12.4.2 N
?(t)(图中标为phi-R,R代?=1,N?=2,N?=3和N?=4时的φ图12.4.5给出的是N
?(t)即是Haar尺度函数,它即是图12.4.1的(a)图。而N?=1时的φ?=2表重建)。显然,N
?(t)即是图12.4.1(b),N?(t)即是图12.4.1(c)。?=3时的φ时的φ它们分别是N1(t),N2(t)和?(t)应是N(t)。注意,φ?(t)只和N?=4时的φ?有关,而和N无关。 N3(t)。显然,图中N4
?=1,?=2,N=2;N?=3,N=3和N?=3,N=5图12.4.6给出的是N N=3;N?有关,而且也和N有关。
?(t)。显然,ψ?(t)不仅和N时的ψ
- 374 -
?=3,N=3,5,7,9时用于分解的φ(t)和ψ(t) 图12.4.4,N
- 376 -
=1, N=3;N=2,N=2;N=3,N=3和N=3,N=5时用图12.4.6,N
?(t) 于重建的ψ
?=1,2,3时N取不同值时H?0(z)/2,H0(z)/2的系数,由于表12.4.1给出了在N
?=2,N=6和8,N?=3,N=5,7和9时的H0(z)的系数过长,故表中没有列入,在N
详细数据可参考文献。
?0(z)及H0(z)的长度差别甚大,由以上讨论可知,按(12.4.15)或(12.4.16)式构造出的H
?(n)的长度决定了h(n)的长度。这样,且N越大,这一差别越明显。由(12.1.15a)式,h01
一对分解滤波器H0(z)和H1(z)的长度将会有着明显的不同。这在一些应用中将会带来不便和麻烦,特别是在语音和图像处理方面。
?0(z)和H(z)长度不同的原因在于对P(y)=P(y)P?0(y)的分解,即(12.4.15)和H00
(12.4.16)式是在假定P0(y)=1,P(y)=P0(y)情况下得到的。若对P(y)作另外形式的分解,即
?k?1+m?m
???P(y)=∑??m??y=P0(y)P0(y),k=l+l,P0(y)≠1 (12.4.17)
m=0??
k?1
?0(y)分别代入(12.3.5)和(12.3.6)式,则可得到保证在双正交条件下且长度然后将P0(y)和P
?(z)和H(z)。Daubechies令 接近的H00
P(y)=A
∏(y?yi)∏(y2?2Re[zi]y+yi)
2
j=1
j=1
J1J2
(12.4.18)
表12.4.1
?(z),H(z)的系数 H00
- 377 -
式中yj(j=1~J1)为P(y)的一阶实根,yi,i=1~J2是P(y)的共扼复根,然后在保证
?(y)系数始终为实数的情况下,考虑yj,y对P(y)和P?(y)的分配。 P0(y)、P00i0
?=4,N=4,即l?=l=2的情况下,有 在N
3
?3+m?m?3+m??2ω?
P(y)=∑??m??y=∑??m???sin2?
?m=0?m=0????
3
m
即
P(z)=?5z3+40z2?131z+208?131z?1+40z?2?5z?3/16
[]
它有两个实根,即z1=0.3289,z2=3.0470,两对共扼复根,即z3,4=0.2841±j0.2432,
- 378 -
z5,6=2.0311±j1.7390。又由于
??
?cos?=?cos?
2?2???
令z=e
jω
ω?
2l
ω?
2l?
?ejω/2+e?jω/2??ejω+2+e?jω?
=?=???24????
?1
42
,则上式等效为z+4z+6+4z
[
2
?0(z),它即属于H也属于H0(z)。+z?2/16。
]
?(z),则H?(z)的长度为7。将余下的z3,4和z5,6赋若将上面分解出的零点z1和z2赋给H00
给H0(z),则H0(z)的长度为9。这样,二者的长度基本相等,即满足了(12.3.9)式,且又具有对称性。
?=5,N=5时,?=4,N=4存在两种分解方式。相应的系数如表12.4.2所示。N当N
?,φ,ψ?,φ,ψ?=N=5时的φ?和ψ如图12.4.7所示,N?和ψ如图12.4.8所示。时的φ
?(t),ψ?,,N=4时用于分解的φ(t),ψ(t)及用于重建的φ?(t) 图12.7.7 N
以上除双正交小波外,Daubechies还构造了接近于正交基的双正交基函数,详细内容见文献,此处不再讨论。
- 379 -
?(t),ψ?,,N=5时用于分解的φ(t),ψ(t)及用于重建的φ?(t) 图12.4.8 N
表11.4.2 具有接近长度的双正交小波对应的滤波器系数
- 380 -
12.5 正交小波包
第十章讨论的多分辨率分析将L2(R)空间逐层进行分解,如将V0分成V1和W1,再将V1分成V2和W2,?,其中V0=V1⊕W1,V1=V2⊕W2,及V0=⊕+Wj。对同一尺度j,Vj
j∈Z
是低频空间,Wj是高频空间,因此,信号x(t)在Vj中的展开系数aj(n)反映了信号的“概貌”,而在Wj中的展开系数dj(n)反映了信号的“细节”,也即x(t)的小波系数。由于这种分解具有恒Q性质,即在高频端可获得很好的时域分辨率而在低频端可获得很好的频域分辨率,因此,这种分解相对均匀滤波器组和短时傅里叶变换有着许多突出的优点,因此获得了广泛的应用。
但这种分解仅是将Vj逐级往下分解。而对Wj不再作分解。将W1和W2相比,显然,W1对应最好的时域分辨率,但是有着最差的频域分辨率。这在既想得到好的时域分辨率又想得到好的频域分辨率的场合是不能满足需要的。当然,在任何情况下,时域-频域分辨率之间都要受到不定原理的制约,但是,我们毕竟可根据工作的需要在二者之间取得最好的折中。例如,在多分辨率分解的基础上,我们可将Wj空间再作分解,如图12.5.1所示。
j=0j=1j=2
j=3
图12.5.1 V0空间的逐级分解
在该图的分解中,任取一组空间进行组合,如果这一组空间:①能将空间V0覆盖;②
- 381 -
相互之间不重合,则称这一组空间中的正交归一基的集合构造了一个小波包(wavelet packet)。显然,小波包的选择不是唯一的,也即对信号分解的方式不是唯一的。如在图12.5.1中,我们可选择
① V31,W31,V32,W32,V33,W33,V34,W34; ② V31,W31,W21,V22,W22; ③ V1,V22,W22
等不同空间来组合,它们都可覆盖V0,相互之间又不重合。如何决定最佳的空间组合及寻找这些空间中的正交归一基便是小波包中的主要研究内容。
图12.5.1的空间分解可用图12.5.2的滤波器组来实现。注意,在实现各级的卷积时,图中滤波器H0,H1的系数要事先翻转,即将hi(n)变成hi(?n),i=0,1 。
图12.5.2 图12.5.1的滤波器组实现
由该图可以看出,基于小波包的信号分解也是用一对滤波器H0(z)和H1(z)来实现的。在第十章的多分辨率分析中,我们详细讨论和论证了在Vj和Wj中分别存在正交归一基φj,k(t)和ψj,k(t),它们和共轭正交镜像滤波器组H0(z)、H1(z)有如下关系:
?t
φ?j?2∞
??t??=2∑h0(k)φ?j?1?k? ??2?k=?∞
(12.5.1a)
- 382 -
?tψ?j?2∞
??t?=φhk2()??j?1?k? ∑1??2?k=?∞
(12.5.1b)
当j=0时,有
φ(t)=2∑h0(k)φ(2t?k)
k=?∞
∞
(12.5.2a) (12.5.2b)
ψ(t)=2∑h1(k)φ(2t?k)
k=?∞
∞
此即二尺度差分方程。式中,h0(k)、h1(k)有如下关系:
h1(k)=(?1)kh0(1?k)
(12.5.3)
在上述的多分辨率分析中,当将Vj分解成Vj+1和Wj+1时,Vj中的正交归一基基φj,k(t)产生了两个正交归一基φj+1,k(t)和ψj+1,k(t),它们分别属于Vj+1和Wj+1,生成办法即是(12.5.1)式的二尺度差分方程。由此我们可以设想,在图12.5.1中,将W1分解生成W21和W22时,
W1中的正交归一基ψ1,k(t)也将会依照二尺度差分方程分别生成W21和W22中的正交归一
基。如果这一结论正确,则图12.5.1中的各个子空间将都存在正交归一基。文献证明了这个一般结论。该结论可由下述定理来描述:
定理12.5
令θj,k(t)是空间Uj中的正交归一基,h0(k),h1(k)是满足(12.5.3)式的
一对共轭正交滤波器,令
θ
0j+1
(t)=
k=?∞
∑h(k)θ(2
∞
?j
t?k)
(12.5.4a)
(12.5.4b)
θ
1j+1
(t)=
k=?∞
∑h(k)θ(2
1
∞
?j
t?k)
则
1
{θ0j+1,k(t),θj+1,k(t)},k∈Z
是Uj中的正交归一基。
该定理的证明见文献。显然,令Uj,Uj分别是θj+1,k(t)和θj+1,k(t)所产生的空间。自然有
- 383 -
1
1
1
U0 j+1⊕Uj+1=Uj+1
(12.5.5)
在图12.5.1中,V0中有正交基φ0,k(t),V1中有正交基φ1,k(t),W1中有正交基ψ1,k(t)。按定理12.5,由ψ1,k(t)可生成W21,W22中的正交基。依次类推,我们可得到图12.5.1中任一子空间中的正交归一基。
对于给定的尺度j,在图12.5.1中,共有2j个子空间。为了讨论的方便,我们将图中的子空间统一标记为Wj,Wj,?,Wj
2p
1
2j?1
。如j=3,共有W3,W3,?,W3个子空间。
2p+1
017
显然,Wj对应每一次剖分的低频部分,而Wj
p
对应其高频部分,p=0,1,?,2
j?1
?1。
令每一个子空间的正交归一基为ψj,k(t)。由(12.5.1) 及(12.5.4)式,有
ψψ
2p
(2
?(j+1)
t)=2
k=?∞
2p+1
∑h(k)ψ(2
p
0∞
p
1
k=?∞
∞
?j
t?k)
(12.5.6a) (12.5.6b) (12.5.7a)
(2
?(j+1)
t)=2
∑h(k)ψ(2
?j
t?k)
及
2h0(k)=2p(2?(j+1)t),ψp(2?jt?k)
2h1(k)=2p+1(2?(j+1)t),ψp(2?jt?k)
(12.5.7b)
由(12.5.5)式,有
2p+1Wj2+p=Wjp 1⊕Wj+1
(12.5.8)
显然,当p=0时,W0=W0即是空间V0,基函数ψp(2?jt?k)=ψ0(t?k)即是V0中的正交归一基φ(t?k),也即尺度函数。由图12.5.1也可看出,Wj即是我们在第十章讨论过的多分辨率分析中的空间Wj,因此Wj中的正交归一基ψ1(2?jt?k)即是Wj中的正交归一基ψj,k(t),也即小波函数。由(12.5.6)式,令j=0,则
- 384 -
1
p0
1
ψψ
2p
(t)=(t)=
2
k=?∞
2p+1
∑h(k)ψ(2t?k)
p
0k=?∞
∞
(12.5.9a) (12.5.9b)
2
∑h(k)ψ(2t?k)
p
1
∞
按照上述思路,只要我们给定了V0中的尺度函数φ(t)及相应的小波函数ψ(t),由(12.5.6),或(12.5.9)式即可递推地求出小波包分解中各个子空间中的基函数ψj,k(t)和ψj,k(t)。
例12.5.1 由(12.5.6)式,有
对Harr小波,h0(0)=h0(1)=
2p
2p+1
111,h1(0)=,h1(1)=? 222
ψ
2p
j+1
(t)=
2∑h0(k)ψjp(t?2jk)
k=0
1
ψ
2p+1j+1
(t)=
2∑h1(k)ψjp(t?2jk)
k=0
1
当j=0时,
(2t)+ψ00(2t?1) ψ10(t)=ψ0
1
(t)=ψ00(2t)?ψ00(2t?1) ψ1
当j=1时,
ψ2(t)=ψ1(2t)+ψ1(2t?1)
1
(t)=ψ10(2t)?ψ10(2t?1) ψ22
(t)=ψ11(2t)+ψ11(2t?1) ψ23
(t)=ψ11(2t)?ψ11(2t?1) ψ2
103
ψ10(t),ψ1(t)的宽度都是2T,而ψ2(t)~ψ2(t)的宽度为4T。j=1,2,3时的ψjp(t)分别
示于图12.5.3a,b和c。
- 385 -
图12.5.3
1
(a) ψ10(t),ψ1(t)
由Harr小波生成的小波包
03
(b) ψ2(t)~ψ2(t)
(c) ψ3(t)~ψ3(t)
07
例12.5.2
令V0空间中的φ(t)为“DB5”小波对应的尺度函数,当j=3时,求出
07
W30~W37中的小波基ψ3(t)~ψ3(t)如图12.5.4所示。
- 386 -
07
图12.5.4 j=3时,由“DB5”小波生成的ψ3(t)~ψ3(t)
上述两例的分解过程可以形象地表为一个二进制的树结构,如图12.5.5所示。图中结点处的数值即为(j,p)。MATLAB中的plottree命令文件可画出此结构图。
令x(t)∈L2(R),则
a0(n)=x(t),φ(t?n)=x(t),ψ0(t?n)
(12.5.10)
是x(t)在空间V0=W0中的“概貌”。我们把它当作一个树状滤波器组的输入信号。在小
- 387 -
波
图12.5.5 j=3时的二进制树结构图
包的分解中,对任意的结点(j,p),则
(0,0)
(1,0)
(2,0)
(2,1)
(2,2)
(1,1)
)
(2,3)
(3,0)(3,1)(3,2)(3,3)(3,4)(3,5)(3,6)(3,7)
djp(n)=x(t),ψjp(t)
为x(t)在该结点(或子空间Wj)处的小波包系数,它是x(t)和基函数ψ(2t?n)作内积的结果。下述定理给出了小波包系数的快速计算方法。
定理12.6
在小波包的分解中,在结点(j+1,p)处的小波包系数由下式给出
2pj+1
pp?j
d
(k)=d(k)?0(2k)=
pj
m=?∞
∑d
∞
∞
pj
(m)h0(m?2k)
(12.5.11a)
d
2p+1j+1
(k)=d(k)?1(2k)=
p
pj
m=?∞
∑d
pj
(m)h1(m?2k)
(12.5.11b)
而在结点(j,p)处的小波包系数dj可由下式重建:
?p?2p+1
?+djp(k)=d2khkd()()j+1j+1(k)?h1(k) 0?2p+1
2p
2p+1
(12.5.12)
式中dj+1(k)和dj+1(k)分别是dj+1(k)和dj+1(k)每两个点插入一个零后所得到的序列。
该定理的证明类似于定理10.6和定理10.7,此处不再讨论。j=2时的分解与重建如
?2p
- 388 -
图12.5.6a和b所示。图中d0(n)即是a0(n)。在图(a)中,当实现各级的卷积时,图中滤波器H0,H1的系数同样要事先翻转,即将hi(n)变成hi(?n),i=0,1 。
02
122232
图12.5.6 基于滤波器组的小波包分解与重建
(a)小波包分解, (b) 小波包重建
例12.5.3
信号x(t)是MATLAB中所给的信号noisdopp.mat,如图12.5.7a的第一
p
个图所示。令j=3,我们可得到对x(t)作小波包分解后的各个子空间的系数dj。由此可看出信号x(t)在各个频带上的时域行为。图中所用小波均是“DB5”正交小波。图a对应
j=1,p=0,图b对应j=2,p=0,1,图c对应j=3,p=0,1,2,3。
图12.5.7信号的小波包分解,
0710123
(a) 原信号x(t),d10和d1;(b) d2,d2,d2及d2;(c) d3。 ~d3
由该图可以看出,当j=1时,d1(n)反映了x(t)的概貌,它相当于多分辨率分析中
1
的V1空间的a1(n)。而d1(n)相当于W1空间的d1(n),它是x(t)中的噪声分量,也即高频
成分;当j=2时,d2(n)仍是信号的概貌,但比d1(n)含有较少的噪声,d2(n),d2(n)及d3(n)也属于噪声分量,当j=3时,反映概貌的d3(n)已几乎不含噪声,d3~d3属于噪声成分,但其幅值已很小。
上述分解是将j=1,2,及3时的分解系数全部求出。实际上,根据小波包的定义,我们只需取部分互不重迭且又能覆盖V0的子空间即可,这即是所谓“最佳小波包”的选择问题。一个最佳小波包的选择取决于三个因素:
(1) 信号本身的性质; (2) 信号分解的目的; (3) “最佳”原则的选择。
显然,一个最佳的小波包应使信号x(t)在其各个子空间中的投影(即dj)尽可能地大。至于说由哪几个空间组成一个最佳的小波包,显然取决于信号x(t)能量随频率的分布。从应用的角度看,分解的目的若是为了去噪,那就应舍弃噪声能量较大的子空间,选择不但信号能量大而且噪声也小的空间;分解的目的若是为了数据的压缩,则应选择信号能量集中的空间。但是,无论何种目的,在决定“最佳小波包”的过程中,我们总要确定一个“代价函数”,从而使在各种小波包选择的可能中,所选择的一种具有最小的代价函数。
至今,人们已提出了代价函数的多种选择方法,如“编码率-失真(R-D)指标”,“范数(Norm)”判据等。若令xi为信号x(t)在某一子空间正交基“Shannon熵判据”,上的投影,则定义
p
3
1
7
0012
E1(x)=?∑xi2lgxi2
i
(12.5.13)
为x的Shannon熵,定义
E2(x)=∑xi
i
p
=xp,1≤p
- 391 -
p
(12.5.14)
为x的l范数。此外,还可定义,
p
E3(x)=∑logxi2
i
(12.5.15)
为x的对数能量熵。这样,对给定的信号x,我们可求出它在每一个子空间中的“熵”或范数,并把它们作为代价函数来决定小波包的选择。文献以j=3为例介绍了一种自底向顶的快速搜索方法。如图12.5.8所示。上方的空间为“母空间”,下方的空间为“子空间”。在每一个空间都标上了由(12.5.13)~(12.5.15)式中任一式求出的代价函数。如果子空间的代价总和小于母空间的代价,这说明这一分解是值得的,因此,该子空间应预以保留。反之,若子空间的代价总和大于母空间的代价,则这一分解是不适当的,应加以放弃,即保留“母空间”。
由于W30加上W31中的代价为3,小于W20的代价11,所以由W20至W30和W31的分解应加以保留;W3加上W3的代价为7,小于W2的代价12,所以这一分解也应保留。同理,
2
3
1
W34加W35的代价为11,小于W22的代价13,应加以保留,但W36加W37的代价为15,大
于W2的代价14,因此这一分解应该放弃。依次往上类推,最后选中的应是子空间W3,W3,
45
。其中,选中的W3,W3W32,W33及W11,由它们构成在所给定意义上最佳的“小波包”
301
和W2由于可由W1所覆盖,故可以放弃。
31
1
w3
w03w23w33w43w53w63w73
图12.5.8
MATLAB中的wavelet Toolbox中有着丰富的有关小波包的命令,其中
Wentropy.m:用来计算所给定信号在各个子空间分解后的熵;
Besttree.m:在指定所用的熵的基础上用来确定最佳小波包的选择。实际上给出的是最佳树的选择;
Wpdec.m:对给定的信号x(t)和熵,确定分解的树结构和分解后的数据结构; Wpcoef.m:由上面的数据结构得到在各子空间的分解系数(即dj); Wprec.m:由上述分解系数重建信号x(t)。
令x(t)仍然为MATLAB中的noisdopp.mat。如图12.5.5(a)的第一个图
p
最佳小波包选择 图中阴影部分为最后所选择的小波包
例12.5.4
所示。用如下三个命令可得到对x(t)的最佳二进制树结构:
[wpt,wpd]=wpdec(x,3,‘wname’);
%:对x作分解,wpt中含树结构,wpd中含数据结构,使用“Shannon”熵; %:‘wname’指定要选的小波。本例中,wname=DB1。3表示j=3。
[t1,d1,e,n]=besttree(wpt,wpd);
%:找最佳树结构,存t1中,d1含新的数据结构,e中含各结点的熵,n为波合并%:的结点索引序号; [t,d]=wp2wtree(t1,d1);
%:抽取最佳树结构t,并得到相应的数据结构; %:画出该树结构。如图12.5.9所示。
Plottree(t);
(3,01,1)
图12.5.9 j=3时的最佳树结构。(a) 最佳树结构, (b) 分解的区间
由该图可以看出,结点(1.1)没有再分解,结点(2.1)也没有再分解,所分解的空间如图(b)所示。因此,我们可选W3,W3,W21及W11作最佳小波包。当然,也可选W20,W21及
1
W11作最佳小波包,但这对应的是j=2的分解,而不是j=3的分解。
- 394 -
标签:胡广书,信号处理,教程