Basic DIF Radix-4 FFT Algorithm

Below is the Fortran code for a simple Decimation-in-Frequency, Radix-4, one butterfly Cooley-Tukey FFT to be followed by an unscrambler.

C   A COOLEY-TUKEY RADIX-4 DIF  FFT PROGRAM
C   COMPLEX INPUT DATA IN ARRAYS X AND Y
C   LENGTH IS  N = 4 ** M
C     C. S. BURRUS, RICE UNIVERSITY, SEPT 1983
C---------------------------------------------------------
    SUBROUTINE  FFT4 (X,Y,N,M)
    REAL X(1), Y(1)
C--------------MAIN FFT LOOPS-----------------------------
    N2 = N
    DO 10 K = 1, M
        N1 = N2
        N2 = N2/4
        E = 6.283185307179586/N1
        A = 0
C--------------------MAIN BUTTERFLIES-------------------
        DO 20 J=1, N2
            B    = A + A
            C    = A + B
            CO1  = COS(A)
            CO2  = COS(B)
            CO3  = COS(C)
            SI1  = SIN(A)
            SI2  = SIN(B)
            SI3  = SIN(C)
            A    = J*E
C----------------BUTTERFLIES WITH SAME W---------------
            DO 30 I=J, N, N1
            I1 = I  + N2
            I2 = I1 + N2
            I3 = I2 + N2
            R1 = X(I ) + X(I2)
            R3 = X(I ) - X(I2)
            S1 = Y(I ) + Y(I2)
            S3 = Y(I ) - Y(I2)
            R2 = X(I1) + X(I3)
            R4 = X(I1) - X(I3)
            S2 = Y(I1) + Y(I3)
            S4 = Y(I1) - Y(I3)
            X(I) = R1 + R2
            R2   = R1 - R2
            R1   = R3 - S4
            R3   = R3 + S4
            Y(I) = S1 + S2
            S2   = S1 - S2
            S1   = S3 + R4
            S3   = S3 - R4
            X(I1) = CO1*R3 + SI1*S3
            Y(I1) = CO1*S3 - SI1*R3
            X(I2) = CO2*R2 + SI2*S2
            Y(I2) = CO2*S2 - SI2*R2
            X(I3) = CO3*R1 + SI3*S1
            Y(I3) = CO3*S1 - SI3*R1
  30            CONTINUE
  20        CONTINUE
  10    CONTINUE
C-----------DIGIT REVERSE COUNTER goes here-----
    RETURN
    END