/*--------------------------------*-C-*---------------------------------* * File: * fftn.c * * multivariate complex Fourier transform, computed in place * using mixed-radix Fast Fourier Transform algorithm. * * Fortran code by: * RC Singleton, Stanford Research Institute, Sept. 1968 * NIST Guide to Available Math Software. * Source for module FFT from package GO. * Retrieved from NETLIB on Wed Jul 5 11:50:07 1995. * translated by f2c (version 19950721) and with lots of cleanup * to make it resemble C by: * MJ Olesen, Queen's University at Kingston, 1995-97 */ /*{{{ Copyright: */ /* * Copyright(c)1995,97 Mark Olesen * Queen's Univ at Kingston (Canada) * * Permission to use, copy, modify, and distribute this software for * any purpose without fee is hereby granted, provided that this * entire notice is included in all copies of any software which is * or includes a copy or modification of this software and in all * copies of the supporting documentation for such software. * * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR * IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S * UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY * KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS * FITNESS FOR ANY PARTICULAR PURPOSE. * * All of which is to say that you can do what you like with this * source code provided you don't try to sell it as your own and you * include an unaltered copy of this message (including the * copyright). * * It is also implicitly understood that bug fixes and improvements * should make their way back to the general Internet community so * that everyone benefits. *----------------------------------------------------------------------*/ /*}}}*/ /*{{{ notes: */ /* * Public: * fft_free * fftn / fftnf * (these are documented in the header file) * * Private: * fftradix / fftradixf * * ----------------------------------------------------------------------* * int fftradix (REAL Re[], REAL Im[], unsigned int nTotal, unsigned int nPass, * unsigned int nSpan, int iSign, unsigned int maxFactors, * unsigned int maxPerm); * * 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. * * * 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: * * MAXFACTORS >= the maximum prime factor of NPASS * MAXPERM >= the number of prime factors of NPASS. In addition, * if the square-free portion K of NPASS has two or more prime * factors, then MAXPERM >= (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. * ----------------------------------------------------------------------*/ /*}}}*/ /*{{{ 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 * * 10 Dec 96 mj olesen * - changes defines to compile *without* float support by default, * use -DFFT_FLOAT to enable. * - shifted some variables to local scope (better hints for optimizer) * - added Michael Steffens * Fortran 90 module * - made it simpler to pass dimensions for 1D FFT. * * 23 Feb 97 Mark Olesen * - removed the GNU Public License (see 21 July 1995 entry), * which should make it clear why I have the right to do so. * - Added copyright notice and submitted to netlib * - Moved documentation for the public functions to the header * files where is will always be available. * ----------------------------------------------------------------------*/ /*}}}*/ #ifndef _FFTN_C #define _FFTN_C /* we use CPP to re-include this same file for double/float cases */ //#if !defined (lint) && !defined (__FILE__) //Error: your compiler is sick! define __FILE__ yourself (a string) //eg, something like -D__FILE__=\"fftn.c\" //#endif #include #include #include #include "fftn.h" /*{{{ defines/constants */ #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 /*}}}*/ /*{{{ static 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 11 static int factor [NFACTOR]; /*}}}*/ /*{{{ fft_free() */ void fft_free (void) { SpaceAlloced = MaxPermAlloced = 0; if (Tmp0) { free (Tmp0); Tmp0 = NULL; } if (Tmp1) { free (Tmp1); Tmp1 = NULL; } if (Tmp2) { free (Tmp2); Tmp2 = NULL; } if (Tmp3) { free (Tmp3); Tmp3 = NULL; } if (Perm) { free (Perm); Perm = NULL; } } /*}}}*/ /* return the number of factors */ static int factorize (unsigned int nPass, int * kt) { int nFactor = 0; int j; unsigned int jj; *kt = 0; /* determine the factors of n */ while ((nPass % 16) == 0) /* factors of 4 */ { factor [nFactor++] = 4; nPass /= 16; } j = 3; jj = 9; /* factors of 3, 5, 7, ... */ do { while ((nPass % jj) == 0) { factor [nFactor++] = j; nPass /= jj; } j += 2; jj = j * j; } while (jj <= nPass); if (nPass <= 4) { *kt = nFactor; factor [nFactor] = nPass; if (nPass != 1) nFactor++; } else { if (nPass - (nPass / 4 << 2) == 0) { factor [nFactor++] = 2; nPass /= 4; } *kt = nFactor; j = 2; do { if ((nPass % j) == 0) { factor [nFactor++] = j; nPass /= j; } j = ((j + 1) / 2 << 1) + 1; } while (j <= nPass); } if (*kt) { j = *kt; do factor [nFactor++] = factor [--j]; while (j); } return nFactor; } /* re-include this source file on the second pass through */ /*{{{ defines for re-including double precision */ #ifdef FFT_NODOUBLE # ifndef FFT_FLOAT # define FFT_FLOAT # endif #else # undef REAL # undef FFTN # undef FFTNS # undef FFTRADIX # undef FFTRADIXS /* defines for double */ # define REAL double # define FFTN fftn # define FFTNS "fftn" # define FFTRADIX fftradix # define FFTRADIXS "fftradix" /* double precision routine */ static int fftradix (double Re[], double Im[], unsigned int nTotal, unsigned int nPass, unsigned int nSpan, int isign, unsigned int maxFactors, unsigned int maxPerm); # include __FILE__ /* include this file again */ #endif /*}}}*/ /*{{{ defines for re-including float precision */ #ifdef FFT_FLOAT # undef REAL # undef FFTN # undef FFTNS # undef FFTRADIX # undef FFTRADIXS /* defines for float */ # 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 */ /* float precision routine */ static int fftradixf (float Re[], float Im[], unsigned int nTotal, unsigned int nPass, unsigned int nSpan, int isign, unsigned int maxFactors, unsigned int maxPerm); # include __FILE__ /* include this file again */ #endif /*}}}*/ #else /* _FFTN_C */ /* * use macros to access the Real/Imaginary parts so that it's possible * to substitute different macros if a complex struct is used */ #ifndef Re_Data # define Re_Data(i) Re[i] # define Im_Data(i) Im[i] #endif /* * */ int FFTN (int ndim, const unsigned int dims [], REAL Re [], REAL Im [], int iSign, double scaling) { unsigned int nTotal; unsigned int maxFactors, maxPerm; /* * tally the number of elements in the data array * and determine the number of dimensions */ nTotal = 1; if (ndim) { if (dims != NULL) { int i; /* number of dimensions was specified */ for (i = 0; i < ndim; i++) { if (dims [i] <= 0) goto Dimension_Error; nTotal *= dims [i]; } } else nTotal *= ndim; } else { int i; /* determine # of dimensions from zero-terminated list */ if (dims == NULL) goto Dimension_Error; for (ndim = i = 0; dims [i]; i++) { if (dims [i] <= 0) goto Dimension_Error; 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? */ if (dims != NULL) { int i; for (maxFactors = maxPerm = 1, i = 0; i < ndim; i++) { if (dims [i] > maxFactors) maxFactors = dims [i]; if (dims [i] > maxPerm) maxPerm = dims [i]; } } else { maxFactors = maxPerm = nTotal; } #else /* use the constants used in the original Fortran code */ maxFactors = 23; maxPerm = 209; #endif /* loop over the dimensions: */ if (dims != NULL) { unsigned int nSpan = 1; int i; for (i = 0; i < ndim; i++) { int ret; nSpan *= dims [i]; ret = FFTRADIX (Re, Im, nTotal, dims [i], nSpan, iSign, maxFactors, maxPerm); /* exit, clean-up already done */ if (ret) return ret; } } else { int ret; ret = FFTRADIX (Re, Im, nTotal, nTotal, nTotal, iSign, maxFactors, maxPerm); /* exit, clean-up already done */ if (ret) return ret; } /* Divide through by the normalizing constant: */ if (scaling && scaling != 1.0) { size_t ist; if (iSign < 0) iSign = -iSign; if (scaling < 0.0) scaling = (scaling < -1.0) ? sqrt (nTotal) : nTotal; scaling = 1.0 / scaling; /* multiply is often faster */ for (ist = 0; ist < nTotal; ist += iSign) { Re_Data (ist) *= scaling; Im_Data (ist) *= scaling; } } return 0; Dimension_Error: fprintf (stderr, "Error: " FFTNS "() - dimension error\n"); fft_free (); /* free-up memory */ return -1; } /*----------------------------------------------------------------------*/ /* * 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 int FFTRADIX (REAL Re [], REAL Im [], unsigned int nTotal, unsigned int nPass, unsigned int nSpan, int iSign, unsigned int maxFactors, unsigned int maxPerm) { int ii, nFactor, kspan, ispan, inc; int j, jc, jf, jj, k, k1, k3, kk, kt, nn, ns, nt; REAL radf; REAL c1, c2, c3, cd; REAL s1, s2, s3, sd; 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 */ #ifndef FFT_RADIX4 REAL s60 = SIN60; /* sin(60 deg) */ REAL s72 = SIN72; /* sin(72 deg) */ REAL c72 = COS72; /* cos(72 deg) */ #endif 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 < maxFactors * sizeof (REAL)) { #ifdef SUN_BROKEN_REALLOC if (!SpaceAlloced) /* first time */ { SpaceAlloced = maxFactors * sizeof (REAL); Tmp0 = malloc (SpaceAlloced); Tmp1 = malloc (SpaceAlloced); Tmp2 = malloc (SpaceAlloced); Tmp3 = malloc (SpaceAlloced); } else { #endif SpaceAlloced = maxFactors * sizeof (REAL); Tmp0 = realloc (Tmp0, SpaceAlloced); Tmp1 = realloc (Tmp1, SpaceAlloced); Tmp2 = realloc (Tmp2, SpaceAlloced); Tmp3 = realloc (Tmp3, SpaceAlloced); #ifdef SUN_BROKEN_REALLOC } #endif } else { /* allow full use of alloc'd space */ maxFactors = SpaceAlloced / sizeof (REAL); } if (MaxPermAlloced < (size_t)maxPerm) { #ifdef SUN_BROKEN_REALLOC if (!MaxPermAlloced) /* first time */ Perm = malloc (maxPerm * sizeof(int)); else #endif Perm = realloc (Perm, maxPerm * sizeof(int)); MaxPermAlloced = maxPerm; } else { /* allow full use of alloc'd space */ maxPerm = MaxPermAlloced; } if (!Tmp0 || !Tmp1 || !Tmp2 || !Tmp3 || !Perm) goto Memory_Error; /* assign pointers */ Rtmp = (REAL *) Tmp0; Itmp = (REAL *) Tmp1; Cos = (REAL *) Tmp2; Sin = (REAL *) Tmp3; /* * Function Body */ inc = iSign; if (iSign < 0) { #ifndef FFT_RADIX4 s60 = -s60; s72 = -s72; #endif 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 */ nFactor = factorize (nPass, &kt); /* test that nFactors is in range */ if (nFactor > NFACTOR) { fprintf (stderr, "Error: " FFTRADIXS "() - exceeded number of factors\n"); goto Memory_Error; } /* compute fourier transform */ for (;;) { sd = radf / (double) kspan; cd = sin (sd); cd = 2.0 * cd * cd; sd = sin (sd + sd); kk = 1; ii++; switch (factor [ii - 1]) { case 2: /* transform for factor of 2 (including rotation factor) */ kspan /= 2; k1 = kspan + 2; do { do { REAL tmpr; REAL tmpi; int k2; k2 = kk + kspan; tmpr = Re_Data (k2); tmpi = Im_Data (k2); Re_Data (k2) = Re_Data (kk) - tmpr; Im_Data (k2) = Im_Data (kk) - tmpi; Re_Data (kk) += tmpr; Im_Data (kk) += tmpi; kk = k2 + kspan; } while (kk <= nn); kk -= nn; } while (kk <= jc); if (kk > kspan) goto Permute_Results; /* exit infinite loop */ do { int k2; c1 = 1.0 - cd; s1 = sd; do { REAL tmp; do { do { REAL tmpr; REAL tmpi; k2 = kk + kspan; tmpr = Re_Data (kk) - Re_Data (k2); tmpi = Im_Data (kk) - Im_Data (k2); Re_Data (kk) += Re_Data (k2); Im_Data (kk) += Im_Data (k2); Re_Data (k2) = c1 * tmpr - s1 * tmpi; Im_Data (k2) = s1 * tmpr + c1 * tmpi; kk = k2 + kspan; } while (kk < nt); k2 = kk - nt; c1 = -c1; kk = k1 - k2; } while (kk > k2); tmp = c1 - (cd * c1 + sd * s1); s1 = sd * c1 - cd * s1 + s1; c1 = 2.0 - (tmp * tmp + s1 * s1); s1 *= c1; c1 *= tmp; kk += jc; } while (kk < k2); k1 += (inc + inc); kk = (k1 - kspan) / 2 + jc; } while (kk <= jc + jc); break; case 4: /* transform for factor of 4 */ ispan = kspan; kspan /= 4; do { c1 = 1.0; s1 = 0.0; do { do { REAL ajm, ajp, akm, akp; REAL bjm, bjp, bkm, bkp; int k2; k1 = kk + kspan; k2 = k1 + kspan; k3 = k2 + kspan; akp = Re_Data (kk) + Re_Data (k2); akm = Re_Data (kk) - Re_Data (k2); ajp = Re_Data (k1) + Re_Data (k3); ajm = Re_Data (k1) - Re_Data (k3); bkp = Im_Data (kk) + Im_Data (k2); bkm = Im_Data (kk) - Im_Data (k2); bjp = Im_Data (k1) + Im_Data (k3); bjm = Im_Data (k1) - Im_Data (k3); Re_Data (kk) = akp + ajp; Im_Data (kk) = bkp + bjp; ajp = akp - ajp; bjp = bkp - bjp; if (iSign < 0) { akp = akm + bjm; bkp = bkm - ajm; akm -= bjm; bkm += ajm; } else { akp = akm - bjm; bkp = bkm + ajm; akm += bjm; bkm -= ajm; } /* avoid useless multiplies */ if (s1 == 0.0) { Re_Data (k1) = akp; Re_Data (k2) = ajp; Re_Data (k3) = akm; Im_Data (k1) = bkp; Im_Data (k2) = bjp; Im_Data (k3) = bkm; } else { Re_Data (k1) = akp * c1 - bkp * s1; Re_Data (k2) = ajp * c2 - bjp * s2; Re_Data (k3) = akm * c3 - bkm * s3; Im_Data (k1) = akp * s1 + bkp * c1; Im_Data (k2) = ajp * s2 + bjp * c2; Im_Data (k3) = akm * s3 + bkm * c3; } kk = k3 + kspan; } while (kk <= nt); c2 = c1 - (cd * c1 + sd * s1); s1 = sd * c1 - cd * s1 + s1; c1 = 2.0 - (c2 * c2 + s1 * s1); s1 *= c1; c1 *= c2; /* values of c2, c3, s2, s3 that will get used next time */ c2 = c1 * c1 - s1 * s1; s2 = 2.0 * c1 * s1; c3 = c2 * c1 - s2 * s1; s3 = c2 * s1 + s2 * c1; kk = kk - nt + jc; } while (kk <= kspan); kk = kk - kspan + inc; } while (kk <= jc); if (kspan == jc) goto Permute_Results; /* exit infinite loop */ break; default: /* transform for odd factors */ #ifdef FFT_RADIX4 fprintf (stderr, "Error: " FFTRADIXS "(): compiled for radix 2/4 only\n"); fft_free (); /* free-up memory */ return -1; break; #else /* FFT_RADIX4 */ ispan = kspan; k = factor [ii - 1]; kspan /= factor [ii - 1]; switch (factor [ii - 1]) { case 3: /* transform for factor of 3 (optional code) */ do { do { REAL aj, tmpr; REAL bj, tmpi; int k2; k1 = kk + kspan; k2 = k1 + kspan; tmpr = Re_Data (kk); tmpi = Im_Data (kk); aj = Re_Data (k1) + Re_Data (k2); bj = Im_Data (k1) + Im_Data (k2); Re_Data (kk) = tmpr + aj; Im_Data (kk) = tmpi + bj; tmpr -= 0.5 * aj; tmpi -= 0.5 * bj; aj = (Re_Data (k1) - Re_Data (k2)) * s60; bj = (Im_Data (k1) - Im_Data (k2)) * s60; Re_Data (k1) = tmpr - bj; Re_Data (k2) = tmpr + bj; Im_Data (k1) = tmpi + aj; Im_Data (k2) = tmpi - aj; kk = k2 + kspan; } while (kk < nn); kk -= nn; } while (kk <= kspan); break; case 5: /* transform for factor of 5 (optional code) */ c2 = c72 * c72 - s72 * s72; s2 = 2.0 * c72 * s72; do { do { REAL aa, aj, ak, ajm, ajp, akm, akp; REAL bb, bj, bk, bjm, bjp, bkm, bkp; int k2, k4; k1 = kk + kspan; k2 = k1 + kspan; k3 = k2 + kspan; k4 = k3 + kspan; akp = Re_Data (k1) + Re_Data (k4); akm = Re_Data (k1) - Re_Data (k4); bkp = Im_Data (k1) + Im_Data (k4); bkm = Im_Data (k1) - Im_Data (k4); ajp = Re_Data (k2) + Re_Data (k3); ajm = Re_Data (k2) - Re_Data (k3); bjp = Im_Data (k2) + Im_Data (k3); bjm = Im_Data (k2) - Im_Data (k3); aa = Re_Data (kk); bb = Im_Data (kk); Re_Data (kk) = aa + akp + ajp; Im_Data (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; Re_Data (k1) = ak - bj; Re_Data (k4) = ak + bj; Im_Data (k1) = bk + aj; Im_Data (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; Re_Data (k2) = ak - bj; Re_Data (k3) = ak + bj; Im_Data (k2) = bk + aj; Im_Data (k3) = bk - aj; kk = k4 + kspan; } while (kk < nn); kk -= nn; } while (kk <= kspan); break; default: k = factor [ii - 1]; if (jf != k) { jf = k; s1 = pi2 / (double) jf; c1 = cos (s1); s1 = sin (s1); if (jf > maxFactors) goto Memory_Error; Cos [jf - 1] = 1.0; Sin [jf - 1] = 0.0; j = 1; do { Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1; Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1; k--; Cos [k - 1] = Cos [j - 1]; Sin [k - 1] = -Sin [j - 1]; j++; } while (j < k); } do { do { REAL aa, ak; REAL bb, bk; int k2; aa = ak = Re_Data (kk); bb = bk = Im_Data (kk); k1 = kk; k2 = kk + ispan; j = 1; k1 += kspan; do { k2 -= kspan; Rtmp [j] = Re_Data (k1) + Re_Data (k2); ak += Rtmp [j]; Itmp [j] = Im_Data (k1) + Im_Data (k2); bk += Itmp [j]; j++; Rtmp [j] = Re_Data (k1) - Re_Data (k2); Itmp [j] = Im_Data (k1) - Im_Data (k2); j++; k1 += kspan; } while (k1 < k2); Re_Data (kk) = ak; Im_Data (kk) = bk; k1 = kk; k2 = kk + ispan; j = 1; do { REAL aj = 0.0; REAL bj = 0.0; k1 += kspan; k2 -= kspan; jj = j; ak = aa; bk = bb; k = 1; do { ak += Rtmp [k] * Cos [jj - 1]; bk += Itmp [k] * Cos [jj - 1]; k++; aj += Rtmp [k] * Sin [jj - 1]; bj += Itmp [k] * Sin [jj - 1]; k++; jj += j; if (jj > jf) jj -= jf; } while (k < jf); k = jf - j; Re_Data (k1) = ak - bj; Im_Data (k1) = bk + aj; Re_Data (k2) = ak + bj; Im_Data (k2) = bk - aj; j++; } while (j < k); kk += ispan; } while (kk <= nn); kk -= nn; } while (kk <= kspan); break; } /* multiply by rotation factor (except for factors of 2 and 4) */ if (ii == nFactor) goto Permute_Results; /* exit infinite loop */ kk = jc + 1; do { c2 = 1.0 - cd; s1 = sd; do { c1 = c2; s2 = s1; kk += kspan; do { REAL tmp; do { REAL ak; ak = Re_Data (kk); Re_Data (kk) = c2 * ak - s2 * Im_Data (kk); Im_Data (kk) = s2 * ak + c2 * Im_Data (kk); kk += ispan; } while (kk <= nt); tmp = s1 * s2; s2 = s1 * c2 + c1 * s2; c2 = c1 * c2 - tmp; kk = kk - nt + kspan; } while (kk <= ispan); c2 = c1 - (cd * c1 + sd * s1); s1 += sd * c1 - cd * s1; c1 = 2.0 - (c2 * c2 + s1 * s1); s1 *= c1; c2 *= c1; kk = kk - ispan + jc; } while (kk <= kspan); kk = kk - kspan + jc + inc; } while (kk <= jc + jc); break; #endif /* FFT_RADIX4 */ } } /* permute the results to normal order -- done in two stages */ /* permutation for square factors of n */ Permute_Results: Perm [0] = ns; if (kt) { int k2; k = kt + kt + 1; if (k > nFactor) k--; Perm [k] = jc; j = 1; do { Perm [j] = Perm [j - 1] / factor [j - 1]; Perm [k - 1] = Perm [k] * factor [j - 1]; j++; k--; } while (j < k); k3 = Perm [k]; kspan = Perm [1]; kk = jc + 1; k2 = kspan + 1; j = 1; if (nPass != nTotal) { /* permutation for multivariate transform */ Permute_Multi: do { do { k = kk + jc; do { /* swap * Re_Data (kk) <> Re_Data (k2) * Im_Data (kk) <> Im_Data (k2) */ REAL tmp; tmp = Re_Data (kk); Re_Data (kk) = Re_Data (k2); Re_Data (k2) = tmp; tmp = Im_Data (kk); Im_Data (kk) = Im_Data (k2); Im_Data (k2) = tmp; kk += inc; k2 += inc; } while (kk < k); kk += (ns - jc); k2 += (ns - jc); } while (kk < nt); k2 = k2 - nt + kspan; kk = kk - nt + jc; } while (k2 < ns); do { do { k2 -= Perm [j - 1]; j++; k2 = Perm [j] + k2; } while (k2 > Perm [j - 1]); j = 1; do { if (kk < k2) goto Permute_Multi; kk += jc; k2 += kspan; } while (k2 < ns); } while (kk < ns); } else { /* permutation for single-variate transform (optional code) */ Permute_Single: do { /* swap * Re_Data (kk) <> Re_Data (k2) * Im_Data (kk) <> Im_Data (k2) */ REAL t; t = Re_Data (kk); Re_Data (kk) = Re_Data (k2); Re_Data (k2) = t; t = Im_Data (kk); Im_Data (kk) = Im_Data (k2); Im_Data (k2) = t; kk += inc; k2 += kspan; } while (k2 < ns); do { do { k2 -= Perm [j - 1]; j++; k2 = Perm [j] + k2; } while (k2 > Perm [j - 1]); j = 1; do { if (kk < k2) goto Permute_Single; kk += inc; k2 += kspan; } while (k2 < ns); } while (kk < ns); } jc = k3; } if ((kt << 1) + 1 >= nFactor) return 0; ispan = Perm [kt]; /* permutation for square-free factors of n */ j = nFactor - kt; factor [j] = 1; do { factor [j - 1] *= factor [j]; j--; } while (j != kt); nn = factor [kt] - 1; kt++; if (nn > maxPerm) goto Memory_Error; j = jj = 0; for (;;) { int k2; k = kt + 1; k2 = factor [kt - 1]; kk = factor [k - 1]; j++; if (j > nn) break; /* exit infinite loop */ jj += kk; while (jj >= k2) { jj -= k2; k2 = kk; kk = factor [k++]; jj += kk; } Perm [j - 1] = jj; } /* determine the permutation cycles of length greater than 1 */ j = 0; for (;;) { do { kk = Perm [j++]; } while (kk < 0); if (kk != j) { do { k = kk; kk = Perm [k - 1]; Perm [k - 1] = -kk; } while (kk != j); k3 = kk; } else { Perm [j - 1] = -j; if (j == nn) break; /* exit infinite loop */ } } maxFactors *= inc; /* reorder a and b, following the permutation cycles */ for (;;) { j = k3 + 1; nt -= ispan; ii = nt - inc + 1; if (nt < 0) break; /* exit infinite loop */ do { do { j--; } while (Perm [j - 1] < 0); jj = jc; do { int k2; if (jj < maxFactors) kspan = jj; else kspan = maxFactors; jj -= kspan; k = Perm [j - 1]; kk = jc * k + ii + jj; k1 = kk + kspan; k2 = 0; do { Rtmp [k2] = Re_Data (k1); Itmp [k2] = Im_Data (k1); k2++; k1 -= inc; } while (k1 != kk); do { k1 = kk + kspan; k2 = k1 - jc * (k + Perm [k - 1]); k = -Perm [k - 1]; do { Re_Data (k1) = Re_Data (k2); Im_Data (k1) = Im_Data (k2); k1 -= inc; k2 -= inc; } while (k1 != kk); kk = k2; } while (k != j); k1 = kk + kspan; k2 = 0; do { Re_Data (k1) = Rtmp [k2]; Im_Data (k1) = Itmp [k2]; k2++; k1 -= inc; } while (k1 != kk); } while (jj); } while (j != 1); } return 0; /* exit point here */ /* alloc or other problem, do some clean-up */ Memory_Error: fprintf (stderr, "Error: " FFTRADIXS "() - insufficient memory.\n"); fft_free (); /* free-up memory */ return -1; } #endif /* _FFTN_C */ /*----------------------- end-of-file (C source) -----------------------*/