Stochastic Maximum Likelihood Direction Finding in the Presence of Nonuniform Noise Fields

Abstract

The maximum likelihood (ML) technique plays an important role in direction-of-arrival (DOA) estimation. In this paper, we employ and design the expectation–conditional maximization either (ECME) algorithm, a generalization of the expectation–maximization algorithm, for solving the ML direction finding problem of stochastic sources, which may be correlated, in unknown nonuniform noise. Unlike alternating maximization, the ECME algorithm updates both the source and noise covariance matrix estimates by explicit formulas, and can guarantee that both estimates are positive semi-definite and definite, respectively. Thus, the ECME algorithm is computationally efficient and operationally stable. Simulation results confirm that the ECME algorithm can efficiently obtain the ML based DOA estimate of each stochastic source.

1. Introduction

Two source signal models are widely used in Cramer–Rao lower bound (CRLB) and maximum likelihood (ML) direction finding, i.e., the deterministic signal model where signals are deterministic and unknown, and the stochastic signal model where signals are Gaussian. For example, various CRLBs using both models have been derived [1,2,3,4,5,6,7]. However, ML direction finding generally involves high-dimensional search algorithms for both models, which causes a significant increase in computational complexity.

In order to reduce the computational complexity, two classic methods have been developed, i.e., alternating maximization (AM)-type [8] and expectation–maximization (EM)-type [9,10,11,12] algorithms. These two methods are first applied under uniform Gaussian noise, which decreases the number of parameters and simplifies the problem. However, the uniform noise model is unrealistic in many situations, and nonuniform noise has been considered in numerous papers [4,13,14,15,16,17,18,19]. In nonuniform noise, the covariance matrix still keeps a diagonal structure, but the diagonal elements are no longer identical, which hinders direction-of-arrival (DOA) estimation. To tackle the problem of direction finding in unknown nonuniform noise, diverse subspace separation approaches based on the subspace technique have been proposed in the literature [13,14,15,16,17,18].

For obtaining ML based solutions, AM- and EM-type algorithms have also been applied to this problem. However, the AM-type algorithms usually require high-dimensional numerical searches due to the noise nonuniformity at each iteration [4,19], which leads to heavy computational burdens. Moreover, when considering Gaussian source signals, the AM algorithm presented in [19] has one severe shortcoming: the source and noise covariance matrix estimates cannot be guaranteed to be positive semi-definite and definite, respectively. To this end, we have designed several computationally efficient EM-type algorithms in [20] that only need low-dimensional (one or two-dimensional) numerical searches at every iteration. In these EM-type algorithms using the stochastic signal model, however, the sources must be uncorrelated. This restricts the use of stochastic ML direction finding in some situations, e.g., multipath conditions. As a consequence, efficient algorithms are urgently needed to address this issue.

In this paper, we employ and design the expectation–conditional maximization either (ECME) algorithm [21], a generalization of the EM algorithm, for solving the ML direction finding problem of stochastic sources, which may be correlated, in unknown nonuniform noise. Unlike the AM algorithm in [19], the ECME algorithm updates both the source and noise covariance matrix estimates by explicit formulas and can guarantee that both estimates are positive semi-definite and definite, respectively. Thus, the ECME algorithm is computationally efficient and operationally stable. Simulation results confirm the effectiveness of the algorithm.

The rest of this paper is outlined as follows: In Section 2, we formulate the stochastic ML direction finding problem in unknown nonuniform noise. In Section 3 and Section 4, we design the ECME algorithm and provide simulation results to show its effectiveness, respectively. Lastly, we conclude this paper in Section 5.

2. Problem Statement

For simplicity, let a uniformly spaced linear array of W sensors receive the plane waves impinging from V $(V<W)$ narrow-band sources of wavelength $\iota$. The distance between any adjacent sensors is $\iota /2$. We denote the direction associated with the vth source by ${\beta }_v\in (0,\pi )$ $\left(\mathrm{radian}\right)$, and write the received signal as

$$\mathbf{r}\left(t\right)={\sum }_{v=1}^V\mathbf{a}\left({\beta }_v\right)k_v\left(t\right)+\mathbf{j}\left(t\right)=\mathbf{A}\left(\mathit{\beta }\right)\mathbf{k}\left(t\right)+\mathbf{j}\left(t\right),$$

(1)

where $\mathbf{a}\left({\beta }_v\right)={[1a_v\cdots a_v^{W-1}]}^T$, $a_v=exp\left(-j\pi cos\left({\beta }_v\right)\right)$, ${[·]}^T$ denotes transposition, $j=\sqrt{-1}$, $k_v\left(t\right)$ is the signal with respect to the vth source, and $\mathbf{j}\left(t\right)$ means nonuniform complex Gaussian noise of zero mean and covariance $\mathbf{Q}$, i.e., $\mathbf{j}\left(t\right)\sim \mathcal{CN}(\mathbf{0},\mathbf{Q})$. Here, $\mathbf{Q}$ is diagonal and expressed as $\mathbf{Q}=\mathrm{diag}\left\{\mathit{\delta }\right\}$, where $\mathit{\delta }={[{\delta }_1\cdots {\delta }_W]}^T>\mathbf{0}$ and $\mathbf{Q}$ is positive definite, i.e., $\mathbf{Q}\succ {\mathbf{0}}_W$ (${\mathbf{0}}_W$ is the $W\times W$ zero matrix). Furthermore, if ${\delta }_1=\cdots ={\delta }_W=\delta >0$, $\mathbf{Q}=\delta {\mathbf{I}}_W$ (${\mathbf{I}}_W$ is the $W\times W$ identity matrix), which makes the noise uniform. In (1), $\mathbf{A}\left(\mathit{\beta }\right)=[\mathbf{a}\left({\beta }_1\right)\cdots \mathbf{a}\left({\beta }_V\right)]$ is the array manifold matrix, $\mathit{\beta }={[{\beta }_1\cdots {\beta }_V]}^T\in \mathbf{\Gamma }$ with $\mathbf{\Gamma }={(0,\pi )}^V$, and $\mathbf{k}\left(t\right)={[k_1\left(t\right)\cdots k_V\left(t\right)]}^T$. For notational convenience, we use $\mathbf{A}$ instead of $\mathbf{A}\left(\mathit{\beta }\right)$ hereafter.

We consider Gaussian source signals, which may be correlated, and have $\mathbf{k}\left(t\right)\sim \mathcal{CN}(\mathbf{0},\mathbf{O})$, where $\mathbf{O}$ is the source covariance matrix and positive semi-definite, i.e., $\mathbf{O}⪰{\mathbf{0}}_V$. Let the sources be uncorrelated with the noise, such that

$$\begin{array}{c}\hfill \mathbf{r}\left(t\right)\sim \mathcal{CN}(\mathbf{0},\mathbf{G}),\mathbf{G}=\mathbf{A}\mathbf{O}{\mathbf{A}}^H+\mathbf{Q}\succ {\mathbf{0}}_W,\end{array}$$

where ${[·]}^H$ is conjugate transposition. On this foundation, the log-likelihood function (LLF) of L statistically independent snapshots can be formulated as

$$\begin{array}{c}\hfill \mathcal{J}(\mathit{\beta },\mathit{\rho },\mathit{\delta })={\sum }_{t=1}^{L}logp\left(\mathbf{r}\left(t\right);\mathit{\beta },\mathit{\rho },\mathit{\delta }\right)=f-L\left(log\left|\mathbf{G}\right|+\mathrm{trace}\left[{\mathbf{G}}^{-1}\widehat{\mathbf{R}}\right]\right),\end{array}$$

(2)

where $|·|$, $\mathrm{trace}[·]$, and ${(·)}^{-1}$ denote determinant, trace, and inversion, respectively. In (2), f is a constant, $\widehat{\mathbf{R}}=(1/L){\sum }_{t=1}^L\mathbf{r}\left(t\right){\mathbf{r}}^H\left(t\right)$ means the covariance matrix of snapshots. Moreover,

$$\begin{array}{c}\hfill \mathit{\rho }=\left({\left[\mathbf{O}\right]}_{1,1},\dots ,{\left[\mathbf{O}\right]}_{V,V},\mathrm{Re}\left\{{\left[\mathbf{O}\right]}_{1,2}\right\},\mathrm{Im}\left\{{\left[\mathbf{O}\right]}_{1,2}\right\},\dots ,\mathrm{Re}\left\{{\left[\mathbf{O}\right]}_{V-1,V}\right\},\mathrm{Im}\left\{{\left[\mathbf{O}\right]}_{V-1,V}\right\}\right),\end{array}$$

where ${\left[\mathbf{O}\right]}_{p,q}$ is the $(p,q)$th element of $\mathbf{O}$, $\mathrm{Re}\left\{a\right\}$ and $\mathrm{Im}\left\{a\right\}$ represent the real part and imaginary part of a, respectively. Consequently, the ML based DOA estimation problem is

$$\begin{array}{c}\hfill \underset{\mathit{\beta }\in \mathbf{\Gamma },\mathbf{O}⪰{\mathbf{0}}_V,\mathit{\delta }>\mathbf{0}}{max}\mathcal{J}(\mathit{\beta },\mathit{\rho },\mathit{\delta }).\end{array}$$

(3)

We assume $\mathcal{R}\left[\mathbf{A}\right]=V$, where $\mathcal{R}\left[\mathbf{A}\right]$ is the rank of $\mathbf{A}$, and can thus eliminate $\mathbf{O}$ in (3) by [22]

$$\begin{array}{ccc}\hfill \widehat{\mathbf{O}}(\mathit{\beta },\mathit{\delta })& =& {\left({\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}\right)}^{-1}{\tilde{\mathbf{A}}}^H\left(\tilde{\mathbf{R}}-{\mathbf{I}}_W\right)\tilde{\mathbf{A}}{\left({\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}\right)}^{-1}\hfill \\ & =& {\left({\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}\right)}^{-1}{\tilde{\mathbf{A}}}^H\tilde{\mathbf{R}}\tilde{\mathbf{A}}{\left({\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}\right)}^{-1}-{\left({\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}\right)}^{-1},\hfill \end{array}$$

(4)

where ${\mathbf{Q}}^{-1/2}=\mathrm{diag}\{1/\sqrt{{\delta }_1},\dots ,1/\sqrt{{\delta }_W}\}$, $\tilde{\mathbf{A}}={\mathbf{Q}}^{-1/2}\mathbf{A}$, and $\tilde{\mathbf{R}}={\mathbf{Q}}^{-1/2}\widehat{\mathbf{R}}{\mathbf{Q}}^{-1/2}$. In other words, $\mathit{\rho }$ can be estimated using the estimates of $\mathit{\beta }$ and $\mathit{\delta }$. On the basis of (4), $\mathbf{G}$ is rewritten as

$$\begin{array}{c}\hfill \mathbf{G}=\mathbf{A}\widehat{\mathbf{O}}(\mathit{\beta },\mathit{\delta }){\mathbf{A}}^H+\mathbf{Q}={\mathbf{Q}}^{1/2}\left({\mathbf{\Pi }}_{\tilde{\mathbf{A}}}\tilde{\mathbf{R}}{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}+{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}^{\perp }\right){\mathbf{Q}}^{1/2},\end{array}$$

where ${\mathbf{\Pi }}_{\tilde{\mathbf{A}}}=\tilde{\mathbf{A}}{\left({\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}\right)}^{-1}{\tilde{\mathbf{A}}}^H$ and ${\mathbf{\Pi }}_{\tilde{\mathbf{A}}}^{\perp }={\mathbf{I}}_W-{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}$. Then, Problem (3) is reduced to [22]

$$\begin{array}{c}\hfill \underset{\mathit{\beta }\in \mathbf{\Gamma },\mathit{\delta }>\mathbf{0}}{min}\mathcal{H}(\mathit{\beta },\mathit{\delta })=log\left|{\mathbf{Q}}^{1/2}\left({\mathbf{\Pi }}_{\tilde{\mathbf{A}}}\tilde{\mathbf{R}}{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}+{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}^{\perp }\right){\mathbf{Q}}^{1/2}\right|+\mathrm{trace}\left[{\left({\mathbf{\Pi }}_{\tilde{\mathbf{A}}}\tilde{\mathbf{R}}{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}+{\mathbf{\Pi }}_{\tilde{\mathbf{A}}}^{\perp }\right)}^{-1}\tilde{\mathbf{R}}\right].\end{array}$$

(5)

In particular, if the noise is uniform Gaussian noise, Problem (5) can be further reduced to [23]

$$\begin{array}{c}\hfill \underset{\mathit{\beta }\in \mathbf{\Gamma }}{min}\mathcal{G}\left(\mathit{\beta }\right)=|\mathbf{A}\widehat{\mathbf{O}}\left(\mathit{\beta }\right){\mathbf{A}}^H+\widehat{\delta }\left(\mathit{\beta }\right){\mathbf{I}}_W|,\end{array}$$

(6)

where

$$\begin{array}{ccc}\hfill \widehat{\delta }\left(\mathit{\beta }\right)& =& \mathrm{trace}\left[\left({\mathbf{I}}_W-\mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}{\mathbf{A}}^H\right)\widehat{\mathbf{R}}\right]/(W-V),\hfill \\ \hfill \widehat{\mathbf{O}}\left(\mathit{\beta }\right)& =& {\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}{\mathbf{A}}^H\left(\widehat{\mathbf{R}}-\widehat{\delta }\left(\mathit{\beta }\right){\mathbf{I}}_W\right)\mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}.\hfill \end{array}$$

Unfortunately, it is very difficult to reduce Problem (5) to some problems with fewer parameters under nonuniform Gaussian noise. Of course, applying gradient-type algorithms to search the solution of Problem (5) is computationally intensive due to the search space of dimension $W+V$ and the complexity of $\mathcal{H}(\mathit{\beta },\mathit{\delta })$.

In fact, when direct maximization over all parameters is intractable, AM can always be utilized. As stated before, the authors in [19] have presented an AM algorithm consisting of two steps at every iteration for Problem (3). Specifically, the first step obtains ${\mathit{\delta }}^{\left(d\right)}$, the estimate of $\mathit{\delta }$ at the dth iteration, by a gradient based algorithm, which is called the “modified inverse iteration algorithm” and satisfies

$$\begin{array}{c}\hfill \mathcal{J}({\mathit{\beta }}^{(d-1)},{\mathit{\rho }}^{(d-1)},{\mathit{\delta }}^{\left(d\right)})\ge \mathcal{J}({\mathit{\beta }}^{(d-1)},{\mathit{\rho }}^{(d-1)},{\mathit{\delta }}^{(d-1)}),\end{array}$$

(7)

where ${[·]}^{\left(0\right)}$ means an initial estimate. Then, the second step simultaneously obtains ${\mathit{\beta }}^{\left(d\right)}$ and ${\mathit{\rho }}^{\left(d\right)}$ by

$$\begin{array}{c}\hfill ({\mathit{\beta }}^{\left(d\right)},{\mathit{\rho }}^{\left(d\right)})=arg\underset{\mathit{\beta }\in \mathbf{\Gamma },\mathbf{O}⪰{\mathbf{0}}_V}{max}\mathcal{J}(\mathit{\beta },\mathit{\rho },{\mathit{\delta }}^{\left(d\right)}),\end{array}$$

(8)

which is solved in a separable manner, i.e.,

$$\begin{array}{ccc}\hfill {\mathit{\beta }}^{\left(d\right)}& =& arg\underset{\mathit{\beta }\in \mathbf{\Gamma }}{min}\mathcal{H}(\mathit{\beta },{\mathit{\delta }}^{\left(d\right)}),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\mathbf{O}}^{\left(d\right)}& =& \widehat{\mathbf{O}}({\mathit{\beta }}^{\left(d\right)},{\mathit{\delta }}^{\left(d\right)}).\hfill \end{array}$$

(10)

However, the AM algorithm has two drawbacks: (1) obtaining ${\mathit{\delta }}^{\left(d\right)}$ and ${\mathit{\beta }}^{\left(d\right)}$ is computationally expensive; (2) ${\mathbf{Q}}^{\left(d\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(d\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(d\right)}⪰{\mathbf{0}}_V$ cannot be guaranteed [24,25]. To efficiently obtain the ML estimate of $\mathit{\beta }$ in (3), we employ and design the ECME algorithm in the next section.

3. ECME Algorithm

Existing EM-type algorithms for stochastic ML direction finding are only applicable to uncorrelated sources [9,12,20], i.e., $\mathbf{O}$ is diagonal. In this section, we employ and design the ECME algorithm [21], a generalization of the EM algorithm, to solve Problem (3) associated with correlated sources.

3.1. Procedure

The sources in (1) may be correlated, so we choose $\mathbf{K}=\left[\mathbf{k}\right(1)\cdots \mathbf{k}(L\left)\right]$ and $\mathbf{J}=\left[\mathbf{j}\right(1)\cdots \mathbf{j}(L\left)\right]$ as augmented data. We express the augmented-data LLF as

$$\begin{array}{cc}\hfill \mathcal{M}(\mathbf{K},\mathbf{J};\mathit{\rho },\mathit{\delta })=& {\sum }_{t=1}^L\left[logp\left(\mathbf{k}\right(t);\mathsf{\rho })+logp\left(\mathbf{j}\right(t);\mathsf{\delta })\right]\\ \hfill =& h-L\left(log\left|\mathbf{O}\right|+\mathrm{trace}\left[{\mathbf{O}}^{-1}{\widehat{\mathbf{N}}}_k\right]\right)+\\ & f-L\left(log\left|\mathbf{Q}\right|+\mathrm{trace}\left[{\mathbf{Q}}^{-1}{\widehat{\mathbf{N}}}_j\right]\right),\end{array}$$

(11)

where h is a constant, ${\widehat{\mathbf{N}}}_k=(1/L){\sum }_{t=1}^L\mathbf{k}\left(t\right){\mathbf{k}}^H\left(t\right)$, and ${\widehat{\mathbf{N}}}_j=(1/L){\sum }_{t=1}^L\mathbf{j}\left(t\right){\mathbf{j}}^H\left(t\right)$. With (11), we first construct the EM algorithm [26], whose expectation and maximization steps at the dth iteration are derived below. Let $\mathcal{E}\{·\}$ and $\mathcal{D}\{·\}$ represent expectation and covariance, respectively.

3.1.1. Expectation Step

Compute the conditional expectation of the augmented-data LLF, i.e.,

$$\begin{array}{cc}\hfill \mathcal{M}\left(\mathit{\rho },\mathit{\delta };{\mathbf{\Omega }}^{(d-1)}\right)=& \mathcal{E}\left\{\mathcal{M}(\mathbf{K},\mathbf{J};\mathit{\rho },\mathit{\delta })\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}\\ \hfill =& h-L\left(log\left|\mathbf{O}\right|+\mathrm{trace}\left[{\mathbf{O}}^{-1}{\widehat{\mathbf{N}}}_k^{\left(d\right)}\right]\right)+\\ & f-L\left(log\left|\mathbf{Q}\right|+\mathrm{trace}\left[{\mathbf{Q}}^{-1}{\widehat{\mathbf{N}}}_j^{\left(d\right)}\right]\right)\end{array}$$

(12)

with ${\mathbf{\Omega }}^{(d-1)}=({\mathit{\beta }}^{(d-1)},{\mathit{\rho }}^{(d-1)},{\mathit{\delta }}^{(d-1)})$, ${\mathit{\beta }}^{(d-1)}\in \mathbf{\Gamma }$, ${\mathbf{O}}^{(d-1)}⪰{\mathbf{0}}_V$, ${\mathit{\delta }}^{(d-1)}>\mathbf{0}$, and $\mathbf{F}=\left[\mathbf{r}\right(1)\cdots \mathbf{r}(L\left)\right]$. Moreover,

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{N}}}_k^{\left(d\right)}& =& \mathcal{E}\left\{{\widehat{\mathbf{N}}}_k\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}\\ & =& {\left[{\mathbf{H}}^{(d-1)}\right]}^H\widehat{\mathbf{R}}{\mathbf{H}}^{(d-1)}+{\mathbf{O}}^{(d-1)}-{\left[{\mathbf{H}}^{(d-1)}\right]}^H{\mathbf{G}}^{(d-1)}{\mathbf{H}}^{(d-1)}⪰{\mathbf{0}}_V,\end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{N}}}_j^{\left(d\right)}& =& \mathcal{E}\left\{{\widehat{\mathbf{N}}}_j\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}\\ & =& {\mathbf{Q}}^{(d-1)}{\left[{\mathbf{G}}^{(d-1)}\right]}^{-1}\widehat{\mathbf{R}}{\left[{\mathbf{G}}^{(d-1)}\right]}^{-1}{\mathbf{Q}}^{(d-1)}+\\ & & {\mathbf{Q}}^{(d-1)}-{\mathbf{Q}}^{(d-1)}{\left[{\mathbf{G}}^{(d-1)}\right]}^{-1}{\mathbf{Q}}^{(d-1)}⪰{\mathbf{0}}_W,\end{array}$$

(14)

where ${\mathbf{H}}^{(d-1)}={\left[{\mathbf{G}}^{(d-1)}\right]}^{-1}{\mathbf{A}}^{(d-1)}{\mathbf{O}}^{(d-1)}$, the conditional distributions of $\mathbf{k}\left(t\right)$ and $\mathbf{j}\left(t\right)$ can be obtained in [27], and

$$\begin{array}{ccc}\hfill \mathcal{E}\left\{\mathbf{k}\left(t\right)\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}& =& {\left[{\mathbf{H}}^{(d-1)}\right]}^H\mathbf{r}\left(t\right),\hfill \\ \hfill \mathcal{D}\left\{\mathbf{k}\left(t\right)\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}& =& {\mathbf{O}}^{(d-1)}-{\left[{\mathbf{H}}^{(d-1)}\right]}^H{\mathbf{G}}^{(d-1)}{\mathbf{H}}^{(d-1)}⪰{\mathbf{0}}_V,\hfill \\ \hfill \mathcal{E}\left\{\mathbf{j}\left(t\right)\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}& =& {\mathbf{Q}}^{(d-1)}{\left[{\mathbf{G}}^{(d-1)}\right]}^{-1}\mathbf{r}\left(t\right),\hfill \\ \hfill \mathcal{D}\left\{\mathbf{j}\left(t\right)\mid \mathbf{F};{\mathbf{\Omega }}^{(d-1)}\right\}& =& {\mathbf{Q}}^{(d-1)}-{\mathbf{Q}}^{(d-1)}{\left[{\mathbf{G}}^{(d-1)}\right]}^{-1}{\mathbf{Q}}^{(d-1)}⪰{\mathbf{0}}_W.\hfill \end{array}$$

3.1.2. Maximization Step

Obtain ${\mathit{\rho }}^{\left(d\right)}$ and ${\mathit{\delta }}^{\left(d\right)}$ by maximizing $\mathcal{M}(\mathit{\rho },\mathit{\delta };{\mathbf{\Omega }}^{(d-1)})$ with respect to $\mathit{\rho }$ and $\mathit{\delta }$, which leads to the two parallel subproblems

$$\begin{array}{c}\hfill \underset{\mathbf{O}⪰{\mathbf{0}}_V}{min}log\left|\mathbf{O}\right|+\mathrm{trace}\left[{\mathbf{O}}^{-1}{\widehat{\mathbf{N}}}_k^{\left(d\right)}\right],\end{array}$$

(15)

$$\begin{array}{c}\hfill \underset{\mathbf{Q}\succ {\mathbf{0}}_W}{min}log\left|\mathbf{Q}\right|+\mathrm{trace}\left[{\mathbf{Q}}^{-1}{\widehat{\mathbf{N}}}_j^{\left(d\right)}\right].\end{array}$$

(16)

${\mathit{\rho }}^{\left(d\right)}$ and ${\mathit{\delta }}^{\left(d\right)}$ are simultaneously obtained by

$$\begin{array}{ccc}\hfill {\mathbf{O}}^{\left(d\right)}& =& {\widehat{\mathbf{N}}}_k^{\left(d\right)}⪰{\mathbf{0}}_V,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\delta }_w^{\left(d\right)}& =& \left\{\begin{array}{cc}{\left[{\widehat{\mathbf{N}}}_j^{\left(d\right)}\right]}_{w,w},\hfill & {\left[{\widehat{\mathbf{N}}}_j^{\left(d\right)}\right]}_{w,w}>0,\hfill \\ {\delta }_w^{(d-1)}/2,\hfill & {\left[{\widehat{\mathbf{N}}}_j^{\left(d\right)}\right]}_{w,w}=0,\hfill \end{array}\right.\forall w.\hfill \end{array}$$

(18)

From (17) and (18), we have the monotonicity of generalized EM algorithms [26], i.e.,

$$\begin{array}{c}\hfill \mathcal{J}\left({\mathit{\beta }}^{(d-1)},{\mathit{\rho }}^{\left(d\right)},{\mathit{\delta }}^{\left(d\right)}\right)\ge \mathcal{J}\left({\mathit{\beta }}^{(d-1)},{\mathit{\rho }}^{(d-1)},{\mathit{\delta }}^{(d-1)}\right).\end{array}$$

(19)

Obviously, ${\mathit{\beta }}^{\left(d\right)}$ is not obtained at the dth iteration of the EM algorithm.

3.1.3. Conditional Maximization Step

In order to obtain ${\mathit{\beta }}^{\left(d\right)}$, we now add a conditional maximization step at this iteration. Considering the following monotonicity:

$$\begin{array}{c}\hfill \mathcal{J}\left({\mathit{\beta }}^{\left(d\right)},{\mathit{\rho }}^{\left(d\right)},{\mathit{\delta }}^{\left(d\right)}\right)\ge \mathcal{J}\left({\mathit{\beta }}^{(d-1)},{\mathit{\rho }}^{\left(d\right)},{\mathit{\delta }}^{\left(d\right)}\right),{\mathit{\beta }}^{\left(d\right)}\in \mathbf{\Gamma },\end{array}$$

(20)

We can design this step as

$$\begin{array}{c}\hfill {\mathit{\beta }}^{\left(d\right)}=arg\underset{\mathit{\beta }\in \mathbf{\Gamma }}{max}\mathcal{J}\left(\mathit{\beta },{\mathit{\rho }}^{\left(d\right)},{\mathit{\delta }}^{\left(d\right)}\right),\end{array}$$

(21)

or use a gradient-type algorithm to obtain ${\mathit{\beta }}^{\left(d\right)}$ based on (20), e.g., Algorithm 1 in the next section. Due to the additional step unrelated to augmented data, the above EM algorithm becomes the ECME algorithm [21].

Algorithm 1 Steepest Descent-Based DOA Estimation

1: $f\left(\mathit{\beta }\right)=-\mathcal{J}\left(\mathit{\beta },{\mathit{\rho }}^{\left(d\right)},{\mathit{\delta }}^{\left(d\right)}\right)/L$, initialize $\mathit{\beta }={\mathit{\beta }}^{(d-1)}\in \mathbf{\Gamma }$. 2: while  ${\parallel \nabla f\left(\mathit{\beta }\right)\parallel }_2>0.001$ do 3:    $t_v=\left\{\begin{array}{cc}-(\pi -{\beta }_v)/f_v^{\prime }\left(\mathit{\beta }\right),\hfill & f_v^{\prime }\left(\mathit{\beta }\right)<0,\hfill \\ {\beta }_v/f_v^{\prime }\left(\mathit{\beta }\right),\hfill & f_v^{\prime }\left(\mathit{\beta }\right)>0,\hfill \\ \infty ,\hfill & f_v^{\prime }\left(\mathit{\beta }\right)=0,\hfill \end{array}\right. \forall v.$ 4:    $t=0.1\times min\{t_1,\cdots ,t_V\}.$ 5:    while $f\left(\mathit{\beta }-t\nabla f\left(\mathit{\beta }\right)\right)>f\left(\mathit{\beta }\right)-0.3t{\parallel \nabla f\left(\mathit{\beta }\right)\parallel }_2^2$ do 6:      $t=0.5t$. 7:    end while 8:    $\mathit{\beta }=\mathit{\beta }-t\nabla f\left(\mathit{\beta }\right)\in \mathbf{\Gamma }$. 9: end while 10: ${\mathit{\beta }}^{\left(d\right)}=\mathit{\beta }$.

3.2. Stability and Complexity

The stable operation of the ECME algorithm requires ${\mathbf{Q}}^{\left(d\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(d\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(d\right)}⪰{\mathbf{0}}_V$ for $d\ge 0$, so we give the following proposition.

Proposition 1.

In the ECME algorithm, ${\mathbf{Q}}^{\left(d\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(d\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(d\right)}⪰{\mathbf{0}}_V$ for $d\ge 1$ if ${\mathbf{Q}}^{\left(0\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(0\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(0\right)}⪰{\mathbf{0}}_V$.

Proof. 

We utilize the mathematical induction method. If ${\mathbf{Q}}^{\left(u\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(u\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(u\right)}⪰{\mathbf{0}}_V$, we have ${\mathbf{G}}^{\left(u\right)}={\mathbf{A}}^{\left(u\right)}{\mathbf{O}}^{\left(u\right)}{\left[{\mathbf{A}}^{\left(u\right)}\right]}^H+{\mathbf{Q}}^{\left(u\right)}\succ {\mathbf{0}}_W$, which leads to ${\widehat{\mathbf{N}}}_j^{(u+1)}⪰{\mathbf{0}}_W$ in (14) and then in (18) ${\delta }_w^{(u+1)}>0,\forall w$, i.e., ${\mathbf{Q}}^{(u+1)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{(u+1)}>\mathbf{0}$). Furthermore, ${\mathbf{O}}^{(u+1)}={\widehat{\mathbf{N}}}_k^{(u+1)}⪰{\mathbf{0}}_V$ is straightforward in (13). The proof is completed. □

Proposition 1 indicates that, when ${\mathbf{Q}}^{\left(0\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(0\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(0\right)}⪰{\mathbf{0}}_V$ in the ECME algorithm, ${\mathit{\rho }}^{\left(d\right)}$ and ${\mathit{\delta }}^{\left(d\right)}$ obtained at the dth iteration are in the parameter spaces, respectively. Hence, the ECME algorithm is operationally stable.

Since ${\mathit{\rho }}^{\left(d\right)}$ and ${\mathit{\delta }}^{\left(d\right)}$ are obtained via the explicit formulas in (17) and (18), the computational complexity of the ECME algorithm is dominated by obtaining ${\mathit{\beta }}^{\left(d\right)}$ in (20). Compared with the AM algorithm in [19], the ECME algorithm is, thus, computationally efficient.

3.3. Limit Point

According to [21,28], we know that the ECME algorithm satisfies certain regularity conditions and always converges to a stationary point of $\mathcal{J}\mathit{(}\mathit{\beta },\mathit{\rho },\mathit{\delta })$. Unfortunately, $\mathcal{J}(\mathit{\beta },\mathit{\rho },\mathit{\delta })$ tends to have multiple stationary points, and the limit point of the ECME algorithm may be an undesirable stationary point. To deal with this issue, we need to provide an accurate initial point. Following the method in [19], we can assume that the noise is uniform and then evaluate $\mathcal{G}\left(\mathit{\beta }\right)$ in (6) on a coarse V-dimensional grid to find a grid point, close to the global minimum of $\mathcal{G}\left(\mathit{\beta }\right)$, as ${\mathit{\beta }}^{\left(0\right)}$ of the ECME algorithm. We can also use the estimate of $\mathit{\beta }$, obtained by a subspace [29] or a sparse representation based [30] algorithm, as ${\mathit{\beta }}^{\left(0\right)}$ due to the higher accuracy of the stochastic ML estimate of $\mathit{\beta }$ [2].

On the boundary of the positive semi-definite region of $\mathit{\rho }$, i.e., the set $ð=\{\rho \mid \mathbf{O}⪰{\mathbf{0}}_V\mathrm{and}\mathcal{R}\left[\mathbf{O}\right]<V\}$, we give the following proposition. Let $\mathcal{N}\left[\mathbf{O}\right]$ denote the null space of $\mathbf{O}$.

Proposition 2.

In the ECME algorithm, $\mathcal{N}\left[{\mathbf{O}}^{\left(d\right)}\right]=\mathcal{N}\left[{\mathbf{O}}^{\left(0\right)}\right]$ for $d\ge 1$ if ${\mathbf{Q}}^{\left(0\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(0\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(0\right)}⪰{\mathbf{0}}_V$.

Proof. 

From Proposition 1, we first know that ${\mathbf{Q}}^{\left(d\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(d\right)}>\mathbf{0}$), ${\mathbf{O}}^{\left(d\right)}⪰{\mathbf{0}}_V$, and ${\mathbf{G}}^{\left(d\right)}\succ {\mathbf{0}}_W$ for $d\ge 0$ due to ${\mathbf{Q}}^{\left(0\right)}\succ {\mathbf{0}}_W$ (or ${\mathit{\delta }}^{\left(0\right)}>\mathbf{0}$) and ${\mathbf{O}}^{\left(0\right)}⪰{\mathbf{0}}_V$. Then, a proof by the mathematical induction method is given.

If ${\mathbf{O}}^{\left(u\right)}\mathbf{v}=\mathbf{0}$, we have ${\mathbf{O}}^{(u+1)}\mathbf{v}={\widehat{\mathbf{N}}}_k^{(u+1)}\mathbf{v}=\mathbf{0}$ in (13) and thus $\mathcal{N}\left[{\mathbf{O}}^{\left(u\right)}\right]\subseteq \mathcal{N}\left[{\mathbf{O}}^{(u+1)}\right]$. Furthermore, if ${\mathbf{O}}^{(u+1)}\mathbf{v}={\widehat{\mathbf{N}}}_k^{(u+1)}\mathbf{v}=\mathbf{0}$, we have ${\mathbf{v}}^H{\widehat{\mathbf{N}}}_k^{(u+1)}\mathbf{v}=0$ and in (13)

$$\begin{array}{ccc}\hfill {\mathbf{v}}^H{\left[{\mathbf{H}}^{\left(u\right)}\right]}^H\widehat{\mathbf{R}}{\mathbf{H}}^{\left(u\right)}\mathbf{v}=0& \Rightarrow & {\left[{\mathbf{H}}^{\left(u\right)}\right]}^H\widehat{\mathbf{R}}{\mathbf{H}}^{\left(u\right)}\mathbf{v}=\mathbf{0},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\mathbf{v}}^H\left({\mathbf{O}}^{\left(u\right)}-{\left[{\mathbf{H}}^{\left(u\right)}\right]}^H{\mathbf{G}}^{\left(u\right)}{\mathbf{H}}^{\left(u\right)}\right)\mathbf{v}=0& \Rightarrow & \left({\mathbf{O}}^{\left(u\right)}-{\left[{\mathbf{H}}^{\left(u\right)}\right]}^H{\mathbf{G}}^{\left(u\right)}{\mathbf{H}}^{\left(u\right)}\right)\mathbf{v}=\mathbf{0}.\hfill \end{array}$$

(23)

In order to proceed, we use the matrix inversion formula [22]

$$\begin{array}{c}\hfill {\mathbf{G}}^{-1}={\mathbf{Q}}^{-1/2}\left[{\mathbf{I}}_W-\tilde{\mathbf{A}}{(\mathbf{O}{\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}+{\mathbf{I}}_V)}^{-1}\mathbf{O}{\tilde{\mathbf{A}}}^H\right]{\mathbf{Q}}^{-1/2}\end{array}$$

(24)

and obtain

$$\begin{array}{c}\hfill \mathbf{O}-{\mathbf{H}}^H\mathbf{G}\mathbf{H}={(\mathbf{O}{\tilde{\mathbf{A}}}^H\tilde{\mathbf{A}}+{\mathbf{I}}_V)}^{-1}\mathbf{O},\end{array}$$

(25)

which suggests $\mathcal{N}\left[\mathbf{O}\right]=\mathcal{N}[\mathbf{O}-{\mathbf{H}}^H\mathbf{G}\mathbf{H}]$ and $\mathcal{N}\left[{\mathbf{O}}^{\left(u\right)}\right]=\mathcal{N}[{\mathbf{O}}^{\left(u\right)}-{\left[{\mathbf{H}}^{\left(u\right)}\right]}^H{\mathbf{G}}^{\left(u\right)}{\mathbf{H}}^{\left(u\right)}]$. Accordingly, ${\mathbf{O}}^{\left(u\right)}\mathbf{v}=\mathbf{0}$ in (22) and (23), leading to $\mathcal{N}\left[{\mathbf{O}}^{\left(u\right)}\right]\supseteq \mathcal{N}\left[{\mathbf{O}}^{(u+1)}\right]$. Finally, by combining $\mathcal{N}\left[{\mathbf{O}}^{\left(u\right)}\right]\subseteq \mathcal{N}\left[{\mathbf{O}}^{(u+1)}\right]$ and $\mathcal{N}\left[{\mathbf{O}}^{\left(u\right)}\right]\supseteq \mathcal{N}\left[{\mathbf{O}}^{(u+1)}\right]$, we obtain $\mathcal{N}\left[{\mathbf{O}}^{\left(u\right)}\right]=\mathcal{N}\left[{\mathbf{O}}^{(u+1)}\right]$. The proof is completed. □

Proposition 2 indicates that if ${\mathit{\rho }}^{\left(0\right)}$ in the ECME algorithm is on the boundary, i.e., ${\mathit{\rho }}^{\left(0\right)}\in ð$ and $\mathcal{N}\left[{\mathbf{O}}^{\left(0\right)}\right]$ is nonempty, the limit point of $\mathit{\rho }$ is also on the boundary. Hence, let $({\mathit{\beta }}^{*},{\mathit{\rho }}^{*},{\mathit{\delta }}^{*})$ denote the solution of Problem (3) and if ${\mathit{\rho }}^{*}\in ð$, we may need to estimate $\mathcal{N}\left[{\mathbf{O}}^{*}\right]$ before implementing the ECME algorithm. Fortunately, $({\mathit{\beta }}^{*},{\mathit{\rho }}^{*},{\mathit{\delta }}^{*})$ is always an interior point of the parameter space (i.e., ${\mathit{\rho }}^{*}\notin ð$, ${\mathbf{O}}^{*}\succ {\mathbf{0}}_V$, and $\mathcal{N}\left[{\mathbf{O}}^{*}\right]$ is empty) in practice even if the true value of $\mathit{\rho }$ is on the boundary. As a result, we can always adopt ${\mathbf{O}}^{\left(0\right)}\succ {\mathbf{0}}_V$ in the ECME algorithm, e.g., the simulation results in Figure 1 related to coherent sources.

4. Simulation Results

Simulation results are provided to confirm the effectiveness of the ECME algorithm, i.e., the ECME algorithm is able to obtain the ML estimate of $\mathit{\beta }$ in (3). We set $V=2$, ${\beta }_1={50}^{\circ }$, $W=6$, ${\beta }_2={100}^{\circ }$, and $\mathit{\delta }={\left[1234210\right]}^T$. Algorithm 1 is used to obtain ${\mathit{\beta }}^{\left(d\right)}$ in (20) and $\parallel {\mathit{\beta }}^{(u+1)}-{\mathit{\beta }}^{\left(u\right)}{\parallel }_2\le 0.{001}^{\circ }$ is adopted as the stopping criterion. The ECME algorithm is given an accurate initial point for obtaining the ML estimate of $\mathit{\beta }$. In Figure 1 and Figure 2, we consider the coherent (or fully correlated) source model with [25]

$$\begin{array}{c}\hfill \mathbf{O}=\left[\begin{array}{cc}\gamma & \gamma \\ \gamma & \gamma \end{array}\right]=\left[\begin{array}{cc}\sqrt{\gamma }& 0\\ 0& \sqrt{\gamma }\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]\left[\begin{array}{cc}\sqrt{\gamma }& 0\\ 0& \sqrt{\gamma }\end{array}\right].\end{array}$$

In Figure 3, we consider the partly correlate source model with

$$\begin{array}{c}\hfill \mathbf{O}=\left[\begin{array}{cc}5& 4\\ 4& 5\end{array}\right]=\left[\begin{array}{cc}\sqrt{5}& 0\\ 0& \sqrt{5}\end{array}\right]\left[\begin{array}{cc}1& 4/5\\ 4/5& 1\end{array}\right]\left[\begin{array}{cc}\sqrt{5}& 0\\ 0& \sqrt{5}\end{array}\right].\end{array}$$

In Figure 1 and Figure 2, we compare the root mean square error (RMSE) performance of the ECME algorithm with the CRLB [4,5]. In addition, we simulate the second space-alternating generalized EM (SAGE) algorithm for uncorrelated sources in [20], and this SAGE algorithm adopts the same simulation settings in [20]. Each RMSE is based on 2000 independent trials and the two algorithms share the same initial point. As expected, the ECME algorithm obtains smaller RMSEs than the SAGE algorithm. More importantly, the ECME algorithm attains the CRLB of $\mathit{\beta }$ when the number of snapshots L or $\gamma$ is large, which coincides with the well-known conclusion that the stochastic CRLB of β can be achieved asymptotically by the stochastic ML estimator of β [2]. Hence, the ECME algorithm is able to obtain the stochastic ML estimate of $\mathit{\beta }$ in (3) given an accurate initial point.

In Figure 3, we compare the ECME algorithm with two subspace-based algorithms, which utilize the state-of-the-art subspace separation approaches in [17,18] and are called “Approach 1+Root-MUSIC” and “Approach 2+Root-MUSIC”, respectively. The three algorithms process the same snapshots of each trial. As expected, the ECME algorithm yields more closely spaced estimates of $({\beta }_1,{\beta }_2)$ centered on $({50}^{\circ },{100}^{\circ })$ since in DOA estimation, the ML technique offers the highest advantage in terms of accuracy.

5. Conclusions

In this paper, we employed and designed the ECME algorithm for stochastic ML direction finding, where sources may be correlated, in unknown nonuniform noise. Theoretical analysis indicated that the ECME algorithm is computationally efficient and operationally stable. Simulation results confirmed that the ECME algorithm can efficiently obtain the ML based DOA estimate of each stochastic source.