© M. Yazdanianasr et al., Published by EDP Sciences, 2025

1 Introduction

Rolling element bearings play a pivotal role in the functioning of rotating equipment. Failure to timely detect any potential damage in bearings can significantly hinder the accurate planning of maintenance interventions. In the worst-case scenario, this can result in catastrophic failures, substantial production losses, and an escalation in overall costs. Therefore, it is crucial to promptly identify any impairment to ensure effective maintenance planning and mitigate the risks. Hence, vibration-based condition monitoring serves as a highly effective approach for diagnosing bearing faults, complementing maintenance strategies aiming at identifying potential machine malfunctions. Nevertheless, signal processing within this context encounters formidable challenges when applied to complex machinery, where multiple concurrent factors contribute to the characterization of the machine's condition. A significant consideration lies in the fact that the signatures emitting from bearing fault signatures are typically obscured amidst a diverse range of noise sources and more dominant components such as gears. Moreover, within the domain of complex machinery, the presence of multiple rolling bearings operating in close proximity is common. Consequently, it is imperative to meticulously explore and identify all potential faults, as some may possess weaker signatures yet still pose a significant risk comparable to faults with more prominent signatures. Examples of these intricate machinery systems are abundant in industrial domains, including wind turbine drive trains, power generation facilities, mining sectors, agricultural enterprises, etc. Hence, there is a need to employ a robust signal processing methodology for effective fault diagnosis within complex systems.

In the realm of bearing diagnostics, the traditional practice of utilizing envelope analysis has long been employed. The fundamental principle of envelope analysis involves the demodulation of signals by subjecting them to a band-pass filter, which restricts the frequency range to a specific interval. Nevertheless, an ongoing question among analysts pertains to the optimal range for applying a band-pass filter to effectively observe the bearing fault signatures. Typically, the vibration signals are filtered around their resonance frequencies, which are stimulated by damaged impulses. These resonance frequencies essentially act as carrier frequencies for the modulation frequencies. The aim is to achieve an optimal filter band that exhibits a high Signal-to-Noise ratio (SNR), resulting in a filtered Squared Envelope Spectrum (SES). This spectral representation effectively amplifies the fault harmonics [1,2]. Recently, researchers have been striving to develop (semi-) automated approaches for the selection of an optimal frequency band, thereby surpassing the reliance on purely subjective engineering knowledge in the context of filtering applications. One of the well-known methods of automated band selection tools is the Fast Kurtogram [3] where the band selection is done based on the maximum kurtosis level. Also, an enhanced version of the Fast Kurtogram (FK), known as the Protrugram, has been introduced [4]. The Protrugram improves upon the FK methodology by extracting the maximum kurtosis from the envelope spectrum rather than the raw signal. The Sparsogram approach [5] operates on evaluating the sparsity levels across various frequency bands using wavelet-packet analysis. On the other hand, the Infogram technique [6] utilizes Negentropy as a feature to identify the impulsive frequency bands within the signal, facilitating the subsequent demodulation. Additionally, the Autogram [7], relies on the concept of maximum kurtosis. However, the Autogram distinguishes itself from the Fast Kurtogram (FK) by computing the maximum kurtosis using the unbiased autocorrelation of the derived squared envelope. It utilized the (undecimated) Maximal Overlap Discrete Wavelet Packet Transform (MODWPT) to split the signal into frequency bands. Additionally it was concluded that multi-band integration is needed by combining filtered Squared Envelope Spectra (SES) with the highest kurtosis at each level, resulting in the Combined Squared Envelope Spectrum (CSES).

Conversely, over the past two decades, the Cyclic Spectral Correlation (CSC) and the Cyclic Spectral Coherence (CSCoh) have emerged as alternatives to SES-based methods [8–10]. A noteworthy advantage of the CSC and CSCoh methods lies in their capability to reveal hidden periodicities associated with second-order cyclostationarity, such as bearing signals that may be obscured by stronger signals and noise. These methods represent the signals in bi-variable maps within the frequency-frequency domain, enabling integration along the spectral axis to obtain either the Enhanced Envelope Spectrum (EES) or the Improved Envelope Spectrum (IES). The bi-variable maps generated by CSC and CSCoh require expert scrutiny to manually select the optimal band for spectral axis integration, leading to the generation of the Improved Envelope Spectrum (IES). To address this limitation, a recent solution known as the IESFOgram using Targeted and Blind Features [11,12] has been proposed as an automated band selection tool for the bi-variable maps produced by the CSC and the CSCoh. The efficacy of this method was evaluated under challenging circumstances like low-speed conditions or in the presence of non-Gaussian noise [13,14] which may exist in complex systems. Furthermore, the IESFOgram Targeted approach has been expanded by incorporating a weighting strategy to emphasize bands with higher Targeted features within the 1/3 binary tree levels, enabling the construction of the Combined Improved Envelope Spectrum (CIES) in [15]. In other targeted attempts, He et al. [16] proposed a method for compound fault diagnosis in gears in which the filtered signal's Log-Envelope Auto-Spectrum (LEAS) along with an indicator called Cyclic Feature Index (CFI) was used for frequency band selection. Moreover, a targeted band selection method named Lagged Information Spectrum (LIS) [17] was proposed in which a correlation measure between the signal and the logarithm of its lagged version was defined to detect compound faults targetedly. Nevertheless, a significant constraint arises when there is a lack of prior knowledge about cyclostationary signals of varying intensities but equal importance, necessitating the adoption of a blind approach for compound fault detection.

Therefore, for such a blind investigation of combined faults, the first objective is not to merge the informative bands but rather to individually explore the potential bands to capture any hidden information they may contain. This paper aims to address the challenge of identifying separate frequency bands associated with diverse fault-induced signatures without prior knowledge. Specifically, the aim is to discern faults with weaker signatures, which tend to be obscured in comparison to stronger faults during combined fault scenarios. The rest of the paper is outlined as follows. In Section 2, the background theory of cyclostationarity is described. In Section 3, the proposed methodology, MIESFOgram, is presented. In the methodology, the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) clustering technique is utilized after the feature extraction step and is described in more detail separately in Section 4. Finally, in Section 5, the methodology is evaluated using both simulated and experimental data. The paper closes in Section 6 with some conclusions.

2 Theory of cyclostationarity

Rotating mechanical components have a tendency to generate cyclical transient signatures that exhibit periodicity when the rotational speed remains constant during signal acquisition. These signals often contain valuable information regarding the condition of machine elements, making signal processing and feature extraction crucial for defect detection and monitoring. In line with the principles of cyclostationary theory, the signals acquired from rotating machinery can be effectively characterized by the presence of the first two orders of cyclostationarity. Signals characterized by first-order cyclostationarity (CS₁) exhibit a periodic behavior in their first-order statistical moment as a function of time (T) that conforms to the conditions specified by equation (1): $$C_{1x}\left(t\right)=\mathbb{𝔼}\left\{x\left(t\right)\right\}=C_{1x}\left(t+T\right).$$(1)

In the given equation, the operator 𝔼 represents the ensemble averaging, while t represents time and T represents the period. In the context of rotating machinery, vibration signals characterized by first-order cyclostationarity (CS₁) exhibit periodic waveforms that are associated with specific components and are phase-locked to the rotor speed. Examples of such phenomena include shaft misalignment and spalling on meshing gears. On the other hand, a second-order cyclostationary (CS₂) signal is defined as a signal whose second-order statistical moment is periodic, specifically demonstrated by the periodicity of its autocorrelation function with a period of T, as described in equation (2): $$C_{2x}\left(\tau ,t\right)=E\left\{x\left(t\right)x{(t-\tau )}^{*}\right\}=C_{2x}\left(\tau ,t+T\right).$$(2)

In the context provided, the variable t represents the continuously sampled time, while τ represents the time delay. The symbol 𝔼 corresponds to the expected value or the ensemble average. Bearing vibration signals are commonly classified as second-order cyclostationary (CS₂) signals, due to their underlying periodicity associated with the rotational speed of the shaft. Signals exhibiting higher orders of cyclostationarity, such as CS_n where n is greater than 2, are not typically considered in the analysis of rotating machinery. This is because CS₁ and CS₂ adequately capture the essential characteristics of the signals generated by such machinery.

The Cyclic Spectral Correlation (CSC) is an analytical technique that effectively characterizes signals exhibiting first-order (CS₁) and second-order (CS₂) cyclostationarity in the frequency-frequency domain. This method is formulated as a distribution function involving two frequency variables: the cyclic frequency α, which is associated with modulations, and the spectral frequency f, which is linked to the carrier signal.

The tool can be described also as the correlation distribution of the carrier and modulation frequencies of the signatures present in the signals, defined in equation (3): $$CSC\left(\alpha ,f\right)={\displaystyle \underset{W\to +\infty }{\text{lim}}}\frac{1}{W}E\left\{F_w\left[x\left(t\right)\right]F_w{[x\left(t+\tau \right)]}^{*}\right\},$$(3)

where F_w [x(t)] indicates the spectrum of the signal x(t) over a finite time duration of W using the Fourier Transform. By normalizing the Cyclic Spectral Correlation, the definition of the Cyclic Spectral Coherence (CSCoh) is derived as in equation (4): $$CSCoh\left(\alpha ,f\right)=\frac{CSC\left(\alpha ,f\right)}{\sqrt{CSC\left(0,f\right)CSC\left(0,f+\alpha \right)}}.$$(4)

Both the Cyclic Spectral Correlation (CSC) and the Cyclic Spectral Coherence (CSCoh) bi-variable maps can be integrated over the spectral frequency axis to derive a spectrum. This integration process yields a one-dimensional spectrum function called Improved Envelope Spectrum (IES) which is solely dependent on the cyclic frequency α as shown in equation (5): $$IES\left(\alpha \right)=\frac{1}{F_2-F_1}{\displaystyle \underset{F_1}{\overset{F_2}{{\displaystyle \int }}}}\text{|}CSCoh(\alpha ,f)|df.$$(5)

The Improved Envelope Spectrum via Feature Optimization-gram (IESFOgram) [11] is a recently proposed band selection tool which can be applied on the bi-variable map, as presented in the scheme depicted in Figure 1. The IESFOgram Targeted feature tries to optimize a Normalized Diagnostic Feature (NDF) based on the cyclic characteristics of interest (e.g. rolling element bearing characteristic fault frequencies/orders) on the spectrum resulting from the integration of the bi-variable map. Moreover, the Blind IESFOgram employs some statistical features on the bi-variable map to choose the optimum band. The method is thought to be general enough to be applied to either the CSC or the CSCoh.

Fig. 1 Schematic diagram of the IESFOgram algorithm.

Fig. 2 Schematic diagram of the MIESFOgram algorithm.

3 Proposed methodology

The utilization of bi-variable maps for diagnostic purposes necessitates a thorough understanding of the map to effectively extract and exploit its embedded information. In contrast, one-dimensional spectrum analysis is more prevalent in both academic and industrial settings, due to its simplicity and widespread application. However, integrating the bi-variable function along its spectral variable yields a one-dimensional spectrum, which can serve as a valuable diagnostic tool in its own right. Nevertheless, it is crucial to acknowledge that diagnostic insights within the integrated spectrum can potentially be obscured by the presence of noise and the overlapping signatures of other components. Therefore, careful analysis and interpretation are required to disentangle and uncover the desired diagnostic information from these potential sources of interference. The process of integrating the specific frequency band that contains the desired signal can significantly improve the spectral characteristics and enhance the detection capability for the target frequencies. This approach can effectively boost the performance and the sensitivity in identifying the frequencies of interest. In practical implementations, multiple faults associated with distinct carrier frequencies and varying levels of insensitivity may manifest. Moreover, fluctuating operating conditions, such as varying speed and load, can influence the identification of carrier frequencies across different frequency bands. Consequently, it becomes imperative to adopt a multifaceted approach encompassing multiple frequency bands in order to comprehensively investigate potential characteristic frequencies rather than limiting the analysis to a single optimal band. To tackle this issue, the current research paper introduces a novel approach known as Multiple Improved Envelope Spectra via Feature Optimization-gram (MIESFOgram). The subsequent sections will outline a detailed, step-by-step explanation of the method's process, presented in Figure 2.

Stage 1: In the initial stage, the bi-variable map is derived from the signal using various estimators, such as the Averaged Cyclic Periodogram (ACP), the Cyclic Modulation Spectrum (CMS), the Fast Spectral Correlation (FSC), or alternative numerical techniques, as referenced in [12]. This process facilitates the extraction of the CSC bi-variable map, as described in equation (3). The normalized version of CSC, denoted as CSCoh, can serve as a substitution for the original CSC representation. Notably, these bi-variable maps manifest in the frequency-frequency domains. Furthermore, in the presence of varying speed conditions, these estimators adhere to a similar principle; however, there are some variations. Specifically, the modulation signals will be expressed in the order domain using angular re-sampling methods, while the bi-variable maps will be represented in the order-frequency domains. These adjustments account for the dynamic nature of speed variations and enhance the applicability of the estimators in such scenarios. The users have the option to select their preferred approach for acquiring the bi-variable map based on methods in [18,19]. These references are indicated as numerical implementations of the CSC in the Order-Frequency domain. Also in the reference [9] Fast Spectral Correlation (FSC) was suggested as a provider of the CSC in the Frequency-Frequency domain with good performance. Additionally, when the method is employed on the order-tracked signal with respect to angle, it yields the Order-Order domain.

Stage 2: The subsequent stage involves segmenting the map based on the spectral axis f using a 1/3-binary tree structure, which is analogous to the one employed in the Fast Kurtogram [3]. The division of the map is accomplished through the utilization of individual trees characterized by decreasing bandwidth (Bw) and incrementally varying center frequencies (Cf). These parameters, in turn, establish the upper and lower limits (f₁ and f₂) required for the integration process described in equation (5). Subsequently, the integration of each band within the map yields a demodulated equivalent spectrum known as the Improved Envelope Spectrum (IES).

Stage 3: During this stage, where there is a lack of prior knowledge regarding the characteristic frequencies of faults, Blind features are employed. This approach draws inspiration from the conventional Blind IESFOgram as depicted in Figure 1. The utilized features encompass statistical parameters, including Kurtosis, Spectral Negentropy, Spectral Flatness, L2/L1 Norm, and Gini index. For a comprehensive understanding of the mathematical formulation of these Blind features, further details can be found in the references provided in [13].

This is the stage in the conventional Blind IESFOgram where the optimal band is chosen by maximizing the Blind feature. However, when multiple fault signatures simultaneously are present across different carrier frequencies throughout the spectral band, relying solely on the selection of a single optimal band becomes insufficient. In such cases, the algorithm will fail to capture the frequency bands that exhibit weak fault signatures across various frequency ranges. Subsequent stages will be undertaken to address and resolve this concern.

Stage 4: The objective of this step is to study the frequency bands (which will be called nodes in this paper) in the 1/3 binary tree that takes values based on the Blind Feature BF(Cf_N, Bw_N) in the vector with the length equal to the number of nodes N and Cf_N and Bw_N being the center frequency and the bandwidth of the node N respectively. Therefore, in this stage, the nodes will be sorted (SBF(Cf_N, Bw_N)) based on the Blind feature from the maximum to the minimum. Afterward, all Improved Envelope Spectra (IESs) for all nodes based on the SBF are stored in the IES_SBF.

Stage 5: The core of the solution is applied in this stage, where the correlation matrix is extracted from the IES_SBF for N nodes. The correlation matrix can have the size of (N × N) in the largest size. However, a smaller correlation matrix with size (M × M) from the initial nodes can sufficiently represent the dominant nodes in the spectral bands. Therefore, a value of a_pq (for p ≠ q) in the correlation matrix A, is calculated as follows: $$a_{pq}=\frac{{{\displaystyle \sum }}_{i=1}^{{\alpha }_{\text{max}}}\left(x_{p,i}-{{\displaystyle \overline{x}}}_p\right)\left(x_{q,i}-{{\displaystyle \overline{x}}}_q\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^{{\alpha }_{\text{max}}}{(x_{p,i}-{{\displaystyle \overline{x}}}_p)}^2}\cdot \sqrt{{{\displaystyle \sum }}_{i=1}^{{\alpha }_{\text{max}}}{(x_{q,i}-{{\displaystyle \overline{x}}}_q)}^2}},$$(6)

where x_p,i represents the ith element of vector x_p and x_q,i represents the ith element of vector x_q . Also, the means of vector x_p and x_q are ${{\displaystyle \overline{x}}}_p=\frac{1}{{\alpha }_{\text{max}}}{\displaystyle {\displaystyle \sum }_{i=1}^{{\alpha }_{\text{max}}}}{x}_{p,i}$ and ${{\displaystyle \overline{x}}}_q=\frac{1}{{\alpha }_{\text{max}}}{\displaystyle {\displaystyle \sum }_{i=1}^{{\alpha }_{\text{max}}}}{x}_{q,i}$. The index i is capped by α_max which corresponds to the length of each vector x_p and x_q. It should be noted that p = q shows the diagonal values which are equal to 1 in the correlation matrix. In this way, the correlation matrix of the Multiple Improved Envelope Spectra via Feature Optimization gram (MIESFOgram) is constructed as follows: $$A=\left(\begin{array}{lllllll}\hfill 1\hfill & a_{12}\hfill & a_{13}\hfill & \cdots \hfill & a_{1q}\hfill & \cdots \hfill & a_{1,M}\hfill \\ a_{21}\hfill & 1\hfill & a_{23}\hfill & \cdots \hfill & a_{2q}\hfill & \cdots \hfill & a_{2,M}\hfill \\ a_{31}\hfill & a_{32}\hfill & 1\hfill & \cdots \hfill & a_{3q}\hfill & \cdots \hfill & a_{3,M}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ a_{p1}\hfill & a_{p2}\hfill & a_{p3}\hfill & \cdots \hfill & 1\hfill & \cdots \hfill & a_{p,M}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ a_{M,1}\hfill & a_{M,2}\hfill & a_{M,3}\hfill & \cdots \hfill & a_{M,q}\hfill & \cdots \hfill & 1\hfill \end{array}\right),$$(7)

where A ∊ ℝ ^M×M and M < N. In the correlation matrix of the MIESFOgram, the diagonal elements always have a value of 1, indicating the correlation of an individual IES signal with itself. The off-diagonal elements represent the similarity between different pairs of IES signals. As a result, it is expected that certain rows and columns in the correlation matrix may exhibit low correlation values with most of the IES spectra, except for the intersections with their corresponding pairs in which they demonstrate high similarity. These specific intersection points within the matrix can provide insights into weaker cyclic signatures that may be less prominent compared to stronger components.

Stage 6: To rank the rows with high to low correlation values, the summation of correlation values in each row of the matrix A is calculated and then sorted in a descending way as equation (9): $$A_p={\displaystyle {\displaystyle \sum }_{q=1}^M}{A}_{pq},$$(8) $$SR=\left(\begin{array}{l}{A}_{n1}\hfill \\ A_{n2}\hfill \\ ⋮\hfill \\ A_{nM}\hfill \end{array}\right),$$(9)

where A_n1 ≥ A_n2 ≥ ....

Consequently, the elements of SR will exhibit a noticeable decrease, indicating the presence of weaker signatures that are ranked lower. These values are clearly distinguishable when there exists a significant disparity between the strong and weak signature components, as observed in combined faults where distinct characteristic frequencies manifest in the IES. The dominance of one intensive frequency tends to obscure the other components within the signature.

Stage 7: In this step, the column SR with the sorted values is passed to the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) to be clustered which is explained in Section 4. In this paper, ε is set to 3, and Minpts is determined to 2, as clusters appear to be relatively dense and packed together in all test cases.

Stage 8: During the concluding phase, within each cluster, the node possessing the highest Blind feature among the BF_(1×N) (Cf_N, Bw_N) vector is identified and reported. This yields multiple IES spectra that reveal distinct components in the bi-variable map without prior knowledge. By comparing the identified peaks within these multiple IES spectra to the potential characteristic frequencies within the machine, one can establish correspondences between the signatures and components in the complex system.

4 Density-based spatial clustering of applications with noise (DBSCAN)

The clustering technique of DBSCAN was introduced by Ester in 1996 [20]. It stands out as a pioneering density-based non-parametric clustering algorithm renowned for its proficiency in processing data exhibiting non-spherical structures. In contrast to traditional methods such as k-means and mean shift algorithms reliant on distance metrics, DBSCAN employs a distinctive approach by identifying high-density regions segregated by low-density zones. Each isolated high-density region is subsequently identified as a unique cluster. In other words, what sets DBSCAN apart is its adaptive determination of the cluster count, dynamically adjusting based on the density distribution within the feature space. This adaptability proves advantageous in scenarios where conventional fixed-cluster assumptions are impractical. Therefore, DBSCAN eliminates the need for a priori specification of the number of clusters in the data, unlike k-means.

Generally, two main parameters are introduced into the DBSCAN clustering algorithm, including ε and Minpts. The parameter ε designates the radius of a neighborhood with a specific point and represents the neighborhood range. Parameter Minpts denotes the minimum number of sample points that must be encompassed within the neighborhood with a specified radius. Generally, the points are identified as three kinds: core points, border points, and noise points. Assuming the sample points of P ={ p_i|i = 1, 2, ..., N }, the number of points whose spatial distance from point p_i is not greater than ε can be described as follows: $$N_{ϵ}\left(p_i\right)=\{p_j\in P|dist\left(p_i,p_j\right)\le ϵ\}.$$(10)

It is noted that the definition of “dist” in equation (10) is considered to be the Euclidean distance metric in this paper. As shown in Figure 3, the sample point p_i is identified as a core point when its ε neighborhood encompasses a minimum of Minpnts sample points. i.e., |N_ε (p_i)| ≥ Minpts. Moreover, a sample point is termed a border point if the number of points is less than the Minpts but it is connected to another core point; otherwise, it is labeled as a noise point.

To set ε, the density of the data and the scale of the clusters should be considered. A smaller ε means tighter clusters, while a larger ε means looser ones. In one-dimensional (1D) data where most points are close together, a relatively small value would be chosen to separate two distinct clusters. Also, defining the Minpts parameter for DBSCAN in 1D data can be challenging, because the concept of density and neighborhood is different in 1D compared to higher dimensions. In this case, the distance between data points should be considered along a single line. As an example, when Minpts is configured as 2, it implies that any data point possessing at least one neighbor within a radius of ε will be identified as a core point. Normally, plotting the 1D data and then visually inspecting it may help to get a sense of the distance between points and the potential gaps between the clusters. This can help estimate an appropriate value for both parameters.

Fig. 3 The principle of DBSCAN clustering algorithm.

5 Experimental applications and results

5.1 Simulated signal of combined faults

Within this section, the proposed methodology will be initially evaluated through the application of a simulated signal. The signal generated utilizing a phenomenological model [21], encompasses two distinct fault types, including an inner race and an outer race defect. The first component of the signal is simulated to replicate the occurrence of an inner race fault (BPFI) based on the geometric properties of a deep groove ball bearing identified as SKF 6208. Similarly, the second component is generated to imitate an outer race fault (BPFO) using the specifications of a double-row self-aligning bearing known as SKF 2206 ETN9. The detail information of the bearing's geometries has been listed in Table 1.

In the simulation of the signal, several factors are taken into account. Firstly, the signal was generated at a speed of 8.42 Hz. Moreover, white Gaussian noise was intentionally introduced, simulating the presence of background interference. The amplitudes of the signal are carefully adjusted to ensure that the Enhanced Envelope Spectrum (EES), which represents the integration over the entire frequency band, exhibits peaks of BPFO harmonics but BPFI signatures are below the noise threshold level. It is important to highlight that all figures in this paper include the depiction of the noise threshold, which is calculated based on three times the Moving Absolute Deviation (MAD) of each spectrum [22]. In addition, diverse carrier frequencies are assigned to introduce multiple excited carrier frequencies within the signal. The comprehensive information regarding the modulation and the corresponding carrier frequencies together with the signal specifications can be found in Table 2.

For the signal analysis process, the IESFOgram Targeted methods are applied to examine both the Ball Pass Frequency of the Inner race (BPFI) and the Ball Pass Frequency of the Outer race (BPFO). As illustrated in Figure 4, a modulation by the BPFO at a frequency of 39.77 Hz is identified at the frequency band with central frequency equal to 38.4 kHz and bandwidth 12.8 kHz. Additionally, Figure 5 reveals the identification of a double modulation by the frequency of 45.66 Hz, which corresponds to the BPFI, and by the shaft rotating speed of 8.42 Hz.

Subsequently, the Blind IESFOgram, employing the Gini index, is applied to the signal without any prior knowledge. The Gini index, renowned for its robust performance in the IESFOgram, serves as a statistical feature. For detailed information on the calculation of this feature, the reader is referred to the work cited in [13]. The results are shown in Figure 6. As illustrated in Figure 6a, the assignment of values to the nodes has been meticulously carried out, taking into consideration the presence of the predefined carrier frequencies. However, due to the dominant presence of the BPFO signature in the simulated signal and the conventional approach employed by the Blind IESFOgram, which selects only one optimal band, the alternative band encompassing the faint signature of BPFI remains concealed and unchosen. Consequently, the peaks that are expected to correspond to the BPFI remain undetectable. It is important to emphasize that, at this stage, the conventional Blind IESFOgram chooses the optimal band by selecting the node with the highest Gini index and then blindly extracts the IES of that band without any prior indication of fault frequencies. However, in order to comprehensively assess the method's performance, the peaks at the characteristic frequencies and their harmonics are shown in the IES in Figure 6b.

To address this issue, a solution has been proposed as described in Section 3, wherein the nodes are organized based on maximizing the Blind Gini index values. Subsequently, the IESs were extracted for each node and the corresponding Sorted Blind features SBF and the IES_SBF were calculated. To construct the MIESFOgram Correlation matrix, the user can select a group of initial IESs, as stated in Stage 5 of the methodology. Figure 7a showcases the resultant MIESFOgram Correlation Matrix for the first 50 nodes. An analysis of the figure reveals that a significant portion of the IESs exhibits strong correlations with one another, indicating their classification within a shared cluster of behavior in the bi-variable map. However, it is evident from the MIESFOgram Correlation Matrix that there exist some IESs that display weak correlations with a significant portion of other IESs, as indicated by the darker-colored lines. Interestingly, these lines intersect at specific regions where the correlation values are notably high, suggesting a notable similarity among the corresponding IESs.

Consequently, in accordance with the explanation provided in Stage 6 of the methodology, the summations of each row in the MIESFOgram Correlation Matrix were computed and arranged in descending order. The graphi graphical representation in Figure 7b demonstrates that cluster 1 primarily consists of nodes exhibiting a high correlation, while cluster 2 encompasses nodes with lower correlation values. Importantly, these two clusters are distinctly separated based on the total mean value. In conclusion, as elucidated in Stage 7 of the methodology, the node exhibiting the highest blind Gini index within each cluster was selected based on SBF as the representative of the prevailing components present in the signals. The outcomes of this selection process are presented in Figure 8. Consequently, the identification of the second weaker component of the BPFI occurs in node 448, in contrast to the previously discovered node 34 which exemplifies the BPFO in the simulated signal.

Fig. 4 Targeted IESFOgram to BPFO and IES for a simulated signal.

Fig. 5 Targeted IESFOgram to BPFI and IES for a simulated signal.

Fig. 6 Blind IESFOgram using Gini index and IES for a simulated signal.

Fig. 7 MIESFOgram Correlation Matrix and the sorted row summations in function to the signal index for a simulated signal.

Fig. 8 IES of the nodes with the highest Gini index in each cluster.

5.2 Application of the methodology on a complex system with combined faults

This section presents the evaluation of the methodology using real experimental data acquired from the KU Leuven Diagnostic Test rig, as depicted in Figure 9. The test rig comprised an electric drive motor, a first housing containing a healthy bearing, a disk for the installation of an unbalanced weight, a second housing with large damage to the outer race of the test bearing, a double row self-aligned bearing SKF 2206 ETN9, and two gearboxes. Notably, the input shaft of the first gearbox incorporates a deep groove ball bearing SKF 6208 exhibiting inner race fault damage. The detailed information of the bearings has been mentioned in Table 1. The test rig was operated at shaft rotating speed equal to f_shaft = 8.42 Hz. Two ICP accelerometers (PCB-model number 352A10), were mounted on the housing of the bearing and on the gearbox casing in close proximity to the faulty bearing.

To assess the efficacy of the methodology under diverse conditions, the healthy bearing was mounted in both housings and the remaining components of the drivetrain, encompassing the gearboxes, were interconnected temporarily to facilitate the comprehensive capture of the healthy signal emanating from the second (test) housing. As is observed from the results of Blind IESFOgram using Gini index, depicted in Figure 10, no dominant region in the IESFOgram map presented in Figure 10a is shown which indicates a healthy state in the signal. There is only a very narrow frequency band occupied with electromagnetic interference which shows itself as the harmonics of the shaft speed presented in Figure 10b. Therefore, by applying MIESFOgram Correlation Matrix, only the first value will be related to EMI which has the least correlation with the rest. By applying the DBSCAN clustering on the sorted row summation, all nodes indicating the healthy state are clustered together and only one node indicating the EMI frequency band is chosen as the noise point separately as shown in Figure 11.

Later, the outer race-damaged bearing was mounted in the test housing to capture a faulty case's signal. As shown in Figure 12, now the Blind IESFOgram map (Fig. 12a) is occupied with the cyclic frequencies related to the BPFO and the optimum band of IES in Figure 12b shows BPFO harmonics. Also, as expected after applying the MISFOgram Correlation Matrix and (Fig. 13a) applying the DBSCAN clustering on the sorted row summation all first 50 nodes are clustered as only one cluster which is illustrated in using the Gini index Figure 13b.

The main objective of the experiment is to identify compound faults mentioned earlier by utilizing the signal obtained from the accelerometer positioned on the top part of the housing of the test bearing exhibiting large damage to the outer race. Additionally the drive train containing the gearbox with the inner race-damaged bearing was connected to the test housing. By examining the signal using the Targeted IESFOgram and inputting the frequencies associated with the BPFI and the BPFO, the presence of the BPFO can be identified, as illustrated in Figure 14. Moreover, the BPFI of the bearing within the gearbox can also be determined, as demonstrated in Figure 15.

Nonetheless, when utilizing the Blind IESFOgram with the Gini index, only a single optimal frequency band was selected, resulting in a clear depiction of the BPFO signature in the spectrum, as shown in Figure 16. However, this approach may not be adequate for scenarios involving combined fault cases. It is important to note that in the blind method, the frequencies of the faults are not known in advance. However, for the purpose of demonstrating the method's performance, the figures include data point's locations corresponding to both faults and their harmonics.

To tackle this problem, the MIESFOgram method was applied and the correlation matrix was extracted as shown in Figure 17a. Then, the summation of values in each of the rows was calculated and sorted. The result is shown in Figure 17b which illustrates two separate clusters.

In conclusion, the nodes characterized by the highest Gini index within each cluster were identified and presented as the multiple Improved Envelope Spectra (Multiple IES), as are depicted in Figure 18.

As previously stated, the MIESFOgram algorithm is versatile and can be applied to bi-variable maps derived from various estimators, regardless of whether they belong to the frequency-frequency domain in steady speed conditions or the order-frequency domain in varying speed scenarios. Consequently, the method was also employed on a signal captured under a randomly varying speed profile. The raw time signal and the random varying speed profile are depicted in Figure 19.

Considering that the signal was acquired under varying speed conditions, an angular resampling procedure was performed. Subsequently, an Order Frequency Cyclic Modulation Spectrum (OFCMS) bi-variable map was extracted. The application of the IESFOgram Blind feature Gini index, as illustrated in Figure 20, successfully identifies the optimal band that contains the BPFO component, similar to the steady speed case. Consequently, the subsequent steps involve the application of the MIESFOgram method, resulting in the extraction of the correlation matrix, as demonstrated in Figure 21.

The presence of darker lines exhibiting high correlation intersections suggests the existence of additional cyclostationary clusters. By performing row summation and subsequent value sorting on the correlation matrix, nodes with the highest Gini index pertaining to each cluster were identified and reported. The outcomes, depicting multiple Improved Envelope Spectra (IESs), are presented in Figure 22. As a result of this analysis, the node containing the signatures of the BPFI was also detected, as illustrated in Figure 22b.

Fig. 9 KU Leuven Diagnostic test rig.

Fig. 10 Blind IESFOgram using the Gini index and IES on a healthy signal.

Fig. 11 MIESFOgram Correlation Matrix and the sorted row summations for a healthy signal.

Fig. 12 Blind IESFOgram using the Gini index and IES for one faulty bearing signal.

Fig. 13 MIESFOgram Correlation Matrix and the sorted row summations for one faulty bearing signal.

Fig. 14 Targeted IESFOgram and IES to BPFO on an experimental signal.

Fig. 15 Targeted IESFOmap and IES to BPFI on an experimental signal.

Fig. 16 Blind IESFOgram sing the Gini index and IES on an experimental signal.

Fig. 17 MIESFOgram Correlation Matrix and the sorted row summations for an experimental signal.

Fig. 18 IES of the nodes with the highest Gini index in each cluster.

Fig. 19 Raw time signal and the speed profile of experimental data captured at random varying speed condition.

Fig. 20 IESFOgram extracted from the OFCMS and Blind IES using the Gini index at a signal captured under varying speed.

Fig. 21 MIESFOgram Correlation Matrix and the sorted row summations under varying speed condition.

Fig. 22 IES of the nodes with the highest Gini index in each cluster.

6 Conclusion

This research paper introduces a novel approach called MIESFOgram for the detection of combined bearing faults. The proposed method aims to identify these faults by presenting multiple Improved Envelope Spectra as the final outcomes. Real-world applications typically involve complex machinery comprising diverse bearing and gear components that operate under varying conditions. In scenarios where multiple faults coexist, certain components with subtle signatures may remain concealed, while others exhibiting severe defects can be detected effectively. The MIESFOgram algorithm is developed by exploiting the inherent cyclo-stationary nature of bearing damages, characterized by a periodic autocorrelation function. Consequently, the correlation distribution of carrier and modulation frequencies in bi-variable maps is utilized as an effective tool for detecting cyclo-stationary signals. By integrating the spectral axis of these bi-variable maps, one can successfully extract the Improved Envelope spectra that capture the modulation frequencies of the carriers within the designated frequency range. In recent studies, approaches have attempted to evaluate various frequency bands by blindly examining their statistical characteristics, aiming to identify the optimal band within the bi-variable maps. However, when dealing with combined faults, it has been observed that a single band is inadequate for effectively distinguishing all cyclo-stationarities with varying intensities. The MIESFOgram algorithm incorporates a blind approach to construct the correlation matrix, focusing on the highest statistical features within a specific frequency band. This methodology facilitates the examination of distinct clusters of similar Improved Envelope Spectra (IESs) and, as a result, enables the identification of multiple IESs that encompass diverse cyclo-stationary patterns present in the signal. Subsequently, the analyst gains access to multiple frequency bands, thereby enabling the examination of peak values and their harmonics in relation to characteristic frequencies that may exist within the machinery. This facilitates the identification of faulty components. Importantly, this method can be utilized with any type of estimator to extract the bi-variable data which provides sufficient information about the signal. This research paper presents an evaluation of the proposed method using both simulated and experimental data. The simulation was conducted to incorporate inner race and outer race defects, deliberately designing one defect to be indistinguishable from the other. Furthermore, the proposed method was applied to experimental data obtained from a complex system comprising faulty bearings, gearboxes, and unbalanced components. The evaluation encompassed both the frequency-frequency domain analysis at a constant speed and the order-frequency domain analysis at varying speeds. The results demonstrate that MIESFOgram successfully discerns the presence of a concealed cluster of faults and provides multiple Improved Envelope Spectra (IESs) instead of focusing solely on a single informative frequency band.

Acknowledgments

The authors gratefully acknowledge the Research Foundation Flanders (FWO) under the ROBUSTIFY research grant no. S006119N.

Funding

This research was funded by the Research Foundation − Flanders (FWO) grant number S006119N. The Article Processing Charges for this article are taken in charge by the French Association of Mechanics (AFM).

Conflict of interest

The authors have nothing to disclose.

Data availability statement

Data associated with this article will be made publicly available in the future.

Author contribution statement

Conceptualization, K.G., M.Y.; Methodology, K.G., M.Y., A.M.; Software, M.Y.; Validation, M.Y.; Formal Analysis, M.Y.; Investigation, K.G., M.Y., A.M.; Resources, K.G.; Data Curation, M.Y.; Writing − Original Draft Preparation, M.Y.; Writing − Review & Editing, M.Y., A.M., K.G.; Visualization, M.Y.; Supervision, K.G.; Project Administration, K.G.; Funding Acquisition, K.G.

References

All Tables

All Figures

	Fig. 1 Schematic diagram of the IESFOgram algorithm.
In the text

	Fig. 2 Schematic diagram of the MIESFOgram algorithm.
In the text

	Fig. 3 The principle of DBSCAN clustering algorithm.
In the text

	Fig. 4 Targeted IESFOgram to BPFO and IES for a simulated signal.
In the text

	Fig. 5 Targeted IESFOgram to BPFI and IES for a simulated signal.
In the text

	Fig. 6 Blind IESFOgram using Gini index and IES for a simulated signal.
In the text

	Fig. 7 MIESFOgram Correlation Matrix and the sorted row summations in function to the signal index for a simulated signal.
In the text

	Fig. 8 IES of the nodes with the highest Gini index in each cluster.
In the text

	Fig. 9 KU Leuven Diagnostic test rig.
In the text

	Fig. 10 Blind IESFOgram using the Gini index and IES on a healthy signal.
In the text

	Fig. 11 MIESFOgram Correlation Matrix and the sorted row summations for a healthy signal.
In the text

	Fig. 12 Blind IESFOgram using the Gini index and IES for one faulty bearing signal.
In the text

	Fig. 13 MIESFOgram Correlation Matrix and the sorted row summations for one faulty bearing signal.
In the text

	Fig. 14 Targeted IESFOgram and IES to BPFO on an experimental signal.
In the text

	Fig. 15 Targeted IESFOmap and IES to BPFI on an experimental signal.
In the text

	Fig. 16 Blind IESFOgram sing the Gini index and IES on an experimental signal.
In the text

	Fig. 17 MIESFOgram Correlation Matrix and the sorted row summations for an experimental signal.
In the text

	Fig. 18 IES of the nodes with the highest Gini index in each cluster.
In the text

	Fig. 19 Raw time signal and the speed profile of experimental data captured at random varying speed condition.
In the text

	Fig. 20 IESFOgram extracted from the OFCMS and Blind IES using the Gini index at a signal captured under varying speed.
In the text

	Fig. 21 MIESFOgram Correlation Matrix and the sorted row summations under varying speed condition.
In the text

	Fig. 22 IES of the nodes with the highest Gini index in each cluster.
In the text