/* An FFT routine */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "complexdefs.h"
#define PI 3.1415926536

/* Globally declared array will hold complex exponentials.
   The array must be initialized with the routine initialize_exps below before
   using Fourier transform.  */

Complex *exps;

void initialize_exps(int N)
{
  int i;

  exps=malloc(N*sizeof(Complex));
  for(i=0;i<N;i++)
    {
      exps[i].Re=cos((2*PI*i)/N);
      exps[i].Im=sin((2*PI*i)/N);
    }
}

/*  If the parameter mode is 1, the FFT takes an array
    of N complex numbers,
    allocated in the calling routine, and replaces
    it by its Discrete Fourier Transform. If mode
    is -1, the inverse transform is computed. 
    N is required to be a power of 2; the routine
    does not check for this!
*/

void fft(Complex *sequence,int N,int mode)
{
  int i,j,x,y,m,base1,base2,arraysize,blocksize,group,numgroups;
  Complex temp,multiplier,negmultiplier;
  /* bit reversal */
  for(j=0;j<N;j++)
    {
      m=N;
      x=j;
      y=0;
      while(m!= 1)
	{
	  y=y<<1;
	  if(x&1)
	    y=y|1;
	  x=x>>1;
	  m=m>>1;
	}
      if(j<y)
	{
	  temp=sequence[j];
	  sequence[j]=sequence[y];
	  sequence[y]=temp;
	}
    }
  /* Now do log2(N) passes of coalescence */
  /* In each pass we have numgroups pairs of subarrays of size
     blocksize, and coalesce them into subarrays of size
     arraysize */
  arraysize=2;
  blocksize=1;
  numgroups=N/2;
  while(arraysize<=N)
    {
      for(group=0;group<numgroups;group++)
	{
	  base1=arraysize*group;
	  base2=base1+blocksize;
	  for(i=0;i<blocksize;i++)
	    {
	      if(mode==1)
		{
		  multiplier = exps[(i*numgroups)%N];
		  negmultiplier=exps[(i*numgroups+N/2)%N];
		}
	      else
		{
		  multiplier=exps[(N-i*numgroups)%N];
		  negmultiplier=exps[(N/2-i*numgroups)%N];
		}
	      temp=Cadd(sequence[base1+i],Cmultiply(negmultiplier,sequence[base2+i]));
	      sequence[base1+i]=Cadd(sequence[base1+i],Cmultiply(multiplier,sequence[base2+i]));
	      sequence[base2+i]=temp;
	    }
	}
      blocksize=blocksize<<1;
      numgroups=numgroups>>1;
      arraysize=arraysize<<1;
    }
  if(mode!=1)
    for(j=0;j<N;j++)
      {
	sequence[j].Re=sequence[j].Re/N;
	sequence[j].Im=sequence[j].Im/N;
      }
}