/*--------------------------------*-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 <olesen@me.QueensU.CA>
* 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[], size_t nTotal, size_t nPass,
* size_t nSpan, int iSign, size_t maxFactors,
* size_t 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 <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.
*
* 8 Jan 95 mj olesen <olesen@me.QueensU.CA>
* - 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 <olesen@me.QueensU.CA>
* - 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 <Michael.Steffens@mbox.muk.uni-hannover.de>
* Fortran 90 module
* - made it simpler to pass dimensions for 1D FFT.
*
* 23 Feb 97 Mark Olesen <olesen@me.QueensU.CA>
* - 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.
*
* 11 Nov 97 Robin Williams <rjrw@ast.leeds.ac.uk>
* - fixed all temporaries to double precision using REALFIX,
* following suggestion from Scott Wurcer
* ----------------------------------------------------------------------*/
/*}}}*/
#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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#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 (int nPass, int * kt)
{
int nFactor = 0;
int j, 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 REALFIX
# undef REAL
# undef FFTN
# undef FFTNS
# undef FFTRADIX
# undef FFTRADIXS
/* defines for double */
# define REAL double
# define REALFIX double
# define FFTN fftn
# define FFTNS "fftn"
# define FFTRADIX fftradix
# define FFTRADIXS "fftradix"
/* double precision routine */
static int
fftradix (double Re[], double Im[],
size_t nTotal, size_t nPass, size_t nSpan, int isign,
int maxFactors, int maxPerm);
# include __FILE__ /* include this file again */
#endif
/*}}}*/
/*{{{ defines for re-including float precision */
#ifdef FFT_FLOAT
# undef REAL
# undef REALFIX
# undef FFTN
# undef FFTNS
# undef FFTRADIX
# undef FFTRADIXS
/* defines for float */
# define REAL float
# define REALFIX double
# 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[],
size_t nTotal, size_t nPass, size_t nSpan, int isign,
int maxFactors, 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 int dims [],
REAL Re [],
REAL Im [],
int iSign,
REAL scaling)
{
size_t nTotal;
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)
{
size_t 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)
{
int i;
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 (i = 0; i < nTotal; i += iSign)
{
Re_Data (i) *= scaling;
Im_Data (i) *= 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 [],
size_t nTotal,
size_t nPass,
size_t nSpan,
int iSign,
int maxFactors,
int maxPerm)
{
int ii, nFactor, kspan, ispan, inc;
int j, jc, jf, jj, k, k1, k3, kk, kt, nn, ns, nt;
REALFIX radf;
REALFIX c1, c2, c3, cd;
REALFIX s1, s2, s3, sd;
REALFIX * Rtmp = NULL; /* temp space for real part*/
REALFIX * Itmp = NULL; /* temp space for imaginary part */
REALFIX * Cos = NULL; /* Cosine values */
REALFIX * Sin = NULL; /* Sine values */
#ifndef FFT_RADIX4
REALFIX s60 = SIN60; /* sin(60 deg) */
REALFIX s72 = SIN72; /* sin(72 deg) */
REALFIX c72 = COS72; /* cos(72 deg) */
#endif
REALFIX 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 (REALFIX))
{
#ifdef SUN_BROKEN_REALLOC
if (!SpaceAlloced) /* first time */
{
SpaceAlloced = maxFactors * sizeof (REALFIX);
Tmp0 = malloc (SpaceAlloced);
Tmp1 = malloc (SpaceAlloced);
Tmp2 = malloc (SpaceAlloced);
Tmp3 = malloc (SpaceAlloced);
}
else
{
#endif
SpaceAlloced = maxFactors * sizeof (REALFIX);
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 (REALFIX);
}
if (MaxPermAlloced < 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 = (REALFIX *) Tmp0;
Itmp = (REALFIX *) Tmp1;
Cos = (REALFIX *) Tmp2;
Sin = (REALFIX *) 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 {
REALFIX tmpr;
REALFIX 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 {
REALFIX tmp;
do {
do {
REALFIX tmpr;
REALFIX 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 {
REALFIX ajm, ajp, akm, akp;
REALFIX 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 {
REALFIX aj, tmpr;
REALFIX 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 {
REALFIX aa, aj, ak, ajm, ajp, akm, akp;
REALFIX 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 {
REALFIX aa, ak;
REALFIX 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 {
REALFIX aj = 0.0;
REALFIX 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 {
REALFIX tmp;
do {
REALFIX 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)
*/
REALFIX 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)
*/
REALFIX 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) -----------------------*/