离散小波变换

长期以来，离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及概率论等领域中都得到了广泛的应用。

各种快速傅氏变换(FFT)和离散小波变换(DWT)算法不断出现，成为数值代数方面最活跃的一个研究领域，而其意义远远超过了算法研究的范围，进而为诸多科技领域的研究打开了一个崭新的局面。

本章分别对FFT 和DWT 的基本算法作了简单介绍，若需在此方面做进一步研究，可参考文献[2]。

1.1 离散小波变换DWT
1.1.1 离散小波变换DWT 及其串行算法
先对一维小波变换作一简单介绍。

设f (x )为一维输入信号，记)2(2)(2/k x x j j jk -=--φφ，
)2(2)(2/k x x j j jk -=--ψψ，这里)(x φ与)(x ψ分别称为定标函数与子波函数，)}({x jk φ与
)}({x jk ψ为二个正交基函数的集合。

记P 0f =f ，在第j 级上的一维离散小波变换
DWT(Discrete Wavelet Transform)通过正交投影P j f 与Q j f 将P j -1f 分解为：
∑∑+=+=-k
k
jk j k jk j k j j j d c f Q f P f P ψφ1
其中：∑
=-=-+1
1
2)(p n j n k j
k c n h c ，∑=-=-+1
1
2)(p n j n k j k c n g d )12,...,1,0,,...,2,1(-==j N k L j ，这里，{h (n )}与{g (n )}分别为低通与高通权系数，它们由基函数)}({x jk φ与)}({x jk
ψ
来确定，p 为权系数
的长度。

}{0
n C 为信号的输入数据，N 为输入信号的长度，L 为所需的级数。

由上式可见，每级一维DWT 与一维卷积计算很相似。

所不同的是：在DWT 中，输出数据下标增加1时，权系数在输入数据的对应点下标增加2，这称为“间隔取样”。

算法 一维离散小波变换串行算法 输入：c 0
=d 0
(c 00
, c 10
,…, c N-10
)
h=(h 0, h 1,…, h L-1)
g=(g 0, g 1,…, g L-1)
输出：c i j
, d i j
(i=0, 1,…, N/2j-1
, j ≥0)
Begin (1)j=0, n=N (2)While (n ≥1) do for i=0 to n-1 do (2.1.1)c i j+1
=0, d i j+1
=0 (2.1.2)for k=0 to L-1 do
j n i k k j i j i j
n i) (k k j i j i d g d d c h c c mod )2(11mod 211,
+++++++=+=
end for
end for
j=j+1, n=n/2
end while
End
显然，算法的时间复杂度为O(N*L)。

在实际应用中，很多情况下采用紧支集小波（Compactly Supported Wavelets ），这时相应的尺度系数和小波系数都是有限长度的，不失一般性设尺度系数只有有限个非零值：
h 1,…,h N ，N 为偶数，同样取小波使其只有有限个非零值：g 1,…,g N 。

为简单起见，设尺度系
数与小波函数都是实数。

对有限长度的输入数据序列：M n x c n n
,,2,1,0
==(其余点的值都看成0)，它的离散小波变换为:
k n Z
n j n j k h c c 21-∈+∑=
k n Z
n j n j k g c d 21-∈+∑=
1,,1,0-=J j
其中J 为实际中要求分解的步数，最多不超过log 2M ，其逆变换为
k n Z
k j
k k n Z
k j
k j n h c h c c 221
-∑∈-∑∈-+
=
1,, J j =
注意到尺度系数和输入系列都是有限长度的序列，上述和实际上都只有有限项。

若完全
按照上述公式计算，在经过J 步分解后，所得到的J +1个序列1,,1,0,-=J j d j
k 和的非
零项的个数之和一般要大于M ，究竟这个项目增加到了多少？下面来分析一下上述计算过程。

j =0时计算过程为
k n M
n n k g x d 211-=∑=
不难看出，的非零值范围为：,12,,0,1,,12-⎥⎥
⎤
⎢⎢⎡-+-
=M N k 即有⎥⎦
⎤
⎢⎣⎡-+=⎥⎥⎤⎢⎢⎡+--
=21212N M M N k 个非零值。

的非零值范围相同。

继续往下分解时，非零项出现的规律相似。

分解多步后非零项的个数可能比输入序列的长度增加较多。

例如，若输入序列长度为100，N =4，则有51项非零，有27项非零，有15项非零，有9项非零，有6项非零，有4项非零，有4项非零。

这样分解到6步后得到的序列的非零项个数的总和为116，超过了输入序列的长度。

在数据压缩等应用中，希望总的长度基本不增加，这样可以提高压缩比、减少存储量并减少实现的难度。

可以采用稍微改变计算公式的方法，使输出序列的非零项总和基本上和输入序列的非零项数相等，并且可以完全重构。

这种方法也相当于把输入序列进行延长（增加非零项），因而称为延拓法。

只需考虑一步分解的情形，下面考虑第一步分解(j =1)。

将输入序列作延拓，若M 为偶数，直接将其按M 为周期延拓；若M 为奇数，首先令。

然后按M +1为周期延拓。

作了这种延拓后再按前述公式计算，相应的变换矩阵已不再是H 和G ，事实上这时的变换矩阵类似于循环矩阵。

例如，当M =8，N =4时矩阵H 变为：
2
1
4
3
432143214321
21430
0000000000000000h h h h h h h h h h h h h h h h h h h h 当M =7，N =4时矩阵H 变为：
1
4
3
321432143211430
00000000000000h h h h h h h h h h h h h h h h h
从上述的矩阵表示可以看出，两种情况下的矩阵内都有完全相同的行，这说明作了重复计算，因而从矩阵中去掉重复的那一行不会减少任何信息量，也就是说，这时我们可以对矩阵进行截短（即去掉一行），使得所得计算结果仍然可以完全恢复原输入信号。

当M =8，N =4时截短后的矩阵为：
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡=43
2
1
4321432
1
2143
000000000000h h h h h h h h h h h h h h h h H 当M =7，N =4时截短后的矩阵为：
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡=3
2
1
4321432
1
143
0000000000h h h h h h h h h h h h h h H 这时的矩阵都只有⎥⎥
⎤
⎢⎢⎡2M 行。

分解过程成为：
向量C 1 和D 1
都只有⎥⎥
⎤⎢⎢⎡2M 个元素。

重构过程为：
可以完全重构。

矩阵H ，G 有等式
H *H +G *G =I
一般情况下，按上述方式保留矩阵的⎥⎥
⎤
⎢⎢⎡2M 行，可以完全恢复原信号。

这种方法的优点是最后的序列的非0元素的个数基本上和输入序列的非0元素个数相同，特别是若输入序列长度为2的幂，则完全相同，而且可以完全重构输入信号。

其代价是得到的变换系数D j
中的一些元素已不再是输入序列的离散小波变换系数，对某些应用可能是不适合的，但在数据压缩等应用领域，这种方法是可行的。

Begin
对所有处理器my_rank(my_rank=0,…, p -1)同时执行如下的算法:
(1)j =0;
(2)while (j <r ) do
将数据)12/,,1,0(-=j j n N n c 按块分配给p 台处理器 将处理器i +1中前L -1个数据发送给处理器i 处理器i 负责1
+j n c 和)12
)1(,,2(1
11
-+=+++j j j n p N
i p N i
n d 的计算 j =j +1 end while
End
这里每一步分解后数据1+j n c 和1+j n d 已经是按块存储在P 台处理器上，因此算法第一步中的数据分配除了j =0时需要数据传送外，其余各步不需要数据传送（数据已经到位）。

因此，按LogP 模型，算法的总的通信时间为：2(L max(o ,g )+l)，远小于计算时间O(N)。

MPI 源程序请参见所附光盘。

[1]. Verlag,1982
[2]. 陈崚．二维正交子波变换的VLSI 并行计算．电子学报,1995,23(02):95-97。