EEG-Pype: An accessible MNE-Python pipeline with graphical user interface for preprocessing and analysis of resting-state electroencephalography data.

Code ↔ Paper

8 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 8 matches

[1] § Design and implementation › Workflow for preprocessing ↔ src/eeg_processing_script.py, lines 2314–2406 · score 0.97 · boundary element model, noise covariance matrix, unit noise gain, weight normalization, source space, spatial filter
[2] § Design and implementation › Quantitative analysis ↔ src/eeg_quantitative_analysis.py, lines 1620–1723 · score 0.85 · concatenated broadband, band power, Spectral variability, relative power, Window lengths, coefficient
[3] § Design and implementation › Workflow for preprocessing ↔ src/eeg_processing_script.py, lines 2999–3049 · score 0.79 · custom parser, configurable header, raw EEG, CNT, FIF, generic
[4] § Design and implementation › Quantitative analysis ↔ src/eeg_quantitative_analysis.py, lines 1941–2070 · score 0.70 · minimum spanning tree, connectivity matrices, graph, weighted, network, efficiency
[5] § Results › 5. Deriving quantitative EEG measures ↔ src/eeg_quantitative_analysis.py, lines 3233–3313 · score 0.68 · approximate entropy, sample entropy, Power spectral density, peak frequency, relative power, Multitaper
[6] § Design and implementation › Quantitative analysis ↔ src/eeg_quantitative_analysis.py, lines 2283–2330 · score 0.65 · joint permutation entropy, volume conduction, JPE
[7] § Design and implementation › Workflow for preprocessing ↔ src/eeg_processing_script.py, lines 3427–3482 · score 0.62 · channel noise, ICA components, heart, ICALabel, muscle, eye
[8] § Design and implementation › Quantitative analysis ↔ src/eeg_quantitative_analysis.py, lines 2213–2254 · score 0.56 · embedding dimension, Permutation entropy

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 3,945 lines · 163 KB · Apache-2.0 · 5 matches

"""@authors: Herman van Dellen en Yorben Lodema."""
import itertools
import json
import logging
import math
import os
import pickle
import random
import re
import time
from collections import defaultdict
from datetime import datetime
from multiprocessing import Pool, cpu_count
from typing import Optional
import mne
import networkx as nx
import numpy as np
import pandas as pd
import psutil
from antropy import lziv_complexity, sample_entropy
from scipy import signal
from scipy.signal import hilbert
from scipy.integrate import trapezoid
from scipy.sparse.csgraph import minimum_spanning_tree
import FreeSimpleGUI as sg
EEG_version = "v4.6.0"
logger = logging.getLogger(__name__)
# Configuration
FOLDER_EXTENSION = "bdf" # Change this to match your folder extension (e.g., 'bdf', 'edf', etc.)
MAX_MEMORY_PERCENT = 70 # Maximum memory usage percentage
MIN_WINDOW_SIZE = 100 # Minimum window size for spectral variability in ms
# Be careful, option to change frequency bands (both those recognized in the epoch file names
# and bands used for power and spectral variability calculations. Don't change the format. You can add additional
# bands in the same format. Connectivity measures will still be calculated on epochs with bands not
# listed here though that might affect the naming of the output columns in the Excel output.
# It is probably advisable to limit the broadband range to something like 20-30 Hz to make sure the
# total power used in relative power and spectral variability contain less EMG noise.
FREQUENCY_BANDS = {
"delta": {"pattern": r"0\.5-4\.0|delta", "range": (0.5, 4.0)},
"theta": {"pattern": r"4\.0-8\.0|theta", "range": (4.0, 8.0)},
"alpha": {"pattern": r"8\.0-13\.0|alpha", "range": (8.0, 13.0)},
"alpha1": {"pattern": r"8\.0-10\.0|alpha1", "range": (8.0, 10.0)},
"alpha2": {"pattern": r"10\.0-13\.0|alpha2", "range": (10.0, 13.0)},
"beta": {"pattern": r"13\.0-30\.0|beta", "range": (13.0, 30.0)},
"beta1": {"pattern": r"13\.0-20\.0|beta1", "range": (13.0, 20.0)},
"beta2": {"pattern": r"20\.0-30\.0|beta2", "range": (20.0, 30.0)},
"broadband": {"pattern": r"0\.5-47|broadband", "range": (0.5, 47.0)},
}
def validate_frequency_bands():
"""Validate FREQUENCY_BANDS configuration."""
if not FREQUENCY_BANDS:
msg = "FREQUENCY_BANDS dictionary is empty"
raise ValueError(msg)
for band_name, band_info in FREQUENCY_BANDS.items():
if "pattern" not in band_info or "range" not in band_info:
msg = f"Band {band_name} missing required keys (pattern, range)"
raise ValueError(msg)
fmin, fmax = band_info["range"]
if not (isinstance(fmin, (int, float)) and isinstance(fmax, (int, float))):
msg = f"Band {band_name} range values must be numeric"
raise TypeError(msg)
if fmin >= fmax:
msg = f"Band {band_name} minimum frequency must be less than maximum"
raise ValueError(msg)
if not isinstance(band_info["pattern"], str):
msg = f"Band {band_name} pattern must be a string"
raise TypeError(msg)
# ---------------------------------------------------------------------------
# Preprocessing pass-band discovery
#
# The preprocessing script (EEG-Pype) writes a batch config as
# <batch_output_subdirectory>/<batch_name>.pkl, holding the band cut-offs under
# tuple keys such as ("cut_off_frequency", "broadband_low"/"broadband_high"),
# and a per-recording sidecar <file_output_subdirectory>/EXPORT_METADATA_NAME
# containing filter metadata recorded from raw.info (base_highpass_hz /
# base_lowpass_hz). For this release, the identity-matched batch config is
# authoritative; sidecars are used only as a fallback when no valid config can
# be found.
#
# Knowing this pass-band matters for relative power: the denominator should not
# extend beyond the range that actually contains signal. Nyquist is a far weaker
# bound (e.g. 125 Hz at fs=250) than the real 0.5-47 Hz pass-band, so validating
# against Nyquist alone lets a user silently include filter roll-off in the
# total-power reference.
# ---------------------------------------------------------------------------
EXPORT_METADATA_NAME = "eeg_pype_export_metadata.json"
PASSBAND_HINT_DEFAULT = "(denominator for rel. power & spectral variability)"
PASSBAND_SEARCH_PARENT_LEVELS = 3 # how far up from the selected folder to look
PASSBAND_MATCH_ABS_TOL = 1e-6
class _RestrictedUnpickler(pickle.Unpickler):
"""Unpickler that refuses to import anything outside a small allowlist.
The batch config is a plain dict of scalars, strings and lists, so nothing
exotic is required. Restricting find_class means a corrupted or untrusted
.pkl in the analysis folder cannot execute arbitrary imports while we are
only trying to read two floats out of it.
"""
_ALLOWED = {
"builtins": {
"dict", "list", "tuple", "set", "frozenset",
"str", "int", "float", "bool", "complex", "bytes",
},
"collections": {"OrderedDict", "defaultdict"},
"numpy": {"dtype", "ndarray"},
"numpy.core.multiarray": {"_reconstruct", "scalar"},
"numpy._core.multiarray": {"_reconstruct", "scalar"},
}
def find_class(self, module, name):
if name in self._ALLOWED.get(module, ()):
return super().find_class(module, name)
msg = f"Blocked unpickling of {module}.{name}"
raise pickle.UnpicklingError(msg)
def _coerce_positive_float(value):
"""Return value as a finite float, or None if it is not usable."""
if isinstance(value, bool) or value is None:
return None
try:
out = float(value)
except (TypeError, ValueError):
return None
if not math.isfinite(out) or out < 0:
return None
return out
def _passband_from_batch_config(path):
"""Read and validate broadband cut-offs from a preprocessing batch config."""
with open(path, "rb") as f:
cfg = _RestrictedUnpickler(f).load()
if not isinstance(cfg, dict):
return None
# EEG-Pype creates:
#
# <batch_name>/<batch_name>.pkl
#
# Require the filename and the batch name stored inside the config to agree.
# This prevents an unrelated or copied .pkl from being selected silently.
expected_batch_name = os.path.splitext(os.path.basename(path))[0]
configured_batch_name = cfg.get("batch_name")
if (
configured_batch_name is not None
and str(configured_batch_name) != expected_batch_name
):
logger.warning(
f"Ignoring batch config {path}: internal batch name "
f"{configured_batch_name!r} does not match filename "
f"{expected_batch_name!r}."
)
return None
# Also check the basename of the stored batch output directory when this
# information is available. Comparing only the basename keeps moved batch
# folders usable while still detecting an internally inconsistent config.
configured_batch_directory = cfg.get("batch_output_subdirectory")
if configured_batch_directory:
configured_directory_name = os.path.basename(
os.path.normpath(str(configured_batch_directory))
)
if configured_directory_name != expected_batch_name:
logger.warning(
f"Ignoring batch config {path}: stored output directory "
f"{configured_directory_name!r} does not match batch name "
f"{expected_batch_name!r}."
)
return None
low = _coerce_positive_float(
cfg.get(("cut_off_frequency", "broadband_low"))
)
high = _coerce_positive_float(
cfg.get(("cut_off_frequency", "broadband_high"))
)
if low is None or high is None or low >= high:
return None
return {
"low": low,
"high": high,
"origin": "config",
"detail": (
f"broadband filter cut-offs in "
f"{os.path.basename(path)}"
),
}
def _passband_from_export_metadata(path):
"""Read the applied filter from one EEG-Pype JSON sidecar."""
with open(path, encoding="utf-8") as f:
meta = json.load(f)
if not isinstance(meta, dict):
return None
# Ignore JSON belonging to another application or schema.
schema_name = meta.get("schema_name")
if schema_name not in (None, "EEG-Pype export metadata"):
return None
# An in-progress or failed recording should not determine the analysis
# pass-band. Missing status is accepted for backwards compatibility.
status = meta.get("status")
if status not in (None, "complete"):
logger.debug(
f"Ignoring pass-band metadata from {path}: "
f"recording status is {status!r}."
)
return None
output = meta.get("output")
if not isinstance(output, dict):
return None
low = _coerce_positive_float(
output.get("base_highpass_hz")
)
high = _coerce_positive_float(
output.get("base_lowpass_hz")
)
if low is None or high is None or low >= high:
return None
return {
"low": low,
"high": high,
"origin": "sidecar",
"detail": f"filter recorded in {os.path.basename(path)}",
}
def _iter_candidate_config_paths(folder_path):
"""Yield only identity-matched batch configs, nearest folder first.
EEG-Pype names the batch configuration after its containing batch folder:
<batch_name>/<batch_name>.pkl
Therefore, do not scan arbitrary .pkl files. At every candidate directory,
only inspect the .pkl whose basename exactly matches that directory.
"""
current = os.path.abspath(folder_path)
for _ in range(PASSBAND_SEARCH_PARENT_LEVELS + 1):
directory_name = os.path.basename(
os.path.normpath(current)
)
if directory_name:
candidate = os.path.join(
current,
f"{directory_name}.pkl",
)
if os.path.isfile(candidate):
yield candidate
parent = os.path.dirname(current)
if parent == current:
return
current = parent
def _iter_candidate_sidecar_paths(folder_path):
"""Yield sidecars in the selected folder and its direct subfolders."""
root = os.path.abspath(folder_path)
yielded_paths = set()
direct = os.path.join(
root,
EXPORT_METADATA_NAME,
)
if os.path.isfile(direct):
yielded_paths.add(direct)
yield direct
try:
entries = sorted(os.listdir(root))
except OSError:
return
for entry in entries:
candidate = os.path.join(
root,
entry,
EXPORT_METADATA_NAME,
)
if (
os.path.isfile(candidate)
and candidate not in yielded_paths
):
yielded_paths.add(candidate)
yield candidate
def _passbands_match(first, second):
"""Return whether two discovered pass-bands are numerically equal."""
return (
math.isclose(
first["low"],
second["low"],
rel_tol=0.0,
abs_tol=PASSBAND_MATCH_ABS_TOL,
)
and math.isclose(
first["high"],
second["high"],
rel_tol=0.0,
abs_tol=PASSBAND_MATCH_ABS_TOL,
)
)
def _discover_consistent_sidecar_passband(folder_path):
"""Return one pass-band only when all usable sidecars agree.
Returns
-------
tuple
``(result, conflict_found)``.
- result is a pass-band dict or None.
- conflict_found is True when valid sidecars contained different
pass-bands. Sidecars are fallback metadata, so a conflict matters
only when no valid identity-matched batch config is available.
"""
sidecar_results = []
for path in _iter_candidate_sidecar_paths(folder_path):
try:
result = _passband_from_export_metadata(path)
except (
OSError,
json.JSONDecodeError,
TypeError,
ValueError,
) as exc:
logger.debug(
f"Could not read pass-band from {path}: {exc}"
)
continue
if result is not None:
sidecar_results.append(
(path, result)
)
if not sidecar_results:
return None, False
reference_path, reference_result = sidecar_results[0]
conflicting_results = [
(path, result)
for path, result in sidecar_results[1:]
if not _passbands_match(
reference_result,
result,
)
]
if conflicting_results:
ranges = {
(
round(result["low"], 9),
round(result["high"], 9),
)
for _, result in sidecar_results
}
formatted_ranges = ", ".join(
f"{low:g}-{high:g} Hz"
for low, high in sorted(ranges)
)
logger.warning(
"No valid identity-matched batch config was found, and "
"the available export metadata sidecars contain different "
f"pass-bands: {formatted_ranges}. Automatic pass-band "
"discovery has been disabled."
)
return None, True
return {
"low": reference_result["low"],
"high": reference_result["high"],
"origin": "sidecar",
"source": reference_path,
"detail": (
f"consistent filters in "
f"{len(sidecar_results)} usable export metadata "
f"sidecar(s)"
),
}, False
def discover_preprocessing_passband(folder_path):
"""Discover the preprocessing pass-band for the selected data folder.
Precedence for this release is deliberately simple and conservative:
1. Use the nearest identity-matched <batch_name>/<batch_name>.pkl config.
This is the authoritative source for the configured broadband export
range.
2. Only when no valid matching config is found, inspect completed
recording-level JSON sidecars as a fallback.
3. Fallback sidecars must agree with one another. If they disagree,
automatic discovery is disabled and the user-entered range is retained.
The config and sidecars are not compared directly because the current JSON
fields are derived from raw.info and do not yet unambiguously describe the
pass-band of every exported broadband epoch file.
"""
if not folder_path or not os.path.isdir(folder_path):
return None
# First try the nearest identity-matched batch config. As soon as a valid
# config is found, return it without scanning recording sidecars. Besides
# reflecting the intended precedence, this avoids potentially opening
# hundreds of JSON files on the GUI thread for a normal batch selection.
for path in _iter_candidate_config_paths(folder_path):
try:
result = _passband_from_batch_config(path)
except (
OSError,
pickle.UnpicklingError,
EOFError,
AttributeError,
TypeError,
ValueError,
) as exc:
logger.debug(
f"Could not read pass-band from {path}: {exc}"
)
continue
if result is not None:
result["source"] = path
logger.info(
f"Detected preprocessing pass-band "
f"{result['low']:g}-{result['high']:g} Hz "
f"from {result['detail']}."
)
return result
# No valid identity-matched config was found. Use sidecars only as a
# fallback, and only when all usable sidecars agree.
sidecar_result, sidecar_conflict = (
_discover_consistent_sidecar_passband(folder_path)
)
if sidecar_conflict:
return None
if sidecar_result is not None:
logger.info(
f"Detected preprocessing pass-band "
f"{sidecar_result['low']:g}-{sidecar_result['high']:g} Hz "
f"from fallback {sidecar_result['detail']}."
)
return sidecar_result
logger.info(
"No unambiguous EEG-Pype preprocessing pass-band was found near "
"the selected folder. The relative-power reference range will be "
"bounded by the PSD and Nyquist limits only."
)
return None
def effective_reference_ceiling(power_fs, psd_method, passband=None):
"""Return the highest usable top edge for the relative-power reference range.
Combines the three independent ceilings: Nyquist, the multitaper PSD cut-off
(the multitaper path computes up to 60 Hz), and the pass-band of the
preprocessed data when it could be discovered. Returns (ceiling, reason).
"""
ceiling = power_fs / 2.0
reason = f"Nyquist ({ceiling:g} Hz)"
if psd_method == "multitaper":
mt_ceiling = min(60.0, power_fs / 2.0)
if mt_ceiling < ceiling:
ceiling = mt_ceiling
reason = f"multitaper PSD cut-off ({ceiling:g} Hz)"
if passband is not None and passband["high"] < ceiling:
ceiling = passband["high"]
reason = f"preprocessing low-pass ({ceiling:g} Hz)"
return ceiling, reason
BATCH_SIZE = 10 # Number of subjects to process in parallel
DEFAULT_THREADS = max(1, int(cpu_count() * 0.7))
class MemoryMonitor:
@staticmethod
def get_memory_usage():
"""Get current memory usage percentage."""
return psutil.Process().memory_percent()
@staticmethod
def check_memory():
"""Check if memory usage is too high."""
return MemoryMonitor.get_memory_usage() > MAX_MEMORY_PERCENT
@staticmethod
def check_concatenation_safety(data_size, num_epochs):
"""Check if concatenation is likely to exceed memory limits."""
try:
# Get system memory info
system_memory = psutil.virtual_memory()
available_memory = system_memory.available
# Calculate estimated memory needed (add 20% buffer)
estimated_memory = data_size * num_epochs * 1.2
# Check if we'll exceed the threshold
memory_threshold = (available_memory * MAX_MEMORY_PERCENT) / 100
if estimated_memory > memory_threshold:
logger.warning(
f"Concatenation may exceed memory limits. "
f"Estimated need: {estimated_memory / 1e9:.2f}GB, "
f"Available: {memory_threshold / 1e9:.2f}GB"
)
return False
return True
except Exception:
logger.exception("Error checking memory for concatenation")
return False
def setup_logging(folder_path):
"""Set up logging for the current analysis run."""
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
log_filename = os.path.join(folder_path, f"eeg_analysis_{timestamp}.log")
# Clear any existing handlers
logging.getLogger().handlers = []
# Configure logging with both file and console output
logging.basicConfig(
level=logging.INFO,
format="%(asctime)s - %(levelname)s - %(message)s",
handlers=[logging.FileHandler(log_filename, mode="w"), logging.StreamHandler()],
)
# Test the logging setup
logger.info("Logging initialized")
logger.info(f"Log file created at: {log_filename}")
return log_filename
def compute_resolution_readout(values):
"""Build a human-readable frequency-resolution summary for the selected PSD method.
Uses the sampling rate and the manually entered epoch length to derive bin
spacing, and for multitaper also the resolution bandwidth and implied taper
count. Never raises — returns a hint string if inputs are incomplete.
"""
method = values.get("-PSD_METHOD-", "Multitaper").lower()
try:
fs = float(values.get("-POWER_FS-", ""))
if fs <= 0:
raise ValueError
except (ValueError, TypeError):
return "Frequency resolution: enter a valid sampling rate"
try:
epoch_sec = float(values.get("-EPOCH_LENGTH-", ""))
if epoch_sec <= 0:
raise ValueError
n_samples = int(round(epoch_sec * fs))
except (ValueError, TypeError):
epoch_sec = None
n_samples = None
if method == "welch":
try:
win_ms = float(values.get("-WELCH_WINDOW-", ""))
if win_ms <= 0:
raise ValueError
except (ValueError, TypeError):
return "Frequency resolution: enter a valid Welch window"
return f"Welch bin spacing: {1000.0 / win_ms:.3f} Hz (= 1000 / {win_ms:.0f} ms)"
if method == "fft":
if epoch_sec is None:
return "FFT bin spacing: enter epoch length to compute"
return (f"FFT bin spacing: {1.0 / epoch_sec:.3f} Hz "
f"(= 1 / {epoch_sec:g} s, fixed by epoch length)")
# Multitaper
if epoch_sec is None:
return "Multitaper: enter epoch length to compute resolution"
bin_hz = fs / n_samples
bw_raw = values.get("-MT_BANDWIDTH-", "").strip()
if bw_raw == "" or bw_raw.lower() == "auto":
bandwidth = 8.0 * fs / n_samples # MNE default
bw_label = f"{bandwidth:.3f} Hz (auto)"
else:
try:
bandwidth = float(bw_raw)
if bandwidth <= 0:
raise ValueError
bw_label = f"{bandwidth:.3f} Hz"
except (ValueError, TypeError):
return "Multitaper: enter a valid bandwidth (or 'Auto')"
n_tapers = int(np.floor(epoch_sec * bandwidth - 1))
taper_note = f"{n_tapers} taper(s)"
if n_tapers < 1:
taper_note += " — TOO LOW, will fail"
elif n_tapers < 3:
taper_note += " — low, little variance benefit"
return (f"Multitaper bandwidth {bw_label}, smoothing ±{bandwidth / 2:.3f} Hz; "
f"bin spacing {bin_hz:.3f} Hz; ~{taper_note}")
def create_gui():
"""Create the GUI layout for EEG analysis settings (Pipeline Layout)."""
suggested_threads = DEFAULT_THREADS
# Adapt to the available screen height. The log box is the one element that
# can give back rows without hiding a setting, so it shrinks on short
# screens (e.g. 1366x768 laptops) and expands again on tall ones. It also
# has expand_y=True, so it reclaims any spare height at run time.
try:
_screen_h = sg.Window.get_screen_size()[1]
except Exception: # noqa: BLE001 - screen probing must never block the GUI
_screen_h = 1080
log_rows = 8 if _screen_h < 900 else 12
# Seed the relative-power reference range from the 'broadband' band definition
# so the GUI default matches the historical behaviour, while still letting the
# user override it independently of the band definitions.
if "broadband" in FREQUENCY_BANDS:
default_total_min, default_total_max = FREQUENCY_BANDS["broadband"]["range"]
else:
default_total_min, default_total_max = 0.5, 47.0
HEADER_BG = "#2C5784"
HEADER_TEXT = "#FFFFFF"
MAIN_BG = "#F0F2F6"
BUTTON_COLOR = ("#FFFFFF", "#2C5784")
STEP_TEXT_COLOR = "#505050"
sg.theme("Default1")
sg.set_options(font=("Helvetica", 10))
# --- Header ---
header = [
[
sg.Text(
"EEG-Pype Quantitative Analysis",
font=("Helvetica", 20, "bold"),
text_color=HEADER_TEXT,
background_color=HEADER_BG,
pad=(10, 5),
)
],
[
sg.Text(
"Author: Yorben Lodema",
font=("Helvetica", 10, "italic"),
text_color=HEADER_TEXT,
background_color=HEADER_BG,
pad=(10, 5),
)
],
]
# --- Step 1: Setup (Inputs & Epochs) ---
step_1_layout = [
[sg.Text("Step 1: Setup", font=("Helvetica", 12, "bold"), text_color=STEP_TEXT_COLOR, background_color=MAIN_BG)],
[
sg.Frame(
"Input Settings",
[
[sg.Text("Select data folder:", font=("Helvetica", 11, "bold"), background_color=MAIN_BG)],
[sg.Input(key="-FOLDER-", size=(25, 1), enable_events=True),
sg.FolderBrowse(button_color=BUTTON_COLOR)],
[
sg.Text("Folder extension:", background_color=MAIN_BG),
sg.Input(FOLDER_EXTENSION, key="-EXTENSION-", size=(8, 1)),
],
[
sg.Text("Processing threads:", background_color=MAIN_BG),
sg.Input(suggested_threads, key="-THREADS-", size=(5, 1)),
],
[
sg.Checkbox(
"Epoch files have headers", key="-HAS_HEADERS-", default=True, background_color=MAIN_BG
)
],
],
background_color=MAIN_BG, expand_x=True
)
],
[
sg.Frame(
"Epoch Selection",
[
[
sg.Checkbox(
"Limit number of epochs", key="-LIMIT_EPOCHS-", default=False, background_color=MAIN_BG
)
],
[
sg.Text("Number of epochs:", background_color=MAIN_BG),
sg.Input("5", key="-MAX_EPOCHS-", size=(5, 1)),
],
[
sg.Text("Random Seed:", background_color=MAIN_BG),
sg.Input("111", key="-SEED-", size=(5, 1)),
],
],
background_color=MAIN_BG, expand_x=True
)
],
# Step 2 continues here: the complexity measures belong to the same
# analysis step as the spectral and connectivity blocks in the next
# column, and are placed here to keep the window short enough for
# smaller screens.
[sg.Text("Step 2: Complexity", font=("Helvetica", 12, "bold"), text_color=STEP_TEXT_COLOR,
background_color=MAIN_BG, pad=((5, 3), (12, 3)))],
[
sg.Frame(
"Complexity Measures",
[
[sg.Checkbox("Calculate JPE/PE", key="-CALC_JPE-", default=False, background_color=MAIN_BG)],
[sg.Text("Time step (tau):", background_color=MAIN_BG), sg.Input("1", key="-JPE_ST-", size=(5, 1))],
[sg.Checkbox("Invert JPE (1-entropy)", key="-INVERT-", default=True, background_color=MAIN_BG)],
[
sg.Checkbox(
"Calculate Sample Entropy", key="-CALC_SAMPEN-", default=False, background_color=MAIN_BG
)
],
[sg.Text("Order (m):", background_color=MAIN_BG), sg.Input("2", key="-SAMPEN_M-", size=(3, 1))],
[
sg.Checkbox(
"Calculate Approx. Entropy", key="-CALC_APEN-", default=False, background_color=MAIN_BG
)
],
[sg.Text("Order (m):", background_color=MAIN_BG), sg.Input("2", key="-APEN_M-", size=(3, 1))],
[
sg.Text("Tolerance (r):", background_color=MAIN_BG),
sg.Input("0.25", key="-APEN_R-", size=(3, 1)),
],
[
sg.Checkbox(
"Calculate Lempel-Ziv Complexity", key="-CALC_LZC-", default=False, background_color=MAIN_BG
)
],
[
sg.Text("Threshold:", background_color=MAIN_BG),
sg.Combo(
["median", "mean"], default_value="median", key="-LZC_THRESH-", size=(8, 1)
),
],
],
background_color=MAIN_BG, expand_x=True
)
],
]
# --- Step 2: Analysis Configuration ---
step_2_layout = [
[sg.Text("Step 2: Spectral & Connectivity", font=("Helvetica", 12, "bold"), text_color=STEP_TEXT_COLOR, background_color=MAIN_BG)],
[
sg.Frame(
"Spectral Analysis",
[
[sg.Text("Sampling rate (Hz):", background_color=MAIN_BG),
sg.Input(key="-POWER_FS-", size=(8, 1), enable_events=True)],
[sg.Text("Epoch length (s):", background_color=MAIN_BG),
sg.Input(key="-EPOCH_LENGTH-", size=(8, 1), enable_events=True),
sg.Text("(for fr. resolution readout)", background_color=MAIN_BG,
font=("Helvetica", 9, "italic"))],
[
sg.Text("PSD Method:", background_color=MAIN_BG),
sg.Combo(
["Multitaper", "Welch", "FFT"], default_value="Multitaper",
key="-PSD_METHOD-", size=(10, 1), enable_events=True
),
],
[sg.HorizontalSeparator(color="#D0D0D0")],
[
sg.Text("Welch window (ms):", background_color=MAIN_BG, size=(18, 1)),
sg.Input("1000", key="-WELCH_WINDOW-", size=(7, 1), enable_events=True),
sg.Text("overlap (%):", background_color=MAIN_BG),
sg.Input("50", key="-WELCH_OVERLAP-", size=(5, 1)),
],
[
sg.Text("Multitaper bandwidth (Hz):", background_color=MAIN_BG, size=(18, 1)),
sg.Input("Auto", key="-MT_BANDWIDTH-", size=(7, 1), enable_events=True),
],
[
sg.Text("", key="-FREQ_RES-", background_color=MAIN_BG,
font=("Helvetica", 10, "italic"), size=(46, 3)),
],
[sg.HorizontalSeparator(color="#D0D0D0")],
[
sg.Checkbox("Calculate power bands", key="-CALC_POWER-",
default=False, background_color=MAIN_BG),
sg.Text("ref. (Hz):", background_color=MAIN_BG),
sg.Input(f"{default_total_min:g}", key="-TOTAL_MIN-", size=(5, 1)),
sg.Text("-", background_color=MAIN_BG),
sg.Input(f"{default_total_max:g}", key="-TOTAL_MAX-", size=(5, 1)),
],
[
sg.Text(PASSBAND_HINT_DEFAULT, key="-PASSBAND_INFO-", background_color=MAIN_BG,
text_color="#606060", font=("Helvetica", 9, "italic"), size=(48, 1)),
],
[
sg.Checkbox("Calculate peak frequency", key="-CALC_PEAK-",
default=False, background_color=MAIN_BG),
sg.Text("range:", background_color=MAIN_BG),
sg.Input("4", key="-PEAK_MIN-", size=(4, 1)),
sg.Text("-", background_color=MAIN_BG),
sg.Input("13", key="-PEAK_MAX-", size=(4, 1)),
],
[
sg.Checkbox("Calc. spectral variability", key="-CALC_SV-",
default=False, background_color=MAIN_BG),
sg.Text("window (ms):", background_color=MAIN_BG),
sg.Input("2000", key="-SV_WINDOW-", size=(6, 1)),
],
],
background_color=MAIN_BG, expand_x=True
)
],
[
sg.Frame(
"Connectivity",
[
[sg.Checkbox("Calculate PLI", key="-CALC_PLI-", default=False, background_color=MAIN_BG)],
[
sg.Checkbox(
"Calc. PLI MST measures", key="-CALC_PLI_MST-", default=False, background_color=MAIN_BG
)
],
[sg.Checkbox("Calculate PLT", key="-CALC_PLT-", default=False, background_color=MAIN_BG)],
[
sg.Text("PLT Threshold (ms):", background_color=MAIN_BG),
sg.Input("30", key="-PLT_THRESH-", size=(4, 1)),
],
[
sg.Checkbox(
"Calc. PLT MST measures", key="-CALC_PLT_MST-", default=False, background_color=MAIN_BG
)
],
[sg.Checkbox("Calculate AEC", key="-CALC_AEC-", default=False, background_color=MAIN_BG)],
[
sg.Checkbox(
"Use orthogonalization (AECc)", key="-USE_AECC-", default=False, background_color=MAIN_BG
)
],
[
sg.Checkbox(
"Concatenate for AEC(c)",
key="-CONCAT_AECC-",
default=False,
background_color=MAIN_BG,
)
],
[
sg.Checkbox(
"Calc. AEC(c) MST measures",
key="-CALC_AEC_MST-",
default=False,
background_color=MAIN_BG,
)
],
[
sg.Checkbox(
"AEC force positive",
key="-AEC_FORCE_POSITIVE-",
default=True,
background_color=MAIN_BG,
)
],
],
background_color=MAIN_BG, expand_x=True
)
],
]
# --- Step 3: Output & Execution ---
step_3_layout = [
[sg.Text("Step 3: Output & Run", font=("Helvetica", 12, "bold"), text_color=STEP_TEXT_COLOR, background_color=MAIN_BG)],
[
sg.Frame(
"Matrix & Metric Export",
[
[
sg.Checkbox(
"Save conn. matrices (Per Epoch)", key="-SAVE_MATRICES-", default=False, background_color=MAIN_BG
)
],
[
sg.Text("Matrix folder:", background_color=MAIN_BG),
sg.Input("connectivity_matrices", key="-MATRIX_FOLDER-", size=(25, 1)),
],
[sg.Checkbox("Save MST matrices (Per Epoch)", key="-SAVE_MST-", default=False, background_color=MAIN_BG)],
[
sg.Text("MST folder:", background_color=MAIN_BG),
sg.Input("mst_matrices", key="-MST_FOLDER-", size=(25, 1)),
],
[
sg.Checkbox(
"Export Per-Epoch metrics to Excel sheet", key="-SAVE_EPOCH_METRICS-", default=True, background_color=MAIN_BG
)
],
[
sg.Checkbox(
"Save channel-level averages",
key="-SAVE_CHANNEL_AVERAGES-",
default=False,
background_color=MAIN_BG,
)
],
],
background_color=MAIN_BG, expand_x=True
)
],
[
sg.Frame(
"PSD Export",
[
[sg.Checkbox("Save averaged PSD (epochs + ch/region)", key="-SAVE_PSD_AVG-", default=False, background_color=MAIN_BG)],
[sg.Checkbox("Save per-channel/region PSD (epoch-avg)", key="-SAVE_PSD_PERCHANNEL-", default=False, background_color=MAIN_BG)],
[sg.Checkbox("Save PSD plots (PNG, MNE-style)", key="-SAVE_PSD_PLOTS-", default=False, background_color=MAIN_BG)],
[
sg.Text("PSD folder:", background_color=MAIN_BG),
sg.Input("psd_output", key="-PSD_FOLDER-", size=(25, 1)),
],
[sg.Text("(PSD uses broadband epochs + sampling rate)",
background_color=MAIN_BG, font=("Helvetica", 10, "italic"))],
],
background_color=MAIN_BG, expand_x=True
)
],
[
sg.Frame(
"Progress & Execution",
[
[
sg.Button(
"Process",
size=(12, 1),
button_color=BUTTON_COLOR,
font=("Helvetica", 12, "bold"),
bind_return_key=True
),
sg.Button(
"Exit",
size=(8, 1),
button_color=(HEADER_TEXT, "#AB4F4F"),
font=("Helvetica", 12),
)
],
[sg.ProgressBar(100, orientation="h", size=(45, 20), key="-PROGRESS-", bar_color=(HEADER_BG, "#FFFFFF"))],
[sg.Text("Log Output:", background_color=MAIN_BG, font=("Helvetica", 9, "bold"))],
[
sg.Multiline(
size=(60, log_rows),
key="-LOG-",
autoscroll=True,
reroute_stdout=True,
disabled=True,
background_color="#FFFFFF",
text_color="#000000",
font=("Courier New", 8),
expand_x=True, # Allow horizontal expansion
expand_y=True # Allow vertical expansion
)
],
],
background_color=MAIN_BG, expand_x=True, expand_y=True
)
],
]
# --- Main Layout Assembly ---
layout = [
[sg.Column(header, background_color=HEADER_BG, expand_x=True)],
[
sg.Column(step_1_layout, background_color=MAIN_BG, vertical_alignment='top', expand_y=True),
sg.VerticalSeparator(color="#D0D0D0"), # Visual divider
sg.Column(step_2_layout, background_color=MAIN_BG, vertical_alignment='top', expand_y=True),
sg.VerticalSeparator(color="#D0D0D0"), # Visual divider
sg.Column(step_3_layout, background_color=MAIN_BG, vertical_alignment='top', expand_y=True, expand_x=True),
],
]
return sg.Window("EEG-Pype", layout, background_color=MAIN_BG, finalize=True, margins=(0, 0), resizable=True)
def create_matrix_folder_structure(base_folder, matrix_folder_name, mst_folder_name=None):
"""Create folder structure with subject subfolders."""
folders = {
"jpe": os.path.join(base_folder, matrix_folder_name, "jpe"),
"pli": os.path.join(base_folder, matrix_folder_name, "pli"),
"plt": os.path.join(base_folder, matrix_folder_name, "plt"),
"aec": os.path.join(base_folder, matrix_folder_name, "aec"),
}
if mst_folder_name:
folders.update(
{
"pli_mst": os.path.join(base_folder, mst_folder_name, "pli_mst"),
"aec_mst": os.path.join(base_folder, mst_folder_name, "aec_mst"),
"plt_mst": os.path.join(base_folder, mst_folder_name, "plt_mst"),
}
)
# Create base folders
for folder in folders.values():
os.makedirs(folder, exist_ok=True)
return folders
def extract_freq_band(condition):
"""
Parse the filename or condition string to identify frequency band.
Based on the FREQUENCY_BANDS config.
"""
for band_name, band_info in FREQUENCY_BANDS.items():
pattern = band_info["pattern"]
# Add Hz to pattern if not already included
if not pattern.endswith("Hz"):
# Wrap in a non-capturing group: '|' has the lowest precedence, so an
# unwrapped pattern like 'a|b' would attach the optional Hz to 'b' only.
search_pattern = f"(?:{pattern})(\\s*Hz)?"
else:
search_pattern = pattern
if re.search(search_pattern, condition, re.IGNORECASE):
return band_name
return "unknown"
def is_broadband_condition(condition):
"""
Check if condition matches broadband pattern from FREQUENCY_BANDS config.
"""
if "broadband" not in FREQUENCY_BANDS:
return False
pattern = FREQUENCY_BANDS["broadband"]["pattern"]
return bool(re.search(pattern, condition, re.IGNORECASE))
def save_connectivity_matrix(matrix, folder_path, subject, freq_band, feature, channel_names, level_type=None):
"""Save connectivity matrix to CSV with proper channel names, prepending subject to the filename."""
subject_folder = os.path.join(folder_path, subject)
os.makedirs(subject_folder, exist_ok=True)
# Construct the latter part of the filename as before
if level_type:
base_filename_part = f"{level_type}_{freq_band}_{feature}"
else:
base_filename_part = f"{freq_band}_{feature}"
# Prepend the subject identifier
filename = f"{subject}_{base_filename_part}.csv"
# Example: "eeg_filename_id_source_delta_pli.csv" or "eeg_filename_id_delta_pli.csv"
filepath = os.path.join(subject_folder, filename)
# Convert matrix to DataFrame with channel names
df = pd.DataFrame(matrix)
df.index = channel_names
df.columns = channel_names
df.to_csv(filepath, float_format='%.5f')
return filepath
def linear_detrend(data):
"""Apply linear detrending to each channel."""
return signal.detrend(data, axis=0, type="linear")
def read_epoch_array(file_path, has_headers, detrend=True):
"""Read one epoch file into a (samples x channels) array, linearly detrended by default.
Single source of truth for epoch reading, used by both the per-epoch loop and the
concatenated-AEC / spectral-variability paths so they all handle files identically:
- has_headers=True -> first row is treated as channel names (pandas default).
- has_headers=False -> the first row is auto-inspected; a non-numeric row is
dropped as a stray header, and all columns are coerced to numeric.
Returns
-------
tuple (np.ndarray, list[str] | None)
The (detrended) data array, and the channel names parsed from the header
when has_headers=True, else None (caller generates generic names).
"""
if has_headers:
df = pd.read_csv(file_path, sep=None, engine="python")
channel_names = df.columns.tolist()
else:
# Detect and drop a stray non-numeric header row, matching the legacy per-epoch logic.
first_row = pd.read_csv(file_path, sep=None, engine="python", header=None, nrows=1)
is_header = False
try:
first_row.astype(float)
except (ValueError, TypeError):
is_header = True
df = pd.read_csv(
file_path, sep=None, engine="python", header=None, skiprows=1 if is_header else 0
)
for col in df.columns:
df[col] = pd.to_numeric(df[col], errors="coerce")
channel_names = None
arr = df.to_numpy()
if detrend:
arr = linear_detrend(arr)
return arr, channel_names
def calculate_PSD(
data: np.ndarray,
fs: float,
method: str = "multitaper",
freq_range: Optional[tuple[float, float]] = None,
**kwargs,
) -> dict[str, np.ndarray]:
"""
Calculate Power Spectral Density (PSD) using specified method.
Parameters
----------
data : np.ndarray
Time series data (samples x channels)
fs : float
Sampling frequency in Hz
method : str
Method to use for PSD calculation ('multitaper', 'welch', 'fft')
freq_range : tuple, optional
Frequency range to return (min_freq, max_freq)
**kwargs : dict
Method-specific parameters:
Welch:
window_length_ms : float (window length in milliseconds)
overlap_percent : float (0 to 100)
Multitaper:
time_bandwidth : float (default 4)
n_tapers : int (optional, computed from time_bandwidth)
Returns
-------
dict
Dictionary containing:
'frequencies' : np.ndarray
Frequency values
'psd' : np.ndarray
Power spectral density (frequencies x channels)
'spectrogram' : np.ndarray, optional
Time-frequency representation (only if compute_spectrogram=True)
"""
if method not in ["multitaper", "welch", "fft"]:
msg = f"Unknown method: {method}"
raise ValueError(msg)
# Input validation
if not isinstance(data, np.ndarray):
msg = "Data must be a numpy array"
raise TypeError(msg)
if data.ndim != 2: # noqa: PLR2004
msg = "Data must be 2D array (samples x channels)"
raise ValueError(msg)
if fs <= 0:
msg = "Sampling frequency must be positive"
raise ValueError(msg)
# Initialize return dictionary
result = {}
# Calculate PSD based on method
if method == "multitaper":
try:
bandwidth = kwargs.get("bandwidth", None)
frequencies, psd = _calculate_multitaper_psd(data, fs, bandwidth=bandwidth)
result["frequencies"] = frequencies
result["psd"] = psd
except Exception:
logger.exception("Error calculating multitaper PSD")
raise
elif method == "welch":
try:
window_length_ms = kwargs.get("window_length_ms", 1000) # Default 1000ms
overlap_percent = kwargs.get("overlap_percent", 50) # Default 50%
frequencies, psd = _calculate_welch_psd(
data, fs, window_length_ms=window_length_ms, overlap_percent=overlap_percent
)
result["frequencies"] = frequencies
result["psd"] = psd
except Exception:
logger.exception("Error calculating Welch PSD")
raise
elif method == "fft":
try:
frequencies, psd = _calculate_fft_psd(data, fs)
result["frequencies"] = frequencies
result["psd"] = psd
except Exception:
logger.exception("Error calculating FFT PSD")
raise
# Apply frequency range if specified
if freq_range is not None:
fmin, fmax = freq_range
if not (0 <= fmin < fmax <= fs / 2):
msg = f"Invalid frequency range: {freq_range}"
raise ValueError(msg)
freq_mask = (result["frequencies"] >= fmin) & (result["frequencies"] <= fmax)
result["frequencies"] = result["frequencies"][freq_mask]
result["psd"] = result["psd"][freq_mask]
return result
def _calculate_welch_psd(
data: np.ndarray, fs: float, window_length_ms: float = 1000, overlap_percent: float = 50
) -> tuple[np.ndarray, np.ndarray]:
"""
Calculate PSD using Welch's method.
Parameters
----------
data : np.ndarray
Time series data (samples x channels)
fs : float
Sampling frequency in Hz
window_length_ms : float
Length of each segment in milliseconds
overlap_percent : float
Overlap between segments in percentage (0-100)
Returns
-------
frequencies : np.ndarray
Frequency values
psd : np.ndarray
Power spectral density (frequencies x channels)
"""
# Convert window length from ms to samples
nperseg = int((window_length_ms / 1000) * fs)
# Convert overlap from percentage to samples
noverlap = int(nperseg * (overlap_percent / 100))
# Initialize array for PSD results
n_channels = data.shape[1]
# Calculate PSD for first channel to get frequency axis
frequencies, temp_psd = signal.welch(
data[:, 0], fs=fs, nperseg=nperseg, noverlap=noverlap, detrend=False, scaling="density"
)
# Initialize PSD array with correct dimensions
psd = np.zeros((len(frequencies), n_channels))
psd[:, 0] = temp_psd
# Calculate for remaining channels
for ch in range(1, n_channels):
_, psd[:, ch] = signal.welch(
data[:, ch],
fs=fs,
nperseg=nperseg,
noverlap=noverlap,
detrend=False,
scaling="density",
)
return frequencies, psd
def _calculate_fft_psd(data: np.ndarray, fs: float) -> tuple[np.ndarray, np.ndarray]:
"""Single-segment PSD via a Hann-windowed FFT over the whole epoch.
Uses the same one-sided density scaling as scipy.signal.welch so absolute
power is consistent across FFT, Welch and multitaper: every bin except DC
(and Nyquist, when present) is doubled to account for the discarded
negative-frequency half.
"""
n_samples = data.shape[0]
n_channels = data.shape[1]
window = signal.windows.hann(n_samples)
# scipy-style density normalisation: 1 / (fs * sum(window**2))
scale = 1.0 / (fs * np.sum(window ** 2))
frequencies = np.fft.rfftfreq(n_samples, d=1 / fs)
psd = np.zeros((len(frequencies), n_channels))
for ch in range(n_channels):
fft_data = np.fft.rfft(data[:, ch] * window)
ch_psd = (np.abs(fft_data) ** 2) * scale
# One-sided correction: double all but DC and (if present) Nyquist.
ch_psd[1:] *= 2
if n_samples % 2 == 0:
ch_psd[-1] /= 2 # Nyquist bin is unique; undo its doubling
psd[:, ch] = ch_psd
return frequencies, psd
def _calculate_multitaper_psd(data: np.ndarray, fs: float, bandwidth=None):
"""Calculate PSD using MNE's multitaper implementation.
Parameters
----------
data : np.ndarray
Time series data (samples x channels).
fs : float
Sampling frequency in Hz.
bandwidth : float or None
Full multitaper resolution bandwidth in Hz (frequencies within
± bandwidth/2 are smoothed together). None uses MNE's default of
8 * fs / n_samples.
"""
fmax = min(60.0, fs / 2.0)
psds, freqs = mne.time_frequency.psd_array_multitaper(
data.T,
sfreq=fs,
fmin=0,
fmax=fmax,
bandwidth=bandwidth,
normalization="full",
n_jobs=1,
verbose=False,
)
return freqs, psds.T
def save_psd_outputs(
mean_psd,
frequencies,
channel_names,
subject,
condition,
freq_band,
level_type,
psd_folder,
save_avg=False,
save_perchannel=False,
save_plots=False,
):
"""Write the epoch-averaged PSD to CSV and/or PNG (broadband conditions only).
Parameters
----------
mean_psd : np.ndarray
Epoch-averaged PSD, shape (n_freqs x n_channels).
frequencies : np.ndarray
Frequency axis (length n_freqs).
channel_names : list[str] | None
Channel / region labels; generic names are generated if missing.
psd_folder : str
Base output folder. A per-subject subfolder is created inside it.
save_avg : bool
Write the channel-averaged PSD (single curve) to CSV.
save_perchannel : bool
Write the per-channel PSD (frequencies x channels) to CSV.
save_plots : bool
Also write MNE-style PNG plots (power in dB vs frequency).
"""
if mean_psd is None or frequencies is None:
return
n_channels = mean_psd.shape[1]
if not channel_names or len(channel_names) != n_channels:
channel_names = [f"Channel_{i + 1}" for i in range(n_channels)]
global_psd = np.nanmean(mean_psd, axis=1) # channel-averaged curve
subject_folder = os.path.join(psd_folder, subject)
os.makedirs(subject_folder, exist_ok=True)
stem = f"{subject}_{level_type}_{freq_band}_psd"
if save_avg:
pd.DataFrame({"frequency_Hz": frequencies, "psd": global_psd}).to_csv(
os.path.join(subject_folder, f"{stem}_global_avg.csv"),
index=False, float_format="%.6e",
)
if save_perchannel:
df_pc = pd.DataFrame(mean_psd, index=frequencies, columns=channel_names)
df_pc.index.name = "frequency_Hz"
df_pc.to_csv(
os.path.join(subject_folder, f"{stem}_per_channel.csv"),
float_format="%.6e",
)
if save_plots:
try:
_plot_psd(frequencies, mean_psd, global_psd, subject_folder, stem,
title=f"{subject} — {condition}")
except Exception:
logger.exception(f"Error plotting PSD for {subject} - {condition}")
def _plot_psd(frequencies, mean_psd, global_psd, out_folder, stem, title=""):
"""Render two MNE-style PSD figures (power in dB vs frequency) as PNGs."""
import matplotlib
matplotlib.use("Agg") # headless / multiprocessing-safe backend
import matplotlib.pyplot as plt
eps = 1e-20 # guard against log10(0)
# Per-channel: faint per-channel traces + bold channel mean.
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(frequencies, 10 * np.log10(mean_psd + eps),
color="0.7", linewidth=0.5, alpha=0.6)
ax.plot(frequencies, 10 * np.log10(global_psd + eps),
color="C0", linewidth=2.0, label="Mean")
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Power (dB)")
ax.set_title(f"{title} — per channel")
ax.legend(loc="upper right")
fig.tight_layout()
fig.savefig(os.path.join(out_folder, f"{stem}_per_channel.png"), dpi=300)
plt.close(fig)
# Channel-averaged: single mean curve.
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(frequencies, 10 * np.log10(global_psd + eps),
color="C0", linewidth=2.0)
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Power (dB)")
ax.set_title(f"{title} — channel-averaged")
fig.tight_layout()
fig.savefig(os.path.join(out_folder, f"{stem}_global_avg.png"), dpi=300)
plt.close(fig)
def calculate_sampen_for_channels(data, m=2):
"""
Calculate Sample Entropy for each channel using antropy.
This version includes robust type casting to match Numba's requirements.
"""
n_channels = data.shape[1]
sampen_values = np.zeros(n_channels)
# 1. Force the order parameter to be a 32-BIT INTEGER (i4).
# This is required by the explicit Numba signature in antropy.
order_m = np.int32(m)
for ch in range(n_channels):
try:
# 2. Force the data to be a C-contiguous array of 64-bit floats (f8[:]).
# Numba can be sensitive to array memory layout.
channel_data = np.ascontiguousarray(data[:, ch], dtype=np.float64)
if np.std(channel_data) == 0:
sampen_values[ch] = 0
continue
sampen_values[ch] = sample_entropy(channel_data, order=order_m)
if ch % 10 == 0:
logger.info(f"Processed SampEn for {ch}/{n_channels} channels")
except Exception:
logger.exception(f"Error calculating SampEn for channel {ch}")
sampen_values[ch] = np.nan
return sampen_values
def calculate_lzc_for_channels(data, threshold="median"):
"""Calculate Lempel-Ziv Complexity for each channel using antropy.
Parameters
----------
data : np.ndarray
Time series data (samples x channels).
threshold : str or float
Binarization strategy. ``"median"`` (default) uses the per-channel
median; ``"mean"`` uses the per-channel mean; a float value is used
directly as a fixed threshold.
Returns
-------
np.ndarray
Normalized LZC value for each channel.
"""
n_channels = data.shape[1]
lzc_values = np.zeros(n_channels)
for ch in range(n_channels):
try:
x = np.ascontiguousarray(data[:, ch], dtype=np.float64)
if len(x) < 2 or np.std(x) == 0:
lzc_values[ch] = 0.0
continue
# Binarize
if threshold == "median":
thr = np.median(x)
elif threshold == "mean":
thr = np.mean(x)
else:
thr = float(threshold)
binary = (x >= thr).astype(np.int32)
# antropy handles LZ76 + normalization internally
lzc_values[ch] = lziv_complexity(binary, normalize=True)
if ch % 10 == 0:
logger.info(f"Processed LZC for {ch}/{n_channels} channels")
except Exception:
logger.exception(f"Error calculating LZC for channel {ch}")
lzc_values[ch] = np.nan
return lzc_values
def calculate_apen_for_channels(data, m=2, r=0.25):
"""Calculate Approximate Entropy for each channel.
Follows Pincus 1995, with optimized implementation using vectorization.
"""
n_channels = data.shape[1]
apen_values = np.zeros(n_channels)
for ch in range(n_channels):
try:
# Get channel data
x = data[:, ch]
# Scale r by standard deviation of the data
r_scaled = r * np.std(x)
# Calculate phi(m) and phi(m+1)
phi_m = _phi_vectorized(x, m, r_scaled)
phi_m_plus_1 = _phi_vectorized(x, m + 1, r_scaled)
# Calculate ApEn
apen_values[ch] = phi_m - phi_m_plus_1
if ch % 10 == 0: # Log progress every 10 channels
logger.info(f"Processed ApEn for {ch}/{n_channels} channels")
except Exception:
logger.exception(f"Error calculating ApEn for channel {ch}")
apen_values[ch] = np.nan
return apen_values
def _phi_vectorized(x, m, r):
"""
Vectorized calculation of Φᵐ(r) following Pincus 1995.
"""
N = len(x)
N_m = N - m + 1
# Create embedding matrix efficiently
# Each row is a pattern of length m
patterns = np.zeros((N_m, m))
for i in range(m):
patterns[:, i] = x[i : i + N_m]
# Calculate distances using broadcasting
# This computes the maximum absolute difference between all pairs of patterns
diff = np.abs(patterns[:, None, :] - patterns[None, :, :])
max_diff = np.max(diff, axis=2)
# Count similar patterns (within tolerance r)
similar_patterns = np.sum(max_diff <= r, axis=1)
# Normalize counts by N_m
C = similar_patterns / N_m
# Calculate Φᵐ(r) with small constant to avoid log(0)
return np.mean(np.log(C + 1e-10))
def calculate_spectral_variability(data_values, fs, window_length=2000, total_power_range=None):
"""Calculate spectral variability per channel from concatenated broadband data.
Uses FREQUENCY_BANDS for band definitions.
- Expects pre-concatenated data with channel means already removed.
- ``total_power_range`` sets the denominator for the relative-power time
series whose coefficient of variation is reported. It defaults to
FREQUENCY_BANDS["broadband"]["range"], and should normally be the same
range passed to calculate_power_bands so the two outputs are comparable.
"""
try:
num_samples, num_channels = data_values.shape
samples_per_window = int(window_length * fs / 1000)
# Require at least 3 windows for a meaningful coefficient of variation
if num_samples < 3 * samples_per_window:
logger.warning(
f"Data length ({num_samples}) too short for meaningful "
f"variability calculation with window length {samples_per_window} samples."
)
return None
# 1) Identify the reference range used for total power
if total_power_range is None:
if "broadband" not in FREQUENCY_BANDS:
logger.error(
"No reference range available for spectral variability "
"('broadband' missing from FREQUENCY_BANDS and no explicit range given)"
)
return None
total_power_range = FREQUENCY_BANDS["broadband"]["range"]
broadband_min, broadband_max = total_power_range
# Prepare output dict of CV values
cv_values = {}
for band_name in FREQUENCY_BANDS:
if band_name.lower() == "broadband":
continue
cv_values[band_name] = np.zeros(num_channels)
# 2) Loop over channels and calculate spectrogram
for channel in range(num_channels):
try:
# Compute spectrogram for this channel
f, t, Sxx = signal.spectrogram(
data_values[:, channel],
fs=fs,
nperseg=samples_per_window,
noverlap=samples_per_window // 2,
detrend="constant",
window="hann",
)
# Create mask for broadband total power
total_mask = (f >= broadband_min) & (f < broadband_max)
if not np.any(total_mask):
logger.error(
f"No spectrogram bins in the reference range "
f"{broadband_min}-{broadband_max} Hz."
)
for band_name in cv_values: # noqa: PLC0206
cv_values[band_name][channel] = np.nan
continue
total_power = np.sum(Sxx[total_mask, :], axis=0) # shape: (time_windows,)
# 3) Loop over the user-defined frequency bands
for band_name, band_info in FREQUENCY_BANDS.items():
if band_name.lower() == "broadband":
continue # skip calculating a separate "broadband" measure
low_freq, high_freq = band_info["range"]
band_mask = (f >= low_freq) & (f < high_freq)
if not np.any(band_mask):
# If no frequencies found in this range, skip
cv_values[band_name][channel] = np.nan
continue
band_power = np.sum(Sxx[band_mask, :], axis=0) # shape: (time_windows,)
# Compute relative power time series
with np.errstate(divide="ignore", invalid="ignore"):
relative_power = np.where(total_power > 0, band_power / total_power, 0)
# Remove NaN / Inf
valid_power = relative_power[np.isfinite(relative_power)]
if len(valid_power) > 0:
# Coefficient of Variation: std / mean
cv_values[band_name][channel] = np.std(valid_power) / np.mean(valid_power)
else:
cv_values[band_name][channel] = np.nan
except Exception:
logger.exception(f"Error processing channel {channel}")
# Fill with NaN for all bands on this channel
for band_name in cv_values: # noqa: PLC0206
cv_values[band_name][channel] = np.nan
return cv_values
except Exception:
logger.exception("Error in spectral variability calculation")
return None
def smooth_spectrum_savgol(power_spectrum: np.ndarray, window_length: int = 5, polyorder: int = 2) -> np.ndarray:
"""
Apply Savitzky-Golay smoothing to a power spectrum.
"""
if window_length % 2 == 0:
raise ValueError("window_length must be an odd integer.")
if polyorder >= window_length:
raise ValueError("polyorder must be less than window_length.")
if len(power_spectrum) < window_length:
# Not enough data points to apply the filter with the given window length.
return power_spectrum
smoothed_spectrum = signal.savgol_filter(power_spectrum, window_length, polyorder)
return smoothed_spectrum
def calculate_avg_peak_frequency(frequencies, psd, freq_range=(4, 13), sg_window_length=5, sg_polyorder=2):
"""
Calculate peak frequency using pre-computed PSD with improved peak detection.
"""
num_channels = psd.shape[1]
peak_frequencies = np.zeros(num_channels)
# Create frequency mask
freq_mask = (frequencies >= freq_range[0]) & (frequencies <= freq_range[1])
freq_range_idx = np.where(freq_mask)[0]
if len(freq_range_idx) == 0:
logger.warning(f"No frequencies found in range {freq_range[0]}-{freq_range[1]} Hz")
return np.full(num_channels, np.nan)
# Get masked frequencies and PSD
frequencies_masked = frequencies[freq_mask]
psd_masked = psd[freq_mask, :]
for channel in range(num_channels):
try:
channel_psd = psd_masked[:, channel]
# Apply Savitzky-Golay smoothing
smoothed_psd = smooth_spectrum_savgol(channel_psd,
window_length=sg_window_length,
polyorder=sg_polyorder)
# Find all peaks (using the basic signal.find_peaks for initial candidates)
peak_indices = signal.find_peaks(smoothed_psd)[0]
if len(peak_indices) == 0:
peak_frequencies[channel] = np.nan
continue
# Calculate peak properties (prominences) for these peaks
# Ensure smoothed_psd is not empty and peak_indices are valid
if smoothed_psd.size > 0 and peak_indices.size > 0 and np.all(peak_indices < len(smoothed_psd)):
peak_props = signal.peak_prominences(smoothed_psd, peak_indices)
prominences = peak_props[0]
else:
peak_frequencies[channel] = np.nan # Not enough data or invalid indices for prominence
continue
# Sort peaks by prominence
if prominences.size > 0:
sorted_peak_indices = peak_indices[np.argsort(-prominences)]
else:
peak_frequencies[channel] = np.nan # No prominences to sort by
continue
if len(sorted_peak_indices) > 0:
# Get the frequency of the most prominent peak
peak_frequencies[channel] = frequencies_masked[sorted_peak_indices[0]]
else:
peak_frequencies[channel] = np.nan
except Exception:
logger.exception(f"Error calculating peak frequency for channel {channel}")
peak_frequencies[channel] = np.nan
return peak_frequencies
def integrate_psd_band(frequencies, psd, frequency_mask):
"""Integrate PSD over a selected frequency range.
Parameters
----------
frequencies : np.ndarray
One-dimensional frequency axis in Hz.
psd : np.ndarray
PSD values with shape frequencies x channels. Values are expected to
be power spectral density, for example µV²/Hz or V²/Hz.
frequency_mask : np.ndarray
Boolean mask selecting the frequencies to integrate.
Returns
-------
np.ndarray
Integrated power for each channel.
"""
selected_frequencies = frequencies[frequency_mask]
selected_psd = psd[frequency_mask, :]
if selected_frequencies.size == 0:
raise ValueError("No frequency bins available for PSD integration.")
if selected_frequencies.size == 1:
# Trapezoidal integration over a single point would return zero.
# Approximate the integral using the frequency-bin width.
if frequencies.size < 2:
raise ValueError(
"At least two PSD frequency bins are required to estimate "
"integrated power."
)
frequency_resolution = float(np.median(np.diff(frequencies)))
return selected_psd[0, :] * frequency_resolution
return trapezoid(
selected_psd,
x=selected_frequencies,
axis=0,
)
def calculate_power_bands(frequencies, psd, total_power_range=None):
"""Calculate integrated absolute and relative power.
The input PSD is numerically integrated over each frequency band. Absolute
power therefore has units of signal squared rather than signal squared per Hz,
and does not scale directly with the PSD frequency resolution.
The denominator used for relative power ("total power") is an explicit
parameter rather than being implicitly tied to the 'broadband' entry, so the
reference range can be chosen independently of the band definitions.
Parameters
----------
frequencies : np.ndarray
1D array of frequency values.
psd : np.ndarray
2D array of power spectral density (frequencies x channels).
total_power_range : tuple (fmin, fmax), optional
Frequency range used as the denominator for relative power. Defaults to
FREQUENCY_BANDS["broadband"]["range"] for backwards compatibility.
Returns
-------
tuple
(powers, channel_powers)
- powers: dict with mean abs/rel power across channels per band
- channel_powers: dict with channel-level arrays (one entry per band)
"""
if total_power_range is None:
total_power_range = FREQUENCY_BANDS["broadband"]["range"]
total_min, total_max = total_power_range
total_mask = (frequencies >= total_min) & (frequencies < total_max)
if not np.any(total_mask):
msg = f"No PSD bins in reference range {total_min}-{total_max} Hz"
raise ValueError(msg)
# Warn if the PSD does not actually reach the requested top edge (e.g. the
# multitaper path truncates at 60 Hz, or fs/2 is lower than requested).
if frequencies.max() < total_max:
logger.warning(
f"Reference range top edge {total_max} Hz exceeds highest PSD bin "
f"({frequencies.max():.2f} Hz); total power is effectively "
f"{total_min}-{frequencies.max():.2f} Hz."
)
# Integrate PSD over the relative-power reference range
total_power = integrate_psd_band(
frequencies,
psd,
total_mask,
)
# Prepare output containers
powers = {}
channel_powers = {}
# Loop over the frequency bands
for band_name, band_info in FREQUENCY_BANDS.items():
# Calculate for all bands including 'broadband'
fmin, fmax = band_info["range"]
band_mask = (frequencies >= fmin) & (frequencies < fmax)
# Integrate PSD over this frequency band
if not np.any(band_mask):
logger.warning(
f"No PSD bins available for frequency band "
f"{band_name} ({fmin}-{fmax} Hz)."
)
abs_power = np.full(psd.shape[1], np.nan)
else:
abs_power = integrate_psd_band(
frequencies,
psd,
band_mask,
)
# Compute relative power. Note that the 'broadband' band is no longer
# forced to 1.0: with a user-chosen reference range its relative power is
# a useful diagnostic (how much of the reference range it covers).
with np.errstate(divide="ignore", invalid="ignore"):
rel_power = np.where(total_power > 0, abs_power / total_power, np.nan)
# Store mean abs/rel power across channels
powers[f"{band_name}_abs_power"] = np.nanmean(abs_power)
powers[f"{band_name}_rel_power"] = np.nanmean(rel_power)
# Also store channel-level arrays
channel_powers[f"{band_name}_abs_power"] = np.nan_to_num(abs_power, nan=np.nan)
channel_powers[f"{band_name}_rel_power"] = np.nan_to_num(rel_power, nan=np.nan)
return powers, channel_powers
def calculate_mst_measures(connectivity_matrix, used_channels=None):
"""
Calculate MST measures from a connectivity matrix with additional error handling for disconnected graphs.
Args:
connectivity_matrix (numpy.ndarray): Square connectivity matrix (e.g., PLI matrix)
used_channels (numpy.ndarray, optional): Boolean array indicating which channels are used.
If None, all channels are considered used.
Returns
-------
tuple: (dict of MST measures, MST matrix, bool indicating success)
"""
# Initialize used_channels if not provided
if used_channels is None:
used_channels = np.ones(len(connectivity_matrix), dtype=bool)
# Get number of total channels (N) and used channels (M)
n_total = len(connectivity_matrix) # N in BrainWave
n_used = np.sum(used_channels) # M in BrainWave
norm_factor = n_used - 1 # (M-1) for initial normalization
# Create MST from connectivity matrix
# Using -connectivity to get maximum spanning tree
mst_matrix = minimum_spanning_tree(-connectivity_matrix).toarray()
# Convert to NetworkX graph to check connectivity
G = nx.from_numpy_array(-mst_matrix)
# Check if the graph is connected
if not nx.is_connected(G):
return None, None, False
mst_matrix = -mst_matrix
G = nx.from_numpy_array(mst_matrix)
measures = {}
try:
# 1. Maximum degree calculation
degrees = defaultdict(float)
for edge in G.edges():
degrees[edge[0]] += 1.0 / norm_factor # norm_factor is (M-1)
degrees[edge[1]] += 1.0 / norm_factor
measures["degree"] = max(degrees.values()) if degrees else 0
# 2. Eccentricity - normalize by (M-1)
eccentricity = nx.eccentricity(G)
normalized_eccentricity = {node: ecc / norm_factor for node, ecc in eccentricity.items() if used_channels[node]}
measures["eccentr"] = np.mean(list(normalized_eccentricity.values()))
# 3. Betweenness centrality
betweenness = nx.betweenness_centrality(G)
measures["betweenness"] = max(list(betweenness.values()))
# 4. Diameter - normalize by (M-1)
raw_diameter = nx.diameter(G)
measures["diameter"] = raw_diameter / norm_factor
# 5. Leaf fraction
leaf_nodes = sum(1 for node, deg in degrees.items() if abs(deg - 1.0 / norm_factor) < 1e-10) # noqa: PLR2004
measures["leaf"] = leaf_nodes / n_used
max_betweenness = max(betweenness.values()) if betweenness else 0
if max_betweenness > 0:
measures["hierarchy"] = leaf_nodes / (2 * max_betweenness * norm_factor)
else:
measures["hierarchy"] = 0
# 6. Kappa (degree divergence)
sum_x = sum((n_total - 1) * deg for node, deg in degrees.items() if used_channels[node])
sum_x2 = sum(((n_total - 1) * deg) ** 2 for node, deg in degrees.items() if used_channels[node])
measures["kappa"] = sum_x2 / sum_x if sum_x > 0 else 0
# 7. Average shortest path (ASP)
paths = dict(nx.all_pairs_shortest_path_length(G))
sum_distances = 0
for i in range(n_total):
if used_channels[i]:
node_distances = 0
for j in range(n_total):
if used_channels[j] and i != j and i in paths and j in paths[i]:
node_distances += paths[i][j]
sum_distances += node_distances
measures["asp"] = sum_distances / (n_used * (n_used - 1)) if n_used > 1 else 0
# 8. Tree efficiency (Teff)
normalized_diam = raw_diameter / norm_factor
measures["teff"] = 1.0 - (normalized_diam * (n_used - 1)) / (n_used - (n_used - 1) * measures["leaf"] + 1.0)
# 9. R (degree correlation)
degree_pairs = []
for edge in G.edges():
i, j = edge[0], edge[1]
if used_channels[i] and used_channels[j]:
degree_pairs.append((degrees[i], degrees[j]))
if degree_pairs:
deg_i, deg_j = zip(*degree_pairs)
deg_i = np.array(deg_i)
deg_j = np.array(deg_j)
mean_i = np.mean(deg_i)
mean_j = np.mean(deg_j)
cov = np.mean((deg_i - mean_i) * (deg_j - mean_j))
var_i = np.mean((deg_i - mean_i) ** 2)
var_j = np.mean((deg_j - mean_j) ** 2)
if var_i * var_j > 0:
measures["r"] = cov / np.sqrt(var_i * var_j)
else:
measures["r"] = 0
else:
measures["r"] = 0
# 10. Mean edge weight
edge_weights = [abs(d.get("weight", 1.0)) for _, _, d in G.edges(data=True)]
measures["mean"] = np.mean(edge_weights) if edge_weights else 0
# 11. Reference value
mst_sum = np.sum(abs(mst_matrix[used_channels][:, used_channels]))
orig_sum = np.sum(connectivity_matrix[used_channels][:, used_channels])
measures["ref"] = mst_sum / orig_sum if orig_sum > 0 else 0
return measures, mst_matrix, True
except Exception:
logger.exception("Error in MST measures calculation")
return None, None, False
def calculate_pli(data):
"""Optimized PLI calculation using vectorization."""
analytic_signal = hilbert(data, axis=0)
phases = np.angle(analytic_signal)
n_channels = data.shape[1]
pli = np.zeros((n_channels, n_channels))
# Vectorized phase difference calculation
for i in range(n_channels):
phase_diffs = phases[:, i : i + 1] - phases[:, i:]
signs = np.sign(np.sin(phase_diffs))
means = np.abs(np.mean(signs, axis=0))
pli[i, i:] = means
pli[i:, i] = means
return pli
def PLT(data, fs, threshold_ms=30):
"""
Calculates the Phase Lag Time (PLT) index.
Expects data in (samples, channels) or (channels, samples).
"""
# 1. Data Standardisation
if hasattr(data, 'values'): # Check if pandas DataFrame
data = data.values
# Ensure shape is (channels x samples) for Hilbert
if data.shape[0] > data.shape[1]:
data = data.T
n_channels = data.shape[0]
n_samples = data.shape[1]
# 2. Calculate Threshold in Samples
min_samples = int((threshold_ms / 1000) * fs)
# 3. Analytic Signal
analytic = hilbert(data, axis=1)
PLT_matrix = np.zeros((n_channels, n_channels))
total_time_seconds = n_samples / fs
for i in range(n_channels):
for j in range(i + 1, n_channels):
# Phase Difference & Sign Logic
phase_diff_complex = analytic[i] / analytic[j]
sign_diff = np.sign(phase_diff_complex.imag)
sign_diff[sign_diff == 0] = 1
# Detect Crossings
crossings = np.where(np.abs(np.diff(sign_diff)) > 1)[0]
# 4. Apply Time-Based Threshold
if crossings.size > 1 and min_samples > 0:
valid_crossings = [crossings[0]]
for k in range(1, len(crossings)):
if crossings[k] - valid_crossings[-1] > min_samples:
valid_crossings.append(crossings[k])
crossings = np.array(valid_crossings)
num_crossings = crossings.size
# 5. Calculate Score
if num_crossings == 0:
PLT_matrix[i, j] = 1.0
else:
avg_t = total_time_seconds / (num_crossings + 1)
PLT_matrix[i, j] = 1 - np.exp(-avg_t)
# Mirror matrix
PLT_matrix = PLT_matrix + PLT_matrix.T
return PLT_matrix
def calculate_aecc(data, orthogonalize=False, force_positive=True):
"""
Calculate amplitude envelope correlation with optional orthogonalization.
Parameters
----------
data : numpy array (time points * channels)
EEG data array
orthogonalize : bool, optional
Whether to perform orthogonalization
force_positive : bool, optional
Whether to force negative correlations to zero
Returns
-------
numpy array (channels * channels)
AEC(c) correlation matrix
"""
def process_correlation(corr):
"""Process correlation based on force_positive setting."""
return max(0.0, corr) if force_positive else corr
n_channels = data.shape[1]
correlation_matrix = np.zeros((n_channels, n_channels))
if orthogonalize:
# Process all channels pairwise
for i in range(n_channels):
for j in range(i + 1, n_channels):
# Orthogonalize in both directions
d_orth_ij = data[:, j] - np.dot(data[:, j], data[:, i]) * data[:, i] / np.dot(data[:, i], data[:, i])
d_orth_ji = data[:, i] - np.dot(data[:, i], data[:, j]) * data[:, j] / np.dot(data[:, j], data[:, j])
# Calculate envelopes
env_i = np.abs(hilbert(data[:, i]))
env_j = np.abs(hilbert(data[:, j]))
env_orth_ij = np.abs(hilbert(d_orth_ij))
env_orth_ji = np.abs(hilbert(d_orth_ji))
# Calculate correlations
corr_ij = process_correlation(np.corrcoef(env_i, env_orth_ij)[0, 1])
corr_ji = process_correlation(np.corrcoef(env_j, env_orth_ji)[0, 1])
# Update correlation matrix
correlation_matrix[i, j] = (corr_ij + corr_ji) / 2
correlation_matrix[j, i] = correlation_matrix[i, j]
# Clean up memory
del env_i, env_j, env_orth_ij, env_orth_ji, d_orth_ij, d_orth_ji
else:
# Simple amplitude envelope correlation without orthogonalization
envs = np.abs(hilbert(data, axis=0))
correlation_matrix = np.corrcoef(envs.T)
if force_positive:
correlation_matrix = np.maximum(correlation_matrix, 0)
del envs
# Zero the diagonal
np.fill_diagonal(correlation_matrix, 0)
return correlation_matrix
def calculate_pe(data, n=4, st=1):
"""Calculate Permutation Entropy for each channel.
Parameters
----------
data : numpy array (time points * channels)
n : int, embedding dimension
st : int, time delay (should scale with sampling frequency)
Returns
-------
numpy array : PE values for each channel
"""
sz = data.shape[0]
combinations = list(itertools.permutations(np.arange(0, n), n))
PEs = []
for ch in range(data.shape[1]):
pattern_counts = np.zeros(len(combinations))
# Step size for moving between patterns should be fixed (e.g., 1)
# Only the sampling interval (st) within patterns should scale with frequency
for i in range(0, sz - n * st, 1):
dat_array = data[i : i + n * st : st, ch]
if len(dat_array) < n:
break
dat_order = dat_array.argsort()
rank = dat_order.argsort()
pattern_idx = combinations.index(tuple(rank))
pattern_counts[pattern_idx] += 1
# Calculate PE
total_patterns = np.sum(pattern_counts)
if total_patterns > 0:
prob = pattern_counts[pattern_counts > 0] / total_patterns
entr = -np.sum(prob * np.log(prob))
pe_norm = entr / np.log(math.factorial(n))
PEs.append(pe_norm)
else:
PEs.append(np.nan)
return np.array(PEs)
def find_mirror_patterns(combinations):
"""Create a lookup dictionary for mirror patterns (assumes 0-based ranks)."""
if not combinations:
return {}
mirrors = {}
n = len(combinations[0]) # Determine embedding dimension from the permutation length
mirror_sum = n - 1
for i, perm1 in enumerate(combinations):
# Optimization: only need to check j > i
for j in range(i + 1, len(combinations)):
perm2 = combinations[j]
# Check if perm2 is the mirror of perm1
is_mirror = True
for k in range(n): # Iterate through elements of the permutations
if perm1[k] + perm2[k] != mirror_sum:
is_mirror = False
break
if is_mirror:
mirrors[i] = j
mirrors[j] = i
return mirrors
def is_volume_conduction(pattern1, pattern2, mirrors):
"""Check for volume conduction."""
return pattern1 == pattern2 or pattern2 == mirrors.get(pattern1, -1)
def calculate_jpe(data, n=4, st=1, invert=True):
"""Calculate joint permutation entropy with corrected time delay handling.
Parameters
----------
data : numpy array (time points * channels)
n : int, embedding dimension
st : int, time delay (should scale with sampling frequency)
invert : bool, whether to return 1-JPE
"""
data = np.asarray(data)
sz = data.shape[0]
combinations = list(itertools.permutations(np.arange(0,n), n))
mirrors = find_mirror_patterns(combinations)
rank_inds = []
for i in range(0, sz - n * st, 1):
dat_array = data[i : i + n * st : st, :]
if dat_array.shape[0] < n: # Safety check
break
dat_order = dat_array.argsort(axis=0)
rank = dat_order.argsort(axis=0)
rank_inds.append([combinations.index(tuple(r)) for r in rank.T])
rank_inds = np.array(rank_inds).T
JPE = np.zeros((data.shape[1], data.shape[1]))
for ch, x in enumerate(rank_inds):
for ind, y in enumerate(rank_inds):
if ind > ch:
pattern_counts = 0
jpe_mat = np.zeros((len(combinations), len(combinations)))
for i, j in zip(x, y):
if not is_volume_conduction(i, j, mirrors):
pattern_counts += 1
jpe_mat[i, j] += 1
if pattern_counts > 0:
jpe_mat = jpe_mat / pattern_counts
prob = jpe_mat[jpe_mat > 0]
entr = -np.sum(prob * np.log(prob))
jpe_norm = entr / np.log(math.factorial(n) * math.factorial(n) - 2 * math.factorial(n))
JPE[ch, ind] = 1 - jpe_norm if invert else jpe_norm
return JPE + JPE.T
def extract_epoch_id(filename):
"""Extract epoch identifier from filename for logging and saving."""
match = re.search(r'Epoch_?(\d+)', os.path.basename(filename), re.IGNORECASE)
return f"Epoch_{match.group(1)}" if match else "Epoch_?"
def parse_epoch_filename(filename):
"""Parse epoch filename to extract components.
Handles both legacy format:
subject_Sensor_level_4.0-8.0 Hz_Epoch_1.txt
subject_Source_level_4.0-8.0 Hz_Epoch_1.txt
And new multi-atlas format:
subject_Source_desikan_cortical_4.0-8.0 Hz_Epoch_1.txt
subject_Source_bna_full_4.0-8.0 Hz_Epoch_1.txt
"""
base_name = filename.split("_")[0]
# Try new method+atlas source format first: _Source_{method}_{atlas}_{region}_
method_atlas_match = re.search(
r"_(Source)_(beamformer|sLORETA|eLORETA|dSPM|MNE)_(desikan|bna|aal2|aal3)_(cortical|full)_",
filename,
re.IGNORECASE,
)
# Also try atlas-only format (no method): _Source_{atlas}_{region}_
atlas_match = re.search(
r"_(Source)_(desikan|bna|aal2|aal3)_(cortical|full)_",
filename,
re.IGNORECASE,
)
if method_atlas_match:
method = method_atlas_match.group(2)
atlas = method_atlas_match.group(3).lower()
region = method_atlas_match.group(4).lower()
level_type = "source"
level_detail = f"source_{method}_{atlas}_{region}"
elif atlas_match:
atlas = atlas_match.group(2).lower()
region = atlas_match.group(3).lower()
level_type = "source"
level_detail = f"source_{atlas}_{region}"
else:
# Legacy format: _Source_level_ or _Sensor_level_
level_match = re.search(r"(Source|Sensor)_level", filename, re.IGNORECASE)
level_type = level_match.group(1).lower() if level_match else "unknown"
level_detail = level_type
# Extract frequency band - try numerical range first
freq_match = re.search(r"(\d+\.?\d*-\d+\.?\d*)\s*Hz", filename)
if freq_match:
freq_band = freq_match.group(1)
else:
# Try to extract broadband text from filename parts
parts = filename.split("_")
for i, part in enumerate(parts):
if part.lower() == "level" and i + 1 < len(parts):
freq_band = parts[i + 1]
break
else:
freq_band = "unknown"
return {
"base_name": base_name,
"level_type": level_type,
"level_detail": level_detail,
"freq_band": freq_band,
"condition": f"{level_detail}_{freq_band}",
}
def process_subject_condition(args):
"""Process a single subject-condition combination and save matrices directly to disk."""
(
subject, condition, epoch_files, invert, calc_jpe, calc_pli, calc_pli_mst,
calc_aec, use_aecc, force_positive, jpe_st, calc_aec_mst, calc_power,
power_fs, calc_peak, peak_min, peak_max, calc_sampen, sampen_m, calc_apen,
apen_m, apen_r, calc_sv, sv_window, save_matrices, save_mst,
save_channel_averages, concat_aecc, has_headers, psd_method,
welch_window_ms, welch_overlap, calc_plt, plt_threshold_ms,
calc_lzc, lzc_threshold, calc_plt_mst, matrix_folders, save_epoch_metrics,
mt_bandwidth,
save_psd_avg, save_psd_perchannel, save_psd_plots, psd_folder,
total_power_range,
) = args
MST_MEASURES = ["degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean"]
try:
channel_results = defaultdict(lambda: defaultdict(list))
jpe_values, pe_values, pli_values, plt_values, aec_values = [], [], [], [], []
pli_mst_values = defaultdict(list)
aec_mst_values = defaultdict(list)
plt_mst_values = defaultdict(list)
power_values = defaultdict(list)
apen_values, sampen_values, lzc_values = [], [], []
sv_values = {}
channel_names = None
# PSD export accumulators (running mean over broadband epochs)
save_psd = save_psd_avg or save_psd_perchannel or save_psd_plots
psd_sum = None
psd_count = 0
psd_frequencies = None
successful_mst_epochs = 0
epoch_level_results = []
logger.info(f"Processing {subject} - {condition} ({len(epoch_files)} epochs)")
# Determine whether anything is actually computed per-epoch. If not (e.g. only
# spectral variability or only concatenated AEC is requested), we still need
# channel_names but can skip re-reading every epoch in this loop.
per_epoch_work = (
calc_jpe or calc_pli or calc_plt or calc_sampen or calc_apen or calc_lzc
or (calc_aec and not concat_aecc)
or ((calc_power or calc_peak or save_psd) and is_broadband_condition(condition))
)
for i, file_path in enumerate(epoch_files):
if MemoryMonitor.check_memory():
logger.warning(f"High memory usage detected while processing {subject}")
time.sleep(1)
try:
epoch_id = extract_epoch_id(file_path)
current_epoch_data = {"Subject": subject, "Condition": condition, "Epoch": epoch_id}
# --- Read Data (centralized in read_epoch_array) ---
data_values, file_channel_names = read_epoch_array(file_path, has_headers)
if channel_names is None:
channel_names = (
file_channel_names
if file_channel_names is not None
else [f"Channel_{j + 1}" for j in range(data_values.shape[1])]
)
# channel_names is now set; if there is no per-epoch work to do, stop
# after the first file rather than re-reading every epoch needlessly.
if not per_epoch_work:
del data_values
break
# Use level_detail (not level_type) so multi-atlas / multi-method source
# files don't collide on disk at the same band + epoch.
file_info = parse_epoch_filename(os.path.basename(file_path))
level_type = file_info["level_detail"]
freq_band = file_info["freq_band"]
# --- Spectral (Power / Peak) ---
need_spectral = (calc_power or calc_peak or save_psd) and is_broadband_condition(condition)
if need_spectral:
try:
# psd_kwargs = {"window_length_ms": welch_window_ms, "overlap_percent": welch_overlap} if psd_method == "welch" else {}
if psd_method == "welch":
psd_kwargs = {"window_length_ms": welch_window_ms, "overlap_percent": welch_overlap}
elif psd_method == "multitaper":
psd_kwargs = {"bandwidth": mt_bandwidth}
else:
psd_kwargs = {}
spectral_data = calculate_PSD(data=data_values, fs=power_fs, method=psd_method, **psd_kwargs)
if calc_power:
powers, channel_powers = calculate_power_bands(
spectral_data["frequencies"],
spectral_data["psd"],
total_power_range=total_power_range,
)
for measure, value in powers.items():
power_values[measure].append(value)
current_epoch_data[measure] = value
if save_channel_averages:
for band_name in FREQUENCY_BANDS:
for ch in range(len(channel_names)):
channel_results[channel_names[ch]][f"{band_name}_abs_power"].append(channel_powers[f"{band_name}_abs_power"][ch])
channel_results[channel_names[ch]][f"{band_name}_rel_power"].append(channel_powers[f"{band_name}_rel_power"][ch])
if calc_peak:
peak_freqs = calculate_avg_peak_frequency(spectral_data["frequencies"], spectral_data["psd"], freq_range=(peak_min, peak_max))
current_epoch_data["peak_frequency"] = np.nanmean(peak_freqs)
power_values["peak_frequency"].append(np.nanmean(peak_freqs))
power_values["channels_without_peak"].append(np.sum(np.isnan(peak_freqs)))
if save_channel_averages:
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["peak_frequency"].append(peak_freqs[ch])
if save_psd:
cur_psd = spectral_data["psd"]
if psd_sum is None:
psd_sum = cur_psd.astype(np.float64).copy()
psd_frequencies = spectral_data["frequencies"]
psd_count = 1
elif cur_psd.shape == psd_sum.shape:
psd_sum += cur_psd
psd_count += 1
else:
logger.warning(f"PSD shape mismatch in {epoch_id}; excluded from PSD average.")
del spectral_data
except Exception:
logger.exception(f"Error in spectral calculations for {epoch_id}")
# --- JPE / PE ---
if calc_jpe:
try:
jpe_matrix = calculate_jpe(data_values, n=4, st=jpe_st, invert=invert)
mask = ~np.eye(jpe_matrix.shape[0], dtype=bool)
mean_jpe = jpe_matrix[mask].mean()
jpe_values.append(mean_jpe)
current_epoch_data["JPE"] = mean_jpe
pe_values_array = calculate_pe(data_values, n=4, st=jpe_st)
pe_values.append(pe_values_array.mean())
current_epoch_data["PE"] = pe_values_array.mean()
if save_channel_averages:
channel_jpe = np.mean(jpe_matrix, axis=1)
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["jpe"].append(channel_jpe[ch])
channel_results[channel_names[ch]]["pe"].append(pe_values_array[ch])
if save_matrices and matrix_folders:
save_connectivity_matrix(jpe_matrix, matrix_folders["jpe"], subject, freq_band, f"jpe_{epoch_id}", channel_names, level_type)
except Exception:
logger.exception(f"Error calculating JPE/PE for {epoch_id}")
# --- SampEn ---
if calc_sampen:
try:
sampen_values_ch = calculate_sampen_for_channels(data_values, m=sampen_m)
sampen_values.append(np.nanmean(sampen_values_ch))
current_epoch_data["SampEn"] = np.nanmean(sampen_values_ch)
if save_channel_averages:
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["sampen"].append(sampen_values_ch[ch])
except Exception:
logger.exception(f"Error calculating SampEn for {epoch_id}")
# --- ApEn ---
if calc_apen:
try:
apen_values_ch = calculate_apen_for_channels(data_values, m=apen_m, r=apen_r)
apen_values.append(np.nanmean(apen_values_ch))
current_epoch_data["ApEn"] = np.nanmean(apen_values_ch)
if save_channel_averages:
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["apen"].append(apen_values_ch[ch])
except Exception:
logger.exception(f"Error calculating ApEn for {epoch_id}")
# --- LZC ---
if calc_lzc:
try:
lzc_values_ch = calculate_lzc_for_channels(data_values, threshold=lzc_threshold)
lzc_values.append(np.nanmean(lzc_values_ch))
current_epoch_data["LZC"] = np.nanmean(lzc_values_ch)
if save_channel_averages:
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["lzc"].append(lzc_values_ch[ch])
except Exception:
logger.exception(f"Error calculating LZC for {epoch_id}")
# --- PLI ---
if calc_pli:
try:
pli_matrix = calculate_pli(data_values)
mask = ~np.eye(pli_matrix.shape[0], dtype=bool)
current_epoch_data["PLI"] = pli_matrix[mask].mean()
pli_values.append(current_epoch_data["PLI"])
if save_channel_averages:
channel_pli = np.mean(pli_matrix, axis=1)
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["pli"].append(channel_pli[ch])
if save_matrices and matrix_folders:
save_connectivity_matrix(pli_matrix, matrix_folders["pli"], subject, freq_band, f"pli_{epoch_id}", channel_names, level_type)
if calc_pli_mst:
mst_measures, mst_matrix, success = calculate_mst_measures(pli_matrix)
if success:
for measure, value in mst_measures.items():
pli_mst_values[measure].append(value)
current_epoch_data[f"PLI_MST_{measure}"] = value
if save_mst and matrix_folders:
save_connectivity_matrix(mst_matrix + mst_matrix.T, matrix_folders["pli_mst"], subject, freq_band, f"pli_mst_{epoch_id}", channel_names, level_type)
except Exception:
logger.exception(f"Error calculating PLI for {epoch_id}")
# --- PLT ---
if calc_plt:
try:
plt_matrix = PLT(data_values, fs=power_fs, threshold_ms=plt_threshold_ms)
mask = ~np.eye(plt_matrix.shape[0], dtype=bool)
current_epoch_data["PLT"] = plt_matrix[mask].mean()
plt_values.append(current_epoch_data["PLT"])
if save_channel_averages:
channel_plt = np.mean(plt_matrix, axis=1)
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["plt"].append(channel_plt[ch])
if save_matrices and matrix_folders:
save_connectivity_matrix(plt_matrix, matrix_folders["plt"], subject, freq_band, f"plt_{epoch_id}", channel_names, level_type)
if calc_plt_mst:
mst_measures, mst_matrix, success = calculate_mst_measures(plt_matrix)
if success:
for measure, value in mst_measures.items():
plt_mst_values[measure].append(value)
current_epoch_data[f"PLT_MST_{measure}"] = value
if save_mst and matrix_folders:
save_connectivity_matrix(mst_matrix + mst_matrix.T, matrix_folders["plt_mst"], subject, freq_band, f"plt_mst_{epoch_id}", channel_names, level_type)
except Exception:
logger.exception(f"Error calculating PLT for {epoch_id}")
# --- AEC (Standard Epoch) ---
if calc_aec and not concat_aecc:
try:
aec_matrix = calculate_aecc(data_values, orthogonalize=use_aecc, force_positive=force_positive)
mask = ~np.eye(aec_matrix.shape[0], dtype=bool)
current_epoch_data["AEC"] = aec_matrix[mask].mean()
aec_values.append(current_epoch_data["AEC"])
if save_channel_averages:
channel_aec = np.mean(aec_matrix, axis=1)
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["aec"].append(channel_aec[ch])
if save_matrices and matrix_folders:
save_connectivity_matrix(aec_matrix, matrix_folders["aec"], subject, freq_band, f"aec_{epoch_id}", channel_names, level_type)
if calc_aec_mst:
mst_measures, mst_matrix, success = calculate_mst_measures(aec_matrix)
if success:
successful_mst_epochs += 1
for measure, value in mst_measures.items():
aec_mst_values[measure].append(value)
current_epoch_data[f"AEC_MST_{measure}"] = value
if save_mst and matrix_folders:
save_connectivity_matrix(mst_matrix + mst_matrix.T, matrix_folders["aec_mst"], subject, freq_band, f"aec_mst_{epoch_id}", channel_names, level_type)
except Exception:
logger.exception(f"Error calculating AEC for {epoch_id}")
if save_epoch_metrics:
epoch_level_results.append(current_epoch_data)
except Exception:
logger.exception(f"Error processing file {os.path.basename(file_path)}")
continue
# --- AEC (Concatenated Edge Case) ---
if calc_aec and concat_aecc and epoch_files:
try:
all_data = []
first_nbytes = None
for aec_epoch_file_path in epoch_files:
try:
arr, _ = read_epoch_array(aec_epoch_file_path, has_headers) # now linearly detrended
if first_nbytes is None:
first_nbytes = arr.nbytes
if not MemoryMonitor.check_concatenation_safety(first_nbytes, len(epoch_files)):
logger.warning("Memory may be insufficient for concatenated AEC; proceeding cautiously")
all_data.append(arr)
except Exception:
logger.exception(f"Error reading {os.path.basename(aec_epoch_file_path)} for concatenated AEC")
if all_data:
data_values = np.concatenate(all_data, axis=0)
del all_data
aec_matrix = calculate_aecc(data_values, orthogonalize=use_aecc, force_positive=force_positive)
mask = ~np.eye(aec_matrix.shape[0], dtype=bool)
mean_aec = aec_matrix[mask].mean()
aec_values.append(mean_aec)
file_info = parse_epoch_filename(os.path.basename(epoch_files[0]))
level_type = file_info["level_detail"]
freq_band = file_info["freq_band"]
concat_epoch_data = {"Subject": subject, "Condition": condition, "Epoch": "Concatenated", "AEC": mean_aec}
if save_channel_averages:
channel_aec = np.mean(aec_matrix, axis=1)
for ch in range(len(channel_names)):
channel_results[channel_names[ch]]["aec"].append(channel_aec[ch])
if save_matrices and matrix_folders:
save_connectivity_matrix(aec_matrix, matrix_folders["aec"], subject, freq_band, "aec_concatenated", channel_names, level_type)
if calc_aec_mst:
mst_measures, mst_matrix, success = calculate_mst_measures(aec_matrix)
if success:
successful_mst_epochs += 1
for measure, value in mst_measures.items():
aec_mst_values[measure].append(value)
concat_epoch_data[f"AEC_MST_{measure}"] = value
if save_mst and matrix_folders:
save_connectivity_matrix(mst_matrix + mst_matrix.T, matrix_folders["aec_mst"], subject, freq_band, "aec_mst_concatenated", channel_names, level_type)
if save_epoch_metrics:
epoch_level_results.append(concat_epoch_data)
except Exception:
logger.exception("Error processing concatenated AEC")
# --- Spectral Variability (computed once on concatenated broadband data) ---
if calc_sv and is_broadband_condition(condition) and epoch_files:
try:
all_data_sv = []
first_nbytes = None
for sv_fp in epoch_files:
try:
sv_arr, _ = read_epoch_array(sv_fp, has_headers) # now linearly detrended
if first_nbytes is None:
first_nbytes = sv_arr.nbytes
if not MemoryMonitor.check_concatenation_safety(first_nbytes, len(epoch_files)):
logger.warning("Memory may be insufficient for concatenated SV; proceeding cautiously")
all_data_sv.append(sv_arr)
except Exception:
logger.exception(f"Error reading {os.path.basename(sv_fp)} for SV")
if all_data_sv:
concatenated_sv = np.concatenate(all_data_sv, axis=0)
del all_data_sv
sv_results = calculate_spectral_variability(
concatenated_sv,
fs=power_fs,
window_length=sv_window,
total_power_range=total_power_range,
)
del concatenated_sv
if sv_results:
for band_name, values in sv_results.items():
sv_values[f"sv_{band_name}"] = np.nanmean(values)
if save_channel_averages:
for ch in range(len(channel_names)):
channel_results[channel_names[ch]][f"sv_{band_name}"] = [values[ch]]
except Exception:
logger.exception("Error in spectral variability calculation")
# --- Save PSD export (epoch-averaged, broadband only) ---
if save_psd and is_broadband_condition(condition) and psd_sum is not None and psd_count > 0:
try:
mean_psd = psd_sum / psd_count
psd_info = parse_epoch_filename(os.path.basename(epoch_files[0]))
save_psd_outputs(
mean_psd,
psd_frequencies,
channel_names,
subject,
condition,
psd_info["freq_band"],
psd_info["level_detail"],
psd_folder,
save_avg=save_psd_avg,
save_perchannel=save_psd_perchannel,
save_plots=save_psd_plots,
)
except Exception:
logger.exception("Error saving PSD export")
# --- Channel-level averages across epochs ---
channel_averages = None
if save_channel_averages and channel_names:
channel_averages = {}
for channel in channel_names:
channel_averages[channel] = {
measure: np.mean(values) for measure, values in channel_results[channel].items()
}
# --- Build whole-brain averaged results dict (Legacy Wide Format) ---
used_epoch_indices = [extract_epoch_id(fp).replace("Epoch_", "") for fp in epoch_files]
results = {
"avg_jpe": np.mean(jpe_values) if jpe_values else np.nan,
"avg_pe": np.mean(pe_values) if pe_values else np.nan,
"avg_pli": np.mean(pli_values) if pli_values else np.nan,
"avg_plt": np.mean(plt_values) if plt_values else np.nan,
"avg_aec": np.mean(aec_values) if aec_values else np.nan,
"avg_sampen": np.mean(sampen_values) if sampen_values else np.nan,
"avg_apen": np.mean(apen_values) if apen_values else np.nan,
"avg_lzc": np.mean(lzc_values) if lzc_values else np.nan,
"n_epochs": len(epoch_files),
"used_epochs": ", ".join(used_epoch_indices),
"channel_names": channel_names if channel_names else [],
"channel_averages": channel_averages,
}
# Average MST values for wide format
for mst_dict, prefix in [(aec_mst_values, "aec"), (pli_mst_values, "pli"), (plt_mst_values, "plt")]:
for measure in MST_MEASURES:
if mst_dict[measure]:
results[f"{prefix}_mst_{measure}"] = mst_dict[measure][0] if (prefix == "aec" and concat_aecc) else np.mean(mst_dict[measure])
results[f"{prefix}_mst_{measure}_valid_epochs"] = 1 if (prefix == "aec" and concat_aecc) else len(mst_dict[measure])
else:
results[f"{prefix}_mst_{measure}"] = np.nan
results[f"{prefix}_mst_{measure}_valid_epochs"] = 0
results["aec_mst_successful_epochs"] = successful_mst_epochs
results["aec_mst_total_epochs"] = 1 if concat_aecc else len(epoch_files)
results["pli_mst_successful_epochs"] = len(pli_mst_values["degree"]) if pli_mst_values["degree"] else 0
results["pli_mst_total_epochs"] = len(epoch_files)
results["plt_mst_successful_epochs"] = len(plt_mst_values["degree"]) if plt_mst_values["degree"] else 0
results["plt_mst_total_epochs"] = len(epoch_files)
if calc_power and power_values:
for band_name in FREQUENCY_BANDS:
results[f"{band_name}_abs_power"] = np.mean(power_values[f"{band_name}_abs_power"]) if power_values[f"{band_name}_abs_power"] else np.nan
results[f"{band_name}_rel_power"] = np.mean(power_values[f"{band_name}_rel_power"]) if power_values[f"{band_name}_rel_power"] else np.nan
if calc_peak:
results["peak_frequency"] = np.mean(power_values["peak_frequency"]) if power_values["peak_frequency"] else np.nan
results["channels_without_peak"] = np.mean(power_values["channels_without_peak"]) if power_values["channels_without_peak"] else np.nan
if calc_sv:
for band_name in FREQUENCY_BANDS:
results[f"sv_{band_name}"] = sv_values.get(f"sv_{band_name}", np.nan)
return subject, condition, results, epoch_level_results
except Exception:
logger.exception(f"Error processing {subject} - {condition}")
return subject, condition, {"n_epochs": 0, "used_epochs": "Error", "channel_names": [], "channel_averages": None}, []
def process_all_subjects(
grouped_files, invert, n_threads, calc_jpe, calc_pli, calc_pli_mst, calc_aec,
use_aecc, force_positive=True, jpe_st=1, calc_aec_mst=False, calc_power=False,
power_fs=256, calc_peak=False, peak_min=3, peak_max=13, calc_sampen=False,
sampen_m=2, calc_apen=False, apen_m=1, apen_r=0.25, calc_lzc=False,
lzc_threshold="median", calc_sv=False, sv_window=1000, save_matrices=False,
save_mst=False, save_channel_averages=False, concat_aecc=False, has_headers=True,
psd_method="multitaper", welch_window_ms=1000, welch_overlap=50, calc_plt=False,
plt_threshold_ms=30, calc_plt_mst=False, limit_epochs=False, max_epochs=5,
random_seed=111, progress_callback=None, matrix_folders=None, save_epoch_metrics=True,
mt_bandwidth=None,
save_psd_avg=False, save_psd_perchannel=False, save_psd_plots=False, psd_folder=None,
total_power_range=None,
):
process_args = []
for subject, conditions in grouped_files.items():
for condition, epoch_files in conditions.items():
files_to_process = epoch_files
if limit_epochs and len(epoch_files) > max_epochs:
sorted_files = sorted(epoch_files)
subject_specific_seed = random_seed + sum(ord(c) for c in str(subject))
rng = random.Random(subject_specific_seed)
files_to_process = sorted(rng.sample(sorted_files, max_epochs))
process_args.append((
subject, condition, files_to_process, invert, calc_jpe, calc_pli, calc_pli_mst,
calc_aec, use_aecc, force_positive, jpe_st, calc_aec_mst, calc_power,
power_fs, calc_peak, peak_min, peak_max, calc_sampen, sampen_m, calc_apen,
apen_m, apen_r, calc_sv, sv_window, save_matrices, save_mst,
save_channel_averages, concat_aecc, has_headers, psd_method,
welch_window_ms, welch_overlap, calc_plt, plt_threshold_ms,
calc_lzc, lzc_threshold, calc_plt_mst, matrix_folders, save_epoch_metrics,
mt_bandwidth,
save_psd_avg, save_psd_perchannel, save_psd_plots, psd_folder,
total_power_range,
))
total_tasks = len(process_args)
results = defaultdict(dict)
all_epoch_results = []
completed = 0
for i in range(0, len(process_args), BATCH_SIZE):
batch = process_args[i : i + BATCH_SIZE]
batch_results = process_batch(batch, n_threads)
for subject, condition, result, epoch_results in batch_results:
results[subject][condition] = result
all_epoch_results.extend(epoch_results)
completed += 1
if progress_callback: progress_callback(completed / total_tasks * 100)
return dict(results), all_epoch_results
def process_batch(batch_args, n_threads):
"""Process a batch of subjects using multiprocessing with fallback."""
try:
with Pool(processes=n_threads, maxtasksperchild=1) as pool:
return list(pool.imap_unordered(process_subject_condition, batch_args))
except Exception:
logger.exception("Pool processing failed, falling back to single thread")
return [process_subject_condition(args) for args in batch_args]
def group_epochs_by_condition(folder_path, folder_ext):
"""Group epoch files by their base name and condition.
Only processes folders containing valid epoch files
Returns a dictionary: {base_name: {condition: [epoch_files]}}.
"""
grouped_files = defaultdict(lambda: defaultdict(list))
# Get immediate subdirectories
try:
subdirs = [
d for d in os.listdir(folder_path) if os.path.isdir(os.path.join(folder_path, d)) and d.endswith(folder_ext)
]
except Exception:
sg.popup_error("Error accessing directory")
return grouped_files
if not subdirs:
sg.popup_error(f"No folders ending with '{folder_ext}' found in the selected directory.")
return grouped_files
for subdir in subdirs:
subdir_path = os.path.join(folder_path, subdir)
# Look for epoch files in this directory
for file in os.listdir(subdir_path):
if file.startswith((".", "._")):
continue
if not file.endswith(".txt"):
continue
# Check if file matches epoch pattern
is_legacy_format = "_level_" in file
is_new_atlas_format = bool(
re.search(r"_(Source|Sensor)_(beamformer|sLORETA|eLORETA|dSPM|MNE)_(desikan|bna|aal2|aal3)_", file, re.IGNORECASE)
or re.search(r"_(Source|Sensor)_(desikan|bna|aal2|aal3)_(cortical|full)_", file, re.IGNORECASE)
)
if (is_legacy_format or is_new_atlas_format) and ("_Epoch_" in file or "_Epoch" in file):
try:
file_info = parse_epoch_filename(file)
full_path = os.path.join(subdir_path, file)
base_name = subdir.replace(folder_ext, "") # Use folder name without extension
condition = file_info["condition"]
grouped_files[base_name][condition].append(full_path)
except Exception:
print(f"Skipping file {file}")
continue
# Print summary of what was found
print("\nFound the following data:")
for base_name, conditions in grouped_files.items():
print(f"\nSubject: {base_name}")
for condition, files in conditions.items():
print(f" {condition}: {len(files)} epochs")
found_bands = set()
unknown_conditions = set()
has_broadband = False
for conditions in grouped_files.values():
for condition in conditions:
band = extract_freq_band(condition)
if band != "unknown":
found_bands.add(band)
if band == "broadband":
has_broadband = True
else:
unknown_conditions.add(condition)
logger.info("Found the following frequency bands in the data:")
for band in sorted(found_bands):
logger.info(f" - {band} ({FREQUENCY_BANDS[band]['pattern']})")
if has_broadband:
logger.info("Broadband epochs are present - spectral calculations will be performed on these epochs")
else:
logger.warning("No broadband epochs found - spectral calculations will be skipped")
if unknown_conditions:
logger.warning("Found conditions with unrecognized frequency bands:")
for cond in sorted(unknown_conditions):
logger.warning(f" - {cond}")
return grouped_files
def save_results_to_excel(
results_dict,
output_path,
invert,
calc_pli_mst,
epoch_level_data=None,
calc_jpe=True,
calc_pli=True,
calc_plt=False,
calc_plt_mst=False,
calc_aec=False,
use_aecc=False,
force_positive=True,
calc_aec_mst=False,
calc_power=False,
power_fs=256,
calc_peak=False,
peak_min=None,
peak_max=None,
calc_sampen=False,
calc_apen=False,
calc_lzc=False,
lzc_threshold="median",
calc_sv=False,
save_channel_averages=False,
concat_aecc=False,
has_headers=True,
sv_window=None,
psd_method="multitaper",
welch_window_ms=None,
welch_overlap=None,
mt_bandwidth=None,
plt_threshold_ms=None,
save_psd_avg=False,
save_psd_perchannel=False,
save_psd_plots=False,
total_power_range=None,
reference_range_source=None,
):
"""
Save results to Excel with organized columns by condition.
- Moves 'used_epochs' to a specific Metadata sheet.
- Ensures epoch numbers are sorted numerically (1, 2, 10) instead of alphabetically (1, 10, 2).
"""
with pd.ExcelWriter(output_path, engine="openpyxl") as writer:
# 1) Gather all unique conditions
all_conditions = set()
for subject_data in results_dict.values():
all_conditions.update(subject_data.keys())
# 2) Build a list of primary columns; start with 'subject'
columns = ["subject"]
for condition in sorted(all_conditions):
# --- Complexity Measures (JPE/PE) ---
if calc_jpe:
measure_name = "jpe_inv" if invert else "jpe"
columns.append(f"{condition}_avg_{measure_name}")
columns.append(f"{condition}_avg_pe")
# --- PLI ---
if calc_pli:
columns.append(f"{condition}_avg_pli")
if calc_pli_mst:
mst_measures = [
"degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean",
]
for mm in mst_measures:
columns.append(f"{condition}_pli_mst_{mm}")
columns.append(f"{condition}_pli_mst_successful_epochs")
columns.append(f"{condition}_pli_mst_total_epochs")
# --- PLT ---
if calc_plt:
columns.append(f"{condition}_avg_plt")
if calc_plt_mst:
mst_measures = [
"degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean",
]
for mm in mst_measures:
columns.append(f"{condition}_plt_mst_{mm}")
columns.append(f"{condition}_plt_mst_successful_epochs")
columns.append(f"{condition}_plt_mst_total_epochs")
# --- AEC ---
if calc_aec:
columns.append(f"{condition}_avg_aec")
if calc_aec_mst:
mst_measures = [
"degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean",
]
for mm in mst_measures:
columns.append(f"{condition}_aec_mst_{mm}")
columns.append(f"{condition}_aec_mst_successful_epochs")
columns.append(f"{condition}_aec_mst_total_epochs")
# --- SampEn, ApEn and LZC ---
if calc_sampen:
columns.append(f"{condition}_avg_sampen")
if calc_apen:
columns.append(f"{condition}_avg_apen")
if calc_lzc:
columns.append(f"{condition}_avg_lzc")
is_broadband_cond = is_broadband_condition(condition)
# --- Power band measures - only for broadband conditions ---
if calc_power and is_broadband_cond:
for band_name in FREQUENCY_BANDS:
if band_name.lower() != "broadband": # Skip broadband
columns.extend([f"{condition}_{band_name}_abs_power", f"{condition}_{band_name}_rel_power"])
# --- Peak frequency ---
if calc_peak and is_broadband_cond:
columns.append(f"{condition}_peak_frequency")
columns.append(f"{condition}_channels_without_peak")
# --- Spectral Variability ---
if calc_sv and is_broadband_cond:
for band_name in FREQUENCY_BANDS:
if band_name.lower() == "broadband":
continue
columns.append(f"{condition}_sv_{band_name}")
# Always add epoch count for each condition
columns.append(f"{condition}_n_epochs")
# 3) Build rows of data
rows = []
for subject, conditions in results_dict.items():
row = {"subject": subject}
for condition in sorted(all_conditions):
data_for_condition = conditions.get(condition, {})
# JPE/PE measures
if calc_jpe:
measure_name = "jpe_inv" if invert else "jpe"
row[f"{condition}_avg_{measure_name}"] = data_for_condition.get("avg_jpe", np.nan)
row[f"{condition}_avg_pe"] = data_for_condition.get("avg_pe", np.nan)
# PLI
if calc_pli:
row[f"{condition}_avg_pli"] = data_for_condition.get("avg_pli", np.nan)
if calc_pli_mst:
mst_measures = [
"degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean",
]
for mm in mst_measures:
row[f"{condition}_pli_mst_{mm}"] = data_for_condition.get(f"pli_mst_{mm}", np.nan)
row[f"{condition}_pli_mst_successful_epochs"] = data_for_condition.get("pli_mst_successful_epochs", np.nan)
row[f"{condition}_pli_mst_total_epochs"] = data_for_condition.get("pli_mst_total_epochs", np.nan)
# PLT
if calc_plt:
row[f"{condition}_avg_plt"] = data_for_condition.get("avg_plt", np.nan)
if calc_plt_mst:
mst_measures = [
"degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean",
]
for mm in mst_measures:
row[f"{condition}_plt_mst_{mm}"] = data_for_condition.get(f"plt_mst_{mm}", np.nan)
row[f"{condition}_plt_mst_successful_epochs"] = data_for_condition.get("plt_mst_successful_epochs", np.nan)
row[f"{condition}_plt_mst_total_epochs"] = data_for_condition.get("plt_mst_total_epochs", np.nan)
# AEC
if calc_aec:
row[f"{condition}_avg_aec"] = data_for_condition.get("avg_aec", np.nan)
if calc_aec_mst:
mst_measures = [
"degree", "eccentr", "betweenness", "kappa", "r",
"diameter", "leaf", "hierarchy", "teff", "asp", "ref", "mean",
]
for mm in mst_measures:
row[f"{condition}_aec_mst_{mm}"] = data_for_condition.get(f"aec_mst_{mm}", np.nan)
row[f"{condition}_aec_mst_successful_epochs"] = data_for_condition.get("aec_mst_successful_epochs", np.nan)
row[f"{condition}_aec_mst_total_epochs"] = data_for_condition.get("aec_mst_total_epochs", np.nan)
# SampEn, ApEn and LZC
if calc_sampen:
row[f"{condition}_avg_sampen"] = data_for_condition.get("avg_sampen", np.nan)
if calc_apen:
row[f"{condition}_avg_apen"] = data_for_condition.get("avg_apen", np.nan)
if calc_lzc:
row[f"{condition}_avg_lzc"] = data_for_condition.get("avg_lzc", np.nan)
# Power & Peak Frequency
is_broadband_cond = is_broadband_condition(condition)
if calc_power and is_broadband_cond:
for band_name in FREQUENCY_BANDS:
if band_name.lower() == "broadband":
continue
abs_key = f"{band_name}_abs_power"
rel_key = f"{band_name}_rel_power"
row[f"{condition}_{band_name}_abs_power"] = data_for_condition.get(abs_key, np.nan)
row[f"{condition}_{band_name}_rel_power"] = data_for_condition.get(rel_key, np.nan)
if calc_peak and is_broadband_cond:
row[f"{condition}_peak_frequency"] = data_for_condition.get("peak_frequency", np.nan)
row[f"{condition}_channels_without_peak"] = data_for_condition.get("channels_without_peak", np.nan)
if calc_sv and is_broadband_cond:
for band_name in FREQUENCY_BANDS:
if band_name.lower() == "broadband":
continue
sv_key = f"sv_{band_name}"
row[f"{condition}_sv_{band_name}"] = data_for_condition.get(sv_key, np.nan)
# n_epochs
row[f"{condition}_n_epochs"] = data_for_condition.get("n_epochs", 0)
rows.append(row)
# 4) Create the DataFrame, reorder columns, and export
df = pd.DataFrame(rows)
df = df[columns] # Force the column order we built above
df.to_excel(writer, sheet_name="Whole Brain Results", index=False)
# Add analysis information sheet
info_data = {
"Parameter": [
"Analysis Date",
"JPE Inversion",
"PLI MST Calculated",
"AEC Type",
"AEC Concatenated Epochs",
"AEC MST Calculated",
"PLT Calculated",
"PLT MST Calculated",
"PLT Threshold (ms)",
"AEC Force Positive",
"Power Bands Calculated",
"Relative Power Reference Range (Hz)",
"Reference Range Source",
"Sampling Frequency (Hz)",
"PSD Method",
"Welch Window Length (ms)",
"Welch Overlap (%)",
"Multitaper Bandwidth (Hz)",
"Peak Frequency Analysis",
"Peak Frequency Range (Hz)",
"Sample Entropy Calculated",
"Approximate Entropy Calculated",
"Lempel-Ziv Complexity Calculated",
"LZC Binarization Threshold",
"Spectral Variability Calculated",
"Spectral Variability Window (ms)",
"Channel Averages Calculated",
"Channel Names Source",
"PSD Averaged CSV Saved",
"PSD Per-Channel CSV Saved",
"PSD Plots Saved",
],
"Value": [
datetime.now().strftime("%Y-%m-%d %H:%M:%S"),
"Yes" if invert else "No",
"Yes" if calc_pli_mst else "No",
"AECc (orthogonalized)" if calc_aec and use_aecc else "AEC" if calc_aec else "Not calculated",
"Yes" if concat_aecc else "No",
"Yes" if calc_aec_mst else "No",
"Yes" if calc_plt else "No",
"Yes" if calc_plt_mst else "No",
str(plt_threshold_ms) if calc_plt else "N/A",
"Yes" if force_positive else "No",
"Yes" if calc_power else "No",
(
f"{total_power_range[0]:g}-{total_power_range[1]:g}"
if (calc_power or calc_sv) and total_power_range is not None
else "N/A"
),
(
reference_range_source
if (calc_power or calc_sv) and reference_range_source
else "N/A"
),
str(power_fs),
psd_method,
str(welch_window_ms) if psd_method == "welch" else "N/A",
str(welch_overlap) if psd_method == "welch" else "N/A",
("Auto (8 x sfreq/n_samples)" if mt_bandwidth is None else str(mt_bandwidth))
if psd_method == "multitaper" else "N/A",
"Yes" if calc_peak else "No",
f"{peak_min}-{peak_max}" if calc_peak else "N/A",
"Yes" if calc_sampen else "No",
"Yes" if calc_apen else "No",
"Yes" if calc_lzc else "No",
str(lzc_threshold) if calc_lzc else "N/A",
"Yes" if calc_sv else "No",
str(sv_window) if calc_sv else "N/A",
"Yes" if save_channel_averages else "No",
"File Headers" if has_headers else "Auto-generated",
"Yes" if save_psd_avg else "No",
"Yes" if save_psd_perchannel else "No",
"Yes" if save_psd_plots else "No",
],
}
info_df = pd.DataFrame(info_data)
info_df.to_excel(writer, sheet_name="Analysis Information", index=False)
# Save channel-level averages if requested
if save_channel_averages:
all_channels = set()
for conditions in results_dict.values():
for result in conditions.values():
if result.get("channel_averages"):
all_channels.update(result["channel_averages"].keys())
channel_rows = []
measures_with_data = set()
for subject, conditions in results_dict.items():
for channel in sorted(all_channels):
row = {"subject": subject, "channel": channel}
for condition in sorted(all_conditions):
if condition in conditions and conditions[condition].get("channel_averages"):
channel_data = conditions[condition]["channel_averages"].get(channel, {})
for measure, value in channel_data.items():
column_name = f"{condition}_{measure}"
row[column_name] = value
if not pd.isna(value):
measures_with_data.add(column_name)
channel_rows.append(row)
if channel_rows:
df_channels = pd.DataFrame(channel_rows)
base_cols = ["subject", "channel"]
data_cols = sorted(measures_with_data, key=lambda x: (x.split("_")[0], x))
column_order = base_cols + data_cols
df_channels = df_channels[column_order]
df_channels.to_excel(writer, sheet_name="Channel Averages", index=False)
# --- Save metadata ---
metadata_rows = []
for subject, conditions in results_dict.items():
for condition in sorted(all_conditions):
if condition in conditions:
# Get the string (e.g., "1, 10, 2")
used_str = conditions[condition].get("used_epochs", "")
# --- Sort the epochs numerically ---
if used_str:
try:
# Split, convert to integers for proper sorting, then back to string
# We use (0, int(x)) for numbers and (1, x) for non-numbers to handle "?" gracefully
parts = [p.strip() for p in used_str.split(',')]
parts.sort(key=lambda x: (0, int(x)) if x.isdigit() else (1, x))
used_str = ", ".join(parts)
except Exception:
pass # Keep original string if logic fails
# ----------------------------------------
meta_row = {
"subject": subject,
"condition": condition,
"n_epochs_used": conditions[condition].get("n_epochs", 0),
"used_epochs": used_str,
}
if "channel_names" in conditions[condition]:
meta_row["n_channels"] = len(conditions[condition]["channel_names"])
meta_row["channels"] = ", ".join(conditions[condition]["channel_names"])
metadata_rows.append(meta_row)
if metadata_rows:
metadata_df = pd.DataFrame(metadata_rows)
desired_cols = ["subject", "condition", "n_epochs_used", "used_epochs", "n_channels", "channels"]
existing_cols = [c for c in desired_cols if c in metadata_df.columns]
metadata_df = metadata_df[existing_cols]
metadata_df.to_excel(writer, sheet_name="Epoch & Channel Info", index=False)
# --- Save Per-Epoch Metrics (Long Format) ---
if epoch_level_data:
df_epochs = pd.DataFrame(epoch_level_data)
# Reorder columns to ensure Subject, Condition, Epoch are first
base_cols = ['Subject', 'Condition', 'Epoch']
existing_base = [c for c in base_cols if c in df_epochs.columns]
other_cols = [c for c in df_epochs.columns if c not in existing_base]
df_epochs = df_epochs[existing_base + sorted(other_cols)]
# Sort rows logically. imap_unordered returns conditions in completion order,
# and within a condition rows follow filesystem order, so the raw frame is
# unordered. Sort by Subject, Condition, then numeric epoch number; any
# non-numeric label (e.g. "Concatenated") sorts last within its condition.
if {'Subject', 'Condition', 'Epoch'}.issubset(df_epochs.columns):
epoch_nums = df_epochs['Epoch'].astype(str).str.extract(r'(\d+)')[0]
df_epochs['_is_special'] = epoch_nums.isna().astype(int)
df_epochs['_epoch_num'] = pd.to_numeric(epoch_nums, errors='coerce')
df_epochs = (
df_epochs
.sort_values(by=['Subject', 'Condition', '_is_special', '_epoch_num'])
.drop(columns=['_is_special', '_epoch_num'])
.reset_index(drop=True)
)
df_epochs.to_excel(writer, sheet_name="Epoch-Level Results", index=False)
logger.info(f"Results saved to: {output_path}")
print(f"\nResults saved to {output_path}")
def main():
"""
Main entry point for the EEG-Pype quantitative analysis GUI.
Handles the event loop, thorough parameter validation, directory creation,
multiprocessing execution, and final result compilation.
"""
window = create_gui()
# Prime the live frequency-resolution readout with the default field values.
initial_event, initial_values = window.read(timeout=0)
if initial_values is not None:
window["-FREQ_RES-"].update(compute_resolution_readout(initial_values))
# Tracks the last folder we probed for a preprocessing pass-band, so that
# typing a path by hand does not re-read the config on every keystroke.
last_passband_folder = None
# --- GUI Event Loop ---
while True:
event, values = window.read()
if event in (sg.WIN_CLOSED, "Exit"):
break
# Live-update the frequency-resolution readout when any relevant field changes.
if event in ("-PSD_METHOD-", "-POWER_FS-", "-EPOCH_LENGTH-", "-WELCH_WINDOW-", "-MT_BANDWIDTH-"):
window["-FREQ_RES-"].update(compute_resolution_readout(values))
window.refresh()
# When a new folder is picked, seed the relative-power reference range
# from the pass-band the preprocessing script recorded, if we can find it.
if event == "-FOLDER-":
candidate = values["-FOLDER-"]
if candidate and candidate != last_passband_folder and os.path.isdir(candidate):
last_passband_folder = candidate
detected = discover_preprocessing_passband(candidate)
if detected is not None:
window["-TOTAL_MIN-"].update(f"{detected['low']:g}")
window["-TOTAL_MAX-"].update(f"{detected['high']:g}")
window["-PASSBAND_INFO-"].update(
f"Detected pass-band: {detected['low']:g}-{detected['high']:g} Hz",
text_color="#1A6B3C",
)
else:
window["-PASSBAND_INFO-"].update(
"No unambiguous preprocessing pass-band found",
text_color="#8A6D1F",
)
window.refresh()
if event == "Process":
# --- 1. Basic Setup & File Path Extraction ---
folder_path = values["-FOLDER-"]
folder_ext = values["-EXTENSION-"].strip()
# Initialize logging immediately for this run
log_file = setup_logging(folder_path)
logger.info("=== Starting new analysis run ===")
logger.info(f"Folder path: {folder_path}")
logger.info(f"Extension: {folder_ext}")
logger.info(f"Processing files with{'out' if not values['-HAS_HEADERS-'] else ''} headers")
if not values["-HAS_HEADERS-"]:
logger.info("Channel names will be auto-generated")
if not folder_path or not folder_ext:
sg.popup_error("Please select a folder and specify the folder extension.")
continue
try:
validate_frequency_bands()
except ValueError:
sg.popup_error("Invalid frequency band configuration in FREQUENCY_BANDS.")
continue
# --- 2. Parameter Extraction & Validation ---
# Threads
try:
n_threads = int(values["-THREADS-"])
if n_threads < 1 or n_threads > cpu_count():
raise ValueError(f"Number of threads must be between 1 and {cpu_count()}")
except ValueError:
sg.popup_error("Invalid number of threads.")
continue
# Matrix Export Configuration
save_matrices = values["-SAVE_MATRICES-"]
matrix_folder = values["-MATRIX_FOLDER-"]
save_mst = values["-SAVE_MST-"]
mst_folder = values["-MST_FOLDER-"]
save_epoch_metrics = values.get("-SAVE_EPOCH_METRICS-", True)
# PSD Export Configuration
save_psd_avg = values["-SAVE_PSD_AVG-"]
save_psd_perchannel = values["-SAVE_PSD_PERCHANNEL-"]
save_psd_plots = values["-SAVE_PSD_PLOTS-"]
psd_folder_name = values["-PSD_FOLDER-"].strip()
save_psd_any = save_psd_avg or save_psd_perchannel or save_psd_plots
if save_psd_any and not psd_folder_name:
sg.popup_error("Please specify a folder name for saving PSD output.")
continue
if save_matrices and not matrix_folder.strip():
sg.popup_error("Please specify a folder name for saving connectivity matrices.")
continue
if save_mst and not mst_folder.strip():
sg.popup_error("Please specify a folder name for saving MST matrices.")
continue
# Epoch Limiting Configuration
limit_epochs = values["-LIMIT_EPOCHS-"]
max_epochs = 5
random_seed = 111
if limit_epochs:
try:
max_epochs = int(values["-MAX_EPOCHS-"])
random_seed = int(values["-SEED-"])
if max_epochs < 1:
raise ValueError("Max epochs must be at least 1")
except ValueError:
sg.popup_error("Invalid Max Epochs or Seed value.")
continue
# Spectral & Welch Parameters
psd_method = values["-PSD_METHOD-"].lower()
welch_window_ms = None
welch_overlap = None
if psd_method == "welch":
try:
welch_window_ms = float(values["-WELCH_WINDOW-"])
welch_overlap = float(values["-WELCH_OVERLAP-"])
if welch_window_ms <= 0:
raise ValueError("Welch window length must be greater than 0")
if not 0 <= welch_overlap < 100:
raise ValueError("Welch overlap must be between 0 and 100")
except ValueError:
sg.popup_error("Invalid Welch parameters.")
continue
# Multitaper bandwidth (Hz). "Auto"/blank -> None (MNE default).
mt_bandwidth = None
if psd_method == "multitaper":
bw_raw = values["-MT_BANDWIDTH-"].strip()
if bw_raw and bw_raw.lower() != "auto":
try:
mt_bandwidth = float(bw_raw)
if mt_bandwidth <= 0:
raise ValueError("Bandwidth must be positive")
except ValueError:
sg.popup_error("Invalid multitaper bandwidth. Use a positive number or 'Auto'.")
continue
# Guardrail: check implied taper count against the entered epoch length.
if mt_bandwidth is not None and (values["-CALC_POWER-"] or values["-CALC_PEAK-"]
or save_psd_avg or save_psd_perchannel or save_psd_plots):
epoch_raw = values["-EPOCH_LENGTH-"].strip()
if epoch_raw:
try:
epoch_sec = float(epoch_raw)
implied_tapers = int(math.floor(epoch_sec * mt_bandwidth - 1))
if implied_tapers < 1:
if sg.popup_yes_no(
f"Bandwidth {mt_bandwidth} Hz with epoch length {epoch_sec} s "
f"yields {implied_tapers} tapers — multitaper will likely fail.\n\n"
"Proceed anyway?",
title="Bandwidth too small",
) != "Yes":
continue
elif implied_tapers < 3:
logger.warning(
f"Multitaper bandwidth {mt_bandwidth} Hz with epoch length "
f"{epoch_sec} s yields only {implied_tapers} taper(s)."
)
except ValueError:
pass # Non-critical; skip guardrail if epoch length is unparseable
power_fs = 256
try:
power_fs = float(values["-POWER_FS-"])
if power_fs <= 0:
raise ValueError("Sampling frequency must be greater than 0")
except ValueError:
sg.popup_error("Invalid sampling frequency value.")
continue
# Relative-power reference range (denominator for rel. power and SV)
total_power_range = (
FREQUENCY_BANDS["broadband"]["range"] if "broadband" in FREQUENCY_BANDS else (0.5, 47.0)
)
passband = None
reference_range_source = None
was_clamped = False
if values["-CALC_POWER-"] or values["-CALC_SV-"]:
# Prefer the pass-band the preprocessing script actually applied
# over the much weaker Nyquist bound.
passband = discover_preprocessing_passband(folder_path)
ceiling, ceiling_reason = effective_reference_ceiling(
power_fs, psd_method, passband
)
try:
total_min = float(values["-TOTAL_MIN-"])
total_max = float(values["-TOTAL_MAX-"])
if total_min < 0:
raise ValueError("Reference minimum must be >= 0 Hz.")
if total_min >= total_max:
raise ValueError("Reference minimum must be below the reference maximum.")
if total_max > power_fs / 2:
raise ValueError(
f"Reference maximum must not exceed Nyquist ({power_fs / 2:g} Hz)."
)
except ValueError as exc:
sg.popup_error(f"Invalid relative-power reference range.\n\n{exc}")
continue
if total_max > ceiling:
if total_min >= ceiling:
sg.popup_error(
"Invalid relative-power reference range.\n\n"
f"The reference minimum is {total_min:g} Hz, but the "
f"highest usable frequency is {ceiling:g} Hz "
f"because of the {ceiling_reason}.\n\n"
"Choose a lower reference minimum."
)
continue
detail = (
passband["detail"]
if passband is not None
else None
)
message = (
f"The reference range you entered ({total_min:g}-{total_max:g} Hz) "
f"extends beyond the {ceiling_reason}.\n\n"
)
if detail:
message += f"Detected from: {detail}.\n\n"
message += (
f"Above {ceiling:g} Hz the data contains only filter roll-off, so "
f"including it dilutes every relative power value.\n\n"
f"Clamp the reference range to {total_min:g}-{ceiling:g} Hz?"
)
choice = sg.popup_yes_no(message, title="Reference range exceeds usable band")
if choice == "Yes":
total_max = ceiling
was_clamped = True
window["-TOTAL_MAX-"].update(f"{total_max:g}")
logger.info(
f"Reference range clamped to {total_min:g}-{total_max:g} Hz "
f"({ceiling_reason})."
)
else:
logger.warning(
f"Reference range {total_min:g}-{total_max:g} Hz exceeds the "
f"{ceiling_reason}; relative power will include out-of-band content."
)
if passband is not None and total_min < passband["low"]:
logger.warning(
f"Reference minimum {total_min:g} Hz is below the preprocessing "
f"high-pass ({passband['low']:g} Hz); the bottom of the reference "
f"range contains attenuated signal only."
)
total_power_range = (total_min, total_max)
# Record where the bound came from, so the Excel metadata makes
# a clamped range traceable rather than silent.
if passband is None:
reference_range_source = (
"user-specified "
"(no unambiguous preprocessing pass-band found)"
)
else:
origin_label = {
"config": "batch config",
"sidecar": "fallback export metadata",
}.get(
passband["origin"],
str(passband["origin"]),
)
reference_range_source = (
f"pass-band {passband['low']:g}-{passband['high']:g} Hz "
f"from {origin_label}"
)
if was_clamped:
reference_range_source += "; range clamped to fit"
# Flag bands that stick out of the denominator.
outside = [
name for name, info in FREQUENCY_BANDS.items()
if name != "broadband"
and (info["range"][0] < total_min or info["range"][1] > total_max)
]
if outside:
logger.warning(
f"Bands partly outside the reference range {total_min:g}-{total_max:g} Hz: "
f"{', '.join(outside)}. Their relative power will be misleading."
)
logger.info(f"Relative power reference range: {total_min:g}-{total_max:g} Hz")
peak_min = peak_max = None
if values["-CALC_PEAK-"]:
try:
peak_min = float(values["-PEAK_MIN-"])
peak_max = float(values["-PEAK_MAX-"])
if peak_min >= peak_max:
raise ValueError("Minimum frequency must be less than maximum")
if peak_min < 0 or peak_max > (power_fs / 2):
raise ValueError(f"Frequency range must be between 0 and {power_fs / 2} Hz")
except ValueError:
sg.popup_error("Invalid peak frequency range.")
continue
# Connectivity Parameters (PLT, JPE)
jpe_st = 1
if values["-CALC_JPE-"]:
try:
jpe_st = int(values["-JPE_ST-"])
if jpe_st < 1:
raise ValueError("Time step must be greater than 0")
except ValueError:
sg.popup_error("Invalid JPE time step value.")
continue
calc_plt = values.get("-CALC_PLT-", False)
plt_threshold_ms = 30
if calc_plt:
try:
plt_threshold_ms = float(values["-PLT_THRESH-"])
if plt_threshold_ms < 0:
raise ValueError("PLT Threshold must be positive")
except ValueError:
sg.popup_error("Invalid PLT Threshold value.")
continue
# Complexity Parameters
sampen_m = 2
if values["-CALC_SAMPEN-"]:
try:
sampen_m = int(values["-SAMPEN_M-"])
if sampen_m < 1:
raise ValueError("Order m must be greater than 0")
except ValueError:
sg.popup_error("Invalid SampEn order parameter.")
continue
apen_m = 1
apen_r = 0.25
if values["-CALC_APEN-"]:
try:
apen_m = int(values["-APEN_M-"])
apen_r = float(values["-APEN_R-"])
if apen_m <= 0 or apen_r <= 0:
raise ValueError("Order m and tolerance r must be greater than 0")
except ValueError:
sg.popup_error("Invalid ApEn parameter.")
continue
sv_window = None
if values["-CALC_SV-"]:
try:
sv_window = int(values["-SV_WINDOW-"])
if sv_window < MIN_WINDOW_SIZE:
raise ValueError(f"Window length must be at least {MIN_WINDOW_SIZE} ms")
except ValueError:
sg.popup_error("Invalid spectral variability window.")
continue
# --- 3. Directory Parsing & Setup ---
grouped_files = group_epochs_by_condition(folder_path, folder_ext)
if not grouped_files:
continue # Grouping function handles its own error popups
# Create destination folders BEFORE processing so workers can save to disk
matrix_folders = None
if save_matrices or save_mst:
matrix_folders = create_matrix_folder_structure(
folder_path, matrix_folder, mst_folder if save_mst else None
)
logger.info("Initialized output directories for per-epoch matrix saving.")
psd_folder = None
if save_psd_any:
psd_folder = os.path.join(folder_path, psd_folder_name)
os.makedirs(psd_folder, exist_ok=True)
logger.info(f"PSD export directory ready: {psd_folder}")
def update_progress(value):
window["-PROGRESS-"].update(value)
window.refresh()
# --- 4. Parallel Execution ---
try:
# Capture both wide-format (results) and long-format (epoch_level_data) output
results, epoch_level_data = process_all_subjects(
grouped_files,
invert=values["-INVERT-"],
n_threads=n_threads,
calc_jpe=values["-CALC_JPE-"],
calc_pli=values["-CALC_PLI-"],
calc_pli_mst=values["-CALC_PLI_MST-"],
calc_aec=values["-CALC_AEC-"],
use_aecc=values["-USE_AECC-"],
force_positive=values["-AEC_FORCE_POSITIVE-"],
concat_aecc=values["-CONCAT_AECC-"],
has_headers=values["-HAS_HEADERS-"],
jpe_st=jpe_st,
calc_aec_mst=values["-CALC_AEC_MST-"],
calc_power=values["-CALC_POWER-"],
power_fs=power_fs,
calc_peak=values["-CALC_PEAK-"],
peak_min=peak_min,
peak_max=peak_max,
calc_sampen=values["-CALC_SAMPEN-"],
sampen_m=sampen_m,
calc_apen=values["-CALC_APEN-"],
apen_m=apen_m,
apen_r=apen_r,
calc_lzc=values["-CALC_LZC-"],
lzc_threshold=values["-LZC_THRESH-"],
calc_sv=values["-CALC_SV-"],
sv_window=sv_window,
save_matrices=save_matrices,
save_mst=save_mst,
save_channel_averages=values["-SAVE_CHANNEL_AVERAGES-"],
psd_method=psd_method,
welch_window_ms=welch_window_ms,
welch_overlap=welch_overlap,
mt_bandwidth=mt_bandwidth,
calc_plt=calc_plt,
plt_threshold_ms=plt_threshold_ms,
calc_plt_mst=values["-CALC_PLT_MST-"],
progress_callback=update_progress,
limit_epochs=limit_epochs,
max_epochs=max_epochs,
random_seed=random_seed,
matrix_folders=matrix_folders, # Routes paths to workers
save_epoch_metrics=save_epoch_metrics, # Tells workers to collect row data
save_psd_avg=save_psd_avg,
save_psd_perchannel=save_psd_perchannel,
save_psd_plots=save_psd_plots,
psd_folder=psd_folder,
total_power_range=total_power_range,
)
# --- 5. Export Results ---
if results:
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
output_path = os.path.join(folder_path, f"EEG_analysis_{timestamp}.xlsx")
try:
save_results_to_excel(
results,
output_path,
values["-INVERT-"],
values["-CALC_PLI_MST-"],
epoch_level_data=epoch_level_data, # Pass long-format data
calc_jpe=values["-CALC_JPE-"],
calc_pli=values["-CALC_PLI-"],
calc_aec=values["-CALC_AEC-"],
use_aecc=values["-USE_AECC-"],
force_positive=values["-AEC_FORCE_POSITIVE-"],
calc_aec_mst=values["-CALC_AEC_MST-"],
calc_power=values["-CALC_POWER-"],
power_fs=power_fs,
calc_peak=values["-CALC_PEAK-"],
peak_min=peak_min,
peak_max=peak_max,
calc_sampen=values["-CALC_SAMPEN-"],
calc_apen=values["-CALC_APEN-"],
calc_lzc=values["-CALC_LZC-"],
lzc_threshold=values["-LZC_THRESH-"],
calc_sv=values["-CALC_SV-"],
sv_window=sv_window if values["-CALC_SV-"] else None,
save_channel_averages=values["-SAVE_CHANNEL_AVERAGES-"],
concat_aecc=values["-CONCAT_AECC-"],
has_headers=values["-HAS_HEADERS-"],
psd_method=psd_method,
welch_window_ms=welch_window_ms,
welch_overlap=welch_overlap,
mt_bandwidth=mt_bandwidth,
calc_plt=calc_plt,
calc_plt_mst=values["-CALC_PLT_MST-"],
plt_threshold_ms=plt_threshold_ms,
save_psd_avg=save_psd_avg,
save_psd_perchannel=save_psd_perchannel,
save_psd_plots=save_psd_plots,
total_power_range=total_power_range,
reference_range_source=reference_range_source,
)
# Generate dynamic success summary
success_msg = f"Analysis complete!\nResults saved to:\n{output_path}"
if save_matrices:
success_msg += f"\n\nConnectivity matrices saved per-epoch in:\n{matrix_folder}"
if save_mst:
success_msg += f"\n\nMST matrices saved per-epoch in:\n{mst_folder}"
logger.info("Analysis run finished successfully.")
sg.popup(success_msg, title="Success")
except Exception:
logger.exception("Error during Excel export")
sg.popup_error("Analysis completed, but an error occurred while saving the Excel file. Check logs.")
else:
logger.warning("No results were generated")
sg.popup_error("Processing finished but no valid results were generated.")
except Exception:
logger.exception("Fatal error during parallel processing")
sg.popup_error("A fatal error occurred during processing. Please check the log file for details.")
finally:
logger.info("--- Run Terminated ---")
window.close()
if __name__ == "__main__":
main()

eeg_quantitative_analysis.py at commit f25ba96, under Apache-2.0 · at the source

Overview

Authors: D. Yorben Lodema¹, Herman J. van Dellen¹, Willem de Haan², Margot van Hest¹, Arjan Hillebrand^3,4,5, Edwin van Dellen^1,6

ORCID iDs: D. Yorben Lodema

Institutions: University Medical Center Utrecht (Netherlands); Amsterdam Neuroscience (Netherlands); Amsterdam University Medical Centers (Netherlands); Amsterdam UMC Location Vrije Universiteit Amsterdam (Netherlands); Vrije Universiteit Amsterdam (Netherlands); Vrije Universiteit Brussel (Belgium); Universitair Ziekenhuis Brussel (Belgium)

Journal: PLoS computational biology, volume 22, issue 3, article e1014043

Dates: received 30 October 2025; accepted 19 February 2026; published online 2 March 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1371/journal.pcbi.1014043 · PMID 41770742 · PMCID PMC12970966 · OpenAlex W7133208668

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), human (organism), methods / tools (subfield)

Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Complexity, Preprocessing, Connectivity, Source localization, Physiology & signal measures

MeSH: Electroencephalography*, Signal Processing, Computer-Assisted*, Software*, User-Computer Interface*, Algorithms, Brain, Computational Biology, Computer Graphics, Humans, Reproducibility of Results (* major topic)

Journal subjects: Software, Research and Analysis Methods, Bioassays and Physiological Analysis, Electrophysiological Techniques, Brain Electrophysiology, Electroencephalography, Biology and Life Sciences, Physiology, Electrophysiology, Neurophysiology, Neuroscience, Brain Mapping, Medicine and Health Sciences, Clinical Medicine, Clinical Neurophysiology, Imaging Techniques, Neuroimaging, Computer and Information Sciences, Software Engineering, Preprocessing, Engineering and Technology, Physical Sciences, Physics, Thermodynamics, Entropy, Computer Software, Open Source Software, Science Policy, Open Science, Mathematics, Discrete Mathematics, Combinatorics, Permutation, Signal Processing, Signal Filtering, Research Design, Quantitative Analysis, Electronics Engineering, Computer Engineering, Man-Computer Interface, Graphical User Interfaces, Computer Architecture, User Interfaces

Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Citations: not cited yet (Europe PMC); 23 references in the paper

Abstract

Processing of electroencephalography (EEG) data requires multiple steps to remove noise and artifacts and select good-quality data. While powerful open-source toolboxes like MNE-Python exist, their command-line nature can pose a barrier for researchers without programming experience. Here, we present EEG-Pype, an open-source (Apache-2.0 licensed) graphical user interface application using MNE-Python functions. EEG-Pype provides an intuitive workflow tailored for preprocessing of resting-state EEG data, including frequency band filtering, independent component analysis and atlas-based beamforming for source-level analysis. The application supports several common raw EEG input file formats and guides users through a comprehensive pipeline focused on manual bad channel and epoch selection. Manual steps are streamlined using MNE-Python’s interactive plots, resulting in a user-friendly experience. Configuration saving and loading allows for batch (re)runs, while a separate log is also saved, improving reproducibility and documentation. Output can be saved after filtering in canonical frequency bands, ready for further analysis. EEG-Pype includes a module for calculating quantitative EEG measures on preprocessed data, including spectral, functional connectivity and network analysis metrics. This software aims to lower the entry barrier for standardized EEG preprocessing, promoting reproducible research practices among neuroscientists and clinicians without requiring programming knowledge. EEG-Pype can be downloaded from: https://github.com/yorbenlodema/EEG-Pype and is not dependent on a specific operating system.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 8 matches between paragraphs and lines of code.

yorbenlodema/EEG-Pype

License: Apache-2.0

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Commit: f25ba96d55174d4a5d21c9016fed933bfd8873b6, 6 August 2026

Languages: Python (5)

Size: 29 files, 5 scripts

Software Heritage: not archived

Found in: “Availability and future directions”

Holds: README, license file, CITATION.cff, environment (Environment.yml, pyproject.toml)

Not found: tests, continuous integration, documentation

Tools: Matplotlib (4 files), MNE-Python (4 files), NumPy (3 files), pandas (3 files), SciPy (2 files), ICLabel (1 file), NetworkX (1 file), specparam (formerly FOOOF) (1 file)

Availability: 1 check, the latest on 30 September 2026: the link answers

30 September 2026: the link answers

7 files

src/__init__.py, Python, 1 line
src/eeg_aperiodic_analysis.py, Python, 3,485 lines
src/eeg_processing_script.py, Python, 5,097 lines, 3 matches
src/eeg_processing_settings.py, Python, 361 lines
src/eeg_quantitative_analysis.py, Python, 3,945 lines, 5 matches
LICENSE, License, 204 lines
README.md, Text, 469 lines

At the source: github.com/yorbenlodema/EEG-Pype

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
5 scripts, each with its path and the digest of its content;
8 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All relevant data are within the manuscript and its Supporting Information files, excluding our test.bdf EEG file due to file size constraints. This file is, and will stay, available on our GitHub page: https://github.com/yorbenlodema/EEG-Pype.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 10 MeSH terms, 20 references.

Cite

This paper

Lodema, D. Y., van Dellen, H. J., de Haan, W., van Hest, M., Hillebrand, A., & van Dellen, E. (2026). EEG-Pype: An accessible MNE-Python pipeline with graphical user interface for preprocessing and analysis of resting-state electroencephalography data. PLoS computational biology, 22(3), e1014043. https://doi.org/10.1371/journal.pcbi.1014043

BibTeX

@article{lodema2026eeg, author = {Lodema, D. Yorben and van Dellen, Herman J. and de Haan, Willem and van Hest, Margot and Hillebrand, Arjan and van Dellen, Edwin}, title = {{EEG-Pype: An accessible MNE-Python pipeline with graphical user interface for preprocessing and analysis of resting-state electroencephalography data}}, journal = {PLoS computational biology}, year = {2026}, month = mar, volume = {22}, number = {3}, pages = {e1014043}, publisher = {PLOS}, issn = {1553-734X}, doi = {10.1371/journal.pcbi.1014043}, url = {https://doi.org/10.1371/journal.pcbi.1014043}, pmid = {41770742}, pmcid = {PMC12970966} }

RIS

TY - JOUR AU - Lodema, D. Yorben AU - van Dellen, Herman J. AU - de Haan, Willem AU - van Hest, Margot AU - Hillebrand, Arjan AU - van Dellen, Edwin TI - EEG-Pype: An accessible MNE-Python pipeline with graphical user interface for preprocessing and analysis of resting-state electroencephalography data T2 - PLoS computational biology J2 - PLoS Comput Biol PY - 2026 DA - 2026/03/02 VL - 22 IS - 3 SP - e1014043 SN - 1553-734X PB - PLOS DO - 10.1371/journal.pcbi.1014043 UR - https://doi.org/10.1371/journal.pcbi.1014043 LA - en ER -

CSL-JSON

{ "id": "10.1371/journal.pcbi.1014043", "type": "article-journal", "title": "EEG-Pype: An accessible MNE-Python pipeline with graphical user interface for preprocessing and analysis of resting-state electroencephalography data", "container-title": "PLoS computational biology", "author": [ { "family": "Lodema", "given": "D. Yorben" }, { "family": "van Dellen", "given": "Herman J." }, { "family": "de Haan", "given": "Willem" }, { "family": "van Hest", "given": "Margot" }, { "family": "Hillebrand", "given": "Arjan" }, { "family": "van Dellen", "given": "Edwin" } ], "container-title-short": "PLoS Comput Biol", "volume": "22", "issue": "3", "page": "e1014043", "DOI": "10.1371/journal.pcbi.1014043", "PMID": "41770742", "PMCID": "PMC12970966", "ISSN": "1553-734X", "publisher": "PLOS", "URL": "https://doi.org/10.1371/journal.pcbi.1014043", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 2 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.3390/s26134019 [code]: NeuroStat: An Open-Source EEG Connectivity Platform for Randomised Controlled Trials.
Journal: Sensors (Basel, Switzerland)
In common: ICLabel, MNE-Python, NetworkX, 4 other tools, EEG, 8 references
[2] doi:10.1093/cercor/bhag113 [code]: Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 4 other tools, methods / tools, EEG, 2 references
[3] doi:10.1002/mds.70348 [code]: Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease.
Journal: Movement disorders : official journal of the Movement Disorder Society
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 4 other tools, EEG, 2 references
[4] doi:10.1162/imag.a.1269 [code]: From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 4 other tools, EEG, 2 references
[5] doi:10.3389/fncom.2026.1786996 [code]: Schumann-anchored golden ratio organization of human neural oscillations.
Journal: Frontiers in computational neuroscience
In common: specparam (formerly FOOOF), MNE-Python, NetworkX, 4 other tools, EEG, 2 references
[6] doi:10.1097/j.pain.0000000000004044 [code]: No effect of rhythmic visual stimulation on experimental pain perception.
Journal: Pain
In common: specparam (formerly FOOOF), ICLabel, MNE-Python, 4 other tools, EEG, 1 reference
[7] doi:10.1162/imag.a.1245 [code]: Towards precision EEG connectomics: Evaluating the benefits of dense sampling.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ICLabel, MNE-Python, pandas, 3 other tools, EEG, 3 references
[8] doi:10.1002/hbm.70628 [code]: EEG Biomarkers for Affective Disorders Diagnosis: An Evaluation and Validation Study.
Journal: Human brain mapping
In common: ICLabel, MNE-Python, pandas, 3 other tools, EEG, 2 references
[9] doi:10.1038/s42003-026-10394-7 [code]: Cognitive load weakens neural speech tracking without altering response timing.
Journal: Communications biology
In common: MNE-Python, NetworkX, pandas, 3 other tools, EEG, 3 references
[10] doi:10.7554/elife.100605 [code]: Age-related changes in ‘cortical’ 1/f dynamics are linked to cardiac activity
Journal: n/a
In common: specparam (formerly FOOOF), MNE-Python, NetworkX, 4 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.