mtspec module

Module with routines for univariate multitaper spectrum estimation (1D). Contains the main MTSpec and MTSine classes where the estimates are made and stored.

See module mtcross for bi-variate problems

Classes

  • MTSine - A class to represent Sine Multitaper estimates

  • MTSpec - A class to represent Thomson’s multitaper estimates

Functions

  • spectrogram - Computes a spectrogram with consecutive multitaper estimates.


class mtspec.MTSine(x, ntap=0, ntimes=0, fact=1.0, dt=1.0)

Bases: object

A class for sine multitaper estimates

Attributes

Parameters
nptsint

number of points of time series

nfftint

number of points of FFT. nfft = 2*npts

Time series
xndarray [npts]

time series

xvarfloat

variance of time series

dtfloat

sampling interval

Frequency vector
nfint

number of unique frequency points of spectral estimate, assuming real time series

freqndarray [nfft]

frequency vector in Hz

dffloat

frequncy sampling interval

Method
ntapint

fixed number of tapers if ntap<0, use kopt

koptndarray [nfft,1]

number of tapers at each frequency

ntimesint

number of max iterations to perform

irealint

0 - real time series 1 - complex time series

Spectral estimates
specndarray [nfft,1]

multitaper estimate

errndarray [nfft,2]

1-std confidence interval of spectral estimate simple dof estimate

Notes

The class is in charge of estimating the adaptive sine multitaper as in Riedel and Sidorenko (1995). This is done by performing a MSE adaptive estimation. First a pilot spectral estimate is used, and S” is estimated, in order to get te number of tapers to use, using (13) of R & S for a min square error spectrum.

Unlike the prolate spheroidal multitapers, the sine multitaper adaptive process introduces a variable resolution and error in the frequency domain. Complete error information is contained in the output variables file as the corridor of 1-standard-deviation errors, and in K, the number of tapers used at each frequency. The errors are estimated in the simplest way, from the number of degrees of freedom (two per taper), not by jack-knifing. The frequency resolution is found from K*fN/Nf where fN is the Nyquist frequency and Nf is the number of frequencies estimated. The adaptive process used is as follows. A quadratic fit to the log PSD within an adaptively determined frequency band is used to find an estimate of the local second derivative of the spectrum. This is used in an equation like R & S (13) for the MSE taper number, with the difference that a parabolic weighting is applied with increasing taper order. Because the FFTs of the tapered series can be found by resampling the FFT of the original time series (doubled in length and padded with zeros) only one FFT is required per series, no matter how many tapers are used. This makes the program fast. Compared with the Thomson multitaper programs, this code is not only fast but simple and short. The spectra associated with the sine tapers are weighted before averaging with a parabolically varying weight. The expression for the optimal number of tapers given by R & S must be modified since it gives an unbounded result near points where S” vanishes, which happens at many points in most spectra. This program restricts the rate of growth of the number of tapers so that a neighboring covering interval estimate is never completely contained in the next such interval.

This method SHOULD not be used for sharp cutoffs or deep valleys, or small sample sizes. Instead use Thomson multitaper in mtspec in this same library.

References

Riedel and Sidorenko, IEEE Tr. Sig. Pr, 43, 188, 1995

Based on Bob Parker psd.f codes. Most of the comments come from his documentation as well.

Modified

September 22 2005

Calls

utils.quick, utils.adapt


class mtspec.MTSpec(x, nw=4, kspec=0, dt=1.0, nfft=0, iadapt=0, vn=None, lamb=None)

Bases: object

A class for univariate Thomson multitaper estimates

Attributes

Parameters
nptsint

number of points of time series

nfftint

number of points of FFT. Dafault adds padding.

nwflaot

time-bandwidth product

kspecint

number of tapers to use

Time series
xndarray [npts]

time series

xvarfloat

variance of time series

dtfloat

sampling interval

Frequency vector
nfint

number of unique frequency points of spectral estimate, assuming real time series

freqndarray [nfft]

frequency vector in Hz

dffloat

frequncy sampling interval

Method
iadaptint

defines methos to use 0 - adaptive multitaper 1 - unweighted, wt =1 for all tapers 2 - wt by the eigenvalue of DPSS

irealint

0 - real time series 1 - complex time series

DPSS tapers and eigenvalues
vnndarray [npts,kspec]

Slepian sequences

lambndarray [kspec]

Eigenvalues of Slepian sequences

Spectral estimates
ykcomplex ndarray [nfft,kspec]

eigencoefficients, fft of tapered series

skndarray [nfft,kspec]

eigenspectra, power spectra for each yk

specndarray [nfft,1]

multitaper estimate

sendarray [nfft,1]

degrees of freedom of estimate

wtndarray [nfft,kspec]

weights for each eigencoefficient at each frequency

Methods

  • rspec : returns the positive frequency of the spectra only

  • reshape : reshape yk’s based on F-test of line components

  • jackspec : estimate 95% confidence interval of multitaper estimate

  • qiinv : the quadratic inverse spectral estimate

  • ftest : F-test of line components in the spectra

  • df_spec : dual-frequency autospectra

References

Based on David J. Thomson’s codes, Alan Chave and Thomson’s Codes and partial codes from EISPACK, Robert L. Parker and Glenn Ierley from Scripps Institution of Oceanography. And my own Fortran90 library.

Notes

The class is in charge of estimating the adaptive weigthed multitaper spectrum, as in Thomson 1982. This is done by estimating the dpss (discrete prolate spheroidal sequences), multiplying each of the kspec tapers with the data series, take the fft, and using the adaptive scheme for a better estimation.

As a by product of the spectrum (spec), all intermediate steps are retained, which can be used for bi-variate analysis, deconvolotuion, returning to the time domain, etc. By-products include the complex information in yk, the eigenspectra sk, the jackknife 95% confidence intervals (spec_ci), the degrees of freedom (se) and the weigths wt(nf,kspec) used. See below for a complete list.

Modified

January 2022 (German Prieto)


df_spec()

Performs the dual-frequency spectrum of a signal with itself.

Returns

df_specndarray complex, 2D (nf,nf)

the complex dual-frequency cross-spectrum. Not normalized

df_cohendarray, 2D (nf,nf)

MSC, dual-freq coherence matrix. Normalized (0.0,1.0)

df_phasendarray, 2D (nf,nf)

the dual-frequency phase

Calls

utils.df_spec


ftest()

Performs the F test for a line component

Computes F-test for single spectral line components at the frequency bins given in the MTSpec class.

Returns

F : ndarray [nfft] vector of f test values, real p : ndarray [nfft] vector with probability of line component

Calls

utils.f_test


jackspec()

code to calculate adaptively weighted jackknifed 95% confidence limits

Returns

spec_cindarray [nfft,2]

real array of jackknife error estimates, with 5 and 95% confidence intervals of the spectrum.

Calls

utils.jackspec

qiinv()

Function to calculate the Quadratic Spectrum using the method developed by Prieto et al. (2007).

The first 2 derivatives of the spectrum are estimated and the bias associated with curvature (2nd derivative) is reduced.

Calculate the Stationary Inverse Theory Spectrum. Basically, compute the spectrum inside the innerband.

This approach is very similar to D.J. Thomson (1990).

Returns

qispecndarray [nfft,0]

the QI spectrum estimate

dsndarray [nfft,0]

the estimate of the first derivative

ddsndarray [nfft,0]

the estimate of the second derivative

References

G. A. Prieto, R. L. Parker, D. J. Thomson, F. L. Vernon, and R. L. Graham (2007), Reducing the bias of multitaper spectrum estimates, Geophys. J. Int., 171, 1269-1281. doi: 10.1111/j.1365-246X.2007.03592.x.

Calls

utils.qiinv


reshape(fcrit=0.95, p=None)

Reshape eigenft’s (yk) around significant spectral lines The “significant” means above fcritical probability (0.95) If probability is large at neighbouring frequencies, I will only remove the largest probability energy.

Returns recalculated yk, sk, spec, wt, and se

Parameters

fcritfloat optional

Probability value over which to reshape, default = 0.95

pndarray optional [nfft]

F-test probabilities to find fcritical If None, it will be calculated

Returns

respecndarray [nfft]

The reshaped PSD estimate

specndarray [nfft]

the PSD without the line components

ykndarray [nfft,kspec]

the eigenft’s without line components

slinendarray [nfft]

the PSD of the line components only

Calls

utils.yk_reshape


rspec(*args)

Returns the spetra at positive frequencies, checking that a real input signal was used.

Parameters

argsndarray

another array to return the positive frequencies. Could be qispec, spec_ci, etc.


mtspec.spectrogram(data, dt, twin, olap=0.5, nw=3.5, kspec=5, fmin=0.0, fmax=-1.0, iadapt=0)

Computes a spectrogram with consecutive multitaper estimates. Returns both Thomson’s multitaper and the Quadratic multitaper estimate

Parameters

dataarray_like (npts,)

Time series or sequence

dtfloat

Sampling interval in seconds of the time series.

twinfloat

Time duration in seconds of each segment for a single multitaper estimate.

olapfloat, optional

Overlap requested for the segment in building the spectrogram. Defaults = 0.5, values must be (0.0 - 0.99). Overlap rounds to the nearest integer point.

nwfloat, optional

Time-bandwidth product for Thomson’s multitaper algorithm. Default = 3.5

kspecint, optional

Number of tapers for avearaging the multitaper estimate. Default = 5

fminfloat, optional

Minimum frequency to estimate the spectrogram, otherwise returns the entire spectrogram matrix. Default = 0.0 Hz

fmaxfloat, optional

Maximum frequency to estimate the spectrogram, otherwise returns the entire spectrogram matrix. Default = 0.5/dt Hz (Nyquist frequency)

iadaptinteger, optional

User defined, determines which method for multitaper averaging to use. Default = 0 0 - Adaptive multitaper 1 - Eigenvalue weights 2 - Constant weighting

Returns

fndarray

Array of sample frequencies.

tndarray

Array of segment times.

Quadndarray

Spectrogram of x using the quadratic multitaper estimate.

MTndarray

Spectrogram of x using Thomson’s multitaper estimate.

By default, the last axis of Quad/MT corresponds to the segment times.

See Also

MTSpec: Multitaper estimate of a time series.

Notes

The code assumes a real input signals and thus mainly returns the positive frequencies. For a complex input signals, code qould require adaptation.

References

Prieto, G.A. (2022). The multitaper spectrum analysis package in Python. Seism. Res. Lett In review.

Examples

To do