A FORTRAN implementation of the first-order Goertzel algorithm with in-order input as given in ((Reference)) and [1] is given below.
|
Skip to content Skip to navigation Skip to collection information
Summary: Fortran programs for efficient DFT, Cooley-Tukey, and Prime Factor Algorithm FFTs.
A FORTRAN implementation of the first-order Goertzel algorithm with in-order input as given in ((Reference)) and [1] is given below.
|
Below is the program for a second order Goertzel algorithm.
|
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
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.
|
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
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
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
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
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
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
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
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
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
More about this module: Metadata | Downloads | Version History
This work is licensed by C. Sidney Burrus under a Creative Commons Attribution License (CC-BY 2.0), and is an Open Educational Resource.
Last edited by Daniel Williamson on Sep 18, 2009 2:45 pm -0500.
"The Fast Fourier Transform (FFT) is a landmark algorithm used in fields ranging from signal processing to high-performance computing. First popularized by two American scientists in 1965, the […]"