1 C ------------------------------------------------------------------------- 2 C | File name: fft.for Sean O'Connor 9-23-87 | 3 C ------------------------------------------------------------------------- 4 5 C ------------------------------------------------------------------------- 6 C | VAX 11/785 VMS 4.5 VAX EXTENDED FORTRAN 77 | 7 C ------------------------------------------------------------------------- 8 9 C 10 C **************************************************************************** 11 C * * 12 C * FFT * 13 C * * 14 C **************************************************************************** 15 C 16 C 17 C DESCRIPTION 18 C 19 C This program does an in-place one dimensional complex discrete Fourier 20 C transform. It uses an FFT algorithm for a number of input points which 21 C is a power of two. 22 C 23 C 24 C INPUT 25 C 26 C N, The number of complex points to be transformed. 27 C N is a power of 2. N must be at least 2; it 28 C must not exceed the dimension of array X, or the 29 C maximum size of the sine and cosine table. (See 30 C below). 31 C 32 C X, The array holding the complex data X( k ) 33 C for k = 0, 1, . . ., N - 1 which is to be 34 C transformed. It is dimensioned in the calling 35 C program. Only the array address is passed to 36 C this subroutine. 37 C 38 C D, The direction of the FFT. D is 1 for a 39 C forward transform and -1 for an inverse 40 C transform. 41 C 42 C 43 C OUTPUT 44 C 45 C X, On output, X holds the complex Fourier coefficients 46 C ^ 47 C X( k ) for k = 0, 1, . . ., N - 1. 48 C 49 C Warning messages to the terminal under these conditions: 50 C 51 C (1) D is not 1 or -1 52 C (2) N is less than 2. 53 C (3) N is greater than the sine table can handle. 54 C (4) N is not a power of 2. 55 C 56 C 57 C CALLING SEQUENCE 58 C 59 C CALL FFT ( N, X, D ) 60 C 61 C 62 C EXAMPLE 63 C 64 C DIMENSION X( 0 : 3 ) ! 4 point FFT 65 C 66 C N = 4 ! 4 point FFT 67 C D = 1 ! Forward transform 68 C 69 C CALL FFT( N, X, D ) ! Do the FFT 70 C 71 C D = -1 ! Do an inverse FFT 72 C 73 C CALL FFT( N, X, D ) 74 C 75 C Suppose X contains the data 76 C 77 C X ( 0 ) = 1 + i 78 C X ( 1 ) = 2 + 2 i 79 C X ( 2 ) = 3 + 3 i 80 C X ( 3 ) = 4 + 4 i 81 C 82 C then the forward FFT returns 83 C 84 C X ( 0 ) = 5 + 5 i 85 C X ( 1 ) = -2 86 C X ( 2 ) = -1 - i 87 C X ( 3 ) = - 2 i 88 C 89 C and the inverse FFT returns the original data. 90 C 91 C 92 C METHOD 93 C ^ 94 C The Fourier coefficients X( k ) are defined for k = 0, . . ., N - 1 95 C by the formula 96 C 97 C N - 1 2 Pi i 98 C ^ 1 ____ - ------ j k 99 C X( k ) = --------- \ N 100 C ----- / X( j ) e 101 C \ / N ---- 102 C v j = 0 103 C 104 C where e is 2.1718 . . ., Pi = 3.14159. . . and i is the square root 105 C of -1. 106 C 107 C This is the formula for the forward transform. The inverse transform 108 C is the same, except that the minus sign in the exponent above is a 109 C plus sign instead. 110 C --- 111 C The normalizing factor 1 / \/ N in front was picked so that the 112 C forward and inverse transforms look the same. 113 C 114 C 115 C The particular type of FFT algorithm we use is a 116 C 117 C complex, 118 C two to the N, 119 C decimation in frequency, 120 C input in order, output bit reversed 121 C 122 C type of FFT. It is described in the book 123 C 124 C THE FAST FOURIER TRANSFORM, 125 C F. Oran Brigham, Prentice Hall, 1974 126 C 127 C 128 C BUGS 129 C 130 C The bit-reversal routine could be made faster by using the 131 C VAX FORTRAN 77 library routines for testing, setting, clearing and 132 C shifting bits: BTEST, IBSET, IBCLR, ISHFT. 133 C 134 C You might want to use double precision complex numbers to get more 135 C accuracy. 136 C 137 C **************************************************************************** 138 C * * 139 C **************************************************************************** 140 141 142 SUBROUTINE FFT( N, X, D ) 143 144 IMPLICIT NONE ! Force declaration of all variables. 145 146 147 C ------------------------------------------------------------------------- 148 C | MACRO DEFINITIONS | 149 C ------------------------------------------------------------------------- 150 151 INTEGER MAX_TABLE_SIZE 152 REAL PI 153 154 PARAMETER ( PI = 3.1415926535897932385 ) 155 156 C We must set the size of the sine and cosine table used in the FFT. 157 C The table is called starting_multiplier. 158 C 159 C Suppose MAX is the largest FFT we would ever want 160 C to perform. Then 161 C 162 C MAX_TABLE_SIZE = LOG ( MAX ) - 1 163 C 2 164 C 165 C is the size of the table. This table can be used by the 166 C FFT for any number of points from 2 up to MAX. 167 C 168 C For example, if MAX_TABLE_SIZE = 14, then we can transform 169 C anywhere from 2 to 2 ** 15 = 32,768 points, using the same 170 C sine and cosine table. 171 C 172 173 PARAMETER ( MAX_TABLE_SIZE = 20 ) 174 175 176 C ------------------------------------------------------------------------- 177 C | CALLING ARGUMENTS | 178 C ------------------------------------------------------------------------- 179 180 INTEGER N ! Number of points in transform 181 182 COMPLEX X( 0 : * ) ! Complex data to be transformed 183 184 C Array X is dimensioned in the calling program. Only its starting 185 C address is passed into this subroutine. 186 187 INTEGER D ! Transform direction: 1 for a forward 188 ! transform, -1 for an inverse transform. 189 190 191 C ------------------------------------------------------------------------- 192 C | LOCAL VARIABLES | 193 C ------------------------------------------------------------------------- 194 195 C A column refers to a column of the FFT signal flow graph. A group is 196 C a collection of butterflies which have the same multiplier W^p. A 197 C butterfly is the basic computational unit of the FFT. 198 199 INTEGER log_of_n, ! Log base 2 of the number of points, N. 200 $ column, ! Current column. 201 $ number_of_groups, ! Number of groups in a column. 202 $ group, ! Current group number. 203 $ start_of_group, ! Array index of first butterfly in group. 204 $ end_of_group, ! Array index of last butterfly in group. 205 $ butterfly_height, ! Height of a butterfly. 206 $ offset, ! Lower branch of a butterfly. 207 $ i, ! Loop counter. Upper branch of butterfly. 208 $ rev_i, ! Bit-reversed array index i. 209 $ REVERSE_BITS ! Bit reversal function. 210 211 212 COMPLEX product, ! Temporary storage in a butterfly calculation. 213 $ temp, ! Same as above. 214 $ increment, ! Difference between current multiplier and next. 215 $ multiplier ! The multiplier W^p for a group of butterflies. 216 217 COMPLEX starting_multiplier( 0 : MAX_TABLE_SIZE ) ! Table of sines 218 ! and cosines. 219 220 REAL normalizer, ! 1/\/N normalizing factor. 221 $ angle ! Current angle in sine and cosine table. 222 223 224 LOGICAL called_already ! .TRUE. if this is the second call or 225 ! later call to this subroutine. 226 227 228 C ------------------------------------------------------------------------- 229 C | SUBROUTINE CODE | 230 C ------------------------------------------------------------------------- 231 232 C 233 C Initialize the variable called_already to .FALSE. at compile time. 234 C Change it to .TRUE. after the first call to this subroutine. 235 C Save the value so it does not get destroyed after returning from 236 C this subroutine. 237 C 238 SAVE called_already 239 DATA called_already /.FALSE./ 240 241 242 C 243 C Compute log base 2 of N. Perturb slightly to avoid roundoff error. 244 C For example, INT( 4.99999999999 ) = 4 (wrong!), but 245 C INT( 4.999999999 + .01 ) = 5 (correct!) 246 C 247 log_of_n = INT( ALOG( FLOAT( N ) ) / ( ALOG( 2.0 ) ) + .01 ) 248 249 250 251 C 252 C Error checking. Print a message to the terminal if any one of these 253 C conditions is true: 254 C 255 C (1) D is not 1 or -1. 256 C (2) N < 2. 257 C (3) log2( N ) > MAX_TABLE_SIZE + 1. 258 C (4) N is not a power of 2. 259 C 260 261 IF ((D .NE. 1) .AND. (D .NE. -1)) THEN 262 263 PRINT *, ' ERROR --- D is not 1 or -1' 264 RETURN 265 266 ELSE IF (N .LT. 2) THEN 267 268 PRINT *, ' ERROR --- N is less than 2. ' 269 RETURN 270 271 ELSE IF ( (log_of_n - 1) .GT. MAX_TABLE_SIZE) THEN 272 273 PRINT *, ' ERROR --- N is larger than the sine table can' 274 PRINT *, ' handle. Increase the parameter MAX_TABLE_SIZE' 275 PRINT *, ' in this subroutine and recompile.' 276 277 RETURN 278 279 ELSE IF ( (2 ** log_of_n) .NE. N ) THEN 280 281 PRINT *, ' ERROR --- N is not a power of 2.' 282 RETURN 283 284 END IF 285 286 287 C 288 C If called_already is false, generate the sine and cosine table for 289 C the first time, then set called_already to .TRUE. Now later calls to 290 C the subroutine don't recreate the sine and cosine table. 291 C 292 293 IF (.NOT. called_already) THEN 294 295 DO i = 0, MAX_TABLE_SIZE - 1 296 297 angle = - PI / ( 2 ** i ) 298 starting_multiplier( i ) = CEXP( CMPLX(0.0 , angle) ) 299 300 END DO 301 302 called_already = .TRUE. 303 304 END IF 305 306 307 308 C 309 C Compute the FFT according to its signal flow graph with 310 C computations arranged into butterflies, groups of butterflies 311 C with the same multiplier, and columns. 312 C 313 butterfly_height = N / 2 314 number_of_groups = 1 315 316 317 DO column = 1, log_of_n 318 319 multiplier = starting_multiplier( column - 1 ) 320 321 C 322 C Modify starting multiplier W^p when computing the inverse transform. 323 C 324 IF ( D .EQ. -1 ) THEN 325 326 multiplier = CONJG( multiplier ) 327 328 END IF 329 330 increment = multiplier 331 332 333 C 334 C Compute the first group. 335 C 336 DO i = 0, butterfly_height - 1 337 338 offset = i + butterfly_height 339 340 temp = X( i ) + X( offset ) 341 X( offset ) = X( i ) - X( offset ) 342 X( i ) = temp 343 344 END DO ! i 345 346 347 C 348 C Now do the remaining groups. 349 C 350 DO group = 2, number_of_groups 351 352 start_of_group = REVERSE_BITS( group - 1, log_of_n ) 353 end_of_group = start_of_group + 354 $ butterfly_height - 1 355 356 357 C Compute all butterflies in a group. 358 359 DO i = start_of_group, end_of_group 360 361 offset = i + butterfly_height 362 product = x( offset ) * multiplier 363 364 temp = X( i ) + product 365 X( offset ) = X( i ) - product 366 X( i ) = temp 367 368 END DO ! i 369 370 multiplier = multiplier * increment 371 372 END DO ! group 373 374 butterfly_height = butterfly_height / 2 375 number_of_groups = number_of_groups * 2 376 377 END DO ! column 378 379 380 C 381 C Normalize the results and permute them into order, two at a time. 382 C 383 normalizer = 1 / SQRT( FLOAT( N ) ) 384 385 386 DO i = 0, N - 1 387 388 rev_i = REVERSE_BITS( i, log_of_n ) 389 390 IF ( rev_i .GT. i ) THEN ! exchange not yet done. 391 392 temp = X( i ) 393 X( i ) = X( rev_i ) * normalizer 394 X( rev_i ) = temp * normalizer 395 396 ELSE IF ( rev_i .EQ. i ) THEN ! exchange with itself. 397 398 X( i ) = X( i ) * normalizer 399 400 END IF 401 402 END DO ! i 403 404 405 RETURN 406 END ! subroutine FFT 407 408 409 410 C **************************************************************************** 411 C * * 412 C * REVERSE_BITS * 413 C * * 414 C **************************************************************************** 415 C 416 C DESCRIPTION 417 C 418 C This function reverses the order of the bits in a number. 419 C 420 C 421 C INPUT 422 C 423 C n The number whose bits are to be reversed. 424 C 425 C num_bits The number of bits in N, counting leading zeros. 426 C 427 C 428 C OUTPUT 429 C 430 C REVERSE_BITS The number N with all its num_bits bits reversed. 431 C 432 C 433 C EXAMPLE 434 C 435 C REVERSE_BITS( (10) , 2 ) = (01) = 1, but 436 C 2 2 437 C 438 C REVERSE_BITS( (10) , 3 ) = REVERSE_BITS( (010) , 3 ) = (010) = (10) 439 C 2 2 2 2 440 C 441 C = 10 442 C 443 C 444 C METHOD 445 C 446 C The algorithm works by a method similar to base conversion. 447 C As the low order bits are broken off, one by one, they are 448 C multiplied by 2 and accumulated to generate the bit-reversed 449 C number. 450 C 451 C 452 C BUGS 453 C 454 C The bit-reversal routine could be made faster by using the 455 C VAX FORTRAN 77 library routines for testing, setting, clearing and 456 C shifting bits: BTEST, IBSET, IBCLR, ISHFT. 457 C 458 C 459 C **************************************************************************** 460 C * * 461 C **************************************************************************** 462 463 INTEGER FUNCTION REVERSE_BITS( n, num_bits ) 464 465 IMPLICIT NONE 466 467 468 C ------------------------------------------------------------------------- 469 C | CALLING ARGUMENTS | 470 C ------------------------------------------------------------------------- 471 472 INTEGER n, num_bits 473 474 475 C ------------------------------------------------------------------------- 476 C | LOCAL VARIABLES | 477 C ------------------------------------------------------------------------- 478 479 INTEGER bit_num, ! Current bit number. 480 $ remainder, ! Bits remaining to be reversed. 481 $ reversed_n, ! Bit-reversed value of n. 482 $ bit, ! The current least significant bit. 483 $ MOD, INT ! Function return types. 484 485 486 C ------------------------------------------------------------------------- 487 C | SUBROUTINE CODE | 488 C ------------------------------------------------------------------------- 489 490 remainder = n ! Load n. 491 reversed_n = 0 ! Clear buffer for bit-reversed n. 492 493 DO bit_num = 1, num_bits 494 495 bit = MOD( remainder, 2 ) ! Next least significant bit. 496 reversed_n = bit + 2 * reversed_n ! Shift left and add to buffer. 497 remainder = INT( remainder / 2 ) ! Remaining bits. 498 499 END DO 500 501 REVERSE_BITS = reversed_n 502 503 RETURN 504 END ! subroutine REVERSE_BITS
