Principle of PFDM
The core process of PFDM is frequency demodulation based on the method of complex demodulation (CDM) [15–17]. CDM is a non-linear time-domain method for time series analysis, which provides amplitude and frequency of non-stationary/unstable oscillatory signal as a continuous function of time. Principle and computer algorithm of CDM have been published previously [16]. Briefly, CDM extracts the time dependent functions of instantaneous frequency through the following four steps: (1) the spectral region of interest (the frequency range of target oscillation) is shifted to zero frequency by forming a product, throughout the record, of the original signal and a complex sinusoid at a reference frequency (Fr, the center frequency of the spectral region of interest), (2) the resultant complex signal is low-pass filtered so that only frequency components around zero remain, (3) the real and imaginary parts of the low-pass filtered signal are converted to a polar form, yielding the instantaneous phase, as a function of time, of the component identified at or near the Fr, and (4) the time series of frequency of the target oscillation is calculated through adding the first-order derivative of the phase to the Fr, since the slope of the phase versus time curve indicates deviation of instantaneous operative frequency from the Fr.
In PFDM, the CDM is customized for analyzing pulse wave signal. Pulse frequency (instantaneous pulse rate) could change widely from less than 0.5 Hz (30 bpm) to more than 3 Hz (180 bpm). Because CDM extracts frequency of oscillatory components within an assigned spectral region (range of CDM filter), an appropriate frequency range needs to be selected so that it covers the possible range of pulse frequency. This is, however, a trade-off with the requirement that the frequency range of CDM filter needs to be narrow enough to avoid influence of harmonics and subharmonics of the fundamental oscillation (pulse wave).
When frequency of pulse wave decreases to a frequency as low as f Hz, the upper limit of the frequency range should be less than 2f Hz (frequency of the 2nd harmonic of the pulse wave). Likewise, when frequency of pulse wave increases to a frequency as high as f Hz, the lower limit of the frequency range should be greater than f/2 Hz (frequency of the 2nd subharmonic of the pulse wave). If the range for CDM filter is expressed as Fr ± Fc, where Fr is the reference frequency (the center frequency of the spectral region of interest) and Fc is the corner frequency of the low-pass filter, then the necessary condition for avoiding the influence of harmonics and subharmonics is
Fc < Fr/3.
For example, when mean pulse rate is 60 bpm, the widest possible range for CDM filter is 60 ± 20 (40 to 80) bpm; if pulse frequency deviates from this range during analyzing period, CDM would not provide an adequate estimation of pulse rate.
To overcome this problem, we devised the algorithm of PFDM incorporating the following three features: (1) overlapping short-segmentation of data, (2) adaptive determination of the Fr, (3) iterative algorithm for optimizing Fr, and (4) stepwise convergence of Fc toward Fr/3 during the iterative process. Briefly, data are first divided into short segments so that pulse frequency within each segment can be expected to remain within the range that can be covered by a feasible CDM filter. In each segment, the mean pulse frequency of the previous segment is used as the initial Fr value (adaptation). The Fr is further optimized through iteration processes, i.e., mean pulse frequency calculated is iteratively used as the Fr in the next iteration until the mean pulse frequency calculated agrees with the Fr used. Finally, to further guard against the possibility that pulse frequency deviates from the range of CDM filter, an Fc of Fr/2 is used for the initial process and it is converged to Fr/3 in the following iteration processes. The effectiveness of each of these processes is demonstrated in the simulation studies.
Simulation data
Simulated pulse wave signal was generated as a sequence of a unit pulse wave that was sampled from actual pulse wave of a healthy subject. The pulse interval was modulated by various source interval functions (modulator functions), which were used as the reference signals for evaluating the performance of PFDM (Fig. 1). The width of each unit pulse wave was also modulated according to the changes in pulse interval as 18.97 × (PI)1/2 ms, where PI is the pulse interval.
Measurement of actual data
We studied 33 subjects (23 males and 10 females); the mean age (standard deviation [SD]) was 34 (7) yr and age range was 22–47 yr. They were recruited from the staffs in the work places of the authors. All subjects gave their written informed consent and the procedure was performed according to the Ethical Guidelines of Nagoya City University Graduate School of Medical Sciences. Continuous pulse wave was recorded from the dorsal side of the wrist with a wireless trans-dermal photoelectric pulse wave sensor system (Prototype C, DENSO, Japan, Fig. 2 and 3), by which the pulse wave data were digitized (14 bit). ECG was recorded simultaneously with a digital Holter recorder (RAC-2102, Nihon Koden, Tokyo, Japan). In order to obtain exact matching of time between the two recordings, several event markers were also recorded simultaneously. Both continuous pulse wave and ECG data were recorded during all-night sleep in a laboratory chamber, because the pulse wave device allowed stable measure of pulse wave signal only during rest at present. To our knowledge, no noninvasive ambulatory devices are currently available for measuring stable pulse wave signal during activities.
Data analysis
For both simulated and actual data analysis, pulse wave data were sampled at 20 Hz. The PFDM was performed with custom-made software written with FORTRAN 95 (Salford Software Ltd, Old Trafford, Manchester, UK). For the PFDM, the data segments had a length of 30 s with overlapping for 10 s at both ends. For the iteration for optimizing Fr, the tolerance for the difference between the mean pulse frequency and Fr was set at 0.001 bpm. Although the instantaneous pulse frequencies were calculated for every sampling point (at 20 Hz), they were averaged over every 500 ms and converted into pulse intervals in order to compare with the source interval functions (for the simulation data) or R-R intervals of ECG (for the actual data).
HRV, by definition, is the beat-to-beat variability of sinus rhythm [1, 18]. Thus, all R-R interval data involving ectopic beat(s) or heart block(s) were excluded from the analysis. The pulse wave data, however, have less information about arrhythmias. An ectopic beat, either supra-ventricular or ventricular, could result in consecutive short pulse intervals, a long pulse interval or even a normal interval depending on the timing and whether it generates a detectable pulse wave or not. Heart blocks may also result in long intervals. To avoid the effects of rhythm disturbances and noises on the analysis of PRV, we excluded all abnormal pulse intervals, which were defined as those deviating 12% or more from the local moving average over the preceding 20 s. The definition of abnormal pulse interval was determined tentatively considering expected range of physiological PRV for the subject population.
The ECG data were analyzed with a Holter scanner (DSC-3100, Nihon Koden, Tokyo, Japan), on which the results of automatic labeling of QRS complexes were reviewed and manually edited for all errors. The ECG were analyzed with a sampling frequency of 125 Hz and, thus, the R-R intervals were measured at a time resolution of 8 ms. The time series data of R-R interval were interpolated along the time axis with a horizontal-step function (R-R interval was considered as constant during each R-R interval) and resampled at 2 Hz.
Spectral analysis
Fast Fourier transformation with a Hanning window was performed for sequential 5-minute segments of both pulse interval and R-R interval data. A segment was excluded from the analysis, if the ratio of valid data points in the segment was <80%. In the segments analyzed, defected parts of data, if any, were interpolated by the horizontal-step function. After correcting for the loss of variance resulting from the window process, power spectral density was integrated over 0.04–0.15 Hz and 0.15–0.45 Hz for assessing the power of low frequency (LF) and high frequency (HF) components, respectively. The power of these components was converted into mean amplitude ([2 × power]1/2) to reduce the skewness from the normal distribution.
Statistical analysis
The agreement between PRV and HRV measures was evaluated from two aspects; (1) agreement when these measures are used for assessing intra-individual variations and (2) agreement when they are used for inter-individual comparisons.
To examine the former agreement, minute-to-minute pulse rate and heart rate and spectral components of PRV and HRV for 5-min segments were compared within each subject. The agreement between corresponding values were evaluated with (a) the difference-against-mean plot and the 'limits of agreement' of Bland and Altman method [19] and (b) intraclass correlation coefficients for 2-way mixed effects analysis of variance with defining segments as the random factor and methods as the fixed factor [20]. The statistical adequacy of the segments as random factor is partly supported by the fact that long-term heart rate fluctuation has the characteristics of random fractal noise [1, 21].
To examine the latter agreement for inter-individual comparisons, the minute-to-minute and 5-min segment values were averaged over the entire recording length for individual subjects. The agreement between the corresponding mean values were evaluated with (a) the Bland and Altman method same as above [19] and (b) intraclass correlation coefficients for 2-way mixed effects analysis of variance with defining subjects as the random factor and methods as the fixed factor [20].
Data were presented as mean (SD) or median (range). Type 1 error level was set at probability (p value) of < 0.05.
