CNR3CDecon.h

#ifndef __CNR3C_DECON_H__
#define __CNR3C_DECON_H__
#include <vector>
#include <math.h>
#include "mspass/utility/AntelopePf.h"
#include "mspass/algorithms/deconvolution/FFTDeconOperator.h"
#include "mspass/algorithms/deconvolution/ShapingWavelet.h"
#include "mspass/algorithms/deconvolution/MTPowerSpectrumEngine.h"
#include "mspass/algorithms/amplitudes.h"
#include "mspass/algorithms/Taper.h"
#include "mspass/algorithms/TimeWindow.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/seismic/TimeSeries.h"
#include "mspass/seismic/Seismogram.h"
namespace mspass::algorithms::deconvolution{
 
class Base3CDecon
{
public:
    virtual ~Base3CDecon() {};
    /*
    virtual void change_parameters(const mspass::BasicMetadata &md)=0;
    virtual void loaddata(mspass::Seismogram& d,const int comp)=0;
    virtual void loadwavelet(const mspass::TimeSeries& w)=0;
    */
    /* \brief Return the ideal output of the deconvolution operator.
 
    All deconvolution operators have a implicit or explicit ideal output
    signal. e.g. for a spiking Wiener filter it is a delta function with or
    without a lag.  For a shaping wavelt it is the time domain version of the
    wavelet. */
    virtual mspass::seismic::Seismogram process()=0;
    virtual mspass::seismic::TimeSeries actual_output()=0;
 
    //virtual mspass::TimeSeries inverse_wavelet() = 0;
    virtual mspass::seismic::TimeSeries inverse_wavelet(double) = 0;
    virtual mspass::utility::Metadata QCMetrics()=0;
};
/* This enum 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.  The main difference is it avoids the messy trampoline class
needed with pybind11 to create wrappers for the python bindings.
This enum will not be exposed to python as it is totally internal to the
CNR3CDecon class.
*/
enum class CNR3C_algorithms{
    generalized_water_level,
    colored_noise_damping,
    undefined
};
//class CNR3CDecon : public mspass::FFTDeconOperator
class CNR3CDecon : public Base3CDecon, public FFTDeconOperator
{
public:
  CNR3CDecon();
  CNR3CDecon(const mspass::utility::AntelopePf& pf);
  CNR3CDecon(const CNR3CDecon& parent);
  ~CNR3CDecon();
  CNR3CDecon& operator=(const CNR3CDecon& parent);
  void change_parameters(const mspass::utility::BasicMetadata& md);
  void loaddata(mspass::seismic::Seismogram& d, const int wcomp,const bool loadnoise=false);
  void loaddata(mspass::seismic::Seismogram& d, const bool loadnoise=false);
  void loadnoise_data(const mspass::seismic::Seismogram& n);
  void loadnoise_data(const mspass::seismic::PowerSpectrum& n);
  void loadwavelet(const mspass::seismic::TimeSeries& w);
  void loadnoise_wavelet(const mspass::seismic::TimeSeries& n);
  void loadnoise_wavelet(const mspass::seismic::PowerSpectrum& n);
  /* These same names are used in ScalarDecon but we don't inherit them
  here because this algorithm is 3C data centric there is a collision
  with the ScalarDecon api because of it.  */
  mspass::seismic::Seismogram process();
 
  /* \brief Return the ideal output of the deconvolution operator.
 
  All deconvolution operators have a implicit or explicit ideal output
  signal. e.g. for a spiking Wiener filter it is a delta function with or
  without a lag.  For a shaping wavelt it is the time domain version of the
  wavelet. */
  mspass::seismic::TimeSeries ideal_output();
  mspass::seismic::TimeSeries actual_output();
 
  mspass::seismic::TimeSeries inverse_wavelet(double tshift);
  mspass::utility::Metadata QCMetrics();
  mspass::seismic::PowerSpectrum wavelet_noise_spectrum()
  {
    return psnoise;
  };
  mspass::seismic::PowerSpectrum data_noise_spectrum()
  {
    return psnoise_data;
  };
  mspass::seismic::PowerSpectrum wavelet_spectrum()
  {
    return pswavelet;
  };
  mspass::seismic::PowerSpectrum data_spectrum()
  {
    return pssignal;
  };
private:
  CNR3C_algorithms algorithm;
  bool taper_data;  //Set false only if none specified
  double operator_dt;   // Data must match this sample interval
  /* Expected time window size in samples
  (computed from processing_window and operator dt)*/
  int winlength;
  double decon_bandwidth_cutoff;
  /* Added Dec 2021 - this defines the upper starting frequency used for
  the EstimateBandwith function.   */
  double fhs;
  /* Defines relative time time window - ignored if length of input is
  consistent with number of samples expected in this window */
  mspass::algorithms::TimeWindow processing_window;
  mspass::algorithms::TimeWindow noise_window;
  MTPowerSpectrumEngine signalengine,waveletengine;
  MTPowerSpectrumEngine dnoise_engine, wnoise_engine;
  ShapingWavelet shapingwavelet;
  /* This contains the noise power spectrum to use for regularization
  of the inverse.  It should normally be created from a longer window
  than the data.
  */
  mspass::seismic::PowerSpectrum psnoise;
  /* This contains the power spectrum of the data used to estimate
     snr-based QC estimates.   It can be the same as the data but
     not necessarily.
  */
  mspass::seismic::PowerSpectrum psnoise_data;
  /* Because of the design of the algorithm we also have to save a power
  spectral estimate for the wavelet and data signals.  We use those when
  automatic bandwidth adjustment is enabled.*/
  mspass::seismic::PowerSpectrum pssignal;
  mspass::seismic::PowerSpectrum pswavelet;
  /* Cached data to be deconvolved - result of loaddata methds*/
  mspass::seismic::Seismogram decondata;
  /* Cached wavelet for deconvolution - result of loadwavelet*/
  mspass::seismic::TimeSeries wavelet;
  /* As the name suggest we allow different tapers for data and wavelet */
  std::shared_ptr<mspass::algorithms::BasicTaper> wavelet_taper;
  std::shared_ptr<mspass::algorithms::BasicTaper> data_taper;
  /* 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, noise_floor;
  /* 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;
  /* this parameter does a similar thing to regularization floor but is
  by the bandwidth estimation algorithm to define the edge of the working
  frequency band. */
  double snr_bandwidth;
  //ComplexArray winv;
  /* winv is in FFTDeconOperator*/
  ComplexArray ao_fft;
  mspass::algorithms::amplitudes::BandwidthData wavelet_bwd;
  mspass::algorithms::amplitudes::BandwidthData signal_bwd;
  /* We cache wavelet snr time series as it is more efficiently computed during
     the process routine and then used in (optional) qc methods */
  std::vector<double> wavelet_snr;
  /* SNR bandbwidth estimates count frequencies with snr above this value */
  double band_snr_floor;
  /* This array stores snr band fractions for each component.*/
  double signal_bandwidth_fraction[3];
  double peak_snr[3];
  double regularization_bandwidth_fraction;
  void read_parameters(const mspass::utility::AntelopePf& pf);
  int TestSeismogramInput(mspass::seismic::Seismogram& d,const int comp,const bool loaddata);
  void compute_gwl_inverse();
  void compute_gdamp_inverse();
  mspass::seismic::PowerSpectrum ThreeCPower(const mspass::seismic::Seismogram& d);
  void update_shaping_wavelet(const mspass::algorithms::amplitudes::BandwidthData& bwd);
};
}  // End namespace
 
#endif