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

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Mjma::deconvolveData(Response *resp,
                                             DoubleArray &data,
                                             int numberOfIntegrations) {
	Math::Restitution::FFT::TransferFunctionPtr tf =
		resp->getTransferFunction(numberOfIntegrations);

	if ( tf == NULL ) {
		setStatus(DeconvolutionFailed, 0);
		return false;
	}

	
	Math::SeismometerResponse::Seismometer5sec paz(Math::Velocity);
	Math::Restitution::FFT::PolesAndZeros seis5sec(paz);

	Math::Restitution::FFT::TransferFunctionPtr cascade =
		*tf / seis5sec;

	// Remove linear trend
	double m,n;
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	return Math::Restitution::transformFFT(data.size(), data.typedData(),
	                                         _stream.fsamp, cascade.get(),
	                                         _config.respTaper, _config.respMinFreq, _config.respMaxFreq);
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Md::deconvolveData(Response* resp,
                                           DoubleArray& data,
                                           int numberOfIntegrations) {
	if ( numberOfIntegrations < -1 )
		return false;

	SEISCOMP_DEBUG("Inside deconvolve function");

	double m, n;
	Math::Restitution::FFT::TransferFunctionPtr tf =
		resp->getTransferFunction(numberOfIntegrations < 0 ? 0 : numberOfIntegrations);

	if ( !tf )
		return false;

	Math::GroundMotion gm;

	if ( numberOfIntegrations < 0 )
		gm = Math::Displacement;
	else
		gm = Math::Velocity;

	Math::Restitution::FFT::TransferFunctionPtr cascade;
	Math::SeismometerResponse::WoodAnderson woodAndersonResp(gm, _config.woodAndersonResponse);
	Math::SeismometerResponse::Seismometer5sec seis5sResp(gm);
	Math::SeismometerResponse::L4C_1Hz l4c1hzResp(gm);

	Math::Restitution::FFT::PolesAndZeros woodAnderson(woodAndersonResp);
	Math::Restitution::FFT::PolesAndZeros seis5sec(seis5sResp);
	Math::Restitution::FFT::PolesAndZeros l4c1hz(l4c1hzResp);

	SEISCOMP_DEBUG("SEISMO = %d", aFile.SEISMO);

	switch ( aFile.SEISMO ) {
		case 1:
			cascade = *tf / woodAnderson;
		break;
		case 2:
			cascade = *tf / seis5sec;
		break;
		case 9:
			SEISCOMP_INFO("%s Applying filter L4C 1Hz to data", AMPTAG);
			cascade = *tf / l4c1hz;
		break;
		default:
			cascade = tf;
			SEISCOMP_INFO("%s No seismometer specified, no signal reconvolution performed", AMPTAG);
			return false;
		break;
	}

	// Remove linear trend
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	return Math::Restitution::transformFFT(data.size(), data.typedData(),
	    _stream.fsamp, cascade.get(), _config.respTaper, _config.respMinFreq,
	    _config.respMaxFreq);
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor::deconvolveData(Response *resp, DoubleArray &data,
                                        int numberOfIntegrations) {
	// Remove linear trend
	double m,n;
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	return resp->deconvolveFFT(data, _stream.fsamp, _config.respTaper,
	                           _config.respMinFreq, _config.respMaxFreq,
	                           numberOfIntegrations);
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor::deconvolveData(Response *resp, DoubleArray &data,
                                        int numberOfIntegrations) {
	// Remove linear trend
	double m,n;
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	bool ret = resp->deconvolveFFT(data, _stream.fsamp, _config.respTaper,
	                               _config.respMinFreq, _config.respMaxFreq,
	                               numberOfIntegrations < 0 ? 0 : numberOfIntegrations);

	if ( !ret )
		return false;

	// If number of integrations are negative, derive data
	while ( numberOfIntegrations < 0 ) {
		Math::Filtering::IIRDifferentiate<double> diff;
		diff.setSamplingFrequency(_stream.fsamp);
		diff.apply(data.size(), data.typedData());
		++numberOfIntegrations;
	}

	return true;
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Mjma::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;

	size_t imax = find_absmax(data.size(), data.typedData(), si1, si2, offset);
	amax = fabs(data[imax] - offset);
	dt->index = imax;

	if ( *_noiseAmplitude == 0. )
		*snr = 1000000.0;
	else
		*snr = amax / *_noiseAmplitude;

	if ( *snr < _config.snrMin ) {
		setStatus(LowSNR, *snr);
		return false;
	}

	*period = -1;

	amplitude->value = amax;

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

	// - convert to micrometer
	amplitude->value *= 1E06;

	// - estimate peak-to-peak from absmax amplitude
	amplitude->value *= 2.0;

	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;
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
void AmplitudeProcessor::prepareData(DoubleArray &data) {
	Sensor *sensor = _streamConfig[_usedComponent].sensor();

	// When using full responses then all information needs to be set up
	// correctly otherwise an error is set
	if ( _enableResponses ) {
		if ( !sensor ) {
			setStatus(MissingResponse, 1);
			return;
		}

		if ( !sensor->response() ) {
			setStatus(MissingResponse, 2);
			return;
		}

		// If the unit cannot be converted into the internal
		// enum (what basically means "unknown") then the deconvolution
		// cannot be correctly. We do not want to assume a unit here
		// to prevent computation errors in case of bad configuration.
		SignalUnit unit;
		if ( !unit.fromString(sensor->unit().c_str()) ) {
			// Invalid unit string
			setStatus(IncompatibleUnit, 2);
			return;
		}

		int intSteps = 0;
		switch ( unit ) {
			case Meter:
				intSteps = -1;
				break;
			case MeterPerSecond:
				break;
			case MeterPerSecondSquared:
				intSteps = 1;
				break;
			default:
				setStatus(IncompatibleUnit, 1);
				return;
		}

		if ( _responseApplied ) return;

		_responseApplied = true;

		if ( !deconvolveData(sensor->response(), _data, intSteps) ) {
			setStatus(DeconvolutionFailed, 0);
			return;
		}
	}
	else {
		// If the sensor is known then check the unit and skip
		// non velocity streams. Otherwise simply use the data
		// to be compatible to the old version. This will be
		// changed in the future and checked more strictly.
		if ( sensor ) {
			SignalUnit unit;
			if ( !unit.fromString(sensor->unit().c_str()) ) {
				// Invalid unit string
				setStatus(IncompatibleUnit, 4);
				return;
			}

			switch ( unit ) {
				case Meter:
					if ( _enableUpdates ) {
						// Updates with differentiation are not yet supported.
						setStatus(IncompatibleUnit, 5);
						return;
					}

					// Derive to m/s
					{
						Math::Filtering::IIRDifferentiate<double> diff;
						diff.setSamplingFrequency(_stream.fsamp);
						diff.apply(data.size(), data.typedData());
					}
					break;

				case MeterPerSecond:
					break;

				default:
					setStatus(IncompatibleUnit, 3);
					return;
			}
		}
	}
}