spectral_connectivity

Tutorials | Documentation | Usage Example | Installation | Developer Installation | Contributing | License | Citation

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.

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

• Tutorial
• Usage Examples
• More Usage Examples

The quickest way to get connectivity measures is the high-level multitaper_connectivity function, which runs the multitaper transform and returns labeled xarray objects (dimensions such as time, frequency, source, target):

from spectral_connectivity import multitaper_connectivity

# time_series has shape (n_time_samples, n_trials, n_signals)
coherence = multitaper_connectivity(
    time_series,
    sampling_frequency=sampling_frequency,
    method="coherence_magnitude",
    time_halfbandwidth_product=3,
)

# Ask for several measures at once (returns an xarray.Dataset)
measures = multitaper_connectivity(
    time_series,
    sampling_frequency=sampling_frequency,
    method=["coherence_magnitude", "imaginary_coherence", "phase_locking_value"],
)

Not sure which method names are valid? list_measures() enumerates every supported measure with its output category and a one-line description, and an unknown method name raises an error that suggests the closest matches:

from spectral_connectivity import list_measures

# Every valid method name, e.g. filter to the directed measures
directed = [m.name for m in list_measures(directed=True)]

See the Cookbook for short, copy-pasteable recipes covering functional and directed connectivity, reading the labeled output, frequency bands, and using your own Fourier coefficients.

Crop, decimate, or aggregate the frequency coordinate without losing labels:

band_connectivity = multitaper_connectivity(
    time_series,
    sampling_frequency=sampling_frequency,
    method="coherence_magnitude",
    frequency_range=(4, 100),
    frequency_bands={"theta": (4, 8), "alpha": (8, 13), "beta": (13, 30)},
)

Named-band means are arithmetic means of the computed score; complex measures use a complex-vector mean and coherence_phase uses a circular mean. Band integration is restricted to power and cross-spectral density, where it has the physical interpretation of band power/covariance. The same reduction is available for an already-computed result via frequency_band_reduce(result, bands, reduction="mean").

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 output's source and target coordinates unless signal_names is passed explicitly. 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.

Not sure what time_halfbandwidth_product, number of tapers, or window duration to use? suggest_parameters picks reasonable values from your sampling rate, signal duration, and desired frequency resolution:

from spectral_connectivity import suggest_parameters

params = suggest_parameters(
    sampling_frequency=1000,
    signal_duration=2.0,  # seconds
    desired_freq_resolution=4.0,  # Hz
)
print(params)  # -> time_halfbandwidth_product, time_window_duration, n_tapers, ...

See also estimate_frequency_resolution, estimate_n_tapers, and Multitaper.summarize_parameters() for related helpers.

For finer control (or if you already have Fourier coefficients from another method), 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],
    # fft_workers=-1,  # optional: parallelize the CPU FFT across all cores
)

# Sets up computing connectivity measures/power from multitaper spectral estimate
# (`from_transform` is the transform-neutral spelling; `from_multitaper` is a
# backward-compatible alias for it.)
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)
# Exact Vidaurre CaCoh is complex and component-resolved. Rich MIC uses the
# same result schema; both include filters, patterns, and group membership.
cacoh = c.canonical_coherency(brain_area_labels, n_components=2)
mic = c.maximized_imaginary_coherency_components(brain_area_labels, n_components=2)

The xarray interfaces preserve these non-pairwise contracts as labeled axes: group measures use source_group/target_group; CaCoh and rich MIC use connection, component, side, and signal plus group-membership metadata; delay exposes a candidate axis; global coherence returns score and vector variables; and group delay returns delay, slope, and correlation variables. Phase slope and group delay are already frequency-reduced, so select their band with connectivity_kwargs={"frequencies_of_interest": (...)}.

The same connectivity core accepts Hann short-time Fourier, Welch, and Morlet wavelet transforms. Morlet coefficients are one-sided, so they support functional measures but deliberately reject Wilson-factorized directed measures:

from spectral_connectivity import (
    Connectivity,
    MorletWavelet,
    ShortTimeFourierTransform,
    Welch,
)

stft = ShortTimeFourierTransform(time_series, sampling_frequency, time_window_duration=1)
welch = Welch(time_series, sampling_frequency, segment_duration=1)
# For a single continuous trial, collect a local time/frequency neighborhood so
# normalized measures are estimated over multiple observations. Invalid wavelet
# edges can be retained, marked NaN, or trimmed.
wavelet = MorletWavelet(
    time_series,
    sampling_frequency,
    frequencies=[4, 8, 16, 32],
    smoothing_time=0.5,
    smoothing_frequency=3,
    smoothing_kernel="hann",
    padding_mode="reflect",
    edge_mode="nan",
)

stft_coherence = Connectivity.from_transform(stft).coherence_magnitude()

Which transform? Use Multitaper for stationary spectra with controlled variance, ShortTimeFourierTransform/Welch for a simple single-window or segment-averaged estimate, and MorletWavelet for time-resolved analysis at specific frequencies. Welch's frequency resolution is 1 / segment_duration Hz, so set segment_duration explicitly for electrophysiology data. MorletWavelet.valid_time_frequency identifies bins with full in-record support; the xarray wrapper carries it as a two-dimensional coordinate.

For DPSS transforms, taper_weighting="uniform" preserves the historical behavior; "eigen" weights by concentration ratio and "adaptive" applies Thomson's frequency- and signal-specific iterative weights (useful when high spectral dynamic range or line noise makes some tapers systematically leakier).

Already have full two-sided FFT coefficients? fourier_connectivity accepts NumPy arrays or labeled DataArrays in the core (time, trial, taper, frequency, signal) layout (with documented shorter forms), preserving frequency, time, and signal coordinates.

Uncertainty for real-valued measures is available through Connectivity.jackknife(...), with automatic variance-stabilizing transformations (log for power, atanh(sqrt(.)) for magnitude-squared coherence, and circular for phase).

See the documentation on ReadTheDocs.

For a canonical reference of connectivity metric value ranges, see Connectivity Metric Ranges.

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

spectral_connectivity requires:

• python
• numpy
• matplotlib
• scipy
• xarray

See pyproject.toml for the authoritative list of dependencies and their minimum versions.

spectral_connectivity supports GPU acceleration using CuPy, which can accelerate large workloads when a suitable CUDA device is available.

There are three ways to enable GPU acceleration:

Option 1: Environment Variable (Shell)

export SPECTRAL_CONNECTIVITY_ENABLE_GPU=true
python your_script.py

Option 2: Environment Variable (Python Script)

import os

# IMPORTANT: Must set BEFORE importing spectral_connectivity
# (Python loads modules once; changing the variable after import has no effect)
os.environ["SPECTRAL_CONNECTIVITY_ENABLE_GPU"] = "true"

from spectral_connectivity import Multitaper, Connectivity

# Verify GPU is active
import spectral_connectivity as sc

backend = sc.get_compute_backend()
print(backend["message"])
# Should print: "Using GPU backend with CuPy on <your GPU name>"

Option 3: Environment Variable (Jupyter Notebook)

# In first cell (before any imports):
%env SPECTRAL_CONNECTIVITY_ENABLE_GPU=true

# In second cell:
from spectral_connectivity import Multitaper, Connectivity
import spectral_connectivity as sc

# Verify GPU is active
backend = sc.get_compute_backend()
print(f"Backend: {backend['backend']}")
print(f"Device: {backend['device_name']}")
# Should show: Backend: gpu, Device: <your GPU name>

# Note: If you already imported spectral_connectivity before setting the
# environment variable, you must restart your kernel for changes to take effect:
# Kernel → Restart & Clear Output, then run cells again

Recommended (conda - auto-detects CUDA version):

conda install -c conda-forge cupy

Recommended pip installation (CUDA 12):

pip install "spectral-connectivity[gpu]"

The package's GPU extra installs cupy-cuda12x>=13.0. For another CUDA version, install the corresponding CuPy 13+ distribution directly instead of the extra.

# Check your CUDA version first: nvidia-smi
pip install "cupy-cuda11x>=13.0"  # CUDA 11.x
pip install "cupy-cuda12x>=13.0"  # CUDA 12.x

See CuPy Installation Guide for detailed instructions and GPU-specific requirements.

Use get_compute_backend() to check if GPU acceleration is enabled:

import spectral_connectivity as sc

backend = sc.get_compute_backend()
print(backend["message"])
# Example output (GPU enabled):
# "Using GPU backend with CuPy on NVIDIA Tesla V100-SXM2-16GB."

# Or if GPU not available:
# "Using CPU backend with NumPy. To enable GPU acceleration:
#   1. Install CuPy: 'conda install -c conda-forge cupy' or
#      'pip install "spectral-connectivity[gpu]"'
#   2. Set environment variable SPECTRAL_CONNECTIVITY_ENABLE_GPU='true' before importing
# See documentation for detailed setup instructions."

# Check all details
for key, value in backend.items():
    print(f"{key}: {value}")

Output fields:

• backend: Either "cpu" or "gpu"
• gpu_enabled: Whether GPU was requested via environment variable
• gpu_available: Whether CuPy is installed and importable
• device_name: Name of compute device (e.g., "CPU" or "GPU (Compute Capability 7.5)")
• message: Human-readable explanation of current configuration

Issue: "GPU support was requested but CuPy is not installed"

Solution: Install CuPy as shown above, ensuring the CUDA version matches your system.

Issue: GPU not being used even after setting environment variable

Possible causes:

1. Environment variable set after importing spectral_connectivity
   • Solution: Set SPECTRAL_CONNECTIVITY_ENABLE_GPU=true before any imports
   • In scripts: Move the os.environ[...] line to the very top, before all spectral_connectivity imports
   • In notebooks: Restart kernel (Kernel → Restart & Clear Output) and set variable in first cell
2. CuPy not installed or CUDA version mismatch
   • Solution: Run python -c "import cupy; print(cupy.__version__)" to verify installation
3. CUDA not available on system
   • Solution: Check CUDA installation with nvidia-smi
4. Environment variable set in script but after import statement
   • Solution: Ensure os.environ['SPECTRAL_CONNECTIVITY_ENABLE_GPU'] = 'true' appears before from spectral_connectivity import ...

Issue: Out of memory errors on GPU

Solution: Use smaller batch sizes or switch back to CPU for very large datasets:

# Remove or unset the environment variable
os.environ.pop("SPECTRAL_CONNECTIVITY_ENABLE_GPU", None)

Issue: Need to check GPU vs CPU performance

Use get_compute_backend() to verify which backend is active, then compare timing:

import time
import spectral_connectivity as sc

backend = sc.get_compute_backend()
print(f"Running on: {backend['backend']}")

start = time.time()
# Your spectral connectivity code here
elapsed = time.time() - start
print(f"Elapsed time: {elapsed:.2f}s")

GPU acceleration is most beneficial for:

• Large datasets (many signals, long recordings, or many trials)
• High frequency resolution (small time windows, many tapers)
• Computing multiple connectivity measures from the same data

For small datasets (< 10 signals, < 1000 time points), CPU may be faster due to GPU transfer overhead.

pip install spectral_connectivity

or

conda install -c edeno spectral_connectivity

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.anaconda.com/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 create the development environment. Type into bash:

# uv: creates .venv from the committed uv.lock with the dev tools installed
uv sync
uv run pytest                     # any command runs in that environment
uvx pre-commit install            # optional: run ruff, codespell and mypy on each commit

# Or conda
conda env create -f environment.yml
conda activate spectral_connectivity
pip install -e .[dev]

uv.lock pins the development environment; uv lock refreshes it after a dependency change, and CI checks that it is current.

This package uses dynamic versioning with Hatch based on git tags. The version is automatically determined from the repository state:

• Tagged releases: 1.2.0
• Development versions: 1.2.0.dev5+g1a2b3c4 (5 commits since tag + git hash)

Releases are published automatically by CI: pushing a v* tag runs the Test, Build, and Publish workflow, which tests, builds, attests, and publishes to PyPI via trusted publishing (after approval in the protected pypi environment). Do not run twine upload by hand — that bypasses the tests, attestations, and approval gate. See CONTRIBUTING.md for the full process and the required one-time setup.

# Update CHANGELOG.md, then tag and push -- CI does the rest.
git tag -a v1.2.0 -m "v1.2.0"
git push origin v1.2.0

Conda packages are published separately and manually:

conda build conda-recipe/ --output-folder ./conda-builds
anaconda upload ./conda-builds/noarch/spectral_connectivity-*.tar.bz2

The version number is automatically extracted from the git tag (without the 'v' prefix).

This package is also available on conda via the edeno channel:

conda install -c edeno spectral_connectivity

Not yet on conda-forge? Help us get there! If you'd like this package on conda-forge for easier installation, please:

• 👍 React to this issue requesting conda-forge support
• Or volunteer to help maintain the conda-forge feedstock

We welcome contributions to spectral_connectivity! Please see our Contributing Guidelines for details on:

• How to report bugs and request features
• Development workflow and coding standards
• Testing requirements
• Code review process

For questions or discussions, please open an issue on GitHub.

This project is licensed under the GPL-3.0 License - see the LICENSE file for details.

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.

• Detection of Directed Connectivities in Dynamic Systems for Different Excitation Signals using Spectral Granger Causality https://doi.org/10.1007/978-3-662-58485-9_11
• Network Path Convergence Shapes Low-Level Processing in the Visual Cortex https://doi.org/10.3389/fnsys.2021.645709
• Subthalamic–Cortical Network Reorganization during Parkinson's Tremor https://doi.org/10.1523/JNEUROSCI.0854-21.2021
• Unifying Pairwise Interactions in Complex Dynamics https://doi.org/10.48550/arXiv.2201.11941
• Phencyclidine-induced psychosis causes hypersynchronization and disruption of connectivity within prefrontal-hippocampal circuits that is rescued by antipsychotic drugs https://doi.org/10.1101/2021.02.03.429582
• The cerebellum regulates fear extinction through thalamo-prefrontal cortex interactions in male mice https://doi.org/10.1038/s41467-023-36943-w