C
C-----------------------------------------------------------------------
C MAIN PROGRAM: FOURSUBT
C AUTHOR:       C. M. RADER
C               MIT LINCOLN LABORATORY, LEXINGTON, MA 02173
C
C INPUT:        FUNCTION IS GENERATED BY THE PROGRAM TO TEST
C               SUBROUTINE FOUREA
C-----------------------------------------------------------------------
C
      COMPLEX W, C(32), D, E, B(32), QB(32), A
C
C SEQUENCE LENGTH N=2**MU; USES MU=5; THUS N=32
C METHOD IS TO COMPUTE DFT OF KNOWN FUNCTION; A**I, I=0,1,...,N-1
C OUTPUT OF PROGRAM IS LARGEST DIFFERENCE BETWEEN DFT COMPUTED TWO WAYS.
C WHEN RUN ON AN IBM 370, THE MAX DIFFERENCE WAS 0.389E-04
C WHEN RUN ON A HONEYWELL 6080N, THE MAX DIFFERENCE WAS 0.238E-06
C A MAX DIFFERENCE LARGE COMPARED TO THE ACCURACY OF THE COMPUTER
C WOULD PROBABLY INDICATE A PROGRAM ERROR.
C
      MU = 5
C
C SET UP MACHINE CONSTANT
C
      IOUTD = I1MACH(2)
      NN = 2**MU
      TPI = 8.*ATAN(1.)
      TPION = TPI/FLOAT(NN)
      W = CMPLX(COS(TPION),-SIN(TPION))
C
C GENERATE A**K AS TEST FUNCTION
C
      A = (.9,.3)
      B(1) = (1.,0.)
      QB(1) = B(1)
      DO 10 K=2,NN
        B(K) = A**(K-1)
        QB(K) = B(K)
  10  CONTINUE
C
C PRINT COMPLEX INPUT SEQUENCE
C
      WRITE (IOUTD,9999)
9999  FORMAT (1H1//////24H COMPLEX INPUT SEQUENCE )
      WRITE (IOUTD,9998) (I,QB(I),I=1,NN)
9998  FORMAT (2(2X, 1H(, I3, 1H), 2E14.6))
C
C B(1) CONTAINS A**0; B(K) CONTAINS A**K-1; ETC
C
C COMPUTE DFT OF B IN CLOSED FORM
C
      D = (1.,0.) - A**NN
      DO 20 K=1,NN
        E = (1.,0.) - A*W**(K-1)
        C(K) = D/E
  20  CONTINUE
C
C DFT OF B IS (1-A**NN)/(1-AW**K)
C
C NOW COMPUTE DFT OF B USING FOUREA
C
      CALL FOUREA(B, NN, -1)
C
C PRINT FOUREA DFT AND THEORETICAL DFT
C
      WRITE (IOUTD,9997)
9997  FORMAT (/12H FOUREA DFT )
      WRITE (IOUTD,9998) (I,B(I),I=1,NN)
      WRITE (IOUTD,9996)
9996  FORMAT (/17H THEORETICAL DFT )
      WRITE (IOUTD,9998) (I,C(I),I=1,NN)
C
C FIND MAX DIFFERENCE BETWEEN B AND C
C
      DD = 0.
      DO 30 I=1,NN
        DE = CABS(C(I)-B(I))
        IF (DD.GT.DE) GO TO 30
        DD = DE
  30  CONTINUE
C
C COMPUTE INVERSE DFT OF C
C
      CALL FOUREA(C, NN, 1)
C
C PRINT INVERSE DFT
C
      WRITE (IOUTD,9995)
9995  FORMAT (1H1//////20H FOUREA INVERSE DFT )
      WRITE (IOUTD,9998) (I,C(I),I=1,NN)
C
C COMPUTE MAX. DIFF. BETWEEN INPUT AND INVERSE OF THEOR. DFT
C
      GG = 0.
      DO 40 I=1,NN
        GG1 = CABS(QB(I)-C(I))
        IF (GG1.LE.GG) GO TO 40
        GG = GG1
  40  CONTINUE
C
C PRINT MAXIMUM DIFFERENCE
C
      WRITE (IOUTD,9994) DD
      WRITE (IOUTD,9993) GG
9994  FORMAT (/42H MAX DIFF BETWEEN THEOR AND FOUREA DFT IS , E12.3)
9993  FORMAT (/51H MAX DIFF BETWEEN ORIGINAL DATA AND INVERSE DFT IS ,
     *    E12.3)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FOUREA
C PERFORMS COOLEY-TUKEY FAST FOURIER TRANSFORM
C-----------------------------------------------------------------------
C
      SUBROUTINE FOUREA(DATA, N, ISI)
C
C THE COOLEY-TUKEY FAST FOURIER TRANSFORM IN ANSI FORTRAN
C
C DATA IS A ONE-DIMENSIONAL COMPLEX ARRAY WHOSE LENGTH, N, IS A
C POWER OF TWO.  ISI IS +1 FOR AN INVERSE TRANSFORM AND -1 FOR A
C FORWARD TRANSFORM.  TRANSFORM VALUES ARE RETURNED IN THE INPUT
C ARRAY, REPLACING THE INPUT.
C TRANSFORM(J)=SUM(DATA(I)*W**((I-1)*(J-1))), WHERE I AND J RUN
C FROM 1 TO N AND W = EXP (ISI*2*PI*SQRT(-1)/N).  PROGRAM ALSO
C COMPUTES INVERSE TRANSFORM, FOR WHICH THE DEFINING EXPRESSION
C IS INVTR(J)=(1/N)*SUM(DATA(I)*W**((I-1)*(J-1))).
C RUNNING TIME IS PROPORTIONAL TO N*LOG2(N), RATHER THAN TO THE
C CLASSICAL N**2.
C AFTER PROGRAM BY BRENNER, JUNE 1967. THIS IS A VERY SHORT VERSION
C OF THE FFT AND IS INTENDED MAINLY FOR DEMONSTRATION. PROGRAMS
C ARE AVAILABLE IN THIS COLLECTION WHICH RUN FASTER AND ARE NOT
C RESTRICTED TO POWERS OF 2 OR TO ONE-DIMENSIONAL ARRAYS.
C SEE -- IEEE TRANS AUDIO (JUNE 1967), SPECIAL ISSUE ON FFT.
C
      COMPLEX DATA(1)
      COMPLEX TEMP, W
      IOUTD = I1MACH(2)
C
C CHECK FOR POWER OF TWO UP TO 15
C
      NN = 1
      DO 10 I=1,15
        M = I
        NN = NN*2
        IF (NN.EQ.N) GO TO 20
  10  CONTINUE
      WRITE (IOUTD,9999)
9999  FORMAT (30H N NOT A POWER OF 2 FOR FOUREA)
      STOP
  20  CONTINUE
C
      PI = 4.*ATAN(1.)
      FN = N
C
C THIS SECTION PUTS DATA IN BIT-REVERSED ORDER
C
      J = 1
      DO 80 I=1,N
C
C AT THIS POINT, I AND J ARE A BIT REVERSED PAIR (EXCEPT FOR THE
C DISPLACEMENT OF +1)
C
        IF (I-J) 30, 40, 40
C
C EXCHANGE DATA(I) WITH DATA(J) IF I.LT.J.
C
  30    TEMP = DATA(J)
        DATA(J) = DATA(I)
        DATA(I) = TEMP
C
C IMPLEMENT J=J+1, BIT-REVERSED COUNTER
C
  40    M = N/2
  50    IF (J-M) 70, 70, 60
  60    J = J - M
        M = (M+1)/2
        GO TO 50
  70    J = J + M
  80  CONTINUE
C
C NOW COMPUTE THE BUTTERFLIES
C
      MMAX = 1
  90  IF (MMAX-N) 100, 130, 130
 100  ISTEP = 2*MMAX
      DO 120 M=1,MMAX
        THETA = PI*FLOAT(ISI*(M-1))/FLOAT(MMAX)
        W = CMPLX(COS(THETA),SIN(THETA))
        DO 110 I=M,N,ISTEP
          J = I + MMAX
          TEMP = W*DATA(J)
          DATA(J) = DATA(I) - TEMP
          DATA(I) = DATA(I) + TEMP
 110    CONTINUE
 120  CONTINUE
      MMAX = ISTEP
      GO TO 90
 130  IF (ISI) 160, 140, 140
C
C FOR INV TRANS -- ISI=1 -- MULTIPLY OUTPUT BY 1/N
C
 140  DO 150 I=1,N
        DATA(I) = DATA(I)/FN
 150  CONTINUE
 160  RETURN
      END