release 0.9.0: waveform wiggle and speed up sync

release_notes.md

@@ -1,3 +1,8 @@
+# 0.9.0
+## 0.9.0 2024-01-17
+-   `neurodsp.utils.sync_timestamps`: uses FFT based correlation to speed up large arrays alignments
+-   `waveforms`: new wiggle plot for multi-trace waveforms
+
 # 0.8.0
 ## 0.8.1 2023-09-21
 - revert to reading channel info from private methods in shank splitting NP2.4 code to get the original channel layout from shank metadata file

setup.py

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

src/neurodsp/utils.py

@@ -2,7 +2,7 @@
 Window generator, front detections, rms
 """
 import numpy as np
-from scipy import interpolate
+import scipy
 
 
 def sync_timestamps(tsa, tsb, tbin=0.1, return_indices=False):
@@ -23,7 +23,7 @@ def sync_timestamps(tsa, tsb, tbin=0.1, return_indices=False):
     def _interp_fcn(tsa, tsb, ib):
         # now compute the bpod/fpga drift and precise time shift
         drift_ppm = np.polyfit(tsa[ib >= 0], tsb[ib[ib >= 0]] - tsa[ib >= 0], 1)[0] * 1e6
-        fcn_a2b = interpolate.interp1d(tsa[ib >= 0], tsb[ib[ib >= 0]], fill_value="extrapolate")
+        fcn_a2b = scipy.interpolate.interp1d(tsa[ib >= 0], tsb[ib[ib >= 0]], fill_value="extrapolate")
         return fcn_a2b, drift_ppm
 
     # assert sorted inputs
@@ -34,8 +34,7 @@ def _interp_fcn(tsa, tsb, ib):
     y = np.zeros_like(x)
     x[np.int32(np.floor((tsa - tmin) / tbin))] = 1
     y[np.int32(np.floor((tsb - tmin) / tbin))] = 1
-    delta_t = (parabolic_max(np.correlate(x, y, mode='full'))[0] - x.shape[0] + 1) * tbin
-
+    delta_t = (parabolic_max(scipy.signal.correlate(x, y, mode='full'))[0] - x.shape[0] + 1) * tbin
     # do a first assignment at a DT threshold
     ib = np.zeros(tsa.shape, dtype=np.int32) - 1
     threshold = tbin

src/neurodsp/waveforms.py

@@ -5,6 +5,7 @@
 """
 import numpy as np
 import pandas as pd
+import matplotlib.pyplot as plt
 
 
 def _validate_arr_in(arr_in):
@@ -215,7 +216,57 @@ def find_tip_trough(arr_peak, arr_peak_real, df):
     return df, arr_peak
 
 
-def plot_peaktiptrough(df, arr, ax, nth_wav=0, plot_grey=True):
+def plot_wiggle(wav, ax=None, scalar=0.3, clip=1.5, **axkwargs):
+    """
+    Displays a multi-trace waveform in a wiggle traces with negative
+    amplitudes filled
+    :param wav: (nchannels, nsamples)
+    :param axkwargs: keyword arguments to feed to ax.set()
+    :return:
+    """
+    if ax is None:
+        fig, ax = plt.subplots()
+    nc, ns = wav.shape
+    vals = np.c_[wav, wav[:, :1] * np.nan].ravel()  # flat view of the 2d array.
+    vect = np.arange(vals.size).astype(
+        np.float32)  # flat index array, for correctly locating zero crossings in the flat view
+    crossing = np.where(np.diff(np.signbit(vals)))[0]  # index before zero crossing
+    # use linear interpolation to find the zero crossing, i.e. y = mx + c.
+    x1 = vals[crossing]
+    x2 = vals[crossing + 1]
+    y1 = vect[crossing]
+    y2 = vect[crossing + 1]
+    m = (y2 - y1) / (x2 - x1)
+    c = y1 - m * x1
+    # tack these values onto the end of the existing data
+    x = np.hstack([vals, np.zeros_like(c)]) * scalar
+    x = np.maximum(np.minimum(x, clip), -clip)
+    y = np.hstack([vect, c])
+    # resort the data
+    order = np.argsort(y)
+    # shift from amplitudes to plotting coordinates
+    x_shift, y = y[order].__divmod__(ns + 1)
+    ax.plot(y, x[order] + x_shift + 1, 'k', linewidth=.5)
+    x[x > 0] = np.nan
+    x = x[order] + x_shift + 1
+    ax.fill(y, x, 'k', aa=True)
+    ax.set(xlim=[0, ns], ylim=[0, nc], xlabel='sample', ylabel='trace')
+    plt.tight_layout()
+    return ax
+
+
+def plot_peaktiptrough(df, arr, ax=None, nth_wav=0, plot_grey=True):
+    """
+    Plots the peak, tip and trough of a waveform to check the feature extraction
+    :param df:
+    :param arr:
+    :param ax:
+    :param nth_wav:
+    :param plot_grey:
+    :return:
+    """
+    if ax is None:
+        fig, ax = plt.subplots()
     if plot_grey:
         ax.plot(arr[nth_wav], c='gray', alpha=0.5)
     # Peak channel