// root.cern.ch/root/htmldoc/guides/spectrum/Spectrum.html#transforms
// 
// This macro demonstrates the Fourier Transform of signal via TSpectrumTransform

Double_t fun_periodic (Double_t* xval, Double_t* par) {
  Double_t x = xval[0];
  Double_t MaxOrder   = par[0];
  Double_t T0         = par[1];
  Double_t AttenuCoef = par[2];

  Double_t result   = 1. + AttenuCoef;

  for (Double_t order = 1 ; order <= MaxOrder ; order++)
    result += TMath::Power (AttenuCoef, order) * cos (2.*M_PI * x / (T0/order)) ;

  return result;
}


void FourierTransform (Int_t MaxOrder = 6, Double_t T0 = 256, Double_t AttenuCoef = 0.5)
{
  Int_t NBins = 256;

  TF1* FunPeriodic = new TF1 ("FunPeriodic", fun_periodic, 0. , (Double_t) NBins, 3);
  FunPeriodic->SetParameters (MaxOrder, T0, AttenuCoef);
  FunPeriodic->SetNpx (100000);

  TH1F* hsignal = new TH1F ("hsignal", "Input signal", NBins, 0., (Double_t) NBins);
  hsignal->FillRandom ("FunPeriodic", 2000000);
  TCanvas* can1 = new TCanvas ("can1", "Transform demo", 10, 10, 1000, 700);
  can1->Divide (1,2);
  can1->cd (1);
  hsignal->SetMinimum (0.);
  hsignal->Draw ();

  Double_t* source = new Double_t [NBins];
  Double_t* destin = new Double_t [2*NBins];
  for (Int_t i = 0; i < NBins; i++)  source[i] = hsignal->GetBinContent (i+1);

  TSpectrumTransform* Transformer1dim = new TSpectrumTransform ( NBins );

  Transformer1dim->SetTransformType (TSpectrumTransform::kTransformFourier, 0);
/*
  void TSpectrumTransform::SetTransformType (
   Int_t  transType,  : Type of transform (eg. kTransformFourier)
   Int_t  degree      : Degree of mixed transform
  )

  Types of transform: Haar, Walsh, Cosine, Sine, Fourier, Hartley, Fourier-Walsh,
                      Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, 
                      Sine-Walsh, Sine-Haar

  Degree: Applies only for Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh,
                           Cosine-Haar, Sine-Walsh, Sine-Haar transforms
*/

  Transformer1dim->SetDirection (Transformer1dim->kTransformForward); // or kTransformInverse

  Transformer1dim->Transform (source, destin);

/*
  Caution: for Fourier, Fourier-Walsh and Fourier-Haar.
  
  For direct  transforms, the destination array should be 2x larger than source
                          to store  Re and Im coefficients

  For inverse transforms, the source      array should be 2x larger than source
                          to supply Re and Im coefficients

  See: root.cern.ch/doc/master/classTSpectrumTransform.html#a3782d288aab50a3354949c0ef2631972
*/

  TH1F* hresult = new TH1F ("hresult", "Fourier Transform of input", 2*NBins , 0., 2. * NBins);
  for (Int_t i = 0; i < 2*NBins; i++)  hresult->SetBinContent (i+1, destin[i]);
  hresult->SetLineColor (kRed);
  can1->cd (2);
  hresult->SetMarkerStyle (20);
  hresult->SetMarkerSize  (0.5);
  hresult->Draw("P");
}