An Enhanced Intrinsic Time-Scale Decomposition Method Based on Adaptive Lévy Noise and Its Application in Bearing Fault Diagnosis

When early failures in rolling bearings occur, we need to be able to extract weak fault characteristic frequencies under the influence of strong noise and then perform fault diagnosis. Therefore, a new method is proposed: complete ensemble intrinsic time-scale decomposition with adaptive Lévy noise (CEITDALN). This method solves the problem of the traditional complete ensemble intrinsic time-scale decomposition with adaptive noise (CEITDAN) method not being able to filter nonwhite noise in measured vibration signal noise. Therefore, in the method proposed in this paper, a noise model in the form of parameter-adjusted noise is used to replace traditional white noise. We used an optimization algorithm to adaptively adjust the model parameters, reducing the impact of nonwhite noise on the feature frequency extraction. The experimental results for the simulation and vibration signals of rolling bearings showed that the CEITDALN method could extract weak fault features more effectively than traditional methods.


Introduction
Rolling bearings are the most commonly used and extremely vulnerable core supporting components in rotating machinery. If they fail, they directly affect the safe and stable operation of the whole unit [1][2][3][4]. Fault detection through analyzing the vibration signals of rolling bearings is currently one of the most common methods used for bearing analysis. By installing an acceleration sensor near the bearing housing, the vibration signal generated by the rotating rolling bearing is collected and analyzed to determine the condition of the bearing. Affected by many environmental factors (such as vibration caused by the non-linearity of bearing stiffness, vibration caused by the corrugation of machining surface, vibration caused by the eccentricity of bearings, etc.) during the signal acquisition process, the fault signal of rotating machinery usually contains a large amount of noise. Meanwhile, the fault information is difficult to identify due to its unobvious early default characteristics. Additionally, the rotating machinery fault signal usually shows non-smooth, non-linear characteristics, because the dynamic signal will be non-smooth when a bearing fault occurs during operation, such as friction, cracks, oil film vortex, oil film oscillation, etc. Therefore, non-smoothness can characterize the presence of faults. Meanwhile, when a fault occurs and develops, the damping, stiffness, and elastic force of the bearing change, and the dynamic signal also exhibits nonlinear characteristics. The traditional time-domain analysis methods and Fourier transform-based analysis methods can capture only the overall statistical characteristics of the signal, but not the time-frequency detail characteristics of the signal, and cannot handle nonlinear, non-stationary signals. With the continuous development of signal processing technology, modern time-frequency analysis methods have made great progress in their ability to extract fault features from non-stationary, non-linear signals [5][6][7].
There are many weak signal processing methods available, and wavelet transform (WT) is often effective in analyzing non-stationary and non-linear signals [8,9]. However, WT has some shortcomings, such as the choice of certain wavelet basis functions. In recent years, some adaptive signal decomposition methods have been applied to the analysis of the vibration signals of bearings, such as empirical mode decomposition (EMD) [10] and local mean decomposition (LMD) [11]. Intrinsic time-scale decomposition (ITD) is another adaptive signal decomposition method that was proposed by Frei and Osorio [12] and first applied in medical diagnosis. ITD uses linear interpolation to construct the baseline signal, instead of averaging the envelope using cubic spline interpolation as in EMD. As a result, ITD has a clear advantage in terms of operational efficiency. Each proper rotational component (PRC) retains almost all of the polar information of the original signal. As a result, ITD is now widely used to extract the fault characteristics of rolling bearings, as the ITD process happens to be a demodulation process. On the basis of the above advantages, scholars introduced the ITD method from the medical field into mechanical signal fault diagnosis. For example, Lin and Chang proposed a rolling bearing fault diagnosis method based on an enhanced kurtosis spectrum and intrinsic time-scale decomposition [13]. Duan and Yao et al. proposed to apply intrinsic time-scale decomposition to the fault diagnosis of gearboxes under variable operating conditions [14]. Xing and Qu et al. proposed a gear fault diagnosis method with variable working conditions based on intrinsic time-scale and singular value decomposition [15]. Zhang and Liu et al. proposed a diesel engine fault diagnosis method with complete ensemble intrinsic timescale decomposition [16]. Hu and Xiang proposed the application of the ensemble intrinsic time-scale decomposition of fan gear [17]. Tong and Cao et al. proposed an improved intrinsic time-scale decomposition with wavelet packet transform and singular-value decomposition to perform fault diagnosis on rolling bearings [18]. Liu and Zhang et al. proposed using intrinsic time-scale decomposition to diagnose diesel-engine faults [19]. Bi and Ma et al. proposed detecting gasoline-engine knock using the improved intrinsic time-scale decomposition of complete integration [20]. Yu and Liu proposed sparse coding based on intrinsic time-scale decomposition to diagnose weak bearing faults [21]. Yuan and Peng proposed smooth intrinsic time-scale decomposition for the fault diagnosis of rolling bearings [22]. Lei and Zhou et al. used intrinsic time-scale decomposition to monitor tool wear during milling [23]. Ma and Zhan et al. proposed complete ensemble intrinsic time-scale decomposition with adaptive noise (CEITDAN), which was applied to the feature extraction of rolling bearings [24]. According to the research findings of previous ITD methods in the process of signal decomposition based on ITD research. Through the analysis of previous studies, their proposed methods improve the problem of the modal blending of ITD methods, expand the application scope of ITD methods, and improve the effect of feature extraction, but there are still certain shortcomings, such as the fact that a certain amount of nonwhite noise components are still retained in the signals after the noise reduction in these methods. They did not find an effective method to filter out the nonwhite noise components in the signal (e.g., industrial frequency interference, colored noise caused by a drive motor, etc.), which can be contained between the feature frequencies and interfere with the extraction of the feature frequencies, making it difficult to locate them accurately. Additionally, the existence of these noises will cause the calculation of feature parameters such as the entropy value to have large errors, thus reducing the fault diagnosis accuracy.
For the current bearing fault signal processing method for use after noise reduction, where there is still a residual nonwhite noise component interference problem, this paper is inspired by the literature [25] and aims to propose a fault diagnosis method for the signal-to-noise ratio (SNR) −25 dB. This paper aims to improve the noise reduction after the nonwhite noise interference fault diagnosis, reducing the problem of fault diagnosis accuracy. Therefore, this article proposes a method that involves using variable-parameter analog noise (i.e., Lévy noise) to replace traditional white noise. This can effectively filter out noise in the decomposition process and improve the influence of the nonwhite noise part of noise on feature frequency extraction; this is called complete ensemble intrinsic time-scale decomposition with adaptive Lévy noise (CEITDALN). Lévy noise is a digital noise model that contains four parameters. Changes in the parameters can cause dramatic changes in the noise model [26]. Therefore, in order to avoid the difficulty of artificially setting parameters, this paper proposes an improved coyote algorithm to adaptively adjust the noise intensity coefficient and Lévy noise parameters in the CEITDALN method. The energy of the maximal center frequency band of the kurtosis spectrum was used as the objective function of the optimization algorithm to achieve the best feature extraction.
The rest of the article is organized as follows: Section 2 introduces related works such as CEITDAN and the theory of the coyote optimization algorithm. Section 3 introduces CEITDALN and the improved coyote optimization algorithm. Section 4 uses analog signals to verify the effectiveness of the proposed method and introduces the engineering application platform. Through the above platform, the effectiveness of the proposed method is further proved. Section 5 concludes this article.

CEITDAN
CEITDAN's decomposition method superimposes a component decomposed by ITD from white noise with a certain signal-to-noise ratio (we add a white noise whose standard deviation is normalized and whose noise amplitude can be adjusted adaptively, thus ensuring that the SNR of each decomposition stage is fixed) onto the original signal and uses the average residue component obtained by the ITD to obtain the first rotation component. It then subtracts the first rotation component from the original signal.
The residual component is obtained, white noise is added to the residual component, and the same operation is performed until the residual component is a monotonic function or there are less than three extreme points. A proper rotation component can accurately define the instantaneous information of the signal and is repeated several times in this manner.
Step 1. The first round of noise addition: ω i n (t)(i = 1, 2, · · · , I) is the first-order proper rotation component decomposed by ITD from white noise with a certain SNR (here, we use the coefficient β 0 = ε 0 std(x)/std ω (i) to control the SNR of each unchanged decomposition stage), where I is the number of noise additions superimposed onto the original signal x(t). The average value is used as the residue component of the method; then, the first-order component of the origin signal is obtained as follows: n is the corresponding order of ITD white noise. β 0 is used as β 0 = ε 0 std(x)/std(ω (i) ), β 0 is the noise-adding amplitude coefficient, ε i is the quotient of standard deviation, L k [A(t)] represents the residue component obtained by the ITD of signal A(t), and A(t) is the signal with white noise. The proper components at this time are as follows: Step 2. The second round of noise addition: The second-order PRCs of white noise ω i (t) are superimposed onto r 1 (t), and the first-order components of the mixed signal are The second-order proper component obtained at this time is as follows: Step 3. TheM-th round of noise addition: The (m − 1)-th-order PRCs of white noise ω i (t) are superimposed onto the remaining terms r m−1 (t). Then, the first-order component of the signal with adaptive white noise is obtained by ITD, and the average value is taken as the m-th-order residue component of the method, as follows: The PRC obtained at this point is: Step 4. Repeat Steps 1 to 3 several times until the residual component is a monotonic function or the number of extreme points in this residual component is fewer than three. The final remaining term is as follows: At this point, the entire CEITDAN decomposition process ends.
The final shifted terms are as follows: The flowchart of this method is shown in Figure 1:

Group Optimization Algorithm Theory
The coyote optimization algorithm (COA) was proposed by Pierezan et al. [27] in 2018. It is a new intelligent optimization algorithm that simulates the phenomena of coyotes living in groups, including their growth, life, and death and being driven away and accepted by groups. It obtains better optimization results in the optimization of benchmark functions. The COA divides a population into several subgroups through random grouping, determines the alpha wolf and cultural trend of the subgroup, and randomly selects two coyotes. These four factors influence coyote growth. It reforms the growth process through coyote social adaptability. The birth of a coyote is affected by two randomly selected parent coyotes and environmental variation. In terms of social adaptability, if a newborn coyote is better than an old and poor coyote, the old coyote dies. Otherwise, the newborn coyote dies. According to a certain probability among subgroups, some coyotes are driven out of the group and accepted by other groups, thereby changing the coyote grouping status through the continuous evolution of this process of growth, death, driving away, and acceptance. It obtains the coyote that is most suitable for the social environment as the best solution to an optimization problem.

Group Optimization Algorithm Theory
The coyote optimization algorithm (COA) was proposed by Pierezan et al. [27] in 2018. It is a new intelligent optimization algorithm that simulates the phenomena of coyotes living in groups, including their growth, life, and death and being driven away and accepted by groups. It obtains better optimization results in the optimization of benchmark functions. The COA divides a population into several subgroups through random grouping, determines the alpha wolf and cultural trend of the subgroup, and randomly selects two coyotes. These four factors influence coyote growth. It reforms the growth process through coyote social adaptability. The birth of a coyote is affected by two randomly selected parent coyotes and environmental variation. In terms of social adaptability, if a newborn coyote is better than an old and poor coyote, the old coyote dies. Otherwise, the newborn coyote dies. According to a certain probability among subgroups, some coyotes are driven out of the group and accepted by other groups, thereby changing the coyote grouping status through the continuous evolution of this process of growth, death, driving away, and acceptance. It obtains the coyote that is most suitable for the social environment as the best solution to an optimization problem.
In the COA, each coyote represents a candidate solution. Each solution vector is composed of the social state factor of the coyote. These status factors include coyote internal factors and external factors. Each state factor represents a decision variable, and D state factors constitute a solution vector containing D decision variables. Every coyote is measured in terms of its social adaptability. The COA is mainly divided into four steps: randomly initializing the coyote group and randomly grouping them, the coyotes growing up in the group, the life and death of the coyotes, and coyotes being driven away and accepted.
Step 1: Parameters such as coyote population p N , number of coyote individuals per group c N , dimension D , and termination condition nfevalMax were set.
Step 2: Randomly initialize the coyote pack; the i -th coyote individual in pack p at time t is defined as: In the COA, each coyote represents a candidate solution. Each solution vector is composed of the social state factor of the coyote. These status factors include coyote internal factors and external factors. Each state factor represents a decision variable, and D state factors constitute a solution vector containing D decision variables. Every coyote is measured in terms of its social adaptability. The COA is mainly divided into four steps: randomly initializing the coyote group and randomly grouping them, the coyotes growing up in the group, the life and death of the coyotes, and coyotes being driven away and accepted.
Step 1: Parameters such as coyote population N p , number of coyote individuals per group N c , dimension D, and termination condition n f eval Max were set.
Step 2: Randomly initialize the coyote pack; the i-th coyote individual in pack p at time t is defined as: x p,t c,j = lb j + r j ub j − lb j (9) where ub j and lb j denote the upper and lower bounds of the j-th dimensional values, respectively, and r j is a randomly generated real number in the range of [0, 1].
Step 3: Evaluation of coyote adaptability: where f it p,t i ∈ R, x p,t i implies in the coyote's adaptation to the environment (cost of the objective function f it p,t i ).
Step 4: Coyotes sometimes break away or are expelled from their original pack, forming a pack shift; the probability of the occurrence of this is defined as P e : Symmetry 2021, 13, 617 6 of 26 Step 5: Find the head wolf al pha p,t in the current pack and calculate the cultural trend cult p,t of the current coyote pack: where O p,t Nc +1 2 ,j denotes the median of the j-th dimensional variable for all coyotes within the p group at moment t when N c is odd.
Step 6: Simulating the birth-at-death practice in genetics, the age (in years) of the coyote was written as age p,t c . The birth of a new coyote ( pup p,t ) was written as a combination of the social status of both parents (randomly selected) plus environmental influences: where m 1 and m 2 are random coyotes from within the p wolf pack, j 1 and j 2 are two random dimensions of the problem, and R j and rand j are both random numbers within [0, 1] generated by uniform probability. The discrete probability (P s ) and the association probability (P a ) affect the cultural diversity of individuals in the coyote pack. P s and P a are defined as: Assuming that ω indicates that the coyotes in the group are less adaptive than the pups and that φ represents the number of coyotes in the current group, if φ is 1 and ω holds-i.e., if the number of coyotes in a group is 1 and the pups are more adaptive than only one coyote-the pups survive and the oldest coyote in the group dies; in the remaining cases, the pups die.
Step 7: Calculate the effect of the head wolf and group cultural trends on the renewal of individuals within the coyote pack corresponding to the current moments δ 1 and δ 2 : where cr 1 and cr 2 represent the random coyotes in the current pack, respectively.
Step 8: All coyote individuals in the coyote pack were updated in turn to obtain new coyote individuals new_x p,t i , the size of the adaptation of the new coyote to the original coyote was selected by merit, and the optimal coyote x p,t+1 i was retained: where r 1 and r 2 are uniform probability generated real numbers in the range of [0, 1] representing the magnitude of the weight of individual coyotes influenced by the cultural trend of al pha wolves versus the group.
Step 9: Update the ages of individual coyotes by simulating the growth of individuals over time.
Step 10: Judge the termination condition, and if it is reached, output the social state of the coyote with the best adaptation ability; otherwise, return to Step 3 to continue.
The above steps show that COA has the following advantages: (a) COA has a good search model and framework. The coyotes are randomly divided into multiple subgroups. After all the coyotes grow up, cultural exchanges are carried out through repelling and acceptance. Compared with algorithms such as particle swarm optimization (PSO), this search model and framework has stronger exploration capabilities. (b) COA guides the growth of coyotes through wolves and cultural trends. The algorithm has a strong local search ability. (c) The emergence of newborn COA coyotes comes from the combined effect of two randomly selected parent coyotes and random variations in the social environment, meaning that the algorithm has a certain global search capability. (d) COA updates each coyote in the group. Compared with the similar structure of SFLA, the COA update method is simple. (e) COA is randomly grouped after initialization. The coyotes in the group are randomly driven out and accepted, so that there is an exchange of information between groups.

Improved CEITDAN Based on Lévy Noise
The core idea of the CEITDAN algorithm is to add groups of white noise. During the decomposition process, the noise component in the noisy signal is cancelled out to achieve the purpose of noise reduction. However, the signal noise obtained by actual engineering is more complicated. There is both white and colored noise. Therefore, when white noise in the noisy signal does not exist as the main noise component, it is still difficult to extract the characteristic frequency after filtering by this method. For this reason, this article proposes a new noise addition. Instead of a single white noise, digital noise is added that can be adjusted to different noise types by changing parameters, improving the filtering effect, named complete ensemble intrinsic timescale decomposition with adaptive Lévy noise.
Lévy noise is also called stable noise [28]. Since the distribution and probability density functions of Lévy noise do not have explicit expressions, the characteristic function of the distribution of Lévy noise is as follows: where α ∈ ( 0, 2] is the characteristic index. The smaller α is, the more obvious the noise impulse and tailing characteristics are. α = 1 is the Cauchy distribution and α = 2 is the Gaussian distribution; β ∈ [−1, 1] is the symmetry parameter that determines the symmetry of the noise distribution. When β = 0, the distribution is symmetrical; σ ∈ [ 0, +∞) is the dispersion coefficient that represents the range of the Lévy distribution-that is, the degree of dispersion. µ ∈ (−∞, +∞) is the position parameter that is used to determine the center of the distribution. In this paper, on the basis of the Chambers-Mallows-Stuck (CMS) method, the Lévy stable random variable can be expressed by the following formula: When α = 1, In the above formula, V is uniformly distributed in the interval (−π/2, π/2) and W obeys an exponential distribution with an average value of 1. C α,β and D α,β,σ can be calculated according to the following formula: When α = 1, The proposed method in this paper is based on the above-mentioned Lévy noise. Therefore, the CEITDAN method was rewritten as follows.
Let the PR component calculated by pooled averaging be PR k . E j (·) is the j-th PR component, No i is the Lévy noise and satisfies a normal distribution, and x(t) is the original signal. The CEITDALN algorithm proceeds as follows: (1) Using ITD to decompose each x i (t) = x(t) + No i , followed by averaging the first residual component, the first residual component, r 1 , is: (2) By eliminating r 1 from the most primitive signal, the first PR signal PR 1 (t) is obtained as: (3) Construct the set residual signal r 1 (t) + ε 1 E 1 (No i (t)), (i = 1, 2, · · · , I) and decompose it to obtain PR 2 as: (4) When k = 2, 3, · · · , K, analogous to the above calculation, find r k (t) and then find PR k+1 : (5) The algorithm terminates when the R K (t)-pole is less than 3. The final decomposition results in: In the following, the time complexity of the proposed method is analyzed for this paper. Firstly, the problem size of the method is L. The time complexity of the method is denoted as T. According to the definition of the Cox-de Boor recursion, 2L additive operations and 3L multiplicative operations are required for all the baseline functions, and L additive operations are required for the debase lining, so the time complexity of the proposed method is where O + () and O × () denote the addition of time complexity and multiplication time complexity, respectively; n is the number of PRCs.

Improved COA Based on GBO Optimization
As the coyote algorithm was proposed recently, there are still the problems mentioned above. Therefore, the gradient-based optimizer (GBO) method proposed in [29] was introduced into the coyote optimization algorithm. The gradient search rule (GSR) in the GBO method can speed up the convergence speed of the coyote optimization algorithm and enhance the exploration trend. The local escaping operator (LEO) can effectively solve the problem of local optimal solutions. Therefore, this article proposes using the GBO method to improve the COA method and named it GBO-COA.
The GBO method comes from a combination of the gradient method and the population optimization algorithm. It explores the search domain by specifying the search Symmetry 2021, 13, 617 9 of 26 direction by Newton's method and using gradient search rules and local escape operators. In the gradient search rule, better search and better positions are achieved by controlling the movement of vectors. The method uses GSR to increase the speed of convergence. Considering that many of the optimization probabilities in Newton's gradient-based method [30] cannot be solved by differential calculations, the numerical gradient method is used instead of the direct function derivation method. The solution of GSR requires the calculation of both the first derivatives according to the Newtonian method in [31] and the Taylor series. Additionally, there is the following central difference formula [32]: According to Newton's method and Formula (31), the new (x n+1 ) can be defined as: In the GBO algorithm, adjacent locations are replaced by two other vectors in the population. The GSR is modified using an adaptive factor ρ 1 to improve the ability of GBO to balance global and local searches.
In Formulas (33) to (35), β min and β max are 0.2 and 1.2, respectively; m is the number of iterations and M is the total number of iterations. In order to balance global and local search capabilities, parameter ρ 1 changes were based on the change in sine function α. This parameter is dynamically adjusted at each iteration and can be set to a larger value at the early iterations in order to increase the population diversity. Then, as the number of iterations increases, the value is gradually reduced to speed up convergence. The generated solution should be able to explore the search space around the corresponding best solution. This helps the GSR method to get rid of any local optima, as it increases the diversity of the population to search for the best solution obtained. Thus, GSR follows as: where ε is a decimal number in the range of [0, 0.1]. x worst and x best are the worst and best solutions obtained in the optimization process, respectively. rand(1 : N) is an Ndimensional random number; r1, r2, r3, and r4 are different integers randomly selected from r4. step is the step size, affected by x best and x m r1 . The GBO algorithm uses the best space and moves the current space x n in the (x best − x) direction as the direction movement (DM). Therefore, an appropriate local search trend was generated in this process that could improve the speed of the GBO algorithm.
where rand is the random number of [0, 1] and ρ 2 is a random parameter that guarantees that each space has non-synchronization length.
In summary, the optimal space x m n is: The GBO algorithm also used LEO operators to improve the efficiency in solving complex problems. This operator can effectively change the position of x m+1 n to produce a solution X m LEO . In the x best , solutions included X1 m n and X2 m n , two random schemes x m r1 and x m r2 , and a randomly generated scheme (x k m ). The generation of solution X m LEO is as follows in Algorithm 1:

END END
In the above formula, f 1 is the uniform random number of [−1, 1]; f 2 is the random number that obeys the normal distribution; the average is 0; the standard deviation is 1; ρ 1 is the probability; and u 1 , u 2 , and u 3 are three random numbers. The formula definition is: where rand is a random number in the range of [0, 1] and µ 1 is a constant in the range of [0, 1].
The traditional COA has shown strong optimization capabilities in solving optimization problems [27]. However, the COA method also suffers from the following deficiencies. For example, the following problems exist in solving complex optimization problems. (a) The growth process of coyotes in COA is influenced by the difference between the calculated group al pha and the group breeding trend, and the convergence of the algorithm is slow when two coyotes are randomly selected in the group. (b) Guided by the group al pha wolf and the group culture trends, it may appear that the group al pha wolf and the group culture trends are locally optimal solutions, causing the algorithm to fall into a local optimum. (c) COA uses a dynamic greedy algorithm. To a certain extent, it speeds up the convergence speed, but the probability of falling into a local optimum increases. The GSR and LEO used by the GBO algorithm can improve the detection ability of COA very well. Considering the advantages and disadvantages of the COA and GBO methods, this paper proposes using the GBO method instead of the COA method to calculate the growth, life, and death of coyotes in order to improve the convergence speed and solve the local optimum stagnation problem. In contrast, the GBO method can produce smooth transitions in the exploration phase, thereby improving the speed of convergence and the ability to avoid local optima. Therefore, the combination of GBO and COA can improve the performance of COA. Figure 2 shows the flowchart of GBO-COA. methods, this paper proposes using the GBO method instead of the COA method to cal culate the growth, life, and death of coyotes in order to improve the convergence speed and solve the local optimum stagnation problem. In contrast, the GBO method can pro duce smooth transitions in the exploration phase, thereby improving the speed of conver gence and the ability to avoid local optima. Therefore, the combination of GBO and COA can improve the performance of COA. Figure 2 shows the flowchart of GBO-COA.

Proposed Method
This paper proposes an adaptive intrinsic time-scale decomposition algorithm based on Lévy noise instead of white noise. Parameters were optimized using the improved coy ote optimization algorithm, avoiding the drawbacks caused by the artificial setting of pa rameters. The method flow is shown in Figure 3.

Proposed Method
This paper proposes an adaptive intrinsic time-scale decomposition algorithm based on Lévy noise instead of white noise. Parameters were optimized using the improved coyote optimization algorithm, avoiding the drawbacks caused by the artificial setting of parameters. The method flow is shown in Figure 3.
CEITDALN includes an improved coyote optimization algorithm to adaptively adjust Lévy noise parameters and add noise intensity parameters. The most important aspect in the adaptive algorithm is defining a suitable objective function [33]. Therefore, the CEITDALN method generates multiple components in the process of decomposing the signal. If mode components are analyzed and decomposed one by one, this wastes a great amount of time. Therefore, traditional time domain parameters and other objective functions are not applicable in this method, and it is necessary to select an objective function that can judge multiple mode components after decomposition. By consulting the related literature, this paper proposes a brand-new objective function: calculating the spectral kurtosis of the input signal [34]. The center frequency band is found, and the energy of the center frequency band of the decomposed component is calculated. Maximal energy is taken as the optimization target. CEITDALN includes an improved coyote optimization algorithm to adaptively adjust Lévy noise parameters and add noise intensity parameters. The most important aspect in the adaptive algorithm is defining a suitable objective function [33]. Therefore, the CEITDALN method generates multiple components in the process of decomposing the signal. If mode components are analyzed and decomposed one by one, this wastes a great amount of time. Therefore, traditional time domain parameters and other objective functions are not applicable in this method, and it is necessary to select an objective function that can judge multiple mode components after decomposition. By consulting the related literature, this paper proposes a brand-new objective function: calculating the spectral kurtosis of the input signal [34]. The center frequency band is found, and the energy of the center frequency band of the decomposed component is calculated. Maximal energy is taken as the optimization target.

Case A: Numerical Simulation Analysis
In this section, the CEEMDAN, CEITDAN and CEITDALN methods are used to decompose the analogue signals of bearing faults to compare the effectiveness of the CEITDALN method proposed in this paper in filtering. In order to objectively evaluate the effectiveness of the CEITDALN method, two metrics-the orthogonality measure and the energy conservation measure-were introduced to quantitatively evaluate the decomposition performance of the different methods [35,36]. Theoretically, the mode components of different decomposition algorithms should all be completely orthogonal. The orthogonality index is equal to zero. However, due to computer errors and environmental disturbances, the orthogonality index cannot be zero [35]. Therefore, in practical analysis, the closer the orthogonality index is to zero, the more desirable the decomposition result is. In addition, the decomposition result should satisfy the law of energy conservation, and the energy conservation exponent should be equal to 1. However, in practice, the same from computer errors can lead to some loss of signal energy, so the value of the energy conservation exponent should be close to 1 in practical applications, indicating less signal loss in the decomposition, while when the ECI is greater than 1, it indicates the presence of spurious components [37]. Here, we compare these three methods using the

Case A: Numerical Simulation Analysis
In this section, the CEEMDAN, CEITDAN and CEITDALN methods are used to decompose the analogue signals of bearing faults to compare the effectiveness of the CEIT-DALN method proposed in this paper in filtering. In order to objectively evaluate the effectiveness of the CEITDALN method, two metrics-the orthogonality measure and the energy conservation measure-were introduced to quantitatively evaluate the decomposition performance of the different methods [35,36]. Theoretically, the mode components of different decomposition algorithms should all be completely orthogonal. The orthogonality index is equal to zero. However, due to computer errors and environmental disturbances, the orthogonality index cannot be zero [35]. Therefore, in practical analysis, the closer the orthogonality index is to zero, the more desirable the decomposition result is. In addition, the decomposition result should satisfy the law of energy conservation, and the energy conservation exponent should be equal to 1. However, in practice, the same from computer errors can lead to some loss of signal energy, so the value of the energy conservation exponent should be close to 1 in practical applications, indicating less signal loss in the decomposition, while when the ECI is greater than 1, it indicates the presence of spurious components [37]. Here, we compare these three methods using the orthogonal index (OI), root mean square error (RMSE), energy conservation index (ECI) and improved SNR (ISNR). The definitions of OI, ECI, and ISNR are as follows.
where x(t) represents the original signal, mod n (t) represents the i-th mode component decomposed by different methods, and r m (t) indicates residual error. This paper also uses the improved SNR to compare the superiority of the method proposed in this paper in the envelope spectrum, calculated as follows: where ISNR is the improved signal-to-noise ratio. NN stands for total regional power, and its calculation formula is as follows: where MM is appropriately selected according to the number of sampling points and the size of the sampling rate. In order to verify the effectiveness of the proposed method, a bearing failure model was used to simulate the generated shock signal. The structure of the simulation signal is as follows [24]: where s(t) is the periodic fault shock; f r is the frequency conversion, and f r = 86 Hz; C A is 1; C is the attenuation coefficient (C = 760); f n is the resonance frequency ( f n = 3000 Hz); and f i is the characteristic frequency of the inner ring fault ( f i = 1/T = 130 Hz). τ i is the slight sliding of the i − th shock with respect to period T, after which the sliding follows a normal distribution with a mean of 0 (standard deviation is 0.5% of the frequency conversion). n(t) is Gaussian white noise with an SNR of −35 dB, and f s is the sampling frequency ( f s = 8192 Hz). The number of analysis data points was 8192, as shown in Figure 4. Figure 5 shows the decomposition result of CEEMDAN, CEITDAN, and CEITDALN.    Figure 5 shows that the CEEMDAN method decomposed a total of 12 mode components, and CEITDAN and CEITDALN decomposed 10 rotation components. Fewer model components can save subsequent analysis time for decomposition results. More accurate mode components show that the number of false components was reduced. Figure 6 shows the kurtosis spectrum required in the fault diagnosis scheme proposed in this paper to extract the center frequency of the signal. We extracted the signal component on the basis of the central frequency for reconstruction. The results in Figure 6 show that the center frequency was consistent with the analog signal, so it could be used for a follow-up analysis. Figure 7 shows the center frequency of the different rotation components in the decomposition result of the CEITDALN method. The center frequencies of the first three components accounted for the largest proportion. Therefore, the first three components were selected for signal reconstruction, and the envelope spectrum was drawn.   Figure 5 shows that the CEEMDAN method decomposed a total of 12 mode components, and CEITDAN and CEITDALN decomposed 10 rotation components. Fewer model components can save subsequent analysis time for decomposition results. More accurate mode components show that the number of false components was reduced. Figure 6 shows the kurtosis spectrum required in the fault diagnosis scheme proposed in this paper to extract the center frequency of the signal. We extracted the signal component on the basis of the central frequency for reconstruction. The results in Figure 6 show that the center frequency was consistent with the analog signal, so it could be used for a followup analysis. Figure 7 shows the center frequency of the different rotation components in the decomposition result of the CEITDALN method. The center frequencies of the first three components accounted for the largest proportion. Therefore, the first three components were selected for signal reconstruction, and the envelope spectrum was drawn. Figures 8-10 are the envelope spectra of the signals obtained by the CEEMDAN, CEITDAN, and CEITDALN methods, respectively. In the above figure, fn-1 is the fault characteristic frequency, fn-2 is twice the fault characteristic frequency, and the rest can be deduced by analogy. The CEEMDAN and CEITDAN methods had similar denoising capabilities to strong noisy signals. However, the difference between the peak value of the fault characteristic frequency position and the noise peak value of the CEITDAN method was slightly larger than that of the CEEMDAN method. The CEITDALN method had the largest difference and the most obvious fault characteristic frequency. At the same time, the filtering effect of the high-frequency part was obvious compared with those of the CEEMDAN and CEITDAN methods. This did not affect the extraction of the fault characteristic frequency because the proposed method has good application prospects in fault diagnosis.     Figures 8-10 are the envelope spectra of the signals obtained by the CEEMDAN, CEITDAN, and CEITDALN methods, respectively. In the above figure, fn-1 is the fault characteristic frequency, fn-2 is twice the fault characteristic frequency, and the rest can be deduced by analogy. The CEEMDAN and CEITDAN methods had similar denoising capabilities to strong noisy signals. However, the difference between the peak value of the fault characteristic frequency position and the noise peak value of the CEITDAN method was slightly larger than that of the CEEMDAN method. The CEITDALN method had the largest difference and the most obvious fault characteristic frequency. At the same time, the filtering effect of the high-frequency part was obvious compared with those of the CEEMDAN and CEITDAN methods. This did not affect the extraction of the fault characteristic frequency because the proposed method has good application prospects in fault diagnosis.      In order to further compare the advantages and disadvantages of the three methods, OI, ECI, and ISNR were compared with the signal output SNR, as shown in Table 1. Table  1 shows that all the indicators of the proposed method were better than those of the other methods. Therefore, the effectiveness of the proposed method is proven. In order to further compare the advantages and disadvantages of the three methods, OI, ECI, and ISNR were compared with the signal output SNR, as shown in Table 1. Table 1 shows that all the indicators of the proposed method were better than those of the other methods. Therefore, the effectiveness of the proposed method is proven.  Figure 11 illustrates the improvement effect of the proposed COA. The method for optimizing the fitness value and the number of iteration steps was used for comparison. The smaller the optimization fitness value is, the faster the method's optimization speed is, and the easier it is to quickly achieve convergence. Therefore, the GBO-COA method proposed in this paper is better than the traditional COA method, as the proposed method performed well in the simulation results and can be used in engineering practice to decompose strong noisy signals and perform fault diagnosis.   Figure 11 illustrates the improvement effect of the proposed COA. The method for optimizing the fitness value and the number of iteration steps was used for comparison. The smaller the optimization fitness value is, the faster the method's optimization speed is, and the easier it is to quickly achieve convergence. Therefore, the GBO-COA method proposed in this paper is better than the traditional COA method, as the proposed method performed well in the simulation results and can be used in engineering practice to decompose strong noisy signals and perform fault diagnosis.

Case B: Experiment Analysis
In this section, we describe the fault signal acquisition and the test equipment used ( Figure 12). The test platform part consists of the test bench and control system, which is mainly composed of the tested bearings, accompanying test bearings, test spindle, bearing outer ring fixture, drive unit, and loading system. During the test, there are four sets of bearings, divided into two sets of bearings under test and two sets of accompanying bearings. Two sets of loading systems apply loads to the test bearings, respectively. The drive unit provides the power for the whole test bench, and the loading system provides radial force for the test bearings, which causes the spindle to drive the bearings to rotate and

Case B: Experiment Analysis
In this section, we describe the fault signal acquisition and the test equipment used ( Figure 12). The test platform part consists of the test bench and control system, which is mainly composed of the tested bearings, accompanying test bearings, test spindle, bearing outer ring fixture, drive unit, and loading system. During the test, there are four sets of bearings, divided into two sets of bearings under test and two sets of accompanying bearings. Two sets of loading systems apply loads to the test bearings, respectively. The drive unit provides the power for the whole test bench, and the loading system provides radial force for the test bearings, which causes the spindle to drive the bearings to rotate and obtains the vibration signal of the bearings during the rotation process through vibration sensors. The drive unit can provide the bearing speed in the range of 1000 rpm~20,000 rpm, which is continuous. The loading system has a loading range of 0~6 KN and a loading accuracy of ±2% and is continuously adjustable. The operating limit temperature of the test bench is 250 • C [38].   Table 2 shows the test bearing parameters. The sampling frequency was 18,400 Hz, and the data length was 8192 points. The bearing rotation speed was 3000 rpm. On the basis of Table 2, the fault characteristic frequency of the inner loop was calculated as Formula (50), and the result was 407 Hz. Due to computer error, the inner loop fault frequency in the envelope spectrum was 408 Hz.
where N is the number of rolling bodies, r f is the inner ring rotation frequency, d is the rolling body diameter, D is the pitch diameter, and α is the contact angle.  Figure 13 shows that the CEEMDAN method decomposed a total of 13 mode components. The CEITDAN and CEITDALN methods decomposed 10 mode components, respectively. This is consistent with the previous part of the analog signal results, because the three types are relatively stable in decomposing strong noise signals. The CEEMDAN method led to the appearance of false modes due to envelope fitting. The CEITDAN and CEITDALN methods showed better results in terms of signal decomposition. Figure 14 shows the kurtosis spectrum required by CEITDALN in the decomposition process to extract the center frequency and optimize the parameters. Figure 15 shows the percentage of the center frequency contained in the different rotation components of the CEITDALN method. Different from the analog signal, the second rotation component accounted for the largest proportion of the center frequency in the actual engineering signal. However, the results of the first three components were consistent with those of the analog signal. Figures 16-18 show the signal envelope spectrum after decomposition by the CEEMDAN, CEITDAN, and CEITDALN methods, respectively. In the above figure, fn-1 is the fault characteristic frequency, fn-2 is twice the fault characteristic frequency, and the rest can be deduced by analogy. Figures 16 and 17 show that, after the signal processing of the CEEMDAN and CEITDAN methods, the fault characteristic frequency could be observed.  Table 2 shows the test bearing parameters. The sampling frequency was 18,400 Hz, and the data length was 8192 points. The bearing rotation speed was 3000 rpm. On the basis of Table 2, the fault characteristic frequency of the inner loop was calculated as Equation (50), and the result was 407 Hz. Due to computer error, the inner loop fault frequency in the envelope spectrum was 408 Hz.
where N is the number of rolling bodies, f r is the inner ring rotation frequency, d is the rolling body diameter, D is the pitch diameter, and α is the contact angle.  Figure 13 shows that the CEEMDAN method decomposed a total of 13 mode components. The CEITDAN and CEITDALN methods decomposed 10 mode components, respectively. This is consistent with the previous part of the analog signal results, because the three types are relatively stable in decomposing strong noise signals. The CEEMDAN method led to the appearance of false modes due to envelope fitting. The CEITDAN and CEITDALN methods showed better results in terms of signal decomposition. Figure 14 shows the kurtosis spectrum required by CEITDALN in the decomposition process to extract the center frequency and optimize the parameters. Figure 15 shows the percentage of the center frequency contained in the different rotation components of the CEITDALN method. Different from the analog signal, the second rotation component accounted for the largest proportion of the center frequency in the actual engineering signal. However, the results of the first three components were consistent with those of the analog signal. Figures 16-18 show the signal envelope spectrum after decomposition by the CEEMDAN, CEITDAN, and CEITDALN methods, respectively. In the above figure, fn-1 is the fault characteristic frequency, fn-2 is twice the fault characteristic frequency, and the rest can be deduced by analogy. Figures 16 and 17 show that, after the signal processing of the CEEMDAN and CEITDAN methods, the fault characteristic frequency could be observed. However, the difference between the peak values of the characteristic frequency and of the noise was small. High-frequency noise was not completely filtered out, and there was some retention. In the CEITDALN method shown in Figure 18, all the high-frequency noise was filtered out and the peak difference was the largest. The extraction of characteristic frequencies was the most obvious. In engineering applications, it is generally only necessary to observe the fault characteristic frequency of onefold frequency. A more intuitive consideration of the filtering method involves the degree of noise filtering. The results shown by the method in this paper can be used in engineering practice.
Symmetry 2021, 13, x FOR PEER REVIEW 20 of 29 However, the difference between the peak values of the characteristic frequency and of the noise was small. High-frequency noise was not completely filtered out, and there was some retention. In the CEITDALN method shown in Figure 18, all the high-frequency noise was filtered out and the peak difference was the largest. The extraction of characteristic frequencies was the most obvious. In engineering applications, it is generally only necessary to observe the fault characteristic frequency of onefold frequency. A more intuitive consideration of the filtering method involves the degree of noise filtering. The results shown by the method in this paper can be used in engineering practice.       In order to objectively compare the three methods, this section still uses the same indicators as those in the analog signal section. However, due to the actual signal, the input SNR and output SNR could not be accurately calculated. Therefore, only four indicators-OI, ECI, ISNR, and RMSE-were used, as shown in Table 3. The results were consistent with those of analog signals, and the proposed method showed a strong applicability in index comparison. Figure 19 shows the comparison between the improved COA and the traditional COA algorithms in optimizing the fitness value. Although the basic value was improved, the improved method proposed in this paper maintained the characteristic of faster convergence compared with the traditional method.  In order to objectively compare the three methods, this section still uses the same indicators as those in the analog signal section. However, due to the actual signal, the input SNR and output SNR could not be accurately calculated. Therefore, only four indicators-OI, ECI, ISNR, and RMSE-were used, as shown in Table 3. The results were consistent with those of analog signals, and the proposed method showed a strong applicability in index comparison. Figure 19 shows the comparison between the improved COA and the traditional COA algorithms in optimizing the fitness value. Although the basic value was improved, the improved method proposed in this paper maintained the characteristic of faster convergence compared with the traditional method.  In order to objectively compare the three methods, this section still uses the same indicators as those in the analog signal section. However, due to the actual signal, the input SNR and output SNR could not be accurately calculated. Therefore, only four indicators-OI, ECI, ISNR, and RMSE-were used, as shown in Table 3. The results were consistent with those of analog signals, and the proposed method showed a strong applicability in index comparison. Figure 19 shows the comparison between the improved COA and the traditional COA algorithms in optimizing the fitness value. Although the basic value was improved, the improved method proposed in this paper maintained the characteristic of faster convergence compared with the traditional method.  In order to more fully reflect the effectiveness of the proposed method in this paper, we use the same experimental platform for troubleshooting different bearing faults, including outer ring faults. The outer ring fault characteristic frequency was calculated as shown in Equation (51). The experimental bearing outer ring failure characteristic frequency was calculated as 292 Hz.
Bearing outer race (BPFL) = (N/2) f r where N is the number of rolling bodies, f r is the inner ring rotation frequency, d is the rolling body diameter, D is the pitch diameter, and α is the contact angle.  In order to more fully reflect the effectiveness of the proposed method in this paper, we use the same experimental platform for troubleshooting different bearing faults, including outer ring faults. The outer ring fault characteristic frequency was calculated as shown in Equation (51). The experimental bearing outer ring failure characteristic frequency was calculated as 292 Hz.
where N is the number of rolling bodies, r f is the inner ring rotation frequency, d is the rolling body diameter, D is the pitch diameter, and α is the contact angle. According to the previous comparison of the CEITDAN and CEEMDAN methods, it can be seen that the proposed method has a better effect, but in order to have more convincing experimental results, this paper introduces wavelet threshold noise reduction to compare with the proposed method, and considering that the fault diagnosis results are more strict on time requirements, it is also necessary to add a time complexity index to evaluate the different methods. The results are as follows. Figure 20 shows the decomposition results of the three methods for the bearing outer ring fault signal. From the figure, it can be seen that the signal decomposition effect of the proposed method in this paper is better than the other two, which is consistent with the performance of the simulated signal and the experimental signal of the inner ring fault in the previous paper. Figure 21 shows the time domain diagram of the signal after wavelet threshold noise reduction, and a subsequent envelope spectrum transformation of the noise reduction signal will be performed to obtain the fault characteristic frequency. Figure 22 shows the center band energy of the cliff spectrum used in the proposed method. It can be seen that the 2nd mode component contains the greatest center band energy, so this component is selected for the signal envelope spectrum transformation to obtain the fault characteristic frequency. Figures 23-26 show the results of the envelope spectrum transformation for different methods. In the above figure, fn-1 is the fault characteristic frequency, fn-2 is twice the fault characteristic frequency, and the rest can be deduced by analogy. In Figures 23-25, it can be seen that the filtering effect is not satisfactory, the According to the previous comparison of the CEITDAN and CEEMDAN methods, it can be seen that the proposed method has a better effect, but in order to have more convincing experimental results, this paper introduces wavelet threshold noise reduction to compare with the proposed method, and considering that the fault diagnosis results are more strict on time requirements, it is also necessary to add a time complexity index to evaluate the different methods. The results are as follows. Figure 20 shows the decomposition results of the three methods for the bearing outer ring fault signal. From the figure, it can be seen that the signal decomposition effect of the proposed method in this paper is better than the other two, which is consistent with the performance of the simulated signal and the experimental signal of the inner ring fault in the previous paper. Figure 21 shows the time domain diagram of the signal after wavelet threshold noise reduction, and a subsequent envelope spectrum transformation of the noise reduction signal will be performed to obtain the fault characteristic frequency. Figure 22 shows the center band energy of the cliff spectrum used in the proposed method. It can be seen that the 2nd mode component contains the greatest center band energy, so this component is selected for the signal envelope spectrum transformation to obtain the fault characteristic frequency. Figures 23-26 show the results of the envelope spectrum transformation for different methods. In the above figure, fn-1 is the fault characteristic frequency, fn-2 is twice the fault characteristic frequency, and the rest can be deduced by analogy. In Figures 23-25, it can be seen that the filtering effect is not satisfactory, the difference between the fault characteristic frequency and the noise is small, and it is difficult to extract the fault characteristic frequency quickly and accurately. However, the gap between the peak of the characteristic frequency and the peak of the noise is large in Figure 26, and the high-frequency noise is completely filtered out cleanly, showing the same better application effect as the previous paper. Table 4 shows a summary table of the comparison indexes of different methods, and it is obvious from the table that the proposed method in this paper has a better effect than other methods in practical applications.
In summary, the method proposed in this paper maintained a good applicability in engineering practice, solving the problem of excessive signal noise and the influence of nonwhite noise on fault diagnosis. In the traditional decomposition method, the frequency bands of nonwhite noise are mostly concentrated in a specific frequency band. Therefore, it is necessary to analyze the decomposition results. We filtered out mode components that were not white noise. In the proposed method, nonwhite noise could be partially offset and preprocessed in the decomposition process, reducing the proportion of nonwhite noise components in the decomposition results. This improves the signal-to-noise ratio and diagnosis accuracy and has a strong practicability and operability.
cult to extract the fault characteristic frequency quickly and accurately. However, the gap between the peak of the characteristic frequency and the peak of the noise is large in Figure 26, and the high-frequency noise is completely filtered out cleanly, showing the same better application effect as the previous paper. Table 4 shows a summary table of the comparison indexes of different methods, and it is obvious from the table that the proposed method in this paper has a better effect than other methods in practical applications.

Conclusions
This study solves the problem of the traditional CEITDAN method not being able to offset the effects of strong and nonwhite noise. The method proposed in this paper was verified by analog and actual engineering signals. This paper compares the proposed method and the traditional method using both qualitative analysis (the number of feature frequencies extracted from the envelope spectrum and the degree of noise interference included) and quantitative analysis (ECI, OI, RMSE and ISNR), and finds that the proposed method outperforms the traditional method in the fault diagnosis process. The results show that this method has a good performance because it has a certain reference and validity for signal processing and fault diagnosis. The CEITDALN method proposed in this paper changes the traditional decomposition algorithm by adding white noise to groups, adding analog noise in the form of adjustable parameter transformation noise. This can better cancel noise in strong noise signals in the process of signal decomposition. The proportion of nonwhite noise in the signal decomposition result was reduced, reducing the impact of nonwhite noise on the accuracy of the decomposition results. The method performs well in both analog and measured signals and can be further used in signal processing and fault diagnosis in actual projects. In the future, we will conduct in-depth research on the endpoint effects of the ITD method itself to solve the problem of curve distortion and further improve the accuracy and visualization of signal decomposition.