Iterative algorithm for solving acoustic source characterization problems under block sparsity constraints

Noise Source Identification (NSI) techniques seek to pinpoint the locations, waveforms, and the associated acoustic variables of the noise sources in industrial, transportation, and environmental scenarios.¹ Many array-based NSI techniques such as the Nearfield Acoustic Holography (NAH), beamforming, deconvolution method, etc., have been suggested in the past.^2,3 NAH that requires measurement in the nearfield of sources is a high-resolution reconstruction method.⁴ NAH methods can be implemented on the basis of spatial Fourier transforms,^5,6 plane or spherical wave expansions,^7–9 boundary element methods,^10,11 etc.

The Equivalent Source Method (ESM),^12–14 also known as the method of wave superposition, is a simple method that is not limited to sources of regular geometry. The idea of ESM lies in that an arbitrary sound field can be approximated as the superposition of waves emanated by an array of discrete simple sources. A small standoff distance is generally maintained between the equivalent sources and the real source to prevent singularity problems when reconstructing the sound field on the source surface. A linear inverse problem can be formulated via ESM. However, for the spatially extended, the number of equivalent sources assumed can be far greater than the number of sensors. In this case, the inverse problem becomes an underdetermined one whose solution is not unique. The classical Least-Squares (LS) solution with pseudoinverse gives only minimum-norm solution, which tends to unduly spread errors in the solution space. A newly emerging technique, Compressive Sensing (CS), tackles this nonuniqueness problem by imposing the sparsity constraint. The sparse solution can be obtained by solving a convex optimization problem.^15–18 Numerical packages such as CVX^19,20 are available for solving the convex optimization problems. However, the CS approaches are useful when the sources are discrete and sparse. For spatially extended “block” sources, the CS approaches could break down because the sparsity condition is not fulfilled. In addition, package CVX is very computationally expansive. As a more efficient alternative, the Wideband Holography (WBH), was suggested by Hald²¹ to solve the inverse problem by using a modified method of steepest descent with pruning. Although WBH is computationally more efficient than the CVX, it primarily only deals with point sources. To handle block-sparse sources, beamforming using Total Variation (TV) norm was proposed by Xenaki et al.²² Their method mainly deals with one-dimensional and coherent block sources.

In the present paper, an iterative CS algorithm is proposed to solve block-sparse inverse problems and to reconstruct sound field. Source localization and signal extraction are accomplished in a unified CS framework. The iterative CS algorithm, also referred to as Compressive Newton's (CNT) hereafter, is based on a modified Newton's method^23–25 to calculate the source amplitudes with the aid of a special pruning procedure. The CNT method converges much faster than the WBH method with a pruning method different from the method by Xenaki et al.,²² which can handle two-dimensional and incoherent block sources. In addition, by using the CNT method, the acoustic variables including sound pressure, particle velocity, sound intensity, and sound power, and the radiation pattern can all be calculated, based on the ESM.

This paper is organized as follows. Section II outlines the basic theory of ESM and underdetermined inverse NSI problem. Section III describes the key concepts of the CNT algorithm. Section IV presents calculations of acoustic variables based on the ESM. Numerical and experimental investigations are presented in Secs. V and VI. Section VII summarizes the strengths as well as limitations of the CNT technique.

In this section, the acoustic source characterization problem is formulated as an inverse problem on the basis of the ESM. The simple sources used as equivalent sources are monopoles.

Figure 1 depicts the idea of the ESM by using which the sound field is represented with a large number of virtual monopole sources. The sound pressure emitted by the actual source is sampled in space and time by an array of microphones. The equivalent sources are assumed to be distributed behind the actual source surface with a standoff distance. The acoustic variables on the reconstruction surface can be calculated by the field propagation of these equivalent sources.

Consider M microphones and N equivalent monopole sources that emit spherical waves. The pressure vector p(ω) received at the microphones can be written in the frequency domain as²⁶ 

where $p∈ℂM$⁠, the matrix $A=[a1a2⋯aN]∈ℂM×N$ is the steering matrix and $an=[e−jk⋅r1n/r1ne−jk⋅r2n/r2n⋯e−jk⋅rMn/rMn]T$⁠, $n=1,…,N,$ denotes the steering vector of the nth source, with the superscript “T” denoting matrix transpose. Here, time dependence $ejωt$ is assumed, with $ω$ and t denoting angular frequency and time, and $j=−1$⁠. The distance, $rmn=‖rn−rm‖$⁠, with r_n and r_m being the position vectors of the nth equivalent source and the mth microphone. The source amplitude vector s(ω) = [s₁(ω)… s_N(ω)]^T$∈ℂN$⁠. The vector n(ω) represents the additive noise that is uncorrelated with the source signals.

Given the pressure vector p and the steering matrix A, solving the problem of Eq. (1) for the unknown source amplitude vector s is a linear inverse problem. For spatially extended sources, the number of equivalent sources must be sufficiently large, or equivalently the search grid must be fine enough, to serve as the dictionary and is in general far greater than the number of microphones. Therefore, the problem is typically underdetermined, i.e., $N≫M$⁠. In this case, the solution is generally not unique unless imposing additional constraints. Although a traditional LS solution using pseudoinverse gives a unique minimum-norm LS solution for the underdetermined problem, the error tends to spread evenly in the solution space. Unlike the LS solution, the CS approach imposes a sparsity constraint to limit the cardinality (nonzero entries) of the solutions.^15–17 The CS problem can be cast into the following constrained optimization problem:

where $‖s‖0$ denotes the 0-norm, or the cardinality, of the vector s and A is the steering matrix or the sensing matrix in the CS context. However, the CS problem above is not convex. A convex approximation of the problem that also provides robustness to noise is generally used.

The parameter ε is a threshold that can be selected with reference to the aforementioned LS solution. Numerous methods can be used to solve this convex optimization problem.^27,28 This problem can be solved numerically by a convex optimization package (CVX).^18–20 It is noteworthy that the preceding CS approach works for only point-sparsity. In order for the nonzero point sources to be resolvable, the number of these sources must be less than the number of microphones (M).²² This requirement is a clear violation for problems of spatially extended sources which should be characterized as problems with block-sparsity.

Although CS problems can be solved by using the package CVX, the computational cost is extremely high to prevent itself from many time-critical applications. The present paper proposes a CNT algorithm to solve CS problems with block-sparsity constraints. The proposed algorithm is based on Newton's (NT) method into which a pruning procedure of sparsity constraints is incorporated.

The NT method, also known as the Newton-Raphson method, is a numerical method that was originally developed to find successively the roots (or zeroes) of a real-valued function.^20–22 The NT method can be extended to vector and complex-valued cases. The NSI problem described as a linear system of equations in Eq. (1) can be solved iteratively by using the NT method. Instead of solving the linear system of equations directly, the problem is formulated as an unconstrained optimization problem with the cost function

to be minimized. In Eq. (4), s is source amplitude vector. The preceding expression is of hermitian form and thus has a unique minimum. The stationary point can be determined by taking the complex gradient of F.²⁹ 

Note that the gradient vector

with the residual vector $r=As−p$ and G being the Jacobian matrix whose entries are defined as $Gkj=∂rk/∂xj$⁠.

It follows from Eq. (5) that the optimization problem of Eq. (4) becomes a root-finding problem of Eq. (5). Let s₀ be a vector in the neighborhood of the root vector. Expand the gradient vector in the neighborhood of the root vector by the vector Taylor series

This leads to

where H is the Hessian matrix

With some manipulations, it can be shown that the ijth-entry of H is $Hij=∂/∂sj*(∑k=1m∂rk/∂si*(rk))=∑k=1m((∂rk/∂si*)(∂rk/∂sj*)+(∂2rk/∂sj*∂si*)rk)$⁠. In matrix form

$where Dij=∑k=1m(∂2rk/∂si*∂sj*)rk$⁠. As inspired by Eq. (8), the solution can be updated as

where the superscript “+” denotes the pseudoinverse operation. This is generally referred to as the Newton-Raphson method. By omitting D, which is comprised of the second-order derivatives, the Newton-Raphson method reduces to the Gauss-Newton method

Another well-known steepest decent method can also be regarded as an approximation of the Newton-Raphson method when the Hessian matrix $H$ is omitted.

The preceding formulation is general in that it is applicable to nonlinear system of equations, $f(s)=p$⁠. In a linear case as examined in the paper, $f(s)=As$⁠, it is straightforward to show that

It follows from Eq. (11) that

By the singular value decomposition (SVD) of A

As a result, both the Newton-Raphson method and the Gauss-Newton method lead to the same update equation in the linear case.

where i is the iteration index. The step size μ is used to control the convergence rate.

In the proposed CNT method, the key step to impose the sparsity constraint is to incorporate a pruning procedure into the Newton's iteration. In the following, a binary mask is employed to promote the sparsity constraints:

where i denotes iteration index, n denotes source number, $smax(i)=max(|sn(i)|,n=1,…,N)$ with $sn(i)$ being the nth entry of the source amplitude vector s(i) at the ith iteration, $sn(i)/|sn(i)|$ being the phase term of s_n(i), and D(i) is the pruning threshold selected to be a fraction of $smax(i)$

The choice of the parameter γ enables migration from point sparsity constraint to block sparsity constraint as γ varies from 1 to 0. In light of the pruning action, the magnitude of the equivalent source amplitude will be reset by the binary mask to a maximum level with phase unchanged if its magnitude is greater than the threshold. Otherwise, the source amplitude will be reset to zero. Figure 2 illustrates the pruning procedure. With the pruning procedure, the source magnitudes become increasingly sparse and clustered such that the block characteristics are enhanced. Hence, the source vector is updated by the Newton's method and pruned with the binary mask in each iteration. The stopping criterion is

The iteration stops if the difference between two successive solution norms is less than 0.01. The proposed CNT algorithm is summarized in the following pseudo-codes and the flowchart of Fig. 3.

By using the CNT algorithm, the complex amplitudes of the source signals are obtained. The acoustic variables including sound pressure (p), particle velocity (u), sound intensity (I), sound power (W), and radiation pattern can then be calculated on the basis of the ESM.^30,31

The sound pressure propagated from equivalent sources can be reconstructed as

where r_mv is the distance of the vth equivalent source and the mth field point on the reconstruction surface, i.e., $rmv=|zm−zv|$⁠, with z_m and z_v being the position vectors of the vth equivalent source and the mth field point on the reconstruction surface.

The normal component of particle velocity can also be obtained

where ρ₀ is the air density, $e=(zm−zv)/rmv$ and n is the unit vector normal to the reconstruction surface.

The microphone array model used in the present work is a random array with sensor locations determined in previous research.³¹ The array size is 0.48 m × 0.4 m, and has 30 microphones. The sampling frequency is 16 kHz. The Fast Fourier Transform (FFT) block size is 2048, processed with Hanning window, and averaged with 50% overlap.

For the simulation, consider two point sources located at (0 m, 0 m, 1 m) and (0.8 m, 0.3 m, 1 m), emanating independent speech signals. The microphone array is placed on the plane, z = 0 m. On the equivalent source plane that covers a 1.5 m × 1.5 m square area, 100 (10 $×$ 10) virtual monopole sources are distributed at a distance 1.05 m from the array. The reconstruction plane is parallel to the microphone array plane at z = 1 m. The reconstruction plane is of the same size as the microphone array plane. On the reconstruction plane, $31×31$ grid points are selected at a distance 0.05 m from the equivalent source plane. The signal-to-noise ratio (SNR) of the test data used in the simulation is 40 dB. It accommodates 100× loss of SNR due to the illposedness in the inversion process. Lower SNR would only result in worse quality of the extracted source signal, which necessitates heavier regularization. Comprehensive investigation on this aspect can be found in the Ref. 25.

The pruning parameter γ = 0.95 is assumed in the CNT algorithm because only point sparsity is involved in this case. In addition, the steepest descent (SD)²⁰ and conjugate gradient (CG)^32,33 method with pruning are also used for benchmarking the CNT method. Figure 4 shows the reconstructed pressure magnitude map. For simplicity, all data have been normalized within 0–1 in linear scale. The results have shown that all methods but the CG method located the sources correctly. Furthermore, we not only compare the proposed CNT method with the CVX package, but also compare with the SD, CG, and FOCUSS algorithms.³⁴ The source speech signals have been well separated by using the CNT algorithm, as summarized in Table I. Perceptual Evaluation of Speech Quality (PESQ)³⁵ is employed as an objective test metric to evaluate speech quality. The mean opinion score (MOS) is calculated, ranging from 1 to 5. MOS signifies the difference in speech quality between the clean source signal and the source signal extracted by the separation algorithms.

Although CVX achieves the highest PESQ, it is extremely time-consuming. Three iterative CS approaches are far more computationally efficient than the CVX method. The processing time required by the proposed CNT method which achieves slightly lower MOS than the CVX is only 1% of that required by the CVX.

The baffled circular rigid piston of radius 0.144 m, which has a farfield closed-form solution, is employed as a spatially extended source model. The rigid piston is driven with unit velocity amplitude at the frequency 3800 Hz, which amounts to $ka=10$⁠. The piston is assumed to be located at 1 m away from the microphone array. There are 100 × 100 virtual sources, covering 0.6 m × 0.6 m square area, distributed 1.05 m away from the microphone array. The reconstruction surface has $31×31$ grid points located at a standoff distance 0.05 m away from the equivalent source surface, covering a 0.3 m $×$ 0.3 m area.

The pruning parameter γ = 0.7 is assumed in the CNT algorithm because of the block sparsity involved in this case. The acoustic variables are reconstructed by using the CNT method on the source surface, as shown in Fig. 5, where the circular piston is indicated by a blue dotted curve. The reconstructed velocity coincides well with the piston area. In particular, particle velocity map has achieved a nearly perfect fit. Figure 5(b) further compares the particle velocity profiles reconstructed using the CNT, CVX, and the desired profile of the piston. We compare the reconstructed particle velocity with ground truth. The error is quantified using the metric $‖x̂−x‖/‖x‖(%)$⁠, where $x̂$ is the reconstructed velocity and $x$ is ground truth data. The reconstruction error of particle velocity profiles for the CNT and the CVX methods are 11.06% and 80.18%, respectively. The velocity profile reconstructed by the CNT method fits quite well the desired profile, whereas the velocity profile reconstructed by the CVX yields a hill instead of a terrace. This result demonstrates that the CNT is able to handle problems with block sparsity in contrast to the CVX that only deals with problems with point sparsity.

The example above is a case of “coherent” block source. In the following, another example of “incoherent” block source is examined by using the same piston driven by unit-magnitude but random-phase velocity distribution. Figure 6 shows the surface acoustic variables reconstructed using the CNT algorithm for the incoherent piston source. The reconstructed particle velocity remains to cover the piston area with constant velocity magnitude. This result demonstrates the capability of the proposed CNT algorithm in imaging incoherent block sources.

Figure 7 shows the learning curve of solution norms obtained with three iterative CS methods. Apparently, the CNT method converges at a faster speed than the CG method to the desired reconstruction within pre-specified tolerance of the solution norm increment [Eq. (17)]. Interestingly, the CG settles to a biased solution norm, which explains its inferior performance as compared to the other two methods.

To verify the proposed algorithm, experiments are undertaken in a 5.4 m × 3.5 m × 2 m anechoic room, as shown in Fig. 8. In the experimental arrangement, 30 channels of 1/4 in. PCB^® condenser microphones are utilized to capture the noise signals. NI PXI-1042Q^® is used as the data acquisition system operated under the sample rate 16 kHz. The origin of the coordinate system is selected to be the bottom-left corner of the rectangular random array holder.

In Fig. 8(a), a midrange loudspeaker broadcasting female speech and a hand cutter are used as two example sources. These sources are located at (0 m, 0 m, 1 m) and (0.7 m, 0.5 m, 1 m), respectively. The hand cutter can be regarded as a point source, while the loudspeaker can be regarded as a coherent block source. In total, $31×31$ equivalent sources covering a 1.5 m × 1.5 m square area 1.05 m away from the array are used in the ESM. The reconstruction surface has identical grid points allocated, with a 0.05 m standoff distance from the equivalent source surface.

The pruning parameter γ = 0.7 is assumed in the CNT algorithm because both point sparsity and block sparsity are involved in this case. Figure 9 shows the result of virtual source amplitude obtained using different methods. All methods demonstrate accurate localization capability. In particular, the virtual sources obtained using the CNT method appear to be distributed in a more concentrated way than the other two methods for the hand cutter. However, there is a drastic difference in the processing time (310 s for the CNT method and 30 842 s for the CVX). That is, the CNT method is computationally more efficient than the CVX by approximately a factor of 100. The acoustic variables and radiation pattern calculated based on the ESM are shown in Fig. 10. All imagery, in particular the velocity map, coincides with the actual source locations very clearly. The sound powers produced by the hand cutter and the loudspeaker are estimated to be 1.05 and 0.077 W, respectively.

In Fig. 8(b), all experimental arrangements are similar to those in Sec. VI A. Instead of a loudspeaker which is better characterized as a coherent source, however, an electrical fan is utilized to serve as an incoherent block source. The sources are located at (0 m, 0 m, 1 m) and (0.8 m, 0.5 m, 1 m), respectively. As before, the hand cutter can be regarded as a monopole source.

Figure 11 shows the result of virtual source amplitudes calculated using three CS methods. The CNT method fits the outline of the fan, whereas the CG method yields a slightly larger area. Note that the CVX method gives a concentrated hotspot for the fan because it deals with only point sparsity. The acoustic variables calculated are shown in Fig. 12. All imagery, especially the particle velocity map, shows the source locations clearly. The choice of the pruning parameter (γ) has a profound impact on the reconstructed velocity distribution. Solutions progressively migrating from point sparsity to block sparsity is promoted as γ decreases from 1 to 0. To better manifest this important point, we calculate the velocity maps with three different pruning settings: without pruning (γ = 1.0), with pruning (γ = 0.9), and with pruning (γ = 0.7). Without pruning, only the hand cutter is identified as a point source, while the reconstructed velocity of the fan is blurred. With light pruning (γ = 0.9), both sources are identified as point sources. Last, with heavy pruning (γ = 0.7), only the fan source is identified as a continuous source, while the hand cutter is nearly missed because of the over-dimmed image. With light pruning (γ = 0.75), both sources are identified correctly. Figure 13 shows acoustic variables calculated by the CVX method. The results indicate that the CVX method fails to “outline” the fan which represents an incoherent block source. The processing times are 327 s for the CNT method and 30 147 s for the CVX method. The sound powers produced by the hand cutter and the fan are estimated to be 0.094 and 1.023 W, respectively.

A noise source characterization technique based on the iterative CS method has been developed. The proposed CNT method incorporates point and block sparsity constraints into the iteration process by using a binary mask, enabling the identification of not only point sources but also two-dimensional incoherent block sources.

Simulated results have shown that the source signals extracted using the CNT method attain good voice quality assessed by PESQ, with little processing time. For the rigid piston example, the CNT method provides very good reconstruction of the sound field in agreement with the theoretical solution.

The experimental results have demonstrated that the proposed CNT method is capable of reconstructing the sound field produced by an incoherent block source. The CNT method requires little computational cost compared to the numerical package CVX in an order of 100. In the post-processing phase, acoustic variables including particle velocity can be calculated on the basis of ESM. In particular, particle velocity map coincides well with the outline of block sources regardless of the coherence of source surface velocity distribution.

The work was supported by the Ministry of Science and Technology (MOST) in Taiwan, Republic of China, under the project number 105-2221-E-007-030-MY3. Thanks also go to Mr. E. Chin-Pu Tsai for his enthusiastic support of the Telecom Electroacoustics Audio laboratory (TEA lab).