**1. Introduction**

Advancements in modern technologies enabled field monitoring of some parameters of human health [1]; for example, heart monitoring is used for o ff-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 di fferent 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 cuto ff 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–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-su fficient mobile sensor. The main contributions of this paper are as follows:


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.

**Figure 1.** Photodetector(PD) waveforms: raw signals obtained from photoplethysmography (PPG) sensor in reflective configuration placed on the index tip (dotted line) and on the wrist (continuous line).

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.

(**a**) 

**Figure 2.** *Cont*.

**Figure 2.** Waveforms of the signals from the PPG detector: (**a**) the raw signal; (**b**) the raw signal from (**a**) filtered only by a band-pass filter; and (**c**) the raw signal from (**a**) filtered first with the limiting section described in the paper and then by the same biquadratic filter as used in (**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):

$$y\_i = \frac{1}{a\_0}(b\_0\mathbf{x}\_i + b\_1\mathbf{x}\_{i-1} + b\_2\mathbf{x}\_{i-2} - a\_1y\_{i-1} - a\_2y\_{i-2}),\tag{1}$$

where *yi* is *i*-th output sample; *xi* is *i*-th input sample; *a*0, *a*1, and *a*2 and *b*0, *b*1, and *b*2 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:

$$\begin{array}{c} \text{if } y\_i = \begin{cases} \quad y\_{C,i} \text{ if } L\_L \le y\_{C,i} \le L\_H\\ \quad L\_L \text{ if } y\_{C,i} < L\_L\\ \quad L\_H \text{ if } y\_{C,i} > L\_H \end{cases} \end{array} \tag{2}$$

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 <sup>o</sup>ff. 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.

**Figure 3.** Block diagram of PPG signal acquisition, preprocessing, and final HR estimation by TDHR algorithm.

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* = [*<sup>x</sup>*1, *x*2, ... , *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):

$$m\_{k,i} = \begin{cases} 0 & \text{x}\_{i-1} > \text{x}\_{i-k-1} \land \text{x}\_{i-1} > \text{x}\_{i+k-1} \\ r+1 & \text{otherwise} \end{cases},\tag{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.

JJJHaving 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*]:

$$\gamma\_k = \sum\_{i=1}^{N} m\_{k,i} \text{ for } k \in \{1, 2, \dots, L\} \tag{4}$$

and the global minimum of γ, λ = arg min(γ*k*). The value of γ is used to obtain the matrix Mr, 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:

$$\sigma\_{i} = \frac{1}{\lambda - 1} \sum\_{k=1}^{\lambda} \left[ \left( m\_{k,i} - \frac{1}{\lambda} \sum\_{k=1}^{\lambda} m\_{k,i} \right)^{2} \right]^{\frac{1}{2}}. \tag{5}$$

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 e fficient 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:

$$
\lambda > \lambda\_{\text{max}} \tag{6}
$$

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 Mr with the matrix Mr' containing 1-bit binary values, *<sup>m</sup>*'*k,i*, as follows:

$$m'\_{k,i} = \begin{cases} 0 & \mathbf{x}\_{i-1} > \mathbf{x}\_{i-k-1} \land \mathbf{x}\_{i-1} > \mathbf{x}\_{i+k-1} \\ 1 & \text{otherwise} \end{cases} . \tag{7}$$

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:

$$s\_i = \sum\_{k=1}^{A} m'\_{k,i}.\tag{8}$$

As the *<sup>m</sup>*'*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*.

$$t'\_i = \frac{t\_{i-1} + t\_i}{2} \tag{9}$$

for which the following simple condition (10) holds for *si* from Equation (8):

$$s\_{i\cdot 2} > 1 \land s\_{i\cdot 1} = 1 \land s\_i = 1 \land s\_{i+1} > 1. \tag{10}$$

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.

**Figure 4.** Sample of (**a**) filtered input PPG signal, (**b**) calculated and normalized values of σ*i* from Equation (5) and *si* from Equation (8), and (**c**) calculated values of γ*k* from Equation (4).

**Figure 5.** An example of special case of a "flat" peak consisting of two equal values at *ti*-1 and *ti*.

**Figure 6.** Comparison of the results of the peak detection algorithm proposed in this paper (**a**) with the AMPD peak detection algorithm (**b**) from [19] on the same waveform. The arrows indicate the peaks not detected by one of the algorithms but detected by the other one. All other peaks were detected at the same positions, by both methods.

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.

**Figure 7.** Elimination of the patient's movement from pulse-signal detection: (**a**) the PPG signal distorted with the patient's movement; the areas in gray are excluded from pulse period calculation due to the detected movement; (**b**) acceleration values read from the accelerometer placed together with the heart rate detector on the wrist, (**c**) and the values of movement *di* from Equation (11); (**d**) signal *di* after threshold, as in Equation (12), and with the eliminated short-term peaks (the eliminated peak is shown with a dotted line).

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:

$$d\_i = \frac{1}{G} \sqrt{(g\_{x,i} - g\_{x,i-1})^2 + \left(g\_{y,i} - g\_{y,i-1}\right)^2 + \left(g\_{z,i} - g\_{z,i-1}\right)^2},\tag{11}$$

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*:

$$V\_i = \begin{cases} 1 & d\_i > H \\ 0 & \text{otherwise} \end{cases} \tag{12}$$

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):

$$IP = flow\left(\frac{N \cdot BPM\_{\text{min}}}{60f\_s}\right) \tag{13}$$

where *fs* is the sampling frequency, and *BPM*min 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:

$$MAE = \frac{1}{W} \sum\_{j=1}^{W} \left( BPM\_{est}(j) - BPM\_{ref}(j) \right). \tag{14}$$

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 1.** Comparison of the accuracy of the presented algorithm to the algorithms from the literature on the large PPG-DaLia dataset as MAE (bpm). The proposed TDHR algorithm was tested in several versions, with various numbers of analyzed samples *N* (for all other algorithms *N* = 256), sampled with a frequency of 32 Hz, with a *BPMmin* = 40 bpm. The measurements were done every 2 s. For the TDHR algorithm, the accuracy was also calculated for valid measurements, and then the percentage of valid measurements was also given.


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*.

**Table 2.** Comparison of the performance versus computational cost of the TDHR algorithm and the algorithms from the literature. The computational cost of the TDHR algorithm was calculated as the number of arithmetical operations and the number of memory bytes needed for algorithm realization in a microcontroller.


As can be seen from Tables 1 and 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 Tables 1 and 2 can help to select *N* to achieve a compromise between the accuracy and computational cost.
