<figure><img src="image_1.png" /></figure>

HAL
open science

# An Approximate Message Passing Approach for DOA Estimation in Phase Noisy Environments

Guillaume Beaumont, Ronan Fablet, Angélique Drémeau

► To cite this version:

Guillaume Beaumont, Ronan Fablet, Angélique Drémeau. An Approximate Message Passing Approach for DOA Estimation in Phase Noisy Environments. LVA/ICA 2018 - international conference on Latent Variable Analysis and Signal Separation, Jun 2018, Surrey, United Kingdom. pp.385-394, 10.1007/978-3-319-93764-9_36. hal-05352375

HAL Id: hal-05352375

https://hal.science/hal-05352375v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An Approximate Message Passing approach for DOA estimation in phase noisy environments

Guillaume Beaumont¹,*, Ronan Fablet², and Angélique Drémeau¹

¹Lab-STICC UMR 6285, CNRS, ENSTA Bretagne, Brest, F-29200, France

²IMT-Atlantique, Lab-STICC UMR 6285, UBL, Brest F-29200, France

**Abstract.** In underwater acoustics, wave propagation can be greatly disrupted by random fluctuations in the ocean environment. In particular, phase measurements of the complex pressure field can be heavily noisy and can defeat conventional direction-of-arrival (DOA) estimation algorithms.

In this paper, we propose a new Bayesian approach to address such phase noise through an informative prior on the measurements. This is combined to a sparse assumption on the directions of arrival to achieve a highly-resolved estimation and integrated into a message-propagation algorithm referred to as the “paSAMP” algorithm (for Phase-Aware Swept Approximate Message Passing). Our algorithm can be seen as an extension of the recent phase-retrieval algorithm “prSAMP” to phase-aware priors.

Experiments on simulated data mimicking real environments demonstrate that paSAMP outperform conventional approaches (e.g. classic beamforming) in terms of DOA estimation. paSAMP also proves to be more robust to additive noise than other variational methods (e.g. based on mean-field approximation).

**Keywords:** DOA estimation, sparse representation, Bayesian estimation, variational Bayesian approximations, message passing algorithms

## 1 Introduction

Common to many applications such as SONAR, RADAR, and telecommunications, direction-of-arrival (DOA) estimation aims at locating one or more sources emitting in some propagation media. Various methods have been proposed to address this problem. They can be distinguished by the assumptions made on the propagating medium and sources.

The beamforming approach [1] constitutes the most famous approach. As it implicitly assumes the noise to be Gaussian and additive, it leads to poor estimation performance for compelex phase perturbations. The so-called “high-resolution” techniques consider additional assumptions over the number or the nature of the sources. This is the case of the well-known MUSIC method [2].

<footnote>* This work has been supported by the DGA MRIS.</footnote>
MUSIC assumes the number of sources to be known and the separability of the sub-spaces where the noise and the signal live. More recently, “compressive” beamforming approaches proposed e.g. in [3] benefit from an explicit sparse model on the sources.

While all the previously cited approaches rely on an additive Gaussian noise model, recent work has focused on the integration of phase-noise models better accounting for complex propagation processes. Such approaches aim to take into account the wave-front distortion occurring when waves travel through fluctuat- ing media. This is of key interest for a wide range of application fields including as underwater acoustics [4, 5] or atmospheric sound propagation [6, 7]. These contributions mainly relate to recent advances in phase recovery (see e.g. [8–11]) and the use of informative priors on the missing phases. In this respect, we can mention the Bayesian approach “paVBEM” based on a mean-field approxima- tion [12].

Here, we further explore a variational Bayesian approach. Knowing that higher-order approximations and associated message-passing algorithms outper- form mean-field approximations for a wide range of inverse problems [13], we propose a novel approach based on the “swept approximate message passing” (SwAMP) algorithm introduced in [14]. Our algorithm is proven to be more ro- bust to additive noise and multiplicative phase noise than previous approaches using phase-aware priors such as the paVBEM approach [12] and those using non-informative phase priors [9].

## 2 PROBLEM STATEMENT

In this section, we recall the Bayesian modeling introduced in [12], which we shall follow throughout of this paper, and introduce the estimation problem we propose to solve.

### 2.1 Observation Model

Our objective is to design an algorithm able to recover the directions of arrival of a few waves, despite a structured phase-noisy environment, exploiting one single temporal snapshot on a uniform linear sensor array. In underwater acoustics, this noise is mainly due to internal waves, changing the local sound-speed (see e.g. [4]). These internal waves and their impact on the acoustic measurements have been studied in different works (see [4, 5]), which leads to a statistical characterization of the phase noise.

In this context, we propose the following observation model: we consider a linear antenna composed of N regularly-spaced sensors and assume that the received signal at sensor n can be expressed as

$$y_n = e^{j\theta_n} \sum_{m=1}^{M} d_{nm}x_m + \omega_n, \quad (1)$$
where $\theta_n$ stands for the phase noise due to the propagation through the fluctuating medium and $\omega_n$ an additive noise. The scalar $d_{nm}$ is the n-th element of the steering vector $\mathbf{d}_m = [e^{j\frac{2\pi}{\lambda}\Delta\sin(\phi_m)} \dots e^{j\frac{2\pi}{\lambda}\Delta^N\sin(\phi_m)}]^T$ where the $\phi_m$'s are some potential angles of arrival, $\Delta$ is the distance between two adjacent sensors and $\lambda$ is the wavelength of the propagation waves.

Within model (1), at each sensor of the antenna, we assume that the received field is a combination of a few waves arriving from different angles $\phi_m$. The DOA estimation problem then consists in identifying the positions of the non-zero coefficients in $\mathbf{x} \triangleq [x_1 \dots x_M]^T$. In underwater acoustics, the phase noise considered in (1) is well-suited to characterize phase perturbations of the wave front in a fluctuating ocean [5], especially in the case of the so-called “partially saturated” propagation regime defined in [4]. This regime focuses on far-field propagation at high frequency with no multipath. In this case, amplitude variations of the measured acoustic field can be neglected regarding the high sensitivity to a structured phase-noise. Note that a similar fluctuation regime has been also identified in atmospheric sound propagation (see [7]).

## 2.2 Bayesian formulation of the problem

We address the estimation of $\mathbf{x}$ from the measurements $\mathbf{y} \triangleq [y_1, \dots, y_N]^T$ in the presence of (unknown) additive noise $\omega \triangleq [\omega_1, \dots, \omega_N]^T$ and multiplicative phase noise $\theta \triangleq [\theta_1, \dots, \theta_N]^T$. To solve this problem, we consider a Bayesian framework and define some prior assumptions on the different variables in (1).

A first assumption is set on the number of sources (i.e. the non-zero coefficients in $\mathbf{x}$) that we suppose to be small in front of the number of sensors. To take into account this sparse hypothesis, we adopt a Bernoulli-Gaussian model $\forall m \in \{1, \dots, M\}$

$$p(x_m) = \rho \mathcal{CN}(x_m; m_x, \sigma_x^2) + (1 - \rho)\delta_0(x_m), \quad (2)$$

where $\rho$ is the Bernoulli parameter and equals the probability for $m_x$ to be non-zero<fnref n="1" />, $\mathcal{CN}(x_m; m_x, \sigma_x^2)$ stands for the circular complex Gaussian distribution with mean $m_x$ and variance $\sigma_x^2$, and $\delta_0(x_m)$ for the Dirac distribution. The Bernoulli-Gaussian model is widely used when considering Bayesian inference methods for sparsity-constrained problems (see e.g. [15, 16]).

Previous studies of the statistical characterization of fluctuation phenomena [4, 5] provide the basis for the definition of a phase-noise prior. In underwater acoustics, [4, 5] exhibited and characterized the existence of a spatial correlation of the measured field all along the antenna. To account for the resulting coherence length, we consider a Markovian model as.

$$p(\theta_n|\theta_{n-1}) = \mathcal{N}(\theta_n; \beta \theta_{n-1}, \sigma_\theta^2), \quad \forall n \in \{2, \dots, N\}, \quad (3)$$

$$p(\theta_1) = \mathcal{N}(\theta_1; 0, \sigma_1^2), \quad (4)$$

<footnote n="1">¹ We assume the Bernoulli parameter to be the same for each $m \in \{1, \dots, M\}$.</footnote>
with $\beta \in \mathbb{R}_+$. Variance $\sigma_{\theta}^2$ is related to the coherence length and accounts for the strength of the fluctuations. As an example, a large $\sigma_{\theta}^2$ models strong fluctuations of the medium and results in a small coherence length, such that the phase noise varies widely from a sensor to the neighboring ones.

We also introduce an additive noise $\omega$ to account for the combination of a large number of random parasitic phenomena. Based on the central limit theorem, we consider with a classic zero-mean Gaussian distribution with variance $\sigma^2$.

Overall, our Bayesian formulation leads to the following Minimum Mean Square Error (MMSE) problem:

$$ \hat{\mathbf{x}} = \underset{\tilde{\mathbf{x}}}{\operatorname{argmin}} \int_{\mathbf{x}} \|\mathbf{x} - \tilde{\mathbf{x}}\|_2^2 p(\mathbf{x}|\mathbf{y})d\mathbf{x} \quad (5) $$

where $p(\mathbf{x}|\mathbf{y}) = \int_{\theta} p(\mathbf{x}, \theta|\mathbf{y})d\theta$.

To solve efficiently this problem, we propose to exploit a variational Bayesian inference strategy, that approximates the posterior joint distribution $p(\mathbf{x}, \theta|\mathbf{y})$ by a distribution having a suitable factorization. In [12], a mean-field approximation $p(\mathbf{x}, \theta|\mathbf{y}) \approx q(\theta) \prod_{m=1}^{M} q(x_m)$ was considered. Here, we address a different type of factorization, called the Bethe approximation, relating to the “approximate message passing” (AMP) algorithms [13]. This approximation exploits higher-order terms which result in better estimation performance [13].

We motivate and detail our approach in the next section.

# 3 The “paSAMP” algorithm

In this section, we motivate and present the novel algorithm proposed to solve problem (5).

## 3.1 Motivation and main principles of the approach

AMP algorithms have been considered for a few years as a serious solution to linear problems under sparsity constraints. First considered in the sole case of i.i.d (sub-)Gaussian matrices, they have been recently extended to random but more general matrices by the “vector approximate message passing” (VAMP) algorithm [17] and to highly correlated matrices by the “swept approximate message passing” (SwAMP) approach [14]. Both methods aim at alleviating the convergence issues of AMP (notably highlighted in [18]) due to its parallel update structure.

AMP, VAMP and SwAMP have been extended to generalized but component-wise measurement models [19, 20, 14]. They have been then successfully applied to the phase recovery task where $\theta_n \sim \mathcal{U}[0, 2\pi]$, $\forall n \in \{1, \dots, N\}$, giving raise to the so-called “prGAMP” [21], “prVAMP” [10] and “prSAMP” [9] algorithms. In particular, the latter was shown to outperform other state-of-the-art algorithms among which the mean-field approximation [8].
<caption>**Algorithm 1 paSAMP Algorithm**</caption>

**Input:** y, D, σ<fnref n="2" />, ρ, σ<fnref n="2" />_x, μ_θ, Σ_θ, T_max
**Define:**
$$g_{out,n} \triangleq \frac{1}{\Sigma_{z_n}} (E_{Z|Y,P}\{z_n|y_n, \mu_{z_n}, \Sigma_{z_n}\} - \mu_{z_n})$$
$$g'_{out,n} \triangleq \frac{1}{\Sigma_{z_n}} \left( \frac{\mathrm{var}_{Z|Y,P}\{z_n|y_n, \mu_{z_n}, \Sigma_{z_n}\}}{\Sigma_{z_n}} - 1 \right)$$
$$g_{in,m} \triangleq E_{X|Y}\{x_m|\mu_{x_m}, \Sigma_{x_m}\}$$
$$g'_{in,m} \triangleq \mathrm{var}_{X|Y}\{x_m|\mu_{x_m}, \Sigma_{x_m}\}$$
1: **while** t < T_max **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;**for** n = 1...N **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\hat{z}_n(t) = \sum_{m=1}^{M} d_{nm} a_m(t)$$
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\Sigma_{z_n}^1(t+1) = \sum_{m=1}^{M} |d_{nm}|^2 v_m(t)$$
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\mu_{z_n}^1(t+1) = \hat{z}_n(t) - \Sigma_{z_n}^1(t) g_{out,n}$$
6: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
7: &nbsp;&nbsp;&nbsp;&nbsp;**for** m = permute[1...M] **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\Sigma_{x_m}(t+1) = \left(-\sum_{n=1}^{N} |d_{nm}|^2 g'_{out,n}\right)^{-1}$$
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\mu_{x_m}(t+1) = a_m(t) + \Sigma_{x_m}(t+1) \sum_{n=1}^{N} d_{nm} g_{out,n}$$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$v_m(t+1) = \Sigma_{x_m}(t+1) g'_{in,m}$$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$a_m(t+1) = g_{in,m}$$
12: &nbsp;&nbsp;&nbsp;&nbsp;**for** n = 1...N **do**
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\Sigma_{z_n}^{m+1}(t+1) = \Sigma_{z_n}^{m}(t+1) + |d_{nm}|^2 (v_m(t+1) - v_m(t))$$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\mu_{z_n}^{m+1}(t+1) = \mu_{z_n}^{m}(t+1) + d_{nm}(a_m(t+1) - a_m(t)) - g_{out,n}(t) (\Sigma_{z_n}^{m+1}(t+1) - \Sigma_{z_n}^{m}(t+1))$$
15: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
16: &nbsp;&nbsp;&nbsp;&nbsp;**update** σ<fnref n="2" /> according to [12]
17: &nbsp;&nbsp;&nbsp;&nbsp;**update** [$\theta_{m_n}$, $\Sigma_{\theta_n}$] according to (14-15)
18: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
19: **end while**
20: **Output:** {$\hat{x}_m = a_m(T_{max})$}_m

The prSAMP algorithm constitutes thus a promising approach for our DOA estimation<fnref n="2" /> problem (5). However, here, the phases $\theta_n$'s are spatially-correlated (as represented in the Markov model). This prevents us from a direct application of prSAMP.

We thus propose an iterative algorithm based on the two following mathematical derivations:

i) the extension of prSAMP to a i.i.d. Gaussian prior on the phases,

ii) the use of a mean-field approximation to estimate the (Gaussian) posterior distribution on the phases.

We detail both aspects in the next two sub-sections. In the following, we refer to the proposed procedure as “paSAMP” for “phase-aware SwAMP algorithm”. The pseudo-code of paSAMP is presented in **Algorithm 1**.

<footnote n="2">² Note in addition that the DOA estimation problem involves a highly-correlated matrix. This further motivates a SwAMP-like approach.</footnote>
## 3.2 Extension of prSAMP to i.i.d. Gaussian phases

AMP algorithms are based on the propagation of two types of messages: the “outgoing” messages and the “ingoing” messages from and to variables’ nodes {$x_m$}$_{m}$$_{={1...M}}$ . These messages are derived here for the prior distributions attached to the considered problem, namely (2) and (3)-(4)<fnref n="3" />. We first focus on the “outgoing messages”.

Considering $z_n \triangleq \sum_{m=1}^{M} d_{nm}x_m$, $\forall n \in \{1, ..., N\}$, we assume that the $z_n$'s follow Gaussian distributions with means $\mu_{z_n}$ and variances $\Sigma_{z_n}$ as linear combinations of $x_m$'s following Bernoulli-Gaussian distributions. By integrating over $\theta_n$ and resorting<fnref n="4" /> to an identification with a Von Mises distribution [22], we can write the moments of the posterior distribution as

$$E_{Z|Y}\{z_n|y_n, \mu_{z_n}, \Sigma_{z_n}\} = \frac{\Sigma_{z_n}}{\sigma^2 + \Sigma_{z_n}} \mathbf{R}_0\left(\frac{1}{\Sigma_{\theta}^z}\right) y_n e^{-j\mu_{\theta} z_n} + \frac{\sigma^2}{\sigma^2 + \Sigma_{z_n}} \mu_{z_n}, \quad (6)$$

$$\mathrm{var}_{Z|Y}\{z_n|y_n, \mu_{z_n}, \Sigma_{z_n}\} = \frac{|\Sigma_{z_n} y_n e^{-j\mu_{\theta} z_n} + \mu_{z_n} \sigma^2|^2}{|\sigma^2 + \Sigma_{z_n}|^2} \mathbf{R}_0\left(\frac{1}{\Sigma_{\theta}^z}\right) + \frac{\Sigma_{z_n} \sigma^2}{\sigma^2 + \Sigma_{z_n}} - E_{Z|Y}\{z_n|y_n, \mu_{z_n}, \Sigma_{z_n}\}^2, \quad (7)$$

with

$$\frac{1}{\Sigma_{\theta}^{z}} = \frac{1}{\alpha} + \frac{1}{\Sigma_{\theta_n}}, \quad \mu_{\theta_n}^{z} = \frac{-\frac{\arg(y_n^* \mu_{z_n})}{\alpha} + \frac{\mu_{\theta_n}}{\Sigma_{\theta_n}}}{\frac{1}{\alpha} + \frac{1}{\Sigma_{\theta_n}}}, \quad \alpha = \frac{\Sigma_{z_n} + \sigma^2}{|y_n| |\mu_{z_n}|},$$

$\mu_{\theta_n}$ (resp. $\Sigma_{\theta_n}$) is the marginalized posterior mean (resp. variance) of the phase noise $\theta_n$ as discussed in the next section, and $\mathbf{R}_0(\cdot) = \frac{I_1(\cdot)}{I_0(\cdot)}$ where $I_n(\cdot)$ is the modified Bessel function of the first kind at order $n$. We refer the reader to our technical report [23] which details the derivations of the computations.

Regarding the “ingoing” messages, which carry the prior information on the {$x_m$}$_{m}$$_{={1...M}}$ , the Bernoulli-Gaussian case has already been considered within the AMP context, in particular in [15]. Similarly to the “outgoing” messages, the moments of the “ingoing” messages resort to intermediary parameters $\mu_{x_m}$ and $\Sigma_{x_m}$ resp. homogeneous to the mean and variance of a Gaussian distribution:

$$E_{X|Y}(x_m|\mu_{x_m}, \Sigma_{x_m}) = \frac{\rho\sqrt{2\pi\nu^2}}{Z_{nor}} e^{-\frac{|m_x - \mu_{x_m}|^2}{2(\sigma^2 + \Sigma_{x_m})}} \gamma, \quad (8)$$

$$\mathrm{var}_{X|Y}(x_m|\mu_{x_m}, \Sigma_{x_m}) = \frac{\rho\sqrt{2\pi\nu^2}}{Z_{nor}} e^{-\frac{|m_x - \mu_{x_m}|^2}{2(\sigma^2 + \Sigma_{x_m})}} |\gamma^2 + \nu^2| - E_{X|Y}(x_m|\mu_{x_m}, \Sigma_{x_m})^2 \quad (9)$$

with

$$Z_{nor} = \rho\sqrt{2\pi\nu^2}e^{-\frac{|m_x - \mu_{x_m}|^2}{2(\sigma^2 + \Sigma_{x_m})}} + (1 - \rho)e^{-\frac{|\mu_{x_m}|^2}{2\Sigma_{x_m}}}, \quad (10)$$

$$\gamma = \frac{\sigma^2 \mu_{x_m} + \Sigma_{x_m} m_x}{\Sigma_{x_m} + \sigma^2}, \quad \nu^2 = \frac{\sigma^2 \Sigma_{x_m}}{\Sigma_{x_m} + \sigma^2}. \quad (11)$$

<footnote n="3">³ We refer the reader to papers [9, 14] for a more general presentation of the approach.</footnote>

<footnote n="4">⁴ We justify and develop this point in the technical report [23].</footnote>
We implement those calculations to paSAMP as new definitions of the two func-
tions $g_{in}$ and $g'_{in}$ defined in the pseudo-code **Algorithm 1**. We remind the reader
that, as an extended implementation of the SwAMP algorithm, the paSAMP al-
gorithm will conserve the structure described in [14] and [9]. For sake of clarity,
we use the notations introduced in [21] except for the scalar $d_{nm}$.

3.3 Mean-field approximation for the phase noise

The above expressions call on the knowledge of the moments of the posterior
distribution on θ. To simplify the latter computation, we propose in this step to
resort to a mean-field approximation. Following a similar reasoning as in [12],
we get

$$
q(\theta) = \mathcal{N}(\theta; \mu_{\theta}, \Sigma_{\theta}), \quad (12)
$$

$$
\text{where} \quad \Sigma_{\theta}^{-1} = \Lambda_{\theta}^{-1} + \operatorname{diag}\left(\frac{2}{\sigma^2} |\eta|\right), \qquad (13)
$$

$$
\mu_{\theta} = \Sigma_{\theta} \left( \operatorname{diag}\left( \frac{2}{\sigma^2} |\eta| \right) \operatorname{arg}(\eta) \right), \quad (14)
$$

with η_n = y_n Σ_{m=1}^{M} d_{nm}^* E_{X|Y}^* \{x_m | \mu_{x_m}, \Sigma_{x_m}\}, the n-th element in η (|η| stands here for the element-wise absolute value of η and * for the complex conjugate), and Λ_{\theta}^{-1} is the precision matrix attached to the prior distribution (4) on θ, i.e.

$$
\Lambda_{\theta}^{-1} =
\begin{pmatrix}
\frac{1}{\sigma_1^2} + \frac{\beta^2}{\sigma_{\theta}^2} - \frac{\beta}{\sigma_{\theta}^2} & 0 & 0 \\
-\frac{\beta}{\sigma_{\theta}^2} & \frac{1+\beta^2}{\sigma_{\theta}^2} & \cdots & 0 \\
0 & \cdots & \cdots & -\frac{\beta}{\sigma_{\theta}^2} \\
0 & 0 & -\frac{\beta}{\sigma_{\theta}^2} & \frac{1}{\sigma_{\theta}^2}
\end{pmatrix}.
\quad (15)
$$

Note that since the distribution $q(\theta)$ is Gaussian, marginals $q(\theta_n)$ used in
the previous “prSAMP-step” of the algorithm come as

$$
q(\theta_n) = \mathcal{N}(\theta_n; \mu_{\theta_n}, \Sigma_{\theta_n}) \quad (16)
$$

where $\mu_{\theta_n}$ (resp. $\Sigma_{\theta_n}$) is the $n$th element in $\mu_{\theta}$ (resp. in the diagonal of $\Sigma_{\theta_n}$).

Finally, following [12], we insert an estimation of the variance $\sigma^2$ of the ad-
ditive noise as a maximization step of an Expectation-Maximization (EM) algo-
rithm. Due to space limitation, we omit here the derivation of the computation,
but we refer again the reader to our technical report [23].

4 Numerical Experiments

In this section, we perform a quantitative and qualitative evaluation of the pro-
posed approach with respect to state-of-the-art algorithms.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Evolution of the (averaged) normalized correlation as a function of the variance $\sigma^2$ for $K = 2$ (left) and $K = 5$ (right), Comparison of the performance of conventionnal (delay-and-sum) beamforming (triangle mark), “prSAMP” (diamond mark), “paVBEM” (circle mark) and “paSAMP” (square mark). Experiments show that “paSAMP” provides better results and successfully integrates the phase noisy observation model.</figcaption></figure>

We consider the problem of the identification of the directions of arrival of $K$ plane waves from an antenna composed of $N = 256$ sensors. We assume that the angles of the $K$ incident waves can be written as $\phi_k = \frac{\pi}{2} + i_k \frac{\pi}{50}$ with $i_k \in [1, 50]$. A set of $M = 50$ steering vectors is defined from a set of angles $\{\phi_i = -\pi + i \frac{\pi}{50}\}_{i \in \{1,...,50\}}$ and the choice of the parameter $\Delta/\lambda = 4$. For each of the $K$ incident waves, the coefficient $x_{i_k}$ is initialized with $m_x = 0.5 + j0.5$, $\rho = K/M$ and $\sigma_x^2 = 0.1$. Finally, we set the following parameters for the phase Markov model (3): $\sigma_0^2 = 10$, $\sigma_\theta^2 = 0.1$ and $\beta = 0.8$. This corresponds to a situation where we have a high uncertainty on the initial value but a physical link between two space-consecutive angle measurements is taken into account.

We compare the performance of the following 4 different algorithms: i) the standard beamforming introduced in [1] (dashed black curve, triangle mark); ii) the prSAMP algorithm proposed in [9] as a solution to the phase retrieval problem (continuous black curve, diamond mark); iii) the paVBEM procedure proposed in [12] exploiting the same prior models (dashed red curve, circle mark); iv) the paSAMP algorithm described in Section 3 (continuous blue curve, square mark). To evaluate the performance of these procedures, we consider the normalized correlation between the ground truth $\mathbf{x}$ and its reconstruction $\hat{\mathbf{x}}$, that is $\frac{|\mathbf{x}^H\hat{\mathbf{x}}|}{\|\mathbf{x}\|\|\hat{\mathbf{x}}\|}$, as a function of the additive noise variance $\sigma^2$. This quantity is averaged over 100 realizations for each point of simulation.

The results achieved by the 4 procedures are presented in Figure 1, resp. for $K = 2$ (left) and $K = 5$ (right) sources. In both cases, we see that the conventional beamforming and the prSAMP algorithm fail to reconstruct $\mathbf{x}$ properly. These results illustrate the benefits of carefully accounting for the phase noise in fluctuating environments. We can also notice the superiority of paSAMP over its mean-field counterpart paVBEM, especially in presence of a strong additive
noise. This comes in the continuity of previous work [9], where prSAMP proved to outperform prVBEM. Finally, it is interesting to compare the performance of both paSAMP and paVBEM algorithms with regard to the number of sources. Both achieve better performance when confronting to $K = 5$ sources than to $K = 2$ sources. As mentioned in [12], this behavior is typical for the phase retrieval problems, where the loss information on the phases can be compensated by a larger number of sources. In addition, we observe that the performance gap between paSAMP and paVBEM tends to increase with the number of sources. This is in accordance with previous work [13] demonstrating the relevance of the Bethe approximation over the mean-field approximation when the signal to recover exhibits a low sparsity (i.e. a high number of non-zero coefficients).

# 5 CONCLUSION

We have presented here a novel AMP algorithm able to perform DOA estimation in a corrupted phase-noisy environment. This approach exploits both a sparsity prior on the sources and a structured prior on the phase noise. Compared to state-of-the-art algorithms, the approach presents a good behaviour illustrating a successful inclusion of the different assumptions. In particular, it outperforms a recent algorithm dealing with the same DOA estimation problem in fluctuating environments. Future work will include further assessment on real data.

# 6 ACKNOWLEDGMENT

The authors thank Boshra Rajaei for sharing her MATLAB implementation of the prSAMP algorithm.

# References

1. Johnson, D.H., Dudgeon, D.E.: Array Signal Processing: Concepts and Techniques. Prentice Hall (1993).

2. Schmidt, R.: Multiple emitter location and signal parameter estimation. IEEE transactions on antennas and propagation. vol. 34, pp. 276–280 (1986).

3. Xenaki, A., Gerstoft, P., Mosegaard, K.: Compressive beamforming. The Journal of the Acoustical Society of America. vol. 136, pp. 260–271 (2014).

4. Dashen, R., Flatté, S. M., Munk, W. H., Watson, K. M., Zachariasen, F.: Sound transmission through a fluctuating ocean. Cambridge University Press (2010).

5. Colosi, J.A.: Sound Propagation through the Stochastic Ocean. Cambridge University Press (2016).

6. Cheinet, S., Ehrhardt, L., Juve, D., Blanc-Benon, P.: Unified modeling of turbulence effects on sound propagation. The Journal of the Acoustical Society of America. vol. 132, pp. 2198–2209 (2012).

7. Ehrhardt, L., Cheinet, S., Juve, D., Blanc-Benon, P.: Evaluating a linearized Euler equations model for strong turbulence effects on sound propagation. The Journal of the Acoustical Society of America. vol. 133, pp. 1922–1933 (2013).
8. Dremeau, A., Krzakala F.: Phase recovery from a Bayesian point of view: the variational approach. Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3661–3665 (2015).

9. Rajaei, B., Gigan, S., Krzakala, F., Daudet, L.: Robust phase retrieval with the swept approximate message passing (prSAMP) algorithm. IPOL. vol. 7, pp. 43–55 (2016).

10. Metzler, C. A., Sharma, M. K., Nagesh, S., Baraniuk, R. G., Cossairt, O. and Veeraraghavan, A. : Coherent Inverse Scattering via Transmission Matrices: Efficient Phase Retrieval Algorithms and a Public Dataset, Proc. IEEE International Conference on Computational Photography, pp. 1–16 (2017).

11. Waldspurger, I., d'Aspremont, A., Mallat, S.: Phase recovery, maxcut and complex semidefinite programming. Mathematical Programming, vol. 149. Springer (2015).

12. Dremeau, A. and Herzet, C.: DOA estimation in structured phase noisy environments: technical report (2016).

13. Krzakala, F., Manoel, A., Tramel, E.W., Zdeborova, L.: Variational free energies for compressed sensing. 2014 Proc. IEEE International Symposium on Information Theory (ISIT), pp. 1499–1503 (2014).

14. Manoel, A., Krzakala, F., Tramel, E., Zdeborova, L.: Swept approximate message passing for sparse estimation. Proc. of the 32nd International Conference on Machine Learning (ICML-15). pp. 1123–1132 (2015).

15. Vila, J., Schniter, P.: Expectation-maximization Bernoulli-Gaussian approximate message passing. Proc. Signals, Systems and Computers (ASILOMAR), pp. 799–803 (2011).

16. Soussen, C., Idier, J., Brie, D., Duan, J.: From Bernoulli-Gaussian deconvolution to sparse signal restoration. IEEE Transactions on Signal Processing (2011).

17. Rangan, S., Schniter, P., Fletcher, A. K. : Vector Approximate Message Passing. Proc. IEEE International Symposium on Information Theory (ISIT), pp. 1588–1592 (2017).

18. Caltagirone, F., Zdeborova, L., Krzakala, F.: On Convergence of Approximate Message Passing. Proc. IEEE International Symposium on Information Theory (ISIT), pp. 1812–1816 (2014).

19. Rangan, S.: Generalized approximate message passing for estimation with random linear mixing. Proc. IEEE International Symposium on Information Theory (ISIT), pp. 2168–2172 (2011).

20. Schniter, P., Rangan, S., Fletcher, A. K : Vector approximate message passing for the generalized linear model. Asilomar Conference on Signals, Systems and Computers, pp. 1525–1529 (2016).

21. Schniter, P., Rangan, S: Compressive phase retrieval via generalized approximate message passing. IEEE Transactions on Signal Processing, vol. 63 pp. 1043–1055 (2015).

22. Mardia, K. V., Jupp, P. E.: Directional statistics. John Wiley & Sons (2009).

23. Beaumont, G., Fablet, R., Dremeau, A.: DOA estimation in fluctuating environments: an approximate message-passing approach. Technical Report (2018). https://hal.archives-ouvertes.fr/hal-01624855v4/document