-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREALFT.f
53 lines (53 loc) · 1.32 KB
/
REALFT.f
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
SUBROUTINE realft(data,n,isign)
INTEGER isign,n
REAL data(n)
CU USES four1
INTEGER i,i1,i2,i3,i4,n2p3
REAL c1,c2,h1i,h1r,h2i,h2r,wis,wrs
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
theta=3.141592653589793d0/dble(n/2)
c1=0.5
if (isign.eq.1) then
c2=-0.5
call four1(data,n/2,+1)
else
c2=0.5
theta=-theta
endif
wpr=-2.0d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.0d0+wpr
wi=wpi
n2p3=n+3
do 11 i=2,n/4
i1=2*i-1
i2=i1+1
i3=n2p3-i2
i4=i3+1
wrs=sngl(wr)
wis=sngl(wi)
h1r=c1*(data(i1)+data(i3))
h1i=c1*(data(i2)-data(i4))
h2r=-c2*(data(i2)+data(i4))
h2i=c2*(data(i1)-data(i3))
data(i1)=h1r+wrs*h2r-wis*h2i
data(i2)=h1i+wrs*h2i+wis*h2r
data(i3)=h1r-wrs*h2r+wis*h2i
data(i4)=-h1i+wrs*h2i+wis*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
11 continue
if (isign.eq.1) then
h1r=data(1)
data(1)=h1r+data(2)
data(2)=h1r-data(2)
else
h1r=data(1)
data(1)=c1*(h1r+data(2))
data(2)=c1*(h1r-data(2))
call four1(data,n/2,-1)
endif
return
END
C (C) Copr. 1986-92 Numerical Recipes Software *!^3#!0Y..