EM and SAGE Algorithms for DOA Estimation in the Presence of Unknown Uniform Noise

Abstract

The existing expectation maximization (EM) and space-alternating generalized EM (SAGE) algorithms are only applied to direction of arrival (DOA) estimation in known noise. In this paper, the two algorithms are designed for DOA estimation in unknown uniform noise. Both the deterministic and random signal models are considered. In addition, a new modified EM (MEM) algorithm applicable to the noise assumption is also proposed. Next, these EM-type algorithms are improved to ensure the stability when the powers of sources are not equal. After being improved, simulation results illustrate that the EM algorithm has similar convergence with the MEM algorithm, the SAGE algorithm outperforms the EM and MEM algorithms for the deterministic signal model, and the SAGE algorithm cannot always outperform the EM and MEM algorithms for the random signal model. Furthermore, simulation results show that processing the same snapshots from the random signal model, the SAGE algorithm for the deterministic signal model can require the fewest computations.

1. Introduction

Direction of arrival (DOA) estimation is an important part of array signal processing and some high-resolution estimation techniques have been developed in the literature [1,2]. In particular, the maximum likelihood (ML) technique plays a critical role since it can offer the highest advantage in terms of both accuracy and spatial resolution. However, ML direction finding problems are non-convex and difficult to obtain their solutions in closed form.

One computationally efficient method to solve ML estimation problems is the classic expectation maximization (EM) algorithm [3,4], which has been employed for ML direction finding [5,6]. Each iteration of the EM algorithm is composed of an expectation step (E-step) and a maximization step (M-step). At the M-step, however, the EM algorithm updates all of the parameter estimates simultaneously, which causes slow convergence. In order to speed up the convergence of the EM algorithm, the space-alternating generalized EM (SAGE) algorithm has been proposed in [7]. References [8,9] show that the SAGE algorithm does yield faster convergence in terms of DOA estimation.

The existing EM and SAGE algorithms are usually derived under known noise [5,6,8,9]. The known noise is without unknown parameters, which may be unrealistic in certain applications. In fact, many seminal works in ML direction finding consider the so-called unknown uniform noise model [1,2,10,11,12], i.e., the noise covariance matrix can be expressed as $\tau {\mathbf{I}}_K$, where $\tau$ is the only unknown noise parameter and ${\mathbf{I}}_K$ is the $K\times K$ identity matrix. Under this noise assumption, a computationally attractive alternating projection algorithm is presented for computing the deterministic ML based DOA estimator in [10]. The authors in [11] investigate the statistical performance of this ML estimator and derive the Cramer–Rao lower bound. Moreover, some statistical properties of both the deterministic and random ML based DOA estimators under unknown uniform noise are compared in [12]. In addition to uniform noise, non-uniform noise has also attracted increasing attention [13,14,15,16]. The non-uniform noise has an arbitrary diagonal covariance matrix and, thus, makes DOA estimation complex. For efficiently computing both the deterministic and random ML based DOA estimators in unknown non-uniform noise, some feasible algorithms have been proposed in [17,18,19,20].

In this paper, we apply and design the EM and SAGE algorithms for DOA estimation in unknown uniform noise. Theoretical analyses indicate that for the deterministic signal model, $\tau$ has little effect on the two algorithms. However, the problem in the M-step of the EM algorithm for the random signal model can be no longer decomposed into parallel subproblems easily when $\tau$ is unknown. To proceed, we divide the M-step into two conditional M-steps (CM-steps) based on the expectation CM (ECM) algorithm [21]. In addition, we propose a new modified EM (MEM) algorithm applicable to the unknown uniform noise assumption. Note that although the EM algorithm in [22] is similar to the MEM algorithm, it is incorrectly derived. In brief, the proposed EM-type algorithms only need low-dimensional numerical searches at each iteration and are easy to perform. However, the proposed EM-type algorithms require accurate initial points, which is generally a computationally expensive task.

Existing simulations using EM-type algorithms always adopt sources of equal power [5,6,8,9,22]. We find, however, that when the powers of sources are unequal, the DOA estimates of multiple sources obtained by these EM-type algorithms tend to be consistent with the true DOA of the source with the largest power. To this end, we improve the proposed EM-type algorithms. After being improved, simulation results illustrate that (1) the EM algorithm has similar convergence with the MEM algorithm, (2) the SAGE algorithm outperforms the EM and MEM algorithms for the deterministic signal model, i.e., the SAGE algorithm converges faster and can avoid the convergence to an undesirable stationary point of the log-likelihood function (LLF) more efficiently, and (3) the SAGE algorithm cannot always outperform the EM and MEM algorithms for the random signal model.

The proposed EM-type algorithms for the deterministic signal model can process snapshots from the random signal model, so we, via simulation, compare these algorithms for both signal models. Simulations show that, under the same snapshots, initial DOA estimates, and stopping criterion, the SAGE algorithm for the deterministic signal model can require the fewest iterations and computations.

The contributions of this paper can be enumerated as follows:

• We apply and design the EM and SAGE algorithms for DOA estimation in unknown uniform noise. In particular, we derive the SAGE algorithm for random ML direction finding, which is not discussed in [7,8,22].
• We propose a new MEM algorithm applicable to the unknown uniform noise assumption.
• We improve these EM-type algorithms to ensure the stability when the powers of sources are not equal.
• Via simulation we show that the EM algorithm has similar convergence with the MEM algorithm and the SAGE algorithm outperforms the EM and MEM algorithms for the deterministic signal model. However, the SAGE algorithm cannot always outperform the EM and MEM algorithms for the random signal model.
• Via simulation we show that processing the same snapshots from the random signal model, the SAGE algorithm for the deterministic signal model can require the fewest iterations and computations.

The rest of this paper is outlined as follows: in Section 2, we formulate both the deterministic and random ML direction finding problems in unknown uniform noise. In Section 3, Section 4 and Section 5, we design the EM, MEM, and SAGE algorithms, respectively. We analyze some convergence properties of these EM-type algorithms in Section 6 and provide simulation results to compare the convergence of these EM-type algorithms in Section 7. Finally, we conclude this paper in Section 8.

2. Signal Model and Problem Statement

An array of K sensors is assumed to receive the plane waves emitted from G narrowband sources, which share the same known center wavelength $\chi$. We use the Cartesian coordinate ${\zeta }_k={\left[x_k{y}_k{z}_k\right]}^T$ and the Spherical coordinate $(1,{\mu }_g,{\eta }_g)$ to locate the kth sensor and the direction of the gth source, respectively. Here, ${[·]}^T$ denotes transpose, ${\mu }_g$ and ${\eta }_g$ denote the elevation and azimuth angles of the gth source, respectively. For convenience, we transform $(1,{\mu }_g,{\eta }_g)$ into the corresponding Cartesian coordinate ${\gamma }_g={[sin\left({\mu }_g\right)cos\left({\eta }_g\right)sin\left({\mu }_g\right)sin\left({\eta }_g\right)cos\left({\mu }_g\right)]}^T$. Let the origin be the reference point such that the signal received at the array is written as

$$\begin{array}{c}\hfill \mathbf{w}\left(t\right)=\sum _{g=1}^G\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)+\mathbf{v}\left(t\right)=\mathbf{B}\left(\xi \right)\mathbf{m}\left(t\right)+\mathbf{v}\left(t\right),\end{array}$$

(1)

where $\mathbf{b}\left({\xi }_g\right)={[e^{𝚥{\varphi }_{1,g}}\cdots e^{𝚥{\varphi }_{K,g}}]}^T$ with DOA ${\xi }_g=({\mu }_g,{\eta }_g)\in \Omega$, ${\varphi }_{k,g}=-\frac{2\pi }{\chi }{\zeta }_k^T{\gamma }_g$, and $𝚥=\sqrt{-1}$, $m_g\left(t\right)$ is the signal of the gth source, and $\mathbf{v}\left(t\right)$ is complex Gaussian noise of zero mean and covariance $\tau {\mathbf{I}}_K$, i.e., $\mathbf{v}\left(t\right)\sim \mathcal{CN}(\mathbf{0},\tau {\mathbf{I}}_K)$ with $\mathbf{0}={[0\cdots 0]}^T$. In (1), $\mathbf{B}(\mathit{\xi })=[\mathbf{b}\left({\xi }_1\right)\cdots \mathbf{b}\left({\xi }_G\right)]$, $\mathit{\xi }=({\xi }_1,\cdots ,{\xi }_G)\in \Omega$ with $\Omega ={\Omega }^G$, and $\mathbf{m}\left(t\right)={[m_1\left(t\right)\cdots m_G\left(t\right)]}^T$.

EM-type algorithms need to define the unavailable complete data. According to the classic EM paradigm for superimposed signals [5,6], we construct L independent snapshots, the incomplete data of the EM algorithm, by 

$$\begin{array}{c}\hfill \mathbf{w}\left(t\right)=\sum _{g=1}^G\left[\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)+{\mathbf{v}}_g\left(t\right)\right]=\sum _{g=1}^G{\mathbf{h}}_g\left(t\right),t=1,2,\cdots ,L,\end{array}$$

(2)

where the ${\mathbf{h}}_g\left(t\right)$’s are the complete data. Moreover, the ${\mathbf{v}}_g\left(t\right)$’s are mutually uncorrelated and ${\mathbf{v}}_g\left(t\right)\sim \mathcal{CN}(\mathbf{0},{\beta }_g\tau {\mathbf{I}}_K)$, where $\beta ={[{\beta }_1\cdots {\beta }_G]}^T>\mathbf{0}$ and ${\mathbf{1}}^T\beta =1$ with $\mathbf{1}={[1\cdots 1]}^T$. Note that the incomplete- and complete-data LLFs require the distributions of the $m_g\left(t\right)$’s, we adopt the following two statistical models separately.

2.1. Deterministic Signal Model

We let the $m_g\left(t\right)$’s be deterministic and unknown [5,6,10,11], which leads to ${\mathbf{h}}_g\left(t\right)\sim \mathcal{CN}\left(\mathbf{b}\left({\xi }_g\right)m_g\left(t\right),{\beta }_g\tau {\mathbf{I}}_K\right)$ and $\mathbf{w}\left(t\right)\sim \mathcal{CN}\left(\mathbf{B}\left(\xi \right)\mathbf{m}\left(t\right),\tau {\mathbf{I}}_K\right)$. Then, the incomplete- and complete-data LLFs are formulated by

$$\begin{array}{ccc}\hfill \mathcal{J}(\mathbf{\Psi },\tau )& =& \sum _{t=1}^{L}logp\left(\mathbf{w}\left(t\right);\xi ,\mathbf{m}\left(t\right),\tau \right)\hfill \\ & =& -LKlog\left(\pi \tau \right)-\frac{1}{\tau }\sum _{t=1}^L{\parallel \mathbf{w}\left(t\right)-\mathbf{B}\left(\xi \right)\mathbf{m}\left(t\right)\parallel }^2,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \mathcal{Z}(\mathbf{\Psi },\tau )& =& \sum _{t=1}^L\sum _{g=1}^{G}logp\left({\mathbf{h}}_g\left(t\right);{\xi }_g,m_g\left(t\right),\tau \right)\hfill \\ & =& -LKGlog\left(\pi \tau \right)-LK\sum _{g=1}^{G}log\left({\beta }_g\right)-\frac{1}{\tau }\sum _{t=1}^L\sum _{g=1}^G\frac{1}{{\beta }_g}{\parallel {\mathbf{h}}_g\left(t\right)-\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)\parallel }^2,\hfill \end{array}$$

(4)

where $\parallel ·\parallel$ denotes Euclidean norm and $\mathbf{M}=\left[\mathbf{m}\right(1)\cdots \mathbf{m}(L\left)\right]$. Note that $\mathbf{\Psi }=(\mathit{\xi },\mathbf{M})$ denotes the signal parameters while $\tau$ is the only noise parameter. Finally, the ML estimation problem is

$$\begin{array}{c}\hfill \underset{\xi \in \Omega ,\mathbf{M},\tau >0}{max}\mathcal{J}(\mathbf{\Psi },\tau ).\end{array}$$

(5)

2.2. Random Signal Model

We assume $m_g\left(t\right)\sim \mathcal{CN}(0,{\rho }_g)$ where ${\rho }_g$ is the power of the gth source. For simplicity, all of the $m_g\left(t\right)$’s and ${\mathbf{v}}_g\left(t\right)$’s are assumed to be mutually uncorrelated [5,6]. Next, we have ${\mathbf{h}}_g\left(t\right)\sim \mathcal{CN}(\mathbf{0},{\mathbf{N}}_g)$ with ${\mathbf{N}}_g={\rho }_g\mathbf{b}\left({\xi }_g\right){\mathbf{b}}^H\left({\xi }_g\right)+{\beta }_g\tau {\mathbf{I}}_K$, where ${[·]}^H$ denotes conjugate transpose, and $\mathbf{w}\left(t\right)\sim \mathcal{CN}(\mathbf{0},{\mathbf{N}}_w)$ with ${\mathbf{N}}_w={\sum }_{g=1}^G{\mathbf{N}}_g$. The incomplete- and complete-data LLFs are formulated by

$$\begin{array}{ccc}\hfill \mathcal{J}(\mathbf{\Psi },\tau )& =& \sum _{t=1}^{L}logp\left(\mathbf{w}\left(t\right);\xi ,\rho ,\tau \right)\hfill \\ & =& -L\left[Klog\left(\pi \right)+log\left(\mathrm{Det}\left({\mathbf{N}}_w\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_w^{-1}{\widehat{\mathbf{P}}}_w\right)\right],\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \mathcal{Z}(\mathbf{\Psi },\tau )& =& {\displaystyle \sum _{t=1}^L}{\displaystyle \sum _{g=1}^G}logp\left({\mathbf{h}}_g\left(t\right);{\xi }_g,{\rho }_g,\tau \right)\hfill \\ & =& -L{\displaystyle \sum _{g=1}^G}\left[Klog\left(\pi \right)+log\left(\mathrm{Det}\left({\mathbf{N}}_g\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_g^{-1}{\widehat{\mathbf{P}}}_g\right)\right],\hfill \end{array}$$

(7)

where ${[·]}^{-1}$, $\mathrm{Det}(·)$, and $\mathrm{Tr}(·)$ denote inversion, determinant, and trace, respectively. Moreover, $\mathbf{\Psi }=(\mathit{\xi },\mathit{\rho })$, ${\widehat{\mathbf{P}}}_w=(1/L){\sum }_{t=1}^L\mathbf{w}\left(t\right){\mathbf{w}}^H\left(t\right)$, ${\widehat{\mathbf{P}}}_g=(1/L){\sum }_{t=1}^L{\mathbf{h}}_g\left(t\right){\mathbf{h}}_g^H\left(t\right)$, and $\mathit{\rho }={[{\rho }_1\cdots {\rho }_G]}^T$. Finally, the ML estimation problem is

$$\begin{array}{c}\hfill \underset{\xi \in \Omega ,\mathit{\rho }\ge \mathbf{0},\tau >0}{max}\mathcal{J}(\mathbf{\Psi },\tau ).\end{array}$$

(8)

3. EM Algorithm

In this section, we design and derive the EM algorithm for solving problems (5) and (8). The E- and M-steps at the rth iteration are introduced below. Let ${[·]}^{\left(r\right)}$, $\mathcal{E}\{·\}$, and $\mathcal{D}\{·\}$ denote an iterative value at the rth iteration, expectation, and covariance, respectively. ${[·]}^{\left(0\right)}$ is an initial estimate.

3.1. Deterministic Signal Model

3.1.1. E-Step

Calculate the conditional expectation of the complete-data LLF in (4)

$$\begin{array}{ccc}\hfill \mathcal{E}\left\{\mathcal{Z}(\mathbf{\Delta })|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}& =& C-L\{KGlog\left(\tau \right)+K\sum _{g=1}^{G}log\left({\beta }_g\right)+\frac{1}{\tau }[u^{\left(r\right)}+\hfill \\ & & \frac{1}{L}\sum _{t=1}^L\sum _{g=1}^G\frac{1}{{\beta }_g}{\parallel {\mathbf{h}}_g^{\left(r\right)}\left(t\right)-\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)\parallel }^2\left]\right\},\hfill \end{array}$$

(9)

where $\mathbf{\Delta }=(\mathbf{\Psi },\mathbf{\tau })=(\mathbf{\xi },\mathbf{M},\mathbf{\tau })$, $\mathbf{W}=\left[\mathbf{w}\right(1)\mathbf{\cdots }\mathbf{w}(L\left)\right]$, and $C=-LKGlog\left(\pi \right)$. In (9), the conditional distribution of ${\mathbf{h}}_g\left(t\right)$ can be derived from [23] and

$$\begin{array}{ccc}\hfill {\mathbf{h}}_g^{\left(r\right)}\left(t\right)& =& \mathcal{E}\left\{{\mathbf{h}}_g\left(t\right)|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}\hfill \\ & =& \mathbf{b}\left({\xi }_g^{(r-1)}\right)m_g^{(r-1)}\left(t\right)+{\beta }_g\left[\mathbf{w}\left(t\right)-\mathbf{B}\left({\xi }^{(r-1)}\right){\mathbf{m}}^{(r-1)}\left(t\right)\right],\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill u^{\left(r\right)}& =& \frac{1}{L}\sum _{t=1}^L\sum _{g=1}^G\frac{1}{{\beta }_g}\mathrm{Tr}\left(\mathcal{D}\left\{{\mathbf{h}}_g\left(t\right)|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}\right)=\frac{1}{L}\sum _{t=1}^L\sum _{g=1}^G\frac{1}{{\beta }_g}\mathrm{Tr}\left({\beta }_g(1-{\beta }_g){\tau }^{(r-1)}{\mathbf{I}}_K\right)\hfill \\ & =& K(G-1){\tau }^{(r-1)}.\hfill \end{array}$$

(11)

3.1.2. M-Step

Update the estimates of $\mathbf{\Psi }$ and $\tau$ by solving

$$\begin{array}{c}\hfill \underset{\xi \in \Omega ,\mathbf{M},\tau >0}{min}KGlog\left(\tau \right)+\frac{1}{\tau }\left[u^{\left(r\right)}+\frac{1}{L}\sum _{t=1}^L\sum _{g=1}^G\frac{1}{{\beta }_g}\parallel {\mathbf{h}}_g^{\left(r\right)}\left(t\right)-\mathbf{b}\left({\xi }_g\right)m_g\left(t\right){\parallel }^2\right].\end{array}$$

(12)

${\mathbf{\Psi }}^{\left(r\right)}=({\mathit{\xi }}^{\left(r\right)},{\mathbf{M}}^{\left(r\right)})$ and ${\tau }^{\left(r\right)}$ are obtained by [6]

$$\begin{array}{ccc}\hfill {\xi }_g^{\left(r\right)}& =& arg\underset{{\xi }_g\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill m_g^{\left(r\right)}\left(t\right)& =& {\mathbf{b}}^H\left({\xi }_g^{\left(r\right)}\right){\mathbf{h}}_g^{\left(r\right)}\left(t\right)/K,\forall g,t,\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill {\tau }^{\left(r\right)}& =& (1-1/G){\tau }^{(r-1)}+(1/K/G){\sum }_{g=1}^G{d}_g^{\left(r\right)}/{\beta }_g,\hfill \end{array}$$

(15)

where ${\mathbf{\Gamma }}_g=\mathbf{b}\left({\xi }_g\right){\mathbf{b}}^H\left({\xi }_g\right)/K$, ${\widehat{\mathbf{P}}}_g^{\left(r\right)}=(1/L){\sum }_{t=1}^L{\mathbf{h}}_g^{\left(r\right)}\left(t\right){\left[{\mathbf{h}}_g^{\left(r\right)}\left(t\right)\right]}^H$, and $d_g^{\left(r\right)}=\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_g^{\left(r\right)}){\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\ge 0$. In (15), ${\tau }^{\left(r\right)}>0$ if ${\tau }^{(r-1)}>0$.

Remark 1.

Note that the ${\mathbf{h}}_g^{\left(r\right)}\left(t\right)$’s in (10), the ${\xi }_g^{\left(r\right)}$’s in (13), and the $m_g^{\left(r\right)}\left(t\right)$’s in (14) are unrelated to ${\tau }^{(r-1)}$, we can omit (15) due to the nuisance parameter τ.

3.2. Random Signal Model

3.2.1. E-Step

Calculate the conditional expectation of the complete-data LLF in (7)

$$\begin{array}{c}\hfill \mathcal{E}\left\{\mathcal{Z}(\mathbf{\Delta })|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}=C-L\sum _{g=1}^G\left[log\left(\mathrm{Det}\left({\mathbf{N}}_g\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_g^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right],\end{array}$$

(16)

where $\mathbf{\Delta }=(\mathbf{\Psi },\tau )=(\mathit{\xi },\mathit{\rho },\tau )$ and

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{P}}}_g^{\left(r\right)}& =& \mathcal{E}\left\{{\widehat{\mathbf{P}}}_g|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}\hfill \\ & =& {\mathbf{N}}_g^{(r-1)}-{\mathbf{N}}_g^{(r-1)}{\left[{\mathbf{N}}_w^{(r-1)}\right]}^{-1}{\mathbf{N}}_g^{(r-1)}+{\mathbf{N}}_g^{(r-1)}{\left[{\mathbf{N}}_w^{(r-1)}\right]}^{-1}{\widehat{\mathbf{P}}}_w{\left[{\mathbf{N}}_w^{(r-1)}\right]}^{-1}{\mathbf{N}}_g^{(r-1)}.\hfill \end{array}$$

(17)

3.2.2. M-Step

Update the estimates of $\mathbf{\Psi }$ and $\tau$ by solving

$$\begin{array}{c}\hfill \underset{\xi \in \Omega ,\rho \ge \mathbf{0},\tau >0}{min}\sum _{g=1}^G\left[log\left(\mathrm{Det}\left({\mathbf{N}}_g\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_g^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right],\end{array}$$

(18)

which is difficult to be decomposed into parallel subproblems due to $\tau$. To obtain ${\mathbf{\Psi }}^{\left(r\right)}$ and ${\tau }^{\left(r\right)}$ easily, we rewrite (18) as

$$\begin{array}{c}\hfill \underset{\xi \in \Omega ,\sigma \ge \mathbf{0},\tau >0}{min}KGlog\left(\tau \right)+\sum _{g=1}^G\left[log\left(\mathrm{Det}\left({\mathbf{Q}}_g\right)\right)+\frac{1}{\tau }\mathrm{Tr}\left({\mathbf{Q}}_g^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right],\end{array}$$

(19)

where ${\mathbf{N}}_g=\tau {\mathbf{Q}}_g$ with ${\mathbf{Q}}_g={\sigma }_g\mathbf{b}\left({\xi }_g\right){\mathbf{b}}^H\left({\xi }_g\right)+{\beta }_g{\mathbf{I}}_K$ and $\sigma ={[{\sigma }_1\cdots {\sigma }_G]}^T$ with ${\sigma }_g={\rho }_g/\tau$. We now divide the M-step into the following two CM-steps based on the ECM algorithm [21], i.e., the algorithm becomes the ECM algorithm. For convenience, we still call it the EM algorithm.

• First CM-step: Estimate $\mathbf{\Psi }$ but hold $\tau ={\tau }^{(r-1)}$ fixed. Then, problem (19) can be decomposed into the G parallel subproblems

  $$\begin{array}{c}\hfill \underset{{\xi }_g\in \Omega ,{\sigma }_g\ge 0}{min}log\left(\mathrm{Det}\left({\mathbf{Q}}_g\right)\right)+\frac{1}{{\tau }^{(r-1)}}\mathrm{Tr}\left({\mathbf{Q}}_g^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g.\end{array}$$

  (20)
  ${\mathbf{\Psi }}^{\left(r\right)}=({\mathit{\xi }}^{\left(r\right)},{\mathit{\rho }}^{\left(r\right)})$ is obtained by [6]

  $$\begin{array}{ccc}\hfill {\xi }_g^{\left(r\right)}& =& arg\underset{{\xi }_g\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g,\hfill \end{array}$$

  (21)

  $$\begin{array}{ccc}\hfill {\rho }_g^{\left(r\right)}& =& {\sigma }_g^{\left(r\right)}{\tau }^{(r-1)}=max\left\{\left(e_g^{\left(r\right)}-{\beta }_g{\tau }^{(r-1)}\right)/K,0\right\},\forall g,\hfill \end{array}$$

  (22)

  where $e_g^{\left(r\right)}=\mathrm{Tr}\left({\mathbf{\Gamma }}_g^{\left(r\right)}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)$ and ${\xi }_g^{\left(r\right)}$ is indeterminate if ${\rho }_g^{\left(r\right)}=0$.
• Second CM-step: Estimate $\tau$ but hold $\mathbf{\Psi }={\mathbf{\Psi }}^{\left(r\right)}$ fixed. Then, problem (19) is simplified to

  $$\begin{array}{c}\hfill \underset{\tau >0}{min}KGlog\left(\tau \right)+\frac{1}{\tau }\sum _{g=1}^G\mathrm{Tr}\left({\left[{\mathbf{Q}}_g^{\left(r\right)}\right]}^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right).\end{array}$$

  (23)
  ${\tau }^{\left(r\right)}$ is obtained by

  $$\begin{array}{c}\hfill {\tau }^{\left(r\right)}=\frac{1}{KG}\sum _{g=1}^G\mathrm{Tr}\left({\left[{\mathbf{Q}}_g^{\left(r\right)}\right]}^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)=\frac{1}{K}{\tau }^{(r-1)}+\frac{1}{KG}\sum _{g=1}^G{d}_g^{\left(r\right)}/{\beta }_g,\end{array}$$

  (24)

  where $d_g^{\left(r\right)}=\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_g^{\left(r\right)}){\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\ge 0$ and ${\tau }^{\left(r\right)}>0$ if ${\tau }^{(r-1)}>0$.

4. MEM Algorithm

In the previous section, $\beta$ is fixed and known. In this section, we regard $\beta$ as a parameter to be estimated and, thus, propose an MEM algorithm applicable to the unknown uniform noise assumption.

To estimate $\tau$ and $\beta$ easily, we introduce ${\tau }_g={\beta }_g\tau$ as the common noise variance of the gth source and have

$$\begin{array}{c}\hfill {\mathbf{v}}_g\left(t\right)\sim \mathcal{CN}(\mathbf{0},{\tau }_g{\mathbf{I}}_K).\end{array}$$

(25)

Clearly, $\tau ={\sum }_{g=1}^G{\tau }_g$ and ${\beta }_g={\tau }_g/\tau$. The E- and M-steps at the rth iteration are introduced below. Let $\mathit{\tau }={[{\tau }_1\cdots {\tau }_G]}^T$.

4.1. Deterministic Signal Model

Based on (25), the complete-data LLF in (4) is rewritten as

$$\begin{array}{c}\hfill \mathcal{Z}(\mathbf{\Delta })=C-LK\sum _{g=1}^{G}log\left({\tau }_g\right)-\sum _{t=1}^L\sum _{g=1}^G\frac{1}{{\tau }_g}{\parallel {\mathbf{h}}_g\left(t\right)-\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)\parallel }^2,\end{array}$$

(26)

where $\mathbf{\Delta }=(\mathbf{\Psi },\mathit{\tau })=(\mathit{\xi },\mathbf{M},\mathit{\tau })$ and ${\mathbf{h}}_g\left(t\right)\sim \mathcal{CN}\left(\mathbf{b}\left({\xi }_g\right)m_g\left(t\right),{\tau }_g{\mathbf{I}}_K\right)$.

4.1.1. E-Step

Calculate the conditional expectation of the complete-data LLF in (26)

$$\begin{array}{ccc}\hfill \mathcal{E}\left\{\mathcal{Z}(\mathbf{\Delta })|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}& =& C-L\sum _{g=1}^G\{Klog\left({\tau }_g\right)+\frac{1}{{\tau }_g}[u_g^{\left(r\right)}+\hfill \\ & & \frac{1}{L}\sum _{t=1}^L{\parallel {\mathbf{h}}_g^{\left(r\right)}\left(t\right)-\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)\parallel }^2]\},\hfill \end{array}$$

(27)

where

$$\begin{array}{ccc}\hfill {\mathbf{h}}_g^{\left(r\right)}\left(t\right)& =& \mathcal{E}\left\{{\mathbf{h}}_g\left(t\right)|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}\hfill \\ & =& \mathbf{b}\left({\xi }_g^{(r-1)}\right)m_g^{(r-1)}\left(t\right)+\left({\tau }_g^{(r-1)}/{\tau }^{(r-1)}\right)\left[\mathbf{w}\left(t\right)-\mathbf{B}\left({\xi }^{(r-1)}\right){\mathbf{m}}^{(r-1)}\left(t\right)\right],\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill u_g^{\left(r\right)}& =& \frac{1}{L}\sum _{t=1}^L\mathrm{Tr}\left(\mathcal{D}\left\{{\mathbf{h}}_g\left(t\right)|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}\right)=K{\tau }_g^{(r-1)}\left(1-{\tau }_g^{(r-1)}/{\tau }^{(r-1)}\right).\hfill \end{array}$$

(29)

4.1.2. M-Step

Update the estimates of $\mathbf{\Psi }$ and $\mathbf{\tau }$ by solving the G parallel subproblems

$$\begin{array}{c}\hfill \underset{{\xi }_g\in \Omega ,{\mathbf{m}}_g,{\tau }_g>0}{min}Klog\left({\tau }_g\right)+\frac{1}{{\tau }_g}\left[u_g^{\left(r\right)}+\frac{1}{L}\sum _{t=1}^L\parallel {\mathbf{h}}_g^{\left(r\right)}\left(t\right)-\mathbf{b}\left({\xi }_g\right)m_g\left(t\right){\parallel }^2\right],\forall g,\end{array}$$

(30)

where ${\mathbf{m}}_g=[m_g\left(1\right)\cdots m_g\left(L\right)]$. ${\mathbf{\Psi }}^{\left(r\right)}=({\mathbf{\xi }}^{\left(r\right)},{\mathbf{M}}^{\left(r\right)})$ and ${\mathbf{\tau }}^{\left(r\right)}$ are obtained by

$$\begin{array}{ccc}\hfill {\xi }_g^{\left(r\right)}& =& arg\underset{{\xi }_g\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g,\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill m_g^{\left(r\right)}\left(t\right)& =& {\mathbf{b}}^H\left({\xi }_g^{\left(r\right)}\right){\mathbf{h}}_g^{\left(r\right)}\left(t\right)/K,\forall g,t,\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill {\tau }_g^{\left(r\right)}& =& {\tau }_g^{(r-1)}\left(1-{\tau }_g^{(r-1)}/{\tau }^{(r-1)}\right)+d_g^{\left(r\right)}/K,\forall g,\hfill \end{array}$$

(33)

where ${\widehat{\mathbf{P}}}_g^{\left(r\right)}=(1/L){\sum }_{t=1}^L{\mathbf{h}}_g^{\left(r\right)}\left(t\right){\left[{\mathbf{h}}_g^{\left(r\right)}\left(t\right)\right]}^H$ and $d_g^{\left(r\right)}=\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_g^{\left(r\right)}){\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\ge 0$. In (33), if ${\mathit{\tau }}^{(r-1)}>\mathbf{0}$, we have $1-{\tau }_g^{(r-1)}/{\tau }^{(r-1)}>0,\forall g$, and then ${\tau }_g^{\left(r\right)}>0,\forall g$, i.e., ${\mathit{\tau }}^{\left(r\right)}>\mathbf{0}$.

4.2. Random Signal Model

Based on (25), the complete-data LLF in (7) is rewritten as

$$\begin{array}{c}\hfill \mathcal{Z}(\mathbf{\Delta })=C-L\sum _{g=1}^G\left[log\left(\mathrm{Det}\left({\mathbf{N}}_g\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_g^{-1}{\widehat{\mathbf{P}}}_g\right)\right],\end{array}$$

(34)

where $\mathbf{\Delta }=(\mathbf{\Psi },\mathit{\tau })=(\mathit{\xi },\mathit{\rho },\mathit{\tau })$ and ${\mathbf{N}}_g={\rho }_g\mathbf{b}\left({\xi }_g\right){\mathbf{b}}^H\left({\xi }_g\right)+{\tau }_g{\mathbf{I}}_K$.

4.2.1. E-Step

Calculate the conditional expectation of the complete-data LLF in (34)

$$\begin{array}{c}\hfill \mathcal{E}\left\{\mathcal{Z}(\mathbf{\Delta })|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}=C-L\sum _{g=1}^G\left[log\left(\mathrm{Det}\left({\mathbf{N}}_g\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_g^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right],\end{array}$$

(35)

where

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{P}}}_g^{\left(r\right)}& =& \mathcal{E}\left\{{\widehat{\mathbf{P}}}_g|\mathbf{W};{\mathbf{\Delta }}^{(r-1)}\right\}\hfill \\ & =& {\mathbf{N}}_g^{(r-1)}-{\mathbf{N}}_g^{(r-1)}{\left[{\mathbf{N}}_w^{(r-1)}\right]}^{-1}{\mathbf{N}}_g^{(r-1)}+{\mathbf{N}}_g^{(r-1)}{\left[{\mathbf{N}}_w^{(r-1)}\right]}^{-1}{\widehat{\mathbf{P}}}_w{\left[{\mathbf{N}}_w^{(r-1)}\right]}^{-1}{\mathbf{N}}_g^{(r-1)}.\hfill \end{array}$$

(36)

4.2.2. M-Step

Update the estimates of $\mathbf{\Psi }$ and $\mathit{\tau }$ by solving the G parallel subproblems

$$\begin{array}{c}\hfill \underset{{\xi }_g\in \Omega ,{\sigma }_g\ge 0,{\tau }_g>0}{min}Klog\left({\tau }_g\right)+log\left(\mathrm{Det}\left({\mathbf{Q}}_g\right)\right)+\frac{1}{{\tau }_g}\mathrm{Tr}\left({\mathbf{Q}}_g^{-1}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g,\end{array}$$

(37)

where ${\mathbf{N}}_g={\tau }_g{\mathbf{Q}}_g$ with ${\mathbf{Q}}_g={\sigma }_g\mathbf{b}\left({\xi }_g\right){\mathbf{b}}^H\left({\xi }_g\right)+{\mathbf{I}}_K$ and ${\sigma }_g={\rho }_g/{\tau }_g$. Since $\mathrm{Det}\left({\mathbf{Q}}_g\right)=K{\sigma }_g+1$ and ${\mathbf{Q}}_g^{-1}={\mathbf{I}}_K-\frac{K{\sigma }_g}{K{\sigma }_g+1}{\mathbf{\Gamma }}_g$, subproblems (37) are rewritten as

$$\begin{array}{c}\hfill \underset{{\xi }_g\in \Omega ,{\sigma }_g\ge 0,{\tau }_g>0}{min}Klog\left({\tau }_g\right)+log(K{\sigma }_g+1)+\frac{1}{{\tau }_g}\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)-\frac{K{\sigma }_g}{{\tau }_g(K{\sigma }_g+1)}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g.\end{array}$$

(38)

To proceed, we first eliminate $\mathit{\sigma }={[{\sigma }_1\cdots {\sigma }_G]}^T$ in (38) [24,25]. Thus, when obtaining ${\mathit{\xi }}^{\left(r\right)}$ and ${\mathit{\tau }}^{\left(r\right)}$, ${\mathit{\sigma }}^{\left(r\right)}$ and ${\mathit{\rho }}^{\left(r\right)}$ are obtained by [6]

$$\begin{array}{ccc}\hfill {\sigma }_g^{\left(r\right)}& =& max\left\{\left(e_g^{\left(r\right)}/{\tau }_g^{\left(r\right)}-1\right)/K,0\right\},\forall g,\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill {\rho }_g^{\left(r\right)}& =& {\sigma }_g^{\left(r\right)}{\tau }_g^{\left(r\right)}=max\left\{\left(e_g^{\left(r\right)}-{\tau }_g^{\left(r\right)}\right)/K,0\right\},\forall g,\hfill \end{array}$$

(40)

where $e_g^{\left(r\right)}=\mathrm{Tr}\left({\mathbf{\Gamma }}_g^{\left(r\right)}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)$. Note that if ${\sigma }_g^{\left(r\right)}=0$, ${\xi }_g^{\left(r\right)}$ is indeterminate and ${\tau }_g^{\left(r\right)}=\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K$ by (38). To obtain ${\mathit{\xi }}^{\left(r\right)}$ and ${\mathit{\tau }}^{\left(r\right)}$, we assume ${\mathit{\sigma }}^{\left(r\right)}>\mathbf{0}$. After eliminating $\sigma$, subproblems (38) are simplified to

$$\begin{array}{c}\hfill \underset{{\xi }_g\in \Omega ,{\tau }_g>0}{min}(K-1)log\left({\tau }_g\right)+log\left(\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right)+\frac{1}{{\tau }_g}\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_g){\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g.\end{array}$$

(41)

Next, we eliminate $\mathit{\tau }$ in (41). Thus, when obtaining ${\mathit{\xi }}^{\left(r\right)}$, ${\mathit{\tau }}^{\left(r\right)}$ is obtained by ${\tau }_g^{\left(r\right)}=d_g^{\left(r\right)}/(K-1),\forall g$, where $d_g^{\left(r\right)}=\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_g^{\left(r\right)}){\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)$. After eliminating $\mathit{\tau }$, subproblems (41) are simplified to

$$\begin{array}{c}\hfill \underset{{\xi }_g\in \Omega }{min}(K-1)log\left(\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)-\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right)+log\left(\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\right),\forall g,\end{array}$$

(42)

where ${\xi }_g\in {\delta }_g^{\left(r\right)}=\left\{{\xi }_g\in \Omega \mid \mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)>\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K\right\}$ due to the fact that when ${\xi }_g^{\left(r\right)}\in {\delta }_g^{\left(r\right)}$, $e_g^{\left(r\right)}=\mathrm{Tr}\left({\mathbf{\Gamma }}_g^{\left(r\right)}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)>\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K$ and

$$\begin{array}{cc}\hfill & {\sigma }_g^{\left(r\right)}=max\left\{\left(e_g^{\left(r\right)}/{\tau }_g^{\left(r\right)}-1\right)/K,0\right\}=\left(e_g^{\left(r\right)}-\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K\right)/d_g^{\left(r\right)}>0,\forall g.\end{array}$$

(43)

Since $(K-1)log\left(\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)-x\right)+log\left(x\right)$ is a monotonically decreasing function of x for $x\ge \mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K$, subproblems (42) are equivalent to

$$\begin{array}{c}\hfill \underset{{\xi }_g\in {\delta }_g^{\left(r\right)}}{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g.\end{array}$$

(44)

Based on the above analysis, ${\mathbf{\Psi }}^{\left(r\right)}=({\mathit{\xi }}^{\left(r\right)},{\mathit{\rho }}^{\left(r\right)})$ and ${\mathit{\tau }}^{\left(r\right)}$ are obtained by

$$\begin{array}{ccc}\hfill {\xi }_g^{\left(r\right)}& =& arg\underset{{\xi }_g\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g,\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill {\tau }_g^{\left(r\right)}& =& \left\{\begin{array}{cc}{d}_g^{\left(r\right)}/(K-1)\ge 0,\hfill & e_g^{\left(r\right)}>\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K,\hfill \\ \mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K\ge 0,\hfill & e_g^{\left(r\right)}\le \mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K,\hfill \end{array}\right.\forall g,\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill {\rho }_g^{\left(r\right)}& =& max\left\{\left(e_g^{\left(r\right)}-\mathrm{Tr}\left({\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)/K\right)/(K-1),0\right\},\forall g,\hfill \end{array}$$

(47)

where we can use a proof by contradiction to verify that ${\mathit{\tau }}^{\left(r\right)}>\mathbf{0}$ if ${\mathit{\tau }}^{(r-1)}>\mathbf{0}$.

5. SAGE Algorithm

In the SAGE algorithm, each iteration consists of G cycles and ${\xi }_q^{\left(r\right)}$ is obtained at the qth cycle of the rth iteration. Let ${[·]}^{(r,q)}$ mean an iterative value at the qth cycle of the rth iteration, ${[·]}^{(r-1)}={[·]}^{(r-1,G)}={[·]}^{(r,0)}$.

At the qth cycle of the rth iteration, the SAGE algorithm first constructs the complete data by [7,8]

$${\mathbf{h}}_g\left(t\right)=\left\{\begin{array}{cc}\mathbf{b}\left({\xi }_g\right)m_g\left(t\right)+\mathbf{v}\left(t\right)\hfill & g=q,\hfill \\ \mathbf{b}\left({\xi }_g\right)m_g\left(t\right)\hfill & g\ne q.\hfill \end{array}\right.$$

(48)

Then, the E- and M-steps at the qth cycle of the rth iteration are introduced below.

5.1. Deterministic Signal Model

Based on (48), we have ${\mathbf{h}}_q\left(t\right)\sim \mathcal{CN}\left(\mathbf{b}\left({\xi }_q\right)m_q\left(t\right),\tau {\mathbf{I}}_K\right)$ and the ${\mathbf{h}}_g\left(t\right)$’s with $g\ne q$ are deterministic. The complete-data LLF is expressed as

$$\begin{array}{ccc}\hfill \mathcal{Z}({\xi }_q,{\mathbf{m}}_q,\tau )& =& \sum _{t=1}^{L}logp\left({\mathbf{h}}_q\left(t\right);{\xi }_q,m_q\left(t\right),\tau \right)\hfill \\ & =& -LKlog\left(\pi \tau \right)-\frac{1}{\tau }\sum _{t=1}^L{\parallel {\mathbf{h}}_q\left(t\right)-\mathbf{b}\left({\xi }_q\right)m_q\left(t\right)\parallel }^2.\hfill \end{array}$$

(49)

5.1.1. E-Step

Calculate the conditional expectation of the complete-data LLF in (49)

$$\begin{array}{c}\hfill \mathcal{E}\left\{\mathcal{Z}({\xi }_q,{\mathbf{m}}_q,\tau )|\mathbf{W};{\mathbf{\Delta }}^{(r,q-1)}\right\}=-LKlog\left(\pi \tau \right)-\frac{1}{\tau }\sum _{t=1}^L{\parallel {\mathbf{h}}_q^{\left(r\right)}\left(t\right)-\mathbf{b}\left({\xi }_q\right)m_q\left(t\right)\parallel }^2,\end{array}$$

(50)

where $\mathbf{\Delta }=(\mathit{\xi },\mathbf{M},\tau )$, ${\mathbf{0}}_K$ is the $K\times K$ zero matrix, and 

$$\begin{array}{ccc}\hfill {\mathbf{h}}_q^{\left(r\right)}\left(t\right)& =& {\mathbf{h}}_q^{(r,q)}\left(t\right)=\mathcal{E}\left\{{\mathbf{h}}_q\left(t\right)|\mathbf{W};{\mathbf{\Delta }}^{(r,q-1)}\right\}\hfill \\ & =& \mathbf{b}\left({\xi }_q^{(r,q-1)}\right)m_q^{(r,q-1)}\left(t\right)+\left[\mathbf{w}\left(t\right)-\mathbf{B}\left({\xi }^{(r,q-1)}\right){\mathbf{m}}^{(r,q-1)}\left(t\right)\right],\hfill \end{array}$$

(51)

$$\begin{array}{ccc}& & \mathcal{D}\left\{{\mathbf{h}}_q\left(t\right)|\mathbf{W};{\mathbf{\Delta }}^{(r,q-1)}\right\}={\mathbf{0}}_K.\hfill \end{array}$$

(52)

5.1.2. M-Step

Update the estimates of ${\xi }_q$, ${\mathbf{m}}_q$, and $\tau$ by solving

$$\begin{array}{c}\hfill \underset{{\xi }_q\in \Omega ,{\mathbf{m}}_q,\tau >0}{min}Klog\left(\tau \right)+\frac{1}{\tau L}\sum _{t=1}^L{\parallel {\mathbf{h}}_q^{\left(r\right)}\left(t\right)-\mathbf{b}\left({\xi }_q\right)m_q\left(t\right)\parallel }^2.\end{array}$$

(53)

${\xi }_q^{\left(r\right)}$, ${\mathbf{m}}_q^{\left(r\right)}$, and ${\tau }^{(r,q)}$ are obtained by

$$\begin{array}{ccc}\hfill {\xi }_q^{\left(r\right)}& =& {\xi }_q^{(r,q)}=arg\underset{{\xi }_q\in \Omega }{max}\mathrm{Tr}\left\{{\mathbf{\Gamma }}_q{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right\},\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill m_q^{\left(r\right)}\left(t\right)& =& m_q^{(r,q)}\left(t\right)={\mathbf{b}}^H\left({\xi }_q^{\left(r\right)}\right){\mathbf{h}}_q^{\left(r\right)}\left(t\right)/K,\forall t,\hfill \end{array}$$

(55)

$$\begin{array}{ccc}\hfill {\tau }^{(r,q)}& =& d_q^{\left(r\right)}/K,\hfill \end{array}$$

(56)

where ${\widehat{\mathbf{P}}}_q^{\left(r\right)}=(1/L){\sum }_{t=1}^L{\mathbf{h}}_q^{\left(r\right)}\left(t\right){\left[{\mathbf{h}}_q^{\left(r\right)}\left(t\right)\right]}^H$ and $d_q^{\left(r\right)}=\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_q^{\left(r\right)}){\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)$.

Moreover, the other parameter estimates are not updated at this cycle and their iterative values are

$$\begin{array}{ccc}\hfill {\xi }_g^{(r,q)}& =& {\xi }_g^{(r,q-1)},\forall g\ne q,\hfill \end{array}$$

(57)

$$\begin{array}{ccc}\hfill m_g^{(r,q)}\left(t\right)& =& m_g^{(r,q-1)}\left(t\right),\forall g\ne q,t.\hfill \end{array}$$

(58)

Remark 2.

Since the ${\mathbf{h}}_q^{\left(r\right)}\left(t\right)$’s in (51), ${\xi }_q^{\left(r\right)}$ in (54), and the $m_q^{\left(r\right)}\left(t\right)$’s in (55) are unrelated to ${\tau }^{(r,q-1)}$, we can omit (56) due to the nuisance parameter τ.

5.2. Random Signal Model

Based on (48), we have ${\mathbf{h}}_q\left(t\right)\sim \mathcal{CN}(\mathbf{0},{\mathbf{N}}_q)$ with ${\mathbf{N}}_q={\rho }_q\mathbf{b}\left({\xi }_q\right){\mathbf{b}}^H\left({\xi }_q\right)+\tau {\mathbf{I}}_K$. However, the distribution of ${\mathbf{h}}_g\left(t\right)$ with $g\ne q$ is only associated with $m_g\left(t\right)$. The complete-data LLF is written as

$$\begin{array}{ccc}\hfill \mathcal{Z}({\xi }_q,\rho ,\tau )& =& \sum _{t=1}^L\sum _{g\ne q}logp\left(m_g\left(t\right);{\rho }_g\right)+\sum _{t=1}^{L}logp\left({\mathbf{h}}_q\left(t\right);{\xi }_q,{\rho }_q,\tau \right)\hfill \\ & =& -L(G-1)log\left(\pi \right)-L\sum _{g\ne q}\left[log\left({\rho }_g\right)+{\widehat{P}}_g/{\rho }_g\right]\hfill \\ & & -LKlog\left(\pi \right)-L\left[log\left(\mathrm{Det}\left({\mathbf{N}}_q\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_q^{-1}{\widehat{\mathbf{P}}}_q\right)\right],\hfill \end{array}$$

(59)

where ${\widehat{P}}_g=(1/L){\sum }_{t=1}^L{\left|m_g\left(t\right)\right|}^2$ and $|·|$ denotes modulus.

5.2.1. E-Step

Calculate the conditional expectation of the complete-data LLF in (59)

$$\begin{array}{ccc}\hfill \mathcal{E}\left\{\mathcal{Z}({\xi }_q,\mathit{\rho },\tau )|\mathbf{W};{\mathbf{\Delta }}^{(r,q-1)}\right\}& =& V-L\sum _{g\ne q}\left[log\left({\rho }_g\right)+{\widehat{P}}_g^{(r,q)}/{\rho }_g\right]\hfill \\ & & -L\left[log\left(\mathrm{Det}\left({\mathbf{N}}_q\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_q^{-1}{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)\right],\hfill \end{array}$$

(60)

where $\mathbf{\Delta }=(\mathit{\xi },\mathit{\rho },\tau )$ and $V=-L(K+G-1)log\left(\pi \right)$. In (60),

$$\begin{array}{ccc}\hfill {\widehat{P}}_g^{(r,q)}& =& \mathcal{E}\left\{{\widehat{P}}_g|\mathbf{W};{\mathbf{\Delta }}^{(r,q-1)}\right\}\hfill \\ & =& {\rho }_g^{(r,q-1)}\left[1-{\mathbf{b}}^H\left({\xi }_g^{(r,q-1)}\right){\mathbf{d}}_g^{(r,q-1)}\right]+{\left[{\mathbf{d}}_g^{(r,q-1)}\right]}^H{\widehat{\mathbf{P}}}_w{\mathbf{d}}_g^{(r,q-1)}\ge 0\hfill \end{array}$$

(61)

with ${\mathbf{d}}_g^{(r,q-1)}={\left[{\mathbf{N}}_w^{(r,q-1)}\right]}^{-1}\mathbf{b}\left({\xi }_g^{(r,q-1)}\right){\rho }_g^{(r,q-1)}$ and

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{P}}}_q^{\left(r\right)}& =& {\widehat{\mathbf{P}}}_q^{(r,q)}=\mathcal{E}\left\{{\widehat{\mathbf{P}}}_q|\mathbf{W};{\mathbf{\Delta }}^{(r,q-1)}\right\}\hfill \\ & =& {\mathbf{N}}_q^{(r,q-1)}-{\mathbf{N}}_q^{(r,q-1)}{\left[{\mathbf{N}}_w^{(r,q-1)}\right]}^{-1}{\mathbf{N}}_q^{(r,q-1)}\hfill \\ & & +{\mathbf{N}}_q^{(r,q-1)}{\left[{\mathbf{N}}_w^{(r,q-1)}\right]}^{-1}{\widehat{\mathbf{P}}}_w{\left[{\mathbf{N}}_w^{(r,q-1)}\right]}^{-1}{\mathbf{N}}_q^{(r,q-1)}.\hfill \end{array}$$

(62)

5.2.2. M-Step

Update the estimates of ${\xi }_q$, $\mathit{\rho }$, and $\tau$ by solving

$$\begin{array}{c}\hfill \underset{{\xi }_q\in \Omega ,\rho \ge \mathbf{0},\tau >0}{min}\sum _{g\ne q}\left[log\left({\rho }_g\right)+{\widehat{P}}_g^{(r,q)}/{\rho }_g\right]+log\left(\mathrm{Det}\left({\mathbf{N}}_q\right)\right)+\mathrm{Tr}\left({\mathbf{N}}_q^{-1}{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right).\end{array}$$

(63)

We, thus, have ${\rho }_g^{(r,q)}={\widehat{P}}_g^{(r,q)},\forall g\ne q$, and ${\xi }_q^{\left(r\right)}$, ${\rho }_q^{(r,q)}$, and ${\tau }^{(r,q)}$ are obtained by solving

$$\begin{array}{c}\hfill \underset{{\xi }_q\in \Omega ,{\sigma }_q\ge 0,\tau >0}{min}Klog\left(\tau \right)+log\left(\mathrm{Det}\left({\mathbf{Q}}_q\right)\right)+\frac{1}{\tau }\mathrm{Tr}\left({\mathbf{Q}}_q^{-1}{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right),\end{array}$$

(64)

where ${\mathbf{N}}_q=\tau {\mathbf{Q}}_q$ with ${\mathbf{Q}}_q={\sigma }_q\mathbf{b}\left({\xi }_q\right){\mathbf{b}}^H\left({\xi }_q\right)+{\mathbf{I}}_K$ and ${\sigma }_q={\rho }_q/\tau$. Following (45)–(47), ${\xi }_q^{\left(r\right)}$, ${\rho }_q^{(r,q)}$, and ${\tau }^{(r,q)}$ are obtained by

$$\begin{array}{ccc}\hfill {\xi }_q^{\left(r\right)}& =& {\xi }_q^{(r,q)}=arg\underset{{\xi }_q\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_q{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right),\hfill \end{array}$$

(65)

$$\begin{array}{ccc}\hfill {\tau }^{(r,q)}& =& \left\{\begin{array}{cc}{d}_q^{\left(r\right)}/(K-1)\ge 0,\hfill & e_q^{\left(r\right)}>\mathrm{Tr}\left({\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)/K,\hfill \\ \mathrm{Tr}\left({\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)/K\ge 0,\hfill & e_q^{\left(r\right)}\le \mathrm{Tr}\left({\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)/K,\hfill \end{array}\right.\hfill \end{array}$$

(66)

$$\begin{array}{ccc}\hfill {\rho }_q^{(r,q)}& =& max\left\{\left(e_q^{\left(r\right)}-\mathrm{Tr}\left({\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)/K\right)/(K-1),0\right\},\hfill \end{array}$$

(67)

where ${\tau }^{(r,q)}=0$ is possible although its probability is very low. For example, if $G=2$, $L=1$, and ${\rho }_1^{(1,1)}=0$ at 1st cycle of the 1st iteration, at the 2nd cycle of the 1st iteration we will have

$$\begin{array}{cc}\hfill & {\mathbf{N}}_w^{(1,1)}={\rho }_2^{(1,1)}\mathbf{b}\left({\xi }_2^{(1,1)}\right){\mathbf{b}}^H\left({\xi }_2^{(1,1)}\right)+{\tau }^{(1,1)}{\mathbf{I}}_K={\mathbf{N}}_2^{(1,1)}\end{array}$$

(68)

and then ${\widehat{\mathbf{P}}}_2^{\left(1\right)}=\mathbf{w}\left(1\right){\mathbf{w}}^H\left(1\right)$ by (62). Furthermore, if $\mathbf{w}\left(1\right)=h\mathbf{b}\left(\overline{\xi }\right)$$(h\ne 0,\overline{\xi }\in \Omega )$, we will obtain ${\xi }_2^{\left(1\right)}=\overline{\xi }$ by (65) and ${\tau }^{(1,2)}=d_2^{\left(1\right)}/(K-1)=0$ by (66).

To avoid ${\tau }^{(r,q)}=0$, we use the following two CM-steps to re-estimate ${\xi }_q$, ${\rho }_q$, and $\tau$ based on problem (64) if ${\tau }^{(r,q)}=0$ in (66).

• First CM-step: Estimate ${\xi }_q$ and ${\rho }_q$ but hold $\tau ={\tau }^{(r,q-1)}$ fixed. Then, problem (64) is simplified to

  $$\begin{array}{c}\hfill \underset{{\xi }_q\in \Omega ,{\sigma }_q\ge 0}{min}log\left(\mathrm{Det}\left({\mathbf{Q}}_q\right)\right)+\frac{1}{{\tau }^{(r,q-1)}}\mathrm{Tr}\left({\mathbf{Q}}_q^{-1}{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right).\end{array}$$

  (69)
  Following (21) and (22), ${\xi }_q^{\left(r\right)}$ and ${\rho }_q^{(r,q)}$ are obtained by

  $$\begin{array}{ccc}\hfill {\xi }_q^{\left(r\right)}& =& {\xi }_q^{(r,q)}=arg\underset{{\xi }_q\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_q{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right),\hfill \end{array}$$

  (70)

  $$\begin{array}{ccc}\hfill {\rho }_q^{(r,q)}& =& {\sigma }_q^{(r,q)}{\tau }^{(r,q-1)}=max\left\{\left(e_q^{\left(r\right)}-{\tau }^{(r,q-1)}\right)/K,0\right\}.\hfill \end{array}$$

  (71)
• Second CM-step: Estimate $\tau$ but hold ${\xi }_q={\xi }_q^{\left(r\right)}$ and ${\sigma }_q={\sigma }_q^{(r,q)}$ fixed. Then, problem (64) is simplified to

  $$\begin{array}{c}\hfill \underset{\tau >0}{min}Klog\left(\tau \right)+\frac{1}{\tau }\mathrm{Tr}\left({\left[{\mathbf{Q}}_q^{\left(r\right)}\right]}^{-1}{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right),\end{array}$$

  (72)

  where ${\mathbf{Q}}_q^{\left(r\right)}={\sigma }_q^{(r,q)}\mathbf{b}\left({\xi }_q^{\left(r\right)}\right){\mathbf{b}}^H\left({\xi }_q^{\left(r\right)}\right)+{\mathbf{I}}_K$. We obtain ${\tau }^{(r,q)}$ by

  $$\begin{array}{cc}\hfill & {\tau }^{(r,q)}=\mathrm{Tr}\left\{{\left[{\mathbf{Q}}_q^{\left(r\right)}\right]}^{-1}{\widehat{\mathbf{P}}}_q^{\left(r\right)}\right\}/K=\left({\tau }^{(r,q-1)}+d_q^{\left(r\right)}\right)/K,\end{array}$$

  (73)

  where $d_q^{\left(r\right)}=\mathrm{Tr}\left(({\mathbf{I}}_K-{\mathbf{\Gamma }}_q^{\left(r\right)}){\widehat{\mathbf{P}}}_q^{\left(r\right)}\right)\ge 0$ and ${\tau }^{(r,q)}>0$ if ${\tau }^{(r,q-1)}>0$.

Moreover, the other parameter estimate(s) is (are) not updated at this cycle and the iterative value(s) is (are)

$$\begin{array}{cc}\hfill & {\xi }_g^{(r,q)}={\xi }_g^{(r,q-1)},\forall g\ne q.\end{array}$$

(74)

6. Properties of the Proposed EM, MEM, and SAGE Algorithms

6.1. Convergence Point

It is easy to verify that the above EM-type algorithms satisfy standard regularity conditions [5,21,26] and always converge to stationary points of $\mathcal{J}(\mathbf{\Psi },\tau )$. Of course, the convergence points of these EM-type algorithms depend on their initial points. To generate appropriate initial points, we can employ the method presented in [10] using the deterministic signal model.

6.2. Complexity and Stability

At the rth iteration, the computational burdens of the above EM-type algorithms lie in solving the G maximization problems

$$\begin{array}{c}\hfill {\xi }_g^{\left(r\right)}=arg\underset{{\xi }_g\in \Omega }{max}\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g.\end{array}$$

(75)

Thus, these EM-type algorithms have almost the same computational complexity at every iteration and if an algorithm has faster convergence, its number of iterations required will be smaller and its computations will be fewer.

However, when the powers of sources are unequal, we have found via simulation that the DOA estimates of multiple sources, updated by (75), tend to be consistent with the true DOA of the source with the largest power. Accordingly, these EM-type algorithms may be unstable. To address this issue, we can reduce the difference between ${\xi }_g^{\left(r\right)}$ and ${\xi }_g^{(r-1)}$ by increasing $\mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)$ rather than maximizing it, i.e.,

$$\begin{array}{c}\hfill \mathrm{Tr}\left({\mathbf{\Gamma }}_g^{\left(r\right)}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)\ge \mathrm{Tr}\left({\mathbf{\Gamma }}_g^{(r-1)}{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right),\forall g,\end{array}$$

(76)

which guarantees the monotonicity of these EM-type algorithms [3]. As a good example, Algorithm 1 in the next section has given excellent simulation results.

Algorithm 1 Gradient ascent with backtracking line search

1: $y\left({\eta }_g\right)=K\times \mathrm{Tr}\left({\mathbf{\Gamma }}_g{\widehat{\mathbf{P}}}_g^{\left(r\right)}\right)$, ${\eta }_g={\eta }_g^{(r-1)}\in (0,\pi )\left(\mathrm{radian}\right)$. 2: while  $|y^{\prime }\left({\eta }_g\right)|>0.001$  do 3:    $s=0.1\times \left\{\begin{array}{cc}(\pi -{\eta }_g)/y^{\prime }\left({\eta }_g\right),\hfill & y^{\prime }\left({\eta }_g\right)>0,\hfill \\ -{\eta }_g/y^{\prime }\left({\eta }_g\right),& y^{\prime }\left({\eta }_g\right)<0.\end{array}\right.$ 4:    while $y\left({\eta }_g+s{y}^{\prime }\left({\eta }_g\right)\right)<y\left({\eta }_g\right)+0.3s{\left(y^{\prime }\left({\eta }_g\right)\right)}^2$ do 5:      $s=0.5s$. 6:    end while 7:    ${\eta }_g={\eta }_g+s{y}^{\prime }\left({\eta }_g\right)\in (0,\pi )\left(\mathrm{radian}\right)$. 8: end while 9: ${\eta }_g^{\left(r\right)}={\eta }_g$.

7. Simulation Results

We give simulation results to compare the proposed EM-type algorithms. For convenience, the array is assumed to be a uniformly spaced linear array, in which ${\mathit{\zeta }}_k={\left[\frac{\chi }{2}(k-1)00\right]}^T$. We set $G=2$, $L=20$, and $\beta =\mathbf{1}/G$ in the EM algorithm. Here, ${\mu }_1={\mu }_2=90°$ is known while ${\eta }_1$ and ${\eta }_2$ need to be estimated. $\mathbf{M}$ in the deterministic signal model is also generated by the independent random numbers $m_g\left(t\right)\sim \mathcal{CN}(0,{\rho }_g)$. All the algorithms adopt the stopping criterion $\parallel {\mathit{\xi }}^{(r+1)}-{\mathit{\xi }}^{\left(r\right)}\parallel \le 0.001°$. Algorithm 1 is designed to search the ${\eta }_g^{\left(r\right)}$’s in (76) [27]. Moreover, ${\mathbf{M}}^{\left(0\right)}={\left[\mathbf{1}\mathbf{1}\right]}^T$, ${\mathit{\rho }}^{\left(0\right)}=\mathbf{1}$, ${\mathit{\tau }}^{\left(0\right)}=\mathbf{1}/G$, and ${\tau }^{\left(0\right)}=1$. In this section, the EM, MEM, and SAGE algorithms are simply written as EM, MEM, and SAGE, respectively.

7.1. Deterministic Signal Model

To compare the convergence of EM, MEM, and SAGE, Figure 1 plots their $\mathcal{J}({\mathbf{\Psi }}^{\left(r\right)},{\tau }^{\left(r\right)})$’s, ${\eta }_1^{\left(r\right)}$’s, and ${\eta }_2^{\left(r\right)}$’s under one trial. In Figure 1, $K=10$, ${\eta }_1=20°$, ${\eta }_2=80°$, ${\rho }_1=-2\mathrm{dB}$, ${\rho }_2=4\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=24°$, and ${\eta }_2^{\left(0\right)}=84°$. It is easy to see that EM, MEM, and SAGE converge to a consistent $({\eta }_1,{\eta }_2)$ estimate given an accurate initial point. Moreover, EM has a similar convergence with MEM while SAGE converges faster than EM and MEM.

Figure 2 and Figure 3 show two scatter plots of $({\eta }_1,{\eta }_2)$ estimates under 200 trials. In Figure 2, $K=10$, ${\eta }_1=25°$, ${\eta }_2=75°$, ${\rho }_1=-4\mathrm{dB}$, ${\rho }_2=2\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=40°$, and ${\eta }_2^{\left(0\right)}=60°$. Moreover, the numbers of desirable points obtained by EM, MEM, and SAGE are 68, 72, and 179, respectively. In Figure 3, $K=10$, ${\eta }_1=70°$, ${\eta }_2=78°$, ${\rho }_1=-2\mathrm{dB}$, ${\rho }_2=4\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=50°$, and ${\eta }_2^{\left(0\right)}=58°$. Moreover, the numbers of desirable points obtained by EM, MEM, and SAGE are 159, 157, and 190, respectively. Figure 2 and Figure 3 imply that given a poor initial point, SAGE can avoid the convergence to an undesirable stationary point of $\mathcal{J}(\mathbf{\Psi },\tau )$ more efficiently than EM and MEM.

Note that in each of Figure 2 and Figure 3, the SAGE algorithm requires the least processing time due to the fastest convergence and thus performs the fewest computations required. Moreover, both sources in Figure 2 are not closely spaced, so it is very difficult to mix up both sources and the desirable points in Figure 2 are centered around the true position $(25°,75°)$. However, both sources in Figure 3 are closely spaced and the desirable points are centered around $(78°,70°)$ or $(70°,78°)$, i.e., these EM-type algorithms are likely to mix up closely spaced sources.

According to these simulations, we can conclude that (1) EM has similar convergence with MEM, and (2) SAGE outperforms EM and MEM.

7.2. Random Signal Model

To compare the convergence of EM, MEM, and SAGE, Figure 4 plots their $\mathcal{J}({\mathbf{\Psi }}^{\left(r\right)},{\tau }^{\left(r\right)})$’s, ${\eta }_1^{\left(r\right)}$’s, and ${\eta }_2^{\left(r\right)}$’s under one trial. In Figure 4, $K=10$, ${\eta }_1=20°$, ${\eta }_2=80°$, ${\rho }_1=-4\mathrm{dB}$, ${\rho }_2=4\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=24°$, and ${\eta }_2^{\left(0\right)}=84°$. We can also observe that given an accurate initial point, EM, MEM, and SAGE converge to a consistent $({\eta }_1,{\eta }_2)$ estimate. Moreover, EM has similar convergence with MEM while SAGE converges faster than EM and MEM.

Figure 5 and Figure 6 show two scatter plots of $({\eta }_1,{\eta }_2)$ estimates under 200 trials. In Figure 5, $K=10$, ${\eta }_1=25°$, ${\eta }_2=75°$, ${\rho }_1=-4\mathrm{dB}$, ${\rho }_2=2\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=40°$, and ${\eta }_2^{\left(0\right)}=60°$. Moreover, the numbers of desirable points obtained by EM, MEM, and SAGE are 185, 186, and 175, respectively. In Figure 6, $K=10$, ${\eta }_1=70°$, ${\eta }_2=78°$, ${\rho }_1=-2\mathrm{dB}$, ${\rho }_2=-1\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=55°$, and ${\eta }_2^{\left(0\right)}=63°$. Moreover, the numbers of desirable points obtained by EM, MEM, and SAGE are 161, 161, and 172, respectively. Figure 5 and Figure 6 imply that EM has similar convergence with MEM but compared to EM and MEM, SAGE is less and more efficient for avoiding the convergence to an undesirable stationary point of $\mathcal{J}(\mathbf{\Psi },\tau )$ in Figure 5 and Figure 6, respectively. Note that in each of Figure 5 and Figure 6, the SAGE algorithm requires the least processing time due to the fastest convergence and thus performs the fewest computations required. In addition, these EM-type algorithms are likely to mix up closely spaced sources, so the desirable points in Figure 5 are centered around $(25°,75°)$ and the desirable points in Figure 6 are centered around $(78°,70°)$ or $(70°,78°)$.

Based on the above figures, we can conclude that (1) EM has similar convergence with MEM, and (2) SAGE cannot always outperform EM and MEM.

7.3. Deterministic and Random Signal Models

Snapshots from the random signal model can be processed by these EM-type algorithms for the deterministic signal model, which means that we can compare these algorithms for both signal models. The above simulation results have shown that EM has similar convergence with MEM, so we only compare EM and SAGE for both signal models in this subsection for simplicity.

Since both signal models have the same DOA parameter $\mathit{\xi }$, the stopping criterion $\parallel {\mathit{\xi }}^{(r+1)}-{\mathit{\xi }}^{\left(r\right)}\parallel \le 0.001°$ is suitable. Figure 7 shows a scatter plot of $({\eta }_1,{\eta }_2)$ estimates under 50 trials. In Figure 7, $K=10$, ${\eta }_1=50°$, ${\eta }_2=100°$, ${\rho }_1=-4\mathrm{dB}$, ${\rho }_2=4\mathrm{dB}$, $\tau =4\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=55°$, and ${\eta }_2^{\left(0\right)}=95°$. We can see that EM and SAGE for both signal models yield close $({\eta }_1,{\eta }_2)$ estimates for each trial.

Based on Figure 7, Figure 8 compares the numbers of iterations required by these algorithms. We can observe that EM for the deterministic signal model generally needs more iterations than EM for the random signal model. The reason is that EM for deterministic signal model needs to update more parameter estimates at each iteration and, thus, has slower convergence than EM for the random signal model. Moreover, SAGE for the deterministic signal model generally needs fewer iterations than SAGE for the random signal model. More importantly, SAGE for the deterministic signal model always requires the fewest iterations for each trial. Thus, we can conclude that SAGE for the deterministic signal model is superior to the other algorithms in the computational cost.

Figure 9 shows the root mean square error (RMSE) performances of DOA estimation obtained by the SAGE algorithm for the deterministic and random signal models. In Figure 9, ${\eta }_1=60°$, ${\eta }_2=120°$, ${\rho }_1=0\mathrm{dB}$, ${\rho }_2=1\mathrm{dB}$, $\tau =3\mathrm{dB}$, ${\eta }_1^{\left(0\right)}=55°$, and ${\eta }_2^{\left(0\right)}=115°$. Each RMSE is computed from 1000 independent realizations. We can observe that as the number of sensors K increases, the SAGE algorithm for each signal model yields small RMSEs, which indicates that increasing the number of sensors K can improve the accuracy of DOA estimation.

8. Conclusions

In this paper, we applied and designed the EM and SAGE algorithms for DOA estimation in unknown uniform noise and proposed a new MEM algorithm applicable to the noise assumption. Next, we improved these EM-type algorithms to ensure the stability when the powers of sources are unequal. After being improved, simulation results illustrated that the EM algorithm has similar convergence with the MEM algorithm, the SAGE algorithm outperforms the EM and MEM algorithms for the deterministic signal model, and the SAGE algorithm cannot always outperform the EM and MEM algorithms for the random signal model. In addition, simulation results indicated that when these EM-type algorithms process the same snapshots from the random signal model, the SAGE algorithm for the deterministic signal model can require the fewest iterations and computations.