/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "scft.h" #include #include /* PARAMETER translations */ #define KEDIM 30 #define MAXM KEDIM #define NDMAX 6 #define ONE 1.e0 #define SPI4 .7071067811865475244008443621048490e0 /* end of PARAMETER translations */ /* COMMON translations */ struct t_csfftc { LOGICAL32 needst; long int mt, nt, mm, ks, ilast, ke[KEDIM]; } csfftc; /* end of COMMON translations */ void /*FUNCTION*/ scft( float a[], char *mode, long m[], long nd, long *ms, float s[]) { long int _d_l, _d_m, _do0, _do1, i, ii, iic, ir, irc, j, jl, k, ksm, l, ma, mmax, msi, ndd, nf[NDMAX + 1], ntot2; float fn; static char cerror[14] = "MODE(K:K) = "; long int *const kee = (long*)((long*)csfftc.ke + -1); /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const A = &a[0] - 1; long *const Kee = &kee[0] - 1; long *const M = &m[0] - 1; long *const Nf = &nf[0] - 1; float *const S = &s[0] - 1; /* end of OFFSET VECTORS */ /*>> 1997-03-31 SCFT Krogh Increased KEDIM, more sine table checks. *>> 1994-10-20 SCFT Krogh Change CERROR to simplify conversion to C. *>> 1994-10-20 SCFT Krogh Changes to use M77CON *>> 1994-04-19 SCFT CLL Edited to make DP & SP files similar. *>> 1989-06-16 FTK Fix error message on MODE. *>> 1989-06-05 WVS Change length of MODE from *(ND) to *(*). *>> 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 complex Fourier transforms in up to 6 * dimensions 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. * * Values for A, MODE, M, ND, and MS must be specified before * calling the subroutine. * * In describing the usage the following notation is used * N(K) = 2**M(K) * WK = EXP(2*PI*I/N(K)), where I = sqrt(-1) and PI = 3.14159... * MA = M(1) + M(2) + ... + M(ND) * NA = 2**MA * MX = MAX(M(1), M(2), ..., M(ND)) * NX = 2**MX * * 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, C, if one is * doing Fourier synthesis. On output, A contains C 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 complex array with * dimension A(N(1),N(2),...,N(ND)). (Of course this is not necessary * as long as one declares enough storage for A, and keeps in mind that * the array is treated in this subroutine as if it had been specified * in this way.) It is assumed that the imaginary part of a complex * number is stored in the cell immediately following its real part. * * MODE A character variable of length at least ND. * The k-th character must be 'A' or 'a' to select Analysis in the * k-th dimension, or 'S' or 's' to select Synthesis in the * k-th dimension. * The subroutine sets A(j1+1, j2+1, ..., jND+1) to the sum * over 0 .le. k1 .le. N(1)-1, 0 .le. k2 .le. N(2)-1, ..., * 0 .le. kND .le. N(ND)-1 of A(k1+1, k2+1, ..., kND+1) * T(1,j1,k1) * * T(2,j2,k2) * ... * T(ND,jND,kND), 0 .le. j1 .le. N(1)-1, ..., * 0 .le. jnd .le. N(ND)-1 * with T(L,j,k) = * WL ** (j*k) if MODE(L:L) = 'S' or 's' * (1/N(L)) * WL ** (-j*k) if MODE(L:L) = 'A' or 'a' * * M() is a vector used to indicate N(k) = 2**M(k), the number of * complex points in the k-th variable. M must be such that MX < 21. * If M(L)=0, no change is made with respect to the L-th dimension. * * ND is 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 ND, and MAXM = the maximum * permitted for M(1), ..., M(ND) * The dimension of KE must be at least as large as MAXM. * The named common CSFFTC is used for communication between this * subroutine and the subroutine SFFT which computes a one * dimensional complex Fourier transform and computes the sine table. * The use of the variables in CSFFTC is contained in the listing * of SFFT. * * NF(1) = 2, NF(K+1) = 2**(1+M(1)+...+M(K)), K = 1,..., ND * ------------------------------------------------------------------ *--S replaces "?": ?CFT, ?FFT, C?FFTC * Both versions use IERM1 * ------------------------------------------------------------------ */ /* Common variables */ /* Note that KEE(1) is equivalent to ILAST. */ /* ------------------------------------------------------------------ * */ ndd = nd; if ((ndd <= 0) || (ndd > NDMAX)) { /* Fatal error, default is to stop in IERM1 */ ierm1( "SCFT", 1, 2, "BAD ND", "ND", nd, '.' ); *ms = -2; return; } Nf[1] = 2; ma = 0; mmax = 0; fn = ONE; for (k = 1; k <= ndd; k++) { csfftc.mm = M[k]; mmax = max( csfftc.mm, mmax ); ma += csfftc.mm; i = ipow(2,csfftc.mm); if ((csfftc.mm < 0) || (csfftc.mm > MAXM)) { /* Fatal error, default is to stop in IERM1 */ ierm1( "SCFT", 2, 2, "Need 0 .le. M(K) .le. 30", "M(K)" , csfftc.mm, '.' ); *ms = -2; return; } Nf[k + 1] = Nf[k]*i; if (mode[k - 1] == 'A' || mode[k - 1] == 'a') { fn /= i; } else if (mode[k - 1] != 'S' && mode[k - 1] != 's') { cerror[12] = mode[k - 1]; ierm1( "SCFT", 3, 2, cerror, "for K =", k, '.' ); *ms = -2; return; } } msi = *ms; csfftc.needst = mmax > msi; if (!csfftc.needst) { /* Check internal parameters to catch certain user errors. */ if (csfftc.mt < KEDIM) { if (mmax <= csfftc.mt + 2) { /* Skip sine table computation if all appears O.K. */ if (csfftc.mt <= 0) goto L_20; if (fabsf( S[csfftc.nt/2] - SPI4 ) <= 1.e-7) goto L_20; } } csfftc.needst = TRUE; ermsg( "DRFT1", 3, 1, "Invalid sine table (re)computed", '.' ); } *ms = mmax; csfftc.mt = mmax - 2; sfft( a, a, s ); if (msi == -1) return; /* All setup for computation now */ L_20: ntot2 = Nf[ndd + 1]; if (fn != ONE) { for (i = 1; i <= ntot2; i++) { A[i] *= fn; } } /* Beginning of loop for computing multiple sum */ for (k = 1; k <= ndd; k++) { if (mode[k - 1] == 'A' || mode[k - 1] == 'a') { irc = 1; iic = 0; } else { irc = 0; iic = 1; } csfftc.mm = M[k]; csfftc.ks = Nf[k]; ksm = csfftc.ks - 1; csfftc.ilast = Nf[k + 1]; for (l = 1; l <= csfftc.mm; l++) { Kee[l + 1] = Kee[l]/2; } for (j = 1, _do0=DOCNT(j,ntot2,_do1 = csfftc.ilast); _do0 > 0; j += _do1, _do0--) { jl = j + ksm; for (i = j; i <= jl; i += 2) { ir = i + irc; ii = i + iic; sfft( &A[ir], &A[ii], s ); } } } return; /* End of loop for computing multiple sum * */ } /* end of function */