RecCalc.cpp

Go to the documentation of this file.

#include "RecCalc.h"
#include "types.h"
#include "convert.h"
#include "miscfunc.h"
#include "IterDecon.h"
#include <iostream>
#include <fstream>
#include <algorithm>
#include <boost/bind.hpp>
#include <boost/cast.hpp>
#include <boost/filesystem.hpp>
#include <numeric>
#include <functional>
#include <cassert>
#include <gsl/gsl_math.h>
#include <sys/types.h>
#include <unistd.h>
#include "statutils.h"
#include "WFunc.h"
#include "CRFSpecElement.h"
#include "SeisTools.h"
#include "Adaptors.h"

using namespace std;

RecCalc::RecCalc(const int myshift, const double myomega, const double mysigma,
    const double myc, bool multicalc, const trfmethod themethod) :
  method(themethod), Spectrum(multicalc), c(myc), omega(myomega),
      sigma(mysigma), shift(myshift), normalize(false)
  {

  }

RecCalc::RecCalc(const RecCalc &Old) :
  method(Old.method), Spectrum(Old.Spectrum), c(Old.c), omega(Old.omega),
      sigma(Old.sigma), shift(Old.shift), normalize(Old.normalize),
      RadComp(Old.RadComp), VerComp(Old.VerComp)
  {

  }

RecCalc& RecCalc::operator=(const RecCalc& source)
  {
    if (this == &source) return *this;
    method = source.method;
    Spectrum = source.Spectrum;
    c = source.c;
    omega = source.omega;
    sigma = source.sigma;
    shift = source.shift;
    normalize = source.normalize;
    RadComp = source.RadComp;
    VerComp = source.VerComp;
    return *this;
  }

RecCalc::~RecCalc()
  {
  }

void RecCalc::SpectralDivision(const SeismicDataComp &RComp,
    const SeismicDataComp &VComp, SeismicDataComp &Receiver)
  {
    const double omegastep = 2. * PI/(RComp.GetDt() * RComp.GetData().size());
    // Make local copies for trimming an mean removal
    ttsdata LocalRad(RComp.GetData());
    ttsdata LocalVer(VComp.GetData());
    // Make both components the same length 
    if (LocalRad.size() > LocalVer.size())
      LocalRad.erase(LocalRad.end() - (LocalRad.size() - LocalVer.size()),
          LocalRad.end());
    else if (LocalVer.size() > LocalRad.size())
      LocalVer.erase(LocalVer.end() - (LocalVer.size() - LocalRad.size()),
          LocalVer.end());
    //remove the mean
    SubMean(LocalVer.begin(), LocalVer.end());
    SubMean(LocalRad.begin(), LocalRad.end());
    //apply a window function for spectral calculation
    ApplyWindow(LocalVer.begin(), LocalVer.end(), LocalVer.begin(), Steep());
    ApplyWindow(LocalRad.begin(), LocalRad.end(), LocalRad.begin(), Steep());
    //assign space for spectra
    tcompdata RadSpec(LocalRad.size()/2+1);
    tcompdata VerSpec(LocalVer.size()/2+1);
    //calculate spetrca
    Spectrum.CalcSpectrum(LocalRad.begin(), LocalRad.end(), RadSpec.begin(),
        RadSpec.end());
    Spectrum.CalcSpectrum(LocalVer.begin(), LocalVer.end(), VerSpec.begin(),
        VerSpec.end());
    tcompdata RecSpec(RadSpec.size(), 0);

    typedef tcompdata::iterator tcomit;

    Receiver.GetData().assign(LocalRad.size(), 0);
    // determine the waterlevel threshold from the maximum vertical spectrum and the waterlevel parameter c
    tcomp cmaxelem = *max_element(VerSpec.begin(), VerSpec.end(),
        gplib::absLess<tcomp,tcomp>());
    cmaxelem *= conj(cmaxelem);
    cmaxelem *= c;
    // the denominator for the gaussian filter
    const double denom = 1./(4.0 * sigma*sigma);
    //calculate the spectrum of the receiver function
    transform(RadSpec.begin(), RadSpec.end(), VerSpec.begin(), RecSpec.begin(),
    CalcSpectralElement(cmaxelem));
    //determine the phaseshift for the desired shift in points
    const dcomp phaseshift(0, shift);
    // apply the gaussian filter
    for (tcomit recit = RecSpec.begin(); recit < RecSpec.end(); ++recit)
      {
        double freq = omegastep * distance(RecSpec.begin(), recit);
        *recit *= exp(- (freq*freq)*denom);
        *recit *= exp( -phaseshift * freq);
      }
    //transform receiver function spectrum back to time domain
    Spectrum.CalcTimeSeries(RecSpec.begin(), RecSpec.end(), Receiver.GetData().begin(), Receiver.GetData().end());
  }

void RecCalc::IterativeDeconvolution(const SeismicDataComp &RComp,
    const SeismicDataComp &VComp, SeismicDataComp &Receiver)
  {
    gplib::rvec Filter(RComp.GetData().size());
    const double dt = RComp.GetDt();
    double normfactor = 1.; // /(sigma * sqrt(2.* acos(-1.)));
    const int shiftpoints = boost::numeric_cast<int>(shift/dt);
    const double variancesq = sigma*dt * sigma*dt;
    const unsigned int maxit = 100;
    const double minimprove = 1e-4; // we want to stop when the improvement is less than 0.1%
    const int maxpoints = Filter.size()/2; //we store it in a signed variable to avoid compiler warning for the loop below 
    //the loop index HAS TO BE SIGNED ! 
    for (int i = 0; i <= maxpoints; ++i)
      {
        Filter(i) = normfactor * exp( -pow((i-shiftpoints),2.) * variancesq);
        Filter(Filter.size() -i-1) = normfactor * exp( -pow((-i-1-shiftpoints),2.) * variancesq);
      }

    gplib::rvec RadFilter(RComp.GetData().size()), VerFilter(VComp.GetData().size());
    copy(RComp.GetData().begin(),RComp.GetData().end(),RadFilter.begin());
    copy(VComp.GetData().begin(),VComp.GetData().end(),VerFilter.begin());
    Convolve(RadFilter,Filter,RadFilter);
    Convolve(VerFilter,Filter,VerFilter);
    gplib::rvec Current(RadFilter), Output(RadFilter);

    IterDecon Decon(VerFilter.size());
    const double power = ublas::prec_inner_prod(RadFilter,RadFilter);
    double error = power;
    double lasterror = 1.;
    double improvement = 1.;
    unsigned int iter = 0;
    while (improvement> minimprove && iter < maxit)
      {
        Decon.AdaptFilter(VerFilter,Current);
        Decon.CalcOutput(VerFilter,Output);

        Current = RadFilter-Output;
        error = ublas::prec_inner_prod(Current,Current) /power;
        improvement = lasterror - error;
        lasterror = error;
        ++iter;
      }
    Convolve(Decon.GetWeights(),Filter,Current); //we re-use current to store the filtered weights, to save another object
    Receiver.GetData().assign(Current.size(),0);
    copy(Current.begin(),Current.end(),Receiver.GetData().begin());
  }

void RecCalc::CalcRecData(const SeismicDataComp &RadComp,
    const SeismicDataComp &VerComp, SeismicDataComp &Receiver)
  {
    Receiver.CopyHeader(RadComp);
    switch (method)
      {
    case iterdecon:
      IterativeDeconvolution(RadComp, VerComp, Receiver);
      break;
    case specdiv:
      SpectralDivision(RadComp, VerComp, Receiver);
      break;
      }
    // set parameters for the timeseries object
    Receiver.SetB(-shift);
    Receiver.SetDt(RadComp.GetDt());
    if (normalize)
      {
        Normalize(Receiver.GetData());
      }
  }

void RecCalc::CalcRecSynth(const std::string &filename, ResPkModel &Model,
    SeismicDataComp &Receiver, const bool cleanfiles)
  {
    SynthPreParallel(filename, Model, Receiver, cleanfiles);
    SynthSafeParallel(filename, Model, Receiver, cleanfiles);
    SynthPostParallel(filename, Model, Receiver, cleanfiles);
  }

void RecCalc::SynthPreParallel(const std::string &filename, ResPkModel &Model,
    SeismicDataComp &Receiver, const bool cleanfiles)
  {
    Model.WriteRunFile(filename);
    Model.WriteData(filename+".mod");
  }
void RecCalc::SynthSafeParallel(const std::string &filename, ResPkModel &Model,
    SeismicDataComp &Receiver, const bool cleanfiles)
  {
    system((filename).c_str());
    //system(("./iter"+filename).c_str());

    RadComp.GetData(filename+".mod_sp.r");
    VerComp.GetData(filename+".mod_sp.z");
    Receiver = RadComp;
    CalcRecData(RadComp, VerComp, Receiver);
  }

void RecCalc::SynthPostParallel(const std::string &filename, ResPkModel &Model,
    SeismicDataComp &Receiver, const bool cleanfiles)
  {
    if (cleanfiles)
      {
        system(("rm "+filename+"*").c_str());
      }
  }

---