A statistical filtering method for denoising of micro-force measurements

Precise measurement of mechanical forces is crucial to efficient micro-manufacturing. The quality of such measurements depends heavily on the properties of the noise inevitably accompanying every measurement process. In the micro-range, the signal-to-noise ratio tends to be very low, and the noise dynamic varies for different frequencies. In result, common denoising methods that assume white noise perform poorly in this setting. In this paper, a novel, easily implementable denoising method based on a local statistic of the measured data’s spectrum is proposed. By testing it on a representative dataset, it is shown that the proposed method is robust and stable. Particularly, it allows for an efficient retrieval of the force signal encountered in micro-milling processes.


Introduction
In recent years, efficient and flexible production of small components has become indispensable to high-precision mechatronics. This increased demand in miniaturization had an immense impact on the manufacturing technology and led to the development of many novel and innovative processing techniques for high-precision micro-parts (cf., e.g., [39]).
Micro-milling is a prominent method for rapid and cheap production of micro-components. Together with other material removal processes, it has become one of the most promising methods for the production of micro-forming tools as well as micro-components (cf., e.g., [12,39]). This is especially true for the manufacturing of components with high surface quality and with textured or structured surfaces (cf. [3,37]). For end products where the desired functionality is determined by these surface properties micro-milling is a crucial micro-production method. Accordingly, in order to satisfy the desired requirements and properties of the final product optimal process parameters have to be achieved and maintained within the manufacturing process (cf. [8,30,35,37]).
To determine optimal process parameters, efficient models and simulations are needed. The predictive power of such models depends heavily on the modeling of the cutting forces with respect to the process parameters and the desired structures. Furthermore, in order to maintain optimal parameters and thereby in order to avoid unsatisfactory manufacturing results careful monitoring of the process is indispensable. In this case, also, the analysis of micro-forces plays an important role in the detection of aberrations from the desired process parameters.
Altogether, with the focus on surface texture and tool wear, it is therefore essential to model and predict microforces involved in micro-milling precisely. Consequently, a thorough knowledge of these forces is necessary to guarantee the high quality of the manufactured micro-components.
Usually, the force models in cutting processes are based on the common assumption that the force is proportional to the cross-sectional area [2,40]. However, the experimental data indicates that this approach is not well-equipped to fully describe effects occurring in micro-processes. Further, conventional methods and models cannot be scaled effectively to the micro-scale, due to the occurring size effects [38]. Consequently, the development of new forcemodels for micro-processes is of uttermost importance for manufacturing technology.
The detailed research on turning processes and model development for complex surface generation is presented in [25]. In particular, there, the authors consider a dynamic cutting force model, a regenerative vibration model, a machining system response model, and a tool profile model.
Further, a two-degrees-of-freedom model which includes dynamical changes caused by forces acting during the micro-process was presented and discussed for turning processes in [11,27]. This model was successfully further developed and extended to the micro-end-ball-milling processes in [28]. In particular, the extended, more general force model is based on the tool-workpiece-interaction. An analytical force model for micro-end-milling including tool run-out and tool wear effects was investigated in a series of papers by Bao et al. see [4][5][6]. In addition, Li et al. [24] proposed a model of a three dimensional cutting force for micro-end-milling operation which includes a new nominal uncut chip thickness algorithm. This model was developed under the consideration of tool run-out and an exact trochoidal trajectory of the tool tip. Similar investigation together with a finite element model for the orthogonal cutting is presented in [1]. Finally, additional investigations of micro-milling forces and their models were carried out in [9,20] and [21].
However, many questions concerning force models for micro-cutting processes are still open. In order to develop new, improved models as well as to qualitatively verify the capability of already existing models to precisely predict the micro-forces encountered in such micro-processes a reliable recovery of these forces from noisy measurements is indispensable.
As mentioned above, unanticipated manufacturing errors of various kinds may occur during micro-processes. In order to correctly and reliably assess the state of such processes monitoring is necessary. However, monitoring methods for miniaturized manufacturing processes are more complex and challenging than their counterparts in the conventional manufacturing. A vast number of methods has been proposed for the monitoring of micro-processes. The most common among them are the measurement of microforces; further, the measurement of audio signals generated by the process; and, finally, the recording and analysis of the process by video (cf. [19,26,44]). Of course, all of these methods may be considered separately as well as in combination of each other.
In particular, in the case of the tool monitoring, one of the most significant problems is the tool wear. In addition to the above mentioned methods, a monitoring and tool wear prediction method based on accelerometer measurements was presented by Stavropoulos et al. in [33]. Another, innovative approach to process monitoring for micro-milling operations has been studied by Kuram and Ozcelik in [23].
Finally, in order to determine the correct machining conditions, chatter and/or vibrations have also to be taken into consideration. To that end and to ensure the maximal efficiency in monitoring, a reliable measurement of micro-forces is required.
Altogether, it can be summarized that accurate measurement micro-forces in micro-milling processes are quintessential for achieving and maintaining optimal results within the manufacturing process.
However, a reliable quantification of micro-forces within measured signals is a highly challenging task. As every other real-world measurement process, the measurement of process forces in micro-range is inevitably accompanied by noise. Therefore, the quality-and subsequently the usefulness--of the acquired data depends on the ability to separate and remove the noise part from the measured data.
Forces encountered in micro-milling may have low magnitude as compared to the noise of the measurement. Therefore, a low signal-to-noise ratio may be expected in these measurements. In particular, this tends to be the case for high-speed micro-milling processes. Consequently denoising methods tailored to handle micro-forces in microprocesses are of fundamental interest to practitioners.

Approach
The main aim of this paper is to provide a fast and reliable method for removing the noise parts from the force signals measured in micro-milling processes. The presented method has very low computational complexity, allowing for easy implementation.
Currently, already a plethora of standard denoising methods for time-dependent signals is currently available to practitioners. These range from classical Wiener-filtering (i.e., Fourier filtering) methods (cf. [41]), more recent Gabor Transform and Short Time Fourier Transform (STFT)-based methods (cf. [7,17] and the references therein) up to wavelet-based methods (cf. [14,16] and the references therein). Further practices include non-Fourier based, nonlinear transforms like the empirical mode decomposition (cf. [10,18]). Finally, we also mention the recently developed, patch-based methods [13].
These general methods have been successfully applied to the problem of denoising force signals in micro-milling. For example, results for denoising based on independent component analysis (ICA), Short Time Fourier Transform, and wavelets are presented by Zhu et al. in [44]. Another approach based on sparse representation of signals in Short Time Fourier Transform has been studied by Zhu et al. in [45]. Further study and denoising of force signals by means of independent component analysis has been discussed in [34,42,43].
However, a typical, underlying assumption of these standard, general denoising methods is that the noise is white and Gaussian. This means that the noise level is constant for all spectral (Fourier) frequencies. This assumption may not be necessary true for measurement of mechanical forces in the micro-range. Exemplary spectra of the three force components where the white noise assumption is invalid are depicted in Fig. 1.
Denoising methods operating in the original domain of the signal are essentially smoothing the signal. Therefore, they are virtually low-pass filters. As such they remove high-frequency features important for the analysis of force models. On the other hand, denoising methods operating in the frequency domain are usually thresholding the spectrum. Referring again to Fig. 1, it is observed that in order to remove the noise completely the threshold has to be chosen high. Therefore, a filter based on hard-thresholding will not preserve the information on the higher harmonics. Consequently, it can be stated that standard denoising methods generally fail to produce high-quality results in the setting depicted in Fig. 1.
Hence, in order to perform the separation reliably and efficiently, the characteristic properties of the noise and of the force signal have to be determined. To that end, typical power spectra of the three components of the force depicted in Fig. 1 are considered. The following, crucial points can be derived from visual inspection of this figure: As already mentioned above, the noise level is dynamic, i.e., the noise-level present in the data is frequency dependent, and thus, is not white. Consequently, the local noise level varies for different frequencies. Further, it is observed that the force parts of the signal consist of peaks or outliers in the power spectrum which are distinctively higher than the surrounding, local noise level. Furthermore, these peaks follow roughly the dynamic properties of the noise, i.e., they are high where the noise is high and are low where the noise is low.
Additionally, the signal exhibits a clear fundamental frequency (which is connected to the rotational frequency f as well as the rotational frequency per cutting edge f k , where k is the number of the cutting edges). This fundamental frequency (and the related peaks in the plots of Fig. 1) is accompanied by higher order harmonics (and their related peaks in the plots of Fig. 1) at integer multiples of the fundamental frequency. In the data depicted in Fig. 1, the fundamental frequency is n = 40, 000min −1 ≈ 666.667s −1 and therefore higher harmonics are expected at multiples of f = 666.667 Hz (rotational frequency).
Previous literature (cf. [2,36]) and research of the authors (cf. [28]) indicates that preservation of the information about the higher harmonics within the filtering process is crucial to the modeling of the underlying forces. Therefore, in the setting described in this paper, an efficient filter must preserve information of higher harmonics of the fundamental frequency.
As observed above, the essential part of the force signal consists of peaks in the spectrum. Therefore, in order to provide efficient denoising, the related filter has to be able to detect peaks for variable levels of noise. Further, the filter has to be fast in order to provide the denoised signal in real time, e.g., for monitoring applications. Finally, in order to be easily implementable, the filter should have low mathematical and computational complexity. We close this section with a rough overview of the proposed method. A detailed description is provided in Section 5. Further, a visualization of the main steps is presented in Fig. 2. It is also remarked that the proposed method operates on individual force components separately.
In a first step, the (Fourier) power spectrum is computed by means of the Fast Fourier Transform (FFT) from the measured force signal.
Then, this data is employed to derive a peak indicator. The peak indicator is based on a gliding squared coefficient of variation (SCV) of the power spectrum. It can be shown that the value of the SCV in the absence of peaks has a value which does not depend on the noise level (i.e., SCV is agnostic with respect to the noise level). Further, where peaks are present the SCV have much higher value. The reasons for this behavior are presented in detail in Section 4.
In the next step, the Fourier spectrum is decimated according to the value of SCV. This means that the Fourier coefficients for all frequencies for which the SCV is smaller than a certain value (chosen by the user) are set to zero, i.e., the Fourier spectrum is thresholded with respect to the SCV.
Finally, the filtered force signal is derived from the decimated spectrum by Inverse Fast Fourier Transform (IFFT).
The gliding SCV is fast and easy to implement, since it can be reduced to two gliding mean filters (cf. Section 4). Consequently, the main computational steps of the approach, namely FFT, SCV, decimation, and IFFT are fast and either already available in standard software packages (FFT and IFFT) or very easy to implement (SCV, decimation). Altogether, the proposed approach fulfills all of the above requirements and provides a simple yet powerful denoising method for the measurement of process forces in micro-regime. Results concerning the performance of the filter are presented in Section 6.

Experimental setup
All cutting processes were carried out on a DMG Sauer Ultrasonic 20 linear machine tool. The machine tool design and a high-speed spindle (maximum rotational speed 42, 000 min −1 ) allow to conduct micro-milling operations. The workpiece material for the cutting processes was hardened cold working steel 1.2379 with a hardness of 60 HRC. The steel was produced through powder metallurgy and exhibits a fine grained micro-structure without primary carbides, making this appropriate for micro-manufacturing.
A Kistler 9119AA2 MiniDyn piezo-electric multicomponent force transducer was used to measure process forces. It allows for the simultaneous acquisition of three orthogonal forces (F x , F y , and F z ) or three orthogonal momentums (MA, MB, and MC). The sensitivity of the force transducer is given by the manufacturer as less than 2 mN and the force measurement range with −2.5 kN to 2.5 kN for all measurement directions. The force transducer was mounted to an adapter plate, attached to the  The measurement chain comprises a Kistler 5019 multichannel charge amplifier as well as a computer equipped with an I/O A/D converter card for data sampling. The sampling rate was 30 kHz for all cutting processes. All considered processes were slot milling processes.

Identification of the cutting processes in the data
A measurement campaign consists of a set of three repeated cutting processes which are carried out for an identical set of parameters. An individual cutting process consists of the measurement of the force for two known depth of cut levels a 1 p and a 2 p with a 1 p > a 2 p . The milling tool is moving through the material with the depth of cut a 1 p and after time t 1 changes to level a 2 p . In sum, all three force components F x , F y , F y are being measured in six individual cutting processes within one campaign. For example, the force measurements in x-direction for one whole campaign are presented in Fig. 4 (light gray line).
The cutting processes are separated in time. However, the boundaries of these processes are not synchronized with the force measurement aperture. Therefore, in the first step, these boundaries have to be recovered from the data measured within the whole campaign.
Because of the high level of noise in the measured data, it is not possible to determine the exact boundaries of the milling processes and a smoothing step is carried out. After smoothing the data, individual cutting processes at the different depths of cut can be identified. In Fig. 4 (black line), the result of smoothing with Daubechies wavelet filter (of order 8) is presented. The boundaries of the individual processes-the sections divided by dark gray, vertical lines in Fig. 4-are then clearly distinguishable. In the subsequent steps, every process is considered separately. Only one force component for one depth of cut is presented exemplarily in Fig. 4, as the other two were processed in the same way. However, it is remarked that the smoothed data is only used in order to identify the boundaries of the individual milling processes. The processed data (black line in Fig. 4) is too over-smoothed to be useful in further processing. For this reason, in the subsequent steps, the original datarestricted to the temporal boundaries of each process-is employed again.

Statistical peak detection
The proposed filter for the individual milling processes is developed in two parts. In this section, a local peakdetection criterion is discussed. This criterion will then be developed into a full denoising method in the next section.
A pure white noise signal, e.g., a vector of length K defined via  u ) are sample vectors of two, normally distributed, stochastically independent random variables with zero mean and variance σ 2 · K 2 . Thus, they can also be considered pure-noise signals. Consequently, the vector 2 ) 2 may be interpreted as a sample vector of a Rayleigh-distributed random variable with the scale parameter σ · √ K/2. The expected value of that random variable is given via √ π/2·σ and the variance by 4−π 2 σ 2 , (cf. [15] p. 173). Hence, the so called squared coefficient of variation (SCV) of that random variable, which is the ratio of the variance to the square of the expected value, is independent of the scale parameter of the Rayleigh distribution and has the value 4 π − 1 ≈ 0.27. An estimator for SCV is given by: and → w : are the standard unbiased estimators for the expected value (mean) and for the variance. Of course, this estimator is by itself not expected to be unbiased. However, due to numerical experiments and the analogue computations for a normal distributed random variable (cf. [32] p. 58.) the error is expected to be negligible.
In other words, for pure-noise signals It should be noted that the assumptions on the distribution properties of → v 1 and → v 2 were only used to derive the peak detector. In order for the filter to be sensible, the only important property is that → w is a sample of a Rayleigh-distributed random variable. As can be seen in Fig. 5, this assumption is fulfilled for high-frequency part of the data (higher than 5 kHz). This frequency was chosen since there the noise level is approximately stationary. Obviously, the spectrum is indeed Rayleigh-distributed there. Therefore, it can be stated that the data (locally) fulfills the assumptions of the filter.
It should be stressed that the detector described above is scale invariant, i.e., the scale parameter of the Rayleighdistribution, respectively, the noise level is not parameters within the detector and therefore have neither to be estimated nor to be known. Further, the filter may be applied to all signals which are (locally) stationary in the power spectrum. Finally, it should also be emphasized that the pointwise absolute-value of the Fourier transform is used as the definition of spectral power of a signal. Another common definition for the spectral power is the square of the absolute-value. Of course in that case the SCV can be employed to construct a scale invariant peak detector, too. However, in our experiments, the definition used here achieved better performance.

Description of the filter
As was mentioned before, it is assumed that all three force components F x , F y , F y are measured and thus also have to be filtered. The data for every of the six individual experimental parts of the data (cf. Section 3.1) is denoised separately. The three force components for the measured, noisy data of a single process are denoted by → F x , → F y , → F z -an index identifying the process is dropped in order to simplify the notation. As the filter treats every component separately the action of the filter is presented only for the force com- x ] of the data → P x . Further, the data is extended symmetrically at the boundaries. For example, for N = 2 and m = 1 the extended subset is given by [p (3) x , p (2) x , p (1) x , p (2) x , p (3) x ]. -The parameter N controls the width of typically detected peaks, while the parameter c described in the last section controls the height of typically detected peaks. Together, these two parameters control the quality of the denoising. -The proposed filter is a gliding filter of width K = 2N + 1. In particular, every peak is usually surrounded by a window of width 2N + 1 where the noise is not filtered. In the given dataset, this window was helpful for the modeling of the micro-forces. However, if such behavior is not desired, additional morphological operations [31] may be applied to narrow or even delete that surrounding window.
-Since it is known that → F x is real-valued, the left half of the FFT is complex conjugate of the right half (if the vector is considered to be a row vector). As such, both halves contain the same information. Therefore, subsequently, the vector → G x could also be cropped to the left half, i.e., it is sufficient to only consider entries with index 1 ≤ m ≤ L, where L := M/2 and · denotes the rounding towards infinity (also known as ceil-function).

Results
In this section, the performance of the filter in the spectral domain and in the time domain based on real data is described. Since in the tested dataset, the timestamps were equidistant, FFT was used. The algorithms were implemented in MATLAB and have been tested for several parameter settings.

Performance in the spectral domain
In Fig. 6 (top row), the spectrum of a milling process' signal is depicted which mainly exhibits peaks at the first and second multiple of the rotational frequency n = 40, 000/min (≈ 666.67 Hz). These peaks are recognizable by visual inspection. The width N of the filter was chosen to be 6 in all three force components while the peak detection parameter c was chosen 10 for the F x force component, 7 for the F y force component, and 15 for the F z force component of the force signal. In the bottom row of the figure, it can be seen that the filter successfully recognized and retained these peaks. At this point, it is once again stressed that the filter is agnostic with respect to the information on the rotational frequency, i.e., that information was not fed to the filter. In Fig. 7, the performance of the filter is shown for signals which exhibit peaks at more than only two multiples of the rotational frequency. In this case, the width N was chosen to be 12 and the peak detection parameter c was chosen to be 5 in all three force components. As before, the top row of the figure shows that the filter is at least as good as choosing the peaks by visual inspection. Furthermore, the bottom row shows that also in this case the filter recovers the correct peaks at the multiples of the rotational frequency. It is worth pointing out that characteristic peaks at about 250 Hz are recognizable in all three force components as well. Although the mechanical interpretation of this frequency is currently still an open question. It is clear that these peaks have different characteristics than the surrounding noise and, therefore, are correctly recognized as 'not-noise' part of the signal.
Notice that in both cases (Figs. 6 and 7), the filter is at least as good as traditional filters. This is especially evident for the force component in x-direction (left row) and z-direction (right row).
The considered processes were slot milling processes with feed velocity in x-direction. As such, the F y component of the cutting force is expected to have much less pronounced peaks then the other two components. This is clearly visible in the middle column of Figs. 6 and 7. In order to remove the noise caused by the machine table (best visible at frequencies 2.5-3 kHz) a filter based on spectral thresholding would also remove all peaks belonging to the signal, as these peaks have significantly smaller magnitude than the surrounding noise. However, the presented filter is clearly able to recover the correct peaks for the signal.
The filter was designed to be scale invariant and therefore to automatically adapt to the level of noise present in  Fig. 6 the spectrum. This fact is depicted in the subplots in the upper right corner of plots of Fig. 7. Firstly, these sub-graphs show that information about the force component is present even at higher multiples of the rotational frequency (in the depicted case at the sixth and seventh multiple). Secondly, it shows that the magnitude of this information is smaller than the magnitude at low multiples (e.g., the first and second multiple). Finally, this shows that the filter works in the way it was designed and successfully retains the force information while removing the noise from the signal, even if the noise level varies within the signal spectrum.
Summarizing, this shows that the presented filter is often as good and in critical cases superior to the established spectral thresholding methods.

Performance in the time domain and comparison with force model
In order to fully evaluate the performance of the filter, the results in the time domain are considered also. In Figs. 8 and 9, the measured and filtered F x (top) and F z (bottom) force components for the same process as in Fig. 6 is depicted. The force component in y-direction is omitted, as the process is a slot milling process in the x-direction. Consequently that force component is of less interest as it is expected to be much less variable than the other two components (a view which is confirmed by the much lower magnitudes of the relevant spectral information, as shown in Fig. 6). It is remarked that the unfiltered signal (gray line in Fig. 8) is relatively smooth. Therefore, smoothing-based filters are expected to perform poorly for that signal. However, this also means that a smoothness-based evaluation of the performance of the filter not is not the right choice. Therefore, instead, the filtered force measurement is compared with an well-established force model.
For conventional cutting processes, a number of cutting force models have been described in the literature. In the case of milling processes, the acting forces are separated into three components, radial force, tangential force, and axial force and are denoted by F r , F t , and F a , respectively. A common approach relies on fact that the force components are proportional to the cross section area of cut A c and can be computed using the specific cutting force components K rc , K tc , and K ac . This results in the forces F r = K rc A c , F t = K tc A c , F a = K ac A c . By applying the coordinate transformation operator, the force components F x , F y , and F z in the Cartesian tool coordinate system are obtained. This cutting force model was first introduced by Weck and Teipel [40]. Altintas extend this model by adding a second term which additionally takes into account friction, cf. [2, p. [35][36][37][38][39][40][41][42][43][44][45][46][47]. The model of Altintas is applied in the remainder of this section. Accordingly, the model for the (instantaneous) cross section area A c (t), the (instantaneous) chip thickness h(t) and the instantaneous angle of immersion φ(t) are given by where the depth of cut a p , the feed velocity v f , rotational frequency n, and the number of cutting edges k are not time varying within a single milling process. However, they may have different values for different milling processes. The quantity t is the time measured in seconds, and φ 0 is the initial angle of immersion. All other quantities have the dimensions given in Figs. 6 and 7. The model of Altintas then gives the following formulas for radial-, tangential-, and axial-force components: F t (t) = K tc A c (t) + K te a p (7) F r (t) = K rc A c (t) + K re a p (8) F a (t) = K ac A c (t) + K ae a p A visualization of the forces in the x, y-plane is presented in Fig. 10. Finally, the forces in the Cartesian tool coordinate system are given via F y (t) = +F t (t) · sin φ(t) − F r (t) · cos φ(t) + c y (11) F z (t) = +F a (t) + c z (12) where the offset parameters c x , c y , and c z are additionally introduced to counter the offset caused by the measurement aperture.
It should be noted that there are ten unknown parameters φ 0 , K tc , K rc , K ac , K te , K re , K ae , c x , c y , c z . These are recovered from the filtered data by least square optimization. The results of the optimization are depicted in Fig. 9 for the F x and F z force component. The component F y is once again omitted for the reasons mentioned above, i.e., the nature of the milling process at hand.
A very good fit between the filtered force component (gray line, top of Fig. 9) and the optimized model (black line, top of Fig. 9) is achieved for the F x force component. This validates the effectiveness of the filter in that case.
The difference between the filtered signal (gray line, bottom of Fig. 9) and the optimized model (black line, bottom of Fig. 9) for the F z direction can be accounted for by the fact that the macro-scale milling models (like the Altintas model) do only describe the chip thickness h in terms of the feed and feed per tooth (which translates to the parameters v f and k in the above definition of h(t)). Therefore, the rotational frequency is dominant in such models in the axial force component F a and therefore also in the F z force component.
Summarizing, the presented filtering is as good as the model-based constraints allow it to be. However, the development of force models for milling processes in the microscale is currently an ongoing research effort and the comparison of the filter to these models is an interesting future research endeavor.

Conclusions
A new approach for signal denoising based on a statistical peak detection methods is presented. It can be applied to a wide range of noise models. Moreover, no assumptions from the milling process are needed to model the signal. The only necessary assumption is that the signal parts exhibit a peaklike behavior in the Fourier spectrum. The filter is controlled by two parameters which have to be chosen in advance by the user. The filtering method has succeeded to determine the fundamental frequency as well as its harmonics from the input signal.
The successful application of the presented filter combined with further cutting experiments and theoretical research will allow the development of novel cutting force modeling techniques which account for the higher multiples of the rotational frequency within the signal.