1C -------------------------------------------------------------------------
2C | File name: fft.for Sean O'Connor 9-23-87 |
3C -------------------------------------------------------------------------
4
5C -------------------------------------------------------------------------
6C | VAX 11/785 VMS 4.5 VAX EXTENDED FORTRAN 77 |
7C -------------------------------------------------------------------------
8
9C
10C ****************************************************************************
11C * *
12C * FFT *
13C * *
14C ****************************************************************************
15C
16C
17C DESCRIPTION
18C
19C This program does an in-place one dimensional complex discrete Fourier
20C transform. It uses an FFT algorithm for a number of input points which
21C is a power of two.
22C
23C
24C INPUT
25C
26C N, The number of complex points to be transformed.
27C N is a power of 2. N must be at least 2; it
28C must not exceed the dimension of array X, or the
29C maximum size of the sine and cosine table. (See
30C below).
31C
32C X, The array holding the complex data X( k )
33C for k = 0, 1, . . ., N - 1 which is to be
34C transformed. It is dimensioned in the calling
35C program. Only the array address is passed to
36C this subroutine.
37C
38C D, The direction of the FFT. D is 1 for a
39C forward transform and -1 for an inverse
40C transform.
41C
42C
43C OUTPUT
44C
45C X, On output, X holds the complex Fourier coefficients
46C ^
47C X( k ) for k = 0, 1, . . ., N - 1.
48C
49C Warning messages to the terminal under these conditions:
50C
51C (1) D is not 1 or -1
52C (2) N is less than 2.
53C (3) N is greater than the sine table can handle.
54C (4) N is not a power of 2.
55C
56C
57C CALLING SEQUENCE
58C
59C CALL FFT ( N, X, D )
60C
61C
62C EXAMPLE
63C
64C DIMENSION X( 0 : 3 ) ! 4 point FFT
65C
66C N = 4 ! 4 point FFT
67C D = 1 ! Forward transform
68C
69C CALL FFT( N, X, D ) ! Do the FFT
70C
71C D = -1 ! Do an inverse FFT
72C
73C CALL FFT( N, X, D )
74C
75C Suppose X contains the data
76C
77C X ( 0 ) = 1 + i
78C X ( 1 ) = 2 + 2 i
79C X ( 2 ) = 3 + 3 i
80C X ( 3 ) = 4 + 4 i
81C
82C then the forward FFT returns
83C
84C X ( 0 ) = 5 + 5 i
85C X ( 1 ) = -2
86C X ( 2 ) = -1 - i
87C X ( 3 ) = - 2 i
88C
89C and the inverse FFT returns the original data.
90C
91C
92C METHOD
93C ^
94C The Fourier coefficients X( k ) are defined for k = 0, . . ., N - 1
95C by the formula
96C
97C N - 1 2 Pi i
98C ^ 1 ____ - ------ j k
99C X( k ) = --------- \ N
100C ----- / X( j ) e
101C \ / N ----
102C v j = 0
103C
104C where e is 2.1718 . . ., Pi = 3.14159. . . and i is the square root
105C of -1.
106C
107C This is the formula for the forward transform. The inverse transform
108C is the same, except that the minus sign in the exponent above is a
109C plus sign instead.
110C ---
111C The normalizing factor 1 / \/ N in front was picked so that the
112C forward and inverse transforms look the same.
113C
114C
115C The particular type of FFT algorithm we use is a
116C
117C complex,
118C two to the N,
119C decimation in frequency,
120C input in order, output bit reversed
121C
122C type of FFT. It is described in the book
123C
124C THE FAST FOURIER TRANSFORM,
125C F. Oran Brigham, Prentice Hall, 1974
126C
127C
128C BUGS
129C
130C The bit-reversal routine could be made faster by using the
131C VAX FORTRAN 77 library routines for testing, setting, clearing and
132C shifting bits: BTEST, IBSET, IBCLR, ISHFT.
133C
134C You might want to use double precision complex numbers to get more
135C accuracy.
136C
137C ****************************************************************************
138C * *
139C ****************************************************************************
140
141
142 SUBROUTINE FFT( N, X, D )
143
144 IMPLICIT NONE ! Force declaration of all variables.
145
146
147C -------------------------------------------------------------------------
148C | MACRO DEFINITIONS |
149C -------------------------------------------------------------------------
150
151 INTEGER MAX_TABLE_SIZE
152 REAL PI
153
154 PARAMETER ( PI = 3.1415926535897932385 )
155
156C We must set the size of the sine and cosine table used in the FFT.
157C The table is called starting_multiplier.
158C
159C Suppose MAX is the largest FFT we would ever want
160C to perform. Then
161C
162C MAX_TABLE_SIZE = LOG ( MAX ) - 1
163C 2
164C
165C is the size of the table. This table can be used by the
166C FFT for any number of points from 2 up to MAX.
167C
168C For example, if MAX_TABLE_SIZE = 14, then we can transform
169C anywhere from 2 to 2 ** 15 = 32,768 points, using the same
170C sine and cosine table.
171C
172
173 PARAMETER ( MAX_TABLE_SIZE = 20 )
174
175
176C -------------------------------------------------------------------------
177C | CALLING ARGUMENTS |
178C -------------------------------------------------------------------------
179
180 INTEGER N ! Number of points in transform
181
182 COMPLEX X( 0 : * ) ! Complex data to be transformed
183
184C Array X is dimensioned in the calling program. Only its starting
185C 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
191C -------------------------------------------------------------------------
192C | LOCAL VARIABLES |
193C -------------------------------------------------------------------------
194
195C A column refers to a column of the FFT signal flow graph. A group is
196C a collection of butterflies which have the same multiplier W^p. A
197C 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
228C -------------------------------------------------------------------------
229C | SUBROUTINE CODE |
230C -------------------------------------------------------------------------
231
232C
233C Initialize the variable called_already to .FALSE. at compile time.
234C Change it to .TRUE. after the first call to this subroutine.
235C Save the value so it does not get destroyed after returning from
236C this subroutine.
237C
238 SAVE called_already
239 DATA called_already /.FALSE./
240
241
242C
243C Compute log base 2 of N. Perturb slightly to avoid roundoff error.
244C For example, INT( 4.99999999999 ) = 4 (wrong!), but
245C INT( 4.999999999 + .01 ) = 5 (correct!)
246C
247 log_of_n = INT( ALOG( FLOAT( N ) ) / ( ALOG( 2.0 ) ) + .01 )
248
249
250
251C
252C Error checking. Print a message to the terminal if any one of these
253C conditions is true:
254C
255C (1) D is not 1 or -1.
256C (2) N < 2.
257C (3) log2( N ) > MAX_TABLE_SIZE + 1.
258C (4) N is not a power of 2.
259C
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
287C
288C If called_already is false, generate the sine and cosine table for
289C the first time, then set called_already to .TRUE. Now later calls to
290C the subroutine don't recreate the sine and cosine table.
291C
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
308C
309C Compute the FFT according to its signal flow graph with
310C computations arranged into butterflies, groups of butterflies
311C with the same multiplier, and columns.
312C
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
321C
322C Modify starting multiplier W^p when computing the inverse transform.
323C
324 IF ( D .EQ. -1 ) THEN
325
326 multiplier = CONJG( multiplier )
327
328 END IF
329
330 increment = multiplier
331
332
333C
334C Compute the first group.
335C
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
347C
348C Now do the remaining groups.
349C
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
357C 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
380C
381C Normalize the results and permute them into order, two at a time.
382C
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
410C ****************************************************************************
411C * *
412C * REVERSE_BITS *
413C * *
414C ****************************************************************************
415C
416C DESCRIPTION
417C
418C This function reverses the order of the bits in a number.
419C
420C
421C INPUT
422C
423C n The number whose bits are to be reversed.
424C
425C num_bits The number of bits in N, counting leading zeros.
426C
427C
428C OUTPUT
429C
430C REVERSE_BITS The number N with all its num_bits bits reversed.
431C
432C
433C EXAMPLE
434C
435C REVERSE_BITS( (10) , 2 ) = (01) = 1, but
436C 2 2
437C
438C REVERSE_BITS( (10) , 3 ) = REVERSE_BITS( (010) , 3 ) = (010) = (10)
439C 2 2 2 2
440C
441C = 10
442C
443C
444C METHOD
445C
446C The algorithm works by a method similar to base conversion.
447C As the low order bits are broken off, one by one, they are
448C multiplied by 2 and accumulated to generate the bit-reversed
449C number.
450C
451C
452C BUGS
453C
454C The bit-reversal routine could be made faster by using the
455C VAX FORTRAN 77 library routines for testing, setting, clearing and
456C shifting bits: BTEST, IBSET, IBCLR, ISHFT.
457C
458C
459C ****************************************************************************
460C * *
461C ****************************************************************************
462
463 INTEGER FUNCTION REVERSE_BITS( n, num_bits )
464
465 IMPLICIT NONE
466
467
468C -------------------------------------------------------------------------
469C | CALLING ARGUMENTS |
470C -------------------------------------------------------------------------
471
472 INTEGER n, num_bits
473
474
475C -------------------------------------------------------------------------
476C | LOCAL VARIABLES |
477C -------------------------------------------------------------------------
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
486C -------------------------------------------------------------------------
487C | SUBROUTINE CODE |
488C -------------------------------------------------------------------------
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
