mtspec.sine_psd

mtspec.multitaper.sine_psd(data, delta, number_of_tapers=None, number_of_iterations=2, degree_of_smoothing=1.0, statistics=False, verbose=False)

Wrapper method for the sine_psd subroutine in the library by German A. Prieto.

The subroutine is in charge of estimating the adaptive sine multitaper as in Riedel and Sidorenko (1995). It outputs the power spectral density (PSD).

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 Riedel and Sidorenko 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 as the corridor of 1-standard-deviation errors, and in 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.

Parameters:

• data – numpy.ndarray Array with the data.
• delta – float Sample spacing of the data.
• number_of_tapers – integer/None, optional Number of tapers to use. If none is given, the library will perform an adaptive taper estimation with a varying number of tapers for each frequency. Defaults to None.
• number_of_iterations – integer, optional Number of iterations to perform. Values less than 2 will be set to 2. Defaults to 2.
• degree_of_smoothing – float, optional Degree of smoothing. Defaults to 1.0.
• statistics – bool, optional Calculates and returns statistics. See the notes in the docstring for further details.
• verbose – bool, optional Passed to the fortran library. Defaults to False.

Returns:

Returns a list with numpy.ndarray. See the note below for details.

Note

This method will at return at least two arrays: The calculated spectrum and the corresponding frequencies. If statistics is True is will also return (in the given order) (multidimensional) arrays containing the 1-std errors (a simple dof estimate) and the number of tapers used for each frequency point.