Implements the Generalized Iterative Method of Wang and Pavlis. More...
#include <GeneralIterDecon.h>
Inheritance diagram for mspass::algorithms::deconvolution::GeneralIterDecon:
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.
53 : ScalarDecon() {
54 const string base_error("GeneralIterDecon AntelopePf contructor: ");
55 stringstream ss; // used for constructing error messages
56 /* The pf used for initializing this object has Antelope Arr section
57 for each algorithm. Since the generalized iterative method is a
58 two-stage algorithm we have a section for the iterative algorithm
59 and a (variable) sectin for the preprocessor algorithm. We use
60 the AntelopePf to parse this instead of raw antelope pfget
61 C calls. */
62 try {
65 IterDeconType dct = parse_for_itertype(mdgiter);
66 this->decon_type = dct;
67 double ts, te;
68 ts = mdgiter.get<double>("full_data_window_start");
69 te = mdgiter.get<double>("full_data_window_end");
70 dwin = TimeWindow(ts, te);
71 ts = mdgiter.get<double>("deconvolution_data_window_start");
72 te = mdgiter.get<double>("deconvolution_data_window_end");
73 fftwin = TimeWindow(ts, te);
74 ts = mdgiter.get<double>("noise_window_start");
75 te = mdgiter.get<double>("noise_window_end");
76 nwin = TimeWindow(ts, te);
77 /* We need to make sure the noise and decon windows are inside the
78 * full_data_window*/
79 if (fftwin.start < dwin.start || fftwin.end > dwin.end) {
80 stringstream ss;
81 ss << base_error << "decon window error" << endl
82 << "Wavelet inversion window is not inside analysis window" << endl
83 << "full_data_window (analysis) range=" << dwin.start << " to "
84 << dwin.end << endl
85 << "decon_window (wavelet inversion) range=" << fftwin.start << " to "
86 << fftwin.end << endl;
87 throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
88 }
89 noise_component = mdgiter.get<int>("noise_component");
90 double target_dt = mdgiter.get<double>("target_sample_interval");
91 int maxns = static_cast<int>((fftwin.end - fftwin.start) / target_dt);
92 ++maxns; // Add one - points not intervals
93 nfft = nextPowerOf2(maxns);
94 /* This should override this even if it was previously set */
95 mdgiter.put("operator_nfft", nfft);
96 this->ScalarDecon::changeparameter(mdgiter);
97 shapingwavelet = ShapingWavelet(mdgiter, nfft);
98 /* Note minus sign here */
99 time_shift = -(dwin.start) / target_dt;
100 AntelopePf mdleaf;
101 /* We make sure the window parameters in each algorithm mactch what
102 is set for this algorithm. Abort if they are not consistent. The
103 test code is a bit repetitious but a necessary evil to allow the message
104 to be clearer. */
105 int n1, n2; // temporaries used below - needed because declrations illegal
106 // inside case
107 switch (decon_type) {
108 case WATER_LEVEL:
109 mdleaf = md.get_branch("water_level");
110 ts = mdleaf.get<double>("deconvolution_data_window_start");
111 te = mdleaf.get<double>("deconvolution_data_window_end");
112 if ((ts != fftwin.start) || (te != fftwin.end)) {
113 ss << base_error
114 << "water level method specification of processing window is not"
115 " consistent with gid parameters"
116 << endl
117 << "water level parameters: deconvolution_data_window_start=" << ts
118 << ", decon_window_end=" << te << endl
119 << "GID parameters: decon_window_start=" << fftwin.start
120 << ", decon_window_end=" << fftwin.end << endl;
121 throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
122 }
123 preprocessor = new WaterLevelDecon(mdleaf);
124 break;
125 case LEAST_SQ:
126 mdleaf = md.get_branch("least_square");
127 ts = mdleaf.get<double>("deconvolution_data_window_start");
128 te = mdleaf.get<double>("deconvolution_data_window_end");
129 if ((ts != fftwin.start) || (te != fftwin.end)) {
130 ss << base_error
131 << "least square method specification of processing window is not"
132 " consistent with gid parameters"
133 << endl
134 << "least square parameters: deconvolution_data_window_start=" << ts
135 << ", decon_window_end=" << te << endl
136 << "GID parameters: decon_window_start=" << fftwin.start
137 << ", decon_window_end=" << fftwin.end << endl;
138 throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
139 }
140 preprocessor = new LeastSquareDecon(mdleaf);
141 break;
142 case MULTI_TAPER:
143 default:
144 mdleaf = md.get_branch("multi_taper");
145 ts = mdleaf.get<double>("deconvolution_data_window_start");
146 te = mdleaf.get<double>("deconvolution_data_window_end");
147 if ((ts != fftwin.start) || (te != fftwin.end)) {
148 ss << base_error
149 << "mulittaper method specification of processing window is not"
150 " consistent with gid parameters"
151 << endl
152 << "multitaper parameters: deconvolution_data_window_start=" << ts
153 << ", decon_window_end=" << te << endl
154 << "GID parameters: decon_window_start=" << fftwin.start
155 << ", decon_window_end=" << fftwin.end << endl;
156 throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
157 }
158 /* Here we also have to test the noise parameters, but the gid
159 window can be different fromt that passed to the mulittaper method.
160 Hence we test only that the mulittaper noise window is within the bounds
161 of the gid noise window */
162 n1 = static_cast<int>((fftwin.end - fftwin.start) / target_dt) + 1;
163 n2 = static_cast<int>((nwin.end - nwin.start) / target_dt) + 1;
164 if (n1 > n2) {
165 ss << base_error << "inconsistent noise window specification" << endl
166 << "multitaper parameters specify taper length=" << n1 << " samples"
167 << endl
168 << "GID noise window parameters define noise_window_start="
169 << nwin.start << " and noise_window_end=" << nwin.end << endl
170 << "The GID window has a length of " << n2 << " samples" << endl
171 << "GID implementation insists multitaper noise window be smaller "
172 "or equal to GID noise window"
173 << endl;
174 throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
175 }
176 preprocessor = new MultiTaperXcorDecon(mdleaf);
177 };
178 /* Because this may evolve we make this a private method to
179 make changes easier to implement. */
180 this->construct_weight_penalty_function(mdgiter);
181 /* Set convergencd parameters from md keys */
182 iter_max = mdgiter.get<int>("maximum_iterations");
183 lw_linf_floor = mdgiter.get<double>("lag_weight_Linf_floor");
184 lw_l2_floor = mdgiter.get<double>("lag_weight_rms_floor");
185 resid_linf_prob =
186 mdgiter.get<double>("residual_noise_rms_probability_floor");
187 resid_l2_tol = mdgiter.get<double>("residual_fractional_improvement_floor");
188 } catch (...) {
189 throw;
190 };
191}
AntelopePf get_branch(const std::string key) const
used for subtrees (nested Arrs in antelope)
Definition AntelopePf.cc:392
References mspass::utility::Metadata::get(), and mspass::utility::AntelopePf::get_branch().
◆ ~GeneralIterDecon()
| mspass::algorithms::deconvolution::GeneralIterDecon::~GeneralIterDecon | ( | ) |
192{ delete preprocessor; }
Member Function Documentation
◆ actual_output()
|
inlinevirtual |
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.
98 {
100 };
virtual mspass::seismic::CoreTimeSeries actual_output()=0
Return the actual output of the deconvolution operator.
References mspass::algorithms::deconvolution::ScalarDecon::actual_output().
◆ changeparameter()
|
inlinevirtual |
Implements mspass::algorithms::deconvolution::BasicDeconOperator.
77 {
78 this->preprocessor->changeparameter(md);
79 };
◆ getresult()
| CoreSeismogram mspass::algorithms::deconvolution::GeneralIterDecon::getresult | ( | ) |
698 {
699 try {
700 string base_error("GeneralIterDecon::getresult: ");
701 CoreSeismogram result(d_all);
702 /* We will make the output the size of the processing window for the
703 iteration. May want to alter this to trim the large lag that would not
704 be allowed due to wavelet duration anyway, BUT for GID method the
705 wavelet should be compact enough that should be a small factor. Hence
706 for now I omit that complexity until proven to be an issue. */
707 result = WindowData(result, dwin);
708 result.u.zero();
709 /* The spike sequences uses the time reference of the data in the
710 private copy r. This is the computed offset in samples to correct
711 lags in the spikes list container to be at correct time in result */
712 double dt0;
713 int delta_col;
714 dt0 = result.t0() - r.t0();
715 delta_col = round(dt0 / r.dt());
716 list<ThreeCSpike>::iterator sptr;
717 int k, resultcol;
718 for (sptr = spikes.begin(); sptr != spikes.end(); ++sptr) {
719 if (((sptr->col) < 0) || ((sptr->col) >= result.npts()))
720 throw MsPASSError(
721 base_error +
722 "Coding error - spike lag is outside output data range",
723 ErrorSeverity::Fatal);
724 resultcol = (sptr->col) - delta_col;
725 for (k = 0; k < 3; ++k)
726 result.u(k, resultcol) = sptr->u[k];
727 }
728 string wtype = this->shapingwavelet.type();
729 if (wtype != "none") {
730 CoreTimeSeries w(this->shapingwavelet.impulse_response());
731 result = sparse_convolve(w, result);
732 }
733 return result;
734 } catch (...) {
735 throw;
736 };
737}
mspass::seismic::CoreTimeSeries impulse_response()
Definition ShapingWavelet.cc:265
◆ ideal_output()
|
inline |
95 {
96 return this->preprocessor->ideal_output();
97 };
◆ inverse_wavelet() [1/2]
|
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.
101 {
103 };
virtual mspass::seismic::CoreTimeSeries inverse_wavelet()=0
Return a FIR represention of the inverse filter.
References mspass::algorithms::deconvolution::ScalarDecon::inverse_wavelet().
◆ inverse_wavelet() [2/2]
|
inlinevirtual |
Implements mspass::algorithms::deconvolution::ScalarDecon.
◆ load() [1/2]
| int mspass::algorithms::deconvolution::GeneralIterDecon::load | ( | const mspass::seismic::CoreSeismogram & | d, |
| mspass::algorithms::TimeWindow | dwin | ||
| ) |
289 {
290 try {
291 dwin = dwin_in;
292 /* First we load the requested window. Note we MUST always make this window
293 a bit larger than the range of desired lags as the iterative algorithm will
294 not allow lags at the edges (defined by a construction parameter
295 wavelet_pad)
296 */
297 d_all = WindowData(draw, dwin);
298 ndwin = d_all.npts();
299 return 0;
300 } catch (...) {
301 throw;
302 };
303}
◆ load() [2/2]
| 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.
320 {
321 try {
322 int iretn, iret;
323 iretn = this->loadnoise(draw, nwin);
324 iret = this->load(draw, dwin);
325 return (iretn + iret);
326 } catch (...) {
327 throw;
328 };
329}
◆ loadnoise()
| int mspass::algorithms::deconvolution::GeneralIterDecon::loadnoise | ( | const mspass::seismic::CoreSeismogram & | d, |
| mspass::algorithms::TimeWindow | nwin | ||
| ) |
305 {
306 try {
307 nwin = nwin_in;
308 n = WindowData(draw, nwin);
309 nnwin = n.npts();
310 double ret = this->compute_resid_linf_floor();
311 if (ret > 0)
312 return 0;
313 else
314 return 1;
315 } catch (...) {
316 throw;
317 };
318}
◆ process()
|
virtual |
Implements mspass::algorithms::deconvolution::ScalarDecon.
480 {
481 const string base_error("GeneralIterDecon::process method: ");
482 try {
483 /* We first have to run the signal processing style deconvolution.
484 This is defined by the base pointer available through the symbol
485 preprocessor. All those altorithms require load methods to be called
486 to initate the computation. A complication is that the multitaper is
487 different and requires a noise signal to also be loaded through loadnoise.
488 That complicates this a bit below, but the flow of the algorithm should
489 still be clear. Outer loop is over the three components were we assemble
490 a full 3c record. Note this is the same algorithm use in trace_decon
491 for anything but this iterative algorithm.
492 */
493 /* d_decon will hold the preprocessor output. We normally expect to
494 derive it by windowing of t_all. We assume WindowData will be
495 successful - constructor should guarantee that. */
496 d_decon = WindowData(d_all, fftwin);
497 dmatrix uwork(d_decon.u);
498 uwork.zero();
499 /* We assume loadnoise has been called previously to set put the
500 right data here. We need a scalar function to pass to the multtitaper
501 algorithm though. */
502 if (decon_type == MULTI_TAPER) {
503 CoreTimeSeries nts(ExtractComponent(n, noise_component));
504 dynamic_cast<MultiTaperXcorDecon *>(preprocessor)->loadnoise(nts.s);
505 }
506 /* For this case of receiver function deconvolution we always get the
507 wavelet from component 2 - assumed here to be Z or L. */
508 CoreTimeSeries srcwavelet(ExtractComponent(d_decon, 2));
509 for (int k = 0; k < 3; ++k) {
510 CoreTimeSeries dcomp(ExtractComponent(d_decon, k));
511 /* Need the qualifier or we get the wrong overloaded
512 * load method */
513 preprocessor->ScalarDecon::load(srcwavelet.s, dcomp.s);
514 preprocessor->process();
515 vector<double> deconout(preprocessor->getresult());
516 int copysize = deconout.size();
518 copysize = d_decon.npts();
519 cblas_dcopy(copysize, &(deconout[0]), 1, uwork.get_address(k, 0), 3);
520 }
521 d_decon.u = uwork;
522 /* The inverse wavelet and the actual output signals are determined in all
523 current algorithms from srcwavelet. Hence, what is now stored will work.
524 If this is extended make sure that condition is satisfied.
525
526 We extract the inverse filter and use a time domain convolution with
527 data in the longer time window. Note for efficiency may want to
528 convert this to a frequency domain convolution if it proves to be
529 a bottleneck */
530 // double dt;
531 // dt=this->shapingwavelet.sample_interval();
532 // TimeSeries winv=this->preprocessor->inverse_wavelet(tshift,d_decon.t0());
533 // DEBUG
534 // cerr << "inverse wavelet tshift="<<tshift/5.0<<"
535 // d_decon.to="<<d_decon.t0()<<endl; TimeSeries
536 // winv=this->preprocessor->inverse_wavelet(tshift/5.0,d_decon.t0());
537 // TimeSeries winv=this->preprocessor->inverse_wavelet(0.0,d_decon.t0());
539 // This is a test - need a more elegant solution if it works. Remove me
540 // when finished with this test
541 // if(d_decon.t0!=0) winv.t0 -= d_decon.t0;
542 // DEBUG
543 /*
544 cerr << "Inverse wavelet"<<endl
545 << winv
546 <<endl;
547*/
548 /* The actual output signal is used in the iterative
549 * recursion of this algorithm. For efficiency it is important
550 * to trim the fir filter. The call to trim does that.*/
551 CoreTimeSeries actual_out(this->preprocessor->actual_output());
552 // DEBUG
553 /*
554 cerr << "Actual output raw"<<endl;
555 cerr << actual_out<<endl;
556*/
557 actual_out = trim(actual_out);
558 actual_o_fir = actual_out.s;
559 actual_o_0 = actual_out.sample_number(0.0);
560 double peak_scale = actual_o_fir[actual_o_0];
561 vector<double>::iterator aoptr;
562 for (aoptr = actual_o_fir.begin(); aoptr != actual_o_fir.end(); ++aoptr)
563 (*aoptr) /= peak_scale;
564 /* This is the size of the inverse wavelet convolution transient
565 we use it to prevent iterations in transient region of the deconvolved
566 data */
567 wavelet_pad = winv.s.size();
568 if (2 * wavelet_pad > ndwin) {
569 stringstream ss;
570 ss << base_error << "Inadequate data window size" << endl
571 << "trimmed FIR filter size for actual output signal=" << wavelet_pad
572 << endl
573 << "Data window length=" << ndwin << endl
574 << "Window size must be larger than two times FIR filter size" << endl;
575 throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
576 }
577 /* These two signals should be trimmed by winv.npts() on both ends to remove
578 sections that are pure edge transients.
579 REMOVE me when that is done*/
580
581 r = sparse_convolve(winv, d_all);
582 /* Replace n by convolution with inverse wavelet to get the levels correct
583 */
584 n = sparse_convolve(winv, n);
585 TimeWindow trimwin;
588 n = WindowData(n, trimwin);
589 // double nfloor;
590 // nfloor=compute_resid_linf_floor();
591 // DEBUG - for debug always print this. Should be a verbose option
592 // cerr << "Computed noise floor="<<nfloor<<endl;
593
594 /* d_all now contains the deconvolved data. Now enter the
595 generalized iterative method recursion */
596 int i, k;
597 lag_weights.clear();
598 vector<double> amps, wamps; // raw and weighted amplitudes
599 amps.reserve(r.npts());
600 wamps.reserve(r.npts());
601 /* We need these iterators repeatedly in the main loop below */
602 vector<double>::iterator amax;
604 lag_weights.push_back(1.0);
605 // DEBUG - temporarily disabled for testing
606 // for(i=0; i<wavelet_pad; ++i) lag_weights[i]=0.0;
607 // for(i=0; i<wavelet_pad; ++i) lag_weights[r.npts()-i-1]=0.0;
608 /* These are initial values of convergence parameters */
609 lw_linf_initial = 1.0;
610 lw_l2_initial = 1.0;
611 resid_linf_initial = Linf(r.u);
612 resid_l2_initial = L2(r.u);
613 iter_count = 0;
614 // DEBUG - remove after testing
615 lw_linf_history.push_back(lw_linf_initial);
616 lw_l2_history.push_back(lw_l2_initial);
617 resid_l2_history.push_back(resid_l2_initial);
618 resid_linf_history.push_back(resid_linf_initial);
619 do {
620 /* Compute the vector of amplitudes and find the maximum */
621 amps = amp3c(r.u);
622 wamps.clear();
623 for (k = 0; k < amps.size(); ++k)
624 wamps.push_back(amps[k] * lag_weights[k]);
625 amax = max_element(wamps.begin(), wamps.end());
626 /* The generic distance algorithm used here returns an integer
627 that would work to access amps[imax] so we can use the same index
628 for the column of the data in d.u. */
629 int imax = distance(wamps.begin(), amax);
630 /* Save the 3c amplitude at this lag to the spike condensed
631 respresentation of the output*/
632 ThreeCSpike spk(r.u, imax);
633 spikes.push_back(spk);
634 // DEBUG
636 << " amps=" << spk.u[0] << ", " << spk.u[1] << ", " << spk.u[2]
637 << endl;
638 /* This private method defines how the lag_weights vector is changed
639 in the vicinity of this spike. The tacit assumption is the weight is
640 made smaller (maybe even zero) at the spike point and a chosen recipe
641 for points in the vicinity of that spike */
642 this->update_lag_weights(imax);
643 /* Subtract the actual output from the data at lag imax. Note
644 We don't test validity of the lag in spk, but depend on dmatrix to throw
645 an exception of the range is invalid */
646 this->update_residual_matrix(spk);
647 ++iter_count;
648 } while (this->has_not_converged());
649 if (iter_count >= iter_max)
650 throw MsPASSError("GeneralIterDecon::process did not converge",
651 ErrorSeverity::Suspect);
652 } catch (...) {
653 throw;
654 };
655}
double endtime() const noexcept
Definition CoreSeismogram.h:497
◆ 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.
738 {
739 cerr << "QCMetrics method not yet implemented" << endl;
740 exit(-1);
741}
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
