Program 13: In-Place, In-Order Prime Factor FFT Algorithm
Summary: Prime Factor FFT Algorithm without unscrambler
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
References
Comments, questions, feedback, criticisms?
Send feedback
- E-mail the author of the module, Program 13: In-Place, In-Order Prime Factor FFT Algorithm
Footer
More about this content: Metadata | Version History
This work is licensed by C. Sidney Burrus under a Creative Commons Attribution License (CC-BY 2.0).
Last edited by C. Sidney Burrus on Sep 2, 2008 2:13 pm GMT-5.