mspass::algorithms::deconvolution::ShapingWavelet Class Reference

Frequency domain shaping wavelet. More...

#include <ShapingWavelet.h>

Public Member Functions
	ShapingWavelet (const mspass::utility::Metadata &md, int npts=0)
	Construct using a limited set of analytic forms for the wavelet.
	ShapingWavelet (mspass::seismic::CoreTimeSeries d, int nfft=0)
	Use a wavelet defined by a TimeSeries object.
	ShapingWavelet (const double fpeak, const double dtin, const int n)
	ShapingWavelet (const int npolelo, const double f3dblo, const int npolehi, const double f3dbhi, const double dtin, const int n)
	ShapingWavelet (const ShapingWavelet &parent)
ShapingWavelet &	operator= (const ShapingWavelet &parent)
void	swap (ShapingWavelet &other) noexcept
ComplexArray *	wavelet ()
mspass::seismic::CoreTimeSeries	impulse_response ()
double	freq_bin_size ()
double	sample_interval ()
std::string	type ()
int	size () const

Detailed Description

Frequency domain shaping wavelet.

Frequency domain based deconvolution methods all use a shaping wavelet for output to avoid ringing. Frequency domain deconvolution methods in this Library contain an instance of this object. In all cases it is hidden behind the interface. A complexity, however, is that all frequency domain methods will call the Metadata driven constructor.

This version currently allows three shaping wavelets: Gaussin, Ricker, and Slepian0. The first two are standard. The last is novel and theoretically can produce an actual output with the smalle posible sidebands

Constructor & Destructor Documentation

ShapingWavelet() [1/6]

mspass::algorithms::deconvolution::ShapingWavelet::ShapingWavelet ( )

inline

: nfft(0), dt(-1), df(-1) {};

ShapingWavelet() [2/6]

mspass::algorithms::deconvolution::ShapingWavelet::ShapingWavelet	(	const mspass::utility::Metadata &	md,
		int	npts = 0
	)

Construct using a limited set of analytic forms for the wavelet.

This constructor is used to create a ricker or gaussian shaping wavelet with parameters defined by parameters passed through the Metadata object.

Parameters

md	- Metadata object with parameters specifying the wavelet.
npts	length of signal to be generated which is the same as the fft size for real valued signals. If set 0 (the default) the constructor will attempt to get npts from md using the keyword "operator_nfft".

                                                             {
  const string base_error("ShapingWavelet object constructor:  ");
  try {
    /* We use this to allow a nfft to be set in md.   A bit error prone
     * we add a special error handler. */
    if (nfftin > 0)
      this->nfft = nfftin;
    else {
      try {
        nfft = GetIntRequired(md, "operator_nfft");
      } catch (MetadataGetError &mderr) {
        throw MsPASSError(base_error +
                              "Called constructor with nfft=0 but parameter "
                              "with key=operator_nfft is not in parameter file",
                          ErrorSeverity::Invalid);
      }
    }
    if (nfft <= 0)
      throw MsPASSError(base_error + "FFT length must be positive",
                        ErrorSeverity::Invalid);
    /* these are workspaces used by gnu's fft algorithm.  Other parts
     * of this library cache them for efficiency, but this one we
     * compute ever time the object is created and then discard it. */
    auto resources = detail::AllocateGSLFFTResources(nfft, base_error);
    string wavelettype = md.get_string("shaping_wavelet_type");
    wavelet_name = wavelettype;
    dt = GetDoubleRequired(md, "shaping_wavelet_dt");
    if (wavelettype == "gaussian") {
      float fpeak = GetDoubleRequired(md, "shaping_wavelet_frequency");
      // construct wavelet and fft
      unique_ptr<double[]> r(gaussian(fpeak, (float)dt, nfft));
      if (!r)
        throw bad_alloc();
      w = ComplexArray(nfft, r.get());
      gsl_fft_complex_forward(w.ptr(), 1, nfft, resources.wavetable.get(),
                              resources.workspace.get());
    }
    /* Note for CNR3CDecon the initial values on construction for
    ricker or butterworth are irrelevant and wasted effort.  We keep
    them in this class because these forms are needed by the family of
    scalar deconvolution algorithms */
    else if (wavelettype == "ricker") {
      float fpeak =
          static_cast<float>(GetDoubleRequired(md, "shaping_wavelet_frequency"));
      // construct wavelet and fft
      unique_ptr<double[]> r(rickerwavelet(fpeak, (float)dt, nfft));
      if (!r)
        throw bad_alloc();
      // DEBUG
      // cerr << "Ricker shaping wavelet"<<endl;
      // for(int k=0;k<nfft;++k) cerr << r[k]<<endl;
      w = ComplexArray(nfft, r.get());
      gsl_fft_complex_forward(w.ptr(), 1, nfft, resources.wavetable.get(),
                              resources.workspace.get());
    } else if (wavelettype == "butterworth") {
      double f3db_lo, f3db_hi;
      f3db_lo = GetDoubleRequired(md, "f3db_lo");
      f3db_hi = GetDoubleRequired(md, "f3db_hi");
      int npoles_lo, npoles_hi;
      npoles_lo = GetIntRequired(md, "npoles_lo");
      npoles_hi = GetIntRequired(md, "npoles_hi");
      Butterworth bwf(true, true, true, npoles_lo, f3db_lo, npoles_hi, f3db_hi,
                      this->dt);
      w = bwf.transfer_function(this->nfft);
    }
    /*   This option requires a package to compute zero phase wavelets of
    some specified type and bandwidth.  Previous used antelope filters which
    colides with open source goal of mspass.   Another implementation of
    TimeInvariantFilter and this could be restored.
    */
    /*
    else if(wavelettype=="filter_response")
    {
        string filterparam=md.get_string("filter");
        TimeInvariantFilter filt(filterparam);
        TimeSeries dtmp(nfft);
        dtmp.ns=nfft;
        dtmp.dt=dt;
        dtmp.live=true;
        dtmp.t0=(-dt*static_cast<double>(nfft/2));
        dtmp.tref=relative;
        int i;
        i=dtmp.sample_number(0.0);
        dtmp.s[i]=1.0;  // Delta function at 0
        filt.zerophase(dtmp);
        // We need to shift the filter response now back to
       // zero to avoid time shifts in output
        dtmp.s=circular_shift(dtmp.s,nfft/2);
        w=ComplexArray(dtmp.s.size(), &(dtmp.s[0]));
        gsl_fft_complex_forward(w.ptr(), 1, nfft, wavetable, workspace);
    }
    */
    else if ((wavelettype == "slepian") || (wavelettype == "Slepian")) {
      double tbp = GetDoubleRequired(md, "time_bandwidth_product");
      double target_pulse_width =
          GetDoubleRequired(md, "target_pulse_width");
      /* Sanity check on pulse width */
      if (target_pulse_width > (nfft / 4) || (target_pulse_width < tbp)) {
        stringstream ss;
        ss << "ShapingWavelet Metadata constructor:  bad input parameters for "
              "Slepian shaping wavelet"
           << endl
           << "Specified target_pulse_width of " << target_pulse_width
           << " samples" << endl
           << "FFT buffer size=" << nfft
           << " and time bandwidth product=" << tbp << endl
           << "Pulse width should be small fraction of buffer size but ge than "
              "tbp"
           << endl;
        throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
      }
      double c = tbp / target_pulse_width;
      int nwsize = round(c * (static_cast<double>(nfft)));
      unique_ptr<double[]> wtmp(slepian0(tbp, nwsize));
      if (!wtmp)
        throw bad_alloc();
      vector<double> work(nfft, 0.0);
      dcopy(nwsize, wtmp.get(), 1, &(work[0]), 1);
      w = ComplexArray(nfft, work);
      gsl_fft_complex_forward(w.ptr(), 1, nfft, resources.wavetable.get(),
                              resources.workspace.get());
    } else if (wavelettype == "none") {
      /* The shaping wavelet is stored in the frequency domain.  The identity
       * filter is therefore one at every frequency bin, not a time-domain
       * delta stored without an FFT. */
      w = ComplexArray(nfft, 1.0);
    } else {
      throw MsPASSError(
          base_error + "illegal value for shaping_wavelet_type=" + wavelettype,
          ErrorSeverity::Invalid);
    }
    df = 1.0 / (dt * ((double)nfft));
  } catch (MsPASSError &err) {
    cerr << base_error
         << "Something threw an unhandled MsPASSError with this message:"
         << endl;
    err.log_error();
    cerr << "This is a nonfatal bug that needs to be fixed" << endl;
    throw err;
  } catch (...) {
    throw;
  }
}

References mspass::utility::Metadata::get_string(), mspass::utility::MsPASSError::log_error(), mspass::algorithms::deconvolution::ComplexArray::ptr(), and mspass::algorithms::Butterworth::transfer_function().

ShapingWavelet() [3/6]

mspass::algorithms::deconvolution::ShapingWavelet::ShapingWavelet	(	mspass::seismic::CoreTimeSeries	d,
		int	nfft = 0
	)

Use a wavelet defined by a TimeSeries object.

This constructor uses the data stored in a TimeSeries object to define the shaping wavelet. Note to assure output is properly time aligned this signal is forced to be symmetric around relative time 0. That is necessary because of the way this, like every fft we are aware of, handles phase. Hence, the user should assure zero time is an appropriate zero reference (e.g. a zero phase filter should the dominant peak at 0 time.) If the input data size is smaller than the buffer size specified the buffer is zero padded.

Parameters

d	TimeSeries specifying wavelet. Note dt will be extracted and stored in this object.
nfft	- buffer size = fft length for frequency domain representation of the wavelet. If ns of d is less than nfft or the time range defined by d is does not exceed -(nfft/2)*dt to (nfft2)*dt the result will be sero padded before computing the fft. if nfft is 0 the d.ns will set the fft length.

                                                         {
  const string base_error("ShapingWavelet TimeSeries constructor:  ");
  wavelet_name = string("data");
  /* Silently handle this default condition that allows calling the
   * constructor without the nfft argument */
  if (nfft <= 0)
    nfft = d.npts();
  if (nfft <= 0)
    throw MsPASSError(base_error + "FFT length must be positive",
                      ErrorSeverity::Invalid);
  this->nfft = nfft;
  dt = d.dt();
  df = 1.0 / (dt * ((double)nfft));
  /* This is prone to an off by one error */
  double t0;
  if (nfft % 2)
    t0 = -(double)(nfft / 2) + 1;
  else
    t0 = -(double)(nfft / 2);
  t0 *= dt;
  /* A basic sanity check on d */
  if (d.time_is_UTC())
    throw MsPASSError(
        base_error +
            "Shaping wavelet must be defined in relative time centered on 0",
        ErrorSeverity::Invalid);
  if ((d.endtime() < t0) || (d.t0() > (-t0)))
    throw MsPASSError(base_error + "Input wavelet time span is illegal\n" +
                          "Wavelet must be centered on 0 reference time",
                      ErrorSeverity::Invalid);
  /* Create a work vector an initialize it to zeros to make insertion
   * algorithm easier */
  vector<double> dwork;
  dwork.reserve(nfft);
  int i;
  for (i = 0; i < nfft; ++i)
    dwork.push_back(0.0);
  /* We loop over
   * we skip through d vector until we insert */
  for (i = 0; i < d.npts(); ++i) {
    double t;
    t = t0 + dt * ((double)i);
    int iw = d.sample_number(t);
    if (t > d.endtime())
      break;
    if ((iw >= 0) && (iw < nfft))
      dwork[i] = d.s[iw];
  }
  auto resources = detail::AllocateGSLFFTResources(nfft, base_error);
  w = ComplexArray(nfft, &(dwork[0]));
  gsl_fft_complex_forward(w.ptr(), 1, nfft, resources.wavetable.get(),
                          resources.workspace.get());
}

References mspass::seismic::BasicTimeSeries::dt(), mspass::seismic::BasicTimeSeries::endtime(), mspass::seismic::BasicTimeSeries::npts(), mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::seismic::CoreTimeSeries::s, mspass::seismic::BasicTimeSeries::sample_number(), mspass::seismic::BasicTimeSeries::t0(), and mspass::seismic::BasicTimeSeries::time_is_UTC().

ShapingWavelet() [4/6]

mspass::algorithms::deconvolution::ShapingWavelet::ShapingWavelet	(	const double	fpeak,
		const double	dtin,
		const int	n
	)

Construct a Ricker wavelet shaping filter with peak frequency fpeak.

    : wavelet_name("ricker") {
  const string caller("ShapingWavelet Ricker constructor");
  if (n <= 0)
    throw MsPASSError(caller + ": FFT length must be positive",
                      ErrorSeverity::Invalid);
  nfft = n;
  dt = dtin;
  df = 1.0 / (dt * static_cast<double>(n));
  unique_ptr<double[]> r(rickerwavelet((float)fpeak, (float)dt, nfft));
  if (!r)
    throw bad_alloc();
  w = ComplexArray(nfft, r.get());
  auto resources = detail::AllocateGSLFFTResources(nfft, caller);
  gsl_fft_complex_forward(w.ptr(), 1, nfft, resources.wavetable.get(),
                          resources.workspace.get());
}

References mspass::algorithms::deconvolution::ComplexArray::ptr().

ShapingWavelet() [5/6]

mspass::algorithms::deconvolution::ShapingWavelet::ShapingWavelet	(	const int	npolelo,
		const double	f3dblo,
		const int	npolehi,
		const double	f3dbhi,
		const double	dtin,
		const int	n
	)

Construct a zero phase Butterworth filter wavelet.

This is the recommended constructor to use for adjustable bandwidth shaping. It is the default for CNR3CDecon.

Parameters

npolelo	is the number of poles for the low corner
f3dblo	is the 3db point for the low corner of the passband
npolehi	is the number of poles for the upper corner filter
f3dbhi	is the 3db point for the high corner of the passband.
dtin	sample interval of the wavelet.
n	number of samples in the wavelet.

    : wavelet_name("butterworth") {
  if (n <= 0)
    throw MsPASSError(
        "ShapingWavelet Butterworth constructor: FFT length must be positive",
        ErrorSeverity::Invalid);
  nfft = n;
  dt = dtin;
  df = 1.0 / (dt * static_cast<double>(n));
  Butterworth bwf(true, true, true, npolelo, f3dblo, npolehi, f3dbhi, dtin);
  w = bwf.transfer_function(nfft);
}

References mspass::algorithms::Butterworth::transfer_function().

ShapingWavelet() [6/6]

mspass::algorithms::deconvolution::ShapingWavelet::ShapingWavelet

(

const ShapingWavelet &

parent

)

Copy constructor.

    : nfft(parent.nfft), w(parent.w), dt(parent.dt), df(parent.df),
      wavelet_name(parent.wavelet_name) {}

Member Function Documentation

freq_bin_size()

double mspass::algorithms::deconvolution::ShapingWavelet::freq_bin_size ( )

inline

Return the frequency bin size.

{ return df; };

impulse_response()

CoreTimeSeries mspass::algorithms::deconvolution::ShapingWavelet::impulse_response

(

)

Return the impulse response of the shaping filter. Expect the result to be symmetric about 0 (i.e. output spans nfft/2 to nfft/2.

                                                {
  try {
    int nfft = w.size();
    const string caller("ShapingWavelet::impulse_response");
    if (nfft <= 0)
      throw MsPASSError(caller + ": shaping wavelet is empty",
                        ErrorSeverity::Invalid);
    auto resources = detail::AllocateGSLFFTResources(nfft, caller);
    /* We need to copy the current shaping wavelet or the inverse fft
     * will make it invalid */
    ComplexArray iwf(w);
    gsl_fft_complex_inverse(iwf.ptr(), 1, nfft, resources.wavetable.get(),
                            resources.workspace.get());
    CoreTimeSeries result(nfft);
 
    result.set_tref(TimeReferenceType::Relative);
    result.set_npts(nfft);
    result.set_dt(dt);
    result.set_t0(dt * (-(double)nfft / 2));
    result.set_live();
    /* Unfold the fft output */
    int k, shift;
    shift = nfft / 2;
    // Need this because constructor fills initially with nfft zeros and
    // we use this approach to unfold the fft output
    result.s.clear();
    for (k = 0; k < nfft; ++k)
      result.s.push_back(iwf[k].real());
    result.s = circular_shift(result.s, shift);
    result.s = normalize<double>(result.s);
    return result;
  } catch (...) {
    throw;
  };
}

References mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::seismic::CoreTimeSeries::s, mspass::seismic::CoreTimeSeries::set_dt(), mspass::seismic::BasicTimeSeries::set_live(), mspass::seismic::CoreTimeSeries::set_npts(), mspass::seismic::CoreTimeSeries::set_t0(), mspass::seismic::BasicTimeSeries::set_tref(), and mspass::algorithms::deconvolution::ComplexArray::size().

operator=()

ShapingWavelet & mspass::algorithms::deconvolution::ShapingWavelet::operator=

(

const ShapingWavelet &

parent

)

Assignment operator.

                                                                      {
  if (this != &parent) {
    nfft = parent.nfft;
    w = parent.w;
    dt = parent.dt;
    df = parent.df;
    wavelet_name = parent.wavelet_name;
  }
  return *this;
}

sample_interval()

double mspass::algorithms::deconvolution::ShapingWavelet::sample_interval ( )

inline

Return the sample interval.

{ return dt; };

size()

int mspass::algorithms::deconvolution::ShapingWavelet::size ( ) const

inline

Return the number of complex frequency samples.

{ return w.size(); };

References mspass::algorithms::deconvolution::ComplexArray::size().

swap()

void mspass::algorithms::deconvolution::ShapingWavelet::swap ( ShapingWavelet &  other )

inlinenoexcept

Exchange complete shaping state without allocation.

                                            {
    std::swap(nfft, other.nfft);
    w.swap(other.w);
    std::swap(dt, other.dt);
    std::swap(df, other.df);
    wavelet_name.swap(other.wavelet_name);
  }

References mspass::algorithms::deconvolution::ComplexArray::swap().

type()

std::string mspass::algorithms::deconvolution::ShapingWavelet::type ( )

inline

Return the shaping wavelet type name.

{ return wavelet_name; };

wavelet()

ComplexArray * mspass::algorithms::deconvolution::ShapingWavelet::wavelet ( )

inline

Return a pointer to the shaping wavelet this object defines in the frequency domain.

{ return &w; };

---

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

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