/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:04 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dfft.h" #include /* PARAMETER translations */ #define HALF .5e0 #define KEDIM 30 #define PI4 .7853981633974483096156608458198757e0 #define SPI4 .7071067811865475244008443621048490e0 /* 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*/ dfft( double ar[], double ai[], double s[]) { LOGICAL32 spcase; long int _d_l, _d_m, _do0, _do1, _do2, _do3, _do4, _do5, i, i1, i2, i3, ij, j, jdif, ji, ji2, jj, jr, k, ksi, l, l1, l4, lj, ll, mtc; double t, t1, t1i, t2, t2i, t3, t3i, theta, ti, tp, tp1, tp1i, tpi, wi1, wi2, wi3, wr1, wr2, wr3; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Ai = &ai[0] - 1; double *const Ar = &ar[0] - 1; long *const Ke = &cdfftc.ke[0] - 1; double *const S = &s[0] - 1; /* end of OFFSET VECTORS */ /*>> 1997-03-28 DFFT Krogh KEDIM increased, removed assigned go to's. *>> 1994-11-11 DFFT Krogh Declared all vars. *>> 1994-10-20 DFFT Krogh Changes to use M77CON *>> 1989-08-15 DFFT FTK -- Parameter constants given to more precision. *>> 1989-04-21 DFFT FTK & CLL * Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * * This subroutine computes one dimensional complex Fourier * transforms using the Cooley-Tukey algorithm. It is used by a * number of subroutines which compute Fourier transforms. * * Programmed by Fred T. Krogh at the Jet Propulsion Laboratory, * Pasadena, Calif. August 1, 1969. * Revised by Krogh at JPL -- January 19, 1988 -- For portability * */ /* Minimum dimensions are AR(ILAST), AI(ILAST), S(NT-1), where ILAST * and NT are defined in the common block below. * * AR and AI give the arrays used to hold the real and imaginary parts of * the Fourier coefficients and data. * S contains the sine table required in the calculations. * * ----------------------------------------------------------------- * Notes on COMMON, PARAMETER's, and local variables * * SPI4 = SIN(PI/4) = .5 * SQRT(2) * PI4 = PI / 4 * THE DIMENSION OF KE MUST BE AS LARGE AS THE MAXIMUM VALUE * PERMITTED FOR MM. * * WHERE * NEEDST= .TRUE. if the sine table must be computed. * MT = base 2 log(NT) * NT = number of entries in the sine table * MM = base 2 log(number of complex fourier coefficients) * KS = distance in memory between successive points. The i-th * coefficient, a(i) is given by AR((I+1)*KS)+AI((I+1)*KS)* * sqrt(-1), i=0, 1, ..., (2**MM)-1. * ILAST = KS * 2**MM * KE(L) = KS * 2**(MM-L), L=1, 2, ..., MM * ----------------------------------------------------------------- *--D replaces "?": ?FFT, C?FFTC * ----------------------------------------------------------------- */ /* Common variables */ /* Note that KEE(1) is equivalent to ILAST. */ /* ----------------------------------------------------------------- * */ if (cdfftc.needst) { /* Compute the sine table */ cdfftc.needst = FALSE; mtc = cdfftc.mt; cdfftc.nt = ipow(2,mtc); if (mtc > 0) { /* SET FOR L=1 */ j = cdfftc.nt; jj = j/2; S[jj] = SPI4; if (mtc > 1) { theta = PI4; for (l = 2; l <= mtc; l++) { theta *= HALF; k = j; j = jj; jj /= 2; /* At this point J = 2**(MT-L+1) and JJ = 2**(MT-L) */ S[jj] = sin( theta ); l1 = cdfftc.nt - jj; S[l1] = cos( theta ); ll = cdfftc.nt - k; if (ll >= j) { for (i = j, _do0=DOCNT(i,ll,_do1 = j); _do0 > 0; i += _do1, _do0--) { i1 = cdfftc.nt - i; i2 = i + jj; S[i2] = S[i]*S[l1] + S[jj]*S[i1]; } } } } } } else { /* Compute the transform * Scramble A * */ ij = 1; ji = 1; l1 = cdfftc.ks; if (cdfftc.mm > 1) { L_30: ij += l1; for (i = 1; i <= cdfftc.mm; i++) { ji += Ke[i]; Ke[i] = -Ke[i]; if (Ke[i] < 0) { if (ij < ji) { /* INTERCHANGE THE IJ-TH COEFFICIENT WITH THE JI-TH */ t = Ar[ij]; Ar[ij] = Ar[ji]; Ar[ji] = t; t = Ai[ij]; Ai[ij] = Ai[ji]; Ai[ji] = t; } goto L_30; } } } /* END OF SCRAMBLING A */ jdif = cdfftc.nt; if ((cdfftc.mm%2) != 0) { /* SPECIAL CASE -- L = 1, MM ODD (RADIX 2 ALGORITHM) */ l1 += l1; for (i = 1, _do2=DOCNT(i,cdfftc.ilast,_do3 = l1); _do2 > 0; i += _do3, _do2--) { ksi = i + cdfftc.ks; t = Ar[i]; Ar[i] = t + Ar[ksi]; Ar[ksi] = t - Ar[ksi]; t = Ai[i]; Ai[i] = t + Ai[ksi]; Ai[ksi] = t - Ai[ksi]; } jdif /= 2; } for (l = 2; l <= cdfftc.mm; l += 2) { l4 = 4*l1; lj = l1/2; spcase = TRUE; j = 0; ji = 0; /* START OF I LOOP (RADIX 4 ALGORITHM) */ L_60: ij = j + 1; for (i = ij, _do4=DOCNT(i,cdfftc.ilast,_do5 = l4); _do4 > 0; i += _do5, _do4--) { i1 = i + l1; i2 = i1 + l1; i3 = i2 + l1; if (spcase) { if (j == 0) { /* SPECIAL CASE -- J = 0 */ t = Ar[i] + Ar[i1]; t1 = Ar[i] - Ar[i1]; ti = Ai[i] + Ai[i1]; t1i = Ai[i] - Ai[i1]; t2 = Ar[i2] + Ar[i3]; t3 = Ar[i2] - Ar[i3]; t2i = Ai[i2] + Ai[i3]; t3i = Ai[i2] - Ai[i3]; goto L_110; } /* SPECIAL CASE -- J = L1 / 2 */ t2 = SPI4*Ar[i2]; t2i = SPI4*Ai[i2]; t3 = SPI4*Ar[i3]; t3i = SPI4*Ai[i3]; tp = t2 - t2i; tpi = t2 + t2i; tp1 = -t3 - t3i; tp1i = t3 - t3i; t = Ar[i] - Ai[i1]; t1 = Ar[i] + Ai[i1]; ti = Ai[i] + Ar[i1]; t1i = Ai[i] - Ar[i1]; } else { /* USUAL CASE -- J .NE. 0 AND J .NE. L1 / 2 * * WRK AND WIK (K = 1, 2, 3) ARE THE REAL AND IMAGINARY PART * RESP. OF EXP(SQRT(-1) * PI * K*(2**(-L-MOD(MM, 2)))*J/KS) * =EXP(SQRT(-1) * PI * K * (J / (4 * L1))) * */ t2 = wr2*Ar[i1] - wi2*Ai[i1]; t2i = wi2*Ar[i1] + wr2*Ai[i1]; t = Ar[i] + t2; t1 = Ar[i] - t2; ti = Ai[i] + t2i; t1i = Ai[i] - t2i; tp = wr1*Ar[i2] - wi1*Ai[i2]; tpi = wi1*Ar[i2] + wr1*Ai[i2]; tp1 = wr3*Ar[i3] - wi3*Ai[i3]; tp1i = wi3*Ar[i3] + wr3*Ai[i3]; } t2 = tp + tp1; t3 = tp - tp1; t2i = tpi + tp1i; t3i = tpi - tp1i; L_110: Ar[i] = t + t2; Ai[i] = ti + t2i; Ar[i1] = t1 - t3i; Ai[i1] = t1i + t3; Ar[i2] = t - t2; Ai[i2] = ti - t2i; Ar[i3] = t1 + t3i; Ai[i3] = t1i - t3; } /* END OF I LOOP * */ if (j <= lj) { if (spcase) { if (j != 0) goto L_130; j = cdfftc.ks; spcase = FALSE; } else { j = l1 - j; /* COMPUTE WR-S AND WI-S FOR J REPLACED BY L1 - J */ t = wi1; wi1 = wr1; wr1 = t; wr2 = -wr2; t = -wi3; wi3 = -wr3; wr3 = t; goto L_60; } } else { j = l1 - j + cdfftc.ks; } if (j < lj) { /* COMPUTE WR-S AND WI-S */ ji += jdif; jr = cdfftc.nt - ji; wr1 = S[jr]; wi1 = S[ji]; ji2 = ji + ji; wi2 = S[ji2]; jr = cdfftc.nt - ji2; wr2 = S[jr]; ji2 += ji; if (ji2 <= cdfftc.nt) { jr = cdfftc.nt - ji2; wr3 = S[jr]; wi3 = S[ji2]; goto L_60; } jr = ji2 - cdfftc.nt; ji2 = cdfftc.nt - jr; wi3 = S[ji2]; wr3 = -S[jr]; goto L_60; } else if (j == lj) { spcase = TRUE; goto L_60; } /* END OF COMPUTING WR-S AND WI-S * */ L_130: l1 = l4; jdif /= 4; } /* END OF L LOOP */ } return; } /* end of function */