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

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 L1-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.


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 1 -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 1 -SRACV algorithm for estimating sparse signals' DOAs. Compared to L 1 -SVD, the L 1 -SRACV algorithm does not need to determine the regularization parameter and has a higher stability. However, the L 1 -SRACV algorithm needs to estimate the noise power. In [7,8], L 1 -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 Sections 4 and 5 respectively.

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 θ " rθ 1 θ 2¨¨¨θK s. The array received signal is given by where A " rapθ 1 q apθ 2 q¨¨¨apθ K qs is the MˆK array manifold matrix whose steering vector is apθ k q " r1, e jφpθ k q ,¨¨¨, e jpM´1qφpθ k q s T , φpθ k q "´2π d λ sinpθ k q, λ denotes the wavelength of the incident signals, p¨q T denotes the transpose, sptq is the Gaussian signal with mean zero, nptq is additive Gaussian white noise, and L is the number of snapshots.

Signal Model Based on Sparse Representation
Suppose that the grid α " rα 1 , α 2 ,¨¨¨, α P s T (P " K) contains all the potential directions in the spatial domain. The sparse representation of the signal can be expressed as follows: where A " rapα 1 q apα 2 q¨¨¨apα P qs is an overcomplete dictionary of size MˆP, X " rx 1 , x 2 ,¨¨¨, x P s T is a sparse signal of size PˆL in the spatial domain, and x p ptq is equal to s k ptq when the kth target locates at α p and is zero otherwise. The array covariance matrix is where p¨q H denotes the conjugate transpose, R X " E XX H ( is the covariance matrix of the sparse signal X where E t¨u denotes the expectation, and R X " diag σ 2 1 , σ 2 2 ,¨¨¨, σ 2 P ( , σ 2 p is the pth spatial signal's power, σ 2 denotes the noise power, and I M denotes an identity matrix of size MˆM.

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 1 -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: where I M 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 where p¨q˚denotes the complex conjugate. The array manifold matrix of the isometric uniform linear array satisfies that J M A˚" AB and we have where the diagonal matrix B " Φ 1´M , and the rotation matrix Φ " diag ! e jφpα 1 q¨¨¨ejφpα P q ) . The real-valued symmetrical matrix can be achieved via a unitary transformation The equation B " Φ 1´M is substituted into Equation (7) and we have where Rep¨q denotes the real part and A " AΦ 1´M 2 " rapα 1 q apα 2 q¨¨¨apα P qs (9) and the real-valued overcomplete dictionary is Ψ " QA " rψpα 1 q ψpα 2 q¨¨¨ψpα P qs If M is even, the real-valued steering vector is and, if M is odd, the real-valued steering vector is 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]: where vecp¨q is the vectorization operation that stacks each column of a matrix one by one, Υ " rυ 1 , υ 2 ,¨¨¨, υ k s P C mˆk , Ξ " rξ 1 , ξ 2 ,¨¨¨, ξ k s P C nˆk , d denotes the KR product, Ξ˚d Υ " rξ1 b υ 1 , ξ2 b υ 2 ,¨¨¨, ξk b υ k s, b denotes Kronecker product, ζ " rζ 1 , ζ 2 ,¨¨¨, ζ k s T , and Z " diagpζq is a diagonal matrix. According to the above property, Equation (14) is applied to Equation (8), and Equation (8) can be vectorized as where r r " vecpR r q is the vector form of R r , and u " diag is a vector which is consist of the diagonal elements of the matrix RepΦ (15) can be considered a new observation model, in which the observation vector r r , the dictionary pΨ d Ψq, 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 pΨ d Ψq is M 2 , and the number of the distinct rows in pΨ d Ψq is 2M´1, which is greater than the matrix Ψ's dimension M. Therefore, the array aperture is extended, and the resolving power is effectively improved.

DOA Estimation Based on a SRACV
According to Equations (5) and (7), R r can also be formulated as and As the number of snapshots is finite in practice, the practical sampling covariance matrix is given byR where ∆R "R Y´RY is the estimate error. The vector ∆r " vecp∆Rq " vecpR Y´RY q satisfies an asymptotic normal distribution [17,18] ∆r " where AsNpµ, σ 2 q denotes an asymptotic normal distribution whose expectation is µ and whose variance is σ 2 . According to the literature [17][18][19], it can be known that and In order to obtain the covariance of ∆r and ∆r˚, firstly we have wherer i , r i , y i ptq denote the ith column ofR Y , the ith column of R Y and the ith element of yptq, respectively. Then, the covariance matrix of ∆r and ∆r˚is given by where covpx, yq denotes the covariance matrix of x and y, and we have covpvecp∆R r1 q, vecp∆R r2 qq " covpG 1 ∆r, From Equation (30), it can be known that vecp∆R r1 q and vecp∆R r2 q are dependent. Obviously, the sum of vecp∆R r1 q and vecp∆R r2 q still satisfies an asymptotic normal distribution wherer r " vecpR r q " vecp 1 2 QpR Y`JMRY J M qQ H q, Asχ 2 pM 2 q denotes an asymptotic chi-square with M 2 degrees of freedoms. Then, the following formula can hold with a high probability q where η " χ 2 pM 2 q, which can be obtained through the function chi2invpq, M 2 q in Matlab. Further, the formula for sparse DOA estimation is given by 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 OpP 3 q in the proposed algorithm, while it requires OpM 3 P 3 q in the L 1 -SRACV algorithm and OpK 3 P 3 q in the L 1 -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 1 -RVSKR method can be considered as the unitary L 1 -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.

Simulation Experiments
In this section, several simulations are performed with MATLAB R2014a and CVX toolbox to verify the performance of the proposed L 1 -RVSKR algorithm. Firstly, the spatial spectra of the L 1 -RVSKR method is compared with L 1 -SRACV, L 1 -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 λ{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 η, and the average value of the M´K smallest eigenvalues ofR Y is used as σ 2 in the L 1 -RVSKR and L 1 -SRACV algorithms.

The Spatial Spectra Comparison
In this simulation, the spatial spectra of the L1-RVSKR algorithm is compared with L1-SRACV, L1-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 L1-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 L1-RVSKR algorithm can distinguish each signal in different SNRs. Several pseudo-peaks appear with the L1-RVSKR and L1-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 L1-RVSKR algorithm left, which is less than that of the L1-SRACV algorithm. Moreover, the L1-SVD and MUSIC algorithms cannot distinguish the signals with low SNR, although they have no pseudo-peaks.
According to the above simulation results, the L1-RVSKR and L1-SRACV algorithms are easy to produce pseudo-peaks because they are realized by L1-norm minimization. The problem can be solved by the idea of a weighted L1-norm constraint [21,22].

The Spatial Spectra Comparison
In this simulation, the spatial spectra of the L 1 -RVSKR algorithm is compared with L 1 -SRACV, L 1 -SVD, and MUSIC. Consider four uncorrelated far-field narrowband signals arriving at the array from directions r´35˝,´10˝, 15˝, 40˝s. 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 1 -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 1 -RVSKR algorithm can distinguish each signal in different SNRs. Several pseudo-peaks appear with the L 1 -RVSKR and L 1 -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 1 -RVSKR algorithm left, which is less than that of the L 1 -SRACV algorithm. Moreover, the L 1 -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 1 -RVSKR and L 1 -SRACV algorithms are easy to produce pseudo-peaks because they are realized by L 1 -norm minimization. The problem can be solved by the idea of a weighted L 1 -norm constraint [21,22].

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: 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 L1-RVSKR is less than those of the other two methods when the SNR is less than −6 dB. It is declared that the proposed L1-RVSKR algorithm can achieve better performance than the other two methods with low SNR.

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: RMSE " where N denotes the number of independent Monte Carlo simulations, andθ kn is the estimated value of θ k in the nth 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 1 -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 1 -RVSKR algorithm can achieve better performance than the other two methods with low SNR.
value of k  in the th n 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 L1-RVSKR is less than those of the other two methods when the SNR is less than −6 dB. It is declared that the proposed L1-RVSKR algorithm can achieve better performance than the other two methods with low SNR.

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 ∆θ`8˝. The angle interval ∆θ 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 1 -SVD, and L 1 -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 1 -RVSKR decreases significantly when the angle interval is smaller than 4˝, and the performance of L 1 -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.

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 8     . The angle interval   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, L1-SVD, and L1-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 L1-RVSKR decreases significantly when the angle interval is smaller than 4 , and the performance of L1-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.

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 SNR 0dB  . 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 L1-RVSKR algorithm is smaller than the other two methods, and the L1-RVSKR algorithm has a better performance with a small number of snapshots.

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 SNR " 0dB. 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 1 -RVSKR algorithm is smaller than the other two methods, and the L 1 -RVSKR algorithm has a better performance with a small number of snapshots.

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 SNR 0dB  . 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 L1-RVSKR algorithm is smaller than the other two methods, and the L1-RVSKR algorithm has a better performance with a small number of snapshots.

The Algorithm Complexity Analysis
In this section, we perform 50 Monte Carlo simulations to verify the low computational complexity of the proposed L 1 -RVSKR algorithm. The average running time of the three algorithms with different numbers of sensors is shown in Table 1 and Figure 6. It is shown that the running time of L 1 -RVSKR is the shortest among the three algorithms. The proposed L 1 -RVSKR algorithm runs much faster than the other two algorithms and has a higher computational efficiency. Therefore, L 1 -RVSKR outperforms the other two algorithms with low computational complexity according to the above analysis and simulation results.

The Algorithm Complexity Analysis
In this section, we perform 50 Monte Carlo simulations to verify the low computational complexity of the proposed L1-RVSKR algorithm. The average running time of the three algorithms with different numbers of sensors is shown in Table 1 and Figure 6. It is shown that the running time of L1-RVSKR is the shortest among the three algorithms. The proposed L1-RVSKR algorithm runs much faster than the other two algorithms and has a higher computational efficiency. Therefore, L1-RVSKR outperforms the other two algorithms with low computational complexity according to the above analysis and simulation results.

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

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 1 -SVD and L 1 -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.