spectral_connectivity.connectivity.Connectivity

class Connectivity(fourier_coefficients: ~numpy.ndarray[tuple[int, ...], ~numpy.dtype[~numpy.complexfloating]], expectation_type: str = 'trials_tapers', frequencies: ~numpy.ndarray[tuple[int, ...], ~numpy.dtype[~numpy.floating]] | None = None, time: ~numpy.ndarray[tuple[int, ...], ~numpy.dtype[~numpy.floating]] | None = None, dtype: ~numpy.dtype[~typing.Any] | None | type[~typing.Any] | ~numpy._typing._dtype_like._SupportsDType[~numpy.dtype[~typing.Any]] | str | tuple[~typing.Any, int] | tuple[~typing.Any, ~typing.SupportsIndex | ~collections.abc.Sequence[~typing.SupportsIndex]] | list[~typing.Any] | ~numpy._typing._dtype_like._DTypeDict | tuple[~typing.Any, ~typing.Any] = <class 'numpy.complex128'>, minimum_phase_tolerance: float = 1e-08, minimum_phase_max_iterations: int = 500, is_one_sided: bool = False, observation_weights: ~numpy.ndarray[tuple[int, ...], ~numpy.dtype[~numpy.floating]] | None = None, observations_are_independent: bool = True, time_bins_are_independent: bool = True, *, _adopt_fourier_coefficients: bool = False)

Bases: object

Compute functional and directed connectivity measures from spectral data.

This class provides a comprehensive suite of connectivity analysis methods based on cross-spectral matrices derived from Fourier-transformed time series. Methods range from basic coherence to advanced Granger causality measures.

Parameters:

• fourier_coefficients (NDArray[complexfloating], shape (n_time_windows, n_trials, n_tapers, n_frequencies, n_signals)) – Complex-valued Fourier coefficients from spectral analysis. Must be two-sided (positive and negative frequencies) for Granger methods. Usually obtained from multitaper or other spectral estimation methods. Validation: Must be 5-dimensional with at least 2 signals and contain only finite values (no NaN/Inf).

• expectation_type ({"trials_tapers", "trials", "tapers", "time",) – “time_trials”, “time_tapers”, “time_trials_tapers”}, default=”trials_tapers” Specifies how to average the cross-spectral matrix: - “trials_tapers”: average over trials and tapers (most common) - “trials”: average over trials only (keep taper dimension) - “tapers”: average over tapers only (keep trial dimension) - “time”: average over time windows - combinations: average over multiple dimensions

• frequencies (NDArray[floating], shape (n_frequencies,), optional) – Frequency values in Hz corresponding to FFT bins. If None, uses normalized frequencies.

• time (NDArray[floating], shape (n_time_windows,), optional) – Time values in seconds for each time window. If None, uses indices.

• dtype (np.dtype, default=complex128) – Data type for internal computations. Should match input precision.

• minimum_phase_tolerance (float, default=1e-8) – Relative convergence tolerance for the Wilson minimum-phase factorization used by the directed measures (spectral Granger, DTF, PDC, and relatives).

• minimum_phase_max_iterations (int, default=500) – Maximum Wilson iterations. Near-singular cross-spectral matrices (highly correlated channels) can need several hundred iterations; if the directed measures return NaN with a non-convergence warning, increase this value. The factorization returns early once every sub-spectrum has converged, so a large ceiling is cheap for well-conditioned data.

• is_one_sided (bool, default=False) – Whether coefficients contain only non-negative frequencies. One-sided transforms are returned without FFT half-spectrum slicing or power doubling and cannot be used by Wilson-factorized directed measures.

• observation_weights (ndarray, optional) – Finite non-negative weights with shape (time, trial, taper, frequency, 1). They are applied to every expectation and shared across signals. Transform constructors supply these automatically when smoothing uses a non-uniform kernel or masks invalid edge estimates.

• observations_are_independent (bool, default=True) – Whether the trial/taper observations are statistically independent. Transform constructors set this False when the observation axis holds correlated estimates (a MorletWavelet smoothing neighborhood, or Welch segments overlapping by more than half). Measures whose finite-sample corrections or null distributions count observations warn in that case, and jackknife refuses to leave out tapers.

• time_bins_are_independent (bool, default=True) – Whether the time bins may be counted as independent observations when the expectation averages over time. Transform constructors set this False when successive time bins are correlated (Multitaper or ShortTimeFourierTransform windows overlapping by more than half, or MorletWavelet samples closer than four wavelet standard deviations). The measures that count observations then warn for expectations that include time; it has no effect otherwise.

n_observations

Number of trial/taper observations reduced by the expectation. This is the raw count, not a weighted effective sample size, and it is not reduced for correlated observations (see observations_are_independent).

Type:

int

See also

spectral_connectivity.transforms.Multitaper

Produce the Fourier coefficients this class consumes.

spectral_connectivity.wrapper.multitaper_connectivity

High-level interface returning labeled xarray results with explicit source/target axes.

Notes

Array orientation: pairwise measures end in an (n_signals, n_signals) pair ((n_groups, n_groups) for group measures, ordered as the returned labels), indexed source first. For the directed measures (the spectral Granger family, directed transfer function, directed coherence, (generalized) partial directed coherence, and direct directed transfer function) result[..., i, j] is the influence of signal i on signal j (i -> j). For the lead/lag measures (directed_phase_lag_index, phase_slope_index, group_delay, delay) and the antisymmetric phase measures (coherence_phase, imaginary_coherency, phase_lag_index, weighted_phase_lag_index) positive result[..., i, j] (above 0.5 for directed_phase_lag_index) means signal i leads signal j.

The labeled wrappers multitaper_connectivity() and fourier_connectivity() use the same order: result.sel(source=a, target=b) is a -> b (or “a leads b”). Prefer them unless you need this lower-level API.

Intermediates shared across measures (the expected cross-spectral matrix, power, phase-lag moments, and the minimum-phase factor with the transfer function, noise covariance, and MVAR coefficients derived from it) are cached on first access. Reassigning fourier_coefficients, expectation_type, or observation_weights automatically invalidates these caches, so reusing an instance for new data is safe (constructing a new instance is still the clearer choice). Call clear_cache() to release them while keeping the instance.

The class supports both CPU (NumPy) and GPU (CuPy) computation depending on the SPECTRAL_CONNECTIVITY_ENABLE_GPU environment variable. For Granger causality measures, minimum phase decomposition [1] is used to estimate transfer functions and noise covariances non-parametrically.

References

[1]

Dhamala, M., Rangarajan, G., and Ding, M. (2008). Analyzing information flow in brain networks with nonparametric Granger causality. NeuroImage 41, 354-362.

[2]

Bastos, A. M., & Schoffelen, J. M. (2016). A tutorial review of functional connectivity analysis methods and their interpretational pitfalls. Frontiers in systems neuroscience, 9, 175.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity
>>> rng = np.random.default_rng(0)
>>> n_times, n_trials, n_tapers, n_freqs, n_signals = 50, 10, 5, 100, 2
>>> # Create complex coefficients with coherence injected at frequency bin 10
>>> phase_diff = np.pi / 4  # 45 degree phase difference
>>> coeffs = (
...     rng.standard_normal((n_times, n_trials, n_tapers, n_freqs, n_signals))
...     + 1j
...     * rng.standard_normal((n_times, n_trials, n_tapers, n_freqs, n_signals))
... )
>>> coeffs[:, :, :, 10, 1] = coeffs[:, :, :, 10, 0] * np.exp(1j * phase_diff)
>>> conn = Connectivity(coeffs, expectation_type="trials_tapers")
>>> coherence = conn.coherence_magnitude()
>>> coherence.shape  # (n_times, non-negative freqs, n_signals, n_signals)
(50, 51, 2, 2)
>>> print(f"Peak coherence: {np.max(coherence[:, 10, 0, 1]):.3f}")
Peak coherence: 1.000

Attributes:

all_frequencies

Return positive and negative frequencies of the transform.

expectation_type

Which dimensions the cross-spectral matrix is averaged over.

fourier_coefficients

Multitaper Fourier coefficients.

frequencies

Return non-negative frequencies of the transform.

is_one_sided

Whether the input contains only non-negative frequencies.

n_observations

Return the raw number of observations averaged by the expectation.

n_signals

Number of signals represented by the Fourier coefficients.

observation_weights

Non-negative weights used when averaging spectral observations.

observations_are_independent

Whether the trial/taper observations are statistically independent.

time_bins_are_independent

Whether the time bins may be counted as independent observations.

Methods

blockwise_spectral_granger_prediction(...)	Return spectral Granger prediction between multichannel groups.
canonical_coherence(group_labels)	Return the historical magnitude-squared canonical correlation.
canonical_coherency(group_labels, *[, rank, ...])	Return exact complex canonical coherency (CaCoh) components.
clear_cache()	Free the intermediates cached for reuse across measures.
coherence_magnitude()	Return the magnitude squared of the complex coherency.
coherence_phase()	Return the phase angle of the complex coherency.
coherency()	Return the complex-valued linear association between time series.
conditional_spectral_granger_prediction()	Return pairwise spectral Granger prediction conditioned on all others.
corrected_imaginary_phase_locking_value()	Return corrected imaginary phase-locking value (ciPLV).
cross_spectral_density()	Return the one-sided cross-spectral density matrix.
debiased_squared_phase_lag_index()	Return square of phase lag index corrected for positive bias.
debiased_squared_weighted_phase_lag_index()	Return square of weighted phase lag index corrected for bias.
delay([frequencies_of_interest, ...])	Find a range of possible delays from the coherence phase.
direct_directed_transfer_function()	Return the direct directed transfer function (dDTF).
directed_coherence()	Return the squared directed coherence (noise-weighted DTF).
directed_phase_lag_index()	Return the directed phase-lag index (dPLI).
directed_transfer_function()	Return transfer function coupling strength normalized by inflow.
from_multitaper(multitaper_instance[, ...])	Construct from a spectral transform; the original name of from_transform.
from_transform(transform[, ...])	Construct from any spectral transform.
generalized_partial_directed_coherence()	Return generalized partial directed coherence.
global_coherence([max_rank, ...])	Find linear combinations that capture the most coherent power.
group_delay([frequencies_of_interest, ...])	Return the average time-delay of a broadband signal.
imaginary_coherence()	Return the normalized imaginary component of the cross-spectrum.
imaginary_coherency()	Return the signed imaginary component of coherency.
jackknife(method, *[, confidence_level, ...])	Estimate uncertainty by leaving out one trial/taper observation.
maximized_imaginary_coherency(group_labels)	Return maximized imaginary coherency (MIC) between signal groups.
maximized_imaginary_coherency_components(...)	Return component-resolved MIC scores, filters, and patterns.
minimum_phase_reconstruction_error()	Return the relative reconstruction error of the Wilson factorization.
multivariate_interaction_measure(group_labels)	Return the multivariate interaction measure (MIM) between groups.
pairwise_phase_consistency()	Return square of phase locking value corrected for bias.
pairwise_spectral_granger_prediction()	Return amount of power at a node explained by other nodes.
partial_coherence([regularization])	Return magnitude-squared coherence conditional on all other signals.
partial_directed_coherence()	Return transfer function coupling strength normalized by outflow.
phase_lag_index()	Return non-parametric synchrony measure mitigating power differences.
phase_locking_value()	Return the cross-spectrum with power scaled to magnitude 1.
phase_slope_index([frequencies_of_interest, ...])	Return weighted average of slopes projected onto imaginary axis.
power()	Return the one-sided power spectral density of the signal.
subset_pairwise_spectral_granger_prediction(pairs)	Return predictive power for a subset of signal pairs.
time_reversed_spectral_granger_prediction()	Return pairwise spectral Granger prediction after time reversal.
weighted_phase_lag_index()	Return weighted average of phase lag index using imaginary coherency magnitudes.

Methods

blockwise_spectral_granger_prediction	Return spectral Granger prediction between multichannel groups.
canonical_coherence	Return the historical magnitude-squared canonical correlation.
canonical_coherency	Return exact complex canonical coherency (CaCoh) components.
clear_cache	Free the intermediates cached for reuse across measures.
coherence_magnitude	Return the magnitude squared of the complex coherency.
coherence_phase	Return the phase angle of the complex coherency.
coherency	Return the complex-valued linear association between time series.
conditional_spectral_granger_prediction	Return pairwise spectral Granger prediction conditioned on all others.
corrected_imaginary_phase_locking_value	Return corrected imaginary phase-locking value (ciPLV).
cross_spectral_density	Return the one-sided cross-spectral density matrix.
debiased_squared_phase_lag_index	Return square of phase lag index corrected for positive bias.
debiased_squared_weighted_phase_lag_index	Return square of weighted phase lag index corrected for bias.
delay	Find a range of possible delays from the coherence phase.
direct_directed_transfer_function	Return the direct directed transfer function (dDTF).
directed_coherence	Return the squared directed coherence (noise-weighted DTF).
directed_phase_lag_index	Return the directed phase-lag index (dPLI).
directed_transfer_function	Return transfer function coupling strength normalized by inflow.
from_multitaper	Construct from a spectral transform; the original name of from_transform.
from_transform	Construct from any spectral transform.
generalized_partial_directed_coherence	Return generalized partial directed coherence.
global_coherence	Find linear combinations that capture the most coherent power.
group_delay	Return the average time-delay of a broadband signal.
imaginary_coherence	Return the normalized imaginary component of the cross-spectrum.
imaginary_coherency	Return the signed imaginary component of coherency.
jackknife	Estimate uncertainty by leaving out one trial/taper observation.
maximized_imaginary_coherency	Return maximized imaginary coherency (MIC) between signal groups.
maximized_imaginary_coherency_components	Return component-resolved MIC scores, filters, and patterns.
minimum_phase_reconstruction_error	Return the relative reconstruction error of the Wilson factorization.
multivariate_interaction_measure	Return the multivariate interaction measure (MIM) between groups.
pairwise_phase_consistency	Return square of phase locking value corrected for bias.
pairwise_spectral_granger_prediction	Return amount of power at a node explained by other nodes.
partial_coherence	Return magnitude-squared coherence conditional on all other signals.
partial_directed_coherence	Return transfer function coupling strength normalized by outflow.
phase_lag_index	Return non-parametric synchrony measure mitigating power differences.
phase_locking_value	Return the cross-spectrum with power scaled to magnitude 1.
phase_slope_index	Return weighted average of slopes projected onto imaginary axis.
power	Return the one-sided power spectral density of the signal.
subset_pairwise_spectral_granger_prediction	Return predictive power for a subset of signal pairs.
time_reversed_spectral_granger_prediction	Return pairwise spectral Granger prediction after time reversal.
weighted_phase_lag_index	Return weighted average of phase lag index using imaginary coherency magnitudes.

Attributes

all_frequencies	Return positive and negative frequencies of the transform.
expectation_type	Which dimensions the cross-spectral matrix is averaged over.
fourier_coefficients	Multitaper Fourier coefficients.
frequencies	Return non-negative frequencies of the transform.
is_one_sided	Whether the input contains only non-negative frequencies.
n_observations	Return the raw number of observations averaged by the expectation.
n_signals	Number of signals represented by the Fourier coefficients.
observation_weights	Non-negative weights used when averaging spectral observations.
observations_are_independent	Whether the trial/taper observations are statistically independent.
time_bins_are_independent	Whether the time bins may be counted as independent observations.

property all_frequencies: ndarray[tuple[int, ...], dtype[floating]]

Return positive and negative frequencies of the transform.

Returns:

All frequency values including negative frequencies.

Return type:

NDArray[floating], shape (n_frequencies,)

blockwise_spectral_granger_prediction(group_labels: ndarray[tuple[int, ...], dtype[integer]]) → tuple[ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[integer]]]

Return spectral Granger prediction between multichannel groups.

Parameters:

group_labels (array-like, shape (n_signals,)) – Label assigning each signal to one non-overlapping group.

Returns:

• blockwise_granger (array) – Shape (..., n_nonnegative_frequencies, n_groups, n_groups). Output [..., i, j] is the influence of group i on group j (i -> j), with groups ordered as in labels. The diagonal is NaN.

• labels (array, shape (n_groups,)) – Sorted unique group labels.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Group "a": two noisy copies of signal 0; "b": two of its 6 ms-delayed copy.
>>> pair = np.stack([leader[3:], leader[:-3]], axis=-1)
>>> signals = np.repeat(pair, 2, axis=-1)  # (time, trials, 4 signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> granger, labels = connectivity.blockwise_spectral_granger_prediction(
...     ["a", "a", "b", "b"]
... )
>>> granger.shape  # (n_time_windows, n_frequencies, n_groups, n_groups)
(1, 501, 2, 2)
>>> labels
array(['a', 'b'], dtype='<U1')
>>> # [..., i, j] is group i -> group j, so a -> b is [..., 0, 1] (10 Hz = bin 20).
>>> bool(granger[0, 20, 0, 1] > 10 * granger[0, 20, 1, 0])
True

canonical_coherence(group_labels: ndarray[tuple[int, ...], dtype[integer]]) → tuple[ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[integer]]]

Return the historical magnitude-squared canonical correlation.

The canonical coherence finds two sets of weights such that the coherence between the linear combination of group1 and the linear combination of group2 is maximized.

Parameters:

group_labels (array-like, shape (n_signals,)) – Links each signal to a group.

Returns:

• canonical_coherence (array) – Shape (n_time_windows, n_nonnegative_frequencies, n_groups, n_groups). The maximal coherence for each group pair (symmetric; NaN diagonal).

• labels (array, shape (n_groups,)) – The sorted unique group labels that correspond to n_groups.

Notes

Range: [0, 1]. Maximal coherence values are bounded like coherence magnitude.

Trials x tapers are the observations. A group pair with more signals than observations has intersecting observation subspaces, which forces its value to 1 for any data; this case warns.

References

[1]

Stephen, E.P. (2015). Characterizing dynamically evolving functional networks in humans with application to speech. Boston University.

See also

canonical_coherency

Exact complex, phase-optimised Vidaurre CaCoh with component filters and patterns.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Group "a": two noisy copies of signal 0; "b": two of its 6 ms-delayed copy.
>>> pair = np.stack([leader[3:], leader[:-3]], axis=-1)
>>> signals = np.repeat(pair, 2, axis=-1)  # (time, trials, 4 signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> coherence, labels = connectivity.canonical_coherence(["a", "a", "b", "b"])
>>> coherence.shape  # (n_time_windows, n_frequencies, n_groups, n_groups)
(1, 501, 2, 2)
>>> labels
array(['a', 'b'], dtype='<U1')
>>> bool(coherence[0, 20, 0, 1] > 0.5)  # strong a-b coupling at 10 Hz (bin 20)
True

canonical_coherency(group_labels: ndarray[tuple[int, ...], dtype[Any]], *, rank: int | None = None, n_components: int = 1, regularization: float = 1e-12) → MultivariateConnectivityResult

Return exact complex canonical coherency (CaCoh) components.

This implements Vidaurre et al.’s phase-optimised CaCoh definition: for each candidate phase, the real projection of the between-group CSD is whitened by the real within-group CSDs, and the phase giving the largest singular value is selected. Scores are complex, with magnitude equal to the maximised coherence and phase encoded using MNE’s magnitude * exp(-1j * phi) convention.

Unlike the historical canonical_coherence(), this method returns component-resolved scores, spatial filters, and Haufe-style patterns. Additional components are extracted by CCA-style deflation in the whitened space: each is sought in the orthogonal complement of the previous components’ whitened directions, so successive component signals are uncorrelated within each group (a_j^T Re(Caa) a_k = 0 for j != k) and every component is invariant to invertible real mixing of a group’s channels.

Parameters:

• group_labels (array-like, shape (n_signals,)) – Label assigning each signal to a group; every unordered pair of groups is one connection.

• rank (int, optional) – Retain at most this many within-group whitening directions per group. None keeps every numerically non-zero direction. Applied identically to both groups of every connection.

• n_components (int, default=1) – Number of coherency components to return. A connection whose smaller group has fewer channels returns NaN for the unavailable components.

• regularization (float, default=1e-12) – Relative diagonal loading used by the whitening decomposition.

Returns:

Complex scores of shape (..., frequency, connection, component) plus real filters and patterns; see the class docstring.

Return type:

MultivariateConnectivityResult

Notes

Range: score magnitudes lie in [0, 1]. The whitening, singular-value decomposition, and phase optimization (every local maximum of a 74-point coarse phase grid is refined by a batched Newton iteration and the best is kept, so near-equal lobes of the phase objective are resolved unless two lie within about 2 * pi / 74 of each other) are vectorized over the time/frequency axes on the active xp backend, so this runs on the GPU when GPU support is enabled.

Phase convention: a spatial filter and its negative span the same direction, so the canonical phase is intrinsically defined only modulo pi. Each filter’s sign is fixed so that the largest-magnitude coefficient of its spatial pattern is positive; with that convention the score is the conjugate of the canonical coherency between the positively oriented components, and for single-channel groups it reduces exactly to the conjugate pairwise coherency.

Relation to mne-connectivity: the first component maximizes the same objective as mne-connectivity’s cacoh and matches it (magnitude, and phase modulo pi) wherever both optimizers reach the global maximum. Components >= 2 differ from mne-connectivity 0.9, which deflates the cross-spectrum with the previous components’ channel-space filters; the whitened-space deflation used here keeps them uncorrelated within each group and invariant to invertible real within-group mixing.

References

[1]

Vidaurre C, et al. (2019) Canonical maximization of coherence: A novel tool for investigation of neuronal interactions between two datasets. NeuroImage 201:116009.

[2]

Haufe S, et al. (2014) On the interpretation of weight vectors of linear models in multivariate neuroimaging. NeuroImage 87:96-110.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Group "a": two noisy copies of signal 0; "b": two of its 6 ms-delayed copy.
>>> pair = np.stack([leader[3:], leader[:-3]], axis=-1)
>>> signals = np.repeat(pair, 2, axis=-1)  # (time, trials, 4 signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> result = connectivity.canonical_coherency(["a", "a", "b", "b"])
>>> result.scores.shape  # (n_time_windows, n_frequencies, n_connections, n_components)
(1, 501, 1, 1)
>>> result.connections  # one row per group pair
array([['a', 'b']], dtype='<U1')
>>> result.filters.shape  # (..., n_connections, n_components, side, n_signals)
(1, 501, 1, 1, 2, 4)
>>> bool(abs(result.scores[0, 20, 0, 0]) > 0.5)  # |CaCoh| at 10 Hz (bin 20)
True

clear_cache() → None

Free the intermediates cached for reuse across measures.

Measures computed on one instance share intermediates such as the expected cross-spectral matrix and the minimum-phase factorization, which can each take n_frequencies * n_signals**2 values for every observation that expectation_type leaves unaveraged (each time window by default). Call this after the last measure that needs them to release the memory while keeping the instance; later measures recompute them and return the same results. Replacing fourier_coefficients, expectation_type, or observation_weights clears the cache automatically.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity
>>> rng = np.random.default_rng(0)
>>> fourier_coefficients = rng.standard_normal((1, 5, 3, 16, 4)) + 0j
>>> connectivity = Connectivity(fourier_coefficients)
>>> coherence = connectivity.coherence_magnitude()
>>> connectivity.clear_cache()
>>> bool(np.array_equal(connectivity.coherence_magnitude(), coherence, equal_nan=True))
True

coherence_magnitude() → ndarray[tuple[int, ...], dtype[floating]]

Return the magnitude squared of the complex coherency.

Note that the squared modulus of coherency (originally a complex quantity) is the magnitude-squared coherence (i.e., the normalized, real component of coherency). This value should be bounded by 0 and 1.

Returns:

magnitude – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Magnitude-squared coherence values (symmetric in the signal pair).

Return type:

array

Notes

Range: [0, 1]. Implementation may produce tiny numerical excursions beyond bounds due to floating-point precision.

References

[1]

Hansson-Sandsten M (2011) Cross-spectrum and coherence function estimation using time-delayed Thomson multitapers. In: 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp 4240-4243.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> coherence = connectivity.coherence_magnitude()
>>> coherence.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> round(float(coherence[0, 20, 0, 1]), 1)  # strong coupling at 10 Hz (bin 20)
0.6

coherence_phase() → ndarray[tuple[int, ...], dtype[floating]]

Return the phase angle of the complex coherency.

Returns:

phase – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Phase angles in radians. Positive [..., i, j] means signal i leads signal j; the result is antisymmetric in the signal pair.

Return type:

array

Notes

Range: [-π, π]. Phase angles in radians for complex coherency. A pure delay tau (seconds) gives a phase of 2 * pi * f * tau (wrapped into [-π, π]). canonical_coherency() reports its phase in the conjugate convention (magnitude * exp(-1j * phi)), so for single-channel groups its angle is the negative of this one.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> phase = connectivity.coherence_phase()
>>> phase.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # Signal 0 leads, so [..., 0, 1] is positive: about 2 * pi * 10 Hz * 6 ms at 10 Hz.
>>> round(float(phase[0, 20, 0, 1]), 1)  # bin 20 is 10 Hz
0.4
>>> round(float(phase[0, 20, 1, 0]), 1)
-0.4

coherency() → ndarray[tuple[int, ...], dtype[complexfloating]]

Return the complex-valued linear association between time series.

Computed in the frequency domain.

Returns:

complex_coherency – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Complex coherency between all signal pairs; Hermitian in the signal pair, with its angle given by coherence_phase().

Return type:

array

Notes

Range: Magnitude \(|C_{xy}(f)|\) is in [0, 1]; phase is in [-π, π]. Values lie in the unit disk of the complex plane.

Phase convention: [..., i, j] is S_ij / sqrt(S_ii S_jj) with S_ij = E[X_i conj(X_j)], so a positive angle means signal i leads signal j. canonical_coherency() uses the conjugate convention (magnitude * exp(-1j * phi)): with single-channel groups its score is conj of this coherency.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> coherency = connectivity.coherency()
>>> coherency.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # At 10 Hz (bin 20) the magnitude is the coupling strength and a positive
>>> # angle at [..., 0, 1] means signal 0 leads signal 1.
>>> round(float(abs(coherency[0, 20, 0, 1])), 1)
0.8
>>> bool(np.angle(coherency[0, 20, 0, 1]) > 0)
True

conditional_spectral_granger_prediction() → ndarray[tuple[int, ...], dtype[floating]]

Return pairwise spectral Granger prediction conditioned on all others.

For each ordered source-target pair, the influence of the source on the target is measured after accounting for every remaining signal, using the frequency-domain conditional measure of Chen, Bressler and Ding (2006): the full model containing every signal and the reduced model omitting the source are each spectrally factorized, and the reduced model’s innovation spectrum for the target is split into the part explained by the target’s own full-model innovations and the remainder attributable to the source. With two signals this reduces to ordinary pairwise spectral Granger.

Returns:

conditional_granger – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Output [..., i, j] is the influence of signal i on signal j (i -> j), conditional on every signal other than i and j. The diagonal is NaN.

Return type:

array

Notes

Non-negativity: spectral Granger is >= 0 by definition. Negative estimates within roundoff of zero – above -100 * eps of the result dtype on the log-ratio scale, about -2e-14 for float64 – are clipped to 0; materially negative bins below that threshold (a degenerate factorization) are returned as NaN; use minimum_phase_reconstruction_error() to diagnose them. Other packages (FieldTrip, MVGC, mne-connectivity) return such values as-is.

Range: [0, ∞). The measure is a log-ratio of a total to an intrinsic innovation spectrum, so it is non-negative up to roundoff; bins where either spectrum is not positive (a degenerate factorization) are returned as NaN with a warning.

Cost: n_signals + 1 minimum-phase factorizations (the full system once, plus one (n_signals - 1)-channel system per source), each shared by every target. The full-system factorization is cached and shared with the other directed measures.

References

[1]

Chen, Y., Bressler, S.L., and Ding, M. (2006). Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. Journal of Neuroscience Methods 150, 228-237.

[2]

Geweke, J.F. (1984). Measures of conditional linear dependence and feedback between time series. Journal of the American Statistical Association 79, 907-915.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> # Chain 0 -> 1 -> 2: each signal is the previous one delayed 3 samples, plus noise.
>>> source = rng.standard_normal((1006, 20))
>>> middle = source[:-3] + 0.5 * rng.standard_normal((1003, 20))
>>> target = middle[:-3] + 0.5 * rng.standard_normal((1000, 20))
>>> signals = np.stack([source[6:], middle[3:], target], axis=-1)  # (time, trials, 3)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> conditional = connectivity.conditional_spectral_granger_prediction()
>>> conditional.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 3, 3)
>>> # [..., i, j] is i -> j. Pairwise Granger sees the indirect 0 -> 2 ([..., 0, 2]),
>>> # but conditioning on signal 1 removes it while keeping 1 -> 2 (10 Hz = bin 20).
>>> pairwise = connectivity.pairwise_spectral_granger_prediction()
>>> bool(pairwise[0, 20, 0, 2] > 0.5), bool(conditional[0, 20, 0, 2] < 0.05)
(True, True)
>>> bool(conditional[0, 20, 1, 2] > 0.5)
True

corrected_imaginary_phase_locking_value() → ndarray[tuple[int, ...], dtype[floating]]

Return corrected imaginary phase-locking value (ciPLV).

ciPLV removes the contribution of zero- and pi-lag phase locking while correcting the imaginary PLV for the reduction in its attainable range.

Returns:

corrected_imaginary_phase_locking_value – Shape (..., n_nonnegative_frequencies, n_signals, n_signals).

Return type:

array

Notes

Range: [0, 1]. Exact zero- or pi-lag locking has a zero numerator and denominator and is defined as zero.

References

[1]

Bruña, R., Maestú, F., and Pereda, E. (2018). Phase locking value revisited: teaching new tricks to an old dog. Journal of Neural Engineering 15, 056011.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> ciplv = connectivity.corrected_imaginary_phase_locking_value()
>>> ciplv.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # The 6 ms lag is not zero-phase, so lagged locking remains (10 Hz = bin 20).
>>> bool(ciplv[0, 20, 0, 1] > 0.2)
True

cross_spectral_density() → ndarray[tuple[int, ...], dtype[complexfloating]]

Return the one-sided cross-spectral density matrix.

The diagonal contains the one-sided power spectral densities returned by power(); off-diagonal entries retain both the amplitude and relative-phase information between signal pairs. Interior positive frequency bins are doubled so that the one-sided result has the same total power as the two-sided spectrum. DC and, for an even FFT length, Nyquist are not doubled.

Returns:

cross_spectral_density – Shape (..., n_nonnegative_frequencies, n_signals, n_signals).

Return type:

array

Notes

The matrix is Hermitian at every time-frequency bin and has physical units of signal squared per Hz when the input signal has physical units. Unlike connectivity measures normalized to [0, 1], its magnitude has no finite upper bound.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> csd = connectivity.cross_spectral_density()
>>> csd.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # The diagonal is the power spectrum; the matrix is Hermitian.
>>> bool(np.allclose(csd[..., 0, 0].real, connectivity.power()[..., 0]))
True
>>> bool(np.allclose(csd[..., 0, 1], np.conj(csd[..., 1, 0])))
True

debiased_squared_phase_lag_index() → ndarray[tuple[int, ...], dtype[floating]]

Return square of phase lag index corrected for positive bias.

The square of the phase lag index corrected for the positive bias induced by using the magnitude of the complex cross-spectrum.

Returns:

phase_lag_index – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Debiased squared phase lag index values (symmetric in the signal pair, so they carry no lead/lag direction).

Return type:

array

Notes

Range: [-1 / (n_observations - 1), 1]. The unbiased finite-sample estimate can be negative when the observed phase consistency is below its null bias; negative values do not represent negative coupling. Pairs whose imaginary cross-spectrum is exactly zero for every observation (in-phase signals, and the diagonal) have no phase lag to estimate and are returned as 0 rather than the lower bound.

References

[1]

Vinck, M., Oostenveld, R., van Wingerden, M., Battaglia, F., and Pennartz, C.M.A. (2011). An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias. NeuroImage 55, 1548-1565.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> debiased_pli = connectivity.debiased_squared_phase_lag_index()
>>> debiased_pli.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> bool(debiased_pli[0, 20, 0, 1] > 0.05)  # lagged coupling at 10 Hz (bin 20)
True

debiased_squared_weighted_phase_lag_index() → ndarray[tuple[int, ...], dtype[floating]]

Return square of weighted phase lag index corrected for bias.

The square of the weighted phase lag index corrected for the positive bias induced by using the magnitude of the complex cross-spectrum.

Returns:

weighted_phase_lag_index – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Debiased squared weighted phase lag index values (symmetric in the signal pair, so they carry no lead/lag direction).

Return type:

array

Notes

Range: [-1, 1]. The debiased finite-sample estimate can be negative when the signed cross-products are dominated by inconsistent phase lags; negative values do not represent negative coupling. Pairs whose imaginary cross-spectrum is zero up to rounding (in-phase signals, the diagonal, and the purely real DC and Nyquist bins) have no phase lag to estimate and are returned as 0.

References

[1]

Vinck, M., Oostenveld, R., van Wingerden, M., Battaglia, F., and Pennartz, C.M.A. (2011). An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias. NeuroImage 55, 1548-1565.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> debiased_wpli = connectivity.debiased_squared_weighted_phase_lag_index()
>>> debiased_wpli.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> bool(debiased_wpli[0, 20, 0, 1] > 0.3)  # lagged coupling at 10 Hz (bin 20)
True

delay(frequencies_of_interest: ndarray[tuple[int, ...], dtype[floating]] | None = None, frequency_resolution: float | None = None, significance_threshold: float = 0.05, n_range: int = 3) → ndarray[tuple[int, ...], dtype[floating]]

Find a range of possible delays from the coherence phase.

The delay (and phase) at each frequency is indistinguishable from 2π phase jumps, but we can look at a range of possible delays and see which one is most likely.

Parameters:

• frequencies_of_interest (array-like, shape (2,), optional) – Frequency band (low, high) to evaluate, in the units of frequencies. Both edges are exclusive: only bins strictly inside the band are returned, and at least one must be. None uses every frequency.

• frequency_resolution (float, optional) – Frequency resolution for independent samples.

• significance_threshold (float, default=0.05) – P-value threshold for significance.

• n_range (int, default=3) – Number of phases to consider.

Returns:

possible_delays – Shape (…, n_frequencies, (n_range * 2) + 1, n_signals, n_signals), where n_frequencies counts only the frequencies inside frequencies_of_interest. Candidate k (index k + n_range) adds k cycles of phase. Array of possible time delays in the reciprocal units of frequencies (seconds for Hz, samples for cycles/sample); positive [..., i, j] means signal i leads signal j. The true delay is the candidate that is consistent (frequency-independent) across the band. Frequencies without significant coherence, and the 0 Hz (DC) bin, are undefined and returned as NaN.

Return type:

array

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> delays = connectivity.delay(frequencies_of_interest=[5, 50], n_range=3)
>>> delays.shape  # (n_time_windows, n_band_freqs, n_candidates, n_signals, n_signals)
(1, 89, 7, 2, 2)
>>> # The zero-wrap candidate (index n_range) recovers the 6 ms lead of signal 0.
>>> round(float(np.nanmedian(delays[0, :, 3, 0, 1])), 3)
0.006

direct_directed_transfer_function() → ndarray[tuple[int, ...], dtype[floating]]

Return the direct directed transfer function (dDTF).

The squared dDTF of Korzeniewska et al. (2003) multiplies the full-frequency DTF, which keeps direct and indirect (cascade) influence, by the squared partial coherence, which is zero between two signals whose relation is fully explained by the others. The product therefore keeps only direct influence:

chi^2_ij(f) = ffDTF^2_ij(f) * kappa^2_ij(f), with ffDTF^2_ij(f) = |H_ij(f)|^2 / sum_f' sum_k |H_ik(f')|^2 (summed over the non-negative frequencies) and kappa^2_ij(f) = |G_ij(f)|^2 / (G_ii(f) G_jj(f)), where G = A^H Sigma^-1 A is the inverse spectral matrix of the MVAR model (A = H^-1, Sigma the innovation covariance). These formulas are written in the transfer function’s native [target, source] indexing (chi^2_ij is j -> i).

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Returns:

direct_directed_transfer_function – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Output [..., i, j] is the direct influence of signal i on signal j (i -> j).

Return type:

array

Notes

Range: [0, 1]. Like directed_transfer_function(), this returns the squared quantity; SCoT and ConnectiviPy report its square root, |ffDTF| * |kappa|. The diagonal is the full-frequency DTF of each signal with itself (kappa_ii = 1).

References

[1]

Korzeniewska, A., Manczak, M., Kaminski, M., Blinowska, K.J., and Kasicki, S. (2003). Determination of information flow direction among brain structures by a modified directed transfer function (dDTF) method. Journal of Neuroscience Methods 125, 195-207.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> # Chain 0 -> 1 -> 2: each signal is the previous one delayed 3 samples, plus noise.
>>> source = rng.standard_normal((1006, 20))
>>> middle = source[:-3] + 0.5 * rng.standard_normal((1003, 20))
>>> target = middle[:-3] + 0.5 * rng.standard_normal((1000, 20))
>>> signals = np.stack([source[6:], middle[3:], target], axis=-1)  # (time, trials, 3)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> ddtf = connectivity.direct_directed_transfer_function()
>>> ddtf.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 3, 3)
>>> # [..., i, j] is i -> j. The direct 1 -> 2 ([..., 1, 2]) far exceeds both the
>>> # reverse 2 -> 1 ([..., 2, 1]) and the indirect 0 -> 2 ([..., 0, 2]) that is
>>> # relayed through signal 1 (10 Hz = bin 20).
>>> direct = ddtf[0, 20, 1, 2]
>>> bool(direct > 10 * ddtf[0, 20, 2, 1]), bool(direct > 10 * ddtf[0, 20, 0, 2])
(True, True)

directed_coherence() → ndarray[tuple[int, ...], dtype[floating]]

Return the squared directed coherence (noise-weighted DTF).

Like the directed transfer function, but the noise variance weights both the numerator and the inflow normalization. The returned value is the squared directed coherence nv_j |H_ij|^2 / sum_k nv_k |H_ik|^2, written in the transfer function’s native [target, source] indexing (H_ij is j -> i), where nv is the per-signal innovation (noise) variance and H is the transfer function. Each target’s values sum to 1 over sources (result.sum(axis=-2) is 1).

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Returns:

directed_coherence – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Squared directed coherence values. Output [..., i, j] is the influence of signal i on signal j (i -> j).

Return type:

array

Notes

Range: [0, 1]. Normalized directional connectivity measure.

Assumption: This measure follows Baccala et al. (1998), which assumes the MVAR innovations are uncorrelated (a diagonal noise covariance). The denominator then equals the signal’s power spectral density S_ii = sum_k nv_k |H_ik|^2. When the estimated innovation covariance has non-negligible off-diagonal terms (common for non-parametrically estimated MVARs), the true PSD S_ii = (H Cov H^H)_ii also contains cross-power between correlated sources that this diagonal formula omits, so the values are approximate. A UserWarning is emitted in that case.

References

[1]

Baccala, L., Sameshima, K., Ballester, G., Do Valle, A., and Timo-Iaria, C. (1998). Studying the interaction between brain structures via directed coherence and Granger causality. Applied Signal Processing 5, 40.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> directed_coherence = connectivity.directed_coherence()
>>> directed_coherence.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # [..., i, j] is i -> j, so 0 -> 1 is [..., 0, 1]. It dominates at 10 Hz (bin 20).
>>> bool(directed_coherence[0, 20, 0, 1] > 10 * directed_coherence[0, 20, 1, 0])
True
>>> bool(np.allclose(directed_coherence.sum(axis=-2), 1))  # each target's inflow
True

directed_phase_lag_index() → ndarray[tuple[int, ...], dtype[floating]]

Return the directed phase-lag index (dPLI).

Values above 0.5 indicate that one signal consistently phase-leads the other; values below 0.5 indicate that it phase-lags. A value of 0.5 represents no preferred phase-lag direction, including in-phase or anti-phase pairs, whose imaginary cross-spectrum is zero up to rounding.

Returns:

directed_phase_lag_index – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). [..., i, j] above 0.5 means signal i leads signal j; below 0.5 means signal i lags signal j.

Return type:

array

Notes

Range: [0, 1]. With the convention H(0) = 0.5, dPLI is (1 + signed_PLI) / 2. Consequently dPLI[i, j] = 1 - dPLI[j, i] and the diagonal is 0.5.

References

[1]

Stam, C.J., and van Straaten, E.C.W. (2012). Go with the flow: use of a directed phase lag index (dPLI) to characterize patterns of phase relations in a large-scale model of brain dynamics. NeuroImage 62, 1415-1428.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> dpli = connectivity.directed_phase_lag_index()
>>> dpli.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # Signal 0 leads signal 1, so [..., 0, 1] > 0.5 > [..., 1, 0] (10 Hz = bin 20).
>>> bool(dpli[0, 20, 0, 1] > 0.5 > dpli[0, 20, 1, 0])
True

directed_transfer_function() → ndarray[tuple[int, ...], dtype[floating]]

Return transfer function coupling strength normalized by inflow.

The transfer function coupling strength normalized by the total influence of other signals on that signal (inflow).

Characterizes the direct and indirect coupling to a node.

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Returns:

directed_transfer_function – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Directed transfer function values. Output [..., i, j] is the influence of signal i on signal j (i -> j).

Return type:

array

Notes

Range: [0, 1] (normalized). Represents proportion of inflow via transfer function; each target’s values sum to 1 over sources (result.sum(axis=-2) is 1).

References

[1]

Kaminski, M., and Blinowska, K.J. (1991). A new method of the description of the information flow in the brain structures. Biological Cybernetics 65, 203-210.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> dtf = connectivity.directed_transfer_function()
>>> dtf.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # [..., i, j] is i -> j, so 0 -> 1 is [..., 0, 1]. It dominates at 10 Hz (bin 20).
>>> bool(dtf[0, 20, 0, 1] > 10 * dtf[0, 20, 1, 0])
True
>>> bool(np.allclose(dtf.sum(axis=-2), 1))  # each target's inflow sums to 1
True

property expectation_type: str

Which dimensions the cross-spectral matrix is averaged over.

Reassigning clears all cached intermediates (see clear_cache()).

property fourier_coefficients: Any

Multitaper Fourier coefficients.

Shape (n_time_windows, n_trials, n_tapers, n_fft_samples, n_signals). The instance owns an immutable snapshot so cached calculations cannot become stale through in-place mutation. This accessor returns a detached copy, marked read-only when the backend supports it. Assign a new array through the setter to replace the coefficients and clear the caches.

property frequencies: ndarray[tuple[int, ...], dtype[floating]]

Return non-negative frequencies of the transform.

Returns:

Non-negative frequency values.

Return type:

NDArray[floating], shape (n_frequencies,)

classmethod from_multitaper(multitaper_instance: SpectralTransform, expectation_type: str = 'trials_tapers', dtype: ~typing.Any = <class 'numpy.complex128'>, minimum_phase_tolerance: float = 1e-08, minimum_phase_max_iterations: int = 500) → Connectivity

Construct from a spectral transform; the original name of from_transform.

Accepts any SpectralTransform.

Parameters:

• multitaper_instance (SpectralTransform) – The transform; transform in from_transform().

• expectation_type – As in from_transform().

• dtype – As in from_transform().

• minimum_phase_tolerance – As in from_transform().

• minimum_phase_max_iterations – As in from_transform().

Returns:

New Connectivity instance.

Return type:

Connectivity

classmethod from_transform(transform: SpectralTransform, expectation_type: str = 'trials_tapers', dtype: ~typing.Any = <class 'numpy.complex128'>, minimum_phase_tolerance: float = 1e-08, minimum_phase_max_iterations: int = 500) → Connectivity

Construct from any spectral transform.

This is the transform-neutral spelling of from_multitaper(); the older method remains fully supported.

Parameters:

• transform (SpectralTransform) – Such as Multitaper, MorletWavelet, or a user-defined class; see SpectralTransform for the required and optional members. fft() must return fresh, unshared storage on each call. Neither the transform nor its caller may subsequently mutate it through any alias: the result is used without copying and marked read-only where the backend supports this.

• expectation_type (str, default="trials_tapers") – How to average the cross-spectral matrix.

• dtype (np.dtype, default=complex128) – Data type for computations.

• minimum_phase_tolerance (float, default=1e-8) – Relative convergence tolerance for the Wilson minimum-phase factorization used by the directed measures.

• minimum_phase_max_iterations (int, default=500) – Maximum Wilson iterations. Increase for near-singular cross-spectral matrices (highly correlated channels) that fail to converge.

Returns:

New Connectivity instance.

Return type:

Connectivity

generalized_partial_directed_coherence() → ndarray[tuple[int, ...], dtype[floating]]

Return generalized partial directed coherence.

The transfer function coupling strength normalized by its strength of coupling to other signals (outflow).

The partial directed coherence tries to regress out the influence of other observed signals, leaving only the direct coupling between two signals.

The generalized partial directed coherence scales the relative strength of coupling by the noise variance.

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Returns:

generalized_partial_directed_coherence – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Generalized partial directed coherence values. Output [..., i, j] is the influence of signal i on signal j (i -> j).

Return type:

array

Notes

Range: [0, 1]. Normalized, scaled by noise variance; each source’s values sum to 1 over targets (result.sum(axis=-1) is 1).

References

[1]

Baccala, L.A., Sameshima, K., and Takahashi, D.Y. (2007). Generalized partial directed coherence. In Digital Signal Processing, 2007 15th International Conference on, (IEEE), pp. 163-166.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> gpdc = connectivity.generalized_partial_directed_coherence()
>>> gpdc.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # [..., i, j] is i -> j, so 0 -> 1 is [..., 0, 1]. It dominates at 10 Hz (bin 20).
>>> bool(gpdc[0, 20, 0, 1] > 10 * gpdc[0, 20, 1, 0])
True
>>> bool(np.allclose(gpdc.sum(axis=-1), 1))  # each source's outflow sums to 1
True

global_coherence(max_rank: int = 1, max_workspace_elements: int = 16000000) → tuple[ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[complexfloating]]]

Find linear combinations that capture the most coherent power.

The linear combinations of signals that capture the most coherent power at each frequency and time window.

This is a frequency domain analog of PCA over signals at a given frequency/time window.

Parameters:

• max_rank (int, default=1) – The number of components to keep (like the number of PC dimensions).

• max_workspace_elements (int, default=16_000_000) – Approximate working-set target, in array elements, for the batched decomposition: frequency bins are processed in chunks sized so the main intermediates stay near this many complex elements (the default ~16M ≈ 256 MB of complex128). It is a soft target, not a hard memory cap — it counts the dominant per-bin intermediates, not the outputs or LAPACK’s internal workspace, and it never goes below one bin per chunk, so actual peak memory is somewhat higher. Lower it to reduce peak memory on a constrained CPU or GPU (at the cost of more, smaller chunks); the default favors speed and does not change the result. Ignored on the per-bin fallback path used for a large decomposition dimension.

Returns:

• global_coherence (ndarray) – Shape (n_time_windows, n_fft_samples, n_components). The fraction of total coherent power captured by each component (eigenvalue of the cross-spectral matrix divided by the sum of all eigenvalues), ordered strongest component first.

• unnormalized_global_coherence (ndarray) – Shape (n_time_windows, n_fft_samples, n_signals, n_components). The global coherence vectors (left singular vectors).

Notes

Frequency axis: unlike every other public measure, which returns only the non-negative frequencies, this method returns all n_fft_samples bins of the (two-sided) transform in FFT order; index it with all_frequencies rather than frequencies.

Range: [0, 1]. Each value is the fraction of total coherent power in that component, so the measure is scale-invariant and the components sum to at most 1.

Algorithm: when the number of estimates (n_trials * n_tapers) is at least n_signals and n_signals is small (<= 64), the components are obtained from an eigendecomposition of the (n_signals, n_signals) cross-spectral matrix A @ Aᴴ rather than a singular value decomposition of A. This is substantially faster but squares the condition number, so for a nearly rank-deficient cross-spectral matrix (near-duplicate channels) the weakest returned components (large max_rank) may lose relative precision. The dominant component(s) — the usual use of this measure — are unaffected. A thin matrix (fewer estimates than signals) uses the economy SVD directly.

References

[1]

Cimenser, A., Purdon, P.L., Pierce, E.T., Walsh, J.L., Salazar-Gomez, A.F., Harrell, P.G., Tavares-Stoeckel, C., Habeeb, K., and Brown, E.N. (2011). Tracking brain states under general anesthesia by using global coherence analysis. Proceedings of the National Academy of Sciences 108, 8832-8837.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> global_coherence, vectors = connectivity.global_coherence(max_rank=1)
>>> # Unlike the pairwise measures, the frequency axis spans all FFT bins.
>>> global_coherence.shape  # (n_time_windows, n_fft_samples, n_components)
(1, 1000, 1)
>>> vectors.shape  # (n_time_windows, n_fft_samples, n_signals, n_components)
(1, 1000, 2, 1)
>>> # One component captures most of the power of the coupled pair (10 Hz = bin 20).
>>> bool(global_coherence[0, 20, 0] > 0.8)
True

group_delay(frequencies_of_interest: ndarray[tuple[int, ...], dtype[floating]] | None = None, frequency_resolution: float | None = None, significance_threshold: float = 0.05) → tuple[ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[floating]]]

Return the average time-delay of a broadband signal.

Parameters:

• frequencies_of_interest (array-like, shape (2,), optional) – Frequency band (low, high) to fit over, in the units of frequencies. Both edges are exclusive: only bins strictly inside the band are used, and at least one must be. None uses every frequency.

• frequency_resolution (float, optional) – Frequency resolution for independent samples.

• significance_threshold (float, default=0.05) – P-value threshold for significance.

Returns:

• delay (array, shape (…, n_signals, n_signals)) – Time delays between signal pairs, in the reciprocal units of frequencies: seconds for Hz, samples for cycles/sample. Positive [..., i, j] means signal i leads signal j. The diagonal is NaN.

• slope (array, shape (…, n_signals, n_signals)) – Slope of the coherence phase vs frequency, in radians per unit of frequencies (delay = slope / (2 * pi)); same sign convention as delay.

• r_value (array, shape (…, n_signals, n_signals)) – Correlation coefficient of the linear phase-frequency fit, with the sign of slope (so [..., j, i] is -[..., i, j]).

Notes

Range: (-∞, ∞). Time delays can be positive or negative.

References

[1]

Gotman, J. (1983). Measurement of small time differences between EEG channels: method and application to epileptic seizure propagation. Electroencephalography and Clinical Neurophysiology 56, 501-514.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> delay, slope, r_value = connectivity.group_delay(frequencies_of_interest=[5, 50])
>>> delay.shape  # (n_time_windows, n_signals, n_signals)
(1, 2, 2)
>>> # Signal 0 leads signal 1 by 6 ms, so [..., 0, 1] is +0.006 s.
>>> round(float(delay[0, 0, 1]), 3), round(float(delay[0, 1, 0]), 3)
(0.006, -0.006)
>>> bool(r_value[0, 0, 1] > 0.9)  # phase is linear in frequency
True

imaginary_coherence() → ndarray[tuple[int, ...], dtype[floating]]

Return the normalized imaginary component of the cross-spectrum.

Projects the cross-spectrum onto the imaginary axis to mitigate the effect of volume-conducted dependencies. Assumes volume-conducted sources arrive at sensors at the same time, resulting in a cross-spectrum with phase angle of 0 (perfectly in-phase) or π (anti-phase) if the sensors are on opposite sides of a dipole source. With the imaginary coherence, in-phase and anti-phase associations are set to zero.

Returns:

imaginary_coherence_magnitude – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Imaginary coherence magnitudes (symmetric in the signal pair; see imaginary_coherency() for the signed, lead/lag-aware version).

Return type:

array

Notes

Range: [0, 1]. Magnitude version of imaginary part of coherency. Raw imaginary component ranges in [-1, 1].

References

[1]

Nolte, G., Bai, O., Wheaton, L., Mari, Z., Vorbach, S., and Hallett, M. (2004). Identifying true brain interaction from EEG data using the imaginary part of coherency. Clinical Neurophysiology 115, 2292-2307.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> imaginary_coherence = connectivity.imaginary_coherence()
>>> imaginary_coherence.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # The 6 ms lag is not zero-phase, so the imaginary part survives (10 Hz = bin 20).
>>> bool(imaginary_coherence[0, 20, 0, 1] > 0.2)
True

imaginary_coherency() → ndarray[tuple[int, ...], dtype[floating]]

Return the signed imaginary component of coherency.

This is the signed counterpart of imaginary_coherence(), which returns its magnitude. The sign is antisymmetric across a signal pair and preserves the pair’s phase-lead/phase-lag orientation.

Returns:

imaginary_coherency – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Positive [..., i, j] means signal i leads signal j.

Return type:

array

Notes

Range: [-1, 1]. The diagonal and pairs involving zero-power signals are undefined and returned as NaN, matching coherency().

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> imaginary_coherency = connectivity.imaginary_coherency()
>>> imaginary_coherency.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # Signal 0 leads, so [..., 0, 1] is positive and [..., 1, 0] negative (10 Hz).
>>> bool(imaginary_coherency[0, 20, 0, 1] > 0 > imaginary_coherency[0, 20, 1, 0])
True

property is_one_sided: bool

Whether the input contains only non-negative frequencies.

jackknife(method: str, *, confidence_level: float = 0.95, transformation: Literal['auto', 'identity', 'log', 'fisher', 'fisher_squared', 'circular'] = 'auto', **method_kwargs: Any) → JackknifeResult

Estimate uncertainty by leaving out one trial/taper observation.

The configured expectation must average trials, tapers, or their combination. For trials_tapers each trial-taper eigencoefficient is treated as one observation. The method is recomputed for every leave-one-out sample, so it supports nonlinear measures without an analytic variance formula (at a computational cost proportional to n_observations).

Parameters:

• method (str) – Name of a public, real-valued connectivity measure to recompute, e.g. "coherence_magnitude". Complex-valued and tuple-valued measures are not supported.

• confidence_level (float, default=0.95) – Two-sided coverage of the interval, in (0, 1). The critical value is Student t with n_observations - 1 degrees of freedom.

• transformation ({"auto", "identity", "log", "fisher",) –

  “fisher_squared”, “circular”} Scale on which the interval is formed. "auto" resolves to:

  • "log" for power;

  • "fisher_squared" (atanh(sqrt(.))) for the magnitude-squared measures in [0, 1]: coherence_magnitude and partial_coherence;

  • "fisher" (atanh(.)) for the magnitudes in [0, 1]: phase_locking_value and imaginary_coherence. For these two measures the Fisher lower bound and bias-corrected estimate are clamped at 0, since a magnitude cannot be negative;

  • "circular" for coherence_phase;

  • "identity" for every other measure.

• **method_kwargs – Keyword arguments forwarded to method on every replicate.

Returns:

Estimate, bias-corrected estimate, standard error, and confidence bounds, each with the measure’s own shape (for a pairwise measure (n_time, n_nonnegative_frequencies, n_signals, n_signals)). For directed measures [..., i, j] is the influence i -> j, the same order as the xarray wrapper’s sel(source=i, target=j).

Return type:

JackknifeResult

Raises:

ValueError – If the expectation averages tapers ("tapers" or "trials_tapers") but observations_are_independent is False: correlated observations (a MorletWavelet smoothing neighborhood, overlapping Welch segments) are not valid leave-one-out units. Recompute the transform with independent observations instead: Welch with segment_overlap <= 0.5, a Multitaper transform, or MorletWavelet without smoothing_time on at least 3 trials. (expectation_type="trials" is accepted but keeps the correlated axis as an output dimension instead of averaging it, a different measure.)

Notes

Circular confidence bounds are wrapped to (-pi, pi], so when the interval crosses +/-pi the lower bound exceeds the upper bound; the interval is then the arc from lower up through pi and on to upper.

If the input Fourier coefficients were produced with Multitaper(taper_weighting="adaptive"), the leave-one-out replicates reuse the full-sample Thomson weights (the adaptive weights are not recomputed for each reduced taper set), so the interval is an approximation in that case.

The interval describes the size of a measure; it is not a test that the measure differs from 0. Monte Carlo coverage of the default ("auto") 95% intervals, with independent complex-Gaussian observations and 5 to 100 observations per dataset:

• coherence_magnitude (fisher_squared): 94-95% when the true |coherency| is 0.3-0.8, but at zero true coherence the interval excludes 0 in 11-15% of datasets, and more observations do not help (the estimated magnitude’s spread shrinks with n at the same rate as its mean). “The interval excludes 0” is therefore not evidence of nonzero coherence. Test that with the exact zero-coherence null instead: spectral_connectivity.statistics.coherence_significance_pvalue() applied to coherency() and n_observations (valid for independent, equally weighted observations), corrected across frequencies and pairs with spectral_connectivity.statistics.adjust_for_multiple_comparisons(). For measures without an analytic null, use a permutation or surrogate test.

• phase_locking_value (fisher): at zero true PLV the interval excludes 0 in 5% of datasets at 5 observations but 14% at 100; at a high true PLV it under-covers (85-95% at a true PLV of 0.82, 76-90% at 0.93).

• coherence_phase (circular): under-covers when coherence is weak (77-88% at a true |coherency| of 0.1, 84-93% at 0.2, 87-94% at 0.3, against 92-95% at 0.6).

maximized_imaginary_coherency(group_labels: ndarray[tuple[int, ...], dtype[integer]], rank: int | None = None, regularization: float = 1e-12) → tuple[ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[integer]]]

Return maximized imaginary coherency (MIC) between signal groups.

Each group’s real within-group cross-spectrum is whitened before the largest singular value of the between-group imaginary cross-spectrum is taken. This makes the result invariant to invertible, static real-valued mixing within either group.

Parameters:

• group_labels (array-like, shape (n_signals,)) – Label assigning each signal to a group.

• rank (int, optional) – Retain at most this many within-group whitening components. None retains every numerically non-zero component independently per bin.

• regularization (float, default=1e-12) – Relative diagonal loading used by the whitening decomposition.

Returns:

• mic (array) – Shape (..., n_nonnegative_frequencies, n_groups, n_groups).

• labels (array, shape (n_groups,)) – Sorted unique group labels.

Notes

Range: [0, 1]. The diagonal is returned as NaN.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Group "a": two noisy copies of signal 0; "b": two of its 6 ms-delayed copy.
>>> pair = np.stack([leader[3:], leader[:-3]], axis=-1)
>>> signals = np.repeat(pair, 2, axis=-1)  # (time, trials, 4 signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> mic, labels = connectivity.maximized_imaginary_coherency(["a", "a", "b", "b"])
>>> mic.shape  # (n_time_windows, n_frequencies, n_groups, n_groups)
(1, 501, 2, 2)
>>> labels
array(['a', 'b'], dtype='<U1')
>>> bool(mic[0, 20, 0, 1] > 0.2)  # lagged a-b coupling at 10 Hz (bin 20)
True

maximized_imaginary_coherency_components(group_labels: ndarray[tuple[int, ...], dtype[Any]], *, rank: int | None = None, n_components: int = 1, regularization: float = 1e-12) → MultivariateConnectivityResult

Return component-resolved MIC scores, filters, and patterns.

The singular vectors of the whitened imaginary between-group CSD are returned in descending singular-value order. Filters map channel data to the components; patterns map the components back to channel space. This is the component-resolved counterpart of the scalar maximized_imaginary_coherency().

Parameters:

• group_labels (array-like, shape (n_signals,)) – Label assigning each signal to a group; every unordered pair of groups is one connection.

• rank (int, optional) – Retain at most this many within-group whitening directions per group. None keeps every numerically non-zero direction.

• n_components (int, default=1) – Number of singular components to return. A connection whose smaller group has fewer channels returns NaN for the unavailable components.

• regularization (float, default=1e-12) – Relative diagonal loading used by the whitening decomposition.

Returns:

Real scores of shape (..., frequency, connection, component) plus filters and patterns; see the class docstring.

Return type:

MultivariateConnectivityResult

Notes

Range: scores lie in [0, 1]. The whitening and singular-value decomposition are vectorized over the time/frequency axes on the active xp backend, so this runs on the GPU when GPU support is enabled.

References

[1]

Ewald A, et al. (2012) Estimating true brain connectivity from EEG/ MEG data invariant to linear and static transformations in sensor space. NeuroImage 60(1):476-488.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Group "a": two noisy copies of signal 0; "b": two of its 6 ms-delayed copy.
>>> pair = np.stack([leader[3:], leader[:-3]], axis=-1)
>>> signals = np.repeat(pair, 2, axis=-1)  # (time, trials, 4 signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> labels = ["a", "a", "b", "b"]
>>> result = connectivity.maximized_imaginary_coherency_components(labels)
>>> result.scores.shape  # (n_time_windows, n_frequencies, n_connections, n_components)
(1, 501, 1, 1)
>>> result.connections  # one row per group pair
array([['a', 'b']], dtype='<U1')
>>> result.patterns.shape  # (..., n_connections, n_components, side, n_signals)
(1, 501, 1, 1, 2, 4)
>>> bool(result.scores[0, 20, 0, 0] > 0.2)  # lagged a-b coupling at 10 Hz (bin 20)
True

minimum_phase_reconstruction_error() → ndarray[tuple[int, ...], dtype[floating]]

Return the relative reconstruction error of the Wilson factorization.

This diagnostic checks how faithfully the cached minimum-phase factor reconstructs the expected cross-spectral matrix. One value is returned per retained batch (normally time-window) dimension. Values near machine precision indicate a faithful factorization; large values suggest that the spectrum is too coarsely resolved for directed-connectivity measures. Non-converged factorizations return NaN.

Returns:

Maximum relative reconstruction error for each sub-spectrum.

Return type:

array, shape (…,)

Notes

A full two-sided spectrum is required, as for the directed measures that use the Wilson factorization.

multivariate_interaction_measure(group_labels: ndarray[tuple[int, ...], dtype[integer]], rank: int | None = None, regularization: float = 1e-12) → tuple[ndarray[tuple[int, ...], dtype[floating]], ndarray[tuple[int, ...], dtype[integer]]]

Return the multivariate interaction measure (MIM) between groups.

MIM sums the squared singular values of the whitened imaginary cross-spectrum, incorporating every phase-lagged interaction component rather than only the strongest component returned by MIC.

Parameters:

• group_labels (array-like, shape (n_signals,)) – Label assigning each signal to a group.

• rank (int, optional) – Retain at most this many within-group whitening components. None retains every numerically non-zero component independently per bin.

• regularization (float, default=1e-12) – Relative diagonal loading used by the whitening decomposition.

Returns:

• mim (array) – Shape (..., n_nonnegative_frequencies, n_groups, n_groups).

• labels (array, shape (n_groups,)) – Sorted unique group labels.

Notes

Range: [0, min(rank_group_1, rank_group_2)]; unlike MIC, MIM can exceed one. The diagonal is returned as NaN.

References

[1]

Ewald, A., Marzetti, L., Zappasodi, F., Meinecke, F.C., and Nolte, G. (2012). Estimating true brain connectivity from EEG/MEG data invariant to linear and static transformations in sensor space. NeuroImage 60, 476-488.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Group "a": two noisy copies of signal 0; "b": two of its 6 ms-delayed copy.
>>> pair = np.stack([leader[3:], leader[:-3]], axis=-1)
>>> signals = np.repeat(pair, 2, axis=-1)  # (time, trials, 4 signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> mim, labels = connectivity.multivariate_interaction_measure(["a", "a", "b", "b"])
>>> mim.shape  # (n_time_windows, n_frequencies, n_groups, n_groups)
(1, 501, 2, 2)
>>> labels
array(['a', 'b'], dtype='<U1')
>>> bool(mim[0, 20, 0, 1] > 0.05)  # lagged a-b interaction at 10 Hz (bin 20)
True

property n_observations: int

Return the raw number of observations averaged by the expectation.

Returns:

Product of the lengths of the averaged observation axes (for the default "trials_tapers" expectation, n_trials * n_tapers). This is a raw count, not an effective number of independent observations: it ignores observation_weights and is not reduced when the transform’s observations are correlated (see observations_are_independent).

Return type:

int

property n_signals: int

Number of signals represented by the Fourier coefficients.

property observation_weights: Any | None

Non-negative weights used when averaging spectral observations.

Weights have shape (time, trial, taper, frequency, 1) and are shared by every signal. A detached, read-only copy is returned so cached expectations cannot be invalidated by mutation.

property observations_are_independent: bool

Whether the trial/taper observations are statistically independent.

False when the transform collected correlated estimates on the observation axis (a smoothing neighborhood of a MorletWavelet, or Welch segments overlapping by more than half). n_observations then overstates the effective sample size, so the measures that rely on it warn and jackknife refuses to treat tapers as leave-one-out units.

pairwise_phase_consistency() → ndarray[tuple[int, ...], dtype[floating]]

Return square of phase locking value corrected for bias.

The square of the phase locking value corrected for the positive bias induced by using the magnitude of the complex cross-spectrum.

Returns:

phase_locking_value – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Pairwise phase consistency values.

Return type:

array

Notes

Range: [-1 / (n_observations - 1), 1]. The unbiased finite-sample estimate can be negative when phase consistency is below its null bias; negative values do not represent negative coupling.

References

[1]

Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P., and Pennartz, C.M.A. (2010). The pairwise phase consistency: A bias-free measure of rhythmic neuronal synchronization. NeuroImage 51, 112-122.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> ppc = connectivity.pairwise_phase_consistency()
>>> ppc.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> bool(ppc[0, 20, 0, 1] > 0.3)  # consistent phase difference at 10 Hz (bin 20)
True

pairwise_spectral_granger_prediction() → ndarray[tuple[int, ...], dtype[floating]]

Return amount of power at a node explained by other nodes.

The amount of power at a node in a frequency explained by (is predictive of) the power at other nodes.

Also known as spectral granger causality.

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Returns:

pairwise_granger – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Spectral Granger prediction values. Output [..., i, j] is the influence of signal i on signal j (i -> j).

Return type:

array

Notes

Non-negativity: spectral Granger is >= 0 by definition. Negative estimates within roundoff of zero – above -100 * eps of the result dtype on the log-ratio scale, about -2e-14 for float64 – are clipped to 0; materially negative bins below that threshold (a degenerate factorization) are returned as NaN; use minimum_phase_reconstruction_error() to diagnose them. Other packages (FieldTrip, MVGC, mne-connectivity) return such values as-is.

Range: [0, ∞). Non-negative values with no finite upper bound.

References

[1]

Geweke, J. (1982). Measurement of Linear Dependence and Feedback Between Multiple Time Series. Journal of the American Statistical Association 77, 304.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> granger = connectivity.pairwise_spectral_granger_prediction()
>>> granger.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # [..., i, j] is i -> j, so 0 -> 1 is [..., 0, 1]. It dominates at 10 Hz (bin 20).
>>> bool(granger[0, 20, 0, 1] > 10 * granger[0, 20, 1, 0])
True

partial_coherence(regularization: float = 1e-12) → ndarray[tuple[int, ...], dtype[floating]]

Return magnitude-squared coherence conditional on all other signals.

Partial coherence is computed by normalizing the off-diagonal elements of the inverse cross-spectral density (the spectral precision matrix). It measures the remaining linear association between each pair after conditioning on every other observed signal.

Parameters:

regularization (float, default=1e-12) – Non-negative relative diagonal loading applied independently to each time-frequency cross-spectral matrix before inversion. The absolute loading is regularization * rms(abs(S)). Increase this value for statistically rank-deficient or ill-conditioned spectra.

Returns:

partial_coherence – Shape (..., n_nonnegative_frequencies, n_signals, n_signals).

Return type:

array

Notes

Range: [0, 1]. The diagonal is undefined and returned as NaN. This undirected measure is distinct from partial directed coherence. Regularization stabilizes inversion but also changes the estimand, so analyses should report a non-default value.

The averaged trials x tapers are the observations. With fewer observations than signals the cross-spectral matrix is rank-deficient and its (regularized) inverse is dominated by the null space; with exactly one null direction every partial coherence is forced to 1 for any data. This case warns.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1000, 20, 1))
>>> # Signals 1 and 2 are noisy copies of signal 0 and share nothing else.
>>> copies = leader + 0.5 * rng.standard_normal((1000, 20, 2))
>>> signals = np.concatenate([leader, copies], axis=-1)  # (time, trials, signals)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> partial = connectivity.partial_coherence()
>>> partial.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 3, 3)
>>> # Signals 1 and 2 are coherent, but not once signal 0 is accounted for (10 Hz).
>>> bool(connectivity.coherence_magnitude()[0, 20, 1, 2] > 0.5)
True
>>> bool(partial[0, 20, 1, 2] < 0.05)
True

partial_directed_coherence() → ndarray[tuple[int, ...], dtype[floating]]

Return transfer function coupling strength normalized by outflow.

The transfer function coupling strength normalized by its strength of coupling to other signals (outflow).

The partial directed coherence tries to regress out the influence of other observed signals, leaving only the direct coupling between two signals.

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Returns:

partial_directed_coherence – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Partial directed coherence values. Output [..., i, j] is the influence of signal i on signal j (i -> j).

Return type:

array

Notes

Range: [0, 1]. Normalized direct coupling measure; each source’s values sum to 1 over targets (result.sum(axis=-1) is 1).

References

[1]

Baccala, L.A., and Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. Biological Cybernetics 84, 463-474.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> pdc = connectivity.partial_directed_coherence()
>>> pdc.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # [..., i, j] is i -> j, so 0 -> 1 is [..., 0, 1]. It dominates at 10 Hz (bin 20).
>>> bool(pdc[0, 20, 0, 1] > 10 * pdc[0, 20, 1, 0])
True
>>> bool(np.allclose(pdc.sum(axis=-1), 1))  # each source's outflow sums to 1
True

phase_lag_index() → ndarray[tuple[int, ...], dtype[floating]]

Return non-parametric synchrony measure mitigating power differences.

A non-parametric synchrony measure designed to mitigate power differences between realizations (tapers, trials) and volume-conduction.

The phase lag index is the average sign of the imaginary component of the cross-spectrum. The imaginary component sets in-phase or anti-phase signals to zero and the sign scales it to have the same magnitude regardless of phase.

Note that this is the signed version of the phase lag index. In order to obtain the unsigned version, as in [1], take the absolute value of this quantity.

Returns:

phase_lag_index – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Signed phase lag index values for all signal pairs. Positive [..., i, j] means signal i leads signal j.

Return type:

array

Notes

Range: [-1, 1] (signed version). For unsigned version (as in [1]), take absolute value to get range [0, 1]. In-phase or anti-phase pairs, whose imaginary cross-spectrum is zero up to rounding, are 0.

References

[1]

Stam, C.J., Nolte, G., and Daffertshofer, A. (2007). Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Human Brain Mapping 28, 1178-1193.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> pli = connectivity.phase_lag_index()
>>> pli.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # Signal 0 leads, so [..., 0, 1] is positive and [..., 1, 0] negative (10 Hz).
>>> bool(pli[0, 20, 0, 1] > 0 > pli[0, 20, 1, 0])
True

phase_locking_value() → ndarray[tuple[int, ...], dtype[floating]]

Return the cross-spectrum with power scaled to magnitude 1.

The phase locking value attempts to mitigate power differences between realizations (tapers or trials) by treating all values of the cross-spectrum as the same power. This has the effect of downweighting high power realizations and upweighting low power realizations.

Returns:

phase_locking_value – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Phase locking values between all signal pairs.

Return type:

array

Notes

Range: [0, 1]. 0 indicates random phases; 1 indicates constant phase difference.

References

[1]

Lachaux, J.-P., Rodriguez, E., Martinerie, J., Varela, F.J., and others (1999). Measuring phase synchrony in brain signals. Human Brain Mapping 8, 194-208.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> plv = connectivity.phase_locking_value()
>>> plv.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> bool(plv[0, 20, 0, 1] > 0.5)  # consistent phase difference at 10 Hz (bin 20)
True

phase_slope_index(frequencies_of_interest: ndarray[tuple[int, ...], dtype[floating]] | None = None, frequency_resolution: float | None = None) → ndarray[tuple[int, ...], dtype[floating]]

Return weighted average of slopes projected onto imaginary axis.

The phase slope index sums the product of the coherency at adjacent frequencies, conj(C(f)) * C(f + df), over the band and projects the result onto the imaginary axis to avoid volume-conduction effects (Nolte et al. 2008). The magnitude of the coherency at each frequency therefore weights the contribution of that frequency step.

Parameters:

• frequencies_of_interest (array-like, shape (2,), optional) – Frequency band (low, high) to sum over, in the units of frequencies. Both edges are exclusive: only bins strictly inside the band contribute, and at least one must. None uses every frequency.

• frequency_resolution (float, optional) – Frequency resolution for independent samples.

Returns:

phase_slope_index – Phase slope index values. Positive [..., i, j] means signal i leads signal j; the result is antisymmetric.

Return type:

array, shape (…, n_signals, n_signals)

Notes

Range: (-∞, ∞). Signed directional measure with no bounds.

References

[1]

Nolte, G., Ziehe, A., Nikulin, V.V., Schlogl, A., Kramer, N., Brismar, T., and Muller, K.-R. (2008). Robustly Estimating the Flow Direction of Information in Complex Physical Systems. Physical Review Letters 100.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> psi = connectivity.phase_slope_index(frequencies_of_interest=[5, 50])
>>> psi.shape  # (n_time_windows, n_signals, n_signals)
(1, 2, 2)
>>> # Signal 0 leads signal 1, so [..., 0, 1] is positive and [..., 1, 0] negative.
>>> bool(psi[0, 0, 1] > 0 > psi[0, 1, 0])
True

power() → ndarray[tuple[int, ...], dtype[floating]]

Return the one-sided power spectral density of the signal.

Only the non-negative frequencies are returned, with the interior positive-frequency bins doubled so that integrating the returned spectrum over frequency recovers the full signal power (the negative frequencies of a real signal carry equal power). The DC bin, and the Nyquist bin for an even FFT length, are not doubled.

For one-sided input (is_one_sided=True) the coefficients are returned as provided: a one-sided transform such as MorletWavelet already scales its coefficients to the one-sided density, and externally supplied one-sided coefficients keep whatever scale they were given.

Returns:

One-sided power spectral density for non-negative frequencies, shape (..., n_nonnegative_frequencies, n_signals).

Return type:

NDArray[floating]

Notes

Range: [0, ∞). Power spectral density is always non-negative with no finite upper bound.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> power = connectivity.power()
>>> power.shape  # (n_time_windows, n_frequencies, n_signals)
(1, 501, 2)
>>> connectivity.frequencies[[0, -1]]  # 0 Hz up to Nyquist in 0.5 Hz steps
array([  0., 250.])

subset_pairwise_spectral_granger_prediction(pairs: Sequence[Sequence[int]] | ndarray[tuple[int, ...], dtype[integer]]) → ndarray[tuple[int, ...], dtype[floating]]

Return predictive power for a subset of signal pairs.

Changed in version 3.0: The output is source first ([..., i, j] is i -> j); 2.x returned the transpose (j -> i).

Parameters:

pairs (array_like, shape (n_pairs, 2)) – Pairs of signal indices. Each pair is estimated in both directions.

Returns:

pairwise_granger – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Spectral Granger prediction for the specified pairs; entries for pairs not requested (and the diagonal) are NaN. Output [..., i, j] is the influence of signal i on signal j (i -> j).

Return type:

array

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz); signal 2 is noise.
>>> signals = np.stack([leader[3:], leader[:-3], np.zeros((1000, 20))], axis=-1)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> granger = connectivity.subset_pairwise_spectral_granger_prediction(pairs=[(0, 1)])
>>> granger.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 3, 3)
>>> # [..., i, j] is i -> j, so 0 -> 1 is [..., 0, 1]. It dominates at 10 Hz (bin 20).
>>> bool(granger[0, 20, 0, 1] > 10 * granger[0, 20, 1, 0])
True
>>> bool(np.isnan(granger[0, 20, 0, 2]))  # pair (0, 2) was not requested
True

property time_bins_are_independent: bool

Whether the time bins may be counted as independent observations.

False when the transform’s successive time bins are correlated (Multitaper or ShortTimeFourierTransform windows overlapping by more than half, or MorletWavelet samples closer than four wavelet standard deviations). It matters only when the expectation averages over time: n_observations then counts the correlated bins and overstates the effective sample size, so the measures that rely on it warn.

time_reversed_spectral_granger_prediction() → ndarray[tuple[int, ...], dtype[floating]]

Return pairwise spectral Granger prediction after time reversal.

For a real stationary process, time reversal transposes the cross-spectral matrix at every frequency. Contrasting this result with pairwise_spectral_granger_prediction() helps identify apparent directionality caused by instantaneous mixing or data asymmetries.

Returns:

time_reversed_granger – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Output [..., i, j] is the influence of signal i on signal j (i -> j) in the time-reversed data.

Return type:

array

Notes

Non-negativity: spectral Granger is >= 0 by definition. Negative estimates within roundoff of zero – above -100 * eps of the result dtype on the log-ratio scale, about -2e-14 for float64 – are clipped to 0; materially negative bins below that threshold (a degenerate factorization) are returned as NaN; use minimum_phase_reconstruction_error() to diagnose them. Other packages (FieldTrip, MVGC, mne-connectivity) return such values as-is.

References

[1]

Winkler, I., Panknin, D., Bartz, D., Müller, K.-R., and Haufe, S. (2016). Validity of time reversal for testing Granger causality. IEEE Transactions on Signal Processing 64, 2746-2760.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> reversed_granger = connectivity.time_reversed_spectral_granger_prediction()
>>> reversed_granger.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # A genuine 0 -> 1 lead flips under time reversal: 1 -> 0 ([..., 1, 0]) dominates.
>>> bool(reversed_granger[0, 20, 1, 0] > 10 * reversed_granger[0, 20, 0, 1])
True

weighted_phase_lag_index() → ndarray[tuple[int, ...], dtype[floating]]

Return weighted average of phase lag index using imaginary coherency magnitudes.

Weighted average of the phase lag index using the imaginary coherency magnitudes as weights.

Note that this is the signed version of the weighted phase lag index (mirroring phase_lag_index()). In order to obtain the unsigned version, as in [1], take the absolute value of this quantity.

Returns:

weighted_phase_lag_index – Shape (..., n_nonnegative_frequencies, n_signals, n_signals). Signed weighted phase lag index values. Positive [..., i, j] means signal i leads signal j.

Return type:

array

Notes

Range: [-1, 1] (signed version). For the unsigned version (as in [1]), take the absolute value to get range [0, 1]. The result is antisymmetric in the signal pair (wpli[..., i, j] = -wpli[..., j, i]).

References

[1]

Vinck, M., Oostenveld, R., van Wingerden, M., Battaglia, F., and Pennartz, C.M.A. (2011). An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias. NeuroImage 55, 1548-1565.

Examples

>>> import numpy as np
>>> from spectral_connectivity import Connectivity, Multitaper
>>> rng = np.random.default_rng(0)
>>> leader = rng.standard_normal((1003, 20))
>>> # Signal 1 is signal 0 delayed by 3 samples (6 ms at 500 Hz), plus noise.
>>> signals = np.stack([leader[3:], leader[:-3]], axis=-1)  # (time, trials, signals)
>>> signals += 0.5 * rng.standard_normal(signals.shape)
>>> multitaper = Multitaper(signals, sampling_frequency=500)
>>> connectivity = Connectivity.from_transform(multitaper)
>>> wpli = connectivity.weighted_phase_lag_index()
>>> wpli.shape  # (n_time_windows, n_frequencies, n_signals, n_signals)
(1, 501, 2, 2)
>>> # Signal 0 leads, so [..., 0, 1] is positive and [..., 1, 0] negative (10 Hz).
>>> bool(wpli[0, 20, 0, 1] > 0 > wpli[0, 20, 1, 0])
True