Basic DIF Split Radix FFT Algorithm

Below is the Fortran code for a simple Decimation-in-Frequency, Split-Radix, one butterfly FFT to be followed by a bit-reversing unscrambler.

C   A DUHAMEL-HOLLMANN SPLIT RADIX  FFT PROGRAM
C   FROM: ELECTRONICS LETTERS, JAN. 5, 1984
C   COMPLEX INPUT DATA IN ARRAYS X AND Y
C   LENGTH IS  N = 2 ** M
C     C. S. BURRUS, RICE UNIVERSITY, MARCH 1984
C
C---------------------------------------------------------
    SUBROUTINE  FFT (X,Y,N,M)
    REAL X(1), Y(1)
C--------------MAIN FFT LOOPS-----------------------------
C
    N1 = N
    N2 = N/2
    IP = 0
    IS = 1
    A  = 6.283185307179586/N
    DO 10 K = 1, M-1
        JD = N1 + N2
        N1 = N2
        N2 = N2/2
        J0 = N1*IP + 1
        IP = 1 - IP
        DO 20 J = J0, N, JD
            JS = 0
        JT = J + N2 - 1
            DO 30 I = J, JT
            JSS= JS*IS
            JS = JS + 1
                C1 = COS(A*JSS)
                C3 = COS(3*A*JSS)
                S1 = -SIN(A*JSS)
                S3 = -SIN(3*A*JSS)
            I1 = I  + N2
            I2 = I1 + N2
            I3 = I2 + N2
            R1    = X(I ) + X(I2)
            R2    = X(I ) - X(I2)
            R3    = X(I1) - X(I3)
            X(I2) = X(I1) + X(I3)
            X(I1) = R1
C
            R1    = Y(I ) + Y(I2)
            R4    = Y(I ) - Y(I2)
            R5    = Y(I1) - Y(I3)
            Y(I2) = Y(I1) + Y(I3)
            Y(I1) = R1
C
            R1    = R2 - R5
            R2    = R2 + R5
            R5    = R4 + R3
            R4    = R4 - R3
C
            X(I)  = C1*R1 + S1*R5
            Y(I)  = C1*R5 - S1*R1
            X(I3) = C3*R2 + S3*R4
            Y(I3) = C3*R4 - S3*R2
  30            CONTINUE
  20        CONTINUE
        IS = IS + IS
  10    CONTINUE
    IP = 1 - IP
    J0 = 2 - IP
    DO 5 I = J0, N-1, 3
       I1 = I + 1
       R1    = X(I) + X(I1)
       X(I1) = X(I) - X(I1)
       X(I)  = R1
       R1    = Y(I) + Y(I1)
       Y(I1) = Y(I) - Y(I1)
       Y(I)  = R1
   5    CONTINUE
    RETURN
    END