/*========================================================================
Filename:	FFT.c
Language:	C

FFT and associated functions.

Modifications:

02/02/91	Gerry Beauregard	Switched to on-the-fly computation of Wkn/N
								factors to reduce memory requirements.
02/23/90	Gerry Beauregard	Created this file
=====================================================================*/

/*====================================================================
Includes
=====================================================================*/

#ifndef		_STDIO_H
#include	<stdio.h>
#endif

#ifndef		_MATH_H
#include	<math.h>
#endif

#ifndef		_COMPLEXMATH_H
#include	"ComplexMath.h"
#endif

#ifndef		_H_FFT
#include	"BeauregardFFT.h"
#endif

/*========================================================================
Prototypes
==========================================================================*/


/*=========================================================================
Globals
============================================================================*/

static Real	WAngleIncrement;	/* Used in computation of Wkn/N factors */
 
/*--FFT---------------------------------------------------------------
Performs a decimation in frequency FFT or IFFT on an array of
complex numbers.  Note that the output is in bit reversed order!
The arguments are:
	X[]			An array of Complex numbers. Used for input and 
				output. 
	LogN		Log to base 2 of array size.
	Direction	The direction of the FFT. 

These constants are useful for setting Direction:
	FORWARD (0)		Do forward FFT.
	INVERSE (1)		Do inverse FFT.		

No return values.  The output array is in the same variable as the input.
-----------------------------------------------------------------------*/
void BeauregardFFT(Complex* X, int LogN, int Direction)
{
	int		Stage;			
	int		FliesPerFFT;	/* Number of butterflies per sub FFT */
	int		Span;			/* Width of the butterfly */
	int		FFTStart;		/* Starting index of sub FFT */
	int		FFTSpacing;		/* Distance between start of sub FFTs */
	int		FlyCount;		/* Counter for number of butterflies */
	int		Top,Bottom;		/* Top and bottom index of a butterfly */
	int		k;				/* General purpose index */
	int		N;				/* Array size */
	int		WIndex;			/* Index for Wkn/N factors */
	int		WIndexStep;		/* Increment for index into WTable */
	Real				XTopRe;	/* = X[Top] */
	Real				XTopIm;	
	Real		XBotRe;	/* = X[Bottom] */
	Real		XBotIm;	
	Real		angle;	/* Angle of Wkn/N factor */
	Real		WRe;	/* Wkn/N factor */ 
	Real		WIm;		

	N = (int)(pow(2.0,(Real)LogN) + 0.001);

	if ( Direction == FORWARD )
		WAngleIncrement = -2.0*MY_PI/(Real)N;
	else
		WAngleIncrement = +2.0*MY_PI/(Real)N;
		
	FliesPerFFT	= N >> 1;	/* N/2 */
	Span 		= N >> 1;
	FFTSpacing	= N;
	WIndexStep	= 1; 	  

	for ( Stage = 0; Stage < LogN; ++Stage )
	{
		FFTStart = 0;
		for (FFTStart = 0; FFTStart < N; FFTStart += FFTSpacing)
		{
			WIndex = 0;
			Top    = FFTStart;
			Bottom = Top + Span;
			for (FlyCount = 0; FlyCount < FliesPerFFT; ++FlyCount)
			{
				angle = WAngleIncrement * WIndex;
				WRe = cos(angle);  
				WIm = sin(angle);  

				XTopRe = X[Top].Re;  
				XTopIm = X[Top].Im;
				XBotRe = X[Bottom].Re;
				XBotIm = X[Bottom].Im;

				X[Top].Re = XTopRe + XBotRe;
				X[Top].Im = XTopIm + XBotIm;
				XBotRe    = XTopRe - XBotRe;
				XBotIm    = XTopIm - XBotIm;
				X[Bottom].Re = XBotRe*WRe - XBotIm*WIm;
				X[Bottom].Im = XBotRe*WIm + XBotIm*WRe;

				++Top;
				++Bottom;
				WIndex += WIndexStep;
			}
		}
		FliesPerFFT		>>= 1;    	/* Divide by 2 by right shift */
		Span		 	>>= 1;
		FFTSpacing 		>>= 1;
		WIndexStep	 	<<= 1;		/* Divide by 2 by left shift */

	}

	/* The algorithm leaves the result in a scrambled order. This unscrambles it */
	BitReverseArray(X,LogN);
	
	/* Divide everything by N for IFFT */
	if (Direction != FORWARD)
		for (k = 0; k < N; k++)
		{
			X[k].Re = X[k].Re/(Real)N;
			X[k].Im = X[k].Im/(Real)N;
		}	
/*	printf("...Finished FFT\n"); */
}			


/*--BitReverseArray-------------------------------------------------------
Takes array elements and puts them in bit-reversed order.
-------------------------------------------------------------------*/
void BitReverseArray(Complex* X, int LogN)
{
	int	Counter;
	int	BitReverseCounter;
	int	N;
	Complex	Temp;
	
	N = (int)(pow(2.0,(Real)LogN) + 0.0001);
	for (Counter = 0; Counter < N; Counter++)
	{
		BitReverseCounter = BitReverse(Counter,LogN);
		if (BitReverseCounter > Counter)
		{
			Temp = X[Counter];
			X[Counter] = X[BitReverseCounter];
			X[BitReverseCounter] = Temp;
		}
	}
}

/*--BitReverse-------------------------------------------------------
Do bit reversal of specified number of places of an int
-------------------------------------------------------------------*/
int BitReverse(int Number, int Places)
{
	int	i;
	int	Reverse;

	Reverse = 0;
	for (i = 0; i < Places; i++)
	{
		Reverse = Reverse << 1;
		Reverse = Reverse | (Number & 0x0001);
		Number = Number >> 1;
	}
	return Reverse;
}


/*--RealToComplex-------------------------------------------------------
Transform a real array to a complex array.
For N entries in Real Array there are N Complex entries
and 2 N slots worth of Reals used.
-------------------------------------------------------------------*/
void RealToComplex(Real* Y, Complex* X, int logN)
{
	int	i;
	int size;
	
	size = 1>>logN;
	
	for (i = 0; i < size; i++)
	{
		X[i].Re = Y[i];
		X[i].Im = 0.0 ;
	}
}

/*--ComplexToReal-------------------------------------------------------
Transform a Complex array to a Real array.
For N entries in Real Array there are N Complex entries
and 2 N slots worth of Reals used.
-------------------------------------------------------------------*/
void ComplexToReal(Complex* X, Real* Y, int logN)
{
	int	i;
	int size;
	
	size = 1>>logN;
	
	for (i = 0; i < size; i++)
	{
		Y[i] = X[i].Re;
	}
}