ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  A real-valued, in place, split-radix IFFT program
c  Hermitian symmetric input and real output in array X
c  Length is N = 2 ** M
c  Decimation-in-frequency, cos/sin in second loop
c  Input order:
c         [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
c
c  This FFT computes
c     x(j) = (1/N) * sum_{k=0}^{N-1} X(k)*exp(2ijk*pi/N)
c
c
c  H.V. Sorensen, Rice University, Nov. 1985
c
c  Reference:  H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus;
c              Real-Valued Fast Fourier Transform Algorithms; IEEE
c              Trans. Acoust., Speech, Signal Process.; Vol ASSP-35,
c              June 1987, pp. 849-863.
c
c  This code was originally named IRVFFT.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE SFFTCB( X, N, M )
c ... Scalar arguments ...
      INTEGER N, M
c ... Array arguments ...
      REAL  X(*)
c ... Local scalars ...
      INTEGER  J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8
      INTEGER  N1, N2, N4, N8
      REAL  SQRT2, TWOPI, XT, R1, T1, T2, T3, T4, T5
      REAL  A, A3, E, CC1, SS1, CC3, SS3
c ... Intrinsic functions ...
      INTRINSIC  SIN, COS
c ... Parameters ...
      PARAMETER ( SQRT2 = 1.4142135623730950488,
     +            TWOPI = 6.2831853071795864769 )
c
c ... Exe. statements ...
c
      IF ( N .EQ. 1 ) RETURN
c
      N2 = 2 * N
      DO 10, K = 1, M-1
         IS = 0
         ID = N2
         N2 = N2 / 2
         N4 = N2 / 4
         N8 = N4 / 2
         E = TWOPI / N2
 17      DO 15, I = IS, N-1, ID
            I1 = I + 1
            I2 = I1 + N4
            I3 = I2 + N4
            I4 = I3 + N4
            T1 = X(I1) - X(I3)
            X(I1) = X(I1) + X(I3)
            X(I2) = 2 * X(I2)
            X(I3) = T1 - 2 * X(I4)
            X(I4) = T1 + 2 * X(I4)
            IF ( N4 .EQ. 1 ) GOTO 15
            I1 = I1 + N8
            I2 = I2 + N8
            I3 = I3 + N8
            I4 = I4 + N8
            T1 = ( X(I2) - X(I1) ) / SQRT2
            T2 = ( X(I4) + X(I3) ) / SQRT2
            X(I1) = X(I1) + X(I2)
            X(I2) = X(I4) - X(I3)
            X(I3) = 2 * ( - T2 - T1 )
            X(I4) = 2 * ( -T2 + T1 )
 15      CONTINUE
         IS = 2 * ID - N2
         ID = 4 * ID
         IF ( IS .LT. N-1 ) GOTO 17
         A = E
         DO 20, J = 2, N8
            A3 = 3 * A
            CC1 = COS(A)
            SS1 = SIN(A)
            CC3 = COS(A3)
            SS3 = SIN(A3)
            A = J * E
            IS = 0
            ID = 2 * N2
 40         DO 30, I = IS, N-1, ID
               I1 = I + J
               I2 = I1 + N4
               I3 = I2 + N4
               I4 = I3 + N4
               I5 = I + N4 - J + 2
               I6 = I5 + N4
               I7 = I6 + N4
               I8 = I7 + N4
               T1 = X(I1) - X(I6)
               X(I1) = X(I1) + X(I6)
               T2 = X(I5) - X(I2)
               X(I5) = X(I2) + X(I5)
               T3 = X(I8) + X(I3)
               X(I6) = X(I8) - X(I3)
               T4 = X(I4) + X(I7)
               X(I2) = X(I4) - X(I7)
               T5 = T1 - T4
               T1 = T1 + T4
               T4 = T2 - T3
               T2 = T2 + T3
               X(I3) = T5 * CC1 + T4 * SS1
               X(I7) = - T4 * CC1 + T5 * SS1
               X(I4) = T1 * CC3 - T2 * SS3
               X(I8) = T2 * CC3 + T1 * SS3
 30         CONTINUE
            IS = 2 * ID - N2
            ID = 4 * ID
            IF ( IS .LT. N-1 ) GOTO 40
 20      CONTINUE
 10   CONTINUE
c
      IS = 1
      ID = 4
 70   DO 60, I0 = IS, N, ID
         I1 = I0 + 1
         R1 = X(I0)
         X(I0) = R1 + X(I1)
         X(I1) = R1 - X(I1)
 60   CONTINUE
      IS = 2 * ID - 1
      ID = 4 * ID
      IF ( IS .LT. N ) GOTO 70
c
 100  J = 1
      N1 = N - 1
      DO 104, I = 1, N1
         IF ( I .GE. J ) GOTO 101
         XT = X(J)
         X(J) = X(I)
         X(I) = XT
 101     K = N / 2
 102     IF ( K .GE. J ) GOTO 103
            J = J - K
            K = K / 2
            GOTO 102
 103     J = J + K
 104  CONTINUE
      XT = 1.0 / FLOAT( N )
      DO 99, I = 1, N
         X(I) = XT * X(I)
 99   CONTINUE
      RETURN
c
c ... End of subroutine SFFTCB ...
c
      END