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

Scalar least-squares deconvolution operator. More...

#include <LeastSquareDecon.h>

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

Collaboration diagram for mspass::algorithms::deconvolution::LeastSquareDecon:

Public Member Functions
	LeastSquareDecon ()
	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 ()
	Construct an empty FFT operator shell.
	FFTDeconOperator (const mspass::utility::Metadata &md)
	Construct from deconvolution Metadata.
	FFTDeconOperator (const FFTDeconOperator &parent)
	Copy constructor.
	~FFTDeconOperator ()
	Release allocated GSL FFT work objects.
FFTDeconOperator &	operator= (const FFTDeconOperator &parent)
	Assign FFT size and lag shift from another operator.
void	changeparameter (const mspass::utility::Metadata &md)
	Reconfigure the FFT operator from Metadata.
void	change_size (const int nfft_new)
	Change the FFT work-buffer length.
void	change_shift (const int shift)
	Set the sample offset used to unwrap deconvolution lag windows.
int	get_size ()
	Return the current FFT work-buffer length.
int	get_shift ()
	Return the current deconvolution lag-window sample shift.
int	operator_size ()
	Return the current FFT work-buffer length.
int	operator_shift ()
	Return the current deconvolution lag-window sample shift.
double	df (const double dt)
	Return frequency spacing for the current FFT length.
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 ()
	Construct an empty scalar deconvolution operator.
	ScalarDecon (const mspass::utility::Metadata &md)
	Construct from Metadata.
	ScalarDecon (const std::vector< double > &d, const std::vector< double > &w)
	Construct with loaded data and source-wavelet vectors.
	ScalarDecon (const ScalarDecon &parent)
	Copy constructor preserving loaded vectors and shaping wavelet.
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)
	Assign loaded wavelet, data, and result vectors from parent.
std::vector< double >	getresult ()
	Return the current deconvolution result vector.
void	changeparameter (const mspass::utility::Metadata &md)
void	change_shaping_wavelet (const ShapingWavelet &nsw)
ShapingWavelet	get_shaping_wavelet () const
virtual mspass::seismic::CoreTimeSeries	output_shaping_wavelet ()
	Return the output shaping wavelet.
mspass::seismic::CoreTimeSeries	ideal_output ()
	Legacy alias for output_shaping_wavelet.
mspass::seismic::CoreTimeSeries	resolution_kernel ()
	Alias for actual_output using inverse-theory terminology.
Public Member Functions inherited from mspass::algorithms::deconvolution::BasicDeconOperator
virtual	~BasicDeconOperator ()

Additional Inherited Members
Protected Member Functions inherited from mspass::algorithms::deconvolution::ScalarDecon
mspass::utility::Metadata	BasicQCMetrics (const std::string &operator_name, const bool processed)
	Build QC Metadata fields common to scalar deconvolution operators.
Protected Attributes inherited from mspass::algorithms::deconvolution::FFTDeconOperator
int	nfft
	Current FFT work-buffer length in samples.
int	sample_shift
	Sample offset used to place zero lag in deconvolution outputs.
gsl_fft_complex_wavetable *	wavetable
	GSL factorization table allocated for nfft.
gsl_fft_complex_workspace *	workspace
	GSL scratch workspace allocated for nfft.
ComplexArray	winv
	Frequency-domain inverse wavelet coefficients for derived operators.
Protected Attributes inherited from mspass::algorithms::deconvolution::ScalarDecon
std::vector< double >	data
	Data vector to be deconvolved by concrete scalar methods.
std::vector< double >	wavelet
	Source-wavelet estimate used by concrete scalar methods.
std::vector< double >	result
	Deconvolved output vector produced by process.
ShapingWavelet	shapingwavelet
	Output shaping wavelet applied by scalar deconvolution methods.

Detailed Description

Scalar least-squares deconvolution operator.

This FFT-based operator estimates a regularized inverse filter from an input wavelet and applies it to scalar data.

Constructor & Destructor Documentation

LeastSquareDecon() [1/4]

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

inline

Construct with default damping.

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

LeastSquareDecon() [2/4]

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

(

const LeastSquareDecon &

parent

)

Copy constructor.

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

LeastSquareDecon() [3/4]

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

(

const mspass::utility::Metadata &

md

)

Construct from Metadata parameters.

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

LeastSquareDecon() [4/4]

mspass::algorithms::deconvolution::LeastSquareDecon::LeastSquareDecon	(	const mspass::utility::Metadata &	md,
		const std::vector< double > &	wavelet,
		const std::vector< double > &	data
	)

Construct from Metadata plus wavelet and data vectors.

                                {
  try {
    this->read_metadata(md);
  } catch (...) {
    throw;
  };
  wavelet = w;
  data = d;
}

References mspass::algorithms::deconvolution::ScalarDecon::data, and mspass::algorithms::deconvolution::ScalarDecon::wavelet.

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 {
    vector<double> wavelet_padded(wavelet);
    if (wavelet_padded.size() < nfft)
      wavelet_padded.resize(nfft, 0.0);
    ComplexArray W(nfft, &(wavelet_padded[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::algorithms::deconvolution::FFTDeconOperator::nfft, mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::algorithms::deconvolution::ScalarDecon::result, mspass::algorithms::deconvolution::ShapingWavelet::sample_interval(), mspass::seismic::CoreTimeSeries::set_t0(), mspass::algorithms::deconvolution::ScalarDecon::shapingwavelet, mspass::algorithms::deconvolution::ComplexArray::size(), mspass::algorithms::deconvolution::ScalarDecon::wavelet, mspass::algorithms::deconvolution::ShapingWavelet::wavelet(), mspass::algorithms::deconvolution::FFTDeconOperator::wavetable, mspass::algorithms::deconvolution::FFTDeconOperator::winv, and mspass::algorithms::deconvolution::FFTDeconOperator::workspace.

changeparameter()

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

virtual

Update operator parameters from Metadata.

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;
  };
}

References inverse_wavelet().

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;
  };
}

References mspass::algorithms::deconvolution::FFTDeconOperator::FourierInverse(), mspass::algorithms::deconvolution::ShapingWavelet::sample_interval(), mspass::algorithms::deconvolution::ScalarDecon::shapingwavelet, mspass::algorithms::deconvolution::ShapingWavelet::wavelet(), and mspass::algorithms::deconvolution::FFTDeconOperator::winv.

process()

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

virtual

Compute the inverse operator and deconvolved output.

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

                               {
 
  const string base_error("LeastSquareDecon::process:  ");
  result.clear();
  const int output_length = data.size();
  if (output_length > nfft)
    throw MsPASSError(base_error + "data vector exceeds padded fft length",
                      ErrorSeverity::Invalid);
  // apply fft to the input trace data
  vector<double> data_padded(data);
  if (data_padded.size() < nfft)
    data_padded.resize(nfft, 0.0);
  ComplexArray d_fft(nfft, &(data_padded[0]));
  gsl_fft_complex_forward(d_fft.ptr(), 1, nfft, wavetable, workspace);
 
  // apply fft to wavelet
  if (wavelet.size() > nfft)
    throw MsPASSError(base_error + "wavelet vector exceeds padded fft length",
                      ErrorSeverity::Invalid);
  vector<double> wavelet_padded(wavelet);
  if (wavelet_padded.size() < nfft)
    wavelet_padded.resize(nfft, 0.0);
  ComplexArray b_fft(nfft, &(wavelet_padded[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();
  if (b_rms == 0.0)
    throw MsPASSError(base_error + "wavelet data vector is all zeros",
                      ErrorSeverity::Invalid);
  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);
  result = ExtractLagWindow(rf_fft, output_length, sample_shift);
}

References mspass::algorithms::deconvolution::ComplexArray::conj(), mspass::algorithms::deconvolution::ScalarDecon::data, mspass::algorithms::deconvolution::ScalarDecon::get_shaping_wavelet(), mspass::algorithms::deconvolution::FFTDeconOperator::nfft, mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::algorithms::deconvolution::ScalarDecon::result, mspass::algorithms::deconvolution::ComplexArray::rms(), mspass::algorithms::deconvolution::FFTDeconOperator::sample_shift, mspass::algorithms::deconvolution::ScalarDecon::wavelet, mspass::algorithms::deconvolution::ShapingWavelet::wavelet(), mspass::algorithms::deconvolution::FFTDeconOperator::wavetable, mspass::algorithms::deconvolution::FFTDeconOperator::winv, and mspass::algorithms::deconvolution::FFTDeconOperator::workspace.

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.

                                     {
  Metadata md(this->BasicQCMetrics("LeastSquareDecon", !result.empty()));
  md.put("decon_operator_nfft", nfft);
  md.put("decon_operator_sample_shift", sample_shift);
  md.put("damping_factor", damp);
  md.put("least_square_damping_factor", damp);
  return md;
}

References mspass::algorithms::deconvolution::ScalarDecon::BasicQCMetrics(), mspass::algorithms::deconvolution::FFTDeconOperator::nfft, mspass::utility::Metadata::put(), mspass::algorithms::deconvolution::ScalarDecon::result, and mspass::algorithms::deconvolution::FFTDeconOperator::sample_shift.

---

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