mspass::algorithms::deconvolution::WaterLevelDecon Class Reference

Inheritance diagram for mspass::algorithms::deconvolution::WaterLevelDecon:

Public Member Functions
	WaterLevelDecon (const WaterLevelDecon &parent)
	WaterLevelDecon (const mspass::utility::Metadata &md)
	WaterLevelDecon (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	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

WaterLevelDecon() [1/3]

mspass::algorithms::deconvolution::WaterLevelDecon::WaterLevelDecon ( )

inline

                    : FFTDeconOperator(), ScalarDecon() {
    this->wlv = 0.1;
    this->regularization_fraction = 0.0;
  };

WaterLevelDecon() [2/3]

mspass::algorithms::deconvolution::WaterLevelDecon::WaterLevelDecon

(

const WaterLevelDecon &

parent

)

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

WaterLevelDecon() [3/3]

mspass::algorithms::deconvolution::WaterLevelDecon::WaterLevelDecon

(

const mspass::utility::Metadata &

md

)

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

Member Function Documentation

actual_output()

CoreTimeSeries mspass::algorithms::deconvolution::WaterLevelDecon::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 (unsigned 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);
    ao = normalize<double>(ao);
    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_tref(TimeReferenceType::Relative);
    result.set_npts(nfft);
    for (int k = 0; k < nfft; ++k)
      result.s[k] = ao[k];
    ;
    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::WaterLevelDecon::changeparameter ( const mspass::utility::Metadata &  md )

virtual

Reimplemented from mspass::algorithms::deconvolution::ScalarDecon.

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

inverse_wavelet() [1/2]

CoreTimeSeries mspass::algorithms::deconvolution::WaterLevelDecon::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::WaterLevelDecon::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.

Output wavelet is always circular shifted with 0 lag at center.

Parameters

t0parent

- time zero of parent seismograms (see above).

Exceptions

SeisppError

is thrown if the tshift value is more than half the length of the data (nfft*dt). Reason is justified 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->FFTDeconOperator::FourierInverse(
        this->winv, *shapingwavelet.wavelet(), dt, t0parent));
  } catch (...) {
    throw;
  };
}

process()

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

virtual

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

                              {
 
  // apply fft to the input trace data
  //  data and wavelet sizes need to be zero padded if the are short
  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);
 
  double b_rms = b_fft.rms();
  if (b_rms == 0.0)
    throw MsPASSError(
        "WaterLevelDecon::process():  wavelet data vector is all zeros");
 
  // water level - count the number of points below water level
  int nunderwater(0);
  for (int i = 0; i < nfft; i++) {
    /* We have to avoid hard zeros because the formula below will
     * yield an NaN when that happens from 0/0 */
    double bamp;
    bamp = abs(b_fft[i]);
    /*WARNING:  this is dependent on implementation detail of ComplexArray
    that defines the data as a FortranComplex64 - real,imag pairs. */
    if (bamp < b_rms * wlv) {
      /* Use 64 bit epsilon as the floor for defining a pure zero */
      if (bamp / b_rms < DBL_EPSILON) {
        *b_fft.ptr(i) = b_rms * wlv;
        *(b_fft.ptr(i) + 1) = b_rms * wlv;
      } else {
        // real part
        *b_fft.ptr(i) = (*b_fft.ptr(i) / abs(b_fft[i])) * b_rms * wlv;
        // imag part
        *(b_fft.ptr(i) + 1) =
            (*(b_fft.ptr(i) + 1) / abs(b_fft[i])) * b_rms * wlv;
      }
      ++nunderwater;
    }
  }
  regularization_fraction = ((double)nunderwater) / ((double)nfft);
  ComplexArray rf_fft;
  rf_fft = d_fft / b_fft;
  /* Make numerator for inverse from zero lag spike */
  double *d0 = new double[nfft];
  for (int k = 0; k < nfft; ++k)
    d0[k] = 0.0;
  d0[0] = 1.0;
  ComplexArray delta0(nfft, d0);
  delete[] d0;
  gsl_fft_complex_forward(delta0.ptr(), 1, nfft, wavetable, workspace);
  winv = delta0 / b_fft;
 
  // apply shaping wavelet to rf estimate
  rf_fft = (*shapingwavelet.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 (unsigned int k = 0; k < data.size() - sample_shift; k++)
      result.push_back(rf_fft[k].real());
  } else {
    for (unsigned int k = 0; k < data.size(); k++)
      result.push_back(rf_fft[k].real());
  }
}

QCMetrics()

Metadata mspass::algorithms::deconvolution::WaterLevelDecon::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/WaterLevelDecon.h
• /home/runner/work/mspass/mspass/cxx/src/lib/algorithms/deconvolution/WaterLevelDecon.cc