DoubleSignal is a vector of double floating points with a sampling frequency, and a a type. More...

#include <DoubleSignal.h>

Inheritance diagram for GeopsyCore::DoubleSignal:

Public Types
enum	SignalType { Waveform = 0, Spectrum = 1, Phase = 2, RealSpectrum = 3, CAmpWaveform = 4, CPhaseWaveform = 5, ArrivalTime = 6, UndefinedSignalType = 7 }
Public Member Functions
bool	abs ()
int	add (const DoubleSignal *sig, double sigStart, double thisStart, double tToCopy, bool checkInvalid=false, double invalidValue=1e99)
int	add (const DoubleSignal *sig)
bool	agc (double width)
double	amplitude (const double *samples, int i) const
double	amplitude2 (const double *samples, int i) const
double	amplitudeAt (double abscissa) const
double	averageAmplitude (int itmin, int itmax) const
void	clipMaxToKeep (double percentage, KeepSignal *keep, double delay) const
Complex	complex (const double *samples, int i) const
ComplexSignal *	complexSpectrum () const
void	constLockProperties () const
bool	convolution (const DoubleSignal s1, const DoubleSignal s2)
void	copyAmplitudeFrom (const DoubleSignal *sig)
void	copyBasicProperties (const DoubleSignal &o)
void	copySamplesFrom (const DoubleSignal *sig)
void	copySamplesFrom (const DoubleSignal *sig, double sigStart, double thisStart, double tToCopy)
void	copySamplesFrom (const DoubleSignal *sig, int isigStart, int ithisStart, int nToCopy)
bool	correlation (const DoubleSignal s1, const DoubleSignal s2, double maxDelay)
double	correlation (const DoubleSignal *sig, double sigStart, double thisStart, double length) const
void	debugPrint () const
bool	decimateAmplitude (double delta)
bool	decimateTime (int factor, const DoubleSignal *sig)
double	deltaF () const
double	deltaT () const
bool	discreteFourierTransform ()
	DoubleSignal ()
	DoubleSignal (int n)
	DoubleSignal (const DoubleSignal &o)
double	duration () const
double	energy (double minF, double maxF) const
void	fastFourierTransform (SignalType st)
bool	filter (const FilterParameters &param)
bool	hv (const DoubleSignal zSig, const DoubleSignal hSig)
double	im (const double *samples, int i) const
double	interpoleAt (double abscissa) const
bool	isComplex () const
bool	isReal () const
bool	isSpectrum () const
void	lockProperties () const
double	maximumAmplitude (int itmin, int itmax) const
int	maximumAmplitudeAt (int itmin, int itmax) const
ComplexSignal *	morletWavelet (const MorletParameters &param) const
bool	normalizedCorrelation (const DoubleSignal s1, const DoubleSignal s2, double maxDelay)
bool	overSample (double factor, const DoubleSignal *sig)
double	phase (const double *samples, int i) const
ComplexSignal *	positiveComplexSpectrum (bool fullSize=false) const
double	power (double excludeBorderRatio=0.0) const
double	re (const double *samples, int i) const
void	restoreType (SignalType st)
void	restoreType (const DoubleSignal *sig) const
double	samplingFrequency () const
SignalType	saveType (SignalType st)
const DoubleSignal *	saveType (SignalType st) const
void	setAmplitude (double *samples, int i, double val) const
void	setComplex (double *samples, int i, const Complex &value) const
void	setDeltaF (double newval)
virtual void	setDeltaT (double newval)
void	setIm (double *samples, int i, double value) const
virtual void	setNSamples (int n)
void	setRe (double *samples, int i, double value) const
void	setType (SignalType t)
bool	setValue (double val)
bool	shift (double dt)
bool	smooth (const SmoothingParameters &param, double minFreq, double maxFreq)
bool	sqrt ()
bool	square ()
void	stalta (double tsta, double tlta, DoubleSignal *ratio) const
void	staltaToKeep (const WindowingParameters::ShortLongTermAverage &param, KeepSignal *keep, double delay) const
bool	stddevClip (double factor)
bool	subtractValue (double val=0.0)
bool	taper (const TimeRange &t, const TaperParameters &param)
bool	taper (const TaperParameters &param)
bool	timeCorrelation (const DoubleSignal s1, const DoubleSignal s2, double maxDelay)
SignalType	type () const
bool	unglitch (double threshold)
void	unlockProperies () const
bool	whiten ()
Static Public Member Functions
static SignalType	type (QString t)
static QString	typeString (SignalType st)
Protected Attributes
double	_deltaT
int	_nyquistIndex
QReadWriteLock	_propertiesMutex
SignalType	_type

Detailed Description

DoubleSignal is a vector of double floating points with a sampling frequency, and a a type.

DoubleSignal has only three additional public data members compared to SignalTemplate<double>:

sampling period (deltaT())
type: Waveform,Spectrum,Phase,RealSpectrum, CAmpWaveform, CPhaseWaveform, ArrivalTime
[internal] the last spectrum index if the number of samples is even

fft can be applied to this object as well as all other basic signal processings.

Member Enumeration Documentation

enum GeopsyCore::DoubleSignal::SignalType

Enumerator:

Waveform
Spectrum
Phase
RealSpectrum
CAmpWaveform
CPhaseWaveform
ArrivalTime
UndefinedSignalType

                  {
    Waveform=0,
    Spectrum=1,
    Phase=2,
    RealSpectrum=3,   // For real spectrum and their fft's
    CAmpWaveform=4,
    CPhaseWaveform=5,
    ArrivalTime=6,
    UndefinedSignalType=7
  };

Constructor & Destructor Documentation

GeopsyCore::DoubleSignal::DoubleSignal ( )

Construct a null signal (null sampling frequency, null number of samples).

References _deltaT, _nyquistIndex, _type, TRACE, and Waveform.

{
  TRACE;
  _type=Waveform;
  _deltaT=0.0;
  _nyquistIndex=0; // Set as an odd number of samples
}

GeopsyCore::DoubleSignal::DoubleSignal ( int n )

Construct a signal with n samples (null sampling frequency).

References _deltaT, GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, _type, TRACE, and Waveform.

                                :
    SignalTemplate<double>(n)
{
  TRACE;
  _type=Waveform;
  _deltaT=0.0;
  _nyquistIndex=_nSamples >> 1;
}

GeopsyCore::DoubleSignal::DoubleSignal ( const DoubleSignal & o )

Construct a a copy of signal o.

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, copyBasicProperties(), and TRACE.

                                                :
    SignalTemplate<double>(o)
{
  TRACE;
  _nyquistIndex=_nSamples >> 1;
  copyBasicProperties(o);
}

Member Function Documentation

bool DoubleSignal::abs ( )

Convert all samples to positive numbers

References LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::StationSignals::setKeep(), and GeopsyCore::SubSignalPool::stalta().

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    for(int i=0;i < _nSamples;i++ )
      thisSamples[ i ]=fabs(thisSamples[ i ] );
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

int GeopsyCore::DoubleSignal::add	(	const DoubleSignal *	sig,
		double	sigStart,
		double	thisStart,
		double	tToCopy,
		bool	checkInvalid = `false`,
		double	invalidValue = `1e99`
	)

Referenced by HVStationSignals::horizontal().

int DoubleSignal::add ( const DoubleSignal * sig )

References deltaT(), QGpCoreTools::endl(), QGpCoreTools::tr(), and TRACE.

{
  TRACE;
  if(_deltaT<=0) {
    App::stream() << tr("DoubleSignal::add: null deltaT, aborting.") << endl;
    return 0;
  }
  if(_deltaT!= sig->deltaT()) {
    App::stream() << tr("DoubleSignal::add: deltaT must be equal, aborting.") << endl;
    return 0;
  }
  SignalType st=saveType(Waveform);
  int nCommon=SignalTemplate<double>::add (sig);
  restoreType(st);
  return nCommon;
}

bool DoubleSignal::agc ( double width )

Automatic Gain Control

width is in seconds, it is the half width of the time window interval used to calculate the energy level.

Return false if not processed due to memory problem

References LOCK_SAMPLES, QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::agc().

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    SignalType st=saveType(Waveform);
    int iWidth=(int) floor(width/_deltaT + 0.5);
    double sum=0.0, val, nwin;
    int i;
    // Large windows (larger than signal length)
    if(2 * iWidth + 1 >= _nSamples) {
      for(i=0; i < _nSamples; ++i) {
        val=thisSamples[ i ];
        sum += val * val;
      }
      if(sum > 0) {
        sum=::sqrt(sum);
        for(i=0; i < _nSamples; ++i) {
          thisSamples[ i ] /= sum;
        }
      }
    } else {
      double * agcSig=new double [ _nSamples ];
      // compute initial window for first sample (only on the right side)
      for(i=0; i < iWidth + 1; ++i) {
        val=thisSamples[ i ];
        sum += val * val;
      }
      nwin=2.0 * iWidth + 1.0;
      agcSig[ 0 ]=(sum <= 0.0) ? 0.0 : thisSamples[ 0 ]/::sqrt(sum/nwin);

      // Next samples, increasing left side of the window
      for(i=1; i <= iWidth; ++i) {
        val=thisSamples[ i + iWidth ];
        sum += val * val;
        ++nwin;
        agcSig[ i ]=(sum <= 0.0) ? 0.0 : thisSamples[ i ]/::sqrt(sum/nwin);
      }

      // Full window available (2*width+1)
      int n=_nSamples - 1 - iWidth;
      for(i=iWidth + 1; i <= n; ++i) {
        val=thisSamples[ i + iWidth ];
        sum += val * val;
        val=thisSamples[ i - iWidth ];
        sum -= val * val;
        agcSig[ i ]=(sum <= 0.0) ? 0.0 : thisSamples[ i ]/::sqrt(sum/nwin);
      }

      // Decreasing right side of the window
      for(i=_nSamples - iWidth; i < _nSamples; ++i) {
        val=thisSamples[ i - iWidth ];
        sum -= val * val;
        --nwin;
        agcSig[ i ]=(sum <= 0.0) ? 0.0 : thisSamples[ i ]/::sqrt(sum/nwin);
      }

      // Copying back to original signal
      memcpy(thisSamples, agcSig, _nSamples * sizeof(double) );
      delete [] agcSig;
    }

    restoreType(st);
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

double GeopsyCore::DoubleSignal::amplitude	(	const double *	samples,
		int	i
	)		const `[inline]`

Returns the amplitude for frequency index i. samples must be a pointer to signal samples locked with LOCK_SAMPLES or CONST_LOCK_SAMPLES.

See also:: re(), im(), complex(), amplitude2(), phase()

References amplitude2(), and sqrt().

Referenced by copyAmplitudeFrom(), GeopsyGui::SignalLayer::drawSignal(), GeopsyGui::SignalLayer::updateGrid(), and GeopsyCore::Signal::writeAscii().

{
  return ::sqrt(amplitude2(samples, i));
}

double GeopsyCore::DoubleSignal::amplitude2	(	const double *	samples,
		int	i
	)		const `[inline]`

Returns the squared amplitude for frequency index i. samples must be a pointer to signal samples locked with LOCK_SAMPLES or CONST_LOCK_SAMPLES.

See also:: re(), im(), complex(), amplitude(), phase()

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and TRACE.

Referenced by ArrayCore::FKStationSignals::absolutePower(), amplitude(), and energy().

{
  TRACE;
  if(i==0) {
    return samples[i]*samples[i];
  } else if(i!=_nyquistIndex || (_nSamples & 0x00000001)) {
    return samples[i]*samples[i] +
           samples[_nSamples-i]*samples[_nSamples-i];
  } else { // For even sequency, sample at nyquist frequency
    return samples[i]*samples[i];
  }
}

double DoubleSignal::amplitudeAt ( double abscissa ) const

Returns the amplitude of the signal at time or frequency abscissa (depends upon signal type). The index is rounded to the lowest sample. No interpolation is perfomed.

Reimplemented in GeopsyCore::Signal.

References CONST_LOCK_SAMPLES, phase(), TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  double val;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    switch (_type) {
    case Waveform:
    case RealSpectrum: {
        int i=(int) (abscissa/_deltaT);
        if(i < 0 || i >= _nSamples)
          val=0.0;
        else
          val=thisSamples[ i ];
      }
      break;
    case Spectrum: {
        int i=(int) (abscissa * _deltaT * (double) _nSamples);
        if(i < 0 || i > (_nSamples >> 1) )
          val=0.0;
        else
          val=amplitude(thisSamples, i);
      }
      break;
    case Phase: {
        int i=(int) (abscissa * _deltaT * (double) _nSamples);
        if(i < 0 || i > (_nSamples >> 1) )
          val=0.0;
        else
          val=phase(thisSamples, i);
      }
      break;
    default:
      val=0.0;
      break;
    }
  UNLOCK_SAMPLES(this)
  else {
    val=0.0;
  }
  return val;
}

double DoubleSignal::averageAmplitude	(	int	itmin,
		int	itmax
	)		const

Returns the average amplitude of the signal between itmin and itmax

Reimplemented in GeopsyCore::Signal.

References CONST_LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  if(itmin < 0) itmin=0;
  if(itmax >= _nSamples) itmax=_nSamples - 1;
  double nSamples=(double) (itmax-itmin+1);

  double sum=0.0;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    switch (_type) {
    case Waveform:
    case RealSpectrum:
      for(int i=itmin; i <= itmax; i++ ) {
        sum += thisSamples[ i ];
      }
      break;
    case CAmpWaveform:
    case CPhaseWaveform:
    case Spectrum:
    case Phase: {
        int nSamples2=_nSamples >> 1;
        int it2max;
        if(itmax > nSamples2) it2max=nSamples2; else it2max=itmax;
        for(int i=itmin; i <= it2max; i++ ) {
          sum += amplitude2(thisSamples, i);
        }
      }
      break;
    default:
      sum=1.0;
      break;
    }
  UNLOCK_SAMPLES(this)
  return sum/nSamples;
}

void DoubleSignal::clipMaxToKeep	(	double	percentage,
		KeepSignal *	keep,
		double	delay
	)		const

Set keep signal according to the amplitude of signal (if above percentage of maxValue()). delay is the delay of this sig relative to the beginning of keep signal. keep and this signals must have the same deltaT().

References CONST_LOCK_SAMPLES, GeopsyCore::KeepSignal::deltaT(), QGpCoreTools::endl(), LOCK_SAMPLES, GeopsyCore::SignalTemplate< sampleType >::nSamples(), QGpCoreTools::tr(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::StationSignals::setKeep().

{
  TRACE;
  if(_deltaT<=0) {
    App::stream() << tr("DoubleSignal::clipMaxToKeep: null deltaT, aborting.") << endl;
    return;
  }
  if(_deltaT!=keep->deltaT()) {
    App::stream() << tr("DoubleSignal::clipMaxToKeep: deltaT must be equal, aborting.") << endl;
    return;
  }
  int thisStart, keepStart;
  int nSamples=commonSamples(delay, keep->nSamples(), thisStart, keepStart);
  CacheProcess cp;
  cp << this << keep;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    LOCK_SAMPLES(int, keepSamples, keep)
      double th=maximum() * percentage * 0.01;
      for(int i=0;i < nSamples;i++ ) {
        if(fabs( thisSamples[ i+thisStart ] ) >= th)
          keepSamples[ i+keepStart ]=0;
      }
    UNLOCK_SAMPLES(keep)
  UNLOCK_SAMPLES(this)
}

Complex GeopsyCore::DoubleSignal::complex	(	const double *	samples,
		int	i
	)		const `[inline]`

Retruns the complex sample with index i. samples must be a pointer to signal samples locked with LOCK_SAMPLES or CONST_LOCK_SAMPLES.

See also:: re(), im(), amplitude2(), amplitude(), phase()

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and TRACE.

Referenced by phase().

{
  TRACE;
  if(i==0) {
    return samples[i];
  } else if(i<_nyquistIndex) {
    return Complex(samples[i], samples[_nSamples-i]);
  } else if(i==_nyquistIndex && !(_nSamples & 0x00000001) ) {
    return samples[i]; // For even sequency, sample at nyquist frequency
  } else {
    return Complex(samples[_nSamples-i], -samples[i]);
  }
}

ComplexSignal * GeopsyCore::DoubleSignal::complexSpectrum ( ) const

Return the complex spectrum for positive and negative frequencies of a real signal.

(from FFTW documentation)) For those who like to think in terms of positive and negative frequencies, this means that the positive frequencies are stored in the first half of the output and the negative frequencies are stored in backwards order in the second half of the output. (The frequency -k/n is the same as the frequency (n-k)/n.)

See also:: re()

References GeopsyCore::SignalTemplate< double >::_nSamples, QGpCoreTools::conjugate(), CONST_LOCK_SAMPLES, GeopsyCore::csig, LOCK_SAMPLES, restoreType(), saveType(), Spectrum, TRACE, and UNLOCK_SAMPLES.

Referenced by convolution().

{
  TRACE;
  int n=_nSamples >> 1;
  // Copying to a full complex representation of the signal
  ComplexSignal * csig=new ComplexSignal(_nSamples);
  LOCK_SAMPLES(Complex, csamples, csig)
    // Assign values to the complex signal spectra from the real2complex representation
    const DoubleSignal * sig=saveType(Spectrum);
    CONST_LOCK_SAMPLES(double, sigSamples, sig)
      csamples[0]=sigSamples[ 0 ];
      for(int i=1; i<n; i++) {
        csamples[i].setRe(sigSamples[i]);              // positive frequency
        csamples[i].setIm(sigSamples[_nSamples-i]);
        csamples[_nSamples-i]=conjugate(csamples[i]);  // negative frequency
      }
      if(_nSamples & 0x00000001) {
        csamples[n].setRe(sigSamples[n]);
        csamples[n].setIm(sigSamples[_nSamples-n]);
      } else {
        csamples[n]=sigSamples[n];
      }
    UNLOCK_SAMPLES(sig)
    restoreType(sig);
  }

void GeopsyCore::DoubleSignal::constLockProperties ( ) const [inline]

{_propertiesMutex.lockForRead();}

bool DoubleSignal::convolution	(	const DoubleSignal *	s1,
		const DoubleSignal *	s2
	)

Calculates the convolution between signal s1 and s2 and stores the result in this signal.

This is the fast alternative based on Fourier's transform.

References complexSpectrum(), CONST_LOCK_SAMPLES, deltaT(), QGpCoreTools::endl(), QGpCoreTools::Complex::im, LOCK_SAMPLES, GeopsyCore::SignalTemplate< sampleType >::nSamples(), QGpCoreTools::Complex::re, QGpCoreTools::Complex::set(), sOut(), QGpCoreTools::Complex::toString(), QGpCoreTools::tr(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::convolution().

{
  TRACE;
  // Samples of this cannot be allocated
  if(isAllocated()) {
    App::stream() << tr("DoubleSignal::convolution: samples of destination signal cannot be allocated, aborting.") << endl;
    return false;
  }
  if(s1->deltaT()<=0) {
    App::stream() << tr("DoubleSignal::convolution: null deltaT, aborting.") << endl;
    return false;
  }
  if(s1->deltaT()!=s2->deltaT()) {
    App::stream() << tr("DoubleSignal::convolution: deltaT must be equal, aborting.") << endl;
    return false;
  }
  int n=s1->nSamples();
  if(n!=s2->nSamples()) {
    App::stream() << tr("DoubleSignal::convolution: number of samples must be equal, aborting.") << endl;
    return false;
  }
  int n2=n << 1;
  setNSamples(n2);
  setDeltaT(s1->deltaT());
  setType(Spectrum);

  Complex * t=new Complex[n];
  CONST_LOCK_SAMPLES(double, s1Samples, s1)
    for(int i=0; i<n; i++) {
      t[i].set(s1Samples[i], 0);
    }
  UNLOCK_SAMPLES(s1)
  FastFourierTransform::instance()->forward(n, t);
  QTextStream sOut(stdout);
  for(int i=0; i<n; i++) {
    sOut << t[i].toString(6) << endl;
  }
  sOut << "end fft complex" << endl;
  ComplexSignal * cspec=s1->complexSpectrum();
  CONST_LOCK_SAMPLES(Complex, cSamples, cspec)
  for(int i=0; i<n; i++) {
    sOut << cSamples[i].toString(6) << endl;
  }
  UNLOCK_SAMPLES(cspec)
  sOut << "end r2c rearranged" << endl;
  /*const DoubleSignal * spec1=s1->saveType(Spectrum);
  CONST_LOCK_SAMPLES(double, spec1Samples, spec1)
    for(int i=0; i<n/2; i++) {
      sOut << complex(spec1Samples, i).toString(6) << endl;
    }
  UNLOCK_SAMPLES(s1)
  const DoubleSignal * spec2=s2->saveType(Spectrum);*/



  bool ret=false;
#if 0
  LOCK_SAMPLES(double, convSamples, (&conv))
    CONST_LOCK_SAMPLES(double, s1Samples, spec1)
      CONST_LOCK_SAMPLES(double, s2Samples, spec2)
        convSamples[0]=s1Samples[0]*s2Samples[0];
        Complex c1, c2;
        for(int i=1; i< conv._nyquistIndex; i++ ) {
          c1.set(s1Samples[ i ], s1Samples[ n2 - i ] );
          c2.set(s2Samples[ i ], -s2Samples[ n2 - i ] );
          c1*=c2;
          convSamples[i]=c1.re();
          convSamples[n2 - i]=c1.im();
        }
        if( !(conv._nSamples & 0x00000001) ) {
          convSamples[conv._nyquistIndex]=s1Samples[conv._nyquistIndex]*s2Samples[conv._nyquistIndex];
        }
        ret=true;
      UNLOCK_SAMPLES(spec2)
    UNLOCK_SAMPLES(spec1)
    conv.fastFourierTransform(Waveform);
    LOCK_SAMPLES(double, thisSamples, this)
      // Negative lags
      int n0=n2-_nyquistIndex;
      for(int i=0; i< _nyquistIndex; i++ ) {
        thisSamples[i]=convSamples[n0+i];
      }
      // Positive lags
      n0=-_nyquistIndex;
      for(int i=_nyquistIndex; i< _nSamples; i++ ) {
        thisSamples[i]=convSamples[n0+i];
      }
    UNLOCK_SAMPLES(this)
  UNLOCK_SAMPLES((&conv))
#endif
  //delete spec1;
  //delete spec2;
  return ret;
}

void DoubleSignal::copyAmplitudeFrom ( const DoubleSignal * sig )

Copy all amplitudes from signal sig. Hence transform a complex spectrum into a real amplitude spectrum. This signal must not be allocated.

References amplitude(), CONST_LOCK_SAMPLES, deltaT(), QGpCoreTools::endl(), LOCK_SAMPLES, GeopsyCore::SignalTemplate< sampleType >::nSamples(), restoreType(), saveType(), QGpCoreTools::tr(), TRACE, type(), and UNLOCK_SAMPLES.

Referenced by HVStationSignals::horizontal(), StructureStationProcessSignals::setProcessed(), HVStationSignals::setProcessed(), and MonoStation::SpectrumStationSignals::setProcessed().

{
  TRACE;
  if(isAllocated()) {
    App::stream() << tr("DoubleSignal::copyAmplitudeFrom: samples of destination signal cannot be allocated, aborting.") << endl;
    return;
  }
  setDeltaT(1.0/(sig->deltaT() * sig->nSamples()) );
  int nSamples2=sig->nSamples() >> 1;
  setNSamples(nSamples2+1);

  const DoubleSignal * sigfft=sig->saveType(Spectrum);
  CONST_LOCK_SAMPLES(double, sigSamples, sigfft)
    LOCK_SAMPLES(double, thisSamples, this)
      for(int i=0;i <= nSamples2;i++ ) {
        thisSamples[ i ]=sig->amplitude(sigSamples, i);
      }
    UNLOCK_SAMPLES(this)
  UNLOCK_SAMPLES(sigfft)
  sig->restoreType(sigfft);

  switch (sig->type()) {
  case Spectrum:
  case Phase:
    _type=RealSpectrum;
    break;
  case CAmpWaveform:
  case CPhaseWaveform:
    _type=Waveform;
  default:
    break;
  }
}

void GeopsyCore::DoubleSignal::copyBasicProperties ( const DoubleSignal & o )

References _deltaT, _type, and TRACE.

Referenced by DoubleSignal().

{
  TRACE;
  _deltaT=o._deltaT;
  _type=o._type;
}

void GeopsyCore::DoubleSignal::copySamplesFrom ( const DoubleSignal * sig ) [inline]

Referenced by TFAResults::compute(), DampingResults::compute(), GeopsyCore::StationProcessSignals::copyOriginalSignal(), GeopsyCore::Signal::copySamplesFrom(), decimateTime(), GeopsyCore::StationSignals::setKeep(), GeopsyCore::SubSignalPool::stalta(), and T0Correlation::value().

{
  SignalTemplate<double>::copySamplesFrom(sig);
}

void DoubleSignal::copySamplesFrom	(	const DoubleSignal *	sig,
		double	sigStart,
		double	thisStart,
		double	tToCopy
	)

References deltaT(), QGpCoreTools::endl(), QGpCoreTools::tr(), and TRACE.

{
  TRACE;
  if(_deltaT<=0) {
    App::stream() << tr("DoubleSignal::copySamplesFrom: null deltaT, aborting.") << endl;
    return;
  }
  if(_deltaT!=sig->deltaT()) {
    App::stream() << tr("DoubleSignal::copySamplesFrom: deltaT must be equal, aborting.") << endl;
    return;
  }
  double fSamp=1.0/_deltaT;
  copySamplesFrom(sig, (int) round(sigStart * fSamp),
                           (int) round(thisStart * fSamp),
                           (int) round(tToCopy * fSamp) );
}

void DoubleSignal::copySamplesFrom	(	const DoubleSignal *	sig,
		int	isigStart,
		int	ithisStart,
		int	nToCopy
	)

References restoreType(), saveType(), and TRACE.

{
  TRACE;
  const DoubleSignal * wavSig=sig->saveType(Waveform);
  SignalType st=saveType(Waveform);
  SignalTemplate<double>::copySamplesFrom(wavSig, isigStart, ithisStart, nToCopy);
  restoreType(st);
  sig->restoreType(wavSig);
}

bool DoubleSignal::correlation	(	const DoubleSignal *	s1,
		const DoubleSignal *	s2,
		double	maxDelay
	)

Calculates the unormalized time correlation between signal s1 and s2 and stores the result in this signal.

This is the fast method based on Fourier's transform.

References TRACE.

Referenced by GeopsyCore::Signal::correlation(), and T0Correlation::value().

{
  TRACE;
  if(!correlationCheck(s1, s2, maxDelay)) {
    return false;
  }
  setNSamples((int)round(2*maxDelay/s1->deltaT())+1);
  setDeltaT(s1->deltaT());

  DoubleSignal * s1t=new DoubleSignal(*s1);
  DoubleSignal * s2t=new DoubleSignal(*s2);
  bool ret=correlationProcess(s1, s2, s1t, s2t);
  delete s1t;
  delete s2t;
  return ret;
}

double GeopsyCore::DoubleSignal::correlation	(	const DoubleSignal *	sig,
		double	sigStart,
		double	thisStart,
		double	length
	)		const

void DoubleSignal::debugPrint ( ) const

References CONST_LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  printf( "Debug samples for signal %s\n", debugName().toAscii().data());
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    for(int i=0;i < _nSamples;i++ ) printf( "%i\t%lg\n", i, thisSamples[ i ] );
  UNLOCK_SAMPLES(this)
}

bool DoubleSignal::decimateAmplitude ( double delta )

Decimate along amplitude scale. delta is the amplitude step used to discretize amplitudes.

References LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::decimateAmplitude().

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    SignalType st=saveType(Waveform); // make sure we are in time domain
    double invDelta=1.0/delta;
    for(int i=0; i<_nSamples; i++) {
      if(thisSamples[i]>0)
        thisSamples[i]=ceil(thisSamples[i] * invDelta) * delta;
      else
        thisSamples[i]=floor(thisSamples[i] * invDelta) * delta;
    }
    restoreType(st);
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

bool DoubleSignal::decimateTime	(	int	factor,
		const DoubleSignal *	sig
	)

Decimate sig along time scale. The result is set to this. An adequat anti-aliasing filter is used before any processing. One sample over factor is kept.

sig is the reference signal. The results are stored in this signal. This signal must not be allocated.

References GeopsyCore::FilterParameters::convolutionWindow(), copySamplesFrom(), QGpCoreTools::endl(), fastFourierTransform(), filter(), LOCK_SAMPLES, GeopsyCore::SignalTemplate< sampleType >::nSamples(), GeopsyCore::FilterParameters::setBand(), GeopsyCore::FilterParameters::setMethod(), GeopsyCore::FilterParameters::setMinimumFrequency(), GeopsyCore::FilterParameters::setWidth(), GeopsyCore::TaperParameters::setWindow(), QGpCoreTools::tr(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::decimateTime().

{
  TRACE;
  if(factor<=1) {
    App::stream() << tr("DoubleSignal::decimateTime: factor cannot be less than 1, aborting.") << endl;
    return false;
  }
  if(isAllocated()) {
    App::stream() << tr("DoubleSignal::decimateTime: samples of destination signal cannot be allocated, aborting.") << endl;
    return false;
  }
  int nResampled=(int) floor(sig->nSamples()/(double)factor);
  setNSamples(nResampled);
  setDeltaT(factor * deltaT());

  // Anti aliasing filter, need to make a copy signal
  DoubleSignal * filtsig=new DoubleSignal(*sig);
  filtsig->copySamplesFrom(sig);
  // Use an anti-aliasing filter: Low Pass with cut-off frequency set to 0.8*Nyquist frequency
  FilterParameters param;
  param.setMinimumFrequency(0.5/_deltaT);
  param.setWidth(0.2);
  param.setBand(FilterParameters::LowPass);
  param.setMethod(FilterParameters::Convolution);
  param.convolutionWindow().setWindow(TaperParameters::Cosine);
  if( !filtsig->filter(param) ) {
    delete filtsig;
    return false;
  }
  filtsig->fastFourierTransform(Waveform);
  bool ret=false;
  CacheProcess cp;
  cp << filtsig << this;
  LOCK_SAMPLES(double, sigSamples, filtsig)
    LOCK_SAMPLES(double, thisSamples, this)
      for(int i=0; i<_nSamples; i++) {
        thisSamples[i]=sigSamples[ i*factor ];
      }
      ret=true;
    UNLOCK_SAMPLES (this)
  UNLOCK_SAMPLES(filtsig)
  delete filtsig;
  return ret;
}

double GeopsyCore::DoubleSignal::deltaF ( ) const [inline]

Referenced by GeopsyCore::SubSignalPool::timeRange(), GeopsyGui::SignalLayer::updateGrid(), and GeopsyGui::SignalLayer::visibleSamples().

{return _deltaT;}

double GeopsyCore::DoubleSignal::deltaT ( ) const [inline]

{return _deltaT;}

bool DoubleSignal::discreteFourierTransform ( )

References QGpCoreTools::cos(), LOCK_SAMPLES, QGpCoreTools::sin(), and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::discreteFourierTransform().

{
  long i,k;
  double arg;
  double cosarg,sinarg;

  int n=nSamples();
  double * x2=new double[n];
  double * y2=new double[n];

  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)

  for(i=0;i<n;i++) {
    x2[i]=0;
    y2[i]=0;
    arg=- 2.0 * M_PI * (double)i/(double)n;
    for(k=0;k<n;k++) {
      cosarg=cos(k * arg);
      sinarg=sin(k * arg);
      x2[i] += (thisSamples[k] * cosarg);
      y2[i] += (thisSamples[k] * sinarg);
    }
  }

  // Compute spectrum amplitude
  // From 0 to n/2 : positive frequencies
  // From n/2 to n-1: negative frequencies in reverse order (antisymetric for real input signal
  int n2=n >> 1;
  double rn=1.0/(_nSamples*_deltaT);
  for(i=0;i<=n2;i++) {
    thisSamples[i]=x2[i]*rn;
  }
  for(i=n2+1;i<n;i++) {
    thisSamples[i]=y2[n-i]*rn;
  }
  ret=true;
  UNLOCK_SAMPLES(this);
  setType(Spectrum);

  delete [] x2;
  delete [] y2;
  return ret;
}

double GeopsyCore::DoubleSignal::duration ( ) const [inline]

Referenced by GeopsyCore::Signal::compare(), GeopsyCore::Signal::header(), GeopsyCore::CitySignal::loadSignals(), GeopsyCore::SignalHeaderObject::setEndTime(), and GeopsyGui::SignalLayer::visibleSamples().

{return _deltaT * _nSamples;}

double DoubleSignal::energy	(	double	minF,
		double	maxF
	)		const

Returns the energy between frequency minF and maxF.

References amplitude2(), CONST_LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by LinearFKActiveStationSignals::normalize().

{
  TRACE;
  double en=0.0;
  const DoubleSignal * sig=saveType(Spectrum); // make sure we are in frequency domain
  CONST_LOCK_SAMPLES(double, sigSamples, sig)
    double fac=_deltaT * (double) _nSamples;
    int minI=(int) floor(minF * fac);
    if(minI<0) minI=0;
    int maxI=(int) ceil(maxF * fac);
    if(maxI>=_nSamples) maxI=_nSamples-1;
    if(minI>=maxI) return 0.0;
    for(int i=minI; i<maxI; i++ ) {
      en+=sig->amplitude2(sigSamples, i);
    }
  UNLOCK_SAMPLES(sig);
  restoreType(sig);
  return en;
}

void DoubleSignal::fastFourierTransform ( SignalType st )

A r2r kind of FFTW_R2HC (r2hc) corresponds to an r2c DFT (see One-Dimensional DFTs of Real Data) but with "halfcomplex" format output, and may sometimes be faster and/or more convenient than the latter. The inverse hc2r transform is of kind FFTW_HC2R. This consists of the non-redundant half of the complex output for a 1d real-input DFT of size n, stored as a sequence of n real numbers (double) in the format:

r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1

Here, rk is the real part of the kth output, and ik is the imaginary part. (Division by 2 is rounded down.) For a halfcomplex array hc[n], the kth component thus has its real part in hc[k] and its imaginary part in hc[n-k], with the exception of k==0 or n/2 (the latter only if n is even)--in these two cases, the imaginary part is zero due to symmetries of the real-input DFT, and is not stored. Thus, the r2hc transform of n real values is a halfcomplex array of length n, and vice versa for hc2r.

References LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::StationSignals::addSignal(), StructureStationSignals::addSignals(), decimateTime(), LinearFKActiveStationSignals::endPreprocess(), GeopsyCore::SubSignalPool::fastFourierTransform(), LinearFKActiveStationSignals::normalize(), StructureStationProcessSignals::setProcessed(), RealTimeStationSignals::setProcessed(), MonoStation::SpectrumStationSignals::setProcessed(), and ArrayCore::FKStationSignals::setProcessed().

{
  TRACE;
  if(isComplex()) {
    switch (st) {
    case Waveform:
    case RealSpectrum: {
        LOCK_SAMPLES(double, thisSamples, this)
          FastFourierTransform::instance()->backward(_nSamples, thisSamples);
          for(int i=0; i<_nSamples; i++) {
            thisSamples[i]*=_deltaT;
          }
        UNLOCK_SAMPLES(this)
        switch (_type) {
        case CAmpWaveform:
        case CPhaseWaveform:
          _type=RealSpectrum;
          break;
        case Spectrum:
        case Phase:
          _type=Waveform;
        default:
          break;
        }
      }
      break;
    case Spectrum:
    case CAmpWaveform:
      if(_type==CPhaseWaveform)
        _type=CAmpWaveform;
      else if(_type==Phase)
        _type=Spectrum;
      break;
    case Phase:
    case CPhaseWaveform:
      if(_type==CAmpWaveform)
        _type=CPhaseWaveform;
      else if(_type==Spectrum)
        _type=Phase;
    case ArrivalTime:
    case UndefinedSignalType:
      break;
    }
  } else if(isReal()) {
    switch (st) {
    case Waveform:
    case RealSpectrum:
    case ArrivalTime:
    case UndefinedSignalType:
      break;
    case Spectrum:
    case Phase:
    case CAmpWaveform:
    case CPhaseWaveform: {
        LOCK_SAMPLES(double, thisSamples, this)
          FastFourierTransform::instance()->forward(_nSamples, thisSamples);
          double rn=1.0/(_nSamples*_deltaT);
          for(int i=0; i < _nSamples; i++ ) {
            thisSamples[i]*=rn;
          }
        UNLOCK_SAMPLES(this)
        switch (_type) {
        case RealSpectrum:
          if(st==Phase || st==CPhaseWaveform)
            _type=CPhaseWaveform;
          else
            _type=CAmpWaveform;
          break;
        case Waveform:
          if(st==Phase || st==CPhaseWaveform)
            _type=Phase;
          else
            _type=Spectrum;
        default:
          break;
        }
        break;
      }
    }
  }
}

bool GeopsyCore::DoubleSignal::filter ( const FilterParameters & param )

Butterworth filter, cosine filter, and Morlet wavelet filter

Return false if not processed due to memory problem

Referenced by DampingResults::compute(), decimateTime(), GeopsyCore::SubSignalPool::filter(), Process::run(), GeopsyCore::StationSignals::setKeep(), and ArrayCore::SPACStationSignals::setProcessed().

{
  TRACE;
  bool ret=true;
  switch (param.method()) {
  case FilterParameters::Butterworth: {
      SignalType ost=saveType(Waveform);
      LOCK_SAMPLES(double, thisSamples, this)
        switch (param.band()) {
        case FilterParameters::BandPass:
          filrec(param.maximumFrequency(), 1, 0, _nSamples-1, 1, thisSamples, thisSamples, param.order(), _deltaT);
          filrec(param.minimumFrequency(), -1, 0, _nSamples-1, 1, thisSamples, thisSamples, param.order(), _deltaT);
          break;
        case FilterParameters::BandReject: {
            double * s=new double [_nSamples];
            int i;
            for(i=0; i<_nSamples; i++) {
              s[i]=thisSamples[i];
            }
            filrec(param.minimumFrequency(), 1, 0, _nSamples-1, 1, thisSamples, thisSamples, param.order(), _deltaT);
            filrec(param.maximumFrequency(), -1, 0, _nSamples-1, 1, s, s, param.order(), _deltaT);
            for(i=0; i<_nSamples; i++) {
              thisSamples[ i ]+=s[i];
            }
            delete [] s;
          }
          break;
        case FilterParameters::LowPass:
          filrec(param.minimumFrequency(), 1, 0, _nSamples-1, 1, thisSamples, thisSamples, param.order(), _deltaT);
          break;
        case FilterParameters::HighPass:
          filrec(param.maximumFrequency(), -1, 0, _nSamples-1, 1, thisSamples, thisSamples, param.order(), _deltaT);
        default:
          break;
        }
      UNLOCK_SAMPLES(this) else ret=false;
      restoreType(ost);
    }
    break;
  case FilterParameters::Convolution: {
      SignalType st=saveType(Spectrum);
      double fac=_deltaT*(double)_nSamples;
      LOCK_SAMPLES(double, thisSamples, this)
        TaperDelegate::ComplexSample samples(_nSamples, thisSamples);
        TaperDelegate window(&param.convolutionWindow());
        switch (param.band()) {
        case FilterParameters::LowPass:
          window.apply(samples, 0, _nyquistIndex, fac*param.minimumFrequency()*(1.0-param.width()), fac*param.minimumFrequency());
          break;
        case FilterParameters::HighPass:
          window.apply(samples, 0, _nyquistIndex, fac*param.minimumFrequency(), fac*param.minimumFrequency()*(1.0+param.width()));
          break;
        case FilterParameters::BandPass:
          window.apply(samples, 0, _nyquistIndex, fac*param.minimumFrequency(), fac*param.maximumFrequency());
          break;
        case FilterParameters::BandReject:
          window.apply(samples, 0, _nyquistIndex, fac*param.minimumFrequency(), fac*param.maximumFrequency());
          break;
        }
      UNLOCK_SAMPLES(this) else ret=false;
      restoreType(st);
    }
    break;
  }
  return ret;
}

bool DoubleSignal::hv	(	const DoubleSignal *	zSig,
		const DoubleSignal *	hSig
	)

Compute the H/V spectra ratio of signals zSig (vertical component), hSig (horizontal component).

References CONST_LOCK_SAMPLES, deltaT(), QGpCoreTools::endl(), LOCK_SAMPLES, GeopsyCore::SignalTemplate< sampleType >::nSamples(), QGpCoreTools::tr(), TRACE, type(), and UNLOCK_SAMPLES.

{
  TRACE;
  // Samples of this cannot be allocated
  if(isAllocated()) {
    App::stream() << tr("DoubleSignal::hv: samples of destination signal cannot be allocated, aborting.") << endl;
    return false;
  }
  // Samples of zSig, h1Sig and h2Sig must be in real spectrum
  if(zSig->type()!=RealSpectrum || hSig->type()!=RealSpectrum) return false;
  _type=RealSpectrum;
  _deltaT=zSig->deltaT();
  setNSamples(zSig->nSamples());
  bool ret=false;
  CacheProcess cp;
  cp << zSig << hSig << this;
  CONST_LOCK_SAMPLES(double, zSamples, zSig)
    CONST_LOCK_SAMPLES(double, hSamples, hSig)
      LOCK_SAMPLES(double, thisSamples, this)
        for(int i=0;i < _nSamples;i++ ) {
          thisSamples[ i ]=hSamples[ i ]/zSamples[ i ];
        }
        ret=true;
      UNLOCK_SAMPLES(this)
    UNLOCK_SAMPLES(hSig)
  UNLOCK_SAMPLES(zSig)
  return ret;
}

double GeopsyCore::DoubleSignal::im	(	const double *	samples,
		int	i
	)		const `[inline]`

Retruns the imaginary part of sample with index i. samples must be a pointer to signal samples locked with LOCK_SAMPLES or CONST_LOCK_SAMPLES.

See also:: re(), complex(), amplitude2(), amplitude(), phase()

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and TRACE.

{
  TRACE;
  if(i==0) {
    return 0.0;
  } else if(i < _nyquistIndex) {
    return samples[ _nSamples -i ];
  } else if(i==_nyquistIndex && !(_nSamples & 0x00000001)) {
    return 0.0; // For even sequency, sample at nyquist frequency
  } else {
    return -samples[ i ];
  }
}

double DoubleSignal::interpoleAt ( double abscissa ) const

Returns the amplitude of the signal at time or frequency abscissa (depends upon signal type). The returned amplitude is interpolated.

References CONST_LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by MonoStation::StationResults::setWindow().

{
  TRACE;
  double val;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    switch (_type) {
    case Waveform:
    case RealSpectrum: {
        double di=abscissa/_deltaT;
        int i=(int) ceil(di);
        if(i < 0 || i >= _nSamples) val=0.0;
        else {
          if(i==0) i=1;
          val=thisSamples[ i ] + (thisSamples[ i ] - thisSamples[ i - 1 ] ) * (di - (double) i);
        }
      }
      break;
    default:
      val=0.0;
      break;
    }
  UNLOCK_SAMPLES(this)
  else {
    val=0.0;
  }
  return val;
}

bool GeopsyCore::DoubleSignal::isComplex ( ) const [inline]

References _type, CAmpWaveform, CPhaseWaveform, Phase, Spectrum, and TRACE.

Referenced by GeopsyCore::SubSignalPool::timeRange(), and GeopsyCore::Signal::writeAscii().

{
  TRACE;
  switch (_type) {
  case Spectrum:
  case Phase:
  case CAmpWaveform:
  case CPhaseWaveform:
    return true;
  default:
    return false;
  }
}

bool GeopsyCore::DoubleSignal::isReal ( ) const [inline]

References _type, RealSpectrum, TRACE, and Waveform.

Referenced by GeopsyCore::Signal::writeAscii().

{
  TRACE;
  switch (_type) {
  case Waveform:
  case RealSpectrum:
    return true;
  default:
    return false;
  }
}

bool GeopsyCore::DoubleSignal::isSpectrum ( ) const [inline]

References _type, Phase, RealSpectrum, Spectrum, and TRACE.

Referenced by GeopsyGui::SignalLayer::drawSignal(), GeopsyGui::PickLayer::mouseMoveEvent(), GeopsyGui::PickLayer::paintData(), GeopsyGui::SignalLayer::properties(), and GeopsyCore::SubSignalPool::timeRange().

{
  TRACE;
  switch (_type) {
  case Spectrum:
  case Phase:
  case RealSpectrum:
    return true;
  default:
    return false;
  }
}

void GeopsyCore::DoubleSignal::lockProperties ( ) const [inline]

{_propertiesMutex.lockForWrite();}

double DoubleSignal::maximumAmplitude	(	int	itmin,
		int	itmax
	)		const

Returns the maximum amplitude of the signal between itmin and itmax

Reimplemented in GeopsyCore::Signal.

References CONST_LOCK_SAMPLES, QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

Referenced by LinearFKActiveStationSignals::normalize().

{
  TRACE;
  if(itmin < 0) itmin=0;
  if(itmax >= _nSamples) itmax=_nSamples - 1;

  double vmax=0.0;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    double v;
    switch (_type) {
    case Waveform:
    case RealSpectrum:
      for(int i=itmin; i <= itmax; i++ ) {
        v=fabs(thisSamples[ i ] );
        if(v > vmax) vmax=v;
      }
      break;
    case CAmpWaveform:
    case CPhaseWaveform:
    case Spectrum:
    case Phase: {
        int nSamples2=_nSamples >> 1;
        int it2max;
        if(itmax > nSamples2) it2max=nSamples2; else it2max=itmax;
        for(int i=itmin; i <= it2max; i++ ) {
          v=amplitude2(thisSamples, i);
          if(v > vmax) vmax=v;
        }
        vmax=::sqrt(vmax);
      }
      break;
    default:
      vmax=1.0;
      break;
    }
  UNLOCK_SAMPLES(this)
  return vmax;
}

int DoubleSignal::maximumAmplitudeAt	(	int	itmin,
		int	itmax
	)		const

Returns the index of the maximum amplitude of the signal between itmin and itmax

References CONST_LOCK_SAMPLES, QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  if(itmin<0) itmin=0;
  if(itmax>=_nSamples) itmax=_nSamples-1;

  double vmax=0.0;
  int maxIndex=itmin;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    double v;
    switch (_type) {
    case Waveform:
    case RealSpectrum:
      for(int i=itmin; i<=itmax; i++) {
        v=fabs(thisSamples[i]);
        if(v>vmax) {
          vmax=v;
          maxIndex=i;
        }
      }
      break;
    case CAmpWaveform:
    case CPhaseWaveform:
    case Spectrum:
    case Phase: {
        int nSamples2=_nSamples >> 1;
        int it2max;
        if(itmax > nSamples2) it2max=nSamples2; else it2max=itmax;
        for(int i=itmin; i <= it2max; i++) {
          v=amplitude2(thisSamples, i);
          if(v>vmax) {
            vmax=v;
            maxIndex=i;
          }
        }
        vmax=::sqrt(vmax);
      }
      break;
    default:
      break;
    }
  UNLOCK_SAMPLES(this)
  return maxIndex;
}

ComplexSignal * GeopsyCore::DoubleSignal::morletWavelet ( const MorletParameters & param ) const

Return the convolution of this signal with Morlet Wavelet in freshly allocated complex signal

Input:

w0 factor > 5.5 to has an admissible wavelet
fmin, central frequency of filter (s=w0/(2*pi*fmin))

References _deltaT, GeopsyCore::SignalTemplate< double >::_nSamples, GeopsyCore::FastFourierTransform::backward(), GeopsyCore::csig, QGpCoreTools::exp(), GeopsyCore::MorletParameters::fi(), GeopsyCore::FastFourierTransform::instance(), LOCK_SAMPLES, GeopsyCore::MorletParameters::m(), GeopsyCore::SignalTemplate< sampleType >::nSamples(), positiveComplexSpectrum(), TRACE, and UNLOCK_SAMPLES.

Referenced by TFAResults::compute(), and GeopsyCore::SubSignalPool::waveletTransform().

{
  TRACE;
  // Number of samples that describe the spectrum for positive part
  int n=_nSamples >> 1;
  ComplexSignal * csig=positiveComplexSpectrum(true);
  LOCK_SAMPLES(Complex, csamples, csig)
    // Make convolution of this complex spectrum with the Morlet function
    // defined only for the positive frequencies.
    //
    // Morlet Wavelet:
    //      Modified Morlet wavelet=1/pow(pi,0.25) * exp (-(f*s - w0)^2*m)
    //      Common value for m is 10
    // s=w0/param.fmin
    // However for efficiency: s=param.w0/param.fmin*df (see computation of a in loop)
    double df=1.0/(_deltaT * (double) _nSamples);
    double s=6.0*df/param.fi();
    double c=pow(M_PI, -0.25);
    csamples[0]=0.0; // DC component
    for(int i=1;i<=n; i++) { // 0 frequency is already set to 0
      double a=s*i-6.0;
      a=a*a;
      if(a<700) {
        a=c*exp(-a*param.m());
        csamples[i]*= a;
      } else {
        csamples[i]=0.0;
      }
    }
    // Negative frequency: convolution gives null contribution
    for(int i=_nSamples-1;i>n; i--) {
      csamples[i]=0.0;
    }
    // Back transform to time domain
    FastFourierTransform::instance()->backward(csig->nSamples(), csamples);
  UNLOCK_SAMPLES(csig)
  return csig;
}

bool DoubleSignal::normalizedCorrelation	(	const DoubleSignal *	s1,
		const DoubleSignal *	s2,
		double	maxDelay
	)

Calculates the normalized time correlation between signal s1 and s2 and stores the result in this signal.

This is a fast alternative based on Fourier's transform.

This implementation is based on "Fast Normalized Cross-Correlation" by J. P. Lewis, downgraded to 1D signals (instead of 2D images). (http://scribblethink.org/Work/nvisionInterface/nip.html)

References CONST_LOCK_SAMPLES, deltaT(), LOCK_SAMPLES, QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::Signal::normalizedCorrelation().

{
  TRACE;
  if(!correlationCheck(s1, s2, maxDelay)) {
    return false;
  }
  setNSamples((int)round(2*maxDelay/s1->deltaT())+1);
  setDeltaT(s1->deltaT());

  DoubleSignal * s1t=new DoubleSignal(*s1);
  DoubleSignal * s2t=new DoubleSignal(*s2);
  if(!correlationProcess(s1, s2, s1t, s2t)) {
    delete s1t;
    delete s2t;
    return false;
  }

  // Computation of normalization factors
  DoubleSignal * s1Norm=new DoubleSignal(_nSamples);
  DoubleSignal * s2Norm=new DoubleSignal(_nSamples);
  CacheProcess cp;
  cp << this << s1Norm << s2Norm;
  correlationNormalization(s1t, s1Norm);
  correlationNormalization(s2t, s2Norm);
  delete s1t;
  delete s2t;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    CONST_LOCK_SAMPLES(double, s1Samples, s1Norm)
      CONST_LOCK_SAMPLES(double, s2Samples, s2Norm)
        for(int i=0; i<_nSamples; i++) {
          double en1=s1Samples[i];
          double en2=s2Samples[_nSamples-i-1];
          if(en1>0.0 && en2>0.0) {
            thisSamples[i]/=::sqrt(en1*en2);
          }
        }
        ret=true;
      UNLOCK_SAMPLES(s2Norm)
    UNLOCK_SAMPLES(s1Norm)
  UNLOCK_SAMPLES(this)
  delete s1Norm;
  delete s2Norm;
  return ret;
}

bool DoubleSignal::overSample	(	double	factor,
		const DoubleSignal *	sig
	)

Divide the sampling period by factor, equivalent to zero padding on the complex spectrum. sig is the reference signal. The results are stored in this signal. This signal must not be allocated.

Organisation of samples:

r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1

Definition of evenNyquistIndex:

if( !(n & 0x00000001) ) _evenNyquistIndex=_nSamples >> 1; else _evenNyquistIndex=0; // Set as an odd number of samples

Its use to get the complex spectrum:

if(i==0 || i==_evenNyquistIndex) return samples[ i ]; else return Complex(samples[ i ], samples[ _nSamples - i ] );

References GeopsyCore::SignalTemplate< sampleType >::_nSamples, _nyquistIndex, CONST_LOCK_SAMPLES, QGpCoreTools::endl(), LOCK_SAMPLES, restoreType(), saveType(), QGpCoreTools::tr(), TRACE, type(), and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::overSample().

{
  TRACE;
  if(factor <=1) {
    App::stream() << tr("DoubleSignal::overSample: factor cannot be less than 1, aborting.") << endl;
    return false;
  }
  if(isAllocated()) {
    App::stream() << tr("DoubleSignal::overSample: samples of destination signal cannot be allocated, aborting.") << endl;
    return false;
  }
  int nOriginalSpectrum=_nSamples >> 1;
  int nResampled=(int) round(_nSamples*factor);
  setNSamples(nResampled);
  setDeltaT(deltaT()/factor);

  bool ret=false;
  const DoubleSignal * sigfft=sig->saveType(Spectrum);
  CONST_LOCK_SAMPLES(double, sigSamples, sigfft)
    LOCK_SAMPLES(double, thisSamples, this)
      thisSamples[0]=sigSamples[0]*factor;
      for(int i=1; i<nOriginalSpectrum; i++) {
        thisSamples[i]=sigSamples[i]*factor;
        thisSamples[_nSamples-i]=sigSamples[sig->_nSamples -i]*factor;
      }
      if(_nSamples & 0x00000001) {
        thisSamples[nOriginalSpectrum]=0.0;
      } else {
        thisSamples[nOriginalSpectrum]=sigSamples[sig->_nyquistIndex]*factor;
      }
      int n=_nSamples-nOriginalSpectrum;
      for(int i=nOriginalSpectrum+1; i<=n; i++) {
        thisSamples[i]=0.0;
      }
      setType(Spectrum);
      fastFourierTransform(sig->type());
      ret=true;
    UNLOCK_SAMPLES(this);
  UNLOCK_SAMPLES(sigfft)
  sig->restoreType(sigfft);
  return ret;
}

double GeopsyCore::DoubleSignal::phase	(	const double *	samples,
		int	i
	)		const `[inline]`

Returns the phase for frequency index i. samples must be a pointer to signal samples locked with LOCK_SAMPLES or CONST_LOCK_SAMPLES.

See also:: re(), im(), complex(), amplitude2(), amplitude()

References complex(), and QGpCoreTools::Complex::phase().

Referenced by GeopsyGui::SignalLayer::drawSignal(), GeopsyGui::SignalLayer::updateGrid(), and GeopsyCore::Signal::writeAscii().

{
  return complex(samples, i).phase();
}

ComplexSignal * GeopsyCore::DoubleSignal::positiveComplexSpectrum ( bool fullSize = false ) const

Return the complex spectrum for positive frequencies. If fullSize is true, reserve a complex signal with full number of samples (both positive and negative frequency) but only positive frequencies are initialized.

(from FFTW documentation)) For those who like to think in terms of positive and negative frequencies, this means that the positive frequencies are stored in the first half of the output and the negative frequencies are stored in backwards order in the second half of the output. (The frequency -k/n is the same as the frequency (n-k)/n.)

See also:: re()

References GeopsyCore::SignalTemplate< double >::_nSamples, CONST_LOCK_SAMPLES, GeopsyCore::csig, LOCK_SAMPLES, restoreType(), saveType(), Spectrum, TRACE, and UNLOCK_SAMPLES.

Referenced by HVStationSignals::horizontal(), and morletWavelet().

{
  TRACE;
  int n=_nSamples >> 1;
  // Copying to a full complex representation of the signal excluding negative frequencies
  ComplexSignal * csig=new ComplexSignal(fullSize ? _nSamples : n+1);
  LOCK_SAMPLES(Complex, csamples, csig)
    // Assign values to the complex signal spectra from the real2complex representation
    const DoubleSignal * sig=saveType(Spectrum);
    CONST_LOCK_SAMPLES(double, sigSamples, sig)
      csamples[0]=sigSamples[ 0 ];
      for(int i=1; i<n; i++) {
        csamples[i].setRe(sigSamples[i]);
        csamples[i].setIm(sigSamples[_nSamples-i]);
      }
      if(_nSamples & 0x00000001) {
        csamples[n].setRe(sigSamples[n]);
        csamples[n].setIm(sigSamples[_nSamples-n]);
      } else {
        csamples[n]=sigSamples[ n ];
      }
    UNLOCK_SAMPLES(sig)
    restoreType(sig);
  UNLOCK_SAMPLES(csig)
  return csig;
}

double DoubleSignal::power ( double excludeBorderRatio = 0.0 ) const

Time domain computation of total power. excludeBorderRatio can for instance 0.1 to exclude 10% of signal on both ends.

References CONST_LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  ASSERT(excludeBorderRatio<1.0 && excludeBorderRatio>=0.0);
  double p=0.0;
  CONST_LOCK_SAMPLES(double, samples, this)
    int n=(int)round(_nSamples*(1.0-excludeBorderRatio));
    for(int i=(int)round(_nSamples*excludeBorderRatio); i<n; i++) {
      p+=samples[i]*samples[i];
    }
  UNLOCK_SAMPLES(this);
  return p;
}

double GeopsyCore::DoubleSignal::re	(	const double *	samples,
		int	i
	)		const `[inline]`

Retruns the real part of sample with index i. samples must be a pointer to signal samples locked with LOCK_SAMPLES or CONST_LOCK_SAMPLES.

Format of inplace r2c fft:

r0, r1, r2, ..., r(n/2), i((n+1)/2-1), ..., i2, i1

Format of r2c fft:

c0, c1, c2, ..., c(n/2-1), c(n/2), c(n/2-1)*, ..., c2*, c1*

Index i can range from 0 to n. Positive frequencies are stored from 0 to n/2, negative frequencies from n/2+1 (lowest) to n (highest frequency closest to 0 Hz).

See also:: im(), complex(), amplitude2(), amplitude(), phase()

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and TRACE.

{
  TRACE;
  if(i<=_nyquistIndex) {
    return samples[i];
  } else {
    return samples[_nSamples-i];
  }
}

void DoubleSignal::restoreType ( SignalType st )

Transform to initial type

See also:: saveType()

References TRACE.

Referenced by complexSpectrum(), DampingResults::compute(), copyAmplitudeFrom(), copySamplesFrom(), filter(), overSample(), positiveComplexSpectrum(), subtractValue(), PtMotionResults::updateSignals(), GeopsyCore::Signal::writeGse(), GeopsyCore::Signal::writeSac(), GeopsyCore::Signal::writeSeg2(), and GeopsyCore::Signal::writeSegySu().

{
  TRACE;
  fastFourierTransform(st);
}

void GeopsyCore::DoubleSignal::restoreType ( const DoubleSignal * sig ) const

double GeopsyCore::DoubleSignal::samplingFrequency ( ) const [inline]

Referenced by GeopsyCore::Signal::compare(), Process::run(), GeopsyCore::SignalHeaderObject::samplingFrequency(), timeCorrelation(), and GeopsyCore::Signal::writeMiniSeed().

{return 1.0/_deltaT;}

DoubleSignal::SignalType DoubleSignal::saveType ( SignalType st )

Transform type into st

References TRACE.

Referenced by complexSpectrum(), DampingResults::compute(), copyAmplitudeFrom(), copySamplesFrom(), filter(), overSample(), positiveComplexSpectrum(), GeopsyCore::SubSignalPool::saveGeopsySignal(), subtractValue(), PtMotionResults::updateSignals(), GeopsyCore::Signal::writeGse(), GeopsyCore::Signal::writeSac(), GeopsyCore::Signal::writeSeg2(), and GeopsyCore::Signal::writeSegySu().

{
  TRACE;
  SignalType oldSt=_type;
  fastFourierTransform(st);
  return oldSt;
}

const DoubleSignal* GeopsyCore::DoubleSignal::saveType ( SignalType st ) const

void GeopsyCore::DoubleSignal::setAmplitude	(	double *	samples,
		int	i,
		double	val
	)		const `[inline]`

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and sqrt().

{
  if(i==0) {
    samples[i]=val;
  } else if(i!=_nyquistIndex || (_nSamples & 0x00000001)) {
    double ampl=::sqrt(samples[i]*samples[i]+
                       samples[_nSamples-i]*samples[_nSamples-i]);
    if(ampl==0.0) return;
    double factor=val/ampl;
    samples[i ]*=factor;
    samples[_nSamples-i]*=factor;
  }  else { // For even sequency, sample at nyquist frequency
    samples[i]=val;
  }
}

void GeopsyCore::DoubleSignal::setComplex	(	double *	samples,
		int	i,
		const Complex &	value
	)		const `[inline]`

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, QGpCoreTools::Complex::im, QGpCoreTools::Complex::re, and TRACE.

{
  TRACE;
  if(i>0 && i<_nyquistIndex) {
    samples[i]=value.re();
    samples[_nSamples-i]=value.im();
  } else {
    samples[i]=value.re();
  }
}

void GeopsyCore::DoubleSignal::setDeltaF ( double newval ) [inline]

Referenced by HVStationSignals::horizontal().

{setDeltaT(newval);}

virtual void GeopsyCore::DoubleSignal::setDeltaT ( double newval ) [inline, virtual]

Reimplemented in GeopsyCore::Signal.

Referenced by TFAResults::compute(), DampingResults::compute(), and GeopsyCore::StationProcessSignals::setHighPassFilter().

{_deltaT=newval;}

void GeopsyCore::DoubleSignal::setIm	(	double *	samples,
		int	i,
		double	value
	)		const `[inline]`

References GeopsyCore::SignalTemplate< double >::_nSamples, and TRACE.

{
  TRACE;
  if(i<_nyquistIndex && i>0) {
    samples[_nSamples-i]=value;
  }
  // Not allowed to change the value for all other cases
}

void GeopsyCore::DoubleSignal::setNSamples ( int n ) [inline, virtual]

Reimplemented from GeopsyCore::SignalTemplate< double >.

Reimplemented in GeopsyCore::Signal, and GeopsyCore::DynamicSignal.

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and TRACE.

{
  TRACE;
  SignalTemplate<double>::setNSamples(n);
  _nyquistIndex=_nSamples >> 1;
}

void GeopsyCore::DoubleSignal::setRe	(	double *	samples,
		int	i,
		double	value
	)		const `[inline]`

References GeopsyCore::SignalTemplate< double >::_nSamples, _nyquistIndex, and TRACE.

{
  TRACE;
  if(i<=_nyquistIndex) {
    samples[i]=value;
  } else {
    samples[_nSamples-i]=value;
  }
}

void GeopsyCore::DoubleSignal::setType ( SignalType t ) [inline]

Referenced by GeopsyCore::Signal::correlation(), HVStationSignals::horizontal(), MatFormat::load(), GeopsyCore::CitySignal::loadSignals(), GeopsyCore::Signal::normalizedCorrelation(), GeopsyCore::Signal::read(), GeopsyCore::GeopsyCoreEngine::showSignal(), GeopsyCore::SubSignalPool::waveletTransform(), and GeopsyCore::Signal::xml_setProperty().

{if(t!=UndefinedSignalType) _type=t;}

bool DoubleSignal::setValue ( double val )

Set value to all samples.

References LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    for(int i=0; i<_nSamples; i++) {
      thisSamples[i]=val;
    }
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

bool DoubleSignal::shift ( double dt )

Shift the signal by dt seconds. It may be any number even not a multiple of the sampling period. Better if dt is not greater than the sampling period. If it is greater prefer a shift based on sample index. This shift is rather slow because it requires two fft transforms to and from frequency domain.

A positive dt shifts the signal towards the left (or the past).

References QGpCoreTools::cos(), QGpCoreTools::Complex::im, LOCK_SAMPLES, QGpCoreTools::Complex::re, QGpCoreTools::sin(), TRACE, and UNLOCK_SAMPLES.

Referenced by Process::run(), GeopsyCore::SubSignalPool::shift(), and T0Correlation::value().

{
  TRACE;
  bool ret=false;
  SignalType st=saveType(Spectrum);
  LOCK_SAMPLES(double, thisSamples, this)
    int n=_nSamples >> 1;
    double domega=2*M_PI/duration()*dt;
    for(int i=1; i<=n; i++) {
      Complex c(thisSamples[ i ], thisSamples[ _nSamples -i ]);
      Complex dphi(cos( i*domega), sin (i*domega) );
      c*=dphi;
      thisSamples[ i ]=c.re();
      thisSamples[ _nSamples -i ]=c.im();
    }
    ret=true;
  UNLOCK_SAMPLES(this)
  restoreType(st);
  return ret;
}

bool DoubleSignal::smooth	(	const SmoothingParameters &	param,
		double	minFreq,
		double	maxFreq
	)

smoothing with the "Konno-Ohmachi" function

(sin(exp*log10(f/fc))/(exp*log10(f/fc))^4

where fc is the frequency around which the smoothing is performed exp determines the exponent 10^(1/exp) is the half-width of the peak

References TRACE.

Referenced by HVStationSignals::horizontal(), StructureStationProcessSignals::setProcessed(), HVStationSignals::setProcessed(), and MonoStation::SpectrumStationSignals::setProcessed().

{
  TRACE;
  switch (param.type()) {
  case SmoothingParameters::NoSmoothing:
    break;
  case SmoothingParameters::KonnoOhmachi:
    return smoothKonno(param.constant(), minFreq, maxFreq);
  case SmoothingParameters::Constant:
    return smoothFixedTriang(param.constant(), minFreq, maxFreq);
  case SmoothingParameters::Proportional:
    return smoothFreqDepTriang(param.constant(), minFreq, maxFreq);
  }
  return true;
}

bool DoubleSignal::sqrt ( )

Take the square root of all samples

References LOCK_SAMPLES, QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

Referenced by amplitude(), HVStationSignals::horizontal(), and setAmplitude().

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    for(int i=0;i < _nSamples;i++ )
      thisSamples[ i ]=::sqrt(thisSamples[ i ] );
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

bool DoubleSignal::square ( )

Take the square of all samples

References LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by HVStationSignals::horizontal().

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    for(int i=0;i < _nSamples;i++ )
      thisSamples[ i ] *= thisSamples[ i ];
     ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

void DoubleSignal::stalta	(	double	tsta,
		double	tlta,
		DoubleSignal *	ratio
	)		const

Calculates the sta/lta ratio and stores it in ratio. while i<nsta (or lsta), sta (or lsta) is set to the average value of the signal.

The sta and lta are always computed over the past samples.

Averages are calculated with arithmetic sums. Use abs() or square() before calling this function.

References CONST_LOCK_SAMPLES, LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::stalta().

{
  TRACE;
  CacheProcess cp;
  cp << this << ratio;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    LOCK_SAMPLES(double, ratioSamples, ratio)
      double fSamp=1.0/_deltaT;
      int nsta=(int) (tsta * fSamp);
      int nlta=(int) (tlta * fSamp);
      if(nsta > _nSamples) nsta=_nSamples;
      if(nlta > _nSamples) nlta=_nSamples;
      double moy=average(0, _nSamples);
      double sta_sum=moy * (double) nsta, lta_sum=moy * (double) nlta;
      double r=tlta/tsta;
      int i;
      // lta_sum=sta_sum => slta=nlta/nsta;
      //printf("i\tsig\tstasum\tltasum\n");
      for(i=0;i < nsta;i++ ) {
        //printf("%i\t%lf\t%lf\t%lf\n",i,_samples[i],sta_sum,lta_sum);
        sta_sum += thisSamples[ i ] - moy;
        lta_sum += thisSamples[ i ] - moy;
        ratioSamples[ i ]=r * sta_sum/lta_sum;
      }
      for( ;i < nlta;i++ ) {
        //printf("%i\t%lf\t%lf\t%lf\n",i,_samples[i],sta_sum,lta_sum);
        sta_sum += thisSamples[ i ] - thisSamples[ i - nsta ];
        lta_sum += thisSamples[ i ] - moy;
        ratioSamples[ i ]=r * sta_sum/lta_sum;
      }
      for( ;i < _nSamples;i++ ) {
        //printf("%i\t%lf\t%lf\t%lf\n",i,_samples[i],sta_sum,lta_sum);
        sta_sum += thisSamples[ i ] - thisSamples[ i - nsta ];
        lta_sum += thisSamples[ i ] - thisSamples[ i - nlta ];
        ratioSamples[ i ]=r * sta_sum/lta_sum;
      }
    UNLOCK_SAMPLES(ratio)
  UNLOCK_SAMPLES(this)
}

void DoubleSignal::staltaToKeep	(	const WindowingParameters::ShortLongTermAverage &	param,
		KeepSignal *	keep,
		double	delay
	)		const

Calculates the sta/lta ratio and test if it is between min and max, set keep accordingly. delay is the delay of this sig relative to the beginning of keep signal. keep and this signals must have the same deltaT().

while i<nsta (or lsta), sta (or lsta) is set to the average value of the signal.

The sta and lta are always computed over the past samples.

Averages are calculated with arithmetic sums. Use abs() or square() before calling this function.

References CONST_LOCK_SAMPLES, GeopsyCore::KeepSignal::deltaT(), QGpCoreTools::endl(), LOCK_SAMPLES, GeopsyCore::WindowingParameters::ShortLongTermAverage::longLength, GeopsyCore::WindowingParameters::ShortLongTermAverage::maximumThreshold, GeopsyCore::WindowingParameters::ShortLongTermAverage::minimumThreshold, GeopsyCore::SignalTemplate< sampleType >::nSamples(), GeopsyCore::WindowingParameters::ShortLongTermAverage::shortLength, QGpCoreTools::tr(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::StationSignals::setKeep().

{
  TRACE;
  if(_deltaT<=0) {
    App::stream() << tr("DoubleSignal::staltaToKeep: null deltaT, aborting.") << endl;
    return;
  }
  if(_deltaT!=keep->deltaT()) {
    App::stream() << tr("DoubleSignal::staltaToKeep: deltaT must be equal, aborting.") << endl;
    return;
  }
  int thisStart, keepStart;
  int nSamples=commonSamples(delay, keep->nSamples(), thisStart, keepStart);
  // OK, calculate keep from keepStart to keepStart+nSamp
  CacheProcess cp;
  cp << this << keep;
  CONST_LOCK_SAMPLES(double, thisSamples, this)
    LOCK_SAMPLES(int, keepSamples, keep)
      int nsta=(int) (param.shortLength/_deltaT);
      int nlta=(int) (param.longLength/_deltaT);
      double aver=average(thisStart, nSamples);
      //printf("nsta %i nlta %i moy %lf\n",nsta, nlta, moy);
      double sta_sum=aver*(double)nsta;
      double lta_sum=aver*(double)nlta;
      double corr=param.shortLength/param.longLength;
      double ratioCorr;
      double minCorr=corr * param.minimumThreshold;
      double maxCorr=corr*param.maximumThreshold;
      int i;
      // lta_sum=sta_sum => slta=nlta/nsta;
      //printf("i\tsig\tstasum\tltasum\n");
      double val;
      int ithis;
      int n=nsta;
      if(nSamples < n)
        n=nSamples;
      for(i=0;i < n;i++ ) {
        //printf("%i\t%lf\t%lf\t%lf\n",i,_samples[i],sta_sum,lta_sum);
        val=thisSamples[i+thisStart]-aver;
        sta_sum+=val;
        lta_sum+=val;
        ratioCorr=sta_sum/lta_sum;
        if(ratioCorr < minCorr || ratioCorr > maxCorr)
          keepSamples[ i + keepStart ]=0;
      }
      n=nlta;
      if(nSamples < n)
        n=nSamples;
      for( ;i < n;i++ ) {
        //printf("%i\t%lf\t%lf\t%lf\n",i,_samples[i],sta_sum,lta_sum);
        ithis=i+thisStart;
        val=thisSamples[ithis];
        sta_sum+=val-thisSamples[ithis-nsta];
        lta_sum+=val-aver;
        ratioCorr=sta_sum/lta_sum;
        if(ratioCorr<minCorr || ratioCorr>maxCorr)
          keepSamples[i+keepStart]=0;
      }
      for(; i<nSamples; i++) {
        //printf("%i\t%lf\t%lf\t%lf\n",i,_samples[i],sta_sum,lta_sum);
        ithis=i+thisStart;
        val=thisSamples[ithis];
        sta_sum+=val-thisSamples[ithis-nsta];
        lta_sum+=val-thisSamples[ithis-nlta];
        ratioCorr=sta_sum/lta_sum;
        if(ratioCorr<minCorr || ratioCorr>maxCorr)
          keepSamples[i+keepStart]=0;
      }
    UNLOCK_SAMPLES(keep)
  UNLOCK_SAMPLES(this)
}

bool DoubleSignal::stddevClip ( double factor )

Clip signal exceeding factor times standard deviation. Re-compute standard deviation and clip again until no more clipping is required. The DC is first removed.

References QGpCoreTools::endl(), LOCK_SAMPLES, QGpCoreTools::sqrt(), QGpCoreTools::tr(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::stddevClip().

{
  TRACE;
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    SignalType st=saveType(Waveform); // make sure we are in time domain
    double sum=average(0, _nSamples)*_nSamples;
    double sum2=average2(0, _nSamples)*(_nSamples-1);
    double invNSamples=1.0/_nSamples;
    double factor2=factor * factor;
    double varianceFactor=factor2 * (sum2 - sum*invNSamples*sum)/(_nSamples-1);
    double stddevFactor=::sqrt(varianceFactor);
    double newStddevFactor;
    for(int ir=0; ir<10; ir++) {
      for(int i=0; i<_nSamples; i++) {
        if(thisSamples[i]>stddevFactor) {
          sum+=stddevFactor - thisSamples[i];
          sum2+=varianceFactor - thisSamples[i]*thisSamples[i];
          thisSamples[i]=stddevFactor;
        } else if(thisSamples[i] < -stddevFactor) {
          sum+=-stddevFactor-thisSamples[i];
          sum2+=varianceFactor-thisSamples[i]*thisSamples[i];
          thisSamples[i]=-stddevFactor;
        }
      }
      if(_nSamples<2) {
        varianceFactor=0.0;
      } else {
        varianceFactor=factor2 *(sum2-sum*invNSamples*sum)/(_nSamples-1);
        if(varianceFactor<0.0) varianceFactor=0.0;
      }
      newStddevFactor=::sqrt(varianceFactor);
      if(newStddevFactor>0.9*stddevFactor) {
        break;
      }
      stddevFactor=newStddevFactor;
    }
    if(newStddevFactor<0.9*stddevFactor) {
      App::stream() << tr("Convergence criterium not reached (less than 10% variantion of stddev).") << endl;
    }
    restoreType(st);
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

bool GeopsyCore::DoubleSignal::subtractValue ( double val = 0.0 )

Subtract value val from all samples. If val is null, the average of all samples is subtracted (DC removal).

References GeopsyCore::SignalTemplate< double >::_nSamples, GeopsyCore::SignalTemplate< double >::average(), LOCK_SAMPLES, restoreType(), saveType(), TRACE, UNLOCK_SAMPLES, and Waveform.

Referenced by LinearFKActiveStationSignals::beginPreprocess(), DampingResults::compute(), RealTimeStationSignals::setProcessed(), StructureStationProcessSignals::setProcessed(), ArrayCore::SPACStationSignals::setProcessed(), MonoStation::SpectrumStationSignals::setProcessed(), ArrayCore::FKStationSignals::setProcessed(), and GeopsyCore::SubSignalPool::subtractValue().

{
  TRACE;
  if(val==0.0) val=average(0, _nSamples);
  bool ret=false;
  if(val!= 0.0) {
    LOCK_SAMPLES(double, thisSamples, this)
      SignalType ost=saveType(Waveform);
      for(int i=0; i < _nSamples; i++ )
        thisSamples[ i ] -= val;
      restoreType(ost);
      ret=true;
    UNLOCK_SAMPLES(this)
  }
  return ret;
}

bool DoubleSignal::taper	(	const TimeRange &	t,
		const TaperParameters &	param
	)

Applies a taper window to the signal over the range t.

Returns false if not processed due to memory problem.

Reimplemented in GeopsyCore::Signal.

References GeopsyCore::TaperDelegate::apply(), LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by LinearFKActiveStationSignals::endPreprocess(), RealTimeStationSignals::setProcessed(), StructureStationProcessSignals::setProcessed(), ArrayCore::SPACStationSignals::setProcessed(), MonoStation::SpectrumStationSignals::setProcessed(), ArrayCore::FKStationSignals::setProcessed(), and LinearFKActiveStationSignals::taper().

{
  TRACE;
  bool ret=false;
  SignalType st=saveType(Waveform);
  LOCK_SAMPLES(double, thisSamples, this)
    TaperDelegate::RealSample samples(thisSamples);
    TaperDelegate window(&param);
    window.apply(samples, _nSamples, t, samplingFrequency());
    ret=true;
  UNLOCK_SAMPLES(this)
  restoreType(st);
  return ret;
}

bool DoubleSignal::taper ( const TaperParameters & param )

Applies a taper window on the complete signal.

Returns false if not processed due to memory problem.

References GeopsyCore::TaperDelegate::apply(), LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  bool ret=false;
  SignalType st=saveType(Waveform);
  LOCK_SAMPLES(double, thisSamples, this)
    TaperDelegate::RealSample samples(thisSamples);
    TaperDelegate window(&param);
    window.apply(samples, 0, _nSamples-1, 0, _nSamples-1);
    ret=true;
  UNLOCK_SAMPLES(this)
  restoreType(st);
  return ret;
}

bool DoubleSignal::timeCorrelation	(	const DoubleSignal *	s1,
		const DoubleSignal *	s2,
		double	maxDelay
	)

Calculates the correlation between signal s1 and s2 in time domain and stores the result in this signal. The correlation is computed with a maximum time delay of maxDelay seconds, assuming that signals are already shifted by dt seconds. s1(t1+dt) is at the same time as s2(t2). The number of samples for the resulting signal is 2*(maxDelay /detlaT())+1. Sampling frequencies of two signals must be equal.

The time window used for correlation is only the overlapping window decreased by the double of maxDelay (md).

          S1      --------------------------------------
                 | dt | md |     overlap           | md |
          S2           ------------------------------------------------------
                  -----------> time scale

          S1      ----------0000000000000000000000000000
                 | dt | md |     overlap           | md |
          S2           0000000000000000000000000000--------------------------

          S1      -----0000000000000000000000000000-----
                 | dt | md |     overlap           | md |
          S2           -----0000000000000000000000000000---------------------

If this signal has the correct number of samples and the correct sampling frequency, the value of this signal are not cleaned and the results are stacked on existing values. If not, the signal is initialized with 0.

This is the slow alternative not based on Fourier's transform. Prefer correlation() for all applications.

References CONST_LOCK_SAMPLES, deltaT(), LOCK_SAMPLES, GeopsyCore::SignalTemplate< sampleType >::nSamples(), samplingFrequency(), QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

{
  TRACE;
  if(!correlationCheck(s1, s2, maxDelay)) {
    return false;
  }

  double sampFreq=s1->samplingFrequency();
  int nmd=(int)round(maxDelay*sampFreq);     // number of samples for maxDelay
  int n1=s1->nSamples();
  int n2=s2->nSamples();

  // Init the result signal if necessary, else stack
  if(nSamples()!=2*nmd+1 || deltaT()!=s1->deltaT()) {
    if(isAllocated()) {
      return false;
    }
    setNSamples(2*nmd+1);
    setDeltaT(s1->deltaT());
    initValues (0.0, 0, 2*nmd+1);
  }

  bool ret=false;
  CONST_LOCK_SAMPLES(double, s1Samples, s1)
    CONST_LOCK_SAMPLES(double, s2Samples, s2)
      LOCK_SAMPLES(double, thisSamples, this)
        thisSamples+=nmd;
        double corr, en1, en2, v1, v2;
        for(int delay=-nmd; delay<=0; delay++) {
          corr=0.0;
          en1=0.0;
          en2=0.0;
          int n=n1;
          if(n>n2-delay) {
            n=n2-delay;
          }
          for(int i=-delay; i<n; i++) {
            v1=s1Samples[i];
            v2=s2Samples[i+delay];
            corr+=v1*v2;
            en1+=v1*v1;
            en2+=v2*v2;
          }
          if(en1>0.0 && en2>0.0) {
            thisSamples[delay]+=corr/::sqrt(en1*en2);
          }
        }
        for(int delay=1; delay<=nmd; delay++) {
          corr=0.0;
          en1=0.0;
          en2=0.0;
          int n=n2;
          if(n>n1+delay) {
            n=n1+delay;
          }
          for(int i=delay; i<n; i++) {
            v1=s1Samples[i-delay];
            v2=s2Samples[i];
            corr+=v1*v2;
            en1+=v1*v1;
            en2+=v2*v2;
          }
          if(en1>0.0 && en2>0.0) {
            thisSamples[delay]+=corr/::sqrt(en1*en2);
          }
        }
        ret=true;
      UNLOCK_SAMPLES(this)
    UNLOCK_SAMPLES(s2)
  UNLOCK_SAMPLES(s1)
  return ret;
}

static SignalType GeopsyCore::DoubleSignal::type ( QString t ) [static]

QString DoubleSignal::typeString ( SignalType st ) [static]

Return signal type as string.

References TRACE.

Referenced by GeopsyCore::SignalProcess::fastFourierTransform(), GeopsyCore::Signal::header(), GeopsyCore::SignalHeaderObject::type(), and GeopsyCore::Signal::xml_writeProperties().

{
  TRACE;
  switch (st) {
  case RealSpectrum: return "RealSpectrum";
  case CAmpWaveform: return "CAmpWaveform";
  case CPhaseWaveform: return "CPhaseWaveform";
  case Spectrum: return "Spectrum";
  case Phase: return "Phase";
  case ArrivalTime: return "ArrivalTime";
  case Waveform: return "Waveform";
  default: return "Undefined";
  }
}

bool DoubleSignal::unglitch ( double threshold )

References LOCK_SAMPLES, QGpCoreTools::sqrt(), TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::unglitch().

{
  TRACE;
  // TODO: not checked for spectrum amplitude (odd or even)
  // would be better to use inline functions
  bool ret=false;
  LOCK_SAMPLES(double, thisSamples, this)
    bool modif=false;
    double mul=threshold/100.;
    if(_type==Waveform) {
      double vm1, v, vp1;
      v=thisSamples[ 0 ];
      vp1=thisSamples[ 1 ];
      if(v > vp1 * mul) {
        thisSamples[ 0 ]=vp1;
        modif=true;
      }
      int i;
      for(i=1; i < _nSamples - 1; i++ ) {
        vm1=v;
        v=vp1;
        vp1=thisSamples[ i + 1 ];
        if(( v > vm1 * mul) && (v > vp1 * mul) ) {
          thisSamples[ i ]=(vm1 + vp1)/2.;
          modif=true;
        }
      }
      vm1=v;
      v=vp1;
      if(v > vm1 * mul) {
        thisSamples[ i ]=vm1;
        modif=true;
      }
    } else if(isSpectrum()) {
      double vm1, v, vp1;
      mul *= mul;
      v=thisSamples[ 0 ] * thisSamples[ 0 ];
      vp1=thisSamples[ 1 ] * thisSamples[ 1 ] + thisSamples[ _nSamples - 1 ] * thisSamples[ _nSamples - 1 ];
      if(v > vp1 * mul) {
        thisSamples[ 0 ]=::sqrt(vp1);
        modif=true;
      }
      int i;
      for(i=1; i < _nSamples/2; i++ ) {
        vm1=v;
        v=vp1;
        vp1=thisSamples[ i + 1 ] * thisSamples[ i + 1 ] + thisSamples[ _nSamples - i - 1 ] * thisSamples[ _nSamples - i - 1 ];
        if(( v > vm1 * mul) && (v > vp1 * mul) ) {
          double corr=( ::sqrt(vm1) + ::sqrt(vp1) )/2./::sqrt(v);
          thisSamples[ i ] *= corr;
          thisSamples[ _nSamples - i ] *= corr;
          modif=true;
        }
      }
      vm1=v;
      v=vp1;
      if(v > vm1 * mul) {
        thisSamples[ i ]=::sqrt(vm1);
        modif=true;
      }
    }
    ret=true;
  UNLOCK_SAMPLES(this)
  return ret;
}

void GeopsyCore::DoubleSignal::unlockProperies ( ) const [inline]

{_propertiesMutex.unlock();}

bool DoubleSignal::whiten ( )

Whiten the spectrum

References LOCK_SAMPLES, TRACE, and UNLOCK_SAMPLES.

Referenced by GeopsyCore::SubSignalPool::whiten().

{
  TRACE;
  bool ret=false;
  SignalType st=saveType(Spectrum); // make sure we are in frequency domain
  LOCK_SAMPLES(double, thisSamples, this)
    for(int i=0; i<_nyquistIndex; i++) {
      setAmplitude(thisSamples, i, 1.0);
    }
    ret=true;
  UNLOCK_SAMPLES(this)
  restoreType(st);
  return ret;
}