from compressive to adaptive sampling of neural and ecg recordings
FROM COMPRESSIVE TO ADAPTIVE SAMPLING OF NEURAL AND ECG RECORDINGS Alexander Singh Alvarado and Jos´e C. Pr´ıncipe University of Florida, Electrical and Computer Engineering Department Gainesville, Florida, 32611 ABSTRACT The miniaturization required for interfacing with the brain demands new methods of transforming neuron responses (spikes) into digital representations. The sparse nature of neural recordings is evident when represented in a shift invariant basis. Although a compressive sensing (CS) framework may seem suitable in reducing the data rates, we show that the time varying sparsity in the signals makes it difficult to apply. Furthermore, we present an adaptive sampling scheme which takes advantage of the local characteristics of the neural spike trains and electrocardiograms (ECG). In contrast to the global constraints imposed in CS our solution is sensitive to the local time structure of the input. The simplicity in the design of the integrate-and-fire (IF) make it a viable solution in current brain machine interfaces (BMI) and ambulatory cardiac monitoring. Index Terms— Adaptive sampling, brain-machine interface, non-uniform sampling, integrate-and-fire model, ECG 1. INTRODUCTION Information in the brain is encoded in the stereotypical output of neurons known as action potentials (AP) that stand out above a noisy background. The goal of BMIs is to extract meaningful information from these neural recordings, in order to take future actions. The first stage in any BMI is to capture the behavior of the hundreds of neurons and reliably transmit APs to the external world. Cortical neurons spike at rates of 5-100 AP per second, the upper limit when they form neural assemblies that are related to behavioral tasks. The algorithms used in conventional BMIs require that the AP be catalogued by neuron and assume the rest of the signal as featureless. It is evident that the regions of interest in neural recordings are localized in time and amplitude. The only information that must be preserved is the shape of each action potential, in order to guarantee discriminability. The recording of the electrocardiogram (ECG) presents a very similar structure and there is a great need to compress ECG for ambulatory, 24 hour monitoring. The sparse events immediately suggest the use of a CS approach. The authors acknowledge NINDS Grant Number: NS053561.
978-1-4577-0539-7/11/$26.00 ©2011 IEEE
633
Although the approach seems reasonable, a number of difficulties arise. First, the CS framework assumes the sparsity of the input is known and constant. This does not hold for neural recordings, since the firing patterns may change drastically. Specially when neurons are involved in neural assemblies. Furthermore, measurements in CS are created based on projections on the entire sparse feature. Although this may reduce the data rates in comparison to conventional analog to digital converters (ADC), the time structure of the signal is being ignored. In contrast we propose an adaptive sampling scheme that can be tuned to the specific features of the input. Instead of imposing global sparseness constraints, we exploit the local structure of the signal. The adaptive sampling scheme proposed is based on the integrate-and-fire model [1]. The sampler is asynchronous and will generate samples based on level crossings of the integral of the signal. Therefore, the resulting sample sets are nonuniformly spaced and clustered in the vicinity of the action potentials, given their prominent amplitude. Recently a variety of adaptive sampling schemes have been implemented, mostly in wireless sensor networks [2]. Other approaches also include information about the specific task in the design of the sampler [3], these are related to active learning in the machine learning context. In this case, the sample locations and values are determined based on the final goal and the available information. Given the abstract nature of the problem, it is difficult to describe a general framework for all adaptive sampling schemes. In this paper we focus on the IF sampler which is a special case of a large class of asynchronous Σ − Δ modulators [4]. These samplers are typically designed such that they oscillate in an autonomous mode. The oscillation frequency is analogous to the clock period in their discrete counterparts. These modulators encode the input into the output sample rate, which is typically hundreds of times the Nyquist rate. Nevertheless, they are attractive since coarse quantization of the amplitude still provides an accurate representation. The input can be recovered using a simple low pass filter. In contrast, the IF sampler is used to sub-sample the input. When used to encode neural recordings or ECG we can still accurately recover the APs and QRS complexes at sub-Nyquist rates. Nevertheless, the recovery algorithms are nonlinear and the compression depends on the signal characteristics.
ICASSP 2011
The paper is organized as follows. We describe both CS and IF sampling schemes and their corresponding recovery algorithms. Both approaches are then compared in terms of the data rates and the recovery accuracy when applied to neural recordings and ECG. Since BMI algorithms are designed on conventional digital systems, we must also consider the effects of quantization under both schemes. In the case of CS sample amplitudes are quantized assuming a constant sampling period while the IF events are quantized in time and consider a coarse amplitude range.
Fig. 1. Bi-phasic integrate-and-fire block diagram. 3. THE BI-PHASIC INTEGRATE-AND-FIRE SAMPLER
2. A CS APPROACH FOR ANALOG SIGNALS Compressive sensing literature [5], was initially framed in terms of discrete finite length signals. The basic framework exploits the fact that given an appropriate representation most signals can be considered sparse. Therefore a relatively small number of measurements is necessary to ensure recovery. It is crucial to define the necessary and sufficient conditions under which a feasible recovery is possible, i.e. constraints on the projection matrix must be imposed. Most authors refer to the sufficient conditions as the coherence measure or the Uniform Uncertainty Principle (UUP) otherwise known as the restricted isometry property (RIP). In recent years, compressive sensing has been extended to the analog signal class providing reconstruction algorithms and hardware implementations. We have identified at least three different trends in the implementation of compressive sampling algorithms for streaming data. The first, consists of using finite length random filters [6]. The measurements are the convolution of the finite response filter and the input. The process can be translated into a linear system of equations y = Φx where y are the measurements, x is the sparse input and Φ is a banded projection matrix for which each row has no more than b nonzero elements, where b is the length of the filter. In this paper we use a similar implementation assuming we have access to the entire sparse signal. The integrate and fire, although a nonlinear system can also be framed in terms of linear constraints but the matrix rows are built from ones with a varying number of nonzero elements which depend on the distance between consecutive samples. On the other hand the measurement vector would only consist of ±θ. Other approaches, such as the random analog demodulator proposed in [7, 8] mix the sparse components directly in the continuous time domain. These are then uniformly sampled presumably at sub-Nyquist rates. A third class of implementations combines traditional frequency based analysis and CS proposed in [9]. Other alternatives that assume sparsity constraints on the input include the work on finite rate of innovation [10]. The reconstruction algorithms separate the problem in two stages the first uses Prony’s method (a.k.a annihilation filter) to determine the location of the sparse components. The second determines the sample values at these locations by least squares.
634
Since the inter-spike activity is mostly noise, an appropriate sampler should exploit the sparse nature of the APs. A similar case occurs in ECG, where the most relevant feature is the QRS complex. The IF model takes advantage of these features, Figure 1 shows the basic components. The continuous input x(t) is convolved with an “averaging function” uk (t), the result is compared against two fixed thresholds, when either of these is reached, a pulse is created at time tk representing the threshold value. The integrator is reset and held at this state depending on the refractory period specified by τ ; the process then repeats. Note that samples are bi-polar, which is one of the differences when compared to other approaches using these type of models as sampling devices [11]. The output pulse train is defined over a discrete set of nonuniformly spaced time stamps. Here we assume uk (t) is constant between two samples. Other averaging functions have been considered in [12]. Therefore, the IF samples satisfy: tk+1 θk = x(t)dt (1) tk +τ
The input is modeled in a shift invariant spaces (SIS), a natural representation for the action potentials. These spaces are well known in approximation theory, nonuniform sampling and compressive sensing. In this case, the input is defined in a M dimensional space: x(t) =
M
ck φ(t − T k)
(2)
k=1
Where φ(t) are known as the generators of the space and T is the underlying period. Which is determined in relation to the Nyquist period. The generator is chosen to be a cubic B(asis)spline defined over a uniform knot sequence. The samples generated by the IF provide a system of equations from which we can determine the sparse coefficients, θ = Sc, where the elements in S are defined by: ti+1 Si,k = φ(t − T k)dt (3) ti +τ
The vector θ contains the threshold values and c are the coefficients. Note that the size and elements of the matrix are signal
dependent and since the basis are compactly supported most of the non zero elements are near the diagonal. In this case we can ensure perfect recovery as long as, one sample falls in the support of each basis. Nevertheless, we are not concerned with full recovery of the signal, but rather the spike regions [12]. We have also presented recovery algorithms in the typical bandlimited space, along with the corresponding error bounds [13].
but then grows slowly, while CS soars to nearly 100 dB. The low SNR reconstruction in the IF is due to the low amplitude details in the AP as seen in Figure 2(a), in order to reproduce these small variations we would need to select the threshold extremely low. In contrast, CS only assumes a sparse signal and does not constrain the amplitude of these elements, therefore the conditions are exactly satisfied. Nevertheless, these small details are in the range of the noise for this application and are not crucial. This trend does not extend to real
4. RESULTS
−5
−5
x 10
x 10 6
4 Amplitude [Volts]
4
Amplitude [Volts]
6
original rec−CS rec−IF
2
0
2 0 −2 −4 −6
−2
original rec−CS rec−IF
−8
−4
0.01
0.011
0.012
0.013 0.014 Time [s]
0.015
0.016
(a) Neural simulator.
Threshold
−9
−6
5
x 10
0 −5
0.02
0.022
0.024 Time [s]
0.026
0.028
(b) Real data
Fig. 2. Data recovered using CS and IF. The underline defines the AP region. require knowledge of the precise timings of the APs for each neuron. Therefore, we are only concerned with the reconstruction accuracy in the vicinity of the AP. The recovery error is measured in terms of the signal to error ratio (SER) (10 log(power signal/power error)) within a window around the AP, where the error is simply given by the difference between both signals. Figure 3(a) shows the average variation of the SER around the three AP regions for both approaches in relation to the number of samples. In the case of the IF the samples were increased by reducing the threshold, keeping the refractory period at zero. Note that the reconstruction performance of the IF outperforms the CS for very few samples
635
120
Average spike SER [dB]
80 60 40 20 0 0
50
IF CS
100 Average spike SER [dB]
We show results both on neural and ECG recordings. The neural data used has been recorded from a behaving animal. This data was provided by the Neuro-prosthetics Research Group at the University of Florida. We also use synthetic data in order to satisfy the sparsity constraints imposed by CS. The synthetic data was created by placing synthetic APs and forcing the inter-spike regions to zero. Figure 2(a) shows the reconstruction for this type of signals using the IF and CS approaches. The insert shows the complete signal consisting of three APs. Since the IF assumes an analog input, the precise timings for the samples are determined based on an oversampled version of the input. The measurements in the CS case, were determined by projections of the input with a random matrix whose elements are instances of a zero mean unit variance gaussian random variable. We use the optimization methods proposed in [5], based on an interior point method to determine the sparse vector. As mentioned earlier BMIs
100
200
300 400 500 Number of samples
600
(a) Neural simulator.
700
IF CS
40
30
20 Nyquist Rate 10
0 0
100
200
300 400 500 Number of samples
600
700
800
(b) Real data.
Fig. 3. Variation in the reconstruction accuracy in relation to the number of samples.
recordings, as seen in Figure 2(b). Here we show a zoomed in version of the last AP, these have amplitudes near 50μV and a bandwidth of 5KHz. Since they are no longer strictly sparse the performance using CS deteriorates as seen in Figure 3(b). The IF outperforms the conventional CS approach, and provides a feasible representation even at sub-Nyquist rates. In this case the Nyquist rate is 10Khz and the signal duration is 37ms which translates into 368 samples. We have also extended our results to ECG recordings, which share similar time structure with neural signals. Here we used the QT database available at PhysioBank [14]. The recordings have been sampled at 250Hz. Figure 4 shows the original and recovered signals, as well as the variation of the reconstruction on the entire window in relation to the number of samples. The entire signal has a duration of approximately 4s therefore the number of samples equivalent to the Nyquist rate would be 1024. These results are comparable to those in [15] based on a CS approach using wavelet decompositions. Although the SER seems to be low, sufficient information of the signal is preserved to make reliable estimates on cardiac features. Since the data must eventually be processed by a digital system we also show the variation in the recovery based on the quantization. In CS the amplitudes of the measurements are quantized, while in the IF the sample time instants are quantized to a specific clock as seen in Figure 4. These results were obtained using the synthetic data. In the case of the IF we require at least 1μs resolution at the receiver side, this does not imply we need to transmit the quantized time. In contrast, each measurement transmitted using the CS approach requires at least 10 bits.
6. REFERENCES
35
Amplitude [mV]
1
Average spike SER [dB]
Original IF
1.2
0.8 0.6 0.4 0.2 0 −0.2
Threshold
−0.4 2
−3
2
x 10
2.2
2.4
2.6
2.8
[2] R.M. Willett, A.M. Martin, and R.D. Nowak, “Adaptive sampling for wireless sensor networks,” Information Theory, 2004. ISIT 2004. Proceedings. International Symposium on, p. 519, jun. 2004.
25
20
15 0
3
[1] Wulfram Gerstner and Werner Kistler, “Spiking neuron models: An introduction,” Cambridge University Press, 2002.
Nyquist Rate
30
500
1000
1500
Number of samples
[3] Rui Castro, Active Learning and Adaptive Sampling for NonParametric Inference, Ph.D. thesis, Rice University, aug 2007.
2000
[4] S. Ouzounov, Engel Roza, J.A. Hegt, G. van der Weide, and A.H.M. van Roermund, “Analysis and design of high-performance asynchronous sigma-delta modulators with a binary quantizer,” Solid-State Circuits, IEEE Journal of, vol. 41, no. 3, pp. 588 – 596, mar. 2006.
0 −2
2
2.2
2.4
2.6
Time [s]
2.8
1
3
2
Time [s]
3
4
Fig. 4. ECG reconstruction and performance at an average sample rate of 50 samp/s. 90
35
Average spike SER [dB]
Average spike SER [dB]
80 30
25
20
15
70
−8
10
−7
10 Clock period [s]
−6
10
−5
10
50 40 30
10 0
5
10 15 20 25 Number of quantization bits
30
[6] J.A. Tropp, M.B. Wakin, M.F. Duarte, D. Baron, and R.G. Baraniuk, “Random filters for compressive sampling and reconstruction,” Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on, vol. 3, pp. III –III, may. 2006. [7] J.N. Laska, S. Kirolos, M.F. Duarte, T.S. Ragheb, R.G. Baraniuk, and Y. Massoud, “Theory and implementation of an analog-to-information converter using random demodulation,” Circuits and Systems, 2007. ISCAS 2007. IEEE International Symposium on, pp. 1959 –1962, may. 2007.
60
20 10 −9 10
[5] E.J. Candes and T. Tao, “Decoding by linear programming,” Information Theory, IEEE Transactions on, vol. 51, no. 12, pp. 4203–4215, dec. 2005.
35
(a) Time quantization of IF samples. (b) Amplitude quantization of CS samples.
Fig. 5. Effect of quantization on recovery methods for both approaches.
[8] J.A. Tropp, J.N. Laska, M.F. Duarte, J.K. Romberg, and R.G. Baraniuk, “Beyond nyquist: Efficient sampling of sparse bandlimited signals,” Information Theory, IEEE Transactions on, vol. 56, no. 1, pp. 520 – 544, jan. 2010. [9] M. Mishali and Y.C. Eldar, “Xampling: Analog data compression,” Data Compression Conference (DCC), 2010, pp. 366 –375, mar. 2010. [10] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” Signal Processing, IEEE Transactions on, vol. 50, no. 6, pp. 1417 –1428, jun. 2002.
5. CONCLUSIONS We have shown that in an ideal noiseless scenario CS provides a feasible solution to neural encoding. Nevertheless, when applied to real recordings the lack of sparsity in the time domain restricts the use of this type of CS method. Furthermore, when dealing with bursting neurons, or signals with varying sparsity modifications are required in order to apply CS. In contrast, the IF can be readily used in streaming data applications, since it is based on local structure. The IF provides compression, while accurately representing the spike regions and the QRS complex, even at sub-Nyquist rates. However, if high accuracy is needed the IF data rates will increase substantially. Since BMI algorithms are based on conventional digital systems quantization is inevitable. We evaluate the effect of amplitude quantization of the measurements in the case of CS as well as time quantization in the IF. Although we presented both approaches independently, it may also be possible to extend the IF encoding scheme in order to take advantage of the optimization techniques developed for CS, this would involve input dependent projection matrices.
636
[11] Aurel A. Lazar, “Time encoding with an integrate-and-fire neuron with a refractory period,” Neurocomputing, vol. 58-60, pp. 53–58, jun 2004. [12] A.S. Alvarado, J.C. Principe, and J.G. Harris, “Stimulus reconstruction from the biphasic integrate-and-fire sampler,” Neural Engineering, 2009. NER ’09. 4th International IEEE/EMBS Conference on, pp. 415 –418, apr. 2009. [13] H.G. Feichtinger, J. C. Principe, J.L. Romero, A. S. Alvarado, and G. Velasco, “Approximate reconstruction of bandlimited functions for the integrate and fire sampler,” Advances in computational mathematics, 2010, [Accepted]. [14] A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. Ch. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley, “PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals,” Circulation, vol. 101, no. 23, pp. e215–e220, 2000 (June 13), Circulation Electronic Pages: http://circ.ahajournals.org/cgi/content/full/101/23/e215. [15] Octavian Adrian Postolache Eduardo Correia Pinheiro and Pedro Silva Girao, “Implementation of compressed sensing in telecardiology sensor networks,” International Journal of Telemedicine and Applications, jul. 2010.
978-1-4577-0539-7/11/$26.00 ©2011 IEEE
633
Although the approach seems reasonable, a number of difficulties arise. First, the CS framework assumes the sparsity of the input is known and constant. This does not hold for neural recordings, since the firing patterns may change drastically. Specially when neurons are involved in neural assemblies. Furthermore, measurements in CS are created based on projections on the entire sparse feature. Although this may reduce the data rates in comparison to conventional analog to digital converters (ADC), the time structure of the signal is being ignored. In contrast we propose an adaptive sampling scheme that can be tuned to the specific features of the input. Instead of imposing global sparseness constraints, we exploit the local structure of the signal. The adaptive sampling scheme proposed is based on the integrate-and-fire model [1]. The sampler is asynchronous and will generate samples based on level crossings of the integral of the signal. Therefore, the resulting sample sets are nonuniformly spaced and clustered in the vicinity of the action potentials, given their prominent amplitude. Recently a variety of adaptive sampling schemes have been implemented, mostly in wireless sensor networks [2]. Other approaches also include information about the specific task in the design of the sampler [3], these are related to active learning in the machine learning context. In this case, the sample locations and values are determined based on the final goal and the available information. Given the abstract nature of the problem, it is difficult to describe a general framework for all adaptive sampling schemes. In this paper we focus on the IF sampler which is a special case of a large class of asynchronous Σ − Δ modulators [4]. These samplers are typically designed such that they oscillate in an autonomous mode. The oscillation frequency is analogous to the clock period in their discrete counterparts. These modulators encode the input into the output sample rate, which is typically hundreds of times the Nyquist rate. Nevertheless, they are attractive since coarse quantization of the amplitude still provides an accurate representation. The input can be recovered using a simple low pass filter. In contrast, the IF sampler is used to sub-sample the input. When used to encode neural recordings or ECG we can still accurately recover the APs and QRS complexes at sub-Nyquist rates. Nevertheless, the recovery algorithms are nonlinear and the compression depends on the signal characteristics.
ICASSP 2011
The paper is organized as follows. We describe both CS and IF sampling schemes and their corresponding recovery algorithms. Both approaches are then compared in terms of the data rates and the recovery accuracy when applied to neural recordings and ECG. Since BMI algorithms are designed on conventional digital systems, we must also consider the effects of quantization under both schemes. In the case of CS sample amplitudes are quantized assuming a constant sampling period while the IF events are quantized in time and consider a coarse amplitude range.
Fig. 1. Bi-phasic integrate-and-fire block diagram. 3. THE BI-PHASIC INTEGRATE-AND-FIRE SAMPLER
2. A CS APPROACH FOR ANALOG SIGNALS Compressive sensing literature [5], was initially framed in terms of discrete finite length signals. The basic framework exploits the fact that given an appropriate representation most signals can be considered sparse. Therefore a relatively small number of measurements is necessary to ensure recovery. It is crucial to define the necessary and sufficient conditions under which a feasible recovery is possible, i.e. constraints on the projection matrix must be imposed. Most authors refer to the sufficient conditions as the coherence measure or the Uniform Uncertainty Principle (UUP) otherwise known as the restricted isometry property (RIP). In recent years, compressive sensing has been extended to the analog signal class providing reconstruction algorithms and hardware implementations. We have identified at least three different trends in the implementation of compressive sampling algorithms for streaming data. The first, consists of using finite length random filters [6]. The measurements are the convolution of the finite response filter and the input. The process can be translated into a linear system of equations y = Φx where y are the measurements, x is the sparse input and Φ is a banded projection matrix for which each row has no more than b nonzero elements, where b is the length of the filter. In this paper we use a similar implementation assuming we have access to the entire sparse signal. The integrate and fire, although a nonlinear system can also be framed in terms of linear constraints but the matrix rows are built from ones with a varying number of nonzero elements which depend on the distance between consecutive samples. On the other hand the measurement vector would only consist of ±θ. Other approaches, such as the random analog demodulator proposed in [7, 8] mix the sparse components directly in the continuous time domain. These are then uniformly sampled presumably at sub-Nyquist rates. A third class of implementations combines traditional frequency based analysis and CS proposed in [9]. Other alternatives that assume sparsity constraints on the input include the work on finite rate of innovation [10]. The reconstruction algorithms separate the problem in two stages the first uses Prony’s method (a.k.a annihilation filter) to determine the location of the sparse components. The second determines the sample values at these locations by least squares.
634
Since the inter-spike activity is mostly noise, an appropriate sampler should exploit the sparse nature of the APs. A similar case occurs in ECG, where the most relevant feature is the QRS complex. The IF model takes advantage of these features, Figure 1 shows the basic components. The continuous input x(t) is convolved with an “averaging function” uk (t), the result is compared against two fixed thresholds, when either of these is reached, a pulse is created at time tk representing the threshold value. The integrator is reset and held at this state depending on the refractory period specified by τ ; the process then repeats. Note that samples are bi-polar, which is one of the differences when compared to other approaches using these type of models as sampling devices [11]. The output pulse train is defined over a discrete set of nonuniformly spaced time stamps. Here we assume uk (t) is constant between two samples. Other averaging functions have been considered in [12]. Therefore, the IF samples satisfy: tk+1 θk = x(t)dt (1) tk +τ
The input is modeled in a shift invariant spaces (SIS), a natural representation for the action potentials. These spaces are well known in approximation theory, nonuniform sampling and compressive sensing. In this case, the input is defined in a M dimensional space: x(t) =
M
ck φ(t − T k)
(2)
k=1
Where φ(t) are known as the generators of the space and T is the underlying period. Which is determined in relation to the Nyquist period. The generator is chosen to be a cubic B(asis)spline defined over a uniform knot sequence. The samples generated by the IF provide a system of equations from which we can determine the sparse coefficients, θ = Sc, where the elements in S are defined by: ti+1 Si,k = φ(t − T k)dt (3) ti +τ
The vector θ contains the threshold values and c are the coefficients. Note that the size and elements of the matrix are signal
dependent and since the basis are compactly supported most of the non zero elements are near the diagonal. In this case we can ensure perfect recovery as long as, one sample falls in the support of each basis. Nevertheless, we are not concerned with full recovery of the signal, but rather the spike regions [12]. We have also presented recovery algorithms in the typical bandlimited space, along with the corresponding error bounds [13].
but then grows slowly, while CS soars to nearly 100 dB. The low SNR reconstruction in the IF is due to the low amplitude details in the AP as seen in Figure 2(a), in order to reproduce these small variations we would need to select the threshold extremely low. In contrast, CS only assumes a sparse signal and does not constrain the amplitude of these elements, therefore the conditions are exactly satisfied. Nevertheless, these small details are in the range of the noise for this application and are not crucial. This trend does not extend to real
4. RESULTS
−5
−5
x 10
x 10 6
4 Amplitude [Volts]
4
Amplitude [Volts]
6
original rec−CS rec−IF
2
0
2 0 −2 −4 −6
−2
original rec−CS rec−IF
−8
−4
0.01
0.011
0.012
0.013 0.014 Time [s]
0.015
0.016
(a) Neural simulator.
Threshold
−9
−6
5
x 10
0 −5
0.02
0.022
0.024 Time [s]
0.026
0.028
(b) Real data
Fig. 2. Data recovered using CS and IF. The underline defines the AP region. require knowledge of the precise timings of the APs for each neuron. Therefore, we are only concerned with the reconstruction accuracy in the vicinity of the AP. The recovery error is measured in terms of the signal to error ratio (SER) (10 log(power signal/power error)) within a window around the AP, where the error is simply given by the difference between both signals. Figure 3(a) shows the average variation of the SER around the three AP regions for both approaches in relation to the number of samples. In the case of the IF the samples were increased by reducing the threshold, keeping the refractory period at zero. Note that the reconstruction performance of the IF outperforms the CS for very few samples
635
120
Average spike SER [dB]
80 60 40 20 0 0
50
IF CS
100 Average spike SER [dB]
We show results both on neural and ECG recordings. The neural data used has been recorded from a behaving animal. This data was provided by the Neuro-prosthetics Research Group at the University of Florida. We also use synthetic data in order to satisfy the sparsity constraints imposed by CS. The synthetic data was created by placing synthetic APs and forcing the inter-spike regions to zero. Figure 2(a) shows the reconstruction for this type of signals using the IF and CS approaches. The insert shows the complete signal consisting of three APs. Since the IF assumes an analog input, the precise timings for the samples are determined based on an oversampled version of the input. The measurements in the CS case, were determined by projections of the input with a random matrix whose elements are instances of a zero mean unit variance gaussian random variable. We use the optimization methods proposed in [5], based on an interior point method to determine the sparse vector. As mentioned earlier BMIs
100
200
300 400 500 Number of samples
600
(a) Neural simulator.
700
IF CS
40
30
20 Nyquist Rate 10
0 0
100
200
300 400 500 Number of samples
600
700
800
(b) Real data.
Fig. 3. Variation in the reconstruction accuracy in relation to the number of samples.
recordings, as seen in Figure 2(b). Here we show a zoomed in version of the last AP, these have amplitudes near 50μV and a bandwidth of 5KHz. Since they are no longer strictly sparse the performance using CS deteriorates as seen in Figure 3(b). The IF outperforms the conventional CS approach, and provides a feasible representation even at sub-Nyquist rates. In this case the Nyquist rate is 10Khz and the signal duration is 37ms which translates into 368 samples. We have also extended our results to ECG recordings, which share similar time structure with neural signals. Here we used the QT database available at PhysioBank [14]. The recordings have been sampled at 250Hz. Figure 4 shows the original and recovered signals, as well as the variation of the reconstruction on the entire window in relation to the number of samples. The entire signal has a duration of approximately 4s therefore the number of samples equivalent to the Nyquist rate would be 1024. These results are comparable to those in [15] based on a CS approach using wavelet decompositions. Although the SER seems to be low, sufficient information of the signal is preserved to make reliable estimates on cardiac features. Since the data must eventually be processed by a digital system we also show the variation in the recovery based on the quantization. In CS the amplitudes of the measurements are quantized, while in the IF the sample time instants are quantized to a specific clock as seen in Figure 4. These results were obtained using the synthetic data. In the case of the IF we require at least 1μs resolution at the receiver side, this does not imply we need to transmit the quantized time. In contrast, each measurement transmitted using the CS approach requires at least 10 bits.
6. REFERENCES
35
Amplitude [mV]
1
Average spike SER [dB]
Original IF
1.2
0.8 0.6 0.4 0.2 0 −0.2
Threshold
−0.4 2
−3
2
x 10
2.2
2.4
2.6
2.8
[2] R.M. Willett, A.M. Martin, and R.D. Nowak, “Adaptive sampling for wireless sensor networks,” Information Theory, 2004. ISIT 2004. Proceedings. International Symposium on, p. 519, jun. 2004.
25
20
15 0
3
[1] Wulfram Gerstner and Werner Kistler, “Spiking neuron models: An introduction,” Cambridge University Press, 2002.
Nyquist Rate
30
500
1000
1500
Number of samples
[3] Rui Castro, Active Learning and Adaptive Sampling for NonParametric Inference, Ph.D. thesis, Rice University, aug 2007.
2000
[4] S. Ouzounov, Engel Roza, J.A. Hegt, G. van der Weide, and A.H.M. van Roermund, “Analysis and design of high-performance asynchronous sigma-delta modulators with a binary quantizer,” Solid-State Circuits, IEEE Journal of, vol. 41, no. 3, pp. 588 – 596, mar. 2006.
0 −2
2
2.2
2.4
2.6
Time [s]
2.8
1
3
2
Time [s]
3
4
Fig. 4. ECG reconstruction and performance at an average sample rate of 50 samp/s. 90
35
Average spike SER [dB]
Average spike SER [dB]
80 30
25
20
15
70
−8
10
−7
10 Clock period [s]
−6
10
−5
10
50 40 30
10 0
5
10 15 20 25 Number of quantization bits
30
[6] J.A. Tropp, M.B. Wakin, M.F. Duarte, D. Baron, and R.G. Baraniuk, “Random filters for compressive sampling and reconstruction,” Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on, vol. 3, pp. III –III, may. 2006. [7] J.N. Laska, S. Kirolos, M.F. Duarte, T.S. Ragheb, R.G. Baraniuk, and Y. Massoud, “Theory and implementation of an analog-to-information converter using random demodulation,” Circuits and Systems, 2007. ISCAS 2007. IEEE International Symposium on, pp. 1959 –1962, may. 2007.
60
20 10 −9 10
[5] E.J. Candes and T. Tao, “Decoding by linear programming,” Information Theory, IEEE Transactions on, vol. 51, no. 12, pp. 4203–4215, dec. 2005.
35
(a) Time quantization of IF samples. (b) Amplitude quantization of CS samples.
Fig. 5. Effect of quantization on recovery methods for both approaches.
[8] J.A. Tropp, J.N. Laska, M.F. Duarte, J.K. Romberg, and R.G. Baraniuk, “Beyond nyquist: Efficient sampling of sparse bandlimited signals,” Information Theory, IEEE Transactions on, vol. 56, no. 1, pp. 520 – 544, jan. 2010. [9] M. Mishali and Y.C. Eldar, “Xampling: Analog data compression,” Data Compression Conference (DCC), 2010, pp. 366 –375, mar. 2010. [10] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” Signal Processing, IEEE Transactions on, vol. 50, no. 6, pp. 1417 –1428, jun. 2002.
5. CONCLUSIONS We have shown that in an ideal noiseless scenario CS provides a feasible solution to neural encoding. Nevertheless, when applied to real recordings the lack of sparsity in the time domain restricts the use of this type of CS method. Furthermore, when dealing with bursting neurons, or signals with varying sparsity modifications are required in order to apply CS. In contrast, the IF can be readily used in streaming data applications, since it is based on local structure. The IF provides compression, while accurately representing the spike regions and the QRS complex, even at sub-Nyquist rates. However, if high accuracy is needed the IF data rates will increase substantially. Since BMI algorithms are based on conventional digital systems quantization is inevitable. We evaluate the effect of amplitude quantization of the measurements in the case of CS as well as time quantization in the IF. Although we presented both approaches independently, it may also be possible to extend the IF encoding scheme in order to take advantage of the optimization techniques developed for CS, this would involve input dependent projection matrices.
636
[11] Aurel A. Lazar, “Time encoding with an integrate-and-fire neuron with a refractory period,” Neurocomputing, vol. 58-60, pp. 53–58, jun 2004. [12] A.S. Alvarado, J.C. Principe, and J.G. Harris, “Stimulus reconstruction from the biphasic integrate-and-fire sampler,” Neural Engineering, 2009. NER ’09. 4th International IEEE/EMBS Conference on, pp. 415 –418, apr. 2009. [13] H.G. Feichtinger, J. C. Principe, J.L. Romero, A. S. Alvarado, and G. Velasco, “Approximate reconstruction of bandlimited functions for the integrate and fire sampler,” Advances in computational mathematics, 2010, [Accepted]. [14] A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. Ch. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley, “PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals,” Circulation, vol. 101, no. 23, pp. e215–e220, 2000 (June 13), Circulation Electronic Pages: http://circ.ahajournals.org/cgi/content/full/101/23/e215. [15] Octavian Adrian Postolache Eduardo Correia Pinheiro and Pedro Silva Girao, “Implementation of compressed sensing in telecardiology sensor networks,” International Journal of Telemedicine and Applications, jul. 2010.
