PROBLEM #1 – $\delta$-insensitive TDM¶
The stimulus of a single-input single-output (SISO) TEM is modeled as
a bandlimited function of the form
u(t)= \sum_{k=1}^{15} u(kT) \frac{\sin \Omega (t-kT)}{\Omega (t-kT)},
where $\Omega = 2 \pi \cdot 25$ Hz and $T=\frac{\pi}{\Omega}$.
Assume that the TEM describes an ASDM. Generate the trigger times in the time interval $[-3T, 15T]$. Parameters:
$C=0.3, \delta=5e-3$, and choose the injected current to be $b = (C*\delta*\Omega/\pi + 1)$ (it is instructive think why this choice of $b$ makes sense).
Implement a threshold insensitive ($\delta$-insensitive) decoding algorithm for ASDM.
Plot the recovery error (difference between the input stimulus and the recovered waveform) and calculate the Signal-to-Noise ratio of the recovery in decibel (dB).
Compare the $\delta$-insensitive recovery result with the $\delta$-sensitive result (implementation of $\delta$-sensitive algorithm can be found in lecture notebook).
The derivation of $\delta$-insensitive TDM decoding is available in Chapter 6. Highly recommend reading through the entire chapter thoroughly as it will massively help with the following chapters and lectures.
Whole-signal signal-to-Noise ratio (SNR) of the recovery can be calculated for signal $u(t)$ and recovered signal $u_{rec}(t)$ as
SNR = 10\log_{10}\left(\frac{mean(u^2)}{mean((u-u_{rec})^2)}\right)
the result will be in dB
Problem 2 – Derivation of TEM/TDM algorithm for HH with dendritic processing¶
In this problem you are asked to derive and implement a TEM/TDM algorithm for the Hodgkin-Huxley neuron (BSG in the figure below) with dendritic processing modeled as a linear filter ($h(t)$ in the figure below).
The encoding circuit consisting of a filter in cascade with a HH neuron:
In the figure above, the external current injected into the neuron is given as
I_{ext}(t) = b + \underbrace{\left( u \ast h \right)(t)}_{v(t)}
where $b$ is the injected current (also written as $I$ in other texts).
Generate the impulse response of $h(t)$ and visualize.
With $b=20$, encode a randomly generated input stimulus (from Problem 1) using a reduced PIF neuron that is equivalent to the HH neuron model.
Derive an algorithm to recover the signal $u(t)$ from the recieved spikes. Writing down the $t$-transform of the encoding circuit shown above in an inner product form.
Find the time decoding machine (TDM) that recovers the signal $u$. Particularly, provide forms for $q_k$ and $[G]_{lk}$. Please write down the important procedures.
Recover the signal $u(t)$ from output spike times of the reduced PIF and show encoding error and SNR (as in Problem #1 above)
Stimulus and Filter¶
Use the stimulus $u(t)$ from Problem #1, and the filter $h$ is
h(t)= 3 \cdot 120 \cdot \mbox{exp}\left(-100 t\right)\left(\frac{(150 t)^2}{2!}-\frac{(150t)^4}{4!}\right) \cdot \mathbb{1}_{t\ge 0}
note that $\mathbb{1}_{t\ge 0}$ is also known as the Heaviside Step function which ensures that the filter $h(t)$ is causal.
You know the filter $h$ and filtered output $v(t) = (u \ast h)(t)$, but you do not know $u(t)$. You can read the spike times, and you want to recover $u(t)$ from the spikes.
Problem 1¶
%matplotlib inline
import numpy as np
np.random.seed(0) # fix random seed
import matplotlib.pyplot as plt
import typing as tp
from scipy import signal
from scipy.linalg import circulant
from scipy.integrate import cumulative_trapezoid as cumtrapz
from compneuro.utils.phase_response import PIF, iPRC
from compneuro.utils.neuron import limit_cycle
from compneuro.utils.signal import spike_detect, spike_detect_local
from compneuro.neurons.hodgkin_huxley import HodgkinHuxley
from compneuro.utils.signal import convolve
plt.rcParams[“figure.figsize”] = [10, 5]
plt.rcParams[“figure.dpi”] = 100
# TODO: the stimulus ut
def ut(t,T, omega):
return np.nan
Part 1. Implement a threshold insensitive ($\delta$-insensitive) decoding algorithm for ASDM.¶
# TODO: implement ASDM
def integrate_step_asdm():
return np.nan
def asdm_encode():
return np.nan
# TODO: endode with ASDM
# T = np.pi/omega
# t = np.arange(-3*T, 15*T, dt)
# u = ut(t,T,omega)
Disclaimer implementation of $\delta$-sensitive algorithm can be found in lecture notebook.
# TODO: Recovery
# functions for computing the matrices, you can (and probably should) add more
def compute_G(t, tk):
# return G
def compute_q(t, tk):
# return q
def asdm_decoder_sen(t, tk):
# TODO: find the sensitive case
return np.nan
def asdm_decoder_insen(t, tk):
q = compute_q(t, tk)
G = compute_G(t, tk)
# TODO: implement the insensitive case
return np.nan
c_insens, u_rec_insens = np.NaN # delta-insensitive recovery coefficients
# TODO: Plot recovery (in the same plot with the original)
Part 2: Plot the recovery SNR¶
# Compute SNR
# the SNR formula given in the question prompt is for whole-signal SNR
# which is a single value
# for SNR vs t plot use the function one below
def SNR_f(u, u_rec):
return 10 * np.log10(u**2 / (u-u_rec)**2)
# TODO: plot the error and SNR as function of time
Part 3: Compare the $\delta$-insensitive recovery result with the $\delta$-sensitive result¶
# TODO: implement delta-sensitive case
c_sens, u_rec_sens = np.NaN # delta-sensitive recovery coefficients
# TODO: Compare between sensitive and insensitive,
# plot both outputs in same plot with original
# (make use of legends and labels, feel free to adjust linewidth)
# plot both SNR in the same plot
Problem 2¶
Part 1. Generate and show the inpulse response of the filter¶
HINT: For this problem, pay attention to singal time vector, filter time vector and their convolution output lengths.
t_filt = np.arange(T_1, T_2, dt)
* np.exp(-a * t_filt)
(a * t_filt) ** 2 / np.math.factorial(2)
– (a * t_filt) ** 4 / np.math.factorial(4)
) # dendritic filter
# TODO: plot the inpulse response
Part 2. Encode a randomly generated input stimulus (from Problem 1) using a reduced PIF neuron that is equivalent to the HH neuron model.¶
# TODO: filter signal
v = np.NaN
# TODO: compute iPRC at given bias
hh = HodgkinHuxley()
# Plot the iPRC
# TODO: Compute PIF from iPRC
pif_spike_time = np.nan
Part 3. Recover the signal $u(t)$ from the recieved spikes.¶
TODO: A. Write down the t-transform of the encoding circuit in an inner product form
TODO: B. Write down $q_k$ and $[G]_{lk}$, as well the important procedures/steps to obtain $u(t)$ given $h$ (in equation form would suffice)
Part 4. Recover the signal $u(t)$ from PIF spike times and show the statistics¶
# TODO:implement PIF decoder
# def compute_G():
# def compute_q():
# TODO: recover u(t)
c_rec, u_rec = np.nan
# TODO: plot recovery result (in same plot as original), and SNR (as fn of time)
# The recovery may be bad at the boundary because all signals are finite,
# you can focus on the middle part of the signal for comparison.