Skip to main content
Engineering LibreTexts

14.3: Appendix 3 - FFT Computer Programs

  • Page ID
    2051
  • [ "article:topic" ]

    Goertzel Algorithm

    A FORTRAN implementation of the first-order Goertzel algorithm with in-order input is given below.

    C----------------------------------------------
    C   GOERTZEL'S  DFT  ALGORITHM
    C   First order, input inorder
    C   C. S. BURRUS,   SEPT 1983
    C---------------------------------------------
        SUBROUTINE DFT(X,Y,A,B,N)
        REAL X(260), Y(260), A(260), B(260)
        Q = 6.283185307179586/N
        DO 20 J=1, N
           C  = COS(Q*(J-1))
           S  = SIN(Q*(J-1))
           AT = X(1)
           BT = Y(1)
           DO 30 I = 2, N
              T  = C*AT - S*BT + X(I)
              BT = C*BT + S*AT + Y(I)
              AT = T
    30     CONTINUE
           A(J) = C*AT - S*BT
           B(J) = C*BT + S*AT
    20  CONTINUE
        RETURN
        END
    NOT_CONVERTED_YET: caption

    First Order Goertzel Algorithm

    Second Order Goertzel Algorithm

    Below is the program for a second order Goertzel algorithm.

    C----------------------------------------------
    C   GOERTZEL'S  DFT  ALGORITHM
    C   Second order, input inorder
    C   C. S. BURRUS,   SEPT 1983
    C---------------------------------------------
        SUBROUTINE DFT(X,Y,A,B,N)
        REAL X(260), Y(260), A(260), B(260)
    C
        Q = 6.283185307179586/N
        DO 20 J = 1, N
           C  = COS(Q*(J-1))
           S  = SIN(Q*(J-1))
           CC = 2*C
           A2 = 0
           B2 = 0
           A1 = X(1)
           B1 = Y(1)
           DO 30 I = 2, N
              T  = A1
              A1 = CC*A1 - A2 + X(I)
              A2 = T
              T  = B1
              B1 = CC*B1 - B2 + Y(I)
              B2 = T
    30     CONTINUE
           A(J) = C*A1 - A2 - S*B1
           B(J) = C*B1 - B2 + S*A1
    20  CONTINUE
    C
        RETURN
        END
    NOT_CONVERTED_YET: caption

    Second Order Goertzel Algorithm

    Second Order Goertzel Algorithm 2

    Second order Goertzel algorithm that calculates two outputs at a time.

    C-------------------------------------------------------
    C GOERTZEL'S  DFT  ALGORITHM,  Second order
    C Input inorder, output by twos;  C.S. Burrus, SEPT 1991
    C-------------------------------------------------------
        SUBROUTINE DFT(X,Y,A,B,N)
        REAL X(260), Y(260), A(260), B(260)
        Q = 6.283185307179586/N
        DO 20 J = 1, N/2 + 1
           C  = COS(Q*(J-1))
           S  = SIN(Q*(J-1))
           CC = 2*C
           A2 = 0
           B2 = 0
           A1 = X(1)
           B1 = Y(1)
           DO 30 I = 2, N
              T  = A1
              A1 = CC*A1 - A2 + X(I)
              A2 = T
              T  = B1
              B1 = CC*B1 - B2 + Y(I)
              B2 = T
    30     CONTINUE
           A2  = C*A1 - A2
           T   = S*B1
           A(J)     = A2 - T
           A(N-J+2) = A2 + T
           B2  = C*B1 - B2
           T   = S*A1
           B(J)     = B2 + T
           B(N-J+2) = B2 - T
    20  CONTINUE
        RETURN
        END

    Figure.  Second Order Goertzel Calculating Two Outputs at a Time

    Basic QFT Algorithm

    A FORTRAN implementation of the basic QFT algorithm is given below to show how the theory is implemented. The program is written for clarity, not to minimize the number of floating point operations.

    C
       SUBROUTINE QDFT(X,Y,XX,YY,NN)
       REAL X(0:260),Y(0:260),XX(0:260),YY(0:260)
       C
       N1 = NN - 1
       N2 = N1/2
       N21 = NN/2
       Q   = 6.283185308/NN
       DO 2 K = 0, N21
          SSX = X(0)
          SSY = Y(0)
          SDX = 0
          SDY = 0
          IF (MOD(NN,2).EQ.0) THEN
             SSX = SSX + COS(3.1426*K)*X(N21)
             SSY = SSY + COS(3.1426*K)*Y(N21)
          ENDIF
          DO 3 N = 1, N2
             SSX = SSX + (X(N) + X(NN-N))*COS(Q*N*K)
             SSY = SSY + (Y(N) + Y(NN-N))*COS(Q*N*K)
             SDX = SDX + (X(N) - X(NN-N))*SIN(Q*N*K)
             SDY = SDY + (Y(N) - Y(NN-N))*SIN(Q*N*K)
    3     CONTINUE
          XX(K) = SSX + SDY
          YY(K) = SSY - SDX
          XX(NN-K) = SSX - SDY
          YY(NN-K) = SSY + SDX
    2  CONTINUE
       RETURN
       END

    NOT_CONVERTED_YET: caption

    Simple QFT Fortran Program

    Basic Radix-2 FFT Algorithm

    Below is the Fortran code for a simple Decimation-in-Frequency, Radix-2, one butterfly Cooley-Tukey FFT followed by a bit-reversing unscrambler.

    C
    C   A COOLEY-TUKEY RADIX-2, DIF  FFT PROGRAM
    C   COMPLEX INPUT DATA IN ARRAYS X AND Y
    C      C. S. BURRUS, RICE UNIVERSITY, SEPT 1983
    C---------------------------------------------------------
        SUBROUTINE FFT (X,Y,N,M)
        REAL X(1), Y(1)
    C--------------MAIN FFT LOOPS-----------------------------
    C
        N2 = N
        DO 10 K = 1, M
            N1 = N2
            N2 = N2/2
            E  = 6.283185307179586/N1
            A  = 0
            DO 20 J = 1, N2
            C = COS (A)
            S = SIN (A)
            A = J*E
            DO 30 I = J, N, N1
                        L = I + N2
                        XT   = X(I) - X(L)
                        X(I) = X(I) + X(L)
                        YT   = Y(I) - Y(L)
                        Y(I) = Y(I) + Y(L)
                        X(L) = C*XT + S*YT
                        Y(L) = C*YT - S*XT
       30           CONTINUE
       20       CONTINUE
       10   CONTINUE
    C
    C------------DIGIT REVERSE COUNTER-----------------
      100   J = 1
        N1 = N - 1
        DO 104 I=1, N1
            IF (I.GE.J) GOXTO 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

    Figure: Radix-2, DIF, One Butterfly Cooley-Tukey FFT

    Basic DIT Radix-2 FFT Algorithm

    Below is the Fortran code for a simple Decimation-in-Time, Radix-2, one butterfly Cooley-Tukey FFT preceeded by a bit-reversing scrambler.

    C
    C   A COOLEY-TUKEY RADIX-2, DIT  FFT PROGRAM
    C   COMPLEX INPUT DATA IN ARRAYS X AND Y
    C      C. S. BURRUS, RICE UNIVERSITY, SEPT 1985
    C
    C---------------------------------------------------------
        SUBROUTINE FFT (X,Y,N,M)
        REAL X(1), Y(1)
    C------------DIGIT 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
    C--------------MAIN FFT LOOPS-----------------------------
    C
        N2 = 1
        DO 10 K = 1, M
            E  = 6.283185307179586/(2*N2)
            A  = 0
            DO 20 J = 1, N2
            C = COS (A)
            S = SIN (A)
            A = J*E
            DO 30 I = J, N, 2*N2
                        L = I + N2
                        XT = C*X(L) + S*Y(L)
                        YT = C*Y(L) - S*X(L)
                        X(L) = X(I) - XT
                        X(I) = X(I) + XT
                        Y(L) = Y(I) - YT
                        Y(I) = Y(I) + YT
       30           CONTINUE
       20       CONTINUE
            N2 = N2+N2
       10   CONTINUE
    C
        RETURN
        END
     

    DIF Radix-2 FFT Algorithm

    Below is the Fortran code for a Decimation-in-Frequency, Radix-2, three butterfly Cooley-Tukey FFT followed by a bit-reversing unscrambler.

    C   A COOLEY-TUKEY RADIX 2, DIF  FFT PROGRAM
    C   THREE-BF, MULT BY  1  AND  J  ARE REMOVED
    C   COMPLEX INPUT DATA IN ARRAYS X AND Y
    C   TABLE LOOK-UP OF W VALUES
    C      C. S. BURRUS, RICE UNIVERSITY, SEPT 1983
    C---------------------------------------------------------
        SUBROUTINE FFT (X,Y,N,M,WR,WI)
        REAL X(1), Y(1), WR(1), WI(1)
    C--------------MAIN FFT LOOPS-----------------------------
    C
        N2 = N
        DO 10 K = 1, M
            N1 = N2
            N2 = N2/2
            JT = N2/2 + 1
            DO 1 I = 1, N, N1
            L = I + N2
            T    = X(I) - X(L)
            X(I) = X(I) + X(L)
            X(L) = T
            T    = Y(I) - Y(L)
            Y(I) = Y(I) + Y(L)
            Y(L) = T
       1        CONTINUE
            IF (K.EQ.M) GOTO 10
            IE  = N/N1
            IA  = 1
            DO 20 J = 2, N2
            IA = IA + IE
            IF (J.EQ.JT) GOTO 50
            C = WR(IA)
            S = WI(IA)
            DO 30 I = J, N, N1
                L = I + N2
                T    = X(I) - X(L)
                X(I) = X(I) + X(L)
                TY   = Y(I) - Y(L)
                Y(I) = Y(I) + Y(L)
                X(L) = C*T + S*TY
                Y(L) = C*TY - S*T
       30       CONTINUE
            GOTO 25
       50       DO 40 I = J, N, N1
                L = I + N2
                T    = X(I) - X(L)
                X(I) = X(I) + X(L)
                TY   = Y(I) - Y(L)
                Y(I) = Y(I) + Y(L)
                X(L) = TY
                Y(L) =-T
       40       CONTINUE
       25       A = J*E
       20       CONTINUE
       10   CONTINUE
    C------------DIGIT REVERSE COUNTER Goes here----------
        RETURN
        END
     

    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

     

    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
     

    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
     

    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
     

    Prime Factor FFT Algorithm

    Below is the Fortran code for a Prime-Factor Algorithm (PFA) FFT allowing factors of the length of 2, 3, 4, 5, and 7. It is followed by an unscrambler.

    C---------------------------------------------------
    C
    C   A PRIME FACTOR FFT PROGRAM WITH GENERAL MODULES
    C   COMPLEX INPUT DATA IN ARRAYS  X AND Y
    C   COMPLEX OUTPUT IN  A AND B
    C   LENGTH  N  WITH  M  FACTORS IN ARRAY  NI
    C     N = NI(1)*NI(2)* ... *NI(M)
    C   UNSCRAMBLING CONSTANT  UNSC
    C     UNSC = N/NI(1) + N/NI(2) +...+ N/NI(M), MOD N
    C      C. S. BURRUS, RICE UNIVERSITY, JAN 1987
    C
    C--------------------------------------------------
    C
        SUBROUTINE PFA(X,Y,N,M,NI,A,B,UNSC)
    C
        INTEGER  NI(4), I(16), UNSC
            REAL X(1), Y(1), A(1), B(1)
    C
        DATA  C31, C32  / -0.86602540,-1.50000000 /
        DATA  C51, C52  /  0.95105652,-1.53884180 /
        DATA  C53, C54  / -0.36327126, 0.55901699 /
        DATA  C55       / -1.25  /
        DATA  C71, C72  / -1.16666667,-0.79015647 /
        DATA  C73, C74  /  0.055854267, 0.7343022 /
        DATA  C75, C76  /  0.44095855,-0.34087293 /
        DATA  C77, C78  /  0.53396936, 0.87484229 /
    C
    C-----------------NESTED LOOPS----------------------
    C
        DO 10 K=1, M
              N1 = NI(K)
          N2 = N/N1
          DO 15 J=1, N, N1
                 IT   = J
                 DO 30 L=1, N1
                    I(L) = IT
                A(L) = X(IT)
                B(L) = Y(IT)
                    IT = IT + N2
                    IF (IT.GT.N)  IT = IT - N
       30        CONTINUE
                 GOTO (20,102,103,104,105,20,107), N1
    C
    C----------------WFTA N=2--------------------------------
    C
      102   R1    = A(1)
        A(1)  = R1 + A(2)
        A(2)  = R1 - A(2)
    C
        R1   = B(1)
        B(1) = R1 + B(2)
        B(2) = R1 - B(2)
    C
        GOTO 20
    C----------------WFTA N=3--------------------------------
    C
      103   R2 = (A(2) - A(3)) * C31
        R1 =  A(2) + A(3)
        A(1)= A(1) + R1
        R1  = A(1) + R1 * C32
    C
        S2 = (B(2) - B(3)) * C31
        S1 =  B(2) + B(3)
        B(1)= B(1) + S1
        S1  = B(1) + S1 * C32
    C
        A(2) = R1 - S2
        A(3) = R1 + S2
        B(2) = S1 + R2
        B(3) = S1 - R2
    C
        GOTO 20
    C
    C----------------WFTA N=4---------------------------------
    C
      104   R1 = A(1) + A(3)
        T1 = A(1) - A(3)
        R2 = A(2) + A(4)
        A(1) = R1 + R2
        A(3) = R1 - R2
    C
        R1 = B(1) + B(3)
        T2 = B(1) - B(3)
        R2 = B(2) + B(4)
        B(1) = R1 + R2
        B(3) = R1 - R2
    C
        R1 = A(2) - A(4)
        R2 = B(2) - B(4)
    C
        A(2) = T1 + R2
        A(4) = T1 - R2
        B(2) = T2 - R1
        B(4) = T2 + R1
    C
        GOTO 20
    C
    C----------------WFTA N=5--------------------------------
    C
      105   R1 = A(2) + A(5)
        R4 = A(2) - A(5)
        R3 = A(3) + A(4)
        R2 = A(3) - A(4)
    C
        T = (R1 - R3) * C54
        R1 = R1 + R3
        A(1) = A(1) + R1
        R1   = A(1) + R1 * C55
    C
        R3 = R1 - T
        R1 = R1 + T
    C
        T = (R4 + R2) * C51
        R4 =  T + R4 * C52
        R2 =  T + R2 * C53
    C
        S1 = B(2) + B(5)
        S4 = B(2) - B(5)
        S3 = B(3) + B(4)
        S2 = B(3) - B(4)
    C
        T = (S1 - S3) * C54
        S1 = S1 + S3
        B(1) = B(1) + S1
        S1   = B(1) + S1 * C55
    C
        S3 =  S1 - T
        S1 =  S1 + T
    C
        T = (S4 + S2) * C51
        S4 =  T + S4 * C52
        S2 =  T + S2 * C53
    C
        A(2) = R1 + S2
        A(5) = R1 - S2
        A(3) = R3 - S4
        A(4) = R3 + S4
    C
        B(2) = S1 - R2
        B(5) = S1 + R2
        B(3) = S3 + R4
        B(4) = S3 - R4
    C
        GOTO 20
    C-----------------WFTA N=7--------------------------
    C
      107   R1 = A(2) + A(7)
        R6 = A(2) - A(7)
        S1 = B(2) + B(7)
        S6 = B(2) - B(7)
        R2 = A(3) + A(6)
        R5 = A(3) - A(6)
        S2 = B(3) + B(6)
        S5 = B(3) - B(6)
        R3 = A(4) + A(5)
        R4 = A(4) - A(5)
        S3 = B(4) + B(5)
        S4 = B(4) - B(5)
    C
        T3 = (R1 - R2) * C74
        T  = (R1 - R3) * C72
        R1 = R1 + R2 + R3
        A(1) = A(1) + R1
        R1   = A(1) + R1 * C71
        R2 =(R3 - R2) * C73
        R3 = R1 - T + R2
        R2 = R1 - R2 - T3
        R1 = R1 + T + T3
        T = (R6 - R5) * C78
        T3 =(R6 + R4) * C76
        R6 =(R6 + R5 - R4) * C75
        R5 =(R5 + R4) * C77
        R4 = R6 - T3 + R5
        R5 = R6 - R5 - T
        R6 = R6 + T3 + T
    C
        T3 = (S1 - S2) * C74
        T  = (S1 - S3) * C72
        S1 =  S1 + S2 + S3
        B(1) = B(1) + S1
        S1   = B(1) + S1 * C71
        S2 =(S3 - S2) * C73
        S3 = S1 - T  + S2
        S2 = S1 - S2 - T3
        S1 = S1 + T  + T3
        T  = (S6 - S5) * C78
        T3 = (S6 + S4) * C76
        S6 = (S6 + S5 - S4) * C75
        S5 = (S5 + S4) * C77
        S4 = S6 - T3 + S5
        S5 = S6 - S5 - T
        S6 = S6 + T3 + T
    C
        A(2) = R3 + S4
        A(7) = R3 - S4
        A(3) = R1 + S6
        A(6) = R1 - S6
        A(4) = R2 - S5
        A(5) = R2 + S5
        B(4) = S2 + R5
        B(5) = S2 - R5
        B(2) = S3 - R4
        B(7) = S3 + R4
        B(3) = S1 - R6
        B(6) = S1 + R6
    C
       20        IT   = J
                 DO 31 L=1, N1
                    I(L) = IT
                X(IT) = A(L)
                Y(IT) = B(L)
                    IT = IT + N2
                    IF (IT.GT.N)  IT = IT - N
       31        CONTINUE
       15     CONTINUE
       10   CONTINUE
    C
    C-----------------UNSCRAMBLING----------------------
    C
        L = 1
        DO 2 K=1, N
           A(K) = X(L)
           B(K) = Y(L)
               L = L + UNSC
               IF (L.GT.N)  L = L - N
        2   CONTINUE
        RETURN
        END
    C
     

    In Place, In Order Prime Factor FFT Algorithm

    Below is the Fortran code for a Prime-Factor Algorithm (PFA) FFT allowing factors of the length of 2, 3, 4, 5, 7, 8, 9, and 16. It is both in-place and in-order, so requires no unscrambler.

    C
    C   A PRIME FACTOR FFT PROGRAM
    C   IN-PLACE AND IN-ORDER
    C   COMPLEX INPUT DATA IN ARRAYS  X AND Y
    C   LENGTH  N  WITH  M  FACTORS IN ARRAY  NI
    C       N = NI(1)*NI(2)*...*NI(M)
    C   REDUCED TEMP STORAGE IN SHORT WFTA MODULES
    C   Has modules 2,3,4,5,7,8,9,16
    C   PROGRAM BY  C. S. BURRUS,  RICE UNIVERSITY
    C                  SEPT 1983
    C----------------------------------------------------
    C
        SUBROUTINE PFA(X,Y,N,M,NI)
        INTEGER  NI(4), I(16), IP(16), LP(16)
            REAL X(1), Y(1)
        DATA  C31, C32  / -0.86602540,-1.50000000 /
        DATA  C51, C52  /  0.95105652,-1.53884180 /
        DATA  C53, C54  / -0.36327126, 0.55901699 /
        DATA  C55       / -1.25  /
        DATA  C71, C72  / -1.16666667,-0.79015647 /
        DATA  C73, C74  /  0.055854267, 0.7343022 /
        DATA  C75, C76  /  0.44095855,-0.34087293 /
        DATA  C77, C78  /  0.53396936, 0.87484229 /
        DATA  C81       /  0.70710678 /
        DATA  C95       / -0.50000000 /
        DATA  C92, C93  /  0.93969262, -0.17364818 /
        DATA  C94, C96  /  0.76604444, -0.34202014 /
        DATA  C97, C98  / -0.98480775, -0.64278761 /
        DATA  C162,C163 /  0.38268343,  1.30656297 /
        DATA  C164,C165 /  0.54119610,  0.92387953 /
    C
    C-----------------NESTED LOOPS----------------------------------
    C
        DO 10 K=1, M
               N1 = NI(K)
           N2 = N/N1
           L  = 1
           N3 = N2 - N1*(N2/N1)
           DO 15 J = 1, N1
              LP(J) = L
                  L = L + N3
              IF (L.GT.N1) L = L - N1
       15      CONTINUE
    C
           DO 20 J=1, N, N1
              IT   = J
              DO 30 L=1, N1
                 I(L) = IT
                 IP(LP(L)) = IT
                 IT = IT + N2
                 IF (IT.GT.N)  IT = IT - N
       30         CONTINUE
              GOTO (20,102,103,104,105,20,107,108,109,
         +              20,20,20,20,20,20,116),N1

    C----------------WFTA N=2--------------------------------
    C
      102   R1       = X(I(1))
        X(I(1))  = R1 + X(I(2))
        X(I(2))  = R1 - X(I(2))
    C
        R1       = Y(I(1))
        Y(IP(1)) = R1 + Y(I(2))
        Y(IP(2)) = R1 - Y(I(2))
    C
        GOTO 20
    C
    C----------------WFTA N=3--------------------------------
    C
      103   R2 = (X(I(2)) - X(I(3))) * C31
        R1 =  X(I(2)) + X(I(3))
        X(I(1))= X(I(1)) + R1
        R1     = X(I(1)) + R1 * C32
    C
        S2 = (Y(I(2)) - Y(I(3))) * C31
        S1 =  Y(I(2)) + Y(I(3))
        Y(I(1))= Y(I(1)) + S1
        S1     = Y(I(1)) + S1 * C32
    C
        X(IP(2)) = R1 - S2
        X(IP(3)) = R1 + S2
        Y(IP(2)) = S1 + R2
        Y(IP(3)) = S1 - R2
    C
        GOTO 20
    C
    C----------------WFTA N=4---------------------------------
    C
      104   R1 = X(I(1)) + X(I(3))
        T1 = X(I(1)) - X(I(3))
        R2 = X(I(2)) + X(I(4))
        X(IP(1)) = R1 + R2
        X(IP(3)) = R1 - R2
    C
        R1 = Y(I(1)) + Y(I(3))
        T2 = Y(I(1)) - Y(I(3))
        R2 = Y(I(2)) + Y(I(4))
        Y(IP(1)) = R1 + R2
        Y(IP(3)) = R1 - R2
    C
        R1 = X(I(2)) - X(I(4))
        R2 = Y(I(2)) - Y(I(4))
    C
        X(IP(2)) = T1 + R2
        X(IP(4)) = T1 - R2
        Y(IP(2)) = T2 - R1
        Y(IP(4)) = T2 + R1
    C
        GOTO 20

    C----------------WFTA N=5--------------------------------
    C
      105   R1 = X(I(2)) + X(I(5))
        R4 = X(I(2)) - X(I(5))
        R3 = X(I(3)) + X(I(4))
        R2 = X(I(3)) - X(I(4))
    C
        T = (R1 - R3) * C54
        R1 = R1 + R3
        X(I(1)) = X(I(1)) + R1
        R1      = X(I(1)) + R1 * C55
    C
        R3 = R1 - T
        R1 = R1 + T
    C
        T = (R4 + R2) * C51
        R4 =  T + R4 * C52
        R2 =  T + R2 * C53
    C
        S1 = Y(I(2)) + Y(I(5))
        S4 = Y(I(2)) - Y(I(5))
        S3 = Y(I(3)) + Y(I(4))
        S2 = Y(I(3)) - Y(I(4))
    C
        T = (S1 - S3) * C54
        S1 = S1 + S3
        Y(I(1)) = Y(I(1)) + S1
        S1      = Y(I(1)) + S1 * C55
    C
        S3 =  S1 - T
        S1 =  S1 + T
    C
        T = (S4 + S2) * C51
        S4 =  T + S4 * C52
        S2 =  T + S2 * C53
    C
        X(IP(2)) = R1 + S2
        X(IP(5)) = R1 - S2
        X(IP(3)) = R3 - S4
        X(IP(4)) = R3 + S4
    C
        Y(IP(2)) = S1 - R2
        Y(IP(5)) = S1 + R2
        Y(IP(3)) = S3 + R4
        Y(IP(4)) = S3 - R4
    C
        GOTO 20

    C-----------------WFTA N=7--------------------------
    C
      107   R1 = X(I(2)) + X(I(7))
        R6 = X(I(2)) - X(I(7))
        S1 = Y(I(2)) + Y(I(7))
        S6 = Y(I(2)) - Y(I(7))
        R2 = X(I(3)) + X(I(6))
        R5 = X(I(3)) - X(I(6))
        S2 = Y(I(3)) + Y(I(6))
        S5 = Y(I(3)) - Y(I(6))
        R3 = X(I(4)) + X(I(5))
        R4 = X(I(4)) - X(I(5))
        S3 = Y(I(4)) + Y(I(5))
        S4 = Y(I(4)) - Y(I(5))
    C
        T3 = (R1 - R2) * C74
        T  = (R1 - R3) * C72
        R1 = R1 + R2 + R3
        X(I(1)) = X(I(1)) + R1
        R1      = X(I(1)) + R1 * C71
        R2 =(R3 - R2) * C73
        R3 = R1 - T + R2
        R2 = R1 - R2 - T3
        R1 = R1 + T + T3
        T = (R6 - R5) * C78
        T3 =(R6 + R4) * C76
        R6 =(R6 + R5 - R4) * C75
        R5 =(R5 + R4) * C77
        R4 = R6 - T3 + R5
        R5 = R6 - R5 - T
        R6 = R6 + T3 + T
    C
        T3 = (S1 - S2) * C74
        T  = (S1 - S3) * C72
        S1 =  S1 + S2 + S3
        Y(I(1)) = Y(I(1)) + S1
        S1      = Y(I(1)) + S1 * C71
        S2 =(S3 - S2) * C73
        S3 = S1 - T  + S2
        S2 = S1 - S2 - T3
        S1 = S1 + T  + T3
        T  = (S6 - S5) * C78
        T3 = (S6 + S4) * C76
        S6 = (S6 + S5 - S4) * C75
        S5 = (S5 + S4) * C77
        S4 = S6 - T3 + S5
        S5 = S6 - S5 - T
        S6 = S6 + T3 + T
    C
        X(IP(2)) = R3 + S4
        X(IP(7)) = R3 - S4
        X(IP(3)) = R1 + S6
        X(IP(6)) = R1 - S6
        X(IP(4)) = R2 - S5
        X(IP(5)) = R2 + S5
        Y(IP(4)) = S2 + R5
        Y(IP(5)) = S2 - R5
        Y(IP(2)) = S3 - R4
        Y(IP(7)) = S3 + R4
        Y(IP(3)) = S1 - R6
        Y(IP(6)) = S1 + R6
    C
        GOTO 20

    C-----------------WFTA N=8--------------------------
    C
      108   R1 = X(I(1)) + X(I(5))
        R2 = X(I(1)) - X(I(5))
        R3 = X(I(2)) + X(I(8))
        R4 = X(I(2)) - X(I(8))
        R5 = X(I(3)) + X(I(7))
        R6 = X(I(3)) - X(I(7))
        R7 = X(I(4)) + X(I(6))
        R8 = X(I(4)) - X(I(6))
        T1 = R1 + R5
        T2 = R1 - R5
        T3 = R3 + R7
        R3 =(R3 - R7) * C81
        X(IP(1)) = T1 + T3
        X(IP(5)) = T1 - T3
        T1 = R2 + R3
        T3 = R2 - R3
        S1 = R4 - R8
        R4 =(R4 + R8) * C81
        S2 = R4 + R6
        S3 = R4 - R6
        R1 = Y(I(1)) + Y(I(5))
        R2 = Y(I(1)) - Y(I(5))
        R3 = Y(I(2)) + Y(I(8))
        R4 = Y(I(2)) - Y(I(8))
        R5 = Y(I(3)) + Y(I(7))
        R6 = Y(I(3)) - Y(I(7))
        R7 = Y(I(4)) + Y(I(6))
        R8 = Y(I(4)) - Y(I(6))
        T4 = R1 + R5
        R1 = R1 - R5
        R5 = R3 + R7
        R3 =(R3 - R7) * C81
        Y(IP(1)) = T4 + R5
        Y(IP(5)) = T4 - R5
        R5 = R2 + R3
        R2 = R2 - R3
        R3 = R4 - R8
        R4 =(R4 + R8) * C81
        R7 = R4 + R6
        R4 = R4 - R6
        X(IP(2)) = T1 + R7
        X(IP(8)) = T1 - R7
        X(IP(3)) = T2 + R3
        X(IP(7)) = T2 - R3
        X(IP(4)) = T3 + R4
        X(IP(6)) = T3 - R4
        Y(IP(2)) = R5 - S2
        Y(IP(8)) = R5 + S2
        Y(IP(3)) = R1 - S1
        Y(IP(7)) = R1 + S1
        Y(IP(4)) = R2 - S3
        Y(IP(6)) = R2 + S3
    C
        GOTO 20

    C-----------------WFTA N=9-----------------------
    C
      109   R1 = X(I(2)) + X(I(9))
        R2 = X(I(2)) - X(I(9))
        R3 = X(I(3)) + X(I(8))
        R4 = X(I(3)) - X(I(8))
        R5 = X(I(4)) + X(I(7))
        T8 =(X(I(4)) - X(I(7))) * C31
        R7 = X(I(5)) + X(I(6))
        R8 = X(I(5)) - X(I(6))
        T0 = X(I(1)) + R5
        T7 = X(I(1)) + R5 * C95
        R5 = R1 + R3 + R7
        X(I(1)) = T0 + R5
        T5 = T0 + R5 * C95
        T3 = (R3 - R7) * C92
        R7 = (R1 - R7) * C93
        R3 = (R1 - R3) * C94
        T1 = T7 + T3 + R3
        T3 = T7 - T3 - R7
        T7 = T7 + R7 - R3
        T6 = (R2 - R4 + R8) * C31
        T4 = (R4 + R8) * C96
        R8 = (R2 - R8) * C97
        R2 = (R2 + R4) * C98
        T2 = T8 + T4 + R2
        T4 = T8 - T4 - R8
        T8 = T8 + R8 - R2
    C
        R1 = Y(I(2)) + Y(I(9))
        R2 = Y(I(2)) - Y(I(9))
        R3 = Y(I(3)) + Y(I(8))
        R4 = Y(I(3)) - Y(I(8))
        R5 = Y(I(4)) + Y(I(7))
        R6 =(Y(I(4)) - Y(I(7))) * C31
        R7 = Y(I(5)) + Y(I(6))
        R8 = Y(I(5)) - Y(I(6))
        T0 = Y(I(1)) + R5
        T9 = Y(I(1)) + R5 * C95
        R5 = R1 + R3 + R7
        Y(I(1)) = T0 + R5
        R5 = T0 + R5 * C95
        T0 = (R3 - R7) * C92
        R7 = (R1 - R7) * C93
        R3 = (R1 - R3) * C94
        R1 = T9 + T0 + R3
        T0 = T9 - T0 - R7
        R7 = T9 + R7 - R3
        R9 = (R2 - R4 + R8) * C31
        R3 = (R4 + R8) * C96
        R8 = (R2 - R8) * C97
        R4 = (R2 + R4) * C98
        R2 = R6 + R3 + R4
        R3 = R6 - R8 - R3
        R8 = R6 + R8 - R4
    C
        X(IP(2)) = T1 - R2
        X(IP(9)) = T1 + R2
        Y(IP(2)) = R1 + T2
        Y(IP(9)) = R1 - T2
        X(IP(3)) = T3 + R3
        X(IP(8)) = T3 - R3
        Y(IP(3)) = T0 - T4
        Y(IP(8)) = T0 + T4
        X(IP(4)) = T5 - R9
        X(IP(7)) = T5 + R9
        Y(IP(4)) = R5 + T6
        Y(IP(7)) = R5 - T6
        X(IP(5)) = T7 - R8
        X(IP(6)) = T7 + R8
        Y(IP(5)) = R7 + T8
        Y(IP(6)) = R7 - T8
    C
        GOTO 20

    C-----------------WFTA N=16------------------------
    C
      116   R1 = X(I(1)) + X(I(9))
        R2 = X(I(1)) - X(I(9))
        R3 = X(I(2)) + X(I(10))
        R4 = X(I(2)) - X(I(10))
        R5 = X(I(3)) + X(I(11))
        R6 = X(I(3)) - X(I(11))
        R7 = X(I(4)) + X(I(12))
        R8 = X(I(4)) - X(I(12))
        R9 = X(I(5)) + X(I(13))
        R10= X(I(5)) - X(I(13))
        R11 = X(I(6)) + X(I(14))
        R12 = X(I(6)) - X(I(14))
        R13 = X(I(7)) + X(I(15))
        R14 = X(I(7)) - X(I(15))
        R15 = X(I(8)) + X(I(16))
        R16 = X(I(8)) - X(I(16))
        T1 = R1 + R9
        T2 = R1 - R9
        T3 = R3 + R11
        T4 = R3 - R11
        T5 = R5 + R13
        T6 = R5 - R13
        T7 = R7 + R15
        T8 = R7 - R15
        R1 = T1 + T5
        R3 = T1 - T5
        R5 = T3 + T7
        R7 = T3 - T7
        X(IP( 1)) = R1 + R5
        X(IP( 9)) = R1 - R5
        T1 = C81 * (T4 + T8)
        T5 = C81 * (T4 - T8)
        R9 = T2 + T5
        R11= T2 - T5
        R13 = T6 + T1
        R15 = T6 - T1
        T1 = R4 + R16
        T2 = R4 - R16
        T3 = C81 * (R6 + R14)
        T4 = C81 * (R6 - R14)
        T5 = R8 + R12
        T6 = R8 - R12
        T7 = C162 * (T2 - T6)
        T2 = C163 * T2 - T7
        T6 = C164 * T6 - T7
        T7 = R2 + T4
        T8 = R2 - T4
        R2 = T7 + T2
        R4 = T7 - T2
        R6 = T8 + T6
        R8 = T8 - T6
        T7 = C165 * (T1 + T5)
        T2 = T7 - C164 * T1
        T4 = T7 - C163 * T5
        T6 = R10 + T3
        T8 = R10 - T3
        R10 = T6 + T2
        R12 = T6 - T2
        R14 = T8 + T4
        R16 = T8 - T4
        R1 = Y(I(1)) + Y(I(9))
        S2 = Y(I(1)) - Y(I(9))
        S3 = Y(I(2)) + Y(I(10))
        S4 = Y(I(2)) - Y(I(10))
        R5 = Y(I(3)) + Y(I(11))
        S6 = Y(I(3)) - Y(I(11))
        S7 = Y(I(4)) + Y(I(12))
        S8 = Y(I(4)) - Y(I(12))
        S9 = Y(I(5)) + Y(I(13))
        S10= Y(I(5)) - Y(I(13))
        S11 = Y(I(6)) + Y(I(14))
        S12 = Y(I(6)) - Y(I(14))
        S13 = Y(I(7)) + Y(I(15))
        S14 = Y(I(7)) - Y(I(15))
        S15 = Y(I(8)) + Y(I(16))
        S16 = Y(I(8)) - Y(I(16))
        T1 = R1 + S9
        T2 = R1 - S9
        T3 = S3 + S11
        T4 = S3 - S11
        T5 = R5 + S13
        T6 = R5 - S13
        T7 = S7 + S15
        T8 = S7 - S15
        R1 = T1 + T5
        S3 = T1 - T5
        R5 = T3 + T7
        S7 = T3 - T7
        Y(IP( 1)) = R1 + R5
        Y(IP( 9)) = R1 - R5
        X(IP( 5)) = R3 + S7
        X(IP(13)) = R3 - S7
        Y(IP( 5)) = S3 - R7
        Y(IP(13)) = S3 + R7
        T1 = C81 * (T4 + T8)
        T5 = C81 * (T4 - T8)
        S9 = T2 + T5
        S11= T2 - T5
        S13 = T6 + T1
        S15 = T6 - T1
        T1 = S4 + S16
        T2 = S4 - S16
        T3 = C81 * (S6 + S14)
        T4 = C81 * (S6 - S14)
        T5 = S8 + S12
        T6 = S8 - S12
        T7 = C162 * (T2 - T6)
        T2 = C163 * T2 - T7
        T6 = C164 * T6 - T7
        T7 = S2 + T4
        T8 = S2 - T4
        S2 = T7 + T2
        S4 = T7 - T2
        S6 = T8 + T6
        S8 = T8 - T6
        T7 = C165 * (T1 + T5)
        T2 = T7 - C164 * T1
        T4 = T7 - C163 * T5
        T6 = S10 + T3
        T8 = S10 - T3
        S10 = T6 + T2
        S12 = T6 - T2
        S14 = T8 + T4
        S16 = T8 - T4
        X(IP( 2)) = R2 + S10
        X(IP(16)) = R2 - S10
        Y(IP( 2)) = S2 - R10
        Y(IP(16)) = S2 + R10
        X(IP( 3)) = R9 + S13
        X(IP(15)) = R9 - S13
        Y(IP( 3)) = S9 - R13
        Y(IP(15)) = S9 + R13
        X(IP( 4)) = R8 - S16
        X(IP(14)) = R8 + S16
        Y(IP( 4)) = S8 + R16
        Y(IP(14)) = S8 - R16
        X(IP( 6)) = R6 + S14
        X(IP(12)) = R6 - S14
        Y(IP( 6)) = S6 - R14
        Y(IP(12)) = S6 + R14
        X(IP( 7)) = R11 - S15
        X(IP(11)) = R11 + S15
        Y(IP( 7)) = S11 + R15
        Y(IP(11)) = S11 - R15
        X(IP( 8)) = R4 - S12
        X(IP(10)) = R4 + S12
        Y(IP( 8)) = S4 + R12
        Y(IP(10)) = S4 - R12
    C
        GOTO 20
    C
       20      CONTINUE
       10   CONTINUE
        RETURN
        END
     

    Contributor

    • ContribEEBurrus