# A Planar Acoustic Field Reconstruction Method Based on Fast Wave Superposition Spectrum and Sparse Sampling

Authors: Yang Zhang; Yu-Jiang Wang; Yu Xiang; Zi-Yu Shi; Shao-Jie Fan; Jing Lu
Journal: Shock and Vibration
 (2022)

---

## Abstract

Based on the fast Fourier wave superposition spectrum method, a new equivalent source method (ESM) with a sparse sampling technique is proposed. First, the equivalent source intensities are expanded on a rectangular virtual surface using a bidirectional Fourier series, resulting in a semi-analytic and half-numerical acoustic pressure expression. The Fourier coefficients result in good sparsity for continuous acoustic pressures from structural vibration sources, and the proposed sparse sampling method can further reduce correlation in the measurement matrix. Better results can be obtained by solving thel1 norm optimization problem. Finally, the method was verified using several examples. The proposed method offers two main advantages compared with the traditional compressive equivalent source method: (1) the unknown source intensity vector is expanded into a bidirectional Fourier series, thereby transforming an unknown source intensity vector into a sparse Fourier coefficient vector, which has better sparsity; (2) the proposed method constructs a random sampling matrix, which is expanded into a sparse sampling matrix by random distribution, thereby improving the reconstruction accuracy of planar near-field acoustic field compared with the traditional random position sampling method reducing correlation in the transfer matrix.


---

## Body

## 1. Introduction

Near-field acoustic holography (NAH) is a powerful technique for identifying and localizing acoustic sources and visualizing acoustic fields [1]. With the development of NAH technology, several new algorithms have been introduced including the spatial Fourier transform (SFT) [2], boundary element method (BEM) [3], and equivalent source method (ESM) [4]. Among them, the ESM has been widely studied and applied, mainly owing to its simple principle and strong adaptability [5]. According to the principle of the ESM, the acoustic field radiated by a vibrating body of arbitrary shape can be approximated as the superposition of the acoustic field generated by a series of virtual equivalent sources distributed inside the structure. The acoustic field is reconstructed using the weight coefficients of each equivalent source based on measured acoustic pressure on a holographic surface [6]. Since noise in measured data cannot be avoided, the least-squares method based on the l2 norm is usually used with the traditional ESM to obtain a stable solution. In practical applications of NAH, due to the measurement conditions and cost limitations, the acoustic reconstruction problem usually involves solving a set of underdetermined equations. The least-squares method based on the l2 norm is used to solve the undetermined system of equations, which makes it difficult to obtain an ideal reconstruction result.In recent years, compressed sensing (CS) has been widely used in signal and image processing [7] and has become a hot topic of research [8, 9]. The novel compressive sampling technique uses the sparsity of the signal to solve the underdetermined system of equations. The sparse solution is easier to obtain because the reconstruction algorithm is based on the l0 or l1 norm instead of the traditional least-squares method based on the l2 norm. Thus, good reconstruction accuracy is guaranteed with fewer measurement points. In 2012, Chardon et al. [10] introduced CS into NAH for the first time and experimentally demonstrated higher accuracy compared with the traditional method when reducing the measurement points. Fernandez and Xenaki [11] proposed the compressive equivalent source method (C-ESM) by combining CS with ESM and analyzed the influence of column correlation in the transfer matrix on the reconstruction result [12]. On the basis of the C-ESM, Hu et al. [13] studied the reconstruction of sparse sound field through the sparse basis obtained by singular value decomposition of the transfer matrix; a fast sparse acoustic field reconstruction method was proposed to combine the Bayesian compressed sensing with the sparse basis function in the following study [14]. In addition to the C-ESM, there are many methods combining compressed sensing with equivalent source method, for example, compressed fused model equivalent source method (CFMESM) [15], compressed velocity-mode equivalent source method (CVMESM) [16], and fused total generalized variation (FTGV) [17]. The particle velocity maps will have sharp peaks in CVMESM [15]. The mode scaling factor was set in CFMESM, and the need for the parameter is a weakness of the method [16]. FTGV combines sparsity in the second derivatives with sparsity in the amplitudes, but FTGV must be solved using the CVX toolbox [17]. In previous studies, sampling points (microphones) are often randomly arranged to satisfy noncorrelation between the observation matrix and the sparse basis. Random sampling can reduce correlation in the measurement matrix through random positioning. Moreover, since previous methods are based on the assumed sparsity of the equivalent source intensities, it is difficult to obtain accurate prior information of the acoustic source when the type of acoustic source is unknown, and the sparsity requirement is not guaranteed. Using corresponding sparse sampling to the transform matrix, we can not only randomly select the sampling position to reduce correlation in the measurement matrix but also use the random coefficient matrix to further reduce correlation.In this study, a planar acoustic field reconstruction method based on the fast wave superposition spectrum and the sparse sampling matrix is proposed. The bidirectional Fourier series expansion of the equivalent source intensity on a rectangular virtual surface is used to establish a semi-analytical and semi-numerical Fourier series expansion form, and the sparsity of the Fourier coefficient vector is analyzed. To reduce the correlation of the measurement matrix further, a sampling matrix with random spatial positioning and random summation coefficients is introduced. Finally, simulations were performed to compare the proposed method and the traditional C-ESM for planar acoustic field reconstruction.


## 2. Planar Acoustic Field Reconstruction Based on Fast Fourier Wave Superposition Spectrum

### 2.1. Fourier Series Expansion of Acoustic Pressure

From the wave superposition method (WSM) [4], the acoustic pressure at any field point r outside the acoustic source can be superposed by a series of virtual equivalent sources in the virtual surface SE inside the acoustic source. Thus, the acoustic pressure can be expressed using the following integral:(1)pr=∫SEqrEGr,rEdSE,where rE is the position vector of the equivalent source, qrE is the strength of the equivalent source at r, and Gr,rE is Green’s function. As shown in Figure 1, the position coordinates of field point r and equivalent source point rE in the rectangular coordinate system are r=x,y,z and rE=xE,yE,zE, respectively. To reconstruct the planar acoustic field, the virtual surface SE is arranged in a rectangular domain of size 2Lx×2Ly.Figure 1
Field point and source point.Expanding the equivalent source intensityqrE in (1) into a Fourier series along the x-axis and y-axis of the selected virtual surface SE, zE is given, and it becomes the following:(2)qrE=∑m=−∞+∞∑n=−∞+∞CmnzEexp−imπxELxexp−inπyELy,xE∈−Lx,Lx,y∈−Ly,Ly,where the Fourier series coefficient is as follows: CmnzE=1/4LxLy∬SEqrEexpimπxE/LxexpinπyE/LydxEdyE.Substituting equation (2) into equation (1), the pressure field becomes(3)pr=∑m=−∞+∞∑n=−∞+∞CmnzE∫−LyLy∫−LxLxGr;xE,yE,zEexp−imπxELxexp−inπyELydxEdyE.In practical calculations, the summation range of the series in equation (3) needs to be truncated; that is,(4)pr=∑m=−mfmf∑n=−nfnfCmnzE∫−LyLy∫−LxLxGr;xE,yE,zEexp−imπxELxexp−inπyELydxEdyE.


### 2.2. Fast Fourier Wave Superposition Spectrum Method

For convenience, the integral term in equation (4) can be expressed as follows:(5)G¯mnr;zE=∫−LyLy∫−LxLxGr;xE,yE,zEexp−imπxELxexp−inπyELydxEdyE.Dividing the integration range−Lx∼Lx and −Ly∼Ly in equation (5) into equal divisions of M and N, respectively, it becomes(6)ΔxE=2LxM,xE=k1ΔxE,k1=0,1,2…M−1,ΔyE=2LyN,yE=k2ΔyE,k2=0,1,2…N−1.Using the trapezoidal formula to numerically solve equation (5), it becomes(7)G¯mnr;zE=4LxLyMN∑k1=0M−1∑k2=0N−1Gr;2k1LxM,2k2LyN,zEexp−i2πMmk1exp−i2πNnk2m=0,1,2…M−1;n=0,1,2…N−1.The two-dimensional discrete Fourier transform (DFT) and inverse DFT are defined as [18] follows:(8)Xm,n=∑k1=0M−1∑k2=0N−1xk1,k2exp−i2πMmk1exp−i2πNnk2,m=0,1,…,M−1;n=0,1,…,N−1xk1,k2=1MN∑m=0M−1∑n=0N−1Xm,nexpi2πMmk1expi2πNnk2,k1=0,1,…,M−1;k2=0,1,…,N−1,where the value of M×N samples is as follows: xk1,k2k1=0,1,…,M−1;k2=0,1,…,N−1.It can be seen from equation (8) that equation (7) is the standard two-dimensional DFT. When M and N are positive integer powers of 2, G¯mnr;zE can be quickly calculated using the discrete fast Fourier transform.It is worth mentioning that when the DFT is used to calculateG¯mnr;zE, the range of m and n is m=0∼M−1 and n=0∼N−1, respectively. However, the summation ranges in equation (4) are m=−mf∼+mf and n=−nf∼+nf, respectively, and only need to be adjusted by the periodicity of G¯mnr;zE with respect to M and N.The fast Fourier wave superposition spectrum of acoustic pressure at any point can be obtained using equations (4) and (7):(9)pr=∑m=−mfmf∑n=−nfnfCmnzEG¯mnr;zE.The number of summation truncation terms requires2mf+1≤M and 2nf+1≤N.For acoustic field reconstruction in NAH, measured holographic acoustic pressure data are substituted into equation (9) and written in matrix form:(10)Ph=KC,where Ph is the column vector of acoustic pressure on the holographic surface at sampling point Mh and C is the column vector composed of Mv=2mf+12nf+1 coefficients CmnzE of the Fourier series of equivalent source intensity on the virtual surface SE. K is the transfer matrix of Mh×Mv between the acoustic pressure on the holographic sampling surface and the Fourier coefficient of the equivalent source intensity on the virtual surface SE. Then, the least-squares method based on the l2 norm is used to solve equation (10):(11)argminPh−KC2+λC2,where λ is the regularization parameter and ·2 is the l2 norm.After solving the Fourier coefficient column vectorC of equation (10), the acoustic pressure of any reconstructed surface in the acoustic field can be obtained as follows:(12)Pr=DC,where D is the transfer matrix between the acoustic pressure on any reconstructed surface and the Fourier coefficient CmnzE of the equivalent source strength on the virtual surface SE.From the above derivation, the equivalent source intensityqrE on the virtual surface SE can be expanded into a Fourier series and the acoustic pressure expression can be obtained using the fast Fourier wave superposition spectrum. Since this formula is semi-analytical and semi-numerical, the solution accuracy will be higher than that of the traditional ESM [19]. To improve the calculation results, the traditional least-squares method based on the l2 norm must satisfy Mh≥Mv and equation (10) must be a set of overdetermined equations.In practical applications of NAH, the number of holographic measuring points is often insufficient to meet the above requirements due to measurement conditions and cost constraints. Therefore, equation (10) is generally a set of underdetermined equations, making it difficult to obtain satisfactory results using the least-squares method based on the l2 norm. However, if the sparsity of the Fourier coefficient vector C of equivalent source intensity is known, sparse sampling can be carried out on the holographic surface and the underdetermined equations can be solved using the l0 or l1 norm methods of CS, which can obtain a better solution result.According to the theory of structural dynamics, the dynamic response (vibration velocity) of a structure under an accidental load excitation is mainly composed of a superposition of low-order modes. Then, the Rayleigh integral shows that the radiated acoustic field is also mainly concentrated in the low-order modes. In other words, the contribution of structural vibration to the acoustic field is mainly concentrated in the low-order modes, whereas the high-order modes are similar to the evanescent modes in the SFT. The higher the order, the faster the attenuation; that is, the higher-orderCmn becomes smaller and smaller. Therefore, provided that the number of summation truncation terms mf and nf of equation (9) is large enough, the acoustic field represented by this formula will contain both low-order modes with strong radiation ability and evanescent modes. Thus, the Fourier coefficient vector Cmn must be a sparse vector with certain sparsity. As the number of summation truncation terms mf and nf increases, the sparsity is strengthened. Therefore, the proposed fast Fourier wave superposition spectrum method can directly solve the acoustic field reconstruction using sparse sampling combined with the l0 or l1 norm of CS.To intuitively explain the sparsity of the Fourier coefficient vectorCmn, a simple analysis of the sparsity of the Fourier coefficient vector Cmn of this kind of continuous vibration acoustic source is presented. The simulation calculations were performed with the vibration of a simply supported plate as the acoustic source and compared with the results of the C-ESM [12]. The dimensions of the simply supported plate were 0.5 m × 0.5 m × 0.03 m (length × width × height). A harmonic excitation force of 1 N was applied at the center of the plate, and the excitation point was located at the center of the simply supported plate. The sampled acoustic pressure on the holographic surface was calculated by the Rayleigh integral method to obtain the theoretical acoustic pressure at the measuring point [20]. The Gaussian white noise with a signal-to-noise ratio (SNR) of 20 dB was added to simulate the actual acoustic pressure. The holographic surface was 0.1 m above the simply supported plate, the reconstruction surface was 0.05 m above the simply supported plate, and the virtual equivalent source surface was arranged on the simply supported plate. Sizes of the holographic surface, reconstruction surface, and virtual equivalent source surface are consistent with those of the simply supported plate. In the C-ESM, 81 measuring points (a 9 × 9 uniform array) were set on the holographic surface and the number of equivalent sources was 441. According to the relationship between the number of summation truncation terms and the sound wavelength, the number of summation and truncation terms in the x and y directions was nf=5 and mf=20, respectively. According to equation (9), the number of integral segments in the x and y directions of the virtual surface was set as N=M=29. The equation was solved by basis pursuit de-noising (BPDN) based on the l1 norm in the SPGL1 toolbox in MATLAB [21].The amplitudes obtained using the C-ESM and the method in this study with a sampling rate of 1000 Hz are presented in Figure2. The virtual source intensity Q obtained by the C-ESM has large amplitudes at most of the source intensity sequence points and the sparsity is poor, as shown in Figure 2(a). The Fourier coefficient vector Cmn of this study has only a small number of large amplitudes and the rest are close to zero, as shown in Figure 2(b), which indicates that the Fourier coefficient vector Cmn has stronger sparsity for continuous structural vibration sources such as the simply supported plate. The larger the number of truncated terms mf and nf, the stronger the sparsity. Therefore, compared with the C-ESM, the vector has better sparsity and can be better solved using the l1 norm.Figure 2
Amplitude obtained using the compressive equivalent source method and the method with a sampling rate of 1000 Hz. (a) Compressive equivalent source method. (b) Method of this study.
(a)(b)


## 2.1. Fourier Series Expansion of Acoustic Pressure

From the wave superposition method (WSM) [4], the acoustic pressure at any field point r outside the acoustic source can be superposed by a series of virtual equivalent sources in the virtual surface SE inside the acoustic source. Thus, the acoustic pressure can be expressed using the following integral:(1)pr=∫SEqrEGr,rEdSE,where rE is the position vector of the equivalent source, qrE is the strength of the equivalent source at r, and Gr,rE is Green’s function. As shown in Figure 1, the position coordinates of field point r and equivalent source point rE in the rectangular coordinate system are r=x,y,z and rE=xE,yE,zE, respectively. To reconstruct the planar acoustic field, the virtual surface SE is arranged in a rectangular domain of size 2Lx×2Ly.Figure 1
Field point and source point.Expanding the equivalent source intensityqrE in (1) into a Fourier series along the x-axis and y-axis of the selected virtual surface SE, zE is given, and it becomes the following:(2)qrE=∑m=−∞+∞∑n=−∞+∞CmnzEexp−imπxELxexp−inπyELy,xE∈−Lx,Lx,y∈−Ly,Ly,where the Fourier series coefficient is as follows: CmnzE=1/4LxLy∬SEqrEexpimπxE/LxexpinπyE/LydxEdyE.Substituting equation (2) into equation (1), the pressure field becomes(3)pr=∑m=−∞+∞∑n=−∞+∞CmnzE∫−LyLy∫−LxLxGr;xE,yE,zEexp−imπxELxexp−inπyELydxEdyE.In practical calculations, the summation range of the series in equation (3) needs to be truncated; that is,(4)pr=∑m=−mfmf∑n=−nfnfCmnzE∫−LyLy∫−LxLxGr;xE,yE,zEexp−imπxELxexp−inπyELydxEdyE.


## 2.2. Fast Fourier Wave Superposition Spectrum Method

For convenience, the integral term in equation (4) can be expressed as follows:(5)G¯mnr;zE=∫−LyLy∫−LxLxGr;xE,yE,zEexp−imπxELxexp−inπyELydxEdyE.Dividing the integration range−Lx∼Lx and −Ly∼Ly in equation (5) into equal divisions of M and N, respectively, it becomes(6)ΔxE=2LxM,xE=k1ΔxE,k1=0,1,2…M−1,ΔyE=2LyN,yE=k2ΔyE,k2=0,1,2…N−1.Using the trapezoidal formula to numerically solve equation (5), it becomes(7)G¯mnr;zE=4LxLyMN∑k1=0M−1∑k2=0N−1Gr;2k1LxM,2k2LyN,zEexp−i2πMmk1exp−i2πNnk2m=0,1,2…M−1;n=0,1,2…N−1.The two-dimensional discrete Fourier transform (DFT) and inverse DFT are defined as [18] follows:(8)Xm,n=∑k1=0M−1∑k2=0N−1xk1,k2exp−i2πMmk1exp−i2πNnk2,m=0,1,…,M−1;n=0,1,…,N−1xk1,k2=1MN∑m=0M−1∑n=0N−1Xm,nexpi2πMmk1expi2πNnk2,k1=0,1,…,M−1;k2=0,1,…,N−1,where the value of M×N samples is as follows: xk1,k2k1=0,1,…,M−1;k2=0,1,…,N−1.It can be seen from equation (8) that equation (7) is the standard two-dimensional DFT. When M and N are positive integer powers of 2, G¯mnr;zE can be quickly calculated using the discrete fast Fourier transform.It is worth mentioning that when the DFT is used to calculateG¯mnr;zE, the range of m and n is m=0∼M−1 and n=0∼N−1, respectively. However, the summation ranges in equation (4) are m=−mf∼+mf and n=−nf∼+nf, respectively, and only need to be adjusted by the periodicity of G¯mnr;zE with respect to M and N.The fast Fourier wave superposition spectrum of acoustic pressure at any point can be obtained using equations (4) and (7):(9)pr=∑m=−mfmf∑n=−nfnfCmnzEG¯mnr;zE.The number of summation truncation terms requires2mf+1≤M and 2nf+1≤N.For acoustic field reconstruction in NAH, measured holographic acoustic pressure data are substituted into equation (9) and written in matrix form:(10)Ph=KC,where Ph is the column vector of acoustic pressure on the holographic surface at sampling point Mh and C is the column vector composed of Mv=2mf+12nf+1 coefficients CmnzE of the Fourier series of equivalent source intensity on the virtual surface SE. K is the transfer matrix of Mh×Mv between the acoustic pressure on the holographic sampling surface and the Fourier coefficient of the equivalent source intensity on the virtual surface SE. Then, the least-squares method based on the l2 norm is used to solve equation (10):(11)argminPh−KC2+λC2,where λ is the regularization parameter and ·2 is the l2 norm.After solving the Fourier coefficient column vectorC of equation (10), the acoustic pressure of any reconstructed surface in the acoustic field can be obtained as follows:(12)Pr=DC,where D is the transfer matrix between the acoustic pressure on any reconstructed surface and the Fourier coefficient CmnzE of the equivalent source strength on the virtual surface SE.From the above derivation, the equivalent source intensityqrE on the virtual surface SE can be expanded into a Fourier series and the acoustic pressure expression can be obtained using the fast Fourier wave superposition spectrum. Since this formula is semi-analytical and semi-numerical, the solution accuracy will be higher than that of the traditional ESM [19]. To improve the calculation results, the traditional least-squares method based on the l2 norm must satisfy Mh≥Mv and equation (10) must be a set of overdetermined equations.In practical applications of NAH, the number of holographic measuring points is often insufficient to meet the above requirements due to measurement conditions and cost constraints. Therefore, equation (10) is generally a set of underdetermined equations, making it difficult to obtain satisfactory results using the least-squares method based on the l2 norm. However, if the sparsity of the Fourier coefficient vector C of equivalent source intensity is known, sparse sampling can be carried out on the holographic surface and the underdetermined equations can be solved using the l0 or l1 norm methods of CS, which can obtain a better solution result.According to the theory of structural dynamics, the dynamic response (vibration velocity) of a structure under an accidental load excitation is mainly composed of a superposition of low-order modes. Then, the Rayleigh integral shows that the radiated acoustic field is also mainly concentrated in the low-order modes. In other words, the contribution of structural vibration to the acoustic field is mainly concentrated in the low-order modes, whereas the high-order modes are similar to the evanescent modes in the SFT. The higher the order, the faster the attenuation; that is, the higher-orderCmn becomes smaller and smaller. Therefore, provided that the number of summation truncation terms mf and nf of equation (9) is large enough, the acoustic field represented by this formula will contain both low-order modes with strong radiation ability and evanescent modes. Thus, the Fourier coefficient vector Cmn must be a sparse vector with certain sparsity. As the number of summation truncation terms mf and nf increases, the sparsity is strengthened. Therefore, the proposed fast Fourier wave superposition spectrum method can directly solve the acoustic field reconstruction using sparse sampling combined with the l0 or l1 norm of CS.To intuitively explain the sparsity of the Fourier coefficient vectorCmn, a simple analysis of the sparsity of the Fourier coefficient vector Cmn of this kind of continuous vibration acoustic source is presented. The simulation calculations were performed with the vibration of a simply supported plate as the acoustic source and compared with the results of the C-ESM [12]. The dimensions of the simply supported plate were 0.5 m × 0.5 m × 0.03 m (length × width × height). A harmonic excitation force of 1 N was applied at the center of the plate, and the excitation point was located at the center of the simply supported plate. The sampled acoustic pressure on the holographic surface was calculated by the Rayleigh integral method to obtain the theoretical acoustic pressure at the measuring point [20]. The Gaussian white noise with a signal-to-noise ratio (SNR) of 20 dB was added to simulate the actual acoustic pressure. The holographic surface was 0.1 m above the simply supported plate, the reconstruction surface was 0.05 m above the simply supported plate, and the virtual equivalent source surface was arranged on the simply supported plate. Sizes of the holographic surface, reconstruction surface, and virtual equivalent source surface are consistent with those of the simply supported plate. In the C-ESM, 81 measuring points (a 9 × 9 uniform array) were set on the holographic surface and the number of equivalent sources was 441. According to the relationship between the number of summation truncation terms and the sound wavelength, the number of summation and truncation terms in the x and y directions was nf=5 and mf=20, respectively. According to equation (9), the number of integral segments in the x and y directions of the virtual surface was set as N=M=29. The equation was solved by basis pursuit de-noising (BPDN) based on the l1 norm in the SPGL1 toolbox in MATLAB [21].The amplitudes obtained using the C-ESM and the method in this study with a sampling rate of 1000 Hz are presented in Figure2. The virtual source intensity Q obtained by the C-ESM has large amplitudes at most of the source intensity sequence points and the sparsity is poor, as shown in Figure 2(a). The Fourier coefficient vector Cmn of this study has only a small number of large amplitudes and the rest are close to zero, as shown in Figure 2(b), which indicates that the Fourier coefficient vector Cmn has stronger sparsity for continuous structural vibration sources such as the simply supported plate. The larger the number of truncated terms mf and nf, the stronger the sparsity. Therefore, compared with the C-ESM, the vector has better sparsity and can be better solved using the l1 norm.Figure 2
Amplitude obtained using the compressive equivalent source method and the method with a sampling rate of 1000 Hz. (a) Compressive equivalent source method. (b) Method of this study.
(a)(b)


## 3. Principle of Compressed Sensing and Sparse Sampling Matrix

### 3.1. Principle of Compressed Sensing

Assuming an N-dimensional signalx∈RN×1 can be linearly represented by a set of basis vectors ψi and defining a matrix composed of basis vectors as ψ=ψ1,ψ1,ψ3…,ψN, then the signal x can be expressed as follows:(13)x=∑i=1Nψiαi=ψα,where α is the decomposition coefficient. If there are only K nonzero values in the decomposition coefficient vector α and K≪N, the sparsity of the coefficient vector α is K and ψ is the sparse basis matrix of signal x. The signal x is sampled M times based on CS through the measurement matrix G, which is not related to the sparse basis matrix, to obtain the measurement value:(14)y=Gψα=ACSα,where ACS is a matrix of M×N. Since the decomposition coefficient vector α has some sparsity, it can be optimally solved using the l0 norm of CS:(15)argminy−ACSα2+λα0.However, the solution of equation (15) belongs to the NP-hard problem, which means that the nondeterministic polynomial (NP) cannot be solved by an exact algorithm, and an effective approximation algorithm for such problems must be sought. To obtain the correct solution, it is necessary to exhaustively enumerate CNK combinations, which will consume a lot of calculation time. Usually, the l1 norm is used instead of the l0 norm to obtain(16)argminy−ACSα2+λα1.


### 3.2. Construction of Sparse Sampling Matrix

In a previous study, the C-ESM reduced correlation in the measurement matrix (transfer matrix) by randomly arranging the sensor positions [12]. In this study, a random sampling matrix is constructed to reduce the correlation. There is a random matrix of Mh′×Mh: Φ=φ1…φj…φMh, where Mh is the number of actual sampling points (microphones). Extending this random matrix to a sampling matrix S of Mh′×Mv and randomly distributing Mh column vectors in Φ to S satisfies Mh≪Mv, that is S=0⋯0⋯φ1⋯0⋯φ2⋯0⋯0⋯φj⋯0⋯φMh. Zero column vectors of matrix S are the positions of no sampling points. Multiplying both sides of the sampling matrix in equation (10), it becomes(17)SPh=SKC.LetP¯h=SPh and H=SK;(18)P¯h=ΗC,where P¯h is the column vector of Mh×1 and Η is the matrix of Mh×Mv. Equation (18) is a set of underdetermined equations that can be solved using the l1 norm minimization method, that is,(19)argminP¯h−ΗC2+λC1.In contrast to the random arrangement of sampling points used in the traditional methods [12], this study proposes a new sparse sampling method. A random matrix is constructed and expanded into the sampling matrix. However, the method of selecting the random sampling matrix is important. Commonly used random sampling matrices are as follows:(1)
Random matrix: construct a random matrix ofMh′×Mh, with values of the array elements evenly distributed between 0 and 1. Measurement matrix P can be obtained by normalizing the column vector.(2)
Gaussian matrix [22]: construct a matrix Φ of Mh′×Mh such that each element of the matrix independently satisfies a Gaussian distribution with mean 0 and variance 1/Mh′ (strongly random, not related to most orthogonal bases).(3)
Bernoulli matrix [23]: construct a matrix Φ of Mh′×Mh such that every element in the matrix obeys a Bernoulli distribution, independently. Compared with the Gaussian matrix, the elements of the Bernoulli matrix are ±1, which are easier to implement and store in practical applications.(4)
Circulant matrix [24]: first generate an Mh dimensional random vector, circulate it MhMh′≤Mh times to construct the remaining Mh′−1 row vector, and finally normalize the column vector to obtain the measurement matrix Φ.(5)
Part Fourier matrix [25]: first generate an orthogonal matrix of Mh×Mh, then randomly select Mh′ row vectors of the orthogonal matrix, and finally normalize the column vectors of Mh′×Mh to obtain the measurement matrix Φ.In this study, the above five random sampling matrices are used, and sparse sampling matrices with random distributions are constructed. Compared with the traditional randomly selected sampling positions for reducing correlation in the measurement matrix, the randomly distributed coefficient matrix ensures that the sparse sampling matrix also has the randomness of the summation coefficients of the acoustic pressure at the measurement points. Therefore, the sparse sampling matrix contains not only the randomness of the position but also the randomness of the summation coefficients of the acoustic pressure at measurement points, which further reduces correlation in the sparse sampling matrix.


## 3.1. Principle of Compressed Sensing

Assuming an N-dimensional signalx∈RN×1 can be linearly represented by a set of basis vectors ψi and defining a matrix composed of basis vectors as ψ=ψ1,ψ1,ψ3…,ψN, then the signal x can be expressed as follows:(13)x=∑i=1Nψiαi=ψα,where α is the decomposition coefficient. If there are only K nonzero values in the decomposition coefficient vector α and K≪N, the sparsity of the coefficient vector α is K and ψ is the sparse basis matrix of signal x. The signal x is sampled M times based on CS through the measurement matrix G, which is not related to the sparse basis matrix, to obtain the measurement value:(14)y=Gψα=ACSα,where ACS is a matrix of M×N. Since the decomposition coefficient vector α has some sparsity, it can be optimally solved using the l0 norm of CS:(15)argminy−ACSα2+λα0.However, the solution of equation (15) belongs to the NP-hard problem, which means that the nondeterministic polynomial (NP) cannot be solved by an exact algorithm, and an effective approximation algorithm for such problems must be sought. To obtain the correct solution, it is necessary to exhaustively enumerate CNK combinations, which will consume a lot of calculation time. Usually, the l1 norm is used instead of the l0 norm to obtain(16)argminy−ACSα2+λα1.


## 3.2. Construction of Sparse Sampling Matrix

In a previous study, the C-ESM reduced correlation in the measurement matrix (transfer matrix) by randomly arranging the sensor positions [12]. In this study, a random sampling matrix is constructed to reduce the correlation. There is a random matrix of Mh′×Mh: Φ=φ1…φj…φMh, where Mh is the number of actual sampling points (microphones). Extending this random matrix to a sampling matrix S of Mh′×Mv and randomly distributing Mh column vectors in Φ to S satisfies Mh≪Mv, that is S=0⋯0⋯φ1⋯0⋯φ2⋯0⋯0⋯φj⋯0⋯φMh. Zero column vectors of matrix S are the positions of no sampling points. Multiplying both sides of the sampling matrix in equation (10), it becomes(17)SPh=SKC.LetP¯h=SPh and H=SK;(18)P¯h=ΗC,where P¯h is the column vector of Mh×1 and Η is the matrix of Mh×Mv. Equation (18) is a set of underdetermined equations that can be solved using the l1 norm minimization method, that is,(19)argminP¯h−ΗC2+λC1.In contrast to the random arrangement of sampling points used in the traditional methods [12], this study proposes a new sparse sampling method. A random matrix is constructed and expanded into the sampling matrix. However, the method of selecting the random sampling matrix is important. Commonly used random sampling matrices are as follows:(1)
Random matrix: construct a random matrix ofMh′×Mh, with values of the array elements evenly distributed between 0 and 1. Measurement matrix P can be obtained by normalizing the column vector.(2)
Gaussian matrix [22]: construct a matrix Φ of Mh′×Mh such that each element of the matrix independently satisfies a Gaussian distribution with mean 0 and variance 1/Mh′ (strongly random, not related to most orthogonal bases).(3)
Bernoulli matrix [23]: construct a matrix Φ of Mh′×Mh such that every element in the matrix obeys a Bernoulli distribution, independently. Compared with the Gaussian matrix, the elements of the Bernoulli matrix are ±1, which are easier to implement and store in practical applications.(4)
Circulant matrix [24]: first generate an Mh dimensional random vector, circulate it MhMh′≤Mh times to construct the remaining Mh′−1 row vector, and finally normalize the column vector to obtain the measurement matrix Φ.(5)
Part Fourier matrix [25]: first generate an orthogonal matrix of Mh×Mh, then randomly select Mh′ row vectors of the orthogonal matrix, and finally normalize the column vectors of Mh′×Mh to obtain the measurement matrix Φ.In this study, the above five random sampling matrices are used, and sparse sampling matrices with random distributions are constructed. Compared with the traditional randomly selected sampling positions for reducing correlation in the measurement matrix, the randomly distributed coefficient matrix ensures that the sparse sampling matrix also has the randomness of the summation coefficients of the acoustic pressure at the measurement points. Therefore, the sparse sampling matrix contains not only the randomness of the position but also the randomness of the summation coefficients of the acoustic pressure at measurement points, which further reduces correlation in the sparse sampling matrix.


## 4. Simulation Analysis

### 4.1. Simply Supported Plate Acoustic Source under Central Excitation

Since the C-ESM has been compared with FTGV, CVMESM, and CFMESM in reference [26], o reduce the amount of computation, the proposed method was only compared with the C-ESM method. To verify the correctness of the proposed method and the accuracy of the acoustic field reconstruction, a simply supported plate was used as a vibration source for simulation calculations. The size of the simply supported steel plate was 0.5 m × 0.5 m, and the thickness was 0.003 m. The plate was driven by a harmonic excitation force of 1 N, and the excitation point was located at the center of the plate. A theoretical radiation acoustic field was calculated using the Rayleigh first integral [20]. The holographic surface was located 0.05 m above the plate, and its dimensions were consistent with the acoustic source surface. The reconstruction surface was located 0.05 m from the holographic surface, and the equivalent source surface was arranged on the acoustic source surface.Dimensions of the reconstruction surface and equivalent source surface were consistent with the holographic surface, and the center of all three surfaces was on thez-axis with the center of the plate. In simulations, the actual measured acoustic pressure on the holographic surface was obtained by adding 30 dB of the Gaussian white noise to the theoretically calculated acoustic pressure. For the C-ESM, the number of virtual equivalent source points was 625 and the number of holographic measurement points was 81. The sampling points were randomly selected from 625 measurement points on the 25 × 25 holographic grid. The distribution of sampling point locations is shown in Figure 3. When the calculations were performed using the method of this study, the microphone array on the measurement surface was consistent with [12]. According to the relationship between the number of summation truncation terms and the sound wavelength, the number along the x and y directions of summation truncation terms was nf=30 and mf=40, respectively. According to equation (9), the sampling number along the x and y directions of the virtual surface was Nf=Mf=29.Figure 3
Distribution of holographic measurement points for compressive equivalent source method and proposed method.The reconstruction error is defined as follows:(20)ξ=P−P¯2P¯2×100%,where P is the reconstructed acoustic pressure and P¯ is the analytic expression of acoustic pressure.Figure4 shows the reconstructed acoustic pressure and theoretical acoustic pressure obtained using the C-ESM or proposed method with five different sparse sampling matrices at frequencies of 500 Hz, 1500 Hz, and 2500 Hz. When the frequency is 500 Hz, all reconstructed pressures match the analytical acoustic pressure, except for the proposed method, which combines part of the Fourier matrix. As the frequency increases, the C-ESM can no longer match the analytical acoustic pressure well, whereas the five different sparse sampling matrices are sufficient for obtaining better results with the proposed method. It can also be seen that the proposed method results better reconstruction effects when the Gaussian matrix, Bernoulli matrix, and circular measurement matrix were used as the random sampling measurement matrix and was in better agreement with the analytic expression of acoustic pressure.Figure 4
Reconstructed and theoretical acoustic pressure nephogram of six methods.Figure4 is only the reconstruction results at three specific frequencies: f = 500 Hz, 1500 Hz, and 2500 Hz. To analyze the reconstruction results at f = 100∼3000 Hz, Figure 5 shows the reconstruction error curves of the C-ESM and the proposed method with five different random sampling measurement matrices in the frequency band from 100 to 3000 Hz. The reconstruction error of the proposed method is lower than that of the C-ESM. Meanwhile, it is difficult to ensure stable reconstruction results in the computational frequency band using either the partial Fourier coefficient matrix or random matrix as the random sampling measurement matrix. The proposed method obtains good results with the Gaussian matrix, Bernoulli matrix, and circular measurement matrix, and the error curves were almost coincident.Figure 5
Acoustic pressure reconstruction error with frequency for various methods at 81 sampling points.


### 4.2. Simple Supported Plate Acoustic Source under Eccentric Excitation

To compare the acoustic field reconstruction of the proposed method and the C-ESM for a simply supported plate acoustic source under eccentric excitation, the parameters of the simply supported plate in Section3.1 were used. The excitation point was located on the top right of the plate (0.375, 0.375), as shown in Figure 6. The holographic surface was set at 0.06 m above the plate, and its size was consistent with the acoustic source surface. The reconstruction surface was 0.04 m above the plate, and the equivalent source surface was 0.01 m above the acoustic source surface. The sizes of the reconstruction surface and the equivalent source surface were consistent with those of the holographic surface.Figure 6
Eccentric excitation point.In simulations, the same measured acoustic pressure on the holographic surface was used as in the previous example of Section3.1, which was obtained by adding 30 dB of the Gaussian white noise to the theoretically calculated acoustic pressure. When the C-ESM is used for calculation, the number of virtual equivalent source points was 484 and the number of holographic measuring points was 64. The number of sampling points was randomly selected from the 484 measuring points on the 22 × 22 holographic grid (intervals of 0.024 m). The method used for selecting the random sampling matrix in this study was the same as that used in the C-ESM.Because of the particularity of central excitation, it is more universal to choose eccentric excitation as random excitation. Figure7 shows the theoretical acoustic pressures of the reconstructed surface at frequencies of 500 Hz, 1500 Hz, and 2000 Hz and acoustic pressure contours reconstructed using 64 sampling points using either the C-ESM or the proposed method with five different sparse sampling matrices. When f = 500 Hz, the reconstructed acoustic pressures of all methods are in good agreement with the analytical acoustic pressure and the error is less than 10%. As the frequency increases, the C-ESM no longer matches the analytical acoustic pressure well, but the proposed method has high accuracy with five different sparse sampling matrices. Using the Gaussian matrix, Bernoulli matrix, and circular measurement matrix as the random sampling measurement matrix, the proposed method can identify richer details of the acoustic field information.Figure 7
Reconstructed and theoretical acoustic pressure nephogram of six methods with 64 sampling points.Figure8 shows the reconstruction error curves of the C-ESM and the proposed method with five different sparse sampling measurement matrices in the frequency band of 100∼2000 Hz. The reconstruction error of the proposed method is lower than that of the C-ESM. Meanwhile, the reconstruction errors using the partial Fourier coefficient matrix or random matrix are higher than the other three methods; the proposed method obtains good results with the Gaussian matrix, Bernoulli matrix, and circular measurement matrix.Figure 8
Acoustic pressure reconstruction error with frequency for various acoustic field reconstruction methods with 64 sampling points.According to Figures5 and 8, the errors of the C-ESM and the proposed method with five different sparse sampling measurement matrices are less than 20% in the 100∼1400 Hz frequency range, and they are all acceptable error percentages. Due to different excitations, the calculation accuracy of five different sparse sampling measurement matrices is also different, the reconstruction results of the Gaussian matrix, Bernoulli matrix, and circular measurement matrix are better than that of the C-ESM, and the reconstruction results of the partial Fourier coefficient matrix or random matrix are similar to that of the C-ESM. However, the error of the C-ESM under central excitation exceeds 20%, the errors of the Gaussian matrix and circular measurement matrix are about 10∼16% in the 1400∼3000 Hz frequency range, and the errors of the Gaussian matrix and circular measurement matrix are acceptable error percentages. Because of the construction characteristic of the Gaussian matrix and circular measurement matrix, satisfactory results could be obtained even in high frequencies. The proposed method with the Gaussian matrix and circular measurement matrix can apply in further studies. In conclusion, the reconstruction results of the Gaussian matrix and circular measurement matrix are the best in five different sparse sampling measurement matrices, and it is further studied in other applications of near-field acoustic holography.Table 1
Comparison of full name and abbreviation in this article.
Full nameAbbreviationBoundary element methodBEMBasis pursuit de-noisingBPDNCompressive equivalent source methodC-ESMCompressed fused model equivalent source methodCFMESMCompressed velocity-mode equivalent source methodCVMESMCompressed sensingCSDiscrete Fourier transformDFTEquivalent source methodESMFused total generalized variationFTGVNear-field acoustic holographyNAHNondeterministic polynomialNPSpatial Fourier transformSFTSignal-to-noise ratioSNRWave superposition methodWSM


## 4.1. Simply Supported Plate Acoustic Source under Central Excitation

Since the C-ESM has been compared with FTGV, CVMESM, and CFMESM in reference [26], o reduce the amount of computation, the proposed method was only compared with the C-ESM method. To verify the correctness of the proposed method and the accuracy of the acoustic field reconstruction, a simply supported plate was used as a vibration source for simulation calculations. The size of the simply supported steel plate was 0.5 m × 0.5 m, and the thickness was 0.003 m. The plate was driven by a harmonic excitation force of 1 N, and the excitation point was located at the center of the plate. A theoretical radiation acoustic field was calculated using the Rayleigh first integral [20]. The holographic surface was located 0.05 m above the plate, and its dimensions were consistent with the acoustic source surface. The reconstruction surface was located 0.05 m from the holographic surface, and the equivalent source surface was arranged on the acoustic source surface.Dimensions of the reconstruction surface and equivalent source surface were consistent with the holographic surface, and the center of all three surfaces was on thez-axis with the center of the plate. In simulations, the actual measured acoustic pressure on the holographic surface was obtained by adding 30 dB of the Gaussian white noise to the theoretically calculated acoustic pressure. For the C-ESM, the number of virtual equivalent source points was 625 and the number of holographic measurement points was 81. The sampling points were randomly selected from 625 measurement points on the 25 × 25 holographic grid. The distribution of sampling point locations is shown in Figure 3. When the calculations were performed using the method of this study, the microphone array on the measurement surface was consistent with [12]. According to the relationship between the number of summation truncation terms and the sound wavelength, the number along the x and y directions of summation truncation terms was nf=30 and mf=40, respectively. According to equation (9), the sampling number along the x and y directions of the virtual surface was Nf=Mf=29.Figure 3
Distribution of holographic measurement points for compressive equivalent source method and proposed method.The reconstruction error is defined as follows:(20)ξ=P−P¯2P¯2×100%,where P is the reconstructed acoustic pressure and P¯ is the analytic expression of acoustic pressure.Figure4 shows the reconstructed acoustic pressure and theoretical acoustic pressure obtained using the C-ESM or proposed method with five different sparse sampling matrices at frequencies of 500 Hz, 1500 Hz, and 2500 Hz. When the frequency is 500 Hz, all reconstructed pressures match the analytical acoustic pressure, except for the proposed method, which combines part of the Fourier matrix. As the frequency increases, the C-ESM can no longer match the analytical acoustic pressure well, whereas the five different sparse sampling matrices are sufficient for obtaining better results with the proposed method. It can also be seen that the proposed method results better reconstruction effects when the Gaussian matrix, Bernoulli matrix, and circular measurement matrix were used as the random sampling measurement matrix and was in better agreement with the analytic expression of acoustic pressure.Figure 4
Reconstructed and theoretical acoustic pressure nephogram of six methods.Figure4 is only the reconstruction results at three specific frequencies: f = 500 Hz, 1500 Hz, and 2500 Hz. To analyze the reconstruction results at f = 100∼3000 Hz, Figure 5 shows the reconstruction error curves of the C-ESM and the proposed method with five different random sampling measurement matrices in the frequency band from 100 to 3000 Hz. The reconstruction error of the proposed method is lower than that of the C-ESM. Meanwhile, it is difficult to ensure stable reconstruction results in the computational frequency band using either the partial Fourier coefficient matrix or random matrix as the random sampling measurement matrix. The proposed method obtains good results with the Gaussian matrix, Bernoulli matrix, and circular measurement matrix, and the error curves were almost coincident.Figure 5
Acoustic pressure reconstruction error with frequency for various methods at 81 sampling points.


## 4.2. Simple Supported Plate Acoustic Source under Eccentric Excitation

To compare the acoustic field reconstruction of the proposed method and the C-ESM for a simply supported plate acoustic source under eccentric excitation, the parameters of the simply supported plate in Section3.1 were used. The excitation point was located on the top right of the plate (0.375, 0.375), as shown in Figure 6. The holographic surface was set at 0.06 m above the plate, and its size was consistent with the acoustic source surface. The reconstruction surface was 0.04 m above the plate, and the equivalent source surface was 0.01 m above the acoustic source surface. The sizes of the reconstruction surface and the equivalent source surface were consistent with those of the holographic surface.Figure 6
Eccentric excitation point.In simulations, the same measured acoustic pressure on the holographic surface was used as in the previous example of Section3.1, which was obtained by adding 30 dB of the Gaussian white noise to the theoretically calculated acoustic pressure. When the C-ESM is used for calculation, the number of virtual equivalent source points was 484 and the number of holographic measuring points was 64. The number of sampling points was randomly selected from the 484 measuring points on the 22 × 22 holographic grid (intervals of 0.024 m). The method used for selecting the random sampling matrix in this study was the same as that used in the C-ESM.Because of the particularity of central excitation, it is more universal to choose eccentric excitation as random excitation. Figure7 shows the theoretical acoustic pressures of the reconstructed surface at frequencies of 500 Hz, 1500 Hz, and 2000 Hz and acoustic pressure contours reconstructed using 64 sampling points using either the C-ESM or the proposed method with five different sparse sampling matrices. When f = 500 Hz, the reconstructed acoustic pressures of all methods are in good agreement with the analytical acoustic pressure and the error is less than 10%. As the frequency increases, the C-ESM no longer matches the analytical acoustic pressure well, but the proposed method has high accuracy with five different sparse sampling matrices. Using the Gaussian matrix, Bernoulli matrix, and circular measurement matrix as the random sampling measurement matrix, the proposed method can identify richer details of the acoustic field information.Figure 7
Reconstructed and theoretical acoustic pressure nephogram of six methods with 64 sampling points.Figure8 shows the reconstruction error curves of the C-ESM and the proposed method with five different sparse sampling measurement matrices in the frequency band of 100∼2000 Hz. The reconstruction error of the proposed method is lower than that of the C-ESM. Meanwhile, the reconstruction errors using the partial Fourier coefficient matrix or random matrix are higher than the other three methods; the proposed method obtains good results with the Gaussian matrix, Bernoulli matrix, and circular measurement matrix.Figure 8
Acoustic pressure reconstruction error with frequency for various acoustic field reconstruction methods with 64 sampling points.According to Figures5 and 8, the errors of the C-ESM and the proposed method with five different sparse sampling measurement matrices are less than 20% in the 100∼1400 Hz frequency range, and they are all acceptable error percentages. Due to different excitations, the calculation accuracy of five different sparse sampling measurement matrices is also different, the reconstruction results of the Gaussian matrix, Bernoulli matrix, and circular measurement matrix are better than that of the C-ESM, and the reconstruction results of the partial Fourier coefficient matrix or random matrix are similar to that of the C-ESM. However, the error of the C-ESM under central excitation exceeds 20%, the errors of the Gaussian matrix and circular measurement matrix are about 10∼16% in the 1400∼3000 Hz frequency range, and the errors of the Gaussian matrix and circular measurement matrix are acceptable error percentages. Because of the construction characteristic of the Gaussian matrix and circular measurement matrix, satisfactory results could be obtained even in high frequencies. The proposed method with the Gaussian matrix and circular measurement matrix can apply in further studies. In conclusion, the reconstruction results of the Gaussian matrix and circular measurement matrix are the best in five different sparse sampling measurement matrices, and it is further studied in other applications of near-field acoustic holography.Table 1
Comparison of full name and abbreviation in this article.
Full nameAbbreviationBoundary element methodBEMBasis pursuit de-noisingBPDNCompressive equivalent source methodC-ESMCompressed fused model equivalent source methodCFMESMCompressed velocity-mode equivalent source methodCVMESMCompressed sensingCSDiscrete Fourier transformDFTEquivalent source methodESMFused total generalized variationFTGVNear-field acoustic holographyNAHNondeterministic polynomialNPSpatial Fourier transformSFTSignal-to-noise ratioSNRWave superposition methodWSM


## 5. Conclusions

(1)
Based on the ESM, the source intensity can be expanded into a bidirectional Fourier series on the rectangular virtual surface, and a semi-analytical and semi-numerical expression of acoustic pressure can be obtained. The proposed method has higher accuracy compared with the traditional numerical discrete source methods. Furthermore, since this method converts the source intensity vector into the sparse Fourier coefficient vector, the vector to be solved also has stronger sparsity and can be more effectively solved using thel1 norm.(2)
Correlation in the measurement matrix is reduced by constructing a random sampling matrix, and the sparse sampling matrix is introduced by expanding it into the sampling matrix. The sparse sampling matrix not only includes the randomness of the sampling position, but also includes the randomness of the sum coefficient of the acoustic pressure at the measuring point, which further reduces correlation in the coefficient matrix.(3)
The acoustic field reconstruction results of the traditional C-ESM and the proposed method for a simply supported plate source were compared. The simulation results show that the reconstruction error of the proposed method is lower than that of the C-ESM, and this method has higher accuracy with the Gaussian matrix, Bernoulli matrix, and circular measurement matrix. In particular, the errors of the Gaussian matrix and circular measurement matrix under central excitation are about less than 16% in the 1400∼3000 Hz frequency range, which are at least 4% lower than the C-ESM.(4)
In this study, center excitation and arbitrary eccentric excitation are used for the excitation of a simply supported plate. The sound field reconstruction of a simply supported plate under other arbitrary excitations can be regarded as the superposition of the excitations in this study. Therefore, this method is of great significance for planar sound sources, but for the sound field reconstruction of nonplanar sound sources, especially the sound field reconstruction of rotating structures, this method can be further studied (Table1).



---

*Source: 6724100-2022-01-24.xml*