Adaptive Beat-to-Beat Heart Rate Estimation in Ballistocardiograms

A ballistocardiograph records the mechanical activity of the heart. We present a novel algorithm for the detection of individual heart beats and beat-to-beat interval lengths in ballistocardiograms (BCGs) from healthy subjects. An automatic training step based on unsupervised learning techniques is used to extract the shape of a single heart beat from the BCG. Using the learned parameters, the occurrence of individual heart beats in the signal is detected. A final refinement step improves the accuracy of the estimated beat-to-beat interval lengths. Compared to many existing algorithms, the new approach offers heart rate estimates on a beat-to-beat basis. The agreement of the proposed algorithm with an ECG reference has been evaluated. A relative beat-to-beat interval error of 1.79% with a coverage of 95.94% was achieved on recordings from 16 subjects.


I. INTRODUCTION
W ORLD-WIDE, cardiovascular diseases are the major cause of death [1]. Unobtrusive long-term monitoring of patient's vital signs shows great promise for the prevention and management of such diseases. Home monitoring, for instance, can help to reduce mortality rates, the amount of time spend in hospitals, and the overall costs of treatment [2]. At the same time, patients do not have to be removed from their familiar environment while adverse developments to their condition can still be detected. Furthermore, unobtrusive vital signs monitoring systems can serve as additional safety measure in the general ward of hospitals. A key requirement for the acceptance of such a system is that it does not reduce the patients' comfort nor increase the burden on the hospital staff.
Ballistocardiography (BCG) allows cardiac activity to be monitored unobtrusively. A ballistocardiograph records the vibrations of the body which are caused by the mechanical activity of the heart. While the basic concept has been known since the late 19th century [3], it has gained renewed interest in recent years. Modern ballistocardiographic systems can be fully integrated into everyday objects such as beds [4]- [8], chairs [9] or weighing scales [10] in ways which make the presence of the system unnoticeable to the patient. Furthermore, wearable BCG system based on accelerometers worn Manuscript  on the subject's chest have been proposed [11]. The main advantage of ballistocardiography over electrocardiography is that, in general, no electrodes, textiles, or the like have to be attached to the patient's body. Therefore, ballistocardiographic systems are well suited for monitoring cardiopulmonary activity at night over longer periods of time. Due to the underlying measurement principle, additional information concerning the general activity level of the patient can also be derived from such a system. Current signal processing algorithms for the estimation of heart rates from BCGs can be divided into two groups: those which only provide heart rates averaged over several seconds and those that detect individual heart beats in the signal.
The latter split the signal into segments. For each segment, the average heart rate is estimated by finding the maximum of the auto-correlation function [12], [13] or the power spectral density [5], [14]. Other approaches average over a number of fiducial points [8] or apply empirical mode decomposition [15]. These algorithms give only coarse estimations of the heart rate. They cannot provide information about beat-to-beat heart rate variability or irregular arrhythmias. Furthermore, they can produce incorrect results when strong harmonics are found in the signal. Beat-to-beat interval information, however, is necessary for advanced applications such as heart rate variability analysis or sleep staging.
Of the algorithms that detect individual heart beats, many work by detecting the largest deflection or other fiducial points in the BCG signal [16]- [20]. In some cases, the signal is pre-processed using a variety of techniques, such as discrete wavelet transform [7], [21]- [23], low-pass filtering [24], wavelet-transform-based de-noising [25], or a Hilbert transform [26]. The algorithms that rely on fiducial points lack robustness on general BCG data, as often dominant deflections are missing from the BCG or change their position with respect to the heart beat.
In [4] and [27], a signal template of a single heart beat is used to detect subsequent beats by calculating the crosscorrelation between the template and the signal. In both works, however, the templates were selected manually, which makes this approach unsuitable for a fully automatic monitoring system. Another approach makes use of a multi-sensor array [28]. The authors apply FFT averaging together with a cepstrum analysis to determine the inter-beat intervals. The method presented in [29] is based on clustering short signal segments centered around minima and maxima of the signal. However, the authors report that their algorithm was not able to find more than 49.2% of all heart beats during their evaluation. The use

Raw signal
Force sensor Time (s) F Fig. 1: Overview of the operating principle of a general bedmounted BCG system. of so-called heart valve components in the BCG signal was suggested in [30]. These components are reportedly related to the closure of the heart valves during the cardiac cycle.
In this article, we present and evaluate a novel method for the estimation of beat-to-beat interval lengths from BCGs: the Beat-to-beat Estimation by Adaptive Training (BEAT) algorithm. The BEAT algorithm works on signals recorded by a single sensor and drops the assumptions of a regularly beating heart. At its core, it is based on unsupervised learning techniques to adapt to the high inter-and intra-subject variability of BCG recordings. The algorithm automatically adapts itself to the given BCG signal using a short (e.g. 30 seconds long) training phase. During this training, a set of parameters describing the properties of an individual heart beat in the given signal is determined. Using these parameters, three independent methods are combined to locate heart beats in the BCG signal. These estimations are further refined to obtain accurate beat-to-beat interval length information. Based on the beat-to-beat interval lengths, a heart rate can be computed for each interval. An initial overview and first results of this method were previously published in [31] and [32].
Section II presents a short background on our acquisition system and general BCG signal characteristics. We continue by discussing the details of our proposed algorithm in Section III. In the final sections, we evaluate our algorithm and discuss the results.

A. Signal Acquisition System
For the purpose of this study, a regular hospital bed was instrumented to allow the sensing of forces acting perpendicular to the surface of the bed (see Fig. 1). This was achieved by attaching four strain gauges to one slat of the bed's slatted frame. Located at the center of the slat, these strain gauges form a full Wheatstone bridge which measures the deformation of the slat. The instrumented slat was installed under the subject's thorax in order to optimally record cardiopulmonary activity. An ordinary mattress was placed on top of the slats. Data were acquired by means of a 12 bit ADC using a sampling rate of 128 Hz. A lead I ECG was simultaneously recorded for reference purposes.

B. Ballistocardiogram
In the 1950s, three axes of measurement for translational BCG recordings were defined [33]: longitudinal (head-foot), transverse (side to side), and dorsoventral (back to chest). Historically, most BCG recording systems were of the longitudinal kind [3], [34]. Weighing scales used in recent studies [10] also measure along this axis. Many recently introduced unobtrusive BCG systems, however, especially those which are bed-based [5]- [8], measure along a combination of the transverse and dorsoventral axes. In these system, the exact axis of measurement, i.e. transversal, dorsoventral, or a combination thereof, depends on the subject's current posture with respect to the sensor. As shown in Fig. 1, the system used in this work belongs to this class of systems. Other examples of systems measuring along the transverse-drosoventral axis include accelerometers worn on the chest [11]. For the sake of brevity, we will deviate from the original 1950s definition and denote modern transverse-dorsoventral BCG systems as transverse in the following discussion.
The distinction between longitudinal and transverse BCGs is important, as the signal morphology differs significantly between the recordings shown in transverse BCG studies and those presented in longitudinal BCG studies. While longitudinal BCGs reportedly have a strong cardiac output related component [34], this is not generally true for transverse BCGs. In fact, the exact combination of sources contributing to the signal components recorded by modern transverse BCG system, such as ours, are not yet well understood. This is partly due to the introduction of difficult to model mechanical components, such as mattresses, into the measurement systems. Hence, the very sources of the measured signal poses interesting research questions for future work in this field.
Furthermore, unlike longitudinal BCGs, transverse BCG recordings are highly variable, much more so than e.g. ECGs. The waveform of a single heart beat in a transverse recording changes depending on the subject and on how the subject is positioned on the bed with respect to the sensor.  shows BCG tracings of two heart beats recorded from two subjects. Each subject was recorded in two different postures. Even tough the transverse BCG seems to come with many strings attached due to its higher morphological variability, its major advantage is that it allows an easy implementation of unobtrusive monitoring systems. Due to this fact, many contributions to the field of ballistocardiography in the recent past have dealt with transverse measurement systems. This is also the reason why we chose to develop an heart beat detection algorithm which specifically addresses the variability that is inherent in these types of BCGs.
Specifically, even if occasionally a characteristic BCG peak seems to coincide with each heart beat, it usually distinguishes itself not enough from other peaks in the signal to be detected reliably. A solution to this problem is to detect signal patterns which repeat themselves with each heart beat and which consist of several peaks instead of a single outstanding peak. Surely, the individual peaks also vary from beat to beat, but the overall pattern is still preserved and can thus be better detected by some algorithm than single peaks.

A. Overview
In the following, we present the Beat-to-beat Estimation by Adaptive Training (BEAT) algorithm which autonomously learns and detects BCG peak patterns corresponding to individual heart beats. An overview of the proposed algorithm is given by the flowchart shown in Fig. 3.
First, the raw BCG signal is pre-processed by applying a second order Butterworth high-pass filter with a 3 dB cut-off frequency of 1 Hz to remove the low-frequency respiratory components. The use of a filter with a cut-off frequency of 1 Hz might appear counter-intuitive since, at 60 bpm, 50 % of the power of the pulse shape will be removed. However, we have found this frequency to be a good trade-off between preserving the pulse shape and suppressing undesired artefacts of the respiratory motion in the filtered signal.
A short segment of the filtered signal is then analyzed to determine the features of the heart beat in the so-called : Parameters used to describe a local maximum in the low-pass filtered BCG signal. The i-th peak is parameterized in terms of its amplitude a max , the amplitude of the local minimum a min between the i-th and the (i + 1)-th peak, the distance d max to this local minimum as well as the distance d min from the local minimum to the local maximum of the (i + 1)-th peak.
"training phase". In the next phase, i.e. the beat detection phase, the remaining signal is scanned for heart beats using the features that were extracted during the training procedure. This results in a list of estimated heart beat locations. In a final step, the estimated heart beat locations are used to produce a refined list of beat-to-beat interval lengths. Whenever subjects enter the bed or change their posture with respect to the BCG sensor, the shape of the heart beat pattern changes. In this case, the training phase is repeated so that the algorithm can adapt to the changed heart beat pattern.

B. Training Phase
During the training phase, a short signal segment with a duration of T train (e.g. T train = 30 s) is analyzed to extract a set of parameters describing the principal morphology of the heart beat in the BCG. First, the pre-processed BCG signal is smoothed by applying a second order Butterworth low-pass filter with a 3 dB cut-off frequency of 10 Hz. In this smoothed signal, maxima and minima are located and parameterized as shown in Fig. 4 in order to obtain a low-dimensional, rough description of the signal.
As outlined before, the presented algorithm tries to detect patterns consisting of several peaks. For this purpose, each peak in the smoothed BCG signal is assigned a feature vector f containing the parameters of the peak itself and of the N −1 consecutive peaks. We set N = 7, as this provides the most accurate results on our dataset. The feature vector assigned to the i-th peak is defined as Each feature vector f i encodes the morphology of a signal segment beginning at the i-th peak and encompassing a total of N peaks. The similarity between two feature vectors can be quantified, for instance, by the Euclidean distance between them. If a group of feature vectors describe a similar signal pattern, they will be located close to each other and thus form a cluster in the feature space. Hence, we identify recurring peak patterns by identifying clusters of similar feature vectors.
Prior to further processing, the feature vectors are normalized. Principal component analysis is applied to reduce the A modified version of the k-means clustering algorithm [35] is then performed to identify clusters of feature vectors. Fig. 5 shows a set of feature vectors forming four clusters in feature space. Each cluster represents a different repeating pattern found in the BCG signal. The described morphology is extracted from each cluster in the form of two parameters: 1) The cluster center c k is the mean of all feature vectors belonging to the k-th cluster.
2) The cluster prototype p k (n) is the signal subsegment from which the feature vector with minimum distance to c k was derived, i.e. a time-domain representation of the signal pattern described by the cluster. In addition to the cluster analysis, a modified version of the heart valve (HV) signal which was presented in [30] is also considered. The heart valve signal is reportedly related to the closure of the heart valves during the cardiac cycle. It is computed from the BCG by first applying a band-pass filter, then squaring the resulting signal, and finally estimating an envelope by low-pass filtering. Instead of using a fixed 20-40 Hz bandpass to extract the heart valve signal, we apply a filter with a narrower, 4 Hz wide, passband and a tunable center frequency. A suitable center frequency is heuristically selected during the training phase.
In combination with the modified heart valve signal, the parameter set (c best , p best (n)) which relates best to the heart beat is identified. This is achieved by using each parameter set separately to detect heart beats in the training segment (see Subsection III-C). The most suitable set of parameters is then determined based on the results of these test runs using a set of heuristics.

C. Detection Phase: Heart Beat Indicators
The following three methods are used to locate individual heart beats in the BCG signal. Each methods makes use of a different parameter that was derived during the training phase. 1) Cross-Correlation: Using the cluster prototype p best (n), further heart beats can be located by computing the crosscorrelation between the prototype signal and the remaining BCG signal. Local maxima in the cross-correlation function are detected and linearly interpolated to obtain an upper envelope of the correlation signal. Clear peaks in this envelope appear at locations where the pattern described by the prototype, i.e. a heart beat, occurs. An exemplary plot of the cross correlation function, its upper envelope, and the indicated heart beat locations are shown Fig. 6.
2) Euclidean Distance: The second method detects heart beat locations by computing the Euclidean distance between the feature vectors of the remaining signal and the cluster center c best . As shown in Fig. 6, the resulting distance function exhibits local minima when a signal segment is similar to the heart beat pattern which is represented by the cluster.
3) Heart Valve Signal: Using the center frequency selected during the training, the heart valve signal is computed. Fig.  6 shows that the maxima in the HV signal coincide with the occurrence of heart beats in the BCG.
D. Detection Phase: Indicator Fusion 1) Reliability Heuristics: Each of the discussed methods produces either a minima or a maxima in its output wherever a heart beat is found in the BCG signal. In order to reasonably merge these outputs, we reduce each output function to a set of indicator pairs. Each indicator pair consists of the time t at which the extreme value occurred and a reliability score w which indicates how reliable this extreme value is as an indicator for a heart beat location. Local maxima in the envelope of the cross-correlation, for example, are evaluated by means of their height relative to their two neighboring minima and maxima. Analogous     to the reliability values for cross-correlation maxima (w c ), reliability heuristics for minima in the distance function (w d ) and for maxima in the heart valve signal (w h ) are computed. By constructing the reliability heuristics based on relative parameters, they are normalized to a range from 0 to 1.
2) Merging Indicators: Ideally, each heart beat will be detected by each method and the resulting indicators would perfectly align. In real world scenarios, however, the indicated locations slightly differ from each other, while spurious indicators and heart beats which are not detected by all of the methods further complicate the situation. We suggest a procedure to coalesce these disparate indicators by identifying and merging indicator triplets which are close to each other and which therefore seem to correspond to the same heart beat.
Each indicator selects up to one partner from each of the other indicator classes. For example, a distance indicator q d should be merged with one cross-correlation indicator q c and one HV signal indicator q h that have a high reliability and that are close to its own position. Possible partner indicators are selected from the interval t d ± 0.33 s. The indicator q d and its partners are merged to form a representative indicator pair which is located at the weighted average of the indicators' locations and whose reliability score is computed as the average of the reliabilities of the indicators.
A representative indicator pair is computed for every indicator, with correlation criterion and HV criterion indicators yielding representative indicators Q c and Q h , respectively. Afterwards, a heart beat score H(t) is computed for every time point t, where H(t) = 0 if no representative indicators are found at time t, otherwise H(t) equals the sum of the reliability scores of the representative indicators located at time t. Figure 7 visualizes this process. If all members of a group select each other, the resulting representative indicators are identical and their reliability scores, therefore, add up. Spurious peaks with low reliabilities, on the other hand, will also form a triplet but no other indicator will select them as partner. Hence, their heart beat scores will be significantly lower than those of triplets that selected each other consensually. The locations of heart beats can then be detected as time points t where the heart beat score H(t) exceeds a certain threshold.

E. Interval Refinement
The accuracy of the beat-to-beat interval lengths is further improved by applying the refinement step proposed in [30]. Given two consecutive heart beat locations, the BCG signal recorded during this interval should represent a characteristic pattern of a single heart beat. We refine each interval individually by determining when the corresponding characteristic pattern first repeats itself using the autocorrelation function.

A. Performance Metrics
The performance of the BEAT algorithm with respect to the estimated heart beat events and the beat-to-beat interval lengths was evaluated using a synchronized lead I ECG as gold standard. After acquisition, the recorded signals were processed off-line by algorithms implemented using the MAT-LAB software package. In the recorded ECGs, the locations of the R peaks were detected using the Pan-Tompkins algorithm [36]. Furthermore, the resulting R-R intervals were computed.
The closest BCG heart beat location within a window of ±0.25 s around each R peak was assigned to the R peak. Unassigned BCG heart beats were considered false positives while R peaks to which no BCG heart beat was assigned were considered false negatives.
However, the number of false positives and false negatives gives no indication on the accuracy of the estimated beat-tobeat intervals. Therefore, for each estimated BCG interval, the absolute and relative error between its length and the length of the corresponding R-R interval was computed (RR abs and RR rel , respectively). Since false negatives and false positives directly affect the lengths of their neighboring beat-to-beat intervals, interval errors are suitable for measuring the overall quality of the heart beat estimation.
Another quality measure is aimed at evaluating the algorithm's adequacy for heart rate monitoring. According to the recommendations of ANSI/AAMI/ISO EC13 (cf. [37]), cardiac monitors are supposed to update their displayed heart rates every 10 seconds. Hence, for each non-overlapping 10 second window, the error between the mean heart rate derived from the BCG during that window and that derived from the ECG is computed (HR 10 ).
For each of the discussed error measures the mean of the absolute errors (x) is reported. Furthermore, the spread of the errors is computed using the 90 th percentile of the errors (P 90 ), i.e. 90 % of the errors lie below this value.
In addition to the accuracy of the estimation, the percentage of the BCG signal that was automatically classified as artefactfree (coverage) is also recorded. Artefacts are detected based on the signal energy in a moving window. Only the signal segments that are artefact-free were considered for the evaluation as the presence of artefacts impedes a reliable heart rate estimation (see Fig. 8).

B. Measurement Scenario
Measurements were performed in the instrumented bed on a group of 16 healthy volunteers (9 male, 7 female, ages: 20-50) who gave their informed written consent to participate in the test. Each test subject was measured for a total duration of 30 minutes switching their positions every 7.5 minutes (left lateral → supine → right lateral → prone).
Segments belonging to the individual postures in bed were identified manually in the resulting recordings and were analyzed separately by the algorithm. Remaining body motion artefacts in the segments were detected automatically by the algorithm and disregarded in the further analysis.
The signal segments were processed by the BEAT algorithm and the results obtained from each subject were combined by averaging them with respect to the segments' lengths. Table I shows the results of the BEAT algorithm for each subject. Furthermore, we also evaluated the performance of the beat-to-beat estimation method proposed in [30] using the same dataset. The results of this comparison are also shown in Table I. The BEAT algorithm was able to derive beat-to-beat heart rates from the BCG signals of all subjects. On average, 95.94 % of each signal was identified as artefact-free and usable for heart rate estimation. The quality of the estimations is rather consistent between all subjects. The average false positive and false negative rates were 0.12 % and 0.41 %, respectively. Further, the mean absolute RR error was 16.61 ms (P 90 : 30.62 ms) while the mean relative error was 1.79 % (P 90 : 3.38 %). This value can be put into perspective by considering that the time interval between two consecutive sampling points in the BCG signal, which is sampled at 128 Hz, is 7.8 ms. This means that the average interval length is off by two samples. The mean error of the heart rates averaged over 10 seconds was 0.46 bpm (P 90 : 1.19 bpm). Figure 9 shows a modified Bland-Altman plot comparing the averaged heart rates to the ECG reference. The plot indicates that there is no systematic bias in the heart rates obtained through the BEAT algorithm. Compared to [30], the BEAT algorithm provides better overall accuracy but slightly lower coverage as shown in Table I. A second comparison was performed using the averaging heart rate estimator presented in [12]. This algorithm provides one heart rate estimate every 30 seconds, hence we averaged  the beat-to-beat heart rates obtained by the BEAT methods using consecutive non-overlapping windows of equal length. Table II shows the results of both algorithms. Compared to the beat-to-beat results in Table I, higher coverages could by achieved since the larger window length makes the algorithms less susceptible to artefacts of short duration. The average heart rate error of the BEAT algorithm (0.5 bpm) is lower than the error achieved by [12] (2.6 bpm). It should be noted, however, that the methods presented in [12] and [30] were developed for the use with an EMFi foil BCG sensor as opposed to the strain gauge sensor used to record the data for this study. Hence, the inferior performance of these methods could be due to the differences in the sensor modalities.

C. Results
We also evaluated the BEAT algorithm's performance with respect to the four postures in which the subjects were recorded. Table III shows the relative beat-to-beat interval error for each posture. While the prone and the left lateral positions give the best results, the overall accuracy does not vary significantly between the different postures. It is not surprising that the supine position produces the worst results as it causes the ribcage and spine of the subjects to be located between the heart and the sensor, thus causing a slight attenuation of the signal.
Further evaluation of the BEAT algorithm was performed with respect to the contribution of the three types of heart beat indicators presented in Section III-C. First, we applied each of the methods separately to detect beat-to-beat interval lengths. Then, we tested combinations of two methods using our proposed merging algorithm. Finally, all three methods were combined. Figure 10 shows the results for each of these tests. By combining multiple methods using the proposed algorithm, the relative interval error can be significantly improved. Combining two methods already provides a significant benefit over using a single method. However, the combination of all three methods provides the best overall results.
The BEAT algorithm was designed to also correctly identify heart beats with irregular beat-to-beat intervals. Figure 11 demonstrates this capability on the recordings taken from subject 11 who seems to exhibit a very pronounced case of respiratory sinus arrhythmia. The BEAT algorithm clearly follows these fluctuations on a beat-to-beat basis. While this is a promising result, further work is necessary to assess the algorithms performance in the presence of arrhythmias which not only affect the timing between beats but also their morphology, such as atrial fibrillations.
While the BEAT algorithms provides accurate heart rate estimations for the majority of the analysed signals, its results can be unreliable if the assumption of a repeating pattern coinciding with each heart beat is not fulfilled. Figure 12 shows a short segment of the BCG signal recorded from subject 09 as well as the computed heart rates and the ECG reference values. The BCG signal in this example exhibits significant irregularities which causes severe errors in the output of the BEAT algorithm. However, these severe irregularities were only observed with subject 09 lying in a supine position. Whether they were caused by external influences or by the subject's physiology could not be determined. We would also like to point out that the dataset used to evaluate the BEAT algorithm presents a best-case situation for heart rate estimation from BCGs. All subjects were healthy and were instructed to move as little as possible between the prescribed turning events. Hence, we expect that, especially, the coverage will be reduced for subjects that are more active and thus cause more motion artefacts. Increasing artefacts are also likely to negatively affect the accuracy of the beat-to-beat estimation.

V. CONCLUSION
A novel algorithm for the estimation of beat-to-beat heart rates from BCGs was presented. The BEAT algorithm uses unsupervised learning techniques to derive information about the heart beat pattern from a short BCG signal segment. This training step is based on the parameterization and collation of local maxima in the BCG as well as on a modified k-means clustering algorithm. Using the parameters derived during the training phase, heart beats are located in the remaining BCG signal by merging the results of three independent indicator functions. The resulting estimated heart beat locations are processed by a refinement step to obtain a list of beat-to-beat interval lengths.
To evaluate the achieved quality of the beat-to-beat heart rate estimations, BCGs and reference ECGs were recorded from 16 test subjects in a laboratory environment. Good agreement between the output of the BEAT algorithm and the gold standard could be observed. We also compared the proposed algorithm to two methods known in the literature with favorable results.
Future work on the BEAT algorithm could focus on finetuning the automatic triggering of the training phase. Furthermore, the performance of the BEAT algorithm on explicit arrhythmia data could be analysed as well as its application to other sensor modalities.