<figure><img src="image_21.png" /><figcaption>Fig. 1. Illustration of a f-k plane obtained within configuration setup exposed in Section IV-A using 240 sensors.</figcaption></figure>

applications in underwater acoustics [22], [23], and for imaging the modal dispersion of surface waves [24]. The first application of CS to estimate modal spectrum in underwater acoustics has been proposed in [25].

When considering a broadband source, wavenumber spectra can be estimated at each of the source frequencies. The concatenation of them results in a frequency-wavenumber ($f-k$) diagram (see Fig. 1), representative of the waveguide dispersion [12], [13], [25]. As a consequence, the $f-k$ diagram constructions suffer from the drawbacks inherited from each wavenumber spectrum estimation. To partially solve this problem, Le Courtois and Bonnel have proposed [26] a postprocessing method to track the wavenumbers in badly resolved $f-k$ diagrams. Tracking is performed using particle filtering (PF) [27]. To do so, the wavenumber spectrum (at a given frequency) is modeled as a dynamic system parametrized by two equations: a system equation and an observation equation. The observation equation is the equation that allows generating a wavenumber spectrum if the (discrete) wavenumbers are known. On the other hand, the system equation is an iterative relation linking the wavenumbers from one frequency to another. This system equation thus benefits from a physical hypothesis on the wavenumber dispersion in the frequency dimension, and such an idea will also be used in our paper. Note that the PF method is particularly interesting when the number of sensors is too small to separate the modes. However, the number of propagating modes must be known a priori to initiate the tracking operation. In addition, the resolution of the method depends on the distribution of the particles and cannot be easily estimated. These last two constraints will be eased in this paper.

This paper introduces a framework relying on the hypothesis that the propagation is sparse (in the wavenumber dimension) and dispersive (in the frequency dimension). A methodology allowing $f-k$ estimation and integrating these two physical hypotheses is presented. To this end, the work focuses on the sparse estimation of the wavenumber spectrum using a Bernoulli-Gaussian model [28]. In this framework, each wavenumber spectrum is a sparse vector. Each wavenumber bin is supposed to be “on” (i.e., there is a mode propagating at this particular

wavenumber value) or “off” (i.e., there is no mode propagating at this particular wavenumber value) according to a Bernoulli distribution. Once a wavenumber bin is “on,” the corresponding mode amplitude is supposed to follow a Gaussian distribution. In this approach, the Bernoulli parameter offers a great opportunity to take into account the dispersive propagation, so that the sparse wavenumber spectra estimated at each frequency are naturally related with each other. In this paper, the Soft Bayesian Pursuit (SoBaP) algorithm [29] is considered to efficiently perform the estimation procedure.

The paper is organized as follows. Section II recalls the physical principles of the modal propagation and the prior knowledge exploited to build up the $f-k$ representation. Section III is dedicated to a mathematical formulation of the problem and the integration of the physical prior. Finally, in Section IV, the method is applied to synthetic simulations and to marine data recorded in the North Sea.

## II. ACOUSTIC PROPAGATION IN DISPERSIVE SHALLOW-WATER ENVIRONMENTS

### A. Received Signal

In shallow-water environments, the acoustic propagation is described by modal theory. The propagation depends on the frequency $f$ and is dispersive. When considering an emitting source $s(f)$ at depth $z_s$, the received signal on a sensor located at distance $r$ and depth $z$ can be written as [30]

$$y(f, r) = Q \frac{s(f)}{\sqrt{r}} \sum_{m=1}^{M(f)} \psi_m(f, z_s) \psi_m(f, z) \frac{e^{-jrk_{rm}(f)}}{\sqrt{k_{rm}(f)}} + w(f, r) \quad (1)$$

where $Q$ is a constant factor, $M(f)$ is the number of propagating modes at frequency $f$, $k_{rm}(f)$ is the horizontal wavenumber of the $m$th mode, and $\psi_m(f, z)$ is the modal depth function of the $m$th mode. The quantity $w(f, r)$ stands for the TFT of the noise attached to the measurements. Note that, for the sake of clarity, we deliberately omit here the dependence in $z$ in the notation of $y(f, r)$ and $w(f, r)$ as we will consider HLAs, i.e., sensors at a constant depth.

For subsequent convenience, we define the amplitude of the $m$th mode as

$$A_m(f) \triangleq \frac{\psi_m(f, z_s) \psi_m(f, z)}{\sqrt{k_{rm}(f)}} \quad (2)$$

so that (1) can finally be reexpressed as

$$y(f, r) = Q \frac{s(f)}{\sqrt{r}} \sum_{m=1}^{M(f)} A_m(f) e^{-jrk_{rm}(f)} + w(f, r). \quad (3)$$

Considering an HLA of regularly spaced sensors aligned with the source (endfire position), the $l$th sensor is at a distance of $r = (l - 1)\Delta_r + r_0$ from the source, where $r_0$ is the distance of the first sensor and $\Delta_r$ is sensor spacing. In this paper, we assume a distant source (so that $r_0 \gg \Delta_r$). Under this assumption, the geometrical attenuation factor $1/\sqrt{r}$ can be considered as constant over the entire sensor array. Similarly, we will neglect
any eventual modal attenuation, so that the modal wavenumber $k_{rm}(f)$ may be considered as real numbers.

As stated in the introduction, considering a monochromatic source, the wavenumber spectrum can be estimated by performing a simple SFT along the HLA. If a broadband source is available, an $f-k$ diagram can be obtained by concatenating the wavenumber spectra at several frequencies [12]. However, this requires a long and dense HLA. The goal of this paper is to estimate $f-k$ diagram in less constrained contexts, in particular, when the number of sensors is low.

## B. Dispersion Relation

In waveguides, the horizontal wavenumbers $k_{rm}$ are linked to their vertical counterparts $k_{zm}$ by the dispersion relation, that is, for a given frequency $f$

$$ \left(\frac{2\pi f}{c}\right)^2 = k_{rm}(f)^2 + k_{zm}(f)^2 \quad (4) $$

where $c$ is the speed of sound.<fnref n="1" />

Discretizing the frequency axis (with $f = \nu\Delta_f$, $\nu \in \mathbb{N}$) and denoting $k_{rm}[\nu] = k_{rm}(\nu\Delta_f)$, the wavenumbers attached to two successive indices are then defined as [25]

$$ k_{rm}[\nu + 1]^2 = k_{rm}[\nu]^2 + (2\nu + 1) \left(\frac{2\pi\Delta_f}{c}\right)^2 + \epsilon[\nu] \quad (5) $$

where $\Delta_f$ is the sampling period, i.e., the step between two frequency bins, and $\epsilon[\nu] = k_{zm}[\nu]^2 - k_{zm}[\nu + 1]^2$. In shallow-water environments, the vertical wavenumbers $k_{zm}$ weakly depend on the frequency [30]; the quantity $\epsilon$ is smaller than the other terms of the equation and can be neglected. This approximation may not be true when the propagation becomes strongly range dependent. This case will not be considered in the remainder of the paper.

Equation (5) expresses the dispersive nature of the wavenumbers in a shallow-water environment. In this paper, we propose to take into account this physical information in the reconstruction procedure of the $f-k$ diagram. Note that (5) was used in [26] and [25] to postprocess badly resolved $f-k$ diagrams. In this paper, (5) will be embedded as a prior information into a Bayesian CS algorithm, allowing a direct estimation of the $f-k$ diagram.

## III. COMPRESSED SENSING

When dealing with compressed sensing, two properties have to be verified:

1) sparsity: the signal to be acquired can be represented with a few nonzero elements in a given representation basis;

2) incoherence: the sensing must be made in a domain “as orthogonal as possible” to the representation basis.

Under these conditions, the theory of compressed sensing guarantees that we can recover the signal from a number of

samples of the order of the number of its nonzero elements in the representation basis.

In other words, “incoherence” expresses the fact that a signal admitting a sparse representation in a particular domain (for example, in the frequency domain) must be acquired in the domain where its representation is spread out (for example, the time domain). In the time-frequency example, the sensing domain is the canonical basis (time) and the sparsity domain is the Fourier basis (frequency).

The reconstruction problem of the $f-k$ diagram is then a perfect application case for compressed sensing: the signal of interest is sparse in the wavenumber domain (only few wavenumbers propagate in shallow-water environments and at low frequencies), and is acquired in the spatial domain. From a mathematical point of view, this means that we acquire the signal in the canonical basis (space, through the sensor array) and the sparsity domain is the spatial Fourier basis (wavenumber). We refer the reader to [31] for a brief introduction to compressed sensing, or to [32] for more insights.

### A. State of the Art

Formally, let $y_{\nu} \in \mathbb{C}^L$ be the signal measured over the $L$ sensors at the frequency index $\nu$. Adopting a discretized matrix formulation, (3) can be reexpressed as

$$ y_{\nu} = Dz_{\nu} + w_{\nu} \quad (6) $$

where $D$ is an $(L \times N)$-dictionary of Fourier discrete atoms, namely whose $(\ell, n)$-element is $d_{nl} \triangleq e^{-j2\pi \frac{n\ell}{N}}$, and $z_{\nu} = [z_{\nu,1}, \dots, z_{\nu,N}]^T \in \mathbb{C}^N$ is the wavenumber spectrum at the frequency index, i.e., the $\nu$th transposed line of the $f-k$ diagram to estimate.

Then, $N$ corresponds to the number of discretized points in the horizontal wavenumber domain, say

$$ \kappa_{rn} \triangleq 2\pi \frac{n}{N\Delta_r} \quad \forall n \in \{1, \dots, N\} \quad (7) $$

and $z_{\nu}$ gathers the amplitudes attached to each of the wavenumber bins.

According to the modal theory in shallow-water environments and at low frequencies, we have $M[\nu] \ll N$ [see (3)], where $M[\nu] = M(\nu\Delta_f)$. In other words, the vector $z_{\nu}$ has few nonzero elements, corresponding to the propagating modal wavenumbers. Formally, we define by $M_{\nu} \subset \{1, \dots, N\}$ the set of cardinality $M[\nu]$ gathering the indices of the propagating wavenumbers. For each of the $M[\nu]$ propagating modes, there exists $n \in M_{\nu}$ such as,

$$ k_{rm}[\nu] = \kappa_{rn} \quad \text{and} \quad A_m[\nu] = z_{\nu,n} \quad (8) $$

where $A_m[\nu] = A_m(\nu\Delta_f)$ defined as in (2).

We note here that the precise estimation of the set $M_{\nu}$ is crucial. From it derives the knowledge of the propagating wavenumbers $k_{rm}[\nu]$ and the modal amplitudes $A_m[\nu]$, of particular interest for source depth estimation and/or environmental inversion.

The sparsity of the vector $z_{\nu}$ constitutes important information on the $f-k$ diagram, which should be taken into account

<footnote n="1">¹Note that quantities $c$ and $k_{zm}(f)$ are depth dependent. Here again, we choose to skip this dependence in the expression as it will be considered as constant in our sensing context. As we will see in the following, the way the dispersion relation will be exploited in the proposed approach allows some deviations from this assumption.</footnote>
in the reconstruction procedure. Several formulations of the corresponding sparse recovery problem can then be considered. In this paper, we focus on the following one:

$$ \hat{\mathbf{z}}_{\nu} = \underset{\mathbf{z}_{\nu}}{\operatorname{argmin}} \|\mathbf{y}_{\nu} - \mathbf{D}\mathbf{z}_{\nu}\|_{2}^{2} + \lambda \|\mathbf{z}_{\nu}\|_{0} \quad (9) $$

where $\|\mathbf{z}_{\nu}\|_{0}$ stands for the $l_{0}$-pseudonorm of $\mathbf{z}_{\nu}$ (counting the number of nonzero elements in $\mathbf{z}_{\nu}$) and $\lambda$ is a parameter specifying the tradeoff between the sparsity constraint and the data-fidelity term $\|\mathbf{y}_{\nu} - \mathbf{D}\mathbf{z}_{\nu}\|_{2}^{2}$. The solution of (9) can be naturally interpreted as the least square (LS) solution (which, in our case, is equivalent to an SFT when $\mathbf{y}_{\nu}$ is uniformly sampled) penalized by the sparsity of the solution.

Solving (9) is an NP-hard problem [33], i.e., it generally requires a combinatorial search over the entire solution space. Therefore, heuristic (but tractable) algorithms have been devised to deal with this problem. We can roughly divide them into three families: 1) the greedy algorithms (e.g., [34] and [35]) which build up the sparse solution by making a succession of greedy decisions; 2) the relaxation-based algorithms (e.g., [36] and [37]) which replace the $l_{0}$-pseudonorm by some $l_{p}$-norm (with $p \in \{0, 1\}$) leading to a relaxed problem efficiently solvable by standard optimization procedures; and 3) the Bayesian algorithms (e.g., [38], [29], and [39]) which express the problem as the solution of a Bayesian inference problem and apply statistical tools to solve it.

The Bayesian approaches are particularly suitable when the noise level is not known. In addition, they offer a simple framework to implement any prior information available on the physical context of the problem. In our case, they present in particular a great opportunity to take into account the dispersion relation given in (5).

## B. Bayesian Formulation

Let us first assume $\mathbf{w}_{\nu}$ in (6) to be a circular Gaussian noise (denoted by $\mathcal{CN}$) with zero mean and variance $\sigma_{w}^{2}$. We suppose then that $\mathbf{z}_{\nu}$ is the realization of a Bernoulli-Gaussian (BG) model [29], that is

$$ \mathbf{z}_{\nu} = \mathbf{s}_{\nu} \odot \mathbf{x}_{\nu} \quad (10) $$

where $\odot$ represents here the term-by-term product, and

$$ p(\mathbf{s}_{\nu}) = \prod_{n=1}^{N} p(s_{\nu,n}), \quad \text{with} \quad p(s_{\nu,n}) = \text{Ber}(p_{\nu,n}) \quad (11) $$

$$ p(\mathbf{x}_{\nu}) = \prod_{n=1}^{N} p(x_{\nu,n}), \quad \text{with} \quad p(x_{\nu,n}) = \mathcal{CN}(0, \sigma_{x}^{2}). \quad (12) $$

The Bernoulli distribution (denoted by Ber in the equation above) has a realization domain on $\{0, 1\}$ and depends on a parameter $p_{\nu,n}$ which represents the probability of being equal to 1. The circular Gaussian distribution put on the variable $\mathbf{x}_{\nu}$ is assumed to be with zero mean and variance $\sigma_{x}^{2}$.

In other words, $\mathbf{s}_{\nu}$ is called the support of the sparse representation. It is such as $\mathbf{s}_{\nu} \in \{0, 1\}^{N}$. Every 1-value in $\mathbf{s}_{\nu}$ indicates a wavenumber bin corresponding to a propagating

modal wavenumber; every 0-value in $\mathbf{s}_{\nu}$ indicates a wavenumber bin corresponding to a wavenumber that does not propagate. The quantity $\mathbf{x}_{\nu}$ stands for the amplitude of sparse representation. Coming back to (8), $\mathbf{s}_{\nu}$ serves then to define the set $\mathcal{M}_{\nu}$ of the wavenumber bins corresponding to the propagating wavenumbers

$$ \mathcal{M}_{\nu} = \{n \in \{1, \dots, N\} | s_{\nu,n} = 1\} $$

while $\mathbf{x}_{\nu}$ can be directly linked to the amplitudes $A_{m}[\nu]$ of the propagating modes. As a consequence, the latter are assumed to be independently and identically distributed according to a Gaussian law, as expressed in (12).

Formally, the BG model (11)-(12) is well suited to modeling situations where $\mathbf{y}_{\nu}$ stems from a sparse process: if $p_{\nu,n} \ll 1$, $\forall n$, only a small number of $s_{\nu,n}$'s will typically be nonzero, i.e., the observations $\mathbf{y}_{\nu}$ will be generated with high probability from a small subset of the columns of $\mathbf{D}$. More rigorously, it has been proved (see [28] and [40]) that, within model (11)-(12), the joint maximum a posteriori estimation problem shares the same set of solutions as the standard sparse recovery problem (9). This connection gives weight to the general use of this model in sparse recovery problems.

In this paper, we choose to resort to a particular algorithm of the sparsity literature, that is the SoBaP algorithm [29]. Exploiting model (6)-(12), SoBaP aims at solving the following marginalized MAP problem:

$$ \hat{\mathbf{s}}_{\nu} = \underset{\mathbf{s}_{\nu}}{\operatorname{argmax}} p(\mathbf{s}_{\nu} | \mathbf{y}_{\nu}) \quad (13) $$

where

$$ p(\mathbf{s}_{\nu} | \mathbf{y}_{\nu}) = \int_{\mathbf{x}_{\nu}} p(\mathbf{s}_{\nu}, \mathbf{x}_{\nu} | \mathbf{y}_{\nu}) d\mathbf{x}_{\nu}. \quad (14) $$

The algorithm relies on a variational approximation [41] which constitutes the building block of the current state-of-the-art algorithms [42].

Once $\mathbf{s}_{\nu}$ is estimated, the amplitudes of the propagated wavenumbers can be computed by a simple pseudoinversion of the dictionary $\mathbf{D}$, restricted to its nonzero columns, say $\mathbf{D}_{\hat{\mathbf{s}}_{\nu}}$

$$ \hat{\mathbf{x}}_{\nu} = \mathbf{D}_{\hat{\mathbf{s}}_{\nu}}^{\dagger} \mathbf{y}_{\nu} \quad (15) $$

where $^{+}$ stands for the pseudoinverse operator.

## C. Incorporating A Priori Information About Dispersive Propagation

In our application context, the BG model presents an additional advantage: relying on the dedicated support variable $\mathbf{s}_{\nu}$, it allows for an easy implementation of a priori information through the Bernoulli parameters $p_{\nu,n}$ [see (11)].

By the physical relation (5), the wavenumbers estimated at the frequency $\nu$ define a priori values for the wavenumbers at the next frequency $\nu + 1$. Consequently, the sparse representation support can be propagated to $\nu + 1$. This can be done by means of the Bernoulli parameters $p_{\nu+1,n}$ which set the probabilities for the elements in $\mathbf{s}_{\nu+1}$ to be set to 1 (i.e., for the wavenumbers
<figure><img src="image_22.png" /><figcaption>Fig. 2. Illustration of the definition of the Bernoulli parameters as expressed in (16) and (17), on two successive rows of the f-k diagram.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 3. Scheme summarizing the overall proposed methodology for f-k reconstruction.</figcaption></figure>

to be chosen). Formally, we fix, for all $n \in \{1, \dots, N\}$

$$p_{\nu+1,n} = \begin{cases} 0.7, & \text{if } n \in \mathcal{I}_{\nu} \\ 0.3, & \text{if } (n-1) \in \mathcal{I}_{\nu} \\ 0.3, & \text{if } (n+1) \in \mathcal{I}_{\nu} \\ M[\nu+1]/N, & \text{otherwise} \end{cases} \quad (16)$$

with

$$\mathcal{I}_{\nu} = \left\{ n \mid n = \left[ \sqrt{n'^2 + (2\nu + 1) \left( \frac{N\Delta_r \Delta_f}{c} \right)^2} \right], \hat{s}_{\nu,n'} = 1 \right\} \quad (17)$$

where [.] is the nearest natural number. We explain these equations hereafter.

Expression (17) is directly derived from (5), where we have replaced the wavenumbers $k_{rm}[\nu+1]$ and $k_{rm}[\nu]$ by their discretized expressions (7), introducing indices $n$ and $n'$ [see (8)]. To make it clear, let us assume that at frequency $\nu$, a mode propagates at the wavenumber bin $n'$, i.e., $\hat{s}_{\nu,n'} = 1$. Then, (5) allows to predict the wavenumber value at frequency $\nu + 1$. The corresponding wavenumber bin is the index $n$ given by (17). As a result, a weight of 0.7 is given to the bin $n$ predicted at frequency $\nu + 1$ (i.e., such that $n \in \mathcal{I}_{\nu}$). Fig. 2 illustrates the procedure.

The round operator [.] in (17) throws off the unknown quantity $\epsilon[\nu]$ in (5) which contains the (weak) frequency dependence of the vertical wavenumbers. It is, however, taken into account in a probabilistic way through the Bernoulli parameters of the nearest neighbors of the indices selected by $\mathcal{I}_{\nu}$. Then, if a wavenumber is found to propagate at frequency $\nu$ and bin $n'$ (i.e., $\hat{s}_{\nu,n'} = 1$), a weight of 0.3 is given to the position $n+1$ and $n-1$ at frequency $\nu+1$ where $n \in \mathcal{I}_{\nu}$.

Note, finally, that the numerical values of 0.7 and 0.3 are arbitrarily chosen. The other positions, which have not been predicted by the dispersion relation are assumed to have the same probability to be chosen, equal to $M[\nu+1]/N$.

Interestingly, we note that, by virtue of the flexibility of the Bayesian framework, the value of $M[\nu+1]$ in (16) does not have to be precisely known, a rough guess (e.g., within an error of two modes) being sufficient to guarantee robust results; here, the knowledge of the number of propagating modes is less essential than in HR and tracking methods. As examples, strong a priori information is necessary for the initialization of particle filtering [25], for the choice of the eigenvalues in MUSIC [20] or the order of autoregressive models as in [16] and [18]. Also, the environmental parameters, such as the sound speed, may not be perfectly known without impacting the performance of the approach (see experiments with real data in Section IV-B). Note that, in cases of strong environmental uncertainties, a simple and straightforward way to ensure robustness will be to extend the neighborhood of significative probabilities [indices $n+1$ and $n-1$ in (16)] to farther wavenumber bins (for example $n+2$ and $n-2$). Doing so, we can compensate for the uncertainty over the propagation environment by a more relaxed probabilistic model. Note that the extreme case will be then used to set all parameters $p_{\nu+1,n}$ to the same value, resulting in the standard sparse case, as originally considered in [29].

Fig. 3 illustrates the entire f-k diagram estimation process. Given a frequency index $\nu$ of the f-k plane, wavenumbers are estimated by solving the problem (13) using SoBaP. Then (16) predicts the value of the Bernoulli parameter at the frequency index $\nu + 1$. The frequency is incremented to index $\nu + 1$ and problem (13) is solved using the updated Bernoulli parameters.

## IV. APPLICATIONS

In this section, we propose two experimental setups to assess the performance of the proposed method: to quantify the performance, we first consider synthetic experiments, before applying the method on real data.

### A. Simulations

The synthetic data are simulated using a Pekeris waveguide [30]. The water column is assumed to be $D = 130$ m deep with a sound speed $c_{water} = 1500$ m/s and a density $\rho_{water} = 1$ kg/m$^3$. The seabed is a semi-infinite fluid layer with a sound speed of $c_{seabed} = 2000$ m/s and a density $\rho_{seabed} = 2$ kg/m$^3$. The array is composed of 240 hydrophones lying on the seabed and with spacing $\Delta_r = 25$ m, resulting in a 6000-m-long antenna. The corresponding spectral resolution of the SFT is then $1.05 \times 10^{-4}$ rad/m. The source emits a broadband signal between 0 and
<figure><img src="image_24.png" /><figcaption>Fig. 4. Reconstruction of the dispersion curves (in decibels) using measurements from 30 sensors and an SNR of 10 dB: (a) LS inversion; (b) OMP; and (c) SoBaP with dispersive *a priori*. The ground truth is represented by red points.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 5. Normalized MSE on the estimated wavenumbers (in decibels) as a function of: (a) the SNR for an array of $L = 120$ sensors; and (b) the number of sensors in the array for an SNR of 10 dB.</figcaption></figure>

50 Hz (corresponding to a white noise or an impulsive signal in the time domain; the latter case was simulated here). The frequency resolution is $\Delta_f = 0.2$ Hz. The source is placed at $D = 130$-m depth to avoid nodes of the modal functions. At 50 Hz, five modes are propagating.

Two sparse representation algorithms are compared for the resolution of the problem: orthogonal matching pursuit (OMP) [43] as used in [24], and SoBaP incorporating *a priori* information on the dispersive propagation, as defined earlier. OMP is stopped when $||y_\nu - D\hat{x}_\nu||_2 < \sqrt{L\sigma_w^2}$. SoBaP is stopped when the Kullback-Leibler divergence between the distribution $p(s, x|y)$ and its variational approximation is smaller than $10^{-4}$. Both algorithms are compared to the least square (LS) approach, performing a simple inverse Fourier transform.

The resulting $f-k$ diagrams are shown in Fig. 4 together with the ground truth, represented by red points. For all three, 30 sensors are used and a signal-to-noise ratio (SNR) of 10 dB is considered. To account for the CS framework, the positions of the sensors are randomly selected from the 240 simulated measurements. At these settings, the $f-k$ diagram obtained by simple LS inversion and plotted in Fig. 4(a) cannot be correctly interpreted to infer the expected wavenumbers: they do not appear clearly. The one obtained by OMP seems to be more

sensitive to noise; some wavenumbers are identified but they have no physical justification. The use of the dispersive *a priori* information in SoBaP suppresses all of the artifacts from the $f-k$ plane. These visual results are to be compared with the ground truth, represented by red points on each $f-k$ diagram.

To quantify more precisely the performance of the algorithms, we compute for each of them the normalized mean square error (MSE) between the estimated wavenumbers (corresponding, for each frequency index $\nu$, to the $M[\nu]$ largest coefficients in $\hat{x}_\nu$) and the true values of the Pekeris model. The error is averaged over the mode number and the frequencies to obtain a global value. The MSE is represented, in Fig. 5(a), as a function of the SNR for a given number of sensors $L = 120$ and in Fig. 5(b), as a function of the number of hydrophones exploited to reconstruct the $f-k$ diagram for a given SNR of 10 dB.

In accordance with previous works [25], [24], both figures illustrate the good behavior of sparse-aware algorithms with regard to naïve LS inversion: OMP and the proposed SoBaP procedure outperform the LS estimation in most practical setups, i.e., noisy measurements and few sensors.

As shown in Fig. 5(a), for an SNR greater than 20 dB, the two sparse-aware methods present similar performances, and an exact reconstruction of the wavenumbers is obtained (note
<figure><img src="image_26.png" /><figcaption>Fig. 6. Reconstruction of the dispersion curves (in decibels) from the North Sea measurements using: (a) 240 sensors with an LS inversion; (b) 20 sensors with an LS inversion; (c) 20 sensors with the OMP algorithm; and (d) 20 sensors with the physics-aware SoBaP algorithm. Red points represent the result of a peak-picking method on the LS inversion on the entire antenna (a).</figcaption></figure>

that the residual MSE is linked to the resolution of the Fourier matrix D, independently of the number of measurements). The proposed approach proves to be more interesting in the noisy cases: we thus observe that the SoBaP procedure exploiting the frequency dependence between the wavenumbers shows a faster decrease of the MSE as the SNR increases. For an SNR of 10 dB, the method reaches an MSE around -35 dB, compared to -20 dB for OMP.

This particular setting is considered in more detail in Fig. 5(b), where the performance at SNR = 10 dB is assessed from a sensor point of view. Within this setting, the three methods achieve similar (poor) performance for very small numbers of sensors (below L = 50). Note that this observation does not contradict the visual interpretation of Fig. 4, as the MSE quantifies the position of the $M[ν]$ largest estimated wavenumbers without any consideration over false alarms. Beyond 50 sensors, the proposed SoBaP procedure outperforms indisputably the other two approaches.

## B. North Sea Data

In this section, we apply our approach to real data. The considered data were acquired in the North Sea during a seismic campaign led by the Compagnie Générale de Géophysique [44], [45]. The source is an air gun that is impulsive with a nearly flat spectrum between 0 and 80 Hz. Measurements are performed by a synthetic antenna of 240 ocean-bottom seismometers resting on the seabed. They are spaced at intervals of 25 m, leading to a total length of 6000 m. The pressure field is sampled at 250 Hz. The environment is assumed to be close to a Pekeris waveguide. Within this assumption, the experimental parameters were estimated [45] as: $c_{water} = 1520$ m/s, $c_{seabed} = 1875$ m/s, and $D = 130$ m. For these data, we evaluated the SNR around 13 dB.

Fig. 6 presents the $f-k$ representations obtained by the inversion of 20 hydrophones with: a simple LS method [Fig. 6(b)], the OMP algorithm [Fig. 6(c)], and the SoBaP algorithm exploiting the frequency dependence between the wavenumbers [Fig. 6(d)]. As a comparison, the $f-k$ diagram obtained with
an LS inversion using the entire antenna is displayed in Fig. 6(a). Moreover, in each figure, we have added in red the wavenumbers selected from the LS inversion on the entire antenna [Fig. 6(a)] using a peak-picking method. These results show the relevance of the proposed approach. We can thus see that sparsity enables an undeniable improvement of the inversion (compared to the simple LS estimate) in particular in terms of noise suppression, but missed (in particular at the highest frequencies) and false detections are still noticeable. Adding prior knowledge on the propagated wavenumbers allows us to recover perfectly the wavenumbers, even at the highest frequencies (around 60 Hz in our case), while reducing the noise of the *f* − *k* diagram.

V. CONCLUSION

This paper deals with the estimation of the wavenumbers
in shallow-water environments. An array processing method
is proposed to reduce the number of the sensors required to
separate accurately the wavenumbers over a wide range of fre-
quencies. The proposed method relies on two robust hypotheses.
On the one hand, in the wavenumber domain, the wavenumber
spectrum is sparse. On the other hand, in the frequency domain,
the wavenumber spectra can be related from one frequency to the
next using a general dispersion relationship. These two hypothe-
ses call for the utilization of a Bayesian compressed-sensing
(CS) algorithm. Indeed, the CS framework is suitable for our
sensing context, while the Bayesian framework enables a nat-
ural implementation of the relation linking the wavenumbers
from one frequency to the next.

The performance of the proposed method is assessed on simu-
lations. Our algorithm makes possible an accurate wavenumber
estimation using a small number of hydrophones, even in noisy
contexts. It outperforms other state-of-the-art methods such as
LS SFT (which does not benefit from any physical hypothe-
sis) and OMP (which benefits only from the sparse hypothesis):
estimated wavenumbers are more accurate, and the estimated
*f* − *k* diagrams are globally less noisy. Also, the method is val-
idated on experimental data recorded in the North Sea. In this
context, the proposed method performs the estimation of the
wavenumbers with a 20-hydrophone HLA.

Because the proposed method is based on CS, it is natu-
rally suitable for cases where the number of measurements
(i.e., hydrophones) is small and, ideally, randomly distributed.
However, it does not necessarily help in reducing the required
range aperture (i.e., the HLA length). As a result, the method
may not show its full strength on a small HLA with constant
hydrophone spacing. However, interesting perspectives arise for
analyzing (existing) synthetic aperture data, where the consid-
ered snapshots can be chosen at will along the source/receiver
track. Also, the potential of CS should be taken into account
when designing new at-sea experiment.

REFERENCES

[1] T. C. Yang, "Effectiveness of mode filtering: A comparison of matched-field and matched-mode processing," *J. Acoust. Soc. Amer.*, vol. 87, 1990.

[2] G. Wilson, R. Koch, and P. Vidmar, "Matched mode localization," *J. Acoust. Soc. Amer.*, vol. 84, pp. 310–320, 1988.

[3] J. Buck, J. Preisig, and K. Wage, “A unified framework for mode filtering and the maximum a posteriori mode filter,” *J. Acoust. Soc. Amer.*, vol. 103, p. 1813, 1998.

[4] T. Neilsen and E. Westwood, "Extraction of acoustic normal mode depth functions using vertical line array data," *J. Acoust. Soc. Amer.*, vol. 111, pp. 748–756, 2002.

[5] G. Potty, J. Miller, J. Lynch, and K. Smith, “Tomographic inversion for sediment parameters in shallow water,” *J. Acoust. Soc. Amer.*, vol. 108, pp. 973–986, 2000.

[6] J. Bonnel, B. Nicolas, J. Mars, and S. Walker, "Estimation of modal group velocities with a single receiver for geoacoustic inversion in shallow water," *J. Acoust. Soc. Amer.*, vol. 128, pp. 719–727, 2010.

[7] J. Bonnel, G. Le Touzé, B. Nicolas, and J. Mars, "Physics-based time-frequency representations for underwater acoustics: Power class utilization with waveguide-invariant approximation," *IEEE Signal Process. Mag.*, vol. 30, no. 6, pp. 120–129, 2013.

[8] C. Gervaise, B. Kinda, J. Bonnel, Y. Stéphan, and S. Vallez, “Passive geoacoustic inversion with a single hydrophone using broadband ship noise,” *J. Acoust. Soc. Amer.*, vol. 131, no. 3, pp. 1999–2010, 2012.

[9] Y. Le Gall and J. Bonnel, “Passive estimation of the waveguide invariant per pair of modes,” *J. Acoust. Soc. Amer.*, vol. 134, no. 2, pp. EL230-EL236, 2013.

[10] G. Frisk and J. Lynch, "Shallow water waveguide characterization using the Hankel transform," *J. Acoust. Soc. Amer.*, vol. 76, p. 205, 1984.

[11] G. V. Frisk et al., "Modal mapping experiment and geoacoustic inversion using sonobuoys," *IEEE J. Ocean. Eng.*, vol. 40, no. 3, pp. 607–620, 2015.

[12] Ö. Yilmaz, *Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data*. no. 10, Tulsa, OK, USA: Society of Exploration Geophysicists, 2001, ch. 1.2.

[13] B. Nicolas, J. Mars, and J. Lacoume, “Geoacoustical parameters estimation with impulsive and boat-noise sources,” *IEEE J. Ocean. Eng.*, vol. 28, no. 3, pp. 494–501, 2003.

[14] S. M. Kay, "Modern spectral estimation, theory and application," in *Signal Processing*. Upper Saddle River, NJ, USA: Prentice-Hall, 1988, ch. 6.

[15] E. C. Shang, H. P. Wang, and Z. Y. Huang, "Waveguide characterization and source localization in shallow water waveguides using the prony method," *J. Acoust. Soc. Amer.*, vol. 83, no. 1, pp. 103–108, 1988.

[16] K. M. Becker and G. V. Frisk, "Evaluation of an autoregressive spectral estimator for modal wave number estimation in range-dependent shallow water waveguides," *J. Acoust. Soc. Amer.*, vol. 120, no. 3, pp. 1423–1434, 2006.

[17] F. D. Philippe, P. Roux, and D. Cassereau, “Iterative high-resolution wavenumber inversion applied to broadband acoustic data,” *IEEE Trans. Ultrason. Ferroeletr. Freq. Control*, vol. 55, no. 10, pp. 2306–2311, 2008.

[18] F. Le Courtois and J. Bonnel, "Autoregressive model for high-resolution wavenumber estimation in a shallow water environment using a broadband source," *J. Acoust. Soc. Amer.*, vol. 135, no. 4, pp. EL199-EL205, 2014.

[19] I.-T. Lu, R. C. Qiu, and J. Kwak, "A high-resolution algorithm for complex spectrum search," *J. Acoust. Soc. Amer.*, vol. 104, no. 1, pp. 288–299, 1998.

[20] S. D. Rajan and S. D. Bhatta, "Evaluation of high-resolution frequency estimation methods for determining frequencies of eigenmodes in shallow water acoustic field," *J. Acoust. Soc. Amer.*, vol. 93, no. 1, pp. 378–389, 1993.

[21] E. J. Candès and M. B. Wakin, "An introduction to compressive sampling," *IEEE Signal Process. Mag.*, vol. 25, no. 2, pp. 21–30, 2008.

[22] N. R. Chapman and I. Barrodale, "Deconvolution of marine seismic data using the 11 norm," *Geophys. J. Int.*, vol. 72, no. 1, pp. 93–100, 1983.

[23] P. Gerstoft, A. Xenaki, and C. F. Mecklenbräuker, "Single and multiple snapshot compressive beamforming," arXiv:1503.02339, 2015.

[24] J. B. Harley and J. M. Moura, "Dispersion curve recovery with orthogonal matching pursuit," *J. Acoust. Soc. Amer.*, vol. 137, no. 1, pp. EL1-EL7, 2015.

[25] F. Le Courtois and J. Bonnel, "Compressed sensing for wideband wavenumber tracking in dispersive shallow water," *J. Acoust. Soc. Amer.*, vol. 138, no. 2, pp. 575–583, 2015.

[26] F. Le Courtois and J. Bonnel, "Wavenumber tracking in a low resolution frequency-wavenumber representation using particle filtering," in *Proc. IEEE Int. Conf. Acoust. Speech Signal Process.*, Florence, Italy, Jun. 2014, pp. 6805–6809.

[27] P. M. Djuric et al., "Particle filtering," *IEEE Signal Process. Mag.*, vol. 20, no. 5, pp. 19–38, Sep. 2003.

[28] C. Soussen, J. Idier, D. Brie, and J. Duan, "From Bernoulli-Gaussian deconvolution to sparse signal restoration," *IEEE Trans. Signal Process.*, vol. 59, no. 10, pp. 4572–4584, 2011.
[29] A. Drémeau, C. Herzet, and L. Daudet, “Boltzmann machine and mean-field approximation for structured sparse decompositions,” *IEEE Trans. Signal Process.*, vol. 60, no. 7, pp. 3425–3438, 2012.

[30] F. Jensen, W. Kuperman, M. Porter, and H. Schmidt, *Computational Ocean Acoustics*. New York, NY, USA: American Institute of Physics, 2nd ed., 2011, ch. 5 and 10.

[31] E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” *IEEE Signal Process. Mag.*, vol. 25, no. 2, pp. 21–30, 2008.

[32] D. L. Donoho, “Compressed sensing,” *IEEE Trans. Inf. Theory*, vol. 52, no. 4, pp. 1289–1306, Apr. 2006.

[33] B. K. Natarajan, “Sparse approximate solutions to linear systems,” *SIAM J. Comput.*, vol. 24, pp. 227–234, Apr. 1995.

[34] S. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” *IEEE Trans. Signal Process.*, vol. 41, pp. 3397–3415, Dec. 1993.

[35] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in *Proc. 27th Annu. Asilomar Conf. Signals Syst. Comput.*, Nov. 1993, pp 40–44.

[36] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” *SIAM J. Sci. Comput.*, vol. 20, pp. 33–61, 1998.

[37] R. Tibshirani, “Regression shrinkage and selection via the LASSO,” *J. Roy. Stat. Soc.*, vol. 58, pp. 267–288, 1996.

[38] B. A. Olshausen and D. J. Field, “Sparse coding with an overcomplete basis set: A strategy employed by v1?” *Vis. Res.*, vol. 37, no. 23, pp. 3311–3325, 1997.

[39] F. Krzakala, M. Mézart, F. Sausset, Y. F. Sun, and L. Zdeborová, “Statistical physics-based reconstruction in compressed sensing,” 2011. [Online]. Available at: http://arxiv.org/abs/1109.4424v2

[40] C. Herzet and A. Drémeau, “Bayesian pursuit algorithms,” in *Proc. 18th Eur. Signal Process. Conf.*, Aug. 2010, pp. 1474–1478.

[41] M. Beal, “Variational algorithms for approximate Bayesian inference,” Ph.D. dissertation, Univ. College London, London, U.K., May 2003.

[42] F. Krzakala, A. Manoel, E. W. Tramel, and L. Zdeborová, “Variational free energies for compressed sensing,” in *Proc. IEEE Int. Symp. Inf. Theory*, Jul. 2014, pp. 1499–1503.

[43] Y. C. Pati, R. Rezaiifar, and P. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” in *Conf. Rec. IEEE 27th Asilomar Conf. Signals Syst. Comput.*, 1993, pp. 40–44.

[44] J. Bonnel, B. Nicolas, and J. I. Mars, “Estimation of modal group velocities with a single receiver for geoacoustic inversion in shallow water,” *J. Acoust. Soc. Amer.*, vol. 128, no. 2, pp. 719–727, 2010.

[45] B. Nicolas, J. I. Mars, and J.-L. Lacoume, “Geoacoustical parameters estimation with impulsive and boat-noise sources,” *IEEE J. Ocean. Eng.*, vol. 28, no. 3, pp. 494–501, 2003.

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

**Angélique Drémeau** received the State Engineering degree from Telecom Bretagne, Brest, France, in 2007 and the M.Sc. and the Ph.D. degrees in signal processing and telecommunications from the Université de Rennes, Rennes, France, in 2007 and 2010, respectively.

She is currently an Associate Professor at ENSTA Bretagne, Lab-STICC (UMR 6285), Brest, France. Her research interests include inverse problems, sparse representation algorithms, and underwater acoustics.

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

**Florent Le Courtois** received the Ph.D. degree in acoustics from Université du Maine, Le Mans, France, in 2012.

From 2008 to 2011, he worked for the SNCF, Paris, France, on the problem of imaging acoustic sources on high-speed trains. From 2012 to 2015, he was employed as a Postdoctoral Fellow to work on the modal estimation using hydrophone arrays at ENSTA Bretagne, Brest, France. In January 2016, he joined Shom, Brest, France, to develop the underwater noise pollution descriptor for the MSFD.

His research interests include array processing, time-frequency analysis, sound-scape, and bioacoustics.

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

**Julien Bonnel** (S’08–M’11) received the Ph.D. degree in signal processing from Grenoble Institut National Polytechnique (Grenoble INP), Grenoble, France, in 2010.

Since 2010, he has been an Associate Professor at Lab-STICC (UMR 6285), ENSTA Bretagne, Brest, France. His research in signal processing and underwater acoustics includes time-frequency analysis, source detection/localization, geoacoustic inversion, acoustical tomography, passive acoustic monitoring, and bioacoustics.

Dr. Bonnel is a member of the Acoustical Society of America.
# Modal estimation in underwater acoustics by
data-driven structured sparse decompositions

Clément Dorffer*, Thomas Paviet-Salomon*, Gilles Le Chenadec* and Angélique Drémeau*

*Lab-STICC & ENSTA-Bretagne, UMR CNRS 6285, Brest, France.

E-mail: clement.dorffer@ensta-bretagne.fr, thomas.paviet-salomon@ensta-bretagne.org,
gilles.le_chenadec@ensta-bretagne.fr, angelique.dremeau@ensta-bretagne.fr

***Abstract***—In underwater acoustics, shallow water environments act as modal dispersive waveguides when considering lowfrequency sources. In this context, propagating signals can be described as a sum of few modal components, each of them propagating according to its own wavenumber. Estimating these wavenumbers is of key interest to understand the propagating environment as well as the emitting source. To solve this problem, we proposed recently a Bayesian approach exploiting a sparsity-inforcing prior. When dealing with broadband sources, this model can be further improved by integrating the particular dependence linking the wavenumbers from one frequency to the other. In this contribution, we propose to resort to a new approach relying on a restricted Boltzmann machine, exploited as a generic structured sparsity-inforcing model. This model, derived from deep Bayesian networks, can indeed be efficiently learned on physically realistic simulated data using well-known and proven algorithms.

***Index Terms***—Underwater modal estimation, structured sparse approximation, Restricted Boltzmann Machine, Bayesian algorithm.

<figure><img src="image_30.png" /><figcaption>Fig. 1. Illustration of a (f-k) diagram obtained with a Pekeris waveguide.</figcaption></figure>

## I. INTRODUCTION

In underwater acoustics, shallow environments behave like dispersive waveguides when considering low-frequency sources. An acoustic field received on an antenna is then classically described by a small set of modes propagating longitudinally according to their horizontal wavenumbers. The knowledge of these modes is of great importance for the characterization of the observation environment and, consequently, for the source localization. Among the different methods used to discriminate these modal components, the frequency-wavenumber (f-k) representation (see Fig. 1) allows a direct observation of the dispersion (i.e., the frequency dependence) of the wavenumbers. Inherently conceivable for a horizontal array of sensors aligned with the source, they are particularly used in geophysics [1]. Recent contributions have focused on the construction of (f-k) diagrams by exploiting less constrained acquisition schemes, allowing their use in underwater acoustics.

Since propagation is described by a small number of modes, the use of sparsity-inforcing models seems appropriate. In fact, some contributions (see e.g., [2], [3]) have proposed the use of the “compressed sensing” paradigm to estimate modal dispersion. However, if these methods prove to be relevant, we argue that they can be further improved by precisely integrating the dispersion relation linking the wavenumbers

from one frequency to another into the estimation process of the (f-k) diagram. To this end, we propose to build on some recent contributions in compressed sensing and machine learning.

On the one hand, research on compressed sensing [4]–[6] and more generally sparse decompositions have underlined the interest of taking into account the structures naturally living in signal representations [7]. On the other hand, restricted Boltzmann machines (RBM) [8] are at the heart of many recent contributions in machine learning that emphasize two precious qualities: i) they have been identified as generic models approximating any distribution over {0, 1}ⁿ [9]; ii) efficient algorithms have been developed to train them making them powerful representational models when large data sets are available [10]–[12].

These parallel researches in finally not so distant fields have recently led to the idea of exploiting RBMs learned from large databases to model the structures of sparse representations [13]–[15]. This paper follows on from these contributions. More particularly, we propose to deal with the (f-k) diagram estimation using a new Bayesian algorithm. Our approach considers a RBM as a model for structured sparse representations and exploits it through a mean-field approximation.

Thanks to DGA Naval Systems for funding.
## II. ACOUSTIC PROPAGATION IN SHALLOW WATER ENVIRONMENTS

In shallow water, acoustic propagation is described by modal theory. According to the latter, when considering a source at depth $z_s$ and frequency $f$, the signal received by an hydrophone placed at depth $z$ and distance $r$ of the source can be formulated as

$$y(f, r, z, z_s) = \sum_{m=1}^{K(f)} A_m(f, z, z_s) e^{-irk_{rm}(f)} + n(f, r, z), \quad (1)$$

where $K(f)$ stands for the number of modal components propagating at frequency $f$, $A_m(f, z, z_s)$ is the modal amplitude associated to the $m$-th component, $k_{rm}(f)$ is the $m$-th wavenumber at frequency $f$ and $n(f, r, z)$ stands for some additive measurement noise.

For a $M$ sensor linear antenna, Eq. (1) can be written as

$$y_f = \mathbf{F}z_f + \mathbf{n}_f \quad (2)$$

where $\mathbf{F}$ is a $(M \times N)$-dictionary of Fourier discrete atoms, $z_f$ contains a few non-zero coefficients indicating the wavenumbers propagating at frequency $f$ and quantifying their amplitudes and $\mathbf{n}_f$ stands for the additive noise.

A common approach to estimate $z_f$ from $y_f$ is to use a simple (inverse) spatial Fourier transform. Although very popular, the method requires sensors that are finely spaced sufficiently to avoid aliasing. As an alternative, more recent contributions [2], [16], [17] have proposed to leverage the sparsity of $z_f$ to design more robust methods.

When dealing with broadband sources, the generative model can be written as follows:

$$[\mathbf{y}_{f1}, \dots, \mathbf{y}_{fP}] = \mathbf{F} \cdot [\mathbf{z}_{f1}, \dots, \mathbf{z}_{fP}] + [\mathbf{n}_{f1}, \dots, \mathbf{n}_{fP}], \quad (3)$$

where $P$ is the number of (discretized) frequencies of the source. The estimation of the $z_f$'s can then be performed frequency per frequency with the same tools as for the mono-frequency case. Stacking these estimates $z_f$ one on top of the other leads to the (f-k) diagram (see Fig. 1) and makes appear some structures linking the wavenumbers propagating from one frequency to the other. These structures are well-known in practice under the physical name of dispersion relation. While some contributions in the literature work on an explicit combination of this relation with sparse approaches [2], [17], we propose here to resort to a Bayesian formulation of the (f-k) diagram estimation problem using a particular choice of prior holding the structure information.

To that end, a vectorized version of Eq. (3) is built by defining $\mathbf{y} = [\mathbf{y}_{f1}^T, \dots, \mathbf{y}_{fP}^T]^T$, $\mathbf{z} = [\mathbf{z}_{f1}^T, \dots, \mathbf{z}_{fP}^T]^T$, $\mathbf{n} = [\mathbf{n}_{f1}^T, \dots, \mathbf{n}_{fP}^T]^T$ and $\mathbf{D}$ as a block-diagonal matrix with $\mathbf{F}$ $P$-times repeated on its diagonal. Adopting this formulation, Eq. 3 can be re-formulated as

$$\mathbf{y} = \mathbf{D}\mathbf{z} + \mathbf{n}, \quad (4)$$

where $\mathbf{y}$ (resp. $\mathbf{z}$) is of dimension $MP$ (resp. $NP$) and $\mathbf{D}$ is a $(MP \times NP)$-dictionary.

The inverse problem we are interested in is then formulated as the estimation of the (f-k) diagram $\mathbf{z}$ from the acoustic signal $\mathbf{y}$ measured over the $M$-sensor antenna at all of the $P$ frequencies.

## III. PROBLEM FORMULATION

According to the modal theory, in shallow water environments and at low frequencies, the vector $\mathbf{z}$ has few non-zero elements, corresponding to the propagating modal wavenumbers, that is $\mathbf{z}$ is assumed to be sparse. Modeling such a property in a Bayesian framework can be realized in different ways. A popular model considers $\mathbf{z}$ as the result of the Hadamard product of a centered circular Gaussian variable<fnref n="1" />, say $\mathbf{x} \in \mathbb{C}^{NP}$ with variance $\sigma_x^2$ and a binary variable, say $\mathbf{s} \in \{0, 1\}^{NP}$. This model explicitly expresses the support $\mathbf{s}$ of the sparse vector $\mathbf{z}$, which is a key component of sparse estimation<fnref n="2" />: we look for a good estimation of $\mathbf{s}$ knowing the observations $\mathbf{y}$.

When no structure is assumed in sparse representations, the Bernoulli law constitutes a natural choice for $\mathbf{s}$. Structures, for their part, can be modeled according to various distributions. In [7], [18]–[20], Boltzmann machines are envisaged as generic models encompassing many well-known models. These models depend on parameters which are difficult to train, due - among others - to the presence of a quadratic term. Instead, more recent works [13]–[15] propose to exploit restricted Boltzmann machines, for which efficient training algorithms exist.

Formally, for $\mathbf{s} \in \{0, 1\}^{NP}$, RBM can be expressed as

$$p(\mathbf{s}) = \sum_{\mathbf{h}} p(\mathbf{s}, \mathbf{h}) \propto \exp(\mathbf{a}^T\mathbf{h} + \mathbf{b}^T\mathbf{s} + \mathbf{s}^T\mathbf{W}\mathbf{h}), \quad (5)$$

where $\mathbf{h}$ is a $L$-dimensional binary hidden variable and $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{W}$ are the RBM parameters.

In [14], the authors propose to implement this model into an OMP-like framework, while in [15], a reweighted $\ell_1$-like procedure is considered. Both procedures are deterministic, integrating RBM as an additional extra block and making hard decision on it. In [13], a Bayesian approach based on an approximate message passing procedure is proposed, leading to a more integrated vision of the restricted Boltzmann machine. However, RBM is not fully exploited in the sense that the support $\mathbf{s}$ is still assumed to follow a Bernoulli law, although informed by the RBM. In this paper, we propose a new Bayesian approach, leveraging both $i)$ the probabilistic nature of the RBM, and $ii)$ the explicit modeling of $\mathbf{s}$ as the commonly visible layer of a RBM (see Eq. (5)).

More particularly, considering model (4)-(5), we are interested in the following marginalized maximum a posteriori estimation problem:

$$(\hat{\mathbf{x}}, \hat{\mathbf{s}}) = \underset{\mathbf{x},\mathbf{s}}{\operatorname{argmax}} \log p(\mathbf{x}, \mathbf{s}|\mathbf{y}) \quad (6)$$

<footnote n="1">¹For a sake of simplicity, we will use the same notation for a random variable and its realizations.</footnote>

<footnote n="2">²We note indeed that once the support is estimated, the corresponding coefficients emerge straightforwardly by least-squares estimation.</footnote>
with $p(\mathbf{x}, \mathbf{s}|\mathbf{y}) = \int p(\mathbf{x}, \mathbf{s}, \mathbf{h}|\mathbf{y})$. In the continuation of [7], we propose to resort to a mean-field approximation.

## IV. DEEP STRUCTURED SOFT BAYESIAN PURSUIT

Mean-field (MF) approximations aim at approximating a posterior distribution, here $p(\mathbf{x}, \mathbf{s}|\mathbf{y})$, by a “simpler” distribution, say $q(\mathbf{x}, \mathbf{s})$, having a easy-to-handle factorization. We consider in this paper the factorization

$$q(\mathbf{x}, \mathbf{s}, \mathbf{h}) = \prod_{i=1}^{NP} q(x_i, s_i) \prod_{j=1}^{L} q(h_j), \quad (7)$$

with $x_i$ (resp. $s_i$) the $i$-th element in $\mathbf{x}$ (resp. $\mathbf{s}$) and $h_j$ the $j$-th element in $\mathbf{h}$. The approximation $q(\mathbf{x}, \mathbf{s}, \mathbf{h})$ is chosen to be as close as possible to $p(\mathbf{x}, \mathbf{s}, \mathbf{h}|\mathbf{y})$ in the sense of the Kullback-Leibler (KL) divergence. In practice, this latter optimization problem can be efficiently solved by an iterative algorithm, called “variational Bayesian expectation maximization” (VBEM) algorithm, insuring a decreasing of the KL divergence at each iteration. We expose here below its particularization to model (4)-(5).

The VBEM procedure successively updates all factors in the MF approximation, namely here the $q(x_i, s_i)$'s and $q(h_j)$'s, according to the following iterative rules:

$$q^{(k+1)}(x_i, s_i) \propto \exp \left( \langle \log p(\mathbf{x}, \mathbf{s}, \mathbf{h}, \mathbf{y}) \mid \begin{array}{l} \Pi_j q^{(k)}(h_j) \\ \Pi_{j>i} q^{(k)}(x_j, s_j) \\ \Pi_{j<i} q^{(k+1)}(x_j, s_j) \end{array} \right)$$

$$q^{(k+1)}(h_i) \propto \exp \left( \langle \log p(\mathbf{x}, \mathbf{s}, \mathbf{h}, \mathbf{y}) \mid \begin{array}{l} \Pi_{j>i} q^{(k)}(h_j) \\ \Pi_{j<i} q^{(k+1)}(h_j) \\ \Pi_j q^{(k+1)}(x_j, s_j) \end{array} \right)$$

where $\langle f(\mathbf{u}) \rangle_{q(u_j)} \triangleq \int_{u_j} q(u_j) f(\mathbf{u}) du_j$ and $k$ is the current iteration number. Developing these update rules according to model (4)-(5) gives $q^{(k+1)}(x_i, s_i) = q^{(k+1)}(x_i|s_i)q^{(k+1)}(s_i)$

$$\text{with } q^{(k+1)}(x_i|s_i) = \mathcal{CN}(m_x^{(k+1)}(s_i), \Sigma_x^{(k+1)}(s_i)) \quad (8)$$

$$m_x^{(k+1)}(s_i) = s_i \frac{\sigma_x^2}{\sigma_n^2 + \sigma_x^2 s_i \mathbf{d}_i^H \mathbf{d}_i} \mathbf{d}_i^H \langle \mathbf{r}_i \rangle^{(k+1)} \quad (9)$$

$$\Sigma_x^{(k+1)}(s_i) = \frac{\sigma_n^2 \sigma_x^2}{\sigma_n^2 + \sigma_x^2 s_i \mathbf{d}_i^H \mathbf{d}_i} \quad (10)$$

$$\langle \mathbf{r}_i \rangle^{(k+1)} = \mathbf{y} - \sum_{j \neq i} q^{(n)}(s_j=1) m_x^{(n)}(s_j=1) \mathbf{d}_j \quad (11)$$

with $n=k+1$ if $j<i$ and $n=k$ if $j>i$

where $\mathcal{CN}(\mu, \Gamma)$ stands for the circular Gaussian distribution with mean $\mu$ and variance $\Gamma$ and

$$q^{(k+1)}(s_i) \propto \exp \left( s_i \left( b_i + \sum_l w_{il} q^{(k)}(h_l = 1) \right) \right) \sqrt{\Sigma_x^{(k+1)}(s_i)} \exp \left( \frac{1}{2} \frac{|m_x^{(k+1)}(s_i)|^2}{\Sigma_x^{(k+1)}(s_i)} \right) \quad (12)$$

$$q^{(k+1)}(h_l) \propto \exp \left( h_l \left( a_l + \sum_i w_{il} q^{(k+1)}(s_i = 1) \right) \right) \quad (13)$$

where we have all along denoted by $\mathbf{d}_i$ the $i$-th column in $\mathbf{D}$ and $w_{il}$ the $(i, l)$-th element in $\mathbf{W}$.

To help the convergence of the algorithm, an estimation of the noise variance $\sigma_n^2$ is implemented in the same way as the one proposed in SoBaP [7], leading to a dependence of $\Sigma_x^{(k+1)}(s_i)$ in the iteration number $(k+1)$ (see Eq. (10)). As stopping criterion, the KL divergence is computed at each iteration, between the target distribution $p(\mathbf{x}, \mathbf{s}, \mathbf{h}|\mathbf{y})$ and its current MF approximate $q(\mathbf{x}, \mathbf{s}, \mathbf{h})$: the algorithm stops when this divergence no longer decreases “sufficiently”. Pseudo-code 1 gives a practical implementation of the above equations. The use of RBMs being a natural bridge towards deep networks, we will refer to the proposed procedure as the “Deep Structured Soft Bayesian Pursuit” (DSSoBaP).

### Pseudo-code 1 DSSoBaP algorithm

**Input:** $\mathbf{y}, \mathbf{D}, \mathbf{W}, \mathbf{a}, \mathbf{b}, \sigma_x^2, \sigma_n^2$

**Initialisation:** $\{q^{(0)}(s_i), m_x^{(0)}(s_i)\}_{s_i \in \{0,1\}, i \in \{1,...,NP\}}, \{q^{(0)}(h_l)\}_{l \in \{1,...,L\}}, k = 0, \text{KL} = 0, \text{KL}_{\text{old}} = \infty$

1. **while** $\text{KL}_{\text{old}} - \text{KL} > 10^{-1}$ **do**

2. **Optional:** update $\sigma_n^2$ according to [7]

3. **for** $i = 1...NP$ **do**

4. update $\langle \mathbf{r}_i \rangle^{(k+1)}$ using (11)

5. **Optional:** if step 2, update $\Sigma_x^{(k+1)}(s_i)$ using (10)

6. update $m_x^{(k+1)}(s_i)$ using (9) for $s_i \in \{0,1\}$

7. update $q^{(k+1)}(s_i)$ using (12) for $s_i \in \{0,1\}$

8. **end for**

9. **for** $l = 1...L$ **do**

10. update $q^{(k+1)}(h_l)$ using (13) for $h_l \in \{0,1\}$

11. **end for**

12. $k \leftarrow k + 1$

13. $\text{KL}_{\text{old}} \leftarrow \text{KL}$

14. compute current KL divergence

15. **if** $\text{KL}_{\text{old}} - \text{KL} \le 10^{-1}$ **then**, $k_{\text{max}} = k$

16. **end if**

17. **end while**

**Output:** $\{q^{(k_{\text{max}})}(s_i = 1), m_x^{(k_{\text{max}})}(s_i = 1)\}_{i \in \{1,...,NP\}}$

Note that DSSoBaP has an algorithmic complexity in the order of $\mathcal{O}(NP(MP+L))$ per iteration. This complexity is reasonable compared to other procedures in the literature, such as e.g., RBM-OMP [14]. This drawback in particular prevented us from using RBM-OMP in our particular problem.

## V. EXPERIMENT

In this section, we confront our approach to synthetic experiments and compare its performance to the SoBaP standard procedure.

### A. Data simulation and training phase

Several models of the ocean have been proposed in the literature of modal theory [21]. One of the most popular is the Pekeris waveguide. In this model, the sea surface is assumed to be perfectly reflective, and the sea bottom is considered as a semi-infinite fluid. Its parameters therefore
<figure><img src="image_31.png" /><figcaption>Fig. 2. (f-k) diagrams obtained for $M/N = 0.1$ and SNR= 3 dB by (left) inverse Fourier transform; (middle) SoBaP; (right) DSSoBaP.</figcaption></figure>

include water layer thickness $D$ (constant over the source-receiver distance), water density $\rho_w$ and celerity $c_w$, ground density $\rho_g$ and celerity $c_g$. Considering this simple model, we propose to simulate the propagation of an acoustic signal in various Pekeris environments and to learn the structures in the corresponding (f-k) diagrams. The considered setup is as follows.

Source and antenna are assumed to lie on the ground at same depth $z_s = z = D$ varying from 100 m to 200 m. We simulate the water celerity $c_w$ from $1450 \text{ m s}^{-1}$ to $1550 \text{ m s}^{-1}$, and the water density $\rho_w$ from $950 \text{ kg m}^{-3}$ to $1050 \text{ kg m}^{-3}$. The studied frequencies range from 0 Hz to 30 Hz with a step-size of 0.5 Hz, i.e., $P = 60$ frequencies are considered, and we set the ground properties $c_g$ and $\rho_g$ from the sediments geoacoustic parameters synthesis in [22]. The simulated antenna is composed of $M = 120$ sensors regularly spaced from 50 m.

With such parameters, we generated a set of 3575 signals – of size $120 \times 60$, i.e., 120 sensors, 60 frequencies—and their theoretical sparse decomposition (f-k) diagrams – of size $120 \times 60$, i.e., $N = 120$ potential wavenumbers for each 60 frequencies – whose sparsity ratios fluctuate between $3.10^{-3}$ and $3.10^{-2}$. We then randomly divided these simulations into a training set (3000 samples) and a testing set (575 samples) and binarized the (f-k) diagrams from the training set so as to obtain their theoretical supports.

These supports were used to train a RBM with $L = 30$ hidden units. Training parameters were fixed to $10^{-2}$ for the learning rate, 50 for the minibatch size, 2000 for the number of epochs and the Persistent Contrastive Divergence algorithm [10] was employed with an $l_1$ regularization so as to inforce the sparsity of the learned features. The obtained RBM is then integrated as prior into DSSoBaP.

## B. Tests and results analysis

To assess the relevance of RBM into the considered inverse problem, we compare the performance of DSSoBaP with that of SoBaP, introduced in [7]. This algorithm is based on the same MF approximation as DSSoBaP but instead of the RBM prior (see Eq. (5)), SoBaP exploits a *i.i.d.* Bernoulli model on the SR support. It has to be noticed that we do not compare

the approach with those proposed in [7], [14], mainly because of their computational cost.

Two figures of merit are considered: i) the estimation of the modal amplitudes (i.e., the non-zero coefficients in the $z_f$'s) is evaluated according to the normalized mean square error (NMSE), ii) the detection of the wavenumbers (i.e., the support of the sparse $z_f$'s) is evaluated according to the True Positives Rate (TPR) which is the proportion of the theoretical support that has been correctly recovered and the False Discovery Rate (FDR) which is the proportion of the estimated support that should not have been activated (i.e., the proportion of false detection within the estimated support). Each metric is valued according to the signal-to-noise ratio (SNR) on the one hand, and the number of sensors on the other hand. For the tests related to the SNR, we fix the number of observations to 1440 by randomly selecting 24 sensors and let the SNR range from 0 dB to 5 dB by adding Gaussian noise $n_f$ with modulated variance $\sigma_n^2$. For the tests related to the number of sensors – which is classical in compressed sensing – we choose to fix the SNR to 3 dB and let the proportion of hydrophones range from 0.1 to 1 by randomly selecting from 12 to 120 sensors of the antenna.

SoBaP and DSSoBaP are run using the same stopping criterium as described in pseudo-code 1. For SoBaP, the Bernoulli parameter is set for each test sample to its exact value that is the theoretical proportion of nonzero coefficients in $z$, while for DSSoBaP the prior is provided by the same trained RBM for all test samples. Finally, the performance metrics are averaged over the entire test dataset.

Fig. 3 illustrates the obtained results. One can see that DSSoBaP performance curves follow the same tendencies as SoBaP while being better in most cases, especially in the worst configurations, i.e., with low SNR and low $M/N$ ratio. In particular, the gap between DSSoBaP and SoBaP TPR curves shows that the RBM-based prior has been more effective than the Bernoulli prior for recovering the theoretical support, while the similar FDR curves show that it also conserves the same sparsity properties of the Bernoulli prior by producing only slightly more false detections. Interestingly, we can observe that the improvement of the detection performance leads to an improvement of the NMSE metric, that is the amplitudes of the modes are consequently better recovered with DSSoBaP.
<figure><img src="image_32.png" /><figcaption>Fig. 3. NMSE (top) and TPR/FDR (bottom) w.r.t. the SNR (left) and the sensor ratio *M/N* (right).</figcaption></figure>

than with SoBaP.

An example of recovered (f-k) diagrams by inverse Fourier transform, SoBaP and DSSoBaP is shown on Fig. 2. The considered test parameters are SNR=3 dB and sensor ratio *M/N* = 0.1. One can see that DSSoBaP is able to detect much elements from the ground truth support than SoBaP. In particular, structures that are characteristic of modal propagation – the line shapes induced by wavenumber propagation along frequencies – are better recognized using the RBM-based prior, thus showing that RBM are, as expected, well suited for learning spatial dependencies in the support density.

## VI. CONCLUSION

In this work, we introduced a new Bayesian sparse decomposition algorithm able to exploit an RBM – learned in a previous step – as a structure model on the support of a sparse representation. The proposed algorithm relies on an MF approximation and allows a natural integration of the RBM. Its performance on a current underwater acoustic problem is very promising compared to its “standard” counterpart which exploits non-structured Bernoulli variables. Moreover, as the training phase is performed on a collection of physically-realistic simulations, it naturally incites the prior to respect some theoretical structures, which might be helpful when dealing with real datasets.

## REFERENCES

[1] L. Amundsen and B. Ursin, “Frequency-wavenumber inversion of acoustic data,” *Geophysics*, vol. 56, no. 7, pp. 1027–1039, 1991.

[2] A. Drémeau, F. Le Courtois, and J. Bonnel, “Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array,” *IEEE Journal of Oceanic Engineering*, vol. 42, no. 4, pp. 914–922, 2017.

[3] J. B. Harley and J. M. Moura, “Dispersion curve recovery with orthogonal matching pursuit,” *The Journal of the Acoustical Society of America*, vol. 137, no. 1, pp. EL1–EL7, 2015.

[4] D. L. Donoho, “Compressed sensing,” *IEEE Transactions on Information Theory*, vol. 52, no. 4, pp. 1289–1306, 2006.

[5] E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” *IEEE Signal Processing Magazine*, vol. 25, no. 2, pp. 21–30, 2008.

[6] Y. C. Eldar and G. Kutyniok, *Compressed sensing: theory and applications*. Cambridge university press, 2012.

[7] A. Drémeau, C. Herzet, and L. Daudet, “Boltzmann machine and mean-field approximation for structured sparse decompositions,” *IEEE Transactions On Signal Processing*, vol. 60, no. 7, p. 3425–3438, 2012.

[8] P. Smolensky, “Information processing in dynamical systems: Foundations of harmony theory,” Colorado Univ at Boulder Dept of Computer Science, Tech. Rep., 1986.

[9] N. Le Roux and Y. Bengio, “Representational power of restricted boltzmann machines and deep belief networks,” *Neural computation*, vol. 20, no. 6, pp. 1631–1649, 2008.

[10] T. Tieleman, “Training restricted boltzmann machines using approximations to the likelihood gradient,” in *Proc. 25th international conference on Machine learning*, 2008, pp. 1064–1071.

[11] G. E. Hinton, “A practical guide to training restricted boltzmann machines,” in *Neural networks: Tricks of the trade*. Springer, 2012, pp. 599–619.

[12] A. Fischer and C. Igel, “Training restricted boltzmann machines: An introduction,” *Pattern Recognition*, vol. 47, no. 1, pp. 25–39, 2014.

[13] E. W. Tramel, A. Manoel, F. Caltagirone, M. Gabrié, and F. Krzakala, “Inferring sparsity: Compressed sensing using generalized restricted boltzmann machines,” in *2016 IEEE Information Theory Workshop (ITW)*. IEEE, 2016, pp. 265–269.

[14] L. F. Polania and K. E. Barner, “Exploiting restricted boltzmann machines and deep belief networks in compressed sensing,” *IEEE Transactions on Signal Processing*, vol. 65, no. 17, pp. 4538–4550, 2017.

[15] Z. Liao, J. Zhang, D. Hu, C. Li, L. Zhu, and Y. Li, “Sparse signal reconstruction with statistical prior information: A data-driven method,” *IEEE Access*, vol. 7, pp. 157037–157045, 2019.

[16] F. Le Courtois and J. Bonnel, “Compressed sensing for wideband wavenumber tracking in dispersive shallow water,” *The Journal of the Acoustical Society of America*, vol. 138, no. 2, pp. 575–583, 2015.

[17] T. Paviet-Salomon, C. Dorffer, J. Bonnel, B. Nicolas, T. Chonavel, and A. Drémeau, “Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics,” *Proc. Int’l Conference on acoustics, speech and signal processing*, 2019.

[18] T. Peleg, Y. Eldar, and M. Elad, “Exploiting statistical dependencies in sparse representations for signal recovery,” *IEEE Transactions on Signal Processing*, vol. 60, no. 5, pp. 2286–2303, 2012.

[19] V. Cevher, M. Duarte, C. Hegde, and R. Baraniuk, “Sparse signal recovery using markov random fields,” *Proc. Neural Information Processing Systems (NeurIPS)*, 2008.

[20] J. Garrigues and B. A. Olshausen, “Learning horizontal connections in a sparse coding model of natural images,” *Proc. Neural Information Processing Systems (NeurIPS)*, 2008.

[21] F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt, *Computational ocean acoustics*. Springer Science & Business Media, 2011.

[22] D. Jackson, “High-frequency ocean environmental acoustic models handbook,” *Tech. report, Applied Physics Laboratory, Univ. of Washington*, 1994.
Reference-less measurement of the
transmission matrix of a highly
scattering material using a DMD and
phase retrieval techniques

Angélique Drémeau,¹ Antoine Liutkus,² David Martina,³
Ori Katz,³,⁴ Christophe Schülke,⁵ Florent Krzakala,¹,⁶
Sylvain Gigan,³,⁶,* and Laurent Daudet⁴,⁵

¹LPS-ENS and CNRS UMR 8550, Paris, F-75005, France.

²Inria, CNRS, Loria UMR 7503, Villers-lès-Nancy, F-54600, France

³Laboratoire Kastler Brossel, Université Pierre et Marie Curie, Ecole Normale Supérieure,
Collège de France, CNRS UMR 8552, Paris, F-75005, France

⁴Institut Langevin, ESPCI and CNRS UMR 7587, Paris, F-75005, France

⁵Paris Diderot University, Sorbonne Paris Cité, Paris, F-75013, France

⁶Sorbonne Universités, UPMC Université Paris 06, F-75005, Paris, France

*sylvain.gigan@lkb.ens.fr

**Abstract:** This paper investigates experimental means of measuring the transmission matrix (TM) of a highly scattering medium, with the simplest optical setup. Spatial light modulation is performed by a digital micromirror device (DMD), allowing high rates and high pixel counts but only binary amplitude modulation. On the sensor side, without a reference beam, the CCD camera provides only intensity measurements. Within this framework, this paper shows that the TM can still be retrieved, through signal processing techniques of phase retrieval. This is experimentally validated on three criteria : quality of prediction, distribution of singular values, and quality of focusing.

© 2015 Optical Society of America

**OCIS codes:** (290.4210) Multiple scattering; (070.6120) Spatial light modulators; (100.5070) Phase retrieval

**References and links**

1. P. Sebbah, *Waves and Imaging Through Complex Media* (Springer, 2001).

2. A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, "Controlling waves in space and time for imaging and focusing in complex media," *Nature Photonics* **6**, 283–292 (2012).

3. M. Cui and C. Yang, "Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation," *Opt. Express* **18**, 3444–3455 (2010).

4. I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, "Focusing and scanning light through a multimode optical fiber using digital phase conjugation," *Opt. Express* **20**, 10583 (2012).

5. I. M. Vellekoop and A. P. Mosk, "Focusing coherent light through opaque strongly scattering media," *Opt. Express* **32**, 2309 (2007).

6. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, "Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media," *Phys. Rev. Lett.* **104**, 100601 (2010).

7. S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, "Image transmission through an opaque material," *Nature Commun.* **1**, 81 (2010).

8. Y. Choi, T. D. Yang, C. Fang-Yen, P. Kang, K. J. Lee, R. R. Dasari, M. S. Feld, and W. Choi, "Overcoming the diffraction limit using multiple light scattering in a highly disordered medium," *Phys. Rev. Lett.* **107**, 023902 (2011).
9. M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park, and W. Choi, "Maximal energy transport through disordered media with the implementation of transmission eigenchannels," Nature Photonics **6**, 581–585 (2012).

10. Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, "Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber," Phys. Rev. Lett. **109** (2012).

11. I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, "High-resolution, lensless endoscope based on digital scanning through a multimode optical fiber," Biomed. Opt. Express **4**, 260–270 (2013).

12. S. Bianchi and R. Di Leonardo, "A multi-mode fiber probe for holographic micromanipulation and microscopy," Lab on a Chip **12**, 635 (2012).

13. T. Čižmár and K. Dholakia, "Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics," Opt. Express **19**, 18871 (2011).

14. J. B. Sampsell, "Dmd display system," (1995). US Patent 5,452,024.

15. D. B. Conkey, A. M. Caravaca-Aguirre, and R. Piestun, "High-speed scattering medium characterization with application to focusing light through turbid media," Opt. Express **20**, 1733–1740 (2012).

16. S. A. Goorden, J. Bertolotti, and A. P. Mosk, "Superpixel-based spatial amplitude and phase modulation using a digital micromirror device," arXiv:1405.3893 [physics] (2014).

17. D. Akbulut, T. J. Huisman, E. G. van Putten, W. L. Vos, and A. P. Mosk, "Focusing light through random photonic media by binary amplitude modulation," Opt. Express **19**, 4017–4029 (2011).

18. D. Kim, W. Choi, M. Kim, J. Moon, K. Seo, S. Ju, and W. Choi, "Implementing transmission eigenchannels of disordered media by a binary-control digital micromirror device," Opt. Commun. **330**, 35–39 (2014).

19. J. W. Tay, J. Liang, and L. V. Wang, "Amplitude-masked photoacoustic wavefront shaping and application in flowmetry," Opt. Lett. **39**, 5499–5502 (2014).

20. X. Zhang and P. Kner, "Binary wavefront optimization using a genetic algorithm," Journal of Optics **16**, 125704 (2014).

21. T. Chaigne, O. Katz, A. C. Boccara, M. Fink, E. Bossy, and S. Gigan, "Controlling light in scattering media non-invasively using the photoacoustic transmission matrix," Nature Photonics **8**, 58–64 (2014).

22. D. Akbulut, T. Strudley, J. Bertolotti, T. Zehender, E. Bakkers, A. Lagendijk, W. Vos, O. Muskens, and A. Mosk, "Measurements on the optical transmission matrices of strongly scattering nanowire layers," in *Proc. Conf. on Lasers and Electro-Optics Europe and International Quantum Electronics Conf. (CLEO EUROPE/IQEC)*, (2013).

23. A. Drémeau and F. Krzakala, "Phase recovery from a bayesian point of view: the variational approach," in *Proceedings of IEEE Trans. Acoust. Speech Signal Process*. (2015).

24. R. Gerchberg and W. Saxton, "A practical algorithm for the determination of phase from image and diffraction plane pictures," Optik **35**, 237–246 (1972).

25. J. R. Fienup, "Phase retrieval algorithms: a comparison," Appl. Opt. **21**, 2758–2769 (1982).

26. E. Candès, T. Strohmer, and V. Voroninski, "Phaselift: exact and stable signal recovery from magnitude measurements via convex programming," Commun. in Pure and Applied Mathematics **66**, 1241–1274 (2013).

27. I. Waldspurger, A. d'Aspremont, and S. Mallat, "Phase recovery, maxcut and complex semidefinite programming," Mathematical Programming Series A- Springer (2013).

28. P. Schniter and S. Rangan, "Compressive phase retrieval via generalized approximate message passing," in *Proceedings of Communication, Control, and Computing (Allerton)* (2012).

29. V. A. Marčenko and L. A. Pastur, "Distribution of eigenvalues for some sets of random matrices." Mathematics of the USSR-Sbornik **1**, 457 – 483 (1967).

30. S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, "Controlling light through optical disordered media: transmission matrix approach," New Journal of Physics **13**, 123021 (2011).

31. A. Liutkus, D. Martina, S. Popoff, G. Chardon, O. Katz, G. Lerosey, S. Gigan, L. Daudet, and I. Carron, "Imaging with nature: Compressive imaging using a multiply scattering medium," Sci. Rep. **4** (2014).

# 1. Introduction

Wave propagation in complex media is a fundamental problem in physics, be it in acoustics, optics, or electromagnetism [1]. In optics, it is particularly relevant for imaging applications. Indeed, when light passes through a multiply scattering medium, such as a biological tissue or a layer of paint, ballistic light is rapidly attenuated, preventing conventional imaging techniques, and random scattering events generate a so-called speckle pattern that is usually considered useless for imaging. Recently, wavefront shaping using spatial light modulators (SLM) has emerged as a unique tool to manipulate multiply scattered coherent light, for focusing or imaging in scattering media [2]. In essence, these methods use the linearity and time-reversal symmetry of the wave propagation, whatever the complexity of the medium, to control the output speckle field, by manipulating the light beam impinging on the scattering sample. Dif-
ferent wavefront shaping approaches rely on digital phase-conjugation [3, 4] or iterative algo-
rithms [5], but it is also possible to measure the so-called transmission matrix (TM) of the
medium [6], which fully describes light propagation through the linear medium, from the mod-
ulator device to the detector. This approach has been particularly efficient for focusing, imag-
ing [7, 8] and for studying the transmission modes of the medium [9]. These methods are not
only valid for scattering material but can also be applied to other complex transmission system,
most notably multimode fibers, turning them into minimal footprint endoscopes [10–13].

A major limitation of most of these techniques for imaging is their speed. Indeed, the wave-
front shaping process must be faster than the stability time of the medium, which can be of
only a few milliseconds in biological tissues. Yet, most of the works reported so far have relied
on phase modulators which are usually slow (few tens of Hertz for liquid crystal modulators).
Micro Electro-Mechanical Systems (MEMS) modulators are much faster, but are usually not
phase-only. As a promising alternative for wave shaping in complex media, Digital Micromir-
ror Device (DMD) technology [14] offers binary amplitude modulators (i.e., “ON” or “OFF”)
operating at > 20kHz, with high pixel counts (10⁶) and low pitch (around 10 microns), all this
at low cost. These binary amplitude modulators have been used as phase modulators, using ap-
propriate diffraction and filtering, e.g. by Lee-type amplitude holography [15, 16], as shown on
Fig. 1(b). While phase control is more effective for wavefront shaping than amplitude control,
some works reported on using DMD as genuine binary amplitude modulators for wavefront
shaping through opaque scattering media, albeit usually yielding lower overall efficiency than
phase modulators for focusing or mode matching [17–19]. The DMD configuration can also be
optimized using genetic algorithms [20] to maximize the intensity enhancement.

For the measurement of a TM, an additional issue lies in accessing the amplitude and phase
of the output field, that in optics usually requires a holographic measurement, i.e., a reference
beam, as shown on Fig. 1(a). This reference beam can either be co-propagating in the medium
[7, 21], or use an external reference arm [8, 22]. The phase and amplitude of the measured
field can then be extracted by simple linear combinations of interference patterns with a phase-
shifted or off-axis reference. This however poses the unavoidable experimental problem of the
interferometric stability of the reference arm.

In this work we report on the full measurement of the complex TM of a multiply scatte-
ring medium, using a DMD binary amplitude modulator as an SLM, with no reference on the
detection side, as shown on Fig. 1(c). This approach combines the high-speed and high pixel
counts allowed by DMD devices, with the simplicity and robustness of a reference-less optical
setup. However, it involves advanced signal processing algorithms for phase retrieval, run on
a sufficiently large number of input-output calibration measurements. Here, we use a Bayesian
phase retrieval algorithm [23] for the estimation of a TM, based on actual noisy experimental
measurements. This is done one output pixel at a time, a process that can easily be paral-
lelized for speed. Although a number of phase estimation techniques could also be used, our
approach provides good results at moderate computational costs. This technique is validated
on experimental data, with three criteria. First, for a given arbitrary input pattern, we compare
the measured output and the predicted output computed with the measured TM. Then, we show
that the distribution of the singular values of the measured TM varies according to random ma-
trix theory. Finally, we demonstrate that single- or multi-point light focusing can be achieved.
Interestingly, our Bayesian formulation allows the use of the same model for the estimation
of the optimal DMD binary input pattern. In addition to being an interesting signal processing
problem, this approach is particularly relevant for real-life applications of the TM approach,
since it allows a simple, fast and robust implementation.
<figure><img src="image_33.png" /><figcaption>Fig. 1. Different experimental approaches for measuring the complex-valued transmission matrix of a scattering medium with a binary DMD amplitude modulator: (a) using a reference arm for retrieving the phase of the output field by off-axis or phase-shifting holography; (b) using the DMD as a spatial phase modulator by displaying amplitude holograms, and using the unmodulated parts of the field as a phase-stable reference; (c) presented approach where only intensity values are measured.</figcaption></figure>

## 2. Experimental setup

Our experimental setup, described in Fig. 2, uses a DMD-array from Texas Instrument (1920 × 1080 tilting micromirrors), driven by the DLP V-9500 VIS module (Vialux). The DMD is made of mirrors that can switch between two angular positions separated by 24°, thus reflecting each pixel either toward a beam dump (pixel “OFF”) or towards the focusing system (pixel “ON”).

Under Matlab, an amplitude mask is computed and loaded on the DMD. The pattern corresponding to the ON pixels is focused on the surface of a thick scattering medium by means of a $f = 100$ mm lens L1 (thus the DMD pixels correspond roughly to incidence angles on the sample). The sample is a ~100 microns thick layer of white paint, which is thick enough in order to considerably mix the light on the other side, producing a complex speckle interfering pattern. This speckle pattern is collected through a microscope objective (L2) and detected on a camera (AVT Pike F-100B). In order to measure the TM, we need to send a large series of input patterns (typically a few times the number of input pixels we wish to control), in a time over which the medium can be considered stationary. For this purpose, we use the “high speed” driver provided with the DMD in order to load all the to-be-projected random amplitude masks to the memory of the DMD driver module, and we trigger the display of each mask via a DAQ card (National Instruments, PCI-6221) and a waveform generator. In the same way, in order to be as fast as possible, we also only consider a subregion on the camera of size of 400 × 400 pixels. The overall acquisition rate is 31 images per second. To monitor the stability of the medium, we periodically measure the correlation of the speckle image corresponding to the same input mask. We therefore quantify the stability of the medium, which is better than 98%
over the total measurement time (typically around 5 minutes).

<figure><img src="image_34.png" /><figcaption>Fig. 2. Experimental scheme: A 532 nm CW laser is expanded through a telescope in order to obtain a homogeneous beam. Through a rectangular mask, it illuminates the DMD which acts as binary amplitude spatial light modulator. The DMD reflects the light in two different directions corresponding to either “ON” (unit transmission) or “OFF” (the light is deviated towards a beam dump). The transmitted pattern is focused by a first lens L1 on the scattering medium – here a white paint layer –, acting as a thick multiply scattering medium. The transmitted speckle pattern is collected by a microscope objective and is observed through a polarizer P on a CCD camera.</figcaption></figure>

### 3. Phase retrieval for optics through complex media

The experimental setup described in the previous section presents different challenges in order to exploit the TM framework: the absence of reference requires a solution exploiting only intensity measurements, and this must be robust to experimental noise. This paper shows that signal processing techniques for phase retrieval can be used, first for the estimation of the TM, and then for getting the optimal DMD configuration in a focalization task.

#### 3.1. Calibration as a phase retrieval problem

If the focalization task - or in a more general view, any front-shaping task using spatial light modulators - can be easily expressed as a phase retrieval problem - knowing the TM, which SLM setup best explains the intensity observations or the desired output -, this is not the case for the estimation of the TM.

We formalize the latter as a calibration problem: given P incoming waves, assumed perfectly known, which model explains at best the observed outputs?

Formally, let $x_{\mu} \in \{0, 1\}^N$ stand for the binary DMD inputs related to the $\mu$-th acquisition, where $N$ is the number of pixels (mirrors) used on the DMD. We assume that the partial observations of the sole moduli of the transmitted waves (the square root of the camera measured intensities), denoted by $y_{\mu} \in \mathbb{R}_{+}^{M}$, obey

$$y_{\mu} = |\mathbf{D}x_{\mu}|, \quad \forall \mu \in \{1, \dots, P\}, \qquad (1)$$
where **D** is the *unknown* complex-valued transmission matrix characterizing the scattering material, and *M* is the number of observed pixels on the camera.

Then, adopting a matrix formulation and conjugating-transposing the system, we get

$$ \mathbf{Y}^{\mathrm{H}} = |\mathbf{X}^{\mathrm{H}}\mathbf{D}^{\mathrm{H}}|, \qquad (2) $$

where $\mathbf{Y} = [\mathbf{y}_1, \dots, \mathbf{y}_P]$, $\mathbf{X} = [\mathbf{x}_1, \dots, \mathbf{x}_P]$ and $^{\mathrm{H}}$ denotes the conjugate-transpose of a matrix/vector. This reveals a “classic” phase retrieval problem: given the matrix of inputs $\mathbf{X}^{\mathrm{H}}$, each column of $\mathbf{Y}^{\mathrm{H}}$ is used to estimate each complex-valued column of $\mathbf{D}^{\mathrm{H}}$.

## 3.2. Phase retrieval techniques

The problem of reconstructing a complex vector given only the magnitude of measurements is a non-convex optimization problem notoriously difficult to solve. Many algorithms have been devised in the literature to deal with this problem, especially in optics where most approaches aim at reconstructing a phase object from an intensity image in a Fourier plane [24, 25]. Here we are interested in a somewhat simpler and more general problem where the transform is completely general and the output points can be treated independently. Amongst recent contributions to this signal processing problem, we can mention [26,27] which approximate the phase recovery problem by relaxed problems enabling the use of standard optimization procedures, and the Bayesian procedures proposed in [23, 28], which circumvent the non-linearity of the modulus through the introduction of hidden variables and resort to variational approximations.

In this paper, we chose to use the so-called *prVBEM* Bayesian approach introduced in [23]. Our choice is motivated by two of its main advantages:

* its generic framework, enabling its use in both applications we are here interesting in, namely the calibration of the TM and the focalization task,

* its reasonable computational complexity, of crucial importance for a complete procedure including the calibration of the TM and its immediate use for focusing.

We will not detail here the derivation of the algorithm. Its particularizations to the two problems of interest, calibration of the TM and focalization, are described resp. in appendix 7.1 and 7.2. For a methodological justification of the approach, we refer the reader to the original paper [23].

# 4. Estimating the TM with intensity-only measurements and binary inputs

In this section, we focus on the calibration stage of the proposed approach. We assess the accuracy of the TM estimated through the phase retrieval technique [23] with regard to two different analyses: its prediction performance, namely its ability to predict the observations, and its nature, expected to satisfy random features.

## 4.1. Prediction performance

To assess the prediction performance of the estimated TM, we adopt a cross-validation-like experimental framework. The setup is as follows. We measure the *M* = 40000 camera pixels stemming from *N* = 900 DMD mirrors, 50% of them being turned on, the others off at each displayed pattern. The operation is repeated randomly *P* = 6000 times. Given this dataset, a row of the TM is then learned from *p* = *αN* calibration measurements, with *α* varying in {1, ..., 6}, and used in a second step to predict the *P* – *p* remaining measurements. This estimation is performed on 50 different rows of the TM.

The *prVBEM* algorithm presents of complexity of order $O(p^2)$. Run for 200 iterations, this leads to a computational time of order 0.6 s per estimated row of the TM, in Matlab, on a Macbook Air with a 1.7Ghz i7 processor. This performance is good in comparison with other state-of-the-art algorithms (see e.g., [23]) and, keeping in mind that rows are independent, sufficient
for making possible a complete procedure “calibration of the TM-focalization” in a reasonable time.

<figure><img src="image_35.png" /><figcaption>Fig. 3. Prediction performance according to (a), the mean-square error (MSE), in log scale, and to (b), the normalized cross-correlation between observation predictions using the estimated TM, and actual measurements of the output moduli (square root of the camera intensity values), as a function of the number of calibration measurements (x-axis is $\alpha$, such that $p = \alpha N$ calibration measurements are used).</figcaption></figure>

The prediction performance of the estimated TM is evaluated according to the mean-square error (MSE) and the normalized cross-correlation between the moduli of the $P-p$ predicted measurements and the actual observed ones, such as, for the estimate $\hat{\mathbf{d}}$ of the current (conjugate) row of the TM (see Eq. (2)),

$$MSE = \frac{\|\mathbf{y}_{P-p} - |\mathbf{X}_{P-p}^H \hat{\mathbf{d}}||_2^2}{\|\mathbf{y}_{P-p}\|_2^2}, \qquad (3)$$

$$Corr = \frac{\mathbf{y}_{P-p}^H |\mathbf{X}_{P-p}^H \hat{\mathbf{d}}|}{\|\mathbf{y}_{P-p}^H \mathbf{y}_{P-p}\|_2 \|\hat{\mathbf{d}}^H \mathbf{X}_{P-p}^H \mathbf{X}_{P-p}^H \hat{\mathbf{d}}\|_2}, \qquad (4)$$

where $\mathbf{y}_{P-p}$ stands for the current (real-valued) row of $\mathbf{Y}$ (see Eq. (2)) restricted to its last $P-p$ elements, the first $p$ ones being used for the estimation of $\hat{\mathbf{d}}$.

These figures of merit give different insights into the estimations achieved by the algorithm: the MSE measures the distance (in an euclidian sense) to the actual observations, while the correlation evaluates their angular difference. Figure 3 show these quantities (resp. in (a) and (b)) averaged over the 50 rows of the TM considered for estimation (red lines), as well as the minimum and maximum (given by whiskers), for an increasing $\alpha = p/N$.

Both curves match admirably: increasing $\alpha$ leads to a decrease of the MSE and an increase of the correlation. Interestingly, we see that for $\alpha \ge 3$, that is, for at least 3 times more real measurements than complex unknowns, the estimation of the TM is accurate enough to predict the observations with an average correlation around 0.95 and an average MSE lower than 0.18 (Fig. 3(a) shows this value in log-scale).

## 4.2. Comparison of singular values to Random matrix theory

Interestingly, we can check that the measured TM presents some characteristics as predicted by random matrix theory. One practical way is to verify that the distribution of its normalized
<figure><img src="image_36.png" /><figcaption>Fig. 4. Density of the normalized singular values for different $\gamma = M/N$. Stamped line: experimental results, continuous line: Marčenko-Pastur law.</figcaption></figure>

singular values obeys the Marčenko-Pastur law [29]. It should be noted that such apparently random signals are the hardest case for phase retrieval, where no specific structure can be taken into account.

In order to reduce the influence of specifics of our experimental setting, we perform the following operations, as in [30]:

i) We normalize over the rows and columns, to attenuate the illumination artifacts: residual illumination “by default” on each pixel of the camera for the rows, and inhomogeneous contribution of each DMD mirror on the entire set of camera pixels for the columns.

ii) Because of the size of the speckle grains, two neighboring DMD mirrors may affect the material in the same way, as well, two pixels of the camera will be potentially correlated. To avoid this effect, we subsample the rows and columns of the matrix.

To draw the empirical spectral density, we then consider the following setup. We subsample the columns of the matrix up to $N = 200$ and leave the number of rows varying, more precisely $M = \gamma N$, with $\gamma \in \{1, ..., 6\}$. These sub-matrices thus constitute partitions of the estimated matrix, randomly picked 100 times to average the resulting densities. Figure 4 compares the experimental curves to the theoretical ones drawn according to the Marčenko-Pastur law. We see that the experiments qualitatively follow the predictions. We remark however that the larger $\gamma$ is, the more chances we have to consider the contributions of neighboring correlated pixels. This partly explains the increasing gap between both curves.

## 5. Focusing with the DMD

Knowing the TM gives a powerful and flexible tool to control light within the scattering medium [30]. In particular, it can be used to compute which DMD input has to be set, in order to display a given arbitrary pattern at the receiver end. In this section, we demonstrate the special case of
focusing light with maximum intensity on a desired pattern (a chosen sparse subset of the output pixels), with the TM measured experimentally as in the section 4. It should be emphasized that we keep the same experimental setup, with the binary DMD as input device. Here, simple inversion methods such as [30] cannot be used, as these require phase-modulated input and amplitude and phase detection.

We propose here to resort to the Bayesian variational approach [23], in a similar way as for the calibration but adapted here to the binary nature of the DMD inputs. The particularization of the algorithm to this case is given in appendix 7.2.

<figure><img src="image_37.png" /><figcaption>Fig. 5. Illustration of light focusing on 3 points. The circles mark the positions of the targets.</figcaption></figure>

We assess the performance of the proposed focusing approach through different experiments. The general setting is as follows. The DMD inputs, here taken of dimension $N = 1600$, are estimated from the desired outputs (see appendix 7.2), focusing on 1 to 4 target points. The procedure exploits the TM, reduced to its rows of interest and previously measured as discussed in section 4.

Figure 5 shows an example of the observed output field, corresponding to the estimated DMD configuration, optimized to focus on 3 points. To quantitatively evaluate the focusing performance, we measure the intensity enhancement factor, as:

$$ \eta \triangleq \frac{I_{\text{foc}}}{I_{\text{back}}}, \qquad (5) $$

where $I_{\text{foc}}$ is the intensity inside the target area after spatial binary amplitude modulation is performed, $I_{\text{back}}$ is the average background intensity. This value is measured for 100 trials, as a function of the number of calibration measurements used to learn the TM.

Two different setups are then considered: the single-point focusing case and the multi-target case.

## 5.1. Focusing on a single point

Figure 6 compares the enhancement factors achieved by two different focusing methods, namely a simple binary phase-conjugation - performing $\hat{x} = [\Re(D^H y) > 0]$ - and the proposed method, in the case where only one target point is focused. Results are presented under a “box” format, where:
<figure><img src="image_38.png" /><figcaption>Fig. 6. Single target experiment. Enhancement factor as a function of the number of measurements used to learn the TM (x-axis is $\alpha$, such that $p = \alpha N$ calibration measurements are used). For the same estimation of the TM, 2 focusing techniques are compared: binary phase conjugation (blue boxes), and the Bayesian technique [23] (red boxes).</figcaption></figure>

* the middle segment stands for the average enhancement $\bar{\eta}$ over the 100 trials,

* the upper and lower bounds of the rectangle define the interval $[\bar{\eta} - \sigma_{\eta} \quad \bar{\eta} + \sigma_{\eta}]$ (where $\sigma_{\eta}$ is the experimentally computed standard deviation), in which lies, under the Gaussian assumption, 68% of the trials,

* the whiskers represent the minimum and maximum values observed over the entire set of trials.

For each experiment point $\alpha \in \{2, \dots, 6\}$, such that $p = \alpha N$ calibration measurements are used to compute the TM, we display the boxes related to the binary phase-conjugation method (blue boxes), and the Bayesian technique (red boxes) [23] particularized in appendix 7.2.

As a first observation, we can see that the general dependency with regard to $\alpha$ noticeably resonates with the curve of the *prVBEM* algorithm in Fig. 3(b): there is a clear gap between the performance achieved for $\alpha = 2$ and for $\alpha = 3$, while, for $\alpha \ge 3$, the intensity enhancement keeps increasing but less significantly.

Interestingly, the *prVBEM* algorithm seems to outperform binary phase-conjugation, with regard to the mean and maximum values measured, but not in a statistically significant manner. Focusing on the most favorable case considered here, namely with $\alpha = 6$, the best intensity enhancement factor lies around 140, to be compared with the ideal expected enhancement given by $1 + \frac{1}{\pi}(\frac{N}{2} - 1) \approx 255$, see [17].

## 5.2. Focusing on multiple points

For this second setup, we are interested in the performance of the *prVBEM* approach in a context of multiple target points. Additionally to the intensity enhancement, we consider here the missed detection rate, defined as the number of trials (expressed in percentage) failing to focus on *at least* one of the multiple target points, i.e., the number of trials for which at least one of the *T* largest intensity peaks in the output image does not match any of the *T* targets.
<figure><img src="image_39.png" /><figcaption>Fig. 7. Multiple target experiment. (a) Average enhancement factor as a function of the number of measurements used to learn the TM (x-axis is $\alpha$, such that $p = \alpha N$ calibration measurements are used), and the number of target points (y-axis). (b) Missed detection rate (same axis as in (a)).</figcaption></figure>

Figure 7 represents these two figures of merit under diagram formats. They present an interesting general symmetry: increasing the number of targets or decreasing the number of calibration points leads to an increase of the missed detections and a decrease of the enhancement factor. The missed detection rate seems however more sensitive to the number of calibration points used to learn the TM: for $\alpha = 2$ and 2 target points, the algorithm fails with a rate approaching 40%, while for $\alpha = 3$ and the same number of targets, we keep a reasonable performance (around 10%). In a more general view, these figures greatly highlight the deep relation between the quality of the calibration and the focusing performance.

## 6. Conclusion

This paper shows that the full *complex-valued* transmission matrix of a strongly scattering material can be estimated, up to a global phase factor on each of its rows, with a simple experimental setup involving only *real-valued* inputs and outputs. In our experiment, the inputs are amplitude modulations on a *binary* DMD, and the output is the field intensity measured on a CCD camera, that gathers a significant amount of measurement noise. Note that no reference arm is used, that would allow interferometric measurements, but that would make the experimental setup more complex and considerably more unstable.

We here resort to Bayesian phase retrieval techniques, and we have shown that, amongst such techniques, a recently proposed variational approach (the so-called *prVBEM* algorithm [23]) allows a precise estimation of the transmission matrix, tractable in computational complexity and scalable for large-size signals, provided that we have a sufficiently large number of input-output calibration signals. Experimental results validate this concept, both in terms of output prediction, distribution of singular values, and in an application of light focusing onto a number of target points in the output plane. It should be emphasized that this estimation of the transmission matrix opens many applications beyond light focusing, may it be for imaging through the scattering material [30, 31], or for obtaining information about the scattering material itself.
## 7. Appendix: particularization of the prVBEM algorithm

In this appendix, we outline the main lines of the prVBEM algorithm particularized to its use for calibration of the TM and focalization. We refer the reader to the original paper [23] for a methodological presentation of the procedure.

### 7.1. prVBEM for the calibration of the TM

To circumvent the non-linearity induced by the absence of phase observations, the procedure introduces new variables modeling, on the one hand, the missing phases of the observations, and on the other hand, some acquisition noise. Thus, recalling that we resort to a conjugate-transposition of the matrix system (see Eq. (2)), each absolute-valued measurement $y_μ$, $μ \in \{1...P\}$, of any row y of Y, is expressed as

$$y_{\mu} = e^{j\theta_{\mu}} (\sum_{i=1}^{N} x_{\mu i} d_{i}^{*} + \omega_{\mu}), \qquad (6)$$

where $\theta_{\mu} \in [0, 2\pi)$ stands for its missing conjugate phase, $x_{\mu i}$ is the ith element of the μth row in X, $d_{i}^{*}$ corresponds to the ith conjugate element in the current estimated row d of D and $\omega_{\mu}$ is an additive noise, assumed centered isotropic Gaussian (denoted $CN$ in the following) with variance $\sigma^2$. We moreover suppose that the probability distributions for the entries of the matrix and for the missing phases are:

$$p(\mathbf{d}) = \prod_{i=1}^{N} p(d_i) \quad \text{with} \quad p(d_i) = \mathcal{CN}(0, \sigma_d^2), \qquad (7)$$

$$\text{and} \quad p(\theta) = \prod_{\mu=1}^{P} p(\theta_{\mu}) \quad \text{with} \quad p(\theta_{\mu}) = \frac{1}{2\pi}. \qquad (8)$$

Under these assumptions, the absence of phases in the observations is naturally taken into account in the model since marginalizing on $\theta_{\mu}$ leads to a distribution on $y_{\mu}$ which only depends on the moduli of $y_{\mu}$ and $\sum_{i=1}^{N} x_{\mu i} d_{i}^{*}$.

Within model (6)-(8), the recovery of the complex row d of the TM can be expressed as the solution of the following marginalized Maximum A Posteriori (MAP) estimation problem

$$\hat{\mathbf{d}} = \underset{\mathbf{d}}{\operatorname{argmax}} p(\mathbf{d}|\mathbf{y}), \qquad (9)$$

$$\text{with} \quad p(\mathbf{d}|\mathbf{y}) = \int_{\theta} p(\mathbf{d}, \theta|\mathbf{y}). \qquad (10)$$

Because of the marginalization on the hidden variables $\theta$, the direct computation of $p(\mathbf{d}|\mathbf{y})$ is however intractable in general. The prVBEM algorithm is based on the computation of a particular (the so-called “mean-field”) approximation $\hat{q}(\mathbf{d}, \theta) = \prod_i q(d_i) \prod_{\mu} q(\theta_{\mu})$ of the posterior joint distribution $p(\mathbf{d}, \theta|\mathbf{y})$. The procedure is then iterative, updating the factors $q(d_i)$ as

$$q(d_i) = \mathcal{CN}(m_i, \Sigma_i), \qquad (11)$$

with

$$m_i = \frac{\sigma_d^2}{\sigma_n^2 + \sigma_d^2 \mathbf{x}_i^H \mathbf{x}_i} \langle \mathbf{r}_i \rangle^H \mathbf{x}_i, \quad \Sigma_i = \frac{\sigma_n^2 \sigma_d^2}{\sigma_n^2 + \sigma_d^2 \mathbf{x}_i^H \mathbf{x}_i}, \quad \langle \mathbf{r}_i \rangle = \bar{\mathbf{y}} - \sum_{k \neq i} m_k^* \mathbf{x}_k, \qquad (12)$$
and

$$ \bar{\mathbf{y}} = \left[ y_{\mu} e^{(j \arg(y_{\mu}^* \langle z_{\mu} \rangle))} \frac{I_1(\frac{2}{\sigma_{\mu}^2} |y_{\mu}^* \langle z_{\mu} \rangle|)}{I_0(\frac{2}{\sigma_{\mu}^2} |y_{\mu}^* \langle z_{\mu} \rangle|)} \right]_{\mu=\{1...P\}} \quad (13) $$

$$ \langle z_{\mu} \rangle = \sum_{i=1}^{N} m_{i}^{*} x_{\mu i}. \quad (14) $$

In the equations above, $I_0$ (resp. $I_1$) stands for the modified Bessel function of the first kind for order 0 (resp. 1).

Coming back to problem (9), an approximation of $p(\mathbf{d}|\mathbf{y})$ simply follows from

$$ p(\mathbf{d}|\mathbf{y}) = \int_{\theta} p(\mathbf{d}, \theta|\mathbf{y}), \quad (15) $$

$$ \approx \int_{\theta} q(\mathbf{d}, \theta), \quad (16) $$

$$ \approx \int_{\theta} \prod_{i} q(d_{i}) \prod_{\mu} q(\theta_{\mu}), \quad (17) $$

$$ \approx \prod_{i} q(d_{i}), \quad (18) $$

resulting in an element-wise estimation of row $\hat{\mathbf{d}}$, such as $\hat{d}_i = m_i$.

## 7.2. prVBEM for focusing

The problem of focusing is expressed as an inverse problem, where, knowing the TM $\mathbf{D}$ and the observation $\mathbf{y}$, we look for the DMD input $\mathbf{x}$ such as described in (1). Adopting a similar modeling as in previous section, we then assume, for all elements $y_{\mu}$ with $\mu \in \{1, ..., M\}$,

$$ y_{\mu} = e^{j\theta_{\mu}} \left( \sum_{i=1}^{N} d_{\mu i} x_{i} + \omega_{\mu} \right), \quad (19) $$

where $\theta_{\mu} \in [0, 2\pi)$ stands for the missing conjugate phase, $d_{\mu i}$ is the $\mu$th element of the $i$th column in $\mathbf{D}$, $x_i \in \{0, 1\}$ corresponds to the state of the $i$th DMD pixel and $\omega_{\mu}$ is an additive noise, assumed centered isotropic Gaussian of variance $\sigma^2$. As previously, we suppose that the elements $\theta_{\mu}$ are independently and uniformly distributed in the interval $[0, 2\pi)$, however, in order to accommodate for binary inputs, we consider here a Bernoulli model for $\mathbf{x}$:

$$ p(\mathbf{x}) = \prod_{i=1}^{N} p(x_i) \quad \text{with} \quad p(x_i) = \text{Ber}(p_i) = \begin{cases} p_i & \text{if } x_i = 1, \\ 1 - p_i & \text{if } x_i = 0, \end{cases} \quad (20) $$

where we set $p_i$ to 0.5, noticing that asymptotically half of the DMD pixels are expected to be “ON” [17].

Then, within model (19)-(20), the recovery of the DMD inputs is expressed as the marginalized MAP estimation problem

$$ \hat{\mathbf{x}} = \underset{\mathbf{x}}{\operatorname{argmax}} p(\mathbf{x}|\mathbf{y}), \quad (21) $$

$$ \text{with} \quad p(\mathbf{x}|\mathbf{y}) = \int_{\theta} p(\mathbf{x}, \theta|\mathbf{y}). \quad (22) $$
The procedure is then similar to the one described previously for the calibration task. The prVBEM algorithm iteratively updates the factors of the mean-field approximation $q(\mathbf{x}, \theta) = \prod_i q(x_i) \prod_\mu q(\theta_\mu)$ of $p(\mathbf{x}, \theta|\mathbf{y})$, which leads, in the particular case of focusing, to the expression

$$q(x_i) = p(x_i) \exp \left( x_i \frac{2 \Re(\mathbf{d}_i^H \langle \mathbf{r}_i \rangle) - \mathbf{d}_i^H \mathbf{d}_i}{\sigma^2} \right), \quad (23)$$

where

$$\langle \mathbf{r}_i \rangle = \bar{\mathbf{y}} - \sum_{k \neq i} q(x_k = 1) \mathbf{d}_k, \quad (24)$$

$$\bar{\mathbf{y}} = \left[ y_{\mu} e^{(j \arg(y_{\mu} \langle z_{\mu} \rangle))} \frac{I_1(\frac{2}{\sigma^2} |y_{\mu}^* \langle z_{\mu} \rangle|)}{I_0(\frac{2}{\sigma^2} |y_{\mu}^* \langle z_{\mu} \rangle|)} \right]_{\mu=\{1...M\}} \quad (25)$$

$$\langle z_{\mu} \rangle = \sum_{i} q(x_i = 1) d_{\mu i}, \quad (26)$$

and $I_0$ (resp. $I_1$) stands for the modified Bessel function of the first kind for order 0 (resp. 1). Then, an approximation of $p(\mathbf{x}|\mathbf{y})$ is given by

$$p(\mathbf{x}|\mathbf{y}) = \int_{\theta} p(\mathbf{x}, \theta|\mathbf{y}), \quad (27)$$

$$\approx \int_{\theta} \prod_{i} q(x_i) \prod_{\mu} q(\theta_{\mu}), \quad (28)$$

$$= \prod_{i} q(x_i). \quad (29)$$

Using this approximation, the problem is easy to solve by a simple thresholding operation, i.e., $\hat{x}_i = 1$ if $q(x_i = 1) > 0.5$ and $\hat{x}_i = 0$ otherwise.

## Acknowledgements

AD is currently working at ENSTA Bretagne, LabSTICC (UMR 6285), 2 rue François Verny, F-29200 Brest, France. OK acknowledges the support of the Marie Curie Intra-European Fellowship for career development (IEF). LD acknowledges a joint research position with the Institut Universitaire de France. This work has been supported in part by the CSI:PSL grant and LABEX WIFI (Laboratory of Excellence within the French Program “Investments for the Future”) under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*, and by the ERC under the European Union’s 7th Framework Programme Grant Agreements 307087-SPARCS and 278025-COMEDIA.
# DOA ESTIMATION IN STRUCTURED PHASE-NOISY ENVIRONMENTS

Angélique Drémeau<sup>®*</sup>, Cédric Herzel<sup>*</sup>

<sup>®</sup> ENSTA Bretagne and Lab-STICC UMR 6285, Brest, F-29200, France

<sup>*</sup> INRIA Centre Rennes-Bretagne Atlantique, Campus universitaire de Beaulieu, F-35000 Rennes, France

## ABSTRACT

In this paper we focus on the problem of estimating the directions of arrival (DOA) of a set of incident plane waves. Unlike many previous works, which assume that the received observations are only affected by additive noise, we consider the setup where some phase noise also corrupts the data (as for example observed in atmospheric sound propagation or underwater acoustics). We propose a new methodology to solve this problem in a Bayesian framework by resorting to a variational mean-field approximation. Our simulation results illustrate the benefits of carefully accounting for the phase noise in the DOA estimation process.

**Index Terms**— Direction-of-arrival estimation, phase noise, variational Bayesian approximation, mean-field approximation

## 1. INTRODUCTION

Estimating the directions of arrival (DOA) of propagating plane waves is at the heart of many applicative domains including sonar, radar and mobile telecommunications. Among the rich literature dealing with this problem, the most popular method is probably conventional beamforming [1] which can roughly be interpreted as a least-square estimator. *Per se*, this approach suffers from a lower limit on the resolution achievable in the DOA estimation process (conditioned by the length of the sensor array). To overcome this issue, so-called “high-resolution” techniques, taking advantage of more prior information on (the number and the nature of) the sources, have been proposed in the literature (see e.g., [2, 3, 4, 5]). In [2], the authors introduced the well-known MUSIC algorithm which benefits from the knowledge of the number of the sources and rely on the assumption that the noise and the signal of interest live in perfectly separable subspaces. More recently, a “compressive” beamforming approach was proposed in [3, 4, 5], where a sparse prior was exploited to address the DOA estimation problem.

The contributions mentioned above assume that the incident plane waves are only corrupted by some additive noise. Unfortunately, when the waves travel through highly fluctuating media, as in the case of e.g., atmospheric sound

propagation [6] or underwater acoustics [7], this model does no longer describe accurately the physics underlying the propagation process. In such cases, a multiplicative phase noise typically corrupts the collected signal, making the corresponding DOA estimation problem quite challenging. We address this problem in the present paper.

Our approach is inspired from the recent standard “high-resolution” DOA methods [3, 4, 5] and some phase retrieval algorithms presented in [8, 9]. More specifically, we model the received signal as a sparse combination of elementary signals (taken from a redundant dictionary) and assume that the latter is corrupted by both additive and phase noise. Our methodology is grounded on a probabilistic Bayesian framework and relies on a variational mean-field approximation. In particular, we show how to nicely incorporate fine noise-phase models in this framework, extending in this respect, the approaches proposed in [8, 9].

## 2. PROBLEM STATEMENT

Our derivations are based on the following formulation of the DOA estimation problem: we consider an antenna composed of $N$ sensors and assume that the collected observation vector $y \in \mathbb{C}^N$ can be expressed as

$$y = PDz + \omega, \quad (1)$$

where $\omega \in \mathbb{C}^N$ and $P = \text{diag}(\{e^{j\theta_n}\}_{n=1}^N) \in \mathbb{C}^{N \times N}$ play respectively the role of an additive and a multiplicative phase noise. Matrix $D = [d_1 \dots d_M] \in \mathbb{C}^{N \times M}$ is made up of the steering vectors $d_i \triangleq [e^{j\frac{2\pi}{\lambda}\Delta\sin(\phi_i)} \dots e^{j\frac{2\pi}{\lambda}\Delta N \sin(\phi_i)}]^T$, where $\phi_i$'s are some potential angles of arrival, $\Delta$ is the distance between two adjacent sensors, and $\lambda$ is the wavelength of the propagation waves.

With this formulation, assuming that $y$ results from the combination of a few waves arriving from different angles $\phi_i$'s, the DOA estimation problem is basically equivalent to identifying the positions of the nonzero coefficients in $z$ (since each column of $D$ corresponds to a particular angle of arrival). In the standard DOA estimation problem, $P$ is assumed to be known with $P = I_N$, where $I_N$ is the $N \times N$ identity matrix. In this case, the model connecting the unknown vector $z$ to the measurements $y$ is linear; finding the

<footnote>*This work has been supported by the DGA/MRIS.</footnote>
position of the nonzero elements in z can then be carried out with standard sparse-representation algorithms, see e.g., [10].

In this paper, we consider the more complex case where $\mathbf{P} = \text{diag}(\{e^{j\theta_n}\}_{n=1}^N)$ and the $\theta_n$'s are unknown. More specifically, we assume that the phases $\theta_n$'s obey the following Markov model:

$$p(\boldsymbol{\theta}) = \prod_{n=2}^{N} p(\theta_n | \theta_{n-1}) p(\theta_1), \quad (2)$$

with $p(\theta_n|\theta_{n-1}) = \mathcal{N}(a \theta_{n-1}, \sigma_\theta^2)$, $\forall n \in \{2, \dots, M\}$, $a \in \mathbb{R}_+$, and $p(\theta_1) = \mathcal{N}(0, \sigma_1^2)$. From a practical point of view, this model allows us to describe spatial fluctuations of the propagation medium all along the antenna; the strength of the fluctuations is related to the value of parameter $\sigma_\theta^2$.

As noted in the introduction, assuming that $\mathbf{P}$ is unknown renders the DOA estimation much more difficult since it introduces uncertainties in the observation model. Before proceeding to the presentation of the proposed methodology to address this issue, we draw some connections with other applicable fields:

* Our work relates first to the phase retrieval problem (e.g., [11]) where the phase information of the observations is completely missing: only intensities or amplitudes are acquired. Formally, both problems share - explicitely or not - the same observation model (1) but differ in the prior distribution they enforce on the phases $\theta$, the absence of phase information being modeled by a non-informative prior, such as a uniform law (see [8, 9]).

* Then we note that, in top of being relevant for the DOA estimation problem with fluctuating media, this model is also of interest in the domain of digital communications where it can be used to characterize the transmission of complex modulation symbols over a channel affected by carrier phase noise, see e.g., [12].

### 3. BAYESIAN FORMULATION OF THE PROBLEM

We address the problem of estimating z from y when the realizations of $\omega$ and $\theta$ are unknown. To that end, we first place this problem into a Bayesian framework by defining suitable additional prior distributions on the unknown quantities.

To account for the sparsity of z, we suppose that, $\forall i \in \{1, \dots, M\}$, $z_i = x_i \cdot s_i$ with

$$p(x_i) = \mathcal{CN}(0, \sigma_x^2), \quad \text{and} \quad p(s_i) = \text{Ber}(p_i). \quad (3)$$

This so-called Bernoulli-Gaussian model has been now largely used in the literature to model sparse priors (see e.g., [13]). Next, we assume that $p(\omega)$ is a zero-mean Gaussian distribution with variance $\sigma^2$. This hypothesis is typically justified by the central limit theorem under the assumption

that the additive noise corrupting the data results from the aggregation of a large number of random parasitic contributions. Finally, the probabilistic model describing the behavior of $\theta$ is given by (2). As mentioned in the previous section, this simple model accounts for local variations of the propagation medium along the sensor network.

Based on this probabilistic model, we propose to look for the solution of the following Minimum Mean Square Error (MMSE) problem

$$\hat{\mathbf{z}} = \arg \min_{\tilde{\mathbf{z}}} \mathbf{E}_{\mathbf{z}|\mathbf{y}} \left[ \|\mathbf{z} - \tilde{\mathbf{z}}\|_2^2 \right], \quad (4)$$

relying on the marginal posterior distribution $p(\mathbf{z}|\mathbf{y}) = \int_{\theta} p(\mathbf{z}, \theta|\mathbf{y}) d\theta$. The computation of this marginalization being an intractable problem, we propose hereafter a practical procedure based on a mean-field approximation of the joint distribution $p(\mathbf{z}, \theta|\mathbf{y})$ to approach the solution of (4).

### 4. THE PROPOSED PROCEDURE

Mean-field approximations aim at approximating a posterior joint distribution by another one constrained to have a “simple” factorization while minimizing some distance with the targeted distribution. In the following, we will look for an approximation of $p(\mathbf{z}, \theta|\mathbf{y})$, say $\hat{q}(\mathbf{z}, \theta)$, obeying the following factorization $\hat{q}(\mathbf{z}, \theta) = \hat{q}(\theta) \prod_{i=1}^M \hat{q}(z_i)$. With this approximation of $p(\mathbf{z}, \theta|\mathbf{y})$, the evaluation of the marginal with respect to $\mathbf{z}$ is simplified since $\int_{\theta} \hat{q}(\mathbf{z}, \theta) = \prod_i \hat{q}(z_i)$ and the solution of the MMSE problem (4) can be approximated component-wise as

$$\hat{z}_i \approx \int_{z_i} z_i \hat{q}(z_i) dz_i.$$

A well-known algorithm to find a mean-field approximation of a target distribution is the so-called “Variational Bayes Expectation-Maximization” (VBEM) algorithm (see e.g., [14]). Particularized to our problem, this procedure searches for a local minimum of the following optimization problem:

$$\hat{q}(\mathbf{z}, \theta) = \arg \min_{q} \int_{\mathbf{z}, \theta} q(\mathbf{z}, \theta) \log \left( \frac{q(\mathbf{z}, \theta)}{p(\mathbf{z}, \theta|\mathbf{y})} \right) d\mathbf{z} \, d\theta \\ \text{subject to } q(\mathbf{z}, \theta) = q(\theta) \prod_{i=1}^{M} q(z_i),$$

by sequentially minimizing the cost function with respect to $q(\theta)$ and $q(z_i)$, $\forall i \in \{1, \dots, M\}$.

We show hereafter that the sequence of distributions generated by the VBEM algorithm, say $\{q^{(k)}(\theta), \{q^{(k)}(z_i)\}_i\}_k$, admit simple parametric expressions. We detail these expressions (at a given iteration of the procedure) in the following subsections. Due to space limitation, we do not provide the technical derivations but refer the reader to our companion report [15] for more details. For the sake of clarity, we also omit the iteration index $^{(k)}$ in the notations.
### 4.1. Update of $q(\theta)$

Under the condition<fnref n="1" /> of small $\sigma_{\theta}^2$, we can express the estimate of $q(\theta)$ as a Gaussian distribution:

$$q(\theta) = \mathcal{N}(\mathbf{m}_{\theta}, \Sigma_{\theta}),$$

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

$$\mathbf{m}_{\theta} = \Sigma_{\theta} \left( \operatorname{diag} \left( \frac{2}{\sigma^2} |\eta| \right) \operatorname{arg}(\eta) \right),$$

and $\Lambda_{\theta}^{-1}$ is the precision matrix attached to the prior distribution on $\theta$ (2), i.e.,

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

Vector $\eta$ is defined as $\eta \triangleq [\eta_1, \dots, \eta_N]^T$, with

$$\eta_n = y_n \sum_i \langle z_i \rangle^* d_{ni}^*,$$

where $d_{ni}^*$ is the conjugate of the nth element of $\mathbf{d}_i$, and $\langle z_i \rangle$ is defined using the current estimate of $q(z_i)$ as

$$\langle z_i \rangle \triangleq \int_{z_i} z_i q(z_i) dz_i.$$

Note that the distribution $q(\theta)$ being Gaussian, the marginals $q(\theta_n)$ come straightforwardly as

$$q(\theta_n) = \mathcal{N}(m_{\theta_n}, \Sigma_{\theta_n}),$$

where $m_{\theta_n}$ (resp. $\Sigma_{\theta_n}$) is the nth element in $\mathbf{m}_{\theta}$ (resp. in the diagonal of $\Sigma_{\theta}$). In practice, estimating these parameters can be efficiently implemented through a Kalman smoother [16] because of the particular structure of the precision matrix (5).

### 4.2. Update of $q(z_i)$

The iterates of the $q(z_i)$'s take the form of a mixture of two Gaussian distributions and only depend on the marginal distributions $q(\theta_n)$ of $q(\theta)$. More precisely, given that

$$q(z_i) = q(x_i, s_i) = q(x_i|s_i) q(s_i),$$

we obtain

$$q(x_i|s_i) = \mathcal{CN}(m_{x_i}(s_i), \Sigma_{x_i}(s_i)), \quad (6)$$

$$q(s_i) \propto \sqrt{\Sigma_{x_i}(s_i)} \exp \left( \frac{m_{x_i}(s_i)^* m_{x_i}(s_i)}{\Sigma_{x_i}(s_i)} \right) p(s_i), \quad (7)$$

where $\propto$ and $^H$ means respectively proportionality and trans- pose complex conjugate,

$$\Sigma_{x_i}(s_i) = \frac{\sigma^2 \sigma_x^2}{\sigma^2 + s_i \sigma_x^2 \mathbf{d}_i^H \mathbf{d}_i}, \quad (8)$$

$$m_{x_i}(s_i) = s_i \frac{\sigma_x^2}{\sigma^2 + s_i \sigma_x^2 \mathbf{d}_i^H \mathbf{d}_i} \mathbf{d}_i^H \langle \mathbf{r}_i \rangle, \quad (9)$$

$$\langle \mathbf{r}_i \rangle = \bar{\mathbf{y}} - \sum_{k \neq i} q(s_k = 1) m_{x_k}(s_k = 1) \mathbf{d}_k, \quad (10)$$

$$\bar{\mathbf{y}} = \left[ y_n e^{-j m_{\theta_n}} \frac{I_1(1/\Sigma_{\theta_n})}{I_0(1/\Sigma_{\theta_n})} \right]_{n=\{1...M\}}, \quad (11)$$

and $I_0$ (resp. $I_1$) stands for the modified Bessel of the first kind of order 0 (resp. 1).

We note that the update equation (6)-(11) share some connections with the phase retrieval algorithm presented in [9]. The latter, relying also on a VBEM algorithm, differs from the proposed procedure in the definition of the prior distributions (2) and (3), respectively replaced by a uniform and a Gaussian distributions. In practice, both procedures share a similar structure. Leaving out the choice made here of a sparse-enforcing prior on $\mathbf{z}$, the main difference lies in the “reconstructed” phases $m_{\theta_n}$ in (11): while their definition relies here on the parameters of the Markov chain through the precision matrix (5), they only depend on the observations in [9] where a non-informative prior is considered.

### 4.3. Noise estimation

As emphasized in [13], the estimation of model parameters can easily be embedded within the VBEM procedure. Among them, the noise variance is of particular interest. Measure of the (mean) discrepancies between the observations and the assumed model, its iterative estimation usually helps the convergence of the algorithm to a proper local minimum, as observed in [13]. Particularized to model (1)-(3), this leads to

$$\begin{aligned} \hat{\sigma}^2 &= \arg \max_{\sigma^2} \int_{\mathbf{z}} q(\mathbf{z}, \theta) \log p(\mathbf{z}, \theta; \mathbf{y}) \sigma^2) \, d\mathbf{z} \, d\theta, \\ &= \frac{1}{N} \left( \mathbf{y}^H \mathbf{y} - 2 \sum_i \Re \left[ q(s_i=1) m_{x_i}(s_i=1) \bar{\mathbf{y}}^H \mathbf{d}_i \right] \right. \\ &\quad + \sum_i \sum_{k \neq i} q(s_i=1) q(s_k=1) m_{x_i}(s_i=1)^* m_{x_k}(s_k=1) \mathbf{d}_i^H \mathbf{d}_k \\ &\quad \left. + \sum_i q(s_i=1) \left( \Sigma_{x_i}(s_i=1) + |m_{x_i}(s_i=1)|^2 \right) \mathbf{d}_i^H \mathbf{d}_i \right). \end{aligned}$$

In the following, we will refer to the proposed procedure as “paVBEM” for “phase-aware VBEM algorithm”.

## 5. EXPERIMENTS

In this section, we assess numerically the effectiveness of the proposed approach. We consider the problem of identifying

<footnote n="1">¹We justify this condition in the technical report [15].</footnote>
<figure><img src="image_40.png" /><figcaption>**Fig. 1.** Evolution of the (averaged) normalized correlation as a function of the variance $\sigma^2$ when $K = 2$.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Fig. 2.** Evolution of the (averaged) normalized correlation as a function of the variance $\sigma^2$ when $K = 5$.</figcaption></figure>

the directions of arrival of $K$ plane waves from $N = 256$ observations ($K$ will be specified later on). 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], \forall k \in \{1, \dots, K\}$. The set of angles $\{\phi_i = -\pi + i \frac{\pi}{50}\}_{i \in \{1, \dots, 50\}}$ together with the choice of the parameter $\Delta/\lambda = 4$ define the columns of the dictionary $\mathbf{D}$ (see section 2). We set the following parameters for the phase Markov model (2): $\sigma_1^2 = 10^6, \sigma_\theta^2 = 1$ and $a = 0.8$. This corresponds to the situation where one has a large uncertainty on the initial value of the phase noise but connections exist between the phase noise on adjacent sensors.

As well as for the classic phase retrieval problem, the solution of (4) can only be found up to a global phase. Hence, computing the mean square distance between the ground truth $\mathbf{z}$ and its reconstruction $\hat{\mathbf{z}}$ does not constitute a relevant figure of merit. Instead, we consider their normalized correlation, $\frac{|\mathbf{z}^H \hat{\mathbf{z}}|}{\|\mathbf{z}\|_2 \|\hat{\mathbf{z}}\|_2}$. This quantity is averaged over 50 realizations for each point of simulation. We compare the performance of the following algorithms: i) the conventional beamforming introduced in [1] (⊼); ii) the so-called *prVBEM* algorithm proposed in [9] as a solution to the phase retrieval problem (⊼); iii) the *paVBEM* procedure described in section 4 (⊼); iv) a relaxed version of *paVBEM* in which the sparsity of $\mathbf{z}$ is not exploited but replaced by a Gaussian prior (⊼).

The performance of these procedures are illustrated in Fig. 1 and 2 as a function of the noise variance $\sigma^2$ for $K = 2$ and $K = 5$, respectively. We see that the beamforming algorithm, which was originally proposed to solve the DOA estimation problem in the standard linear setup ($\mathbf{P} = \mathbf{I}_N$), fails to cope with the presence of fluctuations in the phase $\theta$. The three other algorithms achieve different levels of performance, depending on the power of the additive noise and the number of incident waves. We note that all these procedures are derived from a similar optimization procedure

but consider different degrees of knowledge on $\theta$ and $\mathbf{z}$. In [9], the authors assume that the phase is uniformly distributed and the sparse nature of $\mathbf{z}$ is ignored; in the relaxed version of the proposed procedure (see [15]) the phase model (2) is exploited but the sparsity of $\mathbf{z}$ is not taken into account; finally, as explained previously, the methodology presented in this paper integrates both the phase model (2) and the sparsity of $\mathbf{z}$ in the estimation process. We see from Fig. 1 and 2 that the performance of these algorithms directly relates to the level of information they exploit: the proposed methodology outperforms its relaxed counterpart which, in turn, leads to better performance than the procedure proposed in [9].

We also notice that the procedures achieve better performance when $K = 5$ than $K = 2$. This counter-intuitive behavior is typical for phase-retrieval problems. In fact, it is easy to see that, when no phase information is available (i.e., $\theta_n$ is uniform on $[0, 2\pi]$), only the modulus of $y_n$ provides some information on $\mathbf{z}$. In such a case, the worst situation occurs when $K = 1$ since all the elements of $\mathbf{D}$ has the same modulus (equal to 1) and the observations thus provide no information on $\mathbf{z}$. The setup $K = 2$ is close to this worst case, hence explaining the observed behavior.

## 6. CONCLUSION

We have presented a novel algorithm able to estimate DOA in environments corrupted by phase noise. Our approach relies on a mean-field approximation and exploits two types of priors: on the DOA through a sparse-enforcing distribution and on the phase noise through a Markov model. Our experiments have confronted the proposed approach to conventional beamforming and similar variational approaches handicapped by non-informative priors. In this regard, its good performance tends to prove a successful inclusion of the priors. Future work will include further assessment in underwater acoustics.
7. REFERENCES

[1] D. H. Johnson and D. E. Dudgeon, *Array signal processing: concepts and techniques*, Englewood Cliffs, NJ, 1993.

[2] R. O. Schmidt, "Multiple emitter location and signal parameter estimation," *IEEE Trans. Antennas Propagation*, vol. AP-34, pp. 276-280, 1986.

[3] W. Mantzel, J. Romberg, and K. Sabra, "Compressive matched-field processing," *Journal of Acoustical Society of America*, vol. 132, no. 1, pp. 90-102, 2012.

[4] A. Xenaki, P. Gerstoft, and K. Mosegaard, "Compressive beamforming," *Journal of Acoustical Society of America*, vol. 136, no. 1, pp. 260-271, 2014.

[5] S. Fortunati, R. Grasso, F. Gini, M. S. Greco, and K. Lepage, "Single snapshot DOA estimation by using compressed sensing," *EURASIP Journal on Advances in Signal Processing*, vol. 120, 2014.

[6] S. Cheinet, L. Ehrhardt, D. Juvé, and P. Blanc-Benon, "Unified modeling of turbulence effects on sound propagation," *Journal of Acoustical Society of America*, vol. 132, no. 4, pp. 2198-2209, 2012.

[7] R. Dashen, S.M. Flatté, W.H. Munk, K.M. Watson, and F. Zachariasen, *Sound Transmission Through a Fluctuating Ocean*, Cambridge Monographs on Mechanics. Cambridge University Press, 2010.

[8] P. Schniter and S. Rangan, "Compressive phase retrieval via generalized approximate message passing," in *Communications, Control, and Computing (Allerton)*, October 2012.

[9] A. Drémeau and F. Krzakala, "Phase recovery from a bayesian point of view: the variational approach," in *Proc. IEEE Int'l Conference on Acoustics, Speech and Signal Processing (ICASSP)*, Brisbane, Australia, April 2015, pp. 3661-3665.

[10] S. Foucart and H. Rauhut, *A mathematical introduction to compressive sensing.*, Applied and Numerical Harmonic Analysis. Birkhäuser, 2013.

[11] E. J. Candès, T. Strohmer, and V. Voroninski, "Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming," *Communications on Pure and Applied Mathematics*, vol. 66, no. 8, pp. 1241-1274, 2013.

[12] G. Colavolpe, A. Barbieri, and G. Caire, "Algorithms for iterative decoding in the presence of strong phase noise," *IEEE Journal on selected areas in communications*, vol. 23, no. 9, pp. 1748-1757, 2005.

[13] A. Drémeau, C. Herzet, and L. Daudet, "Boltzmann machine and mean-field approximation for structured sparse decompositions," *IEEE Trans. On Signal Processing*, vol. 60, no. 7, pp. 3425-3438, July 2012.

[14] M. J. Beal and Z. Ghahramani, "The variational bayesian EM algorithm for incomplete data: with application to scoring graphical model structures," *Bayesian Statistics*, vol. 7, pp. 453-463, 2003.

[15] A. Drémeau and C. Herzet, "DOA estimation in structured phase-noisy environments: technical report," Tech. Rep., available at http://arxiv.org:1609.03503, September 2016.

[16] J. M. Mendel, *Lessons in Estimation Theory for Signal Processing Communications and Control*, Prentice Hall Signal Processing Series, Englewood Cliffs, NJ, 1995.