/*
 * Purpose : Perform a one dimensional FFT in place.
 *
 * Parameters : input   [char *] array with input values
 *    npoints [int]  length of array.
 *    dir   [int]  if dir =  1 then forward FFT,
 *     if dir = -1 then inverse FFT.
 *
 *
 *
Purpose : Perform a one dimensional FFT in place.

Parameters:
 input [char *] array with input values
 npoints [int]  length of array.
 dir  [int]  if dir =  1 then forward FFT,
     if dir = -1 then inverse FFT.


*/
#include <math.h>
#include "fft.h"

#define PI M_PI


int fft(double *real,double *imag,int npoints,int dir)

{
	double temp;
	register int i,j,k;
	int index,swapindex;
	register double tr,ti,angle,wr,wi;

	/* Swap the input elements for the decimation in time algorithm. */

	for (index = 1, swapindex = 0; index < npoints; index++)
	{
		k = npoints;
		do 
			k /= 2;  
		while ((swapindex + k) >= npoints);
		swapindex = (swapindex % k) + k;
		if (swapindex <= index) continue;
		temp = real[index];
		real[index] = real[swapindex];
		real[swapindex] = temp;
		temp = imag[index];
		imag[index] = imag[swapindex];
		imag[swapindex] = temp;
	}

	/*
   * Do the butterfly computations.
   *
   *   stage index:   k = 1, 2, 4, 8 . . .
   *     For npoints = 8, for example, there will
   *     be three stages of butterfly computations.
   *   b'fly indices:  i and j will be separated by a distance
   *     dependent on the current stage.
   *     k is used as the separation constant.
   */

	for (k = 1 ; k < npoints ; k *= 2)
	{
		for (index = 0 ; index < k ; index++)
		{
			angle = PI * ((double) (index * -dir)) / ((double) k);
			wr = cos(angle);
			wi = sin(angle);
			for (i = index ; i < npoints ; i += 2*k)
			{
				j = i + k;
				tr = (wr * real[j]) - (wi * imag[j]);
				ti = (wr * imag[j]) + (wi * real[j]);

				real[j] = real[i] - tr;
				imag[j] = imag[i] - ti;
				real[i] += tr;
				imag[i] += ti;
			}
		}
	}

	wr = sqrt(1./npoints);

	/* scale the output by sqrt(1./N) */

	for (i = 0 ; i < npoints ; i++)
	{
		real[i] *= wr;
		imag[i] *= wr;
	}

	return (0);
}
/*+
 reorders elements of an block for a real fft

 programmer:   L v.d. Wal
 rewritten in C by:  Frits v.d. Putten
 date:    sept 1980
     April 1988 (conversion to C)

 ****************************************************************

 input variables :
  c   complex data
  n   number of elements in c (power of 2)

 output variables :
  c   complex data
-*/
/*
shufft(COMPLEX *c,int n,int dir)

{
   register int i, n2;
   register double scale;
   double d1;
   double d2;
   COMPLEX cda1;
   COMPLEX cda2;
   COMPLEX w1;
   COMPLEX w2;

   d1 = c[0].real;
   d2 = c[0].imag;
   c[0].real = d1+d2;
   c[0].imag = d1-d2;
   w1.real = w2.real = cos(PI/n);
   w1.imag = w2.imag = -dir*sin(PI/n);
   n2 = n/2;
   if (dir == 1)
      scale = 0.5;
   else
      scale = 1.0;

   for (i=1; i<=n2; i++)
   {
     cda1.real  = scale * (c[i].real + c[n-i].real);
     cda1.imag  = scale * (c[i].imag - c[n-i].imag);
     cda2.real  = scale * (c[n-i].imag + c[i].imag);
     cda2.imag  = scale * (c[n-i].real - c[i].real);
     c[n-i].real = cda1.real - dir*(cda2.real*w1.real - cda2.imag*w1.imag);
     c[n-i].imag = dir*(cda2.real*w1.imag + cda2.imag*w1.real) - cda1.imag;
     c[i].real  = cda1.real + dir*(cda2.real*w1.real - cda2.imag*w1.imag);
     c[i].imag  = cda1.imag + dir*(cda2.real*w1.imag + cda2.imag*w1.real);
     d1   = w1.real * w2.real - w1.imag * w2.imag;
     d2   = w1.real * w2.imag + w1.imag * w2.real;
     w1.real  = d1;
     w1.imag  = d2;
   }

   return;
}
 */
