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)
	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)
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

                   {
    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 = md.get_int("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);
      }
    }
    /* 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. */
    gsl_fft_complex_wavetable *wavetable;
    gsl_fft_complex_workspace *workspace;
    wavetable = gsl_fft_complex_wavetable_alloc(nfft);
    workspace = gsl_fft_complex_workspace_alloc(nfft);
    string wavelettype = md.get_string("shaping_wavelet_type");
    wavelet_name = wavelettype;
    dt = md.get_double("shaping_wavelet_dt");
    double *r;
    if (wavelettype == "gaussian") {
      float fpeak = md.get_double("shaping_wavelet_frequency");
      // construct wavelet and fft
      r = gaussian(fpeak, (float)dt, nfft);
      w = ComplexArray(nfft, r);
      gsl_fft_complex_forward(w.ptr(), 1, nfft, wavetable, workspace);
      delete[] r;
    }
    /* 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 = (float)md.get_double("shaping_wavelet_frequency");
      // construct wavelet and fft
      r = rickerwavelet(fpeak, (float)dt, nfft);
      // DEBUG
      // cerr << "Ricker shaping wavelet"<<endl;
      // for(int k=0;k<nfft;++k) cerr << r[k]<<endl;
      w = ComplexArray(nfft, r);
      gsl_fft_complex_forward(w.ptr(), 1, nfft, wavetable, workspace);
      delete[] r;
    } else if (wavelettype == "butterworth") {
      double f3db_lo, f3db_hi;
      f3db_lo = md.get_double("f3db_lo");
      f3db_hi = md.get_double("f3db_hi");
      int npoles_lo, npoles_hi;
      npoles_lo = md.get_int("npoles_lo");
      npoles_hi = md.get_int("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 = md.get_double("time_bandwidth_product");
      double target_pulse_width = md.get_double("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)));
      double *wtmp = slepian0(tbp, nwsize);
      double *work = new double[nfft];
      for (int k = 0; k < nfft; ++k)
        work[k] = 0.0;
      dcopy(nwsize, wtmp, 1, work, 1);
      delete[] wtmp;
      w = ComplexArray(nfft, work);
      gsl_fft_complex_forward(w.ptr(), 1, nfft, wavetable, workspace);
      delete[] work;
    } else if (wavelettype == "none") {
      /* Prototype code issued an error in this condition, but we accept it
      here as an option defined by none.  We could do this by putting
            all ones in the w array but using a delta function at zero lage
            avoids scaling issues for little cost - this assumes this
            object is created only occassionally and not millions of times */
      double *r = new double[nfft];
      for (int k = 0; k < nfft; ++k)
        r[k] = 0.0;
      r[0] = 1.0;
      w = ComplexArray(nfft, r);
      delete[] r;
    } else {
      throw MsPASSError(
          base_error + "illegal value for shaping_wavelet_type=" + wavelettype,
          ErrorSeverity::Invalid);
    }
    gsl_fft_complex_wavetable_free(wavetable);
    gsl_fft_complex_workspace_free(workspace);
    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_double(), mspass::utility::Metadata::get_int(), mspass::utility::Metadata::get_string(), mspass::utility::MsPASSError::log_error(), and mspass::algorithms::Butterworth::transfer_function().

ShapingWavelet() [3/6]

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

                                                         {
  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();
  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];
  }
  gsl_fft_complex_wavetable *wavetable;
  gsl_fft_complex_workspace *workspace;
  wavetable = gsl_fft_complex_wavetable_alloc(nfft);
  workspace = gsl_fft_complex_workspace_alloc(nfft);
  w = ComplexArray(nfft, &(dwork[0]));
  gsl_fft_complex_forward(w.ptr(), 1, nfft, wavetable, workspace);
  gsl_fft_complex_wavetable_free(wavetable);
  gsl_fft_complex_workspace_free(workspace);
}

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") {
  nfft = n;
  dt = dtin;
  df = 1.0 / (dt * static_cast<double>(n));
  double *r;
  r = rickerwavelet((float)fpeak, (float)dt, nfft);
  w = ComplexArray(nfft, r);
  gsl_fft_complex_wavetable *wavetable;
  gsl_fft_complex_workspace *workspace;
  wavetable = gsl_fft_complex_wavetable_alloc(nfft);
  workspace = gsl_fft_complex_workspace_alloc(nfft);
  gsl_fft_complex_forward(w.ptr(), 1, nfft, wavetable, workspace);
  gsl_fft_complex_wavetable_free(wavetable);
  gsl_fft_complex_workspace_free(workspace);
  delete[] r;
}

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
nplolehi	is the number of poles for the upper corner filter
f3dbhi	is the 3db point for the high corner of the passband.

    : wavelet_name("butterworth") {
  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

)

                                                           : 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 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();
    gsl_fft_complex_wavetable *wavetable;
    gsl_fft_complex_workspace *workspace;
    wavetable = gsl_fft_complex_wavetable_alloc(nfft);
    workspace = gsl_fft_complex_workspace_alloc(nfft);
    /* 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, wavetable, workspace);
    gsl_fft_complex_wavetable_free(wavetable);
    gsl_fft_complex_workspace_free(workspace);
    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::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().

operator=()

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

(

const ShapingWavelet &

parent

)

                                                                      {
  if (this != &parent) {
    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 dt; };

size()

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

inline

{ return w.size(); };

type()

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

inline

{ 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