This page is optimized for mobile devices, if you would prefer the desktop version just
click here
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 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
Read also:
OpenStax, Fast fourier transforms. OpenStax CNX. Nov 18, 2012 Download for free at http://cnx.org/content/col10550/1.22
Google Play and the Google Play logo are trademarks of Google Inc.