osl_dynamics.analysis.spectral#

Functions to perform spectral analysis.

Module Contents#

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])

Calculates coherences from cross power spectral densities.

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.

_welch_spectrogram(data, sampling_frequency[, ...])

Calculates a spectogram (time-varying power spectral density)

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

Calculate a (cross) power spectrum using Welch's method.

_multitaper_spectrogram(data, sampling_frequency[, ...])

Calculates a spectogram (time-varying power spectral density)

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

Calculate a multitaper.

_regress_welch_spectrogram(data, alpha, ...)

Calculate regression spectra for a single array.

_window_mean(alpha, window, window_length, step_size, ...)

Applies a windowing function to a time series and takes the mean.

Attributes#

_logger

osl_dynamics.analysis.spectral._logger[source]#
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 scipy.signal.morlet2 window to calculate the 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) –

    w parameter to pass to scipy.signal.morlet2.

  • 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.

osl_dynamics.analysis.spectral.welch_spectra(data, sampling_frequency, alpha=None, window_length=None, step_size=None, frequency_range=None, standardize=True, 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?

  • 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). 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, 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?

  • 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). 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 _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 calcualte 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]#

Calculates coherences from cross power spectral densities.

Parameters:

  • cpsd (np.ndarray) – Cross power spectra. Shape is (n_channels, n_channels, n_freq) or (n_modes, 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.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.decomposion.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).

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).

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

osl_dynamics.analysis.spectral._welch_spectrogram(data, sampling_frequency, window_length=None, step_size=None, frequency_range=None, calc_cpsd=True, n_sub_windows=1)[source]#

Calculates a spectogram (time-varying power spectral density) using Welch’s method.

Steps:

  1. Segment data into overlapping windows.

  2. Multiply each window by a Hann tapering function.

  3. Calculate a periodogram for each window.

We use the same scaling for the power spectra as SciPy (i.e. when you use scipy.signal.periodogram(..., scaling="density")).

Parameters:

  • data (np.ndarray) – Data to calculate the spectrogram for. Shape must be (n_samples, n_channels).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • window_length (int, optional) – Number of data points to use when calculating the periodogram. Defaults to 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.

  • calc_cpsd (bool, optional) – Should we calculate cross spectra?

  • 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.

Returns:

  • t (np.ndarray) – Time axis.

  • f (np.ndarray) – Frequency axis.

  • P (np.ndarray) – Spectrogram. 3D numpy array with shape (time, channels * (channels + 1) / 2, freq) if calc_cpsd=True, otherwise it is (time, channels, freq).

osl_dynamics.analysis.spectral._welch(data, sampling_frequency, alpha=None, window_length=None, step_size=None, frequency_range=None, standardize=True, calc_coh=False, keepdims=False)[source]#

Calculate a (cross) power spectrum using Welch’s method.

This function first calculates a spectrogram using _welch_spectrogram then takes the time average.

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

Parameters:

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

  • sampling_frequency (float) – Sampling frequency in Hz.

  • alpha (np.ndarray, optional) – Inferred state probability time course. Must have shape (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?

  • calc_coh (bool, optional) – Should we calculate the coherence between channels?

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

Returns:

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

  • psd (np.ndarray) – Power spectra. Shape is (…, n_channels, n_freq).

  • coh (np.ndarray) – Coherences spectra. Shape is (…, n_channels, n_channels, n_freq). Only returned is calc_coh=True.

osl_dynamics.analysis.spectral._multitaper_spectrogram(data, sampling_frequency, window_length=None, step_size=None, frequency_range=None, calc_cpsd=True, n_sub_windows=1, time_half_bandwidth=4, n_tapers=7)[source]#

Calculates a spectogram (time-varying power spectral density) using a multitaper.

Steps:

  1. Segment data into overlapping windows.

  2. Multiply each window by a number of tapering functions.

  3. Calculate a periodogram for each tapered window.

  4. Average over tapers.

Parameters:

  • data (np.ndarray) – Data to calculate the spectrogram for. Shape must be (n_samples, n_channels).

  • sampling_frequency (float) – Sampling frequency in Hz.

  • window_length (int, optional) – Number of data points to use when calculating the periodogram. Defaults to 2 * sampling_frequency.

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

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

  • calc_cpsd (bool, optional) – Should we calculate cross spectra?

  • 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.

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

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

Returns:

  • t (np.ndarray) – Time axis.

  • f (np.ndarray) – Frequency axis.

  • P (np.ndarray) – Spectrogram. 3D numpy array with shape (time, channels * (channels + 1) / 2, freq) if calc_cpsd=True, otherwise it is (time, channels, freq).

osl_dynamics.analysis.spectral._multitaper(data, sampling_frequency, alpha=None, window_length=None, step_size=None, time_half_bandwidth=4, n_tapers=7, frequency_range=None, standardize=True, calc_coh=True, keepdims=False)[source]#

Calculate a multitaper.

This function first calculates a spectrogram using _multitaper_spectrogram then takes the time average.

Parameters:

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

  • sampling_frequency (float) – Sampling frequency in Hz.

  • alpha (np.ndarray, optional) – Inferred state probability time course. Must have shape (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 (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?

  • calc_coh (bool, optional) – Should we calculate the coherence between channels?

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

Returns:

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

  • psd (np.ndarray) – Power spectra. Shape is (…, n_channels, n_freq).

  • coh (np.ndarray) – Coherences spectra. Shape is (…, n_channels, n_channels, n_freq). Only returned is calc_coh=True.

osl_dynamics.analysis.spectral._regress_welch_spectrogram(data, alpha, sampling_frequency, window_length, step_size, n_sub_windows, frequency_range, calc_cpsd)[source]#

Calculate regression spectra for a single array.

Parameters:

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

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

  • sampling_frequency (float) – Sampling_frequency in Hz.

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

  • 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.

  • calc_cpsd (bool, optional) – Should we calculate the cross spectra?

Returns:

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

  • coefs (np.ndarray) – Regression coefficients. Shape is (n_modes, n_channels * (n_channels + 1) / 2, n_freq) if calc_cpsd=True, otherwise it is (n_modes, n_channels, n_freq).

  • intercept (np.ndarray) – Intercept. Shape is (n_channels * (n_channels + 1) / 2, n_freq) if calc_cpsd=True, otherwise it is (n_channels, n_freq).

osl_dynamics.analysis.spectral._window_mean(alpha, window, window_length, step_size, n_sub_windows)[source]#

Applies a windowing function to a time series and takes the mean.

Parameters:

  • alpha (np.ndarray) – Time series data. Shape is (n_samples, n_modes).

  • window (str) – Name of the windowing function.

  • window_length (int) – Number of data points in a window.

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

  • n_sub_windows (int) – We split the window into a number of sub-windows and average across all sub-windows. This is the number of sub-windows.

Returns:

mean_window_alpha – Mean for each window.

Return type:

np.ndarray