diff --git a/nwbwidgets/analysis/spikes.py b/nwbwidgets/analysis/spikes.py index e7c2c471..7a13a858 100644 --- a/nwbwidgets/analysis/spikes.py +++ b/nwbwidgets/analysis/spikes.py @@ -14,10 +14,16 @@ def compute_smoothed_firing_rate(spike_times, tt, sigma_in_secs): Returns: Gaussian smoothing evaluated at array t """ - if len(spike_times) < 2: + if len(spike_times) < 1: return np.zeros_like(tt) binned_spikes = np.zeros_like(tt) - binned_spikes[np.searchsorted(tt, spike_times)] += 1 + # find bin indices for spike times + spike_idx = np.searchsorted(tt, spike_times, side="right") - 1 + # filter out negative spike idxs, though there shouldn't be any + spike_idx = spike_idx[spike_idx >= 0].astype(int) + # increment binned spike count at bin indices + np.add.at(binned_spikes, spike_idx, 1) + # smooth data dt = np.diff(tt[:2])[0] sigma_in_samps = sigma_in_secs / dt smooth_fr = scipy.ndimage.gaussian_filter1d(binned_spikes, sigma_in_samps) / dt diff --git a/test/test_analysis_spikes.py b/test/test_analysis_spikes.py new file mode 100644 index 00000000..082b79a7 --- /dev/null +++ b/test/test_analysis_spikes.py @@ -0,0 +1,27 @@ +import numpy as np + +from nwbwidgets.analysis.spikes import compute_smoothed_firing_rate + + +def test_compute_smoothed_firing_rate(): + spike_times = np.array([1.0, 2.0, 5.0, 5.5, 7.0, 7.5, 8.0]) + tt = np.arange(10, dtype=float) + expected_binned_spikes = np.array([0.0, 1.0, 1.0, 0.0, 0.0, 2.0, 0.0, 2.0, 1.0, 0.0]) + expected_smoothed_spikes = np.array( + [ + 0.35438556, + 0.64574827, + 0.64991247, + 0.40421249, + 0.55136343, + 0.91486675, + 1.02201074, + 1.14797447, + 0.89644961, + 0.41307621, + ] + ) + binned_spikes = compute_smoothed_firing_rate(spike_times, tt, 0.001) + smoothed_spikes = compute_smoothed_firing_rate(spike_times, tt, 1.0) + np.testing.assert_allclose(expected_binned_spikes, binned_spikes) + np.testing.assert_allclose(expected_smoothed_spikes, smoothed_spikes)