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
|
c
c-----------------------------------------------------------------------
c subroutine: fourea
c performs cooley-tukey fast fourier transform
c-----------------------------------------------------------------------
c
subroutine fourea(data, n, isi)
c
c the cooley-tukey fast fourier transform in ansi fortran
c
c data is a one-dimensional complex array whose length, n, is a
c power of two. isi is +1 for an inverse transform and -1 for a
c forward transform. transform values are returned in the input
c array, replacing the input.
c transform(j)=sum(data(i)*w**((i-1)*(j-1))), where i and j run
c from 1 to n and w = exp (isi*2*pi*sqrt(-1)/n). program also
c computes inverse transform, for which the defining expression
c is invtr(j)=(1/n)*sum(data(i)*w**((i-1)*(j-1))).
c running time is proportional to n*log2(n), rather than to the
c classical n**2.
c after program by brenner, june 1967. this is a very short version
c of the fft and is intended mainly for demonstration. programs
c are available in this collection which run faster and are not
c restricted to powers of 2 or to one-dimensional arrays.
c see -- ieee trans audio (june 1967), special issue on fft.
c
complex data(1)
complex temp, w
ioutd = i1mach(2)
c
c check for power of two up to 15
c
nn = 1
do 10 i=1,15
m = i
nn = nn*2
if (nn.eq.n) go to 20
10 continue
write (ioutd,9999)
9999 format (30h n not a power of 2 for fourea)
stop
20 continue
c
pi = 4.*atan(1.)
fn = n
c
c this section puts data in bit-reversed order
c
j = 1
do 80 i=1,n
c
c at this point, i and j are a bit reversed pair (except for the
c displacement of +1)
c
if (i-j) 30, 40, 40
c
c exchange data(i) with data(j) if i.lt.j.
c
30 temp = data(j)
data(j) = data(i)
data(i) = temp
c
c implement j=j+1, bit-reversed counter
c
40 m = n/2
50 if (j-m) 70, 70, 60
60 j = j - m
m = (m+1)/2
go to 50
70 j = j + m
80 continue
c
c now compute the butterflies
c
mmax = 1
90 if (mmax-n) 100, 130, 130
100 istep = 2*mmax
do 120 m=1,mmax
theta = pi*float(isi*(m-1))/float(mmax)
w = cmplx(cos(theta),sin(theta))
do 110 i=m,n,istep
j = i + mmax
temp = w*data(j)
data(j) = data(i) - temp
data(i) = data(i) + temp
110 continue
120 continue
mmax = istep
go to 90
130 if (isi) 160, 140, 140
c
c for inv trans -- isi=1 -- multiply output by 1/n
c
140 do 150 i=1,n
data(i) = data(i)/fn
150 continue
160 return
end
|
