DOA Estimation Based on Sparse Signal Recovery Utilizing Double-Threshold Sigmoid

This paper proposes a new algorithm based on sparse signal recovery for estimating the direction of arrival (DOA) of multiple sources. The problem model we build is about the sample covariance matrix fitting by unknown source powers. We enhance the sparsity by the double-threshold sigmoid penalty function which can approximate the l 0 norm accurately. Our method can distinguish closely spaced sources and does not need the knowledge of the number of the sources. In addition, our method can also perform well in low SNR. Besides, our method can handle more sources accurately than other methods. Simulations are done to certify the great performance of the proposed method.


Introduction
The estimation of the direction of arrival (DOA) of multiple sources plays a key role in many applications including radar, sonar, and wireless communication.So far, amounts of superresolution algorithms for DOA estimation have been developed.The nonparametric methods include Capon method [1] and subspace-based methods.The traditional subspace-based algorithms, like MUSIC which firstly exploits the orthogonality between the signal space and the noise subspace [2] and ESPRIT which utilizes the rotational invariance of the signal subspace [3], can achieve excellent performance in high SNR, specially, when the snapshots are long.The maximum likelihood (ML) methods including the deterministic maximum likelihood (DML) and stochastic maximum likelihood (SML) possess good statistical properties [4][5][6][7] but require a large number of samples.All the above methods need the knowledge of the number of sources.
Sparse representation of signals and compressed sensing have become a hot topic in many fields [8,9], and the DOA estimation methods based on sparse reconstruction have already been paid more attention by researchers.The wellknown  1 -SVD is a pretty good algorithm [10]; it combines the sparse signal recovery method based on  1 norm with the singular value decomposition (SVD).The  1 -SVD algorithm can handle closely spaced correlated sources when the number of sources is known, while its performance is degraded without knowing the number of sources.In [11], the joint  0 approximation DOA (JLZA-DOA) algorithm is proposed.This algorithm processes with a mixed  2,0 approximation approach; it can acquire high resolution without the knowledge of the number of sources and with only a few snapshots.However, the JLZA-DOA algorithm may fail to make a good estimation in low SNR.
The algorithms we mentioned in previous paragraph are all based on the multiple measurement vectors (MMV) problem [12].Recently, the coprime array technique is proposed [13,14].This technique can enhance the degrees of freedom (DOFs) of array.With vectorizing the sample covariance matrix, it deals with the sparse covariance fitting problem eventually when signals are sparsely represented.That is to say, the MMV problem is transformed into the single measurement vector (SMV) problem.In the SMV case, DOA estimation can be implemented using many techniques, such as the orthogonal matching pursuit (OMP) and the improved smoothed  0 approximation algorithm (ISL0) [15,16].The OMP algorithm enjoys a short computational time.However, it needs to know the sparsity in advance; that is, it requires the knowledge of the number of sources.The ISL0 performs better than OMP, but it costs longer computational time.And its performance is degraded in low SNR.In this paper, we adopt the coprime array technique and propose a new Newton-like algorithm based on doublethreshold sigmoid penalty for handling the sparse covariance fitting problem.We know that the direct  0 norm optimization problem is NP-hard.Many algorithms approximate the  0 norm by  1 norm, but the estimation errors increase when the magnitudes of the nonzero elements to be estimated are greater than one.In addition, the  1 norm method is not robust to noise and takes lots of iteration to converge.Instead of replacing the  0 norm with  1 norm, we utilize the doublethreshold sigmoid penalty function to approach the  0 norm [17].And by adjusting the upper threshold at each step on the iteration, our algorithm can preserve most of the advantages of  0 norm; this contributes to improving the estimation performance.Besides it also can accelerate the speed of convergence by adjusting the upper threshold.Numerical simulations demonstrate that our proposed method can achieve high resolution without the knowledge of the number of sources and perform well in low SNR.

Problem Model
Consider a coprime array with compressed interelement spacing, as illustrated in Figure 1.The two subarrays consist of  and  sensors, respectively, where  =  + 1.  is usually set to /2, where  denotes the wavelength.Consequently, the array sensors are positioned at the two subarrays share the first sensor at zeroth position; denote p = [ 1 ,  2 , . . .,  +−1 ]  as the positions of the array sensors, where   ∈ P, ∀.
The covariance matrix of receive data vector y() can be obtained as where R s = [s()s  ()] and I +−1 represents an identity matrix.
The classical DOA estimation problem can be reformulated as a sparse representation problem.We consider a gird of  equally interval angles Φ = [ 1 , . . .,   ] with  ≫ .Assuming Θ is a subset of Φ, we construct an overcomplete matrix A by collecting the steering vectors corresponding to all the potential source locations.Accordingly, received signal model ( 2) can be represented as where A = [a( 1 ), a( 2 ), . . ., a(  )], x() ∈ C ×1 is a sparse vector, and the nonzero entries of x() are the positions which correspond to the source locations.That is to say, the th Consequently, we can obtain a spatial covariance matrix in terms of A; it takes the following form: where In practice, the spatial covariance matrix is replaced by the sample covariance matrix R = (1/) ∑  =1 y()y  (), where  denotes the sample snapshots.Under the assumption that sources are uncorrelated, R x ∈ C × is a sparse diagonal matrix, and only  entries are nonzero.Denote diag(R x ) = [ 1 ,  2 , . . .,   ]  = b, and it is obvious that b is a  × 1 sparse vector with nonzero entries at positions which correspond to source locations.Moreover, the elements of b are real valued and negative; that is, b ∈ R ×1 + .
To avoid the calculation of complex number and reduce the computational cost, we can separate the real and image part.Then (7) can be reformulated as where z  = Re(z), z  = Im(z), Ã = Re( Ã), and Ã = Im( Ã).
Our proposed method can work on both real valued case and complex valued case; we just take the complex valued case into account in the following discussion.

Proposed Method
3.1.Proposed L0A Algorithm.The problem expressed in (7) can be solved by the  1 regularization method whose model can be written as where ‖b‖ 1 is the penalty term which approximates the  0 norm.ℎ denotes the regularization parameter; it controls the tradeoff between the sparsity of the signal and the residual energy.In the following, we derive a gradient-based method with double-threshold sigmoid (DTHS) penalty.Before deriving our method, we first deal with the problem that the derivation of |  | is not defined at zero, which always happens when handling the sparse problem.Some approximation techniques for this problem are adopted in practice [18].Now we make an approximation to |  | as Obviously, this equation is also held in the complex valued case.
In order to approximate the  0 norm more accurately, some thresholding penalties have been utilized.Then (10) can be rewritten as where Th(⋅) denotes the thresholding function and ) is an approximation of ‖b‖ 0 .Now, we adopt a DTHS function in our method.The ideal one is divided into three parts by the upper and lower threshold points, which is shown in Figure 2   are multiplied by a particular constant, and the values larger than   are set to a particular positive value.However, the derivative in the two threshold points does not exist, so we replace the ideal DTHS function by the flowing function given by which can approximate the ideal DTHS function greatly when  → +∞; we name this function the DTHS function and set  to 50 in this paper.Substituting the thresholding function in ( 12) by ( 13), we obtain a new problem: Figure 3 shows the graph of the DTHS function with four different upper thresholds.For any given  > 0, we have lim Consequently, the function ∑  =1 Th(  ) behaves like ‖‖ 0 .That is to say, the DTHS penalty will approach  0 norm when a smaller   is used.However, while we choose a smaller   , the function    ,  (  ) might have many local minima.As   increases,    ,  (  ) becomes smoother.Then, we can handle our problem by solving a sequence optimization problem.Start solving (14) with a larger   ; subsequently, we reduce   by a small factor  and solve (14) again for   =   .
We can obtain the gradient of (11) as where W is a diagonal matrix: Th  (  ;   ,   ) Instead of calculating the Hessian matrix of (12), we adopt H = Ã Ã + ℎW to approximate it; that is, ∇ 2 b (b, ℎ) ≈ H. Then the iteration process can be accelerated, and we obtain the th quasi-Newton iteration formula of problem (8), which is written as where  denotes the step size.The whole algorithm is summarized in Algorithm 1; we call our method the  0 approximation (L0A) algorithm.
Output.The solution is b () after  iterations.
We initialize b 0 by b 0 = Ã+ z for accelerating the iteration, where Ã+ denotes the Moore-Penrose pseudoinverse of Ã.As discussed above, if  0 is too small, function (15) may have many local minima.This will degrade the estimation performance of our method.So  0 should be initialized by a larger value.To preserve most of the advantages of  0 norm,  0 also should not be chosen to be too large.Consequently, the value of  0 should be selected moderately according to the signal magnitude.Usually, we set  0 = 5 max(abs(b 0 )), where abs(⋅) denotes the absolute value function.As for the lower threshold   , some entries close to zero will be picked into the nonzero components of b () when   > 0. Consequently, we fixed it as 0. According to numerous experiments, we suggest that  should be set from 0.3 to 0.5, and  is set from 0.6 to 0.8.In this paper we set  stop to 0.01.Now, we talk about the selection of the regularization parameter ℎ.Firstly, we have estimated an approximate range by calculating the root mean square error (RMSE) of our proposed method as a function of ℎ, as shown in Figure 4.It demonstrates that the low RMSE can be achieved when 0.01 < ℎ < 0.1 and 0.6 < ℎ < 15.However, if ℎ < 5, there will be false peaks.And it is going to worsen as ℎ decreases; the speed of convergence also becomes slower simultaneously.When ℎ > 10, some true peaks may disappear, and the performance will be more terrible as ℎ increases.Consequently, we set ℎ from 5 to 10.Some more simulations about how ℎ affects the estimation performance are conducted in Section 4.

Applying MUSIC Algorithm. Now we make a brief introduction about how to apply MUSIC algorithm correctly for DOA estimation under the problem model
Ĩ which is obtained by vectorizing R in (4); more details can be seen in [19].It is obvious that the virtual source signal becomes a single snapshot of b.And the rank of the covariance matrix of z, R z = zz  , is one in the noise-free case.Then the problem resembles dealing with fully coherent sources.Consequently, MUSIC will fail to work when multiple sources are impinging the array.As described in [19], the spatial smoothing technique can be used to overcome this problem.Since spatial smoothing requires a continuous set of differences, we construct a new matrix Ã1 of size 2( − 1) + 1 ×  where we have extracted precisely those rows from Ã which correspond to the 2( − 1) + 1 successive differences and also sorted them.This is equivalent to removing the corresponding rows from the observation vector z and sorting them to get a new vector z 1 expressed as We divide this virtual array into (−1)+1 overlapping subarrays and construct a full-rank covariance matrix so that the MUSIC algorithm can be applied for DOA estimation.

Simulation Results
In this section, we illustrate the simulation results of our proposed method.We consider the coprime arrays talked about in Section 2 consisting of 12 array sensors,  = 6 and  = 7.
We assume four uncorrelated sources impinging on the array; their locations are 10 ∘ , 14 ∘ , 50 ∘ , and 61 ∘ .The SNR is set to 5 dB and the number of snapshots is 200.The estimated spectrum is shown in Figure 5 with above conditions, where the grid resolution is 0.5 ∘ .It demonstrates that L0A method performs well in low SNR and can distinguish closely spaced sources.We also consider two correlated sources at 21 ∘ and 30 ∘ based on the same conditions.As it is shown in Figure 6, the estimated results are close to the actual angles with small errors.
By adopting the coprime arrays, higher degrees of freedom can be achieved.Now we consider 25 uncorrelated sources with uniform space between −60 ∘ and 60 ∘ .We set the number of snapshots to 500 and the SNR is 5 dB, and the grid resolution is set to 0.2 ∘ .We compare the performance of L0A and MUSIC method under the same condition.As shown in Figures 7 and 8, the L0A method apparently performs better than MUSIC method in low SNR; it can recognize all the closely spaced sources greatly, while the MUSIC method fails to solve some angles of sources.
Figure 9 shows the RMSEs of the L0A, OMP, ISL0, and MUSIC as a function of the SNR.Twenty uncorrelated sources with 500 snapshots are taken into account in the simulations, and the grid resolution is set to 0.5 ∘ .It manifests that the estimation performance of these algorithms degrades with SNR decreasing.Our proposed L0A enjoys a better performance than others, especially in low SNR.In Figure 10, we make a comparison among their performance of DOA estimation under different snapshots in 5 dB.With the snapshots increasing, the performance becomes better.However, to achieve the same RMSE, our proposed method just needs smaller snapshots.Our simulations were run on the computer with a 3.40 GHz Intel (R) Core (TM) i7-2600 CPU, where the operating system is Microsoft Windows XP with 32 bits.Table 1 shows the computation time for different algorithms.Here, we consider the same case used in Figure 9.We can see that SNR has little effect on the computation time of OMP and MUSIC.As SNR becomes low, the computation time of L0A and ISL0 increases.However, our proposed L0A costs shorter computation time than ISL0.Finally, we show the effect of ℎ on the performance of DOA estimation from Figures 11 and 12.We consider the same case used in Figure 7. Figure 11 shows the case where ℎ = 1.Obviously, there are some false peaks in the spectrum.When ℎ = 15, true peaks disappear at the position of −50 ∘ and 55 ∘ as shown in Figure 12.Consequently, an appropriate value of ℎ is important to the DOA estimation.

Conclusion
In this paper, we propose a new method based on sparse reconstruction for finding the directions of sources impinging on a coprime array.By approximating the  0 norm with the DTHS function and adjusting the upper threshold dynamically, our method achieves an excellent performance.The proposed method not only can perform well without the knowledge of the number of sources but also could work on correlated sources with small bias.Furthermore, it can also distinguish the closely spaced sources.With the high DOFs, this method could resolve much more the number of sources   and have a better performance in low SNR and with smaller snapshots.

Figure 4 :
Figure 4: RMSE as a function of ℎ.
2  + , where  is a small positive constant, and √ 2  +  approaches |  | when  → 0 + .It is set to 0.0001 in this paper.Then