/Users/jmonk/Physics/ForIA/src/AnalysisTools/Wavelet.cxx

Go to the documentation of this file.

#include "ForIA/AnalysisTools/Wavelet.hh"

#include "boost/lexical_cast.hpp"

#include <cmath>

#include <stdexcept>
//#include <iostream>

namespace ForIA{
  Wavelet::Wavelet(): m_doInitialise(true), m_haveSetRangeStart(false), m_haveSetRangeEnd(false),
  m_nSubSamples(0), m_wavelet(0), m_workspace(0){
    
  }
  
  Wavelet::~Wavelet(){
    if(m_workspace) gsl_wavelet_workspace_free(m_workspace);
    if(m_wavelet) gsl_wavelet_free(m_wavelet);
  }
  
  void Wavelet::initialise(){
    
    m_nSegments = m_grid.size();
    
    if(!m_haveSetRangeStart || !m_haveSetRangeEnd){
      throw std::runtime_error("Wavelet: Cannot initialise without first setting co-ordinate range of input sequence");
    }
    
    if( ((m_nSegments-1) & m_nSegments)){ 
      throw std::runtime_error("Number of segments must be radix 2.  Number given is " + 
                               boost::lexical_cast<std::string>(m_nSegments));
      return;
    }
  
    m_range = m_max - m_min;
    
    m_freqMax = 0.5 * (double)m_nSegments / m_range;
    
    m_freqMin = 0.5 / m_range;
    
    m_nSubSamples = 0;
    
    int tmp = m_nSegments;
    
    while(tmp >>= 1) ++m_nSubSamples;
    
//    m_nSubSamples = log2(m_nSegments);
        
    m_workspace = gsl_wavelet_workspace_alloc(m_nSegments);
    
    m_wavelet = gsl_wavelet_alloc(gsl_wavelet_daubechies_centered, 4);
  
    m_bands.clear();
    
    unsigned int kMax=1;
    
    unsigned int index=1;
    
    double freqMin = m_freqMin;
    double freqMax = freqMin * 2.;
    
    unsigned int scale = m_nSegments;
    
    double physScale = m_range;
    
    for(unsigned int j=0; j!= nSubBands(); ++j){
      
      m_bands.push_back(FrequencyBand());
      
      m_bands.back().m_parent = this;
      m_bands.back().m_level = nSubBands() - j;
      m_bands.back().m_lowerFreq = freqMin;
      m_bands.back().m_upperFreq = freqMax;
      m_bands.back().m_physScale = physScale;
      
      m_bands.back().m_scale = scale;
      
      double tau_min = m_min;
      double tau_incr = m_range / (double)kMax;
      double tau_max = tau_min + tau_incr;
      
      for(unsigned int k=0; k != kMax; ++k){
        
        m_bands.back().m_coefficients.push_back(Coefficient());
        m_bands.back().m_coefficients.back().m_arrayPosition=index;
        m_bands.back().m_coefficients.back().m_level = m_bands.back().level();
        m_bands.back().m_coefficients.back().m_index = k;
        m_bands.back().m_coefficients.back().m_upper = tau_max;
        m_bands.back().m_coefficients.back().m_lower = tau_min;
        
        m_bands.back().m_coefficients.back().m_scale = scale;
        
        tau_min = tau_max;
        tau_max += tau_incr;
        
        ++index;
      }
      
      freqMin = freqMax;
      freqMax *= 2.;
      kMax *= 2;
      scale /= 2;
      physScale *= 0.5;
    }    
    
    m_doInitialise = false;
    
    return;
  }
  
  void Wavelet::setRangeStart(double min){
    if(!m_doInitialise){
      throw std::runtime_error("Wavelet: Cannot set start range of co-ordinate after initialisation!");
    }
    
    m_min = min;
    
    m_haveSetRangeStart = true;
    return;
  }
  
  void Wavelet::setRangeEnd(double max){
    if(!m_doInitialise){
      throw std::runtime_error("Wavelet: Cannot set end range of co-ordinate after initialisation!");
    }
    
    m_max = max;
    
    m_haveSetRangeEnd = true;
    return;
  }
  
  unsigned int Wavelet::nSubBands()const{
    return m_nSubSamples;
  }
  
  const Wavelet::FrequencyBand &Wavelet::frequencyBand(unsigned int level)const{
    unsigned int index = nSubBands() - level;
    return m_bands[index];
  }
  
  const Wavelet::Coefficient &Wavelet::coefficient(unsigned int level, unsigned int index)const{
    return frequencyBand(level).coefficient(index);
  }
  
  void Wavelet::setGrid(const vector<double> &grid){
        
    m_grid = grid;
    
    if(m_doInitialise) initialise();
    
    if(m_grid.size() != m_nSegments){
      throw std::runtime_error("Number of segments has changed since initialisation! Was first called with " + 
                               boost::lexical_cast<std::string>(m_nSegments) + 
                               " segments, now being called with " + boost::lexical_cast<std::string>(m_grid.size()));
    }
    
    vector<double> output(grid.begin(), grid.end());
        
    int status = gsl_wavelet_transform_forward(m_wavelet, &(output[0]), 1, m_nSegments, m_workspace);
    
    if(status != GSL_SUCCESS){
      throw std::runtime_error("Wavelet transform failed due to either lack of memory or data that is not radix 2!");
    }
            
    m_smoothing = output[0];
    
    unsigned int index=1;
    
    for(vector<FrequencyBand>::iterator band = m_bands.begin();
        band != m_bands.end(); ++band){
    
      band->m_doEnergy = true;
      band->m_doVar = true;
      band->m_doAv = true;
      for(vector<Coefficient>::iterator coeff = band->m_coefficients.begin();
          coeff != band->m_coefficients.end(); ++coeff){
                
        coeff->m_magnitude = output[index];

        ++index;
      }
             
    }
    
    return;
  }
  
  void Wavelet::invert(vector<double> &sequence)const{
    
    if(sequence.size() != m_nSegments){
      throw std::runtime_error("Number of segments in sequence to be inverted (" + 
                               boost::lexical_cast<string>(sequence.size()) + ") does not match wavelet (" + 
                               boost::lexical_cast<string>(m_nSegments) + ")");
    }
        
    gsl_wavelet_transform_inverse(m_wavelet, &(sequence[0]), 1, m_nSegments, m_workspace);
    
    return;
  }
  
  MessageBox &operator << (MessageBox &box, const Wavelet::FrequencyBand &freq){
    box.setTitle("Wavelet Frequency Band");
    
    box.addEntry("Level", freq.level());
    box.addEntry("Scale factor", freq.scale());
    box.addEntry("Physical scale", freq.physicalScale());
    box.addEntry("Number of Coefficients", freq.nCoefficients());
    box.addEntry("Frequency range", "{ " + boost::lexical_cast<string>(freq.lower_frequency()) + ", " + 
                 boost::lexical_cast<string>(freq.upper_frequency()) + "}");
    
    return box;
  }
  
  std::ostream &operator << (std::ostream &out, const Wavelet::FrequencyBand &freq){
    MessageBox tmp;
    
    tmp<<freq;
    out<<tmp;
    
    return out;
  }
  
}