/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:43 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "drft.h"
#include <stdlib.h>
#include <string.h>
/* PARAMETER translations */
#define HALF .5e0
#define KEDIM 30
#define MAXM (KEDIM + 1)
#define NDMAX 6
#define SPI4 .7071067811865475244008443621048490e0
#define TWO 2.e0
/* end of PARAMETER translations */
/* COMMON translations */
struct t_cdfftc {
LOGICAL32 needst;
long int mt, nt, mm, ks, ilast, ke[KEDIM];
} cdfftc;
/* end of COMMON translations */
void /*FUNCTION*/ drft(
double a[],
byte mode,
long m[],
long nd,
long *ms,
double s[])
{
LOGICAL32 anal;
long int _d_l, _d_m, _do0, _do1, i, i1, i2, ii, ii1, ii2, iib,
iic, ir, ir1, ir2, irb, irc, j, jdif, jj, jl, k, k1, k1l, k1n,
kb, kk, kn2, kn2s, ksm, l, lb, m1, ma, mmax, msi, mu[NDMAX - 1],
nda, ntot2, _i, _r;
double fn, t, ti, tt, tti, wi, wr;
static char msg1[20] = "Bad MODE. MODE = ";
static int _aini = 1;
/* EQUIVALENCE translations */
static long _es0[8];
long int *const kee = (long*)((long*)cdfftc.ke + -1);
long int *const n1 = (long*)((long*)_es0 + 1);
long int *const nf = (long*)_es0;
/* end of EQUIVALENCE translations */
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const A = &a[0] - 1;
long *const Kee = &kee[0] - 1;
long *const M = &m[0] - 1;
long *const Mu = &mu[0] - 1;
long *const Nf = &nf[0] - 1;
double *const S = &s[0] - 1;
/* end of OFFSET VECTORS */
if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
Nf[1] = 1;
_aini = 0;
}
/*>> 1997-03-31 DFFT Krogh Increased KEDIM, more sine table checks.
*>> 1996-01-23 DRFT Krogh Changes to simplify conversion to C.
*>> 1994-10-20 DRFT Krogh Changes to use M77CON
*>> 1990-01-23 DRFT CLL Removed extraneous label 80.
*>> 1989-06-16 FTK Fix error message on MODE.
*>> 1989-05-10 FTK & CLL
*>> 1989-05-07 FTK & CLL
*>> 1989-04-21 FTK & CLL
* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*
* This subroutine computes Fourier transforms of real data in up to 6
* dimensions using the Cooley-Tukey fast Fourier transform.
*
* Programmed by Fred F. Krogh at the Jet Propulsion Laboratory,
* Pasadena, Calif. August 1, 1969.
* Revised for portability by Krogh -- February 2, 1988
*
* Variables in the calling sequence have the following types */
/* Values for A(), MODE, M(), ND, and MS must be specified before
* calling the subroutine.
*
* In describing the usage the following notation is used
* NDA = ND
* N(K) = 2 ** M(K)
* MA = M(1) + M(2) + ... + M(ND)
* NA = 2 ** MA
* MX = MAX(M(1),M(2), ..., M(ND))
* NX = 2**MX
* WK = EXP(2*PI*I/N(K)), WHERE I = SQRT(-1) AND PI = 3.14159...
*
* The usage is as follows
*
* A() on input is an array of function values, X, if one is doing
* Fourier analysis, and is an array of Fourier coefficients, B, if one
* is doing Fourier synthesis. On output, A contains B if doing
* analysis, and contains X if doing synthesis. In our description
* here, we assume that A has been declared by the user as a real array
* with dimension A(N(1), N(2), ..., N(NDA)) when A contains the real
* data, and A is a complex array, C, with dimension C(N(1)/2,N(2),...,
* N(NDA)), when A contains complex Fourier coefficients. (B(k1, k2,
* ..., kNDA) for k1 > N(1)/2 need not be saved since B(N(1)-k1,
* N(2)-k2, ..., N(NDA)-kNDA) = conjugate of B(k1, k2, ..., kNDA),
* where the L-th subscripts are to be interpreted modulo N(L). It is
* assumed that the imaginary part of a complex number is stored
* in the cell immediately following its real part, except that
* if k1 = 0 and all of the remaining ki satisfy ki = 0 or ki = N(i)/2,
* then B(0,k2,...,kNDA) and B(N(1)/2,k2,...,kNDA) are real, and
* real part of C(1,k2+1,...,kNDA+1) = B(0,k2,...,kNDA) and
* imag. part of C(0,k2+1,...,kNDA+1) = B(N(1)/2,k2,...,kNDA).
* If k1 = 0 and kL is the first ki that does not satisfy ki = 0 or
* ki = N(i)/2 then C(1,k2+1,...,kNDA+1) = B(0,k2,...,kNDA) for
* kL .LT. N(L)/2, and C(1,k2+1,...,kNDA+1) = B(N(1)/2,k2,...,kNDA) for
* kL.GT.N(L)/2. Of course the storage required for A can be
* reserved by the user in any way that works.
*
* MODE A character variable to select Analysis or Synthesis.
* 'A' or 'a' selects Analysis, and 'S' or 's' selects Synthesis.
* For Analysis the subroutine computes
* B(k1,k2,...,kNDA) = sum over 0 .le. j1 .le. N(1)-1, 0 .le. j2
* .le. N(2)-1,..., 0 .le. jNDA .le. N(NDA)-1, of X(j1,j2,...,jNDA)
* * T(1,j1,k1) * T(1,j2,k2) * ... * T(NDA,jNDA,kNDA), 0 .le.
* k1 .le. N(1)-1,..., 0 .le. kNDA .le. N(NDA)-1
* with T(L,j,k) = (1/N(L)) * WL ** (-j*k).
* For Synthesis the subroutine computes
* X(j1,j2,...,jNDA) = sum over 0 .le. k1 .le. N(1)-1, 0 .le. k2
* .le. N(2)-1,..., 0 .le. kNDA .le. N(NDA)-1, of B(k1,k2,...,kNDA)
* * T(1,j1,k1) * T(1,j2,k2) * ... * T(NDA,jNDA,kNDA), 0 .le.
* j1 .le. N(1)-1,..., 0 .le. jNDA .le. N(NDA)-1
* with T(L,j,k) = WL**(j*k).
*
* M() is a vector used to indicate N(k) = 2**M(k), the number of
* real points in the k-th variable. M must be such that 0 .le. M(k)
* .le. 21, and if M(1)=0, then M(k) must be 0 for k =1, ..., NDA.
*
* ND gives the dimension of (i.e. the number of subscripts in)
* the array A. ND must satisfy 1 .le. ND .le. 6.
*
* MS gives the state of the sine table. If MS > 0, there are NT =
* 2 ** (MS-2) good entries in the sine table. On the initial call,
* MS must be set to 0, and when the return is made, it will be set
* to MX, which is the value of MS required for computing a
* transform of size NX. If MS = -1, the sine table will be computed
* as for the case MS = 0, and then a return to the user will be made
* with MS set as before, but no transform will be computed. This
* option is useful if the user would like access to the sine table
* before computing the FFT.
*
* S is a vector, S(j) = sin(pi*j/2*NT)), j = 1, 2, ..., NT-1, where
* NT is defined in the description of MS above. S is computed by the
* subroutine if MX .gt. MS. (If S is altered, set MS=0 so that S
* is recomputed.)
* ------------------------------------------------------------------
* Notes on COMMON, PARAMETER's, and local variables
*
* NDMAX = the maximum value for ABS(ND), and MAXM = the maximum
* permitted for M(1), ..., M(ND)
*
* The dimension of KE must be at least as large as MAXM-1.
*
* The named common CDFFTC is used for communication between this
* subroutine and the subroutine DFFT which computes a one
* dimensional complex Fourier transform and computes the sine table.
* The use of the variables in CDFFTC is contained in the listing
* of DFFT.
*
* NF(1) = 1, NF(K+1) = 2**(M(1) + M(2) + ... + M(K)), K = 1,...,NDA
*
* MU is used in computing the index associated with a given index.
* (The indices (k1,k2,...,kNDA) and (j1,j2,...,jNDA) are associated
* if k1 = (N(1)/2)-j1 (modulo (N(1)/2)) and ki = N(i)-j1 (modulo ni)
* i = 2,...,ND.)
*
* ANAL = .TRUE. when doing Fourier analysis, and .FALSE. otherwise
*
* N1 = 2**M(1)
* ------------------------------------------------------------------
*--D replaces "?": ?RFT, ?FFT, C?FFTC */
/* Both versions use IERM1
* ------------------------------------------------------------------ */
/* Common variables */
/* Note that KEE(1) is equivalent to ILAST. */
/* ------------------------------------------------------------------
* */
if (mode == 'A' || mode == 'a')
{
anal = TRUE;
}
else if (mode == 'S' || mode == 's')
{
anal = FALSE;
}
else
{
msg1[18] = mode;
ermsg( "DRFT", 2, 2, msg1, '.' );
*ms = -2;
return;
}
nda = nd;
if ((nda == 0) || (nda > NDMAX))
{
/* Fatal error, default is to stop in IERM1 */
ierm1( "DRFT", 1, 2, "BAD ND", "ND", nd, '.' );
*ms = -2;
return;
}
ma = 0;
mmax = -1;
for (k = 1; k <= nda; k++)
{
cdfftc.mm = M[k];
if (cdfftc.mm < 0)
goto L_400;
if (cdfftc.mm > mmax)
{
/* if M(1)=0, then so must all of the rest. */
if ((mmax == 0) || (cdfftc.mm > MAXM))
goto L_400;
mmax = cdfftc.mm;
}
ma += cdfftc.mm;
Nf[k + 1] = Nf[k]*ipow(2,cdfftc.mm);
}
msi = *ms;
cdfftc.needst = mmax > msi;
if (!cdfftc.needst)
{
/* Check internal parameters to catch certain user errors. */
if (cdfftc.mt < KEDIM)
{
if (mmax <= cdfftc.mt + 2)
{
/* Skip sine table computation if all appears O.K. */
if (cdfftc.mt <= 0)
goto L_20;
if (fabs( S[cdfftc.nt/2] - SPI4 ) <= 1.e-7)
goto L_20;
}
}
cdfftc.needst = TRUE;
ermsg( "DRFT1", 3, 1, "Invalid sine table (re)computed", '.' );
}
*ms = mmax;
cdfftc.mt = mmax - 2;
dfft( a, a, s );
if (msi == -1)
return;
/* All setup for computation now */
L_20:
if (mmax == 0)
return;
ntot2 = Nf[nda + 1];
m1 = M[1];
/* Compute constants associated with N1. */
kn2 = *n1/2;
kn2s = kn2;
if (anal)
kn2s = -kn2s;
k1l = kn2 - 1;
jdif = (4*cdfftc.nt)/ *n1;
if (!anal)
{
/* Set flags for Fourier synthesis */
irc = 0;
iic = 1;
goto L_160;
}
/* Set flags for Fourier analysis */
irc = 1;
iic = 0;
fn = HALF/(double)( ntot2 );
/* Doing Fourier analysis, so multiply by 2**(-M(1)-M(2)-...-M(ND)-1) */
for (i = 1; i <= ntot2; i++)
{
A[i] *= fn;
}
/* Beginning of loop for computing multiple sum
* multi-dimensional complex Fourier synthesis or analysis */
L_100:
for (k = 1; k <= nda; k++)
{
cdfftc.mm = M[k];
cdfftc.ks = Nf[k];
if (k != 1)
goto L_110;
cdfftc.mm -= 1;
cdfftc.ks = 2;
L_110:
ksm = cdfftc.ks - 1;
cdfftc.ilast = Nf[k + 1];
for (l = 1; l <= cdfftc.mm; l++)
{
Kee[l + 1] = Kee[l]/2;
}
for (j = 1, _do0=DOCNT(j,ntot2,_do1 = cdfftc.ilast); _do0 > 0; j += _do1, _do0--)
{
jl = j + ksm;
for (i = j; i <= jl; i += 2)
{
ir = i + irc;
ii = i + iic;
dfft( &A[ir], &A[ii], s );
}
}
}
/* End of loop for doing multi-dimensional complex Fourier transforms
* */
if (!anal)
return;
/* Beginning of calculations relating coefficients of real data
* with coefficients of associated complex data */
L_160:
kb = 0;
l = 1;
lb = 1;
/* L = LB occurs when the coefficients of the real data which are
* computed in the case k1 = 0, are real.
*
* Coefficients with index L + k1 and LB + N1 - k1 (k1.gt.0) are
* associated in the sense described above where MU is dimensioned.
* When K.GT.1, LB = KB - L.
* */
k = 1;
/* Special case -- K1 = N1/4 (L = LB) */
L_170:
if (m1 != 1)
{
i1 = l + kn2;
A[i1] *= TWO;
A[i1 + 1] *= -TWO;
}
/* Special case -- K1 = 0 (L = LB) */
t = A[l] + A[l + 1];
ti = A[l] - A[l + 1];
if (anal)
{
t *= TWO;
ti *= TWO;
}
A[l] = t;
A[l + 1] = ti;
goto L_220;
/* Special case -- K1 = N1/4 (L.NE.LB) */
L_200:
if (m1 != 1)
{
i1 = l + kn2;
i2 = lb + kn2;
t = TWO*A[i2];
A[i2] = TWO*A[i1];
A[i1] = t;
t = -TWO*A[i2 + 1];
A[i2 + 1] = -TWO*A[i1 + 1];
A[i1 + 1] = t;
}
/* SET UP BASE INDICES */
ir = l;
irb = lb;
if (anal)
{
ir = irb;
irb = l;
}
ii = ir + 1;
iib = irb + 1;
/* SPECIAL CASE-- K1 = 0 (L.NE.LB) */
t = A[ir] + A[irb];
tt = A[iib] + A[ii];
ti = A[ii] - A[iib];
tti = A[ir] - A[irb];
A[ir] = t - tt;
A[irb] = t + tt;
A[ii] = ti + tti;
A[iib] = tti - ti;
L_220:
if (m1 <= 2)
goto L_260;
j = 0;
/* USUAL CASE -- K1.NE.0 AND K1.NE.N1/4 */
for (k1 = 2; k1 <= k1l; k1 += 2)
{
k1n = *n1 - k1;
if (anal)
{
ir1 = lb + k1n;
ir2 = l + k1;
}
else
{
ir1 = l + k1;
ir2 = lb + k1n;
}
ii2 = ir2 + 1;
ii1 = ir1 + 1;
j += jdif;
wi = S[j];
jj = cdfftc.nt - j;
wr = S[jj];
t = A[ir1] - A[ir2];
ti = A[ii1] + A[ii2];
tt = t*wi + ti*wr;
tti = t*wr - ti*wi;
t = A[ir1] + A[ir2];
ti = A[ii1] - A[ii2];
A[ir1] = t - tt;
A[ir2] = t + tt;
A[ii1] = tti + ti;
A[ii2] = tti - ti;
if (l != lb)
{
ir1 += kn2s;
ir2 -= kn2s;
ii2 = ir2 + 1;
ii1 = ir1 + 1;
t = A[ir1] - A[ir2];
ti = A[ii1] + A[ii2];
tt = t*wr - ti*wi;
tti = t*wi + ti*wr;
t = A[ir1] + A[ir2];
ti = A[ii1] - A[ii2];
A[ir1] = t - tt;
A[ir2] = t + tt;
A[ii1] = ti - tti;
A[ii2] = -ti - tti;
}
}
L_260:
if (k < 2)
goto L_310;
L_270:
l += *n1;
kk = 1;
if (k > 2)
{
/* Logic used to find index associated with a given index. */
L_280:
if (Mu[kk] == 0)
kb += Nf[kk + 2];
Mu[kk] += Nf[kk + 1];
if (Mu[kk] >= Nf[kk + 2])
{
Mu[kk] = 0;
kk += 1;
kb -= Nf[kk + 1];
if (kk <= (k - 2))
goto L_280;
}
}
L_300:
lb = kb - l;
if (lb >= l)
{
if (lb > l)
goto L_200;
goto L_170;
}
if (kk <= (k - 2))
goto L_270;
L_310:
if (k >= nda)
goto L_330;
for (i = 1; i <= k; i++)
{
Mu[i] = 0;
}
k += 1;
kb = Nf[k + 1] + 2;
l = Nf[k] + 1;
goto L_300;
L_330:
if (!anal)
goto L_100;
return;
/* Fatal error, default is to stop in IERM1 */
L_400:
ermsg( "DRFT", 2, 2, "M() specified improperly", '.' );
*ms = -2;
return;
} /* end of function */