MTPowerSpectrumEngine.h

#ifndef __MTPOWERSPECTRUM_ENGINE_H__
#define __MTPOWERSPECTRUM_ENGINE_H__
 
#include "mspass/algorithms/deconvolution/GSLFFTResources.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/seismic/TimeSeries.h"
#include "mspass/utility/dmatrix.h"
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <cmath>
#include <memory>
#include <string>
#include <vector>
 
namespace mspass::algorithms::deconvolution {
class MTPowerSpectrumEngine {
public:
  MTPowerSpectrumEngine();
  MTPowerSpectrumEngine(const int winsize, const double tbp, const int ntapers,
                        const int nfftin = -1, const double dtin = 1.0);
  MTPowerSpectrumEngine(const MTPowerSpectrumEngine &parent);
  ~MTPowerSpectrumEngine();
  MTPowerSpectrumEngine &operator=(const MTPowerSpectrumEngine &parent);
  mspass::seismic::PowerSpectrum apply(const mspass::seismic::TimeSeries &d);
  std::vector<double> apply(const std::vector<double> &d);
  double df() const { return deltaf; };
  std::vector<double> frequencies();
  int taper_length() const { return taperlen; };
  double time_bandwidth_product() const { return tbp; };
  int number_tapers() const { return ntapers; };
  int fftsize() const { return nfft; };
  double dt() const { return operator_dt; };
  double set_df(double dt) {
    const std::string caller("MTPowerSpectrumEngine::set_df");
    if (!std::isfinite(dt) || dt <= 0.0)
      throw mspass::utility::MsPASSError(
          caller + ": sample interval must be finite and positive",
          mspass::utility::ErrorSeverity::Invalid);
    const int this_nf = this->nf();
    if (this_nf <= 1)
      throw mspass::utility::MsPASSError(
          caller + ": engine fft length is not configured",
          mspass::utility::ErrorSeverity::Invalid);
    const double new_deltaf = 1.0 / (static_cast<double>(this->nfft) * dt);
    if (!std::isfinite(new_deltaf) || new_deltaf <= 0.0)
      throw mspass::utility::MsPASSError(
          caller + ": sample interval produces an invalid frequency spacing",
          mspass::utility::ErrorSeverity::Invalid);
    this->operator_dt = dt;
    this->deltaf = new_deltaf;
    return deltaf;
  };
  int nf() {
    /* this simple formula depends upon integer truncation when used with
    nfft as an odd number.   For reference, this is what prieto uses in
    the python multitaper package:
    if (nfft%2 == 0):
        nf = int(nfft/2 + 1)
    else:
        nf = int((nfft+1)/2)
    they will yield the same result but this is simpler and faster
    */
    return (this->nfft) / 2 + 1;
  };
 
private:
  int taperlen;
  int ntapers;
  int nfft;
  double tbp;
  double operator_dt;
  mspass::utility::dmatrix tapers;
  /* Frequency bin interval of last data processed.*/
  double deltaf;
  gsl_fft_complex_wavetable *wavetable;
  gsl_fft_complex_workspace *workspace;
  friend boost::serialization::access;
  template <class Archive>
  void save(Archive &ar, const unsigned int version) const {
    ar & taperlen;
    ar & ntapers;
    ar & nfft;
    ar & tbp;
    ar & operator_dt;
    ar & tapers;
    ar & deltaf;
  }
  template <class Archive> void load(Archive &ar, const unsigned int version) {
    int loaded_taperlen;
    int loaded_ntapers;
    int loaded_nfft;
    double loaded_tbp;
    double loaded_operator_dt;
    mspass::utility::dmatrix loaded_tapers;
    double loaded_deltaf;
    ar & loaded_taperlen;
    ar & loaded_ntapers;
    ar & loaded_nfft;
    ar & loaded_tbp;
    ar & loaded_operator_dt;
    ar & loaded_tapers;
    ar & loaded_deltaf;
 
    const std::string caller("MTPowerSpectrumEngine serialization load");
    if (loaded_nfft <= 0 || loaded_taperlen <= 0 || loaded_ntapers <= 0)
      throw mspass::utility::MsPASSError(
          caller + ": archived fft, taper, and taper-count lengths must be "
                   "positive",
          mspass::utility::ErrorSeverity::Invalid);
    if (loaded_nfft < loaded_taperlen)
      throw mspass::utility::MsPASSError(
          caller + ": archived fft length is shorter than taper length",
          mspass::utility::ErrorSeverity::Invalid);
    if (!std::isfinite(loaded_tbp) || loaded_tbp <= 0.0 ||
        !std::isfinite(loaded_operator_dt) || loaded_operator_dt <= 0.0 ||
        !std::isfinite(loaded_deltaf) || loaded_deltaf <= 0.0)
      throw mspass::utility::MsPASSError(
          caller + ": archived spectral parameters must be finite and "
                   "positive",
          mspass::utility::ErrorSeverity::Invalid);
    if (!loaded_tapers.storage_is_consistent() ||
        loaded_tapers.rows() != static_cast<size_t>(loaded_ntapers) ||
        loaded_tapers.columns() != static_cast<size_t>(loaded_taperlen))
      throw mspass::utility::MsPASSError(
          caller + ": archived taper matrix dimensions are inconsistent",
          mspass::utility::ErrorSeverity::Invalid);
    auto resources = detail::AllocateGSLFFTResources(loaded_nfft, caller);
 
    loaded_tapers.swap(tapers);
    if (wavetable != nullptr)
      gsl_fft_complex_wavetable_free(wavetable);
    if (workspace != nullptr)
      gsl_fft_complex_workspace_free(workspace);
    taperlen = loaded_taperlen;
    ntapers = loaded_ntapers;
    nfft = loaded_nfft;
    tbp = loaded_tbp;
    operator_dt = loaded_operator_dt;
    deltaf = loaded_deltaf;
    wavetable = resources.wavetable.release();
    workspace = resources.workspace.release();
  }
  BOOST_SERIALIZATION_SPLIT_MEMBER()
};
} // namespace mspass::algorithms::deconvolution
#endif