1. Introduction
Advancements in modern technologies enabled field monitoring of some parameters of human health [
1]; for example, heart monitoring is used for off-hospital monitoring and in fitness and professional sport activities. The most commonly measured value is the heart rate (HR), although advanced applications also use other values, e.g., pulse irregularity, as well as biometric identification or analysis of accurate electrical signals that cause heart contraction, i.e., electrocardiography (ECG) [
2]. Accurate ECG requires connecting electrodes to the patient’s body in several different places, which is inconvenient for the patient, and it can be used only in certain situations. A much more convenient method is measuring the pulse on the wrist by using photoelectric methods. The skin of the wrist is irradiated with single or multicolor light, and then the reflected light is measured. The intensity of the reflected light depends on the absorption of the skin, which depends on the blood volume supplied to the tissues. In this way, the received signal contains information about the current blood supply to the vessels near the measuring device. This method, introduced by Hertzman [
3], is known as photoplethysmography (PPG). Unfortunately, PPG signals obtained from a moving person’s wrist are weak, distorted, and contain noise. The noise level is often higher than a usable PPG signal. Correct analysis of a low-quality PPG signal is a very challenging task and can consume significant processing time, energy, and resources.
Heart rate estimation from wrist PPG is now a popular area of investigation, and many algorithms for such a task have been proposed [
4]. Some of them require significant computing power and memory usage, blocking their application in small portable low-power devices. There are two main approaches used: time-domain sophisticated filtering and frequency-domain processing. Both are often accompanied by a movement sensor (accelerometer) for movement-based artefacts/spectrum removal.
A straightforward approach toward heart rate estimation is peak detection in a periodic signal. One of the simplest possibilities is to use threshold or auto-threshold values in the signal time window [
5,
6]. Another way is to use transforms such as continuous wavelet [
7] or Hilbert [
8]. In [
9], a nonlinear filter bank was used, with variable cutoff frequencies. In [
10], the detection algorithm was based on a time-varying autoregressive model, with a Kalman filter used for autoregressive parameter estimations. Hidden Markov models were used in [
11] for combining structural and statistical knowledge of the signal in a single parametric model. A neural network adaptive whitening filter to model the lower frequencies of the signal is presented in [
12]. Much recent work is concentrated on methods using time-frequency spectra [
13,
14] and deep-learning approaches [
15,
16,
17]. A review of the current state-of-the-art signal-processing techniques for HR estimation from a wrist PPG signal can be found, for example, in [
4,
18].
The proposed solution of time-domain heart rate measurement algorithm (TDHR) consists of three main blocks: signal conditioning, peak detection, and heart-rate-measuring blocks. In this work, the modified Automatic Multiscale-based Peak Detection (AMPD) algorithm from [
19], together with a bandpass filter/limiter, was used for finding the HR from a wrist-based PPG signal. All of the signal processing was done while taking into account the need for low power, low resources, and computing power utilization, necessary for a self-sufficient mobile sensor. The main contributions of this paper are as follows:
A two-stage input-signal-conditioning digital nonlinear filtering block with limiter;
Application of the AMPD algorithm for HR peak detection;
Modification of the AMPD algorithm toward efficient implementation in low-power resource-constrained hardware;
The proposition of a time-domain heart-rate-measuring algorithm with an accelerometer-based false measurement removal.
This paper is organized as follows. In
Section 2, the signal conditioning block is described, which is followed by the peak detection algorithm described in
Section 3 and the heart rate calculation block presented in
Section 4. The proposed algorithm was evaluated by using a multi-hour dataset; the results of the evaluation and the comparison to the other solutions from the literature are presented in
Section 5.
Section 6 contains the details of the implementation in a low-power wearable device. Finally, conclusions are provided in
Section 7.
2. Signal Conditioning
PPG-signal measuring is usually done with an infrared (IR), red, or green light-emitting diode (LED) as the light source and a photodetector (PD) receiving the reflected or transmitted light. Often, instead of one, a few LEDs and/or PDs are used, providing the possibility to choose the best observed signal or to use additional preprocessing, such as, for example, averaging. Two measure modes can be used: reflectance and transmission. In the first case, the LED and PD are placed on the same skin surface, close to each other; the typical spacing between the PD and LED is in the range of 5–15 mm. The PD measures the light reflected from the tissue. In transmission mode, the LED and PD are placed on the opposite sides of a human body part, and the light reaching the PD goes through the whole body part. The commonly used places on a human body for such PPG devices are the fingertip, wrist, earlobe, forehead, torso, ankle, and nose [
18,
20]. The wrist is the most convenient place to mount the monitoring device for everyday life, but this is not optimal regarding the strength of the signal. Due to the relatively large thickness and the presence of bones, only the reflectance method of measurement is practical on the wrist. In
Figure 1, the waveforms of the signal received from the sensor in reflectance mode on the index finger and on the outer wrist, where a watch is usually worn, are presented. A useful signal carrying information about the pulse is present in the form of the peaks with the period of about 30–40 samples superimposed on the curve. As shown in
Figure 1, the amplitude of these peaks for the signal measured at the finger is about 200, while the amplitude measured at the wrist is much lower, about 30–50, while the average signal level (baseline) is about 11,100.
As can be seen from the waveforms from
Figure 1, the signal from the wrist is much weaker; moreover, each signal has a variable offset, and it contains noise and interference. To extract the heart rate in a time-domain, the signal must be preprocessed, so that the peak detection algorithm can find the peaks. The simplest approach is to use a band-pass filter that can help to reduce the noise and eliminate the constant component of the signal. Such a solution is in common use [
4]. However, the signal at the input of the filter can have significant changes of the constant component, as well as signal fluctuations caused by wrist movement, resulting in large peaks and oscillations at the output of the band-pass filter, as can be seen in
Figure 2a,b.
The proposed approach for PPG signal acquisition and processing is presented in
Figure 3. The raw signal obtained from the PD is fed, in the first step, to the band-pass biquad section, with an internal limiter. The output value of the standard biquad section using direct version I is given by Equation (1):
where
yi is
i-th output sample;
xi is
i-th input sample;
a0,
a1, and
a2 and
b0,
b1, and
b2 are, respectively, the denominator and nominator coefficients of the biquad transfer function. The biquad section with the internal limiter works in two steps: First, it calculates candidate
yC,i for the output value, according to Equation (1), and then it uses Equation (2) to update the actual output value:
where
LL and
LH are the limiter’s parameters. These parameters should be selected so that the PPG signal from the heart contractions will remain intact, while distortions resulting from hand movements should be cut off.
Figure 2a shows the PPG signal containing distortions taken from a wrist. For this case, the limiter parameters should therefore cut off these high distortions, while the small sawtooth waveform should not be affected. The signal after the bandpass filter does not contain a constant component, as is presented in
Figure 2b. In the case of the signal from
Figure 2b, good choices for the
LL and
LH parameters’ values could be −150 and 150, respectively.
Figure 2c shows the results of preprocessing for such a limiter, together with the band-pass biquad section. In general, these values for a given PPG system should be selected as about 100–150% of the negative and positive amplitude of the useful signal, respectively.
In the second stage of the proposed preprocessing block, a typical fourth-order band-pass filter built from two biquads was employed. The band-pass of both stages was set to 0.5–2.5 Hz.
The introduction of the limiting section in the form of a single biquad stage significantly reduces the rapid changes of the signal at the input of the band-pass filter and improves its recovery after a rapid change of the constant component in the input signal. An example of distorted signal, together with the accelerometer readouts, is presented later in this paper.
3. Peak Detection
The conditioned signal from the detector was used as an input to the block responsible for finding the peak values of the signal; the heart rate can then be easily calculated from the detected peaks. The detection of peaks is based on the AMPD algorithm [
19]. The AMPD algorithm has the capability to work with noisy periodic and quasi-periodic signals. It needs the input signal to be linearly detrended, but the use of the input filter of band-pass characteristic with the limiter described in the previous section satisfies this requirement. The AMPD algorithm performs well for the filtered PPG signal, but it is computationally expensive, which can be unacceptable for wearable devices. The need to calculate a large matrix with real-valued elements, where moving windows are used, can be avoided due to the modifications of the algorithm proposed in further parts of this paper. This section starts with the detailed description of the AMPD algorithm, and then the proposed modifications are introduced. In this way, the authors wanted to make it easier for the reader to track changes that were applied to the original algorithm, without having to refer to the reference.
The main part of the AMPD algorithm is the Local Maxima Scalogram (LMS) matrix M of elements,
mk,i, which is calculated for the discrete uniformly sampled signal
x = [
x1,
x2, …,
xN] in the analyzed window, where
N is the constant number of samples, and scale
k defines the moving window of varying length,
wk, according to Equation (3):
where
wk = 2
k |
k = 1,2, …,
L,
k is
k-th scale of the signal,
L = ceil(
N/2) − 1, and
r is a uniformly distributed random number of values [0,1]. The values of
mk,i are calculated for every scale
k and
i =
k + 2,
k + 3, …,
N −
k + 1.
Having calculated the LMS matrix, the next step in the AMPD algorithm is to calculate the scale-dependent distribution of zeros in the LMS, by calculating vector
γ = [
γ1,
γ2, …,
γk]:
and the global minimum of
γ,
λ = arg min(
γk). The value of
γ is used to obtain the matrix M
r, which is the matrix M with deleted bottom rows for
k >
λ. The peaks are then found for indexes,
i, for which the column-wise standard deviation,
σi, is equal to zero:
Parameter λ enables the determination of how many rows of matrix M should be used for calculating the matrix Mr. The noisier the signal from the PPG sensor, the larger the value of λ.
In this paper, the authors propose modifying the AMPD algorithm toward a more efficient implementation. From practical observation, it has been inferred that the signal is too noisy, and it is of no use for peak detection and heart rate calculation, when we have the following:
The value of λmax = 17 was found empirically to be a good choice. Therefore, it is practical to assume a priori, that the signals for which the condition (6) is true are of no use and the full matrix M is not used. To save the memory, only the matrix Mr, instead of M, can be calculated and stored together with the vector γ.
Processing and storing floating point values,
mk,i, requires storage space and processing resources. To simplify the processing, the authors propose replacing real-valued elements of matrix M
r with the matrix M
r’ containing 1-bit binary values,
m’
k,i, as follows:
This saves a lot of the device’s memory, requires only integer operations, and results in another simplification: Instead of calculating the column-wise standard deviation
σi from Equation (5), calculating the column-wise summation presented in Equation (8) can be used:
As the
m’
k,i are 1-bit binary values,
si can be calculated with fast integer summation. The indices of peaks
pi can be located by finding all indices,
i, for which
si = 0. The values of
σi and
si, together with the values of
γk were shown in
Figure 4b,c, for an exemplar input signal (
Figure 4a).
The PPG signal from a wrist is weak, as can be seen in
Figure 2c, and when sampled, the situation of a “flat” peak with two equal values for samples
ti-1 and
ti, as shown in an example in
Figure 5, can sometimes occur. Such a peak would neither be detected by the AMPD algorithm nor its modified version proposed in this paper. To improve the peak detection in such a case, an additional rule was introduced: The peak is also detected at time
t’
i.
for which the following simple condition (10) holds for
si from Equation (8):
The proposed simplified peak detection algorithm was compared to the original AMPD algorithm presented in [
19]. For this purpose, the PPG signals from the PPG-DaLia database [
15] were used, and the detected peaks were compared. As the original AMPD peak detection algorithm uses a random variable to fill the LMS matrix, its results are slightly different from run to run, depending on the seed value of the random number generator. The results are very similar, while the modified version can be implemented more efficiently. An example is presented in
Figure 6: the arrows indicate the missing peaks, which were exclusively detected by the other algorithm. It can be seen that only a few peaks are differently detected.
The presented example uses an unfiltered PPG signal to present the robustness of the algorithm. The signal after filtration would be more regular; thus, the differences in the results between the two versions of the algorithm would be even smaller.
4. Heart Rate Calculation
In an ideal situation, where the patient is not moving, the detected peaks from the peak detection algorithm can be directly used to calculate the heart rate. However, the PPG signal can be seriously distorted by the movement of the patient. Each movement of the patient can cause a change in the saturation of the tissues with blood, and a sensor displacement on the wrist, which results in artefacts in the signal received by the optical detector. To mitigate this problem, the accelerometer is used for detecting the movement of the patient’s hand. The movement data are aligned with the detected peaks, and the periods between the peaks which are affected by movements are discarded from the calculation of the pulse period. The example of the elimination of the patient’s movement is presented in
Figure 7.
To find the time periods that will be excluded from period calculations, first the acceleration values from the accelerometer are differenced to obtain the movement indicator,
di, according to the following equation:
where
i is the index of the sample;
gx,i,
gy,i, and
gz,i are the acceleration values read from the accelerometer for axes X, Y, and Z, respectively; and
G is a constant value used for normalization to obtain
di ∈ [0,1] for all
i and for the acceleration values in
gx,i,
gy,i, and
gz,i in the accelerometer’s full measuring range. The movement values,
di, are then compared to the constant,
H, to obtain a digital binary signal
Vi:
Signal Vi is then filtered in time, to eliminate spikes longer than Ts. The values of H and Ts were experimentally set to H = 0.0025 and Ts = 500 ms for the accelerometer range ±2 g.
For the final result of the heart rate, the inverse of the median of the periods between the peaks not affected by the movement is calculated. For a valid result, a minimal number of peaks,
P, is required, and it is calculated as shown in Equation (13):
where
fs is the sampling frequency, and
BPMmin is the minimal heart rate required to be measured by the system.
5. Evaluation on Dataset
For the purpose of evaluating the algorithm, the PPG-DaLia database [
15], containing more than 35 h of data recorded from 15 persons, was used. The database contains signals collected from a PPG sensor, accelerometer, and ECG, where the ECG is used as ground truth. As described in [
15], ground truth heart rate values were obtained from an ECG signal processed by R-peak detection algorithm [
21]. Then, the detection results were manually inspected and corrected, mainly for a few cases where significant motion was observed. The ECG signal was segmented with a shifted window; ground truth heart rate was finally calculated as the mean heart rate within each window. The signals were collected during eight different types of typical daily life activities, under controlled but close-to-real-life conditions. This dataset is the longest one available to the authors. The accuracy of the algorithm was evaluated by using the mean absolute error (
MAE) metric of beats per minute (bpm), calculated by using the sliding window approach of length 8 s, with a 2 s shift, according to the following equation:
where
W is the total number of windows,
BPMest(
w) is the heart rate in bpm for window
j, and
BPMref(
j) is the reference heart rate obtained from the ECG ground truth signal for the same window
j. This evaluation method is commonly used in related work [
14,
15,
22,
23].
The quality of the signal and the possibility of the measurement are checked in the TDHR algorithm. When the measurement is not possible due to not satisfying Equation (6), i.e.,
λ >
λmax or due to movement detected by the accelerometer affecting all of the periods between the detected peaks, the measurement result is indicated as invalid. The evaluation of the algorithm is performed in two ways: (i) When the result is not available, the last valid result is used; (ii) only valid results are used to calculate the performance metric; and then the percentage of the valid samples is also given. The results of the evaluation, together with the performance of the SpaMa [
14], SpaMaPlus [
15], Schaeck2017 [
23], CNN average, and CNN ensemble [
15] algorithms are presented in
Table 1. The TDHR algorithm was evaluated for several lengths of the sliding window:
N = 1024 (32 s),
N = 512 (16 s),
N = 256 (8 s), and
N = 128 (4 s).
Table 2 presents the comparison of the computational cost of several algorithms. For the TDHR algorithm, the number of operations per second was estimated from a manual analysis of the C code as the number of arithmetic operations needed to obtain a single heart-rate result. The TDHR algorithm requires only a few parameters, as opposed to the CNN algorithms, but it needs storage memory for calculating the LMS array; the size of this memory depends on the window length,
N.
As can be seen from
Table 1 and
Table 2, the performance of the proposed TDHR algorithm is similar to
SpaMa,
SpaMaPlus, and
Schaeck2017, while it is worse than
CNN average,
CNN ensemble, and
CNN constrained. The computational costs of
SpaMa,
SpaMaPlus, and
Schaeck2017 are not known, but they can be high, as those algorithms require the calculation of the power spectral density and the analysis of the PPG spectrum. The CNN-based algorithms, even
CNN constrained, require larger computational costs than most versions of the TDHR algorithm. The proposed algorithm uses the mechanism of removing the time periods with body motion registered by the accelerometer, so there can be gaps between valid measurements in the case of long-lasting motions. The presented results in
Table 1 and
Table 2 can help to select
N to achieve a compromise between the accuracy and computational cost.
6. Implementation
The proposed algorithm was implemented in low-power wearable hardware and tested. The hardware consists of the main board, power supply PCB, with external coil for wireless inductive charging, and a LiPo battery of capacity 110 mAh. A block diagram of the hardware is presented in
Figure 8. As the processing unit, a PSoC6 microcontroller from Cypress was used. This is an ultra-low-power microcontroller with dual processor architecture: Arm Cortex M4 and M0+ cores. A BH1790GLC optical sensor from Rohm detects the PPG signal and also drives the four green 527 nm LEDs, equally placed in a circle of diameter of 6 mm, with the optical sensor placed in the center. The detector detects the light emitted by the LEDS and reflected from the patient’s skin. An LSM6DSL accelerometer from STMicroelectronics is mounted on the same board as the optical sensor. The optical sensor and the accelerometer are connected to the microcontroller via an I
2C bus, used for configuration and data transfer.
The software for heart rate measurements runs on the Cortex M4 processing core of the PSoC6 microcontroller. The algorithm presented in this paper was written in C as two tasks running on the FreeRTOS operating system: The first task constantly reads data from the sensors, realizes filtering, and buffers data. The second task starts periodically and executes the peak detection algorithm, using data from the buffer. The calculated heart rate is sent wirelessly, using the Bluetooth Low Energy protocol. For this purpose, a built-in BLE transceiver in PSoC6 device was used. According to the BLE nomenclature, in the proposed solution, the BLE transceiver block was configured to perform a peripheral (a device constrained in resources such as energy and computing power) and server role (a device working as a data source and sending that data to the remote master device). As a data format, standard BLE Heart Rate Profile was used. The values of the measured heart rate are transmitted periodically; the user can receive the transmitted values by connecting any Bluetooth receiver compatible with Heart Rate Profile.
The PPG signal is sampled with 14-bit resolution and with frequency fs = 32 Hz, using a buffer of N = 128… 1024 samples, which enables data from 4 to 32 s to be analyzed, depending on the selected value of N, compromising accuracy versus computational cost and memory usage.
The prototype device was built and installed inside a custom-made 3D-printed case with a rubber strap, as shown in
Figure 9. The size of the printed case (without the rubber strap) is 11.7 × 26 × 46 mm. The cost of the components, including the battery, was about
$60, at retail prices, for 15 pieces. The proposed implementation was capable of continuous measuring of the heart rate for more than 24 h.
7. Discussion and Conclusions
Nowadays, there are many smart watches on the market that are capable of measuring the heart rate. The details regarding the optical part of the measurement, such as the number of sensors or the wavelengths of the LEDs used during measurement, are often revealed. However, the details regarding the algorithms used are not available. Most of the publications focus on measuring the accuracy of the popular devices. In [
24], the authors measured the performance of Apple’s iWatch Apple Watch Sport 42 mm (first generation), during cardiopulmonary exercise test (CPET). They observed MAE from 6.34 to 7.55. It is difficult to exactly compare this to our results, as we use different and longer test conditions. In fact, the details of the PPG algorithms can only be found in the scientific publications, where the authors try to increase the accuracy of the measurement by using novel ideas and powerful techniques.
In this paper, a time-domain algorithm for the real-time detection of the heart rate was presented. The algorithm is aimed at wearable, resource-constrained devices, where battery capacity, processor speed and memory size are constrained. The algorithm processes raw data in a time-domain and requires only a few parameters. The approach is simple, the proposed algorithm has a reasonable accuracy, and it can be implemented in a typical (not DSP) microcontroller.
The algorithm consists of a two-stage input-signal-conditioning block with a limiter, the peak detection block, and period calculation block. The two-stage input conditioning block is built out of two digital bandpass filters, where the first filter has been modified to process data nonlinearly, to provide fast recovery after large signal transients. The use of band-pass filters is very simple to implement; it appeared sufficient and very effective for conditioning the signal; therefore, other methods such as wavelet-based baseline removal would not need to be considered.
The peak detection block enables operation at significantly lower processing and implementation costs, compared to the original AMPD algorithm. The period calculation block uses the median to calculate the heart rate based on the time differences between the peaks, with the use of an accelerometer to exclude the time periods affected by body movement.
In this study, most of the parameters were set experimentally, as we targeted on a simple and economical hardware implementation. It would be interesting to provide a method of automatic and dynamic adjustment of the parameters to further reduce the computational cost, basing on the input signal quality and the movement readouts. This will be a topic of our further research.
The proposed algorithm was evaluated in several variants, for different sliding window lengths, N, providing the possibility to compromise the accuracy versus lower operational costs. The authors mainly used N = 1024, which seems to be a good compromise between the calculation cost and the accuracy. The proposed algorithm was compared to the other algorithms from the literature. The achieved accuracy is comparable to the other algorithms at smaller computational costs. The proposed solution was also implemented in a low-power wrist-wearable device.