Sage Journals: Discover world-class research

Abstract

In contrast with the large diesel engine fuel system operational condition complexity, the chaotic fractal characteristics in the combustion process of the diesel engine fuel system are systematically analyzed. Then, based on chaotic fractal method, an effective fault diagnosis method is further proposed. And under the two kinds of typical failure modes, the two important characteristic parameters variation trend for chaotic system, namely, correlation dimensions and the calculation methods of the maximum Lyapunov exponent are discussed. On the process of chaotic parameters calculation the wavelet noise reduction method is used to handle the original signal to conquer the linear scaling region properly. The results show that the maximum Lyapunov exponent method is a more effective method on the diesel engine fuel system fault diagnosis than the correlation dimension method, which can be further used in health condition assessment for this type fuel engine.

Keywords

the largest Lyapunov exponent correlation dimension fault vibration signal diesel engine fuel system health condition assessment

Introduction

The engine combustion process is a complex process. The process can be defined as several development stages, and the degree of combustion of each stage is different. The prior stage is the preparation and foundation of the follow-up stage, which is also influenced by the development of the later stage. From the acoustic point of view, the vibration characteristics of each stage corresponded to the different stages, and under the influence of prior and later vibration signal the components of vibration signal are increasingly complex.¹ Furthermore, within the filtering effect on delivery channel of engine vibration signal and the dissipation and reverberation effect on delivery channel of elastic wave, the vibration ignition transient reflected form of explosive power is quite complex, so the traditional method has the limitation on the analysis of the engine vibration signal.

The engine combustion process has the character of hierarchy and self-similarity. As an operational cycle of engine, the combustion process of the fuel in the cylinder is in an instant conversion process from chemical energy into heat energy. This is the first level of awareness in the combustion process. The traditional method using vibration signal to study the process of the engine combustion is based on this level. If we see the combustion process as a whole, the whole combustion process can be seen as the composition of several secondary combustion processes occurring at different times and different regions. This is the second level of awareness in the combustion process. Based on this point, the typical four-phase method is introduced to analyze the combustion process. Moreover, each of the inferior combustion process can also be seen as the combined action of several subordinate combustion processes with similar property and shape. This is the third level of awareness in the combustion process. And so forth, the combustion process can be distinguished with more inferior scales. From this point of view, the combustion process is similar with the self-similarity feature of fractal chaos. Therefore, the self-similarity fractal feature of the combustion process which is reflected by the acoustic signal can be used in fault diagnosis.²

As we know, in order to realize the fault diagnosis on complicated system, some critical variable should be separately introduced to be monitored. The variation tendency of the critical variable can be used to analyze the complicated system. But in the real condition, the critical variables in the complicated system are influenced by each other. It is a hard work for us to extract the critical variables separately to implement on monitoring. The fuel engine is such a complicated system. As we know, the effectiveness of the combustion is influenced by many factors such as air admission, oil admission, reduction ratio, injection pressure, fuel supply angle, etc. Even now, no independent variables for the process of combustion are admitted officially. There are two reasons. On the one hand, it is difficult to extract the critical variable from the complicated system. And on the other hand, the accuracy is queried when extracting the critical variable to ignore the influence of the other critical variable. So the traditional linear method to fault diagnosis on the engine is difficult to realize. Then, we can use the fractal method to analyze the vibration signal, which is the combined action of all critical variables. If the time series of the vibration signal is provided with the chaotic character, the chaotic method can be used to analyze the time series.

Due to the nonlinear character of the diesel engine fuel system, especially when estimating the failure of the system, two classes of chaotic methods have been used to diagnose. One method is using Duffing oscillator to weak signal detection. The method based on Duffing oscillator has been used to diagnose the gearbox fault.³ The other method is using the chaotic character extraction of fault systems such as Lyapunov exponent and correlation dimension (CD).^4,5 There are various methods used to calculate the character of chaos system. Tan et al.⁶ show that there are fractal and chaos features in the denoising vehicle engine vibration system and the CD can be used as the efficient diagnosis of potential vibration fault from the vehicle diesel engine. Ma et al.⁷ analyzed the chaotic vibration signal among aero engine; the vibration signals are separated by means of fast independent component analysis (ICA) blind source separation, and the chaos is detected via a power spectrum and Lyapunov exponent based on the characteristics of chaos. Wu et al.⁸ resolved engine valve fault diagnosis problems, and meanwhile three chaotic features such as the CDs, maximum Lyapunov exponent, Kolmogorov entropy, and their statistical features were obtained as the fault features. The results indicate that the classification performance of the single chaotic feature is poorer, and with statistical features in tandem with chaotic features as the feature vector of fault signal, the performance of fault diagnosis is better. Yilmaz and Guler⁹ used both the largest Lyapunov exponent (LLE) obtained by Wolf method and the CD obtained by Grassberger–Procaccia (G–P) method to analyze the Doppler signals. Cai et al.¹⁰ proposed an extending method based on Lyapunov exponent forecasting model which is effective in fault diagnosis of rotating machinery, and the results demonstrate that boundary effects can be controlled effectively and faults can be recognized exactly. Wang et al.¹¹ separately calculated the nonlinear multiparameters of the CD, Lyapunov exponent, complexity, and approximate entropy which are applied in the research of the fault diagnosis of rotating machines. Studies show that the fault type is different, nonlinearity is significantly different, which verifies that the nonlinear feature quantities are effective parameters for fault information, and thus a more effective way is provided for studying the fault diagnosis of complexity rotating machinery.

In this paper, we focus on the F3L912 fuel engine, and the vibration acceleration signal of cylinder cover is monitored. The specified chaotic fractal method is introduced, and the effectiveness of the method is testified by the simulation data. Then, the original signal of the engine is tested in the method. Because of the impact of noise, the linear scaleless range can be found clearly. The wavelet denoise method is introduced to decrease the impact of noise. Two classes of fault mode are introduced to diagnose by calculating the CD and the maximum Lyapunov exponent. The result shows that the CD method cannot be used to fault diagnose, while the maximum Lyapunov exponent method can diagnose the engine effectively.

The sections of the paper are organized as follows. “Introduction” section describes that the vibration signal has the chaotic character from qualitative analysis of the combustion process. This is the foundation of our study. In “Phase space reconstruction on chaotic time series” section, in order to study the chaotic character of the monitored time series, we need to use phase space reconstruction method to reveal the dynamic information. The proper embedded dimension and time delay take direct effect on the quality of phase space reconstruction. In “The experiment data collection and pretreament” section, we introduced a proper noise reduction method. The feature extraction is made on the reconstructed time series and the results are given in “Analyses of the diesel engine fuel system fault vibration signal” section. The final section is the conclusion.

Phase space reconstruction on chaotic time series

As the chaotic time series conceal the dynamic information of chaotic system and the phase space reconstruction method can reveal the dynamic information, we can reconstruct the chaotic oscillator which depicts the trajectory of the chaotic system.¹² That is to say, proper reconstructed chaotic oscillator can regain the chaotic system as accurate as possible. From this point, if the engine is a specific chaotic system, the vibration signal monitored by the engine is also the signal with the characteristic of chaos, and the time series of signal includes all information in the engine. Then, we can reconstruct a new time series including all important information, and the analysis of the reconstructed time series can be used to realize the fault diagnosis.

Let {x(t)|t = 1, 2, … , N} be the observed time series. According to Packard et al.¹³ and Takens et al.,¹⁴ based on the selected embedded dimension m and the time delay τ, we can reconstruct the phase space as follows

X_{i} = {x (i), x (i + τ), \dots, x (i + (m - 1) τ)}

(1)

where i = 1, 2, … , M and M = N – (m – 1) τ.

According to Taken’s theorem,¹⁴ the embedded dimension m and the time delay τ are the basic parameters of the constructed phase space. There are two different points on the two parameters calculation. One point is that the two parameters are irrelevant to each other and the value of the parameters can be calculated separately. The time delay τ can be calculated by mutual information method or autocorrelation method, and the embedded dimension m can be calculated by the false nearest neighbors method or G–P method. The other point is that the two parameters are relevant to each other; the two parameters can be calculated simultaneously by the embedded window method or the C-C method proposed by H.S. Kim. Compared with the embedded window method and the C–C method, the amount of computation in the C–C method is less than the embedded window method, and so it is easier to be implemented on the engineering. Takens et al.¹⁴ assume that we have an infinite noise-free dataset, in which case, we can arbitrarily choose the time delay τ. However, as real datasets are finite and noisy, the choice of the time delay plays an important role in the reconstruction of the oscillator from the scalar time series. If τ is too small, the reconstructed oscillator is compressed along the identity line which is called redundancy. If τ is too large, the attractor dynamics may become causally disconnected which is called irrelevance.

CD

The CD can be obtained by the correlation integral function, i.e.

C (r) = 2 M (M - 1) \sum_{i = 1}^{M} \sum_{j = 1}^{M} H (r - ‖ X_{i} - X_{j} ‖)

(2)

where C(r) denotes the CD of the data and M = N – (m – 1)τ denotes the number of points in the phase space. H(u) denotes Heaviside function and it can be defined as

H (u) = {0; u < 0 1; u \geq 0

(3)

r denotes the radius of phase space ball; ‖X_i – X_j‖ denotes the Euclidean distance of the vector X_i and X_j and it can be defined as

‖ X_{i} - X_{j} ‖ = | \sum_{l = 0}^{m - 1} (x_{i + l τ} - x_{j + l τ})^{2} |^{12}

(4)

The correlation integral function C(r) is defined as the probability that the distance between the conjugate vector (X_i and X_j) is not greater than r, and

C (r) \propto r^{D_{2}} (r \to 0)

(5)

where D₂ denotes the CD and it can be obtained by the following equation

D_{2} = {lim}_{r \to 0} (In (C (r)) In (r))

(6)

In the calculation, we can let m in equation (1) change bigger until D₂ is invariant, namely, the linear curve of the diagram In(C(r)) – In(r). Then, the least square method is used to calculate the slope of the linear curve D₂.

The phase trajectory of Lorenz system is shown in Figure 1. There are three invariants in the typical chaotic system, i.e. x, y, z. Among them, the time series curve of invariant y also includes the chaotic characteristic of the Lorenz system as shown in Figure 2. Therefore, based on the above-mentioned method, we can get the CD by using the time series of invariant y. The process of calculation of the CD is shown in Figure 3. In the figure, the slope of the curve In(C(r)) – In(r) labeled in the linear scaleless range is the CD.

Figure 1.

The phase trajectory of Lorenz system.

Figure 2.

The time series curve of the variable y in the Lorenz system.

Figure 3.

The process of calculation of the CD.

In practice, the fault signal will be modulated by the noise, and the robustness of the maximum Lyapunov exponent calculation method must be considered. Because the Wolf method is not suitable for the noised time series and the real data are not enough to be obtained, the small data method is proposed in this paper. Furthermore, the Lyapunov exponent is sensitive to the time delay τ and embed dimension m, the C–C method is used to calculate both the time delay τ and embed dimension m simultaneously which is more accurate than calculating the two items separately.

The maximum Lyapunov exponent

Rosenstein et al.¹⁵ put forward to the small dataset method to calculate the LLE. This method can make full use of the usable data, besides it has the advantage of fast operation and easily implementation. The concrete procedure is shown as follows:

Step 1

Calculate the average period p of the time series {x(n)|n = 1, 2, … , N} by Fast Fourier Transition method.

The time delay and embed dimension are accordance with C–C method as follows:

As in equation (2), the correlation integral function can be obtained as

C (m, N, r, τ) = 2 M (M - 1) \sum_{i = 1}^{M} \sum_{j = 1}^{M} H (r - ‖ X_{i} - X_{j} ‖)

(7)

As in equation (6), we can obtain the CD

D (m, τ) = {lim}_{r \to 0} (In (C (m, r, τ)) In (r)), In (C (m, r, τ)) = {lim}_{N \to \infty} C (m, N, r, τ)

(8)

Let the time series t be misrelated subsidiary sequence, the length is Int(N/t), for any t, we have

{x (1), x (t + 1), x (2 t + 1), \dots} {x (2), x (t + 2), x (2 t + 2), \dots} {x (t), x (t + t), x (2 t + t), \dots}

(9)

Then, we can get the statistics amount of the subsidiary sequence

S (m, N, r, τ) = 1 t \sum_{l = 1}^{t} (C_{l} (m, N r, r, τ) - [C_{l} (1, N r, r, τ)]^{m})

(10)

where C_l denotes the correlation integral of the lth subsidiary sequence.

We can define the D value

Δ S (m, t) = max [S (m, N, r_{i}, t)] - min [S (m, N, r_{j}, t)], i \neq j

(11)

m ∈ [2,3,4,5], r = [0.5Δ,Δ,1.5Δ,2Δ], Δ denotes the average variance of time series. We can get

{S_{cor} (t) = Δ \overset{\land}{S} (t) + | \overset{\land}{S} (t) | Δ \overset{\land}{S} (t) = 14 \sum_{m = 2}^{5} Δ S (m, N, t) \overset{\land}{S} (t) = 116 \sum_{m = 2}^{5} \sum_{j = 2}^{4} S (m, N, r_{j}, t)

(12)

where Ŝ(t) denotes the average value of the statistics amount of the subsidiary sequence S(m, N, r_j, t), the first minimum of ΔŜ(t) corresponds with τ, the minimum value S_cor(t) is the time delay window.

Step 2

In the phase space reconstruction, the algorithm locates the nearest neighbor of each point on the trajectory. The nearest neighbor X(ĵ) is found by the way of searching for the point, namely minimizing the distance to the particular reference point X(j). This can be expressed as

d_{j} (0) = {min}_{j} ‖ X_{j} - X_{j} ‖, | j - j | > p

(13)

where ‖•‖ denotes the Euclidean norm and p denotes the mean period of the time series.

Step 3

For every reference point X(j ), we can calculate the ith time step dispersion distance d_j(i) of the neighborhood point X(ĵ )

Step 4

Suppose the reference point X( j) and the nearest point X(ĵ ) have the exponent divergence ratio λ₁, so we can get

d_{j} (i) = d_{j} (0) e^{λ_{1} (i \cdot Δ t)}

(15)

Step 5

Equation (15) is fetched In(·) on both sides, we can get that

{Ind}_{j} (i) = {InC}_{j} + λ_{1} (i \cdot Δ t)

(16)

From equation (16) we can see that the curve i ∼ lnd_j(i) satisfies the linear relation, the slope of the curve is λ₁Δt. So fix i, for all j corresponding with lnd_j(i) seek the average and divide to Δt, we can get y(i)

y (i) = 1 q Δ t \sum_{j = 1}^{q} {Ind}_{j} (i)

(17)

where q is held constant of d_j(i).

Step 6

The linear region of the curve i ∼ y(i) is selected, and the slope of the curve is made by the least square method, which is the Largest Lyapunov exponent (LLE).

λ_{1} = \sum_{i} i \cdot y_{i} - \bar{y} \sum_{i} i \sum_{i} i^{2} - \bar{i} \sum_{i} i, \bar{i} = 1 n \sum_{i} i, \bar{y} = 1 n \sum_{i} y (i)

(18)

Then we take the Lorenz system, for example; Figure 4 shows the calculation method diagram of the LLE obtained by small data method. The core problem is to find the linear scaleless range as shown in Figure 5; the slope of the linear is the maximum Lyapunov exponent.

Figure 4.

The largest Lyapunov exponent obtained by small data method.

Figure 5.

The comparison between the real curve and the line obtained by least square method.

The experiment data collection and pretreament

A mechanical test bed in the RCM laboratory of Mechanical Engineering College is used in this research to validate the effectiveness of the proposed method in this paper. The engine investigated is F3L912 type of diesel engine, the vibration signal is collected by piezoelectric acceleration sensor which is installed on the cylinder head, and the rotation rate is collected by photoelectric speed sensor. The signal collection system is the system designed by NI company; the collected software is programmed by Labview software.

According to the need of experiment, two kinds of failures are implemented in the first diesel engine block with the engine idling: the one failure is fuel supply advances angle abnormalities, and the other failure is supplies pressure abnormalities. While the F3L912 type of diesel engine is in the normal circumstances, the fuel supply advance angle is 18°A and the injection pressure is 17 MPa. When the fuel supply advance angles are 12, 14, 16, 20, 22, 24, and 26°A and the injection pressures are 11, 13, 15, 19, and 21 MPa, the changes in the maximum Lyapunov exponents and the CD values can be used to analyze the relationship.

Because fuel system fault mainly affect the engine combustion, fuel injection system fault will result in abnormalities of the instantaneous pressure in the cylinder, and the exception will be led directly to the cylinder head on explosion pressure vibration response is abnormal, so we to cylinder head vibration signals as the main object of acquisition and analysis. Figure 6 shows the vibration acceleration diagram in normal state (the injection pressure is 17 MPa, the fuel supply advance angle is 18°A, and the rotation rate is 900 r/min). Figure 7 shows the vibration acceleration diagram with injection pressure fault (the injection pressure is 11 MPa, the fuel supply advance angle is 18°A, and the rotation rate is 900 r/min).

Figure 6.

The vibration acceleration in normal state.

Figure 7.

The vibration acceleration diagram with injection pressure fault.

Then, we can use the above-mentioned method to achieve the characteristic extraction. From the comparison with Figures 4 and 8, we can see that it is hardly for us to locate the linear scaleless region as shown in Figure 9. The disabled reason of the method is the noise. In fact, although the method can make noise reduction to some extent, the method is effective in white noise. From Figures 8 and 9, we can see that we need a new method to reduce the noise reflection.

Figure 8.

The largest Lyapunov exponent of initial data without noise reduction.

Figure 9.

The linear scaleless range selection.

The frequency extent in the combustion process of the engine is 60 Hz–6 kHz. The low frequency range (≤300 Hz) reflects the engine maximum combustion pressure. The medial frequency range (300 Hz–2 kHz) reflects the rising rate of the engine maximum combustion pressure. The high frequency range (≥2 kHz) reflects the maximum acceleration of the engine combustion pressure. The two kinds of abnormalities will result in the medial frequency range changing. Therefore, we choose the wavelet denoise method to make the noise reduction in the range of the frequency (≤300 Hz and ≥2 kHz).

Figure 10 shows the vibration acceleration diagram in normal state by wavelet denoises. Figure 11 shows the vibration acceleration diagram with injection pressure fault by wavelet denoise. From the comparison with Figures 8 and 12, we can see that the fluctuant scope slow downs in the rising range by the wavelet denoise method which is shown in Figure 13. Therefore, we can obtain the maximum Lyapunov exponent by the mean square least method as shown in Figure 14.

Figure 10.

The vibration acceleration diagram in normal state by wavelet denoise.

Figure 11.

The vibration acceleration diagram with injection pressure fault by wavelet denoise.

Figure 12.

The largest Lyapunov exponent of initial data with wavelet denoise.

Figure 13.

The linear scaleless range selection.

Figure 14.

The comparison between the real curve and the line obtained by least square method.

Analyses of the diesel engine fuel system fault vibration signal

CD analysis

The improved G–P method was introduced to calculate the CD valves with different injection pressures and four different speeds based upon the preliminary phase space reconstruction of the diesel engine fuel system fault vibration signal. The changes in the CD values with the injection pressure are shown in Figure 15. From Figure 15, we can see that (1) in the circumstance of injection pressure normal, the value of CD is added with enhancing the rotating speed. (2) In the same rotating speed, the changes of CD were not clearly between the normal and abnormal.

Figure 15.

The correlation dimension values with the injection pressure.

The same method was introduced to calculate the CD valves with fuel supply advance angle. The changes in the CD values with the injection pressure are shown in Figure 16. From Figure 16, we can see that (1) CD values have the tendency of improvement to the constant value by enhancing the rotating speed on the normal situation. (2) On the same level of rotating speed, the changes of CD were not clearly between the normal and abnormal.

Figure 16.

The correlation dimension values with different fuel supply advance angles.

In conclusion, CD values computed with the original signals cannot distinguish whether the diesel engine fuel system operates properly or not.

The maximum Lyapunov exponents analysis

The smaller dataset method was introduced to calculate the maximum Lyapunov exponents with different injection pressures and four different speeds based upon the preliminary phase space reconstruction of the diesel engine fuel system fault vibration signal. The changes in the maximum Lyapunov exponents with the injection pressure are shown in Figure 17. From Figure 17, we can see that (1) the maximum Lyapunov exponents with different injection pressure are the smallest in the proper situation (17 MPa). (2) In the proper situation (17 MPa), the maximum Lyapunov exponent will decrease with the enhancing of the rotating speed.

Figure 17.

The maximum Lyapunov exponents with different injection pressures.

The same method was introduced to calculate the maximum Lyapunov exponents with fuel supply advance angle. The changes in the maximum Lyapunov exponents with the injection pressure are shown in Figure 18. From Figure 18, we can see that (1) the maximum Lyapunov exponent with different fuel supply advance angles is the smallest on the proper situation (18℃A). (2) The changing tendency of the maximum Lyapunov exponent with different fuel supply advance angles is similar with different rotation speeds.

Figure 18.

The maximum Lyapunov exponents with different fuel supply advance angles.

Numerical analysis

Now we can distinguish between the normal signal and the fault signal by the above-mentioned analysis of the diesel engine fuel system fault vibration signal. Furthermore, we need to distinguish between the two kinds of abnormalities. We can use the numerical analysis to achieve the aim.

From Figure 19, we can see that (1) the engine is in the normal state (the injection pressure is 17 MPa and the fuel supply advance angle is 18℃A). In this state, the mean value of the maximum Lyapunov exponent is 0.00020, and the mean variance is 5.236 × 10⁻¹⁰. (2) The engine is in the fault mode 1 state (the injection pressure is 11 MPa and the fuel supply advance angle is 18℃A). In this state, the mean value of the maximum Lyapunov exponent is 0.00035, and the mean variance is 7.8919 × 10⁻⁶. (3) The engine is in the fault mode 2 state (the injection pressure is 17 MPa and the fuel supply advance angle is 22℃A). In this state, the mean value of the maximum Lyapunov exponent is 0.00025, and the mean variance is 3.12 × 10⁻⁸.

Figure 19.

Changes in the maximum Lyapunov exponents with different fault modes.

Conclusions

When the fuel supply advance angle and supply pressure operate improperly, the maximum Lyapunov exponents computed with the original signals have the tendency of improvement. And then it can be used to distinguish whether the diesel engine fuel system operates properly or not.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Haykin S and Princip J. Making sense of a complex world. IEEE Signal Process Mag 1998; 15: 66–81.

Jiang WL, Zhang SQ and Wang YQ. Fault Information Diagnosis based upon Chaos and Wavelets. Beijing: National Defense Industry Press, 2005.

Li J and Zhao JM. Incipient gearbox fault diagnosis based on the reverse state transformation of the chaotic Duffing oscillator and sampling integral technology. Math Prob Eng 2015; 2015: 1–14.

Wang BC, Ren ZH and Wen BC. Chaotic parameters analysis in fault diagnoses. Mach Tool Hydraul 2010; 38: 144–147.

Zheng HM, Li YA and Chen L. Noise reduction of ship signals based on the local projective algorithm. J Northwest Polytech Univ 2011; 29: 569–574.

Tan JQ, E JQ, Tan Q, et al. Chaos diagnosis system of vibration fault for vehicle diesel engine and its application. J Central South Univ (Sci Technol) 2011; 42: 984–989.

Ma JC, Liu Q and Ma YL. Detection of chaotic vibration signals of aeroengines with blind separation. Acta Aeronaut Astronaut Sin 2009; 30: 81–85.

Wu ZY, Yuan HQ and Li L. Fault diagnosis of engine based on chaotic features and support vector machine. J Mech Strength 2010; 32: 723–728.

Yilmaz D and Guler NF. Analysis of the Doppler signals using largest Lyapunov exponent and correlation dimension in healthy and stenosed internal carotid artery patient. Digital Signal Process 2010; 20: 401–409.

10.

Cai YP, Li AH and Shi LS. Processing method for end effects of EMD based on improved chaos forecasting model and its application in machinery fault diagnosis. J Vib Shock 2011; 30: 46–352.

11.

Wang BC, Ren ZH and Wen BC. Fault diagnoses method of rotating machines based on nonlinear multi-parameters. J Mech Eng 2012; 48: 63–69.

12.

Lv JH, Lu JA and Chen SH. Analysis on chaotic time series and its application. Wuhan University School of Public, 2005.

13.

Packard NH, Crutchfield JP and Farmer JD. Geometry from a time series. Phys Rev Lett 1980; 45: 712–716.

14.

Takens F, Rand DA and Young LS. Dynamical systems and turbulence. Lect Notes Math 1981; 898: 336–341.

15.

Rosenstein MT, Collins JJ and De L. Reconstruction expansion as a geometry-based framework for choosing proper delay times. Phys D Nonlinear Phenom 1994; 73: 82–98.