1 Introduction
Heartrate monitoring is ubiquitous in modern wearable devices such as a smart watch (phan2015smartwatch, ; aarts2006apparatus, ) or Electrocardiogram (ECG) necklace (penders2011low, ). ECG sensors (mukhopadhyay2015wearable, ; gyselinckx2005human, ) attached to these devices monitor the electrical potential caused by the systolic activity of heart and then propagated through cardiac muscles. The recorded electrical data are postprocessed, either locally on the sensor (van2015345, ) or on a device (7367297, ) attached to the sensor to estimate heartrate. QRS pattern identification (see Figure 1) from ECG is fundamental to heartrate estimation. Although QRS detection has achieved significant maturity over time (kohler2002principles, ), recent advances in wearable healthcare (gyselinckx2005human, ; otto2006system, ) have motivated researchers to revisit QRS. This is due to

ECG readings from wearable sensors are contaminated with motion artifacts and baseline drifts; and

devices integrating wearable sensors are constrained in terms of area, power consumption and computational capabilities.
Several approaches have been proposed recently to detect QRS complex from wearable ECG using software techniques. A level crossing approach is proposed in (ravanshad2014level, ), where a modified levelcrossing analogtodigital converter is introduced to convert analog ECG data to a meaningful representation. Based on this representation, an algorithm is proposed to detect the RR intervals. A timefrequency representation (called the STransform) is introduced in (zidelmal2014qrs, ) to isolate QRS complexes in timefrequency domain. Shanon energy is computed on these isolated spectrums to localize Rwaves in time domain. A realtime signal processing approach is proposed in (Karimipour2014153, ), which includes high frequency noise filtering and baseline drift reduction using discrete wavelet transform. In (RKSPL2014, ), a threshold independent approach is proposed based on first derivative of ECG signals. There are also other deterministic approaches for QRS detection. For a summary of these approaches, readers are referred to (Jain2017362, ). Some of the recent works have also addressed low power hardware implementations (van2015345, ; deepu2016ecg, ; ieong2016sub, ; kim2017energy, ). Softwarebased QRS detection techniques require general purpose computing platforms (such as a CPU core). Power consumption of these techniques is usually of the order of W. One advantage of these approaches is the ease of updating an existing algorithm or implementing a new usecase (such as arrhythmia detection), with software updates only (no hardware change is required). Hardwarebased QRS detection techniques require dedicated hardware to implement the chosen algorithm. Although, sub
W power consumption is achieved in hardwarebased techniques, limited flexibility is offered by these designs to implement a new algorithm or a new usecase. For these techniques, the usual approach is to do feature extraction or QRS detection on the sensor with usecases such as arrhythmia detection implemented on the device
(tekeste2017nano, ).In recent years, artificial intelligent approaches have been investigated for QRS detection to address flexibility, while consuming lower power, an important consideration for powerconstrained wearable devices. Artificial neural network has been used to carry out the classification task in
(silpioTSP98, ). Such a network forms the basis for ECG analysis including arrhythmia, myocardial ischemia and other chronic alterations. Knearest neighbor classifier is used for QRS complex detection from ECG in
(Saini2013331, ). The ECG signal is postprocessed using a digital bandpass filter to reduce false detection with gradient of the signal used as a feature for the classifier. Radial basis function is used for QRS complex identification in
(Arbateni2014438, ). Similar to the previous approach, this technique uses filters as a preprocessing step with a centering approach for adjusting the Rpeak positions. Support vector machine is used in
magrans2016myocardialto classify QRS segments. An 11 layer deep convolution neural network based approach is proposed in
(acharya2017automated, ; acharya2017automated2, )to classify ECG segments. All these approaches use the classical neural network approaches with supervised learning. Success of these approaches depends on availability of large amount of hand labeled data, generic enough to be applied to a broad range of subjects with and without cardiac irregularities.
In this work we use spiking neural networks (maass1997networks, )
for heartrate estimation. These are powerful and biologically realistic computation models, inspired by the dynamics of human brain. Spiking neural networks are gaining popularity in solving complex pattern recognition
(buonomano1999neural, ; kasabov2013dynamic, ), function approximation (iannella2001spiking, ) and image classification (diehl2015fast, ; diehl2015unsupervised, ; samadi2017deep, ) tasks. Another reason for widespread success of spiking neural networks is their efficient VLSI implementation (indiveri2006vlsi, ; fusi1998collective, ; hsieh2012vlsi, ). Examples are the large scale neuromorphic computing^{1}^{1}1The term neuromorphic computing was first coined in (mead1990neuromorphic, ). systems such as TrueNorth (akopyan2015truenorth, ), CxQuad (indiveri2015neuromorphic, ), NeuCube (kasabov2014neucube, ), SpiNNaker (khan2008spinnaker, ), NeuroGrid (benjamin2014neurogrid, ) and HICANN (schemmel2010wafer, ), among others. Some of these systems are originally designed for highperformance computing (e.g., TrueNorth and SpiNNaker), while others are designed for lowpower embedded systems (e.g., CxQuad and HICANN). Readers are referred to schuman2017surveyfor a recent survey of neuromorphic hardware. In this work we report energy usage on CxQuad, a low power spiking hardware with 1024 neurons and 64K synapses (
(indiveri2015neuromorphic, )).Our work differentiates from existing studies on ECGbased heartrate estimation by (1) using spiking neural networks, which can be implemented on energy efficient neuromorphic hardware; (2) encoding analog ECG signal directly into spike trains, which are then used to excite the network of spiking neurons; and (3) designing an unsupervised readout, facilitating learning from subjectspecific ECG to estimate heartrate. The overall approach allows personalization and eliminates the need to manually annotate training data. We envision our approach to be integrated inside an ECG sensor node. Analog ECG signal is encoded directly into spikes, which are used to excite a reservoir of recurrently connected spiking neurons. This is inspired by the computation model of Liquid State Machine (maass2002real, ). Neurons in the architecture are interconnected using plastic synapses, with weight updates using Spike Timing Dependent Plasticity (STDP) (rao2001spike, ; brader2007learning, ), an important form of Hebbian Learning. Additionally, homeostatic synaptic scaling (liu2011learning, ; galtier2013biological, ; carlson2013biologically, ) is used to stabilize the plastic mechanism, preventing runaway behaviors. At the readout stage of the architecture, we use particle swarm optimization (eberhart1995new, ) to select contributions from a subset (called the winning set) of the spiking neurons. Cumulative responses from these winning neurons are clustered to infer heartrate using Fuzzy cMeans clustering (bezdek1984fcm, ). To validate our approach we used CARLsim (beyeler2015carlsim, ) spiking neural network simulator with ECG data from inhouse clinical trials and also opensource databases. We compared our approach with (1) deterministic QRS detection technique of (van2015345, ); (2) the neural network based technique of (acharya2017automated, ); and the support vector machine based technique of (magrans2016myocardial, )^{2}^{2}2Original flavors of (acharya2017automated, ; magrans2016myocardial, ) are adapted for QRS detection and heartrate estimation only.. Results demonstrate high accuracy of our approach, with applicability to subjects with and without cardiac irregularities.
Contributions: Following are our novel contributions:

a technique to encode timeseries data to spike train, capturing its spatiotemporal characteristics;

a novel architecture inspired by the computation model of Liquid State Machine;

a novel learningrule with soft winnertakeall implementation, facilitating temporal pattern detection;

an unsupervised readout for heartrate estimation using Fuzzy cMeans;

a technique to improve clustering accuracy by selecting neurons responses using swarm intelligence; and

a reallife medical benchmark for neuromorphic computing community.
2 Overview of our Approach
Our approach is based on the computational model of Liquid State Machine (LSM) (maass2002real, )
, which consists of a preprocessing unit (called the “Liquid”) with recurrently connected spiking neurons and a readout unit to decode/interpret the internal states of the liquid. LSM can implement any nonlinear transformation of the input. Internally, the Liquid integrates temporal information into linearly independent Liquid states, without computing precisely any nonlinear transformation; the readout has no notion of temporal aspects and learns to map Liquid states to the function approximation/classification. Following are some characteristics of LSM that are relevant for the heartrate estimation problem we aim to solve in this work.

Recurrently connected neurons in the Liquid can implement memory of the input spikes at many (potentially infinitely many) preceding time instances. This is required as morphological and certain other dynamic features of the ECG are a result of cardiac changes over a period of time, which needs to be incorporated for personalized heartrate estimation (temporal).

Separating the computation into Liquid and readout allows to instantiate multiple readouts (implementing a large range of functionalities) using the same Liquid. In this work we discuss one such readout instantiation using unsupervised learning algorithm, directed by the statistics of the spike input to the Liquid. Additionally in the context of wearables, the Liquid can be instantiated inside the sensor with multiple readouts implemented on the device (flexibility).

LSM allows realtime computation with fast and robust learning from streaming data. This makes the computational model specially suitable for heartrate estimation, where onsensor subjectspecific learning is desirable.
An overview of our approach is shown in Figure 2. Spikes generated from the encoder (Spike Encoder) are used to excite the Liquid (Spiking Network). Liquid states are used to infer heartrate in the readout unit (HeartRate Decoder). We envision the spike encoder and the spiking network (Liquid) to be integrated inside the sensor, while the readout resides on the wearable device. In this way, spatiotemporal characteristics of the ECG processed by the Liquid can be used to implement multiple usecases (as separate readouts) on a wearable device, depending on user preference and market needs. Examples include arrhythmia detection or ECGbased authentication.
2.1 Spike Encoder
Information in a spiking neural network can be encoded using two techniques – ratebased coding and temporal coding (van2001rate, ). Ratebased coding encodes information as number of spikes within an encoding window without considering the temporal characteristics of the signal. Ratebased encoding has been successfully applied to spatial classification tasks such as handwritten digit recognition (diehl2015unsupervised, ). Temporal coding (van2001novel, ), on the other hand, encodes information as interspike intervals, capturing the spatiotemopral structure of the input signal. Temporal encoding has been successfully applied for timeseries processing such as speech processing (Tavanaei2017191, ) and EEGbased brainmachine interface (CorradiTBME, ). For heartrate estimation using ECG, temporal characteristics around QRS complexes need to encoded as interspike intervals and therefore, temporal coding is adopted in this work.
2.1.1 Overview of the Spike Encoding Circuit Design
The spike encoder encodes input ECG signal to interspike intervals using a combination of threshold modulators, voltage comparator, spike generator and a timer. This is shown in Figure 3. The threshold modulators are implemented using full adders. Such adders can be efficiently implemented using thin film transistors (bahubalindruni2016novel, ). Reference analog voltage comparator and spike generation implementations are discussed in (li2014ultra, ; du2011bio, ). The timer in Figure 3 can be implemented using D flip flops (DFF), clocked using an external crystal oscillator. The DFF implements clock division to generate the trigger intervals for the comparators. In this work, we use 2 ms trigger intervals, i.e. the spike encoder circuit operates at 500 Hz. Accuracypower consumption tradeoff obtained by varying this trigger interval is left as future work. It is to be noted that the comparator circuit can be shared to reduce power consumption. Figure 4 shows the schematic of such an arrangement for a sensor node (balsamo2016hibernus, ) with a maximum propagation delay of . The two analog switches multiplex the signal propagated to the comparator^{3}^{3}3Detailed circuit analysis of the spike encoder is beyond the scope of this work.. In the next section we describe the high level working of the spike encoding circuit.
2.1.2 Working Principle of the Spike Encoder
The working principle of the spike encoder is described considering three scenarios – ECG signal is rising, ECG signal is falling and ECG signal is stable within a time window Delta. The algorithm for the spike encoder is shown as a flowchart in Figure 5. The encoding is based on two thresholds – Lthr and Uthr (). Voltage comparisons are performed at a fixed intervals controlled using the timer. At every interval, the comparator is enabled to compare the ECG voltage and the threshold Uthr. For simplicity, we assume that the ECG signal is rising. The voltage comparison is positive, triggering the following sequence of events.

threshold update: Lthr = Uthr and Uthr = Uthr + Delta

enable spike generator to emit one spike

restart timer and return to wait
If the result of the comparison is negative (meaning the ECG signal is either falling or stable within Delta), a second comparison is performed, where the ECG voltage is compared with the threshold Lthr. If this comparison is positive (in the case when the ECG signal is falling), the following sequence of events are triggered.

threshold update: Uthr = Lthr and Lthr = Lthr  Delta

restart timer and return to wait
Finally, if the result of the second comparison is also negative (in the case when the ECG signal is stable within Delta), thresholds are not updated and timer is restarted. Following are the specific characteristics of this approach.

Threshold updates are performed to track the ECG signal in the upward or downward directions.

No spikes are generated when the ECG signal is falling. This is a design choice we use in this work as the spatiotemporal characteristics of the QRS can be captured using the rising part of the voltage waveform.

No spikes are generated if the ECG signal is stable.
Figure 6(a) shows the spike generation from a segment of the ECG signal. Also highlighted in this figure are the regions of interest i.e., the QRS peaks (shown in the middle figure). As discussed before and this can be seen quite well from this figure, spike temporal coding captures important spatiotemporal characteristics from the input signal in the form of interspike intervals. Figure 6(b) shows a zoomed part of the ECG signal. As discussed before, we use spikes generated from the rising signal only to reduce data transfer without compromising accuracy.
2.2 Spiking NeuronsBased Liquid
This section describes the Liquid, starting with the neuron and synapse model and followed by the Liquid topology. The readout unit is described in Section 2.3.
2.2.1 Neuron and Synapse Model
The Liquid in our approach is built using the Izhikevich neuron (izhikevich2003simple, )
model with the following twodimensional ordinary differential equations
(1) 
(2) 
where is the synaptic current to the neuron. In contrast to simple models such as leaky integrateandfire type neurons, the Izhikevich neuron is more biologically realistic with voltage reset occurring at the peak of the spike (as opposed to occurring at threshold) and an instantaneous reset of the membrane potential when the membrane voltage reaches a cutoff voltage (Equation 2). Additionally, Izhikevich neurons can be efficiently implemented in hardware (demirkol2011low, ). Table 1 reports the parameters used for simulation. The subscript E is used to indicate excitatory neurons i.e., neurons that are glutamatergic, while the subscript I is used for inhibitory type neurons, i.e., neurons with GABAergic synapses. These parameters are chosen to model excitatory neurons as regular spiking neurons and inhibitory neurons as fast spiking neurons izhikevich2003simple .
Class  Parameters  Value 
Neuron  ,  0.02, 0.1 
,  0.2, 0.2  
,  65, 65  
,  8, 2  
Initial synaptic strengths  
ESTDP/ISTDP  ,  0.1, 0.1 
,  20, 20  
,  0.1, 0.1  
,  20, 20  
Homeostasis (Excitatory)  0.1, 10, 35  
Homeostasis (Inhibitory)  0.1, 2, 3.5 
Current in the postsynaptic neuron is computed as
(3) 
where is the total number of presynaptic neurons connected to the postsynaptic neuron, is the current through the synaptic connection between presynaptic neuron and postsynaptic neuron and is the corresponding weight. The synaptic currents are modeled using conductance changes and therefore, decay over time.
2.2.2 Liquid Topology
Connection Probability and Synaptic Delay: Our Liquid topology consists of three layers – input, recurrent and output. The first layer is the input layer, which generates spikes (encoded directly from the input ECG). The second layer is the recurrent layer and consists of recurrently connected neurons, where is the number of excitatory neurons and is the number of inhibitory neurons. Motivated by the anatomy of a mammalian cortex (abeles1991corticonics, ; bock2011network, ), our framework consists of . Neurons in the second layer are interconnected satisfying the following – excitatory neurons of the second layer can connect to any neurons (excitatory as well as inhibitory neurons), but the inhibitory neurons are only connected to other excitatory neurons. This decision is influenced by lack of evidence of recurrence in inhibitory neural assembles in mammalian cortex (abeles1991corticonics, )
. Connection probabilities among the neurons are as follows: excitatory to excitatory neurons have connection probabilities of 0.01; excitatory to inhibitory and inhibitory to excitatory neurons have connection probabilities of 0.1. These probabilities are similar to that used in
(pfeil2013six, ). Initial synaptic strengths are reported in Table 1. Changes in synaptic strength are bounded between 0 and . Synaptic connection delays are selected randomly between 1ms and 2ms.Connection Pruning implementing Soft WinnerTakeAll: Inhibitory to excitatory connections are pruned such that an inhibitory neuron is never connected to the excitatory neuron from which it receives a connection. This creates lateral inhibition in the Liquid topology. As an additional step, we fine tuned excitatory and inhibitory synaptic conductances, such that the lateral inhibition is neither too weak (which means lateral inhibition has no influence to distinguish QRS complexes), nor too strong (which means once a winner is selected, the winner prevents other neurons from firing). This implements a soft version of the winnertakeall strategy.
Output Connections: The third layer in the LSM Liquid consists of neurons which are interpreted in the readout unit. In most LSM architectures, the output neurons reside inside the readout. Neurons in the input layer are connected to the excitatory neurons of the second layer using plastic synapses, which are adjusted using interspike intervals (discussed in the next section). Learning of connections from the recurrent layer to the output layer is discussed in Section 2.3 as part of the readout unit.
2.2.3 Learning using STDP and Homeostatic Plasticity
In our framework, synaptic weight updates are disabled after a time interval interval . The time interval 0 (start of ECG) to is training phase of the spiking neural network. In this phase, weights are updated using spike timing dependent plasticity (STDP) (rao2001spike, ; brader2007learning, ). STDP uses correlation between the spikes (based on spike times) to derive a potential causal (or anticausal) relationship between the pre and postsynaptic neurons. This correlation is then used to influence the weight changes (sjostrom2010spike, )^{4}^{4}4
This is contrary to supervised learning such as using error backpropagation
(bohte2002error, ), where a teacher signal is used to adjust the weights.. Synaptic weights connecting an excitatory neuron to any other neurons are updated using excitatorytype STDP (ESTDP), while synaptic weights connecting an inhibitory neuron to an excitatory neurons are updated using inhibitorytype STDP (ISTDP). In this work, STDP (ESTDP/ISTDP) is realized using an exponential function(4) 
To prevent runaway behavior, homeostatic synaptic scaling (carlson2013biologically, ) is used. The combined weight change equation is given by
(5) 
where is the homeostatic scaling factor, is the STDP scaling factor, is the average firing rate of a neuron over a large period of time, is the predefined firing rate of the neuron (a design choice), (and ) is the prebeforepost (and preafterpost) weight change contributions and is the scalability factor
(6) 
where is the time duration over which the firing rate is being averaged and is a constant factor. Our exploration involves tuning the parameters of these equations.
2.2.4 Implementation Framework of Spiking Neural Network
The spiking neural network (spiking neurons and learning rule) is implemented using CARLsim (beyeler2015carlsim, ), a GPUaccelerated library for simulating spiking neural network models with a high degree of biological detail. We have performed exploration with different number of excitatory and inhibitory neurons. Based on these explorations, the spiking network in our approach is built with 64 excitatory and 16 inhibitory neurons.
2.3 Probabilistic Heartrate Decoder as LSM Readout
The heartrate decoder of Figure 2 is the readout stage of our implementation of LSM computation model. Novelty of this readout is its ability to personalize to a subject (with and without cardiac irregularities) in an unsupervised approach, learning directly from the ECG signal processed using the Liquid. For this purpose we use Fuzzy cMeans clustering on a subset of output neurons (called the winning set). Cluster center selection and winning neurons selection subproblems are jointly solved using our proposed particle swarm optimization technique. This optimization problem is solved once during initialization; subsequently, the readout evaluates cluster membership only using the evaluated centers^{5}^{5}5This approach will be extended in future to periodically carry out the initialization (optimization), to take into account major shift in ECG pattern during lifetime of a device or multiple user of the same device.. Finally, we compute discrete heartrate probability from cluster membership using Poisson Binomial Distribution, which is then used to calculate the expected heartrate.
We start our formulation by defining two intervals – the spike integrate interval (SI) and the heartrate inference interval (HI). Spike responses from a neuron are integrated (i.e., counted) in every SI interval. In the remainder of this work, the term spike response is used to denote the spike count of a neuron in a SI interval, which is the smallest unit of time for our purpose. In our experimental setup, SI = 100 ms and HI = 1 min, allowing the heartrate estimated in a HI interval to be directly interpreted as beats per minute (bpm). Each HI interval is composed of 600 SI intervals. Spike responses from the neurons in a HI interval are clustered into two classes (QRS and noQRS) using the Fuzzy cMean algorithm (bezdek1984fcm, )
, where cluster membership of a response is given by the probability distribution over the cluster. We define probability of success for a neuron response as the probability that the neuron response belong to the QRS cluster. Using the Fuzzy cMean outcomes as Bernoulli trails each with its characteristics probability of success, we compute the discrete probability distribution, which describes a fixed number of successes in the
HI interval. This distribution is then used to compute the expected heartrate for the HI interval. Following sections formally describe these steps.2.3.1 Fuzzy cMeans Clustering of Spike Count
The Fuzzy cMeans (FCM) (bezdek1984fcm, )
is a probabilistic clustering algorithm, where instead of computing the absolute membership of an item belonging to a cluster (hard decisions as in K Means clustering
(macqueen1967some, )), FCM computes the likelihood of an item belonging to a cluster (soft or fuzzy decisions).We define the following notations to ease the problem formulation.
The FCM algorithm minimizes the generalized least square error
(7) 
where is the likelihood of belonging to cluster , is a positive definite weight matrix determined by the type of cluster. In this work we consider hyperspherical cluster with
, the identity matrix. Under this condition, the FCM algorithm minimizes the least square Euclidean norm of Equation
7. As discussed in (bezdek1984fcm, ), at the local minimum of ,(8) 
Due to the cyclic dependency of cluster centers on the likelihood and viceversa, Equation 2.3.1 is solved iteratively for a fixed number of iterations (cannon1986efficient, ).
2.3.2 Dimensionality Reduction And Cluster Center Selection
Unlike the native FCM, which involves selecting the cluster centers for a fixed dimension, our approach involves selecting (1) the optimum set of neurons (dimensionality reduction), together with (2) the cluster centers^{6}^{6}6, while the cluster centers ; solution to Equation 7 is therefore of mixedinteger type.
. The problem of finding the optimum dimension is also referred to as feature selection
(keogh2011curse, ). Standard feature selection techniques such as sequential forward (backward) selection, SFS (SBS) (whitney1971direct, ; marill1963effectiveness, ) are often associated with high computation cost and are prone to being stuck at local minimum. We use particle swarm optimization (PSO) (eberhart1995new, ), an evolutionary computing technique inspired by social behaviors such as bird flocking and fish schooling. Evolutionary computing techniques are efficient in avoiding to stuck at local optima. Additionally, PSO is computationally less expensive with faster convergence compared to its counterparts such as genetic algorithm (GA) or simulated annealing (SA), which makes PSO the ideal candidate for resource and power constrained wearable devices.
In general, PSO finds the optimum solution to a fitness function . Each solution is represented as a particle in the swarm. Each particle has a velocity with which it moves in the search space to find the optimum solution. During the movement, a particle updates its position and velocity according to its own experience (closeness to the optimum) and also experience of its neighbors. We introduce the following new notations related to PSO.
Position and velocity updates are performed according to the following equations.
(9) 
where is the iteration number, are constants and (and ) is the particles own (and neighbors) experience. Figure 7 shows the iterative solution of PSO, where position and velocity are updated until a predefined convergence is achieved.
PSO for dimensionality reduction: We introduce ,
, a weighing factor, one for each output neuron. In order to transform real valued weights to binary values (for on/off decision of a neuron), two binarization functions are possible:

Deterministic Binarization:

Stochastic Binarization:
For simplicity we use deterministic binarization with representing the binarization operation performed on all elements of . Modified neuron spike response is
(10) 
To use PSO for selecting the subset of output neurons, the weights are used as the dimensions of the particles in PSO.
PSO for cluster centers: The number of clusters for the FCM algorithm is 2 (QRS and noQRS) for the heartrate estimation problem. This allows us to rewrite Equation 7 using Equation 10 as
(11) 
In order to use PSO to find the cluster centers there are two options:
Using likelihoods as PSO dimensions results in each particle in the swarm to be of dimensions, causing the algorithm to converge at a slow rate to the optimum solution. We therefore use cluster centers as dimensions.
Combining, we have two cluster centers (each of dimension ) together with output neuron weights as combined PSO dimensions i.e., . Under this assumption, the position of a swarm particle at iteration is
(12)  
For every iteration of the PSO algorithm, a particle’s position is used to determine the set of neurons and the centers of the two clusters using Equation 12. Subsequently, each observation is modified according to Equation 10, which is used to evaluate the fitness function using Equation 2.3.2. When the PSO algorithm converges, Equation 2.3.1 is used to compute the probability of belonging to cluster 1 (QRS cluster). The probability of belonging to cluster 0 (no QRS cluster) is .
2.3.3 Fuzzy cMeans Outcome to Discrete Heartrate Probability Distribution
In this section we use the probabilities to infer heartrate. A naive solution is to use hard assignments as shown below.
(13) 
In simple terms, this equation is interpreted as follows: if the probability of an observation belonging to the QRS cluster is higher than the probability of belonging to the noQRS cluster, then the observation is interpreted as QRS and accounted for heartrate estimation. In this work, we take a different approach. Rather than hard assignments as in Equation 13, we compute the discrete heartrate distribution from the probabilities computed using the FCM algorithm. In other words, we consider all observations (even the ones with low probability of association to the QRS cluster) to contribute towards heartrate estimation. For this purpose, let
be the random variable associated with the observation
with probability of success (i.e., belonging to QRS cluster) . The ’s are independent random variables but are not identically distributed, as ’s can be different. This allows us to consider ’s as the outcome of Bernoulli’s trials, with the sum of these variables (representing the heartrate) forming the Poisson Binomial Distribution (le1960approximation, ) of order .(14) 
We are interested in estimating the probability mass function (pmf) of , which gives the discrete probability distribution of heartrate. Let be the set of unique observations in the HI interval with . Then,
(15) 
where the above sum is evaluated over all possible ways of selecting observations from . The solution to Equation 15 involves equating all the permutations for heartrate , which is intractable for larger and
. We use Discrete Fourier Transform (DFT) of the characteristic function to determine the probability mass function (PMF) of the distribution
(zhang2017algorithm, ; hong2013computing, ). For notational brevity, let be the maximum heartrate of humans. PMF of the distribution be represented as , where(16) 
Characteristic function of the random variable is computed as
(17) 
The above equation can be rewritten as
(18) 
The characteristic function can also be written in terms of the PMF as
(19) 
(20) 
Expressing , Equation 20 can be rewritten as
(21) 
where
(22) 
Observe the similarity of Equation 21 to that of inverse discrete Fourier Transform (IDFT), allowing us to compute the PMFs as the discrete Fourier Transform (DFT) of the characteristic function .
The expected heartrate can be represented as
(23) 
Database  Subjects  Duration (mins)  Sampling Rate (Hz) 

Internal (intdb)  8  10  256 
MIT BIH Supraventricular Arrhythmia Database (svdb)  14  30  128 
MITBIH Long Term Database (ltdb)  7  60  128 
Long Term ST Database (ltstdb)  31  60  250 
3 Results and Analysis
We evaluated our approach on ECG recordings from 8 subjects collected during inhouse clinical trials. To demonstrate the applicability of our approach to wider subjects, three other databases are considered from open source Physiobank recordings (goldberger2000physiobank, ). Details of these databases are reported in Table 2. It is to be noted that subjects from our internal database consists of healthy subjects. Additionally, subjects from svdb have cardiac arrhythmia, subjects from ltdb have normal sinus rhythms and subjects from ltstdb have myocardial ischemia. These databases are chosen to demonstrate the applicability of our approach to subjects with healthy cardiac rhythms and also to the ones with different types of cardiac irregularities.
3.1 Inhouse Clinical Data Collection Setup
Inhouse ECG data were collected as part of the variability study of ECG and PPG heartrate intervals. These data were collected using Osram PPG module (test) and gtec (g.USBamp) biosignal amplifier and acquisition system (reference). All the recruited subjects have no history of cardiac disorders. The experimental setup consists of

Osram module (SFH7060) with three green and one each of red and infrared LEDs.

Texas Instruments (TI) analog front end (TI AFE 4403 EVM; evaluation kit).

gtec (g.USBamp) biosignal amplifier and acquisition system (CEcertified and FDAlisted medical device, safety class: II, conformity class: IIa).

Workstation/laptop installed with gtec and TI data acquisition software.

Osram PPG module mounted on a Velcro band for strapping to subject’s wrist during measurements.

3 pregelled snap electrodes (ECG measurements).
Data collection methods for the public databases are partially described in goldberger2000physiobank for svdb and ltdb and in jager2003long for ltstdb.
3.2 Spiking Neural Network Simulation Setup
Our simulation^{7}^{7}7The simulation setup is used here for accelerating ECG simulations. Our final approach is envisioned to be running on a neuromorphic chip such as the CxQuad indiveri2015neuromorphic for realtime ECGbased estimation. framework consists of

A Python module to convert ECG to spikes.

The core CARLsim simulator (in C++) implementing the Liquid.

A Python readout implementing unsupervised heartrate estimation.
The simulation is carried on a system with the following configuration

CPU: Intel Corei7 8 cores

GPU: Nvidia GeForce GTX970

Memory: 32GB

OS: Ubuntu 14.04 (64 bits)
Heartrate estimated using our approach is compared with the golden heartrate calculated using van2015345 , which uses computation intensive (and power hungry) deterministic QRSdetection techniques to compute the heartrate. In the following subsections, we present our detailed evaluations.
3.3 HeartRate Estimation Accuracy
Figure 8 plots the average heartrate estimation accuracy using our approach for all subjects from four databases. The accuracy number for each subject is measured as Mean Average Percent Error (MAPE) calculated as
(24) 
where is the L1 norm or the absolute difference between the actual heartrate () and estimated heartrate () and is the number of oneminute segments in the given ECG duration. As can be seen from the figure, the MAPE using our approach is less than 10% across all subjects in all databases. For our internal database, the MAPE varies between 0.5% and 2.1%, with an average of 1.2%. In terms of absolute values, these subjects have heartrate differences between 0 and 5.3 beats per minute (bpm), with an average difference of 0.8 bpm.
For svdb database, the MAPE varies between 0.53% and 5.3%, with an average of 2.32%. The heartrate difference of the subjects are in the range of 0.2 bpm to 5.1 bpm, an average of 1.8 bpm. For ltdb, the MAPE varies between 1.53% and 8.07% with an average of 3.54%. Heartrate difference of the subjects in this database are between 1.1 bpm and 5.3 bpm, an average 2.6 bpm. Finally for the ltstdb, the MAPE varies between 0.92% and 6.7% with an average of 3.88%. The actual heartrate difference varies between 0.6 bpm and 6.8 bpm, with an average difference of 3.1 bpm.
3.4 Intuition Behind Misprediction and Possible Solution
As discussed in Section 2.1, we use a fixed value for the threshold gap Delta and the timer interval. While this fixed parameter spike encoding scheme works well for identifying most QRS complexes (high accuracy as seen in the previous subsection), there are some scenarios where our approach leads to misprediction (misdetection or false detection, causing false negative or false positive, respectively). In this subsection, we discuss this and propose a simple solution to improve accuracy. Figure 9(a) shows two waveforms (A and B). In this figure are the fixed time instances when these waveforms are compared with the thresholds. Spikes generated from the two waveforms are shown at the bottom of this figure. As seen from this figure, spike trains of the two waveforms are similar. These causes aliasing effect. In other words, the spiking neural network will learn to treat the two waveforms similarly (both as QRS or both as not QRS). This results in misprediction in some cases as seen from Figure 8.
One technique to overcome this situation is to adapt the timer interval in Figure 3 in response to the input waveform. This is shown in figure 9(b), where the timer interval for waveform A is reduced to better track the fast rising waveform A. Spike trains generated from the waveform are shown in the bottom of this figure. It is to be noted that the timer interval for waveform B is unaltered (same as in Figure 9)(a). This results in the spike train shown at the bottom for this waveform. It can be seen quite clearly that spike trains for waveforms A and B are now different. This will allow our spiking neural network to better distinguish the two spike trains, further improving accuracy. One approach to adapt timer intervals is by setting its clock frequency to be proportional to the highest Fourier component extracted from the waveform. This will be addressed in future.
3.5 Histogram of Heartrate Difference
Figure 10 plots the difference between the estimated and the actual heartrate for representative subjects from the four databases. Each subfigure plots the absolute heartrate difference in bpm (on the xaxis) and the distribution i.e., the number of one minute segments on the yaxis. As seen for the subject from intdb in Figure 10(a), out of 13 one minute ECG segments, 8 have heartrate difference of 0.7 bpm or lower, another 4 segments have heartrate difference between 0.7 bpm and 2.1 bpm. Finally, there is one segment with heartrate difference between 2.8 and 3.5 bpm. Results are similar for other subjects of the same database. Similarly, for the subject from svdb in Figure 10(b), 26 out of 28 one minute segments have heartrate difference of 1.4 bpm or lower, the two remaining segments have heartrate difference of 2.3 bpm and 2.8 bpm, respectively. Results are similar for the two representative subjects from ltdb and ltstdb, which are shown in Figures 10(c) and 10(d), respectively.
3.6 Comparison with other learningbased heartrate estimation techniques
In this subsection, we compare other artificial intelligence techniques applied to the heartrate estimation problem. Specifically, we use support vector machine (SVM) implemented with Scikit Python toolbox (pedregosa2011scikit, )
and convolutional neural network (ConvNet) implemented in tensorflow
(abadi2016tensorflow, ), both of which are supervised learning. These toolboxes are used to implement reduced versions of the approaches presented in (magrans2016myocardial, ) and (acharya2017automated, ), respectively. In contrast to these approaches, ours is an unsupervised approach. Summary of these results are presented in Table 3. In this table, Column 1 reports the method used for heartrate estimation; Column 2 reports the underlying algorithm; Columns 36 report the average heartrate estimation accuracy (in terms of MAPE with respect to using van2015345 ) across all subjects in the four databases. It is to be noted that supervised approaches are sensitive to training set selection. To showcase this, we report four numbers each, for the two supervised approaches. The first number is the accuracy obtained when trained with ECG data from regular subjects (i.e., subjects with no cardiac irregularities). The second number is the accuracy when the supervised approaches are trained with subjects from svdb (with arrhythmia). The third number is the accuracy when trained with subjects from ltstdb (with myocardial ischemia). Finally, the fourth number reports the accuracy when trained with representative subjects with and without cardiac irregularities. The training set is the same for the two supervised approaches, while test set is the same for all three approaches.Methods  Approach  intdb  svdb  ltdb  ltstdb 

SVM (magrans2016myocardial, )  Supervised  2.0/6.9/8.2/5.2  9.8/2.3/8.1/4.4  3.8/9.4/11.9/6.9  11.6/14.1/3.7/8.4 
ConvNet (acharya2017automated, )  Supervised  2.5/5.5/10.0/5.7  11.1/2.1/10.1/4.6  3.6/8.9/8.9/3.9  10.3/21.1/3.7/7.4 
Our Approach  Unsupervised  1.2  2.32  3.54  3.88 
Following observations can be made from results presented in Table 3. SVM achieves better accuracy on ltdb and some subjects from intdb, when the classifier is trained with regular subjects (first numbers in row 2, columns 3 and 5). For svdb and ltstdb, the classifier performs poorly (MAPE of approximately 10% and 12% for these two databases, respectively as seen from first numbers in row 2, columns 4 and 6). On the other hand, when SVM is trained with subjects from arrhythmia database (svdb), accuracy for this database increases (2.3%, second number in row 2, column 4). However, accuracies for the other databases are lower (second numbers in row 2, columns 3, 5, 6). Similarly, accuracy is better for ltstdb, when SVM is trained with subjects of this database (third number in row 2, column 6). In this case also the accuracies for other three databases are low (third numbers in row 2, columns 3, 4, 5). This trend is consistent for the convolutional neural network based classifier. From these trend its can be concluded that supervised training approaches need to be specialized for the particular cardiac condition, in order to achieve a good heartrate estimation accuracy.
A second observation can be made from the above results. If the SVM is trained with representative subjects from all databases, average accuracy across these databases increases. However, many subjects in these databases have over 20% error. This is also true for the ConvNet based approach. These results signify that a pretrained model, even with specialized training can result in high error, when applied to unknown subjects. In contrast to these two approaches, our approach learns from the ECG data itself, eliminating the need for a general training set. Additionally, our approach can be applied to new subjects, while consistently providing high estimation accuracy.
3.7 Data Density and Energy Reduction
Table 4 reports the data density and energy reduction summary for subjects in the four databases. For each row item in the table, we report the ECG duration (in sec) in column 2, average spike firing rate (in Hz) in column 3, the data density (bits per spike) in column 4 and the energy consumption (in J) in column 5. The average spike firing rate is the total number of spikes generated in the spiking neural network averaged over the ECG duration (column 2). Data density measures efficiency of the spike encoder and is defined as
(25) 
where input spikes is the number of spikes generated at the output of our spike encoder (and input to the spiking neural network) and ADC bits is the number of bits per ECG sample transmitted post analog to digital conversion (ADC) for QRS detection in a standard nonspiking manner. Finally, column 5 reports the energy reduction achieved using our approach compared to the approach of (deepu2016ecg, ). The energy consumption of our approach is for a stateoftheart spiking hardware – CxQuad (indiveri2015neuromorphic, ), using estimation methodology similar to (cao2015spiking, ). For fairness of comparison, the energy consumption using deepu2016ecg excludes the ADC energy, while the energy consumption using ours exclude the spike encoder energy. The table reports result for each subject in the intdb (rows 29, identified by intdb.S0x in column 1). Average results of these subjects is reported next (row 10, identified by intdb.avg). Finally, average across the subjects of the three other databases (svdb, ltdb and ltstdb) are reported in the next three rows (rows 1113, identified by svdb.avg, ltdb.avg and ltstdb.avg, respectively).
Subjects  ECG duration (sec)  Avg. Spike Firing Rate (Hz)  Data density (bits/spike)  Energy Reduction 

intdb.S01  570  7.2  40.6  38.7x 
intdb.S02  940  7.8  35.5  35.6x 
intdb.S03  940  8.5  35.5  32.7x 
intdb.S04  1060  5.7  63.8  48.6x 
intdb.S05  1060  7.8  33.8  35.6x 
intdb.S06  960  7.8  56.3  36.6x 
intdb.S07  720  8.9  39.5  31.2x 
intdb.S08  980  8.9  44.5  31.3x 
intdb.avg  896  7.8  43.7  36.2x 
svdb.avg  1790  7.6  52.6  37.1x 
ltdb.avg  9900  7.8  40.7  35.9x 
ltstdb.avg  9900  7.7  57.6  36.5x 
As seen from the table, the average data density (averaged over all subjects) for our internal database (intdb) is 43.7. This result can be interpreted as follows: on average for every 43.7 bits of raw ECG data transmitted from the sensor for heartrate detection using standard QRSdetection techniques, our spike encoder transmits one spike for the same purpose. This compression saves energy and databandwidth.The higher the density, the higher the savings. For the three other databases i.e., svdb, ltdb and ltstdb, our approach achieves data density of 52.6 (12.9 – 113.3 for all subjects), 40.7 (24.4 – 56.5 for all subjects) and 57.6 (21.7 – 108.5 for all subjects), respectively.
The average energy reduction using our heartrate estimation approach (averaged over all subjects) for our internal database is 36.2 lower than (deepu2016ecg, ), signifying the importance of our approach for power constrained wearable devices. For the three other databases, the average energy savings of the subjects are 37.1x (32.9x – 49.5x), 35.9x (33.7x – 38.1x) and 36.5x (29.7x – 46.6x), respectively. These results suggest that our approach can be very well integrated in future wearable devices, providing significant battery life and improving user experience.
Subjects  Accuracy  False Positive  False Negative  ECG Detection Offset  

0ms  50ms  50ms  100ms  100ms  200ms  
intdb.S01  100% (414)  2.4%  0%  100%  0%  0% 
intdb.S02  100% (842)  1.9%  0%  100%  0%  0% 
intdb.S03  99.38% (971)  0.6%  0.6%  57.3%  39.3%  3.4% 
intdb.S04  100% (702)  3.2%  0%  100%  0%  0% 
intdb.S05  99.7% (949)  2.5%  0.3%  84.2%  12.3%  3.3% 
intdb.S06  99.7% (934)  1.39%  0.3%  82.7%  17.3%  0% 
intdb.S07  99.0% (907)  1.0%  1.0%  49.7%  39.7%  10.6% 
intdb.S08  100% (1169)  1.19%  0%  100%  0%  0% 
3.8 Extension of our approach to other usecases: QRS Detection Case Study
Our proposed approach can be easily extended to other usecases by designing a new readout using the same Liquid. One most common usecase is the QRS detection, which is commonly used to detect clinically significant heart irregularities such as arrhythmia. In this section we show how our approach can be instantiated for QRS detection. In this case, the probability distribution at the output of the FCM algorithm can be inferred as QRS based on hard assignments (similar to Equation 13). In other words, if the probability of an observation being part of QRS cluster is higher than the probability of it to be part of noQRS cluster, the observation is labeled as QRS. This can be implemented as a separate readout using the same Liquid. Table 5 reports QRS summary for the eight subjects in the intdb. In column 2 we report the detection accuracy in percentage for the actual number of QRS reported in parenthesis. Column 3 and 4 report respectively, the false positives and false negatives as percentage of the actual number of QRS. Column 57 report the offset between the detected QRS position and the position of actual QRS. Column 5 reports the percentage of the detected QRS with a maximum offset of 50ms, Column 6 that between 50ms to 100ms and Column 7 with offset between 100ms and 200ms. As can be seen from the table, our approach results in QRS detection accuracy of over 99% for all subjects in the database. Additionally, the false positives and negatives are less 4% and 1%, respectively. Finally, over 90% of the detected QRS have offset less than 100ms. It is to be noted that results using the three other databases (svdb, ltdb and ltstdb) are similar to that obtained for the intdb. These results signify the relevance of our approach when applied to QRS detection. Our ongoing work is to investigate other usecases, such as arrhythmia detection and ECGbased human authentication.
3.9 Spiking Neural Network Related Results
This section provides results pertaining operation of the spiking neural network. Figure 11 plots the firing rate of the excitatory and the inhibitory neurons of the spiking neural network Liquid. The average firing rate for the excitatory neurons is 45.24 Hz, while that for inhibitory neurons is 85.02 Hz. These results are for a specific subject selected randomly from the intdb. Results for subjects from other databases are similar. These firing rates are consistent with the firing rates reported in amit1997dynamics . We performed design space exploration using the parameters mentioned in Table 1 to adjust the firing rates of the neurons. Finally, Figure 12 plots the raster graph with spike times for 16 (out of 64) excitatory and 16 inhibitory neurons.
4 Conclusions and Outlook
We presented an endtoend approach for heartrate estimation in wearables with embedded neuromorphic hardware. Starting from encoding an analog electrocardiogram (ECG) signal directly into spike train (temporal encoding), our approach uses a network of recurrently connected spiking neurons (the Liquid) to process these spike trains for inferring heartrate. Finally, using a probabilistic inference (readout) unit at the output of the neural network, heartrate is estimated in a completely unsupervised manner. The readout unit is designed using swarm intelligence to reduce dimensionality and select cluster centers of a Fuzzy cMeans Algorithm. Thorough evaluation using data from inhouse clinical trials as well as with publicly available databases demonstrates the high accuracy of our approach together with significant data compression and low power consumption. Following are the useful insight from this work.
Relevance and Practicality: Wearable electronic devices such as the activity tracker are gaining popularity because of their ability to monitor physical activity, sleep and other behaviors. A basic requirement of these devices is to track heartrate accurately using limited energy for different cardiac conditions. Recent studies show that some of these devices suffer from poor accuracy and large power consumptions. Readers are referred to Evenson2015 ; sasaki2015validation for a review of different activity trackers. Another limitation of these devices is that these devices offer limited flexibility to implement complicated clinically significant usecases, such as detecting arrhythmia. This is due to small form factor of these devices coupled with limited flexibility and energy budgets. We have demonstrated that the endtoend solution has a low energy footprint (see Section 3.7) together with high accuracy (see Section 3.3), using a limited number of neurons and synapses. This approach can be readily embedded in future wearable devices with strict energy and area budgets. The proposed LSM computational model offers flexibility by allowing implementation of clinically significant usecases (as readouts) from the spatiotemporal properties of ECG integrated inside a network of spiking neurons (Liquid). We have demonstrated this using the QRS detection usecase in Section 3.8. In future, we will investigate arrhythmia detection. Additionally, the unsupervised readout is conducive to personalized healthcare, by allowing learning from subjects directly, without requiring costly data annotations to train the network. This allows future wearables to be used seamlessly for subjects with and without cardiac irregularities.
Novelty:
Existing machine learning based ECG processing (QRS detection or heartrate estimation) are primarily implemented using supervised learning. These techniques require good training sets to achieve acceptable accuracy as we demonstrated in Section
3.6. In addition to this, classical supervised approaches cannot be easily generalized to different cardiac irregularities. The power consumption is high due to the requirement of transmitting digitized bits between sensor and devices. To address these limitations, our approach presents three novel contributions: (1) a technique to encode spikes from ECG directly, without requiring to digitize the analog ECG signal and thus achieving over 40x reduction in data density (Section 3.7); (2) a novel learning rule for spiking neural network in a LSM computation model, which can be efficiently implemented on stateoftheart neuromorphic hardware, with over 30x improvement in energy consumption (Section 3.7) over existing approaches; and (3) an unsupervised readout for heartrate estimation casestudy, where Liquid states are selected intelligently using Particle Swarm Optimization to improve Fuzzy cMeans clustering (Section 2.3). Unsupervised learning allows generalization to different cardiac conditions, eliminating the need for data annotation. Additionally, the approach can be readily deployed to subjects with rare cardiac conditions, where ECG data is not always available to train a supervised classifier. An example situation is heartrate observation following Transcatheter Aortic Valve Implantation (TAVI) tamburino2011incidence . TAVI subjects are currently monitored in lab facilities for a longer period after the implantation procedure. Our approach can be readily deployed for these subjects, facilitating early rehabilitation and remote monitoring.Advancing neuromorphic computing: Neuromorphic computing has matured significantly in recent years. This growth is fueled by innovations in materials kuzum2013synaptic ; jo2010nanoscale ; amit1994learning , circuits pfeil2013six ; snider2011synapses ; chicca2014neuromorphic ; walter2015neuromorphic , algorithms beyeler2013categorization ; shrestha2017robust ; weng2013modulation and applications CorradiTBME ; verstraeten2005isolated ; carneiro2013event . Our work falls in the interface between neuromorphic applications and algorithms. Towards this end, Table 6 reports some practical applications for neuromorphic computing based on LSM computation model. The table also indicates the contribution of this work. As can be seen, our proposed application contributes towards a reallife benchmark against which future neuromorphic architectures and algorithms can be evaluated.
Approach  Application Domain  Information Coding  Algorithm 

Verstraeten et al. verstraeten2005isolated  Speech Recognition  Temporal  LSM with supervised readout 
Grzyb et al. grzyb2009facial  Expression Recognition  Spatial Coding  LSM with supervised readout 
Corradi et al. CorradiTBME  EEGbased BMI  Temporal  LSM with supervised readout 
Our work  ECGbased heartrate estimation  Temporal  LSM with unsupervised readout 
Acknowledgments
We would like to acknowledge our longterm collaborators Dr. Federico Corradi and Prof. Giacomo Indiveri from the Institute of Neuroinformatics at University of Zurich for the useful discussions related to this work. This work is supported in parts by EUH2020 grant NeuRAM3 Cube (NEUral computing aRchitectures in Advanced Monolithic 3DVLSI nanotechnologies) and ITEA3 proposal PARTNER (Patientcare Advancement with Responsive Technologies aNd Engagement togetheR).
References
References
 (1) D. Phan, L. Y. Siong, P. N. Pathirana, A. Seneviratne, Smartwatch: Performance evaluation for longterm heart rate monitoring, in: Bioelectronics and Bioinformatics (ISBB), 2015 International Symposium on, IEEE, 2015, pp. 144–147.
 (2) R. M. Aarts, M. Ouwerkerk, Apparatus for monitoring a person’s heart rate and/or heart rate variation; wristwatch comprising the same, uS Patent App. 12/097,548 (Nov. 14 2006).
 (3) J. Penders, J. van de Molengraft, M. Altini, F. Yazicioglu, C. Van Hoof, A lowpower wireless ecg necklace for reliable cardiac activity monitoring onthemove, in: Proceedings of the International Conference of the IEEE Engineering in Medicine and Biology Society, 2011.
 (4) S. C. Mukhopadhyay, Wearable sensors for human activity monitoring: A review, IEEE sensors journal 15 (3) (2015) 1321–1330.
 (5) B. Gyselinckx, C. Van Hoof, J. Ryckaert, R. F. Yazicioglu, P. Fiorini, V. Leonov, Human++: autonomous wireless sensors for body area networks, in: Custom Integrated Circuits Conference, 2005. Proceedings of the IEEE 2005, IEEE, 2005, pp. 13–19.
 (6) N. Van Helleputte, M. Konijnenburg, J. Pettine, D.W. Jee, H. Kim, A. Morgado, R. Van Wegberg, T. Torfs, R. Mohan, A. Breeschoten, et al., A 345 w multisensor biomedical soc with bioimpedance, 3channel ecg, motion artifact reduction, and integrated dsp, IEEE Journal of SolidState Circuits 50 (1) (2015) 230–244.
 (7) Z. Krasauskas, L. Telksnys, Ubiquitous personal heart rate long distance transmission to the treatment centers based on smart mobile phone application, in: 2015 IEEE 3rd Workshop on Advances in Information, Electronic and Electrical Engineering (AIEEE), 2015, pp. 1–4. doi:10.1109/AIEEE.2015.7367297.
 (8) B.U. Kohler, C. Hennig, R. Orglmeister, The principles of software qrs detection, IEEE Engineering in Medicine and Biology Magazine 21 (1) (2002) 42–57.
 (9) C. Otto, A. Milenkovic, C. Sanders, E. Jovanov, System architecture of a wireless body area sensor network for ubiquitous health monitoring, Journal of mobile multimedia 1 (4) (2006) 307–326.
 (10) N. Ravanshad, H. RezaeeDehsorkh, R. Lotfi, Y. Lian, A levelcrossing based qrsdetection algorithm for wearable ecg sensors, IEEE Journal of Biomedical and Health Informatics 18 (1) (2014) 183–192.
 (11) Z. Zidelmal, A. Amirou, D. OuldAbdeslam, A. Moukadem, A. Dieterlen, Qrs detection using stransform and shannon energy, Computer methods and programs in biomedicine 116 (1) (2014) 1–9.

(12)
A. Karimipour, M. R. Homaeinezhad,
Realtime
electrocardiogram pqrst detectionâ€“delineation algorithm based on
qualitysupported analysis of characteristic templates, Computers in Biology
and Medicine 52 (2014) 153 – 165.
doi:https://doi.org/10.1016/j.compbiomed.2014.07.002.
URL http://www.sciencedirect.com/science/article/pii/S001048251400167X  (13) A. G. Ramakrishnan, A. P. Prathosh, T. V. Ananthapadmanabha, Thresholdindependent qrs detection using the dynamic plosion index, IEEE Signal Processing Letters 21 (5) (2014) 554–558. doi:10.1109/LSP.2014.2308591.

(14)
S. Jain, M. Ahirwal, A. Kumar, V. Bajaj, G. Singh,
QRS
detection using adaptive filters: A comparative study, ISA Transactions 66
(2017) 362 – 375.
doi:https://doi.org/10.1016/j.isatra.2016.09.023.
URL http://www.sciencedirect.com/science/article/pii/S001905781630386X  (15) C. Deepu, X. Zhang, C. Heng, Y. Lian, An ecgonchip with qrs detection & lossless compression for low power wireless sensors, sensors 2 (2016) 3.
 (16) C.I. Ieong, P.I. Mak, M.I. Vai, R. P. Martins, Subw qrs detection processor using quadratic spline wavelet transform and maxima modulus pair recognition for powerefficient wireless arrhythmia monitoring, in: Design Automation Conference (ASPDAC), 2016 21st Asia and South Pacific, IEEE, 2016, pp. 21–22.

(17)
J. Kim, P. Mazumder, Energyefficient hardware architecture of selforganizing map (som) for ecg clustering in 65nm cmos, IEEE Transactions on Circuits and Systems II: Express Briefs PP (99) (2017) 1–1.
doi:10.1109/TCSII.2017.2672789.  (18) T. Tekeste, H. Saleh, B. Mohammad, A. Khandoker, M. Elnaggar, A nanowatt ecg feature extraction engine in 65nm technology, IEEE Transactions on Circuits and Systems II: Express Briefs PP (99) (2017) 1–1. doi:10.1109/TCSII.2017.2658670.
 (19) R. Silipo, C. Marchesi, Artificial neural networks for automatic ecg analysis, IEEE Transactions on Signal Processing 46 (5) (1998) 1417–1425. doi:10.1109/78.668803.

(20)
I. Saini, D. Singh, A. Khosla,
QRS
detection using knearest neighbor algorithm (KNN) and evaluation on
standard ECG databases, Journal of Advanced Research 4 (4) (2013) 331 –
344.
doi:https://doi.org/10.1016/j.jare.2012.05.007.
URL http://www.sciencedirect.com/science/article/pii/S209012321200046X 
(21)
K. Arbateni, A. Bennia,
Sigmoidal
radial basis function ANN for QRS complex detection, Neurocomputing 145
(2014) 438 – 450.
doi:https://doi.org/10.1016/j.neucom.2014.05.009.
URL http://www.sciencedirect.com/science/article/pii/S0925231214005724  (22) R. Magrans, P. Gomis, P. Caminal, Myocardial ischemia event detection based on support vector machine model using qrs and st segment features, in: Computing in Cardiology Conference (CinC), 2016, IEEE, 2016, pp. 405–408.

(23)
U. R. Acharya, H. Fujita, O. S. Lih, Y. Hagiwara, J. H. Tan, M. Adam,
Automated
detection of arrhythmias using different intervals of tachycardia ECG
segments with convolutional neural network, Information Sciences 405 (2017)
81 – 90.
doi:http://dx.doi.org/10.1016/j.ins.2017.04.012.
URL http://www.sciencedirect.com/science/article/pii/S0020025517306539 
(24)
U. R. Acharya, H. Fujita, O. S. Lih, M. Adam, J. H. Tan, C. K. Chua,
Automated
detection of coronary artery disease using different durations of ECG
segments with convolutional neural network, KnowledgeBased Systemsdoi:http://dx.doi.org/10.1016/j.knosys.2017.06.003.
URL http://www.sciencedirect.com/science/article/pii/S0950705117302769  (25) W. Maass, Networks of spiking neurons: the third generation of neural network models, Neural networks 10 (9) (1997) 1659–1671.
 (26) D. V. Buonomano, M. Merzenich, A neural network model of temporal code generation and positioninvariant pattern recognition, Neural Computation 11 (1) (1999) 103–116.
 (27) N. Kasabov, K. Dhoble, N. Nuntalid, G. Indiveri, Dynamic evolving spiking neural networks for online spatioand spectrotemporal pattern recognition, Neural Networks 41 (2013) 188–201.
 (28) N. Iannella, A. D. Back, A spiking neural network architecture for nonlinear function approximation, Neural networks 14 (6) (2001) 933–939.
 (29) P. U. Diehl, D. Neil, J. Binas, M. Cook, S.C. Liu, M. Pfeiffer, Fastclassifying, highaccuracy spiking deep networks through weight and threshold balancing, in: Neural Networks (IJCNN), 2015 International Joint Conference on, IEEE, 2015, pp. 1–8.
 (30) P. U. Diehl, M. Cook, Unsupervised learning of digit recognition using spiketimingdependent plasticity, Frontiers in computational neuroscience 9 (0) (2015) 0–0.

(31)
A. Samadi, T. P. Lillicrap, D. B. Tweed,
Deep learning with dynamic
spiking neurons and fixed feedback weights, Neural Computation 29 (3) (2017)
578–602, pMID: 28095195.
arXiv:http://dx.doi.org/10.1162/NECO_a_00929, doi:10.1162/NECO_a_00929.
URL http://dx.doi.org/10.1162/NECO_a_00929  (32) G. Indiveri, E. Chicca, R. Douglas, A vlsi array of lowpower spiking neurons and bistable synapses with spiketiming dependent plasticity, IEEE transactions on neural networks 17 (1) (2006) 211–221.
 (33) S. Fusi, M. Mattia, Collective behavior of networks with linear (vlsi) integrateandfire neurons, Neural Computation 11 (3) (1998) 633–652.
 (34) H.Y. Hsieh, K.T. Tang, Vlsi implementation of a bioinspired olfactory spiking neural network, IEEE transactions on neural networks and learning systems 23 (7) (2012) 1065–1073.
 (35) C. Mead, Neuromorphic electronic systems, Proceedings of the IEEE 78 (10) (1990) 1629–1636.
 (36) F. Akopyan, J. Sawada, A. Cassidy, R. AlvarezIcaza, J. Arthur, P. Merolla, N. Imam, Y. Nakamura, P. Datta, G.J. Nam, et al., Truenorth: Design and tool flow of a 65 mw 1 million neuron programmable neurosynaptic chip, IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems 34 (10) (2015) 1537–1557.
 (37) G. Indiveri, F. Corradi, N. Qiao, Neuromorphic architectures for spiking deep neural networks, in: Electron Devices Meeting (IEDM), 2015 IEEE International, IEEE, 2015, pp. 4–2.
 (38) N. K. Kasabov, Neucube: A spiking neural network architecture for mapping, learning and understanding of spatiotemporal brain data, Neural Networks 52 (2014) 62–76.
 (39) M. M. Khan, D. R. Lester, L. A. Plana, A. Rast, X. Jin, E. Painkras, S. B. Furber, Spinnaker: mapping neural networks onto a massivelyparallel chip multiprocessor, in: Neural Networks, 2008. IJCNN 2008.(IEEE World Congress on Computational Intelligence). IEEE International Joint Conference on, Ieee, 2008, pp. 2849–2856.
 (40) B. V. Benjamin, P. Gao, E. McQuinn, S. Choudhary, A. R. Chandrasekaran, J.M. Bussat, R. AlvarezIcaza, J. V. Arthur, P. A. Merolla, K. Boahen, Neurogrid: A mixedanalogdigital multichip system for largescale neural simulations, Proceedings of the IEEE 102 (5) (2014) 699–716.
 (41) J. Schemmel, D. Briiderle, A. Griibl, M. Hock, K. Meier, S. Millner, A waferscale neuromorphic hardware system for largescale neural modeling, in: Circuits and systems (ISCAS), proceedings of 2010 IEEE international symposium on, IEEE, 2010, pp. 1947–1950.
 (42) C. D. Schuman, T. E. Potok, R. M. Patton, J. D. Birdwell, M. E. Dean, G. S. Rose, J. S. Plank, A survey of neuromorphic computing and neural networks in hardware, arXiv preprint arXiv:1705.06963 0 (0) (2017) 0–0.
 (43) W. Maass, T. Natschläger, H. Markram, Realtime computing without stable states: A new framework for neural computation based on perturbations, Neural computation 14 (11) (2002) 2531–2560.
 (44) R. P. Rao, T. J. Sejnowski, Spiketimingdependent hebbian plasticity as temporal difference learning, Neural computation 13 (10) (2001) 2221–2237.
 (45) J. M. Brader, W. Senn, S. Fusi, Learning realworld stimuli in a neural network with spikedriven synaptic dynamics, Neural computation 19 (11) (2007) 2881–2912.
 (46) J. K. Liu, Learning rule of homeostatic synaptic scaling: Presynaptic dependent or not, Neural computation 23 (12) (2011) 3145–3161.
 (47) M. N. Galtier, G. Wainrib, A biological gradient descent for prediction through a combination of stdp and homeostatic plasticity, Neural computation 25 (11) (2013) 2815–2832.
 (48) K. D. Carlson, M. Richert, N. Dutt, J. L. Krichmar, Biologically plausible models of homeostasis and stdp: stability and learning in spiking neural networks, in: Neural Networks (IJCNN), The 2013 International Joint Conference on, IEEE, 2013, pp. 1–8.
 (49) R. Eberhart, J. Kennedy, A new optimizer using particle swarm theory, in: Micro Machine and Human Science, 1995. MHS’95., Proceedings of the Sixth International Symposium on, IEEE, 1995, pp. 39–43.
 (50) J. C. Bezdek, R. Ehrlich, W. Full, Fcm: The fuzzy cmeans clustering algorithm, Computers & Geosciences 10 (23) (1984) 191–203.
 (51) M. Beyeler, K. D. Carlson, T.S. Chou, N. Dutt, J. L. Krichmar, Carlsim 3: A userfriendly and highly optimized library for the creation of neurobiologically detailed spiking neural networks, in: Neural Networks (IJCNN), 2015 International Joint Conference on, IEEE, 2015, pp. 1–8.
 (52) R. Van Rullen, S. J. Thorpe, Rate coding versus temporal order coding: what the retinal ganglion cells tell the visual cortex, Neural computation 13 (6) (2001) 1255–1283.
 (53) M. C. Van Rossum, A novel spike distance, Neural computation 13 (4) (2001) 751–763.

(54)
A. Tavanaei, A. S. Maida,
A
spiking network that learns to extract spike signatures from speech signals,
Neurocomputing 240 (2017) 191 – 199.
doi:https://doi.org/10.1016/j.neucom.2017.01.088.
URL http://www.sciencedirect.com/science/article/pii/S0925231217303119  (55) F. Corradi, G. Indiveri, A neuromorphic eventbased neural recording system for smart brainmachineinterfaces, IEEE Transactions on Biomedical Circuits and Systems 9 (5) (2015) 699–709. doi:10.1109/TBCAS.2015.2479256.
 (56) P. G. Bahubalindruni, V. G. Tavares, E. Fortunato, R. Martins, P. Barquinha, Novel linear analogadder using aigzo tfts, in: Circuits and Systems (ISCAS), 2016 IEEE International Symposium on, IEEE, 2016, pp. 2098–2101.
 (57) Y. Li, W. Mao, Z. Zhang, Y. Lian, An ultralow voltage comparator with improved comparison time and reduced offset voltage, in: Circuits and Systems (APCCAS), 2014 IEEE Asia Pacific Conference on, IEEE, 2014, pp. 407–410.
 (58) D. Du, K. Odame, A bioinspired ultralowpower spike encoding circuit for speech edge detection, in: Biomedical Circuits and Systems Conference (BioCAS), 2011 IEEE, IEEE, 2011, pp. 289–292.
 (59) D. Balsamo, A. S. Weddell, A. Das, A. R. Arreola, D. Brunelli, B. M. AlHashimi, G. V. Merrett, L. Benini, Hibernus++: a selfcalibrating and adaptive system for transientlypowered embedded devices, IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems 35 (12) (2016) 1968–1980.
 (60) E. M. Izhikevich, Simple model of spiking neurons, IEEE Transactions on neural networks 14 (6) (2003) 1569–1572.
 (61) A. S. Demirkol, S. Ozoguz, A low power vlsi implementation of the izhikevich neuron model, in: New Circuits and Systems Conference (NEWCAS), 2011 IEEE 9th International, IEEE, 2011, pp. 169–172.
 (62) M. Abeles, Corticonics: Neural circuits of the cerebral cortex, Cambridge University Press, 1991.
 (63) D. D. Bock, W.C. A. Lee, A. M. Kerlin, M. L. Andermann, G. Hood, A. W. Wetzel, S. Yurgenson, E. R. Soucy, H. S. Kim, R. C. Reid, Network anatomy and in vivo physiology of visual cortical neurons, Nature 471 (7337) (2011) 177–182.
 (64) T. Pfeil, A. Grübl, S. Jeltsch, E. Müller, P. Müller, M. A. Petrovici, M. Schmuker, D. Brüderle, J. Schemmel, K. Meier, Six networks on a universal neuromorphic computing substrate, Frontiers in neuroscience 7 (0) (2013) 0–0.
 (65) J. Sjöström, W. Gerstner, Spiketiming dependent plasticity, Spiketiming dependent plasticity 35 (0) (2010) 0–0.
 (66) S. M. Bohte, J. N. Kok, H. La Poutre, Errorbackpropagation in temporally encoded networks of spiking neurons, Neurocomputing 48 (1) (2002) 17–37.
 (67) J. MacQueen, et al., Some methods for classification and analysis of multivariate observations, in: Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, Vol. 1, Oakland, CA, USA., 1967, pp. 281–297.
 (68) R. L. Cannon, J. V. Dave, J. C. Bezdek, Efficient implementation of the fuzzy cmeans clustering algorithms, IEEE Transactions on Pattern Analysis and Machine Intelligence (2) (1986) 248–255.

(69)
E. Keogh, A. Mueen, Curse of dimensionality, in: Encyclopedia of Machine Learning, Springer, 2011, pp. 257–258.
 (70) A. W. Whitney, A direct method of nonparametric measurement selection, IEEE Transactions on Computers 100 (9) (1971) 1100–1103.
 (71) T. Marill, D. Green, On the effectiveness of receptors in recognition systems, IEEE transactions on Information Theory 9 (1) (1963) 11–17.
 (72) L. Le Cam, et al., An approximation theorem for the poisson binomial distribution, Pacific J. Math 10 (4) (1960) 1181–1197.
 (73) M. Zhang, Y. Hong, N. Balakrishnan, An algorithm for computing the distribution function of the generalized poissonbinomial distribution, arXiv preprint arXiv:1702.01326 0 (2017) 0–0.
 (74) Y. Hong, On computing the distribution function for the poisson binomial distribution, Computational Statistics & Data Analysis 59 (2013) 41–51.
 (75) A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.K. Peng, H. E. Stanley, Physiobank, physiotoolkit, and physionet, Circulation 101 (23) (2000) e215–e220.
 (76) F. Jager, A. Taddei, G. B. Moody, M. Emdin, G. Antolič, R. Dorn, A. Smrdel, C. Marchesi, R. G. Mark, Longterm st database: a reference for the development and evaluation of automated ischaemia detectors and for the study of the dynamics of myocardial ischaemia, Medical and Biological Engineering and Computing 41 (2) (2003) 172–182.
 (77) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Scikitlearn: Machine learning in python, Journal of Machine Learning Research 12 (Oct) (2011) 2825–2830.
 (78) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., Tensorflow: Largescale machine learning on heterogeneous distributed systems, arXiv preprint arXiv:1603.04467 0 (2016) 0–0.

(79)
Y. Cao, Y. Chen, D. Khosla, Spiking deep convolutional neural networks for energyefficient object recognition, International Journal of Computer Vision 113 (1) (2015) 54–66.
 (80) D. J. Amit, N. Brunel, Dynamics of a recurrent network of spiking neurons before and following learning, Network: Computation in Neural Systems 8 (4) (1997) 373–404.

(81)
K. R. Evenson, M. M. Goto, R. D. Furberg,
Systematic review of the
validity and reliability of consumerwearable activity trackers,
International Journal of Behavioral Nutrition and Physical Activity 12 (1)
(2015) 159.
doi:10.1186/s1296601503141.
URL http://dx.doi.org/10.1186/s1296601503141  (82) J. E. Sasaki, A. Hickey, M. Mavilia, J. Tedesco, D. John, S. K. Keadle, P. S. Freedson, Validation of the fitbit wireless activity tracker for prediction of energy expenditure, Journal of Physical Activity and Health 12 (2) (2015) 149–154.
 (83) C. Tamburino, D. Capodanno, A. Ramondo, A. S. Petronio, F. Ettori, G. Santoro, S. Klugmann, F. Bedogni, F. Maisano, A. Marzocchi, et al., Incidence and predictors of early and late mortality after transcatheter aortic valve implantation in 663 patients with severe aortic stenosis, Circulation 123 (3) (2011) 299–308.
 (84) D. Kuzum, S. Yu, H. P. Wong, Synaptic electronics: materials, devices and applications, Nanotechnology 24 (38) (2013) 382001.
 (85) S. H. Jo, T. Chang, I. Ebong, B. B. Bhadviya, P. Mazumder, W. Lu, Nanoscale memristor device as synapse in neuromorphic systems, Nano letters 10 (4) (2010) 1297–1301.
 (86) D. J. Amit, S. Fusi, Learning in neural networks with material synapses, Neural Computation 6 (5) (1994) 957–982.
 (87) G. Snider, R. Amerson, D. Carter, H. Abdalla, M. S. Qureshi, J. Leveille, M. Versace, H. Ames, S. Patrick, B. Chandler, et al., From synapses to circuitry: Using memristive memory to explore the electronic brain, Computer 44 (2) (2011) 21–28.
 (88) E. Chicca, F. Stefanini, C. Bartolozzi, G. Indiveri, Neuromorphic electronic circuits for building autonomous cognitive systems, Proceedings of the IEEE 102 (9) (2014) 1367–1388.
 (89) F. Walter, F. Röhrbein, A. Knoll, Neuromorphic implementations of neurobiological learning algorithms for spiking neural networks, Neural Networks 72 (2015) 152–167.
 (90) M. Beyeler, N. D. Dutt, J. L. Krichmar, Categorization and decisionmaking in a neurobiologically plausible spiking network using a stdplike learning rule, Neural Networks 48 (2013) 109–124.
 (91) S. B. Shrestha, Q. Song, Robust learning in spikeprop, Neural Networks 86 (2017) 54–68.
 (92) J. Weng, S. Paslaski, J. Daly, C. VanDam, J. Brown, Modulation for emergent networks: Serotonin and dopamine, Neural Networks 41 (2013) 225–239.
 (93) D. Verstraeten, B. Schrauwen, D. Stroobandt, J. Van Campenhout, Isolated word recognition with the liquid state machine: a case study, Information Processing Letters 95 (6) (2005) 521–528.
 (94) J. Carneiro, S.H. Ieng, C. Posch, R. Benosman, Eventbased 3d reconstruction from neuromorphic retinas, Neural Networks 45 (2013) 27–38.
 (95) B. J. Grzyb, E. Chinellato, G. M. Wojcik, W. A. Kaminski, Facial expression recognition based on liquid state machines built of alternative neuron models, in: Neural Networks, 2009. IJCNN 2009. International Joint Conference on, IEEE, 2009, pp. 1011–1017.