/* 

inv_dft.cc: Inverse Discrete Fourier Transformation 

(re)generates y(t) given its dft.

input file:    dft.dat (created by dft.cc)
input file format:   frequency index \t real part \t imaginary part

output file format: same as input 

*/

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

using namespace std;

const int i_max = 1000;
const double PI = acos(-1.0);
 
int main()
{
  double imag, real, input [2][i_max]
, x, y;
  int i = 0, j, k;
 
  ifstream input_file;
  ofstream output_file;
   
  input_file.open("fft.dat");           //get data
  output_file.open("invers_dft.dat");   //save data
   
   
  while (!input_file.eof() && (i < i_max)){
    input_file >> j >> x >> y ;
    input[0][i] = x;
    input[1][i] = y;
    ++i;
  }
   
  int N = i;   // number of sampled values
   
  for (j=0; j < N; ++j)	       // loop over frequency index
     {
       real = 0.0;                   // clear variables
       imag = 0.0;			
       for (k = 0; k < N; ++k)	     // loop for sum
	 {
	   real += input[0][k]*cos(2*PI*k*j/N)+input[1][k]*sin(2*PI*k*j/N);
	   imag += input[1][k]*cos(2*PI*k*j/N)-input[0][k]*sin(2*PI*k*j/N);
	 }
       output_file << j << "\t" << real/(double)N << "\t" << imag/(double)N << endl;
     }

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