add bad channel plot helper and default LF parameters

release_notes.md

@@ -1,8 +1,12 @@
 # 1.4
 
 ## 1.4.0 2024-10-05
-- Optimization of the waveform extractor
-- Refactoring ot the waveform loader with back compability
+- Waveform extraction:
+  - Optimization of the waveform extractor, outputs flattened waveforms
+  - Refactoring ot the waveform loader with back compability
+- Bad channel detector:
+  - The bad channel detector has a plot option to visualize the bad channels and thresholds
+  - The default low-cut filters are set to 300Hz for AP band and 2 Hz for LF band
 
 # 1.3
 

setup.py

@@ -8,7 +8,7 @@
 
 setuptools.setup(
     name="ibl-neuropixel",
-    version="1.3.2",
+    version="1.4.0",
     author="The International Brain Laboratory",
     description="Collection of tools for Neuropixel 1.0 and 2.0 probes data",
     long_description=long_description,

src/ibldsp/plots.py

@@ -1,6 +1,5 @@
 import numpy as np
 import matplotlib.pyplot as plt
-import scipy.signal
 
 
 def show_channels_labels(raw, fs, channel_labels, xfeats, similarity_threshold, psd_hf_threshold=0.02):
@@ -11,9 +10,7 @@ def show_channels_labels(raw, fs, channel_labels, xfeats, similarity_threshold,
     """
     nc, ns = raw.shape
     ns_plot = np.minimum(ns, 3000)
-    vaxis_uv = 75
-    sos_hp = scipy.signal.butter(**{"N": 3, "Wn": 300 / fs * 2, "btype": "highpass"}, output="sos")
-    butt = scipy.signal.sosfiltfilt(sos_hp, raw)
+    vaxis_uv = 250 if fs < 2600 else 75
     fig, ax = plt.subplots(1, 5, figsize=(18, 6), gridspec_kw={'width_ratios': [1, 1, 1, 8, .2]})
     ax[0].plot(xfeats['xcor_hf'], np.arange(nc))
     ax[0].plot(xfeats['xcor_hf'][(iko := channel_labels == 1)], np.arange(nc)[iko], 'k*')
@@ -30,7 +27,7 @@ def show_channels_labels(raw, fs, channel_labels, xfeats, similarity_threshold,
     ax[2].plot([-.75, -.75], [0, nc], 'y--')
     ax[2].set(yticklabels=[], xlabel='LF coherence', ylim=[0, nc], title='c) outside')
     ax[2].sharey(ax[0])
-    im = ax[3].imshow(butt[:, :ns_plot] * 1e6, origin='lower', cmap='PuOr', aspect='auto',
+    im = ax[3].imshow(raw[:, :ns_plot] * 1e6, origin='lower', cmap='PuOr', aspect='auto',
                       vmin=-vaxis_uv, vmax=vaxis_uv, extent=[0, ns_plot / fs * 1e3, 0, nc])
     ax[3].set(yticklabels=[], title='d) Raw data', xlabel='time (ms)', ylim=[0, nc])
     ax[3].grid(False)

src/ibldsp/voltage.py

@@ -713,15 +713,18 @@ def nxcor(x, ref):
     raw = raw - np.mean(raw, axis=-1)[:, np.newaxis]  # removes DC offset
     xcor = channels_similarity(raw)
     fscale, psd = scipy.signal.welch(raw * 1e6, fs=fs)  # units; uV ** 2 / Hz
-    if psd_hf_threshold is None:
-        # the LFP band data is obviously much stronger so auto-adjust the default threshold
-        psd_hf_threshold = 1.4 if fs < 5000 else 0.02
-    sos_hp = scipy.signal.butter(
-        **{"N": 3, "Wn": 300 / fs * 2, "btype": "highpass"}, output="sos"
-    )
+    # auto-detection of the band with which we are working
+    band = 'ap' if fs > 2600 else 'lf'
+    # the LFP band data is obviously much stronger so auto-adjust the default threshold
+    if band == 'ap':
+        psd_hf_threshold = 0.02 if psd_hf_threshold is None else psd_hf_threshold
+        filter_kwargs = {"N": 3, "Wn": 300 / fs * 2, "btype": "highpass"}
+    elif band == 'lf':
+        psd_hf_threshold = 1.4 if psd_hf_threshold is None else psd_hf_threshold
+        filter_kwargs = {"N": 3, "Wn": 1 / fs * 2, "btype": "highpass"}
+    sos_hp = scipy.signal.butter(**filter_kwargs, output="sos")
     hf = scipy.signal.sosfiltfilt(sos_hp, raw)
     xcorf = channels_similarity(hf)
-
     xfeats = {
         "ind": np.arange(nc),
         "rms_raw": utils.rms(raw),  # very similar to the rms avfter butterworth filter