From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/ieee/chap1/radix4.f | 488 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 488 insertions(+) create mode 100644 math/ieee/chap1/radix4.f (limited to 'math/ieee/chap1/radix4.f') diff --git a/math/ieee/chap1/radix4.f b/math/ieee/chap1/radix4.f new file mode 100644 index 00000000..cd703305 --- /dev/null +++ b/math/ieee/chap1/radix4.f @@ -0,0 +1,488 @@ +c +c----------------------------------------------------------------------- +c subroutine: radix4 +c computes forward or inverse complex dft via radix-4 fft. +c uses autogen technique to yield time efficient program. +c----------------------------------------------------------------------- +c + subroutine radix4(mm,iflag,jflag) +c +c input: +c mm = power of 4 (i.e., n = 4**mm complex point transform) +c (mm.ge.2 and mm.le.5) +c +c iflag = 1 on first pass for given n +c = 0 on subsequent passes for given n +c +c jflag = -1 for forward transform +c = +1 for inverse transform +c +c input/output: +c a = array of dimensions 2*n with real and imaginary parts +c of dft input/output in odd, even array components. +c +c for optimal time efficiency, common is used to pass arrays. +c this means that dimensions of arrays a, ix, and t can be +c modified to reflect maximum value of n = 4**mm to be used. note +c that array "ix" is also dimensioned in subroutine "rad4sb". +c +c i.e., a( ) ix( ) t( ) +c +c m =2 32 38 27 +c m<=3 128 144 135 +c m<=4 512 658 567 +c m<=5 2048 2996 2295 +c + dimension a(2048),ix(2996),t(2295) + dimension nfac(11),np(209) + common ntypl,kkp,index,ixc + common /aa/a + common /xx/ix +c +c check for mm<2 or mm>5 +c + if(mm.lt.2.or.mm.gt.5)stop +c +c initialize on first pass """""""""""""""""""""""""""""""""""""""" +c + if(iflag.eq.1) go to 9999 +c +c fast fourier transform start #################################### +c +8885 kspan=2*4**mm + if(jflag.eq.1) go to 8887 +c +c conjugate data for forward transform +c + do 8886 j=2,n2,2 +8886 a(j)=-a(j) + go to 8889 +c +c multiply data by n**(-1) if inverse transform +c +8887 do 8888 j=1,n2,2 + a(j)=a(j)*xp +8888 a(j+1)=a(j+1)*xp +8889 i=3 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8),it +c*********************************************************************** +c +c 8 multiply butterfly +c +1 kk=ix(i) +c +11 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp +c + bjp=bkp-bjp +c + a(k2+1)=(akp+bjp-ajp)*c707 + a(k2)=a(k2+1)+bjp*cm141 +c + bkp=bkm+ajm + akp=akm-bjm +c + ac0=(akp+bkp)*c924 + a(k1+1)=ac0+akp*cm541 + a(k1) =ac0+bkp*cm131 +c + bkm=bkm-ajm + akm=akm+bjm +c + ac0=(akm+bkm)*c383 + a(k3+1)=ac0+akm*c541 + a(k3) =ac0+bkm*cm131 +c + i=i+1 + kk=ix(i) + if (kk) 111,111,11 +111 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c 4 multiply butterfly +c +2 kk=ix(i) +c +22 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp + a(k2)=-bkp+bjp + a(k2+1)=akp-ajp +c + bkp=bkm+ajm +c + a(k1+1)=(bkp+akm-bjm)*c707 + a(k1)=a(k1+1)+bkp*cm141 +c + akm=akm+bjm +c + a(k3+1)=(akm+ajm-bkm)*c707 + a(k3)=a(k3+1)+akm*cm141 +c + i=i+1 + kk=ix(i) + if (kk) 222,222,22 +222 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c 8 multiply butterfly +c +3 kk=ix(i) +c +33 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp +c + ajp=akp-ajp +c + a(k2+1)=(ajp+bjp-bkp)*c707 + a(k2)=a(k2+1)+ajp*cm141 +c + bkp=bkm+ajm + akp=akm-bjm +c + ac0=(akp+bkp)*c383 + a(k1+1)=ac0+akp*c541 + a(k1) =ac0+bkp*cm131 +c + bkm=bkm-ajm + akm=akm+bjm +c + ac0=(akm+bkm)*cm924 + a(k3+1)=ac0+akm*c541 + a(k3) =ac0+bkm*c131 +c + i=i+1 + kk=ix(i) + if (kk) 333,333,33 +333 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c general 9 multiply butterfly +c +4 kk=ix(i) +c +44 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp +c + ajp=akp-ajp + bjp=bkp-bjp +c + j=ix(i+1) +c + ac0=(ajp+bjp)*t(j+8) + a(k2+1)=ac0+ajp*t(j+6) + a(k2) =ac0+bjp*t(j+7) +c + bkp=bkm+ajm + akp=akm-bjm +c + ac0=(akp+bkp)*t(j+5) + a(k1+1)=ac0+akp*t(j+3) + a(k1) =ac0+bkp*t(j+4) +c + bkm=bkm-ajm + akm=akm+bjm +c + ac0=(akm+bkm)*t(j+2) + a(k3+1)=ac0+akm*t(j) + a(k3) =ac0+bkm*t(j+1) +c + i=i+2 + kk=ix(i) + if (kk) 444,444,44 +444 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c 0 multiply butterfly +c +5 kk=ix(i) +c +55 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp + a(k2)=akp-ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp + a(k2+1)=bkp-bjp +c + a(k3+1)=bkm-ajm + a(k1+1)=bkm+ajm + a(k3)=akm+bjm + a(k1)=akm-bjm +c + i=i+1 + kk=ix(i) + if (kk) 555,555,55 +555 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c offset reduced +c +6 kspan=kspan/4 + i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c bit reversal (shuffling) +c +7 ip1=ix(i) +77 ip2=ix(i+1) + t1=a(ip2) + a(ip2)=a(ip1) + a(ip1)=t1 + t1=a(ip2+1) + a(ip2+1)=a(ip1+1) + a(ip1+1)=t1 + i=i+2 + ip1=ix(i) + if (ip1) 777,777,77 +777 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +8 if(jflag.eq.1) go to 888 +c +c conjugate output if forward transform +c + do 88 j=2,n2,2 +88 a(j)=-a(j) +888 return +c +c fast fourier transform ends ##################################### +c +c initialization phase starts. done only once +c +9999 ixc=1 + n=4**mm + xp=n + xp=1./xp + ntot=n + n2=n*2 + nspan=n + n1test=n/16 + n2test=n/8 + n3test=(3*n)/16 + nspan4=nspan/4 + ibase=0 + isn=1 + inc=isn + rad=8.0*atan(1.0) + pi=4.*atan(1.0) + c707=sin(pi/4.) + cm141=-2.*c707 + c383=sin(pi/8.) + c924=cos(pi/8.) + cm924=-c924 + c541=c924-c383 + cm541=-c541 + c131=c924+c383 + cm131=-c131 +10 nt=inc*ntot + ks=inc*nspan + kspan=ks + jc=ks/n + radf=rad*float(jc)*.5 + i=0 +c +c determine the factors of n +c all factors must be 4 for this version +c + m=0 + k=n +15 m=m+1 + nfac(m)=4 + k=k/4 +20 if(k-(k/4)*4.eq.0) go to 15 + kt=1 + if(n.ge.256) kt=2 + kspan0=kspan + ntypl=0 +c +100 ndelta=kspan0/kspan + index=0 + sd=radf/float(kspan) + cd=2.0*sin(sd)**2 + sd=sin(sd+sd) + kk=1 + i=i+1 +c +c transform for a factor of 4 +c + kspan=kspan/4 + ix(ixc)=0 + ix(ixc+1)=6 + ixc=ixc+2 +c +410 c1=1.0 + s1=0.0 +420 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan + if(s1.eq.0.0) go to 460 +430 if(kspan.ne.nspan4) go to 431 + t(ibase+5)=-(s1+c1) + t(ibase+6)=c1 + t(ibase+4)=s1-c1 + t(ibase+8)=-(s2+c2) + t(ibase+9)=c2 + t(ibase+7)=s2-c2 + t(ibase+2)=-(s3+c3) + t(ibase+3)=c3 + t(ibase+1)=s3-c3 + ibase=ibase+9 +c +431 kkp=(kk-1)*2 + if(index.ne.n1test) go to 150 + call rad4sb(1) + go to 5035 +150 if(index.ne.n2test) go to 160 + call rad4sb(2) + go to 5035 +160 if(index.ne.n3test) go to 170 + call rad4sb(3) + go to 5035 +170 call rad4sb(4) +5035 kk=k3+kspan + if(kk.le.nt) go to 420 +440 index=index+ndelta + c2=c1-(cd*c1+sd*s1) + s1=(sd*c1-cd*s1)+s1 + c1=c2 + c2=c1*c1-s1*s1 + s2=c1*s1+c1*s1 + c3=c2*c1-s2*s1 + s3=c2*s1+s2*c1 + kk=kk-nt+jc + if(kk.le.kspan) go to 420 + kk=kk-kspan+inc + if(kk.le.jc) go to 410 + if(kspan.eq.jc) go to 800 + go to 100 +460 kkp=(kk-1)*2 + call rad4sb(5) +5050 kk=k3+kspan + if(kk.le.nt) go to 420 + go to 440 +c +800 ix(ixc)=0 + ix(ixc+1)=7 + ixc=ixc+2 +c +c compute parameters to permute the results to normal order +c done in two steps +c permutation for square factors of n +c + np(1)=ks + k=kt+kt+1 + if(m.lt.k) k=k-1 + j=1 + np(k+1)=jc +810 np(j+1)=np(j)/nfac(j) + np(k)=np(k+1)*nfac(j) + j=j+1 + k=k-1 + if(j.lt.k) go to 810 + k3=np(k+1) + kspan=np(2) + kk=jc+1 + k2=kspan+1 + j=1 +c +c permutation for single variate transform +c +820 kkp=(kk-1)*2 + k2p=(k2-1)*2 + ix(ixc)=kkp+1 + ix(ixc+1)=k2p+1 + ixc=ixc+2 + kk=kk+inc + k2=kspan+k2 + if(k2.lt.ks) go to 820 +830 k2=k2-np(j) + j=j+1 + k2=np(j+1)+k2 + if(k2.gt.np(j)) go to 830 + j=1 +840 if(kk.lt.k2) go to 820 + kk=kk+inc + k2=kspan+k2 + if(k2.lt.ks) go to 840 + if(kk.lt.ks) go to 830 + jc=k3 + ix(ixc)=0 + ix(ixc+1)=8 + go to 8885 + end -- cgit