amplitudes.h

#ifndef _AMPLITUDES_H_
#define _AMPLITUDES_H_
#include <sstream>
#include "mspass/utility/MsPASSError.h"
#include "mspass/seismic/TimeSeries.h"
#include "mspass/seismic/Seismogram.h"
#include "mspass/seismic/Ensemble.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/utility/VectorStatistics.h"
#include "mspass/utility/Metadata.h"
#include "mspass/algorithms/TimeWindow.h"
#include "mspass/algorithms/algorithms.h"
namespace mspass::algorithms::amplitudes{
double PeakAmplitude(const mspass::seismic::CoreTimeSeries& d);
double PeakAmplitude(const mspass::seismic::CoreSeismogram& d);
double RMSAmplitude(const mspass::seismic::CoreTimeSeries& d);
double RMSAmplitude(const mspass::seismic::CoreSeismogram& d);
double PercAmplitude(const mspass::seismic::CoreTimeSeries& d,const double perf);
double PercAmplitude(const mspass::seismic::CoreSeismogram& d,const double perf);
double MADAmplitude(const mspass::seismic::CoreTimeSeries& d);
double MADAmplitude(const mspass::seismic::CoreSeismogram& d);
enum class ScalingMethod
{
  Peak, 
  RMS,  
  ClipPerc, 
  MAD   
};
const std::string scale_factor_key("calib");
template <typename Tdata> double scale(Tdata& d,const ScalingMethod method,
  const double level, const mspass::algorithms::TimeWindow win)
{
  if((method==ScalingMethod::ClipPerc) && (level<=0.0 || level>1.0))
    throw mspass::utility::MsPASSError("scale function:  illegal perf level specified for clip percentage scale - must be between 0 and 1\nData unaltered - may cause downstream problems",
       mspass::utility::ErrorSeverity::Suspect);
  try{
    double newcalib(1.0);
    /* the else condition here should perhaps generate an elog message but
    did not implement to allow this template to be used for CoreTimeSeries
    and CoreSeismogram that do not have an elog attribute.*/
    if(d.is_defined(scale_factor_key))
    {
      newcalib=d.get_double(scale_factor_key);
    }
    /* Handle time windowing. Log window mismatches but silently handle
    cast where the window is invalid - used as a way to override any
    time windowing */
    mspass::algorithms::TimeWindow ampwindow;
    if(win.start>win.end)
    {
      ampwindow.start=d.t0();
      ampwindow.end=d.endtime();
    }
    else if( (fabs(win.start-d.t0())/d.dt()>0.5)
       || (fabs(win.end-d.endtime())/d.dt() > 0.5) )
    {
      std::stringstream ss;
      ss << "Window time range is inconsistent with input data range"<<std::endl
         << "Input data starttime="<<d.t0()<<" and window start time="
         << win.start <<" Difference="<<d.t0()-win.start<<std::endl
         << "Input data endtime="<<d.endtime()<<" and window end time="
         << win.end <<" Difference="<<d.endtime()-win.end<<std::endl
         << "One or the other exceeds 1/sample interval="<<d.dt()<<std::endl
         << "Window for amplitude calculation changed to data range";
      d.elog.log_error("scale",ss.str(),mspass::utility::ErrorSeverity::Complaint);
      ampwindow.start=d.t0();
      ampwindow.end=d.endtime();
    }
    else
    {
      ampwindow=win;
    }
    Tdata windowed_data;
    windowed_data=mspass::algorithms::WindowData(d,ampwindow);
    double amplitude,dscale;
    switch(method)
    {
      case ScalingMethod::Peak:
        amplitude=PeakAmplitude(windowed_data);
        dscale = level/amplitude;
        newcalib /= dscale;
        break;
      case ScalingMethod::ClipPerc:
        amplitude=PercAmplitude(windowed_data,level);
        /* for this scaling we use level as perf and output level is frozen
        to be scaled to order unity*/
        dscale = 1.0/amplitude;
        newcalib /= dscale;
        break;
      case ScalingMethod::MAD:
        amplitude=MADAmplitude(windowed_data);
        dscale = level/amplitude;
        newcalib /= dscale;
        break;
      case ScalingMethod::RMS:
      default:
        amplitude=RMSAmplitude(windowed_data);
        dscale = level/amplitude;
        newcalib /= dscale;
    };
    d *= dscale;
    d.put(scale_factor_key,newcalib);
    return amplitude;
  }catch(...){throw;};
}
template <typename Tdata> std::vector<double> scale_ensemble_members(mspass::seismic::Ensemble<Tdata>& d,
  const ScalingMethod& method, const double level, const mspass::algorithms::TimeWindow win)
{
  if((method==ScalingMethod::ClipPerc) && (level<=0.0 || level>1.0))
    throw mspass::utility::MsPASSError("scale_ensemble_members function:  illegal perf level specified for clip percentage scale - must be between 0 and 1\nData unaltered - may cause downstream problems",
       mspass::utility::ErrorSeverity::Suspect);
  try{
    typename std::vector<Tdata>::iterator dptr;
    std::vector<double> amps;
    amps.reserve(d.member.size());
    for(dptr=d.member.begin();dptr!=d.member.end();++dptr)
    {
      double thisamp;
      thisamp=scale(*dptr,method,level,win);
      amps.push_back(thisamp);
    }
    return amps;
  }catch(...){throw;};
}
template <typename Tdata> double scale_ensemble(mspass::seismic::Ensemble<Tdata>& d,
  const ScalingMethod& method, const double level, const bool use_mean)
{
  if((method==ScalingMethod::ClipPerc) && (level<=0.0 || level>1.0))
    throw mspass::utility::MsPASSError("scale_ensemble function:  illegal perf level specified for clip percentage scale - must be between 0 and 1\nData unaltered - may cause downstream problems",
       mspass::utility::ErrorSeverity::Suspect);
  try{
    double avgamp;   //defined here because the value computed here is returned on success
    typename std::vector<Tdata>::iterator dptr;
    std::vector<double> amps;
    amps.reserve(d.member.size());
    size_t nlive(0);
    for(dptr=d.member.begin();dptr!=d.member.end();++dptr)
    {
      double amplitude;
      if(dptr->dead()) continue;
      switch(method)
      {
        case ScalingMethod::Peak:
          amplitude=PeakAmplitude(*dptr);
          break;
        case ScalingMethod::ClipPerc:
          amplitude=PercAmplitude(*dptr,level);
          break;
        case ScalingMethod::MAD:
          amplitude=MADAmplitude(*dptr);
          break;
        case ScalingMethod::RMS:
        default:
          amplitude=RMSAmplitude(*dptr);
      };
      ++nlive;
      amps.push_back(log(amplitude));
    }
    /*Silently return a 0 if there are no live data members*/
    if(nlive==0) return 0.0;
    mspass::utility::VectorStatistics<double> ampstats(amps);
    if(use_mean)
    {
      avgamp=ampstats.mean();
    }
    else
    {
      avgamp=ampstats.median();
    }
 
    /* restore to a value instead of natural log*/
    avgamp=exp(avgamp);
    double dscale=level/avgamp;
    /* Now scale the data and apply calib */
    for(dptr=d.member.begin();dptr!=d.member.end();++dptr)
    {
      if(dptr->live())
      {
        double calib;
        (*dptr) *= dscale;
        if(dptr->is_defined(scale_factor_key))
        {
          calib=dptr->get_double(scale_factor_key);
        }
        else
        {
          calib=1.0;
        }
        calib/=dscale;
        dptr->put(scale_factor_key,calib);
      }
    }
    return avgamp;
  }catch(...){throw;};
}
class BandwidthData
{
public:
  double low_edge_f;
  double high_edge_f;
  double low_edge_snr;
  double high_edge_snr;
  /* This is the frequency range of the original data */
  double f_range;
  BandwidthData()
  {
    low_edge_f=0.0;
    high_edge_f=0.0;
    low_edge_snr=0.0;
    high_edge_snr=0.0;
    f_range=0.0;
  };
  double bandwidth_fraction() const
  {
    if(f_range<=0.0)
      return 0.0;
    else
      return (high_edge_f-low_edge_f)/f_range;
  };
  double bandwidth() const
  {
    if(f_range<=0.0)
    return 0.0;
    else
    {
      double ratio=high_edge_f/low_edge_f;
      return 20.0*log10(ratio);
    }
  };
};
BandwidthData EstimateBandwidth(const double signal_df,
  const mspass::seismic::PowerSpectrum& s, const mspass::seismic::PowerSpectrum& n,
    const double snr_threshold, const double tbp,const double fhs=-1.0,
     const bool fix_high_edge_to_fhs=false);
mspass::utility::Metadata BandwidthStatistics(const mspass::seismic::PowerSpectrum& s,
  const mspass::seismic::PowerSpectrum& n,
      const BandwidthData& bwd);
} // namespace end
#endif