//******************************************************************************
// Calcul d'une FFT avec le librairie GSL
// Dominique Lefebvre TangenteX.com Novembre 2009
//******************************************************************************
#include <cstdlib>
#include <iostream>
#include <math.h>

#include <gsl\gsl_errno.h>
#include <gsl\gsl_fft_complex.h>

#include <dislin_d.h>

#define REAL(z,i)  ((z)[2*(i)])
#define IMAG(z,i)  ((z)[2*(i)+1])

#define N          128

using namespace std;

//******************************************************************************
//routines de service
//******************************************************************************
void InitGraphDislin(void)
{
    metafl("XWIN");
    disini();
    pagera();
    pagfll(255);
    color("black");
    winfnt("arial");
    return;     
}

//******************************************************************************
// programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
    int i;
    double data[2*N];
    
    // Mise a 0 des tableaux;
    for (i=0;i<N;i++)
     {
        REAL(data,i) = 0.0;
        IMAG(data,i) = 0.0;
     }

    // Initialisation de la librairie graphique
    InitGraphDislin();
    wintit("FFT d'une impulsion");
    graf(0, N-1, 0, 20,-0.5,2,-0.5,0.5);
    titlin("Transformee de Fourier d'une impulsion",1);
    title();
            
    // Initialisation de la forme de l'impulsion 
    REAL(data,0) = 1.0;
    for(i=1;i<=10;i++) 
      REAL(data,i) = REAL(data,N-i) = 1.0;

    // Tracé de l'impulsion
    color("blue");
    rlstrt(0,REAL(data,0));                        // positionnement à l'origine
    for(i=1;i<N;i++)
       rline(i-1,REAL(data,i-1),i,REAL(data,i));     

    // Calcul de la FFT
    gsl_fft_complex_radix2_forward(data,1,N);

    // Tracé de la FFT
    color("red");
    rlstrt(0,REAL(data,0)/sqrt(N));                // positionnement à l'origine
    for(i=1;i<N;i++)
       rline(i-1,REAL(data,i-1)/sqrt(N),i,REAL(data,i)/sqrt(N)); 

    disfin();
    return EXIT_SUCCESS;
}