Figure 1 shows the measurement flow. Each volunteer wore a SENSIO smartwatch recording wrist PPG and ECG during the experimental period. For volunteers in the health management center, 3 rounds of measurements were conducted. For volunteers in the outpatient clinic, 5 rounds of measurements were made. In each round, the participants maintained the sitting position, and ECG was measured in the first minute. Blood pressures were then measured by the sphygmomanometer on the other arm (not wearing the smartwatch) with the cuff aligned at the heart level. A one-minute rest was reserved between 2 adjacent rounds. The wrist PPG signals were continuously recorded throughout the course. In the end, baPWV was measured by the OMRON noninvasive vascular screening device, with the cuffs on 4 limbs in the supine position.

The experiment was approved by the research ethics committee of National Taiwan University Hospital (number 201902087RIPA). All data were collected in accordance with the approved protocol. Importantly, the dataset used in this study did not contain any personally identifiable information, and all records were fully anonymized prior to analysis. Informed consent was obtained from all participants, and the study was conducted in compliance with the ethical standards set forth in the Declaration of Helsinki and relevant national regulations.

The signal-processing flow is indicated in Figure 2. The PPG and ECG, sampled at 256 Hz, were extracted from the first minute of each round in the synchronization phase (Figures S1 A and S1 B in Multimedia Appendix 1). In the preprocessing phase, baseline wandering of signals was corrected by the discrete wavelet transform, and the 60-Hz power interference was suppressed by the notch filter. The amplitude of the whole signal segment was then normalized to [−1, +1]. The R peak of ECG and the valley of PPG signals were detected to calculate cycle length (Figures S1 C and S1 D in Multimedia Appendix 1). The skewness and variation of ECG and PPG cycle lengths were adopted to establish the signal quality index to exclude suboptimal ECG or PPG cycles for feature extraction. The first-order derivative PPG (FDPPG) and the second-order derivative PPG (SDPPG) signals were calculated. The systolic peak, notch, and diastolic peak were marked by the algorithm [19] for each PPG cycle (Figure 3A). The maximal slope (max slope) of the ascending systolic pulse, corresponding to the maximal rate of blood volume change, was identified by the first local maximum in FDPPG (Figure 3B) [20]. The local extrema of the SDPPG in systole are defined as a, b, c, and d points, where points a and c are local maxima and points b and d are local minima (Figure 3C) [21]. Point e is the local maximum around the boundary of systole and diastole in SDPPG. Point f is the first local minimum after point e.

The PPG pulse is regarded as a summation of several component waves, including the forward waves by left ventricular contraction and the distally reflected waves due to aortic elasticity and reservoir property [22]. The pulse decomposition analysis helps segregate the component waves [23]. With proper weighting, the variation of component waves can be reduced [24]. Five Gaussian waves are used for synthesizing the PPG pulse. Given θi=αi, βi, γi corresponding to pulse amplitude, pulse position, and pulse width of the component wave i, and Θ={θ1, θ2, …,θ5}, the summation of the Gaussian waves takes the form of

with

Denote G_i as the component wave described by g(t|θi). Given the boundary constraints, Lαi≤αi≤Uαi, Lβi≤βi≤Uβi, and Lγi≤γi≤Uγi [24], the interior-point method is used to solve the following optimization problem,

where w(n) is the weight to emphasize the informative portion of the PPG pulse sn with length M and is given by

Variables na and nf refer to the position of points a and f. The weight ω is set to 80 for stabilizing the variation of component waves in the sequence with acceptable mean square error between the synthesized waveform and original waveform.

Once the component waves are acquired, the forward wave is generated by combining G₁ and G₂. The systolic wave and diastolic wave are derived by combining G₁ to G₃ and G₄ to G₅, respectively. The respective peaks of the synthesized forward wave, systolic wave, and diastolic wave are named as pf, ps, and pd. In the following, the amplitude and position of feature x in the PPG pulse are indicated by Ax and nx, respectively. The amplitude of feature x in the ith-order derivative PPG is represented by Ax(i). The result of decomposed component waves by weighted pulse decomposition (WPD) is shown in Figure 4.

To assess the quality of WPD, WPD signal quality index, which was defined as mean square error between the PPG pulse, s(n), and the synthesized pulse, G(nTs|Θ), of >2×10-3, was implemented to remove disqualified pulses.

A total of 22 features were derived from the PPG pulse, FDPPG, and SDPPG (Table S1 A in Multimedia Appendix 2). The age index, which has been shown to be correlated with the augmentation index of aortic pressure [21,25],

and its related variant combining only highly correlated components,

were also used. There were 27 features derived from WPD (Table S1 B in Multimedia Appendix 2). The stiffness index (SI) is defined as the time interval between the peaks of systolic and diastolic waves [23] and is denoted by npd-nps. The time intervals of the third or fourth component wave to the forward wave were also calculated. Note that nps and npd were obtained from synthesized systolic wave peak ps and diastolic wave peak pd of WPD as shown in Figure 4 while nsys and ndia were marked as the positions of systolic peak and diastolic peak in PPG as shown in Figure 3.

The ECG-related features were also adopted (Table S1 C in Multimedia Appendix 2). The R peak and T peak of the ECG waveform were identified and marked as nR and nT. Since the R peak occurs earlier than the PPG valley of the same heartbeat, nR is negative in number. The pulse arrival time (PAT) measures the time span between R peak and PPG valley, denoted by -nR. PAT² and Height²/PAT² were included since either linear or nonlinear relationship between BP and pulse transit time has been shown [26]. The time span from R peak to maximum slope, peak of systolic wave, or component wave 2 was also considered.

Basic information (Table S1 D in Multimedia Appendix 2) contains age, height (H), weight, BMI, and lengths from arm to wrist (Law) and finger (Laf). The lengths from heart to brachium (Lb) and from heart to ankle (La) can be approximated by [27]

The length difference between ankle and brachium could be expressed by La-Lb.

Feature normalization is often adopted since the relative change of 2 features could provide additional information than each feature alone. To systematically derive the normalization results, we generate combined features by dividing the value of feature u by value of feature v. The combined features contain not only magnitude-normalized or time-normalized features but also basic information features.

The differences between the estimated results v^j and the measured PWV vj of the jth measurement are shown by the mean absolute error, mean error, SD, and root-mean-square error (RMSE), which are defined as follows.

The correlation coefficients together with P values are also provided. Since some participants have more than 1 measurement, to avoid unbalanced weighting, averaged PWV estimation and averaged PWV measurement are used for the statistical results per participant.

In this study, we used wrist PPG and ECG signals to estimate baPWV. The morphology of wrist PPG signals is quite different from that of finger PPG signals. The conventional approach that used finger PPG morphology features may encounter the problem of feature missing due to much fewer identifiable features of wrist PPG signals. In addition, the multivariable regression model used in prior works may be too simple to describe the complicated hemodynamic state in the vessels. Hence, we resorted to the machine learning algorithm to deal with the estimation. Although the wrist PPG and ECG signals were acquired before the baPWV measurement, they are still related to the vessel condition and stiffness. To further improve and refine the estimation results, hierarchical regression was adopted to shrink the range handled in the submodel. The achieved RMSE and SD by our hierarchical regression models for both men and women are lower than the threshold (150 cm per second) of acceptable accuracy for PWV estimation set by the ARTERY Society [32].

With the WPD and feature imputation techniques developed by us, more than 98% of all ambiguous and missing features of wrist PPG can be identified [19]. From the correlation results (Multimedia Appendix 4), besides age (feature 23) and age square (feature 63), correlation related to SDPPG amplitude of point c (feature 18), point d (feature 19), and point e (feature 20) are still obvious as what has been mentioned in finger PPG [25]. In addition, SI (npd-nps; feature 51), which are often missing in the original wrist PPG pulses, can be computed through the synthesized systolic and diastolic waves in decomposed wrist PPG. According to the feature importance (Multimedia Appendix 6), it still plays an important role for PWV estimation.

The multivariable regression uses only a few features. If significantly high correlations of those features to baPWV do not appear, the performance of estimation will be degraded. However, the machine learning algorithm can help exploit more linear or nonlinear information embedded in the PPG waveform or its component waves and thus is suitable for these applications. Furthermore, the combined features from PPG and ECG morphology, WPD, and basic information supplied more feature information sources that can be selected by the model.

The concept of hierarchical regression is to introduce different models to refine the estimation results. However, the global classifier or regressor must provide sufficiently correct classification to avoid model mismatch. From the hierarchical regression results, it is clear that the inclusion of overlapping zone in local regressors indeed improved the estimation results, as reflected in the improved correlation coefficients (Table 2). However, the determination of optimal range of overlapping zone is still controversial. If the overlapping zone is too wide, the hierarchical regression model would become similar to the general regression model. On the other hand, if the overlapping zone is too narrow, the misclassified data cannot be properly handled. In this study, we recommend the overlapping zone of 400 cm per second of 2 subdivision models because the misclassified data are near the boundary due to good capability of the global classifier and can be appropriately covered by the submodel. We conducted further analysis on the features that were misclassified for those samples not near the decision boundary. The results showed no significant outliers. Additionally, the vote counts for 2 classes across the entire forest were close, indicating low confidence among the trees. The latent properties beyond the observed features should be further studied. On the other hand, we also applied a Kernel Density Estimation–based mutual information analysis [33] to assess the relevance of individual features in male and female datasets. The mutual information values from male features were lower than those from female features, which can also explain the lower classification accuracy in male participants of our dataset.

This study has limitations, which point to the directions for future research. First, the sample size remained small and more older adult people were recruited in the study, which might limit its applicability in younger populations. While the current dataset demonstrates feasibility in estimating PWV using wrist PPG in older individuals, the skewed dataset toward older individuals may have influenced the performance due to age-related vascular characteristics. In future work, we plan to expand the study population by actively recruiting more young participants. The inclusion of younger participants will help balance the age distribution and allow for more robust assessment of the model performance across different age groups. This extension will not only improve the generalizability of the model but also enable a more comprehensive evaluation of age-related vascular changes. Second, the current model adopts machine learning algorithms to exploit linear and nonlinear features within the scope of this dataset. As the dataset grows in size and diversity, other deep learning algorithms, such as Bayesian neural networks or multilayer perceptrons, can be applied, which may offer better uncertainty quantification or modeling capabilities. Third, the feature space used in the current model is relatively high-dimensional, which may hinder its practical deployment on wearable or edge devices with limited computational resources. Feature compression or dimensionality reduction techniques can be considered to decrease model complexity in the future. This optimization will help make the system more suitable for real-time, low-power applications in wearable health care settings. Together, these improvements aim to enhance both the robustness and the applicability of the proposed approach, facilitating its transition toward practical use in diverse and real-world scenarios.