J. Biomedical Science and Engineering, 2013, 6, 11431154 JBiSE http://dx.doi.org/10.4236/jbise.2013.612143 Published Online December 2013 (http://www.scirp.org/journal/jbise/) Application of Van der Pol oscillator screening for peripheral arterial disease in patients with diabetes mellitus Jianxing Wu1*, Chianming Li1,2, Weiling Chen1, Chiahung Lin3, Tainsong Chen1* 1Department of Biomedical Engineering, National Cheng Kung University, Tainan City, Taiwan 2Division of Infectious Diseases, Department of Medicine of Chi Mei Medical Center, Tainan City, Taiwan 3Department of Electrical Engineering, KaoYuan University, Kaohsiung City, Taiwan Email: *jian0218@gmail.com, 235813cmli@gmail.com, lynnchen.k@gmail.com, eechl53@gmail.com, *chents@mail.ncku.edu.tw Received 7 November 2013; revised 4 December 2013; accepted 12 December 2013 Copyright © 2013 Jianxing Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. In accor dance of the Creative Commons Attribution License all Copyrights © 2013 are reserved for SCIRP and the owner of the intellectual property Jianxing Wu et al. All Copyright © 2013 are guarded by law and by SCIRP as a guardian. ABSTRACT This paper proposes a Van der Pol (VDP) oscillator screening for peripheral arterial disease (PAD) in patients with diabetes mellitus. The longterm ele vated blood sugar levels produce a high risk of pe ripheral neuropathy and peripheral vascular disease, especially in the foot of a diabetic. Early detection is important, in order to prevent foot problems for dia betic patients with PAD. Photoplethysmography (PPG) is a noninvasive method for the detection of blood volume changes in peripheral arteries. Because of changes in the resistancecompliance, the rise time and transit time for the PPG signals increase and change in their shape are highly correlated with PAD severity. In this study, the Burg autoregressive (AR) method is used to determine the characteristic fre quencies of the right and leftside PPG signals. For bilateral frequency spectra, the VDP oscillator uses asynchronous signals. This produces damped sinu soidal responses and the oscillation overshoot (OS) is an approximating function only of the damped factor. This index is used to estimate the degree of PAD, in cluding normal the condition and diabetic patients with PAD. The results show that the proposed me thod is efficient and more accurate in the estimation of PAD. Keywords: Van der Pol (VDP) Oscillator; Peripheral Arterial Disease (PAD); Burg Autoregressive Method; Oscillation Overshoot 1. INTRODUCTION Type 2 diabetes is a chronic, metabolic disease that is prevalent in developed and developing countries, associ ated with mortality and morbidity due to increasing risk factor for cardiovascular diseases. Patients with diabetes mellitus incur substantantial medical expenditure to con trol their glycemic concentrations and complications. Peripheral arterial disease (PAD) is highly prevalent in patients with diabetes mellitus. However, it is usually underestimated by patients as well as physicians, espe cially in the foot of diabetics, mostly contributing to the slow progress of arterial atherosclerosis in lower limbs. Manifestations of this peripheral arterial obstruction can vary from a symptom to total occlusion or gangrene of toe or foot [1]. Besides classic presentation of intermit tent claudication, foot ulceration or diabetic foot may develop, signaling a significant degree of arterial occlu sion and high risk for future amputation. PAD is respon sible for 47% of amputations in diabetic patients [1]. In recent years, many imaging techniques, such as Doppler ultrasound, colour duplex ultrasound (CDU), Xray an giography, computed tomographic angiography (CTA), and magnetic resonance angiography (MRA), have been used to diagnose PAD, verifying the presence of a steno sis inside the vasculature of diabetes [24]. Although these methodologies are reliable and highly accurate, a noninvasive and inexpensive technique for early diagno sis and monitoring of PAD remains unsettled in hospitals and primary care settings. Modern ultrasound technique is welldeveloped, compact, fast, and free of irradiation, but they can only be correctly operated by a welltrained technician clinical physician. Therefore, this study pro *Corresponding authors. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1144 poses a noninvasive means of measurement and an automatic diagnostic algorithm to estimate the degree of PAD, purporting for screening and monitoring of lower limb vascular occlusion for highrisk groups like patients of diabetes, dyslipidemia, and hypertension, and the eld erly. Photoplethysmography (PPG) is an optical and nonin vasive technique, widely applied to measure oxygen sa turation, monitor heartbeat, and detect blood volume changes in the vascular bed. The advantages of PPG like compactness fast, and easytooperate [58] take the form of a portable instrument for routine screening. PPG sig nals can be obtained from multiplesite measurements including earlobe, forearm, fingertip, thumb, and toe [9 17]. A typical waveform of PPG is triphasic. Biphasic or monophasic form indicates a degree of arterial obstruc tion. In signal processing, timedomain and frequencydo main analysis has been applied to assist diagnosis of PAD. Previous research reports that frequency parame ters verify the differences between healthy and PAD subjects, using highfrequency and lowfrequency char acteristics [57,15]. The most common parameters of PPG waveforms are transit time, amplitude, and shape changes, which are helpful for evaluation of arterial com pliance single or bilateral measurements of lower limbs. Previous research reports that these timing parameters and shape characteristics reflect the degree of vascular stenosis, contributing to time delay, damping, amplitude reduction, and changes of resistance and compliance pro perties of occlusive peripheral arteries, with variations depending on a subject’s age and measurement sites [11, 13,15]. Therefore, an alternative for applying PPG to evaluate PAD is proposed, including bilateral or multi channel PPG techniques resulting from the findings of significant bilateral asymmetry and increased difference in timing and frequency parameters along with increased PAD severity [1315,18]. Multichannel measurements were used to study the agerelated changes in the char acteristics of PPG shape at the ear, finger, and toe sites. The largest changes with age were seen at the ear and finger sites for the systolic rising edge region, and the finger site for the dicrotic notch region. The smallest changes with age are for the toe site [911]. Accordingly, in this experiment, bilateral measurement at the right and left great toes was designed for PAD estimation. The objective of this study was to construct a fre quencybased Van der Pol (VDP) oscillator to evaluate the degree of PAD in patients with diabetes mellitus, and apply Burg autoregressive (AR) method as a parametric method for timedomain and frequencydomain signal processing [1924], purporting to smooth the spectra and find characteristic frequencies from the PPG signals re corded. Given a bilateral frequency spectrum, a second order VDP system has a step response that is character ized by damped oscillation [2528]. The oscillation over shoot (OS) is an approximating function of the damped factor, and determines the severity of PAD: normal con dition (Nor), mildtomoderate disease (MD), and severe disease (SD). With twentyone subjects, the test results demonstrate that the proposed method is efficient and accurate in the estimation of PAD in the feet of diabetics. 2. FREQUENCY SPECTRA ANALYSIS In signal processing of plethysmographic waveforms, higher frequency characteristics are preferred for fre quency analysis to identify normal pulse waves [57]. On the other hand, low frequency is used to characterize PAD in diabetic patients. Both frequencybased methods, such as Fourier transforms (FTs), and frequency spec trum methods, are the nonparametric techniques that are used to preprocess the PPG signals. They are disadvan tageous for their spectral leakage effect, because of the size of sampling window. Wavelet transform is a para metric method and provides better timefrequency reso lution than nonparametric methods, and would be a promising method to extract features from nonstationary biosignals. However, significant frequencies are extracted at specific wavelet coefficients with different wavelets, and several frequency bands are analyzed using the cas caded lowpass or highpass filters and downsampling/ upsampling operations through trial procedures. In literatures [1720], the Burg autoregressive (AR) method is able to derive the frequency spectra by fitting an AR model of a given specific order. Its advantages over abovementioned methods include smooth spectra, high frequency resolution, stable AR model, and compu tationally efficient. The Burg AR method is also a para metric method for timedomain and frequencydomain signal processes. It could smooth spectra in comparison with frequencybased methods to find the characteristic frequencies with distinguishing peak central/main spectra by fitting an AR model of a given specific order. The frequency spectra of PPG signals is expressed as a linear combination of the previous samples and the re sidual values, resi. For a discrete set of n sampling points, P coefficients are used to approximate the original data, i, 1, 2, 3,,in , where is a discrete frequency spec trum of the PPG signal. The residual value is assumed to be independent of the previous samples and is calculated by [20,21] 1 ,1,2,3,, P ii pip p resi n (1) where i is the frequency spectrum of the PPG signal, P is the order of the AR model, and , and p are the coefficients of the AR model. 1,2, 3,,pP Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1145 The optimal coefficients, p, are used to minimize the square error, , between the original data 2 i n Eres and the approximated data. For the forward and back ward linear prediction, the object is to minimize Fp and Bp, as shown in [22] 2 2 2 0 1 p nn n pitiiit p ipip tip fi with 0 p iit i fi (2) 22 2 0 01 0 np np ni pitiiitp ipi ti B bi with 0 p iit i bi (3) where t, t[p, n], is a linear weighted combination of p previous sampling data, and t, t[0, ntp], is a linear weighted combination of p subsequent sampling data. The sum of the residual energy (SORE) in stage p is Ep = Fp + Bp. The LevinsonDurbin recursion algorithm is used to minimize the SOREs to estimate the model coef ficients. The coefficients of p, p[1, P], are stored in a vector and an inverted order vector . Therefore, T 12 1,,, ,, , pp A 0,,,,, pPp V P p T 21 ,,1 the recursion formula is the following [22] 1pp AV (4) 1 11 0 1 2 10 1 0 np pp pp inp n pp ip i 2 ip bi FB ib i n (5) The concept of Burg AR method adjusts the parameter to minimize the SORE Ep using the final prediction error criterion (PEC) [23,24]. Then update the Ap and update fp(i) and bi(i) as 11 ppp fifi bip (6) 11 ppp bibi fip (7) The optimal coefficients p are chosen to minimize the squared error, whereas the forward and backward linear prediction equations attempt to minimize Fp and Bp. The entries of i represent samples of a discrete data, and P is the integer order of an AR model, used in estimating the power spectral density (PSD), log( i)/max{log( i)} dB/kHz, . This property of the AR model is used to smooth the frequency spectrum, which tends to favor the peaky spectra. Qualitatively different PSDs are observed by normalized process, which different values within the specific ranges of normalized logmagnitude and scaled frequency. Therefore, PSD can be used to identify the characteristic frequencies and magnitudes in each frequency spectra. 1,2, 3,,i In the physiological measurements, PPG signals were obtained using optical sensors (940 nm near infrared, spectral bandwidth: 45 nm) with 1 KHz sampling rate, as shown in Figure 1(a). The dashed line represents a nor mal subject and the solid line represents a diabetic pa tient [14,27,28]. It was found that the Fourier transform based estimation methods has some drawbacks, such as spectral leakage into sidelobes and characteristic spectral broadening, as shown in Figures 1(b) and (c). In order to find characteristic frequencies, the Burg AR method is used to estimate the PSD by fitting an AR model with given specific orders, P = 8, 12, and 16, to the PPG sig nals. The characteristic frequencies of the PPG signals easily found in the frequency spectra between 0 Hz and 500 Hz. The PSDs are shown in Figures 1(d)(f). The peaky spectra for diabetic subjects and the smooth fre quency spectra for normal subjects fall into different frequency bands, respectively. Depending on these bands, frequencybased parameters provide information for the estimation of PAD in diabetic patients. 3. THE ESTIMATION OF PERIPHERAL ARTERIAL DISEASE (PAD) SEVERITY 3.1. Feature Extraction Using the Van der Pol (VDP) Oscillator PAD can develop in the arteries of most visceral organs and extremities. Atherosclerosis in the brain and heart draws much attention for clinical diagnosis and interven tion; by contrast, lower limb PAD is underestimated by the patients, even in diabetics whose progress of vascular atherosclerosis is accelerated by the metabolic disease. PPG is a noninvasive technique capable to optically gen erate a plethysmograph at great toes, index fingers, and ears, and monitor blood pressure, blood oxygen satura tion (SaO2), and blood volume changes in an artery or a vascular bed of tissue. A PPG signal consists of AC (Al ternating Current) and DC (Direct Current) components. The AC component reveals physiological information, such as cardiac synchronous changes in the blood vol ume with each heartbeat and vasomotor activity [810, 15]. Timedomain and frequencydomain parameters have been used to detect the degree of PAD severity for diabetic patients with normal condition (Nor), mildto moderate disease (MD), and severe disease (SD). For the timedomain parameters, ECG Rpeak as a timing refer ence is used to obtain pulse rise time (RT), pulse transit time (PTT) from Rpeak to pulse foot (PTTf), pulse tran it time from Rpeak to pulse peak (PTTp), and PPG am s Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 Copyright © 2013 SciRes. 1146 Figure 1. (a) Timedomain PPG signal analysis for normal subjects and diabetic subjects, (b) and (c) PSD estimation using the frequency spectrum method, (d), (e), and (f) PSD estimation using the Burg method for an AR model with orders P = 8, 12, and 16 for normal subjects and diabetic subjects. plitude (plus foottoplus peak amplitude), as shown in Figure 2(a). The timing parameters, PTT and RT, are prolonged as due to the resistance of peripheral vessels increases, and the amplitude (AMP) and shape of pulse waves are smoothed [11]. Parameters, PTT and AMP, and shape differences vary with ages for individual sites, but the smallest variations with age occur at the great toe sites [915]. Thus, bilateral differences in the timing pa rameters, ∆PTTf, ∆PTTp, and ∆RT, provide time delay information for the estimation of PAD. Bilateral differ ences in frequency spectra are similar and are used to distinguish between the normal subjects and diabetic subjects. OPEN ACCESS This study uses frequencybased parameters to esti mate the degree of PAD. The PSDs are estimated using the Burg AR method. The Van der Pol (VDP) system is used to quantify the features of the right and leftside PSDs. It is an oscillator with nonlinear damping that is defined by a secondorder differential equation. The formulation has the form of an autonomous system with two state variables. The state equations are evaluated as below [2326]: d d rl (8) 2 d1 d lrrl (9) where r( ) and l( ) are the functions with respect to frequency , r( ) is the PSD of the rightside PPG sig nal, except for = 0, and l( ) is the PSD of the left side PPG signal, right and left are denoted by the sub scripts r and l. The parameter, , is a control parameter, and > 0 reflects the degree of nonlinearity of the sys tem [21]. When the term, ( r( ))2, becomes dominant, the VDP system becomes a nonlinear equation with posi tive damping. The dynamics of the system are stable and are restricted to a fixed point. Equation (9) gives the frequency of selfoscillation, as determined by a real parameter, , and demonstrates dissipation or damping. If the parameter, = 0, then Equation (9) reduces to that for a simple harmonic oscil lator. When the PSDs, r( ) and l( ), are different (not symmetrical), l( ) provides the damping in the VDP oscillator that results in selfsustained oscillations. The multiple peaks and amplitudes of these peaks demon strate the nonlinearity of the frequency spectrum. PSD l( ) is used to determine the response of the VDP oscillator. The VDP system demonstrates a damped si nusoidal response for a general secondorder system. The transient response consists of a sinusoidal oscillating waveform with an exponentially decaying amplitude. The sinusoidal frequency is called the damped frequency of oscillation. In order to evaluate the discrete frequency spectrum, the PSD is estimated using the Burg AR method. As suming a set of n points from the PSD, i, , (n = 500 in this study), the continuous VDP system can be modified as a discrete VDP system. Therefore, the discrete VDP system that is proposed for the estimation of PAD, as 1,2, 3,,in
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1147 rl i i i n (10) 21 lr rl ii i (11) where the index, i, is the integer scale of the frequency, . A bifurcation is a fundamental change in the nature of a solution. For different initial conditions, 1,2, 3,,i 1 ri and 1 li , the output, li, given by Equation (11), has a step response that is characterized by a damped oscillation in the frequency domain. It ex hibits a rich variety of nonlinear dynamic behaviors and generates the limit cycle for small values of , and de velops into relaxation oscillations when becomes large. If the VDP oscillators use different PSDs, the oscillations become unstable, as the amplitudes of spectra peaks in crease. The step response of the bilateral PSDs is shown in Figure 2. A secondorder step response is character ized by damped oscillations. From Figure 2(c), the per centage overshoot OS% is given by the following Equa tion [29]: max min min % cc OS c 100% (12) The percentage overshoot OS% defines the amount by which the oscillation overshoots the minimum value cmin = b. The term cmax is determined by curve fitting the function at the maximum value, as [29] 2 max 2 2 min exp1 cossin 1 exp 1(13) (14) cb b cb For the step response, substituting Equations (13) and (14) into Equation (12) gives 2 1 % exp1100%OS b (15) Equation (15) demonstrates that the percentage over shoot is a function only of the index, , and allows the determination of the OS%, given the index, . The in verse of the equation allows a solution for given the OS%. The inverse is given by 22 ln % ln % bOS bOS (16) The relationship between the index, , and the multi ple characteristic frequencies is clearly seen. The higher the value of , the more oscillatory is the response in the frequencydomain from Equations (10), (11), (15), and (a) (b) (c) Figure 2. (a) ECG and PPG signals and pulse landmarks, (b) The bilateral PSD functions for rightside and leftside PPG signals, (c) The damped sinusoidal response. (16). The increase in the resistance of peripheral vessels produces multiple characteristic frequencies, so this study uses the index, , is used to estimate the degree of PAD. The various ranges of the index are obtained from specific subjects, including those with normal con dition (Nor), mildtomoderate disease (MD), and severe disease (SD). 3.2. Physiological Measurement PPG signals were collected from twentyone subjects in a hospital (ChiMei Medical Center, Department of In fectious Disease and Immunology, Tainan City, Taiwan). The subjects, aged from 24 to 65 years, were divided into three groups: Nor, MD and SD, in terms of diabetic dis ease. In a clinical examination, the ABPI (AnkleBra Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1148 chial Pressure Index) was used as an early detection method to decide the degree of PAD. The degrees were categorized by clinical manifestations, ABPI values, and angiographic findings. For preliminary PAD estimation using the APBI, the indices signify ABPI ≥ 0.9 for nor mal subjects and diabetic subjects with MD, ABPI < 0.9 (at least one leg) for diabetic subjects with SD, as shown in Table 1 [2830,32]. Despite the noninvasive nature of the examination to assess death rate and cardiovascular diseases in highrisk patients, the accuracy of the ABPI measurement is reduced in patients with calcified blood vessels, caused by conditions such as diabetes and chro nic renal failure. However, measurement of the ABPI must be repeated several times (>10 minutes) and is of limited use in routine vascular screening in a primary care setting [13]. Bilateral timing parameters could offer a quick as sessment for the screening of PAD in primary care. Bi lateral measurements simultaneously acquire PPG sig nals from the right and left big toes. The absolute timing differences are used to reference one side of the body (right side) with the contralateral side (left side), in or der to calculate the parameters, PTTf, PTTp and RT, so various ranges are obtained for specific groups. Each parameter has a mean value and a specific range between maximum (Max) and minimum (Min) values. A com parison of the bilateral differences demonstrates that these parameters increase as the severity of the disease increases. Therefore, clinical physicians would consider timing differences to be a good reference for PAD as sessment. The three groups are divided into ten normal subjects (No. 1  No. 10), eight MD diabetic subjects (No. 11  No. 18) and three SD diabetic subjects (No. 19  No. 21), as shown in Table 1. For these twentyone subjects, Equations (12) and (16) were used to compute the index, . The indices, , also represent the specific ranges for diabetic patients, non diabetic patients with PAD and normal subjects as a con trol, as shown in Figure 3. Severe PAD causes more asymmetrical PPG signals and the damped oscillation increases as PAD increases, which is caused by the re sistance of the peripheral vessels. Therefore, the index, , also increases in diabetic patients with increased PAD. The index, , is less than 0.64 for normal subjects with an ABPI ≥ 0.9, from 0.64 to 0.65 in diabetic patients with an ABPI ≥ 0.9 and greater than 0.65 for diabetic patients with an ABPI < 0.9. The standard indices, , are established, as shown in Table 1. The sensitivity is greater than 85.7%, and positive predictivity is also greater than 80.0% to quantify the performance of the proposed mea surement method. These specific ranges also confirm the degree of PAD noted by clinical physicians. The trend in the degeneration of PAD is clearly demonstrated. De pending on the severity of the PAD assessment, these specific ranges provide key information for the evalua tion of the degree of PAD. 4. EXPERIMENTAL RESULTS AND DISCUSSION 4.1. Experimental Setup Using the bilateral noninvasive measurement, two opti cal sensors (reflection mode), consisting of light sources, photodetectors, transimpedance amplifiers, and high pass filters, were placed at the right and left great toes. Near infrared (NIR) has large differences in the extinc tion coefficients of deoxyhaemoglobin and oxyhaemo globin. Thus, the light source 940 nm NIR was chosen in this optoelectronic design. The reflection mode, includ ing light source and photodetector, was positioned side by side with 5 mm spacing, and the light is directed down into the skin and is backscattered from the skin adjacent to the photodetector [10]. Two optoelectronic probes (circle shape, diameter: 20.0 mm, height: 8.2 mm) were synchronized using a data acquisition controller. PPG signals are captured at a sampling rate of 1 kHz for 15 minutes [30]. A DAQ card (National Instruments DAQ Card, 16 Channels, 1.25 MS/s) served as an ana logtodigital (A/D) converter between the optical meas urement system and a computer. Locating each pulse foot (PF)pulse foot (PF) interval of PPG signal, 800 sampling data were acquired within a sampling window. In this study, we used the Burg AR method to estimate the PSDs of PPG signals. The suitable AR order can be used to identify the peaky spectra. Figure 4 shows the variation of residual energy versus AR model order. We use model orders, from order 1 to order 30, to calculate the SORE. To obtain the suitable order, Akaike’s final prediction criterion is used to select the model order [22,24]. For considering convergent condition, we con sidered SORE ≤ 10−1 to stop the Burg AR algorithm. Then, the Burg method using model order P = 8 with optimal coefficients was used to estimate the PSDs. The VDP oscillator shows the step responses in the fre quencydomain, as shown in Figure 5. The oscillations of the step responses are smaller than those for normal sub jects. The higher the value of , the more oscillatory the response are important information for the degree of PAD assessment, which may be useful for determining the trends of PAD in routine examination. The proposed algo rithms were developed on a PC AMD Celeron (R) CPU 2.40GHz with 2.39 GHz, 224 MB RAM GHz with 1.75 GB RAM and Matlab software. Portable optical meas urement was used to obtain the PPG signals in the labo ratory and medical center. To demonstrate the effective ness of the proposed method, the database of 21 subjects were used to design the gold standard of PAD assess ment using the VDP oscillator from selected 6 subjects Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 Copyright © 2013 SciRes. 1149 Table 1. The ABPI, bilateral differences in timing parameters [30], and index for PAD estimation. ABPI ≥ 0.9 (10) 0.5 ≤ ABPI < 0.9 (8) ABPI < 0.5 (3) Normal subject Diabetic subject Diabetic subject Subject Category Parameter Normal (Nor) Mildtomoderate disease (MD) Severe disease (SD) PTTf (ms) Mean MinMax 2.58 0.3  7.4 9.2 5.1  23.7 29.2 23.6  34.8 PTTP (ms) Mean MinMax 7.3 0.4  22.3 22.6 14.3  56.5 52 46.2  57.8 RT (ms) Mean MinMax 6.84 1.3  15.6 14.6 3.4  32.3 23.4 11.5  35.3 Index Mean MinMax 0.616 0 .014 0.597  0.636 0.645 0.002 0.640  0.647 0.668 0.012 0.652  0.681 Note: The values of PTTf, PTTp, and RT are absolute values. (1) fRf PTT PTTPTT Lf , (2) Rp Lp PTT PTTPTT , (3) RL RT RTRT , where suffix words R and L are defined right and left legs, (4) ABPI: it is calculated using the highest of the right and left ankle systolic blood pressure divided by the highest of the right and left arm brachial systolic blood pressure [13]. Figure 3. The specific ranges of the index in normal subjects and diabetic patients with PAD. Figure 5. The step responses of the VDP Oscillator for severe disease (SD), mildtomoderate disease (MD), and normal con dition (Nor). and diabetic patients with PAD were tested using physio logical measurement. There were two normal subjects (No. 1 and No. 2), four diabetic patients, two of which were MD subjects (No. 11 and No. 12) and two of which were SD subjects (No. 19 and No. 20). The preliminary diagnosis and ABPI classification were performed by clinical physicians and the results are shown in Table 2. However, the ABPI is significantly lower for normal subjects and diabetic patients with MD (ABPI ≥ 0.9). For bilateral PPG measurement, two optical sensors were placed on the right and left great toes and PPG signals were obtained from the three groups, as shown in Figure 6. The transit time is prolonged for an occluded leg and gradually increases as the disease becomes more severe, as seen in the interval between the dashed lines. For ex ample, the bilateraltiming differences are Figure 4. Variation of residual energy versus AR model order. (three groups: Nor, MD, and SD) were used to test, as shown in Table 2. The result demonstrates the computational efficiency and accurate diagnosis achieved by this study. 4.2. Experimental Results Normal SubjectNo. 1: PTTf = 2.733, PTTp = 1.977, RT = 3.133, In order to verify the proposed method, normal subjects OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1150 Table 2. The experimental results of PAD estimation. Mean values of bilateral differences Subject No. PTTf (ms) PTTp (ms) RT (ms) ABPI RLeg ABPI LLeg Clinical physician decision Bilateral timing difference ANFIS decision (output) Fuzzy logic decision (output) The proposed method (index ) 1 2.7330% 7.05% 1.9970% 16.26% 3.1330% 0.19%1.06501.1138Nor Nor Nor (0.1988) Nor (0.1750) Nor (0.6199) 2 1.1730% 2.78% 5.8700% 4.05% 4.6950% 3.97%1.14041.1333Nor Nor Nor (0.1990) Nor (0.1620) Nor (0.6223) 3 2.3330% 12.07% 3.6660% 5.44% 13.330% 3.97%1.14401.1101Nor Nor Nor (0.1978) *Failure Nor (0.6336) 11 15.321% 16.92% 17.286% 10.50% 1.9640% 4.05%1.17881.2357MD MD MD (0.4000) MD (0.4000) MD (0.6451) 12 6.3850% 13.44% 16.923% 16.76% 17.556% 0.39%1.27351.3076MD MD MD (0.3988) MD (0.3830) MD (0.6478) 13 23.667% 13.62% 33.333% 7.58% 32.833% 16.63%1.0714 1.0446MD *SD MD (0.4000) MD (0.4020) MD (0.6476) 19 33.072% 10.06% 48.500% 14.18% 15.428% 8.57%1.08170.8941SD SD SD (0.6000) SD (0.6000) SD (0.6733) 20 23.606% 13.64% 57.700% 6.05% 35.000% 10.83%1.04420.8945SD SD SD (0.5988) *MD (0.4000) SD (0.6810) Note: ANFIS and Fuzzy Logic: 1 output variable with 3 triangular membership functions, the center values of the functions are 0.2, 0.4, and 0.6 for Nor, MD, and SD. Figure 6. Timedomain PPG signals of the three groups of patients. MD SubjectNo. 11: PTTf = 15.321, PTTp = 17.286, RT = 1.964, SD SubjectNo. 19: PTTf = 33.072, PTTp = 48.500, RT = 15.428. The shape of the PPG signals for a SD subject with a unilateral occluded leg is substantially different. The PPG signals at the right and left sites are asynchronous in the timedomain. According to the standard established in Table 1 [30], the degree of PAD from the preliminary estimate is confirmed, as shown in Table 2. However, the ABPI and bilateral timing differences are not the only parameters that determine classification. Using twenty bilateral PPG signals, the diagnostic re sults for six subjects are shown in Figure 7, where the horizontal axis is the number of PPG signals and the ver tical axis is the index, , for each PAD assessment. The Burg AR method is used to estimate the frequency spec tra from the bilateral PPG signals. Using the bilateral frequency spectra, the VDP oscillator shows the dynamic characteristics to be step damped oscillation responses for diabetic patients with PAD. In this study, the index, , is a parameter that represents the degree of PAD that is Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1151 Figure 7. Different indexes for separating the diabetic PAD from the normal subjects. assessed using the VDP oscillator, for the Nor, MD and SD groups. A comparison of the results shown in Figure 7, demonstrates that the index, , is a significant pa rameter for the separation of diabetic patients with PAD from normal subjects. According to the standard values in Table 1, the values are centered at 0.645 (±0.002) for MD subjects and are greater than 0.65 for SD subjects (0.668 ± 0.012). Within these specific ranges, the index, , can also be used to estimate the degree of PAD. Experimental tests show that the proposed VDP oscil lator can be used to assess the degree of severity and can be used to monitor the trends in PAD degeneration at the borders between Nor and MD (subjects No. 3, No. 6, and No. 12), or between MD and SD (subject No. 21). In this study, the index, , is used to monitor the trends in the degeneration of PAD to allow early detection or screen ing and monitoring. Diabetes mellitus is a chronic dis ease that may require therapy to prevent further problems. In particular, for the feet of diabetics, adequate treatment may improve the risk profile of chronic complications. This research provides a potential method for early de tection, monitoring, and prevention of PAD by early in tervention to control its risk factors. 4.3. Discussions Table 2 shows a comparison of the experimental results using the proposed method, bilateral timing difference, Fuzzy logic decision, and an adaptive network based fuzzy inference system (ANFIS). Using the bilateral timing difference, the timing parameters, PTTf, PTTp and RT, are used to estimate the degrees of PAD sever ity. For example, subjects, No. 2 and No. 12, were con firmed as a healthy subject and a MD diabetic patient, respectively. According to the Table 1, the timing pa rameters, ∆PTTf = 6.3850 and ∆PTTp = 16.923, signify an overlap between Nor (standard: ∆PTTf = 0.3  7.4 and ∆PTTp = 0.4  22.3) and MD (standard: ∆PTTf = 5.1  23.7 and ∆PTTp = 14.3  56.5). If there is overlapping of the ranges of the timing parameters, the inferences are affected by the timing reference (heart rate) and meas urement errors. This demonstrates that timing parameters alone are insufficient for the estimation of the degree of PAD. The clinical physicians decided the possible degree using a combination of the variances in the three timing parameters. However, this examination method required offline analysis to obtain a diagnosis. According to the Table 1, each timing parameter has a mean value and a specific range between the minimum (Min) and maximum (Max) values. A triangular mem bership function is parameterized by its triplet values (Min, Mean, Max) and is decomposed into three fuzzy partitions, defined as Nor, MS and SD, as shown in Fig ure 8. Therefore, the Fuzzy logic decision has three input variables with nine triangular membership functions and 1 output variable with three triangular membership func tions. The fourteen linguistic Fuzzy rules for the three degrees are determined by professional physicians over many examinations, as shown in Table 1. Inference uses the centre of the mean of the maximal defuzzifier. For adaptive and selforganizing applications, the literature [30,31,33,34] also cites the gradient descent method, the leastsquare algorithm and the modified leastsquare al gorithm for training a network structure and network parameters. The ANFIS structure is immediately deter mined after the presentation of each inputoutput pair and the inference rules. Updating of the parameters is then performed after the network structure has been de cided. In this study, the leastsquare algorithm is used to tune network parameters. For the same subjects, the overall results are shown in Table 2. The sensitivity is 76.2% using Fuzzy logic deciion. The ANFIS has more s Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1152 Figure 8. The value of certainty degree versus index, , and bilateral timing parameters for separating the diabetic PAD. than 90% of sensitivity, but the adaptive mechanism is difficult implement in hardware devices, because of needs to assign and update network parameters [35]. The experimental results validate the effectiveness of the proposed VDP oscillator. A comparison with the pro posed method and other examination methods, the results of the 21 patients demonstrate the 85.7% of sensitivity and allow early detection between Nor and MD or be tween MD and SD. However, the trends of the index, , increases as the PAD gradually becomes severity. The proposed screening method provides a finding to evalu ate the orders of PAD for a routine examination. Its ad vantages are summarized as follows. The Burg AR method overcomes spectral leakage and smoothes the spectra to determine the characteristic frequencies. The proposed algorithm for the VDP oscillator is eas ily implemented in hardware devices, such as em bedded system (ES) and fieldprogrammable gate ar ray (FPGA) chips. The proposed VDP oscillator can also be imple mented using analog electronic circuits. The proposed method potentially allows the use of a portable monitor for the estimation of diabetic foot PAD in daily homecare. 5. CONCLUSION The Van der Pol (VDP) oscillator was able to detect pho Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 1153 toplethysgraphic signals to estimate severity of periph eral arterial disease (PAD) in patients with diabetes mel litus. In addition, bilateral noninvasive optical sensors were helpful for acquiring the PPG signals at the right and left great toes, combined with the Burg AR method for determining the frequency spectra of the right and leftside PPG signals. Using the bilateral frequency spec tra, the VDP oscillator shows the step damped oscillation responses for diabetic patients with mild and severe se verity, respectively, with the increasing amplitudes of damped oscillation as PAD severity worsening. The os cillation amplitude is expressed as an index to differen tiate between diabetic patients with PAD and normal subjects. The numerical experiments reveal specific ranges that allow confirmation of the degree, which allows tracking the trends of PAD for screening and monitoring of highrisk patients in primary settings as well as in hos pitals. For an improved diagnostic design, the proposed method can be combined with a fuzzy inference system to derive an index for the diagnosis of PAD. The pro posed method potentially allows for constructing a port able biomonitor and the use of telemedicine in a variety of clinical environments. 6. ACKNOWLEDGEMENTS This work is supported in part by the National Science Council of Tai wan under contract number: NSC992221E244007 (August 1 2010 October 31 2011). The authors would like to thank Dr. ChianMing Li for providing his valuable suggestion and help on experiments. REFERENCES [1] I. S. Muller, M. L. Bartelink, W. J. C. Grauw, H. J. M. Hooden, W. H. E. M. Gerwen and G. E. H. M. Rutten, “Foot Ulceration and Low Limb Amputation in Type 2 Diabetic Patients in Dutch Primary Health Care,” Dia betes Care, Vol. 25, No. 3, 2002, pp. 570574. http://dx.doi.org/10.2337/diacare.25.3.570 [2] J.Y. David, S. A. Jones and D. P. Giddens, “Modern Spectral Analysis Techniques for Blood Flow Velocity and Spectral Measurements with Pulsed Doppler Ultra sound,” IEEE Transactions on Biomedical Engineering, Vol. 38, No. 6, 1991, pp.589596. http://dx.doi.org/10.1109/10.81584 [3] C. Loewe, M. Schoder, T. Rand, U. Hoffmann, J. Sailer, T. Kos, J. Lammer and S. Thurnher, “Peripheral Vascular Occlusive Disease: Evaluation with ContrastEnhanced MovingBed MR Angiography versus Digital Subtraction Angiography in 106 Patients,” American Journal of Ro entgenology, Vol. 179, No. 4, 2002, pp. 10131021. http://dx.doi.org/10.2214/ajr.179.4.1791013 [4] O. Pablo Vasquez, M. Marco Munguia and B. Manders son, “Arteriovenous Fistula Stenosis Detection Using Wavelets and Support Vector Machines,” 31st Annual International Conference of the IEEE EMBS, Mineapolis, 26 September 2009, pp. 12981301. [5] M. Nitzan, A. Babchenko and B. Khonokh, “Very Low Frequency Variability in Arterial Blood Pressure and Blood Volume Pulse,” Medical & Biological Engineering & Computing, Vol. 37, No. 1, 1999, pp. 5458. http://dx.doi.org/10.1007/BF02513266 [6] E. D. Übeyli, D. Cvetkovic and I. Cosic, “AR Spectral Analysis Technique for Human PPG, ECG and EEG Signals,” Journal of Medical Systems, Vol. 32, No. 3, 2008, pp. 201206. http://dx.doi.org/10.1007/s1091600791237 [7] Y. M. Akay, M. Akay, W. Welkowitz, S. Lewkowicz and J. L. Semmlow, “Non Invasive Acoustical Detection of Coronary Artery Disease: A Comparative Study of Signal Processing Methods,” IEEE Transactions on Biomedical Engineering, Vol. 40, No. 6, 1993, pp. 571578. http://dx.doi.org/10.1109/10.237677 [8] J. Allen, “Photoplethysmography and Its Application in Clinical Physiological Measurement,” Physiological Meas urement, Vol. 28, No. 3, 2007, pp. R1R39. http://dx.doi.org/10.1088/09673334/28/3/R01 [9] J. Allen and A. Murray, “AgeRelated Changes in the Characteristics of the Photoplethysmographic Pulse Shape at Various Body Sites,” Physiological Measurement, Vol. 24, No. 2, 2003, pp. 297307. http://dx.doi.org/10.1088/09673334/24/2/306 [10] R. Erts, J. Spigulis, I. Kukulis and M. Ozols, “Bilateral Photoplethysmography Studies of the Leg Arterial Steno sis,” Physiological Measurement, Vol. 26, No. 5, 2005, pp. 865874. http://dx.doi.org/10.1088/09673334/26/5/022 [11] J. Allen and A. Murray, “Variability of Photoplethysmo graphy Peripheral Pulse Measurements at the Ears, Thumbs, and Toes,” IEE Proceedings on Science, Meas urement and Technology, Vol. 147, No. 6, 2000, pp. 403 407. http://dx.doi.org/10.1049/ipsmt:20000846 [12] P. A. Bonham, T. Kelechi, M. Mueller and J. Robison, “Are Toe Pressures Measured by a Portable PhotoPhle thysmograph Equivalent to Standard Laboratory Tests,” Journal of Wound, Ostomy and Continence Nurses Soci ety, Vol. 37, No. 5, 2010, pp. 475486. [13] J. Allen, K. Overbeck, A. F. Nath, A. Murray and G. Stansby, “A Prospective Comparison of Bilateral Pho toplethysmography versus the AnkleBrachial Pressure Index for Detecting and Quantifying Lower Limb Pe ripheral Arterial Disease,” Journal of Vascular Surgery, Vol. 47, No. 4, 2008, pp. 794802. http://dx.doi.org/10.1016/j.jvs.2007.11.057 [14] C.H. Lin, “Assessment of Bilateral Photoplethysmogra phy for Lower Limb Peripheral Vascular Occlusive Dis ease Using Color Relation Analysis Classifier,” Com puter Method and Program in Biomedicine, Vol. 103, No. 3, 2011, pp. 121131. [15] J. Allen, C. P. oates, T. A. Lees and A. Murray, “Pho toplethysmography Detection of Lower Limb Peripheral Arterial Occlusive Disease: A Comparison of Pulse Tim ing, Amplitude and Shape Characteristics,” Physiological Measurement, Vol. 26, No. 5, 2005, pp. 811821. http://dx.doi.org/10.1088/09673334/26/5/018 [16] M. Fallahpour, D. Megias and M. Ghanbari, “Reversible Copyright © 2013 SciRes. OPEN ACCESS
J. X. Wu et al. / J. Biomedical Science and Engineering 6 (2013) 11431154 Copyright © 2013 SciRes. 1154 OPEN ACCESS and Highcapacity Data Hiding in Medical Images,” IET Image Processing, Vol. 5, No. 2, 2011, pp. 190197. http://dx.doi.org/10.1049/ietipr.2009.0226 [17] D. M. Vázquez, J. J. Rubio and J. Pacheco, “A Charac terization Framework for Epileptic Signals,” IET Image Processing, Vol. 6, No. 9, 2012, pp. 12271235. http://dx.doi.org/10.1049/ietipr.2012.0037 [18] H. Heidrich, R. Wenk and P. Hesse, “Frequency of As ymptomatic Peripheral Arterial Disease in Patients En tering the Department of General and Internal Medicine of a General Care Hospital,” Vasa, Vol. 33, No. 2, 2004, pp. 6367. http://dx.doi.org/10.1024/03011526.33.2.63 [19] I. Kauppinen, J. Kauppinen and P. Saarinen, “A Method for Long Extrapolation of Audio Signals,” Journal of the Audio Engineering Society, Vol. 49, No. 12, 2001, pp. 11671180. [20] A. Broadman, F. S. Schlindwein, A. P. Rocha and A. Leite, “A Study on the Optimum Order of Autoregressive Models for Heart Rate Variability,” Physiological Meas urement, Vol. 23, 2002, pp. 324336. [21] K. Roth, I. Kauppinen, P. A. A. Esquef and V. Valimaki, “Frequency Warped Burg’s Method for ARModeling,” IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, 1922 October 2003, pp. 58. [22] C. Collomb, “Linear Prediction and LevinsonDurbin Al gorithm,” 2009. http://www.emptyloop.com/technotes/A%20tutorial%20o n%20linear%20prediction%20and%20LevinsonDurbin.p df [23] N. Kannathal, U. Rajendra Acharya, P. Joseph and E. Y. K. Ng, “Analysis of EEG Signals with and without Re flexology Using FFT and Auto Regressive Modeling Techniques,” Chinese Journal of Medicine, Vol. 1, No. 1, 2006, pp. 1220. [24] M. Akay, J. L. Semmlow, W. Welkowitz, M. D. Bauer and J. B. Kostis, “Detection of Coronary Occlusions Us ing Autoregressive Modeling of Diastolic Heart Sounds,” IEEE Transactions on Biomedical Engineering, Vol. 37, No. 4, 1990, pp. 366373. http://dx.doi.org/10.1109/10.52343 [25] S. B. Waluya and W. T. van Horssen, “On the Periodic Solutions of a Generalized Nonlinear Van der Pol Os cillator,” Journal of Sound and Vibration, Vol. 268, No. 1, 2003, pp. 209215. http://dx.doi.org/10.1016/S0022460X(03)002517 [26] Q. S. Bi, “Dynamical Analysis of Two Coupled Paramet rically Excited Van der Pol Oscillators,” International Journal of Nonlinear Mechanics, Vol. 39, No. 1, 2004, pp. 3354. http://dx.doi.org/10.1016/S00207462(02)001269 [27] G. M. Mahmoud and A. A. M. Farghaly, “Chaos Control of Chaotic Limit Cycles of Real and Complex Van der Pol Oscillators,” Chaos, Solitons and Fractals, Vol. 21, No. 4, 2004, pp. 915924. http://dx.doi.org/10.1016/j.chaos.2003.12.039 [28] Y. Abbas, J. Ann and A. Merna, “A Multistage Adomian Decomposition Method for Solving the Autonomous Van Der Pol System,” Australian Journal of Basic and Ap plied Sciences, Vol. 3, No. 4, 2009, pp. 43974407. [29] N. S. Nise, “Control Systems Engineering,” 4th Edition, John Wiley & Sons, INC., 2004. [30] Y. C. Du and C. H. Lin, “Adaptive NetworkBased Fuzzy Inference System for Assessment of Lower Limb Peri pheral Vascular Occlusive Disease,” Journal of Medical System, Vol. 36, No. 1, 2012, pp. 301310. http://dx.doi.org/10.1007/s1091601094761 [31] C. H. Lin, Y. F. Chen, Y. C. Du, J. X. Wu and T. S. Chen, “Chaos Synchronization Detector Combining Radial Ba sis Network for Estimation of Lower Limb Peripheral Vascular Occlusive Disease,” Lecture Notes in Computer Science, Vol. 6165, 2010, pp. 126136. http://dx.doi.org/10.1007/9783642139239_13 [32] P. A. Bonham, T. Kelechi, M. Mueller and J. Robison, “Are Toe Pressures Measured by a Portable Photophle thysmograph Equivalent to Standard Laboratory Tests?” Journal of Wound Ostomy & Continence Nursing, Vol. 37, No. 5, 2010, pp. 475486. http://dx.doi.org/10.1097/WON.0b013e3181eda0c5 [33] P. Angelov, E. Lughofer and X. Zou, “Evolving Fuzzy Classifiers Using Different Model Architectures,” Fuzzy Sets and Systems, Vol. 159, No. 23, 2008, pp. 31603182. http://dx.doi.org/10.1016/j.fss.2008.06.019 [34] J. J. Rubio, “SOFMLS: Online SelfOrganizing Fuzzy Modified Least Squares Network,” IEEE Transactions on Fuzzy Systems, Vol. 17, No. 6, 2009, pp. 12961309. http://dx.doi.org/10.1109/TFUZZ.2009.2029569 [35] C.H. Lin and G.W. Lin, “FPGA Implementation of Fractal Patterns Classifier for Multiple Cardiac Arrhyth mias Detection,” Journal of Biomedical Science and En gineering, Vol. 5, No. 3, 2012, pp. 120132.
