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

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 coeﬃcients 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 the l 1 norm optimization problem. Finally, the method was veriﬁed using several examples. The proposed method oﬀers 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 coeﬃcient 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-ﬁeld acoustic ﬁeld compared with the traditional random position sampling method reducing correlation in the transfer matrix.


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. e 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 l 2 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. e least-squares method based on the l 2 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]. e novel compressive sampling technique uses the sparsity of the signal to solve the underdetermined system of equations. e sparse solution is easier to obtain because the reconstruction algorithm is based on the l 0 or l 1 norm instead of the traditional leastsquares method based on the l 2 norm. us, 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]. e particle velocity maps will have sharp peaks in CVMESM [15]. e 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. e bidirectional Fourier series expansion of the equivalent source intensity on a rectangular virtual surface is used to establish a semianalytical 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.

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 S E inside the acoustic source. us, the acoustic pressure can be expressed using the following integral: where r E is the position vector of the equivalent source, q(r E ) is the strength of the equivalent source at r, and G(r, r E ) is Green's function. As shown in Figure 1, the position coordinates of field point r and equivalent source point r E in the rectangular coordinate system are r � (x, y, z) and r E � (x E , y E , z E ), respectively. To reconstruct the planar acoustic field, the virtual surface S E is arranged in a rectangular domain of size 2L x × 2L y . Expanding the equivalent source intensity q(r E ) in (1) into a Fourier series along the x-axis and y-axis of the selected virtual surface S E , z E is given, and it becomes the following: where the Fourier series coefficient is as follows: Substituting equation (2) into equation (1), the pressure field becomes

Shock and Vibration
In practical calculations, the summation range of the series in equation (3) needs to be truncated; that is,

Fast Fourier Wave Superposition Spectrum Method.
For convenience, the integral term in equation (4) can be expressed as follows: Dividing the integration range (−L x ∼ L x ) and (−L y ∼ L y ) in equation (5) into equal divisions of M and N, respectively, it becomes Using the trapezoidal formula to numerically solve equation (5), it becomes e two-dimensional discrete Fourier transform (DFT) and inverse DFT are defined as [18] follows:

Shock and Vibration
where the value of M × N samples is as follows: 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 mn (r; z E ) can be quickly calculated using the discrete fast Fourier transform.
It is worth mentioning that when the DFT is used to calculate G mn (r; z E ), 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 � (−m f ∼ + m f ) and n � (−n f ∼ + n f ), respectively, and only need to be adjusted by the periodicity of G mn (r; z E ) with respect to M and N. e fast Fourier wave superposition spectrum of acoustic pressure at any point can be obtained using equations (4) and (7): e number of summation truncation terms requires For acoustic field reconstruction in NAH, measured holographic acoustic pressure data are substituted into equation (9) and written in matrix form: where P h is the column vector of acoustic pressure on the holographic surface at sampling point M h and C is the column vector composed of M v � (2m f + 1)(2n f + 1) coefficients C mn (z E ) of the Fourier series of equivalent source intensity on the virtual surface S E . K is the transfer matrix of M h × M v between the acoustic pressure on the holographic sampling surface and the Fourier coefficient of the equivalent source intensity on the virtual surface S E . en, the least-squares method based on the l 2 norm is used to solve equation (10): where λ is the regularization parameter and ‖ · ‖ 2 is the l 2 norm. After solving the Fourier coefficient column vector C of equation (10), the acoustic pressure of any reconstructed surface in the acoustic field can be obtained as follows: where D is the transfer matrix between the acoustic pressure on any reconstructed surface and the Fourier coefficient C mn (z E ) of the equivalent source strength on the virtual surface S E .
From the above derivation, the equivalent source intensity q(r E ) on the virtual surface S E 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 l 2 norm must satisfy M h ≥ M v 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. erefore, equation (10) is generally a set of underdetermined equations, making it difficult to obtain satisfactory results using the least-squares method based on the l 2 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 l 0 or l 1 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. en, 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. e higher the order, the faster the attenuation; that is, the higher-order C mn becomes smaller and smaller. erefore, provided that the number of summation truncation terms m f and n f 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. us, the Fourier coefficient vector C mn must be a sparse vector with certain sparsity. As the number of summation truncation terms m f and n f increases, the sparsity is strengthened. erefore, the proposed fast Fourier wave superposition spectrum method can directly solve the acoustic field reconstruction using sparse sampling combined with the l 0 or l 1 norm of CS.
To intuitively explain the sparsity of the Fourier coefficient vector C mn , a simple analysis of the sparsity of the Fourier coefficient vector C mn of this kind of continuous vibration acoustic source is presented. e 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]. e 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. e 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].
e Gaussian white noise with a signal-to-noise ratio (SNR) of 20 dB was added to simulate the actual acoustic pressure. e 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 n f � 5 and m f � 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 � 2 9 . e equation was solved by basis pursuit de-noising (BPDN) based on the l 1 norm in the SPGL1 toolbox in MATLAB [21]. e amplitudes obtained using the C-ESM and the method in this study with a sampling rate of 1000 Hz are presented in Figure 2. e 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). e Fourier coefficient vector C mn 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 C mn has stronger sparsity for continuous structural vibration sources such as the simply supported plate. e larger the number of truncated terms m f and n f , the stronger the sparsity. erefore, compared with the C-ESM, the vector has better sparsity and can be better solved using the l 1 norm.

Principle of Compressed Sensing and Sparse
Sampling Matrix

Principle of Compressed Sensing.
Assuming an N-dimensional signal x ∈ R N×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: 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. e 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: where A CS is a matrix of M × N. Since the decomposition coefficient vector α has some sparsity, it can be optimally solved using the l 0 norm of CS: 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 C K N combinations, which will consume a lot of calculation time. Usually, the l 1 norm is used instead of the l 0 norm to obtain

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. ere is a random that is S � [0 · · · 0 · · · φ 1 · · · 0 · · · φ 2 · · · 0 · · · 0 · · · φ j · · · 0 · · · φ M h ]. 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 Let P h � SP h and H � SK; where P h is the column vector of M h × 1 and Η is the matrix of M h × M v . Equation (18) is a set of underdetermined equations that can be solved using the l 1 norm minimization method, that is, 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: 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. erefore, 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.

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. e size of the simply supported steel plate was 0.5 m × 0.5 m, and the thickness was 0.003 m. e 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]. e holographic surface was located 0.05 m above the plate, and its dimensions were consistent with the acoustic source surface. e 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 the z-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. e sampling points were randomly selected from 625 measurement points on the 25 × 25 holographic grid. e 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 n f � 30 and  Shock and Vibration m f � 40, respectively. According to equation (9), the sampling number along the x and y directions of the virtual surface was N f � M f � 2 9 . e reconstruction error is defined as follows: where P is the reconstructed acoustic pressure and P is the analytic expression of acoustic pressure. Figure 4 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 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. e 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. e proposed method obtains good results with the Gaussian matrix, Bernoulli matrix, and circular measurement matrix, and the error curves were almost coincident.

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 Section 3.1 were used. e excitation point was located on the top right of the plate (0.375, 0.375), as shown in Figure 6. e holographic surface was set at 0.06 m above the plate, and its size was consistent with the acoustic source surface. e reconstruction surface was 0.04 m above the plate, and the equivalent source surface was 0.01 m above the acoustic source surface. e sizes of the reconstruction surface and the equivalent source surface were consistent with those of the holographic surface.
In simulations, the same measured acoustic pressure on the holographic surface was used as in the previous example of Section 3.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. e number of sampling points was randomly selected from the 484 measuring points on the 22 × 22 holographic grid (intervals of 0.024 m). e 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. Figure 7 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 8 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. e 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.
According to Figures 5 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. e 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.

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. e 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 the l 1 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. e 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) e acoustic field reconstruction results of the traditional C-ESM and the proposed method for a simply supported plate source were compared. e 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. e sound field reconstruction of a simply supported plate under other arbitrary excitations can be regarded as the superposition of the excitations in this study. erefore, 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 (Table 1).

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.