feat(eeg): add compute_features_from_epoch (PSD bands + 5 statistics)

src/pipelines/eeg_pipeline.py

@@ -15,6 +15,8 @@ from __future__ import annotations
 import mne
 import numpy as np
 from mne.preprocessing import ICA
+from scipy import signal as scipy_signal
+from scipy import stats as scipy_stats
 
 from src.core.logger import get_logger
 
@@ -159,3 +161,57 @@ def remove_artifacts_with_ica(
     )
     ica.apply(out, verbose="ERROR")
     return out
+
+
+EEG_BANDS: dict[str, tuple[float, float]] = {
+    "delta": (1.0, 4.0),
+    "theta": (4.0, 8.0),
+    "alpha": (8.0, 13.0),
+    "beta":  (13.0, 30.0),
+    "gamma": (30.0, 40.0),
+}
+STATS: tuple[str, ...] = ("mean", "std", "var", "skew", "kurtosis")
+
+
+def _band_power(freqs: np.ndarray, psd: np.ndarray, lo: float, hi: float) -> float:
+    """Mean PSD value within the [lo, hi) frequency band."""
+    mask = (freqs >= lo) & (freqs < hi)
+    if not mask.any():
+        return 0.0
+    return float(psd[mask].mean())
+
+
+def compute_features_from_epoch(epoch: np.ndarray, sfreq: float) -> np.ndarray:
+    """Compute PSD-band + statistical features for one epoch.
+
+    Per channel, the feature block is:
+      [psd_delta, psd_theta, psd_alpha, psd_beta, psd_gamma,
+       mean, std, var, skew, kurtosis]
+    Channels are stacked in their input order. The resulting 1-D vector has
+    length ``n_channels * (len(EEG_BANDS) + len(STATS))``.
+
+    PSD is computed with Welch's method (`scipy.signal.welch`) at the
+    epoch's sample rate. Higher moments use `scipy.stats` with default
+    bias correction.
+
+    Args:
+        epoch: A 2-D array shape (n_channels, n_samples).
+        sfreq: Sampling rate in Hz.
+
+    Returns:
+        A 1-D `np.ndarray` of dtype float64.
+    """
+    n_channels, n_samples = epoch.shape
+    nperseg = min(256, n_samples)
+    feats: list[float] = []
+    for ch in range(n_channels):
+        x = epoch[ch]
+        freqs, psd = scipy_signal.welch(x, fs=sfreq, nperseg=nperseg)
+        for _band, (lo, hi) in EEG_BANDS.items():
+            feats.append(_band_power(freqs, psd, lo, hi))
+        feats.append(float(np.mean(x)))
+        feats.append(float(np.std(x)))
+        feats.append(float(np.var(x)))
+        feats.append(float(scipy_stats.skew(x)))
+        feats.append(float(scipy_stats.kurtosis(x)))
+    return np.asarray(feats, dtype=np.float64)

tests/pipelines/test_eeg_pipeline.py

@@ -9,6 +9,7 @@ import pytest
 
 from src.pipelines.eeg_pipeline import (
     bandpass_filter,
+    compute_features_from_epoch,
     is_valid_epoch,
     remove_artifacts_with_ica,
 )
@@ -17,6 +18,10 @@ from src.pipelines.eeg_pipeline import (
 FIXTURE = Path(__file__).parent.parent / "fixtures" / "eeg_sample.fif"
 
 
+EEG_BANDS = ("delta", "theta", "alpha", "beta", "gamma")
+STATS = ("mean", "std", "var", "skew", "kurtosis")
+
+
 class TestIsValidEpoch:
     def test_accepts_2d_finite_array(self) -> None:
         epoch = np.zeros((4, 256), dtype=np.float64)
@@ -175,3 +180,42 @@ class TestRemoveArtifactsWithIca:
         np.testing.assert_allclose(out.get_data(), raw.get_data(), rtol=1e-6, atol=1e-12)
         log_output = buf.getvalue()
         assert "ICA skipped: eog_ch_name='EOG_DOES_NOT_EXIST' not found" in log_output
+
+
+class TestComputeFeaturesFromEpoch:
+    def test_returns_1d_float_array(self) -> None:
+        epoch = np.random.default_rng(0).standard_normal((4, 256))
+        out = compute_features_from_epoch(epoch, sfreq=256.0)
+        assert isinstance(out, np.ndarray)
+        assert out.ndim == 1
+        assert out.dtype == np.float64
+
+    def test_feature_count_matches_contract(self) -> None:
+        """Each channel contributes len(EEG_BANDS) PSD features + len(STATS) stats."""
+        n_channels = 4
+        epoch = np.random.default_rng(0).standard_normal((n_channels, 256))
+        out = compute_features_from_epoch(epoch, sfreq=256.0)
+        expected = n_channels * (len(EEG_BANDS) + len(STATS))
+        assert out.shape == (expected,)
+
+    def test_alpha_band_dominates_for_alpha_signal(self) -> None:
+        """Pure 10 Hz sine on 1 channel should put most PSD power in alpha (8-13 Hz)."""
+        sfreq = 256.0
+        t = np.arange(int(sfreq * 2.0)) / sfreq
+        signal = np.sin(2 * np.pi * 10.0 * t)[None, :]  # (1, n_samples)
+        out = compute_features_from_epoch(signal, sfreq=sfreq)
+        # Layout for n_channels=1: [psd_delta, psd_theta, psd_alpha, psd_beta, psd_gamma, mean, std, var, skew, kurtosis]
+        psd_block = out[: len(EEG_BANDS)]
+        alpha_idx = EEG_BANDS.index("alpha")
+        assert psd_block[alpha_idx] == psd_block.max()
+
+    def test_finite_output(self) -> None:
+        epoch = np.random.default_rng(0).standard_normal((4, 256))
+        out = compute_features_from_epoch(epoch, sfreq=256.0)
+        assert np.all(np.isfinite(out))
+
+    def test_deterministic_for_same_input(self) -> None:
+        epoch = np.random.default_rng(0).standard_normal((4, 256))
+        a = compute_features_from_epoch(epoch, sfreq=256.0)
+        b = compute_features_from_epoch(epoch, sfreq=256.0)
+        np.testing.assert_array_equal(a, b)