/*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 */