小波变换原理推导以及程序

1. 和离散小波变换

长期以来，快速傅氏变换(Fast Fourier Transform)和离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及 概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT)算法 不断出现， 成为数值代数方面最活跃的一个研究领域， 而其意义远远超过了算法研究的范围， 进而为诸多科技领域的研究打开了一个崭新的局面。 本章分别对 FFT 和 DWT 的基本算法作 了简单介绍，若需在此方面做进一步研究，可参考文献[2]。

1.1 离散小波变换 DWT 1.1.1 离散小波变换 DWT 及其串行算法

先对一维小波变换作一简单介绍。 f(x)为一维输入信号，  jk ( x)  2 j / 2  (2 j x  k ) ， 设 记

 jk ( x)  2  j / 2 (2  j x  k ) ，这里  (x) 与  (x) 分别称为定标函数与子波函数， { jk ( x)} 与

记 在第 j 级上的一维离散小波变换 DWT(Discrete { jk ( x)}为二个正交基函数的集合。 P0f=f， Wavelet Transform)通过正交投影 Pjf 与 Qjf 将 Pj-1f 分解为：

Pj 1 f  Pj f  Q j f   c kj  jk   d kj

k k jk

  其中： c kj   h(n)c 2jk 1 n ， d kj   g (n)c 2jk 1 n   n 0 n 0

p 1

p 1

( j  1,2,..., L, k  0,1,..., N 2 j  1) ，这里，{h(n)}

与{g(n)}分别为低通与高通权系数，它们由基函数 { jk ( x)}与 { jk ( x)} 来确定，p 为权系数 

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

每级一维 DWT 与一维卷积计算很相似。所不同的是：在 DWT 中，输出数据下标增加 1 时， 权系数在输入数据的对应点下标增加 2，这称为“间隔取样” 。 算法 22.3 一维离散小波变换串行算法 输入：c0=d0(c00, c10,…, cN-10) h=(h0, h1,…, hL-1) g=(g0, g1,…, gL-1) 输出：cij , dij (i=0, 1,…, N/2j-1, j≥0) Begin (1)j=0, n=N (2)While (n≥1) do (2.1)for i=0 to n-1 do (2.1.1)cij+1=0, dij+1=0

(2.1.2)for k=0 to L-1 do

j cij 1  cij 1  hk c(k  2i) mod n,

di j 1  di j 1  gk d(jk  2i ) mod n

end for end for (2.2)j=j+1, n=n/2 end while End 显然，算法 22.3 的时间复杂度为 O(N*L)。 在实际应用中，很多情况下采用紧支集小波（Compactly Supported Wavelets） ，这时相 应的尺度系数和小波系数都是有限长度的，不失一般性设尺度系数只有有限个非零值： h1,„,hN，N 为偶数，同样取小波使其只有有限个非零值：g1,„,gN。为简单起见，设尺度系

0 数与小波函数都是实数。对有限长度的输入数据序列： cn  xn , n  1,2,, M (其余点的值都

看成 0)，它的离散小波变换为:

c kj 1   c nj hn  2 k

nZ

d kj 1   c nj g n  2 k

nZ

j 

0,1,  , J  1

其中 J 为实际中要求分解的步数，最多不超过 log2M，其逆变换为

j cn 1   ckj hn2k   ckj hn2k kZ kZ

j  J ,,1

注意到尺度系数和输入系列都是有限长度的序列， 上述和实际上都只有有限项。 若完全 按照上述公式计算，在经过 J 步分解后，所得到的 J+1 个序列 d kj , j  0,1,, J  1 和 c k 的

j

非零项的个数之和一般要大于 M，究竟这个项目增加到了多少？下面来分析一下上述计算 过程。 j=0 时计算过程为

c1   x n hn2k k

n 1

M 1 d k   xn g n2k n 1

M

1 N M 不 难 看 出 ， c k 的 非 零 值 范 围 为 ： k    1,,1,0,,    1, 即 有 2 2





k 

N 1  M   M  N  1 1       个非零值。 d k 的非零值范围相同。继续往下分解时，非零项出 2 2 2  

现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加较多。例如，若输入序

1 2 3 4 列长度为 100，N=4，则 d k 有 51 项非零， d k 有 27 项非零， d k 有 15 项非零， d k 有 9 项非

5 6 零， d k 有 6 项非零， d k 有 4 项非零， c k 有 4 项非零。这样分解到 6 步后得到的序列的非

6

零项个数的总和为 116，超过了输入序列的长度。在数据压缩等应用中，希望总的长度基本 不增加，这样可以提高压缩比、减少存储量并减少实现的难度。 可以采用稍微改变计算公式的方法， 使输出序列的非零项总和基本上和输入序列的非零 项数相等，并且可以完全重构。这种方法也相当于把输入序列进行延长（增加非零项） ，因 而称为延拓法。 只需考虑一步分解的情形，下面考虑第一步分解(j=1)。将输入序列作延拓，若 M 为偶 数，直接将其按 M 为周期延拓；若 M 为奇数，首先令 x M 1  0 。然后按 M+1 为周期延拓。 作了这种延拓后再按前述公式计算，相应的变换矩阵已不再是 H 和 G，事实上这时的变换 矩阵类似于循环矩阵。例如，当 M=8，N=4 时矩阵 H 变为：

h3 h1 0 0 h3 h4 h2 0 0 h4 0 h3 h1 0 0 0 h4 h2 0 0 0 0 h3 h1 0 0 0 h4 h2 0 h1 0 0 h3 h1 h2 0 0 h4 h2

当 M=7，N=4 时矩阵 H 变为：

h3 h1 0 0 h3 h4 h2 0 0 h4 0 h3 h1 0 0 0 h4 h2 0 0 0 0 h3 h1 0 0 0 h4 h2 0 h1 0 0 h3 h1

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

h3  h H  1 0  0  h4 h2 0 0 0 h3 h1 0 0 h4 h2 0 0 0 h3 h1 0 0 h4 h2 h1 0 0 h3 h2   0 0  h4  

当 M=7，N=4

时截短后的矩阵为：

h3  h H  1 0  0  h4 h2 0 0 0 h3 h1 0 0 h4 h2 0 0 0 h3 h1 0 0 h4 h2 h1 0 0 h3       

M  这时的矩阵都只有   行。分解过程成为： 2

C 1  HC 0 D1  GC 0

M  向量 C1 和 D1 都只有   个元素。重构过程为： 2

C 0  H * C 1  G * D1

可以完全重构。矩阵 H，G 有等式 H*H+G*G=I

M  一般情况下，按上述方式保留矩阵的   行，可以完全恢复原信号。 2 这种方法的优点是最后的序列的非 0 元素的个数基本上和输入序列的非 0 元素个数相 同，特别是若输入序列长度为 2 的幂，则完全相同，而且可以完全重构输入信号。其代价是 得到的变换系数 Dj 中的一些元素已不再是输入序列的离散小波变换系数，对某些应用可能 是不适合的，但在数据压缩等应用领域，这种方法是可行的。

1.1.2 离散小波变换并行算法

下设输入序列长度 N=2t，不失一般性设尺度系数只有有限个非零值：h0，„，hL-1，L 为偶数，同样取小波使其只有有限个非零值：g0，„，gL-1。为简单起见，我们采用的延拓 0 方法计算。即将有限尺度的序列 cn  xn , (n  0,1,, N  1) 按周期 N 延长，使他成为无限长度的 序列。这时变换公式也称为周期小波变换。变换公式为：

ckj 1   cnj hn2k   hn c j

nZ n 0 L 1 n2k

d kj 1   g n d j

n 0

L 1

n 2 k

j  0,1,  , J  1

其中表示 n+2k 对于模 N/2j 的最小非负剩余。注意这时 ckj 和 d kj 是周期为 N/2j 的 周期序列。其逆变换为

c nj 1 

kZ

 ckj hn2k   ckj g n2 k

kZ

j  J ,,1

j 从变换公式中可以看出， 计算输出点 c kj 1 和 d kj 1 ， 需要输入序列 cn 在 n=2k 附近的值 （一

j 般而言， 远远小于输入序列的长度） 设处理器台数为 p， L 。 将输入数据 cn (n  0,1,, N / 2 j  1)

按块分配给 p 台处理器， 处理器 i 得到数据 c nj (n  i 和 d nj 1 (n  i

N N ,  , (i  1)  1) ， 让处理器 i 负责 c nj 1 j P2 P2 j

N N , , (i  1)  1) 的计算，则不难看出，处理器 i 基本上只要用到局部数据， P2 j 1 P2 j 1

只有 L/2 个点的计算要用到处理器 i+1 中的数据，这时做一步并行数据发送：将处理器 i+1 中前 L-1 个数据发送给处理器 i，则各处理器的计算不再需要数据交换，关于本算法其它描 述可参见文献[1]。 算法 22.4 离散小波变换并行算法 输入：hi(i=0,„, L-1), gi(i=0,„, L-1), ci0(i=0,„, N-1) 输出：cik (i=0,„, N/2k-1,k>0)

Begin 对所有处理器 my_rank(my_rank=0,„, p-1)同时执行如下的算法: (1)j=0; (2)while (j中前 L-1 个数据发送给处理器 i (2.3)处理器 i 负责 c n (2.4)j=j+1 end while End 这里每一步分解后数据 c nj 1 和 d nj 1 已经是按块存储在 P 台处理器上，因此算法第一步 中的数据分配除了 j=0 时需要数据传送外，其余各步不需要数据传送（数据已经到位） 。因 此，按 LogP 模型，算法的总的通信时间为：2(Lmax(o,g)+l)，远小于计算时间 O(N)。 MPI 源程序请参见所附光盘。

j 1

和 d n (n  i

j 1

N N ,, (i  1)  1) 的计算 j 1 p2 p 2 j 1

1.2 小结

本章主要讨论一维 FFT 和 DWT 的简单串、并行算法，二维 FFT 和 DWT 在光学、地震 以及图象信号处理等方面起着重要的作用。限于篇幅，此处不再予以讨论。同样，FFT 和 DWT 的并行算法的更为详尽描述可参见文献[2]和文献[3],专门介绍快速傅氏变换和卷积算 法的著作可参见[4]。另外，二维小波变换的并行计算及相关算法可参见文献[5]，LogP 模型 可参考文献[3]。

参考文献

王能超 著．数值算法设计．华中理工大学出版社,1988.9 陈国良 编著．并行计算——结构·算法·编程．高等教育出版社,1999.10 陈国良 编著．并行算法设计与分析（修订版） ．高等教育出版社，2002.11 Nussbaumer H J. Fast Fourier Transform and Convolution Algorithms.2nded. SpringerVerlag,1982 [5]. 陈崚．二维正交子波变换的 VLSI 并行计算．电子学报,1995,23(02):95-97 [1]. [2]. [3]. [4].

附录 FFT 并行算法的 MPI 源程序

1. 源程序 fft.c

#include #include #include #include

#include

puts(nalize(); }

void shuffle(complex* f, int beginPos, int endPos) {

void evaluate(complex* f, int beginPos, int endPos, const complex* x, complex* y, int leftPos, int rightPos, int totalLength) { int i;

complex temp[2*MAX_N]; int i, j;

for(i = beginPos; i

j = beginPos; for(i = beginPos; i

isPrint = true;

for(i = 1; i EPS) {

for(i = beginPos +1; i

if (isPrint) printf(

1. 和离散小波变换

长期以来，快速傅氏变换(Fast Fourier Transform)和离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及 概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT)算法 不断出现， 成为数值代数方面最活跃的一个研究领域， 而其意义远远超过了算法研究的范围， 进而为诸多科技领域的研究打开了一个崭新的局面。 本章分别对 FFT 和 DWT 的基本算法作 了简单介绍，若需在此方面做进一步研究，可参考文献[2]。

1.1 离散小波变换 DWT 1.1.1 离散小波变换 DWT 及其串行算法

先对一维小波变换作一简单介绍。 f(x)为一维输入信号，  jk ( x)  2 j / 2  (2 j x  k ) ， 设 记

 jk ( x)  2  j / 2 (2  j x  k ) ，这里  (x) 与  (x) 分别称为定标函数与子波函数， { jk ( x)} 与

记 在第 j 级上的一维离散小波变换 DWT(Discrete { jk ( x)}为二个正交基函数的集合。 P0f=f， Wavelet Transform)通过正交投影 Pjf 与 Qjf 将 Pj-1f 分解为：

Pj 1 f  Pj f  Q j f   c kj  jk   d kj

k k jk

  其中： c kj   h(n)c 2jk 1 n ， d kj   g (n)c 2jk 1 n   n 0 n 0

p 1

p 1

( j  1,2,..., L, k  0,1,..., N 2 j  1) ，这里，{h(n)}

与{g(n)}分别为低通与高通权系数，它们由基函数 { jk ( x)}与 { jk ( x)} 来确定，p 为权系数 

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

每级一维 DWT 与一维卷积计算很相似。所不同的是：在 DWT 中，输出数据下标增加 1 时， 权系数在输入数据的对应点下标增加 2，这称为“间隔取样” 。 算法 22.3 一维离散小波变换串行算法 输入：c0=d0(c00, c10,…, cN-10) h=(h0, h1,…, hL-1) g=(g0, g1,…, gL-1) 输出：cij , dij (i=0, 1,…, N/2j-1, j≥0) Begin (1)j=0, n=N (2)While (n≥1) do (2.1)for i=0 to n-1 do (2.1.1)cij+1=0, dij+1=0

(2.1.2)for k=0 to L-1 do

j cij 1  cij 1  hk c(k  2i) mod n,

di j 1  di j 1  gk d(jk  2i ) mod n

end for end for (2.2)j=j+1, n=n/2 end while End 显然，算法 22.3 的时间复杂度为 O(N*L)。 在实际应用中，很多情况下采用紧支集小波（Compactly Supported Wavelets） ，这时相 应的尺度系数和小波系数都是有限长度的，不失一般性设尺度系数只有有限个非零值： h1,„,hN，N 为偶数，同样取小波使其只有有限个非零值：g1,„,gN。为简单起见，设尺度系

0 数与小波函数都是实数。对有限长度的输入数据序列： cn  xn , n  1,2,, M (其余点的值都

看成 0)，它的离散小波变换为:

c kj 1   c nj hn  2 k

nZ

d kj 1   c nj g n  2 k

nZ

j 

0,1,  , J  1

其中 J 为实际中要求分解的步数，最多不超过 log2M，其逆变换为

j cn 1   ckj hn2k   ckj hn2k kZ kZ

j  J ,,1

注意到尺度系数和输入系列都是有限长度的序列， 上述和实际上都只有有限项。 若完全 按照上述公式计算，在经过 J 步分解后，所得到的 J+1 个序列 d kj , j  0,1,, J  1 和 c k 的

j

非零项的个数之和一般要大于 M，究竟这个项目增加到了多少？下面来分析一下上述计算 过程。 j=0 时计算过程为

c1   x n hn2k k

n 1

M 1 d k   xn g n2k n 1

M

1 N M 不 难 看 出 ， c k 的 非 零 值 范 围 为 ： k    1,,1,0,,    1, 即 有 2 2





k 

N 1  M   M  N  1 1       个非零值。 d k 的非零值范围相同。继续往下分解时，非零项出 2 2 2  

现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加较多。例如，若输入序

1 2 3 4 列长度为 100，N=4，则 d k 有 51 项非零， d k 有 27 项非零， d k 有 15 项非零， d k 有 9 项非

5 6 零， d k 有 6 项非零， d k 有 4 项非零， c k 有 4 项非零。这样分解到 6 步后得到的序列的非

6

零项个数的总和为 116，超过了输入序列的长度。在数据压缩等应用中，希望总的长度基本 不增加，这样可以提高压缩比、减少存储量并减少实现的难度。 可以采用稍微改变计算公式的方法， 使输出序列的非零项总和基本上和输入序列的非零 项数相等，并且可以完全重构。这种方法也相当于把输入序列进行延长（增加非零项） ，因 而称为延拓法。 只需考虑一步分解的情形，下面考虑第一步分解(j=1)。将输入序列作延拓，若 M 为偶 数，直接将其按 M 为周期延拓；若 M 为奇数，首先令 x M 1  0 。然后按 M+1 为周期延拓。 作了这种延拓后再按前述公式计算，相应的变换矩阵已不再是 H 和 G，事实上这时的变换 矩阵类似于循环矩阵。例如，当 M=8，N=4 时矩阵 H 变为：

h3 h1 0 0 h3 h4 h2 0 0 h4 0 h3 h1 0 0 0 h4 h2 0 0 0 0 h3 h1 0 0 0 h4 h2 0 h1 0 0 h3 h1 h2 0 0 h4 h2

当 M=7，N=4 时矩阵 H 变为：

h3 h1 0 0 h3 h4 h2 0 0 h4 0 h3 h1 0 0 0 h4 h2 0 0 0 0 h3 h1 0 0 0 h4 h2 0 h1 0 0 h3 h1

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

h3  h H  1 0  0  h4 h2 0 0 0 h3 h1 0 0 h4 h2 0 0 0 h3 h1 0 0 h4 h2 h1 0 0 h3 h2   0 0  h4  

当 M=7，N=4

时截短后的矩阵为：

h3  h H  1 0  0  h4 h2 0 0 0 h3 h1 0 0 h4 h2 0 0 0 h3 h1 0 0 h4 h2 h1 0 0 h3       

M  这时的矩阵都只有   行。分解过程成为： 2

C 1  HC 0 D1  GC 0

M  向量 C1 和 D1 都只有   个元素。重构过程为： 2

C 0  H * C 1  G * D1

可以完全重构。矩阵 H，G 有等式 H*H+G*G=I

M  一般情况下，按上述方式保留矩阵的   行，可以完全恢复原信号。 2 这种方法的优点是最后的序列的非 0 元素的个数基本上和输入序列的非 0 元素个数相 同，特别是若输入序列长度为 2 的幂，则完全相同，而且可以完全重构输入信号。其代价是 得到的变换系数 Dj 中的一些元素已不再是输入序列的离散小波变换系数，对某些应用可能 是不适合的，但在数据压缩等应用领域，这种方法是可行的。

1.1.2 离散小波变换并行算法

下设输入序列长度 N=2t，不失一般性设尺度系数只有有限个非零值：h0，„，hL-1，L 为偶数，同样取小波使其只有有限个非零值：g0，„，gL-1。为简单起见，我们采用的延拓 0 方法计算。即将有限尺度的序列 cn  xn , (n  0,1,, N  1) 按周期 N 延长，使他成为无限长度的 序列。这时变换公式也称为周期小波变换。变换公式为：

ckj 1   cnj hn2k   hn c j

nZ n 0 L 1 n2k

d kj 1   g n d j

n 0

L 1

n 2 k

j  0,1,  , J  1

其中表示 n+2k 对于模 N/2j 的最小非负剩余。注意这时 ckj 和 d kj 是周期为 N/2j 的 周期序列。其逆变换为

c nj 1 

kZ

 ckj hn2k   ckj g n2 k

kZ

j  J ,,1

j 从变换公式中可以看出， 计算输出点 c kj 1 和 d kj 1 ， 需要输入序列 cn 在 n=2k 附近的值 （一

j 般而言， 远远小于输入序列的长度） 设处理器台数为 p， L 。 将输入数据 cn (n  0,1,, N / 2 j  1)

按块分配给 p 台处理器， 处理器 i 得到数据 c nj (n  i 和 d nj 1 (n  i

N N ,  , (i  1)  1) ， 让处理器 i 负责 c nj 1 j P2 P2 j

N N , , (i  1)  1) 的计算，则不难看出，处理器 i 基本上只要用到局部数据， P2 j 1 P2 j 1

只有 L/2 个点的计算要用到处理器 i+1 中的数据，这时做一步并行数据发送：将处理器 i+1 中前 L-1 个数据发送给处理器 i，则各处理器的计算不再需要数据交换，关于本算法其它描 述可参见文献[1]。 算法 22.4 离散小波变换并行算法 输入：hi(i=0,„, L-1), gi(i=0,„, L-1), ci0(i=0,„, N-1) 输出：cik (i=0,„, N/2k-1,k>0)

Begin 对所有处理器 my_rank(my_rank=0,„, p-1)同时执行如下的算法: (1)j=0; (2)while (j中前 L-1 个数据发送给处理器 i (2.3)处理器 i 负责 c n (2.4)j=j+1 end while End 这里每一步分解后数据 c nj 1 和 d nj 1 已经是按块存储在 P 台处理器上，因此算法第一步 中的数据分配除了 j=0 时需要数据传送外，其余各步不需要数据传送（数据已经到位） 。因 此，按 LogP 模型，算法的总的通信时间为：2(Lmax(o,g)+l)，远小于计算时间 O(N)。 MPI 源程序请参见所附光盘。

j 1

和 d n (n  i

j 1

N N ,, (i  1)  1) 的计算 j 1 p2 p 2 j 1

1.2 小结

本章主要讨论一维 FFT 和 DWT 的简单串、并行算法，二维 FFT 和 DWT 在光学、地震 以及图象信号处理等方面起着重要的作用。限于篇幅，此处不再予以讨论。同样，FFT 和 DWT 的并行算法的更为详尽描述可参见文献[2]和文献[3],专门介绍快速傅氏变换和卷积算 法的著作可参见[4]。另外，二维小波变换的并行计算及相关算法可参见文献[5]，LogP 模型 可参考文献[3]。

参考文献

王能超 著．数值算法设计．华中理工大学出版社,1988.9 陈国良 编著．并行计算——结构·算法·编程．高等教育出版社,1999.10 陈国良 编著．并行算法设计与分析（修订版） ．高等教育出版社，2002.11 Nussbaumer H J. Fast Fourier Transform and Convolution Algorithms.2nded. SpringerVerlag,1982 [5]. 陈崚．二维正交子波变换的 VLSI 并行计算．电子学报,1995,23(02):95-97 [1]. [2]. [3]. [4].

附录 FFT 并行算法的 MPI 源程序

1. 源程序 fft.c

#include #include #include #include

#include

puts(nalize(); }

void shuffle(complex* f, int beginPos, int endPos) {

void evaluate(complex* f, int beginPos, int endPos, const complex* x, complex* y, int leftPos, int rightPos, int totalLength) { int i;

complex temp[2*MAX_N]; int i, j;

for(i = beginPos; i

j = beginPos; for(i = beginPos; i

isPrint = true;

for(i = 1; i EPS) {

for(i = beginPos +1; i

if (isPrint) printf(