#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int main (void)
{
  int i;
  const int n = 256;
  double data[2*n];
  
  double omega1=50;
  double omega2=75;
  double phase1=-M_PI*30.0/180.0;
  double phase2=M_PI*90.0/180.0;
  
  double amplitude0=2.0;
  double amplitude1=3.0;
  double amplitude2=1.50;

  double t=0.0;  
  
  gsl_fft_complex_wavetable * wavetable;
  gsl_fft_complex_workspace * workspace;
  
  for (i = 0; i < n; i++)
  {
     double fs =2.0*M_PI;
     double fn = fs/n;
     
    t= fn*i;   
    REAL(data,i) = amplitude0+amplitude1*cos(omega1*t+phase1)+amplitude2*cos(omega2*t+phase2);
    IMAG(data,i) = 0.0;
    char str[64]="123.45";
    //double test =  atof(str);
  }
  
  
  wavetable = gsl_fft_complex_wavetable_alloc (n);
  workspace = gsl_fft_complex_workspace_alloc (n);
  
  gsl_fft_complex_forward (data, 1, n, wavetable, workspace);  

  for (i = 0; i < n; i++)
  {    
   
      double R = sqrt(REAL(data,i)*REAL(data,i) + IMAG(data,i)*IMAG(data,i));
      double amplitude =0;    
      double phase = 0.0;
      phase = atan2(IMAG(data,i),REAL(data,i));    
      
      if(0 == i)
      {
        amplitude = R/n;
      }
      else
      {
         amplitude = R/n*2;
      }    
      printf ("(%d,%f,%f) (%f,%f)\n", i, amplitude,phase,REAL(data,i),IMAG(data,i));
 
  } 
  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);
  return 0;
}