Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 12 additions & 3 deletions meegkit/asr.py
Original file line number Diff line number Diff line change
Expand Up @@ -288,6 +288,15 @@ def _recency_weights(counter, sample_weight):
return weights


def _round_half_away(x):
"""Round half away from zero (MATLAB `round` semantics).

numpy's ``np.round`` rounds half to even (banker's rounding); this
matches MATLAB/Octave which round halves away from zero.
"""
return np.floor(np.abs(x) + 0.5) * np.sign(x)


def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5],
win_len=.5, win_overlap=0.66, min_clean_fraction=0.25,
max_dropout_fraction=0.1, show=False):
Expand Down Expand Up @@ -354,11 +363,11 @@ def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5],
truncate_quant = [0.0220, 0.6000]
step_sizes = [0.01, 0.01]
shape_range = SHAPE_RANGE
max_bad_chans = np.round(X.shape[0] * max_bad_chans)
max_bad_chans = _round_half_away(X.shape[0] * max_bad_chans)

# set data indices
[nc, ns] = X.shape
N = int(win_len * sfreq)
N = int(_round_half_away(win_len * sfreq))
N_raw = win_len * sfreq # non-truncated N, avoids step-size phase drift
offsets = np.round(np.arange(0, ns - N, N_raw * (1 - win_overlap)))
offsets = offsets.astype(int)
Expand Down Expand Up @@ -543,7 +552,7 @@ def asr_calibrate(X, sfreq, cutoff=5, blocksize=100, win_len=0.5,
X, _zf = yulewalk_filter(X, sfreq, ab=None)

# window length for calculating thresholds
N = int(np.round(win_len * sfreq))
N = int(_round_half_away(win_len * sfreq))
if ns < N:
raise ValueError(
f"Calibration data has {ns} samples, shorter than one analysis "
Expand Down
30 changes: 29 additions & 1 deletion tests/test_asr.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
import pytest
from scipy import signal

from meegkit.asr import ASR, asr_calibrate, asr_process, clean_windows
from meegkit.asr import ASR, _round_half_away, asr_calibrate, asr_process, clean_windows
from meegkit.utils.asr import SHAPE_RANGE, fit_eeg_distribution, yulewalk, yulewalk_filter
from meegkit.utils.covariances import block_covariance
from meegkit.utils.matrix import sliding_window
Expand Down Expand Up @@ -142,6 +142,34 @@ def test_yulewalk_filter(n_chans, show=False):
plt.show()


def test_round_half_away():
"""_round_half_away rounds ties away from zero, unlike numpy's banker's rounding."""
# exact-half (tie) behavior: rounds away from zero
assert _round_half_away(2.5) == 3
assert _round_half_away(128.5) == 129
assert _round_half_away(-2.5) == -3
assert _round_half_away(0.5) == 1
assert _round_half_away(1.5) == 2 # np.round gives 2 here too (even)

# divergence from banker's rounding on ties
assert _round_half_away(2.5) != np.round(2.5) # 3 != 2
assert _round_half_away(128.5) != np.round(128.5) # 129 != 128

# non-tie values round normally
assert _round_half_away(2.4) == 2
assert _round_half_away(2.6) == 3
assert _round_half_away(-2.6) == -3

# int() of the result is exact for scalars
assert int(_round_half_away(2.5)) == 3
assert int(_round_half_away(-2.5)) == -3

# array input
ties = np.array([2.5, 128.5, -2.5, 0.5, 1.5])
expected = np.array([3, 129, -3, 1, 2])
np.testing.assert_array_equal(_round_half_away(ties), expected)


def test_asr_functions(show=False, method="riemann"):
"""Test ASR functions (offline use).

Expand Down
Loading