Direction-of-Arrival Estimation of Virtual Array Signals Based on Doppler Effect

Abstract

The estimation accuracy of direction-of-departure (DOD) and direction-of-arrival (DOA) is reduced because of Doppler shifts caused by the high-speed moving sources. In this paper, an improved DOA estimation method which combines the forward-backward spatial smoothing (FBSS) technique with the MUSIC algorithm is proposed for virtual MIMO array signals in high mobility scenarios. Theoretical analysis and experiment results demonstrate that the resolution capability can be significantly improved by using the proposed method compared to the MUSIC algorithm for the moving sources with limited array elements, especially the DOA which can still be accurately estimated when the sources are much closely spaced.

1. Introduction

In the last several decades, the high resolution direction-of-arrival (DOA) estimation methods using antenna arrays have played an important role in various fields, including mobile communications, radar, sonar, and seismology (a few examples of the many possible applications). For example, public safety communications systems can also benefit from DOA methods for search and rescue operations based on estimated DOAs from victim mobile devices [1]. There are less reflector and direct path dominating in high mobility scenarios. Although Doppler diffusion is not outstanding, but the Doppler shift is more serious, which could affect the symbol detection for the wireless digital communication system, and even leads to significant degradation of the communication quality. The system performance can be improved by these methods that include the augment of the transmission power and signal-to-noise ratio or the decrease of communication link distance, but these methods are limited in practical applications. Nevertheless, it is a grand challenge when we take the extension of antenna aperture on the small aircraft into consideration to improve the resolution ratio of the system. Therefore, it is necessary to improve the resolution of multiple sources by combining the virtual method of the extended array elements with existing MIMO technology in high mobility scenarios.

The virtual array was formed to extend the equivalent array aperture by using the conjugate counterpart of the array output, so that it can handle more sources than the original arrays [2]. Antenna arrays can obtain additional freedom by exploiting virtual antenna arrays, which led to narrower beam width and lower side lobes [3, 4]. However, they did not consider the correlation between the signals, so some scholars consider the decorrelation of received signals firstly and then estimate the DOA with estimation algorithm in hand to improve resolution ratio [5–7]. In addition, the DOD was estimated by using an ESPRIT algorithm with two-dimensional direction search method and the DOA was estimated by using a ROOT-MUSIC algorithm with one-dimensional search method are not only complicated to calculate, but need extra angle matching [8]. But it has been shown that the DOD and DOA could be estimated by using one-dimensional search method in [9], which avoid nonlinear search of 2D spectrum peak and iteration calculation; thus the method reduced computational complexity greatly without losing accuracy and at the same time parameters can be paired automatically. And DOA was estimated by using a modified MUSIC algorithm based on virtual array transformation and spatial smoothing technique [10–12], which avoided seeking the peak of the spectrum of traditional MUSIC. Moreover, the new DOA estimation algorithm [13] was proposed, which made full use of weighted spatial smoothing combining the autocorrelation information with cross-correlation information of the submatrix. At the same time, the DOA estimation algorithm of coherent sources was proposed by using the ESPRIT algorithm and the MUSIC algorithm with multiple and invariant features [14, 15]. However, these literature have only done virtual extension for receiver arrays. Few of them considered the problem of combining virtual MIMO arrays with the MUSIC algorithm (using the FBSS technique) in high mobility scenarios. So this work will study the DOD-DOA estimation which can be applied to angle estimation of high-mobility sources by combining the virtual array transformation with the MUSIC algorithm (using FBSS processing), which can choose flexibly the number of virtual array elements and their locations according to the actual environment. Moreover, our proposed approach can improve localization precision effectively in case of the sources with similar angle and limited array elements.

This paper has the following five parts. The system model and the mathematical expression for the virtual array signals are provided in Section 2. The proposed MUSIC algorithm by using FBSS processing with a detailed analysis of its decorrelation effect is investigated in Section 3. Simulation results are presented in Section 4 and conclusions are drawn in Section 5.

2. Virtual Array Model

Array dates or information of virtual location was constructed to achieve the goal of array extension based on the received signals of actual arrays, which is called virtual arrays extension technology [2]. But the number of antenna elements is limited by apertures of platform in high mobility scenarios, so resolution is improved by using the virtual array extension technology.

2.1. Virtual Aperture Extension

Assume that two antennas are installed in the plane with the limited volume, black antenna $s_0(t)$ and $s_1(t)$ express real-time position of plane and white antenna. $s_{-1}(t)$ is the virtual antenna at one moment, at the same time $s_{-1}(t)$ and $s_1(t)$ are the center symmetry about $s_0(t)$ . The transmitter arrays are composed of two uniform linear array elements (ULA), and receiver comprised N ones as shown in Figure 1. Element spacing of transmitter and receiver are $d_t$ and $d_r$ , respectively. Virtual MIMO system model of transmitter and receiver is shown in Figure 1.

Figure 1

Virtual uniform linear array model.

Suppose that there are K sources close to the receiver at a constant V rate in the air. The DOD and DOA of the kth source are represented as ${\varphi }_k$ , ${\theta }_k$ , respectively. The mutually orthogonal narrow-band coding signals are emitted for the transmitted signals which are represented as $\{s_0(t),\dots ,s_M(t)\}$ .

Suppose that virtual array are symmetrical extrapolation values about the position of real elements; namely, the transmitted signals are represented as $s_{-1}(t)=s_1^{\ast }(t)$ . The complex conjugate output of the mth element can be written as

$$\begin{array}{l}{x}_m^{\ast }\left(t\right)={\displaystyle \sum }_{m=\mathrm{1}}^M{s}_m^{\ast }\left(t\right)e^{-j(m-\mathrm{1})\mathrm{2}\pi d_t/\lambda \mathrm{s}\mathrm{i}\mathrm{n}\left({\phi }_m\right)}+n_m^{\ast }\left(t\right)\\ ={\displaystyle \sum }_{m=\mathrm{1}}^M{s}_m\left(t\right)e^{-j(m-\mathrm{1})\mathrm{2}\pi d_t/\lambda \mathrm{s}\mathrm{i}\mathrm{n}\left({\phi }_m\right)}+n_m^{\ast }\left(t\right).\end{array}$$

(1)

Obviously, $x_m^{\ast }(t)$ can be considered as the output signals of the virtual elements. Therefore, by concatenating the actual array output signal and its conjugate components, the new array output vector [2] can be written as

$$\begin{array}{l}\mathit{y}\left(\mathit{t}\right)=\left[\begin{array}{l}\mathit{T}{\mathit{x}}^{\ast }\left(\mathit{t}\right)\\ \mathit{x}\left(\mathit{t}\right)\end{array}\right]=\left(\begin{array}{l}\mathit{T}{\mathit{A}}^{\ast }\\ \mathit{A}\end{array}\right)\mathit{s}\left(\mathit{t}\right)+\left(\begin{array}{l}\mathit{T}{\mathit{n}}^{\ast }\left(\mathit{t}\right)\\ \mathit{n}\left(\mathit{t}\right)\end{array}\right)\\ =\tilde{\mathit{A}}\mathit{s}\left(\mathit{t}\right)+\mathit{v}\left(\mathit{t}\right),\end{array}$$

(2) where $\mathit{T}={\mathit{I}}_{\mathit{M}}\hspace{0.17em}\hspace{0.17em}(M-\mathrm{1}:\mathrm{1},:)$ and   $\tilde{\mathit{A}}=\left[\mathit{a}({\varphi }_{\mathrm{1}}),\mathit{a}({\varphi }_{\mathrm{2}}),\mathit{a}({\varphi }_M)\right]$ is the direction matrix of the virtual array, whose mth column denoted as $\mathit{a}({\varphi }_m)={\left[e^{-j(m-\mathrm{1})\mathrm{2}\pi d_t/\lambda \mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_m)},\dots ,\mathrm{1},\dots ,e^{j(m-\mathrm{1})\mathrm{2}\pi d_t/\lambda \mathrm{s}\mathrm{i}\mathrm{n}({\varphi }_m)}\right]}^T\in \left(\mathrm{1}\le m\le M\right)$ is the corresponding extended array steering vector.

2.2. Received Signal with Doppler Frequency Shift

The Doppler frequency is less than zero when the sources are far away from the receiver. In other words, the frequency of the received signal is less than the frequency of the transmitted signal, and vice versa. Consider a single transmitted source firstly. Delay of the actual transmitted signals arriving at a certain receiver can be written as

$$\begin{array}{cc}& \tau =\frac{d_t\mathrm{s}\mathrm{i}\mathrm{n}\theta }{c}=\frac{vt\mathrm{s}\mathrm{i}\mathrm{n}\theta }{c}.\end{array}$$

(3)

Phase difference of transmitted arrays arriving at receiver can be obtained as $\phi (t)={\omega }_{\mathrm{0}}\tau =\mathrm{2}\pi \left(vt\mathrm{s}\mathrm{i}\mathrm{n}\theta /\lambda \right)$ which is a function of t. When the radial velocity is constant, Doppler is caused with the frequency difference $\phi (t)$ , which can be expressed as

$$\begin{array}{cc}& f_d=\frac{\mathrm{1}}{\mathrm{2}\pi }\frac{d\phi }{dt}=\frac{v}{\lambda }\mathrm{s}\mathrm{i}\mathrm{n}\theta =\frac{v{f}_{\mathrm{0}}}{c}\mathrm{s}\mathrm{i}\mathrm{n}\theta .\end{array}$$

(4) It is straightforward to obtain expression of N received signals as

$$\begin{array}{l}y\left(t\right)={\left[y_{\hspace{0.17em}1}\left(t\right),y_{\hspace{0.17em}2}\left(t\right),\dots ,y_{N\hspace{0.17em}}\left(t\right)\right]}^T\\ =\alpha \mathrm{e}\mathrm{x}\mathrm{p}\left(j2\pi f_{d}t\right){\mathit{b}}_r\left(\theta \right)\times {\mathit{a}}_t^{\mathit{T}}\left(\varphi \right)\times \mathit{s}\left(t\right)+\mathit{n}\left(t\right).\end{array}$$

(5) The array direction vector [8] in the transmitter can be written as

$$\begin{array}{l}{\mathit{a}}_t\left(\varphi \right)=\left[e^{j\mathrm{2}\pi (M-\mathrm{1})(d_t/\lambda )\mathrm{s}\mathrm{i}\mathrm{n}\varphi },\dots ,\mathrm{1},e^{-j\mathrm{2}\pi (d_t/\lambda )\mathrm{s}\mathrm{i}\mathrm{n}\varphi },\right.\\ {\left.\hspace{1em}\hspace{1em}L{e}^{-j\mathrm{2}\pi (M-\mathrm{1})(d_t/\lambda )\mathrm{s}\mathrm{i}\mathrm{n}\varphi }\right]}^T.\end{array}$$

(6) The array direction vector in the receiver can be written as

$$\begin{array}{cc}& {\mathit{b}}_r\left(\theta \right)={\left[\mathrm{1},e^{-j\mathrm{2}\pi (d_r/\lambda )\mathrm{s}\mathrm{i}\mathrm{n}\theta },L{e}^{-j\mathrm{2}\pi (N-\mathrm{1})(d_r/\lambda )\mathrm{s}\mathrm{i}\mathrm{n}\theta }\right]}^T,\end{array}$$

(7) where $\mathit{n}(t)=[{\mathit{n}}_1(t),{\mathit{n}}_2(t),\dots ,{\mathit{n}}_N(t){]}^T\in \left(N\times L\right)$ is the noise matrix, whose each component $\left(N\times \mathrm{1}\right)$ is assumed to be independent, zero mean, complex, and Gaussian distributed. α is the amplitude attenuation factor with the transmission path, ${\mathit{a}}_t({\varphi }_i)$ is the path response of array from the ${\varphi }_i$ direction, and ${\mathit{b}}_r({\theta }_i)$ is the path response of array from the ${\theta }_i$ direction. In the case of multiple targets, the angles of the transmitting arrays and receiving arrays are $({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_K)$ and $({\theta }_1,{\theta }_2,\dots ,{\theta }_K)$ , respectively, $f_{dk}$ is the Doppler shift of the kth target, and received signals at the time t can be expressed as

$$\begin{array}{cc}& \mathit{y}\left(t\right)={\displaystyle \sum }_{k=\mathrm{1}}^K{\alpha }_k\mathrm{e}\mathrm{x}\mathrm{p}\left(j\mathrm{2}\pi f_{dk}t\right){\mathit{b}}_r\left({\theta }_k\right)·{\mathit{a}}_t^T\left({\phi }_k\right)·\mathit{s}\left(t\right)+\mathit{n}\left(t\right).\end{array}$$

(8) And then (8) is discretized and expression of received arrays [16] can be denoted as

$$\begin{array}{l}\mathit{y}={\displaystyle \sum }_{k=1}^K{\alpha }_k\mathrm{e}\mathrm{x}\mathrm{p}\left(j2\pi f_{dk}l{T}_s\right){\mathit{b}}_r\left({\theta }_k\right)·{\mathit{a}}_t^T\left({\phi }_k\right)\mathit{s}+\mathit{n}\left(l{T}_s\right)\\ =\mathit{B}\mathit{\Lambda }{\mathit{A}}^T\mathit{S}+\mathit{N},\end{array}$$

(9) where the size of receive signal y is $N\times L$ and $l\in (1,L)$ is the number of snapshots within a single pulse. Transmitted signal $S=[s_1,s_2,\dots ,s_m{]}^T$ is a $M\times L$ matrix. Defining ${\mathit{a}}_t({\varphi }_k)$ and ${\mathit{b}}_r({\theta }_k)$ as the steering vectors and the $M\times K$ and $N\times K$ size are represented as transmitted arrays manifold matrix $\mathit{A}=[{\mathit{a}}_t({\varphi }_1),{\mathit{a}}_t({\varphi }_2),\dots ,{\mathit{a}}_t({\varphi }_k)]$ and received arrays manifold matrix $\mathit{B}=[{\mathit{b}}_r({\theta }_1),{\mathit{b}}_r({\theta }_2),\dots ,{\mathit{b}}_r({\theta }_k)]$ , respectively. $T_S$ is the duration time of each symbol, and N is a $N\times L$ noise matrix. $\Lambda =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(c)$ is a $K\times K$ diagonal matrix, where $\mathit{c}=[{\alpha }_1\mathrm{e}\mathrm{x}\mathrm{p}(j2\pi f_{d1}l{T}_S),{\alpha }_2\mathrm{e}\mathrm{x}\mathrm{p}(j2\pi f_{d2}l{T}_S),\dots ,{\alpha }_k\mathrm{e}\mathrm{x}\mathrm{p}(j2\pi f_{dk}l{T}_S)]$ .

3. Proposed Algorithm

Received signals are superposition of K sources, which arrive in N received antennas. Assume that there are $K\times N$ matched filter in the receiver. Then, each signal is separately isolated from the received signal by adopting modern digital signal processing technology such as mixing, filtering, and sampling. Because the signal satisfies orthogonality, there are

$$\begin{array}{cc}& \frac{\mathrm{1}}{L}{\displaystyle \sum }_{l=\mathrm{1}}^L{S}_i\left(l\right)S_j^{\ast }\left(l\right)=\left\{\begin{array}{cc}\mathrm{1},\hfill & i=j\hfill \\ \mathrm{0},\hfill & i\ne j.\hfill \end{array}\right.\end{array}$$

(10) Then the $K\times N$ independent output signal ${\mathit{y}}_{\mathit{k},\mathit{n}}$ is obtained by multiplying the separated received signals y by the transmitted signal $S_1^H,S_2^H,\dots ,S_K^H$ . The matched filtering structure [16] is shown in Figure 2.

Figure 2

Matched filtering of the receiver array.

The output of l sample point after the relevant processing of receiver at a time is written as

$$\begin{array}{cc}& \mathbf{Y}=\mathbf{y}{\mathbf{S}}^H=\mathbf{B}\mathit{\Lambda }{\mathbf{A}}^T\mathbf{S}{\mathbf{S}}^H+\mathbf{W}{\mathbf{S}}^H=\mathbf{B}\mathit{\Lambda }{\mathbf{A}}^T+\mathbf{N}\end{array}$$

(11) according to the column vector; the $N\times M$ matrix Y can be turned into $MN\times L$ dimensional vector $\widehat{\mathit{Y}}$

$$\begin{array}{cc}& \widehat{\mathit{Y}}=\left(\mathit{A}\oplus \mathit{B}\right){\mathit{c}}_k+\widehat{\mathit{N}},\end{array}$$

(12) where $\widehat{\mathbf{Y}}=\mathrm{v}\mathrm{e}\mathrm{c}(\mathbf{Y})$ , $\widehat{\mathbf{N}}=\mathrm{v}\mathrm{e}\mathrm{c}(\mathbf{N})$ , the ⊗ is Kronecker product, $\mathrm{v}\mathrm{e}\mathrm{c}(·)$ is matrix which was written as a column vector according to the principle $\widehat{\mathbf{N}}=\mathrm{v}\mathrm{e}\mathrm{c}(\mathbf{N})$ from left to right, and $\mathbf{A}\oplus \mathbf{B}=[{\mathbf{a}}_t({\phi }_{\mathrm{1}})\otimes {\mathbf{b}}_r({\theta }_{\mathrm{1}}),{\mathbf{a}}_t({\phi }_{\mathrm{2}})\otimes {\mathbf{b}}_r({\theta }_{\mathrm{2}}),\dots ,{\mathbf{a}}_t({\phi }_k)\otimes {\mathbf{b}}_r({\theta }_k)]$ . So the total output signal matrix after matched filter is represented as

$$\begin{array}{cc}& \tilde{\mathbf{Y}}={\displaystyle \sum }_{k=\mathrm{1}}^K(\mathbf{A}\oplus \mathbf{B})\mathbf{C}+\tilde{\mathbf{N}},\end{array}$$

(13) where the $C=[{\sum }_{k=1}^k{c}_k]$ is a $1\times L$ matrix, the $\widehat{\mathbf{N}}$ is a $MN\times L$ matrix, the $\mathit{A}\oplus \mathit{B}$ is $MN\times 1$ matrix, and C contains the target scattering coefficient and Doppler shifts information.

First received signals are divided into two subspaces [11, 17]; namely, they are divided into array element $x_0$ with carried information ${\tilde{\mathbf{Y}}}_{\mathbf{1}}=\tilde{\mathbf{Y}}(1:N_r,:)$ and array element $x_1$ with carried information ${\tilde{\mathit{Y}}}_2=\tilde{\mathit{Y}}(N_r:2N_r,:)$ . Then, the covariance matrix of first submatrix ${\tilde{Y}}_{\mathit{1}}$ is expressed as

$$\begin{array}{cc}& {\mathbf{R}}_{\mathrm{11}}=E\left[{\tilde{\mathbf{Y}}}_{\mathrm{1}}{\tilde{\mathbf{Y}}}_{\mathrm{1}}^H\right].\end{array}$$

(14) Mutual covariance of two submatrix can be given as

$$\begin{array}{cc}& {\mathit{R}}_{21}=\frac{{\tilde{\mathit{Y}}}_2{\tilde{\mathit{Y}}}_1}{L}.\end{array}$$

(15) The singular value decomposition (SVD) [18] is performed for the autocorrelation matrix ${\mathit{R}}_{11}$ , namely, $[\mathit{U},\mathit{S},\mathit{V}]=\mathrm{S}\mathrm{V}\mathrm{D}({\mathit{R}}_{11})$ . The reconstructed covariance matrix by the SVD of the ${\mathbf{R}}_{\mathrm{11}}$ is denoted as

$$\begin{array}{cc}& {\mathbf{R}}_{\mathrm{11}}^{\prime }=\mathbf{V}{\mathbf{S}}^{-\mathrm{1}}{\mathbf{U}}^H.\end{array}$$

(16) Mixed covariance by extracting the signal subspace are represented as

$$\begin{array}{cc}& \mathbf{R}={\mathbf{R}}_{\mathrm{21}}{\mathbf{R}}_{\mathrm{11}}^{\prime }.\end{array}$$

(17) Mixed covariance R of expression (17) by using eigenvalue decomposition is expressed as

$$\begin{array}{cc}& \left[\mathit{\vartheta },\mathbf{D}\right]=\mathrm{e}\mathrm{i}\mathrm{g}\left(\mathbf{R}\right),\end{array}$$

(18) where $\mathbf{D}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[e^{-j2\pi d_t{f}_0\mathrm{s}\mathrm{i}\mathrm{n}\varphi /c},0,\dots ,0{]}_{N_r\times N_r}$ . Thus, the DOD [7–9] can be given as

$$\begin{array}{cc}& {\widehat{\phi }}_k=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{c}{\mathrm{2}\pi d{f}_{\mathrm{0}}}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}\left(\mathbf{D}\right)\right)\hspace{1em}k=\mathrm{1,2},\dots ,K.\end{array}$$

(19) Signal subspace ${\vartheta }_s=\vartheta (:,\mathrm{1})$ is extracted from feature vectors ϑ which is decomposed by R and subspace is processed as follows:

$$\begin{array}{cc}& {\widehat{\mathit{\vartheta }}}_s=\frac{{\mathit{\vartheta }}_s(\mathrm{2}:N_r,:)}{{\mathit{\vartheta }}_s(\mathrm{1}:N_r-\mathrm{1},:)}.\end{array}$$

(20) Therefore, the DOA can be expressed as

$$\begin{array}{cc}& {\theta }_k=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{c}{\mathrm{2}\pi d{f}_{\mathrm{0}}}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}\left({\widehat{\mathit{\vartheta }}}_s\right)\right)\hspace{1em}k=\mathrm{1,2},\dots ,K.\end{array}$$

(21) The DOA is estimated by further precise processing with the formula (21)

$$\begin{array}{cc}& {\widehat{\theta }}_k=\frac{\mathrm{1}}{N_r-\mathrm{1}}{\displaystyle \sum }_{k=\mathrm{1}}^{N_r-\mathrm{1}}{\theta }_k.\end{array}$$

(22)

The receiving signals of arrays are coherent through the above algorithm analysis, for which the dimensions of the signal subspace can only reflect direction of arrival of the irrelevant signals. In other words, corresponding characteristic vector of the coherent signals transforms into the noise subspace, resulting in performance being deteriorated. So some direction vectors of coherent signals and noise subspaces are not orthogonal, for which we cannot estimate the DOA of coherent signal by using high resolution feature decomposition algorithm [7, 13]. In the preprocessing of the coherent signals, it is necessary to remove the correlation by using a spatial smoothing [12–14] with reducing dimension method between signals. The modified algorithm can be described as follows.

We consider that the received signals can remove the correlation. Forward-backward spatial smoothing (FBSS) is used to preprocess the data of the received arrays. After the received signal $\widehat{\mathit{Y}}$ performs subspace classification [17], then we divide it into the information ${\widehat{\mathit{Y}}}_1=\widehat{\mathit{Y}}(1:N_r,:)$ carried by $x_0$ and the information ${\widehat{\mathit{Y}}}_2=\widehat{\mathit{Y}}(N_r:2N_r,:)$ carried by $x_1$ , respectively. So covariance matrix of forward subarray of array elements $x_0$ in which the number of snapshots is l can be obtained as

$$\begin{array}{cc}& {\mathbf{R}}_{x^{\prime }}=E\left[{\widehat{\mathbf{Y}}}_{\mathrm{l}}{{\widehat{\mathbf{Y}}}_{\mathrm{l}}}^H\right].\end{array}$$

(23) Define the covariance matrix of the forward spatial smoothing as

$$\begin{array}{cc}& {\mathbf{R}}_x=\frac{\mathrm{1}}{L}{\displaystyle \sum }_{l=\mathrm{1}}^L{\mathbf{R}}_{x^{\prime }}.\end{array}$$

(24) In the same way, we can obtain the covariance matrix of backward spatial smoothing as

$$\begin{array}{cc}& {\mathbf{R}}_b=\mathbf{J}{{\mathbf{R}}_x}^{\mathit{\ast }}\mathbf{J}\mathbf{,}\end{array}$$

(25) where $\mathbf{J}={\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \cdots & \mathrm{1}\\ \mathrm{0}& \cdots & \mathrm{1}& \mathrm{0}\\ ⋮& ⋮& ⋮& ⋮\\ \mathrm{1}& \mathrm{0}& \cdots & \mathrm{0}\end{array}\right]}_{N_r\times N_r}$ . The conjugate of ${\mathbf{R}}_x$ is ${\mathbf{R}}_x^{\ast }$ .

Therefore, we can define forward-backward smoothing covariance matrix [5] as

$$\begin{array}{cc}& {\mathbf{R}}_{xx}=\frac{\mathrm{1}}{\mathrm{2}}\left({\mathbf{R}}_x+{\mathbf{R}}_b\right).\end{array}$$

(26) Overestimate operation on the processed covariance can be obtained as

$$\begin{array}{cc}& {\mathbf{R}}_{\mathrm{11}}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum }_{m=\mathrm{1}}^{\mathrm{2}}{\mathbf{R}}_{xxm}{\left({\mathbf{R}}_{xxm}^H{\mathbf{R}}_{xxm}\right)}^{-\mathrm{1}}{\mathbf{R}}_{xxm}^H,\end{array}$$

(27) where ${\mathbf{R}}_{xxm}={\mathbf{R}}_{xx}(:,m:N_r+m-\mathrm{2})$ .

The steps for the proposed algorithm are given as follows.

Step 1.

Divide the received signal $\widehat{\mathbf{Y}}$ into two subarrays: ${\widehat{\mathbf{Y}}}_{\mathrm{1}}=\widehat{\mathbf{Y}}(\mathrm{1}:N_r,:)$ and ${\widehat{\mathbf{Y}}}_{\mathrm{2}}=\widehat{\mathbf{Y}}(N_r:2N_r,:)$ .

Step 2.

Autocorrelation matrix ${\mathbf{R}}_{\mathbf{11}}$ of the subarray ${\widehat{\mathbf{Y}}}_{\mathbf{1}}$ can be obtained by using the formula (27), and cross-correlation matrix ${\mathbf{R}}_{\mathbf{21}}$ of the two subarrays can be get by using the formula (15).

Step 3.

The mixed covariance of signal subspace R can be obtained by using the formula (17).

Step 4.

The DOA is estimated by using the formula (19) and (22).

4. Simulation Results

The actual number of transmitting array and receiving array are set as $M=2$ and $N=2$ , respectively, and virtual array elements are $M=3$ , element spacing is $d_t=d_r=\lambda /2$ , and $L=\mathrm{1024}$ is code length of transmitted signal. Mutually orthogonal code signals are transmitted for each transmit array element, and the duration time of each code is $Ts=10\hspace{0.17em}\mu s$ . The duration time of the transmitted pulse is 1.024 ms. The amplitude attenuation of K target is equal to 1 [16].

Part  1. We consider there are three unrelated targets $K=3$ existing in the air and the DOA and DOD of three sources are assumed to be located at $\varphi =(-\mathrm{2}{\mathrm{5}}^{°},{\mathrm{5}}^{°},\mathrm{2}{\mathrm{2}}^{°})$ and $\theta =(-\mathrm{2}{\mathrm{2}}^{\circ },\mathrm{1}{\mathrm{8}}^{\circ },{\mathrm{5}}^{\circ })$ , respectively. The Doppler shifts are $\mathrm{100}\hspace{0.17em}\mathrm{H}\mathrm{z}$ , $\mathrm{1550}\hspace{0.17em}\mathrm{H}\mathrm{z}$ , $\mathrm{3000}\hspace{0.17em}\mathrm{H}\mathrm{z}$ , respectively, and snapshot L is 1024. All simulations are averaged over 1000 independent runs. Figure 3 shows DOD-DOA estimation performance of three signals using our proposed method. It can be seen that the proposed algorithm can accurately estimate the angle of the far spaced sources when $\mathrm{S}\mathrm{N}\mathrm{R}=\mathrm{20}\hspace{0.17em}\mathrm{d}\mathrm{B}$ .

Figure 3

The DOD-DOA estimation performance of three sources.

Part  2. The DOA and DOD of three sources are assumed to be located at $\phi =(\mathrm{2}{\mathrm{5}}^{\circ },\mathrm{2}{\mathrm{0}}^{\circ },\mathrm{2}{\mathrm{2}}^{\circ })$ and $\theta =(\mathrm{2}{\mathrm{2}}^{\circ },\mathrm{2}{\mathrm{1}}^{\circ },\mathrm{2}{\mathrm{0}}^{\circ })$ , respectively, and the other parameters are the same as before.

Figure 4 illustrates that the MIMO Array MUSIC algorithm [10] cannot distinguish the location information when sources are closed to each other, whatever there are two or three transmitting array elements. Nevertheless, it is evident from Figure 5 that we propose the modified algorithm have also been identified exactly when the three sources are very close.

Figure 4

DOD-DOA estimation performance of three sources ( $\mathrm{S}\mathrm{N}\mathrm{R}=\mathrm{25}\hspace{0.17em}\mathrm{d}\mathrm{B}$ ).

Figure 5

The DOD-DOA estimation performance of three sources. By using the modified algorithm ( $\mathrm{S}\mathrm{N}\mathrm{R}=\mathrm{25}\hspace{0.17em}\mathrm{d}\mathrm{B}$ ).

Part  3. Define the root mean squared error (RMSE) of the DOA as

$$\begin{array}{cc}& \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\frac{\mathrm{1}}{L}{\displaystyle \sum }_{l=\mathrm{1}}^L\sqrt{\frac{\mathrm{1}}{\mathrm{1000}}{{\displaystyle \sum }}_{n=1}^{\mathrm{1000}}{\left({\widehat{\theta }}_k^n-{\theta }_k\right)}^{\mathrm{2}}},\end{array}$$

(28) where ${\widehat{\theta }}_k^n$ is the estimate of DOA ${\theta }_k$ of the nth Monte Carlo simulation experiment, K is the number of sources.

Figure 6 shows the RMSE of DOA estimation as a function of input SNR, and the condition of simulation is the same as first source of Part  1. We can see that DOA estimation performance of different array spacing is very enormous when there are two real arrays or three virtual arrays. As the signal-to-noise ratio increases, the performance of the RMSE of DOA estimation can achieve the optimal when the array interval is $d_r=\lambda /2$ . So the performance of the system is best when distance between the adjacent antenna elements are placed with half-wavelength.

Figure 6

The RMSE-DOA estimation performance when different interval.

Part  4. We can see from the Figure 7 that estimation accuracy decreases as the Doppler shift increases. The bigger the Doppler frequency shift, the worse the estimation performance. However, virtual three arrays are better than two arrays for the DOA estimation accuracy. Figure 8 shows the more virtual array elements, the higher estimation accuracy when the Doppler shift is $f_d=100$ constantly. In addition, Figure 8(a) shows the RMSE of DOA estimation of the three virtual elements superior to the two actual elements with the value of SNR increasing for $L=1024$ and $f_d=100$ . And Figure 8(b) shows the RMSE of DOA estimation of the three virtual elements superior to the two actual elements with the number of snapshots increasing for $\mathrm{S}\mathrm{N}\mathrm{R}=\mathrm{20}$ and $f_d=\mathrm{100}$ .

Figure 7

The RMSE-DOA estimation performance of the three sources. In case of different Doppler shift.

Figure 8

RMSE-DOA performance under the various numbers of SNR and snapshots.

Part  5. By comparing proposed algorithm with the MIMO Array MUSIC algorithm [10], Figure 9 shows performance of the RMSE of DOA estimation as a function of input SNR for the simulation conditions are similar to Part 1, which in case of two real arrays or three virtual arrays, separately. It can be clearly seen from Figure 9 that as the SNR increases, the performance of the proposed method is much better than the MIMO Array MUSIC one, in addition, which demonstrate that the proposed algorithm can also achieve good performance under the condition of two and three transmitting array elements.

Figure 9

RMSE of DOA estimation versus SNR.

Obviously, Figure 10 also shows our proposed algorithm has higher accuracy than the MIMO Array MUSIC algorithm for estimation precision of three sources with the number of snapshots increasing, in case of the two real arrays or three virtual arrays, respectively. At the same time as the number of snapshots increases, estimation accuracy tends to a constant.

Figure 10

RMSE-DOA versus snapshot curves of three sources ( $\mathrm{S}\mathrm{N}\mathrm{R}=\mathrm{20}$ ).

5. Conclusion

The problem of DOA estimation has been an active research topic in array processing for several decades. In this paper, a new method for virtual array generation and an improved MUSIC algorithm are proposed, which can correctly estimate the two-dimensional angle of sources with limited array elements when the FBSS technique is performed on received signals. Simulation results show that when the virtual array element spacing is equal to the half wavelength, the proposed algorithm can not only estimate the DOA more accurately than the MUSIC algorithm, but also distinguish between closely spaced sources; even the difference of the angle is close to 1 degree.