Direction of arrival (DOA) estimation of multiple emitting sources is an important issue in array processing and has various applications in military, radar, sonar, wireless communications and source localization [1], [2]. A large number of solutions have been proposed to solve this problem during the past years. Usually, these solutions can be categorized into three groups: timedelay based methods, beamforming methods and signal subspace methods. However, the majority of DOA estimation algorithms are developed under certain assumptions: the source signal needs to be statistically stationary and uncorrelated, the number of snapshots is sufficient, and the signalnoise ratio (SNR) is moderately high. Practically, these conditions are barely satisfied, thus these methods achieve the limited estimation accuracy. In order to increase the DOA estimation accuracy, the wellknown subspacebased method of multiple signal classification (MUSIC) algorithm and estimation method of signal parameters via rotational invariance techniques (ESPRIT) have been widely used due to their high estimation accuracies but at the price of high complexity.
Recently, the sparse representation technique of signal has been applied in many areas, such as image processing, wireless channel estimation and biomedical signal processing, which also provides a new idea for DOA estimation based on the fact that the number of sources is in general much smaller than that of potential source points when implementing the array processing algorithms. Several DOA estimation methods based on sparse representation have been proposed [3][16]. In [3], [4], a whiten sparse covariancebased representation model is first presented for the source parameter estimation by applying the global matched filter (GMF). In [5] the most representative sparse recovery algorithm for DOA estimation (
To make use of the second order statistics of the array output, a sparse iterative covariancebased estimation (SPICE) approach for array signal processing by the minimization of a covariance matrix fitting criterion which can be used in both single and multiple measurement cases is proposed in [10]. Another method called
All the above mentioned sparse representation based DOA estimation algorithms assume that the ambient noise is Gaussian distributed. However, the noise in practice often exhibits nonGaussian properties, sometimes accompanied by strong impulsiveness [17]. For example, atmospheric noise (thunderstorms), car ignitions, microwave ovens, office equipments, and other types of naturally occurring or manmade signal sources can result in aggregating noise components that may exhibit high amplitudes for small time intervals. Under investigation, it is found that
An important characteristic of the
The
$ \phi(t)=e^{\{jat\gammat^{\alpha}[1+j\beta {\rm sgn}(t)\varpi(t, \alpha)]\}} $  (1) 
where
Consider the case of
$ X(t)=A(\theta)S(t)+N(t) $  (2) 
where
$ \begin{align*} &X_{M\times1}(t)=[x_1(t), x_2(t), ..., x_M(t)]^T\\ &A_{M\times K}(\theta)=[{a}(\theta_1), {a}(\theta_2), ..., {a}(\theta_K)]\\ &S_{K\times 1}(t)=[s_1(t), s_2(t), ..., s_K(t)]^T\\ &N_{M\times1}(t)=[n_1(t), n_2(t), ..., n_M(t)]^T \end{align*} $ 
where
$ {a}(\theta_n)=\left[1, e^{j\frac{2\pi}{\lambda}d\sin\theta_n}, ..., e^{j\frac{2\pi}{\lambda}(M1)d\sin\theta_n}\right]^T $  (3) 
where
Assume the noise in (2) is zeromean Gaussian white noise with the power of
$ R=E(X(t)X^H(t))=A(\theta)R_sA^H(\theta)+\sigma_n^2 I_M $  (4) 
where the source covariance matrix
Applying the vectorization operator in (4), we have [22]
$ y=\mbox{vect}(R)=B(\theta)\sigma_s+\sigma_n^2\mbox{vect}(I_M) $  (5) 
$ B(\theta)=[a^\ast(\theta_1)\otimes a(\theta_1), ..., a^\ast(\theta_K)\otimes a(\theta_K)] $  (6) 
where
$ y=B(\hat{\theta})\nu+\sigma_n^2\mbox{vect}(I_M). $  (7) 
Hence the DOA estimation can be reduced to the detection of the nonzero elements of
$ \min\\hat{\nu}\_1, \quad {\rm s.t.}~\\hat{y}B(\hat{\theta})\hat{\nu}\sigma_n^2\mbox{vect}(I_M)\_2\leq \varepsilon $  (8) 
where
$ \Delta y=\hat{y}y=\mbox{vect}(\hat{R}R)\sim \mbox{AsN}\left(0_{M^2, 1}, \frac{1}{N}R^{T}\otimes R\right). $  (9) 
Define the weighted matrix
$ W^{\frac{1}{2}}\Delta y\sim \mbox{AsN}\left(0_{M^2, 1}, I_{M^2}\right). $  (10) 
Then from (7), we can further get
$ \left\W^{\frac{1}{2}}\left[\hat{y}B(\hat{\theta})\hat{\nu} \sigma_n^2\mbox{vect}(I_M)\right]\right\_2^2 \sim \mbox{As}\chi^2(M^2) $  (11) 
where
$ \left\W^{\frac{1}{2}}\left[\hat{y}B(\hat{\theta})\hat{\nu} \sigma_n^2\mbox{vect}(I_M)\right]\right\_2^2\leq \varepsilon^2 $ 
with a high probability
$ \min\\hat{\nu}\_1, \quad {\rm s.t.}~ \left\\hat{W}^{\frac{1}{2}}\left[\hat{y}B(\hat{\theta})\hat{\nu} \hat{\sigma}_n^2\mbox{vect}(I_M)\right]\right\_2\leq \varepsilon. $  (12) 
This DOA estimation algorithm based on the sparse representation of the second order statistics covariance vector can be named as SS_SOSCV algorithm.
Ⅳ. DOA ESTIMATION BASED ON SPARSE REPRESENTATION OF FRACTIONAL LOWER ORDER STATISTICS VECTORWhen the noise in (2) is
The FLOM matrix
$ C_{i, k}={\rm E}\{x_i(t)x_k(t)^{p2}x_k^\ast(t)\} $  (13) 
where
$ C=A(\theta)\Lambda_sA^H(\theta)+\xi I_M $  (14) 
where the diagonal matrix
$ y_{\rm FLOM}=\mbox{vect}(C)=B(\theta)\gamma_s+\xi\mbox{vect}(I_M). $  (15) 
The vector
$ y_{\rm FLOM}=B(\hat{\theta})\nu_{\rm FLOM}+\xi\mbox{vect}(I_M). $  (16) 
As with the SS_SOSCV algorithm, the DOA estimation can be resolved by the following convex optimization problem:
$ \begin{align} &\min\\hat{\nu}_{\rm FLOM}\_1 \\ & {\rm s.t.}\quad \left\\hat{W}^{\frac{1}{2}}_{\rm FLOM}\left[ \hat{y}_{\rm FLOM}B(\hat{\theta})\hat{\nu}_{\rm FLOM}\hat{\xi}\mbox{vect}(I_M)\right]\right\_2 \leq\varepsilon \end{align} $  (17) 
where the weighted matrix
$ \hat{W}^{\frac{1}{2}}_{\rm FLOM}=\sqrt{N}\hat{C}^{\frac{T}{2}}\otimes\hat{C}^{\frac{1}{2}} $  (18) 
$ \hat{C}_{i, k}=\frac{1}{N}\sum\limits_{t=1}^N\left\{x_i(t)x_k(t)^{p2}x_k^\ast(t)\right\} $  (19) 
From the FLOM definition, it can be seen that it is limited in range of
$ \begin{equation} \Gamma_{i, k}={\rm E}\left\{x_i^{\langle b\rangle}(t)x_k^{\langle b\rangle}(t)\right\}, \quad 0 <b <\frac{\alpha}{2} \end{equation} $  (20) 
where the PFLOM operation on a complex number
$ \begin{align} z^{\langle b\rangle}= \begin{cases} \frac{z^{b+1}}{z^\ast}, &z\neq 0\\ 0, &z=0 \end{cases} \end{align} $  (21) 
and the conjugate of the
$ \Gamma=A(\theta)\Phi_s A^H(\theta)+\kappa I_M $  (22) 
where the diagonal matrix
$ y_{\rm PFLOM}=\mbox{vect}(\Gamma)=B(\theta)\varphi_s+\kappa \mbox{vect}(I_M) $  (23) 
$ y_{\rm PFLOM}=B(\theta)\nu_{\rm PFLOM}+\kappa \mbox{vect}(I_M). $  (24) 
Likewise, the DOAs can be estimated by solving the following optimization problem:
$ \begin{align} &\min\\hat{\nu}_{\rm PFLOM}\_1 \\ & {\rm s.t.}\ \Big\\hat{W}_{\rm PFLOM}^{\frac{1}{2}}\Big[ \hat{y}_{\rm PFLOM}B(\hat{\theta})\hat{\nu}_{\rm PFLOM} \\ &\qquad\hat{\kappa} \mbox{vect}(I_M)\Big]\Big\_2 \leq \varepsilon \end{align} $  (25) 
where the weighted matrix
$ \hat{W}_{\rm PFLOM}^{\frac{1}{2}}=\sqrt{N}\hat{\Gamma}^{\frac{T}{2}}\otimes \hat{\Gamma}^{\frac{1}{2}} $  (26) 
$ \hat{\Gamma}_{i, k}=\frac{1}{N}\sum\limits_{t=1}^N\left\{x_i^{\langle b\rangle}(t)x_k^{\langle b\rangle}(t)\right\} $  (27) 
Equations (16) and (24) can be unified and expressed as
$ y_{\rm FLOS}=B(\hat{\theta})\nu_{\rm FLOS}+(\xi\kappa)\mbox{vect}(I_M). $  (28) 
Notice that the vector
$ \begin{align} y_{\rm IFLOS}=&\ Jy_{\rm FLOS}\\ =&\ J\left\{B(\hat{\theta})\nu_{\rm FLOS}+(\xi\kappa)\mbox{vect}(I_M)\right\}\\ =&\ D(\hat{\theta})\nu_{\rm FLOS} \end{align} $  (29) 
where
$ J=[J_1, J_2, ..., J_{M1}]^T $  (30) 
where
$ \begin{align} J_m=&\ [\pmb{e}_{(m1)(M+1)+2}, \pmb{e}_{(m1)(M+1)+3}, \\ &\ ..., \pmb{e}_{m(M+1)}] \in\mathbb{R}^{M^2\times M}, \\ &\qquad\qquad\qquad\qquad m=1, ..., M1 \end{align} $  (31) 
$ \min\\nu_{\rm FLOS}\_1, \quad {\rm s.t.}~\left\y_{\rm IFLOS}D(\hat{\theta})\nu_{\rm FLOS}\right\_2 <\varepsilon_I. $  (32) 
Letting
$ \Delta y_{\rm IFLOS}\sim \mbox{AsN}\left(0_{M(M1), 1}, \frac{1}{N} J\left((C\Gamma)^T\otimes (C\Gamma)\right)J^T\right). $  (33) 
Define the weighted matrix
$ {W}_{\rm IFLOS}^{\frac{1}{2}}=\sqrt{N} J^{\frac{1}{2}}\left((C\Gamma)^{\frac{T}{2}}\otimes (C\Gamma)^{\frac{1}{2}}\right)J^{\frac{T}{2}} $  (34) 
and its estimation as
$ W_{\rm IFLOS}^{\frac{1}{2}}\Delta y_{\rm IFLOS}\sim \mbox{AsN}\left(0_{M(M1), 1}, I_{M(M1), 1}\right) $  (35) 
$ \left\W_{\rm IFLOS}^{\frac{1}{2}}\left[\hat{y}_{\rm IFLOS} D(\hat{\theta})\hat{\nu}_{\rm FLOS}\right]\right\_2^2 \sim \mbox{As}\chi^2[M(M1)]. $  (36) 
Therefore, a parameter
$ \left\W_{\rm IFLOS}^{\frac{1}{2}}\left[\hat{y}_{\rm IFLOS} D(\hat{\theta})\hat{\nu}_{\rm FLOS}\right]\right\_2^2\leq \varepsilon_I^2 $ 
with a high probability
$ \varepsilon_I=\sqrt{\chi^2_{pc}(M(M1))} . $ 
Likewise, the DOAs can be estimated by solving the following optimization problem:
$ \begin{align} &\min\\hat{\nu}_{\rm FLOS}\_1\\ & {\rm s.t.}\quad \Big\\hat{W}_{\rm IFLOS}^{\frac{1}{2}}\big[\hat{y}_{\rm IFLOS} D(\hat{\theta})\hat{\nu}_{\rm FLOS}\big]\Big\_2\leq \varepsilon_I. \end{align} $  (37) 
The DOA estimation methods based on (37) by using the FLOM matrix
The main computational costs of the SS_FLOMV and SS_PFLOMV algorithms include the calculation of the FLOM matrix
From the above analysis, the SS_FLOMV and SS_ PFLOMV algorithms' steps can be summarized as following:
Step 1: Obtain the FLOM estimate matrix
Step 2: Get the estimation of the parameter
Step 3: Calculate the weighted matrix
Step 4: Solve the convex optimization problem of (17) or (25) to get the estimation of the vector
Step 5: Estimate the DOAs according the location of nonzero elements in the vector
The SS_IFLOMV and SS_IPFLOMV algorithms' steps are similar to those of SS_FLOMV and SS_PFLOMV algorithms except that Step 2 is applying the elimination operation (29) on the vector
In this section, a series of numerical experiments under different conditions are conducted to compare the performances of the proposed SS_FLOMV, SS_PFLOMV, SS_IFLOMV and SS_IPFLOMV algorithms with those of the FLOM_MUSIC, PFLOM_MUSIC and SS_SOSCV methods. Throughout this section, the convex optimization problems of (12), (17), (25) and (37) are resolved by using the software package CVX [24], the probability
$ { RMSE}=\frac{1}{N_{ok}}\sum\limits_{n=1}^{N_{ok}}\sqrt{\frac{1}{K} \sum\limits_{i=1}^K\left(\bar{\theta}_i(n)\theta_i(n)\right)^2} . $  (38) 
As the characteristic of the
$ {GSNR}=10\lg\frac{\delta_s^2}{\gamma} $  (39) 
where
Example 1: Three sources impinging on array for
Download:


Fig. 1 Normalized spatial spectrum of FLOM_MUSIC algorithm. 
Download:


Fig. 2 Normalized spatial spectrum of PFLOM_MUSIC algorithm. 
Download:


Fig. 3 Normalized spatial spectrum of SS_SOSCV algorithm. 
Download:


Fig. 4 Normalized spatial spectrum of SS_FLOMV algorithm. 
Download:


Fig. 5 Normalized spatial spectrum of SS_PFLOMV algorithm. 
Download:


Fig. 6 Normalized spatial spectrum of SS_IFLOMV algorithm. 
Download:


Fig. 7 Normalized spatial spectrum of SS_IPFLOMV algorithm. 
Example 2: Three sources impinging on array for
Download:


Fig. 8 Probability of resolution versus GSNR. 
Download:


Fig. 9 RMSE of DOA estimation versus GSNR. 
Example 3: Three sources impinging on array for
Download:


Fig. 10 Probability of resolution versus number of snapshots. 
Download:


Fig. 11 RMSE versus number of snapshots. 
Example 4: In this example, the performances of the proposed algorithms versus the characteristic exponent
Download:


Fig. 12 Probability of resolution versus characteristic exponent 
Download:


Fig. 13 RMSE versus characteristic exponent 
Although the FLOM and PFLOM have the good performances in suppressing the
Download:


Fig. 14 Probability of resolution versus characteristic exponent 
Download:


Fig. 15 RMSE versus characteristic exponent 
The new methods based on sparse representation of the fractional lower order statistics vector are proposed for the DOAs estimation under
[1]  H. C. So, "Source localization: algorithms and analysis, " in Handbook of Position Location: Theory, Practice, and Advances, S. A. Zekavat and R. M. Buehrer, Eds. New York, NY, USA: WileyIEEE Press, 2011. http://onlinelibrary.wiley.com/doi/10.1002/9781118104750.ch2/pdf 
[2]  H. Krim and M. Viberg, "Two decades of array signal processing research: the parametric approach, " IEEE Signal Process. Mag., vol. 13, no. 4, pp. 6794, Jul. 1996. http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=526899 
[3]  J. J. Fuchs, "On the use of the global matched filter for DOA estimation in the presence of correlated waveforms, " in Proc. of the 42nd Asilomar Conf. on Signals, Systems and Computers, Pacific Grove, CA, USA, 2008, pp. 269273. http://ieeexplore.ieee.org/document/5074406/ 
[4]  J. J. Fuchs, "Identification of real sinusoids in noise, the global matched filter approach, " in Proc. of the 15th IFCA Symp. System Identification, SaintMalo, France, 2009, pp. 11271132. http://dx.doi.org/10.3182/200907063FR2004.00187 
[5]  D. Malioutov, M. Cetin, and A. S. Willsky, "A sparse signal reconstruction perspective for source localization with sensor arrays, " IEEE Trans. Signal Process., vol. 53, no. 8, pp. 30103022, Aug. 2005. https://academic.oup.com/bioinformatics 
[6]  X. Xu, X. Wei, and Z. Ye, "DOA estimation based on sparse signal recovery utilizing weighted l_{1}norm penalty, " IEEE Signal Process. Lett., vol. 19, no. 3, pp. 155158, Mar. 2012. http://ieeexplore.ieee.org/document/6129390/ 
[7]  N. Hu, Z. F. Ye, D. Y. Xu, and S. H. Cao, "A sparse recovery algorithm for DOA estimation using weighted subspace fitting, " Signal Process., vol. 92, no. 10, pp. 25662570, Oct. 2012. http://www.sciencedirect.com/science/article/pii/S0165168412001119 
[8]  X. H. Wei, Y. B. Yuan, and Q. Ling, "DOA estimation using a greedy block coordinate descent algorithm, " IEEE Trans. Signal Process., vol. 60, no. 12, pp. 63826394, Dec. 2012. http://dl.acm.org/citation.cfm?id=2710152.2710644&coll=DL&dl=GUIDE 
[9]  M. M. Hyder and K. Mahata, "Direction of arrival estimation using a mixed l_{2, 0} norm approximation, " IEEE Trans. Signal Process., vol. 58, no. 9, pp. 46464655. http://ieeexplore.ieee.org/document/5466152/ 
[10]  P. Stoica, P. Babu, and J. Li, "SPICE: a sparse covariancebased estimation method for array processing, " IEEE Trans. Signal Process., vol. 59, no. 12, pp. 629638, Feb. 2011. http://www.ams.org/mathscinetgetitem?mr=2790596 
[11]  J. H. Yin and T. Q. Chen, "Directionofarrival estimation using a sparse representation of array covariance vectors, " IEEE Trans. Signal Process., vol. 59, no. 92, pp. 44894493, Sep. 2011. http://dl.acm.org/citation.cfm?id=2333345 
[12]  Z. Q. He, Q. H. Liu, L. N. Jin, and S. Ouyang, "Low complexity method for DOA estimation using array covariance matrix sparse representation, " Electron. Lett., vol. 49, no. 3, pp. 228230, Jan. 2013. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6457588 
[13]  Z. Q. He, Z. P. Shi, and L. Huang, "Covariance sparsityaware DOA estimation for nonuniform noise, " Digital Signal Process., vol. 28, pp. 7581, May 2014. http://www.sciencedirect.com/science/article/pii/S1051200414000505 
[14]  Y. Tian, X. Y. Sun, and S. S. Zhao, "DOA and power estimation using a sparse representation of second order statistics vector and l_{0}norm approximation, " Signal Process., vol. 105, pp. 98108, Dec. 2014. http://www.sciencedirect.com/science/article/pii/S0165168414002321 
[15]  Y. Tian, X. Y. Sun, and S. S. Zhao, "Sparsereconstructionbased direction of arrival, polarisation and power estimation using a crossdipole array, " IET Radar, Sonar Navig., vol. 9, no. 6, pp. 727731, Jun. 2015. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=7127163 
[16]  Z. H. Sha, Z. M. Liu, Z. T. Huang, and Y. Y. Zhou, "Covariancebased directionofarrival estimation of wideband coherent chirp signals via sparse representation, " Sensors, vol. 13, no. 9, pp. 1149011497, Aug. 2013. http://pubmedcentralcanada.ca/pmcc/articles/PMC3821291/ 
[17]  A. Mahmood, M. Chitre, and M. A. Armand, "PSK communication with passband additive symmetric αstable noise, " IEEE Trans. Commun., vol. 60, no. 10, pp. 29903000, Oct. 2012. http://ieeexplore.ieee.org/document/6253211/ 
[18]  C. L. Nikias and M. Shoa, Signal Processing with AlphaStable Distributions and Applications. New York, NY, USA: Wiely, 1995. 
[19]  T. H. Liu and J. M. Mendel, "A subspacebased direction finding algorithm using fractional lower order statistics, " IEEE Trans. Signal Process., vol. 49, no. 8, pp. 16051613, Aug. 2001. http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=934131 
[20]  H. Belkacemi and S. Marcos, "Robust subspacebased algorithms for joint angle/Doppler estimation in nonGaussian clutter, " Signal Process., vol. 87, no. 7, pp. 15471558, Jul. 2007. http://dl.acm.org/citation.cfm?id=1233044 
[21]  G. Tzagkarakis, P. Tsakalides, and J. L. Starck, "Covariationbased subspaceaugmented MUSIC for joint sparse support recovery in impulsive environments, " Signal Process., vol. 93, no. 5, pp. 13651373, May 2013. http://www.sciencedirect.com/science/article/pii/S0165168412004161 
[22]  W. K. Ma, T. H. Hsieh, and C. Y. Chi, "DOA estimation of quasistationary signals with less sensors than sources and unknown spatial noise covariance: a KhatriRao subspace approach, " IEEE Trans. Signal Process., vol. 58, no. 4, pp. 21682180, Apr. 2010. http://ieeexplore.ieee.org/document/5290056/ 
[23]  B. Ottersten, P. Stoica, and R. Roy, "Covariance matching estimation techniques for array signal processing applications, " Digit. Signal Process., vol. 8, no. 3, pp. 185210. http://kth.divaportal.org/smash/record.jsf?pid=diva2:501168 
[24]  M. Grant, S. Boyd, and Y. Ye, CVX: MATLAB Software for Disciplined Convex Programming. CVX Research, Inc., 2008. 