/* 

dft.cc: Discrete Fourier Transformation 

input file format: real y(t) values separated by whitespace
input file name:   fourier_raw.dat

output file name:     dft.dat
output file format:   frequency index \t real paret \t imaginary part

related programs: inv_dft.cc 

*/


#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;


const int i_max = 1000;		     // max number of input data 
const double PI = acos(-1.0);        //  let the machine define pi


int main()
{
  double imag, real, input[i_max], x;
   int i=0,j,k;

 
   ifstream input_file;
   ofstream output_file;

   input_file.open("fourier_raw.dat");   // get input data
   output_file.open("dft.dat");          // save processed data


   for (i = 0; i < i_max; ++i) // zero out working array
     {
      input[i] = 0.0;
     }

   i = 0; // reset i

   while (!input_file.eof() && (i < i_max)){  // fill sampled data array
     input_file >> j >> x;
     input[i] = x;
     ++i;
   }

  
   int N = i; // number of data points read

   for (j=0; j < N ; ++j)		  // loop over frequency index 
     {
       real = 0.0;
       imag = 0.0;		         //clear variables 
       for (k = 0; k < N; ++k)		 // loop for sums 
	 {
	   real += input[k]*cos((2*PI*k*j)/N);
	   imag += input[k]*sin((2*PI*k*j)/N);
	 }
    
       output_file << j << "\t" << real << "\t" << imag << endl;
     }

   input_file.close();
   output_file.close();
   cout << "data stored in dft.dat" << endl;
   return 0;
}