Jérôme Antoni
email: jerome.antoni@insa-lyon.fr

Thibaut Le Magueresse

Quentin Leclere

Patrice Simard

Quentin Leclère

Sparse acoustical holography from iterated Bayesian focusing

Keywords: Acoustical holography, acoustic imaging, source identification, inverse problem, sparsity AF Aperture Function BIBF Broadband Iterated Bayesian Focusing BF

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The reconstruction of sound sources from acoustic array measurements is known as a difficult inverse problem. Several methods have been proposed over the years with sustained efforts to reach improved performance. An unavoidable limitation arises from the fact that reconstructing a continuous acoustical field from a limited number of discrete and remote measurements, as returned by an array of microphones, has no unique solution in general. This issue is particularly critical in near-field configurations, where the spectrum of the propagation operator is characterized by a fast exponential decay and is therefore hardly invertible. Consequently, sound sources cannot be reliably reconstructed unless additional information is incorporated into the inverse problem in order to reduce the space of possible solutions. So far, Tikhonov regularization has been the prevalent method to address this issue, by selecting solutions with minimal energy. However, this is often at the expense of accuracy of the reconstruction, and especially in term of spatial resolution. When the sources of interest are known to have a sparse representation, much better strategies than Tikhonov regularization are actually feasible.

Compressed sensing and the related concept of sparsity have recently encountered a huge success for solving illposed problems in many scientific fields [START_REF] Starck | An overview of inverse problem regularization using sparsity[END_REF]. Compressed sensing is essentially a new sampling theory which provides a valid alternative to Shannon theory in situations where samples are forced to remain limited. It relies on the premise that there exists a sparse representation of the quantities of interest and furnishes unexpected optimistic bounds on the number of measurements required to reconstruct it without loss. When applied to inverse problems that seek a sparse solution, sparsity intervenes as an alternative to Tikhonov regularization by imposing an apriori on the shape of the expected solution and not only on its energy. This is apt to considerably narrow the space of possible solutions and, in some situations, even shrink it to a unique solution. Technically speaking, the difference is in replacing the 2 -norm penalty used in Tikhonov regularization by a p -norm, with p a power strictly less than 2, typically equal to 1. Whereas Tikhonov regularization tries to minimize the energy (reflected by the 2 -norm) of the solution with the side-effect of spreading the source in space, "sparse regularization" controls the source magnitude by forcing it to take only a few high values and shrinking it towards zero everywhere else.

The above benefits have been recognized early for the localizing of sources, a subject which has nourished an abundant literature in array processing (beamforming, Direction Of Arrival) [START_REF] Fuchs | Linear programming in spectral estimation. application to array processing[END_REF][START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. In these seminal works, sparsity arrives from the assumption that the acoustical field is produced by a few point sources (monopoles) located in the farfield. Remarkable results have been obtained based on minimizing the 1 -norm, which proved superior to those of the state-of-the-art, such as returned by classical high-resolution (parametric) methods (e.g. Capon, MUSIC). Numerous works have followed, which constitute a body of literature known as compressive beamforming [START_REF] Gurbuz | A compressive beamforming method[END_REF]. Applications to acoustics are found in Refs. [START_REF] Edelmann | Beamforming using compressive sensing[END_REF][START_REF] Simard | Acoustic source identification: Experimenting the 1 minimization approach[END_REF][START_REF] Zhong | Compressive sensing beamforming based on covariance for acoustic imaging with noisy measurements[END_REF][START_REF] Xenaki | Compressive beamforming[END_REF][START_REF] Gerstoft | Multiple and single snapshot compressive beamforming[END_REF][START_REF] Xenaki | Grid-free compressive beamforming[END_REF].

Incidentally, compressive beamforming shares some similarities with deconvolution algorithms -originally derived as a post-processing of beamforming results -such as DAMAS [START_REF] Brooks | A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays[END_REF], DAMAS2, SC-DAMAS [START_REF] Yardibi | Sparsity constrained deconvolution approaches for acoustic source mapping[END_REF], CLEAN-SC 2 [START_REF] Sijtsma | CLEAN Based on Spatial Source Coherence[END_REF], NNLS, since they all try to return a final beamforming map made of a few point sources [START_REF] Chu | Comparison of deconvolution methods for the visualization of acoustic sources based on cross-spectral imaging function beamforming[END_REF]. Another connection can be made with covariance matching methods [START_REF] Yardibi | A covariance fitting approach for correlated acoustic source mapping[END_REF], which also enforce sparse solutions as pointed in [START_REF] Yang | A discretization-free sparse and parametric approach for linear array signal processing[END_REF]. This fact has been exploited in Refs. [START_REF] Chu | Robust Bayesian super-resolution approach via sparsity enforcing a priori for near-field aeroacoustic source imaging[END_REF][START_REF] Chu | A robust super-resolution approach with sparsity constraint in acoustic imaging[END_REF] which proposed super-resolution algorithms derived in a Bayesian setting to estimate the power of incoherent point sources in the presence of background noise. In addition to refining spatial resolution, deconvolution and covariance matching methods have been shown to improve the estimation of source powers over a large interval of frequencies; this pushes still further one of the recurrent limits of acoustic imaging techniques. However, this family of techniques critically relies on the assumption of uncorrelated sources, a constraint which is actually relaxed in compressive beamforming [START_REF] Xenaki | Compressive beamforming[END_REF].

Generalized versions of beamforming with sparsity constraints have also been investigated in acoustics independently from the recent literature on compressive beamforming. The so-called Generalized inverse beamforming (GIBF) is truly an inverse problem that tries to recover a limited number of monopoles and dipoles by imposing a 1 -norm penalty [START_REF] Suzuki | L1 generalized inverse beam-forming algorithm resolving coherent/incoherent, distributed and multipole sources[END_REF][START_REF] Zavala | Generalized inverse beamforming with optimized regularization strategy[END_REF]. Strongly inspired by deconvolution algorithms, it has been shown to outperform them in several instances, in particular in the presence of correlated sources. More general configurations have recently been investigated where the acoustical field is parsimoniously represented by a few elements taken in dictionaries different from a set of monopoles or dipoles. Reference [START_REF] Peillot | Acoustic sources joint localization and characterization using compressive sampling to cite this version : Acoustic sources joint localization and characterization using compressive sampling[END_REF] proposes to use a dictionary made of spherical harmonics so as to jointly localize the sources and identify their directivity.

Although compressive beamforming (broadly speaking as including all the aforementioned techniques) has probably reached its maturity, it has so far mostly relied on the far-field approximation of the acoustic propagator. This is an explicit assumption in Refs. [START_REF] Edelmann | Beamforming using compressive sensing[END_REF][START_REF] Simard | Acoustic source identification: Experimenting the 1 minimization approach[END_REF][START_REF] Zhong | Compressive sensing beamforming based on covariance for acoustic imaging with noisy measurements[END_REF][START_REF] Xenaki | Compressive beamforming[END_REF][START_REF] Gerstoft | Multiple and single snapshot compressive beamforming[END_REF][START_REF] Xenaki | Grid-free compressive beamforming[END_REF] concerned with DOA and is often implicit in deconvolution approaches based on beamforming. While GIBF is theoretically not subjected to this constraint, it seems that it has been mostly validated in the far-field as well.

There are few reports to date of the application of sparse inversion to (near-field) acoustical holography (NAH). Reference [START_REF] Chardon | Near-field acoustic holography using sparse regularization and compressive sampling principles[END_REF] was probably first to address the subject. The objective therein was to demonstrate that a substantial reduction in the number of microphones can be achieved (by more than a factor 10) without affecting the reconstruction performance by exploiting the sparsity of the acoustical field. The two fundamental assumptions are that the acoustical field can be represented by a few plane waves (sparsity in the wavenumber domain) and that the microphones are distributed randomly.

Steered by somewhat different objectives, Ref. [START_REF] Hald | Fast wideband acoustical holography[END_REF] recently introduced the concept of wideband acoustical holography (WBH). The motivation in this work is for a method that can bridge the gap between NAH and classical beamforming in order to reconstruct sources over a wide frequency range. On the one hand, NAH can cover the low frequencies but requires the array to be placed at a small distance. On the other hand, beamforming has good performance in a high frequency range but requires the array to be moved at a larger distance. By enforcing sparsity in the Equivalent Source Method (ESM) (i.e. representation by a few point sources), the proposed WBH method is able to reconstruct sources with a remarkable quantification of the sound power over a wide frequency range. A numerical experiment seems to indicate that WBH is also able to consider extended sources such as produced by the vibration of a plate. The algorithm proposed in Ref. [START_REF] Hald | Fast wideband acoustical holography[END_REF] is named "fast" because it uses an Iterative Hard Thresholding algorithm instead of a direct 1 -norm minimization, yet it still depends on several user-dependent parameters.

Following a similar same idea, Ref. [START_REF] Fernandez-Grande | A sparse equivalent source method for near-field acoustic holography[END_REF] also proposed a sparse version of ESM, yet based on direct 1 -norm minimization. In accordance with the results of Ref. [START_REF] Chardon | Near-field acoustic holography using sparse regularization and compressive sampling principles[END_REF], the authors also stress the importance of using a sampling scheme that leads to low column coherence of the sensing matrix, a condition which is better met when the distance between the measurement plane and the source plane is decreased. As in Ref. [START_REF] Hald | Fast wideband acoustical holography[END_REF], the applicability to the reconstruction of extended sources is also demonstrated by placing the equivalent sparse sources behind the actual radiating surface. The investigated frequency range is limited to that of classical NAH, where moderate improvement is observed in the reconstruction of source spatial distribution as compared to 2 regularization. Reference [START_REF] Abusag | Near-field acoustic holography using the method of superposition[END_REF] introduced independently a sparse wave superposition method (with sparsity placed on the so-called "charge points" which actually coincide with the monopoles generating the acoustic field), yet its validation was limited to a numerical simulation. Most recently, the spare ESM has been extended to the case where a limited number of acoustical modes rather than monopoles are recovered in order to explain a source field; in Ref. [START_REF] Bi | Sound field reconstruction using compressed modal equivalent point source method[END_REF] the "radiation" modes are obtained from the eigenvalue decomposition of the power resistance matrix [START_REF] Yong-Fa | Acoustic field reconstruction using source strength density acoustic radiation modes[END_REF], whereas in Ref. [START_REF] Canclini | Dictionary-based equivalent source method for near-field acoustic holography[END_REF] they are computed from a finite element model of the radiating object.

It is noteworthy that Ref. [START_REF] Pereira | Acoustic imaging in closed spaces[END_REF] introduced a sparse ESM before the former publications, which was later named Iterative ESM (iESM). It was compared with GIBF in Ref. [START_REF] Oudompheng | A theoretical and experimental comparison of the iterative equivalent source method and the generalized inverse beamforming[END_REF] and also with deconvolution approaches in Ref. [START_REF] Leclère | A unified formalism for acoustic imaging based on microphone array measurements[END_REF].

Two remarks are worthwhile in light of the above literature review on compressive beamforming and sparse NAH. First, most of the aforementioned works enforce sparsity through an 1 -norm penalty and essentially differ in the optimization algorithm used to minimize it. Although few methods set up the inverse problem with a general pnorm, p ≤ 2, the inverse problem is often solved with p = 1 in order to allow the use of a convex optimization algorithm (e.g Iterated Reweighted Least Squares (IRLS) in [START_REF] Suzuki | L1 generalized inverse beam-forming algorithm resolving coherent/incoherent, distributed and multipole sources[END_REF] and linear programming in [START_REF] Fernandez-Grande | A sparse equivalent source method for near-field acoustic holography[END_REF]). One exception is Ref. [START_REF] Leclère | A unified formalism for acoustic imaging based on microphone array measurements[END_REF] which shows the advantage of using p < 1 and in particular p → 0 by inserting a thresholding step in IRLS. Another exception is WBH which introduces an optimization algorithm similar to Iterative Hard Thresholding, yet without specifying which norm is actually minimized.

Second, a challenging aspect of all the above approaches is to involve a regularization parameter (sometimes implicitly) that balances the effect of the sparsity enforcing penalty with the data fidelity term. Although setting the regularization parameter is known to be extremely critical, especially in NAH, no systematic rule seems to have found a consensus in the sparse case, contrary to state-of-the-art on 2 -norm Tikhonov regularization [START_REF] Hansen | Rank-deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion[END_REF]. The "crossvalidation" criterion is often advocated, yet it is costly to apply since it has no "generalized" equivalent as can be for Tikhonov regularization [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF].

The main contribution of the present paper is to demonstrate that sparse NAH is naturally devised by iterating the Bayesian Focusing (BF) method introduced in Ref. [START_REF] Antoni | A bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF]. BF solves the inverse acoustical problem by considering the sound sources as random variables. They are first vaguely described by a prior probability density function (PDF) which reflects the user expectation before the experiment is run. Next, the posterior PDF is found as the product of the likelihood function and the prior PDF, which provides the probability of the sound sources after taking the measurements. A typical point estimate is then returned by the maximum value of the posterior PDF -the maximum aposteriori (MAP) -i.e. the source configuration with highest probability of occurrence given the observed data. In BF, the variance of the prior PDF -coined the aperture function (AF) -appears as a key quantity; intuitively, it will take high values where the sound sources are expected to radiate from and nil values elsewhere. Therefore, by shrinking the space of solutions to a confined region, the AF improves the quality of the reconstruction in terms of spatial resolution and of source strength quantification over a larger frequency interval than allowed by conventional methods. This is all the more remarkable as the AF is allowed to become very narrow, as in the case of sparse sources. By analogy with optics, it plays the role of a lens that focuses the light onto a point.

In Ref. [START_REF] Antoni | A bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF], the AF was manually tuned by the user. A natural idea is to automate this process by using the current estimate of the source distribution as the AF to apply in the next iteration. It is shown in this paper that iterating BF generally leads to the recovery of sparse solutions. Not only is the common 1 -norm penalty recovered as a particular case, but many other sparsity enforcements -possibly stronger than the 1 -norm -can be devised depending on how the AF is updated. The reason is that when the AF is considered as a random quantity (in the Bayesian setting all unknowns are described by random variables), it then assigns the sound sources with a prior PDF in the form of a "scale mixture of Gaussians" (SMoG), which necessarily promotes sparsity. Based on this finding , the convergence of the iterations can be proved by establishing a formal equivalence with the Expectation-Maximization algorithm. This algorithm will be referred as Iterating Bayesian Focusing (IBF), which is reminiscent of the general concept of "sparse Bayesian learning".

One advantage of IBF is to provide some physical insight into the mechanism of sparsity enforcement. Through the role played by the AF, it explains why promoting sparse solutions not only increases the spatial resolution, but also improves the estimation of the source levels and the source directivity over a larger frequency range.

Another prime advantage of the approach is to inherit the automated regularization of BF [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF]. Contrary to Tikhonov regularization for which well established algorithms exist (e.g. the L-curve, GCV), this point remains an issue in compressive beamforming. Here, by taking a Bayesian perspective, the regularization parameter is jointly inferred in IBF with proved convergence.

The Bayesian framework also easily allows the processing of multiple snapshots and multiple frequencies, which both involve the consideration of group sparsity [START_REF] Huang | The benefit of group sparsity[END_REF]. The multi-snapshot case is addressed here by means of the cross-spectral matrix (CSM), which turns out a sufficient statistic (the approach differs in this aspect from that of Ref. [START_REF] Gerstoft | Multiple and single snapshot compressive beamforming[END_REF]). The advantage of using the CSM is purely practical, since it is often the only quantity recorded by commercial data acquisition systems; indeed, it happens to be sufficient for running most imaging methods while allowing at the same time a substantial reduction of the amount of stored data. The multi-frequency approach is of interest when the sound sources have a constant position in space, independently of their frequency content, thus leading to broadband processing. Both the multi-snapshot and the multi-frequency considerations are shown to improve the estimation of the sound sources.

The rest of the manuscript is organized as follows. Section 2 starts with a review of BF and of the key concept of an AF. Next, Section 3 explains the focusing mechanism produced by the AF and the resulting improvement in the estimation of the source strength over a wide frequency range. Section 4 then introduces the principle of the IBF method and its formalization by means of prior PDFs in the form of an SMoG. Specific cases are detailed relating to the multivariate complex Student-t, K-distribution, Laplacian, Generalized Gaussian and mixture of Gaussians. Section 5 compares these different possibilities in terms of the number of parameters to set and of systematic and random errors. Section 6 extends IBF to the broadband case. Finally, section 7 illustrates the application of IBF to experimental data. In order to ease the reading of the manuscript, all technical parts are moved to the Appendix.

The following notations are used throughout the document. By convention, vectors are denoted by bold lower-case letters and matrices by bold upper-case letters (e.g. u and U, respectively). I stands for the identity matrix. Symbol * stands for the complex conjugate of a complex variable. Symbol H stands for the transpose conjugate of a vector or a matrix. Diag(a q ) is the diagonal matrix with element a q in its q-th diagonal entry. Operator E{X} denotes the expected value of a random variable X and E{X|Y} the expected value of X conditioned to random variable Y. By convention, square brackets are used to denote a PDF (e.g. [X] denotes the PDF of random variable X). Notation [X|Y] stands for the PDF of X conditioned on Y. N C (x; µ, Ω) means that vector x follows a complex Gaussian with mean µ and covariance matrix Ω. The collection of indexed vectors

{p i } N i=1 , p i ∈ C M (resp. {s i } N i=1
) is denoted as P (resp. S), an M × N matrix whose i-th column is p i . The p -norm of a vector x in C N is defined as ||x|| p . = ( N q=1 |x q | p ) 1/p . The position vector in the continuous space is denoted by r. A function of space (e.g. σ) may be noted in three different forms: i) the functional notation, σ(r), where it is viewed as a continuous function of the space variable r, ii) the discrete notation, σ q . = σ(r q ), where it is viewed as a sample at position r q , and iii) the matrix notation, Σ . = Diag(σ q ), which collects all the space samples in a diagonal matrix.

Review of Bayesian Focusing (BF)

Parametrization of the problem

The array is composed of M microphones indexed by the lower-case letter m, located at positions r m , m = 1, ..., M. All measurements are considered at a given frequency f , after application of the Fourier transform on a series of N snapshots indexed by the lower-case letter i (i = 1, ..., N). The Fourier coefficients of the sound pressure measured by microphone m and assigned to snapshot i then reads p m,i , where explicit dependence on frequency f is dropped to simplify the notation. The M Fourier coefficients p m,i , m = 1, ..., M are collected in the column vector p i . The cross-spectral matrix (CSM) of the measurements averaged over N snapshots is defined as

S pp = 1 N N i=1 p i p H i . (1)
In the following, formulations in either continuous or discrete space will be used alternatively, the former because it provides deeper physical insight and the latter because it corresponds to the numerical resolution of the problem.

In the continuous formulation, the source distribution is assumed to be a scalar field -e.g. parietal pressure, normal velocity -that adheres on the source surface Γ s . In a functional notation, it is noted s i (r) where r ∈ Γ s stands for the position vector. In the discrete formulation of the problem, the spatial samples of the source distribution taken at positions r q , q = 1, ..., Q are noted s q,i . = s i (r q ). They are referred to as the "source coefficients" and are collected in the column vector s i , with dimension Q. The power spectrum of the source coefficients at position r q , averaged over N snapshots, is defined as

S qq = 1 N N i=1 |s q,i | 2 .
(2)

Direct problem

The relationship between the source distribution and the radiated pressures measured by the microphones is hereafter considered in the form of a Fredholm integral. In a functional notation, the direct problem therefore reads

p m,i = Γ s G(r m |r)s i (r)dΓ(r) + n m,i , i = 1, ..., N (3)
where G(r m |r) is the Green function between the source point r ∈ Γ s and the microphone location r m subjected to Neumann boundary conditions on Γ s . The quantity n m,i stands for additive noise at microphone m. In anticipation of a numerical resolution of the inverse problem, the discrete version of Eq. ( 3) is considered,

p i = Gs i + n i , i = 1, ..., N, (4)
where G stands for an M × Q matrix whose entry (i, q) is fed with G(r m |r q )∆Γ(r q ) -with ∆Γ(r q ) a small surface element at position r q -and where vector n i stacks the noise terms n m,i .

The aim of the inverse problem is to recover an estimate of the source distribution s i (r) from the observation of the measured pressures P . = {p i } N i=1 returned by the array of microphone. BF provides a solution to this problem that makes intensive use of prior information about the spatial structure of the source field to be recovered.

The Bayesian program

The philosophy in the Bayesian approach is to see all unknowns in the inverse problem as random variables and to infer them through their PDFs. Of interest here is the posterior PDF [s i |p i ] of the vector s i of coefficients given the measurement p i . The knowledge of [s i |p i ] completely characterizes the information that can be gained on s i once the data have been measured. A popular point estimate used in this paper is the MAP estimate, ŝi = Argmax

s i [s i |p i ], (5)
to be interpreted as the value of s i with maximum probability, yet another plausible estimate is the posterior mean

E{s i |p i } = s i [s i |p i ]ds i .
Bayes's rule then makes it possible to express the "inverse probability" [s i |p i ] in terms of the probabilities assigned to the direct problem, that is

[s i |p i ] = [p i |s i ][s i ] [p i ] (6)
where [p i |s i ] is the PDF of observing the data given an instance of s i (the so-called likelihood function), [s i ] is the PDF of the values possibly taken by s i before the data are observed (the so-called prior) and [p i ] = [p i |s i ][s i ]ds i is the "evidence". The Bayesian program is to introduce the PDFs [p i |s i ] and [s i ] according to the specificity of the problem and to deduce the posterior PDF [s i |p i ] by using Eq. ( 6). Since the data are processed in the Fourier domain, it holds from the Central Limit Theorem applied to the Fourier transform that the elements of the noise vector n i rapidly converge to a complex Gaussian [START_REF] Brillinger | Time Series: Data Analysis and Theory[END_REF]. Assuming for simplicity a zero mean and a covariance matrix proportional to the identity (i.e. a homogeneous field),

[p i |s i ] = e -β -2 (p i -Gs i ) H (p i -Gs i ) β 2M π M = N C (p i ; Gs i , β 2 I) (7)
where β 2 stands for the value of the noise power spectrum at frequency f (independently of space according to the previous assumption) and π = 3.14159... is Archimedes' constant.

The specification of the prior PDF [s i ] allows much more flexibility and is actually central in this paper. The assumption initially used in Ref. [START_REF] Antoni | A bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF] -to be shortly abandoned from section 4 -is to resort to a complex Gaussian, mainly because it makes analytic calculation tractable and also because it coincides with classical Tikhonov regularization. Namely,

[s i ] = e -α -2 s H i Σ -2 0 s i α 2Q π Q |Σ 2 0 | = N C (s i ; 0, α 2 Σ 2 0 ) (8)
where α 2 Σ 2 0 = E{s H i s i } stands for the prior covariance matrix of the coefficients, α 2 for their mean power and Σ 2 0 for a matrix of correlation coefficients conventionally normalized such that trace{Σ 2 0 } = Q. The fact that α 2 Σ 2 0 does not depend on index i reflects the assumption of a stationary acoustic field, i.e. with constant statistical properties during the acquisition of the N snapshots. Using Eq. ( 6), one then finds the posterior PDF in the form of a complex Gaussian, [s i |p i ] = N C (s i ; ŝi , Ω), with posterior mean and covariance matrix

ŝi = Σ 2 0 G H (GΣ 2 0 G H + η 2 I) -1 p i , η 2 = β 2 α 2 (9
)

Ω = β -2 G H G + α -2 Σ -2 0 -1 .
In this case, the posterior mean ŝi coincides with the MAP estimate as well as with the Tikhonov regularized leastsquare solution (LS) returned by ŝi = Argmax

s i [p i |s i ][s i ] = Argmin s i (-ln[p i |s i ][s i ]) = Argmin s i ||p i -Gs i || 2 2 + η 2 Σ -1 0 s i 2 2 (10) with Σ -1 0 s i 2 2
.

= s H i Σ -2 0 s i and wherein all terms not depending on s i have been dropped out. The ratio η 2 = β 2 /α 2 in the above equations acts as the Tikhonov regularization parameter which controls the stability of the solution. One benefit of the Bayesian framework is to provide a solution to automatically tune η 2 , as described in Ref. [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF] and resumed in Appendix B.

It is explicit from Eq. ( 9) that the MAP estimate directly depends on the setting of the prior covariance Σ 0 . This aspect is crucial in the present paper, since it is the mechanism by which sparse solutions will be promoted as explained in the next subsections.

Definition of a spatial prior through the aperture function

Although the exact source distribution will remain unknown until the inverse problem is solved, prior information is often available before the experiment is run, even in a vague form. Typically, spatial regions may be identified where sources are expected to radiate from -or not -and a certain degree of spatial correlation may be specified on the expected source field. Such elements of information constitute the "spatial prior".

Of specific concern in this paper is the case where a majority of the source coefficients are nil (or nearly so) whereas a few coefficients are allowed to take large values. This is referred to as a "sparse" source distribution. Although this situation will be properly addressed in section 4 by introducing specific non-Gaussian prior PDFs [s i ], a necessary condition to meet from the onset is to assume a priori that the source distribution is incoherent on the surface Γ s . This will return the finest spatial resolution among all possible spatial priors. It is important to emphasize at this juncture that the assumption of prior incoherence does not imply that the reconstructed sources are a posteriori incoherent as well: correlation among the estimated source coefficients is truly allowed as seen immediately from the structure of Eq. ( 9). Indeed, most methods of the literature (e.g. such as NAH or ESM) happen to assume prior incoherence in an implicit manner when interpreted from a Bayesian perspective.

Assuming time stationarity and a finite-power (i.e. with integrable squared-magnitude) spatial spectrum, the spatial correlation of an incoherent source field has the functional form E s i (r q )s i (r q ) * = α 2 σ 2 0 (r q )δ qq [START_REF] Brooks | A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays[END_REF] where δ qq stands for the Kronecker delta and σ 2 0 (r q ) reflects the intensity of the field at position r q . The continuous function σ 2 0 (r), r ∈ Γ s , arbitrarily normalized to unit volume, Γ s σ 2 0 (r)dΓ(r) = 1, is called the "aperture function" (AF) [START_REF] Antoni | A bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF] for reasons to become clear soon and α 2 the prior "source power". This allows for possible non-homogeneity of the source field with prior local variance E{|s i (r)| 2 } = α 2 σ 2 0 (r). Therefore, spatial regions with high likelihood of radiation are assigned a relatively high value of σ 2 0 (r) whereas regions less likely to radiate are given a value of σ 2 0 (r) close to zero. Based on this model, one has the Q × Q diagonal prior covariance matrix

Σ 2 0 = Diag σ 2 0,q , (12)
where σ 2 0,q . = σ 2 0 (r q ). It is proved in Appendix A that this parameterization of the problem actually minimizes the mean-square error of the estimated source. In a discrete setting, this implies the following probabilistic model for the source coefficients, s q,i = α • σ 0,q • q,i , q,i ∼ N C (0, 1), [START_REF] Sijtsma | CLEAN Based on Spatial Source Coherence[END_REF] where i,q denotes a standardized (i.e. zero-mean and unit variance) complex Gaussian random variable. The latter encodes the fluctuations of the source coefficient in time -hence its dependence on the snapshot index i -whereas the (squared) product α 2 σ 2 0,q returns the source power at spatial position r q independently of i. Since σ 2 0,q is normalized, it is convenient to endow the global variable α with the physical dimension of the sources s q,i .

Eventually, more insight is provided by the functional notation of Eq. ( 9) as allowed by the introduction of the AF σ 2 0 (r). It is seen that the source at position r is then estimated as

ŝi (r) = σ 2 0 (r) M m=1 G(r|r m ) * q m,i (14)
where the coefficients q m,i are the elements of vector

q i = (G + η 2 I) -1 p i with [G] kl = Γ s G(r k |r)G(r l |r) * σ 2 0 (r)dΓ(r). (15)
Formula ( 14) explicitly shows that the number of "degrees of freedom" of the estimated source distribution is equal to the number of microphones. Its physical interpretation is that ŝi (r) is estimated by constructively interfering on surface Γ s the waves radiated by the M microphones as if they were point sources with "strength" q m,i -from the principle of reciprocity, the conjugated Green function G(r|r m ) * appearing in Eq. ( 14) means that the direction of propagation is reversed, from the microphone position r m to the source coordinate r. At last, the estimated source is focused on the region of interest by multiplication of the sum in Eq. ( 14) with the AF σ 2 0 (r).

Focusing mechanism and its link to wideband acoustical holography

The AF σ 2 0 (r) plays a fundamental role in BF. By shrinking the reconstruction of the source distribution on a restrained spatial region where sound is a priori supposed to radiate from, it can considerably improve the estimation performance in terms of spatial resolution, quantification of sound levels, and frequency range. The improvement in spatial resolution is an obvious consequence of using a narrow AF. If S σ denotes the effective area of the surface covered by σ 2 0 (r), then the spatial resolution of ŝ(r) is on the order of S σ /M. The improvement in the quantification of sound levels is perhaps more unexpected. It may be explained as follows: since a narrow AF alleviates the effort of localizing the sources, more "degrees of freedom" are left available to improve the accuracy of the source amplitudes. An extreme example is when a small number -say N s -of point sources are to be estimated with a priori known positions; the AF is then shrunk into a few Dirac deltas at the source positions so that the inverse problem boils down to estimating N s amplitudes from M observations. As long as N s ≤ M the problem is potentially over-determined and can be solved more accurately than if the full spatial distribution coded on Q spatial samples had to be inferred.

Eventually, the use of a narrow AF can also enlarge the frequency range of the reconstructed source distribution. This is particularly compelling since traditional approaches in inverse acoustics are known to be fundamentally limited downward and upward in frequency. The reason may be better explained trough an example.

Example 1. Let us consider the case of a monopole source on a linear domain Γ s = {r : -1 ≤ r ≤ 1} parallel to a uniform linear array (ULA) of M = 11 microphones located at a distance R = 0.35 m -see Fig. 1. A "high" frequency f = 6 kHz is purposely selected that makes the resolution of the inverse problem difficult in this configuration. The signal-to-noise ratio (SNR) is set to 20 dB. The situation without AF is first considered, which is implicitly equivalent to using an uniform AF with unit amplitude over Γ s . The reconstructed source ŝ(r) (whose real and imaginary parts are shown by the black and grey curves, respectively, in Fig. 2(a)) exhibits marked fluctuations over the full domain Γ s and significantly differs from the theoretical monopole. As a result, the pressure field radiated by ŝ(r) (real part shown in Fig. 2(a)) comprises a lot of fictitious interferences. The squared modulus of the radiated pressure field (shown in Fig. 2(d)) -which represents one part of the acoustic energy -evidences a complex radiation of the energy from Γ s towards the array and some hot spots are seen at the positions of the microphones. This can be understood as an artifact of the classical LS solution with Tikhonov regularization when no AF is used. In this case, the regularization term tries to globally minimize the energy of the source while still matching the measurements at the microphones. Therefore, a source distribution is found with a directivity pattern that primarily points towards the microphone positions. As a result, the average acoustic power produced by the reconstructed source is also largely underestimated: in the present case, it is found equal to only 14.4% of the true value.

The situation is next considered where a narrow AF is inserted in the form of a Hann function (see Fig. 2(c)) in order to encourage the reconstruction of the source in a small neighborhood of r = 0. The reconstructed source (real and imaginary parts shown by the black and grey curves, respectively, in Fig. 2(b)) now resembles much more the theoretical monopole and it radiates a pressure field (real part shown in Fig. 2(b)) that correctly reproduces outgoing spherical waves. The squared modulus of the radiated pressure field (shown in Fig. 2(e)) also correctly follows the spherical spreading law of a monopole, without artificial directivity towards the microphones. The average acoustic power is now found equal to 97% of the true value.

The benefit offered by the insertion of a narrow AF can be understood from a different perspective by introducing the concept of a representation basis. From the interpretation following Eq. ( 14), the source estimate ŝ(r) is constructed by interfering waves radiated from the microphones and focused through the AF on the source surface Γ s . Therefore, these waves provide a spatial basis,

{φ k (r) = σ 2 0 (r)G(r|r k ) * , r ∈ Γ s } M k=1
, for representing ŝ(r) (note that the basis is not necessarily complete). Figure 3(a) illustrates how one basis function is produced by the radiation of one microphone. The full family of basis functions relating to Example 1 is illustrated in Fig. 3(b-c). In the absence of an AF -or equivalently when the latter is set uniform -the basis functions reduce to the Green functions, {G(r|r k ) * , r ∈ Γ s } M k=1 , which oscillate in Γ s with wavelength λ ∼ c 0 / f . Because the array is parallel to Γ s , the basis is composed of the same waveform shifted at different locations. The corresponding spatial spectra are displayed in Fig. 4(a), which all have the same magnitude and only differ by a linear phase. It is seen that the spectra have a marked peak around the spatial frequency f /c 0 , which explains why the source distribution in Fig. 2(a) has been reconstructed with so many ripples. On the contrary, the insertion of a narrow AF produces a much richer representation basis, quite close to "Gabor wavelets" as seen in Fig. 3(c). The corresponding spatial spectra (see 4(b)) are close to a filterbank decomposition, where each basis function spans a different frequency band. In addition, the spectral range has been enlarged, which explains why a sharper source can be reconstructed. This is a manifestation of the uncertainty principle stating that the product of the widths of the AF,say ∆r, and of its spatial spectrum, say ∆k, is upper bounded, i.e. ∆r∆k ≤ C for some constant C. There is indeed no limit to the highest spatial frequency that can be reconstructed this way since it is inversely proportional to the AF width which can be tuned arbitrarily narrow -i.e. ∆k ∼ C/∆r. In particular, a monopole can be reconstructed perfectly when the AF ultimately tends to a Dirac delta. This phenomena -although not explained from this point of view -has been coined "wideband holography" in Ref. [START_REF] Hald | Fast wideband acoustical holography[END_REF].

k/2π (m -1 ) k/2π (m -1 ) f c

Sparse Bayesian Focusing

BF, as reviewed in the previous section, involves a manually-tuned AF. If properly designed, the AF can greatly improve the quality of the reconstructed source distribution in terms of localization, quantification and frequency range. The question then arises as whether it is possible to automatically tune the AF, by considering it as another unknown quantity in the inverse problem. An intuitive solution consists in first solving the inverse problem with a large AF and next in using the so-estimated source distribution as a new AF to solve the inverse problem a second time. By iterating the process several times, the initial AF progressively tunes itself so as to focus on the different point sources that constitute the source distribution, as depicted in Fig. 5. This section justifies the strategy from a theoretical point of view in the Bayesian setting. It also demonstrates that it applies with various adaptation rules depending on the choice of the prior PDF of the source coefficients. The principle is to consider the AF as a random quantity that is inferred in an iterative resolution of the inverse problem. How to model the AF as random quantity is first addressed in subsection 4.1. By assigning the AF with a PDF, various source priors can thus be generated, different from the Gaussian one of Eq. ( 8), which promote sparse sources; they are explored in subsections 4.2. Subsection 4.3 then introduces an iterative algorithm -coined IBF -which jointly infers the source coefficients and updates the AF.

Γ s s [0] (r) σ 2 (r) 0 σ 2 (r) [1]
σ 2 (r) [START_REF] Fuchs | Linear programming in spectral estimation. application to array processing[END_REF] s [START_REF] Starck | An overview of inverse problem regularization using sparsity[END_REF] (r) ...

Sparse spatial prior

In the context of this paper, the sources coefficients are assumed to be sparsely distributed in space; this means that only a very elements in vector s i (always the same independently of snapshot i) are expected to have significant values, the other ones being close (or ideally equal) to zero. One way to enforce this property is to assign to s i a prior PDF with heavy tails and very peaked at zero as compared to the Gaussian, which is denoted hereafter as leptokurtic (i.e. with a kurtosis greater than the Gaussian). The use of leptokurtic prior PDFs (sometimes referred to as sparse priors, sparsity enforcing priors or shrinkage prior) within the Bayesian framework for enforcing sparse solutions has been extensively studied in the recent statistical literature and is also known as sparse Bayesian learning1 (see e.g. [START_REF] Ji | Bayesian compressive sensing[END_REF][START_REF] Babacan | Bayesian compressive sensing using laplace priors[END_REF]).

A physical interpretation is given hereafter in terms of SMoGs. Following the same lines as in subsection 2.3, let the source coefficients have the following probabilistic model,

s q,i = α • σ 0,q • τ q • q,i , q,i ∼ N C (0, 1), ( 16
)
where the source power α 2 , the AF σ 2 0,q and the random variables q,i have been defined in Eq. ( 13). Equation ( 16) separates the source coefficients into time-dependent (index i) and space-dependent (index q) variables. All variables are dimensionless, except α which carries the dimension of s q,i . All variables are also seen as random, except σ 2 0,q which is "deterministically" designed by the user to reflect a global (or "initial") AF (e.g. see dotted bell-shaped curves in Fig. 7); it is introduced for the sake of generality, yet it could just as well be ignored to simplify the exposition. As compared to Eq. ( 13), τ q is now an additional non-negative random variable, which reflects the "relative intensity" at the spatial position r q (see Fig. 6). It may also be interpreted as the variance of τ q • q,i ∼ N C (0, τ 2 q ), which is the mechanism by which "local" aperture functions can be produced at position r q . It is noted at this juncture that α cannot be absorbed into τ q because the latter will be shortly seen to act as a "sparsifying" variable, with high degree of sparsity achieved for (very) small values of its scale parameter; this means that the scale parameter of τ q must be left free and the global scaling assigned to α. The stochastic model ( 16) is actually equivalent to considering right out a random AF

σ 2 (r) = σ 2 0 (r)τ 2 (r), (17)
(given here in functional notation with τ q . = τ(r q )) instead of a deterministic one as done in the previous section. The random AF in Eq. ( 17) yields the conditional prior variance E{|s

i (r)| 2 |τ 2 (r)} = α 2 • σ 2 0 (r) • τ 2 (
r) at position r on Γ s , independently of the snapshot index i, which locally and constantly encourages high values of the source coefficients where the relative intensity τ 2 (r) is large and tends to switch them off where the relative intensity is small. As compared to the uniform relative intensity used in BF (see Eq. ( 11)), the resulting spatial distribution of the source coefficients is thus sparser.

It is noteworthy that the proposed probabilistic model allows enforcing sparsity directly on the modulus of the source coefficients. This is advantageous as compared to the situation imposed by some solvers for 1 minimization, implemented for real-value data and whose extension to the complex case is not trivial. Sparsity on the moduluswhich is equivalent to group sparsity on the real and imaginary parts -easily allows switching off a large number of source coefficients. On the contrary, enforcing separate sparsity on the real and imaginary parts -as would be done by mapping a complex number to a two-dimensional real vector -can only switch off a source coefficient when both the real and imaginary happen to vanish simultaneously, which occurs with a much lower probability. In addition, separate sparsity would produce several instances of exactly nil real (resp. imaginary) parts together with non-zero imaginary (resp. real) parts, which reflects unrealistic physical combinations.

Although it is difficult to specify the values of τ 2 q before the inverse problem is solved, it is possible to assign them a hyperprior PDF, τ 2 q , that allows local variations in space. Given a spatial position r q and considering all snapshots i = 1, ..., N observed at this position, the prior PDF of the source coefficients S . = {s q,i } N i=1 then reads

[S] = ∞ 0 S|τ 2 q τ 2 q dτ 2 q = E τ 2 q        N i=1 N C s q,i ; 0, α 2 σ 2 0,q τ 2 q        (18)
where E τ 2 q means that the expected value is taken with respect to the random variable τ 2 q . Inspection of Eq. [START_REF] Chu | A robust super-resolution approach with sparsity constraint in acoustic imaging[END_REF] shows that the prior PDF of the source coefficients is a continuous mixture of Gaussians N C s q,i ; 0, α 2 σ 2 0,q τ 2 q weighted by the "scale" τ 2 q -the so-called SMoG. This is in general no longer a Gaussian, but a leptokurtic PDF. Figure 6 illustrates how a leptokurtic PDF is obtained by mixing Gaussians with random variances distributed according to a given law. It is seen that PDFs with very heavy tails can thus be obtained by allowing the variance to be broadly distributed.

2 2 2 0,q q α σ τ × × 2 q
σ updated aperture function source power relative intensity initial aperture function

0 iii) [s] = [s|τ 2 ][τ 2 ]dτ 2 τ 2 i) [τ 2 ] ii) [s|τ 2 ] s ∫ a) b)
Figure 6: a) Structure of the variance of the sparse spatial prior, where index q relates to spatial position r q . α 2 is the global source power and the initial aperture function σ 2 0,q reflects its expected spatial distribution before the experiment is run. The spatial variance α 2 σ 2 0,q is used in the Gaussian prior PDF of BF; its multiplication by the relative intensity τ 2 q returns the spatial variance α 2 σ 2 0,q τ 2 q used in the SMoG prior PDF of Iterated Bayesian Focusing. The method updates the aperture function as σ 2 q ← σ 2 0,q τ 2 q . b) Principle of construction of an SMoG. i) PDF [τ 2 ] of the variance τ 2 , ii) isoprobabiliy contours of the conditional Gaussian PDF, [s|τ 2 ] = N C s; 0, τ 2 , as a function of the random variable s and of the variance parameter τ 2 , iii) PDF [s] = [s|τ 2 ][τ 2 ]dτ 2 of the random variable s averaged over all values of the variance τ 2 (continuous curve) together with the Gaussian PDF with same variance (dotted curve). It is seen that a leptokurtic PDF is produced, with heavy tails and a peaked mode.

Examples of prior source fields simulated with the SMoGs investigated in the next subsection are displayed in Fig. 6. In order to allow an easy representation, the source domain is limited to a line segment. The initial AF σ 2 0,q has been designed in the form of a Hann function that progressively downsizes the values at the extremities of the domain where sound is less likely to radiate from. In each case, the value of α 2 is tuned so as to impose a unit variance (spatially averaged over the source domain). The source coefficients are displayed as a function of the space position for a given snapshot i. It is seen that the real and imaginary parts are jointly sparse, which means that they both have the same variance proportional to the updated AF σ 2 0,q τ 2 q . In other words, sparsity is enforced on the modulus -and not on the real and imaginary parts, independently -in order to force a majority of source coefficients to take small values and, ideally, to switch off completely.

-1 0 1 -1 0 1 -1 0 1 r q r q a) b) c) d) e) f)
Figure 7: Examples of one-dimensional cuts through complex random fields synthesized with different SMoG PDFs (the real and imaginary parts are displayed by black and gray curves, respectively). a) Complex Gaussian, b) complex Student-t (rooted on the inverse Gamma PDF [τ 2 ] with a = b = 1), c) complex K-distribution (rooted on the Gamma PDF [τ 2 ] with a = b = 0.05), d) complex Laplacian (generalized complex Gaussian with p = 1), e) generalized complex Gaussian with p = 0.1, f) mixture of complex Gaussians with π = 0.1, τ2 1 = 0.1 and τ2 2 = 10. Also shown is the initial AF σ 2 0,q (dotted curve).

Prior PDFs in the SMoG family

Different prior PDFs are obtained depending on the choice of the hyperprior τ 2 q in Eq. ( 18). Although closedform results exist in the literature in the case of real-valued random variables, they must be extended here to the scenario with complex-valued source coefficients and with several snapshots -i.e. to the complex multivariate case. Explicit formula are derived in Appendix C for some particular cases of interest and resumed in Tab. 1. As far as the authors know, some of these formula are given here for the first time.

It is noteworthy that the SMoG model allows expressing the prior PDF [S] by means of the power spectrum of the source coefficients only. Specifically, the quantity

χ 2 q = N S qq α 2 σ 2 0,q , (19)
which involves the ratio of the source power spectrum S qq (see Eq. ( 2)) to the prior spatial variance α 2 σ 2 0,q at position q, is a sufficient statistic. It will be used from now on as a shortcut notation.

Prior source distributions of Tab. 1 are illustrated in Fig. 7 and also displayed in Fig. 11 (a-b) for N = 1 and N = 50 snapshots. It is verified that all PDFs are apt to produce very sparse fields, even though the Laplacian (corresponding to the 1 penalty) seems slightly inferior than the other candidates.

Hyperparameters

Sparse prior PDF [S] = S|τ 2 q τ 2 q dτ 2 q Inverse Gamma shape parameter a > 0 scale parameter b > 0 Multivariate complex Student-t (MCS)

Gamma

Inference in the SMoG family

The price to pay for the introduction of a non-Gaussian prior PDF is that, generally, a closed-form solution no longer exists for estimating the source coefficients. However, it is shown hereafter that a point estimate can still be expressed with the same structure as in Eq. ( 9), yet where Σ 2 0 is replaced by a different AF. This is a direct consequence of the fact that an SMoG is a weighted average of Gaussian priors as expressed in Eq. [START_REF] Chu | A robust super-resolution approach with sparsity constraint in acoustic imaging[END_REF]. To see this, let us introduce the notation ŝi

(Σ 2 ) = Σ 2 G H (GΣ 2 G H + η 2 I) -1 p i ( 20
)
for the classical penalized LS solution given in Eq. ( 9) rooted on a given AF, Σ 2 = Diag(σ 2 q ) (note that, as a particular case, Eq. ( 9) is exactly recovered with Σ 2 = Σ 2 0 ). Then, it is established in Appendix D that the posterior mean under an SMoG prior is

E{s i |P} = ŝi (Σ 2 )[Σ 2 |P]dΣ 2 (21)
where [Σ 2 |P] is the posterior PDF of the AF given the measurements P. Equation ( 21) is recognized as a weighted average of ŝi (Σ 2 ) over realizations of the AF conditioned on the data. More importance is therefore given to solutions whose posterior PDF [Σ 2 |P] is high. In the extreme case where the latter PDF is very peaked, the posterior mean may be evaluated directly from Eq. ( 20) wherein the AF is replaced by its MAP estimate Σ 2 MAP = Diag(σ 2 0,q •Arg max[τ 2 q |P]). Alternatively, it is proved in Appendix D that the MAP estimate under an SMoG prior -which will be considered from now on -is ŝi (Σ 2 ) with AF

Σ 2 = Diag         σ 2 0,q E τ -2 q | Ŝ         (22)
where Ŝ .

= {ŝ i } N i=1 . This means that the MAP is obtained by updating the AF σ 2 0,q by dividing it with the expected value of τ -2 q conditioned on Ŝ. This is expected to provide a "sharper" AF than the initial one. It appears in light of these results that the MAP estimate will generally differ from the posterior mean in the SMoG family -contrary to the Gaussian case of section 2 -which reflects the fact that [s i |P] is not a symetric PDF.

Although inserting Eq. ( 22) into ( 20) returns an implicit equation with no closed-form solution, it can be solved iteratively as shown in the next subsection.

At this juncture, the inverse problem has been formally reformulated as inference in a Bayesian hierarchical model with a structured prior, as depicted in Fig. 8. In short, the essential difference with BF is that the AF is now an unknown quantity to be inferred simultaneously with the source coefficients.

Iterated Bayesian Focusing (IBF)

An algorithm is proposed hereafter to reach the MAP estimated given in Eq. [START_REF] Chardon | Near-field acoustic holography using sparse regularization and compressive sampling principles[END_REF]. The IBF algorithm is described hereafter in terms of the normalized sum of squares

χ 2 q,[k-1] = N S qq,[k-1] α 2 [k-1] σ 2 0,q (23)
which involves the power spectrum

S qq,[k-1] = 1 N N i=1 | ŝq,i,[k-1] | 2 of
the source coefficients and the source power α 2

[k-1] estimated at iteration k -1. The reason for introducing this variable is manifold: beyond greatly simplifying forthcoming notations, it is also a sufficient statistics which plays a central role in the IBF algorithm. As seen from its definition, it also has the physical significance of the "normalized" power of the sources.

Algorithm 1

•
Step 0: Design an initial AF σ 2 0,q and construct the diagonal matrix Σ 2 0 whose q-th diagonal entry is σ 2 0,q .

• Step 1: Set k = 0. Estimate the source power α 2 [0] and the regularization parameter η 2 [0] by using Eqs. (B.3) and (Appendix B) evaluated with Σ 2 0 . Initialize the algorithm with the MAP solution ŝi,[0] of Eq. ( 9) as returned from a Gaussian prior.

Repeat steps 2 to 5 until convergence.

• Step 2: Do k ←-k + 1.
Estimate the relative intensity as

τ2 q,[k] = -       ∂ ∂(χ 2 q ) ln χ 2 q       -1 χ 2 q =χ 2 q,[k-1] , (24)
where [χ 2 q,[k-1] ] is the prior PDF expressed as a function of the normalized sum of squares χ 2 q,[k-1] at iteration k -1 (closed-form expressions of τ2

q,[k] are given by Eqs. (G.1), (G.2), (G.4), (G.5), (G.6), (G.8), and (G.9)).

• Step 2(bis): Set τ2 q,[k] = max τ2 q,[k] ; 1 • max q τ2 q,[k]
where 0 < 1 1.

• Step 3: Update the AF as

σ 2 q,[k] = σ 2 q,0 τ2 q,[k] (25)
and construct the diagonal matrix Σ 2

[k] whose q-th diagonal entry is σ 2 q,[k] .

• Step 4: Estimate the source power and the regularization parameter as

α 2 [k] , η 2 [k] = Argmax α 2 , η 2 |S pp , Σ 2 [k] (26)
by using Eqs. (B.3) and (Appendix B) evaluated with the updated AF.

• Step 5: Estimate the source coefficients as

ŝi,[k] = Σ 2 [k] G H (GΣ 2 [k] G H + η 2 [k] I) -1 p i . ( 27
)
Convergence criterion:

Stop iterations when ||ŝ i,[k] -ŝi,[k-1] ||/||ŝ i,[k-1]
|| is smaller than a given threshold, 0 < 2 < 1.

(Step 2(bis) is optional with the MCS and MMCG priors, yet it is recommended with the MCK, MCL and GMCG priors as explained later in subsection 5.4.) IBF belongs to the family Iterative Re-weighted Least Squares (IRLS) algorithms [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm[END_REF], since at each iteration it solves a weighted least square problem whose solution is given by Eq. ( 9) with Σ 2 0 replaced by it updated version Σ 2

[k] as obtained from Eq. ( 25). Yet, it is more rigorously derived in Appendix E as an Expected-Maximization (EM) algorithm, which makes possible to prove its convergence to the MAP solution -insightful connections between IRLS and EM have been investigated in Ref. [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF]. The proposed IBF algorithm has the following features:

1. It applies generally in the SMoG family, as long a closed-form expression of the prior PDF is available. Consequently, different structures and degrees of sparsity can be promoted by choosing among the various possible prior PDFs of subsection 4.2. 2. The regularization parameter η 2 is iteratively updated together with the source coefficients with proved convergence (see Appendix B). This is to be contrasted with other approaches found in the literature where the regularization parameter is either arbitrarily imposed or tuned afterwards, typically by using cross-validation or the L-curve. The latter approaches require running the optimization as many times as values of the regularization parameter are to be tested, and are therefore computationally more demanding -if not prohibitive in some applications such as that illustrated in subsection 6. In addition, the estimation of the regularization parameter by means of the Bayesian criterion used in Step 4 has been found more efficient than the cross-validation or the L-curve in many inverse acoustic scenarios [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF]. 3. The algorithm applies to complex-valued source coefficients and sparsity is enforced to their modulus (and not to the real an imaginary coefficients independently). Consequently, a majority of sources are forced to switch off. 4. A key step in the algorithm is the update of the AF in Eq. ( 25) rooted on the current estimate of the relative intensity given by Eq. ( 24). Closed-form expressions of the updating rules corresponding to the prior PDFs of Tab. 1 are given in Appendix G. 5. The updating rule involves the consideration of possibly several snapshots and therefore involves the power spectra of the source coefficients -see Eq. ( 23). Accordingly, the IBF algorithm can be formulated in terms of the CSM of the measurements, S pp , only (see Eq. ( 1)). As explained in the introduction, this may be determinant in practice, since some commercial data acquisition system used in acoustic imaging only record the CSM by default. The algorithm is as follows.

Algorithm 2

Step 1 to 4: Unchanged.

Step 5: Estimate the power spectrum of the source coefficients at position r q as

S qq,[k-1] = σ 4 q,[k] g H q (GΣ 2 [k] G H + η 2 [k] I) -1 S pp (GΣ 2 [k] G H + η 2 [k] I) -1 g q (28)
where g q is the q-th column of matrix G.

Example 2. The IBF algorithm is illustrated on Example 1 of section 3 which involved the reconstruction of a monopole source at f = 6 kHz with an ULA placed at a distance R = 0.35 m (see Fig. 2) under an SNR of 20 dB. The MCS prior is used with hyperparameters a = b = 0.01 (comparison between the several priors introduced in subsection 4.2 will be addressed in the next section). To start with, only N = 1 snapshot is recorded. The IBF algorithm is initialized with the BF solution of section 2 with an AF σ 2 q,0 shaped like a Hann window (raised cosine) over the domain Γ s = {r : -1 ≤ r ≤ 1} (with discretizaton step of 0.01 m). The algorithm is stopped after 10 iterations. Figure 9(a-b) displays the evolution of the AF σ 2 q,[k] on linear and dB scales with respect to the iteration number k. It is seen that the initial AF quickly tuned itself into a narrow peak around the actual source position shown by a red dot. At the same time, the dynamic range of the AF improves from 5 dB (for σ 2 q,0 ) to more than 70 dB (for σ 2 q, [START_REF] Gerstoft | Multiple and single snapshot compressive beamforming[END_REF] ). As explained in section 3, the shrinkage of the AF introduces a focusing mechanism that improves the accuracy of the estimated acoustic power. Whereas it was only 17.6% of the true value at the first iteration (note a slight improvement as compared to 14.4% for the uniform window used in section 3), it is 99.9% at the ninth iteration.

Figure 10 displays the spectra of the estimated acoustic power (relatively to the true value) for several working frequencies in the range [0.02; 50] kHz. The monopole is driven by a white Gaussian signal and results are averaged over 100 independent draws of complex Gaussian noise added to the measurements with an SNR of 20 dB. As expected, the BF solution has an upper frequency limit on the order of f max ∼ c 0 /d 2.3 kHz with d = 0.15 m the spacing between the microphones. On the other hand, the IBF solution does not show such a limit. In the present case with a monopole source, it is virtually unlimited in frequency since the AF is allowed to become arbitrarily narrow. The IBF solution is remarkably accurate over the full frequency range as compared to the BF one. It is also found more precise with a 3-sigma confidence interval about half as large as that of the BF solution. In addition, the coefficient of variation (CV) is constant (CV = 0.44%) as compared to that of the BF solution which increases when approaching the upper frequency limit f max .

The same experiment is now run with two incoherent monopole sources driven by white Gaussian signals in order the check the ability of the IBF algorithm to self-tune the AF around several distinct source points. The first monopole is located at r = 0 m and the second one at r = -0.05 m. They are both fed with mutually independent complex Gaussian random noises. The level of the second monopole is set to 4% the level of the first one. In this case, N = 100 snapshots are recorded in order to better estimate the partial power due to each monopole (i.e. to average out incoherent cross-terms). The IBF is run with the same settings as before. Figure 9(c-d) displays the evolution of the AF on linear and dB scales; it is seen to quickly split in two narrow peaks anchored to the two monopole positions. As before, the dynamic range after 10 iterations is greater than 70 dB. The two peaks have different magnitudes -which actually makes it difficult to see the smallest one on the linear scale although its presence is obvious on the dB scaleyet their respective areas converges to exactly the same ratio as the levels of the two monopoles, i.e. 0.04. This might be an issue only for detecting very small sources masked by a dominant source (here below 70 dB of the largest one), a situation which places the performance of IBF much higher than CV or BF in this respect. Despite the two peaks of the updated AF having different magnitudes, the estimated source coefficients are perfectly calibrated to recover the correct strengths of the two monopoles. After convergence, the acoustic power radiated by the first monopole is found equal to 100.1% of theoretical power and that of the second one to 100.7%. This is to be compared to 8.6% and 304.6%, respectively, for the initial estimation at k = 0, which on the one hand still underestimates the largest source but on the other hand largely overestimates the smallest one due to the presence of sidelobes from the dominant source.

Comparison of priors

The use of an SMoG model has been shown to yield several possible candidates for the prior PDF of the source coefficients. The question naturally arises as which one(s) should be selected in practice. Although the exact answer to this question is surely problem dependent, some considerations are possible based on the following arguments.

Number of hyperparameters

Unless the user has a very accurate knowledge of the prior PDF, it is often convenient to adopt a choice that involves the fewest hyperparameters. In this respect, The MCL arrives in first position without any hyperparameter to set. The GMCG arrives in second position with only one hyperparameter, p, which directly reflects the amount of expected sparsity. The MCS and the MCK come next with two hyperparameters, a and b. However, as seen in Eq. (G.1), the effect of hyperparameter in the MCS tends to vanish when the number of snapshots becomes large. Therefore, a convenient choice is to set a = 0 and keep b as the only hyperparameter to tune; another possibility is to set a = b to a small value that promotes sparsity. No similar simplification seems to hold for the MCK, for which the presence of two hyperparameters is the price to make it very flexible. The MMCG comes last with three hyperparameters, τ1 , τ2 and π. It has been checked in simulations that their ratio is really determinant, so that one can set τ2 = 1 for instance from the onset. These findings are resumed in Tab. 2.

Sensitivity to hyperparameters

The hyperparameters of the SMoG model arrive in the last level of the Bayesian hierarchical model (see Fig. 8) and are therefore, in theory, much less influential than hyperparameters in lower levels such as typically the regularization parameter η 2 to which the inverse problem is known to be extremely sensitive (see e.g. Ref. [START_REF] Gelman | Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)[END_REF] for an excellent discussion on the influence of hyperparameters in high-level priors). As explained in the previous subsection, the general recommendation is to tailor the prior PDF to make it sparse. Extensive experiments have shown that, on the one hand, the MCS prior is very stable when a and b are set small as indicated in Tab. 2. On the other hand, the MCK, the MMCG and the GMCG are quite sensitive to the setting of their hyperparameters. This is indeed the price to pay for the PDF to gain flexibility in its shape. The GMCG has yet an advantage in that it is controlled by only one parameter, the power p, with perhaps a clearer meaning than for the MCS and MCK hyperparameters. In this respect, the GMCG prior is a relevant choice if the prior PDF has to be carefully tuned in a simple way, as will be evidenced in the experiments of section 7. Automatic tuning of the power p remains challenging; although theoretically possible in the Bayesian setting [START_REF] Aucejo | On a full Bayesian inference for force reconstruction problems[END_REF], this issue is outside the scope of the present paper.

Bias

The choice of a particular SMoG has a direct impact on the estimation error of the source coefficients. This may be investigated by considering independently the bias and the variance of the MAP estimate. Based on the results of subsection 4.3, the bias can be expressed as (see Appendix I)

s i -E {ŝ i } = I + η -2 G H G • Diag σ 2 0,q E τ -2 q | Ŝ -1 -1 s i . ( 29
)
Therefore, it is seen that a reduced bias is generally achieved by SMoG models which can return high values of the relative intensity τ2 q . As seen in Fig. 11(c-d), this puts the MCS, the MCK with N > a, and the GMCG with p 0 (symbol means close to but greater than), in first position. On the opposite, the MMCG is expected to return the strongest bias as its relative intensity in Eq. (G.9) is bounded upward by τ2 2 . Incidentally, expression [START_REF] Pereira | Acoustic imaging in closed spaces[END_REF] for the bias also indicates that matrix G should be designed such that G H G is as close as possible to a diagonal matrix in order for the bias at position r q not to be influenced by the non-zero values of the source coefficients at other positions r q r q . This condition is actually reminiscent of the "restricted isometry property" [START_REF] Simard | Acoustic source identification: Experimenting the 1 minimization approach[END_REF] and "mutual coherence criterion" [START_REF] Xenaki | Compressive beamforming[END_REF][START_REF] Fernandez-Grande | A sparse equivalent source method for near-field acoustic holography[END_REF] of compressed sensing.

Eventually, it is seen from Eq. ( 29) that the bias is inflated by the regularization parameter η 2 , which is the price to pay for stabilizing an ill-posed inverse problem. Given a non-zero AF, the bias completely vanishes only in the particular case when η 2 is nil.

Variance and rate of convergence

The exact expression of the covariance matrix of the estimated source coefficients is not easy to derive. However, it can be shown within the framework of the EM algorithm that it directly depends on the so-called "missing information matrix" I m [START_REF] Mclachlan | The EM algorithm and extensions, 2nd Edition[END_REF], which measures the discrepancy of information due to the introduction of the unknown relative intensity τ 2 q in the EM algorithm as compared to the case where it would be known for all indices q: the larger the norm of I m , the larger that of the covariance matrix. Besides, the speed of convergence of the EM algorithm is similarly controlled by the norm of I m : the larger the "missing information" and the slower the convergence.

It is proved in Appendix J that the missing information matrix averaged over N snapshots takes the following expression

I m = 1 α 2 N Diag        χ 2 q σ 2 0,q • Var τ -2 q | Ŝ        where Var τ -2 q | Ŝ = ∂ 2 (∂χ 2 q ) 2 ln χ 2 q . ( 30
)
This result shows that the measure of missing information essentially depends on the conditional variance of the inverse of the relative intensity, Var{τ -2 q | Ŝ}, which is the price paid for sparsity. This unveils a connection between the missing information matrix and sparsity: missing information increases when more sparsity is introduced; intuitively, this corresponds to requiring more information to tune the unknown random variables τ 2 q when moving from a smooth LS solution (no sparsity) to a highly structured (highly sparse) solution. This equivalently reflects a balance between high sparsity and reduced estimation variance. Equation ( 30) also evidences that the variance decreases proportionally with the number of snapshots; therefore taking N large should be encouraged in general. The conditional variance of the PDFs introduced in subsection 4.2 are displayed in Fig. 11(e-f) for N = 1 and N = 50 snapshots.

An important requirement for the convergence of the IBF is to verify that the diagonal elements of I m remain bounded for all 0 ≤ χ 2 q < ∞. It is established in Appendix J that this is the case for the MCS and the MMCG, yet not for the GMCG with p < 2 neither for the MCK in general. This has a practical implication for the IBF algorithm based on the GMCG and MCK priors since the AF should be controlled in some way in order to keep it stable when χ 2 q takes very small values. One such strategy is to systematically set to zero the relative intensities τ2 q,[k] that converge below a small threshold and to remove them from the subsequent iterations. Another simple strategy is to replace them by the value of the threshold. This is the rationale behind Step 2(bis) in Algorithm 1.

In conclusion, the considerations of subsections 5.1, 5.3 and 5.4 place the MCS in first position since it 1) allows very sparse solutions while 2) requiring the setting of only one significant hyperparameter and 3) having a reduced bias and 4) an unconditionally bounded variance. 2) (black line), the MCK case with (a = 3N/4, b = 1) (red dotted line), the MCL case (yellow line) and the MMCG case with (π = 0.9, τ2 1 = 0.0001, τ2 2 = 1) (purple line) compared with the multivariate complex Gaussian (MCG) case (dotted black line). Sparse behaviors are essentially recognized by the distance to the Gaussian at the origin (χ 2 q = 0) and towards infinity (χ 2 q → ∞). c) Corresponding relative intensity τ2 q = E{τ -2 q | Ŝ} -1 = -(∂ ln([χ 2 q ]/∂(χ 2 q )) -1 for N = 1 and d) for N = 50, all normalized to 1 at χ 2 q /N = 1. e) Corresponding variance Var{τ -2 q | Ŝ} = ∂ 2 ln[χ 2 q ]/(∂χ 2 q ) 2 for N = 1 and f) for N = 50 with excerpt in the range 0 ≤ χ 2 q /N ≤ 0.06.

Example 3. Experimental comparison of the above priors is carried on in a configuration similar to Examples 1 and 2. The range distance is R = 0.5 m and an ULA with 20 microphones and spacing 0.1 m is used. The SNR is set to 20 dB. The source domain Γ s = {r : -1 ≤ r ≤ 1} (in meters) is discretized with a spatial step of 0.01 m. Initialization of the IBF algorithm is done as in Example 2 with a Hann window for the AF. A threshold set to 2 = 10 -3 times the highest amplitude is imposed with the GMCG and MCK priors in order to avoid possible divergence of the EM algorithm (see discussion of subsection 5.4). The threshold used in the stopping criterion of the iteration is 1 = 10 -3 . The MCS prior is used with a = b = 0.01, the MCK prior with a = b = 0.1 and the MMCG prior with π = 0.99 and τ2 2 /τ 2 1 either set to 10 or to 100. These parameters have been kept constant for all experiments in order to fairly test the robustness of the method in this respect. With these settings, the EM algorithm has been found to converge within a few iterations (about 10 on the average), with the fastest speed attained with the MMCG prior and the slowest one with the MCL prior.

The first experiment investigates the reconstruction of two incoherent monopoles with equal strengths, located at r 1 = -0.05 m and r 2 = 0.1 m, at frequency f = 1 kHz and with N = 1 snapshot. Reconstruction results are superimposed for 4 independent draws of the snapshot. Figure 12(a) displays the results of conventional beamforming and BF, which both suffer from a poor spatial resolution in this case. Apart from an offset due to a different scaling (the amplitude of beamforming returns the source strength whereas BF returns a source density), the two curves are quite similar. Figure 12(b) displays the result of the MCL prior (corresponding to the 1 -norm penalty), which improves the spatial resolution but still has a limited dynamic range. Numerous sidelobes are however present, probably due to the use of a uniform array, and a marked dispersion is observed with respect to snapshot draws. The estimation using the MCK prior is found very similar to the MCL prior in this case (Figure 12(e)). The GMCG ( 0 -norm penalty) and MCS priors in Figure 12(c-d) show the sparsest results, by refining even more the spatial resolution and greatly increasing the dynamic range. At the same time, they are found rather stable with respect to snapshot draws. The MMCG reconstruction shown in Fig. 12(f) shows intermediate sparsity, which have been checked to strongly depend on the value of the ratio τ2 2 /τ 2 1 . In some cases, only one of the two sources is identified. The second experiment repeats the first one, yet using the CSM computed with N = 100 snapshots. The main finding is that results have been considerably stabilized in accordance with the conclusions of subsection 5.4 (they are now almost invariant with respect to draws of the CSM). While the use the CSM has clearly upgraded the MCK estimation by placing it to close to the GMCG (p = 0) and MCS results, it has reduced the spatial resolution of MMCG. This is explained by the fact that the updating rule (G.9) tends to a hard threshold when N increases, which basically clips the source coefficients estimated in the first iteration of the IBF algorithm, whose magnitudes are smaller than γ c as given by Eq. (G.10).

The third experiment tackles a more difficult configuration with a very low frequency f = 100 Hz and with the first monopole strength 20% smaller than the second one (i.e. 14 dB difference in power level). Only results obtained with the CMS computed on 100 snapshots are displayed (see Fig. 14). It is seen that both beamforming and BF fail to detect the smallest source due to poor spatial resolution. Although the MCL can correctly identify the two sources, even better results are achieved by sparser priors such as the GMCG (p = 0) and the MCS (note a small bias on the localization of the small source). Both the MCK and the MMCG priors failed to reconstruct the small source with the given parameter settings. The spatial resolution of MMCG is only slightly better than that of beamforming.

The fourth experiment investigates another extreme case with a very high frequency f = 10 kHz and the same two unequal source strengths (see Fig. Beamforming and BF both fail in this case. The MCL prior clearly improves the situation, yet GMCG (p = 0) and MCS still return the best estimates with the deepest dynamic range achieved by the latter. Although sparse, the MMCG fail to correctly identify the two sources.

Similar conclusions have been obtained in several other experimental configurations not reported in the paper. In general, it turned out that sparser results than allowed by the 1 -norm penalty could be obtained by using the other priors investigated in the paper. However, only the GMCG ( 0 -norm penalty) and the MCS priors showed a constant stability of the estimates. On the contrary, the MMCG prior was found extremely dependent on the setting of its parameters τ1 and τ2 . Overall, results have been found remarkably independent of whether the point sources are incoherent or partially coherent (with the above settings, perfectly coherent sources of equal strengths could be resolved, while sources with 14 dB difference in power could be correctly identified up to 80% coherence). 6. Broadband Iterated Bayesian Focusing (BIBF)

So far, a different AF has been assumed at each frequency f . Nevertheless, it is meaningful to consider instances where the AF is independent of frequency, as would typically occur when the sound source is attached to a specific device (e.g. an obstruction, an opening, etc.). This situation has been considered for example in acoustic imaging in Refs. [START_REF] Wang | Wideband relax and wideband clean for aeroacoustic imaging[END_REF][START_REF] Gemba | Multi-frequency sparse bayesian learning for robust matched field processing[END_REF] and is reminiscent of the concept of "multi-frequency" group sparsity, where the same sparse structure is enforced to the source distribution for a group of frequencies.

Forcing a constant AF in a frequency band is straightforward within the proposed IBF framework. It is also found to yield improved results. The algorithm proceeds as follows (see Appendix K). Algorithm 3

• Define a frequency band B = { f 1 , ..., f B } as a collection of frequency bins (as typically returned by the discrete Fourier transform) where the self-tunning AF -as introduced in Eq. ( 17) -is invariant.

• Apply Algorithm 1 of section 4.4 where steps 0 and 3 are unchanged, step 2 is modified with

χ 2 q,[k-1] = B l=1 N i=1 | ŝq,i,[k-1] ( f l )| 2 α 2 [k-1] σ 2 0,q = N B l=1 S qq,[k-1] ( f l ) α 2 [k-1] σ 2 0,q (31)
substituted for χ 2 q,[k-1] , where ŝq,i,[k-1] ( f l ) stands for the evaluation of ŝq,i,[k-1] at frequency f l and snapshot i, steps 1, 4 and 5 are unchanged but performed for each frequency f k in band B.

Equation ( 31) is the main modification to address the broadband case. The quantity χ 2 q,[k-1] may be interpreted as the normalized energy of the source coefficients -as measured by their power spectra -summed in band B. A resolution in terms of the CSM is also trivially by proceeding as in Algorithm 3. Eventually, it is noteworthy that a different regularization parameter is reasonably assumed at each frequency.

Since BIBF solves the inverse problem at once for all frequencies in band B (say a total B frequencies), its computational cost is obviously increased as compared to one call of BF (for one frequency). However, it is lower than repeating Algorithm 1 successively for each frequency in the band, since only one AF is inferred instead of B. In addition, according to subsection 5.4, BIBF is expected to converge more rapidly than IBF solved at a single frequency since the corresponding missing information matrix (see Eq. ( 30)) is divided by an extra factor B. Overall, this endows BIBF with less computational effort than the repeated call of BF for scanning all frequencies in the band.

Example 4. BIBF is illustrated here for the reconstruction of four uncorrelated monopole sources driven by white noises in a wide frequency interval B ranging from 100 Hz to 20 kHz. The sources are located at r 1 = -0.1, r 2 = 0, r 3 = 0.1 and r 4 = 0.2 m on the source domain Γ s = {r : -1 ≤ r ≤ 1} (in meters) and have decreasing levels in the proportion 1, 0.25, 0.04 and 0.01 (i.e. 0dB, -6dB, -14dB, -20dB). The array is the same as in Example 3 and the range distance is R = 0.1 m. The SNR is set to 20 dB. BIBF is run with the CSM computed on 100 snapshots and a spatial resolution of 10 Hz is used (thus leading to B = 1991 frequency bins). The results of conventional beamforming, BF, IBF and BIBF are compared in Fig. 16. It is again checked that beamforming and BF both suffer from 1) a poor spatial resolution which makes the localization of the sources difficult below 2 kHz and 2) a limited dynamic range which prevents the two smallest sources from being identified. On the one hand, whereas IBF greatly upgrades the results, it is limited downward to 1 kHz for the reconstruction of the strongest source and to 2 kHz for the smallest source. The identification of the two small sources is actually intermittent and the smallest one is missed at most frequencies. On the other hand, BIBF completely fixes these errors. The accuracy in terms of localization and quantification is excellent over the full frequency range and for all sources. This is because a frequency-independent AF is better estimated from an "average" of the map in Fig. 16(c), which in turns improves the estimation of the sources at low frequencies and low SNRs.

Figure 17 displays the estimated hyperparameters α 2 and η 2 as a function of frequency. The prior source power spectrum α 2 ( f ) is constant on the average, as expected from sources driven by white noises. The spectrum of the regularization parameter, η 2 ( f ), shows a hyperbolic decrease because the conditioning of the inverse problem actually improves with frequency. The fact that theses hyperparameters have to estimated repeatedly for all frequencies in the band B justifies the need for the automatic process based on Bayesian regularization given in Appendix B.

Concerning computational time, the BIBF algorithm was actually found faster than repeating the IBF for all frequencies in the scrutinized band. As previoulsy explained, this is because accounting for all frequencies at once makes the convergence of the iterations faster (convergence was achieved after 7 iterations for BIBF while each run of the IBF required about twice as many iterations on the average).

Experimental application

Experiments have been conducted to validate the theoretical results of the previous sections. The experimental setup is displayed in Fig. 18. It consists of an 117 microphone LMS Sound Camera Digital Array placed in a semianechoic room, 30 cm away from "point" sources flush-mounted on a wood plate and spaced apart by a distance of 14 cm. The sources are fed by mutually uncorrelated white-noise signals with a global level difference of about 12 dB. The radiation of each source with the other one switched off was first recorded by an intensity probe in order to provide a reference based on a normative sound power measurement technique.

Since the radiation of the acoustic energy is in the half space, the baffled-source Green functions were used in the inverse problem. In addition, since the position of the sources is expected to be invariant in frequency, the experimental setup provides an excellent example to test the BIBF algorithm. Two priors are chosen to enforce sparsity of the reconstructed sources: i) the MCS prior (with a = b = 0.01 =) because it has been previously found to enforce high sparsity while being simple to use and ii) the GMCG prior with p = 1.3 because it best matches the actual sources which, in this experience, do not reduce to perfect points. The results are displayed by means of the sound power, W(r q ), calculated at each node q of the calculation grid with the formula

f (kHz) f (kHz) a) b) α 2 η 2
W(r q ) = ρc 0 k 2 8π S qq - ρc 0 k 2 l q Im G(r q |r l )S ql , (32)
where operator Im{• • • } extracts the imagy part, k stands for the wavenumber, ρ for the air density, c 0 for the sound speed, and S ql is the estimated cross-power spectrum between the complex volume velocities at positions r q and r l [START_REF] Pereira | Acoustic imaging in enclosed spaces[END_REF].

Examples of sound power maps reconstructed in the third-octave band at 1 kHz by IBF and BIBF and using the two aforementoned priors are shown in Fig. 19(a). It is seen that the benefit of using BIBF instead of IBF is to improve the spatial resolution and thus to better reveal the presence of the smallest source. This is further evidenced by the one dimensional cut of the sound power distribution at y = 0 shown in Fig. 19(b), where the very low frequency 300 Hz is now considered. Besides, it is observed in Fig. 19 that the MCS prior returns a sparser reconstruction than the GMCG prior with p = 1.3. Although this may be advantageous for localization, it does not necessarily reflects to the actual source characteristics as shown hereafter.

Figure 20(a) displays the sound power spectra estimated from BF and IBF with the two different priors. The sound power was evaluated from Eq. ( 32) by integrating over the large circular region depicted in Fig. 19(a), which encloses the two sources. The spectra are also compared with the reference returned by the intensity probe, which imposes the upper limit of 6.5 kHz of the frequency interval. On the one hand, it is seen that BF underestimates the actual sound power, especially in the higher frequency range. This complies with previous numerical simulations, even though the under-estimation is not as important as before since the array now covers a large solid angle from the source point of view. On the other hand, the MCS prior is found to overestimate the sound power. This is not contradictory with the results of the paper since the actual sources used in the experiment are not true monopoles: their supports are not exactly point-wise on surface Γ and their directivities are not exactly omnidirectional. Consequently, the GMCG prior with a controlled degree of sparsity (p = 1.3) seems to better reflect their actual nature. This observation points out the importance of correctly choosing the prior type and, therefore, of having a general algorithm that allows testing several of them.

Finally, Fig. 20(b) displays the estimated sound power of each source compared with the intensity reference obtained with one of the two soures alternatively switched off. For the secondary source this is achieved by integrating over the small circular region depicted in Fig. 19(a) and, for the primary source, by taking the difference with the spectra of Fig. 20(a). Only the results obtained with the GMCG prior are shown, yet compared for IBF and BIBF. It is again checked that BIBF allows a better resolution of the secondary source in the very low frequency range, down to 300 Hz, which is compliant with the result shown in Fig. 19(b).

Conclusion

This paper has demonstrated how sparse acoustical holography is naturally achieved by simply iterating the Bayesian Focusing method previously proposed in Ref. [START_REF] Antoni | A bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF]. The principle consists in using the estimated sources in the current iteration to update the aperture function in the next iteration. Formally, this is equivalent to considering the aperture function as a random quantity (endowed with a prior probability density function) that is estimated conjointly with the sources. Not only does this point of view provide physical insight into the mechanism of sparse acoustical holography -and in particular why it can considerably improve the estimation of source quantification and directivity high in frequency -but it also allows direct extension to different scenarios. For instance, the method is easily formulated in terms of the cross-spectral matrix and adapted to the broadband case (i.e. frequency group sparsity) where source positions are assumed stationary in space. The method also allows complete flexibility in the choice of the sparsity penalty. Several updating rules for the aperture function have been investigated, rooted on different source priors. The multivariate complex Student-t prior appears to be a good candidate, which returns sparser results than the typical 1 -norm approach while at the time offering stable results. However, the Generalized multivariate complex Gaussian may be advantageous in cases where the degree of sparisty needs to be controlled, for instance when the actual sources are not truly pointwise. A last but not least advantage is that Iterated Bayesian Focusing inherits the possibility of automatic setting of the regularization parameter, an issue which has remained critical in sparse methods. Overall, it is believed that Iterated Bayesian Focusing provides a comprehensive view of sparse acoustical holography that might not be shared by other approaches.

Beyond methodological aspects, this paper has also confirmed the experimental results of Ref. [START_REF] Hald | Fast wideband acoustical holography[END_REF] establishing that sparse holography improves the reconstruction of sources not only in terms of localization -as does compressive beamforming in the far-field -but also of quantification and of directivity, in a frequency range that is considerably enlarged as compared to classical methods. It has been found that these improvements can be achieved even with regular arrays, provided that sparser priors than those leading to the standard 1 penalty are used; therefore, using random arrays is not a necessary condition as is in compressed beamforming or ESM. correlation

E s i s H i = ΦE c i c H i Φ H = ΦΩ c Φ H (A.2)
where Ω c stands as the prior covariance matrix of c i . Therefore, by assigning to c i a prior PDF in the form of the complex Gaussian,

[c i ] = N C (c i ; 0, Ω c ), one finds the MAP estimate ĉi = Ω c Φ H G H GΦΩ c Φ H G H + β 2 I -1 p i (A.3)
with η 2 as defined in Eq. ( 9). Equation (A.3) generally applies for any choice of Φ and Ω c . In the case of sparse sources, it make sense to expect a priori an incoherent source field, yet possibly non-homogeneous, as given by Eq. [START_REF] Brooks | A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays[END_REF]. It is proved in Ref. [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF] that the prior covariance matrix Ω c which makes E s i s H i in Eq. (A.2) as close as possible to α 2 Σ 0 -with Σ 2 0 diagonal as defined in Eq. ( 12) -in the mean-square sense is

Ω c = α 2 Φ H Σ -2 0 Φ -1 . (A.4)
With such a choice, it is further proved in Ref. [START_REF] Antoni | A bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF] that the optimal Φ that minimizes the mean-square error (or returns the highest MAP) of the estimated source distribution is

Φ = Σ 2 0 G H ∈ C Q×M . (A.5)
This shows in passing that for Φ to be a basis, one must have M = Q. Substituting expressions (A.4) and (A.5) back into Eqs. (A.1) and (A.3) then returns Eq. ( 9). Comparison of Eqs. (A.3) and ( 9) shows that, with a complex Gaussian prior PDF N C (0, α 2 Σ 2 0 ), the optimal parameterization that enforces an incoherent source field is actually equivalent to using the canonical basis I Q×Q together with the prior covariance matrix α 2 Σ 2 0 on the spatial samples s i of the source seen as the unknowns of the problem. The same parameterization can actually be adopted afterwards with Q > M, provided that Eq. ( 9) remains invertible by properly setting the value of the regularization parameter η 2 .

Although this parameterization turns out somewhat similar to the Equivalent Source Method (ESM) [START_REF] Pereira | A theoretical and experimental comparison of the equivalent source method and a bayesian approach to noise source identification[END_REF][START_REF] Fernandez-Grande | A sparse equivalent source method for near-field acoustic holography[END_REF], it bears several differences which are worth being highlighted.

• First, the present proof is deductive and it establishes that the parameterization of the problem is optimal under the assumption of prior incoherence of the source distribution. Other assumptions on the prior spatial correlation would lead to other parameterizations.

• Second, it emphasizes that the actual dimension of the problem (i.e. the dimension of the spatial basis Φ as given in Eq. (A.5)) is K = M (the number of microphones) and not Q (the dimension of the discretized source surface), which is advantagous since one often has M << Q. For this reason, the optimal basis introduced in Eq. (A.5) leads to a minimal parameterization. This is not explicit in the ESM method.

• Third, it strictly involves Green functions satisfying Neumann boundary conditions on Γ s and not necessarily free-field Green functions (of monopoles) as is commonplace with the ESM.

Appendix B. Bayesian regularization

Following the Bayesian program, the regularization parameter appearing in Eqs. ( 9) together with the source power needed in Eq. ( 23) are estimated as the MAP of the joint PDF [α 2 , η 2 |S pp , Σ 2 ] with S pp the CSM (see Eq. ( 1)) and Σ 2 a given AF. This subsection resumes the results of Ref. [START_REF] Pereira | Empirical bayesian regularization of the inverse acoustic problem[END_REF] where the solution to this problem is addressed in details. The algorithm is as follows.

• Step 1: Compute the eigenvalue decomposition

GΣ 2 G H = M m=1 s 2 m u m u H m (B.1)
where u k are and s 2 k , k = 1, ..., M, are the eigenvectors and eigenvalues of GΣ 2 G H , respectively.

• Step 2: Project the data in the space spanned by the eigenvectors:

< |y k | 2 > . = u H k S pp u k , k = 1, ..., M. (B.2)
• Step 3: Estimate the source power as

α2 (η 2 ) = 1 M M k=1 < |y k | 2 > s 2 k + η 2 .
(B.3)

• Step 4: Estimate the regularization parameter as the minimum of the cost function

J(η 2 ) = 1 M M k=1 ln(s 2 k + η 2 ) + M - 2 N ln α2 (η 2 ), (B.4)
which may be easily found by a simple grid search and eventually substitute back this result into Eq. (B.3).

Appendix C. Prior PDFs in the SMoG family

This appendix provides expressions of some prior PDFs in the SMoG family obtained from different choices of the hyperprior τ 2 q in Eq. ( 18). Expressions are given in terms of the sufficient statistic χ 2 q .

= NS qq /α 2 σ 2 0,q .

• Multivariate complex Student-t (MCS). Let us assume that τ 2 q is distributed like an inverse Gamma, with shape and scale parameters a > 0 and b > 0, respectively:

τ 2 q = b a Γ(a) e -bτ -2 q τ 2(a+1) q . (C.1)
It then holds that

χ 2 q = Γ(N + a)b a (πα 2 σ 2 0,q ) N Γ(a) χ 2 q + b (N+a) , (C.2)
which is recognized as a multivariate complex version of the Student-t distribution. This PDF is made very leptokurtic by setting the value of a small. By substituting the MCS for the prior PDF into the MAP (see Eq. (E.1) in appendix Appendix E), one arrives at an LS formulation with a logarithmic penalty, viz ŝi = Argmax

s i [P|S][S] = Argmin s i (-ln[p i |s i ][S]) = Argmin s i         ||p i -Gs i || 2 2 + β 2 (N + a) ln         α -2 N j=1 s H j Σ -2 0 s j + b                 . (C.3)
The logarithmic penalty is here naturally endowed with a regularization term b > 0 that prevents it from diverging.

• Multivariate complex K-distribution (MCK). let us assume that τ 2 q is distributed like a Gamma, with shape and rate parameters a > 0 and b > 0, respectively:

τ 2 q = b a Γ(a) e -bτ 2 q τ 2(a-1) q . (C.4)
It then holds that

χ 2 q = 2b (N+a) 2 (πα 2 σ 2 0,q ) N Γ(a) χ 2 q (a-N) 2 K N-a 2 bχ 2 q , (C.5)
where K ν (z) is the modified Bessel function of the second kind with order ν (note that K ν (z) = K -ν (z)). This is recognized as a multivariate complex version of the K-distribution [START_REF] Eltoft | Multivariate Scale Mixture of Gaussians Modeling[END_REF]. This PDF is made very leptokurtic by setting the value of a small. The MCK also leads to an equivalent penalized LS problem with a flexible penalty function that can be given various shapes depending on the values of a, b and N.

• Multivariate complex Laplacian (MCL). As a special case when a = N + 1/2 and b = 1/4, the MCK becomes

χ 2 q = 1/2 (πα 2 σ 2 0,q ) N Γ(N) Γ(2N) e - √ χ 2 q (C.6)
which is a multivariate complex version of the Laplacian distribution [START_REF] Eltoft | Multivariate Scale Mixture of Gaussians Modeling[END_REF] as well as a particular case of the generalized multivariate complex Gaussian introduced in the next paragraph. The Laplacian PDF is a popular choice to enforce sparsity of real-value coefficients since it corresponds to a cost function with a penalty term based on the 1 -norm. This is exploited in the Bayesian LASSO [START_REF] Park | The Bayesian Lasso[END_REF][START_REF] Hans | Bayesian lasso regression[END_REF]. It is noteworthy that in the complex case, the Laplacian originates from an SMoG rooted on a different prior than in the real case (i.e. an exponential PDF). The equivalent penalized LS problem reads ŝi = Argmin

s i           ||p i -Gs i || 2 2 + β 2 α         N j=1 s H j Σ -2 0 s j         1 2           (C.7)
where the penalty function happens to mix the 1 and the 2 norms in a way that may not be obvious without using the Bayesian approach.

• Generalized multivariate complex Gaussian (GMCG). The multivariate complex version of a generalized Gaussian PDF reads

χ 2 q = p/2 (πα 2 σ 2 0,q ) N Γ(N) Γ(2N/p) e -(χ 2 q ) p 2 , 0 < p ≤ 2. (C.8)
Since the generalized Gaussian (also known as the exponential power distribution in the statistical literature) is an SMoG, the GMCG is also an SMoG (actually for p ≤ 2). The expression of τ 2 q happens to be quite intricate, yet only the expression of the marginal PDF (C.8) will be needed in the following. Formula (C.8) extends the result of Ref. [START_REF] Novey | A complex generalized gaussian distribution -characterization, generation, and estimation[END_REF] to the multivariate case.

Sparsity is enforced by setting small values of p, strictly smaller than 2 and typically smaller than or equal to 1. As special cases, p = 2 corresponds to the multivariate complex Gaussian (MCG) and p = 1 to the MCL.

The GMCG prior gives rise to the penalized LS problem, ŝi = Argmin

s i           ||p i -Gs i || 2 2 + β 2 α p         N j=1 s H j Σ -2 0 s j         p 2           , (C.9)
where the penalty function mixes the 2 and the p norms.

• Mixture of multivariate complex Gaussians (MMCG). Mixture of multivariate complex Gaussians are obtained from the SMoG model when τ 2 q is a discrete random variable. The simplest situation is when it takes two values, say τ2 1 and τ2 2 , with probabilities π and (1π), respectively, which means that a fraction π of the coefficients are distributed like a complex Gaussian with variance τ2

1 and a fraction 1π like a complex Gaussian with variance τ2

2 . The PDF then reads τ 2 q = πδ(τ 2 q -τ2 1 ) + (1π)δ(τ 2 q -τ2 2 ) (with δ the Dirac delta) and the corresponding SMoG is χ 2 q = π • N C,N (χ 2 q,[k-1] ; 0, πα 2 σ 2 0,q τ2 1 ) + (1π) • N C,N (χ 2 q,[k-1] ; 0, πα 2 σ 2 0,q τ2 2 ) (C.10)

where N C,N (χ 2 q,[k-1] ; 0, σ 2 ) . = σ -2N e -χ 2 q /σ 2 . Very sparse PDF can thus be constructed when τ1 << τ2 and π → 1. The extreme scenario when τ1 → 0 returns the Bernouilli-Gauss model where a fraction π of the coefficients are exactly zero on the average, the other coefficients being distributed like a complex Gaussian with variance τ2

2 . The equivalence with a penalized LS problem is here less obvious.

Appendix E. Convergence of the IBF algorithm

Let us first assume that the regularization parameter η 2 and the prior source power α 2 are known. Then, considering all snapshots together, the Bayesian program is to find the MAP estimate ŝi = Argmax The second line results from the application of Bayes rule [START_REF] Simard | Acoustic source identification: Experimenting the 1 minimization approach[END_REF] where the evidence [P] has been ignored since it does not depend on the unknown quantity ŝi . Next, replacing [P|S] by its multivariate complex Gaussian expression and using Eq. ( 18) to express [S] as an SMoG, it comes ŝi = Argmax

s i          -β -2 ||p i -Gs i || 2 + ln Diag          E τ 2 q          N j=1 s q, j |τ 2 q [τ 2 q ]                            (E.1)
where only functions of the unknown quantity ŝi have been kept in the expression to be maximized, the diagonal operator applies with respect to index q and E τ 2 q means that the expected value is taken with respect to the random variable τ 2 q . At this stage, it is noticed that if the order of the expectation operator and of the logarithm were reversed in the last term of Eq. (E.1), the cost function would be quadratic in ŝi and its optimization would considerably simplify. This is the essence of the EM algorithm, which solves problem (E.1) by iterating between the following two steps [START_REF] Mclachlan | The EM algorithm and extensions, 2nd Edition[END_REF].

E step: Calculate Q(S; S [k-1] ) = E τ 2 q ln N j=1 s q, j |τ 2 q |P, S [k-1] where the expectation is now conditioned to the data p j and the previous estimate of the source coefficients s j,[k-1] for all j = 1, ..., N. = -E{τ -2 q |P, S [k-1] } N j=1 |s q, j | 2 α 2 σ 2 q,0 -E{ln(τ 2 q )|P, S [k-1] } -N ln(πα 2 σ 2 q,0 )

where [s q, j |τ 2 q ] = N C (s q, j ; 0, α 2 σ 2 q,0 , τ 2 q ) has been used (see Eq. ( 18)). By ignoring all quantities which do not depend on s j,[k-1] and inverting the sign of the cost function, the EM algorithm takes the simplified form:

E step: Calculate τ -2 q,[k] = E{τ -2 q |P, S [k-1] } (E.2)
which, in the SMoG, family takes the simplified form E{τ -2 q |χ 2 q,[k-1] } where conditioning is constricted to the normalized sum of squares χ 2 q,[k-1] defined in Eq. ( 19) and evaluated at the previous iteration k -1. M step: Calculate

s i,[k] = Argmin β -2 ||p i -Gs i || 2 + α -2 s H i Diag(σ -2 q,0 τ 2 q , [k] -2 )s i = G H G + η 2 Diag(σ -2 q,0 τ -2 q,[k] ) -1
G H p i (E.3a) = Diag(σ 2 q,0 τ 2 q,[k] )G H (GDiag(σ 2 q,0 τ 2 q,[k] )G H + η 2 I) -1 p i . (E.3b)

This last equation is recognized as being equivalent to Eq. ( 9). Eventually, it is found that source coefficients can

• Multivariate complex Student-t (MCS). The expected relative intensity related to the MCS prior PDF is easily found by application of Eq. ( 24) to Eq. (C.2):

τ2 q,[k] = χ 2 q,[k-1] + b N + a . (G.1)
As seen in Fig. 11(c-d), the estimated relative intensity τ2 q,[k] is a linearly increasing function of the normalized sum of squares χ 2 q,[k-1] , which means that source coefficients with high values are progressively emphasized by the updating rule [START_REF] Abusag | Near-field acoustic holography using the method of superposition[END_REF] in the IBF algorithm. At convergence, only the coefficients with largest moduli will survive. Very sparse solutions are expected when small values are quickly attenuated, which requires b to be small, and large values are amplified, which requires a to be small as well.

• Multivariate complex K-distribution (MCK). Using the properties of the modified Bessel functions, the expected relative intensity related to the MCK is found as

τ2 q,[k] = χ 2 q,[k-1] N -a             1 + bχ 2 q,[k-1] N -a K N-a-1 2 bχ 2 q,[k-1] K N-a 2 bχ 2 q,[k-1]             -1
, N > a, (G.2)

with asymptotic behavior

τ2 q,[k] ∼ χ 2 q,[k-1] N -a , N > a, for small χ 2 q,[k-1] (G.3) and τ2 q,[k] ∼ χ 2 q,[k-1] N -a             1 + bχ 2 q,[k-1] N -a             -1
, N > a, for large χ 2 q,[k-1] .

(G.4) This is a crescent function of the normalized sum of squares χ 2 q,[k-1] , yet with a nonlinear behavior intermediate between the linear and the square root case. It is noteworthy that for large χ 2 q,[k-1] the MCK prior behaves like the MCL given hereafter (see Eq. (G.6)). In practice, the numerical evaluation of the modified Bessel functions may be troublesome at high orders and for small values of χ 2 q,[k-1] . A good approximation is proposed for large orders Na is (see Appendix Appendix H)

τ2 q,[k] ≈ χ 2 q,[k-1] N -a             1 + bχ 2 q,[k-1] N -a        1 -e - √ bχ 2 q,[k-1] N-a-1                    -1
, N > a + 1.

(G.5)

• Generalized multivariate complex Gaussian (GMCG). The expected relative intensity related to the GMCG is readily found as

τ2 q,[k] = 2 p χ 2 q,[k-1] 1-p 2 , 0 < p ≤ 2. (G.6)
This is again a monotonic power function of χ 2 q,[k-1] , where high sparsity can be reached by setting p small. The special case p = 1 returns a square root behavior similar to the asymptote of the MCK. One definite advantage of the GMCG is to involve only one parameter p.

The updating rule of the GMCG is particularly intuitive since it reads, in logarithmic scale, ln σ 2 q,[k] = ln σ 2 q,0 + ln τ2 In accordance with Ref. [START_REF] Mclachlan | The EM algorithm and extensions, 2nd Edition[END_REF], let define the observed information matrix of the incomplete data (i.e. without introducing τ 2 q ),

I i = - ∂ 2 ∂s i ∂s H i ln [S|P] ,
(J.1)

and the expected information matrix of the complete data (i.e. marginalized over τ = {τ 2 q , q = 1, .., Q})

I c,i = -E τ 2
∂ 2 ∂s i ∂s H i ln S|P, τ 2 |P, S , (J.2) both conditioned to the observations P and evaluated for snapshot i. The missing information matrix is then defined as I m,i = I c,i -I i , (J.3) which measures information deficiency when the relative intensity τ 2 q is not directly observed. As explained in Ref. [START_REF] Mclachlan | The EM algorithm and extensions, 2nd Edition[END_REF], an estimate of the covariance matrix of the estimated source coefficients is returned by the inverse of the observed information matrix I i = I c,i -I m,i averaged over all snapshots, whereas the speed of convergence of the EM algorithm is controlled by the largest eigenvalue of J i = I -1 c,i I m,i , both evaluated at the MAP solution ŝi . In both cases, the missing information matrix is seen to act as the main leverage and its norm should therefore be minimized in order to achieve good convergence of the estimates.

The missing information matrix is calculated hereafter from Eq. (J.2) after evaluating the expressions for I c,i and I i . Since S|P, τ 2 q ∝ [P|S] S|τ 2 q where [P|S] and S|τ 2 q are complex Gaussians PDFs as given by Eqs. ( 7) and ( 8), it readily comes that I c,i = β -2 G H G + α -2 Diag σ -2 0,q E{τ -2 q |S} (J.4)

where E{τ -2 q |S} = -∂ ∂(χ 2 q ) ln χ 2 q from Eq. (K.4). Similarly, since [S|P] ∝ [P|S] [S],

I i = β -2 G H G - ∂ 2 ∂s i ∂s H i ln [S] = β -2 G H G -α -2 Diag        σ -2 0,q        ∂ ∂(χ 2 q ) + |s q,i | 2 α 2 σ 2 0,q ∂ 2 (∂χ 2 q ) 2        ln χ 2 q        (J.5)
where [S] is given by Eq. ( 18) and the chain rule, ∂ ∂s q,i = s * q,i α 2 σ 2 0,q ∂ ∂(χ 2 q ) with χ 2 q as defined in Eq. ( 23), has been used. Therefore,

I m,i = Diag        |s q,i | 2 (α 2 σ 2 0,q ) 2 ∂ 2 (∂χ 2 q ) 2 ln χ 2 q        , (J.6)
whose average over snapshot i = 1, ..., N and evaluation at the MAP solution returns the first equality in [START_REF] Oudompheng | A theoretical and experimental comparison of the iterative equivalent source method and the generalized inverse beamforming[END_REF]. Now, considering that [χ 2 q ] ∝ e -χ 2 q /τ 2 q [τ 2 q ]dτ 2 q , it can be checked that ∂ 2 (∂χ 2 q ) 2 ln χ 2 q = ∞ 0 τ -4 q τ 2 q |{s q, j,[k-1] } N j=1 τ 2 q dτ 2 q -∞ 0 τ -4 q τ 2 q |{s q, j,[k-1] } N j=1 τ 2 q dτ 2 q 2 = Var τ -2 q |S , (J.7) which proves the second equality of Eq. [START_REF] Oudompheng | A theoretical and experimental comparison of the iterative equivalent source method and the generalized inverse beamforming[END_REF].

The above results lead to the following particular cases:

• The diagonal elements of the information matrix of the MCS prior are upper bounded provided that b > 0:

[I m ] MCS qq = 1 α 2 σ 2 0,q (N + a) N χ 2 q (χ 2 q + b) 2 ≤ 1 α 2 σ 2 0,q (N + a) N . (J.8)

Figure 1 :Figure 2 :

 12 Figure 1: Scheme of a simple inverse problem comprising an ULA of microphones (black circles) and a monopole source (red dot) on a linear source domain Γ s (black segment).

Figure 3 :

 3 Figure3: a) b) Real part of the pressure field radiated by one microphone in the array (here the central one) as if it was a source ( f = 6 kHz). The black (resp. grey) curve is the real (resp. imaginary) part of the pressure radiated on Γ s as a function of r. b) Real and imaginary parts (black and grey curves, respectively) of the basis functions (as obtained by the mechanism illustrated in (a)) rooted on an uniform aperture function and c) the narrow aperture function of Fig.2(c). The aperture function used in (c) is displayed by the black dotted curves.

Figure 4 :

 4 Figure 4: Amplitudes of spatial spectra of the M basis functions of a) Fig. 3(a) and b) Fig. 3(b). Note that all amplitude spectra are superimposed in (a).

Figure 5 :

 5 Figure 5: Principle of Iterated Bayesian Focusing: the reconstructed source ŝ[k] (r) at the current iteration k is used to refine the aperture function σ 2 [k+1] (r) used in the next iteration.

shape parameter a > 0 scale parameter b > 0 4 (

 04 Multivariate complex K-distribution (MCK) particular case: Multivariate complex Laplace (MCL) when a = N + 1 2 and b = 1 Not required) power 0 < p ≤ 2 Generalized multivariate complex Gaussian (GMCG) particular cases: Multivariate complex Gaussian (MCG) when p = 2 Multivariate complex Laplace (MCL) when p = complex Gaussians (MMCG)

2 Figure 8 :

 28 Figure 8: Directed graph of the hierarchical model for inferring S. Unknown parameters to infer are shown in circles and known quantities in rectangles. Note that P and S can be replaced by the CSM S pp and the source power spectrum S qq .

Figure 9 :

 9 Figure 9: Update of the aperture function σ 2q,[k] with respect to iteration numbers k = 0, ..., 9 for the MCS prior. The zeroth iteration corresponds to the initial aperture function σ 2 q,0 . Configurations with a-b) one monopole at r = 0 m and c-d) two incoherent monopoles at r = 0 m and r = -0.05 m (shown by red dots) with relative levels in the ratio 0.04. The two columns relate to linear and dB scales, respectively.

Figure 10 :

 10 Figure10: Spectra of the acoustic power (relatively to the true value) radiated by a monopole source estimated with Bayesian Focusing (BF) (blue curve) based on the large aperture function σ 2 0 shown in Fig.9and with Iterated Bayesian Focusing (IBF) (black curve). The shaded areas indicate 3-sigma confidence intervals and the dotted curves the coefficients of variation of BF (blue) and IBF (black) amplified by a factor 10. The vertical dotted line corresponds to the cutting frequency f max ∼ c 0 /d. Results are averaged over 100 draws of the independent complex Gaussian noise with 20 dB SNR. (CV = coefficient of variation.)

Figure 11 :

 11 Figure 11: a) Prior PDF [χ 2q ] for N = 1 (in logarithmic scale) and b) for N = 50 in the MCS case with (a = 0.1, b = 0.2) (black line), the MCK case with (a = 3N/4, b = 1) (red dotted line), the MCL case (yellow line) and the MMCG case with (π = 0.9, τ2 1 = 0.0001, τ2 2 = 1) (purple line) compared with the multivariate complex Gaussian (MCG) case (dotted black line). Sparse behaviors are essentially recognized by the distance to the Gaussian at the origin (χ 2 q = 0) and towards infinity (χ 2 q → ∞). c) Corresponding relative intensity τ2q = E{τ -2 q | Ŝ} -1 = -(∂ ln([χ 2 q ]/∂(χ 2 q )) -1 for N = 1 and d) for N = 50, all normalized to 1 at χ 2 q /N = 1. e) Corresponding variance Var{τ -2 q | Ŝ} = ∂ 2 ln[χ 2 q ]/(∂χ 2 q ) 2 for N = 1 and f) for N = 50 with excerpt in the range 0 ≤ χ 2 q /N ≤ 0.06.

Figure 12 :

 12 Figure 12: Estimated source power spectrum S qq at f = 1kHz with N = 1 snapshot in the case of two incoherent monopoles with equal strengths (indicated by blue circles). Results from four independent snapshot realizations are superimposed. The positions of the microphones is marked by black circles in (a). a) Results returned by conventional beamforming (blue) and Bayesian Focusing (red). Results returned by Iterated Bayesian Focusing with b) the MCL prior (GMCG with p = 1), c) the GMCG prior with p = 0, d) the MCS prior (a = b = 0.01), e) the MCK prior (a = b = 0.1) and f) the MMCG prior with π = 0.99, τ2 2 = 1, τ2 1 = 0.01 (red) and τ2 1 = 0.001 (orange). The red/orange curves in (b-f) relate to the source power spectrum S qq and the black/gray curves to the updated aperture function σ 2 q .

Figure 13 :Figure 14 :

 1314 Figure 13: Same configuration as in Fig. 12 with N = 100 snapshots.

Figure 15 :

 15 Figure 15: Same configuration as in Fig. 14 at f = 10 kHz.

Figure 16 :

 16 Figure 16: Estimated source distribution on the 1D domain Γ s = {r : -0.2 ≤ r ≤ 0.3}, as a function of frequency, when the theoretical distribution comprises four monopoles placed at r = 0.2 m, r = 0.1 m, r = 0 m and r = -0.1 m, driven by uncorrelated white Gaussian signals with levels in the proportion 1, 0.25, 0.04 and 0.01, respectively, and SNR of 20 dB. Estimates returned by a) BF with a fixed uniform aperture function, b) conventional beamforming, c) IBF with the MCS prior, and d) broadband IBF with the MCS prior.

Figure 17 :Figure 18 :

 1718 Figure 17: Estimated hyperparameters a) α 2 and b) η 2 , as a function of frequency, returned by broadband IBF with the MCS prior.

Figure 19 :

 19 Figure 19: Sound power maps reconstructed in the 1 kHz third-octave band from iterated Bayesian Focusing using a) the GMCG prior (with p = 1.3) and b) the MCS prior (with a = b = 0.01). Reconstruction from Broadband Iterated Bayesian Focusing MCS using c) the GMCG prior and d) the MCS prior. e) Comparison of the sound power distribution (|x| ≤ 0.5m, y = 0) at 300 Hz estimated from Iterated Bayesian Focusing and Broadband Iterated Bayesian Focusing with GMCG (p = 1.3) and MCS (a = b = 0.01) priors.

2 frequencyFigure 20

 220 Figure 20: a) Overall sound power spectra estimated with Iterated Bayesian Focusing, with the MCS (a = b = 0.01) and GMCG (p = 1.3) priors, compared to Bayesian Focusing and to the intensity probe reference. b) Sound power spectra of each source estimated with Iterated Bayesian Focusing and Broadband Iterated Bayesian Focusing, with the GMCG (p = 1.3) prior.

  P|S] + ln [S]) .

2 q

 2 M step: Calculate s i,[k] = Argmax s i -β -2 ||p i -Gs i || 2 + Q(S; S [k-1] ) .The E step is further worked out asQ(S; S [k-1] ) = ln s q, j |τ 2 q |P, S [k-1]

  q,0 + 1 -p 2 ln S qq,[k-1] + C (G.7)

Table 1 :

 1 Prior PDFs constructed in the SMoG family Hyperprior PDF τ 2

	q

Table 2 :

 2 Hyperparameters in the SMoG family

	Sparse prior PDF	Number of relevant hyperparameters Setting
	Multivariate complex Laplace (MCL)	0
	Generalized multivariate complex Gaussian (GMCG) 1	0 < p ≤ 2
	Multivariate complex Student-t (MCS)	2 reduced to 1	a = 0 and b small or a = b small
	Multivariate complex K-distribution (MCK)	2	a and b
	Mixture of multivariate complex Gaussians (MMCG) 3 reduced to 2	π and τ2 1 /τ 2 2

Note that sparse Bayesian learning should be undertood as a general terminology that embodies the family of Bayesian methods dedicated to enforcing sparsity, yet it is not an algorithm by itself.

Acknowledgments

This work was performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007), and supported by the Labcom P3A (ANR-13-LAB2-0011-01).

Appendix A. Optimal basis for sparse sources

In a general setting and in functional notation, the source distribution is expanded onto a spatial basis {φ k (r)} Q k=1 , M ≤ Q, i.e.

where coefficients c k,i , k = 1, ..., Q, i = 1, ..., N, are the unknowns of the problem. In discrete form, this reads s i = Φc i where [Φ] kq . = φ k (r q ) and vector c i collects coefficients c k,i . This parametrization imposes a priori the spatial Appendix D. Point estimates in the SMoG family of priors

It is proved in this Appendix that, in the SMoG family of priors, the MAP and the posterior mean estimates of s i both take the structure of the penalized LS solution given in Eq. [START_REF] Zavala | Generalized inverse beamforming with optimized regularization strategy[END_REF].

Appendix D.1. Posterior mean

The posterior mean is defined as

where the notation Σ 2 . = Diag(σ 2 q ) has been used,

ds i is the posterior mean ŝi (Σ 2 ) conditioned on Σ 2 as given in Eq. ( 9), one has

MAP estimate

The MAP estimate satisfies the optimality condition ∂ ∂s H i

[S|P] = 0. Using Eq. (D.1), this is equivalent to

where [S|P, Σ 2 ][Σ 2 |P]dΣ 2 = [S|P] has been used. Dividing all terms by β -2 [S|P] and recognizing that

be estimated for each snapshot i, independently, yet by using estimates of η 2 and τ 2 q which are evaluated by using all available snapshots j = 1, ..., N.

Let us now consider the situation where the regularization parameter and the source power are unknown. Then, the goal is to estimate them jointly with the source coefficients, viz ŝi , α2 , η2 = Arg max

The maximization of the above cost function can be significantly simplified by alternating between the maximization of [S|α 2 , η 2 , P] and [α 2 , η 2 |{s j , p j } N j=1 ], which is the principle of the Expectation-Conditional Maximization (ECM) algorithm [START_REF] Mclachlan | The EM algorithm and extensions, 2nd Edition[END_REF]. The convergence of the ECM is further accelerated by replacing the maximization of [α 2 , η 2 |{s j , p j } N j=1 ] by that of its marginal over the source coefficients, that is

(where conditioning on P is eventually to be equivalent to conditioning on S pp ). This is known as "collapsing" [START_REF] Van Dyk | Partially collapsed gibbs samplers[END_REF]. This justifies interlacing the Bayesian regularization algorithm described in Appendix Appendix B with the EM steps.

Appendix F. Estimation of the relative intensity

The conditional expectation (E.2) is evaluated as follows

where it has been recognized that τ 2 q depends directly and only on s q, j,[k-1] . Using Bayes' rule,

Now, exploiting the fact that

in the SMoG model and setting

, one finds

Application of the latter formula to the prior PDFs (C.2), (C.5), (C.8) and (C.10) returns Eqs. (G.1), (G.2), (G.6) and (G.9), respectively. Incidentally, the above result also proves that τ -2 q,[k] depends on χ 2 q,[k-1] only, which justifies the contracted notation τ -2 q,[k] = E{τ -2 q |χ 2 q,[k-1] } instead of the full notation of (E.2).

Appendix G. Updating rules for the aperture function

This appendix provides closed-form expressions for τ2 q,[k] in the case of the prior PDFs introduced in section 4.2 (see Eqs. where C . = ln (2/p) + (1p/2) ln N/α 2 is a global scaling that does not depend on the local position r q . Therefore, the updating is interpreted as the weighted geometric mean between the initial AF σ 2 q,0 and the previous estimation of the source power spectrum S qq,[k-1] . The value of p specifies the balance between the two terms. For p = 1, uniform weights of value 1/2 are assigned to the two terms, whereas for 0 < p < 1, more weight is given to the data-dependent term.

• Multivariate complex Laplacian (MCL). Since the multivariate complex Laplacian is a particular case of the GMCG with p = 1, the expected relative intensity is

with updating rule ln σ 2 q,[k] = 1 2 ln σ 2 q,0 + 1 2 ln S qq,[k-1] + C. This is found to generalize the Bayesian LASSO [START_REF] Park | The Bayesian Lasso[END_REF] to the multivariate complex case.

• Mixture of multivariate complex Gaussians (MMCG). According to Eq. (C.10),

where

. Equation (G.9) has an insightful interpretation: τ-2 q is obtained as the weighted averaged of τ-2 1 and τ-2 2 where weights p 1 and p 2 reflect the likelihood of the data χ 2 q,[k-1] under a Gaussian model with relative intensities τ1 and τ2 , respectively. The MMCG case stands out from the other cases by producing a soft threshold, as shown in Fig. 11(c-d). Values of χ 2 q,[k-1] below a certain quantity γ c are attenuated by a factor nearly equal to τ2 1 whereas values above γ c are amplified by a factor nearly equal to τ2

2 . The value of the threshold may be assessed from the inflection point of τ2

q,[k] as a function of χ 2 q,[k-1] , viz

For large values of N, this tends to a hard threshold where all values below it would be exactly zeroed.

Appendix H. Proof of approximation (G.5)

For N > a + 1, the ratio of the modified Bessel functions in Eq. (G.2) is very well fitted by a function of the form

), where coefficient c = 1/(Na -1) is found so as to match the derivative at χ 2 q,[k-1] = 0.

Appendix I. Derivation of the bias -Eq. (29)

Setting p i = Gs i in Eq. ( 20) and using the matrix inversion lemma, one has

with Σ 2 as given in Eq. [START_REF] Chardon | Near-field acoustic holography using sparse regularization and compressive sampling principles[END_REF]. Now, setting

) in the former equality yields Eq. ( 29)

• The information matrix of the MMCG prior, [I m ] MMCG qq , is found upper bounded by 1

2

) N ).

• The information matrix of the GMCG prior with p < 2 (which includes the MCL) is seen to diverge at χ

(J.9)

• The MCK has unbounded variance,

, N > a, for small χ 2 q,[k-1] , (J.10) which diverges at χ 2 q = 0.

Appendix K. Broadband Iterated Bayesian focusing (BIBF)

The essential difference in BIBF is to assume that the AF is not a function of frequency. Therefore, the relative intensity τ 2 q should be estimated once for all frequencies. The equivalent of Eq. (F.1) then becomes

(with explicit dependence of the source coefficients s q,i,[k-1] ( f l ) on frequency f l ). Using Bayes' rule, this gives

q N i=1 B l=1 s q,i,[k-1]( f l ) |τ 2 q τ 2 q dτ 2 q ∞ 0 N i=1 B l=1 s q,i,[k-1]( f l ) |τ 2 q τ 2 q dτ 2 q (K.2)

where independence of the source coefficients have been assumed with respect to frequencies (which holds true under the assumption of a (time) stationary sound field). Since N i=1 B l=1 s q,i,[k-1]( f l ) |τ 2 q = 1 (α 2 σ 2 0,q τ 2 q ) NB e -B l=1 N i=1 |s q,i,[k-1]( f l ) | 2 τ 2 q α 2 σ 2 0,q (K.3) or, equivalently, χ 2 q,[k-1] ∝ e -χ 2 q,[k-1] with χ 2 q,[k-1] as defined in Eq. ( 31), one eventually finds