/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:48 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "stcst.h" #include #include /* PARAMETER translations */ #define FOUR 4.0e0 #define KEDIM 30 #define MAXMX (KEDIM + 1) #define NDMAX 6 #define ONE 1.0e0 #define SPI4 .7071067811865475244008443621048490e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* 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*/ stcst( float a[], char *tcs, char *mode, long m[], long nd, long *ms, float s[]) { char _c0[2]; long int _d_l, _d_m, _do0, _do1, _do2, _do3, i, i1, ii, ir, itcs[NDMAX], itcsk, j, jdif, jj, jk, k, kdr, kii, kin, kk, kki, kkl, kkn, l, ma, mi, mmax, msi, mu[NDMAX], n, ndd, ndiv, nf[NDMAX + 1], ni, ni1, ni2, ni2i, ntot2; float fn, sum, t, t1, tp, tpi, ts, ts1, wi, wr; static char msg1[14] = "MODE(K:K) = ?"; static char msg2[13] = "TCS(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 Itcs = &itcs[0] - 1; long *const Ke = &csfftc.ke[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; float *const S = &s[0] - 1; /* end of OFFSET VECTORS */ /*>> 1997-03-31 STCST Krogh Increased KEDIM, more sine table checks. *>> 1996-01-23 STCST Krogh Changes to simplify conversion to C. *>> 1994-11-11 STCST Krogh Declared all vars. *>> 1994-10-20 STCST Krogh Changes to use M77CON *>> 1989-06-16 STCST FTK Fix error message on MODE, and TCS. *>> 1989-06-05 WVS Change length of MODE and TCS from (ND) to (*) *>> 1989-05-08 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 trigonometirc (sine-cosine), sine, or * cosine transforms of real data in up to 6 dimensions using the * Cooley-Tukey fast Fourier transform. * * Variables in the calling sequence have the following types */ /* Programmed by Fred F. Krogh at the Jet Propulsion Laboratory, * Pasadena, Calif. August 1, 1969. * Revised for portability by Krogh -- January 29, 1988 * * Values for A, TCS, 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) * MA = M(1) + M(2) + ... + M(ND) * NA = 2 ** MA * * MTCS(K) = M(K) TCS(K:K) = 'T' * = M(K)+1 otherwise * * MX = MAX(MTCS(1), MTCS(2), ..., MTCS(ND)) * NX = 2 ** MX * * T(L,j,k) is defined differently for different values of TCS(L) * * if TCS(L:L) = 'T' and MODE(L:L) = 'S', T(L,j,k) * =1/2 if k = 0 * =(1/2)*(-1)**j if k = 1 * =COS(j*k*PI/N(L)) if k IS EVEN (k = 2, 4, ..., N(L)-2) * =SIN(j*(k-1)*PI/N(L)) if k IS ODD (k = 3, 5, ..., N(L)-1) * and if MODE(L:L) = 'A', T(L,j,k) * = (4/N) * (value of T(L,k,j) defined above) If j<2 * = (2/N) * (value of T(L,k,j) defined above) Otherwise * * if TCS(L:L) = 'C' and MODE(L:L) = 'S', T(L,j,k) * =1/2 if k = 0 * =COS(j*k*PI/N(L)) k = 1, 2, ..., N(L)-1 * =(1/2)*(-1)**j if k = N(L) * and if MODE(L:L) = 'A', T(L,j,k) * = (2/N) * (value of T(L,j,k) defined above) * * if TCS(L:L) = 'S' and MODE(L:L) = 'S', T(L,j,k) * =SIN(j*k*PI/N(L)) k = 0, 1, ..., N(L)-1 * and if MODE(L:L) = 'A', T(L,j,k) * = (2/N) * (value of T(L,j,k) defined above) * * D(L) = N(L) if TCS(L:L) .ne. 'C' * = N(L)+1 if TCS(L:L) = 'C' * * The usage is as follows * * A() on input 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. On output, these are reversed. In either case * A() is a real array with dimension A(D(1), D(2), ..., D(ND)). * * TCS is a character variable of length ND. The k-th character must be * 'T' or 't' to select the general Trigonometric transform, or * 'C' or 'c' to select the Cosine transform, or * 'S' or 's' to select the Sine transform. * See the description of T(L,j,k) and M above. * * MODE A character variable of length 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. * One may be doing analysis, MODE(k:k) = 'A', with respect to one * dimension and synthesis, MODE(k:k) = 'S', with respect to * another. A(j1+1, j2+1, ..., jND+1) is replaced by the sum over * 0 .le. k1 .le. D(1)-1, 0 .le. k2 .le. D(2)-1, ..., 0 .le. kND .le. * D(ND)-1, of A(k1+1, k2+1, ..., kND+1) * T(1, k1, j1) * T(2, k2, j2) * ... * T(ND, kND, jND), 0 .le. j1 .le. D(1)-1, ..., 0 .le. jND .le. * D(ND)-1. */ /* M() is a vector used to indicate N(k) = 2**M(k)). The number of * points in the k-th variable is then given by D(k) (see above). M * must be specified in such a way that 0 < M(k) < 22 * for k = 1, 2, ..., ND. * * 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 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 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 MAXMX = the maximum * permitted for MTCS(1), ..., MTCS(ND) * */ /* NF(1) = 1, NF(K+1) = NF(K) * D(K) K = 1, 2, ..., ND * * MU is used in the process of eliminating transforms with respect * to the first subscript of transforms with TCS(:) = 'S'. * (This is only necessary if ND.GT.1.) * * The dimension of KE must be at least as large as MAXMX-1. * 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. * * The input character variable TCS is mapped to the internal * integer array ITCS() by mapping 'T' to 1, 'C' to 2, 'S' to 3. * * ----------------------------------------------------------------- *--S replaces "?": ?TCST, ?FFT, C?FFTC * Both versions use IERM1 * and need ERFIN, IERV1 * ----------------------------------------------------------------- */ /* 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( "STCST", 1, 2, "BAD ND", "ND", nd, '.' ); *ms = -2; return; } ma = 0; mmax = 0; ndiv = 1; /* Every element in the array A is divided by NDIV before computing * the transform. The value computed for NDIV depends on whether * one is doing analysis or synthesis and on the type of * transform being computed. */ for (k = 1; k <= ndd; k++) { csfftc.mm = M[k]; if ((csfftc.mm < 0) || (csfftc.mm > MAXMX)) goto L_200; ma += csfftc.mm; n = ipow(2,csfftc.mm); if (mode[k - 1] == 'A' || mode[k - 1] == 'a') { ndiv *= n; } else if (mode[k - 1] == 'S' || mode[k - 1] == 's') { ndiv *= 2; } else { msg1[12] = mode[k - 1]; ierm1( "STCST", 2, 2, msg1, "for K =", k, '.' ); *ms = -2; return; } itcsk = (istrstr( "TtCcSs", STR1(_c0,tcs[k - 1]) ) + 1)/2; Itcs[k] = itcsk; if (itcsk == 0) { msg2[11] = tcs[k - 1]; ierm1( "STCST", 3, 2, msg2, "for K =", k, '.' ); return; } Nf[1] = 1; if (itcsk >= 2) { if (itcsk == 2) n += 1; ndiv *= 2; csfftc.mm += 1; } Nf[k + 1] = Nf[k]*n; if (csfftc.mm > mmax) { mmax = csfftc.mm; } } 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_15; if (fabsf( S[csfftc.nt/2] - SPI4 ) <= 1.e-7) goto L_15; } } csfftc.needst = TRUE; ermsg( "STCST", 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_15: ntot2 = Nf[ndd + 1]; fn = ONE/(float)( ndiv ); /* Divide every element of A by NDIV */ for (i = 1; i <= ntot2; i++) { A[i] *= fn; } /* Beginning of loop for computing multiple sum */ for (k = 1; k <= ndd; k++) { itcsk = Itcs[k]; mi = M[k]; csfftc.mm = mi - 1; if (mode[k - 1] == 'A' || mode[k - 1] == 'a') mi = -mi; kdr = Nf[k]; csfftc.ks = kdr + kdr; csfftc.ilast = Nf[k + 1]; if (itcsk == 2) csfftc.ilast -= kdr; for (l = 1; l <= csfftc.mm; l++) { Kee[l + 1] = Kee[l]/2; } i = 1; j = ndd; L_40: for (l = 1; l <= j; l++) { Mu[l] = 0; if ((l != k) && (Itcs[l] > 2)) { /* Skip the part of the array left empty by the sine transform */ Mu[l] = Nf[l]; i += Nf[l]; } } /* Compute indices associated with the current value of I (and K) */ L_60: i1 = i + kdr; ni1 = i + Nf[k + 1]; if (itcsk == 2) ni1 -= kdr; ni = ni1 - kdr; ni2 = (ni1 + i)/2; ni2i = ni2 + kdr; if (itcsk != 1) { /* Doing a cosine or a sine transform -- set MI = 0 and do * calculations not required for sine-cosine transforms */ mi = 0; j = ni; sum = A[i1]; t = A[j]; L_70: jk = j - csfftc.ks; if (jk >= i1) { sum += A[j]; A[j] = A[jk] - A[j]; j = jk; goto L_70; } if (itcsk != 2) { /* Calculations for the sine transform */ A[i] = TWO*A[i1]; A[i1] = -TWO*t; if (csfftc.mm == 0) goto L_90; t = TWO*A[ni2]; A[ni2] = -TWO*A[ni2i]; A[ni2i] = t; goto L_90; } /* Set for cosine transform */ A[i1] = A[ni1]; } if (csfftc.mm == 0) goto L_90; if (mi < 0) goto L_160; /* Begin calculations for the sine-cosine transform */ L_80: A[ni2] *= TWO; A[ni2i] *= TWO; L_90: t = A[i] + A[i1]; A[i1] = A[i] - A[i1]; if (mi < 0) { if (itcsk == 1) { A[i1] *= TWO; t *= TWO; } } A[i] = t; j = 0; jdif = ipow(2,csfftc.mt - csfftc.mm + 1); kkl = Ke[1] - kdr; if (csfftc.mm > 1) { for (kk = csfftc.ks, _do0=DOCNT(kk,kkl,_do1 = csfftc.ks); _do0 > 0; kk += _do1, _do0--) { kki = i + kk; kii = kki + kdr; kkn = ni1 - kk; kin = kkn + kdr; j += jdif; wi = S[j]; jj = csfftc.nt - j; wr = S[jj]; t = A[kki] + A[kkn]; ts = A[kkn] - A[kki]; t1 = A[kin] - A[kii]; ts1 = A[kin] + A[kii]; if (itcsk > 2) { /* The sine-cosine transform must be computed * differently in the case of the sine transform * because the input data is stored differently. */ tp = wr*t - wi*t1; tpi = wi*t + wr*t1; A[kki] = tp - ts1; A[kkn] = -tp - ts1; A[kii] = tpi + ts; A[kin] = tpi - ts; } else { tp = wr*ts1 + wi*ts; tpi = wi*ts1 - wr*ts; A[kki] = t + tp; A[kkn] = t - tp; A[kii] = t1 + tpi; A[kin] = tpi - t1; } } } else if (csfftc.mm == 0) { goto L_120; } /* End of computing sine-cosine transform * */ if (mi < 0) goto L_140; ir = i; ii = ir + kdr; /* Compute a one-dimensional complex Fourier transform */ L_110: sfft( &A[ir], &A[ii], s ); if (mi < 0) goto L_80; if (mi != 0) goto L_140; L_120: if (itcsk == 1) goto L_140; if (itcsk == 3) { A[i] = ZERO; } else { /* Compute first and last elements of the cosine transform */ sum *= FOUR; t = TWO*A[i]; A[ni1] = t - sum; A[i] = t + sum; if (csfftc.mm >= 0) A[ni2] *= TWO; } if (csfftc.mm > 0) { /* Extra calculations required by sine and cosine transform */ j = 0; jdif /= 2; for (kk = kdr, _do2=DOCNT(kk,kkl,_do3 = kdr); _do2 > 0; kk += _do3, _do2--) { kki = i + kk; kkn = ni1 - kk; j += jdif; wi = TWO*S[j]; t = A[kki] + A[kkn]; ts = A[kki] - A[kkn]; if (itcsk != 2) { t /= wi; } else { ts /= wi; } A[kki] = t + ts; A[kkn] = t - ts; } } /* Logic for deciding which one-dimensional transform to do next */ L_140: j = 0; L_150: j += 1; if (j == k) { j += 1; i += Nf[j] - Nf[j - 1]; } if (j > ndd) goto L_170; Mu[j] += Nf[j]; if (Mu[j] >= Nf[j + 1]) goto L_150; i += Nf[1]; j -= 1; if (j != 0) goto L_40; goto L_60; /* Set for Fourier analysis */ L_160: ii = i; ir = ii + kdr; goto L_110; L_170: ; } /* End of K loop */ return; /* Fatal error, default is to stop in IERM1 */ L_200: ierm1( "STCST", 4, 2, "Require 0 .le. max(M(K)) .le. 31", "M", csfftc.mm, '.' ); *ms = -2; return; } /* end of function */