Photoplethysmography sampling frequency: pilot assessment of how low can we go to analyze pulse rate variability with reliability?

Pulse rate variability (PRV) analysis appears as the first alternative to heart rate variability analysis for wearable devices; however, there is a constraint on computational load and energy consumption for the limited system resources available to the devices. Considering that adjustment of the sampling frequency is one of the strategies for reducing computational load and power consumption, this study aimed to investigate the influence of sampling frequency (fs) on PRV analysis and to find the minimum sampling frequency while maintaining reliability. We generated 5000, 2500, 1000, 500, 250, 100, 50, 25, 20, 15, 10, 5 Hz down-sampled photoplethysmograms from 10 kHz-sampled PPGs and derived time- and frequency-domain variables of the PRV. These included AVNN, SDNN, SDSD, RMSSD, NN50, pNN50, total power, VLF, LF, HF, LF/HF, nLF and nHF for each down-sampled signal. Derived variables were compared with heart rate variability of the 10 kHz-sampled electrocardiograms, and then statistically investigated using one-way ANOVA test and Bland–Altman analysis. As a result, significant differences (P  <  0.05) were found for SDNN, SDSD, RMSSD, NN50, pNN50, TP, HF, LF/HF, nLF and nHF, but not for AVNN, VLF and LF. Based on the post hoc tests, it was found that the NN50 and pNN50, SDSD and RMSSD, LF/HF and nHF, SDNN, TP and nLF analysis had significant differences at fs  ⩽  20 Hz, fs  ⩽  15 Hz, fs  ⩽10 Hz; fs  =  5 Hz, respectively. In other words, a significant difference was not seen for any variable if the fs was greater than 25 Hz. Consequently, our pilot study suggests that analysis of variability in the time and frequency domain from pulse rate obtained through PPG may be potentially as reliable as that derived from the analysis of the electrocardiogram, provided that fs  ⩾25 Hz sampling frequency is used.

was not seen for any variable if the f s was greater than 25 Hz. Consequently, our pilot study suggests that analysis of variability in the time and frequency domain from pulse rate obtained through PPG may be potentially as reliable as that derived from the analysis of the electrocardiogram, provided that f s ⩾25 Hz sampling frequency is used.
Keywords: heart rate variability, photoplethysmography, pulse rate variability, sampling frequency (Some figures may appear in colour only in the online journal)

Introduction
Pulse rate variability (PRV) is an analytic method for assaying pulse-to-pulse intervals, and is usually derived from analyzing photoplethysmograms (PPGs). PRV is similar to heart rate variability (HRV) in that it is a method for gauging the autonomic nervous system in terms of analyzing the interbeat interval (IBI) generated from heart activity (Taskforce 1996). However, HRV and PRV analysis are different, because HRV is derived from an electrocardiogram (ECG) that is an electrical activity of the heart, and PRV is derived from PPG that is a mechanical phenomenon of pulsatile blood. Several research studies have shown that the PRV has significant coherence with the HRV in resting condition (Hayano et al 2005, Lu et al 2008, Gil et al 2010. Between PRV and HRV, however, there is significant divergence in dynamic conditions such as exercise (Carrasco et al 1998, McKinley et al 2003, Charlot et al 2009. As such, PRV is regarded as a surrogate for HRV under certain static conditions like resting, but more research is required to generalize it for practical use.
Nevertheless, the scenarios for using PRV have not been thoroughly investigated, even though numerous devices ready to adapt PRV for healthcare use. Representative examples are wristband-type wearable devices and smartphone applications. In both types of devices, PRV is usually applied to the evaluation of emotional status like stress or mood (Lakens 2013, Muaremi et al 2013. Many of recent wearable devices basically contain an optical sensor for PPG measurement, and the obtained PPG is usually used for calculating the pulse rate as a surrogate for the heart rate, with PRV then derived from the pulse rate. In case of a smartphone, camera-based PPG recording has been researched in several studies (Jonathan and Leahy 2010, Pelegris et al 2010, Gregoski et al 2012, Scully et al 2012, Heathers 2013. Although there has been some research on the clinical significance of PRV information from wearable devices or smartphone camera-based systems (Heathers 2013, Lakens 2013, Flatt and Esco 2015, Mahdiani et al 2015 a clear analysis on the utility of these methods has not been made. For camera-based PPG measurements, sampling frequency is one of the most important limitations of its use. Camera-based PPG recording usually depends on video recording or a consecutive screen capturing. However, the frame rate of recorded video is usually 29.97 frames per second (fps). Moreover, according to some manufacturers, data may be sampled at less than a 15 Hz rate. Considering that the sampling frequency for PRV analysis is around 100 Hz in clinical settings (Saeed et al 2011, Garde et al 2014, camera-based PPG sampling frequencies are too low and it could cause critical errors in PRV analysis. For wearable devices, the sampling rate is a key factor in optimizing performance. As the wearable platform has a very limited computing power and battery life, and at the same time, the power consumption of its optical sensor is very high (~1 mA), it is necessary to minimize sampling and computation as much as possible. With this in mind, finding the optimal or minimal frequency to maintain significance for the PRV reads is very important in PRV-based applications. Moreover, such research could be used for evaluating of the clinical significance or availability of the current PRV devices. Furthermore, the research could also provide guidelines for sampling frequency and power source design in implementing PPG devices.
In HRV research, several studies on the effect of the sampling rate have been conducted. Ziemssen et al generated 100 and 200 Hz-sampled ECGs from 500 Hz-sampled ECG, and then derived the HRV for comparison (Ziemssen et al 2008). They found that there were no statistically significant differences in baroreflex sensitivity and frequency bands ranging from VLF to HF using 100 Hz sampling instead of the original 500 Hz one. Only the UHF band (>0.40 Hz) was significantly different. In another study, Hejjel and Roth resampled tachograms with different ratios and compared the obtained HRV parameters (Hejjel and Roth 2004). They proposed a 1 kHz sampling as the optimum sampling rate for the ECG signal to obtain an accurate time-domain HRV parameter without interpolation at 1 kHz. They noticed that by decreasing the ECG sampling interval, the poorest accuracy and precision was for the pNN50 measure. Mahdiani et al generated down-sampled ECGs and investigated the error rate (Mahdiani et al 2015). They argued that a 50 Hz sampling rate could be used for measuring the ECG signal without compromising the accuracy of the calculated time domain in HRV parameters. However, this was not statistically validated, and further research is still required.
Even with the popularity of PPG-based personal healthcare applications, research on finding the proper sampling rate for PPG has not been conducted. The purpose of this study was to investigate changes in PRV results according to a given sampling frequency and to find the minimum sampling frequency, which could maintain the coherence for PRV analysis. For this, we analyzed PRV not only with a high rate-sampled PPG but also with the down-sampled PPG generated systematically by decimation. We then validated the significance of coherence of PRV analysis result between the high-sampled PPG and the down-sampled one.

Data summary
We obtained PPG and ECG simultaneously from 30 volunteers and used the data from 28 subjects, except for duplicated data by operational mistake or distorted data from severe motion artifacts. Finally, the dataset from 28 subjects (20 males, 8 females, age 22.0 ± 2.9 years, height 170.2 ± 7.3 cm, and weight 70.5 ± 15.4 kg) was used for analysis. Every data was recorded at left index finger for 20 min with a 10 kHz sampling frequency in supine position. Subjects were asked to be at rest but not to sleep or to make motion noise during the signal acquisition. Every participant reported not having any cardiovascular symptoms or diseases, impairment or other long-term health problem. Also, all participants were instructed to avoid caffeine, alcohol and smoking a day before and during the experiment as they could affect the autonomic nervous activity.
We generated down-sampled PPGs by decimation to investigate the sampling frequency effect on PRV analysis. Before decimation, we interpolated the 10 kHz-sampled PPG with 15 kHz sampling rate because 10 kHz cannot be divided in frequencies such as 15 Hz. Finally, we generated 11 down-sampled PPGs: 5000 Hz, 2500 Hz, 1000 Hz, 500 Hz, 250 Hz, 100 Hz, 50 Hz, 25 Hz, 15 Hz, 10 Hz and 5 Hz. The MP150, ECG100C, OXY100E and Acknowledge 4.0 (Biopac Systems Inc., Santa Barbara, CA, USA) were used for 16-bit analog and digital conversion, ECG amplification, PPG amplification and data recording, respectively. These modules were connected to record ECG and PPG simultaneously. The TSD124A, transmission-type finger clip optical transducer which includes 910 nm light emitting diode, was also used for PPG recording, and this configuration is possible to record PPG correctly within 18-321 bpm pulse rate range from manufacturer's datasheet.

Data analysis
We analyzed PRV on time-and frequency-domain with 20 min length data and the typical variables-AVNN (average NN interval), SDNN (standard deviation of NN interval), SDSD (standard deviation of successive difference of NN interval), RMSSD (root mean square of successive difference of NN interval), NN50 (number of pairs of successive NNs that differ by more than 50 ms), pNN50 (the proportion of NN50 divided by total number of NNs), TP (total power), VLF (very low frequency power), LF (low frequency power), HF (high frequency power), LF/HF (ratio of LF and HF), nLF (normalized LF), nHF (normalized HF)-were derived (Cardiology 1996). Long-term variables like SDANN (standard deviation of all the 5 min NN interval) and SDNN index (standard deviation of all the 5 min NN interval) were excluded in our experiment because SDANN and SDNN index are used for the 24 h HRV analysis like holter records in most cases (Taskforce 1996).
Unlike the time-domain analysis, frequency-domain analysis of PRV needed preprocessing for IBI. The preprocessing for this research were as follows: deriving the RR tachogram, cubic spline interpolation of tachogram with 1 kHz sampling rate, resampling with 4 Hz sampling rate, linear trend removal, and windowing with Hanning window. The power spectrum density (PSD) was then estimated with magnitude spectrum of fast Fourier transform (FFT), and VLF, LF and HF power were calculated by integrating PSD on the interval (0.0033, 0.04), (0.04, 0.15) and (0.15, 0.4) with respect to frequency (Hz), respectively.
Every PRV analysis was conducted with upper peak of PPG waveform that indicates maximum systolic phase. We used adaptive threshold-based PPG peaks detection algorithm as an initial screening (Shin et al 2009) after finite impulse response (FIR) band pass filtering of 0.5-10 Hz passband. Then every miss-detected peak location was manually corrected by mutual agreement of three experienced researchers using proprietary developed peak correction GUI (graphical user interface) in figure 1. Developed GUI includes waveform panel and tachogram panel to represent PPG waveform and beat-to-beat interval, respectively. Context menu of GUI provides four functions; remove current peak location, find new peak location based on maximum value, find new peak location based on minimum value and export peaks. In peak location correction, first we searched every miss-detected peak location from detected peak location set one by one, and then delete or re-find correct peak locations using graphical interface.
Matlab 2015b (MathWorks Inc., Natick, MA, USA), which was used for the whole procedure of signal processing and data analysis.

Statistical analysis.
Derived PRV were compared with the result of HRV derived from the 10 kHz-sampled ECG, and statistically validated by one-way ANOVA. Post hoc tests were then conducted if the result by sampling frequency was not specified a priori (P < 0.05). As a post hoc test, when the variances could be assumed to be equal, we used the Bonferroni test. Otherwise, we used Tamhane's T2 test. Levene's test (Levene 1960), which is an inferential statistic used to assess the equality of variances for a variable calculated for two or more groups, was used to assess the equality of variances for the derived PRV variables. Bonferroni test, which is appropriate when the variances are equal, uses t-tests to perform pairwise comparisons between group means, but controls overall error rate by setting the error rate for each test to the experimentwise error rate divided by the total number of tests. Hence, the observed significance level is adjusted for the fact that multiple comparisons are being made. Tamhane's T2, which is appropriate when the variances are unequal, is a conservative pairwise comparisons test based on a t-test. All statistical analyses were conducted using IBM SPSS Statistics software, Version 23.0 (IBM Corp., Armonk, NY, USA). The Bland-Altman analysis was also adapted to compare between the derived variability variables of HRV and PRV, according to the sampling frequency. In Bland-Altman analysis, we calculated bias and 95% limits of agreements (LoA) to evaluate the agreement of the down-sampled PRV compared with the 10 kHz HRV.

PPG characteristics
An example of the transition on waveform is shown in figure 2(a) along with power spectral density of PPG in figure 2(b). From the top, the figures respectively show a waveform or a PSD of 10 kHz-sampled PPG, a transition image based on various down-sampled PPG, and a waveform or PSD of 5 Hz-sampled PPG. Figure 2(a) is a truncated PPG waveform (2 s duration). In this figure, we can find that the PPG waveform becomes more simplified as the sampling frequency decreases. The changing of waveform shape started around f s = 25 Hz and the displacement of peak positions was also observed at lower frequencies. In figure 2(b), PSD of PPG has a periodic characteristic, for which power is decreased over frequency. This is a well-known phenomenon that the frequency component of PPG has the standing wave around the heartbeat frequency and the harmonic characteristic. At top of figure 2(a), we can find the fundamental frequency related to the heartbeat of around 1 Hz and four other harmonics. Middle of figure 2(b) shows that the specific frequency range was omitted by the decreasing sampling frequency. For example, the spectrum disappeared at >5 Hz in P 10 and at >2.5 Hz in P 5 . This phenomenon could be explained with Shannon's sampling theorem. Because the sampling frequency should be higher than or equal to twice the signal bandwidth to prevent aliasing, 10 and 5 Hz sampling could not contain the frequency components of over 5 and 2.5 Hz, respectively. As the PSD of PPG was decreased over the frequency change, losing the higher frequency component had less of an effect on the signal obtained. However, in this example, we could confirm that relatively large frequency components (1st to 5th harmonic) started to disappear at f s ⩽ 15 Hz. (a) Upper panel shows the beat-to-beat interval to support for finding detection error; lower panel shows the waveform and detected peak locations; context menu includes functions for editing such as 'delete', 'find max', 'find min', and 'save peaks'. (b) Context menu function will be applied after selecting the interest region by dragging.

Reliability of PRV analysis result over sampling rate
Tables 1 and 2 are the statistical result of HRV and PRV analysis on time-and frequencydomain, respectively. Each row is for each variable of PRV analysis, and the columns are the signal types. The notation 'Signal x ' denotes the signal for the sampling frequency x. For example, ECG 10000 denotes the 10 kHz-sampled ECG. All values in the tables are represented as mean ± SD. In tables 1 and 2, a and b denotes the type of post hoc test used, namely Bonferroni and Tamhane's T2 tests, respectively. To indicate statistical significance for the P-values, an asterisk (*) was used: * P < 0.05, ** P < 0.01 and *** P < 0.001.
From a one-way ANOVA test, significant difference (P < 0.05) was found in SDNN, SDSD, RMSSD, NN50, pNN50, TP, HF, LF/HF, nLF and nHF, but not in AVNN, VLF and LF. Thus, we then conducted post hoc analysis between down-sampled PPGs. Before the post hoc test, we performed the Levene's test to prove equality of variance and found that every variable had the equal variance except on the NN50, pNN50 and LF/HF. Therefore, we performed   Tamhane's T2 test for NN50, pNN50 and LF/HF, and Bonferroni test for other variables as post hoc tests. The post hoc test analysis indicated a significance difference (P < 0.05) between the frequency groups; NN50 and pNN50 were significantly different at f s ⩽ 20 Hz. In addition, SDSD and RMSSD were significantly different at f s ⩽ 15 Hz; also, LF/HF and nHF were significantly different at f s ⩽ 10 Hz. SDNN, TP and nLF were significantly different at f s = 5 Hz, and no significant differences were found for any of the variables at f s ⩾ 25 Hz.
The results could be interpreted that there were no differences in the frequency-domain analysis between reference HRV and PRV when PPG's sampling frequency was over 25 Hz. However, at f s ⩽ 20 Hz, significant differences were found for some of the variables. Overall, the frequency-domain variables showed robust characteristics at a low sampling rate compared with the time-domain variables. Also, the low frequency-related variables such as AVNN, SDNN, VLF and LF were less effected by down-sampling. Figures 3 and 4 graphically represent the mean and standard error of variables on timedomain analysis (figure 3) and frequency-domain analysis (figure 4). Because the results did not show significant differences for f s > 25 Hz, we show the results of HRV analysis for f s ⩾ 100 Hz. As in figure 3, AVNN was not changed significantly despite the sampling frequency decreases. However, SDNN, SDSD and RMSSD gradually increased by downsampling. The NN50 and pNN50 had no specific trend over the sampling frequency; however, values were dramatically increased with statistical significance at f s ⩾ 20.
In figure 4, we can find that variables of frequency analysis changed with a specific trend except for VLF. Values for TP, LF, HF and nHF increased, but LF/HF and nLF decreased by the decreasing sampling frequency. Remarkable differences were found at f s = 5 Hz. (n.u.: null unit) Figures 5 and 6 respectively show the results of the Bland-Altman analysis between the 10 kHz ECG based on HRV and the down-sampled PPG based on PRV. From the time-domain analysis as in figure 5, greater bias was observed by decreasing the sampling frequency. Moreover, the width of LoA was increased with lower sampling frequencies. Tendency of bias and LoA showed different characteristics according to the variable analyzed. Bias was slightly increased in AVNN, abruptly in SDNN, SDSD, and RMSSD, and randomly in NN50 and pNN50. LoA width was abruptly increased at f s > 10 Hz in SDNN, SDSD, and RMSSD, at f s > 20 Hz in pNN50 and NN50, and at f s < 15 Hz in AVNN. Figure 6 shows the Bland-Altman analysis of the result of the frequency-domain analysis. Similar to the result from time-domain analysis, bias and LoA width were increased by the decreasing sampling frequency. In terms of each variable, bias and LoA were abruptly increased at f s < 10 Hz in TP, VLF, LF and HF. However, it was also found that bias and LoA were decreased abruptly at f s < 25 Hz in LF/HF, nLF and nHF.

Reliability change in PRV by decreasing the sampling rate
The results could be discussed with respect to a number of viewpoints related to sampling frequency. First, we found the number of significant differences with respect to various variables time-and frequency-domain variables of the PRV were dramatically higher at f s < 10 Hz. It could be postulated that this phenomenon originates from the spectral characteristic of PPG. As seen above, the power spectrum of PPG consists of the fundamental frequency around the heart rate, which has most of the energy and its harmonics. Considering a normal range for a heart rate of 70 bpm (1.2 Hz), a 10 Hz sampling frequency could cause a distortion of the 5th (6.0 Hz) or higher order harmonic components by Shannon's sampling theorem (Oppenheim and Schafer 2010). Besides, at 5 Hz sampling rate, aliasing could have occurred at a frequency of over 2.5 Hz and this could cause a severe distortion as normally this frequency range is around the 2nd harmonic.
From the PRV analysis by decimation, the changes followed a certain rule: the lower the sampling frequency was, the more irregular the waveform became. In other words, the variation in the peak location would increase by decreasing the sampling frequency and it could cause the variations seen from analysis of PRV. Increased variation could also raise the total frequency power. This was shown in figure 4(a) by the increases in total power at lower f s values, particularly for f s = 5 Hz.
To explain the changes for other variables, we should remember the implicit prerequisite for our simulation conditions that at the very least, the pulsatile shape of PPG waveform needs to be preserved. From this prerequisite, we set 5 Hz as a minimum frequency of simulation and in this condition, every peak location could be varied most within a pulse duration. Any slight variation of peak location could increase total frequency power; particularly, HF power would be more affected as HF is related to the short-term changes on the time domain.
This trend is clearly represented in the time-domain analysis. Significant differences were apparent for SDSD and RMSSD at f s ⩽ 15 Hz and in NN50 and pNN50 at f s ⩽ 20 Hz and these variables are known to reflect high frequency activity derived from parasympathetic activation (Taskforce 1996). Otherwise, variables that are known to reflect low frequency activity-AVNN and SDSD-showed little or no significant changes. This result is apparent with Bland-Altman analysis. In figure 5, bias and LoA varied only within 4 ms in AVNN, even at f s = 5 Hz. For SDNN, bias and LoA were around 15 ms at f s = 10 Hz and around 40 ms at f s = 5 Hz. This was a very rapid variation compared with AVNN. However, it is a relatively tiny variation compared with changes in SDSD or RMSSD. Bias and LoA variation of SDSD or RMSSD were around 50 ms at f s = 10 Hz and around 100 ms at f s = 5 Hz. Furthermore, for NN50 and in pNN50, bias and LoA showed abrupt variations at f s ⩽ 20 Hz.
The results of the frequency-domain analysis could be interpreted in the same way as above. Increasing the irregularity of peak position increases the total power. In this respect, high frequency variation greatly contributes to the total power increases. Figure 4(a) shows that the total power was increased by decreasing the sampling frequency. The VLF, LF, HF and nHF also show an increasing trend at low sampling frequencies. However, we also found that the extent of the increases depended on the variables analyzed: VLF and LF showed a slight increase without significance, but HF and nHF showed abrupt and significant changes. As HF increased to a much higher extent than LF at low sampling frequencies, the variables LF/HF and nLF also decreased by decreasing the sampling frequency. Similar to the timedomain results, Bland-Altman analysis (figure 6) provided another viewpoint of significant change in variables according to the sampling frequency change. In case of VLF, bias and LoA showed little variation (<50 ms 2 of bias) in every sampling frequency tested. However, for TP and HF, bias and LoA were dramatically (~100 ms 2 of bias) increased at f s = 5 Hz. Bias and LoA were also increased for LF at f s = 5 Hz, but the increment was relative small (~150 ms 2 of bias). The variation of LF/HF, nLF and nHF are gradually increased and in every case larger variation was observed around 5 Hz or 10 Hz sampling frequency.

Limitations
The above results demonstrate that decreasing the sampling frequency could degrade reliability of the PRV analysis. However, there are additional conditions such as changes in the heartbeat interval could change the pulse wave characteristics. For example, exercise or certain body positioning like standing could increase the heart rate, and influence the change in pulse rate as well as the shape of the pulse wave. Moreover, various cardiovascular conditions need to be considered for generalizing the relationship between the sampling frequency and the reliability of the PRV analysis; for example, cardiovascular disease (CVD) such as arrhythmia could also change the pulse wave characteristics. The limitation of our research is that we focused only on normal conditions for our test volunteers such as having a normal heart rate, being in the supine position and resting (without sleeping). Therefore, additional analysis and validation in more complex cardiovascular conditions including subjects with CVD or  subjects in various additional behaviors for a more generalized result for setting the lower limit on the frequency sampling rate.

Conclusions
Since the PPG and ECG have different wave shapes, it is natural that they are affected differently by the sampling frequency changes. Considering that the practical uses of PPG have increased in the modern society, the effect of the sampling rate for PPG-based applications deserves attention. In this research, we tried to find the proper sampling frequency to analyze PRV, and the results suggest that the minimal sampling rate from PRV analysis depends on the PRV variables analyzed; namely, f s > 20 Hz for NN50 and pNN50, f s > 15 Hz for NN50 and pNN50, f s > 10 Hz for LF/HF and nHF, and f s > 5 Hz for SDNN, TP, HF and nLF are required. One thing we should consider is that our research was based on the perfect detection of the PPG peak position. In actual situations where the performance of the peak detector could be degraded, the reliability of the analysis result could vary. As such, to apply the result of this research to actual applications, the performance of the peak detector at low frequency should be considered with an adequate sampling frequency. In conclusion, our pilot study suggests that analysis of variability in the time and frequency domain from pulse rate obtained through PPG may be potentially as reliable as that derived from the analysis of the electrocardiogram, provided that appropriate sampling frequency is used. Further studies in different conditions and situations, for the elderly subject or for those with cardiac diseases, are needed to prove the reliability of the methodology.