Improved iterative reweighted L1 norm minimization method for sound source identification

2.1. Compressed sensing reconstruction algorithm

The fundamental concept of compressed sensing theory involves three sequential operations: first representing the original signal as a sparse signal compatible with compressed sensing processing, then performing compressive sampling on the sparse signal using a measurement matrix to acquire measurement data and ultimately recovering the sparse signal from these measurements through reconstruction algorithms. In the context of sound source identification applications, the critical implementation challenge resides in constructing an appropriate sensing matrix, where successful acquisition of this matrix determines the effectiveness of spatial sound field reconstruction and direction-of-arrival estimation.

A planar microphone array measurement model is initially established as shown in Fig. 1, comprising $M$ microphones (denoted by ●). With the array center as the origin, a Cartesian coordinate system is constructed where the Direction of Arrival (DOA) of sound sources is characterized by coordinates (${\theta }_i$, ${\phi }_i$). Here, ${\theta }_i$ represents the elevation angle between the incident direction of the $i$-th sound source and the $z$-axis, while ${\phi }_i$ denotes the azimuth angle between the $x$-axis and the projection of the $i$-th source’s incident direction onto the $x$-$y$ plane, with angular constraints defined as 0° $\le {\theta }_i\le$ 90° and 0° $\le {\phi }_i<$ 360°.

Fig. 1Planar array sampling acoustic signal model

Assuming the target sound source region is discretized into $N$ fixed grid points, the $M$-dimensional vector $\mathbf{P}$ formed by the microphone-received signals can be expressed as:

where, $\mathbf{A}$ is the sensing matrix; $\mathbf{q}={\left[\left(q_1,q_2,\dots ,q_N\right)\right]}^T$ represents the source strength distribution vector.

Let $\mathbf{x}={\left[(x_1,x_2,\dots ,x_M)\right]}^T$ and $\mathbf{y}={\left[(y_1,y_2,\dots ,y_M)\right]}^T$ represent the $M$-dimensional vectors composed of the $x$-axis and $y$-axis coordinates of all microphones, respectively. The sensing matrix ($\mathbf{A}$) can be expressed as:

where, $\mathbf{d}\left(t_{1i},t_{2i}\right)={\left[\exp \left(j2\pi \left(x_1{t}_{1i}+y_1{t}_{2i}\right)\right),\exp \left(j2\pi \left(x_2{t}_{1i}+y_2{t}_{2i}\right)\right),\dots ,\exp \left(j2\pi \left(x_M{t}_{1i}+y_M{t}_{2i}\right)\right)\right]}^T$, $t_{1i}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_i/\lambda$; $t_{2i}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_i/\lambda$, $\lambda$ is the wavelength.

In practice, when $M\times N$, Eq (1) becomes an underdetermined system of linear equations, for which no analytical solution exists. However, if the vector $\mathbf{q}$ possesses sparsity, it can be accurately recovered by solving the following $l_0$-norm minimization problem:

where $\xi$ represents the tolerance for the noise signal ($\mathbf{n}$). Although Eq. (3) is an NP-hard problem and generally difficult to solve, it is equivalent to the following $l_1$-norm minimization problem:

For Eq. (4), it can be solved using the convex optimization toolkit CVX to obtain the source strength corresponding to each fixed grid point, thereby enabling sound source DOA (Direction of Arrival) estimation and source strength quantification.

2.2. Improved iteratively reweighted L1 minimization method

The output results obtained by directly solving Eq. (4) often exhibit certain deviations. IRL1 further reduces these deviations by performing iterative optimization of the sound source distribution based on the L1-norm minimization solution. This method first employs a logarithmic sum penalty function that promotes sparsity more effectively than the $l_1$ norm, specifically $\sum _{n=1}^N\mathrm{l}\mathrm{n}({\left|q_n\right|}^2+ϵ)$, to construct the objective function, leading to the following optimization problem:

where, $ϵ$ is a positive parameter that serves a dual role: on one hand, it ensures the logarithmic function is properly defined, and on the other hand, it acts as a control parameter for the iterative process. By initializing $ϵ$ to a small value (e.g., 1) and gradually reducing it to zero during iterations, it guarantees that the global optimal solution of Eq. (5) converges to a neighborhood near the true solution.

Let ${\mathbf{q}}^{\left(\gamma \right)}=\left[q_1^{\left(\gamma \right)},q_2^{\left(\gamma \right)},\dots ,q_N^{\left(\gamma \right)}\right]$ denote the source strength distribution vector obtained after the $\gamma$-th iteration. In the $\gamma +1$-th iteration, construct a surrogate function $\mathrm{\Omega }\left(q\right)$ for $L\left(q\right)$, satisfying $\mathrm{\Omega }\left(q\right)-L\left(q\right)\ge 0$, where equality holds if and only if $\mathbf{q}={\mathbf{q}}^{\left(\gamma \right)}$. Here:

By removing the constant term in Eq. (6), we obtain:

Correspondingly, Eq. (5) is reformulated into the following surrogate function form:

Let the weighted matrix $\mathbf{W}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\left[w_1{,w}_2,\dots ,w_N\right]}^T\right)$, where the $n$-th weight coefficient $w_n^{\left(\gamma \right)}$ is expressed as follows:

Therefore, the problem is transformed into solving the following equation:

Utilize the CVX toolbox to iteratively solve Eq. (10) for sound source DOA estimation.

Compared to the traditional iteratively reweighted $l_1$-norm minimization algorithm, the proposed algorithm introduces a log-sum penalty function with enhanced sparsity-promoting capability, and simplifies the problem into a form that efficiently solves for the source strength distribution vector by constructing a surrogate function.

3.1. Performance analysis of various algorithms at different frequencies

The numerical simulation encompassed two distinct scenarios: a single sound source and dual coherent sound sources. In the single-source configuration, a point source was positioned at (18°, 72°), while the dual-source scenario involved two equal-intensity coherent point sources located at (54°, 72°) and (54°, 126°) respectively. The measurement array, designed in accordance with compressed sensing theory, adopted a sector-wheel topology comprising 18 microphones strategically distributed across nodal positions within a 0.4 m×0.4 m rectangular region. This array maintained parallelism with the focal plane at a standoff distance of 0.5 m. The focal plane itself was discretized into a 6×21 grid of reconstruction points. To enhance practical relevance, Gaussian white noise was introduced into the acoustic pressure measurements, achieving a SNR of 30 dB through controlled additive noise implementation.

Fig. 2 presents the sound source identification results of the conventional IRL1 algorithm and the enhanced IRL1 method at a source frequency of 1000 Hz. In the figure, “○” denotes actual source positions, while “*” represents identified source locations.

Fig. 2Source localization performance at a sound source frequency of 1000 Hz

a) IRL1

b) Improved IRL1

c) IRL1

d) Improved IRL1

As shown in Fig. 2 for the 1000 Hz acoustic source frequency, the conventional IRL1 algorithm in Fig. 2(a) and 2(c) yields completely erroneous identification results. In contrast, under both single-source (Fig. 2(b)) and dual-source (Fig. 2(d)) conditions, our enhanced IRL1 algorithm successfully localizes all acoustic sources with positional accuracy.

When the source frequency increases to 1500 Hz, Fig. 3(a) and Fig. 3(c) reveal that while the conventional IRL1 algorithm shows improved localization performance compared to lower frequencies, it still fails to correctly identify the single-source position. In dual-source scenarios, it accurately locates the right-side source but exhibits significant positional error for the left-side source. Conversely, our enhanced IRL1 algorithm (Fig. 3(b) and 3(d)) demonstrates robust performance by accurately identifying all true source positions under both test conditions.

Fig. 3Source localization performance at a sound source frequency of 1500 Hz

a) IRL1

b) Improved IRL1

c) IRL1

d) Improved IRL1

3.2. Performance analysis of algorithms under low SNR condition

Maintaining the original simulation conditions, simulation tests were conducted with the frequency set to 2000 Hz and SNR at 5 dB.

From Fig. 4(a), it can be observed that the conventional IRL1 algorithm accurately identifies the left-side acoustic source but still exhibits significant positional error in localizing the right-side source. In Fig. 4(b), the enhanced IRL1 algorithm proposed in this study achieves precise localization of both acoustic sources at their true positions.

Fig. 4Source localization effect diagram of two algorithms at a SNR of 5 dB

a) IRL1

b) Improved IRL1