Fast Fourier Transform Algorithms

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: FASTMAIN  -  FAST FOURIER TRANSFORMS
C AUTHORS:      G. D. BERGLAND AND M. T. DOLAN
C               BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974
C
C INPUT:        THE PROGRAM CALLS ON A RANDOM NUMBER
C               GENERATOR FOR INPUT AND CHECKS DFT AND
C               IDFT WITH A 32-POINT SEQUENCE
C-----------------------------------------------------------------------
C
      DIMENSION X(32), Y(32), B(34)
C
C GENERATE RANDOM NUMBERS AND STORE ARRAY IN B SO
C THE SAME SEQUENCE CAN BE USED IN ALL TESTS.
C NOTE THAT B IS DIMENSIONED TO SIZE N+2.
C
C IW IS A MACHINE DEPENDENT WRITE DEVICE NUMBER
C
      IW = I1MACH(2)
C
      DO 10 I=1,32
        X(I) = UNI(0)
        B(I) = X(I)
  10  CONTINUE
      M = 5
      N = 2**M
      NP1 = N + 1
      NP2 = N + 2
      KNT = 1
C
C TEST FAST-FSST THEN FFA-FFS
C
      WRITE (IW,9999)
  20  WRITE (IW,9998) (B(I),I=1,N)
      IF (KNT.EQ.1) CALL FAST(B, N)
      IF (KNT.EQ.2) CALL FFA(B, N)
      WRITE (IW,9997) (B(I),I=1,NP1,2)
      WRITE (IW,9996) (B(I),I=2,NP2,2)
      IF (KNT.EQ.1) CALL FSST(B, N)
      IF (KNT.EQ.2) CALL FFS(B, N)
      WRITE (IW,9995) (B(I),I=1,N)
      KNT = KNT + 1
      IF (KNT.EQ.3) GO TO 40
C
      WRITE (IW,9994)
      DO 30 I=1,N
        B(I) = X(I)
  30  CONTINUE
      GO TO 20
C
C TEST FFT842 WITH REAL INPUT THEN COMPLEX
C
  40  WRITE (IW,9993)
      DO 50 I=1,N
        B(I) = X(I)
        Y(I) = 0.
  50  CONTINUE
  60  WRITE (IW,9992) (B(I),I=1,N)
      WRITE (IW,9991) (Y(I),I=1,N)
      CALL FFT842(0, N, B, Y)
      WRITE (IW,9997) (B(I),I=1,N)
      WRITE (IW,9996) (Y(I),I=1,N)
      CALL FFT842(1, N, B, Y)
      WRITE (IW,9990) (B(I),I=1,N)
      WRITE (IW,9989) (Y(I),I=1,N)
      KNT = KNT + 1
      IF (KNT.EQ.5) GO TO 80
C
      WRITE (IW,9988)
      DO 70 I=1,N
        B(I) = X(I)
        Y(I) = UNI(0)
  70  CONTINUE
      GO TO 60
C
9999  FORMAT (19H1TEST FAST AND FSST)
9998  FORMAT (20H0REAL INPUT SEQUENCE/(4E17.8))
9997  FORMAT (29H0REAL COMPONENTS OF TRANSFORM/(4E17.8))
9996  FORMAT (29H0IMAG COMPONENTS OF TRANSFORM/(4E17.8))
9995  FORMAT (23H0REAL INVERSE TRANSFORM/(4E17.8))
9994  FORMAT (17H1TEST FFA AND FFS)
9993  FORMAT (37H1TEST FFT842 WITH REAL INPUT SEQUENCE/(4E17.8))
9992  FORMAT (34H0REAL COMPONENTS OF INPUT SEQUENCE/(4E17.8))
9991  FORMAT (34H0IMAG COMPONENTS OF INPUT SEQUENCE/(4E17.8))
9990  FORMAT (37H0REAL COMPONENTS OF INVERSE TRANSFORM/(4E17.8))
9989  FORMAT (37H0IMAG COMPONENTS OF INVERSE TRANSFORM/(4E17.8))
9988  FORMAT (40H1TEST FFT842 WITH COMPLEX INPUT SEQUENCE)
  80  STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FAST
C REPLACES THE REAL VECTOR B(K), FOR K=1,2,...,N,
C WITH ITS FINITE DISCRETE FOURIER TRANSFORM
C-----------------------------------------------------------------------
C
      SUBROUTINE FAST(B, N)
C
C THE DC TERM IS RETURNED IN LOCATION B(1) WITH B(2) SET TO 0.
C THEREAFTER THE JTH HARMONIC IS RETURNED AS A COMPLEX
C NUMBER STORED AS  B(2*J+1) + I B(2*J+2).
C THE N/2 HARMONIC IS RETURNED IN B(N+1) WITH B(N+2) SET TO 0.
C HENCE, B MUST BE DIMENSIONED TO SIZE N+2.
C THE SUBROUTINE IS CALLED AS  FAST(B,N) WHERE N=2**M AND
C B IS THE REAL ARRAY DESCRIBED ABOVE.
C
      DIMENSION B(2)
      COMMON /CONS/ PII, P7, P7TWO, C22, S22, PI2
C
C IW IS A MACHINE DEPENDENT WRITE DEVICE NUMBER
C
      IW = I1MACH(2)
C
      PII = 4.*ATAN(1.)
      PI8 = PII/8.
      P7 = 1./SQRT(2.)
      P7TWO = 2.*P7
      C22 = COS(PI8)
      S22 = SIN(PI8)
      PI2 = 2.*PII
      DO 10 I=1,15
        M = I
        NT = 2**I
        IF (N.EQ.NT) GO TO 20
  10  CONTINUE
      WRITE (IW,9999)
9999  FORMAT (33H N IS NOT A POWER OF TWO FOR FAST)
      STOP
  20  N4POW = M/2
C
C DO A RADIX 2 ITERATION FIRST IF ONE IS REQUIRED.
C
      IF (M-N4POW*2) 40, 40, 30
  30  NN = 2
      INT = N/NN
      CALL FR2TR(INT, B(1), B(INT+1))
      GO TO 50
  40  NN = 1
C
C PERFORM RADIX 4 ITERATIONS.
C
  50  IF (N4POW.EQ.0) GO TO 70
      DO 60 IT=1,N4POW
        NN = NN*4
        INT = N/NN
        CALL FR4TR(INT, NN, B(1), B(INT+1), B(2*INT+1), B(3*INT+1),
     *      B(1), B(INT+1), B(2*INT+1), B(3*INT+1))
  60  CONTINUE
C
C PERFORM IN-PLACE REORDERING.
C
  70  CALL FORD1(M, B)
      CALL FORD2(M, B)
      T = B(2)
      B(2) = 0.
      B(N+1) = T
      B(N+2) = 0.
      DO 80 IT=4,N,2
        B(IT) = -B(IT)
  80  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FSST
C FOURIER SYNTHESIS SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE FSST(B, N)
C
C THIS SUBROUTINE SYNTHESIZES THE REAL VECTOR B(K), FOR
C K=1,2,...,N, FROM THE FOURIER COEFFICIENTS STORED IN THE
C B ARRAY OF SIZE N+2.  THE DC TERM IS IN B(1) WITH B(2) EQUAL
C TO  0.  THE JTH HARMONIC IS STORED AS B(2*J+1) + I B(2*J+2).
C THE N/2 HARMONIC IS IN B(N+1) WITH B(N+2) EQUAL TO 0.
C THE SUBROUTINE IS CALLED AS FSST(B,N) WHERE N=2**M AND
C B IS THE REAL ARRAY DISCUSSED ABOVE.
C
      DIMENSION B(2)
      COMMON /CONST/ PII, P7, P7TWO, C22, S22, PI2
C
C IW IS A MACHINE DEPENDENT WRITE DEVICE NUMBER
C
      IW = I1MACH(2)
C
      PII = 4.*ATAN(1.)
      PI8 = PII/8.
      P7 = 1./SQRT(2.)
      P7TWO = 2.*P7
      C22 = COS(PI8)
      S22 = SIN(PI8)
      PI2 = 2.*PII
      DO 10 I=1,15
        M = I
        NT = 2**I
        IF (N.EQ.NT) GO TO 20
  10  CONTINUE
      WRITE (IW,9999)
9999  FORMAT (33H N IS NOT A POWER OF TWO FOR FSST)
      STOP
  20  B(2) = B(N+1)
      DO 30 I=4,N,2
        B(I) = -B(I)
  30  CONTINUE
C
C SCALE THE INPUT BY N
C
      DO 40 I=1,N
        B(I) = B(I)/FLOAT(N)
  40  CONTINUE
      N4POW = M/2
C
C SCRAMBLE THE INPUTS
C
      CALL FORD2(M, B)
      CALL FORD1(M, B)
C
      IF (N4POW.EQ.0) GO TO 60
      NN = 4*N
      DO 50 IT=1,N4POW
        NN = NN/4
        INT = N/NN
        CALL FR4SYN(INT, NN, B(1), B(INT+1), B(2*INT+1), B(3*INT+1),
     *      B(1), B(INT+1), B(2*INT+1), B(3*INT+1))
  50  CONTINUE
C
C DO A RADIX 2 ITERATION IF ONE IS REQUIRED
C
  60  IF (M-N4POW*2) 80, 80, 70
  70  INT = N/2
      CALL FR2TR(INT, B(1), B(INT+1))
  80  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FR2TR
C RADIX 2 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE FR2TR(INT, B0, B1)
      DIMENSION B0(2), B1(2)
      DO 10 K=1,INT
        T = B0(K) + B1(K)
        B1(K) = B0(K) - B1(K)
        B0(K) = T
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FR4TR
C RADIX 4 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE FR4TR(INT, NN, B0, B1, B2, B3, B4, B5, B6, B7)
      DIMENSION L(15), B0(2), B1(2), B2(2), B3(2), B4(2), B5(2), B6(2),
     *    B7(2)
      COMMON /CONS/ PII, P7, P7TWO, C22, S22, PI2
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
C
C JTHET IS A REVERSED BINARY COUNTER, JR STEPS TWO AT A TIME TO
C LOCATE THE REAL PARTS OF INTERMEDIATE RESULTS, AND JI LOCATES
C THE IMAGINARY PART CORRESPONDING TO JR.
C
      L(1) = NN/4
      DO 40 K=2,15
        IF (L(K-1)-2) 10, 20, 30
  10    L(K-1) = 2
  20    L(K) = 2
        GO TO 40
  30    L(K) = L(K-1)/2
  40  CONTINUE
C
      PIOVN = PII/FLOAT(NN)
      JI = 3
      JL = 2
      JR = 2
C
      DO 120 J1=2,L1,2
      DO 120 J2=J1,L2,L1
      DO 120 J3=J2,L3,L2
      DO 120 J4=J3,L4,L3
      DO 120 J5=J4,L5,L4
      DO 120 J6=J5,L6,L5
      DO 120 J7=J6,L7,L6
      DO 120 J8=J7,L8,L7
      DO 120 J9=J8,L9,L8
      DO 120 J10=J9,L10,L9
      DO 120 J11=J10,L11,L10
      DO 120 J12=J11,L12,L11
      DO 120 J13=J12,L13,L12
      DO 120 J14=J13,L14,L13
      DO 120 JTHET=J14,L15,L14
        TH2 = JTHET - 2
        IF (TH2) 50, 50, 90
  50    DO 60 K=1,INT
          T0 = B0(K) + B2(K)
          T1 = B1(K) + B3(K)
          B2(K) = B0(K) - B2(K)
          B3(K) = B1(K) - B3(K)
          B0(K) = T0 + T1
          B1(K) = T0 - T1
  60    CONTINUE
C
        IF (NN-4) 120, 120, 70
  70    K0 = INT*4 + 1
        KL = K0 + INT - 1
        DO 80 K=K0,KL
          PR = P7*(B1(K)-B3(K))
          PI = P7*(B1(K)+B3(K))
          B3(K) = B2(K) + PI
          B1(K) = PI - B2(K)
          B2(K) = B0(K) - PR
          B0(K) = B0(K) + PR
  80    CONTINUE
        GO TO 120
C
  90    ARG = TH2*PIOVN
        C1 = COS(ARG)
        S1 = SIN(ARG)
        C2 = C1**2 - S1**2
        S2 = C1*S1 + C1*S1
        C3 = C1*C2 - S1*S2
        S3 = C2*S1 + S2*C1
C
        INT4 = INT*4
        J0 = JR*INT4 + 1
        K0 = JI*INT4 + 1
        JLAST = J0 + INT - 1
        DO 100 J=J0,JLAST
          K = K0 + J - J0
          R1 = B1(J)*C1 - B5(K)*S1
          R5 = B1(J)*S1 + B5(K)*C1
          T2 = B2(J)*C2 - B6(K)*S2
          T6 = B2(J)*S2 + B6(K)*C2
          T3 = B3(J)*C3 - B7(K)*S3
          T7 = B3(J)*S3 + B7(K)*C3
          T0 = B0(J) + T2
          T4 = B4(K) + T6
          T2 = B0(J) - T2
          T6 = B4(K) - T6
          T1 = R1 + T3
          T5 = R5 + T7
          T3 = R1 - T3
          T7 = R5 - T7
          B0(J) = T0 + T1
          B7(K) = T4 + T5
          B6(K) = T0 - T1
          B1(J) = T5 - T4
          B2(J) = T2 - T7
          B5(K) = T6 + T3
          B4(K) = T2 + T7
          B3(J) = T3 - T6
 100    CONTINUE
C
        JR = JR + 2
        JI = JI - 2
        IF (JI-JL) 110, 110, 120
 110    JI = 2*JR - 1
        JL = JR
 120  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FR4SYN
C RADIX 4 SYNTHESIS
C-----------------------------------------------------------------------
C
C
      SUBROUTINE FR4SYN(INT, NN, B0, B1, B2, B3, B4, B5, B6, B7)
      DIMENSION L(15), B0(2), B1(2), B2(2), B3(2), B4(2), B5(2), B6(2),
     *    B7(2)
      COMMON /CONST/ PII, P7, P7TWO, C22, S22, PI2
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
C
      L(1) = NN/4
      DO 40 K=2,15
        IF (L(K-1)-2) 10, 20, 30
  10    L(K-1) = 2
  20    L(K) = 2
        GO TO 40
  30    L(K) = L(K-1)/2
  40  CONTINUE
C
      PIOVN = PII/FLOAT(NN)
      JI = 3
      JL = 2
      JR = 2
C
      DO 120 J1=2,L1,2
      DO 120 J2=J1,L2,L1
      DO 120 J3=J2,L3,L2
      DO 120 J4=J3,L4,L3
      DO 120 J5=J4,L5,L4
      DO 120 J6=J5,L6,L5
      DO 120 J7=J6,L7,L6
      DO 120 J8=J7,L8,L7
      DO 120 J9=J8,L9,L8
      DO 120 J10=J9,L10,L9
      DO 120 J11=J10,L11,L10
      DO 120 J12=J11,L12,L11
      DO 120 J13=J12,L13,L12
      DO 120 J14=J13,L14,L13
      DO 120 JTHET=J14,L15,L14
        TH2 = JTHET - 2
        IF (TH2) 50, 50, 90
  50    DO 60 K=1,INT
          T0 = B0(K) + B1(K)
          T1 = B0(K) - B1(K)
          T2 = B2(K)*2.0
          T3 = B3(K)*2.0
          B0(K) = T0 + T2
          B2(K) = T0 - T2
          B1(K) = T1 + T3
          B3(K) = T1 - T3
  60    CONTINUE
C
        IF (NN-4) 120, 120, 70
  70    K0 = INT*4 + 1
        KL = K0 + INT - 1
        DO 80 K=K0,KL
          T2 = B0(K) - B2(K)
          T3 = B1(K) + B3(K)
          B0(K) = (B0(K)+B2(K))*2.0
          B2(K) = (B3(K)-B1(K))*2.0
          B1(K) = (T2+T3)*P7TWO
          B3(K) = (T3-T2)*P7TWO
  80    CONTINUE
        GO TO 120
  90    ARG = TH2*PIOVN
        C1 = COS(ARG)
        S1 = -SIN(ARG)
        C2 = C1**2 - S1**2
        S2 = C1*S1 + C1*S1
        C3 = C1*C2 - S1*S2
        S3 = C2*S1 + S2*C1
C
        INT4 = INT*4
        J0 = JR*INT4 + 1
        K0 = JI*INT4 + 1
        JLAST = J0 + INT - 1
        DO 100 J=J0,JLAST
          K = K0 + J - J0
          T0 = B0(J) + B6(K)
          T1 = B7(K) - B1(J)
          T2 = B0(J) - B6(K)
          T3 = B7(K) + B1(J)
          T4 = B2(J) + B4(K)
          T5 = B5(K) - B3(J)
          T6 = B5(K) + B3(J)
          T7 = B4(K) - B2(J)
          B0(J) = T0 + T4
          B4(K) = T1 + T5
          B1(J) = (T2+T6)*C1 - (T3+T7)*S1
          B5(K) = (T2+T6)*S1 + (T3+T7)*C1
          B2(J) = (T0-T4)*C2 - (T1-T5)*S2
          B6(K) = (T0-T4)*S2 + (T1-T5)*C2
          B3(J) = (T2-T6)*C3 - (T3-T7)*S3
          B7(K) = (T2-T6)*S3 + (T3-T7)*C3
 100    CONTINUE
        JR = JR + 2
        JI = JI - 2
        IF (JI-JL) 110, 110, 120
 110    JI = 2*JR - 1
        JL = JR
 120  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FORD1
C IN-PLACE REORDERING SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE FORD1(M, B)
      DIMENSION B(2)
C
      K = 4
      KL = 2
      N = 2**M
      DO 40 J=4,N,2
        IF (K-J) 20, 20, 10
  10    T = B(J)
        B(J) = B(K)
        B(K) = T
  20    K = K - 2
        IF (K-KL) 30, 30, 40
  30    K = 2*J
        KL = J
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FORD2
C IN-PLACE REORDERING SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE FORD2(M, B)
      DIMENSION L(15), B(2)
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
      N = 2**M
      L(1) = N
      DO 10 K=2,M
        L(K) = L(K-1)/2
  10  CONTINUE
      DO 20 K=M,14
        L(K+1) = 2
  20  CONTINUE
      IJ = 2
      DO 40 J1=2,L1,2
      DO 40 J2=J1,L2,L1
      DO 40 J3=J2,L3,L2
      DO 40 J4=J3,L4,L3
      DO 40 J5=J4,L5,L4
      DO 40 J6=J5,L6,L5
      DO 40 J7=J6,L7,L6
      DO 40 J8=J7,L8,L7
      DO 40 J9=J8,L9,L8
      DO 40 J10=J9,L10,L9
      DO 40 J11=J10,L11,L10
      DO 40 J12=J11,L12,L11
      DO 40 J13=J12,L13,L12
      DO 40 J14=J13,L14,L13
      DO 40 JI=J14,L15,L14
        IF (IJ-JI) 30, 40, 40
  30    T = B(IJ-1)
        B(IJ-1) = B(JI-1)
        B(JI-1) = T
        T = B(IJ)
        B(IJ) = B(JI)
        B(JI) = T
  40    IJ = IJ + 2
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FFA
C FAST FOURIER ANALYSIS SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE FFA(B, NFFT)
C
C THIS SUBROUTINE REPLACES THE REAL VECTOR B(K),  (K=1,2,...,N),
C WITH ITS FINITE DISCRETE FOURIER TRANSFORM.  THE DC TERM IS
C RETURNED IN LOCATION B(1) WITH B(2) SET TO 0.  THEREAFTER, THE
C JTH HARMONIC IS RETURNED AS A COMPLEX NUMBER STORED AS
C B(2*J+1) + I B(2*J+2).  NOTE THAT THE N/2 HARMONIC IS RETURNED
C IN B(N+1) WITH B(N+2) SET TO 0.  HENCE, B MUST BE DIMENSIONED
C TO SIZE N+2.
C SUBROUTINE IS CALLED AS FFA (B,N) WHERE N=2**M AND B IS AN
C N TERM REAL ARRAY.  A REAL-VALUED, RADIX 8  ALGORITHM IS USED
C WITH IN-PLACE REORDERING AND THE TRIG FUNCTIONS ARE COMPUTED AS
C NEEDED.
C
      DIMENSION B(2)
      COMMON /CON/ PII, P7, P7TWO, C22, S22, PI2
C
C IW IS A MACHINE DEPENDENT WRITE DEVICE NUMBER
C
      IW = I1MACH(2)
C
      PII = 4.*ATAN(1.)
      PI8 = PII/8.
      P7 = 1./SQRT(2.)
      P7TWO = 2.*P7
      C22 = COS(PI8)
      S22 = SIN(PI8)
      PI2 = 2.*PII
      N = 1
      DO 10 I=1,15
        M = I
        N = N*2
        IF (N.EQ.NFFT) GO TO 20
  10  CONTINUE
      WRITE (IW,9999)
9999  FORMAT (30H NFFT NOT A POWER OF 2 FOR FFA)
      STOP
  20  CONTINUE
      N8POW = M/3
C
C DO A RADIX 2 OR RADIX 4 ITERATION FIRST IF ONE IS REQUIRED
C
      IF (M-N8POW*3-1) 50, 40, 30
  30  NN = 4
      INT = N/NN
      CALL R4TR(INT, B(1), B(INT+1), B(2*INT+1), B(3*INT+1))
      GO TO 60
  40  NN = 2
      INT = N/NN
      CALL R2TR(INT, B(1), B(INT+1))
      GO TO 60
  50  NN = 1
C
C PERFORM RADIX 8 ITERATIONS
C
  60  IF (N8POW) 90, 90, 70
  70  DO 80 IT=1,N8POW
        NN = NN*8
        INT = N/NN
        CALL R8TR(INT, NN, B(1), B(INT+1), B(2*INT+1), B(3*INT+1),
     *      B(4*INT+1), B(5*INT+1), B(6*INT+1), B(7*INT+1), B(1),
     *      B(INT+1), B(2*INT+1), B(3*INT+1), B(4*INT+1), B(5*INT+1),
     *      B(6*INT+1), B(7*INT+1))
  80  CONTINUE
C
C PERFORM IN-PLACE REORDERING
C
  90  CALL ORD1(M, B)
      CALL ORD2(M, B)
      T = B(2)
      B(2) = 0.
      B(NFFT+1) = T
      B(NFFT+2) = 0.
      DO 100 I=4,NFFT,2
        B(I) = -B(I)
 100  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FFS
C FAST FOURIER SYNTHESIS SUBROUTINE
C RADIX 8-4-2
C-----------------------------------------------------------------------
C
      SUBROUTINE FFS(B, NFFT)
C
C THIS SUBROUTINE SYNTHESIZES THE REAL VECTOR B(K), WHERE
C K=1,2,...,N. THE INITIAL FOURIER COEFFICIENTS ARE PLACED IN
C THE B ARRAY OF SIZE N+2.  THE DC TERM IS IN B(1) WITH
C B(2) EQUAL TO 0.
C THE JTH HARMONIC IS STORED AS B(2*J+1) + I B(2*J+2).
C THE N/2 HARMONIC IS IN B(N+1) WITH B(N+2) EQUAL TO 0.
C THE SUBROUTINE IS CALLED AS FFS(B,N) WHERE N=2**M AND
C B IS THE N TERM REAL ARRAY DISCUSSED ABOVE.
C
      DIMENSION B(2)
      COMMON /CON1/ PII, P7, P7TWO, C22, S22, PI2
C
C IW IS A MACHINE DEPENDENT WRITE DEVICE NUMBER
C
      IW = I1MACH(2)
C
      PII = 4.*ATAN(1.)
      PI8 = PII/8.
      P7 = 1./SQRT(2.)
      P7TWO = 2.*P7
      C22 = COS(PI8)
      S22 = SIN(PI8)
      PI2 = 2.*PII
      N = 1
      DO 10 I=1,15
        M = I
        N = N*2
        IF (N.EQ.NFFT) GO TO 20
  10  CONTINUE
      WRITE (IW,9999)
9999  FORMAT (30H NFFT NOT A POWER OF 2 FOR FFS)
      STOP
  20  CONTINUE
      B(2) = B(NFFT+1)
      DO 30 I=1,NFFT
        B(I) = B(I)/FLOAT(NFFT)
  30  CONTINUE
      DO 40 I=4,NFFT,2
        B(I) = -B(I)
  40  CONTINUE
      N8POW = M/3
C
C REORDER THE INPUT FOURIER COEFFICIENTS
C
      CALL ORD2(M, B)
      CALL ORD1(M, B)
C
      IF (N8POW.EQ.0) GO TO 60
C
C PERFORM THE RADIX 8 ITERATIONS
C
      NN = N
      DO 50 IT=1,N8POW
        INT = N/NN
        CALL R8SYN(INT, NN, B, B(INT+1), B(2*INT+1), B(3*INT+1),
     *      B(4*INT+1), B(5*INT+1), B(6*INT+1), B(7*INT+1), B(1),
     *      B(INT+1), B(2*INT+1), B(3*INT+1), B(4*INT+1), B(5*INT+1),
     *      B(6*INT+1), B(7*INT+1))
        NN = NN/8
  50  CONTINUE
C
C DO A RADIX 2 OR RADIX 4 ITERATION IF ONE IS REQUIRED
C
  60  IF (M-N8POW*3-1) 90, 80, 70
  70  INT = N/4
      CALL R4SYN(INT, B(1), B(INT+1), B(2*INT+1), B(3*INT+1))
      GO TO 90
  80  INT = N/2
      CALL R2TR(INT, B(1), B(INT+1))
  90  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R2TR
C RADIX 2 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
C
      SUBROUTINE R2TR(INT, B0, B1)
      DIMENSION B0(2), B1(2)
      DO 10 K=1,INT
        T = B0(K) + B1(K)
        B1(K) = B0(K) - B1(K)
        B0(K) = T
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R4TR
C RADIX 4 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE R4TR(INT, B0, B1, B2, B3)
      DIMENSION B0(2), B1(2), B2(2), B3(2)
      DO 10 K=1,INT
        R0 = B0(K) + B2(K)
        R1 = B1(K) + B3(K)
        B2(K) = B0(K) - B2(K)
        B3(K) = B1(K) - B3(K)
        B0(K) = R0 + R1
        B1(K) = R0 - R1
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: R8TR
C RADIX 8 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE R8TR(INT, NN, BR0, BR1, BR2, BR3, BR4, BR5, BR6, BR7,
     *    BI0, BI1, BI2, BI3, BI4, BI5, BI6, BI7)
      DIMENSION L(15), BR0(2), BR1(2), BR2(2), BR3(2), BR4(2), BR5(2),
     *    BR6(2), BR7(2), BI0(2), BI1(2), BI2(2), BI3(2), BI4(2),
     *    BI5(2), BI6(2), BI7(2)
      COMMON /CON/ PII, P7, P7TWO, C22, S22, PI2
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
C
C SET UP COUNTERS SUCH THAT JTHET STEPS THROUGH THE ARGUMENTS
C OF W, JR STEPS THROUGH STARTING LOCATIONS FOR THE REAL PART OF THE
C INTERMEDIATE RESULTS AND JI STEPS THROUGH STARTING LOCATIONS
C OF THE IMAGINARY PART OF THE INTERMEDIATE RESULTS.
C
      L(1) = NN/8
      DO 40 K=2,15
        IF (L(K-1)-2) 10, 20, 30
  10    L(K-1) = 2
  20    L(K) = 2
        GO TO 40
  30    L(K) = L(K-1)/2
  40  CONTINUE
      PIOVN = PII/FLOAT(NN)
      JI = 3
      JL = 2
      JR = 2
      DO 120 J1=2,L1,2
      DO 120 J2=J1,L2,L1
      DO 120 J3=J2,L3,L2
      DO 120 J4=J3,L4,L3
      DO 120 J5=J4,L5,L4
      DO 120 J6=J5,L6,L5
      DO 120 J7=J6,L7,L6
      DO 120 J8=J7,L8,L7
      DO 120 J9=J8,L9,L8
      DO 120 J10=J9,L10,L9
      DO 120 J11=J10,L11,L10
      DO 120 J12=J11,L12,L11
      DO 120 J13=J12,L13,L12
      DO 120 J14=J13,L14,L13
      DO 120 JTHET=J14,L15,L14
        TH2 = JTHET - 2
        IF (TH2) 50, 50, 90
  50    DO 60 K=1,INT
          T0 = BR0(K) + BR4(K)
          T1 = BR1(K) + BR5(K)
          T2 = BR2(K) + BR6(K)
          T3 = BR3(K) + BR7(K)
          T4 = BR0(K) - BR4(K)
          T5 = BR1(K) - BR5(K)
          T6 = BR2(K) - BR6(K)
          T7 = BR3(K) - BR7(K)
          BR2(K) = T0 - T2
          BR3(K) = T1 - T3
          T0 = T0 + T2
          T1 = T1 + T3
          BR0(K) = T0 + T1
          BR1(K) = T0 - T1
          PR = P7*(T5-T7)
          PI = P7*(T5+T7)
          BR4(K) = T4 + PR
          BR7(K) = T6 + PI
          BR6(K) = T4 - PR
          BR5(K) = PI - T6
  60    CONTINUE
        IF (NN-8) 120, 120, 70
  70    K0 = INT*8 + 1
        KL = K0 + INT - 1
        DO 80 K=K0,KL
          PR = P7*(BI2(K)-BI6(K))
          PI = P7*(BI2(K)+BI6(K))
          TR0 = BI0(K) + PR
          TI0 = BI4(K) + PI
          TR2 = BI0(K) - PR
          TI2 = BI4(K) - PI
          PR = P7*(BI3(K)-BI7(K))
          PI = P7*(BI3(K)+BI7(K))
          TR1 = BI1(K) + PR
          TI1 = BI5(K) + PI
          TR3 = BI1(K) - PR
          TI3 = BI5(K) - PI
          PR = TR1*C22 - TI1*S22
          PI = TI1*C22 + TR1*S22
          BI0(K) = TR0 + PR
          BI6(K) = TR0 - PR
          BI7(K) = TI0 + PI
          BI1(K) = PI - TI0
          PR = -TR3*S22 - TI3*C22
          PI = TR3*C22 - TI3*S22
          BI2(K) = TR2 + PR
          BI4(K) = TR2 - PR
          BI5(K) = TI2 + PI
          BI3(K) = PI - TI2
  80    CONTINUE
        GO TO 120
  90    ARG = TH2*PIOVN
        C1 = COS(ARG)
        S1 = SIN(ARG)
        C2 = C1**2 - S1**2
        S2 = C1*S1 + C1*S1
        C3 = C1*C2 - S1*S2
        S3 = C2*S1 + S2*C1
        C4 = C2**2 - S2**2
        S4 = C2*S2 + C2*S2
        C5 = C2*C3 - S2*S3
        S5 = C3*S2 + S3*C2
        C6 = C3**2 - S3**2
        S6 = C3*S3 + C3*S3
        C7 = C3*C4 - S3*S4
        S7 = C4*S3 + S4*C3
        INT8 = INT*8
        J0 = JR*INT8 + 1
        K0 = JI*INT8 + 1
        JLAST = J0 + INT - 1
        DO 100 J=J0,JLAST
          K = K0 + J - J0
          TR1 = BR1(J)*C1 - BI1(K)*S1
          TI1 = BR1(J)*S1 + BI1(K)*C1
          TR2 = BR2(J)*C2 - BI2(K)*S2
          TI2 = BR2(J)*S2 + BI2(K)*C2
          TR3 = BR3(J)*C3 - BI3(K)*S3
          TI3 = BR3(J)*S3 + BI3(K)*C3
          TR4 = BR4(J)*C4 - BI4(K)*S4
          TI4 = BR4(J)*S4 + BI4(K)*C4
          TR5 = BR5(J)*C5 - BI5(K)*S5
          TI5 = BR5(J)*S5 + BI5(K)*C5
          TR6 = BR6(J)*C6 - BI6(K)*S6
          TI6 = BR6(J)*S6 + BI6(K)*C6
          TR7 = BR7(J)*C7 - BI7(K)*S7
          TI7 = BR7(J)*S7 + BI7(K)*C7
C
          T0 = BR0(J) + TR4
          T1 = BI0(K) + TI4
          TR4 = BR0(J) - TR4
          TI4 = BI0(K) - TI4
          T2 = TR1 + TR5
          T3 = TI1 + TI5
          TR5 = TR1 - TR5
          TI5 = TI1 - TI5
          T4 = TR2 + TR6
          T5 = TI2 + TI6
          TR6 = TR2 - TR6
          TI6 = TI2 - TI6
          T6 = TR3 + TR7
          T7 = TI3 + TI7
          TR7 = TR3 - TR7
          TI7 = TI3 - TI7
C
          TR0 = T0 + T4
          TI0 = T1 + T5
          TR2 = T0 - T4
          TI2 = T1 - T5
          TR1 = T2 + T6
          TI1 = T3 + T7
          TR3 = T2 - T6
          TI3 = T3 - T7
          T0 = TR4 - TI6
          T1 = TI4 + TR6
          T4 = TR4 + TI6
          T5 = TI4 - TR6
          T2 = TR5 - TI7
          T3 = TI5 + TR7
          T6 = TR5 + TI7
          T7 = TI5 - TR7
          BR0(J) = TR0 + TR1
          BI7(K) = TI0 + TI1
          BI6(K) = TR0 - TR1
          BR1(J) = TI1 - TI0
          BR2(J) = TR2 - TI3
          BI5(K) = TI2 + TR3
          BI4(K) = TR2 + TI3
          BR3(J) = TR3 - TI2
          PR = P7*(T2-T3)
          PI = P7*(T2+T3)
          BR4(J) = T0 + PR
          BI3(K) = T1 + PI
          BI2(K) = T0 - PR
          BR5(J) = PI - T1
          PR = -P7*(T6+T7)
          PI = P7*(T6-T7)
          BR6(J) = T4 + PR
          BI1(K) = T5 + PI
          BI0(K) = T4 - PR
          BR7(J) = PI - T5
 100    CONTINUE
        JR = JR + 2
        JI = JI - 2
        IF (JI-JL) 110, 110, 120
 110    JI = 2*JR - 1
        JL = JR
 120  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R4SYN
C RADIX 4 SYNTHESIS
C-----------------------------------------------------------------------
C
      SUBROUTINE R4SYN(INT, B0, B1, B2, B3)
      DIMENSION B0(2), B1(2), B2(2), B3(2)
      DO 10 K=1,INT
        T0 = B0(K) + B1(K)
        T1 = B0(K) - B1(K)
        T2 = B2(K) + B2(K)
        T3 = B3(K) + B3(K)
        B0(K) = T0 + T2
        B2(K) = T0 - T2
        B1(K) = T1 + T3
        B3(K) = T1 - T3
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R8SYN
C RADIX 8 SYNTHESIS SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE R8SYN(INT, NN, BR0, BR1, BR2, BR3, BR4, BR5, BR6, BR7,
     *    BI0, BI1, BI2, BI3, BI4, BI5, BI6, BI7)
      DIMENSION L(15), BR0(2), BR1(2), BR2(2), BR3(2), BR4(2), BR5(2),
     *    BR6(2), BR7(2), BI0(2), BI1(2), BI2(2), BI3(2), BI4(2),
     *    BI5(2), BI6(2), BI7(2)
      COMMON /CON1/ PII, P7, P7TWO, C22, S22, PI2
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
      L(1) = NN/8
      DO 40 K=2,15
        IF (L(K-1)-2) 10, 20, 30
  10    L(K-1) = 2
  20    L(K) = 2
        GO TO 40
  30    L(K) = L(K-1)/2
  40  CONTINUE
      PIOVN = PII/FLOAT(NN)
      JI = 3
      JL = 2
      JR = 2
C
      DO 120 J1=2,L1,2
      DO 120 J2=J1,L2,L1
      DO 120 J3=J2,L3,L2
      DO 120 J4=J3,L4,L3
      DO 120 J5=J4,L5,L4
      DO 120 J6=J5,L6,L5
      DO 120 J7=J6,L7,L6
      DO 120 J8=J7,L8,L7
      DO 120 J9=J8,L9,L8
      DO 120 J10=J9,L10,L9
      DO 120 J11=J10,L11,L10
      DO 120 J12=J11,L12,L11
      DO 120 J13=J12,L13,L12
      DO 120 J14=J13,L14,L13
      DO 120 JTHET=J14,L15,L14
        TH2 = JTHET - 2
        IF (TH2) 50, 50, 90
  50    DO 60 K=1,INT
          T0 = BR0(K) + BR1(K)
          T1 = BR0(K) - BR1(K)
          T2 = BR2(K) + BR2(K)
          T3 = BR3(K) + BR3(K)
          T4 = BR4(K) + BR6(K)
          T6 = BR7(K) - BR5(K)
          T5 = BR4(K) - BR6(K)
          T7 = BR7(K) + BR5(K)
          PR = P7*(T7+T5)
          PI = P7*(T7-T5)
          TT0 = T0 + T2
          TT1 = T1 + T3
          T2 = T0 - T2
          T3 = T1 - T3
          T4 = T4 + T4
          T5 = PR + PR
          T6 = T6 + T6
          T7 = PI + PI
          BR0(K) = TT0 + T4
          BR1(K) = TT1 + T5
          BR2(K) = T2 + T6
          BR3(K) = T3 + T7
          BR4(K) = TT0 - T4
          BR5(K) = TT1 - T5
          BR6(K) = T2 - T6
          BR7(K) = T3 - T7
  60    CONTINUE
        IF (NN-8) 120, 120, 70
  70    K0 = INT*8 + 1
        KL = K0 + INT - 1
        DO 80 K=K0,KL
          T1 = BI0(K) + BI6(K)
          T2 = BI7(K) - BI1(K)
          T3 = BI0(K) - BI6(K)
          T4 = BI7(K) + BI1(K)
          PR = T3*C22 + T4*S22
          PI = T4*C22 - T3*S22
          T5 = BI2(K) + BI4(K)
          T6 = BI5(K) - BI3(K)
          T7 = BI2(K) - BI4(K)
          T8 = BI5(K) + BI3(K)
          RR = T8*C22 - T7*S22
          RI = -T8*S22 - T7*C22
          BI0(K) = (T1+T5) + (T1+T5)
          BI4(K) = (T2+T6) + (T2+T6)
          BI1(K) = (PR+RR) + (PR+RR)
          BI5(K) = (PI+RI) + (PI+RI)
          T5 = T1 - T5
          T6 = T2 - T6
          BI2(K) = P7TWO*(T6+T5)
          BI6(K) = P7TWO*(T6-T5)
          RR = PR - RR
          RI = PI - RI
          BI3(K) = P7TWO*(RI+RR)
          BI7(K) = P7TWO*(RI-RR)
  80    CONTINUE
        GO TO 120
  90    ARG = TH2*PIOVN
        C1 = COS(ARG)
        S1 = -SIN(ARG)
        C2 = C1**2 - S1**2
        S2 = C1*S1 + C1*S1
        C3 = C1*C2 - S1*S2
        S3 = C2*S1 + S2*C1
        C4 = C2**2 - S2**2
        S4 = C2*S2 + C2*S2
        C5 = C2*C3 - S2*S3
        S5 = C3*S2 + S3*C2
        C6 = C3**2 - S3**2
        S6 = C3*S3 + C3*S3
        C7 = C3*C4 - S3*S4
        S7 = C4*S3 + S4*C3
        INT8 = INT*8
        J0 = JR*INT8 + 1
        K0 = JI*INT8 + 1
        JLAST = J0 + INT - 1
        DO 100 J=J0,JLAST
          K = K0 + J - J0
          TR0 = BR0(J) + BI6(K)
          TI0 = BI7(K) - BR1(J)
          TR1 = BR0(J) - BI6(K)
          TI1 = BI7(K) + BR1(J)
          TR2 = BR2(J) + BI4(K)
          TI2 = BI5(K) - BR3(J)
          TR3 = BI5(K) + BR3(J)
          TI3 = BI4(K) - BR2(J)
          TR4 = BR4(J) + BI2(K)
          TI4 = BI3(K) - BR5(J)
          T0 = BR4(J) - BI2(K)
          T1 = BI3(K) + BR5(J)
          TR5 = P7*(T0+T1)
          TI5 = P7*(T1-T0)
          TR6 = BR6(J) + BI0(K)
          TI6 = BI1(K) - BR7(J)
          T0 = BR6(J) - BI0(K)
          T1 = BI1(K) + BR7(J)
          TR7 = -P7*(T0-T1)
          TI7 = -P7*(T1+T0)
          T0 = TR0 + TR2
          T1 = TI0 + TI2
          T2 = TR1 + TR3
          T3 = TI1 + TI3
          TR2 = TR0 - TR2
          TI2 = TI0 - TI2
          TR3 = TR1 - TR3
          TI3 = TI1 - TI3
          T4 = TR4 + TR6
          T5 = TI4 + TI6
          T6 = TR5 + TR7
          T7 = TI5 + TI7
          TTR6 = TI4 - TI6
          TI6 = TR6 - TR4
          TTR7 = TI5 - TI7
          TI7 = TR7 - TR5
          BR0(J) = T0 + T4
          BI0(K) = T1 + T5
          BR1(J) = C1*(T2+T6) - S1*(T3+T7)
          BI1(K) = C1*(T3+T7) + S1*(T2+T6)
          BR2(J) = C2*(TR2+TTR6) - S2*(TI2+TI6)
          BI2(K) = C2*(TI2+TI6) + S2*(TR2+TTR6)
          BR3(J) = C3*(TR3+TTR7) - S3*(TI3+TI7)
          BI3(K) = C3*(TI3+TI7) + S3*(TR3+TTR7)
          BR4(J) = C4*(T0-T4) - S4*(T1-T5)
          BI4(K) = C4*(T1-T5) + S4*(T0-T4)
          BR5(J) = C5*(T2-T6) - S5*(T3-T7)
          BI5(K) = C5*(T3-T7) + S5*(T2-T6)
          BR6(J) = C6*(TR2-TTR6) - S6*(TI2-TI6)
          BI6(K) = C6*(TI2-TI6) + S6*(TR2-TTR6)
          BR7(J) = C7*(TR3-TTR7) - S7*(TI3-TI7)
          BI7(K) = C7*(TI3-TI7) + S7*(TR3-TTR7)
 100    CONTINUE
        JR = JR + 2
        JI = JI - 2
        IF (JI-JL) 110, 110, 120
 110    JI = 2*JR - 1
        JL = JR
 120  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  ORD1
C IN-PLACE REORDERING SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE ORD1(M, B)
      DIMENSION B(2)
C
      K = 4
      KL = 2
      N = 2**M
      DO 40 J=4,N,2
        IF (K-J) 20, 20, 10
  10    T = B(J)
        B(J) = B(K)
        B(K) = T
  20    K = K - 2
        IF (K-KL) 30, 30, 40
  30    K = 2*J
        KL = J
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  ORD2
C IN-PLACE REORDERING SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE ORD2(M, B)
      DIMENSION L(15), B(2)
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
      N = 2**M
      L(1) = N
      DO 10 K=2,M
        L(K) = L(K-1)/2
  10  CONTINUE
      DO 20 K=M,14
        L(K+1) = 2
  20  CONTINUE
      IJ = 2
      DO 40 J1=2,L1,2
      DO 40 J2=J1,L2,L1
      DO 40 J3=J2,L3,L2
      DO 40 J4=J3,L4,L3
      DO 40 J5=J4,L5,L4
      DO 40 J6=J5,L6,L5
      DO 40 J7=J6,L7,L6
      DO 40 J8=J7,L8,L7
      DO 40 J9=J8,L9,L8
      DO 40 J10=J9,L10,L9
      DO 40 J11=J10,L11,L10
      DO 40 J12=J11,L12,L11
      DO 40 J13=J12,L13,L12
      DO 40 J14=J13,L14,L13
      DO 40 JI=J14,L15,L14
        IF (IJ-JI) 30, 40, 40
  30    T = B(IJ-1)
        B(IJ-1) = B(JI-1)
        B(JI-1) = T
        T = B(IJ)
        B(IJ) = B(JI)
        B(JI) = T
  40    IJ = IJ + 2
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FFT842
C FAST FOURIER TRANSFORM FOR N=2**M
C COMPLEX INPUT
C-----------------------------------------------------------------------
C
      SUBROUTINE FFT842(IN, N, X, Y)
C
C THIS PROGRAM REPLACES THE VECTOR Z=X+IY BY ITS  FINITE
C DISCRETE, COMPLEX FOURIER TRANSFORM IF IN=0.  THE INVERSE TRANSFORM
C IS CALCULATED FOR IN=1.  IT PERFORMS AS MANY BASE
C 8 ITERATIONS AS POSSIBLE AND THEN FINISHES WITH A BASE 4 ITERATION
C OR A BASE 2 ITERATION IF NEEDED.
C
C THE SUBROUTINE IS CALLED AS SUBROUTINE FFT842 (IN,N,X,Y).
C THE INTEGER N (A POWER OF 2), THE N REAL LOCATION ARRAY X, AND
C THE N REAL LOCATION ARRAY Y MUST BE SUPPLIED TO THE SUBROUTINE.
C
      DIMENSION X(2), Y(2), L(15)
      COMMON /CON2/ PI2, P7
      EQUIVALENCE (L15,L(1)), (L14,L(2)), (L13,L(3)), (L12,L(4)),
     *    (L11,L(5)), (L10,L(6)), (L9,L(7)), (L8,L(8)), (L7,L(9)),
     *    (L6,L(10)), (L5,L(11)), (L4,L(12)), (L3,L(13)), (L2,L(14)),
     *    (L1,L(15))
C
C
C IW IS A MACHINE DEPENDENT WRITE DEVICE NUMBER
C
      IW = I1MACH(2)
C
      PI2 = 8.*ATAN(1.)
      P7 = 1./SQRT(2.)
      DO 10 I=1,15
        M = I
        NT = 2**I
        IF (N.EQ.NT) GO TO 20
  10  CONTINUE
      WRITE (IW,9999)
9999  FORMAT (35H N IS NOT A POWER OF TWO FOR FFT842)
      STOP
  20  N2POW = M
      NTHPO = N
      FN = NTHPO
      IF (IN.EQ.1) GO TO 40
      DO 30 I=1,NTHPO
        Y(I) = -Y(I)
  30  CONTINUE
  40  N8POW = N2POW/3
      IF (N8POW.EQ.0) GO TO 60
C
C RADIX 8 PASSES,IF ANY.
C
      DO 50 IPASS=1,N8POW
        NXTLT = 2**(N2POW-3*IPASS)
        LENGT = 8*NXTLT
        CALL R8TX(NXTLT, NTHPO, LENGT, X(1), X(NXTLT+1), X(2*NXTLT+1),
     *      X(3*NXTLT+1), X(4*NXTLT+1), X(5*NXTLT+1), X(6*NXTLT+1),
     *      X(7*NXTLT+1), Y(1), Y(NXTLT+1), Y(2*NXTLT+1), Y(3*NXTLT+1),
     *      Y(4*NXTLT+1), Y(5*NXTLT+1), Y(6*NXTLT+1), Y(7*NXTLT+1))
  50  CONTINUE
C
C IS THERE A FOUR FACTOR LEFT
C
  60  IF (N2POW-3*N8POW-1) 90, 70, 80
C
C GO THROUGH THE BASE 2 ITERATION
C
C
  70  CALL R2TX(NTHPO, X(1), X(2), Y(1), Y(2))
      GO TO 90
C
C GO THROUGH THE BASE 4 ITERATION
C
  80  CALL R4TX(NTHPO, X(1), X(2), X(3), X(4), Y(1), Y(2), Y(3), Y(4))
C
  90  DO 110 J=1,15
        L(J) = 1
        IF (J-N2POW) 100, 100, 110
 100    L(J) = 2**(N2POW+1-J)
 110  CONTINUE
      IJ = 1
      DO 130 J1=1,L1
      DO 130 J2=J1,L2,L1
      DO 130 J3=J2,L3,L2
      DO 130 J4=J3,L4,L3
      DO 130 J5=J4,L5,L4
      DO 130 J6=J5,L6,L5
      DO 130 J7=J6,L7,L6
      DO 130 J8=J7,L8,L7
      DO 130 J9=J8,L9,L8
      DO 130 J10=J9,L10,L9
      DO 130 J11=J10,L11,L10
      DO 130 J12=J11,L12,L11
      DO 130 J13=J12,L13,L12
      DO 130 J14=J13,L14,L13
      DO 130 JI=J14,L15,L14
        IF (IJ-JI) 120, 130, 130
 120    R = X(IJ)
        X(IJ) = X(JI)
        X(JI) = R
        FI = Y(IJ)
        Y(IJ) = Y(JI)
        Y(JI) = FI
 130    IJ = IJ + 1
      IF (IN.EQ.1) GO TO 150
      DO 140 I=1,NTHPO
        Y(I) = -Y(I)
 140  CONTINUE
      GO TO 170
 150  DO 160 I=1,NTHPO
        X(I) = X(I)/FN
        Y(I) = Y(I)/FN
 160  CONTINUE
 170  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R2TX
C RADIX 2 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE R2TX(NTHPO, CR0, CR1, CI0, CI1)
      DIMENSION CR0(2), CR1(2), CI0(2), CI1(2)
      DO 10 K=1,NTHPO,2
        R1 = CR0(K) + CR1(K)
        CR1(K) = CR0(K) - CR1(K)
        CR0(K) = R1
        FI1 = CI0(K) + CI1(K)
        CI1(K) = CI0(K) - CI1(K)
        CI0(K) = FI1
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R4TX
C RADIX 4 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE R4TX(NTHPO, CR0, CR1, CR2, CR3, CI0, CI1, CI2, CI3)
      DIMENSION CR0(2), CR1(2), CR2(2), CR3(2), CI0(2), CI1(2), CI2(2),
     *    CI3(2)
      DO 10 K=1,NTHPO,4
        R1 = CR0(K) + CR2(K)
        R2 = CR0(K) - CR2(K)
        R3 = CR1(K) + CR3(K)
        R4 = CR1(K) - CR3(K)
        FI1 = CI0(K) + CI2(K)
        FI2 = CI0(K) - CI2(K)
        FI3 = CI1(K) + CI3(K)
        FI4 = CI1(K) - CI3(K)
        CR0(K) = R1 + R3
        CI0(K) = FI1 + FI3
        CR1(K) = R1 - R3
        CI1(K) = FI1 - FI3
        CR2(K) = R2 - FI4
        CI2(K) = FI2 + R4
        CR3(K) = R2 + FI4
        CI3(K) = FI2 - R4
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  R8TX
C RADIX 8 ITERATION SUBROUTINE
C-----------------------------------------------------------------------
C
      SUBROUTINE R8TX(NXTLT, NTHPO, LENGT, CR0, CR1, CR2, CR3, CR4,
     *    CR5, CR6, CR7, CI0, CI1, CI2, CI3, CI4, CI5, CI6, CI7)
      DIMENSION CR0(2), CR1(2), CR2(2), CR3(2), CR4(2), CR5(2), CR6(2),
     *    CR7(2), CI1(2), CI2(2), CI3(2), CI4(2), CI5(2), CI6(2),
     *    CI7(2), CI0(2)
      COMMON /CON2/ PI2, P7
C
      SCALE = PI2/FLOAT(LENGT)
      DO 30 J=1,NXTLT
        ARG = FLOAT(J-1)*SCALE
        C1 = COS(ARG)
        S1 = SIN(ARG)
        C2 = C1**2 - S1**2
        S2 = C1*S1 + C1*S1
        C3 = C1*C2 - S1*S2
        S3 = C2*S1 + S2*C1
        C4 = C2**2 - S2**2
        S4 = C2*S2 + C2*S2
        C5 = C2*C3 - S2*S3
        S5 = C3*S2 + S3*C2
        C6 = C3**2 - S3**2
        S6 = C3*S3 + C3*S3
        C7 = C3*C4 - S3*S4
        S7 = C4*S3 + S4*C3
        DO 20 K=J,NTHPO,LENGT
          AR0 = CR0(K) + CR4(K)
          AR1 = CR1(K) + CR5(K)
          AR2 = CR2(K) + CR6(K)
          AR3 = CR3(K) + CR7(K)
          AR4 = CR0(K) - CR4(K)
          AR5 = CR1(K) - CR5(K)
          AR6 = CR2(K) - CR6(K)
          AR7 = CR3(K) - CR7(K)
          AI0 = CI0(K) + CI4(K)
          AI1 = CI1(K) + CI5(K)
          AI2 = CI2(K) + CI6(K)
          AI3 = CI3(K) + CI7(K)
          AI4 = CI0(K) - CI4(K)
          AI5 = CI1(K) - CI5(K)
          AI6 = CI2(K) - CI6(K)
          AI7 = CI3(K) - CI7(K)
          BR0 = AR0 + AR2
          BR1 = AR1 + AR3
          BR2 = AR0 - AR2
          BR3 = AR1 - AR3
          BR4 = AR4 - AI6
          BR5 = AR5 - AI7
          BR6 = AR4 + AI6
          BR7 = AR5 + AI7
          BI0 = AI0 + AI2
          BI1 = AI1 + AI3
          BI2 = AI0 - AI2
          BI3 = AI1 - AI3
          BI4 = AI4 + AR6
          BI5 = AI5 + AR7
          BI6 = AI4 - AR6
          BI7 = AI5 - AR7
          CR0(K) = BR0 + BR1
          CI0(K) = BI0 + BI1
          IF (J.LE.1) GO TO 10
          CR1(K) = C4*(BR0-BR1) - S4*(BI0-BI1)
          CI1(K) = C4*(BI0-BI1) + S4*(BR0-BR1)
          CR2(K) = C2*(BR2-BI3) - S2*(BI2+BR3)
          CI2(K) = C2*(BI2+BR3) + S2*(BR2-BI3)
          CR3(K) = C6*(BR2+BI3) - S6*(BI2-BR3)
          CI3(K) = C6*(BI2-BR3) + S6*(BR2+BI3)
          TR = P7*(BR5-BI5)
          TI = P7*(BR5+BI5)
          CR4(K) = C1*(BR4+TR) - S1*(BI4+TI)
          CI4(K) = C1*(BI4+TI) + S1*(BR4+TR)
          CR5(K) = C5*(BR4-TR) - S5*(BI4-TI)
          CI5(K) = C5*(BI4-TI) + S5*(BR4-TR)
          TR = -P7*(BR7+BI7)
          TI = P7*(BR7-BI7)
          CR6(K) = C3*(BR6+TR) - S3*(BI6+TI)
          CI6(K) = C3*(BI6+TI) + S3*(BR6+TR)
          CR7(K) = C7*(BR6-TR) - S7*(BI6-TI)
          CI7(K) = C7*(BI6-TI) + S7*(BR6-TR)
          GO TO 20
  10      CR1(K) = BR0 - BR1
          CI1(K) = BI0 - BI1
          CR2(K) = BR2 - BI3
          CI2(K) = BI2 + BR3
          CR3(K) = BR2 + BI3
          CI3(K) = BI2 - BR3
          TR = P7*(BR5-BI5)
          TI = P7*(BR5+BI5)
          CR4(K) = BR4 + TR
          CI4(K) = BI4 + TI
          CR5(K) = BR4 - TR
          CI5(K) = BI4 - TI
          TR = -P7*(BR7+BI7)
          TI = P7*(BR7-BI7)
          CR6(K) = BR6 + TR
          CI6(K) = BI6 + TI
          CR7(K) = BR6 - TR
          CI7(K) = BI6 - TI
  20    CONTINUE
  30  CONTINUE
      RETURN
      END
Return to Main Software Page