mspass::algorithms::deconvolution::LeastSquareDecon Class Reference

Inheritance diagram for mspass::algorithms::deconvolution::LeastSquareDecon:

Public Member Functions
	LeastSquareDecon (const LeastSquareDecon &parent)
	LeastSquareDecon (const mspass::utility::Metadata &md)
	LeastSquareDecon (const mspass::utility::Metadata &md, const std::vector< double > &wavelet, const std::vector< double > &data)
void	changeparameter (const mspass::utility::Metadata &md)
void	process ()
mspass::seismic::CoreTimeSeries	actual_output ()
	Return the actual output of the deconvolution operator.
mspass::seismic::CoreTimeSeries	inverse_wavelet (const double t0parent=0.0)
	Return a FIR respresentation of the inverse filter.
mspass::seismic::CoreTimeSeries	inverse_wavelet ()
	Return default FIR represesentation of the inverse filter.
mspass::utility::Metadata	QCMetrics ()
	Return appropriate quality measures.
Public Member Functions inherited from mspass::algorithms::deconvolution::FFTDeconOperator
	FFTDeconOperator (const mspass::utility::Metadata &md)
	FFTDeconOperator (const FFTDeconOperator &parent)
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)
int	get_size ()
int	get_shift ()
int	operator_size ()
int	operator_shift ()
double	df (const double dt)
mspass::seismic::CoreTimeSeries	FourierInverse (const ComplexArray &winv, const ComplexArray &sw, const double dt, const double t0parent)
	Return inverse wavelet for Fourier methods.
Public Member Functions inherited from mspass::algorithms::deconvolution::ScalarDecon
	ScalarDecon (const mspass::utility::Metadata &md)
	ScalarDecon (const std::vector< double > &d, const std::vector< double > &w)
	ScalarDecon (const ScalarDecon &parent)
int	load (const std::vector< double > &wavelet, const std::vector< double > &data)
	Load all data required for decon.
int	loaddata (const std::vector< double > &data)
int	loadwavelet (const std::vector< double > &wavelet)
ScalarDecon &	operator= (const ScalarDecon &parent)
std::vector< double >	getresult ()
void	changeparameter (const mspass::utility::Metadata &md)
void	change_shaping_wavelet (const ShapingWavelet &nsw)
ShapingWavelet	get_shaping_wavelet () const
mspass::seismic::CoreTimeSeries	ideal_output ()

Additional Inherited Members
Protected Attributes inherited from mspass::algorithms::deconvolution::FFTDeconOperator
int	nfft
int	sample_shift
gsl_fft_complex_wavetable *	wavetable
gsl_fft_complex_workspace *	workspace
ComplexArray	winv
Protected Attributes inherited from mspass::algorithms::deconvolution::ScalarDecon
std::vector< double >	data
std::vector< double >	wavelet
std::vector< double >	result
ShapingWavelet	shapingwavelet

Constructor & Destructor Documentation

LeastSquareDecon() [1/3]

mspass::algorithms::deconvolution::LeastSquareDecon::LeastSquareDecon ( )

inline

: FFTDeconOperator(), ScalarDecon() { damp = 1.0; };

LeastSquareDecon() [2/3]

mspass::algorithms::deconvolution::LeastSquareDecon::LeastSquareDecon

(

const LeastSquareDecon &

parent

)

    : FFTDeconOperator(parent), ScalarDecon(parent) {
  damp = parent.damp;
}

LeastSquareDecon() [3/3]

mspass::algorithms::deconvolution::LeastSquareDecon::LeastSquareDecon

(

const mspass::utility::Metadata &

md

)

                                                     : FFTDeconOperator(md) {
  try {
    this->read_metadata(md);
  } catch (...) {
    throw;
  }
}

Member Function Documentation

actual_output()

CoreTimeSeries mspass::algorithms::deconvolution::LeastSquareDecon::actual_output ( )

virtual

Return the actual output of the deconvolution operator.

The actual output is defined as w^-1*w and is compable to resolution kernels in linear inverse theory. Although not required we would normally expect this function to be peaked at 0. Offsets from 0 would imply a bias.

Implements mspass::algorithms::deconvolution::ScalarDecon.

                                               {
  try {
    ComplexArray W(nfft, &(wavelet[0]));
    gsl_fft_complex_forward(W.ptr(), 1, nfft, wavetable, workspace);
    ComplexArray ao_fft;
    ao_fft = winv * W;
    /* We always apply the shaping wavelet - this perhaps should be optional
    but probably better done with a none option for the shaping wavelet */
    ao_fft = (*shapingwavelet.wavelet()) * ao_fft;
    gsl_fft_complex_inverse(ao_fft.ptr(), 1, nfft, wavetable, workspace);
    vector<double> ao;
    ao.reserve(nfft);
    for (int k = 0; k < ao_fft.size(); ++k)
      ao.push_back(ao_fft[k].real());
    /* We always shift this wavelet to the center of the data vector.
    We handle the time through the CoreTimeSeries object. */
    int i0 = nfft / 2;
    ao = circular_shift(ao, i0);
    CoreTimeSeries result(nfft);
    /* Getting dt from here is unquestionably a flaw in the api, but will
    retain for now.   Perhaps should a copy of dt in the ScalarDecon object. */
    double dt = this->shapingwavelet.sample_interval();
 
    result.set_t0(-dt * ((double)i0));
    result.set_dt(dt);
    result.set_live();
    result.set_npts(nfft);
    result.set_tref(TimeReferenceType::Relative);
    for (int k = 0; k < nfft; ++k)
      result.s[k] = ao[k];
    result.s = normalize<double>(result.s);
    return result;
  } catch (...) {
    throw;
  };
}

References mspass::seismic::CoreTimeSeries::s, mspass::seismic::CoreTimeSeries::set_dt(), mspass::seismic::BasicTimeSeries::set_live(), mspass::seismic::CoreTimeSeries::set_npts(), mspass::seismic::CoreTimeSeries::set_t0(), and mspass::seismic::BasicTimeSeries::set_tref().

changeparameter()

void mspass::algorithms::deconvolution::LeastSquareDecon::changeparameter ( const mspass::utility::Metadata &  md )

virtual

Implements mspass::algorithms::deconvolution::BasicDeconOperator.

                                                         {
  try {
    this->read_metadata(md);
  } catch (...) {
    throw;
  };
}

inverse_wavelet() [1/2]

CoreTimeSeries mspass::algorithms::deconvolution::LeastSquareDecon::inverse_wavelet ( )

virtual

Return default FIR represesentation of the inverse filter.

This is an overloaded version of the parameterized method. It is equivalent to this->inverse_wavelet(0.0,0.0);

Implements mspass::algorithms::deconvolution::ScalarDecon.

                                                 {
  try {
    return this->inverse_wavelet(0.0);
  } catch (...) {
    throw;
  };
}

inverse_wavelet() [2/2]

CoreTimeSeries mspass::algorithms::deconvolution::LeastSquareDecon::inverse_wavelet ( const double  t0parent = 0.0 )

virtual

Return a FIR respresentation of the inverse filter.

An inverse filter has an impulse response. For some wavelets this can be respresented by a FIR filter with finite numbers of coefficients. Since this is a Fourier method the best we can do is return the inverse fft of the regularized operator. The output usually needs to be phase shifted to be most useful. For typical seismic source wavelets that are approximately minimum phase the shift can be small, but for zero phase input it should be approximately half the window size. This method also has an optional argument for t0parent. Because this processor was written to be agnostic about a time standard it implicitly assumes time 0 is sample 0 of the input waveforms. If the original data have a nonzero start time this should be passed as t0parent or the output will contain a time shift of t0parent. Note that tshift and t0parent do very different things. tshift is used to apply circular phase shift to the output (e.g. a shift of 10 samples causes the last 10 samples in the wavelet to be wrapped to the first 10 samples). t0parent only changes the time standard so the output has t0 -= parent.t0.

The returned wavelet is always time shifted by nfft/2.

Parameters

t0parent

- time zero of parent seismograms (see above).

Implements mspass::algorithms::deconvolution::ScalarDecon.

                                                                      {
  try {
    /* Getting dt from here is unquestionably a flaw in the api, but will
     *         retain for now.   Perhaps should a copy of dt in the ScalarDecon
     * object. */
    double dt = this->shapingwavelet.sample_interval();
    return (this->FourierInverse(this->winv, *shapingwavelet.wavelet(), dt,
                                 t0parent));
  } catch (...) {
    throw;
  };
}

process()

void mspass::algorithms::deconvolution::LeastSquareDecon::process ( )

virtual

Implements mspass::algorithms::deconvolution::ScalarDecon.

                               {
 
  const string base_error("LeastSquareDecon::process:  ");
  // apply fft to the input trace data
  if (data.size() < nfft)
    for (int i = data.size(); i < nfft; ++i)
      data.push_back(0.0);
  ComplexArray d_fft(nfft, &(data[0]));
  gsl_fft_complex_forward(d_fft.ptr(), 1, nfft, wavetable, workspace);
 
  // apply fft to wavelet
  if (wavelet.size() < nfft)
    for (int i = wavelet.size(); i < nfft; ++i)
      wavelet.push_back(0.0);
  ComplexArray b_fft(nfft, &(wavelet[0]));
  gsl_fft_complex_forward(b_fft.ptr(), 1, nfft, wavetable, workspace);
 
  // deconvolution: RF=conj(B).*D./(conj(B).*B+damp)
  b_fft.conj();
  ComplexArray rf_fft, conj_b_fft(b_fft);
  b_fft.conj();
 
  double b_rms = b_fft.rms();
  ComplexArray denom(conj_b_fft * b_fft);
  double theta(b_rms * damp);
  /* This is like normal equation form for damped inverse so theta as computed
  needs to be squared */
  theta = theta * theta;
  for (int k = 0; k < nfft; ++k) {
    double *ptr;
    ptr = denom.ptr(k);
    /* ptr points to the real part - an oddity of this interface */
    *ptr += theta;
  }
  rf_fft = (conj_b_fft * d_fft) / denom;
  // rf_fft=(conj_b_fft*d_fft)/(conj_b_fft*b_fft+b_rms*damp);
  /* Compute the frequency domain version of the inverse wavelet now
  and save it the object for efficiency. */
  // winv=conj_b_fft/(conj_b_fft*b_fft+b_rms*damp);
  winv = conj_b_fft / denom;
 
  // apply shaping wavelet but only to rf estimate - actual output and
  // inverse_wavelet methods apply it to when needed there for efficiency
  ShapingWavelet sw(this->get_shaping_wavelet());
  rf_fft = (*sw.wavelet()) * rf_fft;
 
  // ifft gets result
  gsl_fft_complex_inverse(rf_fft.ptr(), 1, nfft, wavetable, workspace);
  if (sample_shift > 0) {
    for (int k = sample_shift; k > 0; k--)
      result.push_back(rf_fft[nfft - k].real());
    for (int k = 0; k < data.size() - sample_shift; k++)
      result.push_back(rf_fft[k].real());
  } else if (sample_shift == 0) {
    for (int k = 0; k < data.size(); k++)
      result.push_back(rf_fft[k].real());
  } else {
    throw MsPASSError(base_error + "Coding error - trying to use an illegal "
                                   "negative time shift parameter",
                      ErrorSeverity::Fatal);
  }
}

QCMetrics()

Metadata mspass::algorithms::deconvolution::LeastSquareDecon::QCMetrics ( )

virtual

Return appropriate quality measures.

Each operator commonly has different was to measure the quality of the result. This method should return these in a generic Metadata object.

Implements mspass::algorithms::deconvolution::ScalarDecon.

                                     {
  /* Return only an empty Metadata container.  Done as it is
  easier to maintain the code letting python do this work.
  This also anticipates new metrics being added which would be
  easier in python.*/
  Metadata md;
  return md;
}

---

The documentation for this class was generated from the following files:

• /home/runner/work/mspass/mspass/cxx/include/mspass/algorithms/deconvolution/LeastSquareDecon.h
• /home/runner/work/mspass/mspass/cxx/src/lib/algorithms/deconvolution/LeastSquareDecon.cc