任意长二维DFT的多项式变换算法及其并行算法

任意长二维DFT的多项式变换算法及其并行算法 任意长二维DFT的多项式变换算法及其并行算法 N二维复序列,其二维DFT定义为X k1,k2=N-1 设x n1,n2为N× n1=0N-1 n2=0 x n1,n2 W n1k1N W n2k2N(k 1,k 2=0,1,…,N-1)(1)其中 W N=e-2iπ/N。当N=2 t时,(1)的FPT算法最为简单,此处从略。当N=p c(p为奇素数,c为正整数),可将(1)式按k 2=l modp,即k 2=p u+l(l=0,1,…,p-1;u=0,1,…,N/p-1)分离为如下p个式子X k1,pu+l=N-1 n1=0N/p-1 u=0 x(l)n1,n2 W n1k1N W(pu+l)n2N(2)(k 1=0,1,…,N-1;u=0,1,… ,N/P-1;l=0,1,…,p-1)其中x(l)n1,n2=p-1 n=0 x n1,n2+nN/p W nlp(3) (n 1=0,1,…,N-1;n 2=0,1,…,N/p-1;l=0,1,…,p-1)(3)式是N 2/p个p点一维DFT,如设M u(p)和A d(p)是p点DFT所需的乘、加法量,那么计算x(l)n1,n2共需N 2 M u(p)/p和N 2 Ad(p)/p次乘法和加法。(1)当l=0时,若改变 (2)式的计算次序,如果会(?)X(*,0)n1,n2=W-n2N x(0)n1,n2 (n 1=0,1,…,N-1;n 2=0,1,…,N/p-1)X(0)n1(z)=N/p-1 n2=0 x(*,0)n1,n2 z n2 (n 1=0,1,…,N-1)X(0)k1(z)=N-1 n1=0 X(0)n1(z)z n1k1 modP e(z) (k 1=0,1,…,N-1)(4) 就有X〈k1(Pu+1)〉,pu?X(0)k1(z)mod(z-W pu+1N)(5)(k 1=0,1,…,N-1;u=0,1,…,N/p-1)由于(P e(z),z,N)构成多项式变换,故(4)式是{X(0)n1(Z)}的多项式变换,可用FPT计算,只需(p-1)NA d(N)/p+(p-1)N 2/p次加法(这里第二项是模运算所需的加法量,A d(N)是N点一维FFT所需的加法量)。如设(4)式的计算结果为 X(0)k1(z)=(p-1)N/p-1 m=0 y(0)k1,m z m,(k 1=0,1,…,N-1)则由(5)式有 X〈k1(pu+1)〉,pu=(p-1)N/p-1 m=0 y(0)k1 ,m W m(pu+1)N=N/p-1 m=0[W mNp-2 j=0 W jp y(0)k1,m+jN/p]W umN/p(6)(k 1=0,1,…,N-1;u=0,1,…,N/p-1)(6)式是以k 1为参数的序列y(0)k1,m=W mNp-2 j=0 W jp y(0)k1m+jN/p的N/p点一维DFT(共N个 ),可用基-pFFT计算。因此由方程组(?)和(6)式计算X R1,pu共需M 1=NM u N p+N 2-N A 1=NA d N p+p-1 p NA d(N)+2p-3 p N 2次乘法和加法(其中M u(k),A d(k)表示k点FFT的乘加量,下同)。(2)l?0,改变 (2)式计算次序,若令(?)X(l)n1(z)=N/p-1 n2=0 x(l)n1,n2 z n2 (n 1=0,1,…,N-1)Xlk1(z)=N-1 n1=0 X(l)n1(z)z n1k1 modP e(z) (k 1=0,1,…,N-1)(7)则有X〈k1(pu+1)〉,pu+l?X(l)k1(z)mod(z-W pu+lN )(8)(k 1=0,1,…,N-1;u=0,1,…,N/p-1,l?0)与(4)式一样,(7)式是多项式序列{X(l)n1(z)}的多项式变换,可用相同方法计算,所需加法也相同,如设(7)式的计算结果为X(l)k1(z)=(p-1)N/p-1 m=0 y(l)k1,m z m,(k 1=0,1 ,…,N-1,l?0)则由(8)式有 X〈k1(pu+l)〉,pu+l=(p-1)N/p-1 m=0 y(l)k1,m W m(pu+l)N=N/p-1 m=0 W mlNp-2 j=0 W jlp y(l)k1,m+jN/p W umN/p(9)(k 1=0,1,…,N-1;u=0,1,…,N/p-1;l?0)(9)式是以k 1为参数的序列y(l)k1,m=W mlNp-2 j=0 W jlp y(l)k1,m+jN/p的N/p点DFT(共N个),可用基-pFFT计算。因此用方程组(?)和(9)式计算X k1,p u+l(l?0)的运算量是 M 2=NM u N p+p-1 p N 2 A 2=NA d N p+p-1 p NA d(N)+2 p-3 p N 2综上所述,p c×p c二维DFT的FPT算法为:Step1:计算(3)式中N 2/p个p点DFT。Step2:利用方程组(?)和(5)式计算X k1,pu。Step3:利用方程组(?)和(9)式计算X k1,pu+l(l?0)。Step4:结束。上述FPT算法的运算量为M(N)=N 2 p M u(p)+M 1+(p-1)M 2=p NM u N p+M u(p)+p 2-p+1 p N 2-N A(N)=N 2 p A d(p)+A 1+(p-1)A 2=pNA d N p+(p-1)NA d(N) +A d(p)+2p 2-3p p N 2(10)特别,当p=3,N=p s时 ,N×N二维DFT的FPT算法运算量为M(3 s)=3NM u N 3+3N 2-N=2sN 2+2N A(3 S)=3NA d N 3+2NA d(N)+5N 2=(6s+3)N 2(11)式中 用到了M u(3 S)=2s?3 S-3 S+1,A d(3 S)=2s?3 S(s>1),M u(3)=2,A d(3 )=6。当用行列算法计算3 S×3 S二维DFT时,运算量为M 0(3 S)=2?NM u(3 S)=(4s-2)N 2+2N A 0(3 S)=2?N S A d(3 S)=4sN 2(12)比较(11)与(12)式可知,FPT算法的乘法量约减少一半,加法量有所增加(可通过改进算法使加法量减少)。下面给出上述FPT算法的并行算法。Algorithm?:向量计算机的向量并行算法。Step1:用多道FFT计算(3)式中N个p点DFT;所需运算量是Nmu(p)/p次和NA d(p)/p次向量长度为N的并行乘法和加法。Step2:利用方程组(?)和(5)式向量并行计算X k1,pu:方程(?)的第一式需N/p-1次长度为N的并行乘法;用PFPT计算多项式变换需A d(N)+N次长度为(p-1)N/p的并行加法;(5)式用多道FFT计算,共需M u(N/p)+(p-1)N/p次和A d(N/p)+(p-2)N/p次长度为N的并行乘法和加法。Step3:利用方程组(?)和(9)式向量并行计算X k1,pu+l:这一步除Step2中第一步外,其余相同。Step4:结束。Algorithm?所需的运算量为M=p M u N p+[M u(p)+p 2-p+1]N/p-1 A=pA d N p+[A d(p)+p 2-2p ]N/p(13)次向量长度为N的并行乘法和加法,以及A*=pA d(N)+pN次向量长度为(p-1)N/p的并行加法。加速比约等于N。事实上,上述算法还可加以改进,细节从略。现在来给出任意长二维DFT的FDT算法及其并行算法。当N=N 1?N 2,(N 1,N 2)=1,利用简单映射和孙子定理映射,可将N×N二维DFT(1)式映射成(N 1×N 1)×(N 2×N 2)四维DFT,当应用素因子算法计算这个四维DFT时,运算量为M=N 21 M 2+N 22 M 1;A=N 21 A 2+N 22 A 1(14)当用Winograd嵌套算法计算这个四维DFT时,运算量为M=M 1 M 2;A=N 22 A 1+M 1 A 2式中 M I,A I为N I×N I二维DFT的运算量。当N I=2 t或p c时,N I×N I二维DFT如前述可用FPT算法及其并行算法计算,因此N有双因子时,N×N二维DFT可用F PT算法及其并行算法计算。当N=N 1 N 2…N r,(N I,N j)=1(I?j)时,同样用简单映射和孙子定理映射将(1)式映射为(N 1×N 1)×(N 2×N 2)×…×(N r×N r)2r维DFT,当用素因子算法时,运算量为M=r I=1 M I N 2 N 2i,A=r I=1 A I N 2 N 2i(15)当用Winogral d嵌套算法时,运算量为M=r I=1 M I A=A 1 N 2 N 21+M 1 A 2 N 2 N 21 N 22+…+M 1 M 2…M r-1 A r(16)当N I=2 t或p Ci 时,各个N I×N I均可用前述FPT算法及其并行算法计算。 (原载《国防科技大学学报》1995年第1期,有删节) 作者:蒋增荣