# DAMAS on one million points: Greedy coordinate descent for the Covariance Matrix Fitting method for acoustical source imaging

Gilles Chardon

► To cite this version:

Gilles Chardon. DAMAS on one million points: Greedy coordinate descent for the Covariance Matrix Fitting method for acoustical source imaging. Journal of Sound and Vibration, 2025, pp.119483. 10.1016/j.jsp.2025.119483. hal-05301250

HAL Id: hal-05301250

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

Submitted on 7 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

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

# DAMAS on one million points: Greedy coordinate descent for the Covariance Matrix Fitting method for acoustical source imaging

Gilles Chardon

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

Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France

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

## ARTICLE INFO

**Keywords:**

Source localization

Inverse problems

Deconvolution

## ABSTRACT

A variant of the DAMAS algorithm for acoustical source imaging by deconvolution is proposed, where the coefficient of the power map to be updated is chosen at each iteration, instead of being fixed in advance, or random. As the original DAMAS algorithm, the proposed algorithm converges to a solution of the Covariance Matrix Fitting problem. Experimental evaluation of the algorithm shows faster convergence than the DAMAS algorithm, with reduced memory usage (respectively 40 times faster, and 30 more efficient in memory in a 2D setting). The performances of the proposed method lie between CLEAN-SC, which performances are limited by its resolution, and the more performant but computationally intensive Lawson-Hanson algorithm and gridless method Sliding Frank-Wolfe. The computation complexity in space and time of the proposed method, both linear in function of the size of the discretization and number of microphones, allows imaging on a three dimensional grid of more than one million points, yielding sparser source distributions than with coarser grids in the case of point sources.

## 1. Introduction

The limitations of conventional beamforming for acoustical imaging, in particular in terms of resolution, are known [1], and several methods have been proposed with improved performances. Two such methods, DAMAS (Deconvolution Approach for the Mapping of Acoustical Sources) [2], based on a deconvolution approach, and CMF [3] (Covariance Matrix Fitting), based on a least-squares fit of the Cross Spectral Matrix (CSM), the Sample Covariance Matrix (SCM) of the measurements, have been shown to be equivalent [4], in the sense that the DAMAS algorithm converges, slowly, towards a solution to the optimization problem considered by the CMF method. In the same work, an algorithm based on the Lawson-Hanson algorithm (LH) [5] exploiting a factorization of the matrix involved in the problem was also proposed, with substantial acceleration of the convergence.

Several methods have been proposed since the introduction of DAMAS. In particular, CLEAN-SC [6] is based on an iterative algorithm aiming at “cleaning” the beamforming map. Gridless methods have been recently investigated, based on the Sliding Frank-Wolfe algorithm (SFW) [7] to solve a convex problem formulated in a measure space, or with global optimization [8] (e.g. with genetic algorithms [9]) operating on the parameters of the sources themselves (positions, powers), a nonconvex problem. Recent works also considered learning based methods, such as dense neural nets [10], learned ISTA (Iterated Soft Thresholding Algorithm) [11,12], or transformer based methods [13]. However the present paper will focus on methods that do not necessitate a learning phase.

Owing to its simplicity, DAMAS is still frequently used, as a baseline [10,12,14,15] or for actual use in acoustical imaging [16,17]. In this paper, a minor modification is applied to the DAMAS algorithm, making it competitive with the LH algorithm while avoiding the necessity of solving linear problems. Instead of selecting the coefficient of the power map to be updated by traversing the
grid in a predetermined order, or randomly, the proposed algorithm DAMAS-MI (MI for Maximal Improvement [18]) selects the coefficient to update in order to maximize the improvement of the objective function of the CMF problem. A careful implementation of the algorithm makes this new selection rule efficient in terms of computational complexity, while accelerating the convergence of the algorithm, and reducing its memory demands. Comparison with the LH algorithm on experimental data shows that the DAMAS-MI algorithm reaches lower objective function values earlier that LH, allowing faster results when exact minimization of the objective function is not necessary. Connections with beamforming and CLEAN-PSF [6] are also outlined. While the main focus is the computational efficiency of the proposed algorithm, performances of the proposed DAMAS-MI method are compared to DAMAS, LH, CLEAN-SC and the gridfree method CMF-SFW, complementing the already available results on the performances of DAMAS or CMF [19–21]. Previous methods aiming at improving the efficiency of DAMAS or CMF were mostly concerned with the reduction of the dimension of the problem through screening (using the beamforming map [22], wavelets [23], or neural networks [15]), leaving the algorithms unchanged.

The paper is structured as follows. The model considered in the CMF method, its optimization problem, and existing algorithms to solve it are discussed in Section 2. The proposed DAMAS-MI algorithm is introduced in Section 3, with discussion of its computational complexity and connections to relevant methods. Numerical and experimental results are given in Section 4, with comparisons to CLEAN-SC and CMF-SFW, and an application to a three dimensional region discretized with more that one million points. Concluding remarks are given in Section 5.

## 2. The covariance matrix fitting method

In this section, the signal model, the CMF method (also known as Spectral Estimation Method, SEM [24]) optimization problem, and algorithms from the state of the art are recalled.

### 2.1. Model

The sound radiated by a distribution of sources is measured on an array of M microphones. The sources are assumed to be monopolar and emitting uncorrelated noises. Their positions are assumed to be on a grid of N points in the domain of interest. A short-time Fourier transform is performed on the signals, yielding L snapshots $m_l$, vectors of complex amplitudes at a frequency of interest at time $l$ for each microphone.

These M dimensional vectors can be written as

$$m_l = \sum_{j=1}^{N} g_j s_{jl} + n_l, \qquad (1)$$

where $g_j$ is the vector of Green's functions from the jth source to the microphones, $n_l$ is an additive noise, assumed white in time and space and uncorrelated with the signals $s_{jl}$. Under these assumptions, the theoretical covariance matrix of the measurements $m_l$ is

$$R = \sum_{j=1}^{N} p_j g_j g_j^H + \sigma^2 I \qquad (2)$$

where $p_j$ is the powers of the jth source, $\sigma^2$ is the power of the noise and $I$ the M-dimensional identity matrix. With $G$ the M × N matrix of columns $g_j$ and $p$ the vector collecting the powers $p_j$, the covariance matrix $R$ can be rewritten

$$R = G \text{diag}(p) G^H + \sigma^2 I. \qquad (3)$$

In the CMF method, the powers of the sources are estimated by matching the above theoretical model to the sample covariance matrix, or cross spectral matrix, $C$, estimated from the L snapshots by

$$C = \frac{1}{L} \sum_{j=1}^{L} m_l m_l^H. \qquad (4)$$

More precisely, the powers $p$ are estimated by solving the non-negative least-squares (NNLS) problem

$$\underset{p \in \mathbb{R}_+^N}{\operatorname{argmin}} \frac{1}{2} \| G \operatorname{diag}(p) G^H - C \|_F^2 \qquad (5)$$

with $\| \cdot \|_F$ the Frobenius norm, This form of the CMF problem neglects the effect of the noise. However, this simpler problem streamlines the exposition of the proposed method. Modifications necessary to account for microphone self-noise will be given in Section 3.3.

The CMF problem has several interesting features for source imaging. Firstly, it boils down to a NNLS problem, a well known problem in optimization [5]. Secondly, solutions to the CMF problem are sparse. Indeed, solutions to NNLS problems in the underdetermined setting are sparse, with at least one solution of (5) that has at most $M(M+1)/2$ nonzero coefficients [25]. Finally, the non-negativity constraints allow super-resolution, the ability to separates close sources [26,27].

The objective function of the CMF problem can be rewritten [4]

$$J(p) = \frac{1}{2} \|G \operatorname{diag}(p) G^H - C\|_F^2 \qquad (6)$$
$$ = \frac{1}{2} \mathbf{p}^T \mathbf{A} \mathbf{p} - \mathbf{b}_0^T \mathbf{p} + C \quad (7) $$

where the positive semidefinite matrix A has coefficients

$$ A_{ij} = |\mathbf{g}_i^H \mathbf{g}_j|^2. \quad (8) $$

or

$$ \mathbf{A} = (\mathbf{D}^H \mathbf{D}) \odot (\overline{\mathbf{D}^H \mathbf{D}}), \quad (9) $$

where $\odot$ denotes the Hadamard (component-wise) product and $\overline{\cdot}$ component-wise complex conjugation, and $\mathbf{b}_0$ is defined by its coordinates

$$ \mathbf{b}_{0,i} = \mathbf{g}_i^H \mathbf{C} \mathbf{g}_i, \quad (10) $$

values of a beamforming map computed with unnormalized steering vectors, justifying its notation. C is a constant. The negative gradient $\mathbf{b}(\mathbf{p})$ of J is

$$ \mathbf{b}(\mathbf{p}) = \mathbf{b}_0 - \mathbf{A} \mathbf{p}. \quad (11) $$

## 2.2. Algorithms

Several algorithms can be used to solve the CMF problem. Yardibi et al. [3] used the SeDuMi Matlab package, based on the interior-point algorithm. The LH algorithm is used in Acouil [28], and an optimized version exploiting the structure of the problem was introduced in [4]. It is remarkable that the DAMAS algorithm, introduced by Brooks and Humphreys in 2006 [2], two years before Yardibi et al. introduced the CMF method, converges to a solution to the CMF problem, albeit slowly. Indeed, the DAMAS algorithm is a coordinate descent algorithm for the CMF problem, iterating cyclically through the coefficients. The iterates $\mathbf{p}^k$ are obtained by choosing an index $i_k$, setting $p_j^k = p_j^{k-1}$ for $j \neq i_k$, and

$$ p_{i_k}^k = \max \left( 0, \frac{1}{A_{i_k i_k}} \left( b_{0,i_k} - \sum_{l=0}^{i_k-1} A_{i_k l} p_l^{k-1} - \sum_{l=i_k+1}^{N} A_{i_k l} p_l^{k-1} \right) \right). \quad (12) $$

DAMAS was introduced and is usually used with a deterministic choice of the updates. Using $i_k = k \mod N$ will be referred to as cyclic DAMAS, and iterating from 1 to N, and back from N to 1, and repeating, will be referred to as roundtrip DAMAS. A final variant, Random DAMAS, uses a uniform random choice for $i_k$ at each iteration.

Although the cost of an iteration is low, in $O(N)$, the convergence of DAMAS is known to be slow, in part because every point of the grid if visited, even if the actual solution to the CMF problem is sparse. Moreover, the DAMAS algorithm involves all $N \times N$ coefficients of A explicitly, which necessitates either to store them in memory, or to compute them at each iteration, saving memory but increasing the cost of an iteration from $O(N)$ to $O(MN)$.

## 3. Greedy coordinate descent

The convergence rate of coordinate descent, in particular DAMAS, can be improved by choosing the coordinate to update according to appropriate rules [18]. In particular, we focus on the Maximum Improvement rule that selects the coordinate to update by computing the improvement offered by updating the ith coordinate, and choosing the best improvement.

At a given iteration k, varying the ith coordinate yields

$$ J(\mathbf{p}^{k-1} + \alpha \mathbf{e}_i) = J(\mathbf{p}^{k-1}) + \frac{1}{2} A_{ii} \alpha^2 - b_i(\mathbf{p}^{k-1})\alpha \quad (13) $$

with $\mathbf{e}_i$ the N dimensional vector having value 1 at index i and 0 elsewhere. The minimum of the unconstrained problem is at $\alpha_i^* = \frac{b_i(\mathbf{p}^{k-1})}{A_{ii}}$, yielding an improvement of

$$ \delta_i^* = J(\mathbf{p}^{k-1}) - J(\mathbf{p}^{k-1} + \alpha_i^* \mathbf{e}_i) \quad (14) $$

$$ = \frac{1}{2} \frac{b_i(\mathbf{p}^{k-1})^2}{A_{ii}}. \quad (15) $$

In the case where $p_i^{k-1} + \alpha_i^* < 0$, then the constrained minimum of J is at $\alpha_i^* = -p_i^{k-1}$ (that is, the power of the source i is set at 0), with improvement

$$ \delta_i^* = J(\mathbf{p}^{k-1}) - J(\mathbf{p}^{k-1} + \alpha_i^* \mathbf{e}_i) \quad (16) $$

$$ = -p_i^{k-1} b_i(\mathbf{p}^{k-1}) - \frac{1}{2} A_{ii} (\mathbf{p}_i^{k-1})^2. \quad (17) $$

The improvement obtained by updating the coordinate i is then

$$ \delta_i = \begin{cases} \delta_i^* & \text{if } p_i^{k-1} + \alpha_i^* \ge 0 \\ \delta_i^* & \text{else} \end{cases} \quad (18) $$
Then the i*-th coordinate maximizing the improvement is updated, with $p_{i*}^k = p_{i*}^{k-1} + \frac{b_{i*}(p^{k-1})}{A_{i*}A_{i*}}$ in the first case, $p_{i*}^k = 0$ in the second case, or equivalently,

$$p_{i*}^k = \max(0, p_{i*}^{k-1} + \alpha_{i*}^k), \qquad (19)$$

leaving the other coordinates constant.

The iteration ends with the update of the gradient, with

$$\mathbf{b}(\mathbf{p}^k) = \mathbf{b}(\mathbf{p}^{k-1}) - (\mathbf{p}_{i*}^k - \mathbf{p}_{i*}^{k-1})\mathbf{a}_{i*} \qquad (20)$$

where $\mathbf{a}_{i*}$ is the i*-th column of A.

The algorithm can be stopped after a given number K of iterations, or when a stopping condition is reached. Here, we will stop the iterations when the improvement of the objective function by an iteration falls below a fixed threshold $\epsilon$, i.e. when

$$\frac{J(\mathbf{p}^{k-1}) - J(\mathbf{p}^k)}{J(\mathbf{p}^k)} < \epsilon. \qquad (21)$$

The influence of the choice of $\epsilon$ will be investigated in the numerical experiments, in particular its effect on the resolving power of the method, and overfitting at low SNR.

The algorithm can be simplified and accelerated by normalizing the columns of G, which implies $A_{ii} = 1$ for all i, and rescaling the estimated vector p after stopping the iterations. This version is summarized in Alg. 1, where $\tilde{\epsilon}$ indicates values computed using normalized columns $\tilde{\mathbf{g}}_i$.

**Algorithm 1 DAMAS-MI algorithm**

Initialize the estimated powers: $\tilde{\mathbf{p}}^0 = 0$, Initialize the set of computed columns of A: $A = \emptyset$
Compute the norms $n_i$ of the columns of G
Compute the normalized columns $\tilde{\mathbf{g}}_i = \mathbf{g}_i/n_i$ of the matrix $\tilde{\mathbf{G}}$
Initialize the negative gradient $\tilde{\mathbf{b}}_0 = \tilde{\mathbf{b}}(\mathbf{p}^0)$ with
$$\tilde{\mathbf{b}}_{0,i} = \tilde{\mathbf{g}}_i^H \mathbf{C} \tilde{\mathbf{g}}_i$$

for k=1, ..., K do

Compute

$$\alpha_i = \tilde{b}(\mathbf{p}^{k-1})_i \qquad (22)$$

$$\delta_i^* = \frac{1}{2} \tilde{b}(\mathbf{p}^{k-1})_i^2 \qquad (23)$$

$$\delta_i^\dagger = -\tilde{\mathbf{p}}_i^{k-1} \tilde{b}(\tilde{\mathbf{p}}_i^{k-1})_i - \frac{1}{2}(\tilde{\mathbf{p}}_i^{k-1})^2 \qquad (24)$$

$$\delta_i = \begin{cases} \delta_i^* & \text{if } \tilde{\mathbf{p}}_i^{k-1} + \alpha_i \ge 0 \\ \delta_i^\dagger & \text{else} \end{cases} \qquad (25)$$

Find the index $i^*$ of the maximal $\delta_i$

Update

$$\tilde{\mathbf{p}}_{i*}^k = \max(0, \tilde{\mathbf{p}}_{i*}^{k-1} + \alpha_{i*})$$

Retrieve $\tilde{\mathbf{a}}_{i*}$ from A, or compute it and store it in A

Update

$$\tilde{\mathbf{b}}(\mathbf{p}^k) = \tilde{\mathbf{b}}(\mathbf{p}^{k-1}) - (\tilde{\mathbf{p}}_{i*}^k - \tilde{\mathbf{p}}_{i*}^{k-1})\tilde{\mathbf{a}}_{i*}$$

end for

Rescale the estimated powers with $p_i = \tilde{\mathbf{p}}_i^K/n_i^2$

## 3.1. Computational complexity

The computational complexities of classical DAMAS and the proposed DAMAS-MI are now compared. The computational complexity of the algorithms is composed of the complexity of the initialization, and that of the iterations.

An initialization task common to all algorithms is the computation of the beamforming map $\mathbf{b}_0$, with cost $O(NM^2)$, that will not be considered when comparing the methods. In DAMAS, the matrix A can be precomputed with a cost $O(N^2M)$. This case will be denoted DAMAS offline.

The cost of an iteration of DAMAS, i.e. computation of Eq. (12), is either $O(N)$ when A is precomputed, or $O(NM)$ when the necessary coefficients of A are not stored, but computed when needed (DAMAS online). The total complexity of the DAMAS algorithm is thus $O(N^2M) + O(KN)$ for the offline case, and $O(KNM)$ for the online case. In typical applications, K is much larger than N, making the online DAMAS much less efficient in terms of time complexity. However, online DAMAS is frugal in memory compared to offline DAMAS, avoiding the storage of the large $N \times N$ matrix A, limiting the necessary storage to the iterates $\mathbf{p}^k$ of size N.
<table><caption>**Table 1**
Computational complexity of DAMAS variants.</caption><thead><tr><th>Algorithm</th><th>Space</th><th>Time</th></tr></thead><tbody><tr><td>DAMAS-MI</td><td>O(S<sub>K</sub>N)</td><td>O(S<sub>K</sub>NM) + O(KN)</td></tr><tr><td>DAMAS offline</td><td>O(N<sup>2</sup>)</td><td>O(N<sup>2</sup>M) + O(KN)</td></tr><tr><td>DAMAS online</td><td>O(N)</td><td>O(KNM)</td></tr><tr><td>CLEAN-SC</td><td>O(N)</td><td>O(KNM<sup>2</sup>)</td></tr><tr><td>CLEAN-SC with Eq. (20)</td><td>O(N)</td><td>O(KNM)</td></tr></tbody></table>

In the proposed DAMAS-MI algorithm, the cost of an iteration consists of computing the improvement for each index *i* with Eqs. (15), (17) and (18), for a cost linear in *N*, and updating the gradient. The update of the gradient with Eq. (20) costs either $O(N)$ if $a_{i*}$ is stored in memory, or $O(NM)$ if $a_{i*}$ is computed from Eq. (8). Due to the sparsity of the solution, it is likely that only a few coefficients will be chosen in the iterations. This can be exploited by computing $a_{i*}$ only when needed with a cost $O(NM)$, and storing it for further usage. The number of columns computed and stored in $\Lambda$ up to iteration *K* (or equivalently, the number of unique coefficients updated until iteration *K*) is denoted $S_K$, with the obvious bounds $1 \le S_K \le \min(N, K)$. The total complexity is then $O(S_K NM) + O(KN)$, for a storage space of $O(S_K N)$ for the computed columns of $\Lambda$.

As a conclusion, for a given number of iterations, the complexity in time of the proposed DAMAS-MI algorithm is lower than the complexities of DAMAS (see Table 1).

Note that this discussion does not take into account the convergence rate of the objective function in function of the number of iterations. Numerical results will show that the proposed DAMAS-MI is both more efficient for a given number of iterations, and converges faster than the DAMAS algorithm.

## 3.2. Connections with other methods

It was already shown in [4] that the original DAMAS algorithm is closely connected to the CMF method, as its iterates converge towards a solution to the CMF problem. We give here additional connections, in particular to CLEAN-PSF and CLEAN-SC.

**Beamforming.** The connection with beamforming is obvious when considering the first iteration of the algorithm: a new source is identified by maximizing

$$ \delta_i^* = \frac{b_{0,i}^2}{2A_{ii}} \qquad (26) $$

$$ = \frac{1}{2} \left( \frac{g_i^H C g_i}{\|g_i\|_2^2} \right)^2 \qquad (27) $$

which corresponds to a beamforming map obtained by the formulation IV of steering vectors [29]. The power is then updated with

$$ p_{i*}^1 = \frac{b_{0,i*}}{A_{i*i*}} \qquad (28) $$

$$ = \frac{g_i^H C g_i}{\|g_i\|_2^4}. \qquad (29) $$

with here formulation III of the steering vectors. The combination of formulations IV and III has been shown to be optimal in terms of bias and variance for estimating the position and power of a source [30]. The next iterations proceeds similarly, operating on the beamforming map $b(p^k)$ obtained by removing the contribution of the identified sources in the SCM. We note that after the first iteration, $b(p^k)$ is not necessarily positive, and the location of the maximum of the absolute value of the normalized beamforming map is selected.

**CLEAN-PSF.** CLEAN-PSF [6] operates in a similar way, by selecting a source at the maximum of a beamforming map and removing its contribution with an optional damping parameter $\phi$. However, this algorithm does not have a procedure for diminishing, or removing entirely, powers of sources that have been overestimated at an early iteration. We note that this interpretation of the CLEAN-PSF algorithm suggests to use formulation III to locate the source, and formulation IV to estimate its power, as in beamforming. Additionally, Eq. (20) can be used to update the beamforming map instead of computing it from the cleaned SCM at each iteration, reducing the cost of an iteration from $O(NM^2)$ to $O(NM)$. These remarks remain valid for CLEAN-SC, and will be followed in the numerical experiments.

**Matching pursuit.** The greedy selection of the updated coefficient bears resemblance with the Matching Pursuit algorithm. However here the nonnegativity constraint is considered. We note that Orthogonal Matching Pursuit was proposed for source deconvolution in [20,31], with no nonnegativity constraint, and with the goal of decomposing the beamforming map, and not the covariance matrix.

**Lawson–hanson.** The Lawson–Hanson algorithm [5] finds a minimizer of the CMF problem in a finite number of iterations. Similarly to the other methods, sources are identified one at a time, and powers of all sources are updated. This update may need solving several unconstrained least-squares problems, which comes at high computational cost.
<figure><img src="image_6.png" /><figcaption>Fig. 1. Experimental setup. (a) Acoustical sources. (b) Positions of the array microphones.</figcaption></figure>

### 3.3. Diagonal removal

In experimental settings where self-noise of the microphones is non-negligible, it is usual to use diagonal removal in beamforming, CLEAN-SC, DAMAS or similar methods. The influence of uncorrelated microphone noise in the covariance matrix is the addition of diagonal terms. For better robustness with respect to this noise, the diagonal terms of the covariance matrix can be neglected, with a modified Frobenius norm $||·||_{F,dr}$ where the diagonal is not considered. The obtained criterion is

$$J_{\text{dr}}(\mathbf{p}) = \frac{1}{2} \|\text{Gdiag}(\mathbf{p})\mathbf{G}^H - \mathbf{C}\|_{F,dr}^2 \quad (30)$$

$$= \frac{1}{2} \|\text{Gdiag}(\mathbf{p})\mathbf{G}^H - \mathbf{C}\|_{F}^{2} - \frac{1}{2} \sum_{m=1}^{M} \left( (\text{Gdiag}(\mathbf{p})\mathbf{G}^H)_{mm} - C_{mm} \right)^2 \quad (31)$$

where the first term is $J(\mathbf{p})$, and the second term, denoted $J_d(\mathbf{p})$, cancels the diagonal coefficients in the Frobenius norm. Remarking that $(\text{Gdiag}(\mathbf{p})\mathbf{G}^H)_{mm} = \sum_{i=1}^{N} p_i |g_{i,m}|^2$ and expanding the square, the second term is

$$J_d(\mathbf{p}) = \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} p_i p_j \sum_{m=1}^{M} |g_{i,m} g_{j,m}|^2 - \sum_{i=1}^{N} p_i \sum_{m=1}^{M} |g_{i,m}|^2 C_{mm} + \frac{1}{2} \sum_{m=1}^{M} C_{mm}^2 \quad (32)$$

With the matrix $A_d$ and the vector $b_{0,d}$ with coefficients

$$A_{d,i,j} = \sum_{m=1}^{M} |g_{i,m} g_{j,m}|^2 \quad (33)$$

$$b_{0,d,i} = \sum_{m=1}^{M} |g_{i,m}|^2 C_{mm} \quad (34)$$

$$= g_i^H \text{diag}(C) g_i, \quad (35)$$

where diag(C) is the diagonal matrix with the same diagonal coefficients as C, $J_d(\mathbf{p})$ can be rewritten

$$J_d(\mathbf{p}) = \frac{1}{2} \mathbf{p}^T A_d \mathbf{p} - \mathbf{b}_{0,\mathrm{d}}^T \mathbf{p} + C_d \quad (36)$$

with $C_d$ a constant. Combining with Eq. (7), the criterion $J_{\text{dr}}(\mathbf{p})$ is written

$$J_{\text{dr}} = \frac{1}{2} \mathbf{p}^T A_{\text{dr}} \mathbf{p} - \mathbf{b}_{0,\text{dr}}^T \mathbf{p} + C_{\text{dr}} \quad (37)$$

where the matrix $A_{\text{dr}}$ is given by its coefficients

$$A_{\text{dr},ij} = A_{ij} - A_{d,ij} \quad (38)$$
<figure><img src="image_7.png" /><figcaption>Fig. 2. Comparison of (a) beamforming (formulation IV), (b) LH, (c) method DAMAS-MI, DAMAS in its (d) random, (e) cyclic and (f) roundtrip variants, as well a CLEAN-SC with (g) $\phi = 1$ and (h) $\phi = 0.1$, Powers in dB.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 3. Decay of the objective in function of the iteration count, 2D case.</figcaption></figure>

$$ = |g_i^H g_j|^2 - \sum_{m=1}^{M} |g_{i,m}^2 g_{j,m}|, \qquad (39) $$

the vector $b_{0,dr}$ by

$$ b_{0,dr,i} = b_{0,i} - b_{0,d,i} \qquad (40) $$

$$ = g_i^H (C - \text{diag}(C)) g_i \qquad (41) $$

and $C_{dr} = C - C_d$ is a constant.

As a conclusion, using a diagonal removal formulation consists in using $A_{dr}$ and $b_{0,dr}$ instead of the original $A$. To ensure that $A_{dr,ii} = 1$, the columns of $G$ are normalized by $n_i = \sqrt{\|g_i\|^4 - \sum_{m=1}^M |g_{i,m}^4|}$ in Alg. 1.

## 4. Experimental and numerical results

The proposed method is tested using simulated data and experimental data recorded in an anechoic chamber, and compared to LH and CLEAN-SC in particular. The setup of the experiments is identical to the experiments described in [32].
<figure><img src="image_9.png" /><figcaption>Fig. 4. Decay of the objective in function of computational time, 2D case.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5. Size of the support of the iterates, 2D case.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 6. Number of necessary columns of A in function of the iteration count, 2D case.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 7.** Comparison of estimated sources distribution by (a), (b) beamforming, (c), (d) LH, (e), (f) DAMAS-MI and (g), (h) CLEAN-SC at (a), (c), (e), (g) *F* = 781 Hz and (b), (d), (f), (h) *F* = 5004 Hz. Powers in dB.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Influence of the number of sources. Positions of the sources used in the simulation (a), and beamforming map (b) when all sources are active.</figcaption></figure>

The soundfield is measured by a microphone array with 128 digital MEMS microphone (Invensense - INMP441) with a 26 dBFS sensitivity (1kHz, 94dB SPL) and a flat response in the band [150 Hz – 15 kHz]. The geometry of the array is given in Fig. 1. The antenna counts 128 elements distributed along 16 linear rays. On each ray, the 8 microphones are spaced regularly at 17 cm intervals, the distance between the first microphone and the center of the array follows a pseudo-random distribution.

In the experiment, four sources, baffled broadband omnidirectional loudspeakers (Visaton-BF32 - [150Hz–20kHz]) are used (see Fig. 1).

The microphone signals are sampled at *F*<sub>s</sub> = 50 kHz and analyzed by Short-Term Fourier Transform, with a 2048 samples Hann window (41 ms duration and 75% overlap), yielding 939 snapshots.

Numerical computation are conducted on a laptop with 32 GB of memory and an Intel Core Ultra 9 185H CPU using python code available online [33]. The indicated sources powers are given at a distance of 1 meter to the source.

### 4.1. 2D case, experimental data

The method is first tested on a 2D case, where sources are searched in a planar region at a distance of 4.4 m to the array, of rectangular shape with dimensions 3 m × 1 m, discretized with a step *Δ* of 0.02 m, yielding *N* = 7701 points.

Fig. 2 shows the results, at *F* = 2002 Hz, of (a) beamforming (formulation IV), (b) LH, (c) the proposed method DAMAS-MI, DAMAS in its (d) random, (e) cyclic and (f) roundtrip variants, as well a CLEAN-SC with (g) *φ* = 1 and (h) *φ* = 0.1, with 75000
<figure><img src="image_14.png" /><figcaption>Fig. 9. Influence of the number of sources. (a) Total specific error, (b) Inverse error, (c) Computational time.</figcaption></figure>

iterations for each algorithm (CLEAN-SC is stopped when the maximal value of the beamforming map falls below $10^{-3}$ times its initial value).

The decrease of the objective function

$$ \delta_k = \frac{J(\mathbf{p}^k) - J(\mathbf{p}^*)}{J(\mathbf{p}^*)} \qquad (42) $$

relative to the actual minimum with the exact solution $\mathbf{p}^*$ obtained by the Lawson-Hanson algorithm, is plotted on Fig. 3 in function of the number of iterations, showing that the proposed method is not only more efficient than DAMAS in terms of computational cost per iteration, but also exhibits faster convergence in function of the number of iterations. CLEAN-SC is stopped early, as it is not capable of reducing the power of previously identified sources, with, as expected, slower decay for $\phi = 0.1$.

Fig. 4 shows the decrease of the objective in function of time. While LH converges in finite time, early iterations of DAMAS-MI have lower objectives than LH.
<figure><img src="image_15.png" /><figcaption>Fig. 10. Errors with two sources for varying frequency. (a) Specific error source 1, (b) Specific error source 2, (c) Inverse error.</figcaption></figure>

Fig. 5 shows the evolution of the size of the support of the iterates, with sparse iterates for DAMAS-MI. The number of columns of A to be explicitly computed is given in Fig. 6. As expected, DAMAS-MI needs a reduced number of columns of A compared to the other versions of DAMAS, reducing its memory usage by a factor 30.

Further reconstructions are given on Fig. 7 at two additional frequencies, $F = 781$ Hz and $F = 5004$ Hz, for beamforming, LH, DAMAS-MI and CLEAN-SC with $\phi = 0.1$. Here DAMAS-MI is stopped using $\epsilon = 10^{-5}$. At low frequency, DAMAS-MI, with the value of $\epsilon$ used here, cannot resolve the two central sources. CLEAN-SC cannot resolve the central sources either, and spread the outer sources on curves. At higher frequency, all methods performs well.

## 4.2. Simulations

Simulations are used to assess the resolution, accuracy, and numerical efficiency of the proposed method compared to LH, CLEAN-SC, and CMF-SFW, a gridless method based on the Sliding Frank-Wolfe algorithm which was shown to outperform gridless methods based on CLEAN-SC, HR-CLEAN-SC [34] or global optimization.

The methods are compared using the inverse level error and specific level errors defined in [21]. The specific level error for a given source is defined here as the mean squared error between its power and the sum of the powers estimated in a ball of radius R around its position. The inverse level error is the sum of the powers estimated outside of the union of these balls.
<figure><img src="image_16.png" /><figcaption>Fig. 11. Errors with varying SNR. (a) Total specific error, (b) Inverse error (zero errors for COMET2 and CLEAN-SC with $\phi = 1$ are not plotted).</figcaption></figure>

<table><caption>**Table 2**
Powers of the simulated sources.</caption><tr><td>Source</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td></tr><tr><td>Powers (Pa²)</td><td>1</td><td>0.8</td><td>0.5</td><td>1</td><td>0.7</td><td>0.8</td><td>0.4</td><td>0.3</td></tr></table>

*Number of sources.* In a first set of simulations, the performances are evaluated in function of the number of sources present in the domain of interest, with a maximum of 8 sources. The positions and ranks of the sources are given in Fig. 8, at a distance of 3 meters from the array, with $F = 1642$ Hz. Their powers are given in Table 2. We used here $R = 0.05$ m. The power of the noise is $\sigma^2 = 0.1$ Pa², resulting in an SNR ranging from 0.14 dB (one source) to 7.31 dB (eight sources).

The sum of the specific errors is plotted in Fig. 9-(a) in function of the number of sources. This sum is naturally increasing, as more sources are taken into account. Sharp increases are seen when the third and fifth sources are introduced. Being close to the first and second sources respectively, their introduction implies resolution issues. The inverse level errors are plotted on Fig. 9-(b).

Best performances are obtained by CMF-SFW and LH, with the caveat that the number of sources is assumed to be known and used as the number of iterations for CMF-SFW. Moreover, computation times for LH and CMF-SFW are large (see Fig. 9-(c)), and increasing with the number of sources for CMF-SFW.

Performances of CLEAN-SC are limited by its resolution power, with errors increasing significantly when the third and fifth sources are introduced. DAMAS-MI offers a good compromise between performances and computational time, that can be tuned with the stopping criterion $\epsilon$.

*Resolution.* The resolution of the tested methods is assessed more precisely by simulations with two sources at fixed positions, at varying frequency between 273 Hz and 3248 Hz. The specific errors and the inverse level error are plotted on Fig. 10 for the methods tested above, with $R = 0.02$ m, and a SNR of 2.70 dB. It can be seen that the resolving power of DAMAS-MI can be improved by
<figure><img src="image_17.png" /><figcaption>**Fig. 12.** Errors with pairs of correlated sources. (a) Total specific error, (b) Inverse error (zero errors for COMET2 and CLEAN-SC with $\phi = 1$ are not plotted).</figcaption></figure>

lowering $\epsilon$ (raising the number of iterations), at a cost of longer computations. CLEAN-SC is capable of resolving the sources only at higher frequencies, with slightly better performances obtained at $\phi = 0.1$.

These numerical experiments show that DAMAS-MI has better performances than CLEAN-SC, in particular in resolution power, while maintaining reasonable computation times.

**SNR.** The performances are now evaluated in function of the SNR, ranging from -15 dB to 15 dB, plotted on Fig. 11. To avoid performances issues caused by closely spaced sources, only sources 1, 2, 4 and 6 on Fig. 8 are used. At low SNR, best performances are obtained by DAMAS-MI with $\epsilon = 10^{-6}$. These performances can be explained by the value of $\epsilon$ avoiding overfitting of the model to the data. This observation is supported by the higher errors obtained with $\epsilon = 10^{-8}$, and LH, which returns the exact solution of the optimization problem.

**Source correlation.** In the theoretical model, sources are assumed to be uncorrelated. The robustness of the method to source correlation is assessed, using two pairs of correlated sources (sources 1 and 2, and sources 4 and 6), with correlations ranging from 0 (no correlation) to 1 (full correlation, the sources emit identical signals). As shown on Fig. 12, errors increase with correlation for all methods. However, the CMF based methods (LH and DAMAS-MI) are less affected by source correlation. For CLEAN-SC, the degradation of the performances with correlation can be attributed to the manipulation of the steering vector at each iteration aiming at removing sources correlated to the identified source.

**Noise correlation.** Uncorrelated noise (i.e., diagonal covariance matrix of the noise) is assumed in the model. We test here the robustness of the method to noise correlation, using a covariance function $c(x, y) = \sigma^2 \exp(-\|x - y\|_2/l)$ where $l$ is a correlation length, ranging from 0 to 0.4 m. The errors plotted on Fig. 13 show that the influence of noise correlation on the performances is rather weak for the proposed method.
<figure><img src="image_18.png" /><figcaption>Fig. 13. Errors with correlated noise. (a) Total specific error, (b) Inverse error (zero errors for COMET2 and CLEAN-SC with $\phi = 1$ are not plotted).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14. Three dimensional distribution of sources obtained at $F = 2002$ Hz by DAMAS-MI ((a) front view and (b) top view) and LH ((c) front view and (d) top view).</figcaption></figure>

### 4.3. 3D case, experimental data

The method is now tested on a 3D case, with a region of interest of dimensions 3.4 × 1.6 × 1.6 (in meters), discretized with a step $\Delta$ of 0.02 m, containing 1,121,931 points. Fig. 14 compares the source distributions at $F = 2002$ Hz obtained by LH (computational time 17 min) and DAMAS-MI with $\epsilon = 10^{-3}$ (23 iterations, 4.7 s, including the initialization of the gradient by beamforming Eq. (10),
<figure><img src="image_20.png" /><figcaption>Fig. 15. Decay of the objective in function of computational time, 3D case.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 16. Support of the solution and computed columns of A, 3D case.</figcaption></figure>

2.4 s). Convergence of the objective function is plotted on Fig. 15 for up to 200,000 iterations, showing lower values of the objective function at early iterations for DAMAS-MI compared to LH.

The size of the support of the iterates and the number of computed columns of A are plotted in Fig. 16. Only 719 columns of A are needed, 0.06% of the total number of columns.

The effect of the size of the grid is demonstrated on DAMAS-MI with grid steps $\Delta$ of 0.02 m (1,121,931 points), 0.05 m (75,141) and 0.1 m (10,115). Fig. 17 shows the front view of the source distributions estimated at $F = 6982$ Hz. At coarser grid steps, the sources are spread over several grid points (all sources at $\Delta = 0.1$ m, left source at $\Delta = 0.05$ m), while there are mostly limited to one grid points at $\Delta = 0.02$ m.

This property is shown at frequencies of $F = 6982$ Hz, 7983 Hz, 8984 Hz and 9985 Hz. The estimated power for each grid node is plotted in Fig. 18 in decreasing order, for the 10 largest values. Compared to the coarser grids, the source powers obtained with the finer grid with $\Delta = 0.01$ m are concentrated on the first four most powerful sources, the following sources (due to spreading of the sources, noise, model errors, etc.) having lower powers than with coarser grid steps.

## 5. Conclusion

A modification of the DAMAS algorithm is proposed, based on a greedy selection of the coordinate to be updated at each iteration. The choice of the grid point, and the amount of power added or removed, is consistent with classical beamforming, and are based on the interpretation of DAMAS as a coordinate descent algorithm for the CMF problem.

Simulations and application to experimental measurements showed that convergence is indeed faster than deterministic or random DAMAS, with convergence 40 times faster in the 2D case, and lower memory usage (memory usage divided by 30). Additionally the iterates are sparse, an expected property of at least one of the solution of the CMF problem. The proposed algorithm is complementary to the Lawson-Hanson algorithm and gridless methods, as it offers an approximate, but fast, solution to the CMF problem, with an algorithm similar to established methods from the state of the art and based on elementary mathematical operations (in particular, no matrix inversion is needed).

Comparison to the similar CLEAN-SC algorithm showed that DAMAS-MI has better performances in terms of power estimation, that can be attributed to its better resolving power and better robustness with respect to correlations between sources.
<figure><img src="image_22.png" /><figcaption>Fig. 17. Front view of the source distributions estimated by DAMAS-MI at $F = 6982$ Hz, with (a) $\Delta = 0.1$ m, (b) $\Delta = 0.05$ m and (c) $\Delta = 0.02$ m.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 18. Powers of the 20 most powerful grid nodes, in decreasing orders, for $\Delta = 0.1$ m, $\Delta = 0.05$ m and $\Delta = 0.02$ m, for DAMAS-MI with $\varepsilon = 10^{-3}$ at (a) $F = 6982$ Hz, (b) $F = 7983$ Hz., (c) $F = 8984$ Hz and (d) $F = 9985$ Hz.</figcaption></figure>

The numerical efficiency of the method makes large scale applications possible, e.g. on a three dimensional grid with more than one million points, with running time in the order of the time necessary to compute a beamforming map. This make the use of fine grids possible, enhancing the sparsity of the estimated source distribution, as demonstrated in a case with four point sources.

CRediT authorship contribution statement

Gilles Chardon: Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Conceptualization.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data is available online, Ref. 33.

References

[1] R. Merino-Martínez, P. Sijtsma, M. Snellen, T. Ahlefeldt, J. Antoni, C.J. Bahr, D. Blacodon, D. Ernst, A. Finez, S. Funke, T.F. Geyer, S. Baxter, G. Herold, X. Huang, W.M. Humphreys, Q. Leclère, A. Malgoezar, U. Michel, T. Padois, A. Pereira, C. Picard, E. Sarradj, H. Siller, D.G. Simons, C. Spehr, A review of acoustic imaging methods using phased microphone arrays: Part of the aircraft noise generation and assessment special issue, CEAS Aeronaut. J. 10 (1) (2019) 197–230, http://dx.doi.org/10.1007/s13272-019-00383-4.

[2] T.F. Brooks, W.M. Humphreys, A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays, J. Sound Vib. 294 (4) (2006) 856–879, http://dx.doi.org/10.1016/j.jsv.2005.12.046.

[3] T. Yardibi, J. Li, P. Stoica, L.N. Cattafesta, Sparsity constrained deconvolution approaches for acoustic source mapping, J. Acoust. Soc. Am. 123 (5) (2008) 2631–2642, http://dx.doi.org/10.1121/1.2896754.

[4] G. Chardon, J. Picheral, F. Ollivier, Theoretical analysis of the DAMAS algorithm and efficient implementation of the covariance matrix fitting method for large-scale problems, J. Sound Vib. 508 (2021) 116208, http://dx.doi.org/10.1016/j.jsv.2021.116208.

[5] C.L. Lawson, R.J. Hanson, Solving least squares problems, Soc. Ind. Appl. Math. (1995) http://dx.doi.org/10.1137/1.9781611971217.

[6] P. Sijtsma, CLEAN Based on Spatial Source Coherence, Int. J. Aeroacoustics 6 (4) (2007) 357–374, http://dx.doi.org/10.1260/147547207783359459.

[7] G. Chardon, Gridless covariance matrix fitting methods for three dimensional acoustical source localization, J. Sound Vib. 551 (2023) 117608, http://dx.doi.org/10.1016/j.jsv.2023.117608.

[8] A. Goudarzi, Global, and local optimization beamforming for broadband sources, J. Acoust. Soc. Am. 155 (1) (2024) 262–273, http://dx.doi.org/10.1121/10.0024247.

[9] B. von den Hoff, R. Merino-Martínez, D.G. Simons, M. Snellen, Using global optimization methods for three-dimensional localization and quantification of incoherent acoustic sources, JASA Express Lett. 2 (5) (2022) 054802, http://dx.doi.org/10.1121/10.0010456.

[10] Y. Zhao, Y. He, H. Chen, Z. Zhang, Z. Xu, Three-dimensional grid-free sound source localization method based on deep learning, Appl. Acoust. 227 (2025) 110261, http://dx.doi.org/10.1016/j.apacoust.2024.110261.

[11] C. Kayser, A. Kujawski, E. Sarradj, A fast data-driven method for inverse microphone array signal processing, JASA Express Lett. 3 (4) (2023) 042401, http://dx.doi.org/10.1121/10.0017882.

[12] H. Liang, G. Zhou, X. Tu, A. Jakobsson, X. Ding, Y. Huang, Learning an interpretable end-to-end network for real-time acoustic beamforming, J. Sound Vib. 591 (2024) 118620, http://dx.doi.org/10.1016/j.jsv.2024.118620.

[13] A. Kujawski, E. Sarradj, Fast grid-free strength mapping of multiple sound sources from microphone array data using a Transformer architecture, J. Acoust. Soc. Am. 152 (5) (2022) 2543–2556, http://dx.doi.org/10.1121/10.0015005.

[14] Y. Yang, Y. Yang, Z. Chu, Off-grid deconvolution beamforming for acoustic source identification, Appl. Acoust. 218 (2024) 109909, http://dx.doi.org/10.1016/j.apacoust.2024.109909.

[15] T. Lobato, R. Sottek, M. Vorländer, Deconvolution with neural grid compression: A method to accurately and quickly process beamforming results, J. Acoust. Soc. Am. 153 (4) (2023) 2073, http://dx.doi.org/10.1121/10.0017792.

[16] E.J.G. Arcondoulis, D. Ragni, D. Fiscaletti, R. Merino-Martinez, Y. Liu, Acoustic response of structured and randomized porous blunt trailing edges subject to turbulent boundary layers, J. Acoust. Soc. Am. 156 (2) (2024) 1029–1040, http://dx.doi.org/10.1121/10.0028189.

[17] M. Awasthi, D. Moreau, P. Croaker, P. Dylejko, The sound radiated by tip clearances submerged in a boundary layer, Appl. Acoust. 238 (2025) 110741, http://dx.doi.org/10.1016/j.apacoust.2025.110741.

[18] J. Nutini, M. Schmidt, I.H. Laradji, M. Friedlander, H. Koepe, Coordinate descent converges faster with the Gauss-Southwell rule than random selection, in: International Conference on Machine Learning, ICML’15, 2015, pp. 1632–1641.

[19] T. Padois, A. Berry, Two and three-dimensional sound source localization with beamforming and several deconvolution techniques, Acta Acust. United Acust. 103 (2017) 392–480, http://dx.doi.org/10.3813/AAA.919069.

[20] T.F. Bergh, Deconvolution approach to the mapping of acoustic sources with matching pursuit and matrix factorization, J. Sound Vib. 459 (2019) 114842, http://dx.doi.org/10.1016/j.jsv.2019.07.008.

[21] G. Herold, E. Sarradj, Performance analysis of microphone array methods, J. Sound Vib. 401 (2017) 152–168, http://dx.doi.org/10.1016/j.jsv.2017.04.030.

[22] W. Ma, X. Liu, Compression computational grid based on functional beamforming for acoustic source localization, Appl. Acoust. 134 (2018) 75–87, http://dx.doi.org/10.1016/j.apacoust.2018.01.006.

[23] W. Ma, X. Liu, Improving the efficiency of DAMAS for sound source localization via wavelet compression computational grid, J. Sound Vib. 395 (2017) 341–353, http://dx.doi.org/10.1016/j.jsv.2017.02.005.

[24] D. Blacodon, G. Elias, Level estimation of extended acoustic sources using a parametric method, J. Aircr. 41 (6) (2004) 1360–1369, http://dx.doi.org/10.2514/1.3053.

[25] C. Boyer, A. Chambolle, Y.D. Castro, V. Duval, F. de Gournay, P. Weiss, On representer theorems and convex regularization, SIAM J. Optim. 29 (2) (2019) 1260–1281, http://dx.doi.org/10.1137/18M1200750.

[26] V.I. Morgenshtern, E.J. Candès, Super-resolution of positive sources: The discrete setup, SIAM J. Imaging Sci. 9 (1) (2016) 412–444, http://dx.doi.org/10.1137/15M1016552.

[27] M. Ferreira Da Costa, Y. Chi, Compressed super-resolution of positive sources, IEEE Signal Process. Lett. 28 (2021) 56–60, http://dx.doi.org/10.1109/LSP.2020.3045343.

[28] E. Sarradj, G. Herold, A python framework for microphone array data processing, Appl. Acoust. 116 (2017) 50–58, http://dx.doi.org/10.1016/j.apacoust.2016.09.015.

[29] E. Sarradj, Three-dimensional acoustic source mapping with different beamforming steering vector formulations, Adv. Acoust. Vib. 2012 (2012) e292695, http://dx.doi.org/10.1155/2012/292695.

[30] G. Chardon, Theoretical analysis of beamforming steering vector formulations for acoustic source localization, J. Sound Vib. 517 (2022) 116544, http://dx.doi.org/10.1016/j.jsv.2021.116544.

[31] T. Padois, A. Berry, Orthogonal matching pursuit applied to the deconvolution approach for the mapping of acoustic sources inverse problem, J. Acoust. Soc. Am. 138 (6) (2015) 3678–3685, http://dx.doi.org/10.1121/1.4937609.
[32] G. Chardon, F. Ollivier, J. Picheral, Localization of sparse and coherent sources by orthogonal least squares, J. Acoust. Soc. Am. 146 (2019) 4873–4882, http://dx.doi.org/10.1121/1.5138931.

[33] G. Chardon, DAMAS-MI, 2025, http://dx.doi.org/10.5281/zenodo.17063650, https://zenodo.org/record/17063650 (Accessed 05 September 2025).

[34] P. Sijtsma, R. Merino-Martinez, A.M. Malgoezar, M. Snellen, High-resolution CLEAN-SC: Theory and experimental validation, Int. J. Aeroacoustics 16 (4–5) (2017) 274–298, http://dx.doi.org/10.1177/1475472X17713034.