<jupyter_start><jupyter_text>Parameter analysis: kernel shapeWe analyse how the parameters $\tau_{{\sf half}}$ and $\sigma_{{\sf gauss}}$ affect the shape of the kernel $F = G \ast E$.We use the original processing code and empirically obtain impulse responses to illustrate the shape of the kernels<jupyter_code>import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline
import pandas as pd
mpl.rcParams['savefig.dpi'] = 150 # for print, go to 600
def compute_sig_smooth_diff_exp(signal, sigma_gauss, tau_half):
    """
    Smooth the signal with gaussian kernel, form the discrete differential, and
    convolve it with an exponential kernel.
    
    sigma_gauss - width of the gaussian kernel
    tau_half - half-life of the exponential kernel
    """
    windowsize = int(10 * sigma_gauss) # 5 sigma windowsize
    sig_smooth = pd.Series(signal).rolling(windowsize,  win_type='gaussian').mean(std=sigma_gauss)
    sig_diff = np.diff(sig_smooth)
    sig_diff_ewma = pd.Series(sig_diff).ewm(halflife=tau_half, adjust=False, ignore_na=True).mean()
    return sig_diff_ewma

def compute_sig_smooth_exp(signal, sigma_gauss, tau_half):
    """
    Smooth the signal with gaussian kernel, and convolve it with an exponential kernel 
    (no differential formed).
    sigma_gauss - width of the gaussian kernel
    tau_half - half-life of the exponential kernel
    """
    windowsize = int(10 * sigma_gauss) # 5 sigma windowsize
    sig_smooth = pd.Series(signal).rolling(windowsize, win_type='gaussian').mean(std=sigma_gauss)
    sig_ewma = pd.ewma(sig_smooth, halflife=tau_half, adjust=False, ignore_na=True)
    return sig_ewma


#construct a signal of length 2000, equivalent to 20 seconds at a sampling rate of 100 Hz as in the real data.

sigma_gauss = 30 # 0.3 s
tau_half = 40 # 0.4 s
siglength = 2000 # 20 s

signal = np.zeros(siglength)
signal[500] = 1. # impulse at 5.0 s<jupyter_output><empty_output><jupyter_text>Impulse response and spectral analysis of the kernels<jupyter_code>hls = np.arange(5.,51.,20.)
std = 25.

colors = cm.gray(np.linspace(0,1,4))

f = plt.figure(figsize=(8,5))
gs = mpl.gridspec.GridSpec(2,2, hspace=0.6, wspace=0.4)
axa = f.add_subplot(gs[0,0])
axb = f.add_subplot(gs[0,1])
for i,h in enumerate(hls):
    kernel_gde = compute_sig_smooth_diff_exp(signal, sigma_gauss=std, tau_half=h)
    kernel_gde[np.isnan(kernel_gde)] = 0. # replace nans with 0
    kernel_gde_norm = kernel_gde / np.nanmax(kernel_gde)
    xaxis = np.arange(len(kernel_gde_norm))/100.
    axa.plot(xaxis, kernel_gde_norm, color=colors[i], label="{:.2f}".format(h/100.))
    spec = np.fft.fft(kernel_gde_norm)
    xaxis_spec = np.fft.fftfreq(len(kernel_gde_norm), 0.01)
    axb.semilogx(xaxis_spec[:int(siglength/2.-1)], np.abs(spec)[:int(siglength/2.-1)],
                 color=colors[i], label="{:.2f}".format(h/100.))
axa.set_frame_on(False)
axa.xaxis.set_ticks_position('bottom')
axa.yaxis.set_ticks_position('left')
axa.set_xlim(5,9)
axa.set_ylim(-1.05, 1.05)
axa.set_xlabel('time [s]')
axa.set_title(r"$\sigma_{{\sf smooth}}$=  {:.2f} s".format(std/100.))
axa.legend(frameon=False, title=r"$\tau_{{\sf half}}$ [s]", fontsize=8, loc=(0.8, 0.7))
axa.spines['left'].set_position(('outward',6))
axa.spines['bottom'].set_position(('outward',5))
f.text(0.07, 0.93, "A)", weight='bold', fontsize=12)

axb.set_xlim(0,4.)
axb.set_xlabel('f [Hz]')
axb.set_ylabel("rel. power [a.u]")
miny = axb.get_yticks()[0]
maxy = axb.get_yticks()[-1]
axb.set_yticks([miny, maxy/2, maxy])
axb.set_frame_on(False)
axb.xaxis.set_ticks_position('bottom')
axb.yaxis.set_ticks_position('left')
axb.legend(frameon=False, title=r"$\tau_{{\sf half}}$ [s]", fontsize=8, loc=(0.8, 0.7))
axb.spines['left'].set_position(('outward',6))
axb.spines['bottom'].set_position(('outward',5))
f.text(0.07, 0.45, "C)", weight='bold', fontsize=12)

stds = np.arange(5.,51.,20.)
hl = 25.
axc = f.add_subplot(gs[1,0])
axd = f.add_subplot(gs[1,1])
for i,s in enumerate(stds):
    kernel_gde = compute_sig_smooth_diff_exp(signal, sigma_gauss=s, tau_half=hl)
    kernel_gde[np.isnan(kernel_gde)] = 0. # replace nans with 0
    kernel_gde_norm = kernel_gde / np.nanmax(kernel_gde)
    xaxis = np.arange(len(kernel_gde_norm))/100.
    axc.plot(xaxis, kernel_gde_norm, color=colors[i], label="{:.2f}".format(s/100.))
    spec = np.fft.fft(kernel_gde_norm)
    xaxis_spec = np.fft.fftfreq(len(kernel_gde_norm), 0.01)
    axd.semilogx(xaxis_spec[:int(siglength/2.-1)], np.abs(spec)[:int(siglength/2.-1)],
                 color=colors[i], label="{:.2f}".format(s/100.))
axc.set_frame_on(False)
axc.xaxis.set_ticks_position('bottom')
axc.yaxis.set_ticks_position('left')
axc.set_xlim(5,9)
axc.set_ylim(-1.05, 1.05)
axc.set_xlabel('time [s]')
axc.set_title(r"$\tau_{{\sf half}}$=  {:.2f} s".format(hl/100.))
axc.legend(frameon=False, title=r"$\sigma_{{\sf smooth}}$ [s]", fontsize=8, loc=(0.8, 0.7))
axc.spines['left'].set_position(('outward',6))
axc.spines['bottom'].set_position(('outward',5))
f.text(0.50, 0.93, "B)", weight='bold', fontsize=12)

axd.set_xlim(0,4.)
axd.set_xlabel('f [Hz]')
axd.set_ylabel("rel. power [a.u]")
maxy = axd.get_yticks()[-1]
axd.set_yticks([miny, maxy/2, maxy])
axd.set_frame_on(False)
axd.xaxis.set_ticks_position('bottom')
axd.yaxis.set_ticks_position('left')
axd.set_title(r" ")
axd.legend(frameon=False, title=r"$\sigma_{{\sf smooth}}$ [s]", fontsize=8, loc=(0.8, 0.7))
axd.spines['left'].set_position(('outward',6))
axd.spines['bottom'].set_position(('outward',5))
f.text(0.50, 0.45, "D)", weight='bold', fontsize=12)

f.savefig("Figures/Fig. 11 - Kernel shapes.png", dpi=600)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup