A Fast Estimation Method for Direction of Arrival Using Tripole Vector Antenna

Abstract

The tripole vector antenna comprises three orthogonal dipole antennas, so it could completely capture all the electric field of the incident electromagnetic (EM) wave. Then, the electric field information could be used to estimate the direction of arrival (DOA) of the EM wave if two conditions are satisfied. One is that there exists only one single EM wave in space. The other is that the EM wave is elliptically or circularly polarized. The new estimation method obtains two snapshot vectors through the output of a tripole antenna and computes their cross-product vector. Furthermore, the direction of the cross-product vector is used to estimate the DOA of the EM wave directly. We analyze the statistical characteristics of the DOA estimation error to prove that the new scheme is an asymptotic unbiased estimation. Furthermore, unlike the existing Multiple Signal Classification (MUSIC)-based algorithms, the proposed approach only need one tripole vector antenna instead of an antenna array. Meanwhile, the new method also outperforms existing MUSIC-based algorithms in the term of computational complexity. Finally, the performance and advantages of the proposed method are verified by numerical simulations.

1. Introduction

The incident angle of the electromagnetic (EM) wave relative to the receiver is called direction of arrival (DOA). In radar, communication, radio astronomy, telemetry, track and command, and other fields [1,2,3], users always need to estimate the DOA of EM waves. For example, the DOA information wave can be used in the beam control system of the receiver antenna, so that the main beam of the receiver antenna can track the incident signal all the time. Different DOA estimation methods are adopted in various application fields. Each method has its own advantages and keeps developing continuously.

Nowadays, the most traditional DOA estimation equipment is the interferometer, which estimates the DOA according to the phase difference between EM waves received by adjacent antenna elements. The phase difference is caused by the difference in distance between the radiation source and different antenna elements, and it has a certain geometric relationship with the incident angle. However, there are some shortcomings in DOA estimation of the interferometer. Firstly, the estimation accuracy is not high because the algorithm is vulnerable to noise. Obtaining the phase difference is a key step for interferometer. Common methods include zero-crossing detection, Fast Fourier Transform (FFT) [4], and Hilbert orthogonal transformation [5] and so on. The zero-crossing detection uses the approximation of the function near the zero point, i.e., sinx = x, to accurately determine the zero-crossing position of the sample sequence. These positions can help to calculate the path difference between the signals received by different antenna elements. Unfortunately, noise results in serious uncertainty when estimating the zero point. In addition to noise, the estimation accuracy is also related to the sample frequency of the signal, which requires high performance of the hardware. FFT is another way to estimate the phase difference between signals. The phase difference between the two signals can be easily obtained by FFT and cross-correlation operation on the sample sequence. Compared with zero-crossing detection, the FFT-based algorithm performs much better in anti-noise and anti-harmonic interference, so it could achieve higher accuracy in the phase difference estimation. The second shortcoming for the interferometer is that the estimation accuracy is closely related to the array spacing. The shortest interval in an interferometer array is usually less than half a wavelength; hence, the phase measured in this case is unique [6,7]. However, in broadband direction-finding systems over a range such as 0.1–3 GHz or 2–8 GHz, the physical dimensions of the antenna element are dictated by the lowest operating frequency, whereas the shortest interelement spacing is decided by the highest frequency. This poses the problem of not being able to employ a short spacing of less than half-wavelength [8,9]. However, if the baseline length exceeds the half-wavelength of the pulse signal, the periodic ambiguity will be introduced into the phase measurement and its ambiguity is 2π. The combination of long baseline and short baseline is an effective method to solve the phase ambiguity. Even so, the combination of multiple baselines greatly increases the complexity of the system and requires more accurate calibration of the system. Thirdly, the interferometer can only aim at a single signal in the main lobe, and cannot be used when multiple EM waves arrive at the antenna at the same time. In this case, it means insufficient resolution.

In order to promote the resolution and realize the DOA estimation of multiple incident signals, various methods based on antenna arrays have become the research hotspot in the last decades. So far, several classes of DOA estimation techniques for incoherently distributed sources have been developed, such as the subspace method and the parameter method. Multiple Signal Classification (MUSIC) and Estimating Signal Parameter via Rotational Invariance Techniques (ESPRIT) are the most important subspace methods [10,11,12]. Furthermore, the maximum likelihood (ML) algorithm [13,14] and Covariance Matching Estimator (COMET) algorithm [15,16] are typical parameter estimation approaches, which can obtain better estimation performance than subspace methods. However, their computing load is also heavier. Meanwhile, DOA estimation methods based on antenna arrays are facing two great challenges. First, the complexity of the algorithm is too high to implement from an engineering perspective. In these technologies, most methods are proposed for uniform linear arrays; therefore, the spectral peak search happens only on a two-dimensional plane, which makes it impossible to estimate the elevation angle and azimuth angle simultaneously. Later, researchers further developed these algorithms to make them suitable for planar arrays. In this way, the range of the spectral peak search is extended to three-dimensional space, which means that the DOA of the signal source can completely include both the elevation angle and azimuth angle. However, due to multi-dimensional space spectrum search, a large amount of computational burden has always been a problem that researchers need to solve [17]. In addition to the huge computation of the spectral peak search, this kind of DOA algorithm involves various complex operations such as matrix eigenvalue decomposition, which resulting in difficulties in the engineering implementation of programmable devices, such as FPGA. Another challenge is that these antenna arrays generally have a large number of elements and each element needs a corresponding radio frequency front-end channel, which makes the hardware cost very high.

Different from the popular algorithms mentioned above, we use a tripole vector antenna to make full use of the electric field of the EM wave and develop a fast DOA estimation algorithm for a single signal. The tripole vector antenna is also known as a polarization sensitive antenna, which is different from the traditional single polarization antenna because it can accurately sense the polarization characteristics of EM waves [18]. The tripole vector antenna is proposed by R.T Compton [19,20], and it is composed of three orthogonal dipole antennas. Document [21] proposed a DOA estimation method based on polarization-spatial steering vector. It uses a tripole antenna but needs to know the elliptical susceptibility of the polarized signal. In recent years, polarization-sensitive arrays have become the focus of academic research, where each element is one tripole vector antenna. Researchers try to achieve higher performance through polarization-sensitive arrays in the fields of spatial spectrum estimation and adaptive anti-jamming array. In fact, it is a signal processing method combining spatial and polarimetry resources, which improves the array freedom degree of the receiver. For example, the tripole antenna has been proposed as one possible choice of antenna construction to be used in the low-frequency array (LOFAR) project [22]. The primary field of application for the project is radio-based astrophysics.

Due to the structural characteristics of the tripole vector antenna, we only need to use one tripole vector antenna to complete the DOA estimation of the EM wave if there is only one signal in space. We map the sample data of the tripole vector antenna to the space rectangular coordinate system and calculate the DOA by vector product according to the geometric relationship. The new method has some obvious advantages over the other methods. First of all, the proposed DOA algorithm does not need to solve equations such as the interferometer or search the spectral peak such as array signal processing methods [23]; as a result, the signal processing process is much simpler. Secondly, the main lobe of the tripole vector antenna is wide enough so that any signal impinging from the main lobe area could be captured and no servo system is needed. Obviously, this reduces the cost of hardware greatly. Thirdly, the DOA estimation method based on the tripole vector antenna is suitable for any nonlinearly polarized EM wave without any polarization-matching problem between antennas and EM waves, because it captures all the polarization information of EM waves. This is also the advantage that other DOA estimation methods do not have. Similar to the interferometer, the new method also has the disadvantage of insufficient resolution, but it can be solved by combining with other subspace estimation algorithms.

The rest of this paper is as follows. Section 2 is a brief introduction of the received signal model. Section 3 is the overall description of the method. The method complexity and error analysis are discussed in Section 4. The experimental results of different methods are presented in Section 5. Section 6 concludes the paper.

2. Signal Model

The tripole vector antenna is composed of three identical and mutually perpendicular linearly polarized dipole antennas centered at the same position [19,20,21], as shown in Figure 1. Each linearly polarized dipole antenna can only receive the electric field component in its polarization direction and has an output feed point. The polarization direction of the dipoles in a tripole are respectively parallel to the three coordinate axes of the space rectangular coordinate system, which ensures that no matter which direction the EM wave is incident from, and no matter what polarization the EM wave is, all the energy of the EM wave can be received by the tripole vector antenna. So, there will be no polarization-matching problem. In addition, since all three feed points are concentrated at the same location, phase shifts between units will not occur in a tripole vector antenna.

In the coordinate system shown in Figure 1, we suppose that the DOA of an elliptically polarized EM wave signal is incident from a certain direction. Once the signal reaches the antenna, its electric field will be simultaneously received by the three orthogonal dipoles of the antenna. We name the three poles of the antenna as x, y, z, hence the signal received by the antenna can be expressed as

$$\begin{array}{ll}x\left(n\right)& =E_x\left(n\right)+n_x\left(n\right)\\ & =E_{xm}\cos \left(\omega n{T}_s+{\phi }_x\right)+n_x(n)\end{array}$$

(1)

$$\begin{array}{ll}y\left(n\right)& =E_y\left(n\right)+n_y(n)\\ & =E_{ym}\cos \left(\omega n{T}_s+{\phi }_y\right)+n_y(n)\end{array}$$

(2)

$$\begin{array}{ll}z\left(n\right)& =E_z\left(n\right)+n_z(n)\\ & =E_{zm}\cos \left(\omega n{T}_s+{\phi }_z\right)+n_z(n)\end{array}$$

(3)

where E_xm, E_ym, and E_zm and ϕ_x, ϕ_y, and ϕ_z are the amplitudes and initial phases of the electric field components in the x, y, and z directions. n_x(n), n_y(n), and n_z(n) are independent white Gaussian noise in the x, y, and z directions. n is the sample number, T_s is the sample period, and ω is the angular frequency of the EM wave.

3. DOA Estimation Algorithm

3.1. Principle of Algorithm

For an elliptically polarized EM wave, according to the signal model, in the antenna coordinate system, the electric field of the received signal can be expressed as

$$E\left(n\right)=x\left(n\right)e_x+y\left(n\right)e_y+z\left(n\right)e_z$$

(4)

where $E\left(n\right)$ is the electric field of the received signal, and e_x, e_y, and e_z are unit direction vectors in the x, y, and z directions respectively.

According to the definition of elliptically polarized electromagnetic waves, the electric field vector E(n) is orthogonal to the propagation direction, and its trajectory is elliptical when viewed from the direction of propagation. Obviously, the output vector of the tripole vector antenna is equal to the electric field vector E(n) at the origin of the coordinate. In addition, the plane of the trajectory must be perpendicular to the incident direction of the EM wave. According to the definition of the vector operation, the cross-product of two non-parallel vectors is perpendicular to the plane where they are located. Therefore, we can calculate the direction vector of the incident EM wave through these two snapshots. Of course, this also applies to circular polarization, because circular polarization is only a special case of elliptical polarization.

3.2. Caculation of Direction Vector

In the space rectangular coordinate system as shown in Figure 2, without considering noise, the electric field trajectory is an ellipse in the plane $\alpha$. Two continuous sample snapshots of the tripole vector antenna are taken to describe two electric field vectors p₁ and p₂ at two different moments, which are

$$p_1={[\begin{array}{ccc}{E}_x\left(2k\right)& E_y\left(2k\right)& E_z\left(2k\right)\end{array}{\displaystyle ]}}^T$$

(5)

$$p_2={\left[\begin{array}{ccc}{E}_x\left(2k+1\right)& E_y\left(2k+1\right)& E_z\left(2k+1\right)\end{array}\right]}^T$$

(6)

where symbol “T” means transposition. 2k and 2k + 1 are two continuous sample indexes.

Apparently, as shown in the Figure 2, the angle between p₁ and p₂ is certainly less than π and greater than 0 because of the Nyquist sample theorem. Therefore, these two vectors are not on the same line. According to vector operation principle, the normal vector of the plane α can be represented by the cross product of p₁ and p₂, which is

$$m\left(k\right)=p_1\times p_2={\left[\begin{array}{ccc}{m}_x& m_y& m_z\end{array}\right]}^T$$

(7)

where

$$m_x=E_y\left(2k\right)E_z\left(2k+1\right)-E_z\left(2k\right)E_y\left(2k+1\right)$$

(8)

$$m_y=E_z\left(2k\right)E_x\left(2k+1\right)-E_x\left(2k\right)E_z\left(2k+1\right)$$

(9)

$$m_z=E_x\left(2k\right)E_y\left(2k+1\right)-E_y\left(2k\right)E_x\left(2k+1\right).$$

(10)

Obviously, the vector m is orthogonal to the plane $\alpha$, since it is the cross-product of p₁ and p₂. Then, the direction of the vector m is parallel to the DOA of the EM wave, and that means we can estimate the DOA by computing the angle of vector m. From the engineering point of view, we only consider incident EM waves from the upper hemisphere. Therefore, for incident EM waves from the lower hemisphere, all coordinates of vector m are inverted before the DOA estimation.

3.3. Noise-Free Angle Estimation Algorithm

For convenience, we assume temporarily that the received signal contains no noise. According to Figure 2, the elevation angle θ₀ (θ₀ ∈ [0, π/2]) is the angle between the vector m and the z axis. The azimuth angle ϕ₀ is the angle between the projection of m in the xy plane and the x axis. If this projection is in quadrant 1 or 2, then ϕ₀ ∈ [0, π], and if it is in quadrant 3 or 4, then ϕ₀ ∈ [−π, 0]. From Figure 2, the elevation angle ${\theta }_0$ should satisfy.

$$\cos {\theta }_0=\frac{m_z^{}}{\sqrt{m_x^2+m_y^2+m_z^2}}$$

(11)

The value of Equation (11) is always not negative because we only consider EM waves coming above the xy plane; thus, the estimated elevation angle θ₀ can be calculated by the inverse cosine function, which is

$${\widehat{\theta }}_0=\arccos \left(\frac{m_z}{\sqrt{m_x^2+m_y^2+m_z^2}}\right).$$

(12)

On the one hand, when ${\widehat{\theta }}_0=0$, it is obvious that the vector m is parallel to the z axis, which means it is meaningless to estimate azimuth, so we stop estimating. On the other hand, when ${\widehat{\theta }}_0\ne 0$, the vector m is not parallel to the z axis, hence

$$m_x^2+m_y^2\ne 0.$$

(13)

Similarly, according to the relationship between the projection of vector m onto the xy plane and the x axis in Figure 2, it can be known that

$$\cos {\varphi }_0=\frac{m_x}{\sqrt{m_x^2+m_y^2}}.$$

(14)

Furthermore, according to the definition of azimuth, the estimated value of azimuth is

$${\widehat{\varphi }}_0=\{\begin{array}{c}\arccos \left(\frac{m_x}{\sqrt{m_x^2+m_y^2}}\right)\begin{array}{ccc}& & m_y\ge 0\end{array}\\ -\arccos \left(\frac{m_x}{\sqrt{m_x^2+m_y^2}}\right)\begin{array}{ccc}& & m_y\le 0\end{array}\end{array}.$$

(15)

3.4. Actual Angle Estimation Algorithm

The implementation of the algorithm is shown in Figure 3. In order to decrease the effect of noise, we can take the average value with multiple estimates. Obviously, the calculation of the anti-trigonometric function in Equations (12) and (15) are the most complex calculation steps in the whole estimation algorithm, and these steps are also the most sensitive to noise. Therefore, when calculating the vector m, taking the average value with K sample snapshots is an effective action to reduce the influence of noise, because this operation is before the calculation of anti-trigonometric functions. Then, the estimated value of m is

$$\overline{m}=\frac{1}{K}{\displaystyle \sum _{k=1}^{K}m\left(k\right)}={\left[\begin{array}{ccc}{\overline{m}}_x& {\overline{m}}_y& {\overline{m}}_z\end{array}\right]}^T.$$

(16)

Correspondingly, DOA estimates are modified to

$${\widehat{\theta }}_0=\arccos \left(\frac{{\overline{m}}_z}{\sqrt{{\overline{m}}_x^2+{\overline{m}}_y^2+{\overline{m}}_z^2}}\right)$$

(17)

$${\widehat{\varphi }}_0=\{\begin{array}{c}\arccos \left(\frac{{\overline{m}}_x}{\sqrt{{\overline{m}}_x^2+{\overline{m}}_y^2}}\right)\begin{array}{ccc}& & {\overline{m}}_y\ge 0\end{array}\\ -\arccos \left(\frac{{\overline{m}}_x}{\sqrt{{\overline{m}}_x^2+{\overline{m}}_y^2}}\right)\begin{array}{ccc}& & {\overline{m}}_y\le 0\end{array}\end{array}.$$

(18)

4. Performance Analysis

4.1. Probability Distribution of Estimated Angles

Since there is noise in the snapshot, the estimated angles ${\widehat{\theta }}_0$ and ${\widehat{\varphi }}_0$ are considered as random variables. In order to find out the Probability Density Function (PDF), the PDFs of m_x, m_y, and m_z need to be calculated first.

Considering noise components in received signals, we plug Equations (1)–(3) into Equations (8)–(10) respectively, and obtain m_x, m_y, and m_z with noise. Since the three variables have the same structure, only the PDF of m_x is analyzed. For m_x, which can be expanded as

$$m_x=m_{x0}+P_x+Q_x$$

(19)

where

$$m_{x0}=y\left(2k\right)z\left(2k+1\right)-z\left(2k\right)y\left(2k+1\right)$$

(20)

$$P_x=n_y\left(2k\right)n_z\left(2k+1\right)-n_z\left(2k\right)n_y\left(2k+1\right)$$

(21)

$$\begin{array}{ll}{Q}_x& =n_y\left(2k\right)E_z^{}\left(2k+1\right)+n_z\left(2k+1\right)E_y^{}\left(2k\right)\\ & -n_z\left(2k\right)E_y^{}\left(2k+1\right)-n_y\left(2k+1\right)E_z^{}\left(2k\right)\end{array}.$$

(22)

Obviously, m_x0 is ideal value without noise, but P_x and Q_x have noise. Apparently, noise snapshot n_x, n_y, and n_z are Gaussian random variables with a mean 0 and variance ${\sigma }_x^2$, ${\sigma }_y^2$, and ${\sigma }_z^2$ respectively.

Since Q_x is the sum of multiple independent Gaussian random variables, Q_x is also a random variable of Gaussian distribution. According to Equation (22), it is not difficult to find out the mean and variance are respectively:

$$E\left[Q_x\right]=0$$

(23)

$$\begin{array}{ll}D\left[Q_x\right]=& {\sigma }_{Q_x}^2={\sigma }_y^2{E}_z^2\left(2k+1\right)+{\sigma }_z^2{E}_y^2\left(2k\right)\\ & \begin{array}{ccc}& & +{\sigma }_z^2{E}_y^2\left(2k+1\right)+{\sigma }_y^2{E}_z^2\left(2k\right)\end{array}\end{array}.$$

(24)

Therefore, the PDF of Q_x is:

$$f_{Q_x}\left(Q_x\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{Q_x}^{}}\exp \left(-\frac{Q_x^2}{2{\sigma }_{Q_x}^2}\right).$$

(25)

In order to obtain the PDF of P_x, two related theorems are given in advance.

Theorem 1. 

If two independent random variables x and y obey the Gaussian distribution with mean 0 and variance${\sigma }_x^2$,${\sigma }_y^2$the PDF of the variable z = xy is:

$$f_z\left(z\right)=g\left(z\right)=\frac{K_0\left(\frac{\left|z\right|}{{\sigma }_x^{}{\sigma }_y^{}}\right)}{\pi {\sigma }_x^{}{\sigma }_y^{}}$$

(26)

where the function K₀ is the modified Bessel function of the second kind of order 0.

Theorem 2 

[22].If the PDFs of two independent random variables x and y are f_x(x) and f_y(y) respectively, the PDF of the variable z = x + y is:

$$f_z\left(z\right)={\displaystyle {\int }_{-\infty }^{+\infty }{f}_x\left(t\right)f_y\left(z-t\right)dt}.$$

(27)

Looking back to Equation (21), function P_x is the difference of two independent components, and each independent component is the product of two independent Gaussian variables. Hence, using Theorems 1 and 2, we can get the PDF of the variable P_x:

$$f_{P_x}\left(P_x\right)={\displaystyle {\int }_{-\infty }^{+\infty }g\left(t\right)g\left(t-P_x\right)dt}=\frac{1}{2{\sigma }_z{\sigma }_y}\exp \left(\frac{-\left|P_x\right|}{{\sigma }_z{\sigma }_y}\right).$$

(28)

According to Equation (19), and using Theorem 2 again, we can get the PDF of variable m_x; similarly, we can get PDFs of m_y, m_z.

$$\begin{array}{l}{f}_{m_x}\left(m_x\right)={\displaystyle {\int }_{-\infty }^{+\infty }{f}_{P_x}\left(t\right)f_{Q_x}\left(m_x-m_{x0}-t\right)dt}\\ ={\displaystyle {\int }_{-\infty }^{+\infty }\frac{1}{2\sqrt{2\pi }{\sigma }_{Q_x}{\sigma }_z{\sigma }_y}\exp \left(-\frac{{\left(m_x-t-m_{x0}\right)}^2}{2{\sigma }_{Q_x}^2}-\frac{\left|t\right|}{{\sigma }_z{\sigma }_y}\right)dt}\end{array}$$

(29)

$$\begin{array}{l}{f}_{m_y}\left(m_y\right)={\displaystyle {\int }_{-\infty }^{+\infty }{f}_{P_y}\left(t\right)f_{Q_y}\left(m_y-m_{y0}-t\right)dt}\\ ={\displaystyle {\int }_{-\infty }^{+\infty }\frac{1}{2\sqrt{2\pi }{\sigma }_{Q_y}{\sigma }_z{\sigma }_x}\exp \left(-\frac{{\left(m_y-t-m_{y0}\right)}^2}{2{\sigma }_{Q_y}^2}-\frac{\left|t\right|}{{\sigma }_z{\sigma }_x}\right)dt}\end{array}$$

(30)

$$\begin{array}{l}{f}_{m_z}\left(m_z\right)={\displaystyle {\int }_{-\infty }^{+\infty }{f}_{P_z}\left(t\right)f_{Q_z}\left(m_z-m_{z0}-t\right)dt}\\ ={\displaystyle {\int }_{-\infty }^{+\infty }\frac{1}{2\sqrt{2\pi }{\sigma }_{Q_z}{\sigma }_x{\sigma }_y}\exp \left(-\frac{{\left(m_z-t-m_{z0}\right)}^2}{2{\sigma }_{Q_z}^2}-\frac{\left|t\right|}{{\sigma }_x{\sigma }_y}\right)dt}\end{array}$$

(31)

According to Equation (24), similarly, we can define ${\sigma }_{Q_y}^2$, ${\sigma }_{Q_z}^2$ as follows

$${\sigma }_{Q_y}^2={\sigma }_x^2{E}_z^2\left(2k+1\right)+{\sigma }_z^2{E}_x^2\left(2k\right)+{\sigma }_z^2{E}_x^2\left(2k+1\right)+{\sigma }_x^2{E}_z^2\left(2k\right)$$

(32)

$${\sigma }_{Q_z}^2={\sigma }_x^2{E}_y^2\left(2k+1\right)+{\sigma }_y^2{E}_x^2\left(2k\right)+{\sigma }_y^2{E}_x^2\left(2k+1\right)+{\sigma }_x^2{E}_y^2\left(2k\right).$$

(33)

In order to simplify Equations (29)–(31), we replace exp(−|t|) with a function, which fits well with the original function. This alternative function is

$$h\left(t\right)=a_1{e}^{\frac{-t^2}{b_1}}+a_2{e}^{\frac{-t^2}{b_2}}$$

(34)

where a₁ = 0.4652, b₁ = 0.5524, a₂ = 0.4396, b₂ = 1.911, and max |exp(−|t|) − h(t)| is less than 0.03, which means these two functions fit well.

In this way, Equations (29)–(31) can be replaced as:

$$\begin{array}{ll}{f}_{m_x}\left(m_x\right)& =\frac{a_1}{4{\sigma }_{Q_x}^2{\sigma }_z^2{\sigma }_y^2}\sqrt{\frac{{\sigma }_z^2{\sigma }_y^2+2b_1{\sigma }_{Q_x}^2}{b_1}}\exp \left(-\frac{2b_1{\sigma }_{Q_x}^2{\sigma }_z^2{\sigma }_y^2}{{\sigma }_z^2{\sigma }_y^2+2b_1{\sigma }_{Q_x}^2}{\left(m_x-m_{x0}\right)}^2\right)\\ & +\frac{a_2}{4{\sigma }_{Q_x}^2{\sigma }_z^2{\sigma }_y^2}\sqrt{\frac{{\sigma }_z^2{\sigma }_y^2+2b_2{\sigma }_{Q_x}^2}{b_2}}\exp \left(-\frac{2b_2{\sigma }_{Q_x}^2{\sigma }_z^2{\sigma }_y^2}{{\sigma }_z^2{\sigma }_y^2+2b_2{\sigma }_{Q_x}^2}{\left(m_x-m_{x0}\right)}^2\right)\end{array}$$

(35)

$$\begin{array}{ll}{f}_{m_y}\left(m_y\right)& =\frac{a_1}{4{\sigma }_{Q_y}^2{\sigma }_z^2{\sigma }_x^2}\sqrt{\frac{{\sigma }_z^2{\sigma }_x^2+2b_1{\sigma }_{Q_y}^2}{b_1}}\exp \left(-\frac{2b_1{\sigma }_{Q_y}^2{\sigma }_z^2{\sigma }_x^2}{{\sigma }_z^2{\sigma }_x^2+2b_1{\sigma }_{Q_y}^2}{\left(m_y-m_{y0}\right)}^2\right)\\ & +\frac{a_2}{4{\sigma }_{Q_y}^2{\sigma }_z^2{\sigma }_x^2}\sqrt{\frac{{\sigma }_z^2{\sigma }_x^2+2b_2{\sigma }_{Q_y}^2}{b_2}}\exp \left(-\frac{2b_2{\sigma }_{Q_y}^2{\sigma }_z^2{\sigma }_x^2}{{\sigma }_z^2{\sigma }_x^2+2b_2{\sigma }_{Q_y}^2}{\left(m_y-m_{y0}\right)}^2\right)\end{array}$$

(36)

$$\begin{array}{ll}{f}_{m_z}\left(m_z\right)& =\frac{a_1}{4{\sigma }_{Q_z}^2{\sigma }_x^2{\sigma }_y^2}\sqrt{\frac{{\sigma }_x^2{\sigma }_y^2+2b_1{\sigma }_{Q_z}^2}{b_1}}\exp \left(-\frac{2b_1{\sigma }_{Q_z}^2{\sigma }_x^2{\sigma }_y^2}{{\sigma }_x^2{\sigma }_y^2+2b_1{\sigma }_{Q_z}^2}{\left(m_z-m_{z0}\right)}^2\right)\\ & +\frac{a_2}{4{\sigma }_{Q_z}^2{\sigma }_x^2{\sigma }_y^2}\sqrt{\frac{{\sigma }_x^2{\sigma }_y^2+2b_2{\sigma }_{Q_z}^2}{b_2}}\exp \left(-\frac{2b_2{\sigma }_{Q_z}^2{\sigma }_x^2{\sigma }_y^2}{{\sigma }_x^2{\sigma }_y^2+2b_2{\sigma }_{Q_z}^2}{\left(m_z-m_{z0}\right)}^2\right)\end{array}.$$

(37)

In order to get PDFs of ${\widehat{\theta }}_0$ and ${\widehat{\varphi }}_0$, we quote two more related theorems.

Theorem 3 

[22].If$u=g\left(x,y,z\right)$, and$z\left(u,x,y\right)$is the inverse function of u on z, u has continuous partial derivative to x, y, z, the joint PDF of x, y, z is$f\left(x,y,z\right)$, hence the PDF of u is

$$g\left(u\right)={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }f\left(z\left(u,x,y\right),x,y\right)\left|\frac{\partial z}{\partial u}\right|dx}dy}.$$

(38)

Theorem 4 

[22].If$u=g\left(x,y\right)$, and$y\left(u,x\right)$is the inverse function of u on y, u has continuous partial derivative to x, y, the joint PDF of x, y is$f\left(x,y\right)$, hence the PDF of u is

$$g\left(u\right)={\displaystyle {\int }_{-\infty }^{+\infty }f\left(y\left(u,x\right),x\right)\left|\frac{\partial y}{\partial u}\right|dx}.$$

(39)

If we define the joint PDF of random variables (m_x, m_y, m_z) is f₁(m_x, m_y, m_z), then

$$f_1\left(x,y,z\right)=f_{m_x}\left(x\right)f_{m_y}\left(y\right)f_{m_z}\left(z\right).$$

(40)

Replacing m_x, m_y, and m_z with x, y, and z, according to Equation (12) and Theorem 3, the PDF of the estimated elevation angle ${\widehat{\theta }}_0$ can be obtained as

$$f_{\theta }\left({\widehat{\theta }}_0\right)={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }{f}_1\left(x,y,z\left(x,y,{\widehat{\theta }}_0\right)\right)\left|\frac{\partial \left(x,y,z\right)}{\partial \left({\widehat{\theta }}_0,x,y\right)}\right|dxdy}}$$

(41)

where

$$\left|\frac{\partial z}{\partial {\widehat{\theta }}_0}\right|=\frac{\sqrt{x^2+y^2}}{{\sin }^2{\widehat{\theta }}_0}.$$

(42)

So,

$$f_{\theta }\left({\widehat{\theta }}_0\right)={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }{f}_{m_x}\left(x\right)f_{m_y}\left(y\right)f_{m_z}\left(\sqrt{x_{}^2+y_{}^2}\cot {\widehat{\theta }}_0\right)\frac{\sqrt{x^2+y^2}}{{\sin }^2{\widehat{\theta }}_0}dxdy}}.$$

(43)

In the same way, if the joint PDF of random variables (m_x, m_y) is

$$f_2\left(x,y\right)=f_{m_x}\left(x\right)f_{m_y}\left(y\right),$$

(44)

according to Equation (15) and Theorem 4, the PDF of the estimated azimuth ${\widehat{\varphi }}_0$ can be obtained as

$$f_{\varphi }\left({\widehat{\varphi }}_0\right)={\displaystyle {\int }_{-\infty }^{+\infty }{f}_2\left(x,y\left(x,{\widehat{\varphi }}_0\right)\right)\left|\frac{\partial y}{\partial {\widehat{\varphi }}_0}\right|dx}$$

(45)

where

$$\left|\frac{\partial y}{\partial {\widehat{\varphi }}_0}\right|=\left|x\right|{\sec }^2{\widehat{\varphi }}_0.$$

(46)

Hence

$$f_{\varphi }\left({\widehat{\varphi }}_0\right)={\displaystyle {\int }_{-\infty }^{+\infty }{f}_{m_x}\left(x\right)}{f}_{m_y}\left(x\tan {\widehat{\varphi }}_0\right)\left|x\right|{\sec }^2{\widehat{\varphi }}_{0}dx$$

(47)

after simplification,

$$f_{\varphi }\left({\widehat{\varphi }}_0\right)={\sec }^2{\widehat{\varphi }}_0{\displaystyle \sum _{i=1}^2{\displaystyle \sum _{j=1}^2{q}_{ij}\left({\widehat{\varphi }}_0\right)\left(2-\Phi \left(\frac{-{\mu }_{ij}}{{\sigma }_{ij}}\right)-\Phi \left(\frac{{\mu }_{ij}}{{\sigma }_{ij}}\right)\right)}}$$

(48)

where $\Phi \left(x\right)={\displaystyle {\int }_{-\infty }^x\frac{\exp \left(-t^2\right)}{\sqrt{2\pi }}}dt$ and

$$q_{ij}\left({\widehat{\varphi }}_0\right)=r_{ij}\exp \left(s_{ij}{\left(m_{y0}^{}-{\tan }^{}{\widehat{\varphi }}_0{m}_{x0}^{}\right)}^2\right)\{\begin{array}{c}i=1,2\\ j=1,2\end{array}$$

(49)

$$r_{ij}=\frac{a_i{a}_j}{8\pi {\sigma }_{Q_x}{\sigma }_{Q_y}{\sigma }_z^2{\sigma }_y{\sigma }_x}\sqrt{\pi \left(\frac{2b_i{\sigma }_{Q_x}^2{\sigma }_z{\sigma }_y}{{\sigma }_z{\sigma }_y+2b_i{\sigma }_{Q_x}^2}+\frac{2b_j{\sigma }_{Q_y}^2{\sigma }_z{\sigma }_x}{{\sigma }_z{\sigma }_x+2b_j{\sigma }_{Q_y}^2}{\tan }^2{\widehat{\varphi }}_0\right)}$$

(50)

$${\mu }_{ij}=\frac{b_i{\sigma }_{Q_x}^2{\sigma }_y\left({\sigma }_z{\sigma }_x+2b_j{\sigma }_{Q_y}^2\right)m_{x0}+b_j{\sigma }_{Q_y}^2{\sigma }_x\left({\sigma }_z{\sigma }_y+2b_i{\sigma }_{Q_x}^2\right)m_{y0}\tan {\widehat{\varphi }}_0}{b_i{\sigma }_{Q_x}^2{\sigma }_y\left({\sigma }_z{\sigma }_x+2b_j{\sigma }_{Q_y}^2\right)+b_j{\sigma }_{Q_y}^2{\sigma }_x\left({\sigma }_z{\sigma }_y+2b_i{\sigma }_{Q_x}^2\right){\tan }^2{\widehat{\varphi }}_0}$$

(51)

$${\sigma }_{ij}=\frac{1}{2}\sqrt{\frac{\left({\sigma }_z{\sigma }_y+2b_i{\sigma }_{Q_x}^2\right)\left({\sigma }_z{\sigma }_x+2b_j{\sigma }_{Q_y}^2\right)}{b_i{\sigma }_{Q_x}^2{\sigma }_y{\sigma }_z\left({\sigma }_z{\sigma }_x+2b_j{\sigma }_{Q_y}^2\right)+b_j{\sigma }_{Q_y}^2{\sigma }_x{\sigma }_z\left({\sigma }_z{\sigma }_y+2b_i{\sigma }_{Q_x}^2\right){\tan }^2{\widehat{\varphi }}_0}}$$

(52)

$$s_{ij}=\frac{2b_i{b}_j{\sigma }_{Q_x}^2{\sigma }_{Q_y}^2{\sigma }_z{\sigma }_y{\sigma }_x}{b_i{\sigma }_{Q_x}^2{\sigma }_y\left({\sigma }_z{\sigma }_x+2b_j{\sigma }_{Q_y}^2\right)+b_j{\sigma }_{Q_y}^2{\sigma }_x\left({\sigma }_z{\sigma }_y+2b_i{\sigma }_{Q_x}^2\right){\tan }^2{\widehat{\varphi }}_0}.$$

(53)

Clearly, $\Phi \left(x\right)$ is a cumulative distribution function of the standard normal distribution with an integral. For the sake of simplicity, we use function I(x) to replace $\Phi \left(x\right)$ if max|$\Phi \left(x\right)$ − I(x)| is less than 0.01,

$$\mathrm{I}\left(x\right)=\frac{1}{1+\exp \left(-1.701x\right)}.$$

(54)

Hence,

$$f_{\varphi }\left({\widehat{\varphi }}_0\right)={\sec }^2{\widehat{\varphi }}_0{\displaystyle \sum _{i=1}^2{\displaystyle \sum _{j=1}^2{q}_{ij}\left({\widehat{\varphi }}_0\right)\left(2-\mathrm{I}\left(\frac{-{\mu }_{ij}}{{\sigma }_{ij}}\right)-\mathrm{I}\left(\frac{{\mu }_{ij}}{{\sigma }_{ij}}\right)\right)}}.$$

(55)

4.2. Statistical Characteristics of Angle Estimatior

The statistical characteristics of angle estimators are the main indicators to measure the estimation performance. In order to analyze the statistical characteristics of angle estimates, we first study the distribution characteristics of the vector $\overline{m}$, which is obtained by Equation (16). According to Section 4.1, the estimated $\overline{m}$ is symmetrically distributed in space about the accurate value m₀; as a result, the mean value of $\overline{m}$ should be equal to the accurate value m₀, namely

$$E\left[{\overline{m}}_x\right]=m_{x0}$$

(56)

$$E\left[{\overline{m}}_y\right]=m_{y0}$$

(57)

$$E\left[{\overline{m}}_z\right]=m_{z0}.$$

(58)

Meanwhile, let ${\sigma }_{mx}^2$, ${\sigma }_{my}^2$ and ${\sigma }_{mz}^2$ as the variance of m_x, m_y, and m_z and ${\overline{\sigma }}_{mx}^2$, ${\overline{\sigma }}_{my}^2$ and ${\overline{\sigma }}_{mz}^2$ as the variance of ${\overline{m}}_x$, ${\overline{m}}_y$ and ${\overline{m}}_z$ according to the relationship among variables above, we can obtain

$${\overline{\sigma }}_{mx}^2=\frac{{\sigma }_{mx}^2}{K}$$

(59)

$${\overline{\sigma }}_{my}^2=\frac{{\sigma }_{my}^2}{K}$$

(60)

$${\overline{\sigma }}_{mz}^2=\frac{{\sigma }_{mz}^2}{K}.$$

(61)

Obviously, when the number of snapshots K tends to infinity, ${\overline{\sigma }}_{mx}^2$, ${\overline{\sigma }}_{my}^2$ and ${\overline{\sigma }}_{mz}^2$ are all close to 0.

According to above conclusion, using Chebyshev’s theorem, the following inequalities are established:

$$P\left\{\left|{\overline{m}}_x-m_{x0}\right|\ge {\epsilon }_1\right\}\le \frac{{\overline{\sigma }}_{mx}^2}{{\epsilon }_1^2}$$

(62)

$$P\left\{\left|{\overline{m}}_y-m_{y0}\right|\ge {\epsilon }_2\right\}\le \frac{{\overline{\sigma }}_{my}^2}{{\epsilon }_2^2}$$

(63)

$$P\left\{\left|{\overline{m}}_z-m_{z0}\right|\ge {\epsilon }_3\right\}\le \frac{{\overline{\sigma }}_{mz}^2}{{\epsilon }_3^2}$$

(64)

where ε₁, ε₂, and ε₃ are arbitrary positive numbers.

Clearly, the probability for ${\overline{m}}_x$, ${\overline{m}}_y$ and ${\overline{m}}_z$ deviate from m_x0, m_y0, and m_z0 tends to 0 respectively if the number of snapshots K is very large. Combining Equations (16) and (17), we can get

$$\underset{K\to \infty }{\lim }\arccos \left(\frac{{\overline{m}}_z}{\sqrt{{\overline{m}}_x^2+{\overline{m}}_y^2+{\overline{m}}_z^2}}\right)=\arccos \left(\frac{m_{z0}}{\sqrt{m_{x0}^2+m_{y0}^2+m_{z0}^2}}\right)={\theta }_0$$

(65)

$$\{\begin{array}{c}\underset{K\to \infty }{\lim }\arccos \left(\frac{{\overline{m}}_x}{\sqrt{{\overline{m}}_x^2+{\overline{m}}_y^2}}\right)=\arccos \left(\frac{m_{x0}}{\sqrt{m_{x0}^2+m_{y0}^2}}\right)={\varphi }_0\left({\overline{m}}_y\ge 0\right)\\ \underset{K\to \infty }{\lim }-\arccos \left(\frac{{\overline{m}}_x}{\sqrt{{\overline{m}}_x^2+{\overline{m}}_y^2}}\right)=\arccos \left(\frac{m_{x0}}{\sqrt{m_{x0}^2+m_{y0}^2}}\right)={\varphi }_0\left({\overline{m}}_y\le 0\right)\end{array}.$$

(66)

Therefore, we can conclude that our algorithm is an asymptotic unbiased estimation, and the estimated value will become closer to the true value as the number of snapshots increases.

4.3. Analysis of Algorithm Complexity

MUSIC is a DOA estimation method based on subspace, which uses eigenvalue decomposition to the covariance matrix and searches peaks by spatial spectrum. Its advantage is that D different signals can be estimated at the same time. However, the disadvantage is that the N-element array must be used and results in heavy computation load, where N > D.

The computation load of the interferometer algorithm is relatively small, but the amount of operation increases significantly with the number of array elements. In addition, if the correction for amplitude and phase errors is considered, the amount of operation would further increase.

Compared with above these two algorithms, the proposed algorithm based on a tripole vector antenna has a significant advantage in computation. When K snapshots are used for estimation, the new algorithm here only needs to complete multiplication in 6K − 1 times, addition in 6K − 1 times, division in 2 times, square root in 2 times, and anti-cosine function in 2 times.

In details, the square root operation is calculated by Newton iterative method, and the iterative formula is

$$x_{n+1}=\frac{b}{2x_n}+\frac{x_n}{2}\begin{array}{cc}& x_0=b\end{array}$$

(67)

where b is a radicand. Since A means iteration times and x_A is the square root of b, there are A times of addition, A times of multiplication, and A times of division in one square root operation. On the other hand, we calculate the inverse cosine function by Taylor’s expansion

$$\arccos x=\frac{\pi }{2}-{\displaystyle \sum _{n=0}^M\left(\frac{\left(2n\right)!}{2^{{}_{2n}}{\left(n!\right)}^2}\right)\frac{x^{2n+1}}{2n+1}}=\frac{\pi }{2}-{\displaystyle \sum _{n=0}^M{a}_n{x}^{2n+1}}\left(\left|x\right|\le 1\right).$$

(68)

So, the amount of an anti-cosine function computation consists of (M + 1) (M + 2) times of multiplication, $n+1$ times of addition.

In conclusion, the computation load of the proposed algorithm includes M² + 3M + 1 + 2A + 6K times of multiplication, n + 2A + 6K times of addition, 2 + 2A times of division. The computation load of each algorithm can be seen in

Table 1. The computational complexity of the algorithm. MUSIC: Multiple Signal Classification.

Algorithm	Times of Multiplication	Times of Addition	Times of Division
MUSIC [11]	(2N − 1) (N − D) IQ + NK	(2N − 1) (N − D) IQ + (2N − 1) K	1
Interferometer [9]	8 + A + 4 (M + 1) (M + 2) + 2K + 0.5 (N − 2) Klog2K	A + 1 + 2K + 4 (n + 1) + Klog2K	6 + A
Proposed Algorithm	M2 + 3M + 1 + 2A + 6K	n + 2A + 6K	2 + 2A

.

Where Q is spectral peak search times of azimuth and I is that of elevation.

It can be seen from Table 1 that compared with existing methods, our algorithm of angle estimation based on a single tripole vector antenna has obvious advantages in computational complexity.

5. Simulations

In this section, we present some numerical results to assess the performance of our algorithm. For a comparison, we investigate the performance of MUSIC [11], Interferometer [9], and our proposed algorithm.

In the following simulations, we consider an L-shaped array consisting of three parallel dipoles used for the interferometer, and three triple vector antennas used for both MUSIC and the proposed algorithm. The spectral peak search interval of the MUSIC algorithm is 0.1 degree. In addition, we assume that the DOA of the incident EM wave is ${\theta }_0={\phi }_0={30}^{\circ }$, and the frequency of the EM wave is 1 GHz.

For a statistical comparison, we define the root mean square error (RMSE) as

$$\mathrm{RMSE}=\sqrt{\frac{1}{J}{\displaystyle \sum _{j=1}^J\left[{\left({\varphi }_0-{\widehat{\varphi }}_{0j}\right)}^2+{\left({\theta }_0-{\widehat{\theta }}_{0j}\right)}^2\right]}}$$

(69)

where ${\varphi }_0$ and ${\theta }_0$ are the actual values of azimuth and elevation angle, respectively. In addition, ${\widehat{\varphi }}_{0j}$ and ${\widehat{\theta }}_{0j}$ are the estimated values of azimuth and elevation angle in the j-th Monte Carlo trial where j = 1, 2, ..., J.

We also consider the statistical characteristics of the estimation error for the proposed algorithm. The estimation error is defined as the difference between the estimated value and true value, that is

$$\{\begin{array}{l}\Delta \theta ={\theta }_0-{\widehat{\theta }}_0\\ \Delta \varphi ={\varphi }_0-{\widehat{\varphi }}_0\end{array}.$$

(70)

Their mean and variance are

$$\{\begin{array}{l}E\left[\Delta \theta \right]=\frac{1}{J}{\displaystyle \sum _{j=1}^J\left({\theta }_0-{\widehat{\theta }}_{0j}\right)}\\ D\left[\Delta \theta \right]=E^2\left[\Delta \theta \right]-E{\left[\Delta {\theta }^2\right]}^{}\end{array}$$

(71)

$$\{\begin{array}{l}E\left[\Delta \varphi \right]=\frac{1}{J}{\displaystyle \sum _{j=1}^J\left({\varphi }_0-{\widehat{\varphi }}_{0j}\right)}\\ D\left[\Delta \varphi \right]=E^2\left[\Delta \varphi \right]-E{\left[\Delta {\varphi }^2\right]}^{}\end{array}.$$

(72)

In addition, the signal-to-noise ratio is defined as

$$\mathrm{SNR}=10\lg \frac{P_s}{P_n}$$

(73)

where P_s and P_n represent the effective power of signal and noise, respectively.

5.1. Simulation 1: Statistical Characteristics of Angle Error

In the first simulation, we mainly introduce the statistical characteristics of the estimation error of the proposed algorithm according to Equations (71) and (72). Supposing J = 200 and SNR = 10 dB, the results are shown in Figure 4.

In Figure 4, the mean of the estimation error is almost always 0 even in the case of small number of snapshots, which illustrates that the estimation method can be considered unbiased and the error mean is not sensitive to the number of snapshots. However, the variance of estimation error is relatively large when the number of snapshots is small. Increasing the number of snapshots can reduce the variance. Furthermore, the variance of azimuth is always larger than that of elevation angle.

5.2. Simulation 2: Comparison of RMSE Performance Versus Different SNR

The second simulation illustrates the comparison of the performance among three kinds of DOA estimation algorithms, which are phase interferometer, MUSIC, and the proposed tripole vector antenna method. In the simulation, K = 1000, J = 200, and SNR varies from 0 to 10 dB. The results are shown in Figure 5.

It depicts in Figure 5 that the MUSIC algorithm has the best performance followed by the proposed algorithm here. In addition, due to the polarization mismatch between the EM wave and the antenna for interferometer, the input SNR loss is 3 dB and leads to an estimation error almost twice that of the proposed method.

5.3. Simulation 3: Comparison of RMSE Performance Versus Different Snapshots

The third simulation compares the performance of the three algorithms under different snapshots. In this simulation, SNR = 10 dB, J = 200, and K varies from 100 to 1000. The results are shown in Figure 6.

As illustrated in Figure 6, when the number of snapshots is the same, the performance of MUSIC is the best, followed by the proposed algorithm and the interferometer. This result is very similar to the simulation result in Section 5.2, which shows that increasing number of snapshots is equal to improving the SNR for the proposed algorithm.

5.4. Simulation 4: RMSE for Different DOAs in Degrees

In order to fully demonstrate the estimated performance of the proposed algorithm at various angles, we simulated all angles under the condition of 10 dB and 1000 snapshots, and the results are as follows.

According to Figure 7, estimate errors are different in various incident angles. In the situation of invariant azimuth, as a decline of elevation, RMSE decreases; when the elevation stays invariant, the RMSE varies periodically with the change of azimuth, the RMSE is relatively small when azimuth is kπ/2, and it is relatively large when kπ/2 + π/4. The above simulation results are demonstrated from Equations (17) and (18). According to Equation (17), the elevation angle depends on $\frac{{\overline{m}}_x^2+{\overline{m}}_y^2}{{\overline{m}}_z^2}$; according to Equation (18), azimuth depends on $\frac{{\overline{m}}_x^2}{{\overline{m}}_x^2+{\overline{m}}_y^2}$, which has nothing to do with ${\overline{m}}_z$.

5.5. Simulation 5: Algorithm Complexity

In this simulation, we record the running time of different algorithms in MATLAB, which is obtained by taking the average of the running time of multiple Monte Carlo experiments. The version of MATLAB is 2018b. The CPU of the computer running the experiment is i7-7700HQ with a clock frequency of 2.8 GHz. The simulation results are shown in Figure 8.

Figure 8 presents that the speed of the proposed algorithm and the interferometer algorithm are both much faster than that of the MUSIC algorithm. In detail, the running time of our algorithm is at least 5 orders of magnitude lower than that of MUSIC algorithm. Furthermore, as the number of snapshots increases, the calculation time of the interferometer algorithm is approximately 40 to 200 times that of the proposed algorithm. Obviously, the significant advantage of the proposed tripole antenna algorithm is the low computational complexity.

6. Conclusions

We investigate the DOA estimation with a tripole vector antenna for a single nonlinearly polarized EM wave. Unlike traditional dipole antennas, tripole vector antennas can completely capture the electric field of EM waves, which makes it possible to estimate the DOA directly through the geometric relationship between the signal propagation direction and the tripole vector antenna. In the proposed method, the antenna array is not required so that hardware cost is reduced. In addition, many complex operations are not required, such as spectrum search and matrix eigenvalue decomposition, which means the proposed algorithm has extremely fast operation speed with satisfactory accuracy. In the future, we will modify the new algorithm to enlarge the suitable application fields. In addition, combining the proposed algorithm with other schemes to obtain better performance in DOA estimation will also be considered.