/Users/jmonk/Physics/ForIA/ForIA/AnalysisTools/Wavelet.hh

Go to the documentation of this file.

#ifndef FORIA_WAVELET_HH
#define FORIA_WAVELET_HH

#include "ForIA/IFourMomentum.hh"

#include "ForIA/MessageBox.hh"

#include <gsl/gsl_wavelet.h>

#include <cmath>
#include <vector>
#include <ostream>
//#include <iostream>
namespace ForIA{
 
  using std::vector;
  
  class Wavelet{
    
  public:
    
    Wavelet();
    
    ~Wavelet();
    class Coefficient{
      
    public:
      
      Coefficient(): m_doTranslation(true){}
      
      unsigned int scale()const{return m_scale;};
      
      double translation()const{
        if(m_doTranslation){
          m_translation = 0.5 * (translationMax() + translationMin());
          m_doTranslation = false;
        }
        return m_translation;
      }
      
      double translationMax()const{return m_upper;}
      
      double translationMin()const{return m_lower;}
      
      double magnitude()const{return m_magnitude;};
      
      unsigned int arrayPosition()const{return m_arrayPosition;};
      
      unsigned int level()const{return m_level;};
      
      unsigned int index()const{return m_index;};
      
      struct ByAmplitudeDown{
        bool operator()(const Coefficient &left, const Coefficient &right){
          return fabs(left.magnitude()) > fabs(right.magnitude());
        }
      };
      
    private:
      
      friend class Wavelet;
      
      unsigned int m_scale;
      
      mutable double m_translation;
      
      mutable bool m_doTranslation;
      
      double m_upper;
      
      double m_lower;
      
      double m_magnitude;
      
      unsigned int m_arrayPosition;
      
      unsigned int m_level;
      
      unsigned int m_index;
      
      
    };
    
    class FrequencyBand{
      
    public:
      
      FrequencyBand(): m_doEnergy(true), m_doVar(true), m_doAv(true){}
      
      double lower_frequency()const{return m_lowerFreq;};
      
      double upper_frequency()const{return m_upperFreq;};
      
      double physicalScale()const{return m_physScale;}
      
      unsigned int scale()const{return m_scale;}
      
      double variance()const{
        if(m_doVar){
          m_var=0.;
          for(const_iterator c=begin(); c!= end(); ++c){
            double y = c->magnitude() - average();
            m_var += y*y;
          }
          m_var=sqrt(m_var / (double)nCoefficients());
          m_doVar = false;
        }
        return m_var;
      }
      
      
      double average()const{
        if(m_doAv){
          m_av=0.;
          for(const_iterator c=begin(); c!= end(); ++c){
            m_av += c->magnitude();
          }
          m_av /= (double)nCoefficients();
          m_doAv = false;
        }
        return m_av;
      }
      
      double energy()const{
        if(m_doEnergy){
//          std::cout<<"========"<<std::endl;
          int p=1;
          p = p << level();
          m_inverse = vector<double>(p * nCoefficients(), 0.);
          
          for(const_iterator c=begin(); c != end(); ++c){
//            std::cout<<c->arrayPosition()<<std::endl;
            m_inverse[c->arrayPosition()] = c->magnitude();
          }
          
          m_parent->invert(m_inverse);
          m_energy=0.;
          for(vector<double>::const_iterator seg=m_inverse.begin();
              seg != m_inverse.end(); ++seg){
            m_energy += *seg;
          }
          
          m_doEnergy = false;
        }
        return m_energy;
      }
      
      unsigned int level()const{return m_level;};
      
      unsigned int nCoefficients()const{return m_coefficients.size();}
      
      const Coefficient &coefficient(unsigned int index)const{return m_coefficients[index];}
      
      typedef vector<Coefficient>::const_iterator const_iterator;
      
      typedef vector<Coefficient>::const_reverse_iterator const_reverse_iterator;
      
      const_iterator begin()const{return m_coefficients.begin();}
      
      const_reverse_iterator rbegin()const{return m_coefficients.rbegin();}
      
      const_iterator end()const{return m_coefficients.end();}
      
      const_reverse_iterator rend()const{return m_coefficients.rend();}
            
    private:
      
      friend class Wavelet;
      
      friend std::ostream &operator << (std::ostream &out, const FrequencyBand &freq);
      
      friend MessageBox &operator <<(MessageBox &box, const FrequencyBand &freq);
      
      Wavelet *m_parent;
      
      double m_lowerFreq;
      
      double m_upperFreq;
      
      double m_physScale;
      
      mutable double m_energy;
      
      mutable bool m_doEnergy;
      
      mutable double m_var;
      
      mutable bool m_doVar;
      
      mutable double m_av;
      
      mutable double m_doAv;
      
      unsigned int m_scale;
      
      unsigned int m_level;
      
      vector<Coefficient> m_coefficients;    
      
      mutable vector<double> m_inverse;
    };
    
    
    typedef vector<FrequencyBand>::const_iterator const_iterator;
        
    typedef vector<FrequencyBand>::const_reverse_iterator const_reverse_iterator;
    
    void setGrid(const vector<double> &grid);
    
    void setRangeStart(double min);
    
    void setRangeEnd(double max);
    
    unsigned int nSubBands()const;
    
    void invert(vector<double> &sequence)const;
    
    const_iterator begin()const{return m_bands.begin();}
    
    const_reverse_iterator rbegin()const{return m_bands.rbegin();}
    
    const_iterator end()const{return m_bands.end();}
    
    const_reverse_iterator rend()const{return m_bands.rend();}
    
    const FrequencyBand &frequencyBand(unsigned int level)const;
    
    const Coefficient &coefficient(unsigned int level, unsigned int index)const;
    
    double smoothingCoefficient()const{return m_smoothing;};
    
  private:
    
    void initialise();
    
    bool m_doInitialise;
    
    unsigned int m_nSegments;
    
    vector<double> m_grid;
    
    double m_min;
    
    double m_max;
    
    double m_range;
    
    double m_freqMax;
    
    double m_freqMin;
    
    bool m_haveSetRangeStart;
    
    bool m_haveSetRangeEnd;
    
    int m_nSubSamples;
        
    vector<FrequencyBand> m_bands;
    
    double m_smoothing;
    
    gsl_wavelet *m_wavelet;
    
    gsl_wavelet_workspace *m_workspace;
    
  };
}

#endif