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/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) {
    // std::cout << "Entered FFTDecon serializaton load function" << std::endl;
    ar &this->nfft;
    ar &this->sample_shift;
    // std::cout << "Loading winv vector" << std::endl;
    ar & winv;
    // std::cout << "Creating wavetable " << std::endl;
    this->wavetable = gsl_fft_complex_wavetable_alloc(this->nfft);
    // std::cout << "Creating workspace" << std::endl;
    this->workspace = gsl_fft_complex_workspace_alloc(this->nfft);
    // std::cout << "Exiting load" << std::endl;
  }
  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