mspass::algorithms::deconvolution::GeneralIterDecon Class Reference

Implements the Generalized Iterative Method of Wang and Pavlis. More...

#include <GeneralIterDecon.h>

Inheritance diagram for mspass::algorithms::deconvolution::GeneralIterDecon:

Public Member Functions
	GeneralIterDecon (mspass::utility::AntelopePf &md)
	Create an initialize an operator for subseqquent processing. More...
	GeneralIterDecon (const GeneralIterDecon &parent)
void	changeparameter (const mspass::utility::Metadata &md)
int	load (const mspass::seismic::CoreSeismogram &d, mspass::algorithms::TimeWindow dwin)
int	loadnoise (const mspass::seismic::CoreSeismogram &d, mspass::algorithms::TimeWindow nwin)
int	load (const mspass::seismic::CoreSeismogram &d, mspass::algorithms::TimeWindow dwin, mspass::algorithms::TimeWindow nwin)
	Load all needed data and process. More...
void	process ()
mspass::seismic::CoreSeismogram	getresult ()
mspass::seismic::CoreTimeSeries	ideal_output ()
mspass::seismic::CoreTimeSeries	actual_output ()
mspass::seismic::CoreTimeSeries	inverse_wavelet ()
	Return a FIR represention of the inverse filter. More...
mspass::seismic::CoreTimeSeries	inverse_wavelet (double t0parent)
mspass::utility::Metadata	QCMetrics ()
	Return appropriate quality measures. More...
Public Member Functions inherited from mspass::algorithms::deconvolution::ScalarDecon
	ScalarDecon (const mspass::utility::Metadata &md)
	ScalarDecon (const std::vector< double > &d, const std::vector< double > &w)
	ScalarDecon (const ScalarDecon &parent)
int	load (const std::vector< double > &wavelet, const std::vector< double > &data)
	Load all data required for decon. More...
int	loaddata (const std::vector< double > &data)
int	loadwavelet (const std::vector< double > &wavelet)
ScalarDecon &	operator= (const ScalarDecon &parent)
std::vector< double >	getresult ()
void	changeparameter (const mspass::utility::Metadata &md)
void	change_shaping_wavelet (const ShapingWavelet &nsw)
mspass::seismic::CoreTimeSeries	ideal_output ()

Additional Inherited Members
Protected Attributes inherited from mspass::algorithms::deconvolution::ScalarDecon
std::vector< double >	data
std::vector< double >	wavelet
std::vector< double >	result
ShapingWavelet	shapingwavelet

Detailed Description

Implements the Generalized Iterative Method of Wang and Pavlis.

This class is an extension of the idea o the generalized iterative method described in the paper by Wang and Pavlis. Their original paper worked on scalar seismograms. This implemenation differs by iterating on the vector data effectively reducing a 3c seismogram to a sequence of 3D vectors at computed lags. The iterative reduction then iterates on the data matrix. In addition, Wang's original implementation did the iteration by repeated forward and inverse fourier transforms. This algorithm works completely in the time domain.

The class retains a concept Wang used in the original implementation. Many OOP books stress the idea of construction is initialization. That is not pruddent fo this case since construction involves some work that we do not want to repeat for each new seismogram. Hence, construction only initializes parameters required to run the algorithm. Calculation is triggered by calling one of the two overloaded load methods. The user should be warned, however, that we did not build ina test to assure the noise and signal windows are consistent. Anyone wanting to use this code outside of applications developed by the authors must recognize That the assure consistency they need to call loadnoise before load OR call only the load method that includes a signal and noise window.

Constructor & Destructor Documentation

GeneralIterDecon()

mspass::algorithms::deconvolution::GeneralIterDecon::GeneralIterDecon

(

mspass::utility::AntelopePf &

md

)

Create an initialize an operator for subseqquent processing.

This is the primary constructor for this object. The PfStyleMetadata object is expected to contain a suite of parameters used to define how this operator will behave along with what inverse filter is to be constructed ans applied before the iterative algorithm is applied. The pf is required to have multiple Arr sections the components.

                                                         : ScalarDecon()
{
    const string base_error("GeneralIterDecon AntelopePf contructor:  ");
    stringstream ss;   // used for constructing error messages
    /* The pf used for initializing this object has Antelope Arr section
    for each algorithm.   Since the generalized iterative method is a
    two-stage algorithm we have a section for the iterative algorithm
    and a (variable) sectin for the preprocessor algorithm.  We use
    the AntelopePf to parse this instead of raw antelope pfget
    C calls. */
    try {
        AntelopePf md=mdtoplevel.get_branch("deconvolution_operator_type");
        AntelopePf mdgiter=md.get_branch("generalized_iterative_deconvolution");
        IterDeconType dct = parse_for_itertype(mdgiter);
        this->decon_type=dct;
        double ts,te;
        ts=mdgiter.get<double>("full_data_window_start");
        te=mdgiter.get<double>("full_data_window_end");
        dwin=TimeWindow(ts,te);
        ts=mdgiter.get<double>("deconvolution_data_window_start");
        te=mdgiter.get<double>("deconvolution_data_window_end");
        fftwin=TimeWindow(ts,te);
        ts=mdgiter.get<double>("noise_window_start");
        te=mdgiter.get<double>("noise_window_end");
        nwin=TimeWindow(ts,te);
        /* We need to make sure the noise and decon windows are inside the full_data_window*/
        if(fftwin.start<dwin.start || fftwin.end>dwin.end)
        {
            stringstream ss;
            ss << base_error << "decon window error"<<endl
               << "Wavelet inversion window is not inside analysis window"<<endl
               << "full_data_window (analysis) range="<<dwin.start<<" to "<<dwin.end<<endl
               << "decon_window (wavelet inversion) range="<<fftwin.start<<" to "<<fftwin.end<<endl;
            throw MsPASSError(ss.str(),ErrorSeverity::Invalid);
        }
        noise_component=mdgiter.get<int>("noise_component");
        double target_dt=mdgiter.get<double>("target_sample_interval");
        int maxns=static_cast<int>((fftwin.end-fftwin.start)/target_dt);
        ++maxns;   // Add one - points not intervals
        nfft=nextPowerOf2(maxns);
        /* This should override this even if it was previously set */
        mdgiter.put("operator_nfft",nfft);
        this->ScalarDecon::changeparameter(mdgiter);
        shapingwavelet=ShapingWavelet(mdgiter,nfft);
        /* Note minus sign here */
        time_shift = -(dwin.start)/target_dt;
        AntelopePf mdleaf;
        /* We make sure the window parameters in each algorithm mactch what
        is set for this algorithm.  Abort if they are not consistent.  The
        test code is a bit repetitious but a necessary evil to allow the message
        to be clearer. */
        int n1,n2; //temporaries used below -  needed because declrations illegal inside case
        switch(decon_type)
        {
        case WATER_LEVEL:
            mdleaf=md.get_branch("water_level");
            ts=mdleaf.get<double>("deconvolution_data_window_start");
            te=mdleaf.get<double>("deconvolution_data_window_end");
            if((ts!=fftwin.start) || (te!=fftwin.end))
            {
                ss << base_error
                   <<"water level method specification of processing window is not"
                   " consistent with gid parameters"<<endl
                   << "water level parameters:  deconvolution_data_window_start="<<ts
                   << ", decon_window_end="<<te<<endl
                   << "GID parameters: decon_window_start="<<fftwin.start
                   << ", decon_window_end="<<fftwin.end<<endl;
                throw MsPASSError(ss.str(),ErrorSeverity::Invalid);
            }
            preprocessor=new WaterLevelDecon(mdleaf);
            break;
        case LEAST_SQ:
            mdleaf=md.get_branch("least_square");
            ts=mdleaf.get<double>("deconvolution_data_window_start");
            te=mdleaf.get<double>("deconvolution_data_window_end");
            if((ts!=fftwin.start) || (te!=fftwin.end))
            {
                ss << base_error
                   <<"least square method specification of processing window is not"
                   " consistent with gid parameters"<<endl
                   << "least square parameters:  deconvolution_data_window_start="<<ts
                   << ", decon_window_end="<<te<<endl
                   << "GID parameters: decon_window_start="<<fftwin.start
                   << ", decon_window_end="<<fftwin.end<<endl;
                throw MsPASSError(ss.str(),ErrorSeverity::Invalid);
            }
            preprocessor=new LeastSquareDecon(mdleaf);
            break;
        case MULTI_TAPER:
        default:
            mdleaf=md.get_branch("multi_taper");
            ts=mdleaf.get<double>("deconvolution_data_window_start");
            te=mdleaf.get<double>("deconvolution_data_window_end");
            if((ts!=fftwin.start) || (te!=fftwin.end))
            {
                ss << base_error
                   <<"mulittaper method specification of processing window is not"
                   " consistent with gid parameters"<<endl
                   << "multitaper parameters:  deconvolution_data_window_start="<<ts
                   << ", decon_window_end="<<te<<endl
                   << "GID parameters: decon_window_start="<<fftwin.start
                   << ", decon_window_end="<<fftwin.end<<endl;
                throw MsPASSError(ss.str(),ErrorSeverity::Invalid);
            }
            /* Here we also have to test the noise parameters, but the gid
            window can be different fromt that passed to the mulittaper method.
            Hence we test only that the mulittaper noise window is within the bounds
            of the gid noise window */
            n1=static_cast<int>((fftwin.end-fftwin.start)/target_dt)+1;
            n2=static_cast<int>((nwin.end-nwin.start)/target_dt)+1;
            if(n1>n2)
            {
                ss << base_error << "inconsistent noise window specification"<<endl
                   << "multitaper parameters specify taper length="<<n1<<" samples"<<endl
                   << "GID noise window parameters define noise_window_start="<<nwin.start
                   << " and noise_window_end="<<nwin.end<<endl
                   << "The GID window has a length of "<<n2<<" samples"<<endl
                   << "GID implementation insists multitaper noise window be smaller or equal to GID noise window"
                   <<endl;
                throw MsPASSError(ss.str(),ErrorSeverity::Invalid);
            }
            preprocessor=new MultiTaperXcorDecon(mdleaf);
        };
        /* Because this may evolve we make this a private method to
        make changes easier to implement. */
        this->construct_weight_penalty_function(mdgiter);
        /* Set convergencd parameters from md keys */
        iter_max=mdgiter.get<int>("maximum_iterations");
        lw_linf_floor=mdgiter.get<double>("lag_weight_Linf_floor");
        lw_l2_floor=mdgiter.get<double>("lag_weight_rms_floor");
        resid_linf_prob=mdgiter.get<double>("residual_noise_rms_probability_floor");
        resid_l2_tol=mdgiter.get<double>("residual_fractional_improvement_floor");
    } catch(...) {
        throw;
    };
}

References mspass::utility::Metadata::get(), and mspass::utility::AntelopePf::get_branch().

Member Function Documentation

actual_output()

mspass::seismic::CoreTimeSeries mspass::algorithms::deconvolution::GeneralIterDecon::actual_output ( )

inlinevirtual

\brif 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.

                                                {
        return this->preprocessor->actual_output();
    };

References mspass::algorithms::deconvolution::ScalarDecon::actual_output().

inverse_wavelet()

mspass::seismic::CoreTimeSeries mspass::algorithms::deconvolution::GeneralIterDecon::inverse_wavelet ( )

inlinevirtual

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->preprocessor->inverse_wavelet();
    };

References mspass::algorithms::deconvolution::ScalarDecon::inverse_wavelet().

load()

int mspass::algorithms::deconvolution::GeneralIterDecon::load	(	const mspass::seismic::CoreSeismogram &	d,
		mspass::algorithms::TimeWindow	dwin,
		mspass::algorithms::TimeWindow	nwin
	)

Load all needed data and process.

This method is little more than a call to loadnoise followed immediately by a call to load that is assumed to initiate the computation.

{
    try {
        int iretn, iret;
        iretn=this->loadnoise(draw,nwin);
        iret=this->load(draw,dwin);
        return (iretn+iret);
    } catch(...) {
        throw;
    };
}

QCMetrics()

Metadata mspass::algorithms::deconvolution::GeneralIterDecon::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.

{
    cerr << "QCMetrics method not yet implemented"<<endl;
    exit(-1);
}

---

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

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