Increasing accuracy of pulse arrival time estimation in low frequency recordings

Objective. Wearable devices that measure vital signals using photoplethysmography are becoming more commonplace. To reduce battery consumption, computational complexity, memory footprint or transmission bandwidth, companies of commercial wearable technologies are often looking to minimize the sampling frequency of the measured vital signals. One such vital signal of interest is the pulse arrival time (PAT), which is an indicator of blood pressure. To leverage this non-invasive and non-intrusive measurement data for use in clinical decision making, the accuracy of obtained PAT-parameters needs to increase in lower sampling frequency recordings. The aim of this paper is to develop a new strategy to estimate PAT at sampling frequencies up to 25 Hertz. Approach. The method applies template matching to leverage the random nature of sampling time and expected change in the PAT. Main results. The algorithm was tested on a publicly available dataset from 22 healthy volunteers, under sitting, walking and running conditions. The method significantly reduces both the mean and the standard deviation of the error when going to lower sampling frequencies by an average of 16.6% and 20.2%, respectively. Looking only at the sitting position, this reduction is even larger, increasing to an average of 22.2% and 48.8%, respectively. Significance. This new method shows promise in allowing more accurate estimation of PAT even in lower frequency recordings.


Introduction
Due to the prevalence of wearable devices such as smartwatches, an abundance of continuous data on heart activity is available in the form of photoplethysmography (PPG) signals.PPG is an optical measurement of changes in blood volume in the tissue underneath the sensor location and reflects heart activity non-intrusively and non-invasively (Allen 2007).From PPG, several vital signals can be determined, such as heart rate and heart rate variability (Biswas et al 2019, Galli et al 2022).In combination with an electrocardiogram (ECG), PPG can be used to calculate the pulse arrival time (PAT), which is widely used to measure blood pressure non-invasively and non-intrusively over longer periods of time (Zhou et al 2023).
PAT is defined as the time between the R-peak of a pulse in the ECG signal and the start of the same pulse in the PPG signal (Peng et al 2021), and represents the time between the depolarization of the ventricles and the arrival detection at the site of the PPG measuring device.PAT can be used for tracking blood pressure (Deshmukh et al 2022) and has also been shown to be able to detect mechanical alternans (van Duijvenboden et al 2019); and while the pulse transit time more directly translates to blood pressure, PAT is often used as a substitute (Finnegan et al 2021).
PPG is often embedded in wearable devices, such as smartwatches, to enable continuous monitoring of cardiac activity and blood pressure through PAT.One of the facets of wearable design is battery consumption, which must be minimized in order to improve the usability of the device.To reduce the frequency of having to charge the device, its battery size may be increased at the cost of bigger devices, which might negatively affect their non-obtrusiveness.An alternative solution is to optimize the load of signal acquisition and processing methods that are executed by the software and hardware.One possible way to do this is to reduce the sampling frequency at which the PPG and ECG signals are acquired (Kim and Baek 2023).
Medical applications require accurate measurement of the vital signs, such as PAT, and a low sampling frequency can lead to a decrease in accuracy.In most medical applications a higher sampling frequency is therefore preferred (Peláez-Coca et al 2022, Wang et al 2022).Hence, there is an increased interest in technologies that can extract health parameters with high accuracy in low sampling frequency recordings.
Recently, diastole-patching has been proposed to identify the start of the pulse in the PPG signal (Vardoulis et al 2013).Diastole-patching obtains a so-called 'patch' of a pulse, by taking the part of the pulse centered around the diastolic minimum, up to the maximum of the first derivative of the pulse.This patch is then aligned to a subsequent pulse.At the location where the patch best matches the following pulse, the exact delay between the placed patch and the previous pulse becomes the newly found start of the pulse.Diastole-patching operates under the assumption that the morphology and timing of the heart beat and the resulting PPG wave are comparable.Diastole-patching was originally applied to determine pulse wave velocity, interbeat interval, and the start of a pulse for PAT calculation (Bachler et al 2013, Hemon andPhillips 2016).However, the algorithm was only tested on synthetic data with a sampling frequency of 1 kHz, and on samples of the MIMIC Database (Moody and Mark 1996), which contains PPG recorded at a sampling frequency of 125 Hz.
The application of diastole-patching with lower sampling frequency recordings is hampered by a high sensitivity to noise due to the limited number of samples in the patch.Indeed, a small number of samples entails that even a minor disturbance to one of the values within the patch has a relatively large effect on the subsequent matching.The low sampling frequency means a low resolution in the time domain, as the time between two consecutive samples is increased.As a consequence, in lower frequency environments, the limitation of aligning two patches over integer samples leads to a limited time resolution, and hence inaccuracy, in the resulting PAT.
The aim of this paper is to develop a new strategy to estimate the PAT at sampling frequencies as low as 25 Hz.In this work, we built on diastole-patching by matching a whole PPG pulse to a template.This template is created by aligning all pulses within a segment of the signal and taking the average.Each individual pulse is then aligned to this template through a parameterized transformation.By taking a whole pulse into account and allowing for subsample alignment, as well as using a template instead of just the subsequent pulse, we hypothesize that the proposed method will allow for more accurate PAT estimation when going to lower sampling frequencies, which are commonly used in commercial products.

Data and preprocessing
The method has been developed and tested on the pulse transit time PPG Dataset (Mehrgardt et al 2022), publicly available on PhysioNet (Goldberger et al 2000).The dataset consists of data gathered from 22 healthy participants (of which 6 female) in an age range between 20 and 53 (with a mean of 28.52 years old), who performed three activities (sitting, walking, and running) in random order, for a total of 1 h per participant.During these activities, participants had both PPG and ECG recorded.Every signal was measured at a sampling frequency of 500 Hz.The signals were synchronized and the R-peaks in the ECG were manually annotated.
The PAT was calculated after filtering the original PPG signal using a band-pass filter between 0.07 and 35 Hz, as described by Peláez-Coca et al (2022).The PAT is defined as the time between the R-peak in the ECG and the start of the corresponding pulse in the PPG signal.The location of the R-peak was taken from the manual annotation, while the corresponding start of the PPG pulse was determined as the minimum value before the rising flank of the pulse, also known as the diastolic minimum.The rising flank corresponds to the maximum value of the first derivative of the PPG signal.The PAT obtained from the filtered signal recorded at 500 Hz is used as the reference PAT ref .
Next, the original PPG signal was downsampled to 250, 100, 50 and 25 Hz by applying an anti-aliasing filter and then decimation.This downsampled signal mimics the signal acquired at a lower sampling frequency, and is used as the input of the developed algorithm.Figure 1 shows an overview of the developed algorithm.

Obtaining pulses
The first step of the proposed method is to partition the PPG signal into individual pulses.Here, each pulse is defined as the PPG signal between two consecutive R-peaks in the ECG signal.The method could also be applied through the use of a rudimentary peak finder, as long as the location within the pulse is consistent.
The starts of these pulses in the PPG signal can then be used in combination with the location of the corresponding R-peaks to determine the PAT of the lower frequency signal.

Template creation
The PPG signal is split into 5 min, non-overlapping segments.For each of the segments, a template is computed from the pulses within the segment.The segment length was chosen as a trade-off between the inclusion of enough data to create the template while minimizing the influences of non-stationary effects.
The template was computed by averaging 80 % of the pulses in the segment.The 20 % of pulses with a higher average difference from the mean of all the pulses is excluded from the template in order to filter erroneous pulses.Since pulses have varying lengths, the tail of the template has less individual values to average.Because of this, the values of the template are calculated only if there are at least 11 values to average.

Calculating the parameters
To determine the PAT, we need to determine the delay between the R-peak of the ECG and the onset of the PPG pulse.To detect the onset of the PPG pulse in an accurate way, we propose to first align each individual PPG pulse in the segment, with the template pulse.The time shift needed in this alignment can then be used to determine an accurate PAT.In the alignment between an individual PPG pulse and template, a transformation based on three parameters is applied, which includes a slope adjustment θ 1 , a baseline shift θ 2 , and a time shift θ 3 .A graphical representation of each of these parameters can be seen in panel D of figure 1.
Figure 1.Overview of the method.In panel (A) the PPG signal with a low sampling frequency is represented.The signal is segmented in order to separate single pulses (panel B).The start of these pulses is determined and used to compute the 'raw' PAT (PAT raw ) (Panel G).In panel (C) the pulses from panel (B) are used to create the template (black dotted line).Then, in (D) the parameters of the transformation applied to each original pulse are estimated by matching the pulses (in gray) with the template to yield the adjusted pulse (in blue).The parameters of the transformation are θ 1 to adjust the slope, θ 2 to shift vertically, and θ 3 to shift horizontally, and their effects are shown individually in panel D. After obtaining the optimal parameters opt  q , the pulse is adjusted (panel E).In panel E, the gray line represents the original pulse, the blue line represents the adjusted pulse, and the black dashed line represents the template.The obtained adjusted pulse is used to calculate the adjusted PAT (PAT adj ) (panel F), which is then used in combination with PAT raw to compute the final PAT (PAT est ) (panel H).
In order to calculate the parameters θ we propose a parameterized model of the PPG pulse s(t) and template pulse S(t).In this parameterized model, the line segment s n (t) between two consecutive samples n − 1 and n is described as a linear function s n (t) = a n t + b n .This series of line segments is compared to a similar parameterization of the template S: S n (t) = A n t + B n .
The estimation of the transformation parameters  q can now be defined as an optimization problem, in which the cost function is defined as the sum of the squared area between the template and each transformed pulse: where ε is the error to be minimized, T s is the sampling interval, S n (t) is the line segment of the template, s n (t) is the line segment of the individual pulse, and | ( ) is a correction function applied to s n (t) in order to match the template segment S n (t).The goal of this correction function is to adjust for differences in the slope and baseline between the individual pulse and the template pulse.As such, we propose to use a first-order correcting function f (t) = θ 1 t + θ 2 .The use of θ 3 allows a subsample shift that is applied to the series of line segments s n (t).Effectively, this allows to perform a time alignment between the individual pulses and the template that has a time resolution that is better than the sampling frequency.The error function is globally minimized for [ ] , , 1 2 3   q q q q = with the optimal parameter opt  q defined as Because of the definition of s(t), S(t) and f (t) as a series of linear segments, the final system of equations for the optimal value of  q is a linear system, with three equations and three variables, and has one unique solution opt  q .
The optimized set of parameters opt  q defines the optimal correction to be applied to the pulse s(t) in order to align it with the template pulse.

Adjusting the pulses
As mentioned before, the correction function applies a slope adjustment and baseline shift to each individual pulse with f (t) = θ 1,opt t + θ 2,opt , where t = 0 is the start of the pulse.The time shift is applied by delaying the samples in time by the factor θ 3,opt .This means the updated pulse is given as where ( ) s t new is the adjusted PPG wave, while ( ) s t is the raw PPG wave.Even after the correction, there can still be a considerable difference between the adjusted pulse and the template.This difference indicates either a pulse that does not fit the standard shape, or that the segment is influenced by too many artifacts, causing the template to not be fully representative of the morphology of the PPG pulses in the segment.Therefore, the error ε obtained from equation (1) after having adjusted the pulses can be used as an indication of how likely the new estimate is correct.
To reduce the effect of an possibly erroneous adjustment on the final PAT, the original PAT raw is considered as well.The final estimated PAT is calculated as where PAT adj and PAT raw are the PAT calculated using the adjusted and the raw pulses, respectively, and η is the error factor defined as where ε is the error of the current pulse, max e is the maximum error for the current segment, and x is a scaling factor for the error.This scaling factor is introduced, as it often occurs that the maximum error is an order of magnitude larger than the other errors.This factor has been empirically set to be equal to 5, as it provides a tradeoff between the difference in magnitudes of error and bias towards the newly obtained PAT.
The final estimated PAT is thus a weighted average of the newly calculated PAT adj and PAT raw , where the weight is dependent on the remaining error ε.

Methodology for evaluation
To determine the accuracy of the PAT at lower sampling frequencies with respect to the reference (i.e. the PAT determined at high sampling frequency), the root mean squared error (RMSE) between the reference and the calculated PAT was determined for each subject in the dataset.The normality of distributions is checked with the Shapiro-Wilk test.The difference between distributions is tested by using the one-sided Wilcoxon signed-rank test.To check whether the proposed method reduces the error in PAT estimation, the alternative hypothesis states that the error of the final estimated PAT (PAT est ) is lower than the error of the raw PAT (PAT raw ).

Results
The results obtained on the signals acquired while sitting, walking, and running are shown in figures 2(a)-(c), respectively.To determine whether the method significantly improves the results, the one-sided Wilcoxon signed rank test was used, where the alternative hypothesis states that the error of the estimated PAT is lower than the raw PAT, as described in section 2. The Wilcoxon signed-rank test was used as all distributions were not normally distributed (Shapiro-Wilk test, p < 0.01), except the adjusted error for the sitting activity at a sampling frequency of 25 Hz, which was not significantly different from a normal distribution (p = 0.243).

Discussion
The sampling frequency is an important parameter that influences the accuracy of the estimation of crucial vital signs such as the PAT.In this work, we proposed a method to improve the acquisition accuracy of the PAT when the sampling frequency is as low as 25 Hz.
The results on the raw data (i.e.before applying our proposed method) show that the error in PAT increases when the sampling frequency lowers.This is due to a loss of information caused by the lower sampling frequency, as well as the lower resolution in time.The walking and running activities also show an increase in error in the PAT, which is in accordance with the increased prevalence and severity of movement artifacts.
The same behavior as the raw data can also be seen when applying our proposed method; the error in PAT increases when the sampling frequency lowers.This is to be expected, as there still is a loss of information due to the lower sampling frequency.However, this increase in error is smaller than for the raw data.Furthermore, the differing types of activity have a strong effect on the obtained error.
By comparing different activities, both the mean and the standard deviation of the error during the sitting and the running datasets decreased for all sampling frequencies.On average, the mean and standard deviation for the sitting activity decreased by 22.2 % and 48.8 %, respectively.For the walking dataset the mean decreased by 2.9 %, while the standard deviation increased by 4.2 %.Lastly, for the running activity, the mean and standard deviation decreased by 24.6 and 16.0%.Combining this, the method improves the error between the reference PAT and the estimated PAT by an average of 16.6 and 20.2 %, for the mean and the standard deviation respectively.
As can be seen in the results, the PAT calculated on the dataset acquired during walking shows a much larger RMSE than the other two datasets.This is likely due to the lower signal quality of the measurements during this activity.One cause of this could be the rhythmic effect of walking having close to the same frequency as the heart beat.The artifacts caused by the movement are therefore consistent in timing and shape, to the point where outlier removal does not function as expected.Another explanation could be that the activity caused an increase in movement between the sensor and the skin, further reducing the acquired signal quality.Applying the method described in this paper therefore also marginally improves the estimation of the PAT in this case and has less significance than the other two activities.
With our method we can improve the accuracy in PAT estimation for low sampling frequency recordings significantly, possibly allowing PAT analysis in common settings through unobtrusive monitoring.In case a high accuracy in PAT estimation is required, like in critical situations, a high sampling frequency could still be warranted, as suggested by other researchers (Béres andHejjel 2021, Peláez-Coca et al 2022).
One of the main limitations of the proposed method is the calculation of the reference template.The current approach to create the template performs adequately under the assumption that the window is clean of artifacts.The performance of the proposed method depends on the quality of the template and the data that is used to create the template.The performance reduces when artifacts such as motion are present in the data.An example of this can be seen in the results for the walking activity.The method still improves the estimation, but not as strongly as for the other two activities in the dataset.This limitation also imposes the assumption that the underlying heart activity, and hence the shape of the resulting pulse, is relatively constant within the window segment.
One possible way to mitigate this limitation is by adjusting the window size depending on the specific application.A larger window means that more pulses can be included to create the template, leading to higher quality of the template, but a larger window also means the template cannot account for possibly present disruptions or slow changes in the signal.Additionally, using a different method for a more robust template creation under certain conditions, such as the methods described by Papini et al (2017) or Margarito et al (2016) could improve the robustness of the proposed method, allowing it to be used in those circumstances.
Finally, it should be noted that the dataset used for this work contains only 22 subjects.While the results of this work are shown to be statistically significant, the findings can be solidified by applying the method to more subjects, with both healthy and diseased heart rhythm behavior, to observe the differences in results.
Another avenue to pursue in the future is the translation to blood pressure estimation, and the effect and behavior of the pre-ejection period during measurement, as well as the effect it has on this method.
The proposed method shows promise in allowing an accurate estimation of PAT even in low sampling frequency recordings.Our strategy opens new possibilities for continuous monitoring of health status through wearable devices.One application could be monitoring during sleep, as the improved accuracy during unobtrusive monitoring could allow more accurate estimation of vital information.In situations like these, our strategy can support increasing the accuracy of the acquisition of the PAT, especially when a lower sampling frequency is being used.

Figure 2 .
Figure 2. Box plots of the error distributions for different sampling frequencies, for the different activities.The whiskers extend to the 5th and 95th percentile, respectively.Both the error from the raw ('PAT raw ') and the estimated ('PAT est ') signal are shown.Significance of the one-sided Wilcoxon rank test is indicated where ' * ' is p < 0.05, ' ** ' is p < 0.01 and ' *** ' is p < 0.001.The blue diamond indicates the mean value.