Instantaneous frequency from Hilbert-Huang transformation of digital volume pulse as indicator of diabetes and arterial stiffness in upper-middle-aged subjects

To investigate the value of decomposed short-time digital volume pulse (DVP) signals in discerning systemic vascular anomaly in diabetic patients, demographic and anthropometric parameters, serum lipid profile, fasting blood glucose and glycated hemoglobin (HbA1c) levels were obtained from 29 healthy adults (Group 1) and 29 age-matched type 2 diabetes mellitus patients (Group 2). Six-second DVP signals from right index finger acquired through photoplethysmography were decomposed using ensemble empirical mode decomposition. Using one intrinsic mode function (IMF5), stiffness index (SI) and instantaneous energy of maximal energy (fEmax) were obtained. Other indicators of arterial stiffness, including electrocardiogram-pulse wave velocity of foot (ECG-PWVfoot), crest time (CT) and crest time ratio (CTR), were obtained from the testing subjects for comparison. The mean body weight, body mass index, waist circumference, HbA1c and fasting blood sugar levels were higher in Group 2 than those in Group 1, whereas values of systolic and diastolic blood pressure were lower in Group 2 than those in Group 1. SI and fEmax were significantly higher in Group 2 than those in Group 1. Moreover, fEmax was positively associated with HbA1c concentration, CT and SI in Group 2 (p < 0.05) but not in Group 1. When all subjects were considered, fEmax was highly significantly associated with HbA1c and fasting blood sugar levels, and SI (all p < 0.001). After Hilbert-Huang transformation, short-time DVP signals could give significant information on arterial stiffness and vascular anomaly in diabetic patients.


could give significant information on arterial stiffness and vascular anomaly in diabetic patients.
Not only is cardiovascular disease a major killer in the developed world, but it is also a non-communicable disease that poses an ever-increasing health threat in the developing countries 1 . Besides, the gruesome disabilities for the survivors result in staggering social and economic burden worldwide 2 . Accordingly, prevention and early detection of the disease in high-risk patients such as those with diabetes have become important issues in preventive medicine and public health 3 .
Cost-effective and non-invasive approaches to early detection of cardiovascular disease have been widely investigated. Parameters including pulse wave velocity (PWV) 4 , cardio-ankle vascular index (CAVI) 5 , and ankle-brachial blood pressure index (ABI) 6 have been validated as non-invasive assessment tools for vascular health. Despite their non-invasiveness, unfavorable factors including the cost of the equipment, the need for technical assistance during measurement, the relatively long time for data acquisition as well as their hospital-based settings hinder their popularity as screening tools. With the advance of artificial intelligence and computational technologies, it is now possible to analyze physiological signals using portable electronic equipment with built-in computational programs to produce affordable devices for daily use 7 .
Arterial waveform analysis previously focused on direct observation of changes in waveform signals with the time domain as in the computation of stiffness index (SI) and crest time (CT) 8 . Digital volume pulse (DVP) acquired through photoplethysmography (PPG) has gained much popularity in assessing vascular health because of its non-invasiveness, convenience, and low cost 9 . After detrending and eliminating background noises using ensemble empirical mode decomposition (EEMD) (i.e., a component of HHT), DVP signals from the finger have been previously shown to successfully differentiate between normal and diabetic subjects using the multiscale cross-approximate entropy approach 10 .
Since conventional Fourier analysis only displays data as sine and cosine functions, it cannot help in the computation of stiffness index (SI) that requires discernible systolic and diastolic digital volume pulse (DVP) waveforms which are often difficult to obtain in subjects with systemic diseases such as diabetes 11 . Ensemble empirical mode decomposition (EEMD), which is the first component of Hilbert-Huang transformation (HHT), separates acquired physiological signals into a set of distinct physiological information known as intrinsic mode functions (IMFs) 12 . It has been shown that IMF5 best represents DVP waveforms through which SI can be calculated 13 . It is our hypothesis that physiological information hidden in the inseparable signal of an IMF as reflected in the energy-frequency width spectrum can be revealed through Hilbert-Huang spectrum analysis, which is the second component of HHT. Therefore, the present study aimed at assessing the significance of instantaneous frequency of maximal energy (f Emax ) derived from marginal Hilbert spectrum and instantaneous frequency in vascular health through comparing the index in diabetes subjects with that in healthy volunteers. Established indices of arterial stiffness including SI, electrocardiogram-pulse wave velocity of foot (ECG-PWV foot ) 14 , crest time (CT) and crest time ratio (CTR) 11 were also acquired from the testing subjects for comparison.

Results
Baseline characteristics of study subjects. Diabetic group (Group 2) had significantly higher mean body weight, body mass index (BMI), and waist circumference (all p < 0.05) as well as concentrations of glycated hemoglobin (HbA1c) and fasting blood sugar (both p < 0.001) than those in the non-diabetic group (Group 1) ( Table 1). On the other hand, the systolic (SBP) and diastolic blood pressure (DBP) were lower in Group 2 than those in Group 1 (both p < 0.05).

Comparison of computational parameters between diabetic and non-diabetic patients.
There was no significant difference in ECG-PWV foot , CT and CTR between subjects without diabetes (Group 1) and those with the disease (Group 2) (p = 0.082, 0.287, and 0.896, respectively) (Table 1). However, stiffness index (SI) and instantaneous energy of maximal energy (f Emax ) were significantly higher in the diabetic patients (Group 2) than those in the non-diabetic subjects (Group 1) (p < 0.05 and <0.001, respectively) ( Table 1, Fig. 1). Correlations between f Emax and risk factors of atherosclerosis. Despite the lack of significant correlation between f Emax and conventional risk factors of atherosclerosis in the non-diabetic subjects (Group 1), f Emax was found to be positively associated with plasma concentration of glycated hemoglobin (p = 0.011) and stiffness index (SI) (p = 0.029) but negatively correlated with crest time (CT) (p = 0.015) in the diabetic patients (Group 2) ( Table 2). When all testing subjects were taken into account, f Emax was shown to be highly significantly associated with glycated hemoglobin and fasting blood sugar levels as well as stiffness index (all p < 0.001) but negatively correlated with CT (p = 0.017).

Discussion
Not only is the present study the first to apply EEMD in acquiring only six-second detrended DVP signals for assessing arterial stiffness, but it also unveiled the clinical significance of Hilbert-Huang spectrum-derived instantaneous energy of maximal energy (f Emax ) in differentiating non-diabetic middle-aged and elderly subjects from those with diabetes by showing highly significant positive associations of f Emax with key diabetes parameters and stiffness index. Despite currently available non-invasive equipment for the assessment of vascular health, the cost of the devices and the requirement for designated personnel for operation restrict their being used as routine screening tools. For instance, the use of the air pressure sensing system (APSS) mandates the use of an inflatable pressure cuff and full assistance for measurement. By contrast, DVP acquired through photoplethysmography (PPG) is a parameter that can be acquired using simple portable devices. We have previously shown that discrepancy in DVP waveform amplitude between two upper limbs can differentiate young from aged subjects and patients with good from those with poor diabetic control as well as identify arteriosclerosis risks 15 . Using the multiscale entropy approach, DVP waveform amplitudes have also been demonstrated to differentiate healthy from diabetic subjects 16,17 .
Nevertheless, the major obstacle to the use of DVP is the noises that notably interfere with data interpretation. Ensemble empirical mode decomposition (EEMD), one of the two components of the Hilbert-Huang transformation, is an adaptive time-frequency data analysis method found to be useful for extracting signals from data produced in a variety of non-stationary and nonlinear processes [18][19][20] . By applying EEMD, the current study demonstrated that DVP signals can be detrended and decomposed into IMFs for arterial stiffness computation (i.e., IMF5). Moreover, the current study demonstrated that the other component of the Hilbert-Huang transformation, Hilbert-Huang spectrum, can be utilized to compute the instantaneous frequency of maximal energy (f Emax ) that exhibited significant positive association with the parameters of acute (i.e., fasting blood sugar) and chronic (i.e., HbA1c) blood sugar control as well as the conventional stiffness indices (i.e., CT and SI).
On the other hand, the question arises regarding the cause of increase in f Emax in the diabetic subjects. One plausible explanation would be the production of diabetes-associated advanced glycation end-products (AGEs) that causes collagen cross-linking in the arterial medial layer, thereby contributing to arterial stiffness 21 . The results of the present study further support the link between diabetes and arterial stiffness. Additionally, diabetes has been well-documented as a risk factor for the formation of atherosclerotic plaques that can promote turbulence in the laminar blood flow 22 . Therefore, it is rational to propose that the turbulence thus generated may partly account for the increase in f Emax in the diabetic patients in the current study. The lack of significant difference in conventional parameters for assessing arterial stiffness (i.e., ECG-PWV foot , CT, CTR) between non-diabetic and diabetic subjects in the present study may further highlight the unique characteristic of f Emax in reflecting diabetes-associated vascular dysfunction. This study has its limitations. Firstly, the sample size of the present study is relatively small. On the other hand, the highly significant associations of f Emax with the parameters of sugar control and arterial stiffness highlight the significance of our findings. Secondly, other clinical conditions associated with atherosclerosis such as hypertension, smoking, and hyperlipidemia have not been separately examined to compare their impacts on f Emax . Finally, there were several inherent limitations interfering with proper acquisition of DVP signal with photoplethysmography sensors. To obtain reliable data, potential interferences including skin pigmentation, tissue characteristics, blood flow in the measured area, involuntary vibrations of the subjects being examined as well as low temperature-induced peripheral vessel constriction should be minimized. Therefore, all measurements were performed in a noise-free, humidity-controlled room with the temperature maintained at 26 ± 1 °C.
In conclusion, the results of the present study showed that Hilbert-Huang transformation not only enabled the assessment of arterial stiffness from merely 6 seconds of decomposed digital volume pulse signals, but it can also produce an instantaneous frequency of maximal energy that correlated positively with the parameters of blood sugar control and conventional arterial stiffness indices. The findings shed light on the possibility of adopting this simple parameter as an index of arterial stiffness and diabetes control. The short time required for measurement and the simplicity of equipment underscore its potential widespread use in portable devices.

Methods
Study population. From July 2009 to October 2010, 63 middle-to-old aged men and women with and without diabetes mellitus type 2 were recruited. All diabetic patients were enrolled from the outpatient clinic for diabetes care of the Hualien Hospital. The diagnosis of diabetes was based on either a fasting glucose concentration higher than 126 mg/ dL or a glycosylated hemoglobin (HbA1c) level >6.5% 23 . All patients underwent regular treatment as well as follow-up in the clinic for over two years. The criteria for enrollment of non-diabetic subjects included a negative history for diabetes mellitus, an HbA1c level less than 6.5% and a fasting glucose level lower than 126 mg/dL. Subjects with history of atherosclerosis-associated complications, such as stroke, angina, myocardial infarction and peripheral vascular diseases within three months of the present study were excluded. Of the 63 subjects recruited, five were excluded because of either inadequate length of follow-up for diabetic patients or history of known cardiovascular diseases. As a result, totally 58 subjects were recruited who were divided into two groups, including 29 non-diabetic subjects (15 females,14 males) recruited from a routine annual physical check-up program at the same hospital (Group 1) and 29 patients with diabetes mellitus type 2 (15 females,14 males) (Group 2). The protocol and procedures of this study were in accordance with the principles of the Declaration of Helsinki. The study was approved by the Institutional Review Board of Hualien Hospital. Informed consents were obtained from all study participants.  Acquisition of anthropometric, serum biochemical, and hemodynamic parameters. All measurements and procedures were performed in the morning (i.e., 8:30-10:30 a.m.). Demographic (i.e., age, gender) and anthropometric (i.e., body weight, body height, waist circumference, body mass index) parameters, serum lipid profile (i.e., high-density lipoprotein cholesterol, low-density lipoprotein cholesterol, total cholesterol, and triglyceride), and parameters of glucose control (i.e., fasting blood glucose and glycated hemoglobin [HbA1c] levels) were obtained at the hospital. Resting blood pressure was measured once over the left arm of the supine participants using an automated oscillometric device (Microlife BP3AG1, Taiwan) with a cuff of appropriate size. Cholesterol and triglycerides concentrations were determined from blood samples obtained after overnight fasting for 12 hours. The participants were asked to refrain from caffeine-containing beverages and theophylline-containing medications for at least 12 hours before each hospital visit. Additionally, to minimize potential errors in infrared sensor readings arising from involuntary vibrations of the examinees and a decreased environmental temperature that may cause constriction of peripheral vessels, all participants received blood sampling before data acquisition and were allowed to rest in a supine position for 10 minutes in a quiet room with temperature maintained at 26 ± 1 °C. The PPG sensor was simultaneously applied to left index fingertip of each participant for the acquisition of data for 30 minutes 16,17 . From all testing subjects, 6 seconds of DVP signals were acquired after the start of recording for 5 minutes.

Protocol of measurement of DVP and other computational parameters.
The protocol for acquisition of DVP signals with photoplethysmography was described in our previous study 10 . Briefly, a self-developed six-channel electrocardiography (ECG)-PWV-based equipment was used to obtain 3000 successive DVP signals within 6 seconds from the right index finger at a frequency of 500 Hz. Infrared sensors were put on the points of reference simultaneously to acquire data. After being processed through an analog-to-digital converter USB-6009 DAQ, National Instruments, Austin, TX) with a sampling frequency of 500 Hz, the digitized signals were stored in a computer 24 . The DVP signals thus obtained were subject to two-staged Hilbert-Huang transformation (HHT) that included ensemble empirical mode decomposition (EEMD) and Hilbert-Huang spectrum. For comparison, ECG-PWV foot , crest time (CT) and crest time ratio (CTR) from each subject were acquired. The computation of ECG-PWV foot was conducted as previously described 14 . The time difference (ΔT) between the R peak of an electrocardiogram and the foot point of a DVP waveform during the same cardiac cycle was obtained. The distance from the sternal notch to the foot was the sum of the shortest distance from the sternal notch to medial patella, from medial patella to medial malleolus, and from medial malleolus to the tip of the second toe of each foot. ECG-PWV foot was calculated by dividing the distance with ΔT and taking the average from both sides 14 . Crest time (CT) is defined as time from foot point to peak of a pulse wave, while crest time ratio (CTR) is the ratio of CT to the time from one foot point to the next foot point (i.e., cycle time) 11 .

Hilbert-Huang transformation (HHT) of DVP signals.
In general, DVP signals comprise noise-free and noise components: in which y(t) is the acquired data, and s(t) and n(t) are the desired signal and white noise, respectively. As shown in (1), all data are combinations of signal and noise. To enhance the accuracy of measurements, the ensemble approach is adopted in which the data are collected by separate observations that may contain different noises. To generalize this ensemble concept, noise is introduced to the single data set y(t) as if different observations were actually being made as an analog to a physical experiment which could be repeated many times. Under such circumstances, the ith "artificial" observation will be i i In the case of a single observation, each multiple observation ensemble is simulated through the addition of not arbitrary but different realizations of white noise ω i (t) to that observation as shown in (2). EEMD decomposes the signal into different IMFs, each of which is a mono-component function 9 . Hilbert transformation is then applied to calculate the instantaneous frequencies of the original signal. After the identification of all local maxima and minima of the signal, the upper and lower envelopes are created through curve fitting. Mean values of the upper and lower envelopes of the signal, m 11 (t), are then obtained. Therefore, the difference between the signal y i (t) and its envelopes m 11 (t), which is denoted as h 11 (t), can be computed i 11 11 The approximation nature of the curve fitting method necessitates further procession of h 11 (t) (i.e., treating h 11 (t) as the signal itself with continual repetitions of the process) until satisfaction of the following two conditions 9 : (1) The number of extrema and the number of zero-crossings are either equal or differ by at most one, and (2) At all times, the mean value between the envelope defined by local maxima and that defined by local minima is zero.
Through iteration (for a total of k times), the difference between the signal and the mean envelope values, h 1k (t), is computed as where m 1k (t) is the mean envelope value after the kth iteration, and h 1(k−1) (t) is the difference between the signal and the mean envelope value at the (k − 1)th iteration. The function h 1k (t) is then denoted as the first IMF component: After removing IMF 1 (t) from the original signal y(t), the residue can be obtained The residue r 1 (t) can then be treated as the new signal, and the above iteration process is repeated to extract the rest of the IMFs from the signal y i (t) as: The signal decomposition process continues until r n (t) becomes a monotonic function, from which no further IMFs are extractable. The signal y i (t) is decomposed through the summation of (6) and (7) into a number of IMFs that are distinct components of the signal. Hence, the signal y(t) can be defined as: 1 2 n In this way, the data are decomposed into n-empirical IMF modes and a residue, r n (t), which can either mean a constant or trend in (8). The concept of EEMD is to add white noise, which populates the whole time-frequency space uniformly with the components of different scales separated by a filter bank. The EEMD process is further explained as follows 9 : Step 1. Addition of a white noise series to the targeted data; Step 2. Decomposition of the data with added white noise into IMFs; Step 3. Repetition of step 1 and step 2, but with different white noise series each time; Step 4. Acquisition of the (ensemble) means of corresponding IMFs of the decompositions as the final result. When the number in the ensemble approaches infinity, the result of EEMD is obtained in which IMF i (t) + αr k (t) is the kth realization of the ith IMF in the signal with added noise, α is the standard deviation of the added noise, and r k (t) is the residual after extraction of the first k IMF components. The number of trials in the ensemble (N) needs to be large.
In the present study, six-second DVP signals were obtained from each testing subject with 3000 successive samplings at a rate of 500 Hz. Hence, α was set to be 0.2, and N was equal to 3000 for easy computing 20,25 . Using the Matlab package, the original 3000 samplings within the six-second DVP signals from a testing subject (i.e., y(t) in Eq. (1)) were subject to EEMD to obtain Eqs (2)(3)(4)(5)(6)(7)(8). The process generated 8 IMFs and residual r 8 (t) as illustrated in Fig. 2 (Left panel).
According to the definition of stiffness index (SI) 26 , it can be computed from IMF5(t) in Fig. 2 as follows: DVP Hilbert-Huang transformation (HHT) is signal analysis in the time-frequency domain by combining EEMD with Hilbert transformation 16,19,20 . Unlike Fourier analysis that produces a series of sine and cosine functions of fixed amplitudes to represent each frequency constituent in the signal, the Hilbert-Huang spectrum approach is based on the computation of instantaneous frequency from Hilbert transformation of the signal. Generally, the definition of Hilbert transform for the signal IMFn(t) is In theory, any complex signal z(t) can be considered to be the sum of its real part IMFn(t) and imaginary part Im(t) so that To be expressed in a polar coordinate system, Eq. (11) can be rewritten as 1 denote the instantaneous amplitude and phase of the analytic complex signal, respectively. The instantaneous frequency ω(t) of the signal can be derived from the instantaneous phase θ(t) as where Re{z(t)} represents the real part of the complex signal z(t). To ensure that the instantaneous frequency obtained from Eq. (15) is physically meaningful, the instantaneous phase θ(t) must be a mono-component function (i.e., a single-valued function over time).
As shown in (13) and (14), the instantaneous phase θ(t) is derived from IMFn(t) and its Hilbert transform. Therefore, IMFn(t) must also be a monocomponent function. Eq. (13) allows the amplitude and the instantaneous frequency to be presented in a 3-D plot, in which the amplitude is represented by the height in the time-frequency plane. This time-frequency distribution thus generated is known as "Hilbert-Huang spectrum" H(ω, t): On the other hand, the marginal spectrum, h(ω), can be defined as (Fig. 2, right panel). including coefficient of variation, t test, and correlation. Data were expressed as mean ± standard deviation (SD). The continuous variables between the two groups were compared using two-tailed t-test. Correlations between instantaneous frequency (f Emax ) and risk factors was analyzed with Pearson correlation test. A p value < 0.05 was considered statistically significant.

Data Availability
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.