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

Scalar water-level deconvolution operator. More...

#include <WaterLevelDecon.h>

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

Collaboration diagram for mspass::algorithms::deconvolution::WaterLevelDecon:

Public Member Functions
	WaterLevelDecon ()
	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 ()
	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	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 water-level deconvolution operator.

This FFT-based operator stabilizes spectral division by flooring the wavelet power spectrum at a configured water level.

Constructor & Destructor Documentation

WaterLevelDecon() [1/4]

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

inline

Construct with default water level and QC state.

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

WaterLevelDecon() [2/4]

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

(

const WaterLevelDecon &

parent

)

Copy constructor.

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

WaterLevelDecon() [3/4]

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

(

const mspass::utility::Metadata &

md

)

Construct from Metadata parameters.

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

WaterLevelDecon() [4/4]

mspass::algorithms::deconvolution::WaterLevelDecon::WaterLevelDecon	(	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::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 {
    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 (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::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::WaterLevelDecon::changeparameter ( const mspass::utility::Metadata &  md )

virtual

Update operator parameters from Metadata.

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

References inverse_wavelet().

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

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::WaterLevelDecon::process ( )

virtual

Compute the water-level inverse operator and deconvolved output.

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

                              {
 
  result.clear();
  const string base_error("WaterLevelDecon::process():  ");
  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
  //  data and wavelet sizes need to be zero padded if the are short
  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);
 
  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 */
  vector<double> d0(nfft, 0.0);
  d0[0] = 1.0;
  ComplexArray delta0(nfft, 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);
  result = ExtractLagWindow(rf_fft, output_length, sample_shift);
}

References mspass::algorithms::deconvolution::ScalarDecon::data, 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::shapingwavelet, 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::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.

                                    {
  Metadata md(this->BasicQCMetrics("WaterLevelDecon", !result.empty()));
  md.put("decon_operator_nfft", nfft);
  md.put("decon_operator_sample_shift", sample_shift);
  md.put("water_level", wlv);
  md.put("water_level_regularization_fraction", regularization_fraction);
  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/WaterLevelDecon.h
• /home/runner/work/mspass/mspass/cxx/src/lib/algorithms/deconvolution/WaterLevelDecon.cc