diff options
Diffstat (limited to 'fft/fftn.c')
-rw-r--r-- | fft/fftn.c | 1046 |
1 files changed, 1046 insertions, 0 deletions
diff --git a/fft/fftn.c b/fft/fftn.c new file mode 100644 index 0000000..6ee7c5a --- /dev/null +++ b/fft/fftn.c @@ -0,0 +1,1046 @@ +/*--------------------------------*-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 int +fftradix (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 int +fftradixf (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 11 +static int factor [NFACTOR]; + +void +fft_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 */ + +/* + * + */ +int +FFTN (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 int +FFTRADIX (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; + 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 { + k2 = kk + kspan; + ak = Re [k2]; + bk = Im [k2]; + Re [k2] = Re [kk] - ak; + Im [k2] = Im [kk] - bk; + Re [kk] += ak; + Im [kk] += bk; + kk = k2 + kspan; + } while (kk <= nn); + kk -= nn; + } while (kk <= jc); + if (kk > kspan) + goto Permute_Results_Label; /* exit infinite loop */ + do { + c1 = 1.0 - cd; + s1 = sd; + do { + do { + do { + k2 = kk + kspan; + ak = Re [kk] - Re [k2]; + bk = Im [kk] - Im [k2]; + Re [kk] += Re [k2]; + Im [kk] += Im [k2]; + Re [k2] = c1 * ak - s1 * bk; + Im [k2] = s1 * ak + c1 * bk; + kk = k2 + kspan; + } while (kk < nt); + k2 = kk - nt; + c1 = -c1; + kk = k1 - k2; + } while (kk > k2); + ak = c1 - (cd * c1 + sd * s1); + s1 = sd * c1 - cd * s1 + s1; + c1 = 2.0 - (ak * ak + s1 * s1); + s1 *= c1; + c1 *= ak; + 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 { + k1 = kk + kspan; + k2 = k1 + kspan; + k3 = k2 + kspan; + akp = Re [kk] + Re [k2]; + akm = Re [kk] - Re [k2]; + ajp = Re [k1] + Re [k3]; + ajm = Re [k1] - Re [k3]; + bkp = Im [kk] + Im [k2]; + bkm = Im [kk] - Im [k2]; + bjp = Im [k1] + Im [k3]; + bjm = Im [k1] - Im [k3]; + Re [kk] = akp + ajp; + Im [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 [k1] = akp; + Re [k2] = ajp; + Re [k3] = akm; + Im [k1] = bkp; + Im [k2] = bjp; + Im [k3] = bkm; + } else { + Re [k1] = akp * c1 - bkp * s1; + Re [k2] = ajp * c2 - bjp * s2; + Re [k3] = akm * c3 - bkm * s3; + Im [k1] = akp * s1 + bkp * c1; + Im [k2] = ajp * s2 + bjp * c2; + Im [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_Label; /* exit infinite loop */ + break; + + default: + /* transform for odd factors */ +#ifdef FFT_RADIX4 + fputs ("Error: " FFTRADIXS "(): compiled for radix 2/4 only\n", stderr); + fft_free (); /* free-up memory */ + return -1; + break; +#else /* FFT_RADIX4 */ + k = factor [ii - 1]; + ispan = kspan; + kspan /= k; + + switch (k) { + case 3: /* transform for factor of 3 (optional code) */ + do { + do { + k1 = kk + kspan; + k2 = k1 + kspan; + ak = Re [kk]; + bk = Im [kk]; + aj = Re [k1] + Re [k2]; + bj = Im [k1] + Im [k2]; + Re [kk] = ak + aj; + Im [kk] = bk + bj; + ak -= 0.5 * aj; + bk -= 0.5 * bj; + aj = (Re [k1] - Re [k2]) * s60; + bj = (Im [k1] - Im [k2]) * s60; + Re [k1] = ak - bj; + Re [k2] = ak + bj; + Im [k1] = bk + aj; + Im [k2] = bk - 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 { + k1 = kk + kspan; + k2 = k1 + kspan; + k3 = k2 + kspan; + k4 = k3 + kspan; + akp = Re [k1] + Re [k4]; + akm = Re [k1] - Re [k4]; + bkp = Im [k1] + Im [k4]; + bkm = Im [k1] - Im [k4]; + ajp = Re [k2] + Re [k3]; + ajm = Re [k2] - Re [k3]; + bjp = Im [k2] + Im [k3]; + bjm = Im [k2] - Im [k3]; + aa = Re [kk]; + bb = Im [kk]; + Re [kk] = aa + akp + ajp; + Im [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 [k1] = ak - bj; + Re [k4] = ak + bj; + Im [k1] = bk + aj; + Im [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 [k2] = ak - bj; + Re [k3] = ak + bj; + Im [k2] = bk + aj; + Im [k3] = bk - aj; + kk = k4 + kspan; + } while (kk < nn); + kk -= nn; + } while (kk <= kspan); + break; + + default: + if (k != jf) { + jf = k; + s1 = pi2 / (double) k; + c1 = cos(s1); + s1 = sin(s1); + if (jf > max_factors) + goto Memory_Error_Label; + 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 { + k1 = kk; + k2 = kk + ispan; + ak = aa = Re [kk]; + bk = bb = Im [kk]; + j = 1; + k1 += kspan; + do { + k2 -= kspan; + j++; + Rtmp [j - 1] = Re [k1] + Re [k2]; + ak += Rtmp [j - 1]; + Itmp [j - 1] = Im [k1] + Im [k2]; + bk += Itmp [j - 1]; + j++; + Rtmp [j - 1] = Re [k1] - Re [k2]; + Itmp [j - 1] = Im [k1] - Im [k2]; + k1 += kspan; + } while (k1 < k2); + Re [kk] = ak; + Im [kk] = bk; + k1 = kk; + k2 = kk + ispan; + j = 1; + do { + k1 += kspan; + k2 -= kspan; + jj = j; + ak = aa; + bk = bb; + aj = 0.0; + bj = 0.0; + k = 1; + do { + k++; + ak += Rtmp [k - 1] * Cos [jj - 1]; + bk += Itmp [k - 1] * Cos [jj - 1]; + k++; + aj += Rtmp [k - 1] * Sin [jj - 1]; + bj += Itmp [k - 1] * Sin [jj - 1]; + jj += j; + if (jj > jf) { + jj -= jf; + } + } while (k < jf); + k = jf - j; + Re [k1] = ak - bj; + Im [k1] = bk + aj; + Re [k2] = ak + bj; + Im [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 == mfactor) + goto Permute_Results_Label; /* exit infinite loop */ + kk = jc + 1; + do { + c2 = 1.0 - cd; + s1 = sd; + do { + c1 = c2; + s2 = s1; + kk += kspan; + do { + do { + ak = Re [kk]; + Re [kk] = c2 * ak - s2 * Im [kk]; + Im [kk] = s2 * ak + c2 * Im [kk]; + kk += ispan; + } while (kk <= nt); + ak = s1 * s2; + s2 = s1 * c2 + c1 * s2; + c2 = c1 * c2 - ak; + 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_Label: + Perm [0] = ns; + if (kt) { + k = kt + kt + 1; + if (mfactor < k) + k--; + j = 1; + Perm [k] = jc; + 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_Label: + do { + do { + k = kk + jc; + do { + /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ + ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; + bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; + 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_Label; + kk += jc; + k2 += kspan; + } while (k2 < ns); + } while (kk < ns); + } else { +/* permutation for single-variate transform (optional code) */ +Permute_Single_Label: + do { + /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */ + ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak; + bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk; + 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_Label; + kk += inc; + k2 += kspan; + } while (k2 < ns); + } while (kk < ns); + } + jc = k3; + } + + if ((kt << 1) + 1 >= mfactor) + return 0; + ispan = Perm [kt]; + /* permutation for square-free factors of n */ + j = mfactor - kt; + factor [j] = 1; + do { + factor [j - 1] *= factor [j]; + j--; + } while (j != kt); + kt++; + nn = factor [kt - 1] - 1; + if (nn > max_perm) + goto Memory_Error_Label; + j = jj = 0; + for (;;) { + 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; + k++; + kk = factor [k - 1]; + jj += kk; + } + Perm [j - 1] = jj; + } + /* determine the permutation cycles of length greater than 1 */ + j = 0; + for (;;) { + do { + j++; + kk = Perm [j - 1]; + } 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 */ + } + } + max_factors *= 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 { + kspan = jj; + if (jj > max_factors) { + kspan = max_factors; + } + jj -= kspan; + k = Perm [j - 1]; + kk = jc * k + ii + jj; + k1 = kk + kspan; + k2 = 0; + do { + k2++; + Rtmp [k2 - 1] = Re [k1]; + Itmp [k2 - 1] = Im [k1]; + k1 -= inc; + } while (k1 != kk); + do { + k1 = kk + kspan; + k2 = k1 - jc * (k + Perm [k - 1]); + k = -Perm [k - 1]; + do { + Re [k1] = Re [k2]; + Im [k1] = Im [k2]; + k1 -= inc; + k2 -= inc; + } while (k1 != kk); + kk = k2; + } while (k != j); + k1 = kk + kspan; + k2 = 0; + do { + k2++; + Re [k1] = Rtmp [k2 - 1]; + Im [k1] = Itmp [k2 - 1]; + 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_Label: + fputs ("Error: " FFTRADIXS "() - insufficient memory.\n", stderr); + fft_free (); /* free-up memory */ + return -1; +} +#endif /* _FFTN_C */ +/* ---------------------- end-of-file (c source) ---------------------- */ |