Refactor artifacts detection + add sturation detection by samuelgarcia · Pull Request #4297 · SpikeInterface/spikeinterface

samuelgarcia · 2026-01-06T16:00:18Z

@JoeZiminski @oliche

This refactor the artifact detection and blanking

with a main multi method function : detect_artifact_periods()
return periods with the new base_period_dtype instead list of list
SilencedPeriodsRecording also use this new dtype
DetectArtifactAndSilentPeriodsRecording (replace SilencedArtifactsRecording) do detect and silence at the same
time to be used in preprocess dict (@chrishalcrow this is for you)
refactor the method of Pierre based on envleope
incorporate the method of @oliche and @JoeZiminski to detect saturation

Could also be done in this PR:

add a margin_ms option in SilentPeriodsRecording to ewxtend the zeroing with a mragin : do we need this ?

To be done in 2 futures PR:

handle margin and taper on border of SilentPeriod (more tricky that it appear when period is at border...)
rewrite RemoveArtifactsRecording based on SilencedPeriodsRecording and the new dtype

yger · 2026-01-12T14:38:15Z

@samuelgarcia do not forget, in this PR, to change DetectThresholdCrossing(recording, ...) to DetectThresholdCrossing(envelope, ...) in the detec_period_artefacts_by_envelope

samuelgarcia · 2026-01-21T10:03:45Z

Hi @oliche @JoeZiminski
I copy paste and clean a bit your PR #4301 into this one.
Have a look and tell if this is OK.

I think we do not need to scaled traces at each chunk we could only could the unscaled threhold in the init and then we could avoid the costly convertion to float32 or float64.

Also I put the unit in uV to be more SI friendly.

… compatibility

samuelgarcia · 2026-01-21T13:22:02Z

@alejoe91 @chrishalcrow ready to review on my side

src/spikeinterface/preprocessing/detect_artifacts.py

samuelgarcia · 2026-01-22T10:49:51Z

src/spikeinterface/preprocessing/silence_periods.py

+            self, recording, periods=artifacts, mode=mode, noise_levels=noise_levels, seed=seed, **noise_levels_kwargs
+        )
+        # note self._kwargs["periods"] is done by SilencedPeriodsRecording and so the computaion is done once
+


add detect_artifact_method detect_artifact_kwargs

chrishalcrow · 2026-01-22T16:05:20Z

src/spikeinterface/preprocessing/detect_artifacts.py

+      * 'saturation' using detect_artifact_periods_by_envelope()
+      * 'envelope' using detect_saturation_periods()


Suggested change

* 'saturation' using detect_artifact_periods_by_envelope()

* 'envelope' using detect_saturation_periods()

* 'saturation' using detect_saturation_periods()

* 'envelope' using detect_artifact_periods_by_envelope()

chrishalcrow · 2026-01-22T16:06:35Z

src/spikeinterface/preprocessing/detect_artifacts.py

+# from spikeinterface.core.core_tools import define_function_handling_dict_from_class
+# from spikeinterface.preprocessing.silence_periods import SilencedPeriodsRecording


Suggested change

# from spikeinterface.core.core_tools import define_function_handling_dict_from_class

# from spikeinterface.preprocessing.silence_periods import SilencedPeriodsRecording

src/spikeinterface/preprocessing/detect_artifacts.py

chrishalcrow · 2026-01-22T16:44:12Z

Hello, trying this out on real data with saturation. I get this output

artifact_periods=array([(0, 637721, 638769), (0, 638955, 639206), (0, 923219, 923602),
       (0, 923605, 923606), (0, 923607, 923835)],

Are these periods inclusive? And if they are, should the last pair be merged?

chrishalcrow · 2026-01-22T16:47:22Z

src/spikeinterface/preprocessing/detect_artifacts.py

+
+def detect_saturation_periods(
+    recording,
+    saturation_threshold_uV,  # 1200 uV


Is it assumed that saturation_threshold_uV is symmetric around 0uV? E.g. you'll get saturation at $\pm 1200$uV? Is this always the case?

I think that's how saturation works, at least for int16. Maybe we can make a symmetric arguments?

Makes sense to me, @oliche? (we have a plus-minus saturation and voltage derivative threshold?

JoeZiminski · 2026-01-22T22:51:27Z

Hey @samuelgarcia thanks for this! will check this out tomorrow. I accidentally did not push two tidy-up commits to #4301, do you mind if I push to your branch to add some changes? I think it is mostly superficial stuff like docstrings etc.

samuelgarcia · 2026-01-23T16:59:41Z

Hello, trying this out on real data with saturation. I get this output
artifact_periods=array([(0, 637721, 638769), (0, 638955, 639206), (0, 923219, 923602),
       (0, 923605, 923606), (0, 923607, 923835)],
Are these periods inclusive? And if they are, should the last pair be merged?

Yes, the last pair should be merged.
Let us check. the _collapse_events() is done for this purpose.
@JoeZiminski can you look at this ?

samuelgarcia · 2026-01-23T17:00:10Z

Hey @samuelgarcia thanks for this! will check this out tomorrow. I accidentally did not push two tidy-up commits to #4301, do you mind if I push to your branch to add some changes? I think it is mostly superficial stuff like docstrings etc.

yes of course. push anything you want

for more information, see https://pre-commit.ci

JoeZiminski · 2026-01-30T13:47:39Z

Hey @samuelgarcia thanks for this, I added the content of those commits I missed. I ran into some issues with the (now renamed) uV_per_second_threshold . This is a threshold for the first derivative in units uV/s so it can be a very small number (e.g. IBL default in volts, converted to uV is 1e-2; for example, with the test data the threshold is around 0.0004 which goes negative after the round trip to int and back. Please correct me if wrong, the dynamic range of NP probes pm 5mV at int16 the bit resolution is 0.153 μV. I think this might be why @chrishalcrow is seeing those missed samples in his data.

If possible I'd prefer to cast and scale the data to avoid such issues, then we always know we are solid whatever the user inputs. I get 0.047298s to cast, offset and scale 10 million samples so I don't think this will be a bottleneck. WDYT?

samuelgarcia · 2026-02-02T12:36:31Z

src/spikeinterface/preprocessing/tests/test_detect_artifacts.py

    max_ = np.max(np.r_[data_seg_1.flatten(), data_seg_2.flatten()])
    gain =  max_ / 2**15
-    offset = 0
+    offset = 50


@JoeZiminski : why did you put 50 as offset here ?

Apologies forgot to say, just to ensure it is tested (if it was zero and was accidentally not applied in the function, the test would not flag it)

(that was just a quick test, see update below)

samuelgarcia · 2026-02-02T12:39:05Z

Please correct me if wrong, the dynamic range of NP probes pm 5mV at int16 the bit resolution is 0.153 μV. I think this might be why @chrishalcrow is seeing those missed samples in his data.

NP1 is 12bits, so I think that the lsb is 5000 / 2**11 = 2.4uV
Am I wrong ?

JoeZiminski · 2026-02-03T00:15:33Z

Hey @samuelgarcia I didn't realize the bit depth of the ADC is not 16-bit! I see here ADC resolution is given as 10 bits. In any case, the resolution from the gains/offsets used to scale the data won't be able to represent the values used for the uV_per_second_threshold?

oliche · 2026-02-03T15:42:25Z

@JoeZiminski the sample to volts conversion is a probe dependent variable. To get a feel of the complexity of the matter, here are the Neuropixel possible configs for all probes released to date: https://docs.google.com/spreadsheets/d/1UZycmxwwaJGCJ4AknyXrHkDXWOUO24TtAnBbe1mQ9u4/edit?gid=0#gid=0

As you can see there are potentially multiple conversion factors depending on the probe vintage.

It seems the SpikeInterface canon is to leave this up to the user. Our design in https://github.com/int-brain-lab/ibl-neuropixel has been to read it from the metadata, making sure we are compatible with the oldest probe versions we observed at IBL.

JoeZiminski · 2026-02-03T15:53:23Z

Thanks @oliche, in this case @samuelgarcia I think we should definitely scale the data rather than the parameters, as we cannot be sure the gains / offsets used to scale the data will not destroy information when converting float-to-int e.g. uV_per_second_threshold. We could play around with applying the division by fs at another point but it seems a little brittle and not worth the very minor speed-up.

alejoe91 · 2026-02-05T16:32:25Z

src/spikeinterface/preprocessing/detect_artifacts.py

+        # 0.98 is empirically determined as the true saturating point is
+        # slightly lower than the documented saturation point of the probe
+        self.saturation_threshold_unscaled = (thresh - offsets) / gains * 0.98
+        self.uV_per_sec_threshold = (uV_per_sec_threshold - offsets) / gains


It should be easy to transform this in uv_per_sample_threshold as:

self.uv_per_sample_threshold = (uV_per_sec_threshold * fs - offsets) / gains

alejoe91 · 2026-02-05T16:38:23Z

src/spikeinterface/preprocessing/detect_artifacts.py

+        self,
+        recording,
+        saturation_threshold_uV,
+        uV_per_sec_threshold,


I think it would be better to define uV_per_ms_threshold (or uv_per_sample as gfix does). The value in per seconds is very hard to grasp!

JoeZiminski · 2026-02-06T15:58:48Z

Hey @samuelgarcia @alejoe91 thanks for the feedback, I'm running into the below issue calculating the derivative on int16 data. I don't think this is an artefact of the test data but would be good to check, @oliche too.

See below for the test data in float and int. When taking the derivative, information is lost due to overflow.

The above is the data converted to int, each timeseries is the data from one channel. In many instances, the data after the saturation is negative in int space, resulting in overflow

e.g. by definition the saturation point is at the max of the int16 range, 32767. Often, when the saturation goes back to baseline, the datapoint is negative in int space e.g. -10000. When you take the diff e.g. -10000 - 32767 you get overflow and the information is lost. It seems like this would happen in real data as well?

for more information, see https://pre-commit.ci

oliche · 2026-02-07T13:05:04Z

@JoeZiminski I think the test data is not realistic. Usually the maximum int recorded on neuropixel is 512 on NP1 and 8192 on NP2 so you wouldn't have the issues above as the max diff would be 8192 * 2 (2 ** 14).

However in general I think it is a poor practice to compute in int16 for this kind of issues and patchy support of the numpy suite. Users may unwittingly create weird situations, exactly like this one: #4175.

Anyways for this specific saturation detection problem, you just have to fix your test data by staying in the maxint ranges above.

JoeZiminski · 2026-02-07T13:39:47Z

Cheers @oliche! That's good to know for NP probes. For other providers like Nexus, run their ADC at 16 bit-depth (e.g. here), in this case it might be a problem? That being said I don't know how often these probes saturate and whether people are likely to run this function for non-NP1 probes. But yes I think the easiest solution to avoid these worries is to perform computation in float.

Refactor aftifacts detection.

5610518

alejoe91 added the preprocessing Related to preprocessing module label Jan 6, 2026

alejoe91 changed the title ~~Refactor aftifacts detection.~~ Refactor artifacts detection Jan 7, 2026

JoeZiminski mentioned this pull request Jan 15, 2026

Compute saturation intervals ala IBL #4301

Draft

Integrate the saturation detector

ef3e820

samuelgarcia added 2 commits January 21, 2026 14:11

Progagate new periods dtype to SilencedPeriodsRecording with backward…

a71dade

… compatibility

oups

0d709e6

samuelgarcia changed the title ~~Refactor artifacts detection~~ Refactor artifacts detection + add sturation detection Jan 21, 2026

samuelgarcia marked this pull request as ready for review January 21, 2026 13:21

yger reviewed Jan 21, 2026

View reviewed changes

src/spikeinterface/preprocessing/detect_artifacts.py Show resolved Hide resolved

oups

7ab75c0

samuelgarcia commented Jan 22, 2026

View reviewed changes

chrishalcrow reviewed Jan 22, 2026

View reviewed changes

src/spikeinterface/preprocessing/detect_artifacts.py Outdated Show resolved Hide resolved

chrishalcrow reviewed Jan 22, 2026

View reviewed changes

JoeZiminski and others added 3 commits January 29, 2026 15:40

Applying extras from other PR, adding voltage_per_sec_threshold.

7c7446e

Rename uV_per_sec_threshold.

86e6924

[pre-commit.ci] auto fixes from pre-commit.com hooks

8b17890

for more information, see https://pre-commit.ci

samuelgarcia commented Feb 2, 2026

View reviewed changes

samuelgarcia added this to the 0.104.0 milestone Feb 3, 2026

alejoe91 reviewed Feb 5, 2026

View reviewed changes

JoeZiminski mentioned this pull request Feb 6, 2026

Suggestion to always cast the data to float and scale on load #4371

Open

JoeZiminski and others added 2 commits February 6, 2026 16:01

Updates testing converting to int.

4d507ad

[pre-commit.ci] auto fixes from pre-commit.com hooks

6546166

for more information, see https://pre-commit.ci

		* 'saturation' using detect_artifact_periods_by_envelope()
		* 'envelope' using detect_saturation_periods()

		# from spikeinterface.core.core_tools import define_function_handling_dict_from_class
		# from spikeinterface.preprocessing.silence_periods import SilencedPeriodsRecording

Conversation

samuelgarcia commented Jan 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

yger commented Jan 12, 2026

Uh oh!

samuelgarcia commented Jan 21, 2026

Uh oh!

samuelgarcia commented Jan 21, 2026

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

chrishalcrow commented Jan 22, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JoeZiminski commented Jan 22, 2026

Uh oh!

samuelgarcia commented Jan 23, 2026

Uh oh!

samuelgarcia commented Jan 23, 2026

Uh oh!

JoeZiminski commented Jan 30, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

samuelgarcia commented Feb 2, 2026

Uh oh!

JoeZiminski commented Feb 3, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

oliche commented Feb 3, 2026

Uh oh!

JoeZiminski commented Feb 3, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JoeZiminski commented Feb 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

oliche commented Feb 7, 2026

Uh oh!

JoeZiminski commented Feb 7, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

6 participants

samuelgarcia commented Jan 6, 2026 •

edited

Loading

JoeZiminski commented Feb 3, 2026 •

edited

Loading

JoeZiminski commented Feb 6, 2026 •

edited

Loading