Singular spectrum modal synthesis algorithm for separating unsteady components in spectrum aliasing environment

Abstract

The invention discloses a singular spectrum modal synthesis algorithm for separating unsteady state components in a spectrum aliasing environment, which comprises the following steps of: reconstructing all elements in an 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 aiming at the left matrix, the singular value matrix and the right matrix, calculating a singular value of the left matrix, the singular value matrix and the right matrix; according to the singular spectrum modal comprehensive algorithm, fast singular spectrum decomposition is carried out, m basic components yi are divided into # imgabs0 # groups, yi divided into the same group is summed, and an intrinsic mode function component of # imgabs1 # original signals x. According to the singular spectrum modal comprehensive algorithm, fast Fourier transform is utilized for fast singular spectrum decomposition, a large amount of matrix calculation in singular spectrum analysis is replaced, the calculation time is greatly shortened, and the calculation efficiency is improved. While the calculation speed is ensured, non-steady-state intrinsic mode function components overlapped in a frequency domain can be separated, and non-steady-state components in voice and mechanical vibration signals can be separated.

Description

Technical Field

[0001] The present invention belongs to the technical field of non-steady-state signal processing, and in particular relates to a singular spectrum mode ensemble (SME) algorithm for separating non-steady-state components in a spectrum aliasing environment. Background Art

[0002] There are non-steady-state components in speech and mechanical vibration signals, and the characteristics of these components play a very important role in speech quality analysis and mechanical fault diagnosis. These non-steady-state features cannot be extracted using traditional short-time analysis methods such as 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 such as 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 Due to the limitation of decomposition ability, it is unable to fully separate the overlapping components on the spectrum. Therefore, it is often unable to separate the non-steady-state components when analyzing real speech and mechanical vibration signals.

[0003] Singular Spectrum Analysis (SSA) technique (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 matrix to separate trends, amplitude modulation and noise from time series data. Studies have shown that high-rank Hankel matrix is ​​crucial for separating overlapping components in signals because it allows higher frequency domain decomposition resolution. However, high-rank Hankel matrix causes SSA to consume huge computational resources and cause over-decomposition. How to ensure computational speed and adaptively combine over-decomposed components while using high-rank Hankel matrix has become a key issue. Summary of the invention

[0004] In view of the deficiencies of the prior art, an object of the present invention is to provide a high-speed, accurate and fully adaptive singular spectrum mode ensemble (SME) algorithm for separating non-steady-state components in a spectrum aliasing environment.

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

[0006] A singular spectrum modal synthesis algorithm for separating non-stationary components in a spectrum aliasing environment comprises the following steps:

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

[0008]

[0009] Among them, x k is the kth element in the original signal x, k = 1, 2, ..., N, N is the number of sampling points of the original signal x, 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, and the original signal x is set to be a signal composed of multiple intrinsic mode functions;

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

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

[0013]

[0014] Among them, r j (t) is the functional expression of the jth intrinsic mode function, j=1, 2, ..., J, and J is the number of intrinsic mode functions contained in the original signal x.

[0015] In the above technical solution, Among them, a j (t) is the amplitude modulation function in the jth eigenmode function, f i (τ) is the frequency modulation function in the jth eigenmode 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 ​​the matrix [U] m×m , the singular value matrix is ​​the matrix [S] m×n , the right matrix is ​​the matrix [V T ] n×n ;

[0017] S3, performing fast singular spectrum decomposition on the left matrix, singular value matrix and right matrix, including: S3-1 to S3-4,

[0018] S3-1, calculate m vectors U according to the left matrix and the right matrix 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 ) means selecting the left matrix [U] m×m All elements of the ith column of T i，： ) represents the right matrix [VT ] n×n All elements of the i-th row are calculated using the Fast Fourier Transform;

[0023] S3-2, in each of the vectors U i (k) and each of the vectors V i The end position of (k) is filled with 0 for each element in it, so that the vector U i (k) and vector V i The length of (k) is N, and the vectors are obtained respectively. and vector

[0024] S3-3, m vectors obtained according to S3-2 and m vectors Calculate and obtain m Y i (k), each Y i The calculation formula of (k) is:

[0025] S3-4, for each Y i (k) Calculate and obtain the basic component y i :

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

[0027] Among them, IFFT is the inverse fast Fourier transform calculation, σ i represents the singular value matrix [S] m×n The singular value of the i-th row and i-th column in , ⊙ is the vector multiplication symbol, L is a vector containing N elements, and the l-th element L in L l for:

[0028]

[0029] S4, m basic components y i Divide groups, and divide the y into the same group i Sum, we get The intrinsic mode function components of the original signal x, where m basic components y i Divide The specific steps for grouping are as follows:

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

[0031] S4-2, calculate y Δ and g+ 1 99% of the occupied bandwidth, each group y Δ and g+1 After calculating the 99% occupied bandwidth, two occupied bandwidth 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 ), we get (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 <γ, generate a new group and set y g+ 1 as the first element of the new group and y g+1 Assign to y Δ ;

[0038] When (P OBW ) g ≥γ, y g+1 Add to y g In the group, compare y Δ and g+1 99% of the occupied bandwidth, assign an element with a wider 99% occupied bandwidth to y Δ ;

[0039] S4-6, when g is less than m-1, increase the value of g by 1 and add y obtained in S4-5 to the value of Δ Substitute into S4-2 and execute S4-2 to S4-6; when g equals m-1, end the calculation.

[0040] The characteristics and beneficial effects of the present invention are:

[0041] (1) The fast singular spectrum decomposition (FSSD) in the singular spectrum modal synthesis algorithm of the present invention uses fast Fourier transform to replace a large number of matrix calculations in singular spectrum analysis, greatly shortening the calculation time.

[0042] (2) The present invention proposes an adaptive grouping method based on the proportion of singular values ​​and bandwidth intersection, which adaptively groups a large number of over-decomposed elementary components (elementary component matrix) and solves the problem of over-decomposition.

[0043] (3) The singular spectrum modal synthesis (SME) algorithm proposed in the present invention can separate the overlapping non-steady-state intrinsic mode function (IMF) components in the frequency domain while ensuring the calculation speed, and can separate the non-steady-state components in speech and mechanical vibration signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 The instantaneous time-frequency spectrum of the real component of the simulated signal 1 (left) and the instantaneous time-frequency spectrum calculated by the singular spectrum modal synthesis algorithm (right);

[0045] Figure 2 The instantaneous time-frequency spectrum of the real component of the simulated signal 2 (left) and the instantaneous time-frequency spectrum calculated by the singular spectrum modal synthesis algorithm (right);

[0046] Figure 3 The spectrogram of speech signal 1 (top) and the instantaneous time-frequency spectrum (bottom) calculated using the singular spectrum modal synthesis algorithm;

[0047] Figure 4 The acoustic spectrum of the rolling bearing fault signal (top) and the instantaneous time-frequency spectrum calculated using the singular spectrum modal synthesis algorithm (bottom). DETAILED DESCRIPTION

[0048] The singular spectrum modal synthesis (SME) algorithm for separating non-steady-state components in a spectrum aliasing environment of the present invention is described in detail below with reference to the accompanying drawings.

[0049] Example 1

[0050] A singular spectrum modal synthesis algorithm for separating non-stationary components in a spectrum aliasing environment comprises the following steps:

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

[0052]

[0053] Among them, x k is the kth element in the original signal x, k = 1, 2, ..., N, N is the number of sampling points of the original signal x, 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 the vibration signal. The vibration signal may 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). The function expression x(t) of the original signal x is:

[0055]

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

[0057] r j (t) is:

[0058]

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

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

[0061] S3, performing fast singular spectrum decomposition (FSSD) on the left matrix, the singular value matrix and the right matrix, the specific steps are as follows:

[0062] S3-1, calculate m vectors U according to the left matrix and the right matrix i (k) and m vectors V i(k), k represents the phase, 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 ) means selecting the left matrix [U] m×m All elements of the i-th column of are calculated using the Fast Fourier Transform (FFT), FFT(V T i，： ) represents the right matrix [V T ] n×n All elements of the i-th row are calculated 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 end position of (k) is filled with 0 for each element in it, so that the vector U i (k) and vector V i The length of (k) is N, and the vectors are obtained respectively. and vector Wherein, i = 1, 2, ..., m;

[0067] S3-3, m vectors obtained according to S3-2 and m vectors Calculate and obtain m Y i (k), each Y i The calculation formula of (k) is:

[0068]

[0069] Among them, i=1, 2,…, m.

[0070] S3-4, for each Y i (k) Calculate and obtain the elementary component matrix y i :

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

[0072] Among them, IFFT is the inverse fast Fourier transform calculation, σ i =S i，i , σ i represents the singular value matrix [S] m×n The singular value of the i-th row and i-th column in , i = 1, 2, ..., m, ⊙ is the vector multiplication symbol, L is a vector containing N elements, and the l-th element L in L l for:

[0073]

[0074] S4, m basic components y i Divide groups, and divide the y into the same group i Sum, we get The intrinsic mode function (IMF) components of the original signal x ( The m basic components y are adaptively determined by the grouping algorithm of S4. i Divide The specific steps for grouping are as follows:

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

[0076] S4-2, calculate y Δ and g+1 99% occupied bandwidth (OBW, bandwidth that accounts for 99% of the signal power), each group y Δ and g+1 After calculating the 99% occupied bandwidth, two occupied bandwidth 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 ), we get (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 <γ, generate a new group and set y g+1 As the first element of the new group and y g+1 Assign to y Δ ;

[0083] When (P OBW ) g ≥γ, y g+1 Add to y g In the group, compare y Δ and g+ 1 has a bandwidth of 99%, and an element with a wider bandwidth of 99% is assigned to y Δ ;

[0084] S4-6, when g is less than m-1, increase the value of g by 1 and add y obtained in S4-5 to the value of Δ Substitute into S4-2 and execute S4-2 to S4-6; when g equals m-1, end the calculation.

[0085] Example 2

[0086] A singular spectrum modal synthesis algorithm for separating non-steady-state components in a spectrum aliasing environment, based on Example 1, the original signal x is analog signal 1, analog signal 2, voice signal 1 or rolling bearing fault signal. Analog signal 1 and analog signal 2 are obtained by computer simulation respectively, analog signal 1 is a signal composed of a harmonic signal and a triangular frequency modulation signal, and the instantaneous time-frequency spectrum of all its components is as follows Figure 1 (Left) As shown. Analog signal 2 is composed of a harmonic signal and a second-order frequency sweep signal. The instantaneous time spectrum of all its components is shown as Figure 2 (Left) shown.

[0087] The rolling bearing fault signal is collected by the vibration sensor. The rolling bearing fault signal is subjected to short-time Fourier transform (Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.) (short-time Fourier transform, STFT), and the spectrogram of the rolling bearing fault signal is obtained as shown in the figure. Figure 4 (above) as shown;

[0088] The speech signal 1 is obtained by microphone acquisition, and its content is the vowel / a / . The speech signal 1 is subjected to short-time Fourier transform (Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.) (short-time Fourier transform, STFT), and the obtained spectrogram of the speech signal 1 is as follows Figure 3 (above) as shown;

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

[0090] The speech signal 1 is used as the original signal x to calculate the singular spectrum modal synthesis algorithm in Example 1 (excluding S4), and the time consumption is compared with the traditional singular spectrum analysis (Singular Spectrum Analysis, SSA) technology (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 is reduced from 36.28s to 0.25s, and the calculation time is greatly shortened.

[0091] Table 1

[0092]

[0093] The decomposition results of the singular spectrum modal synthesis algorithm of the present invention are visualized and analyzed. The Hilbert transform of each intrinsic mode function (IMF) component is calculated separately (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.), and the instantaneous frequency sequence ω and instantaneous amplitude sequence a of each intrinsic mode function (IMF) component are obtained. Then the instantaneous frequency sequence ω of all intrinsic mode function (IMF) components is plotted on the graph, and the color depth and size of each point are determined by the instantaneous energy sequence e = a 2 Finally, the instantaneous time spectrum diagram is obtained.

[0094] In this embodiment, the instantaneous time spectrum of all real components of the analog signal 1 is as follows: Figure 1 As shown in (left), the instantaneous time-frequency diagram of the intrinsic mode function (IMF) component is calculated by using the singular spectrum mode synthesis algorithm for the analog signal 1 as shown in Figure 1 (right) shows; the instantaneous time spectrum of all real components of analog signal 2 is as follows Figure 2 As shown in (left), the instantaneous time-frequency diagram of the intrinsic mode function (IMF) component is calculated by using the singular spectrum mode synthesis algorithm for the analog signal 2 as shown in Figure 2 (right) shown.

[0095] Will Figure 1 (Left) and Figure 1 (right) for comparison, Figure 2 (Left) and Figure 2 (right) Through the two groups of comparison, it can be seen that even if the components of the analog signal overlap in frequency, the singular spectrum modal synthesis algorithm of the present invention can still separate the harmonics and non-steady-state components.

[0096] The instantaneous time-frequency diagram of the intrinsic mode function (IMF) component is calculated by using the singular spectrum mode synthesis algorithm for the speech signal 1 as shown in the figure: Figure 3 (below) Figure 3(above) By comparison, it can be seen that the singular spectrum modal synthesis algorithm of the present invention can not only separate the harmonic components in the speech signal 1, but also separate the non-steady-state components. The frequency fluctuation period of the non-steady-state components separated from the speech signal 1 is about 0.06s, which is equivalent to the fluctuation period of the non-steady-state components corresponding to the bubble sound in the speech signal 1.

[0097] The instantaneous time-frequency diagram of the intrinsic mode function (IMF) component is calculated by using the singular spectrum mode synthesis algorithm for the rolling bearing fault signal as shown in the figure: Figure 4 (below) Figure 4 (above) From the comparison, it can be seen that the singular spectrum modal synthesis algorithm of the present invention can separate the non-steady-state component corresponding to the fault component, and the fluctuation frequency of the time-frequency curve of the separated non-steady-state component is 106 Hz, which is equivalent to the fault characteristic frequency of the rolling bearing.

[0098] The present invention is described above by way of example. It should be noted that, without departing from the core of the present invention, any simple deformation, modification or other equivalent replacement that can be made by those skilled in the art without inventive effort falls within the protection scope of the present invention.

Claims

1. A singular spectrum modal synthesis algorithm for separating non-steady-state components in a spectrum aliasing environment, characterized in that: The following steps are involved: S1, reconstruct all elements in the original signal x into a high-rank Hankel matrix to obtain a high-rank Hankel matrix X: Among them, x k is the kth element in the original signal x, k = 1, 2, ..., N, N is the number of sampling points of the original signal x, 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, and the original signal x is set to be a signal composed of multiple intrinsic mode functions; 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 ​​the matrix [U] m×m , the singular value matrix is ​​the matrix [S] m×n , the right matrix is ​​the matrix [V T ] n×n ; S3, performing fast singular spectrum decomposition on the left matrix, singular value matrix and right matrix, including: S3-1 to S3-4, S3-1, calculate m vectors U according to the left matrix and the right matrix i (k) and m vectors V i (k), each vector U i (k) is expressed as follows: i (k) = FFT(U ：，i ), each vector V i (k) is expressed as follows: V i (k) = FFT(V T i，： ), where i = 1, 2, ..., m, FFT(U ：，i ) means selecting the left matrix [U] m×m All elements of the ith column of are calculated using the Fast Fourier Transform, FFT(V T i，： ) represents the right matrix [V T ] n×n All elements of the i-th row are calculated using the Fast Fourier Transform; S3-2, in each of the vectors U i (k) and each of the vectors V i The end position of (k) is filled with 0 for each element in it, so that the vector U i (k) and vector V i The length of (k) is N, and the vectors are obtained respectively. and vector S3-3, m vectors obtained according to S3-2 and m vectors Calculate and obtain m Y i (k), each Y i The calculation formula of (k) is: S3-4, for each Y i (k) Calculate and obtain the basic component y i :y i =IFFT(Y i (k))⊙(σ i L); Among them, IFFT is the inverse fast Fourier transform calculation, σ i represents the singular value matrix [S] m×n The singular value of the i-th row and i-th column in , ⊙ is the vector multiplication symbol, L is a vector containing N elements, and the l-th element L in L l for: S4, m basic components y i Divide groups, and divide the y into the same group i Sum, we get The intrinsic mode function components of the original signal x, where m basic components y i Divide The specific steps for grouping are as follows: S4-1, set the discriminant parameter γ∈(0,1), assign y1 to y Δ , create the first group, take y1 as the first element of the first group, and set g = 1; S4-2, calculate y Δ and g+1 99% of the occupied bandwidth, each group y Δ and g+1 After calculating the 99% occupied bandwidth, two occupied bandwidth ranges are obtained: (s Δ , e Δ ) and (s g+1 , e g+1 ); S4-3, calculate (s Δ , e Δ ) and (s g+1 , e g+1 ), we get (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 <γ, generate a new group and set y g+1 As the first element of the new group and y g+1 Assign to y Δ ; When (P OBW ) g ≥γ, y g+1 Add to y g In the group, compare y Δ and g+1 99% of the occupied bandwidth, assign an element with a wider 99% occupied bandwidth to y Δ ; S4-6, when g is less than m-1, increase the value of g by 1 and add y obtained in S4-5 to the value of Δ Substitute into S4-2 and execute S4-2 to S4-6; when g equals m-1, end the calculation.

2. The singular spectrum modal synthesis algorithm according to claim 1, characterized in that: In S1, the vibration signal is a voice signal or a mechanical vibration signal.

3. The singular spectrum modal synthesis algorithm according to claim 1, characterized in that: In S1, the function expression x(t) of the original signal x is: Among them, r j (t) is the functional expression of the jth intrinsic mode function, j=1, 2, ..., J, and J is the number of intrinsic mode functions contained in the original signal x.

4. The singular spectrum modal synthesis algorithm according to claim 3, characterized in that: Among them, a j (t) is the amplitude modulation function in the jth eigenmode function, f j (τ) is the frequency modulation function in the jth eigenmode function.

Images