Reconstructing Classes of Non-bandlimited Signals from Time Encoded Information

TL;DR

This paper explores reconstructing non-bandlimited signals from time-encoded samples using filtering and non-uniform sampling, introducing new algorithms and conditions for perfect recovery, even with noisy data.

Contribution

It introduces a novel framework for reconstructing non-bandlimited signals from time encoding, focusing on exponential and polynomial spline filters, with algorithms for perfect recovery.

Findings

Sampling by timing is equivalent to non-uniform sampling.

Fundamental properties of filters are locally preserved in non-uniform sampling.

Proposed algorithms achieve perfect reconstruction under certain conditions.

Abstract

We investigate time encoding as an alternative method to classical sampling, and address the problem of reconstructing classes of non-bandlimited signals from time-based samples. We consider a sampling mechanism based on first filtering the input, before obtaining the timing information using a time encoding machine. Within this framework, we show that sampling by timing is equivalent to a non-uniform sampling problem, where the reconstruction of the input depends on the characteristics of the filter and on its non-uniform shifts. The classes of filters we focus on are exponential and polynomial splines, and we show that their fundamental properties are locally preserved in the context of non-uniform sampling. Leveraging these properties, we then derive sufficient conditions and propose novel algorithms for perfect reconstruction of classes of non-bandlimited signals such as: streams of Diracs, sequences of pulses and piecewise constant signals. Next, we extend these methods to operate with arbitrary filters, and also present simulation results on synthetic noisy data.

Figures

Tables

Table 1. TABLE I: Choice of hyperparameters of the crossing TEM for reconstructing bursts of K 𝐾 K Diracs.

$K$	$M$	$P$	$L$	$f_s^{min}$	Ideal sample density (samples/s)	Sample density (samples/s)
1	1	$2$	$2$ sec	$\frac{5}{2L}=1.25$ Hz	$\frac{2}{L+ϵ}$	$2f_s=2.5$
2	2	$2$	$2$ sec	$\frac{7}{2L}=1.75$ Hz	$\frac{4}{L+ϵ}$	$2M{f}_s=7$
4	4	$2$	$2$ sec	$\frac{7}{2L}=1.75$ Hz	$\frac{8}{L+ϵ}$	$2M{f}_s=14$

Table 2. TABLE II: Choice of hyperparameters for estimating a burst of 2 2 2 Diracs, using an M 𝑀 M -channel integrate-and-fire system.

$M$	$P$	$L$	$C_T^{max}$	${\omega }_0$	${\Delta }_b$	Average number of samples per burst
$2$	$2$	$2$	$0.1218$	$-\frac{\pi }{3}$	$0.2$	44.8
$2$	$2$	$2$	$0.0947$	$-\frac{\pi }{3}$	$0.3$	56.21
$1$	$4$	$4$	$0.0114$	$-\frac{\pi }{3}$	$0.2$	202.24
$2$	$2$	$2$	$0.1286$	$-\frac{\pi }{2}$	$0.2$	43.25
$2$	$4$	$4$	$0.0133$	$-\frac{\pi }{3}$	$0.2$	413.3

Table 3. TABLE III: Effect of noise on the estimation of a piecewise constant signal, from spikes obtained using the integrate-and-fire TEM. The error ϵ t subscript italic-ϵ 𝑡 \epsilon_{t} is the average absolute difference between the true and estimated locations, ϵ A subscript italic-ϵ 𝐴 \epsilon_{A} is the relative error of the estimated amplitudes of the input discontinuities and SER is the signal-to-error-ratio, for amplitude estimation.

SNR(dB)	$\sigma$	${ϵ}_t$	${ϵ}_A$	SER(dB)
$43.33$	0.01	$2.61\times {10}^{-4}$	$6.21\times {10}^{-5}$	$43.11$
$29.63$	0.05	$0.0015$	$2.1509\times {10}^{-4}$	$24.59$
$23.38$	0.1	$0.0042$	$0.0026$	$22.83$

Table 4. TABLE IV: Probability of perfect reconstruction of a burst of 2 Diracs, with random Gaussian amplitudes.

${\Delta }_b={\tau }_2-{\tau }_1$	$C_T$	P{perfect estimation}
$0.05$	$0.01$	$0.950$
$0.1$	$0.01$	$0.959$
$0.5$	$0.001$	$0.959$
$1.5$	$0.001$	$0.978$

Equations

$$y_n=y(t_n)=\int x(\tau )\phi (\tau -t_n)d\tau =⟨x(t),\phi (t-t_n)⟩.$$

$$y_n=y(t_n)=\pm C_T={\int }_{t_{n-1}}^{t_n}f(\tau )d\tau ,$$

$$f(t)=\int x(\alpha )\phi (\alpha -t)d\alpha ,\text{ for }t\in [t_{n-1},t_n].$$

$$y_1=y(t_1)=\pm C_T={\int }_{{\tau }_1}^{t_1}f(\tau )d\tau .$$

$$y(t_n)=⟨x(t),(\phi \ast q_{{\theta }_n})(t-t_{n-1})⟩,$$

$$q_{{\theta }_n}(t)=\{\begin{array}{ll}1,& 0\le t\le {\theta }_n,\\ 0,& otherwise.\end{array}$$

$$\begin{array}{rl}y(t_n)& ={\int }_{t_{n-1}}^{t_n}f(\tau )d\tau ={\int }_{t_{n-1}}^{t_n}{\int }_{-\infty }^{\infty }x(t)\phi (t-\tau )dtd\tau \\ & \stackrel{\text{mbox}(a)}{=}{\int }_{-\infty }^{\infty }x(t){\int }_{t_{n-1}}^{t_n}\phi (t-\tau )d\tau dt\\ & \stackrel{\text{mbox}(b)}{=}{\int }_{-\infty }^{\infty }x(t){\int }_{t-t_n}^{t-t_{n-1}}\phi (\tau )d\tau dt\\ & \stackrel{\text{mbox}(c)}{=}{\int }_{-\infty }^{\infty }x(t){\int }_{t-t_n}^{t-t_{n-1}}\phi (\tau )q_{{\theta }_n}(t-t_{n-1}-\tau )d\tau dt\\ & \stackrel{\text{mbox}(d)}{=}{\int }_{-\infty }^{\infty }x(t)(\phi \ast q_{{\theta }_n})(t-t_{n-1})dt\\ & =⟨x(t),(\phi \ast q_{{\theta }_n})(t-t_{n-1})⟩.\end{array}$$

$$y(t_1)\stackrel{\text{mbox}(a)}{=}{\int }_{{\tau }_1}^{t_1}f(\tau )d\tau =⟨x(t),(\phi \ast q_{{\theta }_1})(t-{\tau }_1)⟩,$$

$${\beta }_P(t)=\underset{P+1\text{ times}}{\underbrace{{\beta }_0(t)\ast {\beta }_0(t)....\ast {\beta }_0(t)}},$$

$${\beta }_0(t)=\{\begin{array}{ll}1,& -1\le t\le 0,\\ 0,& \text{otherwise}.\end{array}$$

$$\sum _{n\in \mathbb{Z}}{c}_{m,n}{\beta }_P(t-n)=t^m,$$

$${\beta }_1(t)=\{\begin{array}{ll}-t,& -1\le t\le 0,\\ 2+t,& -2\le t<-1,\\ 0,& \text{otherwise}.\end{array}$$

$$\sum _{n=0}^{N-1}{c}_{m,n}^I{\beta }_1(t-t_n)=t^m,$$

$${\phi }_1(t)=\{\begin{array}{ll}{e}^{-{\alpha }_{0}t},& -1\le t\le 0,\\ 0,& \text{otherwise}.\end{array}$$

$${\phi }_P(t)={\phi }_{{\alpha }_0}(t)\ast {\phi }_{{\alpha }_1}(t)....\ast {\phi }_{{\alpha }_{P-1}}(t).$$

$$\sum _{n\in \mathbb{Z}}{c}_{m,n}\phi (t-n)=e^{-{\alpha }_{m}t},$$

$${\phi }_2(t)=\{\begin{array}{ll}\frac{e^{c_1-c_0}}{c_1-c_0}{e}^{-{\alpha }_{0}t}+\frac{e^{-c_1+c_0}}{c_0-c_1}{e}^{-{\alpha }_{1}t},& -L\le t<-\frac{L}{2},\\ \frac{1}{c_0-c_1}{e}^{-{\alpha }_{0}t}+\frac{1}{c_1-c_0}{e}^{-{\alpha }_{1}t},& -\frac{L}{2}\le t\le 0,\\ 0,& \text{otherwise},\end{array}$$

$$\sum _{n=0}^{N-1}{c}_{m,n}^I{\phi }_2(t-t_n)=e^{-{\alpha }_{m}t},$$

$$({\phi }_P\ast g)(t)=\sum _{m=0}^{P-1}{a}_m{G}_m{e}^{{\alpha }_{m}t},$$

$$\sum _{n=0}^{N-1}{c}_{m,n}^I({\phi }_P\ast g)(t-t_n)=e^{-{\alpha }_{m}t},$$

$$\text{vspace}-0.5emx(t)=x_1\delta (t-{\tau }_1).$$

$$\sum _{n=1}^2{c}_{m,n}^I\phi (t-t_n)=e^{j{\omega }_{m}t},\text{for }m\in \{0,1\}.\text{vspace}-0.5em$$

$$\text{vspace}-0.5em\begin{array}{rl}{s}_m& =\sum _{n=1}^2{c}_{m,n}^{I}y(t_n)\stackrel{\text{mbox}(a)}{=}\sum _{n=1}^2{c}_{m,n}^I⟨x(t),\phi (t-t_n)⟩\\ & \stackrel{\text{mbox}(b)}{=}{\int }_{-\infty }^{\infty }x(t)\sum _{n=1}^2{c}_{m,n}^I\phi (t-t_n)dt\\ & \stackrel{\text{mbox}(c)}{=}{\int }_{-\infty }^{\infty }{x}_1\delta (t-{\tau }_1)\sum _{n=1}^2{c}_{m,n}^I\phi (t-t_n)dt\\ & \stackrel{\text{mbox}(d)}{=}{\int }_I{x}_1\delta (t-{\tau }_1)e^{j{\omega }_{m}t}dt=x_1{e}^{j{\omega }_m{\tau }_1}=b_1{u}_1^m,\end{array}\text{vspace}-1.5em$$

$$\text{vspace}-0.7emx(t)=\sum _k{x}_k\delta (t-{\tau }_k).\text{vspace}-0.1em$$

$$s_m=c_{m,n}^{I}y(t_n)+c_{m,n+1}^{I}y(t_{n+1})=x_2{e}^{j{\omega }_m{\tau }_2}.$$

$$x(t)=\sum _b\sum _{k=1}^K{x}_{b,k}\delta (t-{\tau }_{b,k}),$$

$$\text{vspace}-0.5em0<C_T<\frac{A_{min}}{3}{\int }_0^{\frac{L}{2}}\phi (-t)dt,$$

$${\int }_0^{\frac{L}{2}}\phi (-t)dt={\int }_{{\tau }_1}^{{\tau }_1+\frac{L}{2}}\phi ({\tau }_1-t)dt\stackrel{\text{mbox}(a)}{=}\frac{1}{{\omega }_0^2}[1-\cos ({\omega }_0\frac{L}{2})],$$

$$3C_T<A_{\min }{\int }_0^{\frac{L}{2}}\phi (-t)\stackrel{\text{mbox}(b)}{<}{\int }_{{\tau }_1}^{{\tau }_1+\frac{L}{2}}{x}_1\phi ({\tau }_1-t)dt,$$

$$\begin{array}{rl}3C_T& ={\int }_{{\tau }_1}^{t_3}f(t)dt\stackrel{\text{mbox}(c)}{=}{\int }_{{\tau }_1}^{t_3}{x}_1\phi ({\tau }_1-t)dt,\end{array}$$

Full text

I Introduction

Sampling plays a fundamental role in signal processing and communications, achieving the conversion of continuous time phenomena into discrete sequences [1]. From the Whittaker-Shannon theorem [2], to recent theories in compressed sensing [3, 4], super-resolution [5] and finite rate of innovation [6, 7, 8, 9, 10], sampling theory has provided precise answers on when a faithful conversion of a continuous waveform into a discrete sequence is possible. These methods are generally based on recording the signal amplitude at specified times, which lead to uniform sampling if the samples are evenly spaced, and non-uniform sampling otherwise.

In this paper, we concentrate on an alternative method to classical sampling, which encodes the input into a sequence of non-uniformly spaced time events or spikes. In other words, rather than recording the value of the signal at preset times, one records the instants when the signal crosses a pre-defined threshold or triggers a pre-defined event. ††Some of the work in this paper was, in part, presented at the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Brighton, UK, May 2019 [11], and the International Conference on Sampling Theory and Applications (SampTA), Bordeaux, France, July 2019 [12].Acquisition models inspired by this mechanism include zero-crossing detectors [13], delta-modulation schemes[14], as well as the time encoding machine (TEM) introduced in [15]. This latter model is of particular interest, as it mimics the integrate-and-fire mechanism of neurons in the human brain. Biological neurons use time encoding to represent sensory information as action potentials [16, 17, 18], which allows them to process information very efficiently. In the same manner, sampling inspired by the brain could lead to very simple and highly efficient devices, ranging from analog to digital converters [15], to neuromorphic computing or event-based vision sensors, which record only changes in the input intensity, leading to low power consumption and fewer storage requirements [19].

At the same time, time-encoding methods extend theories of traditional sampling, and this makes this topic intriguing also from a research perspective. Within the study of time encoding, the key problem that arises is to find methods to retrieve the input signal from its timing information, and hence the key questions to pursue are the following. 1) Is time encoding invertible, and which classes of signals can be uniquely represented using timing information? 2) What algorithms allow perfect retrieval of these signals from their time-encoded samples?

To address these questions, several authors have provided ways to sample and reconstruct bandlimited signals [20, 21, 22, 23, 24, 25]. These initial results on time-encoding machines have also been extended to functions that belong to shift-invariant spaces [26, 27], typically by connecting time encoding with the problem of non-uniform sampling [28, 29, 30]. Time encoding theory has also been generalized to the case of non-bandlimited signals in [31], however in the context of studying the dynamics of populations of neurons, by leveraging stochastic assumptions on the firing parameters.

In this paper, we show that it is possible to perfectly reconstruct particular classes of continuous-time signals which are neither bandlimited nor belong to shift-invariant subspaces, from samples obtained using a time encoding mechanism. The signals we focus on are infinite streams of Diracs, sequences of pulses, as well as piecewise constant signals. Sampling and reconstructing pulses is of significant relevance to many real-world applications. For example, time-of-flight cameras probe the 3D scene with pulses of light and reconstruct the scene by measuring their round trip time. In applications which require reduced computational power and speed, e.g. robots mapping their surroundings, time-of-flight technology may benefit from a time encoding framework which would significantly lower the sampling rate. Signals consisting of a stream of pulses appear in many other applications, including: ultrawideband communications [32], ECG acquisition and compression [33], radio-astronomy [34], image processing[35], ultrasound imaging [8] and processing of neuronal signals [36].

At the same time, time encoding principles have already been integrated in bio-inspired technologies such as dynamic vision sensors (DVS) [19], which have many real-world applications, ranging from robotics to autonomous driving as well as low-power surveillance. In a DVS camera, each pixel only records changes in the input at the time instants they occur, by taking a time derivative of the signal. Hence, at the local pixel-level, this is equivalent to time encoding of piecewise constant signals, which is studied in this paper.

Motivated by these real-world applications, the time encoding strategy we propose is based on filtering the input signal before extracting the timing information using a crossing or an integrate-and-fire TEM. The filter may be used to reduce noise, or may model the distortion introduced by the acquisition device, for example the optics in a time-of-flight scanner or the photoreceptors in a DVS camera. In order to develop a framework for exact reconstruction, we initially focus on two classes of compact-support filters (sampling kernels): exponential and polynomials splines. Please note that exponential splines are very useful since they can be used to model any convolution operator with rational transfer function as for example, simple RC circuits [37, 6]. Our first main contribution is to prove that exponential (polynomial) splines locally preserve their exponential (polynomial) reproducing properties in the context of time-based sampling. Specifically, we show that within intervals where there are no knots of at least $N$ non-uniformly shifted kernels, we can locally reproduce exponentials (polynomials) of degree $N$. The second aspect of our contribution is to leverage these properties to address the problem of reconstructing some classes of non-bandlimited signals from timing information. We initially develop our reconstruction framework for the case of one Dirac, where we show how a linear combination of its non-uniform samples leads to a sequence of signal moments, which can then be annihilated using Prony’s method [38], in order to retrieve the free parameters of the input. Furthermore, we extend this method to reconstruct infinite streams and bursts of Diracs, sequences of pulses as well as piecewise constant signals, for which we can achieve local reconstruction given the compact support of the filter. Finally, we depart from the ideal case, and present a universal reconstruction strategy that works with timing-based samples taken by arbitrary kernels.

This paper is organized as follows. In Section II-A, we describe the principles of time encoding, with two exemplary cases. Then, in Section II-B we show that sampling kernels which reproduce exponentials or polynomials preserve this property locally, when sampling is based on timing information. Furthermore, in Section III we present methods for the reconstruction of non-bandlimited signals from their timing information obtained using a crossing TEM. We first propose a method for estimation of a single Dirac, and extend this to retrieve streams of Diracs and bursts of Diracs. Then, in Section IV we demonstrate the perfect retrieval of classes of non-bandlimited signals from timing information, obtained using an integrate-and-fire TEM. These estimation methods are then extended in Section V to the case of arbitrary sampling kernels. Here we also present results for the case of noisy signals. Finally, we highlight the high efficiency of sampling based on timing information in Section VI, and present concluding remarks in Section VII. Please note that the code to reproduce our simulations is available online [39].

II Time Encoding Mechanisms

II-A Acquisition Models

In this section, we introduce the time encoding machines considered in this paper: the crossing TEM and the integrate-and-fire TEM. Specifically, we show how these TEMs map a real signal $x(t)$ to a strictly increasing sequence of times $\{t_{n}\}$ [27]. We also show that although no measure of the amplitude of the signal is recorded, time encoding is equivalent to a non-uniform sampling problem.

II-A1 Crossing Time Encoding Machine

The crossing time encoding strategy is inspired by the A/D conversion scheme in e.g. [40, 27], and is depicted in Fig. 1. It consists of a compact-support filter $\varphi(-t)$, and a comparator with a sinusoidal reference $g(t)$. The output of the acquisition device is the sequence $\{t_{n}\}$, corresponding to the time instants when the filtered input signal crosses the reference, i.e. when $y(t_{n})-g(t_{n})=0$. Moreover, since the shape of the test function $g(t)$ is known, we can retrieve the amplitudes of the output samples, given by $y_{n}=y(t_{n})=g(t_{n})$. Hence, decoding the input signal is equivalent to a non-uniform sampling problem, where we aim to reconstruct $x(t)$ from the non-uniform samples given by:

[TABLE]

In Fig. 2 we depict the time encoded information of an input signal of 3 Diracs, obtained using the TEM in Fig. 1.

II-A2 Time Encoding based on an Integrate-and-fire System

The operating principle of this time encoding strategy is similar to the one in [20], and is depicted in Fig. 3. The signal is first filtered with a compact-support filter with impulse response $\varphi(-t)$, before being passed to an integrator. When the output of the integrator reaches the positive trigger mark $C_{T}$, the time encoding machine outputs a spike and the integrated signal $y(t)$ is reset to zero. Similarly, a spike is generated and $y(t)$ resets to zero, when the integrator reaches the negative trigger mark $-C_{T}$. The time instants when the integrator reaches the threshold $\pm C_{T}$ are recorded in the sequence $\{t_{n}\}$. Then, we can compute the output sample $y(t_{n})$ at each spike $t_{n}$ as:

[TABLE]

where $n\geq 2$ and $f(t)$ is defined as:

[TABLE]

Similarly, assuming that the input signal $x(t)=0$, for $t<\tau_{1}$, and that the filter $\varphi(-t)$ is causal, then the first output sample is given by:

[TABLE]

Hence, time encoding with an integrate-and-fire model is equivalent to a non-uniform sampling problem, where we aim to estimate the input $x(t)$ from the non-uniform samples $y(t_{n})$. In Fig. 4 we depict the time encoding of an input signal, obtained using the device in Fig. 3, for $C_{T}=0.15$.

Furthermore, leveraging the results in [41], we can show that the non-uniform output samples we obtain using the acquisition model in Fig. 3 are the same as those obtained by filtering the input with the modified kernel $(\varphi*q_{\theta_{n}})(t)$:

[TABLE]

where $\theta_{n}=t_{n}-t_{n-1}$ and $q_{\theta_{n}}(t)$ is defined as:

[TABLE]

We can prove Eq. (5) by re-writing Eq. (2) as follows:

[TABLE]

In the derivations above, $(a)$ holds since we assume both the input $x(t)$ and the filter $\varphi(t)$ have compact support, and $(b)$ follows from a change of variable. Moreover, $(c)$ follows from the fact that $q_{\theta_{n}}(t-t_{n-1}-\tau)=1$ for $\tau\in[t-t_{n},t-t_{n-1}]$ and $(d)$ holds since $q_{\theta_{n}}(t-t_{n-1}-\tau)=0$ for $\tau\notin[t-t_{n},t-t_{n-1}]$, as defined in Eq. (6).

Finally, the first output sample can be computed as:

[TABLE]

where $\theta_{1}=t_{1}-\tau_{1}$, and $(a)$ follows from Eq. (4).

We conclude this subsection by making the following remark. We observe that from the timing sequence $\{t_{n}\}$, we can either recover $y(t_{n})=\langle x(t),\varphi(t-t_{n})\rangle$ for the case of the crossing TEM or $y(t_{n})=\langle x(t),(\varphi*q_{\theta_{n}})(t-t_{n-1})\rangle$ for the integrate-and-fire model. This means that in both cases, the reconstruction of $x(t)$ will depend on the proper choice of the sampling kernel $\varphi(t)$ and on its non-uniform shifts $\varphi(t-t_{n})$.

In what follows we focus on two families of kernels, polynomial and exponential splines [6, 42, 43], and show that some of their fundamental properties are preserved in the case of non-uniform shifts.

II-B Sampling Kernels

The sampling kernels $\varphi(t)$, that we consider in this paper are all anti-causal since they are the time reversed versions of causal filters.

II-B1 Polynomial splines

A B-spline $\beta_{P}(t)$ of order $P$ is computed as the $(P+1)$-fold convolution of the box function $\beta_{0}(t)$ [42]:

[TABLE]

where the anti-causal version of $\beta_{0}(t)$ is defined as:

[TABLE]

The B-spline of order $P$ satisfies the Strang-Fix conditions [44] and hence, together with its uniform shifts, it can reproduce polynomials of maximum degree $P$:

[TABLE]

where $m\in\{0,1,...,P\}$, and for a proper choice of the coefficients $c_{m,n}$.

For instance, the first-order B-spline satisfies Eq. (9) for $P=1$, which means it can reproduce constant and linear polynomials, and is defined as:

[TABLE]

The first order B-spline has two continuous regions, each of which is a linear polynomial: $\beta_{1}^{A}(t)=-t$, for $t\in(-1,0)$ and $\beta_{1}^{B}(t)=2+t$, for $t\in(-2,-1)$. Using this observation, it is possible to show that the first-order B-spline, together with its non-uniformly shifted versions can locally reproduce polynomials of maximum degree $1$. In other words, it is possible to prove that within a time interval $I$ where the shifted kernels $\beta_{1}(t-t_{n})$ have no knots, the following equation holds:

[TABLE]

where $N\geq 2$, $m\in\{0,1\}$, $t\in I$ and $\{t_{n}\}$ are non-uniform.

The proof can be outlined by setting $N=2$ for simplicity. Then, let $I$ be an interval where there are no knots of $\beta_{1}(t-t_{0})$ and $\beta_{1}(t-t_{1})$, with $I\subset(t_{1}-1,t_{0})$. Furthermore, let $v_{0}(t)=\beta_{1}(t-t_{0})=-t+t_{0}$ for $t\in I$ and $v_{1}(t)=\beta_{1}(t-t_{1})=-t+t_{1}$ for $t\in I$. In the vector space of linear polynomials in $I$ which is a two-dimensional space, the elements $v_{0}(t)$ and $v_{1}(t)$ are linearly independent and so form a basis of the space, provided $t_{0}\neq t_{1}$. Therefore, using a linear combination of the two functions, we can uniquely represent any vector in this space, including the vector $t$. In other words, we can determine the unique coefficients $c_{1,0}^{I}=\frac{t_{1}}{t_{0}-t_{1}}$ and $c_{1,1}^{I}=\frac{t_{0}}{t_{1}-t_{0}}$ that ensure $c_{1,0}^{I}v_{0}(t)+c_{1,1}^{I}v_{1}(t)=t$, for $t\in I$. Similarly, we find the unique coefficients $c_{0,0}^{I}=\frac{1}{t_{0}-t_{1}}$ and $c_{0,1}^{I}=\frac{1}{t_{1}-t_{0}}$ such that $c_{0,0}^{I}v_{0}(t)+c_{0,1}^{I}v_{1}(t)=1$, for $t\in I$. Hence, Eq. (10) is satisfied in the knot-free interval $I$ for $N=2$.

In the same manner, one can show that reproduction of constant and linear polynomials is achieved on any interval spanned by knot-free regions of at least two non-uniformly shifted B-splines. Lastly but importantly, for different knot-free intervals, the solution to Eq. (10) differs, and this fact is highlighted in Fig. 5. Here, we depict two non-uniform shifts of the first-order B-spline, namely $\beta_{1}(t-2)$ and $\beta_{1}(t-2.625)$. The shifted kernel $\beta_{1}(t-2)$ has knots at $t=0$, $t=1$ and $t=2$, whilst $\beta_{1}(t-2.625)$ has knots at $t=0.625$, $t=1.625$ and $t=2.625$. As a result, reproduction of polynomials is possible within the knot-free regions $I_{1}=(0.625,1)$ and $I_{2}=(1,1.625)$, however with a different linear combination of the B-splines overlapping these regions, i.e. with $c_{m,n}^{I_{1}}\neq c_{m,n}^{I_{2}}$.

One can extend this result to the case of higher order polynomials by using B-splines of order $P>1$. This is due to the fact that polynomial splines are piecewise polynomial functions of degree $P$. Hence, in any interval $I$ that contains $P+1$ knot-free shifted versions of splines, it is possible to reproduce polynomials up to degree $P$.

II-B2 Exponential splines

The anti-causal version of the E-spline of first-order is defined as:

[TABLE]

where $\alpha_{0}$ can be either real or complex.

As with polynomial splines, E-splines of order $P$ are obtained from the convolution of first-order E-splines [43]:

[TABLE]

An E-spline of order $P$ has compact support and can reproduce $P$ different exponentials of the form $e^{-\alpha_{m}t}$ [43]:

[TABLE]

where $m=0,1,...,P$, and for a suitable choice of the coefficients $c_{m,n}$.

For example, the E-spline of order $P=2$ of support of arbitrary length $L$ is defined as:

[TABLE]

where $\alpha_{i}\in\mathbb{C}$ (if $\mathchar 572\relax\{\alpha_{i}\}=0$ then $\varphi_{2}(t)\in\mathbb{R}$), and where $c_{i}=\alpha_{i}\frac{L}{2}$ for $i=0,1$ in order to ensure continuity of $\varphi_{2}(t)$. Throughout the remainder of the paper, we assume for simplicity that $L=2$.

The second-order E-spline can reproduce the exponentials $e^{-\alpha_{0}t}$ and $e^{-\alpha_{1}t}$. In fact, we notice that within each of its knot-free regions, the function $\varphi_{2}(t)$ can be expressed as a linear combination of the exponentials $e^{-\alpha_{0}t}$ and $e^{-\alpha_{1}t}$. This observation helps us prove that within any time interval $I$ which contains knot-free regions of non-uniformly shifted first-order E-splines, we can reproduce two exponentials:

[TABLE]

where $N\geq 2$, $m\in\{0,1\}$, $t\in I$ and $\{t_{n}\}$ are non-uniform.

For example, let $I$ be an interval which contains knot-free regions of $\varphi_{2}(t-t_{0})$ and $\varphi_{2}(t-t_{1})$, with $I\subset(t_{1}-L,t_{0}-\frac{L}{2})$. Moreover, let $v_{0}(t)=\varphi_{2}(t-t_{0})$ for $t\in I$ and $v_{1}(t)=\varphi_{2}(t-t_{1})$ for $t\in I$. The elements $v_{0}(t)$ and $v_{1}(t)$ are linear combinations of $e^{-\alpha_{0}t}$ and $e^{-\alpha_{1}t}$, and therefore belong to the vector space spanned by these two exponentials. Moreover, $v_{0}(t)$ and $v_{1}(t)$ are linearly independent and so, form a basis of that vector space, since $t_{1}\neq t_{0}$. Hence, using a linear combination of $v_{0}$ and $v_{1}$, we can uniquely represent any vector in this space, including $e^{-\alpha_{0}t}$ and $e^{-\alpha_{1}t}$. Therefore, in the interval $I$ where there are no knots, we can find unique coefficients $c_{m,0}^{I}$ and $c_{m,1}^{I}$ such that Eq. (13) holds for $m\in\{0,1\}$.

Similarly, reproduction of two different exponentials is possible on any time interval spanned by knot-free regions of at least two shifted E-splines. Note that for different intervals $I_{1}$ and $I_{2}$, the solution to Eq. (13) differs, i.e. $c_{m,n}^{I_{1}}\neq c_{m,n}^{I_{2}}$. This is highlighted in Fig. 6, where exponential reproduction is possible in the regions $I_{1}$ and $I_{2}$, but using a different linear combination of the E-splines that overlap these regions.

By using the same argument we can prove similar results for the general case of an E-spline of order $P$ and support of length $L$ which can reproduce $P$ different exponentials. Specifically, within an interval $I$ containing knot-free regions of at least $P$ non-uniformly shifted E-splines, we can reproduce $P$ different exponentials, such that Eq. (13) holds for $N\geq P$ and $m\in\{0,1,...,P-1\}$. This is due to the fact that any knot-free interval of an E-spline of order $P$ is a linear combination of $P$ different exponentials.

Finally, let us consider the kernel $(\varphi_{P}*g)(t)$, where $\varphi_{P}(t)$ is a $P$-order E-spline which can reproduce the exponentials $e^{\alpha_{m}t}$, for $m=0,1,...,P-1$. Furthermore, let us assume that $g(t)$ has compact support $[-\epsilon,\epsilon]$. The support of $\varphi_{P}(t)$ is $[-L,0]$ and its knots are located at instants $(-L+n\frac{L}{P})$ with $n\in\mathbb{N}$. Then, in the knot-free interval $(\frac{-L}{P},0)$ we can compactly represent $\varphi_{P}(t)=\sum_{m=0}^{P-1}a_{m}e^{\alpha_{m}t}$, for some coefficients $a_{m}$.

If the length of the support of $g(t)$ satisfies $2\epsilon\leq\frac{L}{P}$ and $\int g(t)e^{-\alpha_{m}t}dt$ exists, then $(\varphi_{P}*g)(t)$ is given by:

[TABLE]

where $t\in(\epsilon-\frac{L}{P},-\epsilon)$ and $G_{m}=\int_{-\epsilon}^{\epsilon}g(t)e^{-\alpha_{m}t}dt$.

Therefore, in the interval $(\epsilon-\frac{L}{P},-\epsilon)$, $(\varphi_{P}*g)(t)$ is a linear combination of $P$ exponentials. As a result, within $I=(t_{N-1}+\epsilon-\frac{L}{P},t_{1}-\epsilon)$, $(\varphi_{P}*g)(t)$ and its non-uniform shifts can reproduce $P$ exponentials, as follows:

[TABLE]

where $N\geq P$, and $m\in\{0,1,...,P-1\}$.

III Perfect Recovery of Signals from Timing Information obtained with a Crossing TEM

In the previous section, we showed how time encoding maps the input signal to a sequence of non-uniform samples, which depend on the signal and non-uniform shifts of the sampling kernel. In what follows we assume that the sampling kernel $\varphi(t)$ is a second-order exponential reproducing spline, such that a linear combination of its non-uniformly shifted versions can reproduce two different exponentials, as described in Section II-B2. Moreover, $\varphi(t)$ has compact support of length $L$, with $\varphi(t)=0$ for $t\notin[-L,0]$ and the two frequencies that this kernel can reproduce are $\alpha_{0}=j\omega_{0}$ and $\alpha_{1}=-\alpha_{0}$, which ensures that $\varphi(t)$ is a real-valued function.

Under these assumptions, we study the problem of reconstructing different classes of non-bandlimited signals, from timing information obtained using the crossing TEM in Fig. 1. Specifically we present a method for estimation of an input Dirac. Here we show that two output spikes are sufficient to retrieve the input, provided they are located suitably close to the Dirac, which is guaranteed by imposing conditions on the frequency and amplitude of the comparator’s sinusoidal reference signal. We then extend this to retrieval of streams of Diracs and bursts of Diracs. While the reconstruction method proposed to retrieve one Dirac might not be unique, it has the advantage that it naturally generalizes to multiple Diracs. We note that similar results could be proved using polynomial splines, but we omit these proofs to keep the focus of the paper on E-splines.

III-A Estimation of an Input Dirac

Let us consider a single input Dirac of the form:

[TABLE]

Proposition 1**.**

Let the sampling kernel $\varphi(t)$ be a second-order E-spline of support of length $L$, defined as in Eq. (12), with $\omega_{1}=-\omega_{0}$ and $0<\omega_{0}\leq\frac{\pi}{L}$. The filter $\varphi(t)$ and its non-uniform shifts can reproduce the exponentials $e^{j\omega_{0}t}$ and $e^{j\omega_{1}t}$ as in Eq. (13). In addition, suppose that the reference signal $g(t)=A\cos(w_{s}t)$ of the comparator in Fig. 1 has amplitude $A>|x_{1}|$ and period $T_{s}<\frac{2L}{5}$. Then, the timing information $\{t_{1},t_{2},...,t_{N}\}$ provided by the comparator TEM is a sufficient representation of an input Dirac as in Eq. (16).

Proof.

From the timing information $\{t_{1},t_{2}\}$, we can retrieve the non-uniform output samples $y(t_{1})$ and $y(t_{2})$, as described in Eq. (1). In what follows we show that we can find a linear combination of the samples $y(t_{1})$ and $y(t_{2})$ to get $c_{m,1}y(t_{1})+c_{m,2}y(t_{2})=x_{1}e^{j\omega_{m}\tau_{1}}$, for $m=0,1$, from which we can retrieve the input parameters $x_{1}$ and $\tau_{1}$.

For simplicity, suppose that the amplitude of the input Dirac satisfies $x_{1}>0$. In addition, the hypothesis that $\varphi(t)$ reproduces $e^{\pm j\omega_{0}t}$ with $0<\omega_{0}\leq\frac{\pi}{L}$ means that $0\leq\varphi(t)<1$, for $t\in[-\frac{L}{2},0]$. Then, since $0<x_{1}<A$, the output $y(t)=x_{1}\varphi(\tau_{1}-t)$ of the crossing TEM satisfies $0\leq y(t)<A=\max(g(t))$, for $t\in[\tau_{1},\tau_{1}+\frac{L}{2}]$. Since we assume $\frac{5T_{s}}{4}<\frac{L}{2}$, this means that $0\leq y(t)<A=\max(g(t))$, for $t\in[\tau_{1},\tau_{1}+\frac{5T_{s}}{4}]$.

Let us then define the continuous function $h(t)=g(t)-y(t)$. Using Bolzano’s intermediate value theorem [45] and the fact that $0\leq y(t)<\max(g(t))$, we show that within the interval $(\tau_{1},\tau_{1}+\frac{5T_{s}}{4}]$, the signal $h(t)$ crosses zero at least twice. In other words, $\exists t_{1},t_{2}\in(\tau_{1},\tau_{1}+\frac{5T_{s}}{4}]$ such that $h(t_{1})=h(t_{2})=0$. For example, if we assume $h(\tau_{1})=g(\tau_{1})>0$, then $g(\tau_{1}+\frac{T_{s}}{2})<0$ and since $y(t)\geq 0$, we get $h(\tau_{1}+\frac{T_{s}}{2})=g(\tau_{1}+\frac{T_{s}}{2})-y(\tau_{1}+\frac{T_{s}}{2})<0$. Then, Bolzano’s intermediate value theorem states that $\exists t_{1}\in(\tau_{1},\tau_{1}+\frac{T_{s}}{2}]$ such that $h(t_{1})=0$.

Using the same argument one can then show that $\exists t_{2}\in(\tau_{1}+\frac{T_{s}}{2},\tau_{1}+\frac{5T_{s}}{4}]$ such that $h(t_{2})=0$. This follows from the assumption that $g(\tau_{1})>0$, which implies that $\exists\epsilon\in[0,\frac{T_{s}}{2}]$ such that $g(\tau_{1}+\frac{3T_{s}}{4}+\epsilon)=\cos(\tau_{1}+\frac{3T_{s}}{4}+\epsilon)=\max(g(t))=A$, as highlighted in Fig. 7. At the same time, we showed that $0\leq y(t)<A$ for $t\in[\tau_{1},\tau_{1}+\frac{5T_{s}}{4}]$ and hence, we get $h(\tau_{1}+\frac{3T_{s}}{4}+\epsilon)=g(\tau_{1}+\frac{3T_{s}}{4}+\epsilon)-y(\tau_{1}+\frac{3T_{s}}{4}+\epsilon)>0$. Since $h(\tau_{1}+\frac{T_{s}}{2})<0$ and $h(\tau_{1}+\frac{3T_{s}}{4}+\epsilon)>0$, Bolzano’s intermediate value theorem guarantees that $\exists t_{2}\in(\tau_{1}+\frac{T_{s}}{2},\tau_{1}+\frac{3T_{s}}{4}+\epsilon]$ such that $h(t_{2})=0$, for $\epsilon\in[0,\frac{T_{s}}{2}]$. Therefore, at the maximum value of $\epsilon$, we have proved that the second output spike satisfies $t_{2}\in[\tau_{1}+\frac{T_{s}}{2},\tau_{1}+\frac{5T_{s}}{4}]$.

Hence, since we assume $\frac{L}{2}>\frac{5T_{s}}{4}$, we obtain the inequality $t_{2}\leq\tau_{1}+\frac{5T_{s}}{4}<\tau_{1}+\frac{L}{2}$. This guarantees that in a region $(\tau_{1},\tau_{1}+\frac{L}{2})$ following a Dirac at $\tau_{1}$, there are at least 2 output samples, namely $y(t_{1})$ and $y(t_{2})$, as depicted in Fig. 8.

Then, in the interval $I=(t_{2}-\frac{L}{2},t_{1})$, which does not contain knots of either $\varphi(t-t_{1})$ or $\varphi(t-t_{2})$, we can reproduce two exponentials as described in Section II-B2. Specifically, we can find coefficients $c_{m,n}^{I}$ such that:

[TABLE]

We then define the signal moments $s_{m}$ as follows:

[TABLE]

where $b_{1}\vcentcolon=x_{1}e^{j\omega_{0}\tau_{1}}$, $u_{1}\vcentcolon=e^{j\lambda\tau_{1}}$, the frequencies $\omega_{m}=\omega_{0}+\lambda m$, for $m\in\{0,1\}$, and $\lambda=-2\omega_{0}$.

In these derivations, $(a)$ follows from Eq. (1), $(b)$ from the linearity of the inner product, and $(c)$ from Eq. (16). Moreover, $(d)$ holds since $t_{1},t_{2}\in(\tau_{1},\tau_{1}+\frac{L}{2})$, which means $\tau_{1}\in I$, and from the local exponential reproduction property of $\varphi(t)$ in the region $I$, as given in Eq. (17).

The unknowns $\{b_{1},u_{1}\}$ can be uniquely retrieved from the signal moments, as follows: $b_{1}=s_{0}$ and $u_{1}=\frac{s_{1}}{s_{0}}$. More generally, the parameters $\{b_{1},u_{1}\}$ can also be found using the annihilating filter method [46], also known as Prony’s method [38] (see Appendix A). Then, we get the Dirac’s amplitude and location, using $b_{1}=x_{1}e^{j\omega_{0}\tau_{1}}$ and $u_{1}=e^{j\lambda\tau_{1}}$. ∎

III-B Estimation of a Stream of Diracs

Let us now consider the case of a stream of Diracs:

[TABLE]

Proposition 2**.**

Let the sampling kernel $\varphi(t)$ be a second-order E-spline of support of length $L$, defined as in Eq. (12), with $\omega_{1}=-\omega_{0}$ and $0<\omega_{0}\leq\frac{\pi}{L}$. The filter $\varphi(t)$ and its non-uniform shifts can reproduce the exponentials $e^{j\omega_{0}t}$ and $e^{j\omega_{1}t}$ as in Eq. (13). In addition, suppose that the reference signal $g(t)=A\cos(w_{s}t)$ of the comparator has amplitude $A>|x_{k}|$, $\forall k$ and period $T_{s}<\frac{2L}{5}$, and that the minimum spacing between consecutive Diracs is larger than $L$. Then, the timing information $\{t_{1},t_{2},...,t_{N}\}$ provided by the device shown in Fig. 1 is a sufficient representation of a stream of Diracs as in Eq. (19).

Proof.

The input stream of Diracs can be sequentially estimated as follows. The first Dirac $x_{1}\delta(t-\tau_{1})$ can be uniquely estimated using the first two non-zero samples $y(t_{1})$ and $y(t_{2})$, as presented in Section III-A. Once we know $\tau_{1}$, we retrieve the first two non-zero samples $y(t_{n})$ and $y(t_{n+1})$ located after $\tau_{1}+L$, and use these to estimate the second Dirac $x_{2}\delta(t-\tau_{2})$. We then sequentially retrieve the next Dirac using the first two non-zero samples located after $\tau_{2}+L$, as illustrated in Fig. 9. In what follows we show that once $x_{1}\delta(t-\tau_{1})$ has been estimated, we can use $y(t_{n})$ and $y(t_{n+1})$ to estimate the second Dirac in the stream. Since we assume that the separation between input Diracs is larger than the length $L$ of the kernel’s support, then the location $\tau_{2}$ of the second Dirac satisfies $\tau_{1}+L<\tau_{2}<t_{n}$. Moreover, provided the period of the comparator’s signal satisfies $T_{s}<\frac{2L}{5}$, Bolzano’s intermediate value theorem [45] guarantees that $y(t_{n}),y(t_{n+1})\in(\tau_{2},\tau_{2}+\frac{L}{2})$, as previously outlined in Section III-A. Then, the interval $I=(t_{n+1}-\frac{L}{2},t_{n})$ contains no knots of either $\varphi(t-t_{n})$ or $\varphi(t-t_{n+1})$, and perfect exponential reproduction can be achieved. Hence we can compute the signal moments using similar derivations as in Eq. (18):

[TABLE]

Finally, we can estimate $x_{2}$ and $\tau_{2}$ from $s_{m}$, using Prony’s method. Once the second Dirac has been estimated, we use subsequent non-uniform output samples after $\tau_{2}+L$ in order to sequentially retrieve the next Diracs. ∎

The time encoding of the stream of Diracs is depicted in Fig. 9. Here, the filter is a second-order E-spline, of support length $L=2$, which can reproduce the exponentials $e^{\pm j\frac{\pi}{3}t}$. The frequency of the comparator’s test signal is $f_{s}=1.26>\frac{5}{2L}$ and the separation between Diracs is at least $L=2$. The amplitudes and locations of the estimated Diracs are exact to numerical precision.

III-C Multi-channel Estimation of Bursts of Diracs

Let us now consider a sequence of bursts of $K$ Diracs:

[TABLE]

where the amplitudes $x_{b,k}$ in the same burst $b$ have the same sign and satisfy $|x_{b,k}|<A_{max}$.

Proposition 3**.**

Let us consider a system of $M\geq K$ TEM devices as in Fig. 1.The filter $\varphi(t)$ of the $m^{th}$ TEM is a second-order E-spline whose support has length $L$, and which can reproduce two different exponentials, $e^{j\omega_{m_{0}}t}$ and $e^{j\omega_{m_{1}}t}$, with $\omega_{m_{0}}=\omega_{0}+\lambda m$, $\lambda=\frac{-2\omega_{0}}{2M-1}$, $0<\omega_{0}\leq\frac{\pi}{L}$, $\omega_{m_{1}}=-\omega_{m_{0}}$, and $m=0,1,...,M-1$. Furthermore, suppose the reference signal $g(t)=A\cos(w_{s}t)$ has amplitude $A>KA_{max}$ and period $T_{s}<\frac{2L}{7}$. In addition, let us assume the spacing between consecutive bursts is larger than $L$, and the maximum separation between the last and first Dirac in any burst $b$ satisfies $\tau_{b,K}-\tau_{b,1}<\frac{T_{s}}{2}$. Then, the timing information $t_{1,m},t_{2,m},...,t_{N,m}$ for $m=0,1,...,M-1$ provided by $M$ devices as in Fig. 1 is a sufficient representation of bursts of $K$ Diracs as in Eq. (20).

Proof.

See Appendix B. ∎

We summarize the results in this section by showing possible choices of the hyperparameters of the crossing TEM, and how they influence the density of output samples. This relationship is presented in Table I, for the case of a sequence of bursts of $K$ Diracs. Here, $M$ is the minimum number of channels, $P$ is the order of the sampling kernel for each channel, $L$ is the length of the support of the sampling kernel and $f_{s}^{min}$ the minimum frequency of the comparator’s sinusoidal reference. The table shows both the average sample density, as well as the ideal sample density111The ideal sample density is computed under the assumption that two samples are necessary to reconstruct one Dirac and that on average there are $K$ Diracs in an interval $L+\epsilon$ with $\epsilon>0$. required for perfect estimation of each of the bursts of $K$ Diracs.

We conclude this section by making the observation that the number of redundant samples of the crossing TEM is large, since samples are recorded even when the input is zero. In this case, output samples are recorded at the time instants when the sinusoidal reference signal crosses zero. In what follows, we aim to use the same decoding framework, however with a more efficient acquisition device, the integrate-and-fire TEM.

IV Perfect Recovery from Timing Information obtained with an Integrate-and-fire TEM

We now shift our focus on the integrate-and-fire TEM in Fig. 3. In particular, we show how to perfectly estimate an input Dirac, and extend this method to streams and bursts of Diracs, streams of pulses as well as piecewise constant signals. The retrieval of these signals from their timing information is perfect, provided the threshold of the trigger comparator is small enough to ensure a sufficient density of output samples. As it will become evident in Section VI, an important feature of the integrate-and-fire model is that it can be more efficient than the comparator or a system based on uniform sampling, in the case of input signals with a small number of Diracs, because it leads to a smaller number of samples.

IV-A Estimation of an Input Dirac

Proposition 4**.**

Let the sampling kernel $\varphi(t)$ be a second-order E-spline of support of length $L$, defined as in Eq. (12), with $\omega_{1}=-\omega_{0}$ and $0<\omega_{0}\leq\frac{\pi}{L}$. The filter $\varphi(t)$ and its non-uniform shifts can reproduce the exponentials $e^{j\omega_{0}t}$ and $e^{j\omega_{1}t}$ as in Eq. (13). In addition, suppose that the trigger mark of the comparator satisfies:

[TABLE]

where $A_{min}$ is the absolute minimum amplitude of the Dirac. Then, the timing information $\{t_{1},t_{2},...,t_{N}\}$ provided by the integrate-and-fire TEM in Fig. 3 is a sufficient representation of an input Dirac as in Eq. (16).

Proof.

We will prove that the upper bound on $C_{T}$ guarantees that the integrated filtered input $y(t)=x_{1}\varphi(\tau_{1}-t)$ reaches the trigger mark at least three times in the interval $(\tau_{1},\tau_{1}+\frac{L}{2})$. We will then show how we can use the second and third output samples $y(t_{2})$ and $y(t_{3})$ to perfectly estimate the input Dirac, given that the integrated filtered input has no discontinuities in the interval $(\tau_{1},\tau_{1}+\frac{L}{2})$.

First, we note that:

[TABLE]

where $(a)$ follows from Eq. (12), given $\alpha_{0}=-j\omega_{0}$, $\alpha_{1}=-j\omega_{1}$ and $\omega_{1}=-\omega_{0}$.

Then, we assume for simplicity that the Dirac’s amplitude satisfies $x_{1}>0$ and re-write the upper bound on $C_{T}$ as:

[TABLE]

where $(b)$ follows from Eq. (21).

Furthermore, from Eq. (2) and (4), we know that:

[TABLE]

where $(c)$ follows from Eq. (3) and given the input signal is $x(t)=x_{1}\delta(t-\tau_{1})$.

Then, from Eq. (22) and Eq. (23), we obtain the inequality:

[TABLE]

Using the hypothesis $\omega_{1}=-\omega_{0}$, together with Eq. (12), we obtain $\varphi(\tau_{1}-t)=\frac{\sin(\omega_{0}(t-\tau_{1}))}{\omega_{0}}$, for $t\in(\tau_{1},\tau_{1}+\frac{L}{2})$. Given the assumption $0<\omega_{0}\leq\frac{\pi}{L}$, we get $0\leq\omega_{0}(t-\tau_{1})\leq\frac{\pi}{2}$ and therefore, $\sin(\omega_{0}(t-\tau_{1}))>0$ for $t\in(\tau_{1},\tau_{1}+\frac{L}{2})$. Hence, since $\varphi(\tau_{1}-t)$ is positive in the range $(\tau_{1},\tau_{1}+\frac{L}{2})$ and using Eq. (24), we get that $t_{3}<\tau_{1}+\frac{L}{2}$.

As a result, the locations of the first non-uniform output samples satisfy $t_{1},t_{2},t_{3}\in(\tau_{1},\tau_{1}+\frac{L}{2})$, and can be computed using Eq. (8) and Eq. (7) as follows:

[TABLE]

[TABLE]

[TABLE]

for $\theta_{1}=t_{1}-\tau_{1}$, $\theta_{2}=t_{2}-t_{1}$ and $\theta_{3}=t_{3}-t_{2}$.

Furthermore, since $\varphi(t)$ is a second-order E-spline which can reproduce the exponentials $e^{j\omega_{0}t}$ and $e^{j\omega_{1}t}$ as in Eq. (12), and given the definition of $q_{\theta_{n}}(t)$ in Eq. (6), we have that:

[TABLE]

for $t\in(t_{1}-\frac{L}{2},t_{1})$.

Similarly:

[TABLE]

for $t\in(t_{2}-\frac{L}{2},t_{2})$, and

[TABLE]

for $t\in(t_{3}-\frac{L}{2},t_{3})$.

The shifted kernel $(\varphi*q_{\theta_{1}})(t-\tau_{1})$ depends on the Dirac’s location $\tau_{1}$, and hence its shape cannot be determined a-priori. On the other hand, the shifted kernels $(\varphi*q_{\theta_{2}})(t-t_{1})$ and $(\varphi*q_{\theta_{3}})(t-t_{2})$ are independent of $\tau_{1}$ and can be written as a linear combination of the exponentials $e^{j\omega_{0}t}$ and $e^{j\omega_{1}t}$, for $t\in(t_{3}-\frac{L}{2},t_{1})$. Therefore, in the interval $I=(t_{3}-\frac{L}{2},t_{1})$, where there are no knots of either the shifted kernel $(\varphi*q_{\theta_{2}})(t-t_{1})$ or $(\varphi*q_{\theta_{3}})(t-t_{2})$, we can use the proof in Section II-B2 to find the unique coefficients $c_{m,2}^{I}$ and $c_{m,3}^{I}$ such that:

[TABLE]

for $m\in\{0,1\}$ and $t\in(t_{3}-\frac{L}{2},t_{1})$.

Then, we can define the signal moments as:

[TABLE]

In the derivations above, $(d)$ follows from Eq. (16), (25) and (26), and $(e)$ follows from $\tau_{1}\in(t_{3}-\frac{L}{2},t_{1})$ which is true given Eq. (24), and since the property in Eq. (27) holds within $(t_{3}-\frac{L}{2},t_{1})$. Finally, using Prony’s method we can uniquely estimate parameters $x_{1}$ and $\tau_{1}$, from the two signal moments $s_{m}$ given by Eq. (28), for $m\in\{0,1\}$ and $\omega_{1}=-\omega_{0}$.

∎

IV-B Estimation of a Stream of Diracs

Proposition 5**.**

Let the sampling kernel $\varphi(t)$ be a second-order E-spline of support of length $L$, defined as in Eq. (12), with $\omega_{1}=-\omega_{0}$ and $0<\omega_{0}\leq\frac{\pi}{L}$. The filter $\varphi(t)$ and its non-uniform shifts can reproduce the exponentials $e^{j\omega_{0}t}$ and $e^{j\omega_{1}t}$ as in Eq. (13). In addition, assume that the minimum separation between consecutive Diracs is $L$ and the trigger mark of the comparator satisfies:

[TABLE]

where $A_{\min}$ is the absolute minimum amplitude of any Dirac in the input signal.

Then, the timing information $\{t_{1},t_{2},...,t_{N}\}$ provided by the integrate-and-fire TEM in Fig. 3 is a sufficient representation of a stream of Diracs as in Eq. (19).

Proof.

The first Dirac $\delta_{1}=x_{1}\delta(t-\tau_{1})$ can be correctly estimated using the method in Section IV-A, since Eq. (29) satisfies the requirements of Proposition 4. Then, suppose we aim to estimate the second Dirac in the input signal, and let us assume for simplicity that its amplitude satisfies $x_{2}>0$. Moreover, let us denote the output spike locations in the interval $(\tau_{1},\tau_{1}+L)$ with $t_{1},t_{2},...,t_{n-1}$, and the time information after $\tau_{1}+L$ with $t_{n},t_{n+1},...,t_{N}$. Then, given the hypothesis that the minimum separation between consecutive Diracs is $L$, the location of the second Dirac must satisfy $\tau_{2}\in(\tau_{1}+L,t_{n})$. We also have that:

[TABLE]

where $(a)$ follows from Eq. (12), for $\omega_{1}=-\omega_{0}$.

This shows the upper bound in Eq. (29) is equivalent to:

[TABLE]

Furthermore, we have that:

[TABLE]

where $(b)$ follows from Eq. (2), and $(c)$ holds since $t_{n-1}$ and $t_{n}$ are consecutive output spikes, and $t_{n}>\tau_{2}>t_{n-1}$.

As a result, Eq. (30) and (31) give the following inequality:

[TABLE]

As shown in Section IV-A, the sampling kernel satisfies $\varphi(t)>0$ for $x_{2}>0$, within the interval $(\tau_{2},\tau_{2}+\frac{L}{2})$. This means that the inequality in Eq. (32) is equivalent to $t_{n+2}<\tau_{2}+\frac{L}{2}$, which guarantees that the output samples $y_{n}$, $y_{n+1}$ and $y_{n+2}$ occur in the time interval $(\tau_{2},\tau_{2}+\frac{L}{2})$. Using the model of Fig. 3, we compute these non-uniform output samples as:

[TABLE]

[TABLE]

[TABLE]

The sample $y(t_{n})$ contains information of both $\delta_{1}$ and $\delta_{2}$, and hence cannot be used for estimation of the latter Dirac. On the other hand, since $t_{n+1},t_{n+2}\in(\tau_{2},\tau_{2}+\frac{L}{2})$, we can use the samples $y_{n+1}$ sand $y_{n+2}$ to compute the signal moments as in Section IV-A:

[TABLE]

Once $\delta_{2}$ is estimated from $s_{m}$ using Prony’s method, we use the non-uniform output samples after $\tau_{2}+L$, in order to sequentially retrieve the next Diracs in the input signal. ∎

The sampling and reconstruction of a stream of $K=3$ Diracs of minimum absolute amplitude $A_{\min}=1$ are depicted in Fig. 10. Here, the filter is a second-order E-spline, of support of length $L=2$, which can reproduce the exponentials $e^{\pm j\frac{\pi}{3}t}$, and the comparator’s trigger mark is $C_{T}=0.11$, which satisfies Eq. (29). Fig. 10(b) shows the filtered input and the output of the integrator. The amplitudes and locations of the estimated Diracs are exact to numerical precision. Finally, in Fig. 10(c) we observe that there are no output spikes in a region where the input signal is constant (zero), which leads to small average density of samples.

IV-C Estimation of a Stream of Pulses

Let us now consider a stream of pulses of the form $x(t)*g(t)$, where $x(t)$ is defined in Eq. (19) and the support of $g(t)$ is $[-\epsilon,\epsilon]$. Filtering this signal with the second-order E-spline $\varphi(t)$ is equivalent to filtering the stream of Diracs $x(t)$ with the kernel $(\varphi*g)(t)$. As a case in point, let us consider the cosine-squared pulse $g(t)=\cos^{2}(t)$, and assume that $2\epsilon<\frac{L}{2}$, where $L$ is the length of the filter’s support. In addition, suppose we want to estimate the first pulse $(x_{1}*g)(t)$ in the stream $x(t)*g(t)$ and denote its timing information with $t_{1},t_{2},...,t_{N}$. The first three output samples can be computed as follows:

[TABLE]

[TABLE]

[TABLE]

where $x_{1}(t)=x_{1}\delta(t-\tau_{1})$ is the first Dirac in the stream $x(t)$ with $x_{1}>0$, $\theta_{1}=t_{1}-\tau_{1}$, $\theta_{2}=t_{2}-t_{1}$ and $\theta_{3}=t_{3}-t_{2}$.

Assuming $2\epsilon<\frac{L}{2}$ we can leverage the results in Eq. (14) to show that in the interval $(\epsilon-\frac{L}{2},-\epsilon)$, we get $(\varphi*g)(t)=Ae^{\alpha_{0}t}+Be^{\alpha_{1}t}$, for some constants $A$ and $B$. Then, in the interval $(t_{2}-\frac{L}{2}+\epsilon,t_{1}-\epsilon)$, the function $(\varphi*g*q_{\theta_{2}})(t-t_{1})$ can also be expressed as a linear combination of the exponentials $e^{\alpha_{0}t}$ and $e^{\alpha_{1}t}$. Similarly, $(\varphi*g*q_{\theta_{3}})(t-t_{2})$ is a linear combination of the same exponentials in the interval $(t_{3}-\frac{L}{2}+\epsilon,t_{2}-\epsilon)$. As a result, in the knot-free interval $I=(t_{3}-\frac{L}{2}+\epsilon,t_{1}-\epsilon)$, we can perfectly reproduce two exponentials as in Eq. (15), using the shifted kernels $(\varphi*g*q_{\theta_{n+1}})(t-t_{n})$, for $n=1,2$. We can then compute two signal moments as in Eq. (28), and retrieve the amplitude and location of the first Dirac $x_{1}\delta(t-\tau_{1})$ in the stream $x(t)$ using Prony’s method.

In order for these derivations to hold we need to ensure that $\tau_{1}\in I$, or in other words that $t_{1}>\tau_{1}+\epsilon$ and $t_{3}<\tau_{1}+\frac{L}{2}-\epsilon$. Since the filtered input corresponding to the first pulse satisfies $x_{1}(\varphi*g)(\tau_{1}-t)>0$, for $t\in(\tau_{1}-\epsilon,\tau_{1}+L+\epsilon)$ and $x_{1}(\varphi*g)(-t+\tau_{1})=0$ otherwise, the condition $t_{1}>\tau_{1}+\epsilon$ holds provided the trigger mark of the comparator satisfies:

[TABLE]

Using the same reasoning as in Section IV-B, the condition $t_{3}<\tau_{1}+\frac{L}{2}-\epsilon$ holds provided:

[TABLE]

where $A_{min}$ is the minimum amplitude of the Diracs in $x(t)$.

We also note that in order for Eq. (35) and (36) to be simultaneously satisfied, we need to impose additional constraints on $\epsilon$, such that:

[TABLE]

Finally, once the first pulse centered around $\tau_{1}$ has been estimated, and assuming a minimum separation between consecutive pulses of at least $L+2\epsilon$ (which is the length of the support of $(\varphi*g)(t)$), we can use subsequent samples after $\tau_{1}+L+2\epsilon$ to retrieve the next pulse $(x_{2}*g)(t)$ in the input.

The sampling and perfect retrieval of a stream of cosine-squared pulses are depicted in Fig. 11, for $C_{T}=0.8$.

IV-D Multi-channel Estimation of Bursts of Diracs

Let us now consider the estimation of a sequence of bursts of Diracs as in Eq. (20). This problem is equivalent to the estimation of a stream of Diracs, however, this time, the $K$ Diracs can be arbitrarily close to each other. Therefore, the estimation of a burst of $K$ Diracs involves retrieving a larger number of moments (at least $2K$) to accurately retrieve the Diracs. We employ a multi-channel scheme of $M\geq K$ different acquisition devices, each of which will help us compute 2 different signal moments. We will show that it is sufficient to record 2 output samples for each channel in order to perfectly reconstruct each burst in the input signal, and that the trigger mark of the threshold detector can be adjusted to ensure these output samples are located suitably close to the input burst of Diracs. Specifically, we need to ensure that the 2 samples we use for estimation have contributions from all the $K$ Diracs, and hence, occur after the last Dirac in the burst.

Proposition 6**.**

Let us consider a system of $M\geq K$ TEM devices as in Fig. 3. The filter $\varphi(t)$ of the $m^{th}$ TEM is a second-order E-spline whose support has length $L$, and which can reproduce two different exponentials, $e^{j\omega_{m_{0}}t}$ and $e^{j\omega_{m_{1}}t}$, with $\omega_{m_{0}}=\omega_{0}+\lambda m$, $\lambda=\frac{-2\omega_{0}}{2M-1}$, $0<\omega_{0}\leq\frac{\pi}{L}$, $\omega_{m_{1}}=-\omega_{m_{0}}$, and $m=0,1,...,M-1$. Moreover, let us assume the spacing between consecutive bursts is larger than $L$, and the maximum separation between the last and first Dirac in any burst $b$ satisfies $\tau_{b,K}-\tau_{b,1}<\frac{L}{2}$. In addition, suppose that the comparator’s trigger mark $C_{T}$ satisfies the following conditions for each device $m$ and burst $b$:

[TABLE]

[TABLE]

where $A_{\max}$ and $A_{\min}$ are the absolute maximum and minimum amplitudes of the input, respectively.

Then, the timing information $t_{1,m},t_{2,m},...,t_{N,m}$ for $m=0,1,...,M-1$ provided by $M$ devices as in Fig. 3 is a sufficient representation of bursts of $K$ Diracs as in Eq. (20).

Proof.

See Appendix C. ∎

Even though we considered the sampling of bursts of Diracs using a multi-channel system, it is possible under slightly more restrictive conditions, to achieve the same using a single TEM device. Therefore, for the sake of completeness, we state the following result without proof:

Proposition 7**.**

Let us consider the integrate-and-fire TEM in Fig. 3. Let the sampling kernel $\varphi_{P}(t)$ be an E-spline of order $P\geq 2K$ and support of length $L$, which can reproduce $P$ different exponentials $e^{j\omega_{m}t}$, with $\omega_{m}=\omega_{0}+m\lambda$, $m=0,1,...,P-1$, and $0<\omega_{0}\leq\frac{\pi}{L}$. In addition, setting $P$ even and $\lambda=\frac{-\pi}{P}$ ensures $\varphi(t)$ is a real-valued function. In this setting, let us assume the spacing between bursts is larger than $L$, and the separation between the last and first Diracs within any burst $b$ satisfies $\tau_{b,K}-\tau_{b,1}<\frac{L}{P}$. In addition, suppose the trigger mark of the comparator $C_{T}$ satisfies:

[TABLE]

[TABLE]

where $\Delta_{b}=\max(\tau_{b,K}-\tau_{b,1})$.

Then, the timing information $\{t_{1},t_{2},...,t_{N}\}$ provided by the integrate-and-fire TEM in Fig. 3 is a sufficient representation of a sequence of bursts of $K$ Diracs as in Eq. (20).

IV-E Estimation of Piecewise Constant Signals

Let us now consider a piecewise constant signal $x(t)$, and assume that we filter this with the derivative of an E-spline $\varphi(t)$ of order $P\geq 2$, obtained using Eq. (11). Filtering $x(t)$ with $\frac{d\varphi(t)}{dt}$ ensures that in a region where the input is constant, there are no output spikes, since $\frac{d\varphi(t)}{dt}$ has average value equal to zero. This leads to energy-efficient sampling of the piecewise constant signal, resulting in a small average number of output spikes. In this setting, the filtered input is given by:

[TABLE]

This shows that filtering a piecewise constant signal $x(t)$ with $\frac{d\varphi(t)}{dt}$ is equivalent to filtering the stream of Diracs corresponding to the discontinuities of the piecewise constant signal with the E-spline $\varphi(t)$. The discontinuities of $\frac{dx(t)}{dt}$ can be estimated from the output spikes, by extending the results of Proposition 5 to the case of a $P$-order E-spline $\varphi_{P}(t)$, with $P\geq 2$. In this case, the E-spline $\varphi_{P}(t)$ of support of length $L$ can reproduce $P\geq 2$ different complex exponentials $e^{j\omega_{m}t}$, with $\omega_{m}=\omega_{0}+\lambda m$. and $m=0,1,...,P-1$. Moreover, choosing $\lambda=\frac{-2\omega_{0}}{P-1}$ and $P$ even ensures the kernel $\varphi_{P}(t)$ is a real-valued function. As before, the separation between consecutive Diracs must be larger than $L$ and the trigger mark of the comparator satisfies:

[TABLE]

Suppose we want to estimate the $k^{th}$ discontinuity in $\frac{dx(t)}{dt}$, of amplitude $z_{k}$ and located at $\tau_{k}$, and let us denote the locations of the first output spikes after $\tau_{k}$ with $t_{n},t_{n+1},...t_{N}$. Then, using a similar proof as in Section IV-B, we can show that the constraint in Eq. (41) guarantees that $\tau_{k}\in I=(t_{n+P}-\frac{L}{P},t_{n})$. Then, we can compute the following signal moments:

[TABLE]

In these derivations, $(a)$ follows from Eq. (7), and $(b)$ holds given $\tau_{k}\in(t_{n+P}-\frac{L}{P},t_{n})$, and the fact that none of the kernels $(\varphi_{P}*q_{\theta_{n+i}})(\tau_{k}-t_{n+i-1})$ have any discontinuities in $(t_{n+P}-\frac{L}{P},t_{n})$, for $i=1,2,...,P$. As before, we can use Prony’s method to estimate $z_{k}$ and $\tau_{k}$ from the signal moments $s_{m}$. Finally, we can retrieve the piecewise constant signal $x(t)$ once we have estimated its discontinuities $\frac{dx(t)}{dt}$.

The sampling and reconstruction of a piecewise constant signal are depicted in Fig. 12. The filter is the derivative of the fourth-order E-spline, with support length $L=4$, as seen in Fig. 12(b), the separation between input discontinuities is larger than the length of the kernel’s support as depicted in Fig. 12(a), and the comparator’s trigger mark is $C_{T}=0.001$. The estimation of the input is exact to numerical precision.

Similarly, the results of Propositions 6 and 7 can be extended to the case of a piecewise constant signal $x(t)$, where the discontinuities $\frac{dx(t)}{dt}$ are bursts of arbitrarily close Diracs, as in Eq. (20). For example, in Fig. 13, we show the time encoding and perfect decoding of a piecewise constant signal, with two arbitrarily close discontinuities.

We conclude this section by summarizing possible choices of hyperparameters in our sampling framework based on an integrate-and-fire system. Specifically, let us consider the case of streams of bursts of $K$ Diracs, and discuss the relationship between the sampling kernel and the trigger mark of the comparator, and how these parameters determine the density of output samples. The sampling kernel is assumed to be the E-spline given in Eq. (11) of order $P$, and support length $L=P$. Furthermore, the conditions of the trigger mark $C_{T}$ ensure that the output samples used for reconstruction are located sufficiently close to the burst of Diracs, and in a region where the filtered input is continuous. As a result, these conditions depend on the separation $\Delta_{b}$ between the Diracs, as well as on the location of the knots of the sampling kernel, which in turn depends on the length of the support of this kernel. Setting $C_{T}$ to its maximum theoretical value ensures that the number of samples is minimised.

The choice of the hyperparameters of the integrate-and-fire TEM is summarised in Table II. Here, a burst of $K=2$ Diracs was time encoded using an $M$-channel system. The amplitudes of each of the Diracs was chosen uniformly at random in the interval $[1,2]$ and the trigger mark $C_{T}$ computed using Eq. (38). The results were averaged over 100 different experiments.

Finally, we make the remark that only some of the output samples are used for input reconstruction. For online reconstruction applications, these are the only samples that need to be stored. This is depicted in Fig. 14, where we highlight in red the samples used for reconstruction of each burst of Diracs, of one of the two channels. Only the second and third output samples, located at $t_{2}^{(1)}$ and $t_{3}^{(1)}$ need to be recorded and used to retrieve the first burst of 2 Diracs. Once the first burst has been estimated, we record the second and third output samples after $\tau_{2}+L$, located at $t_{2}^{(2)}$ and $t_{3}^{(2)}$ in order to retrieve the next burst.

V Generalized Time-based Sampling

To highlight the potential practical implications of the methods developed in the previous sections, we present here extensions of our framework to deal with arbitrary kernels and the noisy scenario, and show that reliable input reconstruction can be achieved also in these scenarios.

V-A Sampling with Arbitrary Kernels

In the previous sections we have presented methods for perfect retrieval of certain classes of non-bandlimited signals from timing information. We have seen that these methods require the sampling kernel $\varphi(t)$ to locally reproduce exponentials, in order to be able to map this problem to Prony’s method. In reality, however, the sampling kernel may not have the exponential reproducing property as in Eq. (13). Let us now consider an arbitrary kernel $\tilde{\varphi}(t)$, and find a linear combination of its non-uniform shifted versions that gives the best approximation of $P$ exponentials $f(t)=e^{j\omega_{m}t}$ within an interval $I$, for $\omega_{m}=\omega_{0}+\lambda m$, $m=0,1,...,P-1$, and $\lambda=\frac{-2\omega_{0}}{P-1}$. In other words, we want to find the optimal coefficients $c_{m,n}^{I}$ such that:

[TABLE]

for $t\in I$ and $n=1,2,...,N$, with $N$ being the number of kernels $\tilde{\varphi}(t-t_{n})$ overlapping $I$.

We find the coefficients $c_{m,n}$ using the least-squares approximation method described in [47]. The coefficients are computed using the orthogonal projection of $f(t)$ onto the space spanned by the non-uniform shifts $\tilde{\varphi}(t-t_{n})$, such that:

[TABLE]

for $t\in I$ and $n=1,2,..N$.

Furthermore, Eq. (43) is equivalent to:

[TABLE]

which represents a system of $N$ equations from which we can determine the $N$ coefficients $c_{m,k}^{I}$, for each $m=0,1,...,P-1$.

We then use the calculated coefficients $c_{m,k}^{I}$ to compute the signal moments as in Section IV. Finally, the estimation of the input can be further refined using the Cadzow iterative algorithm in order to increase the accuracy of the signal moments, before applying Prony’s method [48, 49].

The sampling and reconstruction of bursts of 2 Diracs are depicted in Fig. 16. We use the multi-channel estimation method presented in Section IV-D, where the filter of each channel is a third order B-spline $\beta_{3}(t)$, such that the modified kernel $(\beta_{3}*q_{\theta_{n}})(t)$ in Eq. (5) cannot reproduce exponentials. Moreover, we aim to approximately reproduce 4 different exponentials for each channel, and hence we require a number of 4 non-uniform samples, as discussed in Section II-B. In Fig. 15, we depict the approximate exponential reproduction in Eq. (42), within the interval $I=(0.82,1.4)$ overlapping the first burst of Diracs. The mean-squared error of the exponential reproduction within this interval is $9.47\times 10^{-13}$. Finally, the estimation of the input is close to exact, as depicted in Fig. 16(c). In particular, the mean-squared error in the time locations of the Diracs is $2.44\times 10^{-4}$, and the mean-squared error in the amplitudes of the Diracs is $2.04\times 10^{-10}$.

V-B Robustness of the Integrate-and-fire TEM to Noise

In many practical circumstances, the input signal is corrupted by noise, which is typically assumed to be white, additive Gaussian noise. When this happens, the non-uniform times $\{t_{n}\}$ change which means that the sequence of moments $s_{m}$ is also corrupted, and perfect reconstruction may no longer be possible. Nevertheless, if the noise has average value equal to 0, it is in part removed by the integrator in the TEM, as a result of the averaging effect of the integral.

In Fig. 17 we show the reconstruction of a piecewise constant signal corrupted by white, additive Gaussian noise, using the method in Section IV-E. The filter is the derivative of a fourth-order E-spline with support length $L=4$ which can reproduce the exponentials $e^{\pm j\frac{\pi}{3}t}$ and $e^{\pm j\frac{\pi}{6}t}$, the trigger mark of the comparator is $C_{T}=0.001$, the standard deviation of the noise is $\sigma=0.1$ (SNR$=21.56$dB), and the separation between consecutive discontinuities of the input is larger than $L$. The reconstruction of the input from noisy samples is very accurate. A quantitative analysis of the effect of noise on the retrieval of this piecewise constant signal is presented in Table III. The table shows the error of the estimated locations and the relative error of the estimated amplitudes of the discontinuities in the input signal, averaged over 10000 experiments.

In Fig. 18 we show the reconstruction errors, for the case of a stream of Diracs, for different SNR values, averaged over 1000 experiments. Here, the input signal is corrupted by white, additive Gaussian noise, and the sampling kernel is a second-order E-spline whose support has length $L=2$, defined as in Eq. (12), for $\alpha_{0}=-\alpha_{1}=j\frac{\pi}{3}$. When SNR$=20dB$, the amplitude mean-squared error is $1.17\times 10^{-3}$ and the mean-squared error in time locations is $2.5\times 10^{-3}$. Finally, in Fig. 19 we depict the estimation of an input stream of Diracs corrupted by noise, from its timing information, when SNR$=10dB$.

V-C Time Encoding and Decoding of Bursts of Diracs of Arbitrary Signs

In Section IV-D we presented sufficient conditions for perfect retrieval of bursts of Diracs defined as in Eq. (20). These conditions rely on various assumptions, including that the amplitudes $x_{b,k}$ of the Diracs $x_{b,k}\delta(t-\tau_{b,k})$ in the same burst $b$, have the same sign. In reality, this assumption may not always hold, and in this section we empirically show that the reconstruction framework presented in this paper usually performs well also when the amplitudes of the Diracs in the same burst have opposite signs. We consider the estimation of a burst of 2 Diracs from its time encoded information using a 2-channel approach. We assume that the filter of each channel is a second-order E-spline defined as in Eq. (12) with support length $L=2$. We denote the output information of channel 1 with $t_{1,1},t_{2,1},...,t_{N,1}$, and that of channel 2 with $t_{1,2},t_{2,2},...,t_{N,2}$. We assume that the amplitudes of these Diracs are distributed as Gaussian variables, of mean $\mu=0$ and variance $\sigma=1$.

The decoding scheme presented in Section IV-D showed that we can reliably use the samples $y(t_{3,1}),y(t_{4,1})$ of the first channel and $y(t_{3,2}),y(t_{4,2})$ of the second channel, in order to perfectly retrieve an input burst of 2 Diracs of same sign. The sufficient conditions on the trigger mark of the integrator given in Eq. (37) and (38) ensure that these samples are located after the second Dirac at $\tau_{2}$. Here, we choose $C_{T}$ below its minimum theoretical value given in Eq. (37), in order to ensure a sufficient number of output samples is obtained, even when the two Diracs in the input have opposite signs and are located closely to each other. However, when lowering $C_{T}$, the samples $y(t_{3,1}),y(t_{4,1})$ and $y(t_{3,2}),y(t_{4,2})$ are not guaranteed to occur after $\tau_{2}$, and hence, may not be reliably used for reconstruction of both Diracs. Therefore, we adjust the reconstruction scheme as follows. Using $y(t_{2,1}),y(t_{3,1})$ of the first channel and $y(t_{2,2}),y(t_{3,2})$ of the second channel we compute the signal moments $s_{m}$ as described in Appendix C and then build matrix $S$ as in Appendix A. If the rank of the matrix $S$ is 1, then $t_{3,1}<\tau_{2}$ and $t_{3,2}<\tau_{2}$. Hence, we can use $y(t_{2,1}),y(t_{3,1})$ to estimate the first Dirac $x_{1}\delta(t-\tau_{1})$. Once the first Dirac has been estimated, we remove its contribution from the output spikes, and use the next non-zero samples in order to estimate the second Dirac $x_{2}\delta(t-\tau_{2})$. Otherwise, if the rank of matrix $S$ is 2, then at least for one of the channels $i$, we get $t_{i,3}>\tau_{2}$. As a result of the similarity between the sampling kernels of the two channels, it is likely that $t_{i,3}>\tau_{2}$ for both $i=1$ and $i=2$. In other words, the samples $y(t_{4,1}),y(t_{5,1})$ and $y(t_{4,2}),y(t_{5,2})$ are likely to have contributions from both Diracs and hence, we can use the method in Section IV-D to estimate the burst. In Table IV we show the probability of correct estimation of the 2 Diracs, against different values of $\Delta_{b}$ and trigger mark $C_{T}$, averaged over 1000 experiments. The results show that we still achieve perfect reconstruction in most cases. The reconstruction typically fails when the number of samples required for estimation is not achieved (for example, when the amplitudes of the Diracs are very small or when they have similar magnitudes, but opposite signs).

VI Density of Non-uniform Samples Obtained with an Integrate-and-fire TEM

In the previous sections, we have presented techniques for estimation of non-bandlimited signals from timing information. We have seen that perfect estimation can be achieved using simple algorithms, and physically realisable kernels. In this section we outline the fact that in many settings sampling based on timing using our integrate-and-fire system is an efficient way to acquire signals, resulting in a smaller density of samples, compared to classical sampling.

As a case in point we consider the retrieval of bursts of $K$ Diracs, described in Section IV-D. We have seen that perfect reconstruction from timing information can be achieved, provided the separation between consecutive bursts is at least $L$, and that the Diracs within any burst are sufficiently close. In particular, let us denote the maximum separation between the last and first Dirac within a burst with $\Delta=\max(\tau_{K}-\tau_{1})<\frac{L}{2}$, which can be determined according to Eq. (37) and (38). Moreover, let us assume the input is sufficiently sparse, such that the average separation between consecutive bursts is $L+S$, with $S>0$. Under these assumptions, the results in [6] show that in order to retrieve the $K$ Diracs from uniform samples, we need at least $2K$ samples within the interval $L-\Delta$ following the burst of Diracs. As a result, the uniform sampling period must satisfy $T\leq\frac{L-\Delta}{2K}$. Then, the number of uniform samples we record within an interval of length $L+S$ is $\frac{L+S}{T}=\frac{2K(L+S)}{L-\Delta}$. On the other hand, in the case of time encoding using the integrate-and-fire TEM in Fig. 3, the results in Section IV-D show that we need to record 4 output samples for each of the $K$ channels (or equivalently, $4K$ samples for the case of single-channel sampling), for each burst of $K$ Diracs. We note that Eq. (38) shows that in many situations, the TEM outputs more than 4 spikes per channel. Nevertheless, these samples can be discarded since they are not used in estimation. For example, one way to stop recording spikes once we have obtained 4 non-zero samples, is to increase the trigger mark $C_{T}$ of the comparator in Fig. 3, for a duration of $L-\Delta$.

Moreover, when the input is constant (zero), the integrate-and-fire TEM does not fire, and hence there are no output samples. Therefore, in an interval of size $L+S$, the number of stored samples from a $K$-Dirac burst is $4K$, $\forall S$.

Furthermore, $\frac{2K(L+S)}{L-\Delta}>4K$ for $S\geq L-2\Delta>0$ and $\forall K$, which shows that the average number of non-uniform spikes required for the retrieval of $K$ Diracs is lower than the number of uniform samples required to estimate the same number of free input parameters, when the input is sufficiently sparse.

VII Conclusions

In this work we established time encoding as an alternative sampling method for some classes of signals that are neither bandlimited, nor belong to shift-invariant subspaces. The proposed sampling scheme is based on first filtering the input signal, before retrieving the timing information using a crossing or integrate-and-fire TEM. We demonstrated sufficient conditions for the exact recovery of streams of Diracs, streams of pulses and piecewise constant signals, from their time-based samples. Central to our reconstruction methods is the use of specific filters that we proved can locally reproduce polynomials or exponentials. We further highlighted the potential of this new framework by showing that it is resilient to noise and that it can handle non-ideal filters. Finally, the diverse applications of previous results of finite rate of innovation theory [33, 34, 8, 32] also serve as evidence for the potential for real-world applications of the theoretical framework developed in this paper.

Appendix A Prony’s method

One way to solve the problem of estimating the parameters $\{b_{k},u_{k}\}_{k=1}^{K}$ from the sequence $s_{m}=\sum_{k=1}^{K}b_{k}u_{k}^{m}$ is given by the annihilating filter method, also referred to as Prony’s method [38]. The name of this approach comes from the observation that if we filter $s_{m}$ with a filter which has zeros at $\{u_{k}\}_{k=1}^{K}$, the output is zero, or in other words, this filter annihilates the sequence $s_{m}$.

The z-transform of the annihilating filter satisfies:

[TABLE]

which evaluates to zero when $z=u_{k}$.

Filtering the sequence $s_{m}$ with $h_{m}$ corresponds to the convolution of these sequences:

[TABLE]

where $(a)$ holds since $z=u_{k}$ gives $H(z)=0$ in Eq. (44).

Eq. (45) can be written in matricial form as follows:

[TABLE]

It can be shown that provided $\{b_{k}\}_{k=1}^{K}$ are non-zero and $\{u_{k}\}_{k=1}^{K}$ are distinct, matrix S has full row rank $K$, which means the solution h given by Eq. (46) is unique. Moreover, the solution h can be obtained by performing a singular value decomposition of S, where h is the singular vector corresponding to the zero singular value.

Then, once the coefficients $h_{m}$ of the polynomial $H(z)$ are known, the parameters $\{u_{k}\}_{k=1}^{K}$ are obtained from the roots of this filter. Finally, once $\{u_{k}\}_{k=1}^{K}$ are found, the parameters $\{b_{k}\}_{k=1}^{K}$ can be computed from the linear system of $K$ equations given by $s_{m}=\sum_{k=1}^{K}b_{k}u_{k}^{m}$, with $m=0,1,...,K-1$.

Appendix B

B-A Proof of Proposition 3

For simplicity, let us assume the number of devices equals the number of Diracs in a burst, i.e. $M=K$. Suppose we want to estimate the Diracs in the first burst, located at $\tau_{1,1},....,\tau_{1,K}$. Moreover, assume for simplicity that their amplitudes satisfy $x_{1,1},...,x_{1,K}>0$. In addition, let us consider the output of the $m^{th}$ TEM device, and denote its timing information with $\{t_{1},t_{2},...,t_{N}\}$.

Since we assume all the amplitudes in the first burst satisfy $0<x_{1,k}<A_{max}$, and since $0\leq\varphi(t)<1$, we get $0\leq y(t)$ and $y(t)=\sum_{k=1}^{K}x_{k}\varphi(\tau_{k}-t)<KA_{max}<A=\max(g(t))$.

Then, Bolzano’s intermediate value theorem [45] guarantees that the $m^{th}$ TEM outputs at most one sample in the interval $(\tau_{1,1},\tau_{1,K})$, given the assumption $\tau_{1,K}-\tau_{1,1}<\frac{T_{s}}{2}$, and the fact that $0\leq y(t)<\max(g(t))$. At the same time, this theorem also guarantees that the filtered input $y(t)$ crosses the sinusoidal reference signal in at least 3 points, within the window $(\tau_{1,1},\tau_{1,1}+\frac{7T_{s}}{4})$, such that $t_{3}-\tau_{1,1}\leq\frac{7T_{s}}{4}$. Moreover, the assumption $T_{s}\leq\frac{2L}{7}$ ensures that $t_{3}-\tau_{1,1}\leq\frac{L}{2}$. Hence, whilst the spike at $t_{1}$ may occur before $\tau_{1,K}$, the second and third spikes satisfy $t_{2},t_{3}\in(\tau_{1,K},\tau_{1,1}+\frac{L}{2})$, which means that $\tau_{1,1},\tau_{1,2},...\tau_{1,K}\in(t_{3}-\frac{L}{2},t_{2})$.

Since in the interval $I=(t_{3}-\frac{L}{2},t_{2})$ there are no knots of either $\varphi(t-t_{2})$ or $\varphi(t-t_{3})$, we can compute the following signal moments for the $m^{th}$ channel:

[TABLE]

where $i\in\{0,1\}$, and $\omega_{m_{0}}=\omega_{0}+\lambda m$ and $\omega_{m_{1}}=-\omega_{m_{0}}$.

In the derivations above, $(a)$ follows from Eq. (1), $(b)$ from the linearity of the inner product, and $(c)$ from the local exponential reproduction property of the sampling kernel described in Eq. (13), for $N=2$. Moreover, $(d)$ follows from Eq. (20), and given that $\tau_{1,1},\tau_{1,2},...\tau_{1,K}\in(t_{3}-\frac{L}{2},t_{1})$.

By using the same approach on each of the $K$ channels, we can retrieve $2K$ different moments and, due to the specific choice of exponents, the $2K$ moments can be expressed as:

[TABLE]

where $b_{k}=e^{j\omega_{0}\tau_{1,k}}$, $u_{k}=e^{j\lambda\tau_{1,k}}$, and $p=0,1,...,2M-1$.

We can then apply Prony’s method on $s_{p}$ to retrieve the $K$ amplitudes and the $K$ locations of the Diracs. Finally, we use subsequent output samples, located after $\tau_{1,K}+L$ to retrieve the free parameters of the Diracs in the second burst, and we reiterate the process for the following bursts.

The sampling and reconstruction of a sequence of bursts of 2 Diracs are depicted in Fig. 20. Here, the sampling kernel is a second-order E-spline for each channel, of support of length $L=2$, shown in Fig. 20(c) and (d). The first channel’s kernel reproduces the exponentials $e^{\pm j\frac{\pi}{3}t}$, whereas the second kernel reproduces $e^{\pm j\frac{\pi}{9}t}$. Moreover, the comparator’s reference signal has frequency $f_{s}=1.76>\frac{7}{2L}$, and the separation between consecutive bursts of Diracs is at least $L$. The amplitudes and locations of the estimated Diracs are exact.

Appendix C

C-A Proof of Proposition 6

The input stream of bursts of Diracs can be sequentially estimated as follows. We estimate the first burst using the first four non-zero samples of each channel and the methods presented below. We then retrieve the second burst using the first four non-zero samples of each channel located after $\tau_{1,K}+L$, where $\tau_{1,K}$ denotes the estimated location of the last Dirac in the first burst, and $L$ is the length of the kernel’s support. We then use the first non-zero samples located after $\tau_{2,K}+L$ to estimate the third burst, and repeat this procedure to estimate the subsequent bursts of Diracs.

Let us assume we want to retrieve burst $b$ and denote with $t_{n},t_{n+1},t_{n+2},t_{n+3}$ the first four output spikes located after $\tau_{b-1,K}+L$. Then we have that $t_{n}>\tau_{b,1}>t_{n-1}$, where $\tau_{b,1}$ is the location of the first Dirac in the $b^{th}$ burst. Furthermore, let us assume for simplicity that the Diracs in the $b^{th}$ burst satisfy $x_{b,1},...,x_{b,K}>0$, as depicted in Fig. 21.

In what follows, we show that the samples $y(t_{n+2})$ and $y(t_{n+3})$ can be reliably used to estimate the $b^{th}$ burst.

We first prove that the following conditions hold:

[TABLE]

and:

[TABLE]

We note that since we assume $x_{b,1},...,x_{b,K}>0$, the filtered input defined in Eq. (3) satisfies $f(\tau)>0$, and hence the condition in Eq. (47) is equivalent to:

[TABLE]

The left-hand side of this inequality can be expressed as:

[TABLE]

where $(a)$ holds given Eq. (2) and $(b)$ since $t_{n}>\tau_{b,1}>t_{n-1}$.

The right-hand side of Eq. (49) can be re-written as:

[TABLE]

which proves the inequality in Eq. (47).

In the derivations above, $(c)$ follows from the definition in Eq. (3) and since we assume $A_{max}>x_{b,1},...,x_{b,K}>0$. In addition, $(e)$ follows from Eq. (37) and $(f)$ from Eq. (50). Finally, condition $(d)$ follows from:

[TABLE]

where $(h)$ follows from the definition of $\varphi(\tau_{b,k}-\tau)$ in Eq. (12) for $\tau\in[\tau_{b,k},\tau_{b,K}]$ with $\tau_{b,K}<\tau_{b,k}+\frac{L}{2}$, and from the hypothesis that $\varphi(\tau)$ reproduces the exponentials $e^{\pm j\omega_{m_{0}}\tau}$. Moreover, $(i)$ follows from the hypothesis that $0<\omega_{m_{0}}\leq\frac{\pi}{L}$ which is equivalent to $0<\frac{\omega_{m_{0}}L}{2}\leq\frac{\pi}{2}$, and from the assumption that $\tau_{b,K}-\tau_{b,k}<\frac{L}{2}$, which means that $0<\omega_{m_{0}}(\tau_{b,K}-\tau_{b,k})<\frac{\pi}{2}$, and hence $1-\cos(\omega_{m_{0}}(\tau_{b,K}-\tau_{b,1}))>1-\cos(\omega_{m_{0}}(\tau_{b,K}-\tau_{b,k}))$ $\forall k=2,...,K$.

In Fig. 21 we notice that in some cases (third burst of 2 Diracs), the spike $t_{n}$ may occur in the interval $(\tau_{b,1},\tau_{b,K})$. Nevertheless, the condition in Eq. (49) ensures that the sample at $t_{n+1}$ happens after $\tau_{b,K}$.

Similarly, since $f(\tau)>0$, Eq. (48) is equivalent to:

[TABLE]

where the left-hand side can be expressed as:

[TABLE]

where $(a)$ follows from Eq. (3), $(b)$ follows from the definition of $\varphi(\tau_{b,k}-\tau)$ in Eq. (12) for $\tau\in(\tau_{b,k},\tau_{b,1}+\frac{L}{2})$, and $(c)$ follows from the hypothesis that $0<\omega_{m_{0}}\leq\frac{\pi}{L}$ which is equivalent to $0<\frac{\omega_{m_{0}}L}{2}\leq\frac{\pi}{2}$, and since $\tau_{b,k}-\tau_{b,1}<\frac{L}{2}$ $\forall k=2,...,K$. Moreover, $(d)$ holds since we assume $x_{b,1},...,x_{b,K}>0$, and $(e)$ follows from Eq. (38).

Finally, the right-hand side of Eq. (51) is equivalent to:

[TABLE]

hence proving the result in Eq. (48).

In these derivations, $(f)$ follows from Eq. (2), $(g)$ holds since $t_{n}>\tau_{b,1}>t_{n-1}$ and $(h)$ follows from Eq. (52).

The conditions in Eq. (47) and (48) ensure that the output samples $y(t_{n+2})$ and $y(t_{n+3})$ have contributions only from all the Diracs in the $b^{th}$ burst. These samples can be computed using Eq. (5) and (20) for each channel $m$, as follows:

[TABLE]

Similarly, we can write $y(t_{n+3})$ as:

[TABLE]

For each channel $m$, the signal $(\varphi*q_{\theta_{n+2}})(t-t_{n+1})$ is a linear combination of the exponentials $e^{j\omega_{m_{0}}t}$ and $e^{j\omega_{m_{1}}t}$, for $t\in(t_{n+2}-\frac{L}{2},t_{n+1})$, given Eq. (12) and Eq. (6). Similarly, $(\varphi*q_{\theta_{n+3}})(t-t_{n+2})$ is a linear combination of the exponentials $e^{j\omega_{m_{0}}t}$ and $e^{j\omega_{m_{1}}t}$, for $t\in(t_{n+3}-\frac{L}{2},t_{n+1})$. Therefore, in the interval $(t_{n+3}-\frac{L}{2},t_{n+1})$, where there are no knots of either $(\varphi*q_{\theta_{n+2}})(t-t_{n+1})$ or $(\varphi*q_{\theta_{n+3}})(t-t_{n+2})$, we use the proof in Section II-B2 to find unique $c_{m_{i},2}$ and $c_{m_{i},3}$ such that:

[TABLE]

for $i\in\{0,1\}$, $t\in[t_{n+3}-\frac{L}{2},t_{n+1}]$, $m_{0}=m$ and $m_{1}=2K-1-m$ (which ensures $\omega_{m_{1}}=-\omega_{m_{0}}$).

Then, for each channel $m$ we can compute the signal moments as before:

[TABLE]

where $i\in\{0,1\}$, $m_{0}=m$ and $m_{1}=2K-1-m$.

In the derivations above, (a) follows from Eq. (53) and (54), and (b) from $\tau_{b,1},...,\tau_{b,K}\in(t_{n+3}-\frac{L}{2},t_{n+1})$ and the fact that Eq. (55) holds within this interval. We can then uniquely retrieve the $2K$ input parameters of the $b^{th}$ burst from the $2K$ signal moments $s_{m_{i}}$ of all channels, using Prony’s method.

Finally, we make the observation that the inequalities in Eq. (37) and Eq. (38) impose additional constraints on the maximum separation between the Diracs in a burst $b$, namely on $\tau_{b,K}-\tau_{b,1}$. Specifically, we need to impose:

[TABLE]

which may give different constraints on the Dirac separation according to the filter characteristics, $A_{\max}$ and $A_{\min}$.