Appendix 3: FFT Computer Programs 1.2 2008/09/02 16:34:54 GMT-5 2008/09/04 21:12:15.384 GMT-5 C. Sidney Burrus csb@rice.edu C. Sidney Burrus csb@rice.edu Daniel Collins Williamson dwilliamson1285@gmail.com FFT FFT programs FORTRAN PFA Fortran programs for efficient DFT, Cooley-Tukey, and Prime Factor Algorithm FFTs.

Goertzel Algorithm A FORTRAN implementation of the first-order Goertzel algorithm with in-order input as given in () and 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 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 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 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

Burrus, C. S. and Parks, T. W. DFT/FFT and Convolution Algorithms John Wiley & Sons 1985 New York