[comp.lang.c] /**

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;
**