PPGTempStitch: A MATLAB Toolbox for Augmenting Annotated Photoplethsmogram Signals

An annotated photoplethysmogram (PPG) is required when evaluating PPG algorithms that have been developed to detect the onset and systolic peaks of PPG waveforms. However, few publicly accessible PPG datasets exist in which the onset and systolic peaks of the waveforms are annotated. Therefore, this study developed a MATLAB toolbox that stitches predetermined annotated PPGs in a random manner to generate a long, annotated PPG signal. With this toolbox, any combination of four annotated PPG templates that represent regular, irregular, fast rhythm, and noisy PPG waveforms can be stitched together to generate a long, annotated PPG. Furthermore, this toolbox can simulate real-life PPG signals by introducing different noise levels and PPG waveforms. The toolbox can implement two stitching methods: one based on the systolic peak and the other on the onset. Additionally, cubic spline interpolation is used to smooth the waveform around the stitching point, and a skewness index is used as a signal quality index to select the final signal output based on the stitching method used. The developed toolbox is free and open-source software, and a graphical user interface is provided. The method of synthesizing by stitching introduced in this paper is a data augmentation strategy that can help researchers significantly increase the size and diversity of annotated PPG signals available for training and testing different feature extraction algorithms.


Introduction
Photoplethysmography is a technology that optically detects changes in the blood volume of microvascular tissue beds. This technology, which can obtain a wealth of information about the cardiovascular system, has received extensive attention from scientists with different backgrounds, as it is non-invasive and can be continuously monitored. It is also used in devices like smartphones and wearable devices for health monitoring and primary health checks [1,2].
When analyzing photoplethysmogram (PPG) signals, researchers usually need to extract the features in the time domain, such as the locations and amplitudes of the PPG systolic peaks and the onsets [12]. Therefore, the accuracy of the feature extraction algorithm will affect the accuracy of the analysis results. To evaluate the performance of PPG time domain feature extraction algorithms, a large number of PPG signals with different sampling frequencies, noise levels, morphologies, and heart rhythms are required. However, few publicly accessible annotated PPG signals exist. The MIMIC database [13] contains 64 annotated PPG signal records, the length of which is one hour, and only the onsets of the PPG waveforms are annotated. Additionally, the IEEE TBME Respiratory Rate Benchmark dataset has 42 annotated PPG signal records, each eight minutes long [14], but only the systolic peaks are annotated in these records. Therefore, this study aimed to develop a toolbox to generate PPGs due to the high interest in research related to PPG signal analysis and app development.
Some mathematical models have been used to generate PPGs. In [15], a nonlinear model was used to generate reference PPG to clean PPGs before extracting features. Another method used the two Gaussian functions to model a single pulse to generate PPGs [16]; this method used the average of two autoregressive moving average models to approximate the parameters of the pulse model. Similarly, another method using a single pulse modeled by a log-normal and two Gaussian functions was applied in [17]. The beat-to-beat intervals were extracted from the RR interval in the electrocardiogram signal, and each PPG pulse was then connected according to the RR interval. Furthermore, the authors' previous work utilized a Two Gaussian Functions model based on circular motions to generate regular and irregular PPGs [18,19]. However, these methods [15][16][17][18][19] were focused on the morphology of the PPG, and no annotations were provided.
In this study, a MATLAB toolbox called "PPGTempStitch" was developed to stitch short-time annotated PPG templates together to generate a new long-time PPG signal for any sampling frequency. Stitching technology has been widely used in image processing [20][21][22]. To our knowledge, this present study is the second attempt to apply the stitching idea to time-series PPG signals. The first attempt to stitch time-series PPG signals was discussed in [23]. The "PPGTempStitch" toolbox generated by this study is an opensource package designed to use existing routines in MATLAB 2020a for reproducibility.

Materials and Methods
To evaluate the performance of PPGTempStitch when stitching different PPG types together, four types of normalized PPG templates were used for testing: regular PPG, irregular PPG, fast-rhythm PPG, and noisy PPG. Regular PPG occurs when the PPG rhythm is regular and the mean heart rate is within 59-99 beats per minute. Irregular PPG occurs when the mean heart rate is within 59-99 and the PPG rhythm is irregular. Fast-rhythm PPG occurs when the PPG rhythm is irregular but the mean heart rate is higher than 100. Lastly, noisy PPG occurs when the PPG segments contain noise. Figure 1 shows the dataset that was used in this study, these PPG templates came from different subjects in the MIMIC III database. The templates, which have already been normalized, were recorded at a 125 Hz sampling frequency. Figure 2 shows the flow chart for the steps used to stitch two PPGs (x 1 and x 2 ). These steps were normalization, stitching, interpolation, and method selection.

Normalization
PPGs may have different amplitudes. To enable comparison of different PPGs, they are usually normalized in pre-processing. In the present study, min-max and z-score normalization were both utilized to scale the two PPG segments before they were stitched together; this ensured that the amplitudes of the PPGs were consistent.
Min-max normalization was used to scale PPGs in the 0-1 range, as follows: where x is an original PPG signal, and y is the normalized PPG signal. The other normalization method, called the z-score, was used to center PPGs with a mean of 0 and scaled them to have a standard deviation of 1, as follows: where x is an original PPG signal; µ and σ are the mean and standard deviation of x, respectively; and y is the normalized PPG signal. In the current study's toolbox, the two normalization methods could not be used at the same time. The normalized PPGs were named "y 1 " and "y 2 ", which correspond to the first PPG (x 1 ) and second PPG (x 2 ), respectively.

Figure 1.
Four types of PPG templates. The regular, irregular, fast-rhythm, and noisy templates are (a-d), respectively. The "*" and "+" in (a-d) refer to the annotated systolic peak and onset, respectively.

Figure 2.
Flowchart for stitching two PPGs together. The "Peak" refers to the systolic peak.

Stitching
A major challenge for stitching two PPG signals together is choosing the stitching point, as the amplitude of the end point of the first PPG (y 1 ) is always different from that of the starting point of the second PPG (y 2 ). In this case, two signals were stitched together based on the systolic peak or the onset-the main two features of the PPG signal in the time domain. The difference between these two methods of selecting the stitching point was as follows: • Based on the systolic peak. For the first PPG, the stitching point was the last systolic peak, while the stitching point of the second PPG was its first systolic peak. • Based on the onset. For the first PPG, the stitching point was the last onset, while the stitching point of the second PPG was its first onset.
The two PPGs were then aligned in time based on the stitching point. The segment after the last onset of the first PPG and the segment before the first onset of the second PPG were discarded. Figure 3 compares stitching based on the systolic peak with stitching based on the onset.

Interpolation
The value of the stitching point in the first PPG (y 1 ) may be different from that in the second PPG (y 2 ). Therefore, to smooth the stitching point, cubic spline data interpolation was used [24]. Cubic spline interpolation involves a spline in which each piece is a thirddegree polynomial specified by its values and first derivatives at the end points of the corresponding domain interval.
The interpolation involved a total of 20 samples: the stitching point, nine samples before the stitching point, and 10 samples after the stitching point. To express these more clearly, the PPG segments involved in the interpolation were named "y t " and those that were used after the interpolation were named "z".
In this study, the first five and last five samples of y t were used to fit the interpolation function, and the values of the 10 other samples in the middle of y t were replaced by the values of the samples in z.

Stitching Beat Selection
As skewness is the optimal signal quality index (SQI) of PPGs [23], the skewness of the previous, stitched, and next beats were calculated. One PPG beat began with the onset of a beat and ended with the onset of the following beat. For the stitching method based on the systolic peak, the stitching beat was the beat where the stitching point was located. For the stitching method based on the onset, the stitching beat was the beat following the stitching point. The previous beat was the beat before the stitching beat, while the next beat was the beat following the stitching beat. To examine a PPG beat, the skewness was applied as follows: where µ and σ are the mean and standard deviation of the beat, respectively; z(i) is the value of the samples; and N is the number of samples in the beat. The skewness index was calculated as follows: where k previous , k stitching , and k next are the skewness of the previous, stitched, and next beats, respectively. In this study, "s 1 " was the skewness index of the stitching method based on the systolic peak, and "s 2 " was the skewness index of the stitching method based on the onset. The skewness index was used to evaluate the quality of stitching, and the stitching method with a higher skewness index was considered better than that with a lower skewness index. To improve the signal quality, the result of the stitching method with a higher skewness index was chosen as the final output in this study. Table 1 shows the skewness of the aforementioned three beats and the skewness index in different combinations. Figure 4 provides three examples of stitching results. Figure 4a,b are the results of stitching a regular PPG with another regular PPG based on the systolic peak and the onset, respectively; Figure 4c,d are the results of stitching a regular PPG with an irregular PPG based on the systolic peak and the onset, respectively; and Figure 4e,f are the results of stitching a fast-rhythm PPG with a noisy PPG based on the systolic peak and the onset, respectively. The skewness indices in Figure 4a,c,f are greater than those in Figure 4b,d,e, respectively. The differences in the signal quality of the PPG waveforms between Figure 4a,b and between Figure 4c,d are not very obvious, but the signal quality of the PPG waveform in Figure 4f is obviously better than that in Figure 4e. The length of the stitched PPG was less than the sum of the lengths of two PPGs. Table 1. Signal quality indices of the stitching results ("stitching beat" refers to a beat that has been stitched, "previous beat" to the beat prior to the stitching beat, and "next beat" to the beat after the stitching beat. The "signal quality index" refers to skewness of the beat). PPG templates were augmented to generate new PPG signals (see Figure 4), which can be helpful when the training PPG sample size is small (e.g., 10 PPG recordings of subjects with a specific disease). The introduced strategy stitched templates in a random fashion, similar to combining different images using the data augmentation strategies that are utilized by 2D image processing to generate new images [25]. The math behind data augmentation transformations in 2D signal processing is usually simple (e.g., scaling, translating). The present study opened this area of augmentation for 1D signal processing using PPG signals, and the results were impressive. To the authors' knowledge, this was the first study of 1D signal analysis to discuss stitching as an augmentation step to generate annotated time-series health data, such as PPGs. (e,f) are the results of stitching a noisy PPG with an irregular PPG based on the systolic peak and the onset, respectively. The stitching points are denoted as "*" in (a,c,e) and as "+" in (b,d,f). The cyan, red, and magenta parts are the previous, stitching, and next beats, respectively. The numbers above the curve represent the skewness of the beat. The "s 1 " and "s 2 " refer to the skewness index in the stitching methods based on the systolic peak and the onset, respectively.

Annotations
Annotations are the labeled PPG waveform, such as a systolic peak, a diastolic peak, an onset, and a dicrotic notch in the time domain [26,27]. For the toolbox developed in this study, the onset and systolic peak were supported. Additionally, the annotations could be merged after the PPGs were stitched together. For the first PPG, the annotations before the stitching point were reserved, while those after the stitching point were discarded. For the second PPG, the annotations before the stitching point were discarded, while those after the stitching point were merged with the reserved annotations in the first PPG according to their positions relative to the stitching point.
Since the interpolation step changed the values of the stitching point, the previous four samples, and the following five samples, the annotation at the stitching point needed to be corrected. For the stitching method based on the systolic peak, the systolic peak at the stitching point was corrected to the maximum value of the 10 samples after interpolation. For the stitching method based on the onset, the onset at the stitching point was corrected to the minimum value of the 10 samples after the interpolation.

Onset and Systolic Peak Detection Algorithms
To obtain the optimal onset and systolic peak extraction algorithm, three pulse wave extraction algorithms were compared, as follows: • Method I [28] used signal derivative, Hilbert transform, amplitude thresholding, and slope-reversal based approaches. The steps were as follows: -PPGs were filtered using a sixth-order Butterworth low-pass filter with a cutoff frequency of 15 Hz.

-
The Hilbert transform was applied to the second derivative of the PPG signal.

-
In the Hilbert transform data, the region greater than 50% amplitude was selected. The slope reversal points within these areas were determined as the maximum peaks. • Method II [29] first extracted the initial peaks, and classified them as true and false peaks. For each false peak, an algorithm was used to relocate the true peak. The steps were as follows: -Systolic peak detection: * The missing data points or data points with values greater than 20 times the PPG's median waveform height were defined as outliers. The outlier data points that lasted for less than 0.2 s were linearly interpolated, and the result was denoted as w0.
* w0 was linearly detrended, and its power spectrum density was computed using a fast Fourier transform. The maximum power spectrum in the 0.8-3.0 Hz range was the heartbeats' frequency, while the average beat-to-beat interval was its reciprocal value.

*
The waveform w0 was smoothed using a center median filter with a window size set to one-fifth of the estimated beat-to-beat interval, followed by a center moving-average filter with the same window size to generate w1. * w2 was generated by the smoothed w1 passed through a third-order, lowpass Butterworth filter with a cutoff frequency of one-and-a-half times the estimated heartbeat frequency.

*
The baseline b of the PPG waveform was generated by applying a center moving-average filter to waveform w2, with a window size set to 1.5 of the estimated beat-to-beat interval.
* To extract the systolic peaks, an initial set of peaks on w0 was used to find a local maximum in each time interval where the filtered waveform w2 was above the baseline b.
* All peaks on w0 were sorted by amplitude in increasing order and selected the amplitude value at the 2/3 length position. Each initial peak needed to satisfy two conditions: it needed to be greater than half of the selected amplitude value, and each peak-to-peak interval could not deviate from the median interval by more than two times, referred to as MAD. MAD was calculated as MAD = median[|t − median(t)|], where t is the peak-to-peak intervals. The PPG peaks that did not satisfy these two conditions were identified as potential false peaks, while the remaining were identified as true peaks.
* An attempt was made to relocate the putative false peaks and identify the locations of true ones. Specifically, starting from each interval's left-boundary position, median(t) seconds were added, and this point was considered as the expected position of the next peak. Next, the area around the expected position was searched to identify a local maximum on the smoothed waveform w1 within a window size of a length set to MAD. If the local maximum was located at either end of the window, the window size increased by MAD seconds, and the process was repeated. Then, the equivalent maximum was on w0 and labeled as the next peak. Starting from this newly discovered peak position, the procedure above was repeated until the end of the interval was reached and multiple consecutive missed peaks were recovered.

-Onset detection:
The onset corresponding to each true peak was detected in three steps, as follows: * Ranges were located where w0 was below both the filtered waveform w2 and the baseline b.
* If multiple ranges were identified, they were ranked based on their lengths. Then, the top two ranges' rightmost one was selected.

*
The minimum position on the waveform w0 in the selected range was considered the onset.
• Method III extracted the systolic peaks based on a block-based method [30], and the local minimum between two successive peaks was defined as the onset.

-Systolic peak detection:
* A second-order Butterworth 0.5-8 Hz bandpass filter was applied, and then the filter's output was clipped by keeping the signal above 0 to produce the x[n] signal.

*
The filtered signal was squared to emphasize the large difference between the systolic wave and diastolic wave in x[n]. This squaring step produced the y[n] signal. * Blocks of interest were generated based on two event-related moving averages and an offset threshold, as follows: where MA peak is the first moving average used to emphasize the systolic peak area, MA beat is the second moving average used to emphasize the beat area to be used as a threshold for the first moving average,ȳ is the mean of the squaring result, and W 1 , W 2 , and β are parameters. The blocks of interest were generated by comparing the MA peak signal with THR 1 , and the blocks wider than or equal to THR 2 were classified as the systolic peak area. The optimized parameters (W 1 , W 2 , β = 111 ms, 667 ms, 2) were used. * The maximum value in each block was considered the systolic peak.
-Onset detection: This study did not discuss the onset detection algorithm.

Performance Evaluation
Algorithms were evaluated using two statistical measures: Positive predictivity (PP) = TP TP + FP × 100%, where true positive (TP) is the number of features detected as features, false negative (FN) is the number of features that were not detected, and false positive (FP) is the number of non-features that were detected as features. In this study, each feature (systolic peak and onset) obtained by the algorithms was defined as a true positive if it deviated from the corresponding annotation within ±10 ms [31]. Figure 5 shows the main dialogue of the graphical user interface (GUI). Four types of PPG templates were supported in the GUI. These templates were extracted from different subjects in MIMIC III database. To synthesize an annotated PPG, users could complete all their work in the main menu. First, the signal length and sampling frequency of the synthetic PPG were modified. The PPG signal needed be longer than five seconds, the length of all the PPG templates, and so the stitched PPG was resampled to the required sampling rate when the needed sampling frequency was not 125 Hz. The second step was modifying the percentages of the different PPG types. The sum of all the percentages should be 100. When the percentage of a certain PPG multiplied by the signal length is less than 5 s, this signal type may not appear in the final synthesized PPG. The third step was selecting the normalization method. The toolbox provided two commonly used normalization methods: Min-Max and z-score normalization. When the "Generator" button was pressed, the GUI attempted to generate the synthetic PPG and showed it at the bottom of the dialogue. The toolbox generated PPGs using the following steps:

1.
Determine the proportion of beat types. The synthesized signal's length (n s ) was divided by the template's length 5 s to obtain the number of templates n at . Then, n at was multiplied by the percentage of each template to get the results n atp . The integer part of n atp was used as the number of each template. The four templates represented four beat types: regular, irregular, fast rhythm, and noisy.

2.
Randomly arrange the determined number of templates. Then, the method proposed in this study was used to stitch the determined templates in sequence one by one. When stitching was based on the systolic peak, the stitching point of the first PPG was the last annotated systolic peak, and the stitching point of the second PPG was the first annotated systolic peak. When stitching was based on the onset, the stitching point of the first PPG was the last annotated onset, and the stitching point of the second PPG was the first annotated onset.

3.
Adjust the signal length. If the length of the synthesized PPG obtained by Step 2 was less than the required length n, the difference between the lengths was defined as the new n s . Then, Steps 1 and 2 were repeated until the synthesized PPG was longer than n, at which point the first n seconds of the signal were considered as the signal output s p .

4.
Re-sample the signal. If the desired sampling frequency was not 125 Hz, s p was resampled to the desired frequency to generate the final output s end . The annotations were mapped to s end according to the ratio of the required sampling rate to 125. Each annotated systolic peak was corrected to the maximum position in the 0.01 s window centered on the peak. Likewise, the onset was corrected to the minimum position in the 0.01 s window centered on the onset.
After the signals are synthesized, users of the toolbox can press the "Save" button to save the synthetic PPG to a comma-separated values file (.csv), a text file (.txt), a generic data file (.dat), or an Excel spreadsheet file (.xls, .xlsm, or .xlsx). If more noise levels and waveform types are required, users can also add new PPG templates with annotations to the "Manage Templates" tab.
To maintain data consistency, users can only add five-second PPG templates, and the added templates will be resampled at 125 Hz when saving. When the signal length is short, the percentage of each template in the synthesized PPG obtained by this toolbox is not exactly the same as that set by the user.

Performance Evaluation Results
To evaluate a feature extraction algorithm, data containing different conditions are required. Typically, a large number of annotated PPGs are used to evaluate the algorithm's performance under different conditions. However, this toolbox could generate PPG signals under different conditions, thereby providing a solution to the current shortage of annotated PPGs.
To compare the three feature extraction algorithms, 10 PPG records were synthesized using this toolbox. These records contain different conditions. These records contained different conditions: records 1-5 were used to compare these methods' accuracy in different signal compositions at the same sampling rate and length, records 5-8 were used to compare these accuracies with different sampling rates and lengths under the same signal composition (e.g., the same proportion of beat types), and records 1, 7, and 10 compared these methods under unique conditions. Table 2 compares the application of the three feature extraction algorithms on different synthetic PPGs. Since Method III did not discuss the onset detection algorithm, its onset accuracy was not discussed. Record 1 was stitched only by the regular template. In this case, all three methods got SE = 100% and PP = 100% in systolic peak detection, and no assessment about which algorithm was best could be made. For onset detection, Method I got SE = 100% and PP = 100%. However, Method II got SE = 61.1% and PP = 61.1%. Figure 6a shows a segment of Record 1. Method II failed to detect the second and fifth onsets, and the deviation between the result of Method II and the annotation was greater than 10 ms. When the signal quality decreased, the accuracy of Methods I and II decreased, while Method III still achieved high accuracy.
Comparisons of the three algorithms' average accuracy revealed that Record 6 had the lowest accuracy. Figure 6b shows a segment of Record 6. Both Methods I and II missed the second systolic peak. Method 1 only selected the region greater than 50% amplitude for Hilbert transform data. This may have been due to the fact that the amplitude of the systolic peak may become lower than the threshold when noisy or irregular events occur. Method II judged the second peak as a false peak that was caused by the peak-to-peak interval threshold.
For the sixth systolic peak, the offset between the detection result of Method I and the annotation was greater than 10 ms. The results of Method III were not exactly the annotation location, but the offset was less than 10 ms. After comparing the average accuracy of 10 records, Method III was found to be the optimal method for detecting the systolic peak. Likewise, the accuracy of Method I was better than Method II for onset detection. Furthermore, by comparing the accuracy of these 10 records, the changes in sampling frequency were found to have no great influence on the three methods' accuracy.  Method I [28] Method II [29] Method III [30] Sampling

Discussion
After comparing the average of the three algorithms' accuracy on different records, the following ratio is recommended to test the feature extraction algorithm: the proportions of regular, irregular, fast rhythm, and noisy beats should be 50%, 20%, 20%, and 10%, respectively. Users of the toolbox can synthesize more PPGs with different sampling rates and different waveforms to evaluate the feature extraction algorithms. Because the order of the templates is random, the toolbox's PPG signal output may differ even if the same parameters (e.g., scale, sampling rate, signal length) are used.
To ensure the authenticity of the synthetic PPGs, all PPG templates came from real PPGs in the MIMIC III database. Noise was not added to the PPG templates. If noise is required, users can generate PPGs and then add noise. The advantage of this toolbox is that it can generate PPGs with annotations; previous studies [15][16][17][18][19] only generated the PPGs without annotation. Thus, by using this toolbox, users can generate PPGs to test and compare the developed feature extraction algorithms. Several PPG features exist in the time domain-onset, systolic peak, dicrotic notch, and diastolic peak. This study focused only on the main features of a PPG waveform: systolic peak and onset. However, other annotated features can also be generated using the same method, as discussed in [32,33].
One advantage of this toolbox is that the augmented PPGs are generated based on human data. Augmenting PPGs under the same condition can mislead the evaluation of the feature extraction algorithms. However, the combination of different data types covers different conditions to simulate real-life situations, reducing bias in performance assessment. One limitation of this toolbox is that only four types of PPG templates are included by default. If more PPG morphologies are needed, users can add their templates in the "ManageTemplates" tab of the toolbox.

Conclusions
PPGTempStitch, a new MATLAB toolbox for generating annotated PPG signals with any sampling frequency and a time length longer than 5 s, was described herein. The generated PPGs can contain different noise, waveform, and heart rhythm levels to simulate real-life conditions. Users can generate different annotated PPG types by adding new PPG templates. The developed toolbox is free and open-source software. Hopefully, the user-friendly toolbox will make PPG research easier, especially when evaluating PPG feature extraction algorithms.
Author Contributions: M.E. designed the study. Q.T., Z.C., C.M., R.W. and M.E. conceived the study, provided directions, feedback, and/or revised the manuscript. M.E. led the investigation and drafted the manuscript for submission with revisions and feedback from the contributing authors. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.