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
|
c
c-----------------------------------------------------------------------
c subroutine: ffa
c fast fourier analysis subroutine
c-----------------------------------------------------------------------
c
subroutine ffa(b, nfft)
c
c this subroutine replaces the real vector b(k), (k=1,2,...,n),
c with its finite discrete fourier transform. the dc term is
c returned in location b(1) with b(2) set to 0. thereafter, the
c jth harmonic is returned as a complex number stored as
c b(2*j+1) + i b(2*j+2). note that the n/2 harmonic is returned
c in b(n+1) with b(n+2) set to 0. hence, b must be dimensioned
c to size n+2.
c subroutine is called as ffa (b,n) where n=2**m and b is an
c n term real array. a real-valued, radix 8 algorithm is used
c with in-place reordering and the trig functions are computed as
c needed.
c
dimension b(2)
common /con/ pii, p7, p7two, c22, s22, pi2
c
c iw is a machine dependent write device number
c
iw = i1mach(2)
c
pii = 4.*atan(1.)
pi8 = pii/8.
p7 = 1./sqrt(2.)
p7two = 2.*p7
c22 = cos(pi8)
s22 = sin(pi8)
pi2 = 2.*pii
n = 1
do 10 i=1,15
m = i
n = n*2
if (n.eq.nfft) go to 20
10 continue
write (iw,9999)
9999 format (30h nfft not a power of 2 for ffa)
stop
20 continue
n8pow = m/3
c
c do a radix 2 or radix 4 iteration first if one is required
c
if (m-n8pow*3-1) 50, 40, 30
30 nn = 4
int = n/nn
call r4tr(int, b(1), b(int+1), b(2*int+1), b(3*int+1))
go to 60
40 nn = 2
int = n/nn
call r2tr(int, b(1), b(int+1))
go to 60
50 nn = 1
c
c perform radix 8 iterations
c
60 if (n8pow) 90, 90, 70
70 do 80 it=1,n8pow
nn = nn*8
int = n/nn
call r8tr(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1),
* b(4*int+1), b(5*int+1), b(6*int+1), b(7*int+1), b(1),
* b(int+1), b(2*int+1), b(3*int+1), b(4*int+1), b(5*int+1),
* b(6*int+1), b(7*int+1))
80 continue
c
c perform in-place reordering
c
90 call ord1(m, b)
call ord2(m, b)
t = b(2)
b(2) = 0.
b(nfft+1) = t
b(nfft+2) = 0.
do 100 i=4,nfft,2
b(i) = -b(i)
100 continue
return
end
|
