mspass::algorithms::deconvolution::MTPowerSpectrumEngine Class Reference

Multittaper power spectral estimator. More...

#include <MTPowerSpectrumEngine.h>

Public Member Functions
	MTPowerSpectrumEngine ()
	MTPowerSpectrumEngine (const int winsize, const double tbp, const int ntapers, const int nfftin=-1, const double dtin=1.0)
	construct with full definition.
	MTPowerSpectrumEngine (const MTPowerSpectrumEngine &parent)
	~MTPowerSpectrumEngine ()
MTPowerSpectrumEngine &	operator= (const MTPowerSpectrumEngine &parent)
mspass::seismic::PowerSpectrum	apply (const mspass::seismic::TimeSeries &d)
std::vector< double >	apply (const std::vector< double > &d)
	Low level processing of vector of data.
double	df () const
std::vector< double >	frequencies ()
int	taper_length () const
double	time_bandwidth_product () const
int	number_tapers () const
int	fftsize () const
double	dt () const
double	set_df (double dt)
	Putter equivalent of df.
int	nf ()

Detailed Description

Multittaper power spectral estimator.

The multitaper method uses averages of spectra windowed by Slepian functions. This class can be used to compute power spectra. For efficiency the design has constructors that build the Slepian functions and cache them in a private area. We use this model because computing spectra on a large data set in parallel will usually be done with a fixed time window. The expected use is that normally the engine is created once and passed as an argument to functions using it in a map operator.

This class uses the apply model for processing. It accepts raw vector or TimeSeries data. The former assumes the sample interval is 1 while the second scales the spectrum to have units of 1/Hz.

Constructor & Destructor Documentation

MTPowerSpectrumEngine() [1/3]

mspass::algorithms::deconvolution::MTPowerSpectrumEngine::MTPowerSpectrumEngine

(

)

Default constructor. Do not use as it produces a null object that is no functional.

                                             {
  taperlen = 0;
  ntapers = 0;
  nfft = 0;
  tbp = 0.0;
  deltaf = 1.0;
  operator_dt = 1.0;
  wavetable = NULL;
  workspace = NULL;
}

MTPowerSpectrumEngine() [2/3]

mspass::algorithms::deconvolution::MTPowerSpectrumEngine::MTPowerSpectrumEngine	(	const int	winsize,
		const double	tbp,
		const int	ntapers,
		const int	nfftin = -1,
		const double	dtin = 1.0
	)

construct with full definition.

This should be the normal constructor used to create this object. It creates and caches the Slepian tapers that are used on calls the apply method.

Parameters

winsize	is the length of time windows in samples the operator will be designed to compute.
tbp	is the time bandwidth product to use for the operator.
ntapers	is the number of tapers to actually use for the operator. Note the maximum ntapers is always int(tbp*2). If ntapers is more than 2*tbp a mesage will be posted to cerr and ntapers set to tbp*2.
nfftin	is the size of the fft workspace to use for computation. When less than the winsize (the default forces this) set to 2*winsize+1.
dtin	sets the operator sample interval stored in the object and used to compute frequency bin size from fft length.

                                                                {
  taperlen = winsize;
  tbp = tbpin;
  ntapers = ntpin;
  if (nfftin < winsize) {
    nfft = nextPowerOf2(winsize);
    if (nfft == winsize)
      nfft = 2 * winsize;
  } else
    nfft = nfftin;
  operator_dt = dtin;
  this->set_df(dtin);
  int nseq = static_cast<int>(2.0 * tbp);
  if (ntapers > nseq) {
    cerr << "MTPowerSpectrumEngine (WARNING):  requested number of tapers="
         << ntpin << endl
         << "is inconsistent with requested time time bandwidth product ="
         << tbp << endl
         << "Automatically reset number tapers to max allowed=" << nseq << endl;
    ntapers = nseq;
  }
  int seql(0);
  int sequ = ntapers - 1;
  double *work = new double[ntapers * taperlen];
  dpss_calc(taperlen, tbp, seql, sequ, work);
  tapers = dmatrix(ntapers, taperlen);
  int i, ii, j;
  for (i = 0, ii = 0; i < ntapers; ++i) {
    for (j = 0; j < taperlen; ++j) {
      tapers(i, j) = work[ii];
      ++ii;
    }
  }
  delete[] work;
  /* To be consistent with Prieto we use this algorithm to convert to
  what he calls the "positive standard".   That means we assure the
  center point is positive.
  */
  for (i = 0; i < ntapers; ++i) {
    int lh; // matches Prieto algorithm name - see multitaper module
    if (taperlen % 2)
      lh = static_cast<int>((taperlen + 1) / 2);
    else
      lh = static_cast<int>(taperlen / 2);
    if (tapers(i, lh) < 0.0) {
      for (j = 0; j < taperlen; ++j)
        tapers(i, j) = -tapers(i, j);
    }
  }
  wavetable = gsl_fft_complex_wavetable_alloc(nfft);
  workspace = gsl_fft_complex_workspace_alloc(nfft);
}

References set_df().

MTPowerSpectrumEngine() [3/3]

mspass::algorithms::deconvolution::MTPowerSpectrumEngine::MTPowerSpectrumEngine

(

const MTPowerSpectrumEngine &

parent

)

Standard copy constructor

    : tapers(parent.tapers) {
  taperlen = parent.taperlen;
  ntapers = parent.ntapers;
  nfft = parent.nfft;
  tbp = parent.tbp;
  operator_dt = parent.operator_dt;
  deltaf = parent.deltaf;
  wavetable = gsl_fft_complex_wavetable_alloc(nfft);
  workspace = gsl_fft_complex_workspace_alloc(nfft);
}

~MTPowerSpectrumEngine()

mspass::algorithms::deconvolution::MTPowerSpectrumEngine::~MTPowerSpectrumEngine

(

)

Destructor. Not trivial as it has to delete the fft workspace and cached tapers.

                                              {
  if (wavetable != NULL)
    gsl_fft_complex_wavetable_free(wavetable);
  if (workspace != NULL)
    gsl_fft_complex_workspace_free(workspace);
}

Member Function Documentation

apply() [1/2]

PowerSpectrum mspass::algorithms::deconvolution::MTPowerSpectrumEngine::apply

(

const mspass::seismic::TimeSeries &

d

)

\process a TimeSeries.

This is one of two methods for applying the multiaper algorithm to data. This one uses dt and data length to set the Rayleigh bin size (df). If the input data vector length is not the same as the operator length an elog complaint is posted to parent. Short data are processed but should be considered suspect unless the sizes differ by only a tiny fraction (e.g. and off by one error from rounding). Long data will be truncated on the right (i.e. sample 0 will be the start of the window used). The data return will be scaled to psd in units if 1/Hz.

Parameters

parent

is the data to process

Returns

vector containing estimated power spwecrum

                                                              {
  try {
    int k;
    /* Used to test for operator sample interval against data sample interval.
    We don't use a epsilon comparison as slippery clock data sometime shave
    sample rates small percentage difference from nominal.*/
    const double DT_FRACTION_TOLERANCE(0.001);
    const string algorithm("MTPowerSpectrumEngine");
    /* We need to define this here to allow posting problems to elog.*/
    PowerSpectrum result;
    int dsize = d.npts();
    vector<double> work;
    work.reserve(this->nfft);
    double dtfrac = fabs(d.dt() - this->operator_dt) / this->operator_dt;
    if (dtfrac > DT_FRACTION_TOLERANCE) {
      stringstream ss;
      ss << "Date sample interval=" << d.dt()
         << " does not match operator sample interval=" << this->operator_dt
         << endl
         << "Cannot proceed.  Returning a null result";
      result.elog.log_error("MTPowerSpectrumEngine::apply", ss.str(),
                            ErrorSeverity::Invalid);
      result.kill();
      return result;
    }
 
    if (dsize < taperlen) {
      stringstream ss;
      ss << "Received data window of length=" << d.npts() << " samples" << endl
         << "Operator length=" << taperlen << endl
         << "Results may be unreliable" << endl;
      result.elog.log_error(algorithm, string(ss.str()),
                            ErrorSeverity::Suspect);
      for (k = 0; k < taperlen; ++k)
        work.push_back(0.0);
      for (k = 0; k < dsize; ++k)
        work[k] = d.s[k];
    } else {
      if (dsize > taperlen) {
        stringstream ss;
        ss << "Received data window of length=" << d.npts() << " samples"
           << endl
           << "Operator length=" << taperlen << endl
           << "Results may be unreliable because data will be truncated to "
              "taper length"
           << endl;
        result.elog.log_error(algorithm, ss.str(), ErrorSeverity::Suspect);
      }
      for (k = 0; k < taperlen; ++k)
        work.push_back(d.s[k]);
    }
    /* intentionally omit try catch here because the above logic assures Sizes
    must match here. This overloaded method will throw an exception in that
    case. Note in this implementation the result returned by apply is scaled to
    assumed properly scaled to power spectrum and normalized for multitapers.*/
    vector<double> spec(this->apply(work));
 
    result = PowerSpectrum(dynamic_cast<const Metadata &>(d), spec, deltaf,
                           string("Multitaper"), 0.0, d.dt(), d.npts());
    /* We post these to metadata for the generic PowerSpectrum object. */
    result.put<double>("time_bandwidth_product", tbp);
    result.put<long>("number_tapers", ntapers);
    return result;
  } catch (...) {
    throw;
  };
}

References apply(), mspass::seismic::BasicTimeSeries::dt(), mspass::seismic::PowerSpectrum::elog, mspass::seismic::BasicSpectrum::kill(), mspass::utility::ErrorLogger::log_error(), mspass::seismic::BasicTimeSeries::npts(), and mspass::seismic::CoreTimeSeries::s.

apply() [2/2]

std::vector< double > mspass::algorithms::deconvolution::MTPowerSpectrumEngine::apply

(

const std::vector< double > &

d

)

Low level processing of vector of data.

This is lower level function that processes a raw vector of data. Since it does not know the sample interval it cannot compute the rayleigh bin size so if callers need that feature they must do that (simple) calculation themselves. Unlike the TimeSeries method this one will throw an exception if the input data size does not match the operator size. It returns power spectral density assuming a sample rate of 1. i.e. it scales to correct for the gsl fft scaling by of the forward transform by N.

Parameters

d	is the vector of data to process. d.size() must this->taperlen() value.

Returns

vector containing estimated power spectrum (usual convention with 0 containing 0 frequency value)

Exceptions

throw

a MsPASSError if the size of d does not match operator length

df()

double mspass::algorithms::deconvolution::MTPowerSpectrumEngine::df ( ) const

inline

Return the frquency bin size defined for this operator.

{ return deltaf; };

dt()

double mspass::algorithms::deconvolution::MTPowerSpectrumEngine::dt ( ) const

inline

Retrieve the internally cached required data sample interval.

{ return operator_dt; };

fftsize()

int mspass::algorithms::deconvolution::MTPowerSpectrumEngine::fftsize ( ) const

inline

Return size of fft used by this operator - usually not the same as taper length.

{ return nfft; };

frequencies()

vector< double > mspass::algorithms::deconvolution::MTPowerSpectrumEngine::frequencies

(

)

Return and std::vector of all frequencies for spectral estimates this operator computes.

                                                  {
  vector<double> f;
  /* If taperlen is odd this still works according to gsl documentation.*/
  for (int i = 0; i < this->nf(); ++i) {
    /* Here we assume i=0 frequency is 0 */
    f.push_back(deltaf * ((double)i));
  }
  return f;
}

References nf().

nf()

int mspass::algorithms::deconvolution::MTPowerSpectrumEngine::nf ( )

inline

Return tne number of frequency bins in estimates the operator will compute.

           {
    /* this simple formula depends upon integer truncation when used with
    nfft as an odd number.   For reference, this is what prieto uses in
    the python multitaper package:
    if (nfft%2 == 0):
        nf = int(nfft/2 + 1)
    else:
        nf = int((nfft+1)/2)
    they will yield the same result but this is simpler and faster
    */
    return (this->nfft) / 2 + 1;
  };

number_tapers()

int mspass::algorithms::deconvolution::MTPowerSpectrumEngine::number_tapers ( ) const

inline

Return number of tapers used by this engine.

{ return ntapers; };

operator=()

MTPowerSpectrumEngine & mspass::algorithms::deconvolution::MTPowerSpectrumEngine::operator=

(

const MTPowerSpectrumEngine &

parent

)

Standard assignment operator.

                                                                    {
  if (&parent != this) {
    taperlen = parent.taperlen;
    ntapers = parent.ntapers;
    nfft = parent.nfft;
    tbp = parent.tbp;
    operator_dt = parent.operator_dt;
    deltaf = parent.deltaf;
    tapers = parent.tapers;
    wavetable = gsl_fft_complex_wavetable_alloc(nfft);
    workspace = gsl_fft_complex_workspace_alloc(nfft);
  }
  return *this;
}

set_df()

double mspass::algorithms::deconvolution::MTPowerSpectrumEngine::set_df ( double  dt )

inline

Putter equivalent of df.

The computation of the Rayleigh bin size is complicated a bit by the folding properties of fft algorithms that have to handle odd and even length inputs differently. This algorithm uses the internally set nfft value to set the frequency bin size for even or odd nfft and the input sample interval. NOTE POSSIBLE CONFUSION that input is time sample interval NOT the actual frquency bin size. The reason is that the odd/even issue makes df dependent on if the fft size is even or odd. We include this method as a convenience as that is an implementation detail for the fft algorithm.

Note also this method sets not just df but the internally stored sample interval (symbol operator_dt in the source code.)

Parameters

dt	is the data sample interval (time domain)

Returns

computed df

                           {
    this->operator_dt = dt;
    int this_nf = this->nf();
    double fny = 1.0 / (2.0 * dt);
    this->deltaf = fny / static_cast<double>(this_nf - 1);
    return deltaf;
  };

References dt(), and nf().

taper_length()

int mspass::algorithms::deconvolution::MTPowerSpectrumEngine::taper_length ( ) const

inline

Retrieve the taper length.

{ return taperlen; };

time_bandwidth_product()

double mspass::algorithms::deconvolution::MTPowerSpectrumEngine::time_bandwidth_product ( ) const

inline

Retrieve time-bandwidth product.

{ return tbp; };

---

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

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