/*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 "drft.h" #include #include /* 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 */