Multivariate Dynamic Mode Decomposition and Its Application to Bearing Fault Diagnosis

In practical engineering applications, the multivariate signal contains more fault feature information than the single-channel signal. How to realize synchronous extraction of fault features from the multivariate signal is of great significance in fault diagnosis of rotary machinery. Dynamic mode decomposition (DMD) has attracted much attention due to its excellent dynamic feature extraction ability. However, DMD lacks mode aliasing property when dealing with the multivariate signal, which may lead to the loss of critical fault feature information. Cater to this problem, this article proposed a multivariate DMD (MDMD) algorithm that is the multivariate extension of DMD. First, snapshot tensors are defined to convert multivariate signals to tensor format. Then, the MDMD algorithm is proposed by introducing tensor operations into the original DMD algorithm, where a tensor low tubal rank component extraction framework is constructed to enable simultaneous extraction of bearing fault features from the multivariate signal, to enhance the robustness and effectiveness of the algorithm. Finally, both numerical simulations and experiments verify that the proposed MDMD has higher fault diagnosis accuracy than other multivariate signal-processing methods.


I. INTRODUCTION
F AULT diagnosis is essential to ensure the stable operation of machinery [1], [2]. As a key component of rotating machinery, rolling bearings are prone to failure due to heavy loads and extreme temperatures during operation [3], [4], [5]. Once the bearing fails, the maintenance cost of the equipment will increase due to unexpected downtime and even cause catastrophic accidents [6], [7], [8]. Therefore, early fault diagnosis for rolling bearings has attracted much attention from researchers.
The early fault of the bearing usually behaves as a single failure, and the periodic impulses will be generated when the rollers pass the damaged part [9]. Due to the complex transmission chain of the mechanical system, the periodic impulses are easy to be coupled with the periodic components generated by other components (such as shaft misalignment, gear meshing, etc.), which makes the collected signals show obvious nonlinear characteristics [10]. Meanwhile, the fault features are often submerged in strong noise, which brings great difficulty to fault feature extraction. Therefore, it is the main task of bearing fault diagnosis to accurately extract fault feature components from noise signals. In recent years, many feature extraction methods have been proposed to improve the accuracy of fault diagnosis. Signal filtering methods, such as spectral kurtosis [12] and empirical wavelet transform [13], can extract fault features by identifying fault-related This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ frequency bands. In addition to the above methods, the cyclostationary analysis methods [14], [15] can display the energy spectrum that is composed of spectral frequency and cycle frequency. These methods can effectively extract fault characteristic frequencies through envelope analysis techniques.
Signal decomposition methods, such as empirical mode decomposition [16], local mean decomposition [17], variational mode decomposition [18], and dynamic mode decomposition (DMD) [19], decompose the signal into fault-related modes to extract fault features. The above methods provide a rich theoretical basis for bearing fault diagnosis. However, considering the practical application, the location of the fault cannot be known in advance, and the prearrangement position of the sensor may not be optimal, which may lead to weak fault features contained in the collected signal. Therefore, it is a common operation in practical engineering applications to use multiple sensors to collect multichannel signals to increase the collection of fault-related information. However, the single-channel fault diagnosis method cannot efficiently utilize the collected signals due to the uneven energy distribution of fault features among multiple channels. Hence, it is of great significance to develop a high-accuracy and robust multivariate signal-processing method for bearing fault diagnosis. How to realize information fusion between channels, extract fault features synchronously, and improve the robustness of the algorithm are the key issues of multivariate signal-processing methods.
To solve the above problems, researchers extended some mode decomposition algorithms to the field of multivariate signal processing. Multivariate empirical mode decomposition (MEMD) [20] realizes information fusion between channels by constructing a multidimensional space using a uniform projection-based approach. Then the local mean is estimated to fulfill the simultaneous decomposition of multivariate signals at different frequency scales. With the merits, MEMD and its improved methods [21], [22], [23] have been widely employed in mechanical fault diagnosis. However, as the multivariate extension of the EMD method, MEMD inherits the limitations of single-channel EMD. In addition to mode mixing and endpoint effects, the lack of noise robustness also restricts its engineering applications. Therefore, how to enhance the noise robustness of MEMD has aroused the interest of researchers. Lv et al. [24] improved the noise robustness of MEMD using nonlocal means denoising approach. They verified that MEMD can effectively extract bearing multivariate fault features. Hao et al. [25] proposed a joint denoising framework for MEMD. The denoising performance of MEMD is improved by using a subspace projection scheme. Multivariate variational mode decomposition (MVMD) [26] extracts multivariate modulation oscillations synchronously from the Fourier domain and thus directly has mode alignment property. In other words, the mode alignment property guarantees the accuracy of MVMD in fault diagnosis when all channels contain common fault information. However, the performance of MVMD depends on the selection of two key parameters. Song et al. [27] proposed a self-adaptive MVMD method for bearing fault diagnosis. They determined the mode parameters and initial center frequencies based on the convergence trend of MVMD. Lu et al. [28] proposed an adaptive parameter selection method for MVMD based on the difference in correlation coefficients.
The above research provides valuable ideas for solving the key problems of multivariate signal processing. In recent years, DMD has been widely used in various fields [29], [30], [31], [32] due to its complete theoretical foundation and excellent dynamics feature extraction capability. DMD characterizes the dynamics of the original system by constructing a proxy dynamical system matrix and decomposing it into single-frequency DMD modes to extract the dynamical features hidden in the signal. However, the standard single-channel DMD lacks the mode alignment property. The energy-based DMD mode selection strategy cannot accurately locate fault features due to the uneven energy distribution of bearing fault features between channels. To further solve the problems of information fusion and simultaneous extraction of bearing fault features, a novel multivariate DMD (MDMD) algorithm is proposed in our research here. First, we introduce the tensor product into the operation of MDMD and realize the information fusion between channels by using the characteristic of circular convolution of the tensor product in the time-domain signal. Then, a robust tensor low-rank component extraction framework is constructed in MDMD for the simultaneous extraction of bearing fault features from the multivariate signal, which enhances the noise robustness of the proposed method. In our previous work [33], we have verified that the fault features are low-rank distributed in the dynamical system matrix. As the multivariate extension of the traditional DMD algorithm, multivariate fault features are also lowrank distributed in the dynamic system tensor (DST). Finally, numerical simulations verify the mode alignment property of the proposed approach, while the mode alignment rate (MAR) of the proposed method is higher than that of MEMD and MVMD. The subsequent experiments further verified the conclusion. The innovations of this article are summarized as follows.
1) A novel MDMD algorithm is proposed. The tensor singular value decomposition (TSVD) and tensor product are constructed during the operation of MDMD. It is worth noting that the proposed MDMD has mode alignment property due to the nature of circular convolution operations of the tensor product.
2) A novel robust low-rank tensor component extraction framework is constructed in MDMD for multivariate fault feature extraction. The proposed framework not only enhances the noise robustness of the algorithm, but also avoids the problem of multivariate DMD mode selection.
3) A bearing fault diagnosis method for multivariate signals based on MDMD is proposed. The snapshot tensor is defined to convert the multivariate bearing fault signal into a tensor format. Then the proposed method provides a robust fault feature extraction for multivariate fault signals. Numerical simulation and experimental analysis verify the effectiveness of MDMD in bearing fault diagnosis.
The remaining structure of this article is organized as follows. Section II provides the theoretical background of DMD. Section III describes the proposed MDMD algorithm. Section IV illustrates the mode alignment property of the proposed MDMD algorithm. Section V verifies the effectiveness of the proposed multivariate signal-processing method in bearing fault diagnosis, and the conclusions are drawn in Section VI.

II. DYNAMIC MODE DECOMPOSITION
In this section, the basic idea of standard DMD is introduced so that it can be extended to the field of tensor computing later. The DMD algorithm provides a spatiotemporal decomposition of 1-D time series into DMD modes so that the low-rank structure in a complex system can be extracted by combining specific modes. A sampling window divides the 1-D time series into snapshots as [x 1 , x 2 , . . . , x n ], and the time interval between the two snapshots is t, DMD assumes that the nonlinear system is approximately linear in time t, and the actual dynamic system can be characterized by a high-dimension dynamical system matrix A as The optimal dynamical system matrix A is given by where X † is the pseudoinverse of X . Then, an singular value decomposition (SVD)-based pseudoinverse solving method is introduced below. Assume that X = U V T , where is the diagonal singular value matrix. The pseudoinverse of X can be calculated as follows: To increase the efficiency of the algorithm, it is more convenient to computeÃ by projecting a dynamical system matrix A on UÃ AssumingÃ = W is the eigendecomposition ofÃ, where and W are the eigenvalue matrix and the eigenvector matrix, respectively. The DMD modes are defined as follows: The 1-D time series can be reconstructed by where b = † x 1 is the initial conditions of the reconstruction process. The pseudo-code of standard DMD is shown in Algorithm 1.

III. MULTIVARIATE DYNAMIC MODE DECOMPOSITION
In this section, the MDMD algorithm is proposed to process multichannel or multivariate signals synchronously. First, the tensor product is introduced as the basis for the multivariate DMD. Second, the matrix operations in standard DMD are rewritten as tensor operations. In particular, the detailed algorithms for the computation of tensor pseudoinverse and DMD modes are discussed. Afterward, a tensor low tubal rank component extraction framework is introduced to fulfill the automatic identification of DMD modes. Finally, the algorithm of MDMD is given at the end.

A. Notations and Preliminaries
In this article, tensors are represented by Euler script letters, for example, X ; matrices are represented by capital letters, for example, X ; vectors are represented by lowercase letters, for example, x; and scalars in tensor are denoted by x i jk , where x i jk denote X (i, j, k) in MATLAB notation. The inner product between X and Y is defined as ⟨X, Y ⟩ = T r (X * Y ). The inner product between X and Y is defined Some norms of the tensor are used below. The ℓ 1 -norm of the tensor is ∥X ∥ 1 = i jk |x i jk |; the infinity norm of the tensor is ∥X ∥ ∞ = max i jk |x i jk |; the Frobenius norm of the tensor is ∥X ∥ F = ( i jk |x i jk | 2 ) 1/2 ; the tensor nuclear norm Definition 1 (TSVD [34]): Let X ∈ R n 1 ×n 2 ×n 3 . Then it can be factorized as where U ∈ R n 1 ×n 1 ×n 3 , V ∈ R n 2 ×n 2 ×n 3 are orthogonal, and S ∈ R n 1 ×n 2 ×n 3 is an f-diagonal tensor. Definition 2 (Tensor Product [35]): Let X ∈ R n 1 ×n 2 ×n 3 and Y ∈ R n 2 ×l×n 3 . Then the tensor product X * Y is defined as where fold(X ), unfold(X ), and bcirc(X ) are defined as It can be seen from Definition 2 that the tensor product in the time domain is actually the circular convolution between the layers, while it is equivalent to the matrix multiplication in the Fourier domain. Benefiting from this property, the tensor product operations are performed in the Fourier domain to improve the operational efficiency of the MDMD algorithm.

B. Multivariate DMD
Before formulating the MDMD algorithm, a multivariate signal needs to be constructed as snapshot tensors. Assume a set of n-channel signal S = [s 1 (t), s 2 (t), . . . , s n (t)] ∈ R n×m . First, the snapshot matrices are constructed by using a sampling window to synchronously intercept the data of different channels. Then the snapshot tensors are obtained by horizontally stacking the snapshot matrices. It is worth noting that signals from the same channel should be placed in the same forward slice to activate the information fusion property of tensor operations. Therefore, proper rotation of the snapshot matrix is necessary. The detailed steps are shown in Fig. 1. Similar to standard DMD, the two snapshot tensors X and Y are constructed with a delay of t. The DST of MDMD is defined as Next, a TSVD-based tensor pseudo-inverse framework is constructed, and the detailed algorithm is shown in

Algorithm 2 Calculation Steps of the Tensor Pseudoinverse
Input: tensor X ∈ R n 1 ×n 2 ×n 3 . 1 Algorithm 2. It can be observed from Algorithm 2 that all operations are done in the frequency domain. It is the convolution in the time domain that can be converted into the product in the frequency domain that the operation can speed up the efficiency of the algorithm. Meanwhile, the layer-by-layer inversion of the diagonal matrix S in the frequency domain simplifies the steps of the conventional tensor unfolding.
The projected DST B can be calculated as The multivariate dynamic modes are defined as the eigenvector of B. To calculate MDMD modes, we also use the properties of the tensor product to perform the eigendecomposition of frontal slices in the Fourier domain. Thus, the tensor B can be factorized as where E is an f-diagonal tensor and W can be viewed as a tensor whose frontal slices are composed of eigenvectors.
Thus, the MDMD modes can be calculated as The procedure for computing the MDMD modes is given in Algorithm 3.
To this end, we have completed the formulation of MDMD and realized the synchronous decomposition for multichannel signals. It is worth noting that MDMD behaves as a slice operation on tensors in the Fourier domain, which may cause a misunderstanding that the proposed method does not involve data fusion between multichannel signals. But the fact is that the tensor product in MDMD can be viewed as the circular convolution of the tensor slices. Hence, the proposed MDMD is able to extract fault feature information from the multichannel signals, synchronously.
Another point to highlight is the identification of multivariate DMD modes. In standard DMD, the low-rank dynamical model is constructed through specific combinations of DMD modes. Compared with the single-channel DMD modes, multivariate DMD modes are obtained in a matrix form, which increases the difficulty of bearing fault feature identification. In this article, a tensor low tubal rank component extraction framework for MDMD is proposed, which is aimed at directly extracting fault feature components from DST and avoiding the selection of multivariate DMD modes. The tensor low tubal rank component extraction problem in MDMD can be written as where B L represents the low tubal rank component in DST and B N represents the noise that is sparsely distributed in tensors. The joint optimization problem is known as the tensor robust principal component analysis [36] model which can be solved via the alternating direction method of multipliers (ADMMs). The augmented Lagrangian function of (15) can be written as

Algorithm 4 Tensor Low Tubal Rank Component Extraction
Input: B ∈ R n 1 ×n 1 ×n 3 , parameter λ. Initialization: where Y is the Lagrange multiplier. B L and B N can be updated by iteratively minimizing (16), while fixing the others.
Step 1: Update B L k+1 Step 2: Update B L N +1 Step Step 4: Update µ k+1 The initial parameter settings and iteration-stopping conditions are presented in Algorithm 4. The regularization parameter λ controls the process of low-rank approximation, which is related to the noise in DST. The error ε is the stopping condition of the iteration, which is related to the accuracy and running time of the algorithm. The recommended settings for the above parameters can be referred to [36].
To sum up, the framework of the MDMD algorithm is given in Algorithms 2 and 3, which extends DMD to multivariate signal processing. Algorithm 4 is designed to extract low tubal rank components from DST directly, which increases the noise robustness while avoiding the mode selection problem. It is worth noting that the low tubal rank component extraction is proposed for fault diagnosis of rotating machinery, which has a prior that the fault component is low rank [28]. The whole algorithm of the proposed MDMD is given in Algorithm 5.

IV. NUMERICAL SIMULATION
In this section, the mode alignment property of MDMD is discussed. Then, the comparative studies with representative multivariate signal-processing methods (MEMD and MVMD) are proposed.

A. Research on the Mode Alignment Property of MDMD
The mode alignment of a multichannel signal refers to the alignment of joint oscillations across different channels [26]. In other words, similar vibration frequencies are present in each channel of a multichannel signal. For bearing fault diagnosis, the mode alignment property of signal-processing methods is of fundamental importance to ensure physically meaningful decisions. However, the arrangement of the sensors and noise interference may cause energy imbalance between channels during the acquisition of the signal. Taking the signal collected by the multichannel sensor as an example, the energy levels of fault components between each channel are different, but the energy levels of noise interference between channels are the same. Therefore, the bearing outer-race fault simulation model [37] is introduced to simulate the multichannel fault signal where x 1 (t), x 2 (t), x 3 (t) represent the signal collected from the three channels. The amplitudes of the three channels are set to [1], [2], [4] to simulate energy imbalance. The parameters are shown in Table I.
The simulated multichannel signal can be expressed as is the clean signal and N is the white Gaussian noise (WGN) with the signal-to-noise ratio (SNR) = -5 dB. The waveform and the spectrum of the simulated signal are shown in Fig. 2. It can be observed from the spectrum that only the fault characteristic frequency (FCF) can be identified in Channel 1, while the FCF in Channels 2 and 3 are masked by noise.
First, MDMD is employed to process the above-simulated signal. It can be observed from Fig. 3 that the FCF( f o ) can be identified in all three channels. It is worth noting that, although the SNR of the multichannel signal is −5 dB, the SNR of Channel 3 is much smaller than −5 dB. It may cause a misunderstanding that more WGN is added to Channel 3 than the other two channels. But the fact is that the amplitude of the   vibration signal collected by Channel 3 is lower than the other two channels, which is consistent with the actual situation. Afterward, DMD is employed to process three single-channel signals. Fig. 4 is the result of standard DMD. Obviously, the interference frequency and noise still exist in Channels 2 and 3, which makes it difficult to identify the FCF and its multiples.
Next, a novel FCF detection indicator is proposed to evaluate the mode alignment property of MDMD, which is based on the idea of outlier detection in statistics. In the processed signal, the FCFs are inevitably mixed into the interference frequencies and even submerged by the noise. But it should be noted that the FCFs are sparsely distributed in the frequency domain, and the amplitude of the noise is relatively low. From this, it can be considered that the FCFs and the interference frequencies are the outliers of the frequency-domain sequence obtained by the discrete Fourier transform. The frequency detection indicator is defined as k = mean (y) + 3σ (27) where mean(y) and σ denote the mean and standard deviation of the frequency-domain sequence from [0, 4 f o ], respectively. Then, the MAR and the mode alignment error (MAE) are defined as where e i is the weight of the ith FCF MAR which is set to [0.5, 0.3, 0.2], n i is the number of the ith FCF detected in all channels, a is the number of channels, n is the total number of the detected FCFs, and m is the total number of the detected frequency. Therefore, the larger MAR indicates that more FCFs can be identified in the frequency domain, while the smaller MAE indicates fewer interfering frequencies exist in the frequency domain.
To visualize the indicators, the frequency detection threshold k is indicated by a red dotted line. The MAR of MDMD and standard DMD are 0.93 and 0.80, respectively. The MAE of MDMD and standard DMD are 0.11 and 0.63, respectively. Compared to the results of standard DMD, MDMD has better mode alignment properties and less interference frequencies. Benefiting from the circular convolution in the process of the tensor product, MDMD can fulfill the information fusion between channels, which enables MDMD to extract the common fault impulses in the multichannel signal.

B. Performance Analysis
In the above simulations, the mode alignment property of MDMD has been illustrated. Next, a comparative analysis with representative multivariate signal-processing methods (MEMD and MVMD) is presented below.
The performance of MDMD under strong noise is tested. The simulated signal with SNR = −10 dB is employed in this case. It can be seen from Fig. 5 that there are no obvious periodic impulses in the time-domain diagram, and corresponding FCFs in the spectrum are buried by strong noise. Fig. 6 is the envelope spectrum (ES) of the result obtained by MDMD. It can be seen from Fig. 6 that all FCFs in the three channels are extracted. Meanwhile, the noise mixed in the signal is almost removed due to the robust approach in MDMD. Fig. 7 shows the ES of the top six order intrinsic mode functions (IMFs) of MEMD. Although the FCFs can be identified from IMF2 of Channels 1 and 2, noise still exists in Channel 3, which leads to the MAR of MEMD only reaching 0.83. Meanwhile, a slight modal aliasing phenomenon appears in IMF2 and IMF3 of Channel 1, resulting in a lower amplitude of MEMD than the other two methods. Fig. 8 shows the ES of the result obtained by MVMD. Similar to the results of MEMD, the FCFs in Channel 3 are buried by noise. Hence, the simulation results verify that the proposed MDMD is superior to MEMD and MVMD under strong noise. There are two reasons for this: 1) the tensor operation in MDMD involves circular convolution between multichannel data, which realizes information fusion and synchronous extraction between channels and 2) the tensor low-rank approximation framework further enhances the noise robustness of MDMD.      fault diagnosis, the collected signal contains not only WGN, but also interference caused by other components. Therefore, two sets of actual multichannel signals are employed to further verify the effectiveness of the proposed method.

A. Bearing Inner-Race Fault Diagnosis
In this section, a multichannel fault signal from Case Western Reserve University (CWRU) is used to verify the effectiveness of MDMD. As shown in Fig. 9, the test rig consists of a motor, a torque sensor, and a dynamometer. The fault bearing is SKF 6205 that is located at the driving end of the motor. The speed of the test bearing is 1750 r/min, and the theoretical FCF of the inner-race is 158.20 Hz. The multichannel signals are collected by three sensors with a sampling frequency of 12 kHz, and the sensors' location is shown in Fig. 9(b). The waveform and the spectrum are shown in Fig. 10. As can be seen that the fault in the three channels shows obvious uneven energy distribution.
First, the proposed MDMD is employed to process the above multichannel signal. In this case, the size of snapshot matrices is 76 × 3 and the size of snapshot tensors is 76 × 150 × 3. The result of MDMD is shown in Fig. 11. It can be seen that all three channels clearly show the inner-race FCF and its multiples, and the MAR of MDMD is 0.9. Moreover, it is worth noting that although the fault impact still exhibits For comparison, two classic multichannel signal-processing methods (MEMD and MVMD [27]) are employed in this case. The result of MEMD is shown in Fig. 12. Although FCF can be identified from all three channels, there are still interference frequencies in the ES that affect fault diagnosis. The MAE of MEMD is 0.59, greater than 0.27 of MDMD. The results show that although MEMD can extract fault features synchronously, the noise robustness needs to be further improved. The ES of MVMD is shown in Fig. 13. The MAR of MVMD is only 0.5, far less than 0.9 of MDMD. The results show that although MVMD can extract fault features synchronously, its performance is greatly affected by the interference frequency in the case of uneven energy distribution. In this case, the runtime of MDMD is 1.4194 s, while the runtimes of MEMD and MVMD are 43.5722 and 1.0974 s, respectively. It should be noted that although the computational efficiency of MDMD is slightly lower than that of MVMD in this dataset when the parameter selection problem is ignored, the feature extraction performance of MDMD is much better than that of MVMD. The above results demonstrate that MDMD has satisfactory computational efficiency. The above algorithms are run on a desktop with a CPU of AMD Ryzen 2700, 3.20 GHz. The above comparison verifies the mode alignment property of MDMD, and MDMD outperforms MEMD and MVMD when multichannel signals have uneven energy distribution.

B. Bearing Outer-Race Fault Diagnosis
In this section, a set of multichannel signals with an outerrace fault is used to verify the noise robustness of the proposed approach. As shown in Fig. 14, the test rig is composed of an ac motor, a support shaft, a gearbox, and a magnetic powder brake. The fault bearing is SKF 6202 that is located at the support shaft [see Fig. 14(b)]. The speed of the test bearing is 900 r/min, and the theoretical FCF is 45.60 Hz. The sampling frequency is 12 800 Hz. The time-domain waveform and the   frequency spectrum are shown in Fig. 15. It can be seen from Fig. 15 that the FCFs are buried by strong noise.
First, MDMD is employed to process the above signal. In this case, the size of snapshot matrices is 288 × 3 and the size of snapshot tensors is 288 × 40 × 3. Fig. 16 is the result of MDMD. Obviously, MDMD can efficiently extract the FCFs in all channels. Meanwhile, it can be seen from Fig. 16 that MDMD almost eliminates noise and interference frequencies in the ES. Then, MEMD and MVMD are used to process the above signal, respectively. The result of MEMD and MVMD are shown in Figs. 17 and 18. The MAR of MDMD, MEMD, and MVMD are 1, 0.8, and 0.7, respectively. It can be seen from the results that MDMD is more noise robust than MEMD and MVMD due to the tensor low tubal rank component extraction framework of MDMD. In addition, the result verifies MAR indicator can accurately identify FCFs that are buried by noise. Moreover, MDMD has high computing efficiency. In this case, the runtime of MDMD is 0.7356 s, while the runtimes of MEMD and MVMD are 53.2031 and 7.7086 s, respectively. The comparison analysis verifies the noise robustness of MDMD, and its superiority   over other multivariate signal-processing algorithms in dealing with the multichannel signal containing strong noise.

VI. CONCLUSION
In this article, a novel multivariate signal-processing-based fault diagnosis approach is proposed, where DMD has been developed into multivariate scope. This article provides a new insight into multivariate signal processing, indicating an inspiration for the formulation of the proposed MDMD algorithm by tensor operations. The tensor low tubal rank component extraction framework is proposed to enhance the noise robustness of MDMD. Numerical simulation demonstrates the mode alignment property and noise robustness of MDMD. The subsequent experiments verify that MDMD can effectively extract multivariate fault features from multichannel bearing fault signals. The main novelties and contributions of the conducted research can be summarized as follows.
1) A novel multivariate extension algorithm of DMD which is MDMD is proposed in the article. This novel multivariate algorithm has inherent mode alignment property due to the tensor operation and is of great significance in fault feature extraction of multivariate signals. 2) A robust low tubal rank tensor component extraction framework is constructed in MDMD for the synchronous extraction of fault features hidden in the multivariate signal. In addition, the constructed framework enhances the noise robustness of the proposed MDMD algorithm. 3) An FCF detection indicator called MAR is proposed to measure the mode alignment property of multivariate signal-processing algorithms quantitatively. Both simulation and experiment verify the effectiveness of MAR in identifying FCF when the spectrum is interfered with by strong noises. 4) A novel multivariate signal-processing-based bearing fault diagnosis approach is proposed. The snapshot tensor is defined to convert the multivariate fault signal into the tensor format. The experimental results illustrate the superiority of MDMD over the traditional multivariate signal-processing algorithms in fault diagnosis of rolling bearing. It should be noted that compound fault diagnosis is very difficult in the field of fault diagnosis of rotary machinery, since MDMD achieves satisfactory results during fault feature extraction of multivariate signals of rolling bearing, it could offer a novel solution for compound fault diagnosis. Thus, in our future work, MDMD will be employed to diagnose the compound fault of rotating machinery. More tensor lowrank approximation frameworks of MDMD will be explored to enhance the performance of the proposed method.