Filter-dependant voltage shift in BioCAM recordings

[b-grimaud]

Describe the bug
Raw traces from BioCAM recording have an unexpected + 2048 applied to the voltage measurements, i.e. traces are all shifted 2048 µV higher than they actually are. This issue disappears when the recording is filtered through SpikeInterface.

To Reproduce
Using either version of BRW recordings available on GIN :

from pathlib import Path
import matplotlib.pyplot as plt
import spikeinterface as si
import spikeinterface.extractors as se
import spikeinterface.preprocessing as spre

path = Path(r'/home/lumin/Documents/ephy/biocam_hw3.0_fw1.6.brw') # or biocam_hw3.0_fw1.7.0.12_raw.brw

recording = se.read_biocam(path, mea_pitch=81, electrode_width=21)
# recording = spre.bandpass_filter(recording=recording,\
#                                   freq_min=300, freq_max=3000)

traces = recording.get_traces(channel_ids=['35','26','15','23','8'],\
                               end_frame=4096) # Picking a few random channels
plt.plot(traces)

from pathlib import Path
import matplotlib.pyplot as plt
import spikeinterface as si
import spikeinterface.extractors as se
import spikeinterface.preprocessing as spre

path = Path(r'/home/lumin/Documents/ephy/biocam_hw3.0_fw1.6.brw') # or biocam_hw3.0_fw1.7.0.12_raw.brw

recording = se.read_biocam(path, mea_pitch=81, electrode_width=21)
recording = spre.bandpass_filter(recording=recording,\
                                  freq_min=300, freq_max=3000)

traces = recording.get_traces(channel_ids=['35','26','15','23','8'],\
                               end_frame=4096) # Picking a few random channels
plt.plot(traces)

Expected behaviour
Raw baseline should be around zero.

Environment:

• OS: Ubuntu
• Python version : 3.10.10
• Neo version : 0.13.0.dev0
• NumPy version : 1.24.2

Additional context
I first reported the problem in a SpikeInterface issue a while ago, and detailed it a bit more in this recent PR.

Briefly, looking through @mhhennig 's read function in HerdingSpikes, signal is either read as data - 2048 or 2048 - data depending on file version and signal inversion.
There also seem to be some kind of mask to filter out certain values ? Not sure about that : d[np.abs(d) > 1500] = 0.

On the other hand, neo reads signal as either data or 4096 - data, depending on signal inversion.
Furthermore, signal inversion itself isn't actually taken into account :

    if format_100:
        if signal_inv == 1:
            read_function = readHDF5t_100
        elif signal_inv == 1:
            read_function = readHDF5t_100_i
        else:
            raise Exception("Unknown signal inversion")
    else:
        if signal_inv == 1:
            read_function = readHDF5t_101
        elif signal_inv == 1:
            read_function = readHDF5t_101_i
        else:
            raise Exception("Unknown signal inversion")

By quickly running through random data point in the HDF5 file, it looks to me that the - 2048 is necessary :

On biocam_hw3.0_fw1.6.brw :

>>> f['3BData']['Raw'][0]
2043
>>> f['3BData']['Raw'][758]
2045
>>> f['3BData']['Raw'][92475]
2037
>>> 

On biocam_hw3.0_fw1.7.0.12_raw.brw :

>>> f['Well_A1']['Raw'][0]
2048
>>> f['Well_A1']['Raw'][725]
2046
>>> f['Well_A1']['Raw'][98712]
2050
>>> 

[alejoe91]

Thanks for writing this up @BptGrm !

The problem is gone after filtering because the filtering removes the offset.

Is this information about the offset clear in the brw format documentation?
Feel free to make a PR to fix the behavior.

I think I ported Matthias' code a while back so not sure where the wrong scaling came from!
We generally do not deal with biocam systems, so this might have gone under the radar.

[alejoe91]

For reference, the NEO code was ported from here: https://github.com/SpikeInterface/spikeextractors/blob/master/spikeextractors/extractors/biocamrecordingextractor/biocamrecordingextractor.py

[b-grimaud]

The problem is gone after filtering because the filtering removes the offset.

Makes sense !

Is this information about the offset clear in the brw format documentation?

Couldn't find any explicit mention of it the docs for both BRW v3.x and v4.x. The only code example provided deals with event-base data. 3Brain shared some code with me a while ago, but all of the IO is done by their own DLLs.

Feel free to make a PR to fix the behavior.

Will do, still not sure about the "mask" d[np.abs(d) > 1500] = 0 or if the data type is important .flatten('C').astype(ctypes.c_short).
I know detection in HerdingSpikes is done in C so that might be why it was there originally.

I can wait for the previous PR #1326 to merge so that I can do a single one myself.

For reference, the NEO code was ported from here:

Maybe there's a clue in here ?

    def write_recording(recording, save_path):
        # Convert to uV:
        # AnalogValue = MVOffset + DigitalValue * ADCCountsToMV
        # Where ADCCountsToMV is defined as:
        # ADCCountsToMV = SignalInversion * ((MaxVolt - MinVolt) / 2^BitDepth)
        # And MVOffset as:
        # MVOffset = SignalInversion * MinVolt
        # conversion back
        # DigitalValue = (AnalogValue - MVOffset)/ADCCountsToMV
        # we center at 2048

[b-grimaud]

See #1348