meegkit.utils#

Utility functions.

`auditory`	Tools for working with auditory signals.
`coherence`	Signal coherence tools.
`covariances`	Covariance calculation.
`denoise`	Denoising utilities.
`matrix`	Matrix operation utility functions.
`sig`	Audio and signal processing tools.
`stats`	Statistics utilities.

Auditory#

Tools for working with auditory signals.

class meegkit.utils.auditory.AuditoryFilterbank(sfreq, b=1.019, order=1, q=9.26449, min_bw=24.7)#

Bases: GammatoneFilterbank

Special case of Gammatone filterbank with preset center frequencies.

Methods

filter(X)

Filter X along its last dimension.

__init__(sfreq, b=1.019, order=1, q=9.26449, min_bw=24.7)#

class meegkit.utils.auditory.GammatoneFilterbank(sfreq, cf, b=1.019, order=1, q=9.26449, min_bw=24.7)#

Bases: object

Gammatone Filterbank.

This class computes the filter coefficients for a bank of Gammatone filters. These filters were defined by Patterson and Holdworth for simulating the cochlea, and originally implemented in [1].

Parameters:

sfreq (float) – Sampling frequency of the signals to filter.
cf (array_like) – Center frequencies of the filterbank.
b (float) – beta of the gammatone filters (default=1.019).
order (int) – Order (default=1).
q (float) – Q-value of the ERB (default=9.26449).
min_bw (float) – Minimum bandwidth of an ERB.

Notes

The python was adapted from Alexandre Chabot-Leclerc’s pambox, and Jason Heeris’ gammatone toolbox: achabotl/pambox detly/gammatone

References

[1]

Slaney, M. (1993). An efficient implementation of the Patterson-Holdsworth auditory filter bank. Apple Computer, Perception Group, Tech. Rep, 35(8).

Methods

filter(X)

Filter X along its last dimension.

__init__(sfreq, cf, b=1.019, order=1, q=9.26449, min_bw=24.7)#

filter(X)#

Filter X along its last dimension.

Process an input waveform with a gammatone filter bank. This function takes a single sound vector, and returns an array of filter outputs, one channel per row.

Parameters:: X (ndarray, shape=(n_chans, n_times)) – Signal to filter.
Returns:: Filtered signals with shape (M, N), where M is the number of channels, and N is the input signal’s number of samples.
Return type:: ndarray

meegkit.utils.auditory.erb2hz(erb)#: Convert ERB-rate values to the corresponding frequencies in Hz.

meegkit.utils.auditory.erb_bandwidth(fc)#

Bandwidth of an Equivalent Rectangular Bandwidth (ERB).

Parameters:: fc (ndarray) – Center frequency, or center frequencies, of the filter.
Returns:: Equivalent rectangular bandwidth of the filter(s).
Return type:: ndarray or float

meegkit.utils.auditory.erbspace(flow, fhigh, n)#: Generate n equidistantly spaced points on ERB scale.

meegkit.utils.auditory.hz2erb(f)#

Convert frequencies to the corresponding ERB-rates.

Notes

There is a round-off error in the Glasberg & Moore paper, as 1000 / (24.7 * 4.37) * log(10) = 21.332 and not 21.4 as is stated.

Coherence#

Signal coherence tools.

Compute 2D, 1D, and 0D bicoherence, polycoherence, bispectrum, and polyspectrum.

Bicoherence is a measure of the degree of phase coupling between different frequency components in a signal. It’s essentially a normalized form of the bispectrum, which itself is a Fourier transform of the third-order cumulant of a time series: - 2D bicoherence is the most common form, where one looks at a

two-dimensional representation of the phase coupling between different frequencies. It’s a function of two frequency variables.

1D bicoherence would imply a slice or a specific condition in the 2D bicoherence, reducing it to a function of a single frequency variable. It simplifies the analysis by looking at the relationship between a particular frequency and its harmonics or other relationships
0D bicoherence would imply a single value representing some average or overall measure of phase coupling in the signal. It’s a highly condensed summary, which might represent the average bicoherence over all frequency pairs, for instance.

meegkit.utils.coherence.compute_spectrogram(X, sfreq, **kwargs)#

Compute the complex spectrogram of X.

Simple wrapper around scipy.signal.spectrogram.

Returns:

f (ndarray, shape=(n_freqs,)) – Positive frequencies.
t (ndarray, shape=(n_timesteps,)) – Time axis.
S (ndarray, shape=(n_chans, n_timesteps, n_freqs)) – Complex spectrogram (one-sided).

meegkit.utils.coherence.cross_coherence(x1, x2, sfreq, norm=2, **kwargs)#

Compute the bispectral cross-coherence between two signals of same length.

Code adapted from [2].

Parameters:

x1 (array-like, shape=([n_channels, ]n_samples)) – First signal.
x2 (array-like, shape=([n_channels, ]n_samples)) – Second signal.
sfreq (float) – Sampling sfreq.
norm (int | None) – Norm (default=2). If None, return bispectrum.

Returns:

f1 (array-like) – Frequency axis.
B (array-like) – Bicoherence between s1 and s2.

References

[1]

http://wiki.fusenet.eu/wiki/Bicoherence

[2]

https://stackoverflow.com/a/36725871

meegkit.utils.coherence.norm_spectrum(spec, P1, P2, P12, time_axis=-2)#

Compute bicoherence from bispectrum.

For formula see [1].

Parameters:

spec (ndarray, shape=(n_chans, n_freqs)) – Polyspectrum.
P1 (array-like, shape=(n_chans, n_times, n_freqs_f1)) – Spectrum evaluated at f1.
P2 (array-like, shape=(n_chans, n_times, n_freqs_f2)) – Spectrum evaluated at f2.
P12 (array-like, shape=(n_chans, n_times, n_freqs)) – Spectrum evaluated at f1 + f2.
time_axis (int) – Time axis.

Returns:

coh – Polycoherence.

Return type:

ndarray, shape=(n_chans,)

References

[1]

https://en.wikipedia.org/wiki/Bicoherence

meegkit.utils.coherence.plot_polycoherence(freq1, freq2, bicoh, ax=None)#: Plot polycoherence.

meegkit.utils.coherence.plot_polycoherence_1d(freq, coh)#: Plot polycoherence for fixed frequencies.

meegkit.utils.coherence.plot_signal(t, signal, ax=None)#: Plot signal and spectrum.

meegkit.utils.coherence.polycoherence_0d(X, sfreq, freqs, norm=2, synthetic=None, **kwargs)#

Polycoherence between freqs and sum of freqs.

Parameters:

X (ndarray, shape=(n_channels, n_samples)) – Input data.
sfreq (float) – Sampling rate.
freqs (list[float]) – Fixed frequencies.
norm (int | None) – Norm (default=2).
synthetic (tuple(float, float, float)) – Used for synthetic signal for some frequencies (freq, amplitude, phase), freq must coincide with the first fixed frequency.
**kwargs (dict) – Additional parameters passed to scipy.signal.spectrogram. Important parameters are nperseg, noverlap, nfft.

Returns:

B – Polycoherence

Return type:

ndarray, shape=(n_channels,)

meegkit.utils.coherence.polycoherence_1d(X, sfreq, f2, norm=2, synthetic=None, **kwargs)#

1D polycoherence as a function of f1 and at least one fixed frequency f2.

Parameters:

X (ndarray) – Input data.
sfreq (float) – Sampling rate
f2 (list[float]) – Fixed frequencies.
norm (int | None) – Norm (default=2).
synthetic (tuple(float, float, float)) – Used for synthetic signal for some frequencies (freq, amplitude, phase), freq must coincide with the first fixed frequency.
**kwargs – Additional parameters passed to scipy.signal.spectrogram. Important parameters are nperseg, noverlap, nfft.

Returns:

freq (ndarray, shape=(n_freqs_f1,)) – Frequencies
B (ndarray, shape=(n_channels, n_freqs_f1)) – 1D polycoherence.

meegkit.utils.coherence.polycoherence_1d_sum(X, sfreq, fsum, *ofreqs, norm=2, synthetic=None, **kwargs)#

1D polycoherence with fixed frequency sum fsum as a function of f1.

Returns polycoherence for frequencies ranging from 0 up to fsum.

Parameters:

X (ndarray) – Input data.
sfreq (float) – Sampling rate.
fsum (float) – Fixed frequency sum.
ofreqs (list[float]) – Fixed frequencies.
norm (int or None) – If 2 - return polycoherence, n0 = n1 = n2 = 2 (default)
synthetic (tuple(float, float, float) | None) – Used for synthetic signal for some frequencies (freq, amplitude, phase), freq must coincide with the first fixed frequency.

Returns:

freq (ndarray, shape=(n_freqs,)) – Frequencies.
B (ndarray, shape=(n_channels, n_freqs)) – Polycoherence for f1+f2=fsum.

meegkit.utils.coherence.polycoherence_2d(X, sfreq, ofreqs=None, norm=2, flim1=None, flim2=None, synthetic=None, **kwargs)#

2D polycoherence between freqs and their sum as a function of f1 and f2.

2D bicoherence is the most common form, where one looks at a two-dimensional representation of the phase coupling between different frequencies. It is a function of two frequency variables.

2D polycoherence as a function of f1 and f2, ofreqs are additional fixed frequencies.

Parameters:

X (ndarray) – Input data.
sfreq (float) – Sampling rate.
ofreqs (list[float]) – Fixed frequencies.
norm (int or None) – If 2 - return polycoherence (default), else return polyspectrum.
flim1 (tuple | None) – Frequency limits for f1. If None, it is set to (0, nyquist / 2)
flim2 (tuple | None) – Frequency limits for f2.

Returns:

freq1 (ndarray, shape=(n_freqs_f1,)) – Frequencies for f1.
freq2 (ndarray, shape=(n_freqs_f2,)) – Frequencies for f2.
B (ndarray, shape=([n_chans, ]n_freqs_f1, n_freqs_f2)) – Polycoherence.

meegkit.utils.coherence.synthetic_signal(t, synthetic)#

Create a synthetic complex signal spectrum.

Parameters:

t (array-like) – Time.
synthetic (list[tuple(float, float, float)]) – List of tuples with (freq, amplitude, phase).

Return type:

Complex signal.

Covariances#

Covariance calculation.

meegkit.utils.covariances.block_covariance(data, window=128, overlap=0.5, padding=True, estimator='cov')#

Compute blockwise covariance.

Parameters:

data (array, shape=(n_chans, n_samples)) – Input data (must be 2D)
window (int) – Window size.
overlap (float) – Overlap between successive windows.

Returns:

cov – Block covariance.

Return type:

array, shape=(n_blocks, n_chans, n_chans)

meegkit.utils.covariances.convmtx(V, n)#

Generate a convolution matrix.

CONVMTX(V,N) returns the convolution matrix for vector V. If V is a column vector and X is a column vector of length N, then CONVMTX(V,N)*X is the same as CONV(V,X). If R is a row vector and X is a row vector of length N, then X*CONVMTX(R,N) is the same as CONV(R,X).

Given a vector V of length N, an N+n-1 by n convolution matrix is generated of the following form:

    |  V(0)  0      0     ...      0    |
    |  V(1) V(0)    0     ...      0    |
    |  V(2) V(1)   V(0)   ...      0    |
X = |   .    .      .              .    |
    |   .    .      .              .    |
    |   .    .      .              .    |
    |  V(N) V(N-1) V(N-2) ...  V(N-n+1) |
    |   0   V(N)   V(N-1) ...  V(N-n+2) |
    |   .    .      .              .    |
    |   .    .      .              .    |
    |   0    0      0     ...    V(N)   |

That is, V is assumed to be causal, and zero-valued after N.

Parameters:

V (array, shape=(N,) or(N, 1) or (1, N)) – Input vector.
n (int)

Returns:

Return type:

array, shape=(N * n - 1, n)

Examples

Generate a simple convolution matrix:

>>> h = [1, 2, 1]
>>> convmtx(h,7)
np.array(
    [[1. 2. 1. 0. 0. 0.]
     [0. 1. 2. 1. 0. 0.]
     [0. 0. 1. 2. 1. 0.]
     [0. 0. 0. 1. 2. 1.]]
)

meegkit.utils.covariances.cov_lags(X, Y, shifts=None)#

Empirical covariance of the joint array [X, Y] with lags.

Parameters:

X (array, shape=(n_times, n_chans_x[, n_trials])) – Time shifted data.
Y (array, shape=(n_times, n_chans_y[, n_trials])) – Reference data.
shifts (array, shape=(n_shifts,)) – Positive lag means X is delayed relative to Y.

Returns:

C (array, shape=(n_chans_x + n_chans_y, n_chans_x + n_chans_y, n_shifts)) – Covariance matrix (3D if n_shifts > 1).
tw (float) – Total weight.
m (int) – Number of columns in X.

See also

relshift, tscov, tsxcov

meegkit.utils.covariances.nonlinear_eigenspace(L, k, alpha=1)#

Nonlinear eigenvalue problem: total energy minimization.

This example is motivated in [1] and was adapted from the manopt toolbox in Matlab.

TODO : check this

Parameters:

L (array, shape=(n_channels, n_channels)) – Discrete Laplacian operator: the covariance matrix.
alpha (float) – Given constant for optimization problem.
k (int) – Determines how many eigenvalues are returned.

Returns:

Xsol (array, shape=(n_channels, n_channels)) – Eigenvectors.
S0 (array) – Eigenvalues.

References

[1]

“A Riemannian Newton Algorithm for Nonlinear Eigenvalue Problems”, Zhi Zhao, Zheng-Jian Bai, and Xiao-Qing Jin, SIAM Journal on Matrix Analysis and Applications, 36(2), 752-774, 2015.

meegkit.utils.covariances.pca(cov, max_comps=None, thresh=0)#

PCA from covariance.

Parameters:

cov (array, shape=(n_chans, n_chans)) – Covariance matrix.
max_comps (int | None) – Maximum number of components to retain after decomposition. None (the default) keeps all suprathreshold components (see thresh).
thresh (float) – Discard components below this threshold.

Returns:

V (array, shape=(max_comps, max_comps)) – Eigenvectors (matrix of PCA components).
d (array, shape=(max_comps,)) – PCA eigenvalues

meegkit.utils.covariances.regcov(Cxy, Cyy, keep=None, threshold=0)#

Compute regression matrix from cross covariance.

Parameters:

Cxy (array) – Cross-covariance matrix between data and regressor.
Cyy (array) – Covariance matrix of regressor.
keep (array) – Number of regressor PCs to keep (default=None, which keeps all).
threshold (float) – Eigenvalue threshold for discarding regressor PCs (default=0, which keeps all components).

Returns:

R – Matrix to apply to regressor to best model data.

Return type:

array

meegkit.utils.covariances.tscov(X, shifts=None, weights=None, assume_centered=True)#

Time shift covariance.

This function calculates, for each pair [X[i], X[j]] of columns of X, the cross-covariance matrix between the time-shifted versions of X[i].

Parameters:

X (array, shape=(n_times, n_chans[, n_trials])) – Data, can be 1D, 2D or 3D.
shifts (array, shape=(n_shifts,) | None) – Array of time shifts.
weights (array, shape=(n_times,) or) –

shape=(n_times, 1, n_trials), or
shape=(n_times, n_chans, n_trials)

Weights, 1D (if X is 1D or 2D) or 2D (if X is 3D). The weights are not shifted.
assume_centered (bool) – If False, remove the mean of X before computing the covariance (default=True).

Returns:

C (array, shape=(n_chans * n_shifts, n_chans * n_shifts)) – Covariance matrix. This matrix is made up of a (n_times, n_times) matrix of submatrices of dimensions (n_shifts, n_shifts).
tw (array) – Total weight (C/tw is the normalized covariance).

meegkit.utils.covariances.tsxcov(X, Y, shifts=None, weights=None, assume_centered=True)#

Calculate cross-covariance of X and time-shifted Y.

This function calculates, for each pair of columns (Xi, Yj) of X and Y, the scalar products between Xi and time-shifted versions of Yj.

Output is a 2D matrix with dimensions .

Parameters:

X (arrays, shape=(n_times, n_chans[, n_trials])) – Data to cross correlate. X can be 1D, 2D or 3D.
Y (arrays, shape=(n_times, n_chans[, n_trials])) – Data to cross correlate. X can be 1D, 2D or 3D.
shifts (array) – Time shifts.
weights (array) – The weights that are applied to X. 1D (if X is 1D or 2D) or 2D (if X is 3D).
assume_centered (bool) – If False, remove the mean of X before computing the covariance (default=True).

Returns:

C (array, shape=(n_chans_x, n_chans_y * n_shifts)) – Cross-covariance matrix.
tw (total weight)

Denoising#

Denoising utilities.

meegkit.utils.denoise.demean(X, weights=None, return_mean=False, inplace=False)#

Remove weighted mean over rows (samples).

Parameters:

X (array, shape=(n_samples, n_channels[, n_trials])) – Data.
weights (array, shape=(n_samples))
return_mean (bool) – If True, also return signal mean (default=False).
inplace (bool) – If True, save the resulting array in X (default=False).

Returns:

X (array, shape=(n_samples, n_channels[, n_trials])) – Centered data.
mn (array) – Mean value.

meegkit.utils.denoise.find_outlier_samples(X, toobig1, toobig2=[])#: Find outlier trials using an absolute threshold.

meegkit.utils.denoise.find_outlier_trials(X, thresh=None, show=True)#

Find outlier trials.

For example thresh=2 rejects trials that deviate from the mean by more than twice the average deviation from the mean.

Parameters:

X (ndarray, shape=(n_times, n_chans[, n_trials])) – Data array.
thresh (float or array of floats) – Keep trials less than thresh from mean.
show (bool) – If True (default), plot trial deviations before and after.

Returns:

bads (list of int) – Indices of trials to reject.
d (array) – Relative deviations from mean.

meegkit.utils.denoise.mean_over_trials(X, weights=None)#: Compute weighted mean over trials (last dimension).

meegkit.utils.denoise.wpwr(X, weights=None)#: Weighted power.

Matrix manipulation#

Matrix operation utility functions.

meegkit.utils.matrix.fold(X, epoch_size)#: Fold 2D (n_times, n_channels) X into 3D (n_times, n_chans, n_trials).

meegkit.utils.matrix.matmul3d(X, mixin)#

Apply mixing matrix to each trial of X.

This is a simple wrapper or np.einsum.

Parameters:

X (array, shape=(n_samples, n_chans, n_trials))
mixing (array, shape=(n_chans, n_components))

Returns:

proj – Projection.

Return type:

array, shape=(n_samples, n_components, n_trials)

meegkit.utils.matrix.multishift(X, shifts, fill_value=0, axis=0, keep_dims=False, reshape=False, solution='full')#

Apply several shifts along specified axis.

If shifts has multiple values, the output will contain one shift per page. Shifted data are padded with fill_value.

Parameters:

X (array, shape=(n_samples[, n_chans][, n_trials])) – Array to shift.
shifts (array, shape=(n_shifts,)) – Array of shifts. Positive shifts mean that X is ‘delayed’ in time (i.e. y[shift] = X[0]). Conversely, a negative shift means that X is ‘advanced’ (i.e. y[0] =
fill_value (float | np.nan) – Value to pad output axis by.
axis (int, optional) – The axis along which elements are shifted.
keep_dims (bool) – If True, keep singleton dimensions in output.
reshape (bool) – If True, concatenate channels and lags, yielding an array of shape (n_samples, n_chans*n_shifts[, n_trials])
solution ({'valid', 'full'}) – If valid, the output’s is cropped along axis by n_shifts in order to remove edge artifacts. If full, the output has the same size as X.

Returns:

y – Shifted array.

Return type:

array, shape=(n_samples[, n_chans][, n_trials], n_shifts)

See also

relshift, shift, shiftnd

meegkit.utils.matrix.multismooth(X, smooths, axis=0, keep_dims=False)#

Apply several shifts along specified axis.

If shifts has multiple values, the output will contain one shift per trial. Shifted data are padded with fill_value.

Parameters:

X (array, shape=(n_samples[, n_epochs][, n_trials])) – Array to shift.
smooths (array) – Array of smoothing values (in samples).
axis (int, optional) – The axis along which elements are shifted (default=0).
keep_dims (bool) – If True, keep singleton dimensions in output.

Returns:

y – Shifted array.

Return type:

array, shape=(n_samples[, n_epochs][, n_trials], n_shifts)

See also

multishift, smooth

meegkit.utils.matrix.normcol(X, weights=None, return_norm=False)#

Normalize each column so that its weighted mean square value is 1.

If X is 3D, pages are concatenated vertically before calculating the norm.

Weight should be either a column vector, or a matrix (2D or 3D) of same size as X.

Parameters:

X (array) – X to normalize.
weights (array) – Weights.
return_norm (bool) – If True, also return norm vector.

Returns:

X_norm (array) – Normalized X.
norm (array) – Norm.

meegkit.utils.matrix.relshift(X, ref, shifts, fill_value=0, axis=0)#

Create shifted versions of X relative to ref with padding.

ref is replicated to have the same shape as X and padded accordingly.

Parameters:

X (array, shape=(n_samples[, n_epochs][, n_trials])) – Array to shift.
ref (array, shape=(n_samples[, n_epochs][, n_trials])) – Reference array against which X is shifted.
shifts (array | int) – Array of shifts. Positive shifts mean that X is ‘delayed’ in time (i.e. y[shift] = X[0]). Conversely, a negative shift means that X is ‘advanced’ (i.e. y[0] =
fill_value (float) – Value to pad output axis by.
axis (int) – The axis along which elements are shifted.

Returns:

y (array, shape=(n_samples[, n_epochs][, n_trials], n_shifts)) – Shifted array.
y_ref (array, shape=(n_samples[, n_epochs][, n_trials], n_shifts)) – Reference array, repeated to match y.shape. Padding matches that of y.

See also

multishift, shift, shiftnd

meegkit.utils.matrix.shift(X, shift, fill_value=0, axis=0)#

Shift array along its first, second or last dimension.

Output is padded by fill_value.

Parameters:

X (array, shape=(n_samples[, n_epochs][, n_trials])) – Multidimensional input array.
shift (int) – The number of places by which elements are shifted along axis. Positive shifts mean that X is ‘delayed’ in time (i.e. y[shift] = X[0]). Conversely, a negative shift means that X is ‘advanced’ (i.e. y[0] = X[shift]).
fill_value (float) – Value to pad output axis by.
axis (int, optional) – The axis along which elements are shifted.

Returns:

y – Output array, with the same shape as X.

Return type:

array

See also

relshift, multishift, shiftnd

meegkit.utils.matrix.shiftnd(X, shift, fill_value=0, axis=None)#

Roll array elements along a given axis with padding.

Elements off the end of the array are treated as zeros. This function is slower than function:shift, so prefer the latter when possible.

Parameters:

X (array) – Multidimensional input array.
shift (int) – The number of places by which elements are shifted along axis.
fill_value (float) – Value to pad output axis by.
axis (int, optional) – The axis along which elements are shifted. By default, the array is flattened before shifting, after which the original shape is restored.

Returns:

y – Output array, with the same shape as X.

Return type:

array, (n_samples, [n_epochs, ][n_trials, ])

See also

np.roll: Elements that roll off one end come back on the other.
np.rollaxis: Roll the specified axis backwards, until it lies in a given position.

Examples

>>> x = np.arange(10)
>>> shiftnd(x, 2)
array([0, 0, 0, 1, 2, 3, 4, 5, 6, 7])
>>> x2 = np.reshape(x, (2,5))
>>> x2
array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])
>>> shiftnd(x2, 1)
array([[0, 0, 1, 2, 3],
       [4, 5, 6, 7, 8]])
>>> shiftnd(x2, -2)
array([[2, 3, 4, 5, 6],
       [7, 8, 9, 0, 0]])
>>> shiftnd(x2, 1, axis=0)
array([[0, 0, 0, 0, 0],
       [0, 1, 2, 3, 4]])
>>> shiftnd(x2, -1, axis=0)
array([[5, 6, 7, 8, 9],
       [0, 0, 0, 0, 0]])
>>> shiftnd(x2, 1, axis=1)
array([[0, 0, 1, 2, 3],
       [0, 5, 6, 7, 8]])

meegkit.utils.matrix.sliding_window(data, window, step=1, axis=-1, copy=True)#

Calculate a sliding window over a signal.

Parameters:

data (array) – The array to be slided over.
window (int) – The sliding window size.
step (int) – The sliding window stepsize (default=1).
axis (int) – The axis to slide over (defaults=-1).
copy (bool) – Return strided array as copy to avoid sideffects when manipulating the output array.

Returns:

data – A matrix whose last dimension corresponds to the window size, and the second-to-last dimension corresponds to the number of slices.

Return type:

array, shape=(…, n_windows, window_size)

Notes

Be wary of setting copy to False as undesired side effects with the output values may occur.

Examples

>>> a = numpy.array([1, 2, 3, 4, 5])
>>> sliding_window(a, size=3)
array([[1, 2, 3],
       [2, 3, 4],
       [3, 4, 5]])

>>> sliding_window(a, size=3, stepsize=2)
array([[1, 2, 3],
       [3, 4, 5]])

meegkit.utils.matrix.theshapeof(X)#: Return the shape of X.

meegkit.utils.matrix.unfold(X)#: Unfold 3D X into 2D (concatenate trials).

meegkit.utils.matrix.unsqueeze(X)#: Append singleton dimensions to an array.

meegkit.utils.matrix.widen_mask(mask, widen=4, axis=0)#

Widen each ‘bad’ section of a binary mask by n cells.

Parameters:

mask (array) – Masking array. If an element of the mask is True (or 1), the corresponding element of the associated array is masked (marked as invalid).
widen (int) – Number of cells to widen mask by.
axis (int) – Axis to operate on.

Returns:

out – Widened mask, of same shape as input mask.

Return type:

array

Examples

>>> test = widen_mask(np.array([False, False, False, True, False], 1)
>>> print(test)
[False False False True True]]

Signal#

Audio and signal processing tools.

meegkit.utils.sig.gaussfilt(data, srate, f, fwhm, n_harm=1, shift=0, return_empvals=False, show=False)#

Narrow-band filter via frequency-domain Gaussian.

Empirical frequency and FWHM depend on the sampling rate and the number of time points, and may thus be slightly different from the requested values.

Parameters:

data (ndarray) – EEG data, shape=(n_samples, n_channels[, …])
srate (int) – Sampling rate in Hz.
f (float) – Break frequency of filter.
fhwm (float) – Standard deviation of filter, defined as full-width at half-maximum in Hz.
n_harm (int) – Number of harmonics of the frequency to consider.
shift (int) – Amount shift peak frequency by (only useful when considering harmonics, otherwise leave to 0).
return_empvals (bool) – Return empirical values (default: False).
show (bool) – Set to True to show the frequency-domain filter shape.

Returns:

filtdat (ndarray) – Filtered data.
empVals (float) – The empirical frequency and FWHM.

meegkit.utils.sig.hilbert_envelope(x)#

Calculate the Hilbert envelope of a signal.

Parameters:: x (array) – Signal on which to calculate the hilbert envelope. The calculation is done along the last axis (i.e. axis=-1).
Return type:: ndarray

meegkit.utils.sig.lowpass_env_filtering(x, cutoff=150.0, n=1, sfreq=22050)#

Low-pass filters a signal using a Butterworth filter.

Parameters:

x (ndarray)
cutoff (float, optional) – Cut-off frequency of the low-pass filter, in Hz. The default is 150 Hz.
n (int, optional) – Order of the low-pass filter. The default is 1.
sfreq (float, optional) – Sampling frequency of the signal to filter. The default is 22050 Hz.

Returns:

Low-pass filtered signal.

Return type:

ndarray

meegkit.utils.sig.modulation_index(phase, amp, n_bins=18)#

Compute the Modulation Index (MI) between two signals.

MI is a measure of the amount of phase-amplitude coupling. Phase angles are expected to be in radians [1]. MI is derived from the Kullbach-Leibner distance, a measure for the disparity of two distributions, which is also returned here. MI is recommend the modulation index for noisy and short data epochs with unknown forms of coupling [2].

Parameters:

phase (array) – Phase vector, in radians.
amp (array) – Amplitude vector.
n_bins (int) – Number of bins in which to discretize phase (default=18 bins, giving a 20-degree resolution).

Returns:

MI (array) – Tort’s Modulation index.
KL (array) – Kullbach-Leibner distance.

Examples

>> phas = np.random.rand(100, 1) * 2 * np.pi - np.pi >> ampl = rng.standard_normal((100, 1) * 30 + 100 >> MI, KL = modulation_index(phas, ampl)

Notes

Phase and amplitude can be derived directly from any time series through the analytic signal: >> analytic_signal = hilbert(filtered_data) >> phase = np.phase(analytic_signal) >> amplitude = np.abs(analytic_signal)

MI can be subjected to permutation testing to assess significance. For permutation testing, the observed coupling value is compared to a distribution of shuffled coupling values. Shuffled coupling values are constructed by calculating the coupling value between the original phase time series and a permuted amplitude time series. The permuted amplitude time series can be constructed by cutting the amplitude time series at a random time point and reversing the order of both parts [2]. The observed coupling value is standardized to the distribution of the shuffled coupling values according to the following formula:

MI_z = (MI_observed − µ_MI_shuffled) / σ_MI_shuffled

where μ denotes the mean and σ the standard deviation. Only when the observed phase-locking value is larger than 95% of shuffled values, it is defined as significant. See [2] for details.

References

[1]

Tort AB, Komorowski R, Eichenbaum H, Kopell N. Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J Neurophysiol. 2010 Aug;104(2):1195-210. doi: 10.1152/jn.00106.2010.

[2] (1,2,3)

Huelsemann, M. J., Naumann, E., & Rasch, B. (2018). Quantification of Phase-Amplitude Coupling in Neuronal Oscillations: Comparison of Phase-Locking Value, Mean Vector Length, and Modulation Index. bioRxiv, 290361.

meegkit.utils.sig.prony(h, nb, na)#

Prony’s method for time-domain IIR filter design.

[B,A] = PRONY(H, NB, NA) finds a filter with numerator order NB, denominator order NA, and having the impulse response in vector H. The IIR filter coefficients are returned in length NB+1 and NA+1 row vectors B and A, ordered in descending powers of Z. H may be real or complex.

If the largest order specified is greater than the length of H, H is padded with zeros.

Parameters:

h (array) – Impulse response.
nb (int) – Numerator order.
na (int) – Denominator order.

References

T.W. Parks and C.S. Burrus, Digital Filter Design, John Wiley and Sons, 1987, p226.

Notes

meegkit.utils.sig.slope_sum(x, w: int, axis=0)#

Slope sum function.

The discrete version of the Teager-Kaiser operator is computed according to:

y[n] = \sum_{i=k-w}^k (y_i -y_{i-1})

Parameters:

x (array, shape=(n_samples[, n_channels][, n_trials])) – Input data.
w (int) – Window.

References

https://ieeexplore.ieee.org/document/8527395

meegkit.utils.sig.smooth(x, window_len, window='square', axis=0, align='left')#

Smooth a signal using a window with requested size along a given axis.

This method is based on the convolution of a scaled window with the signal. The signal is prepared by introducing reflected copies of the signal (with the window size) in both ends so that transient parts are minimized in the beginning and end part of the output signal.

Parameters:

x (array) – The input signal.
window_len (float) – The dimension of the smoothing window (in samples). Can be fractionary.
window (str) – The type of window from ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’ flat window will produce a moving average smoothing.
axis (int) – Axis along which smoothing will be applied (default=0).
align ({'left' | 'center'}) – If left (default), the convolution if computed in a causal way by shifting the output of a normal convolution by the kernel size. If center, the center of the impulse is used as the location where the convolution is summed.

Returns:

y – The smoothed signal.

Return type:

array

Examples

>> t = linspace(-2, 2, 0.1) >> x = sin(t) + randn(len(t)) * 0.1 >> y = smooth(x, 2)

See also

np.hanning, np.hamming, np.bartlett, np.blackman, np.convolve, scipy.signal.lfilter

Notes

length(output) != length(input), to correct this, we return : >> y[(window_len / 2 - 1):-(window_len / 2)] # noqa instead of just y.

meegkit.utils.sig.spectral_envelope(x, sfreq, lowpass=32)#

Compute envelope with convolution.

Notes

The signal is first padded to avoid edge effects. To align the envelope with the input signal, we return : >> y[(window_len / 2 - 1):-(window_len / 2)] # noqa

meegkit.utils.sig.stmcb(x, u_in=None, q=None, p=None, niter=5, a_in=None)#

Compute linear model via Steiglitz-McBride iteration.

[B,A] = stmcb(H,NB,NA) finds the coefficients of the system B(z)/A(z) with approximate impulse response H, NA poles and NB zeros.

[B,A] = stmcb(H,NB,NA,N) uses N iterations. N defaults to 5.

[B,A] = stmcb(H,NB,NA,N,Ai) uses the vector Ai as the initial guess at the denominator coefficients. If you don’t specify Ai, STMCB uses [B,Ai] = PRONY(H,0,NA) as the initial conditions.

[B,A] = STMCB(X,Y,NB,NA,N,Ai) finds the system coefficients B and A of the system which, given X as input, has Y as output. N and Ai are again optional with default values of N = 5, [B,Ai] = PRONY(Y,0,NA). Y and X must be the same length.

Parameters:

x (array)
u_in (array)
q (int)
p (int)
n_iter (int)
a_in (array)

Returns:

b (array) – Filter coefficients (denominator).
a (array) – Filter coefficients (numerator).

Examples

Approximate the impulse response of a Butterworth filter with a system of lower order:

>>> [b, a] = butter(6, 0.2)                # Butterworth filter design
>>> h = filter(b, a, [1, zeros(1,100)])    # Filter data using above filter
>>> freqz(b, a, 128)                       # Frequency response
>>> [bb, aa] = stmcb(h, 4, 4)
>>> plt.plot(freqz(bb, aa, 128))

References

meegkit.utils.sig.teager_kaiser(x, m=1, M=1, axis=0)#

Mean Teager-Kaiser energy operator.

The discrete version of the Teager-Kaiser operator is computed according to:

y[n] = x[n] ** {2 / m} - (x[n - M] * x[n + M]) ** {1 / m}

with m the exponent parameter and M the lag parameter which both are usually equal to 1 for a conventional operator. The Teaser-Kaiser operator can be used to track amplitude modulations (AM) and/or frequency modulations (FM).

Parameters:

x (array, shape=(n_samples[, n_channels][, n_trials])) – Input data.
m (int) – Exponent parameter.
M (int) – Lag parameter.
axis (int) – Axis to compute metric on.

Returns:

Instantaneous energy.

Return type:

array, shape=(n_samples - 2 * M[, n_channels][, n_trials])

References

Adapted from the TKEO function in R library seewave.

Examples

>>> x = np.array([1,  3, 12, 25, 10])
>>> tk_energy = teager_kaiser(x)

Statistics#

Statistics utilities.

meegkit.utils.stats.bootstrap_ci(X, n_bootstrap=2000, ci=(5, 95), axis=-1)#

Confidence intervals from non-parametric bootstrap resampling.

Bootstrap is computed over the chosen axis, with replacement. By default, the resampling is performed over the last dimension (trials), but this can also be done over channels (axis=1).

Parameters:

X (array, shape=(n_times, n_chans[, n_trials])) – Input data.
n_bootstrap (int) – Number of bootstrap resamplings (default=2000). For publication quality results, it is usually recommended to have more than 5000 iterations.
ci (length-2 tuple) – Confidence interval (default=(5, 95)).
axis (int) – Axis to operate on.

Returns:

ci_low, ci_up – Confidence intervals.

Return type:

arrays, shape=(n_times, n_chans)

meegkit.utils.stats.bootstrap_snr(epochs, n_bootstrap=2000, baseline=None, window=None)#

Get bootstrap SNR from non-parametric bootstrap resampling.

The lower bound SNR value can serve as a minimum threshold of signal quality. Parks et al. (2016) suggest using a threshold of 3 dB when assessing whether to remove subjects from ERP studies.

Parameters:

epochs (mne.Epochs instance) – Epochs instance to compute SNR from.
n_bootstrap (int) – Number of bootstrap iterations (should be > 10000 for publication quality).
baseline (tuple or list of length 2, or None) – The time interval to apply rescaling / baseline correction. If None do not apply it. If baseline is (bmin, bmax) the interval is between bmin (s) and bmax (s). If bmin is None the beginning of the data is used and if bmax is None then bmax is set to the end of the interval. If baseline is (None, None) the entire time interval is used. If baseline is None, no correction is applied.
window (tuple or list of length 2, or None) – The time interval used to compute the numerator of the SNR. If window is (bmin, bmax) the interval is between bmin (s) and bmax (s). If bmin is None the beginning of the data is used and if bmax is None then bmax is set to the end of the interval. If window is (None, None) the entire time interval is used. If None use all positive times.

Returns:

ERP (length-3 tuple) – Mean bootstrap ERP, with 90% confidence intervals.
SNR (length-3 tuple) – Mean bootstrap SNR, with 90% confidence intervals.

References

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4751267/

meegkit.utils.stats.cronbach(epochs, K=None, n_bootstrap=2000, tmin=None, tmax=None)#

Compute reliability of ERP.

Internal reliability of the ERN and Pe as a function of increasing trial numbers can be quantified with Cronbach’s alpha:

$$alpha = frac{K}{K-1} left(1-frac{sum_{i=1}^K
sigma^2_{Y_i}}{ sigma^2_X}right)$$

Hinton, Brownlow, McMurray, and Cozens (2004) have suggested that Cronbach’s alpha exceeding .90 indicates excellent internal reliability, between .70 and .90 indicates high internal reliability, from .50 to .70 indicates moderate internal reliability, and below .50 is low.

Parameters:

epochs (mne.Epochs | ndarray, shape=(n_trials, n_chans, n_samples)) – Epochs to compute alpha from.
K (int) – Number of trials to use for alpha computation.
n_bootstrap (int) – Number of bootstrap resamplings.
tmin (float) – Start time of epoch.
tmax (float) – End of epoch.

Returns:

alpha (array, shape=(n_chans,)) – Cronbach alpha value
bounds (length-2 tuple) – Lower and higher bound of CI.

meegkit.utils.stats.rms(X, axis=0)#: Root-mean-square along given axis.

meegkit.utils.stats.robust_mean(X, axis=0, percentile=[5, 95])#: Do robust mean based on JR Kings implementation.

meegkit.utils.stats.rolling_corr(X, y, window=None, sfreq=1, step=1, axis=0)#

Calculate rolling correlation between some data and a reference signal.

Parameters:

X (array, shape=(n_times, n_chans[, n_epochs])) – Test signal. First dimension should be time.
y (array, shape=(n_times[, n_epochs])) – Reference signal.
window (int) – Number of timepoints for to include for each correlation calculation.
sfreq (int) – Sampling frequency (default=1).
step (int) – If > 1, only compute correlations every step samples.

Returns:

corr (array, shape=(n_times - window, n_channels[, n_epochs])) – Rolling window correlation.
t_corr (array) – Corresponding time vector.

meegkit.utils.stats.snr_spectrum(X, freqs, n_avg=1, n_harm=1, skipbins=1)#

Compute Signal-to-Noise-corrected spectrum.

The implementation tries to replicate examples in [1; 2; 3]_.

Parameters:

X (ndarray , shape=(n_freqs, n_chans,[ n_trials,])) – One-sided power spectral density estimate, specified as a real-valued, nonnegative array. The power spectral density must be expressed in linear units, not decibels.
freqs (array, shape=(n_freqs,)) – Frequency bins.
n_avg (int) – Number of neighbour bins to estimate noise over. Make sure that this value doesn’t overlap with neighbouring target frequencies.
n_harm (int) – Compute SNR at each frequency bin as a pooled RMS over this bin and n_harm harmonics (see references below).
skipbins (int) – Number of bins skipped to estimate noise of neighbouring bins.

Returns:

SNR – Signal-to-Noise-corrected spectrum.

Return type:

ndarray, shape=(n_freqs, n_chans, n_trials) or (n_freqs, n_chans)

References

[1]

Cottereau, B. R., McKee, S. P., Ales, J. M., & Norcia, A. M. (2011). Disparity-tuned population responses from human visual cortex. The Journal of Neuroscience, 31(3), 954-965.

[2]

Cottereau, B. R., McKee, S. P., & Norcia, A. M. (2014). Dynamics and cortical distribution of neural responses to 2D and 3D motion in human. Journal of neurophysiology, 111(3), 533-543.

[3]

de Heering, A., & Rossion, B. (2015). Rapid categorization of natural face images in the infant right hemisphere. Elife, 4.

TRCA utilities#

TRCA utils.

meegkit.utils.trca.bandpass(eeg, sfreq, Wp, Ws)#

Filter bank design for decomposing EEG data into sub-band components.

Parameters:

eeg (np.array, shape=(n_samples, n_chans[, n_trials])) – Training data.
sfreq (int) – Sampling frequency of the data.
Wp (2-tuple) – Passband for Chebyshev filter.
Ws (2-tuple) – Stopband for Chebyshev filter.

Returns:

y – Sub-band components decomposed by a filter bank.

Return type:

np.array, shape=(n_trials, n_chans, n_samples)

See also

scipy.signal.cheb1ord: Chebyshev type I filter order selection.

meegkit.utils.trca.itr(n, p, t)#

Compute information transfer rate (ITR).

Definition in [1].

Parameters:

n (int) – Number of targets.
p (float) – Target identification accuracy (0 <= p <= 1).
t (float) – Average time for a selection (s).

Returns:

itr – Information transfer rate [bits/min]

Return type:

float

References

[1]

M. Cheng, X. Gao, S. Gao, and D. Xu, “Design and Implementation of a Brain-Computer Interface With High Transfer Rates”, IEEE Trans. Biomed. Eng. 49, 1181-1186, 2002.

meegkit.utils.trca.normfit(data, ci=0.95)#

Compute the mean, std and confidence interval for them.

Parameters:

data (array, shape=()) – Input data.
ci (float) – Confidence interval (default=0.95).

Returns:

m (float) – Mean.
sigma (float) – Standard deviation
[m - h, m + h] (list) – Confidence interval of the mean.
[sigmaCI_lower, sigmaCI_upper] (list) – Confidence interval of the std.

meegkit.utils.trca.round_half_up(num, decimals=0)#

Round half up round the last decimal of the number.

The rules are: from 0 to 4 rounds down from 5 to 9 rounds up

Parameters:

num (float) – Number to round
decimals (number of decimals)

Return type:

num rounded

meegkit.utils.trca.schaefer_strimmer_cov(X)#

Schaefer-Strimmer covariance estimator.

Shrinkage estimator described in [1]:

\[\hat{\Sigma} = (1 - \gamma)\Sigma_{scm} + \gamma T\]

where $T$ is the diagonal target matrix:

\[T_{i,j} = \{ \Sigma_{scm}^{ii} \text{if} i = j, 0 \text{otherwise} \}\]

Note that the optimal $\gamma$ is estimated by the authors’ method.

Parameters:: X (array, shape=(n_chans, n_samples)) – Signal matrix.
Returns:: cov – Schaefer-Strimmer shrinkage covariance matrix.
Return type:: array, shape=(n_chans, n_chans)

References

[1]

Schafer, J., and K. Strimmer. 2005. A shrinkage approach to large-scale covariance estimation and implications for functional genomics. Statist. Appl. Genet. Mol. Biol. 4:32.