AnalysisWaveform.h

#ifndef ANALYSIS_WAVEFORM_HH
#define ANALYSIS_WAVEFORM_HH

class FFTWComplex; 
#include "TGraphAligned.h" 
#include "Math/Interpolator.h"
#include "FFTWindow.h" 

class AnalysisWaveform 
{


  public: 

    static void enableDebug(bool enable);  

    static void allowEvenToUnevenConversion(bool allow); 

    enum InterpolationType
    {
      AKIMA, 
      SPARSE_YEN,    
      REGULARIZED_SPARSE_YEN    
    }; 

    enum EvenEvaluationType
    {
      EVAL_LINEAR, 
      EVAL_AKIMA
    }; 

    static InterpolationType defaultInterpolationType; 

    struct InterpolationOptions 
    {

      InterpolationOptions() 
        : max_distance(8), regularization_order(0), mu(1e-3) {; } 

      int max_distance;  
      int regularization_order; 
      double mu; 
    }; 


    class PowerCalculationOptions
    {
      public: 
      PowerCalculationOptions() 
        : method(FFT), window_size(64)  { type = &FFTtools::GAUSSIAN_WINDOW; } 

      enum method
      {
        FFT, 
        BARTLETT,
        WELCH
      } method;  


      int window_size; 
      const FFTtools::FFTWindowType *type; 
    }; 

    static InterpolationOptions defaultInterpolationOptions; 

    /* Constructors */ 

    AnalysisWaveform(int Nt, const double * x, const double * y, double nominal_dt = 1./2.6, InterpolationType type = defaultInterpolationType, InterpolationOptions * opt  = &defaultInterpolationOptions);

    AnalysisWaveform(int Nt, const double * y, double dt, double t0);  

    static AnalysisWaveform * makeWf(const TGraph * g, bool even = true, AnalysisWaveform * replace = 0); 


    AnalysisWaveform(int Nt, const FFTWComplex * f, double df, double t0);  

    AnalysisWaveform(int Nt = 260, double dt=1./2.6, double t0=0);

    AnalysisWaveform(const AnalysisWaveform & other);  

    /*  Computes the RMS, either in time or frequency domain to avoid an FFT*/ 
    double getRMS() const;  

    AnalysisWaveform & operator=(const AnalysisWaveform & other); 


    static AnalysisWaveform * correlation(const AnalysisWaveform * A, const AnalysisWaveform * B, int npadfreq = 0, double scale =1 , int window_normalize = 0); 

    static AnalysisWaveform * convolution(const AnalysisWaveform * A, const AnalysisWaveform * B, int npadfreq = 0, double scale =1 ); 

    /* Enables or disables (check) and nagging on attemptoin correlation
     * without zero padding. Useful if you actually want to do a circular
     * correlation or if the check is too expensive for your tastes.
     *
     * Note that this only affects the current thread. 
     **/
    static void setCorrelationNag(bool nag); 


    AnalysisWaveform * autoCorrelation(int npadtime =1, int npadfreq = 0, double scale = 1); 

    bool checkIfPaddedInTime() const; 

    virtual ~AnalysisWaveform(); 

    /* Constant accessors. If you coerce the compiler into allowing modification, coupled waveforms won't be updated */ 

    const TGraphAligned * uneven() const ;

    const TGraphAligned * even()  const; 

    const FFTWComplex * freq() const; 

    const TGraphAligned * power() const; 

    const TGraphAligned * powerdB() const; 

    const TGraphAligned * phase() const; 

    const TGraphAligned * groupDelay() const; 

    const TGraphAligned * hilbertEnvelope() const; 

    const AnalysisWaveform * hilbertTransform() const; 

    int Nfreq() const;

    int Neven() const { return even()->GetN(); } 

    int Nuneven() const { return uneven()->GetN(); } 

    double deltaT() const { return dt ; } 

    double deltaF() const { return df ; }

    void forceEvenSize(int size); 

    // drawers since drawing is non-const (and we don't care about silly things like axes for constness)
    
    void drawEven(const char * opt = "", int color=-1) const; 

    void drawHilbertEnvelope(const char * opt = "", int color=-1) const; 

    void drawUneven(const char * opt = "", int color=-1) const; 

    void drawPower(const char * opt = "", int color=-1) const; 


    void drawPowerdB(const char * opt = "", int color=-1) const; 

    void drawPhase(const char * opt = "", int color=-1) const; 

    double evalEven(double t, EvenEvaluationType = EVAL_LINEAR) const; //TODO: add additional evaluation methods other than linear interpolation. indeed, would be best to eval multiple points at same time

    void evalEven(int N, const double * __restrict t , double * __restrict v, EvenEvaluationType = EVAL_LINEAR) const; 


    FFTWComplex * updateFreq(); 
    void updateFreq(int new_N, const FFTWComplex * new_freq, double new_df = 0) ;

    TGraphAligned * updateEven(); 

    void updateEven(const TGraph * replace_even); 

    TGraphAligned * updateUneven(); 

    void updateUneven(const TGraph * replace_uneven); 


    void padEven(int factor, int where = 1); 

    void padFreq(int factor); 


    void padFreqAdd(int npad);


    void setPowerCalculationOptions(PowerCalculationOptions & opt); 


    static void basisChange(AnalysisWaveform * __restrict hpol_or_lcp, AnalysisWaveform * __restrict vpol_or_rcp); 

    static void sumDifference(AnalysisWaveform * __restrict a, AnalysisWaveform * __restrict b); 


    void zeroMeanEven(); 

    void setTitle(const char * title); 
    /* Sets the line and marker color */ 
    void setColor(int c); 
    /* Sets the line width */ 
    void setWidth(int w);  

    /* Sets range user for x axis of frequency plots */ 
    void setFreqDisplayRange(double low, double high); 

    /* Sets range user for x axis of time plots */ 
    void setTimeDisplayRange(double low, double high); 

  private: 
    void calculateEvenFromUneven() const; 
    void calculateEvenFromFreq() const;  
    void calculateFreqFromEven() const;  
    void calculateUnevenFromEven() const;

    /* Storage vars */ 
    mutable TGraphAligned g_uneven; 
    mutable TGraphAligned g_hilbert_envelope; 
    mutable TGraphAligned g_even; 
    mutable TGraphAligned g_power; 
    mutable TGraphAligned g_power_db;
    mutable TGraphAligned g_phase;  
    mutable TGraphAligned g_group_delay; 


    mutable double dt;  
    mutable double df; 
    mutable int fft_len; 
    mutable FFTWComplex * fft; 

    mutable ROOT::Math::Interpolator theEvenAkimaInterpolator; 
    const ROOT::Math::Interpolator * evenAkimaInterpolator() const; 

    InterpolationType interpolation_type; 
    InterpolationOptions interpolation_options; 
    PowerCalculationOptions power_options; 

    /* State vars */ 
    mutable bool just_padded;  
    mutable bool must_update_uneven; 
    mutable bool must_update_freq; 
    mutable bool must_update_even; 
    mutable bool uneven_equals_even; 
    mutable bool power_dirty; 
    mutable bool power_db_dirty; 
    mutable bool phase_dirty; 
    mutable bool group_delay_dirty;
    mutable bool hilbert_envelope_dirty; 
    mutable bool hilbert_dirty; 
    mutable bool even_akima_interpolator_dirty; 
    mutable AnalysisWaveform * hilbert_transform;
    mutable int force_even_size; 
    unsigned uid; 
    void nameGraphs(); 
}; 



#endif