Basic DIF Radix-4 FFT Algorithm Below is the Fortran code for a Decimation-in-Frequency, Radix-4, three butterfly Cooley-Tukey FFT followed by a bit-reversing unscrambler. Twiddle factors are precalculated and stored in arrays WR and WI.
C
C   A COOLEY-TUKEY RADIX-4 DIF  FFT PROGRAM
C   THREE BF, MULTIPLICATIONS BY  1, J, ETC. ARE REMOVED
C   COMPLEX INPUT DATA IN ARRAYS X AND Y
C   LENGTH IS  N = 4 ** M
C   TABLE LOOKUP OF W VALUES
C
C     C. S. BURRUS, RICE UNIVERSITY,  SEPT 1983
C
C---------------------------------------------------------
C
    SUBROUTINE  FFT4 (X,Y,N,M,WR,WI)
    REAL X(1), Y(1), WR(1), WI(1)
    DATA C21 / 0.707106778 /
C
C--------------MAIN FFT LOOPS-----------------------------
C
    N2 = N
    DO 10 K = 1, M
        N1 = N2
        N2 = N2/4
        JT = N2/2 + 1
C---------------SPECIAL BUTTERFLY FOR W = 1---------------
        DO 1 I = 1, 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)
C
            X(I) = R1 + R2
            X(I2)= R1 - R2
            X(I3)= R3 - S4
            X(I1)= R3 + S4
C
            Y(I) = S1 + S2
            Y(I2)= S1 - S2
            Y(I3)= S3 + R4
            Y(I1)= S3 - R4
C
   1        CONTINUE
        IF (K.EQ.M) GOTO 10
        IE = N/N1
        IA1 = 1
 
C--------------GENERAL BUTTERFLY-----------------
        DO 20 J = 2, N2
        IA1  = IA1 + IE
        IF (J.EQ.JT) GOTO 50
        IA2  = IA1 + IA1 - 1
            IA3  = IA2 + IA1 - 1
            CO1  = WR(IA1)
            CO2  = WR(IA2)
            CO3  = WR(IA3)
            SI1  = WI(IA1)
            SI2  = WI(IA2)
            SI3  = WI(IA3)
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)
C
            X(I) = R1 + R2
            R2   = R1 - R2
            R1   = R3 - S4
            R3   = R3 + S4
C
            Y(I) = S1 + S2
            S2   = S1 - S2
            S1   = S3 + R4
            S3   = S3 - R4
C
            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
        GOTO 20
C------------------SPECIAL BUTTERFLY FOR  W = J-----------
  50            DO 40 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)
C
            X(I) = R1 + R2
            Y(I2)=-R1 + R2
            R1   = R3 - S4
            R3   = R3 + S4
C
            Y(I) = S1 + S2
            X(I2)= S1 - S2
            S1   = S3 + R4
            S3   = S3 - R4
C
            X(I1) = (S3 + R3)*C21
            Y(I1) = (S3 - R3)*C21
            X(I3) = (S1 - R1)*C21
            Y(I3) =-(S1 + R1)*C21
  40            CONTINUE
  20        CONTINUE
  10    CONTINUE
C-----------DIGIT REVERSE COUNTER----------
  100   J = 1
    N1 = N - 1
    DO 104 I = 1, N1
        IF (I.GE.J) GOTO 101
        R1   = X(J)
        X(J) = X(I)
        X(I) = R1
        R1   = Y(J)
        Y(J) = Y(I)
        Y(I) = R1
 101        K = N/4
 102        IF (K*3.GE.J) GOTO 103
                J = J - K*3
                K = K/4
                GOTO 102
 103        J = J + K
 104    CONTINUE
    RETURN
    END