diff --git a/meegkit/asr.py b/meegkit/asr.py index 58ae53b5..05914263 100755 --- a/meegkit/asr.py +++ b/meegkit/asr.py @@ -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): @@ -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) @@ -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 " diff --git a/tests/test_asr.py b/tests/test_asr.py index 3aa9304c..842e6586 100644 --- a/tests/test_asr.py +++ b/tests/test_asr.py @@ -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 @@ -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).