Skip to main content
Engineering LibreTexts

14.3: Appendix 3 - FFT Computer Programs

  1. Last updated
  2. Save as PDF
  • Page ID
    2051

    • \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

    \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

    \( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)

    ( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)

    \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

    \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)

    \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

    \( \newcommand{\Span}{\mathrm{span}}\)

    \( \newcommand{\id}{\mathrm{id}}\)

    \( \newcommand{\Span}{\mathrm{span}}\)

    \( \newcommand{\kernel}{\mathrm{null}\,}\)

    \( \newcommand{\range}{\mathrm{range}\,}\)

    \( \newcommand{\RealPart}{\mathrm{Re}}\)

    \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

    \( \newcommand{\Argument}{\mathrm{Arg}}\)

    \( \newcommand{\norm}[1]{\| #1 \|}\)

    \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

    \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)

    \( \newcommand{\vectorA}[1]{\vec{#1}}      % arrow\)

    \( \newcommand{\vectorAt}[1]{\vec{\text{#1}}}      % arrow\)

    \( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

    \( \newcommand{\vectorC}[1]{\textbf{#1}} \)

    \( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)

    \( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)

    \( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)

    \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

    \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

    \(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)

    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

    This page titled 14.3: Appendix 3 - FFT Computer Programs is shared under a CC BY license and was authored, remixed, and/or curated by C. Sidney Burrus.

    1. Back to top
    • Was this article helpful?

    Recommended articles

    1. Article type
      Section or Page
      Author
      C. Sidney Burrus
      License
      CC BY
      OER program or Publisher
      OpenStax CNX
      Show TOC
      no
    2. Tags
      This page has no tags.