CN110361691B - Implementation method of coherent source DOA estimation FPGA based on non-uniform array - Google Patents
Implementation method of coherent source DOA estimation FPGA based on non-uniform array Download PDFInfo
- Publication number
- CN110361691B CN110361691B CN201910669180.5A CN201910669180A CN110361691B CN 110361691 B CN110361691 B CN 110361691B CN 201910669180 A CN201910669180 A CN 201910669180A CN 110361691 B CN110361691 B CN 110361691B
- Authority
- CN
- China
- Prior art keywords
- module
- matrix
- calculation
- value
- characteristic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 20
- 230000001427 coherent effect Effects 0.000 title claims abstract description 18
- 239000011159 matrix material Substances 0.000 claims abstract description 209
- 238000004364 calculation method Methods 0.000 claims abstract description 114
- 238000001228 spectrum Methods 0.000 claims abstract description 28
- 238000012545 processing Methods 0.000 claims abstract description 26
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 25
- 238000005070 sampling Methods 0.000 claims abstract description 6
- 230000007704 transition Effects 0.000 claims description 63
- 239000013598 vector Substances 0.000 claims description 56
- 238000004140 cleaning Methods 0.000 claims description 42
- 238000009499 grossing Methods 0.000 claims description 22
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 9
- 230000003044 adaptive effect Effects 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 230000003595 spectral effect Effects 0.000 abstract description 7
- 238000013461 design Methods 0.000 abstract description 4
- 230000006870 function Effects 0.000 description 14
- 238000010586 diagram Methods 0.000 description 9
- 238000004088 simulation Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 101100322581 Caenorhabditis elegans add-1 gene Proteins 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D30/00—Reducing energy consumption in communication networks
- Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks
Abstract
The invention belongs to the technical field of array signal processing, and particularly relates to a coherent source DOA estimation FPGA implementation method based on a non-uniform array. The system comprises a self-adaptive covariance matrix calculation module, a multistage pipeline feature decomposition module, a double-threshold information source number estimation module and a spectrum peak search module; the input of the self-adaptive covariance matrix calculation module is sampling data, and the output of the self-adaptive covariance matrix calculation module is connected with the input of the multistage pipelining characteristic decomposition module; the output of the multi-stage pipelining characteristic decomposition module is connected with the input of the double-threshold information source number estimation module, and the output of the double-threshold information source number estimation module is connected with the spectral peak search module; the DOA estimation of the central symmetry non-uniform array coherent source can be realized, different data formats can be adopted, the design is carried out according to the requirement, and the calculation of covariance matrix and parallel Jacobian calculation scheme according to the snapshot number is adopted, so that the calculation time is shortened; the design idea of adopting the pipeline structure optimizes the resources and the speed to the maximum degree on the premise of ensuring the operation precision.
Description
Technical Field
The invention belongs to the technical field of array signal processing, and particularly relates to a coherent source DOA estimation FPGA implementation method based on a non-uniform array.
Background
The estimation of the direction of the wave is one of important branches of array signal processing, and has wide application prospect in the fields of radar, communication, sonar and the like. The uniform array has the defects of small array aperture, low direction-finding precision, poor resolution and incapability of meeting the minimum array element spacing smaller than the half-wavelength non-fuzzy direction-finding condition of an incident information source due to the physical size of an antenna in practical application, and is difficult to be suitable for the modern complex electromagnetic environment. Compared with a uniform array with the same array element number, the non-uniform array can have larger aperture, higher resolution and higher degree of freedom, the reasonable array structure can avoid direction finding ambiguity and improve direction finding performance, and the array aperture can be enlarged, the angle finding precision can be improved, decoherence can be carried out by utilizing space smoothing, and the decoherence capacity of the non-uniform array can be improved by adopting the symmetrical array element arrangement or array element interval combination arrangement array.
The MUSIC algorithm is mainly divided into four parts of covariance calculation, characteristic decomposition, information source number estimation and spectrum peak search, a large amount of calculation is needed for the characteristic decomposition and the spectrum peak search, the time required by DOA estimation is difficult to meet the requirement of real-time direction finding at ms level by adopting a DSP, the parallelism of DOA estimation is high by adopting an FPGA, the DOA estimation can be realized by adopting running water, the time required by DOA estimation is greatly reduced, and the us-level DOA estimation is realized. In the existing research of the MUSIC algorithm, a covariance calculation part adopts a parallel-execution grouping scheme, a control circuit is complex, snapshot variable calculation cannot be realized, a corresponding decorrelation module is not designed, and the consideration of non-centrosymmetric array real number is less. The feature decomposition part comprises a large number of nonlinear operations, the existing serial Jacobi operation is used for realizing the time consumption, and the parallel Jacobi algorithm is used for realizing the feature decomposition, so that the time consumption of the feature decomposition is reduced. The spectral peak search section generally obtains a steering vector by using a memory circuit, calculates a spatial spectral function by directly reading the value, and consumes a large amount of memory resources.
Disclosure of Invention
The invention aims to provide a method for realizing coherent source DOA estimation FPGA based on a non-uniform array, which is suitable for a non-uniform array which can be smoothly decohered, and the feasibility is verified by hardware design realized on an FPGA platform.
The purpose of the invention is realized in the following way:
a method for realizing coherent source DOA estimation FPGA based on a non-uniform array comprises the following steps:
step 1: inputting the received signals into an adaptive covariance matrix calculation module, and calculating to obtain an adaptive covariance matrix of the received signals;
step 2: inputting the obtained self-adaptive covariance matrix into a multi-stage running water characteristic decomposition module to decompose the multi-running water characteristic values, so as to obtain characteristic values and characteristic vectors of the covariance matrix;
step 3: inputting the obtained characteristic value and the characteristic vector into a double-threshold information source number estimation module, judging the information source number and obtaining a noise subspace;
step 4: and outputting the double-threshold information source number estimation module, inputting the output of the double-threshold information source number estimation module into a spectrum peak searching module, and calculating to obtain a stored value obtained after the spectrum function judges and compares all extreme points, namely the direction in which the information source is located.
The self-adaptive covariance matrix calculation module in the step 1 comprises an autocorrelation matrix calculation module, a smoothing processing module and a real-realization processing module;
the autocorrelation matrix calculation module adopts multipath parallel calculation to obtain the calculation result of each snapshot, and adds up the calculation result of each snapshot to obtain the triangle element value on the autocorrelation matrix, and the autocorrelation matrix calculation formula of the received data is as follows:
wherein L is the number of shots,
the calculation formula of (2) is as follows:
wherein 
(p, q=1,., 6); the autocorrelation matrix calculation module comprises a complex multiplier and an accumulator; the input of the complex multiplier is corresponding sampling data, an accumulator accumulates the single snapshot result, and the input snapshot number L is used for controlling to obtain the calculation result of the autocorrelation matrix;
the smoothing processing module performs summation operation on the data of the corresponding position of the autocorrelation matrix according to the space smoothing principle, and performs summation operation on the autocorrelation matrix obtained by calculation
Corresponding toThe elements of the positions are summed, and the calculation formula is:
wherein rpq (p, q=1, 2,., 6) is an autocorrelation matrix
Elements of the corresponding position in (a);
the real number processing module performs a negation operation on the smoothed partial result, converts the covariance matrix of the signal from a complex domain to a real number domain, and obtains ru=λu according to the characteristic decomposition property, wherein λ is a characteristic value of the covariance matrix, R is the covariance matrix of the signal and is the complex matrix, and the representation is r=r r +iR i ，R r Representing the real part of matrix R, R i The imaginary part of the matrix R, u is the eigenvector of the covariance matrix, u is denoted u=u r +iu i ，u r Representing the real part of the vector u, u i The imaginary part representing the vector u, is taken into the above equation:
matrix is formed
According to->
The conversion is to convert the complex matrix into a real matrix, solve the opposite number of the imaginary part of the triangular element on the smoothed matrix and output the calculation result of the triangular element on the covariance matrix to be decomposed.
The multistage pipelining feature decomposition module in step 2 is used for calculating feature values and feature vectors of a covariance matrix, and comprises a multistage cleaning module, wherein the input of the multistage cleaning module is the covariance matrix, and the output of the multistage cleaning module is data of main diagonal lines in a feature value transition matrix and a feature vector matrix; and each stage of cleaning module completes cleaning of all upper triangle elements in the eigenvalue transition matrix and updating of the eigenvector transition matrix.
The dual-threshold information source number estimation module in the step 3 is used for judging the information source number and obtaining a noise subspace, and comprises a characteristic value ordering module, a threshold judgment module and a noise subspace solving module; the characteristic value ordering module is used for ordering the input characteristic values, the threshold value judging module is used for judging the number of the information sources, and the noise subspace solving module is used for processing and screening the input characteristic vector matrix through the ordered and information source number estimation result to obtain a noise subspace matrix.
The spectrum peak searching module in the step 4 comprises a guide vector calculating module, a space spectrum function calculating module, an extremum checking module and a peak value judging module; the guiding vector calculating module inputs the information source frequency, the array element position coordinates and the sin (theta) value stored in the ROM and is used for calculating the guiding vector value, and the calculating formula is as follows:
wherein θ is the search angle, d i F is the frequency of an incident information source, and c is the speed of light; the spatial spectrum function calculation module is input into a noise subspace and a steering vector, and the spectrum function calculation formula is as follows:
wherein Q is the number of sources, v i An ith column vector in the noise subspace, and a (theta) is a guide vector; the extremum checking module is used for judging the minimum value point in the spectrum function and outputting the extremum point to the peak value judging module, the peak value judging module judges the extremum point and the stored value, if the extremum point is smaller than the stored value, the stored value is updated, otherwise, the stored value is kept unchanged, and the stored value obtained after judging and comparing all the extremum points is the direction of the information source.
Each stage of cleaning module comprises an input control module, a single cleaning module and an output control module; the input control module is used for controlling the position adjustment of the characteristic value transition matrix elements and the input of the characteristic vector transition matrix; the single cleaning module cleans elements at a certain position in the eigenvalue transition matrix in parallel, and calculates to obtain an updated eigenvalue transition matrix and an updated eigenvector transition matrix; the output control module is used for transmitting the operation result to the input control module after single cleaning and controlling the output of data after the primary cleaning is finished.
The single cleaning module comprises an input buffer module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module and an output buffer module; the rotation matrix calculation module is used for calculating the rotation angle of a certain position element in the eigenvalue transition matrix and transmitting the rotation matrix to the eigenvalue transition matrix calculation module and the eigenvector transition matrix calculation module, the eigenvalue transition matrix calculation module multiplies the input cached eigenvalue transition matrix with the rotation matrix to obtain an updated eigenvalue transition matrix and transmits the updated eigenvalue transition matrix to the output cache module, the eigenvector transition matrix calculation module multiplies the input cached eigenvector transition matrix with the rotation matrix to obtain an updated eigenvector transition matrix and transmits the updated eigenvector transition matrix to the output cache module, the rotation matrix calculation module comprises 5 identical calculation units and is used for calculating sine and cosine values of rotation angles, and the calculation formula of the rotation angles is as follows:
taking the first set S1 as an example, the rotation matrix is constructed using the obtained rotation angles as follows:
the eigenvalue matrix update calculation formula is therefore:
R'＝WR'W T
the feature vector moment update calculation formula is:
V'＝WV'
wherein the initial value of R' is a real symmetric matrix
The initial value of V' is the identity matrix.
The characteristic value sequencing module comprises N identical comparison units, wherein N is the number of characteristic values, N paths of parallel are adopted to realize comparison between every two characteristic values, one port of each comparison unit inputs a fixed characteristic value a, the other port sequentially inputs different characteristic values b, and if a is less than b, a result +1 is output; otherwise, if a=b, comparing the position marks of the characteristic values a and b, and outputting a result +1 before the position of a is the position of b; the position of a is not before the position of b, and the output result remains unchanged.
The invention has the beneficial effects that:
(1) The self-adaptive covariance calculation module calculates according to the snapshot number, can realize the calculation of the autocorrelation matrix of different snapshot sampling data, is beneficial to the realization of running water and shortens the operation time;
(2) In the invention, the DOA estimation of the coherent information source is realized by adopting the space smoothing module, and the real number module converts complex numbers into real numbers, so that the calculated amount of the multistage pipelining characteristic decomposition module is reduced;
(3) The multistage characteristic value decomposition module adopts multistage cleaning, so that the realization of running water is facilitated, the operation time is shortened, and the resource and the speed are well optimized by adopting a parallel Jacobian calculation method in each stage of cleaning;
(4) The sorting mode in the dual-threshold information source number estimation module adopts a grouping parallel and serial data input comparison mode to sort, so that the sorting method has the advantages of increasing the operation period compared with parallel full sorting, reducing the occupation of resources and well optimizing the resources and the speed.
Drawings
FIG. 1 is a schematic diagram of a centrally symmetric non-uniform array structure;
FIG. 2 is a general block diagram of the present invention;
FIG. 3 is a schematic diagram of an adaptive covariance matrix calculation module;
FIG. 4 is a schematic diagram of an autocorrelation matrix calculation base module;
FIG. 5 is a schematic diagram of a multi-stage pipeline feature decomposition module;
FIG. 6 is a schematic diagram of a dual threshold source number estimation module;
FIG. 7 is a schematic diagram of a spectral peak search module;
fig. 8 is a model sim simulation result diagram.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
The implementation method of the coherent source DOA estimation FPGA based on the non-uniform array comprises a self-adaptive covariance matrix calculation module, a multi-stage flow characteristic decomposition module, a double-threshold source number estimation module and a spectrum peak search module;
the self-adaptive covariance matrix calculation module comprises an autocorrelation matrix calculation module, a smoothing processing module and a real-realization processing module; the autocorrelation matrix calculation module adopts multipath parallel calculation to obtain the calculation result of each snapshot, and the calculation result of each snapshot is accumulated to obtain the triangular element value on the autocorrelation matrix; the smoothing processing module performs summation operation on the data of the corresponding position of the autocorrelation matrix according to a space smoothing principle; and the real number processing module performs negation operation on the smoothed partial result.
The multistage pipelining characteristic decomposition module is used for calculating characteristic values and characteristic vectors of a covariance matrix and comprises a multistage cleaning module, wherein the input of the multistage cleaning module is the covariance matrix, and the output of the multistage cleaning module is data of main diagonal lines in a characteristic value transition matrix and a characteristic vector matrix; each stage of cleaning module completes cleaning of all upper triangle elements in the characteristic value transition matrix and updating of the characteristic vector transition matrix;
each stage of cleaning module comprises an input control module, a single cleaning module and an output control module; the input control module is used for controlling the position adjustment of the characteristic value transition matrix elements and the input of the characteristic vector transition matrix; the single cleaning module cleans elements at specific positions in the characteristic value transition matrix in parallel, and calculates to obtain an updated characteristic value transition matrix and an updated characteristic vector transition matrix; the output control module is used for transmitting the operation result to the input control module after single cleaning and controlling the output of data after the primary cleaning is finished.
The single cleaning module comprises an input buffer module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module and an output buffer module; the rotation matrix calculation module is used for calculating the rotation angle of a specific position element in the eigenvalue transition matrix, transmitting the rotation matrix to the eigenvalue transition matrix calculation module and the eigenvector transition matrix calculation module, multiplying the eigenvalue transition matrix of the input buffer by the rotation matrix to obtain an updated eigenvalue transition matrix, transmitting the updated eigenvalue transition matrix to the output buffer module, and multiplying the eigenvector transition matrix of the input buffer by the rotation matrix to obtain an updated eigenvector transition matrix and transmitting the updated eigenvector transition matrix to the output buffer module.
The double-threshold information source number estimation module is used for judging the information source number and obtaining a noise subspace, and comprises a characteristic value ordering module, a threshold judgment module and a noise subspace solving module; the characteristic value ordering module is used for ordering the input characteristic values, the threshold value judging module is used for judging the number of the information sources, and the noise subspace solving module is used for processing and screening the input characteristic vector matrix through the ordered estimation results of the number of the information sources to obtain a noise subspace matrix;
the characteristic value sequencing module comprises N identical comparison units, wherein N is the number of characteristic values, N paths of parallel are adopted to realize the comparison between every two characteristic values, one port of each comparison unit inputs a fixed characteristic value a, the other port sequentially inputs different characteristic values b, and if a is smaller than b, a result +1 is output; otherwise, if a=b, comparing the position marks of the characteristic values a and b, and outputting a result +1 before the position of a is the position of b; the position of a is not before the position of b, and the output result remains unchanged.
The spectrum peak searching module comprises a guide vector calculating module, a space spectrum function calculating module, an extremum checking module and a peak value judging module; the guiding vector calculating module inputs the information source frequency, the array element position coordinates and the sin (theta) value stored in the ROM and is used for calculating the guiding vector value, and the calculating formula is as follows:
where θ is the search angle, d i F is the frequency of an incident information source, and c is the speed of light; the spatial spectrum function calculation module is input into a noise subspace and a steering vector, and the spectrum function calculation formula is as follows:
wherein Q is the number of sources, v i An ith column vector in the noise subspace, and a (theta) is a guide vector; the extremum checking module is used for judging the minimum value point in the spectrum function and outputting the extremum point to the peak value judging module, the peak value judging module judges the extremum point and the stored value, if the extremum point is smaller than the stored value, the stored value is updated, otherwise, the stored value is kept unchanged, and the stored value obtained after judging and comparing all the extremum points is the direction of the information source.
DOA estimation is performed by adopting a centrosymmetric non-uniform array shown in figure 1, wherein the array element number is 6, the incident information source number is 2, the incident information sources are related or coherent, and the array element spacing is integer multiple of half wavelength, namely [ d ] 1 ,d 2 ,d 3 ]＝λ 0 /2·[k 1 ,k 2 ,k 3], wherein k1 ,k 2 ,k 3 Is a positive integer and is mutually equal to each other; the angular ambiguity can be effectively restrained by utilizing the array element position design, and in addition, the array structure can adopt a forward and backward space smoothing method to carry out decorrelation, and the smoothing mode is shown as an arrow in the attached figure 1.
An embodiment of a DOA estimation hardware implementation device based on a non-uniform array coherent source, as shown in fig. 2, comprises a self-adaptive covariance matrix calculation module, a multi-stage feature decomposition module, a dual-threshold source number estimation module and a spectrum peak search module.
(1) Self-adaptive covariance matrix calculation module
As shown in fig. 3, the adaptive covariance matrix calculation module comprises an autocorrelation matrix calculation module, a smoothing processing module and a real-realization processing module; the external M paths of sampling data are input into an autocorrelation matrix calculation module according to the snapshot number, and after the received data autocorrelation matrix with the snapshot number of L is calculated, the smoothing processing module and the real number processing module can process the autocorrelation matrix. (1.1) autocorrelation matrix calculation Module
The autocorrelation matrix of the received data is calculated as r=e { XX H And X is an input received data matrix, and since the calculation result R is a hermitian matrix, only the values of the upper 21 triangular elements of the matrix R need to be calculated, and the autocorrelation matrix is calculated as follows:
wherein L is the snapshot number, x i x i H The formula of (2) is as follows:
wherein 
(p, q=1,) 6. R is not used when the forward and backward space smoothing is adopted as shown in figure 1 16 The autocorrelation matrix calculation module therefore comprises 20 identical calculation units mult_accum1 to mult_accum20, each calculation unit comprising 1 complex multiplier cmult and 2 accumulators accumulate1 and accumulate2; the input of the complex multiplier is corresponding sampling data, the accumulator accumulates the single snapshot result, and the calculation result of the autocorrelation matrix is obtained by controlling the input snapshot number L.
(1.2) smoothing processing Module
The forward and backward smoothing is performed in the manner shown in fig. 1, i.e. the autocorrelation matrix obtained by the calculation is used
The elements of the corresponding positions in the table are summed, and the calculation formula is as follows:
wherein rpq (p, q=1, 2,., 6) is an autocorrelation matrix
Elements of corresponding positions in (b) due to r pq Complex, so that every two numbers add requiring two adders; as can be seen from the above, the same elements exist, so that 11 complex additions are required to realize smoothing, namely 22 adders add_1 to add_22 are required to calculate; the smoothing module comprises 22 adders and performs parallel calculation to obtain a smoothed matrix +.>
(1.3) real-ization processing Module
The real number module converts the covariance matrix of the signal from a complex number domain to a real number domain, so that the subsequent feature decomposition is convenient to realize, the subsequent calculation is simplified, ru=lambdau can be obtained according to the property of the feature decomposition, wherein R is the covariance matrix of the signal, lambdais the feature value of the covariance matrix, u is the feature vector of the covariance matrix, R is a complex number matrix and can be expressed as R=R r +iR i ，R r Representing the real part of matrix R, R i Representing the imaginary part of matrix R, u may be expressed as u=u r +iu i ，u r Representing the real part of the vector u, u i The imaginary part representing the vector u, the carry-over equation is available:
thus smoothing the 5 x 5 matrix
According to->
The complex matrix can be converted into a real matrix by conversion, and the number of the characteristic values obtained by carrying out characteristic decomposition on the complex matrix is 10 and is approximately equal to each other; the number of the feature vectors is 10, and the feature vectors are related in pairs; therefore, the real number processing module comprises 11 opposite number computing modules not_1-not_11, and obtains opposite numbers for the imaginary parts of the triangular elements on the smoothed matrix and outputs the computing results of the triangular elements on the covariance matrix to be decomposed.
(2) Multistage running water characteristic decomposition module
As shown in figure 5, the multistage running water feature decomposition module comprises 3 stages of cleaning modules, each stage of cleaning module completes Jacobian rotation of triangle elements on the whole covariance matrix, obtains a corresponding feature vector matrix, obtains non-diagonal elements after three stages of cleaning to be basically zeroed, obtains feature values and feature vectors of the covariance matrix,
in each level of cleaning, the upper triangle of the covariance matrix is divided into 9 groups, each group contains 5 elements, the cleaning of each group of elements is single cleaning, serial operation processing is adopted between groups, the 5 elements in each group are processed by parallel jacobian rotation, and the specific distribution rule of each group is shown in table 1:
table 1 distribution rules of each group
| S1 | S2 | S3 | S4 | S5 | S6 | S7 | S8 | S9 |
| (1,2) | (1,4) | (1,6) | (1,8) | (1,10) | (1,9) | (1,7) | (1,5) | (1,3) |
| (3,4) | (2,6) | (4,8) | (6,10) | (8,9) | (10,7) | (9,5) | (7,3) | (5,2) |
| (5,6) | (3,8) | (2,10) | (4,9) | (6,7) | (8,5) | (10,3) | (9,2) | (7,4) |
| (7,8) | (5,10) | (3,9) | (2,7) | (4,5) | (6,3) | (8,2) | (10,4) | (9,6) |
| (9,10) | (7,9) | (5,7) | (3,5) | (2,3) | (4,2) | (6,4) | (8,6) | (10,8) |
The single cleaning module comprises an input buffer module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module,
the rotation matrix calculating module comprises 5 identical calculating units and is used for calculating sine and cosine values of a rotation angle, and the calculation formula of the rotation angle is as follows:
taking the first set S1 as an example, the rotation matrix shown below is constructed using the obtained rotation angles:
the eigenvalue matrix update calculation formula is therefore:
R'＝WR'W T
the feature vector moment update calculation formula is:
V'＝WV'
wherein the initial value of R' is a real symmetric matrix
The initial value of V' is the identity matrix.
(3) Dual-threshold source number estimation module
As shown in fig. 6, the source number estimation module includes a feature value ordering module, a threshold value judging module, and a noise subspace solving module; storing 10 eigenvalues in memory data_r, and storing r 1 ～r 10 Ordering in order of from big to small, and sorting the ordered result num 1 ～num a Stored in memory sort_result; the characteristic value ordering module sort comprises 10 identical comparison units, and is exemplified by sort_one_1, and one port input r 1 Another port inputs r k (k=1, 2., (a), 10), r 1 And r k Comparing the sizes of (3) if r 1 ＜r k Result=result+1, the initial value of result is 0; if r 1 ＝r k Comparing the position marks, namely the subscript values, and if 1 is less than k, result=result+1; the result value remains unchanged under other conditions. r is (r) 1 The result value after the comparison with 10 feature values is the position coordinate of the result value in the ordered feature values and is stored in num 1 Is a kind of medium.
Threshold value judgment will be performed
Respectively withThe threshold values th1 and th2 are compared, and if the threshold value th1 is larger than the threshold value th2, the source is determined to be the source, if the threshold value th1 is smaller than the threshold value th2, the source is likely to be the source, and if the threshold value th2 is smaller than the threshold value th2, the source is not the source.
(4) Spectral peak search module
As shown in fig. 7, the spectral peak search includes a steering vector calculation module, a spatial spectral function calculation module, an extremum inspection module and a peak value judgment module; firstly, calculating the value of a guide vector, and simplifying the value by utilizing the position coordinate setting of an array element, wherein the calculation formula is as follows:
wherein ki (i=1, 2,.,. 5) is a positive integer, f is the incident source frequency, f 0 For the highest frequency of the incident source, θ is the search angle, the search step is 0.1 °, and since the incident angle ranges from 0 ° to 60 °, the values of sin (-0.5 °) -sin (60.5 °) are stored in the ROM in search steps;
then, calculating a spatial spectrum function, simplifying the inversion operation, and equivalently, the position of a peak of the search spectrum function is equivalent to the position of a search trough, and as the noise subspace matrix is real, the square is simplified to the absolute value operation, and the calculation formula of the spectrum function is as follows:
the calculated spatial spectrum value P (theta i ) Respectively with P (theta) i-1) and P(θi+1 ) By comparison, if P (θ i )＜P(θ i-1 ) And P (theta) i )＜P(θ i+1 ) At the same time satisfy then P (θ) i ) Is an extreme point; comparing the extreme point with the smaller two values of the extreme point stored in the memory, and if the value is smaller than the value stored in the register, then P (θ i ) If not, waiting for the input of the next extreme point, and after all search angle calculation is completed, storing the corresponding reference number of the peak value stored in the registerAnd converting to obtain the DOA estimated angle.
FIG. 8 is a diagram of simulation results of a model, in which a Qurtus platform is used to implement FPGA implementation on a board of model Stratix10x, simulation is performed using model sim, the incident source is 2, the incident angle is 15 degrees, 10 degrees, the simulation results show that the value of aoa_1 is 150, the value of aoa_2 is 100, the product of the result and the search step is the correct DOA estimation angle, the required period is 3415 periods, 21.4us is required when 160MHz clock is used, pipelined calculation is used, the second set of results only requires 905 periods, and 5.7us is required when 160MHz clock is used. The number of occupied resources is shown in table 2:
TABLE 2 number of resources occupied
Claims (7)
1. The implementation method of the coherent source DOA estimation FPGA based on the non-uniform array is characterized by comprising the following steps:
step 1: inputting the received signals into an adaptive covariance matrix calculation module, and calculating to obtain an adaptive covariance matrix of the received signals;
the self-adaptive covariance matrix calculation module comprises an autocorrelation matrix calculation module, a smoothing processing module and a real processing module;
the autocorrelation matrix calculation module adopts multipath parallel calculation to obtain the calculation result of each snapshot, adds up the calculation result of each snapshot to obtain the triangle element value on the autocorrelation matrix, and the autocorrelation matrix calculation formula of the received data is as follows:
wherein the method comprises the steps ofL is the snapshot number;
the calculation formula of (2) is as follows:
wherein ,
the autocorrelation matrix calculation module comprises a complex multiplier and an accumulator; the input of the complex multiplier is corresponding sampling data, an accumulator accumulates the single snapshot result, and the input snapshot number L is used for controlling to obtain the calculation result of the autocorrelation matrix;
the smoothing processing module performs a smoothing operation on the autocorrelation matrix according to the spatial smoothing principle
The elements of the corresponding positions in the table are summed, and the calculation formula is as follows:
wherein ,rpq Is an autocorrelation matrix
Elements of the corresponding position in (a);
the real-number processing module converts the covariance matrix of the signal from a complex number domain to a real number domain; deriving ru=λu from the property of the eigen decomposition, where λ is the eigenvalue of the covariance matrix, R is the covariance matrix of the signal and is a complex matrix, denoted r=r r +iR i ，R r Representing the real part of matrix R, R i The imaginary part of the matrix R, u is the eigenvector of the covariance matrix, u is denoted u=u r +iu i ，u r Representing vectorsThe real part of u, u i An imaginary part representing the vector u;
matrix is formed
According to->
Converting the complex matrix into a real matrix, solving the opposite number of the imaginary part of the triangular element on the smoothed matrix, and outputting the calculation result of the triangular element on the covariance matrix to be decomposed;
step 2: inputting the obtained self-adaptive covariance matrix into a multi-stage running water characteristic decomposition module to decompose the multi-running water characteristic values, so as to obtain characteristic values and characteristic vectors of the covariance matrix;
step 3: inputting the obtained characteristic value and the characteristic vector into a double-threshold information source number estimation module, judging the information source number and obtaining a noise subspace;
step 4: and outputting the double-threshold information source number estimation module, inputting the output of the double-threshold information source number estimation module into a spectrum peak searching module, and calculating to obtain a stored value obtained after the spectrum function judges and compares all extreme points, namely the direction in which the information source is located.
2. The non-uniform array-based coherent source DOA estimation FPGA implementation method of claim 1, wherein: the multistage pipelining feature decomposition module in step 2 is used for calculating feature values and feature vectors of a covariance matrix, and comprises a multistage cleaning module, wherein the input of the multistage cleaning module is the covariance matrix, and the output of the multistage cleaning module is data of main diagonal lines in a feature value transition matrix and a feature vector matrix; and each stage of cleaning module completes cleaning of all upper triangle elements in the eigenvalue transition matrix and updating of the eigenvector transition matrix.
3. The non-uniform array-based coherent source DOA estimation FPGA implementation method of claim 1, wherein: the dual-threshold information source number estimation module in the step 3 is used for judging the information source number and obtaining a noise subspace, and comprises a characteristic value ordering module, a threshold judgment module and a noise subspace solving module; the characteristic value ordering module is used for ordering the input characteristic values, the threshold value judging module is used for judging the number of the information sources, and the noise subspace solving module is used for processing and screening the input characteristic vector matrix through the ordered and information source number estimation result to obtain a noise subspace matrix.
4. The non-uniform array-based coherent source DOA estimation FPGA implementation method of claim 1, wherein: the spectrum peak searching module in the step 4 comprises a guide vector calculating module, a space spectrum function calculating module, an extremum checking module and a peak value judging module; the guiding vector calculating module inputs the information source frequency, the array element position coordinates and the sin (theta) value stored in the ROM and is used for calculating the guiding vector value, and the calculating formula is as follows:
wherein θ is the search angle, d i F is the frequency of an incident information source, and c is the speed of light; the spatial spectrum function calculation module is input into a noise subspace and a steering vector, and the spectrum function calculation formula is as follows:
wherein Q is the number of sources, v i An ith column vector in the noise subspace, and a (theta) is a guide vector; the extremum checking module is used for judging the minimum value point in the spectrum function and outputting the extremum point to the peak value judging module, the peak value judging module judges the extremum point and the stored value, if the extremum point is smaller than the stored value, the stored value is updated, otherwise, the stored value is kept unchanged, and the obtained extremum points are judged and comparedThe stored value is the direction in which the source is located.
5. The non-uniform array-based coherent source DOA estimation FPGA implementation method of claim 1, wherein: each stage of cleaning module comprises an input control module, a single cleaning module and an output control module; the input control module is used for controlling the position adjustment of the characteristic value transition matrix elements and the input of the characteristic vector transition matrix; the single cleaning module cleans elements at a certain position in the eigenvalue transition matrix in parallel, and calculates to obtain an updated eigenvalue transition matrix and an updated eigenvector transition matrix; the output control module is used for transmitting the operation result to the input control module after single cleaning and controlling the output of data after the primary cleaning is finished.
6. The non-uniform array-based coherent source DOA estimation FPGA implementation method of claim 5, wherein the method comprises the steps of: the single cleaning module comprises an input buffer module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module and an output buffer module; the rotation matrix calculation module is used for calculating the rotation angle of a certain position element in the eigenvalue transition matrix and transmitting the rotation matrix to the eigenvalue transition matrix calculation module and the eigenvector transition matrix calculation module, the eigenvalue transition matrix calculation module multiplies the input cached eigenvalue transition matrix with the rotation matrix to obtain an updated eigenvalue transition matrix and transmits the updated eigenvalue transition matrix to the output cache module, the eigenvector transition matrix calculation module multiplies the input cached eigenvector transition matrix with the rotation matrix to obtain an updated eigenvector transition matrix and transmits the updated eigenvector transition matrix to the output cache module, the rotation matrix calculation module comprises 5 identical calculation units and is used for calculating sine and cosine values of rotation angles, and the calculation formula of the rotation angles is as follows:
taking the first set S1 as an example, the rotation matrix is constructed using the obtained rotation angles as follows:
the eigenvalue matrix update calculation formula is therefore:
R′＝WR′W T
the feature vector moment update calculation formula is:
V′＝WV′
wherein the initial value of R' is a real symmetric matrix
The initial value of V' is the identity matrix.
7. The method for implementing the non-uniform array-based coherent source DOA estimation FPGA according to claim 3, wherein the method comprises the following steps: the characteristic value sequencing module comprises N identical comparison units, wherein N is the number of characteristic values, N paths of parallel are adopted to realize comparison between every two characteristic values, one port of each comparison unit inputs a fixed characteristic value a, the other port sequentially inputs different characteristic values b, and if a is less than b, a result +1 is output; otherwise, if a=b, comparing the position marks of the characteristic values a and b, and outputting a result +1 before the position of a is the position of b; the position of a is not before the position of b, and the output result remains unchanged.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910669180.5A CN110361691B (en) | 2019-07-24 | 2019-07-24 | Implementation method of coherent source DOA estimation FPGA based on non-uniform array |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910669180.5A CN110361691B (en) | 2019-07-24 | 2019-07-24 | Implementation method of coherent source DOA estimation FPGA based on non-uniform array |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110361691A CN110361691A (en) | 2019-10-22 |
| CN110361691B true CN110361691B (en) | 2023-05-30 |
Family
ID=68220117
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910669180.5A Active CN110361691B (en) | 2019-07-24 | 2019-07-24 | Implementation method of coherent source DOA estimation FPGA based on non-uniform array |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110361691B (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112543025B (en) * | 2020-12-08 | 2023-03-14 | 成都天奥信息科技有限公司 | High-speed serial AD sampling and data processing system and method based on matrixing |
| CN114237548B (en) * | 2021-11-22 | 2023-07-18 | 南京大学 | Method and system for complex point multiplication operation based on nonvolatile memory array |
| CN114325579B (en) * | 2022-03-09 | 2022-07-15 | 网络通信与安全紫金山实验室 | Positioning parameter estimation method, positioning parameter estimation device, positioning parameter estimation apparatus, storage medium, and program product |
| CN114545324B (en) * | 2022-04-24 | 2022-07-12 | 南京宇安防务科技有限公司 | Rapid direction finding method suitable for non-uniform array |
| CN117192471A (en) * | 2023-09-06 | 2023-12-08 | 无锡芯光互连技术研究院有限公司 | HLS-based two-dimensional DOA estimation method, device and storage medium |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104052701A (en) * | 2014-06-03 | 2014-09-17 | 哈尔滨工程大学 | Intra-pulse modulation characteristic real-time extraction and classification system based on FPGA |
| CN105703083A (en) * | 2016-04-26 | 2016-06-22 | 深圳前海智讯中联科技有限公司 | Multi-beam selection smart antenna array and system having antenna array |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101272168B (en) * | 2007-03-23 | 2012-08-15 | 中国科学院声学研究所 | Signal sources estimation method and its DOA estimation method |
-
2019
- 2019-07-24 CN CN201910669180.5A patent/CN110361691B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104052701A (en) * | 2014-06-03 | 2014-09-17 | 哈尔滨工程大学 | Intra-pulse modulation characteristic real-time extraction and classification system based on FPGA |
| CN105703083A (en) * | 2016-04-26 | 2016-06-22 | 深圳前海智讯中联科技有限公司 | Multi-beam selection smart antenna array and system having antenna array |
Non-Patent Citations (1)
| Title |
|---|
| 二维高精度MUSIC算法的高速实现;张多利等;《合肥工业大学学报（自然科学版）》;20180331;第41卷(第3期);论文摘要，第1部分 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110361691A (en) | 2019-10-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110361691B (en) | Implementation method of coherent source DOA estimation FPGA based on non-uniform array | |
| CN106980106B (en) | Sparse DOA estimation method under array element mutual coupling | |
| CN107589399B (en) | Estimation method of direction of arrival of co-prime array based on singular value decomposition of multi-sampling virtual signal | |
| CN106019215B (en) | Nested array direction of arrival angle method of estimation based on fourth-order cumulant | |
| CN110222307B (en) | Parallel implementation method for eigenvalue decomposition of real symmetric matrix based on FPGA | |
| CN104793176A (en) | FPGA (field programmable gate array) based DOA (direction of arrival) estimation fast-implementing method | |
| CN111413668B (en) | DOA estimation method based on DFT enhancement in large-scale array | |
| CN103516643A (en) | MIMO detecting preprocessing device and method | |
| CN114488034A (en) | Passive detection and interference reconnaissance integrated device and method | |
| Li et al. | Hardware acceleration of MUSIC algorithm for sparse arrays and uniform linear arrays | |
| CN103837858A (en) | Far field direction of arrival estimation method applied to plane array and system thereof | |
| Chen et al. | Fast algorithm for DOA estimation with partial covariance matrix and without eigendecomposition | |
| Feng et al. | An off-grid iterative reweighted approach to one-bit direction of arrival estimation | |
| Huang et al. | An efficient FPGA implementation for 2-D MUSIC algorithm | |
| CN105608057A (en) | FPGA realization module and FPGA realization method for signal subspace decomposition by time-sharing multiplexing of hardware resources | |
| CN113203997B (en) | FPGA-based radar super-resolution direction finding method, system and application | |
| CN109507634B (en) | Blind far-field signal direction-of-arrival estimation method based on propagation operator under any sensor array | |
| CN108008665B (en) | Large-scale circular array real-time beam former based on single-chip FPGA and beam forming calculation method | |
| CN109683128B (en) | Single-snapshot direction finding method under impact noise environment | |
| Jiménez | Smart cities, open innovation and open government: towards | |
| CN111865385B (en) | Two-dimensional planar array digital beam forming method based on FPGA | |
| CN115079090A (en) | Multi-array non-circular source rapid positioning method based on ESPRIT and weighted dimension reduction search | |
| CN114217265A (en) | Minimum variance distortionless response-based source DOA estimation method and system | |
| CN108415005A (en) | A kind of passive location delay time estimation method and device | |
| CN114488064A (en) | Distance-speed joint estimation method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |