MultiTaperXcorDecon.h

#ifndef __SIMPLE_MULTITAPER_DECON_H__
#define __SIMPLE_MULTITAPER_DECON_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 MultiTaperXcorDecon : public FFTDeconOperator, public ScalarDecon {
public:
  MultiTaperXcorDecon() : FFTDeconOperator(), ScalarDecon() {};
  MultiTaperXcorDecon(const mspass::utility::Metadata &md,
                      const std::vector<double> &noise,
                      const std::vector<double> &wavelet,
                      const std::vector<double> &data);
  MultiTaperXcorDecon(const mspass::utility::Metadata &md);
  MultiTaperXcorDecon(const MultiTaperXcorDecon &parent);
  ~MultiTaperXcorDecon() {};
  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();
  mspass::utility::Metadata QCMetrics();
  int get_taperlen() { return taperlen; };
  int get_number_tapers() { return nseq; };
  double get_time_bandwidth_product() { return nw; };
 
private:
  /* noise data vector */
  std::vector<double> noise;
  double nw, damp;
  int nseq; // number of tapers
  unsigned int taperlen;
  mspass::utility::dmatrix tapers;
  /* 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 */
  ComplexArray ao_fft;
  int read_metadata(const mspass::utility::Metadata &md, bool refresh);
  /* Returns a tapered data in container of ComplexArray objects*/
  std::vector<ComplexArray> taper_data(const std::vector<double> &signal);
  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 & ao_fft;
  }
};
} // namespace mspass::algorithms::deconvolution
#endif