小波变换详解

第10章 小波变换与JPEG 2000编码
之小波变换
虽然基于DCT 的JPEG 标准的压缩效果已经很不错，但在较高压缩比时会出现明显的马赛克现象，且不能渐进传输。

为了适应网络发展的需要，JPEG 于2000年底推出了采用DWT (Discrete Wavelet Transform 离散小波变换)的JPEG 2000标准。

小波变换是1980年代中期发展起来的一种时频分析方法，比DCT 这样的傅立叶变换的性能更优越，被广泛应用于调和分析、语音处理、图像分割、石油勘探和雷达探测等等方面，也被应用于音频、图像和视频的压缩编码。

本章先介绍小波变换的来龙去脉，然后分别介绍连续小波变换、离散小波变换、Haar 小波变换和整数小波变换，最后介绍JPEG 2000的编码算法和标准。

10.1 小波变换
小波变换(wavelet transform)是傅立叶变换的发展，中间经历了窗口傅立叶变换。

原始数据一般是时间或空间信号，在时空上有最大分辨率。

时空信号经傅立叶变换后得到频率信号，在频域上有最大分辨率，但其本身并不包含时空定位信息。

窗口傅立叶变换通过对时空信号进行分段或分块进行时空-频谱分析，但由于其窗口的大小是固定的，不适用于频率波动大的非平稳信号。

而小波变换可以根据频率的高低自动调节窗口大小，是一种自适应的时频分析方法，具有多分辨分析功能。

本节先讨论小波变换与（窗口）傅立叶变换的关系，然后依次介绍连续小波变换、离散小波变换、Haar 小波变换和第二代小波变换（整数小波变换）。

10.1.1 傅立叶变换与小波变换
傅立叶变换(Fourier transform)是法国科学家Joseph Fourier 发表
于1822年的他在用无穷三角级数求解热传导偏微分方程时所提出的一种数学方法，它可将时空信号变换成频率信号。

鉴于傅立叶变换不含时空定位信息，（1971年的诺贝尔物理学奖
获得者）匈牙利人Dennis Gabor 于1946年提出窗口傅立叶变换（window Fourier transform ）。

可以用于时频分析，但是窗口大小是固定的。

1984年法国的物理学家Jean Morlet 和A. Grossman ，在进行石油勘探的地震数据处理分析时，又提出了具有可变窗口的自适应时频分析方法——小波变换（wavelet transform ）。

Joseph Fourier
● 傅立叶变换
傅立叶变换(Fourier transform)是1807年法国科学家Joseph Fourier 在研究热力学问题时所提出来的一种全新的数学方法，当时曾受到数学界的嘲笑与抵制，后来却得到工程技术领域的广泛应用，并成为分析数学的一个分支——傅立叶分析。

原始的多媒体数据一般为时空信号，在时空上有最大分辨率，并可利用时空上的相关性进行数据压缩。

Fourier 变换可将时空域中的多媒体信号映射到频率域来研究，即更符合人类感觉特征，也可以利用信号在频率域中的冗余进行数据压缩。

Fourier 变换所得的频率信号，在频率域上有最大分辨率，但其本身并不包含时空定位信息。

时空信号：
f (t ), t ∈(-∞, ∞) （一维时间信号，参见图10-1）
f (x , y ), x , y ∈(-∞, ∞) （二维空间信号）
图10-1 音频信号的时间波形图
Fourier 变换，F (w )为频率信号：
⎰∞
∞--=dt e t f w F jwt )()( （参见图10-2） dy dx e y x f v u F vy ux j ⎰⎰∞
∞-∞∞-+-=)(),(),(
图10-2 音频信号的频率图
● 窗口傅立叶变换
虽然基于Fourier 变换的频谱分析，在需要信号分析及数据处理的物理、电子、化学、生物、医学、军事、语音、图像、视频等众多科学研究与工程技术的广阔领域得到了非常广
泛和深入应用，但对既需要频谱分析又要求时空定位的应用，如雷达探测、语音识别、图像处理、地震数据分析等等，Fourier 分析技术就显得力不从心了。

为了弥补Fourier 变换不能时空定位的不足，工程技术领域长期以来一直采用D.Gabor 开发的窗口Fourier 变换（短时Fourier 变换），来对时空信号进行分段或分块的时空-频谱分析（时频分析）。

窗口Fourier 变换：（参见图10-4）
⎰∞
∞---=dt e t g t f w F jwt g )()(),(ττ 其中，g 为窗口函数（参见图10-3）。

图10-3 音频处理中常用的几种窗口函数
图10-4 音频信号的三维频谱图
虽然窗口Fourier 变换能部分解决Fourier 变换时空定位问题，但由于窗口的大小是固定的，对频率波动不大的平稳信号还可以，但对音频、图像等突变定信号就成问题了。

本来对高频信号应该用较小窗口，以提高分析精度；而对低频信号应该用较大窗口，以避免丢失低频信息；而窗口Fourier 变换则不论频率的高低，都统一用同样宽度的窗口来进行变换，所以分析结果的精度不够或效果不好。

迫切需要一种更好的时频分析方法。

小波变换
近二十年来发展起来的小波(wavelet)分析正是这样一种时频分析方法，具有多分辨分析
功能，被誉为数学显微镜。

它是继一百多年前发明傅立叶分析之后的又一个重大突破，对许多古老的自然学科和新兴的高新技术应用学科都产生了强烈冲击，并迅速应用到图像处理和语音分析等众多领域。

1）函数展开与积分变换
小波分析是傅立叶分析的发展，是分析数学的一个新分枝，高等数学中的微积分（数学分析）就是分析数学的基础。

与幂级数、三角级数或傅立叶级数等一样，小波分析研究用一组简单函数，如{x n }、{sin nx , cos nx }等，来表示任意函数，如
∑∞
==0
)(n n n x a x f （幂级数）
∑∑∞-∞=∞==⎪⎭⎫ ⎝⎛++=n x l n j n n n n e c l x n b l x n a a x f πππ10sin cos 2)(（三角级数/傅立叶级数）
其中
1 ,sin cos ),(2
1 ),(21-=+=+=-=-j j e jb a c jb a c j n n n n n n θθθ 被表示的函数的全体构成一个函数空间（一种函数的集合），而表示这些函数的函数族{x n }与{sin nx , cos nx }等则为函数空间的基底。

函数展开式中的系数为该函数在函数空间中相对于此基底的坐标，对应于函数空间的一个点。

这相当于将函数从原来的域变到新的域，如三角级数将时空域的函数变换到频率域。

为了求得展开式的系数，需要对原函数求微积分，如幂级数中的
!
)0()(n f a n n = 三角级数中的
⎰⎰--==
l l n l l n dx l
x n x f l b dx l x n x f l a ππsin )(1 ,cos )(1 和傅立叶级数中的 ⎰--=l l x l n j n dx e x f l c π)(1
若f (x )不是以2 l 为周期的函数，在上式中改记x 为t 、
w l n =π，并让∞→l ，则得Fourier 变换：
⎰∞
∞--=dt e t f w F jwt )()( 这是一种复变函数的广义积分，也是一种积分变换。

2）小波的发展
自从近两百年前Joseph Fourier 在研究热力学问题提出Fourier 分析以后，长期以来许多数学家一直在寻找更广泛函数空间的性能更好的基底函数族，工程技术领域也一直在寻找更好的时频分析方法，但收获甚微。

1984年法国的年轻的地球物理学家Jean Morlet 在进行石油勘探的地震数据处理分析时与法国理论物理学家A.Grossman 一起提出了小波变换(wavelet transform, WT)的概念并定义了小波函数的伸缩平移系：
⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧⎪⎭⎫ ⎝⎛-a b x a ψ||1， 但并没有受到学术界的重视。

直到1986年法国大数学家Yves Meyer 构造出平方可积空间L 2的规范正交基——二进制伸缩平移系：
⎭
⎬⎫⎩⎨⎧-=--)2(2)(2,k x x j j k j ψψ
小波才得到数学界的认可。

1987年正在读硕士的Stephane Mallat 将自己熟悉的图像处理的塔式算法引入小波分析，提出多分辨分析的概念和构造正交小波的快速算法——Mallat 算法。

1988年法国女科学家Inrid Daubechies 构造出具有紧支集的正交小波基——Daubechies 小波。

1990年美籍华裔数学家崔锦泰和武汉大学的数学教授王建忠又构造出基于样条函数的单正交小波函数——样条小波。

1992年Daubechies 在美国费城举行的CBMS-NFN 应用数学大会上作了著名的《小波十讲Ten Lectures on Wavelets 》报告，掀起了学习与应用小波的高潮。

1994年Wim Swelden 提出了一种不依赖于Fourier 变换的新的小波构造方法——提升模式(lifting scheme)，也叫第二代小波或整数小波变换。

3）连续小波变换
连续小波变换(CWT = Continuous wavelet transform)的定义为：
⎰∞
∞-⎪⎭
⎫ ⎝⎛-=dx a b x x f a b a W f ψ)(||1
),( 其中，a 为缩放因子（对应于频率信息），b 为平移因子（对应于时空信息），)(x ψ为小波函数（又叫基本小波或母小波），)(x ψ表示)(x ψ的复共轭。

连续小波变换的过程可参见图10-5。

图10-5 连续小波变换的过程
小波变换的特点有：（参见图10-6）
⏹ 时频局域性、多分辨分析、数学显微镜
⏹ 自适应窗口滤波：低频宽、高频窄
⏹ 适用于去噪、滤波、边缘检测等
图10-6 窗口傅立叶变换与小波变换的时频特征
如同三角函数sin x 和cos x 及e -jx 可以缩放构成函数空间的基底{sin nx , cos nx }及{ e -jwx }一样，母小波也可以缩放和平移而构成函数空间的基底：
⎭⎬⎫⎩⎨⎧-=--)2(2)(2,k x x j j k j ψψ及⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧⎪⎭⎫ ⎝⎛-a b x a ψ||1 与傅立叶变换不同，小波变换的结果有两个参数，多了一个可以表示时空位置信息的平移因子，所以其图示为一个二维曲面。

图10-7/8是Mallat 构造的一组典型数据的曲线及其连续小波变换曲面的二维与三维图示：
图10-7 Mallat 数据及其连续小波变换的二维图示
图10-8 Mallat 数据及其连续小波变换的三维图示
4）小波函数
小波变换与傅立叶变换比较，它们的变换核不同：傅立叶变换的变换核为固定的虚指数函数（复三角函数）e -jwx ,而小波变换的变换核为任意的母小波)(x ψ。

前者是固定的，而后者是可选的，实际上母小波有无穷多种，只要)(x ψ满足下列条件即可。

小波函数需满足的条件：
绝对可积且平方可积，即2
1L L ⋂∈ψ
⏹ 正负部分相抵，即0)(=⎰∞
∞-dx x ψ（即0)0(ˆ=ψ
） ⏹ 满足允许条件(admissible condition)，即∞<⎰∞∞-ωω
ωψd 2)(ˆ（广义积分收敛） 其中)(ˆωψ
为)(x ψ的傅立叶变换 常见的小波函数有：
⏹ Haar 小波(Alfred Haar ，1910年)：⎪⎩
⎪⎨⎧<≤-<≤=其他 ,015.0 ,15.00 ,1)(x x x ψ，参见图10-9。

图10-9 Haar 小波函数及其Fourier 变换
⏹ 墨西哥草帽(Mexican hat)小波：2222)(x e dx d x -=ψ，参见图10-10。

图10-10 墨西哥草帽小波函数及其Fourier 变换
⏹ Morlet 小波(Jean Morlet ，1984年)：5,)(2 2
≥⋅=-C e e x x x jC ψ，参见图10-11。

图10-11 Morlet 小波函数(C=5)及其Fourier 变换
除了Haar 小波外，其他紧支集小波都不是初等函数，有的小波函数是用导数/积分或微分方程/积分方程来定义，有的小波用其傅立叶变换定义，有的小波甚至没有解析表达式，而只是一些数字解，很多小波为复函数，所以不太直观。

可以把小波与三角函数中正弦波加以比较（参见图10-12）。

图10-12 小波与正弦波
离散小波变换
将连续小波变换的缩放因子a 离散化，得到二进小波变换；再将其平移因子b 也离散化，就得到离散小波变换。

1) 二进小波变换与滤波器
为了适应数字信号处理，需要将小波变换离散化。

可以先进行缩放因子的离散：若小波函数ψ满足
1|)2
(ˆ|2=∑∈Z k k ωψ，
则称ψ为基本二进小波。

在连续小波变换中，若ψ为基本二进小波，则令a = 2k ，得到二进小波变换：
⎰∞∞-⎪⎭
⎫ ⎝⎛-=dx b x x f b f W k k k 2)(21)(2ψ 为了构造基本二进小波，可设φ满足：
∑∞==122|)2(ˆ||)(ˆ|j j
ωψωφ
可推出1|)0(ˆ|2=φ
，则φ大体上相当于一个低通滤波器，因此，φ(2x )的通道比φ(x )的宽，可设φ满足如下的双尺度方程：
∑∈-=Z
n n n x h x )2(2)(φφ
其Fourier 变换为：
⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛=2ˆ2)(ˆωφωωφH ，其中∑∈-=Z
n n i n e h H ωω)( 为低通滤波器。

由1|)0(ˆ|2=φ
，可得H (0) = 1 即∑h n = 1。

若设
22|)(|1|)(|ωωH G -=，其中∑∈-=Z
n n i n e g G ωω)(
则G 为高通滤波器，有
∑∈-=Z
n n n x g x )2(2)(φψ
其Fourier 变换为：
()()ωφωωψ
ˆ)2(ˆG = 因0)0(ˆ=ψ
且1|)0(ˆ|2=φ，得G (0) = 0 即∑g n = 0。

例如（B2滤波器），若取φ为二次B 样条，则
)33(8
1)2(cos )(222ωωωω
ωωi i i i e e e e H ---+++== 可得h n = h 1- n ，h 0 = h 1 = 3/8 = 0.375，h -1 = h 2= 1/8 = 0.125，其余h n = 0；
因G 不唯一，可令)()(ωω--=G G ，g n = -g 1- n ，解得-g 0 = g 1 = 0.5798，-g -1 = g 2 = 0.0869，-g -2 = g 3 = 0.0061，其余g n =0。

又例如（B3滤波器），若取φ为中心三次B 样条，则
)464(16
1)2(cos )(224ωωωωω
ωi i i i e e e e H --++++== 可得h n = h - n ，h 0 = 3/8 = 0.375，h -1 = h 1= 1/4 = 0.25，h -2 = h 2= 1/16 = 0.0625，其余h n = 0；
似上例可得g n = -g - n ，-g -1 = g 1 = 0.59261，-g -2 = g 2 = 0.10872，-g -3 = g 3 = 0.00008，其余g n 为0。

2) 离散小波变换
下面再将二进小波变换中的平移因子也离散化：令b = n 2k ，则可得离散小波变换：
⎰∞
∞----=dx n x x f n f W k k k )2()(2)(2/2ψ 可以用前面所讲的滤波器系数改写成如下循环形式：
∑∑---=-=--k j k k
j k k n f S h n f S k n f S g n f W j j j j )2
()()
2()(12212211
其中，)()(02n f n f S =，D = Wf 为差——高频部分，A = Sf 为剩余——低频部分，h k 与g k 为上面讲过的滤波器)()(ωωG H 与之系数。

可以写出正反离散小波变换的具体算法如下：
⏹ 正变换（分解）（保存f S J 2和所有f W j 2）
j = 0; )()(02n f n f S =;
while ( j < J ) {
∑∑-=-=++k j
k k
j k k n f S h n f S k n f S g n f W j j j j )2()()
2()(222211
j++;
}
⏹ 逆变换（重构）（利用正变换所保存下来的f S J 2和所有f W j 2）
j = J; while ( j > 0 ) {
∑∑-----+-=-k
j k k j k k n f W g k n f S h n f S j j j )2()2()(121221
j--;
}
)()(02n f S n f =
图10-13 Mallat数据的离散小波变换
说明：
图形的横纵坐标分别表示时间（平移因子）和变换结果S f与W f的值。

小波分解可以无限进行下去，J是自己指定的最大分解次数，一般为8~10。

求和符号中k∈Z，无上下限，但具体计算时，由于只有有限个h k、g k不为0，所以实际上是有限的。

逆变换中h与g上的一杠表示复数的共轭，对于实h与g，则共轭与不共轭相同。

求S f与W f都涉及到对所有的样本求和，不可能只处理一个样本。

3) 小波分解
执行离散小波变换的有效方法是使用滤波器。

该方法是Mallat在1988年开发的，叫做Mallat算法。

这种方法实际上是一种信号的分解方法，在数字信号处理中称为双通道子带编码。

用滤波器执行离散小波变换的概念如图10-14所示。

图10-14 双通道滤波过程
图中，S表示原始的输入信号，通过两个互补的滤波器产生A和D两个信号，A表示信号的近似值(approximations)，D表示信号的细节值(detail)。

在许多应用中，信号的低频部分是最重要的，而高频部分起一个“添加剂”的作用。

犹如声音那样，把高频分量去掉之后，听起来声音确实是变了，但还能够听清楚说的是什么内容。

相反，如果把低频部分去掉，听起
来就莫名其妙。

在小波分析中，近似值是大的缩放因子产生的系数，表示信号的低频分量。

而细节值是小的缩放因子产生的系数，表示信号的高频分量。

由此可见，离散小波变换可以被表示成由低通滤波器和高通滤波器组成的一棵树。

原始信号通过这样的一对滤波器进行的分解叫做一级分解。

信号的分解过程可以叠代，也就是说可进行多级分解。

如果对信号的高频分量不再分解，而对低频分量连续进行分解，就得到许多分辨率较低的低频分量，形成如图10-15所示的一棵比较大的树。

这种树叫做小波分解树(wavelet decomposition tree)。

分解级数的多少取决于要被分析的数据和用户的需要。

小波分解树表示只对信号的低频分量进行连续分解。

图10-15 小波分解树
随便要提及的是，在使用滤波器对真实的数字信号进行变换时，得到的数据将是原始数据的两倍。

例如，如果原始信号的数据样本为1000个，通过滤波之后每一个通道的数据均为1000个，总共为2000个。

于是,根据尼奎斯特(Nyquist)采样定理就提出了降采样(downsampling)的方法，即在每个通道中每两个样本数据取一个，得到的离散小波变换的系数(coefficient)分别用cD和cA表示，如图10-16所示。

图中的符号表示降采样。

图10-16 降采样过程
4) 小波重构
离散小波变换可以用来分析或者叫做分解信号，这个过程叫做分解或者叫做分析。

把分解的系数还原成原始信号的过程叫做小波重构(wavelet reconstruction)或者叫做合成(synthesis)，数学上叫做逆离散小波变换(inverse discrete wavelet transform，IDWT)。

在使用滤波器做小波变换时包含滤波和降采样两个过程，在小波重构时要包含升采样(upsampling)和滤波过程。

小波重构的方法如图10-17所示，图中的符号表示升采样。

图10-17 小波重构方法
升采样是在两个样本数据之间插入“0”，目的是把信号的分量加长。

升采样的过程如图10-18所示。

图10-18 升采样的方法
图10-19是对某数据进行离散小波变换后结果。

图10-19 某数据的离散小波变换。