MISS026%ECNCDC.BITNET@wiscvm.wisc.EDU (GREENY) (07/21/87)
HEADER: TITLE:**FAST FOURIER TRANSFORM; DATE:**05/18/1985; DESCRIPTION:*"PERFORMS FAST FOURIER TRANSFORM USING METHOD DESCRIBED **BY E. O. BRIGHAM. FOR DETAILS OF THE METHOD, REFER **TO BRIGHAM'S BOOK. THE FAST FOURIER TRANSFORM";* KEYWORDS: *FOURIER, TRANSFORM; FILENAME:*FFT.C; WARNINGS:* "THIS PROGRAM IS SELF-CONTAINED, ALL THAT IS NEEDED IS A MANNER OF GETTING THE DATA INTO THE ARRAY REAL_DATA (& IMAG_DATA, IF APPLICABLE). THE TRANSFORMED DATA WILL RESIDE IN THESE TWO ARRAYS UPON RETURN WITH THE ORIGINAL DATA BEING DESTROYED." AUTHORS:*JIM PISANO; COMPILERS:*DESMET; REFERENCES:*AUTHORS:*E. O. BRIGHAM; **TITLE:**"THE FAST FOURIER TRANSFORM"; **CITATION:*""; *ENDREF */ /**FILE NAME FFT.C **PROGRAM NAME FFT() ... FAST FOURIER TRANSFORM * **PERFORM FAST FOURIER TRANSFORM USING METHOD DESCRIBED BY E. O. BRIGHAM. * FOR DETAILS OF THE METHOD, REFER TO BRIGHAM'S BOOK * **TRANSLATED TO C FROM FORTRAN BY * ***JIM PISANO ***P.O. BOX 3134 ***UNIVERSITY STATION ***CHARLOTTESVILLE, VA 22903 * * THIS PROGRAM IS IN THE PUBLIC DOMAIN & MAY BE USED BY ANYONE FOR COMMERCIAL * OR NON-COMMERCIAL PURPOSES. * * REAL_DATA ... PTR. TO REAL PART OF DATA TO BE TRANSFORMED * IMAG_DATA ... PTR. TO IMAG " " " " " " * INV ..... SWITCH TO FLAG NORMAL OR INVERSE TRANSFORM * N_PTS ... NUMBER OF REAL DATA POINTS * NU ...... LOGARITHM IN BASE 2 OF N_PTS E.G. NU = 5 IF N_PTS = 32. * * THIS PROGRAM IS SELF-CONTAINED, ALL THAT IS NEEDED IS A MANNER OF GETTING * THE DATA INTO THE ARRAY REAL_DATA (& IMAG_DATA, IF APPLICABLE). THE * TRANSFORMED DATA WILL RESIDE IN THESE TWO ARRAYS UPON RETURN WITH THE * ORIGINAL DATA BEING DESTROYED. */ #INCLUDE*<STDIO.H> #INCLUDE*<MATH.H>** /* DECLARE MATH FUNCTIONS TO USE */ #DEFINE**PI**3.1419527 FFT( REAL_DATA, IMAG_DATA, N_PTS, NU, INV ) DOUBLE *REAL_DATA, *IMAG_DATA; INT N_PTS, NU, INV; * *INT N2, J, J1, L, I, IB, K, K1, K2; *INT SGN; *DOUBLE TR, TI, ARG, NU1;*/* INTERMEDIATE VALUES IN CALCS. */ *DOUBLE C, S;* /* COSINE & SINE COMPONENTS OF FOURIER TRANS. */ *N2 = N_PTS / 2; *NU1 = NU - 1.0; *K = 0; /* * SIGN CHANGE FOR INVERSE TRANSFORM */ *SGN = INV ? -1 : 1;** /* * CALCULATE THE COMPONETS OF THE FOURIER SERIES OF THE FUNCTION */ *FOR( L = 0; L != NU; L++ ) ** **DO *** ***FOR( I = 0; I != N2; I++ ) ***** ****J = K / ( POW( 2.0, NU1 ) ); ****IB = BIT_SWAP( J, NU ); ****ARG = 2.0 * PI * IB / N_PTS; ****C = COS( ARG ); ****S = SGN * SIN( ARG ); ****K1 = K; ****K2 = K1 + N2; ****TR = *(REAL_DATA+K2) * C + *(IMAG_DATA+K2) * S; ****TI = *(IMAG_DATA+K2) * C - *(REAL_DATA+K2) * S; *****(REAL_DATA+K2) = *(REAL_DATA+K1) - TR; *****(IMAG_DATA+K2) = *(IMAG_DATA+K1) - TI; *****(REAL_DATA+K1) = *(REAL_DATA+K1) + TR; *****(IMAG_DATA+K1) = *(IMAG_DATA+K1) + TI; ****K++; **** ***K += N2; *** WHILE( K < N_PTS - 1); **K = 0; **NU1 -= 1.0; **N2 /= 2; *** *FOR( K = 0; K != N_PTS; K++ ) ** **IB = BIT_SWAP( K, NU ); **IF( IB > K) *** ***SWAP( (REAL_DATA+K), (REAL_DATA+IB) ); ***SWAP( (IMAG_DATA+K), (IMAG_DATA+IB) ); *** ** /* * IF CALCULATING THE INVERSE TRANSFORM, MUST DIVIDE THE DATA BY THE NUMBER OF * DATA POINTS. */ *IF( INV ) **FOR( K = 0; K != N_PTS; K++) *** ****(REAL_DATA+K) /= N_PTS; ****(IMAG_DATA+K) /= N_PTS; *** * /* * BIT SWAPING ROUTINE IN WHICH THE BIT PATTERN OF THE INTEGER I IS REORDERED. * SEE BRIGHAM'S BOOK FOR DETAILS */ BIT_SWAP( I, NU ) INT I, NU; * *INT IB, I1, I2; *IB = 0; *FOR( I1 = 0; I1 != NU; I1++ ) ** **I2 = I / 2; **IB = IB * 2 + I - 2 * I2; **I = I2; ** *RETURN( IB ); * /* * SIMPLE EXCHANGE ROUTINE WHERE *X1 & *X2 ARE SWAPPED */ SWAP( X1, X2 )* DOUBLE *X1, *X2; * *INT *TEMP_IY; *DOUBLE *TEMP_X; **TEMP_X = *X1; **X1 = *X2; **X2 = *TEMP_X; **