concord-metrics-spectral

Frequency-domain metrics: Welch PSD and band power

Role in the system Implements spectral analysis metrics. Both classes use scipy.signal.welch internally and return MetricResult containers with freq_axis populated.

Package Structure

concord-metrics-spectral/
  src/concord_metrics_spectral/
    __init__.py     exports WelchPSD, BandPower
    welch.py        WelchPSD metric
    band_power.py   BandPower metric

WelchPSD

Computes the Power Spectral Density (PSD) of each channel using Welch's method. Welch's method reduces variance by averaging the periodograms of overlapping segments.

What is PSD? PSD tells you how much power the signal contains at each frequency. The x-axis is frequency (Hz), the y-axis is power density (V²/Hz). A peak at 10 Hz means the signal has strong alpha-band activity.

classWelchPSD(window_s=4.0, overlap=0.5, fmin=0.0, fmax=None, window="hann")

Parameter	Default	Description
window_s	4.0	Length of each FFT window in seconds. Longer = better frequency resolution.
overlap	0.5	Fractional overlap between FFT windows (0–1). 0.5 = 50% overlap.
fmin	0.0	Minimum frequency to return (Hz).
fmax	None	Maximum frequency to return (Hz). If None, returns up to Nyquist (fs/2).
window	"hann"	Window function name (passed to scipy). "hann" is the standard choice for spectral analysis.

name: "welch_psd"

units: "V^2/Hz"

Requires: min_duration_s = window_s

Output

Field	Value
data.shape	(n_channels, n_freqs)
freq_axis	Array of frequency values in Hz, length n_freqs
time_axis	None (PSD is a global estimate over the whole recording)

Frequency resolution

Frequency resolution = 1 / window_s. With the default 4-second window, resolution is 0.25 Hz. This means frequencies are reported at 0, 0.25, 0.5, ... Hz.

How it works internally

import scipy.signal as ss

# For each channel:
freqs, psd = ss.welch(
    data[ch],
    fs=fs,
    window=window,
    nperseg=int(window_s * fs),
    noverlap=int(window_s * fs * overlap),
)
# Then mask to [fmin, fmax]

BandPower

Computes the total power in each of the standard EEG frequency bands for every channel. Internally runs WelchPSD then integrates the PSD over each band using the trapezoidal rule.

classBandPower(bands=None, window_s=4.0, overlap=0.5)

Parameter	Default	Description
bands	None	Dict mapping band name → (fmin, fmax). If None, uses standard EEG bands (see below).
window_s	4.0	Passed to WelchPSD.
overlap	0.5	Passed to WelchPSD.

name: "band_power"

units: "V^2"

Default frequency bands

Band	Range	Clinical significance
delta	0.5–4 Hz	Deep sleep, pathological slow waves, post-ictal suppression
theta	4–8 Hz	Drowsiness, hippocampal activity, memory encoding
alpha	8–13 Hz	Relaxed wakefulness, posterior dominant rhythm
beta	13–30 Hz	Active cognition, motor activity, anxiolytic drug effects
gamma	30–80 Hz	High-frequency oscillations, active processing, seizure fast ripples
high_gamma	80–150 Hz	High-frequency oscillations (HFOs), very-high-frequency seizure activity

Output

Field	Value
data.shape	(n_channels, n_bands)
freq_axis	None (bands summarize ranges)
metadata["bands"]	List of band names in axis-1 order, e.g. ["delta", "theta", "alpha", "beta", "gamma", "high_gamma"]

Example Usage

from concord_io import read_bids_ieeg
from concord_metrics_spectral import WelchPSD, BandPower

rec = read_bids_ieeg("sub-HUP117_run-01_ieeg.edf")

# PSD
psd = WelchPSD(window_s=4.0, fmin=1.0, fmax=150.0).compute(rec)
print(psd.data.shape)        # (72, n_freqs)
print(psd.freq_axis[:5])     # [1.0, 1.25, 1.5, 1.75, 2.0]  @ 4s window

# Band power
bp = BandPower().compute(rec)
print(bp.data.shape)         # (72, 6)
print(bp.metadata["bands"])  # ["delta", "theta", "alpha", "beta", "gamma", "high_gamma"]
gamma_power = bp.data[:, 4] # gamma band for all channels