| << Chapter < Page | Chapter >> Page > |
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, 1984C COMPLEX INPUT DATA IN ARRAYS X AND Y
C LENGTH IS N = 2 ** MC C. S. BURRUS, RICE UNIVERSITY, MARCH 1984
CC---------------------------------------------------------
SUBROUTINE FFT (X,Y,N,M)REAL X(1), Y(1)
C--------------MAIN FFT LOOPS-----------------------------C
N1 = NN2 = N/2
IP = 0IS = 1
A = 6.283185307179586/NDO 10 K = 1, M-1
JD = N1 + N2N1 = N2
N2 = N2/2J0 = N1*IP + 1
IP = 1 - IPDO 20 J = J0, N, JD
JS = 0JT = J + N2 - 1
DO 30 I = J, JTJSS= JS*IS
JS = JS + 1C1 = COS(A*JSS)
C3 = COS(3*A*JSS)S1 = -SIN(A*JSS)
S3 = -SIN(3*A*JSS)I1 = I + N2
I2 = I1 + N2I3 = 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) = R1C
R1 = Y(I ) + Y(I2)R4 = Y(I ) - Y(I2)
R5 = Y(I1) - Y(I3)Y(I2) = Y(I1) + Y(I3)
Y(I1) = R1C
R1 = R2 - R5R2 = R2 + R5
R5 = R4 + R3R4 = R4 - R3
CX(I) = C1*R1 + S1*R5
Y(I) = C1*R5 - S1*R1X(I3) = C3*R2 + S3*R4
Y(I3) = C3*R4 - S3*R230 CONTINUE
20 CONTINUEIS = IS + IS
10 CONTINUEIP = 1 - IP
J0 = 2 - IPDO 5 I = J0, N-1, 3
I1 = I + 1R1 = 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) = R15 CONTINUE
RETURNEND
Below is the Fortran code for a simple Decimation-in-Frequency, Split-Radix, two butterfly FFT to be followed by a bit-reversing unscrambler. Twiddlefactors are precalculated and stored in arrays WR and WI.
C--------------------------------------------------------------C
C A DUHAMEL-HOLLMAN SPLIT RADIX FFT CC REF: ELECTRONICS LETTERS, JAN. 5, 1984 C
C COMPLEX INPUT AND OUTPUT DATA IN ARRAYS X AND Y CC LENGTH IS N = 2 ** M, OUTPUT IN BIT-REVERSED ORDER C
C TWO BUTTERFLIES TO REMOVE MULTS BY UNITY CC SPECIAL LAST TWO STAGES C
C TABLE LOOK-UP OF SINE AND COSINE VALUES CC C.S. BURRUS, RICE UNIV. APRIL 1985 C
C--------------------------------------------------------------CC
SUBROUTINE FFT(X,Y,N,M,WR,WI)REAL X(1),Y(1),WR(1),WI(1)
C81= 0.707106778N2 = 2*N
DO 10 K = 1, M-3IS = 1
ID = N2N2 = N2/2
N4 = N2/42 DO 1 I0 = IS, N-1, ID
I1 = I0 + N4I2 = I1 + N4
I3 = I2 + N4R1 = 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 - R2R1 = 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 + R2Y(I3) = R1 + R2
1 CONTINUEIS = 2*ID - N2 + 1
ID = 4*IDIF (IS.LT.N) GOTO 2
IE = N/N2IA1 = 1
DO 20 J = 2, N4IA1 = IA1 + IE
IA3 = 3*IA1 - 2CC1 = WR(IA1)
SS1 = WI(IA1)CC3 = WR(IA3)
SS3 = WI(IA3)IS = J
ID = 2*N240 DO 30 I0 = IS, N-1, ID
I1 = I0 + N4I2 = I1 + N4
I3 = I2 + N4C
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)
CS3 = R1 - S2
R1 = R1 + S2S2 = R2 - S1
R2 = R2 + S1X(I2) = R1*CC1 - S2*SS1
Y(I2) =-S2*CC1 - R1*SS1X(I3) = S3*CC3 + R2*SS3
Y(I3) = R2*CC3 - S3*SS330 CONTINUE
IS = 2*ID - N2 + JID = 4*ID
IF (IS.LT.N) GOTO 4020 CONTINUE
10 CONTINUEC
IS = 1ID = 32
50 DO 60 I = IS, N, IDI0 = I + 8
DO 15 J = 1, 2R1 = 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 - R2R1 = 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 - R2Y(I0+2) = S3 - R4
Y(I0+3) = S3 + R4X(I0+2) = R3 + S4
X(I0+3) = R3 - S4I0 = I0 + 4
15 CONTINUE60 CONTINUE
IS = 2*ID - 15ID = 4*ID
IF (IS.LT.N) GOTO 50C
IS = 1ID = 16
55 DO 65 I0 = IS, N, IDR1 = 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 - R3R1 = R1 + R3
R3 = R2 - R4R2 = R2 + R4
X(I0) = R1 + R2X(I0+1) = R1 - R2
CR1 = 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 + S3S3 = R2 - R4
R2 = R2 + R4Y(I0) = R1 + R2
Y(I0+1) = R1 - R2X(I0+2) = T1 + S3
X(I0+3) = T1 - S3Y(I0+2) = T2 - R3
Y(I0+3) = T2 + R3C
R1 = (R6 - R8)*C81R6 = (R6 + R8)*C81
R2 = (S6 - S8)*C81S6 = (S6 + S8)*C81
CT1 = R5 - R1
R5 = R5 + R1R8 = R7 - R6
R7 = R7 + R6T2 = S5 - R2
S5 = S5 + R2S8 = S7 - S6
S7 = S7 + S6X(I0+4) = R5 + S7
X(I0+7) = R5 - S7X(I0+5) = T1 + S8
X(I0+6) = T1 - S8Y(I0+4) = S5 - R7
Y(I0+7) = S5 + R7Y(I0+5) = T2 - R8
Y(I0+6) = T2 + R865 CONTINUE
IS = 2*ID - 7ID = 4*ID
IF (IS.LT.N) GOTO 55C
C------------BIT REVERSE COUNTER-----------------C
100 J = 1N1 = N - 1
DO 104 I=1, N1IF (I.GE.J) GOTO 101
XT = X(J)X(J) = X(I)
X(I) = XTXT = Y(J)
Y(J) = Y(I)Y(I) = XT
101 K = N/2102 IF (K.GE.J) GOTO 103
J = J - KK = K/2
GOTO 102103 J = J + K
104 CONTINUERETURN
END
Notification Switch
Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|