A comparison of three heart rate detection algorithms over ballistocardiogram signals
AA comparison of three heart rate detection algorithms over ballistocardiogram signals
Ibrahim Sadek , Bessam Abdulrazak AMI-Lab, Computer Science Department, Faculty of Science, University of Sherbrooke, Sherbrooke, QC, Canada
Abstract
Heart rate (HR) detection from ballistocardiogram (BCG) signals is challenging for various reasons. For example, BCG signals’ morphology can vary between and within-subjects. Also, it differs from one sensor to another. Hence, it is essential to evaluate HR detection algorithms across different datasets and under different experimental setups. This paper investigated the suitability of three algorithms (i.e., MODWT-MRA, CWT, and template matching) for HR detection across three independent BCG datasets. The first two datasets (Datset1 and DataSet2) were obtained using a microbend fiber optic (MFOS) sensor, while the last one (DataSet3) was obtained using a fiber Bragg grating (FBG) sensor. DataSet1 was collected from 10 OSA patients during an in-lab PSG study, Datset2 was obtained from 50 subjects in a sitting position, and DataSet3 was gathered from 10 subjects in a sleeping position. The CWT with derivative of Gaussian (Gaus2) provided superior results than the MODWT-MAR, CWT (frequency B-spline — Fbsp-2-1-1), and CWT (Shannon — Shan1.5-1.0) for DataSet1 and DataSet2. That said, a BCG template was constructed from DataSet1. Then, it was applied for HR detection across DataSet2. The template matching method achieved slightly superior results than CWT — Gaus2 for DataSet2. Furthermore, it has proved useful for HR detection across DataSet3 despite that BCG signals were obtained from a different sensor and under different conditions. Overall, the time required to analyze a 30-second BCG signal was in a millisecond resolution for the three proposed methods. The MODWT-MRA had the highest performance, with an average time of 4.92 ms.
Keywords : Ballistocardiography; Mobile health; Homecare; Heart rate; Wavelet transform; Template matching Introduction
Remote monitoring of vital signs, that is, body temperature, heart rate (HR), blood oxygen saturation, respiratory rate, and blood pressure, has attracted public health attention due to rapidly emerging infectious diseases, e.g., coronavirus disease [1]. Changes in vital signs are critical in assessing the severity and prognosis of epidemic diseases. Specifically, these changes are significant signs of deteriorating patient health and thus present an Corresponding author. E-mail address: [email protected] (I. Sadek) pportunity for early detection and intervention. In hospital practice, nursing staff and doctors rely on intermittent vital signs usually measured every 8-h shift. Hence, early deterioration indicators could be overlooked, particularly at night, when deterioration could progress undetected until the next morning [2]. Bed-embedded ballistocardiogram (BCG) sensors have presented encouraging results for detecting particular vital signs, namely HR and respiratory rate [3,4]. Additionally, these sensors have been implemented for diagnosing severe sleep disorders, that is, sleep apnea [5 – ’ s medical conditions. There will always be a compromise between the continuity of data and patient convenience. Although BCG sensors can help alleviate some shortcomings of wearable sensors, they are highly prone to motion artifacts, e.g., body movements. Furthermore, they can only be practical for observing patients in a single bed setting. That said, these sensors are not designed to deliver spot readings to substitute nurse observations. However, they are intended for monitoring trends in vital signs, taking into account their capacity to collect longitudinal data [2]. Various signal processing and machine learning algorithms have been suggested to scrutinize BCG signals, considering the multiple factors that affect the signal quality. These algorithms ’ objective is to automatically identify the “ J ” peak of the “ I-J-K ” complex [3]. Under controlled conditions, if the subject sleeps on the bed without movement, this peak can be detected using a classical peak detector. Nonetheless, this is not conceivable in real-life scenarios. The sensor location is another aspect that has to be underlined when we consider the signal quality. Ideally, the closer the sensor is to the chest and abdomen region, the better the signal quality. Hence, the sensor ’ s desired location is under the upper part of the body, wherein it can be placed under the bed sheet or the mattress. In real-life scenarios, we cannot predict subjects ’ sleep positions, and thus, unless the bed is covered entirely by pressure sensors, the signal quality can be highly degraded. Still, this arrangement will increase the deployment ’ s cost. Besides, a BCG signal structure can vary from one sensor to another and between and within patients. These restrictions should always be considered when a system is designed for analyzing BCG signals. Fast Fourier transform (FFT), Hilbert transform, template matching, autocorrelation, cepstrum, wavelet transform, and empirical mode decomposition, among others, have been implemented for automatic HR detection from BCG signals [3,11,12]. oreover, convolutional neural networks (CNNs) have been employed to segment the “I -J- K” complexes a nd detect HR in BCG signals [13 – Related Work
The wavelet transform (WT) aims at decomposing the signal into “approximation” and “detail” components. Thus, the component ( or the sum of multiple “detail” components ), including the most similar pulse-like peaks, can be adopted to locate the J-peaks of the BCG signal. Zhu et al. [16,17] applied the “ à trous ” WT to raw sensor data acquired from a liquid pressure sensor under a pillow. The signals were obtained from 13 healthy subjects during sleep for about 2-h. Motion artifacts caused by head and body movements were reduced by a threshold method using the raw signals’ standard deviation (S
D). The “Cohen -Daubechies-
Feauveau” th and 5 th “detail” components were realigned in the signal phases, and their amplitudes were summed to estimate the BCG signal. Finally, the J-peaks were detected using a modified Pan-Tompkins algorithm [18] after noise reduction with a soft threshold method. Jin et al. [19] employed a translation-invariant WT based on adaptive threshold wavelet shrinkage for signal denoising. The Symlet wavelet of order 8 (sym8) was chosen because it was closer in shape to the BCG signal. The signal was collected from a healthy subject, but there was no information on the data acquisition process. Then, -peaks were detected using a pseudo-period detection approach by locating the signal’s largest swings [20]. Postolache et al. [21] designed a framework for measuring HR from two EMFi sensors embedded in the seat and backrest of a wheelchair. BCG signals were collected from 8 subjects seated in the chair over 15 min. First, signals were denoised using discrete stationary WT combined with a soft threshold method. Second, the denoised signals were decomposed via a discrete WT — Daubechies (db5) wavelet function was utilized, and the cardiac signal was reconstructed by summing the 8 th , 9 th , and 10 th “detail” components. At last, a time-domain peak detection algorithm was used to detect J-peaks. A similar approach was introduced by Pino et al. [22], in which BCG signals were acquired from two EMFi sensors embedded in the seat and backrest of a typical chair. Raw sensor data were gathered from 19 subjects in a lab for over 1 min. Likewise, 35 subjects were in a hospital waiting area for over 2 min. Daubechies (db6) wavelet function was used for the decomposition, and the cardiac signal was reconstructed by summing the 4 th to 7 th “detail” components . J-peaks were detected using a customized peak detector algorithm. Gilaberte et al. [23] proposed to use CWT to detect HR from subjects standing on a bathroom scale. Six subjects participated in the study, and data were recorded over 10-second in different days and conditions (i.e., before and after meals). Subjects were instructed not to talk or move to eliminate noise. The cardiac signal was located using Daubechies (db10) wavelet function at different scale ranges. The authors suggested that two ranges must be explored in the case of very different HR values. Alvarado-Serrano et al. [24] implemented CWT with B-splines to detect HR using data from subjects in a sitting position. A piezo-electric sensor was fixed to a typical chair seat, and raw sensor data was gathered from 7 subjects for about 100-second. The 5 th scale of the CWT was defined as the optimal scale for HR detection. J-peaks were detected through learning and decision stages. In these stages, several experimental parameters were determined that could limit their use in another dataset. Table 1 presents a summary of wavelet-based methods used in the literature to detect HR from BCG signals. On the other hand, Shin et al. [25] proposed to use a TM approach for BCG beat detection. BCG signals were recorded using three sensors, that is, air-mattress, loadcell, and EMFi-film. A BCG template was constructed for each sensor using ensemble averaging for valid BCG cycles centered at J-peak points. Five records of 30-second were gathered for each sensor, and the matching was conducted using the correlation coefficient function. Paalasmaa et al. [26] presented a method for interbeat intervals from BCG signals acquired with a piezo-electric sensor. A BCG template was created using a clustering-based method. Then, interbeat intervals were detected via a customized cross correlation-based approach. The BCG template was continually being updated based on the detected interbeat intervals. Raw sensor data were recorded overnight from 40 patients in a sleep clinic and 20 subjects at their homes. Nevertheless, only 46 overnight recordings were used in the study. Cathelain et al. [27] ntroduced a similar approach to [25]. However, the matching was achieved using dynamic time wrapping. In this study, a Murata SCA11H BCG sensor was deployed, and data were acquired from 10 healthy subjects over 20 to 50 minutes long naps. The initial BCG template was updated with the newly detected J-peaks to alleviate the BCG signal's shape variabilities. So far, we discussed how the wavelet transform and template matching methods were implemented for HR detection. Next, we provide detailed information about our proposed method. Table 1 Summary of wavelet-based approaches used to detect HR from BCG signals. “ CDF ” : Cohen-Daubechies-Feauveau; “ sym ” : Symlet; “ db ” : Daubechies; “ D ” : detail-component ; “s”: seconds; “min”: minutes. Authors Subjects Sensor Environment Acquisition time Wavelet Wavelet function Cardiac signal Zhu et al. [16,17]
13 Liquid pressure sensor Sleep lab 120 min DWT CDF 9/7 Sum(D4:D5)
Jin et al. [19]
1 N/A Lab N/A DWT sym8 N/A
Postolache et al. [21]
8 EMFi sensor Lab 15 min DWT db5 Sum(D8:D10)
Pino et al. [22]
35 EMFi sensor Lab and hospital 2 min DWT db6 Sum(D4:D7)
Gilaberte et al. [23]
6 Strain gauges Lab 10 s CWT db10 N/A
Alvarado-Serrano et al. [24]
7 Piezo- electric sensor lab 100 s CWT B-splines Scale 5 Methodology
In this section, we describe the proposed framework. First, in subsection 3.1, we demonstrate how we collected the two BCG datasets. Second, we introduce the MODWT, CWT, and TM methods in subsections 3.2, 3.3, and 3.4, respectively. Third, in subsection 3.5, we demonstrate how we detect HR through each method.
Experimental Setup and Data Collection
The first dataset was acquired by a microbend fiber optic sensor (MFOS) from 10 sleep apnea patients. The patients underwent polysomnography (PSG) at the sleep lab of Khoo Teck Puat Hospital, Singapore (elapsed time: hours).
The MFOS was placed under the patient’s mattress in the upper part of the bed.
The PSG electrocardiogram (ECG) signals were used as a gold-standard to assess the proposed HR detection methods. For more details about the MFOS and dataset, readers are referred to [5]. The second dataset was collected in a realistic setting by an MFOS from 50 participants on a massage chair (elapsed time: hours). The MFOS was installed on the chair's headrest, and BCG signals were transmitted wirelessly to a computer by Bluetooth.
The study aimed to evaluate the participants’ stress levels at arious time points. The participants underwent a sequence of stress-induced activities, rest (no-activity), and relaxation therapy [28,29]. The continuity of contact was a significant issue in this study; if the participants had lifted or relocated their heads, we could not have recorded the BCG signals.
We manually discarded participants’ data with artifacts severe enough to degrade BCG signal quality in light of this issue. Therefore, we could only analyze data from 39 participants. ECG signals were simultaneously recorded, and they were used as a reference for HR detection. In both situations, obtained BCG signals were preprocessed to separate motion artifacts and no-activity intervals. The preprocessing step was carried out using a sliding time window of 30-seconds with an overlap of 15-seconds. The SD of each time window was computed. Then, the median absolute deviation (MAD) of the SDs was calculated. Time windows with SD greater than 2 times the MAD were regarded as motion artifacts. Furthermore, time windows with SD less than a fixed value (10 mV) were regarded as no-activity intervals and discarded from further analysis. No-activity implies that no pressure force was applied to the MFOS. The remaining time windows were considered informative signals wherein BCG signals could be extracted. Chebyshev Type I band-pass filter was applied to the informative signals (artifact-free) to obtain BCG signals. The band-pass filter consisted of high- and low-pass filters as follows: 1) second-order Chebyshev high-pass filter with a maximum pass-band ripple of 0.5 dB, and a critical frequency of 2.5Hz was applied, 2) fourth-order Chebyshev low-pass filter with a maximum pass-band ripple of 0.5dB and critical frequency of 5Hz.
Maximal Overlap Discrete Wavelet Transform
Wavelet transform is particularly appropriate for analyzing nonstationary signals, such as in the case of BCG signals - In other words, it is a logical approach for analyzing signals with localized features whose locations are unknown a priori . Fourier transform, on the other hand, performs best for stationary signals. Unlike the DWT, the MODWT skips the downsampling after filtering the signal. The reason is that it gains other features, e.g., invariant to time-shifting, the ability to analyze any time series with arbitrary sample size, and increased resolution at a coarser scale. Besides, it generates a more asymptotically efficient wavelet variance estimator than the DWT [30,31]. The MODWT decomposes a signal into a number of “details” and a single “smooth . ” The “details” describe variations at a particular time scale, whereas the “smooth” describes the low -frequency variations. Given a time series 𝑋 𝑡 of 𝑁 samples, the level 𝐽 MODWT is a transform consisting of 𝐽 + 1 vector, that is, 𝑊̃ , ⋯ , 𝑊̃ 𝐽 and 𝑉̃ 𝐽 . All these vectors have a dimension 𝑁 . The vector 𝑊̃ 𝑗 comprises wavelet coefficients linked to changes on the scale 𝜏 𝑗 = 2 𝑗−1 , whereas the 𝑉̃ 𝐽 comprises the MODWT scaling coefficients linked to averages on the scale 𝜆 𝐽 = 2 𝐽 [32]. The 𝑊̃ 𝑗 and 𝑉̃ 𝑗 can be constructed by filtering 𝑋 𝑡 as follows: 𝑊̃ 𝑗,𝑡 = ∑ ℎ̃ 𝑗,𝑙∘ 𝑋 𝑡−𝑙 𝑚𝑜𝑑 𝑁𝑁−1𝑙=0 , (1) 𝑉̃ 𝑗,𝑡 = ∑ 𝑔̃ 𝑗,𝑙∘ 𝑋 𝑡−𝑙 𝑚𝑜𝑑 𝑁𝑁−1𝑙=0 , (2) 𝑡 = 0, ⋯ , 𝑁 − 1 and 𝑗 = 1, 2, ⋯ , 𝐿 , where ℎ̃ 𝑗,𝑙∘ and 𝑔̃ 𝑗,𝑙∘ are the 𝑗 th level MODWT wavelet and scaling filters (high- and low-pass filters) obtained by periodizing ℎ̃ 𝑗,𝑙 and 𝑔̃ 𝑗,𝑙 to length 𝑁 . These filters can be defined by renormalizing the DWT wavelet and scaling filters such as ℎ̃ 𝑗,𝑙 = ℎ 𝑗,𝑙/2 𝑗/2 and 𝑔̃ 𝑗,𝑙 = 𝑔 𝑗,𝑙/2 𝑗/2 . The multiresolution analysis (MRA) of the MODWT breaks up a signal into high- pass filtered “detail” components and a low - pass filtered “smooth” component. The MRA of the MODWT can be expressed as follows: 𝑋 = ∑ 𝐷 𝑗 + 𝑆 𝑗𝐿𝑗=1 , (3) 𝑋 = ∑ 𝐷 𝑗 + 𝑆 𝑗𝐿𝑗=1 , (4) 𝑆 𝑗,𝑡 = ∑ 𝑔̃ 𝑗,𝑙∘ 𝑁−1𝑙=0 𝑉̃ 𝑗,𝑡+𝑙 𝑚𝑜𝑑 𝑁 , (5) Where 𝐷 𝑗 is the wavelet “detail” at decomposition 𝑗 , and 𝑆 𝑗 is the wavelet “smooth” at decomposition 𝑗 . Figure 1 shows an example of the MODWT multiresolution analysis for a 10-second BCG signal. The 4 th level smooth coefficient (S4) can be used to represent the J-Peaks of the BCG signal. We briefly discuss the CWT in the next subsection. Continuous Wavelet Transform
The continuous wavelet transform is a time-frequency (more correctly, a time-scale) transform that is a useful tool for examining nonstationary signals. The CWT is a generalization of the short-time Fourier transform (STFT) commonly used to analyze nonstationary signals at multiple scales [33]. In a similar way to the STFT, the CWT applies an analysis window, that is, a wavelet, to extract segments from a signal. In contrast to the STFT, the wavelet is not only translated but dilated and contracted to take into account the scale of the activity under consideration. The wavelet’s dilation and contraction serve two purposes, that is, increasing the CWT’s sensi tivity o long- and short-time scale events, respectively. Given a continuous input signal 𝑥(𝑡) , the CWT can be defined as follows:
Figure 1 MODWT multiresolution analysis for a 10-second BCG signal. Wavelet Biorthogonal 3.9 (bior3.9) with 4 decomposition levels were opted to analyze the BCG signal. The maximum peaks of the 4 th level smooth coefficient (S4) correspond to the J-Peaks. The amplitude was normalized (z-score) for better visualization. 𝐶(𝑎, 𝜏) = ∫ 1√2 𝜓 (𝑡 − 𝜏𝑎 ) 𝑥(𝑡)𝑑𝑡, (6)
Figure 2 An example of a BCG signal (top), scalogram (center), and coefficients line at scale 20 (bottom). Gauss2 wavelet was designated for analyzing the BCG signal. The dashed black line on the scalogram was the scale (i.e., scale 20) where the J-Peaks were detected. The amplitude was normalized (z-score) for better visualization.
Where 𝜓(𝑡) is the mother wavelet, 𝑎 is a scale, 𝑏 is a shift parameter; 𝐶(𝑎, 𝜏) is a bivariate function obtained by mapping 𝑥(𝑡) and a wavelet scaled by 𝑎 at a given time 𝜏 . The localized correlation in time is determined over an integral starting with 𝑡 = 𝜏 and ending duration 𝑡 = 𝜏 + 𝐿 , where 𝐿 is the wavelet’s duration . It is noteworthy that short-term events (high-frequency signal components) such as spikes and transients can be determined when the wavelet is contracted (𝑎 < 1) , whereas long-time events (low-frequency signal components) such as baseline oscillations can be determined when the wavelet is stretched (𝑎 > 1) [33,34]. The result of the CWT can be shown in a graph known as a scalogram. It can be created by estimating the correlation between a signal and wavelets with different scales and then plotting how the correlation of each wavelet changes over a given period [33]. Figure 2 shows an example of the CWT for a 10-second BCG signal. Gauss2 wavelet was opted to analyze the signal, and the wavelet coefficient at scale 20 (scales 1 to 30) was used to detect the J-Peaks. .4. Template Matching
The template matching approach's main challenge is to choose the prototype/template and the similarity measure. The prototype (that is, cardiac cycle) was constructed from the first dataset due to the close contact between the MFOS mat and the participants. However, in the second dataset, the BCG signal ’s morphology was primarily affected by the massage chair's frequent movement. We specified the prototype as follows. Firstly, we manually selected high-quality BCG segments with a size of 30-second (1500-samples).
Figure 3
An ensemble averaging of BCG signals. The “I -J-
K” represents the ejection phase of the cardiac cycle.
The term “high - quality” implies that the segment does not include any signs of motion artifacts . Furthermore, cardiac cycles are easily identified. Secondly, we divided each 30-second segment into equal slices of 1-second (50-samples) with an overlap of 0.5-second. The redundancy created by the overlapped slices enabled us to accurately detect cardiac cycles, considering the relatively small sampling frequency of the MFOS. Thirdly, we only considered slices containing the “I -J- K” complex of the BCG signal.
The remaining slices were discarded from our analysis. The minimum peak distance (MPD) was 0.3 seconds, and it was selected using experimental observation. Several peak distances were evaluated, i.e., ranging from 0.2 to 0.7 seconds with a step size of 0.05 seconds (Figure 4). The MPD was appointed by examining the effects of two measures, i.e., the precision (Prec) and mean absolute rror (MAE), on HR detection across the 10 patients. Detected HR values were classified into correct and incorrect detections for each MPD. Then, the Prec was calculated to provides “ a rough estimate of how large portion of the detect HR values were correct ” , that is, “ how correct the detected HR values were ” [26,35,36]. It was calculated as follows: 𝑃𝑟𝑒𝑐 = 𝑐𝑜𝑟𝑟𝑒𝑐𝑡/(𝑐𝑜𝑟𝑟𝑒𝑐𝑡 + 𝑖𝑛𝑐𝑜𝑟𝑟𝑒𝑐𝑡) . The average MAE (i.e., between true and correctly predicted HR values) in beats per minute (BPM) tended to increase with increasing the distance. In addition, the average precision tended to decrease with increasing the distance. Therefore, the 0.3-second interval was assigned as an optimal interval to strike a balance between lower MAE (5.02) and higher precision (68.91%).
Figure 4 HR detection performance metrics (i.e., MAE and Prec) against the minimum peak distances.
Finally, the prototype was constructed by an ensemble averaging the valid slices centered at J-peak points (Figure 3). For each cardiac cycle, a candidate J-peak was detected by finding the maximum peak of the cross-correlation function (CCF) between the template and the BCG signal. The CCF is defined by calculating the correlation coefficients between the samples of the template ( 𝑥 ) and the BCG signal shifted by 𝑘 , ( 𝑦(𝑛 − 𝑘) ) [37]. The formula is as follows: 𝜌 𝑥𝑦 (𝑘) = 1𝑁 ∑ (𝑥(𝑛) − 𝑥̅) . (𝑦(𝑛 − 𝑘) − 𝑦̅) 𝑁−1𝑛=0 √( 1𝑁 ∑ (𝑥(𝑛) − 𝑥̅) ) . ( 1𝑁 ∑ (𝑦(𝑛) − 𝑦̅) ) , (7) Both signals were supposed to have 𝑁 samples each (that is, 50 samples). J-peaks with minimal intervals of 0.3 seconds were only deemed to be heartbeats. Heart Rate Detection
The HR was computed on a 30-second window and sliding the window by 15 seconds. The time window choice was based on previous studies [5,29,38], taking into account the sampling rate of the sensor (i.e., 50Hz). Regarding the MODWT-MRA, the Biorthogonal 3.9 (bior3.9) wavelet was appointed to determine cardiac cycles' shape. The bior3.9 wavelet proved to be the most suitable to characterize the profile of cardiac cycles across different avelets, precisely Daubechies 1 (db1), Symlet 2 (sym2), Coiflets 1 (coif1), and Reverse Biorthogonal 3.1 (rbior3.1) [12]. BCG signals were analyzed using 4 decomposition levels, and the 4 th level “ smooth/approximation ” component was employed for J-peaks detection [5,29,38]. The periodicity of the smooth coefficient reflected the same periodicity as the HR (Figure 1). At last, J-peaks were localized using a peak detector. For the CWT, Gaussian Derivative (GausP), Frequency B-Spline (FbspM-B-C), and Shannon (ShanB-C) wavelets were tested for HR detection. 𝑃 is an order-dependent normalization constant, 𝑀 is the spline order, 𝐵 is the bandwidth, and 𝐶 is the center frequency. For each wavelet, BCG signals were analyzed at different scales using the scalogram (Figure 2), and then the scale reflecting the same periodicity as the HR was designated for J-peaks detection [23]. The designated parameters and scales of the three wavelets are given in Table 2. For the rest of the paper, Gaus2, Fbsp2-1-1, and Shan1.5-1.0 will be used to refer to the CWT wavelets. Table 2 Parameters and scales of the CWT for HR detection; 𝑃 is an order-dependent normalization constant, 𝑀 is the spline order, 𝐵 is the bandwidth, and 𝐶 is the center frequency. Wavelet Parameters Scales GausP
𝑃 = 2
Range(1, 30), HR scale: 20
FbspM-B-C
𝑀 = 2, 𝐵 = 1, 𝐶 = 1
Range(1, 100), HR scale: 45
ShanB-C
𝐵 = 1.5, 𝐶 = 1
Range(1, 100), HR scale: 75 On the other hand, BCG signals obtained from the first dataset were used to construct a BCG template (training phase). The created template was then employed to detect HR in the second dataset, as outlined in subsection 3.4. The HR value at a time 𝑡 𝑛 , at which the nth maximum occurred, was defined as follows: 𝐻𝑅 𝑛 = 60𝑡 𝑛 − 𝑡 𝑛−1 , (8) Where 𝑡 𝑛 is the time at 𝑛 𝑡ℎ local maxima and 𝑡 𝑛−1 is the time at (𝑛 − 1) 𝑡ℎ local maxima in the designated MODWT component or CWT scale. Results and Discussion
This section presents the results of the three proposed methods (that is, MODWT-MRA, CWT, and TM) for HR detection across two independent BCG datasets (DataSet1 and DataSet2). For each method, the BPM error between the reference ECG and the MFOS was evaluated separately using the MAE and mean absolute percentage error (MAPE), and root mean square error (RMSE). All figures were generated using Python (Matplotlib, Plotline, Bokeh, and Seaborn).
Performance Evaluation of Heart Rate Detection — DataSet1
As presented in Table 3, the MODWT-MRA achieved overall MAE, MAPE, and RMSE, such as 4.71 (1.02), 7.61% (1.57), and 5.59 (0.96), respectively. The smallest and largest error metrics values were (3.13, 4.55%, 4.05) and (7.16, 9.99%, 7.83) for patients 1 and 8. Likewise, both patients had the highest and lowest precision, such as 96.53% and 30.77%. The past medical history of patient 8 indicated hypertension and dyslipidemia. Furthermore, this patient had severe obstructive sleep apnea (OSA) with an apnea-hypopnea index (AHI) of 78.2 [5]. The BCG signal’s quality was poorly affected by these health problems. Thus, reported Prec was low compared with other patients. Furthermore, this low Prec value triggered a large variability in the overall Prec, with a mean and standard deviation (SD) of 80.22% (18.03).
Table 3 HR detection performance metrics for DataSet1 using MODWT-MRA.
Metrics Patients Ids Mean (SD) 1 2 3 4 5 6 7 8 9 10 MAE
MAPE (%)
RMSE
Prec (%) can be more susceptible to patients’ comorbidities. As a result, cardiac cycles were not appropriately captured for various time intervals, triggering a total Prec of 69.57% (24.58). able 4 HR detection performance metrics for DataSet1 using CWT-Gaus2, CWT-Fbsp2-1-1, and CWT-Shan1.5-1.0. M e t h o d Metrics Patients Ids Mean (SD) 1 2 3 4 5 6 7 8 9 10 G a u s MAE 2.67 4.33 5.58 5.42 3.35 4.81 5.92 6.70 4.18 4.10 4.71 (1.16) MAPE (%) 3.88 6.51 9.52 9.30 5.25 8.61 10.61 9.17 6.86 6.11 7.58 (2.06) RMSE 3.51 5.25 6.44 6.34 4.22 5.84 6.74 7.40 5.03 5.04 5.58 (1.13) Prec (%) 98.35 89.75 73.11 72.63 96.45 58.45 68.43 45.98 93.57 91.62 78.83 (16.93) F b s p - - MAE 2.27 4.34 6.26 5.90 4.48 5.63 6.76 5.49 5.29 3.97 5.04 (1.24) MAPE (%) 3.32 6.64 10.13 9.45 7.15 9.48 11.48 7.23 8.69 6.00 7.96 (2.25) RMSE 3.07 5.22 7.04 6.75 5.32 6.48 7.51 6.38 6.17 4.85 5.88 (1.23) Prec (%) 99.18 91.11 47.22 45.73 93.33 39.67 31.75 76.22 76.28 95.21 69.57 (24.58) S h a n . - . MAE 3.36 4.46 5.76 5.61 3.60 4.81 5.50 6.57 4.97 4.18 4.88 (0.96) MAPE (%) 4.90 6.78 9.74 9.56 5.65 8.60 9.89 9.03 8.14 6.21 7.85 (1.73) RMSE 4.26 5.39 6.62 6.51 4.47 5.84 6.40 7.34 5.83 5.10 5.78 (0.93) Prec (%) 95.36 83.37 63.97 70.44 94.54 58.51 70.68 41.26 82.33 89.82 75.03 (16.47) Based on these results, Gaus2 seems to provide information about cardiac cycles more accurately than other wavelets. The total Prec accomplished by Gaus2, i.e., 78.83% (16.93), was slightly inferior to MODWT-MRA, i.e., 80.22% (18.03). Nonetheless, Gaus2 generated more favorable results with respect to the error measures (Table 4). The HR absolute error of each wavelet method is represented as histograms in Figure 5. It is clear from the figure that the HR detection performance of Gaus2 and MODWT-MRA is comparable. Moreover, Figure 6 shows the Bland-Altman plot of HR for the Gaus2 function across DataSet1. The limits of agreement (LoA) were computed as described in [39,40] given the fact that multiple observations per individual are available. The upper and lower LoA values are 10.95 and -11.17 BPM ( 𝑟 𝑚𝑟 = 0.38, 𝑃 < .001 ); 𝑟 𝑚𝑟 is the “ repeated measures correlation ” described in [41]. Heretofore, we examined the effectiveness of each wavelet function across DataSet1. After that, we discuss the performance of the wavelet functions and the template matching method across DataSet2. Figure 5 Distribution of the HR absolute error as histograms for each wavelet function across DataSet1.
Figure 6 Bland-Altman plot of the Gaus2 method across DataSet1.
Markers’ color was randomly assigned for each patient. .2.
Performance Evaluation of Heart Rate Detection — DataSet2
This particular dataset is challenging because BCG signals were acquired in a noisy environment. The signal quality was affected to a large degree by the massage chair’s movement and loss of contact with the MFOS. Regarding the wavelet-based functions, overall error measures (MAE, MAPE, and RMSE) for a) MODWT-MRA, b) Gaus2, c) Fbsp-2-1-1, and d) Shan1.5-1.0 were as follows: a) b) c) d) Figure 7 Overall performance measures (MAE, MAPE, RMSE) of the HR detection across DataSet2 using the 4 wavelet-based functions, i.e., MODWT-MRA, Gaus2, Fbsp2-1-1, and Shan1.5-1.0.
Overall, Gaus2 and Shan1.5-1.0 performed slightly better than MODW-MRA and Fbsp2-1-1. Akin to DataSet1, the error measures' largest values were scored by Fpsp2-1-1, such as 9.58, 16.68%, and 9.73. Figure 7 illustrates the HR detection error of each wavelet-based function across all participants. That said, Gaus2 scored the largest Prec value, i.e., 81.14% (14.17), whereas the Prec values for MODWT-MRA, Fbsp2-1-1, and Shan1.5-1.0 were as follows 77.12% (18.48), 76.24% (23.38), and 76.02% (14.44). The maximum, minimum, and overall Prec values for each wavelet-based function are specified in Table 5. Up to this point, the provided results demonstrate the superiority of Gaus2 for HR detection across two independent datasets. In the following, we discuss the findings of the template matching method and compare them with Gaus2. able 5 The maximum, minimum, and total values of the precision for the 4 wavelet-based functions (i.e., MODWT-MRA, Gaus2, Fbsp2-1-1, and Shan1.5-1.0) and the template matching approach.
Methods Precision (%) Minimum Maximum Mean (SD) MODWT-MRA
Gaus2
Fbsp2-1-1
Shan1.5-1.0
Template Matching
Figure 8 Overall performance measures (MAE, MAPE, RMSE) of the HR detection across DataSet2 using the template matching approach.
Overall, the TM achieved MAE, MAPE, and RMSE such as 4.74 (0.65), 7.46% (1.40), and 5.67 (0.68). As illustrated in Figure 8, the smallest values are 2.59, 4.04%, and 3.17 for participant 21, while the largest values are 5.96, 11.40%, and 6.84 for participant 13. Although the obtained results are reasonable, the total Prec, i.e., 72.83% (14.80), is not as good as Gaus2 (Table 5). Still, this fairly small Prec value is expected given the fact that the template was generated from a different dataset. Moreover, the BCG signals in DataSet2 were heavily corrupted by head movement artifacts. Figure 9 shows the Bland-Altman plot of HR for the TM approach across DataSet2. The upper and lower LoA values are 12.12 and -10.04 BPM ( 𝑟 𝑚𝑟 = 0.34, 𝑃 < .001 ). Similarly, it can be seen from he plot the relatively small number of HR points (i.e., Prec) in contrast to Gaus2. Additionally, 𝑟 𝑚𝑟 was 0.38; however, for Gaus2 it was 0.37. Figure 9 Bland-Altman plot of the TM approach across DataSet2 . Markers’ color was randomly assigned for each subject.
It is important to highlight that our findings for the Gaus2 and TM methods provided accepted overall aggregated results that were less than 10% MAPE. As stated in [42], an error rate of ±10% can be regarded as an accurate threshold for medical ECG monitors. To this end, preferring one method over the other will depend on the application requirements. To illustrate, the TM method produced a little better result in detecting HR than Gaus2 (Figure 10 and Figure 11). Yet, the total Prec was smaller than Gaus2. Thus, it will be more appropriate to use the TM in a well-controlled environment where motion artifacts can be minimized, despite that template selection is time-consuming and labor-intensive. On the other hand, the Gaus2 method seems less susceptible to motion artifacts. Hence, it can be more practical in real-life situations. Still, selecting an optimal wavelet function and scale requires prior knowledge about the BCG signal morphology. These two parameters will vary from one specific sensor to another. Besides, it should be pointed out that the HR detection results can differ significantly from one scale to another.
Figure 10 HR distribution for the reference ECG, Gau2, and TM methods for participant 2 (DataSet2). Time-windows were included in the diagram if they had been evaluated by both methods.
In other words, a particular scale and/or wavelet can only provide adequate results for individual cases while the opposite happens for other cases. A situation like this occurred, for example, with Fbsp2-1-1 in DataSet1 (Section 4.2). In addition to the results mentioned above, Gaus2 and TM methods were applied to a recently published dataset [43] and our findings are presented in Figure 14 and Figure 15 (Appendix). The processing time is another point that requires attention. The proposed methods were implemented using Python 3.9 on a computer with Windows 10 Pro installed (Intel(R) Core(TM) i7-4510U CPU @ 2.00GHz — MRA method was applied using the “ wmtsa-python ” library , while CWT-based methods were applied using the “ Scaleogram ” library . Basically, the average time needed to analyze a 30-second BCG signal was in a millisecond (ms) resolution for the 5 methods. Yet, the MODWT-MRA took less time compared to other methods, that is, 4.92 ms. Shan1.5-1.0 and Fbsp2-1-1 required more time to analyze a BCG signal, i.e., 46.3 and 47.5 ms. This performance is expected because a large number of scales were used, specifically 100. The time required for the TM was 14.4 ms. The improved performance for the MODWT- MRA occurred because the “ wmtsa-python ” library is written in Python and Cython. However, the “ Scaleogram ” library uses an adaptive convolution algorithm selection, that is, the scale processing loop switches to FFT-based convolution when the complexity is better in 𝑁 ∗ 𝑙𝑜𝑔2(𝑁) . The MODW-MRA is going to be more efficient for the applications that require real-time processing of the data considering its improved performance. https://gitlab.com/Cimatoribus/wmtsa-python https://github.com/alsauve/scaleogram Figure 11 Boxplots with p-values for Gaus2 and TM methods vs. the reference ECG across DataSet2. Table 6 The average time (in milliseconds) for analyzing a 30-second BCG signal using MODWT-MRA, Gaus2, Fbsp-2-1-1, Shan1.5-1.0, and TM methods.
Method Time (millisecond) MODWT-MRA
Gaus2
Fbsp2-1-1
Shan1.5-1.0 TM The sensor arrays were placed under the subjects’ a) head, b) chest, c) chest and abdomen, and d) under hip. The experiment was split into two phases: 10 minutes — supine sleeping and 10 minutes — side sleeping. The ideal sensor’s location was under the chest and abdomen. Similarly, optimal results were obtained by averaging the 6 sensors’ signals in the time domain [44]. As a result, the fused signal was employed to serve our purpose. We down-sampled the FBG signals to 50Hz so that cardiac cycles can match up to the BCG template. As given in Table 7, reasonable results were obtained. The total MAE, MAPE, and RMSE were , , and . able 7 HR detection performance metrics for DataSet3 using the TM approach. Metrics Subjects Ids Mean (SD) 1 2 3 4 5 6 7 8 9 10 MAE
MAPE (%)
RMSE
Prec (%)
Figure 12 HR detection results using the TM approach (subject 1, DataSet3). The top figure shows a 30-second BCG signal with the J-peaks annotated by up-pointing triangles. The bottom figure shows the corresponding ECG signal with the R-peaks annotated by up-pointing triangles.
Also, the total Prec was acceptable, that is, . This Prec value is fairly similar to a previous work in which the fusion was performed in the frequency domain using the cepstrum, and the total Prec reported was 84% [45]. Figure 12 demonstrates the performance of the TM approach for J-peaks detection across DataSet3. The top section of the figure shows a 30-second BCG signal with the J-peaks marked by up-pointing triangles. Besides, the bottom section of the figures shows the equivalent ECG signal with the R-peaks labeled by up-pointing riangles. Figure 13 shows the repeated measures correlation (rmcorr) plot for HR detection across DataSet3 using the TM method. Across the 10 subjects, Rmcorr and p-value were as follows: 𝑟 𝑚𝑟 = 0.39 𝑎𝑛𝑑 𝑃 < .001 . In summary, these results verify the potential of using a BCG template from a particular dataset to detect HR in a different dataset and/or under different conditions. The three methods (MODWT-MRA, CWT, and TM) described thus far have provided consistent results for HR detection from BCG signals. Furthermore, the average processing time seems reasonable enough to support real-time data analysis. Moreover, the total Prec values achieved by the three methods were fairly reasonable, considering that BCG signals were recorded in non-restricted environments; in other words, subjects ’ movements were allowed. Figure 13 Repeated measures correlation (rmcorr) coefficient plot [41] for HR detection across DataSet3 using the TM method. Conclusion
This paper examined the usefulness of three methods (i.e., MODWT-MRA, CWT, and TM) for HR detection from BCG signals. Four evaluation metrics (i.e., MAE, MAPE, RMSE, and Prec) were employed to assess the proposed methods’ performance for HR detection. The MODWT -MRA and CWT were applied to two datasets denoted as DataSet1 and DataSet2. The first dataset was obtained from 10 OSA patients during an overnight sleep study in which the MFOS was placed under the bed mattress. The second dataset was obtained from 50 subjects sitting in a massage recliner chair where the MFOS was placed on the chair's headrest. We excluded BCG signals from 11 subjects due to low-quality data. Hence, signals from 39 subjects were only considered in our analysis. The bior3.9 wavelet function (4 levels) was selected for the MODWT-MRA, and the 4 th level smooth coefficient was assigned for HR detection. That said, three wavelet functions were selected for the CWT, that is, Gaus2 (20 th scale of 30), Fbsp2-1-1 (45 th scale of 100), and Shan1.5-1.0 (75 th scale of 100). Across the two datasets, the CWT-Gaus2 method chieved more favorable outcomes than the other wavelets. For DataSet1, the total MAE, MAPE, RMSE, and Prec were 4.71 (1.16), 7.58% (2.06), and 5.58 (1.13), and 78.83% (16.93); however, for Dataset2, the total MAE, MAPE, RMSE, and Prec were 4.95 (1.17), 7.57% (1.71), and 5.77 (1.08), and 81.14% (14.17). For the TM method, a BCG template was created from DataSet1, that is, a training set, and then it was applied for HR across DataSet2, that is, a test set. The designated template contributed to satisfactory results despite the different experimental setup of DataSet2. The MAE, MAPE, RMSE, and Prec were 4.74 (0.65), 7.46% (1.40), and 5.67 (0.68), and 72.83% (14.80). The TM method achieved slightly better results than the CWT-Gaus2 but with a smaller Prec value. We double-checked the effectiveness of the TM method across a third dataset denoted as DataSet3. The dataset was collected from 10 subjects using a different BCG sensor, that is, FBG. Although BCG signals were generated from a different sensor, the TM method provided reasonably good results such as 3.43 (1.20), 5.51% (2.20), and 4.58 (1.19), and 80.88% (13.96). Finally, the average time required to analyze a 30-second BCG signal was in a millisecond’s resolution for the three methods. However, the MODWT-MRA had the highest performance, with an average time of 4.92 ms. In summation, the provided results showed the suitability of the wavelet-based methods and TM methods for HR detection across different datasets and/or under different experimental setups. Declaration of Competing Interest
None. eferences [1] T. Negishi, S. Abe, T. Matsui, H. Liu, M. Kurosawa, T. Kirimoto, G. Sun, Contactless vital signs measurement system using RGB-thermal image sensors and its clinical screening test on patients with seasonal influenza, Sensors (Switzerland). 20 (2020). https://doi.org/10.3390/s20082171. [2] M.J.M.M. Breteler, E.J. KleinJan, D.A.J.J. Dohmen, L.P.H.H. Leenen, R. van Hillegersberg, J.P. Ruurda, K. van Loon, T.J. Blokhuis, C.J. Kalkman, Vital signs monitoring with wearable sensors in high-risk surgical patients a clinical validation study, Anesthesiology. 132 (2020) 424 – D. Cimr, F. Studnička, Automatic detection of breathing disorder from ballistocardiography signals,
Knowledge-Based Syst. 188 (2020) 104973. https://doi.org/10.1016/j.knosys.2019.104973. [7] D. Cimr, F. Studnicka, H. Fujita, H. Tomaskova, R. Cimler, J. Kuhnova, J. Slegr, Computer aided detection of breathing disorder from ballistocardiography signal using convolutional neural network, Inf. Sci. (Ny). 541 (2020) 207 – – – – – – – – –
6. https://doi.org/10.1109/TBME.1985.325532. [19] J. Jin, X. Wang, S. Li, Y. Wu, A novel heart rate detection algorithm in ballistocardiogram based on wavelet transform, Second Int. Work. Knowl. Discov. Data Min. (2009) 76 –
79. https://doi.org/10.1109/WKDD.2009.98. 20] P. Smrcka, M. Jirina, Z. Trefny, K. Hana, New methods for precise detection of systolic complexes in the signal acquired from quantitative seismocardiograph, 2005 IEEE Int. Work. Intell. Signal Process. - Proc. (2005) 375 – – – Biol. Soc. EMBC’10. 2010 (2010) 2557– – th Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. EMBS’08 - "Personalized Healthc. through Technol. (2008) 1144 – – –
4. https://doi.org/10.22489/CinC.2019.145. [28] I. Sadek, J. Biswas, Z. Yongwei, Z. Haihong, J. Maniyeri, C. Zhihao, T.J. Teng, N.S. Huat, M. Mokhtari, Sensor data quality processing for vital signs with opportunistic ambient sensing, in: 2016 38th Annu. Int. Conf. IEEE Eng. Med. Biol. Soc., IEEE, 2016: pp. 2484 – – – – – — wavelets; multiscale activity in physiological signals, Biomed. Signal Image Process. Spring 2005. (2006) 1 – – – – – – Altman limits of agreement with multiple observations per individual, Stat. Methods Med. Res. 22 (2013) 630 – – – ppendix Table 8 Rmcorr, upper LoA, and lower LoA for wavelet-based methods and template matching across the three datasets.
Methods Rmcorr Upper LoA Lower LoA DataSet1 MODWT-MRA
Gaus2
Fbsp2-1-1
Shan1.5-1.0
Methods DataSet2 MODWT-MRA
Gaus2
Fbsp2-1-1
Shan1.5-1.0 TM Methods DataSet3 TM
Figure 14
Overall performance measures (MAE, MAPE, RMSE) of the HR detection across DataSet4 [43] using Gaus2 and TM methods. DataSet4 was obtained from 40 subjects using 4 EMFi sensors placed underneath the bed mattress. The TM method failed to analyze BCG signals from subjects X1001 and X1005 because the signal’s morphology was quite different compared to the MFOS. The HR was computed by fusing the 4 EMFi signals using a pairwise maximum operation. Average fusion was also examined. Nevertheless, the fused signal was distorted and the main features of a typical BCG signal (i.e., “I -J- K” omplexes) were missed. Overall, the MAE, MAPE, RMSE were (2.15 (2.30), 2.91% (2.79), and 2.62 (2.39)) and (4.07 (2.31), 5.79% (2.87), and 4.76 (2.26)) for Gaus2 and TM methods, respectively. The Gaus2 method achieved close results to the reference ECG compared to the TM. Furthermore, the total Prec attained by Gaus2, i.e., 93.65% (16.36) was better than the TM, i.e., 84.08 (20.61). Rmcorr and p-value were as follows: 𝑟 𝑚𝑟 = 0.65 𝑎𝑛𝑑 𝑃 < .001 and 𝑟 𝑚𝑟 = 0.35 𝑎𝑛𝑑 𝑃 < .001 for Gaus2 and TM. Figure 15 Superimposed histograms of the ECG — HR and Gaus2 ——