algorithms.spectral

Module: algorithms.spectral

Spectral transforms are used in order to estimate the frequency-domain representation of time-series. Several methods can be used and this module contains implementations of several algorithms for the calculation of spectral transforms.

Functions

nitime.algorithms.spectral.dpss_windows(N, NW, Kmax, interp_from=None, interp_kind='linear')

Returns the Discrete Prolate Spheroidal Sequences of orders [0,Kmax-1] for a given frequency-spacing multiple NW and sequence length N.

Parameters:

N : int

sequence length

NW : float, unitless

standardized half bandwidth corresponding to 2NW = BW/f0 = BW*N*dt but with dt taken as 1

Kmax : int

number of DPSS windows to return is Kmax (orders 0 through Kmax-1)

interp_from : int (optional)

The dpss can be calculated using interpolation from a set of dpss with the same NW and Kmax, but shorter N. This is the length of this shorter set of dpss windows.

interp_kind : str (optional)

This input variable is passed to scipy.interpolate.interp1d and specifies the kind of interpolation as a string (‘linear’, ‘nearest’, ‘zero’, ‘slinear’, ‘quadratic, ‘cubic’) or as an integer specifying the order of the spline interpolator to use.

Returns:

v, e : tuple,

v is an array of DPSS windows shaped (Kmax, N) e are the eigenvalues

Notes

Tridiagonal form of DPSS calculation from:

Slepian, D. Prolate spheroidal wave functions, Fourier analysis, and uncertainty V: The discrete case. Bell System Technical Journal, Volume 57 (1978), 1371430

nitime.algorithms.spectral.freq_response(b, a=1.0, n_freqs=1024, sides='onesided')

Returns the frequency response of the IIR or FIR filter described by beta and alpha coefficients.

Parameters:

b : beta sequence (moving average component)

a : alpha sequence (autoregressive component)

n_freqs : size of frequency grid

sides : {‘onesided’, ‘twosided’}

compute frequencies between [-PI,PI), or from [0, PI]

Returns:

fgrid, H(e^jw)

Notes

For a description of the linear constant-coefficient difference equation, see http://en.wikipedia.org/wiki/Z-transform

nitime.algorithms.spectral.get_spectra(time_series, method=None)

Compute the spectra of an n-tuple of time series and all of the pairwise cross-spectra.

Parameters:

time_series : float array

The time-series, where time is the last dimension

method : dict, optional

contains: this_method:’welch’

indicates that mlab.psd() will be used in order to calculate the psd/csd, in which case, additional optional inputs (and default values) are:

NFFT=64

Fs=2pi

detrend=mlab.detrend_none

window=mlab.window_hanning

n_overlap=0

this_method:’periodogram_csd’

indicates that periodogram() will be used in order to calculate the psd/csd, in which case, additional optional inputs (and default values) are:

Skx=None

Sky=None

N=None

sides=’onesided’

normalize=True

Fs=2pi

this_method:’multi_taper_csd’

indicates that multi_taper_psd() used in order to calculate psd/csd, in which case additional optional inputs (and default values) are:

BW=0.01

Fs=2pi

sides = ‘onesided’

Returns:

f : float array

The central frequencies for the frequency bands for which the spectra are estimated

fxy : float array

A semi-filled matrix with the cross-spectra of the signals. The csd of signal i and signal j is in f[j][i], but not in f[i][j] (which will be filled with zeros). For i=j fxy[i][j] is the psd of signal i.

nitime.algorithms.spectral.get_spectra_bi(x, y, method=None)

Computes the spectra of two timeseries and the cross-spectrum between them

Parameters:

x,y : float arrays

Time-series data

method : dict, optional

See get_spectra() documentation for details

Returns:

f : float array

The central frequencies for the frequency bands for which the spectra are estimated

fxx : float array

The psd of the first signal

fyy : float array

The psd of the second signal

fxy : float array

The cross-spectral density of the two signals

nitime.algorithms.spectral.mtm_cross_spectrum(tx, ty, weights, sides='twosided')

The cross-spectrum between two tapered time-series, derived from a multi-taper spectral estimation.

Parameters:

tx, ty : ndarray (K, ..., N)

The complex DFTs of the tapered sequence

weights : ndarray, or 2-tuple or list

Weights can be specified as a length-2 list of weights for spectra tx and ty respectively. Alternatively, if tx is ty and this function is computing the spectral density function of a single sequence, the weights can be given as an ndarray of weights for the spectrum. Weights may be

  • scalars, if the shape of the array is (K, ..., 1)
  • vectors, with the shape of the array being the same as tx or ty

sides : str in {‘onesided’, ‘twosided’}

For the symmetric spectra of a real sequence, optionally combine half of the frequencies and scale the duplicate frequencies in the range (0, F_nyquist).

Notes

spectral densities are always computed as

S_{xy}^{mt}(f) = \frac{\sum_k
[d_k^x(f)s_k^x(f)][d_k^y(f)(s_k^y(f))^{*}]}{[\sum_k
d_k^x(f)^2]^{\frac{1}{2}}[\sum_k d_k^y(f)^2]^{\frac{1}{2}}}

nitime.algorithms.spectral.multi_taper_csd(s, Fs=6.283185307179586, NW=None, BW=None, low_bias=True, adaptive=False, sides='default', NFFT=None)

Returns an estimate of the Cross Spectral Density (CSD) function between all (N choose 2) pairs of timeseries in s, using the multitaper method. If the NW product, or the BW and Fs in Hz are not specified by the user, a bandwidth of 4 times the fundamental frequency, corresponding to NW = 4 will be used.

Parameters:

s : ndarray

An array of sampled random processes, where the time axis is assumed to be on the last axis. If ndim > 2, the number of time series to compare will still be taken as prod(s.shape[:-1])

Fs : float, Sampling rate of the signal

NW : float

The normalized half-bandwidth of the data tapers, indicating a multiple of the fundamental frequency of the DFT (Fs/N). Common choices are n/2, for n >= 4. This parameter is unitless and more MATLAB compatible. As an alternative, set the BW parameter in Hz. See Notes on bandwidth.

BW : float

The sampling-relative bandwidth of the data tapers, in Hz.

adaptive : {True, False}

Use adaptive weighting to combine spectra

low_bias : {True, False}

Rather than use 2NW tapers, only use the tapers that have better than 90% spectral concentration within the bandwidth (still using a maximum of 2NW tapers)

sides : str (optional) [ ‘default’ | ‘onesided’ | ‘twosided’ ]

This determines which sides of the spectrum to return. For complex-valued inputs, the default is two-sided, for real-valued inputs, default is one-sided Indicates whether to return a one-sided or two-sided

Returns:

(freqs, csd_est) : ndarrays

The estimatated CSD and the frequency points vector. The CSD{i,j}(f) are returned in a square “matrix” of vectors holding Sij(f). For an input array of (M,N), the output is (M,M,N)

Notes

The bandwidth of the windowing function will determine the number tapers to use. This parameters represents trade-off between frequency resolution (lower main lobe BW for the taper) and variance reduction (higher BW and number of averaged estimates). Typically, the number of tapers is calculated as 2x the bandwidth-to-fundamental-frequency ratio, as these eigenfunctions have the best energy concentration.

nitime.algorithms.spectral.multi_taper_psd(s, Fs=6.283185307179586, NW=None, BW=None, adaptive=False, jackknife=True, low_bias=True, sides='default', NFFT=None)

Returns an estimate of the PSD function of s using the multitaper method. If the NW product, or the BW and Fs in Hz are not specified by the user, a bandwidth of 4 times the fundamental frequency, corresponding to NW = 4 will be used.

Parameters:

s : ndarray

An array of sampled random processes, where the time axis is assumed to be on the last axis

Fs : float

Sampling rate of the signal

NW : float

The normalized half-bandwidth of the data tapers, indicating a multiple of the fundamental frequency of the DFT (Fs/N). Common choices are n/2, for n >= 4. This parameter is unitless and more MATLAB compatible. As an alternative, set the BW parameter in Hz. See Notes on bandwidth.

BW : float

The sampling-relative bandwidth of the data tapers, in Hz.

adaptive : {True/False}

Use an adaptive weighting routine to combine the PSD estimates of different tapers.

jackknife : {True/False}

Use the jackknife method to make an estimate of the PSD variance at each point.

low_bias : {True/False}

Rather than use 2NW tapers, only use the tapers that have better than 90% spectral concentration within the bandwidth (still using a maximum of 2NW tapers)

sides : str (optional) [ ‘default’ | ‘onesided’ | ‘twosided’ ]

This determines which sides of the spectrum to return. For complex-valued inputs, the default is two-sided, for real-valued inputs, default is one-sided Indicates whether to return a one-sided or two-sided

Returns:

(freqs, psd_est, var_or_nu) : ndarrays

The first two arrays are the frequency points vector and the estimated PSD. The last returned array differs depending on whether the jackknife was used. It is either

  • The jackknife estimated variance of the log-psd, OR
  • The degrees of freedom in a chi2 model of how the estimated PSD is distributed about the true log-PSD (this is either 2*floor(2*NW), or calculated from adaptive weights)

Notes

The bandwidth of the windowing function will determine the number tapers to use. This parameters represents trade-off between frequency resolution (lower main lobe BW for the taper) and variance reduction (higher BW and number of averaged estimates). Typically, the number of tapers is calculated as 2x the bandwidth-to-fundamental-frequency ratio, as these eigenfunctions have the best energy concentration.

nitime.algorithms.spectral.periodogram(s, Fs=6.283185307179586, Sk=None, N=None, sides='default', normalize=True)

Takes an N-point periodogram estimate of the PSD function. The number of points N, or a precomputed FFT Sk may be provided. By default, the PSD function returned is normalized so that the integral of the PSD is equal to the mean squared amplitude (mean energy) of s (see Notes).

Parameters:

s : ndarray

Signal(s) for which to estimate the PSD, time dimension in the last axis

Fs : float (optional)

The sampling rate. Defaults to 2*pi

Sk : ndarray (optional)

Precomputed FFT of s

N : int (optional)

Indicates an N-point FFT where N != s.shape[-1]

sides : str (optional) [ ‘default’ | ‘onesided’ | ‘twosided’ ]

This determines which sides of the spectrum to return. For complex-valued inputs, the default is two-sided, for real-valued inputs, default is one-sided Indicates whether to return a one-sided or two-sided

PSD normalize : boolean (optional, default=True) Normalizes the PSD

Returns:

(f, psd) : tuple

f: The central frequencies for the frequency bands PSD estimate for each row of s

nitime.algorithms.spectral.periodogram_csd(s, Fs=6.283185307179586, Sk=None, NFFT=None, sides='default', normalize=True)

Takes an N-point periodogram estimate of all the cross spectral density functions between rows of s.

The number of points N, or a precomputed FFT Sk may be provided. By default, the CSD function returned is normalized so that the integral of the PSD is equal to the mean squared amplitude (mean energy) of s (see Notes).

s : ndarray
Signals for which to estimate the CSD, time dimension in the last axis
Fs : float (optional)
The sampling rate. Defaults to 2*pi
Sk : ndarray (optional)
Precomputed FFT of rows of s
NFFT : int (optional)
Indicates an N-point FFT where N != s.shape[-1]
sides : str (optional) [ ‘default’ | ‘onesided’ | ‘twosided’ ]
This determines which sides of the spectrum to return. For complex-valued inputs, the default is two-sided, for real-valued inputs, default is one-sided Indicates whether to return a one-sided or two-sided
normalize : boolean (optional)
Normalizes the PSD
Returns:

freqs, csd_est : ndarrays

The estimated CSD and the frequency points vector. The CSD{i,j}(f) are returned in a square “matrix” of vectors holding Sij(f). For an input array that is reshaped to (M,N), the output is (M,M,N)

nitime.algorithms.spectral.tapered_spectra(s, tapers, NFFT=None, low_bias=True)

Compute the tapered spectra of the rows of s.

Parameters:

s : ndarray, (n_arr, n_pts)

An array whose rows are timeseries.

tapers : ndarray or container

Either the precomputed DPSS tapers, or the pair of parameters (NW, K) needed to compute K tapers of length n_pts.

NFFT : int

Number of FFT bins to compute

low_bias : Boolean

If compute DPSS, automatically select tapers corresponding to > 90% energy concentration.

Returns:

t_spectra : ndarray, shaped (n_arr, K, NFFT)

The FFT of the tapered sequences in s. First dimension is squeezed out if n_arr is 1.

eigvals : ndarray

The eigenvalues are also returned if DPSS are calculated here.