A beamforming interference suppression method based on partial correlation orthogonal matching pursuit
By constructing a pure spatial spectrum peak polynomial through a partial correlation orthogonal matched tracking beamforming method, and optimizing array errors, the problem of reduced suppression capability under the conditions of desired signal self-cancellation and joint error is solved, thus achieving efficient interference suppression and desired signal enhancement.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- HARBIN INST OF TECH
- Filing Date
- 2026-04-20
- Publication Date
- 2026-07-14
AI Technical Summary
Existing beamforming interference suppression methods have reduced suppression capabilities under scenarios involving desired signal self-cancellation and joint errors, and sparse signal recovery methods have high computational overhead, limiting practical deployment and application.
A beamforming method based on partial correlation orthogonal matched tracking is adopted. By constructing a pure spatial spectrum peak polynomial, iteratively estimating the angle of arrival of the incident signal, optimizing the array amplitude and phase error and position error, generating beamforming weights, and avoiding the self-cancellation effect of the desired signal.
It improves interference suppression capability, enhances the reception of desired signals, has stronger robustness and lower computational complexity, and is suitable for communication in non-ideal environments.
Smart Images

Figure CN122394616A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of array signal processing technology, specifically relating to a beamforming interference suppression method based on partial correlation orthogonal matched tracking. Background Technology
[0002] Adaptive beamforming, as an important research branch in the field of array signal processing, can filter out various interferences in the spatial domain based on antenna arrays and provide beam gain for the desired signal. It has been widely used in wireless communication, radar, sonar, medical imaging and other fields.
[0003] While it has been theoretically proven that minimum variance distortionless response beamformers can maximize the output signal-to-interference-plus-noise ratio (SNR) after interference suppression, these beamformers are highly susceptible to unavoidable non-ideal factors in real-world operating environments, such as limited snapshots, beam pointing errors, array amplitude and phase errors, and positional errors, leading to a sharp decline in output signal quality. To address this, robust beamforming algorithms based on diagonal loading, feature space, uncertainty sets, and covariance matrix reconstruction have emerged. However, the performance degradation of beamformers caused by the self-cancellation of desired signals has not been completely resolved, and most robust measures only address single non-ideal factors, which are significantly affected in scenarios with combined errors. In recent years, sparse signal recovery methods have gradually been applied in beamforming, demonstrating unique advantages in overcoming the self-cancellation effect of desired signals; however, the high computational cost and accurate dictionary matrix required for sparse signal recovery limit its practical deployment and application in mismatch scenarios. Summary of the Invention
[0004] The purpose of this invention is to address the problem that the suppression capability of existing beamforming interference suppression methods is reduced due to the self-cancellation of the desired signal and the influence of joint errors. Therefore, this invention proposes a beamforming interference suppression method based on partial correlation orthogonal matched tracking.
[0005] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:
[0006] According to one aspect of the present invention, a beamforming interference suppression method based on partial correlation orthogonal matched tracking is provided, the method specifically comprising the following steps:
[0007] Step 1: Use N receiving antennas to receive analog signals from free space. The analog signals received by each receiving antenna are processed by their own antenna radio frequency processing links to form N digital sequences.
[0008] Step 2: Integrate the N digital sequences formed in Step 1. Specifically, summarize the data from the N receiving antennas at L consecutive sampling points into an N-row, L-column received data matrix X.
[0009] Step 3: Based on the partial correlation orthogonal matching pursuit criterion and combined with the polynomial root-finding method, estimate the angle of arrival of the incident signal in space to obtain the estimated angle vector Θ. f ;
[0010] Step 4, if Θ f If the number of elements in the matrix is greater than or equal to 1, proceed to step five; otherwise, let the projection matrix P under the joint error be... Ae Given an N x N identity matrix, an N x N amplitude-phase matrix G, and a 1 x N zero vector p, continue with step seventeen F.
[0011] Step 5: Initialize loop variable u = 1, initialize loop variable c = 1;
[0012] Step 6: Initialize the amplitude and phase matrix G as an N-row, N-column identity matrix and the position error vector p as a 1-row, N-column vector of all zeros;
[0013] Set power threshold t p = N / C, where Let X be the average power, C be a preset constant, and set the iteration stopping threshold ε3;
[0014] Step 7: For the estimated angle vector Θ f The angles in the middle are precisely estimated and redundancy is removed;
[0015] Step 8: Determine whether to exit the loop controlled by loop variable c; specifically:
[0016] Step 8A: Determine the vector Θ fc With vector Θ f Do they have the same number of elements?
[0017] If they are different, do not exit the loop and proceed to step nine; otherwise, continue to step eightB.
[0018] Step 8B: Let vector Θ fc With vector Θ f Let d be the distance vector between them, and let Θ be the k-th element of the distance vector d. fc The k-th element and Θ f The absolute value of the difference between the k-th elements; k = 1, 2, ..., card(Θ) f );
[0019] Step 8C: Determine whether the mean of all elements in the distance vector d is less than the threshold ε3;
[0020] If yes, then exit the loop controlled by loop variable c and continue to step ten; otherwise, do not exit the loop controlled by loop variable c and continue to step nine.
[0021] Step 9: Determine if the loop variable c is equal to the set iteration number iter2;
[0022] If yes, continue to step ten; otherwise, set c = c + 1 and return to step seven.
[0023] Step 10: Initialize the relative power vector as P r Vector P r The number of elements and the current angle vector Θ f The number of elements is the same, and vector P r All elements are 0;
[0024] Step 11: Calculate the relative power vector P r And according to the relative power vector P r Adjust the current angle vector Θ f The element order in the vector is used to obtain the angle vector Θ after the element order is adjusted. f ;
[0025] Step 12: Based on the relative power vector P r and the angle vector Θ after element order adjustment f Reconstruct the original transmitted signal S;
[0026] Step 12A: Construct a copy of the estimated angle vector Θ fc equal to the vector Θ after element order adjustment f Construct a relative power vector replica P rc Equal to the relative power vector P r ;
[0027] Step 12B, if Θ fc If the set is empty, then stop reconstructing the original transmitted signal and let the amplitude-phase matrix copy G... e p is an N x N identity matrix and a copy of the position error vector. e Given a vector of 1 row and N columns containing only zeros, proceed to step fourteen.
[0028] Otherwise, continue reconstructing the original transmitted signal, retaining a copy of the currently estimated angle vector Θ. fc The first N k N elements, where N k N is a positive integer. k Satisfy: P rc The first N k The sum of the elements and P rc The ratio of the sum of all elements of P is greater than or equal to ε4, and P rc The first N k – The sum of one element and P rc The ratio of the sum of all its elements is less than ε4;
[0029] Step 12C: Construct the reconstructed angle steering vector matrix A under joint error. r ;
[0030] Step 12D: Reconstruct the original transmitted signal S. The reconstruction method is as follows:
[0031]
[0032] Step 13: Update the amplitude and phase matrix copy G e and position error vector copy p e The specific update method is as follows:
[0033] Step 13A: Initialize the amplitude-phase matrix replica G e p is an N x N identity matrix and a copy of the position error vector. e Create a 1-row, N-column vector of all zeros; and initialize the loop variable q = 1;
[0034] Step 13B: Let x be the q-th row of X;
[0035] Step 13C: Calculate vector w r :
[0036]
[0037] Step 13D: Calculate vector u r u r The kth element is w r The kth element and The product of θ, where θ r, k Represents vector Θ fc The k-th element, k = 1, 2, ..., card(Θ) fc );
[0038] Step 13E: Take vector u respectively r The phase of each element in the vector is used to obtain the vector a. r ;
[0039] Step 13F: Transfer matrix G e The q-th diagonal element is updated to vector u r The first N in k The mean of the moduli of each element;
[0040] Step 13 G, if N k If N is greater than 1, then proceed to step thirteen H. k If the value is 1, then proceed to step thirteen (I).
[0041] Step 13H, for vector p e Matrix G e Update:
[0042] Let Θ fc The sine value of the midpoint angle is the independent variable, a. r The first N k Performing linear regression with each element as the dependent variable yields the slope k. r and intercept b r and vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ;
[0043] Then proceed to step thirteen (J);
[0044] Step Thirteen I: Preserve vector p e Let b remain unchanged. r Let vector a r The first element will be the matrix G. e The q-th diagonal element is multiplied by the original value. Then proceed to step thirteen (J).
[0045] Step 13J: Determine if q is equal to N;
[0046] If yes, proceed to step thirteen K; otherwise, set q = q + 1 and return to step thirteen B.
[0047] Step 13 K, Obtain the vector |p e The largest element p in | max Calculate the vector |diag(G) e The largest element G in () – 1| max Where, diag(·) represents the extraction matrix G e The main diagonal elements are used to form a column vector, and |·| represents the vector diag(G) e Take the absolute value of each element in 1.
[0048] Step 13 L, if p max If it is greater than the preset threshold ε5, then let p e If it is a vector of all zeros, otherwise keep p. e Unchanged; if G max If it is greater than the preset threshold ε6, then let G e It is an N x N identity matrix, otherwise G e Remain unchanged;
[0049] Then proceed to step fourteen;
[0050] Step 14: Let vector Δp = |p – p e|, vector Δg = |diag(G e ) – diag(G)|;
[0051] Step 15: Determine whether to exit the loop controlled by loop variable u, specifically:
[0052] Step 15A: Update p = p e G = G e ;
[0053] Step 15B: If the mean of all elements in vector Δp is greater than or equal to ε3 or the mean of all elements in vector Δg is greater than or equal to ε3, then do not exit the loop and proceed to step 16; otherwise, exit the loop and continue to step 17.
[0054] Step 16: Determine if the loop variable u is equal to the preset maximum number of iterations iter3;
[0055] If yes, continue to step seventeen; otherwise, set u = u + 1 and c = 1, and return to step seven.
[0056] Step 17: Within the desired signal range Θ s The search for the desired signal angle involves the following steps:
[0057] Step 17A: If the angle vector Θ has been estimated f If the number of elements in the vector is greater than or equal to 1, then from the estimated angle vector Θ f Delete the region Θ located in the desired signal range s The elements in the vector form the residual angle vector Θ. d Then proceed to step seventeen-b;
[0058] If the angle vector Θ has been estimated f If the number of elements in the matrix is 0, then let the projection matrix P under the joint error be... Ae Given an N x N identity matrix, proceed to step seventeen (F).
[0059] Step 17B: Based on the remaining angle vector Θ d Update the projection matrix P under joint error Ae And obtain the candidate angle set Θ ce ;
[0060] Step 17C: In the candidate angle set Θ ce Only the region containing the desired signal Θ is retained. s The angle in the set is obtained by removing the remaining elements, resulting in the updated set of candidate angles Θ. ce ;
[0061] Step 17.D. If the updated candidate angle set Θ ceIf the number of elements located within the desired signal interval is greater than or equal to 1, then proceed to step 17E; otherwise, proceed to step 17F.
[0062] Step 17E: Based on the updated candidate angle set Θ ce Calculate the estimated angle of arrival θ of the desired signal s =θ opt2 Then proceed to step eighteen;
[0063] Step 17F: Let the estimated angle of arrival of the desired signal be θ. s For the desired signal interval Θ s The midpoint, then proceed to step eighteen;
[0064] Step 18: Calculate the beamforming weights. The specific calculation method is as follows:
[0065]
[0066] Among them, P b (θ s Let be an N x N diagonal matrix, and let the nth diagonal element be 0. p n Let n be the nth element of vector p, where n = 1, 2, ..., N;
[0067] Step 19: Output the beamformed signal s to complete interference suppression and desired signal reception.
[0068] .
[0069] According to another aspect of the present invention, a beamforming interference suppression method based on partial correlation orthogonal matched tracking is provided, the method specifically including the following steps:
[0070] Step 1: Use N receiving antennas to receive analog signals from free space. The analog signals received by each receiving antenna are processed by their own antenna radio frequency processing links to form N digital sequences.
[0071] Step 2: Integrate the N digital sequences formed in Step 1. Specifically, summarize the data from the N receiving antennas at L consecutive sampling points into an N-row, L-column received data matrix X.
[0072] Step 3: Based on the partial correlation orthogonal matching pursuit criterion and combined with the polynomial root-finding method, estimate the angle of arrival of the incident signal in space to obtain the estimated angle vector Θ. f ;
[0073] Step 4, if Θ f If the number of elements in matrix P is greater than or equal to 1, then proceed to steps five through twelve; otherwise, let matrix P... AdGiven an N x N identity matrix, proceed directly to step 10F;
[0074] Step 5: Initialize the loop variable c = 1;
[0075] Step 6: Set the power threshold t p = N / C, where Let X be the average power, C be a preset constant, and set the iteration stopping threshold ε3;
[0076] Step 7: Perform precise estimation and redundancy removal on the angle of arrival of the incident signal; specifically:
[0077] Step 7A: Set the current vector Θ f The number of elements in the array is denoted as N. f = card(Θ f Let the already estimated angle vector copy Θ fc Equal to Θ f ;
[0078] Step 7B: Initialize the loop variable d = 1, and the current refinement angle index index = 1;
[0079] Step 7C: Let the angle to be refined be θ p For Θ f Let the index of the element be the remaining angle vector Θ. d For Θ f Remove θ p The vector after;
[0080] Step 7D: Construct the remaining angle guidance vector matrix A d ;
[0081] Step 7E: Constructing and maintaining the remaining angle guidance vector matrix A d The corresponding projection matrix P Ad ;
[0082] Step 7F: Based on matrix P Ad Generate the coefficient vector v of a 2N-2 degree polynomial. pd ;
[0083] Step 7 G, based on matrix B d Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. qd ;
[0084] Step 7H, based on matrix B d Generate the coefficient vector v of a 2N-2 degree polynomial. rd ;
[0085] Step 7.1: Based on matrix P AdΓ1 generates the coefficient vector v of a 2N-2 degree polynomial. sd ;
[0086] Step 7J, using v pd and v qd Convolution minus v rd and v sd The convolution yields the coefficient vector v of the final 4N – 4th degree pure spatial spectral peak polynomial. cd ;
[0087] Step 7K: Transfer vector v cd The nth element is considered as the coefficient of the 4N – 3 – n power in the peak polynomial of the pure space spectrum. The roots of the peak polynomial of the pure space spectrum are solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1, and the absolute values of each subtraction are taken. Roots whose absolute values are greater than the threshold ε1 are then removed. The roots of the remaining polynomials are used to form a set S. Rd , n = 1, 2,..., 4N – 3;
[0088] Step 7 L, for set S Rd Each element in the algorithm is processed individually to obtain the candidate angle set Θ. cd ;
[0089] Step 7 M: Delete the candidate angle set Θ cd The redundant angles in the equation are used to obtain the updated candidate angle set Θ. cd ;
[0090] Step 7 N: Update the candidate angle set Θ cd Calculate the function value f1(θ) for each angle in the equation, and output the maximum function value f. max2 and the maximum function value f max2 The corresponding angle θ opt2 ;
[0091] Step 7 O, f max2 Is it less than or equal to the power threshold t? p ;
[0092] If so, then delete Θ. f The index-th element forms the new estimated angle vector Θ. f ;
[0093] Otherwise, vector Θ f Replace the index-th element with θ opt2 And let index = index + 1;
[0094] Step 7 P: Determine if d equals N f ;
[0095] If not, let d = d + 1, and return to step seven C;
[0096] If so, then output the final angle vector Θ. f And continue with step eight;
[0097] Step 8: Determine whether to exit the loop controlled by loop variable c; specifically:
[0098] Step 8A: Determine the vector Θ fc With vector Θ f Do they have the same number of elements?
[0099] If they are different, do not exit the loop and proceed to step nine; otherwise, continue to step eightB.
[0100] Step 8B: Let vector Θ fc With vector Θ f Let d be the distance vector between them, and let Θ be the k-th element of the distance vector d. fc The k-th element and Θ f The absolute value of the difference between the k-th elements; k = 1, 2, ..., card(Θ) f );
[0101] Step 8C: Determine whether the mean of all elements in the distance vector d is less than the threshold ε3;
[0102] If yes, then exit the loop controlled by loop variable c and continue to step ten; otherwise, do not exit the loop controlled by loop variable c and continue to step nine.
[0103] Step 9: Determine if the loop variable c is equal to the set iteration number iter2;
[0104] If yes, continue to step ten; otherwise, set c = c + 1 and return to step seven.
[0105] Step 10: Within the desired signal range Θ s The search for the desired signal angle involves the following steps:
[0106] Step 10A: From the estimated angle vector Θ f Delete the region Θ located in the desired signal range s The elements in the vector form the residual angle vector Θ. d ;
[0107] Step 10B: Repeat steps 7D to 7L to obtain the candidate angle set Θ. cd And update matrix P Ad ;
[0108] Step 10C: In the candidate angle set Θ cd Only the region containing the desired signal Θ is retained. s From the given angles, remove the remaining elements to obtain the updated candidate angle set Θ. cd ;
[0109] Step 10D: If the updated candidate angle set Θ cd If the number of elements in the array is greater than or equal to 1, proceed to step 10E; otherwise, proceed to step 10F.
[0110] Step 10E: Execute step 7N to obtain the estimated angle of arrival (θ) of the desired signal. s = θ opt2 Then proceed to step eleven;
[0111] Step 10F: Let the estimated angle of arrival of the desired signal be θ. s For the desired signal interval Θ s At the midpoint, proceed to step eleven;
[0112] Step 11: Calculate beamforming weights:
[0113]
[0114] Step 12: Output the beamformed signal s to complete interference suppression and desired signal reception.
[0115] .
[0116] The beneficial effects of this invention are:
[0117] This invention first constructs a pure spatial spectrum peak polynomial and selects the solution that maximizes the pure spatial spectrum value as the angle of arrival estimate. Second, iteratively obtains the angle of arrival for all incident signals and continuously refines its values. Then, the original transmitted signal, array amplitude and phase errors, and array position errors are alternately optimized using the received signal and the estimated angles of arrival for each incident signal. Finally, beamforming weights are generated based on the estimated angles of arrival, array amplitude and phase errors, and array position errors, thus avoiding the self-cancellation effect of the desired signal and achieving interference suppression and desired signal enhancement. This invention incorporates the correlation between array steering vectors, possessing a stronger ability to distinguish between neighboring interference and related interference than the original orthogonal matched pursuit algorithm, improving interference suppression capabilities. It also achieves lower computational complexity while maintaining performance comparable to sparse Bayesian algorithms. Furthermore, this invention considers the simultaneous existence of array amplitude and phase errors and position errors, exhibiting stronger robustness in joint error scenarios and providing support for reliable communication in non-ideal environments. Attached Figure Description
[0118] Figure 1This is a flowchart of a beamforming interference suppression method based on partial correlation orthogonal matched tracking according to the present invention;
[0119] Figure 2 This is a simulation performance diagram of a beamforming interference suppression method based on partial correlation orthogonal matched tracking according to a specific embodiment of the present invention.
[0120] The interference signal has two angles of arrival (ACOs), uniformly distributed within the ranges of [-45°, -35°] and [15°, 25°], respectively. The desired signal ACO is uniformly distributed within the range of [-20°, -10°]. The array amplitude error is uniformly distributed within the range of [-0.5dB, 0.5dB]. The array phase error is uniformly distributed within the range of [-5°, 5°]. The array position error is uniformly distributed within the range of [-0.1dB, 0.5dB]. , 0.1 Evenly distributed within the range, The nominal array spacing, λ represents the wavelength of the electromagnetic wave used in communication;
[0121] Figure 3 This is a simulation performance diagram of a beamforming interference suppression method based on partial correlation orthogonal matched tracking according to a specific embodiment of the present invention.
[0122] The interference consists of two elements, with interference angles of arrival uniformly distributed within the ranges of [-45°, -35°] and [15°, 25°], respectively; the desired signal angle of arrival is uniformly distributed within the range of [-20°, -10°]; the array amplitude error is uniformly distributed within the range of [-0.5dB, 0.5dB]; the array phase error is uniformly distributed within the range of [-5°, 5°]; and the array position error is negligible.
[0123] Figure 4 This is a simulation performance diagram of a beamforming interference suppression method based on partial correlation orthogonal matched tracking according to a fifth specific embodiment of the present invention.
[0124] The number of interferences is 2, and the angles of arrival of the interferences are uniformly distributed within the ranges of [-45°, -35°] and [15°, 25°], respectively; the angle of arrival of the desired signal is uniformly distributed within the range of [-20°, -10°]; the array position error is within... Uniformly distributed within the range; array amplitude and phase errors are negligible;
[0125] Figure 5 The simulation performance diagram of the beamforming interference suppression method based on partial correlation orthogonal matched tracking of the present invention is shown in the case that the array amplitude and phase errors and position errors can be ignored in the sixth specific implementation method.
[0126] The number of interferences is 2, and the angles of arrival of the interferences are uniformly distributed in the ranges of [-45°, -35°] and [15°, 25°], respectively; the angle of arrival of the desired signal is uniformly distributed in the range of [-20°, -10°]. Detailed Implementation
[0127] This invention provides a multi-antenna receiver equipped with N = 8 receiving antennas, all uniformly arranged in a straight line with a spacing of λ / 2, and each antenna having its own dedicated RF chain. For cases where the sum of the number of interference signals and the desired signal number is not greater than the number of receiving antennas, this invention proposes a beamforming interference suppression method, which is described in detail below.
[0128] Specific Implementation Method 1: Combination Figure 1 This embodiment describes a beamforming interference suppression method based on partial correlation orthogonal matched tracking. When both array amplitude and phase errors and position errors are non-negligible, the method specifically includes the following steps:
[0129] Step 1: Receive analog signals from free space using N receiving antennas. Each receiving antenna receives an analog signal containing an incident angle located in the interval Θ. s = [ - 5°, The expected signal in [+ 5°] and the interference signals of unknown number, unknown incident angle range and unknown power are processed by the analog signals received by each receiving antenna through its own antenna radio frequency processing link to form N digital sequences.
[0130] in, The angle of arrival of the a priori expected signal shall not exceed 5° from the actual expected signal angle of arrival.
[0131] Step 2: Integrate the N digital sequences formed in Step 1. Specifically, summarize the data from the N receiving antennas at L consecutive sampling points into an N-row, L-column received data matrix X (in this invention, L is 50).
[0132] Step 3: Based on the partial correlation orthogonal matching pursuit criterion and combined with the polynomial root-finding method, estimate the angle of arrival of possible incident signals in space, and obtain the estimated angle vector Θ. f Specifically:
[0133] Step 3A, Initialization: The estimated angle vector Θ f An empty vector, an estimated angle steering vector matrix A is an empty matrix, and a projection matrix P A Given an N x N identity matrix, with loop variable i = 1;
[0134] Step 3B: Based on matrix P AGenerate the coefficient vector v of a 2N-2 degree polynomial. p For example, when the highest degree term of the polynomial is 6, the polynomial takes the form of: , , … The coefficient vectors v are respectively p The element value in;
[0135] Step 3C: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix BΓ1. q ;
[0136] Where, the intermediate variable matrix B = P A XX H P A The superscript "H" indicates the conjugate transpose. Matrix Γ1 is an N-row N-column diagonal matrix, and the nth diagonal element of the diagonal matrix Γ1 is n - 1, where n = 1, 2, ..., N.
[0137] Step 3D: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix B. r ;
[0138] Step 3E: Based on matrix P A Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. s ;
[0139] Step 3F: Calculate the coefficient vector v p and coefficient vector v q Convolution, then calculate the coefficient vector v r and coefficient vector v s The convolution, using the coefficient vector v p and coefficient vector v q The convolution result minus the coefficient vector v r and coefficient vector v s The convolution result yields the coefficient vector v of the 4N – 4th order pure spatial spectral peak polynomial. c ;
[0140] In this invention, the peak polynomial of the pure space spectrum is defined as:
[0141]
[0142] in, Let a(z) = [1, z, z] represent the independent variable of the peak polynomial of the pure space spectrum. 2 , …, z N-1 ] T ;
[0143] Step 3 G: Convert the coefficient vector vc The nth element in the equation is considered as the coefficient of the peak polynomial of the pure space spectrum to the power of 4N – 3 – n. Then, each root of the peak polynomial of the pure space spectrum is solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1 to obtain the subtraction results.
[0144] Next, calculate the absolute value of each difference result, and remove the roots corresponding to the absolute values greater than the threshold ε1 (in this invention, the value of ε1 is 0.7). Use the remaining roots to form the set S of polynomial roots. R ;
[0145] Where n = 1, 2, ..., 4N – 3;
[0146] Step 3H, for set S R Each element in the algorithm is processed individually to obtain the candidate angle set Θ. c ;
[0147] Step 3I, if Θ f If the set is empty, proceed to step three K; otherwise, proceed to step three J.
[0148] Step 3J: Delete the candidate angle set Θ c The redundant angles in the calculation are as follows:
[0149] Step 3 J1: Initialize the loop variable b = 1;
[0150] Step 3 J2, for vector Θ f The b-th element in Θ c Find the four elements closest to it, and move the four found elements from Θ. c Delete them to form a new candidate angle set Θ c ;
[0151] Step 3 J3, if b equals vector Θ f If the number of elements in the array is determined, proceed to step three K; otherwise, set b = b + 1 and return to step three J2.
[0152] Step 3K: Construct a function f1(θ) with respect to the angle:
[0153]
[0154] Where f1(θ) is the arrival angle quality metric based on the partial correlation coefficient, also known as the pure space spectrum expression; a(θ) represents the steering vector of the receiving array for angle θ;
[0155]
[0156] Where θ represents the angle. The base of the natural logarithm. The superscript T indicates the imaginary unit;
[0157] Step 3 L: For the candidate angle set Θ c Calculate the function value f1(θ) for each angle in the equation, and output the calculated maximum function value f. max1 Minimum function value f min1 and the maximum function value f max1 The corresponding angle θ opt1 ;
[0158] Step 3 M, if f max1 / f min1 If the value is greater than or equal to ε2 (ε2 is 2 in this invention), then proceed to step three N; otherwise, jump directly to step three Q.
[0159] Step 3 N, θ opt1 Add vector Θ f Update the estimated angle steering vector matrix A, where the k-th column of A is Θ. f The guiding vector of the kth element, k = 1, 2, ..., card(Θ) f ), where card(Θ) f ) represents the vector Θ f The number of elements in;
[0160] Step 3: Update the projection matrix P A The specific update formula is as follows:
[0161]
[0162] Among them, I N Represents an N x N identity matrix, where the superscript -1 indicates the inverse of the matrix;
[0163] Step 3P: If i equals N, then execute Step 3Q; otherwise, let i = i + 1 and return to Step 3B.
[0164] Step 3 Q: Output the estimated angle vector Θ f ;
[0165] Step 4, if Θ f If the number of elements in the matrix is greater than or equal to 1, proceed to step five; otherwise, let the projection matrix P under the joint error be... Ae Given an N x N identity matrix, an N x N amplitude-phase matrix G, and a 1 x N zero vector p, continue with step seventeen F.
[0166] Step 5: Initialize loop variable u = 1, initialize loop variable c = 1;
[0167] Step 6: Initialize the amplitude and phase matrix G as an N-row, N-column identity matrix and the position error vector p as a 1-row, N-column vector of all zeros;
[0168] Set power threshold t p = N / C, where Let X be the average power, C be a preset constant (in this invention, C is 50), and set the iteration stopping threshold ε3 (in this invention, threshold ε3 is 10). -6 );
[0169] Step 7: For the estimated angle vector Θ f The angle in the middle is precisely estimated and redundancy is removed; specifically:
[0170] Step 7A: Set the current vector Θ f The number of elements in the array is denoted as N. f = card(Θ f Let the already estimated angle vector copy Θ fc Equal to Θ f ;
[0171] Step 7B: Initialize the loop variable d = 1, and the current refinement angle index index = 1;
[0172] Step 7C: Let the angle to be refined be θ p For Θ f Let the index of the element be the remaining angle vector Θ. d For Θ f Remove θ p The vector after;
[0173] Step 7D: Construct the residual angle steering vector matrix A under the joint error. e Matrix A e The number of rows is N and the number of columns is card(Θ). d ), let Θ d The kth element is θ d,k Then A e The k-th column is represented as:
[0174]
[0175] Among them, P b (θ d,k Let P be an N x N diagonal matrix. b (θ d,k The nth diagonal element is p n Let n be the nth element of vector p, where n = 1, 2, ..., N, k = 1, 2, ..., card(Θ).d );
[0176] Step 7E: Construct the projection matrix P under joint error. Ae :
[0177]
[0178] Step 7F: Based on matrix G H P Ae G generates the coefficient vector v of a 2N-degree polynomial. pe ;
[0179] Step 7 G, based on matrix G H B e GΓ generates the coefficient vector v of a 2N-degree polynomial. qe ;
[0180] Wherein, matrix B e = P Ae XX H P Ae , matrix Γ = Γ1 + Diag(p), Diag(·) means to transform vector p into a diagonal matrix, and the nth diagonal element of the diagonal matrix Diag(p) is the nth element of vector p, n = 1, 2, ..., N;
[0181] Step 7H, based on matrix G H B e G generates the coefficient vector v of a 2N-degree polynomial. re ;
[0182] Step 7.1: Based on matrix G H P Ae GΓ generates the coefficient vector v of a 2N-degree polynomial. se ;
[0183] Step 7J, using v pe and v qe Convolution minus v re and v se The convolution yields the coefficient vector v of the final 4N – 4th degree pure spatial spectral peak polynomial. ce ;
[0184] Step 7K: Transfer vector v ce The nth element is considered as the coefficient of the 4N – 3 – n power in the peak polynomial of the pure space spectrum. The roots of the peak polynomial of the pure space spectrum are solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1, and the absolute values of each subtraction are taken. Roots whose absolute values are greater than the threshold ε1 are then removed. The roots of the remaining polynomials are used to form a set S.Re , n = 1, 2,..., 4N – 3;
[0185] Step 7 L, for set S Re Each element in the equation is processed separately (i.e., steps H1 to H3 are performed) to obtain the candidate angle set Θ. ce ;
[0186] Step 7 M: Delete the candidate angle set Θ ce The redundant angle (the specific steps are the same as steps three J1 to three J3, only Θ needs to be removed) f Replace with Θ d Θ c Replace with Θ ce ), thus obtaining the updated candidate angle set Θ ce ;
[0187] Step 7N: Construct the function f2(θ) of the angle when constructing the joint error:
[0188]
[0189] f2(θ) is the arrival angle quality measure based on the partial correlation coefficient when there is joint error, also known as the pure space spectrum expression when there is joint error;
[0190] Step 7.0: Update the candidate angle set Θ ce Calculate the function value of f2(θ) for each angle in the equation, and output the maximum function value f. max2 and the maximum function value f max2 The corresponding angle θ opt2 ;
[0191] Step 7 P, for angle θ opt2 The following corrections have been made:
[0192] Step 7 P1, based on matrix G H P Ae G generates the coefficient vector v of an N(N – 1) term power function. pl Specifically:
[0193] Step 7 P1a: Initialize the coefficient vector v pl It is a vector consisting of 1 row and N(N – 1) + 1 columns, all containing zeros;
[0194] Step 7 P1b: Initialize loop variable uu = 1, loop variable vv = 1, counter count = 0;
[0195] Step 7 P1c: If uu = vv, then v pl The N(N-1)+1th element is updated to v.pl The N(N-1)+1th element and matrix G H P Ae The sum of the values in row uu and column v of G is used; otherwise, count = count + 1, and v is set to... pl The count-th element is assigned to matrix G. H P Ae The value of the element in row uu and column vv of G;
[0196] Step 7 P1d: Determine if vv is equal to N. If yes, proceed to step 7 P1e; otherwise, set vv = vv + 1 and return to step 7 P1c.
[0197] Step 7 P1e: Determine if uu is equal to N;
[0198] If yes, proceed to step seven P1f; otherwise, set uu = uu + 1 and vv = 1, and return to proceed to step seven P1c.
[0199] Step 7 P1f: Output the coefficient vector v of the N(N – 1) term power function. pl ;
[0200] Step 7 P2, based on matrix G H B e GΓ generates the coefficient vector v of an N(N – 1) term power function. ql ;
[0201] v ql The generation method and v pl The same applies; only the matrix G in step seven, P1c, needs to be changed. H P Ae G is replaced with matrix G. H B e GΓ;
[0202] Step 7 P3, based on matrix G H B e G generates the coefficient vector v of an N(N – 1) term power function. rl ;
[0203] v rl The generation method and v pl The same applies; only the matrix G in step seven, P1c, needs to be changed. H P Ae G is replaced with matrix G. H B e G;
[0204] Step 7 P4, based on matrix G H P Ae GΓ generates the coefficient vector v of an N(N – 1) term power function.sl ;
[0205] v sl The generation method and v pl The same applies; only the matrix G in step seven, P1c, needs to be changed. H P Ae G is replaced with matrix G. H P Ae GΓ;
[0206] Step 7, P5: Construct the power vector t1 of the N(N – 1) term power function. The dimension of the power vector t1 is 1 row and N(N – 1) + 1 column. The specific construction method is as follows:
[0207] Step 7 P5a: Construct a position vector p with array position error. v Vector p v The nth element is the nth element of vector p plus n – 1, where n = 1, 2, ..., N;
[0208] Step 7 (P5b): Initialize the loop variable = 1, loop variable f = 1, counter count = 0;
[0209] Step 7 P5c, if If the value is f, then the N(N – 1) + 1th element of the exponentiation vector t1 is assigned the value 0; otherwise, let count = count + 1, and assign the countth element of the exponentiation vector t1 the value p. v The f-th element and p v The The difference between the elements;
[0210] Step 7 P5d: Determine if f is equal to N;
[0211] If yes, proceed to step seven P5e; otherwise, let f = f + 1 and return to proceed to step seven P5c.
[0212] Step 7 P5e, Judgment Is it equal to N?
[0213] If so, proceed to step seven (P5f); otherwise, let... = + 1, f = 1, and return to execute step seven P5c;
[0214] Step 7 P5f: Obtain the maximum value in the power vector t1, and add the obtained maximum value to each element in the power vector t1 to obtain the output power vector t1.
[0215] Step 7, P6: According to v pl and vql Construct N(N – 1)(N 2 – N + 3) / 2 term power function coefficient vector v b Specifically:
[0216] Step 7 (P6a): Initialize the coefficient vector v b For a single row N(N – 1)(N 2 – A vector of all zeros consisting of N + 3) / 2 + 1 columns;
[0217] Step 7 (P6b): Initialize loop variable g = 1; Initialize loop variable h = g; Initialize counter count = 0;
[0218] Step 7 (P6c): Let count = count + 1;
[0219] Step 7 P6d: If g = h, then let v b The count-th element is v pl The g-th element and v ql The product of h elements; otherwise, let v b The count-th element is v pl The g-th element and v ql The product of the h-th elements plus v pl The h-th element and v ql The product of the g-th elements, i.e., v b The count-th element is: v pl (g)×v ql (h)+
[0220] v pl (h)×v ql (g), v pl (g) represents v pl The g-th element, v ql (h) represents v ql The h-th element, v pl (h) represents v pl The h-th element, v ql (g) represents v ql The g-th element;
[0221] Step 7 P6e: Determine if h is equal to N(N – 1) + 1;
[0222] If yes, proceed to step seven P6f; otherwise, set h = h + 1 and return to step seven P6c.
[0223] Step 7 P6f: Determine if g is equal to N(N – 1) + 1. If yes, execute step 7 P6g; otherwise, let g = g + 1 and h = g, and return to execute step 7 P6c.
[0224] Step 7 P6g, output according to v pl and v ql Construct N(N – 1)(N) 2 – N + 3) / 2 term power function coefficient vector v b ;
[0225] Step 7 (P7): According to v rl and v sl Construct N(N – 1)(N 2 – N + 3) / 2 term power function coefficient vector v a ;
[0226] The specific construction method is the same as in step seven, P6, except that v in step seven, P6d needs to be changed. pl Replace with v rl 1. In step seven, P6d, v ql Replace with v sl ;
[0227] Step 7, P8: Construct N(N-1)(N) based on the power vector t1 of the N(N-1) term power function. 2 – N + 3) / 2 power function power vector t2, power vector t2 has a dimension of 1 row N(N – 1)(N 2 – N + 3) / 2 + 1 columns, the specific construction method is as follows:
[0228] Step 7 (P8a): Initialize loop variable l = 1; Initialize loop variable m = l; Initialize counter count = 0;
[0229] Step 7 (P8b): Let count = count + 1;
[0230] Step 7 P8c: Assign the count-th element of t2 to the sum of the l-th element and the m-th element of t1;
[0231] Step 7 P8d: Determine if m is equal to N(N – 1) + 1;
[0232] If yes, proceed to step seven P8e; otherwise, set m = m + 1 and return to step seven P8b.
[0233] Step 7 P8e: Determine if l is equal to N(N – 1) + 1;
[0234] If yes, proceed to step seven P8f; otherwise, let l = l + 1, m = l, and return to proceed to step seven P8b.
[0235] Step 7 P8f, Output N(N – 1)(N 2 – N + 3) / 2 power vectors t2;
[0236] Step 7, P9: According to v b v a The peak power function f(x) of the pure spatial spectrum under joint error and the derivative function g(x) of f(x) are constructed as follows:
[0237]
[0238] Where x represents the independent variable of the function, Indicates v a The nth element, Indicates v b Let t2(n) represent the nth element of t2, where n = 1, 2, ..., N(N – 1)(N 2 – N + 3) / 2 + 1; N’ = N(N – 1)(N 2 – N + 3) / 2 + 1;
[0239] Step 7, page 10: Use Newton's method to iteratively find the maximum points of the peak power function f(x) of the pure space spectrum. The specific iterative steps are as follows:
[0240] Step 7 (P10a): Initialize the iteration variable r = 1, and the initial iteration value x. (0) = ;
[0241] Step 7 P10b, x (r) = x (r – 1) – f[x (r – 1) ] / g[x (r – 1) ];
[0242] Step 7 P10c: Determine whether r has reached the maximum number of iterations iter1 (in this invention, the value is 3);
[0243] If yes, proceed to step seven P10d; otherwise, set r = r + 1 and return to step seven P10b.
[0244] Step 7 P10d, Output x (iter1) As the local maximum point of f(x);
[0245] Step 7, page 11, regarding x(iter1) Perform the processing (i.e., steps three H1 to three H3), and use the processing result as the updated θ. opt2 ;
[0246] Step 7 Q, f max2 Is it less than or equal to the power threshold t? p ;
[0247] If so, then delete Θ. f The index-th element forms the new estimated angle vector Θ. f ;
[0248] Otherwise, vector Θ f Replace the index-th element with θ opt2 And let index = index + 1;
[0249] Step 7 R: Determine if d equals N f ;
[0250] If not, let d = d + 1, and return to step seven C;
[0251] If so, then output the final angle vector Θ. f And continue with step eight;
[0252] Step 8: Determine if it is possible to exit the loop controlled by loop variable c; Specifically:
[0253] Step 8A: Determine the vector Θ fc With vector Θ f Do they have the same number of elements?
[0254] If they are different, do not exit the loop and proceed to step nine; otherwise, continue to step eightB.
[0255] Step 8B: Let vector Θ fc With vector Θ f Let d be the distance vector between them, and let Θ be the k-th element of the distance vector d. fc The k-th element and Θ f The absolute value of the difference between the k-th elements; k = 1, 2, ..., card(Θ) f );
[0256] Step 8C: Determine whether the mean of all elements in the distance vector d is less than the threshold ε3;
[0257] If so, the loop controlled by loop variable c can be exited, and step ten can be executed; otherwise, the loop controlled by loop variable c can not be exited, and step nine can be executed.
[0258] Step 9: Determine if the loop variable c is equal to the set iteration number iter2;
[0259] If yes, continue to step ten; otherwise, set c = c + 1 and return to step seven.
[0260] The iteration count iter2 can be set to balance the system's requirements for real-time performance and reliability. Increasing the value of iter2 can improve system reliability, but it will sacrifice real-time performance, and the reliability gain follows the law of diminishing marginal returns.
[0261] Step 10: Initialize the relative power vector as P r Vector P r The number of elements and the current angle vector Θ f The number of elements is the same, and vector P r All elements are 0;
[0262] Step 11: Calculate the relative power vector P r And according to the relative power vector P r Adjust the current angle vector Θ f The element order in the vector is used to obtain the angle vector Θ after the element order is adjusted. f Specifically:
[0263] Step 11A: Initialize the loop variable p = 1;
[0264] Step 11B: Let the relative power angle to be calculated be θ. w Θ is the angle vector f The p-th element, let the remaining angle vector Θ d For Θ f Remove θ w The vector after;
[0265] Step 11C: Construct the residual angle steering vector matrix A under the joint error. e The construction method is the same as in step seven D;
[0266] Step 11D: Construct the projection matrix P under joint error. Ae The construction method is the same as in step seven E;
[0267] Step 11E: Construct the function f2(θ) of the angle when constructing the joint error. The specific expression has been given in step 7N, where B e = P Ae XX H P Ae ;
[0268] Step 11F, θ w Substitute the values into the function f2(θ), and then use the resulting function value as the vector P.r The p-th element;
[0269] Step 11 G: Determine if p is equal to the angle vector Θ f The number of elements in;
[0270] If yes, continue to step 11H; otherwise, let p = p + 1 and return to step 11B.
[0271] Step 11H: Convert the relative power vector P r The elements in the array are sorted from largest to smallest, and the changes in the positions of the elements before and after the sorting are recorded. The current angle vector Θ is then used. f The elements in the array are rearranged according to the recorded positional changes to obtain the angle vector after the element order is adjusted.
[0272] Step 12: Based on the relative power vector P r and the angle vector Θ after element order adjustment f Reconstruct the original transmitted signal S;
[0273] Step 12A: Construct a copy of the estimated angle vector Θ fc equal to the vector Θ after element order adjustment f Construct a relative power vector replica P rc Equal to the relative power vector P r ;
[0274] Step 12B, if Θ fc If the set is empty, then stop reconstructing the original transmitted signal and let the amplitude-phase matrix copy G... e p is an N x N identity matrix and a copy of the position error vector. e Given a vector of 1 row and N columns containing only zeros, proceed to step fourteen.
[0275] Otherwise, continue reconstructing the original transmitted signal, retaining a copy of the currently estimated angle vector Θ. fc The first N k N elements, that is, deleting the remaining elements, where N k N is a positive integer. k Satisfy: P rc The first N k The sum of the elements and P rc The ratio of the sum of all elements of P is greater than or equal to ε4 (0 < ε4 ≤ 1), and P rc The first N k – The sum of one element and P rc The ratio of the sum of all its elements is less than ε4;
[0276] Step 12C: Construct the reconstructed angle steering vector matrix A under joint error. rThe construction method is the same as in step seven D, except that the vector Θ is used. d Replace with vector Θ fc ;
[0277] Step 12D: Reconstruct the original transmitted signal S. The reconstruction method is as follows:
[0278]
[0279] Step 13: Update the amplitude and phase matrix copy G e and position error vector copy p e The specific update method is as follows:
[0280] Step 13A: Initialize the amplitude-phase matrix replica G e p is an N x N identity matrix and a copy of the position error vector. e Create a 1-row, N-column vector of all zeros; and initialize the loop variable q = 1;
[0281] Step 13B: Let x be the q-th row of X;
[0282] Step 13C: Calculate vector w r :
[0283]
[0284] Step 13D: Calculate vector u r u r The kth element is w r The kth element and The product of θ, where θ r, k Represents vector Θ fc The k-th element, k = 1, 2, ..., card(Θ) fc );
[0285] Step 13E: Take vector u respectively r The phase of each element in the vector is used to obtain the vector a. r ;
[0286] Step 13F: Transfer matrix G e The q-th diagonal element is updated to vector u r The first N in k The mean of the moduli of each element;
[0287] Step 13 G, if N k If N is greater than 1, then proceed to step thirteen H. k If the value is 1, then proceed to step thirteen (I).
[0288] Step 13H, for vector p e Matrix G e Update:
[0289] Let Θ fc The sine value of the midpoint angle is the independent variable, a. r The first N k Performing linear regression with each element as the dependent variable yields the slope k. r and intercept b r and vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ;
[0290] Then proceed to step thirteen (J);
[0291] Step Thirteen I: Preserve vector p e Let b remain unchanged. r Let vector a r The first element will be the matrix G. e The q-th diagonal element is multiplied by the original value. Then proceed to step thirteen (J).
[0292] Step 13J: Determine if q is equal to N;
[0293] If yes, proceed to step thirteen K; otherwise, set q = q + 1 and return to step thirteen B.
[0294] Step 13 K, for vector p e Take the absolute value of each element in the vector |p. e |, to obtain the vector |p e The largest element p in | max Calculate the vector |diag(G) e The largest element G in () – 1| max Where, diag(·) represents the extraction matrix G e The main diagonal elements are combined to form a column vector, diag(G e ) – 1 indicates that each element in the resulting column vector is subtracted from 1, and the resulting vector is formed using the results of these subtractions. |·| indicates that for the vector diag(G) e Take the absolute value of each element in -1;
[0295] Step 13 L, if p max If the value is greater than the preset threshold ε5 (0 < ε5 ≤ 1), then let p e If it is a vector of all zeros, otherwise keep p. e Unchanged; if G max If it is greater than the preset threshold ε6 (0 < ε6 ≤ 1), then let G...e It is an N x N identity matrix, otherwise G e Remain unchanged;
[0296] Then proceed to step fourteen;
[0297] Step 14: Let vector Δp = |p – p e |, vector Δg = |diag(G e ) – diag(G)|;
[0298] Step 15: Determine if it is possible to exit the loop controlled by the loop variable u, specifically:
[0299] Step 15A: Update p = p e G = G e ;
[0300] Step 15B: If the mean of all elements in vector Δp is greater than or equal to ε3 or the mean of all elements in vector Δg is greater than or equal to ε3, then the loop cannot be exited and step 16 should be executed; otherwise, the loop can be exited and step 17 should be executed.
[0301] Step 16: Determine if the loop variable u is equal to the preset maximum number of iterations iter3;
[0302] If yes, continue to step seventeen; otherwise, set u = u + 1 and c = 1, and return to step seven.
[0303] iter3 can be used to balance the system's requirements for real-time performance and reliability. Adding iter3 can improve system reliability, but it will sacrifice real-time performance, and the reliability gain follows the law of diminishing marginal returns.
[0304] Step 17: Within the desired signal range Θ s The search for the desired signal angle involves the following steps:
[0305] Step 17A: If the angle vector Θ has been estimated f If the number of elements in the vector is greater than or equal to 1, then from the estimated angle vector Θ f Delete the region Θ located in the desired signal range s The elements in the vector form the residual angle vector Θ. d Then proceed to step seventeen-b;
[0306] If the angle vector Θ has been estimated f If the number of elements in the matrix is 0, then let the projection matrix P under the joint error be... Ae Given an N x N identity matrix, proceed to step seventeen (F).
[0307] Step 17B: Based on the remaining angle vector Θ dUpdate the projection matrix P under joint error Ae And obtain the candidate angle set Θ ce Specifically, this is achieved by repeating steps seven D to seven M;
[0308] Step 17C: In the candidate angle set Θ ce Only the region containing the desired signal Θ is retained. s The angle in the set is obtained by removing the remaining elements, resulting in the updated set of candidate angles Θ. ce ;
[0309] Step 17.D. If the updated candidate angle set Θ ce If the number of elements located within the desired signal interval is greater than or equal to 1, then proceed to step 17E; otherwise, proceed to step 17F.
[0310] Step 17E: Based on the updated candidate angle set Θ ce Calculate the estimated angle of arrival θ of the desired signal s =θ opt2 Specifically, this is achieved through steps seven N to seven P, followed by step eighteen.
[0311] Step 17F: Let the estimated angle of arrival of the desired signal be θ. s For the desired signal interval Θ s The midpoint, then proceed to step eighteen;
[0312] Step 18: Calculate the beamforming weights. The specific calculation method is as follows:
[0313]
[0314] Among them, P b (θ s Let be an N x N diagonal matrix, and let the nth diagonal element be 0. p n Let n be the nth element of vector p, where n = 1, 2, ..., N;
[0315] Step 19: Output the beamformed signal s to complete interference suppression and desired signal reception.
[0316]
[0317] This invention proposes a beamforming interference suppression method based on partial correlation orthogonal matched pursuit. The partial correlation coefficient between candidate atoms and the current residual is used as a new atom selection criterion for the orthogonal matched pursuit algorithm. Polynomial root finding is used instead of peak search to reduce computational complexity. Finally, alternating optimization and linear regression are combined to deal with joint errors and improve robustness, thereby ensuring reliable communication in non-ideal working environments.
[0318] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that the specific process of step three (B) is as follows:
[0319] Step 3B1: Initialize the coefficient vector v p It is a vector consisting of 1 row, 2N – 1 columns, all zeros;
[0320] Step 3B2: Initialize the loop variable ii = 1;
[0321] Step 3B3: Initialize the loop variable jj = 1;
[0322] Step 3B4: Convert the coefficient vector v p The i-th element plus the projection matrix P A The element in the jjth row and N – ii + jjth column will be the coefficient vector v. p The ii-th element is updated to the summation result;
[0323] Step 3B5: Determine if jj is equal to ii. If not, set jj = jj + 1 and return to step 3B4. If yes, execute step 3B6.
[0324] Step 3B6: Determine if ii is equal to N. If not, set ii = ii + 1 and return to step 3B3. If yes, execute step 3B7.
[0325] Step 3B7: Set the loop variable ii = N + 1;
[0326] Step 3B8: Set the loop variable jj = 1;
[0327] Step 3B9: Convert the coefficient vector v p The i-th element plus matrix P A The element in the ii – N + jjth row and jjth column will be the coefficient vector v. p The ii-th element is updated to the summation result;
[0328] Step 3B10: Determine if jj is equal to 2N – ii. If not, set jj = jj + 1 and return to step 3B9. If yes, execute step 3B11.
[0329] Step 3B11: Determine if ii is equal to 2N – 1. If not, set ii = ii + 1 and return to step 3B8. If yes, execute step 3B12.
[0330] Step 3B12: Output the coefficient vector v of the 2N – 2 degree polynomial. p .
[0331] The other steps and parameters are the same as in Specific Implementation Method 1.
[0332] It should be noted that v q The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix BΓ1. r The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix B. s The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix P A Γ1. v pe The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix G H P Ae G;v qe The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix G H B e GΓ;v re The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix G H B e G;v se The generation method and v p The process is the same; only the matrix P in steps B4 and B9 needs to be changed. A Replace with matrix G H P Ae GΓ.
[0333] Specific Implementation Method Three: This implementation method differs from Specific Implementation Method Two in that the specific process of step threeH is as follows:
[0334] For set S R For each element in the sequence, execute steps 3H1 to 3H3 respectively;
[0335] Step 3H1: Take the phase of the element and divide the result of taking the phase by -π;
[0336] The phase value ranges from [-π, π), where π is the mathematical constant pi.
[0337] Step 3H2: Calculate the arcsine function value of the result obtained in Step 3H1;
[0338] Step 3H3: Multiply the result obtained in Step 3H2 by 180, and then divide the result by π to obtain the candidate angle corresponding to the current element.
[0339] The other steps and parameters are the same as in Specific Implementation Method Two.
[0340] Using the method of this embodiment, set S R After processing each element in the algorithm, a candidate angle set is formed using all the candidate angles obtained.
[0341] Specific Implementation Method Four: This implementation method differs from Specific Implementation Method Three in that, when the array amplitude and phase errors are not ignored and the position errors are ignored, step thirteen H is replaced with:
[0342] Fixed slope k r = 0, take a r The first N k The mean of the elements is used as the intercept b. r and vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ;
[0343] Then proceed to step thirteen.
[0344] The other steps and parameters are the same as in Specific Implementation Method 3.
[0345] Specific Implementation Method Five: This implementation method differs from Specific Implementation Method Three in that, when both array position error and amplitude / phase error are ignored, step thirteen H is replaced with:
[0346] Let Θ fc The sine value of the midpoint angle is the independent variable, a. r The first N k Each element is a dependent variable, and the element "0" is added to both the independent and dependent variables. Then, linear regression is performed to obtain the slope k. r and intercept b r , will vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ;
[0347] Then proceed to step thirteen.
[0348] The other steps and parameters are the same as in Specific Implementation Method 3.
[0349] Specific Implementation Method Six: The beamforming interference suppression method based on partial correlation orthogonal matched tracking described in this implementation method, when both array amplitude and phase errors and position errors are negligible, specifically includes the following steps:
[0350] Step 1: Receive analog signals from free space using N receiving antennas. Each receiving antenna receives an analog signal containing an incident angle located in the interval Θ. s = [ - 5°, The expected signal in [+ 5°] and the interference signals of unknown number, unknown incident angle range and unknown power are processed by the analog signals received by each receiving antenna through its own antenna radio frequency processing link to form N digital sequences.
[0351] in, The angle of arrival of the a priori expected signal shall not exceed 5° from the actual expected signal angle of arrival.
[0352] Step 2: Integrate the N digital sequences formed in Step 1. Specifically, summarize the data from the N receiving antennas at L consecutive sampling points into an N-row, L-column received data matrix X (in this invention, L is 50).
[0353] Step 3: Based on the partial correlation orthogonal matching pursuit criterion and combined with the polynomial root-finding method, estimate the angle of arrival of possible incident signals in space, and obtain the estimated angle vector Θ. f Specifically:
[0354] Step 3A, Initialization: The estimated angle vector Θ f An empty vector, an estimated angle steering vector matrix A is an empty matrix, and a projection matrix P A Given an N x N identity matrix, with loop variable i = 1;
[0355] Step 3B: Based on matrix P A Generate the coefficient vector v of a 2N-2 degree polynomial. p For example, when the highest degree term of the polynomial is 6, the polynomial takes the form of: , , … The coefficient vectors v are respectively p The element value in;
[0356] Step 3C: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix BΓ1. q ;
[0357] Where, the intermediate variable matrix B = P A XX H PA The superscript "H" indicates the conjugate transpose. Matrix Γ1 is an N-row N-column diagonal matrix, and the nth diagonal element of the diagonal matrix Γ1 is n - 1, where n = 1, 2, ..., N.
[0358] Step 3D: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix B. r ;
[0359] Step 3E: Based on matrix P A Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. s ;
[0360] Step 3F: Calculate the coefficient vector v p and coefficient vector v q Convolution, then calculate the coefficient vector v r and coefficient vector v s The convolution, using the coefficient vector v p and coefficient vector v q The convolution result minus the coefficient vector v r and coefficient vector v s The convolution result yields the coefficient vector v of the 4N – 4th order pure spatial spectral peak polynomial. c ;
[0361] In this invention, the peak polynomial of the pure space spectrum is defined as:
[0362]
[0363] in, Let a(z) = [1, z, z] represent the independent variable of the peak polynomial of the pure space spectrum. 2 , …, z N-1 ] T ;
[0364] Step 3 G: Convert the coefficient vector v c The nth element in the equation is considered as the coefficient of the peak polynomial of the pure space spectrum to the power of 4N – 3 – n. Then, each root of the peak polynomial of the pure space spectrum is solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1 to obtain the subtraction results.
[0365] Next, calculate the absolute value of each difference result, and remove the roots corresponding to the absolute values greater than the threshold ε1 (in this invention, the value of ε1 is 0.7). Use the remaining roots to form the set S of polynomial roots. R ;
[0366] Where n = 1, 2, ..., 4N – 3;
[0367] Step 3H, for set SR Each element in the algorithm is processed individually to obtain the candidate angle set Θ. c ;
[0368] Step 3I, if Θ f If the set is empty, proceed to step three K; otherwise, proceed to step three J.
[0369] Step 3J: Delete the candidate angle set Θ c The redundant angles in the calculation are as follows:
[0370] Step 3 J1: Initialize the loop variable b = 1;
[0371] Step 3 J2, for vector Θ f The b-th element in Θ c Find the four elements closest to it, and move the four found elements from Θ. c Delete them to form a new candidate angle set Θ c ;
[0372] Step 3 J3, if b equals vector Θ f If the number of elements in the array is determined, proceed to step three K; otherwise, set b = b + 1 and return to step three J2.
[0373] Step 3K: Construct a function f1(θ) with respect to the angle:
[0374]
[0375] Where f1(θ) is the arrival angle quality metric based on the partial correlation coefficient, also known as the pure space spectrum expression; a(θ) represents the steering vector of the receiving array for angle θ;
[0376]
[0377] Where θ represents the angle. The base of the natural logarithm. The superscript T indicates the imaginary unit;
[0378] Step 3 L: For the candidate angle set Θ c Calculate the function value f1(θ) for each angle in the equation, and output the calculated maximum function value f. max1 Minimum function value f min1 and the maximum function value f max1 The corresponding angle θ opt1 ;
[0379] Step 3 M, if f max1 / f min1If the value is greater than or equal to ε2 (ε2 is 2 in this invention), then proceed to step three N; otherwise, jump directly to step three Q.
[0380] Step 3 N, θ opt1 Add vector Θ f Update the estimated angle steering vector matrix A, where the k-th column of A is Θ. f The guiding vector of the kth element, k = 1, 2, ..., card(Θ) f ), where card(Θ) f ) represents the vector Θ f The number of elements in;
[0381] Step 3: Update the projection matrix P A The specific update formula is as follows:
[0382]
[0383] Among them, I N Represents an N x N identity matrix, where the superscript -1 indicates the inverse of the matrix;
[0384] Step 3P: If i equals N, then execute Step 3Q; otherwise, let i = i + 1 and return to Step 3B.
[0385] Step 3 Q: Output the estimated angle vector Θ f ;
[0386] Step 4, if Θ f If the number of elements in matrix P is greater than or equal to 1, then proceed to steps five through twelve; otherwise, let matrix P... Ad Given an N x N identity matrix, proceed directly to step 10F;
[0387] Step 5: Initialize the loop variable c = 1;
[0388] Step 6: Set the power threshold t p = N / C, where Let X be the average power, C be a preset constant (in this invention, C is 50), and set the iteration stopping threshold ε3 (in this invention, threshold ε3 is 10). -6 );
[0389] Step 7: For the estimated angle vector Θ f The angle in the middle is precisely estimated and redundancy is removed; specifically:
[0390] Step 7A: Set the current vector Θ f The number of elements in the array is denoted as N. f = card(Θ fLet the already estimated angle vector copy Θ fc Equal to Θ f ;
[0391] Step 7B: Initialize the loop variable d = 1, and the current refinement angle index index = 1;
[0392] Step 7C: Let the angle to be refined be θ p For Θ f Let the index of the element be the remaining angle vector Θ. d For Θ f Remove θ p The vector after;
[0393] Step 7D: Construct the remaining angle guidance vector matrix A d The construction process is the same as step three, N, except that Θ is removed. f Replace with Θ d ;
[0394] Step 7E: Constructing and maintaining the remaining angle guidance vector matrix A d The corresponding projection matrix P Ad The construction method is the same as in step 3O, except that the estimated angle steering vector matrix A is replaced with the remaining angle steering vector matrix A. d ;
[0395] Step 7F: Based on matrix P Ad Generate the coefficient vector v of a 2N-2 degree polynomial. pd The generation process is the same as in step 3B, except that matrix P needs to be... A Replace with P Ad ;
[0396] Step 7 G, based on matrix B d Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. qd Among them, the intermediate variable matrix B d = P Ad XX H P Ad The generation process is the same as in step 3C, except that matrix B needs to be replaced with B. d ;
[0397] Step 7H, based on matrix B d Generate the coefficient vector v of a 2N-2 degree polynomial. rd The generation process is the same as in step 3D, except that matrix B needs to be replaced with B. d ;
[0398] Step 7.1: Based on matrix P Ad Γ1 generates the coefficient vector v of a 2N-2 degree polynomial.sd The generation process is the same as step three (E), only the matrix P needs to be changed. A Replace with P Ad ;
[0399] Step 7J, using v pd and v qd Convolution minus v rd and v sd The convolution yields the coefficient vector v of the final 4N – 4th degree pure spatial spectral peak polynomial. cd ;
[0400] Step 7K: Transfer vector v cd The nth element is considered as the coefficient of the 4N – 3 – n power in the peak polynomial of the pure space spectrum. The roots of the peak polynomial of the pure space spectrum are solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1, and the absolute values of each subtraction are taken. Roots whose absolute values are greater than the threshold ε1 are then removed. The roots of the remaining polynomials are used to form a set S. Rd , n = 1, 2,..., 4N – 3;
[0401] Step 7 L, for set S Rd Each element in the equation is processed separately (i.e., steps H1 to H3 are performed) to obtain the candidate angle set Θ. cd ;
[0402] Step 7 M: Delete the candidate angle set Θ cd The redundant angle (the specific steps are the same as steps three J1 to three J3, only Θ needs to be removed) f Replace with Θ d Θ c Replace with Θ cd ), thus obtaining the updated candidate angle set Θ cd ;
[0403] Step 7 N: Update the candidate angle set Θ cd Calculate the function value f1(θ) for each angle in the equation, and output the maximum function value f. max2 and the maximum function value f max2 The corresponding angle θ opt2 ;
[0404] Step 7 O, f max2 Is it less than or equal to the power threshold t? p ;
[0405] If so, then delete Θ. f The index-th element forms the new estimated angle vector Θ. f ;
[0406] Otherwise, vector Θ f Replace the index-th element with θ opt2 And let index = index + 1;
[0407] Step 7 P: Determine if d equals N f ;
[0408] If not, let d = d + 1, and return to step seven C;
[0409] If so, then output the final angle vector Θ. f And continue with step eight;
[0410] Step 8: Determine if it is possible to exit the loop controlled by loop variable c; specifically:
[0411] Step 8A: Determine the vector Θ fc With vector Θ f Do they have the same number of elements?
[0412] If they are different, do not exit the loop and proceed to step nine; otherwise, continue to step eightB.
[0413] Step 8B: Let vector Θ fc With vector Θ f Let d be the distance vector between them, and let Θ be the k-th element of the distance vector d. fc The k-th element and Θ f The absolute value of the difference between the k-th elements; k = 1, 2, ..., card(Θ) f );
[0414] Step 8C: Determine whether the mean of all elements in the distance vector d is less than the threshold ε3;
[0415] If so, the loop controlled by loop variable c can be exited, and step ten can be executed; otherwise, the loop controlled by loop variable c can not be exited, and step nine can be executed.
[0416] Step 9: Determine if the loop variable c is equal to the set iteration number iter2;
[0417] If yes, continue to step ten; otherwise, set c = c + 1 and return to step seven.
[0418] The iteration count iter2 can be set to balance the system's requirements for real-time performance and reliability. Increasing the value of iter2 can improve system reliability, but it will sacrifice real-time performance, and the reliability gain follows the law of diminishing marginal returns.
[0419] Step 10: Within the desired signal range Θ sThe search for the desired signal angle involves the following steps:
[0420] Step 10A: From the estimated angle vector Θ f Delete the region Θ located in the desired signal range s The elements in the vector form the residual angle vector Θ. d ;
[0421] Step 10B: Repeat steps 7D to 7L to obtain the candidate angle set Θ. cd And update matrix P Ad ;
[0422] Step 10C: In the candidate angle set Θ cd Only the region containing the desired signal Θ is retained. s From the given angles, remove the remaining elements to obtain the updated candidate angle set Θ. cd ;
[0423] Step 10D: If the updated candidate angle set Θ cd If the number of elements in the array is greater than or equal to 1, proceed to step 10E; otherwise, proceed to step 10F.
[0424] Step 10E: Execute step 7N to obtain the estimated angle of arrival (θ) of the desired signal. s = θ opt2 Then proceed to step eleven;
[0425] Step 10F: Let the estimated angle of arrival of the desired signal be θ. s For the desired signal interval Θ s At the midpoint, proceed to step eleven;
[0426] Step 11: Calculate beamforming weights:
[0427]
[0428] Step 12: Output the beamformed signal s to complete interference suppression and desired signal reception.
[0429]
[0430] The above examples of the present invention are merely illustrative of the computational model and process of the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is impossible to exhaustively list all possible implementations here. Any obvious variations or modifications derived from the technical solutions of the present invention are still within the scope of protection of the present invention.
Claims
1. A beamforming interference suppression method based on partial correlation orthogonal matched tracking, characterized in that, When the array amplitude and phase errors and position errors are not ignored, the method specifically includes the following steps: Step 1: Use N receiving antennas to receive analog signals from free space. The analog signals received by each receiving antenna are processed by their own antenna radio frequency processing links to form N digital sequences. Step 2: Integrate the N digital sequences formed in Step 1. Specifically, summarize the data from the N receiving antennas at L consecutive sampling points into an N-row, L-column received data matrix X. Step 3: Based on the partial correlation orthogonal matching pursuit criterion and combined with the polynomial root-finding method, estimate the angle of arrival of the incident signal in space to obtain the estimated angle vector Θ. f ; Step 4, if Θ f If the number of elements in the matrix is greater than or equal to 1, proceed to step five; otherwise, let the projection matrix P under the joint error be... Ae Given an N x N identity matrix, an N x N amplitude-phase matrix G, and a 1 x N zero vector p, proceed to step seventeen F. Step 5: Initialize loop variable u = 1, initialize loop variable c = 1; Step 6: Initialize the amplitude and phase matrix G as an N-row, N-column identity matrix and the position error vector p as a 1-row, N-column vector of all zeros; Set power threshold t p = N / C, where Let X be the average power, C be a preset constant, and set the iteration stopping threshold ε3; Step 7: For the estimated angle vector Θ f The angles in the middle are precisely estimated and redundancy is removed; Step 8: Determine whether to exit the loop controlled by loop variable c; specifically: Step 8A: Determine the vector Θ fc With vector Θ f Do they have the same number of elements? If they are different, do not exit the loop and proceed to step nine; otherwise, continue to step eightB. Step 8B: Let vector Θ fc With vector Θ f Let d be the distance vector between them, and let Θ be the k-th element of the distance vector d. fc The k-th element and Θ f The absolute value of the difference between the k-th elements; k = 1, 2, ..., card(Θ) f ); Step 8C: Determine whether the mean of all elements in the distance vector d is less than the threshold ε3; If yes, then exit the loop controlled by loop variable c and continue to step ten; otherwise, do not exit the loop controlled by loop variable c and continue to step nine. Step 9: Determine if the loop variable c is equal to the set iteration number iter2; If yes, continue to step ten; otherwise, set c = c + 1 and return to step seven. Step 10: Initialize the relative power vector as P r Vector P r The number of elements and the current angle vector Θ f The number of elements is the same, and vector P r All elements are 0; Step 11: Calculate the relative power vector P r And according to the relative power vector P r Adjust the current angle vector Θ f The element order in the vector is used to obtain the angle vector Θ after the element order is adjusted. f ; Step 12: Based on the relative power vector P r and the angle vector Θ after element order adjustment f Reconstruct the original transmitted signal S; Step 12A: Construct a copy of the estimated angle vector Θ fc equal to the vector Θ after element order adjustment f Construct a relative power vector replica P rc Equal to the relative power vector P r ; Step 12B, if Θ fc If the set is empty, then stop reconstructing the original transmitted signal and let the amplitude-phase matrix copy G... e p is an N x N identity matrix and a copy of the position error vector. e Given a vector of 1 row and N columns containing only zeros, proceed to step fourteen. Otherwise, continue reconstructing the original transmitted signal, retaining a copy of the currently estimated angle vector Θ. fc The first N k There are N elements, where N k N is a positive integer. k Satisfy: P rc The first N k The sum of the elements and P rc The ratio of the sum of all elements of P is greater than or equal to ε4, and P rc The first N k – The sum of one element and P rc The ratio of the sum of all its elements is less than ε4; Step 12C: Construct the reconstructed angle steering vector matrix A under joint error. r ; Step 12D: Reconstruct the original transmitted signal S. The reconstruction method is as follows: Step 13: Update the amplitude and phase matrix copy G e and position error vector copy p e The specific update method is as follows: Step 13A: Initialize the amplitude-phase matrix replica G e p is an N x N identity matrix and a copy of the position error vector. e Create a 1-row, N-column vector of all zeros; and initialize the loop variable q = 1; Step 13B: Let x be the q-th row of X; Step 13C: Calculate vector w r : Step 13D: Calculate vector u r u r The kth element is w r The kth element and The product of θ, where θ r, k Represents vector Θ fc The k-th element, k = 1, 2, ..., card(Θ) fc ); Step 13E: Take vector u respectively r The phase of each element in the vector is used to obtain the vector a. r ; Step 13F: Transfer matrix G e The q-th diagonal element is updated to vector u r The first N in k The mean of the moduli of each element; Step 13 G, if N k If N is greater than 1, then proceed to step thirteen H. k If the value is 1, then proceed to step thirteen (I). Step 13H, for vector p e Matrix G e Update: Let Θ fc The sine value of the midpoint angle is the independent variable, a. r The first N k Performing linear regression with each element as the dependent variable yields the slope k. r and intercept b r and vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ; Then proceed to step thirteen; Step Thirteen I: Preserve vector p e Let b remain unchanged. r Let vector a r The first element will be the matrix G. e The q-th diagonal element is multiplied by the original value. Then proceed to step thirteen (J). Step 13J: Determine if q is equal to N; If yes, proceed to step thirteen K; otherwise, set q = q + 1 and return to step thirteen B. Step 13 K, Obtain the vector |p e The largest element p in | max Calculate the vector |diag(G) e The largest element G in () – 1| max Where, diag(·) represents the extraction matrix G e The elements of the main diagonal are used to form a column vector, and |·| represents the vector diag(G) e Take the absolute value of each element in 1. Step 13 L, if p max If it is greater than the preset threshold ε5, then let p e If it is a vector of all zeros, otherwise keep p. e Unchanged; if G max If it is greater than the preset threshold ε6, then let G e It is an N x N identity matrix, otherwise G e Remain unchanged; Then proceed to step fourteen; Step 14: Let vector Δp = |p – p e |, vector Δg = |diag(G e ) – diag(G)|; Step 15: Determine whether to exit the loop controlled by loop variable u, specifically: Step 15A: Update p = p e G = G e ; Step 15B: If the mean of all elements in vector Δp is greater than or equal to ε3 or the mean of all elements in vector Δg is greater than or equal to ε3, then do not exit the loop and proceed to step 16; otherwise, exit the loop and continue to step 17. Step 16: Determine if the loop variable u is equal to the preset maximum number of iterations iter3; If yes, continue to step seventeen; otherwise, set u = u + 1 and c = 1, and return to step seven. Step 17: Within the desired signal range Θ s The search for the desired signal angle involves the following steps: Step 17A: If the angle vector Θ has been estimated f If the number of elements in the vector is greater than or equal to 1, then from the estimated angle vector Θ f Delete the region Θ located in the desired signal range s The elements in the vector form the residual angle vector Θ. d Then proceed to step seventeen (B). If the angle vector Θ has been estimated f If the number of elements in the matrix is 0, then let the projection matrix P under the joint error be... Ae Given an N x N identity matrix, proceed to step seventeen (F). Step 17B: Based on the remaining angle vector Θ d Update the projection matrix P under joint error Ae And obtain the candidate angle set Θ ce ; Step 17C: In the candidate angle set Θ ce Only the region containing the desired signal Θ is retained. s The angle in the set is obtained by removing the remaining elements, resulting in the updated set of candidate angles Θ. ce ; Step 17.D. If the updated candidate angle set Θ ce If the number of elements located within the desired signal interval is greater than or equal to 1, then proceed to step 17E; otherwise, proceed to step 17F. Step 17E: Based on the updated candidate angle set Θ ce Calculate the estimated angle of arrival θ of the desired signal s = θ opt2 Then proceed to step eighteen; Step 17F: Let the estimated angle of arrival of the desired signal be θ. s For the desired signal interval Θ s The midpoint, then proceed to step eighteen; Step 18: Calculate the beamforming weights. The specific calculation method is as follows: Among them, P b (θ s Let be an N x N diagonal matrix, and let the nth diagonal element be 0. p n Let n be the nth element of vector p, where n = 1, 2, ..., N; Step 19: Output the beamformed signal s to complete interference suppression and desired signal reception. 。 2. The beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 1, characterized in that, The specific process of step three is as follows: Step 3A, Initialization: The estimated angle vector Θ f An empty vector, an estimated angle steering vector matrix A is an empty matrix, and a projection matrix P A Given an N x N identity matrix, with loop variable i = 1; Step 3B: Based on matrix P A Generate the coefficient vector v of a 2N-2 degree polynomial. p ; Step 3C: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix BΓ1. q ; Where, the intermediate variable matrix B = P A XX H P A The superscript "H" indicates the conjugate transpose. Matrix Γ1 is an N-row N-column diagonal matrix, and the nth diagonal element of the diagonal matrix Γ1 is n - 1, where n = 1, 2, ..., N. Step 3D: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix B. r ; Step 3E: Based on matrix P A Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. s ; Step 3F: Calculate the coefficient vector v p and coefficient vector v q Convolution, then calculate the coefficient vector v r and coefficient vector v s The convolution, using the coefficient vector v p and coefficient vector v q The convolution result minus the coefficient vector v r and coefficient vector v s The convolution result yields the coefficient vector v of the 4N – 4th order pure spatial spectral peak polynomial. c ; Step 3 G: Convert the coefficient vector v c The nth element in the equation is considered as the coefficient of the peak polynomial of the pure space spectrum to the power of 4N – 3 – n. Then, each root of the peak polynomial of the pure space spectrum is solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1 to obtain the subtraction results. Next, calculate the absolute value of each difference result, remove the roots corresponding to the absolute values greater than the threshold ε1, and use the remaining roots to form the set S of polynomial roots. R ; Where n = 1, 2, ..., 4N – 3; Step 3H, for set S R Each element in the algorithm is processed individually to obtain the candidate angle set Θ. c ; Step 3I, if Θ f If the set is empty, proceed to step three K; otherwise, proceed to step three J. Step 3J: Delete the candidate angle set Θ c The redundant angles in the calculation are as follows: Step 3 J1: Initialize the loop variable b = 1; Step 3 J2, for vector Θ f The b-th element in Θ c Find the four elements closest to it, and move the four found elements from Θ. c Delete them to form a new candidate angle set Θ c ; Step 3 J3, if b equals vector Θ f If the number of elements in the array is determined, proceed to step three K; otherwise, set b = b + 1 and return to step three J2. Step 3K: Construct a function f1(θ) with respect to the angle: Where a(θ) represents the guide vector of the receiving array with respect to angle θ; Where θ represents the angle. The base of the natural logarithm. The superscript T indicates the imaginary unit; Step 3 L: For the candidate angle set Θ c Calculate the function value f1(θ) for each angle in the equation, and output the calculated maximum function value f. max1 Minimum function value f min1 and the maximum function value f max1 The corresponding angle θ opt1 ; Step 3 M, if f max1 / f min1 If the value is greater than or equal to ε2, then proceed to step three N; otherwise, jump directly to step three Q. Step 3 N, θ opt1 Add vector Θ f Update the estimated angle steering vector matrix A, where the k-th column of A is Θ. f The guiding vector of the kth element, k = 1, 2, ..., card(Θ) f ), where card(Θ f ) represents the vector Θ f The number of elements in; Step 3: Update the projection matrix P A The specific update formula is as follows: Among them, I N Represents an N x N identity matrix, where the superscript -1 indicates the inverse of the matrix; Step 3P: If i equals N, then execute Step 3Q; otherwise, let i = i + 1 and return to Step 3B. Step 3 Q: Output the estimated angle vector Θ f .
3. The beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 2, characterized in that, The specific process of step 3B is as follows: Step 3B1: Initialize the coefficient vector v p It is a vector consisting of 1 row, 2N – 1 columns, all zeros; Step 3B2: Initialize the loop variable ii = 1; Step 3B3: Initialize the loop variable jj = 1; Step 3B4: Convert the coefficient vector v p The i-th element plus the projection matrix P A The element in the jjth row and N – ii + jjth column will be the coefficient vector v. p The ii-th element is updated to the summation result; Step 3B5: Determine if jj is equal to ii. If not, set jj = jj + 1 and return to step 3B4. If yes, execute step 3B6. Step 3B6: Determine if ii is equal to N. If not, set ii = ii + 1 and return to step 3B3. If yes, execute step 3B7. Step 3B7: Set the loop variable ii = N + 1; Step 3B8: Set the loop variable jj = 1; Step 3B9: Convert the coefficient vector v p The i-th element plus matrix P A The element in the ii – N + jjth row and jjth column will be the coefficient vector v. p The ii-th element is updated to the summation result; Step 3B10: Determine if jj is equal to 2N – ii. If not, set jj = jj + 1 and return to step 3B9. If yes, execute step 3B11. Step 3B11: Determine if ii is equal to 2N – 1. If not, set ii = ii + 1 and return to step 3B8. If yes, execute step 3B12. Step 3B12: Output the coefficient vector v of the 2N – 2 degree polynomial. p .
4. The beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 3, characterized in that, The specific process of step three H is as follows: For set S R For each element in the sequence, execute steps 3H1 to 3H3 respectively; Step 3H1: Take the phase of the element and divide the result of taking the phase by -π; Step 3H2: Calculate the arcsine function value of the result obtained in Step 3H1; Step 3H3: Multiply the result obtained in Step 3H2 by 180, and then divide the result by π to obtain the candidate angle corresponding to the current element.
5. The beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 4, characterized in that, The specific process of step seven is as follows: Step 7A: Set the current vector Θ f The number of elements in the array is denoted as N. f = card(Θ f Let the already estimated angle vector copy Θ fc Equal to Θ f ; Step 7B: Initialize the loop variable d = 1, and the current refinement angle index index = 1; Step 7C: Let the angle to be refined be θ p For Θ f Let the index of the element be the remaining angle vector Θ. d For Θ f Remove θ p The vector after; Step 7D: Construct the residual angle steering vector matrix A under the joint error. e Matrix A e The number of rows is N and the number of columns is card(Θ). d ), let Θ d The kth element is θ d,k Then A e The k-th column is represented as: Among them, P b (θ d,k Let P be an N x N diagonal matrix. b (θ d,k The nth diagonal element is p n Let n be the nth element of vector p, where n = 1, 2, ..., N, k = 1, 2, ..., card(Θ). d ); Step 7E: Construct the projection matrix P under joint error. Ae : Step 7F: Based on matrix G H P Ae G generates the coefficient vector v of a 2N-degree polynomial. pe ; Step 7 G, based on matrix G H B e GΓ generates the coefficient vector v of a 2N-degree polynomial. qe ; Wherein, matrix B e = P Ae XX H P Ae , matrix Γ = Γ1 + Diag(p), Diag(·) means to transform vector p into a diagonal matrix, and the nth diagonal element of the diagonal matrix Diag(p) is the nth element of vector p, n = 1, 2, ..., N; Step 7H, based on matrix G H B e G generates the coefficient vector v of a 2N-degree polynomial. re ; Step 7.1: Based on matrix G H P Ae GΓ generates the coefficient vector v of a 2N-degree polynomial. se ; Step 7J, using v pe and v qe Convolution minus v re and v se The convolution yields the coefficient vector v of the final 4N – 4th degree pure spatial spectral peak polynomial. ce ; Step 7K: Transfer vector v ce The nth element is considered as the coefficient of the 4N – 3 – n power in the peak polynomial of the pure space spectrum. The roots of the peak polynomial of the pure space spectrum are solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1, and the absolute values of each subtraction are taken. Roots whose absolute values are greater than the threshold ε1 are then removed. The roots of the remaining polynomials are used to form a set S. Re , n = 1, 2,..., 4N – 3; Step 7 L, for set S Re Each element in the algorithm is processed individually to obtain the candidate angle set Θ. ce ; Step 7 M: Delete the candidate angle set Θ ce The redundant angles in the equation are used to obtain the updated candidate angle set Θ. ce ; Step 7N: Construct the function f2(θ) of the angle when constructing the joint error: Step 7.0: Update the candidate angle set Θ ce Calculate the function value of f2(θ) for each angle in the equation, and output the maximum function value f. max2 and the maximum function value f max2 The corresponding angle θ opt2 ; Step 7 P, for angle θ opt2 The following corrections have been made: Step 7 P1, based on matrix G H P Ae G generates the coefficient vector v of an N(N – 1) term power function. pl Specifically: Step 7 P1a: Initialize the coefficient vector v pl It is a vector consisting of 1 row and N(N – 1) + 1 columns, all containing zeros; Step 7 P1b: Initialize loop variable uu = 1, loop variable vv = 1, counter count = 0; Step 7 P1c: If uu = vv, then v pl The N(N-1)+1th element is updated to v. pl The N(N-1)+1th element and matrix G H P Ae The sum of the values in row uu and column v of G is used; otherwise, count = count + 1, and v is set to... pl The count-th element is assigned to matrix G. H P Ae The value of the element in row uu and column vv of G; Step 7 P1d: Determine if vv is equal to N. If yes, proceed to step 7 P1e; otherwise, set vv = vv + 1 and return to step 7 P1c. Step 7 P1e: Determine if uu is equal to N; If yes, proceed to step seven P1f; otherwise, set uu = uu + 1 and vv = 1, and return to proceed to step seven P1c. Step 7 P1f: Output the coefficient vector v of the N(N – 1) term power function. pl ; Step 7 P2, based on matrix G H B e GΓ generates the coefficient vector v of an N(N – 1) term power function. ql ; Step 7 P3, based on matrix G H B e G generates the coefficient vector v of an N(N – 1) term power function. rl ; Step 7 P4, based on matrix G H P Ae GΓ generates the coefficient vector v of an N(N – 1) term power function. sl ; Step 7, P5: Construct the power vector t1 of the N(N – 1) term power function. The dimension of the power vector t1 is 1 row and N(N – 1) + 1 column. The specific construction method is as follows: Step 7 P5a: Construct a position vector p with array position error. v Vector p v The nth element is the nth element of vector p plus n – 1, where n = 1, 2, ..., N; Step 7 (P5b): Initialize the loop variable = 1, loop variable f = 1, counter count = 0; Step 7 P5c, if If the value is f, then the N(N – 1) + 1th element of the exponentiation vector t1 is assigned the value 0; otherwise, let count = count + 1, and assign the countth element of the exponentiation vector t1 the value p. v The f-th element and p v The The difference between the elements; Step 7 P5d: Determine if f is equal to N; If yes, proceed to step seven P5e; otherwise, let f = f + 1 and return to proceed to step seven P5c. Step 7 P5e, Judgment Is it equal to N? If so, proceed to step seven (P5f); otherwise, let... = + 1, f = 1, and return to execute step seven P5c; Step 7 P5f: Obtain the maximum value in the power vector t1, and add the obtained maximum value to each element in the power vector t1 to obtain the output power vector t1. Step 7, P6: According to v pl and v ql Construct N(N – 1)(N 2 – N + 3) / 2 term power function coefficient vector v b Specifically: Step 7 (P6a): Initialize the coefficient vector v b For a single row N(N – 1)(N 2 – A vector of all zeros consisting of N + 3) / 2 + 1 columns; Step 7 (P6b): Initialize loop variable g = 1; Initialize loop variable h = g; Initialize counter count = 0; Step 7 (P6c): Let count = count + 1; Step 7 P6d: If g = h, then let v b The count-th element is v pl The g-th element and v ql The product of h elements; otherwise, let v b The count-th element is v pl The g-th element and v ql The product of the h-th elements plus v pl The h-th element and v ql The product of the g-th elements; Step 7 P6e: Determine if h is equal to N(N – 1) + 1; If yes, proceed to step seven P6f; otherwise, set h = h + 1 and return to step seven P6c. Step 7 P6f: Determine if g is equal to N(N – 1) + 1. If yes, execute step 7 P6g; otherwise, let g = g + 1 and h = g, and return to execute step 7 P6c. Step 7 P6g, output according to v pl and v ql Construct N(N – 1)(N) 2 – N + 3) / 2 term power function coefficient vector v b ; Step 7 (P7): According to v rl and v sl Construct N(N – 1)(N 2 – N + 3) / 2 term power function coefficient vector v a ; Step 7, P8: Construct N(N-1)(N) based on the power vector t1 of the N(N-1) term power function. 2 – N + 3) / 2 power function power vector t2, power vector t2 has a dimension of 1 row N(N – 1)(N 2 – N + 3) / 2 + 1 columns, the specific construction method is as follows: Step 7 (P8a): Initialize loop variable l = 1; Initialize loop variable m = l; Initialize counter count = 0; Step 7 (P8b): Let count = count + 1; Step 7 P8c: Assign the count-th element of t2 to the sum of the l-th element and the m-th element of t1; Step 7 P8d: Determine if m is equal to N(N – 1) + 1; If yes, proceed to step seven P8e; otherwise, set m = m + 1 and return to step seven P8b. Step 7 P8e: Determine if l is equal to N(N – 1) + 1; If yes, proceed to step seven P8f; otherwise, let l = l + 1, m = l, and return to proceed to step seven P8b. Step 7 P8f, Output N(N – 1)(N 2 – N + 3) / 2 power vectors t2; Step 7, P9: According to v b v a The peak power function f(x) of the pure spatial spectrum under joint error and the derivative function g(x) of f(x) are constructed as follows: Where x represents the independent variable of the function, Indicates v a The nth element, Indicates v b Let t2(n) represent the nth element of t2, where n = 1, 2, ..., N(N – 1)(N 2 – N + 3) / 2 + 1; N’ = N(N – 1)(N 2 – N +3) / 2 + 1; Step 7, page 10: Use Newton's method to iteratively find the maximum points of the peak power function f(x) of the pure space spectrum. The specific iterative steps are as follows: Step 7 (P10a): Initialize the iteration variable r = 1, and the initial iteration value x. (0) = ; Step 7 (P10b): Calculate x (r) = x (r – 1) – f[x (r – 1) ] / g[x (r – 1) ]; Step 7 P10c: Determine if r has reached the maximum number of iterations iter1; If yes, proceed to step seven P10d; otherwise, set r = r + 1 and return to step seven P10b. Step 7 P10d, Output x (iter1) As the local maximum point of f(x); Step 7, page 11, regarding x (iter1) Perform processing, and use the processing result as the updated θ. opt2 ; Step 7 Q, f max2 Is it less than or equal to the power threshold t? p ; If so, then delete Θ. f The index-th element forms the new estimated angle vector Θ. f ; Otherwise, vector Θ f Replace the index-th element with θ opt2 And let index = index + 1; Step 7 R: Determine if d equals N f ; If not, let d = d + 1, and return to step seven C; If so, then output the final angle vector Θ. f Then proceed to step eight.
6. The beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 5, characterized in that, The specific process of step eleven is as follows: Step 11A: Initialize the loop variable p = 1; Step 11B: Let the relative power angle to be calculated be θ. w Θ is the angle vector f The p-th element, let the remaining angle vector Θ d For Θ f Remove θ w The vector after; Step 11C: Construct the residual angle steering vector matrix A under the joint error. e ; Step 11D: Construct the projection matrix P under joint error. Ae ; Step 11E: Construct the function f2(θ) of the angle when calculating the joint error; Step 11F, θ w Substitute the values into the function f2(θ), and then use the resulting function value as the vector P. r The p-th element; Step 11 G: Determine if p is equal to the angle vector Θ f The number of elements in; If yes, continue to step 11H; otherwise, let p = p + 1 and return to step 11B. Step 11H: Convert the relative power vector P r The elements in the array are sorted from largest to smallest, and the changes in the positions of the elements before and after the sorting are recorded. The current angle vector Θ is then used. f The elements in the array are rearranged according to their recorded positional changes to obtain the angle vector after the element order is adjusted.
7. The beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 6, characterized in that, When both array amplitude and phase errors and position errors are ignored, step thirteen H is replaced with: Fixed slope k r = 0, take a r The first N k The mean of the elements is used as the intercept b. r and vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ; Then proceed to step thirteen.
8. A beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 6, characterized in that, When both array position error and amplitude / phase error are ignored, step thirteen H is replaced with: Let Θ fc The sine value of the midpoint angle is the independent variable, a. r The first N k Each element is a dependent variable, and the element "0" is added to both the independent and dependent variables. Then, linear regression is performed to obtain the slope k. r and intercept b r , will vector p e The q-th element is updated to -k r / π, will matrix G e The q-th diagonal element is multiplied by the original value. ; Then proceed to step thirteen.
9. A beamforming interference suppression method based on partial correlation orthogonal matched tracking, characterized in that, When the array amplitude and phase errors and position errors are ignored, the method specifically includes the following steps: Step 1: Use N receiving antennas to receive analog signals from free space. The analog signals received by each receiving antenna are processed by their own antenna radio frequency processing links to form N digital sequences. Step 2: Integrate the N digital sequences formed in Step 1. Specifically, summarize the data from the N receiving antennas at L consecutive sampling points into an N-row, L-column received data matrix X. Step 3: Based on the partial correlation orthogonal matching pursuit criterion and combined with the polynomial root-finding method, estimate the angle of arrival of the incident signal in space to obtain the estimated angle vector Θ. f ; Step 4, if Θ f If the number of elements in matrix P is greater than or equal to 1, then proceed to steps five through twelve; otherwise, let matrix P... Ad Given an N x N identity matrix, proceed directly to step 10F; Step 5: Initialize the loop variable c = 1; Step 6: Set the power threshold t p = N / C, where Let X be the average power, C be a preset constant, and set the iteration stopping threshold ε3; Step 7: Perform precise estimation and redundancy removal on the angle of arrival of the incident signal; specifically: Step 7A: Set the current vector Θ f The number of elements in the array is denoted as N. f = card(Θ f Let the already estimated angle vector copy Θ fc Equal to Θ f ; Step 7B: Initialize the loop variable d = 1, and the current refinement angle index index = 1; Step 7C: Let the angle to be refined be θ p For Θ f Let the index of the element be the remaining angle vector Θ. d For Θ f Remove θ p The vector after; Step 7D: Construct the remaining angle guidance vector matrix A d ; Step 7E: Constructing and maintaining the residual angle guidance vector matrix A d The corresponding projection matrix P Ad ; Step 7F: Based on matrix P Ad Generate the coefficient vector v of a 2N-2 degree polynomial. pd ; Step 7 G, based on matrix B d Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. qd ; Step 7H, based on matrix B d Generate the coefficient vector v of a 2N-2 degree polynomial. rd ; Step 7.1: Based on matrix P Ad Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. sd ; Step 7J, using v pd and v qd Convolution minus v rd and v sd The convolution yields the coefficient vector v of the final 4N – 4th degree pure spatial spectral peak polynomial. cd ; Step 7K: Transfer vector v cd The nth element is considered as the coefficient of the 4N – 3 – n power in the peak polynomial of the pure space spectrum. The roots of the peak polynomial of the pure space spectrum are solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1, and the absolute values of each subtraction are taken. Roots whose absolute values are greater than the threshold ε1 are then removed. The roots of the remaining polynomials are used to form a set S. Rd , n = 1, 2,..., 4N – 3; Step 7 L, for set S Rd Each element in the algorithm is processed individually to obtain the candidate angle set Θ. cd ; Step 7 M: Delete the candidate angle set Θ cd The redundant angles in the equation are used to obtain the updated candidate angle set Θ. cd ; Step 7 N: Update the candidate angle set Θ cd Calculate the function value f1(θ) for each angle in the equation, and output the maximum function value f. max2 and the maximum function value f max2 The corresponding angle θ opt2 ; Step 7 O, f max2 Is it less than or equal to the power threshold t? p ; If so, then delete Θ. f The index-th element forms the new estimated angle vector Θ. f ; Otherwise, vector Θ f Replace the index-th element with θ opt2 And let index = index + 1; Step 7 P: Determine if d equals N f ; If not, let d = d + 1, and return to step seven C; If so, then output the final angle vector Θ. f And continue with step eight; Step 8: Determine whether to exit the loop controlled by loop variable c; specifically: Step 8A: Determine the vector Θ fc With vector Θ f Do they have the same number of elements? If they are different, do not exit the loop and proceed to step nine; otherwise, continue to step eightB. Step 8B: Let vector Θ fc With vector Θ f Let d be the distance vector between them, and let Θ be the k-th element of the distance vector d. fc The k-th element and Θ f The absolute value of the difference between the k-th elements; k = 1, 2, ..., card(Θ) f ); Step 8C: Determine whether the mean of all elements in the distance vector d is less than the threshold ε3; If yes, then exit the loop controlled by loop variable c and continue to step ten; otherwise, do not exit the loop controlled by loop variable c and continue to step nine. Step 9: Determine if the loop variable c is equal to the set iteration number iter2; If yes, continue to step ten; otherwise, set c = c + 1 and return to step seven. Step 10: Within the desired signal range Θ s The search for the desired signal angle involves the following steps: Step 10A: From the estimated angle vector Θ f Delete the region Θ located in the desired signal range s The elements in the vector form the residual angle vector Θ. d ; Step 10B: Repeat steps 7D to 7L to obtain the candidate angle set Θ. cd And update matrix P Ad ; Step 10C: In the candidate angle set Θ cd Only the region containing the desired signal Θ is retained. s From the given angles, remove the remaining elements to obtain the updated candidate angle set Θ. cd ; Step 10D: If the updated candidate angle set Θ cd If the number of elements in the array is greater than or equal to 1, proceed to step 10E; otherwise, proceed to step 10F. Step 10E: Execute step 7N to obtain the estimated angle of arrival (θ) of the desired signal. s = θ opt2 Then proceed to step eleven; Step 10F: Let the estimated angle of arrival of the desired signal be θ. s For the desired signal interval Θ s Find the midpoint, then proceed to step eleven; Step 11: Calculate beamforming weights: Step 12: Output the beamformed signal s to complete interference suppression and desired signal reception. 。 10. A beamforming interference suppression method based on partial correlation orthogonal matched tracking according to claim 9, characterized in that, The specific process of step three is as follows: Step 3A, Initialization: The estimated angle vector Θ f An empty vector, an estimated angle steering vector matrix A is an empty matrix, and a projection matrix P A Given an N x N identity matrix, with loop variable i = 1; Step 3B: Based on matrix P A Generate the coefficient vector v of a 2N-2 degree polynomial. p ; Step 3C: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix BΓ1. q ; Where, the intermediate variable matrix B = P A XX H P A The superscript "H" indicates the conjugate transpose. Matrix Γ1 is an N-row N-column diagonal matrix, and the nth diagonal element of the diagonal matrix Γ1 is n - 1, where n = 1, 2, ..., N. Step 3D: Generate the coefficient vector v of the 2N – 2 degree polynomial based on matrix B. r ; Step 3E: Based on matrix P A Γ1 generates the coefficient vector v of a 2N-2 degree polynomial. s ; Step 3F: Calculate the coefficient vector v p and coefficient vector v q Convolution, then calculate the coefficient vector v r and coefficient vector v s The convolution, using the coefficient vector v p and coefficient vector v q The convolution result minus the coefficient vector v r and coefficient vector v s The convolution result yields the coefficient vector v of the 4N – 4th order pure spatial spectral peak polynomial. c ; Step 3 G: Convert the coefficient vector v c The nth element in the equation is considered as the coefficient of the peak polynomial of the pure space spectrum to the power of 4N – 3 – n. Then, each root of the peak polynomial of the pure space spectrum is solved, and the modulus of each root is calculated. The modulus of each root is subtracted from 1 to obtain the subtraction results. Next, calculate the absolute value of each difference result, remove the roots corresponding to the absolute values greater than the threshold ε1, and use the remaining roots to form the set S of polynomial roots. R ; Where n = 1, 2, ..., 4N – 3; Step 3H, for set S R Each element in the algorithm is processed individually to obtain the candidate angle set Θ. c ; Step 3I, if Θ f If the set is empty, proceed to step three K; otherwise, proceed to step three J. Step 3J: Delete the candidate angle set Θ c The redundant angles in the calculation are as follows: Step 3 J1: Initialize the loop variable b = 1; Step 3 J2, for vector Θ f The b-th element in Θ c Find the four elements closest to it, and move the four found elements from Θ. c Delete them to form a new candidate angle set Θ c ; Step 3 J3, if b equals vector Θ f If the number of elements in the array is determined, proceed to step three K; otherwise, set b = b + 1 and return to step three J2. Step 3K: Construct a function f1(θ) with respect to the angle: Where a(θ) represents the guide vector of the receiving array with respect to angle θ; Where θ represents the angle. The base of the natural logarithm. The superscript T indicates the imaginary unit; Step 3 L: For the candidate angle set Θ c Calculate the function value f1(θ) for each angle in the equation, and output the calculated maximum function value f. max1 Minimum function value f min1 and the maximum function value f max1 The corresponding angle θ opt1 ; Step 3 M, if f max1 / f min1 If the value is greater than or equal to ε2, then proceed to step three N; otherwise, jump directly to step three Q. Step 3 N, θ opt1 Add vector Θ f Update the estimated angle steering vector matrix A, where the k-th column of A is Θ. f The guiding vector of the kth element, k = 1, 2, ..., card(Θ) f ), where card(Θ f ) represents the vector Θ f The number of elements in; Step 3: Update the projection matrix P A The specific update formula is as follows: Among them, I N Represents an N x N identity matrix, where the superscript -1 indicates the inverse of the matrix; Step 3P: If i equals N, then execute Step 3Q; otherwise, let i = i + 1 and return to Step 3B. Step 3 Q: Output the estimated angle vector Θ f .