Filter and Sampling Rate Optimization for PPG-Based Detection of Autonomic Dysfunction: An ECG-guided Approach

In this study, Fig. 1 shows the environment and the instruments for signal acquisition. The physiological signal acquisition instruments, the PowerLab data acquisition device from ADInstruments [39], were used to acquire the ECG and PPG signals. This device is engineered for precise, consistent and reliable data acquisition and is connected to computers via USB links. With the support of the LabChart software package from ADInstruments [31], it can immediately display the acquired signals and store them as files for later analysis and reading. In this study, both signal sampling rates for ECG and PPG were set to 1000 Hz in the PowerLab system. It is worth noting that the resulting PPG signal with a 1000 Hz sampling rate is treated as the original signal to generate a downscaled PPG signal with a lower sampling rate. A lower sampling rate signal generation scheme is described below to obtain the lower sampling rates for the later experimental procedure in this study.

To simulate lower sampling rates from the original PPG signal (sampled at 1 kHz), we applied a uniform down-sampling method with a parameterized sampling period (per) defined in milliseconds.

Let:

PPGraw[n] be the raw PPG signal sampled at 1 kHz, where n is the signal length;

The pseudocode below describes the process: # Inputs: # PPGraw: raw PPG signal sampled at 1 kHz # per: sampling period in milliseconds # Output: # S: down-sampled PPG signal j = 240 # Number of seconds (4 minutes) S = [] for i in range(0, j * (1000 // per)): S.append(PPGraw[i * per])

This procedure was implemented in Python, using basic list indexing for down-sampling. The parameter per directly determines the sampling rate (e.g., per = 1 ms = 1000 Hz, per = 10 ms = 100 Hz, etc.).

In this study, the following sampling periods were used to generate lower-rate signals:

Sampling_periods = {2, 4, 5, 8, 10, 20, 25, 40, 50}, corresponding to frequencies between 20 Hz and 500 Hz. The lowest assumed sampling rate is 20 Hz for the application of video-based PPG [22].

After the acquisition of ECG and PPG signals, the signal processing flow, shown in Fig. 3, consists of four stages, namely preprocessing, feature finding, signal decomposition and feature calculation, and performance measurement.

Fig. 3.

Signal processing flow chart for ECG and PPG.

It should be noted that PPG signals are more susceptible to interference from baseline wander noise than ECG signals [8]. As in Stage 1 of Fig. 3, a baseline wander noise removal algorithm is developed to remove the baseline wander noise that was present in the raw ECG and PPG signals during preprocessing. To find the features of the signals after removing the baseline wander noise, the Pan-Tomkins algorithm was applied to localize the R peaks (marked as R_i in Fig. 2) in the cardiac cycles of the denoised ECG (Stage 2). At the same time, different filters were used to remove the PPG signals to localize the S peaks (marked as S_i in Fig. 2) in each cardiac cycle (Stage 2). Then, the two time domain features RRIi and SSIi (see Fig. 2) are calculated for ECG and PPG to determine the temporal length of each cardiac cycle in Stage 3. Finally, performance measurements of the variation of interbeat intervals and the associated HRV parameters were performed. The acquired parameters are the RRIV and HRV of the ECG signal at 1000 Hz sampling rate, and the SSIV and associated HRV of the PPG signal at different down-sampling rates in Stage 4.

To remove the raw baseline wander in the ECG and PPG data, a baseline wander removal scheme based on the work of Chazal et al. [32] was considered. In two steps, the scheme uses a median filter [33] to filter out the core signals of ECG and PPG from the raw data. Chazal’s method generates a baseline drift signal and uses it to offset the baseline wander with the raw data. In this way, a set of clear ECG or PPG signals is generated. The scheme is briefly described below.

To remove the baseline wander from the ECG and PPG signals, we applied a two-stage median filtering process as follows:

Step 1: A median filter is applied to the raw signal using a shorter window size (1^st_WinSize), resulting in the first filtered output.

Due to the different characteristics of ECG and PPG signals, different window sizes are used:

For ECG:

1^st_WinSize corresponds to 5 Hz (QRS wave)

2^nd_WinSize corresponds to 1.67 Hz (P and T waves)

This approach effectively suppresses the low-frequency drift while preserving the main features of the signal. These values were initially taken from the literature [40], [41] and then tested and fine-tuned to achieve the best results.

For the detection of R peaks, a real-time QRS detection algorithm proposed by Pan and Tompkins is used to determine the temporal locations of all R peaks [34]. Since the Pan-Tompkins method is very efficient and powerful in detecting R peaks in baseline wander-free ECG signals, no further signal processing is required before detecting R peaks. Fig. 4 illustrates the process of extracting QRS from raw signals. After squaring the signal, all negative values were eliminated. Based on the above results for detecting the temporal location of each R peak, the interval of each pair of neighboring R peaks (i.e., RRI) can be determined, and thus, the RRIV variation (i.e., RRIV) can also be obtained.

Fig. 4.

The process of extracting QRS from raw signals.

Since PPG signals are susceptible to motion artifacts, it is conceivable that PPG signals require careful signal processing to accurately detect PPG S peaks. In Table 1, four classical filters (Butterworth, Bessel, Chebyshev, and Elliptic), three smoothing filters (Savizky-Golay, average, and periodic moving average), and three wavelet-based filtering methods (Daubechies db6 (DWT), Coiflet C3 (DWT), and Morlet (CWT)) are tried to find the most appropriate PPG S peak for each cardiac cycle. The classical filters (F01 to F04) have some outstanding features in signal filtering, but disadvantages in roll-off rate, pass/stopband ripple, and phase response [10]–[13]. Due to their low implementation complexity and low latency, they were used to determine the PPG S peak in this study. Then, the three smoothing filters (F11∼F13) are considered since the PPG signal is not similar to the embedded ECG with the sharp QRS complex. These filters can provide different experimental aspects in detecting the PPG S peak. The wavelet-based methods (F21 to F23) include discrete and continuous wavelet transforms. For the DWT, Daubechies db6 and Coiflet C3 were used to perform a test for PPG signal processing, and Morlet was chosen for the CWT. Since these WT techniques have been used for ECG signal processing [35], it seems interesting to apply these methods to the PPG for investigation. The parameter settings required for the above filters are given in the description in Table 1.

To detect S peaks in the filtered PPG signal, a Python-coded peak search program was developed using a specific threshold value. Fig. 5 shows the results of the search for the ECG and PPG peaks.

Fig. 5.

Detected ECG and PPG peaks.

After peak detection processing, as shown in Fig. 5, the continuous cardiac cycles in the ECG or PPG are decomposed and stored in an array P(•), which records a series of peaks. Then the interval between each pair of neighboring peaks, i.e., an RRI or SSI, can be calculated using (1). (1) Inter_Beat_Intervali=Pi+1−Pi Inter\_Beat\_Interva{l_i} = P\left( {i + 1} \right) - P\left( i \right) where i ∈ Z⁺, the set of positive integers.

Based on (1), let IB_Interval be the set of all Inter_Beat_Interval_i, e.g., {Inter_Beat_Interval_i | i ∈ Z⁺}. In this study, the method for calculating the interbeat variation in the ECG or PPG (i.e., the RRIV or SSIV) is shown in (2). (2) Variation=maxIB_Interval−minIB_IntervalmeanIB_Interval*100% Variation = {{\max \left( {IB\_Interval} \right) - \min \left( {IB\_Interval} \right)} \over {{\rm{mean}}\left( {IB\_Interval} \right)}}*100\% where the functions max (), min (), and mean () generate the maximum, minimum, and mean values for IB_Interval.

To evaluate the effectiveness of SSI as an alternative to RRI, we used correlation and accuracy as performance metrics in this study. First, the correlation coefficients could be used to examine the trend of positive, negative, or no correlation between two discrete data sets. Let {r₁,r₂,...,r_n} and {s₁,s₂,...,s_n} be two data sets consisting of RRIs and SSIs, respectively. To measure the strength of a linear relationship between ECG and PPG parameters based on RRI and SSI, the Pearson correlation coefficient R_rs 3) is used. Through the correlation analysis, we can examine the closeness of the RRIs to the different SSIs obtained by the aforementioned preprocessing schemes (Table 1). (3) Rrs=n∑risi−∑ri∑sin∑ri2−∑ri2n∑si2−∑si2 {R_{rs}} = {{n\sum {{r_i}{s_i}} - \sum {{r_i}} \sum {{s_i}} } \over {\sqrt {n\sum {r_i^2} - {{\left( {\sum {{r_i}} } \right)}^2}} \sqrt {n\sum {s_i^2 - {{\left( {\sum {{s_i}} } \right)}^2}} } }}

Second, we used the average absolute errors to evaluate the accuracy of the different SSIs from the PPG compared to the RRI from the ECG. After obtaining the RRI_i and SSI_i sets, i.e., {r₁,r₂,...,r_n} and {s₁,s₂,...,s_n}, the average absolute error (Abs_Err) was calculated. For the performance evaluation of the interbeat variation, we evaluate the accuracy of the different SSIVs with respect to the RRIV using Abs_Err.

To investigate the applicability of the PPG-based solution for telemedicine and clinical applications, we include two HRV time domain parameters, SDNN and NN50, in addition to the interbeat interval and interbeat variation in the PPG in the accuracy evaluation [3]. In [36], the definition of SDNN is the standard deviation of all normal-to-normal intervals, and NN50 is the number of pairs of adjacent NN intervals that differ by more than 50 ms in the entire recording. In our study, the inclusion of SDNN and NN50 was essential for evaluating the validity of PPG-based detection of cardiac autonomic function.

To examine the use of PPG in screening for ANS disorders, the following four aspects of the examination are highlighted. The ANS regulates vital involuntary functions such as heart rate, blood pressure, and respiration through its sympathetic and parasympathetic branches. ANS dysfunction is associated with conditions such as diabetic neuropathy, cardiovascular disease, and postural orthostatic tachycardia syndrome (POTS). Due to its close connection to cardiovascular regulation, the ANS can be effectively monitored by non-invasive methods such as PPG:

In this study, the RRI generated by the Pan-Tompkins algorithm is used as a criterion for comparison with the SSIs generated by the filters with different down-scaled sampling rates. To calculate a series of intervals in the ECG and PPG, (1) is applied to determine the RRI and different SSIs. Then, (2) is required to obtain the interbeat variation for the ECG and PPG series. Then, SDNN and NN50 were determined by RRI and the different SSIs to evaluate the validity of PPG-based detection of the cardiac autonomic function. To calculate the correlation coefficients, see (3) to evaluate the accuracy of using PPG as a substitute for ECG.

Since the physiological signals of the participants are recorded for four minutes in the resting phase and the deep breathing phase, each 4-minute signal unit is divided into the first and last 30-second signal units and three middle 1-minute signal units for analysis. For the sake of signal stability, the three units in the middle were used to calculate the evaluation parameters.

To assess the validity of the interbeat interval in the PPG using the ECG evaluation criteria, Table 2 and Table 3 show the correlation coefficient and absolute error rate between the RRI at 1000 Hz and an SSI at different sampling rates.

Fig. 6 shows that the IIR filter bank performs better than the other filter banks in both the resting phase and the deep breathing phase. The four IIR filters are very similar in both phases, with the exception that Butterworth (F01) is slightly better than the other three. At the same time, each correlation coefficient (CC) value in the four IIR filters of the two phases decreases linearly and slowly when the different sampling rates are examined. In Fig. 6, we can see that the IIR filter bank has very high CC values in both phases even at low sampling rates. This shows that the filter bank is suitable for wearables in healthcare.

Fig. 6.

HeatMap of correlation coefficients between the RRI at 1 kHz and the SSIs at various rates in resting (a), breathing phase (b).

For the other filters, F11 performs best in both phases at different sampling rates in the FIR filter bank, and its CC value is close to the CC values of the IIR filters. F23 outperforms the other two filters in the WT filter bank and its effectiveness is close to the IIR and F11 filters.

Fig. 8 shows the absolute error, measured in microseconds, for all filters. For the IIR filter bank, it can be seen that the absolute errors of the filters are very low. Especially, at low sampling rates, the errors in both phases remain relatively low. This means that their accuracy is good enough for wearable healthcare applications. For the FIR filter bank, F11 has the lowest absolute error, which is close to that of the IIF filter bank. For F12 and F13, the accuracy becomes very poor in both phases below 25 Hz and 200 Hz, respectively. Therefore, they are not suitable for healthcare applications, especially at lower sampling rates. For the WT filter bank, F23 performs better than the other two filters in both phases in terms of absolute error.

At the same time, the anomaly of F21 in absolute error at 1000 Hz can be seen in the table. This could be due to the unsuitable selection of the mother wavelets.

In Table 2, three parameters related to heart rhythm variability, SSIV, SDNN, and NN50, were used to evaluate the accuracy of using PPG instead of ECG. Since the calculation of the absolute errors of these parameters is performed using (4), the units of measurement for the absolute errors are expressed in percent (%) for SSIV and NN50 and in ms for SDNN.

The evaluation tests showed the average absolute errors of SSIV, SDNN, and NN50 from the filters F01, F11, and F23 at different sampling rates, related to the test phase in Table 2. The F01 IIR filter bank produces the most accurate values for SSIV, SDNN, and NN50 compared to all other filters, and F01 in particular achieves the best accuracy. In addition, F11 and F23 showed very low average absolute errors in their respective filter bank, as can be seen in Fig. 8.

Fig. 7.

Evaluation of the best SSIV assessment range of normal cardiac autonomic function with the ROC curve, (a) – the resting phase; (b) – the deep breathing phase.

Fig. 8.

HeatMap of average absolute errors between the RRI at 1000 Hz and the SSIs at different sampling rates in resting (a) and breathing phases (b).

In addition, F23 shows greater superiority than F11 as the accuracy values are better at the lower sampling rates of Table 2. This indicates that the WT filters have good potential for use in wearable healthcare applications.

Table 3 shows the RRIV assessment range for normal cardiac autonomic function during the resting and deep breathing phases of an outpatient neurological examination [31]. Due to the superiority of acquiring PPG signals over ECG, the idea of replacing the ECG data in Table 3 with the PPG data based on the previous analysis is very promising for remote healthcare monitoring.

To obtain a case of the most appropriate SSIV assessment range for normal cardiac autonomic function, it is necessary to find the best SSIV assessment range based on a set of optimal sensitivity and specificity values of the statistical test results. To give an example, the SSIV result with a 1000 Hz sampling rate from the Butterworth filter (F01) was used with the previous correlation and accuracy analysis.

To replace RRIV with SSIV in Table 3, the RRIV and SSIV values of all participants, totaling 244 cases, were screened and confirmed by the neurologist in our research team to determine whether their cardiac ANS was normal or abnormal. Subsequently, the diagnostic results are treated with RRIV as the ground truth and the SSIV results are the predictive values. Table 4 shows a confusion matrix based on the RRIV ground truth and SSIV substitution value, with the SSIV assessment range values for all age groups taken from Table 3. In this case, the calculated sensitivity and specificity, which are 86.27 % and 98.45 % for the resting phase and 84.62 % and 98.96 % for the deep breathing phase, are shown in the column ‘RRIV source’ in Table 5.

To obtain the most appropriate normal assessment range of the SSIV, an inward/outward scaling strategy using the ROC curve is attempted [37]. For a scale-out increment of 1, the scope of the RRIV assessment range (Table 3) is extended by decreasing the lower limit of the range by 1 and increasing the upper limit of the range by 1. For example, the normal assessment range in the resting phase for the 20 to 29 year olds is 12 ∼ 46 ms; if the normal assessment range is in a scale-out increment of 2, it becomes 10 ∼ 48 ms; and if the range is decreased inward by 4, it becomes 16 ∼ 42 ms. Apply the same inward/outward scaling strategy for other age groups. In Table 5, the term ‘Out+2’ means that the normal assessment range is in a scale-out increment of 2, and ‘In-4’ means that the range decreases inward by 4. Consequently, Table 5 shows seven possible prediction models with sensitivity and specificity for different increments and decrements. Note that the term ‘RRIV source’ denotes an increment or decrement of 0.

As a screening tool for daily life, the best normal assessment range of the SSIV should have the characteristic of high sensitivity and low (1 - specificity), resulting in most true positives being positive and few true negatives being positive. Based on the seven SSIV prediction models in Table 5, the ROC curves of the resting and deep breathing phases are outlined in Fig. 5. In Fig. 7, the sensitivity and specificity of the best cut-off point are 96 % and 92 % in the resting phase, and 94 % and 96 % in the deep breathing phase, respectively. It can be seen that the SSIV assessment range of ‘In-1’ in Table 5 comes closest to the best cut points for both phases.

In addition, Fig. 7 shows an evaluation index, namely the area under the curve (AUC). In the ROC curves, the AUC indicates the discrimination capacity of an SSIV prediction model. In the case of Fig. 7, the AUC of the resting and deep breathing phases are 0.982 and 0.979, respectively. If the value of the AUC is greater than 0.9, the SSIV prediction model has excellent discrimination. As a screening tool, this result is very accurate and sufficient. The most appropriate SSIV assessment range is shown in Table 6.