| << Chapter < Page | Chapter >> Page > |
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 YC LENGTH IS N = 4 ** M
C C. S. BURRUS, RICE UNIVERSITY, SEPT 1983C---------------------------------------------------------
SUBROUTINE FFT4 (X,Y,N,M)REAL X(1), Y(1)
C--------------MAIN FFT LOOPS-----------------------------N2 = N
DO 10 K = 1, MN1 = N2
N2 = N2/4E = 6.283185307179586/N1
A = 0C--------------------MAIN BUTTERFLIES-------------------
DO 20 J=1, N2B = A + A
C = A + BCO1 = 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 + N2I2 = I1 + N2
I3 = I2 + N2R1 = 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 - R2R1 = R3 - S4
R3 = R3 + S4Y(I) = S1 + S2
S2 = S1 - S2S1 = S3 + R4
S3 = S3 - R4X(I1) = CO1*R3 + SI1*S3
Y(I1) = CO1*S3 - SI1*R3X(I2) = CO2*R2 + SI2*S2
Y(I2) = CO2*S2 - SI2*R2X(I3) = CO3*R1 + SI3*S1
Y(I3) = CO3*S1 - SI3*R130 CONTINUE
20 CONTINUE10 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 PROGRAMC THREE BF, MULTIPLICATIONS BY 1, J, ETC. ARE REMOVED
C COMPLEX INPUT DATA IN ARRAYS X AND YC LENGTH IS N = 4 ** M
C TABLE LOOKUP OF W VALUESC
C C. S. BURRUS, RICE UNIVERSITY, SEPT 1983C
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 = NDO 10 K = 1, M
N1 = N2N2 = N2/4
JT = N2/2 + 1C---------------SPECIAL BUTTERFLY FOR W = 1---------------
DO 1 I = 1, N, N1I1 = I + N2
I2 = I1 + N2I3 = 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)
CX(I) = R1 + R2
X(I2)= R1 - R2X(I3)= R3 - S4
X(I1)= R3 + S4C
Y(I) = S1 + S2Y(I2)= S1 - S2
Y(I3)= S3 + R4Y(I1)= S3 - R4
C1 CONTINUE
IF (K.EQ.M) GOTO 10IE = N/N1
IA1 = 1C--------------GENERAL BUTTERFLY-----------------DO 20 J = 2, N2
IA1 = IA1 + IEIF (J.EQ.JT) GOTO 50
IA2 = IA1 + IA1 - 1IA3 = 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 + N2I2 = I1 + N2
I3 = I2 + N2R1 = 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)
CX(I) = R1 + R2
R2 = R1 - R2R1 = R3 - S4
R3 = R3 + S4C
Y(I) = S1 + S2S2 = S1 - S2
S1 = S3 + R4S3 = S3 - R4
CX(I1) = CO1*R3 + SI1*S3
Y(I1) = CO1*S3 - SI1*R3X(I2) = CO2*R2 + SI2*S2
Y(I2) = CO2*S2 - SI2*R2X(I3) = CO3*R1 + SI3*S1
Y(I3) = CO3*S1 - SI3*R130 CONTINUE
GOTO 20C------------------SPECIAL BUTTERFLY FOR W = J-----------
50 DO 40 I = J, N, N1I1 = I + N2
I2 = I1 + N2I3 = 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)
CX(I) = R1 + R2
Y(I2)=-R1 + R2R1 = R3 - S4
R3 = R3 + S4C
Y(I) = S1 + S2X(I2)= S1 - S2
S1 = S3 + R4S3 = S3 - R4
CX(I1) = (S3 + R3)*C21
Y(I1) = (S3 - R3)*C21X(I3) = (S1 - R1)*C21
Y(I3) =-(S1 + R1)*C2140 CONTINUE
20 CONTINUE10 CONTINUE
C-----------DIGIT REVERSE COUNTER----------100 J = 1
N1 = N - 1DO 104 I = 1, N1
IF (I.GE.J) GOTO 101R1 = X(J)
X(J) = X(I)X(I) = R1
R1 = Y(J)Y(J) = Y(I)
Y(I) = R1101 K = N/4
102 IF (K*3.GE.J) GOTO 103J = J - K*3
K = K/4GOTO 102
103 J = J + K104 CONTINUE
RETURNEND
Notification Switch
Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|