C   These two programs set up the complex exponential table (PREFFT) and
C   compute the discrete-time Fourier series of an array of complex data
C   samples using a decimation-in-frequency fast Fourier transform (FFT)
C   algorithm.
C
      SUBROUTINE PREFFT (N,MODE,NEXP,W)
C
C   Input Parameters:
C
C     N     - Number of data samples to be processed (integer-must be a
C             power of two)
C     MODE  - Set to 0 for discrete-time Fourier series (Eq. 2.C.1) or
C             1 for inverse (Eq. 2.C.2)
C
C   Output Parameters:
C
C     NEXP  - Indicates power-of-2 exponent such that N=2**NEXP .
C             Will be set to -1 to indicate error condition if N
C             is not a power of 2 (this integer used by sub. FFT)
C     W     - Complex exponential array
C
C   Notes:
C
C     External array W must be dimensioned .GE. N by calling program.
C
      COMPLEX W(1),C1,C2
      NEXP=1
5     NT=2**NEXP
      IF (NT .GE. N)  GO TO 10
      NEXP=NEXP+1
      GO TO 5
10    IF (NT .EQ. N)  GO TO 15
      NEXP=-1                 ! Error:  N is not a power of 2
      RETURN
15    S=8.*ATAN(1.)/FLOAT(NT)
      C1=CMPLX(COS(S),-SIN(S))
      IF (MODE .NE. 0)  C1=CONJG(C1)
      C2=(1.,0.)
      DO 20 K=1,NT
        W(K)=C2
20      C2=C2*C1
      RETURN
      END


      SUBROUTINE FFT (N,MODE,T,NEXP,W,X)
C
C   Input Parameters:
C
C     N,MODE,NEXP,W - See parameter list for subroutine PREFFT
C     T             - Sample interval in seconds
C     X             - Array of N complex data samples, X(1) to X(N)
C
C   Output Parameters:
C
C     X - N complex transform values replace original data samples
C         indexed from k=1 to k=N, representing the frequencies
C         (k-1)/NT hertz
C
C   Notes:
C
C     External array X must be dimensioned .GE. N by calling program.
C
      COMPLEX X(1),W(1),C1,C2
      MM=1
      LL=N
      DO 70 K=1,NEXP
        NN=LL/2
        JJ=MM+1
        DO 40 I=1,N,LL
          KK=I+NN
          C1=X(I)+X(KK)
          X(KK)=X(I)-X(KK)
40        X(I)=C1
        IF (NN .EQ. 1) GO TO 70
        DO 60 J=2,NN
          C2=W(JJ)
          DO 50 I=J,N,LL
            KK=I+NN
            C1=X(I)+X(KK)
            X(KK)=(X(I)-X(KK))*C2
50          X(I)=C1
60        JJ=JJ+MM
        LL=NN
        MM=MM*2
70      CONTINUE
      NV2=N/2
      NM1=N-1
      J=1
      DO 90 I=1,NM1
        IF (I .GE. J)  GO TO 80
        C1=X(J)
        X(J)=X(I)
        X(I)=C1
80      K=NV2
85      IF (K .GE. J)  GO TO 90
        J=J-K
        K=K/2
        GO TO 85
90      J=J+K
      IF (MODE .EQ. 0)  S=T
      IF (MODE .NE. 0)  S=1./(T*FLOAT(N))
      DO 100 I=1,N
100     X(I)=X(I)*S
      RETURN
      END
