Relative vibration identification of cutter and workpiece based on improved bidimensional empirical mode decomposition

In the process of cutting, the relative vibration between the cutter and the workpiece has an important effect on the surface topography. In this study, the bidimensional empirical mode decomposition (BEMD) method is used to identify such effect. According to Riesz transform theory, a type of isotropic monogenic signal is proposed. The boundary data is extended on the basis of a similarity principle that deals with serious boundary effect problem. The decomposition examples show that the improved BEMD can effectively solve the problem of boundary effect and decompose the original machined surface topography at multiple scales. The characteristic surface topography representing the relative vibration between the cutter and the workpiece through feature identification is selected. In addition, the spatial spectrum analysis of the extracted profile is carried out. The decimal part of the frequency ratio that has an important effect on the shape of the contour can be accurately identified through contour extraction and spatial spectrum analysis. The decomposition results of simulation and experimental surface morphology demonstrate the validity of the improved BEMD algorithm in realizing the relative vibration identification between the cutter and the workpiece.


Introduction
In ultra-precision machining, the relative dynamic displacement between the cutter and the workpiece is inevitable. Thus, changing the relative motion between the cutter and the workpiece is constantly reflected on a machined surface. The source responsible for tool-tip vibration cannot be determined exactly on the basis of existing research [1]. Vibration ripples and shape errors are formed on machined surface, which will deteriorate machining quality, aggravate the cutter wear, and reduce productivity. The weak vibration in ultra-precision machining has great influence on the surface quality of the workpiece, given that surface roughness is only a few nanometers. The surface characteristics reflecting the machining process will be reproduced on the surface. Likewise, the changes of cutting parameters and vibration in the process will be reflected. This study provides a new method for identifying the relative vibration between the cutter and the workpiece. The cutting surface can be regarded as a sensor, by which all the information of the machining process is perceived and recorded. The relative vibration between the cutter and the workpiece can be identified according to surface morphology data. This study focuses on relative vibration without considering other factors, such as the influence of impurities in the work material and the fabrication error of the diamond tool.
Researchers analyzed the main factors that affect the surface topography of ultra-precision machining, such as brittle materials [2], cutting parameters [3], tool geometries [4], and size effect [5]. Relative vibration plays an important role in the surface topography of ultra-precision machining [6,7] and is the main reason of the ripple error [8]. Cheung and Lee [9] conducted the analysis of power spectral density of surface profile and found that surface profile was mainly affected by the relative vibration of the cutter and the workpiece along the cutting trajectory. An et al. [10] established the dynamic equation of the vertical spindle supported by aerostatic bearings to study the rotation characteristics of the spindle. The influence of the impact force caused by interrupted cutting, cutting parameters, and the moment of inertia of the spindle on the deviation angle of the spindle and the frequency and amplitude of the ripple error of the cutting surface are analyzed qualitatively. Yang et al. [11] and Miao et al. [12] used Fourier and wavelet transform to analyze the contour values of machined surface along the cutting direction performed by spatial frequency domain. They also analyzed the relationship between vibration frequency and spindle natural frequency through finite element and modal test. Huang et al. [13] established the dynamic model of spindle of ultra-precision turning machine and determined the periodic vibration phenomenon of frequency aliasing of spindle under the action of dynamic unbalance. He and Zong [14] discussed the effect of the vibration between the diamond tool and workpiece along the radial direction and provided a summary of the theoretical models. He et al. [15] also innovatively established an accurate 3D surface topography model for the diamond-turning process considering kinematic factors, dynamic factors, and material defects. In terms of dynamic factors, he concluded that multifrequency vibrations determine the complication of surface topography and irregularities of surface bulge segments. Chen et al. [16] presented the multimode frequency vibration of the machine tool and its influences on the surface generation in flycutting. Zhang et al. [17] theoretically and experimentally studied the dynamic characteristics of spindle imbalance induced forced vibration and its effect on surface generation in diamond turning and concluded that it would yield a significant effect upon surface topography. Tian et al. [18] theoretically and experimentally analyzed surface generation in ultra-precision single-point diamond turning considering the basic machining parameters and the relative vibration along the cutting and feeding directions. Gao et al. [19] and Chen et al. [20] investigated the sources of frequency domain error of ultra-precision spindle and its impact on surface quality based on airinduced vibration of aero-static bearing. He and Zong [21] revealed the mechanism of the influence of multifrequency vibration on the optical performance of diamond-turned workpiece and adopted the well-established two-step process technology including optimization of the cutting parameters and the strict balance of the spindle to reduce the effect of vibration, which is experimentally proven to be sufficiently effective.
The topography of ultra-precision machined surface is composed of numerous waves with different scales superimposed on each other. To identify the relative vibration between the cutter and the workpiece, multi-scale decomposition of cutting surface is needed. At present, multi-scale decomposition of machined surface topography is widely adopted by wavelet analysis method. The original signal is decomposed into scale space to observe and analyze the surface morphology characteristics of the workpiece at different scales [22]. However, the wavelet base and decomposition layer number should be given in advance when surface morphology is decomposed by wavelet analysis method. Under the same conditions, different wavelet basis functions and decomposition layer number have exhibited great influence on the extraction of morphology features [23,24].
Scholars started to extend this treatment method to bidimensional space. For instance, Nunes [25] processed bidimensional signals as a whole and extended empirical mode decomposition (EMD) method to bidimensional space. This study is the first to implement the bidimensional empirical mode decomposition (BEMD) method and apply it to texture analysis and natural image processing. The BEMD is developed on the basis of the one-dimensional (1D) EMD, which can decompose 2D image signals at multiple scales [26,27]. Although similar to the wavelet method, BEMD does not need to consider the problem on wave basis function, which is difficult to determine in the wavelet decomposition process. Zhou and Li [28] proposed fast bidimensional ensemble EMD to denoise scheme for digital speckle pattern interferometry fringes and bidimensional intrinsic mode function (BIMF) energy estimation to reduce speckle noise.
The relative vibration of the cutter and the workpiece with a single frequency can be identified by the spatial spectrum combined with the bidimensional section profile analysis. However, when two or more different frequencies of relative vibration exist between the cutter and the workpiece, such relative vibration is difficult to determine by using spectrum analysis method. Therefore, the BEMD and spatial spectrum analysis methods are innovatively combined in this study. Experimental surfaces of ultraprecision machining, which are detected through dynamic interferometer, are decomposed to obtain BIMFs and trend term of each layer. The boundary effect problems caused by the traditional BEMD decomposition can be reduced by adopting the improved method proposed in this study. The BEMD can effectively decompose the original machined surface topography at multiple scales, and then the characteristic surface topography representing the relative vibration between the cutter and the workpiece through feature identification are selected. The relative vibration can be accurately identified by using contour extraction and spatial spectrum analysis from selected BIMFs. The results are consistent with the simulation conditions, which prove the feasibility of the improved BEMD algorithm in realizing the relative vibration identification between the cutter and the workpiece.

Relative vibration between the cutter and the workpiece
The relative vibration displacement between the cutter and the workpiece can be expressed by sinusoidal vibration or the superposition of multiple sinusoidal harmonics [14,15], which can be assumed as: where z v , A v , and f v are the relative displacement, amplitude, and frequency of mono-frequency vibration, respectively.
where φ is the spindle rotation angular, and f s denotes the spindle rotation frequency. Equation (2) is substituted into Eq. (1) to obtain the following: The forming process of machined surface can be regarded as the process of sampling dynamic cutter path. Under the action of dynamic displacement, the machined surface topography is directly determined by f v and f s . The frequency ratio f r between the two is where I is the integral part of the frequency ratio, and D is the decimal part of the frequency ratio in the range of -0.5 to 0.5.
The frequency ratio f r represents the number of times the cutter experiences the relative vibration period during the rotation of the spindle. When the cutter vibration and the spindle rotation frequency are an integral multiple of the relationship, which means D ¼ 0, the machining vibration is the same for each revolution. When D≠0, the initial phase angle of vibration at each turn is different, which will lead to obvious ripples in the feed direction.
The phase shift f after one revolution can be therefore expressed as [29] f ¼ 2πD: (5) When the relative vibration occurs between the tool and the workpiece, the wavelength of vibration pattern formed by the vibration on the workpiece surface along the feed direction is l n , which can be expressed as a function of f and the feed per revolution S fn .
Thus, the feed per revolution S fn can be further expressed as where f is the feed speed of hydrostatic guide, and ω is the spindle speed. By substituting Eqs. (5) and (7) into Eq. (6), as the following is obtained: Therefore, the spatial frequency can be deduced as 3 Bidimensional empirical mode decomposition

Traditional bidimensional empirical mode decomposition
Traditional BEMD algorithm has the advantage of multiscale and multi-resolution in bidimensional signal processing and is mostly used for image and texture analyses. The BEMD is used to decompose the signal into a series of 2D intrinsic mode functions with different scale characteristics and residual term representing the overall trend of the data. The original sample data is Fðx, yÞ, where x is the row sampling point and y is the column sampling point. The main steps of traditional BEMD algorithm are adopted from Ref. [30].

Riesz transform
To separate the surface topography errors of each frequency band, the frequency characteristics of the signal must be first determined. To obtain the frequency characteristics, constructing the corresponding bidimensional analytical signal is necessary. However, in the processing of bidimensional signals, the bidimensional analytic function constructed by the bidimensional Hilbert transform lacks multi-directional extension. On the basis of theory of Riesz transform and vector field, a type of isotropic monogenic signal is proposed. The monogenic signal is an extension of the 1D Hilbert transform in the bidimensional Euclidean space. The spatial domain of Riesz transform can be expressed as R x ðx, yÞ, R y ðx, yÞ where R x ðx, yÞ is performed with a partial differentiation to x, and R y ðx, yÞ is performed with a partial differentiation to y. For the bidimensional input signal Fðx, yÞ, its monogenic signal F M ðx, yÞ is defined as [31] where * stands for convolution. The local amplitude l a and phase l p of the bidimensional analytical signal spectrum information can be obtained through the following formula: , p 2 ½0, πÞ: (13) where i stands for the complex number and i 2 ¼ -1, and p is the phase angle. By further calculating the local phase, the local frequencies l f of physical significance are obtained: 3

.3 Boundary effect solution
The original sample data shown in the orange part of Fig. 1 is Fðx, yÞ. Figure 1 describes the schematic diagram of boundary data continuation. The solid line represents the original data, whereas the dotted line represents the extended boundary data. The first step is to extract the boundary data of the original data and mark it as Parts A and C. The second step is to extend the left and right boundary data (marked as Part A) of the original data, whereas the dotted line (marked as Part B) is the boundary data after left and right side continuation. Similarly, the upper and lower boundary data are extended (the third step), wherein extended data can be obtained. Traditional BEMD algorithm is combined with Riesz transform, and the boundary data are extended to reduce the boundary effect problems. The flowchart for the improved algorithm is shown as Fig. 2. The upper and lower envelope surfaces of the extended data are obtained by using the algorithm. The intrinsic mode functions and residual term are obtained by decomposing the original data. The results before and after the algorithm improvement will be compared and analyzed in Section 5.3. The comparison of the results indicates that the boundary effect problems can be reduced obviously, which is of great importance in evaluating surface accuracy after decomposition.

Experimental setup
The BEMD method is used to decompose the original cutting surface at multiple scales, and the surface information with different frequency characteristics is extracted from the original surface topography. After finishing cutting, the 3D surface topography of the workpiece is tested online. Fizcam 2000-100 (caliber 100 mm) dynamic interferometer from 4D company is adopted. The reflection mirror installed on the test station can be used to measure directly and reduce the deformation caused by the movement of the workpiece. Figure 3 shows the dynamic tester and the optical path of the in situ measurement. The optical path of the machined surface is directly reflected back to the interferometer by a reflector mounted on the upper part of the measured piece.
The site of experimental test is shown in Fig. 3, whereas the experiment parameters are shown in Table 1. After each processing, the machined surface is measured. The surface topography obtained after processing of the workpiece with different shapes is shown in Fig. 4.  The characteristics of surface topography caused by relative vibration are selected by applying feature identification. Then, the identification of relative vibration between the cutter and the workpiece is realized by frequency domain analysis of feature surface topography. When the cutter vibrates relative to the workpiece, the surface of the ultra-precision machined workpiece will produce a texture pattern different from the feed trace. Therefore, this feature can be used to identify the decomposed BIMF layers. Figures 5 and 6 are BIMFs and trend term of each layer obtained by BEMD method for surface topography of square and circular workpiece, respectively. As shown from Figs. 5(a) and 6(a), the first BIMF is high-frequency error, which is formed by theoretical cutting traces. Figures 5(b) and 6(b) describe that the surface topography shown by the second BIMF is the medium frequency error. Figures 5(c) and 6(c) show that the surface topography shown by the third BIMF is low-frequency error. The contour space is much larger than the feed rate, which is in line with the characteristics of the impact of the relative vibration between the cutter and the workpiece on the surface topography. Therefore, the surface topography shown from the second and the third BIMF is the characteristic surface topography of the cutting surface. Figures 5(d) and 6(d) describe that the surface topography shown from the 4th BIMF is residual term.
Through the above decomposition examples, the BEMD can effectively decompose the original machined surface topography at multiple scales. Then, the characteristic surface topography representing the relative vibration between the cutter and the workpiece through feature identification can be selected.

Feature analysis of bidimensional section profile
The surface topography is simulated and analyzed under different relative vibration conditions. Likewise, the bidimensional section profile of the cutting surface is analyzed along the feed direction by using the discrete Fourier transformation (DFT) method. The shape of the cutting surface profile is related to the cutting angle, and phase difference exists between different section profile shapes. The spatial frequency can be deduced as Eq. (9). Therefore, for the given spindle speed and feed rate, D can be obtained by using the DFT method of the feed direction contour. Table 2 shows two surface simulation schemes with different D values, and the relative vibration frequency f v is different. The I values are both 6, whereas the D values are 0.3 and 0.1, respectively.
Relative vibration identification of a single frequency is shown in Figs. 7 and 8. Figure 7 shows that single relative vibration amplitude of 10 nm and vibration frequency of 29.4 Hz are added into the simulation conditions. Case 1 in Table 2 shows Cutting parameters. Figure 8 shows that single-relative vibration amplitude of 10 nm and vibration frequency of 37.8 Hz are added into the simulation conditions. The feed direction section profiles and spatial spectrum of the simulated surfaces are analyzed. Figures 7(a) and 8(a) describe the section profile curve in feed direction. When D is non-zero, the fluctuation of   Figure 7(b) describes that the spatial frequency of the section profile in the feed direction is 14.4 mm -1 . The given feeding rate is 6 mm/min, and spindle speed is 280 r/min. According to Eq. (9), D = 14Â 6/280 = 0.3; thus D is consistent with the simulation conditions. Figure 8(b) shows that the spatial frequency of the section contour in the feed direction is 4.8 mm -1 . According to Eq. (9), D = 4.8Â 6/280&0.1, which is consistent with the simulation conditions. The relative vibration identification of the surface topography formed by the relative vibration of the cutter and the workpiece with a single frequency can be realized by the spatial spectrum combined with the bidimensional section profile analysis.
The relative vibration identification of two frequencies is shown in Figs. 9 and 10. As shown in Figs. 9(a) and 9(b), two vibration sources with vibration frequencies of 29.4 and 37.8 Hz and vibration amplitudes of 10 nm are added into the simulation surface. Figure 9(a) shows the section profile curve along the feed direction, and Fig. 9(b) shows its spatial spectrum with numerous obvious frequency components in the spatial spectrum. According to the vibration amplitude, 4.8 and 14.4 mm -1 are the main frequency components in the feed direction caused by the vibration source. Other frequency components, such as 33.6 and 47.6 mm -1 , are false frequencies or multiples of the main frequency components introduced due to the sampling length, and the existence of the radius of the diamond tool can also produce these frequency components.
As shown in Figs. 10(a) and 10(b), two vibration sources with vibration frequencies of 29.4 and 37.8 Hz and amplitudes of 10 and 2 nm, respectively, are added into the simulation surface. Numerous obvious frequency components in the spatial spectrum are observed. According to the vibration amplitude, 14.4 mm -1 is the main frequency component caused by the vibration source, corresponding to the vibration source with the vibration frequency of 29.4 Hz. According to the simulation conditions, the vibration source with a vibration frequency of 37.8 Hz and an amplitude of 2 nm is also added into the simulation    surface. However, due to its small vibration amplitude, such vibration source is even smaller than the corresponding vibration amplitude of the false frequency introduced due to the sampling length and the existence of the radius of the diamond tool or the multiple frequencies of the main frequency components such as 33.6 and 47.6 mm -1 . In this case, under the influence of vibration amplitude, the surface topography features formed by the relatively small vibration amplitude will be overwhelmed by the surface topography features formed by the relatively large vibration amplitude. Thus, the accurate vibration frequency cannot be obtained. The relative vibration identification of the surface topography formed by the relative vibration of the cutter   and the workpiece with a single frequency can be realized by the spatial spectrum combined with the bidimensional section profile analysis. However, when two or more different frequencies of relative vibration exist between the cutter and the workpiece, identifying the relative vibration between the cutter and the workpiece by directly analyzing the spatial spectrum of section profile curve is difficult. Therefore, the machined surface topography needs to be processed to separate and extract the characteristic surface topography that can reflect the relative vibration.

Identification of relative vibration
According to the analysis in Section 5.2, conventional spectral analysis method cannot identify multiple vibration sources. This section aims at verifying the effectiveness of the proposed method to identify the relative vibration between the cutter and the workpiece by simulating the cutting surface. Different D values have an important effect on the shape of the contour. The proposed method can only identify the value of D. However, the value of I cannot be identified with this method, given that the experimental surface shape of ultraprecision fly-cutting machining is asymmetrical.
The simulation parameters are adopted for the machined surface topography, as shown in Table 3. According to the analysis in Section 5.1, the second and third BIMFs of the experimental surface topography to the characteristics accord with the relative vibration between the cutter and the workpiece on the surface topography. They are considered characteristic surface topographies. DFT is used to analyze the section profile curve along the feed direction to obtain the bidimensional contour curve and the corresponding spatial spectrum diagram. Figure 11 is decomposed by the classical BEMD method, and the boundary effect problem is serious. By using the method proposed in Section 3.3, the boundary data can be extended on the basis of the similarity principle that can deal with serious boundary effect problem. The comparison results are shown in Fig. 12. Figure 13 shows the bidimensional contour curve of the second BIMF and the corresponding spatial spectrum. Figure 13(b) shows that the spatial frequency of the section profile curve in the feeding direction is 14.4 mm -1 (approximately 14 mm -1 ). The given feeding rate is 6 mm/min and the spindle speed is 280 r/min. According to Eq. (9), D = 14Â 6/280 = 0.3, which is consistent with the simulation condition Case 1. Figure 14 shows the bidimensional contour curve of the third BIMF and the corresponding spatial spectrum. Figure 14(b) shows that the spatial frequency of the section profile curve in the feed direction is 4.8 mm -1 (approximately 5 mm -1 ). The given feeding rate is 6 mm/min and the spindle speed is 280 r/min. According to Eq. (9), D = 4.8Â 6/280&0.1, which is consistent with the simulation condition Case 2.

Conclusions
The surface micromorphology of ultra-precision machined surface that contains abundant surface feature information is composed of numerous waves with different scales superimposed on each other. The vibration in the process will be reflected in the surface characteristics. This study provides a new method for the identification of the relative vibration between the cutter and the workpiece, which can be identified according to surface morphology data. The following conclusions are drawn in this study.
1) BEMD is applied to the identification of relative vibration stored in the surface topography in ultraprecision machining. Although the vibration signal is extremely weak, it can be accurately identified by the method introduced in this study.
2) According to Riesz transform and vector field, a type of isotropic monogenic signal is proposed to be combined with BEMD. The monogenic signal is an extension of the 1D Hilbert transform in the bidimensional Euclidean space.
3) The similarity principle is used and the serious boundary effect problem in the decomposition process of surface morphology can be overcome by extending boundary data. 4) Decimal part of the frequency ratio which has an  important effect on the shape of the contour can be accurately identified through contour extraction and spatial spectrum analysis. The decomposition results of simulation and experimental surface morphology have demonstrated the validity of the improved BEMD algorithm.