Using high-dimensional features for high-accuracy pulse diagnosis

1.   Introduction

Pulse taking has been used as one of the most important diagnostic methods in traditional Chinese medicine (TCM). TCM doctors use their index, middle, and ring fingers to take pulses at a patient’s wrist. Based on the tactile sensations on the fingertips, doctors can determine patients’ physiological conditions. However, accurate pulse diagnosis often requires rich clinical experience ^[1]. In addition, a doctor’s finger tactile sensation might be affected by many subjective and objective factors. Therefore, modern computer-aided pulse diagnostic methods have been developed to help doctors to rapidly and precisely determine a patient’s physiological conditions.

Most computer-aided pulse diagnostic methods first extract pulse features and then use certain classifiers to identify diseases or pulse types. The most commonly used classifiers are support vector machine (SVM) and artificial neural network (ANN). For example, Wang et al. ^[2] used multiscale sample entropy (Multi-SampEn) to extract 4 pulse features and then used SVM to classify healthy persons and diabetic patients. Guo et al. ^[3] used the wavelet packet transform (WPT) to extract 7 pulse features and then used SVM to hierarchically discriminate cholecystitis patients and nephrotic patients out of the rest of the population. Wang et al. ^[4] designed a multichannel sensor fusion device to create a 1572-dimension fusion feature vector and then used SVM to classify healthy persons and diabetic patients. Zhang et al. ^[5] applied the Jin’s pulse diagnosis method to extract 26 features and then used a cubic SVM algorithm to classify healthy individuals and lung cancer patients.

Other than using pulse features only, some studies integrated parameters from the pulse-taking devices to enhance classification performance. For example, Zhang et al. ^[6] used 5 Doppler ultrasonic diagnostic parameters and 16 pulse features extracted from WPT and then used SVM to discriminate the cholecystitis group from the healthy group. Chen et al. ^[7] used 4 Doppler ultrasonic diagnostic parameters and 2 pulse features extracted from an auto-regressive model and then used SVM to classify appendicitis patients and healthy persons. SVM was developed mainly for binary classification. If there were more than two groups, SVM usually compared each group alone against the whole set of other groups ^[3,8,9].

Artificial neural network (ANN) also has been used in the medical fields to help identify diseases. For example, Tang et al. ^[10] used a visual scale to extract 8 pulse features and then used an ANN to differentiate essential hypertension from normotension. Du and Stephanus ^[11] generated 4 features from photoplethysmography signals and then used an ANN to classify the degree of arteriovenous fistula stenosis in hemodialysis patients.

Contrary to the above hand-designed feature extractors, advanced machine learning techniques can extract features automatically ^[12]. Machine learning techniques can extract informative abstraction and features directly from the input pulse images. For example, Li et al. ^[13] used a convolution neural network (CNN) to classify 5 cardiovascular diseases and healthy persons. Although CNN has the advantage of excluding additional feature extraction tools, it fails to reveal the underlying physical meaning or implications hidden in pulse discrimination.

In TCM, pulses are classified into 28 single pulse types based on four main elements, pulse depth, pulse rate, pulse shape, and pulse strength ^[1,14]. Pulse depth describes the vertical position of a pulse. Pulse rate describes the number of pulses per unit time. Pulse shape describes the width and length of a pulse. Pulse strength describes the forcefulness of a pulse. If the normal pulse, i.e., the pulse of a healthy person, is included, totally, there are 29 single pulse types. A patient’s pulse might be a composition of several single pulse types, which is called a complex pulse ^[15]. In this study, a complex pulse is represented in the form of (single pulse A + single pulse B). Complex pulses might carry more physiological information than single pulses. For example, the (floating + rough) pulse might be related to a cold; the (slow + slippery) pulse might be related to asthma.

Some researchers have attempted to classify multiple TCM pulse types. For example, Wang and Cheng ^[16] extracted 13 features directly from the time-domain pulse waveforms and then used Bayesian networks to classify 5 TCM pulse types. Xu et al. ^[17] extracted 4 pulse features from the time-domain pulse waveforms and then used the Lempel-Ziv analysis to classify 7 TCM pulse types. Xu et al. ^[18] extracted 17 pulse features from the time-domain pulse waveforms and then used a fuzzy ANN to classify 16 TCM pulse types. Shu and Sun ^[19] used the gamma density function to obtain 7 pulse features to classify 13 TCM pulse types.

Most of the above-mentioned pulse classification methods used low-dimensional feature vectors to classify certain diseases or pule types. However, a feature generation method might only be effective for extracting certain features but not for others. For example, time domain features could not reveal the same pulse complexity as approximate entropy (ApEn) ^[20] and Multi-SampEn ^[21]. Some important pulse information is even too subtle to be extracted. To understand the underlying physical meaning or implications hidden in pulse discrimination, this study attempted to generate a high-dimensional feature vector from the time, spatial, and frequency domains to cover as much pulse information as possible to increase pulse classification accuracy.

2.   Data collection

In this study, pulse signals were measured using an ANSWatch wrist monitor (Taiwan Scientific Corp., Taipei, Taiwan). The measurement results given by ANSWatch include pulse types. Ten different wrist pulse types were collected, including floating pulse, rapid pulse, moderate pulse, normal pulse, sunken pulse, (rapid + stirred) pulse, (long + replete) pulse, (rapid + sunken) pulse, (normal + replete) pulse, (sunken + rapid + stirred) pulse. In TCM, floating pulse and sunken pulse are recognized by the depth of a pulse. Rapid pulse and moderate pulse are recognized by the rate of a pulse. Stirred pulse and long pulse are recognized by the shape of a pulse. Replete pulse is recognized by the strength of a pulse. The sampling rate was 500 Hz. Fifty 6-second data sets were collected for each pulse type. Since there were 10 pulse types, in all total, there were 500 data sets. The device used a pressure sensor to measure pulse signals, and they were converted to digital signals using a Micro Control Unit (MCU), ranging from 0 to 4095 digits. The digit is proportional to the pressure measured by the device sensor. However, because the actual pressure is unknown, “pressure index” is used as the unit of the vertical axis coordinate, as shown in Figure 1.

Figure 1.  Ten pulse types collected in this study.

DownLoad: Full-Size Img PowerPoint

It is essential to record high quality pulse signals for precise computer-aided diagnosis. However, a subject’s movement or respiration might easily cause baseline drift in pulse signals. In this study, to achieve high classification accuracy, Zhang et al.’s iterative sliding window algorithm was used to remove baseline drift ^[5]. First, a cubic spline was created to fit the local minimums of the pulse signals, as shown in Figure 2(a). Then, the cubic spline was subtracted from the original pulse signals to remove the baseline drift, as shown in Figure 2(b). After the baseline drift was removed, the following feature extraction methods were applied to extract pulse features.

Figure 2.  (a) Before baseline removal; (b) after baseline removal.

DownLoad: Full-Size Img PowerPoint

2.1.   Waveform shape

Seven features were extracted directly from the time series waveforms, including mean, standard deviation, variance, root mean square, average of the pulse intervals, standard deviation of the pulse intervals, and average of the peak amplitudes, and they are corresponding to “Mean”, “SDV”, “Variance”, “RMS”, “Avg. of the pulse interval”, “SDV of the pulse interval”, and “Avg. of the peak amp.” in Table 1, respectively.

Table 1.  Features of the 500 data sets and the loadings of the first 3 PCs.

Features		Original data set				Principal component
		1	2		500	PC1	PC2	PC3
Waveform shape	Mean	585.992	559.883	…	551.754	-0.151	0.197	-0.041
	SDV	550.465	545.270	…	538.055	-0.109	0.017	0.098
	Variance	303012	297320	…	289503	-0.133	0.043	0.197
	RMS	803.923	781.464	…	770.600	-0.148	0.049	0.209
	Avg. of the pulse interval	0.720	0.741	…	0.743	-0.129	-0.015	0.252
	SDV of the pulse interval	0.029	0.072	…	0.053	-0.143	0.204	-0.068
	Avg. of the peak amp.	1832.130	1824.517	…	1780.895	-0.129	0.183	-0.124
Power spectral density	Main freq.	1.465	1.465	…	1.465	-0.135	0.199	-0.113
	Main freq. PSD	260585	226554	…	238887	-0.135	0.212	-0.112
	BER1 (0-0.488HZ)	28.609	23.850	…	22.926	-0.134	0.215	-0.113
	BER2 (0.488-1.465HZ)	42.456	44.513	…	47.382	-0.133	0.216	-0.109
	BER3 (1.465-2.441HZ)	9.975	12.375	…	11.247	-0.133	0.212	-0.111
	BER4 (2.441-3.418HZ)	11.778	9.472	…	10.190	-0.009	-0.076	0.004
	BER5 (3.418-4.395HZ)	5.184	7.071	…	5.596	0.178	0.043	-0.050
	BER6 (4.395-5.371HZ)	0.903	1.241	…	1.275	-0.016	0.063	-0.012
	BER7 (5.371-6.348HZ)	0.634	0.708	…	0.853	0.141	-0.193	0.122
	BER8 (6.348-7.324HZ)	0.337	0.533	…	0.380	-0.095	0.086	-0.181
	BER9 (7.324-8.301HZ)	0.087	0.179	…	0.109	-0.123	0.139	-0.051
	BER10 (8.301-9.277HZ)	0.028	0.041	…	0.035	-0.042	0.058	-0.138
	BER11 (9.277-10.254HZ)	0.009	0.018	…	0.007	-0.102	0.174	-0.021
	Mean freq.	1.611	1.714	…	1.677	-0.013	0.001	0.268
Gamma density function	$\rm{ \mathsf{ α} }$	1.813	1.210	…	1.786	-0.045	0.132	0.057
	$\rm{ \mathsf{ β} }$	3.850	2.059	…	4.625	-0.019	0.184	0.094
	A	2563.632	2221.687	…	2549.752	-0.044	0.184	0.048
	B	865.356	991.690	…	711.243	-0.024	0.153	0.206
	$\Delta$	7	10	…	8	-0.118	0.148	-0.057
HHT	IMF1 (inst. amp.)	83.463	49.466	…	5.519	0.200	0.049	-0.054
	IMF2 (inst. amp.)	237.164	257.480	…	281.664	0.149	0.178	0.152
	IMF3 (inst. amp.)	354.514	342.040	…	511.438	0.166	0.188	0.056
	IMF4 (inst. amp.)	442.211	444.300	…	443.832	0.171	0.174	0.046
	IMF1 (inst. freq.)	0.077	0.008	…	0.046	0.174	0.164	0.029
	IMF2 (inst. freq.)	0.010	0.023	…	0.023	0.200	0.049	-0.054
	IMF3 (inst. freq.)	0.005	0.011	…	0.004	0.174	0.164	0.029
	IMF4 (inst. freq.)	0.012	0.009	…	0.004	0.173	0.167	0.028
	IMF1 (power)	1286.833	465.873	…	12.844	0.128	0.185	0.139
	IMF2 (power)	4571.933	5682.005	…	5538.122	0.174	0.169	0.037
	IMF3 (power)	11795.380	9420.495	…	16683.370	0.124	0.183	0.147
	IMF4 (power)	11237.530	11440.120	…	8996.956	-0.051	0.113	0.006
ApEn	ApEn	0.056	0.057	…	0.055	0.028	0.103	0.096
	ApEn (IMF1)	0.087	0.110	…	0.024	0.110	0.118	0.091
	ApEn (IMF2)	0.102	0.110	…	0.103	0.109	0.155	0.168
	ApEn (IMF3)	0.039	0.048	…	0.039	0.022	0.113	0.167
	ApEn (IMF4)	0.033	0.032	…	0.027	0.075	0.110	0.191
Multi-SampEn	$\tau$=1	0.035	0.034	…	0.031	0.017	0.041	0.116
	$\tau$=5	0.122	0.122	…	0.109	0.008	0.024	0.160
	$\tau$=10	0.159	0.1558	…	0.157	0.016	0.079	0.157
	$\tau$=15	0.193	0.191	…	0.191	0.030	0.030	0.149
	$\tau$=20	0.216	0.215	…	0.223	0.142	0.083	-0.122
	$\tau$=25	0.237	0.233	…	0.248	0.126	0.029	-0.237
	$\tau$=30	0.259	0.247	…	0.261	0.160	-0.005	-0.052
WT	power (cA4)	1.85E+09	1.77E+09	…	1.54E+09	0.093	-0.091	0.142
	power (cD1)	1570.188	1549.465	…	1434.384	-0.066	0.000	0.039
	power (cD2)	6239.389	6607.213	…	5080.174	-0.099	-0.027	0.155
	power (cD3)	85840	95480	…	74108	-0.084	0.031	0.083
	power (cD4)	1263707	1381416	…	1090375	-0.080	0.011	0.117
WPT	power C(4, 0)	1.85E+09	1.77E+09	…	1.54E+09	0.107	0.099	-0.115
	power C(4, 1)	1263707	1381416	…	1090375	0.069	0.031	-0.236
	power C(4, 2)	82527.030	92128.26	…	71890.45	0.140	-0.017	-0.013
	power C(4, 3)	3813.841	3635.693	…	2255.475	0.050	-0.104	0.186
	power C(4, 4)	5628.319	5902.739	…	4668.485	0.199	0.047	-0.088
	power C(4, 5)	327.595	413.234	…	209.324	0.193	0.001	-0.081
	power C(4, 6)	265.937	137.527	…	94.082	0.192	0.012	-0.050
	power C(4, 7)	224.138	201.812	…	113.029	0.200	0.028	-0.086
	power C(4, 8)	754.763	502.210	…	377.318	-0.009	-0.076	0.004
	power C(4, 9)	141.005	172.995	…	100.334	-0.051	0.113	0.006
	power C(4, 10)	146.798	167.296	…	148.572	0.192	0.019	-0.080
	power C(4, 11)	142.433	127.221	…	117.892	0.015	0.006	0.039
	power C(4, 12)	151.214	153.350	…	117.711	-0.011	-0.001	0.025
	power C(4, 13)	143.662	145.584	…	151.531	-0.005	-0.006	-0.029
	power C(4, 14)	139.631	188.992	…	168.549	-0.013	0.000	0.054
	power C(4, 15)	189.569	185.274	…	268.314	0.012	-0.037	-0.011

---

 | Show Table

DownLoad: CSV

2.2.   Power spectral density (PSD)

Fourteen features were extracted from PSD ^[22]. PSD is a measure of the distribution of the power contents of signals in the frequency domain. Band energy ratio (BER) is derived from the frequency spectrum to show the energy distribution in a specific range. Since the human pulse energy distribution is concentrated below 10 Hz, BER in a particular band with respect to the total energy in the range of 0 to 10 Hz is calculated as follows.

where ${E}_{n}$ is the energy in the nth band, and ${E}_{t}$ is the total energy in the interval [0-10] Hz. In this study, signals in the interval [0-10] Hz were divided into 11 bands, as shown in Figure 3. Thus, 11 BERs were extracted, and they are corresponding to “BER1” to “BER11” in Table 1.

Figure 3.  Divide the interval [0, 10] Hz frequency into 11 bands.

DownLoad: Full-Size Img PowerPoint

The main frequency is the frequency with the largest energy. After the baseline drift is removed, the main frequency is the frequency of the pulse signals. The main frequency of the pulse signals and its PSD were extracted as 2 features, and they are corresponding to “Main freq.” and “Main freq. PSD” in Table 1, respectively. For example, in Figure 3, the main frequency and its PSD are 1.4648 and 2.62 × 10⁵, respectively.

Mean frequency, which is corresponding to “Mean freq.” in Table 1, is the sum of the product of the frequency and the corresponding PSD divided by the total PSD in the interval [0-10] Hz:

where, ${E}_{j}$ is the energy at frequency ${f}_{j}$, and l is the data length in the interval [0-10] Hz.

2.3.   Gamma density function

Five pulse features were extracted from a signal using a gamma density function. Shu and Sun ^[19] decomposed a time series pulse waveform into a forward wave and a backward wave using a gamma density function:

where A is the amplitude of the forward wave, B is the amplitude of the backward wave, α is the shape parameter of the pulse waveform, β is the rate parameter of the pulse waveform, and Δ is the phase shift or time delay between two waves. The best curve fitting was determined based on minimizing the difference between the recorded waveform and the gamma density function $F\left(t\right|\alpha, \beta, \Delta, A, B)$.

For example, in Figure 4, the red dots are an actual floating pulse, and the solid line is the corresponding fitted gamma density function $F\left(t\right|\alpha, \beta, \Delta, A, B)$ . In this study, particle swarm optimization (PSO) was used to acquire the optimized parameters. The average correlation of the fitted gamma density functions and the actual pulse waveforms was 0.99. After a fitted gamma density function $F\left(t\right|\alpha, \beta, \Delta, A, B)$ was found, parameters α, β, Δ, А, and В were taken as 5 pulse features, as shown in Table 1.

Figure 4.  A gamma density function (solid line) was used to fit a floating pulse (red dot line) using the PSO algorithm.

DownLoad: Full-Size Img PowerPoint

2.4.   Hilbert-huang transform (HHT)

Twelve features were extracted using HHT ^[23]. HHT is the result of the empirical mode decomposition (EMD) and Hilbert transform (HT) analysis. HHT was designed to analyze nonlinear and nonstationary data. HHT can provide a meaningful time-frequency-energy description of a time series. EMD decomposes signals into several intrinsic mode functions (IMFs) ^[24]. Figure 5 shows an example of using EMD to decompose a floating pulse into 7 IMFs and one residue. After obtaining all IMFs, average instantaneous amplitude ${\stackrel{-}{h}}_{n}$, average instantaneous frequency ${\stackrel{-}{\omega }}_{n}$, and power ${P}_{n}$ can be calculated as follows ^[8].

Figure 5.  A floating pulse is decomposed into 7 IMFs and one residue.

DownLoad: Full-Size Img PowerPoint

where l is the data length, N is the number of IMFs, ${a}_{n}\left(t\right)$ is the instantaneous amplitude, and ${f}_{n}\left(t\right)$ is the instantaneous frequency. In this study, the first 4 IMFs were extracted. For each IMF, the average instantaneous amplitude, average instantaneous frequency, and power were calculated, and they are corresponding to “IMF1 (inst. amp.)” to “IM4 (inst. amp.)”, “IMF1 (inst. freq.)” to “IM4 (inst. freq.)”, and “IMF1 (power)” to “IM4 (power)” in Table 1, respectively. Therefore, 12 features were obtained using HHT.

2.5.   Approximate entropy (ApEn)

Five features were derived using ApEn ^[20]. ApEn measures the complexity of a time series. Higher entropy values indicate a system exhibiting a greater degree of complex dynamics. Prior research found that the ApEn value of a normal health condition was higher than that of an abnormal health condition ^[25,26]. In this study, the ApEn of the original pulse signals and the first four IMFs were calculated to obtain 5 ApEn features, and they are corresponding to “ApEn” and “ApEn (IMF1)” to “ApEn (IMF4)” in Table 1, respectively.

2.6.   Multiscale sample entropy (Multi-SampEn)

Sample entropy (SampEn) is a modification of ApEn ^[27] and has the advantage of being less dependent on the time series length. However, a single-scale SampEn might not be sufficient to separate healthy and pathologic individuals. Costa et al. ^[21] calculated entropy using multiple time scales, called Multi-SampEn, to robustly separate the time series of healthy subjects and patients with severe heart disease. Compared with the traditional entropy method, Multi-SampEn can yield a higher complexity and more meaningful measure.

For the original signals with length l, a consecutive coarse-grained time series which represented the original time series on different time scales is constructed by

where τ is a scale factor. The original time series is divided into non-overlapping windows of length τ, and the data points inside each window are averaged. By acquiring the SampEn of the new coarse-grained time series, the Multi-SampEn of the signals at scale τ can be obtained. In this study, Multi-SampEn with 7 time scales ($\tau$ = 1, 5, 10, 15, 20, 25, and 30) were calculated for each data set. Therefore, 7 features were derived using Multi-SampEn, as shown in Table 1.

2.7.   Wavelet transform (WT)

Five features were extracted using WT ^[28]. WT has been widely used in various fields such as signal processing and data compression. Using WT, a signal is split into two coefficient groups, a low-frequency approximation group (cA_x) and a high-frequency detail group (cD_x), where x is the level number. The decomposition process is iterated, with the successive approximation groups being decomposed in turn until a predefined level is reached.

In this study, a Daubechies wavelet of order 2 (db2) was used as the filter, and the decomposition process was carried out until level 4 ^[6,29]. The result included 4 detail coefficients (cD₁-cD₄) and 1 approximation coefficient (cA₄). After calculating the wavelet coefficients, the wavelet powers were calculated as follows:

where LcD_i is the data length of cD_i, LcA₄ is the data length of cA₄, k is the sample data, and i is the level number. The wavelet powers of the approximation coefficient (cA₄) and 4 detail coefficients (cD₁-cD₄) are corresponding to “power (cA₄)” and “power (cD₁)” to “power (cD₄)” in Table 1, respectively.

2.8.   Wavelet packet transform (WPT)

Sixteen features were extracted using WPT ^[6,30]. Contrary to WT, WPT splits both the detail coefficients and the approximation coefficients. Thus, WPT provids a more sophisticated analysis of the signals. In this study, a Daubechies wavelet of order 2 (db2) was used as the filter, and the decomposition process was carried out until level 4. Therefore, the signals were decomposed into 16 sub-bands. The wavelet packet power of each sub-band was calculated as follows.

where Lc(4, i) is the length of the i-th sub-band at level 4, and k is the sample data. The wavelet packet power of the 16 sub-bands are corresponding to “power C(4, 0)” to “power C(4, 15)” in Table 1, respectively.

2.9.   High-dimensional feature vector

In this study, 500 pulse data sets were collected. For each original data set, a 71-dimensional feature vector was created using the above-mentioned feature extraction methods. Table 1 shows the final pulse features and the corresponding values of the original 500 dataset and the loadings of each feature element in the first 3 Principal Components (PCs).

2.10.   Dimension reduction

Although a high dimensional feature vector contains more pulse information, it might also contain unwanted noises or redundant information, which might affect classification accuracy and increase computational complexity. To improve classification performance and understand the underlying physical meaning or implications hidden in pulse classification, PCA was applied to extract the most representative information.

Figure 6 shows that the first 3PCs covers about 52% of the total accumulative variances (PC1 gives 32.7478%, PC2 gives 12.9812%, and PC3 gives 6.3234%). That means that the first 3 PCs explain about 52% of the variability in the original 71 variables. If the complexity of the data is reduced by using the first 3 PCs, 48% of data information will be lost. The last three columns in Table 1 show the coefficients (or loadings) of each feature element in the first 3 PCs. The larger the absolute value of the loading, the more important the corresponding feature is in calculating the component. The loadings in the first 3 PCs show that the first component (PC1) has large associations with HHT features and WPT features, which indicates that PC1 measures the energy and local instantaneous characteristics of the pulses. The second component (PC2) has large associations with PSD features around the main frequency, which indicates that PC2 measures the characteristics of the main frequency. The third component (PC3) has large associations with mean frequency, waveform features and entropy features, which indicates that PC3 measures the complexity of the waveforms.

Figure 6.  PCA scree plot.

DownLoad: Full-Size Img PowerPoint

3.   Artificial neural network (ANN)

After extracting all features, a basic ANN was used to classify the 10 pulse types. Although the first 3 PCs covered more than 50% of the total variances, because some pulse information was subtle, five ANN models using the first 8, 13, 21, 27, and 71 PCs, which accounted for 70%, 80%, 90%, 95%, and 100% of the total variances, were trained, respectively, using the scaled conjugate gradient back-propagation algorithm to update the weights and bias values of the ANN.

In this study, the number of nodes in the ANN input layer was based on the feature dimension, and the numbers of nodes in the hidden layer and the output layer were 10 each. The maximum number of epochs was 1000, the learning rate was 0.01, and the minimum gradient was e^-6. The sigmoid activation function was used for the hidden layer, and the softmax function was used for the output layer. The training was terminated when the average error was less than 10^-4. This study randomly selected 70% of the data for training, 15% for validation, and 15% for test.

Table 2 shows that Model 4, using the first 27 PCs, has the highest overall average classification accuracy of 98.2%. Figure 7 shows the distances of the clusters of the 10 pulse types in the 27-dimensional PC space, based on the Ward’s hierarchical clustering method. It reveals that pulses 7 and 9, and pulses 3 and 4 have the shortest inter-cluster distances. It might affect their discriminability. However, since ANN is a nonlinear regression method, the distance between two pulses or two clusters might not be the only factor in determining pulse discriminability.

Table 2.  Average classification accuracy of each model.

ANN model	Model 1	Model 2	Model 3	Model 4	Model 5
No. of PCs	8	13	21	27	71
Accumulative variances (%)	70	80	90	95	100
Classification accuracy (%)	92.8	94.6	97.4	98.2	96.6

---

 | Show Table

DownLoad: CSV

Figure 7.  Clustering of the 10 pulses in the 27-dimensional PC  .

DownLoad: Full-Size Img PowerPoint

Table 3 shows the confusion matrix of Model 4. Each row corresponds to the predicted pulse (Output class) and each column corresponds to the true pulse (Target class). The diagonal cells correspond to the results that are correctly classified. The far right column of the table is the precision rates, which represent the percentages of all the predicted pulses in a row that are correctly classified. The row at the bottom of the table is the recall rates, which represent the percentages of all the true pulses in a column that are correctly predicted. The cell in the bottom right of the table is the overall accuracy. The result shows that the precision and recall rates differ for each pulse type.

Table 3.  Confusion matrix of Model 4.

Pulses		Target class
		1	2	3	4	5	6	7	8	9	10	Precision
	1	48	0	0	0	0	0	0	0	0	0	100%
	2	0	50	0	0	1	0	0	0	0	0	98%
	3	0	0	49	0	0	0	0	0	0	0	100%
	4	0	0	1	50	0	0	0	0	0	0	98%
	5	0	0	0	0	49	0	0	0	0	0	100%
	6	0	0	0	0	0	50	0	0	0	0	100%
	7	1	0	0	0	0	0	47	0	1	0	95.9%
	8	0	0	0	0	0	0	0	50	0	1	98%
	9	0	0	0	0	0	0	3	0	49	0	94.2%
	10	1	0	0	0	0	0	0	0	0	49	98%
	Recall	96%	100%	98%	100%	98%	100%	94%	100%	98%	98%	98.2%

---

 | Show Table

DownLoad: CSV

The F-measure ( = $2\bullet \frac{precision\bullet recall}{precision+recall})$, a balance between precision and recall, was used to represent the classification rate for each pulse type. Table 4 shows that pulses 7 ((long+ replete) pulse) and 9 ((normal + replete) pulse) have the lowest classification rates in average. Most models confuse pulses 7 and 9. It might be because pulses 7 and 9 both contain replete pulse features and have the shortest average inter-cluster distance in the 27-dimensional PC space. It indicates that replete pulse dominates long pulse and normal pulse. It also indicates that the strength features of a pulse dominates the shape features of a pulse.

Table 4.  F-measures of all pulse types in the five models.

	Pulse Type	1	2	3	4	5	6	7	8	9	10
	Model 1	98	90.7	91.9	91.1	93.1	98.1	88.9	97.2	83.7	95.1
	Model 2	97.2	92.9	98	97.2	96	97.1	85.2	99	84.4	100
	Model 3	100	97.1	100	98	99	100	92.5	99	89.8	99
	Model 4	98	99	99	99	99	100	95	99	96.1	98
	Model 5	92.4	93.9	99	97	98.1	97	99	99	96	96
	Avg.	97.1	94.7	97.6	96.5	97.0	98.4	92.1	98.6	90.0	97.6

---

 | Show Table

DownLoad: CSV

However, on the other hand, pulses 2 (rapid pulse), 6 ((rapid + stirred) pulse), 8 ((rapid + sunken) pulse), and 10 (sunken + rapid + stirred) contain rapid pulse and are neighboring each other in the 27-dimensional PC space, their discrimination rates are very high, especially for pulse 6. It indicates that stirred pulse and sunken pulse are very distinct, and they dominate rapid pulse. It also indicate that the shape features and the depth features of a pulse are very distinct, and they dominate the rate features.

4.   Discussions

Table 5 shows a comparison between the proposed method and the prior TCM pulse classification methods. Most prior pulse classification methods used low-dimensional feature vectors. Therefore, some important but subtle pulse information might be ignored. The proposed high-dimensional feature classification method could retain more pulse information and thus achieved a higher classification accuracy.

Table 5.  Comparison between the proposed method and prior methods.

Research	Data	Feature Dimension	Classes	Classifier	Accuracy
Wang and Cheng [16]	407	13	5 pulse types	Bayesian networks	84.20%
Xu et al. [31]	900	n/a	6 pulse types	wavelet network	83%
Xu et al. [17]	n/a	4	7 pulse types	Lempel-Ziv complexity analysis	97.10%
Shu and Sun [19]	n/a	4	13 pulse types	n/a	n/a
Xu et al. [18]	320	17	16 pulse types	FNN	90.25%
The proposed method	500	71 (w/o PCA)	10 pulse types	ANN	96.60%
	500	27 (w. PCA)	10 pulse types	ANN	98.20%

---

 | Show Table

DownLoad: CSV

There are some limitations in this study. First, TCM doctors generally use their index, middle, and ring fingers to take patients’ pulses at three different locations (Cun, Guan, and Chi) on the wrists, using three different pressures (light touch, moderate touch, and heavy touch). However, in this study, pulses at the Guan location on the left wrist were taken only. In the future, pulses on both wrists with different pressing forces will be taken to include more pulse types in pulse classification.

5.   Conclusions

Most prior pulse classification research used low-dimensional feature vectors to classify few known diseases or pule types. However, a feature generation method might only be effective for extracting certain features but not for others. In order to increase pulse classification accuracy, this research proposed a novel high-dimensional feature extraction method to extract as much important and subtle pulse information as possible from the time, spatial, and frequency domains. Eight different feature generation methods were applied to construct a 71-dimensional feature vector. Extracting high-dimensional features can also help to understand the underlying physical meaning or implications hidden in pulse discrimination.

ANN results show that PCA accounted for 95% of the total variances achieved the highest accuracy of 98.2% in pulse classification. The results also showed that pulse energy, local instantaneous characteristics, main frequency, and waveform complexity are the major factors determining pulse discriminability. Some pulse features also outperform or dominate other pulse features. For example, the strength features dominate the shape features, and the shape features and the depth features dominate the rate features.

Acknowledgments

The authors would like to thank the Ministry of Science and Technology of Taiwan for financially supporting this research under contract MOST 108-2221-E-002-161-MY2.

Conflict of interest

The authors declare no conflict of interest.