DIF Split Radix FFT Algorithm Below is the Fortran code for a simple Decimation-in-Frequency, Split-Radix, two butterfly FFT to be followed by a bit-reversing unscrambler. Twiddle factors are precalculated and stored in arrays WR and WI.
C--------------------------------------------------------------C
C       A DUHAMEL-HOLLMAN SPLIT RADIX FFT                      C
C       REF: ELECTRONICS LETTERS, JAN. 5, 1984                 C
C       COMPLEX INPUT AND OUTPUT DATA IN ARRAYS X AND Y        C
C       LENGTH IS N = 2 ** M,  OUTPUT IN BIT-REVERSED ORDER    C
C   TWO BUTTERFLIES TO REMOVE MULTS BY UNITY               C
C       SPECIAL LAST TWO STAGES                                C
C   TABLE LOOK-UP OF SINE AND COSINE VALUES            C
C       C.S. BURRUS,       RICE UNIV.       APRIL 1985         C
C--------------------------------------------------------------C
C
        SUBROUTINE FFT(X,Y,N,M,WR,WI)
        REAL X(1),Y(1),WR(1),WI(1)
    C81= 0.707106778
        N2 = 2*N
        DO  10 K = 1, M-3
        IS  = 1
        ID  = N2
        N2 = N2/2
        N4 = N2/4
 2      DO 1 I0 = IS, N-1, ID
        I1 = I0 + N4
        I2 = I1 + N4
            I3 = I2 + N4
        R1    = X(I0) - X(I2)
        X(I0) = X(I0) + X(I2)
        R2    = Y(I1) - Y(I3)
        Y(I1) = Y(I1) + Y(I3)
        X(I2) = R1 + R2
        R2    = R1 - R2
        R1    = X(I1) - X(I3)
        X(I1) = X(I1) + X(I3)
            X(I3) = R2
            R2    = Y(I0) - Y(I2)
                Y(I0) = Y(I0) + Y(I2)
            Y(I2) =-R1 + R2
            Y(I3) = R1 + R2
  1     CONTINUE
        IS = 2*ID - N2 + 1
        ID = 4*ID
            IF (IS.LT.N) GOTO 2
        IE  = N/N2
            IA1 = 1
            DO  20 J = 2, N4
                IA1 = IA1 + IE
                IA3 = 3*IA1 - 2
                CC1 = WR(IA1)
                SS1 = WI(IA1)
                CC3 = WR(IA3)
                SS3 = WI(IA3)
                IS  = J
                ID  = 2*N2
 40             DO 30 I0 = IS, N-1, ID
                    I1 = I0 + N4
                    I2 = I1 + N4
                    I3 = I2 + N4
C
                    R1    = X(I0) - X(I2)
                    X(I0) = X(I0) + X(I2)
                    R2    = X(I1) - X(I3)
                    X(I1) = X(I1) + X(I3)
                    S1    = Y(I0) - Y(I2)
                    Y(I0) = Y(I0) + Y(I2)
                    S2    = Y(I1) - Y(I3)
                    Y(I1) = Y(I1) + Y(I3)
C
                    S3    = R1 - S2
                    R1    = R1 + S2
                    S2    = R2 - S1
                    R2    = R2 + S1
                    X(I2) = R1*CC1 - S2*SS1
                    Y(I2) =-S2*CC1 - R1*SS1
                    X(I3) = S3*CC3 + R2*SS3
                    Y(I3) = R2*CC3 - S3*SS3
 30             CONTINUE
                IS = 2*ID - N2 + J
                ID = 4*ID
                IF (IS.LT.N) GOTO 40
  20        CONTINUE
  10    CONTINUE
C
        IS =  1
        ID = 32
  50    DO 60 I = IS, N, ID
            I0    = I + 8
            DO 15 J = 1, 2
               R1 = X(I0)   + X(I0+2)
               R3 = X(I0)   - X(I0+2)
               R2 = X(I0+1) + X(I0+3)
               R4 = X(I0+1) - X(I0+3)
               X(I0)   = R1 + R2
               X(I0+1) = R1 - R2
               R1 = Y(I0)   + Y(I0+2)
               S3 = Y(I0)   - Y(I0+2)
               R2 = Y(I0+1) + Y(I0+3)
               S4 = Y(I0+1) - Y(I0+3)
               Y(I0)   = R1 + R2
               Y(I0+1) = R1 - R2
               Y(I0+2) = S3 - R4
               Y(I0+3) = S3 + R4
               X(I0+2) = R3 + S4
               X(I0+3) = R3 - S4
               I0 = I0 + 4
  15        CONTINUE
  60    CONTINUE
            IS = 2*ID - 15
            ID = 4*ID
        IF (IS.LT.N) GOTO 50
C
        IS =  1
        ID = 16
  55    DO 65 I0 = IS, N, ID
                R1 = X(I0)   + X(I0+4)
            R5 = X(I0)   - X(I0+4)
            R2 = X(I0+1) + X(I0+5)
            R6 = X(I0+1) - X(I0+5)
            R3 = X(I0+2) + X(I0+6)
            R7 = X(I0+2) - X(I0+6)
            R4 = X(I0+3) + X(I0+7)
            R8 = X(I0+3) - X(I0+7)
            T1 = R1 - R3
            R1 = R1 + R3
            R3 = R2 - R4
            R2 = R2 + R4
            X(I0)   = R1 + R2
            X(I0+1) = R1 - R2
C
            R1 = Y(I0)   + Y(I0+4)
            S5 = Y(I0)   - Y(I0+4)
            R2 = Y(I0+1) + Y(I0+5)
            S6 = Y(I0+1) - Y(I0+5)
            S3 = Y(I0+2) + Y(I0+6)
            S7 = Y(I0+2) - Y(I0+6)
            R4 = Y(I0+3) + Y(I0+7)
            S8 = Y(I0+3) - Y(I0+7)
            T2 = R1 - S3
            R1 = R1 + S3
            S3 = R2 - R4
            R2 = R2 + R4
            Y(I0)   = R1 + R2
            Y(I0+1) = R1 - R2
            X(I0+2) = T1 + S3
            X(I0+3) = T1 - S3
            Y(I0+2) = T2 - R3
            Y(I0+3) = T2 + R3
C
            R1 = (R6 - R8)*C81
            R6 = (R6 + R8)*C81
            R2 = (S6 - S8)*C81
            S6 = (S6 + S8)*C81
C
            T1 = R5 - R1
            R5 = R5 + R1
            R8 = R7 - R6
            R7 = R7 + R6
            T2 = S5 - R2
            S5 = S5 + R2
            S8 = S7 - S6
            S7 = S7 + S6
            X(I0+4) = R5 + S7
            X(I0+7) = R5 - S7
            X(I0+5) = T1 + S8
            X(I0+6) = T1 - S8
            Y(I0+4) = S5 - R7
            Y(I0+7) = S5 + R7
            Y(I0+5) = T2 - R8
            Y(I0+6) = T2 + R8
  65    CONTINUE
            IS = 2*ID - 7
            ID = 4*ID
        IF (IS.LT.N) GOTO 55
C
C------------BIT REVERSE COUNTER-----------------
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
            XT   = Y(J)
            Y(J) = Y(I)
            Y(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
        RETURN
        END