TimeSeriesComponent.cpp

Go to the documentation of this file.

#include <algorithm>
#include <functional>
#include <boost/bind.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include "TimeSeriesComponent.h"
#include "FatalException.h"
#ifdef HAVEGSL
#include <gsl/gsl_spline.h>
#endif

namespace gplib
  {

    TimeSeriesComponent::TimeSeriesComponent()
      {
        samplerate = 0.0;
      }

    TimeSeriesComponent::~TimeSeriesComponent()
      {
      }

    TimeSeriesComponent::TimeSeriesComponent(const TimeSeriesComponent& source) :
      data(source.data), samplerate(source.samplerate), name(source.name),
          starttime(source.starttime)
      {
      }

    void TimeSeriesComponent::ShiftStart(const int npts)
      {
        //calculate the time_shift with microsecond precision
        boost::posix_time::time_duration time_shift =
            boost::posix_time::microseconds(boost::numeric_cast<unsigned int>(
                npts * 1000000 / samplerate));
        if (npts < 0) // if we want to remove something
          {
            if (static_cast<size_t> (abs(npts)) > data.size())
              throw FatalException("Trying to shift by too many points. ");
            // erase the data points
            data.erase(data.begin(), data.begin() + abs(npts));
            // correct the starttime with microsecond precision
            starttime = starttime + time_shift;
          }
        else //we want to add something
          {
            data.insert(data.begin(), npts, 0);
            starttime = starttime - time_shift;
          }
      }

    void TimeSeriesComponent::ShiftEnd(const int npts)
      {
        if (npts < 0)
          {
            if (static_cast<size_t> (abs(npts)) > data.size())
              throw FatalException("Trying to cut too many points. ");
            data.erase(data.end() + npts, data.end());
          }
        else
          {
            data.insert(data.end(), npts, 0);
          }
      }

    TimeSeriesComponent& TimeSeriesComponent::operator=(
        const TimeSeriesComponent& source)
      {
        if (this != &source)
          {
            this->samplerate = source.samplerate;
            this->name = source.name;
            this->starttime = source.starttime;
            this->data.assign(source.data.size(), 0); //assign and copy clears old data
            std::copy(source.data.begin(), source.data.end(),
                this->data.begin());
          }
        return *this;
      }

    //we want to make the dependency on gsl optional, so we only have this function
    //if gsl exists
#ifdef HAVEGSL
    void TimeSeriesComponent::Resample(const double newdt)
      {
        const size_t oldlength = data.size();
        double *oldtime = new double[oldlength];
        const double dt = GetDt();
        double start = 0.0;
        //we need time information for interpolation, the real starttime is irrelevant for this
        for (size_t i = 0; i < oldlength; ++i)
          {
            oldtime[i] = start + i * dt;
          }
        //initialize the gsl interpolation routines
        gsl_interp_accel *acc = gsl_interp_accel_alloc();
        gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, oldlength);
        gsl_spline_init(spline, oldtime, &data[0], oldlength);
        // calculate the lenght of the resampled timeseries
        const size_t newlength = boost::numeric_cast<size_t>(oldlength * (dt
            / newdt));
        // clear old data
        data.clear();
        //create new timeseries
        for (size_t i = 0; i < newlength; ++i)
          {
            data.push_back(gsl_spline_eval(spline, start + i * newdt, acc));
          }
        gsl_spline_free(spline);
        gsl_interp_accel_free(acc);
        samplerate = 1. / newdt;

      }
#endif //HAVEGSL

    TimeSeriesComponent& TimeSeriesComponent::operator*=(const double factor)
      {
        std::transform(data.begin(), data.end(), data.begin(), boost::bind(
            std::multiplies<double>(), factor, _1));
        return *this;
      }

    TimeSeriesComponent& TimeSeriesComponent::operator/=(const double numerator)
      {
        const double factor = 1. / numerator;
        operator*=(factor);
        return *this;
      }

    TimeSeriesComponent& TimeSeriesComponent::operator+=(const double shift)
      {
        std::transform(data.begin(), data.end(), data.begin(), boost::bind(
            std::plus<double>(), shift, _1));
        return *this;
      }

    TimeSeriesComponent& TimeSeriesComponent::operator+=(
        const TimeSeriesComponent &other)
      {
        std::transform(data.begin(), data.end(), other.data.begin(),
            data.begin(), std::plus<double>());
        return *this;
      }
    TimeSeriesComponent& TimeSeriesComponent::operator-=(const double shift)
      {
        operator+=(-shift);
        return *this;
      }
    //for this operator it is more efficient to reimplement it instead of calling += with -other
    TimeSeriesComponent& TimeSeriesComponent::operator-=(
        const TimeSeriesComponent &other)
      {
        std::transform(data.begin(), data.end(), other.data.begin(),
            data.begin(), std::minus<double>());
        return *this;
      }
  }