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