spectral_connectivity#

Test, Build, and Publish DOI Binder status PyPI version Anaconda-Server Badge Documentation Status codecov

What is spectral_connectivity?#

spectral_connectivity is a Python software package that computes multitaper spectral estimates and frequency-domain brain connectivity measures such as coherence, spectral granger causality, and the phase lag index using the multitaper Fourier transform. Although there are other Python packages that do this (see nitime and MNE-Python), spectral_connectivity has several differences:

  • it is designed to handle multiple time series at once

  • it caches frequently computed quantities such as the cross-spectral matrix and minimum-phase-decomposition, so that connectivity measures that use the same processing steps can be more quickly computed. Call Connectivity.clear_cache() to release these intermediates when memory is tight.

  • it decouples the time-frequency transform and the connectivity measures so that if you already have a preferred way of computing Fourier coefficients (i.e. from a wavelet transform), you can use that instead.

  • it implements the non-parametric version of the spectral granger causality in Python.

  • it implements the canonical coherence, which can efficiently summarize brain-area level coherences from multielectrode recordings.

  • easier user interface for the multitaper fourier transform

  • core transforms and connectivity calculations support GPU acceleration when cupy is installed and SPECTRAL_CONNECTIVITY_ENABLE_GPU=true is set before importing the package. Public results are returned as NumPy arrays.

Tutorials#

See the following notebooks for more information on how to use the package:

Usage Example#

The high-level multitaper_connectivity function runs the multitaper transform and returns a labeled xarray object:

from spectral_connectivity import multitaper_connectivity

coherence = multitaper_connectivity(
    time_series,
    sampling_frequency=sampling_frequency,
    method="coherence_magnitude",
    time_halfbandwidth_product=3,
)

The wrapper also provides coordinate-aware frequency cropping, decimation, and named-band reduction through frequency_range, frequency_decimation, and frequency_bands. Phase bands use a circular mean, complex results use a complex-vector mean, and integration is limited to power/cross-spectral density.

Use fourier_connectivity when coefficients were computed elsewhere. It accepts NumPy arrays or labeled DataArrays, preserves their time/frequency/signal coordinates, and uses the same result and provenance contract.

time_series may also be an xarray.DataArray. For DataArray inputs, dimension names define axis roles; positions do not. Common dimension names are inferred and transposed automatically; for domain-specific names, pass time_dim, trial_dim, and signal_dim explicitly. Ambiguous dimensions raise instead of falling back to axis position; when a single unrecognized dimension is left for the one remaining role, it is assigned by elimination and a warning names the assumed mapping. Numeric time coordinates are interpreted as elapsed seconds and numeric sample coordinates as sample numbers, and are used to label output window centers. When sampling_frequency is given it is checked against the time index; when it is omitted, a numeric elapsed-seconds time coordinate infers it (a sample index cannot, having no time scale). Inference also requires enough coordinate precision to resolve the rate reliably; pass sampling_frequency explicitly for low-precision or large-offset time coordinates. A 1-D index on the signal dimension is preserved—including its label type—as the result’s source and target coordinates unless signal_names is supplied. Signal labels must be unique, non-missing, NetCDF-compatible scalar strings, real numbers, datetimes, or timedeltas; integer labels must fit the signed 32-bit range for portable NetCDF3 serialization.

Datetime, timedelta, and object-valued time coordinates are not yet supported. Convert them to numeric elapsed seconds before calling multitaper_connectivity, for example:

da = da.assign_coords(time=(da.time - da.time[0]) / np.timedelta64(1, "s"))

datetime and timedelta signal labels remain valid.

A dask-backed DataArray is rejected; materialize it first with DataArray.compute() (or .load()) and pass the result.

For directed measures, result.sel(source="a", target="b") means influence from a to b (for phase measures, positive means a leads b). The directed-transfer-function family is available by name as an opt-in method. The lower-level Connectivity methods return plain arrays in one of two orders: for the Granger and directed-transfer-function families result[..., i, j] is j -> i, while for directed_phase_lag_index, phase_slope_index, delay, and group_delay a positive result[..., i, j] means i leads j. list_measures() reports each measure’s array_orientation, value range, units, and interpretation; see Connectivity Metric Ranges.

For finer control, use the Multitaper and Connectivity classes directly:

from spectral_connectivity import Multitaper, Connectivity

# Compute multitaper spectral estimate
m = Multitaper(
    time_series=signals,
    sampling_frequency=sampling_frequency,
    time_halfbandwidth_product=time_halfbandwidth_product,
    time_window_duration=0.060,
    time_window_step=0.060,
    start_time=time[0],
)

# Sets up computing connectivity measures/power from multitaper spectral estimate
# (`from_multitaper` is a backward-compatible alias for the transform-neutral
# `from_transform`.)
c = Connectivity.from_transform(m)

# Here are a couple of examples
power = c.power()  # spectral power
coherence = c.coherence_magnitude()
weighted_phase_lag_index = c.weighted_phase_lag_index()
canonical_coherence = c.canonical_coherence(brain_area_labels)
cacoh = c.canonical_coherency(brain_area_labels, n_components=2)
mic = c.maximized_imaginary_coherency_components(brain_area_labels, n_components=2)

The xarray wrappers retain nonstandard scientific shapes instead of flattening them: group-pair matrices have source_group/target_group; rich CaCoh/MIC results include connection, component, side, signal, filters, patterns, and group membership; delay has a candidate axis; and global coherence and group delay return multi-variable Datasets. Phase slope and group delay have no frequency dimension because their band has already been reduced.

ShortTimeFourierTransform, Welch, and MorletWavelet provide alternative spectral transforms with the same coefficient interface. Morlet output is positive-frequency-only and is therefore limited to functional measures; directed Wilson factorization requires a full two-sided FFT. For single-trial Morlet data, smoothing_time and smoothing_frequency collect a local time/frequency neighborhood; smoothing_kernel selects boxcar or Hann weights. padding_mode controls convolution boundaries, while edge_mode retains, masks, or trims estimates without full wavelet support. The strict valid_time_frequency mask is also exposed by the xarray wrapper. Welch resolves 1 / segment_duration Hz, so set segment_duration explicitly for electrophysiology data. Multitaper DPSS coefficients support uniform (historical), eigenvalue, or Thomson adaptive taper weighting through taper_weighting.

For uncertainty estimates, Connectivity.jackknife(method) recomputes a real-valued measure while leaving out each trial/taper observation and returns a bias-corrected estimate, standard error, and confidence interval, applying an automatic variance-stabilizing transformation (log for power, atanh(sqrt(.)) for magnitude-squared coherence, and circular for phase).

Citation#

For citation, please use the following:

Denovellis, E.L., Myroshnychenko, M., Sarmashghi, M., and Stephen, E.P. (2022). Spectral Connectivity: a python package for computing multitaper spectral estimates and frequency-domain brain connectivity measures on the CPU and GPU. JOSS 7, 4840. 10.21105/joss.04840.

Implemented Measures#

Functional

  1. coherency

  2. cross_spectral_density

  3. coherence_magnitude and coherence_phase

  4. imaginary_coherence and signed imaginary_coherency

  5. partial_coherence

  6. canonical_coherence and exact complex canonical_coherency (CaCoh)

  7. maximized_imaginary_coherency (score-only or component-resolved) and multivariate_interaction_measure

  8. phase_locking_value and corrected_imaginary_phase_locking_value

  9. phase_lag_index, directed_phase_lag_index, and weighted_phase_lag_index

  10. debiased_squared_phase_lag_index

  11. debiased_squared_weighted_phase_lag_index

  12. pairwise_phase_consistency

  13. global_coherence

Directed

  1. directed_transfer_function

  2. directed_coherence

  3. partial_directed_coherence

  4. generalized_partial_directed_coherence

  5. direct_directed_transfer_function

  6. group_delay

  7. pairwise_spectral_granger_prediction

  8. conditional_spectral_granger_prediction

  9. blockwise_spectral_granger_prediction

  10. time_reversed_spectral_granger_prediction

Package Dependencies#

spectral_connectivity requires:

  • python

  • numpy

  • matplotlib

  • scipy

  • xarray

See the repository’s pyproject.toml for the authoritative dependency list.

Installation#

pip install spectral_connectivity

or

conda install -c edeno spectral_connectivity

Developer Installation#

If you want to make contributions to this library, please use this installation.

  1. Install miniconda (or anaconda) if it isn’t already installed. Type into bash (or install from the anaconda website):

wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
bash miniconda.sh -b -p $HOME/miniconda
export PATH="$HOME/miniconda/bin:$PATH"
hash -r

  1. Clone the repository to your local machine (.../spectral_connectivity) and install the anaconda environment for the repository. Type into bash:

conda env create -f environment.yml
conda activate spectral_connectivity
pip install -e .

Recent publications and pre-prints that used this software#