Direction of Arrival Estimation Based on the Multistage Nested Wiener Filter

Abstract

A novel direction of arrival (DOA) estimation technique based on data level and order recursive Multistage Nested Wiener Filters (MSNWF) which is used in adaptive beamforming for subarray signal is proposed in this paper. The two subarrays using the same array geometry are used to form a signal whose phase relative to the reference signal is a function of the DOA. The DOA is estimated by calculating the phase-shift between the reference signal and its phase-shifted version. The performance of this DOA estimation technique is significantly improved due to the application of order recursive MSNWF for the rejection of interference signals. The computation of the proposed method is simple, and the number of detectable signal sources could exceed the number of antenna elements.

1. Introduction

In the last two decades, smart antenna has been widely used in many applications such as radar, sonar, and wireless communication systems [1]. It is also utilized in tracking [2, 3], localization [4, 5], intelligent transportation [6], ultra-wideband wireless sensor networks [7], array calibration [8], scatter cluster model [9], and antijamming [10]. For example, the multiple input multiple output (MIMO) radar utilizes multiple sensor array antennas to simultaneously transmit and receive diverse waveforms, which estimates the signal parameters to locate and track the target [2]. Distributed sensor networks have been used for enhancing signal to noise ratios for space-time localization and tracking of remote objects using phased array antennas [4]. Radio frequency identification (RFID) is widely used for electronically identifying, locating, and tracking products, animals, and vehicles as a very valuable business and technology tool [5]. Vehicular ad hoc networks (VANETs) could be a benefit to the traffic safety and efficiency [6]. The performance of array processing algorithms is improved by the sensor array location error calibration, which made the algorithms insensitive to the model uncertainties and deterministic signals with unknown waveforms [8]. The performance of the wireless communication system is evaluated based on scatter cluster models by estimating the corresponding parameters [9]. Array sensor and subarray adaptive beamforming techniques obtain the best antijamming performance widely used in GNSS receivers [10], active radar, and sonar [11]. In these sensor networks implication systems and scenarios, direction of arrival (DOA) is an important parameter that is needed to be estimated to determine the direction of the located and tracked target or the position of the sensor nodes.

Considerable research efforts have been made in the DOA estimation and various array signal process techniques for DOA estimation have been proposed [12–18]. The most commonly used DOA estimation techniques include (1) spectrum based methods, such as Bartlett [14] and Capon [15]; (2) subspace-based algorithm, such as multiple signal classification (MUSIC) [16]; (3) parametric methods, such as estimation of signal parameters via rotational invariance technique (ESPRIT) [17]. In Capon techniques, the DOAs are determined by finding the directions in which their antenna response vectors lead to peaks in the spectrum formed by the covariance matrix of the observation vectors. Thus, the capacity of this DOA estimation technique is less than the number of antenna elements bounded by the covariance matrix of the observation vectors. In MUSIC techniques, the DOAs of target signals are determined by finding the directions in which their antenna response vectors lead to peaks in the MUSIC spectrum formed by the eigenvectors of the noise subspace. Thus, the capacity of this DOA estimation is equal to the rank of the reciprocal subspace of the selected noise subspace and is also less than the number of the antenna elements. In ESPRIT techniques, two virtual subarrays structures are proposed to obtain two signal subspaces. The eigenvectors of the relevant signal subspaces are rotated for the DOAs of the signals. As a result, the capacity of DOA estimation using ESPRIT is bounded by the number of subarrays. The application of the above techniques is limited to cases where the number of signal sources is less than that of antenna elements. These techniques require subspace estimation, eigendecomposition, and the computation of covariance matrix inversion which leads to high computational complexity, and they are thereby limited to the applications where fast DOA estimation is required. Furthermore, in the presence of interference, these techniques need to estimate the DOAs of all the target signals and interference, which decreases the accuracy of DOA estimation.

The application of adaptive beamforming in DOA estimation has become the research focus on interference existence [18]. In [18], Wang et al. developed a new structure of DOA estimation based on subarray beamforming. This technique has clear advantage on the DOA estimation when interference exists, but it still needs the computation of matrix inversion which is not easy to be applied to a practical system. Based on this structure, a DOA estimation technique based on Multistage Nested Wiener Filter (MSNWF) [19–25] is proposed in [26]. In [26], Yu stated an original MSNWF algorithm [27] to estimate the DOAs, which used a filter and blocking matrix to avoid the calculation of covariance matrix inversion. In this technique, however, it cannot calculate the coefficients of Wiener filter in backward recursion if the forward recursion does not finish the calculation of the match filter and blocking matrix. And the mean squared error (MSE) can not be determined when adding a new stage, except for the last stage.

In this paper, a data level order recursive MSNWF DOA estimation technique that uses a reference signal is proposed in detail. This DOA estimation technique uses two subarray adaptive beamformers based on the data level order recursive MSNWF to construct the same array geometry for forming the phase-shift and rejecting interference at same time. The DOAs of the target signals are estimated from the phase-shifts by using reference signal after the rejection of interference. Therefore, the performance of DOA estimation is significantly improved. This technique can be widely used for the implementation of hardware systems such as wireless communication system, active radar, sonar, and space-time adaptive process (STAP) systems [28, 29].

The advantages of the data level order recursive MSNWF DOA estimation are as follows. (1) Since the use of data level order recursive MSNWF in this DOA estimation technique realizes the subspace eigendecomposition, computation of inversion of covariance matrix becomes unnecessary and thus reduces the complexity of computation; the data level MSNWF DOA estimation technique can be easily applied in hardware platform. (2) An orthogonal basis for the Krylov subspace spanned cross correlation vector and covariance matrix improves the computational efficiency when calculating the weight vector of the match filter. And the order recursion could update the weight vector of the match filter and the MSE at new stage.

The paper is organized as follows. In Section 2, the signal model is described. In Section 3, the structure of the data level order recursive MSNWF DOA estimation system, MSNWF based adaptive beamforming including data level recursion of match filters and order recursion, and DOA calculation of the proposed method are presented. Design examples and simulation results are given in Section 5, and conclusions are drawn in Section 6.

2. Signal Model

Consider a uniform linear array (ULA) system that uses M elements with adjacent element spacing d, deployed at a base station. Assume that K narrowband signals and P unknown interference sources are received at the ULA with different DOAs ${\theta }_k$ , $k=\mathrm{1,2},\dots ,K+P$ .

Using complex envelope representation, the received signals can be expressed by

$$\begin{array}{cc}& \mathbf{x}\left(t\right)={\displaystyle \sum }_{k=\mathrm{1}}^{K+P}\mathbf{a}\left({\theta }_k\right)s_k\left(t\right)+\mathbf{n}\left(t\right),\end{array}$$

(1)

where $s_k\left(t\right)$ denotes the kth signal component, $k=\mathrm{1,2},\dots ,K$ denotes the target components, and $k=K+\mathrm{1},K+\mathrm{2},\dots ,K+P$ are interference components. The $\mathbf{a}\left({\theta }_k\right)$ in (1) denotes the steering vector of the array in direction ${\theta }_k$ , which is given by

$$\begin{array}{cc}& \mathbf{a}\left({\theta }_k\right)={\left[\mathrm{1},e^{-j\mathrm{2}\pi d/\sin \left({\theta }_k\right)},\dots ,e^{-j\mathrm{2}\pi d(M-\mathrm{1})\sin \left({\theta }_k\right)}\right]}^T\end{array}$$

(2)

and $\mathbf{n}\left(t\right)$ denotes the noise vector with zero mean and cross covariance

$$\begin{array}{cc}& E\left[\mathbf{n}\left(t_{\mathrm{1}}\right){\mathbf{n}}^H\left(t_{\mathrm{2}}\right)\right]={\sigma }^{\mathrm{2}}\delta \left(t_{\mathrm{1}}-t_{\mathrm{2}}\right)\mathbf{I},\end{array}$$

(3)

where I is the identity matrix.

Suppose that the received vector $\mathbf{x}\left(t\right)$ is sampled at n, $n=\mathrm{1,2},\dots ,L$ , and the received signal can be expressed by (4) in the matrix notation. Consider

$$\begin{array}{cc}& \mathbf{X}=\mathbf{A}\left(\theta \right)\mathbf{S}+\mathbf{N},\end{array}$$

(4)

where X and N are $M\times L$ matrices:

$$\begin{array}{cc}& \mathbf{X}=\left[\mathbf{x}\left(\mathrm{1}\right),\mathbf{x}\left(\mathrm{2}\right),\dots ,\mathbf{x}\left(L\right)\right],\\ \mathbf{N}=\left[\mathbf{n}\left(\mathrm{1}\right),\mathbf{n}\left(\mathrm{2}\right),\dots ,\mathbf{n}\left(L\right)\right].\end{array}$$

(5)

$\mathbf{A}\left(\theta \right)$ is a $M\times K$ matrix as follows:

$$\begin{array}{cc}& \mathbf{A}\left(\theta \right)=\left[\mathbf{a}\left({\theta }_{\mathrm{1}}\right),\mathbf{a}\left({\theta }_{\mathrm{2}}\right),\dots ,\mathbf{a}\left({\theta }_K\right)\right].\end{array}$$

(6)

And S is a $K\times L$ matrix

$$\begin{array}{cc}& \mathbf{S}=\left[\mathbf{s}\left(\mathrm{1}\right),\mathbf{s}\left(\mathrm{2}\right),\dots ,\mathbf{s}\left(L\right)\right].\end{array}$$

(7)

3. MSNWF DOA Estimation

Compared with the SBDOA estimation technique stated in [8], the proposed MSNWF DOA estimation technique in this paper uses the same uniform linear antenna array at the receiving end and the geometry of the array is similar to that used in ESPRIT techniques. The antenna array is decomposed into two equal-sized subarrays, where the two subarrays are used in conjunction with two subarray MSNWF adaptive beamformers to obtain an optimal estimation of a phase-shift reference signal whose phase relative to that of the reference signal is a function of the target DOA. The target DOA is then computed from the estimated phase-shift between the reference signal ${\mathbf{r}}_k$ and the phase-shifted reference signal $e^{j{\varphi }_k}{\mathbf{r}}_k$ . In order to avoid the inversion computation of covariance matrix when getting the optimal weight vector of the beamformer, the two beamformers as in Figure 1 in [18] are replaced with Multistage Nested Wiener Filters. The block diagram of the MSNWF DOA estimation system is illustrated in Figure 1.

Figure 1

Block diagram of the MSNWF DOA estimation system.

3.1. Subarray Signal Formation

Consider that the array is composed of a ULA of M element as a receiver and decomposed into two sets of $M-\mathrm{1}$ element virtual subarrays, A and B . The downconverted baseband signal received by the mth, $m=\mathrm{1,2},\dots ,M$ element of the antenna array is expressed by

$$\begin{array}{cc}& x_m\left(n\right)={\displaystyle \sum }_{k=\mathrm{1}}^{K+P}{e}^{j\left(m-\mathrm{1}\right){\varphi }_k}{s}_k\left(n\right).\end{array}$$

(8)

The vectors of the A and B are given by

$$\begin{array}{l}{\mathbf{y}}_A={[x_1(n),x_2(n),\dots ,x_{M-1}(n)]}^T,\hfill \\ {\mathbf{y}}_B={[x_2(n),x_3(n),\dots ,x_M(n)]}^T,\hfill \end{array}$$

(9)

respectively. Let

$$\begin{array}{cc}& \mathbf{b}\left({\theta }_k\right)={\left[\mathrm{1},e^{j{\varphi }_k},\dots ,e^{j{\varphi }_k\left(M-\mathrm{2}\right)}\right]}^T\end{array}$$

(10)

and then the subarray signals ${\mathbf{y}}_A$ and ${\mathbf{y}}_B$ can be written as

$$\begin{array}{l}{\mathbf{y}}_A(n)=\sum _{k=1}^{K+P}\mathbf{b}({\theta }_k)s_k(n)+{\mathbf{n}}_A(n),\hfill \\ {\mathbf{y}}_B(n)=\sum _{k=1}^{K+P}{e}^{j{\phi }_k}\mathbf{b}({\theta }_k)s_k(n)+{\mathbf{n}}_B(n),\hfill \end{array}$$

(11)

where vectors ${\mathbf{n}}_A\left(n\right)$ and ${\mathbf{n}}_B\left(n\right)$ are the background noise at the subarray, respectively. The phase-shift factor between the kth components of signals ${\mathbf{y}}_A\left(n\right)$ and ${\mathbf{y}}_B\left(n\right)$ which forms the kth signal is given by

$$\begin{array}{cc}& e^{j{\varphi }_k}=e^{-j\mathrm{2}\pi d\sin \left({\theta }_k\right)/\lambda }.\end{array}$$

(12)

Sampling ${\mathbf{y}}_A\left(n\right)$ and ${\mathbf{y}}_B\left(n\right)$ obtains

$$\begin{array}{l}{\text{Y}}_A=[y_A(1),{\mathbf{y}}_A(2),\dots ,{\mathbf{y}}_A(L)],\hfill \\ {\text{Y}}_B=[{\mathbf{y}}_B(1),{\mathbf{y}}_B(2),\dots ,{\mathbf{y}}_B(L)].\hfill \end{array}$$

(13)

3.2. Recursion Algorithm of MSNWF

3.2.1. Data Level Recursion of Match Filters

In the Wiener filter, the estimation of the desired signal $d_{\mathrm{0}}(n)$ from an observation vector ${\mathit{x}}_0(n)$ is optimal in the minimum mean square error (MMSE) sense. The weight vector ${\mathit{w}}_{\mathit{X}0}$ of the Wiener filter can be obtained via solving the following Wiener-Hopf equations

$$\begin{array}{cc}& {\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{w}}_{\mathit{x}\mathrm{0}}={\mathit{r}}_{\mathit{x}\mathrm{0},\mathit{d}\mathrm{0}},\end{array}$$

(14)

where ${\mathit{R}}_{\mathit{x}0}$ is the covariance matrix of observation vector ${\mathbf{x}}_{\mathrm{0}}(n)$ and ${\mathit{r}}_{\mathit{x}0}$ is the cross correlation vector between the observation vector ${\mathbf{x}}_{\mathrm{0}}\left(n\right)$ and the desired signal $d_{\mathrm{0}}\left(n\right)$ . The covariance matrix ${\mathit{R}}_{\mathit{x}0}$ cannot be readily estimated, if ${\mathbf{x}}_{\mathrm{0}}\left(n\right)$ is of high dimension. Based on this, Goldstein and Reed proposed that if the observation signal ${\mathbf{x}}_{\mathrm{0}}\left(n\right)$ is prefiltered by a full-rank matrix $\mathit{T}\in {\mathbb{C}}^{M\times M}$ to get a new observation signal ${\mathit{z}}_1\left(n\right)=\mathit{T}{\mathit{x}}_0\left(n\right)$ , then the weight factor ${\mathit{w}}_{\mathit{z}1}$ of Wiener filter is used to estimate the desired signal $d_0\left(n\right)$ from ${\mathit{z}}_1\left(n\right)$ results in the same MSE [19–21].

The assumed full-rank prefilter matrix can be chosen as

$$\begin{array}{cc}& {\mathit{T}}_{\mathrm{1}}=\left[\begin{array}{c}{\mathit{h}}_{\mathrm{1}}^H\\ {\mathit{B}}_{\mathrm{1}}\end{array}\right],\end{array}$$

(15)

where H is the complex conjugate transpose operator. Thus,

$$\begin{array}{cc}& {\mathit{z}}_{\mathrm{1}}\left(n\right)=\left[\begin{array}{c}{\mathit{h}}_{\mathrm{1}}^H{\mathit{x}}_{\mathrm{0}}\left(n\right)\\ {\mathit{B}}_{\mathrm{1}}{\mathit{x}}_{\mathrm{0}}\left(n\right)\end{array}\right]=\left[\begin{array}{c}{d}_{\mathrm{1}}\left(n\right)\\ {\mathit{x}}_{\mathrm{1}}\left(n\right)\end{array}\right],\end{array}$$

(16)

where ${\mathit{B}}_1$ is referring to the blocking matrix, ${\mathit{B}}_1{\mathit{h}}_1=0$ and ${\mathit{h}}_1={\mathit{r}}_{\mathit{x}0\mathit{d}0}/{‖{\mathit{r}}_{\mathit{x}0\mathit{d}0}‖}_2$ .

The solution of the Wiener-Hopf equations relative to the transformed system is

$$\begin{array}{cc}& {\mathit{w}}_{\mathit{z}\mathrm{1}}={\mathit{R}}_{\mathit{z}\mathrm{1}}^{-\mathrm{1}}{\mathit{r}}_{\mathit{z}\mathrm{1},\mathit{d}\mathrm{0}}={\alpha }_{\mathrm{1}}\left[\begin{array}{c}\mathrm{1}\\ -{\mathit{R}}_{\mathit{x}\mathrm{1}}^{-\mathrm{1}}{\mathit{r}}_{\mathit{x}\mathrm{1},\mathit{d}\mathrm{1}}\end{array}\right],\end{array}$$

(17)

where ${\mathit{R}}_{\mathit{z}1}$ is the covariance matrix of the new observation signal ${\mathit{z}}_1\left(n\right)$ , ${\alpha }_1={‖{\mathit{r}}_{\mathit{x}0,\mathit{d}0}‖}_2/({\sigma }_{d1}^2-{\mathit{r}}_{\mathit{x}1,\mathit{d}1}^H{\mathit{R}}_{\mathit{x}1}^{-1}{\mathit{r}}_{\mathit{x}1,\mathit{d}1})$ , ${\mathit{r}}_{\mathit{z}1,\mathit{d}1}$ is the cross correlation between ${\mathit{d}}_1\left(n\right)$ and ${\mathit{z}}_1\left(n\right)$ , and ${\mathit{R}}_{\mathit{x}1}={\mathit{B}}_1{\mathit{R}}_{\mathit{x}0}{\mathit{B}}_1^H$ , ${\sigma }_{d1}^2={\mathit{h}}_1^H{\mathit{R}}_{\mathit{x}0}{\mathit{h}}_1$ , ${\mathit{r}}_{\mathit{x}1,\mathit{d}1}={\mathit{B}}_1{\mathit{R}}_{\mathit{x}0}{\mathit{h}}_1$ .

This process produces a new vector Wiener filter, which estimates the signal $d_1\left(n\right)$ from the observation vector ${\mathit{x}}_1\left(n\right)$ , and a scalar Wiener filter is followed. Repeating this process, a nested structure can be obtained, which is defined as the original MSNWF [19–21].

In the original MSNWF, the new desired signal $d_i\left(n\right)$ at the output of the ith stage can be expressed as

$$\begin{array}{cc}& d_i\left(n\right)={\mathit{x}}_{\mathrm{0}}\left(n\right)\left({\displaystyle \prod }_{k=\mathrm{1}}^{i-\mathrm{1}}{\mathit{B}}_k^H\right){\mathit{h}}_i={\mathit{x}}_{i-\mathrm{1}}\left(n\right){\mathit{t}}_i.\end{array}$$

(18)

According to (18), a filter ${\mathit{t}}_i$ is used to replace the ith stage Wiener filter as Figure 2 shows, which could be simply the cross correlation between the new observation ${\mathit{x}}_i\left(n\right)$ and the new desired signal $d_i\left(n\right)$ .

Figure 2

First stage of the original MSNWF.

The new observation vector in Figure 3 is expressed as

$$\begin{array}{cc}& \mathit{d}\left(n\right)={\left[d_{\mathrm{1}}\left(n\right),d_{\mathrm{2}}\left(n\right),\dots ,d_N\left(n\right)\right]}^T\end{array}$$

(19)

which proved that this observation vector has a tridiagonal covariance matrix [21].

Figure 3

Match filter bank structure of MSNWF.

Therefore, the new desired signal $d_i\left(n\right)$ can be seen that it is the output of an N length filter ${\mathit{t}}_i$

$$\begin{array}{cc}& {\mathit{t}}_i=\left({\displaystyle \prod }_{k=\mathrm{1}}^{i-\mathrm{1}}{\mathit{B}}_k^H\right){\mathit{h}}_i.\end{array}$$

(20)

The filter ${\mathbf{t}}_i$ is used to recover all the information of ${\mathit{x}}_{i-1}\left(n\right)$ via $d_{i-1}\left(n\right)$ . The output $d_i\left(n\right)$ is gotten by the filter ${\mathit{t}}_{i+1}$ ; thus, $d_i\left(n\right)$ is correlated with $d_{i-1}\left(n\right)$ and $d_{i+1}\left(n\right)$ . However, $d_{i+1}\left(n\right)$ is from the blocking matrix ${\mathbf{B}}_{i+\mathrm{1}}$ , which is not correlated with $d_{i-1}\left(n\right)$ . Therefore, $d_i\left(n\right)$ is only correlated with its two neighbors. And it is also required to be maximally correlated with $d_{i-1}\left(n\right)$ . Considering the orthogonality conditions, the maximal correlation results in an optimization problem [25] as follows:

$$\begin{array}{cc}& \hspace{0.17em}{\mathit{t}}_i=\arg \underset{\mathit{t}}{\max }\hspace{0.17em}E\left[d_i\left(n\right)d_{i-\mathrm{1}}^{*}\left(n\right)\right],\\ \hspace{0.17em}\text{s}.\text{t}.:{\mathit{t}}^H\mathit{t}=\mathrm{1},\hspace{1em}\hspace{1em}{\mathit{t}}^H{\mathit{t}}_k=\mathrm{0},\hspace{1em}k=\mathrm{1,2},\dots ,i-\mathrm{1}.\end{array}$$

(21)

Using Lagrange multipliers, the solution of (21) is

$$\begin{array}{cc}& {\mathit{t}}_i=\frac{\left({\prod }_{k=\mathrm{1}}^{i-\mathrm{1}}{\mathit{P}}_k\right){\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{t}}_{i-\mathrm{1}}}{{‖\left({\prod }_{k=\mathrm{1}}^{i-\mathrm{1}}{\mathit{P}}_k\right){\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{t}}_{i-\mathrm{1}}‖}_{\mathrm{2}}},\end{array}$$

(22)

where ${\mathit{P}}_k={\mathit{I}}_N-{\mathit{t}}_i{\mathit{t}}_i^H$ .

Herein, if ${\mathit{B}}_i$ is assumed to be equal to ${\mathit{P}}_i$ , the filters ${\mathit{t}}_i$ are an orthonormal basis for the Krylov subspace generated by ${\mathit{r}}_{\mathit{x}0,\mathit{d}0}$ and ${\mathit{R}}_{\mathit{x}0}$ [22]. Therefore, the result of recursion of the MSNWF can be obtained without ${\mathit{B}}_i$ .

At ith stage, let

$$\begin{array}{cc}& {\mathit{u}}_i={\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{t}}_{i-\mathrm{1}}.\end{array}$$

(23)

The filters ${\mathit{t}}_i$ of the recursion are computed as

$$\begin{array}{cc}& {\mathit{t}}_i={\mathit{u}}_i-\left({\mathit{t}}_{i-\mathrm{1}}^H{\mathit{u}}_i\right){\mathit{t}}_{i-\mathrm{1}}-\left({\mathit{t}}_{i-\mathrm{2}}^H{\mathit{u}}_i\right){\mathit{t}}_{i-\mathrm{2}}.\end{array}$$

(24)

In the recursion calculation process, the filters ${\mathit{t}}_i$ are calculated which does not need ${\mathit{B}}_i$ and the inversion of covariance matrix, and this reduces the complexity of computation. The calculation of ${\mathit{t}}_i$ only needs the last two members, which also reduces the complexity of computation.

3.2.2. Order Recursion

At the stage $\left(M-1\right)$ of the MSNWF, the orthogonal basis composed by the matcher filters ${\mathit{t}}_i$ is expressed as

$$\begin{array}{cc}& {\mathit{T}}^{\left(M-\mathrm{1}\right)}=\left[{\mathit{t}}_{\mathrm{1}},{\mathit{t}}_{\mathrm{2}},\dots ,{\mathit{t}}_{M-\mathrm{1}}\right].\end{array}$$

(25)

The new observation vector obtained from the recursion calculation can be written as

$$\begin{array}{lll}{\text{d}}^{(M-1)}(n)\hfill & =\hfill & [d_1(n),d_2(n),\dots ,d_{M-1}(n)]\hfill \\ \hfill & =\hfill & [{\mathbf{t}}_1^H{\mathbf{x}}_0(n),{\mathbf{t}}_2^H{\mathbf{x}}_0(n),\dots ,{\mathbf{t}}_{M-1}^H{\mathbf{x}}_0(n)]\hfill \\ \hfill & =\hfill & {\left({\mathbf{T}}^{(M-1)}\right)}^H{\mathbf{x}}_0(n).\hfill \end{array}$$

(26)

The covariance matrix can be written as

$$\begin{array}{cc}& {\mathit{R}}_{\mathit{d}}^{\left(M-\mathrm{1}\right)}={\left({\mathit{T}}^{\left(M-\mathrm{1}\right)}\right)}^H{\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{T}}^{\left(M-\mathrm{1}\right)}.\end{array}$$

(27)

The recursion coefficients are the components of Wiener filter coefficients as (28) which is used to estimate $d_0\left(n\right)$ from ${\mathit{d}}^{\left(M-1\right)}\left(n\right)$

$$\begin{array}{cc}& {\mathit{w}}_{\mathit{d}}^{\left(M-\mathrm{1}\right)}={\left({\mathit{R}}_{\mathit{d}}^{\left(M-\mathrm{1}\right)}\right)}^{-\mathrm{1}}{\mathit{r}}_{\mathit{d},\mathit{d}\mathbf{0}}^{\left(M-\mathrm{1}\right)}={\left({\mathit{R}}_{\mathit{d}}^{\left(M-\mathrm{1}\right)}\right)}^{-\mathrm{1}}{\left({\mathit{T}}^{\left(M-\mathrm{1}\right)}\right)}^H{\mathit{r}}_{\mathit{x}\mathbf{0},\mathit{d}\mathbf{0}}.\end{array}$$

(28)

Then, the coefficients of MSNWF can be expressed as

$$\begin{array}{cc}& {\mathit{w}}_{\mathrm{0}}^{\left(M-\mathrm{1}\right)}={\mathit{T}}^{\left(M-\mathrm{1}\right)}{\mathit{w}}_{\mathit{d}}^{\left(M-\mathrm{1}\right)}.\end{array}$$

(29)

The MSE of the coefficients is

$$\begin{array}{cc}& {\text{MSE}}^{\left(M-\mathrm{1}\right)}={\sigma }_{d_{\mathrm{0}}}^{\mathrm{2}}-{\mathit{r}}_{\mathit{x}\mathbf{0},\mathit{d}\mathbf{0}}{\mathit{w}}_{\mathrm{0}}^{\left(M-\mathrm{1}\right)}\end{array}$$

(30)

which is updated with ${\mathit{w}}_{\mathit{d}}^{\left(M-1\right)}$ and the ${\text{MSE}}^{\left(M-1\right)}$ at stage $\left(M-1\right)$ and with ${\mathit{w}}_{\mathit{d}}^{\left(M-2\right)}$ and ${\text{MSE}}^{\left(M-2\right)}$ from the $\left(M-2\right)$ stage.

According to (19) and its property, the tridiagonal covariance matrix can be rewritten as

$$\begin{array}{cc}& {\mathit{R}}_{\mathit{d}}^{\left(M-\mathrm{1}\right)}={\left({\mathit{T}}^{\left(M-\mathrm{1}\right)}\right)}^H{\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{T}}^{\left(M-\mathrm{1}\right)}=\left[\begin{array}{cc}{\mathit{R}}_{\mathrm{1,1}}& {\mathit{R}}_{\mathrm{1,2}}\\ {\mathit{R}}_{\mathrm{2,1}}& r_{M-\mathrm{1},M-\mathrm{1}}\end{array}\right],\end{array}$$

(31)

where

$$\begin{array}{cc}& {\mathit{R}}_{\mathrm{1,1}}={\left({\mathit{T}}^{\left(M-\mathrm{2}\right)}\right)}^H{\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{T}}^{\left(M-\mathrm{2}\right)},\\ {\mathit{R}}_{\mathrm{1,2}}={\left[\begin{array}{cc}{\mathbf{0}}^T& r_{M-\mathrm{2},M-\mathrm{1}}\end{array}\right]}^T,\\ {\mathit{R}}_{\mathrm{2,1}}=\left[\begin{array}{cc}{\mathbf{0}}^T& r_{M-\mathrm{2},M-\mathrm{1}}^{*}\end{array}\right].\end{array}$$

(32)

The cross correlation vector between the new observation vector and desired signal $d_0\left(n\right)$ is

$$\begin{array}{cc}& {\mathit{r}}_{\mathit{d},\mathit{d}\mathbf{0}}^{\left(M-\mathrm{1}\right)}={\left({\mathit{T}}^{\left(M-\mathrm{1}\right)}\right)}^H{\mathit{r}}_{\mathit{x}\mathbf{0},\mathit{d}\mathbf{0}}=\left[\begin{array}{c}{‖{\mathit{r}}_{\mathit{x}\mathbf{0},\mathit{d}\mathbf{0}}‖}_{\mathrm{2}}\\ \mathit{0}\end{array}\right].\end{array}$$

(33)

Given ${\mathit{R}}_{\mathit{d}}^{\left(M-2\right)}$ from stage $\left(M-2\right)$ , the new elements of ${\mathit{R}}_{\mathit{d}}^{\left(M-1\right)}$ are calculated as

$$\begin{array}{cc}& r_{M-\mathrm{1},M-\mathrm{1}}={\mathit{t}}_{M-\mathrm{1}}^H{\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{t}}_{M-\mathrm{1}},\\ r_{M-\mathrm{2},M-\mathrm{1}}={\mathit{t}}_{M-\mathrm{2}}^H{\mathit{R}}_{\mathit{x}\mathrm{0}}{\mathit{t}}_{M-\mathrm{1}}.\end{array}$$

(34)

According to (26), the (34) can be rewritten as

$$\begin{array}{cc}& r_{M-\mathrm{1},M-\mathrm{1}}={\displaystyle \sum }_{n=\mathrm{0}}^{L-\mathrm{1}}{d}_{M-\mathrm{1}}^{*}\left(n\right)d_{M-\mathrm{1}}\left(n\right),\\ r_{M-\mathrm{2},M-\mathrm{1}}={\displaystyle \sum }_{n=\mathrm{0}}^{L-\mathrm{1}}{d}_{M-\mathrm{2}}^{*}\left(n\right)d_{M-\mathrm{1}}\left(n\right).\end{array}$$

(35)

Consider the property that only the first element of the cross correlation vector $r_{d,d0}^{\left(M-1\right)}$ is not equal to 0. Therefore, only the first column of the inverse of ${\mathit{R}}_{\mathit{d}}^{\left(M-1\right)}$ is needed to calculate the recursion coefficients via (28).

Let the inverse of ${\mathit{R}}_{\mathit{d}}^{\left(M-1\right)}$ be noted as

$$\begin{array}{lll}{\mathbf{C}}^{(M-1)}\hfill & =\hfill & {\left({\mathbf{R}}_{\mathbf{d}}^{(M-1)}\right)}^{-1}\hfill \\ \hfill & =\hfill & \left[{\mathbf{c}}_1^{(M-1)},{\mathbf{c}}_2^{(M-1)},\dots ,{\mathbf{c}}_{M-1}^{(M-1)}\right]\hfill \\ \hfill & =\hfill & \left[\begin{array}{cc}{\mathbf{C}}^{(M-2)}& \mathbf{0}\\ {\mathbf{0}}^T& 0\end{array}\right]+{\beta }_{M-1}^{-1}{\mathbf{b}}^{(M-1)}{\left({\mathbf{b}}^{(M-1)}\right)}^H.\hfill \end{array}$$

(36)

The various quantities in (36) are defined as in the following equation:

$$\begin{array}{cc}& {\mathit{b}}^{\left(M-\mathrm{1}\right)}=\left[\begin{array}{c}{r}_{M-\mathrm{2},M-\mathrm{1}}{\mathit{c}}_{M-\mathrm{2}}^{\left(M-\mathrm{2}\right)}\\ \mathrm{1}\end{array}\right],\\ {\beta }_{M-\mathrm{1}}=r_{M-\mathrm{1},M-\mathrm{1}}-{\left|r_{M-\mathrm{1},M-\mathrm{1}}\right|}^{\mathrm{2}}{c}_{M-\mathrm{2},M-\mathrm{2}}^{\left(M-\mathrm{2}\right)},\end{array}$$

(37)

where $c_{M-2,M-2}^{\left(M-2\right)}$ is the last element of the last column ${\mathit{c}}_{M-2}^{\left(M-2\right)}$ from the previous stage.

Therefore, the column vector ${\mathit{c}}_1^{\left(M-1\right)}$ of ${\mathit{C}}^{\left(M-1\right)}$ can be calculated as

$$\begin{array}{cc}& {\mathit{c}}_{\mathrm{1}}^{\left(M-\mathrm{1}\right)}=\left[\begin{array}{c}{\mathit{c}}_{\mathrm{1}}^{\left(M-\mathrm{2}\right)}\\ 0\end{array}\right]+{\beta }_{M-\mathrm{1}}^{-\mathrm{1}}{\left(c_{\mathrm{1},M-\mathrm{2}}^{\left(M-\mathrm{2}\right)}\right)}^{*}\left[\begin{array}{c}{\left|r_{M-\mathrm{2},M-\mathrm{1}}\right|}^{\mathrm{2}}{\mathit{c}}_{M-\mathrm{2}}^{\left(M-\mathrm{2}\right)}\\ -r_{M-\mathrm{2},M-\mathrm{1}}^{*}\end{array}\right],\end{array}$$

(38)

where $c_{1,M-2}^{\left(M-2\right)}$ is the first element of ${\mathit{c}}_{M-2}^{\left(M-2\right)}$ .

It can be seen from (38) that ${\mathit{c}}_1^{\left(M-1\right)}$ depends on ${\mathit{c}}_1^{\left(M-2\right)}$ from stage $\left(M-2\right)$ and the new elements $r_{M-2,M-1}$ and $r_{M-1,M-1}$ generated from the covariance matrix at stage $\left(M-1\right)$ . And the coefficients of the Wiener filter ${\mathit{w}}_{\mathit{d}}^{\left(M-1\right)}$ are also depending on that. Moreover, the last column vector ${\mathit{c}}_{M-1}^{\left(M-1\right)}$ of ${\mathit{C}}^{\left(M-1\right)}$ depends on the last column ${\mathit{c}}_{M-2}^{\left(M-2\right)}$ from stage $\left(M-2\right)$ . According to (36), the last column vector ${\mathit{c}}_{M-1}^{\left(M-1\right)}$ can be updated as

$$\begin{array}{cc}& {\mathit{c}}_{M-\mathrm{1}}^{\left(M-\mathrm{1}\right)}={\beta }_{M-\mathrm{1}}^{-\mathrm{1}}\left[\begin{array}{c}-r_{M-\mathrm{2},M-\mathrm{1}}{\mathit{c}}_{M-\mathrm{2}}^{\left(M-\mathrm{2}\right)}\\ \mathrm{1}\end{array}\right].\end{array}$$

(39)

It can be seen from (39) that the last column vector ${\mathit{c}}_{M-1}^{\left(M-1\right)}$ is only depending on the last column vector ${\mathit{c}}_{M-2}^{\left(M-2\right)}$ from stage $\left(M-2\right)$ and the new element $r_{M-2,M-1}$ generated from the covariance matrix at stage $\left(M-1\right)$ .

According to (38) and (39), it can be seen that, in recursive calculation process, only ${\mathit{c}}_1^{\left(M-1\right)}$ and ${\mathit{c}}_{M-1}^{\left(M-1\right)}$ are needed to be updated at each stage. This avoids the calculation of the inversion of covariance matrix, which also reduces the complexity of computation.

As for the MSE expressed as in (30), it can be simple and can be updated with $c_{\mathrm{1,1}}^{\left(M-1\right)}$ generated from the covariance matrix at stage $\left(M-1\right)$ as follows:

$$\begin{array}{cc}& {\text{MSE}}^{\left(M-\mathrm{1}\right)}={\sigma }_{d_{\mathrm{0}}}^{\mathrm{2}}-{‖{\mathit{r}}_{\mathit{x}\mathbf{0},\mathit{d}\mathbf{0}}‖}_{\mathrm{2}}^{\mathrm{2}}{c}_{\mathrm{1,1}}^{\left(M-\mathrm{1}\right)}.\end{array}$$

(40)

According to recursive algorithm about the calculation of the coefficients of the match filters and the nest order, the data level order recursive MSNWF DOA estimation structure can be drawn as in Figure 4.

Figure 4

The structure of data level order recursive MSNWF.

4. MSNWF DOA Estimation System

4.1. Calculation of Weight Vector

In the MSNWF DOA system, the optimal estimation of the phase-shifted reference signal $e^{j{\varphi }_k}{\mathit{r}}_k$ in the minimum mean square error sense can be obtained at the output of the adaptive beamformer B , which uses the adaptive beamforming weights obtained from the adaptive beamformer A with the MSNWF structure.

In the adaptive beamformer B , consider the case where the phase-shifted reference signal $e^{j{\varphi }_k}{\mathit{r}}_k$ is the desired signal, and the output of the adaptive beamformer B can be used to estimate the desired signal. Since the phase-shifted $e^{j{\varphi }_k}$ is unknown, both the phase-shifted reference signal and the weight vector of the adaptive beamformer B are not available. However, the weight vector of the adaptive beamformer B can be obtained from the optimal weights of the adaptive beamformer A .

In the adaptive beamformer A , the desired signal and observation vector can be given by

$$\begin{array}{cc}& d_{A\mathrm{0}}\left(n\right)=r_k\left(n\right),\hspace{1em}\hspace{1em}{\mathit{x}}_{A\mathrm{0}}\left(n\right)={\mathit{y}}_A\left(n\right).\end{array}$$

(41)

The optimal weight vector of adaptive beamformer A can be readily obtained according to (42).

The flow diagram of calculation of weight vectors in adaptive beamformer A is as follows:

$$\begin{array}{l}{\mathbf{r}}_{{\mathbf{x}}_{A}0,{\mathbf{d}}_{A}0}=E[{\mathbf{x}}_{A0}(n)d_{A0}^{*}(n)]\hfill \\ r_{0,1,B}=0,c_{1,A}^{(1)}=r_{1,1,A}^{-1}\hfill \\ {\text{MSE}}_A^{(1)}={\sigma }_{d_{0,A}}^2-{\Vert {\mathbf{r}}_{{\mathbf{x}}_A\mathbf{0},{\mathbf{d}}_A\mathbf{0}}\Vert }_2^2{c}_{1,A}^{(1)}\hfill \\ \text{FOR}\hspace{0.17em}i=2,3,\dots ,M-1\hfill \\ {\mathbf{t}}_{m,A}=\sum _{n=0}^{L-1}{d}_{m-1,A}^{*}(n){\mathbf{x}}_{m-1,A}(n)\hfill \\ d_{i,A}(n)={\mathbf{t}}_{m,A}^H{\mathbf{x}}_{m-1,A}(n)\hfill \\ {\mathbf{x}}_{m,A}(n)={\mathbf{x}}_{m-1,A}(n)-d_{i,A}(n){\mathbf{t}}_{m,A}\hfill \\ r_{m-1,m,A}=\sum _{n=0}^{L-1}{d}_{m-1,A}^{*}(n)d_{m,A}(n)\hfill \\ r_{m,m,A}=\sum _{n=0}^{L-1}{d}_{m,A}^{*}(n)d_{m,A}(n)\hfill \\ {\beta }_{i,A}=r_{i,i,A}-{\left|r_{i,i,A}\right|}^2{c}_{1,i-1,A}^{(i-1)}\hfill \\ {\mathbf{c}}_{1,A}^{(i)}=\left[\begin{array}{c}{\mathbf{c}}_{1,A}^{(i-1)}\\ 0\end{array}\right]+{\beta }_{i,A}^{-1}{\left(c_{i,1,A}^{(i-1)}\right)}^{*}\left[\begin{array}{c}{\left|r_{i-1,i,A}\right|}^2{\mathbf{c}}_{i,A}^{(i-1)}\\ -r_{i-1,i,A}^{*}\end{array}\right]\hfill \\ {\mathbf{c}}_{i,A}^{(i)}={\beta }_{i,A}^{-1}\left[\begin{array}{c}-r_{i-1,i,A}{\mathbf{c}}_{i,A}^{(i-1)}\\ 1\end{array}\right]\hfill \\ {\text{MSE}}_A^{(i)}={\sigma }_{d_{A0}}^2-{\Vert {\mathbf{r}}_{{\mathbf{x}}_A\mathbf{0},{\mathbf{d}}_A\mathbf{0}}\Vert }_2^2{c}_{1,1,A}^{(i)}\hfill \\ \text{END}\hfill \\ {\mathbf{T}}_A^{(M-1)}=[{\mathbf{t}}_{1,A},{\mathbf{t}}_{2,A},\dots ,{\mathbf{t}}_{M-1,A}]\hfill \\ {\mathbf{w}}_{A0}^{(M-1)}={\mathbf{T}}_A^{(M-1)}{\mathbf{c}}_{1,A}^{(M-1)}.\hfill \end{array}$$

(42)

In the adaptive beamformer B , the phase-shifted desired signal and observation vector can be given by

$$\begin{array}{cc}& d_{B\mathrm{0}}\left(n\right)=e^{j{\varphi }_k}{r}_k.\end{array}$$

(43)

And the optimal weight vector of adaptive beamformer B can be obtained according to (42) as shown in (44).

The flow diagram of the calculation of weight vector in adaptive beamformer B is as follows:

$$\begin{array}{l}{\mathbf{r}}_{{\mathbf{x}}_{B}0,d_{B}0}=E[{\mathbf{x}}_{B0}(n)d_{B0}^{*}(n)]\hfill \\ r_{0,1,B}=0,c_{1,B}^{(1)}=r_{1,1,B}^{-1}\hfill \\ {\text{MSE}}_B^{(1)}={\sigma }_{d_{0,B}}^2-{\Vert {\mathbf{r}}_{{\mathbf{x}}_B\mathbf{0},{\mathbf{d}}_B\mathbf{0}}\Vert }_2^2{c}_{1,B}^{(1)}\hfill \\ \text{FOR}\hspace{0.17em}i=2,3,\dots ,M-1\hfill \\ {\mathbf{t}}_{m,B}=\sum _{n=0}^{L-1}{d}_{m-1,B}^{*}(n){\mathbf{x}}_{m-1,B}(n)\hfill \\ d_{i,B}(n)={\mathbf{t}}_{m,B}^H{\mathbf{x}}_{m-1,B}(n)\hfill \\ {\mathbf{x}}_{m,B}(n)={\mathbf{x}}_{m-1,B}(n)-d_{i,B}(n){\mathbf{t}}_{m,B}\hfill \\ r_{m-1,m,B}=\sum _{n=0}^{L-1}{d}_{m-1,B}^{*}(n)d_{m,B}(n)\hfill \\ r_{m,m,B}=\sum _{n=0}^{L-1}{d}_{m,.B}^{*}(n)d_{m,B}(n)\hfill \\ r_{m,m,B}=\sum _{n=0}^{L-1}{d}_{m,B}^{*}(n)d_{m,B}(n)\hfill \\ {\beta }_{i,B}=r_{i,i,B}-{\left|r_{i,i,B}\right|}^2{c}_{1,i-1,B}^{(i-1)}\hfill \\ {\mathbf{c}}_{1,B}^{(i)}=\left[\begin{array}{c}{\mathbf{c}}_{1,B}^{(i-1)}\\ 0\end{array}\right]+{\beta }_{i,B}^{-1}{\left(c_{i,1,B}^{(i-1)}\right)}^{*}\left[\begin{array}{c}{\left|r_{i-1,i,B}\right|}^2{\mathbf{c}}_{i,B}^{(i-1)}\\ -r_{i-1,i,B}^{*}\end{array}\right]\hfill \\ {\mathbf{c}}_{i,B}^{(i)}={\beta }_{i,B}^{-1}\left[\begin{array}{c}-r_{i-1,i,B}{\mathbf{c}}_{i,B}^{(i-1)}\\ 1\end{array}\right]\hfill \\ {\text{MSE}}_B^{(i)}={\sigma }_{d_{B0}}^2-{\Vert {\mathbf{r}}_{{\mathbf{x}}_B\mathbf{0},{\mathbf{d}}_B\mathbf{0}}\Vert }_2^2{c}_{1,1,B}^{(i)}\hfill \\ \text{END}\hfill \\ {\mathbf{T}}_B^{(M-1)}=[{\mathbf{t}}_{1,B},{\mathbf{t}}_{2,B},\cdots ,{\mathbf{t}}_{M-1,B}]\hfill \\ {\mathbf{w}}_{B0}^{(M-1)}={\mathbf{T}}_B^{(M-1)}{\mathbf{c}}_{1,B}^{(M-1)}.\hfill \end{array}$$

(44)

Substituting (11) and (13) into (44), we have

$$\begin{array}{cc}& {\mathit{w}}_{A\mathrm{0}}^{\left(M-\mathrm{1}\right)}={\mathit{w}}_{B\mathrm{0}}^{\left(M-\mathrm{1}\right)}.\end{array}$$

(45)

Therefore, the weight vector ${\mathit{w}}_{B0}^{\left(M-1\right)}$ can be obtained by calculating the optimal weight of the adaptive beamformer A .

4.2. Calculation of DOA

The adaptive beamformer B based on the structure of data level order recursive MSNWF can be simplified to a single stage Wiener filter in virtue of obtaining its weight from the adaptive beamformer A . Let ${\widehat{\mathit{r}}}_k\left(n\right)={({\mathit{w}}_{B0}^{\left(M-1\right)})}^H{\mathit{y}}_B$ denote the output signal of beamformer B . Let

$$\begin{array}{cc}& {\widehat{\mathit{r}}}_k={\left[{\widehat{\mathit{r}}}_k\left(\mathrm{1}\right),{\widehat{\mathit{r}}}_k\left(\mathrm{2}\right),\dots ,{\widehat{\mathit{r}}}_k\left(L\right)\right]}^T.\end{array}$$

(46)

Thus, ${\widehat{\mathit{r}}}_k$ is an optimal estimation of the phase-shifted reference signal   $e^{j{\varphi }_k}{\mathit{r}}_k$ in the MMSE sense, which can be written as

$$\begin{array}{cc}& {\widehat{\mathit{r}}}_k=e^{j{\varphi }_k}{\mathit{r}}_k+{\mathit{N}}_k.\end{array}$$

(47)

Let ${\widehat{\varphi }}_k$ denote an estimation of ${\varphi }_k$ , which can be calculated by using the least square method such that the square error between the two signal vectors ${\widehat{\mathit{r}}}_k$ and ${\mathit{r}}_k$ is minimized

$$\begin{array}{cc}& \underset{\min }{{\widehat{\varphi }}_k}{‖{\widehat{\mathit{r}}}_k-e^{j{\widehat{\varphi }}_k}{\mathit{r}}_k‖}_{\mathrm{2}}.\end{array}$$

(48)

In [18], Wang et al. give the optimum solution of ${\widehat{\varphi }}_k$

$$\begin{array}{cc}& {\widehat{\varphi }}_k\hspace{0.17em}\hspace{0.17em}=\arg \left({\widehat{\mathit{r}}}_k{\mathit{r}}_k^H\right).\end{array}$$

(49)

According to (12), an estimation of the target DOA can be obtained then as

$$\begin{array}{cc}& {\widehat{\theta }}_k\hspace{0.17em}\hspace{0.17em}=\arcsin \left(\frac{-\lambda {\widehat{\varphi }}_k}{\mathrm{2}\pi d}\right).\end{array}$$

(50)

5. Simulation Results

In this section, the performance of the proposed method, including the resolution, capacity, and accuracy of the data level order recursive MSNWF DOA techniques, will be evaluated through numerical simulations. In Sections 5.1 and 5.2, the resolution and the capacity of the DOA estimation using the data level order recursive MSNWF DOA techniques will be illustrated and compared with other techniques, such as MUSIC, ESPRIT, SBDOA, and original MSNWF DOA estimation techniques. In Sections 5.3 and 5.4, the effects of snapshot length and stage of data level order recursive MSNWF on the estimation accuracy will be investigated, respectively.

5.1. Resolution of DOA Estimation

Assume that a ULA of 10 elements, with a spacing of $d=\lambda /2$ deployed at the receiver, was employed in the simulations, to deal with a case where the DOAs of three signals and two interference signals are closely distributed. Further, assume that the DOAs of the target related signal components are at −2°, 0°, and 2°. The DOAs of the interference components are at −4° and 4°. The background noise power spectral density ratio of the received signal is set to 10 dB. Snapshot length is fixed at 100 and the stage of MSNWF is set to 5. One thousand simulation runs were performed. These simulation results are illustrated in Figure 5.

Figure 5

Comparison of the resolution of DOA estimation for signal sources that are closely distributed.

The histograms of the resolution of DOA estimation obtained for these five techniques are shown in Figures 5(a)–5(e). The histogram depicts the number of occurrences estimated DOA as a function of DOA degrees. In Figure 5(a), the histogram of MUSIC technique shows two peak values which deviate from the DOAs of the target signals. In Figure 5(b), although the histogram of ESPIRT technique shows three peak values, the peak values deviate from the DOAs of the target signals. It is seen that the MUSIC technique or ESPRIT technique cannot offer the desired results when the DOAs of target signals are very close. Correspondingly, in Figures 5(c), 5(d), and 5(e), the histogram shows three peak values, indicating that, using the SBDOA, original MSNWF DOA, and the data level order recursive MSNWF DOA techniques, all three DOAs are successfully estimated. Therefore, it proved that the data level order recursive MSNWF DOA technique could obtain a better resolution than MUSIC and ESPRIT techniques. However, the SBDOA requires $O(M^3)$ operations, the original MSNWF DOA technique needs $O(2M^2+9M)$ operations, and the data level order recursive MSNWF DOA technique demands $O(M^2+11M)$ operations, which significantly reduce the complexity of computation. In addition, if the recursive order is enough, the resolution of data level order recursive MSNWF DOA technique will be well as that of SBDOA estimation and it is proved in Section 5.4. The resolution and accuracy of data level order recursive MSNWF are better than the original MSNWF which is due to the update of the MSE at each stage.

5.2. Capacity of DOA Estimation

This simulation deals with a case where the number of target signals and interference is larger than that of antenna elements. The simulation conditions are kept the same as those in Section 5.1 except for the number of signal components considered. The DOAs of 9 target signal components are set from −40° to 40° with interval 10°, and the DOAs of 6 interference components are set from −25° to 25° with interval 10°. The simulation results are shown in Figure 6.

Figure 6

Comparison of the capacity of DOA estimation when the number of signal and interference sources exceeds the number of antenna elements.

Histograms of the obtained estimated DOAs are shown in Figures 6(a)–6(e). In Figures 6(a) and 6(b), the histograms show the deviated peak values and demonstrate that these two techniques cannot provide acceptable DOA estimation when the number of antenna elements is less than the total number of target signals and interference. In contrast, in Figures 6(c), 6(d), and 6(e), the histograms show that all 9 target DOAs are successfully estimated when using the SBDOA, original MSNWF, and data level order recursive MASNWF DOA techniques. As can be seen, the successful probability of DOA estimation in data level order recursive MSNWF DOA technique is the same as that in the SBDOA and original MSNWF DOA estimation techniques.

5.3. Effects of Snapshot Length of MSNWF on DOA Estimation Accuracy

In the simulation of snapshot length effects, the snapshot length for adaptive beamformer A and DOA calculation are set to different values, such as 50, 100, 200, 500, and 1000, and the stages of both original MSNWF and data level order recursive MSNWF are set to 5. The DOA of the target signal is fixed at 0°, and the DOAs of the interference are set from −90° to 90° with interval 10° except 0°. The root mean square error (RMSE) of the estimated target DOA averaged over one thousand simulation runs at different SNR conditions; the RMSE of the estimated target DOA and the snapshot length are illustrated in Figure 7.

Figure 7

RMSE of the estimated DOA for different snapshot length L and the SNR.

As can be seen in Figure 7, both the original MSNWF and the data level order recursive MSNWF DOA techniques lead to a RMSE of less than 5°,when using a small snapshot length such as 50. The simulation results show that when the snapshot length is 500, the data level order recursive MSNWF DOA estimation method will have estimation accuracy similar to that of the SBDOA technique. However, the RMSE of the data level order recursive MSNWF DOA technique is better than that of original MSNWF DOA technique under various snapshot lengths, which is due to the update of MSE at each stage to obtain the optimal weight vector. The RMSE obviously decreases as the snapshot length increases, such as the RMSE which will be less than 1° when using one thousand snapshot length of the signal. This demonstrates that the fast DOA tracking can be implemented by using the data level order recursive MSNWF DOA technique and that the estimation accuracy will be improved when using more sample data. And the simulation also proved that the capacity of data level order recursive MSNWF DOA estimation technique can be larger than the number of the sensor elements.

5.4. Effects of the Stage of MSNWF on DOA Estimation Accuracy

In the simulation of the stage effect, both stages of original MSNWF and data level order recursive MSNWF for adaptive beamformer A are set to the same values, such as 3, 5, and 9. The snapshot length is set to 200. And other simulation conditions are kept the same as those in Section 5.3. The RMSE of the estimated target DOA averaged over one thousand simulations runs at different SNR conditions. The RMSE of the estimated target DOA with different stages of MSNWF and SNR is demonstrated in Figure 8.

Figure 8

RMSE of the estimated DOA for different MSNWF stages and SNR.

As can be seen from Figure 8, both the original MSNWF and the data level order recursive MSNWF DOA techniques lead to a RMSE of less than 3°, when using different stages of MSNWF, and the RMSE decreases as the MSNWF stage increases. However, the RMSE of the data level order recursive MSNWF DOA estimation technique is better than that of original MSNWF DOA estimation technique under various stages, which is mainly due to the update of MSE at each stage.

Moreover, in the same simulation conditions, the RMSE of SBDOA estimation technique is less than 1.5°. In contrast, the RMSE of data level order recursive MSNWF DOA estimation technique is almost equal to that of SBDOA when using 9 stages. However, the original MSNWF DOA estimation technique requires more stages to obtain similar estimation accuracy.

6. Conclusion

A novel DOA estimation method based on data level order recursive MSNWF has been proposed in this paper. In this technique, two subarray adaptive beamformers based on the MSNWF are used to form the phase-shift and reject interference at the same time. The DOAs of target signals are estimated from the phase-shift by using reference signal after interference rejection. Therefore, the performance of DOA estimation such as resolution, capacity, and accuracy is significantly improved. And the complexity of computation is also significantly reduced by avoiding the calculation of covariance matrix inversion when getting the optimal weight vector of the beamformer. This technique can be widely used for the implementation of hardware systems such as wireless communication system, active radar, sonar, and STAP systems. Numerical simulations demonstrating the effectiveness and advantage of this technique are presented.