B(t) = (cos(log[t+tau] - phi) + 1) / 2,  if phi-pi < log[t+tau] < phi+pi

import numpy as np
import matplotlib.pyplot as plt

#%% Helper functions
def get_temporal_kernel(t, kpeaks, delay, weights):
    '''
    Inputs:
        t: [length,1] array
        kpeaks: [1,2] array
        delay: [1,2] array
        weights: [2,] array
    Output:
        kernel: [length,] array
    '''
    log_t = np.log(t+1.3-delay) - np.log(kpeaks)  # Use 1.3 not 0.3 here to compensate for the extra 1-point offset in line 65 of `fitfuns.py`.
    log_t_pi = np.clip(np.pi*log_t/get_temporal_kernel.db2, -np.pi, np.pi)
    log_t_pi[np.isnan(log_t_pi)] = -np.pi
    basis = (np.cos(log_t_pi) + 1) / 2.0
    basis_norm = basis / np.linalg.norm(basis, ord=2, axis=0)
    kernel = basis_norm @ weights
    return kernel
get_temporal_kernel.db2 = 2.0 * np.diff(np.log([100.3, 200.3]))  # 2*db in `fitfuns.py`

#%% Plot kernels with different parameters
weights = np.array([[30.0, -20.0], [30.0, -1.0], [15.0, -20.0]])
kpeaks = np.array([[3.0, 5.0], [3.0, 30.0], [20.0, 40.0]])
delays = np.array([[.0, 0.0], [0.0, 60.0], [20.0, 60.0]])

t = np.expand_dims(np.arange(0,150,1), axis=1)  # milliseconds

fig, axes = plt.subplots(3, 3, figsize=(10, 7))
ri = ci = 0

for ci in range(weights.shape[0]):
    kernel = get_temporal_kernel(t, np.array([[9.67, 20.03]]), np.array([[0.0, 1.0]]), weights[ci,:])
    idx = np.abs(kernel) > 0.0
    axes[ri, ci].plot(-t[idx]/1000, kernel[idx])
    axes[ri, ci].set_ylim([-3.5, 10.0])
    axes[ri, ci].text(0.05, 0.90, 'weights={}'.format(weights[ci,:]), horizontalalignment='left', verticalalignment='top', transform=axes[ri, ci].transAxes)
axes[0, 0].set_ylabel('effect of weights')
ri += 1

# kpeaks parameters controll the spread of both peaks, the second peak must have a bigger spread
for ci in range(kpeaks.shape[0]):
    kernel = get_temporal_kernel(t, kpeaks[[ci],:], np.array([[0.0, 1.0]]), np.array([30.0, -20.0]))
    idx = np.abs(kernel) > 0.0
    axes[ri, ci].plot(-t[idx]/1000, kernel[idx])
    axes[ri, ci].set_xlim([-0.15, 0.005])
    axes[ri, ci].text(0.05, 0.90, 'kpeaks={}'.format(kpeaks[ci,:]), horizontalalignment='left', verticalalignment='top', transform=axes[ri, ci].transAxes)
axes[1, 0].set_ylabel('effects of kpeaks')
ri += 1

for ci in range(delays.shape[0]):
    kernel = get_temporal_kernel(t, np.array([[9.67, 20.03]]), delays[[ci],:], np.array([30.0, -20.0]))
    idx = np.abs(kernel) > 0.0
    axes[ri, ci].plot(-t[idx]/1000, kernel[idx])
    axes[ri, ci].set_xlim([-0.125, 0.001])
    axes[ri, ci].text(0.05, 0.90, 'delays={}'.format(delays[ci,:]), horizontalalignment='left', verticalalignment='top', transform=axes[ri, ci].transAxes)
axes[2, 0].set_ylabel('effects of delays')

fig.show()