Two-dimensional multi-frequency synchronous grid-free compressive beamforming

Abstract

Two-dimensional grid-free compressive beamforming with rectangular planar microphone arrays has attracted considerable attention due to its wide source identification region (the hemispherical space in front of arrays), high accuracy (overcoming basis mismatch), and robust performance. However, the uniform distribution of microphones tends to result in spatial aliasing, i.e., the existing single-frequency methods fail to achieve accurate localization for the sources with wavelengths shorter than twice the minimum spacing between microphones. This paper proposes a two-dimensional multi-frequency synchronous grid-free compressive beamforming method, aimed at mitigating spatial aliasing and extending its applicability beyond the sampling limit imposed by the array geometry. First, a two-dimensional multi-frequency synchronous atomic norm minimization problem is formulated and then converted into an equivalent dual problem. Subsequently, a spatial spectral function is constructed from the solution of the dual problem to extract source positions and strengths. Simulation and experiment results show that when the spatial Nyquist criterion is satisfied, both the proposed multi-frequency synchronous method and existing single-frequency methods can accurately identify sources, with the proposed method reducing strength quantification error by over 80% under strong interference. When the spatial Nyquist criterion is violated, the proposed method can effectively suppress spatial aliasing inherent in the existing single-frequency methods and accurately identify common engineering noise sources with high-frequency components, such as electric motor noise, gear whine, and bearing noise.

1. Introduction

Compressive beamforming (CB)^1,2 enables clear and intuitive acoustic imaging, shows strong anti-interference capability and applicability to both coherent and incoherent sources, and thus has great application value in engineering fields such as aerospace,^3–6 automobile,^7,8 wind power generation^9–11 and so on.

Based on the mathematical models, CB can be categorized into four types: on-fixed-grid methods, off-fixed-grid methods, dynamic-grid methods, and grid-free methods. ¹² On-fixed-grid CB discretizes the imaging region into a set of fixed grid points, constructs a sensing matrix using the transfer functions between the grid points and microphones, and then formulates an underdetermined equation system relating the measured pressures to the sensing matrix and the source distribution. Sparse constraints and sparse recovery algorithms are applied to solve this equation system. When acoustic sources fall on the preset grid points, on-fixed-grid CB performs well. However, when sources deviate from the grid points, these methods localize sources to the nearby grid points, resulting in significant performance degradation. This is known as basis mismatch problem. ² The other three types of methods are established to overcome this problem. Off-fixed-grid CB also discretizes the imaging area but mitigates basis mismatch through initial on-grid localization and off-grid offset compensation.^13–15 Dynamic-grid CB discretizes the imaging region into dynamic grids and drives these grids to move toward true source positions, which performs well with an appropriate initial grid setting. ¹⁶ Different from the above two types of CB methods, grid-free CB treats the source distribution as a continuum, thereby avoiding basis mismatch at the root. ¹⁷

Many researchers have made efforts to develop two-dimensional grid-free CB methods based on planar microphone arrays, which enables the identification of the sources lying in the hemispherical space in front of the arrays. In 2017, Yang et al. proposed a two-dimensional grid-free CB method based on atomic norm minimization (ANM), ¹⁷ and later introduced iterative reweighted ANM to improve spatial resolution.^18,19 In 2019, Yang et al. exploited multiple-snapshot measurements to enhance the performance ²⁰ and developed an efficient solver based on the alternating direction method of multipliers (ADMM). ²¹ Tang et al. solved the ANM problem in the dual domain by defining a spatial spectral function related to the dual polynomial. ²² In 2020, Liu et al. proposed an iterative Vandermonde decomposition and shrinkage thresholding-based two-dimensional grid-free CB method, ²³ which solves its model using accelerated proximal mapping and alternating projections. This method avoids performance degradation caused by inaccurate SNR estimation.

The core of the above two-dimensional grid-free CB methods is to formulate and solve the ANM model, for which the construction of a two-level Toeplitz matrix is essential. Therefore, the existing two-dimensional grid-free CB methods can be only applied to regular array configurations, such as uniform or sparse rectangular microphone arrays. ²⁴ For a given number of microphones and equivalent array diameters, regular microphone arrays have a larger minimum spacing between microphones. According to the spatial Nyquist criterion,^25,26 when the minimum spacing between microphones exceeds half the wavelength, DOA estimation may become ambiguous, i.e., high-frequency spatial aliasing occurs. The existing two-dimensional grid-free CB methods process signals at a single frequency, and therefore, their upper applicable frequency limit is constrained by the spatial Nyquist criterion.

Spatial aliasing is frequency-dependent, meaning that aliasing artifacts appear at different locations across frequencies, while the true source locations remain frequency-independent. Therefore, by establishing a multi-frequency model in which the frequency separation (measured in wavelength) exceeds twice the minimum microphone spacing, the intersection of the source distribution sparsity profile across frequencies contains exclusively the true DOAs, enabling aliasing-free DOA estimation. ²⁷ Based on that, some researchers developed a joint sparse multi-frequency model with one-dimensional linear arrays^28–30 to integrate acoustic information at different frequencies into a unified optimization framework, achieving spatial aliasing suppression in one-dimensional line spectral estimation. Inspired by these reports, this paper proposes a two-dimensional multi-frequency synchronous grid-free CB method. The proposed method aims to suppress high-frequency spatial aliasing caused by the use of regular microphone arrays and effectively identify sources whose wavelengths are smaller than twice the minimum microphone spacing, thereby extending its applicability to identifying engineering noise containing high-frequency components, such as electric motor noise, gear whine, and bearing noise. The proposed method first formulates a two-dimensional multi-frequency model and defines the atomic norm. Then, it constructs an infinite-dimensional ANM problem and converts it to a dual problem with finite-dimensional variables for solving. Finally, a spatial spectral function related to the dual polynomial is used to recover the primal solution, from which DOAs and source strengths are extracted.

The main contributions of this article are as follows:

● A two-dimensional multi-frequency synchronous grid-free CB model is established based on atomic norm minimization for source identification in two-dimensional space.

● A multi-frequency spatial spectral function is proposed to extract source information.

● Compared with the existing two-dimensional grid-free CB methods, the proposed method can effectively suppress spatial aliasing and accurately identify the sources with wavelengths shorter than twice the minimum spacing between microphones due to its multi-frequency synchronous model. Compared with the existing multi-frequency CB methods, the proposed method is the two-dimensional extension of these methods, and can achieve two-dimensional source localization in continuous space.

The remainder of this paper is organized as follows: Section 2 introduces the theory of the proposed two-dimensional multi-frequency synchronous grid-free CB, Section 3 evaluates the performance of the proposed method via simulations and Monte Carlo trials, Section 4 validates the correctness of simulations and the effectiveness of the proposed method via loudspeaker identification experiments, and Section 5 summarizes the paper.

2. Theory of two-dimensional multi-frequency synchronous grid-free compressive beamforming

2.1 Two-dimensional multi-frequency synchronous model

Figure 1 shows the measurement model. Sources are located in the far field and radiate plane waves. ${\Omega }_{\mathrm{S}i}=\left({\theta }_{\mathrm{S}i},{\varphi }_{\mathrm{S}i}\right)$ denotes the DOA of the $i\hspace{0.33em}\text{th}$ source, where ${\theta }_{\mathrm{S}i}\in \left[0°,90°\right]$ and ${\varphi }_{\mathrm{S}i}\in [0°,360°)$ represent elevation and azimuth angles, respectively. A rectangular microphone array, which contains A rows and B columns, i.e., A total of AB microphones are used. The symbol “●” denotes microphones. $\left(a,b\right)\in \{0,1,\cdots ,\left(A-1\right)\}\times \{0,1,\cdots ,\left(B-1\right)\}$ denotes the index of microphones. $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$ represent the distances between the adjacent microphones in the x and y directions, respectively. The $\left(0,0\right)\hspace{0.33em}\text{th}$ microphone is located at the origin. $i=1,2,\cdots ,I$ is the index of sources, and I is the number of sources. $l=1,2,\cdots ,L$ denotes the index of frequencies, and L is the number of frequencies. $\mathbb{C}$ denotes the set of complex numbers. ${\mathbf{s}}_i=\left[s_{i,1},s_{i,2},\cdots ,s_{i,L}\right]\in {\mathbb{C}}^{1\times L}$ is a row vector consisting of the pressures at all frequencies radiated by the $i\hspace{0.33em}\text{th}$ source at the origin. Let $t_1\equiv \sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi \mathrm{\Delta }x$ and $t_2\equiv \sin \hspace{0.17em}\theta \hspace{0.17em}\sin \hspace{0.17em}\varphi \mathrm{\Delta }y$ . The continuous expression of multi-frequency signal $\mathbf{x}\left(t_1,t_2\right)\in {\mathbb{C}}^{1\times L}$ is given by

$\mathbf{x}\left(t_1,t_2\right)={\displaystyle \sum }_{i=1}^I{\mathbf{s}}_i\mathrm{\delta }\left(t_1-t_{1i},t_2-t_{2i}\right),$

(1)

where $t_{1i}\equiv \sin \hspace{0.17em}{\theta }_{\mathrm{S}i}\hspace{0.17em}\cos \hspace{0.17em}{\varphi }_{\mathrm{S}i}\mathrm{\Delta }x$ , $t_{2i}\equiv \sin \hspace{0.17em}{\theta }_{\mathrm{S}i}\hspace{0.17em}\sin \hspace{0.17em}{\varphi }_{\mathrm{S}i}\mathrm{\Delta }y$ , and $\mathrm{\delta }(·)$ is the Dirac Delta function. ${\mathbf{p}}_{a,b}=\left[p_{a,b,1},p_{a,b,2},\cdots ,p_{a,b,L}\right]\in {\mathbb{C}}^{1\times L}$ assembles the pressures at the $\left(a,b\right)\hspace{0.33em}\text{th}$ microphone across all frequencies, and it can be written as

${\mathbf{p}}_{a,b}=\underset{\mathbb{T}}{\iint }\mathbf{x}(t_1,t_2)\circ {\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{\mathbf{f}}{c}(t_{1}a+t_{2}b)}\mathrm{d}{t}_1\mathrm{d}{t}_2=\underset{i=1}{\overset{I}{}}[{\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{f_1}{c}(t_{1i}a+t_{2i}b)},{\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{f_2}{c}(t_{1i}a+t_{2i}b)},\cdots ,{\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{f_L}{c}(t_{1i}a+t_{2i}b)}]{\mathbf{S}}_i,$

(2)

where $\mathbb{T}=\left\{\left(t_1,t_2\right)\right|{\left(t_1/\mathrm{\Delta }x\right)}^2+{\left(t_2/\mathrm{\Delta }y\right)}^2\le 1\},$ the symbol “ $\circ$ ” denotes the Hadamard product, $c$ is the sound speed and $\mathrm{j}=\sqrt{-1}$ is the imaginary unit. The center frequencies of all subbands are collected in a row vector $\mathbf{f}=\left[f_1,f_2,\cdots ,f_L\right]\in {\mathbb{R}}^{1\times L}$ . $\mathbb{R}$ denotes the set of real numbers. ${\mathbf{S}}_i=\text{Diag}\left({\mathbf{s}}_i\right)\in {\mathbb{C}}^{L\times L}$ , where $\text{Diag}(·)$ forms a diagonal matrix with the vector in parentheses as the diagonal. We construct $\mathbf{D}\left(t_{1i},t_{2i},\mathbf{f}\right)\in {\mathbb{C}}^{AB\times L}$ ,

$\mathbf{D}(t_{1i},t_{2i},\mathbf{f})=\mathbf{D}(t_{1i},\mathbf{f})\odot \mathbf{D}(t_{2i},\mathbf{f})=\left[\begin{array}{c}1\\ {\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{\mathbf{f}}{c}{t}_{1i}}\\ ⋮\\ {\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{\mathbf{f}}{c}{t}_{1i}(A-1)}\end{array}\right]\odot \left[\begin{array}{c}1\\ {\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{\mathbf{f}}{c}{t}_{2i}}\\ ⋮\\ {\mathrm{e}}^{\mathrm{j}2\mathrm{\pi }\frac{\mathbf{f}}{c}{t}_{2i}(B-1)}\end{array}\right],$

(3)

where the symbol “ $\odot$ ” represents the Khatri-Rao product. For brevity, $\mathbf{D}\left(t_{1i},t_{2i},\mathbf{f}\right)$ is expressed as ${\mathbf{D}}_i$ , ${\mathbf{D}}_i=\left[{\mathbf{d}}_{i,1},{\mathbf{d}}_{i,2},\cdots ,{\mathbf{d}}_{i,L}\right]$ , and ${\mathbf{d}}_{i,l}\in {\mathbb{C}}^{AB\times 1}$ is the $l\hspace{0.33em}\text{th}$ column of ${\mathbf{D}}_i$ . Let $\mathbf{P}={\left[{\mathbf{p}}_{0,0}^{\mathrm{T}},{\mathbf{p}}_{0,1}^{\mathrm{T}},\cdots ,{\mathbf{p}}_{A-1,B-1}^{\mathrm{T}}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{AB\times L}$ , and ${(·)}^{\mathrm{T}}$ denotes the transpose. $\mathbf{P}$ can be expressed as

$\mathbf{P}={\displaystyle \sum }_{i=1}^I{\mathbf{D}}_i{\mathbf{S}}_i=\mathbf{D}\mathbf{S},$

(4)

where $\mathbf{D}=\left[{\mathbf{D}}_1,{\mathbf{D}}_2,\cdots ,{\mathbf{D}}_I\right]\in {\mathbb{C}}^{AB\times LI}$ and $\mathbf{S}={\left[{\mathbf{S}}_1^{\mathrm{T}},{\mathbf{S}}_2^{\mathrm{T}},\cdots ,{\mathbf{S}}_I^{\mathrm{T}}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{LI\times L}$ . When considering the noise inference, the measured signal is

${\mathbf{P}}^{\#}=\mathbf{P}+\mathbf{N},$

(5)

where $\mathbf{N}\in {\mathbb{C}}^{AB\times L}$ is the noise matrix. The signal-to-noise (SNR) ratio is defined as $\text{SNR}=20\hspace{0.17em}{\log }_{10}\left({\Vert \mathbf{P}\Vert }_{\mathrm{F}}/{\Vert \mathbf{N}\Vert }_{\mathrm{F}}\right)$ , where ${\Vert ·\Vert }_{\mathrm{F}}$ denotes the Frobenius norm.

Figure 1.

Measurement model.

2.2 Primal problem

In the continuous spatial domain, the atomic norm of $\mathbf{S}$ is defined as

${\Vert \mathbf{S}\Vert }_{\mathcal{A}}={\displaystyle \sum }_{i=1}^I\text{tr}\left(\left[\begin{array}{cccc}\left|s_{i,1}\right|& 0& \cdots & 0\\ 0& \left|s_{i,2}\right|& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \left|s_{i,L}\right|\end{array}\right]\right),$

(6)

where $\text{tr}(·)$ denotes the trace of a matrix, $|·|$ denotes to take the absolute value. ${\Vert \mathbf{S}\Vert }_{\mathcal{A}}$ measures the sparsity of source distribution. The DOA estimation in the continuous spatial domain can be achieved by solving the following convex optimization problem:

$\underset{\mathbf{x}}{\min }\hspace{0.33em}{\Vert \mathbf{S}\Vert }_{\mathcal{A}}\hspace{1em}\mathrm{s}.\mathrm{t}.\hspace{0.17em}{\Vert {\mathbf{P}}^{\#}-\mathbf{D}\mathbf{S}\Vert }_{\mathrm{F}}\le \epsilon ,$

(7)

where $\epsilon$ is the upper bound of the Frobenius norm of noise. Since $\mathbf{x}\left(t_1,t_2\right)$ is a variable in the continuous domain, the primal problem is an infinite-dimensional optimization problem and cannot be directly solved. In the next section, we will solve Eq. (7) by solving its corresponding dual problem.

2.3 Dual problem

Firstly, the Lagrange function of Eq. (7) is constructed,

$\mathcal{L}\left(\mathbf{S},\mathbf{C},\beta \right)={\Vert \mathbf{S}\Vert }_{\mathcal{A}}+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\left({\mathbf{P}}^{\#}-\mathbf{D}\mathbf{S}-\mathbf{N}\right)\right)\right)+\beta \left({\Vert \mathbf{N}\Vert }_{\mathrm{F}}^2-{\epsilon }^2\right),$

(8)

where $\beta \in {\mathbb{R}}^{+}$ and $\mathbf{C}=\left[{\mathbf{c}}_1,{\mathbf{c}}_2,\cdots ,{\mathbf{c}}_{\mathrm{L}}\right]\in {\mathbb{C}}^{AB\times L}$ are Lagrangian multipliers, ${\mathbb{R}}^{+}$ denotes the set of positive real numbers, ${(·)}^{\mathrm{H}}$ denotes the Hermitian transpose, and $(·)$ takes the real parts of complex numbers. The dual function $g\left(\mathbf{C},\beta \right)$ is the maximum of the lower bound of Eq. (8),

$g\left(\mathbf{C},\beta \right)=\underset{\mathbf{S}}{\inf }\hspace{0.17em}\mathcal{L}\left(\mathbf{S},\mathbf{C},\beta \right)$

$=\underset{\mathbf{S}}{\inf }\hspace{0.33em}\left({\Vert \mathbf{S}\Vert }_{\mathcal{A}}+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\left({\mathbf{P}}^{\#}-\mathbf{D}\mathbf{S}-\mathbf{N}\right)\right)\right)+\beta \left({\Vert \mathbf{N}\Vert }_{\mathrm{F}}^2-{\epsilon }^2\right)\right)=\underset{\mathbf{S}}{\inf }\hspace{0.33em}\left({\Vert \mathbf{S}\Vert }_{\mathcal{A}}+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}-{\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}-{\mathbf{C}}^{\mathrm{H}}\mathbf{N}\right)\right)+\beta \left(\text{tr}\left(\mathbf{N}{\mathbf{N}}^{\mathrm{H}}\right)-{\epsilon }^2\right)\right)$

$=\underset{\mathbf{S}}{\inf }\hspace{0.17em}\left({\Vert \mathit{S}\Vert }_{\mathcal{A}}-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}\right)\right)\right)+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{*}-{\mathbf{C}}^{\mathrm{H}}\mathbf{N}\right)\right)+\beta \left(\text{tr}\left(\mathbf{N}{\mathbf{N}}^{\mathrm{H}}\right)-{\epsilon }^2\right),$

(9)

where $\inf (·)$ denotes the greatest lower bound.

Then, we minimize $g\left(\mathbf{C},\beta \right)$ by taking its partial derivative with respect to the noise matrix and setting the derivative to 0, i.e.,

$\frac{\partial g\left(\mathbf{C},\beta \right)}{\partial \mathbf{N}}=\frac{\partial \left(-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{N}\right)\right)+\beta \text{tr}\left(\mathbf{N}{\mathbf{N}}^{\mathrm{H}}\right)\right)}{\partial \mathbf{N}}=-\mathbf{C}+2\beta \mathbf{N}=0.$

(10)

Therefore, the optimal noise matrix is ${\mathbf{N}}_{\text{opt}}=\mathbf{C}/\left(2\beta \right)$ . Now, we calculate the dual function at ${\mathbf{N}}_{\text{opt}}$ and maximize it with respect to $\beta$ ,

$\frac{{\partial g\left(\mathbf{C},\beta \right)|}_{{\mathbf{N}}_{\text{opt}}}}{\partial \beta }=\frac{{\partial \beta \left(\text{tr}\left(\mathbf{N}{\mathbf{N}}^{\mathrm{H}}\right)-{\epsilon }^2\right)|}_{{\mathbf{N}}_{\text{opt}}}}{\partial \beta }={\left(\text{tr}\left(\mathbf{N}{\mathbf{N}}^{\mathrm{H}}\right)-{\epsilon }^2\right)|}_{{\mathbf{N}}_{\text{opt}}}=\frac{\text{tr}\left(\mathbf{C}{\mathbf{C}}^{\mathrm{H}}\right)}{4{\beta }^2}-{\epsilon }^2=0,$

(11)

thus the optimal dual variable is ${\beta }_{\text{opt}}=\sqrt{\text{tr}\left(\mathbf{C}{\mathbf{C}}^{\mathrm{H}}\right)/\left(4{\epsilon }^2\right)}={\Vert \mathbf{C}\Vert }_{\mathrm{F}}/\left(2\epsilon \right)$ . The dual function at ${\mathbf{N}}_{\text{opt}}$ and ${\beta }_{\text{opt}}$ is

${g\left(\mathbf{C}\right)|}_{{\mathbf{N}}_{\text{opt}},{\beta }_{\text{opt}}}$

$={\underset{\mathbf{S}}{\inf }\hspace{0.17em}\left({\Vert \mathbf{S}\Vert }_{\mathcal{A}}-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}\right)\right)\right)+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}\right)\right)-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{N}\right)\right)+\beta \left(\text{tr}\left(\mathbf{N}{\mathbf{N}}^{\mathrm{H}}\right)-{\epsilon }^2\right)|}_{{\mathbf{N}}_{\text{opt}},{\beta }_{\text{opt}}}$

$={\underset{\mathbf{S}}{\inf }\hspace{0.17em}\left({\Vert \mathbf{S}\Vert }_{\mathcal{A}}-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}\right)\right)\right)+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}\right)\right)-\text{Re}\left(\text{tr}\left(\frac{{\mathbf{C}}^{\mathrm{H}}\mathbf{C}}{2\beta }\right)\right)+\beta \left(\text{tr}\left(\frac{{\mathbf{C}}^{\mathrm{H}}\mathbf{C}}{4{\beta }^2}\right)-{\epsilon }^2\right)|}_{{\beta }_{\text{opt}}}={\underset{\mathbf{S}}{\inf }\hspace{0.17em}\left({\Vert \mathbf{S}\Vert }_{\mathcal{A}}-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}\right)\right)\right)+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}\right)\right)-\frac{{\Vert \mathbf{C}\Vert }_{\mathrm{F}}^2}{2\beta }|}_{{\beta }_{\text{opt}}}$

$=\underset{\mathbf{S}}{\inf }\hspace{0.17em}\left({\Vert \mathbf{S}\Vert }_{\mathcal{A}}-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}\right)\right)\right)+\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}\right)\right)-\epsilon {\Vert \mathbf{C}\Vert }_{\mathrm{F}}.$

(12)

With respect to the first part in Eq. (12), we have

${\Vert \mathbf{S}\Vert }_{\mathcal{A}}-\text{Re}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}\mathbf{D}\mathbf{S}\right)\right)={\Vert \mathbf{S}\Vert }_{\mathcal{A}}-{\sum }_{i=1}^I\text{tr}\left(\text{Re}\left({\left({\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}\right)}^{\mathrm{H}}{\mathbf{S}}_i\right)\right)$

$={\Vert \mathbf{S}\Vert }_{\mathcal{A}}-{\displaystyle \sum }_{i=1}^I\text{tr}\left(\left|{\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}\right|\circ \left|{\mathbf{S}}_i\right|\circ {\cos \hspace{0.17em}\mathbf{\phi }}_i\right)$

$={\displaystyle \sum }_{i=1}^I\text{tr}\left(\left|{\mathbf{S}}_i\right|-\left|{\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}\right|\circ \left|{\mathbf{S}}_i\right|\circ {\cos \hspace{0.17em}\mathbf{\phi }}_i\right)$

$={\sum }_{i=1}^I\text{tr}\left(\left|{\mathbf{S}}_i\right|\circ \left({\mathbf{I}}_L-\left|{\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}\right|\circ {\cos \hspace{0.17em}\mathbf{\phi }}_i\right)\right)$

$\ge {\sum }_{i=1}^I\text{tr}\left(\left|{\mathbf{S}}_i\right|\circ \left({\mathbf{I}}_L-\left|{\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}\right|\right)\right)$

$={\displaystyle \sum }_{i=1}^I{\displaystyle \sum }_{l=1}^L\left(\left|s_{i,l}\right|\left(1-\left|{\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l\right|\right)\right)$

(13)

where $\left|\mathbf{X}\right|$ represents the modulus operation of each element in the matrix $\mathbf{X}$ , ${\mathbf{I}}_L$ denotes a $L\times L$ unit matrix, and ${\mathbf{\phi }}_i$ is the element-wise included angle between ${\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}$ and ${\mathbf{S}}_i$ . When $\left|{\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l\right|\le 1$ , Eq. (13) is nonnegative, implying that its lower bound is 0. When $\left|{\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l\right|>1$ , the lower bound is $-\infty$ . Therefore, the dual function in Eq. (12) can be rewritten as

${g\left(\mathbf{C}\right)|}_{{\mathbf{N}}_{\text{opt}},{\beta }_{\text{opt}}}=\{\begin{array}{l}\mathrm{Re}\left(\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}\right)\right)\right)-\epsilon {\Vert \mathbf{C}\Vert }_{\mathrm{F}},\\ -\infty ,\end{array}\begin{array}{l}\left|{\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l\right|\le 1\\ \text{otherwise}\end{array}.$

(14)

It is known that ${\mathbf{d}}_{i,l}^{\mathrm{H}}\left({\mathbf{c}}_1{\mathbf{c}}_1^{\mathrm{H}}+{\mathbf{c}}_2{\mathbf{c}}_2^{\mathrm{H}}+\cdots +{\mathbf{c}}_L{\mathbf{c}}_L^{\mathrm{H}}\right){\mathbf{d}}_{i,l}$ is the $l\hspace{0.33em}\text{th}$ diagonal element of ${\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}{\mathbf{D}}_i$ , and

${\left|{\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l\right|}^2={\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l{\mathit{c}}_l^{\mathrm{H}}{\mathbf{d}}_{i,l}\le {\mathbf{d}}_{i,l}^{\mathrm{H}}({\mathbf{c}}_1{\mathbf{c}}_1^{\mathrm{H}}+{\mathbf{c}}_2{\mathbf{c}}_2^{\mathrm{H}}+\cdots +{\mathbf{c}}_L{\mathbf{c}}_L^{\mathrm{H}}){\mathbf{d}}_{i,l}.$

(15)

Therefore, the constraint $\left|{\mathbf{d}}_{i,l}^{\mathrm{H}}{\mathbf{c}}_l\right|\le 1$ in Eq. (14) can be replaced by ${\mathbf{d}}_{i,l}^{\mathrm{H}}\left({\mathbf{c}}_1{\mathbf{c}}_1^{\mathrm{H}}+{\mathbf{c}}_2{\mathbf{c}}_2^{\mathrm{H}}+\cdots +{\mathbf{c}}_L{\mathbf{c}}_L^{\mathrm{H}}\right){\mathbf{d}}_{i,l}\le 1$ .

We construct a Hermitian matrix $\mathbf{H}\in {\mathbb{C}}^{AB\times AB}$ , which satisfies $\mathbf{d}{\left(t_1,t_2,f_l\right)}^H\mathbf{H}\mathbf{d}\left(t_1,t_2,f_l\right)=1$ . Since

$\mathbf{d}{\left(t_1,t_2,f_l\right)}^{\mathrm{H}}\mathbf{H}\mathbf{d}\left(t_1,t_2,f_l\right)=\text{sum}\left\{\left(\mathbf{d}{\left(t_1,t_2,f_l\right)}^{*}\mathbf{d}{\left(t_1,t_2,f_l\right)}^{\mathrm{T}}\right)\circ \mathbf{H}\right\}=\text{sum}\left\{\mathbf{G}\circ \mathbf{H}\right\},$

(16)

where $\text{sum}\{·\}$ takes the sum of all elements in vectors or matrices, ${(·)}^{*}$ denotes the conjugate, $\mathbf{G}=\mathbf{d}{\left(t_1,t_2,f_l\right)}^{*}\mathbf{d}{\left(t_1,t_2,f_l\right)}^{\mathrm{T}}\in {\mathbb{C}}^{AB\times AB}$ is a two-fold Toeplitz matrix with the dimension of $AB\times AB$ , whose elements are all in the exponential form,

$\mathbf{G}=\left[\begin{array}{cccc}{\mathbf{G}}_0& {\mathbf{G}}_1^{\mathrm{H}}& \cdots & \hspace{0.17em}{\mathbf{G}}_{A-1}^{\mathrm{H}}\\ {\mathbf{G}}_1& {\mathbf{G}}_0& \cdots & {\mathbf{G}}_{A-2}^{\mathrm{H}}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{G}}_{A-1}& {\mathbf{G}}_{A-2}& \cdots & {\mathbf{G}}_0\end{array}\right].$

(17)

Each block ${\mathbf{G}}_m,m=0,1,\cdots ,\left(A-1\right)$ is a Toeplitz matrix with the dimension of $B\times B$ . All diagonal elements of ${\mathbf{G}}_0$ are 1. Accordingly, we partition $\mathbf{H}$ into blocks, i.e.,

$\mathbf{H}=\left[\begin{array}{cccc}{\mathbf{H}}_{0,0}& {\mathbf{H}}_{0,1}& \cdots & \hspace{0.17em}{\mathbf{H}}_{0,A-1}\\ {\mathbf{H}}_{0,1}^{\mathrm{H}}& {\mathbf{H}}_{1,1}& \cdots & {\mathbf{H}}_{1,A-1}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{H}}_{0,A-1}^{\mathrm{H}}& {\mathbf{H}}_{1,A-1}^{\mathrm{H}}& \cdots & {\mathbf{H}}_{A-1,A-1}\end{array}\right].$

(18)

Let ${\mathbf{J}}_m={\displaystyle \sum }_{j=0}^{A-1-m}{\mathbf{H}}_{j,j+m},m=0,1,\cdots ,\left(A-1\right)$ . Since $\text{sum}\left\{\mathbf{G}\circ \mathbf{H}\right\}=1$ , ${\mathbf{J}}_m$ satisfies that $\text{sum}\left\{\text{diag}\left({\mathbf{J}}_0,0\right)\right\}=1$ , $\text{sum}\left\{\text{diag}\left({\mathbf{J}}_0,n\right)\right\}=0,n=1,2,\cdots ,\left(B-1\right)$ , and $\text{sum}\{\text{diag}\left({\mathbf{J}}_{\mathrm{m}},n)\right\}=0,m=1,2,\cdots ,\left(A-1\right),n=-\left(B-1\right),\cdots ,-1,0,1,\cdots ,\left(B-1\right).$ $\text{diag}\left(·,n\right)$ takes the $n\hspace{0.33em}\text{th}$ diagonal of matrices and forms a column vector. $n$ represents the offset of diagonals. When $n=0$ , it refers to the principal diagonal. When $n>0$ , it denotes the $n\hspace{0.33em}\text{th}$ minor diagonal above the principal diagonal, and when $n<0$ , it is the $\left|n\right|\hspace{0.33em}\text{th}$ minor diagonal below the principal diagonal.

For an arbitrary $i$ and $l$ , the diagonal elements of ${\mathbf{D}}_i^{\mathrm{H}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}{\mathbf{D}}_i$ are constrained to be no larger than 1, i.e.,

$\mathbf{d}{\left(t_1,t_2,f_l\right)}^{\mathrm{H}}\mathbf{C}{\mathbf{C}}^{\mathrm{H}}\mathbf{d}\left(t_1,t_2,f_l\right)\le \mathbf{d}{\left(t_1,t_2,f_l\right)}^{\mathrm{H}}\mathbf{H}\mathbf{d}\left(t_1,t_2,f_l\right).$

(19)

According to the Schur complement condition, we have

$\left[\begin{array}{ll}\mathbf{H}& \mathbf{C}\\ {\mathbf{C}}^{\mathrm{H}}& {\mathbf{I}}_{\mathrm{L}}\end{array}\right]\succcurlyeq 0,$

(20)

where $\succcurlyeq 0$ denotes that the matrix is positive semidefinite.

Overall, the primal minimization problem (Eq. (7)) is equivalent to the maximization of the dual function (Eq (14).), and it can be rewritten as

$\left\{\widehat{\mathbf{C}},\widehat{\mathbf{H}}\right\}=\arg \hspace{0.17em}\max \hspace{0.17em}\left(\text{tr}\left({\mathbf{C}}^{\mathrm{H}}{\mathbf{P}}^{\#}\right)\right)-\epsilon {\Vert \mathbf{C}\Vert }_{\mathrm{F}}$

$\mathrm{s}.\mathrm{t}.\hspace{0.17em}\left[\begin{array}{ll}\mathbf{H}& \mathbf{C}\\ {\mathbf{C}}^{\mathrm{H}}& {\mathbf{I}}_{\mathrm{L}}\end{array}\right]\succcurlyeq 0,\mathbf{H}=\left[\begin{array}{cccc}{\mathbf{H}}_{0,0}& {\mathbf{H}}_{0,1}& \cdots & \hspace{0.17em}{\mathbf{H}}_{0,A-1}\\ {\mathbf{H}}_{0,1}^{\mathrm{H}}& {\mathbf{H}}_{1,1}& \cdots & {\mathbf{H}}_{1,A-1}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{H}}_{0,A-1}^{\mathrm{H}}& {\mathbf{H}}_{1,A-1}^{\mathrm{H}}& \cdots & {\mathbf{H}}_{A-1,A-1}\end{array}\right]$

${\mathit{J}}_m={\displaystyle \sum }_{j=0}^{A-1-m}{\mathbf{H}}_{j,j+m},m=0,1,\cdots ,\left(A-1\right)$

$\text{sum}\left\{\text{diag}\left({\mathbf{J}}_0,0\right)\right\}=1\text{sum}\left\{\text{diag}\left({\mathbf{J}}_0,n\right)\right\}=0|n=1,2,\cdots ,\left(B-1\right)$

$\text{sum}\left\{\text{diag}\left({\mathbf{J}}_m,n\right)\right\}=0|m=1,2,\cdots ,\left(A-1\right),n=-\left(B-1\right),\cdots ,-1,0,1,\cdots ,\left(B-1\right)$

(21)

This positive semidefinite programming problem can be solved by the SDPT3 solver in the CVX toolbox.

2.4 DOA and strength estimation

Define the spatial spectral function as

$y\left(\theta ,\varphi \right)={\displaystyle \sum }_{l=1}^L\mathbf{d}{\left(t_1,t_2,f_l\right)}^{\mathrm{H}}{\widehat{\mathbf{c}}}_l{\widehat{\mathbf{c}}}_l^{\mathrm{H}}\mathbf{d}\left(t_1,t_2,f_l\right),$

(22)

where ${\widehat{\mathbf{c}}}_l$ is the $l\hspace{0.33em}\text{th}$ column of $\widehat{\mathbf{C}}$ . According to the relationship between the primal problem and the dual problem, the peak positions of $y\left(\theta ,\varphi \right)$ correspond to source DOAs. The number of the peaks whose amplitudes are larger than $0.7\hspace{0.17em}\max \hspace{0.17em}\{y\left(\theta ,\varphi \right)\}$ in $y\left(\theta ,\varphi \right)$ is denoted by $\widehat{k}$ . The locations of these $\widehat{k}$ peaks are the estimated DOAs $\left\{{\widehat{\mathrm{\Omega }}}_{Si}=\left({\widehat{\theta }}_{Si},{\widehat{\varphi }}_{Si}\right)|i=1,2,\cdots ,\widehat{k}\right\}$ .

With the estimated DOAs, we obtain $t_{1i}\equiv \sin \hspace{0.17em}{\widehat{\theta }}_{Si}\hspace{0.17em}\cos \hspace{0.17em}{\widehat{\varphi }}_{Si}\mathrm{\Delta }x$ and $t_{2i}\equiv \sin \hspace{0.17em}{\widehat{\theta }}_{Si}\hspace{0.17em}\cos \hspace{0.17em}{\widehat{\varphi }}_{Si}\mathrm{\Delta }y$ . Then, we construct the sensing matrix $\mathbf{A}=\left[\mathbf{D}\left(t_{11},t_{21},\mathbf{f}\right),\mathbf{D}\left(t_{12},t_{22},\mathbf{f}\right),\cdots ,\mathbf{D}\left(t_{1\widehat{k}},t_{2\widehat{k}},\mathbf{f}\right)\right]\in {\mathbb{C}}^{AB\times \widehat{k}L}$ and quantify the strengths using the least-square method,

$\widehat{\mathbf{S}}={\mathbf{A}}^{+}{\mathbf{P}}^{\#},$

(23)

where ${(\hspace{0.17em})}^{+}$ denotes the pseudo-inverse.

3. Simulation

A rectangular microphone array consisting of 8 rows and 8 columns (a total of 64 microphones) is used, with a spacing of $\mathrm{\Delta }x=\mathrm{\Delta }y=0.035\hspace{0.17em}\mathrm{m}$ between adjacent microphones. The upper frequency limit for aliasing-free sampling is approximately 4900 Hz. With reference to the existing two-dimensional single-frequency grid-free compressive beamforming method in Refs. ¹⁷ , this section analyzes the performance of the proposed method when the source frequency is below 4900 Hz (satisfying the spatial Nyquist criterion and free of spatial aliasing) and exceeds 4900 Hz (spatial aliasing occurs). We conduct D Monte-Carlo trials, $\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}$ denotes the DOA estimation error of the $i\hspace{0.33em}\text{th}$ source in the $d\hspace{0.33em}\text{th}$ trial,

$\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}=\frac{180}{\mathrm{\pi }}\arccos \left(\cos \hspace{0.17em}{\widehat{\theta }}_{\mathrm{S}i,d}\hspace{0.17em}\cos \hspace{0.17em}{\theta }_{\mathrm{S}i,d}+\cos \left({\widehat{\varphi }}_{\mathrm{S}i,d}-{\varphi }_{\mathrm{S}i,d}\right)\sin \hspace{0.17em}{\widehat{\theta }}_{\mathrm{S}i,d}\hspace{0.17em}\sin \hspace{0.17em}{\theta }_{\mathrm{S}i,d}\right).$

(24)

We also define the probability of accurate identification (P), the DOA estimation root mean square error (RMSE), and the strength quantification mean error (ME) as

$P=\frac{\text{size}\left(\left[\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}|\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}\le 5°\right]\right)}{ID},$

(25)

$\text{RMSE}=\frac{{\Vert \left[\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}|\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}\le 5°\right]\Vert }_2}{\sqrt{\text{size}\left(\left[\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}|\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}\le 5°\right]\right)}},$

(26)

$\text{ME}=\frac{\sum \left|\left[20\text{lo}{\mathrm{g}}_{10}\left(\left|{\widehat{s}}_{i,d}\right|/\left|s_{i,d}\right|\right)|\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}\le 5°\right]\right|}{\text{size}\left(\left[\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}|\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}\le 5°\right]\right)},$

(27)

where ${\Vert ·\Vert }_2$ denotes ${\mathcal{l}}_2$ norm, $\text{size}(·)$ represents the size of vectors in parentheses, and $\left[\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}|\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}\le 5°\right]$ denotes a column vector consisting of all $\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{S}i,d}$ which satisfy the condition after the vertical bar.

3.1 Aliasing-free case

Assume three sources. They radiate sound waves with a frequency of 3000 Hz. The DOAs ( ${\mathrm{\Omega }}_{\mathrm{S}i}=\left({\theta }_{\mathrm{S}i},{\varphi }_{\mathrm{S}i}\right)$ ) are $\left(30°,200°\right)$ , $\left(50°,90°\right)$ , and $\left(80°,150°\right)$ . The strengths are 100 dB, 94 dB, and 97 dB (referenced to $2.0\times 10^{-5}\hspace{0.17em}\text{Pa}$ ). When simulating the pressures measured by the microphone array, white Gaussian noise with an SNR of 20 dB is added to the clean signal. The proposed method and the method in Ref. ¹⁷ are used to process the pressures at 3000 Hz. Figure 2(a) shows the spatial spectral function of the proposed method. Figure 2(b) and 2(c) present the source identification results. The spatial spectral function exhibits the peaks with nearly unit amplitudes at locations corresponding to the true source DOAs. Both of the two methods accurately localize the sources. RMSE are $0.59°$ and $0.84°$ , respectively. ME are 0.12 dB and 1.49 dB. These two methods achieve comparable performance in DOA estimation, but the proposed method demonstrates a slightly better performance in strength quantification.

Figure 2.

Identification results of the proposed method (single-frequency modelling) and the method in Ref. ¹⁷ under aliasing-free condition.

We vary SNR and perform Monte Carlo trials. The frequency is set to 3000 Hz. The SNR ranges from 5 dB to 30 dB with a step of 5 dB. At each SNR, 1000 trials are considered. In each trial, four sources are assumed, whose DOAs are randomly generated while satisfying the minimum separation condition. ¹⁸ The strengths are set to unit amplitude with random phases. Figure 3 shows P, RMSE, ME, and the average runtime for the two methods. As shown in Figure 3(a), the proposed method achieves a higher probability of accurate identification than the method in Ref. ¹⁷ across different SNRs. From Figure 3(b) and 3(c), it can be observed that both methods yield comparable DOA estimation accuracy, whereas the proposed method significantly outperforms the method in Ref. ¹⁷ in terms of strength quantification. The advantage becomes more pronounced as the SNR decreases. At SNRs of 10 dB and 5 dB, the strength quantification error is reduced by over 80% and 90%, respectively. In Figure 3(d), the runtimes of the two methods are nearly identical, since they operate on matrices of the same dimension and employ the same solver.

Figure 3.

Performance comparison between the proposed method (single-frequency modelling) and the method in Ref. ¹⁷ under aliasing-free condition at different SNRs.

3.2 Spatial aliasing case

In the case described in Section 3.1, the source frequency is changed to 7000 Hz, while all other conditions remain unchanged. The results are shown in Figure 4. Figure 4(a) presents the spatial spectral function, which produces peaks with high amplitudes at the locations where there is no source due to spatial aliasing. The imaging maps are shown in Figure 4(b) and 4(c). The proposed method identifies some spurious sources, whereas the method in Ref. ¹⁷ could only localize a part of sources. The presence of these spurious sources adversely affects the subsequent strength quantification.

Figure 4.

Identification results of the proposed method (single-frequency modelling) and the method in Ref. ¹⁷ under spatial aliasing.

Assume that the three sources contain spectral components not only at 7000 Hz but also at 7500 Hz and 8000 Hz. The source strengths at the frequencies are [100, 94, 94] dB, [94, 100, 94] dB, and [97, 97, 94] dB (referenced to $2.0\times 10^{-5}\hspace{0.17em}\text{Pa}$ ). Other simulation conditions remain unchanged. The proposed method is used to post the pressures at the three frequencies synchronously. Figure 5(a) and 5(b) present the spatial spectral function and the imaging map, respectively. Compared with Figure 4(a) and 4(b), we find that spurious peaks are effectively suppressed, and sources are accurately identified.

Figure 5.

Identification results of the proposed method under spatial aliasing with multi-frequency synchronous modeling.

SNR estimation is required by CVX toolbox to solve the dual problem. To investigate its influence on the performance of the proposed method, we underestimate and overestimate SNR to 10 dB and 30 dB (true SNR is 20 dB), and examine the performance based on the case in Figure 5. The results under inaccurate SNR estimation are shown in Figure 6. Comparing Figure 5(b) and Figure 6, we observe that underestimation of SNR has negligible effect on both DOA estimation and strength quantification, whereas overestimation slightly increases DOA estimation error, with minimal impact on strength quantification.

Figure 6.

Identification results of the proposed method under spatial aliasing with multi-frequency synchronous modeling and inaccurate SNR estimation.

We vary frequency and perform Monte Carlo trials. The fundamental frequencies are set to 5000 Hz, 6000 Hz, 7000 Hz, 8000 Hz, and 9000 Hz, respectively. For each fundamental frequency, four additional frequency components are defined, each increasing by 100 Hz from the fundamental one (e.g., when the fundamental frequency is 5000 Hz, the frequency components are [5000, 5100, 5200, 5300, 5400] Hz). The method in Ref. ¹⁷ only processes the pressures at the fundamental frequencies, whereas the proposed method utilizes both fundamental and additional frequency components. At each frequency, 1000 trials are conducted with four sources. DOAs are randomly generated while satisfying the minimum separation condition. Each source is assigned a unit amplitude with a random phase, and white Gaussian noise is added to achieve an SNR of 20 dB. Figure 7 shows P, RMSE, ME, and the average runtime for the two methods. Obviously, when frequency increases, the probability of accurately identifying sources of the single-frequency method in Ref. ¹⁷ drops sharply due to spatial aliasing. In contrast, the proposed method maintains a high probability and accuracy in source identification. A slight performance degradation is observed for the proposed method at higher frequencies, which may be attributed to the reduced relative frequency spacing that causes overlapping peaks at aliasing positions. Since multi-frequency synchronous processing increases the dimensionality of matrices, the runtime correspondingly rises.

Figure 7.

Performance comparison between the proposed method (multi-frequency synchronous modelling) and the method in Ref. ¹⁷ under spatial aliasing at different frequencies.

4. Experiments

To validate the effectiveness of the proposed method and the correctness of the simulation results, a loudspeaker source identification experiment is conducted in a semi-anechoic chamber using a $4\times 4$ rectangular array. The experimental setup is shown in Figure 8. The array is placed at a height of 1.1 m above the ground, and Type 4958 microphones from Brüel & Kjær company are used. The spacing between two adjacent microphones is $\mathrm{\Delta }x=\mathrm{\Delta }y=0.07\hspace{0.33em}\mathrm{m}$ , therefore the upper frequency limit of this array is approximately 2450 Hz. A Cartesian coordinate system is established with the (0, 0) th microphone as the origin. A loudspeaker driven by white noise is placed at (−2.24, 0, 5) m. The DOAs of the loudspeaker source and its ground-reflected image source are about $\left(24.1°,180°\right)$ and $\left(32.1°,224.5°\right)$ , respectively. Brüel & Kjær PULSE Type 3660C data acquisition system is used to simultaneously record the signals from all microphones with a sampling rate of 16384 Hz. When performing FFT, the Hanning weighting is used, and the frequency resolution is set to 1 Hz.

Figure 8.

Experimental layout.

Figure 9(a) shows the spatial spectral function of the proposed method when processing the pressures at 6500 Hz and 7000 Hz. It can be observed that many spurious peaks caused by spatial aliasing are effectively suppressed. Figure 9(b) and 9(c) present the imaging results of the proposed method and the method in Ref. ¹⁷ at 6500 Hz. The method in Ref. ¹⁷ regards the positions where spatial aliasing occurs as sources. The proposed method successfully locates the sources close to their true positions, though there is a weak spurious source within the dynamic range. The experiment draws a similar conclusion to the simulations. In summary, the proposed two-dimensional multi-frequency synchronous grid-free CB method effectively suppresses spatial aliasing and accurately identifies sources whose wavelengths are shorter than twice the minimum microphone spacing.

Figure 9.

Experimental results of the proposed method (multi-frequency synchronous modelling) and the method in Ref. ¹⁷ under spatial aliasing.

5. Conclusion

To overcome high-frequency spatial aliasing problem of the two-dimensional single-frequency grid-free CB method, this paper proposes a two-dimensional multi-frequency synchronous grid-free CB method. The approach formulates a two-dimensional multi-frequency synchronous atomic norm minimization problem, then converts it into a dual problem for solving, and finally extracts DOAs and strengths by using a spatial spectral function related to the dual polynomial. Simulations and experiments demonstrate that when the spatial Nyquist criterion is satisfied, both the proposed method and the existing single-frequency method can achieve accurate source identification, while the proposed method reduces strength quantification error by over 80% under strong noise interference. When the spatial Nyquist criterion is violated, the proposed method effectively suppresses spatial aliasing and accurately identifies sources whose wavelengths are shorter than twice the minimum microphone spacing. The dual problem is currently solved using the CVX toolbox, which suffers from low computational efficiency. Future work will focus on developing a more efficient solver to reduce computational cost, as well as extending the method to non-uniform array geometries.