fix sfreq estimation for snirf files (#13184)

Author: drammock

doc/changes/devel/13184.bugfix.rst

@@ -0,0 +1 @@
+Fix bug with sampling frequency estimation in snirf files, by `Daniel McCloy`_.

doc/changes/devel/13184.newfeature.rst

@@ -0,0 +1 @@
+New argument ``sfreq`` to :func:`mne.io.read_raw_snirf`, to allow overriding the sampling frequency estimated from (possibly jittered) sampling periods in the file, by `Daniel McCloy`_.

mne/io/snirf/_snirf.py

@@ -14,14 +14,22 @@
 from ..._freesurfer import get_mni_fiducials
 from ...annotations import Annotations
 from ...transforms import _frame_to_str, apply_trans
-from ...utils import _check_fname, _import_h5py, fill_doc, logger, verbose, warn
+from ...utils import (
+    _check_fname,
+    _import_h5py,
+    _validate_type,
+    fill_doc,
+    logger,
+    verbose,
+    warn,
+)
 from ..base import BaseRaw
 from ..nirx.nirx import _convert_fnirs_to_head
 
 
 @fill_doc
 def read_raw_snirf(
-    fname, optode_frame="unknown", preload=False, verbose=None
+    fname, optode_frame="unknown", *, sfreq=None, preload=False, verbose=None
 ) -> "RawSNIRF":
     """Reader for a continuous wave SNIRF data.
 
@@ -41,6 +49,11 @@ def read_raw_snirf(
         in which case the positions are not modified. If a known coordinate
         frame is provided (head, meg, mri), then the positions are transformed
         in to the Neuromag head coordinate frame (head).
+    sfreq : float | None
+        The nominal sampling frequency at which the data were acquired. If ``None``,
+        will be estimated from the time data in the file.
+
+        .. versionadded:: 1.10
     %(preload)s
     %(verbose)s
 
@@ -54,7 +67,7 @@ def read_raw_snirf(
     --------
     mne.io.Raw : Documentation of attributes and methods of RawSNIRF.
     """
-    return RawSNIRF(fname, optode_frame, preload, verbose)
+    return RawSNIRF(fname, optode_frame, sfreq=sfreq, preload=preload, verbose=verbose)
 
 
 def _open(fname):
@@ -74,6 +87,11 @@ class RawSNIRF(BaseRaw):
         in which case the positions are not modified. If a known coordinate
         frame is provided (head, meg, mri), then the positions are transformed
         in to the Neuromag head coordinate frame (head).
+    sfreq : float | None
+        The nominal sampling frequency at which the data were acquired. If ``None``,
+        will be estimated from the time data in the file.
+
+        .. versionadded:: 1.10
     %(preload)s
     %(verbose)s
 
@@ -83,7 +101,9 @@ class RawSNIRF(BaseRaw):
     """
 
     @verbose
-    def __init__(self, fname, optode_frame="unknown", preload=False, verbose=None):
+    def __init__(
+        self, fname, optode_frame="unknown", *, sfreq=None, preload=False, verbose=None
+    ):
         # Must be here due to circular import error
         from ...preprocessing.nirs import _validate_nirs_info
 
@@ -120,7 +140,7 @@ def __init__(self, fname, optode_frame="unknown", preload=False, verbose=None):
 
             last_samps = dat.get("/nirs/data1/dataTimeSeries").shape[0] - 1
 
-            sampling_rate = _extract_sampling_rate(dat)
+            sampling_rate = _extract_sampling_rate(dat, sfreq)
 
             if sampling_rate == 0:
                 warn("Unable to extract sample rate from SNIRF file.")
@@ -531,49 +551,48 @@ def _get_lengthunit_scaling(length_unit):
         )
 
 
-def _extract_sampling_rate(dat):
+def _extract_sampling_rate(dat, user_sfreq):
     """Extract the sample rate from the time field."""
     # This is a workaround to provide support for Artinis data.
     # It allows for a 1% variation in the sampling times relative
     # to the average sampling rate of the file.
     MAXIMUM_ALLOWED_SAMPLING_JITTER_PERCENTAGE = 1.0
 
+    _validate_type(user_sfreq, ("numeric", None), "sfreq")
     time_data = np.array(dat.get("nirs/data1/time"))
-    sampling_rate = 0
-    if len(time_data) == 2:
-        # specified as onset, samplerate
-        sampling_rate = 1.0 / (time_data[1] - time_data[0])
+    time_unit = _get_metadata_str(dat, "TimeUnit")
+    time_unit_scaling = _get_timeunit_scaling(time_unit)  # always 1 (s) or 1000 (ms)
+    if len(time_data) == 2:  # special-cased in the snirf standard as (onset, period)
+        onset, period = time_data
+        file_sfreq = time_unit_scaling / period
     else:
-        # specified as time points
+        onset = time_data[0]
         periods = np.diff(time_data)
-        uniq_periods = np.unique(periods.round(decimals=4))
-        if uniq_periods.size == 1:
-            # Uniformly sampled data
-            sampling_rate = 1.0 / uniq_periods.item()
-        else:
-            # Hopefully uniformly sampled data with some precision issues.
-            # This is a workaround to provide support for Artinis data.
-            mean_period = np.mean(periods)
-            sampling_rate = 1.0 / mean_period
-            ideal_times = np.linspace(time_data[0], time_data[-1], time_data.size)
-            max_jitter = np.max(np.abs(time_data - ideal_times))
-            percent_jitter = 100.0 * max_jitter / mean_period
-            msg = (
-                f"Found jitter of {percent_jitter:3f}% in sample times. Sampling "
-                f"rate has been set to {sampling_rate:1f}."
+        sfreqs = time_unit_scaling / periods
+        file_sfreq = sfreqs.mean()  # our best estimate, likely including some jitter
+    if user_sfreq is not None:
+        logger.info(f"Setting sampling frequency to user-supplied value: {user_sfreq}")
+        if not np.allclose(file_sfreq, user_sfreq, rtol=0.01, atol=0):
+            warn(
+                f"User-supplied sampling frequency ({user_sfreq} Hz) differs by "
+                f"{(user_sfreq - file_sfreq) / file_sfreq:.1%} from the frequency "
+                f"estimated from data in the file ({file_sfreq} Hz)."
             )
-            if percent_jitter > MAXIMUM_ALLOWED_SAMPLING_JITTER_PERCENTAGE:
-                warn(
-                    f"{msg} Note that MNE-Python does not currently support SNIRF "
-                    "files with non-uniformly-sampled data."
-                )
-            else:
-                logger.info(msg)
-    time_unit = _get_metadata_str(dat, "TimeUnit")
-    time_unit_scaling = _get_timeunit_scaling(time_unit)
-    sampling_rate *= time_unit_scaling
-
-    return sampling_rate
+    sfreq = user_sfreq or file_sfreq  # user-passed value overrides value from file
+    # estimate jitter
+    if len(time_data) > 2:
+        ideal_times = onset + np.arange(len(time_data)) / sfreq
+        max_jitter = np.max(np.abs(time_data - ideal_times))
+        percent_jitter = 100.0 * max_jitter / periods.mean()
+        msg = f"Found jitter of {percent_jitter:3f}% in sample times."
+        if percent_jitter > MAXIMUM_ALLOWED_SAMPLING_JITTER_PERCENTAGE:
+            warn(
+                f"{msg} Note that MNE-Python does not currently support SNIRF "
+                "files with non-uniformly-sampled data."
+            )
+        else:
+            logger.info(msg)
+    return sfreq
 
 
 def _get_metadata_str(dat, field):

mne/io/snirf/tests/test_snirf.py

@@ -4,6 +4,7 @@
 
 import datetime
 import shutil
+from contextlib import nullcontext
 
 import numpy as np
 import pytest
@@ -422,7 +423,7 @@ def test_snirf_kernel_hb():
     assert raw.copy().pick("hbo")._data.shape == (180, 14)
     assert raw.copy().pick("hbr")._data.shape == (180, 14)
 
-    assert_allclose(raw.info["sfreq"], 8.257638)
+    assert_allclose(raw.info["sfreq"], 8.256495)
 
     bad_nans = np.isnan(raw.get_data()).any(axis=1)
     assert np.sum(bad_nans) == 20
@@ -434,6 +435,23 @@ def test_snirf_kernel_hb():
     assert raw.annotations.description[1] == "StartIti"
 
 
+@requires_testing_data
+@pytest.mark.parametrize(
+    "sfreq,context",
+    (
+        [8.2, nullcontext()],  # sfreq estimated from file is 8.256495
+        [22, pytest.warns(RuntimeWarning, match="User-supplied sampling frequency")],
+    ),
+)
+def test_user_set_sfreq(sfreq, context):
+    """Test manually setting sfreq."""
+    with context:
+        # both sfreqs are far enough from true rate to yield >1% jitter
+        with pytest.warns(RuntimeWarning, match=r"jitter of \d+\.\d*% in sample times"):
+            raw = read_raw_snirf(kernel_hb, preload=False, sfreq=sfreq)
+    assert raw.info["sfreq"] == sfreq
+
+
 @requires_testing_data
 @pytest.mark.parametrize(
     "fname, boundary_decimal, test_scaling, test_rank",
@@ -535,8 +553,8 @@ def test_sample_rate_jitter(tmp_path):
     with h5py.File(new_file, "r+") as f:
         orig_time = np.array(f.get("nirs/data1/time"))
         acceptable_time_jitter = orig_time.copy()
-        average_time_diff = np.mean(np.diff(orig_time))
-        acceptable_time_jitter[-1] += 0.0099 * average_time_diff
+        mean_period = np.mean(np.diff(orig_time))
+        acceptable_time_jitter[-1] += 0.0099 * mean_period
         del f["nirs/data1/time"]
         f.flush()
         f.create_dataset("nirs/data1/time", data=acceptable_time_jitter)
@@ -545,11 +563,11 @@ def test_sample_rate_jitter(tmp_path):
     lines = "\n".join(line for line in log.getvalue().splitlines() if "jitter" in line)
     assert "Found jitter of 0.9" in lines
 
-    # Add jitter of 1.01%, which is greater than allowed tolerance
+    # Add jitter of 1.02%, which is greater than allowed tolerance
     with h5py.File(new_file, "r+") as f:
         unacceptable_time_jitter = orig_time
         unacceptable_time_jitter[-1] = unacceptable_time_jitter[-1] + (
-            0.0101 * average_time_diff
+            0.0102 * mean_period
         )
         del f["nirs/data1/time"]
         f.flush()