MultiTaperSpecDivDecon.h

#ifndef __PAVLIS_MULTITAPER_DECON_H__
#define __PAVLIS_MULTITAPER_DECON_H__
#include "mspass/algorithms/deconvolution/ComplexArray.h"
#include "mspass/algorithms/deconvolution/FFTDeconOperator.h"
#include "mspass/algorithms/deconvolution/ScalarDecon.h"
#include "mspass/algorithms/deconvolution/ShapingWavelet.h"
#include "mspass/seismic/CoreTimeSeries.h"
#include "mspass/utility/Metadata.h"
#include "mspass/utility/dmatrix.h"
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/vector.hpp>
#include <vector>
namespace mspass::algorithms::deconvolution {
class MultiTaperSpecDivDecon : public FFTDeconOperator, public ScalarDecon {
public:
  MultiTaperSpecDivDecon() : FFTDeconOperator(), ScalarDecon() {};
  MultiTaperSpecDivDecon(const mspass::utility::Metadata &md,
                         const std::vector<double> &noise,
                         const std::vector<double> &wavelet,
                         const std::vector<double> &data);
  MultiTaperSpecDivDecon(const mspass::utility::Metadata &md);
  MultiTaperSpecDivDecon(const MultiTaperSpecDivDecon &parent);
  ~MultiTaperSpecDivDecon() {};
  void changeparameter(const mspass::utility::Metadata &md) {
    this->read_metadata(md, true);
  };
  int loadnoise(const std::vector<double> &noise);
  int load(const std::vector<double> &w, const std::vector<double> &d,
           const std::vector<double> &n);
  void process();
  mspass::seismic::CoreTimeSeries actual_output();
  mspass::seismic::CoreTimeSeries inverse_wavelet(const double t0parent = 0.0);
  mspass::seismic::CoreTimeSeries inverse_wavelet();
  std::vector<mspass::seismic::CoreTimeSeries>
  all_inverse_wavelets(const double t0parent = 0.0);
  std::vector<mspass::seismic::CoreTimeSeries>
  all_rfestimates(const double t0parent = 0.0);
  std::vector<mspass::seismic::CoreTimeSeries>
  all_actual_outputs(const double t0parent = 0.0);
  mspass::utility::Metadata QCMetrics();
  int get_taperlen() { return taperlen; };
  int get_number_tapers() { return nseq; };
  double get_time_bandwidth_product() { return nw; };
 
private:
  std::vector<double> noise;
  double nw, damp;
  int nseq; // number of tapers
  unsigned int taperlen;
  mspass::utility::dmatrix tapers;
  /* With this algorithm we have to keep frequeny domain
   * representations of the inverse for each taper.  The
   * xcor version of this alorithm doesn't need that because
   * the averaging is different. */
  std::vector<ComplexArray> winv_taper;
  /* We also cache the actual output fft because the cost is
   * small compared to need to recompute it when requested.
   * This is a feature added for the GID method that adds
   * an inefficiency for straight application */
  std::vector<ComplexArray> ao_fft;
  /* This contains the rf estimates from each of the tapers.   They
  are averaged to produce final result, but we keep them for the
  option of bootstrap errors. */
  std::vector<ComplexArray> rfestimates;
  /* Private methods */
  /* Returns a tapered data in container of ComplexArray objects*/
  std::vector<ComplexArray> taper_data(const std::vector<double> &signal);
  int read_metadata(const mspass::utility::Metadata &md, bool refresh);
  int apply();
  friend boost::serialization::access;
  template <class Archive>
  void serialize(Archive &ar, const unsigned int version) {
    ar &boost::serialization::base_object<ScalarDecon>(*this);
    ar &boost::serialization::base_object<FFTDeconOperator>(*this);
    ar & noise;
    ar & nw;
    ar & damp;
    ar & nseq;
    ar & taperlen;
    ar & tapers;
    ar & winv_taper;
    ar & ao_fft;
    ar & rfestimates;
  }
};
} // namespace mspass::algorithms::deconvolution
#endif