-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfst.f
154 lines (126 loc) · 3.95 KB
/
fst.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
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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
SUBROUTINE FST (VECTOR, N, Z, NN)
INTEGER N, NN, INV
LOGICAL INVERSE
REAL*8 VECTOR(1:N-1)
COMPLEX*16 Z(0:NN-1)
C
INVERSE = .TRUE.
CALL FAST_SINE (VECTOR, N, INVERSE, Z, NN)
C
RETURN
END
C
SUBROUTINE FSTINV (VECTOR, N, Z, NN)
INTEGER N, NN, INV
LOGICAL INVERSE
REAL*8 VECTOR(1:N-1)
COMPLEX*16 Z(0:NN-1)
C
INVERSE = .FALSE.
CALL FAST_SINE (VECTOR, N, INVERSE, Z, NN)
C
RETURN
END
C
SUBROUTINE FAST_SINE (VECTOR, N, INVERSE, Z, NN)
INTEGER N, K, KK
LOGICAL INVERSE
REAL*8 VECTOR(1:N-1)
COMPLEX*16 Z(0:nn-1)
C**********************************************************
C IF WE EXTEND THE INPUT VECTOR TO THE 2*N VECTOR
C XE = [0,X(1),...,X(N-1),0,-X(N-1),-X(N-2),...,-X(1)]
C THEN
C SINE TRANSFORM = [ i/2 FFT(XE) ](1:N-1)
C**********************************************************
KK = 0
Z(KK) = DCMPLX(0D0,0D0)
KK = KK+1
DO K = 1,N-1
Z(KK) = DCMPLX(VECTOR(K),0D0)
KK = KK+1
END DO
Z(KK) = DCMPLX(0D0,0D0)
KK = KK+1
DO K = N-1,1,-1
Z(KK) = DCMPLX((-1)*VECTOR(K),0D0)
KK = KK+1
END DO
C*************************************************************
C APPLY FAST FOURIER TRANSFORM AND EXTRACT SINE TRANSFORM
C*************************************************************
CALL FFT(Z,2*N)
do i=0,(2*n)-1
Z(i) = Z(i) * (0d0,5d-1)
enddo
do i=1,n-1
VECTOR(i) = Z(i)
enddo
C****************************************************
C INVERSE SINE TRANSFORM = (2/N)SINE TRANSFORM
C****************************************************
IF (INVERSE) THEN
do i=1,n-1
VECTOR(i) = (2 * VECTOR(i)) / N
enddo
END IF
RETURN
END
C================================================================
C COMPLEX FAST FOURIER TRANSFORM
C================================================================
SUBROUTINE FFT(X,M)
COMPLEX*16 X(0:M)
INTEGER M
C BIT REVERSE VARIABLES
INTEGER INDEX, RINDEX, Q, S, LOGCOUNTER
COMPLEX*16 TEMP
C UPDATE VARIABLES
INTEGER LEVEL, COLS, HALFLEVEL, J, K
COMPLEX*16 OMEGA, CTEMP
REAL*8 PI
PI = 3.141592653589793238462643d0
C*********************************************
C BIT REVERSE PERMUTATION OF THE VECTOR
C*********************************************
DO INDEX = 0,M-1
C BIT REVERSE THE INDEX
LOGCOUNTER = M
RINDEX = 0
Q = INDEX
DO while (LOGCOUNTER .gt. 1)
s = Q/2
RINDEX = (2*RINDEX)+ (Q - (2*S))
Q = S
LOGCOUNTER = LOGCOUNTER/2
END DO
C SWAP ENTRIES
IF (RINDEX .gt. INDEX) THEN
TEMP = X(INDEX)
X(INDEX) = X(RINDEX)
X(RINDEX) = TEMP
END IF
END DO
C******************************
C PERFORM BUTTERFLY UPDATE
C******************************
HALFLEVEL = 1
LEVEL = 2
DO while (LEVEL .le. M)
COLS = M/LEVEL
C****************************************************************
C APPLY BUTTERFLY UPDATE TO KTH ROW OF X(LEVEL:LEVEL_COLUMNS)
C****************************************************************
DO K = 0,COLS-1
DO J = 0,HALFLEVEL-1
OMEGA = DCMPLX(COS(2*PI*J/LEVEL),(-1)*SIN(2*PI*J/LEVEL))
CTEMP = OMEGA * X((K*LEVEL)+ J +HALFLEVEL)
X((K*LEVEL)+ J +HALFLEVEL) = X((K*LEVEL)+ J) - CTEMP
X((K*LEVEL)+ J ) = X((K*LEVEL)+ J) + CTEMP
END DO
END DO
HALFLEVEL = LEVEL
LEVEL = LEVEL * 2
END DO
return
end