algorithms

Module: algorithms

This module contains implementations of algorithms for time series analysis. These algorithms include:

1. Spectral estimation: calculate the spectra of time-series and cross-spectra between time-series.

get_spectra(), get_spectra_bi(), periodogram(), periodogram_csd(), DPSS_windows(), multi_taper_psd(), multi_taper_csd(), mtm_cross_spectrum()

2. Coherency: calculate the pairwise correlation between time-series in the frequency domain and related quantities.

coherency(), coherence(), coherence_regularized(), coherency_regularized(), coherency_bavg(), coherence_bavg(), coherence_partial(), coherence_partial_bavg(), coherency_phase_spectrum(), coherency_phase_delay(), coherency_phase_delay_bavg(), correlation_spectrum()

3. Event-related analysis: calculate the correlation between time-series and external events.

freq_domain_xcorr(), freq_domain_xcorr_zscored(), fir()

4. Cached coherency: A set of special functions for quickly calculating coherency in large data-sets, where the calculation is done over only a subset of the adjacency matrix edges and intermediate calculations are cached, in order to save calculation time.

cache_fft(), cache_to_psd(), cache_to_phase(), cache_to_relative_phase(), cache_to_coherency().

5. Wavelet transforms: Calculate wavelet transforms of time-series data.

wmorlet(), wfmorlet_fft(), wlogmorlet(), wlogmorlet_fft()

6. Filtering: Filter a signal in the frequency domain.

boxcar_filter()

The algorithms in this library are the functional form of the algorithms, which accept as inputs numpy array and produce numpy array outputs. Therfore, they can be used on any type of data which can be represented in numpy arrays. See also nitime.analysis for simplified analysis interfaces, using the data containers implemented in nitime.timeseries

Functions

nitime.algorithms.DPSS_windows(N, NW, Kmax)

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

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.LD_AR_est(s, order, Nfreqs, sxx=None, sides='onesided', system=False)

Finds the parameters for an autoregressive model of order norder of the process s. Using these parameters, an estimate of the PSD is calculated from [-PI,PI) in Nfreqs, or [0,PI] in {N/2+1}freqs. Uses the Levinson-Durbin recursion method, and a baised estimate of sxx (unless sxx is provided).

The model for the autoregressive process takes this convention: s[n] = a1*s[n-1] + a2*s[n-2] + ... aP*s[n-P] + v[n]

where v[n] is a zero-mean white noise process with variance=sigma_v

Parameters :

s : ndarray

The sampled autoregressive random process

order : int

The order P of the AR system

Nfreqs : int

The number of spacings on the frequency grid from [-PI,PI). If sides==’onesided’, Nfreqs/2+1 frequencies are computed from [0,PI]

sxx : ndarray (optional)

An optional, possibly unbiased estimate of the autocovariance of s

sides : str (optional)

Indicates whether to return a one-sided or two-sided PSD

system : bool (optional)

If True, return the AR system parameters, sigma_v and a{k}

Returns :

(w, ar_psd) :

w : Array of normalized frequences from [-.5, .5) or [0,.5]

ar_psd : A PSD estimate computed by sigma_v / |1-a(f)|**2 , where

a(f) = DTFT(ak)

nitime.algorithms.boxcar_filter(time_series, lb=0, ub=1, n_iterations=2)

Filters data into a frequency range.

For each of the two bounds, a low-passed version is created by convolving with a box-car and then the low-passed version for the upper bound is added to the low-passed version for the lower bound subtracted from the signal, resulting in a band-passed version

Parameters :

time_series: float array :

the signal

ub : float, optional

The cut-off frequency for the low-pass filtering as a proportion of the sampling rate. Default to 1

lb : float, optional

The cut-off frequency for the high-pass filtering as a proportion of the sampling rate. Default to 0

n_iterations: int, optional :

how many rounds of smoothing to do. Default to 2.

Returns :

float array: :

The signal, filtered

nitime.algorithms.cache_fft(time_series, ij, lb=0, ub=None, method=None, prefer_speed_over_memory=False, scale_by_freq=True)

compute and cache the windowed FFTs of the time_series, in such a way that computing the psd and csd of any combination of them can be done quickly.

Parameters :

time_series : float array

An ndarray with time-series, where time is the last dimension

ij: list of tuples :

Each tuple in this variable should contain a pair of indices of the form (i,j). The resulting cache will contain the fft of time-series in the rows indexed by the unique elements of the union of i and j

lb,ub: float :

Define a frequency band of interest, for which the fft will be cached

method: dict, optional :

See get_spectra() for details on how this is used. For this set of functions, ‘this_method’ has to be ‘welch’

Returns :

freqs, cache :

where: cache =

{‘FFT_slices’:FFT_slices,’FFT_conj_slices’:FFT_conj_slices, ‘norm_val’:norm_val}

Notes

• For these functions, only the Welch windowed periodogram (‘welch’) is available.
• Detrending the input is not an option here, in order to save time on an empty function call.

nitime.algorithms.cache_to_coherency(cache, ij)

From a set of cached spectra, calculate the coherency relationships

Parameters :

cache: dict :

the return value from cache_fft()

ij: list :

a list of (i,j) tuples, the pairs of indices for which the cross-coherency is to be calculated

Returns :

Cxy: dict :

coherence values between the time-series ij. Indexing into this dict takes the form Cxy[i,j] in order to extract the coherency between time-series i and time-series j in the original input to cache_fft()

nitime.algorithms.cache_to_phase(cache, ij)

From a set of cached set of windowed fft’s, calculate the frequency-band dependent phase for each of the channels in ij. Note that this returns the absolute phases of the time-series, not the relative phases between them. In order to get relative phases, use cache_to_relative_phase

Parameters :

cache : dict

The return value of cache_fft()

ij: list :

A list of tuples of the form (i,j) for all the indices for which to calculate the phases

Returns :

Phase : dict

The individual phases, keys are all the i and j in ij, such that Phase[i] gives you the phase for the time-series i in the input to cache_fft()

nitime.algorithms.cache_to_psd(cache, ij)

From a set of cached windowed fft, calculate the psd

Parameters :

cache : dict

Return value from cache_fft()

ij : list

A list of tuples of the form (i,j).

Returns :

Pxx : dict

The phases for the intersection of (time_series[i],time_series[j]). The keys are the intersection of i,j values in the parameter ij

nitime.algorithms.cache_to_relative_phase(cache, ij)

From a set of cached set of windowed fft’s, calculate the frequency-band dependent relative phase for the combinations ij.

Parameters :

cache: dict :

The return value from cache_fft()

ij: list :

A list of tuples of the form (i,j), all the pairs of indices for which to calculate the relative phases

Returns :

Phi_xy : dict

The relative phases between the time-series i and j. Such that Phi_xy[i,j] is the phase from time_series[i] to time_series[j].

nitime.algorithms.coherence(time_series, csd_method=None)

Compute the coherence between the spectra of an n-tuple of time_series.

Parameters of this function are in the time domain.

Parameters :

time_series: float array :

an array of different time series with time as the last dimension

csd_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

c : float array

This is a symmetric matrix with the coherencys of the signals. The coherency of signal i and signal j is in f[i][j].

See also

coherence_calculate()

Notes

This is an implementation of equation (2) of [Sun2005]:

Coh_{xy}(\lambda) = |{R_{xy}(\lambda)}|^2 = \frac{|{f_{xy}(\lambda)}|^2}{f_{xx}(\lambda) \cdot f_{yy}(\lambda)}

[Sun2005]

F.T. Sun and L.M. Miller and M. D’Esposito(2005). Measuring temporal dynamics of functional networks using phase spectrum of fMRI data. Neuroimage, 28: 227-37.

nitime.algorithms.coherence_bavg(time_series, lb=0, ub=None, csd_method=None)

Compute the band-averaged coherence between the spectra of two time series.

Input to this function is in the time domain.

Parameters :

time_series : float array

An array of time series, time as the last dimension.

lb, ub: float, optional :

The upper and lower bound on the frequency band to be used in averaging defaults to 1,max(f)

csd_method: dict, optional. :

See get_spectra() documentation for details

Returns :

c : float

This is an upper-diagonal array, where c[i][j] is the band-averaged coherency between time_series[i] and time_series[j]

nitime.algorithms.coherence_bavg_calculate(fxy, fxx, fyy)

Compute the band-averaged coherency between the spectra of two time series. input to this function is in the frequency domain

Parameters :

fxy : float array

The cross-spectrum of the time series

fyy,fxx : float array

The spectra of the signals

Returns :

float : :

the band-averaged coherence

nitime.algorithms.coherence_calculate(fxy, fxx, fyy)

Compute the coherence between the spectra of two time series.

Parameters of this function are in the frequency domain.

Parameters :

fxy : array

The cross-spectrum of the time series

fyy,fxx : array

The spectra of the signals

Returns :

float : a frequency-band-dependent measure of the linear association

between the two time series

See also

coherence()

nitime.algorithms.coherence_partial(time_series, r, csd_method=None)

Compute the band-specific partial coherence between the spectra of two time series.

The partial coherence is the part of the coherence between x and y, which cannot be attributed to a common cause, r.

Input to this function is in the time domain.

Parameters :

time_series: float array :

An array of time-series, with time as the last dimension.

r: float array :

This array represents the temporal sequence of the common cause to be partialed out, sampled at the same rate as time_series

csd_method: dict, optional :

See get_spectra() documentation for details

Returns :

f: array, :

The mid-frequencies of the frequency bands in the spectral decomposition

c: float array :

The frequency dependent partial coherence between time_series i and time_series j in c[i][j] and in c[j][i], with r partialed out

Notes

This is an implementation of equation (2) of [Sun2004]:

Coh_{xy|r} = \frac{|{R_{xy}(\lambda) - R_{xr}(\lambda) R_{ry}(\lambda)}|^2}{(1-|{R_{xr}}|^2)(1-|{R_{ry}}|^2)}

[Sun2004]

F.T. Sun and L.M. Miller and M. D’Esposito(2004). Measuring

interregional functional connectivity using coherence and partial coherence analyses of fMRI data Neuroimage, 21: 647-58.

nitime.algorithms.coherence_partial_bavg(time_series, r, csd_method=None, lb=0, ub=None)

Compute the band-averaged partial coherence between the spectra of two time series. See coherence_partial().

Input to this function is in the time domain.

Parameters :

time_series : float array

Time series data with the time on the last dimension

r : float array

Cause to be partialed out

csd_method: dict, optional :

See get_spectra() for details

lb: float, optional :

The lower bound frequency (in Hz) of the range over which the average is calculated. Default: 0

ub: float, optional :

The upper bound frequency (in Hz) of the range over which the average is calculated. Defaults to the Nyquist frequency

Returns :

c: float :

the band-averaged coherency

nitime.algorithms.coherence_partial_bavg_calculate(f, fxy, fxx, fyy, fxr, fry, frr)

Compute the band-averaged partial coherence between the spectra of two time series.

Input to this function is in the frequency domain.

Parameters :

f: the frequencies :

fxy : float array

The cross-spectrum of the time series

fyy,fxx : float array

The spectra of the signals

fxr,fry : float array

The cross-spectra of the signals with the event

Returns :

float :

the band-averaged coherency

See also

coherency, coherence, coherence_partial, coherency_bavg

nitime.algorithms.coherence_partial_calculate(fxy, fxx, fyy, fxr, fry, frr)

Compute the band-specific partial coherence between the spectra of two time series. See partial_coherence().

Input to this function is in the frequency domain.

Parameters :

fxy : float array

The cross-spectrum of the time series

fyy,fxx : float array

The spectra of the signals

fxr,fry : float array

The cross-spectra of the signals with the event

Returns :

float :

the band-averaged coherency

nitime.algorithms.coherence_regularized(time_series, epsilon, alpha, csd_method=None)

Same as coherence, except regularized in order to overcome numerical imprecisions

Parameters :

time_series: n-d float array :

The time series data for which the regularized coherence is calculated

epsilon: float :

Small regularization parameter. Should be much smaller than any meaningful value of coherence you might encounter

alpha: float :

large regularization parameter. Should be much larger than any meaningful value of coherence you might encounter (preferably much larger than 1).

csd_method: dict, optional. :

See get_spectra() documentation for details

Returns :

frequencies, coherence :

Notes

The regularization scheme is as follows:

C_{x,y} = \frac{(\alpha f_{xx} + \epsilon)^2}{\alpha^{2}((f_{xx}+\epsilon)(f_{yy}+\epsilon))}

nitime.algorithms.coherence_reqularized_calculate(fxy, fxx, fyy, epsilon, alpha)

A regularized version of the calculation of coherence, which is more robust to numerical noise than the standard calculation.

Input to this function is in the frequency domain

Parameters :

fxy, fxx, fyy: float arrays :

The cross- and power-spectral densities of the two signals x and y

epsilon: float :

First regularization parameter. Should be much smaller than any meaningful value of coherence you might encounter

alpha: float :

Second regularization parameter. Should be much larger than any meaningful value of coherence you might encounter (preferably much larger than 1)

Returns :

float array :

The coherence values

nitime.algorithms.coherency(time_series, csd_method=None)

Compute the coherency between the spectra of n-tuple of time series. Input to this function is in the time domain

Parameters :

time_series: n*t float array :

an array of n different time series of length t each

csd_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

c : float array

This is a symmetric matrix with the coherencys of the signals. The coherency of signal i and signal j is in f[i][j]. Note that f[i][j] = f[j][i].conj()

Notes

This is an implementation of equation (1) of [Sun2005]:

R_{xy} (\lambda) = \frac{f_{xy}(\lambda)} {\sqrt{f_{xx} (\lambda) \cdot f_{yy}(\lambda)}}

[Sun2005]

F.T. Sun and L.M. Miller and M. D’Esposito(2005). Measuring temporal dynamics of functional networks using phase spectrum of fMRI data. Neuroimage, 28: 227-37.

nitime.algorithms.coherency_bavg(time_series, lb=0, ub=None, csd_method=None)

Compute the band-averaged coherency between the spectra of two time series.

Input to this function is in the time domain.

Parameters :

time_series: n*t float array :

an array of n different time series of length t each

lb, ub: float, optional :

the upper and lower bound on the frequency band to be used in averaging defaults to 1,max(f)

csd_method: dict, optional. :

See get_spectra() documentation for details

Returns :

c: float array :

This is an upper-diagonal array, where c[i][j] is the band-averaged coherency between time_series[i] and time_series[j]

Notes

This is an implementation of equation (A4) of [Sun2005]:

\bar{Coh_{xy}} (\bar{\lambda}) = \frac{\left|{\sum_\lambda{\hat{f_{xy}}}}\right|^2} {\sum_\lambda{\hat{f_{xx}}}\cdot sum_\lambda{\hat{f_{yy}}}}

[Sun2005]

F.T. Sun and L.M. Miller and M. D’Esposito(2005). Measuring temporal dynamics of functional networks using phase spectrum of fMRI data. Neuroimage, 28: 227-37.

nitime.algorithms.coherency_bavg_calculate(fxy, fxx, fyy)

Compute the band-averaged coherency between the spectra of two time series.

Input to this function is in the frequency domain.

Parameters :

fxy : float array

The cross-spectrum of the time series

fyy,fxx : float array

The spectra of the signals

Returns :

float :

the band-averaged coherency

Notes

This is an implementation of equation (A4) of [Sun2005]:

\bar{Coh_{xy}} (\bar{\lambda}) = \frac{\left|{\sum_\lambda{\hat{f_{xy}}}}\right|^2} {\sum_\lambda{\hat{f_{xx}}}\cdot sum_\lambda{\hat{f_{yy}}}}

[Sun2005]

F.T. Sun and L.M. Miller and M. D’Esposito(2005). Measuring temporal dynamics of functional networks using phase spectrum of fMRI data. Neuroimage, 28: 227-37.

nitime.algorithms.coherency_calculate(fxy, fxx, fyy)

Compute the coherency between the spectra of two time series.

Input to this function is in the frequency domain.

Parameters :

fxy : float array

The cross-spectrum of the time series

fyy,fxx : float array

The spectra of the signals

Returns :

complex array :

the frequency-band-dependent coherency

See also

coherency()

nitime.algorithms.coherency_phase_delay(time_series, lb=0, ub=None, csd_method=None)

The temporal delay calculated from the coherency phase spectrum.

Parameters :

time_series: float array :

The time-series data for which the delay is calculated.

lb, ub: float :

Frequency boundaries (in Hz), for the domain over which the delays are calculated. Defaults to 0-max(f)

csd_method : dict, optional.

See get_spectra()

Returns :

f : float array

The mid-frequencies for the frequency bands over which the calculation is done.

p : float array

Pairwise temporal delays between time-series (in seconds).

nitime.algorithms.coherency_phase_delay_bavg(time_series, lb=0, ub=None, csd_method=None)

Band-averaged phase delay between time-series

Parameters :

time_series: float array :

The time-series data

lb,ub : float, optional

Lower and upper bounds on the frequency range over which the phase delay is averaged

Returns :

p : float array

The pairwise band-averaged phase-delays between the time-series.

nitime.algorithms.coherency_phase_delay_bavg_calculate(f, fxy)

Compute the band-averaged phase delay between the spectra of two signals

Parameters :

f: float array :

The frequencies

fxy : float array

The cross-spectrum of the time series

Returns :

float :

the phase delay (in sec)

nitime.algorithms.coherency_phase_delay_calculate(f, fxy)

Compute the phase delay between the spectra of two signals. The input to this function is in the frequency domain.

Parameters :

f: float array :

The frequencies

fxy : float array

The cross-spectrum of the time series

Returns :

float array :

the phase delay (in sec) for each frequency band.

nitime.algorithms.coherency_phase_spectrum(time_series, csd_method=None)

Compute the phase spectrum of the cross-spectrum between two time series.

The parameters of this function are in the time domain.

Parameters :

time_series: n*t float array :

The time series, with t, time, as the last dimension :

Returns :

f: mid frequencies of the bands :

p: an array with the pairwise phase spectrum between the time :

series, where p[i][j] is the phase spectrum between time series[i] and :

time_series[j] :

Notes

This is an implementation of equation (3) of Sun et al. (2005) [Sun2005]:

\phi(\lambda) = arg [R_{xy} (\lambda)] = arg [f_{xy} (\lambda)]

[Sun2005]

F.T. Sun and L.M. Miller and M. D’Esposito(2005). Measuring temporal dynamics of functional networks using phase spectrum of fMRI data. Neuroimage, 28: 227-37.

nitime.algorithms.coherency_phase_spectrum_calculate(fxy)

Compute the phase spectrum of the cross-spectrum between two time series.

The parameters of this function are in the frequency domain.

Parameters :

fxy : float array

The cross-spectrum of the time series

Returns :

float :

a frequency-band-dependent measure of the phase between the two time-series

Notes

This is an implementation of equation (3) of Sun et al. (2005) [Sun2005]:

\phi(\lambda) = arg [R_{xy} (\lambda)] = arg [f_{xy} (\lambda)]

[Sun2005]

F.T. Sun and L.M. Miller and M. D’Esposito(2005). Measuring temporal dynamics of functional networks using phase spectrum of fMRI data. Neuroimage, 28: 227-37.

nitime.algorithms.coherency_regularized(time_series, epsilon, alpha, csd_method=None)

Compute a regularized measure of the coherence.

Regularization may be needed in order to overcome numerical imprecisions

Parameters :

time_series: float array :

The time series data for which the regularized coherence is calculated. Time as the last dimension.

epsilon: float :

Small regularization parameter. Should be much smaller than any meaningful value of coherence you might encounter

alpha: float :

Large regularization parameter. Should be much larger than any meaningful value of coherence you might encounter (preferably much larger than 1).

csd_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

c: float array :

This is a symmetric matrix with the coherencys of the signals. The coherency of signal i and signal j is in f[i][j]. Note that f[i][j] = f[j][i].conj()

Notes

The regularization scheme is as follows:

Coh_{xy}^R = \frac{(\alpha f_{xx} + \epsilon) ^2}{\alpha^{2}(f_{xx}+\epsilon)(f_{yy}+\epsilon)}

nitime.algorithms.coherency_reqularized_calculate(fxy, fxx, fyy, epsilon, alpha)

A regularized version of the calculation of coherency, which is more robust to numerical noise than the standard calculation

Input to this function is in the frequency domain.

Parameters :

fxy, fxx, fyy: float arrays :

The cross- and power-spectral densities of the two signals x and y

epsilon: float :

First regularization parameter. Should be much smaller than any meaningful value of coherence you might encounter

alpha: float :

Second regularization parameter. Should be much larger than any meaningful value of coherence you might encounter (preferably much larger than 1).

Returns :

float array :

The coherence values

nitime.algorithms.correlation_spectrum(x1, x2, Fs=2, norm=False)

Calculate the spectral decomposition of the correlation.

Parameters :

x1,x2: ndarray :

Two arrays to be correlated. Same dimensions

Fs: float, optional :

Sampling rate in Hz. If provided, an array of frequencies will be returned.Defaults to 2

norm: bool, optional :

When this is true, the spectrum is normalized to sum to 1

Returns :

ccn: ndarray :

The spectral decomposition of the correlation

f: ndarray :

ndarray with the frequencies

Notes

This method is described in full in [Cordes2000]

[Cordes2000]

D Cordes, V M Haughton, K Arfanakis, G J Wendt, P A Turski,

C H Moritz, M A Quigley, M E Meyerand (2000). Mapping functionally related regions of brain with functional connectivity MR imaging. AJNR American journal of neuroradiology 21:1636-44

nitime.algorithms.fir(timeseries, design)

Calculate the FIR (finite impulse response) HRF, according to [Burock2000]

Parameters :

timeseries : float array

timeseries data

design : int array

This is a design matrix. It has to have shape = (number of TRS, number of conditions * length of HRF)

The form of the matrix is:

A B C ...

where A is a (number of TRs) x (length of HRF) matrix with a unity matrix placed with its top left corner placed in each TR in which event of type A occured in the design. B is the equivalent for events of type B, etc.

Returns :

HRF: float array :

HRF is a numpy array of 1X(length of HRF * number of conditions) with the HRFs for the different conditions concatenated. This is an estimate of the linear filters between the time-series and the events described in design.

Notes

Implements equation 4 in[Burock2000]_:

\hat{h} = (X^T X)^{-1} X^T y

[Burock2000]

M.A. Burock and A.M.Dale (2000). Estimation and Detection of Event-Related fMRI Signals with Temporally Correlated Noise: A Statistically Efficient and Unbiased Approach. Human Brain Mapping, 11:249-260

nitime.algorithms.freq_domain_xcorr(tseries, events, t_before, t_after, Fs=1)

Calculates the event related timeseries, using a cross-correlation in the frequency domain.

Parameters :

tseries: float array :

Time series data with time as the last dimension

events: float array :

An array with time-resolved events, at the same sampling rate as tseries

t_before: float :

Time before the event to include

t_after: float :

Time after the event to include

Fs: float :

Sampling rate of the time-series (in Hz)

Returns :

xcorr: float array :

The correlation function between the tseries and the events. Can be interperted as a linear filter from events to responses (the time-series) of an LTI.

nitime.algorithms.freq_domain_xcorr_zscored(tseries, events, t_before, t_after, Fs=1)

Calculates the z-scored event related timeseries, using a cross-correlation in the frequency domain.

Parameters :

tseries: float array :

Time series data with time as the last dimension

events: float array :

An array with time-resolved events, at the same sampling rate as tseries

t_before: float :

Time before the event to include

t_after: float :

Time after the event to include

Fs: float :

Sampling rate of the time-series (in Hz)

Returns :

xcorr: float array :

The correlation function between the tseries and the events. Can be :

interperted as a linear filter from events to responses (the time-series) :

of an LTI. Because it is normalized to its own mean and variance, it can be :

interperted as measuring statistical significance relative to all :

time-shifted versions of the events. :

nitime.algorithms.gauss_white_noise(npts)

Gaussian white noise.

XXX - incomplete.

nitime.algorithms.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.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.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 tapered complex spectra, with K tapers

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)y_k^x(f)][d_k^y(f)(y_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.multi_taper_csd(s, Fs=6.2831853071795862, BW=None, low_bias=True, adaptive=False, sides='default')

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 :

BW: float, 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).

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)

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

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 :

BW: float, 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).

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, ssigma_or_nu) : ndarrays

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

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

nitime.algorithms.my_freqz(b, a=1.0, Nfreqs=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)

Nfreqs : 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#Linear_constant-coefficient_difference_equation

nitime.algorithms.periodogram(s, Fs=6.2831853071795862, 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

Notes

setting dw = 2*PI/N, then the integral from -PI, PI (or 0,PI) of PSD/(2PI) will be nearly equal to sxx(0), where sxx is the autocovariance function of s(n). By definition, sxx(0) = E{s(n)s*(n)} ~ (s*s.conj()).mean()

nitime.algorithms.periodogram_csd(s, Sk=None, N=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).

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 that is reshaped to (M,N), the output is (M,M,N)

Notes

setting dw = 2*PI/N, then the integral from -PI, PI (or 0,PI) of PSD/(2PI) will be nearly equal to sxy(0), where sxx is the crosscovariance function of s1(n), s2(n). By definition, sxy(0) = E{s1(n)s2*(n)} ~ (s1*s2.conj()).mean()

nitime.algorithms.wfmorlet_fft(f0, sd, samplingrate, ns=5, nt=None)

returns a complex morlet wavelet in the frequency domain

Parameters :

f0 : center frequency

sd : standard deviation of center frequency sampling_rate : samplingrate ns : window length in number of stanard deviations nt : window length in number of sample points

nitime.algorithms.wlogmorlet(f0, sd, sampling_rate, ns=5, normed='area')

returns a complex log morlet wavelet in the time domain

Parameters :

f0 : center frequency

sd : standard deviation of frequency sampling_rate : samplingrate ns : window length in number of stanard deviations

nitime.algorithms.wlogmorlet_fft(f0, sd, sampling_rate, ns=5, nt=None)

returns a complex log morlet wavelet in the frequency domain

Parameters :

f0 : center frequency

sd : standard deviation sampling_rate : samplingrate ns : window length in number of stanard deviations nt : window length in number of sample points

nitime.algorithms.wmorlet(f0, sd, sampling_rate, ns=5, normed='area')

returns a complex morlet wavelet in the time domain

Parameters :

f0 : center frequency

sd : standard deviation of frequency sampling_rate : samplingrate ns : window length in number of stanard deviations

nitime.algorithms.yule_AR_est(s, order, Nfreqs, sxx=None, sides='onesided', system=False)

Finds the parameters for an autoregressive model of order norder of the process s. Using these parameters, an estimate of the PSD is calculated from [-PI,PI) in Nfreqs, or [0,PI] in {N/2+1}freqs. Uses the basic Yule Walker system of equations, and a baised estimate of sxx (unless sxx is provided).

The model for the autoregressive process takes this convention: s[n] = a1*s[n-1] + a2*s[n-2] + ... aP*s[n-P] + v[n]

where v[n] is a zero-mean white noise process with variance=sigma_v

Parameters :

s : ndarray

The sampled autoregressive random process

order : int

The order P of the AR system

Nfreqs : int

The number of spacings on the frequency grid from [-PI,PI). If sides==’onesided’, Nfreqs/2+1 frequencies are computed from [0,PI]

sxx : ndarray (optional)

An optional, possibly unbiased estimate of the autocovariance of s

sides : str (optional)

Indicates whether to return a one-sided or two-sided PSD

system : bool (optional)

If True, return the AR system parameters, sigma_v and a{k}

Returns :

(w, ar_psd) :

w : Array of normalized frequences from [-.5, .5) or [0,.5]

ar_psd : A PSD estimate computed by sigma_v / |1-a(f)|**2 , where

a(f) = DTFT(ak)