10快速傅氏变换和离散小波变换.DOC

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

2、叶变换是 20 世纪 60 年代是计算复杂性研究的主要里程碑之一，1965 年Cooley 和 Tukey 所研究的计算离散傅里叶变换(Discrete Fourier Test)的快速傅氏变换(FFT)将计算量从 (n2)下降至 (nlogn)，推进了 FFT 更深层、更广法的研究与应用。FFT 算法有很多版本，但大体上可分为两类：迭代法和递归法，本节仅讨论迭代法，递归法可参见文献1、2。1.1.1 串行 FFT 迭代算法假定 a0,a1, ,an-1 为一个有限长的输入序列， b0, b1, ,bn-1为离散傅里叶变换的结果序列，则有： ，其中 W ，实际上，上式可)1,.20(10nkW

3、akbnmn ie2写成矩阵 W 和向量 a 的乘积形式： 1210)(1)1(2)(0)(412101210 nnnnn aawwb 一般的 n 阶矩阵和 n 维向量相乘，计算时间复杂度为 n2，但由于 W 是一种特殊矩阵，故可以降低计算量。FFT 的计算流图如 图 22.1 所示，其串行算法如下：算法 22.1 单处理器上的 FFT 迭代算法输入：a=(a 0,a1, ,an-1)输出：b=(b 0,b1, ,bn-1)Begin(1)for k=0 to n-1 do ck=akend for(2)for h=logn-1 downto 0 do (2.1) p=2h(2.2) q=n/

4、p(2.3) z=wq/2(2.4) for k=0 to n-1 doif (k mod p=k mod2p) then (i)ck = ck + ck +p(ii)ck +p=( ck - ck +p)z k modp /* ck 不用(i)计算的新值 */end ifend forend for(3)for k=1 to n-1 do br(k) = ck /* r(k)为 k 的位反 */end forEnda0a1a2a3a4a5a6a7b0b4b2b6b1b5b3b7000000221230图 1.1 n=4 时的 FFT 蝶式变换图显然, FFT 算法的计算复杂度为 O(nlog

5、n)。1.1.2 并行 FFT 算法设 P 为处理器的个数，一种并行 FFT 实现时是将 n 维向量 a 划分成 p 个连续的 m 维子向量，这里 ，第 i 个子向量中下标为 im, , (i+1)m-1，其元素被分配至第 ipnm/号处理器（i=0,1, , p-1） 。由 图 1.1 可以看到，FFT 算法由 logn 步构成，依次以 2logn-1、2 logn-2、2、1 为下标跨度做蝶式计算，我们称下标跨度为 2h 的计算为第 h 步（h=logn-1 , logn-2, , 1, 0） 。并行计算可分两阶段执行：第一阶段，第 logn-1 步至第 logm步，由于下标跨度 h m，

6、各处理器之间需要通信；第二阶段，第 logm-1 步至第 0 步各处理器之间不需要通信。具体并行算法框架描述如下：算法 22.2 FFT 并行算法输入：a=(a 0,a1, ,an-1)输出：b=(b 0,b1, ,bn-1)Begin对所有处理器 my_rank(my_rank=0, p-1)同时执行如下的算法:(1)for h=logp-1 downto 0 do /* 第一阶段，第 logn-1 步至第 logm 步各处理器之间需要通信*/(1.1) t=2i, ,l=2(i+logm) ,q=n/l , z=wq/2 , j= j+1 ,v=2j /*开始 j=0*/(1.2)if (

7、my_rank mod t)=(my_rank mod 2t) then /*本处理器的数据作为变换的前项数据 */(i) tt= my_rank+p/v(ii)接收由 tt 号处理器发来的数据块，并将接收的数据块存于cmy_rank*m+n/v到 cmy_rank*m+n/v+m之中(iii)for k=0 to m-1 do ck=ck+ck+n/vck+n/v=( ck- ck+n/v)*z(my_rank*m+k) mod lend for(iv)将存于 cmy_rank*m+n/v到 cmy_rank*m+n/v+m之中的数据块发送到 tt 号处理器else /*本处理器的数据作为变

8、换的后项数据*/(v)将存于之中的数据块发送到 my_rank-p/v 号处理器(vi)接收由 my_rank-p/v 号处理器发来的数据块存于 c 中end if end for(2)for i=logm-1 downto 0 do /*第二阶段，第 logm-1 步至第 0 步各处理器之间不需要通信*/(2.1) l=2i ,q=n/l , z=wq/2 (2.2)for k=0 to m-1 do if (k mod l)=(k mod 2l) then ck=ck+ck+lck+l=( ck- ck+l)*z(my_rank*m+k) mod lend if end forend fo

9、rEnd由于各处理器对其局部存储器中的 m 维子向量做变换计算，计算时间为 ；点对nmlog点通信 次，每次通信量为 m，通信时间为 ，因此快速傅里叶变换的plog2 ptwslog)(2并行计算时间为 。ptnTwslog)(2logMPI 源程序请参见章末附录。1.2 离散小波变换 DWT1.2.1 离散小波变换 DWT 及其串行算法先对一维小波变换作一简单介绍。设 f(x)为一维输入信号，记， ，这里 与 分别称为定标函)2()(/kxxjjjk 2)(/kxjjj )(x数与子波函数， 与 为二个正交基函数的集合。记 P0f=f，在第 级上的一jkjk j维离散小波变换 DWT(Dis

10、crete Wavelet Transform)通过正交投影 Pjf 与 Qjf 将 Pj-1f 分解为：kkjkjjjj dcfQPf 1其中： ， ，这里，h( n)与102)(pnjnkjkchc02)(pjnjcgd )12,.0,.,1( jNLg(n)分别为低通与高通权系数，它们由基函数 与 来确定，p 为权系数)xjk(jk的长度。 为信号的输入数据，N 为输入信号的长度，L 为所需的级数。由上式可见，0nC每级一维 DWT 与一维卷积计算很相似。所不同的是：在 DWT 中，输出数据下标增加 1 时，权系数在输入数据的对应点下标增加 2，这称为“间隔取样” 。算法 22.3 一维

11、离散小波变换串行算法输入：c 0=d0(c00, c10, cN-10)h=(h0, h1, hL-1)g=(g0, g1, gL-1)输出：c ij , 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 jnikjijij ni)(kjiji dgdchmo)2(1mo21, end forend for(2.2)j=j+1, n=n/2end whileEnd显然，算法 22.3 的时间

12、复杂度为 O(N*L)。在实际应用中，很多情况下采用紧支集小波（Compactly Supported Wavelets） ，这时相应的尺度系数和小波系数都是有限长度的，不失一般性设尺度系数只有有限个非零值：h1,hN，N 为偶数，同样取小波使其只有有限个非零值：g 1,gN。为简单起见，设尺度系数与小波函数都是实数。对有限长度的输入数据序列： (其余点的值都Mnxc,20看成 0)，它的离散小波变换为 : knZjjkhc21knjjkgd21,0Jj其中 J 为实际中要求分解的步数，最多不超过 log2M，其逆变换为knZkjnZkjjnhchc2211,Jj注意到尺度系数和输入系列都是有

13、限长度的序列，上述和实际上都只有有限项。若完全按照上述公式计算，在经过 J 步分解后，所得到的 J+1 个序列 和1,0,Jjdkjkc的非零项的个数之和一般要大于 M，究竟这个项目增加到了多少？下面来分析一下上述计算过程。j=0 时计算过程为 knkhxc21knMkgd21不难看出， 1kc的非零值范围为： 即有,1,0,N个非零值。 k的非零值范围相同。继续往下分解时，非零项22NMNk出现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加较多。例如，若输入序列长度为 100，N=4 ，则 1kd有 51 项非零， 2kd有 27 项非零， 3kd有 15 项非零， 4kd有9

14、项非零， 5kd有 6 项非零， 6有 4 项非零， 6c有 4 项非零。这样分解到 6 步后得到的序列的非零项个数的总和为 116，超过了输入序列的长度。在数据压缩等应用中，希望总的长度基本不增加，这样可以提高压缩比、减少存储量并减少实现的难度。可以采用稍微改变计算公式的方法，使输出序列的非零项总和基本上和输入序列的非零项数相等，并且可以完全重构。这种方法也相当于把输入序列进行延长（增加非零项） ，因而称为延拓法。只需考虑一步分解的情形，下面考虑第一步分解(j=1)。将输入序列作延拓，若 M 为偶数，直接将其按 M 为周期延拓；若 M 为奇数，首先令 01Mx。然后按 M+1 为周期延拓。作

15、了这种延拓后再按前述公式计算，相应的变换矩阵已不再是 H 和 G，事实上这时的变换矩阵类似于循环矩阵。例如，当 M=8，N =4 时矩阵 H 变为： 2143 432143214321 2143 0000hhhhh当 M=7，N =4 时矩阵 H 变为：143 32143214321 143 0000hhhh从上述的矩阵表示可以看出，两种情况下的矩阵内都有完全相同的行，这说明作了重复计算，因而从矩阵中去掉重复的那一行不会减少任何信息量，也就是说，这时我们可以对矩阵进行截短（即去掉一行） ，使得所得计算结果仍然可以完全恢复原输入信号。当M=8，N=4 时截短后的矩阵为： 432143214321

16、 2143000hhhH当 M=7，N =4 时截短后的矩阵为： 32143214321 143000hhhH这时的矩阵都只有 行。分解过程成为：201CGD向量 C1 和 D1 都只有 个元素。重构过程为：2M110*CH可以完全重构。矩阵 H，G 有等式H*H+G*G=I一般情况下，按上述方式保留矩阵的 行，可以完全恢复原信号。2M这种方法的优点是最后的序列的非 0 元素的个数基本上和输入序列的非 0 元素个数相同，特别是若输入序列长度为 2 的幂，则完全相同，而且可以完全重构输入信号。其代价是得到的变换系数 Dj 中的一些元素已不再是输入序列的离散小波变换系数，对某些应用可能是不适合的，

17、但在数据压缩等应用领域，这种方法是可行的。1.2.2 离散小波变换并行算法下设输入序列长度 N=2t，不失一般性设尺度系数只有有限个非零值：h 0，h L-1，L为偶数，同样取小波使其只有有限个非零值：g 0，g L-1。为简单起见，我们采用的延拓方法计算。即将有限尺度的序列 按周期 N 延长，使他成为无限长度的)1,(,0nxc序列。这时变换公式也称为周期小波变换。变换公式为：10221LnkjkZnjjk chc102nkjjkdg,Jj其中 表示 n+2k 对于模 N/2j 的最小非负剩余。注意这时 和 jkd是周期为 N/2jjc的周期序列。其逆变换为 knZkjnZkjjngchc2

18、21,Jj从变换公式中可以看出，计算输出点 和 ，需要输入序列 在 n=2k 附近的值1jkcjdjc（一般而言，L 远远小于输入序列的长度） 。设处理器台数为 p，将输入数据按块分配给 p 台处理器，处理器 i 得到数据)12/,10(jjnNc，让处理器 i 负责 和 的计算，(,jjjnPii 1jnc )12)(,2(11 jjjn PNiPid则不难看出，处理器 i 基本上只要用到局部数据，只有 L/2 个点的计算要用到处理器 i+1中的数据，这时做一步并行数据发送：将处理器 i+1 中前 L-1 个数据发送给处理器 i，则各处理器的计算不再需要数据交换，关于本算法其它描述可参见文献

19、1。算法 22.4 离散小波变换并行算法输入：h i(i=0, L-1), gi(i=0, L-1), ci0(i=0, N-1) 输出：c ik (i=0, N/2k-1,k0)Begin对所有处理器 my_rank(my_rank=0, p-1)同时执行如下的算法:(1)j=0;(2)while (j#include #include #include #include “mpi.h“#define MAX_PROCESSOR_NUM 12#define MAX_N 50#define PI 3.141592653#define EPS 10E-8#define V_TAG 99#defi

20、ne P_TAG 100#define Q_TAG 101#define R_TAG 102#define S_TAG 103#define S_TAG2 104void evaluate(complex* f, int beginPos, int endPos, const complex* x, complex *y, int leftPos,int rightPos, int totalLength);void shuffle(complex* f, int beginPos,int endPos);void print(const complex* f, int fLength);vo

21、id readDoubleComplex(FILE *f,complex int main(int argc, char *argv)complex pMAX_N, qMAX_N, s2*MAX_N, r2*MAX_N;complex w2*MAX_N;complex temp;int variableNum;MPI_Status status;int rank, size;int i, j, k, n;int wLength;int everageLength;int moreLength;int startPos;int stopPos;FILE *fin;MPI_Init( MPI_Ge

22、t_rank(MPI_COMM_WORLD, MPI_Get_size(MPI_COMM_WORLD, if(rank = 0)fin = fopen(“dataIn.txt“, “r“);if (fin = NULL)puts(“Not find input data file“);puts(“Please create a file “dataIn.txt“);puts(“ “);puts(“2“);puts(“1.0 2“);puts(“2.0 -1“);exit(-1);readDoubleComplex(fin, variableNum);if (variableNum MAX_N)

23、puts(“variableNum out of range!“);exit(-1);for(i = 0; i (cos(i*2*PI/wLength),sin(i*2*PI/wLength);everageLength = wLength / size;moreLength = wLength % size;startPos = moreLength + rank * everageLength;stopPos = startPos + everageLength - 1;if(rank = 0)startPos = 0;stopPos = moreLength+everageLength

24、- 1;evaluate(p, 0, variableNum - 1, w, s,startPos, stopPos, wLength);evaluate(q, 0, variableNum - 1, w, r, startPos, stopPos, wLength);for(i = startPos; i 0)MPI_Send(s+startPos), everageLength,MPI_DOUBLE_COMPLEX, 0, S_TAG, MPI_COMM_WORLD);MPI_Recv(s,wLength, MPI_DOUBLE_COMPLEX,0, S_ TAG2,MPI_ COMM_W

25、ORLD,elsefor(i = 1; i 0)MPI_Send(r+startPos), everageLength,MPI_DOUBLE_COMPLEX,0, R_TAG,MPI_COMM_WORLD);elsefor(i = 1; i size; i +)MPI_Recv(r+moreLength+i*everageLength),everageLength, MPI_DOUBLE_COMPLEX,i, R_TAG,MPI_COMM_WORLD,puts(“nAfter FFT r(t)=p(t)q(t)“);printf(“r(t) = “);print(r, wLength - 1);puts(“);