_extract_sampling_rate function for io.read_raw_snirf incorrectly estimates sample rate due to floating point error.

[kdarti]

Description of the problem

I have a snirf file with a sample rate of 50hz, the /nirs/data1/time array is correctly incremented by 0.02 per sample, as expected.

However, when mne.io_read_raw_snirf() is used to read the file, the sample rate is estimated to be not 50Hz:

This is due to floating errors introduced in this section:

mne-python/mne/io/snirf/_snirf.py

Lines 542 to 558 in 415e7f6

	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}."
	)

More precisely, by .np.diff() on line 541.

Easiest fix would be to use Decimal to avoid these floating point errors. using float64 should also work, at least in this case.

Steps to reproduce

import mne

mne.io.read_raw_snirf(r"WR11 Baseline SNIRF.snirf")

Link to data

WR11 Baseline SNIRF.zip

Expected results

Sample rate of 50Hz

Actual results

49.999996Hz

Additional information

Platform Windows-10-10.0.19041-SP0
Python 3.8.3 (default, Jul 2 2020, 17:30:36) [MSC v.1916 64 bit (AMD64)]
Executable C:\Users\kdahlslatt\Anaconda3\python.exe
CPU Intel64 Family 6 Model 158 Stepping 10, GenuineIntel (12 cores)
Memory 15.7 GB

Core
├☑ mne 1.6.1 (latest release)
├☑ numpy 1.24.3 (MKL 2023.1-Product with 6 threads)
├☑ scipy 1.10.1
├☑ matplotlib 3.7.1 (backend=Qt5Agg)
C:\Users\kdahlslatt\Anaconda3\lib\site-packages\paramiko\transport.py:219: CryptographyDeprecationWarning: Blowfish has been deprecated
"class": algorithms.Blowfish,
├☑ pooch 1.8.0
└☑ jinja2 3.1.2

Numerical (optional)
├☑ sklearn 1.3.0
├☑ numba 0.57.1
├☑ nibabel 3.2.1
├☑ nilearn 0.9.1
├☑ pandas 1.5.3
└☐ unavailable dipy, openmeeg, cupy

Visualization (optional)
├☑ pyvista 0.37.0 (OpenGL 4.5.0 - Build 30.0.101.1404 via Intel(R) UHD Graphics 630)
├☑ vtk 9.2.5
├☑ qtpy 2.2.0 (PyQt5=5.15.2)
├☑ ipympl 0.9.3
├☑ pyqtgraph 0.11.0
├☑ ipywidgets 8.0.4
└☐ unavailable pyvistaqt, mne-qt-browser, trame_client, trame_server, trame_vtk, trame_vuetify

Ecosystem (optional)
└☐ unavailable mne-bids, mne-nirs, mne-features, mne-connectivity, mne-icalabel, mne-bids-pipeline

[larsoner]

@kdarti agreed ideally we'd get 50 here, are you up for trying some fix?

[drammock]

IIRC, @rob-luke and I wrote this code; if you're up for making a PR yourself to fix this, feel free to ping me for reviews!

[kdarti]

I tried simply using np.float64, but there are still floating point errors. The option of using Decimal package isn't great because it is very very slow (in a quick test, 750x slower with default precision). The remaining option would be to upscale to int and then later downscale, but I don't really have time to do this unfortunately.