Time-Shift Multiscale Entropy Analysis of Physiological Signals

Abstract

Measures of predictability in physiological signals using entropy measures have been widely applied in many areas of research. Multiscale entropy expresses different levels of either approximate entropy or sample entropy by means of multiple factors for generating multiple time series, enabling the capture of more useful information than using a scalar value produced by the two entropy methods. This paper presents the use of different time shifts on various intervals of time series to discover different entropy patterns of the time series. Examples and experimental results using white noise, $1/f$ noise, photoplethysmography, and electromyography signals suggest the validity and better performance of the proposed time-shift multiscale entropy analysis of physiological signals than the multiscale entropy.

1. Introduction

The notion of the approximate entropy (ApEn) [1], which quantifies irregularity or predictability in scalar time series has been increasingly applied in various scientific domains of signal processing. The modified version of ApEn known as the sample entropy (SampEn) [2] was introduced to remove self-matching in ApEn, and reported to produce results better than the use of ApEn in several studies of time-series analysis [3,4,5].

Consider a scalar time series X of length N taken at regular intervals: $X=(x_1,x_2,\cdots ,x_N)$, and a given embedding dimension m, a set of newly reconstructed time series from X, denoted as Y, can be established as $Y=(y_1,y_2,\cdots ,y_{N-m+1})$, where $y_i=(x_i,x_{i+1},\cdots ,x_{i+m-1})$, $i=1,2,\cdots ,N-m+1$. Given a positive tolerance value r, the probability of vector $y_i$ being similar to vector $y_j$ is computed as

$$C_i^m\left(r\right)=\frac{1}{N-m+1}\sum _{j=1}^{N-m+1}\theta \left[d(y_i,y_j)\right],$$

(1)

where $\theta \left(d(y_i,y_j)\right)$ is the step function defined as

$$\theta \left[d(y_i,y_j)\right]=\left\{\begin{array}{cc}1\hfill & \hfill :d(y_i,y_j)\le r,\\ 0\hfill & \hfill :d(y_i,y_j)>r.\end{array}\right.$$

(2)

The distance between the two vectors can be obtained by

$$d(y_i,y_j)=\underset{k}{\max }\left(\right|x_{i+k-1}-x_{j+k-1}\left|\right),k=1,2,\cdots ,m.$$

(3)

The probabilities of all vectors being similar to one another are computed as

$$C^m\left(r\right)=\frac{1}{N-m+1}\sum _{i=1}^{N-m+1}\log \left[C_i^m\left(r\right)\right].$$

(4)

Approximate entropy, denoted as ApEn, is defined as

$$\mathrm{ApEn}=C^m\left(r\right)-C^{m+1}\left(r\right).$$

(5)

For the calculation of ApEn, m = 2 or 3, and r = 0.1 to 0.25 $\times \sigma$, where $\sigma$ is standard deviation of the time series, were typically suggested [6]. Let $B_i^m\left(r\right)$ be defined as

$$B_i^m\left(r\right)=\frac{1}{N-m-1}\sum _{j=1}^{N-m}\theta \left[d(y_i,y_j)\right],i\ne j.$$

(6)

Thus, $B^m\left(r\right)$ is given by

$$B^m\left(r\right)=\frac{1}{N-m}\sum _{i=1}^{N-m}{B}_i^m\left(r\right).$$

(7)

The formulation of SampEn is expressed as

$$\mathrm{SampEn}=-\log \left[\frac{B^{m+1}\left(r\right)}{B^m\left(r\right)}\right].$$

(8)

In particular, the multiscale entropy (MSE) [7] was developed to measure entropy such as SampEn at different scales by averaging non-overlapping time points of the original time series in order to better reveal patterns of predictability or regularity in the time series. The use of MSE analysis of time series is able to explain the inconsistency encountered with single-scale analysis in the increase and decrease in entropy values of certain physiological signals, and have been found useful in many applications, most recently such as [8,9,10,11]. The MSE works by applying a “coarse-graining” process to the original time series to generate several time series of different scales, and then computing either SampEn or ApEn for all coarse-grained time series, which are plotted as a function of the scale factor. Consider the time series X of length N taken at regular intervals: $X=(x_1,x_2,\cdots ,x_N)$. For a scale factor $\tau$, a new time series $X^{\tau }$ is created by the MSE as follows [7]:

$$X_j^{\tau }=\frac{1}{\tau }\sum _{i=(j-1)\tau +1}^{j\tau }{x}_i,1\le j\le N/\tau .$$

(9)

2. Time-Shift Multiscale Entropy

The proposed time-shift multiscale entropy, denoted as TSME, is motivated by a popular approach for computing the fractal dimension of irregular time series, known as the Higuchi’s fractal dimension (HFD) [12]. The HFD computes the “mean length" of the curve of a time series by constructing a set of new time series that has the property of a fractal curve over all time scales as each time series can be considered a reduced scale form of the whole. A set of new time series constructed from the original time series by the HFD, which are utilized in this study, are based on the consideration of the phase distribution. This phase distribution can reveal strong effects of the irregularity of time series [13]. In fact, it has been reported that the HFD is a stable numerical approach for time series analysis, including stationary, nonstationary, deterministic, and stochastic signals [14,15,16,17]. Therefore, the time-shift multiscale entropy can be expected to be applied to these types of signals for studying irregularity in time series.

The new time series generated by the HFD are constructed as follows by once again considering the time series X of length N : $X=(x_1,x_2,\cdots ,x_N)$. Let $\beta$ and k be positive integers, where $\beta =1,2,\cdots ,k$, then k new time series can be generated using the following equation [12]:

$$X_k^{\beta }=(x_{\beta },x_{\beta +k},x_{\beta +2k},\cdots ,x_{\beta +\lfloor \frac{N-\beta }{k}\rfloor k}),$$

(10)

where $\lfloor \frac{N-\beta }{k}\rfloor$ is the “floor” function that rounds $\frac{N-\beta }{k}$ to the largest integer not exceeding $\frac{N-\beta }{k}$.

From Equation (10), $\beta$ and k indicate the initial time point and time interval, respectively, that is, for a given time interval k, k new time series are constructed using k time shifts. For example [18], let k = 3, and N = 100, three new time series are generated as: $X_{k=3}^{\beta =1}:(x_1,x_4,x_7,\cdots ,x_{100})$, $X_{k=3}^{\beta =2}:(x_2,x_5,x_8,\cdots ,x_{98})$, and $X_{k=3}^{\beta =3}:(x_3,x_6,x_9,\cdots ,x_{99})$.

The proposed TSME method works by constructing k time-shift series for a given time interval k, then computing either SampEn or ApEn for all time-shift time series, denoted as $TSM{E}_k^{\beta }$, $\beta$ = 1, ..., k. The TSME for each k, denoted as $TSM{E}_k$, k = 1, ..., $k_{max}$, is defined as the average value of all $TSM{E}_k^{\beta }$, that is,

$$TSM{E}_k=\frac{1}{k}\sum _{\beta =1}^{k}TSM{E}_k^{\beta }.$$

(11)

The procedure for computing TSME is described as follows:

• Given a scalar time series X, dimension m, tolerance r, and $k_{max}$.
• Set k = 1.
• Using Equation (10) to construct k time-shift time series from X : $X_k^{\beta }$, $\beta =1,\cdots ,k$.
• For each $X_k^{\beta }$, $\beta =1,\cdots ,k$, compute $TSM{E}_k^{\beta }$.
• Compute $TSM{E}_k$ using Equation (11).
• Set $k=k+1$.
• Repeat steps 3–6 until $k=k_{max}$.

3. Results and Discussion

3.1. Analysis of Signals with Known Properties

In signal processing, white noise and $1/f$ (pink) noise are signals with known statistical properties, which are widely found as good approximation of many real-world data and generate mathematically tractable models. Time series of white noise are sequences of serially uncorrelated random variables with zero mean and finite variance, having equal intensity at different frequencies and constant power spectral density [19]. Time series of $1/f$ noise is a signal with a frequency spectrum such that the power spectral density is inversely proportional to the frequency of the signal; it is an intermediate between the white noise and random walk noise with no correlation between increments [20]. Signals of white noise (mean = 0 and variance = 1), and $1/f$ noise, both consist of 10,000 samples, were used in this analysis. Another type of signal of known property is the Lorenz system [21], consisting of x (convection velocity), y (temperature difference), and z (temperature gradient) components, which are well known to be chaotic as shown in Figure 1. These three signals, each having the length of 4000 samples, behave irregularly with time: the x and y components fluctuate around positive and negative values, and the z component oscillates around the range from about 10 to 40.

To carry out the TSME analysis, we set m = 2 and r = $0.15\times \sigma$, where $\sigma$ is the standard deviation of the original time series, to compute ApEn and SampEn. With $k_{max}$ = 10, the ApEn-based $TSM{E}_k^{\beta }$ tends to slightly decrease with increasing k from 2.4 for k = 1 to 1.4 for k = 10, while the SampEn-based $TSM{E}_k^{\beta }$ remains fairly constant around the value of 2.4. Figure 2 shows the $TSM{E}_k$ of the white noise and $1/f$ noise, respectively, where the SampEn-based $TSM{E}_k$ for both signals are fairly constant. Both SampEn-based $TSM{E}_k^{\beta }$ and SampEn-based $TSM{E}_k$ meet the expectation for the analysis of randomly generated time series. The trend in ApEn-based $TSM{E}_k^{\beta }$ and ApEn-based $TSM{E}_k$ should be due to the bias of self-matching in the computation of ApEn [2].

Figure 3 shows that $TSM{E}_k^{\beta }$ using either ApEn or SampEn can distinguish white noise and $1/f$ noise from the chaotic time series of the Lorenz system. Similarly, Figure 4 shows that $TSM{E}_k$ using either ApEn or SampEn can distinguish white noise and $1/f$ noise from the chaotic time series. Particularly, $TSM{E}_k$ can also separate white noise from $1/f$ noise. $TSM{E}_k^{\beta }$ and $TSM{E}_k$ values of white noise and $1/f$ noise are higher than those of the three chaotic time series as shown in Figure 3 and Figure 4, respectively, which suggest the reliability of the TSME analysis, as chaotic signals are deterministic and therefore more predictable than noise data.

3.2. Analysis of PPG and EMG Signals

To illustrate the performance of the proposed TSME and comparison with the MSE, experiments using the photoplethysmography (PPG) and electromyography (EMG) signals are presented. The pulses of the index fingers of the left hands of 43 elderly participants and the middle-age caregiver were synchronously measured with a PPG sensor, and studied in [22] for automated assessment of therapeutic communication for cognitive stimulation for people with cognitive decline. The EMG signals were obtained from the Physical Action Data Set [23], where the channel measured on the right bicep of the participants with the Delsys EMG wireless apparatus was used in this study, for the classification of normal and aggressive human physical actions. Figure 5 shows the first 5000 samples of synchronized PPG signals of two elderly participants and the caregiver. Figure 6 shows the first 5000 samples of the EMG signals of six normal actions: bowing, clapping, handshaking, hugging, jumping, and running. Figure 7 shows the first 5000 samples of the EMG signals of six aggressive actions: elbowing, front kicking, hammering, kneeing, pulling, and punching.

For the MSE analysis, we set $\tau$ = 20 to compute $X_j^{\tau }$ expressed in Equation (9) as the MSE feature of 20 scale factors. For the TSME, we set $k_{max}$ = 20 to compute $TSM{E}_k$ expression in Equation (11) as the TSME feature of 20 phase shifts. Only SampEn was used in this analysis for computing both MSE and TSME, with m = 2 and r = $0.15\times \sigma$, where again $\sigma$ is the standard deviation of the original time series. Both PPG and EMG datasets were used for pattern classification.

The first 5000 samples of the synchronized, de-trended PPG signals of 43 elderly participants and the middle-age caregiver were used to classify the PPG signals of the elderly from the caregiver. MSE and TSME features using SampEn were extracted from these signals, and classified using the linear discriminant analysis (LDA) [24,25]. The leave-one-out (LOO) cross validation was applied to test the accuracy of LDA-based classification of the two types of features. Figure 8 shows the areas under the receiver operating characteristic (ROC) curves [26], denoted as AUC (area under curve), obtained from the MSE and TSME features extracted from the PPG signals, where the AUC of the MSE = 0.74, and the AUC of the TSME = 0.84 (the higher the value of the AUC, the better the performance).

Table 1. Sensitivity (SEN), specificity (SPE), and leave-one-out (LOO) cross validation obtained from linear discriminant analysis of synchronized photoplethysmography (PPG) signals of elderly participants and caregiver, using multiscale entropy (MSE) and time-shift multiscale entropy (TSME).

Feature	SEN (%)	SPE (%)	LOO (%)
MSE	51.16	93.02	44.19
TSME	62.50	97.50	62.50

shows the sensitivity (the rate of elderly features that are correctly identified), specificity (the rate of caregiver features that are correctly identified), and accuracy rates obtained from the LDA using the MSE and TSME features. The TSME feature yielded better results than the MSE feature in terms of accuracy and receiver operating characteristics.

The first 5000 samples of the EMG right bicep signals of the four subjects in the Physical Action Data Set [23] were also used for the differentiating normal from aggressive actions, using LDA analysis with MSE and TSME features. The experiments were carried out on each subject performing ten normal and ten aggressive physical actions that represented human activities. The ten normal actions are: (1) bowing; (2) clapping; (3) handshaking; (4) hugging; (5) jumping; (6) running; (7) seating; (8) standing; (9) walking; and (10) waving. The ten aggressive actions include: (1) elbowing; (2) front kicking; (3) hammering; (4) heading; (5) kneeing; (6) pulling; (7) punching; (8) pushing; (9) side kicking; and (10) slapping. Figure 9 shows the ROC areas under curves obtained from the MSE and TSME features, where the AUC of the MSE = 0.79, and the AUC of the TSME = 0.84, showing the better performance of the TSME.

Table 2. Sensitivity (SEN), specificity (SPE), and leave-one-out (LOO) cross validation obtained from linear discriminant analysis of electromyography (EMG) signals of normal and aggressive actions, using MSE and TSME.

Feature	SEN (%)	SPE (%)	LOO (%)
MSE	67.50	77.50	50.00
TSME	82.50	72.50	61.25

shows the sensitivity (the rate of aggressive action features that are correctly identified), specificity (the rate of normal action features that are correctly identified), and accuracy rates obtained from the LDA using the MSE and TSME features. Once again, the TSME feature yielded better results than the MSE feature in general. The specificity (77.50%) obtained from the MSE is 5% higher than the specificity (72.50%) obtained from the TSME, while the sensitivity (67.50%) obtained from the MSE is 15% lower than the sensitivity (82.50%) obtained from the TSME, and the accuracy (50%) of the MSE is 11.25% lower than the accuracy (61.25%) of the TSME.

4. Conclusions

A method for computing SampEn or ApEn on multiple time shifts has been presented and discussed. The use of SampEn is more theoretically sound for computing TSME, as it is adopted for computing MSE. The construction of the new time series for computing the entropy profiles introduced in this study was known to be able to provide stable time scale and indices corresponding to the characteristics of irregular time series, including short time series, by taking into account self-similarity across the characteristic time scale [12]. While MSE uses the averaging of time series on several interval scales, the proposed TSME applies time shifting in time series that is based on the calculation of the “mean length” of the curve of a time series implemented in the Higuchi’s fractal dimension.

The examples using noise and chaotic time series suggest the validity of the TSME, and classification results using the physiological data illustrate the better performance of the TSME over the MSE. The computational time required for computing the TSME is higher than for computing the MSE, particularly when $k_{max}$ is large. Being similar to the specification of the $\tau$ parameter used in the MSE, an optimal selection for the $k_{max}$ used in the TSME is still an ad hoc choice reported in many studies [17], which needs further investigation. Future development of the TSME for multivariate multiscale entropy is worth pursuing, as it has been formulated based on the concept of the MSE [27,28]. An extension of the differential Shannon entropy rate of time series using kernel density estimators for selecting the order and bandwidth parameters [29] may have implications for improving the performance of the TSME.

As another issue for the implementation of the TSME regarding a good selection of the tolerance parameter r in the computation of ApEn or SampEn, the quadratic sample entropy (QSE) [30] was introduced as a complementary stochastic approach to measuring complexity by adding $\log \left(2r\right)$ to the deterministic approach ApEn or SampEn so that r can be optimally varied for individual time series. In studying neonatal heart rate data with transient decelerations, it was shown that, for the stochastic approach, the entropy tends to be converging as r approaches 0, while for the deterministic approach (SampEn), the entropy is diverging [30]. The QSE also found potential in the application to the analysis of very short physiological time series for the automated detection of atrial fibrillation in implanted ventricular devices [31]. Therefore, the utilization of the QSE in the computation of the TSME for the classification of physiological signals is worth investigating in future study.

The Matlab code for computing the TSME is available at the author’s personal homepage: https://sites.google.com/site/professortuanpham/codes.