subroutine fft(a,b,ntot,n,nspan,isn) c Multivariate complex fourier transform, computed in place c using mixed-radix fast fourier transform algorithm. c By R. C. Singleton, Stanford Research Institute, Sept. 1968 c Arrays A and B originally hold the real and imaginary c components of the data, and return the real and c imaginary components of the resulting fourier coefficients. c Multivariate data is indexed according to the fortran c array element successor function, without limit c on the number of implied multiple subscripts. c the subroutine is called once for each variate. c the calls for a multivariate transform may be in any order. c NTOT is the total number of complex data values. c N is the dimension of the current variable. c NSPAN/N is the spacing of consecutive data values c while indexing the current variable. c The sign of ISN determines the sign of the complex c exponential, and the magnitude of ISN is normally one. c A tri-variate transform with A(N1,N2,N3), B(N1,N2,N3) c is computed by c call fft(a,b,n1*n2*n3,n1,n1,1) c call fft(a,b,n1*n2*n3,n2,n1*n2,1) c call fft(a,b,n1*n2*n3,n3,n1*n2*n3,1) c For a single-variate transform, c NTOT = N = NSPAN = (number of complex data values), e.g. c call fft(a,b,n,n,n,1) c The data can alternatively be stored in a single complex array C c in standard fortran fashion, i.e. alternating real and imaginary c parts. Then with most fortran compilers, the complex array C can c be equivalenced to a real array A, the magnitude of ISN changed c to two to give correct indexing increment, and A(1) and A(2) used c to pass the initial addresses for the sequences of real and c imaginary values, e.g. c complex c(ntot) c real a(2*ntot) c equivalence (c(1),a(1)) c call fft(a(1),a(2),ntot,n,nspan,2) c Arrays AT(MAXF), CK(MAXF), BT(MAXF), SK(MAXF), and NP(MAXP) c are used for temporary storage. If the available storage c is insufficient, the program is terminated by a stop. c MAXF must be .ge. the maximum prime factor of N. c MAXP must be .gt. the number of prime factors of N. c in addition, if the square-free portion K of N has two or c more prime factors, then MAXP must be .ge. K-1. dimension a(1),b(1) c Array storage in NFAC for a maximum of 15 prime factors of N. c If N has more than one square-free factor, the product of the c square-free factors must be .le. 210 dimension nfac(11),np(209) c Array storage for maximum prime factor of 23 dimension at(23),ck(23),bt(23),sk(23) equivalence (i,ii) c The following two constants should agree with the array dimensions. maxp=209 C C Date: Wed, 9 Aug 1995 09:38:49 -0400 C From: ldm@apollo.numis.nwu.edu maxf=23 C if(n .lt. 2) return inc=isn c72=0.30901699437494742 s72=0.95105651629515357 s120=0.86602540378443865 rad=6.2831853071796 if(isn .ge. 0) go to 10 s72=-s72 s120=-s120 rad=-rad inc=-inc 10 nt=inc*ntot ks=inc*nspan kspan=ks nn=nt-inc jc=ks/n radf=rad*float(jc)*0.5 i=0 jf=0 c Determine the factors of N m=0 k=n go to 20 15 m=m+1 nfac(m)=4 k=k/16 20 if(k-(k/16)*16 .eq. 0) go to 15 j=3 jj=9 go to 30 25 m=m+1 nfac(m)=j k=k/jj 30 if(mod(k,jj) .eq. 0) go to 25 j=j+2 jj=j**2 if(jj .le. k) go to 30 if(k .gt. 4) go to 40 kt=m nfac(m+1)=k if(k .ne. 1) m=m+1 go to 80 40 if(k-(k/4)*4 .ne. 0) go to 50 m=m+1 nfac(m)=2 k=k/4 50 kt=m j=2 60 if(mod(k,j) .ne. 0) go to 70 m=m+1 nfac(m)=j k=k/j 70 j=((j+1)/2)*2+1 if(j .le. k) go to 60 80 if(kt .eq. 0) go to 100 j=kt 90 m=m+1 nfac(m)=nfac(j) j=j-1 if(j .ne. 0) go to 90 c Compute Fourier transform 100 sd=radf/float(kspan) cd=2.0*sin(sd)**2 sd=sin(sd+sd) kk=1 i=i+1 if(nfac(i) .ne. 2) go to 400 c Transform for factor of 2 (including rotation factor) kspan=kspan/2 k1=kspan+2 210 k2=kk+kspan ak=a(k2) bk=b(k2) a(k2)=a(kk)-ak b(k2)=b(kk)-bk a(kk)=a(kk)+ak b(kk)=b(kk)+bk kk=k2+kspan if(kk .le. nn) go to 210 kk=kk-nn if(kk .le. jc) go to 210 if(kk .gt. kspan) go to 800 220 c1=1.0-cd s1=sd 230 k2=kk+kspan ak=a(kk)-a(k2) bk=b(kk)-b(k2) a(kk)=a(kk)+a(k2) b(kk)=b(kk)+b(k2) a(k2)=c1*ak-s1*bk b(k2)=s1*ak+c1*bk kk=k2+kspan if(kk .lt. nt) go to 230 k2=kk-nt c1=-c1 kk=k1-k2 if(kk .gt. k2) go to 230 ak=c1-(cd*c1+sd*s1) s1=(sd*c1-cd*s1)+s1 c1=2.0-(ak**2+s1**2) s1=c1*s1 c1=c1*ak kk=kk+jc if(kk .lt. k2) go to 230 k1=k1+inc+inc kk=(k1-kspan)/2+jc if(kk .le. jc+jc) go to 220 go to 100 c Transform for factor of 3 (optional code) 320 k1=kk+kspan k2=k1+kspan ak=a(kk) bk=b(kk) aj=a(k1)+a(k2) bj=b(k1)+b(k2) a(kk)=ak+aj b(kk)=bk+bj ak=-0.5*aj+ak bk=-0.5*bj+bk aj=(a(k1)-a(k2))*s120 bj=(b(k1)-b(k2))*s120 a(k1)=ak-bj b(k1)=bk+aj a(k2)=ak+bj b(k2)=bk-aj kk=k2+kspan if(kk .lt. nn) go to 320 kk=kk-nn if(kk .le. kspan) go to 320 go to 700 c Transform for factor of 4 400 if(nfac(i) .ne. 4) go to 600 kspnn=kspan kspan=kspan/4 410 c1=1.0 s1=0 420 k1=kk+kspan k2=k1+kspan k3=k2+kspan akp=a(kk)+a(k2) akm=a(kk)-a(k2) ajp=a(k1)+a(k3) ajm=a(k1)-a(k3) a(kk)=akp+ajp ajp=akp-ajp bkp=b(kk)+b(k2) bkm=b(kk)-b(k2) bjp=b(k1)+b(k3) bjm=b(k1)-b(k3) b(kk)=bkp+bjp bjp=bkp-bjp if(isn .lt. 0) go to 450 akp=akm-bjm akm=akm+bjm bkp=bkm+ajm bkm=bkm-ajm if(s1 .eq. 0) go to 460 430 a(k1)=akp*c1-bkp*s1 b(k1)=akp*s1+bkp*c1 a(k2)=ajp*c2-bjp*s2 b(k2)=ajp*s2+bjp*c2 a(k3)=akm*c3-bkm*s3 b(k3)=akm*s3+bkm*c3 kk=k3+kspan if(kk .le. nt) go to 420 440 c2=c1-(cd*c1+sd*s1) s1=(sd*c1-cd*s1)+s1 c1=2.0-(c2**2+s1**2) s1=c1*s1 c1=c1*c2 c2=c1**2-s1**2 s2=2.0*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 450 akp=akm+bjm akm=akm-bjm bkp=bkm-ajm bkm=bkm+ajm if(s1 .ne. 0) go to 430 460 a(k1)=akp b(k1)=bkp a(k2)=ajp b(k2)=bjp a(k3)=akm b(k3)=bkm kk=k3+kspan if(kk .le. nt) go to 420 go to 440 c Transform for factor of 5 (optional code) 510 c2=c72**2-s72**2 s2=2.0*c72*s72 520 k1=kk+kspan k2=k1+kspan k3=k2+kspan k4=k3+kspan akp=a(k1)+a(k4) akm=a(k1)-a(k4) bkp=b(k1)+b(k4) bkm=b(k1)-b(k4) ajp=a(k2)+a(k3) ajm=a(k2)-a(k3) bjp=b(k2)+b(k3) bjm=b(k2)-b(k3) aa=a(kk) bb=b(kk) a(kk)=aa+akp+ajp b(kk)=bb+bkp+bjp ak=akp*c72+ajp*c2+aa bk=bkp*c72+bjp*c2+bb aj=akm*s72+ajm*s2 bj=bkm*s72+bjm*s2 a(k1)=ak-bj a(k4)=ak+bj b(k1)=bk+aj b(k4)=bk-aj ak=akp*c2+ajp*c72+aa bk=bkp*c2+bjp*c72+bb aj=akm*s2-ajm*s72 bj=bkm*s2-bjm*s72 a(k2)=ak-bj a(k3)=ak+bj b(k2)=bk+aj b(k3)=bk-aj kk=k4+kspan if(kk .lt. nn) go to 520 kk=kk-nn if(kk .le. kspan) go to 520 go to 700 c Transform for odd factors 600 k=nfac(i) kspnn=kspan kspan=kspan/k if(k .eq. 3) go to 320 if(k .eq. 5) go to 510 if(k .eq. jf) go to 640 jf=k s1=rad/float(k) c1=cos(s1) s1=sin(s1) if(jf .gt. maxf) go to 998 ck(jf)=1.0 sk(jf)=0.0 j=1 630 ck(j)=ck(k)*c1+sk(k)*s1 sk(j)=ck(k)*s1-sk(k)*c1 k=k-1 ck(k)=ck(j) sk(k)=-sk(j) j=j+1 if(j .lt. k) go to 630 640 k1=kk k2=kk+kspnn aa=a(kk) bb=b(kk) ak=aa bk=bb j=1 k1=k1+kspan 650 k2=k2-kspan j=j+1 at(j)=a(k1)+a(k2) ak=at(j)+ak bt(j)=b(k1)+b(k2) bk=bt(j)+bk j=j+1 at(j)=a(k1)-a(k2) bt(j)=b(k1)-b(k2) k1=k1+kspan if(k1 .lt. k2) go to 650 a(kk)=ak b(kk)=bk k1=kk k2=kk+kspnn j=1 660 k1=k1+kspan k2=k2-kspan jj=j ak=aa bk=bb aj=0.0 bj=0.0 k=1 670 k=k+1 ak=at(k)*ck(jj)+ak bk=bt(k)*ck(jj)+bk k=k+1 aj=at(k)*sk(jj)+aj bj=bt(k)*sk(jj)+bj jj=jj+j if(jj .gt. jf) jj=jj-jf if(k .lt. jf) go to 670 k=jf-j a(k1)=ak-bj b(k1)=bk+aj a(k2)=ak+bj b(k2)=bk-aj j=j+1 if(j .lt. k) go to 660 kk=kk+kspnn if(kk .le. nn) go to 640 kk=kk-nn if(kk .le. kspan) go to 640 c Multiply by rotation factor (except for factors of 2 and 4) 700 if(i .eq. m) go to 800 kk=jc+1 710 c2=1.0-cd s1=sd 720 c1=c2 s2=s1 kk=kk+kspan 730 ak=a(kk) a(kk)=c2*ak-s2*b(kk) b(kk)=s2*ak+c2*b(kk) kk=kk+kspnn if(kk .le. nt) go to 730 ak=s1*s2 s2=s1*c2+c1*s2 c2=c1*c2-ak kk=kk-nt+kspan if(kk .le. kspnn) go to 730 c2=c1-(cd*c1+sd*s1) s1=s1+(sd*c1-cd*s1) c1=2.0-(c2**2+s1**2) s1=c1*s1 c2=c1*c2 kk=kk-kspnn+jc if(kk .le. kspan) go to 720 kk=kk-kspan+jc+inc if(kk .le. jc+jc) go to 710 go to 100 c Permute the results to normal order---done in two stages c permutation for square factors of n 800 np(1)=ks if(kt .eq. 0) go to 890 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 if(n .ne. ntot) go to 850 c permutation for single-variate transform (optional code) 820 ak=a(kk) a(kk)=a(k2) a(k2)=ak bk=b(kk) b(kk)=b(k2) b(k2)=bk 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 go to 890 c permutation for multivariate transform 850 k=kk+jc 860 ak=a(kk) a(kk)=a(k2) a(k2)=ak bk=b(kk) b(kk)=b(k2) b(k2)=bk kk=kk+inc k2=k2+inc if(kk .lt. k) go to 860 kk=kk+ks-jc k2=k2+ks-jc if(kk .lt. nt) go to 850 k2=k2-nt+kspan kk=kk-nt+jc if(k2 .lt. ks) go to 850 870 k2=k2-np(j) j=j+1 k2=np(j+1)+k2 if(k2 .gt. np(j)) go to 870 j=1 880 if(kk .lt. k2) go to 850 kk=kk+jc k2=kspan+k2 if(k2 .lt. ks) go to 880 if(kk .lt. ks) go to 870 jc=k3 890 if(2*kt+1 .ge. m) return kspnn=np(kt+1) c permutation for square-free factors of n j=m-kt nfac(j+1)=1 900 nfac(j)=nfac(j)*nfac(j+1) j=j-1 if(j .ne. kt) go to 900 kt=kt+1 nn=nfac(kt)-1 if(nn .gt. maxp) go to 998 jj=0 j=0 go to 906 902 jj=jj-k2 k2=kk k=k+1 kk=nfac(k) 904 jj=kk+jj if(jj .ge. k2) go to 902 np(j)=jj 906 k2=nfac(kt) k=kt+1 kk=nfac(k) j=j+1 if(j .le. nn) go to 904 c Determine the permutation cycles of length greater than 1 j=0 go to 914 910 k=kk kk=np(k) np(k)=-kk if(kk .ne. j) go to 910 k3=kk 914 j=j+1 kk=np(j) if(kk .lt. 0) go to 914 if(kk .ne. j) go to 910 np(j)=-j if(j .ne. nn) go to 914 maxf=inc*maxf c Reorder A and B, following the permutation cycles go to 950 924 j=j-1 if(np(j) .lt. 0) go to 924 jj=jc 926 kspan=jj if(jj .gt. maxf) kspan=maxf jj=jj-kspan k=np(j) kk=jc*k+ii+jj k1=kk+kspan k2=0 928 k2=k2+1 at(k2)=a(k1) bt(k2)=b(k1) k1=k1-inc if(k1 .ne. kk) go to 928 932 k1=kk+kspan k2=k1-jc*(k+np(k)) k=-np(k) 936 a(k1)=a(k2) b(k1)=b(k2) k1=k1-inc k2=k2-inc if(k1 .ne. kk) go to 936 kk=k2 if(k .ne. j) go to 932 k1=kk+kspan k2=0 940 k2=k2+1 a(k1)=at(k2) b(k1)=bt(k2) k1=k1-inc if(k1 .ne. kk) go to 940 if(jj .ne. 0) go to 926 if(j .ne. 1) go to 924 950 j=k3+1 nt=nt-kspnn ii=nt-inc+1 if(nt .ge. 0) go to 924 return c Error finish, insufficient array storage 998 isn=0 print 999 stop 999 format("array bounds exceeded within subroutine fft") end