1. Introduction
Direction of arrival (DOA) estimation is a significant problem in array signal processing and must be studied. DOA estimation is used extensively in radar, sonar, navigation, communication, and other fields [1,2]. Theoretically, in subspace class algorithms, such as multiple signal classification (MUSIC), DOA estimation exceeds the Rayleigh limit to achieve so-called super resolution [3,4]. However, the accuracy of the DOA estimation in engineering implementations is directly related to the accuracy of the covariance matrix eigenvalue decomposition [5,6]. Concurrently, the computational complexity and poor real-time performance of the algorithm limit its widespread application in practical engineering [7,8]. The coordinate rotation digital computer (CORDIC) algorithm is commonly used to construct orthogonal matrices for eigenvalue decomposition [9,10]. A previous study [11] used a double-precision floating-point algorithm instead of the CORDIC IP core to achieve high-precision DOA estimation. Because this method guarantees accuracy only when implemented serially, the parallel processing advantage of the algorithm is lost. An algorithm for fast eigenvalue decomposition has also been studied [12,13]. However, its implementation precision was not considered.
This study focuses on developing a method to obtain high accuracy in the calculation of the covariance matrix eigenvector in DOA estimation. Based on an analysis of the basic principles and characteristics of solving for the eigenvector of a complex matrix in DOA estimation, an efficient and highly accurate DOA estimation algorithm is proposed. The effectiveness of the proposed algorithm was verified using simulation results.
2. Principle of DOA Estimation
Multiple signal classification (MUSIC): The MUSIC algorithm realizes DOA estimation by searching for the peak of the array spatial spectral function [TeX:] $$P(\theta)$$ [14,15].
First, the covariance matrix was estimated using the snap data received by the array. For L snaps, the signal was considered a stationary Gaussian random process; its statistical characteristics remain the same over time, and the signal direction vector was linearly independent, thus making it a complex positive definite Hermite matrix.
Next, R was decomposed by eigenvalues to obtain the signal subspace [TeX:] $$U_{S}$$ and noise subspace [TeX:] $$U_{N}.$$ Subsequently, a peak search was conducted according to the following formula to determine the angle corresponding to the maximum value:
where [TeX:] $$U_{N}$$ is the noise subspace, [TeX:] $$a(\theta)$$ is the signal-steering vector, and H is the conjugate transpose.
According to Eq. (1), the estimation accuracy of the signal subspace [TeX:] $$U_{S}$$ and noise subspace [TeX:] $$U_{N}$$ directly affects the DOA estimation accuracy.
3. Analysis for Complex Matrix Eigenvalue Solving Algorithm
Compared with real symmetric matrices, the processing of complex Hermite matrices is more involved and requires complex covariance matrices to be transformed into real matrices. A typical Hermite matrix is represented as [12]:
where [TeX:] $$B^{T}=B, D^{T}=-D \in R^{n \times n}.$$
The eigenvalue decomposition of matrix A is performed, and a second-order primary subarray is constructed.
where K and L are subscripts of the matrix elements, i is the complex factor [TeX:] $$\sqrt{-1}, \text { and }-$$ indicates conjugation.
Subsequently, a unitary diagonal matrix is constructed.
The imaginary part of [TeX:] $$A_{K L} \text { is } 0, \text { and } \tilde{A}_{K L}$$ is the second-order real matrix:
Then, [TeX:] $$\tilde{A}_{K L}$$ can be diagonalized by Jacobi rotation [12] as
where [TeX:] $$U_{K L}^{H}=T_{K L} J_{K L}, J_{K L}=\left[\begin{array}{cc} C & -S \\ S & c \end{array}\right], s^{2}+c^{2}=1,$$
Further,
Thus, A is diagonalized and the eigenvectors are obtained through ergodic K and L.
Eq. (7) shows that numerous operations, such as squaring and division, are used in the process of constructing the orthogonal unit [TeX:] $$J_{K L},$$ which can be implemented only by serial calculations, thus resulting in the loss of the advantage of the availability of hardware parallel processing.
Let
From Eq. (7), we get
From Eqs. (9) and (10), given the value of [TeX:] $$\theta, \cos (\theta) \text { and } \sin (\theta),$$ that is, the values of c and s, can be obtained in parallel. Therefore, solving for [TeX:] $$\theta$$ is key to efficiently constructing a rotated orthogonal matrix [TeX:] $$J_{K L}$$ and improving the efficiency of the algorithm. The calculation of A, B, and C involves trigonometric function calculations, which are a key technology in algorithm implementation.
4. Parallel High-Precision Eigenvalue Decomposition Algorithm
The CORDIC algorithm is widely used to solve trigonometric functions. It is a unified algorithm that can calculate various transcendent functions. The underlying concept is to approach the target through coordinate rotation [13]. In engineering, to facilitate algorithm implementation, the tangent value of each rotation angle should be an integer multiple of 2. The principle of using fixed-point CORDIC to convert multiple complex algorithms into an iterative operation that requires only shift and subtraction operations was detailed in a previous study [12].
In practical applications, both the relative spatial positions between signals and the antenna front change rapidly; therefore, real-time estimation of DOA is an important issue. Maximizing parallelism in engineering implementations is an important method for ensuring the real-time performance of algorithms. Therefore, based on the analysis presented in Section 2, a parallel algorithm for the complex Hermite matrix eigenvalue decomposition is presented herein. The process is illustrated in Fig. 1.
Based on the CORDIC algorithm analysis, the decomposition accuracy of the matrix eigenvalues depends on the accuracy of the rotation matrix [TeX:] $$J_{K L},$$ that is, the calculation accuracy of the associated trigonometric functions. The CORDIC algorithm assumes that the tangent value of each rotation angle is an integer multiple of 2. Therefore, complex calculations are replaced by multiples of [TeX:] $$2^{-i} \text { and } \pm 1.$$ Thus, to ensure the accuracy of the eigenvector solution, a single-precision floating point was used to complete the trigonometric function solution.
Parallel algorithm for Hermite matrix eigenvalue decomposition.
Single-precision floating-point storage structure.
The single-precision floating-point storage structure according to the IEEE745 standard is shown in Fig. 2.
Single-precision floating-point numbers consist of sign, order, and tail bits, which represent the following values.
Diagram of DOA estimation algorithm.
Based on the storage structure of single-precision floating-point numbers, the multiplication of a single-precision floating-point number and [TeX:] $$2^{-i}$$ is equivalent to the order e minus i. The multiplication of a single-precision floating-point number and [TeX:] $$\pm 1$$ can only change the sign bit. Thus, the calculations based on fixed-point numbers, described above, were performed by shifting and negating, which includes bitwise inversion and plus 1. Therefore, the complexity of single-precision floating-point calculations did not increase compared with fixed-point calculations.
The DOA estimation algorithm is shown in Fig. 3. The basic principles are described in Section 2. The parallel algorithm for Hermite matrix eigenvalue decomposition in Fig. 3 is shown in Fig. 1. The upper limit of iteration time can be set as required. The upper limit of the iterations in this study was chosen to be 10.
5. Algorithm Simulation Analysis
The input conditions for the simulations were as follows. A 4-element array antenna and single-frequency signal were used. The power ratio of the signal-to-noise, elevation of the signal, and azimuth were set at values of 50 dB, [TeX:] $$50 \mathrm{~dB}, 50^{\circ}, \text { and } 10^{\circ},$$ respectively. The performances of the double-precision floating-point calculation, fixed-point calculation, and the proposed algorithm were compared via simulations. The number of iterations for CORDIC of fixed-point and this study was set at 16.
The results of 100 Monte Carlo simulations indicate that the average error in the eigenvector calculation using the proposed algorithm is only 0.000001 compared to that in the double-precision floating-point calculation [11], whereas that of the fixed-point algorithm [10] is 0.005. Some results from eigenvector calculations using the three methods are listed in Table 1 for comparison.
Comparison of eigenvector results
The accuracy of the algorithm in this study was equivalent to that of the double-precision floating-point algorithm [11]. However, the double-precision floating-point algorithm can only complete eigenvalue diagonalization in serial mode and requires the square and division of double-precision floating point numbers. The FPGA implementation requires 300 clock-cycles, whereas the algorithm in this study only requires 70 clock-cycles to complete, 230 clock-cycles fewer than those required for similar calculations presented in the literature [11]. The DOA algorithm must traverse all six non-diagonal elements of the complex matrix, and the clock-cycles used in this algorithm are 1,380 fewer than those in [11]. The entire eigenvalue decomposition process was iterated 16 times, and the clock-cycles used in this algorithm were 22,080 less than those in [11]. The amount of computation required by this algorithm is equivalent to that required by the fixed-point algorithm [10].
The spatial spectrum 3D-plot based on fixed DOA estimation is shown in Fig. 4, and the spatial spectral isobars based on fixed DOA estimation are shown in Fig. 5. The spatial spectrum 3D-plot based on the proposed algorithm is shown in Fig. 6, and the spatial spectral isobars based on the proposed algorithm are shown in Fig. 7.
Spatial spectrum 3D based on fixed DOA estimation.
Spatial spectral isobars based on fixed DOA estimation.
From the simulation results, the estimated elevation and azimuth angles obtained using the proposed algorithm were equal to the set values, and the decomposition accuracy based on the algorithm satisfied the DOA estimation application requirements. The divergence of the DOA estimation results occurred in fixed-point calculations owing to the large implementation errors.
Spatial spectrum 3D based on the proposed algorithm.
Spatial spectral isobars based on the proposed algorithm.
The double-precision floating-point algorithm proposed in [11] can only be implemented serially and uses calculations of squares of double-precision floating-point numbers, which requires at least 300 clock-cycles to diagonalize the matrix using the field-programmable gate array. The proposed algorithm can be completed in 70 clock cycles, and its computational efficiency is 76.67 % higher than that of the double-precision floating-point algorithm proposed in [11], which significantly improves the real-time performance of DOA estimation, while ensuring accuracy.
6. Conclusion
Although DOA estimation of the MUSIC algorithm can theoretically exceed the Rayleigh limit and achieve the supposed super-resolution, its application is limited in practical engineering owing to the real-time performance and accuracy of the algorithm. Ensuring the efficiency and accuracy of an algorithm in engineering is a key problem. Based on the analysis of the characteristics of the MUSIC algorithm, complex matrix eigenvalue decomposition, and CORDIC algorithm, a high-accuracy algorithm for DOA estimation with high parallelism was presented herein. The algorithm uses parallel construction of an orthogonal unit matrix and the CORDIC algorithm based on a single-precision floating-point number to ensure the efficiency and accuracy of the algorithm, respectively. Monte Carlo simulation verified that the average error of the feature vector of the algorithm was only 0.000001 compared with the double-precision floating-point calculation, whereas that of the fixed-point algorithm was 0.005. The divergence in the DOA estimation results appeared because of large error in the fixed-point calculation. The diagonalization of double-precision floating-point computation matrices requires at least 300 clock cycles, whereas the proposed algorithm performed a similar calculation in only 70 clock cycles, which signifies a notable improvement in the real-time performance and accuracy of the algorithm.