The Bivariate Empirical Mode Decomposition and Its Contribution to Grinding Chatter Detection

: Grinding chatter reduces the long-term reliability of grinding machines. Detecting the negative effects of chatter requires improved chatter detection techniques. The vibration signals collected from grinders are mainly nonstationary, nonlinear and multidimensional. Hence, bivariate empirical mode decomposition (BEMD) has been investigated as a multiple signal processing method. In this paper, a feature vector extraction method based on BEMD and Hilbert transform was applied to the problem of grinding chatter. The effectiveness of this method was tested and validated with a simulated chatter signal produced by a vibration signal generator. The extraction criterion of true intrinsic mode functions (IMFs) was also investigated, as well as a method for selecting the most ideal number of projection directions using the BEMD algorithm. Moreover, real-time variance and instantaneous energy were employed as chatter feature vectors for improving the prediction of chatter. Furthermore, the combination of BEMD and Hilbert transform was validated by experimental data collected from a computer numerical control (CNC) guideway grinder. The results reveal the good behavior of BEMD in terms of processing nonstationary and nonlinear signals, and indicating the synchronous characteristics of multiple signals. Extracted chatter feature vectors were demonstrated to be reliable predictors of early grinding chatter.


Introduction
Grinders, which are often high-precision, large-tonnage machines, are essential in the equipment manufacturing industry. Monitoring their condition and diagnosing their faults ensures their safe operation, which has great practical and economic value [1,2]. Most notably, grinders can enter a "chatter" state during machining operations, which causes a series of negative effects [3] such as machined surface undulations, reduced tool life, poor surface finish, noise, and reduced productivity.
Friction, regeneration and mode coupling are three well-known mechanisms that cause grinding chatter. Typically, regeneration chatter is the most common type of chatter encountered in machining processes. It occurs when oscillations in the grinder cause undulations in the machined surface that are regenerated with subsequent passes of the grinder. As numerous experiments have shown, the amplitude of vibration signals increases substantially when a stable grinding state turns into a chatter state [4]. Additionally, there is a transitional phase in this process, which contains considerable information about grinding status [5]. Thus, we can discover the mechanism of grinding chatter and predict it by analyzing transition-phase signals. It is worth noting that the transitional phase occurs within a very short time and is easily affected by many random factors. In order to alleviate the problem of chatter, reliable chatter detection and identification methods are essential so that effective chatter suppression can be applied in a timely manner [6].
Vibration detection in grinding is usually performed by using accelerometers, acoustic emission (AE) sensors or dynamometers. A considerable amount of research work describes these different methods. For example, Tansel et al. used acceleration signals from a lathe tailstock to detect chatter during turning operations, and used an s-transformation to extract a damping index as a descriptive feature of chatter [7]. Yao et al. proposed two-dimensional feature vectors for chatter detection based on the standard deviation of wavelet transform of drilling machining, which allowed rapid chatter identification [8]. Furthermore, Devillez et al. equipped a lathe with eddy current sensors to monitor vibration. These were processed by a fuzzy method to classify the tests [9]. Moreover, many researchers have predicted cutting chatter by using dynamic forces. Gradisek et al. employed a coarse-grained entropy rate as a chatter index in a cutting process, which exhibited a drastic drop at the onset of chatter [10]. Nari et al. presented the usage of permutation entropy, a computer-based, rapid measurement of chatter onset using sound signals recorded by a unidirectional microphone [11]. Yang et al. proposed a method using feature vectors incorporating signal variance and a one-step self-correlation function as a chatter indicator for detecting cutting state [12]. Acoustic emission sensors were used by Li et al. to detect chatter; the specific peculiarities of acoustic emission sensors can effectively reject external disturbances [13].
In summary, wavelet transform [14,15], s-transformation [16], singular value decomposition (SVD) [17], and artificial intelligence [18], are the most popular methods used to process vibration signals. These analysis methods are mostly based on the theory of Fourier transform. However, grinder vibration signals are mostly nonstationary and nonlinear, such that these conventional methods cannot effectively extract signal features [19]. Subsequently, empirical mode decomposition (EMD), proposed by Huang et al., has been demonstrated as an appropriate approach for overcoming the shortcomings of conventional techniques [20]. The EMD approach can decompose signals into several intrinsic mode functions (IMFs) which distribute from high frequency to low frequency. All IMFs derived from grinding signals are characterized as zero-mean oscillations. However, in practice, grinder vibration signals are usually multidimensional signals, yet current methods are only useful with one-dimensional, real-value time series [21]. They are unfavorable for executing information fusion functions and accurately reflecting real-world situations when processing multiple signals.
In order to overcome the aforementioned problems, Rilling et al. proposed an extension to bivariate time series, which generalizes the rationale underlying EMD to the bivariate framework, namely, bivariate empirical mode decomposition (BEMD) [22]. Since being proposed, it has been applied to image processing such as image segmentation [23], image fusion [24], image compression [25] and image watermarking [26]. At the third session of the HHT International Conference, BEMD was reported to have been successfully applied to wind turbine condition monitoring [27]. It has been shown that the BEMD technique has inherited all the merits of the one-dimensional EMD method. However, in contrast to EMD, BEMD aims at decomposing a complex-valued signal into a collection of zero-mean rotating components, which represent fast to slow rotations in the signal. Another greater advantage is that BEMD is more powerful in extracting weak features from nonstable and nonlinear signals [27]. The authors of this paper have made a comparison between the use of EMD and BEMD for extracting features from grinding chatter signals, and demonstrated the superior performance of BEMD [28]. These findings can be summarized as follows: (1) EMD is initially applied to a one-dimensional signal and extracts zero-mean oscillating components, whereas BEMD is applied to a bivariate signal and extracts zero-mean rotating components; (2) BEMD has computational efficiency due to the simultaneous processing of complex-value signals, and it only computes the upper envelope using the maximum points. Meanwhile, EMD can only decompose signals one by one and obtains both upper and lower envelopes by connecting the extreme points; (3) the number of IMFs derived from signals by EMD are unequal. They cannot reveal any synchronous Appl. Sci. 2017, 7, 145 3 of 17 characteristics and phase shifting, nor can EMD extract an information fusion function. In contrast, the number of IMFs derived by BEMD is normally the same, and BEMD can extract an information fusion function and preserve phase differences; (4) BEMD facilitates the establishment of purified shaft vibration orbits and fully guarantees the correctness of results, which EMD cannot.
The objective of this paper is to assess BEMD's performance in processing grinding chatter signals. The remainder of the paper is organized as follows. In Section 2, the algorithms of BEMD and Hilbert transform are reviewed, as well are the extraction criteria of true IMFs. Moreover, real-time variance and instantaneous energy are presented as feature vectors for grinding chatter. In Section 3, a simulated chatter signal is generated from a chatter signal generator (Simulink, MATLAB, MathWorks Company, Natick, MA, USA), then processed by BEMD and Hilbert transform. Moreover, a selection method for the most ideal number of projection directions is developed. In Section 4, the benefits of combining BEMD and Hilbert transform are experimentally validated by processing chatter signals collected from a real grinder (KD4020X16). Section 5 concludes the paper.

The BEMD Algorithm
The main concept of BEMD is to extract a series of IMFs from multiple signals. This method assumes that multiple signals can be described as the sum of a fast rotation component superimposed on a slower rotation component. Each of the IMFs must satisfy the following conditions: 1.
The number of extrema and the number of zero crossings must be equal, or differ at most by one.

2.
The mean value of the envelope at any point defined by the local maximum points and the envelope defined by the local minima must be zero.
Moreover, the BEMD process requires projecting a multiple signal on a set of directions and then applying the sifting process of standard EMD to the projected components. This gives a bivariate signal, s(t) = x(t) + iy(t), and a set of projection directions, ϕ m = 2mπ/N, 1 ≤ m ≤ N. The fundamental process of BEMD can be depicted as follows: 1.
Project the signal s(t) on directions ϕ m : 2. Extract all partial maximum points of p ϕ m (t):{(t m j , p m j )}, where j indicates the number of individual partial maxima.

3.
Interpolate the set of points {(t m j , e iϕ m p m j )} by cubic spline interpolation to obtain the partial envelope curve in direction ϕ m named e ϕ m (t).

4.
Repeat steps 1-3 until the envelope curves in all N projections are obtained.

5.
Compute the mean of all envelope curves: 6. Subtract the mean m(t) from s(t) to obtain g(t):

7.
Test if g(t) is an IMF, if not, replace s(t) with g(t) and repeat the procedure from Step 1 until g(t) is an IMF. If yes, record the obtained IMF and remove it from g(t), i.e., c 1 (t) = g(t), r 1 (t) = s(t) − c 1 (t).

8.
Consider r 1 (t) as the original signal and repeat the above procedure until the second IMF, c 2 (t), is obtained, and the residual component r 2 (t) = r 1 (t) − c 2 (t). 9.
Iterate the previous process until all IMFs are obtained. Referring to the sifting process of BEMD, the term s(t) can be expressed by the procedure: where h k (t) denotes the kth complex-valued IMF and r n (t) denotes a non-zero mean low-degree polynomial residue. Generally, the residue is mostly recognized as the vibration trend of the signal.

Extraction Criterion of True Intrinsic Mode Functions (IMFs)
The majority of vibration signals will be decomposed into excessive IMFs, which are described as spurious vibratory components that directly affect researchers' capacity to analyze signal features. These spurious IMFs also result in enormous interferences, reflecting the peculiarities of the chatter and allowing diagnosis of the mechanical faults of the grinder. One of the major deficiencies is that the estimation of the authenticity of IMF heavily depends on the user's experience. Thus, a reliable method is required to identify and remove spurious components, which is of great importance in extracting the actual vibration mode and the corresponding features of the time-frequency domain.
Spurious IMFs may be caused by several factors: (1) IMFs defined by BEMD algorithms are based on numerical analysis only and do not consider the vibration characteristics of system in the physical sense; (2) the stopping criterion of sifting results in a phenomenon of excessive decomposition and end effects are not well processed; (3) vibration signals mixed with noise or pulse interference will produce a series of spurious components with high frequency in the decomposition process. Since IMFs are recognized as the orthogonal expression of the signals, true IMFs have a higher correlation with original signals than spurious IMFs. Mathematically, a correlation coefficient is used as a statistical indicator to reflect the degree of correlation between variables. It is based on the deviation of two variables and their respective average mean, computed by the method of product comments, as shown in Equation (5).
where Cov(X, Y) denotes the covariance between the original signal and the IMF, E(X) and D(X) denote the expectation and variance of the original signal, respectively, and E(Y) and D(Y) denote the expectation and variance of the kth IMF, respectively. Therefore, it is reasonable to employ the correlation coefficient as an index to exactly distinguish the spurious components from true IMFs and then classify them as a part of residual [16,29]. In view of this, if the correlation coefficient of each IMF, ρ k (k = 1, 2, ... , n), is computed, the true IMFs can be extracted by the following method (Table 1).
Else, Estimate kth IMF, and r n = r n + c k .
The variable λ in Table 1 is a fixed threshold that is generally adopted as a ratio of the maximum correlation coefficient, wherein η is a ratio coefficient larger than 1.

Investigation of Projection Direction Number
Remember that the number of projection directions must be selected in advance because the BEMD algorithm is sensitive to it. Moreover, a large number of directions may be interesting insofar that it reduces the dependence of the final decomposition with respect to rotations of the spatial coordinates. However, excessive directions will reduce the computational efficiency of BEMD. Thus, projection direction number is further investigated by assessing the indicators of computational time and stability degree (SD) which change with the direction number.
The stability degree can be defined as, where X = {X 1 , X 2 , X 3 , ..., X N } T (X j ∈ R, j = 1, 2, ..., N, T is a transpose operator), N represents the sampling points, and µ represents the mean of the signal. Moreover, the SD is an ideal standard for expressing the fluctuation of signal amplitude. When there is a transient mutation in a signal, the stability degree will significantly improve. The purpose of this paper is to find out the most ideal number of projection directions, so that the precise signal features can be obtained as soon as possible.

Hilbert Transform
Through the BEMD method, the chatter signal can be decomposed into a number of complex-valued IMF components. Thereby, the Hilbert transform can be used to obtain the Hilbert and marginal spectrum distributions [30]. The Hilbert transform is defined as, where P is the Cauchy principal value. Define X(t) and Y(t) from c(t) as follows: X k (t) = Re(c k (t)) + ix k (t) and Y k (t) = Im(c k (t)) + iy k (t) (9) so that the amplitude and phase functions can be defined: where θ k (t) represents the phase of the real-part IMF, and ϕ k (t) represents the phase of the imaginary-part IMF. Then, the instantaneous frequency can be expressed as follows: Thus, the Hilbert spectrum can be obtained: Appl. Sci. 2017, 7, 145 6 of 17 The Hilbert spectrum represents the signals as a function that combines time and frequency. If H(ω,δ,t) is integrated over time, the corresponding marginal spectrum is obtained: The marginal spectrum expresses the amplitude of each frequency in space, and represents the accumulated amplitude in a statistical sense.
Although the Hilbert transform can analyze the instantaneous frequency, its applicability is restricted by the form of the original signal. For poorly formed signals, such as very irregular signals, the results produced by Hilbert transform differ significantly from the actual physical status. However, after BEMD decomposition, the IMFs have superior characteristics with a small bandwidth, reasonable symmetry, and smooth envelope changes, which are suitable for the Hilbert transform [31]. Thus, the BEMD method and Hilbert transform are highly connected.

Feature Vectors for Grinding Chatter-Real-Time Variance and Instantaneous Energy
Another major objective of this paper is to accurately detect and identify the onset of grinding chatter before chatter makes damage the workpiece, so that suitable techniques for chatter suppression can be employed as fast as possible. According to the experimental results, the amplitude of chatter signals exhibits a substantial expansion when a stable grinding state changes into a chatter state. Consequently, early grinding chatter can be detected by comparing changes in the time-domain statistical parameters of the signal mathematically. Therefore, appropriate feature vectors should be defined where real-time variance and instantaneous energy are recognized ideal feature vectors.
The real-time variance (R v ) indicates the degree of deviation from the mean chatter signal, which in a sense reflects the trend of signal oscillation. The real-time variance can be described as: where X = {X 1 , X 2 , X 3 , ..., X N } T (X j ∈ R, j = 1, 2, ..., N, T is a transpose operator), N represents the sampling points, and X represents the mean of the signal. In practice, any internal fluctuation may result in a vibration of instantaneous energy, which means that the change in instantaneous energy has a direct relation to abnormal system operation. Thus, according to Equations (10) and (11), instantaneous energy can be defined as follows: where α k (t) represents the instantaneous amplitude of the kth IMF. By these definitions, the R v and E of each IMF can be calculated out in real time.

Application of a Feature Vector Extraction Method for Simulating Chatter Signals
In order to evaluate whether BEMD and Hilbert transform can be a reliable solution for detecting chatter, an experimental program based on BEMD and Hilbert transform was designed ( Figure 1). The process consists of the following steps: (1) signal acquisition and pre-processing, such as signal denoising and reconstruction; (2) applying the BEMD and extraction criterion of true IMFs to obtain all true IMFs; (3) applying the Hilbert transform to true IMFs to obtain Hilbert and marginal spectrum distributions; (4) calculating the real-time variance and instantaneous energy of IMFs, and then superposing and normalizing the feature vectors, respectively; and (5) predicting grinding chatter by contrasting the changes of the feature vectors.

The Chatter Vibration Signal Generator
For the sake of demonstration, a BEMD is adopted to process a simulated chatter signal. According to the mechanism of the chatter and the characteristics of frequency and time domains in the grinding process, a chatter vibration signal generator was established through the Simulink module of MATLAB (2010a, MathWorks Company, Natick, MA, USA). The sampling frequency was set as 800 Hz, as shown in Figure 2. The principle of this chatter generation system can be described as follows. A vibration signal is generated by two sine wave generators and one white noise generator. The frequency of the first sine signal consists of 50 rad/s and 120 rad/s, and the other one consists of 40 rad/s and 80 rad/s. The phase of the second sine signal is shifted by 4 rad. Thus, the two sine signals and white noise can be expressed as:

The Chatter Vibration Signal Generator
For the sake of demonstration, a BEMD is adopted to process a simulated chatter signal. According to the mechanism of the chatter and the characteristics of frequency and time domains in the grinding process, a chatter vibration signal generator was established through the Simulink module of MATLAB (2010a, MathWorks Company, Natick, MA, USA). The sampling frequency was set as 800 Hz, as shown in Figure 2.

The Chatter Vibration Signal Generator
For the sake of demonstration, a BEMD is adopted to process a simulated chatter signal. According to the mechanism of the chatter and the characteristics of frequency and time domains in the grinding process, a chatter vibration signal generator was established through the Simulink module of MATLAB (2010a, MathWorks Company, Natick, MA, USA). The sampling frequency was set as 800 Hz, as shown in Figure 2. The principle of this chatter generation system can be described as follows. A vibration signal is generated by two sine wave generators and one white noise generator. The frequency of the first sine signal consists of 50 rad/s and 120 rad/s, and the other one consists of 40 rad/s and 80 rad/s. The phase of the second sine signal is shifted by 4 rad. Thus, the two sine signals and white noise can be expressed as: The principle of this chatter generation system can be described as follows. A vibration signal is generated by two sine wave generators and one white noise generator. The frequency of the first sine signal consists of 50 rad/s and 120 rad/s, and the other one consists of 40 rad/s and 80 rad/s. The phase of the second sine signal is shifted by 4 rad. Thus, the two sine signals and white noise can be expressed as: x(t) = 4 sin(120t) + 4 sin(50t) + 0.3 and y(t) = 2 sin(80t + 4) This system relies on the clock module to control current simulation time. When t ≤ 1 s, the output is a harmonic vibration signal, in order to simulate stable grinding. When 1 < t ≤ 1.6 s, the output is a harmonic vibration signal multiplied with a slant signal, so as to simulate grinding chatter. Furthermore, when 1.6 < t ≤ 2.5 s, the output is the original harmonic vibration signal multiplied by a gain coefficient, in order to simulate stable chattering. For the sake of simplicity, this chatter generation system can be expressed in a mathematical form: The chatter signal generated by this system is shown in Figure 3. It is clearly seen that signal amplitude is small during stable grinding, but significantly increases after 1 s. Then after 1.6 s, the amplitude becomes steady when the grinder reaches stable chatter. Hence the trend and distribution are similar to experimental chatter images in [7,32], meaning that this system simulates chatter signals well.
This system relies on the clock module to control current simulation time. When t ≤ 1 s, the output is a harmonic vibration signal, in order to simulate stable grinding. When 1 < t ≤ 1.6 s, the output is a harmonic vibration signal multiplied with a slant signal, so as to simulate grinding chatter. Furthermore, when 1.6 < t ≤ 2.5 s, the output is the original harmonic vibration signal multiplied by a gain coefficient, in order to simulate stable chattering. For the sake of simplicity, this chatter generation system can be expressed in a mathematical form: The chatter signal generated by this system is shown in Figure 3. It is clearly seen that signal amplitude is small during stable grinding, but significantly increases after 1 s. Then after 1.6 s, the amplitude becomes steady when the grinder reaches stable chatter. Hence the trend and distribution are similar to experimental chatter images in [7,32], meaning that this system simulates chatter signals well.

The Application of BEMD to Simulated Chatter Signal
Decomposing this chatter signal by BEMD sets 64 projection directions and 10 iterations. The extraction criterion of true IMFs introduced in Section 2.1.2 can be applied to the decomposed IMFs, and the correlation coefficients of each IMF are shown in Table 2. Compared with the data in Table 2, only the first two IMFs, which have relatively high correlation with the original signals, should be reserved. The other three IMFs have to be removed and classified as a part of the residual. The true IMFs derived from the simulated chatter signal are shown in Figure 4.
In Figure 4, it is clearly seen that the phase shifting and synchronization information between the real part and imaginary components are well preserved and detected by BEMD. The true IMFs (i.e., IMF1 and IMF2) also clearly reveal the amplitude and frequency components of the signals which were initially set in the chatter signal generator.

The Application of BEMD to Simulated Chatter Signal
Decomposing this chatter signal by BEMD sets 64 projection directions and 10 iterations. The extraction criterion of true IMFs introduced in Section 2.1.2 can be applied to the decomposed IMFs, and the correlation coefficients of each IMF are shown in Table 2. Compared with the data in Table 2, only the first two IMFs, which have relatively high correlation with the original signals, should be reserved. The other three IMFs have to be removed and classified as a part of the residual. The true IMFs derived from the simulated chatter signal are shown in Figure 4.
In Figure 4, it is clearly seen that the phase shifting and synchronization information between the real part and imaginary components are well preserved and detected by BEMD. The true IMFs (i.e., IMF1 and IMF2) also clearly reveal the amplitude and frequency components of the signals which were initially set in the chatter signal generator.

The Application of Hilbert Transform to Simulated Chatter Signals
By employing the Hilbert transform to true IMFs, the Hilbert and marginal spectra in  In Figure 5, it can be clearly seen that the main amplitudes of the real IMFs are distributed around 8 and 19 Hz, and the imaginary IMFs are around 6.4 and 12.8 Hz. However, note that after 1 s, the amplitude increases and reaches a maximum after 1.6 s due to the onset of chatter. That is to say, chatter results in increased amplitude. In Figure 6, the frequency components are accurately revealed; the peak amplitude values correspond with the frequency components (40, 50, 80 and 120 rad/s, i.e., 6.4, 8, 12.8 and 19 Hz) as set in the chatter generation system. It is clear that the amplitude at low frequency is larger than at high frequency, because the main frequency of vibration becomes low at the onset of grinding chatter. Consequently, the BEMD and Hilbert transform can be successfully combined to extract the actual features from the chatter signals.

The Application of Hilbert Transform to Simulated Chatter Signals
By employing the Hilbert transform to true IMFs, the Hilbert and marginal spectra in

The Application of Hilbert Transform to Simulated Chatter Signals
By employing the Hilbert transform to true IMFs, the Hilbert and marginal spectra in Figures 5 and 6 are obtained. In Figure 5, it can be clearly seen that the main amplitudes of the real IMFs are distributed around 8 and 19 Hz, and the imaginary IMFs are around 6.4 and 12.8 Hz. However, note that after 1 s, the amplitude increases and reaches a maximum after 1.6 s due to the onset of chatter. That is to say, chatter results in increased amplitude. In Figure 6, the frequency components are accurately revealed; the peak amplitude values correspond with the frequency components (40, 50, 80 and 120 rad/s, i.e., 6.4, 8, 12.8 and 19 Hz) as set in the chatter generation system. It is clear that the amplitude at low frequency is larger than at high frequency, because the main frequency of vibration becomes low at the onset of grinding chatter. Consequently, the BEMD and Hilbert transform can be successfully combined to extract the actual features from the chatter signals.

Selection of the Ideal Number of Projection Directions
According to the simulated signal, variation in the correlation coefficient, frequency and calculation time, along with the change in direction number, can be calculated as shown in Figure 7. The trends of correlation coefficient and frequency variation tend to be stable when the maximum projection number is set to 32. Higher numbers are not discussed in this paper. In Figure 5, it can be clearly seen that the main amplitudes of the real IMFs are distributed around 8 and 19 Hz, and the imaginary IMFs are around 6.4 and 12.8 Hz. However, note that after 1 s, the amplitude increases and reaches a maximum after 1.6 s due to the onset of chatter. That is to say, chatter results in increased amplitude. In Figure 6, the frequency components are accurately revealed; the peak amplitude values correspond with the frequency components (40, 50, 80 and 120 rad/s, i.e., 6.4, 8, 12.8 and 19 Hz) as set in the chatter generation system. It is clear that the amplitude at low frequency is larger than at high frequency, because the main frequency of vibration becomes low at the onset of grinding chatter. Consequently, the BEMD and Hilbert transform can be successfully combined to extract the actual features from the chatter signals.

Selection of the Ideal Number of Projection Directions
According to the simulated signal, variation in the correlation coefficient, frequency and calculation time, along with the change in direction number, can be calculated as shown in Figure 7. The trends of correlation coefficient and frequency variation tend to be stable when the maximum projection number is set to 32. Higher numbers are not discussed in this paper.

Selection of the Ideal Number of Projection Directions
According to the simulated signal, variation in the correlation coefficient, frequency and calculation time, along with the change in direction number, can be calculated as shown in Figure 7. The trends of correlation coefficient and frequency variation tend to be stable when the maximum projection number is set to 32. Higher numbers are not discussed in this paper. The degree of stability of the correlation coefficient and frequency for the true IMFs is shown in Figure 8.
In light of the indicators from the above figures, it is clearly seen that the correlation coefficient and frequency components of the signals are almost stable when the direction number is greater than 16. On the other hand, the computational efficiency of the BEMD algorithm gradually reduces with increasing projection number. In general, the ideal number of projection directions is 16, which provides a basis for subsequent BEMD operations. The degree of stability of the correlation coefficient and frequency for the true IMFs is shown in Figure 8.

Chatter Feature Vector Extraction Based on BEMD and Hilbert Transform
By BEMD decomposition and extraction of true IMFs, the chatter signal is decomposed into a series of true IMFs and a residue. According to Equations (15) and (16), the real-time variance and In light of the indicators from the above figures, it is clearly seen that the correlation coefficient and frequency components of the signals are almost stable when the direction number is greater than 16. On the other hand, the computational efficiency of the BEMD algorithm gradually reduces with increasing projection number. In general, the ideal number of projection directions is 16, which provides a basis for subsequent BEMD operations.

Chatter Feature Vector Extraction Based on BEMD and Hilbert Transform
By BEMD decomposition and extraction of true IMFs, the chatter signal is decomposed into a series of true IMFs and a residue. According to Equations (15) and (16), the real-time variance and instantaneous energy of each true IMF can be calculated and then superposed, respectively, to construct synthetic chatter feature vectors for chatter detection and identification. For the sake of comparison and analysis, the chatter feature vectors should be normalized to keep their value between 0 and 1 [33]. The superposed and normalized real-time variance and instantaneous energy of IMFs are shown in Figure 9.

Chatter Feature Vector Extraction Based on BEMD and Hilbert Transform
By BEMD decomposition and extraction of true IMFs, the chatter signal is decomposed into a series of true IMFs and a residue. According to Equations (15) and (16), the real-time variance and instantaneous energy of each true IMF can be calculated and then superposed, respectively, to construct synthetic chatter feature vectors for chatter detection and identification. For the sake of comparison and analysis, the chatter feature vectors should be normalized to keep their value between 0 and 1 [33]. The superposed and normalized real-time variance and instantaneous energy of IMFs are shown in Figure 9. In order to achieve good recognition in the above figure, the maximum and mean values of realtime variance and instantaneous energy in various grinding states can be determined by removing the interference of endpoints ( Figure 10). In order to achieve good recognition in the above figure, the maximum and mean values of real-time variance and instantaneous energy in various grinding states can be determined by removing the interference of endpoints ( Figure 10). From Figures 9 and 10, it can be clearly seen that real-time variance and instantaneous energy in various grinding states exhibit different behaviors. The amplitude of these feature vectors is almost steady in the stable grinding stage, and the maximum and mean values are both very close. While the amplitude drastically increases once the signal enters the transition stage, the mean value of two feature vectors increases 1.4 and 2.5 times respectively, and the maximum increases 1.8 and 3.7 times. From Figures 9 and 10, it can be clearly seen that real-time variance and instantaneous energy in various grinding states exhibit different behaviors. The amplitude of these feature vectors is almost steady in the stable grinding stage, and the maximum and mean values are both very close. While the amplitude drastically increases once the signal enters the transition stage, the mean value of two feature vectors increases 1.4 and 2.5 times respectively, and the maximum increases 1.8 and 3.7 times. Moreover, when the signal is in a chatter stage, the instantaneous energy fluctuates within a certain range that the mean value and maximum are close. However, as real-time variance increases continuously with time, the mean and maximum increase 1.9 and 1.63 times when compared to the transition stage. Though it is difficult to distinguish the transition and chatter states exactly, it is feasible to clearly determine the onset of grinding chatter. Therefore, these two feature vectors can be applied as early chatter detectors.

Grinding Experiments
In order to further validate the feasibility of BEMD for grinding chatter detection, an experimental method was designed to acquire various vibration signals from different grinding parameters for a CNC guideway grinder (KD4020X16, Hangzhou Hangji Machine Tool Co., Ltd., Hangzhou, China). An integral electronic piezoelectric (IEPE) accelerometer sensor with a supporting dynamic signal test and analysis system (TST5912) were used to collect grinding vibration signals ( Figure 11). From Figures 9 and 10, it can be clearly seen that real-time variance and instantaneous energy in various grinding states exhibit different behaviors. The amplitude of these feature vectors is almost steady in the stable grinding stage, and the maximum and mean values are both very close. While the amplitude drastically increases once the signal enters the transition stage, the mean value of two feature vectors increases 1.4 and 2.5 times respectively, and the maximum increases 1.8 and 3.7 times. Moreover, when the signal is in a chatter stage, the instantaneous energy fluctuates within a certain range that the mean value and maximum are close. However, as real-time variance increases continuously with time, the mean and maximum increase 1.9 and 1.63 times when compared to the transition stage. Though it is difficult to distinguish the transition and chatter states exactly, it is feasible to clearly determine the onset of grinding chatter. Therefore, these two feature vectors can be applied as early chatter detectors.

Grinding Experiments
In order to further validate the feasibility of BEMD for grinding chatter detection, an experimental method was designed to acquire various vibration signals from different grinding parameters for a CNC guideway grinder (KD4020X16, Hangzhou Hangji Machine Tool Co., Ltd., Hangzhou, China). An integral electronic piezoelectric (IEPE) accelerometer sensor with a supporting dynamic signal test and analysis system (TST5912) were used to collect grinding vibration signals ( Figure 11). In practice, the grinder is more sensitive to rotational speed, feeding speed and grinding depth of the grinding wheel, which actually contribute to imbalances in grinding vibration. Thus, this experiment was carried out with the following steps:


Keeping the workpiece feeding speed and grinding wheel depth constant while gradually increasing the rotational speed; Stable Transition Chatter Stable Transition Chatter Figure 11. (a) CNC guideway grinder KD4020X16; (b) TST5912 analysis system.
In practice, the grinder is more sensitive to rotational speed, feeding speed and grinding depth of the grinding wheel, which actually contribute to imbalances in grinding vibration. Thus, this experiment was carried out with the following steps:

•
Keeping the workpiece feeding speed and grinding wheel depth constant while gradually increasing the rotational speed; • Keeping the feeding speed and rotational speed constant while gradually increasing the grinding depth; • Keeping the rotational speed and depth of grinding constant while gradually increasing the feeding speed.
Thus, the parameters are listed in Table 3. According to the experimental conditions shown in Table 3, eight piezoelectric acceleration sensors were used to test vibration acceleration signals from the grinding wheel spindle, motor and machine column in various directions [34]. The position of the sensors and their corresponding sensitivities are given in Table 4.
As per the above steps, 80 groups of grinding signals were obtained, of which 45 were in a stable grinding state, while 35 were in a chatter state. Based upon practical experience and considerable experimental data, it has been demonstrated that the chatter of wheel spindles occurs more in the xand z-directions than in the y-direction. Thus, for convenience of this presentation, this paper selects parts of xand z-direction chatter data as an original complex-valued signal. However, in the signal acquisition stage, the disturbing noise will inevitably be mixed in grinding vibration signal duo to the influence of many random factors, directly affecting detection accuracy of grinding chatter. It is reasonable to carry out pre-processing to eliminate noise for signals. Moreover, the wavelet de-noising technique, which has the characteristics of time domain localization, multi-resolution and flexibility of selecting wavelet basis, has been widely employed in signal processing [35,36]. Thus, in this paper, the noise is removed from the original complex-valued signal based on wavelet de-noising. Meanwhile, we resample the complex-valued signal based on Nyquist sampling theory [37]. The new, reconstructed complex-valued signal is shown in Figure 12.
From Figure 12, it is found that grinding chatter emerges at about 5-14 s and the transitional phase lasts almost 9 s. It is clearly seen that the amplitude of the vibration signal rapidly expands when the grinder starts chattering, then the amplitude becomes steady when the grinder enters stable chatter. However, the signal vibrates more markedly compared with the stable grinding state. Moreover, note that the abnormal signals derived from the joints of the workpiece will not greatly affect the overall signal analysis.


Keeping the feeding speed and rotational speed constant while gradually increasing the grinding depth;  Keeping the rotational speed and depth of grinding constant while gradually increasing the feeding speed.
Thus, the parameters are listed in Table 3. According to the experimental conditions shown in Table 3, eight piezoelectric acceleration sensors were used to test vibration acceleration signals from the grinding wheel spindle, motor and machine column in various directions [34]. The position of the sensors and their corresponding sensitivities are given in Table 4. Position As per the above steps, 80 groups of grinding signals were obtained, of which 45 were in a stable grinding state, while 35 were in a chatter state. Based upon practical experience and considerable experimental data, it has been demonstrated that the chatter of wheel spindles occurs more in the x-and z-directions than in the y-direction. Thus, for convenience of this presentation, this paper selects parts of x-and z-direction chatter data as an original complex-valued signal. However, in the signal acquisition stage, the disturbing noise will inevitably be mixed in grinding vibration signal duo to the influence of many random factors, directly affecting detection accuracy of grinding chatter. It is reasonable to carry out pre-processing to eliminate noise for signals. Moreover, the wavelet de-noising technique, which has the characteristics of time domain localization, multi-resolution and flexibility of selecting wavelet basis, has been widely employed in signal processing [35,36]. Thus, in this paper, the noise is removed from the original complex-valued signal based on wavelet de-noising. Meanwhile, we resample the complex-valued signal based on Nyquist sampling theory [37]. The new, reconstructed complex-valued signal is shown in Figure 12. From Figure 12, it is found that grinding chatter emerges at about 5-14 s and the transitional phase lasts almost 9 s. It is clearly seen that the amplitude of the vibration signal rapidly expands when the grinder starts chattering, then the amplitude becomes steady when the grinder enters stable chatter. However, the signal vibrates more markedly compared with the stable grinding state.

The Application of BEMD and Hilbert Transform to Experimental Chatter Signals
As previously mentioned, BEMD demonstrates its powerful capability in terms of processing the two-dimensional simulation chatter signal and extracting true IMFs from it. Moreover, the simulation results well illustrate the phase information and synchronization between real and imaginary parts of the IMFs, which is of significant benefit to applying the BEMD method to real-world chatter. Thus, the BEMD is applicable for decomposing the experimental signal above, which sets 16 projection directions (according to the selection of the ideal projection direction number, Section 3.4) and sets 10 iterations. Then, the extraction criterion based on the correlation coefficient can be applied to the extracted IMFs. The correlation coefficients of each IMF are shown in Table 5. From Table 5, it is clearly seen that the first three IMFs, which have relatively high correlation with the original signals, should be reserved. The other four IMFs have to be removed and classified as a part of the residual. So, the first three IMFs can be recognized as the true IMFs, which contain the main frequency components of the signal. The eventual decomposition result is shown in Figure 13.
imaginary parts of the IMFs, which is of significant benefit to applying the BEMD method to real-world chatter. Thus, the BEMD is applicable for decomposing the experimental signal above, which sets 16 projection directions (according to the selection of the ideal projection direction number, Section 3.4) and sets 10 iterations. Then, the extraction criterion based on the correlation coefficient can be applied to the extracted IMFs. The correlation coefficients of each IMF are shown in Table 5.  Table 5, it is clearly seen that the first three IMFs, which have relatively high correlation with the original signals, should be reserved. The other four IMFs have to be removed and classified as a part of the residual. So, the first three IMFs can be recognized as the true IMFs, which contain the main frequency components of the signal. The eventual decomposition result is shown in Figure 13. Then, the Hilbert transform is performed on each true IMF to obtain the marginal spectrum, as shown in Figure 14. Then, the Hilbert transform is performed on each true IMF to obtain the marginal spectrum, as shown in Figure 14. In Figure 13, the marginal spectrum of both x-and z-directions shows the same frequency components, of about 300, 580, and 1200 Hz. Yet the z-direction IMFs have larger amplitudes relative to the x-direction, as signals tend to vibrate more in the z-direction according to previous studies.

Chatter Feature Vector Extraction from Experimental Chatter Signals
Based on the same principle and method in Section 3.5, the feature vectors can also be extracted from the experimental true IMFs. The superimposed and normalized real-time variance and instantaneous energy vectors with time are shown in Figure 15. In Figure 13, the marginal spectrum of both xand z-directions shows the same frequency components, of about 300, 580, and 1200 Hz. Yet the z-direction IMFs have larger amplitudes relative to the x-direction, as signals tend to vibrate more in the z-direction according to previous studies.

Chatter Feature Vector Extraction from Experimental Chatter Signals
Based on the same principle and method in Section 3.5, the feature vectors can also be extracted from the experimental true IMFs. The superimposed and normalized real-time variance and instantaneous energy vectors with time are shown in Figure 15. In Figure 13, the marginal spectrum of both x-and z-directions shows the same frequency components, of about 300, 580, and 1200 Hz. Yet the z-direction IMFs have larger amplitudes relative to the x-direction, as signals tend to vibrate more in the z-direction according to previous studies.

Chatter Feature Vector Extraction from Experimental Chatter Signals
Based on the same principle and method in Section 3.5, the feature vectors can also be extracted from the experimental true IMFs. The superimposed and normalized real-time variance and instantaneous energy vectors with time are shown in Figure 15. We can see that the behavior of these two feature vectors exhibits the same change regulation as the simulation chatter signal described previously in Section 3.5. However, in order to detect and identify the onset of grinding chatter in time, it is crucial to extract feature vectors during the early stage of chatter development. It shows that the onset time of chatter is nearly 5 s, and that appropriate suppression methods should be carried out to overcome the grinding chatter. We can see that the behavior of these two feature vectors exhibits the same change regulation as the simulation chatter signal described previously in Section 3.5. However, in order to detect and identify the onset of grinding chatter in time, it is crucial to extract feature vectors during the early stage of chatter development. It shows that the onset time of chatter is nearly 5 s, and that appropriate suppression methods should be carried out to overcome the grinding chatter.

Conclusions
In this paper, bivariate empirical mode decomposition (BEMD) and Hilbert transform were demonstrated as reliable methods for processing simulated and experimental chatter signals. Moreover, the extraction criterion of intrinsic mode functions (IMFs) based on correlation coefficients successfully distinguished true IMFs from spurious components. The method of selecting the ideal number of projection directions was investigated. Additionally, real-time variance and instantaneous energy were recognized as effective predictors of chatter. Accordingly, the following conclusions can be made:

•
The BEMD and Hilbert transform methods were combined to decompose a simulated chatter signal derived from a grinding vibration signal generator. This was validated with experimental data collected from a CNC guideway grinder (KD4020X16). The results illustrate the good behavior of BEMD in terms of processing nonstationary and nonlinear signals, and providing good indicators of signal phase shifting and synchronization information. Meanwhile, Hilbert and marginal spectra accurately revealed the actual characteristics of signals.

•
The extraction criterion of true IMFs based on correlation coefficients is a reliable technique that successfully identifies and estimates the spurious components and reserves the main frequency bands of great import for extracting the actual vibration mode and corresponding features of the time and frequency domains.

•
The most ideal number of projection directions was investigated by assessing the indicators of computational time and stability degree (SD) which change along with direction number. According to the results, the ideal number of projection directions is 16. • Real-time variance and instantaneous energy were proposed as appropriate feature vectors for chatter detection and identification. These can be applied as predictors of chatter in order for chatter suppression to be applied in a timely manner.