Program 13: In-Place, In-Order Prime Factor FFT Algorithm
m17389
Program 13: In-Place, In-Order Prime Factor FFT Algorithm
1.2
2008/05/28 13:03:23 GMT-5
2009/04/04 09:37:30.485 GMT-5
C.
Sidney
Burrus
C. Sidney Burrus
csb@rice.edu
C.
Sidney
Burrus
C. Sidney Burrus
csb@rice.edu
Daniel
Collins
Williamson
Daniel Williamson
dcwill@cnx.org
C.
Sidney
Burrus
C. Sidney Burrus
csb@rice.edu
FFT
PFA
Prime Factor Algorithm
Winograd
Mathematics and Statistics
Science and Technology
Prime Factor FFT Algorithm without unscrambler
en
In Place, In Order Prime Factor FFT Algorithm
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 both
in-place and in-order, so requires no unscrambler.
C
C A PRIME FACTOR FFT PROGRAM
C IN-PLACE AND IN-ORDER
C COMPLEX INPUT DATA IN ARRAYS X AND Y
C 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,16
C PROGRAM BY C. S. BURRUS, RICE UNIVERSITY
C SEPT 1983
C----------------------------------------------------
C
SUBROUTINE 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 /
C
C-----------------NESTED LOOPS----------------------------------
C
DO 10 K=1, M
N1 = NI(K)
N2 = N/N1
L = 1
N3 = N2 - N1*(N2/N1)
DO 15 J = 1, N1
LP(J) = L
L = L + N3
IF (L.GT.N1) L = L - N1
15 CONTINUE
C
DO 20 J=1, N, N1
IT = J
DO 30 L=1, N1
I(L) = IT
IP(LP(L)) = IT
IT = IT + N2
IF (IT.GT.N) IT = IT - N
30 CONTINUE
GOTO (20,102,103,104,105,20,107,108,109,
+ 20,20,20,20,20,20,116),N1
C----------------WFTA N=2--------------------------------
C
102 R1 = X(I(1))
X(I(1)) = R1 + X(I(2))
X(I(2)) = R1 - X(I(2))
C
R1 = Y(I(1))
Y(IP(1)) = R1 + Y(I(2))
Y(IP(2)) = R1 - Y(I(2))
C
GOTO 20
C
C----------------WFTA N=3--------------------------------
C
103 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 * C32
C
S2 = (Y(I(2)) - Y(I(3))) * C31
S1 = Y(I(2)) + Y(I(3))
Y(I(1))= Y(I(1)) + S1
S1 = Y(I(1)) + S1 * C32
C
X(IP(2)) = R1 - S2
X(IP(3)) = R1 + S2
Y(IP(2)) = S1 + R2
Y(IP(3)) = S1 - R2
C
GOTO 20
C
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 - R2
C
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 - R2
C
R1 = X(I(2)) - X(I(4))
R2 = Y(I(2)) - Y(I(4))
C
X(IP(2)) = T1 + R2
X(IP(4)) = T1 - R2
Y(IP(2)) = T2 - R1
Y(IP(4)) = T2 + R1
C
GOTO 20
C----------------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))
C
T = (R1 - R3) * C54
R1 = R1 + R3
X(I(1)) = X(I(1)) + R1
R1 = X(I(1)) + R1 * C55
C
R3 = R1 - T
R1 = R1 + T
C
T = (R4 + R2) * C51
R4 = T + R4 * C52
R2 = T + R2 * C53
C
S1 = 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) * C54
S1 = S1 + S3
Y(I(1)) = Y(I(1)) + S1
S1 = Y(I(1)) + S1 * C55
C
S3 = S1 - T
S1 = S1 + T
C
T = (S4 + S2) * C51
S4 = T + S4 * C52
S2 = T + S2 * C53
C
X(IP(2)) = R1 + S2
X(IP(5)) = R1 - S2
X(IP(3)) = R3 - S4
X(IP(4)) = R3 + S4
C
Y(IP(2)) = S1 - R2
Y(IP(5)) = S1 + R2
Y(IP(3)) = S3 + R4
Y(IP(4)) = S3 - R4
C
GOTO 20
C-----------------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))
C
T3 = (R1 - R2) * C74
T = (R1 - R3) * C72
R1 = R1 + R2 + R3
X(I(1)) = X(I(1)) + R1
R1 = X(I(1)) + R1 * C71
R2 =(R3 - R2) * C73
R3 = R1 - T + R2
R2 = R1 - R2 - T3
R1 = R1 + T + T3
T = (R6 - R5) * C78
T3 =(R6 + R4) * C76
R6 =(R6 + R5 - R4) * C75
R5 =(R5 + R4) * C77
R4 = R6 - T3 + R5
R5 = R6 - R5 - T
R6 = R6 + T3 + T
C
T3 = (S1 - S2) * C74
T = (S1 - S3) * C72
S1 = S1 + S2 + S3
Y(I(1)) = Y(I(1)) + S1
S1 = Y(I(1)) + S1 * C71
S2 =(S3 - S2) * C73
S3 = S1 - T + S2
S2 = S1 - S2 - T3
S1 = S1 + T + T3
T = (S6 - S5) * C78
T3 = (S6 + S4) * C76
S6 = (S6 + S5 - S4) * C75
S5 = (S5 + S4) * C77
S4 = S6 - T3 + S5
S5 = S6 - S5 - T
S6 = S6 + T3 + T
C
X(IP(2)) = R3 + S4
X(IP(7)) = R3 - S4
X(IP(3)) = R1 + S6
X(IP(6)) = R1 - S6
X(IP(4)) = R2 - S5
X(IP(5)) = R2 + S5
Y(IP(4)) = S2 + R5
Y(IP(5)) = S2 - R5
Y(IP(2)) = S3 - R4
Y(IP(7)) = S3 + R4
Y(IP(3)) = S1 - R6
Y(IP(6)) = S1 + R6
C
GOTO 20
C-----------------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 + R5
T2 = R1 - R5
T3 = R3 + R7
R3 =(R3 - R7) * C81
X(IP(1)) = T1 + T3
X(IP(5)) = T1 - T3
T1 = R2 + R3
T3 = R2 - R3
S1 = R4 - R8
R4 =(R4 + R8) * C81
S2 = R4 + R6
S3 = 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 + R5
R1 = R1 - R5
R5 = R3 + R7
R3 =(R3 - R7) * C81
Y(IP(1)) = T4 + R5
Y(IP(5)) = T4 - R5
R5 = R2 + R3
R2 = R2 - R3
R3 = R4 - R8
R4 =(R4 + R8) * C81
R7 = R4 + R6
R4 = R4 - R6
X(IP(2)) = T1 + R7
X(IP(8)) = T1 - R7
X(IP(3)) = T2 + R3
X(IP(7)) = T2 - R3
X(IP(4)) = T3 + R4
X(IP(6)) = T3 - R4
Y(IP(2)) = R5 - S2
Y(IP(8)) = R5 + S2
Y(IP(3)) = R1 - S1
Y(IP(7)) = R1 + S1
Y(IP(4)) = R2 - S3
Y(IP(6)) = R2 + S3
C
GOTO 20
C-----------------WFTA N=9-----------------------
C
109 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))) * C31
R7 = X(I(5)) + X(I(6))
R8 = X(I(5)) - X(I(6))
T0 = X(I(1)) + R5
T7 = X(I(1)) + R5 * C95
R5 = R1 + R3 + R7
X(I(1)) = T0 + R5
T5 = T0 + R5 * C95
T3 = (R3 - R7) * C92
R7 = (R1 - R7) * C93
R3 = (R1 - R3) * C94
T1 = T7 + T3 + R3
T3 = T7 - T3 - R7
T7 = T7 + R7 - R3
T6 = (R2 - R4 + R8) * C31
T4 = (R4 + R8) * C96
R8 = (R2 - R8) * C97
R2 = (R2 + R4) * C98
T2 = T8 + T4 + R2
T4 = T8 - T4 - R8
T8 = T8 + R8 - R2
C
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)) + R5
T9 = Y(I(1)) + R5 * C95
R5 = R1 + R3 + R7
Y(I(1)) = T0 + R5
R5 = T0 + R5 * C95
T0 = (R3 - R7) * C92
R7 = (R1 - R7) * C93
R3 = (R1 - R3) * C94
R1 = T9 + T0 + R3
T0 = T9 - T0 - R7
R7 = T9 + R7 - R3
R9 = (R2 - R4 + R8) * C31
R3 = (R4 + R8) * C96
R8 = (R2 - R8) * C97
R4 = (R2 + R4) * C98
R2 = R6 + R3 + R4
R3 = R6 - R8 - R3
R8 = R6 + R8 - R4
C
X(IP(2)) = T1 - R2
X(IP(9)) = T1 + R2
Y(IP(2)) = R1 + T2
Y(IP(9)) = R1 - T2
X(IP(3)) = T3 + R3
X(IP(8)) = T3 - R3
Y(IP(3)) = T0 - T4
Y(IP(8)) = T0 + T4
X(IP(4)) = T5 - R9
X(IP(7)) = T5 + R9
Y(IP(4)) = R5 + T6
Y(IP(7)) = R5 - T6
X(IP(5)) = T7 - R8
X(IP(6)) = T7 + R8
Y(IP(5)) = R7 + T8
Y(IP(6)) = R7 - T8
C
GOTO 20
C-----------------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 + R9
T2 = R1 - R9
T3 = R3 + R11
T4 = R3 - R11
T5 = R5 + R13
T6 = R5 - R13
T7 = R7 + R15
T8 = R7 - R15
R1 = T1 + T5
R3 = T1 - T5
R5 = T3 + T7
R7 = T3 - T7
X(IP( 1)) = R1 + R5
X(IP( 9)) = R1 - R5
T1 = C81 * (T4 + T8)
T5 = C81 * (T4 - T8)
R9 = T2 + T5
R11= T2 - T5
R13 = T6 + T1
R15 = T6 - T1
T1 = R4 + R16
T2 = R4 - R16
T3 = C81 * (R6 + R14)
T4 = C81 * (R6 - R14)
T5 = R8 + R12
T6 = R8 - R12
T7 = C162 * (T2 - T6)
T2 = C163 * T2 - T7
T6 = C164 * T6 - T7
T7 = R2 + T4
T8 = R2 - T4
R2 = T7 + T2
R4 = T7 - T2
R6 = T8 + T6
R8 = T8 - T6
T7 = C165 * (T1 + T5)
T2 = T7 - C164 * T1
T4 = T7 - C163 * T5
T6 = R10 + T3
T8 = R10 - T3
R10 = T6 + T2
R12 = T6 - T2
R14 = T8 + T4
R16 = T8 - T4
R1 = 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 - S9
T3 = S3 + S11
T4 = S3 - S11
T5 = R5 + S13
T6 = R5 - S13
T7 = S7 + S15
T8 = S7 - S15
R1 = T1 + T5
S3 = T1 - T5
R5 = T3 + T7
S7 = T3 - T7
Y(IP( 1)) = R1 + R5
Y(IP( 9)) = R1 - R5
X(IP( 5)) = R3 + S7
X(IP(13)) = R3 - S7
Y(IP( 5)) = S3 - R7
Y(IP(13)) = S3 + R7
T1 = C81 * (T4 + T8)
T5 = C81 * (T4 - T8)
S9 = T2 + T5
S11= T2 - T5
S13 = T6 + T1
S15 = T6 - T1
T1 = S4 + S16
T2 = S4 - S16
T3 = C81 * (S6 + S14)
T4 = C81 * (S6 - S14)
T5 = S8 + S12
T6 = S8 - S12
T7 = C162 * (T2 - T6)
T2 = C163 * T2 - T7
T6 = C164 * T6 - T7
T7 = S2 + T4
T8 = S2 - T4
S2 = T7 + T2
S4 = T7 - T2
S6 = T8 + T6
S8 = T8 - T6
T7 = C165 * (T1 + T5)
T2 = T7 - C164 * T1
T4 = T7 - C163 * T5
T6 = S10 + T3
T8 = S10 - T3
S10 = T6 + T2
S12 = T6 - T2
S14 = T8 + T4
S16 = T8 - T4
X(IP( 2)) = R2 + S10
X(IP(16)) = R2 - S10
Y(IP( 2)) = S2 - R10
Y(IP(16)) = S2 + R10
X(IP( 3)) = R9 + S13
X(IP(15)) = R9 - S13
Y(IP( 3)) = S9 - R13
Y(IP(15)) = S9 + R13
X(IP( 4)) = R8 - S16
X(IP(14)) = R8 + S16
Y(IP( 4)) = S8 + R16
Y(IP(14)) = S8 - R16
X(IP( 6)) = R6 + S14
X(IP(12)) = R6 - S14
Y(IP( 6)) = S6 - R14
Y(IP(12)) = S6 + R14
X(IP( 7)) = R11 - S15
X(IP(11)) = R11 + S15
Y(IP( 7)) = S11 + R15
Y(IP(11)) = S11 - R15
X(IP( 8)) = R4 - S12
X(IP(10)) = R4 + S12
Y(IP( 8)) = S4 + R12
Y(IP(10)) = S4 - R12
C
GOTO 20
C
20 CONTINUE
10 CONTINUE
RETURN
END