Program 11: Split-Radix, DIF, Two-Butterfly, FFT
Summary: Very efficient Split-Radix DIF, Two Butterfly FFT
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
References
Comments, questions, feedback, criticisms?
Send feedback
- E-mail the author of the module, Program 11: Split-Radix, DIF, Two-Butterfly, FFT
Footer
More about this content: Metadata | Version History
This work is licensed by C. Sidney Burrus under a Creative Commons Attribution License (CC-BY 2.0).
Last edited by C. Sidney Burrus on Sep 2, 2008 2:10 pm GMT-5.