1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
c
c-----------------------------------------------------------------------
c subroutine: fft842
c fast fourier transform for n=2**m
c complex input
c-----------------------------------------------------------------------
c
subroutine fft842(in, n, x, y)
c
c this program replaces the vector z=x+iy by its finite
c discrete, complex fourier transform if in=0. the inverse transform
c is calculated for in=1. it performs as many base
c 8 iterations as possible and then finishes with a base 4 iteration
c or a base 2 iteration if needed.
c
c the subroutine is called as subroutine fft842 (in,n,x,y).
c the integer n (a power of 2), the n real location array x, and
c the n real location array y must be supplied to the subroutine.
c
dimension x(2), y(2), l(15)
common /con2/ pi2, p7
equivalence (l15,l(1)), (l14,l(2)), (l13,l(3)), (l12,l(4)),
* (l11,l(5)), (l10,l(6)), (l9,l(7)), (l8,l(8)), (l7,l(9)),
* (l6,l(10)), (l5,l(11)), (l4,l(12)), (l3,l(13)), (l2,l(14)),
* (l1,l(15))
c
c
c iw is a machine dependent write device number
c
iw = i1mach(2)
c
pi2 = 8.*atan(1.)
p7 = 1./sqrt(2.)
do 10 i=1,15
m = i
nt = 2**i
if (n.eq.nt) go to 20
10 continue
write (iw,9999)
9999 format (35h n is not a power of two for fft842)
stop
20 n2pow = m
nthpo = n
fn = nthpo
if (in.eq.1) go to 40
do 30 i=1,nthpo
y(i) = -y(i)
30 continue
40 n8pow = n2pow/3
if (n8pow.eq.0) go to 60
c
c radix 8 passes,if any.
c
do 50 ipass=1,n8pow
nxtlt = 2**(n2pow-3*ipass)
lengt = 8*nxtlt
call r8tx(nxtlt, nthpo, lengt, x(1), x(nxtlt+1), x(2*nxtlt+1),
* x(3*nxtlt+1), x(4*nxtlt+1), x(5*nxtlt+1), x(6*nxtlt+1),
* x(7*nxtlt+1), y(1), y(nxtlt+1), y(2*nxtlt+1), y(3*nxtlt+1),
* y(4*nxtlt+1), y(5*nxtlt+1), y(6*nxtlt+1), y(7*nxtlt+1))
50 continue
c
c is there a four factor left
c
60 if (n2pow-3*n8pow-1) 90, 70, 80
c
c go through the base 2 iteration
c
c
70 call r2tx(nthpo, x(1), x(2), y(1), y(2))
go to 90
c
c go through the base 4 iteration
c
80 call r4tx(nthpo, x(1), x(2), x(3), x(4), y(1), y(2), y(3), y(4))
c
90 do 110 j=1,15
l(j) = 1
if (j-n2pow) 100, 100, 110
100 l(j) = 2**(n2pow+1-j)
110 continue
ij = 1
do 130 j1=1,l1
do 130 j2=j1,l2,l1
do 130 j3=j2,l3,l2
do 130 j4=j3,l4,l3
do 130 j5=j4,l5,l4
do 130 j6=j5,l6,l5
do 130 j7=j6,l7,l6
do 130 j8=j7,l8,l7
do 130 j9=j8,l9,l8
do 130 j10=j9,l10,l9
do 130 j11=j10,l11,l10
do 130 j12=j11,l12,l11
do 130 j13=j12,l13,l12
do 130 j14=j13,l14,l13
do 130 ji=j14,l15,l14
if (ij-ji) 120, 130, 130
120 r = x(ij)
x(ij) = x(ji)
x(ji) = r
fi = y(ij)
y(ij) = y(ji)
y(ji) = fi
130 ij = ij + 1
if (in.eq.1) go to 150
do 140 i=1,nthpo
y(i) = -y(i)
140 continue
go to 170
150 do 160 i=1,nthpo
x(i) = x(i)/fn
y(i) = y(i)/fn
160 continue
170 return
end
|
