#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
const double PI = 3.141592653589793;  // An good approximation to pi
void DirectDFT(double(*)[2], double(*)[2], const 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 direct DFT.  Prepare to wait...\n");

  // Take the direct DFT of functionVals by calling the `DirectDFT' function
  DirectDFT(functionVals, DFT, N);

  // Let the user know you are done.
  printf("Finished direct DFT.\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 )
    {
      // Print this "large" Fourier coefficient.  
      printf("The Frequency %d has Fourier coefficient %f + %fi.\n",i,DFT[i][0],DFT[i][1]);
    }
  }
  
  // Quit the program
  return 0;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// This Function computes the Fourier transform of the Data via direct summation.  It is the slow O(N^2)-time //
// Algorithm.  It takes three 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, ..., N-1, will be stored here.   //
// 3.  The size of the Data array (i.e., the number of equally spaced samples from [0, 2pi] we are using.     //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void DirectDFT( double(*Data)[2], double(*DFT)[2], const int Size )
{
  // Variables for looping
  int j,k = 0;
  const double step = 2*PI / Size;  // The step size we use.
  double temp;  // A variable used for calculations.

  // Compute <f,e^ikx>_Size for each k one-by-one
  for( k =0; k < Size; k++)
  {
    // Compute <f,e^ikx>_Size for this k.  Begin by initializing the temp and DFT variables we need.
    temp = step*k;  // This is now 2*pi*k / N
    DFT[k][0] = 0;  // Real part of kth Fourier coefficient sum starts at 0
    DFT[k][1] = 0;  // Imaginary part of kth Fourier coefficient sum starts at 0
    
    // Compute the kth Fourier coefficient
    for(j = 0; j < Size; j++)
    {
      // Add f(2*pi*j/N)*exp(2*pi*i*j*k/N) to the current sum using Euler's formula
      DFT[k][0] += (Data[j][0]*cos(temp*j) - Data[j][1]*sin(temp*j));
      DFT[k][1] += (Data[j][0]*sin(temp*j) + Data[j][1]*cos(temp*j));
    }
    
    // Rescale the Fourier coefficient by 1/Size
    DFT[k][0] /= Size;  
    DFT[k][1] /= Size;  
  }

  // Return control back to the calling function
  return;
}