#ifndef FSAPPIRRSMPL_SEEN
#define FSAPPIRRSMPL_SEEN


// approximation par serie de Fourier avec echantillonnage irregulier
//------------------------------------------------------------------
//   *****   Guy Le Meur -- LAL-Orsay mars 2002 ****************
//-------------------------------------------------------------------



#include <math.h>
#include <iostream.h>


#include "machdefs.h"    // Definitions specifiques SOPHYA
#include "nbmath.h"
#include "timing.h"
#include "array.h"
#include "fftservintf.h"
#include "fftpserver.h"

#include "toeplitzMatrix.h"


////////////////////////////////////////////////////////////////////
// Un signal donne , suppose a bande de frequences limitee, de longueur 
// finie est periodise. Ce signal periodise est suppose developpe 
// en serie de Fourier finie: l'approximation consiste a rechercher 
// les coefficients de cette serie de Fourier. Cela est fait en utilisant 
// le fait que la matrice du systeme a ecrire est Toeplitz. On utilise
// une classe ToeplitzMatrix, qui utilise les FFT pour la resolution
// L'utilisation standard comporte les etapes suivantes : 
//
//   constructeur --> definit l'echantillonnage et l'amplitude
//                    des abscisses (pour periodisation)
//
//   methode approximateSignal(nbFreq, signal) : definit la bande
//                    de frequences et les valeurs du signal
//
//   recuperation de valeurs approximees ou interpolees : 
//   methodes:  . restaureSignal() (signal "debruite" aux valeurs 
//                                              d'echantionnage) 
//              . restaureRegularlySampledSignal() (signal recalcule
//                          sur un echantillonnage regulier quelconque)
//              . computeSignalOnASampling()  (signal recalcule
//                          sur un echantillonnage irregulier quelconque)
//
//  on peut aussi recuperer les valeurs des coefficients du developpement
//  en serie de Fourier (methodes complexFourierCoef() et coeffCosSin()
/////////////////////////////////////////////////////////////////////


class FSApproximationIrregularSampling
{

private:
  // verouiller le clonage
FSApproximationIrregularSampling(const FSApproximationIrregularSampling&) {}
FSApproximationIrregularSampling &operator = (const FSApproximationIrregularSampling&) {return *this;}


public:
FSApproximationIrregularSampling();

FSApproximationIrregularSampling(TVector<double>& sampling, double offset, double range);

~FSApproximationIrregularSampling(); 

void approximateSignal(int M, const TVector<double>& signal);

void restaureRegularlySampledSignal(int n, TVector<double>& solution) const;


// le vecteur d'abscissses est suppose ordonne et appartenir a l'intervalle 
// de definition du signal
 void computeSignalOnASampling(const TVector<double>& abscisses, TVector<double>& solution ) const;

double estimationConditionnement() const;

void  samplingValues(TVector<double>& sv) const; 

void restaureSignal(TVector<double>& sol) const; 


inline void sampledSignal(TVector<double>& signal) const 
{
  signal = signal_;  
  reshapeSignalInUsersFrame(samplingValues_, signal);
}
inline const TVector<double>& weights() const { return poids_;}

// coefficients complexes ck, pour k=0,M
inline void complexFourierCoef(TVector<complex<double> >& coef) const
{
  int n= coefFourier_.Size();
  coef.ReSize(n);
  coef = coefFourier_;
  coef *=  normeSignal_;
}
// terme constant, puis cos, sin, cos, sin......
void coeffCosSin(TVector<double>& coef) const;

private:

inline void initFFT()
{
  fftIntfPtr_ =new FFTPackServer;
  fftIntfPtr_->setNormalize(false);
    
}
inline double valeursSerie(double u) const
{
  complex<double> somme =complex<double>(0.,0.);
  for (int j=1; j<=M_; j++)
    {
      double angle = 2.*M_PI*j*u;
      complex<double> expon = complex<double>(cos(angle), sin(angle));
      somme += coefFourier_(j)*expon;
    }
  return coefFourier_(0).real()+ 2*somme.real();
}

inline void calculeExponentiellesFourier()
{
  int j;
  int n = samplingValues_.Size();
  exponFourier_.ReSize(n);
  for (j=0; j<n; j++)
    {
      double angle=-2.*M_PI*samplingValues_(j);
      exponFourier_(j) = complex<double>(cos(angle),sin(angle));
    }
}

void matchToSamplingReference(TVector<double>& sampling) const;
void resizeSamplingIn_0_1(double offset, double range); 

void reshapeSignalInUsersFrame(const TVector<double>& abscisses, TVector<double>& resultat) const;
void reshapeSignalInUsersFrame(TVector<double>& resultat) const;
 
void makeToeplitzMatrix(int M);
void makeRHS(TVector<complex<double> >& coefSolution);
void makeSamplingVector(const TVector<double>& sampling, double offset, double range);
void makeSignalVector(const TVector<double>& signal);
void computeWeights();
void NormSignal();


FFTServerInterface* fftIntfPtr_;
Toeplitz tptz_;
TVector<double> samplingValues_;
TVector<double> poids_;
TVector<double> signal_;
TVector<complex<double> > exponFourier_;
TVector<complex<double> > coefFourier_;
double samplingOffset_;
double samplingRange_;
double normeSignal_;
double delta_;
int nokdelta_;
int M_;
};


#endif