# DOA Estimation for a Mixture of Uncorrelated and Coherent Sources Based on Hierarchical Sparse Bayesian Inference with a Gauss-Exp-Chi2 Prior

Authors: Pinjiao Zhao; Weijian Si; Guobing Hu; Liwei Wang
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

Direction of arrival (DOA) estimation algorithms based on sparse Bayesian inference (SBI) can effectively estimate coherent sources without recurring to extra decorrelation techniques, and their estimation performance is highly dependent on the selection of sparse prior. Specifically, the specified sparse prior is expected to concentrate its mass on the zero and distribute with heavy tails; otherwise, these algorithms may suffer from performance degradation. In this paper, we introduce a new sparse-encouraging prior, referred to as “Gauss-Exp-Chi2” prior, and develop an efficient DOA estimation algorithm for a mixture of uncorrelated and coherent sources under a hierarchical SBI framework. The Gauss-Exp-Chi2 prior distribution exhibits a sharp peak at the origin and heavy tails, and this property makes it an appropriate prior to encourage sparse solutions. A three-layer hierarchical sparse Bayesian model is established. Then, by exploiting variational Bayesian approximation, the model parameters are estimated by alternately updating until Kullback-Leibler (KL) divergence between the true posterior and the variational approximation becomes zero. By constructing the source power spectra with the estimated model parameters, the number and locations of the highest peaks are extracted to obtain source number and DOA estimates. In addition, some implementation details for algorithm optimization are discussed and the Cramér-Rao bound (CRB) of DOA estimation is derived. Simulation results demonstrate the effectiveness and favorable performance of the proposed algorithm as compared with the state-of-the-art sparse Bayesian algorithms.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation has been a crucial issue in various application areas involving radar, wireless communication, and navigation [1–3]. Multiple signal classification (MUSIC) [4] and estimation of signal parameter via rotational invariance technique (ESPRIT) [5], known as the two most classical subspace-based DOA estimation algorithms, have been proposed to resolve uncorrelated sources. However, in multipath propagation environments, sources from an identical target may undergo reflection from various surfaces, and thus, the received sources may be a mixture of uncorrelated and coherent sources. In such environments, the aforementioned subspace-based algorithms would suffer from serious performance deterioration or even fail [6].In order to solve the aforementioned problem, several preprocessing techniques are developed for decorrelation. In the related studies, these preprocessing techniques are mainly classified into two categories: spatial smoothing (SS) [7] and matrix reconstruction (MR) [8]. Theoretically, the SS technique is implemented by dividing the whole array into multiple subarrays to combat the rank deficiency, while the MR technique is performed by rearranging the rank-deficient matrix into Hankel matrix, Toeplitz matrix, or other matrices to restore the rank. By combing the SS or MR techniques with subspace-based algorithms, several DOA estimation algorithms have been developed to handle the scenarios where uncorrelated and coherent sources coexist [9–13]. More specifically, the oblique projection spatial smoothing (OPSS) [9] and moduli property spatial smoothing (MPSS) [10] are the typical SS-based algorithms, and the oblique projection Toeplitz matrix reconstruction (OPT-MR) [11], symmetric non-Toeplitz matrix reconstruction (SNT-MR) [12], and real-valued Hankel matrix reconstruction (RVH-MR) [13] are the popular MR-based algorithms.Unlike traditional subspace-based DOA estimation algorithms, the emerging sparse source reconstruction (SSR) algorithms [14–19], including matching pursuit (MP) algorithm [14], lp-norm optimization algorithms [15, 16], and sparse Bayesian inference (SBI) algorithms [17–19], provide a new perspective for DOA estimation. Since SSR-based algorithms realize DOA estimation via sparse source reconstruction, instead of calculating the covariance matrix, they can resolve the coherent sources directly without extra preprocessing techniques. Among these SSR-based algorithms, both the MP and the lp-norm optimization algorithms are restricted to solve an optimization problem and are relied on the point estimate, so that these algorithms ignore the uncertainty of the underlying source in the process of source reconstruction. By contrast, SBI-based algorithms specially consider the uncertainty of the underlying source and estimate the source via choosing an appropriate prior, which yield favorable reconstruction performance [20]. Moreover, they can provide good estimation performance in the case of low SNR or small number of snapshots [21]. In general, SBI-based algorithms first specify a sparsity-encouraging prior to the unknown source model and then the model parameters are estimated via Bayesian inference. Since the exact Bayesian inference is intractable, two mainstream approximation inference algorithms were presented to estimate the model parameters, in which one is evidence procedure [22] and the other is variational Bayesian approximation [23]. For the evidence procedure, some unknown hyperparameters with respect to the hierarchical prior model are estimated iteratively by maximizing the evidence. For the variational Bayesian approximation, the posterior distribution is approximated as the product of several tractable distributions, and the model parameters keep updating to minimize the Kullback-Leibler (KL) divergence between the true posterior and the variational approximation, which has attractive computational efficiency along with high estimation performance [24]. Note that these approximation inference algorithms operate under the premise that an appropriate sparse prior has been imposed on the source model for the purpose of encouraging sparse solutions. Many sparsity-encouraging prior models have been investigated in the SBI framework [25–27]. In [25], Bayesian compressed sensing (BCS) was proposed with a two-layer hierarchical Gaussian-inverse-gamma prior (or Student’s-t prior), where the first layer is a Gaussian probability density function (pdf) and the second layer is a gamma pdf. Babacan et al. [26] proposed an equivalent Laplace prior that is generated by a Gaussian prior and an Exponential prior. In [27], a normal product (NP) prior was developed with two algorithms: NP-0 (using one-layer source model) and NP-1 (using two-layer source model). However, these priors are concentrated near the origin with relatively light tails, which may cause overshrinkage of the incident sources [28].In this paper, we develop a new sparse-encouraging prior (called Gauss-Exp-Chi2 prior) whose pdf distribution exhibits a sharp peak at the origin and heavy tails. With the proposed prior, the DOA estimation for a mixture of uncorrelated and coherent sources is performed under the hierarchical SBI framework using a uniform linear array (ULA). A three-layer hierarchical Bayesian model is established based on the Gauss-Exp-Chi2 prior. Subsequently, according to the variational Bayesian approximation, the model parameters (including the mean and variance of sparse sources and hyperparameters) keep alternately updating until the KL divergence between the true posterior and the variational approximation tends to be zero. By exploiting the estimated model parameters, the source power spectra is constructed, from which the number and locations of the highest peaks are extracted to obtain source number and DOA estimates. Simulation results show that the proposed algorithm has superior estimation performance. Now we briefly summarize the contributions of this work as follows:
(i)
To encourage sparse solutions, we develop a new sparse-encouraging prior, called Gauss-Exp-Chi2 prior, whose pdf distribution has a sharp peak at the origin and heavy tails.
(ii)
By constructing the source powers of all the potential directions in the angular space, both source number and DOA estimates are obtained.
(iii)
Variational approximations are adopted for the estimation of the hierarchical sparse Bayesian model parameters.
(iv)
Several implementation details for algorithm optimization including Woodbury matrix identity for dimension-reduction, pruning a basis function and the third kind Bessel function approximation are discussed, and the CRB of DOA estimation is derived.The remainder of this paper is organized as follows. The DOA estimation model for mixed sources is formulated in Section2. Section 3 presents the Gauss-Exp-Chi2 prior and DOA estimation algorithm for a mixture of uncorrelated and coherent sources within the hierarchical SBI framework. The algorithm optimization and CRB of DOA estimation are discussed in Section 4. Section 5 presents the simulation results of the proposed algorithm. Conclusions are drawn in Section 6.Notations.
Vectors and matrices are denoted by lowercase and uppercase boldface letters, respectively.⋅T, ⋅H, ⋅−1, and ⋅ represent transpose, conjugate transpose, inverse, and the statistical expectation, respectively. I is an identity matrix.⋅2 denotes the Euclidean norm, and diag⋅ denotes a diagonal matrix. Additionally, ∫abvdv represents the integral of v from a to b.


## 2. Problem Formulation

Consider a total ofK far-field narrowband sources impinging on the uniform linear array (ULA) consisting of M omnidirectional sensors with the interspacing between adjacent sensors d being a half of the carrier wavelength λ, that is, d=λ/2. Then, the array output vector at the tth snapshot can be expressed as
(1)xt=∑k=1Kaθkskt+nt,t=1,2,…,T,where T is the number of snapshots, aθk=1,e−jπsinθk,…,e−jπM−1sinθkT is the steering vector corresponding to the direction θk of the kth impinging source skt, and nt denotes the noise vector. Without loss of generality, the impinging source sktk=1K is a mixture of uncorrelated and coherent sources. Note that there exists power fading among the coherent source, and fading coefficients are used to describe the degree of power fading [5]. Specifically, consider K far-field narrowband sources sktk=1K impinging on the ULA, in which the first Ku sources s1t,s2t,…,sKut are uncorrelated, and the last Kcsources sKu+1t,sKu+2t,…,sKtare coherent. Then, the uncorrelated source set is denoted as sktk=1Ku and the coherent source set is denoted as sktk=Ku+1K=1,ρ1,ρ2,…,ρKc−1TsKu+1t with ρ1,ρ2,…,ρKc−1are the fading coefficients.Divide the entire angular space intoJ sampling grids θ⌢=θ⌢1,θ⌢1,…,θ⌢JT, where J represents the grid number and generally satisfies J≫M>K. Assume that θ⌢jk∈θ⌢1,θ⌢1,…,θ⌢JT is the nearest grid to true direction θj∈θ1,θ2,…,θKT; thus, we have hjt=sjktδj−jk. Thus, xt can be rewritten in a sparse form
(2)xt=A⌢ht+nt,t=1,2,…,T,where A⌢=αθ⌢1,αθ⌢2,…,αθ⌢J and ht=h1t,h2t,…,hJtT. Due to the fact that ht has K non-zero elements in J elements, it is a sparse vector. In the case of multiple snapshots, the sparse sources at all the snapshots share the same support, and the array output matrix of T snapshots can be represented by
(3)X=A⌢H+N,where X=x1,x2,…,xT, H=h1,h2,…,hT, and N=n1n2…nT. The goal of this paper is to provide the DOA estimation under the coexistence of uncorrelated and coherent sources from a sparse Bayesian perspective.


## 3. DOA Estimation

In this section, a DOA estimation algorithm for a mixture of uncorrelated and coherent sources is proposed within the hierarchical SBI framework. A Gauss-Exp-Chi2 prior is developed to encourage sparse solutions, and then the parameters of three-layer hierarchical Bayesian model are estimated via variational Bayesian approximations. By constructing source power spectra, the source number and DOA estimation are obtained.


### 3.1. Bayesian Model

In the Bayesian model, the pdf of a joint distributionpX,H,α,η with respect to all the unknown and observed quantities is required to be known. In this paper, the joint distribution can be expressed as
(4)pX,H,α,η=pXH,σ2pHαpατ,ηpηv,where α and η are referred to as the hyperparameters; τ≥0 and v≥0 are, respectively, referred to as the rate parameter and shape parameter. It is assumed that the components nt,t=1,2,…,T are the independently zero-mean stationary Gaussian noise with known variance σ2. Thus, the pdf of the noise matrix N is given by
(5)pNσ2=∏t=1TNnt0,σ2I.Combining (3) and (5), the Gaussian likelihood model is obtained as follows:
(6)pXH,σ2=∏t=1TNxtA⌢ht,σ2I.From a Bayesian perspective, the pdf distribution of an assigned prior is appealing to exhibit a sharp peak at the origin and heavy tails, which favors strong shrinkage of noise sources and avoids overshrinkage of the interest sources. This property is generally considered as a desirable property for enforcing sparsity and variable selection [24]. Some typical sparse priors, such as Gaussian-inverse-gamma prior and Laplace prior, are widely used in the relevant research [25, 26], which, however, are concentrated near the origin with relatively light tails [28]. To alleviate this problem, we here develop a three-layer hierarchical prior, referred to as Gauss-Exp-Chi2 prior, for H. In the first layer of prior, we adopt a zero-mean Gaussian prior
(7)pHα=∏t=1TNht0,Λ−1=∏t=1T∏j=1JNhjt0,αj,with being Λ−1=diagα=diagα1,α2,…,αJ. In the second layer of prior, an exponential hyperprior is imposed on α
(8)pατ,η=∏j=1JExpαjτη,where Exp⋅ denotes the exponential distribution. In the third layer of prior, a chi-square (Chi2) hyperprior is considered over η, that is,
(9)pηv=χ2ηv=Γv2−12−v/2ηv/2−1exp−η2,where Γ⋅ denotes the gamma function.The first two layers (7) and (8) of the proposed prior result in a generative Laplace (Gaussian-Exponential) distribution [26], and the last layer (i.e., a Chi2 distribution) is embedded to obtain the proposed three-layer Gauss-Exp-Chi2 prior. Thus, the proposed prior has more free parameters to control the degree of sparseness as compared with the Laplace prior [29].Based on the above analysis, the directed graph of the sparse Bayesian model is shown in Figure1, where arrows are used to point to the generative model. Note that the first five blocks from the left (corresponding to the variables v, η, τ, α1,…,αJ, and h1t,…,hJt) depict the three-layer hierarchical prior mentioned above.Figure 1
Directed graph of the sparse Bayesian model.Sincepαj∣τv∝ταj+1/2−v/2+1,j=1,2,…,J, it is obvious that most αj for j=1,2,…,J are favored to zero, which leads to the fact that most hjt are favored to zero. Thus, the proposed three-layer hierarchical Gauss-Exp-Chi2 prior is a sparsity-encouraging prior, which can help to improve the performance of source reconstruction [30]. By combining the three layers of the proposed prior, the generated prior can be computed via
(10)pH;τ,v=∏t=1T∏j=1Jphjt;τ,v,where
(11)phjt;τ,v=∫0+∞∫0+∞phjtαjpαjτ,ηpηvdαjdη.Note thatpαj∣τv=∫0+∞pαj∣τηpη∣vdη; thus, phjt;τ,v can be computed by marginalizing phjt∣αjpαj∣τv over αj (see Appendix for detail).
(12)pH;τ,v=∏t=1T∏j=1J−τπΓv/2+1Γv/2+1/2Γv/2Uchfv2+12,12,hj2tτ,where Uchf⋅,⋅,⋅is the confluent hypergeometric function defined by
(13)Uchfa,a−b+1,c=1Γa∫0+∞ta−11+t−bexp−ctdt.The reason for choosing the Gauss-Exp-Chi2 prior is twofold: its pdf distribution has a sharp spike at the origin and heavy tails; it forms in a conjugate manner since the exponential and Chi2 distributions belonging to the exponential distribution family are chosen as the 2nd layer and the 3rd layer of this prior, which significantly simplifies the form of posterior distributions [24]. The tail to spike behavior for Gauss-Exp-Chi2 prior, Laplace prior, Student’s-t prior, and Gaussian prior is illustrated in Figure 2, wherein the parameter selection of these priors is according to the standard derivation [25, 26]. It is observed that the proposed Gauss-Exp-Chi2 prior is a sparsity-encouraging prior which has sharper peak and heavier tails than the existing priors.Figure 2
Four kinds of pdf curves with the standard derivation: (a) compares four distributions at origin and (b) compares four distributions in tail.
(a)
(b)


### 3.2. Variational Bayesian Approximations

It is well known that Bayesian inference operates on the basis of the posterior distributionpH,α,η∣X=pH,α,η,X/pX. However, this posterior distribution is intractable for the reason that pX=∫∫∫pH,α,η,XdHdαdη cannot be calculated analytically. In order to approximate this posterior distribution, we resort to mean-field variational Bayesian approximation [23], whose task is to seek several analytically approximate distributions that minimize the KL divergence between the posterior pH,α,η∣X and its approximation distributionqH,α,η. According to [23], qH,α,η can be factorized as the product of three variational distributions qH, qα, and qη; thus, we have
(14)pH,α,ηX≈qH,α,η=qHqαqη.Denote a setζ=H,α,η, the distribution of each variable ζk is expressed as lnqζk∝<lnqζk,X>qζ/ζk, where ζ/ζk represents the subset excluding ζk. More specifically, the posterior approximation of H, α, and η in this paper are, respectively, given by
(15)lnqH∝<lnpX,H,α,η>qα,(16)lnqα∝<lnpX,H,α,η>qHqη,(17)lnqη∝<lnpX,H,α,η>qα.(1)Update of qH. Keeping the terms of qH only depend on H, we have
(18)lnqH∝<lnpXHpHα>qα∝∑t=1Tσ−2hHtΑ⌢HΑ⌢xt−12σ−2hHtΛHht,with Λ=diag1/α1,1/α2,…,1/αJ. By gathering together the similar terms, we have
(19)lnqH∝∑t=1T−12hHtσ−2Α⌢HΑ⌢+<Λ>ht+σ−2∑t=1ThHtΑ⌢Hxt⇒qH∝∑t=1Texp−12hHtΣ−1ht+hHtΣ−1μt,which implies that qH is the product of multiple Gaussian distributions Nht∣μtΣ with the mean μt and covariance Σ being
(20)μt=σ−2Σ−1Α⌢Hxt,(21)Σ=σ−2Α⌢HΑ⌢+<Λ>−1.(2)Update of qα. Similarly, according to (16), qα is derived as
(22)lnqα∝<lnpHαpατ,η>qHqη∝<ln∏t=1T∏j=1J12παjexp−hj2t2αj∏j=1Jτηexp−τηαj>qHqη.Denotinghj2=∑t=1Thj2t and <hj2>=μj2+Σjj, (22) can be further simplified to
(23)lnqα∝∑j=1J<lnαj−1/2exp−hj22αj−τηαj>qHqη⇒qα∝∏j=1Jαj−1/2exp−<hj2>2αj−τ<η>αj.The posterior distribution ofqα can be approximately equivalent to the product of a series of generalized inverse Gaussian (GIG) distributions, that is, α∼∏j=1JGIGαj∣1/2<hj2>2τ<η>, and thus, we have
(24)<αjn>=<hj2>2τ<η>n/2κn−1/22τ<η><hj2>κ−1/22τ<η><hj2>,where κp⋅ is the third kind Bessel function with order p. The case of n=1 in (24) gives the calculation of <αj>, which can be used in (25). The case of n=−1 gives <αj−1> used for the evaluation of Σ in (21).(3)Update of qη. For qη, we have
(25)lnqη∝<lnpατ,ηpηv>qα∝<ln∏j=1Jτηexp−τηαj2−v/2Γv/2ηv/2−1exp−η2>qα⇒qη∝ηJ+v/2−1exp−τ∑j=1J<αj>η−η2.Thus, it can be known thatqηis a gamma distribution, that is, η∼GammaJ+v/2,τ∑j=1J<αj>+1/2 with the following mean
(26)<η>=J+v/2τ∑j=1J<αj>+1/2.The estimation of model parameters is implemented by alternately updating the meanμ, variance Σ, and hyperparameters <αj> and <η> to minimize the KL divergence between the true posterior and the variational approximation, and the main steps are summarized in Algorithm 1.Algorithm 1: The main steps of the model parameter estimation.
Input: Α⌢, X
Output:μ, Σ
1.Initialize: <αj>, <αj−1>, <η>, τ, v, σ2, ε, imax
2.while <αj>i−<αj>i−1/<αj>i−1>ε and i<imax do
3.  updateΣi using (21)
4.  updateμti using (20)
5.  update<αj>i and <αj−1>i using (24)
6.  update<η>i using (26)
7.i=i+1
8.end while


### 3.3. Source Power and DOA Estimation

The DOAs of the impinging sources are estimated via the following two steps: (1) form spatial spectra using the estimated source powers of all the potential directions and (2) extract the number and the corresponding locations of the peaks beyond a power threshold from the spectra. LetU^=μ1,μ2,…,μT and Σ^ be the final estimates of U and Σ (the mean and variance with respect to H) and consider H row by row, then we have hj∼NU^j,Σ^jjI. As outlined in [18], let Poj be the source power of direction θ⌢j; then we have
(27)P^oj=1TU^j22+TΣ^jj=1TU^j22+Σ^jj.Therefore, the source powers of all the potential directions in the angular space are estimated by calculatingP^o=P^o1,P^o2,…,P^oJT. Suppose that there are K peaks exceeding the power threshold Pthres and the corresponding grid indices are jkjk=1,2,…,K. Then, the estimated source number and DOAs will be K and θ⌢jkjk=1,2,…,K.


## 3.1. Bayesian Model

In the Bayesian model, the pdf of a joint distributionpX,H,α,η with respect to all the unknown and observed quantities is required to be known. In this paper, the joint distribution can be expressed as
(4)pX,H,α,η=pXH,σ2pHαpατ,ηpηv,where α and η are referred to as the hyperparameters; τ≥0 and v≥0 are, respectively, referred to as the rate parameter and shape parameter. It is assumed that the components nt,t=1,2,…,T are the independently zero-mean stationary Gaussian noise with known variance σ2. Thus, the pdf of the noise matrix N is given by
(5)pNσ2=∏t=1TNnt0,σ2I.Combining (3) and (5), the Gaussian likelihood model is obtained as follows:
(6)pXH,σ2=∏t=1TNxtA⌢ht,σ2I.From a Bayesian perspective, the pdf distribution of an assigned prior is appealing to exhibit a sharp peak at the origin and heavy tails, which favors strong shrinkage of noise sources and avoids overshrinkage of the interest sources. This property is generally considered as a desirable property for enforcing sparsity and variable selection [24]. Some typical sparse priors, such as Gaussian-inverse-gamma prior and Laplace prior, are widely used in the relevant research [25, 26], which, however, are concentrated near the origin with relatively light tails [28]. To alleviate this problem, we here develop a three-layer hierarchical prior, referred to as Gauss-Exp-Chi2 prior, for H. In the first layer of prior, we adopt a zero-mean Gaussian prior
(7)pHα=∏t=1TNht0,Λ−1=∏t=1T∏j=1JNhjt0,αj,with being Λ−1=diagα=diagα1,α2,…,αJ. In the second layer of prior, an exponential hyperprior is imposed on α
(8)pατ,η=∏j=1JExpαjτη,where Exp⋅ denotes the exponential distribution. In the third layer of prior, a chi-square (Chi2) hyperprior is considered over η, that is,
(9)pηv=χ2ηv=Γv2−12−v/2ηv/2−1exp−η2,where Γ⋅ denotes the gamma function.The first two layers (7) and (8) of the proposed prior result in a generative Laplace (Gaussian-Exponential) distribution [26], and the last layer (i.e., a Chi2 distribution) is embedded to obtain the proposed three-layer Gauss-Exp-Chi2 prior. Thus, the proposed prior has more free parameters to control the degree of sparseness as compared with the Laplace prior [29].Based on the above analysis, the directed graph of the sparse Bayesian model is shown in Figure1, where arrows are used to point to the generative model. Note that the first five blocks from the left (corresponding to the variables v, η, τ, α1,…,αJ, and h1t,…,hJt) depict the three-layer hierarchical prior mentioned above.Figure 1
Directed graph of the sparse Bayesian model.Sincepαj∣τv∝ταj+1/2−v/2+1,j=1,2,…,J, it is obvious that most αj for j=1,2,…,J are favored to zero, which leads to the fact that most hjt are favored to zero. Thus, the proposed three-layer hierarchical Gauss-Exp-Chi2 prior is a sparsity-encouraging prior, which can help to improve the performance of source reconstruction [30]. By combining the three layers of the proposed prior, the generated prior can be computed via
(10)pH;τ,v=∏t=1T∏j=1Jphjt;τ,v,where
(11)phjt;τ,v=∫0+∞∫0+∞phjtαjpαjτ,ηpηvdαjdη.Note thatpαj∣τv=∫0+∞pαj∣τηpη∣vdη; thus, phjt;τ,v can be computed by marginalizing phjt∣αjpαj∣τv over αj (see Appendix for detail).
(12)pH;τ,v=∏t=1T∏j=1J−τπΓv/2+1Γv/2+1/2Γv/2Uchfv2+12,12,hj2tτ,where Uchf⋅,⋅,⋅is the confluent hypergeometric function defined by
(13)Uchfa,a−b+1,c=1Γa∫0+∞ta−11+t−bexp−ctdt.The reason for choosing the Gauss-Exp-Chi2 prior is twofold: its pdf distribution has a sharp spike at the origin and heavy tails; it forms in a conjugate manner since the exponential and Chi2 distributions belonging to the exponential distribution family are chosen as the 2nd layer and the 3rd layer of this prior, which significantly simplifies the form of posterior distributions [24]. The tail to spike behavior for Gauss-Exp-Chi2 prior, Laplace prior, Student’s-t prior, and Gaussian prior is illustrated in Figure 2, wherein the parameter selection of these priors is according to the standard derivation [25, 26]. It is observed that the proposed Gauss-Exp-Chi2 prior is a sparsity-encouraging prior which has sharper peak and heavier tails than the existing priors.Figure 2
Four kinds of pdf curves with the standard derivation: (a) compares four distributions at origin and (b) compares four distributions in tail.
(a)
(b)


## 3.2. Variational Bayesian Approximations

It is well known that Bayesian inference operates on the basis of the posterior distributionpH,α,η∣X=pH,α,η,X/pX. However, this posterior distribution is intractable for the reason that pX=∫∫∫pH,α,η,XdHdαdη cannot be calculated analytically. In order to approximate this posterior distribution, we resort to mean-field variational Bayesian approximation [23], whose task is to seek several analytically approximate distributions that minimize the KL divergence between the posterior pH,α,η∣X and its approximation distributionqH,α,η. According to [23], qH,α,η can be factorized as the product of three variational distributions qH, qα, and qη; thus, we have
(14)pH,α,ηX≈qH,α,η=qHqαqη.Denote a setζ=H,α,η, the distribution of each variable ζk is expressed as lnqζk∝<lnqζk,X>qζ/ζk, where ζ/ζk represents the subset excluding ζk. More specifically, the posterior approximation of H, α, and η in this paper are, respectively, given by
(15)lnqH∝<lnpX,H,α,η>qα,(16)lnqα∝<lnpX,H,α,η>qHqη,(17)lnqη∝<lnpX,H,α,η>qα.(1)Update of qH. Keeping the terms of qH only depend on H, we have
(18)lnqH∝<lnpXHpHα>qα∝∑t=1Tσ−2hHtΑ⌢HΑ⌢xt−12σ−2hHtΛHht,with Λ=diag1/α1,1/α2,…,1/αJ. By gathering together the similar terms, we have
(19)lnqH∝∑t=1T−12hHtσ−2Α⌢HΑ⌢+<Λ>ht+σ−2∑t=1ThHtΑ⌢Hxt⇒qH∝∑t=1Texp−12hHtΣ−1ht+hHtΣ−1μt,which implies that qH is the product of multiple Gaussian distributions Nht∣μtΣ with the mean μt and covariance Σ being
(20)μt=σ−2Σ−1Α⌢Hxt,(21)Σ=σ−2Α⌢HΑ⌢+<Λ>−1.(2)Update of qα. Similarly, according to (16), qα is derived as
(22)lnqα∝<lnpHαpατ,η>qHqη∝<ln∏t=1T∏j=1J12παjexp−hj2t2αj∏j=1Jτηexp−τηαj>qHqη.Denotinghj2=∑t=1Thj2t and <hj2>=μj2+Σjj, (22) can be further simplified to
(23)lnqα∝∑j=1J<lnαj−1/2exp−hj22αj−τηαj>qHqη⇒qα∝∏j=1Jαj−1/2exp−<hj2>2αj−τ<η>αj.The posterior distribution ofqα can be approximately equivalent to the product of a series of generalized inverse Gaussian (GIG) distributions, that is, α∼∏j=1JGIGαj∣1/2<hj2>2τ<η>, and thus, we have
(24)<αjn>=<hj2>2τ<η>n/2κn−1/22τ<η><hj2>κ−1/22τ<η><hj2>,where κp⋅ is the third kind Bessel function with order p. The case of n=1 in (24) gives the calculation of <αj>, which can be used in (25). The case of n=−1 gives <αj−1> used for the evaluation of Σ in (21).(3)Update of qη. For qη, we have
(25)lnqη∝<lnpατ,ηpηv>qα∝<ln∏j=1Jτηexp−τηαj2−v/2Γv/2ηv/2−1exp−η2>qα⇒qη∝ηJ+v/2−1exp−τ∑j=1J<αj>η−η2.Thus, it can be known thatqηis a gamma distribution, that is, η∼GammaJ+v/2,τ∑j=1J<αj>+1/2 with the following mean
(26)<η>=J+v/2τ∑j=1J<αj>+1/2.The estimation of model parameters is implemented by alternately updating the meanμ, variance Σ, and hyperparameters <αj> and <η> to minimize the KL divergence between the true posterior and the variational approximation, and the main steps are summarized in Algorithm 1.Algorithm 1: The main steps of the model parameter estimation.
Input: Α⌢, X
Output:μ, Σ
1.Initialize: <αj>, <αj−1>, <η>, τ, v, σ2, ε, imax
2.while <αj>i−<αj>i−1/<αj>i−1>ε and i<imax do
3.  updateΣi using (21)
4.  updateμti using (20)
5.  update<αj>i and <αj−1>i using (24)
6.  update<η>i using (26)
7.i=i+1
8.end while


## 3.3. Source Power and DOA Estimation

The DOAs of the impinging sources are estimated via the following two steps: (1) form spatial spectra using the estimated source powers of all the potential directions and (2) extract the number and the corresponding locations of the peaks beyond a power threshold from the spectra. LetU^=μ1,μ2,…,μT and Σ^ be the final estimates of U and Σ (the mean and variance with respect to H) and consider H row by row, then we have hj∼NU^j,Σ^jjI. As outlined in [18], let Poj be the source power of direction θ⌢j; then we have
(27)P^oj=1TU^j22+TΣ^jj=1TU^j22+Σ^jj.Therefore, the source powers of all the potential directions in the angular space are estimated by calculatingP^o=P^o1,P^o2,…,P^oJT. Suppose that there are K peaks exceeding the power threshold Pthres and the corresponding grid indices are jkjk=1,2,…,K. Then, the estimated source number and DOAs will be K and θ⌢jkjk=1,2,…,K.


## 4. Discussions

### 4.1. Algorithm Optimization

For alleviating the computational complexity and speeding up the update efficiency, some implementation details are adopted for algorithm optimization in this subsection.


#### 4.1.1. Woodbury Matrix Identity for Dimension-Reduction

At each iteration, it is inevitable to calculate aJ×J matrix inversion when updating Σ according to (21), which requires OJ3 computations. To reduce the computational complexity, we adopt the Woodbury matrix identity, and then (21) can be rewritten as Σ=<Λ−1>−<Λ−1>A⌢HG−1A⌢<Λ−1> with an M×M matrix G=σ−2I+A⌢<Λ−1>A⌢H. Thus, the computational complexity is reduced to OM3.


#### 4.1.2. Pruning a Basis Function

In order to speed up the update efficiency,<αj> and the corresponding aθ⌢j would be pruned from the original sparse Bayesian model when <αj> is smaller than a certain threshold. This procedure is referred to as “pruning a basis function”, and similar approaches have been used in [29, 30].


#### 4.1.3. The Third Kind Bessel Function Approximation

It is known that (24) involves the third kind Bessel function, and this integral process is relatively complicated. According to [31], the third kind Bessel function of m can be approximated as κpm∼1/2Γpm/2−p, if p>0, and κpm∼1/2Γ−pm/2p, if p<0. Then, (24) can be simplified to
(28)<αjn>=<hj2>2τ<η>1/2,n=1,<hj2>2τ<η>−1/22Γ3/2Γ1/212τ<η><hj2>,n=−1.


### 4.2. Cramér-Rao Bound (CRB)

To evaluate the DOA estimation performance, the CRB of the DOA estimation for a mixture of uncorrelated and coherent sources is derived in this subsection. Consider the sparse form of the array output vectorxt mentioned in (2), the log-likelihood function of x1,x2,…,xT is firstly construct as
(29)lnℒx1,x2,…,xT=const−MTlnσ−σ−1∑t=1Tx∗t−h∗tA⌢∗xt−A⌢ht.Lethrt and hit be the real part and imaginary part of ht, that is, hrt≜Reht and hit≜Imht. Then, the partial derivatives of (29) with respect to hrt, hit, σ, and θ⌢ are, respectively, given as
(30)∂lnℒ∂hrt=2σ−1ReA⌢∗nt,∂lnℒ∂hit=2σ−1ImA⌢∗nt,∂lnℒ∂σ=MTσ−1+σ−2∑t=1Tn∗tnt,∂lnℒ∂θ⌢=2σ−1∑t=1TReHd∗tAD∗nt,where Hd∗t=diagh1t,h2t,…,hJt and AD=daθ⌢1/dθ⌢1,daθ⌢2/dθ⌢2,…,daθ⌢J/dθ⌢J. Subsequently, define
(31)ℵ=E∂lnℒ∂θ⌢∂lnℒ∂θ⌢T=2σ−1∑t=1TReHd∗tAD∗ADHdt,Ξrm=ReΞm=E∂lnℒ∂hrt∂lnℒ∂θ⌢T=2σ−1ReA⌢∗ADHdt,Ξim=ImΞm=E∂lnℒ∂hit∂lnℒ∂θ⌢T=2σ−1ImA⌢∗ADHdt,Qr=E∂lnℒ∂hrt∂lnℒ∂hrkT−1=E∂lnℒ∂hit∂lnℒ∂hikT−1=2σ−1ReA⌢∗A⌢δkt−1,Qi=−E∂lnℒ∂hrt∂lnℒ∂hikT−1=2σ−1ImA⌢∗A⌢δkt−1.Then, the Fisher information matrix (FIM) [32] can be obtained as
(32)FIM=ℵ−Ξr1Ξi1Ξr2Ξi2⋯ΞrMΞiM=Qr−Qi0⋯0QiQr−Qi⋱⋮0Qi⋱⋱0⋮⋱⋱Qr−Qi0⋯0QiQrΞr1Ξi1Ξr2Ξi2⋮ΞrMΞiM.It is noted that the following equality holds(33)ΞrTΞiTQr−QiQiQrΞrΞi=ReΞ∗QΞ,with Q=2σ−1A⌢∗A⌢−1 and Ξ=2σ−1A⌢∗ADHdt. Finally, by combing (32) and (33), we have
(34)CRBθ⌢=FIM−1=0.5σ∑t=1TReHdtAD∗I−A⌢A⌢∗A⌢−1A⌢∗ADHdt−1.


## 4.1. Algorithm Optimization

For alleviating the computational complexity and speeding up the update efficiency, some implementation details are adopted for algorithm optimization in this subsection.


### 4.1.1. Woodbury Matrix Identity for Dimension-Reduction

At each iteration, it is inevitable to calculate aJ×J matrix inversion when updating Σ according to (21), which requires OJ3 computations. To reduce the computational complexity, we adopt the Woodbury matrix identity, and then (21) can be rewritten as Σ=<Λ−1>−<Λ−1>A⌢HG−1A⌢<Λ−1> with an M×M matrix G=σ−2I+A⌢<Λ−1>A⌢H. Thus, the computational complexity is reduced to OM3.


### 4.1.2. Pruning a Basis Function

In order to speed up the update efficiency,<αj> and the corresponding aθ⌢j would be pruned from the original sparse Bayesian model when <αj> is smaller than a certain threshold. This procedure is referred to as “pruning a basis function”, and similar approaches have been used in [29, 30].


### 4.1.3. The Third Kind Bessel Function Approximation

It is known that (24) involves the third kind Bessel function, and this integral process is relatively complicated. According to [31], the third kind Bessel function of m can be approximated as κpm∼1/2Γpm/2−p, if p>0, and κpm∼1/2Γ−pm/2p, if p<0. Then, (24) can be simplified to
(28)<αjn>=<hj2>2τ<η>1/2,n=1,<hj2>2τ<η>−1/22Γ3/2Γ1/212τ<η><hj2>,n=−1.


## 4.1.1. Woodbury Matrix Identity for Dimension-Reduction

At each iteration, it is inevitable to calculate aJ×J matrix inversion when updating Σ according to (21), which requires OJ3 computations. To reduce the computational complexity, we adopt the Woodbury matrix identity, and then (21) can be rewritten as Σ=<Λ−1>−<Λ−1>A⌢HG−1A⌢<Λ−1> with an M×M matrix G=σ−2I+A⌢<Λ−1>A⌢H. Thus, the computational complexity is reduced to OM3.


## 4.1.2. Pruning a Basis Function

In order to speed up the update efficiency,<αj> and the corresponding aθ⌢j would be pruned from the original sparse Bayesian model when <αj> is smaller than a certain threshold. This procedure is referred to as “pruning a basis function”, and similar approaches have been used in [29, 30].


## 4.1.3. The Third Kind Bessel Function Approximation

It is known that (24) involves the third kind Bessel function, and this integral process is relatively complicated. According to [31], the third kind Bessel function of m can be approximated as κpm∼1/2Γpm/2−p, if p>0, and κpm∼1/2Γ−pm/2p, if p<0. Then, (24) can be simplified to
(28)<αjn>=<hj2>2τ<η>1/2,n=1,<hj2>2τ<η>−1/22Γ3/2Γ1/212τ<η><hj2>,n=−1.


## 4.2. Cramér-Rao Bound (CRB)

To evaluate the DOA estimation performance, the CRB of the DOA estimation for a mixture of uncorrelated and coherent sources is derived in this subsection. Consider the sparse form of the array output vectorxt mentioned in (2), the log-likelihood function of x1,x2,…,xT is firstly construct as
(29)lnℒx1,x2,…,xT=const−MTlnσ−σ−1∑t=1Tx∗t−h∗tA⌢∗xt−A⌢ht.Lethrt and hit be the real part and imaginary part of ht, that is, hrt≜Reht and hit≜Imht. Then, the partial derivatives of (29) with respect to hrt, hit, σ, and θ⌢ are, respectively, given as
(30)∂lnℒ∂hrt=2σ−1ReA⌢∗nt,∂lnℒ∂hit=2σ−1ImA⌢∗nt,∂lnℒ∂σ=MTσ−1+σ−2∑t=1Tn∗tnt,∂lnℒ∂θ⌢=2σ−1∑t=1TReHd∗tAD∗nt,where Hd∗t=diagh1t,h2t,…,hJt and AD=daθ⌢1/dθ⌢1,daθ⌢2/dθ⌢2,…,daθ⌢J/dθ⌢J. Subsequently, define
(31)ℵ=E∂lnℒ∂θ⌢∂lnℒ∂θ⌢T=2σ−1∑t=1TReHd∗tAD∗ADHdt,Ξrm=ReΞm=E∂lnℒ∂hrt∂lnℒ∂θ⌢T=2σ−1ReA⌢∗ADHdt,Ξim=ImΞm=E∂lnℒ∂hit∂lnℒ∂θ⌢T=2σ−1ImA⌢∗ADHdt,Qr=E∂lnℒ∂hrt∂lnℒ∂hrkT−1=E∂lnℒ∂hit∂lnℒ∂hikT−1=2σ−1ReA⌢∗A⌢δkt−1,Qi=−E∂lnℒ∂hrt∂lnℒ∂hikT−1=2σ−1ImA⌢∗A⌢δkt−1.Then, the Fisher information matrix (FIM) [32] can be obtained as
(32)FIM=ℵ−Ξr1Ξi1Ξr2Ξi2⋯ΞrMΞiM=Qr−Qi0⋯0QiQr−Qi⋱⋮0Qi⋱⋱0⋮⋱⋱Qr−Qi0⋯0QiQrΞr1Ξi1Ξr2Ξi2⋮ΞrMΞiM.It is noted that the following equality holds(33)ΞrTΞiTQr−QiQiQrΞrΞi=ReΞ∗QΞ,with Q=2σ−1A⌢∗A⌢−1 and Ξ=2σ−1A⌢∗ADHdt. Finally, by combing (32) and (33), we have
(34)CRBθ⌢=FIM−1=0.5σ∑t=1TReHdtAD∗I−A⌢A⌢∗A⌢−1A⌢∗ADHdt−1.


## 5. Simulation Results

In this section, several simulations are presented to illustrate the performance of the proposed algorithm as compared with the BCS [25] and NP-1 [27] algorithms. Consider a 9-sensor ULA with sensor interspacing being a half of the carrier wavelength. In the proposed algorithm, hyperparameters α and η are initialized as 1/T∑t=1TA⌢Hxt and 0.1; the values of rate parameter τ, shape parameter v, and power threshold Pthres are, respectively, set to be τ=1.5, v=1, and Pthres=−10dB. Two hundred independent Monte Carlo trials are conducted for the following simulations, and the success rate and root mean squared error (RMSE) are two significant performance metrics for evaluating the DOA estimation performance. The success rate is defined as the rate between the number of successful estimates and the total number of Monte Carlo trials, where a successful estimate is declared if the estimation error is within a certain small Euclidean distance of the true directions. The RMSE is defined by
(35)RMSE=1200K∑i=1200∑k=1Kθ⌢k,i−θk2,where θ⌢k,i is the estimates of θk in the ith Monte Carlo trial.In the first simulation, we evaluate the effectiveness of the proposed algorithm. Assume that two uncorrelated sources from−32.2°,−11.5° and two coherent sources from 6.7°,23.4° with fading coefficients 1,0.2891−0.7567j impinge on this ULA. The grid number, SNR, and number of snapshots are, respectively, set to be 361, 0 dB, and 200. Figure 3 plots the source power spectrum, which is a function of DOA.Figure 3
Spectra for the proposed algorithm with the fixed SNR 0 dB and number of snapshots 200.As can be seen from Figure3, the power peaks in the vicinity of the true DOAs can be distinguished from other peaks by their relatively strong powers and the biases between these peaks and the true DOAs are slight, which demonstrate the effectiveness of the proposed algorithm.The second simulation compares the estimation performance of three algorithms, including the proposed algorithm, BCS, and NP-1 algorithms. Consider one uncorrelated source from19.2° and two coherent sources from −25.5°,33.3° with fading coefficients being 1,−0.5280+0.6010j. The grid number and SNR are, respectively, fixed at 361 and 0 dB. The success rate versus number of snapshots and SNR are depicted in Figures 4 and 5, whereas the RMSE versus number of snapshots and SNR are shown in Figures 6 and 7. The results from Figures 4 and 5 illustrate that all the three algorithms are able to estimate the DOA correctly under the condition that the SNR or number of snapshots are relatively high, but for low SNR or small number of snapshots, the proposed algorithm has a higher success rate than the BCS and NP-1 algorithms. Figures 6 and 7 show that the estimation accuracy of the three algorithms tends to improve with the increase of SNR or the number of snapshots, and the proposed algorithm yields the best estimation accuracy as compared to the other two algorithms. This is because the proposed Gauss-Exp-Chi2 prior is a three-layer sparsity-encouraging prior, which can help to improve the accuracy of source reconstruction.Figure 4
Success rates versus number of snapshots with the fixed SNR 0 dB.Figure 5
Success rates versus SNR with the fixed number of snapshots 200.Figure 6
RMSE versus number of snapshots with the fixed SNR 0 dB.Figure 7
RMSE versus SNR with the fixed number of snapshots 200.The third simulation investigates the RMSE versus grid number. The simulation settings are the same as those of the second simulation, except that grid number in this simulation is ranged from 61 to 361. Figure8 plots the RMSE curves versus grid number with fixed SNR 10 dB and number of snapshots 500.Figure 8
RMSE versus grid number with the fixed SNR 10 dB and number of snapshots 500.It is observed that the RMSE curves of the three algorithms rapidly decrease as grid number increases from 61 to 181, and these curves tend to slowly decrease with the increase of grid number when grid number is beyond 181. In addition, the proposed algorithm has minimum estimation errors among the three algorithms. Note that the three algorithms recover sparse sources from a Bayesian perspective; thus, they are not restricted to the restricted isometry property. This implies that the grid number can continue to increase to further improve the estimation accuracy and the computational complexity would increase accordingly. Therefore, for the purpose of balancing estimation accuracy and computational complexity, the reasonable range of grid number is from 181 to 361.In the last simulation, we test the angular resolution of the three algorithms. Consider two coherent sources from 10.1° and10.1°+Δθ with the angular separation Δθ ranging from 3° to 10°. The fading coefficients are 1,0.4469−0.7696j, and the grid number is set to be 361. Figure 9 depicts the RMSE versus angular separation with the fixed SNR 10 dB and number of snapshots 500. The results from Figure 9 show that both the BCS and NP-1 algorithms fail to work when the angular separation is less than 4°, whereas the proposed algorithm can still provide accurate DOA estimation as long as the angular separation is no less than 3°. Moreover, the proposed algorithm has smaller RMSE than the BCS and NP-1 algorithms at the same angular separation. Thus, we can say that the proposed algorithm has a higher angular resolution than the BCS and NP-1 algorithms.Figure 9
RMSE versus angular separation with the fixed SNR 10 dB and number of snapshots 500 for coherent sources.


## 6. Conclusion

In this paper, we develop a DOA estimation algorithm for a mixture of uncorrelated and coherent sources using SBI. In the Bayesian framework, a Gauss-Exp-Chi2 prior is introduced to promote sparse solutions, and the corresponding three-layer hierarchical Bayesian model is established. Then, the model parameters are estimated via variational Bayesian approximations. Finally, we form the source power spectra by using the estimated model parameters, from which the number and the locations of the highest peaks are extracted to achieve source number and DOA estimates. Simulation results show that the proposed algorithm can effectively estimate the DOAs of mixed sources and has better estimation performance than the state-of-the-art BCS algorithm and NP-1 algorithm in terms of estimation accuracy, success rate, and angular resolution.



---

*Source: 3505918-2018-07-10.xml*