CNRDeconEngine.h

#ifndef __CNR_DECON_ENGINE_H__
#define __CNR_DECON_ENGINE_H__
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include <boost/serialization/version.hpp>
 
#include "mspass/algorithms/Taper.h"
#include "mspass/algorithms/deconvolution/FFTDeconOperator.h"
#include "mspass/algorithms/deconvolution/MTPowerSpectrumEngine.h"
#include "mspass/algorithms/deconvolution/ShapingWavelet.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/seismic/Seismogram.h"
#include "mspass/seismic/TimeSeries.h"
#include "mspass/utility/AntelopePf.h"
 
namespace mspass::algorithms::deconvolution {
/* This enum is file scope to intentionally exclude it from python wrappers.
It is used internally to define the algorithm the processor is to run.
I (glp) chose that approach over the inheritance approach used in the scalar
methods as an artistic choice.   It is a matter of opinion which approach
is better.  This makes one symbol do multiple things with changes done in
the parameter setup as opposed to having to select the right symbolic name
to construct.  Anyway, this enum defines algorithms that can be chosen for
processing.
*/
enum class CNR3C_algorithms { generalized_water_level, colored_noise_damping };
class CNRDeconEngine : public FFTDeconOperator {
public:
  CNRDeconEngine();
  /* design note - delete when finished.
 
     The constructor uses the pf to initialize operator properties.
     Important in that is the multitaper engine that has a significant
     initialization cost.   the initialize_inverse_operator is a messy
     way to define load data and then compute the inverse stored
     internally.   Doing it that way, however, allows the same code
     to be used for both single station and array decon.   In single
     station use every datum has to call initiallize_inverse_operator
     while with an array decon the operator is define once for an
     ensemble and applied to all members.   THE ASSUMPTION is all the
     grungy work to assure all that gets done correction is handled
     in python.
     */
  CNRDeconEngine(const mspass::utility::AntelopePf &pf);
  CNRDeconEngine(const CNRDeconEngine &parent);
  void changeparameter(const mspass::utility::Metadata &md);
  void
  initialize_inverse_operator(const mspass::seismic::TimeSeries &wavelet,
                              const mspass::seismic::TimeSeries &noise_data);
  void initialize_inverse_operator(
      const mspass::seismic::TimeSeries &wavelet,
      const mspass::seismic::PowerSpectrum &noise_spectrum);
  virtual ~CNRDeconEngine() {};
  mspass::seismic::Seismogram
  process(const mspass::seismic::Seismogram &d,
          const mspass::seismic::PowerSpectrum &psnoise, const double fl,
          const double fh);
  mspass::seismic::Seismogram
  process(const mspass::seismic::Seismogram &d,
          const mspass::seismic::PowerSpectrum &psnoise);
  double get_operator_dt() const { return this->operator_dt; }
  mspass::seismic::PowerSpectrum
  compute_noise_spectrum(const mspass::seismic::TimeSeries &d2use);
  mspass::seismic::PowerSpectrum
  compute_noise_spectrum(const mspass::seismic::Seismogram &d2use);
  mspass::seismic::TimeSeries output_shaping_wavelet() {
    return this->ideal_output();
  }
  mspass::seismic::TimeSeries ideal_output();
  mspass::seismic::TimeSeries
  actual_output(const mspass::seismic::TimeSeries &wavelet);
  mspass::seismic::TimeSeries
  resolution_kernel(const mspass::seismic::TimeSeries &wavelet) {
    return this->actual_output(wavelet);
  }
  mspass::seismic::TimeSeries
  inverse_wavelet(const mspass::seismic::TimeSeries &wavelet,
                  const double t0shift);
  mspass::utility::Metadata QCMetrics();
  CNRDeconEngine &operator=(const CNRDeconEngine &parent);
 
private:
  CNR3C_algorithms algorithm;
  /* For the colored noise damping algorithm the damper is frequency dependent.
     The same issue in water level that requires a floor on the water level
     applies to damping.   We use noise_floor to create a lower bound on
     damper values.   Note the damping constant at each frequency is
     damp*noise except where noise is below noise_floor defined relative to
     maximum noise value where it is set to n_peak*noise_floor*damp. */
  double damp;
  double noise_floor;
  /* SNR bandbwidth estimates count frequencies with snr above this value */
  double band_snr_floor;
  double operator_dt; // Data must match this sample interval
  int shaping_wavelet_number_poles;
  mspass::algorithms::deconvolution::ShapingWavelet shapingwavelet;
  /* PF-defined shaping state.  Legacy process(d,psnoise,fl,fh) may replace
   * shapingwavelet for one call, but configured processing restores this. */
  mspass::algorithms::deconvolution::ShapingWavelet configured_shapingwavelet;
  /* Expected time window size in samples.   When signal lengths
  match this value the slepian tapers are not recomputed.  When there
  is a mismatch it will change.  That means this can change dynamically
  when run on multiple data objects. */
  int winlength;
  mspass::algorithms::deconvolution::MTPowerSpectrumEngine signal_engine;
  mspass::algorithms::deconvolution::MTPowerSpectrumEngine noise_engine;
 
  /* This algorithm uses a mix of damping and water level.   Above this floor,
  which acts a bit like a water level, no regularization is done.  If
  snr is less than this value we regularize with damp*noise_amplitude.
  Note the noise_floor parameter puts a lower bound on the frequency dependent
  regularization.   If noise amplitude (not power) is less than noise_floor
  the floor is set like a water level as noise_max*noise_level.*/
  double snr_regularization_floor;
  /* These are QC metrics computed by process method.  Saved to allow them
  to be use in QCmetrics method. */
  double regularization_bandwidth_fraction;
  double peak_snr[3];
  double signal_bandwidth_fraction[3];
  mspass::algorithms::deconvolution::ComplexArray winv;
  /* This is the lag from sample 0 for the time defines as 0 for the
  wavelet used to compute the inverse.  It is needed to resolve time
  in processing and the actual_output method.*/
  int winv_t0_lag;
  /*** Private methods *****/
  mspass::seismic::Seismogram
  process_with_current_shaping(const mspass::seismic::Seismogram &d,
                               const mspass::seismic::PowerSpectrum &psnoise);
  void update_shaping_wavelet(const double fl, const double fh);
  /* These are two algorithms for computing inverse operator in the frequency
   * domain*/
  void compute_winv(const mspass::seismic::TimeSeries &wavelet,
                    const mspass::seismic::PowerSpectrum &psnoise);
  void compute_gwl_inverse(const mspass::seismic::TimeSeries &wavelet,
                           const mspass::seismic::PowerSpectrum &psnoise);
  void compute_gdamp_inverse(const mspass::seismic::TimeSeries &wavelet,
                             const mspass::seismic::PowerSpectrum &psnoise);
  friend boost::serialization::access;
  template <class Archive>
  void serialize(Archive &ar, const unsigned int version) {
    // std::cout <<"Entered serialize function"<<std::endl;
    ar &boost::serialization::base_object<FFTDeconOperator>(*this);
    // std::cout << "Serializing first group of simple parameters"<<std::endl;
    ar & algorithm;
    ar & damp;
    ar & noise_floor;
    ar & band_snr_floor;
    ar & operator_dt;
    ar & shaping_wavelet_number_poles;
    // std::cout << "Serializing shapingwavelet"<<std::endl;
    ar & shapingwavelet;
    if (version > 0) {
      ar & configured_shapingwavelet;
    } else if (Archive::is_loading::value) {
      configured_shapingwavelet = shapingwavelet;
    }
    ar & winlength;
    // std::cout<<"Serializing power spectrum engine objects"<<std::endl;
    ar & signal_engine;
    ar & noise_engine;
    ar & snr_regularization_floor;
    // std::cout << "Serializin winv vector"<<std::endl;
    ar & winv;
    ar & winv_t0_lag;
    // std::cout<<"Serializing final block of parameters"<<std::endl;
    ar & regularization_bandwidth_fraction;
    /* These fixed length arrays caused probems - seg faults.
     * Apparently boost doesn't handle that corectly.  There may
     * be a more concise way to do this but this should always work. */
    // std::cout << "Entering block for 3 component arrays"<<std::endl;
    for (auto k = 0; k < 3; ++k) {
      // std::cout << "k="<<k<<std::endl;
      ar &peak_snr[k];
      // std::cout << "bandwidth_fraction"<<std::endl;
      ar &signal_bandwidth_fraction[k];
    }
    // std::cout << "Exiting serialize function"<<std::endl;
  }
};
} // namespace mspass::algorithms::deconvolution
BOOST_CLASS_VERSION(mspass::algorithms::deconvolution::CNRDeconEngine, 1)
#endif