A novel Fourier-based deconvolution algorithm with improved efficiency and convergence

Abstract

Various deconvolution algorithms for acoustic source are developed to improve spatial resolution and suppress sidelobe of the conventional beamforming. To improve the computational efficiency and solution convergence of deconvolution, this paper proposes a Fourier-based improved fast iterative shrinkage thresholding algorithm. Simulations and experiments show that Fourier-based improved fast iterative shrinkage thresholding algorithm can achieve excellent acoustic identification performance, with high computational efficiency and good convergence. For Fourier-based improved fast iterative shrinkage thresholding algorithm, the larger the weight coefficient, the narrower the mainlobe width, and the better the convergence, but the spurious source also increases. The recommended weight coefficient for the array described herein is 3. In addition, like other Fourier-based deconvolution algorithms, Fourier-based improved fast iterative shrinkage thresholding algorithm using irregular focus grid can obtain better acoustic source identification performance than using the conventional regular focus grid.

Introduction

Array-based beamforming has become an indispensable acoustic source identification technology in some industries such as aviation aircraft,^1–4 rail train, ⁵ wind turbine,^6–9 mining industry, ¹⁰ and automobile^11–13 over decades. Conventional beamforming (CB) has the advantages of fast measurement speed, high computational efficiency, and suitable for medium-to-long distance measurement. However, due to the insufficient spatial resolution and lack of spurious sources attenuation, the acoustic source identification accuracy is low.^11,14–17 To overcome the above disadvantages, many deconvolution algorithms are developed, such as deconvolution approach for the mapping of acoustic sources (DAMAS),^18,19 non-negative least square (NNLS), ²⁰ Richardson–Lucy (RL), ²⁰ etc. The core idea of these algorithms is to obtain the real information of the acoustic source through the deconvolution operation, because the output of CB can be approximated as the spatial convolution of the sound source distribution and the point spread function (PSF). However, due to disadvantages like time consuming and slow convergence, the practical application of the above algorithms is somewhat limited. To improve the computational efficiency and convergence of these algorithms, the spatial convolution is transformed into the product in wave-number domain based on fast Fourier transform (FFT). The corresponding algorithms are named DAMAS2, ²¹ FFT-NNLS, ²⁰ FFT-RL, ²⁰ etc., respectively. All the above Fourier-based deconvolution algorithms are built on the assumption of a shift-invariant PSF. It results in poor identification performance of the acoustic source away from the center of the calculation plane. Xenaki et al. ²² proposed a specific coordinate transformation that suppresses the spatial shift-variant of PSF to expand the effective angle of the source identification for the Fourier-based deconvolution algorithm. Ehrenfried and Koop ²⁰ and Chu and Yang ²³ compared the computational efficiency and acoustic source identification accuracy of several Fourier-based deconvolution algorithms. The results show that FFT-NNLS has excellent comprehensive performance.

Recently, Lylloff et al. ²⁴ proposed a Fourier-based fast iterative shrinkage thresholding algorithm (FFT-FISTA). Compared with FFT-NNLS, FFT-FISTA enjoys higher computational efficiency and faster convergence. Like other Fourier-based deconvolution algorithms, FFT-FISTA with irregular focus point grid can overcome the deficiency of failing to accurately identify the sound source far from the center of the calculation plane. ²⁵ FFT-FISTA is derived from FISTA and FISTA is originally proposed to solve the linear inverse problem in image processing. ²⁶ Bhotto et al. ²⁷ introduced a positive definite weight matrix in the gradient function minimization of FISTA and proposed improved FISTA (IFISTA) to enhance the convergence and the image reconstruction accuracy of FISTA. Inspired by Lylloff et al., ²⁴ Chu et al., ²⁵ Beck and Teboulle, ²⁶ and Bhotto et al., ²⁷ to further improve the computational efficiency and the convergence of FFT-FISTA, a Fourier-based IFISTA (FFT-IFISTA) for acoustic source identification is proposed in this paper.

In this paper, a FFT-IFISTA is proposed to enhance the calculating efficiency and convergence of deconvolution. In addition, to achieve a better identification performance, the irregular focus grid is applied. Results of the simulation and experiment indicate that the proposed FFT-IFISTA can acquire a better acoustic identification performance with higher calculating efficiency and better convergence than other methods such as FFT-NNLS and FFT-FISTA.

Methods

Figure 1 shows the layout of beamforming measurement. Vector $\mathbf{r}$ ndicates the position of a focus grid point. Vector ${\mathbf{r}}_m$ indicates the position of the $m\text{th}$ microphone, where $m=1,2,3,\dots ,M$ is the serial number of the microphone. $N_r$ and $N_c$ are the row and column numbers of grid points, respectively. $N_{\mathit{rm}}$ and $N_{\mathit{cm}}$ represent the row and column indexes of center grid point. The output of CB based on cross-spectral imaging function ²³ is $$\{\begin{array}{c}b\left(\mathbf{r}\right)={\sum }_{{\mathbf{r}}^{\prime }}q\left({\mathbf{r}}^{\prime }\right)p\mathit{sf}\left(\mathbf{r}|{\mathbf{r}}^{\prime }\right)\\ \mathit{psf}\left(\mathbf{r}|{\mathbf{r}}^{\prime }\right)=\left|{\mathbf{r}}^{\prime }\right|\frac{{\mathbf{v}}^{\text{T}}\left(\mathbf{r}\right)\left[{\mathbf{v}}^{*}\left({\mathbf{r}}^{\prime }\right){\mathbf{v}}^{\text{T}}\left({\mathbf{r}}^{\prime }\right)\right]{\mathbf{v}}^{*}\left(\mathbf{r}\right)}{\sqrt{{\mathbf{w}}^{\text{T}}\left(\mathbf{r}\right)\mathbf{1}{\mathbf{w}}^{*}\left(\mathbf{r}\right)}}\end{array}$$ (1)where ${\mathbf{r}}^{\text{′}}$ is the position of the assumed point source; $q\left({\mathbf{r}}^{\text{′}}\right)$ is its sound pressure contribution at the center of the array; $\mathit{psf}\left(\mathbf{r}|{\mathbf{r}}^{\text{′}}\right)$ is the PSF, which expresses the beamforming output in the $\mathbf{r}$ position of the unit-strength point source in the ${\mathbf{r}}^{\text{′}}$ position; $\mathbf{1}$ is a $M\times M$ matrix with all elements equal to 1; $\mathbf{v}\left(\mathbf{r}\right)={\left[v_1\left(\mathbf{r}\right),v_2\left(\mathbf{r}\right),\dots v_m\left(\mathbf{r}\right),\dots v_M\left(\mathbf{r}\right)\right]}^{\text{T}}$ is the steering vector; and $\mathbf{w}\left(\mathbf{r}\right)\equiv {\left[{\left|v_1\left(\mathbf{r}\right)\right|}^2,{\left|v_2\left(\mathbf{r}\right)\right|}^2,\dots {\left|v_m\left(\mathbf{r}\right)\right|}^2\dots {\left|v_M\left(\mathbf{r}\right)\right|}^2\right]}^{\text{T}}$ . The superscripts “T” and “*” represent the transposition and the conjugation respectively. $v_m\left(\mathbf{r}\right)$ is defined by $$v_m\left(\mathbf{r}\right)=\frac{\exp \left(-\text{i}k\left|\mathbf{r}-{\mathbf{r}}_m\right|\right)}{\left|\mathbf{r}-{\mathbf{r}}_m\right|}$$ (2)where $k=2\pi f/c$ is the wave number, $f$ is the frequency, and $c$ is the sound speed, $\text{i}=\sqrt{-1}$ .

Figure 1.

Layout of beamforming measurement.

Furthermore, the error function $\text{φ}$ is constructed between the PSF, the sound source distribution, and the output of CB, as shown in equation (3) $$\text{φ}=\frac{1}{2}{||\mathbf{Aq}-\mathbf{b}||}_2^2$$ (3)where ${\Vert \Vert }_2$ represents the ℓ₂ norm, $\mathbf{A}=\left[\mathit{psf}\left(\mathbf{r}|{\mathbf{r}}^{\text{′}}\right)\right]$ is the known $N\times N$ PSF matrix, $N=N_r\times N_c$ is the total number of the focus grid points, $\mathbf{q}=\left[q({\mathbf{r}}^{\text{′}})\right]$ is the unknown $N\times 1$ source distribution vector and each element is nonnegative, and $\mathbf{b}=\left[b\left(\mathbf{r}\right)\right]$ is the known $N\times 1$ output vector of CB. In the following, the real parts of elements in $\mathbf{A}$ and $\mathbf{b}$ are taken.

To obtain $\mathbf{q}$ , IFISTA utilizes a fast projected gradient descent algorithm ²² as shown in equation (4) $${\mathbf{q}}_{l+1}=\rho +\left({\mathbf{q}}_l-\frac{1}{L}\nabla \text{φ}\left({\mathbf{q}}_l\right)\right)$$ (4)where ${\rho }_{+}$ is the Euclidean projection of $\mathbf{q}$ onto the nonnegative quadrant, ${\mathbf{q}}_l$ and ${\mathbf{q}}_{l+1}$ represent the source pressure contribution after $l\text{th}$ and $(l+1)\text{th}$ iterations, $L$ is the Lipschitz constant and equals to the largest eigenvalue of ${\mathbf{A}}^{\mathrm{T}}\mathbf{A}$ , ²⁸ and $\nabla \text{φ}$ is the Lipschitz continuous gradient and its expression is $$\nabla \text{φ}\left({\mathbf{q}}_l\right)={\mathbf{W}}_n{\mathbf{A}}^{\text{T}}\left(\mathbf{A}{\mathbf{q}}_l-\mathbf{b}\right)$$ (5)where ${\mathbf{W}}_n$ is the positive definite weight matrix and the expression is given by $${\mathbf{W}}_n={\displaystyle \sum }_{j=1}^n\left(\begin{array}{c}n\\ j\end{array}\right){\left(-\frac{1}{L}\right)}^{j-1}{\left({\mathbf{A}}^{\text{T}}\mathbf{A}\right)}^{j-1}$$ (6)

Due to the introduction of positive definite weight matrix ${\mathbf{W}}_n$ , $n$ is weight coefficient. IFISTA cannot conduct the Fourier transform based on the zero-boundary condition like FFT-NNLS and FFT-FISTA to improve the computational efficiency. In this paper, Fourier transform based on the periodic boundary condition^29,30 is applied in the IFISTA. The PSF matrix $\mathbf{A}$ is a block circulant with circulant blocks (BCCB) normal matrix under the shift-invariant PSF assumption, that is ${\mathbf{A}}^{\text{H}}\mathbf{A}=\mathbf{A}{\mathbf{A}}^{\text{H}}$ , where the superscript “H” represents the transposition conjugation. The BCCB matrix $\mathbf{A}$ has a unitary spectral decomposition, whose expression is given by $$\mathbf{A}={\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}$$ (7)where $\mathbf{F}$ is the two-dimensional unitary discrete Fourier transform matrix. ${\mathbf{F}}^{\text{H}}={\mathbf{F}}^{-1}$ , ${\mathbf{F}}^{-1}$ is the inverse matrix of $\mathbf{F}$ , $\mathbf{\Lambda }$ is the eigenvalue matrix of PSF matrix $\mathbf{A}$ .

For PSF matrix $\mathbf{A}$ is real, that is $${\mathbf{A}}^{\text{T}}\mathbf{A}={\mathbf{A}}^{\text{H}}\mathbf{A}={\left({\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}\right)}^{\text{H}}\left({\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}\right)={\mathbf{F}}^{\text{H}}{\mathbf{\Lambda }}^{\text{H}}\mathbf{F}{\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}$$ (8)where $\mathbf{F}{\mathbf{F}}^{\text{H}}=\mathbf{I}$ , $\mathbf{I}$ is the unit matrix and equation (8) can be expressed as $${\mathbf{A}}^{\text{T}}\mathbf{A}={\mathbf{F}}^{\text{H}}{\mathbf{\Lambda }}^{\text{H}}\mathbf{\Lambda }\mathbf{F}={\mathbf{F}}^{\text{H}}\mathbf{\Delta F}$$ (9)where $\mathbf{\Delta }$ is the eigenvalue matrix of ${\mathbf{A}}^{\text{T}}\mathbf{A}$ and equals to ${\mathbf{\Lambda }}^{\text{H}}\mathbf{\Lambda }$ .

Matrix $\mathbf{F}$ has a very convenient property. The following equivalence relationship is established $$\left\{\begin{array}{l}\sqrt{N}F{\mathbf{q}}_l\leftrightarrow F\left({\mathbf{Q}}_l\right)\\ \frac{1}{\sqrt{N}}{\mathbf{F}}^{\text{H}}{\mathbf{q}}_l\leftrightarrow F^{-1}\left({\mathbf{Q}}_l\right)\\ {\mathbf{F}}^{\text{H}}F{\mathbf{q}}_l\leftrightarrow F^{-1}\left(F\left({\mathbf{Q}}_l\right)\right)\end{array}\right.$$ (10) $\text{“}F\text{”}$ and $\text{“}{F}^{-1}\text{” }$ are Fourier transform and its inverse transform, respectively, $\text{“}\leftrightarrow \text{”}$ represents equivalence relationship, and ${\mathbf{Q}}_l$ is a $N_r\times N_c$ matrix reshaped from the vector ${\mathbf{q}}_l$ . From equation (7) $$\mathbf{FA}\left(:,1\right)=\mathbf{\Lambda }\mathbf{F}\left(:,1\right)$$ (11) $\text{“}:\text{”}$ indicates that all rows are fetched, for $\mathbf{F}\left(:,1\right)={\left[1,1,\cdots ,1\right]}^{\text{T}}/\sqrt{N}$ , that is $$\mathbf{\Lambda }\mathbf{F}\left(:,1\right)=\frac{1}{\sqrt{N}}\widehat{\mathbf{\Lambda }}$$ (12) $\widehat{\mathbf{\Lambda }}$ is a column vector reshaped from $\mathbf{A}$ . Combining equations (10), (11), and (12) $$\widehat{\mathbf{\Lambda }}=\sqrt{N}\mathbf{F}\mathbf{A}\left(:,1\right)\leftrightarrow F\left(\tilde{\mathbf{A}}\right)$$ (13) $\tilde{\mathbf{A}}$ is a $N_r\times N_c$ matrix reshaped from the first column of the matrix $\mathbf{A}$ , which can be obtained by moving up and left the elements in the PSF matrix of the acoustic source at the central focus point $N_{\mathit{rm}}-1$ , $N_{\mathit{cm}}-1$ times, respectively. Combining equations (10) to (13) $$\mathbf{A}{\mathbf{q}}_l={\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}{\mathbf{q}}_l\leftrightarrow F^{-1}\left(F\left({\mathbf{Q}}_l\right)○F\left(\tilde{\mathbf{A}}\right)\right)$$ (14) $\text{“}○\text{”}$ represents the Hadamard product operation.

Substituting equation (9) into equation (6) $${\mathbf{W}}_n=\underset{j=1}{\overset{n}{}}\left(\begin{array}{c}n\\ j\end{array}\right){(-\frac{1}{L})}^{j-1}{\left({\mathbf{F}}^{\text{H}}\mathbf{\Delta F}\right)}^{j-1}={\mathbf{F}}^{\text{H}}\left(\underset{j=1}{\overset{n}{}}\left(\begin{array}{c}n\\ j\end{array}\right){(-\frac{1}{L})}^{j-1}{\mathbf{\Delta }}^{j-1}\right)\mathbf{F}={\mathbf{F}}^{\text{H}}\mathbf{\Theta }\mathbf{F}$$ (15)where $\mathbf{\Theta }$ is the eigenvalue matrix of ${\mathbf{W}}_n$ $$\mathbf{\Theta }=\underset{j=1}{\overset{n}{}}\left(\begin{array}{c}n\\ j\end{array}\right){(-\frac{1}{L})}^{j-1}{\mathbf{\Delta }}^{j-1}$$ (16)

Substituting equations (7) and (15) into equation (5) $${\mathbf{W}}_n{\mathbf{A}}^{\text{T}}\left(\mathbf{A}{\mathbf{q}}_l-\mathbf{b}\right)={\mathbf{F}}^{\text{H}}\mathbf{\Theta }\mathbf{F}{\mathbf{F}}^{\text{H}}{\mathbf{\Lambda }}^{\text{H}}\mathbf{F}\left({\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}{\mathbf{q}}_l-\mathbf{b}\right)={\mathbf{F}}^{\text{H}}\left(\mathbf{\Theta }{\mathbf{\Lambda }}^{\text{H}}\right)\mathbf{F}\left({\mathbf{F}}^{\text{H}}\mathbf{\Lambda }\mathbf{F}{\mathbf{q}}_l-\mathbf{b}\right)$$ (17)

Under the periodic boundary condition, the Fourier transform of equation (17) can be expressed as $${\mathbf{W}}_n{\mathbf{A}}^{\text{T}}\left(\mathbf{A}{\mathbf{q}}_l-\mathbf{b}\right)\leftrightarrow F^{-1}[\mathbf{R}○F[F^{-1}\left[F\left({\mathbf{Q}}_l\right)○F\left(\tilde{\mathbf{A}}\right)\right]-\mathbf{B}]]$$ (18)where $\mathbf{B}$ is a $N_r\times N_c$ matrix reshaped from the vector $\mathbf{b}$ , and $$\mathbf{R}=\left({\displaystyle \sum }_{j=1}^n\left(\begin{array}{c}n\\ j\end{array}\right){\left(-\frac{1}{L}\right)}^{j-1}{\left({\mathbf{\Lambda }}^{\text{H}}\mathbf{\Lambda }\right)}^{j-1}\right){\mathbf{\Lambda }}^{\text{H}}\leftrightarrow \mathbf{R}={\displaystyle \sum }_{j=1}^n\left(\begin{array}{c}n\\ j\end{array}\right){\left(-\frac{1}{L}\right)}^{j-1}F{\left(\tilde{\mathbf{A}}\right)}^{○(j-1)}○{\left({\left(F\left(\tilde{\mathbf{A}}\right)\right)}^{\text{H}}\right)}^{○j}$$ (19) $\text{“}{(\cdot )}^{○(\cdot )}\text{”}$ represents the Hadamard exponentiation.

Hence, FFT-IFISTA converts equation (3) into a Fourier-based minimization equation $$\text{φ}=\frac{1}{2}\Vert F^{-1}\left[F\left({\mathbf{Q}}_l\right)○F\left(\tilde{\mathbf{A}}\right)\right]-\mathbf{B}{]\Vert }_{\text{Fro}}^2$$ (20)where ${\Vert \Vert }_{\text{Fro}}$ denotes the Frobenius norm. The iteration search in the negative gradient direction is applied to obtain $\mathbf{Q}$ . Set ${\mathbf{Q}}_0=\mathbf{0}$ , ${\mathbf{Y}}_1={\mathbf{Q}}_0$ , $t_1=1$ . Before the stopping criterion is reached, the $l\text{th }$ iteration is as follows

${\mathbf{Q}}_l={\rho }_{+}\left({\mathbf{Y}}_l-\frac{1}{L}{F}^{-1}\left[\mathbf{R}○F\left[F^{-1}\left[F\left({\mathbf{Y}}_l\right)○F\left(\tilde{\mathbf{A}}\right)\right]-\mathbf{B}\right]\right]\right)$

$t_{l+1}=\frac{1}{2}\left(1+\sqrt{1+4t_l^2}\right)$

${\mathbf{Y}}_{l+1}={\mathbf{Q}}_l+\frac{t_l-1}{t_{l+1}}\left({\mathbf{Q}}_l-{\mathbf{Q}}_{l-1}\right)$

Results and discussions

Simulations

As shown in Figure 1, the 36-channel Brüel & Kjær sector microphone array with 0.65 m diameter is used in the simulation. Two monopole sources are located at (−0.2, 0.2, 1) m and (0.2, 0.2, 1) m, with a sound pressure contribution of 100 dB. As shown in Figure 2(a), the calculation plane is $1\text{m}\times 1\text{ m}$ with $51\times 51$ focus grid points whose grid space is 0.02 m.

Figure 2.

Focus grid. (a) Regular mesh planar and (b) irregular mesh planar.

Figure 3(a) and (b) shows the results of CB at 2000 and 6000 Hz, respectively. The outputs are normalized based on the maximum value and the display dynamic range of the image is 15 dB. The cross “+” marks the location of the sound source. The results show that the acoustic sources can be localized. However, acoustic sources are fused due to the low spatial resolution of CB at 2000 Hz and lots of spurious sources with high amplitude appear at 6000 Hz.

Figure 3.

CB results at different frequencies. (a) 2000 Hz and (b) 6000 Hz.

Figure 4 shows the 1000 iterations imaging results of FFT-NNLS, FFT-FISTA, and FFT-IFISTA with weight coefficient from 1 to 6 at 2000 Hz. Compared to Figure 3(a), each deconvolution algorithm can effectively attenuate the sidelobe, narrow the mainlobe, and separate the fused acoustic sources. The identified position of the acoustic source is consistent with the true position. However, the above simulations are all based on the regular focus grid as shown in Figure 2(a) whose PSF shift-variant is significant, which deteriorates the identification performance of these Fourier-based deconvolution algorithms. That is, the PSF used in deconvolution is different from the theoretical one at the acoustic source position, which resulted in the shapes of the acoustic sources being irregular and spurious sources appearing near the mainlobe. In short, the identification performance is not good enough.

Figure 4.

Comparison of different deconvolution algorithms with regular focus grid at 2000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; (d) FFT-IFISTA, n = 2; (e) FFT-IFISTA, n = 3; (f) FFT-IFISTA, n = 4; (g) FFT-IFISTA, n = 5; and (h) FFT-IFISTA, n = 6.

To suppress the shift-variant of PSF, an irregular focus grid introduced in Xenaki et al. ²² and Chu et al. ²⁵ as shown in Figure 2(b) is applied. For convenience comparison, the imaging region size of the irregular focus grid is taken the same as that of the regular focus grid.

Figure 5 shows the identification results of the deconvolution algorithm with the irregular focus grid. The frequency is 2000 Hz and the number of iterations is 1000. Compared with the identification results of the regular focus grid, the deconvolution algorithm with irregular focus grid can effectively narrow the mainlobe, improve the spatial resolution, eliminate the spurious source, and round the shape of acoustic source. Figure 5 shows that FFT-NNLS has the widest mainlobe. The mainlobe width of FFT-IFISTA with the value of weight coefficient 1 is similar to that of FFT-FISTA. As the weight coefficient increases, the mainlobe width of FFT-IFISTA gradually decreases. It means that the increment of weight coefficient can make the identification performance of FFT-IFISTA better than the other two algorithms. However, for the array used in this paper, when the weight coefficient increment exceeds 4, the decrease of the mainlobe width is not obvious.

Figure 5.

Comparison of different deconvolution algorithms with irregular focus grid at 2000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; (d) FFT-IFISTA, n = 2; (e) FFT-IFISTA, n = 3; (f) FFT-IFISTA, n = 4; (g) FFT-IFISTA, n = 5; and (h) FFT-IFISTA, n = 6.

Figure 6 shows the simulation acoustic source mapping of deconvolution algorithm with the irregular focus grid at 6000 Hz, and the number of iterations is 1000. Compared with Figure 3(b), the deconvolution algorithms can effectively attenuate spurious sources of CB. After 1000 iterations, a few spurious sources appear in the identification results of FFT-NNLS and FFT-FISTA. While FFT-IFISTA with the weight coefficient 1 and 2 can effectively eliminate the spurious source and achieve a better acoustic source identification. When the weight coefficient is 3, spurious sources appear in Figure 6(e). Figure 6(e) to (h) shows that the larger the weight coefficient is, the more spurious sources appear, which makes the acoustic source identification of FFT-IFISTA worse. Therefore, it is recommended to use a smaller weight coefficient to avoid spurious sources at high frequency. According to the identification results at different frequencies, for the array used in this paper, the recommended weight coefficient should not be larger than 3.

Figure 6.

Comparison of different deconvolution algorithms with irregular focus grid at 6000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; (d) FFT-IFISTA, n = 2; (e) FFT-IFISTA, n = 3; (f) FFT-IFISTA, n = 4; (g) FFT-IFISTA, n = 5; and (h) FFT-IFISTA, n = 6.

Figures 7 and 8 show the above-mentioned deconvolution algorithm results with the signal-to-noise ratio (SNR) of 10 dB at 2000 and 6000 Hz, respectively, and the number of iterations is 1000. Compared to Figure 5, the acoustic source in Figure 7 is deformed by the influence of noise, but still it can be accurately identified by these deconvolution algorithms. Comparing Figure 8(a) and (b) with Figure 6(a) and (b), spurious sources of FFT-NNLS and FFT-FISTA are increased. A similar phenomenon of FFT-IFISTA occurs in Figure 8(f) to (h) when the weight coefficient is greater than 3. However, compared with Figure 6(c) to (e), no more spurious sources appear in Figure 8(c) to (e), which indicates that FFT-IFISTA can effectively suppress noise interference when the weight coefficient is not greater than 3. To sum up, FFT-IFISTA with smaller weight coefficient is more robust to noise than FFT-NNLS and FFT-FISTA.

Figure 7.

Comparison of different deconvolution algorithms with irregular focus grid and 10 dB SNR at 2000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; (d) FFT-IFISTA, n = 2; (e) FFT-IFISTA, n = 3; (f) FFT-IFISTA, n = 4; (g) FFT-IFISTA, n = 5; and (h) FFT-IFISTA, n = 6.

Figure 8.

Comparison of different deconvolution algorithms with irregular focus grid and 10 dB SNR at 6000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; (d) FFT-IFISTA, n = 2; (e) FFT-IFISTA, n = 3; (f) FFT-IFISTA, n = 4; (g) FFT-IFISTA, n = 5; and (h) FFT-IFISTA, n = 6.

To analyze the causes of sidelobes in FFT-IFISTA when the weight coefficient is large, a fixed spurious sound source is selected in the figure and the sidelobe level comparison of FFT-IFISTA with the weight coefficients from 3 to 6 at 6000 Hz is given in Figure 9. The display dynamic range of the image in this figure is 40 dB. Comparing Figure 9(a) with Figure 6(e), the selected spurious source in Figure 9(a) is −17.11 dB, which does not appear in Figure 6(e) of 15 dB display dynamic range. With the increase of the weight coefficient, the intensity of the selected spurious source in Figure 9(a) to (d) increases from −17.11 to −8.48 dB gradually. The reason is that due to the increase of the weight coefficient n, the gradient $\nabla \text{φ}$ in equation (5) is increased, which improves the convergence of FFT-IFISTA and increases the intensity of the mainlobe. What is more, the intensity of the spurious sound source is also increased. Within a certain dynamic range, the number of spurious sound sources increases.

Figure 9.

Sidelobe level comparison of FFT-IFISTA with different weight coefficient at 6000 Hz. (a) FFT-IFISTA, n = 3; (b) FFT-IFISTA, n = 4; (c) FFT-IFISTA, n = 5; and (d) FFT-IFISTA, n = 6.

Furthermore, the convergence and computational efficiency of the above deconvolution algorithms are compared. The convergence can be quantified by the standard deviation function between the theoretical acoustic source contribution and the reconstructed one obtained by iteration ²⁰ ; the standard deviation function is $$\text{σ}=\frac{1}{\sqrt{N}}{\Vert \mathbf{q}}_l-{\mathbf{q}}_{e}2\Vert$$ (21)where ${\mathbf{q}}_l$ is the reconstructed source contribution after $l\text{th}$ iterations and ${\mathbf{q}}_e$ is the exact one.

Figure 10 shows the convergence performance of FFT-NNLS, FFT-FISTA, and FFT-IFISTA at 2000 Hz. The black dash-dotted line and the blue dotted line represent the convergence curves of FFT-NNLS and FFT-FISTA, respectively. Pink dash-dotted line, black solid line and red solid line, pink dash-dotted line, black solid line, red dotted line, blue dash-dotted line, and pink dotted line represent the convergence curves of FFT-IFISTA with the weight coefficients 1, 3, and 5, respectively. In Figure 10(a), at the beginning of the iteration, the descent speed of FFT-NNLS standard deviation is the fastest and tends to be stable after 100 iterations, while the curves of FFT-IFISTA and FFT-FISTA are still decreasing. Comparing the curves of the FFT-FISTA and FFT-IFISTA with different weight coefficients, it can be found that the descent speed of FFT-IFISTA with the weight coefficient 1 is slightly faster than that of FFT-FISTA, and the larger the weight coefficient is, the faster the standard deviation of FFT-IFISTA decreases. By contrast of the stable standard deviation of FFT-NNLS, FFT-FISTA, and FFT-IFISTA, it is found that the larger the weight coefficient of FFT-IFISTA is, the less the iterations number of stable standard deviation requires, and the standard deviation is smallest. From the convergence comparison of three algorithms, conclusion can be drawn that the convergence of FFT-IFISTA is the best, followed by FFT-FISTA, and the last is FFT-NNLS. Furthermore, increasing the weight coefficient appropriately can improve the convergence of FFT-IFISTA. Figure 10(b) describes the mainlobe source power cross-sectional plot of the deconvolution map after 1000 iterations. It shows that the source power estimation of FFT-NNLS is the smallest, followed by FFT-FISTA. The larger the weight coefficient, the more accurate FFT-IFISTA source strength estimate.

Figure 10.

Convergence comparison of three deconvolution algorithms. (a) The relation between standard deviation and iterative times and (b) cross-section of maps at y = 0.2 after 1000 iterations. FFT: fast Fourier transform; FFT-FISTA: Fourier-based fast iterative shrinkage thresholding algorithm; FFT-IFISTA: Fourier-based improved fast iterative shrinkage thresholding algorithm; NNLS: non-negative least square.

To compare the computational efficiency of FFT-NNLS, FFT-FISTA, and FFT-IFISTA, time-consuming curves of three algorithms with different iterations are given in Figure 11. Figure 11(a) describes the relation between time consuming and iterative times. It is worth mentioning that, for FFT-IFISTA, the value of the weight coefficient does not affect its computational efficiency, so the figure only gives the result of the weight coefficient 6. Within five iterations, computational efficiency of FFT-IFISTA is slightly lower than that of FFT-FISTA, after that, its computational efficiency is higher than that of the others. At 5000 iterations, the consuming time of FFT-IFISTA is 3.11 s, FFT-FISTA is 8.35 s, and FFT-NNLS is 12.71 s. Figure 11(b) describes the relation between standard deviation and time. As can be seen from Figure 10(a), since the standard deviation of FFT-NNLS converges to a large value and remains stable, Figure 11(b) does not include the relationship between FFT-NNLS standard deviation and time. It shows that FFT-IFISTA takes less computational time than FFT-FISTA to achieve the same standard deviation and the larger the weight coefficient, the less time it takes. From the relation between standard deviation and computing time, it more intuitively shows that FFT-IFISTA converges faster than FFT-FISTA.

Figure 11.

Computational efficiency comparison of three deconvolution algorithms. (a) The relation between time consuming and iterative times and (b) the relation between standard deviation and time. FFT: fast Fourier transform; FFT-FISTA: Fourier-based fast iterative shrinkage thresholding algorithm; FFT-IFISTA: Fourier-based improved fast iterative shrinkage thresholding algorithm; NNLS: non-negative least square.

In summary, FFT-IFISTA has the highest computation efficiency and it takes less time to achieve a satisfactory acoustic source identification performance than FFT-NNLS and FFT-FISTA.

Experiments

Figure 12 shows the experimental configuration. The 36-channel Brüel & Kjær sector microphone array with 0.65 m diameter is used in the experiment and the distance between array and loudspeaker source plane is 1 m. The loudspeakers are located near (−0.2, 0.2, 1) m and (0.2, 0.2, 1) m.

Figure 12.

Experimental configuration.

Acoustic source mapping of CB at 2000 and 6000 Hz is shown in Figure 13(a) and (b), respectively. Source mainlobes fuse together due to the poor spatial resolution at 2000 Hz, and many spurious sources appear in the mapping due to the high sidelobe level at 6000 Hz. The experimental results are consistent with the simulation results.

Figure 13.

CB experimental results. (a) 2000 Hz and (b) 6000 Hz.

Figures 14 and 15 show the acoustics source mapping of three deconvolution algorithms using an irregular focus grid at 2000 and 6000 Hz, respectively. The number of iterations is 1000. According to the recommendation, in the following experimental results of FFT-IFISTA, the weight coefficients 1 and 3 are used.

Figure 14.

Comparison of different deconvolution algorithms with irregular focus grid at 2000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; and (d) FFT-IFISTA, n = 3.

Figure 15.

Comparison of different deconvolution algorithms with irregular focus grid at 6000 Hz. (a) FFT-NNLS; (b) FFT-FISTA; (c) FFT-IFISTA, n = 1; and (d) FFT-IFISTA, n = 3.

Comparing Figures 14 and 15 to Figure 12, three deconvolution algorithms can significantly narrow the mainlobe, improve the spatial resolution, and attenuate the spurious source. From the comparison of all mappings in Figure 14, it can be found that the mainlobe of FFT-IFISTA with the weight coefficient 3 is narrowest, followed by FFT-IFISTA with the weight coefficient 1, next is FFT-FISTA, and the last is FFT-NNLS. From the comparison of all mappings in Figure 15, for FFT-IFISTA, when the weight coefficient is taken as 1, no spurious sources appear in the identification results. When 3 is taken, a few spurious sources appear, resulting in a slightly poorer acoustic source identification. However, it is still better than that of the other two algorithms.

To sum up, when weight coefficient does not exceed 3, the identification performance of FFT-IFISTA is better than FFT-NNLS and FFT-FISTA, which is consistent with the simulation results.

Conclusions

This paper proposes FFT-IFISTA-based deconvolution algorithms for acoustic source identification. Like other deconvolution algorithms such as FFT-NNLS and FFT-FISTA, FFT-IFISTA with irregular focus grid can achieve better acoustic source identification performance than that with regular focus grid. Compared with FFT-NNLS and FFT-FISTA, FFT-IFISTA has a higher computational efficiency and better convergence. By selecting the appropriate weight coefficient, FFT-IFISTA can achieve excellent performance with the advantages of narrow mainlobe and spurious sources attenuation. Synthesizing various aspects of performance, for the array used in this paper, the recommended weight coefficient is 3.