C---------------------------------------------------------------C
C                                                               C
C   Real-valued, in-place, Cooley-Tukey radix-2 FFT             C
C   Real input and output data in array X                       C
C   Length is N = 2 ** M                                        C
C   Decimation-in-time, cos/sin in innermost loop               C
C   Output in order:                                            C
C       [Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1)]     C
C                                                               C
C   H.V. Sorensen,  Rice University,  January 1985              C
C                                                               C
C   For reference see:                                          C
C   H.V. Sorensen, D.L. Jones, M.T. Heideman, and C.S. Burrus,  C
C   Real-Valued Fast Fourier Transform Algorithms, IEEE Trans.  C
C   Acoust., Speech, Signal Processing, ASSP-35 (1987), pp.     C
C   849-862                                                     C
C                                                               C
C---------------------------------------------------------------C
C
      SUBROUTINE RFFT( X, N, M )
      INTEGER  N, M, J, N1, I, K, N2, N4, I1, I2, I3, I4
      REAL  X(*)
      REAL  XT, A, E, CC, SS, T1, T2, TWOPI
      PARAMETER  ( TWOPI = 6.2831853071795864769 )
C
C-----Digit reverse counter-------------------------------------C
      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
C
C-----Length two butterflies------------------------------------C
      DO 60, I = 1, N, 2
         XT = X(I)
         X(I) = XT + X(I+1)
         X(I+1) = XT - X(I+1)
 60   CONTINUE
C
C-----Other butterflies-----------------------------------------C
      N2 = 1
      DO 10, K = 2, M
         N4 = N2
         N2 = 2 * N4
         N1 = 2 * N2
         E = TWOPI / N1
         DO 20, I = 1, N, N1
            XT = X(I)
            X(I) = XT + X(I+N2)
            X(I+N2) = XT - X(I+N2)
            X(I+N4+N2) = - X(I+N4+N2)
            A = E
C--Note that in the first run through, N4=1, so the next statement
C--is 'DO 30, J = 1, 0', which might cause problems on some compilers
C--         IF ( N4 .EQ. 1 ) GOTO 20
            DO 30, J = 1, N4-1
               I1 = I + J
               I2 = I - J + N2
               I3 = I + J + N2
               I4 = I - J + N1
               CC = COS( A )
               SS = SIN( A )
               A = A + E
               T1 = X(I3) * CC + X(I4) * SS
               T2 = X(I3) * SS - X(I4) * CC
               X(I4) = X(I2) - T2
               X(I3) = - X(I2) - T2
               X(I2) = X(I1) - T1
               X(I1) = X(I1) + T1
 30         CONTINUE
 20      CONTINUE
 10   CONTINUE
C
      RETURN
      END