Signal Filtering Using Neuromorphic Measurements

1. Introduction

In contemporary computing, the majority of tasks necessitate some form of digital signal processing [1]. With the escalating computational capabilities of digital processing systems, there is a concurrent surge in power consumption. This is further exacerbated by the rapid advancements in the field of artificial intelligence. Neuromorphic sampling, or time encoding, is an alternative to traditional digital encoding that transforms an analog signal into a sequence of low-power time-based pulses, often referred to as a spike train. Neuromorphic sampling draws its inspiration from neuroscience and introduces a paradigm shift by significantly reducing power consumption during encoding and transmission [2,3]. Despite these advantages, as of now, there exists no equivalent of digital signal processing operations tailored to neuromorphic sampling. This unexplored territory holds the promise of groundbreaking developments in low-power and efficient signal processing.

In this paper, we address the problem of filtering a signal using its neuromorphic measurements, thus extending the principle of digital signal filtering for the case of neuromorphic sampling. In the literature, the concept of spike train filtering predominantly refers to convolving the filter function with a sequence of Diracs centered in the spike times [4]. This is not an operation on the analog signal that resulted in those neuromorphic samples. Moreover, in the context of neuromorphic measurements generated from multidimensional signals, such as those generated by event-cameras [5], filtering also refers to performing multidimensional spatial convolution [6,7]. The case of filtering the analog signal via its neuromorphic measurements was proposed in [8], but the process also involves signal reconstruction. A filtering approach of neuromorphic signals without reconstruction was first studied in [9].

Therefore, conventionally, to process the signal underlying a sequence of spike measurements, the signal is recovered in the analog domain, followed by filtering and neuromorphic sampling. There are a number of drawbacks with this approach. First, this method does not exploit the power consumption advantage of neuromorphic measurements. Second, this approach is heavier computationally due to the complexity of signal reconstruction from neuromorphic samples [10,11]. Third, reconstruction is only possible if the input satisfies some restrictive smoothness or sparsity constraints. However, for conventional sampling, the process of digital filtering is independent of the charateristics of the signal that generated the measurements.

In this paper, we derive a direct mapping between the time encoded inputs and outputs of an analog filter. The proposed mapping forms the basis for a practical filtering algorithm of the underlying signal corresponding to some given neuromorphic measurements, without direct access to the signal. We introduce theoretical guarantees and error bounds for the measurements generated with the proposed algorithm. Throrugh numerical simulations, we demonstrate the performance of our method in terms of speed, but also reduced restrictions on the input signal, in comparison with the existing conventional method.

This paper is structured as follows. Section 2 presents a brief review of the time encoding model used in this paper. Section 3 introduces the proposed problem. Section 4 describes the proposed filtering method. Numerical results are in Section 5. Section 6 presents the concluding remarks.

2. Time Encoding

The time encoding machine (TEM) is a conceptualization of neuromorphic sampling that maps an input $u\left(t\right)$ into a sequence of samples ${t_k}_{k\in }$. The particular TEM considered here is the integrate-and-fire (IF) TEM, which is inspired from neuroscience. Consequently, sequence ${t_k}_{k\in }$ is called a spike train, where the spike refers to the firing of an action potential, representing the information transmission method in the mammalian cortex. Previously, the IF model was used for system identification of biological neurons [12,13], to perform machine learning tasks [14,15] but also for input reconstruction [11,16,17,18,19]. The IF model adds an input $g\left(t\right)$ with bias parameter b, and subsequently integrates the output to generate a strictly increasing function $y\left(t\right)$. When $y\left(t\right)$ crosses threshold $\delta$, the integrator is reset and the IF generates an output spike time $t_k$. The IF TEM is described by the following equations

$${\int }_{t_k}^{t_{k+1}}u\left(s\right)ds=\delta -b{t}_{k+1}-t_k,k{\in }_{+}^{*}.$$

(1)

Without reducing the generality, it is assumed that $t_0=0$. A common assumption is that the input is bounded by $c{\in }_{+}$, such that $u\left(t\right)\le c<b$. This bound enables deriving the following density guarantees [2]

$$\frac{\delta }{b+c}\le t_{k+1}-t_k\le \frac{\delta }{b-c}.$$

(2)

Signal $u\left(t\right)$ is, in the general case, not recoverable from ${\left\{t_k\right\}}_{k\in \mathbb{Z}}$. To ensure it can be reconstructed, we require imposing restrictive assumptions. A common assumption is that $u\left(t\right)$ belongs to $P{W}_{\mathsf{\Omega }}$, the space of functions bandlimited to $\mathsf{\Omega }\mathrm{rad}/\mathrm{s}$ that are also square integrable, i.e., $u\in L^2$. If this assumption is satisfied, then $u\left(t\right)$ can be recovered from ${\left\{t_k\right\}}_{k\in \mathbb{Z}}$ if

$$\frac{\delta }{b-c}<\frac{\pi }{\mathsf{\Omega }}.$$

(3)

Further generalizations were introduced for the case where $u\left(t\right)$ is a function in a shift-invariant space [11,17] or in a space with finite rate of innovation [18,19]. However, if $u\left(t\right)$ does not belong to one of the classes above, or if it is bandlimited and does not satisfy (3), then the conventional theory does not allow any processing of the signal $u\left(t\right)$ via its samples $t_k$. The same is not true for conventional digital signals, that can be processed even when they’re not sampled at Nyquist rate. We will show how some types of processing such as filtering is still doable even when (3) is not true.

3. Problem Statement

Here we formulate the proposed signal filtering problem as follows. Assume the neuron input is continuous, i.e., $u\in C$. To satisfy the neuron encoding requirement in Section 2, we assume the input is bounded, such that $\left|u\left(t\right)\right|\le c<b,\forall t\in \mathbb{R}$. Furthermore, we assume that the input is absolutely integrable $u\in L^1\left(\mathbb{R}\right)$ and square integrable $u\in L^2\left(\mathbb{R}\right)$. Following the same idea as in digital filtering, we do not impose any general conditions on the bandwidth, smoothness, or sparsity of the analog signal in order to compute the filter output.

The filter is assumed linear, with impulse response $g\left(t\right)$ that is continuous, $g\in C$, and absolutely integrable, $g\in L^1\left(\mathbb{R}\right)$. The output of the filter then satisfies

$$\begin{array}{c}yt={\int }_{-\infty }^{\infty }u\tau gt-\tau d\tau ,\\ yt\le {\int }_{-\infty }^{\infty }u\tau ·gt-\tau d\tau \le c{g}_{L^1}\le c,\end{array}$$

(4)

where the last inequality assumes that $g_{L^1}\le 1$, which is introduced to ensure that $y\left(t\right)\le c$, which in turn allows $y\left(t\right)$ to be sampled by the same neuron. According to the properties of the convolution operator, we also have that $y\in L^2\left(\mathbb{R}\right)\cap L^1\left(\mathbb{R}\right)$, where $L^2\left(\mathbb{R}\right)$ denotes the space of square integrable functions.

Let ${\left\{t_k^u\right\}}_{k\in \mathbb{Z}}$ and ${\left\{t_k^y\right\}}_{k\in \mathbb{Z}}$ be the neuromorphic samples of signals u and y, respectively, computed using an IF neuron with parameters $\delta ,b$. The proposed problem is to compute $t_k^y$ knowing $t_k^u$, sampling parameters $\delta ,b$ and filter $g\left(t\right)$. This problem, illustrated in Figure 1, is inspired from digital signal processing, where a digital filter is applied directly to the samples of a signal. The conventional way to address this problem would be to recover $u\left(t\right)$ from $t_k^u$, apply filter $g\left(t\right)$ in the analog domain, and subsequently sample the output $y\left(t\right)$ with the same IF model to get $t_k^y$. We will refer to this as the indirect method for filtering. The first step of recovery, however, is not possible unless we impose some further restrictive conditions on $u\left(t\right)$ such as being bandlimited [20], shift-invariant [11,17], or having a finite rate of innovation [16,18]. Therefore, the proposed problem is not solvable it its full generality using conventional approaches.

However, if we replace neuromorphic sampling by conventional uniform sampling, this problem leads to the widely used operation of digital filtering. The operation itself does not require any special conditions on the analog signal that generated these samples. Therefore, an equivalent of this solution for the case of neuromorphic sampling is highly desirable.

4. The Proposed Neuromorphic Direct Filtering Method

In this section we describe the proposed direct filtering method. To compute $t_k^y$ from $t_k^u$ we need to create an analytical link between the integrals of the underlying analog signals $u\left(t\right)$ and $y\left(t\right)$ due to the integral operator in (6). To this end, we define the following auxiliary functions

$$U\left(t\right)\triangleq {\int }_0^{t}u\left(\tau \right)d\tau ,Y\left(t\right)\triangleq {\int }_0^{t}y\left(\tau \right)d\tau .$$

(5)

We note that U satisfies $\left|U\left(t\right)\right|\le {∥u∥}_{L^1}$. Using Young’s convolution inequality, we get that $\left|Y\left(t\right)\right|\le {∥y∥}_{L^1}\le {∥u∥}_{L^1}·{∥g∥}_{L^1}$. Using these functions, we derive the equivalent of the t-transform equations (1) as

$$U\left(t_k^u\right)=k\delta -b{t}_k^u,Y\left(t_k^y\right)=k\delta -b{t}_k^y,\forall k\in \mathbb{Z}.$$

(6)

Assuming we know $Y\left(t\right)$, the target spike train $t_k^y$ satisfies $f_k\left(t_k^y\right)=t_k^y$, where

$$f_k\left(t\right)\triangleq \frac{1}{b}k\delta -Y\left(t\right).$$

(7)

The following result shows that $t_k^y$ can be uniquely computed using $f_k\left(t\right)$.

Lemma 1

(Exact Output Samples Computation) Function $f_k\left(t\right)$ has a unique fixed point $t=t_k^y$. Furthermore, let $t_{k,m}^y$ be computed recursively such as $t_{k,0}^y\in$ is arbitrary and

$$t_{k,m+1}^y=f_k{t}_{k,m}^y.$$

(8)

Then ${lim}_{m\to \infty }{t}_{k,m}^y=t_k^y$ and $t_{k,m}^y-t_k^y<t_k^y-t_{k,0}^y{\frac{c}{b}}^m$.

Proof. 

Assume, by contradiction, that $\exists \overline{t}\ne t_k^y$ such that $f_k\left(\overline{t}\right)=\overline{t}$. It follows that

$$b\overline{t}+Y\left(\overline{t}\right)={\int }_0^{\overline{t}}y\left(\tau \right)+bd\tau =k\delta .$$

On the other hand, we know that $y\left(t\right)\le c<b$ due to (4), and thus ${\int }_0^{\overline{t}}y\left(\tau \right)+bd\tau$ is a strictly increasing function of $\overline{t}$, which ensures that ${\int }_0^{\overline{t}}y\left(\tau \right)+bd\tau =k\delta$ has a unique solution. Using (6), we get that $\overline{t}=t_k^y$, which invalidates our initial assumption and proves uniqueness.

From (6), it follows that $\forall k\in \mathbb{Z},f_k\left(t_k^y\right)=t_k^y.$ The following holds

$$\left|f_k^{{}^{\prime }}\left(t\right)\right|=\left|-\frac{y\left(t\right)}{b}\right|\le \frac{c}{b},$$

and thus

$$\forall \zeta \ge 0,\forall t\in \left[t_k^y-\zeta ,t_k^y+\zeta \right],f_k\left(t\right)\in \left[t_k^y-\zeta \frac{c}{b},t_k^y+\zeta \frac{c}{b}\right].$$

(9)

Let $\zeta =t_k^y-t_{k,0}^y$ and $t=t_{k,0}^y$. It follows that

$$f_k\left(t_{k,0}^y\right)-t_k^y=t_{k,1}^y-t_k^y\le t_k^y-t_{k,0}^y\frac{c}{b}.$$

(10)

Similarly, by choosing $\zeta =t_k^y-t_{k,1}^y$ and $t=t_{k,1}^y$ we get

$$t_k^y-t_{k,2}^y\le t_k^y-t_{k,1}^y\frac{c}{b}\le t_k^y-t_{k,0}^y{\frac{c}{b}}^2,$$

and the process continues recursively, which completes the proof via $c<b$. □

Therefore $t_k^y$ can be computed by solving the fixed-point equation $f_k\left(t\right)=t$. This equation requires knowing $Y\left(t\right)$, which satisfies

$$\begin{array}{cc}\hfill Y\left(t\right)=& {\int }_0^t{\int }_{-\infty }^{\infty }g\left(\tau \right)u\left(s-\tau \right)d\tau ds={\int }_{-\infty }^{\infty }g\left(\tau \right){\int }_0^{t}u\left(s-\tau \right)dsd\tau \hfill \\ \hfill =& {\int }_{-\infty }^{\infty }g\left(\tau \right)[U(t-\tau )-U(-\tau )]d\tau ={\int }_{-\infty }^{\infty }U\left(\tau \right)[g(t-\tau )-g(-\tau )]d\tau ,\hfill \end{array}$$

where the last equality uses the variable change $\tau \to t-\tau$. In reality, however, $Y\left(t\right)$ is unknown, since it could only be precisely computed using $U\left(t\right)$. Given that we only know $t_k^u$, and don’t impose any smoothness or sparsity conditions on $u\left(t\right)$, we don’t have access to $U\left(t\right)$, but only to its samples $U{t}_k^u$ via (6). In the following, we show how $Y\left(t\right)$, $f_k\left(t\right)$, and subsequently $t_k^y$ can be estimated using a piecewise constant approximation of $U\left(t\right)$ at points $t_k^u$. Let ${\tilde{f}}_k:\mathbb{R}\to \mathbb{R},$ be defined by

$${\tilde{f}}_k\left(t\right)=\frac{1}{b}\left(k\delta -\tilde{Y}\left(t\right)\right),$$

(11)

where

$$\tilde{Y}\left(t\right)={\int }_{-\infty }^{\infty }{I}_{1}U\left(\tau \right)\left[g\left(t-\tau \right)-g\left(-\tau \right)\right]d\tau ,$$

(12)

where $I_{1}U\left(t\right)$ is the piecewise constant interpolant to U at points ${\left\{t_k^u\right\}}_{k\in \mathbb{Z}}$, such that $I_{1}U\left(t\right)=U{t}_k^u$ for $t\in \left[t_k^u,t_{k+1}^u\right)$. The next proposition derives some properties of $\tilde{Y}\left(t\right)$.

Proposition 1.

Function $\tilde{Y}\left(t\right)$ is continuous and satisfies

$$\tilde{Y}\left(t\right)=\sum _{k\in \mathbb{Z}}U\left(t_k^u\right)\left(G(t-t_k^u)+G(-t_k^u)\right)-U\left(t_k^u\right)\left(G(t-t_{k+1}^u)-G(-t_{k+1}^u)\right),$$

(13)

where $G\left(t\right)\triangleq {\int }_0^{t}g\left(s\right)ds$.

Proof. 

Using (12), we get

$$\begin{array}{cc}\hfill \tilde{Y}\left(t\right)& =\sum _{k\in }{\int }_{t_k^u}^{t_{k+1}^u}{I}_{1}U\left(\tau \right)\left[g\left(t-\tau \right)-g\left(-\tau \right)\right]d\tau \hfill \\ \hfill & =\sum _{k\in }U{t}_k^u{\int }_{t_k^u}^{t_{k+1}^u}g(t-\tau )d\tau -{\int }_{t_k^u}^{t_{k+1}^u}g(-\tau )d\tau \hfill \\ \hfill & =\sum _{k\in }U{t}_k^u{\int }_{t-t_{k+1}^u}^{t-t_k^u}g\left(\tau \right)d\tau -{\int }_{-t_{k+1}^u}^{-t_k^u}g\left(\tau \right)d\tau ,\hfill \end{array}$$

(14)

which proves (13). It follows that $\tilde{Y}\left(t\right)$ is a linear combination of continuous functions, thus it is itself continuous. □

Proposition Section 4 shows that, unlike $Y\left(t\right)$ and $f_k\left(t\right)$, function $\tilde{Y}\left(t\right)$ and, consequently, also ${\tilde{f}}_k\left(t\right)$ are fully known from the IF parameters and input samples $t_k^u$. The remaining challenge is to show that the fixed point equation ${\tilde{f}}_k\left(t\right)=t$ can be solved and to provide an error bound for estimating $t_k^y$. This challenge is addressed rigorously in the next theorem. Moreover, the result allows computing recursively a sequence of estimations ${\tilde{t}}_{k,m}^y$ that converges to a vicinity of $t_k^y$.

Theorem 1

(Estimating Output Samples from Input Samples) Let $u\left(t\right)$ be a signal satisfying $u\in L^1\left(\mathbb{R}\right)\cap L^2\left(\mathbb{R}\right)\cap C\left(\mathbb{R}\right),\left|u\left(t\right)\right|\le c<b$. Furthermore, let $g\left(t\right)$ be the impulse response of a filter satisfying $g\in L^1\left(\mathbb{R}\right)\cap C\left(\mathbb{R}\right),{\parallel g\parallel }_{L^1}\le 1$, and let $y\left(t\right)$ be the filter output in response to input $u\left(t\right)$. Signals $u\left(t\right)$ and $y\left(t\right)$, sampled with an IF neuron with parameters $\delta ,b$, generate values ${\left\{t_k^u\right\}}_{k\in \mathbb{Z}}$ and ${\left\{t_k^y\right\}}_{k\in \mathbb{Z}}$, respectively. Then the following hold true

(a) 

$\forall k\in \mathbb{Z},\exists {\tilde{t}}_k^y\in \left[t_k^y-\frac{2\delta c}{{b-c}^2},t_k^y+\frac{2\delta c}{{b-c}^2}\right]$such that${\tilde{f}}_k\left({\tilde{t}}_k^y\right)={\tilde{t}}_k^y,$where${\tilde{f}}_k\left(t\right)$satisfies (11).

(b) 

Let$k,m$be a sequence defined recursively as$k,m+1={\tilde{f}}_k\left(k,m\right)$, where$k,0\in \mathbb{R}$. Then

$$k,m-t_k^y\le t_k^y-k,0·{\left(\frac{c}{b}\right)}^m+\frac{2\delta c}{b(b-c)}·\sum _{i=1}^m{\left(\frac{c}{b}\right)}^{i-1}$$

(15)

(c) 

For$k,m$defined above,$\exists m_0\in$such that$k,m\in \left[t_k^y-\frac{2\delta c}{{b-c}^2},t_k^y+\frac{2\delta c}{{b-c}^2}\right],\forall m>m_0$.

Proof. 

(a) Function $\tilde{Y}$ satisfies

$$\begin{array}{c}\left|Y\left(t\right)-\tilde{Y}\left(t\right)\right|\le {\int }_{-\infty }^{\infty }\left|U\left(\tau \right)-I_{1}U\left(\tau \right)\right|·\left|g\left(t-\tau \right)-g\left(-\tau \right)\right|d\tau \\ \le \Delta ·\underset{\tau \in \mathbb{R}}{sup}\left|U^{\prime }\left(\tau \right)\right|·{\int }_{-\infty }^{\infty }\left|g\left(t-\tau \right)-g\left(-\tau \right)\right|d\tau \le E,\end{array}$$

(16)

where $E=2\Delta c{∥g∥}_{L^1},\Delta =su{p}_{k\in \mathbb{Z}}\left(t_{k+1}^u-t_k^u\right).$ From (9), the following holds

$$\forall \zeta \ge 0,\forall t\in \left[t_k^y-\zeta ,t_k^y+\zeta \right],f_k\left(t\right)\in \left[t_k^y-\zeta \frac{c}{b},t_k^y+\zeta \frac{c}{b}\right].$$

(17)

Using that $\left|f_k\left(t\right)-{\tilde{f}}_k\left(t\right)\right|\le \frac{E}{b},\forall t\in \mathbb{R}$,

$$\forall \zeta \ge 0,\forall t\in \left[t_k^y-\zeta ,t_k^y+\zeta \right],\widetilde{f_k}\left(t\right)\in \left[t_k^y-\zeta \frac{c}{b}-\frac{E}{b},t_k^y+\zeta \frac{c}{b}+\frac{E}{b}\right].$$

(18)

Unlike in the case of Lemma Section 4, in this case, applying $\widetilde{f_k}\left(t\right)$ recursively does not guarantee exact computation of $t_k^y$. However, we observe that by picking $\zeta =\frac{E}{b-c}$, we get identical intervals for t and $\widetilde{f_k}\left(t\right)$:

$$\forall t\in \left[t_k^y-\frac{E}{b-c},t_k^y+\frac{E}{b-c}\right],\widetilde{f_k}\left(t\right)\in \left[t_k^y-\frac{E}{b-c},t_k^y+\frac{E}{b-c}\right].$$

(19)

This observation is very useful as it enables applying Brouwer fixed-point theorem which states that for any continuous function $f:\mathbb{S}\to \mathbb{S}$, where $\mathbb{S}$ is a nonempty compact convex set, there is a point $t_0$ such that $f\left(t_0\right)=t_0$. Given that $\tilde{Y}\left(t\right)$ is continuous due to Proposition Section 4, it follows that ${\tilde{f}}_k$ is also continuous. By applying Brouwer’s fixed point theorem for $f\left(t\right)={\tilde{f}}_k\left(t\right)$ and $\mathbb{S}=\left[t_k^y-\frac{E}{b-c},t_k^y+\frac{E}{b-c}\right]$, it follows that $f_k\left(t\right)$ has a fixed point in $\mathbb{S}$. We recall that $E=2\Delta c{g}_{L^1}$, which, using (2), leads to

$$E\le \frac{2\delta c}{b-c}{g}_{L^1}\le \frac{2\delta c}{b-c}.$$

(20)

It follows that $\mathbb{S}\subseteq \left[t_k^y-\frac{2\delta c}{{b-c}^2},t_k^y+\frac{2\delta c}{{b-c}^2}\right]$, which yields the required result.

(b) We approach this proof using mathematical induction. We select $t=k,0$ in (18). Using (20), it follows that

$$k,0-t_k^y\le \zeta \frac{c}{b}+\frac{2\delta c}{b(b-c)},\mathrm{s}.\mathrm{t}.k,0\in t_k^y-\zeta ,t_k^y+\zeta .$$

(21)

We note that $k,0\in t_k^y-\zeta ,t_k^y+\zeta$ is always true for $\zeta =t_k^y-k,0$, which yields

$$k,0-t_k^y\le t_k^y-k,0\frac{c}{b}+\frac{2\delta c}{b(b-c)}.$$

(22)

This demonstrates that (15) is true for $m=1$. To finalise the induction, we assume (15) to be true, and show it’s true for $m+1$ as follows.

$$k,m+1-t_k^y=k,m-t_k^y\le \zeta \frac{c}{b}+\frac{2\delta c}{b(b-c)}.$$

(23)

Finally, as before, we use that $\zeta =t_k^y-k,m$ guarantees $k,m\in t_k^y-\zeta ,t_k^y+\zeta$. We also use that $\zeta$ is bounded by (15), which, when substituted in (23), leads to the desired result via (20).

(c) Equation (15) can be expanded into

$$\begin{array}{c}\hfill {\tilde{t}}_{k,m}^y\le t_k^y+t_k^y-k,0{\left(\frac{c}{b}\right)}^m+\frac{2\delta c}{{(b-c)}^2}·\left(1-{\left(\frac{c}{b}\right)}^m\right),\\ \hfill {\tilde{t}}_{k,m}^y\ge t_k^y-t_k^y-k,0{\left(\frac{c}{b}\right)}^m-\frac{2\delta c}{{(b-c)}^2}·\left(1-{\left(\frac{c}{b}\right)}^m\right).\end{array}$$

The required result follows from ${lim}_{m\to \infty }{\frac{c}{b}}^m=0$ via (20). □

Theorem Section 4 shows one can construct a sequence $k,m$ that approximates $t_k^y$ with an error $\frac{2\delta c}{{b-c}^2}$. We note that this error is dependent only on the neuron parameters, and thus can be made arbitrarily small by changing the IF model. In practice, we are using a finite sequence of input measurements ${\left\{t_l^u\right\}}_{l=0,\cdots ,L}$ to approximate the output samples ${\left\{t_k^y\right\}}_{k=0,\cdots ,K}$ that satisfy $t_0^u\le t_k^y\le t_L^u,\forall k=0,\cdots ,K$. The proposed recovery steps are summarized in Algorithm 1.

• Algorithm 6.1.
• Step 1. Set $k=1$ and $0=0$. While $k-1<t_L^u$,
•    Step 1a. Calculate

  $$\tilde{Y}\left({\tilde{t}}_{k-1}^y\right)=\sum _{l=0}^{L}U\left(t_l^u\right)\left(G({\tilde{t}}_{k-1}^y-t_l^u)-G(-t_l^u)\right)-U\left(t_l^u\right)\left(G({\tilde{t}}_{k-1}^y-t_{l+1}^u)+G(-t_{l+1}^u)\right),$$

  (24)

     where $G\left(t\right)={\int }_0^{t}g\left(\tau \right)d\tau ,U\left(t_l^u\right)=l\delta -b{t}_l^u,l=0,\cdots ,L;$
•    Step 1b. $k,0={\tilde{t}}_{k-1}^y$.
•    Step 1c.  Compute $k,m+1=k,m$ for $m=1,\cdots ,M$, where ${\tilde{f}}_k\left(t\right)=\frac{1}{b}\left(\delta -\tilde{Y}\left(t\right)\right)$.
•    Step 1d. ${\tilde{t}}_k^y=k,M$.
•    Step 1e. $k=k+1$;
• Step 2. Set $K=k-2$.

We note that, in practice, convergence was achieved in Algorithm 1 for $M\le 4$ iterations in all examples we evaluated. In the next section, we evaluate numerically the proposed algorithm.

5. Numerical Study

In order to evaluate the performance of the proposed method we define the following approximation errors

• Output inter spike time error between the output of the time-encoded filter and the time-encoded output of the conventional filter:

  $$e_j^t\left(k\right)=t_{k+1}^y-t_k^y-{\tilde{t}}_{k+1}^{y_j}-{\tilde{t}}_k^{y_j},\forall k=1,\cdots ,K,$$

  where ${\tilde{t}}_k^{y_j}$ is the output spike train prediction with Algorithm 1 ($j=1$) and the conventional method, involving the reconstruction of $u\left(t\right)$ ($j=2$).
• Output error between the decoded output of the time-encoded filter and the output of the conventional digital filter:

  $$e_j^y\left(t\right)=y\left(t\right)-{\tilde{y}}_j\left(t\right),\forall t\in [0,2],j=1,2.$$

We compute the errors as percentage of the true quantities $t_{k+1}^y-t_k^y$ and $y\left(t\right)$, respectively, by defining the following normalized error metrics

$$\begin{array}{c}\hfill {\mathrm{ERR}}_j^t=100·\frac{{e_j^t}_{{\ell }^2}}{\sqrt{{\sum }_{k=1}^K{t_{k+1}^{y_j}-t_k^{y_j}}^2}}(\%),{\mathrm{ERR}}_j^y=100·\frac{{e_j^y}_{L^2}}{y_{L^2}}(\%).\end{array}$$

Furthermore, we denote by ${\mathrm{Time}}_1$ and ${\mathrm{Time}}_2$ the computation times of the output spike train with the proposed and indirect methods, respectively. In our numerical examples, we discretised the time using a simulation time step of $\epsilon =1·{10}^{-4}\mathrm{s}$. The parameters used, errors and computing times are summarised in Table 1.

5.1. Low-Pass Filtering of a Bandlimited Signal

The filter selected has the following transfer function

$$G\left(s\right)=\frac{{\omega }_c}{s+{\omega }_c},$$

(25)

where ${\omega }_c=10\mathrm{rad}/\mathrm{s}$. Furthermore, the input is generated as

$$u\left(t\right)=\sum _{k=1}^N{c}_k·\frac{sin\left(\mathsf{\Omega }\right(t-k\pi /\mathsf{\Omega }\left)\right)}{\mathsf{\Omega }(t-k\pi /\mathsf{\Omega })},t\in \left[0,2\right],$$

(26)

where $N=191$ and $c_k$ are random coefficients drawn from the uniform distribution on $(-1,1)$, and $\mathsf{\Omega }=100\mathrm{rad}/2$. We selected $c_i=0,i\in 1,2,3,N-2,N-1,N$ in order to avoid boundary errors for input recovery. The input and output of system $G\left(s\right)$ are sampled with an integrate and fire neuron with parameters $\delta =0.05,b=6$. Signals $u\left(t\right),y\left(t\right)$ and corresponding spike times are depicted in Figure 2 (a) and Figure 2 (b), respectively.

We estimated spike times $t_k^y$ in two different ways, via the indirect method presented in Section 3, and the proposed direct method in Section 4. The results are depicted in fig:spikepred, and the recovery errors are ${\mathrm{ERR}}_1^t=0.49\%$ for the proposed method and ${\mathrm{ERR}}_2^t=0.75\%$ for the indirect method. Furthermore, the filter output errors are ${\mathrm{ERR}}_1^y=0.82\%$ and ${\mathrm{ERR}}_2^y=3.35\%$. However, the computing times are $0.62\mathrm{s}$ for the proposed method and $2.93\mathrm{s}$ for the indirect method. We note that a major contribution in the error induced by the conventional method is due to boundary errors, which are a known artefact of reconstruction [10,21].

5.2. Low-Pass Filtering of Uniform White Noise

To show that the proposed method has no limitations in terms of input bandwidth, as indicated by the bounds in Theorem Section 4, here we choose $u\left(t\right)$ to be drawn from the random uniform distribution on $[-1.5,1.5]$. Given that in this case the filter $G\left(s\right)$ reduces significantly the noise amplitude, we amplify the output by 10 by using the new filter $G\left(s\right)=10·\frac{{\omega }_c}{s+{\omega }_c}$, where ${\omega }_c=10\mathrm{rad}/\mathrm{s}$. The neuron parameters are $\delta =0.1,b=6$. The input and output signals, as well as the output recovered with both methods are in Figure 4a,b. The proposed prediction method does not need any knowledge of the input bandwidth. As an alternative, we reconstructed the input with the conventional indirect method, where the recovery bandwidth is $\mathsf{\Omega }=100\mathrm{rad}/\mathrm{s}$. The value of $\mathsf{\Omega }$ is well above the cuttoff frequency of the filter ${\omega }_c$, which would lead to a accurate recovery of the filter output. However, the random input $u\left(t\right)$ cannot be recovered for any $\mathsf{\Omega }>0$.

The output spike train prediction with both methods is in Figure 5a. The corresponding errors are ${\mathrm{ERR}}_1^t=0.6\%,{\mathrm{ERR}}_2^t=1.06\%,{\mathrm{ERR}}_1^y=25.22\%,$ and ${\mathrm{ERR}}_2^y=26.68\%$, also summarised in Table 1. The computing times are $0.32\mathrm{s}$ for the proposed method and $0.81\mathrm{s}$ for the indirect method.

5.3. The Effect of Spike Density on Performance

During our experiments, we noticed that the neuron threshold $\delta$ has an important role in deciding the computation time of both prediction routines. To illustrate this, we repeated the experiment in Section 5.2 using the same settings where $\delta$ was iteratively assigned 10 values uniformly spaced in interval $[0.02,0.1]$. As before, we measured the error on the predicted spike train, the reconstructed filter output and computation time for each method, illustrated in Figure 5b, Figure 6a and Figure 6b, respectively.

5.4. Output Spikes Prediction for Neurons Sampling Below Nyquist Rate

In the examples presented so far, the filter output is sampled at a rate that allows input reconstruction. In this example, we use a uniform white noise input bounded in $[-1.5,1.5]$. The input is filtered using a bandpass filter computed using the Daubechies orthogonal wavelet 30, depicted in Figure 7. We note that the effective bandwidth is around $1,200\mathrm{rad}/\mathrm{s}$. When sampled with an IF with parameters $\delta =0.04,b=6$, the filter output samples satisfy $max{t}_{k+1}^y-t_k^y=0.0074>\frac{\pi }{\mathsf{\Omega }}=0.0026$, which violates the condition for perfect recovery in (3). Therefore, in this example, we will only evaluate the error using ${\mathrm{ERR}}_j^t$, as ${\mathrm{ERR}}_j^y$ involves the reconstruction of the filter output from IF measurements.

We note that, as in Section 5.2 and Section 5.3, here the input is not bandlimited, and its recovery cannot be guaranteed by (3). In this case, the choice of the bandwidth parameter $\mathsf{\Omega }$ used in recovery is not determined by any existing theoretical result. We run the conventional method by changing the recovery bandwidth across 11 uniformly distributed values in interval $500\mathrm{rad}/\mathrm{s},1500\mathrm{rad}/\mathrm{s}$. The comparative results are depicted in Figure 8a,b. The results show that the error introduced by the conventional method ${\mathrm{ERR}}_2^t$ changes quite significantly with $\mathsf{\Omega }$. An error around $15\%$ higher than the proposed method is achieved for $\mathsf{\Omega }\in 600,1300$, but the difference becomes significantly larger for $\mathsf{\Omega }$ outside this interval. Meanwhile, the computation time of the proposed method is more than one order of magnitude lower than the conventional method (see Figure 8b) for all values of $\mathsf{\Omega }$.

6. Conclusions

In this work we introduced a new method to filter an analog signal via its neuromorphic measurements. Unlike existing approaches, the method does not require imposing smoothness type assumptions on the analog input and filter output, such as a limited bandwidth. We introduced recovery guarantees, showing that it is possible to approximate the output spike train with arbitrary accuracy for an appropriate choice of the sampling model. We compared the proposed method numerically against the conventional solution to this problem, which involves reconstruction of the analog signal. The results show the accuracy of the proposed method is comparable to that of the conventional approach. However, the computing time was smaller for the proposed method in all examples, ranging from 2-3 times up to more than one order of magnitude smaller.

Conceptually, the proposed method has the advantage that it does not depend on the characteristics of the analog signal, and therefore it is not restricted to satisfy any reconstruction guarantees. As demonstrated numerically, the method works well in the case of random inputs, as well as when the input and output of the filter are sampled below Nyquist. Moreover, given it bypasses input reconstruction, the proposed method is not affected by known artefacts of recovery methods such as boundary errors. This work has the potential to lead to the development of neuromorphic filters, that would facilitate a faster transition towards a power efficient computing infrastructure.