MTPowerSpectrumEngine.h

#ifndef __MTPOWERSPECTRUM_ENGINE_H__
#define __MTPOWERSPECTRUM_ENGINE_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 <memory>
#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() { return operator_dt; };
  double set_df(double dt) {
    this->operator_dt = dt;
    int this_nf = this->nf();
    double fny = 1.0 / (2.0 * dt);
    this->deltaf = fny / static_cast<double>(this_nf - 1);
    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) {
    ar & taperlen;
    ar & ntapers;
    ar & nfft;
    ar & tbp;
    ar & operator_dt;
    ar & tapers;
    ar & deltaf;
    this->wavetable = gsl_fft_complex_wavetable_alloc(this->nfft);
    this->workspace = gsl_fft_complex_workspace_alloc(this->nfft);
  }
  BOOST_SERIALIZATION_SPLIT_MEMBER()
};
} // namespace mspass::algorithms::deconvolution
#endif