Bayesian Compressive Sensing Approaches for Direction of Arrival Estimation with Mutual Coupling Effects

TL;DR

This paper introduces two Bayesian compressive sensing methods for accurately estimating the direction of arrival of signals with mutual coupling effects, especially effective in challenging endfire regions.

Contribution

It develops a Bayesian compressive sensing Kalman filter and a Gibbs sampling approach with explicit sparsity enforcement, advancing DOA estimation accuracy under mutual coupling.

Findings

Achieved up to 86.36% reduction in RMSE with the Gibbs sampling method.

Proposed methods outperform traditional relevance vector machine approaches.

Enhanced DOA estimation accuracy in endfire regions of the array.

Abstract

The problem of estimating the dynamic direction of arrival of far field signals impinging on a uniform linear array, with mutual coupling effects, is addressed. This work proposes two novel approaches able to provide accurate solutions, including at the endfire regions of the array. Firstly, a Bayesian compressive sensing Kalman filter is developed, which accounts for the predicted estimated signals rather than using the traditional sparse prior. The posterior probability density function of the received source signals and the expression for the related marginal likelihood function are derived theoretically. Next, a Gibbs sampling based approach with indicator variables in the sparsity prior is developed. This allows sparsity to be explicitly enforced in different ways, including when an angle is too far from the previous estimate. The proposed approaches are validated and evaluated over different test scenarios and compared to the traditional relevance vector machine based method. An improved accuracy in terms of average root mean square error values is achieved (up to 73.39% for the modified relevance vector machine based approach and 86.36% for the Gibbs sampling based approach). The proposed approaches prove to be particularly useful for direction of arrival estimation when the angle of arrival moves into the endfire region of the array.

Tables

Table 1. TABLE I: Performance summary for the endfire region example.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	19.88	0.76
Modified RVM	6.46	0.98
Gibbs	3.82	17.83 excluding burn-in
		107.84 including burn-in

Table 2. TABLE II: Performance summary for the endfire region example with reduced adjacent antenna separations.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	2.21	0.98
Modified RVM	0.86	1.25
Gibbs	3.02	16.81 excluding burn-in
		101.65 including burn-in

Table 3. TABLE III: Performance summary for the non-endfire region example.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	2.91	0.73
Modified RVM	1.85	0.78
Gibbs	1.46	13.08 excluding burn-in
		79.57 including burn-in

Table 4. TABLE IV: Performance summary for the random initial DOA example.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	12.10	0.74
Modified RVM	3.22	0.86
Gibbs	2.89	13.88 excluding burn-in
		84.19 including burn-in

Table 5. TABLE V: Performance summary for the increasing DOA followed by decreasing DOA example.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	4.49	0.69
Modified RVM with	4.45	0.77
constant $+1^{\circ }$ DOA change
Modified RVM with	4.08	0.79
constant $-1^{\circ }$ DOA change
Gibbs	1.41	16.12 excluding burn-in
		97.59 including burn-in

Table 6. TABLE VI: Performance summary for the increasing DOA with an assumed decrease in DOA of 3 ∘ superscript 3 3^{\circ} example.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	10.70	0.66
Modified RVM	8.07	0.77
Gibbs	1.46	13.99 excluding burn-in
		84.72 including burn-in

Table 7. TABLE VII: Performance summary for the random changes in DOA with an assumed increase in DOA of 3 ∘ superscript 3 3^{\circ} example.

	Average $RMSE$	Average Computation
Method	(degrees)	Time (seconds)
RVM	6.89	0.64
Modified RVM	6.37	0.86
Gibbs	1.46	22.26 excluding burn-in
		134.56including burn-in

Equations

$$\text{a}(\Omega ,\theta )=[1,e^{-j{\mu }_2\Omega \cos \theta },\dots ,e^{-j{\mu }_M\Omega \cos \theta }{]}^T,$$

$${\text{y}}_k={\text{A}}_{st}{\text{x}}_k+{\text{n}}_k,$$

$$\text{A}={\text{M}}_{MC}{\text{A}}_{st}.$$

$$\textbf{M}_{MC}=\left(\begin{array}[]{cccccc}1&m_{2}&\ldots&m_{D-1}&\ldots&m_{M}\\ m_{2}&1&m_{2}&\ldots&\ddots&\vdots\\ \vdots&m_{2}&1&m_{2}&\ldots&m_{D-1}\\ m_{D-1}&\ldots&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ldots&m_{2}&1&m_{2}\\ m_{M}&\ldots&m_{D-1}&\ldots&m_{2}&1\\ \end{array}\right).\vspace*{.5cm}$$

$${\text{y}}_k$$

$${\tilde{\text{y}}}_k$$

$$\displaystyle\left[\begin{array}[]{c}\mathcal{R}(\textbf{y}_{k})\\ \mathcal{I}(\textbf{y}_{k})\\ \end{array}\right]$$

$${\tilde{\text{x}}}_{k,opt}=\max \mathcal{P}({\tilde{\text{x}}}_k,{\sigma }_k^2,{\text{p}}_k|{\tilde{\text{y}}}_k,{\tilde{\text{x}}}_p),$$

$$\mathcal{P}(\tilde{\textbf{y}}_{k}|\tilde{\textbf{x}}_{k},\sigma_{k}^{2})=(2\pi\sigma_{k}^{2})^{-M}\exp\Big{\{}-\frac{1}{2\sigma^{2}}||\tilde{\textbf{y}}_{k}-\tilde{\textbf{A}}\tilde{\textbf{x}}_{k}||_{2}^{2}\Big{\}}.$$

$$\mathcal{P}({\tilde{\text{x}}}_k|{\text{p}}_k,{\tilde{\text{x}}}_p)$$

$$\mathcal{P}({\text{p}}_k)=\prod _{n=1}^{2N}G(p_{k,n}|{\beta }_1,{\beta }_2).$$

$$\mathcal{P}({\sigma }_k^2)=G({\sigma }_k^{-2}|{\beta }_3,{\beta }_4),$$

$$\mathcal{P}({\tilde{\text{x}}}_k,{\sigma }_k^2,{\text{p}}_k|{\tilde{\text{y}}}_k,{\tilde{\text{x}}}_p)=\mathcal{P}({\tilde{\text{x}}}_k|{\tilde{\text{y}}}_k,{\sigma }_k^2,{\text{p}}_k,{\tilde{\text{x}}}_p)\mathcal{P}({\sigma }_k^2,{\text{p}}_k|{\tilde{\text{y}}}_k,{\tilde{\text{x}}}_p)$$

$$\mathcal{P}({\tilde{\text{x}}}_k|{\tilde{\text{y}}}_k,{\sigma }_k^2,{\text{p}}_k,{\tilde{\text{x}}}_p)=\frac{\mathcal{P}({\tilde{\text{y}}}_k|{\tilde{\text{x}}}_k,{\sigma }_k^2)\mathcal{P}({\tilde{\text{x}}}_k|{\text{p}}_k,{\tilde{\text{x}}}_p)}{\mathcal{P}({\tilde{\text{y}}}_k|{\sigma }_k^2,{\text{p}}_k,{\tilde{\text{x}}}_p)}$$

$$\displaystyle=(2\pi)^{-N}|\boldsymbol{\Sigma}|^{-1/2}\exp\Bigg{\{}-\frac{1}{2}(\tilde{\textbf{x}}_{k}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\tilde{\textbf{x}}_{k}-\boldsymbol{\mu})\Bigg{\}},$$

$$\mathbf{\Sigma }=({\sigma }_k^{-2}{\tilde{\text{A}}}^T\tilde{\text{A}}+{\text{P}}_k{)}^{-1}$$

$$\mathbf{\mu }=\mathbf{\Sigma }({\sigma }_k^{-2}{\tilde{\text{A}}}^T{\tilde{\text{y}}}_k+{\text{P}}_k{\tilde{\text{x}}}_p),$$

$$\mathcal{P}({\sigma }_k^2,{\text{p}}_k|{\tilde{\text{y}}}_k,{\tilde{\text{x}}}_p)\approx \mathcal{P}({\tilde{\text{y}}}_k|{\sigma }^2,\text{p},{\tilde{\text{x}}}_p)\mathcal{P}({\sigma }_k^2)\mathcal{P}({\text{p}}_k)\mathcal{P}({\tilde{\text{x}}}_p),$$

$$\mathcal{L}({\sigma }_k^2,{\text{p}}_k)$$

$$p_{k,n}^{new}=\frac{{\gamma }_n}{{\mu }_n^2+{\tilde{x}}_{p,n}^2-{\tilde{x}}_{p,n}{\mu }_n},$$

$${\sigma }_{k,new}^2=\frac{||{\tilde{\text{y}}}_k-\tilde{\text{A}}\mathbf{\mu }|{|}_2^2}{2M-\sum _n{\gamma }_n}.$$

$$\tilde{\textbf{x}}_{k,opt}=\Big{(}\frac{\tilde{\textbf{A}}^{T}\tilde{\textbf{A}}}{\sigma_{k,opt}^{2}}+\textbf{P}_{k,opt}\Big{)}^{-1}\Big{(}\frac{\tilde{\textbf{A}}^{T}\tilde{\textbf{y}}_{k}}{\sigma_{k,opt}^{2}}+\textbf{P}_{k,opt}\tilde{\textbf{x}}_{p}\Big{)}$$

$$x_{k,opt,n}={\tilde{x}}_{k,opt,n}+j{\tilde{x}}_{k,opt,N+n},$$

$${\tilde{\text{x}}}_{k|k-1}$$

$${\tilde{\text{y}}}_{k|k-1}$$

$${\tilde{\text{x}}}_k={\tilde{\text{x}}}_{k|k-1}+{\text{K}}_k{\tilde{\text{y}}}_{e,k}$$

$${\text{K}}_k={\mathbf{\Sigma }}_{k|k-1}{\tilde{\text{A}}}^T$$

$$\mathcal{L}({\sigma }_k^2,{\text{p}}_k)$$

$$\mathcal{P}({\tilde{\text{x}}}_k|{\text{p}}_k,{\tilde{\text{z}}}_k)=\prod _{n=1}^{2N}(1-{\tilde{z}}_{k,n}){\delta }_0+{\tilde{z}}_{k,n}\mathcal{N}({\tilde{x}}_{k,n}|0,p_{k,n}),$$

$$z_{k,n}=\{\begin{array}{ll}{z}_{k,n}^1& \text{ if }|n-i|\le j,\\ z_{k,n}^2& \text{ if }|n-i|>j,\end{array}$$

$$z_{k,n}^1$$

$$z_{k,n}^2$$

$$\mathcal{P}(z_{k,n}^1|{\beta }_5,{\beta }_6)$$

Full text

I Introduction

Direction of arrival (DOA) estimation is the process of determining which direction a signal impinging on an array has arrived from. Commonly used methods of solving this problem are: MUSIC [1, 2], ESPRIT [3, 4, 5, 6] and the maximum likelihood DOA estimator [7, 8, 9]. However, these methods have some disadvantages, in particular they require knowledge of the number of signals present beforehand and evaluation of a covariance matrix of the array output (adding computational complexity).

Compressive Sensing (CS) theory says that when certain conditions are met it is possible to recover signals from fewer measurements than used by traditional methods [10, 11]. Hence, CS can be applied to the problem of DOA estimation [12, 13, 14, 15] by splitting the angular region into $N$ potential DOAs, where only $L<<N$ of the DOAs have an impinging signal (alternatively $N-L$ of the angular directions have a zero-valued impinging signal present). These DOAs are then estimated by finding the minimum number of DOAs with a non-zero valued impinging signal that still give an acceptable estimate of the array output.

The problem can also be converted into a probabilistic form and solved via Bayesian compressive sensing (BCS) [16], implemented with a relevance vector machine (RVM) [17, 18, 19]. Such a method has been used to solve the problem of static DOA estimation [20, 21], where a belief of having a sparse received signal is made and the most likely values found.

The Kalman filter (KF) can be used to track dynamic DOAs, with the angular range narrowed to focus in more closely on the DOA estimate from the previous iteration [22]. However, this prevents directly working with the measured array signals and introduces an additional stage of having to reevaluate the steering vector of the array at each iteration of the KF. Hierarchical KFs have been used to track dynamic sparse signals [23, 24], where the predicted mean of the signals at each iteration is taken as the estimate from the previous iteration and the hyperparameters are estimated using BCS, hence the term Bayesian compressive sensing Kalman Filter (BCSKF).

However, a problem remains when a BCSKF is applied to dynamic DOA estimation with a uniform linear array (ULA). The estimation accuracy can be reduced when the DOA approaches the endfire region of the array, i.e. when the impinging signal arrives parallel to or almost parallel to the array. This can be particularly problematic when there is a lot of noise present.

An additional challenge to address when considering the DOA estimation problem is that of mutual coupling. One way of modeling the mutual coupling effects is to use a mutual coupling matrix [25, 26]. In [25] the mutual coupling matrix is found using two methods: minimum mean-square matching and the mutual impedance method. The method in [26] applies a symmetric Toeplitz matrix, where only antennas within a set separation of each other can cause mutual coupling effects. In this work the method in [26] is used to ensure mutual coupling effects are included in the signal model.

The contributions of this paper are: i) A BCSKF with a modified RVM, where the traditional sparsity prior is replaced with a belief that the estimated signals will instead match predicted signal values, is proposed. The result of this new prior is that a new posterior distribution and marginal likelihood have been derived. Initial results for this method using a signal model without mutual coupling have been reported in [27]. ii) A Gibbs sampling approach is proposed. In this approach zero valued signals can be explicitly enforced when there is too large a change in DOA in order to alleviate the estimation accuracy problem for the endfire region of the array. iii) A comprehensive performance evaluation is provided, with the proposed methods being compared to a BCSKF using the traditional RVM approach. Significant improvements in terms of the average root mean square error ($RMSE$) values are observed (up to 73.39$\%$ for the BCSKF with modified RVM and up to 86.36$\%$ for the Gibbs sampling approach).

The remainder of this paper is structured in the following manner: Section II gives details of the proposed estimation methods, including the array model with mutual coupling effects (II-A), the modified RVM framework for BCS (II-B), the BCSKF (II-C) and the Gibbs sampling implementation (II-D). In Section III an evaluation of the effectiveness of the proposed approaches is presented and conclusions are drawn in Section IV.

II Proposed Estimation Methods

II-A Array Model

A narrowband ULA structure consisting of $M$ omnidirectional antennas, with identical responses is shown in Figure 1. Here, a plane-wave signal mode is assumed, i.e. the signal impinges upon the array from the far field and the angle of arrival is limited to $0^{\circ}\leq\theta\leq 180^{\circ}$. The distance from the first antenna to subsequent antennas is denoted as $d_{m}$ for $m=1,2,\ldots,M$, with $d_{1}=0$, i.e. the distance from the first antenna to itself. Note, these values are multiples of a uniform adjacent antenna separation of $\Delta d$.

The steering vector of the array is given by

[TABLE]

where $\Omega=\omega T_{s}$ is the normalised frequency with $T_{s}$ being the sampling period, $\mu_{m}=\frac{d_{m}}{cT_{s}}$ for $m=1,2,\ldots,M$, $c$ gives the wave propagation speed and $\{\cdot\}^{T}$ denotes the transpose operation.

The array output, $\textbf{y}_{k}$, at time snapshot $k$ is then given by

[TABLE]

where $\textbf{x}_{k}=[x_{k,1},x_{k,2},...,x_{k,N}]^{T}\in\mathbb{C}^{N\times 1}$ gives the received source signals, $\textbf{n}_{k}=[n_{k,1},n_{k,2},...,n_{k,M}]^{T}\in\mathbb{C}^{M\times 1}$ is a noise term, given by a zero mean multivariate Gaussian random variable and $\textbf{A}_{st}=[\textbf{a}(\Omega,\theta_{1}),\textbf{a}(\Omega,\theta_{2}),...,\textbf{a}(\Omega,\theta_{N})]\in\mathbb{C}^{M\times N}$ is the matrix containing the steering vectors for each angle of interest. Note, $N$ is the number of points in the grid of potential DOAs the angular region has been split into. However, only $L<<N$ of these angular directions will have an impinging signal present.

In practice there will also be mutual coupling effects present, which alter the pattern of an individual antenna as compared to if it was being used on its’ own. As a result (2) has to be altered to account for this fact. A mutual coupling matrix is used to achieve this [26], by giving the true steering vector matrix as

[TABLE]

Here $\textbf{M}_{MC}\in\mathbb{C}^{M\times M}$ is the mutual coupling matrix given by

[TABLE]

In (4) the mutual coupling coefficients are given by $m_{i}~{}=~{}\rho_{i}\exp\{j\phi_{i}\}$ for $i=2,...,D-1,D,...,M$, where $\rho_{i}$ and $\phi_{i}$ give the amplitude and phase, respectively. The variable $D$ places a limit on the separation between antennas above which there will be no mutual coupling effects. In other words when $i>D$, then $\rho_{i}=0$. This then gives the following:

[TABLE]

Equation (5) can then be split into real and imaginary components (given by $\mathcal{R}(\cdot)$ and $\mathcal{I}(\cdot)$, respectively) as follows

[TABLE]

The difference between $\textbf{y}_{k}$ and $\tilde{\textbf{y}}_{k}$ is that $\textbf{y}_{k}$ has been split into its real and imaginary components in $\tilde{\textbf{y}}_{k}$. As a result the dimensions of $\tilde{\textbf{y}}_{k}$ are increased. A similar relationship exists between A and $\tilde{\textbf{A}}$, $\textbf{x}_{k}$ and $\tilde{\textbf{x}}_{k}$ and $\textbf{n}_{k}$ and $\tilde{\textbf{n}}_{k}$.

II-B Modified Relevance Vector machine for DOA Estimation

The aim is to now find a solution for $\tilde{\textbf{x}}_{k}$ which gives the closest possible match to a predicted set of signal values. To achieve this one can follow a modified RVM framework [27], by evaluating the following

[TABLE]

where $\sigma_{k}^{2}$ is the variance of the Gaussian noise $\textbf{n}_{k}$, $\textbf{p}_{k}=[p_{k,1},p_{k,2},...,p_{k,2N}]^{T}$ contains the hyperparameters that are to be estimated and $\tilde{\textbf{x}}_{p}=[\mathcal{R}(\textbf{x}_{p})^{T},\mathcal{I}(\textbf{x}_{p})^{T}]^{T}=[\mathcal{R}(x_{p,1}),...,\mathcal{R}(x_{p,N}),\mathcal{I}(x_{p,1}),...,\mathcal{I}(x_{p,N})]^{T}$ holds the predicted values of $\tilde{\textbf{x}}_{k}$.

From (II-A) it is possible to find:

[TABLE]

The traditional RVM would now apply a belief that $\tilde{\textbf{x}}_{k}$ is sparse. However, here this is changed to a belief that $\tilde{\textbf{x}}_{k}$ will match the predicted signals $\tilde{\textbf{x}}_{p}$:

[TABLE]

Note, when $\tilde{\textbf{x}}_{p}=[0,0,...,0]^{T}$ then (17) reverts to the hierarchical prior used in the traditional RVM [16, 17] and $|\textbf{P}_{k}|$ indicates the determinant of $\textbf{P}_{k}$, where $\textbf{P}_{k}=\text{diag}(\textbf{p}_{k})$.

It is also necessary to define the hyperparameters over $\textbf{p}_{k}$ and $\sigma_{k}^{2}$. There are various possibilities for the structuring of the priors on $\textbf{p}_{k}$, which represent mixing parameters in a scale mixture of normals representation of the marginal distribution of $\textbf{x}_{k}$, which will here be in the Student-t family, see e.g. [28]. One possibility would be to treat the complex components of $\textbf{x}_{k}$ as complex Student-t distributed, as detailed in [29, 30]. However, this work treats the real and imaginary components of $\textbf{x}_{k}$ as independent Student-t distributed random variables, and hence there are independent Gamma priors for the mixing variables $p_{k,n}$ over all real and imaginary components of $\textbf{x}_{k}$:

[TABLE]

A Gamma prior can also be used for $\sigma_{k}^{2}$

[TABLE]

where $\beta_{1},\beta_{2},\beta_{3}$ and $\beta_{4}$ are scale and shape priors.

It is known that

[TABLE]

and

[TABLE]

where the covariance matrix and the mean are given by

[TABLE]

and

[TABLE]

respectively. Note, the maximum of (II-B) is the posterior mean $\boldsymbol{\mu}$. For a derivation of (II-B) please see Appendix A.

Similarly to [17], the probability $\mathcal{P}(\sigma_{k}^{2},\textbf{p}_{k}|\tilde{\textbf{y}}_{k},\tilde{\textbf{x}}_{p})$ can be represented in the following form:

[TABLE]

where $\mathcal{P}(\tilde{\textbf{x}}_{p})$ is constant as fixed values are used and the second two terms on the right of are constant if $\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=1\times 10^{-4}$ as in [17]. Therefore, maximising $\mathcal{P}(\sigma_{k}^{2},\textbf{p}_{k}|\tilde{\textbf{y}}_{k},\tilde{\textbf{x}}_{p})$ is roughly equivalent to maximising $\mathcal{P}(\tilde{\textbf{y}}_{k}|\sigma_{k}^{2},\textbf{p}_{k},\tilde{\textbf{x}}_{p})$. This can be achieved by a type 2 maximisation of its logarithm, which is given by (please see Appendix B):

[TABLE]

where $\textbf{B}=(\sigma_{k}^{2}\textbf{I}+\tilde{\textbf{A}}\textbf{P}_{k}^{-1}\tilde{\textbf{A}}^{T})^{-1}$ and $\textbf{C}=\textbf{P}_{k}-\textbf{P}_{k}^{T}\boldsymbol{\Sigma}\textbf{P}_{k}$.

This is now differentiated with respect to $p_{k,n}$ and $\sigma_{k}^{-2}$ to obtain the update expressions

[TABLE]

where $\gamma_{n}=1-p_{k,n}\Sigma_{nn}$, $\Sigma_{nn}$ is the $n^{th}$ diagonal element of $\boldsymbol{\Sigma}$ and

[TABLE]

For the derivation of (26) and (27) please see Appendix C.

The maximisation is then achieved by iteratively finding $\boldsymbol{\Sigma}$ and $\boldsymbol{\mu}$, followed by $p_{k,n}^{new}$ for $n=1,...,N$ and $\sigma_{k,new}^{2}$ until a convergence criterion is met [16, 17]. In other words, the new estimates for the noise variance and precision hyperparameters found from (27) and (26) are then used in (22) and (23) to find new estimates of the covariance matrix and mean of the distribution in (II-B). Note that when $\tilde{\textbf{x}}_{p}=[0,0,...,0]^{T}$ the update expressions match those used by the traditional RVM.

The final estimate of the received signals is then given by

[TABLE]

where $\sigma^{2}_{k,opt}$ and $\textbf{P}_{k,opt}=\text{diag}(p_{k,opt,1},p_{k,opt,2},...,p_{k,opt,2N})$ are the result of optimising the noise estimate and hyperparameters, respectively. Now $\tilde{\textbf{x}}_{k,opt}$ can be used to reconstruct the estimated signals as

[TABLE]

where $n=1,2,...,N$.

The thresholding scheme in [20] can then be applied to keep the $\tilde{L}$ most significant signals. To do this find the total energy content of the estimated received signals and then sort them. A threshold value, $\eta$, is then defined as a percentage of the energy content that is to be retained. Starting with the most significant estimated signal, the estimated signals are summed until the threshold is reached and the remaining signals are then set to be equal to 0. The remaining non-zero valued signals then give the DOA estimates and $\tilde{L}$ is an estimate of the number of far field signals impinging on the array.

II-C Bayesian Compressive Sensing Kalman Filter

In order to track the changes in the DOA estimates at each time snapshot the modified RVM based DOA estimation procedure detailed above is combined with a Bayesian KF, giving a BCSKF for DOA estimation [27]. The signal model described above is again used along with the prediction

[TABLE]

and update steps

[TABLE]

of the BCSKF. Here, $k|k-1$ indicates prediction at time instance $k$ given the previous measurements and $\boldsymbol{\Delta}\textbf{x}$ is determined by the assumed DOA change. Note, $\boldsymbol{\Delta}\textbf{x}$ is fixed by the predetermined constant motion rather than being a random noise term. For example, if the angular range is sampled every $1^{\circ}$ and the DOA is assumed to increase by $2^{\circ}$ then $\boldsymbol{\Delta}\textbf{x}$ will be selected to increase the index of the non-zero valued entries in $\tilde{\textbf{x}}_{k-1|k-1}$ by two to give the index of the non-zero valued entries in $\tilde{\textbf{x}}_{k|k-1}$.

At each time snapshot it is necessary to estimate the noise variance and hyperparameters in order to evaluate the prediction and update steps of the BCSKF. This is done by considering the log likelihood function given by

[TABLE]

which can be optimised by following the procedure described in Section II-B. In other words we apply the modified RVM framework to $\tilde{\textbf{y}}_{e,k}$, using the KF prediction $\tilde{\textbf{x}}_{k|k-1}$ as the expected estimate values $\tilde{\textbf{x}}_{p}$.

It is worth noting that the continued accuracy of the proposed BCSKF relies on the accuracy of the initial estimate and the parameter values selected. If the initial estimate (made using the framework described in Section II-B and $\tilde{\textbf{x}}_{p}=[0,0,...,0]^{T}$) of the received signals is accurate and sparse, then the priors that are enforced will ensure this continues to be the case. However, an inaccurate initial DOA estimate or poorly matched expected DOA change can lead to the introduction of inaccuracies in subsequent estimates. Similarly, if the initial estimate of the received signals is not sparse then subsequent estimates are likely to not be sparse. As a result, care should be taken when choosing the initial parameter values and determining the likely DOA change.

II-D Gibbs Sampling for DOA Estimation

The method described in the previous sections based on a BCSKF with a modified RVM required the use of prior knowledge of the predicted change in DOA. However, in practice this may not always be known, making it important to have an alternative method that can still give improved accuracy for the endfire region.

This work proposes using a sparsity prior which is given as a combination of a point mass concentrated at zero (Dirac delta function) and a zero mean Gaussian distribution, [31, 32, 33], giving

[TABLE]

where $\tilde{\textbf{z}}_{k}=[\textbf{z}_{k}^{T},\textbf{z}_{k}^{T}]^{T}$ and $\textbf{z}_{k}=[z_{k,1},z_{k,2},...,z_{k,N}]^{T}$.

Note, $\tilde{z}_{k,n}$ is the indicator variable for $\tilde{x}_{k,n}$ and determines which of the two components in (33) is selected. When $\tilde{z}_{k,n}~{}=~{}0$, the value of $\tilde{x}_{k,n}$ is determined solely by the point mass concentrated at zero. As a result, $\tilde{x}_{k,n}=0$ and sparsity is explicitly introduced. Alternatively, when $\tilde{z}_{k,n}=1$ the value of $\tilde{x}_{k,n}$ is determined by the Gaussian distribution allowing a non-zero valued estimate. The repetition of $\textbf{z}_{k}$ in $\tilde{\textbf{z}}_{k}$ means that the same indicator variable is used for both the real and imaginary parts of each entry in $\textbf{x}_{k}$.

This indicator value can also be used to address the endfire accuracy problem by selecting the value of $z_{k,n}=0$ if $|n-i|>j$. Here $i$ is the index of the closet non-zero valued estimate from the previous time snapshot and $j$ defines a maximum allowed change in the DOA estimate. Only $n=1,2,...,N$ is considered to get the entries for $\textbf{z}_{k}$, with $\tilde{\textbf{z}}_{k}$ then being found as previously stated.

This leaves the following

[TABLE]

where $z_{k,n}^{1}$ and $z_{n}^{2}$ are defined by the following Beta distributions

[TABLE]

In order to enforce zero-valued estimates when $|n-i|>j$, it is necessary to select $\beta_{5}^{2}$ and $\beta_{6}^{2}$ to ensure a zero-valued $z_{k,n}$ is preferred. However, when $|n-i|\leq j$ it is necessary to choose $\beta_{5}^{1}$ and $\beta_{6}^{1}$ so that the chances of $z_{k,n}=0$ and $z_{k,n}=1$ are equal. This can be achieved by

[TABLE]

where $\beta_{5}=\beta_{6}=1$.

The posterior distribution of $\tilde{\textbf{x}}_{k}$ can be written as [33]

[TABLE]

Now also define $\tilde{\textbf{A}}_{n}$ as being the entries in $\tilde{\textbf{A}}$ relating to the index $n$ and $\tilde{\textbf{A}}_{-n}$ are the entries of $\tilde{\textbf{A}}$ excluding the entries relating to index $n$ (and similarly for $\tilde{\textbf{x}}_{k}$). Then as per Appendix D this gives

[TABLE]

where $p_{0}=1/\sigma_{k}^{2}$ and $\tilde{\textbf{y}}_{k,n}=\tilde{\textbf{y}}_{k}-\tilde{\textbf{A}}_{-n}\tilde{\textbf{x}}_{k,-n}$.

There are two further posterior distributions that have to be considered. That is the distributions for $p_{k,n}$ and $p_{0}$ which are given by

[TABLE]

and

[TABLE]

respectively. Note, in (42) $\textbf{x}_{k-1,nj}$ gives the entries within $\textbf{x}_{k-1}$ that have an index within the distance $j$ of index $n$. By using x rather than $\tilde{\textbf{x}}$ to find $\textbf{x}_{k-1,nj}$ it guarantees the same value of $||\textbf{x}_{k-1,nj}||_{2}^{2}$ and $||\textbf{x}_{k-1,nj}||_{0}$ for both the real and imaginary components.

As a result the Gibbs sampling steps are as detailed below:

Sample $\tilde{x}_{k,n}$ from $\mathcal{P}(\tilde{x}_{k,n}|\tilde{\textbf{y}}_{k},\tilde{\textbf{x}}_{k,-n},\sigma_{k}^{2},\textbf{p}_{k},\tilde{\textbf{z}}_{k})$. 2. 2.

Sample $p_{k,n}$ from $\mathcal{P}(p_{k,n}|\textbf{x}_{k-1})$. 3. 3.

$\textbf{if }n\leq N$ then Sample $z_{k,n}^{1}$ from $\mathcal{P}(z_{k,n}^{1}|\beta_{5},\beta_{6})$, $\textbf{else }z_{k,n}^{1}=z_{k,n-N}^{1}$. 4. 4.

$\textbf{if }n\leq N$ then Sample $z_{k,n}^{2}$ from $\mathcal{P}(z_{k,n}^{2}|\beta_{5},\beta_{6})$, $\textbf{else }z_{k,n}^{2}=z_{k,n-N}^{2}$. 5. 5.

Sample $p_{0}$ from $\mathcal{P}(p_{0}|\tilde{\textbf{y}}_{k},\tilde{\textbf{x}}_{k})$.

These steps are done for each of the $T$ iterations of the Gibbs sampler, where the first $T_{BI}$ iterations are the burn-in iterations. The final estimate of the received array signals is then given by the mean values of the final $T-T_{BI}$ iterations [32]. The DOA estimate can then be found using the previously described thresholding scheme (see II-B), with the remaining non-zero valued estimates corresponding to the DOA estimates.

Note, the performance of this method will again heavily depend on the accuracy of the first estimate. As a result, it is possible to use the traditional BCS DOA estimation method (Section II-B with $\tilde{\textbf{x}}_{p}=[0,0,...,0]^{T}$) to ensure an as accurate as possible intial estimate at the first time snapshot. The proposed Gibbs sampling based method can then be used at the subsequent time snapshots to get the next DOA estimate.

III Performance Evaluation

In this section a comparison in performance of the proposed methods and the traditional RVM based BCSKF method will be made over five example scenarios, under the same test conditions. Firstly, an example is considered where the initial DOA starts outside of the endfire region and then moves into it. Secondly, an example is given where the DOA remains out of the endfire region. In the third scenario the initial DOAs and the signal values are randomly generated. Then the evaluation will also consider the scenario where there is a mismatch between the actual and assumed change in DOA. Finally, the evaluation will consider a random change in DOA at each time snapshot. This means that $\Delta\textbf{x}$ which is selected for the modified RVM based BCSKF will not be a true reflection of how the DOA actually changes for the last two examples.

Note, the term traditional RVM based BCSKF method means the entries of $\textbf{P}_{k}$ in the prediction step of the BCSKF are found using the RVM optimisation method as detailed in [16, 17]. In other words this is the method detailed in Section II-B with $\tilde{\textbf{x}}_{p}=\boldsymbol{0}$. All of the examples are implemented in Matlab on a computer with an Intel Xeon CPU E3-1271 (3.60GHz) and 16GB of RAM.

The performance of each method will be measured using the $RMSE$ in the DOA estimate. This is given by

[TABLE]

where $\theta_{l}$ is the actual DOA, $\hat{\theta}_{l}$ is the estimated DOA and $Q$ is the number of Monte Carlo simulations carried out, with $Q=100$ being used in each case. This gives a measure of the estimation accuracy and the computation time will be used as a measure of the complexity of each method.

For the Gibbs sampling method a burn-in period of 250 iterations is used and then 50 further iterations used to find the final estimate of the received array signals. When a distance of $j=5^{\circ}$ is exceeded a zero-valued estimate of the received signals is enforced in order to alleviate the endfire accuracy problem.

For all the design examples considered the selection of $\sigma_{k}^{2}~{}=~{}0.4$ as the noise variance is used, with an initial estimate of the noise variance given by $\sigma_{k,0}^{2}=0.1$. The array geometry being used is that of a ULA with $M=20$ antennas and an adjacent antenna separation of $\frac{\lambda}{2}$, where $\lambda$ is the wavelength of the signal of interest. This gives an array aperture of $9.5\lambda$. For the mutual coupling matrix a value of $D=3$ is selected, meaning that a separation of $1.5\lambda$ or greater gives negligible mutual coupling effects. The values $\rho_{1}=0.65$, $\rho_{2}=0.25$, $\phi_{1}=\pi/7$ and $\phi_{2}=\pi/10$ are then also used. Finally, in each example a single narrowband signal impinging on the array is considered, meaning $L=1$.

III-A Endfire Region

For this example the initial DOA of the signal is $\theta=20^{\circ}$, which then decreases by $1^{\circ}$ at each time snapshot. The signal value at each snapshot is set to be 1. Table I summarises the performance of the three methods for this example, with the $RMSE$ values at each time snapshot being shown in Figure 2.

Here it can be seen that there has been a significant decrease in the average $RMSE$ values for both the modified RVM method (67.51$\%$ improvement) and the Gibbs sampling based method (80.78$\%$ improvement). Overall this suggests that a more accurate estimate of the DOA is possible. It is worth noting that there has still been an increase in the RMSE for the modified RVM based approach in the endfire region of the angular range. However, this has come much later on the than for the traditional RVM based approach (indicating a degradation in performance for a smaller angular range) and the maximum $RMSE$ value reached is lower (indicating the degradation is less severe).

These improvements have come at the cost of an increased computation time in both instances. For the modified RVM method this increase is insignificant as the average computation time is still less than one second. The increase for the Gibbs sampling based method is larger, illustrating an increase in computational complexity. However, it is worth remembering that this increase has resulted in a more accurate DOA estimate being achieved without prior knowledge about what the change in DOA will be.

In this instance the results suggest that one of the two proposed methods should be used when the estimated DOA approaches the endfire region of the array. If the change in DOA is known in advance and computational complexity is a primary concern then the modified RVM based method is the most suitable (a more accurate estimate can be achieved without a large increase in computation time). However, when this information is not available, or computational complexity is not a concern, it is possible to get a significant improvement in accuracy (at the cost of computation time) using the Gibbs sampling based method.

In the previous simulation an adjacent antenna separation of $\lambda/2$ is used as it is known that this is the largest separation that can be used while still avoiding a degraded performance due to the introduction of grating globes [34]. However, an example of what the relative performance of the methods is when a smaller adjacent antenna separation will now be considered.

In this instance an adjacent antenna separation of $\lambda/4$ is selected. As the array aperture is kept constant (to allow a fair comparison between adjacent antenna separation sizes) this means the number of antennas is given by $M=39$. This also means a value of $D=9$ is required to keep the same distance limits on mutual coupling occurring. The values of $\rho_{i}$ and $\phi_{i}$ are then selected to be uniformly spread over the range of 0.65 to 0.25 and $\pi/7$ to $\pi/10$, respectively. The remaining parameters are kept constant and the same test scenario as for the previous example is used.

The performances of each of the methods in this instance are summarised in Table II and Figure 3, respectively. Here it can be seen that the larger number of antennas used has resulted in a lower average $RMSE$ values for all three of the methods. In this instance only the modified RVM method has performed better that the traditional RVM based method when comparing average $RMSE$ values (decrease in average $RMSE$ of 61.09$\%$). However, by looking at the maximum $RMSE$ values it can be seen that the largest estimation error possible with the traditional RVM based method is larger than that for the Gibbs sampling based method ($16.42^{\circ}$ as compared to $3.71^{\circ}$).

It is worth noting that such an array configuration is unlikely to be used in practice. This is due to the costs associated with the number of antennas required. As a result, the remaining examples will stick to the adjacent antenna separation of $\lambda/2$ and associated parameters previously defined.

III-B Non-Endfire Region

For this example the initial DOA is $\theta=100^{\circ}$ with the DOA increasing by $1^{\circ}$ at each time snapshot, with the signal value remaining constant at -1. The performance of the three methods is summarised in Table III, with the $RMSE$ values illustrated in Figure 4.

In this instance it can be seen that there has not been as large an increase in $RMSE$ for the traditional RVM method, as the DOA does not enter the endfire region. However, both the modified RVM and Gibbs sampling based methods have managed to achieve improvements in average $RMSE$ values of 36.42$\%$ and $50.00\%$, respectively. For the Gibbs sampling based method this comes at the expenses of an increase in computation time, whereas the time for the modified RVM based method is comparable to the traditional RVM based method. As with the previous test scenario this would suggest that the modified RVM based method should be used when the expected DOA change information is available and the Gibbs sampling based method when this is not the case, or when computational complexity is not a major concern.

III-C Random Initial DOA

Next consider the case where the initial DOA is randomly chosen from the entire angular range and increased by $1^{\circ}$ at each time snapshot. The signal value is randomly selected as $\pm 1$ for each simulation and remains constant as the DOA changes.

Table IV and Figure 5 summarise the performance of the various methods in this instance. Again, it can be seen that the modified RVM based method has offered improvements in terms of $RMSE$(73.39$\%$), without a significant increase in computation time. The Gibbs sampling based method has also given an estimation accuracy improvement (76.12$\%$) and has even outperformed the modified RVM based method, without prior knowledge of how the DOA was going to change. However, this has come at the expense of an increased computation time.

III-D Mismatched Actual and Assumed DOA Change

This subsection compares the performances of the estimation methods for two situations where the actual change in DOA is not known. First, consider the case where there is an initial DOA of $\theta=100^{\circ}$ which increases by $1^{\circ}$ for 9 time snapshots before decreasing by $1^{\circ}$ for the remaining time snapshots. In this instance, assume a constant signal value of 1 throughout.

The performance comparison is now made between the traditional RVM based method, the modified RVM based method with the assumed DOA change set to a constant increase of $1^{\circ}$, the modified RVM based method with the assumed DOA change set to a constant decrease of $1^{\circ}$ and the Gibbs sampling based method. The performances for each are summarised in Figure 6 and Table V, respectively.

In this instance the average $RMSE$ values suggest a comparable performance in terms of estimation accuracy between the traditional RVM based methods and the two modified RVM based examples. This can be explained by the fact that for both of the modified RVM based examples, the assumed DOA change does not match the actual DOA changes for the entire time range which means the same improvements as for the previous scenario can no longer be guaranteed. Figure 6 highlights this in the results for the two modified RVM examples. It demonstrates that with an assumed increasing DOA the modified RVM offers some initial improvements, while the performance is significantly degraded when the DOA starts to decrease again. On the other hand, the example with an assumed decreasing DOA performs worse than the traditional RVM based method initially and then offers significant improvements when the actual DOA also starts to decrease.

It can also be seen that for the Gibbs sampling based method there has been an improvement in DOA estimation accuracy. In terms of average $RMSE$ values this is a decrease of 68.60$\%$, which has been achieved without any knowledge of how the DOA was going to change. However, there is again an increase in the computation time.

To illustrate how a larger mismatch between actual and assumed DOA changes effects the performance of the modified RVM based method now consider an example where the actual DOA is increasing by $1^{\circ}$ in each snapshot, while the assumed change is a decrease of $3^{\circ}$. Here, the initial DOA is $100^{\circ}$, with a constant signal value of -1. The $RMSE$ values for the methods are shown in Figure 7 and summarised in Table VI along with the computation times.

Here it can be seen that the Gibbs sampling based method has offered an $86.36\%$ improvement in average $RMSE$ as compared to the traditional RVM based method. There has again been a significant increase in the computational complexity.

For the modified RVM based method the average $RMSE$ values suggests that there has been an improvement in estimation accuracy. However, this is smaller than when the actual and assumed DOA changes match. It is also unlikely that this improvement would be obtained in every scenario the the method could be applied to. From looking at Figure 7 we can see that the traditional and modified RVM based methods are showing comparable performance for the just over half of the time frame considered. This is the relative performance that would be expected in the majority of cases.

III-E Random Changes in Direction of Arrival

Finally, consider the example where the initial signal value is assumed to be equal to 1 and the initial DOA is chosen to be $100^{\circ}$. The actual DOA is then allowed to randomly change by up to $\pm 3^{\circ}$ for each time snapshot. For the modified RVM method assume that the actual DOA change is an increase of $3^{\circ}$. This gives the results as summarised in Table VII and Figure 8.

It can be seen that the Gibbs sampling based method has again outperformed the modified RVM based method in terms of estimation accuracy, due to the fact that no prior knowledge of how the DOA will change is required. As compared to the traditional RVM based method there has been an improvement in $RMSE$ of 78.81$\%$. However, as is expected this is at the cost of computation time.

It is also worth noting that the average $RMSE$ values suggest that the modified RVM and traditional RVM have offered a comparable performance. This is due to the fact that the assumption of how the DOA will change is not valid, meaning the modified RVM no longer offers any improvements. Therefore, in this situation the Gibbs sampling based method would be the best to use, assuming computational complexity is not the main motivating factor.

IV Conclusions

This paper has proposed two novel approaches for the estimation of a dynamic direction of arrival using uniform linear arrays with mutual coupling. The first approach is a Bayesian compressed Kalman filter with a modified relevance vector machine, where the traditional sparsity assumption is replaced by an assumption that the estimated signals will instead match predicted signal values. This results in the derivation of a new posterior probability density function of the received signals and the expression for the related marginal likelihood function. The second proposed approach is a Gibbs sampling approach, where sparsity is explicitly enforced if there is a large difference between the previous direction of arrival estimate and the angle currently being considered. The proposed approaches will be particularly useful when applied to the problem of dynamic direction of arrival estimation in the endfire region of antenna arrays. Such problems can arise in numerous application areas such as in communications and surveillance.

An extensive performance evaluation is provided and shows that both of the proposed approaches outperform the traditional relevance vector machine based Bayesian compressive sensing Kalman filter in terms of mean root mean square error values, by up to 73.39$\%$ for the modified relevance vector machine based method and 86.36$\%$ for the Gibbs sampling based method.

Appendix

IV-A Derivation of Posterior Distribution

Bayes’ rule gives

[TABLE]

where $\mathcal{P}(\tilde{\textbf{y}}_{k}|\tilde{\textbf{x}}_{k},\sigma_{k}^{2})$ and $\mathcal{P}(\tilde{\textbf{x}}_{k}|\textbf{p}_{k},\tilde{\textbf{x}}_{p})$ are known from (16) and (17), respectively.

Now following the method suggested in [17] carry out the multiplication on the right hand side of (IV-A), collect terms in $\tilde{\textbf{x}}_{k}$ in the exponential and complete the square.

[TABLE]

where $\boldsymbol{\Sigma}$ and $\boldsymbol{\mu}$ are given by (22) and (23), respectively. This then gives the posterior distribution as (II-B), with the remaining exponential terms

[TABLE]

IV-B Derivation of Marginal Likelihood

From (IV-A) the following is known:

[TABLE]

meaning the term in the exponential will be (47) where

[TABLE]

Therefore the exponential term is given by

[TABLE]

The term outside of the exponential is given by

[TABLE]

This gives the marginal likelihood as

[TABLE]

where B and C are defined as in Section II-B. The log likelihood is then given by

[TABLE]

Using the Woodbury matrix inversion identity gives

[TABLE]

which means

[TABLE]

Also, we know that $\textbf{P}_{k}^{T}=\textbf{P}_{k}$ as $\textbf{P}_{k}$ is a real valued diagonal matrix. This means

[TABLE]

which then gives the log likelihood function in (25).

IV-C Derivation of Update Expressions for Modified RVM

Firstly, differentiating (25) with respect to $p_{k,n}$ gives

[TABLE]

and equating to zero gives

[TABLE]

which leads to (26).

Now collect the terms with $\sigma_{k}^{2}$ in to give

[TABLE]

and then define $\tau=\sigma_{k}^{-2}$ giving

[TABLE]

Now differentiate (60) with respect to $\tau$ and equate to zero to give

[TABLE]

where $\text{tr}(\cdot)$ indicates the trace. As $\text{tr}(\boldsymbol{\Sigma}\tilde{\textbf{A}}^{T}\tilde{\textbf{A}})$ can be written as $\tau^{-1}\sum\limits_{n}\gamma_{n}$ giving

[TABLE]

which in turn gives (27).

IV-D Derivation of (38), (39), (40) and (41)

From (II-D) it is known that

[TABLE]

If we then combine the exponential terms in the second term in (IV-D) we get

[TABLE]

Completing the square gives

[TABLE]

where $\hat{p}_{k,n}$ and $\hat{\mu}_{k,n}$ are given by (39) and (40), respectively.

In order to complete the expression given in (38) it is now necessary to get a new indicator variable, $\hat{z}_{k,n}$ for the new posterior distribution for $\tilde{x}_{k,n}$. To do this, assume that

[TABLE]

Thus giving us (41), allowing us to write the posterior distribution for $x_{k,n}$ in the form given in (38).

Acknowledgments

We appreciate the support of the UK Engineering and Physical Sciences Research Council (EPSRC) via the project Bayesian Tracking and Reasoning over Time (BTaRoT) grant EP/K021516/1. We acknowledge the anonymous reviewers’ suggestions that have helped improve this work and would like to thank the associate editor for handling the review of our paper.