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
|
C
SUBROUTINE FFTRC (A,N,X,IWK,WK)
C SPECIFICATIONS FOR ARGUMENTS
INTEGER N,IWK(1)
REAL*8 A(N),WK(1)
COMPLEX*16 X(1)
C SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER ND2P1,ND2,I,MTWO,M,IMAX,ND4,NP2,K,NMK,J
REAL*8 RPI,ZERO,ONE,HALF,THETA,TP,G(2),B(2),Z(2),AI,
1 AR
COMPLEX*16 XIMAG,ALPH,BETA,GAM,S1,ZD
EQUIVALENCE (GAM,G(1)),(ALPH,B(1)),(Z(1),AR),(Z(2),AI),
1 (ZD,Z(1))
DATA ZERO/0.0D0/,HALF/0.5D0/,ONE/1.0D0/,IMAX/24/
DATA RPI/3.141592653589793D0/
C FIRST EXECUTABLE STATEMENT
IF (N .NE. 2) GO TO 5
C N EQUAL TO 2
ZD = DCMPLX(A(1),A(2))
THETA = AR
TP = AI
X(2) = DCMPLX(THETA-TP,ZERO)
X(1) = DCMPLX(THETA+TP,ZERO)
GO TO 9005
5 CONTINUE
C N GREATER THAN 2
ND2 = N/2
ND2P1 = ND2+1
C MOVE A TO X
J = 1
DO 6 I=1,ND2
X(I) = DCMPLX(A(J),A(J+1))
J = J+2
6 CONTINUE
C COMPUTE THE CENTER COEFFICIENT
GAM = DCMPLX(ZERO,ZERO)
DO 10 I=1,ND2
GAM = GAM + X(I)
10 CONTINUE
TP = G(1)-G(2)
GAM = DCMPLX(TP,ZERO)
C DETERMINE THE SMALLEST M SUCH THAT
C N IS LESS THAN OR EQUAL TO 2**M
MTWO = 2
M = 1
DO 15 I=1,IMAX
IF (ND2 .LE. MTWO) GO TO 20
MTWO = MTWO+MTWO
M = M+1
15 CONTINUE
20 IF (ND2 .EQ. MTWO) GO TO 25
C N IS NOT A POWER OF TWO, CALL FFTCC
CALL FFTCC (X,ND2,IWK,WK)
GO TO 30
C N IS A POWER OF TWO, CALL FFT2C
25 CALL FFT2C (X,M,IWK)
30 ALPH = X(1)
X(1) = B(1) + B(2)
ND4 = (ND2+1)/2
IF (ND4 .LT. 2) GO TO 40
NP2 = ND2 + 2
THETA = RPI/ND2
TP = THETA
XIMAG = DCMPLX(ZERO,ONE)
C DECOMPOSE THE COMPLEX VECTOR X
C INTO THE COMPONENTS OF THE TRANSFORM
C OF THE INPUT DATA.
DO 35 K = 2,ND4
NMK = NP2 - K
S1 = DCONJG(X(NMK))
ALPH = X(K) + S1
BETA = XIMAG*(S1-X(K))
S1 = DCMPLX(DCOS(THETA),DSIN(THETA))
X(K) = (ALPH+BETA*S1)*HALF
X(NMK) = DCONJG(ALPH-BETA*S1)*HALF
THETA = THETA + TP
35 CONTINUE
40 CONTINUE
X(ND2P1) = GAM
9005 RETURN
END
|
