Merge pull request #3844 from h-mayorquin/metric_system_units_docs

Add documentation to handle (physical) units of recording channels

Author: alejoe91

doc/how_to/index.rst

@@ -17,4 +17,5 @@ Guides on how to solve specific, short problems in SpikeInterface. Learn how to.
     drift_with_lfp
     auto_curation_training
     auto_curation_prediction
+    physical_units
     customize_a_plot

doc/how_to/physical_units.rst

@@ -0,0 +1,107 @@
+Working with physical units in SpikeInterface recordings
+========================================================
+
+In neurophysiology recordings, data is often stored in raw ADC (Analog-to-Digital Converter) integer values but needs to be analyzed in physical units.
+For extracellular recordings, this is typically microvolts (µV), but some recording devices may use different physical units.
+SpikeInterface provides tools to handle both situations.
+
+It's important to note that **most spike sorters work fine on raw digital (ADC) units** and scaling is not needed. Going a step further, some sorters, such as Kilosort 3, require their input to be in raw ADC units.
+The specific behavior however depends on the spike sorter, so it's important to understand the specific input requirements on a case per case basis.
+
+Many preprocessing tools are also linear transformations, and if the ADC is implemented as a linear transformation which is fairly common, then the overall effect can be preserved.
+That is, **preprocessing steps can often be applied either before or after unit conversion without affecting the outcome.**. That being said, there are rough edges to this approach.
+preprocessing algorithms like filtering, whitening, centering, interpolation and common reference require casting to float within the pipeline. We advise users to experiment
+with different approaches to find the best one for their specific use case.
+
+
+Therefore, **it is usually safe to work in raw ADC integer values unless a specific tool or analysis requires physical units**.
+If you are interested in visualizations, comparability across devices, or outputs with interpretable physical scales (e.g., microvolts), converting to physical units is recommended.
+Otherwise, remaining in raw units can simplify processing and preserve performance.
+
+Understanding Physical Units
+----------------------------
+
+Most recording devices store data in ADC units (integers) to save space and preserve the raw data.
+To convert these values to physical units, two parameters are needed:
+
+* **gain**: A multiplicative factor to scale the raw values
+* **offset**: An additive factor to shift the values
+
+The conversion formula is:
+
+.. code-block:: text
+
+    physical_value = raw_value * gain + offset
+
+
+Converting to Physical Units
+----------------------------
+
+SpikeInterface provides two preprocessing classes for converting recordings to physical units. Both wrap the
+``RecordingExtractor`` class and ensures that the data is returned in physical units when calling `get_traces <https://spikeinterface.readthedocs.io/en/stable/api.html#spikeinterface.core.BaseRecording.get_traces>`_
+
+1. ``scale_to_uV``: The primary function for extracellular recordings. SpikeInterface is centered around
+    extracellular recordings, and this function is designed to convert the data to microvolts (µV).
+2. ``scale_to_physical_units``: A general function for any physical unit conversion. This will allow you to extract the data in any
+    physical unit, not just microvolts. This is useful for other types of recordings, such as force measurements in Newtons but should be
+    handled with care.
+
+For most users working with extracellular recordings, ``scale_to_uV`` is the recommended choice if they want to work in physical units:
+
+.. code-block:: python
+
+    from spikeinterface.extractors import read_intan
+    from spikeinterface.preprocessing import scale_to_uV
+
+    # Load recording (data is in ADC units)
+    recording = read_intan("path/to/file.rhs")
+
+    # Convert to microvolts
+    recording_uv = scale_to_uV(recording)
+
+For recordings with non-standard units (e.g., force measurements in Newtons), use ``scale_to_physical_units``:
+
+.. code-block:: python
+
+    from spikeinterface.preprocessing import scale_to_physical_units
+
+    # Convert to physical units (whatever they may be)
+    recording_physical = scale_to_physical_units(recording)
+
+Both preprocessors automatically:
+
+1. Detect the appropriate gain and offset from the recording properties
+2. Apply the conversion to all channels
+3. Update the recording properties to reflect that data is now in physical units
+
+Setting Custom Physical Units
+-----------------------------
+
+While most extractors automatically set the appropriate ``gain_to_uV`` and ``offset_to_uV`` values,
+there might be cases where you want to set custom physical units. In these cases, you can set
+the following properties:
+
+* ``physical_unit``: The target physical unit (e.g., 'uV', 'mV', 'N')
+* ``gain_to_unit``: The gain to convert from raw values to the target unit
+* ``offset_to_unit``: The offset to convert from raw values to the target unit
+
+You need to set these properties for every channel, which allows for the case when there are different gains and offsets on different channels. Here's an example:
+
+.. code-block:: python
+
+    # Set custom physical units
+    num_channels = recording.get_num_channels()
+    values = ["volts"] * num_channels
+    recording.set_property(key='physical_unit', values=values)
+
+    gain_values = [0.001] * num_channels  # Convert from ADC to volts
+    recording.set_property(key='gain_to_unit', values=gain_values)  # Convert to volts
+
+    offset_values = [0] * num_channels  # No offset
+    recording.set_property(key='offset_to_unit', values=offset_values)  # No offset
+
+    # Apply the conversion using scale_to_physical_units
+    recording_physical = scale_to_physical_units(recording)
+
+This approach gives you full control over the unit conversion process while maintaining
+compatibility with SpikeInterface's preprocessing pipeline.

src/spikeinterface/core/baserecording.py

@@ -20,7 +20,15 @@ class BaseRecording(BaseRecordingSnippets):
     """
 
     _main_annotations = BaseRecordingSnippets._main_annotations + ["is_filtered"]
-    _main_properties = ["group", "location", "gain_to_uV", "offset_to_uV"]
+    _main_properties = [
+        "group",
+        "location",
+        "gain_to_uV",
+        "offset_to_uV",
+        "gain_to_physical_unit",
+        "offset_to_physical_unit",
+        "physical_unit",
+    ]
     _main_features = []  # recording do not handle features
 
     _skip_properties = [

src/spikeinterface/preprocessing/preprocessinglist.py

@@ -25,8 +25,8 @@
     CenterRecording,
     center,
 )
+from .scale import scale_to_uV, scale_to_physical_units
 
-from .scale import scale_to_uV
 
 from .whiten import WhitenRecording, whiten, compute_whitening_matrix
 from .rectify import RectifyRecording, rectify

src/spikeinterface/preprocessing/scale.py

@@ -1,6 +1,70 @@
 from __future__ import annotations
 
+import numpy as np
+
 from spikeinterface.preprocessing.basepreprocessor import BasePreprocessor
+from spikeinterface.preprocessing.normalize_scale import ScaleRecording
+from spikeinterface.core.core_tools import define_function_handling_dict_from_class
+
+
+class ScaleToPhysicalUnits(ScaleRecording):
+    """
+    Scale raw traces to their physical units using gain_to_physical_unit and offset_to_physical_unit.
+
+    This preprocessor uses the channel-specific gain and offset information
+    stored in the recording extractor to convert the raw traces to their physical units.
+    Most commonly this will be microvolts (µV) for voltage recordings, but some extractors
+    might use different physical units (e.g., Newtons for force measurements).
+
+    Parameters
+    ----------
+    recording : BaseRecording
+        The recording extractor to be scaled. The recording extractor must
+        have gain_to_physical_unit and offset_to_physical_unit properties set.
+
+    Returns
+    -------
+    ScaleToPhysicalUnits
+        The recording with traces scaled to physical units.
+
+    Raises
+    ------
+    ValueError
+        If the recording extractor does not have gain_to_physical_unit and offset_to_physical_unit properties.
+    """
+
+    name = "recording_in_physical_units"
+
+    def __init__(self, recording):
+        if "gain_to_physical_unit" not in recording.get_property_keys():
+            error_msg = (
+                "Recording must have 'gain_to_physical_unit' property to convert to physical units. \n"
+                "Set the gain using `recording.set_property(key='gain_to_physical_unit', values=values)`."
+            )
+            raise ValueError(error_msg)
+        if "offset_to_physical_unit" not in recording.get_property_keys():
+            error_msg = (
+                "Recording must have 'offset_to_physical_unit' property to convert to physical units. \n"
+                "Set the offset using `recording.set_property(key='offset_to_physical_unit', values=values)`."
+            )
+            raise ValueError(error_msg)
+
+        gain = recording.get_property("gain_to_physical_unit")
+        offset = recording.get_property("offset_to_physical_unit")
+
+        # Initialize parent ScaleRecording with the gain and offset values
+        ScaleRecording.__init__(self, recording, gain=gain, offset=offset, dtype="float32")
+
+        # Reset gain/offset since data is now in physical units
+        self.set_property(key="gain_to_physical_unit", values=np.ones(recording.get_num_channels(), dtype="float32"))
+        self.set_property(key="offset_to_physical_unit", values=np.zeros(recording.get_num_channels(), dtype="float32"))
+
+        # Also reset channel gains and offsets
+        self.set_channel_gains(gains=1.0)
+        self.set_channel_offsets(offsets=0.0)
+
+
+scale_to_physical_units = define_function_handling_dict_from_class(ScaleToPhysicalUnits, name="scale_to_physical_units")
 
 
 def scale_to_uV(recording: BasePreprocessor) -> BasePreprocessor:

src/spikeinterface/preprocessing/tests/test_scaling.py

@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
-from spikeinterface.core.testing_tools import generate_recording
-from spikeinterface.preprocessing import scale_to_uV, CenterRecording
+from spikeinterface.core.generate import generate_recording
+from spikeinterface.preprocessing import scale_to_uV, CenterRecording, scale_to_physical_units
 
 
 def test_scale_to_uV():
@@ -70,6 +70,48 @@ def test_scaling_in_preprocessing_chain():
     np.testing.assert_allclose(traces_scaled_with_preprocessor, traces_scaled_with_preprocessor_and_argument)
 
 
+def test_scale_to_physical_units():
+    # Create a sample recording extractor with fake physical unit gains and offsets
+    num_channels = 4
+    sampling_frequency = 30_000.0
+    durations = [1.0]  # seconds
+    recording = generate_recording(
+        num_channels=num_channels,
+        durations=durations,
+        sampling_frequency=sampling_frequency,
+    )
+
+    rng = np.random.default_rng(0)
+    gains = rng.random(size=(num_channels)).astype(np.float32)
+    offsets = rng.random(size=(num_channels)).astype(np.float32)
+
+    # Set physical unit gains/offsets instead of regular gains/offsets
+    recording.set_property("gain_to_physical_unit", gains)
+    recording.set_property("offset_to_physical_unit", offsets)
+
+    # Apply the preprocessor
+    scaled_recording = scale_to_physical_units(recording=recording)
+
+    # Get raw traces and apply scaling manually
+    raw_traces = recording.get_traces(segment_index=0)
+    expected_traces = raw_traces * gains + offsets
+
+    # Get scaled traces
+    scaled_traces = scaled_recording.get_traces(segment_index=0)
+
+    # Check if the traces are scaled correctly
+    np.testing.assert_allclose(scaled_traces, expected_traces)
+
+    # Test for the error when recording doesn't have physical unit properties
+    recording_no_gains = generate_recording(
+        num_channels=num_channels,
+        durations=durations,
+        sampling_frequency=sampling_frequency,
+    )
+    with pytest.raises(ValueError):
+        scale_to_physical_units(recording_no_gains)
+
+
 if __name__ == "__main__":
     test_scale_to_uV()
     test_scaling_in_preprocessing_chain()