Wideband DOA Estimation Utilizing Hierarchical Prior Based on Variational Bayesian Inference

1. Introduction

DOA estimation is an important research focus that has developed rapidly during the past three decades in array signal processing, and its applications are spread across military and civilian fields such as radar, sonar, wireless detection, mobile communication, and biomedicine [1]. It is a remarkable fact that most current researches are based on narrowband signal, while the study for wideband-signal DOA estimation is of equal importance since the common assumption that signals incident from different angles occupy the same frequency band does not always hold in practice [2,3]. In addition, compared to narrowband signals, band signals are characterized by high resolution, which means higher estimation accuracy and better robustness to correlated sources.

Traditional DOA estimation methods for wideband signals can be roughly divided into two categories, one is incoherent signal subspace methods (ISSM) [4,5] and the other is coherent signal subspace methods (CSSM) [6,7,8]. The common point between them is to convert no-spectrum-aliasing wideband signals into narrowband signals for processing, while the most obvious difference is that CSSM can handle coherent sources with lower operational complexity. It is worth noting that there are two methods based on orthogonal testing between ISSM and CSSM, though some literature classifies them as special ISSM category. One is the test of orthogonality of projected subspaces (TOPS), the other is the test of orthogonality of frequency subspaces (TOFS). TOPS and TOFS are with different principles such that they place different emphasis on DOA estimation problems. Their common features are that they can solve the problem of seeking the optimal reference frequency point, and find a balance between the signal to noise ratio (SNR) and the computational complexity. However, all of them need to face the following basic limitations: (ⅰ) Most of them need pre-estimation procedure. (ⅱ) The relatively higher snapshot accumulation is necessary for accurate angular estimation. (ⅲ) A priori of the number of sources is required. (ⅳ) Most of their performance is constrained by the SNR.

For aforementioned limitations, many algorithms based on compressed sensing and sparse representation have been proposed and gradually developed [9,10]. These methods possess better performance than traditional ones such as locating coherent sources well, depending less on high SNR and snapshots and operating with higher efficiency. According to the current research, the methods based on sparse representation for wideband-signal DOA estimation can be broadly classified into three categories: the first one is based on basis pursuit, such as matching pursuit (MP) [11] and orthogonal matching pursuit (OMP) [12]; the second one is based on convex optimization, such as JLZA [13], $l_1$-SVD [14] and W-SpSF [15]; the final one is based on Bayesian compressed sensing [16,17]. The methods’ comprehensive performance of sparse Bayesian learning (SBL) is reported to be better than that based on basis pursuit or convex optimization, though this conclusion was drawn by narrowband signals [18]. There exist several methods based on sparse Bayesian learning, but they all adopted single measure vector model (SMV) [17] or generalized SMV model [16]. In other words, these methods aim at single-snapshot models or vectorized multiple-snapshot models, which would result in a simple decoupling between snapshots. In fact, a reasonable use of temporal correlation not only represents the time-domain coherence of sources, but also improves the algorithm performance. In conclusion, it has been shown that the MMV model has better properties than the SMV model (e.g., stronger joint sparsity and better sparse recovery) [19,20].

In view of this, a novel method, based on the block sparse Bayesian model for wideband-signal DOA estimation, is developed in this paper. The method can be considered as one extension of the SBL based on the SMV model to the MMV model, or as the other extension of the SBL based on narrowband DOA estimation to wideband-signal one. On the basis of SBL, the method takes full advantage of the MMV model by transforming the MMV model into block-sparse model. Not only that, the tactfully design for proposed algorithm allows the joint sparsity between snapshot and spectrum to be exploited simultaneously. In addition, compared to the traditional Gaussian prior, the employed hierarchical structure of the prior can greatly improve the induced sparsity of sparse Bayesian learning.

The rest of this work is organized as follows. In Section 2, the data model is introduced and transformed. In Section 3, the likelihood and prior are presented for variational Bayesian inference and computational complexity of various approaches are compared. The performance of our algorithms is evaluated in Section 4 and the conclusions are drawn in Section 5.

2. Data model

Assume $K$ uncorrelated far-field wideband sources from different directions $\theta =\{{\theta }_1,\mathrm{...},{\theta }_K\}$ impinge on a linear array with $N$ sensors. Without loss of generality, corresponding received data is modeled as

$$x(t)=As(t)+n(t)$$

(1)

where $A=[a({\theta }_1),\mathrm{...},a({\theta }_K)]\in {\coprod }^{N\times K}$ is array manifold matrix with its $k-th$ entry is $a({\theta }_k)={[\exp (-j{v}_1),\mathrm{...}\exp (-j{v}_N)]}^T$, $v_n=2\pi f{d}_n\sin ({\theta }_k)/c$, $n=1,\mathrm{...},N$. $d_n$ denotes the location of the $n-th$ sensor relative to reference one, $c$ is propagation velocity. $s(t)=[s_1(t),\mathrm{...},s_K(t)]$ and $n(t)=[n_1(t),\mathrm{...},n_N(t)]$ are signal waveform vector and additive white Gaussian noise (AWGN) vector, respectively. Then the time-domain data could be converted into frequency-domain data by filter bank (FB) or discrete Fourier transform (DFT).

$$x(f_j)=A(f_j)s(f_j)+n(f_j)$$

(2)

where $f_j$ emphasize that wideband signal is divided into $J$ sub-bands (frequency bins) and their individual center frequency is $f_j$, $j=1,\mathrm{...},J$. In order to ensure the data’s independence in the same sub-band between different snapshots, the total observed time satisfies $T\gg 1/(BL)$, where $B$ is the bandwidth and $L$ is the number of snapshots. On the basis of this, the frequency-domain model of wideband array signal is entirely the same to the time-domain model of narrowband array signal for model structure. In other words, frequency-domain data can be directly treated as time-domain data. After matched filter and snapshots accumulation, (2) could be casted to

$$X_j=A_j{S}_j+N_j$$

(3)

where $A_j=A(f_j)$, $S_j\in {\coprod }^{K\times L}$ denotes the complex amplitude matrix for $L$ snapshots. $N_j\in {\coprod }^{N\times L}$ is the noise matrix. In order to let the model be available for sparse-recovery methods, (3) needs to be converted to following sparse form by sparse representation.

$$X_j={\overline{A}}_j{P}_j+N_j$$

(4)

where ${\overline{A}}_j=[a_j({\overline{\theta }}_1),\mathrm{...},a_j({\overline{\theta }}_M)]\in {\coprod }^{N\times M}$ is the extended manifold matrix, ${\overline{\theta }}_m\in \left\{{\overline{\theta }}_{m=1}^M\right\}$, which is produced by discrete sampling. $P\in {\coprod }^{M\times L}$ is solution matrix and its every row represents a potential source. Therefore, the problem of DOA estimation is casted sparse recovery problem of (4). In terms of [21], (4) can be rewritten as follows.

$$y_j={\Phi }_j{p}_j+n_j$$

(5)

where $y_j=vec(X_j{}^T)\in {\coprod }^{NL\times 1}$, ${\Phi }_j={\overline{A}}_j\otimes I_L\in {\coprod }^{NL\times ML}$,$p_j=vec(P_j{}^T)\in {\coprod }^{ML\times 1}$, $n_j=vec(N_j{}^T)\in {\coprod }^{NL\times 1}$. (5) is block-sparse model since $p_j$ is segmentally sparse (i.e., many segments only contain zeros). As for $p_j$, we just care about the locations of nonzero elements rather than the concrete values, and different $p_j$ indicate the same locations of sources. Therefore, different $p_j$ can be unified as $p$. Consider all the frequency bins, (5) could be extended as

$$y=\Psi p+\overline{n}$$

(6)

where $\overline{y}={[y_1^T,\mathrm{...},y_J^T]}^T\in {\coprod }^{NLJ\times 1}$, $\Psi ={[{\Phi }_1^T,\mathrm{...},{\Phi }_J^T]}^T\in {\coprod }^{NLJ\times ML}$,$\overline{n}={[n_1^T,\mathrm{...},n_J^T]}^T\in {\coprod }^{NLJ\times 1}$.

3. Proposed approach

3.1. Bayesian Model

It is necessary and vital to construct appropriate Bayesian model for sparse Bayesian learning. Priors with hierarchy structure are imposed on hidden variables so as to favouring sparsity more, and the details are as below. As for observed variable $y$, the likelihood is

$$p(y|p;\delta )~CN(\Psi p,{\delta }^{-1}{I}_{NLJ})$$

(7)

Impose Gaussian distribution prior on hidden variable $p$ such that

$$p(p;\gamma )~CN(0,\Sigma )$$

(8)

where $\gamma ={[{\gamma }_{{}_1}^{-1},\mathrm{...},{\gamma }_M^{-1}]}^T$, $\Sigma =diag(\gamma )\otimes I_L\in {\measuredangle }^{ML\times 1}$. Since inverse Gamma distribution is conjugate to Gaussian distribution, the Gamma distribution prior is adopted for each element of $\gamma$.

$$p(\gamma ;a,b)={\displaystyle \prod _{m=1}^M\frac{b^a}{\Gamma (a)}}{\gamma }_m^{a-1}{e}^{-b{\gamma }_m}$$

(9)

where $\Gamma (a)={\displaystyle {\int }_0^{\infty }{x}^{a-1}\exp (-x)dx}$. $a$ is the shape parameter and $b$ is the scale parameter. Similarly, assume $\delta$ obeys Gamma distribution, which yields

$$p(\delta ;c,d)=\frac{c^d}{\Gamma (c)}{\delta }^{c-1}{e}^{-d\delta }$$

(10)

where $c$ and $d$ are corresponding shape and scale parameters, respectively. The directed acyclic graph for representing Bayesian model is as shown in Figure 1.

3.2. Variational Bayesian inference

In the purpose of deriving the iterative algorithm based on SBL, Bayesian inference is adopted. Unluckily, the closed-form solution for posterior cannot be directly obtained due to limitations of mathematics. But the approximate solution would be gained by variational Bayesian inference [21]. The posterior is factorized as below.

$$p(p,\gamma ,\delta |y;a,b,c,d)\approx q(p,\gamma ,\delta )=q(p)q(\gamma )q(\delta )$$

(11)

where $q(p)$, $q(\gamma )$ and $q(\delta )$ are separable marginal distributions of $p$, $\gamma$ and $\delta$. Each of their logarithmic forms can be solved by the others. As for $Inq(p)$, it satisfies

$$Inq(p)={〈Inp(y|p,\delta )p(p;\gamma )〉}_{q(\gamma )q(\delta )}+const$$

(12)

In terms of (7), (8) and (12), $q(p)$ can be solved to obey the Gaussian distribution, whose mean and variance are as below.

$${\mu }_p=〈\delta 〉{\Sigma }_p{\Psi }^{H}y$$

(13)

$${\Sigma }_p={(〈\delta 〉{\Psi }^H\Psi +{〈\Sigma 〉}^{-1})}^{-1}$$

(14)

As usual, $NJ<M$ holds when dense sampling is adopted for high precision estimation, so (13) and (14) could be equivalently transformed into following formulas for less computational complexity.

$${\mu }_p=〈\Sigma 〉{\Psi }^H{({〈\delta 〉}^{-1}{I}_{NLJ}+\Psi 〈\Sigma 〉{\Psi }^H)}^{-1}y$$

(15)

$${\Sigma }_p=〈\Sigma 〉-〈\Sigma 〉{\Psi }^H{({〈\delta 〉}^{-1}{I}_{NLJ}+\Psi 〈\Sigma 〉{\Psi }^H)}^{-1}\Psi 〈\Sigma 〉$$

(16)

Likewise, $q(\gamma )$ satisfies

$$Inq(\gamma )={〈Inp(p;\gamma )p(\gamma ;a,b)〉}_{q(p)q(\delta )}+const$$

(17)

Utilizing (8), (9) and (17), $q(\gamma )$ is identified as a Gamma distribution, whose shape parameter,$m-th$scale parameter and mean are as below.

$$\overline{a}=a+\frac{1}{2}$$

(18)

$${\overline{b}}_m=b+\frac{1}{2}〈p_m^H{p}_m〉=b+\frac{1}{2}{\mu }_{p_m}^H(I_L+diag(diag({\Sigma }_{p_m}))){\mu }_{p_m}$$

(19)

$$〈{\gamma }_m〉=\frac{\overline{a}}{{\overline{b}}_m}$$

(20)

where $p_m,{\mu }_{p_m}\in {\coprod }^{L\times 1}$ is respectively the $m-th$ entry of $p={[p_1^T,\mathrm{...},p_M^T]}^T$, ${\mu }_p={[{\mu }_{p_1}^T,\mathrm{...},{\mu }_{p_M}^T]}^T$, $m=1,2,\mathrm{...},M$. ${\Sigma }_{p_m}={\Sigma }_p([(m-1)L+1:mL],[(m-1)L+1:mL])\in {\coprod }^{L\times L}$is the sub-matrix of ${\Sigma }_p$ according to MATLAB notations.

Similarly, $q(\delta )$ satisfies

$$Inq(\delta )={〈Inp(y|p,\delta )p(\delta ;c,d)〉}_{q(p)q(\gamma )}+const$$

(21)

With (7), (10) and (21), $q(\delta )$ can be also solved as a Gamma distribution and its shape parameter $\overline{c}$, scale parameter $\overline{d}$ and mean are as below.

$$\overline{c}=c+\frac{NLJ}{2}$$

(22)

$$\overline{d}=d+\frac{1}{2}〈{(y-\Psi p)}^H(y-\Psi p)〉=d+\frac{1}{2}{y}^{H}y-real({\mu }_p^H{\Psi }^{H}y)+\frac{1}{2}〈p^H{\Psi }^H\Psi p〉$$

(23)

$$〈\delta 〉=\frac{\overline{c}}{\overline{d}}$$

(24)

For solving the final term of (23), ${\Psi }^H\Psi$ needs to be divided into two parts, i.e., ${({\Psi }^H\Psi )}_{\Lambda }=diag(diag({\Psi }^H\Psi ))$,${({\Psi }^H\Psi )}_{\overline{\Lambda }}={\Psi }^H\Psi -diag(diag({\Psi }^H\Psi ))$. Therefore, (23) can be rewritten as

$$\overline{d}=d+\frac{1}{2}{y}^{H}y-real({\mu }_p^H{\Psi }^{H}y)+\frac{1}{2}{‖{({\Psi }^H\Psi )}_{\overline{\Lambda }}\odot ({\mu }_p^H{\mu }_p)‖}_{1，1}+\frac{1}{2}{‖{({\Psi }^H\Psi )}_{\Lambda }\odot ({\mu }_p^H{\mu }_p+{\Sigma }_p)‖}_{1，1}$$

(25)

So far, the derivation of our proposed iterative algorithm is completed. Utilizing (15), (16), (20) and (24), our iterative algorithm is completed. The specific steps are as follows.

(1) Initialization. First iterative number $k=0$, $a=b=c=d={10}^{-6}$ (ensure uninformative distribution), $p^{(0)}={({\Psi }^H\Psi )}^{-1}{\Psi }^{H}y$. Preset error tolerance$\epsilon$and maximal iterative number $k_{\max }$.

(2) Repetition. Compute hyperparameters with (12), (13), (17) and (22). Iterate alternately until convergence. Regard ${\mu }_p$ as $p^{(k)}{}^{}$ and record $k=k+1$.

(3) Check. If ${‖(p^{(k+1)}-p^{(k)})/p^{(k)}{}^{}‖}_2\le \epsilon$ or $k=k_{\max }$ is satisfied, terminate the algorithm. Otherwise, return to (ii).

(4) Output. Obtain final $p^{(k)}$ and calculate corresponding DOA.

3.3. Computational Complexity

The computational complexity of our algorithm is $O(M{N}^2{L}^3{J}^2)$. When generally dense sampling is adopted, the computational complexity is less than $O(N^2{M}^2)$. The computational complexity of W-SpSF [15], W-SBL [16], $l_1$-SVD [14], JLZA [13] is respectively $O(J^3{M}^3)$, $O(J{M}^3)$,$O(K^3{M}^3)$, $O(M^3+N{M}^2+LNM)$. Generally, our proposed algorithm possesses minimum computational complexity because $M\gg N,L,J,K$ holds in many cases. Of course, the complexity of proposed algorithm is little larger than ISSM or CSSM with many limitations.

4. Numerical Simulation

In this section, the superior performance of our proposed method would be verified by following simulation results. Before that, it is necessary to introduce the Root Mean Square Error (RMSE). RMSE is defined as

$$RMSE=\sqrt{\frac{1}{M_{c}K}{{\displaystyle \sum _{m_c=1}^{M_c}{\displaystyle \sum _{k=1}^K({\widehat{\theta }}_{m_c,k}-{\theta }_k)}}}^2}$$

(26)

where $M_c$ is the Monte Carlo number and ${\widehat{\theta }}_{m_c,k}$ is the angle estimation value for the $k-th$ source in the $m_c-th$ trial, $m_c=1,\mathrm{...},M_c$, $k=1,\mathrm{...},K$.

In order to analyze the performance of our proposed method intuitively, ISSM [4], W-SpSF [15], W-SBL [16], $l_1$-SVD [14] and JLZA [13] are introduced as comparison. The basic experimental conditions are as follows. $K=3$ uncorrelated sources, number of sensors $N=8$, central frequency $f_0=400Hz$, bandwidth $BW=200Hz$, $J=5$, $M_c=300$, grid interval $1^{\circ }$, number of grids $M=180$.Other conditions are indicated in specific experiments.

Experiment 1 tests angular resolution of various algorithms. Other conditions: DOA set $\{-5^{\circ },{35}^{\circ },{60}^{\circ }\}$ (ensure on grid), $SNR=20dB$, number of snapshots $T=20$. As shown in Figure 2, the estimated values of proposed method and ISSM are most precise. Note that proposed method possesses smoother peak values, which expresses small fluctuation with respect to frequency and proves better joint sparsity performance of proposed

Method than ISSM. The results can be explained by the fact that SBL is able to convergence at sparse solutions and possesses robustness to ordinary condition change itself, and proposed method has the ability to enhance sparsity more and stabilize the estimated values of DOA between different frequencies. In a word, proposed method remains the advantages of sparse Bayesian learning and shows better estimation performance than ordinary methods based on SBL.

Experiment 2 examines estimation precision of various algorithms. Other conditions: DOA set $\{-{5.5}^{\circ },{35.1}^{\circ },{60.2}^{\circ }\}$ (ensure off grid). In Figure 3 and Figure 4, four conclusions can be drawn (1) the RMSE performance of proposed method is most excellent; (2) proposed method depend less on number of snapshots than others; (3) the estimation performance of proposed method improves most with SNR increasing; (4) proposed method achieves perfect estimation performance.

5. Conclusions

In this paper, we further expand the use and performance of sparse Bayesian learning with respect to wideband signal by transforming MMV model into block-sparse model. The main contributions are as below: the hierarchical Bayesian prior framework is constructed to favour sparsity; the corresponding iterative process is derived and completed; the whole method is proposed and applied to solve underdetermined DOA estimation problems. In some way, our proposed mehtod fully succeeds the advantages of sparse Bayesian learning and gains more excellent estimation performance and sparse recovery ability. Simulation results of two experiments are presented to prove the superior of our method, such as favouring sparsity, joint sparsity, high estimation accuracy, stable estimated values and strong adaptability.