FDS-MOMEDA: optimization-blind deconvolution in finite high-dimensional spaces for extracting pulse signal in rolling bearing fault diagnosis

Rolling bearing fault diagnosis is crucial for ensuring the safe and reliable operation of mechanical equipment. Detecting faults directly from measurement signals is challenging due to severe noise and interference. Blind deconvolution (BD), as a preferred method, effectively recovers periodic pulses from the measured vibration signals of faulty bearings. This study introduces a simulated annealing-based BD approach to enhance the pulse signal components reflecting faults in vibration signals measured on rolling bearings. This method iteratively searches for the optimal coordinates in a high-dimensional orthogonal optimization space, where the optimal coordinates reflect the combination of the inverse filter coefficients. Compared to the generalized spherical optimization space used in the ‘Optimization-Blind Deconvolution’ method in previous works, the proposed finite high-dimensional optimization space helps overcome the problem of inverse filter coefficient convergence, allowing for the design of inverse filters without limit of its shape. To better accommodate the cyclostationarity characteristics of bearing signal measured in reality, the proposed method employs a target vector that allows for uncertainty in pulse occurrence instants, thus overcomes challenges introduced by pseudo-periodic phenomena resulting from bearing slippage. Numerical simulations and experimental results on real bearing vibration signals confirm that the proposed method can design more flexible filters to enhance pulse-like patterns in signals, effectively utilize limited filter resources. Its capacity to tolerate inaccurate fault period estimates, high background noise, and pulse randomness enables it to effectively address vibration measurement signals in real-world scenarios.


Introduction
Rolling bearings are components widely used in numerous rotating machinery, providing reliable and stable support for Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.such system.Their operational condition significantly affects the performance of rotating machinery.Due to prolonged operation in harsh environments such as heavy loads and high temperatures, key components of bearings are prone to damage.These unexpected failures can lead to unplanned downtime or even catastrophic casualties [1,2].Therefore, early identification of failure types and assessment of their severity are of paramount importance.
Various condition monitoring methods have emerged to identify the type of faults, mainly including vibration-based methods [3][4][5], stator current analysis [6], thermal imaging [7], acoustic noise analysis [8], and multi-sensor fusion [9].Among them, vibration analysis is considered the preferred method for predictive maintenance because it typically only requires non-intrusive collection and analysis of vibration signals using vibration sensors, making it reliable for detecting transient faults and for routine maintenance [3,10].
Over the past few decades, vibration-based rotating machinery fault feature extraction techniques have been extensively researched.A range of technologies has been developed, such as time-frequency analysis [11], spectral kurtosis [12][13][14], cyclostationary methods [15], decomposition-based methods [16], morphological filtering [17], stochastic resonance [18], blind source separation [19], and machine learning methods [20].However, regardless of the method used, identifying early-stage faults from measured signals can be challenging, as the pulse patterns caused by faults are often masked by background noise and other interferences.This situation is further exacerbated by the propagation effects of unknown transmission paths.Therefore, effectively extracting weak fault information from measured vibration signals is crucial for accurate fault diagnosis.
When rotating components pass through a faulty area, impacts occur, leading to the generation of pulse components in the measured signal.Consequently, a series of sparse pulses are considered crucial indicators of mechanical faults in rotating machinery [21].When the rotational speed is constant, these sparse pulses exhibit noticeable periodicity or cyclostationarity.In order to extract these pulse components from complex background noise, blind deconvolution (BD) [22] methods are regarded as effective techniques for enhancing pulse features.BD holds unique advantages in bearing fault diagnosis because the transmission process of bearing fault source signals involves convolution mixing with noise, and BD theory can reliably enhance embedded periodic patterns from noise-modulated data by designing filters [23,24].
In the context of mechanical fault diagnosis, the fundamental principle of BD generally involves designing an inverse filter to recover repetitive transient pulses caused by bearing defects from measured signals by maximizing a certain criterion of the filtered signal.Finite impulse response (FIR) filters offer outstanding performance advantages and find wide applications in digital signal processing.In current BD methods, FIR filters are commonly employed as inverse filters.Based on different criteria, many BD methods such as minimum entropy deconvolution (MED) [25][26][27][28], maximum correlation kurtosis deconvolution (MCKD) [29], optimal MED adjustment (OMEDA) [30], multipoint-OMEDA (MOMEDA) [30], and maximum second-order cyclostationarity BD (CYCBD) [31] have been proposed as powerful tools for extracting rolling bearing fault features.
MED method was originally developed for processing seismic signals and is a classical BD approach.It iteratively solves filter coefficients by maximizing the kurtosis of the filtered signal.However, MED tends to enhance only individual pulses in the signal.For rotating machinery, enhancing a series of similar periodic pulses is more desirable.The filter coefficients of OMEDA aim to maximize the D-norm of the filtered signal and can be directly solved using matrix operations.However, similar to kurtosis, the D-norm is easily influenced by signal extrema, leading OMEDA to also tend to enhance single pulses rather than a series of periodic pulses.
In response to these issues, some BD methods aimed at recovering periodic pulses related to faults from the measured signal have been proposed.Among them, MCKD achieves BD by maximizing the correlation kurtosis of the filtered signal, where correlation kurtosis is a time-domain feature parameter used to quantify the intensity of pulses within specific time intervals.MOMEDA is a non-iterative algorithm for recovering periodic pulses by maximizing multi-D-norm of the filtered signal.CYCBD achieves signal recovery by maximizing cyclostationarity indicators of the filtered signal.
Recently [32][33][34], proposed a method for solving the inverse filter of BD using optimization algorithms.The results indicate that the optimization-based BD method outperforms traditional deconvolution methods in rolling bearing fault detection.Particularly, in cases where the prior knowledge of fault frequencies is inaccurate, the optimization-based BD method can provide better pulse enhancement effects.This is particularly significant for BD methods such as MCKD, MOMEDA, and CYCBD, which aim to recover periodic pulses related to faults from the original vibration signal.These methods often rely on theoretical fault characteristic frequencies calculated from bearing geometry as prior knowledge in bearing fault diagnosis.
It has been noted in [35] that the calculation of theoretical fault characteristic frequencies is idealized, as the occurrence of pulses actually exhibits randomness, resulting in signals that are essentially pseudo-periodic.This theoretical assumption considers no slippage in the kinematic relationship of bearings, but in reality, there are always some slippage and assembly issues, leading to deviations from the fault periods by even as much as 1%-2%.Additionally, theoretical fault characteristic frequencies are linearly correlated with bearing rotational speed, and inaccurate speed estimation can also result in inaccuracies in calculated fault characteristic frequencies.In such cases, optimization-based deconvolution methods mitigate the negative impact introduced by inaccurate estimation of fault characteristic frequencies on BD effectiveness [33], thus making them more suitable for practical diagnostic scenarios.
In [32][33][34]'s method, optimization takes place in a proposed generalized spherical optimization space aimed at reducing the search space for optimal coordinates and speeding up the search process.However, this optimization space inevitably introduces the convergence phenomenon of filters.This issue will be discussed in detail later.Essentially, this phenomenon manifests as the initial coefficients of the filter having relatively large values, but subsequent coefficients tend to approach zero.The potential consequence is that filters designed using this method have only a few coefficients that are relatively effective, resulting in a significant waste of filter resources, making it almost impossible to design a longer filter.This conflicts with the original intention of the BD method, as one of the main advantages of BD is to design filters with as much flexibility in shape as possible [36].
To address the aforementioned issues, this study proposes an optimization-deconvolution method for enhancing periodic pulse components in measured vibration signals from faulty bearings.In this method, the optimization is conducted in a finite high-dimensional orthogonal space, avoiding the filter convergence issue introduced by the generalized spherical optimization space in previous works, thus enabling more flexible filter design.This method allows for the effective design of longer filters and better utilization of filter resources.Conversely, the high efficiency in utilizing filter resources implies that this method only requires the design of shorter filters to achieve similar effects as previous methods when processing the same signals.In the proposed method, the optimization process is implemented in a finite high-dimensional space (FDS) based on the simulated annealing (SA) [37] algorithm, and the deconvolution is formulated with the multi-D-norm used in the MOMEDA method as the optimization criterion, thus named FDS-MOMEDA.
To summarize the novelty, this work firstly proposes a feasible approach to solving the inverse filter coefficients combination in the BD using optimization algorithms.Apart from providing a new insight into such a 'optimization-BD' method, the second innovation of this work lies in the consideration of the pulse randomness brought by the cyclostationarity of bearing vibration signals during the implementation of FDS-MOMEDA.This makes it more adaptable to the characteristics of real-world vibration signals compared to some BD methods currently widely used in mechanical vibration signal processing.Unlike MCKD and MOMEDA methods that strictly adhere to fault characteristic frequencies as prior knowledge, FDS-MOMEDA fully recognizes the phenomenon of random variations in pulse occurrence intervals, accepts fault characteristic frequencies only as suggestions rather than fully trusting them, and employs a pulse position indicating vector (target vector) that tolerates the uncertainty of pulse occurrence instants, thereby effectively adapting to the pseudo-periodicity in real vibration signals.
It should be noted that the FDS proposed in this study is an optimization space applied to optimization-BD process, rather than a specific optimization algorithm.Although in this study, the optimization of FDS-MOMEDA is conducted based on SA, the method itself is transferable, allowing for the application of various optimization algorithms.This work utilizes SA as the optimization tool to clarify this optimization space, while the discussion about other various optimization algorithms is beyond the scope.In external files, implementation of FDS-MOMEDA based on SA is provided as an example.

BD
When a part of the bearing experiences mechanisms such as spalling or pitting faults, each rotation brings the rolling elements into contact with the faulty area.Local faults on the rolling elements introduce impacts, thereby exciting high-frequency resonances in the entire structure between the bearing and the sensor.The response to a single resonance will exhibit a pulse-exponential decay pattern and will continue to occur as the bearing operates.The unit pulse decay I(t) occurring at t 0 can be defined by the following (1), where B is the decay parameter determining the decay rate of the pulse, and f n is the resonant frequency of the bearing system.Figure 1 illustrates a typical time-domain waveform of an exponential decay pattern.
From figure 1, it can be observed that the pulse reaches a high energy level at the moment of occurrence and then decays with oscillations.With each rotation of the bearing, a pulse component similar to the one shown in the graph is generated, forming an approximately periodic pulse chains.
The vibration signal sources include a (pseudo) periodic pulse signal chains b generated by defect faults, a random pulse signal d generated by the vibrating bearing cage, discrete harmonic signal h generated by gear meshing, and Gaussian white noise interference signal n.The excitation generated by the vibration signal sources is synthesized into the measured signal of the sensor, denoted as x, after passing through different and complex transmission paths, as shown in figure 2. Let g(i) represent the transmission path influence of signal source i (expressed in the form of FIR filters), the measured signal can be modeled as (2): BD is used to estimate the signal of interest from mixed signals without needing prior knowledge of the specific values of the convolution kernel.The process of separating the desired fault signal from the measured signal can be seen as solving for an inverse filter, which involves filtering out other components from the measured signal using this inverse filter.Selecting an optimal FIR filter f of length L as the design target for the inverse filter, the BD process applied to bearing fault diagnosis can be described as: Here, y represents the optimal estimation of the fault signal b.The deconvolution is the process of maximizing a predefined criterion for the obtained signal.In the context of mechanical fault diagnosis, this criterion often reflects the richness of fault pulse information in the obtained signal.McDonald and Zhao [30] proposed a deconvolution method called multipoint optimal MED adjustment (MOMEDA).MOMEDA employs the multi-D-norm as this criterion.The multi-D-norm is an indicator that describes the degree of prominence of periodic pulse patterns in a signal and is expressed as [30]: In (4), t represents a constant vector indicating the positions of the pulses in signal y.For example, in a signal of length L, if the pulses occur at the n 1 , n 2 ,…,n m sample positions, with a total of m pulses, it can be expressed as: In ( 5), the positions assigned the value of 1 correspond to the n 1 , n 2 ,..,n m sample positions.
The optimal solution for the inverse filter in MOMEDA can be represented as the maximization of multi-D-norm via a set of filters: Defining the matrix for subsequent filter-solving: Here, X 0 is a constant matrix associated with the measured signal: If the autocorrelation function matrix X 0 X T 0 is invertible, the process of solving for the optimal inverse filter is as follows: The recovered fault signal is then given by: In MOMEDA, the target vector t defines the positions and weights of the pulses to be deconvolved, allowing the algorithm to focus specifically on deconvolving periodic pulse sequences.This characteristic is highly suitable for the nature of bearing faults and has been proven effective in addressing issues related to periodic or quasi-periodic pulses.
An important point to note is that, in ( 5), the instants of pulse occurrences are considered as a series of positions in the target vector t with values of 1, reflecting the early model of bearing fault vibration signals [38], which assumed the occurrence of impulses to be strictly periodic and deterministic.However, the accuracy of this model compared to real bearing signals has been then questioned, as discussed in detail in [39,40].Later [41], recognized the randomness of impulse occurrences and pointed out the presence of some slip in the rolling elements, resulting in uncertainty in the arrival times of impulses.Randall et al [42] modeled this phenomenon, specifically modeling the arrival time T i of the ith impulse as: where T d represents the (average) fault period, and δ i denotes the jitter, thus forming the cyclostationary (CS) model.Although this uncertainty ranges from 1% to 2%, it fundamentally alters the statistical structure of the vibration signal.Subsequently [43], proposed a more accurate model where the uncertainty (variance of jitter) increases with the number of cycles.Specifically, it is formulated as: Thus, a pseudo cyclostationary (PCS) model was constructed.In both of these two models in ( 11) and ( 12), impulses are no longer considered strictly periodic but pseudo-periodic, thus providing a more accurate description of actual bearing fault vibration signals.
In this scenario, it is necessary to refine the target vector t in equation ( 5) to more accurately reflect the faults.Here, for simplification, the case of fewer pulses in the pulse chains is considered, and the uncertainty increasing is assumed to be smaller with increasing number of cycles, i.e. a simplified perfect CS model is used to characterize impulse chains.In this case, the jitter δ i describes the deviation of the actual time-point T i of the ith impulse from the ideal time-point given in equation ( 5).This deviation's objective existence can be described by constructing normal distributions centered at each ideal timepoints.For instance, if T i follows a normal distribution with a standard deviation of σ = 0.02T d relative to the ideal timepoint, this model will suspect with 68.27% confidence level that the actual pulse time-point T i will fall within the range of (i−0.02)Td to (i + 0.02)T d .Using this method, the constructed target vector t will allow for pulse jitter under cyclostationary conditions, better fitting the characteristics of actual measured bearing vibration signals.A typical target vector is illustrated in figure 3, where in this example, with the average fault period T d is 1000 sample points.
The target vector characterizes the pseudo-periodic properties of bearing pulse CS models, allowing for deviations of individual pulse occurrence times relative to ideal times, better adapting to actual bearing vibration pulse signals.

Optimization-BD
The methods for solving inverse filters in deconvolution algorithms vary, but they essentially involve optimization problems: the variables to be optimized are the coefficients of the inverse filter, and the optimization objective is a sparse indicator of the filtered signal.In the context of bearing fault diagnosis, the optimization objective should be closely related to the characteristics of the fault pulse signal of interest, such as the multi-D-norm mentioned in the previous section.
In the context of constructing FIR filters, the order of the filter may need to be sufficiently long to achieve satisfactory results.Consequently, this directly expands the search space of the optimization process, as it always aims at designing the filter coefficients.The dimensionality of the optimizationdeconvolution problem increases with the length of the filter.Maximizing the objective function during the deconvolution process constitutes an unconstrained optimization process, with the search space covering the entire L-dimensional Euclidean space (where L is the length of the filter).
Cheng et al [32][33][34] argue that, in order to shorten the search time and accelerate the search speed, it is necessary to narrow down the search space.Cheng et al [32][33][34] proposed an optimization-BD method that employs an L−1 dimensional hyperspherical search space.This method uses a generalized spherical coordinate transformation to project the coefficients of the inverse filter onto a point on the L−1 dimensional hypersphere with a radius of r, thereby narrowing down the search space.The vector of coefficients for the inverse filter can be represented as: where L is the length of the inverse filter.θ s (s = 1, 2,.., L−1) are the angle parameters of the generalized spherical coordinate transformation, with values ranging from [−π, π].This method transforms the process of solving inverse filter coefficients into an optimization process of finding the vector −1 .In this case, the search space of the deconvolution algorithm is: Based on this, the optimization-deconvolution problem becomes finding an optimal point within the aforementioned search space (14), where the coordinates of this point, projected according to (13), represent the final coefficients of the inverse filter.
However, if it is assumed that the probability density of the values of each angle in (13) follows a uniform distribution, it is possible to obtain the probability distribution of the nth filter coefficient.Taking the first 6 filter coefficients as an example, estimating their probability distribution using the Monte Carlo method yields the results in figure 4.
It can be observed that as the length of the filter increases, the absolute value of the nth filter coefficient tends towards 0. Reviewing the constructed generalized spherical coordinate system (13), it can be noted that the nth coefficient is obtained by multiplying a constant radius r by n trigonometric functions.Since the values of these trigonometric functions always lie within the range [−1, 1], with increasing order, the absolute value of the nth coefficient of the inverse filter tends to converge rapidly towards 0. This potential issue arises because, regardless of how long the filter is designed, a majority of the filter coefficients may converge to values close to zero.Furthermore, the coefficients of the filter designed in the [32][33][34]'s method are eventually bounded because each filter coefficient is obtained by multiplying a fixed radius r (the constant generalized spherical radius) by multiple trigonometric functions (each trigonometric function value lies between −1 and 1).Consequently, each filter coefficient is effectively constrained to the interval [−r, r].

FDS as optimization Space
To address the above-mentioned issues, this study proposes a new optimization approach that utilizes a finite-dimensional orthogonal space as the optimization space.For an inverse filter of length L, a L-dimensional orthogonal space is defined, with lower and upper bounds specified for mth dimension, denoted as [b m , u m ].This construction yields an L-dimensional space: When there are precise boundaries for the upper and lower bounds in each dimension, this space is called a bounded finitedimensional orthogonal space (BFDS); otherwise, it is termed a finite-dimensional orthogonal space (FDS).Based on this space, the search process for the coefficients of the inverse filter becomes finding a point within this finite-dimensional space, where the coordinates of this point project onto each dimension and represent the coefficients of the inverse filter.The dimensionality of the search space is determined by the requirement of the length of the inverse filter, as illustrated in figure 5.
In [32][33][34], the filter coefficients are constrained to the range [−r, r].In contrast, BFDS allows for the adjustment of upper and lower bounds to avoid restricting the coefficients to too narrow of a range.Although there is no universal criterion for selecting the size of the upper and lower bounds, this study suggests choosing a narrow range centered around zero to enhance efficiency.Only when it is not possible to find a combination of filter coefficients that meet the task requirements through experimentation should the consideration of enlarging the upper and lower bounds to expand the search range be made.
Although specifying the lower and upper bounds [b m , u m ] for each dimension can limit the filter coefficients, having bounds (B) is not a mandatory requirement.In such cases, [b m , u m ] = [−∞, +∞], and the search space is referred to as FDS.In this scenario, only the specification of initial values for the filter (e.g.setting all filter coefficients to zero or having each coefficient follow a centered normal distribution around zero) is required.The algorithm still allows for the search for the optimal filter in the vicinity of these initial values.
Similar to the generalized spherical optimization space in [32][33][34], FDS is eventually an optimization space that reflects a new optimization-BD method.Once the optimization space is determined, there are always various options for the specific optimization methods to be used.This work utilizes SA to clarify FDS, while the discussion about other various optimization algorithms will not be carried out.
Taking SA as optimization tool, the implementation principle of FDS-MOMEDA is explained below.
SA [37] is a method designed based on the principles of solid-state annealing.Solid-state annealing is a metallurgical heat treatment process involving heating metal to a certain temperature and holding it for a period, followed by cooling at an appropriate rate.During the heating stage, the particles within the solid undergo increased thermal motion due to the elevated temperature, leading to an increase in internal energy and disorder in molecular arrangement.As the cooling process proceeds slowly, the particles gradually become ordered, reaching equilibrium states at each temperature, ultimately returning to the ground state at room temperature, with their internal energy minimized.At low temperatures, the internal energy of the molecules within the solid is low, causing only slight vibrations about their original positions as the molecules return to an ordered arrangement.
SA is a probabilistic optimization algorithm that simulates the behavior described above.In the aforementioned context, 'temperature' reflects the degree of intensity of perturbation in the current state.Initially, the system has a high temperature, and there is a greater likelihood of the current solution position jumping to a new position.As the temperature decreases, the intensity of perturbation decreases, and there is a greater  likelihood of the current solution position staying in place.At each temperature, if the new solution is better than the current solution, it is accepted; otherwise, based on the Metropolis criterion, the algorithm decides whether to accept the new solution with a certain probability.This acceptance of non-better solutions reduces the possibility of SA getting trapped in local optima and makes it more likely to find the global optimum.
Based on this, the basic process of implementing FDS-MOMEDA using SA can be described as follows: • Select a (B)FDS as the optimization space for SA.
• Randomly find a position within the (B)FDS as the initial position and specify the initial temperature (reflecting the degree of intensity of perturbation from the current position to the new position).• Iterate cooling.With each cooling step, the intensity of perturbation from the current solution position to the new position decreases.For each perturbed new position, decide whether to accept this new position or keep the previous position.If the new solution is better than the current solution, accept the new solution; otherwise, based on the Metropolis criterion, decide whether to accept the new solution with a certain probability.The criterion for determining whether it is better is whether the multi-D-norm of the filtered signal obtained by the newly designed filter is higher than that of the current solution.• If the number of iterations reaches the specified number (or the temperature has decreased to the specified minimum temperature), terminate the iteration, with this position as the optimal solution found by SA in the (B)FDS.• The projection of the optimal solution on each dimension represents each coefficient of the optimal filter.
Although visualizing high-dimensional spaces is quite challenging, describing a three-dimensional BFDS as an example can help.In this case, the goal is to design a filter of length 3. The optimal position found by SA in this 3D BFDS corresponds to a 3D coordinate, and the projection of each dimension of the coordinate represents the value of one of the three filter coefficients.For example, if the optimal position is (x, y, z), then the combination of FIR filter coefficients obtained by solving is [x, y, z], as shown in the figure 6.

Numerical simulation
To validate whether this method effectively finds the occurrence of pulses in the signal, simulated bearing fault signals were first constructed for experimentation.Here, the simulated signals were constructed according to [34]'s recommendations, comprising 2500 sampling points.Among them, (a) represents pulse chains triggered by bearing fault, (b) represents random pulses, (c) represents discrete harmonic signal, (d) represents Gaussian noise signal, and (e) represents the superimposed signal, as shown in figure 7.
The methods proposed by [32][33][34], BFDS-MOMEDA, and FDS-MOMEDA were employed to design filters of lengths 25, 50, 100, and 200.In the case of BFDS, the lower and upper bounds for each dimension were set to −1 and 1, meaning that each coefficient of the designed filter falls within the interval [−1, 1].Similarly, when using the method proposed by [32][33][34], the radius of the hypersphere was also set to 1, implying that each coefficient of the designed filter also falls within the interval [−1, 1].
The inverse filters in figure 8 were obtained using the three methods.
The previous speculation was that in the method proposed by [32][33][34], as the index of the FIR filter coefficients increases, the coefficients tend to take smaller values.It can be observed that as the index of the filter coefficient increases, the coefficient corresponding to that index tends towards 0. This phenomenon is consistent with this speculation.For filter design, this characteristic is always a limitation because the full potential of the filter cannot be realized within the given filter length.Furthermore [36], has pointed out that one of the significant advantages of the BD method compared to other filter design methods is the flexibility in designing FIR filters, which is not constrained by bandwidth, center frequency, or even fixed shapes.However, this convergence phenomenon weakens this advantage of the BD method.
In contrast, BFDS and FDS methods have designed more flexible-shaped filters, fully utilizing the given filter length.In the BFDS-MOMEDA method, each filter coefficient is constrained within the given upper and lower bounds, but there is no convergence issue.Each position of the filter is fully utilized within this range.FDS-MOMEDA has further surpassed the limitations of upper and lower bounds, designing even more flexible filter shapes.
This study considers the method proposed by [32][33][34] as an very efficient and novel BD method, but unfortunately, the issue of filter coefficient convergence is objectively present.In comparison, the (B)FDS-MOMEDA method can overcome this problem.Especially when designing long filters, the proposed method performs better as it can more fully utilize filter resources.
The results in figure 9 shows the pulse recovery signals obtained by the three methods at different filter lengths.From figure 9, it can be observed that in the method proposed by [32][33][34], the characteristics of pulses in all enhanced signals are relatively indistinct.This is because the filters designed using the [32][33][34]'s method are always constrained by shape and cannot achieve satisfactory performance.For BFDS-MOMEDA, due to the relatively small chosen upper and lower bounds (to facilitate comparison, they are consistent with the bounds in the [32][33][34]), the pulse components in the signal are more significant but still not satisfactory.In contrast, regardless of whether shorter or longer filters are used as the design target, the filters obtained using the FDS-MOMEDA method achieve quite impressive results, with the pulse components in the signal being quite noticeable.This is because FDS-MOMEDA designs filters without being constrained by shape and can more flexibly handle complex signals.This indicates that FDS-MOMEDA can be preferred as a solution for recovering fault pulse signals.
It is important to note that longer filters are not always necessarily better, as longer filters can also result in poor handling of the signal front end.This phenomenon can be observed from the results obtained by designing a filter of length 200 using the FDS-MOMEDA method (where the front end exhibits a less pronounced pulse feature).This is because using FIR filters, when processing signals, due to the convolution operation, the original signal front end contains fewer data points than the filter length, resulting in incomplete convolution.
Furthermore, there is a noticeable difference in the amplitude of the signals obtained by the three methods.However, this phenomenon should not be used as a performance metric for evaluating BD methods.When applying BD to process mechanical fault signals, we are more concerned about whether the signals generated by BD exhibit good pulse patterns rather than the overall signal amplitude.As long as good pulse patterns can be observed, it can serve as evidence of suspected mechanical faults.

Case 1: experiment on Case Western Reserve University (CWRU) dataset
The vibration data from the CWRU bearing data center [44] has become one of the standard references for testing bearing fault detection methods [45].The experiment utilized SKF 6205-2RS JEM deep groove ball bearings, with singlepoint faults introduced using electrical discharge machining.Table 1 documents the geometric parameters of the bearings.
It should be noted that the CWRU dataset is relatively characterized by 'clear and simple' bearing fault vibration signal data, with the bearings primarily subjected to radial loads.This implies that envelope analysis can diagnose nominal bearing faults in most inner and outer race fault scenarios of the CWRU dataset.The reason for selecting it for investigation in the first Case is that this dataset has been extensively utilized and considered one of the benchmark datasets in previously reported   works.In this section, FDS-MOMEDA will primarily be compared with the MOMEDA in generalized spherical optimization space to validate its superiority in overcoming filter convergence issues.
The experiment involved vibration data collected from two types of faulty bearings in the CWRU dataset, sampled at a rate of 12 kHz.The faults occurred on the outer race (recorded as X234_DE_time in dataset) and inner race (X209_DE_time), with speeds of 1796 and 1797 rpm, respectively.Additionally, reference signal was obtained from bearings without faults at a speed of 1796 rpm (X097_DE_time).
Based on the provided parameters in table 1, the theoretical fault characteristic frequencies of the bearings can be calculated.Several studies have reported methods for calculating theoretical fault frequencies based on bearing geometric parameters and rotational speeds [4,35,46].In summary, for different types of bearing faults, geometric formulae are used to compute the fault characteristic frequencies of the bearings, as follows.
The frequency of the rolling element passing through the outer ring fault point, that is, the characteristic frequency of the outer ring fault f o (also commonly named as BPFO, ball pass frequency outer race): where f r is the rotation frequency of bearing.The frequency of the rolling element passing through the inner ring fault point, that is, the characteristic frequency of the inner ring fault f i (also commonly named as BPFI, ball pass frequency inner race): In ( 16) and ( 17), the diameter of the rolling element is d, the pitch diameter of the bearing is D, the number of rolling elements is Z, and the contact angle is α.
From ( 16) and ( 17), the theoretical fault characteristic frequencies BPFO and BPFI for the outer race and inner race can be calculated respectively.Nevertheless, CWRU dataset is a dataset where fault features are relatively evident, thus fault characteristic frequencies can also be identified from the envelope spectrum.
The time-domain waveforms for faults in the outer race and inner race are depicted in figure 10.Here, it is worth noting that the prevailing sparse pulses in figure 10 do not reflect the pulse chains corresponding to the outer and inner race fault because the characteristic frequency of the outer race fault corresponds to a period of approximately 112 samples, and for the inner race fault, it corresponds to a period of approximately 74 samples.However, the significant pulse frequency in figure 10 clearly deviates from these fault frequencies.In other words, the pattern of fault pulse chains is masked within the signal.
Filters of lengths 100, 300, and 500 were designed separately, resulting in the outcomes in figure 11.
It can be observed that the method proposed by [32][33][34] consistently fails to effectively identify the pulse chains  reflecting the faults.In contrast, with an increase in filter length, the FDS-MOMEDA method produces signals with increasingly prominent periodic pulse components.For the outer race fault, when the filter length is 300 and 500, the obtained signals already exhibit observable fault pulse chains.Similarly, for the inner race fault, when the filter length is 500, the obtained signals also clearly reveal the fault pulse chains.
Envelope spectra were generated for the signals obtained by processing the outer race and inner race fault signals using 500-length filters designed by FDS-MOMEDA, yielding the results in figure 12.
From the envelope spectra, it is evident that there are prominent peaks at the BPFO and BPFI frequencies and their harmonics.This indicates that the portion of the signal reflecting the pulse components prevails, while other information is effectively suppressed.It can be concluded that FDS-MOMEDA effectively extracts the components in the original signal reflecting the fault pulses.
However, similar to other BD methods that rely on prior knowledge of fault periods, the filters constructed by FDS-MOMEDA, while effectively enhancing fault-related pulse modes, also result in some loss of additional information for diagnosis.This is manifested in the envelope spectra where spectral lines directly related to nominal faults are prominently enhanced, but additional phenomena useful for diagnosing looseness faults, such as modulation sidebands, are suppressed.Such phenomenon may not occur in some non-faultcycle prior BD methods (e.g.MED).However, regardless, this disadvantage is somewhat of a double-edged sword, as it indeed provides a fairly satisfactory reflection of nominal faults.
The experimental results confirm the effectiveness of the FDS-MOMEDA method and its superiority over the method proposed by [32][33][34].In scenarios where longer filters are required to extract fault pulse components in BD, FDS-MOMEDA can more effectively recover the signal components reflecting the fault pulses from the original signal, thus providing a solid recommendation for bearing fault diagnosis.

Case 2: experiment on Safran surveillance dataset
While FDS-MOMEDA showed promising performance on the CWRU dataset in Case 1, it remains uncertain whether it is applicable to the complex and dynamic operating conditions of bearing detection scenarios in the real world.
Case 2 investigates the implementation of FDS-MOMEDA to diagnose bearing faults in an aircraft engine, using vibration data provided by Safran Group.The data originates from Exercise 2 of the 'Safran contest' [47,48], with details available in [48].Case 2 focuses particularly on signals collected by the Acc2 accelerometer in this dataset, used for diagnosing outer race faults in bearings located at position L5.Compared to the CWRU dataset in Case 1, the complexity of this dataset is significantly higher.Antoni et al [48] has noted that due to the presence of numerous discrete components related to shaft rotation, gear meshing, their numerous harmonics, and their interactions in the form of sidebands, the characteristics of bearing faults are completely obscured, with the RMS value of the masking components being at least an order of magnitude higher than that of the bearing fault signals.
Another major difference between the Safran dataset and the CWRU dataset is that the bearings in the CWRU dataset are driven by motors mounted on test rigs, lacking mechanisms such as gears to convert torque into the axial loads borne by the bearings, which almost exclusively bear pure radial loads.In contrast, the Safran dataset is collected from a real aircraft engine gearbox, featuring a complex transmission structure where the bearings must also withstand axial loads.In such scenarios, when certain rolling elements encounter restrictions in the clearance of the cages, they are forced to skid, resulting in random slipping of the entire cage.Moreover, as bearing damage increases, more cage slippage occurs.Thus, there is ample reason to suspect that this slippage may lead to pseudoperiodic phenomena in pulses, potentially posing a significant obstacle in the diagnostic process.
Some BD methods, such as the MCKD, strictly require prior knowledge of fault periods, which may become unsuitable for such a scenario.In contrast, in FDS-MOMEDA, the target vector used to indicate pulse positions takes this into account and adopts a target vector construction method as described in Chapter 2, allowing for the presence of jittered pulses deviating from ideal time-points.In this section, FDS-MOMEDA will be compared with widely-used MCKD, to highlight its advantages in dealing with cyclostationary pseudo-periodic signals.
In the Safran dataset, a sub-signal from 0.6 s to 0.7 s collected by the Acc2 accelerometer is selected for analysis.At this time, the bearing is operating at a constant speed, with the shaft rotational frequency (f L4 ) at approximately 181 Hz, and the BPFO of the bearing with an outer race fault at position L5 is 7.8855 * f L4 .
Figure 13 illustrates the sub-signal and the signals obtained after filtering with the MCKD and FDS-MOMEDA methods, along with their spectra.From the figure, it is evident that the original signal is complex, with pulse features completely submerged in other signal components, making it difficult to effectively extract components highly related to faults from the spectrum.At this point, the MCKD method, which should have served to recover pulses, also fails to handle this signal well due to its strict reliance on pulse cycles as prior knowledge.In the signal extracted by MCKD, the pulse features are not significant, with many pulse intervals still masked by other components.Additionally, there are significant in-balances in the intensities of pulses in different time intervals, as marked in figure 13.
In contrast, FDS-MOMEDA effectively enhances the pulse components in the signal, with pulse-like components clearly visible in the filtered time-domain waveform.From the spectrum, the characteristic frequencies of the outer race fault and many harmonics are prominent, providing ample evidence for locating the outer race fault.
The target vector used in FDS-MOMEDA to indicate the occurrence times of faults takes into account the pseudoperiodic characteristics brought by signal cyclostationarity, accepting fault characteristic frequencies only as suggestions rather than blindly trusting them.The results in Case 2 demonstrate that FDS-MOMEDA, as a method respecting the uncertainty of pulse occurrence times in fault-bearing vibration signals, can more effectively handle real-world bearing vibration fault signals than the widely-used MCKD method.

Inaccurate fault feature frequency estimation
When using BD method to extract periodic fault impulses, the prior knowledge of fault frequencies is often required, and FDS-MOMEDA is not an exception.In the experiments of the fourth section, the theoretical fault characteristic frequencies  of the bearings were used as such prior knowledge.The calculation of theoretical fault characteristic frequencies assumes no slip in the kinematic relationship of the bearing, but in reality, there are always some slip, assembly issues, and deviations in speed estimation.Consequently, each occurrence of a fault pulse carries uncertainty, and the average value of the pulse period also varies.In such cases, it is necessary to evaluate the tolerance of the proposed method to inaccuracies in fault (average) frequency estimation.
Taking the signal in Case 2 as an example, if the theoretical fault characteristic period is assumed to be an unbiased estimate of the true fault period, deviations of 98%, 99%, 99.5%, 100.5%, 101%, and 102% of the period corresponding to BPFO can be used as prior knowledge with bias.Subsequently, FDS-MOMEDA was employed to design filters, yielding the results shown in figure 14.
From figure 14, it can be observed that when using period estimates with slight deviations, the spectral lines highly correlated with fault components and harmonic can still be observed in the spectra.When using the unbiased period estimate of 98% as prior knowledge, the intensity of the third harmonic component displayed in the spectrum is very high, causing other harmonic components to become less prominent.However, these harmonic components are eventually still distinguishable compared to the surrounding regions.This phenomenon confirms that FDS-MOMEDA allows for frequency estimates with slight deviations and can effectively diagnose nominal faults under such circumstances.
Such performance is crucial in practical fault diagnosis, as pointed out by [33], where deviations in fault frequency are inevitable in real-world applications, and theoretical fault characteristic frequencies are unlikely to fully reflect the true fault frequencies.FDS-MOMEDA's tolerance regarding to that allows it to recover pulse chains using approximate fault frequencies, and be more competitive in dealing with the cyclostationarity of measured signals in the real world.The results obtained by FDS-MOMEDA using BPFO and BPFI as prior fault frequency knowledge to process signal from a healthy bearing.

Misdiagnosis of BD
Using the BD method to recover the pulse components from bearing vibration signals raises a potential concern: the possibility of generating pseudo-pulse signals in the output signal when processing signals from healthy bearings.This occurrence could lead to misdiagnosis of faults, mistakenly identifying healthy bearings as faulty.Therefore, it is necessary to verify whether FDS-MOMEDA produces such pseudo-pulses when processing signals from healthy bearings.
Although the signals in Case 2 appear to be more representative of real-world scenarios, it is regrettable that the dataset does not provide signals collected from healthy bearings as a baseline.The CWRU dataset in Case 1 provides a vibration signal (X097_DE_time) collected from a healthy bearing.By estimating the fault frequencies as the BPFO and BPFI, and employing the FDS-MOMEDA method to process these signals, the results are shown in figure 15.From figure 15, it can be observed that when FDS-MOMEDA processes signals collected from healthy bearings, there are no significant periodic pulse components in the resulting signal.This phenomenon indicates that FDS-MOMEDA only produces noticeable pulse patterns in the resulting signal when fault-related pulse components are indeed present in the vibration signal.This further confirms an important prerequisite for applying FDS-MOMEDA to fault diagnosis in bearings: FDS-MOMEDA does not generate pseudo-pulses and cause misdiagnosis of faults when faults have not occurred.

Reason for using shorter signal as the subject
It is worth noting that in this study, whether it is numerical simulation signals or real signals in the two cases, they are all studied with relatively short durations (or fewer sampling points) as the objects of FDS-MOMEDA.For example, the numerical simulation signals in section 2 consist of only 2500 sampling points, the two faulty signals in Case 1 are only 1 s long (12 000 sampling points), and the faulty signals in Case 2 are only 0.1 s long (5000 sampling points).
There are two reasons for this approach.Firstly, from the perspective of the pseudo-periodicity of cyclostationary signals, there is less risk of bearing rotational frequency fluctuations within a shorter duration.Although the FDS-MOMEDA method does not entirely rely on theoretical fault characteristic frequencies and allows for the phenomenon of pulse occurrence time jitter, it still accepts theoretical fault characteristic frequencies as important information, which is indicated in the target vector as the ideal occurrence time of pulses.However, the calculation of theoretical fault characteristic frequencies is strictly linearly correlated with the bearing rotational frequency, as shown in equations ( 16) and (17).If there are significant fluctuations in the bearing rotational frequency, it will cause the signal to no longer exhibit pseudoperiodicity and become completely non-stationary, rendering the FDS-MOMEDA method inapplicable.This potential concern applies not only to the validity of FDS-MOMEDA but also to other BD methods that rely on fault cycles as prior knowledge, including MOMEDA and MCKD.
Another reason is related to the target vector established in this study, which models the pulse chains based on the arrival time of the ith pulse, T i = iT d + δ i [42] (where T d is the average fault period, and δ i represents jitter).However [43], has pointed out that a more accurate model should consider the uncertainty (variance of jitter) increasing with the number of pulse cycles.Here, when the number of pulse cycles is not very large, it is approximated that the increase in the variance of jitter is relatively small.Therefore, this simplification poses a potential limitation when dealing with a larger number of pulse cycles.Undoubtedly, when dealing with longer signals, the increase in this uncertainty must be further considered, thereby demanding a more accurate method for constructing the target vector.This could be one of the main issues to address in the future.

Conclusion
This study proposes an optimization-BD method for fault diagnosis using rolling bearing vibration signals, termed FDS-MOMEDA method, aiming to recover periodic pulse components reflecting faults in the signal.This method employs SA algorithm in a high-dimensional space to search for FIR filter coefficients, maximizing the multi-D-norm metric of the signal to highlight the pulse chains in the signal.FDS-MOMEDA has addressed the issue of filter coefficient convergence in previous optimization-BD methods and resolved it by defining a new optimization space.Numerical simulations demonstrate that the FDS-MOMEDA method can design FIR filters with more flexible shapes, thus achieving better pulse enhancement effects.Experimental validation on signals collected from real faulty bearings confirms that FDS-MOMEDA can effectively enhance the fault pulse modes of both outer and inner race faults hidden in complex noise and suppress other information.FDS-MOMEDA respect the fact that the impulses chains hidden in measured signal is eventually cyclostationary, and the occurrence instants of pulses always exhibit randomness.It tolerates slight deviations in fault frequency estimation, and its capability in this regard surpasses that of the MCKD method.Furthermore, it has been discussed that FDS-MOMEDA does not produce pseudo-pulses when processing signals from healthy bearings, thereby avoiding misdiagnosis.As a novel optimization-BD method, FDS-MOMEDA demonstrates more reliable performance than previous similar methods and can be effectively applied to fault diagnosis tasks in rolling bearings.
The characteristics of FDS-MOMEDA suggest its suitability for dealing with phenomena encountered in real-world measured vibration signals, including the existence of pseudoperiodic pulses and complex noise.It performs well in extracting fault pulse features, indicating its potential to play a significant role in future equipment fault diagnosis and condition monitoring based on vibration measurements.

Figure 1 .
Figure 1.Typical time-domain waveform of an exponential decay pattern.

Figure 2 .
Figure 2. Modeling of measured vibration signal on bearing.

Figure 4 .
Figure 4. Monte Carlo simulation result of the first 6 filter coefficients.

Figure 6 .
Figure 6.A schematic diagram of filter optimization using SA on a 3D BFDS.

Figure 7 .
Figure 7.The simulated vibration signal of a faulty bearing and its components.

Figure 8 .
Figure 8.The FIR filters of different lengths designed by the three methods.

Figure 9 .
Figure 9.The pulse recovery signals obtained by the three methods using FIR filters of different lengths.

Figure 10 .
Figure 10.The two types of bearing fault signals collected from the CWRU dataset.

Figure 11 .
Figure 11.The pulse recovery signals obtained when designing filters of different lengths using the two methods, (a) outer race fault, (b) inner race fault.

Figure 12 .
Figure 12.Envelope spectra of recovered signals, employing FDS-MOMEDA to design filters with length of 500.

Figure 13 .
Figure 13.Sub-signal and the signals obtained after filtering with the MCKD and FDS-MOMEDA methods and their spectra.

Figure 14 .
Figure 14.Pulse recovery signals' spectra obtained by FDS-MOMEDA under different degrees of fault period biases.

Figure 15 .
Figure 15.The results obtained by FDS-MOMEDA using BPFO and BPFI as prior fault frequency knowledge to process signal from a healthy bearing.

Table 1 .
Geometric parameters of the bearings.