#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
const double PI = 3.141592653589793;  // An good approximation to pi
void FFT(double(*)[2], double(*)[2], const int, int, int );  // A function for computing the DFT via direct computation.

int main( void )
{
  const int N = 32768;  // This is the number of equally spaced points in [0, 2pi]
  int i = 0;  // This is a variable used for counting through loops.
 
  // These are function values.  More specifically, functionVals[j][0] holds the ``real'' part of f(2*pi*j/N)
  // and functionVals[j][1] holds the ``imaginary'' part of f(2*pi*j/N).  The imaginary part is always zero
  // when f: [0, 2pi] -> R is a real values function.
  double functionVals[N][2];  
  double DFT[N][2]; // This array will hold the Discrete Fourier Transform of functionVals.  The jth frequency
  // coefficient is stored in DFT[j]:  DFT[j][0] will have the real part, and DFT[j][1] the imaginary part.
  
  // We will now "load up" our function values using the test function f(x) = cos(1234*x)
  for(i = 0; i < N; i++)
  {
    // Save both the real and imaginary parts of the function
    functionVals[i][0] = cos(1234*2*PI*i / N);
    functionVals[i][1] = 0;
  }
  
  // Let the user know you are starting.
  printf("Beginning FFT...\n");

  // Take the DFT of functionVals by calling the FFT function
  FFT(functionVals, DFT, N, 1, 0);

  // Let the user know you are done.
  printf("Finished FFT.\n");
  
  // Print out the frequencies that have Fourier coefficients whose magnitudes are bigger than .0001
  for( i = 0; i < N; i++)
  {
    // Compute the magnitude of the ith Fourier coefficient and check to see if its bigger than .0001
    if( sqrt(DFT[i][0]*DFT[i][0] + DFT[i][1]*DFT[i][1]) > .0001*N )
    {
      // Print this "large" Fourier coefficient.  
      printf("The Frequency %d has Fourier coefficient %f + %fi.\n",i,DFT[i][0] / N,DFT[i][1] / N);
    }
  }
  
  // Quit the program
  return 0;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// This Function computes the Fourier transform of the Data via recurrsive FFT.  It is a fast O(n log N)-time //
// Algorithm.  NOTE:  N must be a power of 2!!!  It takes five inputs:                                        //
// 1.  The Data which represents samples from a periodic function, f: [0, 2pi] -> C.  Data[j] = f(2 pi j/Size)//
// 2.  The Discrete Fourier Transform of the Data, <f,e^ikx>_Size for k = 0, ..., Size-1, will be stored here.//
// 3.  The size of the entire data array (i.e., the number of equally spaced samples from [0, 2pi] we have.   //
// 4.  The step size we are using the in sub-array of the Data that we are currently FFT-ing.                 //
// 5.  The start position in the sub-array of the Data array that we are currently FFT-ing.                   //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void FFT( double(*Data)[2], double(*DFT)[2], const int N, const int step, const int start )
{
  // Variables for looping
  const int size = N / step; // Size of the sub-array of Data that we are FFTing now.
  double NDFT[size][2];  // For computational purposes...
  const int sizeDiv2 = size / 2;  // Half the current size of the array
  const double h = -2*PI / size;  // The step size we use.
  double fsp;  // A variable used for calculations.
  int k, kmod, oddindex, evenindex, index;  // Variables for looping
  const int twostep = 2*step;  // Twice the current step size

  // Base case for recurrsion
  if( size == 1 ) 
  {
    // FFT of size 1 returns the number back again.
    DFT[start][0] = Data[start][0];
    DFT[start][1] = Data[start][1];
  }
  // We need to compute FFTs of the even and odd sub-arrays of the current sub-array
  else
  { 
    // Take the FFTs of the even and odd sub-arrays
    FFT(Data, DFT, N, twostep, start);
    FFT(Data, DFT, N, twostep, start + step);

    // Peice the even and odd sub-FFTs together again...
    for(k = 0; k < size; k++)
    {
      kmod = k % sizeDiv2;
      oddindex = start + step + twostep*kmod;
      evenindex = start + twostep*kmod;
      fsp = h*k;
      NDFT[k][0] = DFT[evenindex][0] + cos(fsp)*DFT[oddindex][0] - sin(fsp)*DFT[oddindex][1];
      NDFT[k][1] = DFT[evenindex][1] + sin(fsp)*DFT[oddindex][0] + cos(fsp)*DFT[oddindex][1];
    }
    // Copy DFT data into the anaswer array.
    for(k = 0; k < size; k++)
    {
      index = start + k*step;
      //DFT[index][0] = DFT[evenindex][0] + cos(fsp)*DFT[oddindex][0] - sin(fsp)*DFT[oddindex][1];
      //DFT[index][1] = DFT[evenindex][1] + sin(fsp)*DFT[oddindex][0] + cos(fsp)*DFT[oddindex][1];
      DFT[index][0] = NDFT[k][0];
      DFT[index][1] = NDFT[k][1];
    }
  }
  
  return;
}