Maximum a posteriori detection of heartbeats from a chest-worn accelerometer

The monitoring of cardiac parameters has a variety of applications in medicine. One example is the assessment of acute-illness severity, where heartrate (HR) is measured to compute an early warning score for patient deterioration (Royal College of Physicians 2017). Another is the monitoring of sleep using cardiac parameters, as an alternative to polysomnography (Bakker et al 2021). For these applications, precise measurement of instantaneous heart rate (IHR) is paramount. IHR is derived from heartbeats, specifically the inter-beat intervals (IBIs) and beat localization times; it is computed as the inverse of the IBIs and can change at every heartbeat. A common characteristic of these applications is the need for prolonged monitoring of patients in bed, either continuously or night-after-night, which is why the measurement must be unobtrusive, be robust against differences in wearing conditions and postures and allow for easy re-attachment of the sensor.

The gold standard for cardiac measurements is electrocardiography, or ECG, in which the QRS complexes are the basis of measurement. The R-peaks are used to localize heartbeats and the intervals between R-peaks (R–R intervals) are used to measure IBIs (Clifford et al 2006) . Alternative detection methods have been developed, which may or may not require body contact. Among the contactless methods are millimeter wave radar (Chen et al 2022) and remote photoplethysmography (PPG) (de Haan and Jeanne 2013). An often used method that requires body contact is PPG, where the blood volume pulse is measured by means of a sensor on the finger or on the wrist (Hoog Antink et al 2021).

Another approach is ballistocardiography, which measures cardiac vibrations through the structure upon which the body resides. The vibrations can be measured using sensors installed on the bed, either on top of or below the mattress. In a comparable approach, seismocardiography (SCG), cardiac vibrations are measured on the chest wall. SCG has long been studied (Agress and Fields 1959) and has regained attention since micro-electromechanical systems such as accelerometers, have become widely available at low cost (Taebi et al 2019, Rai et al 2021). SCG measurement may performed by means of an accelerometer inside a device worn on the chest, as an adhesive patch or through a band. The advantages of using a chest-worn accelerometer are multiple; modern accelerometers have a small footprint, are relatively low in cost, have low power consumption, and only require mechanical contact with the chest at a single point. There is no need for skin contact, neither electrical nor optical, which means that biocompatibility issues leading to allergy, or skin irritations can be avoided, and skin preparation is not required. These properties make it suitable for prolonged wearing, or nightly re-application.

Many studies use SCG to investigate cardiac functioning, specifically cardiac time intervals, such as the pre-ejection period and the left ventricular ejection time (Khosrow-Khavar et al 2017, Dehkordi et al 2019). Some studies also look at the IBIs (Wahlström et al 2017, Cocconcelli et al 2020, Centracchio et al 2023, Parlato et al 2023, Milena et al 2023). In one study, 66 participants equipped with an array of 32 IMUs, were measured in supine posture for one minute, and IBIs were computed by means of a Hidden Markov Model (Wahlström et al 2017). Other studies used templates. In one study, 13 participants wore an accelerometer and rested in a chair for four minutes and templates were automatically calibrated for each participant (Cocconcelli et al 2020). In another study, 100 participants wore an accelerometer and a gyroscope and were measured in supine position for 6 min and 48 s in average, after which participant-specific templates were automatically matched (Centracchio et al 2023, Parlato et al 2023). In yet another study, 21 subjects equipped with an IMU were measured for two minutes in supine posture, and heart beats were detected by envelope detection using a Hilbert transformation (Milena et al 2023).

A common characteristic of all these studies is that the recording times were relatively short and that subjects followed a protocol in which they were asked to not move and remain lying in the same posture. However, during protocol-free conditions, such as sleep, subjects are free to move their body and limbs, and take any posture they prefer. A method for estimating the heart rhythm should therefore be robust against these conditions and be accurate at the same time.

Freedom of movement brings additional challenges to the processing of the data. Global template approaches, for example, cannot be used, as the shape of the waveform may change with posture. This problem has been addressed in the past by exploiting local periodicity, rather than global templates (Brüser et al 2013). Local periodicity may be viewed as maximally localized templates, because two neighbor heartbeats are each other's templates. This method was evaluated on ballistocardiography signals, from which it was able to accurately estimate IBIs, but the accuracy of beat localization remained unclear. From literature we know, that, at least in supine posture, the SCG signal is locally periodic (Dehkordi et al 2019), which is why a local periodicity approach may be suitable for SCG as well.

Here, we evaluate the detection of heartbeats from a chest-worn accelerometer. We evaluate a previously suggested approach based on local periodicity (Brüser et al 2013) and we introduce a new method, which extends the previous method with maximum a posteriori estimation (MAP). We compare the two methods on beat localization, IBI estimation, and on the derived metrics IHR and HR. We use overnight sleep recordings from a population that is diverse in gender, body-mass index, age, with various sleep disorders.

2.1. Data collection

Data were collected at Sleep Medicine Center Kempenhaeghe, Heeze, the Netherlands, as part of the SOMNIA project (van Gilst et al 2019). Recruitment was done among individuals scheduled for overnight polysomnography as part of the standard diagnostic process. Participants undergoing CPAP titration, with intellectual disabilities, or with an earlier diagnosis of atrial fibrillation or other rhythm disorders were excluded. Next to the polysomnographic montage, the participants wore a prototype chest-worn device including a 3D accelerometer (ADXL355, Analog Devices, Wilmington MA, USA), see figure 1. The device had a cable that was connected to a remote recorder, and staff were instructed to mount the sensor device on the thoracic respiratory belt, such that the cable was at the top right. The raw accelerometer signals (${\boldsymbol{x}},$ ${\boldsymbol{y}},$ and ${\boldsymbol{z}},$ each sampled at 250 Hz) were recorded in the recorder and were later downloaded into a computer for analysis. ECG was recorded by the PSG system, using a modified lead II montage sampled at 512 Hz.

Zoom In Zoom Out Reset image size

The SOMNIA study was reviewed by the medical ethical committee of the Maxima Medical Center (Eindhoven, the Netherlands. File no: N16.074). All participants provided consent prior to study initiation. The protocol for data analysis was approved by the Medical Ethical Committee of the Kempenhaeghe hospital and by the Internal Committee of Biomedical Experiments of Philips Research. A total of 147 participants was included (69 female, 78 male), with ages ranging from 17 to 80 years (mean 46.8, SD 15.9). BMI ranged from 18 to 42 kg m⁻² (mean 26, SD 4). The three sleep disorders with the highest prevalence among the included individuals were Sleep Disordered Breathing (71), Insomnia (58), and Movement Disorders (31). Recording durations ranged from 6.7 to 9.8 h (mean 8.3, SD 0.5).

2.2. Data preprocessing and transformation

Prior to detecting heartbeats, we preprocess the data and transform it into spectral and temporal representations, as is shown in figure 2. First, the three components of the acceleration signal $a$ are summed and bandpass filtered with cutoff points of 3 Hz and 25 Hz, such that a one-dimensional SCG signal $s$ is obtained. The 3 Hz cutoff point is chosen to suppress lower-frequency components such as respiration. The summation reduces the sensitivity for changes in the orientation of the sensor with respect to the heart. Even though the instructions in the study prescribe a fixed orientation, in applications in which the sensor does not have a cable, users may mount the sensor in a different orientation. Especially in free-living applications where users re-apply the sensor nightly, there can be quite some variation. Furthermore, the orientation with respect to the heart may vary within a night due to posture changes and arm movements. Finally, the one-dimensional SCG signal s is used for the spectral and temporal transformations.

Zoom In Zoom Out Reset image size

Next to the transformations, the SCG signal $s$ is partitioned into segments. A segment is an interval defined by consecutive samples for which $\left|s\right|\leqslant {s}_{{th}},$ with ${s}_{{th}}$ defining a threshold (empirically set to 0.06 g), such that large motion artifacts and posture changes are excluded from segments. Heartbeat detection is only performed within segments. In time intervals external to segments, no heartbeat detection is performed, as motion artifacts mask any heartbeats. An example of partitioning into segments may be found in figure 5(a).

The spectral representation ${D}_{S}$ has been introduced in Brüser et al (2013) and comprises local periodicity spectra:

$$\begin{eqnarray}&&{D}_{S}={\{S}_{1}\left({\tau }^{{\prime} }\right),{S}_{2}\left({\tau }^{{\prime} }\right),\ldots \},\end{eqnarray} \tag{ 1 }$$

with spectrum ${S}_{j}\left({\tau }^{{\prime} }\right)$ representing local periodicity as a function of ${\tau }^{{\prime} },$ around a center time ${t}_{m}^{j}=\left(j-1\right){T},$ with $T$ = 200 ms. The local periodicity expresses the similarity of two intervals from $s$ as a function of the period time ${\tau }^{{\prime} },$ with one interval ranging from ${t}_{m}^{j}-{\tau }^{{\prime} }$ to $t,$ and one interval ranging from ${t}_{m}^{j}$ to ${t}_{m}^{j}+{\tau }^{{\prime} },$ see figure 3. Due to the prominence of the heartbeats in $s,$ high periodicity typically occurs when ${\tau }^{{\prime} }$ is equal to the time distance between the last heartbeat before ${t}_{m}^{j}$ and the first heartbeat after ${t}_{m}^{j}$ (which together constitute an IBI). As a typical IBI is much longer than $T,$ multiple spectra will have their peak at the same ${\tau }^{{\prime} }.$ When the periodicity over time is plotted on a gray scale, as is shown in figure 4(c), IBIs will become visible as horizontal bars, with height and duration equal the IBI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The temporal representation ${D}_{T}$ expresses the locally normalized magnitude of $s,$ with separate components for the positive and negative parts, see figure 4(d). The normalization for each sample $n$ is relative to the maximum of the product of the absolute value of $s$ and a Gaussian window function centered at the sample $n,$ with a standard deviation of 1.7 s (corresponding to an IBI at a heartrate of 35 beats per minute). If we denominate the normalized positive component ${f}_{p}$ and the normalized negative component ${f}_{n},$ then we get:

$$\begin{eqnarray}&&{D}_{T}={(f}_{p}\left(t\right),{f}_{n}\left(t\right))\end{eqnarray} \tag{ 2 }$$

2.3. MAP estimation

Core to our method is a latent-variable model. This describes the stochastic relation between observed data ${D}_{S}$ and ${D}_{T},$ and a latent random variable $\theta ={(t}_{1},{\boldsymbol{\tau }}$), describing the heartbeats with the starting time ${t}_{1}$ and the IBIs ${\boldsymbol{\tau }}=\left({\tau }_{1},{\tau }_{2},\ldots \right).$ The starting time ${t}_{1}$ is the localization time of the first beat and the IBIs are equal to the difference between the localization times of beats $i$ and $i+1:$ ${\tau }_{i}={t}_{i+1}-{t}_{i}.$ Given knowledge of ${t}_{1},$ the most probable value of ${\boldsymbol{\tau }}$ can be determined by means of MAP estimation:

$$\begin{eqnarray}&&\hat{{\boldsymbol{\tau }}}=\mathrm{argma}{x}_{{\boldsymbol{\tau }}}\,P\left({\boldsymbol{\tau }}\,\right|\,{D}_{S},{D}_{T},{t}_{1}),\end{eqnarray} \tag{ 3 }$$

which can be rewritten by using Bayes' rule:

$$\begin{eqnarray}&&\hat{{\boldsymbol{\tau }}}={\mathrm{argmax}}_{{\boldsymbol{\tau }}}\,P\left({D}_{S},{D}_{T}\,\right|{\,t}_{1},{\boldsymbol{\tau }})P\left({\boldsymbol{\tau }}\,\right|\,{t}_{1}),\end{eqnarray} \tag{ 4 }$$

with likelihood ${\rm{P}}\left({D}_{S},{D}_{T}\,\right|{{t}}_{1},{\boldsymbol{\tau }})$ and prior $P\left({\boldsymbol{\tau }}\,\right|\,{t}_{1}).$

2.3.1. Computing the likelihood

First, we factorize the likelihood per IBI $i:$

$$\begin{eqnarray}&&P\left({D}_{S},{D}_{T}\right|{\boldsymbol{\tau }})\cong \displaystyle \prod _{i}P\left({D}_{S}^{i},{D}_{T}^{i}\right|{t}_{i},{\tau }_{i}),\end{eqnarray} \tag{ 5 }$$

with ${D}_{S}^{i}$ and ${D}_{T}^{i}$ the parts of the data in a time interval that is spanned by ${\tau }_{i}$ (from ${t}_{i}$ to ${t}_{i}+{\tau }_{i}$). Generally, ${D}_{S}^{i},{D}_{T}^{i}$ depend on ${\tau }_{j}$ iff $i=j.$ However, as the last beat of ${\tau }_{j}$ is the same as the first beat of ${\tau }_{j+1},$ the factorization is a (good) approximation. Next, we model the likelihood for each IBI as the product of the spectral likelihood $P\left({D}_{S}^{i}|{{{\rm{t}}}_{{\rm{i}}},\tau }_{i}\right)$ and the temporal likelihood $P\left({D}_{T}^{i}|{{{\rm{t}}}_{{\rm{i}}},\tau }_{i}\right):$

$$\begin{eqnarray}&&P\left({D}_{S}^{i},{D}_{T}^{i}|{{t}_{i},\tau }_{i}\right)\cong P\left({D}_{S}^{i}|{{t}_{i},\tau }_{i}\right)\,P\left({D}_{T}^{i}|{{t}_{i},\tau }_{i}\right).\end{eqnarray} \tag{ 6 }$$

The temporal likelihood is high for data of which the acceleration has a high magnitude at the occurrence of a heartbeat. Therefore, we interpret the locally normalized magnitude (see equation (2)) as the (unnormalized) temporal likelihood:

$$\begin{eqnarray}&&P\left({D}_{T}^{i}|{{t}_{i},\tau }_{i}\right)={c}_{T}{f}_{a}({t}_{i}),\end{eqnarray} \tag{ 7 }$$

with ${c}_{T}$ a normalizing constant, ${t}_{i}$ the localization time of heartbeat $i,$ and ${f}_{a}$ either ${f}_{p}$ or ${f}_{n},$ depending on whether we choose to localize beats at the positive or negative component of $s.$ The choice for either ${f}_{p}$ or ${f}_{n}$ will be described in section C3. Figure 4(d) shows temporally transformed data ${D}_{T}$ of high likelihood, as the peaks in ${f}_{n}(t)$ coincide with the heartbeats, that are shown as dashed lines in figure 4(e).

The spectral likelihood is high for data of which the spectral peaks match with the IBIs, as is the case in figure 4(c), where the dark horizontal bars (the spectral peaks) match in height and duration with the IBIs in figure 4(e). To compute the spectral likelihood, we define ${D}_{S}^{i}$ as a subset of spectra with their center times within the interval spanned by IBI $i:$

$$\begin{eqnarray}&&{D}_{S}^{i}=\{{S}_{n\left(i\right)}\left({\tau }^{{\prime} }\right),\ldots ,{S}_{n\left(i\right)+N\left(i\right)}\left({\tau }^{{\prime} }\right)\},\end{eqnarray} \tag{ 8 }$$

with $n\left(i\right)$ the index of the first spectrum with ${t}_{m}^{n(i)}\geqslant {t}_{i},$ and $n(i)+N\left(i\right)$ the index of the last spectrum with ${t}_{m}^{n\left(i\right)+N(i)}\lt {t}_{i}+{\tau }_{i}.$ We interpret the geometric mean of the spectral values in ${D}_{S}^{i}$ as the (unnormalized) spectral likelihood:

$$\begin{eqnarray}&&P\left({D}_{S}^{i}|{{t}_{i},\tau }_{i}\right)={c}_{S}\sqrt[\left.N(i\right)]{{\prod }_{\left.j=n(i\right)}^{n\left(i\right)+N\left(i\right)}{S}_{j}\left({\tau }_{i}\right)},\end{eqnarray} \tag{ 9 }$$

with ${{\rm{c}}}_{S}$ a normalizing constant. $P\left({D}_{S}^{i}|{{{\rm{t}}}_{{\rm{i}}},\tau }_{i}\right)$ will be high if all ${S}_{j}({\tau }^{{\prime} })$ in ${D}_{T}^{i}$ have a high peak at ${\tau }^{{\prime} }={\tau }_{i}.$

2.3.2. Computing the prior

The prior $P\left({\boldsymbol{\tau }}\right|{t}_{1})$ in equation (4) is computed for each recording individually, by pre-inspecting ${D}_{S}.$ First, for each of the spectra ${S}_{j}\left({\tau }^{{\prime} }\right),$ the values ${v}_{j}$ of the highest peaks (${{{v}_{j}=\max }_{{\tau }^{{\prime} }}\{S}_{j}\left({\tau }^{{\prime} }\right)\}$) and the period times ${p}_{j}$ of the highest peaks (${p}_{j}={{\mathrm{argmax}}_{{\tau }^{^{\prime} }}\{S}_{j}\left({\tau }^{{\prime} }\right)\}$) are computed. Next, the period times ${p}_{j}$ of high peaks (${v}_{j}\geqslant {v}_{{th}},$ with ${v}_{{th}}$ experimentally set to 0.70), are put in a sequence, which is then median filtered (kernel size = 151) to remove outliers. The median-filtered period times, ${p}_{j}^{{\prime} },$ are put in a new sequence (at their original positions $j$), and the period times for the positions $k$ for which no sufficiently strong peak was found (${v}_{k}\lt {v}_{{th}}$), are linearly interpolated. The resulting sequence $({p}_{1}^{{\prime} },{p}_{2}^{{\prime} },\ldots )$ gives a period time for each time point ${t}_{m}^{j},$ which are used for the computation of the prior. The prior is then computed as follows:

$$\begin{eqnarray}&&P\left({\boldsymbol{\tau }}\right|{t}_{1})=\displaystyle \prod _{i}P({{t}_{i},\tau }_{i}),\end{eqnarray} \tag{ 10 }$$

in which each factor $P({\tau }_{i})$ follows a Gaussian distribution with mean ${\mu }_{i}$ and standard deviation $\sigma :$

$$\begin{eqnarray}&&P\left({{t}_{i},\tau }_{i}\right)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{{\tau }_{i}-{\mu }_{i}}{\sigma }\right)}^{2}}.\end{eqnarray} \tag{ 11 }$$

The mean ${\mu }_{i}$ is the geometric mean over the periods ${p}_{j}^{{\prime} }$ with $j$ from $n\left(i\right)$ to $n\left(i\right)+N\left(i\right),$ as in equation (9), and $\sigma$ is the standard deviation of the period times ${p}_{j}.$

2.3.3. Finding the maximum

With (6), (7), and (11), the expression for the MAP estimate $\hat{{\boldsymbol{\tau }}}$ (equation (5)) can be rewritten as follows:

$$\begin{eqnarray}&&\hat{{\boldsymbol{\tau }}}={{\rm{argmax}}}_{{\boldsymbol{\tau }}}\displaystyle \prod _{i}P\left({D}_{T}^{i}\,|{{t}_{i},\tau }_{i}\right)\,P\left({D}_{S}^{i}\,|{{t}_{i},\tau }_{i}\right)\,P\left({{t}_{i},\tau }_{i}\right),\end{eqnarray} \tag{ 12 }$$

where the expressions for the factors in the right-hand side are given by equations (7), (9), and (11).

As there may be many heartbeats in one recording (typically ~30.000), finding the maximum becomes intractable. Fortunately, the problem can be reformulated into a Markov Decision Process, which breaks down the problem into steps (see supplemental materials). In each step, a current localization time ${t}_{i}$ is advanced to a next localization time ${t}_{i+1},$ with ${t}_{i+1}={t}_{i}+{\tau }_{i},$ by deciding for a value of ${\tau }_{i}.$ In each step, a reward is obtained, equal to the $\mathrm{log}$ of the corresponding factor in (12):

$$\begin{eqnarray}&&R\left({t}_{i},{\tau }_{i}\right)=\mathrm{log}{\rm{}}\{P\left({D}_{T}^{i},|,{{t}_{i},\tau }_{i}\right)\,P\left({D}_{S}^{i},|,{{t}_{i},\tau }_{i}\right)\,P\left({{t}_{i},\tau }_{i}\right)\}.\end{eqnarray} \tag{ 13 }$$

The value $V({t}_{i})$ of a localization time ${t}_{i}$ is the maximum of the sum of rewards that may be obtained from ${t}_{i}$ onwards:

$$\begin{eqnarray}&&V({t}_{i})={\max }_{{\tau }_{i}}\{R({t}_{i},{\tau }_{i})+V({t}_{i}+{\tau }_{i})\},\end{eqnarray} \tag{ 14 }$$

which is known as the Bellman equation (Bellman 1952) and provides guidance on decisions for ${\tau }_{i}$ such that an optimal solution for equation (12) is found. During each step, the process aims to decide for a value of ${\tau }_{i}$ that maximizes $V({t}_{i}),$ for which it inspects $C$ IBI candidates, being the ${\tau }_{i}$ for the $C$ highest peaks in $R({t}_{i},{\tau }_{i}).$ For each IBI candidate, equation (14) is evaluated recursively, up to a finite horizon, or recursion depth, $H.$ The IBI candidate ${\tau }_{i}$ for which $R({t}_{i},{\tau }_{i})+V({t}_{i}+{\tau }_{i})$ is maximal is then selected, after which the localization point is advanced to ${t}_{i+1},$ until all data have been processed. We chose $C=3$ and $H=8.$

The decision process needs a starting point ${t}_{1}.$ To this end, the data are partitioned into segments, in which a segment is an interval in which $\left|s\right|\leqslant {s}_{{th}},$ with $s$ the SCG signal and ${s}_{{th}}$ a threshold (empirically set to 0.06 g). We only process segments with a duration $\geqslant$ 15 s as in shorter segments the impact of body and sensor settling may be too prominent. For each segment, a starting point ${t}_{1}$ is searched such that $R({t}_{1},{\tau }_{1})$ is maximal. This is done for both polarities of the SCG signal, so with ${f}_{a}={f}_{p}$ and ${f}_{a}={f}_{n},$ see (7). The polarity ($p$ or $n$) with the highest reward is then chosen and the decision process commences from ${t}_{1}$ in forward direction as described above, and in reverse direction, in a similar manner, but then with ${t}_{i+1}={t}_{i}-{\tau }_{i}.$

2.4. Evaluation

2.4.1. Operating point selection

After estimating the IBIs and localization times, the rewards $R({t}_{i},{\tau }_{i})$ for each of the beats $({t}_{i},{\tau }_{i})$ are compared against a threshold ${R}_{{th}}$ and beats that have their reward below the threshold are discarded. The thresholding enables the selection of an operating point; decreasing ${R}_{{th}}$ decreases the number of detected beats and increases the mean posterior probability of the remaining beats.

2.4.2. Synchronization and pairing with ECG beats

We compute the localization times for the ECG signal, through the detection of R-peaks, using a nonlinear transformation in combination with first-order Gaussian differentiation (Kathirvel et al 2011). The IBIs are then computed as the time-differences between the localization times and IBIs that exceed a threshold ${T}_{{th}}$ are discarded (we chose ${T}_{{th}}$ = 1714 ms, corresponding to 35 beats per minute).

Using the localization times of the SCG and ECG beats, we synchronize the two acquisitions, globally and locally. The global synchronization aims at reducing overall clock offset and drift and was performed by minimizing the localization error between all matched ECG and SCG beats within each recording. An ECG beat and SCG beat become a matched pair if they are each other's nearest neighbors in time, and their distance is less than 250 ms. The localization error, with an example shown in figure 5, may be positive or negative, and may exhibit jumps (caused by changes in the SCG signal morphology) and trends (caused mainly by clock drift). Using a manual procedure, the offset and drift of the SCG signal are adapted such that the localization error trends disappeared, and the mean localization error becomes zero (globally).

Zoom In Zoom Out Reset image size

In addition to the global synchronization, we also performed a local optimization aimed at estimating a lower limit for the localization error. For each segment within a recording, the timing offset of the SCG signal is adjusted starting from the global clock offset, until the mean localization error within the segment becomes zero. No drift adjustment is performed for the local synchronization.

2.4.3. Performance metrics and statistics

For performance on beat detection, localization, and IBI estimation, we synchronize over segments without body motion as well as over full recordings. For both cases, we compute the same performance metrics. First, we label all beats; ECG and SCG beats that are in a pair get the label true positive (TP), ECG beats that are not in a pair get the label false negative (FN), and SCG beats that are not in a pair get the label False Positive (FP). From these labels we compute the metrics sensitivity = TP/(TP + FN) and positive predictive value (PPV) = TP/(TP + FP). To get meaningful Sensitivity and PPV metrics, in case of synchronization over segments, we only consider segments with a duration of at least 5 min. In case of synchronization over full recordings, all processed segments are included and the Sensitivity and PPV are computed over the full recording.

For the TP beats, we further define the IBI error as the difference between the SCG and the ECG IBIs and the localization error as the difference between the SCG and ECG localization times. Finally, we perform factor analysis for the RMSE IBI error. For the numeric factors (BMI and Age), we perform Spearman rank correlation, while we perform a Mann–Whitney U test for the categorical factors (Sex and presence of the three most prevalent disorders in the study).

For performance of HR and IHR, we synchronize over full recordings. For the IHR, for both SCG and ECG, a stepwise constant function is computed, with levels equal to the inverse of the IBIs and level changes at the localization times. The two functions are then resampled to 10 Hz, after which the error is computed as the difference between the two, thereby discarding gaps. For the HR, for both SCG and ECG, we take 15 s intervals over which we compute the HR as the inverse of the average of the IBIs, thereby discarding intervals with less than 5 IBIs. The error is then computed as the difference between the SCG and ECG HR.

As our method comprises the extension of an existing method (Brüser et al 2013) with MAP, we also compute the results for the existing method (without MAP). Furthermore, we tune the operating point of the method with MAP, such that its Sensitivity is close to, but not lower than that of the method without MAP. We compute the statistics for both methods and compare them by means of a Wilcoxon signed-rank test. As we have more than one statistic on the same data, we use a Bonferroni-corrected significance level of $\alpha /m$ ($\alpha \,=$ 0.001, $m\,=$ 9) to prevent false discoveries.

The overnight recordings of 147 participants were preprocessed and segmented, and the median coverage was 79.9% (Q1 = 65.7%, Q3 = 86.2%). This coverage is the combined result of segmentation, with median coverage of 96.6% (Q1 = 94.7%, Q3 = 97.6%), and thresholding within each segment, as described in the section on operating point selection. We obtained performance metrics for the methods with and without MAP, see table 1. All errors were significantly lower with MAP than without. The median RMS IBI error decreased by a factor 2.9 from 25.2 to 8.8 ms (the median IBI was 938 ms). The median RMS localization error decreased by a factor 1.3 from 73.2 to 58.1 ms. Although significant, this reduction factor (1.3) is rather limited when compared to the reduction factor for the RMS IBI error (2.9).

Figure 5 shows how the offset in the localization error can change within one recording and explains the limited decrease of the localization error in case of synchronization over full recordings. Figure 6 shows two short intervals, one from a segment before a posture change and another from a segment after it. In the first segment, the method localized on an earlier positive peak in the SCG signal, while it localized to a later negative peak in the second. The difference in localization caused a change in the offset of the localization error, which could not be compensated for when synchronizing over full recordings. Further examples are given in the supplemental materials.

Zoom In Zoom Out Reset image size

When synchronizing over segments, on the other hand, the offset in the localization error could be compensated for, and the median RMS localization error decreased from 58.1 ms (as in table 1, with MAP) to 7.90 ms (as in table 2, with MAP). This decrease shows that the higher localization error in table 1 was only due to changes in offset and that the MAP method consistently localizes at the same point (phase) in the cardiac cycle. Again, all errors were significantly lower with MAP than without. The IBI localization error decreased by a factor 2.8 from 15.8 to 5.6 ms, while the RMS localization error decreased by a factor 5.6 from 44.6 to 7.9 ms. As only segments longer than 5 min were included in table 2, the coverage was 77.9% (Q1 = 67.0%, Q3 = 84.9%).

The IBI error for the method with MAP, when aggregated over all recordings (N = 3.5 × 10⁶ beats), had 95% limits of agreement of −10.6 and 11.1 ms, while the quantiles Q1, Q2, and Q3 were −2.8, 0.0, and 2.8 ms respectively. Mean coverage was 74.6%, over 1236.8 recorded hours.

We investigated the influence of various factors on the performance. The RMS IBI error increased (p < 0.001) due to the presence of sleep disordered breathing (median from 5.2 to 6.2 ms, +19%), movement disorders (median from 5.4 to 6.4 ms, +19%), and sex (median in males 5.4 ms, median in females 5.9 ms, +9%). Furthermore, positive rank correlations for the RMS IBI error (p < 0.001) were found for BMI (0.12) and age (0.22). These numbers relate to the method with MAP in combination with synchronization over segments.

We computed the derived metrics IHR and HR, of which Bland–Altman density plots are shown in figure 7. For the method with MAP, the median error and 95% limits of agreement were 0.00, −1.20 and 1.25 bpm for the IHR and 0.00, −1.65 and 1.03 bpm for the HR. If we define a measurement as being accurate if it is within a range around the reference of either ±5 bpm or ±10% (whichever is larger), then 98.8% of the IHR measurements and 100.0% of the HR measurements were accurate. These percentages were obtained at median coverages of 80.2% and 79.7% respectively (see tables SI and SII of the supplemental materials). The Spearman correlation coefficient was 0.992 for the IHR (p < 0.001) and 0.997 for the HR (p < 0.001). Further metrics and a comparison between the methods with MAP and without MAP may be found in tables SI and SII in the supplemental materials.

Zoom In Zoom Out Reset image size

We evaluated heartbeat detection from a chest-worn accelerometer in a set of 147 overnight recordings among a diverse population with a variety of sleep disorders. We determined performance for the primary metrics beat detection, beat localization, and IBI estimation, and for the derived metrics HR and IHR. IBI estimation was accurate with an MAE of 3.5 ms, which is less than the sampling period of the SCG signal (4.0 ms) and less than 0.5% of the median IBI time (938 ms). For the derived metrics IHR and HR, we found that 98.8% and 100.0% of all measurements were accurate according to the American National Standard 'Cardiac monitors, heart rate meters, and alarms' (Association for the Advancement of Medical Instrumentation 2002). It should be noted that these percentages were obtained over the covered parts of the recordings (median around 80%).

We have found that adding MAP to the detection method significantly improved the performance. Without reducing sensitivity, the RMS IBI error decreased by a factor 2.8, while the RMS localization error decreased by a factor of 5.6. A possible reason for the improvement is the larger scope of the algorithm; when estimating one IBI (in the Markov Decision Process), all data up to 8 IBIs ahead contribute to the choice for that single IBI.

Table 3 compares our results with other studies on seismocardiography, which reported IBI performance. It should be emphasized that a direct benchmark is not possible, since studies used different protocols, and ultimately, different measuring setups, accelerometers and even performance metrics, so this serves merely as a high-level comparison. Our study shows that it is possible to achieve a high level of accuracy in IBI, HR and IHR estimation also on longer recordings, including changes of posture and body movements, achieving errors at least comparable with the best studies reported even if those were obtained in shorter recordings with subjects laying still in fixed postures (mostly supine).

Table 3. Comparison of results against other studies.

	Year	Study	Sensor	N	Duration	Posture	Results	Method
A	2017	(Wahlström et al 2017)	Accelerometer	66	1 min	Supine	IBI MAE = 6.10 $+/-$ 31.89 ms	Hidden Markov Model
B	2019	(Kathirvel et al 2011)	IMU	25	7 min	Supine	HR LoA = −5.06 to 5.66 bpm	Autocorrelated Differential Algorithm
C	2020	(Cocconcelli et al 2020)	Accelerometer	20	50 min	Supine	IBI RMSE = 4.2 ms	Automatic template calibration and matching
				13	4 min	Resting in chair	IBI RMSE = 6.2 ms
D	2021	(Association for the Advancement of Medical Instrumentation 2002)	Gyroscope at head	19	5 min	Supine	HR MAE ≈ 2.5 bpm	Envelope based (Hilbert transform)
E	2021	(Schipper et al 2021)	Accelerometer	30	40 min	All	HR LoA = −19.3 to 17.8 bpm	Envelope based
F	2023	(Centracchio et al 2023)	IMU	21	2 min	Supine	LoA = −26.43 to 26.43 ms	Envelope based (Hilbert transform)
G	2023	This study	Accelerometer	147	8.3 h	All	IBI RMSE = 5.6 ms
							IBI MAE = 3.5 ms
							LoA:
							- IBI: −10.6 to 11.1 ms	Maximum A posterior estimation
							- HR: −1.6 to 1.0 bpm
							- IHR: −1.2 to 1.3 bpm

As for the localization error, we found that it was substantially higher when the synchronization was performed over whole recordings (MAE = 45.9 ms, table 1), than over segments without body movements (MAE = 4.7 ms, table 2). The body movements, and perhaps more importantly, posture changes, cause the MAP search algorithm to use different localization points in the cardiac cycle. These cause a change in the offset in the localization error, which is removed by synchronization within a stable posture period. In these intervals, the algorithm consistently assigns localization points to the same point in the cardiac cycle, achieving a relatively small localization error. Note that these intervals were detected by the algorithm itself from the acceleration data, without manual annotations of any kind.

The influence of demographics and sleep disorder prevalence on the performance was limited. Spearman correlation coefficients for age and BMI were only 0.22 and 0.12, while the largest increase in RMS IBI error was due to the occurrence of sleep disordered breathing (from 5.2 to 6.2 ms, +19%), such that still a high level could be maintained. Only a small increase (+9%) in IBI error was found for female over male participants. This shows that the method is robust against body movements and posture changes but also against the presence of sleep disorders and different demographic factors, such as gender.

As stated, the MAP estimation was beneficial for the robustness and accuracy of the method. The probabilistic approach with the latent variable model, Bayes' rule, and the MAP estimation, provided a structured framework for estimation. Separate models could be created for the likelihood and the prior, and the product of the two automatically yielded an objective function. The chain rule of probability allowed for the reformulation of the MAP problem into a Markov Decision Process problem, of which an optimum solution is described by the Bellman equation. The optimum can be approximated by a tree-search algorithm, of which the width and depth can be adjusted to trade accuracy against run time.

A limitation of this study is that subjects with prior diagnosis of cardiac rhythm disorders, most notably atrial fibrillation, were excluded from the study. Especially during sleep monitoring, atrial fibrillation is highly prevalent among patients with sleep-disordered breathing (Linz et al 2022). During cardiac arrhythmias, local periodicity, upon which the method relies, may be absent and therefore results may be worse. Future work should include patients with cardiac arrhythmias and annotation of cardiac rhythms should be made, such that performance under such conditions can be investigated. Another limitation is that the computation of the prior was done by inspection of the full recording; we did not investigate the effect of a constrained look-ahead, which is required for low-latency applications. Furthermore, the Markov Decision Process may be improved; instead of exhaustive search, dynamic programming or reinforcement learning may be applied. Finally, further research may be conducted in the mapping of the three accelerometer axes into the one-dimensional SCG signal. The summation that we used in the present study, reduces sensitivity to changes in the orientation of the sensor with respect to the heart and maintains periodicity, as any linear combination of signals with period T is also periodic with T. However, noise may not be uniformly distributed over the axes, so alternative methods may be evaluated in future work, such as a vector sum which may increase SNR while preserving amplitude periodicity.

The possibility to accurately and robustly estimate HR and IHR, as presented in this paper, together with our earlier published work on respiration (Schipper et al 2021, 2023), may enable cardio-respiratory applications to work with a single sensor. The sensor can be small and only needs mechanical attachment at a single point on the chest, without the need for skin contact. One example of such an application is the detection of patient deterioration using respiration rate and HR (Jacobs et al 2021). Another application may be in positional therapy devices for obstructive sleep apnea, which use a chest worn tactile feedback system to entrain the patient to avoid supine sleeping, thereby reducing the apnea burden (Berry et al 2019). The addition of our method may enable the estimation of residual apnea burden, quantified in the form of the apnea-hypopnea index by using IBIs and respiratory effort from an accelerometer in the therapy device itself (Xie et al 2023).

In conclusion, we propose a method that detects heartbeats from a chest-worn accelerometer, demonstrating robustness and high accuracy during overnight recordings with unrestricted postures and body movements, in a diverse population and in the presence of disorders. This enables applications, that directly measure respiration, to monitor cardiac activity robustly and accurately without the need for an additional sensor.

This work has been performed in the IMPULS framework of the Eindhoven MedTech Innovation Center (e/MTIC, incorporating Eindhoven University of Technology, Philips Research, and Sleep Medicine Center Kempenhaeghe). The funders had no role in the study design, decision to publish, or preparation of the manuscript. Fons Schipper, Angela Grassi, Pedro Fonseca, and Jan Brouwer are employed by Philips Research. The employer had no influence on the study and on the decision to publish. Ruud van Sloun is employed by both Philips Research and by the Eindhoven University of Technology. The other authors declare no competing interests.

The authors would like to thank Stan Hullegie, Bertram Hoondert (Kempenhaeghe), Ad Rommers, and Henry van Vugt (Philips) for their help in data acquisition.

The data cannot be made publicly available upon publication because they are owned by a third party and the terms of use prevent public distribution.