Summary: Radix-4, DIF, Three Butterfly FFT
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