* ---------------------------------------------------------------- Here is a nice split-radix, N-dimensional, fast-fourier transform by R. C. Singleton (Stanford Research Institute, Sept. 1968) which has been converted into C code (7/26/95 by John Beale). My very casual tests show this code is significantly faster than the Numerical Recipes fourn() routine (25 vs. 36 seconds for a (1024x1024) floating point matrix). Also, this code is freely redistributable! - 7/26/95 jpb * ---------------------------------------------------------------- Started a minor clean-up of the Fortran-66 code to make it closer to Fortran-77 and added comments to help mark-up the logical structure and then re-did the f2c conversion. There is still a Fortran look to some of the code especially with incrementing a counter and then subtracting 1 for an array subscript, but most compilers should optimize around that and it's best not to fiddle with too many things. I've cleaned-up most of the goto/label spaghetti so that the resulting code is actually recognizable C compared to what f2c produced. Added the wrapper function fftn() to remove the hassle of using fftradix() directly. You MUST call fft_free() to free-up allocated memory. Note that the file re-includes itself with different defines so that the both float and double precision version may compiled together. Don't want the double version? define FFT_NODOUBLE Don't want the float version? define FFT_NOFLOAT Don't want odd radix transforms? define FFT_RADIX4 The suffix `f' indicates a float vs. double version. Made fftradix() a local (static) function since fftn() should be used instead. See INSTALL 1 Aug July 1995 Mark Olesen TODO: - someone should take a look at getting `max_factors' and `max_perm' calculated correctly in fftn() * ---------------------------------------------------------------- /*--------------------------------*-C-*---------------------------------* * File: * fftn.c * * Public: * fft_free (); * fftn / fftnf (); * * Private: * fftradix / fftradixf (); * * Descript: * multivariate complex Fourier transform, computed in place * using mixed-radix Fast Fourier Transform algorithm. * * Fortran code by: * RC Singleton, Stanford Research Institute, Sept. 1968 * * translated by f2c (version 19950721). * * Revisions: * 26 July 95 John Beale * - added maxf and maxp as parameters to fftradix() * * 28 July 95 Mark Olesen * - cleaned-up the Fortran 66 goto spaghetti, only 3 labels remain. * * - added fft_free() to provide some measure of control over * allocation/deallocation. * * - added fftn() wrapper for multidimensional FFTs * * - use -DFFT_NOFLOAT or -DFFT_NODOUBLE to avoid compiling that * precision. Note suffix `f' on the function names indicates * float precision. * * - revised documentation * * 31 July 95 Mark Olesen * - added GNU Public License * - more cleanup * - define SUN_BROKEN_REALLOC to use malloc() instead of realloc() * on the first pass through, apparently needed for old libc * - removed #error directive in favour of some code that simply * won't compile (generate an error that way) * * 1 Aug 95 Mark Olesen * - define FFT_RADIX4 to only have radix 2, radix 4 transforms * - made fftradix /fftradixf () static scope, just use fftn() * instead. If you have good ideas about fixing the factors * in fftn() please do so. * * 8 Jan 95 mj olesen * - fixed typo's, including one that broke scaling for scaling by * total number of matrix elements or the square root of same * - removed unnecessary casts from allocations * ======================================================================* * NIST Guide to Available Math Software. * Source for module FFT from package GO. * Retrieved from NETLIB on Wed Jul 5 11:50:07 1995. * ======================================================================* * *-----------------------------------------------------------------------* * * int fftn (int ndim, const int dims[], REAL Re[], REAL Im[], * int iSign, double scaling); * * NDIM = the total number dimensions * DIMS = a vector of array sizes * if NDIM is zero then DIMS must be zero-terminated * * RE and IM hold the real and imaginary components of the data, and return * the resulting real and imaginary Fourier coefficients. Multidimensional * data *must* be allocated contiguously. There is no limit on the number * of dimensions. * * ISIGN = the sign of the complex exponential (ie, forward or inverse FFT) * the magnitude of ISIGN (normally 1) is used to determine the * correct indexing increment (see below). * * SCALING = normalizing constant by which the final result is *divided* * if SCALING == -1, normalize by total dimension of the transform * if SCALING < -1, normalize by the square-root of the total dimension * * example: * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] * * int dims[3] = {n1,n2,n3} * fftn (3, dims, Re, Im, 1, scaling); * *-----------------------------------------------------------------------* * int fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass, * size_t nSpan, int iSign, size_t max_factors, * size_t max_perm); * * RE, IM - see above documentation * * Although there is no limit on the number of dimensions, fftradix() must * be called once for each dimension, but the calls may be in any order. * * NTOTAL = the total number of complex data values * NPASS = the dimension of the current variable * NSPAN/NPASS = the spacing of consecutive data values while indexing the * current variable * ISIGN - see above documentation * * example: * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3] * * fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp); * fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp); * fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp); * * single-variate transform, * NTOTAL = N = NSPAN = (number of complex data values), * * fftradix (Re, Im, n, n, n, 1, maxf, maxp); * * The data can also be stored in a single array with alternating real and * imaginary parts, the magnitude of ISIGN is changed to 2 to give correct * indexing increment, and data [0] and data [1] used to pass the initial * addresses for the sequences of real and imaginary values, * * example: * REAL data [2*NTOTAL]; * fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp); * * for temporary allocation: * * MAX_FACTORS >= the maximum prime factor of NPASS * MAX_PERM >= the number of prime factors of NPASS. In addition, * if the square-free portion K of NPASS has two or more prime * factors, then MAX_PERM >= (K-1) * * storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS * has more than one square-free factor, the product of the square-free * factors must be <= 210 array storage for maximum prime factor of 23 the * following two constants should agree with the array dimensions. * *-----------------------------------------------------------------------* * * void fft_free (void); * * free-up allocated temporary storage after finished all the Fourier * transforms. * *----------------------------------------------------------------------*/