RadioTraceUtilities.cc

Go to the documentation of this file.

#include "RadioTraceUtilities.h"

#include <utl/TraceAlgorithm.h>
#include <utl/Math.h>


namespace utl {


  void
  RadioTraceUtilities::Noisefinder(const utl::Trace<double>& channeltrace,
                                   double& noiseRMS, double& noiseMean,
                                   double noiseWindowStart, double noiseWindowStop)
  {
    const unsigned int start = noiseWindowStart / channeltrace.GetBinning();
    const unsigned int stop = noiseWindowStop / channeltrace.GetBinning();
    noiseRMS = TraceAlgorithm::RootMeanSquare(channeltrace, start, stop);
    noiseMean = TraceAlgorithm::Mean(channeltrace, start, stop);
  }


  void
  RadioTraceUtilities::Pulsefinder(const utl::TraceV3D& stationtrace,
                                   double& peakAmplitude, double& peakTime, double& peakTimeError,
                                   const double signalSearchWindowStart, const double signalSearchWindowStop,
                                   unsigned int& sample)
  {
    peakAmplitude = 0;
    const unsigned int start = signalSearchWindowStart / stationtrace.GetBinning();
    const unsigned int stop = signalSearchWindowStop / stationtrace.GetBinning();

    for (unsigned int i = start; i < stop; ++i) {
      if (stationtrace[i].GetMag() > peakAmplitude) {
        peakAmplitude = stationtrace[i].GetMag();
        peakTime = i * stationtrace.GetBinning();
        sample = i;
      }
    }
    /* Here only the uncertainty due to the finite binnig is considered,
     * the uncertainty due to noise and from the GPS system is added later. */
    peakTimeError = stationtrace.GetBinning() / std::sqrt(12);
  }


  void
  RadioTraceUtilities::Pulsefinder(const utl::Trace<double>& channeltrace,
                                   double& peakAmplitude, double& peakTime, double& peakTimeError,
                                   const double signalSearchWindowStart, const double signalSearchWindowStop,
                                   unsigned int& sample)
  {
    peakAmplitude = 0;
    const unsigned int start = signalSearchWindowStart / channeltrace.GetBinning();
    const unsigned int stop = signalSearchWindowStop / channeltrace.GetBinning();

    // if trace is an envelope, the abs is not necessary but also does not harm
    for (unsigned int i = start; i < stop; ++i) {
      if (std::abs(channeltrace[i]) > peakAmplitude) {
        peakAmplitude = std::abs(channeltrace[i]);
        peakTime = i * channeltrace.GetBinning();
        sample = i;
      }
    }
    /* Here only the uncertainty due to the finite binnig is considered,
     * the uncertainty due to noise and from the GPS system is added later. */
    peakTimeError = channeltrace.GetBinning() / std::sqrt(12);
  }


  void
  RadioTraceUtilities::PulseFWHMIntegrator(const utl::Trace<double>& channeltrace,
                                           unsigned int sample, const double peakAmplitude,
                                           double& signalTimeFWHM, double& integratedSignal,
                                           double& signalWindowStart, double& signalWindowStop)
  {
    signalTimeFWHM = 0;
    integratedSignal = 0;
    unsigned int tmp = sample;

    while (sample < channeltrace.GetSize() && channeltrace[sample] > 0.5 * peakAmplitude) {
      ++signalTimeFWHM;
      integratedSignal += channeltrace[sample];
      ++sample;
    }
    signalWindowStop = sample * channeltrace.GetBinning();
    sample = tmp;

    while (sample > 0 && channeltrace[sample] > 0.5 * peakAmplitude) {
      ++signalTimeFWHM;
      integratedSignal += channeltrace[sample];
      --sample;
    }
    signalWindowStart = sample * channeltrace.GetBinning();
    signalTimeFWHM = signalTimeFWHM * channeltrace.GetBinning();
  }


  void
  RadioTraceUtilities::PulseFixedWindowIntegrator(const utl::Trace<double>& channeltrace,
                                                  const unsigned int sample, const double integrationTime,
                                                  double& integratedSignal, double& signalWindowStart,
                                                  double& signalWindowStop, const bool usePower)
  {
    integratedSignal = 0;
    signalWindowStart = sample * channeltrace.GetBinning() - integrationTime / 2;
    signalWindowStop = sample * channeltrace.GetBinning() + integrationTime / 2;
    const unsigned int start = std::max(0.0, sample - integrationTime / 2 / channeltrace.GetBinning());
    const unsigned int stop =
      std::min(channeltrace.GetSize() * 1., sample + integrationTime / 2 / channeltrace.GetBinning());

    for (unsigned int i = start; i <= stop; ++i) {
      if (usePower) {
        integratedSignal += Sqr(channeltrace[i]);
      } else {
        integratedSignal += std::abs(channeltrace[i]);
      }
    }
    integratedSignal *= channeltrace.GetBinning();
  }


  void
  RadioTraceUtilities::PulseSlidingWindowIntegrator(const utl::Trace<double>& channeltrace,
                                                    double integrationTime, double& integratedSignal,
                                                    double& signalWindowStart, double& signalWindowStop,
                                                    double signalSearchWindowStart,
                                                    double signalSearchWindowStop,
                                                    int& startSample, int& stopSample)
  {
    const double binning = channeltrace.GetBinning();
    integratedSignal = 0;
    int count = 0;

    for (int windowStartPos = signalSearchWindowStart / binning;
         windowStartPos <= (signalSearchWindowStop - integrationTime) / binning; ++windowStartPos) {
      double newIntegratedSignal = 0;
      count = 0;
      for (int integratePos = 0; integratePos < integrationTime / binning; ++integratePos) {
        newIntegratedSignal += Sqr(channeltrace[windowStartPos + integratePos]);
        ++count;
      }
      newIntegratedSignal /= count;
      if (newIntegratedSignal > integratedSignal) {
        integratedSignal = newIntegratedSignal;
        signalWindowStart = windowStartPos * binning;
        startSample = windowStartPos;
        stopSample = windowStartPos + integrationTime / binning;
      }
    }

    signalWindowStop = signalWindowStart + integrationTime;
  }

}