An Extension of the Virtual Rotating Array Method Using Arbitrary Microphone Configurations for the Localization of Rotating Sound Sources

The characterization of rotating aeroacoustic sources using microphone array methods has been proven to be a useful tool. One technique to identify rotating sources is the virtual rotating array method. The method interpolates the pressure time data signals between the microphones in a stationary array to compensate the motion of the rotating sources. One major drawback of the method is the requirement of ring array geometries that are centred around the rotating axis. This contribution extends the virtual rotating array method to arbitrary microphone configurations. Two different ways to interpolate the time signals between the microphone locations are proposed. The first method constructs a mesh between the microphone positions using Delaunay-triangulation and interpolates over the mesh faces using piecewise linear functions. The second one is a meshless technique which is based on radial basis function interpolation. The methods are tested on synthetic array data from a benchmark test case as well as on experimental data obtained with a spiral array and a five-bladed fan.


Introduction
Beamforming and other array methods for acoustic source localisation have been extensively investigated for multichannel microphone measurements. The characterisation of rotating aeroacoustic sources using microphone array methods has been proven to be a useful tool. Various methods for identifying rotating noise source have been proposed. They can be classified as methods in time domain or in frequency domain. Rotating beamforming in time domain with a rotating focus point has been applied using the Rotating Source Identifier (ROSI) [1] for wind turbine and helicopter blades. The method was later combined with deconvolution algorithms as a hybrid time and frequency domain method [2] and applied to rotating sources [3]. An improved way to calculate the point spread function for the deconvolution of rotating sources was derived by Debrouwere et al. [4].
Rotating beamforming in the frequency domain was fist applied for induct measurements. The method solves the Green's function for a rotating acoustic source in a duct [5,6]. The free field solution of the Green's function for a rotating monopole in spherical coordinates was deduced by Pannert and Maier [7]. This method is known as modal decomposition (MD). Another method in the frequency domain to determine the location of acoustic sources is the virtual rotating array (VRA) method introduced by Herold and Sarradj [8]. It compensates the motion of the source by virtually rotating the array microphone positions. Since both methods work directly in the frequency domain the construction of the cross-spectral matrix and the application of subsequent beamforming, deconvolution and inverse techniques are easily possible.
Various comparisons of the three methods have been conducted [9][10][11]. The resulting delay and sum beamforming output for source localisation have been found to be almost equal for all the methods. Though, the computational cost of ROSI was found to be two to three orders of magnitude higher than that of the VRA method.
The major drawback of the VRA method is the restrictive microphone arrangement. It relies on evenly spaced microphone placements in one or multiple concentric rings around the center of the rotation. Geometries with uneven microphone placements been studied, but still rely on positioning on a ring [12]. To overcome this drawback, a new interpolation method for the virtual rotating array for arbitrary two dimensional and three dimensional microphone placements is proposed here.
The remainder of the paper is organised as follows. In Section 2, the beamforming methods and the implementation of the virtual rotation are described. Application of the methods on simulated benchmark data is presented in Section 3. The test of the methods with experimental data is discussed in Section 4. Section 5 concludes the paper.

Beamforming and Deconvolution Methods
The virtual rotating array method is based on beamforming in the frequency domain. The beamforming output (the squared sound pressures at one of the N focus points x s ) in the frequency with the cross-spectral matrix C of the microphone signals, the steering vector h and the hermitian of the steering vector h H . Here, the latter is chosen according to formulation III from Sarradj [13]: with the number of microphones M, the wavenumber k, r s,m and r s,l indicating the distances between the sources implicitly assumed stationary and the microphone locations and r s,0 between the the sources and the array centre location. This formulation has been shown theoretically to yield the correct source level [13]. In this contribution the (DAMAS) [14] deconvolution is used. The DAMAS algorithm assumes that the beamforming output b consists of the actual source distribution being convoluted by the point spread function from each focus point to all focus points P. The correct source distribution b can be revealed by solving the system of equations: The companding convex optimization problem can be written as: with · 2 denoting the L2-Norm and arg min being the arguments of the minima function which returns the set of points b for which b − Pb 2 attains the smallest value. The problem can be solved using a non negative least squares (NNLS) algorithm [15]. Since beamforming and deconvolution in this form only works for stationary sound sources, the motion of the rotating source has to be compensated.

Virtual Rotating Array
The motion compensation with a virtual rotating array (VRA) [8] tracks the angular position of the rotating source and rotates the virtual microphone positions at the same rate as the source. The signal at each virtual microphone positions is interpolated linearly from the two neighbouring microphones.
The VRA method with a ring array uses only the two neighbouring microphones to interpolate the signal at the virtual microphone positions. Since every point on the ring has exactly 2 neighbouring microphones this approach is suitable.
When using arbitrary microphone geometries the neighbouring microphones cannot be deduced from the virtual microphone position, which makes this interpolation method insufficient. For the extension of this algorithm to 2D and 3D array geometries an interpolation process that is based on the microphone positions has to be applied. Two methods are proposed here for the interpolation over irregularly distributed microphone positions: mesh based triangulation and radial basis function approximation.

Mesh Based Triangulation
For the first method, a two-or three-dimensional mesh needs to be constructed between all neighbouring microphones. A common way to derive such a mesh is the Delaunay triangulation [16]. The triangulation for the microphone positions x k is constructed in such a way that all circumcircles of all triangles have empty interiors, i.e., no other microphone position lies inside the circle. For this contribution all Delaunay triangulations have been calculated with the qhull package which is based on the quickhull algorithm by Barber et al. [17]. Figure 1 shows the Delaunay triangulation of a 2-D spiral microphone geometry. The array consists of 64 microphones which are ordered in a single sunflower spiral [18]. The spiral parameter H which defines the radial distribution of the microphones is chosen to be 1.0 and the parameter V defining the azimuthal distribution of points is set to 5.0. The parameters are chosen according to Sarradj [19]. The spiral array has an aperture of d ring = 1 m.
Since the virtual rotation is carried out in the angular domain the triangulation is also transferred in polar coordinates. The transferred Delauney triangulation is shown in Figure 1b. The angle of interest is set in the domain φ ∈ [−π, π]. To achieve angular periodicity of the interpolation, the mesh is extended to both sides with the repeating microphone position points in angular direction. The interpolation between the points in the Delaunay triangulation is achieved by barycentric interpolation between neighbouring points in the cylindrical-coordinate space. For a 2D Microphone geometry each point in the plane has three vertices of the triangle it is part of for the interpolation, for 3D geometries each point inside the triangulation is interpolated by the four tetrahedron vertices.
This algorithm is limited to interpolation inside the convex hull of the original array. Therefore the aperture of the virtual array has to be the same or smaller than the microphone array. This strict limitation can be overcome with the meshless algorithm which interpolates for the whole domain.

Radial Basis Functions
The meshless method for the virtual rotation uses radial basis function interpolation. A radial basis function (RBF) is an n-dimensional radially symmetric function Φ : R n → R which depends only on the Euclidean distance between the center of the RBF x k and the point of evaluation x. The Euclidean distance is defined as r = ( x − x k 2 ). There are various RBF's available in the literature that can be used for interpolation [20]. The most commonly used RBF's are listed in Table 1, with being a scaling factor.
The pressure time signal is reconstructed at any point x by a linear combination of the RBF's centred at the original microphone positions x i multiplied by the weight coefficients w i : and the weight coefficients are determined by solving the linear system with p 0 · · · p M being the discrete pressure values at the microphone positions. This linear system can be solved iteratively and results in new time signals p(x) at any virtual microphone position. This offers the possibility to extrapolate outside of the array aperture. Although, this only works for points close to the microphone positions. For points further away from the array, the extrapolation is not valid and strongly dependent on the chosen basis function.

Setup
The data used for evaluation are a part of a benchmark series for microphone array measurements. The benchmarks were put together to provide an evaluation standard for array method algorithms and their individual implementations [21]. The benchmark that is used in this contribution is the B11 benchmark [22]. The data of this benchmark comprise two sub-cases, each featuring simulated rotating monopole sources. The calculations of this contribution focus on the sub-case B. This test case consists of three rotating monopoles with a rotational rate that slightly varies over time. Source 1 and 2 are on a radius of 0.25 m with an angle difference ∆ϕ = 40 • . The source power of the second source is 3 dB below the power of the first source. The third source rotates on a radius of 0.125 m and has a 3 dB lower strength than the second source. The rotational rate is tracked by a tacho signal which provides one trigger pulse per revolution. The rotational angle is calculated from the tacho signal.
In order to evaluate the new virtual rotating array method the experimental setup in the Benchmark B11 has been extended. The microphone geometry was changed from an equidistant ring array to a planar spiral array. The number of microphones is kept at 64. The parameters of the spiral are the same as described in Section 2.3. Figure 2 shows the B11 benchmark setup as well as the extended case with the spiral array geometry and Table 2 shows the evaluation parameters for the calculation of the beamforming results.

Results
Two different interpolation methods are shown in this contribution. The modified test case with the spiral array geometry with Delaunay triangulation and barycentric interpolation as well as radial basis function interpolation with cubic basis functions. They are compared to the beamforming results without motion compensation and the VRA method using a ring array and linear interpolation.
The Figures 3-6 show the octave band source maps for three different center frequencies. The beamforming output was calculated using conventional beamforming in the frequency domain.
For the beamforming result without motion compensation the sources are 'smeared' over their trajectory along the ring. The octave band around 2500 Hz shows only one ring for all sources, while to two higher frequency bands show the separate third source at the 0.125 m radius. At 10 kHz, the effect of the beamforming sidelobes is clearly visible. Figure 4 shows the beamforming result for the B11 benchmark setup with a ring array and linear interpolation. The Beamforming results with a spiral array are shown in Figure 5 using motion compensation with Delaunay triangulation and in Figure 6 using motion compensation with radial basis function interpolation.
For all three cases the sources appear at the positions given in the benchmark B11, which indicates that the virtual rotation of the microphone position is feasible for reconstructing rotating source positions. The methods using the spiral array are able to interpolate the pressure values at the virtual microphone positions in a similar way. The width of the beamforming lobes depends strongly on the frequency. For the lowest octave at 2500 Hz the source positions are correct but the source separation of the three single sources is not possible because at this frequency the distance of the sources is below the Rayleigh criterion [23]. The source maps for the 5 kHz octave band are almost equal for the Delaunay and RBF interpolation. The three sources are at the correct positions and distinguishable. For very high frequencies around the 10 kHz band the sources are also displayed correctly but the RBF interpolation produces slightly more artefacts than the Delauney triangulation outside the source regions.
A quantitative comparison of the side lobe levels and the beam width for source 1 is given in Table 3. The beam with is defined by the diameter of the 3 dB level difference distance from source 1's maximum level. The dynamic is the level difference between source 1's maximum level and the level of the highest side lobe. In comparison the interpolation method using Delauney triangulation leads to a slightly better performances than the radial basis function interpolation in beam width as well as dynamic for this test case. The resulting beamforming output using the ring array produces less artefacts at higher frequencies than the spiral array geometry, but produce higher side lobes levels at all octave bands. Those differences are related to the array geometries. The ring array uses the full array aperture and therefore produces a smaller main lobe than the spiral array. The spiral array has a higher dynamic range than the ring array and produces less side lobes. The ring array is less sensitive to artefacts at high frequencies due to the smaller distances between the microphone positions compared to the spiral array.

Computation Time
The computation of the beamforming results was done on a Intel I5-7200U CPU. All results have been calculated with the recent version of the open-source software Acoular [24]. The calculation time for the beamforming without motion compensation was 90 s. Using the motion compensation with the ring array and linear interpolation increased the calculation time to 126 s. For the motion compensation with the Delaunay triangulation the calculation time was 190 s and for the radial basis function interpolation 545 s. While the delaunay triangulation is only 100 s slower than the beamforming without motion compensation, the RBF interpolation takes about five times longer. This is mainly because the triangulation is calculated only once but the linear system from the RBF interpolation has to be solved for every time step.

Setup
In order to test the different interpolation techniques acoustic measurements of a five-bladed fan have been conducted. The axial fan has a diameter of 800 mm and a rated rotational speed of 1020 revolutions per minute.
The array parameters of the spiral array that is used for the measurement are the same as described in Section 2.3. The array aperture of the array is d spiral = 1.5 m. The type of microphones used for both arrays is GRAS 40PL-1 Short CCP and the distance between the array center and the fan was 0.74 m. In addition to the collected sound pressures, the rotational speed of the fan was tracked using a laser trigger. The trigger-per-revolution signal was synchronised to the pressure measurement to take into account small variations in the rotational speed of the fan. The measurement time was 40 s with a sampling frequency of 51,200 Hz.
The beamforming output was calculated using conventional beamforming in the frequency domain as well as DAMAS deconvolution. Diagonal removal was applied to the cross spectral matrix. Diagonal removal removes the main diagonal elements of the cross spectral matrix to neglect the contribution of noise contamination which is incoherent for all the array microphones, e.g., self-induced flow noise [25]. The measurement was conducted in the large anechoic chamber at TU Berlin. Figure 7 shows a photograph of the experimental setup and Figure 8 shows a sketch of the setups side view as well as the front view of the fan. The array measurement and data processing parameters are summarised in Table 4.

Results
The Figures 9 and 10 show the source distribution in terms of the sound pressure level at the array center for three different one-third octave bands that are calculated with conventional beamforming. The displayed dynamic range of the sound pressure level is set to 15 dB. For all three frequency bands the source positions at the tip of the rotor blade are the same for both interpolation algorithms. Previous studies on fan noise localisation have already shown that these source locations are plausible [26,27]. At high flow rate coefficients, as in this setup without static pressure differences between suction and the pressure side of the fan, the main noise sources are at the tips of the rotor blade. At 2.5 and 5 kHz both methods produce a similar beamforming output. The dynamic at the center of the grid is higher with RBF interpolation for 2.5 kHz. At the highest frequency band, sources at the center of the grid are visible. The spatial resolution of the source maps calculated with conventional beamforming is degraded by the influence of the point spread function of the microphone array. The absolute values of the conventional beamforming are overestimated because of the convoluted map. DAMAS deconvolution is applied to predict a correct source strength and enhance the localization. The source maps obtained with DAMAS deconvolution are shown in the Figures 11 and 12. The dynamic range of the maps is increased to 30 dB. For the 2.5 kHz band the distribution of the sources along the 5 blades is visible for both methods. The accuracy of the source positions at the blade tips seems more plausible using the barycentric interpolation. At 5 kHz the sources at the blade tips are reconstructed equally well with both interpolation methods. At the 10 kHz band sources in the center of the focus grid remain visible after deconvolution. The effect is stronger for the radial basis function interpolation. The localisation of the sources at the blade tips also is in better accordance with previous results using barycentric interpolation at this frequency.
In order to compare the quantitative results of the interpolation methods one-third octave band spectra were computed by integrating the sound pressure over the blade areas as well as over the whole focus area. The spectra are shown in Figure 13a using barycentric and in Figure 13b using radial basis function interpolation. Both spectra are based on the DAMAS deconvolution results. The contributions of the individual blades are shown in blue. The sum of the 5 blades is shown in green. The sound pressure level over the whole focus area is shown in black. For both interpolation algorithms the individual blades equally contribute to the overall sound pressure level. Moreover, the blades act as the main noise source, except for very low and very high frequencies. For high frequencies this stems from the sources at the center of the calculation area. For the lower frequencies the differences in the sound pressure level occur due to the insufficient spatial resolution of the source localisation.

Summary
An extension to the original virtual rotating array method has been proposed to overcome the limitation to use only ring arrays. The developed method is able to virtually rotate microphone array geometries with arbitrary microphone arrangements. Two different interpolation methods have been shown. The first one is based on Delaunay triangulation in combination with barycentric interpolation and the second one on radial basis function interpolation. The methods have been tested on simulated benchmark data and have been compared with the original VRA method using a ring array. The location of the simulated sources could be recovered with both methods, but the interpolation bases on Delauney triangulation archived favourable results. Measurements with a spiral array have been performed on a fan. The experimental data were evaluated with the two interpolation methods. Beamforming in frequency domain as well as DAMAS deconvolution has been applied to generate source maps and spectra. This study has shown that both methods have been able to compensate the motion of rotating sources and reconstruct the original source positions and levels with similar quality.