Adaptive variable sampling T-SSD algorithm

Abstract

Addressing the challenges of non-unique decomposition outcomes and prolonged decomposition durations in the fault feature adaptive extraction algorithm based on tensor decomposition, this paper presents a novel algorithm called the adaptive variable sampling tensor singular spectrum decomposition (T-SSD) algorithm. The proposed approach centers on decomposing multichannel time series with adaptive sampling frequency, leveraging tensor singular value decomposition. Initially, the embedding dimension and the number of resampling points were optimized by power spectral density analysis and adaptive sampling algorithm. Subsequently, a third-order tensor is constructed based on the principle of tensor-tensor-preserving order multiplication, combining trajectory tensor construction and embedding dimension. Finally, the signal decomposition and reconstruction of multichannel component signals are achieved through adaptive sampling, interpolation, and complementary back steps. Experimental signal analysis indicates that this algorithm can better extract fault characteristics from signals compared to common signal processing algorithms. In comparison with the traditional T-SSD algorithm, this method significantly improves decomposition efficiency, particularly for low-frequency components. It effectively tackles the efficiency challenge caused by data redundancy, enabling the organic fusion and adaptive decomposition of multichannel signals.

Keywords

Tensor decomposition multichannel signals adaptive decomposition feature extraction third-order tensor

Introduction

Due to the environment in which mechanical equipment works, coupled with factors such as improper human operation equipment aging and other reasons, rotating machinery inevitably experiences failures. These failures often pose significant risks to property and personnel safety during operation. Bearings are the most basic components of rotating machinery. Therefore, the fault diagnosis of bearings has always been the top priority of scientific research. Currently, the field of fault diagnosis is mainly focused in two directions: one involves establishing a kinetic model of the equipment from the perspective of dynamics to analyze and research the failure mechanism; the second involves exploring signal processing algorithms for noise cancellation and fault feature extraction from vibration signals generated during the operation of the equipment. The existing signal processing algorithms mainly include empirical modal decomposition (EMD),¹ variational modal decomposition (VMD),² local mean decomposition³ and singular value decomposition (SVD). Matrix SVD is also suitable for analyzing noise-colored vibration signals. SVD is a form of principal component analysis.⁴ Compared with the sine and cosine of Fourier analysis, the advantage of the SVD is that the derived components are not necessarily harmonic functions, which makes it more suitable for analyzing non-stationary time series. By using SVD to decompose a matrix, singular values that represent the most basic characteristics of the matrix are obtained. The noise-containing signal trajectory matrix can be decomposed into multiple interpretable subspace components, comprising the useful fault characteristics components and useless interference components. The singular values of the components where the noise is located are discarded and only the useful signal is reconstructed, thus obtaining the purpose of removing the noise to obtain the fault characteristic frequency. The use of the SVD method has significant advantages when dealing with noise-containing signals.^5,6

Currently, numerous scholars are devoted to the research of SVD in fault diagnosis. For example, Li et al.⁷ proposed a correlated SVD method for bearing fault diagnosis that combines reconstruction component determination, optimal Hankel matrix structure selection, and envelope power spectrum analysis. In Golafshan and Sanliturk,⁸ an SVD and Hankel matrix-based de-noising algorithm is successfully used for ball bearing localized fault detection in both time and frequency domains. This approach helps eliminate background noise and improve the reliability of the fault detection. Qiao and Pan⁹ proposed an adaptive and singular value selection method based on correlation coefficient for fault feature detection, which achieves good results in detecting weak fault signals of rolling bearings. Wang and Xin¹⁰ proposed a weighted low-rank sparse model for bearing signal analysis that utilizes the low-rank physical structure of the periodic pulse characteristics. Zhang et al.¹¹ proposed a clustered low-rank method for fault diagnosis of engine bearing signals, which effectively eliminates the pathological phenomenon of overlapping of singular values and realizes the simultaneous detection of strong and weak features.

Although matrix SVD has certain advantages in dealing with noisy vibration signals, the algorithm still has certain limitations in dealing with multichannel signals. When a mechanical system fails, multiple parts and measurement points within the system may contain critical fault information. The use of multi-dimensional sensors to collect vibration signals from different locations of the system can effectively avoid the loss of local information.

Tensor as the most intuitive expression of higher-order data, not only can maintain the internal structure characteristics of the data itself, but also in the subsequent computational decomposition process closest to the essential properties of the data, so that the higher-order model will be more accurate in the data expression.¹² Vectors are first-order tensors and matrices are second-order tensors. If the signal is single-channel, the signal is usually expressed in the form of a vector and matrix, such as the time-domain waveform, spectrum, and time-frequency diagram of the signal. According to the mathematical significance of tensors, tensor decomposition is expected to overcome the limitations of matrix SVD when processing multichannel signals. Currently, a growing number of scholars have dedicated their efforts to the study of tensor decomposition, and it has been utilized in various fields.^13–19 The mainstream tensor decomposition methods in use are Canonical Polyadic (CP) decomposition,²⁰ Tucker decomposition,²¹ and higher-order SVD (HOSVD).^22–24 These methods are continually being optimized and refined. Since signals inevitably appear to be missing in the process of storing or transmitting, Zhou et al.²⁵ put forward an iterative HOSVD image missing data recovery method for efficient recovery of missing data in images. This method makes full use of the multiple constraints within the image and between the image and the image, and greatly improves the recovery accuracy. Combining the advantages of HOSVD and quaternion, Miao et al.²⁶ proposed a multi-focus color image fusion method and a color image denoising method. Karim et al.²⁷ proposed a new variant of the Levenberg-Marquardt algorithm used for Tensor CP decomposition with an emphasis on image compression and reconstruction. Li et al.²⁸ presented joint adaptive threshold HOSVD and the rearrangement of pixels in tensors method for image denoising, which exploits more prior knowledge to characterize the image features. Mao et al.²⁹ proposed a new self-supervised deep tensor domain-adversarial anomaly detection model, named Tensor-DAAD. Tensor-DAAD can find early fault occurrence in an earlier location with a much lower false alarm rate. The SVD of matrices has been extended to higher order in many ways,³⁰ but since the essence of the HOSVD algorithm is still to utilize the idea of matrices to deal with signals, there are problems of non-pseudo-diagonal core tensors and non-unique decomposition results in dealing with multichannel signals. To solve the above problems, Kilmer et al.³¹ defined a new tensor multiplication rule and the corresponding T-SVD. Based on the principle of tensor-preserving multiplication proposed by Kilmer, Huang and Cui³² developed a tensor singular spectrum decomposition (T-SSD) algorithm based on T-SVD, which can better extract the weak fault quantization features hidden in the original multichannel signal. Xu et al.³³ used the tensor core norm to represent the low-rank structural information of tensor data and developed a new low-rank tensor-based semi-supervised classifier, which can perform semi-supervised analysis of large multi-sensor signals when the number of labels is limited case, this method achieves high-accuracy fault classification. However, since tensors are usually used as high-order expressions, the length of data involved in T-SSD calculations is very sensitive to its computational efficiency, therefore, how to compress the computing time of T-SSD on the premise of retaining the advantage of tensor decomposition to extract weak fault features is a very challenging topic.

The core contribution of this paper is to propose the adaptive variable sampling T-SSD algorithm, which addresses the inefficiencies found in the traditional T-SSD method. To implement this innovation, the algorithm first evaluates the features of signals within the decomposed bandwidth, which informs the subsequent construction of a third-order tensor tailored to the specific characteristics of the signal. The signals are then resampled using an appropriate sampling frequency, and a tensor with suitable dimensions is built for decomposition. Subsequently, the component signals obtained from the iteration process are thoroughly assessed. The resampled signal components are reintegrated into the original signal’s components of equal length through interpolation. This step ensures that subsequent iterations are not affected by differences in array length. This approach fundamentally solves the problem of low computational efficiency (or even freezing) in traditional T-SSD, especially concerning the decomposition of low-frequency signal components.

The paper’s structure is outlined as follows: the section “Signal processing based on T-SVD” introduces tensor multiplication, the rule of T-SVD, and the traditional T-SSD. Section “The proposed method” presents a novel adaptive variable sampling T-SSD algorithm. This algorithm is developed based on T-SSD, adaptive sampling and the principles of interpolation complementation. It realizes resampling different data points when decomposing different frequency bands. Following components acquisition, the interpolation complementary back algorithm is executed. Section “Verification using experimental signals” the advantages and accuracy of the proposed method are verified by the experimental signal processing results. Section “Conclusion” provides a conclusion summarizing the key findings.

Signal processing based on T-SVD

Definition of notation and new tensor multiplication

$A \in R^{I_{1} \times I_{2} \times I_{3} \times \dots \times I_{N}}$ denotes an $N$ -dimensional tensor, which $I_{1}, I_{2}, \dots I_{N}$ correspond to the dimensions in each order, respectively. For a third-order tensor, it can be observed in different ways. Different features of the data are obtained by using different ways of observation. Figure 1 shows the slices of three different modes.

Figure 1.

Three different mode slices of a third-order tensor.

In tensor analysis, unfold and fold are frequently used operations of a tensor, where a third-order tensor is unfolded to be analyzed in a second-order tensor, and the second-order tensor is collapsed back to a third-order tensor after analysis to observe its features, for example, the positive expansion of the tensor $A \in R^{I_{1} \times I_{2} \times I_{3}}$ is I₃ frontal slice matrices of size I₁ × I₂ of A, (A(:, ; :,1),…, A(:, : ,I₃)), arranged up and down according to the order of the third dimension, it is spliced into a matrix of $I_{1} I_{3} \times I_{2}$ , as in Equation (1):

unfold (A, 1) = [\begin{matrix} A :; :; 1 \\ A :; :; 2 \\ \dots \\ A (:; :; I_{3}) \end{matrix}] \in R^{(I_{1} * I_{3}) \times I_{2}}

(1)

The second parameter $i$ in $unfold (A, i)$ is the specified model of the unfolding tensor, $i$ = 1 is the front unfolding model, $i$ = 2 is the side unfolding model, and $i$ = 3 is the horizontal unfolding model. For example, three unfolding schematics for a 3 × 3 × 3 tensor A are shown in Figure 2.

Figure 2.

Different mode expansions of a third-order tensor.

The multiplication of the different expansion modes of the tensor and the matrix is defined as the n-mode product. The n-mode product does not apply to the multiplication of two third-order tensors. Thus, a new multiplication rule of the tensor is proposed and its properties are described below.

Definition 1: (Tensor multiplication).³¹ It is known that the third-order tensor X is of size $I_{1} \times I_{2} \times I_{3}$ and the third-order tensor Y is of size $I_{2} \times l \times I_{3}$ . The product of X and Y is a tensor of size $I_{1} \times l \times I_{3}$ , the rule of multiplication is defined as follows:

X * Y = fold (circ (unfold (X, 1)) \cdot unfold (Y, 1)

(2)

For example, $X \in R^{n 1 \times n 2 \times 3}$ , $Y \in R^{n 2 \times l \times 3}$ , then by Equation (2):

X * Y = fold ([\begin{matrix} X_{1} & X_{3} & X_{2} \\ X_{2} & X_{1} & X_{3} \\ X_{3} & X_{2} & X_{1} \end{matrix}] [\begin{matrix} Y_{1} \\ Y_{2} \\ Y_{3} \end{matrix}], 1) \in R^{n 1 \times l \times 3}

The unfold and fold in Equation (2) denote the tensor unfolding and folding operations respectively, they are inverse operations to each other; $circ (\cdot)$ denotes the cyclic unfolding operation. If frontal faces of tensor X are $X_{1} = X (:, :, 1)$ ,…, $X_{I_{3}} = X (:, :, I_{3})$ , then:

circ (unfold (X, 1)) = [\begin{matrix} \begin{matrix} X_{1} & X_{I_{3}} \end{matrix} & \begin{matrix} X_{I_{3 - 1}} & \dots \end{matrix} & \begin{matrix} \dots & X_{2} \end{matrix} \\ \begin{matrix} X_{2} & X_{1} \\ \dots & \dots \end{matrix} & \begin{matrix} X_{I_{3}} & \dots \\ \dots & \dots \end{matrix} & \begin{matrix} \dots & X_{3} \\ \dots & \dots \end{matrix} \\ \begin{matrix} X_{I_{3}} & X_{I_{3 - 1}} \end{matrix} & \begin{matrix} \dots & \dots \end{matrix} & \begin{matrix} X_{2} & X_{1} \end{matrix} \end{matrix}]

(3)

From Kilmer et al.,³¹ it is known that the discrete Fourier transform (DFT) can diagonalize cyclic matrices, and similarly, block-cyclic matrices can be diagonalized by DFT. Assuming that the size of the tensor X is I₁ × I₂ × I₃ and $F_{I_{3}}$ is a $I_{3} \times I_{3}$ DFT matrix, then:

\begin{matrix} (F_{I_{3}}^{*} \otimes I_{I_{1}}) \cdot circ (unfold (X, 1)) \cdot (F_{I_{2}}^{*} \otimes I_{I_{2}}) \\ = [\begin{matrix} D_{1} & \begin{matrix}  \end{matrix} \\ D_{2} & \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} & \begin{matrix}  \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \end{matrix} & \begin{matrix} D_{I_{3}} \end{matrix} \end{matrix} \end{matrix}] \end{matrix}

(4)

where “⊗” refers to the Kronecker product, $F^{*}$ refers to the conjugate transpose of $F,$ and “·” refers to the standard matrix product.

Theorem: D_i is a face frontal of the tensor D, which is obtained by applying the fast fourier transform (FFT) to each fiber along, and a fiber is a column vector obtained by keeping $i, j$ fixed and taking $X (i, j, :)$ .

X * (Y * Z) = (X * Y) * Z

(5)

Definition 2 (Unit Tensor): The unit tensor $I_{nnl}$ is a tensor in which only the first frontal slice X , (;, : ,1) is a unit matrix of size $n \times n$ , and the rest of the elements are zero, that is:

Unfold (I_{nnl}, 1) = [\begin{matrix} I \\ 0 \\ 0 \\ \dots \\ 0 \end{matrix}]

(6)

Definition 3 (Inverse of a tensor): Just like the inverse of a matrix, if the inverse of a tensor $X$ of size $n \times n \times n$ is $Y$ , it is satisfied: X*Y = I and Y*X = I.

Definition 4 (Transpose of a tensor): If X is a tensor of $I_{1} \times I_{2} \times I_{3}$ , then $X^{T}$ is a tensor of $I_{2} \times I_{1} \times I_{3}$ , which is given by:

X^{T} = fold (unfold (X, 1, [1, I_{3} : - 1 : 2]), 1)

(7)

That is, $X^{T}$ is the front unfolding mode of $X$ and then the reversal of the order of the slices from the second to I₃ before folding.

Transposition of a tensor $X \in R^{n 1 \times n 2 \times n 3}$ can be obtained from Equation (7).

X^{T} = fold ([\begin{matrix} X_{1}^{T} \\ X_{n 3}^{T} \\ \begin{matrix} X_{n 3 - 1}^{T} \\ \begin{matrix} ⋮ \\ X_{2}^{T} \end{matrix} \end{matrix} \end{matrix}] 1)

(8)

Tensor SVD

If $X$ is a real tensor of size $I_{1} \times I_{2} \times I_{3}$ , then $X$ can be decomposed into the product of three tensors (as shown in Figure 3):

X = U * S * V^{T}

(9)

where $U$ is the orthogonal tensor of size $I_{1} \times I_{1} \times I_{3}$ , $V$ is the orthogonal tensor of size $I_{2} \times I_{2} \times I_{3}$ , $S$ is the f-diagonal tensor of size $I_{1} \times I_{2} \times I_{3}$ .

Figure 3.

Schematic diagram of third-order tensor T-SVD.

Table 1 shows the steps of T-SVD, and Figure 3 shows the corresponding decomposition mode.

Table 1.

T-SVD methods.

Algorithm1: T-SVD methods

Input: tensor

X ϵ R^{I_{1} \times I_{2} \times I_{3}}

of I₁×I₂×I₃.
Outputs:

U ϵ R^{I_{1} \times I_{1} \times I_{3}}

S ϵ R^{I_{1} \times I_{2} \times I_{3}}

and

V ϵ R^{I_{2} \times I_{2} \times I_{3}}

.
Process:D = fft(X,[ ],3);
for i = 1…I₃
[u,s,v] = svd(D( :,:,i ));

\tilde{U}

( :,:,i ) = u;

\tilde{V}

( :,:,i )=v;

\tilde{S}

( :,:,i )=s;
U=ifft(

\tilde{U},

[ ],3); V=ifft(

\tilde{V},

[ ],3); S=ifft(

\tilde{S},

[ ],3 );
end

T-SVD: tensor singular value decomposition.

Proof: SVD is applied to each matrix $D_{i}$ on the diagonal in Equation (4). That is, $D_{i} = U_{i} * S_{i} * V_{i}^{T}$ , then it has:

\begin{matrix} [\begin{matrix} D_{1} & \begin{matrix}  \end{matrix} \\ D_{2} & \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} & \begin{matrix}  \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \end{matrix} & \begin{matrix} D_{I_{3}} \end{matrix} \end{matrix} \end{matrix}] = [\begin{matrix} U_{1} & \begin{matrix}  \end{matrix} \\ U_{2} & \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} & \begin{matrix}  \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \end{matrix} & \begin{matrix} U_{I_{3}} \end{matrix} \end{matrix} \end{matrix}] \\ [\begin{matrix} S_{1} & \begin{matrix}  \end{matrix} \\ S_{2} & \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} & \begin{matrix}  \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \end{matrix} & \begin{matrix} S_{I_{3}} \end{matrix} \end{matrix} \end{matrix}] [\begin{matrix} V_{1}^{T} & \begin{matrix}  \end{matrix} \\ V_{2}^{T} & \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} & \begin{matrix}  \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \end{matrix} & \begin{matrix} V_{I_{3}}^{T} \end{matrix} \end{matrix} \end{matrix}] \end{matrix}

(10)

Simultaneous left multiplication by $(F_{I_{3}}^{*} \otimes I_{I_{1}})$ and right multiplication by $(F_{I_{3}} \otimes I_{I_{2}})$ on both sides of the Equation (10) yields:

\begin{array}{l} (F_{I_{3}}^{*} \otimes I_{I_{1}}) [\begin{matrix} D_{1} \\ D_{2} \\ ⋱ \\ D_{I_{3}} \end{matrix}] (F_{I_{3}} \otimes I_{I_{2}}) \\ = (F_{I_{3}}^{*} \otimes I_{I_{1}}) [\begin{matrix} U_{1} \\ U_{2} \\ ⋱ \\ U_{I_{3}} \end{matrix}] (F_{I_{3}} \otimes I_{I_{1}}) \\ (F_{n 3}^{*} \otimes I_{I_{1}}) [\begin{matrix} S_{1} & \begin{matrix}  \end{matrix} \\ S_{2} & \begin{matrix}  \end{matrix} \\ \begin{matrix}  \end{matrix} & \begin{matrix}  \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \end{matrix} & \begin{matrix} S_{I_{3}} \end{matrix} \end{matrix} \end{matrix}] \\ (F_{I_{3}} \otimes I_{I_{2}}) (F_{I_{3}}^{*} \otimes I_{I_{2}}) [\begin{matrix} V_{1}^{T} \\ V_{2}^{T} \\ ⋱ \\ V_{I_{3}}^{T} \end{matrix}] (F_{I_{3}} \otimes I_{I_{2}}) \end{array}

(11)

It follows from Equations (4) and (11):

\begin{matrix} circ (unfold (X), 1) = circ (unfold (U), 1) circ (unfold (S), 1) \\ circ (unfold (V^{T}), 1) \end{matrix}

(12)

i.e.: $X = U * S * V^{T}$

Construction of tensor models

A single-channel vibration signal acquired from a single sensor is represented as a vector. Therefore, to capture multichannel signals, it is essential to employ multiple sensors. Subsequently, these multichannel signals are combined and organized into a higher-order tensor, facilitating the representation of signals in higher dimensions. The T-SSD adopts the method of constructing the trajectory matrix from the matrix singular spectrum decomposition in Bonizzi et al.³⁴ and extends it to higher order. Taking two-channel signals as an example, assuming that the one-dimensional signals obtained by the two sensors are discrete time series, denoted as $x (n)$ and $y (n)$ , each of length $N$ . In this paper, x(n) and y(n) are constructed as a third-order trajectory tensor $X \in R^{I_{1} \times I_{2} \times I_{3}}$ , where n = 1, 2, 3, …, N; $I_{1} = N$ , $I_{3} = 2$ . The tensor $X \in R^{I_{1} \times I_{2} \times I_{3}}$ of the frontal slice is denoted as $X (:, :, i)$ , when $I_{2} = N$ , the third-order trajectory tensor is formulated as:

{\begin{matrix} X (:, :, 1) = [\begin{matrix} x (1) & x (N) & \begin{matrix} x (N - 1) & \dots & x (2) \end{matrix} \\ x (2) & x (1) & \begin{matrix} x (N) & \dots & x (3) \end{matrix} \\ \begin{matrix} ⋮ \\ x (N) \end{matrix} & \begin{matrix} ⋱ \\ x (N - 1) \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \\ ⋱ \end{matrix} & \begin{matrix} ⋱ \\ x (2) \end{matrix} & \begin{matrix} ⋱ \\ x (1) \end{matrix} \end{matrix} \end{matrix}] \\ X (:, :, 2) = [\begin{matrix} y (1) & y (N) & \begin{matrix} y (N - 1) & \dots & y (2) \end{matrix} \\ y (2) & y (1) & \begin{matrix} y (N) & \dots & y (3) \end{matrix} \\ \begin{matrix} ⋮ \\ y (N) \end{matrix} & \begin{matrix} ⋱ \\ y (N - 1) \end{matrix} & \begin{matrix} \begin{matrix} ⋱ \\ ⋱ \end{matrix} & \begin{matrix} ⋱ \\ y (2) \end{matrix} & \begin{matrix} ⋱ \\ y (1) \end{matrix} \end{matrix} \end{matrix}] \end{matrix}

(13)

The number of frontal slices coincides with the number of sensors. The expression of the Equation (13) is for the relation when $I_{2} = N$ . When upscaling multiple one-dimensional signals to higher-order tensors the size of the second dimension $I_{2}$ of each frontal slice should be updated according to the characteristics of the current iteration component, instead of a fixed value, which can reduce the decomposition error and achieve better decomposition effect. To realize the signal decomposition to achieve the establishment of the adaptive size of the higher-order tensor, the power spectral density (PSD) method is used here to determine the size of $I_{2}$ , and $I_{2} = kk \times f_{n} / f_{\max}$ is defined as the solution equation for $I_{2}$ , where, $f_{n}$ is the sampling frequency of the time series, and the two parameters $f_{\max}$ and $kk$ need to be solved by the following steps. Determine the maximum peak $f_{\max}$ and solve the PSD of the time series $x (n)$ and $y (n)$ to obtain the maximum peak frequencies of the signals of the two channels as $(f_{\max}, x)$ and $(f_{\max}, y)$ . Calculate the maximum peak frequencies of the two channels as the maximum peak of the whole trajectory tensor model $f_{\max}$ = $\max (f_{\max}, x f_{\max}, y)$ . Since only the spatial parameter $I_{2}$ in Equation (11) is a non-deterministic value, it is necessary to quantitatively analyze the effect of $I_{2}$ on the decomposition results and seek the optimal empirical value to prepare for the adaptive decomposition of multichannel signals.

A simulated signal is formed by integrating a harmonic signal (component 1) and an impulse signal (component 2).

X (t) = H \sin (2 π f_{c} t) e^{- pt} + \cos (2 π f_{a} t)

(14)

where $p$ =3000, $f_{a}$ = 3000, $f_{c}$ = 200, $H$ = 0.1. Figure 4 is the time-domain diagram of the simulated signal.

Figure 4.

Time-domain diagram of simulated signals.

The third-order tensor is constructed from the simulated signal, in which the $kk$ in $I_{2} = kk \times f_{n} / f_{\max}$ is increased from 0.2 to 2. Then, the third-order tensor is processed by the proposed method in this paper.

The components obtained by the proposed method are compared with those of the original signal. Correlation coefficients between the components obtained by the proposed method and those of the original signal are calculated for different $kk$ values, as shown in Figure 5.

Figure 5.

Correlation coefficient of decomposed components and the original signal at different $kk$ values: (a) Component 1 and (b) Component 2.

As can be seen in Figure 5, the correlation coefficients are bigger when the $kk$ values of component 1 are 0.8–1.6 and the $kk$ values of component 2 are 0.6–1.2. To obtain more general conclusions, the paper extends this analysis to include signals at various frequencies. Firstly, third-order tensors are established from simulated signals obtained by changing the frequency of harmonic components. Secondly, the correlation coefficients and the root mean square error (RMSE) between the decomposed components and the original simulation signal are calculated. Finally, the optimal range of $kk$ is determined according to the correlation coefficients and RMSE. The correlation coefficients and RMSE corresponding to the different $kk$ values are shown in Figures 6 and 7.

Figure 6.

Correlation coefficients of the components obtained by the adaptive variable sampling T-SSD algorithm and original signal’s components: (a) Component 1 and (b) Component 2.

Figure 7.

RMSE of the components obtained by the adaptive variable sampling T-SSD algorithm and the original signal’s components: (a) Component 1 and (b) Component 2.

From Figures 6 and 7, it can be obtained that the two components exhibit higher correlation coefficients and lower RMSE when kk ranges from 0.8 to 1.2. Therefore, kk is most effective within this range. In this study, kk is set to 1.2 in the algorithm.

Adaptive selection of multichannel component signals

After building the third-order tensor $X \in R^{I_{1} \times I_{2} \times I_{3}}$ based on the bandwidth of the current iterative decomposition, a T-SVD of $X$ is performed. Using the pseudo-diagonal properties of the core tensor obtained by T-SVD, the eigen component corresponding to each singular value can be obtained. The time series signal is reconstructed from the singular values and the corresponding eigen component is obtained from the decomposition. During the reconstruction, it can be noticed that the elements of each cross-diagonal in Equation (13) are the same. Therefore, this step simply averages along each diagonal element. When the reconstructed time series contains noise and other interfering signals, the eigen components corresponding to the singular values obtained from the T-SVD are not all useful signal components. Therefore, it becomes necessary to denoise the time series by “discarding” the corresponding singular values. This process aims to obtain clearer and more useful components free from noise and interferential components.

In determining the current iteration range, the matrix singular spectrum method proposed in Kilmer et al.³¹ is used, that is, an iterative method is used to extract the component signals $g^{(j)} (n)$ and $h^{(j)} (n)$ one by one sequentially. According to the proportion in the power spectrum of the component signals, decompose iteratively from large to small, which ensures that the frequency components of the current component signals are centered in the frequency band $[f_{\max} - Δ f, f_{\max} + Δ f]$ , where $Δ f$ denotes the half bandwidth of the main peak in the PSD of the time series. Thus, a subset $I_{j}$ is obtained based on all the eigen groups that have a prominent dominant frequency in the range of the spectrum $[f_{\max} - Δ f, f_{\max} + Δ f]$ and the one that predominantly contributes to the energy of the dominant peak among the selected modal component $(I_{j} = {i_{1}, \dots, i_{p}})$ . Its corresponding component signal is subsequently reconstructed along the cross-diagonal averaging method. The use of such an iterative approach avoids the identification of large singular value components corresponding to noisy and interfering signals in the signal, and the energy of the signal is used as an indicator to measure the components of the components, which obviously reduces the interference of random noise in the signal. As a result, signal components describing other frequency band components outside the range of $[f_{\max} - Δ f, f_{\max} + Δ f]$ are automatically discarded at the current iteration and will be efficiently identified in subsequent iterations. The purpose of this operation is to include in the frequency band selected in each iteration the frequency of the shock oscillation signal that describes the majority of the energy in the current time series. In addition, the width $Δ f$ of the main peak in the PSD is related to the oscillation time of the signal, since this quantity is unknown, it must be estimated from the PSD of the input signals $x (n)$ and $y (n)$ . Throughout the iterative process, to achieve automatic termination of the whole decomposition process, the termination condition uses the normalized mean square error (NMSE)³⁴ metric of the multichannel signal:

{\begin{matrix} {NMSE}_{x}^{(k)} = \frac{\sum_{i = 1}^{N} (v^{(k + 1)} {(i)}^{2})}{\sum_{i = 1}^{N} {(x (i))}^{2}} \\ {NMSE}_{y}^{(k)} = \frac{\sum_{i = 1}^{N} (w^{(k + 1)} {(i)}^{2})}{\sum_{i = 1}^{N} {(y (i))}^{2}} \\ NMS E^{(k)} = \max ({NMSE}_{x}^{(k)}, {NMSE}_{y}^{(k)}) \end{matrix}

(15)

where, $v^{(k + 1)} (n) = v^{(k)} (n) - {\tilde{g}}^{(k)} (n)$ and $w^{(k + 1)} (n) = w^{(k)} (n) - {\tilde{h}}^{(k)} (n)$ are the residual signals of the two-channel signals after extracting $k$ components. The end condition of this iterative decomposition is that the largest NMSE indicator value in the two-channel signals is less than a set threshold, otherwise, the residual signals will be further decomposed until the termination condition is satisfied, and the decomposition is finished to obtain all the components of the signals.

{\begin{matrix} x (n) = \sum_{k = 1}^{m} {\tilde{g}}^{(k)} (n) + v^{(m + 1)} (n) \\ y (n) = \sum_{k = 1}^{m} {\tilde{h}}^{(k)} (n) + w^{(m + 1)} (n) \end{matrix}

(16)

where, $m$ is the number of component signals obtained in the iterative process, $v^{(m + 1)} (n)$ and $w^{(m + 1)} (n)$ are the two-channel signals, respectively. It can be seen that the energy of the residual signal decreases in each iteration. The aforementioned steps are illustrated in Figure 8.

Figure 8.

Steps of signal decomposition based on T-SVD.

The proposed method

The algorithm described in Figure 8 operates with the frequency of the extracted components selected according to the PSD. The algorithm maintains fixed sampling frequency and sampling point numbers. If a higher sampling frequency persists during the extraction of components with lower frequencies, it will generate redundant data, thus inflating the dimensions for subsequent tensor construction. As a result, the decomposition time is prolonged and the decomposition efficiency is reduced. To solve the above deficiencies, this study proposes the adaptive variable sampling T-SSD algorithm. The optimal parameters are selected according to the center frequency of the selected component range to perform adaptive sampling and interpolation complementing, which solves the problem of low efficiency and even stuck caused by tensor data redundancy.

Adaptive down sampling

Based on the results of the comparison between the center frequency of the band range delineated in the current iteration and the set threshold, it is determined whether adaptive sampling is used in the current decomposition process, where the number of sampling points changes with the change of the center frequency of the component, and the number of sampling points $m$ and the downsampling coefficient $n$ are satisfied at the same time after adaptive sampling.

n = \frac{l}{m}

(17)

where $l$ is the total length of the original signal.

m = Ns * ff_1 (in_3)

(18)

The signal component length $l^{'}$ after adaptive sampling is:

l^{'} = \frac{l}{n} = \frac{Ns * l * ff_1 (in_3)}{f_{n}}

(19)

where $ff_1 (in_3)$ is the center frequency of the band range delineated in the current iteration, $f_{n}$ is the sampling frequency of the original signal, and $Ns$ is the ratio of the new sampling frequency to the center frequency of the current band. $Ns$ should be chosen to be as small as possible while ensuring that the signal is not distorted, thereby reducing the amount of computation. For this reason, a harmonic signal is constructed to determine the value of $Ns$ .

A simulated signal $y 1 = 2 \sin (50 t) + 5 \sin (10 t)$ is created with a sampling frequency of 25,600. Different multiples of component frequencies were used as sampling frequencies for adaptive sampling. Calculate the RMSE between the interpolated supplementary signal and the original signal. By calculating the RMSE, the smallest case is selected to determine the relevant parameter $Ns$ of adaptive sampling. The sampling frequency of the original signal is changed using MATLAB (Natick, Massachusetts, USA) toolbox resampling function. Using different sampling frequencies will lead to different data points collected during the same period. To adaptively adjust the data collection points, the signal is resampled accordingly. The RMSE between the interpolated supplementary signal and the original signal is presented in Table 2.

Table 2.

Multiplier $Ns$ and corresponding RMSE.

$Ns$	20	30	40	50	60
RMSE	0.0044	0.0036	0.0031	0.0027	0.0025

RMSE: root mean square error.

Based on the values of $Ns$ and RMSE, an interpolated fit is performed as shown in Figure 9. From Figure 9, it can be seen that before 50 times, the RMSE decreases significantly with the increase of $Ns$ , and the obvious inflection point of the image occurs at $Ns$ = 50. Taking into account the amount of computation at the later stage of decomposition of the signal and the size of the RMSE, $Ns$ is set to 50, so that the signal after adaptive sampling is closest to the original signal while the number of data points is minimal.

Figure 9.

The sampling frequencies are different multiples $Ns$ of the largest component frequency in the signal power spectrum and the corresponding RMSE.

The threshold $m$ of adaptive sampling is calculated by:

m = f_{n} / 50

(20)

where, $f_{n}$ is the original sampling frequency, $f_{n} =$ 51,200 Hz. Through the calculation of Equation (20), $m$ =1024, then $m$ is adjusted to 1000 to preserve the information contained in the original signal to a greater extent and minimize the loss of information.

When the value of $Ns$ is determined to be 50, and then the downsampling coefficient $n$ in adaptive sampling is calculated by:

n = \frac{Ns * f f_{1 (i n_{3})}}{f_{n}}

(21)

Table 3 gives the steps of adaptive sampling.

Table 3.

Adaptive sampling.

Algorithm2: Adaptive sampling

Input: A signal of length l and sampling frequency

f_{n}

.
Output: A signal of length

\frac{Ns * l * ff_1 (in_3)}{f_{n}}

.
Process: if ff_1(in3_1)<=m
v2_1=downsample(v2_1, factor);
v2_2=downsample(v2_2, factor);
end

Interpolation complement

After adaptive sampling, the iterative decomposition often yields component signals length of shorter length compared to the original signal. This discrepancy can impact the calculation of the residual (residual signal = original signal – iterative component). Hence, when the decomposition and extraction of component signals in adaptive sampling, it becomes necessary to complementarily interpolate back to the length of the original signal before proceeding with subsequent iterative loops. The interpolation complement algorithm in this paper needs to first judge whether the length of the obtained component matches that of the original signal. This assessment helps determine if the current component has undergone adaptive sampling. Subsequently, the algorithm decides whether the current component needs interpolation complementation. In MATLAB, the main interpolation methods for one-dimensional arrays include $linear$ , $pchip$ , $spline$ , and $nearest$ . To select a specific interpolation complement method, we also perform adaptive sampling on simulated signal $y 2 = 2 \sin (500 x) + 5 \sin (10 x) + 50 \sin (8 x) + 100 \sin (500 x)$ . After adaptive sampling of the signal $y 2$ , different interpolation methods are used to fill the data points to be consistent with the original signal. The specific interpolation method is determined by calculating the RMSE between the interpolated supplementary signal and the original signal. Table 4 shows the RMSE values for the different interpolation methods along with their corresponding time consuming. From Table 4, it can be seen that the spline interpolation yields the minimum RMSE. Therefore, the algorithm in this paper adopts spline interpolation as the interpolation method. Table 5 gives the flow of the interpolation complement. Figure 10 shows the main steps of adaptive sampling.

Table 4.

RMSE of interpolation methods and corresponding time consuming.

Methodologies	RMSE	Time consuming
Linear	2.5 × 10⁻³	3.61 s
Pchip	2.18 × 10⁻⁴	3.67 s
Spline	1.74 × 10⁻⁸	3.75 s
Nearest	4.1 × 10⁻¹	3.80 s

RMSE: root mean square error.

Table 5.

Interpolation complement.

Algorithm 3: Interpolation complement

Input: A component of length

(Ns * l * ff_1 (in_3)) / f_{n}

after adaptive sampling.
Output: A component with the same length as the original signal after interpolation complement.
Process ls1=length(RR1), ls2=length(RR2);
if ls1<l|ls2<l
czx=t;
czy1=RR1(k1,:);
czy2=RR2(k1,:);
czx=downsample(t,factor);
clear RR1 RR2;
RR1(k1,:)=interp1( ls1,czy1,czx," interpolation method ");
RR2(k1,:)=interp1( ls2,czy2,czx," interpolation method ");
end

Figure 10.

Main steps of the adaptive sampling algorithm.

Adaptive variable sampling T-SSD algorithm

This section proposes an adaptive variable sampling T-SSD algorithm combining adaptive sampling and interpolated supplementary algorithm. The flowchart of the proposed method is shown in Figure 11, and its main steps are as follows.

Step 1: Determine the range of signals to be decomposed based on the PSD.

Step 2: According to the signal characteristics in the current iteration range, the adaptive downsampling step is carried out to reduce unnecessary redundant data.

Step 3: The higher-order tensors with appropriate dimension size are established and T-SVD decomposition is carried out. Select the corresponding singular values in the range to reconstruct the signal and get the component signals.

Step 4: By comparing the component length obtained by decomposition with the length of the original signal, the implementation of the interpolated supplementary algorithm is selected adaptively to ensure the smooth progress of the subsequent iteration.

Step 5: The decomposition reaches the termination condition, obtains a series of components, and the decomposition is complete.

Figure 11.

Flowchart of the adaptive variable sampling T-SSD algorithm.

Verification using experimental signals

To verify the validity of the adaptive variable sampling T-SSD algorithm, experimental data from a gearbox are used. A fault is located on a bearing on one of the gearbox’s three output shafts. Six acceleration sensors are installed on the experimental gearbox, and the relative positions of these sensors are shown in Figure 12.

Figure 12.

The relative position of the sensor on the experimental gearbox.

Figure 13 shows the faulty bearing with a 1 mm hole machined on the surface of the outer ring of the bearing using electrical impulses to simulate the point of failure.

Figure 13.

Bearing with outer ring failure.

Case 1

The following analysis of the experimental signal is conducted to evaluate the performance of the proposed method of this paper in extracting the fault characteristic frequency, the gearbox speed is set to 1683 rpm. The location of the faulty bearing is shown in Figure 12. Vibration acceleration sensor 1 and sensor 2, corresponding to the horizontal and vertical directions of the bearing respectively, are used with a sampling frequency of $f_{n}$ = 51,200 Hz to acquire vibration acceleration signals of a length of 100,000. The bearing failure characteristic frequency is 434.688 Hz. In order to increase the difficulty of extracting the characteristic frequency, several random interference signals are added to the experimental signal. Figure 14 shows the time-domain diagram, spectrum diagram and envelope spectrum of the experimental signal.

Figure 14.

Time-domain waveform and frequency and envelope spectra of the experimental signal.

The experimental signals are iteratively decomposed by both the traditional T-SSD and the adaptive variable sampling T-SSD algorithm to obtain the component signals, as shown in Figure 15. It can be seen from Figure 15 that the fault characteristic frequency and its harmonics can be found in component 2. This indicates that the proposed method has the same feature extraction capabilities as the traditional T-SSD. However, in terms of decomposition efficiency, the traditional T-SSD algorithm requires 842.46 s, while the method proposed in this paper takes only 171.05 s, achieving an efficiency improvement of nearly five times over the traditional T-SSD.

Figure 15.

The components containing fault characteristic frequency of two algorithms: (a) Component 2 envelope spectrum obtained by traditional T-SSD and (b) Component 2 envelope spectrum obtained of adaptive variable sampling T-SSD algorithm.

Case 2

New fault vibration signals at a speed of 400 rpm are then collected from the experimental gearbox. The time-domain waveform and envelope spectrum of the experimental signal are shown in Figures 16 and 17. At this speed, the fault characteristic frequency is 98.56 Hz. Since the proposed method is a type of signal processing, in order to demonstrate its ability to explore weak fault characteristic frequency in signals, it is compared with traditional signal processing methods (EMD, VMD) using the same set of signals. EMD³⁵ is firstly used to decompose the vibration signal into multiple intrinsic mode functions (IMFs), and the components strongly correlated with the original signal are selected to reconstruct the signal. The time-domain diagram and envelope diagram of the reconstructed signal are shown in Figures 18 and 19.

Figure 16.

Time-domain waveform of experimental signal.

Figure 17.

Envelope spectrum of experimental signal.

Figure 18.

Time-domain waveform after EMD and reconstruction.

Figure 19.

Envelope spectrum after EMD and reconstruction.

In the VMD² algorithm, the number of decomposition levels k, and regularization parameter α greatly influence the decomposition results. Thus, the whale optimization algorithm is used to optimize the parameters k to obtain the optimal decomposition effect. The result shows that when k is 7 and α is 2000, the result of VMD is the optimal solution. The IMFs whose kurtosis value is greater than four obtained by VMD of the original signal are reconstructed. The time-domain waveform and envelope spectrum of the reconstructed signal are shown in Figures 20 and 21.

Figure 20.

Signal time-domain waveform after VMD and reconstruction.

Figure 21.

Signal envelope spectrum after VMD and reconstruction.

It is evident from Figures 18 to 21 that when EMD and VMD algorithms are used to process the early fault vibration signals, the fault characteristic frequency is drowned by overwhelmed by the interference components.

Finally, the adaptive variable sampling T-SSD algorithm is used to process the same signals, and the results are shown in Figures 22 and 23. Component 2 extracted by the adaptive variable sampling T-SSD algorithm contains the fault characteristic frequency and its harmonics, which are relatively more obvious. This indicates that the adaptive variable sampling T-SSD algorithm is more effective at extracting early faults than the EMD and VMD algorithms.

Figure 22.

Time-domain waveform of component 2 obtained after decomposition by adaptive variable sampling T-SSD algorithm

Figure 23.

Envelope spectrum of component 2 obtained after decomposition by adaptive variable sampling T-SSD algorithm.

Efficiency analysis

The above findings underscore the significant enhancement in decomposition effect and efficiency achieved by the method proposed in this paper. To fully prove the superiority of the proposed method in terms of efficiency, three sets of vibration signals were processed using the proposed method and traditional T-SSD. The time consuming of the two decomposition algorithms and the improved efficiency of the adaptive variable sampling T-SSD algorithm are shown in Table 6. The formula used to calculate the efficiency improvement is: Efficiency improvement (%) = 100 × (time of traditional T-SSD − time of proposed method)/time of traditional T-SSD. It can be seen from Table 6 that when decomposing the same data, the adaptive variable sampling T-SSD algorithm has a faster decomposition speed. The decomposition efficiency was improved by 78%, 85.8%, and 95.7%, with an average improvement of 86%.

Table 6.

Comparison of time consumption of traditional T-SSD and the proposed method.

Case	Traditional T-SSD	The proposed method	Efficiency improvement	Average efficiency improvement
1	525.36 s	115.71 s	78%	86%
2	245.50 s	35.50 s	85.5%
3	1231.25 s	52.50 s	95.7%

T-SSD: tensor singular spectrum decomposition.

Error analysis

As can be seen from the section “efficiency analysis,” the proposed method in this paper has obvious advantages in decomposition efficiency. To further assess the accuracy of decomposition results, the RMSE values between components obtained from the traditional T-SSD algorithm and the adaptive variable sampling T-SSD algorithm are analyzed using the same data. These RMSE values of the first 10 components from the traditional T-SSD algorithm and the adaptive variable sampling T-SSD algorithm are presented in Table 7. From the results shown in Table 7, it is evident that the RMSEs of the first 10 components are very small, which proves the accuracy of the adaptive variable sampling T-SSD algorithm.

Table 7.

RMSE of the first 10 components obtained by two algorithms.

Component	RMSE
1	0
2	0
3	0.0405
4	0.0848
5	0.0768
6	0.0519
7	0.0687
8	0.0616
9	0.0624
10	0.0505

RMSE: root mean square error.

Conclusion

An adaptive variable sampling T-SSD algorithm was proposed in this paper to effectively tackle the problem of low decomposition efficiency and potential stagnation caused by data redundancy during signal processing in traditional T-SSD. Several conclusions can be drawn as follows.

The analysis results of experimental data demonstrate that the proposed adaptive variable sampling T-SSD algorithm can greatly improve signal decomposition efficiency. Specifically, when applied to the same set of experimental signals, the proposed algorithm enhances the signal decomposition efficiency by an average of 86% relative to the traditional T-SSD algorithm.

The analysis results of experimental data demonstrate that the proposed algorithm excels at extracting fault features. When analyzing early fault signals from bearings, the fault characteristic frequency and its harmonics can be clearly seen in the component envelope spectrum obtained by the adaptive variable sampling T-SSD algorithm. Furthermore, comparison analysis with traditional signal processing algorithms. (EMD, VMD) confirms that the proposed method provides superior performance in extracting weak fault features.

Footnotes

Appendix

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Natural Science Foundation of China under Grant Nos. 52105109, 52305115, and 52161135101. Natural Science Foundation of Inner Mongolia No. 2023QN05031, and the China Postdoctoral Science Foundation No. 2023M741938.

ORCID iDs

Feibin Zhang

Jinfeng Huang

References

Dou

. Fault feature extraction for gearboxes using empirical mode decomposition. Adv Mater Res 2012; 383–390: 1376–1380.

Dragomiretskiy

Zosso

. Variational mode decomposition. IEEE Trans Signal Process 2014; 62(3): 531–544.

Chen

Zhang

Yang

. Local mean decomposition method and its application to roller bearing fault diagnosis. China Mech Eng 2009; 20: 2711–2717.

Muruganatham

Sanjith

Kumar

, et al. Inner race bearing fault on detection using singular spectrum analysis. In: International conference communication control and computing technologies, Ramanathapuram, India, 2010.

Zhang

Zhao

Xiao

. Fault diagnosis of gear based on singular value decomposition and RBF neural network. In: 2017 2nd International conference on frontiers of sensors technologies (ICFST), China, 2017. IEEE.

Luo

Zhang

. Rolling bearing incipient fault detection based on a multi-resolution singular value decomposition. Appl Sci 2019; 9(20): 4465.

Liu

, et al. Correlated SVD and its application in bearing fault diagnosis. IEEE Trans Neural Netw Learn Syst 2021; 34: 355–365.

Golafshan

Sanliturk

. SVD and Hankel matrix based de-noising approach for ball bearing fault detection and its assessment using artificial faults. Mech Syst Signal Process 2016; 70: 36–50.

Qiao

Pan

. SVD principle analysis and fault diagnosis for bearings based on the correlation coefficient. Meas Sci Technol 2015; 26(8): 085014.

10.

Wang

Liao

Ding

, et al. Periodical sparse low-rank matrix estimation algorithm for fault detection of rolling bearings. ISA Trans 2020; 101(2): 366–378.

11.

Zhang

Chen

, et al. A clustering low-rank approach for aero-enging bearing fault detection. In: 2019 IEEE international instrumentation and measurement technology conference (I2MTC), 2019. IEEE.

12.

Yang

. Research and application of machine learning algorithm based tensor representation. Beijing: China Agricultural University, 2014.

13.

Huang

Zhang

Safaei

, et al. The flexible tensor singular value decomposition and its applications in multisensor signal fusion processing. Mech Syst Signal Process 2024; 220: 111662.

14.

Vasilescu

Terzopoulos

. Multilinear image analysis for face recognition. In: Proceedings of the international conference on pattern recognition (ICPR 2002), 2, Quebec City, Canada, Aug, 2002, pp. 511–514.

15.

Vasilescu

Terzopoulos

. Multilinear subspace analysis of image ensembles. In: Proceedings of the 2003 IEEE computer society conference on computer vision and pattern recognition (CVPR 2003), 2003, vol. 2, p. 128.

16.

Berkant

Eldén

. Handwritten digit classification using higher order singular value decomposition. Pattern Recognit 2007; 40(3): 993–1003.

17.

Kolda

Bader

. Higher-order web link analysis using multilinear algebra. In: Proceedings of the 5th IEEE international conference on data mining (ICDM 2005), Houston, TX, IEEE Computer Society (2005), pp. 242–249.

18.

Das

Kyal

. Efficient multichannel EEG compression by optimal tensor truncation. Biomed Signal Process Control 2021; 68: 102749.

19.

Wang

. Multidimensional denoising of rolling element bearings with compound fault based on tensor factorization. In: International conference on machinery, materials and computer Atlantis press, US, 2017, pp. 350–356.

20.

Carroll

Chang

. Analysis of individual differences in multidimensional scaling via an n-way generalization of “Eckart-Young” decomposition. Psychometrika 1970; 35(3): 283–319.

21.

Tucker

. Some mathematical notes on three-mode factor analysis. Psychometrika 1966; 31(3): 279–311.

22.

Lathauwer

. Blind source separation by higher-order singular value decomposition. Signal Process VII, Theor Appl 1994; 1: 175–178.

23.

Lieven

Bart

Joos

. A multilinear singular value decomposition. SIAM J Matrix Anal Appl 2000; 21(4): 1253–1278.

24.

Omberg

Golub

Alter

. A tensor higher-order singular value decomposition for integrative analysis of DNA microarray data from different studies. Proc Natl Acad Sci 2007; 104(47): 18371–18376.

25.

Zhou

Qiu

Liu

, et al. Iterative tensor high-order singular value decomposition method for image restoration. Appl Res Comput 2013, 30(11): 3488–3491.

26.

Miao

Kou

Cheng

, et al. Quaternion higher-order singular value decomposition and its applications in color image processing. Inform Fusion 2023; 92: 139–153.

27.

Karim

Dulal

Navasca

. Modified Levenberg-Marquardt algorithm for tensor CP decomposition in image compression. arXiv preprint arXiv:2401.04670, 2024.

28.

Pan

, et al. Adaptive thresholding HOSVD with rearrangement of tensors for image denoising. Multimedia Tools Appl 2020; 79: 19575–19593.

29.

Mao

Chen

Zhang

, et al. Tensor-DAAD: when tensor meets online early fault detection with transfer learning. Measurement 2023; 208: 112478.

30.

Kolda

. Tensor decompositions and applications. Siam Review 2009; 51(3): 455–500.

31.

Kilmer

Martin

Perrone

. A third-order generalization of the matrix SVD as a product of third-order tensors. Tech Report TR-2008-4, Tufts University, Department of Computer Science, october 2008.

32.

Huang

Cui

. Tensor singular spectrum decomposition: multi-sensor denoising algorithm and application. IEEE Trans Instrum Meas 2023; 72: 1–15.

33.

Wang

Huang

, et al. Semi-supervised multi-sensor information fusion tailored graph embedded low-rank tensor learning machine under extremely low labeled rate. Inf Fusion 2024; 105: 102222.

34.

Bonizzi

Karel

JMH

Meste

, et al. Singular spectrum decomposition: a new method for time series decomposition. Adv Adaptive Data Anal 2014; 6: 107–109.

35.

Yang

. Adaptive compound fault diagnosis of rolling bearings. J Vibr Shock 2016; 35(10): 145–150.