s2_corr_inline.hpp

/* s2_corr_inline.hpp */
/* $Id: s2_corr_inline.hpp 22125 2013-05-08 16:35:52Z urij $ */
#ifndef __s2_corr_inline_hpp
#define __s2_corr_inline_hpp

#ifdef GE_BUILD
#include <s2_corr.hpp>
#else
#include <s2proc/s2_corr.hpp>
#endif

#include <math.h>

template <template <class> class Functor>
bool S2Corr::traceBaseOp(Trace &tr, const DVector &bases)
{
  Functor<Trace::AmplT> functor;
  
  // number of bases should be > 0 and < than trace lenght
  if (bases.size() < 1 || tr.size() < 1 || bases.size() > tr.size()) 
    {
      printf ("S2Corr->traceBaseOp: error: %zu bases provided for %zu samplse in trace.\n", 
          bases.size(), tr.size());
      return false;
    }
  
  AvgBaseSep<Trace::AmplT> sep(&tr[0], tr.size(), bases.size());
  size_t bn = 0;
  Trace::AmplT res = 0.;
  
  // apply weights for trace
  for (size_t i=0; i<bases.size(); i++)
    {
      Trace::AmplT *trb = sep.get(i, bn); 
      for (size_t j=0; j<bn; j++)
    {
      res = functor(trb[j], bases[i]);
#if 0
      printf ("S2Corr->traceBaseOp: base %d, sample %d : %g op %g = %g\n", 
          i, j, trb[j], bases[i], res);
#endif  
      trb[j] = res;
    }
    }
  
  return true;
}

class CDFT2D;

template <typename FFT_T>
bool S2Corr::tracesCCF(FFT_T &fft, 
               const Trace &tr1, 
               const Trace &tr2, 
               Trace &ccf,
               bool norma)
{
  if(tr1.size() != tr2.size() || tr1.size() == 0 || tr2.size() == 0)
    {
      printf("S2Corr::tracesCCF: error: incorrect size of traces tr1.size() = %zu tr2.size = %zu\n", 
         tr1.size(), tr2.size()); 
      return false;
    }
      
  double En1 = 0.;
  double En2 = 0.;
  for(size_t n=0; n<tr1.size(); ++n)
    {
      En1 += tr1[n]*tr1[n];
      En2 += tr2[n]*tr2[n];
    }
  
  float E = 1e-7;
  double norm_coef = 0;
  
  if(En1 > E && En2 > E)
    norm_coef = 1/sqrt(En1*En2);
  
  fft.setRotateFlag(false);
  fft.setDirection(1);

  Trace sp1_re(tr1.header(), tr1.begin(), tr1.end());
  Trace sp1_im(tr1.header(), tr1.size(), 0.);
  fft.calc(sp1_re, sp1_im);
    
  Trace sp2_re(tr2.header(), tr2.begin(), tr2.end());
  Trace sp2_im(tr2.header(), tr2.size(), 0.);
  fft.calc(sp2_re, sp2_im);
    
  // get spectr of CCF
  Trace resultRe(sp1_re.header(), sp1_re.size(), 0.);
  Trace resultIm(sp1_im.header(), sp1_im.size(), 0.);
  
  for(unsigned int i=0; i<sp1_re.size(); ++i)
    {
      //multiplying (a + ib) * (c - id) = (a*c + d*b) + i(b*c - a*d)
      // a = sp1_re[i]
      // b = sp1_im[i]
      // c = sp2_re[i]
      // d = sp2_im[i]

      resultRe[i] = sp1_re[i]*sp2_re[i] + sp1_im[i]*sp2_im[i];
      resultIm[i] = sp1_im[i]*sp2_re[i] - sp1_re[i]*sp2_im[i];
    }
    
  fft.setDirection(-1); // invers fft

  if(!fft.calc(resultRe, resultIm))
    return false;
  
  size_t half = resultRe.size()/2;
  size_t add = resultRe.size()%2;
  
  //printf("tracesCCF::INFO: result t0 = %g dt = %g \n", resultRe.header().t0, resultRe.header().timestep);
  ccf = Trace(tr1.header(), resultRe.begin(), resultRe.end());
  
  for (size_t i=half+add; i<ccf.size(); i++)
    ccf[i-half-add] = resultRe[i];
  
  for (size_t i=0; i<half+add; i++)
    ccf[i+half] = resultRe[i];
  
  if(norma)
    for (size_t i=0; i<ccf.size(); i++)
      ccf[i]*=norm_coef;
  
  // set proper zero time
  ccf.header().t0 = (-ccf.header().timestep)*(half + add);
  
  return true;
}

template <typename FFT_T>
bool S2Corr::traceACF(FFT_T &fft, 
              const Trace &tr, 
              Trace &acf, 
              bool norma)
{

  if(0 == tr.size())
    {printf("traceACF::ERROR: Empty traces tr.size() = %zu\n", tr.size()); return false;}
  
  fft.setRotateFlag(false);
  fft.setDirection(1);
  
  Trace sp_re = tr;
  Trace sp_im = tr;
  std::fill (sp_im.begin(), sp_im.end(), 0);
  fft.calc(sp_re, sp_im);
     
  // get spectr of ACF
  Trace resultRe = sp_re, resultIm = sp_im;
  
  std::fill (resultRe.begin(), resultRe.end(), 0);
  std::fill (resultIm.begin(), resultIm.end(), 0);

  for(unsigned int i=0; i<sp_re.size(); ++i)
    {
      //multiplying (a + ib) * (a - ib) = (a*a + b*b) + i(b*a - a*b)
      // a = sp_re[i]
      // b = sp_im[i]
     
      resultRe[i] = sp_re[i]*sp_re[i] + sp_im[i]*sp_im[i];
      resultIm[i] = sp_im[i]*sp_re[i] - sp_re[i]*sp_im[i];
    }
    
  fft.setDirection(-1); // invers fft
  
  if(!fft.calc(resultRe, resultIm))
    {
      return false;
    }
  
  size_t half = resultRe.size()/2;
  size_t add = resultRe.size()%2;
        
  //for (size_t i=0; i<half; i++)
  //std::swap<Trace::AmplT>(resultRe[i], resultRe[half + i]);

  //Trace::AmplT first_sample = resultRe[0];
    
  acf = resultRe;
  
  for (size_t i=half+add; i<acf.size(); i++)
    acf[i-half-add] = resultRe[i];
  
  for (size_t i=0; i<half+add; i++)
    acf[i+half] = resultRe[i];
  
  if(norma)
    {
      float norm_coef;
      if(acf[half+add] > 1e-7)
    norm_coef = 1./acf[half+add];
      else
    norm_coef = 0.;
      for (size_t i=0; i<acf.size(); i++)
    acf[i]*=norm_coef;
    }
  
  // set proper zero time
  acf.header().t0 = (-acf.header().timestep)*(half + add);

  return true;
}


template <typename FFT_T>
bool S2Corr::SpectrWidth(FFT_T &fft, 
             const Trace &tr, 
             float &sw_val)
{
  if(0 == tr.size())
    {printf("SpectrWidth::ERROR: Empty traces tr.size() = %zu\n", tr.size()); return false;}
  
  float dt = tr.header().timestep;
  
  fft.setRotateFlag(false);
  fft.setDirection(1);
  
  Trace sp_re = tr;
  Trace sp_im = tr;
  std::fill (sp_im.begin(), sp_im.end(), 0);
  fft.calc(sp_re, sp_im);
     
  // get spectr of trace
  Trace resultRe = sp_re, resultIm = sp_im; // for correct .header() values
  
  float width_val = 0;
  float norm_val = 0;

  std::fill (resultRe.begin(), resultRe.end(), 0);
  std::fill (resultIm.begin(), resultIm.end(), 0);

  for(unsigned int i=0; i<sp_re.size(); ++i)
    {
      //multiplying (a + ib) * (a - ib) = (a*a + b*b) + i(b*a - a*b) = a*a + b*b
      // a = sp_re[i]
      // b = sp_im[i]
      resultRe[i] = sp_re[i]*sp_re[i] + sp_im[i]*sp_im[i];
      
      width_val += (i+1)*(i+1)*resultRe[i];
      norm_val += resultRe[i];
    }
  if(0. == norm_val)
    {
      sw_val = 0.;
      return true;
    }
  
  sw_val = width_val/norm_val;
  sw_val /= 4.*dt*dt;
  
  return true;
}

#endif /* s2_corr_inline.hpp */