Bayesian learning scheme for sparse DOA estimation based on maximum-a-posteriori of hyperparameters

Received Jul 22, 2020 Revised Oct 23, 2020 Accepted Dec 18, 2020 In this paper, the problem of direction of arrival estimation is addressed by employing Bayesian learning technique in sparse domain. This paper deals with the inference of sparse Bayesian learning (SBL) for both single measurement vector (SMV) and multiple measurement vector (MMV) and its applicability to estimate the arriving signal’s direction at the receiving antenna array; particularly considered to be a uniform linear array. We also derive the hyperparameter updating equations by maximizing the posterior of hyperparameters and exhibit the results for nonzero hyperprior scalars. The results presented in this paper, shows that the resolution and speed of the proposed algorithm is comparatively improved with almost zero failure rate and minimum mean square error of signal’s direction estimate.


INTRODUCTION
Direction of arrival (DOA) estimation is a well-known problem in the field of array signal processing where, the angle of arrival (direction) of the signal at the receiver needs to be estimated from the knowledge of the received signal itself. This problem of direction of arrival estimation has attracted many modern researchers because of its large range of applications in certain fields such as RADAR, SONAR, seismology, wireless mobile communication and others. To solve this problem of signal's direction estimation, an array of antennas with linear or non-linear structure having uniform or non-uniform antenna spacing can be used at the receiver. The signals which are generated from the far-field sources arrive at a particular direction and impinge on the antenna array. The received signals from the antenna array of 'M' sensors form the under-sampled observed signal samples. These observed signal samples contain the direction information and hence they are processed to estimate the signal source direction [1]. From past two decades, many algorithms were derived to solve the problem of DOA estimation. These algorithms can be broadly classified into: i) Conventional methods, ii) Subspace methods, iii) Sparse methods.
The standard MUSIC algorithm proposed in [2,3], decomposes signal and noise subspaces along with multiple signal classification methodology to estimate the number of signal sources as well as the spatial spectrum of the received signal. The drawback of this algorithm is that it fails for coherent signal sources. In [4], an improved and modified MUSIC algorithm is proposed by employing matrix decomposition to address the case of coherent signal sources but the performance of this algorithm deteriorates for low SNR region. In [5,6], the performance of all these subspace based standard DOA estimation algorithms are analyzed and found that these techniques offer good speed and less complexity but suffer from low resolution, high sensitivity towards correlated signal sources and high MSE.
In recent years, after the emerge of sparse signal representation several algorithms were derived as solutions to the DOA estimation problem by considering and representing the problem as sparse signal recovery problem by utilizing the sparse nature of to be estimated signal [7]. The sparse based algorithm proposed in [8,9] is based on the orthogonal matching pursuit (OMP) technique, which suffers from low performance and requires the knowledge of source number. In [10], a convex relaxation with l1 norm penalty is applied to mitigate the perturbation effects. In [11,12], an l1-svd based algorithm along with re-weighted l1 minimization is proposed to minimize the complexity of DOA estimation algorithm but suffers with low performance for coherent and closely spaced signal sources. The compressive sensing based DOA estimation algorithms proposed in [13][14][15][16] are based on the simple least squares minimization method which offers good MSE, resolution but suffers from high complexity, the performance of BP and OMP depends on the arraysteering matrix [9], and it degrades for highly correlated array-steering matrix in the case of DOA problem. The scaling/shrinkage operations in convex relaxation may reduce variance for increase in sparsity or vice versa.
In the most recent years, Bayesian methods like maximum a posteriori (MAP) [17,18], maximum likelihood (ML estimation) [18], iterative reweighted l1 and l2 algorithms were applied to solve the DOA estimation problem. These Bayesian algorithms suffer from high MSE, even though true priors are used. In [17], MAP only guarantees maximization of product of likelihood and the prior of the unknown sparse signal. ML estimate in [18], also maximizes only the likelihood function by assuming prior of unknown to be equally likely to occur.
Sparse Bayesian learning (SBL) with relevance vector machine proposed by Tipping in [19] and rerepresented by Wipf in [20] for linear regression/sparse signal recovery problem led a broader way with higher performance results in the research of sparse signal recovery. In this paper, we present the detail inference of Sparse Bayesian learning and its applicability to DOA problem using on-grid approach. We also derive the updating equations for hyperparameters by maximizing the posterior of hyperparameters for nonzero hyperprior scalars.
Further, the paper is organized as: Section 2 describes the signal model used for DOA estimation for a uniform linear array. Section 3 describes the basics and inference of the Sparse Bayesian Learning technique. Section 4 describes about updating of the hyperparameters of SBL estimate by proposing a method of maximum-a-posterior of the hyperparameters. Section 5 summarizes the proposed algorithm. In section 6, the results and performance analysis of the proposed algorithm are presented. Finally, section 7 concludes the paper.

SIGNAL MODEL FOR SPARSE DOA ESTIMATION
Consider 'D' number of arriving signal sources s(n)=[s1(n), s2(n)…. sD(n)] T impinging on the uniform linear array of 'M' sensors with a uniform spacing of d ≤ λ/2, where λ is the wavelength of the arriving signals. Let y(n)=[y1(n), y2(n)…. yM(n)] T be Mx1 observed signal samples received by 'M' antenna array sensors. For simplicity, assuming a single snapshot (single measurement vector) i.e, n=1, the problem of direction of arrival estimation can be modeled as in (1).
where A is MxN array steering matrix given by (2) and a(θi) represents the atom for a particular direction angle θi. For searching the entire angle space for DOA a particular grid of 'N' values of angles are considered. Each atom is a vector of Mx1 antenna array steering vector given in (3).
where, β=2π/λ, x(n)=[x1(n), x2(n)…. xN(n)] T is Nx1 signal vector that needs to be estimated to find the source signal directions in presence of antenna array noise vector w(n)=[w1(n), w2(n)…. wM(n)] T of Mx1 size. The estimated x(n) values are the estimation of signal power s(n) and is related by (4). The model in (1) turns out to be the problem of sparse signal recovery from the under-sampled measurements y [13].

SPARSE BAYESIAN LEARNING INFERENCE
Consider a single snapshot case, DOA estimation problem as in (1). The w(n) are independent noise samples which is assumed to be zero-mean Gaussian random process with noise variance σ 2 . By Bayes's theorem [19], the posterior of unknown x by knowing the observed antenna array received signal y can be expresses as in (5).
where, P(y/x) is the likelihood of observed data for the estimated unknown parameter 'x' which is assumed as P(y/x) = (y/Ax, σ 2 ), where the notation (. ) specifies a guassian distribution over y with mean Ax and variance σ 2 . This assumption is due to another assumption of independence of samples y. Thus the likelihood function of y is given in (6).
The prior of unknown 'x' is also assumed to be as zero mean Gaussian prior distribution over x with variance [19,20]. The Gaussian prior of a single sample of x (i.e, xi) is given in (7).
The overall Gaussian prior for all i = 1 to N is given in (8).
To define prior of unknown x, we require another parameter γ which is variance of unknown x. Thus the γ can be called as a vector of hyperparameters of unknown 'x' [21,22]. Hence, to completely define all the distributions, the hyperparameters γ and noise variance σ 2 needs to be estimated which can be done by defining the hyperpriors of γ and σ 2 as in (9) and (10).
We have chosen gamma distribution because the hyperparameters γ and σ 2 are scale parameters [23,24] where: is the gamma function and a, b, c, d are all hyperprior parameters. After defining the likelihood and priors (5) can be re-written as [19].
In (14) we get one more result for prior of the observed array received signal vector.
with Σ y = (σ 2 I + AΓ −1 A T ) as the prior covariance of observed array received signal vector. Solving (16) and (17) for the known values of hyperparameters γ and σ 2 results in mean and covariance of posterior of unknown x respectively. The posterior mean of unknown x is itself the estimation of the unknown x i.e, ̂= , plotting this ̂ estimate with respect to the on-grid search angle points gives the DOA peaks and hence the arriving signal source's direction can be estimated [25]. In practical situations, the hyperparameters γ and σ 2 will be unknown and there cannot be any closed form expressions obtained for them [26]. Hence, an iterative estimation of hyperparameters γ and σ 2 has to be done.

MAXIMUM A POSTERIOR OF HYPERPARAMETERS
To iteratively estimate the hyperparameters like variance of prior of unknown and variance of noise σ 2 , we maximize the probability function P(γ,σ 2 / ) given by (19). P(γ,σ 2 /y)= P(y/γ,σ 2 )P(γ,σ 2 ) The hyperparameters γ, σ 2 is mutually independent with each other and also the probability of known measured array received signal vector is a constant. Thus, maximizing (19) is equivalent to maximize (20) with respect to γ,σ 2 . P(γ,σ 2 /y)∝P(y/γ,σ 2 )P(γ)P(σ 2 ) As in practice, we assume uniform hyperpriors over a logarithmic scale with the derivatives of the hyperpriors terms goes to zero, we choose to maximize the logarithmic quantity of (20) with respect to log γ and logσ 2 . The logarithm of (20) is given by (21).

The variance of unknown 'x'
The differentiation of (23) with respect to logγ i gives: Setting this (24) to zero and Mackay [27], leads to the update in (25).

The noise variance
The differentiation of (23) with respect to logσ 2 we get: where: γ i −1 Σ xii ) and setting derivative in (26) to zero and re-arranging the terms we get σ 2 update as in (27).

THE PROPOSED ALGORITHM
For a single snapshot case, in DOA estimation, initialize the noise variance σ 2 and the prior variance γ of unknown x to a value (i.e, usually taken as 1). Using (17) will give the covariance estimate and later using (16) will give the first iterative estimate of μ. Before performing the 2 nd iterative estimation of Σ x and μ, let us update the hyperparameters γ and σ 2 using (26) and (27), where the parameters/variables in those equations represent the values of 1 st iteration. Now using these new updated values of γ and σ 2 , estimate the 2 nd iteration values of Σ x and μ. Repeat these steps until a particular stopping criterion is achieved. In this iterative process, some elements of μ vector tend to become very minimum value (i.e, less than a preset threshold), equating these elements to zero, results in sparsity of the solution.
For 'L' number of multiple snapshot/multiple measurement vector (MMV) case also, same procedure can be utilized except that the prior mean of unknown 'x' (i.e, μ) is a matrix, in which each row corresponds to a particular on-grid search point of angle of arrival. Each of these rows of μ for MMV case should be taken as absolute mean square values of all the elements of that particular row. This μ estimate obtained at the final iteration of the proposed algorithm is plotted versus the search grid of angle of arrival. The plot showing peaks corresponding to the particular value of angle of arrival on x-axis, indicates the estimate of direction of arrival. The proposed DOA estimation algorithm based on sparse Bayesian learningmaximum a posterior of hyperparameters (SBL-MAP-H) for MMV case is summarized in Table 1.

RESULTS AND DISCUSSION
In this section, the experimental results of the proposed algorithm is presented for different conditions of various algorithmic parameters. For the simulation of the algorithm, MATLAB R2013a platform has been utilized. The simulation results of the proposed algorithm are compared with standard DOA estimation algorithms like MUSIC [2,3,28], MVDR [29,30] and the recent standard algorithm l1-SVD [11,12]. Considering a uniform linear array (ULA) of M=100 number of array elements with an interelement spacing of λ/2, where λ stands for wavelength of the received signal assumed to be as 1m. Let us assume a single signal source transmitted from a far-field with a direction of 0 0 with respect to the vertical normal axis having an angular frequency of 20π r/s is impinging on the ULA. The proposed algorithm considers a set of on-grid points for searching the direction of the arriving signal with a 0.5 0 step-size. The proposed algorithm is simulated for L=500 number of snapshots in a noisy environment with SNR of 30dB. The hyperparameter updating depends on the hyperprior parameters (a,b,c,d). These parameters highly influence the DOA estimation results as shown in Figure 1. For abcd-parameter values equal to zero, the DOA estimation peak is less steep when compared to the estimation peak obtained for abcd-parameters equal to 0.4. It is also tested with various other values of a,b,c,d and found that for all 0<a,b,c,d<0.5 gives steepest estimation peaks containing maximum peak only at the actual angle of arrival of the received signal and completely flat response for any other grid points. The very high value set for a,b,c,d increases the sparsity in the estimated results and in some cases with weak signal strength, the actual true DOAs also may not contain the estimation peaks. Hence the range of 0<a,b,c,d<0.5 is the optimized option for DOA estimation application. All the next analysis considers abcd-parameters as 0.4. Considering M=10, number of search grid points as N=361, L=100, number of signal sources D=3 with actual true DOAs as -10 0 , 10 0 , 64 0 with corresponding angular frequencies of 20 π, 40 π, 60 π r/s respectively and a noisy environment with SNR 0dB. Figure 2 shows the DOA estimation peaks for the proposed algorithm as well as various standard DOA estimation algorithms. It can be observed that though the value of SNR is very less (i.e, the worst noisy environment), the proposed algorithm shows sharp DOA estimation peaks at the actual true DOAs.
For the same parametric conditions, considering a single snapshot case with L=1 also gives sharp DOA estimation peaks indicating the actual true DOAs with 100% success rate as shown in Figure 3. Figure 4 indicates the estimation case for the number of array elements M=100, which shows that the proposed algorithm performance is almost similar to that for M=10 with respect to mean square error. In the case of very closely spaced source signals with actual true DOA of 10 0 and 11 0 , the proposed algorithm still produces steeper peaks with clear distinguished DOA peaks as compared to other standard algorithms as shown in Figure 5. Figure 6 indicate the case of very closely spaced two coherent signal sources with an angular frequency of 20 π r/s and located at 0 0 and 1 0 . The result in Figure 6 exhibits the high-resolution performance of the proposed algorithm compared to the other algorithms. As the proposed algorithm employs the probability of the measured antenna array signal by knowing the prior of unknowns, good resolution, even for coherent signal sources are obtained. The effect of array sensor noise added up with the received signal for hyperprior parameters a,b,c,d=0 is as shown in Figure 7. Considering L=50, M=100 and a single source signal with actual true DOA of 0 0 , the DOA peak becomes more steeper along with decrease in mean square error for the improvement in SNR value.    Figure 8 shows the effect of number of snapshots 'L' on the DOA estimation peaks for SNR=10 dB, M=100 and actual true DOA of 0 0 . As the number of snapshots increases, the DOA peaks become more steeper with better performance. For the case of increase in number of array elements in the ULA, the DOA estimation performance increases along with the increase in estimation success rate as shown in Figure 9. For a single source arriving at actual true DOA of 0 0 with L=50 and M=100, the performance analysis of various standard algorithms compared with the proposed algorithm with respect to mean square error v/s signal to noise ratio. As seen in Figure 10, the proposed algorithm shows less MSE for all the range of SNR, when compared with other standard DOA estimation algorithms. It is true that the proposed algorithm also exhibits very least failure rate with respect to the SNR as shown in Figure 11.
The execution time consumed by the proposed algorithm is quiet more when compared to MUSIC and other subspace-based algorithms. However, the proposed algorithm's execution time is comparatively less with respect to l1-SVD algorithm and almost becomes equal as the number of snapshots increases as shown in Figure 12. The probability of success rate for measuring the resolution performance of the proposed algorithm is as shown in Figure 13. By maintaining constant antenna array elements in the case of closely angular spaced coherent signal sources, as the SNR of the received antenna array signal increases, the resolution/probability of success rate of estimation of the proposed algorithm also increases and produces almost 100% resolution from SNR range of 0 dB onwards.

CONCLUSION
In this paper, a sparse Bayesian learning approach based on maximum a posterior of hyperparameters for DOA estimation is designed and tested for various DOA estimation conditions and parameters. The proposed algorithm exhibits good mean square error and success rate for all the parametric conditions like low SNR range, less array size and a smaller number of snapshots. It also results in good resolution for very closely spaced signal sources. From the result section, it is observed that the proposed algorithm achieves better MSE v/s SNR performance when compared with other standard DOA estimation algorithms in low to medium SNR range. The only exception of the proposed algorithm is that more computation time with increased complexity for larger number of snapshots. As the proposed algorithm shows good performance for single or a very few snapshot cases with manageable computation time, it can be comparatively of more important solution for the problem of DOA estimation in array signal processing field.