The Real-Valued Sparse Direction of Arrival (DOA) Estimation Based on the Khatri-Rao Product

Abstract

There is a problem that complex operation which leads to a heavy calculation burden is required when the direction of arrival (DOA) of a sparse signal is estimated by using the array covariance matrix. The solution of the multiple measurement vectors (MMV) model is difficult. In this paper, a real-valued sparse DOA estimation algorithm based on the Khatri-Rao (KR) product called the L₁-RVSKR is proposed. The proposed algorithm is based on the sparse representation of the array covariance matrix. The array covariance matrix is transformed to a real-valued matrix via a unitary transformation so that a real-valued sparse model is achieved. The real-valued sparse model is vectorized for transforming to a single measurement vector (SMV) model, and a new virtual overcomplete dictionary is constructed according to the KR product’s property. Finally, the sparse DOA estimation is solved by utilizing the idea of a sparse representation of array covariance vectors (SRACV). The simulation results demonstrate the superior performance and the low computational complexity of the proposed algorithm.

1. Introduction

The direction of arrival (DOA) estimation is an important area of research in array signal processing. Currently, there are many DOA estimation algorithms with good performance, such as the multiple signal classification (MUSIC) algorithm, the estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm, etc. Most of these traditional DOA estimation algorithms require the source number as a priori information and a large number of snapshots to guarantee the estimation precision. In recent years, the DOA estimation that utilizes the idea of sparse representation has become many scholars’ research hotspot [1,2]. The target signals can be regarded to be sparse in a spatial domain, and their DOAs can be estimated according to the array received data and a redundant dictionary. The DOA estimation algorithm by utilizing the idea of sparse representation mainly divides into two kinds [3,4]. One is the sparse model based on the array received data, and the other one is the sparse model based on the array covariance matrix. The singular value decomposition (SVD) of the array received data is used to propose a L₁-SVD algorithm for DOA estimation [5]. Literature [6] introduces the idea of a sparse representation of array covariance vectors (SRACV) to propose a method called the L₁-SRACV algorithm for estimating sparse signals’ DOAs. Compared to L₁-SVD, the L₁-SRACV algorithm does not need to determine the regularization parameter and has a higher stability. However, the L₁-SRACV algorithm needs to estimate the noise power. In [7,8], L₁-SRACV is improved such that the noise power is unnecessary and the algorithm’s robustness is enhanced. However, the above algorithms are all the resolution problems of the multiple measurement vectors (MMV) model and refer to complex operations and thus a heavy calculation burden. In this paper, the sparse model based on the array covariance matrix is transformed to a real-valued sparse model via a unitary transformation so that the amount of calculation is reduced by at least four times. Then, the real-valued sparse model is vectorized, and the MMV model is transformed to a single measurement vector (SMV) one whose calculation complexity is lower. Recently, the Khatri-Rao (KR) approach for DOA estimation increases rapidly. It can extend the array aperture effectively, increase the degree of freedom, and deal with the underdetermined DOA estimation in cases where the number of array elements is less than the source number in some conditions [9,10,11]. According to the property of the KR product and the real-valued overcomplete dictionary achieved by a unitary transformation, this paper constructs a new real-valued dictionary and uses the idea of a SRACV to estimate the target signals’ DOAs. The simulations demonstrate that the proposed algorithm has a good estimation performance and a low amount of calculation.

In this paper, the italic letters, the bold italic capital letters, and the bold italic lower-case letters denote variables, matrices, and vectors, respectively. The remainder of this paper is organized as follows. The sparse DOA estimation model based on the array covariance matrix is introduced in Section 2. The proposed algorithm and the detailed formula derivation of the proposed algorithm is given in Section 3. The simulation results and conclusions are given in Section 4 and Section 5, respectively.

2. Signal Model for DOA Estimation

2.1. Input Signal Model

Suppose a uniform linear array (ULA) which has $M$ elements, and the element spacing is $d$. There are $K$ uncorrelated far-field narrowband signals impinging on the array from directions $\theta =[{\theta }_1{\theta }_2\cdots {\theta }_K]$. The array received signal is given by

$$y(t)=As(t)+n(t),t=1,2,\cdots ,L$$

(1)

where $A=[a({\theta }_1)a({\theta }_2)\cdots a({\theta }_K)]$ is the $M\times K$ array manifold matrix whose steering vector is $a({\theta }_k)={[1,e^{j\varphi ({\theta }_k)},\cdots ,e^{j(M-1)\varphi ({\theta }_k)}]}^{\mathrm{T}}$, $\varphi ({\theta }_k)=-2\pi \frac{d}{\lambda }\sin ({\theta }_k)$, $\lambda$ denotes the wavelength of the incident signals, ${(·)}^{\mathrm{T}}$ denotes the transpose, $s(t)$ is the Gaussian signal with mean zero, $n(t)$ is additive Gaussian white noise, and $L$ is the number of snapshots.

2.2. Signal Model Based on Sparse Representation

Suppose that the grid $\alpha ={[{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_P]}^{\mathrm{T}}$ ($P\gg K$) contains all the potential directions in the spatial domain. The sparse representation of the signal can be expressed as follows:

$$Y=\overline{A}X+N$$

(2)

where $\overline{A}=[a({\alpha }_1)a({\alpha }_2)\cdots a({\alpha }_P)]$ is an overcomplete dictionary of size $M\times P$, $X={[x_1,x_2,\cdots ,x_P]}^{\mathrm{T}}$ is a sparse signal of size $P\times L$ in the spatial domain, and $x_p(t)$ is equal to $s_k(t)$ when the $k\text{th}$ target locates at ${\alpha }_p$ and is zero otherwise.

The array covariance matrix is

$$R_{\mathrm{Y}}=\mathrm{E}\{Y{Y}^{\mathrm{H}}\}=\overline{A}{R}_{\mathrm{X}}{\overline{A}}^{\mathrm{H}}+{\sigma }^2{I}_M$$

(3)

where ${(·)}^{\mathrm{H}}$ denotes the conjugate transpose, $R_{\mathrm{X}}=\mathrm{E}\left\{X{X}^{\mathrm{H}}\right\}$ is the covariance matrix of the sparse signal $X$ where $\text{E}\{·\}$ denotes the expectation, and $R_{\mathrm{X}}=\mathrm{diag}\{{\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_P^2\}$, ${\sigma }_p^2$ is the $p\text{th}$ spatial signal’s power, ${\sigma }^2$ denotes the noise power, and $I_M$ denotes an identity matrix of size $M\times M$.

3. The Real-Valued Sparse DOA Estimation Based on the KR Product

3.1. The Real-Valued Sparse Model for DOA Estimation

In this section, a real-valued sparse DOA estimation algorithm based on the KR product called the L₁-RVSKR is proposed. As multiplication operation occupies the most proportion in the calculation and a complex multiplication requires four real number multiplications, the complex data of the array can be transformed to real data in order to reduce the computational complexity [12,13,14,15]. Define a unitary matrix as follows:

$$Q=\{\begin{array}{ll}\frac{1}{\sqrt{2}}\left[\begin{array}{lll}I{}_{\frac{M}{2}}\hfill & & \hfill J_{\frac{M}{2}}\\ jJ{}_{\frac{M}{2}}\hfill & & \hfill -j{I}_{\frac{M}{2}}\end{array}\right]\hfill & \mathrm{if}M\mathrm{is}\mathrm{even}\hfill \\ \frac{1}{\sqrt{2}}\left[\begin{array}{lllll}I{}_{\frac{M-1}{2}}\hfill & & \hfill 0_{\frac{M-1}{2}\times 1}\hfill & & \hfill J_{\frac{M-1}{2}}\\ 0_{\frac{M-1}{2}\times 1}^{\mathrm{T}}\hfill & & \hfill \sqrt{2}\hfill & & \hfill 0_{\frac{M-1}{2}\times 1}^{\mathrm{T}}\\ jJ{}_{\frac{M-1}{2}}\hfill & & \hfill 0_{\frac{M-1}{2}\times 1}\hfill & & \hfill -j{I}_{\frac{M-1}{2}}\end{array}\right]\hfill & \mathrm{if}M\mathrm{is}\mathrm{odd}\hfill \end{array}$$

(4)

where $I_{\frac{M}{2}}$, $J_{\frac{M}{2}}$, and $0_{\frac{M-1}{2}\times 1}$ denote the $\frac{M}{2}\times \frac{M}{2}$ identity matrix, the $\frac{M}{2}\times \frac{M}{2}$ permutation matrix, and the $\frac{M-1}{2}\times 1$ null matrix, respectively.

Theorem 1. 

For any $M\times M$ Hermitian persymmetric matrix $D$, $QD{Q}^{\mathrm{H}}$ is real and symmetric.

However, as the number of snapshots in practice is finite, the practical sampling covariance matrix is Hermitian but generally not persymmetric. Therefore, the persymmetrized estimator of the practical sampling covariance matrix is

$$\begin{array}{ll}\hfill R& =\frac{1}{2}(R_{\mathrm{Y}}+J_M{R}_{\mathrm{Y}}^{*}{J}_M)\hfill \\ & =\frac{1}{2}(\overline{A}{R}_{\mathrm{X}}{\overline{A}}^{\mathrm{H}}+J_M{(\overline{A}{R}_{\mathrm{X}}{\overline{A}}^{\mathrm{H}})}^{*}{J}_M)+{\sigma }^2{I}_M\hfill \\ & =\frac{1}{2}(\overline{A}{R}_{\mathrm{X}}{\overline{A}}^{\mathrm{H}}+(J_M{\overline{A}}^{*})R_{\mathrm{X}}^{*}{(J_M{\overline{A}}^{*})}^{\mathrm{H}})+{\sigma }^2{I}_M\hfill \end{array}$$

(5)

where ${(·)}^{*}$ denotes the complex conjugate.

The array manifold matrix of the isometric uniform linear array satisfies that $J_M{\overline{A}}^{*}=\overline{A}B$ and we have

$$\begin{array}{cc}\hfill R& =\frac{1}{2}(\overline{A}{R}_{\mathrm{X}}{\overline{A}}^{\mathrm{H}}+\overline{A}B{R}_{\mathrm{X}}^{*}{B}^{\mathrm{H}}{\overline{A}}^{\mathrm{H}})+{\sigma }^2{I}_M\hfill \\ & =\frac{1}{2}\overline{A}(R_{\mathrm{X}}+B{R}_{\mathrm{X}}^{*}{B}^{\mathrm{H}}){\overline{A}}^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \end{array}$$

(6)

where the diagonal matrix $B={\Phi }^{1-M}$, and the rotation matrix $\Phi =\mathrm{diag}\{e^{j\varphi ({\alpha }_1)}\cdots e^{j\varphi ({\alpha }_P)}\}$.

The real-valued symmetrical matrix can be achieved via a unitary transformation

$$\begin{array}{cc}\hfill R_{\mathrm{r}}& =QR{Q}^{\mathrm{H}}\hfill \\ & =\frac{1}{2}Q\overline{A}(R_{\mathrm{X}}+B{R}_{\mathrm{X}}^{*}{B}^{\mathrm{H}}){\overline{A}}^{\mathrm{H}}{Q}^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \end{array}$$

(7)

The equation $B={\Phi }^{1-M}$ is substituted into Equation (7) and we have

$$\begin{array}{cc}\hfill R_{\mathrm{r}}& =\frac{1}{2}Q\overline{A}(R_{\mathrm{X}}+{\Phi }^{1-M}{R}_{\mathrm{X}}^{*}{\Phi }^{M-1}){\overline{A}}^{\mathrm{H}}{Q}^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \\ & =\frac{1}{2}Q\overline{A}{\Phi }^{\frac{1-M}{2}}{\Phi }^{\frac{M-1}{2}}(R_{\mathrm{X}}+{\Phi }^{1-M}{R}_{\mathrm{X}}^{*}{\Phi }^{M-1}){\Phi }^{\frac{1-M}{2}}{\Phi }^{\frac{M-1}{2}}{\overline{A}}^{\mathrm{H}}{Q}^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \\ & =\frac{1}{2}Q\underset{\_}{A}{\Phi }^{\frac{M-1}{2}}(R_{\mathrm{X}}+{\Phi }^{1-M}{R}_{\mathrm{X}}^{*}{\Phi }^{M-1}){\Phi }^{\frac{1-M}{2}}{\underset{\_}{A}}^{\mathrm{H}}{Q}^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \\ & =\frac{1}{2}Q\underset{\_}{A}({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}}+{({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}})}^{*}){\underset{\_}{A}}^{\mathrm{H}}{Q}^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \\ & =\frac{1}{2}\Psi ({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}}+{({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}})}^{*}){\Psi }^{\mathrm{H}}+{\sigma }^2{I}_M\hfill \\ & =\Psi \mathrm{Re}({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}}){\Psi }^{\mathrm{T}}+{\sigma }^2{I}_M\hfill \end{array}$$

(8)

where $\mathrm{Re}(·)$ denotes the real part and

$$\underset{\_}{A}=\overline{A}{\Phi }^{\frac{1-M}{2}}=[\underset{\_}{a}({\alpha }_1)\underset{\_}{a}({\alpha }_2)\cdots \underset{\_}{a}({\alpha }_P)]$$

(9)

$$\underset{\_}{a}({\alpha }_p)=\{\begin{array}{ll}{[e^{-j\frac{M-1}{2}\varphi ({\alpha }_p)},\cdots ,e^{-j\frac{1}{2}\varphi ({\alpha }_p)},e^{j\frac{1}{2}\varphi ({\alpha }_p)},\cdots ,e^{j\frac{M-1}{2}\varphi ({\alpha }_p)}]}^{\mathrm{T}}\hfill & \mathrm{if}M\mathrm{is}\text{ even}\hfill \\ {[e^{-j\frac{M-1}{2}\varphi ({\alpha }_p)},\cdots ,1,\cdots ,e^{j\frac{M-1}{2}\varphi ({\alpha }_p)}]}^{\mathrm{T}}\hfill & \mathrm{if}M\mathrm{is}\mathrm{odd}\hfill \end{array}$$

(10)

and the real-valued overcomplete dictionary is

$$\Psi =Q\underset{\_}{A}=[\psi ({\alpha }_1)\psi ({\alpha }_2)\cdots \psi ({\alpha }_P)]$$

(11)

If $M$ is even, the real-valued steering vector is

$$\psi ({\alpha }_p)=\sqrt{2}{\left[\cos (\frac{M-1}{2}\varphi ({\alpha }_p)),\cdots ,\cos (\frac{1}{2}\varphi ({\alpha }_p)),\sin (\frac{1}{2}\varphi ({\alpha }_p)),\cdots ,\sin (\frac{M-1}{2}\varphi ({\alpha }_p))\right]}^{\mathrm{T}}$$

(12)

and, if $M$ is odd, the real-valued steering vector is

$$\psi ({\alpha }_p)=\sqrt{2}{\left[\cos (\frac{M-1}{2}\varphi ({\alpha }_p)),\cdots ,\cos (\varphi ({\alpha }_p)),1,\sin (\varphi ({\alpha }_p)),\cdots ,\sin (\frac{M-1}{2}\varphi ({\alpha }_p))\right]}^{\mathrm{T}}$$

(13)

Then, an important property of the KR product is used with Equation (8). The property of the KR product can be expressed as follows [9]:

$$\mathrm{vec}(\Upsilon {\rm Z}{\Xi }^{\mathrm{H}})=({\Xi }^{*}\odot \Upsilon )\zeta$$

(14)

where $\text{vec}(·)$ is the vectorization operation that stacks each column of a matrix one by one, $\Upsilon =[{\upsilon }_1,{\upsilon }_2,\cdots ,{\upsilon }_k]\in {ℂ}^{m\times k}$, $\Xi =[{\xi }_1,{\xi }_2,\cdots ,{\xi }_k]\in {ℂ}^{n\times k}$, $\odot$ denotes the KR product, ${\Xi }^{*}\odot \Upsilon =[{\xi }_1^{*}\otimes {\upsilon }_1,{\xi }_2^{*}\otimes {\upsilon }_2,\cdots ,{\xi }_k^{*}\otimes {\upsilon }_k]$, $\otimes$ denotes Kronecker product, $\zeta ={[{\zeta }_1,{\zeta }_2,\cdots ,{\zeta }_k]}^{\mathrm{T}}$, and ${\rm Z}=\text{diag}(\zeta )$ is a diagonal matrix.

According to the above property, Equation (14) is applied to Equation (8), and Equation (8) can be vectorized as

$$\begin{array}{ll}\hfill r_{\mathrm{r}}& =\mathrm{vec}(R_{\mathrm{r}})\hfill \\ & =(\Psi \odot \Psi )·\text{diag}\{\mathrm{Re}({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}})\}+{\sigma }^2\mathrm{vec}(I_M)\hfill \\ & =(\Psi \odot \Psi )·u+{\sigma }^2\mathrm{vec}(I_M)\hfill \end{array}$$

(15)

where $r_{\mathrm{r}}=\mathrm{vec}(R_{\mathrm{r}})$ is the vector form of $R_{\mathrm{r}}$, and $u=\text{diag}\{\mathrm{Re}({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}})\}$ is a vector which is consist of the diagonal elements of the matrix $\mathrm{Re}({\Phi }^{\frac{M-1}{2}}{R}_{\mathrm{X}}{\Phi }^{\frac{1-M}{2}})$. Equation (15) can be considered a new observation model, in which the observation vector $r_{\mathrm{r}}$, the dictionary $(\Psi \odot \Psi )$, and the sparse vector $u$ are all real-valued. Furthermore, the MMV model is transformed to a SMV model by using the KR product’s property. Therefore, the computational complexity can be effectively decreased. Moreover, the dictionary is constructed via the KR product so that a new virtual array is produced [16]. The dimension of the produced virtual array $(\Psi \odot \Psi )$ is $M^2$, and the number of the distinct rows in $(\Psi \odot \Psi )$ is $2M-1$, which is greater than the matrix $\Psi$’s dimension $M$. Therefore, the array aperture is extended, and the resolving power is effectively improved.

3.2. DOA Estimation Based on a SRACV

According to Equations (5) and (7), $R_{\mathrm{r}}$ can also be formulated as

$$\begin{array}{ll}\hfill R_{\mathrm{r}}& =\frac{1}{2}Q(R_{\mathrm{Y}}+J_M{R}_{\mathrm{Y}}^{*}{J}_M)Q^{\mathrm{H}}\hfill \\ & =\frac{1}{2}Q{R}_{\mathrm{Y}}{Q}^{\mathrm{H}}+\frac{1}{2}Q{J}_M{R}_{\mathrm{Y}}^{*}{J}_M{Q}^{\mathrm{H}}\hfill \\ & =R_{\mathrm{r}1}+R_{\mathrm{r}2}\hfill \end{array}$$

(16)

where $R_{\mathrm{r}1}=\frac{1}{2}Q{R}_{\mathrm{Y}}{Q}^{\mathrm{H}}$ and $R_{\mathrm{r}2}=\frac{1}{2}Q{J}_M{R}_{\mathrm{Y}}^{*}{J}_M{Q}^{\mathrm{H}}$.

The vector form of Equation (16) can be expressed as

$$r_{\mathrm{r}}=\mathrm{vec}(R_{\mathrm{r}})=\mathrm{vec}(R_{\mathrm{r}1})+\mathrm{vec}(R_{\mathrm{r}2})$$

(17)

where

$$\mathrm{vec}(R_{\mathrm{r}1})=\frac{1}{2}\mathrm{vec}(Q{R}_{\mathrm{Y}}{Q}^{\mathrm{H}})=\frac{1}{2}(Q^{*}\otimes Q)\mathrm{vec}(R_{\mathrm{Y}})=G_1\mathrm{vec}(R_{\mathrm{Y}})$$

(18)

$$G_1=\frac{1}{2}(Q^{*}\otimes Q)$$

(19)

$$\begin{array}{cc}\hfill \mathrm{vec}(R_{\mathrm{r}2})& =\frac{1}{2}\mathrm{vec}(Q{J}_M{R}_{\mathrm{Y}}^{*}{J}_M{Q}^{\mathrm{H}})\hfill \\ & =\frac{1}{2}(Q^{*}\otimes Q)\mathrm{vec}(J_M{R}_{\mathrm{Y}}^{*}{J}_M)\hfill \\ & =\frac{1}{2}(Q^{*}\otimes Q)(J_M\otimes J_M)\mathrm{vec}(R_{\mathrm{Y}}^{*})\hfill \\ & =G_2\mathrm{vec}(R_{\mathrm{Y}}^{*})\hfill \end{array}$$

(20)

and

$$G_2=\frac{1}{2}(Q^{*}\otimes Q)(J_M\otimes J_M)$$

(21)

As the number of snapshots is finite in practice, the practical sampling covariance matrix is given by

$${\widehat{R}}_{\mathrm{Y}}=\frac{1}{L}{\displaystyle \sum _{t=1}^{L}Y{Y}^{\mathrm{H}}}=R_{\mathrm{Y}}+\Delta R$$

(22)

where $\Delta R={\widehat{R}}_{\mathrm{Y}}-R_{\mathrm{Y}}$ is the estimate error. The vector $\Delta r=\text{vec}(\Delta R)=\text{vec}({\widehat{R}}_{\mathrm{Y}}-R_{\mathrm{Y}})$ satisfies an asymptotic normal distribution [17,18]

$$\Delta r~\mathrm{AsN}(0_{M^2\times 1},\frac{1}{L}{R}_{\mathrm{Y}}^{\mathrm{T}}\otimes R_{\mathrm{Y}})$$

(23)

where $\text{AsN}(\mu ,{\sigma }^2)$ denotes an asymptotic normal distribution whose expectation is $\mu$ and whose variance is ${\sigma }^2$. According to the literature [17,18,19], it can be known that

$$\mathrm{vec}(\Delta R_{\mathrm{r}1})~\text{AsN}(0_{M^2\times 1},G_1(\frac{1}{L}{R}_{\mathrm{Y}}^{\mathrm{T}}\otimes R_{\mathrm{Y}})G_1^{\mathrm{H}})$$

(24)

and

$$\mathrm{vec}(\Delta R_{\mathrm{r}2})~\text{AsN}(0_{M^2\times 1},G_2(\frac{1}{L}{R}_{\mathrm{Y}}\otimes R_{\mathrm{Y}}^{*})G_2^{\mathrm{H}})$$

(25)

where

$$\mathrm{vec}(\Delta R_{\mathrm{r}1})=\mathrm{vec}(\frac{1}{2}Q{\widehat{R}}_{\mathrm{Y}}{Q}^{\mathrm{H}}-\frac{1}{2}Q{R}_{\mathrm{Y}}{Q}^{\mathrm{H}})=G_1\Delta r$$

(26)

and

$$\mathrm{vec}(\Delta R_{\mathrm{r}2})=\mathrm{vec}(\frac{1}{2}Q{J}_M{\widehat{R}}_{\mathrm{Y}}^{*}{J}_M{Q}^{\mathrm{H}}-\frac{1}{2}Q{J}_M{R}_{\mathrm{Y}}^{*}{J}_M{Q}^{\mathrm{H}})=G_2\Delta r^{*}$$

(27)

In order to obtain the covariance of $\Delta r$ and $\Delta r^{*}$, firstly we have

$$\begin{array}{l}\mathrm{E}\{{\widehat{r}}_i{({\widehat{r}}_j^{*})}^{\mathrm{H}}\}=\frac{1}{L^2}{\displaystyle \sum _{t_1=1}^L{\displaystyle \sum _{t_2=1}^L\mathrm{E}\{y(t_1)y_i^{*}(t_1){(y(t_2)y_j^{*}(t_2))}^{\mathrm{T}}\}}}\hfill \\ =\frac{1}{L^2}{\displaystyle \sum _{t_1=1}^L{\displaystyle \sum _{t_2=1,t_2\ne t_1}^L\mathrm{E}\{y(t_1)y_i^{*}(t_1)\}\mathrm{E}\{{(y(t_2)y_j^{*}(t_2))}^{\mathrm{T}}\}}}+\frac{1}{L^2}{\displaystyle \sum _{t_1=1}^L\mathrm{E}\{y(t_1)y_i^{*}(t_1)y_j^{*}(t_1)y^{\mathrm{T}}(t_1)\}}\hfill \\ =\frac{L^2-L}{L^2}{r}_i{r}_j^{\mathrm{T}}+\frac{1}{L}(r_i{r}_j^{\mathrm{T}}+r_j{r}_i^{\mathrm{T}})\hfill \\ =(1-\frac{1}{L})r_i{r}_j^{\mathrm{T}}+\frac{1}{L}(r_i{r}_j^{\mathrm{T}}+r_j{r}_i^{\mathrm{T}})\\ =r_i{r}_j^{\mathrm{T}}+\frac{1}{L}{r}_j{r}_i^{\mathrm{T}}\hfill \end{array}$$

(28)

where ${\widehat{r}}_i$, $r_i$, $y_i(t)$ denote the $i\text{th}$ column of ${\widehat{R}}_{\mathrm{Y}}$, the $i\text{th}$ column of $R_{\mathrm{Y}}$ and the $i\text{th}$ element of $y(t)$, respectively.

Then, the covariance matrix of $\Delta r$ and $\Delta r^{*}$ is given by

$$\begin{array}{cc}C\hfill & =\mathrm{cov}(\Delta r,\Delta r^{*})\triangleq \mathrm{E}\{(\widehat{r}-r){({\widehat{r}}^{*}-r^{*})}^{\mathrm{H}}\}=\mathrm{E}\{\widehat{r}{({\widehat{r}}^{*})}^{\mathrm{H}}\}-r{r}^{\mathrm{T}}\hfill \\ & =\left[\begin{array}{lllll}\frac{1}{L}{r}_1{r}_1^{\mathrm{T}}\hfill & & \hfill \frac{1}{L}{r}_2{r}_1^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill \frac{1}{L}{r}_M{r}_1^{\mathrm{T}}\\ \frac{1}{L}{r}_1{r}_2^{\mathrm{T}}& & \hfill \frac{1}{L}{r}_2{r}_2^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill \frac{1}{L}{r}_M{r}_2^{\mathrm{T}}\\ \hfill ⋮\hfill & & \hfill ⋮\hfill & & \hfill ⋮\\ \frac{1}{L}{r}_1{r}_M^{\mathrm{T}}& & \hfill \frac{1}{L}{r}_2{r}_M^{\mathrm{T}}\hfill & \hfill \cdots \hfill & \hfill \frac{1}{L}{r}_M{r}_M^{\mathrm{T}}\end{array}\right]\hfill \end{array}$$

(29)

where $\mathrm{cov}(x,y)$ denotes the covariance matrix of $x$ and $y$, and we have

$$\mathrm{cov}(\mathrm{vec}(\Delta R_{\mathrm{r}1}),\mathrm{vec}(\Delta R_{\mathrm{r}2}))=\mathrm{cov}(G_1\Delta r,G_2\Delta r^{*})=G_1\mathrm{cov}(\Delta r,\Delta r^{*})G_2^{\mathrm{H}}=G_{1}C{G}_2^{\mathrm{H}}$$

(30)

From Equation (30), it can be known that $\mathrm{vec}(\Delta R_{\mathrm{r}1})$ and $\mathrm{vec}(\Delta R_{\mathrm{r}2})$ are dependent. Obviously, the sum of $\mathrm{vec}(\Delta R_{\mathrm{r}1})$ and $\mathrm{vec}(\Delta R_{\mathrm{r}2})$ still satisfies an asymptotic normal distribution

$$\text{vec}(\Delta R_{\mathrm{r}})~\text{AsN}(0_{M^2\times 1},G_1(\frac{1}{L}{R}_{\mathrm{Y}}^{\mathrm{T}}\otimes R_{\mathrm{Y}})G_1^{\mathrm{H}}+G_2(\frac{1}{L}{R}_{\mathrm{Y}}\otimes R_{\mathrm{Y}}^{*})G_2^{\mathrm{H}}+2G_{1}C{G}_2^{\mathrm{H}})$$

(31)

Define $W=G_1(\frac{1}{L}{R}_{\mathrm{Y}}^{\mathrm{T}}\otimes R_{\mathrm{Y}})G_1^{\mathrm{H}}+G_2(\frac{1}{L}{R}_{\mathrm{Y}}\otimes R_{\mathrm{Y}}^{*})G_2^{\mathrm{H}}+2G_{1}C{G}_2^{\mathrm{H}}$, and we have

$$W^{-\frac{1}{2}}\text{vec}(\Delta R_{\mathrm{r}})~\text{AsN}(0_{M^2\times 1},I_{M^2})$$

(32)

and

$${\Vert W^{-\frac{1}{2}}\left[{\widehat{r}}_{\mathrm{r}}-(\Psi \odot \Psi )·u-{\sigma }^2\mathrm{vec}(I_M)\right]\Vert }_2^2~\mathrm{As}{\chi }^2(M^2)$$

(33)

where ${\widehat{r}}_{\mathrm{r}}=\text{vec}({\widehat{R}}_{\mathrm{r}})=\text{vec}(\frac{1}{2}Q({\widehat{R}}_{\mathrm{Y}}+J_M{\widehat{R}}_{\mathrm{Y}}^{*}{J}_M)Q^{\mathrm{H}})$, $\mathrm{As}{\chi }^2(M^2)$ denotes an asymptotic chi-square with $M^2$ degrees of freedoms. Then, the following formula can hold with a high probability $q$

$${\Vert W^{-\frac{1}{2}}\left[{\widehat{r}}_{\mathrm{r}}-(\Psi \odot \Psi )·u-{\sigma }^2\mathrm{vec}(I_M)\right]\Vert }_2^2\le \eta$$

(34)

where $\eta ={\chi }^2(M^2)$, which can be obtained through the function $\mathrm{chi}2\mathrm{inv}(q,M^2)$ in Matlab. Further, the formula for sparse DOA estimation is given by

$$\min {\Vert u\Vert }_{1}s.t.{\Vert W^{-\frac{1}{2}}\left[{\widehat{r}}_{\mathrm{r}}-(\Psi \odot \Psi )·u-{\sigma }^2\mathrm{vec}(I_M)\right]\Vert }_2\le \sqrt{\eta }$$

(35)

It is obvious that only the matrix $W$ is complex in Equation (35). Therefore, the operations before $W$ are all real number operations for a low calculation complexity. The computational complexity of the operations before $W$ is reduced by at least four times. In fact, solving Equation (35) through a second-order cone programming (SOCP) framework occupies most of the computational complexity in the proposed algorithm. Solving Equation (35) requires $\mathrm{O}(P^3)$ in the proposed algorithm, while it requires $\mathrm{O}(M^3{P}^3)$ in the L₁-SRACV algorithm and $\mathrm{O}(K^3{P}^3)$ in the L₁-SVD algorithm [5,6,20]. It is obvious that the computational complexity of the proposed algorithm is much lower than those of the other two algorithms.

The proposed L₁-RVSKR method can be considered as the unitary L₁-SRACV algorithm. The array covariance matrix is transformed to a real one and vectorized. By this way, a real-valued SMV model is obtained, and the calculation complexity is decreased. Moreover, the KR product is used to construct a new virtual dictionary that increases the degree of freedom.

4. Simulation Experiments

In this section, several simulations are performed with MATLAB R2014a and CVX toolbox to verify the performance of the proposed L₁-RVSKR algorithm. Firstly, the spatial spectra of the L₁-RVSKR method is compared with L₁-SRACV, L₁-SVD, and MUSIC. Then, several experiments are performed to compare the estimation performance of different algorithms versus the signal-to-noise ratio (SNR), the angle interval, and the number of snapshots. Finally, the low computational complexity of the proposed algorithm is verified and analyzed by comparing the running time with the other two algorithms. We consider a ULA with $\lambda /2$ spacing in the following simulations, and the number of sensors is 8 except that shown in Figure 1 and the last simulation. The grid is divided in the range of $-90°$ to $90°$ spacing $1°$. The probability is set to be $q=0.999$ to achieve the parameter $\eta$, and the average value of the $M-K$ smallest eigenvalues of ${\widehat{R}}_{\mathrm{Y}}$ is used as ${\sigma }^2$ in the L₁-RVSKR and L₁-SRACV algorithms.

4.1. The Spatial Spectra Comparison

In this simulation, the spatial spectra of the L₁-RVSKR algorithm is compared with L₁-SRACV, L₁-SVD, and MUSIC. Consider four uncorrelated far-field narrowband signals arriving at the array from directions $[-35°,-10°,15°,40°]$. The number of sensors is 5, and the number of snapshots is 500. Figure 1 shows the spatial spectra of different methods with SNR = 0 dB. It can be seen from Figure 1 that the spatial resolution of the L₁-RVSKR algorithm is better than those of the other algorithms, and it has no pseudo-peaks in such a simulation condition. In Figure 2, we consider two uncorrelated far-field narrowband signals arriving at the array from directions $8°$ and $17°$. The number of sensors is 8, and the number of snapshots is 300. Figure 2 shows that the L₁-RVSKR algorithm can distinguish each signal in different SNRs. Several pseudo-peaks appear with the L₁-RVSKR and L₁-SRACV algorithms when SNR is high, but the pseudo-peaks are lower than that of the real estimated DOA. With the decrease of the SNR, there is only one pseudo-peak of the L₁-RVSKR algorithm left, which is less than that of the L₁-SRACV algorithm. Moreover, the L₁-SVD and MUSIC algorithms cannot distinguish the signals with low SNR, although they have no pseudo-peaks.

According to the above simulation results, the L₁-RVSKR and L₁-SRACV algorithms are easy to produce pseudo-peaks because they are realized by L₁-norm minimization. The problem can be solved by the idea of a weighted L₁-norm constraint [21,22].

4.2. The Estimation Performance versus SNR

In this section, the estimation performance of different methods versus the SNR is compared by an experiment. Firstly, we define the root mean square error (RMSE) as follows:

$$\mathrm{RMSE}=\sqrt{{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^K\frac{{({\widehat{\theta }}_{kn}-{\theta }_k)}^2}{KN}}}}$$

(36)

where $N$ denotes the number of independent Monte Carlo simulations, and ${\widehat{\theta }}_{kn}$ is the estimated value of ${\theta }_k$ in the $n\text{th}$ simulation.

Similarly, we consider two uncorrelated far-field narrowband signals from directions $8°$ and $17°$. The number of snapshots is 300. The SNR varies from −10 dB to 10 dB with 2 dB steps. For each SNR, 50 Monte Carlo simulations are performed to verify the performance of the proposed method. The simulation result is shown in Figure 3. It is shown that the RMSEs of the three algorithms decrease with the increase of the SNR. The RMSE of L₁-RVSKR is less than those of the other two methods when the SNR is less than −6 dB. It is declared that the proposed L₁-RVSKR algorithm can achieve better performance than the other two methods with low SNR.

4.3. The Angle Resolution Capability

In this section, the angle resolution capabilities of different algorithms are compared. We suppose that there are two uncorrelated signals whose directions are $8°$ and $\Delta \theta +8°$. The angle interval $\Delta \theta$ varies from $2°$ to $18°$ with a step of $2°$. For each angle interval, 50 Monte Carlo simulations are performed to compare the angle resolution capability of the proposed method with MUSIC, L₁-SVD, and L₁-SRACV. The SNR is set to 0 dB, and the number of snapshots is 300. Figure 4 shows the comparison of the angular resolution of different methods. It can be seen from Figure 4 that the RMSEs of the four methods decrease with the increase of the angle interval. The RMSE of L₁-RVSKR decreases significantly when the angle interval is smaller than $4°$, and the performance of L₁-RVSKR is superior to the other methods when the angle interval is smaller than $6°$. This simulation illustrates that the resolving ability of the proposed method is superior to those of the other three methods.

4.4. The Estimation Performance versus the Number of Snapshots

In this section, we analyze the relationship of the RMSE of the DOA estimation and the number of snapshots in the case of $\text{SNR}=0\text{dB}$. Consider two uncorrelated far-field narrowband signals impinging on the array from directions $8°$ and $17°$. The number of snapshots varies from 100 to 350 with a step of 50. For each number of snapshots, 50 Monte Carlo simulations are performed to verify the performance of the proposed method. It can be seen from Figure 5 that the RMSEs of the three methods decrease as the increase of the number of snapshots. It is illustrated in this simulation that the RMSE of the L₁-RVSKR algorithm is smaller than the other two methods, and the L₁-RVSKR algorithm has a better performance with a small number of snapshots.

4.5. The Algorithm Complexity Analysis

In this section, we perform 50 Monte Carlo simulations to verify the low computational complexity of the proposed L₁-RVSKR algorithm. The average running time of the three algorithms with different numbers of sensors is shown in

Table 1. The running time versus the number of sensors.

Number of Sensors	L1-RVSKR	L1-SRACV	L1-SVD
6	3.2773 s	36.9106 s	12.9263 s
8	3.4291 s	89.6315 s	12.2859 s
10	5.8132 s	183.1054 s	12.3405 s

and Figure 6. It is shown that the running time of L₁-RVSKR is the shortest among the three algorithms. The proposed L₁-RVSKR algorithm runs much faster than the other two algorithms and has a higher computational efficiency. Therefore, L₁-RVSKR outperforms the other two algorithms with low computational complexity according to the above analysis and simulation results.

5. Conclusions

In this paper, we propose a real-valued sparse DOA estimation algorithm by using the KR product. The sparse model of the array covariance matrix is transformed to a real-valued one via a unitary transformation. Therefore, the calculated amount is reduced. Moreover, the vectorization is made to transform the real-valued MMV model to a SMV one. The array aperture is extended, and the estimation accuracy is improved by using the KR product. Finally, the idea of a SRACV is utilized for DOA estimation. The simulation results show that the proposed method can achieve better performance than L₁-SVD and L₁-SRACV with a low SNR, a small angle interval, and a small number of snapshots. The proposed method also improves the computational efficiency. However, the algorithm cannot deal with an underdetermined DOA estimation because the sparse model used is a model of SMV. We will work on this aspect in the future.