spectral_connectivity.wrapper.multitaper_connectivity#

Compute connectivity measures with multitaper spectral estimation.

This is the main high-level function for connectivity analysis. It performs multitaper spectral analysis on the input time series and computes the requested connectivity measures, returning results as labeled xarray objects.

Changed in version 3.0: For pairwise_spectral_granger_prediction and subset_pairwise_spectral_granger_prediction, sel(source=a, target=b) is a -> b; 2.x returned b -> a. The 2.x wrapper rejected the transfer-function measures, and phase_slope_index, group_delay, and delay are unchanged.

Parameters:

time_series (NDArray[floating] or xarray.DataArray,) – shape (n_times, n_trials, n_channels) or (n_times, n_channels) Time series data. For multiple trials, trials are averaged in spectral domain. For a DataArray, common time/trial/signal dimension names are inferred and transposed automatically; use time_dim, trial_dim, and signal_dim for domain-specific names. Ambiguous names raise rather than falling back to dimension position, though a single unrecognized dimension left for the one remaining role is assigned by elimination with a warning (a spectral name such as frequency or band is rejected instead, since it marks an already-transformed input). A dask-backed DataArray is rejected (materialize it first with DataArray.compute()). A numeric time index is interpreted as elapsed seconds (a sample index as sample numbers) and used to label output window centers. When sampling_frequency is given it is validated against the index; when it is omitted, an elapsed-seconds time coordinate infers it (a sample index cannot, having no time scale). Datetime, timedelta, and object-valued time coordinates are not yet supported and must first be converted to numeric elapsed seconds. When signal_names is omitted, labels from a 1-D index coordinate on the signal dimension are carried to the output’s source and target coordinates without changing their type; if such labels are present but unusable a warning is issued and default string labels are used.
sampling_frequency (float, optional) – Sampling rate in Hz of the time series data. Required for array input; for a DataArray it may be omitted and inferred from a sufficiently precise numeric elapsed-seconds time coordinate. Pass it explicitly when a low-precision or large-offset coordinate cannot resolve the rate reliably.
time_window_duration (float, optional) – Duration of sliding window in seconds for time-resolved analysis. If None, analyzes entire time series (no time resolution).
method (str or list of str, optional) – Connectivity method(s) to compute. If None, computes the default set of real-valued measures that fit the xarray/NetCDF interface (see DEFAULT_METHODS) — not every measure. coherency is left out of the default because complex arrays are not portably serializable across all supported xarray versions and NetCDF engines, but it can be requested by name. Every directed measure is opt-in by name, including pairwise_spectral_granger_prediction and the directed-transfer-function family (directed_transfer_function, directed_coherence, partial_directed_coherence, generalized_partial_directed_coherence, direct_directed_transfer_function); the spectral Granger and transfer-function measures factorize the spectrum and dominate the cost (see the Notes on directed orientation). Measures with nonstandard layouts, including global_coherence, phase_slope_index, group_delay, delay, canonical_coherence, and blockwise spectral Granger, are available by name and return labeled DataArrays or Datasets with their component, group, candidate-delay, or frequency-reduced dimensions. Examples: “coherence_magnitude”, “imaginary_coherence”, “phase_locking_value”.
signal_names (sequence of scalar, optional) – Scalar, non-missing, unique xarray-compatible coordinate labels for signal channels. Integer labels must fit the signed 32-bit range for portable NetCDF3 serialization. Nested or structured labels are not supported. If None, uses the DataArray signal index when available, otherwise stringified indices.
squeeze (bool, default=False) – Only honored when a single method (a string) is requested, so the result is a DataArray. If there are exactly 2 channels, reduce a pairwise measure to the single ordered pair (first source, last target), returning a (time, frequency) array whose selected source and target are retained as scalar coordinates – so the pair (and, for directed measures, the direction) is still recorded. With more than 2 channels a warning is issued and the full matrix is returned; for power (no target axis) squeeze is a no-op. For multi-measure requests (which return a Dataset, whose variables can have incompatible axes such as power’s), squeeze is ignored with a warning.
connectivity_kwargs (dict, optional) – Extra keyword arguments for the measure (a Connectivity method), e.g. pairs for subset_pairwise_spectral_granger_prediction or n_components for canonical_coherency; passed to every requested measure. Transform settings do not go here (see **kwargs).
group_labels (sequence, optional) – One label per signal naming the group it belongs to; labels are scalars such as integers or area names ("CA1"), and missing values (None, NaN) are rejected. Required by the group measures (canonical_coherence, canonical_coherency, maximized_imaginary_coherency, multivariate_interaction_measure, blockwise_spectral_granger_prediction and maximized_imaginary_coherency_components); an error is raised if no requested measure takes it (labels are passed only to the measures that do).
frequency_range ((float, float), optional) – Inclusive frequency interval retained in the labeled result.
frequency_decimation (int, default=1) – Keep every Nth frequency bin after applying frequency_range.
frequency_bands (mapping of str to (float, float), optional) – Reduce the selected bins into named, inclusive bands. With frequency_reduction="mean", scores are averaged, complex measures use a complex vector mean, and coherence_phase uses a circular mean.
frequency_reduction ({"mean", "integral"}, default="mean") – Band reduction. Integration is restricted to power and cross_spectral_density, where it yields band power/covariance.
time_dim (hashable, optional) – DataArray dimension containing time samples. Common names such as "time" and "sample" are inferred automatically.
trial_dim (hashable, optional) – DataArray dimension containing trials or epochs. Required for a 3-D DataArray when its role cannot be inferred unambiguously.
signal_dim (hashable, optional) – DataArray dimension containing signals or channels. Common names such as "signal" and "channel" are inferred automatically.
**kwargs – Extra keyword arguments for the transform (Multitaper), e.g. time_halfbandwidth_product, n_tapers, time_window_step, taper_weighting (or fft_workers=-1 to parallelize the CPU FFT across all cores). Measure settings do not go here (see connectivity_kwargs).

Returns:

result – A plain single-quantity method returns a DataArray. Component-resolved and multi-quantity methods return a Dataset even when requested alone; multiple methods are merged into one Dataset without flattening their semantic dimensions.

Return type:

xarray.DataArray or xarray.Dataset

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> # Generate coupled oscillator data
>>> t = np.arange(0, 1, 1/500)  # 500 Hz, 1 second
>>> sig1 = np.sin(2*np.pi*10*t) + 0.1*rng.standard_normal(len(t))
>>> sig2 = np.sin(2*np.pi*10*t + np.pi/4) + 0.1*rng.standard_normal(len(t))
>>> # Shape (n_time, n_channels); a single trial of 2 signals. The 2-D form
>>> # is promoted to a single-trial 3-D array internally.
>>> data = np.stack([sig1, sig2], axis=-1)  # (500, 2)
>>>
>>> # Compute coherence
>>> coherence = multitaper_connectivity(
...     data, sampling_frequency=500,
...     method="coherence_magnitude",
...     signal_names=["Signal_1", "Signal_2"]
... )
>>> coherence.dims
('time', 'frequency', 'source', 'target')

>>> # Compute multiple measures
>>> measures = multitaper_connectivity(
...     data, sampling_frequency=500,
...     method=["coherence_magnitude", "imaginary_coherence"]
... )
>>> list(measures.data_vars)
['coherence_magnitude', 'imaginary_coherence']

>>> # An xarray.DataArray labels axes by dimension name and can supply the
>>> # sampling rate and channel labels itself (no sampling_frequency needed).
>>> import xarray as xr
>>> da = xr.DataArray(
...     data,
...     dims=("time", "channel"),
...     coords={"time": t, "channel": ["Signal_1", "Signal_2"]},
... )
>>> coherence = multitaper_connectivity(da, method="coherence_magnitude")
>>> coherence.coords["source"].values.tolist()
['Signal_1', 'Signal_2']

Notes

Uses multitaper spectral estimation for robust power spectral density estimation before computing connectivity measures. This provides better spectral estimates than single-taper methods, especially for short time series.

For directed measures (e.g. pairwise_spectral_granger_prediction) the source and target axes are oriented so that result.sel(source=a, target=b) is the influence from a to b, the same order as the underlying Connectivity arrays, where result[..., i, j] is i -> j. Signed undirected phase measures (coherence_phase, imaginary_coherency, phase_lag_index, weighted_phase_lag_index) are positive at sel(source=a, target=b) when a leads b.

Every variable has long_name and units attrs ("1" for dimensionless scores, "rad" for phase, "s" for delay; spectral densities are "(<units>)^2/Hz" when an input DataArray states its units, and "(<units>)^2" once integrated over a band). Non-index coordinates on an input DataArray’s signal dimension (e.g. region) are carried as source_<name>/target_<name>.

Real-valued results write with any NetCDF engine (booleans are stored as 0/1). Complex results (coherency, cross_spectral_density, canonical_coherency, and the global-coherence vectors) need an engine that stores complex data, e.g. result.to_netcdf("result.h5", engine="h5netcdf", invalid_netcdf=True), or netCDF4 >= 1.7 with engine="netcdf4", auto_complex=True (open with the same option).

The result records provenance as NetCDF-safe attributes so a saved file is self-describing:

mt_* – the Multitaper transform parameters.
measure and measure_kwargs_json – the measure name and a canonical, JSON-normalized representation of its keyword arguments.
arg_<key> / arg_<key>_json – each measure keyword argument individually for quick inspection; a scalar is stored as-is under arg_<key>, while a structured or non-finite value is stored as a JSON string under arg_<key>_json (parse with json.loads; measure_kwargs_json is the canonical record).
package, package_version, backend, expectation_type – software provenance.
input_attrs_json – a canonical, JSON-normalized record of attributes carried over from an input xarray.DataArray (e.g. subject or session metadata). Keeping the complete mapping in one record preserves arbitrary keys without collisions or invalid NetCDF attribute names.

JSON records are canonical for scalar, numpy, mapping, and sequence values. A value outside those kinds is recorded best-effort via its repr, which may embed a memory address and is therefore not guaranteed reproducible across runs.

References

[1]

Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096.

[2]

Percival, D. B., & Walden, A. T. (1993). Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques.