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

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.

@alejoe91 @chrishalcrow ready to review on my side

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

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

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.

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 ?

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

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 mV. 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?

Please correct me if wrong, the dynamic range of NP probes pm 5mV at int16 the bit resolution is 0.153 mV. 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 ?

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?

@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.

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.

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?

@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.

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.

If you think you may exceed half the the range of int16 (ie higher than 2 ** 14 or lower than -2 ** 14), then you can check for the max(abs()) values. If they are above or below, you may throw a warning and only perform the absolute saturation detection and skip the one based on the derivative ?

That would prevent too many false positives.

Thanks @oliche, yes in this case I would definitely suggest @samuelgarcia we perform operations in float here to avoid having to worry about the unpredictable bit-depth of the data.