Generalized Multivariate Symplectic Sparsest United Decomposition for Rolling Bearing Fault Diagnosis

Abstract

The non-stationary characteristics of the vibration signals of rolling bearings will be aggravated under variable speed conditions. Meanwhile, multichannel signals can provide a more comprehensive characterization of state information, providing multiple sources of information that facilitate information fusion and enhancement. However, traditional adaptive signal decomposition methods generally assume that the frequency information is constant and stationary, and it is difficult to achieve a unified decomposition when dealing with multichannel time-varying signals. Therefore, the intention of this paper is to propose a multichannel signal adaptive decomposition method applicable to variable speed conditions. Specifically, this paper takes advantage of the strong adaptability and robustness of symplectic geometric mode decomposition (SGMD). To improve its applicability to multichannel time-varying signals at variable rotational speeds, a generalized multivariate symplectic sparsest united decomposition (GMSSUD) method is proposed. In GMSSUD, firstly, the completely adaptive projection (CAP) method is employed to achieve a unified representation of the multichannel signals. Then, the generalized demodulation method is introduced to stabilize the signal and subsequently reduce the noise through component screening and reconstruction. Finally, with the new proposed operator as the optimization objective, the constructed sparse filter parameters are optimized to achieve the frequency band segmentation. The analysis results demonstrate that the GMSSUD method possesses higher decomposition precision for multichannel signals with variable speeds and also has a stronger diagnosis ability for variable-speed bearing faults.

1. Introduction

Rolling bearings are key objects in the field of mechanical fault diagnosis. Classical monitoring and diagnosis methods rely on the vibration signals captured by acceleration sensors which are fitted at the measurement points as a medium for sensing the machine conditions [1,2]. With the advancement of technology, the use of multi-sensor data for diagnostic analysis is becoming a trend [3,4].

The well-known empirical mode decomposition (EMD) [5] starts with processing single-channel signals. Taking empirical mode decomposition as the beginning, a series of adaptive signal decomposition methods such as local characteristic-scale decomposition (LMD) [6], variational mode decomposition (VMD) [7], and symplectic geometric mode decomposition (SGMD) [8] have been derived. All these methods can be used to diagnose bearing faults, but EMD and LMD are limited by mode aliasing and obtaining extreme value points, and VMD is limited by parameters. In contrast, SGMD is less disturbed by these factors. Methods based on various parameters and statistical indicators in the time and frequency domains [9] and the commonly used spectral analysis techniques, such as power spectrum, envelope spectrum, and fast kurtosis spectrum analyses [10], can also be used to extract the fault characteristics in the signal, among which the latest product envelope spectrum generated based on spectral correlation has been applied in bearing fault diagnosis [11]. In addition to these traditional signal decomposition algorithms, more intelligent data-driven machine learning methods can often achieve end-to-end processing of data. In recent years, autoencoders, linear discriminant analysis [12], and deep learning represented by convolutional neural networks have emerged [13], and understanding how to reasonably use them for fault diagnosis is also a major development direction of research.

Although single-channel signals are easy to acquire and process, multichannel signals can decrease the uncertainty of the collected signals and can more accurately and comprehensively represent the state information of the device to be measured [14]. References [15,16] improved the EMD and VMD methods according to the characteristics of multichannel signals and introduced the corresponding multichannel signal processing versions MEMD and MVMD, respectively. Among them, MEMD adopts the Hammersley projection method to obtain uniformly distributed direction vectors in the hypersphere. Further, considering the differences in signals in each channel, reference [17] proposed a completely adaptive projection method, which obtained non-uniformly distributed direction vectors depending on the signal power imbalance of each channel, thus obtaining more realistic projection signals. Reference [18] proposed MESMD according to a new multivariate soft screening stopping criterion that adaptively controlled the number of screening iterations of the components. References [19,20] realized the reasonable fusion of multichannel information. Among them, the 2K-FTSSD method [21] has been successfully applied to the decomposition of multichannel signals, which retains the relationship between different channel signals and has good feature extraction ability.

However, in actual working conditions, the equipment is often under the condition of variable-speed operation; during this time, the signal is no longer stable, but becomes a time-varying signal with a changing instantaneous frequency. The diagnostic effect of the above methods will become worse. The time–frequency analysis method that links time information with frequency information in the signal is suitable for signal analysis under variable speeds [22]. With the help of the phase function derived from the frequency curve obtained by the time–frequency analysis method, theoretically generalized demodulation can straighten the time–frequency curve and make it parallel to the time axis. Reference [23] introduced the generalized demodulation method to design a Vold–Kalman filter and decomposed a variable-speed bearing vibration signal into components containing speed information and fault feature information. Therefore, generalized demodulation can be regarded as a tool that is appropriate for removing the variable-speed condition, so as to stabilize the signal and reduce the interference in the signal frequency caused by the change in speed.

Although the above method for variable-speed bearing fault diagnosis has achieved good results, it entails cumbersome single-channel independent processing and necessary subsequent summation when applied to multichannel signals, lacking applicability to multichannel signals with variable speeds. And the general multichannel signal processing methods are susceptible to non-stationary features under variable speed conditions; in this case, the limitations of the method itself are further magnified. At the same time, there are few academic studies on multichannel signal processing under variable speed conditions. Therefore, this paper makes an attempt to process multichannel signals collected under variable speed conditions. This paper exploits the fact that SGMD has the advantages of decomposing a time series into multiple SGMCs with independent modes, not requiring subjective self-defined parameters, and eliminating noise. On the basis of its theory, CAP and generalized demodulation methods are combined, and the GMSSUD method is proposed. This method utilizes the processing capability of the CAP method for multichannel signals and generalized demodulation technology for variable-speed signals, which is suitable for the fault diagnosis of multichannel bearing signals under the condition of variable speed.

In GMSSUD, the CAP method is employed for adaptive projection to obtain the unified representation of multichannel signals first, and the feature information of different channels is aggregated and fused into a single signal. Next, the signal is stabilized by applying the generalized demodulation method to eliminate the interference of variable speed conditions on the subsequent links. Then, the signal after stabilization is denoised through the screening and reconstruction of symplectic geometric atoms. Finally, the proposed regularized singular local linear operator is taken as the objective function of the constructed sparse filter parameter optimization problem to realize the frequency band segmentation, and inverse generalized demodulation is performed on the components located in different frequency bands to obtain the generalized multivariate symplectic sparsest mode components (GMSSMCs). The analysis results indicate that, in comparison with MVMD and MEMD, GMSSUD has better decomposition performance for multichannel signals under variable speeds. Using envelope order spectrum analysis, the GMSSUD method can effectively extract the fault characteristic order and realize accurate rolling bearing fault diagnosis under actual variable speed conditions. The innovations are as follows:

• Taking the proposed regularized singular local linear operator as the objective function, the signal decomposition problem is replaced by the filter optimization problem. While the number of components is determined adaptively, the component signal is constrained into an AM–FM signal with physical significance.
• With the advantage of adaptive acquisition of projection vectors by CAP, the unified characterization from a multichannel signal to a single-channel signal is realized. Meanwhile, the ability of generalized demodulation to convert a time-varying instantaneous frequency into a constant frequency is exploited to achieve the stabilization of variable-speed signals.
• The generalized multivariate symplectic sparsest united decomposition method is proposed. Compared with the MVMD and MEMD methods, the GMSSUD method has better decomposition accuracy for multichannel time-varying signals.

The following is the arrangement of the remaining sections. Section 2 presents the theory related to the GMSSUD method. Section 3 compares the performance of the GMSSUD method through simulation analysis. Section 4 applies the GMSSUD method to practical fault detection. The research conclusions are organized and summarized in Section 5.

2. Theoretical Background

2.1. Completely Adaptive Projection

The completely adaptive projection (CAP) method can effectively solve the problem that the projection vector generated by traditional hypersphere uniform sampling methods, such as the uniform angle sampling method and Hammersley uniform sampling method, lacks adaptability. The fundamental procedure of the CAP method is described below.

Suppose that a multichannel signal can be represented as a set $\left\{y_1\left(t\right),y_2\left(t\right),{\cdots ,y}_d\left(t\right),{\cdots ,y}_D\left(t\right)\right\}$ ($t=\mathrm{1,2},\cdots ,T^{\prime }$, $D$ is the number of channels). Then, this signal is represented by a multivariate signal as $Y\left(t\right)=\left[y_1\left(t\right),y_2\left(t\right),{\cdots ,y}_d\left(t\right),{\cdots ,y}_D\left(t\right)\right]$.

(1)

A series of projection vectors ${\left\{P^i\right\}}_{i=1}^I$ ($I$ is the number of projections, $P^i=\left[p_1^i,p_2^i,\cdots ,p_D^i\right]$) are obtained by sampling on the $\left(D-1\right)$ dimensional sphere using the Hammersley uniform sampling method. The covariance matrix $E\left\{{Y\left(t\right)}^{T}Y\left(t\right)\right\}$ ($E$ is a statistical expectation operator, $T$ is the transpose of a matrix) is constructed, and then the covariance matrix $E\left\{{Y\left(t\right)}^{T}Y\left(t\right)\right\}$ is decomposed by the eigenvalues. The eigenvector $\vartheta$ corresponding to the largest eigenvalue of $E\left\{{Y\left(t\right)}^{T}Y\left(t\right)\right\}$ is defined as the principal direction vector and ${\vartheta }_o=-\vartheta$ as the opposite direction vector. The Euclidean distance between each uniform projection vector $P^i$ and the principal direction vector $\vartheta$ is calculated, and the half of the uniform projection vector ${\left\{P_{\vartheta }^i\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}$ close to the principal direction vector $\vartheta$ is repositioned by Equation (1).

$${\left\{{\tilde{P}}_{\vartheta }^i\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}={\left\{{\displaystyle \frac{P_{\vartheta }^i+\mathrm{a}\vartheta }{\left|P_{\vartheta }^i+\mathrm{a}\vartheta \right|}}\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}$$

(1)

(2)

The half of the uniform projection vector ${\left\{P_{{\vartheta }_o}^i\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}$ close to the principal opposite vector ${\vartheta }_o$ is repositioned using Equation (2).

$${\left\{{\tilde{P}}_{{\vartheta }_o}^i\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}={\left\{{\displaystyle \frac{P_{{\vartheta }_o}^i+\mathrm{a}{\vartheta }_o}{\left|P_{{\vartheta }_o}^i+\mathrm{a}{\vartheta }_o\right|}}\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}$$

(2)

$\vartheta$ and ${\vartheta }_o$ in Equations (1) and (2) are the aggregation directions of the uniform projection vector ${\left\{P_{{\vartheta }_o}^i\right\}}_{i=1}^I$, and $a$ is the aggregation degree. The degree of aggregation is related to the power imbalance between the multichannel signals.

(3)

The power profile of the multichannel signal can be expressed by calculating the power square of each channel of the multichannel signal, and the Gini index can reflect the imbalance in the data well. Therefore, the degree of aggregation $a$ can be determined by calculating the Gini index of the power square vector of each channel of the multivariate signal. $a$ is expressed by Equation (3) as

$$a=1-{\displaystyle \frac{2}{D-1}}\left(D-{\displaystyle \frac{\sum _{d=1}^{D}d\sum _t\frac{1}{T^{\prime }}{\left|y_d\left(t\right)\right|}^4}{\sum _{d=1}^D\sum _t\frac{1}{T^{\prime }}{\left|y_d\left(t\right)\right|}^4}}\right)$$

(3)

A completely adaptive projection vector for the multivariate signal itself can be obtained through the above steps. The projection strategy for CAP is shown in Figure 1.

$${\left\{\widetilde{P^i}\right\}}_{i=1}^I=\left\{{\left\{P_{\vartheta }^i\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}};{\left\{P_{{\vartheta }_o}^i\right\}}_{i=1}^{{\displaystyle \frac{I}{2}}}\right\}$$

(4)

The red and yellow points in Figure 1 represent the principal component direction vector and its inverse direction vector, respectively. When $a=0$, there is a uniform distribution of projection vectors, and the projection degenerates into the Hammersley uniform projection. When $a=1$, the projection vectors are clustered towards the principal component direction vector and its opposite direction vector. CAP can adaptively obtain $a$ according to the features of the signal itself.

2.2. Generalized Demodulation

The essence of the generalized demodulation method is a generalized Fourier transform; assuming the existence of a real signal $x(t)$, the generalized Fourier transform formula is as follows:

$$X_G(f)={\int }_{-\infty }^{\infty }x(t)e^{-j2\pi (ft+s_0(t))}dt$$

(5)

where $s_0(t)$ is defined as the phase function, which is the key to demodulation. Similar to the formula of the inverse Fourier transform, the formula of the inverse transform is

$$x\left(t\right)={\int }_{-\infty }^{\infty }{X}_G\left(f\right)e^{j2\pi \left(ft+s_0\left(t\right)\right)}df=e^{j2\pi s_0(t)}{\int }_{-\infty }^{\infty }{X}_G(f)e^{j2\pi ft}df$$

(6)

Assuming that $X_G(f)=\delta (f-f_0)$, according to Equation (6), there is ${x\left(t\right)=e}^{j2\pi (f_{0}t+s_0(t))}$. Now, the instantaneous frequency of $x\left(t\right)$ is $f_0+s_0^{\prime }(t)$. Assuming that $f_0$ is the initial or fixed part of the instantaneous frequency, $s_0^{\prime }(t)$ is the variable part, and $s_0(t)$ is used for generalized demodulation, the frequency of $x\left(t\right)$ can be stabilized at a fixed value $f_0$, and the time-varying signal can be stabilized.

2.3. Monocomponent AM–FM Annihilating Operator

According to [24], it is clear that the amplitude and frequency modulation (AM–FM) component of the signal can be eliminated by a second-order differential operator; that is, it meets ${{\rm T}}_1\times A\left(t\right)\cos \left(\mathrm{\varnothing }\left(t\right)\right)=0$, where the symbol $\times$ represents the multiplication operation. The operator is shown in Equation (7), where $D^{″}={\displaystyle \frac{d^2}{{dt}^2}},D^{\prime }={\displaystyle \frac{d}{dt}}$, in turn, indicates that the second and first differential operations are performed for the AM–FM components.

$${{\rm T}}_1=D″+P\left(t\right)\times D^{\prime }+Q\left(t\right)$$

(7)

Furthermore, by resolving the ordinary differential equation formed by the aforementioned equation, the two parameters $P\left(t\right)$ and $Q\left(t\right)$ in the operator ${{\rm T}}_1$ are determined, as shown in Equation (8) below, and thereby the monocomponent AM–FM annihilating operator ${{\rm T}}_2$ is formed. The monocomponent AM–FM annihilating operator is actually a kind of AM–FM operator with application significance, because it gives the exact definition of the unknown term in the operator ${{\rm T}}_1$. Annihilation is a description of the function of the operator, which means that when the operator ${{\rm T}}_2$ is multiplied by the AM–FM component, the result is 0.

$$\left\{\begin{array}{l}{{\rm T}}_2=D″+({\displaystyle \frac{-2A^{\prime }\left(t\right)}{A\left(t\right)}}){\times D}^{\prime }+{{\mathrm{\varnothing }}^{\prime }}^2\left(t\right)+({\displaystyle \frac{2A^{\prime }\left(t\right)}{A\left(t\right)}}{)}^2\\ P\left(t\right)={\displaystyle \frac{-2A^{\prime }\left(t\right)}{A\left(t\right)}}\\ Q\left(t\right)={{\mathrm{\varnothing }}^{\prime }}^2\left(t\right)+({\displaystyle \frac{2A^{\prime }\left(t\right)}{A\left(t\right)}}{)}^2\end{array}\right.$$

(8)

where $A(t)$ and $\mathrm{\varnothing }\left(t\right)$ are the amplitude and phase of the AM–FM signals, respectively, and $\mathrm{\varnothing }\prime \left(t\right)$ is the angular frequency.

2.4. Principle of GMSSUD Method

The GMSSUD method first projects the multichannel signals onto the hypersphere to adaptively obtain the projection vectors using CAP, resulting in a unified characterization of the multichannel signals. Next, the generalized Fourier transform is utilized to stabilize the signal. Then, it selects the effective symplectic geometry atoms for merging by the singular local linear operator to realize the pre-noise reduction of the stabilized signal. Subsequently, the parameters of the constructed sparse filter are optimized with the regularized singular local linear operator as the optimization objective, allowing for the segmentation of the frequency bands. Finally, the final signal components are obtained by inverse generalized demodulation.

The overall flow of the GMSSUD method can be divided into four parts, namely multichannel signal projection, signal stabilization, symplectic geometric pre-noise reduction, and frequency band segmentation, and the algorithm flow is shown in Figure 2, with the following complete process:

Step 1: Multichannel signal projection.

(1)

For the multichannel signal $S(n)=[s_1(n),s_2(n),\dots s_l(n)]$, $J$ projection vectors ${\left\{{\widehat{Y}}^{{\theta }_j}\right\}}_{j=1}^J$ are obtained by the CAP projection strategy according to the characteristics of $S(n)$, and then the projection signal ${\left\{p^{{\theta }_j}(n)={\widehat{Y}}^{{\theta }_j}\cdot S(n{)}^T\right\}}_{j=1}^J$ of $S(n)$ is calculated.

(2)

The unified signal $m(n)$ of the multichannel signal is obtained by averaging the $J$ projected signal $p^{{\theta }_j}(n)$.

$$m(n)={\displaystyle \frac{1}{J}}\sum _{j=1}^J{p}^{{\theta }_j}(n)$$

(9)

Step 2: Signal stabilization.

(1)

The unified representation signal $m(n)$ is transformed by STFT. According to the spectrum distribution of the time–frequency (TF) spectrum, the instantaneous frequency (IA) changing spectral line $s_0^{\prime }(t)$ of the signal is obtained by the ridge extraction method.

(2)

The Hilbert transform is applied to $m(n)$ to construct the analytic signal $h(n)=m(n)+jH[m(n)]$.

(3)

The corresponding phase function $s_0(n)$ is obtained according to $s_0^{\prime }(n)$, and the generalized demodulation signal $Z(n)$ is constructed to realize the stabilization of the time-varying signal, where $Z(n)=h(n)e^{-j2\pi s_0(n)}$.

Step 3: Symplectic geometric pre-noise reduction.

The real part $Z_R(n)$ of the signal $Z(n)$ is denoised, where $Z_R(t)=\sum _{k=1}^K{d}_k(t)$.

(1)

For the stabilized signal $Z_R(n)$, pre-noise reduction is performed, so that $r_1\left(n\right)=Z_R\left(n\right)$, and the phase space trajectory matrix $X$ is constructed.

$$X=\left(\begin{array}{c}{r}_1\left(1\right)\\ ⋮\\ r_1\left(m\right)\end{array}\begin{array}{ccc}{r}_1\left(1+\tau \right)& \dots & r_1\left[1+\left(d-1\right)\tau \right]\\ ⋮& \ddots & ⋮\\ r_1\left(m+\tau \right)& \cdots & r_1\left[m+\left(d-1\right)\tau \right]\end{array}\right)$$

(10)

$N$ is the data length, $d$ is the embedding dimension, usually set to ${\displaystyle \frac{N}{3}}$, and $\tau$ is the delay length, $m=N-\left(d-1\right)\tau$. By determining the embedding dimension $d$ and the delay length $\tau$ through PSD, the reconstruction matrix $X$ can be obtained, subsequently obtaining the symplectic geometric atom library after diagonal averaging $Y=\left[y_1\left(n\right),y_2\left(n\right),\cdots ,y_k\left(n\right),\cdots \right]$.

(2)

The singular local linear operator $T$ is computed for each initial single component.

$$T\left[y_k\left(n\right)\right]={{\rm T}}_2(y_k\left(n\right))+{‖y_k\left(n\right)‖}_2^2$$

(11)

where the singular local linear operator $T$ consists of two parts: the operator ${{\rm T}}_2$ is the monocomponent AM–FM annihilating operator defined in Section 2.3, and ${\Vert \Vert }_2^2$ is a 2-norm term. ${‖y_k\left(n\right)‖}_2^2$ is used to evaluate the symplectic geometry atoms, and the larger ${‖y_k\left(n\right)‖}_2^2$ is, the higher the energy of the symplectic geometry atoms and the lesser the residue, ordering $Y$ from largest to smallest by the singular local linear operator of $y_k\left(n\right)$. The sorted atom is defined as $y_k^{\prime }\left(n\right)$ in order to distinguish it from the original symplectic geometric atom. And $Y^{\prime }$ represents the set of all these atoms.

$$Y^{\prime }=\left[{y_1}^{\prime }\left(n\right),{y_2}^{\prime }\left(n\right),\cdots ,{y_k}^{\prime }\left(n\right),\cdots \right]$$

(12)

(3)

$u={\displaystyle \frac{T\left({y_k}^{\prime }\left(n\right)\right)}{T\left(r_1\left(n\right)\right)}}$ is constructed as the reconstruction threshold index, the symplectic geometry atoms with $u>0.001$ are screened for reconstruction, and part of the symplectic geometry atom matrix $Y^{″}=\left[{y_1}^{\prime }\left(n\right),{y_2}^{\prime }\left(n\right),\cdots ,{y_m}^{\prime }\left(n\right)\right]$ containing the significant modes of the original signal is obtained. The remaining ineffective components do not engage in the reconstruction, thereby decreasing the calculation volume and increasing efficiency.

(4)

${y_1}^{\prime }\left(n\right)$ is merged with other symplectic geometry atoms in turn, and the singular local linear operator value is recalculated. If it is reduced, the atoms are merged, and each component ${y_D}_i\left(n\right)$ after symplectic geometric pre-noise reduction is obtained; then, all the components are added to form the denoising signal $y_D\left(n\right)$.

Step 4: Frequency band segmentation.

(1)

The filter $\chi \left(k|\lambda \right)$ is constructed as shown in Equation (13), $\lambda =\left[\omega ,{\omega }_b,{\omega }_c\right]$, and illustrated in Figure 3.

$$\chi \left(k|\lambda \right)=\left\{\begin{array}{ll}\sin \omega \left[k-{\omega }_c+{\omega }_b+{\displaystyle \frac{\pi }{2\omega }}\right],& {\omega }_c-{\omega }_b-\pi /(2\omega )\le k<{\omega }_c-{\omega }_b\\ 1,& {\omega }_c-{\omega }_b\le k\le {\omega }_c+{\omega }_b\\ \cos \omega \left(k-{\omega }_c-{\omega }_b\right),& {\omega }_c+{\omega }_b<k\le {\omega }_c+{\omega }_b+\pi /(2\omega )\\ 0,& others.\end{array}\right.$$

(13)

(2)

A genetic algorithm is used to obtain the optimal filtering parameters $\lambda =\left[\omega ,{\omega }_b,{\omega }_c\right]$ of the denoising signal $y_D\left(n\right)$ to solve the following optimization problem $P1$:

$$P1:Minimize{‖T\left({\widehat{y}}_D\left(n\right)\right)‖}_2^2+{\beta ‖D\prime \left(y_D\left(n\right)-{\widehat{y}}_D\left(n\right)\right)‖}_2^2$$

(14)

In this optimization problem, ${\widehat{y}}_D\left(n\right)=ifft\left[\chi \left(k|{\lambda }_0\right)fft\left(y_D\left(n\right)\right)\right]$, $T\left({\widehat{y}}_D\left(n\right)\right)$ is the singular local linear operator calculation on ${\widehat{y}}_D\left(n\right)$. By minimizing $T\left({\widehat{y}}_D\left(n\right)\right)$, the filtered signal is constrained to be a local narrowband signal, so as to achieve the purpose of adaptive frequency band segmentation. ${\beta ‖D\left(y_D\left(n\right)-{\widehat{y}}_D\left(n\right)\right)‖}_2^2$ is used to regularize the optimization objective function, $D^{\prime }$ is the differential operation that regulates $\left(y_D\left(n\right)-{\widehat{y}}_D\left(n\right)\right)$, and the weight $\beta$ can control the strength of the regularization term, usually set to $0.5{\displaystyle \frac{{‖T({\widehat{y}}_D\left(n\right))‖}_2}{{‖{\widehat{y}}_D\left(n\right)‖}_2}}$.

(3)

Let $y_D\left(n\right)=y_D\left(n\right)-{\widehat{y}}_D\left(n\right)$ and repeat the previous sub-step (2) to build the optimal filtering parameter matrix $F=\left[{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_i,\cdots \right]$ until ${\displaystyle \frac{{‖{\widehat{y}}_D\left(n\right)‖}_2^2}{{‖y_D\left(n\right)‖}_2^2}}\le 1\times {10}^{-6}$, where ${\lambda }_1,$ ${\lambda }_2$, and ${\lambda }_i$ all refer to a certain set of filter parameters in the filter parameter matrix. Finally, the signal after symplectic geometric noise reduction is filtered by the filter bank constructed in $F$ to achieve segmentation of the signal frequency band, and components $y_F\left(n\right)=\left[y_{f1}\left(n\right),y_{f2}\left(n\right),\cdots ,y_{fi}\left(n\right),\cdots \right]$ are obtained, where $y_{fi}\left(n\right)=ifft\left[\chi \left(k|{\lambda }_{\mathrm{i}}\right)fft\left(y_D\left(n\right)\right)\right]$.

(4)

According to the phase function $s_0(n)$ obtained in Step 2, inverse generalized demodulation is performed on the components of different frequency bands obtained after filtering. In this step, the variable speed characteristics in the original signal are restored to each $GMSSMC$, which is convenient for the subsequent envelope order spectrum processing of ${GMSSMC}_s$. The process of obtaining ${GMSSMC}_s$ is as follows:

$${GMSSMC}_s=y_F(n)e^{j2\pi s_0(n)}$$

(15)

3. Numerical Simulation

3.1. Decomposition Performance Comparison

For verifying the decomposition performance of the GMSSUD method in analyzing non-stationary signals, simulation signals $S_1\left(t\right)$, $S_2\left(t\right)$, and $S_3\left(t\right)$ [25] of the three channels are constructed, as shown in Equation (16). Specifically, the first fraction of Equation (16) is executed to form $S_1\left(t\right)$, and similarly to form the remaining two channel signals $S_2\left(t\right)$ and $S_3\left(t\right)$, which together form the three channel signals. Then, it is decomposed using GMSSUD, MESMD, MVMD, and MEMD. Simulation signals are randomly mixed as an attenuation time-varying signal $x_1\left(t\right)$, time-varying signal $x_2\left(t\right)$, and AM–FM signal $x_3\left(t\right)$, as shown in Equation (16), and the simulation signal is presented in Figure 4b.

$$\left\{\begin{array}{l}{x}_1\left(t\right)=1.5e^{-0.8t}\sin (2\pi (-{24t}^2+60t+60))\\ x_2\left(t\right)=2\cos (2\pi ({16t}^2+20t))\\ x_3\left(t\right)=\sin (24\pi t)\cos (500\pi t)\end{array}\right.$$

(16)

$$\left\{\begin{array}{c}{S}_1\left(t\right)=0.3x_1\left(t\right)+0.3x_2\left(t\right)+0.6x_3\left(t\right)\\ S_2\left(t\right)=0.2x_1\left(t\right)+0.5x_2\left(t\right)+0.3x_3\left(t\right)\\ S_3\left(t\right)=0.8x_1\left(t\right)+0.2x_2\left(t\right)+0.7x_3\left(t\right)\end{array}\right.$$

(17)

Because the three multichannel signal decomposition methods, MESMD, MEMD, and MVMD, all have the property of modal alignment, the components with the same ordinal number in the three channels are summed, and the final components of each method are shown in Figure 5, Figure 6, Figure 7 and Figure 8. The GMSSUD decomposition results are basically consistent with the $x_1\left(t\right)$, $x_2\left(t\right)$, and $x_3\left(t\right)$ components in Figure 4a, and the separation of the different modes in the multichannel signals is achieved. The decomposition result of MESMD is also close to the original signal component, but the error between them increases. In Figure 7, only the AM–FM component $x_3\left(t\right)$ of the simulation signal is completely decomposed, i.e., $M\mathrm{I}\mathrm{M}{F}_2$, and the amplitude of the remaining components shows abnormal fluctuations. Meanwhile, in Figure 8, the effect of decomposition is further deteriorated, except that ${mimf}_1$, which represents the AM–FM component, ${mimf}_2$, and ${mimf}_3$ all show obvious mode aliasing. Furthermore, in the quantitative analysis of the decomposition results of the different methods, the signal energy–error ratio (SER) index [17] is used to evaluate the decomposition accuracy of the components. The calculation equation of SER is as follows:

$$SE{R}_i=\sqrt{{\displaystyle \frac{{\sum }_t{(X_i\left(t\right)-Co{m}_i)}^2}{{\sum }_t{(X_i\left(t\right))}^2}}}$$

(18)

$X_i\left(t\right)$ is the original component in the simulation signal; it is actually the sum of the same original signal components in the three channels, for example, $X_1\left(t\right)=0.3x_1\left(t\right)+0.2x_1\left(t\right)+0.8x_1\left(t\right)$. $Co{m}_i$ is the component obtained by decomposition, and $SE{R}_i$ is the SER of the i-th component. It is easy to see that the smaller this index is, the higher the precision of the component. The SER indices calculated for each component are shown in

Table 1. The SERs of the four components ${GMSSMC}_{\mathrm{s}}$, ${MESMC}_{\mathrm{s}}$, ${MIMF}_{\mathrm{s}}$, and ${mimf}_{\mathrm{s}}$.

Method	Com1	Com2	Com3
GMSSUD	$5.99{\times 10}^{-4}$	$9.17{\times 10}^{-3}$	$6.09{\times 10}^{-4}$
MESMD	$2.38{\times 10}^{-3}$	$6.16{\times 10}^{-3}$	0.0562
MVMD	0.724	$1.98{\times 10}^{-4}$	1.58
MEMD	$4.66{\times 10}^{-4}$	0.481	1.10

Com1, Com2, and Com3 sequentially represent the first, second, and third components obtained by each method’s decomposition.

. And the processing time spent in the decomposition of the four methods is also recorded in

Table 2. The processing times of the four methods.

Method	GMSSUD	MESMD	MVMD	MEMD
Processing time (s)	4.88	0.514	2.61	1.426

. It can be seen from Table 1 and Table 2 that the SER values of the components of GMSSUD are almost the smallest, and the processing time is the longest. Therefore, the decomposition accuracy of GMSSUD is higher than that of the other three methods. But the decomposition efficiency is the lowest. The reason for the superiority of the decomposition results of GMSSUD over those of MVMD and MEMD is that GMSSUD is less influenced by modal aliasing, while the superiority over MESMD is due to the proposed new operator, which is capable of constraining the signal components into AM–FM signals.

In order to track the changing trend of the time–frequency information of the multichannel simulation signals, the differences between ${GMSSMC}_s$, $M\mathrm{E}SM{C}_s$, $M\mathrm{I}\mathrm{M}{F}_s$, and ${mimf}_s$ and the original components are visually compared, and their time–frequency spectrograms are drawn, as shown in Figure 9 below.

Comparing the four time–frequency spectrograms, the time–frequency spectrogram of ${GMSSMC}_S$ in Figure 9b is closest to the time–frequency spectrogram of the simulation signals in the ideal state in Figure 9a, and the change trend of their respective spectral lines is basically the same. In the time–frequency spectrogram of $M\mathrm{E}SM{C}_s$ and ${MIMF}_S$, although every spectral line still maintains the general trend of the spectral lines in Figure 9a, the frequency variation rule with time in the two spectral lines below is destroyed. And there is little gap between the analysis results of MESMD and MVMD. Furthermore, the time–frequency spectral lines of ${mimf}_s$ tend to be chaotic and overlap significantly, which proves the mode aliasing of the $mimfs$ component from the side.

Through the above simulation comparative analysis, compared with MESMD, MVMD, and MEMD, the GMSSUD method can effectively separate the time-varying, attenuation time-varying, and AM–FM components in multichannel signals, and the decomposition accuracy is higher, indicating that GMSSUD has a better feature extraction effect.

3.2. Multichannel Variable-Speed Fault Simulation Signal Analysis

Next, the feasibility of GMSSUD in variable-speed fault diagnosis is explored and further simulations are analyzed. To simulate a multichannel rolling bearing fault signal collected under variable speed conditions, the variable-speed rolling bearing fault simulation signal model defined in reference [26] is extended, and the following three original signals are defined:

$$\left\{\begin{array}{l}{x}_1\left(t\right)={\displaystyle \sum _{i=1}^{46}}0.8e^{-800\left(t-T_i\right)}\times \sin \left(4000(t-T_i\right))\times u(t-T_i)\\ x_2\left(t\right)=0.4\cos (2\pi t\times f_r(t)-25)\\ x_3\left(t\right)=0.2\sin (2\pi t\times {2f}_r(t)+15)\\ f_r\left(t\right)=8t+12\\ T_1=\left(1+\mu \right)\times {\displaystyle \frac{1}{n\times f_r(0)}}\\ T_i=\left(1+\mu \right)\times {\displaystyle \frac{1}{n\times f_r\left(T_{i-1}\right)}},i=\mathrm{2,3},\cdots 46\end{array}\right.$$

(19)

$x_1\left(t\right)$ is an outer ring fault model to simulate the impact caused by the fault in the process of variable-speed rotation of the bearing, which is composed of multiple attenuated pulse oscillation signals. The time interval between each pulse varies with time, and according to $T_1$ and $T_i$, the time node at which each fault pulse occurs is determined. $T_1$ is the time node of the first pulse, $T_i$, and so on. $\mu$ is the slip coefficient, set to 0.01; $f_r\left(t\right)$ is the rotation frequency of the variable-speed rotation; and $n$ is the fixed multiple relationship between the outer ring fault frequency and the rotation frequency $f_r\left(t\right)$—here, $n=3$. $u(t)$ is the unit step function to control the occurrence of pulses. The three signal components $x_1\left(t\right)$, $x_2\left(t\right)$, and $x_3\left(t\right)$ are illustrated in Figure 10a.

To simulate the multichannel vibration signals obtained in the acquisition state, the actual multichannel signals are assumed to be a random mixture of signals $x_1\left(t\right)$, $x_2\left(t\right)$, and $x_3\left(t\right)$ in the analysis, and the following random matrix A is used to mix the three original simulation signals [25]:

$$A=\left[\begin{array}{ccc}0.3399& 0.4613& 0.3427\\ 0.2679& 0.2035& 0.8749\\ 0.2585& 0.9970& 0.9514\end{array}\right]$$

(20)

At the same time, since there is noise in the signal collected by the sensor, white Gaussian noise $S_{noise}$ is also introduced into the multichannel analog signals, and the SNR of $S_{noise}$ is set to −10 here. Then, the multichannel simulation signals are obtained, as shown in Figure 9b. The procedure is described below.

$$\left[\begin{array}{c}{S}_1(t)\\ S_2(t)\\ S_3(t)\end{array}\right]=A\times \left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ x_3\left(t\right)\end{array}\right]+S_{noise}$$

(21)

The signal of each channel is directly analyzed by the envelope order spectrum, and the results are shown in Figure 11. In the envelope order subspectrum, the red dashed line is the marked order line at which the order is indicated by the number to its left. The blue circles indicate the points in the envelope order spectrum that are in the labeled order. According to the definition equation of the simulation signal, the corresponding characteristic order is the multiple relationship between the characteristic frequency of the outer ring fault and the rotation frequency, which is 3. There are no obvious spectral peaks at the fault characteristic order in each spectrum; therefore, the single envelope order spectrum analysis cannot extract the fault characteristic information.

Then, the three multivariate signal decomposition methods mentioned above are used for decomposing the simulated signals of the three channels in turn. As in the previous simulation, among the components obtained by the MESMD, MVMD, and MEMD methods, the components with the same ordinal number in the three channels are summed, and the final components of each method are shown in Figure 12, Figure 13, Figure 14 and Figure 15. The processing time spent in the decomposition of the four methods is also recorded in

Table 3. The processing times of the four methods.

Method	GMSSUD	MESMD	MVMD	MEMD
Processing time (s)	27.02	3.28	13.42	8.66

. The MESMD method has the lowest processing time, while the GMSSUD method still has the lowest decomposition efficiency.

The characteristic frequency intensity coefficient (CFIC) [27] is commonly used to evaluate the fault information contained in the signal and its components. Referring to reference [21], the calculation equation for the CFIC is modified in this paper, and the characteristic frequency order intensity coefficient (CFOIC) is proposed. The calculation formula is shown in Equation (22), where $E_o$ is the amplitude of the subspectrum of the envelope order; $f_b$ is the fault characteristic order; $f_1$ and $f_J$ are the start and end orders in the order spectrum, respectively; $Num$ is the total number of orders in the order spectrum; $K$ represents the maximum multiple of the fault order included in the calculation, which is set to 6; ${\sum }_{f_1}^{f_J}{E}_o(f_j)$ refers to traversing all spectral lines from the lowest order to the highest order, calculating the sum of their amplitudes; and $\mathrm{m}\mathrm{a}\mathrm{x}[E_o(i\times f_b-0.02f_b,i\times f_b+0.02f_b)$ refers to the search for the maximum spectral line amplitude in the range of $0.02f_b$, with the fault characteristic order and its multiple as the center.

From this, the CFOIC for each component of the four methods is calculated, as shown in

Table 4. The CFOICs of the four components ${GMSSMC}_{\mathrm{s}}$, ${MESMC}_{\mathrm{s}}$, ${MIMF}_{\mathrm{s}}$, and ${mimf}_{\mathrm{s}}$.

Method	Com1	Com2	Com3	Com4	Com5	Com6	Com7
GMSSUD	0	0.319	1.03	2.04	0.0083	0.0154	0.00085
MESMD	0.051	0.041	0.991	0.104	0.285	0.125	0.224
MVMD	0	0.06	0.0058	0.0019	0.0011	0.000077	0.008
MEMD	0.045	0.039	0.26	0.14	0.088	0.026	0

, where the largest CFOIC value is marked in red, and the components corresponding to the four methods are ${GMSSMC}_4$, ${MESMC}_3$, ${nIMF}_2$, and ${nimf}_3$. Then, these four components are selected for envelope order spectrum analysis. Figure 16 presents the corresponding results, where the red dotted line is used to mark the order of fault characteristics from one to six times.

$$CFOIC={\displaystyle \frac{Num}{K}}\times {\displaystyle \frac{\sum _{i=1}^K\mathrm{m}\mathrm{a}\mathrm{x}[E_o(i\times f_b-0.02f_b,i\times f_b+0.02f_b)]}{{\sum }_{f_1}^{f_J}{E}_o(f_j)}}$$

(22)

In Figure 16a,b, the fault characteristic order (FCO) with a value of 3 and its multiples can be visually observed. But obviously the order spectrum of ${GMSSUC}_4$ is cleaner. In the envelope order spectra of ${MIMF}_2$, there is no fault characteristic order, only the non-obvious multiple fault characteristic order. This cannot provide a sufficient basis for diagnostic analysis. Similarly, in the spectrum of ${nimf}_3$, there is more interference from other spectral lines. Therefore, compared with the other two decomposition methods, GMSSUD has better decomposition performance and can effectively extract the outer ring fault characteristic. At this time, due to the addition of noise, compared with the other three methods, the main reason is that a pre-noise reduction link is added in the process of the GMSSUD method, so as to have better anti-noise performance.

4. Experiment

In order to investigate the validity of GMSSUD in actual rolling bearing fault diagnosis under variable speed conditions, a bearing test bench, as shown in Figure 17, was built for experimental verification. The test bench was composed of a drive motor, experimental bearing, support bearing, coupling, acceleration sensors, and a B&K Pulse signal acquisition system. The variable-speed rotation of the motor was controlled by a PLC.

Three acceleration sensors were mounted on the top and sides of the experimental bearing housing by magnetic adsorption. In the case that the sampling theorem was satisfied and the amount of data processed was taken into account, we set the sampling rate to 65,536 Hz. When the B&K Pulse signal acquisition system is working normally, the bearing vibration data from the three channels of the three sensors can be recorded at the same time. The experimental bearing type was 6206, and

Table 5. Specification parameters of the experimental bearing.

Parameter	Pitch Circle Diameter (mm)	Contact Angle $(°$)	Ball Diameter (mm)	Number of Balls
Value	46	0	9.42	9

presents the parameters of the bearing. The outer ring fault was set artificially by EDM cutting, and the cutting depth was 0.2 mm. The fault parts are also shown in Figure 17. Based on the empirical formula, the multiple relationship between the outer ring fault characteristic frequency and the rotation frequency of the experimental bearing could be calculated, that is, the fault characteristic order (FCO)—here, it was 3.58.

(a)

Outer race fault

The motor speed changed approximately linearly and became stable after reaching the maximum speed of 1560 rpm for 36 s. The signal collection started from the motor startup to the end after keeping 1560 r/min running for 3 s, and a total of 39 s vibration signal data was collected. The data from 16 s to 25 s were intercepted for analysis. The display screen circled in yellow in Figure 17 can provide the current rpm in real time, and with the use of a PLC to adjust the speed of the test bench, the change law of the speed can be manually set. The speed curve obtained by fitting the actual speed is shown in Figure 18.

The vibration signals of the three channels with a total of 10 s are shown in Figure 19. The experimental signals of the three channels were analyzed by the envelope order spectra, as shown in Figure 20. Among them, no spectrum diagram shows obvious fault order information. The GMSSUD, MESMD, MVMD, and MEMD methods were adopted for decomposing the above signals. The processing time spent in the decomposition of the four methods is recorded in

Table 6. The processing times of the four methods.

Method	GMSSUD	MESMD	MVMD	MEMD
Processing time (s)	41.42	6.54	26.74	13.66

. Then, the CFOIC of the components obtained by the four methods was calculated successively, as shown in

Table 7. The CFOICs of the four components ${GMSSMC}_{\mathrm{s}}$, ${MESMC}_{\mathrm{s}}$, ${MIMF}_{\mathrm{s}}$, and ${mimf}_{\mathrm{s}}$.

Method	Com1	Com2	Com3	Com4	Com5	Com6	Com7
GMSSUD	0.0421	2.1517	0.0566	0.0029	0.0085	0.0123	0.038
MESMD	0.0480	0.0417	0.0329	0.0606	0.0836	0.117	0.7349
MVMD	0.4293	0.0162	0.0049	0.002	0.0011	0.001	0.0024
MEMD	0.0531	0.0717	0.0616	0.0125	0.0175	0.0038	0

. The value of $K$ in the CFOIC equation was set to 5. For brevity, only the time-domain waveforms of the first four components with large CFOIC values are shown in Figure 21. Finally, the component with the largest CFOIC value in each method was analyzed by the envelope order spectra, and the results are illustrated in Figure 22, where subfigures (a), (b), (c), and (d) correspond to the envelope order spectra of components ${GMSSMC}_2$, ${MESMC}_7$, ${MIMF}_1$, and ${mimf}_2$, respectively. The labels to the left of each red dotted line in the envelope order spectra represent the order value of the corresponding light blue circle mark point, where RFO represents the order of rotation frequency, 2RFO represents the order of double rotation frequency, and so on. FCO represents the fault characteristic order, specifically the outer ring fault characteristic order with a value of 3.58, and 2FCO is twice its order.

Figure 22a contains an obvious fault characteristic order (FCO) and rotation frequency order (RFO). In addition, there is rich order information of multiple frequency doubling, indicating that the GMSSUD method extracts fault characteristics in the signal. Figure 22b also shows an abundant RFO and FCO. Although the extraction effect is not as good as that in Figure 22a, the existence of the outer ring fault can also be judged. However, in Figure 22c, only the non-obvious RFO is shown. Similarly, the amplitude of the FCO in Figure 22d is similar to the average amplitude of the entire 0–5 order, which is interfered by the irrelevant order, so its credibility as fault characteristic information is doubtful. The above analysis indicates that compared with MESMD, MVMD, and MEMD methods, the GMSSUD method can accurately extract the outer fault characteristic information in the vibration signal under variable speed conditions. Since the acquired signal is mainly affected by noise and the changing frequency rotation, the increased time-varying characteristics of the signal will interfere with the decomposition effect of the MESMD, MVMD, and MEMD methods, while GMSSUD performs generalized demodulation first and then decomposition, which can reduce the interference of time-varying characteristics in the decomposition process.

(b)

Inner race fault

The experimental bearing in the test bench was replaced with an inner ring fault bearing under the same fault processing conditions; the bearing model and parameters were unchanged, and the other experimental conditions were also unchanged.

The motor speed changed approximately linearly and became stable after reaching the maximum speed of 3094 rpm. The data from 15 s to 24 s were intercepted for analysis. The speed curve obtained by fitting the actual speed is shown in Figure 23. The vibration signals of the three channels with a total of 10 s are shown in Figure 24. The remaining process was consistent with the experimental analysis of the outer ring fault. Figure 25 shows the envelope order spectra of the experimental signals for the three channels, from which very weak frequency conversion orders can be observed. The processing time spent in the decomposition of the four methods is recorded in

Table 8. The processing times of the four methods.

Method	GMSSUD	MESMD	MVMD	MEMD
Processing time (s)	45.93	7.21	25.89	12.52

. The CFOIC of the components obtained by the four methods was calculated, as shown in

Table 9. The CFOICs of the four components ${GMSSMC}_{\mathrm{s}}$, ${MESMC}_{\mathrm{s}}$, ${MIMF}_{\mathrm{s}}$, and ${mimf}_{\mathrm{s}}$.

Method	Com1	Com2	Com3	Com4	Com5	Com6	Com7
GMSSUD	0.1212	0.8013	2.9452	1.1722	0.0819	0.3320	0.1375
MESMD	0.0891	0.0666	0.1034	0.1619	0.4570	0.8757	1.1910
MVMD	0.2411	0.0502	0.0126	0.0053	0.0016	0.0015	0.0011
MEMD	0.1404	0.1628	0.1711	0.5573	1.3163	0.8576	0

. The value of $K$ in the CFOIC equation was set to 3. The decomposition results and the envelope order spectra of the components with the largest CFOIC value are shown in Figure 26 and Figure 27, respectively, where subfigures (a), (b), (c), and (d) correspond to the envelope order spectra of components ${GMSSMC}_3$, ${MESMC}_7$, ${MIMF}_1$, and ${mimf}_5$, respectively. The labels to the left of each red dotted line in the envelope order spectra represent the order value of the corresponding light blue circle mark point, where RFO represents the order of rotation frequency, 2RFO represents the order of double rotation frequency, and so on. FCO represents the fault characteristic order, specifically the inner ring fault characteristic order with a value of 5.442, and 2FCO is twice its order.

Figure 27a contains obvious order information, and the RFO forms a side frequency order similar to the side frequency band on both sides of the FCO. Different from the outer ring fault experiment, the effect of the envelope order spectrum on ${MESMC}_7$ is similar to that of ${MIMF}_1$ and ${mimf}_5$, and the obvious rotation frequency order and fault characteristic order cannot be extracted. In the same case, the difficulty of inner ring fault detection is generally higher than that of outer ring fault detection. Therefore, the fault feature extraction effect of the GMSSUD and MESMD methods is worse, but the GMSSUD method can still extract fault features. The increase in the trend of speed change further affects the stationarity of the collected signal, so the effect of the MVMD and MEMD methods is also worse. This indicates that compared with the MESMD, MVMD, and MEMD methods, the GMSSUD method can accurately extract the inner fault characteristic information in the vibration signal under variable speed conditions and can effectively realize the fault diagnosis of rolling bearings when combined with the envelope order spectra.

5. Conclusions

In the simulation analysis, the GMSSUD method had higher decomposition accuracy and fault feature extraction ability, and in the two experiments, its cooperation with the envelope order spectra enabled it to diagnose the faults of the bearing outer ring and inner ring, which indicates that GMSSUD is more suitable for dealing with signals of variable frequency compared with the other three methods. The biggest reason is that generalized demodulation technology is employed to preprocess the signal before the decomposition, which reduces the interference of time-varying features in the decomposition process. In general, the GMSSUD method combined with envelope order spectrum technology has application potential in multichannel bearing fault diagnosis under variable speed conditions. The specific conclusions are as follows:

(1)

The CAP method is adopted to process multichannel signals to obtain a unified representation signal. This method can obtain the projection vectors adaptively according to the characteristics of the multichannel signals and realize the feature enhancement by fusing the projection signals from different channels.

(2)

Generalized demodulation is adopted to stabilize the variable-speed signal. This method can straighten the time–frequency curve and make it parallel to the time axis, which removes the interference in the instantaneous frequency caused by the change in speed.

(3)

The problem of signal decomposition is transformed into the problem of sparse filter parameter optimization, taking the regularized singular local linear operator as the optimization objective, and by optimizing the parameters of the constructed sparse filter, the decomposition result is constrained to be an amplitude–frequency modulation. The decomposition results have better physical significance.

In spite of this, there are still some shortcomings of GMSSUD. Firstly, the design of the variable speed conditions in this paper was relatively simple, especially in the experimental signal. So, when the variable speed condition becomes complex, understanding how to ensure the reliability of the extracted time–frequency curve will be important. Secondly, the processing time of GMSSUD is obviously longer than that of the comparison methods, and the low decomposition efficiency will restrict the application of this method. So, in the subsequent work, these two problems need to be resolved. For long algorithm processing times, we plan to improve the whole process of the algorithm to shorten the consumption time of each process and replace the algorithm with better performance optimization. However, the problem of insufficient accuracy of the extracted time–frequency curves still needs to be studied to find a solution.