diff --git a/code/chirpdetection.py b/code/chirpdetection.py index c7e1be8..c4a04e7 100755 --- a/code/chirpdetection.py +++ b/code/chirpdetection.py @@ -589,16 +589,16 @@ def main(datapath: str, plot: str) -> None: raw_time = np.arange(data.raw.shape[0]) / data.raw_rate # good chirp times for data: 2022-06-02-10_00 - window_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate - window_duration_index = 60 * data.raw_rate + # window_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate + # window_duration_index = 60 * data.raw_rate # t0 = 0 # dt = data.raw.shape[0] # window_start_seconds = (23495 + ((28336-23495)/3)) * data.raw_rate # window_duration_seconds = (28336 - 23495) * data.raw_rate - # window_start_index = 0 - # window_duration_index = data.raw.shape[0] + window_start_index = 0 + window_duration_index = data.raw.shape[0] # generate starting points of rolling window window_start_indices = np.arange( diff --git a/code/modules/datahandling.py b/code/modules/datahandling.py index 9460b5b..a1e9f18 100644 --- a/code/modules/datahandling.py +++ b/code/modules/datahandling.py @@ -1,9 +1,10 @@ import numpy as np from typing import List, Any from scipy.ndimage import gaussian_filter1d +from scipy.stats import gamma, norm -def norm(data): +def scale01(data): """ Normalize data to [0, 1] @@ -209,6 +210,117 @@ def flatten(list: List[List[Any]]) -> List: return [item for sublist in list for item in sublist] +def causal_kde1d(spikes, time, width, shape=2): + """ + causalkde computes a kernel density estimate using a causal kernel (i.e. exponential or gamma distribution). + A shape of 1 turns the gamma distribution into an exponential. + + Parameters + ---------- + spikes : array-like + spike times + time : array-like + sampling time + width : float + kernel width + shape : int, optional + shape of gamma distribution, by default 1 + + Returns + ------- + rate : array-like + instantaneous firing rate + """ + + # compute dt + dt = time[1] - time[0] + + # time on which to compute kernel: + tmax = 10 * width + + # kernel not wider than time + if 2 * tmax > time[-1] - time[0]: + tmax = 0.5 * (time[-1] - time[0]) + + # kernel time + ktime = np.arange(-tmax, tmax, dt) + + # gamma kernel centered in ktime: + kernel = gamma.pdf( + x=ktime, + a=shape, + loc=0, + scale=width, + ) + + # indices of spikes in time array: + indices = np.asarray((spikes - time[0]) / dt, dtype=int) + + # binary spike train: + brate = np.zeros(len(time)) + brate[indices[(indices >= 0) & (indices < len(time))]] = 1.0 + + # convolution with kernel: + rate = np.convolve(brate, kernel, mode="same") + + return rate + + +def acausal_kde1d(spikes, time, width): + """ + causalkde computes a kernel density estimate using a causal kernel (i.e. exponential or gamma distribution). + A shape of 1 turns the gamma distribution into an exponential. + + Parameters + ---------- + spikes : array-like + spike times + time : array-like + sampling time + width : float + kernel width + shape : int, optional + shape of gamma distribution, by default 1 + + Returns + ------- + rate : array-like + instantaneous firing rate + """ + + # compute dt + dt = time[1] - time[0] + + # time on which to compute kernel: + tmax = 10 * width + + # kernel not wider than time + if 2 * tmax > time[-1] - time[0]: + tmax = 0.5 * (time[-1] - time[0]) + + # kernel time + ktime = np.arange(-tmax, tmax, dt) + + # gamma kernel centered in ktime: + kernel = norm.pdf( + x=ktime, + loc=0, + scale=width, + ) + + # indices of spikes in time array: + indices = np.asarray((spikes - time[0]) / dt, dtype=int) + + # binary spike train: + brate = np.zeros(len(time)) + brate[indices[(indices >= 0) & (indices < len(time))]] = 1.0 + + # convolution with kernel: + rate = np.convolve(brate, kernel, mode="same") + + return rate + + if __name__ == "__main__": timestamps = [