Rotating Machinery Diagnosing in Non-Stationary Conditions with Empirical Mode Decomposition-Based Wavelet Leaders Multifractal Spectra

Diagnosing the condition of rotating machines by non-invasive methods is based on the analysis of dynamic signals from sensors mounted on the machine—such as vibration, velocity, or acceleration sensors; torque meters; force sensors; pressure sensors; etc. The article presents a new method combining the empirical mode decomposition algorithm with wavelet leader multifractal formalism applied to diagnosing damages of rotating machines in non-stationary conditions. The development of damage causes an increase in the level of multifractality of the signal. The multifractal spectrum obtained as a result of the algorithm changes its shape. Diagnosis is based on the classification of the features of this spectrum. The method is effective in relation to faults causing impulse responses in the dynamic signal registered by the sensors. The method has been illustrated with examples of vibration signals of rotating machines recorded on a laboratory stand, as well as on real objects.


Introduction
When properly processed, dynamic signals recorded during the operation of rotating machines are a valuable source of information about the condition of machines. Algorithms are still being developed, which will quickly, unequivocally, and automatically determine the operational state of the machine. One of the most commonly used methods in diagnostics are vibration sensors due to the ease of assembly, the possibility of obtaining information without disassembling the machine, and the wealth of information obtained about the condition of the machine in various frequency ranges. In the case of machine operation with constant speed and load, the most effective method of analysis is the frequency spectrum [1]. However, due to the fluctuating rotational speed of the machine, often due to a variable load, the frequency spectrum becomes blurred. The real challenge is to diagnose machines in non-stationary conditions. In order to avoid resampling in the analysis of dynamic signals, time-frequency methods are used [2]. The authors present a systematic review of over 20 time-frequency methods used to detect machine damage. One of the simplest is the short-time Fourier transform STFT [3]. The improved version of the spectrogram is the Wigner-Ville analysis, the use of which in diagnosing damage to gears is described in [4,5]. The use of wavelet analysis in diagnosing damage to rotating machines can be found, among others in [6,7]. Other papers [8][9][10] present a method of detecting and diagnosing gear damage on the basis of a reference model based on the signal averaging technique. The proposed algorithm first establishes an autoregressive (AR) model of the gear vibration signal in its base state, and then diagnoses the state based on the residual signal. The use of the synchronous method (order-tracking analysis) allows to synchronize the diagnostic signal with the rotation of the diagnosed machine, as well as to eliminate the influence of background noise from other sources [11,12]. The article [13] proposes a method of informative selection of frequency bands. It utilizes the approach of non-negative matrix factorization applied to time-frequency signal representation. For the analysis of non-stationary signals, the empirical mode decomposition EMD algorithm is also used, which decomposes the signal into modes in an empirical manner [14,15]. Since the classic method can mix two or three modes, various modifications of the EMD method are developed [16][17][18].
In recent years, more and more often, efforts are being made to automate the diagnosis process, and qualified diagnosticians are to be replaced by diagnostic algorithms. So-called data-driven statistical methods can also be successfully used in diagnosing machines in non-stationary conditions. Determining the characteristics of a dynamic signal and its classification is carried out by machine learning algorithms [19][20][21][22]. More and more papers are devoted to the diagnostics of machines with the methods of deep learning [23][24][25][26][27][28][29]. Methods that allow detecting failures at an early stage, especially in non-stationary conditions that often occur in real objects, are still being searched and developed.
The fractal theory has unique advantages in dealing with transient, nonlinear, and non-stationary signals. It is a new trend in solving practical problems involving the study of short-term signals in a transient process with strong nonlinearity, based on the fractal theory in order to diagnose damage to mechanical equipment. The impact of damage to the outer race of a rolling bearing on the width of the multifractal spectrum is presented in [30]. Detrended fluctuation analysis MF-DFA is a commonly used algorithm for multifractal analysis. The method is often used to diagnose damage to rolling bearings [31][32][33] and gears [34][35][36]. In [37], the MF-DFA algorithm was used for the analysis of frictional vibrations, where the ensemble empirical mode decomposition EEMD was used to denoise the signal. The adaptive MF-DFA algorithm is presented in [38]. The disadvantage of the MF-DFA algorithm is its sensitivity to analysis parameters. These disadvantages are not present in the method of wavelet leaders (WLMF). The wavelet leaders algorithm was used to analyze the images [39] and diagnose damage to rotating machines [40].
The article proposes a novel EMD-WLMF method consisting in the multifractal analysis of the first empirical component (IMF1) of the vibration signal using WLMF algorithm. As a result, a multifractal spectrum is obtained, the parameters of which are used to distinguish the state of the machine. Since the tested damage causes impulse disturbances of the vibration signal, the signal is first filtered out using the empirical decomposition of the signal EMD. Only the first intrinsic mode function (IMF1) containing high-frequency information is subject to multifractal analysis. The method has been illustrated with examples of signals recorded on a laboratory stand as well as on real objects. The presented examples, especially the analysis of signals recorded in the vehicle, show that the combination of EMD and WLMF methods allows the detection of damage at an early stage in conditions unfavorable to diagnosis-i.e., at floating rotational speed and load.
The article is organized as follows. Section 2 describes the methods of signal analysis used in the study, i.e., the WLMF algorithm, as well as determining the first intrinsic mode function (IMF1) using the EMD method. In Section 3, the described algorithm for fault diagnosis was used on the basis of signals recorded on a laboratory stand and on real objects under operating conditions. The discussion is presented in Section 4.

Multifractal Formalism
Multifractal analysis, which is based on estimated signal scaling exponents, is a popular statistical tool for studying time series. The mathematical formalism is based on the increments of their values, the measure of which are the pointwise Holder exponents h of the time function x(t) at the point t 0 , determined by the supremum of all exponents satisfying, for the constant C > 0, the condition where P n (t − t 0 ) is a polynomial of the order n < h. The pointwise Holder exponent describes the regularity of the function. The greater the local regularity of the time series, the higher the singularity exponent. The multifractal spectrum D(h) represents the histogram of the Holder exponents and determines the fractal dimensions of the subsets of singularities with a given exponent value [41][42][43].
In the time-frequency approach to signal analysis, the WTMM (wavelet transform modulus maxima) method was initially used, which is based on a continuous wavelet transform (CWT). It consists in determining the chains of local, the structure function and its scaling exponents, as well as the Legendre transform. Practical implementations of such an algorithm showed a number of faults that make it impossible to conduct research for some types of real signals [44].
Another approach to the problem of estimating local scaling exponents as a method of studying the regularity of time series and their multifractality is related to the multifractal detrended fluctuation analysis (MF-DFA) method. MF-DFA enables to study the observed signals in terms of their multifractality, provides a more stable approach to multifractal formalism than the WTMM method [31,36,[45][46][47].
The formalism in the time-frequency domain used in the paper allows for the estimation of multifractal parameters using wavelet leaders (WLMF), which are representatives of local Holder exponents of the signal. In the case of wavelet coefficients centered around zero values, it is difficult to guarantee numerical stability. This problem does not arise in the case of the WLMF method, based on the wavelet coefficients obtained as a result of the discrete wavelet transform (DWT). The next steps of the algorithm include the selection of coefficients called wavelet leaders, determination of the structure function with scaling exponents, and the multifractal spectrum. The algorithm shows low computational costs, numerical stability, and high versatility in terms of real signals [39,40].
The scaling exponents of the structural function do not depend on the selection of the wavelet, provided that the number of zero wavelet moments is two times greater than the largest exponent of the signal holder. The mother wavelet must be orthogonal or biorthogonal. In the work, the Daubechies (db2) wavelet mother was arbitrarily selected.
For the coefficients (2) of the discrete wavelet transform (DWT) of the function x(t) and the basic wavelet with a compact support ψ 0 (t) wavelet leaders (3), for the set of the largest coefficients d x (j , k ) ≡ d λ in the vicinity of 3λ, are defined by the dependence on any scale where j, k are integers and 3λ := 3λ j,k = λ j,k−1 ∪ λ j,k ∪ λ j,k+1 and λ := λ j,k = k2 j , (k It can be shown that the Holder exponents are the scaling exponents of wavelet leaders: L x (j, k) ∼ 2 jh . Moreover, the structure function (4) defined for wavelet leaders is described by the power dependence, the exponent of which is the multifractal scaling exponent ζ(q) : R → R .
where q is the order of the structure function and n j is the number of intervals of the multiresolution analysis. The function obtained by the Legendre transformation of the multifractal scaling exponent ζ(q), under mild conditions of signal regularity, is the upper limit for the multifractal spectrum (5) of the tested signal The description of dynamic properties of systems is successfully carried out on the basis of the parameters of multifractal spectra of representative time series.
The following were selected as parameters of the multifractal spectrum D(h) related to singularities h representing the local scaling of the measure in different places of the time series: • multifractality level, representing the heterogeneity of the signal under study, The procedure of the method is shown in Figure 1.
It can be shown that the Holder exponents are the scaling exponents of wavelet leaders: ( , )~2 ℎ . Moreover, the structure function (4) defined for wavelet leaders is described by the power dependence, the exponent of which is the multifractal scaling exponent ( ): → .
where is the order of the structure function and is the number of intervals of the multiresolution analysis.
The function obtained by the Legendre transformation of the multifractal scaling exponent ( ), under mild conditions of signal regularity, is the upper limit for the multifractal spectrum (5) of the tested signal The description of dynamic properties of systems is successfully carried out on the basis of the parameters of multifractal spectra of representative time series.
The following were selected as parameters of the multifractal spectrum (ℎ) related to singularities ℎ representing the local scaling of the measure in different places of the time series: The procedure of the method is shown in Figure 1.  The effect of the algorithm is illustrated in Figure 2 for a harmonic signal and various types of noise, examined in terms of multifractality. Figure 2a shows sample waveforms of signals, Figure 2b shows their multifractal spectra with the characteristic points. The harmonic signal, three types of random signals varying in the level of multifractality, and the sinusoidal signal disturbed by the multifractal signal were analyzed. The level of multifractality can be measured by the width of the multifractal spectrum. Thus, both the The effect of the algorithm is illustrated in Figure 2 for a harmonic signal and various types of noise, examined in terms of multifractality. Figure 2a shows sample waveforms of signals, Figure 2b shows their multifractal spectra with the characteristic points. The harmonic signal, three types of random signals varying in the level of multifractality, and the sinusoidal signal disturbed by the multifractal signal were analyzed. The level of multifractality can be measured by the width of the multifractal spectrum. Thus, both the harmonic signal and the white noise cannot be treated as multifractal signals. Multifractal signals are characterized by impulse disturbances, such as harmonic signals with noise.  Figure 2b shows the characteristic points of the multifractal spectrum, which were used to define the characteristics of the spectrum:  Figure 2b shows the characteristic points of the multifractal spectrum, which were used to define the characteristics of the spectrum: • Spectrum asymmetry Sensors 2021, 21, 7677 6 of 21

•
The singularity exponent with the greatest fractal dimension Figure 2c shows the multifractal spectra determined when selecting the Daubechies mother wavelet of orders 1 to 6. As the selection of the wavelet order has little effect on the parameters of the spectrum, the Daubechies wavelet of order 2 ('db2') was arbitrarily selected. The selection of Daubechies wavelets of different orders should yield consistent results when used repeatedly.

Empirical Mode Decomposition
The EMD is an iterative numerical approximation algorithm designed to extract intristic mode functions IMFs from signals adaptively by cubic spline interpolation according to the local characteristic time-scale. It decomposes a signal into IMFs via iterative sifting.
Any signal can be approximated by a superposition of a series of IMFs as written by where c i (t) is the i-th IMF, and r n (t) is the residual signal which represents the slowly varying or constant trend of the signal.
Applying the Hilbert transform [14] to each IMF in Equation (9), we can construct the corresponding analytic signal. Then we can compute the instantaneous frequency via the derivative of instantaneous phase relative to time. By expressing the analytic signals in polar coordinate form, and taking the real part, we can obtain the Hilbert amplitude spectrum and Hilbert energy spectrum where δ(·) is the Dirac delta function. The EMD is the core of the Hilbert-Huang transform. For a real signal x(t), the EMD procedure is as follows.

2.
Extract the i-th IMF: Find the local minima and the local maxima of h ij (t). c.
Interpolate the local minima and the local maxima with cubic spline to construct the lower and the upper envelopes of h ij (t). d.
Compute the instantaneous mean of the lower and upper envelopes m ij (t) e. Let If h ij (t) satisfies the stop criteria for IMF sifting, then set the i-th IMF c i (t) = h ij (t). Otherwise, let j = j + 1, return to step 2b.

3.
Let If r i (t) satisfies the stop criteria for EMD, then set r i (t) as the residue, and terminate the EMD process. Otherwise, let i = i + 1, return to Step 2.
To illustrate the method, a simulated signal was used, which shows the vibrations of a gear transmission operating at a variable rotational speed [9].
where m (0,1, . . . ,M) is the meshing harmonic number, A m the amplitude at the m-th harmonic frequency f m (i.e., f m = m·N· f s (t), where N is the tooth number and f s is the shaft rotation frequency), t time vector (with sampling time dt), β m the initial phase, z(t) the impact-induced resonant vibration, a m (t) and θ m (t) are the modified amplitude and phase modulation functions at the m-th harmonic, respectively, φ(t) the white noise.
The following model parameters were adopted for the simulation: is a random walk around the value 1. z(t) is described as a convolution of the resonant response to impact during one rotation and the pulse train is the sequence of pulses occurring once per revolution.
Signal-to-noise-ratio, SNR, is calculated as where s i is i-th amplitude of the signal and w i i-th amplitude of the noise. For a simulated signal, SNR = 4.6 dB. The main components of signal power are visible impacts. The signal includes both quasi-harmonic components and impulse responses related to the damage of the gear. The signal presented by the Formula (12) is shown in Figure 3a, together with the following intristic mode functions from 1 to 5, which are the result of decomposition. Figure 3b shows the multifractal spectra for waveforms illustrated in Figure 3a.
The first empirical mode reveals the high-frequency impacts. This high-frequency mode increases in energy as the wear or damage progresses. The next modes show the resonance in the vibration signal and rotating frequency. The first intristic mode function will be used for further analysis. The spectrum for IMF1 is broader and more regular than the spectrum of the raw signal and the other IMF, while the IMF2 is a monofractal.
where is i-th amplitude of the signal and i-th amplitude of the noise. For a simulated signal, SNR = 4.6 dB. The main components of signal power are visible impacts.
The signal includes both quasi-harmonic components and impulse responses related to the damage of the gear. The signal presented by the Formula (12) is shown in Figure 3a, together with the following intristic mode functions from 1 to 5, which are the result of decomposition. Figure 3b shows the multifractal spectra for waveforms illustrated in   The first empirical mode reveals the high-frequency impacts. This high-frequency mode increases in energy as the wear or damage progresses. The next modes show the resonance in the vibration signal and rotating frequency. The first intristic mode function will be used for further analysis. The spectrum for IMF1 is broader and more regular than the spectrum of the raw signal and the other IMF, while the IMF2 is a monofractal. Figure 4 shows a diagram of the EMD-WLMF method. From a dynamic signal using EMD analysis is extracting the first IMF. Then, using the WLMF algorithm, the signal parameters h and D(h) are determined, on the basis of which a multifractal spectrum can be drawn. Based on the characteristics of this spectrum, the operational state of the rotating machine is classified.  Figure 4 shows a diagram of the EMD-WLMF method. From a dynamic signal using EMD analysis is extracting the first IMF. Then, using the WLMF algorithm, the signal parameters h and D(h) are determined, on the basis of which a multifractal spectrum can be drawn. Based on the characteristics of this spectrum, the operational state of the rotating machine is classified. The first empirical mode reveals the high-frequency impacts. This high-frequen mode increases in energy as the wear or damage progresses. The next modes show t resonance in the vibration signal and rotating frequency. The first intristic mode functi will be used for further analysis. The spectrum for IMF1 is broader and more regular th the spectrum of the raw signal and the other IMF, while the IMF2 is a monofractal. Figure 4 shows a diagram of the EMD-WLMF method. From a dynamic signal usi EMD analysis is extracting the first IMF. Then, using the WLMF algorithm, the signal p rameters h and D(h) are determined, on the basis of which a multifractal spectrum can drawn. Based on the characteristics of this spectrum, the operational state of the rotati machine is classified.

Application of the EMD-WLMF Algorithm in Diagnostics of Rotating Machines
The method presented in Section 2.3 was tested on the vibration signals recorded for the gear transmission on a laboratory stand and on a real object, which is a passenger car while driving.

Gear Transmission Vibration Signal Analysis on a Laboratory Stand
Measurements were carried out on a demonstration stand ( Figure 5). The influence of assembly errors and gear teeth wear on vibrations was investigated. The electric motor enables speed control in the range of 100-3000 rpm (without load). The load is pressure regulated with an overflow valve up to 5 MPa. Vibration acceleration was measured with an accelerometer bolted to the gear bearing housing in the vertical direction. The optimal backlash was set to 0.1 mm.

Application of the EMD-WLMF Algorithm in Diagnostics of Rotating Mach
The method presented in Section 2.3 was tested on the vibration signals reco the gear transmission on a laboratory stand and on a real object, which is a passe while driving.

Gear Transmission Vibration Signal Analysis on a Laboratory Stand
Measurements were carried out on a demonstration stand ( Figure 5). The i of assembly errors and gear teeth wear on vibrations was investigated. The electr enables speed control in the range of 100-3000 rpm (without load). The load is regulated with an overflow valve up to 5 MPa. Vibration acceleration was measu an accelerometer bolted to the gear bearing housing in the vertical direction. The backlash was set to 0.1 mm.
Compared to previous papers [48], the experiment was extended with addit termediate operating states and the analysis method was modified. Measurements were carried out for the following five states:  fault-free (new gears, the optimal backlash, parallel shaft location);  new gears and misalignment by an angle up to 1/3° (two cases);  new gears and increased backlash +0.1 mm up to +0.3 mm (three cases);  worn teeth (three cases);  worn teeth and increased backlash +0.3 mm (two cases).
Vibration acceleration signals were recorded for a rotational speed of about 1 and a load of 12%-pressure 0.6 MPa. Each sample included a time series with a n = 10,000, recorded at a sampling frequency of 10 kHz.
The representative waveforms of the vibration acceleration signal for each o states are shown in Figure 6. The first IMF is marked on them, which is taken into in further analysis of the signal. The signals presented in the graph represent a si degree of damage, and their characteristic features are visible. Compared to previous papers [48], the experiment was extended with additional intermediate operating states and the analysis method was modified.
Measurements were carried out for the following five states: • fault-free (new gears, the optimal backlash, parallel shaft location); • new gears and misalignment by an angle up to 1/3 • (two cases); • new gears and increased backlash +0.1 mm up to +0.3 mm (three cases); • worn teeth (three cases); • worn teeth and increased backlash +0.3 mm (two cases).
Vibration acceleration signals were recorded for a rotational speed of about 1050 rpm and a load of 12%-pressure 0.6 MPa. Each sample included a time series with a length of n = 10,000, recorded at a sampling frequency of 10 kHz.
The representative waveforms of the vibration acceleration signal for each of the five states are shown in Figure 6. The first IMF is marked on them, which is taken into account in further analysis of the signal. The signals presented in the graph represent a significant degree of damage, and their characteristic features are visible. Sensors 2021, 21, x FOR PEER REVIEW 10 of 21 Twenty records for each condition were used for further analysis, each 1 s long (10,000 samples) for a total of 220 records. Each record was subjected to the EMD-WLMF analysis. Figure 7 shows the multifractal spectra for various operational states of the gears. The examples show up to two cases per state to increase the visibility of the image.  Twenty records for each condition were used for further analysis, each 1 s long (10,000 samples) for a total of 220 records. Each record was subjected to the EMD-WLMF analysis. Figure 7 shows the multifractal spectra for various operational states of the gears. The examples show up to two cases per state to increase the visibility of the image. Twenty records for each condition were used for further analysis, each 1 s long (10,000 samples) for a total of 220 records. Each record was subjected to the EMD-WLMF analysis. Figure 7 shows the multifractal spectra for various operational states of the gears. The examples show up to two cases per state to increase the visibility of the image.  The least multifractal character is for the fault-free state. As the damage develops, the width of the spectrum increases. The spectral peak is also shifted. For the base state, the spectrum asymmetry index is close to zero. For the tested damages, the absolute value of the asymmetry index increases, and for the wear condition of the gears it has the opposite sign than for the increased backlash. For mixed states of damage, the spectrum is characterized by a large width and indices of asymmetry varying within wide limits.
Three parameters were selected for the diagnosis of the gear condition: the width of the multifractal spectrum, the asymmetry index, and the singularity of the greatest dimension. The list of these parameters is presented in the scatter plot (Figure 8a). The least multifractal character is for the fault-free state. As the damage develops, the width of the spectrum increases. The spectral peak is also shifted. For the base state, the spectrum asymmetry index is close to zero. For the tested damages, the absolute value of the asymmetry index increases, and for the wear condition of the gears it has the opposite sign than for the increased backlash. For mixed states of damage, the spectrum is characterized by a large width and indices of asymmetry varying within wide limits.
Three parameters were selected for the diagnosis of the gear condition: the width of the multifractal spectrum, the asymmetry index, and the singularity of the greatest dimension. The list of these parameters is presented in the scatter plot (Figure 8a). As the damage develops, the level of multifractality of the signal increases and the spectrum is wider. The asymmetry index for the 'base state' is close to zero, and not much As the damage develops, the level of multifractality of the signal increases and the spectrum is wider. The asymmetry index for the 'base state' is close to zero, and not much higher for the 'misalignment' state. The states of increased backlash and wear are characterized by the opposite sign of the spectrum asymmetry index and an increase in its absolute value. The spectrum shift indicator also changes for different operating states.
Any known machine learning method-e.g., k-nearest neighbours (k-NN), support vector machines (SVM), or neural network (NN)-can be used to classify the damage [25]. The fault classification was performed for a newly generated test set of data (240 records) with support vector machines (SVM). The confusion matrix is shown in Figure 8b. The average true positive rate for the test set is 95.16%. The faulty recognition of the state was for 'fault-free' and initial 'misalignment'. The remaining states found to be in error are 'increased backlash' and 'worn teeth and increased backlash'. The tests performed on the laboratory stand brought the desired effect. The influence of damage development on the parameters of the multifractal spectrum was presented. However, the stand lacks additional disturbances that accompany the operation of real machines.

Analysis of the Transmission Vibration Signal in a Car
The tests were carried out on a Punto five-speed gearbox (Figure 9a). The experiment was carried out during road tests for various rotational speeds and loads. The experiment used the recorded signals of acceleration of vibrations of the gearbox housing and additional signals from the crankshaft position and throttle position sensor.
'increased backlash' and 'worn teeth and increased backlash'. The tests performed on the laboratory stand brought the desired effect. The influence of damage development on the parameters of the multifractal spectrum was presented. However, the stand lacks additional disturbances that accompany the operation of real machines.

Analysis of the Transmission Vibration Signal in a Car
The tests were carried out on a Punto five-speed gearbox (Figure 9a). The experiment was carried out during road tests for various rotational speeds and loads. The experiment used the recorded signals of acceleration of vibrations of the gearbox housing and additional signals from the crankshaft position and throttle position sensor.
The Signals of a duration of about 1 min were recorded while driving at the most preset speed possible. Additional signals enabled signal synchronization and measurement of rotational speed and load.
The active experiment consisted in registering signals under simulated mechanical damage to the gearbox, mimicking those often diagnosed in car repair shops. Such damage includes, among others, wear of the gearbox teeth.
In order to investigate the influence of gearbox teeth wear on the vibration signal, an active experiment was performed for the following conditions:   Signals of a duration of about 1 min were recorded while driving at the most preset speed possible. Additional signals enabled signal synchronization and measurement of rotational speed and load.
The active experiment consisted in registering signals under simulated mechanical damage to the gearbox, mimicking those often diagnosed in car repair shops. Such damage includes, among others, wear of the gearbox teeth.
In order to investigate the influence of gearbox teeth wear on the vibration signal, an active experiment was performed for the following conditions:   The EMD-WLMF algorithm was tested for 80 records, 20 records for each state. Fi ure 11a presents sample multifractal spectra for four states, and Figure 11b presents mu tifractal spectra for the same states using the only WLMF algorithm. The EMD-WLMF algorithm was tested for 80 records, 20 records for each state. Figure 11a presents sample multifractal spectra for four states, and Figure 11b presents multifractal spectra for the same states using the only WLMF algorithm. By comparing the multifractal spectra for the first IMF of the vibration signal and the raw of the vibration signal, it can be concluded that the introduction of the EMD enables the identification of tooth wear at an early stage. In Figure 10b, the spectrum at initial wear almost does not differ from the base state and the spectrum for a new pair of gears. The EMD-WLMF analysis showed that the spectrum for the new pair of gears ('after repair') is the most symmetrical. The spectrum for the 'initial state' is narrower, but shows a shift and asymmetry. As the wear develops, the spectrum becomes wider, shifted along the h(q) axis, and shows asymmetry. The illustration of spectrum features for all 80 records is shown in Figure 12a.
The scatter plot for more records confirms that the spectrum width increases as the wear develops. The asymmetry coefficient for all states has a large spread, while for the new pair of teeth ('after repair') the spectrum is the most symmetrical. As the wear of the teeth increases, the spectrum shifts to the right. The average trend of damage development was marked. For a new pair of teeth ('after repair'), the development of the damage may follow a 'different path' on the scatter plot, but the increasing trend of the width indicators and spectrum shift should be analogous.
The fault classification was performed for a newly generated test set of data (120 records) with support vector machines (SVM). The confusion matrix is shown in Figure 12b. The average true positive rate for the test set is 94.17%. The faulty recognition of the state was for 'fault-free' and 'initial wear'. Due to the dispersion of the data, the diagnosis can be confirmed on the basis of the average trend shown in Figure 12a. By comparing the multifractal spectra for the first IMF of the vibration signal and the raw of the vibration signal, it can be concluded that the introduction of the EMD enables the identification of tooth wear at an early stage. In Figure 10b, the spectrum at initial wear almost does not differ from the base state and the spectrum for a new pair of gears. The EMD-WLMF analysis showed that the spectrum for the new pair of gears ('after repair') is the most symmetrical. The spectrum for the 'initial state' is narrower, but shows a shift and asymmetry. As the wear develops, the spectrum becomes wider, shifted along the h(q) axis, and shows asymmetry. The illustration of spectrum features for all 80 records is shown in Figure 12a.
The scatter plot for more records confirms that the spectrum width increases as the wear develops. The asymmetry coefficient for all states has a large spread, while for the new pair of teeth ('after repair') the spectrum is the most symmetrical. As the wear of the teeth increases, the spectrum shifts to the right. The average trend of damage development was marked. For a new pair of teeth ('after repair'), the development of the damage may follow a 'different path' on the scatter plot, but the increasing trend of the width indicators and spectrum shift should be analogous.
The fault classification was performed for a newly generated test set of data (120 records) with support vector machines (SVM). The confusion matrix is shown in Figure 12b. The average true positive rate for the test set is 94.17%. The faulty recognition of the state was for 'fault-free' and 'initial wear'. Due to the dispersion of the data, the diagnosis can be confirmed on the basis of the average trend shown in Figure 12a.

Analysis of the Vibration Signal of the Internal Combustion Engine Head
The next example shows the analysis of the vibration signal of the engine head signal recorded under operating conditions.
The tests were carried out on the drive system of the Fiat Punto car, with a fourcylinder 1.2 spark-ignition engine. The experiment was carried out during road tests for various rotational speeds and loads. The following dynamic signals were recorded:

Analysis of the Vibration Signal of the Internal Combustion Engine Head
The next example shows the analysis of the vibration signal of the engine head signal recorded under operating conditions.
The tests were carried out on the drive system of the Fiat Punto car, with a fourcylinder 1.2 spark-ignition engine. The experiment was carried out during road tests for various rotational speeds and loads. The following dynamic signals were recorded: Signals of approximately 1 min duration were recorded while driving at a constant speed. Additional signals made it possible to identify engine cycles, injection and ignition timing, and valve timing.
A leak in the piston-cylinder system may be caused by burnout of the exhaust valve. During the tests, such damage was simulated by cutting the valve plug at a length of approx. 3 mm (defect 1) and 6 mm (defect 2)- Figure 13.  Signals of approximately 1 min duration were recorded while driving at a constant speed. Additional signals made it possible to identify engine cycles, injection and ignition timing, and valve timing.
A leak in the piston-cylinder system may be caused by burnout of the exhaust valve. During the tests, such damage was simulated by cutting the valve plug at a length of approx. 3 mm (defect 1) and 6 mm (defect 2)- Figure 13. The process of generating vibrations in an internal combustion engine is very complex. The measured vibrations are a combination of periodic waves related to the operation of rotating elements and responses to impulse excitations related to reciprocating motion, as well as excitations caused by gas pressure. Strong transients in the vibroacoustic signal come from the work of the intake and exhaust valves, injectors, the combustion process, and piston strokes against the cylinder sleeve. Additional impulse excitation can be a response to mechanical failures not detected by the on-board diagnostic (OBD) system [49,50].
The time series of vibration accelerations after determining the first IMF for four engine operation cycles and three operating states are shown in Figure 14. Figure 15 presents exemplary multifractal spectra for three operating states using the WLMF and EMD-WLMF algorithms. Due to the presence of transients in the signal, a large share of responses to impulse excitations-such as valve closing, ignition, piston hitting the cylinder, etc.-of the multifractal spectrum after EMD application does not differ much from the spectrum without EMD filtration. As a result of the damage development, the spectrum becomes wider and its maximum value for q = 0 is shifted to the right. The process of generating vibrations in an internal combustion engine is very complex. The measured vibrations are a combination of periodic waves related to the operation of rotating elements and responses to impulse excitations related to reciprocating motion, as well as excitations caused by gas pressure. Strong transients in the vibroacoustic signal come from the work of the intake and exhaust valves, injectors, the combustion process, and piston strokes against the cylinder sleeve. Additional impulse excitation can be a response to mechanical failures not detected by the on-board diagnostic (OBD) system [49,50].
The time series of vibration accelerations after determining the first IMF for four engine operation cycles and three operating states are shown in Figure 14.   Figure 15 presents exemplary multifractal spectra for three operating states using the WLMF and EMD-WLMF algorithms. Due to the presence of transients in the signal, a large share of responses to impulse excitations-such as valve closing, ignition, piston hitting the cylinder, etc.-of the multifractal spectrum after EMD application does not differ much from the spectrum without EMD filtration. As a result of the damage development, the spectrum becomes wider and its maximum value for q = 0 is shifted to the right.   Figure 16 shows the scatter plot for the three operational states, each with 20 records. As the damage develops, the values of the 'width' and 'shift' parameters increase. The 'asymmetry' parameter is characterized by the largest dispersion of values, although it also increases taking into account the average values for the entire set. The average trend for damage development is marked.  Figure 16 shows the scatter plot for the three operational states, each with 20 records. As the damage develops, the values of the 'width' and 'shift' parameters increase. The 'asymmetry' parameter is characterized by the largest dispersion of values, although it also increases taking into account the average values for the entire set. The average trend for damage development is marked.
The fault classification was performed for a newly generated test set of data (120 records) with support vector machines (SVM). The confusion matrix is shown in Figure 16b. The average true positive rate for the test set is 92.5%. The faulty recognition of the state was for 'fault-free' and 'defect 1'. Due to the dispersion of the data, the diagnosis can be confirmed on the basis of the average trend shown in Figure 16a.
The fault classification was performed for a newly generated test set of data (120 records) with support vector machines (SVM). The confusion matrix is shown in Figure 16b. The average true positive rate for the test set is 92.5%. The faulty recognition of the state was for 'fault-free' and 'defect 1'. Due to the dispersion of the data, the diagnosis can be confirmed on the basis of the average trend shown in Figure 16a.

Discussion
The article uses the method of wavelet leaders' multifractal analysis dedicated to the first IMF of the vibration signal of rotating machines. The presented method can be applied to dynamic signals in non-stationary conditions due to the use of the EMD and wavelet analysis methods. Previous works [32] used the EMD analysis and the MF-DFA algorithm to diagnose the damage. However, the WLMF algorithm shows lower computational costs, numerical stability, and high versatility in terms of real signals compared to MF-DFA. The method corresponds to the need to find a relatively fast algorithm for diagnosing the condition of machines in non-stationary conditions. It is a statistical method. It can be an alternative to methods based on a set of signal features. Failure classification can be done using one of the machine learning methods, such as SVM, k-NN, or a neural network based on the features of the multifractal spectrum.

Discussion
The article uses the method of wavelet leaders' multifractal analysis dedicated to the first IMF of the vibration signal of rotating machines. The presented method can be applied to dynamic signals in non-stationary conditions due to the use of the EMD and wavelet analysis methods. Previous works [32] used the EMD analysis and the MF-DFA algorithm to diagnose the damage. However, the WLMF algorithm shows lower computational costs, numerical stability, and high versatility in terms of real signals compared to MF-DFA. The method corresponds to the need to find a relatively fast algorithm for diagnosing the condition of machines in non-stationary conditions. It is a statistical method. It can be an alternative to methods based on a set of signal features. Failure classification can be done using one of the machine learning methods, such as SVM, k-NN, or a neural network based on the features of the multifractal spectrum.
The method has been tested for signals recorded on a laboratory stand and for real objects. Damage was distinguished based on three features of the multifractal spectrum, named here: width, asymmetry, and shift. The tests performed on the laboratory stand brought a good effect. Multifractal analysis of only the first IMF allowed to filter out harmonics that do not show multifractality features. The development of damage to the gear transmission was reflected in the examined features of the multifractal spectrum. However, the stand lacks additional disturbances that accompany the operation of real machines.
The next two investigated cases are the diagnosis of mechanical damage to the vehicle while driving, such as gearbox wear and burnout of the exhaust valves of the internal combustion engine. The signals were recorded in non-stationary conditions. Due to the fluctuations in rotational speed and load, which could affect the dynamics of the signal, the dispersion of the features of the multifractal spectrum was much greater than that observed in the tests on the laboratory stand. Due to the dispersion of the features of the multifractal spectrum, the average value of these features can be used to diagnose the damage. By comparing the multifractal spectra for the first IMF of the vibration signal and the raw of the vibration signal of the gearbox, it can be concluded that the introduction of the EMD enables the identification of tooth wear at an early stage. Due to the presence of transients of the engine head vibrations signal, a large share of responses to impulse excitations-such as valve closing, ignition, piston hitting the cylinder, etc.-the multifractal spectrum after EMD application does not differ much from the spectrum without EMD filtration.
Like any method, this one also has its limitations. The main limitation is the multifractality of the signal, so the signal must contain impulse excitations caused by damage, as in a toothed gear. If the damage causes an increase in the harmonic component, the change in the multifractal spectrum will be hardly noticeable and then other methods-e.g., order tracking or the AR reference model-will bring better results.