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. More...
	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 ()

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/3]

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() [2/3]

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() [3/3]

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

Member Function Documentation

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);
        /* old API
        result.tref=TimeReferenceType::Relative;
        result.ns=nfft;
        result.dt=dt;
        result.t0 = dt*(-(double)nfft/2);
        result.live=true;
        */
 
        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);
        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().

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