CN110361691A - Coherent DOA based on nonuniform noise estimates FPGA implementation method - Google Patents
Coherent DOA based on nonuniform noise estimates FPGA implementation method Download PDFInfo
- Publication number
- CN110361691A CN110361691A CN201910669180.5A CN201910669180A CN110361691A CN 110361691 A CN110361691 A CN 110361691A CN 201910669180 A CN201910669180 A CN 201910669180A CN 110361691 A CN110361691 A CN 110361691A
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 22
- 230000001427 coherent effect Effects 0.000 title claims abstract description 21
- 239000011159 matrix material Substances 0.000 claims abstract description 211
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 27
- 238000012545 processing Methods 0.000 claims abstract description 23
- 238000001228 spectrum Methods 0.000 claims abstract description 23
- 230000003044 adaptive effect Effects 0.000 claims abstract description 10
- 238000004364 calculation method Methods 0.000 claims description 127
- 230000007704 transition Effects 0.000 claims description 65
- 239000013598 vector Substances 0.000 claims description 42
- 238000004140 cleaning Methods 0.000 claims description 41
- 238000009499 grossing Methods 0.000 claims description 20
- 230000003595 spectral effect Effects 0.000 claims description 12
- 238000012163 sequencing technique Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000013139 quantization Methods 0.000 claims description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 abstract description 6
- 238000013461 design Methods 0.000 abstract description 4
- 230000009977 dual effect Effects 0.000 abstract description 4
- 238000005516 engineering process Methods 0.000 abstract 1
- 230000008520 organization Effects 0.000 abstract 1
- 230000006870 function Effects 0.000 description 14
- 238000010586 diagram Methods 0.000 description 9
- 238000004088 simulation Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 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
- 238000005259 measurement Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
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 array signal processing technologies, and in particular to a kind of coherent DOA estimation FPGA implementation method based on nonuniform noise.Including adaptive covariance matrix computing module, multistage flowing water feature decomposition module, dual threshold Sources number estimation module and spectrum peak search module；The input of adaptive covariance matrix computing module is sampled data, and output is connected with multistage flowing water feature decomposition module input；Multistage flowing water feature decomposition module output is connected with the input of dual threshold Sources number estimation module, and the output of dual threshold Sources number estimation module is connected with spectrum peak search module；The DOA estimation of central symmetry nonuniform noise coherent may be implemented in the present invention, can use different data formats, be designed on demand, calculates covariance matrix and parallel Jacobi numerical procedure using by number of snapshots, shortens operation time；Using the design philosophy of pipeline organization under the premise of guaranteeing operational precision, resource and speed are optimized to the maximum extent.
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 arrival is one of important branches of array signal processing, and has wide application prospects in a plurality of fields such as radar, communication, sonar and the like. The uniform array has the defects of small array aperture, low direction finding precision, poor resolution and no fuzzy direction finding condition that the minimum array element distance is smaller than the half wavelength of an incident information source due to the limitation of the physical size of an antenna in practical application, and is difficult to be applied to the modern complicated electromagnetic environment. The non-uniform array can possess bigger aperture than the even array of the same array element number, higher resolution and higher degree of freedom, and reasonable setting array structure can avoid the direction finding fuzzy and can improve the direction finding performance, adopts symmetrical array element setting or array element interval combination to set up the array and not only can enlarge the array aperture, improves the angle finding precision and can also utilize the space smoothness to solve coherently, improves the coherent ability of solution of non-uniform array.
The MUSIC algorithm is mainly divided into four parts, namely covariance calculation, characteristic decomposition, information source number estimation and spectral peak search, a large amount of calculation is needed for characteristic decomposition and spectral peak search, the time needed by the adopted DSP to realize DOA estimation hardly meets the requirement of real-time direction finding at the ms level, the DOA estimation parallelism is high by adopting the FPGA, the running water can be adopted to realize, the time needed by the DOA estimation is greatly reduced, and the us-level DOA estimation is realized. In the existing research of the MUSIC algorithm, a grouping scheme which is executed in parallel is adopted in a covariance calculation part, so that a control circuit is complex, fast-beat-number variable calculation cannot be realized, a corresponding coherent module solving module design is not provided, and the consideration on non-centrosymmetric array real number is less. The characteristic decomposition part comprises a large amount of nonlinear operation, the existing method adopts serial Jacobi operation to realize long time consumption, adopts parallel Jacobi algorithm to realize characteristic decomposition, and reduces the time consumption of characteristic decomposition. The spectral peak searching part usually adopts a memory circuit to obtain the guide vector, and calculates the spatial spectrum function through directly read values, thereby consuming a large amount of storage resources.
Disclosure of Invention
The invention aims to provide a coherent information source DOA estimation FPGA implementation method based on a non-uniform array, which is suitable for the non-uniform array capable of smoothly resolving coherence, and the hardware design implemented on an FPGA platform verifies the feasibility.
The purpose of the invention is realized as follows:
a coherent source DOA estimation FPGA realization method based on a non-uniform array comprises the following steps:
step 1: inputting the received signal into a self-adaptive covariance matrix calculation module, and calculating to obtain a self-adaptive covariance matrix of the received signal;
step 2: inputting the obtained adaptive covariance matrix into a multi-stage flow characteristic decomposition module, and performing multi-flow characteristic value decomposition to obtain a characteristic value and a characteristic vector of the covariance matrix;
and step 3: inputting the obtained characteristic values and the characteristic vectors into a dual-threshold signal source number estimation module, judging the number of signal sources and obtaining a noise subspace;
and 4, step 4: and inputting the output of the dual-threshold information source number estimation module into a spectral peak search module, and calculating to obtain a spectral function, and judging and comparing all extreme points to obtain a stored value, namely the direction of the information source.
The adaptive covariance matrix calculation module in the step 1 comprises an autocorrelation matrix calculation module, a smoothing processing module and a real quantization processing module;
the autocorrelation matrix calculation module adopts multi-path parallel calculation to obtain the calculation result of each snapshot, and accumulates the calculation result of each snapshot to obtain the triangular 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 fast beats,the calculation formula of (A) 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, the accumulator accumulates the result of single snapshot and controls by using the input snapshot number L to obtain the calculation result of the autocorrelation matrix;
the smoothing processing module carries out summation operation on data of corresponding positions of the autocorrelation matrix according to the space smoothing principle, and carries out summation operation on the autocorrelation matrix obtained by the calculationThe elements at the corresponding positions are summed, and the calculation formula is:
wherein r ispq(p, q ═ 1, 2.., 6) is an autocorrelation matrixElements of the corresponding position in;
and the real-number processing module is used for carrying out inverse operation on the smoothed partial result, converting the covariance matrix of the signal from a complex number domain to a real number domain, and obtaining Ru ═ λ u according to the characteristic decomposition property, wherein λ is the characteristic value of the covariance matrix, and R is the covariance matrix of the signal and is a complex matrix represented by R ═ Rr+iRi，RrRepresenting the real part of a matrix R, RiDenotes the imaginary part of the matrix R, u is the eigenvector of the covariance matrix, and u is expressed as u ═ ur+iui，urRepresenting the real part of the vector u, uiRepresenting the imaginary part of vector u, taken into the above equation:
will matrixAccording toAnd the conversion is to convert the complex matrix into a real matrix, solve the inverse 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 flow characteristic decomposition module is used for calculating the eigenvalue and the eigenvector of the 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 the data and the eigenvector matrix of the main diagonal line in the eigenvalue transition matrix; and each level of cleaning module finishes cleaning all upper triangular elements in the characteristic value transition matrix and updating the characteristic vector transition matrix.
The dual-threshold signal source number estimation module is used for judging the number of signal sources and obtaining a noise subspace and comprises a characteristic value sequencing module, a threshold judgment module and a noise subspace solving module; the eigenvalue sorting module is used for sorting the input eigenvalues, the threshold value judging module is used for judging the number of the information sources, and the noise subspace solving module processes and screens the input eigenvector matrix according to the sorted and information source number estimation results 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 extreme value checking module and a peak value judging module; the guiding vector calculation module inputs information source frequency, array element position coordinates and sin (theta) value stored in a ROM to calculate the value of the guiding vector, and the calculation formula is as follows:
where θ is the search angle, diIs the position coordinate of the array element, f is the frequency of the incident information source, and c is the speed of light; the input of the spatial spectrum function calculation module is a noise subspace and a guide vector, and the spectrum function calculation formula is as follows:
wherein Q is the number of sources, viIs the ith column vector in the noise subspace, and a (theta) is the guide vector; the extreme value checking module is used for judging minimum value points in the spectrum function and outputting the extreme value points to the peak value judging module, the peak value judging module judges the extreme value points and the stored values, if the minimum value points are smaller than the stored values, the stored values are updated, otherwise, the stored values are kept unchanged, and the stored values obtained after all the extreme value points are judged and compared are the direction of the information source.
Each level 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 the updated eigenvalue transition matrix and eigenvector transition matrix are obtained through calculation; 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 finishing the first-stage cleaning.
The single cleaning module comprises an input cache module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module and an output cache 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 same calculation units and is used for calculating the sine and cosine values of the rotation angle, and the calculation formula of the rotation angle is as follows:
taking the first group S1 as an example, the rotation matrix is constructed by obtaining the rotation angles and is expressed as follows:
therefore, the eigenvalue matrix update calculation formula is:
R'＝WR'WT
the eigenvector moment update calculation formula is:
V'＝WV'
wherein the initial value of R' is a real symmetric matrixThe 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 are adopted to realize the comparison between every two characteristic values in parallel, one port of each comparison unit inputs a fixed characteristic value a, the other port inputs different characteristic values b in sequence, and if a is smaller than b, a result +1 is output; otherwise, if a is equal to b, comparing the position marks of the characteristic values a and b, and outputting a result +1 when the position of a is before the position of b; the position of a is not before the position of b, and the output result is kept unchanged.
The invention has the beneficial effects that:
(1) the self-adaptive covariance calculation module calculates according to the snapshot number, can realize the autocorrelation matrix calculation of different snapshot sampling data, is beneficial to realizing pipelining, and shortens the operation time;
(2) according to the method, a spatial smoothing module is adopted to realize DOA estimation of a coherent information source, a real number module converts complex numbers into real numbers, and the calculation amount of a multi-stage flow characteristic decomposition module is reduced;
(3) the multistage eigenvalue decomposition module adopts multistage cleaning, is beneficial to the realization of running water, shortens the operation time, adopts a parallel Jacobian calculation method in each stage of cleaning, and well optimizes resources and speed;
(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, increases the parallel full-sorting operation period, reduces the resource occupation, and well optimizes the resource and the speed.
Drawings
FIG. 1 is a schematic diagram of a centrosymmetric 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 block;
FIG. 4 is a schematic diagram of a basic module for autocorrelation matrix calculation;
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 diagram of modelsim simulation results.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
A coherent information source DOA estimation FPGA implementation method based on a non-uniform array comprises a self-adaptive covariance matrix calculation module, a multi-stage flow characteristic decomposition module, a dual-threshold information source number estimation module and a spectral peak search module;
the self-adaptive covariance matrix calculation module comprises a self-correlation matrix calculation module, a smoothing processing module and a real number processing module; the autocorrelation matrix calculation module adopts multi-path parallel calculation to obtain a calculation result of each snapshot, and the calculation results of each snapshot are accumulated to obtain a triangular element value on the autocorrelation matrix; the smoothing processing module carries out summation operation on data at corresponding positions of the autocorrelation matrix according to a spatial smoothing principle; and the real number processing module performs negation operation on the smoothed partial result.
The multi-stage flow characteristic decomposition module is used for calculating the eigenvalue and the eigenvector of the covariance matrix, and comprises a multi-stage cleaning module, wherein the input of the multi-stage cleaning module is the covariance matrix, and the output of the multi-stage cleaning module is the data of a main diagonal line in the eigenvalue transition matrix and the eigenvector matrix; each level of cleaning module finishes cleaning all upper triangular elements in the characteristic value transition matrix and updating the characteristic vector transition matrix;
each level 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 eigenvalue transition matrix in parallel, and the updated eigenvalue transition matrix and eigenvector transition matrix are obtained through calculation; 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 finishing the first-stage cleaning.
The single cleaning module comprises an input cache module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module and an output cache module; the rotation matrix calculation module is used for calculating the rotation angle of a specific 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 eigenvalue transition matrix of the input cache by the rotation matrix to obtain an updated eigenvalue transition matrix and transmits the updated eigenvalue transition matrix to the output cache module, and the eigenvector transition matrix calculation module multiplies the eigenvector transition matrix of the input cache by the rotation matrix to obtain an updated eigenvector transition matrix and transmits the updated eigenvector transition matrix to the output cache module.
The dual-threshold signal source number estimation module is used for judging the number of signal sources and obtaining a noise subspace and comprises a characteristic value sequencing module, a threshold judgment module and a noise subspace solving module; the eigenvalue sorting module is used for sorting the input eigenvalues, the threshold value judging module is used for judging the number of the information sources, and the noise subspace solving module processes and screens the input eigenvector matrix according to the sorted and information source number estimation results 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 are adopted to realize comparison between every two characteristic values in parallel, one port of each comparison unit inputs a fixed characteristic value a, the other port inputs different characteristic values b in sequence, and if a is smaller than b, a result +1 is output; otherwise, if a is equal to b, comparing the position marks of the characteristic values a and b, and outputting a result +1 when the position of a is before the position of b; the position of a is not before the position of b, and the output result is kept unchanged.
The spectrum peak searching module comprises a guide vector calculating module, a space spectrum function calculating module, an extreme value checking module and a peak value judging module; the guiding vector calculation module inputs information source frequency, array element position coordinates and sin (theta) value stored in a ROM to calculate the value of the guiding vector, and the calculation formula is as follows:
where θ is the search angle, diIs the position coordinate of the array element, f is the frequency of the incident information source, and c is the speed of light; the input of the spatial spectrum function calculation module is a noise subspace and a guide vector, and the spectrum function calculation formula is as follows:
wherein Q is the number of sources, viIs the ith column vector in the noise subspace, and a (theta) is the guide vector; the extreme value checking module is used for judging minimum value points in the spectrum function and outputting the extreme value points to the peak value judging module, the peak value judging module judges the extreme value points and the stored values, if the minimum value points are smaller than the stored values, the stored values are updated, otherwise, the stored values are kept unchanged, and the stored values obtained after all the extreme value points are judged and compared are the direction of the information source.
Adopting the centrosymmetric non-uniform array shown in figure 1 to perform DOA estimation, 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 half-waveInteger multiples of length, i.e. [ d ]1,d2,d3]＝λ0/2·[k1,k2,k3]Wherein k is1,k2,k3Is a positive integer and is relatively prime in pairs; the array element position design can effectively inhibit the angle measurement blurring, and in addition, the array structure can be decorrelated by adopting a forward and backward space smoothing method, wherein the smoothing mode is shown as an arrow in the attached figure 1.
An embodiment of a hardware implementation apparatus based on non-uniform array coherent source DOA estimation is shown in fig. 2, and includes an adaptive covariance matrix calculation module, a multi-stage feature decomposition module, a dual-threshold source number estimation module, and a spectral peak search module.
(1) Adaptive covariance matrix calculation module
As shown in fig. 3, the adaptive covariance matrix calculation module includes an autocorrelation matrix calculation module, a smoothing processing module, and a real quantization processing module; and the external M paths of sampling data are input into an autocorrelation matrix calculation module according to the fast beat number, and after a received data autocorrelation matrix with the fast beat number of L is obtained through calculation, the autocorrelation matrix is processed by a smoothing processing module and a real number processing module. (1.1) autocorrelation matrix calculation Module
The autocorrelation matrix of the received data is calculated as R ═ E { XXHAnd X is an input received data matrix, and since the calculation result R is a hermitian matrix, only values of 21 upper triangular elements of the matrix R need to be calculated, and the autocorrelation matrix calculation formula is in the following form:
wherein L is the number of fast beats, xixi HIs calculated as follows:
wherein(p, q ═ 1.., 6). R is not used when employing the forward and backward spatial smoothing as in FIG. 116The autocorrelation matrix calculation module therefore comprises 20 identical calculation units mult _ accum1 to mult _ accum20, each comprising 1 complex multiplier cmult and 2 accumulators accumulate1 and accumulate 2; the input of the complex multiplier is corresponding sampling data, the accumulator accumulates the result of single snapshot and controls the input snapshot number L to obtain the calculation result of the autocorrelation matrix.
(1.2) smoothing Module
Smoothing the forward and backward directions in the manner as shown in FIG. 1, i.e. on the autocorrelation matrix obtained by the above calculationThe elements at the corresponding positions are summed, and the calculation formula is:
wherein r ispq(p, q ═ 1, 2.., 6) is an autocorrelation matrixElement of middle corresponding position, due to rpqIs complex, so that two adders are needed for adding every two numbers; as can be seen from the above formula, the same elements exist, so that 11 complex additions are required for realizing the smoothing, that is, 22 adders add _1 to add _22 are required for calculation; the smoothing module comprises 22 adders for parallel computing to obtain a smoothed matrix
(1.3) real number 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 implementation of subsequent feature decomposition is facilitated, the subsequent calculation is simplified, and Ru is equal to λ u according to the property of the feature decomposition, whereinR is a covariance matrix of the signal, λ is an eigenvalue of the covariance matrix, u is an eigenvector of the covariance matrix, and R is a complex matrix which can be expressed as R ═ Rr+iRi，RrRepresenting the real part of a matrix R, RiRepresenting the imaginary part of the matrix R, u may be expressed as u-ur+iui，urRepresenting the real part of the vector u, uiRepresenting the imaginary part of vector u, and substituting the equation above yields:
thus smoothing the 5 × 5 matrixAccording toThe conversion can convert the complex matrix into a real matrix, and the number of eigenvalues obtained by performing characteristic decomposition on the real matrix is 10 and the eigenvalues are approximately equal to each other; the number of the feature vectors is also 10, and the feature vectors are related pairwise; therefore, the real number processing module comprises 11 inverse number calculation modules not _1 to not _11, and is used for solving the inverse 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.
(2) Multi-stage flow characteristic decomposition module
As shown in fig. 5, the multi-stage flow characteristic decomposition module includes 3 stages of cleaning modules, each stage of cleaning module completes the Jacobi rotation of the triangular elements on the whole covariance matrix, and obtains the corresponding eigenvector matrix, obtains the basic zero of the off-diagonal elements after three stages of cleaning, obtains the eigenvalue and eigenvector of the covariance matrix,
so the element split of the last triangle of covariance matrix is 9 groups in every grade cleans, and every group contains 5 elements, and the cleaning of every group element is for the single cleaning, adopts serial arithmetic processing between group and the group, and 5 elements in every group adopt the rotation processing of parallel jacobian, and concrete every group distribution rule is as shown in table 1:
TABLE 1 distribution rule 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 cache module, a rotation matrix calculation module, an eigenvalue transition matrix calculation module and an eigenvector transition matrix calculation module,
the rotation matrix calculation module comprises 5 same calculation units and is used for calculating the sine and cosine values of the rotation angle, and the calculation formula of the rotation angle is as follows:
taking the first group S1 as an example, a rotation matrix as shown below is constructed by obtaining the rotation angles:
therefore, the eigenvalue matrix update calculation formula is:
R'＝WR'WT
the eigenvector moment update calculation formula is:
V'＝WV'
wherein the initial value of R' is a real symmetric matrixThe initial value of V' is the identity matrix.
(3) Dual-threshold information source number estimation module
As shown in fig. 6, the information source number estimation module includes a eigenvalue sorting module, a threshold judgment module, and a noise subspace solving module; storing 10 characteristic values in a memory data _ r, and r1～r10Sorting according to the sequence from big to small and sorting the result num1～numaStored in the memory sort _ result; the feature value sorting module sort comprises 10 identical comparison units, which is described by taking sort _ one _1 as an example, and one port inputs r1The other port inputs rk(k 1, 2.., 10), mixing r with r1And rkIs compared with the magnitude of (a), if r1＜rkIf the result is greater than the preset value, then result is equal to result +1, and the initial value of result is 0; if r1＝rkThen carry out index or subscript value on the positionComparing, if 1 is less than k, result is equal to result + 1; the value of result remains unchanged under other conditions. r is1The result value after the comparison with the 10 characteristic values is the position coordinate of the result value in the ordered characteristic values and is stored in num1In (1).
When the threshold value is judged, theCompared with the threshold th1 and the threshold th2, if the threshold th1 is larger than the threshold th, the signal source is determined to be the signal source if the threshold th2 is larger than the threshold th1, and the signal source is not determined to be the signal source if the threshold th2 is smaller than the threshold th 2.
(4) Spectral peak searching module
As shown in fig. 7, the spectral peak search includes a guide vector calculation module, a spatial spectrum function calculation module, an extremum checking module and a peak value judgment module; firstly, the value of the guide vector is calculated, and the value is simplified by utilizing the position coordinate setting of the array element, and the calculation formula is as follows:
wherein k isi(i ═ 1, 2.., 5) is a positive integer, f is the incident source frequency0The maximum frequency of an incident information source is, theta is a search angle, the search step is 0.1 degrees, and because the incident angle range is 0-60 degrees, the values of sin (-0.5 degrees) to sin (60.5 degrees) are stored in the ROM according to the search step;
then, a spatial spectrum function is calculated, the reciprocal calculation is simplified, the position of the peak of the search spectrum function is equivalent to the position of the search trough, since the noise subspace matrix is real number, the square is simplified to the absolute value calculation, and the spectrum function calculation formula is as follows:
the calculated spatial spectrum value P (theta)i) Respectively with P (theta)i-1) And P (theta)i+1) For comparison, if P (θ)i)＜P(θi-1) And P (theta)i)＜P(θi+1) While satisfying the criterion P (theta)i) Is an extreme point; comparing the extreme point with the smaller two values of the extreme points stored in the memory, and if the value is smaller than the value stored in the register, then P (theta) is calculatedi) If not, the input of the next extreme point is waited, and after all the search angles are calculated, the labels corresponding to the peak values stored in the register are converted to obtain the DOA estimation angle.
Fig. 8 is a model simulation result diagram, wherein a Quartus platform is used to complete FPGA implementation on a Stratix10x board, and the modelsim is used to perform simulation, an incident signal source is 2, an incident angle is 15 ° and 10 °, it can be seen from the simulation result 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, when a 160MHz clock is used, 21.4us is needed, running water calculation is adopted, the second group of results only needs 905 periods, and when a 160MHz clock is adopted, 5.7us is needed. The number of occupied resources is shown in table 2:
TABLE 2 number of occupied resources
Claims (8)
1. A coherent source DOA estimation FPGA implementation method based on a non-uniform array is characterized by comprising the following steps:
step 1: inputting the received signal into a self-adaptive covariance matrix calculation module, and calculating to obtain a self-adaptive covariance matrix of the received signal;
step 2: inputting the obtained adaptive covariance matrix into a multi-stage flow characteristic decomposition module, and performing multi-flow characteristic value decomposition to obtain a characteristic value and a characteristic vector of the covariance matrix;
and step 3: inputting the obtained characteristic values and the characteristic vectors into a dual-threshold signal source number estimation module, judging the number of signal sources and obtaining a noise subspace;
and 4, step 4: and inputting the output of the dual-threshold information source number estimation module into a spectral peak search module, and calculating to obtain a spectral function, and judging and comparing all extreme points to obtain a stored value, namely the direction of the information source.
2. The non-uniform array based coherent source DOA estimation FPGA realization method as recited in claim 1, further comprising: the adaptive covariance matrix calculation module in the step 1 comprises an autocorrelation matrix calculation module, a smoothing processing module and a real quantization processing module;
the autocorrelation matrix calculation module adopts multi-path parallel calculation to obtain the calculation result of each snapshot, and accumulates the calculation result of each snapshot to obtain the triangular 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 fast beats,the calculation formula of (A) 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, the accumulator accumulates the result of single snapshot and controls by using the input snapshot number L to obtain the calculation result of the autocorrelation matrix;
the smoothing processing module carries out summation operation on data of corresponding positions of the autocorrelation matrix according to the space smoothing principle, and carries out summation operation on the autocorrelation matrix obtained by the calculationThe elements at the corresponding positions are summed, and the calculation formula is:
wherein rpq(p, q ═ 1, 2.., 6) is an autocorrelation matrixElements of the corresponding position in;
and the real-number processing module is used for carrying out inverse operation on the smoothed partial result, converting the covariance matrix of the signal from a complex number domain to a real number domain, and obtaining Ru ═ λ u according to the characteristic decomposition property, wherein λ is the characteristic value of the covariance matrix, and R is the covariance matrix of the signal and is a complex matrix represented by R ═ Rr+iRi，RrRepresenting the real part of a matrix R, RiDenotes the imaginary part of the matrix R, u is the eigenvector of the covariance matrix, and u is expressed as u ═ ur+iui，urRepresenting the real part of the vector u, uiRepresenting the imaginary part of vector u, taken into the above equation:
will matrixAccording toAnd the conversion is to convert the complex matrix into a real matrix, solve the inverse 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.
3. The non-uniform array based coherent source DOA estimation FPGA realization method as recited in claim 2, further comprising: the multistage flow characteristic decomposition module is used for calculating the eigenvalue and the eigenvector of the 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 the data and the eigenvector matrix of the main diagonal line in the eigenvalue transition matrix; and each level of cleaning module finishes cleaning all upper triangular elements in the characteristic value transition matrix and updating the characteristic vector transition matrix.
4. The non-uniform array based coherent source DOA estimation FPGA realization method as recited in claim 1, further comprising: the dual-threshold signal source number estimation module is used for judging the number of signal sources and obtaining a noise subspace and comprises a characteristic value sequencing module, a threshold judgment module and a noise subspace solving module; the eigenvalue sorting module is used for sorting the input eigenvalues, the threshold value judging module is used for judging the number of the information sources, and the noise subspace solving module processes and screens the input eigenvector matrix according to the sorted and information source number estimation results to obtain a noise subspace matrix.
5. The non-uniform array based coherent source DOA estimation FPGA realization method as recited in claim 1, further comprising: the spectrum peak searching module in the step 4 comprises a guide vector calculating module, a space spectrum function calculating module, an extreme value checking module and a peak value judging module; the guiding vector calculation module inputs information source frequency, array element position coordinates and sin (theta) value stored in a ROM to calculate the value of the guiding vector, and the calculation formula is as follows:
where θ is the search angle, diIs the position coordinate of the array element, f is the frequency of the incident information source, and c is the speed of light; the input of the spatial spectrum function calculation module is a noise subspace and a guide vector, and the spectrum function calculation formula is as follows:
wherein Q is the number of sources, viIs the ith column vector in the noise subspace, and a (theta) is the guide vector; the extreme value checking module is used for judging minimum value points in the spectrum function and outputting the extreme value points to the peak value judging module, the peak value judging module judges the extreme value points and the stored values, if the minimum value points are smaller than the stored values, the stored values are updated, otherwise, the stored values are kept unchanged, and the stored values obtained after all the extreme value points are judged and compared are the direction of the information source.
6. The non-uniform array based coherent source DOA estimation FPGA implementation method of claim 3, wherein: each level 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 the updated eigenvalue transition matrix and eigenvector transition matrix are obtained through calculation; 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 finishing the first-stage cleaning.
7. The non-uniform array based coherent source DOA estimation FPGA realization method as recited in claim 6, further comprising: the single cleaning module comprises an input cache module, a rotation matrix calculation module, a characteristic value transition matrix calculation module, a characteristic vector transition matrix calculation module and an output cache 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 same calculation units and is used for calculating the sine and cosine values of the rotation angle, and the calculation formula of the rotation angle is as follows:
taking the first group S1 as an example, the rotation matrix is constructed by obtaining the rotation angles and is expressed as follows:
therefore, the eigenvalue matrix update calculation formula is:
R'＝WR'WT
the eigenvector moment update calculation formula is:
V'＝WV'
wherein the initial value of R' is a real symmetric matrixThe initial value of V' is the identity matrix.
8. The non-uniform array based coherent source DOA estimation FPGA implementation method as recited in claim 4, wherein: the characteristic value sequencing module comprises N identical comparison units, wherein N is the number of characteristic values, N paths are adopted to realize the comparison between every two characteristic values in parallel, one port of each comparison unit inputs a fixed characteristic value a, the other port inputs different characteristic values b in sequence, and if a is smaller than b, a result +1 is output; otherwise, if a is equal to b, comparing the position marks of the characteristic values a and b, and outputting a result +1 when the position of a is before the position of b; the position of a is not before the position of b, and the output result is kept 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 true CN110361691A (en) | 2019-10-22 |
| CN110361691B 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) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112543025A (en) * | 2020-12-08 | 2021-03-23 | 成都天奥信息科技有限公司 | High-speed serial AD sampling and data processing system and method based on matrixing |
| CN114237548A (en) * | 2021-11-22 | 2022-03-25 | 南京大学 | Method and system for complex dot product operation based on nonvolatile memory array |
| CN114325579A (en) * | 2022-03-09 | 2022-04-12 | 网络通信与安全紫金山实验室 | Positioning parameter estimation method, apparatus, device, storage medium and program product |
| CN114545324A (en) * | 2022-04-24 | 2022-05-27 | 南京宇安防务科技有限公司 | 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 (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080231505A1 (en) * | 2007-03-23 | 2008-09-25 | Weiqing Zhu | Method of Source Number Estimation and Its Application in Method of Direction of Arrival Estimation |
| 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 |
-
2019
- 2019-07-24 CN CN201910669180.5A patent/CN110361691B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080231505A1 (en) * | 2007-03-23 | 2008-09-25 | Weiqing Zhu | Method of Source Number Estimation and Its Application in Method of Direction of Arrival Estimation |
| 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算法的高速实现", 《合肥工业大学学报（自然科学版）》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112543025A (en) * | 2020-12-08 | 2021-03-23 | 成都天奥信息科技有限公司 | High-speed serial AD sampling and data processing system and method based on matrixing |
| CN112543025B (en) * | 2020-12-08 | 2023-03-14 | 成都天奥信息科技有限公司 | High-speed serial AD sampling and data processing system and method based on matrixing |
| CN114237548A (en) * | 2021-11-22 | 2022-03-25 | 南京大学 | Method and system for complex dot product operation based on nonvolatile memory array |
| CN114237548B (en) * | 2021-11-22 | 2023-07-18 | 南京大学 | Method and system for complex point multiplication operation based on nonvolatile memory array |
| CN114325579A (en) * | 2022-03-09 | 2022-04-12 | 网络通信与安全紫金山实验室 | Positioning parameter estimation method, apparatus, device, storage medium and program product |
| CN114545324A (en) * | 2022-04-24 | 2022-05-27 | 南京宇安防务科技有限公司 | Rapid direction finding method suitable for non-uniform array |
| 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 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110361691B (en) | 2023-05-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110361691A (en) | Coherent DOA based on nonuniform noise estimates FPGA implementation method | |
| CN104749553B (en) | Direction of arrival angle method of estimation based on rapid sparse Bayesian learning | |
| CN106019215B (en) | Nested array direction of arrival angle method of estimation based on fourth-order cumulant | |
| CN107589399B (en) | Estimation method of direction of arrival of co-prime array based on singular value decomposition of multi-sampling virtual signal | |
| Shi et al. | Accelerating parallel Jacobi method for matrix eigenvalue computation in DOA estimation algorithm | |
| CN110222307B (en) | Parallel implementation method for eigenvalue decomposition of real symmetric matrix based on FPGA | |
| CN106842113A (en) | Direction of arrival high accuracy method of estimation in the case of height 1 bit quantization of sampling | |
| CN106155627B (en) | Low overhead iteration trigonometric device based on T_CORDIC algorithm | |
| CN104793176A (en) | FPGA (field programmable gate array) based DOA (direction of arrival) estimation fast-implementing method | |
| CN110275131B (en) | DOA tracking method and device based on virtual differential array | |
| Jing et al. | An improved fast Root-MUSIC algorithm for DOA estimation | |
| CN109116295A (en) | The passive direction finding algorithm of baseline is chosen based on phased array | |
| Li et al. | Hardware acceleration of MUSIC algorithm for sparse arrays and uniform linear arrays | |
| CN109581277B (en) | A kind of four-dimensional antenna array DOA estimation method based on compressive sensing theory | |
| Chen et al. | Fast algorithm for DOA estimation with partial covariance matrix and without eigendecomposition | |
| 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 | |
| CN108614234B (en) | Direction-of-arrival estimation method based on multi-sampling snapshot co-prime array received signal fast Fourier inverse transformation | |
| Yan et al. | Computationally efficient direction of arrival estimation with unknown number of signals | |
| CN106507952B (en) | A kind of quick spatial spectrum computational methods based on circle battle array | |
| CN113203997B (en) | FPGA-based radar super-resolution direction finding method, system and application | |
| Pasupuleti et al. | Low complex & high accuracy computation approximations to enable on-device RNN applications | |
| CN109685400A (en) | Time-lag power system stability method of discrimination based on time integral IGD | |
| CN212569855U (en) | Hardware implementation device for activating function | |
| CN108008665B (en) | Large-scale circular array real-time beam former based on single-chip FPGA and beam forming calculation 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 |