| << Chapter < Page | Chapter >> Page > |
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 isfollowed by an unscrambler.
C---------------------------------------------------
CC A PRIME FACTOR FFT PROGRAM WITH GENERAL MODULES
C COMPLEX INPUT DATA IN ARRAYS X AND YC COMPLEX OUTPUT IN A AND B
C LENGTH N WITH M FACTORS IN ARRAY NIC N = NI(1)*NI(2)* ... *NI(M)
C UNSCRAMBLING CONSTANT UNSCC UNSC = N/NI(1) + N/NI(2) +...+ N/NI(M), MOD N
C C. S. BURRUS, RICE UNIVERSITY, JAN 1987C
C--------------------------------------------------C
SUBROUTINE PFA(X,Y,N,M,NI,A,B,UNSC)C
INTEGER NI(4), I(16), UNSCREAL X(1), Y(1), A(1), B(1)
CDATA 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, MN1 = NI(K)
N2 = N/N1DO 15 J=1, N, N1
IT = JDO 30 L=1, N1
I(L) = ITA(L) = X(IT)
B(L) = Y(IT)IT = IT + N2
IF (IT.GT.N) IT = IT - N30 CONTINUE
GOTO (20,102,103,104,105,20,107), N1C
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 20C----------------WFTA N=3--------------------------------
C103 R2 = (A(2) - A(3)) * C31
R1 = A(2) + A(3)A(1)= A(1) + R1
R1 = A(1) + R1 * C32C
S2 = (B(2) - B(3)) * C31S1 = B(2) + B(3)
B(1)= B(1) + S1S1 = B(1) + S1 * C32
CA(2) = R1 - S2
A(3) = R1 + S2B(2) = S1 + R2
B(3) = S1 - R2C
GOTO 20C
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 - R2C
R1 = B(1) + B(3)T2 = B(1) - B(3)
R2 = B(2) + B(4)B(1) = R1 + R2
B(3) = R1 - R2C
R1 = A(2) - A(4)R2 = B(2) - B(4)
CA(2) = T1 + R2
A(4) = T1 - R2B(2) = T2 - R1
B(4) = T2 + R1C
GOTO 20C
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)
CT = (R1 - R3) * C54
R1 = R1 + R3A(1) = A(1) + R1
R1 = A(1) + R1 * C55C
R3 = R1 - TR1 = R1 + T
CT = (R4 + R2) * C51
R4 = T + R4 * C52R2 = T + R2 * C53
CS1 = B(2) + B(5)
S4 = B(2) - B(5)S3 = B(3) + B(4)
S2 = B(3) - B(4)C
T = (S1 - S3) * C54S1 = S1 + S3
B(1) = B(1) + S1S1 = B(1) + S1 * C55
CS3 = S1 - T
S1 = S1 + TC
T = (S4 + S2) * C51S4 = T + S4 * C52
S2 = T + S2 * C53C
A(2) = R1 + S2A(5) = R1 - S2
A(3) = R3 - S4A(4) = R3 + S4
CB(2) = S1 - R2
B(5) = S1 + R2B(3) = S3 + R4
B(4) = S3 - R4C
GOTO 20C-----------------WFTA N=7--------------------------
C107 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) * C74T = (R1 - R3) * C72
R1 = R1 + R2 + R3A(1) = A(1) + R1
R1 = A(1) + R1 * C71R2 =(R3 - R2) * C73
R3 = R1 - T + R2R2 = R1 - R2 - T3
R1 = R1 + T + T3T = (R6 - R5) * C78
T3 =(R6 + R4) * C76R6 =(R6 + R5 - R4) * C75
R5 =(R5 + R4) * C77R4 = R6 - T3 + R5
R5 = R6 - R5 - TR6 = R6 + T3 + T
CT3 = (S1 - S2) * C74
T = (S1 - S3) * C72S1 = S1 + S2 + S3
B(1) = B(1) + S1S1 = B(1) + S1 * C71
S2 =(S3 - S2) * C73S3 = S1 - T + S2
S2 = S1 - S2 - T3S1 = S1 + T + T3
T = (S6 - S5) * C78T3 = (S6 + S4) * C76
S6 = (S6 + S5 - S4) * C75S5 = (S5 + S4) * C77
S4 = S6 - T3 + S5S5 = S6 - S5 - T
S6 = S6 + T3 + TC
A(2) = R3 + S4A(7) = R3 - S4
A(3) = R1 + S6A(6) = R1 - S6
A(4) = R2 - S5A(5) = R2 + S5
B(4) = S2 + R5B(5) = S2 - R5
B(2) = S3 - R4B(7) = S3 + R4
B(3) = S1 - R6B(6) = S1 + R6
C20 IT = J
DO 31 L=1, N1I(L) = IT
X(IT) = A(L)Y(IT) = B(L)
IT = IT + N2IF (IT.GT.N) IT = IT - N
31 CONTINUE15 CONTINUE
10 CONTINUEC
C-----------------UNSCRAMBLING----------------------C
L = 1DO 2 K=1, N
A(K) = X(L)B(K) = Y(L)
L = L + UNSCIF (L.GT.N) L = L - N
2 CONTINUERETURN
ENDC
Notification Switch
Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|