1*46439007SCharles.ForsythFFTs: module{ 2*46439007SCharles.Forsyth PATH: con "/dis/math/ffts.dis"; 3*46439007SCharles.Forsyth 4*46439007SCharles.Forsyth ffts: fn(a,b:array of real, ntot,n,nspan,isn:int); 5*46439007SCharles.Forsyth}; 6*46439007SCharles.Forsyth 7*46439007SCharles.Forsyth# multivariate complex fourier transform, computed in place 8*46439007SCharles.Forsyth# using mixed-radix fast fourier transform algorithm. 9*46439007SCharles.Forsyth# arrays a and b originally hold the real and imaginary 10*46439007SCharles.Forsyth# components of the data, and return the real and 11*46439007SCharles.Forsyth# imaginary components of the resulting fourier coefficients. 12*46439007SCharles.Forsyth# multivariate data is indexed according to the fortran 13*46439007SCharles.Forsyth# array element successor function, without limit 14*46439007SCharles.Forsyth# on the number of implied multiple subscripts. 15*46439007SCharles.Forsyth# the subroutine is called once for each variate. 16*46439007SCharles.Forsyth# the calls for a multivariate transform may be in any order. 17*46439007SCharles.Forsyth# ntot is the total number of complex data values. 18*46439007SCharles.Forsyth# n is the dimension of the current variable. 19*46439007SCharles.Forsyth# nspan/n is the spacing of consecutive data values 20*46439007SCharles.Forsyth# while indexing the current variable. 21*46439007SCharles.Forsyth# the sign of isn determines the sign of the complex 22*46439007SCharles.Forsyth# exponential, and the magnitude of isn is normally one. 23*46439007SCharles.Forsyth# univariate transform: 24*46439007SCharles.Forsyth# ffts(a,b,n,n,n,1) 25*46439007SCharles.Forsyth# trivariate transform with a(n1,n2,n3), b(n1,n2,n3): 26*46439007SCharles.Forsyth# ffts(a,b,n1*n2*n3,n1,n1,1) 27*46439007SCharles.Forsyth# ffts(a,b,n1*n2*n3,n2,n1*n2,1) 28*46439007SCharles.Forsyth# ffts(a,b,n1*n2*n3,n3,n1*n2*n3,1) 29*46439007SCharles.Forsyth# the data can alternatively be stored in a single vector c 30*46439007SCharles.Forsyth# alternating real and imaginary parts. the magnitude of isn changed 31*46439007SCharles.Forsyth# to two to give correct indexing increment, and a[0:] and a[1:] used 32*46439007SCharles.Forsyth# for a and b 33