Issue preprocessing: filtered traces after saving show artifacts

[letiziasignorelli]

Hi, I'm having issues when preprocessing my recording and saving with recording.save function: the trace is not the same before and after saving, it seems that the saving function is introducing artifacts. This happens when applying a bandpass filter in my case.

I'm using spikeinterface version 0.101.2

import pathlib, os
import spikeinterface.extractors as se
import spikeinterface.preprocessing as spre
import spikeinterface.full as si
import matplotlib.pyplot as plt

selected_file_path = pathlib.Path("test.bin")

freq_min = 5
freq_max = 200

recording = se.read_axona(selected_file_path)
recordingf = spre.bandpass_filter(recording, freq_min=freq_min, freq_max=freq_max)

# Save
n_cpus = os.cpu_count()
n_jobs = n_cpus - 4
job_kwargs = dict(n_jobs=n_jobs, progress_bar=True)
save_filename = "save_filter_chunk_1s" 
saved_recording = recordingf.save(folder = save_filename, name = save_filename, **job_kwargs)

# Compare recording
reloaded_recording = si.load_extractor(save_filename)
original_trace = recordingf.get_traces()[:,0]
saved_trace = saved_recording.get_traces()[:,0]
reloaded_trace = reloaded_recording.get_traces()[:,0]

delta = original_trace - saved_trace
delta1 = original_trace - reloaded_trace
delta2 = saved_trace - reloaded_trace

fig, ax = plt.subplots(3, 1)
fig.set_figwidth(20)
fig.set_figheight(8)
ax[0].plot(delta)
ax[1].plot(delta1)
ax[2].plot(delta2)

fig, ax = plt.subplots(3, 1)
fig.set_figwidth(20)
fig.set_figheight(8)
ax[0].plot(original_trace)
ax[1].plot(saved_trace)
ax[2].plot(reloaded_trace)

I then tried with chunk_duration = recording.get_duration() since it seemed to have these artifacts every 1 second, and it seemed to work, but since chunk_duration default value is 1s maybe there's something in the saving after filtering that is not working properly?

save_filename = "save_filter_chunk_duration" 
saved_recording = recordingf.save(folder = save_filename, name = save_filename, chunk_duration = recording.get_duration(), **job_kwargs)

# Compare recordings
reloaded_recording = si.load_extractor(save_filename)
original_trace = recordingf.get_traces()[:,0]
saved_trace = saved_recording.get_traces()[:,0]
reloaded_trace = reloaded_recording.get_traces()[:,0]

delta = original_trace - saved_trace
delta1 = original_trace - reloaded_trace
delta2 = saved_trace - reloaded_trace

fig, ax = plt.subplots(3, 1)
fig.set_figwidth(20)
fig.set_figheight(8)
ax[0].plot(delta)
ax[1].plot(delta1)
ax[2].plot(delta2)

Finally, I also tried to just load and save the recording without applying any preprocessing and it does not show this issue.
For reference I uploaded my scripts and data used here for your testing: https://drive.google.com/drive/folders/1h0clLnAmt0-R2TO4yP1rAJwXZW5WgrKC

[chrishalcrow]

Hi @letiziasignorelli , thanks for raising this. I've reproduced from your data and also reproduced for my own data.

Problem has something to do with freq_min being low, so is to do with the highpass_filter part of the bandpass_filter. The artefacts get smaller and smaller as freq_min gets larger and larger. I remember @samuelgarcia discussing the difficulty of parallelising highpass_filter. Maybe we approximately do this at the borders, and the approximation breaks down as 1/freq_min approaches the chunk size?

[zm711]

Yeah Sam had said that low frequency are poorly handled with our parallelization. So you can improve this slightly by increasing chunk size but @samuelgarcia had said we can never truly fix it. Maybe we should document this better since it has now come up in an actual issue.

[letiziasignorelli]

I see, thank you for your answer. I think document this would be helpful. In my case I'm building a pipeline where I need both high and low frequencies for different purposes, and the spikeinterface objects are great to do handle both. Knowing that for the low frequencies there's this problem it's quite crucial

[zm711]

I'll add it to our agenda for a maintenance meeting. It's worth a very careful discussion on our end to make things as clear as possible for end-users.

[chrishalcrow]

Also worth noting that you realised the solution to the problem: to increase the chunk size relative to 1/(freq_min). We'll discuss, then explain why this happens and what the solution is in the docs.

Thanks again for reporting this @letiziasignorelli - I had no idea about this until this discussion!

[yger]

Indeed, for low pass filtering, you need to considerably increase the filter parameter margin_ms that should be used to avoid border effect while filtering. In fact, the problem is that for low pass filtering, this param needs to be increase so much that it makes the chunking rather useless. Each time you'll load a chunk, you'll end up loading much more data before and after to properly estimate the filtering..
For low pass filtering, you can try GaussianFilterRecording, that will perform some sort of low-pass filtering and adjust margin_ms automatically to avoid border effects