| << 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, 7, 8, 9, and 16. It is bothin-place and in-order, so requires no unscrambler.
C
C A PRIME FACTOR FFT PROGRAMC IN-PLACE AND IN-ORDER
C COMPLEX INPUT DATA IN ARRAYS X AND YC LENGTH N WITH M FACTORS IN ARRAY NI
C N = NI(1)*NI(2)*...*NI(M)C REDUCED TEMP STORAGE IN SHORT WFTA MODULES
C Has modules 2,3,4,5,7,8,9,16C PROGRAM BY C. S. BURRUS, RICE UNIVERSITY
C SEPT 1983C----------------------------------------------------
CSUBROUTINE PFA(X,Y,N,M,NI)
INTEGER NI(4), I(16), IP(16), LP(16)REAL X(1), Y(1)
DATA 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 /DATA C81 / 0.70710678 /
DATA C95 / -0.50000000 /DATA C92, C93 / 0.93969262, -0.17364818 /
DATA C94, C96 / 0.76604444, -0.34202014 /DATA C97, C98 / -0.98480775, -0.64278761 /
DATA C162,C163 / 0.38268343, 1.30656297 /DATA C164,C165 / 0.54119610, 0.92387953 /
CC-----------------NESTED LOOPS----------------------------------
CDO 10 K=1, M
N1 = NI(K)N2 = N/N1
L = 1N3 = N2 - N1*(N2/N1)
DO 15 J = 1, N1LP(J) = L
L = L + N3IF (L.GT.N1) L = L - N1
15 CONTINUEC
DO 20 J=1, N, N1IT = J
DO 30 L=1, N1I(L) = IT
IP(LP(L)) = ITIT = IT + N2
IF (IT.GT.N) IT = IT - N30 CONTINUE
GOTO (20,102,103,104,105,20,107,108,109,+ 20,20,20,20,20,20,116),N1C----------------WFTA N=2--------------------------------
C102 R1 = X(I(1))
X(I(1)) = R1 + X(I(2))X(I(2)) = R1 - X(I(2))
CR1 = Y(I(1))
Y(IP(1)) = R1 + Y(I(2))Y(IP(2)) = R1 - Y(I(2))
CGOTO 20
CC----------------WFTA N=3--------------------------------
C103 R2 = (X(I(2)) - X(I(3))) * C31
R1 = X(I(2)) + X(I(3))X(I(1))= X(I(1)) + R1
R1 = X(I(1)) + R1 * C32C
S2 = (Y(I(2)) - Y(I(3))) * C31S1 = Y(I(2)) + Y(I(3))
Y(I(1))= Y(I(1)) + S1S1 = Y(I(1)) + S1 * C32
CX(IP(2)) = R1 - S2
X(IP(3)) = R1 + S2Y(IP(2)) = S1 + R2
Y(IP(3)) = S1 - R2C
GOTO 20C
C----------------WFTA N=4---------------------------------C
104 R1 = X(I(1)) + X(I(3))T1 = X(I(1)) - X(I(3))
R2 = X(I(2)) + X(I(4))X(IP(1)) = R1 + R2
X(IP(3)) = R1 - R2C
R1 = Y(I(1)) + Y(I(3))T2 = Y(I(1)) - Y(I(3))
R2 = Y(I(2)) + Y(I(4))Y(IP(1)) = R1 + R2
Y(IP(3)) = R1 - R2C
R1 = X(I(2)) - X(I(4))R2 = Y(I(2)) - Y(I(4))
CX(IP(2)) = T1 + R2
X(IP(4)) = T1 - R2Y(IP(2)) = T2 - R1
Y(IP(4)) = T2 + R1C
GOTO 20C----------------WFTA N=5--------------------------------C
105 R1 = X(I(2)) + X(I(5))R4 = X(I(2)) - X(I(5))
R3 = X(I(3)) + X(I(4))R2 = X(I(3)) - X(I(4))
CT = (R1 - R3) * C54
R1 = R1 + R3X(I(1)) = X(I(1)) + R1
R1 = X(I(1)) + R1 * C55C
R3 = R1 - TR1 = R1 + T
CT = (R4 + R2) * C51
R4 = T + R4 * C52R2 = T + R2 * C53
CS1 = Y(I(2)) + Y(I(5))
S4 = Y(I(2)) - Y(I(5))S3 = Y(I(3)) + Y(I(4))
S2 = Y(I(3)) - Y(I(4))C
T = (S1 - S3) * C54S1 = S1 + S3
Y(I(1)) = Y(I(1)) + S1S1 = Y(I(1)) + S1 * C55
CS3 = S1 - T
S1 = S1 + TC
T = (S4 + S2) * C51S4 = T + S4 * C52
S2 = T + S2 * C53C
X(IP(2)) = R1 + S2X(IP(5)) = R1 - S2
X(IP(3)) = R3 - S4X(IP(4)) = R3 + S4
CY(IP(2)) = S1 - R2
Y(IP(5)) = S1 + R2Y(IP(3)) = S3 + R4
Y(IP(4)) = S3 - R4C
GOTO 20C-----------------WFTA N=7--------------------------C
107 R1 = X(I(2)) + X(I(7))R6 = X(I(2)) - X(I(7))
S1 = Y(I(2)) + Y(I(7))S6 = Y(I(2)) - Y(I(7))
R2 = X(I(3)) + X(I(6))R5 = X(I(3)) - X(I(6))
S2 = Y(I(3)) + Y(I(6))S5 = Y(I(3)) - Y(I(6))
R3 = X(I(4)) + X(I(5))R4 = X(I(4)) - X(I(5))
S3 = Y(I(4)) + Y(I(5))S4 = Y(I(4)) - Y(I(5))
CT3 = (R1 - R2) * C74
T = (R1 - R3) * C72R1 = R1 + R2 + R3
X(I(1)) = X(I(1)) + R1R1 = X(I(1)) + R1 * C71
R2 =(R3 - R2) * C73R3 = R1 - T + R2
R2 = R1 - R2 - T3R1 = R1 + T + T3
T = (R6 - R5) * C78T3 =(R6 + R4) * C76
R6 =(R6 + R5 - R4) * C75R5 =(R5 + R4) * C77
R4 = R6 - T3 + R5R5 = R6 - R5 - T
R6 = R6 + T3 + TC
T3 = (S1 - S2) * C74T = (S1 - S3) * C72
S1 = S1 + S2 + S3Y(I(1)) = Y(I(1)) + S1
S1 = Y(I(1)) + S1 * C71S2 =(S3 - S2) * C73
S3 = S1 - T + S2S2 = S1 - S2 - T3
S1 = S1 + T + T3T = (S6 - S5) * C78
T3 = (S6 + S4) * C76S6 = (S6 + S5 - S4) * C75
S5 = (S5 + S4) * C77S4 = S6 - T3 + S5
S5 = S6 - S5 - TS6 = S6 + T3 + T
CX(IP(2)) = R3 + S4
X(IP(7)) = R3 - S4X(IP(3)) = R1 + S6
X(IP(6)) = R1 - S6X(IP(4)) = R2 - S5
X(IP(5)) = R2 + S5Y(IP(4)) = S2 + R5
Y(IP(5)) = S2 - R5Y(IP(2)) = S3 - R4
Y(IP(7)) = S3 + R4Y(IP(3)) = S1 - R6
Y(IP(6)) = S1 + R6C
GOTO 20C-----------------WFTA N=8--------------------------C
108 R1 = X(I(1)) + X(I(5))R2 = X(I(1)) - X(I(5))
R3 = X(I(2)) + X(I(8))R4 = X(I(2)) - X(I(8))
R5 = X(I(3)) + X(I(7))R6 = X(I(3)) - X(I(7))
R7 = X(I(4)) + X(I(6))R8 = X(I(4)) - X(I(6))
T1 = R1 + R5T2 = R1 - R5
T3 = R3 + R7R3 =(R3 - R7) * C81
X(IP(1)) = T1 + T3X(IP(5)) = T1 - T3
T1 = R2 + R3T3 = R2 - R3
S1 = R4 - R8R4 =(R4 + R8) * C81
S2 = R4 + R6S3 = R4 - R6
R1 = Y(I(1)) + Y(I(5))R2 = Y(I(1)) - Y(I(5))
R3 = Y(I(2)) + Y(I(8))R4 = Y(I(2)) - Y(I(8))
R5 = Y(I(3)) + Y(I(7))R6 = Y(I(3)) - Y(I(7))
R7 = Y(I(4)) + Y(I(6))R8 = Y(I(4)) - Y(I(6))
T4 = R1 + R5R1 = R1 - R5
R5 = R3 + R7R3 =(R3 - R7) * C81
Y(IP(1)) = T4 + R5Y(IP(5)) = T4 - R5
R5 = R2 + R3R2 = R2 - R3
R3 = R4 - R8R4 =(R4 + R8) * C81
R7 = R4 + R6R4 = R4 - R6
X(IP(2)) = T1 + R7X(IP(8)) = T1 - R7
X(IP(3)) = T2 + R3X(IP(7)) = T2 - R3
X(IP(4)) = T3 + R4X(IP(6)) = T3 - R4
Y(IP(2)) = R5 - S2Y(IP(8)) = R5 + S2
Y(IP(3)) = R1 - S1Y(IP(7)) = R1 + S1
Y(IP(4)) = R2 - S3Y(IP(6)) = R2 + S3
CGOTO 20C-----------------WFTA N=9-----------------------
C109 R1 = X(I(2)) + X(I(9))
R2 = X(I(2)) - X(I(9))R3 = X(I(3)) + X(I(8))
R4 = X(I(3)) - X(I(8))R5 = X(I(4)) + X(I(7))
T8 =(X(I(4)) - X(I(7))) * C31R7 = X(I(5)) + X(I(6))
R8 = X(I(5)) - X(I(6))T0 = X(I(1)) + R5
T7 = X(I(1)) + R5 * C95R5 = R1 + R3 + R7
X(I(1)) = T0 + R5T5 = T0 + R5 * C95
T3 = (R3 - R7) * C92R7 = (R1 - R7) * C93
R3 = (R1 - R3) * C94T1 = T7 + T3 + R3
T3 = T7 - T3 - R7T7 = T7 + R7 - R3
T6 = (R2 - R4 + R8) * C31T4 = (R4 + R8) * C96
R8 = (R2 - R8) * C97R2 = (R2 + R4) * C98
T2 = T8 + T4 + R2T4 = T8 - T4 - R8
T8 = T8 + R8 - R2C
R1 = Y(I(2)) + Y(I(9))R2 = Y(I(2)) - Y(I(9))
R3 = Y(I(3)) + Y(I(8))R4 = Y(I(3)) - Y(I(8))
R5 = Y(I(4)) + Y(I(7))R6 =(Y(I(4)) - Y(I(7))) * C31
R7 = Y(I(5)) + Y(I(6))R8 = Y(I(5)) - Y(I(6))
T0 = Y(I(1)) + R5T9 = Y(I(1)) + R5 * C95
R5 = R1 + R3 + R7Y(I(1)) = T0 + R5
R5 = T0 + R5 * C95T0 = (R3 - R7) * C92
R7 = (R1 - R7) * C93R3 = (R1 - R3) * C94
R1 = T9 + T0 + R3T0 = T9 - T0 - R7
R7 = T9 + R7 - R3R9 = (R2 - R4 + R8) * C31
R3 = (R4 + R8) * C96R8 = (R2 - R8) * C97
R4 = (R2 + R4) * C98R2 = R6 + R3 + R4
R3 = R6 - R8 - R3R8 = R6 + R8 - R4
CX(IP(2)) = T1 - R2
X(IP(9)) = T1 + R2Y(IP(2)) = R1 + T2
Y(IP(9)) = R1 - T2X(IP(3)) = T3 + R3
X(IP(8)) = T3 - R3Y(IP(3)) = T0 - T4
Y(IP(8)) = T0 + T4X(IP(4)) = T5 - R9
X(IP(7)) = T5 + R9Y(IP(4)) = R5 + T6
Y(IP(7)) = R5 - T6X(IP(5)) = T7 - R8
X(IP(6)) = T7 + R8Y(IP(5)) = R7 + T8
Y(IP(6)) = R7 - T8C
GOTO 20C-----------------WFTA N=16------------------------C
116 R1 = X(I(1)) + X(I(9))R2 = X(I(1)) - X(I(9))
R3 = X(I(2)) + X(I(10))R4 = X(I(2)) - X(I(10))
R5 = X(I(3)) + X(I(11))R6 = X(I(3)) - X(I(11))
R7 = X(I(4)) + X(I(12))R8 = X(I(4)) - X(I(12))
R9 = X(I(5)) + X(I(13))R10= X(I(5)) - X(I(13))
R11 = X(I(6)) + X(I(14))R12 = X(I(6)) - X(I(14))
R13 = X(I(7)) + X(I(15))R14 = X(I(7)) - X(I(15))
R15 = X(I(8)) + X(I(16))R16 = X(I(8)) - X(I(16))
T1 = R1 + R9T2 = R1 - R9T3 = R3 + R11
T4 = R3 - R11T5 = R5 + R13
T6 = R5 - R13T7 = R7 + R15
T8 = R7 - R15R1 = T1 + T5
R3 = T1 - T5R5 = T3 + T7
R7 = T3 - T7X(IP( 1)) = R1 + R5
X(IP( 9)) = R1 - R5T1 = C81 * (T4 + T8)
T5 = C81 * (T4 - T8)R9 = T2 + T5
R11= T2 - T5R13 = T6 + T1
R15 = T6 - T1T1 = R4 + R16
T2 = R4 - R16T3 = C81 * (R6 + R14)
T4 = C81 * (R6 - R14)T5 = R8 + R12
T6 = R8 - R12T7 = C162 * (T2 - T6)
T2 = C163 * T2 - T7T6 = C164 * T6 - T7
T7 = R2 + T4T8 = R2 - T4
R2 = T7 + T2R4 = T7 - T2
R6 = T8 + T6R8 = T8 - T6
T7 = C165 * (T1 + T5)T2 = T7 - C164 * T1
T4 = T7 - C163 * T5T6 = R10 + T3
T8 = R10 - T3R10 = T6 + T2
R12 = T6 - T2R14 = T8 + T4
R16 = T8 - T4R1 = Y(I(1)) + Y(I(9))
S2 = Y(I(1)) - Y(I(9))S3 = Y(I(2)) + Y(I(10))
S4 = Y(I(2)) - Y(I(10))R5 = Y(I(3)) + Y(I(11))
S6 = Y(I(3)) - Y(I(11))S7 = Y(I(4)) + Y(I(12))
S8 = Y(I(4)) - Y(I(12))S9 = Y(I(5)) + Y(I(13))
S10= Y(I(5)) - Y(I(13))S11 = Y(I(6)) + Y(I(14))
S12 = Y(I(6)) - Y(I(14))S13 = Y(I(7)) + Y(I(15))
S14 = Y(I(7)) - Y(I(15))S15 = Y(I(8)) + Y(I(16))
S16 = Y(I(8)) - Y(I(16))T1 = R1 + S9
T2 = R1 - S9T3 = S3 + S11
T4 = S3 - S11T5 = R5 + S13
T6 = R5 - S13T7 = S7 + S15
T8 = S7 - S15R1 = T1 + T5
S3 = T1 - T5R5 = T3 + T7
S7 = T3 - T7Y(IP( 1)) = R1 + R5
Y(IP( 9)) = R1 - R5X(IP( 5)) = R3 + S7
X(IP(13)) = R3 - S7Y(IP( 5)) = S3 - R7
Y(IP(13)) = S3 + R7T1 = C81 * (T4 + T8)
T5 = C81 * (T4 - T8)S9 = T2 + T5
S11= T2 - T5S13 = T6 + T1
S15 = T6 - T1T1 = S4 + S16
T2 = S4 - S16T3 = C81 * (S6 + S14)
T4 = C81 * (S6 - S14)T5 = S8 + S12
T6 = S8 - S12T7 = C162 * (T2 - T6)
T2 = C163 * T2 - T7T6 = C164 * T6 - T7
T7 = S2 + T4T8 = S2 - T4
S2 = T7 + T2S4 = T7 - T2
S6 = T8 + T6S8 = T8 - T6
T7 = C165 * (T1 + T5)T2 = T7 - C164 * T1
T4 = T7 - C163 * T5T6 = S10 + T3
T8 = S10 - T3S10 = T6 + T2
S12 = T6 - T2S14 = T8 + T4
S16 = T8 - T4X(IP( 2)) = R2 + S10
X(IP(16)) = R2 - S10Y(IP( 2)) = S2 - R10
Y(IP(16)) = S2 + R10X(IP( 3)) = R9 + S13
X(IP(15)) = R9 - S13Y(IP( 3)) = S9 - R13
Y(IP(15)) = S9 + R13X(IP( 4)) = R8 - S16
X(IP(14)) = R8 + S16Y(IP( 4)) = S8 + R16
Y(IP(14)) = S8 - R16X(IP( 6)) = R6 + S14
X(IP(12)) = R6 - S14Y(IP( 6)) = S6 - R14
Y(IP(12)) = S6 + R14X(IP( 7)) = R11 - S15
X(IP(11)) = R11 + S15Y(IP( 7)) = S11 + R15
Y(IP(11)) = S11 - R15X(IP( 8)) = R4 - S12
X(IP(10)) = R4 + S12Y(IP( 8)) = S4 + R12
Y(IP(10)) = S4 - R12C
GOTO 20C
20 CONTINUE10 CONTINUE
RETURNEND
Notification Switch
Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|