FFTtools.h

#ifndef FFTTOOLS_H
#define FFTTOOLS_H



// If the following is uncommented, in and out arrays for FFTW are allocated contiguously. 
#define FFTTOOLS_ALLOCATE_CONTIGUOUS 
//#define FFTTOOLS_THREAD_SAFE 
//


// You may choose to define FFTTOOLS_COMPAT_MODE
// to get FFTtools to behave like it did several years ago (where it was a class with static methods). 
// This might help with legacy code. 


// c++ libraries thingies
#include <map>
#include <algorithm>
#include <iostream>
#include <fstream>

// ROOT
#include "TGraph.h"
#include "TROOT.h"
#include "TMath.h"
#include "TString.h"
#include "TSystem.h"
#include "TMath.h"
#include "Math/Interpolator.h"
#include "Math/InterpolationTypes.h"

#ifdef FFTTOOLS_COMPAT_MODE
#define FFTTOOLS_METHOD static
#else
#define FFTTOOLS_METHOD  
#include "FFTWindow.h" 
#endif 


//My includes
#include "FFTWComplex.h"
#include "RFSignal.h"

class TRandom; 
class TH2; 

// FFTW
#include <complex>

#ifdef FFTTOOLS_COMPAT_MODE
class 
#else
namespace 
#endif 
FFTtools
{
    
#ifdef FFTTOOLS_COMPAT_MODE
  public: 
#endif
  

   FFTTOOLS_METHOD TGraph *getInterpolatedGraph(const TGraph *grIn, Double_t deltaT);


   FFTTOOLS_METHOD TGraph *getInterpolatedGraphFreqDom(const TGraph *grIn, Double_t deltaT);


   FFTTOOLS_METHOD TGraph *getConvolution(const TGraph *grA, const TGraph *grB);

   FFTTOOLS_METHOD RFSignal *getConvolution(RFSignal *grA, RFSignal *grB);


   FFTTOOLS_METHOD double getAbs(FFTWComplex &theNum);


   FFTTOOLS_METHOD double *doInvFFT(int length, const FFTWComplex *theInput);

   FFTTOOLS_METHOD FFTWComplex *doFFT(int length,const double *theInput);    

   

   FFTTOOLS_METHOD void doFFT(int length, const double * properly_aligned_input, FFTWComplex * properly_aligned_output); 

   FFTTOOLS_METHOD void doInvFFTClobber(int length, FFTWComplex * properly_aligned_input_that_will_likely_be_clobbered, double * properly_aligned_output); 

   FFTTOOLS_METHOD void doInvFFTNoClobber(int length, const FFTWComplex * properly_aligned_input, double * properly_aligned_output); 


   FFTTOOLS_METHOD TGraph *convertMagnitudeToTimeDomain(const TGraph *inputMag);    


   FFTTOOLS_METHOD double *getCorrelation(const TGraph *gr1, const TGraph *gr2,int firstIndex,int lastIndex);

   FFTTOOLS_METHOD double *getCorrelation(int length,const float *oldY1, const float *oldY2);

   FFTTOOLS_METHOD double *getCorrelation(int length,const double *oldY1,const  double *oldY2);



   FFTTOOLS_METHOD TGraph *correlateAndAverage(Int_t numGraphs, TGraph **grPtrPtr);


  FFTTOOLS_METHOD TGraph *correlateAndAveragePredeterminedZone(Int_t numGraphs, TGraph **grPtrPtr, Int_t *correlationBinPtr, Int_t binWiggleRoom);
  

   FFTTOOLS_METHOD TGraph *interpolateCorrelateAndAverage(Double_t deltaTInt,
                        Int_t numGraphs,
                        TGraph **grPtrPtr);
  

  FFTTOOLS_METHOD  Double_t *combineValuesUsingFFTs(Int_t numArrays, Double_t **thePtrPtr, Int_t eachLength);

  //Higher level functions that take and return TGraphs

  FFTTOOLS_METHOD  TGraph *makePowerSpectrum(const TGraph *grWave);

  FFTTOOLS_METHOD  TGraph *makePowerSpectrumPeriodogram(const TGraph *grWave);

  FFTTOOLS_METHOD  TGraph *makePowerSpectrumVoltsSeconds(const TGraph *grWave);

  FFTTOOLS_METHOD  TGraph *makePowerSpectrumVoltsSecondsBartlett(const TGraph *grWave);

  FFTTOOLS_METHOD  TGraph *makePSVSBartlettPaddedOverlap(const TGraph *grWave, Int_t padFactor=4, Int_t numFreqs=64);


  FFTTOOLS_METHOD  TGraph *makePowerSpectrumVoltsSecondsdB(const TGraph *grWave);

  FFTTOOLS_METHOD  TGraph *makePowerSpectrumMilliVoltsNanoSeconds(const TGraph *grWave);
  FFTTOOLS_METHOD  TGraph *makePowerSpectrumMilliVoltsNanoSecondsdB(const TGraph *grWave);

  FFTTOOLS_METHOD  TGraph *makePowerSpectrumVoltsSecondsPadded(const TGraph *grWave, Int_t padFactor=4);

  FFTTOOLS_METHOD  TGraph *makePowerSpectrumVoltsSecondsPaddeddB(const TGraph *grWave, Int_t padFactor=4);
  

  FFTTOOLS_METHOD  TGraph *makeRawPowerSpectrum(const TGraph *grWave);


  FFTTOOLS_METHOD  TGraph *getCorrelationGraph(const TGraph *gr1, const TGraph *gr2, Int_t *zeroOffset=0);




  FFTTOOLS_METHOD  TGraph *getNormalisedCorrelationGraph(const TGraph *gr1, const TGraph *gr2, Int_t *zeroOffset=0);


  FFTTOOLS_METHOD  TGraph *getNormalisedCorrelationGraphTimeDomain(const TGraph *gr1, const TGraph *gr2, Int_t *zeroOffset=0, Int_t useDtRange=0, Double_t dtMin=-1000, Double_t dtMax=1000);


  FFTTOOLS_METHOD  TGraph *getInterpolatedCorrelationGraph(const TGraph *grIn1, const TGraph *grIn2, Double_t deltaT);

  FFTTOOLS_METHOD  TGraph *makeInverseInverseSpectrum(const TGraph *grWave);
  

  FFTTOOLS_METHOD  TGraph *combineGraphsUsingFFTs(Int_t numGraphs, TGraph **grPtr,const double *theWeights=0);
  

   FFTTOOLS_METHOD TGraph *getBoxCar(const TGraph *grWave, Int_t halfWidth);

   FFTTOOLS_METHOD TGraph *getHilbertTransform(const TGraph *grWave);


   FFTTOOLS_METHOD double *  getHilbertTransform(int N, const double * y) ; 



   FFTTOOLS_METHOD TGraph *getHilbertEnvelope(const TGraph *grWave);

  //Utility functions (not necessarily FFT related but they'll live here for now


   FFTTOOLS_METHOD Double_t sumPower(const TGraph *gr,Int_t firstBin=-1,Int_t lastBin=-1);

  FFTTOOLS_METHOD  Double_t integratePower(const TGraph *gr,Int_t firstBin=-1,Int_t lastBin=-1);

  FFTTOOLS_METHOD  Double_t sumVoltageSquared(const TGraph *gr,Int_t firstBin,Int_t lastBin); 

  FFTTOOLS_METHOD  Double_t integrateVoltageSquared(const TGraph *gr,Int_t firstBin=-1,Int_t lastBin=-1);

  FFTTOOLS_METHOD  Int_t getPeakBin(const TGraph *gr); 

  FFTTOOLS_METHOD  Int_t getPeakBin(const TGraph *gr, Int_t firstBin, Int_t lastBin);  

  FFTTOOLS_METHOD Double_t getPeakXvalue(const TGraph *gr); 

  FFTTOOLS_METHOD  Double_t getPeakVal(const TGraph *gr, int *index=0);

  FFTTOOLS_METHOD Double_t getPeakVal(const TGraph *gr, Int_t firstBin, Int_t lastBin, int *index=0);

  FFTTOOLS_METHOD  Double_t getPeakSqVal(const TGraph *gr, int *index=0);

  FFTTOOLS_METHOD  void getPeakRmsSqVal(const TGraph *gr, Double_t &peak, Double_t &rms, Int_t *index=0);
  

  FFTTOOLS_METHOD  void getPeakRmsRectified(const TGraph *gr, Double_t &peak, Double_t &rms, Int_t *index=0);

  //Graph returning utility funcs

  FFTTOOLS_METHOD  TGraph *getSimplePowerEnvelopeGraph(const TGraph *gr);

  FFTTOOLS_METHOD  TGraph *smoothFFT(const TGraph *gr,Int_t factor) ;

  FFTTOOLS_METHOD  TGraph *subtractGraphs(const TGraph *grA, const TGraph *grB);

   FFTTOOLS_METHOD TGraph *divideGraphs(const TGraph *grA, const TGraph *grB);

  

  FFTTOOLS_METHOD  TGraph *translateGraph(const TGraph *grWave, const Double_t deltaT);


  FFTTOOLS_METHOD  TGraph *ratioSubtractOneGraphs(const TGraph *grA, const TGraph *grB) ;

  FFTTOOLS_METHOD  TGraph *dbGraphs(const TGraph *grA, const TGraph *grB);

  FFTTOOLS_METHOD  TGraph *padWave(const TGraph *grA, Int_t padFactor);

  FFTTOOLS_METHOD  TGraph *padWaveToLength(const TGraph *grA, Int_t newLength);

  FFTTOOLS_METHOD  TGraph *rectifyWave(const TGraph *gr, Int_t makeNeg=0);
  
  //Window functions

  FFTTOOLS_METHOD  Double_t  bartlettWindow(Int_t j, Int_t n);

  FFTTOOLS_METHOD  Double_t welchWindow(Int_t j, Int_t n);




  FFTTOOLS_METHOD  void takeDerivative(Int_t numPoints, const Double_t *inputX, const Double_t *inputY, Double_t *outputX, Double_t *outputY);


  FFTTOOLS_METHOD  TGraph *getDerviative(const TGraph *grIn);

  FFTTOOLS_METHOD  TGraph *getDerivative(const TGraph *grIn);


  FFTTOOLS_METHOD  TGraph *simplePassBandFilter(const TGraph *grWave, Double_t minFreq, Double_t maxFreq);


  FFTTOOLS_METHOD  TGraph *simpleNotchFilter(const TGraph *grWave, Double_t minFreq, Double_t maxFreq);



  FFTTOOLS_METHOD  TGraph *multipleSimpleNotchFilters(const TGraph *grWave, Int_t numNotches, const Double_t minFreq[], const Double_t maxFreq[]);


  FFTTOOLS_METHOD  TGraph *cropWave(const TGraph *grWave, Double_t minTime, Double_t maxTime);

  

  FFTTOOLS_METHOD  Double_t getWaveformSNR(const TGraph *gr); 

  FFTTOOLS_METHOD  Double_t getWaveformSNR(const TGraph *gr,Double_t &peakToPeak,Double_t &rms);


  FFTTOOLS_METHOD  Double_t getWaveformPeak(const TGraph *gr);
  FFTTOOLS_METHOD Double_t getEnvelopeSNR(const TGraph *gr);
  FFTTOOLS_METHOD Double_t getEnvelopeSNR(const TGraph *gr,Double_t &peakToPeak,Double_t &rms,Double_t &timeOfPeak);



   FFTTOOLS_METHOD Double_t simpleInterploate(Double_t x1, Double_t y1, Double_t x2, Double_t y2,Double_t x);
  
   FFTTOOLS_METHOD double * FFTCorrelation(int waveformlength, const FFTWComplex * A, const FFTWComplex * B, FFTWComplex * work = 0, 
                           int min_i = 0, int max_i =0, int order=1);  


   FFTTOOLS_METHOD void inPlaceShift(int N, double *x); 


   FFTTOOLS_METHOD void rotate(TGraph * g, int rot); 

   FFTTOOLS_METHOD void polySubtract(TGraph *g, int order=1); 


      
   enum DirectConvolveEdgeBehavior
   {
        ZEROES_OUTSIDE, 
        REPEAT_OUTSIDE
   }; 


    FFTTOOLS_METHOD double * directConvolve(int N, const double *x, int M, const double * h, double *y = 0, int delay = 0,  DirectConvolveEdgeBehavior edge_behavior = ZEROES_OUTSIDE); 


   FFTTOOLS_METHOD void wrap(size_t N, double * vals, double period = 360); 

   FFTTOOLS_METHOD void wrap(size_t N, double * vals, double period, double center); 

   FFTTOOLS_METHOD void wrap(size_t N, float * vals, float period, float center); 

   FFTTOOLS_METHOD void wrap(size_t N, float * vals, float period = 360); 


   FFTTOOLS_METHOD inline int fast_floor(double val) { return (int) val - (val < (int) val); }

   FFTTOOLS_METHOD inline double wrap(double val, double period, double center)
          { return val - period * fast_floor((val-center+period/2)/period); }

   FFTTOOLS_METHOD inline double wrap(double val, double period = 360)
          { return wrap(val,period,period/2); } 


   FFTTOOLS_METHOD void unwrap(size_t N, double * vals, double period = 360); 


   FFTTOOLS_METHOD void unwrap(size_t N, float * vals, float period = 360); 


   FFTTOOLS_METHOD double getDt(const TGraph * g, int realN = 0);  

   FFTTOOLS_METHOD double evalEvenGraph(const TGraph * g, double x); 



#ifndef FFTTOOLS_COMPAT_MODE

   void applyWindow(TGraph *g, const FFTWindowType *w); 
#endif

   
   FFTTOOLS_METHOD double randomRayleigh(double sigma=1, TRandom * rng = 0); 

   FFTTOOLS_METHOD double sinc(double x, double eps = 0); 
 
   FFTTOOLS_METHOD int loadWisdom(const char * file); 
   FFTTOOLS_METHOD int saveWisdom(const char * file); 

#ifndef FFTTOOLS_COMPAT_MODE

   TGraph * welchPeriodogram(const TGraph * gin, int segment_size, double overlap_fraction = 0.5, const FFTWindowType * window = &GAUSSIAN_WINDOW , bool truncate_extra = true, TGraph * gout = 0); 



   
   FFTTOOLS_METHOD double *  lombScarglePeriodogramSlow(int N, const double *x, const double * y, int nfreqs, const double * freqs, double * answer = 0); 

   FFTTOOLS_METHOD TGraph * lombScarglePeriodogram(const TGraph * g, double dt = 0, double oversample_factor  = 4 , 
                       double high_factor = 2, TGraph * replaceme = 0, int extirpolation_factor =4)  ; 
   FFTTOOLS_METHOD TGraph * lombScarglePeriodogram(int N, double dt, const double * __restrict x, 
                                   const double * __restrict y, double oversample_factor  = 4 , 
                                   double high_factor = 2, TGraph * replaceme = 0, int extirpolation_factor = 4)  ; 


#endif

//   TH2 * getPowerVsTimeUsingLombScargle(const TGraph * g, int nbins, double sampling_dt = 0, double oversample_factor = 4, double high_factor = 2, TH2 * useme = 0); 

   /* Compute Stokes parameters from hpol / vpol and their hilbert transforms 
    *
    * This computes either the average, as well as optionally the instantaneous values.  
    *
    *
    * @param N number of samples
    * @param hpol hpol waveform (evenly sampled)
    * @param hpol_hat hilbert transform of hpol waveform (evenly sampled)
    * @param vpol vpol waveform (evenly sampled)
    * @param vpol_hat hilbert transform of vpol waveform (evenly sampled)
    * @param Iavg pointer to where Stokes I integral average will be stored (or null to not store). 
    * @param Qavg pointer to where Stokes Q integral average will be stored (or null to not store). 
    * @param Uavg pointer to where Stokes U integral average will be stored (or null to not store).
    * @param Vavg pointer to where Stokes V integral average will be stored (or null to not store). 
    * @param Iins pointer to where Stokes I instantaneous value will be stored (or null to not store). Must be of length N unless only_store_max_instantaneous is true 
    * @param Qins pointer to where Stokes Q instantaneous value will be stored (or null to not store).  Must be of length N unless only_store_max_instantaneous is true
    * @param Uins pointer to where Stokes U instantaneous value will be stored (or null to not store). Must be of length N unless only_store_max_instantaneous is true
    * @param Vins pointer to where Stokes V instantaneous value will be stored (or null to not store).  Must be of length N unless only_store_max_instantaneous is true
    * @param only_max_instantaneous  if true,  only returns the maximum of the instantaneous quantities (at the maximum I in the range) 
    */ 
   FFTTOOLS_METHOD void stokesParameters(int N, const double *  hpol, const double *  hpol_hat, const double *  vpol, const double *  vpol_hat, 
                         double *  Iavg = 0, double *  Qavg = 0, double *  Uavg = 0, double *  Vavg = 0, 
                         double *  Iins = 0, double *  Qins = 0, double *  Uins = 0, double *  Vins = 0 , 
                         bool only_max_instantaneous = true ); 



   FFTTOOLS_METHOD void dftAtFreq(const TGraph * g, double freq, double * phase, double * amp = 0, double * real = 0, double * imag = 0); 


   FFTTOOLS_METHOD void dftAtFreqAndMultiples(const TGraph * g, double freq, int nmultiples, double * phase, double * amp = 0, double * real = 0, double * imag = 0); 


   FFTTOOLS_METHOD FFTWComplex * makeMinimumPhase(int N, const double * G,  double mindb = -100); 


   FFTTOOLS_METHOD TGraph * makeMinimumPhase(const TGraph *g, double mindb=-100); 


   FFTTOOLS_METHOD TGraph * makeFixedDelay(const TGraph *g, double delay = 1); 


   FFTTOOLS_METHOD double checkCausality(int N, const FFTWComplex * signal) ; 


   FFTTOOLS_METHOD double * rmsEnvelope(int N, double W, const double * x, const double * y, double * out = 0); 

   FFTTOOLS_METHOD double * peakEnvelope(int N, double min_distance, const double * x, const double *y, double * out = 0); 
#ifdef FFTTOOLS_COMPAT_MODE
};
#else
}
#endif 
   
#endif //FFTTOOLS_H