Timely online chatter detection in end milling process

Mechanical Systems and Signal Processing 75 (2016) 668–688

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Timely online chatter detection in end milling process Yang Fu a,1, Yun Zhang a,n, Huamin Zhou a,1, Dequn Li a,1, Hongqi Liu b,2, Haiyu Qiao a,1, Xiaoqiang Wang a,1 a State Key Laboratory of Material Processing and Die & Mold Technology, Huazhong University of Science and Technology, Wuhan, 430074, PR China b National NC System Engineering Research Center, Huazhong University of Science and Technology, Wuhan, 430074, PR China

a r t i c l e in f o

abstract

Article history: Received 14 February 2015 Received in revised form 4 January 2016 Accepted 8 January 2016 Available online 29 January 2016

Chatter is one of the most unexpected and uncontrollable phenomenon during the milling operation. It is very important to develop an effective monitoring method to identify the chatter as soon as possible, while existing methods still cannot detect it before the workpiece has been damaged. This paper proposes an energy aggregation characteristicbased Hilbert–Huang transform method for online chatter detection. The measured vibration signal is firstly decomposed into a series of intrinsic mode functions (IMFs) using ensemble empirical mode decomposition. Feature IMFs are then selected according to the majority energy rule. Subsequently Hilbert spectral analysis is applied on these feature IMFs to calculate the Hilbert time/frequency spectrum. Two indicators are proposed to quantify the spectrum and thresholds are automatically calculated using Gaussian mixed model. Milling experiments prove the proposed method to be effective in protecting the workpiece from severe chatter damage within acceptable time complexity. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Chatter detection Hilbert–Huang transform Energy aggregation Online

1. Introduction Chatter is a severe self-excited vibration in end milling process which brings the system to instability, resulting in great damage to the surface finish integrity of the workpiece, shortened life of the cutter and machine tool, and harsh noises [1–3]. Due to the time-verifying properties and complicated structures of the cutting system, chatter has been one of the most unexpected and uncontrollable phenomenon during the milling operation. With the rapid development of high speed milling, this problem becomes more and more prominent. Since 1907 when Taylor described machine chatter as the “most obscure and delicate of all problems facing the machinist” [4], regenerative chatter has been analytically studied by many researchers, including Altintas [3], Budak [5], Tlusty [6] etc. Because of the tight coupling and time-verifying properties of the whole cutting system, analytic method, like stability lobe analysis [3], still could not identically model the real cutting system and perfectly prevent the occurrence of chatter [1,7]. Therefore the online signal-based cutting state monitoring and chatter detection methods [8,9], which can avoid the complicated modeling of the whole machining system, are very

n

Corresponding author. Tel.: þ86 27 87543492; fax: þ 86 27 87554405. E-mail addresses: [email protected] (Y. Fu), [email protected] (Y. Zhang), [email protected] (H. Zhou), [email protected] (D. Li), [email protected] (H. Liu), [email protected] (H. Qiao), [email protected] (X. Wang). 1 Tel.: þ 86 27 87543492; fax: þ86 27 87554405. 2 Tel.: þ 86 27 87542613; fax: þ 86 27 87554405. http://dx.doi.org/10.1016/j.ymssp.2016.01.003 0888-3270/& 2016 Elsevier Ltd. All rights reserved.

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

669

important for practical machining processes to ensure the safety of the machining system and the workpiece [10], and that is to be focused on in this paper. Propelled by the advancement of computer and sensor technologies [11], online monitoring methods have become more and more powerful [1]. In the past few years, many researchers has conducted a lot of work to acquire chatter-concerned signals and many sensors and signals have been applied in chatter detection, including cutting forces [12,13], vibration [8,14,15], servo current [16], sound [17,18], acoustic emission [15] etc., which have provided substantial information about the machining process. Among these acquired signals, vibration signal stands out for its comprehensive performances in price, installation, signal quality and information usability. Recently, Quintana et al. [1] reviewed the state of the art in chatter research and they pointed out that an efficient signal processing algorithm is another critical issue to be dealt with to identify the chatter before it is completely developed. Because vibration signals are typical non-stationary signals, time frequency analysis (TFA) spectral analysis method should be engaged to process the measured signals. TFA methods, like wavelet method, Hilbert–Huang transform, STFT, WignerVille method et al., map the one-dimensional signal into a two-dimensional time–frequency plane which can simultaneously provide both the frequency distribution and the associated time information. Wavelet transform and Hilbert–Huang transform are two most used TFA methods in recent years. The wavelet transform decomposes the signal into a series of sub-signals and the original signal can be perfectly reconstructed using the wavelet coefficients. These sub-signals then can be used for feature extraction. Kuljanic et al. proposed a multisensory system for milling chatter detection using the wavelet method [8,15]. In their works, the signal was decomposed by db8 mother wavelet and the energy ratio of aperiodic component to periodic component was used as one of the indicators. Yao et al. identified chatter by combining the wavelet transform and the support vector machine [19]. Wavelet transform has been proved to be a powerful time frequency analysis method which has an integrated description in both time domain and frequency domain. However, the choice of mother wavelets is still empirical and there are no standards or general methods on how to appropriately make the selection [20]. Different wavelets would lead to totally different performance. That is why the wavelet transform method remains a manual method and cannot be widely used in online monitoring applications [20]. Hilbert–Huang transform (HHT) is another important TFA method invented by Huang et al. in 1998 [21]. It has been widely used in the field of non-stationary and non-linear signal processing and feature extraction [22–24], including geology [23], oceanography [25], fault diagnosis [24] etc. HHT uses the time-adaptive empirical mode decomposition (EMD) method to decompose the signal into a set of complete and almost orthogonal intrinsic mode functions (IMFs), and then Hilbert spectral analysis (HSA) is applied to obtain instantaneous frequency array. Unlike the Fourier transformation method, HHT does not involve the concept of frequency resolution or time resolution, but introduces the concept of instantaneous frequency. Consequently, a uniform high resolution is obtained in the full frequency range [22]. Carbajo et al. developed a frequency track method using a windowed Hilbert–Huang transform and validated it to be reliable by experiments [26]. Peng et al. presented a comparison study between the Hilbert–Huang transform and wavelet transform in the application of fault diagnosis in rolling bearing [24]. Recently, Lei presented a review on the application of empirical mode decomposition in fault diagnosis of rotation [27]. HHT is self-adaptive and does not involve too much manual intervention which makes it to be a more suitable method for online chatter signal processing. After spectral preprocessing, the following task is to model the relationship between cutting signals and chatter, which has been conducted using different feature extracting methods [28]. Schmitz proposed a method for chatter recognition using statistical evaluations of the milling sound variance with a synchronously sampled signal [18]. Zhang et al. came up with a hybrid approach for chatter monitoring based on the hidden Markov model and artificial neural networks [29]. Liu et al. proposed a method for turning chatter identification using support vector machine based on energy and kurtosis index of intrinsic mode functions [16]. In general, most of them just constructed discriminative models and checked the accuracy, which is not enough for practical online applications. In summary, signal selection, signal processing and feature extraction are the most critical steps in chatter detection, and all these steps have been studied extensively. Unfortunately, the existing methods still cannot identify the chatter before the workpiece has been damaged which should be the core objective of chatter detection, as pointed out by Quintana et al. [1]. This is probably caused by their feature constructing strategy. They generally monitored a certain natural frequency component or frequency band obtained from the modal analysis, while the noisy machining environment and the time-varying property of the system mode will always damage the performance. Chatter is a phenomenon reflecting changes of frequency and energy distribution in machining process [30]. From this perspective, features, which can reflect the energy distribution transform, may achieve some improvements on the timeliness of the detection. Another problem is the need of manual intervention in the model constructing and adjusting which makes the existing methods impractical. Especially when wavelet method is involved, experts will always be needed to set and adjust those model parameters. The preparation of the training samples of machining states also need a lot of meticulous manual work. The last deficiency is a lack of time complexity analysis of the detection algorithm which is a key factor for an online monitoring system. Considering above deficiencies, an energy aggregation characteristic-based Hilbert–Huang transform method is proposed to characterize the progressive change of the vibration energy distribution in this paper. Ensemble empirical mode decomposition is introduced to enhance the narrow band property of the intrinsic mode functions. Features will be extracted from the modified Hilbert–Huang spectrum to indicate the cutting state. Then Gaussian mixture model will be employed to model the feature space, and an automatic threshold determination method is also provided. Experiment

670

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 1. The experimental platform.

Table 1 Configurations of the cutting tool. Terminology

Values

Tool cutting edge inclination (helix angle) (degree) Tool diameter (mm) Tool blade length (mm) Number of teeth Tool length (mm) Tool overhang

30 10 25 3 75 50

evidence will be presented to prove the performance of the proposed detection method in protecting the workpiece from the chatter damage. Time complexity analysis also demonstrates the feasibility. The experimental setup and the method details will be presented in Section 2. Section 3 will discuss the results and conclusions are drawn in Section 4.

2. Methodology 2.1. Experimental set-up In order to support the proposed method for chatter detection, a series of end milling experiments were conducted. Experiments were performed on a TC500 Drilling & Tapping Center, which is able to operate at a high milling speed of about 15,000 rpm. No cutting fluid was used. The workpiece was an Al6061 brick claimed on the vice. Accelerometers (PCB 356A15 3D 2–5 kHz 75%) were mounted on the spindle housing to measure the vibrations during milling. LMS SCADAS Lab was employed to sample vibration signals and transmit them to a laptop. Cutting experiments were conducted by straight cut. The signal acquisition setup is described in Fig. 1. 2.2. Test configurations As a preliminary research, we chose the workpiece to be AL6061, which was widely used in cutting vibration research. The tool holder was examined by dynamic balance test to ensure the safety of the experimentation. Uncoated M2Al tools with three teeth were used in the experiment, and the detailed characteristics parameters are listed in Table 1. Impact testing was conducted and results were illustrated in Fig. 2. It was identified that 3324 Hz is one of the natural frequency of the cutting system. 2.3. Chatter tests Chatter is closely related to the spindle speed and depth of cut according to Altintas [3], so in order to activate the phenomenon, we purposely varied the spindle from 3000 rpm to 10,000 rpm and the depth of cut from 1 mm to 10 mm. For

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

671

Fig. 2. Impact testing.

Table 2 Cutting conditions for the chatter tests. (I for idling moving state, S for stable cutting, C for Chatter). No.

Position

Cutting speed (m/min)

Feed rate (mm/min)

Axial depth (mm)

Radial depth (mm)

Revolution speed (rpm)

Note

1 2 3 4 5 6 7 8 9

R-y L-y L-y L-y L-y L-y L-y L-y L-y

94.2 314.0 314.0 314.0 314.0 314.0 314.0 282.6 251.2

300 300 600 200 300 400 400 400 400

10 10 10 10 10 4 6 6 6

1 0.5 1 1 1 1 1 1 1

3000 10,000 10,000 10,000 10,000 10,000 10,000 9000 8000

I–S I–S S–C I–C I–S–C I–C I–C I–C I–C

end milling, the radial cutting depth greatly influences the cutting forces, so we also conducted some experiments using different radial depth of cut. Because this paper is focus on the detection of chatter, we picked out nine representative results and the parameters are summarized in Table 2. There is an interesting phenomenon in the No. 3 case. When the cutting was setup, it turned out to be stable for a while, and suddenly transformed to severe chatter. This makes it to be a perfect material for this research. The measured signal and the surface topography are illustrated in Fig. 3. And the other eight are displayed in Fig. 4.

2.4. Energy aggregation characteristic-based Hilbert–Huang transform 2.4.1. Energy aggregation property of vibration signal during chatter development Regenerative chatter is a well-known self-excited vibration caused by the interaction of the workpiece and the cutter, represented by changes of frequency and energy distribution in milling process. Siddhpura et al. stated that modulated chip thickness due to vibration affects cutting forces dynamically, which in turn, increase vibration amplitudes [30]. When the modulation frequency of the cutting forces gets close to one of the natural frequencies of the cutting system, the system will become unstable and violently resonate at the certain frequency, which is the so-called chatter [1]. The vibration extracts energy to start and grow from the cutting system. Therefore, the progressive aggregation of energy at a certain frequency is the essential characteristics of chatter onset. In the frequency spectrum, chatter is a small peak at the beginning and then transforms to be dominant in a very short time. When chatter goes to be severe, most of the vibration energy will be absorbed into the certain frequency, therefore the vibration signal will perform narrow-band which can be seen as the

672

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 3. Measured vibration signal (a) and the corresponding surface topography of case 3.

Fig. 4. Measured vibration signals in Case 1 (a), 2 (b), 4 (c), 5 (d), 6 (e), 7 (f), 8 (g) and 9 (h). The stairs in (a), (b), (c), (e), (f), (g), (h) and the first one in (d), represent the change that the cutter and the workpiece got contacted. And the second stair in (d) indicates a cutting state transition from stable state to chatter.

frequency feature of chatter. If the signal is pure, a modulated waveform will be detected whose high frequency component is just the chatter part. Fig. 5 illustrates a measured vibration segment and five frequency spectrums at different moments extracted from the No. 4 case. When the process is stable at 0.082 s, the frequency is flat and the energy is scattered in a wide band. After a short budding stage from 0.78 s to 1.02 s, the energy of the main chatter frequency peak increase sharply in only 0.14 s, from 1.02 s to 1.16 s. From the aspect of energy, the whole transition can be regarded as a convergence process of the signal energy to a certain narrow frequency band where the chatter frequency located. The energy aggregation is a more general description of the chatter phenomenon from the energy perspective. This property will be used in this paper to detect the milling chatter with an energy aggregation based Hilbert–Huang transform which can adaptively decompose the measured vibration signal into narrow band according to energy distribution and effectively quantify the aggregation degree to indicate the chatter evolution.

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

673

Fig. 5. Illustration of the energy aggregation phenomenon during chatter development. (Sample length: 0.0125 s with 256 point).

2.4.2. Hilbert–Huang transform and ensemble empirical mode decomposition improvement In order to characterize the energy aggregation transition, Hilbert–Huang transform (HHT) method is employed in the research. HHT method includes two steps, empirical mode decomposition (EMD) and Hilbert spectral analysis (HSA). The decomposing results of the EMD are a series of complete and almost orthogonal components called intrinsic mode functions (IMFs), expressed as xðt Þ ¼

N X

ci ðt Þ þ r n ðt Þ

ð1Þ

i¼1

where x(t) is the measured signal, ci(t) is the IMFs and rn(t) is the residual. Then Hilbert transform is conducted on the IMFs to calculate the Hilbert spectrum. Z þ1 1 ci ðτÞ H ðci ðt ÞÞ ¼ P dτ ð2Þ π 1 t τ where P represents Cauchy principal value integral. Then the instantaneous frequency can be calculated as f i ðt Þ ¼

1d H ðci ðt ÞÞ arctan 2 dt c i ðt Þ

ð3Þ

The Hilbert–Huang spectrum then can be denoted as 8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < c2 ðt Þ þ H 2 ðc ðt ÞÞ; f ¼ f ðt Þ i i i H i ðt; f Þ ¼ : : 0; f af ðt Þ

ð4Þ

i

The marginal spectrum can be calculated as follow if necessary. Z þ1 H ðt; f Þdt H ðf Þ ¼

ð5Þ

1

The theory and detailed derivation of the method can be referred in [31]. However, there exists a critical defect that the first decomposed IMF usually involves a wide range of frequencies which makes the followed HSA to be meaningless. To avoid the mode mixing problem, ensemble empirical mode decomposition (EEMD) is introduced to conduct the decomposition, instead of EMD. EEMD is a noise-assisted data analysis method. Finite white noise is added to the investigated signal and EMD is applied on the noisy signal. By averaging the IMFs for enough cycles of the above steps, the mode mixing problem can be eliminated automatically. EEMD method is a significant improvement of the EMD [32]. The performance of the EEMD applied on the vibration signal will be discussed in the Section 3. To apply EEMD method, three important parameters, noise level, ensemble number and stoppage criteria, are necessary. Noise level is the amplitude of the white noise added to the signal before applying EMD. It is recommended to be about 0.2 standard deviation of the amplitude of the signal [32]. Through signal tests with different noise level (0.1, 0.2 and 0.4), it seems that a larger noise level performs better in reducing the mode mixing problem for the vibration signals. In this work, the noise level of 0.4 is chosen. In order to reduce the contribution of the added noise to the decomposed results, a large ensemble number, 100, is used [32]. Stoppage criteria determines the number of sifting steps to produce an IMF. A direct

674

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 6. A typical energy limited Hilbert–Huang spectrum of cutting process changing from stable stage to chatter. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

way is to examine the sifting result according to the definition of the IMF given by Huang [21] or using the SD rule [21]. However, when applying EEMD with these rules, the decomposed results and the sifting number perform quite unstable [32], which is a disaster for practical applications. Wu [32] stated that using a fixed sifting number (S-number criteria [40]) is a proper way to solve this problem, which is proved to be quite effective. This S-number rule is adopted in this work and the sifting number is set to 10. Besides, a fixed sifting number is important to control the computational complexity of the whole algorithm for practical applications.

2.4.3. Principal energy based selection Focusing on the goal of indicating energy aggregation, a principal energy based selection is conducted after the EEMD step. The divergence degree of the energy can be expressed using the minimum number of IMFs which contains the majority signal energy. Using ci ; i ¼ 1; 2; …; n; r, which is in descending order according to the total IMF energy, to represent the intrinsic mode functions (IMFs) and the residual, then the normal energy NEi of each IMF can be expressed as X X X X  ð6Þ NE ¼ normalize c21 ; c22 ; …; c2n ; c2r where n is the number of IMF and cr represents the residual. As ci are in descending order, the NEi are also descending sorted. Introducing an energy limit coefficient (ELC) Elim, the top m IMFs are selected according to the following criteria: min m;

s:t:

m X

NEj 4 Elim

ð7Þ

j¼1

Because the IMF is narrow-band, if m is quite small, even one, this represents the majority energy of the vibration system are gathering in a quite narrow frequency band, which meets the above aggregation property description when chatter occurs. The selected m IMFs are named feature IMFs, and then input to the Hilbert spectrum analyzer to calculate the instantaneous frequency and the related amplitude and time. Therefore an energy limited Hilbert–Huang spectrum (EL-HHS) is obtained. Some indicators will be extracted to build the identification system which will be given in Section 2.5. Fig. 6 is a typical energy limited Hilbert–Huang spectrum (EL-HHS) of a milling process developing from stable stage to chatter. It should be stressed that the different colors represent different IMFs rather than aptitudes in the traditional time/ frequency spectrum. The total normalized energy of all the frequency spots is limited by the energy limit coefficient (ELC). The figure shows the change of energy and frequency distribution over time. Before 2.9 s, the majority energy distributes in three feature IMFs, represented in blue, green and red spots respectively. Three colored spots scatter in a broad frequency band. This means that the cutting energy is dispersed and milling process is in a relatively stable stage. Around 2.94 s, a sudden change occurs. The number of feature IMFs falls to only two and then single one and the frequency spots are gathering around 3200 Hz. The energy appears to be attracted to a central frequency and the milling process goes into a prophase of chatter. The prophase lasts about 30 ms and then the frequency distribution gets aggregated and goes out of stable stage. Chatter is under way. 2.5. Chatter indicators Chatter is well known to be a random and chaotic phenomenon. Although the energy limited Hilbert–Huang spectrum (EL-HHS) can intuitively reveal the energy aggregation progress during machining, a precise quantitative characterization of the spectrum is still needed to automatically detect the chatter onset. An exact and robust detection system is very important for practical application. From the above point of chatter identification based on the energy aggregation property, two major indicators are introduced and several auxiliary indicators are recommended to enhance the performance.

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

675

2.5.1. Normalized energy ratio of the most powerful IMF, NER The normalized energy ratio of the most powerful IMF indicates the aggregation degree of the energy. When the cutting is stable, vibration energy distribution is quite uniform and the NER is small. However when chatter begins to develop, NER will increase to nearly 1 in a very short time which is a sign that the system energy is being attracted to a certain frequency and severe resonance is around the corner. NER is a very good indicator to imply the chatter when there is no much external impulse disturbance.

2.5.2. Coefficient of variation of the energy limited Hilbert–Huang spectrum, CV The coefficient of variation (CV) is defined as the ratio of the mean value μ to the stand variance σ of all the frequency points. It is signal-independent and widely used in probability theory and statistics to measure the degree of dispersion [33]. Coefficient of variation of the spectrum implies the degree of frequency scatters around the centric position and represents the aggregation degree from the viewpoint of probability. A proper setting of the CV threshold can help to reduce the interference of the impulse disturbance and enhance the robustness of the detection system. CV ¼

μ σ

ð8Þ

2.5.3. Auxiliary indicators Chatter is a kind of resonance between the cutting system and the workpiece. However the cutting system may also show some slight resonance appearance when it receives a certain external excitation especially when the system vibration amplitude is quite small. As the proposed method is mainly based on the chatter frequency characteristic, this kind of interference will result in confused detection. Activating parameter (AP) is used to activate the detecting method when the actual machining begins and eliminate the perplexity when the lathe is idly moving. Another useful indicator is the central frequency (CF), the mean value of the frequency scatters on the frequency dimension. The CF represents the average frequency of the signal and if the source signal is mono-component or near monocomponent, CF is the dominant frequency. When chatter occurs, the resonance frequency is just the mean value of the frequency points, which is represented by CF. 2.6. Cutting state pattern recognition using Gaussian mixture model Thresholds need be determined to identify the cutting states with the calculated indicators. Gaussian Mixture Model (GMM) is a parametric probability density estimator represented by a weighted sum of Gaussian component densities. GMMs are commonly used as a parametric model of the probability distribution of continuous measurements or features in pattern recognition system [34]. Using the Expectation–Maximization (EM) algorithm from a well-trained prior model, GMM parameters can be estimated. Using the well estimated GMM, threshold can be obtained. A Gaussian mixture model is a weighted sum of M component Gaussian densities, expressed as M   X   ωi g xjμi ; Σi p xjθ ¼

ð9Þ

i¼1

where x is a D-dimensional continuous-valued data vector (i.e. measurement of features), ωi ; i ¼ 1; 2; …; M are the mixture   weights, and g xjμi ; Σi ; i ¼ 1; 2; …; M are the component Gaussian densities, given as Eq. (10).    0  1   1 1 x  g xjμi ; Σi ¼ exp  μ Σ x  μ ð10Þ   i i i 1=2 2 ð2π ÞD=2 Σi  M P ωi ¼ 1: where μi is the mean vector, and Σi is the covariance matrix. The mixture weights satisfy the constrain that i¼1 After the number of components is determined, the complete GMM model is parameterized by the mean vectors, covariance matrices and mixture weights from all component densities. These parameters are collectively represented as θ in Eq. (9).

θ ¼ fωi ; μi ; Σi g;

i ¼ 1; 2; …; M

ð11Þ

A common way to solve the GMM model is Expectation–Maximization (EM) algorithm. Introducing a hidden variable zi to represent the i-th Gaussian component, we can build a condition distribution of the train sample and the Gaussian component using Bayes formula       p z ¼ k; xj jθ ωk g xj jθk ¼ M ð12Þ γ jk ¼ p z ¼ kjxj ; θ ¼ M  P   P  p z ¼ k; xj jθ ωk g xj jθk j¼1

j¼1

676

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 7. Rolling time window sampling strategy.

By maximizing the log likelihood function, we can obtain M P M 1 X ωk ¼ γ ; M j ¼ 1 jk

μk ¼

j¼1

M P

γ jk xj

M P j¼1

γ jk

;

Σk ¼

j¼1



γ jk xj  μk M P j¼1

γ jk

2 :

ð13Þ

The EM algorithm can be summarized as (1) (2) (3) (4)

Initialize θ; Calculate the generation probability γjk using Eq. (12); Calculate the model parameters θ using Eq. (13); Repeat step 2 and 3 until the log likelihood function converge.

With all the above definition, a GMM model is complete and can be used to model the structure represented by the unlabeled train data. In this work, we will employ the GMM method to model different cutting states, and find proper threshold for the online chatter detection. 2.7. Sampling strategy for online chatter detection The efficiency of the real-time system greatly depends on the sample length, the time complexity of the monitoring system and the CPU computing power among which the sample length affects the cycle time cost and the time resolution of the system. In the research, we employ a rolling time window sample strategy. A fixed time length window is set when the application starts. Then a short measured signal is inserted to the window from the end, and the same length signal from the beginning of the window is abandoned. The length of the time window is called frame length and the length of the new segment is called frame shift. This is a quite flexible sample strategy. A large frame length and small frame shift can be set if the CPU is powerful enough and a high time resolution can be obtained. A short frame length and relatively large frame shift is proper when the CPU cannot afford the time complexity, and the time resolution will also reduce. Fig. 7 illustrates the sample principle. 2.8. Parameters summary In order to clarify the presented method, here we summarize all the related parameters in Table 3. The parameters are classified according to processing section they belong to. Some symbols will be proposed in the table and the following derivation will keep to the specified symbols. Finally, a detailed algorithm flow chart of the complete processing procedural is illustrated in Fig. 8.

3. Results and discussion 3.1. Results of the mode decomposition The narrow-band property of the intrinsic mode functions (IMF) is an essential hypothesis of the Hilbert Huang method. It is very important to ensure the narrow-band quality of the decomposing results. The performance of the ensemble empirical mode decomposition (EEMD) has been proved by using simulated signal in many literatures [32,35,36], and in this section, we will focus on the performance of the EEMD in improving the narrow-band property of the IMF when decomposing the measured vibration signal. The mode mixing defect [32] is the main problem affecting the narrow-band performance. According to the existing research, this phenomenon is widespread, especially in the first IMF component [32]. The worst situation is that the first IMF will cover a wide range of frequency, resulting that the energy aggregation representation is negated and the following Hilbert spectral analysis will also be meaningless. The EEMD method is employed here to solve this problem. Fig. 9

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

677

Table 3 Parameters used in the proposed model. Terminology

Symbol Components

Time complexity Notes

Sampling frequency Frame length Frame shift Length of the signal Sift number

Fs tw ts n NS

Sampling Sampling Sampling Sampling EEMD

/ O(n) O(n) / NS*O(n)

Number of IMFs Number of ensemble Noise ratio

NM NE ε

EEMD EEMD EEMD

NM*NS*O(n) NE*NM*NS*O(n) O(n)

Number of feature IMF Energy limitation

m

O(n)

/ Activating parameter Energy ratio

/ AP NER

Energy-based selection Energy-based Selection HSA Identification Identification

Coefficient of variation Central frequency

CV

Identification

O(n)

CF

Identification

O(n)

Elim

/ O(n) O(n) O(n)

Sampling frequency used form measuring the signal. The time length of the rolling window. The update length of the inserted signal. Number of the points of the signal in the rolling window. Sift times in the EMD algorithm to terminate an extracting cycle and determine an IMF. The number of IMFs to be obtained when applying EMD. The cycle times of add noise and apply EMD in the EEMD algorithm. Ratio of the standard deviation of the added white noise and that of the input signal. The number of the feature IMF selected from the EEMD result according to the settled principal energy limitation. Energy limit coefficient used in the IMF selection. Hilbert spectral analysis applied on the selected feature IMFs. AP is used to distinguish the empty moving and the actual cutting. NER physically indicates the aggregation degree of the signal energy along the frequency axis. CV probabilistically implies the degree of frequency scatters around the centric position. The mean value of the frequency scatters which represents the resonance frequency when chatter occurs.

Fig. 8. A complete chatter detection algorithm flow chart.

illustrates decomposing results of a vibration signal segment extracted from the measured vibration signal when the cutting is relatively stable. It can be seen that the main frequency components are decomposed into three different IMFs both by the EMD (blue line) and the EEMD (red line). But there exists an obvious mode mixing problem in the IMF 1 of the EMD. Its frequency range is quite wide and cover the main frequency in the IMF 2 and 3. The IMF1 obtained using EEMD, represented in red line, effectively depresses the extension of the frequency band. More obvious depressing effect can be found in the IMF 3. Fig. 10 is another result showing the depressing effect when chatter occurs. There is no obvious difference in the performance of the IMF 1. Because when chatter occurs, the first IMF which is the most powerful component, absorbs most of the vibration energy to a central frequency. This resonance phenomenon itself depresses the other frequency components. But it should be noted that the depressing effect is clear for the IMF2 and IMF3 where the EEMD results in red lines are beneath the blue ones and occupy a smaller frequency range. Thus the mode mixing problem is lessened. In addition, in

678

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 9. Vibration signals and its EMD and EEMD decomposition results when the cutting is stable. (a) The original vibration signal and its associated FFT spectrum. (b)–(d) the decomposed IMFs from 1 to 3 and the associated FFT spectrum among which the blue line is the results of the EMD and the red line is the results of the EEMD (noise ratio¼ 0.4, NE¼ 100). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Vibration signals and its EMD and EEMD decomposition results when chatter occurs. (a) the original vibration signal and its associated FFT spectrum. (b)–(d) the decomposed IMFs from 1 to 3 and the associated FFT spectrum among which the blue line is the results of the EMD and the red line is the results of the EEMD (noise ratio¼0.4, NE¼ 100). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

other application where the first IMF is not the most powerful one, a sorting according to the IMF energy will be needed before the following processing and similar result can be obtained. It can be concluded that the EEMD method can greatly enhance the narrow-band property of each IMF and effectively tackle the mode mixing problem. This improvement mainly owes to the noise-assistant strategy. It can effectively suppress the signal noise which is the one of the main pollutions.

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

679

Fig. 11. Indicators calculated from a segment extracted from the Case 3. The measured vibration signal (a) and the calculated indicator curves (b) are displayed in the left part. And indicators of three sections are illustrated in the right part. The three sections, named in S1, S2, and S3 relatively represent stable cutting stage, transition stage and chatter. (Fs¼ 12,800 Hz, Elim ¼ 0.9).

3.2. Chatter indicator analysis Normalized energy ratio (NER) and coefficient of variation (CV) are the two main indicators in this research as described in Section 2.5. NER is the normalized ratio of the most powerful IMF and CV is the ratio of mean value to the standard variance of the frequency scatters. Fig. 11 shows the results of NER, CV and central frequency (CF) calculated from a signal segment extracted from the No. 3 case. When the cutting is stable, the vibration energy is dispersive, and the related spectrum shows a lot of decentralized frequency scatters. The NER is low around 0.7 and the CV is high. In this case the CF has no physical meaning. When the cutting turns to be unstable, the energy is gradually absorbed to a central frequency, and the dispersion degree begins to decrease. The disappearance of the background noise scatters and a frequency band occurs ease. The disappearance of the background noise scatters is another sign, and a frequency band gradually forms. The NER begin to grow and the CV decreases. The CF reflects the convergent tendency. When sever chatter occurs, the majority energy gathers around a certain natural frequency, equal to CF, and violently resonance occurs. The NER is very close to 1, such as 0.99 or 0.999, and the CV is relatively quite low. The central frequency (3417 Hz) is also close to one of the natural frequencies of the cutting system, 3324 Hz, obtained from impact test, with considering the signal noise, riding effect and dynamic mode variation. Chatter is underway. The transition can also be found from the two indicators lines showed Fig. 11(b). As chatter develops, the NER line increases quickly to near to 1 and keep steady, and there is a climbing process corresponding to the chatter development. The CV line presents a much more evident drop when the cutting state transform from stable state to chatter. This is mainly caused by the decrease of the number of the feature IMFs. The principle energy based selection filters out those low energy components which can be seen as the background noise compared with the principal components. Because CV is an arithmetic statistical value of the frequency scatters in the spectrum without considering the associated amplitude, the decrease of the number of the feature IMF will remove most of the low energy scatters and the mean value and standard variance will visibly shift, resulting in the jump of the CV, which make CV quite sensitive to the state change. With the NER quickly increasing and the CV jump, it can be confirmed that the cutting system is degenerating to chatter. And if there is no efficient external assistant, violate resonance will be inevitable.

680

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 12. GMM pattern recognition model of the Case 3. (a) feature space; (b) the modeled two Gaussian distribution; (c) the NER view; (d) the CV view.

3.3. Cutting state modeling results In order to obtain a credible threshold to identify the cutting state, GMM method is used to model the feature space constructed by NER and CV. The number of Gaussian components is set according to the signal property, and then the EM method is used to solve the model. The result of the No. 3 case is illustrated in Fig. 12.

 0:3667 0:3498 μ1 ¼ ð 0:9764 0:1060 Þ; Σ 1 ¼ expð  3Þ   0:3498 0:9614

 5:3844 1:6240 μ2 ¼ ð 0:7614 0:4382 Þ; Σ 2 ¼ expð  3Þ  : ð14Þ  1:6240 1:6685 The two solved Gaussian distributions are relatively parameterized as Eq. (14). μ represents the distribution center and Σ the distribution extension. It can be seen from Fig. 12(b) that the two Gaussian components perfectly characterized the two different cutting conditions. The stable cutting spread in a wide range in the feature space and in the GMM it is a broad and flat peak. On the opposite way, the chatter scatters in a quite small range and represented by a cliffy peak in the GMM plane. The GMM method is insensitive to the noise, like the top left outliers, which makes it to be a robust method. The two projection views showed in (c) and (d) provide clearer information about the two components. According to the NER view showed in (c), the two components will be separated around 0.92. Before 0.92, the probability that the sample belong to stable component is larger than the chatter one. When NER is larger than 0.92, it is more likely to indicate the chatter onset. The similar property can be seen in the CV view in (d). The major difference is a wider separation margin. Similar properties can also be found in the other cases. Using the well solved GMM model, thresholds for the online chatter detection can be automatically determined by solving the Eq. (15) in different projection dimensions.     ω1 g xjμ1 ; Σ1 ¼ ω2 g xjμ2 ; Σ2 ð15Þ The GMM spectrum of the other eight cases in Table 2 are illustrated in Fig. 13. Different Gaussian components represent different machining states. It can be seen from (a), (b) and (d) that a flat Gaussian component can be found located in small NER and large CV space, representing the stable cutting state. The chatter state in (c)–(h) is represented by a Gaussian component located in large NER and small CV space. These two states have a clear boundary around (NER ¼0.9, CV ¼0.3). Although the idling moving state presents some uncertainty, it can be easily excluded using an activating parameter to tell the system when cutting starts. From the results, it can be found that the proposed method can exactly identify chatter from stable and idling moving states. From the above results and discussions, initial thresholds with NER¼0.9, CV¼0.3, is a proper set to distinguish the chatter and stable cutting state at the test running. This may result in a little time delay but does not have significant impact

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

681

Fig. 13. Weighted Gaussian components for: (a) Case 1; (b) Case 2; (c) Case 4; (d) Case 5; (e) Case 6; (f) Case 7; (g) Case 8; (h) Case 9. In order to better show the two components, Regularization with 0.001 scale is added when solving the model.

on the detection performance. After the machine has run for a while, these thresholds can be refined based on the GMM results using the measured vibration signals according to the Eq. (15). 3.4. The performance of the presented method in preventing the chatter damage 3.4.1. Performance on the measured signal Using the proposed method, the feature curve can be obtained and illustrated in Fig. 14. Fig. 14(a) represents a stable cut in a low speed at 3000 rpm, the feature curve is noisy but stabilized in the stable region defined by the thresholds; Fig. 14(b) is another stable cutting signal, but the spindle is much higher at 10,000 rpm. The relative signal to noise ratio is much higher, and the NER line and the CV line is much smoother. Fig. 14(c) displays a signal which fall to chatter at the very beginning of the cutting. The NER curve rapidly increases to near 1, and the CV line fall to around 0.2, which is located in the chatter region defined by the thresholds and the chatter is discovered at the budding stage; Fig. 14(d) illustrates a signal which experiences idling moving, stable cutting and chatter. The chatter is exactly identified just when the chatter amplitude begins to increase. Fig. 14(e), (f), (g) and (h) is another four cases for different cutting depth and cutting speed. The eight cases prove the proposed method to be effective in identifying the stable and chatter cutting states using the measured vibration signal.

682

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 14. The feature curves for Case 1 (a), 2 (b), 4 (c), 5 (d), 6 (e), 7 (f), 8 (g) and 9 (h). The measured signal is plotted using the left y axis, and the feature lines are plotted using the right y axis.

Fig. 15. The algorithm performance in protecting the workpiece from the chatter damage in Case 3. (a) the vibration signal and the related NER and CV curve; (b) the surface topography. (feed speed¼ 10 mm/s).

3.4.2. Performance on protecting the workpiece Although the above results show that the proposed method can precisely identify the chatter from the stable cutting, further evidence is still needed to prove whether the proposed detection algorithm can protect the workpiece from the chatter damage. We employ the algorithm to the No. 3 case and results are illustrated in Fig. 15. From Fig. 15(a) chatter can be detected at 1.252s (NER¼0.9225, CV¼0.1562) just when the resonance aptitude starts to grow. Fig. 15(b) shows the associated surface topography, and we can find that before the capturing moment, 1.252 s, the surface is smooth and after the detecting moment, chatter mark turns to be obvious in less than 50 ms, denoting that the

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

683

Fig. 16. Comparisons with wavelet features proposed by Yao. Because the feature scales of T1 and T2 are quite different, they are displayed in different y axis. T1 using the right y axis and T2 the left.

Fig. 17. An example to explain the difference between the HHT and the wavelet. (a) simulated signal; (b) FFT spectrum; (c) Hilbert Huang spectrum; (d) wavelet spectrum.

chatter is more and more violent. The proposed chatter detecting algorithm is proved to be effective in identifying the chatter occurrence at its budding stage and protect the workpiece from the further chatter damage. 3.5. Comparison with wavelet method A comparison is made with wavelet method proposed by Yao [19]. In their work, they proposed two features, standard variation (T1) and energy ratio (T2) of the sub-signal where the chatter frequency located which is obtained by using wavelet package decomposition. The four features, namely NER, CV, T1 and T2, for No. 3 case from 1 s to 1.5 s, are illustrated in Fig. 16. Around the detecting point 1.252 s time predicted by the proposed method, the wavelet features present quite unsatisfying performance. T1 begins to grow after 1.3 s and T2 even perform no more different from the stable cutting around 1.252 s, and even all the time along. A possible reason may be the blurry time border as illustrated in Fig. 17. A numerical example can be used to explain the difference between the HHT and the wavelet method, as illustrated in Fig. 17. Fig. 17(a) shows a simulated sinusoidal signal with two different frequency components. A frequency change from 0.16 to 0.08 occurred at sample 500. The FFT spectrum in Fig. 17(b) clearly shows the two components, but it cannot distinguish the duration of the two frequency component and the changing moment. The HHT spectrum (c) and wavelet scalogram (d) can both tell the frequency change in the time dimension, but there exists obvious difference. The frequency resolution of the wavelet scalogram is clearly lower compared with the HHT spectrum. And the resolution performance of the scalogram varies from the high frequency component to the low frequency component. Besides, the time border is quite vague which greatly influence the timeliness of the detection. On the contrary, the HHT spectrum obviously outperforms both in frequency resolution and time resolution. That is because the HHT method does not involve the concept of the

684

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Fig. 18. HHS results of the measured vibration signals in case 1 using Cao's method. (a) HHS of the raw vibration signals; (b) HHS of the [33] wavelet as proposed in Cao's work; (c) HHS of the [32] wavelet; (d) HHS of the [36] wavelet.

Fig. 19. The energy limited Hilbert Huang spectrum of the vibration signals in Case 1.

frequency resolution and time resolution, but represents the instantaneous frequency and can theoretically obtain good accuracy in both time and frequency resolution, getting rid of the limitation of Heisenberg–Gabor inequality [24]. As for the online cutting state monitoring, a high time resolution is very important, besides a high frequency resolution, to help the system to make timely decision. Another comparison is conducted with the HHT method with WPT pre-denoising proposed by Cao [37]. Db33 mother wavelet is used in the decomposition of the No. 3 case. Results are illustrated in Fig. 18. The four Hilbert Huang spectrums (HHS) all reveal some differences between the stable cutting and chatter. In Cao's work, they stated that the original HHS cannot distinguish the chatter from stable cutting, and the dispersion degree of the frequency spots in the HHS with WPT pre-denoising will be more and more concentrated as long as the chatter develops. However, from the results showed in Fig. 18, Cao's method seems to present no improvement in enhancing the performance of the standard HHT method. On the contrary, the distribution of the frequency spots is more dispersive after the WPT denoising applied, as showed in Fig. 18(b). However another HHS obtained from the [32] wavelet, which is the most powerful sub-signal, seems to reveal a better performance. This is different from the proposed [33] wavelet in the paper, which indicates that the proper choice of wavelet coefficients varies in different cases. The choice of mother wavelet is another problem which is still empirical and there are no efficient rules on how to appropriately make the selection. Different wavelets would lead to a discrepancy in performance [20]. This leads the wavelet transform to remain a popular manual method and cannot be widely used in online monitoring applications. Besides, the choice is hard to make even with a priori knowledge of the modes of the machine because of the non-determinacy and the time-varying characteristics of the cutting progress and the uncertain performance of different mother wavelet [38]. An appropriate combination of the decomposed signals other than a single selection may be more adaptive for online chatter detection which needs further verifications.

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

685

Fig. 20. Time complexity of the proposed method. (Fs¼ 20,480 Hz, NE¼ 100, NS¼ 10, NM ¼8, f c ¼ 41  100  10  8  20480FLOPS ¼ 6:717GFLOPS. Testing platform: Windows 7  64, Intel(R) Core(TM) i7-4790 CPU @ 3.60 GHz (4 C, HT, 3.6 GHz, 4  256 kB L2, 8MB L3, 100 MHz FSB), 16 g DDR3).

Fig. 19 shows the energy limited Hilbert Huang spectrum (EL-HHS) of the measured vibration obtained using the proposed method in the same case. The color definition is the same as in Fig. 6. A clear change can be found around 1.25 s and after that the frequency point centralized to a band around 3200 Hz. Chatter is underway and can be duly identified at 1.252 s. The proposed energy aggregation characteristic-based Hilbert Huang transform method can physically characterize the cutting process based on the energy distribution on frequency over time. The progress, which is self-adaptive, is more suitable for milling chatter detection. 3.6. Time complexity analysis It is critical to consider the time and space complexity for a real-time system. The responding speed mostly depends on the ratio of the processing time and the sampling time regardless of the system delay. For a chatter detection method, it is very important to acquire a reasonable response speed to protect the workpiece from the chatter damage as early as possible, especially in finish machining. The computer space needed is no longer a problem in nowadays. Next we will analyze the time complexity of the total chatter detection algorithm and the selection of the sampling parameters selection in detail. The main procedure of the EEMD part includes ensemble cycle, memory allocation, extreme identification and spline procedure. These parts have been discussed and summarized by Wang [39] as T EEMD ¼ NE  NS  ð41  NMÞ  n The selection part includes energy calculating, energy based descending sort and IMF selection, calculated as   T selection ¼ m  ðn þ n  1Þ þ n  log 2 n þ m ¼ log2 n þ2m n

ð16Þ

ð17Þ

and T HSA ¼ Oðn log nÞ T detection ¼ n þ n þ 2n þ 1 þ 2 ¼ 3n þ3

ð18Þ

As we can see that the coefficient NE  NS  ð41  NMÞ is usually quite large, resulting that the T EEMD is usually much larger than the others, we can obtain T  T EEMD ¼ NE  NS  ð41  NMÞ  n

ð19Þ

Considering the rolling window parameters, frame length ðt w ¼ n=FsÞ and frame shift (ts), using f donating the floating point operations per second of the processor, then the processing time of a signal segment in the rolling window can be expressed as t cpu ¼

T NE  NS  ð41  NMÞ  Fs ¼ tw f f

ð20Þ

686

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

Table 4 A list of the complexity parameter set. Parameter

Symbol Recommended

Notes

Sampling frequency

Fs

/

Frame length

tw

/

Fs should be selected according to the frequency characteristic following sampling theorem, 5–10 times the highest useful frequency Long enough to assure the signal quality and eliminate the effect of energy leakage and the average effect of Fourier transform. Appropriate to perform the function of the EEMD, and not to large which will weaken the noise effect. Author recommender log 2(n) is usually recommended, here we use a fixed number to stabilize the calculating time cost without any loss of the information. Time shift should be large than the time resolution to ensure the system instantaneity.

Number of ensemble NE

100

Sift number Number of IMFs

NS NM

10 8

Frame shift

ts

t s Z Rt ¼

fc tw fp

The actual time length of the update signal is ts, and for a real-time system, the processing time should be at most equal to the updating sampling time, that is t cpu r t s r t w

ð21Þ

when tw ¼ts, we can obtain the critical floating-point operations per second (FLOPS) fc, f c ¼ NE  NS  ð41  NMÞ  Fs

ð22Þ

For a certain processor with fp 4fc processing power, the highest time resolution Rt can be expressed as t s Z Rt ¼

fc tw fp

ð23Þ

The selection method of these parameters are summarized in Table 4. Numerical experiments are conducted to validate these relationships. The measured signal is extracted from No. 5 case which involves abundant types of signals. For every length of the rolling window, many different samples are used to test the CPU time running on a single core, and plot the curve of the mean value. The calculating result is shown in Fig. 20. The computing power of the tested CPU fp is 11.63GFLOPS per core obtained by aggregating arithmetic performance test using SiSoftware Sandra Professional Home 2009 and the criticality is 6.717GFLOPS. The gradient of the CPU time line is around 0.5345, close to the theoretic gradient f c =f p ¼ 6:71=11:63 ¼ 0:5769. When the gradient is smaller than 1, there will be a feasible domain, colored in orange, which is formed by the sampling time line in blue and the upper fluctuating CPU time bound in bright green. There is an unstable region when the length of the rolling window is too short, mainly caused by the neglect of those nonlinear terms. Meanwhile, the bottom bound of the shift time increases and the feasible region enlarges along with the increasing of the frame length, which is in agreement with the previous theoretic analysis. For a practical application, the feature map can be obtained using some benchmark tests. Then an appropriate sample length should be decided to ensure the signal quality and the efficiency of the signal processing method. After that, the frame shift can be determined according to the feature map. The improvement of the processor will always increase the resolution. As proved that the proposed algorithm is theoretically and practically feasible and realistic for the online chatter detection application.

4. Conclusions Focusing on the online chatter detection problem, this paper presents a comprehensive solution by proposing an energy aggregation characteristic-based Hilbert–Huang transform method. Firstly, the measured vibration signal is decomposed into a series of instinct mode functions (IMFs) using ensemble empirical mode decomposition. Then, a principle energy based selection is conducted on the IMFs according to the energy limitation coefficient to acquire the majority energy and feature IMFs. Subsequently, Hilbert spectral analysis is applied on the feature IMFs to calculate the time/frequency spectrum. Afterwards normalized energy ratio (NER) and coefficient of variation (CV) are calculated to quantize the spectrum characteristic for online detection system. Experiments are performed to validate the efficiency of the proposed method and numerical experiments are conducted to verify the instantaneity of the detection algorithm. The numerical and experimental results showed that: (a) As a violate resonance phenomenon, vibration energy, aggregating on a certain natural frequency of the cutting system, is a critical signal characteristic reflecting the chatter development; (b) The time–frequency spectrum obtained from the proposed energy aggregation characteristic-based Hilbert–Huang transform method can dynamically reveal the change from stable stage to chatter, and the transition stage can be

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

687

identified by the indicators at its early beginning so that the workpiece can be successfully protected from serious vibration damage; (c) The detection is based on the physical essence of chatter, which makes it easier to transplant to any other machining scenario without too much manual work compared with other statistical and fitting methods; (d) The proposed method is self-adaptive and robust in processing the unstable and nonlinear chatter signal without too much manual intervention; (e) The time complexity of the proposed method is linear, and the calculating quantity can be adjusted to match the CPU computing power with a proper parameter set. Future improvements of the method may account for a better decomposition method which can self-adaptively divide the signal into stricter narrow-band components with an appropriate time complexity.

Acknowledgments The authors would like to acknowledge financial support from the National Program on Key Basic Research Project (Grant no. 2013CB035805, 2012CB025903), National Natural Science Foundation Council of China (Grant no. 51125021).

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

G. Quintana, J. Ciurana, Y. Altintas, M. Weck, Chatter in machining processes: a review, Int. J. Mach. Tools Manuf. 51 (2011) 363–376. Y. Altintas, M. Weck, Chatter stability of metal cutting and grinding, CIRP Ann.-Manuf. Technol. 53 (2004) 619–642. Y. Altintas, Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design, Cambridge University Press, UK, 2012. F.W. Taylor, On the art of cutting metals, Trans. ASME 28 (1907). E. Budak, Y. Altintas, Analytical prediction of chatter stability in milling – Part 1: general formulation, J. Dyn. Syst. Meas. Control – Trans. ASME 120 (1998) 22–30. J. Tlusty, M. Polacek, The stability of machine tools against self excited vibrations in machining, ASME Int. Res. Prod. Eng. (1963) 465–474. J.V. Abellan-Nebot, F.R. Subirón, A review of machining monitoring systems based on artificial intelligence process models, Int. J. Adv. Manuf. Technol. 47 (2010) 237–257. E. Kuljanic, G. Totis, M. Sortino, Development of an intelligent multisensor chatter detection system in milling, Mech. Syst. Signal Process. 23 (2009) 1704–1718. T.L. Schmitz, K. Medicus, B. Dutterer, Exploring once-per-revolution audio signal variance as a chatter indicator, Mach. Sci. Technol. 6 (2002) 215–233. M. Lamraoui, M. Thomas, M. El Badaoui, E. Girardin, Indicators for monitoring chatter in milling based on instantaneous angular speeds, Mech. Syst. Signal Process. 44 (2014) 72–85. R. Teti, K. Jemielniak, G. O'Donnell, D. Dornfeld, Advanced monitoring of machining operations, CIRP Ann.-Manuf. Technol. 59 (2010) 717–739. H. Moradi, M.R. Movahhedy, G. Vossoughi, Dynamics of regenerative chatter and internal resonance in milling process with structural and cutting force nonlinearities, J. Sound Vib. 331 (2012) 3844–3865. Y. Altintas, M. Eynian, H. Onozuka, H.R. Cao, B. Li, Z.J. He, Identification of dynamic cutting force coefficients and chatter stability with process damping, CIRP Ann.-Manuf. Technol. 57 (2008) 371–374. L. Vela-Martinez, J.C. Jauregui-Correa, J. Alvarez-Ramirez, Characterization of machining chattering dynamics: An R/S scaling analysis approach, Int. J. Mach. Tools Manuf. 49 (2009) 832–842. E. Kuljanic, M. Sortino, G. Totis, Multisensor approaches for chatter detection in milling, J. Sound Vib. 312 (2008) 672–693. H.Q. Liu, Q.H. Chen, B. Li, X.Y. Mao, K.M. Mao, F.Y. Peng, On-line chatter detection using servo motor current signal in turning, Sci. China Technol. Sci. 54 (2011) 3119–3129. N.C. Tsai, D.C. Chen, R.M. Lee, Chatter prevention for milling process by acoustic signal feedback, Int. J. Adv. Manuf. Technol. 47 (2010) 1013–1021. T.L. Schmitz, Chatter recognition by a statistical evaluation of the synchronously sampled audio signal, J. Sound Vib. 262 (2003) 721–730. Z.H. Yao, D.Q. Mei, Z.C. Chen, On-line chatter detection and identification based on wavelet and support vector machine, J. Mater. Process. Technol. 210 (2010) 713–719. Z.K. Peng, F.L. Chu, Application of the wavelet transform in machine condition monitoring and fault diagnostics: a review with bibliography, Mech. Syst. Signal Process. 18 (2004) 199–221. N.E. Huang, Z. Shen, S.R. Long, M. Wu, H.H. Shih, Q.N. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. Lond. Ser. A – Math. Phys. Eng. Sci. 454 (1998) 903–995. Z.S. Yang, Z.H. Yu, C. Xie, Y.F. Huang, Application of Hilbert–Huang transform to acoustic emission signal for burn feature extraction in surface grinding process, Measurement 47 (2014) 14–21. N.E. Huang, Z.H. Wu, A review on Hilbert–Huang transform: method and its applications to geophysical studies, Rev. Geophys. 46 (2008). Z.K. Peng, P.W. Tse, F.L. Chu, A comparison study of improved Hilbert–Huang transform and wavelet transform: Application to fault diagnosis for rolling bearing, Mech. Syst. Signal Process. 19 (2005) 974–988. M. Datig, T. Schlurmann, Performance and limitations of the Hilbert–Huang transformation (HHT) with an application to irregular water waves, Ocean Eng. 31 (2004) 1783–1834. E.S. Carbajo, R.S. Carbajo, C. Mc Goldrick, B. Basu, ASDAH: An automated structural change detection algorithm based on the Hilbert–Huang transform, Mech. Syst. Signal Process. 47 (2014) 78–93. Y.G. Lei, J. Lin, Z.J. He, M.J. Zuo, A review on empirical mode decomposition in fault diagnosis of rotating machinery, Mech. Syst. Signal Process. 35 (2013) 108–126. B. Sick, On-line and indirect tool wear monitoring in turning with artificial neural networks: a review of more than a decade of research, Mech. Syst. Signal Process. 16 (2002) 487–546. C.L. Zhang, X. Yue, Y.T. Jiang, W. Zheng, A hybrid approach of ANN and HMM for cutting chatter monitoring, Adv. Mater. Res. 97 (2010) 3225–3232. M. Siddhpura, R. Paurobally, A review of chatter vibration research in turning, Int. J. Mach. Tools Manuf. 61 (2012) 27–47. N.E. Huang, S.S. Shen, Hilbert–Huang Transform and its Applications, World Scientific, USA, 2005. Z. Wu, N.E. Huang, Ensemble empirical mode decomposition: a noise-assisted data analysis method, Adv. Adapt. Data Anal. 1 (2009) 1–41. A.M. Al-Ghamd, D. Mba, A comparative experimental study on the use of acoustic emission and vibration analysis for bearing defect identification and estimation of defect size, Mech. Syst. Signal Process. 20 (2006) 1537–1571.

688

Y. Fu et al. / Mechanical Systems and Signal Processing 75 (2016) 668–688

[34] P. Paalanen, J.K. Kamarainen, J. Iloen, H. Kalviainen, Feature representation and discrimination based on Gaussian mixture model probability densities – practices and algorithms, Pattern Recogn. 39 (2006) 1346–1358. [35] J.A. Zhang, R.Q. Yan, R.X. Gao, Z.H. Feng, Performance enhancement of ensemble empirical mode decomposition, Mech. Syst. Signal Process. 24 (2010) 2104–2123. [36] J.-R. Yeh, J.-S. Shieh, N.E. Huang, Complementary ensemble empirical mode decomposition: a novel noise enhanced data analysis method, Adv. Adapt. Data Anal. 2 (2010) 135–156. [37] H.R. Cao, Y.G. Lei, Z.J. He, Chatter identification in end milling process using wavelet packets and Hilbert–Huang transform, Int. J. Mach. Tools Manuf. 69 (2013) 11–19. [38] K. Zhu, Y.S. Wong, G.S. Hong, Wavelet analysis of sensor signals for tool condition monitoring: a review and some new results, Int. J. Mach. Tools Manuf. 49 (2009) 537–553. [39] Y.H. Wang, C.H. Yeh, H.W.V. Young, K. Hu, M.T. Lo, On the computational complexity of the empirical mode decomposition algorithm, Physica A 400 (2014) 159–167. [40] N.E. Huang, M.-L.C. Wu, S.R. Long, S.S. Shen, W. Qu, P. Gloersen, K.L. Fan, A confidence limit for the empirical mode decomposition and Hilbert spectral analysis, Proc. R. Soc. Lond. A: Math., Phys. Eng. Sci., R. Soc. (2003) 2317–2345.