Year: 2020

Abstract:
In many animal sensory pathways, the transformation from external stimuli to spike trains is essentially deterministic. In this context, a new mathematical framework for coding and reconstruction, based on a biologically plausible model of the spiking neuron, is presented. The framework considers encoding of a signal through spike trains generated by an ensemble of neurons via a standard convolve-thenthreshold mechanism, albeit with a wide variety of convolution kernels. Neurons are distinguished by their convolution kernels and threshold values. Reconstruction is posited as a convex optimization minimizing energy. Formal conditions under which perfect and approximate reconstruction of the signal from the spike trains is possible are then identified. Coding experiments on a large audio dataset are presented to demonstrate the strength of the framework.

text: In many animal sensory pathways, the transformation from external stimuli to spike trains is essentially deterministic. In this context, a new mathematical framework for coding and reconstruction, based on a biologically plausible model of the spiking neuron, is presented. The framework considers encoding of a signal through spike trains generated by an ensemble of neurons via a standard convolve-thenthreshold mechanism, albeit with a wide variety of convolution kernels. Neurons are distinguished by their convolution kernels and threshold values. Reconstruction is posited as a convex optimization minimizing energy. Formal conditions under which perfect and approximate reconstruction of the signal from the spike trains is possible are then identified. Coding experiments on a large audio dataset are presented to demonstrate the strength of the framework.

heading: 1 INTRODUCTION
text: In biological systems, sensory stimuli is communicated to the brain primarily via ensembles of discrete events that are spatiotemporally compact electrical disturbances generated by neurons, otherwise known as spikes. Spike train representation of signals, when sparse, are not only intrinsically energy efficient, but can also facilitate downstream computation(6; 10). In their seminal work, Olshausen and Field (13) showed how efficient codes can arise from learning sparse representations of natural stimulus statistics, resulting in striking similarities with observed biological receptive fields. (19) developed a biophysically motivated spiking neural network which for the first time predicted the full diversity of V1 simple cell receptive field shapes when trained on natural images. Although these results signify substantial progress, an effective end to end signal processing framework that deterministically represents signals via spike train ensembles is yet to be laid out. Here we present a new framework for coding and reconstruction leveraging a biologically plausible coding mechanism which is a superset of the standard leaky integrate-and-fire neuron model (5).
Our proposed framework identifies reconstruction guarantees for a very general class of signals—those with finite rate of innovation (18)—as shown in our perfect and approximate reconstruction theorems. Most other classes, e.g. bandlimited signals, are subsets of this class. The proposed technique first formulates reconstruction as an optimization that minimizes the energy of the reconstructed signal subject to consistency with the spike train, and then solves it in closed form. We then identify a general class of signals for which reconstruction is provably perfect under certain ideal conditions. Subsequently, we present a mathematical bound on the error of an approximate reconstruction when the model deviates from those ideal conditions. Finally, we present simulation experiments coding for a large dataset of audio signals that demonstrate the efficacy of the framework. In a separate set of experiments on a smaller subset of audio signals we compare our framework with existing sparse coding algorithms viz matching pursuit and orthogonal matching pursuit, establishing the strength of our technique.
The remainder of the paper is structured as follows. In Sections 2 and 3 we introduce the coding and decoding frameworks. Section 4 identifies the class of signals for which perfect reconstruction is achievable if certain ideal conditions are met. In Section 5 we discuss how in practice those ideal conditions can be approached and provide a mathematical bound for approximate reconstruction. Simulation results are presented in Section 6. We conclude in Section 8.

heading: 2 CODING
text: The general class of deterministic mappings (i.e., the set of all nonlinear operators) from continuous time signals to spike trains is difficult to characterize because the space of all spike trains does not lend itself to a natural topology that is universally embraced. The result is that simple characterizations, such as the set of all continuous operators, can not be posited in a manner that has general consensus. To resolve this issue, we take a cue from biological systems. In most animal sensory pathways, external stimulus passes through a series of transformations before being turned into spike trains(17). For example, visual signal in the retina is processed by multiple layers of non-spiking horizontal, amacrine and bipolar cells, before being converted into spike trains by the retinal ganglion cells. Accordingly, we can consider the set of transformations that pass via an intermediate continuous time signal which is then transformed into a spike train through a stereotyped mapping where spikes mark threshold crossings. The complexity of the operator now lies in the mapping from the continuous time input signal to the continuous time intermediate signal. Since any time invariant, continuous, nonlinear operator with fading memory can be approximated by a finite Volterra series operator(2), this general class of nonlinear operators from continuous time signals to spike trains can be modeled as the composition of a finite Volterra series operator and a neuronal thresholding operation to generate a spike train. Here, the simplest subclass of these transformations is considered: the case where the Volterra series operator has a single causal, bounded-time, linear term, the output of which is composed with a thresholding operation of a potentially time varying threshold. The overall operator from the input signal to the spike train remains nonlinear due to the thresholding operation. The code generated by an ensemble of such transformations, corresponding to an ensemble of spike trains, is explored.
Formally, we assume the input signal X(t) to be a bounded square integrable function over the compact interval [0, T ] for some T ∈ R+, i.e., we are interested in the class of input signals F = {X(t)|X(t) ∈ L2[0, T ]}. Since the framework involves signal snippets of arbitrary length, this choice of T is without loss of generalization. We assume an ensemble of convolution kernels K = {Kj |j ∈ Z+, j ≤ n}, consisting of n kernels Kj , j = 1, . . . , n. We assume that Kj(t) is a continuous function on a bounded time interval [0, T ], i.e. ∀j ∈ {1, . . . , n},Kj(t) ∈ C[0, T ], T ∈ R+. Finally, we assume that Kj has a time varying threshold denoted by T j(t). The ensemble of convolution kernels K encodes a given input signal X(t) into a sequence of spikes {(ti,Kji)}, where the ith spike is produced by the jthi kernel Kji at time ti if and only if:∫ X(τ)Kji(ti − τ)dτ = T ji(ti) In our experiments a specific threshold function is assumed in which the time varying threshold T j(t) of the jth kernel remains constant at Cj until that kernel produces a spike, at which time an after-hyperpolarization potential (ahp) is introduced to raise the threshold to a high value M j Cj , which then drops back linearly to its original value within a refractory period δj . Stated formally,
T j(t) = C j , t− δj > tjl (t)
M j − (t−t j l (t))(M j−Cj) δj , t− δj ≤ tjl (t) (1)
Where tjl (t) denotes the time of the last spike generated by K j prior to time t.

heading: 3 DECODING
text: How rich is the coding mechanism just described? We can investigate this question formally by positing a decoding module. The objective of the decoding module is to reconstruct the original signal from the encoded ensemble of spike trains. It is worthwhile to mention that to be able to communicate signals properly by our proposed framework, the decoding module needs to be designed in a manner so that it can operate solely on the spike train data handed over by the encoding module, without explicit access to the input signal itself. Considering the prospect of the invertibility of the coding scheme, we seek a signal that satisfies the same set of constraints as the original signal when generating all spikes apropos the set of kernels in ensemble K. Recognizing that such a signal might not be unique, we choose the reconstructed signal as the one with minimum L2-norm.
Formally, the reconstruction (denoted by X∗(t)) of the input signal X(t) is formulated to be the solution to the optimization problem:
X∗(t) = argmin X̃ ||X̃(t)||22
s.t. ∫ X̃(τ)Kji(ti − τ)dτ = T ji(ti); 1 ≤ i ≤ N (2)
where {(ti,Kji)|i ∈ {1, ..., N}} is the set of all spikes generated by the encoder. The choice of L2 minimization as the objective of the reconstruction problem—which is the linchpin of our framework, as demonstrated in the theorems—can only be weakly justified at the current juncture. The perfect reconstruction theorem that follows provides the strong justification. As it stands, the L2 minimization objective is in congruence with the dictum of energy efficiency in biological systems. The assumption is that, of all signals, the one with the minimum energy that is consistent with the spike trains is desirable. Additionally, an L2 minimization in the objective of (2) reduces the convex optimization problem to a solvable linear system of equations as shown in Lemmas 1 and 3. Later we shall show that L2-minimization has the surprising benefit of recovering the original signal perfectly under certain conditions.

heading: 4 SIGNAL CLASS FOR PERFECT RECONSTRUCTION
text: To establish the effectiveness of the described coding-decoding model, we have to evaluate the accuracy of reconstruction over a class of input signals. We observe that in general the encoding of square integrable signals into spike trains is not a one-to-one map; the same set of spikes can be generated by different signals so as to result in the same convolved values at the spike times. Naturally, with a finite and fixed ensemble of kernelsK, one cannot achieve perfect reconstruction for the general class of signals F as defined in Section 2. We now restrict ourselves to a subset G of the original class F defined as G = {X(t)|X(t) ∈ F , X(t) = ∑N p=1 αpK
jp(tp − t), jp ∈ {1, ..., n}, αp ∈ R, tp ∈ R+, N ∈ Z+} and address the question of reconstruction accuracy. Essentially G consists of all linear combinations of arbitrarily shifted kernel functions. N is bounded above by the total number of spikes that the ensemble K can generate over [0, T ]. In the parlance of signal processing, G constitutes Finite rate of Innovation signals (18). For the class G the perfect reconstruction theorem is presented below. The theorem is proved with the help of three lemmas. Perfect Reconstruction Theorem: Let X(t) ∈ G be an input signal. Then for appropriately chosen time-varying thresholds of the kernels, the reconstruction, X∗(t), resulting from the proposed codingdecoding framework is accurate with respect to the L2 metric, i.e., ||X∗(t)−X(t)||2 = 0. Lemma 1: The solution X∗(t) to the reconstruction problem given by (2) can be written as:
X∗(t) = N∑ i=1 αiK ji(ti − t) (3)
where the coefficients αi ∈ R can be solved from a system of linear equations. Proof: An approach analogous to the Representer Theorem (15), splitting a putative solution to (2) into its within the span of the kernels component and a remnant orthogonal component, results in equation (3). In essence, the reconstructed signal X∗(t) becomes a summation of the kernels, shifted to their respective times of generation of spikes, scaled by appropriate coefficients. Plugging (3) into the constraints (2) gives:
∀1≤i≤N ; ∫ N∑
k=1
αkK jk(tk − t)Kji(ti − t)dτ = T ji(ti)
Setting bi = T ji(ti) and Pik = ∫ Kjk(tk − τ)Kji(ti − τ)dτ results in:
∀1≤i≤N ; N∑ k=1 Pikαk = bi (4)
Equation (4) defines a system of N equations in N unknowns of the form:
Pα = T (5) where α = 〈α1, ..., αN 〉T , T = 〈T j1(t1), ..., T jN (tN )〉T and P is an N ×N matrix with elements Pik = ∫ Kjk(tk − τ)Kji(ti − τ)dτ . Clearly P is the Gramian Matrix of the shifted kernels
{Kji(ti − t)|i ∈ 1, 2, ..., N} in the Hilbert space with the standard inner product. It is well known that P is invertible if and only if {Kji(ti − t)|i ∈ 1, 2, ..., N} is a linearly independent set. If P is invertible α has a unique solution. If, on the other hand, P is not invertible, α has multiple solutions. However, as the next lemma shows, every such solution leads to the same reconstruction X∗(t), and hence any value of α that satisfies 5 can be chosen. We note in passing that in our experiments we have used the least square solution. 2 Import: The goal of the optimization problem is to find the best object in the feasible set. However, the application of the Representer Theorem converts the constraints into a determined system of unknowns and equations, turning the focus onto the feasible set, effectively changing the optimization problem into a solvable system that results in a closed form solution for the αi’s. This implies that instead of solving (2), we can solve for the reconstruction from X∗(t) = ∑N i=1 αiK
ji(ti− t), where αi is the i-th element of α = P−1T . Here, P−1 represents either the inverse or the Moore-Penrose inverse, as the case may be. Lemma 2: Let equation 5 resulting from the optimization problem 2 have multiple solutions. Consider any two different solutions for α, namely α1 and α2, and hence the corresponding reconstructions are given by X1(t) = ∑N i=1 α1iK ji(ti − t) and X2(t) = ∑N i=1 α2iK
ji(ti − t), respectively. Then X1 = X2. Proof: The proof of this lemma follows from the existence of a unique function in the Hilbert Space spanned by {Kji(ti − t)|i ∈ 1, 2, ..., N} that satisfies the constraint of equation 2. The details of the proof is furnished in the appendix A. Import: Lemma 2 essentially establishes the uniqueness of solution to the optimization problem formulated in 2 as any solution to equation 5. The proof follows from the fact that the reconstruction is in the span of the shifted kernels {Kji(ti − t)|i ∈ 1, 2, ..., N} and the inner products of the reconstruction with each of Kji(ti − t) is given (by the spike constraints of 2). Such a reconstruction must be unique in the subspace S. Lemma 3: Let X∗(t) be the reconstruction of an input signal X(t) and {(ti,Kji)}Ni=1 be the set of spikes generated. Then, for any arbitrary signal X̃(t) within the span of {Kji(ti − t)|i ∈ {1, 2, ..., N}}, i.e., the set of shifted kernels at respective spike times, given by X̃(t) =∑N i=1 aiK
ji(ti − t) the following inequality holds: ||X(t)−X∗(t)|| ≤ ||X(t)− X̃(t)|| Proof:
||X(t)− X̃(t)|| = ||X(t)−X∗(t)︸ ︷︷ ︸ A +X∗(t)− X̃(t)︸ ︷︷ ︸ B || First, 〈A,Kji(ti − t)〉 = 〈X(t),Kji(ti − t)〉 − 〈X∗(t),Kji(ti − t)〉,∀i ∈ {1, 2, .., N} = T ji(ti)− T ji(ti) = 0 (Using the constraints in (2) & (2))
Second, 〈A,B〉 = 〈A, N∑ i=1 (αi − ai)Kji(ti − t)〉 (By Lemma 1 X∗(t) = N∑ i=1 αiK ji(ti − t))
= N∑ i=1 (αi − ai)〈A,Kji(ti − t)〉 = 0
=⇒ ||X(t)− X̃(t)||2 = ||A+B||2 = ||A||2 + ||B||2 ≥ ||A||2 = ||X(t)−X∗(t)||2
=⇒ ||X(t)− X̃(t)|| ≥ ||X(t)−X∗(t)|| 2
Import: The implication of the above lemma is quite remarkable. The objective defined in (2) chooses a signal with minimum energy satisfying the constraints, deemed the reconstructed signal. However as the lemma demonstrates, this signal also has the minimum error with respect to the input signal in the span of the shifted kernels. This signifies that our choice of the objective in the decoding module not only draws from biologically motivated energy optimization principles, but also performs optimally in terms of reconstructing the original input signal within the span of the appropriately shifted spike generating kernels.
Corollary: An important consequence of Lemma 3 is that additional spikes in the system do not worsen the reconstruction. For a given input signal X(t) if S1 and S2 are two sets of spike trains where S1 ⊂ S2, the second a superset of the first, then Lemma 3 implies that the reconstruction due to S2 is at least as good as the reconstruction due to S1 because the reconstruction due to S1 is in the span of the shifted kernel functions of S2 as S1 ⊂ S2. This immediately leads to the conclusion that
for a given input signal the more kernels we add to the ensemble the better the reconstruction. Proof of the Theorem: The proof of the theorem follows directly from Lemma 3. Since the input signalX(t) ∈ G, letX(t) be given by: X(t) = ∑N p=1 αpK
jp(tp−t) (αp ∈ R, tp ∈ R+, N ∈ Z+) Assume that the time varying thresholds of the kernels in our kernel ensemble K are set in such a manner that the following conditions are satisfied: 〈X(t),Kjp(tp− t)〉 = T jp(tp) ∀p ∈ {1, ..., N} i.e., each of the kernels Kjp at the very least produces a spike at time tp against X(t) (regardless of other spikes at other times). Clearly then X(t) lies in the span of the appropriately shifted response functions of the spike generating kernels. Applying Lemma 3 it follows that: ||X(t)−X∗(t)||2 ≤ ||X(t)−X(t)||2 = 0 2 Import: In addition to demonstrating the potency of the coding-decoding scheme, this theorem frames Barlow’s efficient coding hypothesis (1)—that the coding strategy of sensory neurons be adapted to the statistics of the stimuli—in mathematically concrete terms. Going by the theorem, the spike based encoding necessitates the signals to be in the span of the encoding kernels for perfect reconstruction. Inverting the argument, kernels must learn to adapt to the basis elements that generate the signal corpora for superior reconstruction.

heading: 5 APPROXIMATE RECONSTRUCTION AND THE EFFECT OF AHP
text: The perfect reconstruction theorem stipulates the conditions under which exact recovery of a signal is feasible in the proposed framework. At first glance, it may seem challenging to meet these conditions for an arbitrary class of natural signals. The concern stems from two difficulties: firstly, given a fixed set of kernels, the input signal may not lie in the span of their arbitrary shifts, and secondly, we may not be able to generate spikes at the desired locations as postulated in the proof of the theorem.
To address these issues, we observe that our decoding model is a continuous transformation from the space of spike trains to L2-functions, in the sense that small changes in spike times or a slight mismatch of the spiking kernels from the components of the signal, bring about only small changes in the reconstruction. In what follows, we furnish an Approximate Reconstruction Theorem (C) that provides a bound on the reconstruction error under such deviations. To address the first problem, it is important to choose kernel functions appropriately so that they can represent the input signals reasonably well. One can leverage biological knowledge; for example, it is well-known that auditory filters are effectively modeled using gammatones (14). Hence our experiments in Section 6 on auditory signals were coded using gammatone kernels. Not surprisingly, the reconstructions were excellent. To alleviate the second problem, we observe that spikes can be produced reasonably close to the desired locations by setting a low baseline threshold and a small refractory period of the after-hyperpolarization potential (ahp) for each kernel, a technique that is guaranteed to give good results as is confirmed by our experiments in Section 6. The following lemma formalizes the notion of how lowering the threshold and the refractory period of a kernel helps in generating spikes at the desired locations. The lemma is followed by the Approximate Reconstruction Theorem. Lemma 4: Let X(t) be an input signal. Let Kp be a kernel for which we want to generate a spike at time tp. Let the inner product 〈X(t),Kp(tp − t)〉 = Ip. Then, if the baseline threshold of the kernel Kp is Cp ≤ Ip and the absolute refractory period is δ as modeled in Equation 1, the kernel Kp must produce a spike in the interval [tp − δ, tp] according to the threshold model defined in Equation 1. Proof: The proof of this lemma follows directly from the intermediate value theorem and is detailed in appendix B. Approximate Reconstruction Theorem: Let the following assumptions be true: • X(t), the input signal to the proposed framework, can be written as a linear combination of some component functions as- X(t) = ∑N i=1 αifpi(ti − t) where αi are bounded real coefficients, the component functions fpi(t) are chosen from a possibly infinite setH = {fi(t)|i ∈ Z+, ||fi(t)|| = 1} of functions of unit L2 norm with compact support, and the corresponding ti ∈ R+ are chosen to be bounded arbitrary time shifts of the component functions so that the overall signal has a compact support in [0, T ] for some T ∈ R+ and thus the input signal still belongs to the same class of signals F , as defined in section 2. • There is at least one kernel Kji from the bag of encoding kernels K, such that the L2-distance of fpi(t) from K
ji(t) is bounded. Formally, ∃ δ ∈ R+ s.t. ||fpi(t)−Kji(t)||2 < δ ∀ i ∈ {1, ..., N}. •When X(t) is encoded by the proposed framework, each one of these kernels Kji produce a spike at time t′i at threshold Ti such that |ti − t′i| < ∆ ∀i, for some ∆ ∈ R+. • Each kernel Kj ∈ K satisfies a Lipschitz type condition as follows:
∃ C∈R s.t. ||Kj(t)−Kj(t−∆t)||2 ≤ C|∆t|, ∀∆t ∈ R, ∀j. • And lastly the shifted component functions satisfy a frame bound type of condition as follows:∑ k 6=i〈fpi(t− ti), fpk(t− tk)〉 ≤ η ∀ i ∈ {1, ..., N} Then, reconstruction X∗(t), resulting from the proposed framework, has a bounded noise to signal ratio. Specifically, the following inequality is satisfied:
||X(t)−X∗(t)||22/||X(t)||22 ≤ (δ + C∆)(1+xmax)/(1−η)
where xmax is a positive number ∈ [0, N − 1] that depends on the maximum overlap of the support of component functions fpi(t− ti). Proof: A detailed proof of this theorem is provided in the appendix C.

heading: 6 EXPERIMENTS ON REAL SIGNALS
text: The proposed framework is general enough to apply to any class of signals. However, since the computational resources necessary to code and reconstruct video signals (function of three variablesx, y, t) would be sufficiently larger than audio signals (function of only one variable t), to demonstrate that the proposed framework can indeed be adopted in real engineering applications as a novel encoding scheme, we ran experiments on a repository of audio signals.

heading: 6.1 DATASET
text: We chose the Freesound Dataset Kaggle 2018 (or FSDKaggle2018 for short), an audio dataset of natural sounds posted on Kaggle referred in (7), containing 18,873 audio files annotated with labels from Google’s AudioSet Ontology (9). For the purpose of the experiments, we ignored the labels and only focused on the sound data, since we were only interested in encoding and decoding the input signals. All audio samples in this dataset are provided as uncompressed PCM 16bit, 44.1kHz, mono audio files, with each file consisting of sound snippets of duration ranging between 300ms to 30s. In the experiment, we ran our proposed methodology over at least 1000 randomly chosen sound snippets from the samples in the dataset. For ease of computation, we kept the length of the input audio snippets to be relatively small (ideally of size less than 50ms), splitting longer signals. This choice of considering small snippets as input made the computation feasible on limited resource machines within reasonable time bounds by reducing the size of the P -matrix referred to in Equation 4. This choice is without loss of generalisation since for encoding signals of greater length, reconstruction using this framework can be done piece-wise: splitting a longer signal into smaller pieces, reconstructing piece-wise and finally stitching the reconstructed pieces together.

heading: 6.2 SET OF KERNELS
text: The proposed encoding technique is operational on a set of kernels, as stated in Equation 2. The first order of business was therefore the choice of a suitable set of kernels for our experiments. Since gammatone filters are widely used as a reasonable model of cochlear filters in auditory systems (14), and mathematically, are fairly simple to represent—atn−1e−2πbt cos(2πft+φ)—in our experiments we chose a set of gammatone filters as our kernels (Figure 1). The implementation of the filterbank is similar to (16), and we used up to 2000 gammatone kernels whose center frequencies were uniformly spaced on the ERB scale between 20 Hz to 20 kHz. In all experiments, the kernels were normalized, and the baseline thresholds and the ahp parameters were kept the same across all kernels.

heading: 6.3 RESULTS
text: Following the assertion of Lemma 4, in all experiments, the baseline threshold and the absolute refractory period were kept low enough so that for each sound snippet near perfect reconstructions could be obtained at a high spike rate. A typical value of the refractory period was ≈ 5ms and the baseline threshold value was kept as low as 10−3. As a consequence of the corollary to Lemma 3, additional spikes did not hurt reconstruction. Experiments were conducted with varying number of kernels. Once a reconstruction at a high spike rate was attained, a greedy technique that removed spikes in order of their impact on the reconstruction was instituted to get a compressed code for each snippet. Reconstructions were then recomputed with the fewer spikes as constraints. We
0.00 0.05 0.10 0.15 0.20 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.00 0.05 0.10 0.000 0.025 0.050 0.075 0.00 0.05
Time in ms
(a) (b) (c) (d) (e)
Figure 1: Five sample gammatone filters used as kernels with center frequencies located at approximately (a) 82 Hz, (b) 118 Hz, (c) 158 Hz, (d) 203 and (e) 253 Hz, respectively.
Amplitude 0
5
10
15
20
25
Ti m e in m
s
0 500 1000 1500 Kernel Indexes
0 200 400 600 800 Kernel Indexes
0 100 200 300 400 500 Kernel Indexes
(a) (b) (c) (d) (e) (f) (g)
Figure 2: Reconstruction of a sample snippet in an experiment with 2000 kernels. (a) the input snippet, extended with zero padding to accommodate future spikes; (b) Spike trains of all kernels obtained at a low threshold and refractory period displayed as a raster plot with time (y-axis) and index of the gammatone kernels in increasing order of center frequencies (x-axis), and (c) the resulting reconstruction. This is an almost perfect reconstruction with a 32.7DB SNR at 1146 kHz spike rate of the ensemble. Subsequently spikes were deleted greedily (see text) to obtain reconstructions at lower spike rates. (d) Resulting spike pattern of the ensemble at 88.4kHz and (e) resulting reconstruction with SNR 19.7DB. Likewise, (f) spike pattern of the ensemble at 17.64kHz and (g) resulting reconstruction with SNR of 9DB. Time scale on left apply to all plots. It is noteworthy that in (d)&(f) spikes of the higher frequency kernels ended up deleted in the culling process.
should emphasize here that soon as spikes are removed to get a compressed representation of the signal, signals are no longer encoded via a simple spike train representation which ideally should communicate only spike times and not their corresponding threshold values. In other words, a compressed signal representation in this scheme needs to communicate both the spike times and the threshold values because once spikes are culled the decoder cannot infer the threshold values from the spike times and the given threshold function of the neurons in equation 1. In that sense a compressed representation of a signal in this approach can be realized through marked spike trains that carry both time as well as threshold information rather than a true spike train based representation. Figure 2 demonstrates this process applied to a sample sound snippet through several stages of removal of spikes. This process was repeated over 1000 randomly chosen sound snippets from the dataset. Figure 3 displays the complete results of the experiment with ≈ 2000 gammatone kernels. As the
figure demonstrates, at high spike rates nearly perfect reconstructions were obtained consistently, and even though lowering the spike rate gradually increased noise, reasonable reconstructions could be obtained at low spike rates (≈ 15DB at 25kHz on average). Since each reconstruction is calculated by solving a system of linear equations involving a P -matrix whose dimension is O(N2) where N is the number of spikes under consideration, computation is fairly time consuming, and the choice of parameters, such as the length of the input snippets, the number of kernels or the threshold parameters were made to ensure feasibility of computation with available resources while maintaining efficacy of the overall reconstruction process.
Since the proposed framework approximates a signal with a sparse linear combination of shifted kernels (as shown by signal class G in 4) and therefore has similarities to compressed sensing, another set of experiments were designed to compare the proposed framework with existing sparse coding techniques viz, Convolutional Matching Pursuit and Convolutional Orthogonal Matching Pursuit. Since the sparse coding techniques are computationally intensive over continuous signals, this set of experiments were restricted to only 10 gammatone kernels (for CMP and COMP all possible shifts of these 10 kernels were considered) and the experiments were run over 30 sound snippets. Our technique was applied as before, starting at a high spike rate with spikes culled gradually to achieve better compression. The results of comparison of average SNR values obtained by the techniques are shown in figure 4. As is evident, our technique in its simplest form does slightly better than COMP up until ≈ 50kHz beyond which COMP performs better. Our technique was ≈ 1.2 times faster than COMP on an i7 hexa-core processor for these experiments and should naturally scale much better than COMP since the proposed technique does not involve repeated computation of inner products. The implementation details of the experiment can be found in our simulation code available at: http://bitbucket.org/crystalonix/oldsensorycoding.git.

heading: 7 RELATION TO PRIOR WORK
text: The problem of representing continuous time signals using ensembles of spike trains has a rich history both in the neuromorphic computing community as well as in computational neuroscience. Most such work rely on classical Nyquist–Shannon sampling theory wherein signals are assumed to be band-limited and reconstruction is realized through sinc filters, albeit via the spike trains. Among existing spike based coding techniques, (4) has explored the spike generating mechanism of the neuron as an oversampling, noise shaping analog-to-digital converter, and (11) represents signals via time encoding machines. Likewise, an image encoding technique has been discussed in (8) using an
integrate-threshold-reset framework that results in spike trains, which leverages differential pulse-code modulation (DPCM) at its core. In our case the input signals considered are elements of L2(R) with finite rate of innovation, the reconstruction error of which tends to zero as the signal approaches the span of appropriately chosen kernel functions which are again a generic class of continuous functions. Our work differs from existing approaches in that using our scheme signal reconstruction is realized via a sparse set of idealized spikes, whereas in the former case signals need to be sampled at a rate higher than the Nyquist rate and reconstruction implicitly relies on sinc interpolation. Since the class of signals G considered in our analysis takes the form X(t) = ∑N p=1 αpK
jp(tp − t), our problem formulation is comparable to that of convolution sparse coding or compressive sensing deconvolution, which, in general, is a hard problem and hence is solved under certain relaxed criteria (3) or by using certain greedy heuristics (12). Our framework provides a corresponding approximate solution to the general problem leveraging a biological thresholding scheme to produce spikes simultaneously at a high rate and then gradually removing unimportant spikes. The proposed technique, therefore, is a novel alternative to existing solutions.

heading: 8 CONCLUSION
text: We have proposed a framework that codes for continuous time signals using an ensemble of spike trains in a manner that is very different from the pulse-density paradigm. The framework applies to all finite rate of innovation signals, which is a very large class that includes bandlimited signals. Although approximate reconstruction is computationally more expensive than interpolation with a sinc kernel (as in Nyquist Shannon), it is feasible, unlike in the case of compressed sensing where the generic case is NP-hard. Fortuitously, the system of linear equations is best solved using the conjugate gradient method since P is a symmetric positive semidefinite matrix. The excellent reconstruction results we have obtained with 2000 kernels—with no parameter tuning—is a testament to the potential of the technique. The human cochlear nerve, in comparison, contains axons of ≈ 50, 000 spiral ganglion cells (corresponding, therefore, to 50,000 kernels). As our theorems show, reconstruction with such a large set of kernels is guaranteed to be even better, albeit at a higher computational cost.

heading: A PROOF LEMMA 2:
text: Lemma 2: Let equation 5 resulting from the optimization problem 2 have multiple solutions. Consider any two different solutions for α, namely α1 and α2, and hence the corresponding reconstructions are given by X1(t) = ∑N i=1 α1iK ji(ti − t) and X2(t) = ∑N i=1 α2iK
ji(ti − t), respectively. Then X1 = X2. Proof: Let S be the subspace of L2-functions spanned by {Kji(ti − t)|i ∈ 1, 2, ..., N} with the standard inner product (by assumption each of {Kji(ti − t)|i ∈ 1, 2, ..., N} are L2-functions and hence S is a subspace of the larger space of all L2-functions). Clearly S is a Hilbert space with dim(S) ≤ N . Hence there exists {e1, ..., eM}, an orthonormal basis of S (where M ≤ N ). Assume that the hypothesis is false, i.e. X1 6= X2. This implies that ∃ ais and bis such that X1(t) = ∑M i=1 aiei and X2(t) = ∑M i=1 biei where not all ais are same as the corresponding bis. =⇒ ∃ k such that ak 6= bk =⇒ 〈X1(t), ek〉 = ak 6= bk = 〈X2(t), ek〉 But ek ∈ span({Kji(ti−t)|i ∈ 1, 2, ..., N}) =⇒ ∃{c1, ..., cN} such that ek = ∑N i=1 ciK
ji(ti−t) =⇒ 〈X1(t), ek〉 = ∑N i=1 ci〈X1(t),Kji(ti − t)〉 = ∑N i=1 ciT
ji(ti) Now, since X1(t), X2(t) are both solutions to the optimization problem 2 =⇒ 〈X1(t), ek〉 = ∑N i=1 ci〈X2(t),Kji(ti − t) = 〈X2(t), ek〉 =⇒ 〈X1(t), ek〉 = 〈X2(t), ek〉 -which contradicts the hypothesis. 2

heading: B PROOF OF LEMMA 4
text: Lemma 4: Let X(t) be an input signal. Let Kp be a kernel for which we want to generate a spike at time tp. Let the inner product 〈X(t),Kp(tp − t)〉 = Ip. Then, if the baseline threshold of the kernel Kp is Cp ≤ Ip and the absolute refractory period is δ as modeled in Equation 1, the kernel Kp must produce a spike in the interval [tp−δ, tp] according to the threshold model defined in Equation 1.
Proof: The lemma is easily proved by contradiction. Assume that prior to and including time tp, the last spike produced by kernel Kp was at time tl. Also assume that tl < tp − δ so that there is no spike in the interval [tp − δ, tp]. Then by Equation 1 the threshold of kernel Kp at time tp is T p(tp) = Cp. But, 〈X(t),Kp(tp − t)〉 = Ip ≥ Cp. Furthermore, since Kp(t) is a continuous function, the convolution C(t) = ∫ X(τ)Kp(t− τ)dτ varies continuously. By Equation 1, as the ahp kicks up the threshold to an arbitrarily high value Mp at time tl and falls linearly to the value Cp before time tp, by the intermediate value theorem, the threshold must be crossed by the convolution C(t) after time tl and before time tp. Since tl was the last spike of kernel Kp before time tp, this is a contradiction. Hence, tl ≮ tp − δ, implying that there must a spike in the interval [tp − δ, tp]. 2

heading: C PROOF OF APPROXIMATE RECONSTRUCTION THEOREM
text: Approximate Reconstruction Theorem: Let the following assumptions be true:
• X(t), the input signal to the proposed framework, can be written as a linear combination of some component functions as below:
X(t) = N∑ i=1 αifpi(ti − t)
where αi are bounded real coefficients, the component functions fpi(t) are chosen from a possibly infinite set H = {fi(t)|i ∈ Z+, ||fi(t)|| = 1} of functions of unit L2 norm with compact support, and the corresponding ti ∈ R+ are chosen to be bounded arbitrary time shifts of the component functions so that overall signal has a compact support in [0, T ] for some T ∈ R+ and thus the input signal still belongs to the same class of signals F , as defined in section 2.
• There is at least one kernel Kji from the bag of encoding kernels K, such that the L2distance of fpi(t) from K
ji(t) is bounded. Formally, ∃ δ ∈ R+ s.t. ||fpi(t)−Kji(t)||2 < δ ∀ i ∈ {1, ..., N}.
• When X(t) is encoded by the proposed framework, each one of these kernels Kji produce a spike at time t′i at threshold Ti such that |ti − t′i| < ∆ ∀i, for some ∆ ∈ R+.
• Each kernel Kj ∈ K satisfies a Lipschitz type of condition as follows: ∃ C∈R s.t. ||Kj(t)−Kj(t−∆t)||2 ≤ C|∆t|, ∀∆t ∈ R, ∀j.
• And lastly the shifted component functions satisfy a frame bound type of condition as follows:∑ k 6=i〈fpi(t− ti), fpk(t− tk)〉 ≤ η ∀ i ∈ {1, ..., N}
Then, reconstruction X∗(t), resulting from the proposed framework, has a bounded noise to signal ratio. Specifically, the following inequality is satisfied:
||X(t)−X∗(t)||22/||X(t)||22 ≤ (δ + C∆)(1+xmax)/(1−η)
where xmax is a positive number ∈ [0, N − 1] that depends on the maximum overlap of the support of component functions fpi(t− ti).
Proof of the Theorem: By hypothesis each kernel Kji produces a spike at time t′i ∀i ∈ {1, ..., N} . Let us call these spikes as fitting spikes. But the coding model might generate some other spikes against X(t) too. Other than the set of fitting spikes {(t′i,Kji)|i ∈ {1, ..., N}}, let {(t̃k,K j̃k)|k ∈ {1, ...,M}} denote those extra set of spikes that the coding model produces for input X(t) against the kernel bag K and call these extra spikes as spurious spikes. Here, M is the number of spurious spikes. By Lemma1 X∗(t) can be represented as below: X∗(t) = ∑N i=1 αiK ji(t′i − t) + ∑M k=1 α̃kK j̃k(t̃k − t)
where αi and α̃k are real coefficients whose values can be formulated again from Lemma1. Let Ti be the thresholds at which kernel Kji produced the spike at time t′i as given in the hypothesis. Hence for generation of the fitting spikes the following condition must be satisfied:
〈X(t),Kji(t′i − t)〉 = Ti∀i ∈ {1, 2, ..., N} (6)
Consider a hypothetical signal Xhyp(t) defined by the equations below:
Xhyp(t) = N∑ i=1 aiK ji(t′i − t), ai ∈ R
s.t.〈Xhyp(t),Kji(t′i − t)〉 = Ti,∀i (7)
Clearly this hypothetical signal Xhyp(t) can be deemed as if it is the reconstructed signal where we are only considering the fitting spikes and ignoring all spurious spikes. Since, Xhyp(t) lies in the span of the shifted kernels used in reconstruction of X(t) using Lemma 3 we may now write:
||X(t)−Xhyp(t)|| ≥ ||X(t)−X∗(t)|| (8) ||X(t)−Xhyp(t)||22 = 〈X(t)−Xhyp(t), X(t)−Xhyp(t)〉 = 〈X(t)−Xhyp(t), X(t)〉 − 〈X(t)−Xhyp(t), Xhyp(t)〉 = 〈X(t)−Xhyp(t), X(t)〉 − ΣNi=1ai〈X(t)−Xhyp(t),Kji(t− t′i)〉 = ||X(t)||22 − 〈X(t), Xhyp(t)〉 (Since by construction〈Xhyp(t),Kji(t− t′i)〉 = Ti∀i ∈ {1...N}) = ΣNi=1Σ
N k=1αiαk〈fi(t− ti), fk(t− tk)〉 − ΣNi=1ΣNk=1αiak〈fi(t− ti),Kjk(t− t′k))〉
= αTFα− αTFKa (9) (denote a = [a1, a2, ..., aN ]T , α = [α1, α2, ..., αN ]T , F = [Fik]N×N , an N ×N matrix, where Fik = 〈fi(t− ti), fk(t− tk)〉 and FK = [(FK)ik]NXN where (FK)ik = 〈fi(t− ti),Kjk(t− t′k)〉)
But using the results of Lemma1 a can be written as:
a = P−1T where P = [Pik]NXN , Pik = 〈Kji(t− t′i),Kjk(t− t′k)〉 And, T = [Ti]N×1 where Ti = 〈X(t),Kji(t− t′i)〉 = ΣNk=1αk〈fk(t− tk),Kji(t− t′i)〉 = FTKα =⇒ a = P−1FTKα
Plugging this expression of a in equations 9 we get,
||X(t)−Xhyp(t)||22 = αTFα− αTFKP−1FTKα (10) But,(FK)ik = 〈fi(t− ti),Kjk(t− t′k)〉 = 〈Kji(t− t′i),Kjk(t− t′k)〉
− 〈Kji(t− t′i)− fi(t− ti),Kjk(t− t′k)〉 = (P )ik − (EK)ik (11)
(denoting EK = [(EK)ik]N×N , where (EK)ik = 〈Kji(t− t′i)− fi(t− ti),Kjk(t− t′k)〉) Also, (F )ik = 〈fi(t− ti), fk(t− tk)〉 = 〈fi(t− ti)−Kji(t− t′i) +Kji(t− t′i),
fk(t− tk)−Kjk(t− t′k) +Kjk(t− t′k)〉 = (E)ik − (EK)ik − (EK)ki + (P )ik (12)
Combining 10, 11 and 12 we get,
||X(t)−Xhyp(t)||22 = αTFα− αTFKP−1FTKα = αTEα− αTEKα− αTETKα+ αTPα − αTPα+ αTEKα+ αTETKα− αTEKP−1ETKα = αTEα− αTEKP−1ETKα ≤ αTEα
(Since, P is an SPD matrix, αTEKP−1ETKα > 0) (13)
We seek for a bound for the above expression. For that we observe the following:
|(E)ik| = |〈fi(t− ti)−Kji(t− t ′ i), fk(t− tk)−Kjk(t− t ′ k)〉|
= ||fi(t− ti)−Kji(t− t ′ i)||2||fk(t− tk)−Kjk(t− t ′
k)||2.xik (where xik ∈ [0, 1]. We also note that xik is close to 0 when there is not much overlap in the support of the two components and their corresponding fitting kernels.)
≤ (||(fi(t− ti)−Kji(t− ti)||+
||Kji(t− ti)−Kji(t− t ′
i))||). (||fk(t− tk)−Kjk(t− tk)||
+ ||Kjk(t− tk)−Kjk(t− t ′
k)||).xik =⇒ (E)ik ≤ xik.(δ + C∆)2 (14)
Using Gershgorin circle theorem, the maximum eigen value of E : Λmax(E) ≤ maxi((E)ii + Σk 6=i|(E)ik|)
≤ (δ + C∆)2(xmax + 1) (Using 14) (15) (where xmax ∈ [0, N − 1] is a positive number that depends on the maximum overlap of the supports of the component signals and their fitting kernels.) Similarly, the minimum eigen value of F is: Λmin(F ) = mini((F )ii − Σi 6=k|〈fpi(t− ti), fpk(t− tk)〉|) ≥ 1− η (16) (By assumption Σi6=k| < fpi(t− ti), fpk(t− tk) > | ≤ η ) Combining the results from 13, 15 and 16 we get: ||X(t)−Xhyp(t)||2/||X(t)||2 ≤ αTEα/αTFα ≤ Λmax(E)/Λmin(F ) ≤ (δ + C∆)2(xmax + 1)/(1− η) (17)
Finally using 8 we conclude, ||X(t)−X∗(t)||2/||X(t)||2 ≤ ||X(t)−Xhyp(t)||2/||X(t)||2
≤ (δ + C∆)2(xmax + 1)/(1− η)

Review summaries:
The authors describe a method for representing a continuous signal by a pulse code, in a manner inspired by auditory processing in the brain. The resulting framework is somewhat like matching pursuit except that filters are run a single time in a causal manner to find the spike times (which would be faster than MP), and then a N*N least squares problem is solved (which makes it slower). The authors claim that their method will perfectly reconstruct signals of finite innovation rate, however there appear to be mathematical errors in the proof.