C++ DoubleArray::slice方法代码示例

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor::computeNoise(const DoubleArray &data, int i1, int i2, double *offset, double *amplitude) {
	// compute offset and rms within the time window
	if(i1<0) i1=0;
	if(i2<0) return false;
	if(i2>(int)data.size()) i2=(int)data.size();

	// If noise window is zero return an amplitude and offset of zero as well.
	if ( i2-i1 == 0 ) {
		*amplitude = 0;
		*offset = 0;
		return true;
	}

	DoubleArrayPtr d = static_cast<DoubleArray*>(data.slice(i1, i2));

	double ofs, amp;

	// compute pre-arrival offset
	ofs = d->median();
	// compute rms after removing offset
	amp = 2 * d->rms(ofs);

	if ( offset ) *offset = ofs;
	if ( amplitude ) *amplitude = amp;

	return true;
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Md::computeAmplitude(const DoubleArray& data, size_t i1,
                                             size_t i2, size_t si1, size_t si2,
                                             double offset, AmplitudeIndex* dt,
                                             AmplitudeValue* amplitude,
                                             double* period, double* snr) {

	double amax, Imax, ofs_sig, amp_sig;
	DoubleArrayPtr d;

    if ( *snr < aFile.SNR_MIN )
        SEISCOMP_DEBUG("%s computed SNR is under configured SNR MIN", AMPTAG);

	if ( _computeAbsMax ) {
		size_t imax = find_absmax(data.size(), data.typedData(), si1, si2, offset);
		amax = fabs(data[imax] - offset);
		dt->index = imax;
	}
	else {
		int lmin, lmax;
		find_minmax(lmin, lmax, data.size(), data.typedData(), si1, si2, offset);
		amax = data[lmax] - data[lmin];
		dt->index = (lmin + lmax) * 0.5;
		dt->begin = lmin - dt->index;
		dt->end = lmax - dt->index;
	}

	Imax = dt->index;

	SEISCOMP_DEBUG("%s Amplitude max: %.2f", AMPTAG, amax);

	//! searching for Coda second by second through the end of the window
	//! if snrMin config is not 0 (config file or waveform review window)
	//! TODO: elevate accuracy by using a nanometers scale (maybe)
	if ( _config.snrMin != 0 ) {

		unsigned int i = si1;
		bool hasEndSignal = false;
		double calculatedSnr = -1;

		for (i = (int) Imax; i < i2; i = i + 1 * (int) _stream.fsamp) {

			int window_end = i + 1 * (int) _stream.fsamp;
			d = static_cast<DoubleArray*>(data.slice(i, window_end));

			//! computes pre-arrival offset
			ofs_sig = d->median();

			//! computes rms after removing offset
			amp_sig = 2 * d->rms(ofs_sig);

			if ( amp_sig / *_noiseAmplitude <= _config.snrMin ) {
				SEISCOMP_DEBUG("%s End of signal found! (%.2f <= %.2f)", AMPTAG,
				    (amp_sig / *_noiseAmplitude), _config.snrMin);
				hasEndSignal = true;
				calculatedSnr = amp_sig / *_noiseAmplitude;
				break;
			}
		}

		if ( !hasEndSignal ) {
			SEISCOMP_ERROR("%s SNR stayed over configured SNR_MIN! (%.2f > %.2f), "
				"skipping magnitude calculation for this station", AMPTAG,
			    calculatedSnr, _config.snrMin);
			return false;
		}

		dt->index = i;
	}
	else dt->index = Imax;

	//amplitude->value = 2 * amp_sig; //! actually it would have to be max. peak-to-peak
	amplitude->value = amp_sig;

	if ( _streamConfig[_usedComponent].gain != 0.0 )
		amplitude->value /= _streamConfig[_usedComponent].gain;
	else {
		setStatus(MissingGain, 0.0);
		return false;
	}

	// Convert m/s to nm/s
	amplitude->value *= 1.E09;

	*period = dt->index - i1 + (_config.signalBegin * _stream.fsamp);

	SEISCOMP_DEBUG("%s calculated event amplitude = %.2f", AMPTAG, amplitude->value);
	SEISCOMP_DEBUG("%s calculated signal end at %.2f ms from P phase", AMPTAG, *period);

	return true;
}