/*--------------------------------*-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 <olesen@me.queensu.ca> *	- 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 <olesen@me.queensu.ca> *	- 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 <olesen@me.queensu.ca> *	- 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. * * ======================================================================* * 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. * *----------------------------------------------------------------------*/#ifndef _FFTN_C#define _FFTN_C#include <stdio.h>#include <stdlib.h>#include <math.h>#include "fftn.h"/* double precision routine */static intfftradix (double Re[], double Im[],      size_t nTotal, size_t nPass, size_t nSpan, int isign,      int max_factors, int max_perm);/* float precision routine */static intfftradixf (float Re[], float Im[],       size_t nTotal, size_t nPass, size_t nSpan, int isign,       int max_factors, int max_perm);/* parameters for memory management */static size_t SpaceAlloced = 0;static size_t MaxPermAlloced = 0;/* temp space, (void *) since both float and double routines use it */static void *Tmp0 = NULL;	/* temp space for real part */static void *Tmp1 = NULL;	/* temp space for imaginary part */static void *Tmp2 = NULL;	/* temp space for Cosine values */static void *Tmp3 = NULL;	/* temp space for Sine values */static int  *Perm = NULL;	/* Permutation vector */#define NFACTOR	11static int factor [NFACTOR];voidfft_free (void){   SpaceAlloced = MaxPermAlloced = 0;   if (Tmp0 != NULL)	{ free (Tmp0);	Tmp0 = NULL; }   if (Tmp1 != NULL)	{ free (Tmp1);	Tmp1 = NULL; }   if (Tmp2 != NULL)	{ free (Tmp2);	Tmp2 = NULL; }   if (Tmp3 != NULL)	{ free (Tmp3);	Tmp3 = NULL; }   if (Perm != NULL)	{ free (Perm);	Perm = NULL; }}#ifndef M_PI# define M_PI	3.14159265358979323846264338327950288#endif#ifndef SIN60# define SIN60	0.86602540378443865	/* sin(60 deg) */# define COS72	0.30901699437494742	/* cos(72 deg) */# define SIN72	0.95105651629515357	/* sin(72 deg) */#endif/* re-include this source file on the second pass through */#undef REAL#undef FFTN#undef FFTNS#undef FFTRADIX#undef FFTRADIXS#ifndef FFT_NOFLOAT# define REAL		float# define FFTN		fftnf		/* trailing 'f' for float */# define FFTNS		"fftnf"		/* name for error message */# define FFTRADIX	fftradixf	/* trailing 'f' for float */# define FFTRADIXS	"fftradixf"	/* name for error message */# include "fftn.c"			/* include this file again */#endif#undef REAL#undef FFTN#undef FFTNS#undef FFTRADIX#undef FFTRADIXS#ifndef FFT_NODOUBLE# define REAL		double# define FFTN		fftn# define FFTNS		"fftn"# define FFTRADIX	fftradix# define FFTRADIXS	"fftradix"# include "fftn.c"	/* include this file again */#endif#if defined (FFT_NOFLOAT) && defined (FFT_NODOUBLE) && !defined (lint)Error: cannot have both -DFFT_NOFLOAT and -DFFT_NODOUBLE#endif#else	/* _FFTN_C *//* * */intFFTN (int ndim, const int dims[],      REAL Re [],      REAL Im [],      int iSign,      double scaling){   size_t nSpan, nPass, nTotal;   int ret, i, max_factors, max_perm;   /*    * tally the number of elements in the data array    * and determine the number of dimensions    */   nTotal = 1;   if (ndim && dims [0])     {    for (i = 0; i < ndim; i++)      {         if (dims [i] <= 0)           {          fputs ("Error: " FFTNS "() - dimension error\n", stderr);          fft_free ();	/* free-up memory */          return -1;           }         nTotal *= dims [i];      }     }   else     {    ndim = 0;    for (i = 0; dims [i]; i++)      {         if (dims [i] <= 0)           {          fputs ("Error: " FFTNS "() - dimension error\n", stderr);          fft_free ();	/* free-up memory */          return -1;           }         nTotal *= dims [i];         ndim++;      }     }   /* determine maximum number of factors and permuations */#if 1   /*    * follow John Beale's example, just use the largest dimension and don't    * worry about excess allocation.  May be someone else will do it?    */   max_factors = max_perm = 1;   for (i = 0; i < ndim; i++)     {    nSpan = dims [i];    if (nSpan > max_factors) max_factors = nSpan;    if (nSpan > max_perm) max_perm = nSpan;     }#else   /* use the constants used in the original Fortran code */   max_factors = 23;   max_perm = 209;#endif   /* loop over the dimensions: */   nPass = 1;   for (i = 0; i < ndim; i++)     {    nSpan = dims [i];    nPass *= nSpan;    ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign,            max_factors, max_perm);    /* exit, clean-up already done */    if (ret)      return ret;     }   /* Divide through by the normalizing constant: */   if (scaling && scaling != 1.0)     {    if (iSign < 0) iSign = -iSign;    if (scaling < 0.0)      {         scaling = nTotal;         if (scaling < -1.0)           scaling = sqrt (scaling);      }    scaling = 1.0 / scaling;	/* multiply is often faster */    for (i = 0; i < nTotal; i += iSign)      {         Re [i] *= scaling;         Im [i] *= scaling;      }     }   return 0;}/* * singleton's mixed radix routine * * could move allocation out to fftn(), but leave it here so that it's * possible to make this a standalone function */static intFFTRADIX (REAL Re[],      REAL Im[],      size_t nTotal,      size_t nPass,      size_t nSpan,      int iSign,      int max_factors,      int max_perm){   int ii, mfactor, kspan, ispan, inc;   int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt;   REAL radf;   REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp;   REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp;   REAL *Rtmp = NULL;	/* temp space for real part*/   REAL *Itmp = NULL;	/* temp space for imaginary part */   REAL *Cos = NULL;	/* Cosine values */   REAL *Sin = NULL;	/* Sine values */   REAL s60 = SIN60;		/* sin(60 deg) */   REAL c72 = COS72;		/* cos(72 deg) */   REAL s72 = SIN72;		/* sin(72 deg) */   REAL pi2 = M_PI;		/* use PI first, 2 PI later */   /* gcc complains about k3 being uninitialized, but I can't find out where    * or why ... it looks okay to me.    *    * initialize to make gcc happy    */   k3 = 0;   /* gcc complains about c2, c3, s2,s3 being uninitialized, but they're    * only used for the radix 4 case and only AFTER the (s1 == 0.0) pass    * through the loop at which point they will have been calculated.    *    * initialize to make gcc happy    */   c2 = c3 = s2 = s3 = 0.0;   /* Parameter adjustments, was fortran so fix zero-offset */   Re--;   Im--;   if (nPass < 2)     return 0;   /*  allocate storage */   if (SpaceAlloced < max_factors * sizeof (REAL))     {#ifdef SUN_BROKEN_REALLOC    if (!SpaceAlloced)	/* first time */      {         SpaceAlloced = max_factors * sizeof (REAL);         Tmp0 = (void *) malloc (SpaceAlloced);         Tmp1 = (void *) malloc (SpaceAlloced);         Tmp2 = (void *) malloc (SpaceAlloced);         Tmp3 = (void *) malloc (SpaceAlloced);      }    else      {#endif         SpaceAlloced = max_factors * sizeof (REAL);         Tmp0 = (void *) realloc (Tmp0, SpaceAlloced);         Tmp1 = (void *) realloc (Tmp1, SpaceAlloced);         Tmp2 = (void *) realloc (Tmp2, SpaceAlloced);         Tmp3 = (void *) realloc (Tmp3, SpaceAlloced);#ifdef SUN_BROKEN_REALLOC      }#endif     }   else     {    /* allow full use of alloc'd space */    max_factors = SpaceAlloced / sizeof (REAL);     }   if (MaxPermAlloced < max_perm)     {#ifdef SUN_BROKEN_REALLOC    if (!MaxPermAlloced)	/* first time */      Perm = (int *) malloc (max_perm * sizeof(int));    else#endif      Perm = (int *) realloc (Perm, max_perm * sizeof(int));    MaxPermAlloced = max_perm;     }   else     {    /* allow full use of alloc'd space */    max_perm = MaxPermAlloced;     }   if (Tmp0 == NULL || Tmp1 == NULL || Tmp2 == NULL || Tmp3 == NULL       || Perm == NULL)     goto Memory_Error_Label;   /* assign pointers */   Rtmp = (REAL *) Tmp0;   Itmp = (REAL *) Tmp1;   Cos  = (REAL *) Tmp2;   Sin  = (REAL *) Tmp3;   /*    * Function Body    */   inc = iSign;   if (iSign < 0) {      s72 = -s72;      s60 = -s60;      pi2 = -pi2;      inc = -inc;		/* absolute value */   }   /* adjust for strange increments */   nt = inc * nTotal;   ns = inc * nSpan;   kspan = ns;   nn = nt - inc;   jc = ns / nPass;   radf = pi2 * (double) jc;   pi2 *= 2.0;			/* use 2 PI from here on */   ii = 0;   jf = 0;   /*  determine the factors of n */   mfactor = 0;   k = nPass;   while (k % 16 == 0) {      mfactor++;      factor [mfactor - 1] = 4;      k /= 16;   }   j = 3;   jj = 9;   do {      while (k % jj == 0) {     mfactor++;     factor [mfactor - 1] = j;     k /= jj;      }      j += 2;      jj = j * j;   } while (jj <= k);   if (k <= 4) {      kt = mfactor;      factor [mfactor] = k;      if (k != 1)    mfactor++;   } else {      if (k - (k / 4 << 2) == 0) {     mfactor++;     factor [mfactor - 1] = 2;     k /= 4;      }      kt = mfactor;      j = 2;      do {     if (k % j == 0) {        mfactor++;        factor [mfactor - 1] = j;        k /= j;     }     j = ((j + 1) / 2 << 1) + 1;      } while (j <= k);   }   if (kt) {      j = kt;      do {     mfactor++;     factor [mfactor - 1] = factor [j - 1];     j--;      } while (j);   }   /* test that mfactors is in range */   if (mfactor > NFACTOR)     {    fputs ("Error: " FFTRADIXS "() - exceeded number of factors\n", stderr);    goto Memory_Error_Label;      }   /* compute fourier transform */   for (;;) {      sd = radf / (double) kspan;