mspass::algorithms::deconvolution::NoiseStableDecon Class Reference

Noise-aware stable FFT inverse used by NS-GID. More...

#include <NoiseStableDecon.h>

Inheritance diagram for mspass::algorithms::deconvolution::NoiseStableDecon:

Collaboration diagram for mspass::algorithms::deconvolution::NoiseStableDecon:

Public Member Functions
	NoiseStableDecon (const mspass::utility::Metadata &md)
	Construct from metadata defining FFT and NS-GID parameters.
	NoiseStableDecon (const mspass::utility::Metadata &md, const std::vector< double > &wavelet, const std::vector< double > &data)
	Construct and load wavelet and data vectors.
void	changeparameter (const mspass::utility::Metadata &md)
	Update operator parameters from a Metadata container.
void	process ()
	Compute the noise-stable inverse filter and result.
void	loadnoise (const std::vector< double > &noise)
	Load a noise waveform used to form frequency-dependent damping.
void	loadnoise (const mspass::seismic::CoreTimeSeries &noise)
	Load a noise waveform from a core time series.
void	loadnoise (const mspass::seismic::PowerSpectrum &noise_spectrum)
	Load a precomputed noise power spectrum.
mspass::seismic::CoreTimeSeries	actual_output ()
	Return the actual output of the deconvolution operator.
mspass::seismic::CoreTimeSeries	inverse_wavelet (const double t0parent)
	Return the inverse filter impulse response.
mspass::seismic::CoreTimeSeries	inverse_wavelet ()
	Return a FIR represention of the inverse filter.
mspass::utility::Metadata	QCMetrics ()
	Return appropriate quality measures.
double	max_gain () const
double	gain_cap () const
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

Noise-aware stable FFT inverse used by NS-GID.

This scalar operator computes

G(f)=B(f) conj(S(f))/(|S(f)|^2 + mu(f))

with frequency-dependent damping from an optional noise spectrum and an explicit gain cap. The minimum damping parameter is scaled by peak wavelet spectral power. It is intended as the inverse operator used to form the GID peak-picking function. As a standalone scalar deconvolution operator it applies one stable inverse filter and returns the finite-bandwidth result; it does not run sparse iteration and does not apply the GID output shaping wavelet.

Constructor & Destructor Documentation

NoiseStableDecon() [1/3]

mspass::algorithms::deconvolution::NoiseStableDecon::NoiseStableDecon

(

)

    : FFTDeconOperator(), ScalarDecon(), mu_min(3.0e-3), alpha(1.0),
      noise_floor(1.0e-12), gain_max(1.0e3), snr_taper_low(1.0),
      snr_taper_high(3.0), use_reliability_taper(false),
      noise_vector_loaded(false), noise_spectrum_loaded(false),
      processed(false), max_gain_actual(0.0), noise_amplification(0.0),
      effective_bandwidth_fraction(0.0) {}

NoiseStableDecon() [2/3]

mspass::algorithms::deconvolution::NoiseStableDecon::NoiseStableDecon

(

const mspass::utility::Metadata &

md

)

Construct from metadata defining FFT and NS-GID parameters.

Parameters

md	metadata containing FFT sizing, shaping-wavelet, and noise-stable inverse parameters.

    : FFTDeconOperator(md), noise_vector_loaded(false),
      noise_spectrum_loaded(false), processed(false) {
  this->read_metadata(md);
}

NoiseStableDecon() [3/3]

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

Construct and load wavelet and data vectors.

Parameters

md	metadata containing FFT sizing, shaping-wavelet, and noise-stable inverse parameters.
wavelet	source wavelet samples.
data	data samples to deconvolve with the source wavelet.

    : FFTDeconOperator(md), noise_vector_loaded(false),
      noise_spectrum_loaded(false), processed(false) {
  this->read_metadata(md);
  this->wavelet = wavelet;
  this->data = data;
}

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

Member Function Documentation

actual_output()

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

                                               {
  if (!processed || winv.size() <= 0)
    throw MsPASSError(
        "NoiseStableDecon::actual_output: process must be called first",
        ErrorSeverity::Invalid);
  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 = winv * W;
  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());
  int i0 = nfft / 2;
  ao = circular_shift(ao, i0);
  CoreTimeSeries result(nfft);
  double dt = shapingwavelet.sample_interval();
  result.set_t0(-dt * static_cast<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;
}

References mspass::algorithms::deconvolution::FFTDeconOperator::nfft, mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::algorithms::deconvolution::ScalarDecon::result, mspass::algorithms::deconvolution::ShapingWavelet::sample_interval(), mspass::algorithms::deconvolution::ScalarDecon::shapingwavelet, mspass::algorithms::deconvolution::ComplexArray::size(), mspass::algorithms::deconvolution::ScalarDecon::wavelet, mspass::algorithms::deconvolution::FFTDeconOperator::wavetable, mspass::algorithms::deconvolution::FFTDeconOperator::winv, and mspass::algorithms::deconvolution::FFTDeconOperator::workspace.

changeparameter()

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

virtual

Update operator parameters from a Metadata container.

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

                                                         {
  this->read_metadata(md);
  result.clear();
  winv = ComplexArray();
  processed = false;
}

References mspass::algorithms::deconvolution::ScalarDecon::result, and mspass::algorithms::deconvolution::FFTDeconOperator::winv.

gain_cap()

double mspass::algorithms::deconvolution::NoiseStableDecon::gain_cap ( ) const

inline

Return the configured hard cap on inverse-filter gain.

{ return gain_max; }

inverse_wavelet() [1/2]

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

virtual

Return a FIR represention of the inverse filter.

After any deconvolution is computed one can sometimes produce a finite impulse response (FIR) respresentation of the inverse filter.

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

                                                 {
  return this->inverse_wavelet(0.0);
}

References inverse_wavelet().

inverse_wavelet() [2/2]

CoreTimeSeries mspass::algorithms::deconvolution::NoiseStableDecon::inverse_wavelet ( const double  t0parent )

virtual

Return the inverse filter impulse response.

Parameters

t0parent

parent waveform start time used to set the output time standard.

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

                                                                      {
  if (!processed || winv.size() <= 0)
    throw MsPASSError(
        "NoiseStableDecon::inverse_wavelet: process must be called first",
        ErrorSeverity::Invalid);
  double dt = shapingwavelet.sample_interval();
  ComplexArray no_shaping(nfft);
  for (int k = 0; k < nfft; ++k) {
    double *ptr = no_shaping.ptr(k);
    ptr[0] = 1.0;
    ptr[1] = 0.0;
  }
  return this->FourierInverse(winv, no_shaping, dt, t0parent);
}

References mspass::algorithms::deconvolution::FFTDeconOperator::FourierInverse(), mspass::algorithms::deconvolution::FFTDeconOperator::nfft, mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::algorithms::deconvolution::ShapingWavelet::sample_interval(), mspass::algorithms::deconvolution::ScalarDecon::shapingwavelet, mspass::algorithms::deconvolution::ComplexArray::size(), and mspass::algorithms::deconvolution::FFTDeconOperator::winv.

loadnoise() [1/3]

void mspass::algorithms::deconvolution::NoiseStableDecon::loadnoise

(

const mspass::seismic::CoreTimeSeries &

noise

)

Load a noise waveform from a core time series.

Parameters

noise

scalar noise time series.

                                                               {
  this->loadnoise(noise_in.s);
}

References loadnoise(), and mspass::seismic::CoreTimeSeries::s.

loadnoise() [2/3]

void mspass::algorithms::deconvolution::NoiseStableDecon::loadnoise

(

const mspass::seismic::PowerSpectrum &

noise_spectrum

)

Load a precomputed noise power spectrum.

Parameters

noise_spectrum

spectrum that covers DC and describes noise power.

                                                                       {
  ValidatePowerSpectrumCoversDC(noise_spectrum_in,
                                "NoiseStableDecon::loadnoise");
  noise_spectrum = noise_spectrum_in;
  noise_spectrum_loaded = true;
  noise_vector_loaded = false;
  result.clear();
  winv = ComplexArray();
  processed = false;
}

References mspass::algorithms::deconvolution::ScalarDecon::result, and mspass::algorithms::deconvolution::FFTDeconOperator::winv.

loadnoise() [3/3]

void mspass::algorithms::deconvolution::NoiseStableDecon::loadnoise

(

const std::vector< double > &

noise

)

Load a noise waveform used to form frequency-dependent damping.

Parameters

noise

scalar noise samples.

                                                               {
  if (noise.empty())
    throw MsPASSError("NoiseStableDecon::loadnoise: noise vector cannot be "
                      "empty",
                      ErrorSeverity::Invalid);
  this->noise = noise;
  noise_vector_loaded = true;
  noise_spectrum_loaded = false;
  result.clear();
  winv = ComplexArray();
  processed = false;
}

References mspass::algorithms::deconvolution::ScalarDecon::result, and mspass::algorithms::deconvolution::FFTDeconOperator::winv.

max_gain()

double mspass::algorithms::deconvolution::NoiseStableDecon::max_gain ( ) const

inline

Return the largest gain reached by the processed inverse filter.

{ return max_gain_actual; }

process()

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

virtual

Compute the noise-stable inverse filter and result.

The method requires loaded wavelet and data vectors and rebuilds cached output, inverse-filter, and QC state.

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

                               {
  const string base_error("NoiseStableDecon::process: ");
  result.clear();
  winv = ComplexArray();
  processed = false;
  const int output_length = data.size();
  if (output_length <= 0)
    throw MsPASSError(base_error + "no data loaded", ErrorSeverity::Invalid);
  if (wavelet.empty())
    throw MsPASSError(base_error + "no wavelet loaded", ErrorSeverity::Invalid);
  if (output_length > nfft || wavelet.size() > nfft)
    throw MsPASSError(base_error + "loaded vectors exceed padded FFT length",
                      ErrorSeverity::Invalid);
 
  vector<double> data_padded(data);
  data_padded.resize(nfft, 0.0);
  ComplexArray d_fft(nfft, &(data_padded[0]));
  gsl_fft_complex_forward(d_fft.ptr(), 1, nfft, wavetable, workspace);
 
  vector<double> wavelet_padded(wavelet);
  wavelet_padded.resize(nfft, 0.0);
  ComplexArray s_fft(nfft, &(wavelet_padded[0]));
  gsl_fft_complex_forward(s_fft.ptr(), 1, nfft, wavetable, workspace);
  if (s_fft.rms() <= 0.0)
    throw MsPASSError(base_error + "wavelet vector is all zeros",
                      ErrorSeverity::Invalid);
 
  const double dt = shapingwavelet.sample_interval();
  vector<double> pn = this->noise_power_spectrum(dt);
  winv = ComplexArray(nfft);
  max_gain_actual = 0.0;
  noise_amplification = 0.0;
  double spectral_power_peak = 0.0;
  for (int k = 0; k < nfft; ++k)
    spectral_power_peak = max(spectral_power_peak, norm(s_fft[k]));
  const double mu_floor = mu_min * spectral_power_peak;
  int usable_bins = 0;
  for (int k = 0; k < nfft; ++k) {
    Complex64 s = s_fft[k];
    double absS2 = norm(s);
    double absS = sqrt(absS2);
    double snr = absS2 / (pn[k] + noise_floor);
    double mu_noise = alpha * absS2 / (snr + 1.0e-12);
    double mu_gain = max(0.0, absS / gain_max - absS2);
    double mu = max({mu_floor, mu_noise, mu_gain});
    double b = this->reliability_taper(snr);
    Complex64 g(0.0, 0.0);
    if ((absS2 + mu) > 0.0)
      g = b * conj(s) / (absS2 + mu);
    double gain = abs(g);
    if (gain > gain_max) {
      g *= gain_max / gain;
      gain = gain_max;
    }
    double *ptr = winv.ptr(k);
    ptr[0] = g.real();
    ptr[1] = g.imag();
    max_gain_actual = max(max_gain_actual, gain);
    noise_amplification += gain * gain * pn[k];
    if (b > 0.0 && gain > 0.0)
      ++usable_bins;
  }
  noise_amplification = sqrt(noise_amplification / static_cast<double>(nfft));
  effective_bandwidth_fraction =
      static_cast<double>(usable_bins) / static_cast<double>(nfft);
 
  ComplexArray rf_fft = winv * d_fft;
  gsl_fft_complex_inverse(rf_fft.ptr(), 1, nfft, wavetable, workspace);
  result = ExtractLagWindow(rf_fft, output_length, sample_shift);
  processed = true;
}

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::ShapingWavelet::sample_interval(), mspass::algorithms::deconvolution::FFTDeconOperator::sample_shift, mspass::algorithms::deconvolution::ScalarDecon::shapingwavelet, mspass::algorithms::deconvolution::ScalarDecon::wavelet, mspass::algorithms::deconvolution::FFTDeconOperator::wavetable, mspass::algorithms::deconvolution::FFTDeconOperator::winv, and mspass::algorithms::deconvolution::FFTDeconOperator::workspace.

QCMetrics()

Metadata mspass::algorithms::deconvolution::NoiseStableDecon::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("NoiseStableDecon", processed));
  md.put("decon_operator_nfft", nfft);
  md.put("decon_operator_sample_shift", sample_shift);
  md.put("ns_gid_gain_max_requested", gain_max);
  md.put("ns_gid_gain_max_actual", max_gain_actual);
  md.put("ns_gid_mu_min", mu_min);
  md.put("ns_gid_alpha", alpha);
  md.put("ns_gid_noise_amplification", noise_amplification);
  md.put("ns_gid_effective_bandwidth_fraction",
         effective_bandwidth_fraction);
  md.put("ns_gid_operator_nfft", nfft);
  md.put("ns_gid_use_reliability_taper", use_reliability_taper);
  return md;
}

References mspass::algorithms::deconvolution::ScalarDecon::BasicQCMetrics(), mspass::algorithms::deconvolution::FFTDeconOperator::nfft, mspass::utility::Metadata::put(), 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/NoiseStableDecon.h
• /home/runner/work/mspass/mspass/cxx/src/lib/algorithms/deconvolution/NoiseStableDecon.cc