Implemented sd_ratio as quality metric

doc/modules/qualitymetrics.rst

@@ -40,6 +40,7 @@ For more details about each metric and it's availability and use within SpikeInt
   qualitymetrics/snr
   qualitymetrics/noise_cutoff
   qualitymetrics/silhouette_score
+  qualitymetrics/sd_ratio
   qualitymetrics/synchrony
 
 

doc/modules/qualitymetrics/references.rst

@@ -21,6 +21,8 @@ References
 
 .. [Llobet] Llobet Victor, Wyngaard Aurélien and Barbour Boris. “Automatic post-processing and merging of multiple spike-sorting analyses with Lussac“. BioRxiv (2022).
 
+.. [Pouzat] Pouzat Christophe, Mazor Ofer and Laurent Gilles. “Using noise signature to optimize spike-sorting and to assess neuronal classification quality“. Journal of Neuroscience Methods (2002).
+
 .. [Rousseeuw] Peter J Rousseeuw. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of computational and applied mathematics, 20(C):53–65, 1987.
 
 .. [Schmitzer-Torbert]  Schmitzer-Torbert, Neil, and A. David Redish. “Neuronal Activity in the Rodent Dorsal Striatum in Sequential Navigation: Separation of Spatial and Reward Responses on the Multiple T Task.” Journal of neurophysiology 91.5 (2004): 2259–2272. Web.

doc/modules/qualitymetrics/sd_ratio.rst

@@ -0,0 +1,38 @@
+Standard Deviation (SD) ratio (:code:`sd_ratio`)
+==============================================
+
+Calculation
+-----------
+
+All spikes from the same neuron should have the same shape. This means that at the peak of the spike, the standard deviation of the voltage should be the same as that of noise. If spikes from multiple neurons are grouped in a single unit, chances are the standard deviation of spikes amplitude is going to increase.
+
+This metric, first described [Pouzat]_ then adapted by Wyngaard, Llobet & Barbour (paper in writing), returns the ratio between both standard deviations:
+
+.. math::
+	S = \frac{\sigma_{\mathrm{unit}}}{\sigma_{\mathrm{noise}}}
+
+To remove the effect of drift on spikes amplitude, :math:`\sigma_{\mathrm{unit}}` is computed by subtracting each spike amplitude, and dividing the resulting standard deviation by :math:`\sqrt{2}`.
+Also to remove the effect of bursts (which can have lower amplitudes), you can specify a censored period (by default 4.0 ms) where spikes happening less than this period after another spike will not be considered.
+
+
+Expectation and use
+-------------------
+
+For a unit representing a single neuron, this metric should return a value close to one. However for unit that are contaminated, the value can be significantly higher.
+
+
+Example code
+------------
+
+.. code-block:: python
+
+	import spikeinterface.qualitymetrics as sqm
+
+	sd_ratio = sqm.compute_sd_ratio(wvf_extractor, censored_period_ms=4.0)
+
+
+Literature
+----------
+
+Introduced by [Pouzat]_ (2002).
+Expanded by Wyngaard, Llobet and Barbour (in writing).

src/spikeinterface/qualitymetrics/misc_metrics.py

@@ -13,8 +13,8 @@
 import numpy as np
 import warnings
 
-from ..postprocessing import correlogram_for_one_segment
-from ..core import get_noise_levels
+from ..postprocessing import compute_spike_amplitudes, correlogram_for_one_segment
+from ..core import WaveformExtractor, get_noise_levels
 from ..core.template_tools import (
     get_template_extremum_channel,
     get_template_extremum_amplitude,
@@ -1365,3 +1365,68 @@ def _compute_rp_violations_numba(nb_rp_violations, spike_trains, spike_clusters,
             spike_train = spike_trains[spike_clusters == i]
             n_v = _compute_nb_violations_numba(spike_train, t_r)
             nb_rp_violations[i] += n_v
+
+
+def compute_sd_ratio(
+    wvf_extractor: WaveformExtractor,
+    censored_period_ms: float = 4.0,
+    correct_for_drift: bool = True,
+    unit_ids=None,
+    **kwargs,
+):
+    """
+    Computes the SD (Standard Deviation) of each unit spikes amplitude, and compare it to that of noise.
+
+    Parameters
+    ----------
+    waveform_extractor : WaveformExtractor
+        The waveform extractor object.
+    censored_period_ms : float
+        The censored period in milliseconds. This is to remove any potential bursts that could affect the SD.
+    correct_for_drift: bool
+        If True, will subtract the amplitudes sequentiially to significantly reduce the impact of drift.
+    unit_ids : list or None
+        The list of unit ids to compute this metric. If None, all units are used.
+    **kwargs:
+        Keyword arguments for computing spike amplitudes and extremum channel.
+    TODO: Possibly remove spikes when computing noise?
+    TODO: Take jitter into account.
+
+    Returns
+    -------
+    num_spikes : dict
+        The number of spikes, across all segments, for each unit ID.
+    """
+
+    from ..curation.curation_tools import _find_duplicated_spikes_keep_first_iterative
+
+    censored_period = int(round(censored_period_ms * 1e-3 * wvf_extractor.sampling_frequency))
+    if unit_ids is None:
+        unit_ids = wvf_extractor.unit_ids
+
+    spikes_amplitude = compute_spike_amplitudes(wvf_extractor, outputs="by_unit", return_scaled=True, **kwargs)
+    noise_levels = get_noise_levels(wvf_extractor.recording, return_scaled=True, method="std")
+    best_channels = get_template_extremum_channel(wvf_extractor, outputs="index", **kwargs)
+
+    sd_ratio = []
+    for segment_index in range(len(spikes_amplitude)):
+        sd_ratio.append({})
+
+        for unit_id, spk_amp in spikes_amplitude[segment_index].items():
+            if unit_id not in unit_ids:
+                continue
+
+            spike_train = wvf_extractor.sorting.get_unit_spike_train(unit_id, segment_index=segment_index).astype(
+                np.int64
+            )
+            censored_indices = _find_duplicated_spikes_keep_first_iterative(spike_train, censored_period)
+            spk_amp = np.delete(spk_amp, censored_indices)
+
+            if correct_for_drift:
+                unit_std = np.std(np.diff(spk_amp)) / np.sqrt(2)
+            else:
+                unit_std = np.std(spk_amp)
+
+            sd_ratio[segment_index][unit_id] = unit_std / noise_levels[best_channels[unit_id]]
+
+    return sd_ratio

src/spikeinterface/qualitymetrics/quality_metric_list.py

@@ -14,6 +14,7 @@
     compute_synchrony_metrics,
     compute_firing_ranges,
     compute_amplitude_cv_metrics,
+    compute_sd_ratio,
 )
 
 from .pca_metrics import (
@@ -46,4 +47,5 @@
     "synchrony": compute_synchrony_metrics,
     "firing_range": compute_firing_ranges,
     "drift": compute_drift_metrics,
+    "sd_ratio": compute_sd_ratio,
 }

src/spikeinterface/qualitymetrics/tests/test_metrics_functions.py

@@ -34,6 +34,7 @@
     compute_synchrony_metrics,
     compute_firing_ranges,
     compute_amplitude_cv_metrics,
+    compute_sd_ratio,
 )
 
 
@@ -70,7 +71,7 @@ def _simulated_data():
 
 
 def _waveform_extractor_simple():
-    recording, sorting = toy_example(duration=50, seed=10)
+    recording, sorting = toy_example(duration=80, seed=10)
     recording = recording.save(folder=cache_folder / "rec1")
     sorting = sorting.save(folder=cache_folder / "sort1")
     folder = cache_folder / "waveform_folder1"
@@ -379,6 +380,14 @@ def test_calculate_drift_metrics(waveform_extractor_simple):
     # assert np.allclose(list(drift_mads_gt.values()), list(drift_mads.values()), rtol=0.05)
 
 
+def test_calculate_sd_ratio(waveform_extractor_simple):
+    sd_ratio = compute_sd_ratio(waveform_extractor_simple)
+
+    assert len(sd_ratio) == waveform_extractor_simple.get_num_segments()
+    assert np.all(list(sd_ratio[0].keys()) == waveform_extractor_simple.unit_ids)
+    assert np.allclose(np.array(list(sd_ratio[0].values())), 1, atol=0.5, rtol=0)
+
+
 if __name__ == "__main__":
     sim_data = _simulated_data()
     we = _waveform_extractor_simple()
@@ -390,5 +399,6 @@ def test_calculate_drift_metrics(waveform_extractor_simple):
     # test_calculate_sliding_rp_violations(we)
     # test_calculate_drift_metrics(we)
     # test_synchrony_metrics(we)
-    test_calculate_firing_range(we)
+    # test_calculate_firing_range(we)
     # test_calculate_amplitude_cv_metrics(we)
+    test_calculate_sd_ratio(we)