FFTDeconOperator.h

#ifndef __FFT_DECON_OPERATOR_H__
#define __FFT_DECON_OPERATOR_H__
#include "mspass/algorithms/TimeWindow.h"
#include "mspass/algorithms/deconvolution/ComplexArray.h"
#include "mspass/algorithms/deconvolution/GSLFFTResources.h"
#include "mspass/seismic/CoreTimeSeries.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/utility/Metadata.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 <string>
namespace mspass::algorithms::deconvolution {
class FFTDeconOperator {
public:
  FFTDeconOperator();
  FFTDeconOperator(const mspass::utility::Metadata &md);
  FFTDeconOperator(const FFTDeconOperator &parent);
  ~FFTDeconOperator();
  FFTDeconOperator &operator=(const FFTDeconOperator &parent);
  void changeparameter(const mspass::utility::Metadata &md);
  void change_size(const int nfft_new);
  void change_shift(const int shift) { sample_shift = shift; };
  int get_size() { return nfft; };
  int get_shift() { return sample_shift; };
  int operator_size() { return static_cast<int>(nfft); };
  int operator_shift() { return sample_shift; };
  double df(const double dt) {
    double period;
    period = static_cast<double>(nfft) * dt;
    return 1.0 / period;
  };
  mspass::seismic::CoreTimeSeries FourierInverse(const ComplexArray &winv,
                                                 const ComplexArray &sw,
                                                 const double dt,
                                                 const double t0parent);
 
protected:
  int nfft;
  int sample_shift;
  gsl_fft_complex_wavetable *wavetable;
  gsl_fft_complex_workspace *workspace;
  ComplexArray winv;
 
private:
  friend boost::serialization::access;
  template <class Archive>
  void save(Archive &ar, const unsigned int version) const {
    // std::cout << "Entered FFTDecon serialization save function"<<std::endl;
    ar & nfft;
    ar & sample_shift;
    ar & winv;
    // std::cout << "Exiting save function" << std::endl;
  }
  template <class Archive> void load(Archive &ar, const unsigned int version) {
    int loaded_nfft;
    int loaded_sample_shift;
    ComplexArray loaded_winv;
    ar & loaded_nfft;
    ar & loaded_sample_shift;
    ar & loaded_winv;
    const std::string caller("FFTDeconOperator serialization load");
    if (loaded_nfft <= 0)
      throw mspass::utility::MsPASSError(
          caller + ": archived fft length must be positive",
          mspass::utility::ErrorSeverity::Invalid);
    if (loaded_sample_shift < 0 || loaded_sample_shift > loaded_nfft)
      throw mspass::utility::MsPASSError(
          caller + ": archived sample shift is outside the fft buffer",
          mspass::utility::ErrorSeverity::Invalid);
    /* An operator may be serialized before its inverse is initialized. */
    if (loaded_winv.size() != 0 && loaded_winv.size() != loaded_nfft)
      throw mspass::utility::MsPASSError(
          caller + ": archived inverse size does not match fft length",
          mspass::utility::ErrorSeverity::Invalid);
    auto resources = detail::AllocateGSLFFTResources(loaded_nfft, caller);
 
    /* Everything above may throw.  Swaps, frees, scalar stores, and releases
     * below form the no-throw commit into either a default or initialized
     * target. */
    winv.swap(loaded_winv);
    if (wavetable != nullptr)
      gsl_fft_complex_wavetable_free(wavetable);
    if (workspace != nullptr)
      gsl_fft_complex_workspace_free(workspace);
    nfft = loaded_nfft;
    sample_shift = loaded_sample_shift;
    wavetable = resources.wavetable.release();
    workspace = resources.workspace.release();
  }
  BOOST_SERIALIZATION_SPLIT_MEMBER()
};
 
/* This helper is best referenced here */
 
std::vector<double> circular_shift(const std::vector<double> &d, const int i0);
int ComputeFFTLength(const mspass::algorithms::TimeWindow w, const double dt);
void ValidateWindowDuration(const mspass::algorithms::TimeWindow w,
                            const std::string &window_name,
                            const std::string &caller);
int ComputeFFTLength(const mspass::utility::Metadata &md);
int ComputeDeconSampleShift(const mspass::utility::Metadata &md);
void ValidatePowerSpectrumCoversDC(
    const mspass::seismic::PowerSpectrum &spectrum,
    const std::string &caller);
std::vector<double> ExtractLagWindow(ComplexArray &fft_buffer,
                                     const int output_length,
                                     const int sample_shift);
extern "C" {
unsigned int nextPowerOf2(unsigned int n);
}
} // namespace mspass::algorithms::deconvolution
#endif