osl_dynamics.analysis.spectral#

Functions to perform spectral analysis.

Functions#

wavelet(data, sampling_frequency[, w, standardize, ...])

Calculate a wavelet transform.

welch_spectra(data, sampling_frequency[, alpha, ...])

Calculates spectra for inferred states using Welch's method.

multitaper_spectra(data, sampling_frequency[, alpha, ...])

Calculates spectra for inferred states using a multitaper.

regression_spectra(data, alpha, sampling_frequency, ...)

Calculates mode-specific spectra by regressing a spectrogram with alpha.

rescale_regression_coefs(psd, alpha, window_length, ...)

Rescales regression coefficients with maximum regressor values.

coherence_spectra(cpsd[, keepdims])

Calculate coherence from cross power spectral densities.

partial_coherence_spectra(cpsd[, keepdims])

Calculate partial coherence from cross power spectral densities.

partial_directed_coherence_spectra(f, cpsd, ...[, ...])

Calculate partial directed coherence from cross power spectra.

decompose_spectra(coherences, n_components[, ...])

Uses non-negative matrix factorization to decompose spectra.

mar_spectra(coeffs, covs, sampling_frequency[, n_freq])

Calculates cross power spectral densities from MAR model parameters.

autocorr_to_spectra(autocorr_func, sampling_frequency)

Calculates spectra from the autocorrelation function.

get_frequency_args_range(frequencies, frequency_range)

Get min/max indices of a range in a frequency axis.

Module Contents#

osl_dynamics.analysis.spectral.wavelet(data, sampling_frequency, w=5, standardize=True, time_range=None, frequency_range=None, df=0.5)[source]#

Calculate a wavelet transform.

This function uses a Morlet2 window to calculate a wavelet transform.

Parameters:

  • data (np.ndarray) – 1D time series data. Shape must be (n_samples,).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • w (float, optional) – Number of oscillations in the Gaussian envelope.

  • standardize (bool, optional) – Should we standardize the data before calculating the wavelet?

  • time_range (list, optional) – Start time and end time to plot in seconds. Default is the full time axis of the data.

  • frequency_range (list of length 2, optional) – Start and end frequency to plot in Hz. Default is [1, sampling_frequency / 2].

  • df (float, optional) – Frequency resolution in Hz.

Returns:

  • t (np.ndarray) – 1D numpy array for the time axis.

  • f (np.ndarray) – 1D numpy array for the frequency axis.

  • wt (np.ndarray) – 2D numpy array (frequency, time) containing the wavelet transform.

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

osl_dynamics.analysis.spectral.welch_spectra(data, sampling_frequency, alpha=None, window_length=None, step_size=None, frequency_range=None, standardize=True, averaging='mean', calc_cpsd=False, calc_coh=True, return_weights=False, keepdims=False, n_jobs=1)[source]#

Calculates spectra for inferred states using Welch’s method.

The scaling for the power spectra calculated by this function matches SciPy (scipy.signal.welch).

Parameters:

  • data (np.ndarray or list) – Time series data. Must have shape (n_sessions, n_samples, n_channels) or (n_samples, n_channels).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • alpha (np.ndarray or list, optional) – Inferred state probability time course. Must have shape (n_sessions, n_samples, n_states) or (n_samples, n_states).

  • window_length (int, optional) – Length of the data segment to use to calculate spectra. If None, we use 2 * sampling_frequency.

  • step_size (int, optional) – Step size for shifting the window. Defaults to window_length // 2 (50% overlap).

  • frequency_range (list, optional) – Minimum and maximum frequency to keep.

  • standardize (bool, optional) – Should we standardize the data before calculating the spectra?

  • averaging (str, optional) – Method used to average periodograms. Must be 'mean' or 'median'.

  • calc_cpsd (bool, optional) – Should we return the cross spectra for psd? If True, we force calc_coh to False.

  • calc_coh (bool, optional) – Should we also return the coherence spectra?

  • return_weights (bool, optional) – Should we return the weights for session-specific spectra? This is useful for calculating a group average.

  • keepdims (bool, optional) – Should we enforce a (n_sessions, n_states, …) array is returned for psd and coh? If False, we remove any dimension of length 1.

  • n_jobs (int, optional) – Number of parallel jobs.

Returns:

  • f (np.ndarray) – Frequencies of the power spectra and coherences. Shape is (n_freq,).

  • psd (np.ndarray) – Power spectra for each session and state. Shape is (n_sessions, n_states, n_channels, n_freq) if calc_cpsd=False, otherwise the shape is (n_sessions, n_states, n_channels, n_channels, n_freq). Any axis of length 1 is removed if keepdims=False.

  • coh (np.ndarray) – Coherences for each session and state. Shape is (n_sessions, n_states, n_channels, n_channels, n_freq). Any axis of length 1 is removed if keepdims=False. Only returned is calc_coh=True.

  • w (np.ndarray) – Weighting for session-specific spectra. Only returned if return_weights=True. Shape is (n_sessions,).

osl_dynamics.analysis.spectral.multitaper_spectra(data, sampling_frequency, alpha=None, window_length=None, step_size=None, time_half_bandwidth=4, n_tapers=7, frequency_range=None, standardize=True, averaging='mean', calc_cpsd=False, calc_coh=True, return_weights=False, keepdims=False, n_jobs=1)[source]#

Calculates spectra for inferred states using a multitaper.

Parameters:

  • data (np.ndarray or list) – Time series data. Must have shape (n_sessions, n_samples, n_channels) or (n_samples, n_channels).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • alpha (np.ndarray or list, optional) – Inferred state probability time course. Must have shape (n_sessions, n_samples, n_states) or (n_samples, n_states).

  • window_length (int, optional) – Length of the data segment to use to calculate spectra.

  • step_size (int, optional) – Step size for shifting the window. Defaults to window_length (no overlap).

  • time_half_bandwidth (float, optional) – Parameter to control the resolution of the spectra.

  • n_tapers (int, optional) – Number of tapers.

  • frequency_range (list, optional) – Minimum and maximum frequency to keep.

  • standardize (bool, optional) – Should we standardize the data before calculating the spectra?

  • averaging (str, optional) – Method used to average periodograms. Must be 'mean' or 'median'.

  • calc_cpsd (bool, optional) – Should we return the cross spectra for psd? If True, we force calc_coh to False.

  • calc_coh (bool, optional) – Should we also return the coherence spectra?

  • return_weights (bool, optional) – Should we return the weights for session-specific spectra? This is useful for calculating a group average.

  • keepdims (bool, optional) – Should we enforce a (n_sessions, n_states, …) array is returned for psd and coh? If False, we remove any dimension of length 1.

  • n_jobs (int, optional) – Number of parallel jobs.

Returns:

  • f (np.ndarray) – Frequencies of the power spectra and coherences. Shape is (n_freq,).

  • psd (np.ndarray) – Power spectra for each session and state. Shape is (n_sessions, n_states, n_channels, n_freq) if calc_cpsd=False, otherwise the shape is (n_sessions, n_states, n_channels, n_channels, n_freq). Any axis of length 1 is removed if keepdims=False.

  • coh (np.ndarray) – Coherences for each session and state. Shape is (n_sessions, n_states, n_channels, n_channels, n_freq). Any axis of length 1 is removed if keepdims=False. Only returned is calc_coh=True.

  • w (np.ndarray) – Weighting for session-specific spectra. Only returned if return_weights=True. Shape is (n_sessions,).

osl_dynamics.analysis.spectral.regression_spectra(data, alpha, sampling_frequency, window_length, step_size, n_sub_windows=1, frequency_range=None, standardize=True, calc_coh=True, return_weights=False, return_coef_int=False, keepdims=False, n_jobs=1)[source]#

Calculates mode-specific spectra by regressing a spectrogram with alpha.

We use osl_dynamics.analysis.spectral._welch_spectrogram() to calculate the spectrogram.

Parameters:

  • data (np.ndarray or list) – Data to calculate a spectrogram for. Shape must be (n_sessions, n_samples, n_channels) or (n_samples, n_channels).

  • alpha (np.ndarray or list) – Inferred mode mixing factors. Shape must be (n_sessions, n_samples, n_modes) or (n_samples, n_modes).

  • sampling_frequency (float) – Sampling_frequency in Hz.

  • window_length (int) – Number samples to use in the window to calculate a spectrum.

  • step_size (int) – Step size for shifting the window.

  • n_sub_windows (int, optional) – We split the window into a number of sub-windows and average the spectra for each sub-window. window_length must be divisible by n_sub_windows.

  • frequency_range (list, optional) – Minimum and maximum frequency to keep.

  • standardize (bool, optional) – Should we standardize the data before calculating the spectrogram?

  • calc_coh (bool, optional) – Should we also return the coherence spectra?

  • return_weights (bool, optional) – Should we return the weights for session-specific spectra? Useful for calculating the group average spectra.

  • return_coef_int (bool, optional) – Should we return the regression coefficients and intercept separately for the mode PSDs?

  • keepdims (bool, optional) – Should we enforce a (n_sessions, n_states, …) array is returned for psd and coh? If False, we remove any dimension of length 1.

  • n_jobs (int, optional) – Number of parallel jobs.

Returns:

  • f (np.ndarray) – Frequency axis. Shape is (n_freq).

  • psd (np.ndarray) – Mode PSDs. A numpy array with shape (n_sessions, 2, n_modes, n_channels, n_freq) where axis=1 is the coefficients/intercept if return_coef_int=True, otherwise shape is (n_sessions, n_modes, n_channels, n_freq). Any axis of length 1 is removed if keepdims=False.

  • coh (np.ndarray) – Mode coherences. Shape is (n_sessions, n_modes, n_channels, n_channels, n_freq). Any axis of length 1 is removed if keepdims=False.

  • w (np.ndarray) – Weight for each session-specific PSD. Shape is (n_sessions,). Only returned if return_weights=True.

osl_dynamics.analysis.spectral.rescale_regression_coefs(psd, alpha, window_length, step_size, n_sub_windows=1)[source]#

Rescales regression coefficients with maximum regressor values.

Parameters:

  • psd (np.ndarray) – Mode PSDs. Shape must be (n_sessions, 2, n_modes, n_channels, n_freq).

  • alpha (np.ndarray or list) – Inferred mode mixing coefficients. Shape must be (n_sessions, n_samples, n_modes) or (n_samples, n_modes).

  • window_length (int) – Number samples used as the window to calculate a PSD.

  • step_size (int) – Step size used for shifting the window.

  • n_sub_windows (int, optional) – Number of sub-windows used to average the spectra. window_length must be divisible by n_sub_windows.

Returns:

psd – Mode PSDs rescaled with maximum regressor values. Shape is identical to the input.

Return type:

np.ndarray

osl_dynamics.analysis.spectral.coherence_spectra(cpsd, keepdims=False)[source]#

Calculate coherence from cross power spectral densities.

We calculate the coherence as:

\[C_{xy}(f) = \frac{|P_{xy}(f)|}{\sqrt{P_{xx}(f) P_{yy}(f)}}\]

where \(x,y\) are different channels.

Parameters:

  • cpsd (np.ndarray) – Cross power spectra. Shape is (…, n_channels, n_channels, n_freq).

  • keepdims (bool, optional) – Should we squeeze any axis of length 1?

Returns:

coh – Coherence spectra. Shape is (…, n_channels, n_channels, n_freq).

Return type:

np.ndarray

osl_dynamics.analysis.spectral.partial_coherence_spectra(cpsd, keepdims=False)[source]#

Calculate partial coherence from cross power spectral densities.

We calculate the partial coherence as:

\[PC_{xy}(f) = \ \frac{|P^{-1}_{xy}(f)|}{\sqrt{P^{-1}_{xx}(f) P^{-1}_{yy}(f)}}\]

where \(x,y\) are different channels and \(P^{-1}(f)\) is the inverse of the cross spectral density matrix.

Parameters:

  • cpsd (np.ndarray) – Cross power spectra. Shape is (…, n_channels, n_channels, n_freq).

  • keepdims (bool, optional) – Should we squeeze any axis of length 1?

Returns:

pcoh – Partial coherence spectra. Shape is (…, n_channels, n_channels, n_freq).

Return type:

np.ndarray

osl_dynamics.analysis.spectral.partial_directed_coherence_spectra(f, cpsd, sampling_frequency, n_iterations=50, tol=1e-07, keepdims=False)[source]#

Calculate partial directed coherence from cross power spectra.

We calculate the partial directed coherence as:

\[PDC_{ij}(f) = \frac{|A_{ij}(f)|}{\sqrt{\sum_m |A_{mj}(f)|^2}}\]

where \(i,j\) are different channels and \(A(f) = H^{-1}(f)\) is the inverse of the transfer function, which is calculated by factorizing the cross spectral density matrix.

By default this function will use all cores available.

Parameters:

  • f (np.ndarray) – Frequency axis. Shape is (n_freq,).

  • cpsd (np.ndarray) – Cross power spectra. Shape is (…, n_channels, n_channels, n_freq).

  • sampling_frequency (float) – Sampling frequency (Hz).

  • n_iterations (int, optional) – Number of iterations in the Wilson factorization of the cross spectra.

  • tol (float, optional) – Tolerance for the Wilson factorization of the cross spectra.

  • keepdims (bool, optional) – Should we squeeze any axis of length 1?

Returns:

pdc – Partial directed coherence spectra. Shape is (…, n_channels, n_channels, n_freq).

Return type:

np.ndarray

osl_dynamics.analysis.spectral.decompose_spectra(coherences, n_components, max_iter=50000, random_state=None, verbose=0)[source]#

Uses non-negative matrix factorization to decompose spectra.

Parameters:

  • coherences (np.ndarray) – Coherences spectra. Shape must be (…, n_channels, n_channels, n_freq,).

  • n_components (int) – Number of spectral components to fit.

  • max_iter (int, optional) – Maximum number of iterations in sklearn.decomposition.non_negative_factorization.

  • random_state (int, optional) – Seed for the random number generator.

  • verbose (int, optional) – Show verbose? (1) yes, (0) no.

Returns:

components – Spectral components. Shape is (n_components, n_freq).

Return type:

np.ndarray

osl_dynamics.analysis.spectral.mar_spectra(coeffs, covs, sampling_frequency, n_freq=512)[source]#

Calculates cross power spectral densities from MAR model parameters.

Parameters:

  • coeffs (np.ndarray) – MAR coefficients. Shape must be (n_modes, n_channels, n_channels) or (n_channels, n_channels) or (n_channels,).

  • covs (np.ndarray) – MAR covariances. Shape must be (n_modes, n_channels, n_channels) or (n_channels, n_channels).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • n_freq (int, optional) – Number of frequency bins in the cross power spectral density to calculate.

Returns:

  • f (np.ndarray) – Frequency axis. Shape is (n_freq,).

  • P (np.ndarray) – Cross power spectral densities. Shape is (n_freq, n_modes, n_channels, n_channels) or (n_freq, n_channels, n_channels).

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

osl_dynamics.analysis.spectral.autocorr_to_spectra(autocorr_func, sampling_frequency, nfft=64, frequency_range=None)[source]#

Calculates spectra from the autocorrelation function.

The power spectrum of each mode is calculated as the Fourier transform of the auto-correlation function. Coherences are calculated from the power spectra.

Parameters:

  • autocorr_func (np.ndarray) – Autocorrelation functions. Shape must be (n_channels, n_channels, n_lags) or (n_modes, n_channels, n_channels, n_lags).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • nfft (int, optional) – Number of data points in the FFT. The auto-correlation function will only have 2 * (n_embeddings + 2) - 1 data points. We pad the auto-correlation function with zeros to have nfft data points if the number of data points in the auto-correlation function is less than nfft.

  • frequency_range (list, optional) – Minimum and maximum frequency to keep (Hz).

Returns:

  • f (np.ndarray) – Frequencies of the power spectra and coherences. Shape is (n_freq,).

  • psd (np.ndarray) – Power (or cross) spectra calculated for each mode. Shape is (n_channels, n_channels, n_freq) or (n_modes, n_channels, n_channels, n_freq).

  • coh (np.ndarray) – Coherences calculated for each mode. Shape is (n_channels, n_channels, n_freq) or (n_modes, n_channels, n_channels, n_freq).

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

osl_dynamics.analysis.spectral.get_frequency_args_range(frequencies, frequency_range)[source]#

Get min/max indices of a range in a frequency axis.

Parameters:

  • frequencies (np.ndarray) – 1D numpy array containing a Frequency axis.

  • frequency_range (list of len 2) – Minimum and maximum frequency.

Returns:

args_range – Indices for minimum and maximum frequency.

Return type:

list of len 2