/* Fast Fourier Transform */
/* Adopted from Tao Pang's book */
/* An example of the fast Fourier transform subroutine with   
   n = 2**m.  ar and ai are the real and imaginary part of   
   data in the input and corresponding Fourier coefficients   
   in the output.   */ 

#include <stdio.h>
#include <math.h>
 
void fft (double* ar, double* ai, int n,int m)
{
 int i,j,k,l,n1,n2,l1,l2;
 double pi,a1,a2,q,v,u;
 
 pi = 4*atan(1);
 n2 = n/2;
 
 n1 = (int)pow(2,m);
 if (n1 != n)
 {
  printf("Indices do not match\n");
//  exit(1);
 }
 
/* Rearrange the data to the bit reversed order */
 
 l = 1;
 for (k = 0; k < n-1; ++k)
 {
  if (k < l-1)
  {
   a1 = ar[l-1];
   a2 = ai[l-1];
   ar[l-1] = ar[k];
   ar[k] = a1;
   ai[l-1] = ai[k];
   ai[k] = a2;
  }
  j = n2;
  while (j < l)
  {
   l = l-j;
   j = j/2;
  }
  l = l+j;
 }
 
/* Perform additions at all levels with reordered data */
 
 l2 = 1;
 for (l = 0; l < m; ++l)
 {
  q = 0;
  l1 = l2;
  l2 = 2*l1;
  for (k = 0; k < l1; ++k)
  {
   u =  cos(q);
   v = -sin(q);
   q =  q+pi/l1;
   for (j = k; j < n; j = j+l2)
   {
    i = j+l1;
    a1 = ar[i]*u-ai[i]*v;
    a2 = ar[i]*v+ai[i]*u;
    ar[i] = ar[j]-a1;
    ar[j] = ar[j]+a1;
    ai[i] = ai[j]-a2;
    ai[j] = ai[j]+a2;
   }
  }
   }
}