# Bandlimited signal reconstruction from
leaky integrate-and-fire encoding using POCS

Nguyen T. Thao, Dominik Rzepka and Marek Miśkowicz N. T. Thao is with the
Department of Electrical Engineering, The City College of New York, CUNY, New
York, USA, email<EMAIL_ADDRESS>Rzepka and M. Miśkowicz are with
the Department of Measurement and Electronics, AGH University of Science and
Technology, Kraków, Poland, emails<EMAIL_ADDRESS><EMAIL_ADDRESS>Rzepka and M. Miśkowicz were supported by the Polish National Center of
Science under grants DEC-2017/27/B/ST7/03082 and DEC-2018/31/B/ST7/03874,
respectively.

###### Abstract

Leaky integrate-and-fire (LIF) encoding is a model of neuron transfer function
in biology that has recently attracted the attention of the signal processing
and neuromorphic computing communities as a technique of event-based sampling
for data acquisition. While LIF enables the implementation of analog-circuit
signal samplers of lower complexity and higher accuracy simultaneously, the
core difficulty of this technique is the retrieval of an input from its LIF-
encoded output. This theoretically requires to perform the pseudo-inversion of
a linear but time-varying operator of virtually infinite size. In the context
of bandlimited inputs to allow finite-rate sampling, we show two fundamental
contributions of the method of projection onto convex sets (POCS) to this
problem: (i) single iterations of the POCS method can be used to
deterministically improve input estimates from any other reconstruction
method; (ii) the iteration limit of the POCS method is the pseudo-inverse of
the above mentioned operator in all conditions, whether reconstruction is
unique or not and whether the encoding is corrupted by noise or not. The
algorithms available until now converge only under particular situations of
unique of reconstruction.

###### Index Terms:

integrate and fire, leakage, bandlimited signals, nonuniform sampling, event-
based sampling, time-encoding machine, time quantization, pseudo-inverse,
POCS.

##  I Introduction

Integrate-and-fire (IF) encoding is a biologically-inspired model for mapping
a continuous-time stimulus to a spike train. As opposed to traditional pulse-
code modulation (PCM) in data conversion, which outputs a sequence of
amplitude values, the input information that an IF encoder provides is in the
timing of its output spikes. It specifically fires an impulse whenever the
integral of the stimulus reaches a given threshold (see Fig. 1). Time encoding
is attracting more and more interest in data acquisition [1, 2, 3, 4, 5, 6, 7,
8, 7, 9, 10] as the downscaling of semiconductor integration is increasing
time precision while resulting in less amplitude accuracy [11, 12, 13]. This
is also part of the trend on event-based sampling [14] with more general
objectives such as having acquisition activity dependent on input activity for
power efficiency. What has made time encoding a difficult orientation is the
non-trivial digital postprocessing it requires to convert the time information
back to an explicit description of the input. This direction of research took
off with the pioneering work of [15] which introduced an algorithm of perfect
reconstruction of bandlimited inputs encoded by asynchronous Sigma-Delta
modulators (ASDM).

The purpose of the present paper is to revisit the problem of bandlimited
input reconstruction from IF outputs with two simultaneous goals: (i) approach
IF encoding from the more fundamental perspective of generalized nonuniform
sampling and (ii) propose a technique of input reconstruction of arbitrary
accuracy that is consistent with this approach. We will perform our analysis
on the more general version of leaky integrate-and-fire (LIF). In data
acquisition, leakage is an artifact that is omnipresent in circuit integrators
[16, 17, 18]. Meanwhile, in neuromorphic systems, leakage is a factor that is
included in the IF model of a neuron to reflect the firing dynamics of the
biological nervous systems with higher fidelity [Stein, 1965]. Interest in the
processing of LIF encoded signals is therefore expected to grow with the
development of neuromorphic hardware [19, 20, 21, 22].

Figure 1: Encoding process of LIF (bipolar configuration with $c=0$).

The earliest results of perfect reconstruction of bandlimited inputs from
time-encoding machines (TEM) are found in [15, 23, 24]. TEM is a generic term
that usually refers to the use of either ASDM, IF or LIF. Although TEM’s are
both nonlinear and time-varying, the fundamental fact used by these works is
that a given time-encoded output does give the knowledge of some linear (but
time-varying) transformation of the input, which depends on the
characteristics of the encoder. In this article, we call the function of this
transformation the sampling operator. Thus, signal reconstruction boils down
to inverting this operator. Although the parameters of this operator are
derived in a different manner in ASDM and LIF, the generic form of the
operator is the same in ASDM and IF. Thus, working with LIF covers the most
general case of TEM’s. The idea of inverting the sampling operator was
presented in [24] in the form of pseudo-inversion of a matrix. However, as the
matrix size is determined by the total number of signal samples, the
computational requirements for this pseudo-inversion become prohibitive in
practical applications, and more particularly in streaming signal processing.
In the case of ASDM, a practical method was later introduced in [25] where
this inversion is performed locally in small blocks of input. Splitting
bandlimited signals into blocks however creates ill-posed inversion situations
[26] with little mathematics and analytical control on the block boundary
effects, and the only available methods of block merging are purely empirical.
Meanwhile, the original works of [15, 23] presented more analytically
motivated algorithms of input estimation, based on iterating contractive
mappings. These methods however require certain sufficient (but not necessary)
conditions on the density of the firing instants to theoretically guarantee
uniqueness of reconstruction. This does not cover the more general situation
of pseudo-inversion when reconstruction is not unique. Moreover, we will see
that a strict oversampling rate of firing instants is not sufficient to meet
these required conditions in the presence of leakage.

After giving the exact description of LIF in Section II together with our
signal assumptions and notation, we start our revisit of bandlimited input
reconstruction in LIF in Section III by briefly reviewing the basics and prior
works on one-step input estimation. This includes a deterministic technique
proposed in [27] and the classic statistical approach of neural research [28,
29, 30]. In Section IV, we give the new interpretation of an LIF output as
providing inner-product values of the input with known functions that we call
the kernel sampling functions. We say that an estimate is consistent with the
input when it yields the same inner-product values [31]. We show that any
estimate (including those of Section III) that is not consistent can be
automatically and deterministically improved by an orthogonal projection. By
alternating this projection with bandlimitation, one obtains an iterative way
to improve an estimate following the classic method of projection onto convex
sets (POCS) [32, 33]. This method was previously used to reconstruct
bandlimited inputs from various signal acquisition encoders in [34, 35, 36].
After connecting the sampling operator with the kernel inner-products, we show
in Section V that the POCS iteration actually tends exactly to the pseudo-
inverse solution of this operator (under proper space norm settings and
assuming a zero initial estimate) in all situations, whether reconstruction is
unique or not, and whether the sampling is noisy or not. In the noisy case,
not only does this guarantee the stability of the iteration, but this also
gives a theoretical characterization of the deviation of the POCS limit.
Section V is also the opportunity to see the LIF encoding of bandlimited
inputs as a mere extension of the pioneering work by Yen (1956) [37] on the
nonuniform point sampling of bandlimited signals. In this process, we point
the reader to a unified generalized view on nonuniform sampling. As the POCS
method is theoretically defined on continuous-time signals, we finally present
in Section VI a rigorous discretization of the POCS iteration that does not
involve Nyquist-rate resampling.

##  II Signal and system settings

###  II-A LIF system

For a given input signal $x(t)$, the bipolar version of LIF outputs an impulse
train sequence of the type

$s(t)=\textstyle\sum\limits_{n=1}^{N}\epsilon_{n}\delta(t-t_{n})$

where $(t_{n})_{1\leq n\leq N}$ is an increasing sequence of positive instants
and $(\epsilon_{n})_{1\leq n\leq N}$ is a sequence of signs. These two
sequences are recursively defined by

$\displaystyle
t_{n}=\min\Big{\\{}t>t_{n-1}:\Big{|}\int_{t_{n-1}}^{t}e^{-\alpha(t-s)}\big{(}x(s)+c\big{)}\,{\mathrm{d}}s\Big{|}=\theta\Big{\\}}$
(1)

and

$\int_{t_{n-1}}^{\raisebox{0.90417pt}{$\scriptstyle
t_{n}$}}e^{-\alpha(t_{n}-s)}\big{(}x(s)+c\big{)}\,{\mathrm{d}}s=\epsilon_{n}\theta$
(2)

starting from $t_{0}=0$, where $\alpha\geq 0$, $c\geq 0$ and $\theta>0$ are
known constants. Fig. 1 gives an illustration of this process with $c=0$.
Unipolar LIF is the particular case where $x(t)+c\geq 0$. Since all integral
values are non-negative in this case, the absolute-value function can then be
removed from (1) and $\epsilon_{n}=+1$ for all $n$ in (2).

For convenience, we will use the following notation throughout the paper,

$\displaystyle{\sf N}:=\big{\\{}1,{\cdots},N\\},$ $\displaystyle\forall
t\in{\mathbb{R}},\qquad\displaystyle
A_{\tau}u(t):=\int_{\tau}^{t}e^{-\alpha(t-s)}u(s)\,{\mathrm{d}}s,\quad\tau\in{\mathbb{R}},$
(3) $\displaystyle x^{c}(t):=x(t)+c.$

Then, (1) and (2) respectively take the simple forms of

$\displaystyle
t_{n}=\min\big{\\{}t>t_{n-1}:|A_{t_{n-1}}x^{c}(t)|=\theta\big{\\}},\qquad
n\in{\sf N}$ (4) $A_{t_{n-1}}x^{c}(t_{n})=\epsilon_{n}\theta,\qquad n\in{\sf
N}.$ (5)

Qualitatively speaking, LIF consists in detecting after each firing instant
$t_{n}$ the next crossing of $A_{t_{n}}x^{c}(t)$ with the levels $\pm\theta$.
LIF can thus be seen as a generalized version of level crossing sampling (LCS)
[38, 39, 35], where the thresholded signal is a dynamically changing
integrated version of the input.

###  II-B Signal setting

All signals in this paper are assumed to be in the space $L^{2}(D)$ of real
square-integrable functions on $D$, where $D$ is either ${\mathbb{R}}$ or
$[0,T]$. This space is equipped with the inner-product

$\langle u,v\rangle:=\int_{D}u(t)v(t)\,{\mathrm{d}}t,$ (6)

which induces the norm $\|u\|:=\langle u,u\rangle^{1\\!/2}$. Fourier transform
is defined in both cases of domain $D$, with the difference that the frequency
values are discrete and multiple of $2\pi/T$ when $D=[0,T]$ from the Fourier
series of $T$-periodic signals.

For convenience, we assume that all bandlimited signals have Nyquist period 1
(which can always be achieved by a change of time unit). The baseband is then
comprised of the frequencies $\omega\in[-\pi,\pi]$. We call ${\cal B}$ the
space of such signals. In the finite time case, we assume that $T$ is a
multiple of the Nyquist period, and hence an integer. We denote by
$\widetilde{u}(t)$ the bandlimited version of any signal $u(t)\in L^{2}(D)$.
In the case $D={\mathbb{R}}$, we have

$\widetilde{u}(t)=\varphi(t)*u(t)$

where $\varphi(t)$ is the sinc function $\sin(\pi t)/(\pi t)$ and $*$ is the
convolution operation. In the case $D=[0,T]$, we also have this expression
assuming that $\varphi(t)$ is the periodic sinc function (or Dirichlet kernel)
and $*$ is the circular convolution operation over $[0,T]$ (we assume that $T$
is an odd integer for the proper formation of the Dirichlet kernel).

The bounded domain $D=[0,T]$ is used in practice to designate the time window
of signal acquisition. Given the decay of the sinc function, the $T$-periodic
bandlimited signals asymptotically match the bandlimited signals of
$L^{2}({\mathbb{R}})$ as one looks at the time instants of $[0,T]$ that get
away from the boundaries. But given that $T$ is in practice virtually infinite
compared to the considered Nyquist period, the boundary effects are typically
neglected. Meanwhile, $D=[0,T]$ theoretically makes ${\cal B}$ of finite
dimension (as a matter of fact, equal to $T$ given our assumptions), which
allows us to define the situations of critical sampling and oversampling from
nonuniform samples.

##  III One-step bandlimited estimation

###  III-A Basic linear reconstruction

In one-step bandlimited estimation, one looks at $x^{c}(t)$ as the input to
the LIF encoder (in other words, one thinks of the constant component $c$ as
part of the input). Using the output of the encoder, the goal is then to
construct a bandlimited signal $u(t)$ that minimizes $\|u-x^{c}\|^{2}$, which
we call the mean-squared error (MSE) of $u(t)$. A basic approach is to start
searching for an estimate $u(t)$ that yields the same output as $x^{c}(t)$
through the same LIF encoder. With the sequences $(t_{n})_{n\in{\sf N}}$ and
$(\epsilon_{n})_{n\in{\sf N}}$ obtained from the encoding of $x(t)$, this
amounts to requiring $u(t)$ to satisfy both (4) and (5) recursively. Note that
(4) amounts to inequality constraints while (5) is an equality. The latter
amounts to providing for each $n\in{\sf N}$ a sample of an affine
transformation of $x(t)$. In the framework of generalized nonuniform sampling,
we will simply be interested in estimates $u(t)$ that are consistent with the
equality constraint of (5), i.e., such that

$\forall n\in{\sf N},\qquad A_{t_{n-1}}u(t_{n})=\epsilon_{n}\theta.$ (7)

The crudest way to obtain such an estimate is to take the signal

$u_{0}(t):={\textstyle\sum\limits_{n\in{\sf
N}}}\epsilon_{n}\theta\,\delta(t-t_{n}).$ (8)

Indeed, if one thinks of $\int_{a}^{b}v(s){\mathrm{d}}s$ as the integral of
$v(s)$ from $a^{+}$ to $b^{+}$, then $u_{0}(t)$ is easily seen to satisfy (7).
But as $u_{0}(t)$ is not bandlimited, one adopts as final approximation of
$x(t)$ its bandlimited version

$\widetilde{u}_{0}(t)=\varphi(t)*u_{0}(t)={\textstyle\sum\limits_{n\in{\sf
N}}}\epsilon_{n}\theta\,\varphi(t-t_{n}).$

In this process, $\widetilde{u}_{0}(t)$ is likely to lose consistency. But we
will keep this simple reconstruction as reference.

###  III-B Relation to the prior work of [27]

With $D={\mathbb{R}}$, the work of [27] can be interpreted as starting from
$u_{0}(t)$ as an initial estimate but proposing a different bandlimited
transformation of it as final approximation of $x(t)$. The authors consider
the global operator

$\forall t\in{\mathbb{R}},\qquad
Au(t):=A_{-\infty}u(t)=\int_{-\infty}^{t}e^{-\alpha(t-s)}u(s)\,{\mathrm{d}}s$

to take advantage of the following features:

  * •

$A$ commutes with bandlimitation since it is a convolution operator,

  * •

$A$ is invertible of inverse $A^{-1}v(t)=v^{\prime}(t)+\alpha\,v(t)$,

  * •

an estimate $u(t)$ that is 0 for all $t\leq 0$ can be shown to satisfy (7) if
and only if $Au(t_{n})=Au_{0}(t_{n})$ for all $n\in{\sf N}$.

Then, they propose to estimate $x(t)$ with a signal $\hat{u}_{0}(t)$ such that
$A\hat{u}_{0}(t)$ is a bandlimited approximation of $Au_{0}(t)$. A downside of
their approach, however, is that they specifically take the Nyquist-rate
aliased version of $Au_{0}(t)$

$A\hat{u}_{0}(t)={\textstyle\sum\limits_{k\in{\mathbb{Z}}}}Au_{0}(k)\,\varphi(t-k)$
(9)

under our setting of Nyquist period 1. While this enables the simple
derivation of $\hat{u}_{0}(t)$ as
$\sum_{k\in{\mathbb{Z}}}Au_{0}(k)\,\psi(t-k)$ where
$\psi(t):=A^{-1}\varphi(t)=\varphi^{\prime}(t)+\alpha\,\varphi(t)$, the
aliasing in (9) is expected to be significant as $Au_{0}(t)$ is typically non-
bandlimited given its jumps at the instants $(t_{n})_{n\in{\sf N}}$. In
contrast, $A\widetilde{u}_{0}(t)$ yields the aliasing-free relation

$\displaystyle A\widetilde{u}_{0}(t)$ $\displaystyle=\varphi(t)*Au_{0}(t)$

since $A$ commutes with bandlimitation. Not only is $\widetilde{u}_{0}(t)$
substantially simpler to implement than $\hat{u}_{0}(t)$, but the numerical
experiments of Section III-D also confirm its higher accuracy as an estimate
of $x^{c}(t)$.

###  III-C Optimal data-driven convolutional reconstruction

Note that the bandlimited estimate $\widetilde{u}_{0}(t)$ is of the form

$u(t)=f(t)*s(t)$

where $f(t)$ is specifically chosen to be $\theta\,\varphi(t)$. A basic
direction is to look for better reconstructions in this larger family of
signals. This has been a basic approach of neural research to estimate the
stimulus of a neuron from its output spike train. As the exact transfer
function of a neuron is unknown, this direction of research has been
considered based on statistics of effective input-output pairs of the neuron
[28] [29, §2.3]. This method has been explicitly applied on LIF in [30,
§4.3.3] as demonstration. Assuming that $x(t)$ is a random process and $s(t)$
is the resulting response of the LIF encoder, the goal is to find the function
$f_{\scriptscriptstyle\mathrm{opt}}(t)$ that minimizes

$E\Big{(}\big{\|}f(t)*s(t)-x^{c}(t)\big{\|}^{2}\Big{)}.$ (10)

This is a Wiener filtering problem. Calling $U(\omega)$ the Fourier transform
of $u(t)$, it is known (see for example [40, (7.3.2)]) that

$F_{\scriptscriptstyle\mathrm{opt}}(\omega)=\frac{E\big{(}S(\omega)^{*}X^{c}(\omega)\big{)}}{E\big{(}|S(\omega)|^{2}\big{)}},$
(11)

for all $\omega$ in the baseband, and
$F_{\scriptscriptstyle\mathrm{opt}}(\omega)=0$ otherwise since
$f_{\scriptscriptstyle\mathrm{opt}}(t)$ must obviously be bandlimited. This
leads to the definition of the new estimate

$u_{\scriptscriptstyle\mathrm{opt}}(t):=f_{\scriptscriptstyle\mathrm{opt}}(t)*s(t)={\textstyle\sum\limits_{n\in{\sf
N}}}\epsilon_{n}\,f_{\scriptscriptstyle\mathrm{opt}}(t-t_{n}).$ (12)

The MSE of $u_{\scriptscriptstyle\mathrm{opt}}(t)$ is of course expected to be
smaller than that of $\widetilde{u}_{0}(t)$. This is at the price of having
the input-output statistics of the system. But it will be interesting to know
for reference how much MSE reduction can be achieved with this method, as done
in the next section.

When $c>0$, note that $x^{c}(t)$ cannot be in $L^{2}({\mathbb{R}})$. In this
case, (10) and (11) are at least well defined with $D=[0,T]$. With the
theoretical setting of $D={\mathbb{R}}$, (10) and (11) can still be used
provided that the constant components are removed from the analysis, thus
making $u_{\scriptscriptstyle\mathrm{opt}}(t)$ optimal only up to a constant
component.

###  III-D Numerical experiments

Figure 2: MSE of $\hat{u}_{0}(t)$, $\widetilde{u}_{0}(t)$ and
$u_{\scriptscriptstyle\mathrm{opt}}(t)$ for various leakage values $\alpha$ in
the unipolar and bipolar cases.

We show in Fig. 2 the MSE results we got numerically for the bandlimited
reconstruction estimates $\hat{u}_{0}(t)$, $\widetilde{u}_{0}(t)$ and
$\widetilde{u}_{\scriptscriptstyle\mathrm{opt}}(t)$ with various leakage
values $\alpha$ in the unipolar and bipolar cases. The numerical difficulty is
to eliminate the boundary effects so that the true MSE from steady state can
be extracted. We do this by working with bandlimited inputs that are periodic
over 61 Nyquist periods. This enables the exact implementation of the sinc
filter without boundary deviations. We adjust the threshold $\theta$ so that
the number of firing instants per Nyquist period is in average equal to 1.5.
We call this number the oversampling ratio. The MSE of each estimate is
averaged over 100 zero-mean bandlimited inputs whose Nyquist-rate samples are
drawn randomly and uniformly in $[-0.7,0.7]$. We use $c=1$ and $c=0$ in the
unipolar and the bipolar cases, respectively. We take as 0 dB reference the
averaged value of $\|x\|^{2}$. We derive the optimal filter
$f_{\scriptscriptstyle\mathrm{opt}}(t)$ of linear reconstruction based on the
statistics of 1000 bandlimited input thus randomly drawn. Given our periodic
input setting, it is not possible to test the estimate $\hat{u}_{0}(t)$ in the
absolute leakage-free case $\alpha=0$ since the operator $A$ does not
converge. To simulate zero leakage, we set $\alpha$ to 0.03, which is close
enough to 0 while allowing numerically a reliable convergence of $A$. Our
observations are as follows:

  1. 1.

Accurate signal reconstruction is more difficult as leakage increases.

  2. 2.

$\widetilde{u}_{0}(t)$ is a better estimate than $\hat{u}_{0}(t)$ in all
cases. The improvement is particularly strong in the leakage-free unipolar
case.

  3. 3.

As expected, the estimate $u_{\scriptscriptstyle\mathrm{opt}}(t)$ is in all
cases as good as or better than $\widetilde{u}_{0}(t)$.

  4. 4.

The experiment tends to show that $\widetilde{u}_{0}(t)$ is optimal as a
convolutional reconstruction with the leakage-free unipolar configuration.
Indeed, not only $\widetilde{u}_{0}(t)$ is close to
$u_{\scriptscriptstyle\mathrm{opt}}(t)$ in MSE in the unipolar case with
$\alpha=0.03$, but $F_{\scriptscriptstyle\mathrm{opt}}(\omega)$ also appears
to be numerically close to the constant $\theta$ in $[-\pi,\pi]$. We attribute
its remaining difference with the ideal constant to boundary effects together
with a leakage that is not exactly 0.

  5. 5.

The estimates $\hat{u}_{0}(t)$ and $\widetilde{u}_{0}(t)$ are basically
useless in the unipolar case with leakage $\alpha\geq 1.5$ as their MSE’s are
of the same order of (if not larger than) the input’s squared sum.

##  IV Reconstruction by POCS

The signal estimation methods presented until now remain of limited accuracy
as can be seen in the results of Fig. 1. This can be explained by a limited
use of the analytical information contained in (5). An early observation is
that little is done in these methods to make the bandlimited estimates
consistent with this information. This was already mentioned concerning
$\widetilde{u}_{0}(t)$. Meanwhile, the estimate
$u_{\scriptscriptstyle\mathrm{opt}}(t)$ only focused on a convolutional
improvement of $\widetilde{u}_{0}(t)$ while the LIF encoder is neither linear
nor time invariant. We show in this section that there is a systematic way to
improve any estimate that is not consistent. The POCS method is an algorithm
that combines this improvement mechanism with the bandlimitation requirement,
to converge to an estimate that is both bandlimited and consistent with
iterative MSE reductions.

###  IV-A LIF encoder as generalized sampler

The starting point is to rigorously reformulate the signal problem in the
Hilbert space $L^{2}(D)$. This leads us to look at $x(t)$ as the input to the
LIF encoder while thinking of $c$ as a parameter of the encoder. In this
context, the MSE of an estimate $u(t)$ will refer to $\|u-x\|^{2}$. Next, (5)
is restated as

$\forall n\in{\sf N},\quad A_{t_{n-1}}x(t_{n})=\theta_{n}:=\epsilon_{n}\theta-
A_{t_{n-1}}c(t_{n})$ (13)

where $c(t):=c$ for all $t$. Let us introduce the notation

$I_{n}:=[t_{n-1},t_{n})\quad\mbox{and}\quad\Delta t_{n}:=t_{n}-t_{n-1},\qquad
n\in{\sf N}.$

From (3), (13) can be rewritten as

$\forall n\in{\sf N},\qquad\langle h_{n},x\rangle=\theta_{n}$ (14)

where $\langle\cdot,\cdot\rangle$ was defined in (6),

$h_{n}(t):=e^{-\alpha(t_{n}-t)}\,1_{I_{n}}(t)$ (15)

and $1_{I_{n}}(t)$ is the indicator function of the time interval $I_{n}$. In
(14), $A_{t_{n-1}}x(t_{n})$ is better seen as a linear functional mapping
$x(t)$ into a single scalar value $\theta_{n}$ that fulfills the role of
sample of $x(t)$, according to the generalized framework of sampling [31]. The
function $h_{n}(t)$ is the associated sampling kernel function. The sample
value $\theta_{n}$ is explicitly given by

$\theta_{n}=\epsilon_{n}\theta-c\int_{t_{n-1}}^{t_{n}}e^{-\alpha(t_{n}-s)}{\mathrm{d}}s=\epsilon_{n}\theta-\mbox{\footnotesize$\displaystyle\frac{c}{\alpha}$}(1-e^{-\alpha\Delta
t_{n}}).$

Now, another contribution of this presentation is the orthogonality that
results from $\langle\cdot,\cdot\rangle$. Induced by this inner-product, the
norm $\|\cdot\|$ satisfies the Pythagorean theorem which, in a stronger form,
lies in the equivalence

$\langle
u,v\rangle=0\qquad\Leftrightarrow\qquad\|u+v\|^{2}=\|u\|^{2}+\|v\|^{2}$ (16)

for any $u(t),v(t)\in L^{2}(D)$. This will be the key to MSE reductions as
seen in the next section.

###  IV-B Set theoretic estimation

Another way to state the property of (14) is to say that $x(t)$ belongs to the
set

${\cal C}({\bm{\theta}}):=\Big{\\{}u(t)\in L^{2}(D):\forall n\in{\sf
N},~{}\langle h_{n},u\rangle=\theta_{n}\Big{\\}}$ (17)

where ${\bm{\theta}}$ symbolizes the sequence $(\theta_{n})_{n\in{\sf N}}$. We
call the elements of ${\cal C}({\bm{\theta}})$ the estimates of $L^{2}(D)$
that are consistent with the sampling of $x(t)$ in the sense of (14). We will
use the shorter notation ${\cal C}:={\cal C}({\bm{\theta}})$ when there is no
ambiguity about the considered sequence ${\bm{\theta}}$, but will keep the
full notation ${\cal C}({\bm{\theta}})$ whenever the dependence with
${\bm{\theta}}$ needs to be highlighted. We alluded in Section III-A to the
difficulty to find a bandlimited signal that is simultaneously consistent.
Whenever an estimate $u(t)$ is in $L^{2}(D)$ but not in ${\cal C}$, there is
in fact a systematic procedure to reduce its MSE. This is due to the
outstanding property that ${\cal C}$ is an affine subspace of $L^{2}(D)$
(i.e., a translated linear subspace) which is moreover closed (since it is of
finite codimension). As a generalization from closed linear subspaces, every
signal $u(t)\in L^{2}(D)$ has an orthogonal projection $P_{\cal C}u(t)$ onto
${\cal C}$ in the sense of $\langle\cdot,\cdot\rangle$, which is the unique
element of ${\cal C}$ such that

$\forall v(t)\in{\cal C},\qquad\big{\langle}u-P_{\cal C}u,\,P_{\cal
C}u-v\big{\rangle}=0.$ (18)

Due to (16), this is equivalent to

$\forall v(t)\in{\cal C},\qquad\|u-v\|^{2}=\|u-P_{\cal C}u\|^{2}+\|P_{\cal
C}u-v\|^{2}.$ (19)

This gives the alternative characterization that $P_{\cal C}u(t)$ is the
element $v(t)$ of ${\cal C}$ that minimizes $\|u-v\|$, i.e., that is closest
to $u(t)$ in the MSE sense. As $x(t)\in{\cal C}$, (19) also implies with
$v(t)=x(t)$ that

$\|P_{\cal C}u-x\|\leq\|u-x\|$

with a strict inequality when $u(t)\notin{\cal C}$ since $\|u-P_{\cal C}u\|>0$
in (19). In other words, $P_{\cal C}u(t)$ is a better estimate of $x(t)$ than
$u(t)$.

###  IV-C Projection implementation

It remains to find the explicit expression of $P_{\cal C}$. What makes its
derivation easy is the outstanding property that $(h_{n}(t))_{n\in{\sf N}}$ is
an orthogonal family of functions in $L^{2}(D)$ since their time supports do
not overlap.

###### Proposition IV.1

$P_{{\cal C}({\bm{\theta}})}u(t)=u(t)+{\textstyle\sum\limits_{n\in{\sf
N}}}\frac{\theta_{n}-\langle h_{n},u\rangle}{\|h_{n}\|^{2}}\,h_{n}(t).$ (20)

###### Proof:

Let $p(t)$ be the right hand side of (20). Since $\langle
h_{m},h_{n}\rangle=0$ for any distinct $m,n\in{\sf N}$, then $\langle
h_{m},p\rangle=\langle h_{m},u\rangle+(\theta_{m}-\langle
h_{m},u\rangle)=\theta_{m}$ for any $m\in{\sf N}$. So $p(t)\in{\cal C}$. For
any $v(t)\in{\cal C}$, $\langle h_{n},p-v\rangle=\theta_{n}-\theta_{n}=0$ for
all $n\in{\sf N}$. Since $u(t)-p(t)=\sum_{n\in{\sf N}}\alpha_{n}h_{n}(t)$ for
some coefficients $\alpha_{n}$, then $\langle u-p,p-v\rangle=0$. Thus
$p(t)=P_{\cal C}u(t)$. ∎

###  IV-D POCS algorithm

We can see from (20), that the correcting term of $P_{\cal C}u(t)$ is a
piecewise exponential function, which is non-bandlimited. By considering the
filtered version

$\widetilde{P}_{\cal C}u(t):=\varphi(t)*P_{\cal C}u(t),$ (21)

we further reduce the error of $P_{\cal C}u(t)$ with $x(t)$. We thus obtain
the error reductions

$\|\widetilde{P}_{\cal C}u-x\|\leq\|P_{\cal C}u-x\|\leq\|u-x\|$ (22)

where at least one of the two inequalities is strict when $u(t)\notin{\cal
B}\cap{\cal C}$, i.e., when $u(t)$ is not simultaneously bandlimited and
consistent (note that $u(t)\in{\cal C}\backslash{\cal B}$ implies that
$P_{\cal C}u(t)=u(t)\notin{\cal B}$). For repetitive improvements, this
suggests the use of the algorithm

$u^{\scriptscriptstyle(k+1)}(t):=\widetilde{P}_{\cal
C}\,u^{\scriptscriptstyle(k)}(t),\qquad k\geq 0$ (23)

with $u^{\scriptscriptstyle(0)}(t):=u(t)$. As a common mathematical notation,
one writes

$u^{\scriptscriptstyle(k)}(t):={\widetilde{P}_{\cal C}}^{k}u(t),\qquad k\geq
0.$

From the set theoretic viewpoint, $\varphi(t)*v(t)$ is nothing but the
orthogonal projection of $v(t)$ onto ${\cal B}$, which is a closed linear
subspace of $L^{2}(D)$. Thus,

$\widetilde{P}_{\cal C}=P_{\cal B}P_{\cal C}.$

Thus, the estimates $u^{\scriptscriptstyle(k)}(t)$ are obtained by alternating
orthogonal projections between ${\cal C}$ and ${\cal B}$. This is a particular
case of the method of projection onto convex sets (POCS) [32, 33]. We know
that $\|u^{\scriptscriptstyle(k)}-x\|$ is strictly decreasing with $k$ as long
as $u^{\scriptscriptstyle(k)}(t)$ is not in ${\cal B}\cap{\cal C}$. In fact,
it is known that $u^{\scriptscriptstyle(k)}(t)$ must eventually converge to an
element of ${\cal B}\cap{\cal C}$. Moreover, as ${\cal B}$ and ${\cal C}$ are
convex sets that are more specifically affine spaces,

$u^{\scriptscriptstyle(\infty)}(t):=\lim_{k\rightarrow\infty}\widetilde{P}_{\cal
C}^{k}u(t)=P_{{\cal B}\cap{\cal C}}\,u(t)$ (24)

where the limit is in the sense of $L^{2}(D)$. Thus,
$u^{\scriptscriptstyle(\infty)}(t)$ is the element of ${\cal B}\cap{\cal C}$
that is closest to the initial guess $u(t)$ in the MSE sense. Finally, it
follows from (20) and (21) that

$\forall u(t)\in{\cal B},~{}~{}\widetilde{P}_{{\cal
C}({\bm{\theta}})}u(t)=u(t)+{\textstyle\sum\limits_{n\in{\sf
N}}}\frac{\theta_{n}-\langle
h_{n},u\rangle}{\|h_{n}\|^{2}}\,\widetilde{h}_{n}(t).$ (25)

###  IV-E Original POCS configuration in nonuniform sampling

The POCS method was originally introduced in nonuniform sampling in [41] with
a different implementation. Its extended version to LIF would be as follows.
We can write

${\cal B}\cap{\cal C}=\textstyle\bigcap\limits_{n\in{\sf N}}\widetilde{\cal
C}_{n}\quad\mbox{where}\quad\widetilde{\cal C}_{n}:=\big{\\{}v(t)\in{\cal
B}:\langle\widetilde{h}_{n},v\rangle=\theta_{n}\big{\\}}$

after noticing that $\langle h_{n},v\rangle=\langle\widetilde{h}_{n},v\rangle$
for all $v\in{\cal B}$, since $\langle h_{n}-\widetilde{h}_{n},v\rangle=0$ by
orthogonality in the frequency domain. So, another way to converge to an
element of ${\cal B}\cap{\cal C}$ is to iterate

$u^{\scriptscriptstyle(k+1)}(t):=P_{\widetilde{\cal C}_{N}}P_{\widetilde{\cal
C}_{N-1}}\cdots P_{\widetilde{\cal C}_{2}}P_{\widetilde{\cal
C}_{1}}u^{\scriptscriptstyle(k)}(t),\qquad k\geq 0.$ (26)

We have observed experimentally that this iteration has a similar convergence
rate to (23). However, its implementation is substantially more complicated.
For each $n\in{\sf N}$, it can be derived that

$\forall u(t)\in{\cal B},\qquad P_{\widetilde{\cal
C}_{n}}u(t)=u(t)+\frac{\theta_{n}-\langle\widetilde{h}_{n},u\rangle}{\|\widetilde{h}_{n}\|^{2}}\,\widetilde{h}_{n}(t).$

While the correcting action of $\widetilde{P}_{\cal C}$ in (25) is performed
by inclusion of all the estimation errors
$(\theta_{n}{-}\langle\widetilde{h}_{n},u\rangle)_{n\in{\sf N}}$ in one shot,
each elementary projection $P_{\widetilde{C}_{n}}$ in (26) needs to wait for
the previous projection $P_{\widetilde{C}_{n-1}}$ to be completed before it
can proceed. In that regard, (23) has some connection to the parallel version
of the POCS method [42].

###  IV-F Lazar’s algorithm

For reference, we will need to compare the POCS method with the first
iterative algorithm of bandlimited reconstruction for LIF proposed by Lazar in
[23]. Under our setting and notation, it consists in the iteration of

$\displaystyle
u^{\scriptscriptstyle(k+1)}(t)=Lu^{\scriptscriptstyle(k)}(t)\quad$ where
$\displaystyle Lu(t):=u(t)+{\textstyle\sum\limits_{n\in{\sf
N}}}\big{(}\theta_{n}-\langle h_{n},u\rangle\big{)}\,\varphi(t-\tau_{n})$

and $\tau_{n}:=(t_{n-1}+t_{n})/2$. Its application in the leakage-free case
was previously analyzed in [15] in the context of ASDM encoding. When ${\sf
N}={\mathbb{Z}}$ and the supremum of $(\Delta t_{n})_{n\in{\mathbb{Z}}}$ is
less than 1, it was shown that $u^{\scriptscriptstyle(k)}(t)$ tends to $x(t)$,
based on mathematical results of [43]. This convergence is shown in [23] to
hold in the presence of leakage under some stricter intricate condition on
$(\Delta t_{n})_{n\in{\mathbb{Z}}}$. An immediate advantage of the POCS method
is that it converges unconditionally when ${\sf N}$ is finite, which is always
the case in practice. This is whether bandlimited consistent reconstruction is
unique or not. When ${\sf N}$ is infinite, it is guaranteed to converge to
$x(t)$ at least whenever the conditions of Lazar’s algorithm are satisfied, as
explained in [36] in the absence of leakage. With a similar reasoning, one
arrives at the same conclusion in the leaky case. What remains is to compare
the numerical performances of these two methods, which is done in the next
section.

###  IV-G Numerical experiments

Figure 3: MSE performance of various iterative algorithms.

Under the same experimental conditions as in Fig. 2, we compare in Fig. 3 the
MSE of estimates from three iterative methods:

  * •

$L^{k}z(t)$, Lazar’s algorithm iterates starting from the zero function

$z(t):=0,$

  * •

${\widetilde{P}_{\cal C}}^{k}z(t)$, the POCS iterates starting from $z(t)$,

  * •

${\widetilde{P}_{\cal C}}^{k}u_{\scriptscriptstyle\mathrm{opt}}^{-c}(t)$, the
POCS iterates starting from
$u_{\scriptscriptstyle\mathrm{opt}}^{-c}(t):=u_{\scriptscriptstyle\mathrm{opt}}(t)-c$,
remembering that $u_{\scriptscriptstyle\mathrm{opt}}(t)$ is an estimate of
$x^{c}(t)=x(t)+c$.

As a fourth experiment, we have added the MSE of ${\widetilde{P}_{{\cal
C}^{\prime}}}^{k}z(t)$ which differs from ${\widetilde{P}_{\cal C}}^{k}z(t)$
only by the density of firing instants, with an oversampling ratio of $2$
instead of $1.5$. Our observations are as follows.

  1. 1.

$L^{k}z(t)$ is outperformed by ${\widetilde{P}_{\cal C}}^{k}z(t)$, except in
the leakage-free unipolar case where it yields a similar MSE. Moreover,
$L^{k}z(t)$ diverges with $k$ in the leaky cases. These are evidently cases
where the sufficient condition for convergence derived in [23] is not met.

  2. 2.

By applying ${\widetilde{P}_{\cal C}}^{k}$ to
$u_{\scriptscriptstyle\mathrm{opt}}^{-c}(t)$, one sees an example where the
POCS method can be used to improve an estimate that is provided by another
method. At the same time, one obtains a head start in MSE reduction compared
to ${\widetilde{P}_{\cal C}}^{k}z(t)$, which is particularly substantial in
the leaky unipolar cases. This however requires the availability of input-
output encoding statistics.

  3. 3.

By comparing ${\widetilde{P}_{{\cal C}^{\prime}}}^{k}z(t)$ with
${\widetilde{P}_{\cal C}}^{k}z(t)$, we see by how much the POCS convergence is
accelerated by increasing the sampling density by a factor of 1.33. The
acceleration is particularly dramatic with the leakage-free unipolar
configuration, where the number of iterations necessary to reach approximately
an MSE of -35dB is reduced by a factor of 4.

##  V Sampling pseudo-inversion

Until now, we have mostly focused on the algorithmic aspects of the POCS
method with iterative MSE decrease and guaranteed convergence. But a remaining
question is to characterize the limit of the POCS iteration from a signal
theoretic viewpoint. We first show that the sampling information of (5)
amounts to providing the image of the input $x(t)$ through a known linear
operator that we call $H$. We then show that the POCS solution amounts to
applying the pseudo-inverse of $H$ on this image, under some suitable
definition of the space norms. This view is actually not specific to LIF
encoding and is characteristic of a whole generic family of nonuniform
sampling problems, including the original basic problem of [37]. The
particularly powerful contribution of the pseudo-inversion property is the
effectiveness of the POCS method under sampling noise.

###  V-A Sampling operator

We introduced the POCS method by formalizing the sample values
$A_{t_{n-1}}x(t_{n})$ as inner-products $\langle h_{n},x\rangle$ in (14).
Thinking of ${\bm{\theta}}=(\theta_{n})_{n\in{\sf N}}$ as an element of
${\mathbb{R}}^{N}$, (14) can be more formally put in the form

$Hx(t)={\bm{\theta}}$ (27)

where $H$ is the linear operator

$H:\begin{array}[t]{ccl}{\cal B}&\rightarrow&{\mathbb{R}}^{N}\\\
u(t)&\mapsto&\big{(}\langle h_{n},u\rangle\big{)}_{n\in{\sf N}}\end{array}.$
(28)

As a natural reflex, one wishes to retrieve $x(t)$ from (27) by “inverting”
$H$. However, $H$ is not invertible in the most general case. In this
situation, one typically resort to the pseudo-inverse $H^{\dagger}$ of $H$ to
obtain the estimate

$x^{\dagger}(t):=H^{\dagger}{\bm{\theta}}.$

Such a principle of signal reconstruction from LIF encoding was first
mentioned in [24, §4.3]. In this work, the signals of ${\cal B}$ are
specifically represented by their discrete components in the sinc basis. This
allows to view $H$ as a matrix, for which the notion of pseudo-inverse is
commonly understood. This presentation however has conceptual limitations when
${\cal B}$ is of infinite dimension, which makes $H$ an infinitely fat matrix.
Whenever ${\bm{\theta}}$ is in the range of $H$, denoted by $\mathrm{ran}(H)$,
$H^{\dagger}{\bm{\theta}}$ has in fact a direct definition with the operator
presentation of (28) [44], which is

$H^{\dagger}{\bm{\theta}}\mbox{ is the minimal-norm element of
}H^{-1}({\bm{\theta}})$ (29)

where

$H^{-1}({\bm{\theta}}):=\big{\\{}u(t)\in{\cal
B}:Hu(t)={\bm{\theta}}\big{\\}}.$

Now, it follows from (17) that

$H^{-1}\big{(}{\bm{\theta}}\big{)}={\cal B}\cap{{\cal C}({\bm{\theta}})}.$
(30)

Then, $H^{\dagger}{\bm{\theta}}$ is the element of ${\cal B}\cap{\cal
C}({\bm{\theta}})$ that is closest to the zero signal $z(t)$. Therefore,
$H^{\dagger}{\bm{\theta}}$ is also the orthogonal projection of $z(t)$ onto
this set. It results from (24) that

$H^{\dagger}{\bm{\theta}}=P_{{\cal B}\cap{{\cal
C}({\bm{\theta}})}}\,z(t)=\lim_{k\rightarrow\infty}{\widetilde{P}_{{\cal
C}({\bm{\theta}})}}^{k}z(t).$ (31)

Hence, the POCS iteration started from 0 tends to the pseudo-inverse of $H$
applied to ${\bm{\theta}}$. This is mathematically true regardless of the
dimension of ${\cal B}$, whether it is finite, virtually infinite, or
theoretically infinite.

###  V-B Generalized problem of nonuniform sampling

The problem solved in the previous section is nothing but a generalization of
the basic nonuniform sampling problem of the pioneering paper by Yen [37].
Given the samples $x_{n}:=x(t_{n})$ of a bandlimited signal $x(t)$ at $N$
given instants $(t_{n})_{n\in{\sf N}}$, this article raised the question of
finding the minimal-norm bandlimited signal that interpolates these samples.
With a formalism similar to the previous section, the solution is
theoretically $S^{\dagger}(x_{n})_{n\in{\sf N}}$ where

$S:\begin{array}[t]{ccl}{\cal B}&\rightarrow&{\mathbb{R}}^{N}\\\
u(t)&\mapsto&\big{(}u(t_{n})\big{)}_{n\in{\sf N}}\end{array}.$ (32)

In their present forms, the operators $S$ and $H$ may not appear to be
similar. However, they are actually both of the form

$F:\begin{array}[t]{ccl}{\cal B}&\rightarrow&{\mathbb{R}}^{N}\\\
u(t)&\mapsto&\big{(}\langle f_{n},u\rangle\big{)}_{n\in{\sf N}}\end{array}$
(33)

for some family $\big{(}f_{n}(t)\big{)}_{n\in{\sf N}}$ of functions in ${\cal
B}$. Specifically,

  * •

$S=F$ with $f_{n}(t)=\varphi_{n}(t):=\varphi(t-t_{n})$ for all $n\in{\sf N}$,

  * •

$H=F$ with $f_{n}(t)=\widetilde{h}_{n}(t)$ for all $n\in{\sf N}$.

Indeed, for any $u(t)\in{\cal B}$,

$u(t_{n})=\langle\varphi_{n},u\rangle\qquad\mbox{and}\qquad\langle
h_{n},u\rangle=\langle\widetilde{h}_{n},u\rangle.$ (34)

This is justified as follows. With the even property of $\varphi(t)$, we have
$\langle\varphi_{n},u\rangle=$ $(\varphi*u)(t_{n})=u(t_{n})$. Meanwhile,
$\langle h_{n}-\widetilde{h}_{n},u\rangle$ must be 0 since
$h_{n}(t)-\widetilde{h}_{n}(t)$ has no in-band component. Both
$(\varphi(t{-}t_{n}))_{n\in{\sf N}}$ and $(\widetilde{h}_{n}(t))_{n\in{\sf
N}}$ are families of ${\cal B}$.

The outstanding result of [37] was that $S^{\dagger}(x_{n})_{n\in{\sf N}}$
must be a linear combination of $(\varphi(t-t_{n}))_{n\in{\sf N}}$. While this
was found by means of Lagrange multipliers, this is in fact a general but
elementary result of linear algebra, which we state here for any operator $F$
of the type (33).

###### Proposition V.1

For any ${\mathbf{v}}\in\mathrm{ran}(F)$, $F^{\dagger}{\mathbf{v}}$ is the
unique linear combination $u(t)$ of $(f_{n}(t))_{n\in{\sf N}}$ such that
$Fu(t)={\mathbf{v}}$.

###### Proof:

Let ${\cal V}_{f}$ be the linear span of $(f_{n}(t))_{n\in{\sf N}}$. Note
first that $Sv(t)=0$ if and only if $v(t)\in{\cal V}_{f}^{\perp}$. Next, the
function $x^{\dagger}(t):=F^{\dagger}{\mathbf{v}}$ is the minimal-norm element
of ${\cal B}$ such that $Sx^{\dagger}(t)={\mathbf{v}}$. Let $u(t)+v(t)$ be the
unique decomposition of $x^{\dagger}(t)$ in ${\cal V}_{f}{\oplus}{\cal
V}_{f}^{\perp}$. Since $Sv(t)=0$, then $Su(t)=Sx^{\dagger}(t)={\mathbf{v}}$.
We cannot have $v(t)\neq 0$ otherwise
$\|u\|^{2}<\|u\|^{2}+\|v\|^{2}=\|x^{\dagger}\|^{2}$ which is not possible by
construction of $x^{\dagger}(t)$. So $x^{\dagger}(t)=u(t)\in{\cal V}_{f}$. Let
$w(t)\in{\cal V}_{f}$. If $Sw(t)={\mathbf{v}}$, then
$S(w(t){-}x^{\dagger}(t))=0$, and hence $w(t){-}x^{\dagger}(t)\in{\cal
V}_{f}^{\perp}$. This is possible only when $w(t)=x^{\dagger}(t)$. ∎

We thus know that $H^{\dagger}{\bm{\theta}}$ must be a linear combination of
$(\widetilde{h}_{n}(t))_{n\in{\sf N}}$. Note that this result was already
observable from (31) given the expression of (25) and the identity $z(t)=0$.

###  V-C POCS iteration under sampling noise

The relation (31) of the POCS iteration to the pseudo-inverse of $H$ was valid
under the assumption that the identity $\langle h_{n},x\rangle=\theta_{n}$ of
(6) is exact for every $n\in{\sf N}$. In practice, however, the sequence
$(t_{n})_{n\in{\sf N}}$ output by the LIF encoder includes errors with respect
to the true crossing instants of $|A_{t_{n-1}}x^{c}(t)|$, due to noise or
other factors. We will see in Section V-F the special case of errors due to
time quantization. Since the sequences $(h_{n}(t))_{n\in{\sf N}}$ and
$(\theta_{n})_{n\in{\sf N}}$ used in the reconstruction are explicitly
calculated from the deviated sequence $(t_{n})_{n\in{\sf N}}$ via (15) and
(13), we can only write that

$\langle h_{n},x\rangle=\theta_{n}^{0}:=\theta_{n}-e_{n},\qquad n\in{\sf N}$
(35)

where $e_{n}$ is some unknown deviation. In other words, the acquired inner-
product value $\theta_{n}$ is thought of as the true inner-product value
$\theta_{n}^{0}$ plus an error $e_{n}$. The question is whether the POCS
iteration still remains convergent, and if so, what is the meaning of its
limit? To answer this question, let us introduce in ${\mathbb{R}}^{N}$ the
weighted inner-product

$\langle{\mathbf{u}},{\mathbf{v}}\rangle_{N}:={\textstyle\sum\limits_{n\in{\sf
N}}}\frac{u_{n}v_{n}}{\|h_{n}\|^{2}},\qquad{\mathbf{u}},{\mathbf{v}}\in{\mathbb{R}}^{N}$

where $u_{n}$ and $v_{n}$ are the $n$th components of ${\mathbf{u}}$ and
${\mathbf{v}}$, respectively. We first have the following result.

###### Proposition V.2

For any ${\bm{\theta}}=(\theta_{n})_{n\in{\sf N}}\in{\mathbb{R}}^{N}$,

$\widetilde{P}_{{\cal C}({\bm{\theta}})}=\widetilde{P}_{{\cal
C}({\bm{\bar{\theta}}})}$ (36)

where, for any vector ${\mathbf{v}}\in{\mathbb{R}}^{N}$, $\mathbf{\bar{v}}$
designates the orthogonal projection of ${\mathbf{v}}$ onto $\mathrm{ran}(H)$
with respect to $\langle\cdot,\cdot\rangle_{N}$.

###### Proof:

Let $D_{\bm{\theta}}u(t):=P_{{\cal C}({\bm{\bar{\theta}}})}u(t)-P_{{\cal
C}({\bm{\theta}})}u(t)$ for any given $u(t)\in L^{2}(D)$. It is clear from
(20) that

$D_{\bm{\theta}}u(t)={\textstyle\sum\limits_{n\in{\sf
N}}}\mbox{\footnotesize$\displaystyle\frac{\bar{\theta}_{n}-\theta_{n}}{\|h_{n}\|^{2}}$}\,h_{n}(t).$

Then, for any $w(t)\in{\cal B}$,

$\big{\langle}D_{\bm{\theta}}u,w\big{\rangle}={\textstyle\sum\limits_{n\in{\sf
N}}}\mbox{\footnotesize$\displaystyle\frac{(\bar{\theta}_{n}-\theta_{n})\langle
h_{n},w\rangle}{\|h_{n}\|^{2}}$}=\big{\langle}{\bm{\bar{\theta}}}-{\bm{\theta}},Hw\big{\rangle}_{N}.$

This is zero since ${\bm{\bar{\theta}}}-{\bm{\theta}}$ is by construction
orthogonal to $\mathrm{ran}(H)$ in the sense of
$\langle\cdot,\cdot\rangle_{N}$. Since this is true for all $w(t)\in{\cal B}$,
then $D_{\bm{\theta}}u(t)\in{\cal B}^{\perp}$. Thus, $P_{\cal
B}D_{\bm{\theta}}u(t)=0$, which implies (36). ∎

Since ${\bm{\bar{\theta}}}\in\mathrm{ran}(H)$, the identity (31) applies to
${\bm{\bar{\theta}}}$. It follows that

$\lim_{k\rightarrow\infty}\widetilde{P}_{{\cal
C}({\bm{\theta}})}^{k}z(t)=\lim_{k\rightarrow\infty}\widetilde{P}_{{\cal
C}({\bm{\bar{\theta}}})}^{k}z(t)=H^{\dagger}{\bm{\bar{\theta}}}.$ (37)

Since ${\bm{\theta}}={\bm{\theta}}^{0}-{\mathbf{e}}$ and
${\bm{\theta}}^{0}\in\mathrm{ran}(H)$ according to (35), then

${\bm{\bar{\theta}}}={\bm{\theta}}^{0}-\mathbf{\bar{e}}.$

Now, the remaining error $\mathbf{\bar{e}}$ is such that

$\|\mathbf{\bar{e}}\|_{N}^{2}\leq\|\mathbf{\bar{e}}\|_{N}^{2}+\|{\mathbf{e}}-\mathbf{\bar{e}}\|_{N}^{2}=\|{\mathbf{e}}\|_{N}^{2}$
(38)

by the Pythagorean theorem. So, not only the POCS iteration keeps on
converging, but it has also a filtering action on the sampling errors, which
is all the greater that $\mathrm{ran}(H)$ is small compared to
${\mathbb{R}}^{N}$. Meanwhile, any error that lies in $\mathrm{ran}(H)$ such
as $\mathbf{\bar{e}}$ is irreversible because undistinguishable from the
samples of an actual input deviation in ${\cal B}$.

###  V-D Complete pseudo-inversion of sampling operator

The characterization (29) defined $H^{\dagger}$ only in $\mathrm{ran}(H)$. Let
us extend the definition of $H^{\dagger}$ to ${\mathbb{R}}^{N}$ by taking
$H^{\dagger}{\bm{\theta}}:=H^{\dagger}{\bm{\bar{\theta}}}$ for all
${\bm{\theta}}\in{\mathbb{R}}^{N}$ (which is trivially true when
${\bm{\theta}}\in\mathrm{ran}(H)$). In this case, we can write from (37) that

$\forall{\bm{\theta}}\in{\mathbb{R}}^{N},\qquad\lim_{k\rightarrow\infty}{\widetilde{P}_{{\cal
C}({\bm{\theta}})}}^{k}z(t)=H^{\dagger}{\bm{\theta}}.$ (39)

With (29) and this extension, $H^{\dagger}{\bm{\theta}}$ can be equivalently
characterized for all ${\bm{\theta}}\in{\mathbb{R}}^{N}$ as

$\displaystyle H^{\dagger}{\bm{\theta}}\mbox{ is the minimal-norm element of
}H^{-1}({\bm{\bar{\theta}}})$ where
$\displaystyle{\bm{\bar{\theta}}}:=\operatorname*{argmin}\limits_{{\mathbf{u}}\in\mathrm{ran}(H)}\|{\mathbf{u}}-{\bm{\theta}}\|_{N},$

since ${\bm{\bar{\theta}}}$ is also by orthogonal projection the element of
$\mathrm{ran}(H)$ that is closest to ${\bm{\theta}}$ with respect to
$\|\cdot\|_{N}$. This description of $H^{\dagger}$ precisely makes it the
pseudo-inverse of $H$ as a mapping between the inner-product spaces
$\big{(}{\cal B},\langle\cdot,\cdot\rangle\big{)}$ and
$\big{(}{\mathbb{R}}^{N},\langle\cdot,\cdot\rangle_{N}\big{)}$ [44]. In this
setting, the POCS iteration completely solves the pseudo-inversion of $H$.

Working with the weighted inner-product $\langle\cdot,\cdot\rangle_{N}$
actually amounts to dealing with the canonical inner-product of
${\mathbb{R}}^{N}$ but using the normalized operator $\hat{H}$ such that

$\displaystyle\forall u(t)\in{\cal B},\quad$
$\displaystyle\hat{H}u(t):=\big{(}\langle\hat{h}_{n},u\rangle\big{)}_{n\in{\sf
N}}$ where $\displaystyle\hat{h}_{n}(t):=h_{n}(t)/\|h_{n}\|,$
$\displaystyle\quad n\in{\sf N}.$

Indeed, it is first clear that $Hx(t)={\bm{\theta}}$ if and only if
$\hat{H}x(t)={\bm{\hat{\theta}}}:=\big{(}\theta_{n}/\|h_{n}\|\big{)}_{n\in{\sf
N}}$. Then, as $\|{\bm{\theta}}\|_{N}=\|{\bm{\hat{\theta}}}\|$ where
$\|\cdot\|$ designates here the Euclidean norm of ${\mathbb{R}}^{N}$, one can
verify that $H^{\dagger}{\bm{\theta}}=\hat{H}^{\dagger}{\bm{\hat{\theta}}}$
where the pseudo-inverse $\hat{H}^{\dagger}$ is defined with respect to this
Euclidean norm. Meanwhile, note that $(\hat{h}_{n}(t))_{n\in{\sf N}}$ is an
orthonormal family of $L^{2}(D)$. For notation simplicity, however, we will
keep working with $H$ and its own definition of pseudo-inverse.

###  V-E Perfect reconstruction and noise filtering

When the sampling data of (14) is free of noise, the POCS iteration leads to
perfect reconstruction whenever ${\cal B}\cap C({\bm{\theta}})$ is a
singleton. As this set coincides with $H^{-1}({\bm{\theta}})$ according to
(30), this is realized if and only if $H$ is injective, i.e., is a one-to-one
mapping. This can theoretically happen only when $\dim({\cal B})\leq N$. This
is the case where $D=[0,T]$ with $T\leq N$, since $\dim({\cal B})=T$ given our
setting of Nyquist period 1 in Section II-B. In other terms, the oversampling
ratio $\rho:=N/T$ must be larger than or equal to 1. As
$(\tilde{h}_{n}(t))_{n\in{\sf N}}$ can be shown to be systematically
overcomplete in ${\cal B}$ in this case, this condition is also sufficient for
the injectivity of $H$.

In the presence of sampling noise, the error reduction of (38) is effective
when $\mathbf{\bar{e}}\neq{\mathbf{e}}$, which requires that $\mathrm{ran}(H)$
be a proper subspace of ${\mathbb{R}}^{N}$, and hence that
$(h_{n}(t))_{n\in{\sf N}}$ be linearly dependent. It can be shown that this
happens only when $\mathrm{rank}(H)=\dim(B)<N$. This is the more specific case
of strict oversampling $\rho>1$. The higher the oversampling ratio is, the
smaller $\mathrm{rank}(H)$ is compared to $N$, and the stronger is the noise
reduction.

###  V-F Time quantization

Regardless of the type of analog-to-digital conversion, the values acquired
from the input must eventually be quantized. In the present case of LIF
encoding, it is the time sequence $(t_{n})_{n\in{\sf N}}$ that must be
quantized. This is a fundamental source of sampling deviation as formulated in
(35). We show in Fig. 4 the MSE of the iterates $\widetilde{P}_{\cal
C}^{k}z(t)$ resulting from this artifact with the leakage-free unipolar
configuration of LIF in the same experimental conditions as in Fig. 2 with
however an oversampling ratio of 2. We perform the experiments with various
quantization resolutions of $(t_{n})_{n\in{\sf N}}$. We say that the
resolution is $b$-bit when the quantization step size is the Nyquist period
divided by $2^{b}$. We represent in mixed lines the MSE obtained from the
pseudo-inverse reconstruction $x^{\dagger}(t)=H^{\dagger}{\bm{\theta}}$. We
observe that the MSE of $\widetilde{P}_{\cal C}^{k}z(t)$ converges to that of
$x^{\dagger}(t)$ for each of the bit resolutions. We also observe the
satisfactory result that the MSE of pseudo-inversion is approximately reduced
by 6dB per bit of resolution.

Figure 4: MSE of $\widetilde{P}_{\cal C}^{k}z(t)$ versus the iteration number
$k$, with the leakage-free unipolar configuration and various resolutions of
time quantization. The MSE of $x^{\dagger}(t)=H^{\dagger}{\bm{\theta}}$ is
represented as a constant in mixed lines.

##  VI Discretization of the POCS iteration

The heavy part in the iteration of (23) is the continuous-time inner-product
$\langle h_{n},u\rangle$ that appears in (25). As the argument $u(t)$ is
bandlimited, this inner-product can be computed in discrete-time using the
Nyquist-rate samples of $u(t)$. But this makes (25) a complicated operation as
it requires a uniform sampling to obtain the coefficients of
$(\widetilde{h}_{n}(t))_{n\in{\sf N}}$ which are themselves not uniformly
distributed in time and not at the same rate as the Nyquist clock. Together
with intricate algebra, this requires in practice the implementation of
complicated buffering systems. There are however techniques to tackle these
issues, which we present in this section.

###  VI-A Case of zero initial estimate

We present first a convenient and rigorous way to discretize the POCS
iteration for $z^{\scriptscriptstyle(k)}(t)={\widetilde{P}_{\cal C}}^{k}z(t)$.
It follows from (23) and (25) that

$z^{\scriptscriptstyle(k+1)}(t)=z^{\scriptscriptstyle(k)}(t)+{\textstyle\sum\limits_{n\in{\sf
N}}}r_{n}^{\scriptscriptstyle(k)}\,\widetilde{h}_{n}(t)$ (40)

where

$r_{n}^{\scriptscriptstyle(k)}:=\frac{\theta_{n}-\langle
h_{n},z^{\scriptscriptstyle(k)}\rangle}{\|h_{n}\|^{2}},\qquad k\geq
0,~{}n\in{\sf N}.$ (41)

As $z^{\scriptscriptstyle(0)}(t)=z(t)=0$, then ${\widetilde{P}_{\cal
C}}^{k}z(t)$ must be of the form

${\widetilde{P}_{\cal C}}^{k}z(t)={\textstyle\sum\limits_{n\in{\sf
N}}}c_{n}^{\scriptscriptstyle(k)}\,\widetilde{h}_{n}(t)$ (42)

where the coefficients $c_{n}^{\scriptscriptstyle(k)}$ satisfy for each
$n\in{\sf N}$ the recursive relation

$c_{n}^{\scriptscriptstyle(k+1)}=c_{n}^{\scriptscriptstyle(k)}+r_{n}^{\scriptscriptstyle(k)}$

starting from $c_{n}^{\scriptscriptstyle(0)}:=0$. For each $m\in{\sf N}$, (40)
implies that

$\langle h_{m},z^{\scriptscriptstyle(k+1)}\rangle=\langle
h_{m},z^{\scriptscriptstyle(k)}\rangle+{\textstyle\sum\limits_{n\in{\sf
N}}}r_{n}^{\scriptscriptstyle(k)}\,\langle h_{m},\widetilde{h}_{n}\rangle.$

It then follows from (41) that

$r_{m}^{\scriptscriptstyle(k+1)}=r_{m}^{\scriptscriptstyle(k)}-{\textstyle\sum\limits_{n\in{\sf
N}}}r_{n}^{\scriptscriptstyle(k)}\frac{\langle
h_{m},\widetilde{h}_{n}\rangle}{\|h_{m}\|^{2}},\qquad m\in{\sf N}.$ (43)

The vectors
${\mathbf{r}}^{\scriptscriptstyle(k)}=(r_{n}^{\scriptscriptstyle(k)})_{n\in{\sf
N}}$ and
${\mathbf{c}}^{\scriptscriptstyle(k)}=(c_{n}^{\scriptscriptstyle(k)})_{n\in{\sf
N}}$ seen as $N$-dimensional column vectors then satisfy the linear system

$\displaystyle{\mathbf{r}}^{\scriptscriptstyle(k+1)}$
$\displaystyle={\mathbf{r}}^{\scriptscriptstyle(k)}-{\mathbf{H}}\,{\mathbf{r}}^{\scriptscriptstyle(k)}$
(44a) $\displaystyle{\mathbf{c}}^{\scriptscriptstyle(k+1)}$
$\displaystyle={\mathbf{c}}^{\scriptscriptstyle(k)}+{\mathbf{r}}^{\scriptscriptstyle(k)}$
(44b)

starting from

${\mathbf{r}}^{\scriptscriptstyle(0)}:=\big{(}\theta_{n}/\|h_{n}\|^{2}\big{)}_{n\in{\sf
N}}\qquad\mbox{and}\qquad{\mathbf{c}}^{\scriptscriptstyle(0)}:=0,$ (45)

where ${\mathbf{H}}$ is the $N{\times}N$ matrix

${\mathbf{H}}:=\Big{[}\mbox{\footnotesize$\displaystyle\frac{\langle
h_{m},\widetilde{h}_{n}\rangle}{\|h_{m}\|^{2}}$}\Big{]}_{(m,n)\in{\sf
N}^{2}}.$

The discrete-time system (44) is synchronous with the firing instants
$(t_{n})_{n\in{\sf N}}$. To obtain the estimate ${\widetilde{P}_{\cal
C}}^{p}z(t)$ at a targeted iteration number $p$, one simply needs to iterate
(44) $p$ times, and perform the discrete-time to continuous-time conversion of
(42) only once at $k=p$.

###  VI-B Precomputation of inner-products

In the system (44), one does not escape from having inner-products, since
$\langle h_{m},\widetilde{h}_{n}\rangle$ appears in the entries of
${\mathbf{H}}$. But they are iterate-independent. So, not only they can be
calculated before the iteration is started, but they can also be determined
analytically without Nyquist-rate resampling. Defining

$T_{n}^{m}:=\int_{t_{n}}^{t_{m}}{\mathrm{d}}t=t_{m}-t_{n},$ (46)

we prove in Appendix -A the following two results.

###### Proposition VI.1

For any given function $f(t)$, let $f_{n}(t):=f(t-t_{n})$. Then,

$\big{\langle}h_{m},f_{n}\big{\rangle}=g(T^{m}_{n})-e^{-\alpha\Delta
t_{m}}g(T^{m-1}_{n})$ (47)

where

$g(t):=\displaystyle\int_{0}^{t}e^{\alpha(s-t)}f(s){\mathrm{d}}s.$ (48)

###### Corollary VI.2

$\displaystyle\langle h_{m},\widetilde{h}_{n}\rangle$
$\displaystyle=e^{-\alpha\Delta t_{n}}\big{(}g(T^{m}_{n-1})-e^{-\alpha\Delta
t_{m}}g(T^{m-1}_{n-1})\big{)}$
$\displaystyle\quad-\big{(}g(T^{m}_{n})-e^{-\alpha\Delta
t_{m}}g(T^{m-1}_{n})\big{)}$ (49)

where

$\displaystyle g(t)$
$\displaystyle:=\frac{1}{\alpha}\int_{0}^{t}\mathrm{sinh}(\alpha(t-s))\,\varphi(s)\,{\mathrm{d}}s.$
(50)

The above function $g(t)$ only depends on the leakage coefficient $\alpha$. It
can be precomputed and numerically stored in a lookup table. The computation
of (49) in terms of $(t_{n})_{n\in{\sf N}}$ then becomes simple, and only
needs to be performed once before the iteration. For the complete
determination of ${\mathbf{H}}$, we need to derive $\|h_{m}\|^{2}$ that yields

$\|h_{m}\|^{2}=\int_{t_{m-1}}^{t_{m}}e^{-2\alpha(t_{m}-t)}{\mathrm{d}}t={\textstyle\frac{1}{2\alpha}}\big{(}1-e^{-2\alpha\Delta
t_{m}}\big{)}.$ (51)

###  VI-C Case of nonzero initial estimate

When one is more interested in the iterates ${\widetilde{P}_{\cal C}}^{k}u(t)$
for some nonzero initial estimate $u(t)$, the simple implementation of Section
VI-A can be used up to some space translation. This is specifically performed
as follows.

###### Proposition VI.3

For any $u(t)\in{\cal B}$ and $k\geq 0$,

${\widetilde{P}_{{\cal C}({\bm{\theta}})}}^{k}u(t)=u(t)+\widetilde{P}_{{\cal
C}({\bm{\hat{\theta}}})}^{k}z(t)$ (52)

where

${\bm{\hat{\theta}}}:={\bm{\theta}}-Hu(t).$ (53)

We prove this in Appendix -B. Note as a straightforward consequence of (39)
that

$\lim_{k\rightarrow\infty}{\widetilde{P}_{{\cal
C}({\bm{\theta}})}}^{k}u(t)=u(t)+H^{\dagger}{\bm{\hat{\theta}}}.$

However, a price to pay for starting the iteration from $u(t)$ is the required
derivation of $Hu(t)=(\langle h_{n},u\rangle)_{n\in{\sf N}}$ involved in (53).
Its computation complexity depends on $u(t)$. If for example one chooses

$u(t)={\textstyle\sum\limits_{n\in{\sf N}}}\epsilon_{n}\,f(t-t_{n})$

for some function $f(t)$, like in (12), then

$\langle h_{m},u\rangle=\textstyle\sum\limits_{n\in{\sf
N}}\epsilon_{n}\,\langle h_{m},f_{n}\rangle$

where $f_{n}(t):=f(t-t_{n})$. Then, the inner-products $\langle
h_{m},f_{n}\rangle$ are given in (47).

###  VI-D Future research

An important aspect of the implementation is yet to be investigated. The
matrix ${\mathbf{H}}$ involved in the discrete-time implementation of the POCS
algorithm in (44) has a size $N$ that is virtually infinite compared to the
practical windows of signal processing. A truncation of ${\mathbf{H}}$
automatically needs to be considered for the realistic implementation of (44a)
with expected performance degradations. Given the decay of the inner-products
$\langle h_{m},\widetilde{h}_{n}\rangle$ as $|m{-}n|$ increases, this is
similar to the problem of FIR filter windowing. The difficulty here is that
the filter is both time-varying and iterated. Some empirical truncation
experiments have been performed in [36] in the leakage-free case of signal
reconstruction from ASDM outputs. But, beyond a needed generalization to leaky
integration, this study deserves deeper theoretical investigations to be
considered for future research.

##  VII Conclusion

While LIF is commonly viewed as a neuron-inspired system that transforms
continuous-time signals into spike trains, we approached the problem of
bandlimited signal reconstruction from LIF outputs from the more general and
abstract perspective of nonuniform sampling. In this view, each output spike
gives the knowledge of an inner-product of the input with a kernel function
that can be determined from the characteristics of the encoder. With the
complete output of the LIF encoder, one thus has access to the transformation
the input by a known linear operator, which we call the sampling operator.
Signal reconstruction then consists in the pseudo-inversion of this operator.
This approach covers all situations including perfect reconstruction when it
is unique, minimal-norm reconstruction when the set of solutions is infinite,
and reconstruction under sampling errors such as deviations from time
quantization. On the practical side, this inversion can be systematically
implemented by the POCS method which, in the present case of linear pseudo-
inversion, consists in iterating orthogonal projections onto affine subspaces.
A single POCS iteration can actually be used to improve in MSE any estimate,
including signals provided by other reconstruction methods. The MSE reduction
is effective whenever the estimate is either not bandlimited or not consistent
with the considered output spike train. Although the POCS method is defined on
continuous-time signals, the POCS iteration can be rigorously implemented in
discrete-time in synchronization with the firing instants and without Nyquist-
rate resampling, up to the use of a lookup table.

The present version of the POCS algorithm only relied on the orthogonality of
the sampling kernel functions in $L^{2}(D)$, and thus can be extended to any
other nonuniform sampling scheme that carries this feature. For example, LIF
with a refractory period [45] fits in this context, with just a small
modification of the kernel functions.

A point of investigation left for future research is the practical
implementation of the POCS method by approximate time-varying FIR filters.

###  -A Proofs of Proposition VI.1 and Corollary VI.2

####  -A1 Preliminary

###### Lemma .1

For any function $f(t)$,

$\displaystyle(h_{m}*f)(t)=e^{-\alpha\Delta
t_{m}}\,\ell(t{-}t_{m-1})-\ell(t{-}t_{m})$ (54)

where

$\ell(t):=\displaystyle\int_{0}^{t}e^{\alpha(t-s)}f(s){\mathrm{d}}s.$

###### Proof:

It follows from (15) that

$\displaystyle(h_{m}*f)(t)$
$\displaystyle=\int_{t_{m-1}}^{t_{m}}e^{\alpha(s-t_{m})}f(t-s){\mathrm{d}}s$
$\displaystyle=\int^{t-t_{m-1}}_{t-t_{m}}e^{\alpha(t-s-
t_{m})}f(s){\mathrm{d}}s$ $\displaystyle=e^{-\alpha\Delta
t_{m}}\int_{0}^{t-t_{m-1}}e^{\alpha(t-t_{m-1}-s)}f(s){\mathrm{d}}s$
$\displaystyle\quad-\int_{0}^{t-t_{m}}e^{\alpha(t-t_{m}-s)}f(s){\mathrm{d}}s.$

∎

####  -A2 Proof of Proposition VI.1

We have

$\langle
h_{m},f_{n}\rangle=(h_{m}*\hat{f})(t_{n})\quad\mbox{where}\quad\hat{f}(t):=f(-t).$

By applying (54) and (46), we obtain

$\displaystyle\big{\langle}h_{m},f_{n}\big{\rangle}$
$\displaystyle=e^{-\alpha\Delta t_{m}}\ell(T^{n}_{m-1})-\ell(T^{n}_{m})$ (55)
$\displaystyle\mbox{where}\qquad\qquad\ell(t)$
$\displaystyle:=\displaystyle\int_{0}^{t}e^{\alpha(t-s)}\hat{f}(s){\mathrm{d}}s$
$\displaystyle=-\int_{0}^{-t}e^{\alpha(t+s)}f(s){\mathrm{d}}s=-g(-t)\qquad$

and $g(t)$ is given in (48). Then, (55) leads to (47).

####  -A3 Proof of Corollary VI.2

From (54), we have

$\displaystyle\widetilde{h}_{n}(t)=(h_{n}*\varphi)(t)=e^{-\alpha\Delta
t_{n}}\ell_{n-1}(t)-\ell_{n}(t)$
$\displaystyle\mbox{where}\quad\ell_{n}(t):=\ell(t-t_{n})\quad\mbox{and}\quad\ell(t):=\displaystyle\int_{0}^{t}e^{\alpha(t-s)}\varphi(s){\mathrm{d}}s,$

so that $\big{\langle}h_{m},\widetilde{h}_{n}\big{\rangle}=e^{-\alpha\Delta
t_{n}}\big{\langle}h_{m},\ell_{n-1}\big{\rangle}-\big{\langle}h_{m},\ell_{n}\big{\rangle}.$
It then follows from (47) that

$\displaystyle\big{\langle}h_{m},\widetilde{h}_{n}\big{\rangle}=e^{-\alpha\Delta
t_{n}}\big{(}g(T^{m}_{n-1})-e^{-\alpha\Delta t_{m}}g(T^{m-1}_{n-1})\big{)}$
$\displaystyle\hskip 71.13188pt-\big{(}g(T^{m}_{n})-e^{-\alpha\Delta
t_{m}}g(T^{m-1}_{n})\big{)}$ $\displaystyle\mbox{where}\qquad
g(t):=\displaystyle\int_{0}^{t}e^{\alpha(s-t)}\ell(s){\mathrm{d}}s$
$\displaystyle=\int_{0}^{t}e^{\alpha(2s-t)}\displaystyle\int_{0}^{s}e^{-\alpha\tau}\varphi(\tau){\mathrm{d}}\tau{\mathrm{d}}s$
$\displaystyle\textstyle=\left[\frac{e^{\alpha(2s-t)}}{2\alpha}\int\limits_{0}^{s}e^{-\alpha\tau}\varphi(\tau)\,{\mathrm{d}}\tau\right]_{s=0}^{t}\\!\\!\\!-\int\limits_{0}^{t}\frac{e^{\alpha(2s-t)}}{2\alpha}e^{-\alpha
s}\varphi(s)\,{\mathrm{d}}s$
$\displaystyle=\mbox{\footnotesize$\displaystyle\frac{1}{2\alpha}$}\int_{0}^{t}e^{\alpha(t-\tau)}\,\varphi(\tau)\,{\mathrm{d}}\tau-\mbox{\footnotesize$\displaystyle\frac{1}{2\alpha}$}\int_{0}^{t}e^{-\alpha(t-s)}\,\varphi(s)\,{\mathrm{d}}s$
$\displaystyle=\mbox{\footnotesize$\displaystyle\frac{1}{\alpha}$}\int_{0}^{t}\mathrm{sinh}(\alpha(t-s))\,\varphi(s)\,{\mathrm{d}}s.$

###  -B Proof of Proposition VI.3

For the given bandlimited function $u(t)$, consider the space translation

$Tv(t):=v(t)-u(t).$

Since $u(t)\in{\cal B}$,

$P_{\cal B}Tv(t)=P_{\cal B}v(t)-u(t)=TP_{\cal B}v(t).$ (56)

With $\widehat{\cal C}:=T({\cal C})$ and $\hat{v}(t):=Tv(t)$, it is easy to
verify that

$P_{\widehat{\cal C}}\hat{v}(t)-\hat{v}(t)=P_{\cal C}v(t)-v(t).$

Thus, $P_{\widehat{\cal C}}\hat{v}(t)=P_{\cal C}v(t)-v(t)+\hat{v}(t)=P_{\cal
C}v(t)-u(t)$, which implies that

$P_{\widehat{\cal C}}Tv(t)=TP_{\cal C}v(t).$ (57)

By combining (56) and (57), we get

$\widetilde{P}_{\widehat{\cal C}}T=P_{\cal B}P_{\widehat{\cal C}}T=P_{\cal
B}TP_{\cal C}=TP_{\cal B}P_{\cal C}=T\widetilde{P}_{\cal C}.$

By induction, one obtains more generally

$\widetilde{P}_{\widehat{\cal C}}^{k}T=T{\widetilde{P}_{\cal C}}^{k},\qquad
k\geq 0.$

Then

$\displaystyle\widetilde{P}_{\widehat{\cal C}}^{k}z(t)$
$\displaystyle=\widetilde{P}_{\widehat{\cal C}}^{k}Tu(t)=T{\widetilde{P}_{\cal
C}}^{k}u(t)={\widetilde{P}_{\cal C}}^{k}u(t)-u(t).$ (58)

Next,

$\displaystyle v(t)\in\widehat{\cal C}=T({\cal C})\quad$
$\displaystyle\Leftrightarrow\quad v(t)+u(t)\in{\cal C}$
$\displaystyle\Leftrightarrow\quad\forall n\in{\sf N},\quad\langle
h_{n},v+u\rangle=\theta_{n}$ $\displaystyle\Leftrightarrow\quad\forall
n\in{\sf N},\quad\langle h_{n},v\rangle=\hat{\theta}_{n}$

where $\hat{\theta}_{n}$ is the $n$th component of ${\bm{\hat{\theta}}}$.
Thus, $\widehat{\cal C}={\cal C}({\bm{\hat{\theta}}})$. So (58) implies (52).

## References

  * [1] Y. C. Yoon, “LIF and simplified SRM neurons encode signals into spikes via a form of asynchronous pulse sigma–delta modulation,” IEEE transactions on neural networks and learning systems, vol. 28, no. 5, pp. 1192–1205, 2016.
  * [2] I. Kianpour, B. Hussain, H. S. Mendonca, and V. G. Tavares, “System-level study on impulse-radio integration-and-fire (irif) transceiver,” AEU-International Journal of Electronics and Communications, vol. 112, p. 152896, 2019.
  * [3] P. Mart\́mathrm{i}nez-Nuevo, H.-Y. Lai, and A. V. Oppenheim, “Delta-ramp encoder for amplitude sampling and its interpretation as time encoding,” IEEE Transactions on Signal Processing, vol. 67, no. 10, pp. 2516–2527, 2019.
  * [4] K. Adam, A. Scholefield, and M. Vetterli, “Sampling and reconstruction of bandlimited signals with multi-channel time encoding,” IEEE Transactions on Signal Processing, vol. 68, pp. 1105–1119, 2020.
  * [5] B. A. Moser and M. Lunglmayr, “On quasi-isometry of threshold-based sampling,” IEEE Transactions on Signal Processing, vol. 67, no. 14, pp. 3832–3841, 2019.
  * [6] R. Alexandru and P. L. Dragotti, “Reconstructing classes of non-bandlimited signals from time encoded information,” IEEE Transactions on Signal Processing, vol. 68, pp. 747–763, 2020.
  * [7] K. Adam, A. Scholefield, and M. Vetterli, “Asynchrony increases efficiency: Time encoding of videos and low-rank signals,” IEEE Transactions on Signal Processing, pp. 1–1, 2021.
  * [8] M. Hilton, R. Alexandru, and P. L. Dragotti, “Time encoding using the hyperbolic secant kernel,” in 2020 28th European Signal Processing Conference (EUSIPCO), pp. 2304–2308, IEEE, 2021.
  * [9] H. Namman, S. Mulleti, and Y. C. Eldar, “FRI-TEM: Time encoding sampling of finite-rate-of-innovation signals,” arXiv preprint arXiv:2106.05564, 2021.
  * [10] D. Florescu, F. Krahmer, and A. Bhandari, “Event-driven modulo sampling,” in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5435–5439, IEEE, 2021.
  * [11] R. B. Staszewski and P. T. Balsara, All-digital frequency synthesizer in deep-submicron CMOS.  John Wiley & Sons, 2006.
  * [12] C. Taillefer, “Analog-to-digital conversion via time-mode signal processing,” PhD. dissertation, McGill University, Montréal, Aug. 2007.
  * [13] P. Toledo, R. Rubino, F. Musolino, and P. Crovetti, “Re-thinking analog integrated circuits in digital terms: A new design concept for the iot era,” IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 68, no. 3, pp. 816–822, 2021.
  * [14] M. Miśkowicz, ed., Event-Based Control and Signal Processing.  Embedded Systems, Boca Raton, FL, USA: CRC Press, 2018.
  * [15] A. Lazar and L. T. Tóth, “Perfect recovery and sensitivity analysis of time encoded bandlimited signals,” IEEE Trans. Circ. and Syst.-I, vol. 51, pp. 2060–2073, Oct. 2004.
  * [16] R. Stata, “Operational integrators,” Analog Dialogue, vol. 1, no. 1, pp. 125–129, 1967.
  * [17] S. Norsworthy, S. Norsworthy, R. Schreier, G. Temes, T. Gc, and I. C. . S. Society, Delta-Sigma Data Converters: Theory, Design, and Simulation.  Wiley, 1996.
  * [18] M. Nahvi and J. Edminister, Schaum’s Outline of Electric Circuits, seventh edition.  McGraw-Hill Education, 2017.
  * [19] M. Davies, N. Srinivasa, T.-H. Lin, G. Chinya, Y. Cao, S. H. Choday, G. Dimou, P. Joshi, N. Imam, S. Jain, Y. Liao, C.-K. Lin, A. Lines, R. Liu, D. Mathaikutty, S. McCoy, A. Paul, J. Tse, G. Venkataramanan, Y.-H. Weng, A. Wild, Y. Yang, and H. Wang, “Loihi: A neuromorphic manycore processor with on-chip learning,” IEEE Micro, vol. 38, no. 1, pp. 82–99, 2018.
  * [20] M. Bouvier, A. Valentian, T. Mesquida, F. Rummens, M. Reyboz, E. Vianello, and E. Beigne, “Spiking neural networks hardware implementations and challenges: A survey,” J. Emerg. Technol. Comput. Syst., vol. 15, apr 2019.
  * [21] A. R. Young, M. E. Dean, J. S. Plank, and G. S. Rose, “A review of spiking neuromorphic hardware communication systems,” IEEE Access, vol. 7, pp. 135606–135620, 2019.
  * [22] S. Furber and P. Bogdan, “Spinnaker-a spiking neural network architecture,” 2020\.
  * [23] A. A. Lazar, “Multichannel time encoding with integrate-and-fire neurons,” Neurocomputing, vol. 65-66, pp. 401–407, 2005.  Computational Neuroscience: Trends in Research 2005.
  * [24] D. Wei, “Time-based analog-to-digital converters,” PhD. dissertation, Dept. of Electrical and Computer Eng., Univ. of Florida, 2005.
  * [25] A. A. Lazar, E. K. Simonyi, and L. T. Toth, “An overcomplete stitching algorithm for time decoding machines,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 55, pp. 2619–2630, Oct 2008.
  * [26] T. Strohmer, “Numerical analysis of the nonuniform sampling problem,” J. Comput. Appl. Math., vol. 122, pp. 297–316, Oct. 2000.
  * [27] H. G. Feichtinger, J. C. Príncipe, J. L. Romero, A. Singh Alvarado, and G. A. Velasco, “Approximate reconstruction of bandlimited functions for the integrate and fire sampler,” Advances in Computational Mathematics, vol. 36, pp. 67–78, Jan 2012.
  * [28] W. Bialek, F. Rieke, R. D. R. Van Steveninck, and D. Warland, “Reading a neural code,” Science, vol. 252, no. 5014, pp. 1854–1857, 1991.
  * [29] F. Rieke, Spikes: Exploring the Neural Code.  Bradford book, MIT Press, 1999.
  * [30] C. Eliasmith and C. Anderson, Neural Engineering: Computation, Representation, and Dynamics in Neurobiological Systems.  A Bradford book, MIT Press, 2003.
  * [31] Y. C. Eldar and T. Werther, “General framework for consistent sampling in hilbert spaces,” International Journal of Wavelets, Multiresolution and Information Processing, vol. 3, no. 04, pp. 497–509, 2005.
  * [32] P. L. Combettes, “The foundations of set theoretic estimation,” Proceedings of the IEEE, vol. 81, pp. 182–208, Feb 1993.
  * [33] H. H. Bauschke and J. M. Borwein, “On projection algorithms for solving convex feasibility problems,” SIAM Rev., vol. 38, no. 3, pp. 367–426, 1996.
  * [34] N. Thao and M. Vetterli, “Deterministic analysis of oversampled a/d conversion and decoding improvement based on consistent estimates,” IEEE Transactions on Signal Processing, vol. 42, no. 3, pp. 519–531, 1994.
  * [35] D. Rzepka, M. Miśkowicz, D. Kościelnik, and N. T. Thao, “Reconstruction of signals from level-crossing samples using implicit information,” IEEE Access, vol. 6, pp. 35001–35011, 2018.
  * [36] N. T. Thao and D. Rzepka, “Time encoding of bandlimited signals: Reconstruction by pseudo-inversion and time-varying multiplierless FIR filtering,” IEEE Transactions on Signal Processing, vol. 69, pp. 341–356, 2021.
  * [37] J. L. Yen, “On nonuniform sampling of bandwidth-limited signals,” IRE Trans. Circ. Theory, vol. CT-3, pp. 251–257, Dec. 1956.
  * [38] J. W. Mark and T. Todd, “A nonuniform sampling approach to data compression,” Communications, IEEE Transactions on, vol. 29, pp. 24–32, Jan 1981.
  * [39] N. Sayiner, H. V. Sorensen, and T. R. Viswanathan, “A level-crossing sampling scheme for a/d conversion,” Circuits and Systems II: Analog and Digital Signal Processing, IEEE Transactions on, vol. 43, no. 4, pp. 335–339, 1996.
  * [40] W. Woyczynski, A First Course in Statistics for Signal Analysis.  SpringerLink : Bücher, Birkhäuser Boston, 2010.
  * [41] S. jen Yeh and H. Stark, “Iterative and one-step reconstruction from nonuniform samples by convex projections,” J. Opt. Soc. Am. A, vol. 7, pp. 491–499, Mar 1990.
  * [42] P. Combettes, “Hilbertian convex feasibility problem: Convergence of projection methods,” Applied Mathematics and Optimization, vol. 35, no. 3, pp. 311–330, 1997.
  * [43] H. G. Feichtinger and K. Gröchenig, “Theory and practice of irregular sampling,” in Wavelets: Mathematics and Applications (J. Benedetto, ed.), pp. 318–324, Boca Raton: CRC Press, 1994.
  * [44] D. G. Luenberger, Optimization by vector space methods.  John Wiley & Sons, Inc., New York-London-Sydney, 1969.
  * [45] A. A. Lazar, “Time encoding with an integrate-and-fire neuron with a refractory period,” Neurocomputing, vol. 58-60, pp. 53–58, 2004.  Computational Neuroscience: Trends in Research 2004.