From 4646440b54a5673fe6c0ce05f9cac8d19329c11b Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 15:39:20 -0700 Subject: [PATCH 1/9] Add MATLAB unittest for numdim --- tests/matlab/test_numdim.m | 52 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 tests/matlab/test_numdim.m diff --git a/tests/matlab/test_numdim.m b/tests/matlab/test_numdim.m new file mode 100644 index 00000000..72da3dc8 --- /dev/null +++ b/tests/matlab/test_numdim.m @@ -0,0 +1,52 @@ +function tests = test_numdim +% TEST_NUMDIM Tests for numdim.m (effective number of sources via +% eigenvalue-entropy of the channel second-order matrix A*A'/100). +% +% Expectations are closed-form / release-invariant: +% - single channel -> lambda == 1 +% - orthogonal equal-energy -> lambda == nchan (uniform eigenvalues) +% - identical channels (rank-1) -> lambda == 1 (one dominant eigenvalue) +% - full-rank -> real, 1 < lambda < nchan +tests = functiontests(localfunctions); +end + +function testSingleChannel(testCase) +% One channel -> one (normalized) eigenvalue == 1 -> entropy 0 -> lambda 1. +a = rand(1, 50) + 0.5; % nonzero energy +verifyEqual(testCase, numdim(a), 1, 'AbsTol', 1e-10); +end + +function testTwoChannelOrthogonal(testCase) +% A*A' = 2*I -> equal eigenvalues -> lambda == nchan == 2. +A = [1 1; 1 -1]; +verifyEqual(testCase, numdim(A), 2, 'AbsTol', 1e-10); +end + +function testOrthogonalEqualEnergy(testCase) +% Hadamard rows are orthogonal with equal norm -> A*A' = 4*I -> lambda == 4. +A = hadamard(4); +verifyEqual(testCase, numdim(A), 4, 'AbsTol', 1e-9); +end + +function testRankDeficientApproxOne(testCase) +% Identical channels -> rank-1 -> one dominant eigenvalue -> ~1 effective dim +% (eig yields tiny nonzero eigenvalues, so the result is finite, not NaN). +A = ones(3, 10); +verifyEqual(testCase, numdim(A), 1, 'AbsTol', 1e-6); +end + +function testFullRankBounds(testCase) +% Full-rank random data: real scalar strictly between 1 and nchan. +rng(42); +A = randn(5, 200); +v = numdim(A); +verifyTrue(testCase, isreal(v)); +verifyGreaterThan(testCase, v, 1); +verifyLessThan(testCase, v, 5); +end + +function testNoArgShowsHelp(testCase) +% nargin<1 branch: prints help and returns without error. +out = evalc('numdim()'); +verifyTrue(testCase, ~isempty(out)); +end From ae60496b9c30522e726fe232b9cf62a9b7f24d7d Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 15:48:36 -0700 Subject: [PATCH 2/9] Add numdim parity tests --- tests/test_numdim.py | 69 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 tests/test_numdim.py diff --git a/tests/test_numdim.py b/tests/test_numdim.py new file mode 100644 index 00000000..427f75b3 --- /dev/null +++ b/tests/test_numdim.py @@ -0,0 +1,69 @@ +"""Parity tests for numdim - effective number of sources (EEGLAB numdim.m). + +numdim estimates a lower bound on the number of discrete sources via the +eigenvalue entropy of the channel second-order matrix ``A @ A.T / 100``. The +measure is scale-invariant (eigenvalues are normalized), so expectations below +are closed-form / property-based. + +Expected values match the MATLAB reference tests in tests/matlab/test_numdim.m +(confirmed on MATLAB R2025a). Regenerate ground truth with +``matlab -batch "addpath('.../functions/miscfunc'); numdim()"``. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from eegprep.functions.miscfunc.numdim import numdim + +pytestmark = pytest.mark.parity + + +def test_single_channel_is_one(): + # One channel -> one normalized eigenvalue -> entropy 0 -> lambda == 1. + rng = np.random.default_rng(0) + a = rng.random((1, 50)) + 0.5 + np.testing.assert_allclose(numdim(a), 1.0, atol=1e-10) + + +def test_two_channel_orthogonal_is_nchan(): + # A @ A.T = 2*I -> equal eigenvalues -> lambda == nchan == 2. + a = np.array([[1.0, 1.0], [1.0, -1.0]]) + np.testing.assert_allclose(numdim(a), 2.0, atol=1e-10) + + +def test_hadamard_equal_energy_is_nchan(): + # Hadamard(4): orthogonal equal-norm rows -> A @ A.T = 4*I -> lambda == 4. + h4 = np.array( + [ + [1, 1, 1, 1], + [1, -1, 1, -1], + [1, 1, -1, -1], + [1, -1, -1, 1], + ], + dtype=float, + ) + np.testing.assert_allclose(numdim(h4), 4.0, atol=1e-9) + + +def test_rank_deficient_is_approx_one(): + # Identical channels -> rank-1 -> one dominant eigenvalue -> ~1 effective dim. + a = np.ones((3, 10)) + np.testing.assert_allclose(numdim(a), 1.0, atol=1e-6) + + +def test_full_rank_matches_matlab(): + # Full-rank deterministic integer matrix: identical in numpy and MATLAB, so + # pin to the exact MATLAB R2025a value (confirmed via matlab -batch). The + # result is real and strictly between 1 and nchan (=3). + a = np.array( + [ + [2.0, 1.0, 0.0, 1.0], + [1.0, 3.0, 1.0, 0.0], + [0.0, 1.0, 2.0, 1.0], + ] + ) + v = numdim(a) + assert np.isreal(v) + np.testing.assert_allclose(v, 2.147746217856, rtol=1e-6) From 624e14551b61b738d0038cd8c9923d696fc9a0fc Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 15:52:59 -0700 Subject: [PATCH 3/9] Add numdim Python implementation --- src/eegprep/functions/miscfunc/numdim.py | 34 ++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/eegprep/functions/miscfunc/numdim.py diff --git a/src/eegprep/functions/miscfunc/numdim.py b/src/eegprep/functions/miscfunc/numdim.py new file mode 100644 index 00000000..029b167b --- /dev/null +++ b/src/eegprep/functions/miscfunc/numdim.py @@ -0,0 +1,34 @@ +"""EEGLAB-compatible estimate of the effective number of sources (numdim).""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + + +def numdim(data: Any) -> float: + """Estimate a lower bound on the number of discrete sources in the data. + + Ports EEGLAB's ``numdim`` (Wackermann's measure): the effective + dimensionality is the exponential of the Shannon entropy of the normalized + eigenvalue spectrum of the channel second-order matrix ``data @ data.T / 100``. + The measure is invariant to overall scaling of the data. + + Args: + data: 2-D array shaped ``(channels, samples)``. + + Returns: + The estimated effective number of sources, a real scalar. + """ + a = np.asarray(data, dtype=float).T # MATLAB: a = a'; + b = a.T @ a / 100.0 # MATLAB: b = a'*a/100; + eigenvalues = np.linalg.eigvals(b) # MATLAB: [v d] = eig(b); + weights = (eigenvalues / np.sum(eigenvalues)).astype(complex) + # Complex log mirrors MATLAB's real(exp(...)): tiny/negative eigenvalues + # from finite precision contribute ~0 instead of producing NaN. + lambda_ = np.exp(-np.sum(weights * np.log(weights))) + return float(np.real(lambda_)) + + +__all__ = ["numdim"] From 8b037eca9f53a3c39b896dac84dfe5edead20076 Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 15:58:15 -0700 Subject: [PATCH 4/9] Add MATLAB unittest for eeg_ms2f --- tests/matlab/test_eeg_ms2f.m | 66 ++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 tests/matlab/test_eeg_ms2f.m diff --git a/tests/matlab/test_eeg_ms2f.m b/tests/matlab/test_eeg_ms2f.m new file mode 100644 index 00000000..db08756a --- /dev/null +++ b/tests/matlab/test_eeg_ms2f.m @@ -0,0 +1,66 @@ +function tests = test_eeg_ms2f +% TEST_EEG_MS2F Tests for eeg_ms2f.m (epoch latency in ms -> nearest 1-based +% epoch frame). outf = 1 + round((pnts-1)*(ms/1000 - xmin)/(xmax - xmin)). +% Closed-form, release-invariant expectations. +tests = functiontests(localfunctions); +end + +function eeg = makeEEG(xmin, xmax, pnts) +eeg = struct('xmin', xmin, 'xmax', xmax, 'pnts', pnts); +end + +function testFirstFrameAtXmin(testCase) +% ms at xmin -> first frame (1-based). +EEG = makeEEG(0, 1, 1001); +verifyEqual(testCase, eeg_ms2f(EEG, 0), 1); +end + +function testLastFrameAtXmax(testCase) +% ms at xmax -> last frame == pnts. +EEG = makeEEG(0, 1, 1001); +verifyEqual(testCase, eeg_ms2f(EEG, 1000), 1001); +end + +function testMidpoint(testCase) +% 500 ms -> 0.5 s -> frame 1 + 1000*0.5 = 501. +EEG = makeEEG(0, 1, 1001); +verifyEqual(testCase, eeg_ms2f(EEG, 500), 501); +end + +function testRoundsToNearest(testCase) +% 499.4 ms -> 1 + round(499.4) = 500. +EEG = makeEEG(0, 1, 1001); +verifyEqual(testCase, eeg_ms2f(EEG, 499.4), 500); +end + +function testEpochCenterNegativeXmin(testCase) +% Epoched data xmin<0: ms=0 -> centre frame. +EEG = makeEEG(-1, 1, 2001); +verifyEqual(testCase, eeg_ms2f(EEG, 0), 1001); +end + +function testBelowRangeErrors(testCase) +% ms below xmin -> error('time out of range'). +EEG = makeEEG(0, 1, 1001); +threw = false; +try + eeg_ms2f(EEG, -1); +catch ME + threw = true; + verifyTrue(testCase, contains(ME.message, 'out of range')); +end +verifyTrue(testCase, threw); +end + +function testAboveRangeErrors(testCase) +% ms above xmax -> error('time out of range'). +EEG = makeEEG(0, 1, 1001); +threw = false; +try + eeg_ms2f(EEG, 2000); +catch ME + threw = true; + verifyTrue(testCase, contains(ME.message, 'out of range')); +end +verifyTrue(testCase, threw); +end From e3ce99709898a9ae93096dc1f59f376b47ee0751 Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 16:00:35 -0700 Subject: [PATCH 5/9] Add eeg_ms2f parity tests --- tests/test_eeg_ms2f.py | 53 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 tests/test_eeg_ms2f.py diff --git a/tests/test_eeg_ms2f.py b/tests/test_eeg_ms2f.py new file mode 100644 index 00000000..ce7eaa12 --- /dev/null +++ b/tests/test_eeg_ms2f.py @@ -0,0 +1,53 @@ +"""Parity tests for eeg_ms2f - epoch latency (ms) to nearest epoch frame. + +Ports EEGLAB eeg_ms2f.m: ``outf = 1 + round((pnts-1)*(ms/1000 - xmin)/(xmax - xmin))`` +with an out-of-range error. Frame numbers are 1-based to match EEGLAB. Expected +values are closed-form and mirror tests/matlab/test_eeg_ms2f.m. +""" + +from __future__ import annotations + +import pytest + +from eegprep.functions.miscfunc.eeg_ms2f import eeg_ms2f + +pytestmark = pytest.mark.parity + + +def _eeg(xmin, xmax, pnts): + return {"xmin": xmin, "xmax": xmax, "pnts": pnts} + + +def test_first_frame_at_xmin(): + # ms at xmin -> first frame (1-based). + assert eeg_ms2f(_eeg(0, 1, 1001), 0) == 1 + + +def test_last_frame_at_xmax(): + # ms at xmax -> last frame == pnts. + assert eeg_ms2f(_eeg(0, 1, 1001), 1000) == 1001 + + +def test_midpoint(): + # 500 ms -> 0.5 s -> frame 1 + 1000*0.5 = 501. + assert eeg_ms2f(_eeg(0, 1, 1001), 500) == 501 + + +def test_rounds_to_nearest(): + # 499.4 ms -> 1 + round(499.4) = 500. + assert eeg_ms2f(_eeg(0, 1, 1001), 499.4) == 500 + + +def test_epoch_center_negative_xmin(): + # Epoched data xmin<0: ms=0 -> centre frame. + assert eeg_ms2f(_eeg(-1, 1, 2001), 0) == 1001 + + +def test_below_range_raises(): + with pytest.raises(ValueError, match="out of range"): + eeg_ms2f(_eeg(0, 1, 1001), -1) + + +def test_above_range_raises(): + with pytest.raises(ValueError, match="out of range"): + eeg_ms2f(_eeg(0, 1, 1001), 2000) From 1675f25b3ee5f49a568bbf063c53fde744965fe6 Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 16:03:03 -0700 Subject: [PATCH 6/9] Add eeg_ms2f Python implementation --- src/eegprep/functions/miscfunc/eeg_ms2f.py | 40 ++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 src/eegprep/functions/miscfunc/eeg_ms2f.py diff --git a/src/eegprep/functions/miscfunc/eeg_ms2f.py b/src/eegprep/functions/miscfunc/eeg_ms2f.py new file mode 100644 index 00000000..8b488666 --- /dev/null +++ b/src/eegprep/functions/miscfunc/eeg_ms2f.py @@ -0,0 +1,40 @@ +"""EEGLAB-compatible conversion of epoch latency (ms) to nearest epoch frame.""" + +from __future__ import annotations + +import math +from typing import Any + + +def eeg_ms2f(EEG: dict[str, Any], ms: float) -> int: + """Convert an epoch latency in milliseconds to the nearest epoch frame. + + Ports EEGLAB's eeg_ms2f.m. The returned frame number is **1-based** to match + EEGLAB; convert to a 0-based index at the point of use. + + Args: + EEG: EEG structure with ``xmin``/``xmax`` (epoch limits, seconds) and + ``pnts`` (samples per epoch). + ms: Epoch latency in milliseconds. + + Returns: + The nearest 1-based epoch frame number. + + Raises: + ValueError: If the latency falls outside ``[xmin, xmax]``. + """ + seconds = ms / 1000.0 + xmin = EEG["xmin"] + xmax = EEG["xmax"] + if seconds < xmin or seconds > xmax: + raise ValueError("time out of range") + frac = (EEG["pnts"] - 1) * (seconds - xmin) / (xmax - xmin) + return 1 + _round_half_away_from_zero(frac) + + +def _round_half_away_from_zero(value: float) -> int: + """Round to the nearest integer with ties going away from zero (MATLAB ``round``).""" + return math.floor(value + 0.5) if value >= 0 else math.ceil(value - 0.5) + + +__all__ = ["eeg_ms2f"] From fa25f6a97b902a9ea3b8e43fdfac8ac772d54b7b Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 16:05:41 -0700 Subject: [PATCH 7/9] Add MATLAB unittest for spher --- tests/matlab/test_spher.m | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 tests/matlab/test_spher.m diff --git a/tests/matlab/test_spher.m b/tests/matlab/test_spher.m new file mode 100644 index 00000000..4048be58 --- /dev/null +++ b/tests/matlab/test_spher.m @@ -0,0 +1,34 @@ +function tests = test_spher +% TEST_SPHER Tests for spher.m (sphering matrix: 2*inv(sqrtm(cov(data')))). +% Properties: symmetric; whitens the channel covariance so that +% S*C*S' = 4*I (since the factor is 2). Single channel has a closed form. +tests = functiontests(localfunctions); +end + +function testSingleChannelClosedForm(testCase) +% 1 channel: cov(data') = var(data) (N-1 normalized) = 2.5 -> 2/sqrt(2.5). +data = [1 2 3 4 5]; +expected = 2 / sqrt(2.5); +verifyEqual(testCase, spher(data), expected, 'RelTol', 1e-10); +end + +function testSymmetric(testCase) +% Sphering matrix of a symmetric PSD covariance is symmetric. +data = [2 1 0 1 3 2; 0 2 1 3 1 0; 1 0 3 1 2 1]; +S = spher(data); +verifyLessThan(testCase, norm(S - S.'), 1e-9); +end + +function testWhitensCovarianceToFourI(testCase) +% S = 2*inv(sqrtm(C)) -> S*C*S' = 4*I. +data = [2 1 0 1 3 2; 0 2 1 3 1 0; 1 0 3 1 2 1]; +S = spher(data); +C = cov(data'); +verifyLessThan(testCase, norm(S * C * S.' - 4 * eye(3)), 1e-6); +end + +function testNoArgShowsHelp(testCase) +% nargin<1 branch: prints help and returns without error. +out = evalc('spher()'); +verifyTrue(testCase, ~isempty(out)); +end From e3727f1115a21c3eb3824633e4afddc36d8b2bad Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 16:07:45 -0700 Subject: [PATCH 8/9] Add spher parity tests --- tests/test_spher.py | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 tests/test_spher.py diff --git a/tests/test_spher.py b/tests/test_spher.py new file mode 100644 index 00000000..4e3219e0 --- /dev/null +++ b/tests/test_spher.py @@ -0,0 +1,44 @@ +"""Parity tests for spher - sphering matrix (EEGLAB spher.m). + +``spher(data) = 2 * inv(sqrtm(cov(data.T)))``. The result is symmetric and +whitens the channel covariance so that ``S @ C @ S.T == 4*I`` (the factor is 2). +Single-channel data has the closed form ``2 / sqrt(var(data))``. Expectations +mirror tests/matlab/test_spher.m (confirmed on MATLAB R2025a). +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from eegprep.functions.sigprocfunc.spher import spher + +pytestmark = pytest.mark.parity + +# Deterministic full-rank 3-channel x 6-sample data (matches the MATLAB test). +DATA = np.array( + [ + [2.0, 1.0, 0.0, 1.0, 3.0, 2.0], + [0.0, 2.0, 1.0, 3.0, 1.0, 0.0], + [1.0, 0.0, 3.0, 1.0, 2.0, 1.0], + ] +) + + +def test_single_channel_closed_form(): + # 1 channel: cov(data.T) = var(data) (N-1) = 2.5 -> 2/sqrt(2.5). + data = np.array([[1.0, 2.0, 3.0, 4.0, 5.0]]) + np.testing.assert_allclose(np.squeeze(spher(data)), 2 / np.sqrt(2.5), rtol=1e-10) + + +def test_symmetric(): + # Sphering matrix of a symmetric PSD covariance is symmetric. + s = spher(DATA) + np.testing.assert_allclose(s, s.T, atol=1e-9) + + +def test_whitens_covariance_to_four_i(): + # S = 2*inv(sqrtm(C)) -> S @ C @ S.T = 4*I. + s = spher(DATA) + c = np.cov(DATA) # rowvar=True, ddof=1 == MATLAB cov(data.T) + np.testing.assert_allclose(s @ c @ s.T, 4 * np.eye(3), atol=1e-6) From 52b430912d44bbf65c99366bf8ef98aa7990e1be Mon Sep 17 00:00:00 2001 From: innaamogolonova Date: Sat, 27 Jun 2026 16:09:01 -0700 Subject: [PATCH 9/9] Add spher Python implementation --- src/eegprep/functions/sigprocfunc/spher.py | 34 ++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/eegprep/functions/sigprocfunc/spher.py diff --git a/src/eegprep/functions/sigprocfunc/spher.py b/src/eegprep/functions/sigprocfunc/spher.py new file mode 100644 index 00000000..676ebe9d --- /dev/null +++ b/src/eegprep/functions/sigprocfunc/spher.py @@ -0,0 +1,34 @@ +"""EEGLAB-compatible sphering (whitening) matrix (spher).""" + +from __future__ import annotations + +from typing import Any + +import numpy as np +from scipy.linalg import sqrtm + + +def spher(data: Any) -> np.ndarray: + """Return the sphering (whitening) matrix for the given data. + + Ports EEGLAB's spher.m: ``2 * inv(sqrtm(cov(data.T)))``. The result is a + symmetric matrix ``S`` that whitens the channel covariance ``C = cov(data.T)`` + so that ``S @ C @ S.T == 4 * I``. + + Args: + data: 2-D array shaped ``(channels, samples)``. + + Returns: + The ``(channels, channels)`` sphering matrix. + """ + array = np.asarray(data, dtype=float) + # np.cov (rowvar=True, ddof=1) matches MATLAB cov(data'); atleast_2d keeps + # the single-channel case as a 1x1 matrix. + covariance = np.atleast_2d(np.cov(array)) + sphere = 2.0 * np.linalg.inv(sqrtm(covariance)) + # sqrtm returns a complex dtype with ~0 imaginary part for real PSD input; + # MATLAB's result is real, so drop the negligible imaginary component. + return np.real_if_close(sphere) + + +__all__ = ["spher"]