Singular Spectral Mode Decomposition Method for Separating Unsteady Components in Spectral Aliasing Environments

Abstract

This invention discloses a singular spectral mode decomposition method for separating non-steady-state components in a spectral aliasing environment, comprising the following steps: reconstructing all elements in the original signal x into a high-rank Hankel matrix X; performing singular value decomposition on the high-rank Hankel matrix X to obtain a left matrix, a singular value matrix, and a right matrix; and performing fast singular spectral decomposition on the left matrix, the singular value matrix, and the right matrix to separate the m basic components y. i Divide into groups, and group y into the same group i Summing yields the eigenmode function components of the original signal x. In the singular spectrum mode decomposition method of this invention, fast singular spectrum decomposition utilizes fast Fourier transform to replace a large number of matrix calculations in singular spectrum analysis, greatly shortening the calculation time. While ensuring calculation speed, it can separate non-steady-state eigenmode function components that overlap in the frequency domain, and can separate non-steady-state components in speech and mechanical vibration signals.

Description

Technical Field

[0001] This invention belongs to the field of nonsteady-state signal processing technology, specifically relating to a Singular Spectrum Mode Ensemble (SME) method for separating nonsteady-state components in a spectral aliasing environment. Background Technology

[0002] Non-steady-state components exist in speech and mechanical vibration signals, and the characteristics of these components are crucial for speech quality analysis and mechanical fault diagnosis. These non-steady-state features cannot be extracted using traditional short-time analysis methods such as the short-time Fourier transform (STFT) (Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.). Mode decomposition methods include Empirical Mode Decomposition (EMD) (NE Huang, et al., “The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences, vol. 454, no. 1971, pp. 903–995, 1998), Empirical Wavelet Transform (EWT) (J. Gilles, “Empirical wavelet transform,” IEEE transactions on signal processing, vol. 61, no. 16, pp. 3999–4010, 2013.), and Variational Mode Decomposition (VMD) (K. Dragomiretskiy and D. Zosso, “Variational mode decomposition,” IEEE transactions on signal processing). (Processing, Vol. 62, No. 3, pp. 531–544, 2013.) Due to limitations in decomposition capabilities, it is impossible to fully separate overlapping components on the spectrum. Therefore, it is often impossible to separate non-steady-state components when analyzing real speech and mechanical vibration signals.

[0003] Singular Spectrum Analysis (SSA) (R. Vautard and M. Ghil, “Singular spectrum analysis in nonlinear dynamics, with applications topaleoclimatic time series,” Physica D: Nonlinear Phenomena, vol. 35, no. 3, pp. 395–424, 1989) is based on Singular Value Decomposition (SVD) and Hankel matrices to separate trends, amplitude modulation, and noise from time series data. Studies have shown that high-rank Hankel matrices are crucial for separating overlapping components in signals because they allow for higher frequency domain resolution. However, high-rank Hankel matrices lead to a huge computational burden for SSA and can cause over-decomposition. A key challenge is how to maintain computational speed while using high-rank Hankel matrices and adaptively combine over-decomposed components. Summary of the Invention

[0004] To address the shortcomings of existing technologies, the purpose of this invention is to provide a high-speed, accurate, and fully adaptive Singular Spectrum Mode Ensemble (SME) method for separating non-steady-state components in spectral aliasing environments.

[0005] The present invention is achieved through the following technical solution.

[0006] A singular spectral mode decomposition method for separating unsteady components in a spectral aliasing environment includes the following steps:

[0007] S1, reconstruct all elements of the original signal x into a high-rank Hankel matrix, obtaining the high-rank Hankel matrix X:

[0008]

[0009] Where, x k Let be the k-th element in the original signal x, where k = 1, 2, ..., N, N is the number of sampling points in the original signal x, and m is the rank of the high-rank Hankel matrix. n = N - m + 1;

[0010] The original signal x is a one-dimensional time series data obtained by sampling the vibration signal. The original signal x is defined as a signal composed of multiple intrinsic mode functions.

[0011] In S1, the vibration signal is either a voice signal or a mechanical vibration signal.

[0012] In S1, the functional expression x(t) of the original signal x is:

[0013]

[0014] Where, r j (t) is the functional expression of the j-th intrinsic mode function, j = 1, 2, ..., J, where J is the number of intrinsic mode functions contained in the original signal x.

[0015] In the above technical solutions, Among them, a j (t) is the amplitude modulation function in the j-th intrinsic mode function, f j (τ) is the frequency modulation function in the j-th intrinsic mode function.

[0016] S2, Perform singular value decomposition on the high-rank Hankel matrix X to obtain a left matrix, a singular value matrix, and a right matrix, where the left matrix is ​​matrix [U]. m×m The singular value matrix is ​​matrix [S]. m×n The right matrix is ​​matrix [V] T ] n×n ;

[0017] S3, for the left matrix, singular value matrix, and right matrix, performs fast singular spectral decomposition, including: S3-1 to S3-4.

[0018] S3-1, based on the left and right matrices, calculate m vectors U. i (k) and m vectors V i (k), each vector U i (k) is represented as follows:

[0019] U i (k)=FFT(U ：，i )

[0020] Each vector V i (k) is represented as follows:

[0021] V i (k)=FFT(V T i，： )

[0022] Where i = 1, 2, ..., m, FFT(U ：，i ) indicates selecting the left matrix [U] m×m All elements in the i-th column are computed using the Fast Fourier Transform (FFT(V)). T i,： ) represents the right matrix [VT ] n×n All elements in the i-th row are computed using the Fast Fourier Transform;

[0023] S3-2, in each of the vectors U i (k) and each of the vectors V i The last position of (k) is padded with 0s for each of its elements, so that vector U i (k) and vector V i The length of (k) is N, and vectors are obtained sequentially. sum vector

[0024] S3-3, the m vectors obtained from S3-2 and m vectors Calculate and obtain m Y values i (k), each Y i The formula for calculating (k) is:

[0025] S3-4, for each Y i (k) is used to calculate and obtain the basic component y. i :

[0026] y i =IFFT(Y i (k))⊙(σ i L)

[0027] Where IFFT is the inverse fast Fourier transform calculation, σ i Represents the singular value matrix [S] m×n Let L be the singular value in the i-th row and i-th column of a vector, ⊙ be the symbol for vector multiplication, and L be a vector with N elements, where L is the l-th element of L. l for:

[0028]

[0029] S4, which divides m basic components y i Divided into Grouping y into the same group i Summing, we get The intrinsic mode function components of an original signal x, where m fundamental components y i Divided into The specific steps for grouping are as follows:

[0030] S4-1, Set the discrimination parameter γ∈(0,1), and assign y1 to y Δ Create the first group, set y1 as the first element of the first group, and set g = 1;

[0031] S4-2, Calculate y Δ and y g+1 99% of the bandwidth is used, per group y Δ and y g+1 After calculating 99% of the bandwidth used, two bandwidth usage ranges are obtained: (s) Δ e Δ ) and (s g+1 e g+1 );

[0032] S4-3, Calculate (s Δ e Δ ) and (s g+1 e g+1 The intersection of ) yields (s Δ，g+1 e Δ，g+1 );

[0033] (s Δ，g+1 e Δ，g+1 )=(s Δ e Δ )∩(s g+1 e g+1 )

[0034] S4-4, Calculate the bandwidth intersection ratio (P) OBW ) g :

[0035]

[0036] S4-5, (P OBW ) g Compare with the discriminant parameter γ:

[0037] When (P) OBW ) g When <γ, generate a new group, and y g+1 As the first element of the new grouping and y g+1 Assign a value to y Δ ;

[0038] When (P) OBW ) g ≥γ, y g+1 Add to y g In the group, compare y Δ and y g+1 99% of the bandwidth is used, so assign a value to y that has a wider bandwidth. Δ ;

[0039] S4-6, When g is less than m-1, increment the value of g by one, and then use the value of y obtained in S4-5. Δ Substitute into S4-2 and execute S4-2 to S4-6; when g equals m-1, the calculation ends.

[0040] The features and beneficial effects of this invention are as follows:

[0041] (1) In the singular spectrum mode decomposition method of the present invention, the fast singular spectrum decomposition (FSSD) uses the fast Fourier transform to replace a large number of matrix calculations in singular spectrum analysis, which greatly shortens the calculation time.

[0042] (2) This invention proposes an adaptive grouping method based on the intersection ratio of singular values ​​and bandwidth, which adaptively groups a large number of over-decomposed elementary component matrices, thus solving the problem of over-decomposition.

[0043] (3) The Singular Spectrum Mode Synthesis (SME) algorithm proposed in this invention can separate non-steady intrinsic mode function (IMF) components that overlap in the frequency domain while ensuring computational speed, and can separate non-steady components in speech and mechanical vibration signals. Attached Figure Description

[0044] Figure 1 The instantaneous time spectrum of the true component of analog signal 1 (left) and the instantaneous time spectrum calculated using the singular spectral mode decomposition method (right);

[0045] Figure 2 The instantaneous time spectrum of the true component of the analog signal (left) and the instantaneous time spectrum calculated using the singular spectral mode decomposition method (right);

[0046] Figure 3 The spectrogram of speech signal 1 (top) and the instantaneous time-frequency spectrum obtained by singular spectral mode decomposition (bottom) are shown.

[0047] Figure 4 The image shows the acoustic spectrum of the rolling bearing fault signal (top) and the instantaneous time spectrum calculated using the singular spectrum mode decomposition method (bottom). Detailed Implementation

[0048] The following detailed description, with reference to the accompanying drawings, describes the singular spectral mode decomposition method of the present invention for separating unsteady components in a spectral aliasing environment.

[0049] Example 1

[0050] A singular spectral mode decomposition method for separating unsteady components in a spectral aliasing environment includes the following steps:

[0051] S1, reconstruct all elements of the original signal x into a high-rank Hankel matrix, obtaining the high-rank Hankel matrix X:

[0052]

[0053] Where, x k Let be the k-th element in the original signal x, where k = 1, 2, ..., N, N is the number of sampling points in the original signal x, and m is the rank of the high-rank Hankel matrix. To maximize m, n = N - m + 1; in this embodiment, N = 2049. n = N - m + 1 = 1026.

[0054] The original signal x is a one-dimensional time series data obtained by sampling a vibration signal. The vibration signal can be a speech signal or a mechanical vibration signal. In this embodiment, the one-dimensional time series data uses an analog signal, a speech signal, or a mechanical vibration signal. The original signal x is set to be a signal composed of J intrinsic mode functions (IMFs), and the functional expression x(t) of the original signal x is:

[0055]

[0056] Where, r j (t) is the functional expression of the j-th intrinsic mode function (IMF), j = 1, 2, ..., J, where J is the number of intrinsic mode functions contained in the original signal x, and t represents time;

[0057] r j (t) is:

[0058] r j (t)=a j (t)cos(2π∫0 t f j (τ)dτ)

[0059] Among them, a j (t) is the amplitude modulation function in the j-th intrinsic mode function (IMF), f j (τ) is the frequency modulation function in the j-th intrinsic mode function (IMF);

[0060] S2, Perform singular value decomposition (SVD) on the high-rank Hankel matrix X to obtain the left matrix, the singular value matrix, and the right matrix, where the left matrix is ​​matrix [U]. m×m The singular value matrix is ​​matrix [S]. m×n The right matrix is ​​matrix [V] T ] n×n ;

[0061] S3. For the left matrix, singular value matrix, and right matrix, perform Fast Singular Spectrum Decomposition (FSSD). The specific steps are as follows:

[0062] S3-1, based on the left and right matrices, calculate m vectors U. i (k) and m vectors V i (k), where k represents the phase and vector U i (k) and vector V i (k) is represented as follows:

[0063] U i (k)=FFT(U ：，i )

[0064] V i (k)=FFT(V T i，： )

[0065] Where i = 1, 2, ..., m, FFT(U ：， i) indicates selecting the left matrix [U] m×m All elements in the i-th column are computed using the Fast Fourier Transform (FFT), FFT(V T i，： ) represents the right matrix [V T ] n×n All elements in the i-th row are computed using the Fast Fourier Transform (FFT);

[0066] S3-2, in each of the vectors U i (k) and each of the vectors V i The last position of (k) is padded with 0s for each of its elements, so that vector U i (k) and vector V i The length of (k) is N, and vectors are obtained sequentially. sum vector Where i = 1, 2, ..., m;

[0067] S3-3, the m vectors obtained from S3-2 and m vectors Calculate and obtain m Y values i (k), each Y i The formula for calculating (k) is:

[0068]

[0069] Where i = 1, 2, ..., m.

[0070] S3-4, for each Y i (k) is used for calculation to obtain the elementary component matrix y. i :

[0071] y i=IFFT(Y i (k))⊙(σ i L)

[0072] Where IFFT is the inverse fast Fourier transform calculation, σ i =S i，i , σ i Represents the singular value matrix [S] m×n Let L be the singular value in the i-th row and i-th column of a vector, where i = 1, 2, ..., m. Let ⊙ be the symbol for vector multiplication by position. Let L be a vector with N elements, and let L be the l-th element of L. l for:

[0073]

[0074] S4, which divides m basic components y i Divided into Grouping y into the same group i Summing, we get The intrinsic mode function (IMF) components of the original signal x (Adaptively determined by the grouping algorithm of S4), where m basic components y i Divided into The specific steps for grouping are as follows:

[0075] S4-1, Set the discrimination parameter γ∈(0,1), in this embodiment γ=0.5, and assign y1 to y Δ Create the first group, set y1 as the first element of the first group, and set g = 1;

[0076] S4-2, Calculate y Δ and y g+1 99% occupied bandwidth (OBW, bandwidth accounting for 99% of signal power), per group y Δ and y g+1 After calculating 99% of the bandwidth used, two bandwidth usage ranges are obtained: (s) Δ ,e Δ ) and (s g+1 e g+1 );

[0077] S4-3, Calculate (s Δ e Δ ) and (s g+1 e g+1 The intersection of ) yields (s Δ，g+1 e Δ，g+1 );

[0078] (s Δ，g+1 e Δg+1 )=(s Δ eΔ )∩(s g+1 e g+1 )

[0079] S4-4, Calculate the bandwidth intersection ratio (P) OBW ) g :

[0080]

[0081] S4-5, (P OBW ) g Compare with the discriminant parameter γ:

[0082] When (P) OBW ) g When <γ, generate a new group, and y g+1 As the first element of the new grouping and y g+1 Assign a value to y Δ ;

[0083] When (P) OBW ) g ≥γ, y g+1 Add to y g In the group, compare y Δ and y g+1 99% of the bandwidth is used, so assign a value to y that has a wider bandwidth. Δ ;

[0084] S4-6, When g is less than m-1, increment the value of g by one, and then use the value of y obtained in S4-5. Δ Substitute into S4-2 and execute S4-2 to S4-6; when g equals m-1, the calculation ends.

[0085] Example 2

[0086] A singular spectral mode decomposition method for separating unsteady components in a spectral aliasing environment is proposed. Based on Example 1, the original signal x is analog signal 1, analog signal 2, speech signal 1, or a rolling bearing fault signal. Analog signal 1 and analog signal 2 are obtained by computer simulation. Analog signal 1 is a signal composed of a harmonic signal and a triangular frequency modulated signal, and the instantaneous time spectrum of all its components is as follows: Figure 1 As shown on the left, analog signal 2 is a signal composed of a harmonic signal and a second-order swept frequency signal. The instantaneous frequency spectrum of all its components is shown below. Figure 2 As shown on the left.

[0087] The rolling bearing fault signal was acquired by a vibration sensor. A short-time Fourier transform (STFT) was applied to the rolling bearing fault signal (Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.). The resulting spectrogram of the rolling bearing fault signal is shown below. Figure 4 (Above)

[0088] Speech signal 1 was acquired by a microphone and contains the vowel / a / . A short-time Fourier transform (STFT) was applied to speech signal 1 (Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.). The resulting spectrogram of speech signal 1 is shown below. Figure 3 (Above)

[0089] The number of sampling points for analog signal 1, analog signal 2, voice signal 1, and rolling bearing fault signal is 2049.

[0090] The speech signal 1 was used as the original signal x, and the singular spectral mode decomposition method in Example 1 was used for calculation (excluding S4). The time consumption was compared with that of traditional Singular Spectrum Analysis (SSA) (R. Vautard and M. Ghil, "Singular spectrum analysis in nonlinear dynamics, with applications topaleoclimatic time series," Physica D: Nonlinear Phenomena, vol. 35, no. 3, pp. 395–424, 1989). The time consumption of each step and the total time consumption are shown in Table 1. The total time consumption was reduced from 36.28s to 0.25s, and the calculation time was greatly shortened.

[0091] Table 1

[0092]

[0093] The decomposition results of the singular spectrum mode decomposition method of this invention are visualized and analyzed. The Hilbert transform of each obtained intrinsic mode function (IMF) component is calculated (NE Huang, et al., “The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary timeseries analysis,” Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences, vol. 454, no. 1971, pp. 903–995, 1998.), yielding the instantaneous frequency sequence ω and instantaneous amplitude sequence a for each IMF component. Then, the instantaneous frequency sequence ω of all IMF components is plotted on a graph, with the color intensity and magnitude of each point determined by the instantaneous energy sequence e = a. 2 The decision was made, and the instantaneous time spectrum was finally obtained.

[0094] In this embodiment, the instantaneous time spectrum of all true components of analog signal 1 is as follows: Figure 1 As shown on the left, the instantaneous time-frequency diagram of the intrinsic mode function (IMF) components obtained by using the singular spectral mode decomposition method for analog signal 1 is as follows. Figure 1 As shown on the right; the instantaneous time spectrum of all real components of analog signal 2 is as follows. Figure 2 As shown on the left, the instantaneous time-frequency diagram of the intrinsic mode function (IMF) components obtained by using the singular spectral mode decomposition method on analog signal 2 is as follows. Figure 2 (As shown on the right).

[0095] Will Figure 1 (left) and Figure 1 (Right) For comparison Figure 2 (left) and Figure 2 (Right) By comparison, it can be seen from the two sets of comparisons that even if the components of the analog signal overlap in frequency, the singular spectrum mode decomposition method of the present invention can still separate the harmonic and non-steady-state components.

[0096] The instantaneous time-frequency plot of the intrinsic mode function (IMF) components of speech signal 1 was obtained by using the singular spectral mode decomposition method, as shown below. Figure 3 (Below) as shown, with Figure 3(Above) As can be seen from the comparison, the singular spectrum mode decomposition method of the present invention can not only separate the harmonic components in speech signal 1, but also separate the non-steady-state components. The frequency fluctuation period of the non-steady-state component separated from speech signal 1 is about 0.06s, which is comparable to the fluctuation period of the non-steady-state component corresponding to the bubbly sound in speech signal 1.

[0097] The instantaneous time-frequency plots of the intrinsic mode function (IMF) components were obtained by using the singular spectrum mode decomposition method to calculate the rolling bearing fault signal, as shown below. Figure 4 (Below) as shown, with Figure 4 (Above) As can be seen from the comparison, the singular spectrum mode decomposition method of the present invention can separate the unsteady components corresponding to the fault components. The time-frequency fluctuation frequency of the separated unsteady components is 106Hz, which is comparable to the fault characteristic frequency of the rolling bearing.

[0098] The present invention has been described above by way of example. It should be noted that any simple modifications, alterations or other equivalent substitutions that can be made by those skilled in the art without creative effort without departing from the core of the present invention fall within the protection scope of the present invention.

Claims

1. A singular spectral mode decomposition method for separating unsteady components in a spectral aliasing environment, characterized in that, Includes the following steps: S1, reconstruct all elements of the original signal x into a high-rank Hankel matrix, obtaining the high-rank Hankel matrix X: Where, x k Let be the k-th element in the original signal x, where k = 1, 2, ..., N, N is the number of sampling points in the original signal x, and m is the rank of the high-rank Hankel matrix. n = N - m + 1; The original signal x is a one-dimensional time series data obtained by sampling the vibration signal. The original signal x is defined as a signal composed of multiple intrinsic mode functions; the vibration signal is a speech signal or a mechanical vibration signal. S2, Perform singular value decomposition on the high-rank Hankel matrix X to obtain a left matrix, a singular value matrix, and a right matrix, where the left matrix is ​​matrix [U]. m×m The singular value matrix is ​​matrix [S]. m×n The right matrix is ​​matrix [V] T ] n×n ; S3, for the left matrix, singular value matrix, and right matrix, performs fast singular spectral decomposition, including: S3-1 to S3-4. S3-1, based on the left and right matrices, calculate m vectors U. i (k) and m vectors V i (k), each vector U i (k) is represented as follows: U i (k)=FFT(U ：，i ), each vector V i (k) is represented as follows: V i (k)=FFT(V T i,： ), where i = 1, 2, ..., m, FFT(U ：，i ) indicates selecting the left matrix [U] m×m All elements in the i-th column are computed using the Fast Fourier Transform (FFT(V)). T i，： ) represents the right matrix [V T ] n×n All elements in the i-th row are computed using the Fast Fourier Transform; S3-2, in each of the vectors U i (k) and each of the vectors V i The last position of (k) is padded with 0s for each of its elements, so that vector U i (k) and vector V i The length of (k) is N, and vectors are obtained sequentially. sum vector S3-3, the m vectors obtained from S3-2 and m vectors Calculate and obtain m Y values i (k), each Y i The formula for calculating (k) is: S3-4, for each Y i (k) is used to calculate and obtain the basic component y. i y i =IFFT(Y i (k))⊙(σ i L); Where IFFT is the inverse fast Fourier transform calculation, σ i Represents the singular value matrix [S] m×n Let L be the singular value in the i-th row and i-th column of a vector, ⊙ be the symbol for vector multiplication, and L be a vector with N elements, where L is the l-th element of L. l for: S4, which divides m basic components y i Divided into Grouping y into the same group i Summing, we get The intrinsic mode function components of an original signal x, where m fundamental components y i Divided into The specific steps for grouping are as follows: S4-1, Set the discrimination parameter γ∈(0,1), and assign y1 to y Δ Create the first group, set y1 as the first element of the first group, and set g = 1; S4-2, Calculate y Δ and y g+1 99% of the bandwidth is used, per group y Δ and y g+1 After calculating 99% of the bandwidth used, two bandwidth usage ranges are obtained: (s) Δ e Δ ) and (s g+1 e g+1 ); S4-3, Calculate (s Δ e Δ ) and (s g+1 e g+1 The intersection of ) yields (s Δ,g+1 e Δ,g+1 ); (s Δ,g+1 ,e Δ，g+1 )=(s Δ ,e Δ )∩(s g+1 ,e g+1 ) S4-4, Calculate the bandwidth intersection ratio (P) OBW ) g : S4-5, (P OBW ) g Compare with the discriminant parameter γ: When (P) OBW ) g When <γ, generate a new group, and y g+1 As the first element of the new grouping and y g+1 Assign a value to y Δ ; When (P) OBW ) g ≥γ, y g+1 Add to y g In the group, compare y Δ and y g+1 99% of the bandwidth is used, so assign a value to y that has a wider bandwidth. Δ ; S4-6, When g is less than m-1, increment the value of g by one, and then use the value of y obtained in S4-5. Δ Substitute into S4-2 and execute S4-2 to S4-6; when g equals m-1, the calculation ends.

2. The singular spectral mode decomposition method according to claim 1, characterized in that, In S1, the functional expression x(t) of the original signal x is: Where, r j (t) is the functional expression of the j-th intrinsic mode function, j = 1, 2, ..., J, where J is the number of intrinsic mode functions contained in the original signal x.

3. The singular spectral mode decomposition method according to claim 2, characterized in that, Among them, a j (t) is the amplitude modulation function in the j-th intrinsic mode function, f j (τ) is the frequency modulation function in the j-th intrinsic mode function.

Images