/*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
#include "fcrt.h"
#include "drft1.h"
#include
#include
/* PARAMETER translations */
#define HALF .5e0
#define KEDIM 30
#define MMAX (KEDIM + 1)
#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*/ drft1(
double a[],
byte mode,
long m,
long *ms,
double s[])
{
LOGICAL32 anal;
long int i, ii, ii1, ii2, ir, ir1, ir2, j, jdif, jj, k1, k1n,
kn2, l, ma, msi, n1p;
double fn, t, ti, tt, tti, wi, wr;
static char msg1[20] = "Bad MODE. MODE = ";
long int *const kee = (long*)((long*)cdfftc.ke + -1);
long int *const n1 = (long*)&cdfftc.ilast;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const A = &a[0] - 1;
long *const Kee = &kee[0] - 1;
double *const S = &s[0] - 1;
/* end of OFFSET VECTORS */
/*>> 1997-03-31 DFFT1 Krogh Increased KEDIM, more sine table checks.
*>> 1996-01-23 DRFT1 Krogh Changes to simplify conversion to C.
*>> 1994-11-11 DRFT1 Krogh Declared all vars.
*>> 1994-10-20 DRFT1 Krogh Changes to use M77CON
*>> 1989-05-07 DRFT1 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 using
* the Cooley-Tukey fast Fourier transform.
*
* Variables in the calling sequence have the following types */
/* Programmed by Fred T. Krogh at the Jet Propulsion Laboratory,
* Pasadena, Calif. August 1, 1969.
* Revised for portability by Krogh -- January 22, 1988
*
* In describing the usage the following notation is used
* N = 2 ** M
* W = EXP(2*PI*I/N), where I = SQRT(-1) and PI = 3.14159...
*
* The usage is as follows
*
* A() is an array of function values if one is doing Fourier analysis,
* and is an array of Fourier coefficients if one is doing Fourier
* synthesis. In our description here we assume that A is a real
* array with dimension A(N) when A contains the real data, X, and
* that A is a complex array with dimension A(N/2) when A contains
* complex Fourier coefficients, C. (C(k) for k > N/2 need not be
* saved since for 0 < k < N, C(N-k) = conjugate of C(k). It is
* assumed that the imaginary part of a complex number is stored
* in the cell immediately following its real part, except that
* A(1) = C(0), and A(2) = C(N/2). This is possible since these
* values of C are real and by doing this both X and C require the
* same storage in A. Of course the storage required for A can be
* reserved by the user in any way that works.
*
* MODE Selects Synthesis or Analysis.
* If MODE = 'A' or 'a', do Fourier analysis, which amounts to setting
* C(k) = sum for j=0 to N-1 of X(j)*T(M,j,k), for k = 0, N/2
* with T(M,j,k) = (1/N) * W ** (-j*k).
* If MODE = 'S' or 's', do Fourier synthesis, which amounts to setting
* X(j) = sum for k=0 to N-1 of C(k)*T(M,j,k), for j = 0, N - 1
* with T(M,j,k) = W ** (j*k)
* (Recall that C(N-k) is the conjugate of C(k).)
*
* M is used to indicate N = 2**M, the number of real points. The
* number of points must satisfy 1 .le. N .le. 2**21.
* M = 0 gives an immediate return.
*
* 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 M, which is the value of MS required for computing a
* transform of size N. 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.
* On detected errors the error message subrs are called and
* execution stops. If the user overrides the stop to cause
* continuation, then this subr will return with MS = -2.
*
* 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 M .gt. MS. (If S is altered, set MS=0 so that S
* is recomputed.)
*
* ------------------------------------------------------------------
* Notes on COMMON, PARAMETER's, and local variables
*
* MMAX is the largest value allowed for M
* The dimension of KE must be at least as large as MMAX-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.
*
* ANAL = .TRUE. when doing Fourier analysis, and .false. otherwise.
*
* N1 = 2 ** M
* ------------------------------------------------------------------
*--D replaces "?": ?RFT1, ?FFT, C?FFTC
* Both versions use ERMSG, IERM1
* and need ERFIN, IERV1
* ------------------------------------------------------------------ */
/* 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( "DRFT1", 2, 2, msg1, '.' );
*ms = -2;
return;
}
ma = m;
if (ma > MMAX || ma < 0)
{
/* Fatal error, default is to stop in IERM1 */
ierm1( "DRFT1", 1, 2, "Require 0 .le. M .le. 31", "M", m,
'.' );
*ms = -2;
return;
}
msi = *ms;
cdfftc.needst = ma > msi;
if (!cdfftc.needst)
{
/* Check internal parameters to catch certain user errors. */
if (cdfftc.mt < KEDIM)
{
if (ma <= cdfftc.mt + 2)
{
/* Skip sine table computation if all appears O.K. */
if (cdfftc.mt <= 0)
goto L_10;
if (fabs( S[cdfftc.nt/2] - SPI4 ) <= 1.e-7)
goto L_10;
}
}
cdfftc.needst = TRUE;
ermsg( "DRFT1", 3, 1, "Invalid sine table (re)computed", '.' );
}
*ms = ma;
cdfftc.mt = ma - 2;
dfft( a, a, s );
/* Return if user requested it. */
if (msi == -1)
return;
/* All setup for computation now */
L_10:
if (ma != 0)
{
cdfftc.mm = ma - 1;
*n1 = ipow(2,ma);
n1p = *n1 + 2;
kn2 = *n1/2;
jdif = (4*cdfftc.nt)/ *n1;
cdfftc.ks = 2;
if (anal)
{
/* Set flags for Fourier analysis */
ir = 2;
ii = 1;
fn = HALF/(double)( *n1 );
/* Doing Fourier analysis, so multiply by 2 ** M */
for (i = 1; i <= *n1; i++)
{
A[i] *= fn;
}
}
else
{
/* Set flags for Fourier synthesis */
ir = 1;
ii = 2;
goto L_50;
}
/* Compute complex Fourier transform */
L_30:
for (l = 1; l <= cdfftc.mm; l++)
{
Kee[l + 1] = Kee[l]/2;
}
dfft( &A[ir], &A[ii], s );
/* End of computing complex transform
* */
if (!anal)
goto L_70;
/* Beginning of calculations relating coefficients of real data
* with coefficients of associated complex data
*
* Special case -- K1 = 0 */
L_50:
t = A[1] + A[2];
ti = A[1] - A[2];
if (anal)
{
t *= TWO;
ti *= TWO;
}
A[1] = t;
A[2] = ti;
if (cdfftc.mm > 0)
{
/* Special kase -- K1 = N1 / 4 */
A[kn2 + 1] *= TWO;
A[kn2 + 2] *= -TWO;
if (cdfftc.mm > 1)
{
j = 0;
for (k1 = 3; k1 <= kn2; k1 += 2)
{
k1n = n1p - k1;
if (anal)
{
ir1 = k1n;
ir2 = k1;
}
else
{
ir1 = k1;
ir2 = 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 (!anal)
goto L_30;
}
L_70:
return;
} /* end of function */