Alternating projection-based model order unknown meshless line spectrum estimation method

Abstract

The invention discloses a model order unknown meshless line spectrum estimation method based on alternating projection, and belongs to the technical field of line spectrum estimation. According to the method, the problem that joint estimation of the signal frequency and the model order cannot be realized due to the fact that the model order needs to serve as prior information in an existing rapid Vandermonde decomposition and shrinkage threshold algorithm based on alternating projection is solved. According to the method, an alternating projection method is developed under the framework of accelerating near-end gradient descent, gradient descent and near-end mapping are iteratively executed, and near-end mapping of atomic norms has no analytical solution, so that near-end mapping is approximated through alternating projection, and a Toeplitz structure in the atomic norms is approximated. And moreover, the sparsity of the signal is fully utilized, and model order updating is carried out through convergence check, so that the joint estimation of the signal frequency and the model order can be realized while the estimation precision is ensured. The method can be applied to line spectrum estimation.

Description

Technical Field

[0001] The present invention belongs to the technical field of line spectrum estimation, and in particular relates to a gridless line spectrum estimation method with unknown model order based on alternating projection. Background Art

[0002] Line spectrum estimation refers to the frequency estimation and model order selection of multiple superimposed complex exponential signals from a small amount of sample data affected by noise or quantized by low-precision quantizers. Line spectrum estimation problems appear in various application fields, including range and direction estimation in radar, channel estimation in wireless communications, and DOA estimation in array signal processing.

[0003] In the single-measurement vector line spectrum estimation problem, the sample It is obtained by superposing K complex sinusoidal original signals and can be expressed as:

[0004]

[0005] Among them, f k represents the frequency, c k represents the complex amplitude and w represents the noise.

[0006] Classical subspace methods can be used to estimate signal frequencies from measurement data, but they require sufficient signal observation snapshots to approximate the statistics well, and require the model order as prior information. Compared with subspace methods, gridless methods show certain advantages in line spectrum estimation, including single measurement vector applicability and super-resolution. Gridless line spectrum estimation methods need to solve the atomic l0 norm minimization problem. The atomic l0 norm directly exploits sparsity and has no resolution limit, but it is NP-hard to calculate. To make it tractable, most works convexly relax the atomic l0 norm to the atomic norm. Atomic norm minimization needs to be converted into a semidefinite program and solved by the interior point method. The interior point method is too computationally complex, making it difficult to use gridless methods in practical applications. Fast Vandermonde decomposition and shrinkage thresholding algorithms based on alternating projections provide a low-complexity solution and guarantee the solution accuracy to a certain extent. However, this type of method requires the model order as prior information.

[0007] In summary, the existing fast Vandermonde decomposition and shrinkage threshold algorithm based on alternating projection is an accelerated algorithm for gridless line spectrum estimation. Although it reduces the computational complexity of the solution to a certain extent and ensures that the solution accuracy is not affected, it requires the model order as prior information and cannot achieve joint estimation of signal frequency and model order. Summary of the invention

[0008] The purpose of the present invention is to solve the problem that the existing fast Vandermonde decomposition and shrinkage threshold algorithm based on alternating projection requires the model order as prior information and cannot realize the joint estimation of signal frequency and model order, and proposes a gridless line spectrum estimation method with unknown model order based on alternating projection.

[0009] The technical solution adopted by the present invention to solve the above technical problems is: a gridless line spectrum estimation method with unknown model order based on alternating projection, the method specifically comprising the following steps:

[0010] Step 1: Construct the sample autocorrelation matrix R based on sample y yy ;

[0011] Step 2: Set the intermediate variable R0 = R yy , intermediate variable y0=y, intermediate variable v0=tr(R0) / N, and define auxiliary variable Z0 according to R0, y0 and v0, tr(R0) is the trace of R0;

[0012] Step 3: Initialize the model order K = 1, and initialize the number of iterations i = 0;

[0013] Step 4: Calculate the intermediate variable y g =y i -δ(y i -y0), intermediate variable R g =R i , the intermediate variable v g =v i , δ is the moving step length;

[0014] Step 5: According to the model order K, R g Perform truncated eigenvalue decomposition, that is, (Σ K ,U K )=eigs(R g ,K),U K is the matrix composed of eigenvectors obtained by truncated eigenvalue decomposition, Σ K is the vector of eigenvalues ​​obtained by truncated eigenvalue decomposition;

[0015] According to Σ K and U K Calculate the intermediate variable matrix According to the matrix Calculating the Matrix And according to the matrix Z g Get the matrix

[0016] Step 6: According to R i+1 , R i , Z i+1 and Z iCalculate the combined residual

[0017] If i=0, set i=i+1 and return to step 4;

[0018] If i>0, execute step 7;

[0019] Step 7: Compare (NK)λ(K) and η, where λ(K) is the Kth eigenvalue obtained by truncated eigenvalue decomposition;

[0020] If (NK)λ(K)≥η, execute step eight;

[0021] If (NK)λ(K)<η, then execute step nine;

[0022] Step 8. Compare and Size:

[0023] like Then continue to judge whether it satisfies If satisfied The iteration ends and Z i+1 As the line spectrum estimation result The model order K at the last iteration is used as the signal model order estimation result like Then let i=i+1, and return to step 4;

[0024] like Then let K=K+1, i=i+1, and return to step 4;

[0025] Step 9: Compare With the threshold ε tol Size:

[0026] like The iteration ends and Z i+1 As the line spectrum estimation result The model order K at the last iteration is used as the signal model order estimation result

[0027] like Then let K=K-1, i=i+1, and return to execute step 4.

[0028] Furthermore, the specific process of step one is:

[0029] Step 1: Calculate the autocorrelation function of sample y, and record the first, second, ..., Nth elements in the calculation result of the autocorrelation function as r(0), r(1), ..., r(N-1), respectively, where N is the signal length;

[0030] Step 2: Construct the autocorrelation matrix R based on the autocorrelation function calculation results yy :

[0031]

[0032] in, is plural, r * (1) is the transpose of r(1), r * (N-1) is the transpose of r(N-1), r * (N-2) is the transpose of r(N-2), and the autocorrelation matrix R yy The dimension is N×N.

[0033] Furthermore, the combined residual for:

[0034]

[0035] Among them, ||·|| F is the F-norm.

[0036] Furthermore, the auxiliary variable Z0 is:

[0037]

[0038] in, is the conjugate transpose of y0.

[0039] Furthermore, the intermediate variable matrix

[0040] Furthermore, according to the matrix Z g Get the matrix Z i+1 The method used is to g Method of projection onto semidefinite sets.

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

[0042] The present invention develops an alternating projection method under the framework of accelerating proximal gradient descent. The method of the present invention iteratively performs gradient descent and proximal mapping. Since the proximal mapping of the atomic norm has no analytical solution, the present invention chooses to approximate the proximal mapping through alternating projection to approximate the Toeplitz structure in the atomic norm. In addition, the sparsity of the signal is fully utilized, and it is proposed to update the model order through convergence check, so that the method of the present invention can realize the joint estimation of signal frequency and model order while ensuring the estimation accuracy, and further reduce the computational complexity and improve the estimation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 is a comparison chart of the solution accuracy between the method of the present invention and the original gridless line spectrum estimation method;

[0044] Figure 2 is a comparison chart of the success rates of the method of the present invention and the original gridless line spectrum estimation method;

[0045] Figure 3 It is a running time comparison chart of the method of the present invention and the original gridless line spectrum estimation method. DETAILED DESCRIPTION

[0046] Specific implementation method 1: Combination Figure 1 The present embodiment is described as follows. The present embodiment describes a gridless line spectrum estimation method for an unknown model order based on alternating projection, and the method specifically comprises the following steps:

[0047] Step 1: Construct the sample autocorrelation matrix R based on sample y yy ;

[0048] Step 2: Set the intermediate variable R0 = R yy , intermediate variable y0=y, intermediate variable v0=tr(R0) / N, and define auxiliary variable Z0 according to R0, y0 and v0, tr(R0) is the trace of R0;

[0049] Step 3: Initialize the model order K = 1, and initialize the number of iterations i = 0;

[0050] Step 4: Under the framework of the accelerated proximal gradient descent algorithm, y moves along the gradient descent with an appropriate small step size δ to calculate the intermediate variable y g =y i -δ(y i -y0), intermediate variable R g =R i , the intermediate variable v g =v i , δ is the moving step length;

[0051] Step 5: According to the model order K, R g Perform truncated eigenvalue decomposition (i.e. retain the first K largest eigenvalues ​​and the eigenvectors corresponding to the eigenvalues ​​after eigenvalue decomposition), that is (Σ K ,U K )=eigs(R g ,K), eigs(R g ,K) indicates R g Perform truncated eigenvalue decomposition, U K is the matrix composed of eigenvectors obtained by truncated eigenvalue decomposition, Σ K is the vector of eigenvalues ​​obtained by truncated eigenvalue decomposition;

[0052] According to Σ K and U K Calculate the intermediate variable matrix According to the matrix Calculating the Matrix And according to the matrix Z g Get the matrix

[0053] Step 6: According to R i+1 , R i , Z i+1 and Z i Calculate the combined residual

[0054] If i=0, set i=i+1 and return to step 4;

[0055] If i>0, execute step 7;

[0056] Step 7. Compare (NK)λ(K) with η (set according to actual conditions), where λ(K) is the Kth eigenvalue obtained by truncated eigenvalue decomposition (also the Kth eigenvalue for R g Perform eigenvalue decomposition to obtain the Kth largest eigenvalue);

[0057] If (NK)λ(K)≥η, execute step eight;

[0058] If (NK)λ(K)<η, then execute step nine;

[0059] Step 8. Compare and Size:

[0060] like Then continue to judge whether it satisfies If satisfied The iteration ends and Z i+1 As the line spectrum estimation result The model order K at the last iteration is used as the signal model order estimation result like Then let i=i+1, and return to step 4;

[0061] like Then let K=K+1, i=i+1, and return to step 4;

[0062] Step 9: Compare With the threshold ε tol Size:

[0063] like The iteration ends and Z i+1 As the line spectrum estimation result The model order K at the last iteration is used as the signal model order estimation result

[0064] like Then let K=K-1, i=i+1, and return to execute step 4.

[0065] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the specific process of step 1 is as follows:

[0066] Step 1: Calculate the autocorrelation function of sample y (i.e., calculate the convolution of sample y with itself), and record the first, second, ..., Nth elements in the calculation result of the autocorrelation function as r(0), r(1), ..., r(N-1), respectively, where N is the signal length;

[0067] Step 2: Construct the autocorrelation matrix R based on the autocorrelation function calculation results yy :

[0068]

[0069] in, is plural, r * (1) is the transpose of r(1), r * (N-1) is the transpose of r(N-1), r * (N-2) is the transpose of r(N-2), and the autocorrelation matrix R yy The dimension is N×N.

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

[0071] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that the combined residual for:

[0072]

[0073] Among them, ||·|| F is the F-norm.

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

[0075] Specific implementation method 4: This implementation method is different from any one of the specific implementation methods 1 to 3 in that the auxiliary variable Z0 is:

[0076]

[0077] in, is the conjugate transpose of y0.

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

[0079] Specific implementation method 5: This implementation method is different from the specific implementation methods 1 to 4 in that the intermediate variable matrix

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

[0081] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the matrix Z g Get the matrix Z i+1 The method used is to g Method of projection onto semidefinite sets.

[0082] The other steps and parameters are the same as those in Specific Implementations 1 to 5.

[0083] The method of the present invention is compared with the original gridless line spectrum estimation method solved by SDP.

[0084] Step 1: Signal Generation

[0085] A sample y of length N is randomly generated; the model order is set to K; and the frequency interval Δf>1 / N.

[0086] Step 2: Evaluation Metrics

[0087] The evaluation indicators are normalized frequency mean square error, running time and success rate. The normalized frequency mean square error is defined as:

[0088]

[0089] When the estimated value of the signal model order is the same as the actual signal model order, it is considered successful. Multiple repeated experiments are performed under the same conditions, and the success rate is defined as the ratio of the number of successes to the number of repeated experiments.

[0090] Experiment 1: Let the signal length N = 15, the model order K = 5, repeat the experiment 100 times at different signal-to-noise ratios, and compare the mean square error and success rate of the normalized frequency.

[0091] The comparison of the normalized frequency mean square error between the original semidefinite programming method, i.e., the interior point method, and the method of the present invention is as follows: Figure 1 It can be seen that under different signal-to-noise ratio conditions, the normalized frequency mean square error of the method of the present invention and the original method is almost the same, which indicates that the method of the present invention has almost no effect on the solution accuracy.

[0092] The comparison of the success rate of model order estimation between the original semidefinite programming method and the method of the present invention is as follows: Figure 2 As shown. Figure 2 It can be seen that the success rate of model order estimation of the method of the present invention is slightly higher than that of the original method.

[0093] Experiment 2: Let the model order K = 5, the signal-to-noise ratio SNR = 10dB, repeat the experiment 100 times, and compare the running time required for different signal lengths N, such as Figure 3 shown.

[0094] It can be seen from this experiment that compared with the original semidefinite programming method, the method of the present invention significantly reduces the computational complexity of gridless line spectrum estimation. Combining the two experiments for analysis, the method of the present invention solves the problem that the existing alternating projection gridless method requires prior information, and its success rate of model order estimation is slightly higher than that of the original semidefinite programming method. At the same time, the method of the present invention has made outstanding contributions in reducing computational complexity.

[0095] 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 gridless line spectrum estimation method for unknown model order based on alternating projection, characterized in that: The method specifically comprises the following steps: Step 1: Construct the sample autocorrelation matrix R based on sample y yy ; Step 2: Set the intermediate variable R0 = R yy , intermediate variable y0=y, intermediate variable v0=tr(R0) / N, and define auxiliary variable Z0 according to R0, y0 and v0, tr(R0) is the trace of R0; Step 3: Initialize the model order K = 1, and initialize the number of iterations i = 0; Step 4: Calculate the intermediate variable y g =y i -δ(y i -y0), intermediate variable R g =R i , the intermediate variable v g =v i , δ is the moving step length; Step 5: According to the model order K, R g Perform truncated eigenvalue decomposition, that is, (Σ K ,U K )=eigs(R g ,K),U K is the matrix composed of eigenvectors obtained by truncated eigenvalue decomposition, Σ K is the vector of eigenvalues ​​obtained by truncated eigenvalue decomposition; According to Σ K and U K Calculate the intermediate variable matrix According to the matrix Calculating the Matrix And according to the matrix Z g Get the matrix Step 6: According to R i+1 , R i , Z i+1 and Z i Calculate the combined residual If i=0, set i=i+1 and return to step 4; If i>0, execute step 7; Step 7: Compare (NK)λ(K) and η, where λ(K) is the Kth eigenvalue obtained by truncated eigenvalue decomposition; If (NK)λ(K)≥η, execute step eight; If (NK)λ(K)<η, then execute step nine; Step 8. Compare and Size: like Then continue to judge whether it satisfies If satisfied The iteration ends and Z i+1 As the line spectrum estimation result The model order K at the last iteration is used as the signal model order estimation result like Then let i=i+1, and return to step 4; like Then let K=K+1, i=i+1, and return to step 4; Step 9: Compare With the threshold ε tol Size: like The iteration ends and Z i+1 As the line spectrum estimation result The model order K at the last iteration is used as the signal model order estimation result like Then let K=K-1, i=i+1, and return to execute step 4.

2. According to claim 1, a gridless line spectrum estimation method with unknown model order based on alternating projection is characterized in that: The specific process of step one is: Step 1: Calculate the autocorrelation function of sample y, and record the first, second, ..., Nth elements in the calculation result of the autocorrelation function as r(0), r(1), ..., r(N-1), respectively, where N is the signal length; Step 2: Construct the autocorrelation matrix R based on the autocorrelation function calculation results yy : in, is plural, r * (1) is the transpose of r(1), r * (N-1) is the transpose of r(N-1), r * (N-2) is the transpose of r(N-2), and the autocorrelation matrix R yy The dimension is N×N.

3. The gridless line spectrum estimation method with unknown model order based on alternating projection according to claim 2, characterized in that: The combined residual for: Among them, ||·|| F is the F-norm.

4. The gridless line spectrum estimation method for unknown model order based on alternating projection according to claim 3, characterized in that: The auxiliary variable Z0 is: in, is the conjugate transpose of y0.

5. The gridless line spectrum estimation method with unknown model order based on alternating projection according to claim 4, characterized in that:The intermediate variable matrix 6. The gridless line spectrum estimation method with unknown model order based on alternating projection according to claim 5, characterized in that: According to the matrix Z g Get the matrix Z i+1 The method used is to g Method of projection onto semidefinite sets.

Images