Map EDF digital endpoints to their physical endpoints

doc/source/authors.rst

@@ -106,6 +106,7 @@ and may not be the current affiliation of a contributor.
 * Aditya Singh (github)
 * Stanley C.
 * Axel Coiffard-Donaty
+* Anton Soloviev
 
 1. Centre de Recherche en Neuroscience de Lyon, CNRS UMR5292 - INSERM U1028 - Université Claude Bernard Lyon 1
 2. Unité de Neuroscience, Information et Complexité, CNRS UPR 3293, Gif-sur-Yvette, France

neo/rawio/edfrawio.py

@@ -105,8 +105,8 @@ def _parse_header(self):
             dtype = "int16"  # assume general int16 based on edf documentation
             units = sig_dict["dimension"]
             physical_range = sig_dict["physical_max"] - sig_dict["physical_min"]
-            # number of digital steps resolved (+1 to account for '0')
-            digital_range = sig_dict["digital_max"] - sig_dict["digital_min"] + 1
+            # The digital endpoints correspond to the physical endpoints.
+            digital_range = sig_dict["digital_max"] - sig_dict["digital_min"]
             gain = physical_range / digital_range
             offset = -1 * sig_dict["digital_min"] * gain + sig_dict["physical_min"]
 

neo/test/iotest/test_edfio.py

@@ -67,6 +67,8 @@ def test_compare_data(self):
         """
         Compare data from AnalogSignal with plain data stored in text file
         """
+        import pyedflib
+
         with EDFIO(self.filename) as io:
             plain_data = np.loadtxt(io.filename.replace(".edf", ".txt"), dtype=np.int16)
             seg = io.read_segment(lazy=True)
@@ -78,16 +80,11 @@ def test_compare_data(self):
                 ana_data = anasig.load(magnitude_mode="raw")
                 np.testing.assert_array_equal(ana_data.magnitude, plain_data[:, aidx : aidx + 1])
 
-                # comparing floating data to original values * gain factor
-                ch_head = io.edf_reader.getSignalHeader(aidx)
-                physical_range = ch_head["physical_max"] - ch_head["physical_min"]
-                # number of digital values used (+1 to account for '0' value)
-                digital_range = ch_head["digital_max"] - ch_head["digital_min"] + 1
-
-                gain = physical_range / digital_range
+                # Compare calibrated values with pyedflib's physical readout.
                 ana_data = anasig.load(magnitude_mode="rescaled")
-                rescaled_data = plain_data[:, aidx : aidx + 1] * gain
-                np.testing.assert_array_equal(ana_data.magnitude, rescaled_data)
+                with pyedflib.EdfReader(io.filename) as reader:
+                    rescaled_data = reader.readSignal(aidx)
+                np.testing.assert_allclose(ana_data.magnitude[:, 0], rescaled_data, rtol=1e-6, atol=1e-5)
 
 
 if __name__ == "__main__":

neo/test/rawiotest/test_edfrawio.py

@@ -3,10 +3,68 @@
 """
 
 import unittest
+import tempfile
+from pathlib import Path
 
-from neo.rawio.edfrawio import EDFRawIO
+import numpy as np
+
+from neo.rawio.edfrawio import EDFRawIO, HAS_PYEDF
 from neo.test.rawiotest.common_rawio_test import BaseTestRawIO
 
+if HAS_PYEDF:
+    import pyedflib
+
+
+@unittest.skipUnless(HAS_PYEDF, "requires pyedflib")
+class TestEDFCalibration(unittest.TestCase):
+    def test_digital_endpoints_map_to_physical_endpoints(self):
+        headers = [
+            dict(
+                label="identity",
+                dimension="uV",
+                sample_frequency=8,
+                digital_min=-1000,
+                digital_max=1000,
+                physical_min=-1000,
+                physical_max=1000,
+            ),
+            dict(
+                label="offset",
+                dimension="uV",
+                sample_frequency=8,
+                digital_min=0,
+                digital_max=100,
+                physical_min=-10,
+                physical_max=30,
+            ),
+            dict(
+                label="full_range",
+                dimension="uV",
+                sample_frequency=8,
+                digital_min=-32768,
+                digital_max=32767,
+                physical_min=-200,
+                physical_max=200,
+            ),
+        ]
+        digital = [np.array([h["digital_min"], h["digital_max"]] * 4, dtype=np.int32) for h in headers]
+        expected = np.array([[h["physical_min"], h["physical_max"]] * 4 for h in headers]).T
+        with tempfile.TemporaryDirectory() as folder:
+            filename = str(Path(folder) / "calibration.edf")
+            with pyedflib.EdfWriter(filename, len(headers), file_type=pyedflib.FILETYPE_EDFPLUS) as writer:
+                writer.setSignalHeaders(headers)
+                writer.writeSamples(digital, digital=True)
+            with EDFRawIO(filename) as reader:
+                reader.parse_header()
+                for channels in (None, [2, 0], [1]):
+                    raw = reader.get_analogsignal_chunk(channel_indexes=channels)
+                    selected = slice(None) if channels is None else channels
+                    np.testing.assert_array_equal(raw, np.asarray(digital).T[:, selected])
+                    for dtype in ("float32", "float64"):
+                        with self.subTest(channels=channels, dtype=dtype):
+                            physical = reader.rescale_signal_raw_to_float(raw, dtype=dtype, channel_indexes=channels)
+                            np.testing.assert_allclose(physical, expected[:, selected], rtol=1e-6, atol=1e-5)
+
 
 class TestExampleRawIO(
     BaseTestRawIO,