amplitudes.h

#ifndef _AMPLITUDES_H_
#define _AMPLITUDES_H_
#include "mspass/algorithms/TimeWindow.h"
#include "mspass/algorithms/algorithms.h"
#include "mspass/seismic/Ensemble.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/seismic/Seismogram.h"
#include "mspass/seismic/TimeSeries.h"
#include "mspass/utility/Metadata.h"
#include "mspass/utility/MsPASSError.h"
#include "mspass/utility/VectorStatistics.h"
#include <algorithm>
#include <cmath>
#include <sstream>
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);
    }
    /* A reversed window is the public signal for full-record scaling.
    Otherwise measure only the physical intersection with the datum. */
    mspass::algorithms::TimeWindow ampwindow;
    if (win.start > win.end) {
      ampwindow.start = d.t0();
      ampwindow.end = d.endtime();
    } else {
      ampwindow.start = std::max(win.start, d.t0());
      ampwindow.end = std::min(win.end, d.endtime());
      if (ampwindow.start > ampwindow.end) {
        std::stringstream ss;
        ss << "scale:  amplitude measurement window [" << win.start << ", "
           << win.end << "] has no intersection with data "
           << "range [" << d.t0() << ", " << d.endtime() << "]";
        throw mspass::utility::MsPASSError(
            ss.str(), mspass::utility::ErrorSeverity::Invalid);
      }
    }
    Tdata windowed_data;
    windowed_data = mspass::algorithms::WindowData(d, ampwindow);
    double amplitude, dscale;
    switch (method) {
    case ScalingMethod::Peak:
      amplitude = PeakAmplitude(windowed_data);
      break;
    case ScalingMethod::ClipPerc:
      amplitude = PercAmplitude(windowed_data, level);
      break;
    case ScalingMethod::MAD:
      amplitude = MADAmplitude(windowed_data);
      break;
    case ScalingMethod::RMS:
    default:
      amplitude = RMSAmplitude(windowed_data);
    };
    /* needed to handle case with a vector of all 0s*/
    if (amplitude > 0.0) {
      dscale = level / amplitude;
      newcalib /= dscale;
      d *= dscale;
      d.put(scale_factor_key, newcalib);
    } else {
      std::stringstream ss;
      ss << "Data array is all 0s and cannot be scaled";
      d.elog.log_error("scale", ss.str(),
                       mspass::utility::ErrorSeverity::Complaint);
      /* This may not be necessary but it assures this value is always set on
      return even if it means nothing*/
      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());
    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);
      };
      if ((amplitude > 0.0) && std::isfinite(amplitude))
        amps.push_back(std::log(amplitude));
    }
    /* Silently return 0 without mutation when no member can define a gain. */
    if (amps.empty())
      return 0.0;
    if (amps.size() == 1) {
      avgamp = amps.front();
    } else {
      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 = std::exp(avgamp);
    /* A core ensemble has no error-log channel, so a nonpositive result is a
    silent no-op consistent with the no-eligible-member case above. */
    if (avgamp <= 0.0) {
      return 0.0;
    }
    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;
  };
}
template <class T> std::vector<T> normalize(const std::vector<T> &d) {
  size_t N = d.size();
  std::vector<T> result;
  result.reserve(N);
  double d_nrm(0.0);
  ;
  for (size_t i = 0; i < N; ++i) {
    d_nrm += (d[i] * d[i]);
    result.push_back(d[i]);
  }
  d_nrm = sqrt(d_nrm);
  for (size_t i = 0; i < N; ++i)
    result[i] /= d_nrm;
  return result;
}
class BandwidthData {
public:
  double low_edge_f;
  double high_edge_f;
  double low_edge_snr;
  double high_edge_snr;
  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 {
    /* All these conditionals are necessary for handling unexpected
     * values.   0 is effectively and error return.  Without these
     * the function can return NaN or generate floating point exceptions. */
    if ((f_range <= 0.0) || (high_edge_f <= low_edge_f) ||
        (high_edge_f < 0.0) || (low_edge_f <= 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 mspass::algorithms::amplitudes
#endif