10快速傅氏变换和离散小波变换 长期以来,快速傅氏变换( Fast fourier transform)和离散小波变换( Discrete Wavelet ansform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及 概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT算法 不断出现,成为数值代数方面最活跃的一个研究领域,而其意义远远超过了算法研究的范围, 进而为诸多科技领域的研究打开了一个崭新的局面。本章分别对FFT和DWT的基本算法作 了简单介绍,若需在此方面做进一步研究,可参考文献(2]。 1.1快速傅里叶变换FFT 离散傅里叶变换是20世纪60年代是计算复杂性研究的主要里程碑之一,1965年 Cooley 和 Tukey所研究的计算高散傅里叶变换 Discrete Fourier Test))的快速傅氏变换(FF1将计算量 从O(m2)下降至O( nlogn,推进了FFT更深层、更广法的研究与应用。FT算法有很多版本, 但大体上可分为两类:迭代法和递归法,本节仅讨论迭代法,递归法可参见文献[1l、[2]。 1.1.1串行FFT选代算法 假定a[0]d1l…anl]为一个有限长的输入序列,b[0],b1l,…b{-1为离散傅里叶 变换的结果序列,则有:研]=∑k(k=012-,n-1),其中Wn=en,实际上,上式可 写成矩阵W和向量a的乘积形式: 10w1w2 =0w2w4 12(m-1) 0(n-)2(n-1)…w(n-Xn-)‖a 般的n阶矩阵和n维向量相乘,计算时间复杂度为n2,但由于W是一种特殊矩阵 故可以降低计算量。FFT的计算流图如图1.1所示,其串行算法如下: 算法221单处理器上的FFI迭代算法 输入:a=(aa,…,an1) 输 (I)for k=0 to n-l Ck-ak (2)for h=logn-I downtoN do (2.3)=nm2
1. 10 快速傅氏变换和离散小波变换 长期以来，快速傅氏变换(Fast Fourier Transform)和离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及 概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT)算法 不断出现，成为数值代数方面最活跃的一个研究领域，而其意义远远超过了算法研究的范围， 进而为诸多科技领域的研究打开了一个崭新的局面。本章分别对 FFT 和 DWT 的基本算法作 了简单介绍，若需在此方面做进一步研究，可参考文献[2]。 1.1 快速傅里叶变换 FFT 离散傅里叶变换是 20 世纪 60 年代是计算复杂性研究的主要里程碑之一，1965 年 Cooley 和 Tukey 所研究的计算离散傅里叶变换(Discrete Fourier Test)的快速傅氏变换(FFT)将计算量 从 О(n 2 )下降至 О(nlogn)，推进了 FFT 更深层、更广法的研究与应用。FFT 算法有很多版本， 但大体上可分为两类：迭代法和递归法，本节仅讨论迭代法，递归法可参见文献[1]、[2]。 1.1.1 串行 FFT 迭代算法 假定 a[0],a[1], …,a[n-1] 为一个有限长的输入序列，b[0], b[1], …,b[n-1]为离散傅里叶 变换的结果序列，则有： [ ] [ ] ( 0,1,2,..., 1) 1 0 = = − − = b k a k W k n n m km n ，其中 W n i n e 2 = ，实际上，上式可 写成矩阵 W 和向量 a 的乘积形式： = − − − − − − − − 1 2 1 0 0 ( 1) 2( 1) ( 1)( 1) 0 2 4 2( 1) 0 1 2 1 0 0 0 0 1 2 1 0 n n n n n n n n a a a a w w w w w w w w w w w w w w w w b b b b 一般的 n 阶矩阵和 n 维向量相乘，计算时间复杂度为 n 2，但由于 W 是一种特殊矩阵， 故可以降低计算量。FFT 的计算流图如图 1.1 所示，其串行算法如下： 算法 22.1 单处理器上的 FFT 迭代算法 输入：a=(a0,a1, …,an-1) 输出：b=(b0,b1, …,bn-1) Begin (1)for k=0 to n-1 do ck=ak end for (2)for h=logn-1 downto 0 do (2.1) p=2h (2.2) q=n/p (2.3) z=w q/2
(2. 4)for k=0 to n-l do if(k mod P=k mod2p) then (Ick=Ck+ Ck+p (i)lx+2=(ck-ck+p)kmo/*c不用()计算的新值* (3)for k-l to n-l do bn(6)=Ck/r(k)为k的位反* end fo 图1.1n=4时的FFT蝶式变换图 显然,FFT算法的计算复杂度为 O(nlogn)。 1.12并行FFT算法 设P为处理器的个数,一种并行FFT实现时是将n维向量a划分成p个连续的m维子 向量,这里m=「n/p,第i个子向量中下标为ⅸ×m,…,(计+1)×m-1,其元素被分配至第i 号处理器(=0,1,…,P1)。由图1.1可以看到,FFT算法由logn步构成,依次以21、 22、…、2、1为下标跨度做蝶式计算,我们称下标跨度为2"的计算为第h步(h=logn-l, logn-2,…,1,0)。并行计算可分两阶段执行:第一阶段,第logn-1步至第logm步,由于下 标跨度h≥m,各处理器之间需要通信:第二阶段,第logm-1步至第0步各处理器之间不 需要通信。具体并行算法框架描述如下: 算法222FFT并行算法 输入:a=(anan,…am1) 输出:b=(bb,…bn-) 对所有处理器 my rank( my rank=0,…,p-1)同时执行如下的算法 (I)for h=logp-I downton do
(2.4) for k=0 to n-1 do if (k mod p=k mod2p) then (i)ck = ck + ck +p (ii)ck +p=( ck - ck +p)z k modp /* ck 不用(i)计算的新值 */ end if end for end for (3)for k=1 to n-1 do br(k) = ck /* r(k)为 k 的位反 */ end for End a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 b 0 b 4 b 2 b 6 b 1 b 5 b 3 b 7 0 0 0 0 0 0 2 2 1 2 3 0 图 1.1 n=4 时的 FFT 蝶式变换图 显然, FFT 算法的计算复杂度为 O(nlogn)。 1.1.2 并行 FFT 算法 设 P 为处理器的个数，一种并行 FFT 实现时是将 n 维向量 a 划分成 p 个连续的 m 维子 向量，这里 m = n / p ，第 i 个子向量中下标为 i×m, …, (i+1)×m-1，其元素被分配至第 i 号处理器（i=0,1, …, p-1）。由图 1.1 可以看到，FFT 算法由 logn 步构成，依次以 2 logn-1、 2 logn-2、…、2、1 为下标跨度做蝶式计算，我们称下标跨度为 2 h 的计算为第 h 步（h=logn-1, logn-2, …, 1, 0）。并行计算可分两阶段执行：第一阶段，第 logn-1 步至第 logm 步，由于下 标跨度 h≥ m，各处理器之间需要通信；第二阶段，第 logm-1 步至第 0 步各处理器之间不 需要通信。具体并行算法框架描述如下： 算法 22.2 FFT 并行算法 输入：a=(a0,a1, …,an-1) 输出：b=(b0,b1, …,bn-1) Begin 对所有处理器 my_rank(my_rank=0,…, p-1)同时执行如下的算法: (1)for h=logp-1 downto 0 do
/*第一阶段,第logn-l步至第logm步各处理器之间需要通信* 1.1)F=2,=2g),q=m,=n92,j=+1,1=2/开始户=0* (.2)if((my rank mod n=(my rank mod 2n)) then /*本处理器的数据作为变换的前项数据* (i)接收由号处理器发来的数据块,并将接收的数据块存于 cmy rank*mm/到 cmy rank*m+ny+m之中 (infor k-=0 to m-I do c{灯]=c[幻]+c[+m c[k+m/=(c灯]-c[k+m)*吗m end for (iv)将存于c[ my rank*m+m小到 cmy rank*m+m+m]之中的数据块发送 到t号处理器 else/*本处理器的数据作为变换的后项数据* (V)将存于之中的数据块发送到 my rank-p/号处理器 ()接收由 my rank-p/v号处理器发来的数据块存于c中 (2) for i=logn- I downto o do/*第二阶段,第logm-1步至第0步各处理器之间 不需要通信* (2.1)=2 (2. 2)for k=0 to m-l do if((k mod /=(k mod 2D)) ck=ck+ck+/ c[+/=(c-c[+门)* my rank.m+h) mod I end if nd f 由于各处理器对其局部存储器中的m维子向量做变换计算,计算时间为mogn;点对 点通信2logp次,每次通信量为m,通信时间为2(;+mr)logp,因此快速傅里叶变换的并 行计算时间为Tp= m logn+2(y+mr) log po MPI源程序请参见章末附录。 12离散小波变换DWT 121离散小波变换DWT及其串行算法 先对一维小波变换作一简单介绍。设(x)为一维输入信号,记中k(x)=2-122-1x-k), (x)=212(2x-k),这里x)与v(x)分别称为定标函数与子波函数,(x)}与
/* 第一阶段，第 logn-1 步至第 logm 步各处理器之间需要通信*/ (1.1) t=2i, ,l=2(i+logm) ,q=n/l , z=w q/2 , j= j+1 ,v=2j /*开始 j=0*/ (1.2)if ((my_rank mod t)=(my_rank mod 2t)) then /*本处理器的数据作为变换的前项数据*/ (i) tt= my_rank+p/v (ii)接收由 tt 号处理器发来的数据块，并将接收的数据块存于 c[my_rank*m+n/v]到 c[my_rank*m+n/v+m]之中 (iii)for k=0 to m-1 do c[k]=c[k]+c[k+n/v] c[k+n/v]=( c[k]- c[k+n/v])*z (my_rank*m+k) mod l end for (iv)将存于 c[my_rank*m+n/v]到 c[my_rank*m+n/v+m]之中的数据块发送 到 tt 号处理器 else /*本处理器的数据作为变换的后项数据*/ (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=w q/2 (2.2)for k=0 to m-1 do if ((k mod l)=(k mod 2l)) then c[k]=c[k]+c[k+l] c[k+l]=( c[k]- c[k+l])*z (my_rank*m+k) mod l end if end for end for End 由于各处理器对其局部存储器中的 m 维子向量做变换计算，计算时间为 mlog n ；点对 点通信 2log p 次，每次通信量为 m，通信时间为 2(ts + mtw)log p ，因此快速傅里叶变换的并 行计算时间为 Tp = mlog n + 2(ts + mt w )log p 。 MPI 源程序请参见章末附录。 1.2 离散小波变换 DWT 1.2.1 离散小波变换 DWT 及其串行算法 先对一维小波变换作一简单介绍。设 f(x)为一维输入信号，记 ( ) 2 (2 ) / 2 x x k j j jk = − − − ， ( ) 2 (2 ) / 2 x x k j j jk = − − − ，这里 (x) 与 (x) 分别称为定标函数与子波函数， { (x)} jk 与
{vA(x)}为二个正交基函数的集合。记Pf千,在第j级上的一维离散小波变换 DWT(Discrete Wavelet transform)通过正交投影Pf与Qf将P分解为: 其中:c=2hnl2k+n,d=2g(nl2kn(=12,-,Lk=0.1-,N/2-1),这里,{h(m) 与{g(m)}分别为低通与高通权系数,它们由基函数钟9(x)与{(x)来确定,P为权系数 的长度。{C}为信号的输入数据,N为输入信号的长度,L为所需的级数。由上式可见 每级一维DWT与一维卷积计算很相似。所不同的是:在DWT中,输出数据下标增加1时 权系数在输入数据的对应点下标增加2,这称为“间隔取样”。 算法223一维高散小波变换串行算法 输入:c=d(co,c1…,cN h1,…,h- g=(80g1,…,.g 输出:c,d(i=0.1,…,MV,j≥0) (1)=0.n=N (2) While(n≥1)do (2. 1.2)for k=0 to L-Ido h k+ 2i)mmd n? d+=d++ k +2i)md n end while 显然,算法22.3的时间复杂度为O∧N*)。 在实际应用中,很多情况下采用紧支集小波( Compactly Supported Wavelets),这时相 应的尺度系数和小波系数都是有限长度的,不失一般性设尺度系数只有有限个非零值 h,…,hw,N为偶数,同样取小波使其只有有限个非零值:g1;…gN。为简单起见,设尺度系 数与小波函数都是实数。对有限长度的输入数据序列:c0=xn,n=12,…,M(其余点的值都 看成0),它的离散小波变换为 =∑cn
{ (x)} jk 为二个正交基函数的集合。记P0f=f，在第 j 级上的一维离散小波变换 DWT(Discrete Wavelet Transform)通过正交投影 Pjf 与 Qjf 将 Pj-1f 分解为： − = + = + k k jk j jk k j Pj 1 f Pj f Qj f ck d 其中： = − = − + 1 0 1 2 ( ) p n j k n j k c h n c ， = − = − + 1 0 1 2 ( ) p n j k n j k d g n c ( = 1,2,..., , = 0,1,..., 2 −1) j j L k N ，这里，{h(n)} 与{g(n)}分别为低通与高通权系数，它们由基函数 { (x)} jk 与 { (x)} jk 来确定，p 为权系数 的长度。 { } 0 Cn 为信号的输入数据，N 为输入信号的长度，L 为所需的级数。由上式可见， 每级一维 DWT 与一维卷积计算很相似。所不同的是：在 DWT 中，输出数据下标增加 1 时， 权系数在输入数据的对应点下标增加 2，这称为“间隔取样”。 算法 22.3 一维离散小波变换串行算法 输入：c 0=d0 (c0 0 , c1 0 ,…, cN-1 0 ) h=(h0, h1,…, hL-1) g=(g0, g1,…, gL-1) 输出：ci j , di j (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)ci j+1=0, di j+1=0 (2.1.2)for k=0 to L-1 do j k k i n j i j i j k ( k i) n j i j ci c h c d d g d( 2 )mod 1 1 2 mod 1 1 , + + + + + + = + = + 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。为简单起见，设尺度系 数与小波函数都是实数。对有限长度的输入数据序列： cn xn ,n 1,2, ,M 0 = = (其余点的值都 看成 0)，它的离散小波变换为: n k n Z j n j ck c h 2 1 − + = n k n Z j n j dk c g 2 1 − + =
其中J为实际中要求分解的步数,最多不超过log2M,其逆变换为 ch-= I ckhm-2k+I ckhm-2k j=J,…l1 注意到尺度系数和输入系列都是有限长度的序列,上述和实际上都只有有限项。若完全 按照上述公式计算,在经过J步分解后,所得到的+1个序列d1,j=01…,J-1和ck的 非零项的个数之和一般要大于M,究竟这个项目增加到了多少?下面来分析一下上述计算 过程。 =0时计算过程为 不难看出 以的非零值范围为:A-+10-121即有 k=-~-1+「M1-「+N-个非零值。d的非零值范围相同。继续往下分解时,非零项出 现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加较多。例如,若输入序 列长度为10,0N=4,则d有51项非零,d有27项非零,d有15项非零,d有9项非 零,d有6项非零,有4项非零,c有4项非零。这样分解到6步后得到的序列的非 零项个数的总和为116,超过了输入序列的长度。在数据压缩等应用中,希望总的长度基本 不增加,这样可以提高压缩比、减少存储量并减少实现的难度 可以采用稍微改变计算公式的方法,使输出序列的非零项总和基本上和输入序列的非零 项数相等,并且可以完全重构。这种方法也相当于把输入序列进行延长(增加非零项),因 而称为延拓法。 只需考虑一步分解的情形,下面考虑第一步分解(=1)。将输入序列作延拓,若M为偶 数,直接将其按M为周期延拓:若M为奇数,首先令xM+1=0。然后按M+l为周期延拓。 作了这种延拓后再按前述公式计算,相应的变换矩阵已不再是H和G,事实上这时的变换 矩阵类似于循环矩阵。例如,当M=8,N4时矩阵H变为: h2h40000Ah2 hh2h3h20000 0000h2h3h h2h40000h2 当M=7,N=4时矩阵H变为:
j = 0,1, , J −1 其中 J 为实际中要求分解的步数，最多不超过 log2M，其逆变换为 n k k Z j n k k k Z j k j cn c h 2 c h 2 1 − − − = + j = J , ,1 注意到尺度系数和输入系列都是有限长度的序列，上述和实际上都只有有限项。若完全 按照上述公式计算，在经过 J 步分解后，所得到的 J+1 个序列 d , j = 0,1, , J −1 j k 和 j k c 的 非零项的个数之和一般要大于 M，究竟这个项目增加到了多少？下面来分析一下上述计算 过程。 j=0 时计算过程为 n k M n ck xnh 2 1 1 − = = n k M n dk xn g 2 1 1 − = = 不 难 看 出 ， 1 k c 的 非 零 值 范 围 为 ： 1, 2 1, , 1,0, , 2 − = − + − N M k 即 有 + − = = − − + 2 1 2 1 2 N M M N k 个非零值。 1 k d 的非零值范围相同。继续往下分解时，非零项出 现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加较多。例如，若输入序 列长度为 100，N=4，则 1 k d 有 51 项非零， 2 k d 有 27 项非零， 3 k d 有 15 项非零， 4 k d 有 9 项非 零， 5 k d 有 6 项非零， 6 k d 有 4 项非零， 6 k c 有 4 项非零。这样分解到 6 步后得到的序列的非 零项个数的总和为 116，超过了输入序列的长度。在数据压缩等应用中，希望总的长度基本 不增加，这样可以提高压缩比、减少存储量并减少实现的难度。 可以采用稍微改变计算公式的方法，使输出序列的非零项总和基本上和输入序列的非零 项数相等，并且可以完全重构。这种方法也相当于把输入序列进行延长（增加非零项），因 而称为延拓法。 只需考虑一步分解的情形，下面考虑第一步分解(j=1)。将输入序列作延拓，若 M 为偶 数，直接将其按 M 为周期延拓；若 M 为奇数，首先令 xM +1 = 0 。然后按 M+1 为周期延拓。 作了这种延拓后再按前述公式计算，相应的变换矩阵已不再是 H 和 G，事实上这时的变换 矩阵类似于循环矩阵。例如，当 M=8，N=4 时矩阵 H 变为： 3 4 1 2 1 2 3 4 1 2 3 4 1 2 3 4 3 4 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h h h h h h h h h h h h h h h h h h h h 当 M=7，N=4 时矩阵 H 变为：
h3h40000h hh2h3h4000 00h1h2h3h0 0000h1h2h 从上述的矩阵表示可以看出,两种情况下的矩阵内都有完全相同的行,这说明作了重复 计算,因而从矩阵中去掉重复的那一行不会减少任何信息量,也就是说,这时我们可以对矩 阵进行截短(即去掉一行),使得所得计算结果仍然可以完全恢复原输入信号。当M=8,N=4 时截短后的矩阵为: h3h200001h2 h h, h, h 0 0 0内h2h3h 当M7,N=4时截短后的矩阵为: 「h3h00004 H=4鸟2鸟000 00h1h2h3h20 000 这时的矩阵都只有 分解过程成为 C= HO D=GCO 向量C和D都只有个元素。重构过程为 H*C+G*D 可以完全重构。矩阵H,G有等式 H*H+G*G=I 一般情况下,按上述方式保留矩阵的行,可以完全恢复原信号。 这种方法的优点是最后的序列的非0元素的个数基本上和输入序列的非0元素个数相 同,特别是若输入序列长度为2的幂,则完全相同,而且可以完全重构输入信号。其代价是 得到的变换系数D中的一些元素已不再是输入序列的离散小波变换系数,对某些应用可能 是不适合的,但在数据压缩等应用领域,这种方法是可行的。 122离散小波变换并行算法 下设输入序列长度N=2,不失一般性设尺度系数只有有限个非零值:ho,…,hL1,L 为偶数,同样取小波使其只有有限个非零值:8o,…,g-1。为简单起见,我们采用的延拓 方法计算。即将有限尺度的序列c0=xn,(n=01…,N-1)按周期N延长,使他成为无限长度的 序列。这时变换公式也称为周期小波变换。变换公式为
3 4 1 1 2 3 1 2 3 4 1 2 3 4 3 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h h h h h h h h h h h h h h h h h 从上述的矩阵表示可以看出，两种情况下的矩阵内都有完全相同的行，这说明作了重复 计算，因而从矩阵中去掉重复的那一行不会减少任何信息量，也就是说，这时我们可以对矩 阵进行截短（即去掉一行），使得所得计算结果仍然可以完全恢复原输入信号。当 M=8，N=4 时截短后的矩阵为： = 1 2 3 4 1 2 3 4 1 2 3 4 3 4 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h h h h h h h h h h h h h h h h H 当 M=7，N=4 时截短后的矩阵为： = 1 2 3 1 2 3 4 1 2 3 4 3 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h h h h h h h h h h h h h h H 这时的矩阵都只有 2 M 行。分解过程成为： 1 0 C = HC 1 0 D = GC 向量 C 1 和 D1 都只有 2 M 个元素。重构过程为： 0 1 1 C = H *C + G * D 可以完全重构。矩阵 H，G 有等式 H*H+G*G=I 一般情况下，按上述方式保留矩阵的 2 M 行，可以完全恢复原信号。 这种方法的优点是最后的序列的非 0 元素的个数基本上和输入序列的非 0 元素个数相 同，特别是若输入序列长度为 2 的幂，则完全相同，而且可以完全重构输入信号。其代价是 得到的变换系数 D j 中的一些元素已不再是输入序列的离散小波变换系数，对某些应用可能 是不适合的，但在数据压缩等应用领域，这种方法是可行的。 1.2.2 离散小波变换并行算法 下设输入序列长度 N=2t，不失一般性设尺度系数只有有限个非零值：h0，…，hL-1，L 为偶数，同样取小波使其只有有限个非零值：g0，…，gL-1。为简单起见，我们采用的延拓 方法计算。即将有限尺度的序列 ,( 0,1, , 1) cn 0 = xn n = N − 按周期 N 延长，使他成为无限长度的 序列。这时变换公式也称为周期小波变换。变换公式为： − = + − + = = 1 0 2 2 1 L n n k j n k n n Z j n j k c c h h c
j=0,1,…,J-1 其中表示叶+2k对于模M2的最小非负剩余。注意这时c和d是周期为M2的 周期序列。其逆变换为 从变换公式中可以看出,计算输出点c和d,需要输入序列e在m=2k附近的值( 般而言,L远远小于输入序列的长度)设处理器台数为p,将输入数据ch(n=01…,N/2/-1) 按块分配给p台处理器,处理器i得到数据m…(+1B,-,让处理器1负数 和d"(n=1P2…,(+1)p2m-D的计算,则不难看出,处理器;基本上只要用到局部数据, 只有L2个点的计算要用到处理器计1中的数据,这时做一步并行数据发送:将处理器+1 中前L-1个数据发送给处理器i,则各处理器的计算不再需要数据交换,关于本算法其它描 述可参见文献[]l 算法224离散小波变换并行算法 输入:h(i 输出:c(=0;…,M2-1,k>0 对所有处理器 my rank( my rank=0…,p-1)同时执行如下的算法: (1)= (2)while (<r)de (2.1)将数据c(n=0,…,N/2J-)按块分配给p台处理器 22)将处理器计+1中前L-1个数据发送给处理器i (23)处理器/负责cn和d(=、 ,…,(i+1) 1)的计算 end while 这里每一步分解后数据cn+和d〃艹已经是按块存储在P台处理器上,因此算法第一步 中的数据分配除了产=0时需要数据传送外,其余各步不需要数据传送(数据已经到位)。因 此,按LogP模型,算法的总的通信时间为:2(Lmax(o.g)+1),远小于计算时间ON。 MPI源程序请参见所附光盘
− = + + = 1 0 2 1 L n n k j n j dk g d j = 0,1, , J −1 其中表示 n+2k 对于模 N/2j 的最小非负剩余。注意这时 j k c 和 j dk 是周期为 N/2j 的 周期序列。其逆变换为 n k k Z j n k k k Z j k j cn c h 2 c g 2 1 − − − = + j = J , ,1 从变换公式中可以看出，计算输出点 j+1 k c 和 j+1 k d ，需要输入序列 j n c 在 n=2k 附近的值（一 般而言，L 远远小于输入序列的长度）。设处理器台数为 p，将输入数据 ( = 0,1, , / 2 −1) j j cn n N 按块分配给 p 台处理器，处理器 i 得到数据 1) 2 , ,( 1) 2 ( = + − j j j n P N i P N c n i ，让处理器 i 负责 j+1 n c 和 1) 2 , ,( 1) 2 ( 1 1 1 = + − + + + j j j n P N i P N d n i 的计算，则不难看出，处理器 i 基本上只要用到局部数据， 只有 L/2 个点的计算要用到处理器 i+1 中的数据，这时做一步并行数据发送：将处理器 i+1 中前 L-1 个数据发送给处理器 i，则各处理器的计算不再需要数据交换，关于本算法其它描 述可参见文献[1]。 算法 22.4 离散小波变换并行算法 输入：hi(i=0,…, L-1), gi(i=0,…, L-1), ci 0 (i=0,…, N-1) 输出：ci k (i=0,…, N/2 k -1,k>0) Begin 对所有处理器 my_rank(my_rank=0,…, p-1)同时执行如下的算法: (1)j=0; (2)while (j<r) do (2.1)将数据 ( = 0,1, , / 2 −1) j j cn n N 按块分配给 p 台处理器 (2.2)将处理器 i+1 中前 L-1 个数据发送给处理器 i (2.3)处理器 i 负责 j+1 n c 和 1) 2 , ,( 1) 2 ( 1 1 1 = + − + + + j j j n p N i p N d n i 的计算 (2.4)j=j+1 end while End 这里每一步分解后数据 j+1 n c 和 j+1 dn 已经是按块存储在 P 台处理器上，因此算法第一步 中的数据分配除了 j=0 时需要数据传送外，其余各步不需要数据传送（数据已经到位）。因 此，按 LogP 模型，算法的总的通信时间为：2(Lmax(o,g)+l)，远小于计算时间 O(N)。 MPI 源程序请参见所附光盘
13小结 本章主要讨论一维FFT和DWT的简单串、并行算法,二维FFT和DWT在光学、地震 以及图象信号处理等方面起着重要的作用。限于篇幅,此处不再予以讨论。同样,FFT和 DWT的并行算法的更为详尽描述可参见文献[2]和文献[3],专门介绍快速傅氏变换和卷积算 法的著作可参见[4]。另外,二维小波变换的并行计算及相关算法可参见文献[5],LogP模型 可参考文献[3] 参考文献 [1]王能超著.数值算法设计.华中理工大学出版社,19889 [2].陈国良编著.并行计算——结构·算法·编程.高等教育出版社,199910 [3].陈国良编著.并行算法设计与分析(修订版).高等教育出版社,200.11 [4]. Nussbaumer H J Fast Fourier Transform and Convolution Algorithms.2nded. Springer Verlag. 1982 [S].陈崚.二维正交子波变换的LSI并行计算.电子学报,1995,23(02):95-97 附录FFT并行算法的MP源程序 1.源程序ftc #include #include *f, int begi #include"mpi.h void print( const complex+f, #define MAX PROCESSOR NUM 12 #define MAX N void read Double Complex(FILE *f, #define Pl 3.141592653 complex &z) #define V TAg int main(int argc, char *argv) Q TAG complex p[MAX NI, q[MAX NI R TAG #define S TAG complex w[2"MAX N] #define S TAG2 complex temp void evaluate(complex*x, int rank size. Int rightPos, int totallLength); Int w
1.3 小结 本章主要讨论一维 FFT 和 DWT 的简单串、并行算法，二维 FFT 和 DWT 在光学、地震 以及图象信号处理等方面起着重要的作用。限于篇幅，此处不再予以讨论。同样，FFT 和 DWT 的并行算法的更为详尽描述可参见文献[2]和文献[3],专门介绍快速傅氏变换和卷积算 法的著作可参见[4]。另外，二维小波变换的并行计算及相关算法可参见文献[5]，LogP 模型 可参考文献[3]。 参考文献 [1]. 王能超 著．数值算法设计．华中理工大学出版社,1988.9 [2]. 陈国良 编著．并行计算——结构·算法·编程．高等教育出版社,1999.10 [3]. 陈国良 编著．并行算法设计与分析（修订版）．高等教育出版社，2002.11 [4]. Nussbaumer H J. Fast Fourier Transform and Convolution Algorithms.2nded. SpringerVerlag,1982 [5]. 陈崚．二维正交子波变换的 VLSI 并行计算．电子学报,1995,23(02):95-97 附录 FFT 并行算法的 MPI 源程序 1. 源程序 fft.c #include #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 #define P_TAG 100 #define Q_TAG 101 #define R_TAG 102 #define S_TAG 103 #define S_TAG2 104 void 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); void readDoubleComplex(FILE *f, complex &z); int main(int argc, char *argv[]) { complex p[MAX_N], q[MAX_N], s[2*MAX_N], r[2*MAX_N]; complex w[2*MAX_N]; complex temp; int variableNum; MPI_Status status; int rank, size; int i, j, k, n; int wLength;
int everagelength; printf("p(t)=) print(p, variableNum); In Int stop Pos: print(q, variableNum) FILE率fin for(i=1; i< size; 1++) MPI Get rank(MPI COMM WORLD MPI Send(&variableNum, I &rank) MPI INTL V TAG MPI COMM WORLD) MPI Get size(MPI COMM WORLD MPI DOUBLE COMPLEX, 1, P TAG PI COMM WORLD) MPI Send(g, variableNum fin= fopen("dataIn. txt,"T), MPI DOUBLE COMPLEX. if (fin==NULL) TAG MPI COMM WORLD) puts("Not find input data file"); puts("<example for dataIn. txt ") puts("2"); MPI Recv(&variableNum, I, MPI INT,O, puts("1.02") V TAG, MPI COMM WORLD, puts("2.0-1"), &status) exit( 1) MPI Recv(p, variableN MPI DOUBLE COMPLEX.O P TAG. PI COMM WORLD read Double Complex(fin, variableNum); MPI Recv(a, variableNum if((var um l)(variable Num MPI DOUBLE COMPLEX. O MAX N) Q TAG, MPI COMM WORLD, { &status) puts("variableNum out of range! " exit( 1) lEngth =2* variableNum for(i=0, 1< variableNum; 1++) readDouble Complex(fin, p); for(i=0; i< variableNum; 1++) read Complex(fin, q: sin(i 2"PI/lEngth) fclose(fin) puts("Read from data file l'dataIn. txt\"); verageLength lEngth/size;
int everageLength; int moreLength; int startPos; int stopPos; FILE *fin; MPI_Init(&argc, &argv); MPI_Get_rank(MPI_COMM_WORLD, &rank); MPI_Get_size(MPI_COMM_WORLD, &size); 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)) { 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 =lEngth size, S TAG2 MPI COMM WORLD stop Pos=start Pos everageLength-1; if(rank ==0) for(int i=1; i0) s0=so*rowLength*1.0) MPI Send((r+startPos), everageLength if(rank >0) MPI DOUBLE COMPLEX. O R TAG MPI COMM WORLD), everagelength MPI DOUBLE COMPLEX, 0, S TAG, MPI COMM WORLD); MPI DOUBLE COMPLEX. O S TAG2, MPI COMM WORLD MPI Recv((r+morelength+i MPI DOUBLE COMPLEX r(i= l; i< size: 1++) MPI COMM WORLD, &status) MPI Recv((s+moreLength+i* erage Length), everageLength MPI DOUBLE COMPLEX puts("nAfter FFTr(tFp(t)q(t"); MPI COMM WORLD print(r, wLength-1); rinf("Use prossor size=%d\n",size) for(i=1; i<size,i++ MPI Send(s, lEngth MPI DOUBLE COMPLEX, i
moreLength = wLength % size; startPos = moreLength + rank * everageLength; stopPos = startPos + everageLength - 1; if(rank == 0) { startPos = 0; stopPos = moreLength+everageLength - 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_WORLD, &status); } else { for(i = 1; i 0) { MPI_Send((r+startPos), everageLength, MPI_DOUBLE_COMPLEX,0, R_TAG, MPI_COMM_WORLD); } else { for(i = 1; i < size; i ++) { MPI_Recv((r+moreLength+i* everageLength), everageLength, MPI_DOUBLE_COMPLEX, i, R_TAG, MPI_COMM_WORLD, &status); } puts("\nAfter FFT r(t)=p(t)q(t)"); printf("r(t) = "); print(r, wLength - 1); puts(""); printf("Use prossor size = %d\n",size); } MPI_Finalize(); }
