1C -------------------------------------------------------------------------
2C | File name: fftd.for Sean O'Connor 9-22-87 |
3C | 7-01-17 |
4C -------------------------------------------------------------------------
5
6C -------------------------------------------------------------------------
7C | VAX 11/785 VMS 4.5 VAX EXTENDED FORTRAN 77 |
8C | MacBook Pro Ubuntu LTS 16.04 Linux gfortran |
9C -------------------------------------------------------------------------
10
11
12C ****************************************************************************
13C * *
14C * FFTD *
15C * *
16C ****************************************************************************
17C
18C
19C DESCRIPTION
20C
21C This program tests the subroutine FFT in the file FFT.FOR It does a
22C forward discrete Fourier transform (DFT). It then does an inverse
23C transform of the forward transform to check the roundoff error.
24C
25C
26C INPUT
27C
28C The number of complex points in the transform.
29C Whether to do a forward (1) or inverse (-1) transform.
30C The data points (complex numbers).
31C
32C
33C OUTPUT
34C
35C The DFT of the input.
36C The inverse DFT of the first DFT.
37C
38C
39C EXAMPLE
40C
41C $ ./FFTD
42C
43C Enter number of complex points in transform: 4
44C
45C Enter 1 (forward) or -1 (inverse) transform: 1
46C
47C Enter point no. 0 (Real, Imag) > (1,1)
48C Enter point no. 1 (Real, Imag) > (2,2)
49C Enter point no. 2 (Real, Imag) > (3,3)
50C Enter point no. 3 (Real, Imag) > (4,4)
51C
52C num_points = 4
53C direction = 1
54C
55C Point No. 0 = ( 1.0000000000000 1.0000000000000 )
56C Point No. 1 = ( 2.0000000000000 2.0000000000000 )
57C Point No. 2 = ( 3.0000000000000 3.0000000000000 )
58C Point No. 3 = ( 4.0000000000000 4.0000000000000 )
59C
60C The transform of the input is . . .
61C
62C Point No. 0= ( 5.0000000000000 5.0000000000000 )
63C Point No. 1= ( -2.0000000000000 0.00000000000000E+000)
64C Point No. 2= ( -1.0000000000000 -1.0000000000000 )
65C Point No. 3= (-0.59604644775391E-007 -2.0000000000000 )
66C
67C The inverse transform of the transform above is . . .
68C
69C Point No. 0= ( 1.0000000000000 1.0000000000000 )
70C Point No. 1= ( 2.0000000000000 2.0000000000000 )
71C Point No. 2= ( 3.0000000000000 3.0000000000000 )
72C Point No. 3= ( 4.0000000000000 4.0000000000000 )
73C
74C
75C INSTALLATION
76C
77C To compile, type
78C
79C $ gfortran FFTD.FOR FORTRAN FFT.FOR -o FFT
80C
81C To run, type
82C
83C $ ./FFTD
84C
85C
86C ****************************************************************************
87C * *
88C ****************************************************************************
89
90 PROGRAM fftd
91
92C -------------------------------------------------------------------------
93C | MACRO DEFINITIONS |
94C -------------------------------------------------------------------------
95
96 PARAMETER ( MAX_NUM_POINTS = 255 ) ! Maximum number of points in FFT.
97
98
99C -------------------------------------------------------------------------
100C | LOCAL VARIABLES |
101C -------------------------------------------------------------------------
102
103 INTEGER i, ! Loop counter
104 $ num_points, ! Number of points in the transform
105 $ direction ! Direction of transform 1 for forward,
106 ! -1 for inverse
107
108 COMPLEX X( 0 : MAX_NUM_POINTS ) ! Array of complex data points
109 ! to be transformed.
110
111
112C -------------------------------------------------------------------------
113C | PROGRAM CODE |
114C -------------------------------------------------------------------------
115
116C
117C Read the number of points in the FFT
118C
119 WRITE( 6, 100 )
120100 FORMAT ( ' Enter number of complex points in transform: ' )
121
122 READ ( 5 , * ) num_points ! Free format input
123
124
125
126C
127C Decide whether to do the forward or inverse FFT
128C
129
130 WRITE( 6, 200 )
131200 FORMAT ( //, ' Enter 1 for a forward transform ', /,
132 $ ' or -1 for an inverse transform: ' )
133
134 READ ( 5, * ) direction ! Free format input
135
136
137C
138C Check to see that the number of points in the FFT
139C is valid. If not, halt.
140C
141
142 IF ( (num_points .LE. 0) .OR.
143 $ (num_points .GT. MAX_NUM_POINTS) ) THEN
144
145 WRITE( 6, 300 )
146300 FORMAT( ' The number of points is out of range ')
147 STOP
148
149 END IF
150
151
152
153C
154C Read in the complex data points, one by one
155C
156 DO i = 0, num_points - 1
157
158 WRITE( 6, 400 ) i
159400 FORMAT ( ' Enter point no. ', I3,' (Real, Imag) > ')
160 READ ( 5 , * ) X( i ) ! Free format input
161
162 END DO
163
164
165C
166C Print out the points just read in.
167C
168
169 WRITE( 6, 500 ) num_points
170500 FORMAT ( //, ' num_points = ', I10, / )
171
172 WRITE( 6, 600 ) direction
173600 FORMAT ( ' direction = ', I10, // )
174
175 DO i = 0, num_points - 1
176
177 WRITE( 6 , 700 ) i , X ( i )
178700 FORMAT ( ' Point no. ', I3, ' = ( ', 2G22.14E3, ') ' )
179
180 END DO
181
182
183C
184C Call the FFT subroutine.
185C
186 CALL FFT( num_points, X, direction )
187
188
189C
190C Print the transform results.
191C
192 WRITE( 6, 800 )
193800 FORMAT ( //, ' The transform of the input data is . . . ', // )
194
195
196 DO i = 0, num_points - 1
197
198 WRITE( 6, 900 ) i , X( i )
199900 FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' )
200
201 END DO
202
203
204C
205C Call the FFT subroutine with -direction to do the inverse transform.
206C
207 CALL FFT( num_points, X, -direction )
208
209
210
211C
212C Print the inverse transform coefficients.
213C
214 WRITE ( 6, 1000 )
2151000 FORMAT ( //,' The inverse transform ',
216 $ 'of the transform above is . . .', // )
217
218
219 DO i = 0, num_points - 1
220
221 WRITE ( 6, 1100 ) i , X( i )
2221100 FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' )
223
224 END DO
225
226
227
228C
229C Test the error checking in the FFT routine.
230C
231
232C Number of points is too large. Currently, MAX_TABLE_SIZE in
233C subroutine FFT is 20, so we can transform up to 2 ** 21 points.
234
235 PRINT *, ' Calling FFT with too many points.'
236
237 CALL FFT( 2 ** 22, X, 1 )
238
239
240C Number of points is less than 2.
241
242 PRINT *, ' Calling FFT with 1 point.'
243
244 CALL FFT( 1, X, 1 )
245
246
247C Calling FFT with 2^ 15 + 1 points (not a power of 2)
248
249 PRINT *, ' Calling FFT with 2 ^ 15 + 1 points'
250
251 CALL FFT( 2 * 15 + 1, X, 1 )
252
253
254C Direction is not 1 or -1.
255
256 PRINT *, ' Calling FFT with direction = 2'
257
258 CALL FFT( 16, X, 2)
259
260
261 STOP
262 END ! of program fftd
