Frequency estimation method for data-driven subspace containing color noise

Abstract

The invention discloses a method for estimating the frequency of a data-driven subspace containing color noise, and belongs to the field of digital signals. According to the method, the problem that the accuracy is low when the frequency of the signal under the color noise is estimated by adopting an existing subspace frequency estimation method is solved. The invention provides a subspace frequency estimation method under color noise for the problem of color noise in signals, a nonlinear discrete time system is used for describing the color noise, a Hankel matrix is constructed according to the number of sample groups, past step length and future step length, an extended input and output model of the nonlinear discrete time system is constructed, and the frequency of the color noise is estimated. And extracting a system matrix through a gap measurement method, further extracting a frequency characteristic matrix A, solving characteristic values of the frequency characteristic matrix A, and determining a periodic signal frequency. The method provided by the invention can be applied to the frequency estimation of the data-driven subspace containing the color noise.

Description

Technical Field

[0001] The invention belongs to the field of digital signals, and in particular relates to a data-driven subspace frequency estimation method containing colored noise. Background Art

[0002] With the rapid development of science and technology, more and more industrial systems are under high-intensity and full-load operation conditions, generating massive amounts of operational data. In order to obtain system or equipment information in a timely manner, sensors are usually installed to collect data, such as drone blade vibration sensors, high-voltage power grid voltage sensors, motor bearing vibration sensors, etc. The safety and reliability of the system are analyzed by analyzing the data of the core components of the equipment. However, in complex environments, the collected data has the characteristics of periodic oscillation and contains colored noise. Therefore, considering the signal problem under colored noise, accurately estimating the frequency of the periodic oscillation signal under colored noise is of great significance for the fault diagnosis and safe operation of the system.

[0003] At present, the frequency estimation methods of periodic oscillation signals can be divided into non-parametric estimation and parameter estimation. The commonly used non-parametric estimation method is the fast Fourier transform method for frequency estimation. Specifically, the time domain and frequency characteristics of the signal cannot be linked in the fast Fourier transform method. When the amount of data is limited, the spectrum leakage problem cannot be avoided, resulting in inaccurate frequency estimation. However, as an important branch of parameter estimation methods, subspace-assisted frequency estimation methods have received increasing attention in the past few decades. This method can effectively estimate the main frequency of the signal more accurately, link the time domain with the frequency domain, and avoid the spectrum leakage problem. However, the existing subspace frequency estimation method is still not suitable for the frequency estimation of periodic oscillation signals under colored noise. The accuracy of using the existing subspace frequency estimation method to estimate the frequency of periodic oscillation signals under colored noise is low. Therefore, combined with the problem of colored noise, designing a subspace frequency estimation method to estimate the frequency of periodic oscillation signals is a problem that needs to be solved urgently. Summary of the invention

[0004] The purpose of the present invention is to solve the problem of low accuracy when the frequency of a signal under colored noise is estimated by using the existing subspace frequency estimation method, and to propose a data driven subspace frequency estimation method containing colored noise.

[0005] The technical solution adopted by the present invention to solve the above technical problems is: a data-driven subspace frequency estimation method containing color noise, the method specifically comprising the following steps:

[0006] Step 1: construct the Hankel matrix according to the signals at each sampling moment;

[0007] Step 2: Perform LQ decomposition on the Hankel matrix constructed in step 1, and calculate the parameter matrix of the system containing colored noise according to the LQ decomposition result;

[0008] Step 3: Extract the extended observable matrix Γ based on the parameter matrix calculated in step 2 sf ;

[0009] Step 4: According to the extended observable matrix Γ sf Calculate the frequency characteristic matrix A, and then estimate the signal frequency based on the frequency characteristic matrix A.

[0010] The beneficial effects of the present invention are:

[0011] Aiming at the problem of colored noise in signals, the present invention proposes a subspace frequency estimation method under colored noise, uses a nonlinear discrete time system to describe the colored noise, constructs a Hankel matrix according to the number of sample groups, past step lengths and future step lengths, constructs an extended input-output model of the nonlinear discrete time system, extracts the system matrix through a gap measurement method, and further extracts a frequency characteristic matrix A, solves its eigenvalues ​​and determines the frequency of the periodic signal.

[0012] The method for estimating the fundamental frequency of a periodic signal designed by the present invention takes into account the influence of colored noise in the signal, avoids the disadvantage of spectrum leakage in the previous fast Fourier transform method, and improves the accuracy of frequency estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 is a flow chart of a data driven subspace frequency estimation method containing colored noise of the present invention;

[0014] Figure 2 This is the result diagram of system order determination. DETAILED DESCRIPTION

[0015] Specific implementation method 1: Combination Figure 1 The present embodiment is described as follows. The present embodiment describes a method for estimating frequency of a data-driven subspace containing color noise, and the method specifically comprises the following steps:

[0016] Step 1: construct the Hankel matrix according to the signals at each sampling moment;

[0017] Step 2: Perform LQ decomposition on the Hankel matrix constructed in step 1, and calculate the parameter matrix of the system containing colored noise according to the LQ decomposition result;

[0018] Step 3: Extract the extended observable matrix Γ based on the parameter matrix calculated in step 2 sf , the matrix Γ sf Used to analyze the observability of the system;

[0019] Step 4: According to the extended observable matrix Γ sf Calculate the frequency characteristic matrix A, and then estimate the signal frequency based on the frequency characteristic matrix A.

[0020] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the signal at each sampling time is:

[0021]

[0022] Among them, n f represents the number of cosine components, α i is the amplitude of the i-th cosine component, ω i is the frequency of the i-th cosine component, is the phase of the i-th cosine component, t s is the time interval between two adjacent sampling moments, y k is the signal at the kth sampling moment, v k is the colored noise at the kth sampling moment.

[0023] The other steps and parameters are the same as those in the first embodiment.

[0024] Specific implementation method three: This implementation method is different from specific implementation method one or two in that the color noise v at the kth sampling moment k for:

[0025] v k =h(x n,k ,n k )=h C (x n,k )x n,k +h D (x n,k ) k (2)

[0026] Among them, x n,k is the state of generating colored noise at the kth sampling moment (the state at the initial moment is obtained by random initialization), x n,k+1 is the state where color noise is generated at the k+1th sampling moment, x n,k+1 =f(x n,k ,n k ) = f A (x n,k )x n,k +f B (x n,k ) k , n k is the Gaussian white noise at the kth sampling moment, f A (x n,k ), fB (x n,k ),h C (x n,k ) and h D (x n,k ) is the color noise system parameter matrix, and has the characteristics of time variation;

[0027]

[0028] h C (x n,k )=[sin(k) cos(k)] (5)

[0029] h D (x n,k )=[cos(k)] (6)

[0030] The other steps and parameters are the same as those in the first or second embodiment.

[0031] Specific implementation method 4: This implementation method is different from one of the specific implementation methods 1 to 3 in that the specific process of step 1 is as follows:

[0032] Step 1: Construct the Hankel matrix based on the signals at each sampling time and

[0033]

[0034] Where R represents the real number domain, represents the Hankel matrix The number of rows is s p , the number of columns is N, represents the Hankel matrix The number of rows is s f , the number of columns is N, N is the number of given sample groups, and is the intermediate variable of the Hankel matrix;

[0035] Step 1 and 2: Construct the Hankel matrix based on the signals at each sampling time and

[0036]

[0037] Where R represents the real number domain, represents the Hankel matrix The number of rows is s p , the number of columns is N, represents the Hankel matrix The number of rows is s f, the number of columns is N, N is the number of given sample groups, and is the intermediate variable of the Hankel matrix.

[0038] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0039] In the present invention, the value of t is any positive integer.

[0040] Specific implementation method 5: This implementation method is different from the specific implementation methods 1 to 4 in that the intermediate variable of the Hankel matrix and They are:

[0041]

[0042]

[0043] Among them, s p is the past step length, s f For future step length, For the ks p The signal at the sampling time, y k-1 is the signal at the k-1th sampling moment, For the ks p +1 sampling time signal, y k is the signal at the kth sampling moment, For the ks p +N-1 sampling time signal, y k+N-2 is the signal at the k+N-2th sampling moment, For the ks f +1 sampling time signal, y k+1 is the signal at the k+1th sampling moment, is the k+sth f The signal at the sampling time, y k+N-1 is the signal at the k+N-1th sampling moment, is the k+sth f +N-2 signal sampling times.

[0044] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0045] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the intermediate variable of the Hankel matrix and They are:

[0046]

[0047]

[0048] in, is the k+ts p The signal at the sampling time, y k+t-1 is the signal at the k+t-1th sampling time, is the k+ts p +1 sampling time signal, y k+t is the signal at the k+tth sampling moment, is the k+ts p +N-1 sampling time signal, y k+t+N-2 is the signal at the k+t+N-2th sampling time, y k+t is the signal at the k+tth sampling moment, is the k+t+s f -1 sampling time signal, y k+t+1 is the signal at the k+t+1th sampling time, is the k+t+s f The signal at the sampling time, y k+t+N-1 is the signal at the k+t+N-1th sampling time, is the k+t+s f +N-2 signal sampling times.

[0049] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0050] Specific implementation method 7: This implementation method is different from the specific implementation methods 1 to 6 in that the specific process of step 2 is as follows:

[0051] For the Hankel matrix and Do LQ decomposition:

[0052]

[0053] Among them, L 11 , L 21 and L 22 are all sub-blocks of the lower triangular matrix obtained by LQ decomposition, Q1 and Q2 are unitary matrix sub-blocks, and They are the Hankel matrices constructed from past and future steps respectively;

[0054] For the matrix L 21 Q1 performs singular value decomposition:

[0055]

[0056] Wherein, U is a left singular matrix, U1 and U2 are sub-block matrices of U; V is a right singular matrix, V1 and V2 are sub-block matrices of V; Σ is a diagonal matrix, Σ1 and Σ2 are sub-block matrices of Σ, and Σ2≈0;

[0057] Then the parameter matrix Γ of the system containing colored noise at the kth sampling moment is sf,k for:

[0058]

[0059] For the Hankel matrix and Do LQ decomposition:

[0060]

[0061] Among them, L′ 11 , L′ 21 and L′ 22 They are all sub-blocks of the lower triangular matrix obtained by LQ decomposition, and Q′1 and Q′2 are unitary matrix sub-blocks;

[0062] For the matrix L′ 21 Q′1 performs singular value decomposition:

[0063]

[0064] Wherein, U′ is a left singular matrix, U′1 and U′2 are sub-block matrices of U′; V′ is a right singular matrix, V1′ and V2′ are sub-block matrices of V′; Σ′ is a diagonal matrix, Σ′1 and Σ′2 are sub-block matrices of Σ′, and Σ′2≈0;

[0065] Then the parameter matrix Γ of the system containing colored noise at the k+tth sampling time is sf,k+t for:

[0066]

[0067] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0068] Specific implementation eight: This implementation differs from specific implementations one to seven in that the specific process of step three is as follows:

[0069] Step 31: Initialize i=1;

[0070] Step 3.2: Parameter matrix Γ sf,k and Γ sf,k+t Perform singular value decomposition on the first 2i columns of the subspace:

[0071]

[0072] in, is the parameter matrix Γ sf,k The first 2i columns of the subspace, is the parameter matrix Γ sf,k+t The first 2i columns of the subspace, and For The sub-block matrix of the left unitary matrix obtained by singular value decomposition, and For The sub-block matrix of the left unitary matrix obtained by singular value decomposition, Σ ki and Σ (k+t)i is a diagonal matrix, For The right unitary matrix obtained by singular value decomposition is For The right unitary matrix obtained by singular value decomposition;

[0073] Step 3. Calculate Γ sf,k and Γ sf,k+t The gap measurement result of i (Γ sf,k ,Γ sf,k+t )for:

[0074]

[0075] Among them, the superscript T represents the transpose of the matrix;

[0076] If gap i (Γ sf,k ,Γ sf,k+t ) is less than or equal to the set threshold ε, then set i=i+1 and return to step 32;

[0077] If gap i (Γ sf,k ,Γ sf,k+t ) is greater than the set threshold ε, then n f =i-1, the extracted extended observable matrix Γ sf =Γ sf,k (:,1:2n f ), that is, extract the matrix Γ sf,k The first 2n f List.

[0078] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0079] The threshold ε in this embodiment is an arbitrarily small positive number.

[0080] Specific implementation method 9: This implementation method is different from any one of the specific implementation methods 1 to 8 in that the extended observable matrix Γ sf Calculate the frequency feature matrix A, specifically:

[0081]

[0082] in, represents the pseudo-inverse of the matrix, Γ sf (2:s f ,:) represents the matrix Γ sf From the 2nd to the sth row f OK, Representation Matrix From row 1 to row s f -1 line.

[0083] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0084] Specific implementation method ten: This implementation method is different from any one of specific implementation methods one to nine in that the signal frequency is estimated according to the frequency characteristic matrix A, specifically:

[0085]

[0086] Among them, λ n is the eigenvalue of the frequency characteristic matrix A, Re(λ n ) represents the eigenvalue λ n The real part of is the frequency estimation result of the periodic signal, is the frequency estimation result of the i-th cosine component, i=1,2,...,n f .

[0087] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0088] Example

[0089] This embodiment proposes a data-driven subspace frequency estimation method containing color noise, and the method is specifically as follows:

[0090] Step 1: Obtain the motor bearing vibration sensor signal at each sampling moment, and the signal collected by the vibration sensor is the vibration velocity or vibration displacement:

[0091]

[0092] Among them, n f represents the number of cosine components, α i is the amplitude of the i-th cosine component, ω i is the frequency of the i-th cosine component, is the phase of the i-th cosine component, t s is the time interval between two adjacent sampling moments, y k is the signal at the kth sampling moment, v k is the colored noise at the kth sampling moment.

[0093] v k =h(x n,k ,n k )=h C (x n,k )x n,k +h D (x n,k ) k (2)

[0094] Among them, x n,k is the state of generating colored noise at the kth sampling moment (the state at the initial moment is obtained by random initialization), x n,k+1 is the state where color noise is generated at the k+1th sampling moment, x n,k+1 =f(x n,k ,n k ) = f A (x n,k )x n,k +f B (x n,k ) k , n k is the Gaussian white noise at the kth sampling moment, f A (x n,k ), f B (x n,k ),h C (x n,k ) and h D (x n,k ) is the color noise system parameter matrix, and has the characteristics of time variation;

[0095]

[0096] h C (x n,k )=[sin(k) cos(k)] (5)

[0097] h D (x n,k )=[cos(k)] (6)

[0098] Then, the Hankel matrix is ​​constructed according to the signals obtained at each sampling moment; specifically:

[0099] Step 1: Construct the Hankel matrix based on the signals at each sampling time and

[0100]

[0101] Where R represents the real number domain, represents the Hankel matrix The number of rows is s p , the number of columns is N, represents the Hankel matrix The number of rows is s f , the number of columns is N, N is the number of given sample groups, and is the intermediate variable of the Hankel matrix;

[0102]

[0103] Among them, s p is the past step length, s f For future step length, For the ks p The signal at the sampling time, y k-1 is the signal at the k-1th sampling moment, For the ks p +1 sampling time signal, y k is the signal at the kth sampling moment, For the ks p +N-1 sampling time signal, y k+N-2 is the signal at the k+N-2th sampling moment, For the ks f +1 sampling time signal, y k+1 is the signal at the k+1th sampling moment, is the k+sth f The signal at the sampling time, y k+N-1 is the signal at the k+N-1th sampling moment, is the k+sth f +N-2 signal sampling times.

[0104] Step 1 and 2: Construct the Hankel matrix based on the signals at each sampling time and

[0105]

[0106] Where R represents the real number domain, represents the Hankel matrix The number of rows is s p , the number of columns is N, represents the Hankel matrix The number of rows is sf , the number of columns is N, N is the number of given sample groups, and is the intermediate variable of the Hankel matrix;

[0107]

[0108] in, is the k+ts p The signal at the sampling time, y k+t-1 is the signal at the k+t-1th sampling time, is the k+ts p +1 sampling time signal, y k+t is the signal at the k+tth sampling moment, is the k+ts p +N-1 sampling time signal, y k+t+N-2 is the signal at the k+t+N-2th sampling time, y k+t is the signal at the k+tth sampling moment, is the k+t+s f -1 sampling time signal, y k+t+1 is the signal at the k+t+1th sampling time, is the k+t+s f The signal at the sampling time, y k+t+N-1 is the signal at the k+t+N-1th sampling time, is the k+t+s f +N-2 signal sampling times.

[0109] Step 2: Perform LQ decomposition on the Hankel matrix constructed in step 1, and calculate the parameter matrix of the system containing colored noise according to the LQ decomposition result;

[0110] For the Hankel matrix and Do LQ decomposition:

[0111]

[0112] Among them, L 11 , L 21 and L 22 are all sub-blocks of the lower triangular matrix obtained by LQ decomposition, Q1 and Q2 are unitary matrix sub-blocks, and They are the Hankel matrices constructed from past and future steps respectively;

[0113] For the matrix L 21 Q1 performs singular value decomposition:

[0114]

[0115] Wherein, U is a left singular matrix, U1 and U2 are sub-block matrices of U; V is a right singular matrix, V1 and V2 are sub-block matrices of V; Σ is a diagonal matrix, Σ1 and Σ2 are sub-block matrices of Σ, and Σ2≈0;

[0116] Then the parameter matrix Γ of the system containing colored noise at the kth sampling moment is sf,k for:

[0117]

[0118] For the Hankel matrix and Do LQ decomposition:

[0119]

[0120] Among them, L′ 11 , L′ 21 and L′ 22 They are all sub-blocks of the lower triangular matrix obtained by LQ decomposition, and Q′1 and Q′2 are unitary matrix sub-blocks;

[0121] For the matrix L′ 21 Q′1 performs singular value decomposition:

[0122]

[0123] Wherein, U′ is a left singular matrix, U′1 and U′2 are sub-block matrices of U′; V′ is a right singular matrix, V1′ and V2′ are sub-block matrices of V′; Σ′ is a diagonal matrix, Σ′1 and Σ′2 are sub-block matrices of Σ′, and Σ′2≈0;

[0124] Then the parameter matrix Γ of the system containing colored noise at the k+tth sampling time is sf,k+t for:

[0125]

[0126] Step 3: Extract the extended observable matrix Γ based on the parameter matrix calculated in step 2 sf , the matrix Γ sf Used to analyze the observability of the system; the specific process of step three is:

[0127] Step 31: Initialize i=1;

[0128] Step 3.2: Parameter matrix Γ sf,k and Γ sf,k+t Perform singular value decomposition on the first 2i columns of the subspace:

[0129]

[0130] in, is the parameter matrix Γ sf,k The first 2i columns of the subspace, is the parameter matrix Γ sf,k+t The first 2i columns of the subspace, and For The sub-block matrix of the left unitary matrix obtained by singular value decomposition, and For The sub-block matrix of the left unitary matrix obtained by singular value decomposition, Σ ki and Σ (k+t)i is a diagonal matrix, For The right unitary matrix obtained by singular value decomposition is For The right unitary matrix obtained by singular value decomposition;

[0131] Step 3. Calculate Γ sf,k and Γ sf,k+t The gap measurement result of i (Γ sf,k ,Γ sf,k+t )for:

[0132]

[0133] Among them, the superscript T represents the transpose of the matrix;

[0134] If gap i (Γ sf,k ,Γ sf,k+t ) is less than or equal to the set threshold ε, then set i=i+1 and return to step 32;

[0135] If gap i (Γ sf,k ,Γ sf,k+t ) is greater than the set threshold ε, then n f =i-1, the extracted extended observable matrix Γ sf =Γ sf,k (:,1:2n f ), that is, extract the matrix Γ sf,k The first 2n f List.

[0136] Step 4: According to the extended observable matrix Γ sf Calculate the frequency feature matrix A:

[0137]

[0138] in, represents the pseudo-inverse of the matrix, Γ sf (2:s f ,:) represents the matrix Γ sf From the 2nd to the sth row f OK, Representation Matrix From row 1 to row s f -1 line.

[0139] Then estimate the signal frequency based on the frequency characteristic matrix A:

[0140]

[0141] Among them, λ n is the eigenvalue of the frequency characteristic matrix A, Re(λ n ) represents the eigenvalue λ n The real part of is the frequency estimation result of the periodic signal, is the frequency estimation result of the i-th cosine component, i=1,2,...,n f .

[0142] The effectiveness of the method of the present invention will be illustrated below in conjunction with specific simulation results:

[0143] According to formula (1), a set of periodic signal data is constructed, with frequency ω = 2πf, f = [46, 42, 30, 18, 8], amplitude α = [30, 25, 15, 10, 5], and phase The nonlinear state space parameters of the color noise are set as:

[0144]

[0145] h C (x n,k )=[sin(k)cos(k)]

[0146] h D (x n,k )=[cos(k)]

[0147] Therefore, the nonlinear state space model of colored noise is:

[0148]

[0149] v k =[sin(k)cos(k)]x n,k +[cos(k)]n k

[0150] The initial state is x n,0 =[0.10.1] T, n k is Gaussian white noise, and its energy is 1×10 -4 .

[0151] Figure 2 The result diagram of the system order calculated by the method of the present invention, matrix Γ sf By the expansion matrix Γ sf,k The simulation experiment was repeated 50 times to avoid the accidental results of the system order. The experimental results show that when the gap metric is about 0, the system order is determined to be 10. That is, the number of frequencies n in the signal f = 5. In this case, the data length N is set to 800, and the data length s p and the future length s f They are set to 70 and 30 respectively, and the frequency component estimation results are shown in Table 1:

[0152] Table 1 Frequency estimation results

[0153]

[0154]

[0155] The simulation results show that the frequency estimation method of the present invention has an extremely high estimation accuracy.

[0156] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. A data-driven subspace frequency estimation method containing colored noise, characterized in that: The method specifically comprises the following steps: Step 1: construct the Hankel matrix according to the signals at each sampling moment; Step 2: Perform LQ decomposition on the Hankel matrix constructed in step 1, and calculate the parameter matrix of the system containing colored noise according to the LQ decomposition result; Step 3: Extract the extended observable matrix Γ based on the parameter matrix calculated in step 2 sf ; Step 4: According to the extended observable matrix Γ sf Calculate the frequency characteristic matrix A, and then estimate the signal frequency based on the frequency characteristic matrix A.

2. The method for estimating frequency of a data-driven subspace containing color noise according to claim 1, characterized in that: The signals at each sampling moment are: Among them, n f represents the number of cosine components, α i is the amplitude of the i-th cosine component, ω i is the frequency of the i-th cosine component, is the phase of the i-th cosine component, t s is the time interval between two adjacent sampling moments, y k is the signal at the kth sampling moment, v k is the colored noise at the kth sampling moment.

3. The method for estimating frequency of a data-driven subspace containing color noise according to claim 2, characterized in that: The colored noise v at the kth sampling moment k for: v k =h C (x n,k )x n,k +h D (x n,k )n k (2) Among them, x n,k is the state of generating colored noise at the kth sampling moment, x n,k+1 is the state where color noise is generated at the k+1th sampling moment, x n,k+1 =f A (x n,k )x n,k +f B (x n,k ) k , n k is the Gaussian white noise at the kth sampling moment, f A (x n,k ), f B (x n,k ),h C (x n,k ) and h D (x n,k ) is the color noise system parameter matrix; h C (x n,k )=[sin(k) cos(k)] (5)。 h D (x n,k )=[cos(k)] (6)。 4. The method for estimating frequency of a data-driven subspace containing color noise according to claim 3, characterized in that: The specific process of step one is: Step 1: Construct the Hankel matrix based on the signals at each sampling time and Where R represents the real number domain, represents the Hankel matrix The number of rows is s p , the number of columns is N, represents the Hankel matrix The number of rows is s f , the number of columns is N, N is the given number of sample groups, and is the intermediate variable of the Hankel matrix; Step 1 and 2: Construct the Hankel matrix based on the signals at each sampling time and Where R represents the real number domain, represents the Hankel matrix The number of rows is s p , the number of columns is N, represents the Hankel matrix The number of rows is s f , the number of columns is N, and is the intermediate variable of the Hankel matrix.

5. The method for estimating frequency of a data-driven subspace containing color noise according to claim 4, characterized in that: The intermediate variable of the Hankel matrix and They are: Among them, s p is the past step length, s f For future step length, For the ks p The signal at the sampling time, y k-1 is the signal at the k-1th sampling moment, For the ks p +1 sampling time signal, y k is the signal at the kth sampling moment, For the ks p +N-1 sampling time signal, y k+N-2 is the signal at the k+N-2th sampling moment, For the ks f +1 sampling time signal, y k+1 is the signal at the k+1th sampling moment, is the k+sth f The signal at the sampling time, y k+N-1 is the signal at the k+N-1th sampling moment, is the k+sth f +N-2 signal sampling times.

6. The method for estimating frequency of a data-driven subspace containing color noise according to claim 5, characterized in that: The intermediate variable of the Hankel matrix and They are: in, is the k+ts p The signal at the sampling time, y k+t-1 is the signal at the k+t-1th sampling time, is the k+ts p +1 sampling time signal, y k+t is the signal at the k+tth sampling moment, is the k+ts p +N-1 sampling time signal, y k+t+N-2 is the signal at the k+t+N-2th sampling time, y k+t is the signal at the k+tth sampling moment, is the k+t+s f -1 sampling time signal, y k+t+1 is the signal at the k+t+1th sampling time, is the k+t+s f The signal at the sampling time, y k+t+N-1 is the signal at the k+t+N-1th sampling time, is the k+t+s f +N-2 signal sampling times.

7. The method for estimating frequency of a data-driven subspace containing color noise according to claim 6, characterized in that: The specific process of step 2 is as follows: For the Hankel matrix and Do LQ decomposition: Among them, L 11 , L 21 and L 22 They are all sub-blocks of the lower triangular matrix obtained by LQ decomposition, and Q1 and Q2 are unitary matrix sub-blocks; For the matrix L 21 Q1 performs singular value decomposition: Among them, U is the left singular matrix, U1 and U2 are the sub-block matrices of U; V is the right singular matrix, V1 and V2 are the sub-block matrices of V; Σ is a diagonal matrix, Σ1 and Σ2 are the sub-block matrices of Σ; Then the parameter matrix Γ of the system containing colored noise at the kth sampling moment is sf,k for: For the Hankel matrix and Do LQ decomposition: Among them, L′ 11 , L′ 21 and L′ 22 They are all sub-blocks of the lower triangular matrix obtained by LQ decomposition, and Q′1 and Q′2 are unitary matrix sub-blocks; For the matrix L′ 21 Q′1 performs singular value decomposition: Among them, U′ is a left singular matrix, U′1 and U′2 are sub-block matrices of U′; V′ is a right singular matrix, V1′ and V2′ are sub-block matrices of V′; Σ′ is a diagonal matrix, Σ′1 and Σ′2 are sub-block matrices of Σ′; Then the parameter matrix Γ of the system containing colored noise at the k+tth sampling time is sf,k+t for: C sf,k+t =U′1Σ′1 12 (26)。 8. The method for estimating frequency of a data-driven subspace containing color noise according to claim 7, characterized in that: The specific process of step three is: Step 31: Initialize i=1; Step 3.2: Parameter matrix Γ sf,k and Γ sf,k+t Perform singular value decomposition on the first 2i columns of the subspace: in, is the parameter matrix Γ sf,k The first 2i columns of the subspace, is the parameter matrix Γ sf,k+t The first 2i columns of the subspace, and For The sub-block matrix of the left unitary matrix obtained by singular value decomposition, and For The sub-block matrix of the left unitary matrix obtained by singular value decomposition, Σ ki and Σ (k+t)i is a diagonal matrix, For The right unitary matrix obtained by singular value decomposition is For The right unitary matrix obtained by singular value decomposition; Step 3. Calculate Γ sf,k and Γ sf,k+t The gap measurement result of i (Γ sf,k ,Γ sf,k+t )for: Among them, the superscript T represents the transpose of the matrix; If gap i (Γ sf,k ,Γ sf,k+t ) is less than or equal to the set threshold ε, then set i=i+1 and return to step 32; If gap i (Γ sf,k ,Γ sf,k+t ) is greater than the set threshold ε, then n f =i-1, the extracted extended observable matrix Γ sf =Γ sf,k (:,1:2n f ).

9. The method for estimating frequency of a data-driven subspace containing color noise according to claim 8, characterized in that: According to the extended observable matrix Γ sf Calculate the frequency feature matrix A, specifically: in, represents the pseudo-inverse of the matrix, Γ sf (2:s f ,:) represents the matrix Γ sf From the 2nd to the sth row f OK, Representation Matrix From row 1 to row s f -1 line.

10. The method for estimating frequency of a data-driven subspace containing color noise according to claim 9, characterized in that: The signal frequency is estimated according to the frequency characteristic matrix A, specifically: Among them, λ n is the eigenvalue of the frequency characteristic matrix A, Re(λ n ) represents the eigenvalue λ n The real part of is the frequency estimation result of the periodic signal, is the frequency estimation result of the i-th cosine component, i=1,2,…,n f .

Images