An Inverse Microphone Array Method for the Estimation of a Rotating Source Directivity

: Microphone arrays methods are useful for determining the location and magnitude of rotating acoustic sources. This work presents an approach to calculating a discrete directivity pattern of a rotating sound source using inverse microphone array methods. The proposed method is divided into three consecutive steps. Firstly, a virtual rotating array method that compensates for motion of the source is employed in order to calculate the cross-spectral matrix. Secondly, the source locations are determined by a covariance matrix ﬁtting approach. Finally, the sound source directivity is calculated using the inverse method SODIX on a reduced focus grid. Experimental validation and synthetic data from a simulation are used for the veriﬁcation of the method. For this purpose, a rotating parametric loudspeaker array with a controllable steering pattern is designed. Five different directivity patterns of the rotating source are compared. The proposed method compensates for source motion and is able to reconstruct the location as well the directivity pattern of the rotating beam source.


Introduction
Beamforming techniques for the source localization of rotating noise sources have been widely applied in various industry applications, such as axially blowing fans, wind turbines and rotor blades. The estimation of source positions and strengths with microphone array methods has been studied in numerous cases. Beamforming in the time domain with a rotating focus point has been introduced as ROSI and was applied to several test cases [1,2]. Frequency domain methods have been used to model propagation in a duct [3,4] and under free field conditions [5,6]. Common to all those techniques is the assumption of monopole sound sources. Although such beamforming techniques work well in locating sources with uniform directivities, they can perform poorly when used to reconstruct directional noise sources. Since most aerodynamic sound sources have non-uniform directivity patterns, beamforming algorithms for rotating sources using monopole transfer functions may result in the incorrect estimation of positions and strengths of sources. Sound propagation with non-uniform directivity has been paid much attention to in the field of aeroacoustic measurements over the recent years. A dipole directivity in the transfer function was introduced for beamforming algorithms [7] and applied to rotating sources using dipole sources with unknown orientations [8]. A correction method for ROSI with dipole sources was implemented in [9].
Inverse beamforming methods have also been used for non-uniform transfer functions. A generalized inverse beamforming technique with L 1 regularization and different multipole directivities has been implemented [10]. Pan et al. [11] used orthogonal beamforming in combination with the general inverse beamforming approach with several non-uniform transfer functions.
Another way to introduce the directivity of the sources is the inverse method source directivity modeling in the cross-spectral matrix (SODIX). This method is an extension of the spectral estimation method (SEM) from Blacodon and Elias [12], which was initially introduced by Funke and Michel [13,14]. It is mainly used to evaluate the directional sound radiation of turbofan engines. An adaption of the algorithm to a new solver and regularization was introduced by Oertwig [15,16]. In the case of SODIX, the directivity of the sources is not governed by the transfer function. Instead, the solution of the inverse problem directly gives the source strength towards the direction of every microphone in the array.
In this paper, a combination of the virtual rotating array method and the inverse method SODIX is developed and applied to rotating sources with a non-uniform directivity. Considering that SODIX is based on the cross-spectral matrix, motion compensation has to be applied before it can be used on moving sources. Therefore, the VRA method is used on the time data to achieve a static source distribution before CSM calculation. The method is validated by numerical simulations and experiments.
The remainder of the paper is organized as follows: In Section 2, the implementation of the virtual rotation, CMF and SODIX are described. The experimental setup of the rotating source and the microphone array measurements as well as the procedure to obtain simulated measurement data are presented in Section 3. To demonstrate the performance of the methods, the source distributions and directivities are calculated and discussed for experimental and simulated data in Section 4. Section 5 concludes the paper.

Virtual Rotating Array Method
In order to use microphone array methods in the frequency domain, the motion of the rotating sources needs to be compensated. The so-called virtual rotating array (VRA) method [6] tracks the angular position of the rotating source and rotates the virtual microphone positions accordingly.
In this work, a planar spiral array is used for source localization. The same array geometry was tested with different interpolation methods by Jekosch [17]. A grid-based technique and a meshless technique were developed for arbitrary planar array configuration. For the source localization with this array geometry, the grid-based interpolation provided favorable results and is therefore used in this work. The grid-based interpolation is achieved by using a Delaunay triangulation inside the convex hull of the array coordinates in cylindrical coordinates. The sound pressure signal at the virtual rotated microphone positions is then interpolated from the three neighboring microphones in the plane.

Inverse Methods Based on CSM Modeling
The inverse microphone array methods are based on the cross-spectral matrix, which can be calculated from the measured pressure signals p m of M microphones via The cross-spectral matrix can be computed using Welch's method [18] in which the signal is divided into K blocks that are then transformed into the frequency domain by means of an FFT. For each discrete frequency, the complex-valued sound pressure values p m are then averaged over all blocks which results in a matrix C that contains the cross-spectra for every frequency of the FFT.
The spectral estimation method (SEM) [12] and covariance matrix fitting (CMF) [19] both assume that the pressure signal p m at the m-th microphone is caused by acoustic sources q j at the j-th focus grid point: The aim of these methods is to model a cross-spectral matrix C mod with a known sound propagation function a jm from the grid of possible sources to the microphones and the unknown source pressure amplitudes q j at the grid points j. The source matrix D j is a diagonal matrix containing the source strengths D j = diag(q j ): for monopole Green's functions, the sound propagation a jm is described by with r jm being the distances from the microphones to the grid points and c 0 being the speed of sound. To minimize the difference between the modeled and the measured CSM, a cost function F(D) can be established. The cost function F CMF for the optimization problem is given by The problem formulation for SODIX extends the monopole assumption to discrete directivities in the direction of every microphone. The matrix D jm contains the source strengths for every focus grid pointing in the direction of each microphone.
This cost function needs to be solved iteratively by a minimization algorithm. An analytical derivative of the cost function was provided by Funke [20]. The present paper uses the L-BFGS-B gradient solver [21] as proposed by Oertwig [15], as it already fulfills the bound constraints D jm ≥ 0. The derivative of the SODIX cost function with the bound constraint D jm ≥ 0 is given by As an inverse problem, the stability of SODIX can be improved with regularization methods. The regularization increases the stability of the optimization by reducing the ill-posedness. Therefore, a penalty function in the form of a matrix p-norm regularization factor D jm p with a regularization parameter α is added to the SODIX cost function: The derivative of the p-norm x p can be calculated via with sgn(x) being the sign function. This leads to a SODIX derivative with an additional regularization term: and the regularization term can be reduced when applying L 1 regularization: The sound localization is calculated with the model for the covariance matrix fitting described in Equation (6). The SODIX modeling described in Equations (9) and (12) is used for the calculation of the discrete source directivities.

Rotating Loudspeaker Array
With a focus on non-uniform rotating sources, a parametric loudspeaker array was constructed for the measurement. It consists of eight Visaton BF37 speakers that are mounted on a rotating arm. The rotating arm is 800 mm long, includes a power bank and two micro controllers and serves as a housing for the speakers. The micro controllers handle real-time audio output and provide access to the test rig via Wi-Fi.
The wireless communication is enabled by a micro controller (ESP32) that controls the audio signal settings. Another micro controller (Nucleo G431KB) serves as the digital signal processor. Communication between the micro controllers is implemented with the serial peripheral interface. The digital audio signal is transmitted to the eight audio amplifiers via the serial audio interface in the form of time division multiplexing. All eight speakers can be activated and delayed individually. The distance between the speakers is 40 mm and the steering angle of the main lobe can be changed in a 10 degree resolution by shifting the signal between the speakers by 1 sample. A controlled AC motor with 1.05 kW power, which is mounted to a height-adjustable test stand frame, is used to ensure a rotation with a constant rotational speed. The motor can rotate the speaker array with up to 400 revolutions per minute. The parametric loudspeaker array on the rotating beam is shown in Figure 1.

Array Measurement
The acoustic measurements were conducted in the large anechoic chamber at TU Berlin. For the array measurement, a planar array with 64 channels is used. The microphones are ordered in a sunflower spiral [22] and the parameters of the sunflower spiral are chosen to be H = 1.0 and V = 5.0 according to Sarradj [23] while the array aperture of the spiral is d = 1.5 m. The distribution of the microphone sensors is shown in Figure 2b). GRAS 40PL-1 Short CCP microphones are used for the array. The rotational axis of the speaker array was aligned to the center of the array. The distance between the microphone array plane and the speakers is 0.991 m. In addition to the sound pressures, the rotational speed of the speaker array was tracked using a laser trigger. One trigger is logged per revolution and the trigger signal was synchronized to the pressure data recording. The measurement time was 40 seconds at a sampling frequency of 51,200 Hz. Figure 3 shows a photograph of the experimental setup. The array measurement and data processing parameters are summarized in Table 1.   Five different cases are measured for comparison. The first case is measured without rotation and uses the same white noise signal for all eight speakers. In the remaining four cases, the speaker array rotates with an average speed of 205 revolutions per minute. For each case, the main radiation direction is changed by delaying the signal between the speakers. This way, the radiation can be changed in both directions by starting to delay the signals at either the first or last speaker in the line. The five cases are summarized in Table 2. The different radiation angles result in a change in directivity maxima towards the array plane. Figure 2a shows the side view of the setup including the main radiation direction for the four cases, and Figure 2b shows the sensor arrangement in the microphone array with the same radiation maxima.

Simulation
In order to compare to the measurement data with an ideal synthetic setup, test case 2 was simulated using the same parameters as those in the measurement. The rotating speaker array is represented by eight monopole sources on circular trajectories with the same radius as that in the experiment. All sources emit the same white noise signal with an RMS source strength of 1 Pa. The distance between the sources is 40 mm, and the rotational speed amounts to 205 revolutions per minute. The microphone array geometry and the distance between the array and the source plane are kept the same as those in the measurement. The sampling rate and measurement duration are also left unchanged. However, the simulated data do not contain any additional noise from aerodynamic sources or the motor. The simulated case is referred to as test case 6.

Results
To obtain the sound source distribution and directivities from the measurement data, the linear two-dimensional VRA method was applied. The cross-spectral matrix was calculated from the interpolated pressure data using a block size of 1024 samples and a 50% overlapping von Hann windowing function. Diagonal removal of the cross-spectral matrix was applied to reduce noise. The source locations were calculated with the CMF method using a LASSO solver on the full calculation grid consisting of 1681 points. L 1 regularization was applied to the measurement data. To correctly tune the value of the regularization parameter α, the Bayesian information criterion [24] was used in order to archive a sparse source distribution while simultaneously preserving all of the significant components.
Since the resulting source distribution is sparse and only shows significant sources at the area of the line source, the calculation grid for the SODIX method can be reduced to 60 grid points. This extenuates the ill-posedness of the optimization problem to the extent that regularization can be omitted entirely. The SODIX method was used with a fixed number of 100 iterations per discrete frequency. The calculation parameters can be found in Table 3. The calculations and simulation of test case 6 were performed with the open-source beamforming library Acoular [25]. The dynamic range was set to 15 dB for the source maps and 30 dB for the directivity distributions. Because the speakers of the loudspeaker array emit coherent white noise, the source distribution along the speakers does not appear as a perfect line source as would be expected for incoherent noise sources.  Figures 4-9 show the source distribution in terms of the sound pressure level at the array center and the source directivities towards the array microphones at the octave band around 5 kHz. On the left hand side of Figures 4-9, the rotating speaker is depicted in angular alignment with the laser trigger point. The directivities towards the microphones are calculated for the sources in the reduced source map that correspond to the speakers in the line array. The discrete directivities are linearly interpolated in between the microphone positions and plotted over the convex hull of the array. The blue dotted lines indicate the theoretical maximum of the speaker array. Figure 4 shows test case 1 with a non-rotating speaker array. There is no delay between the speakers, which results in a 0 degree main lobe direction. Given that the center of the speaker array is at a vertical position of −225 mm, the source locations are found at the correct positions, and the directivity distribution also has its maximum at microphones at this height. Since there are no aerodynamic and motor noises, the solution of the directivity distribution is sparse. Test case 2 also uses the same signal on every speaker, but the speakers are rotated at 205 revolutions per minute. The resulting source distributions and directions are shown in Figure 5. Here, the source localization plot shows artifacts from the virtual rotation, but the strongest source is still located at the correct position. The directivity has its maximum at the same height, but it spreads over more microphone positions.
Test cases 3-5 are shown in Figures 6-8. In these test cases, the main lobe of the speaker array is steered by delaying the speakers. All three cases show the source at the correct location with some minor errors due to the virtual rotation. The directivities differ according to the main radiation direction of the speaker array. For the −20 degree direction, the maximum of the source directivity is close to the theoretical maximum of −585 mm at the array plane. For the positive 20 degree, directivity appears on the top half of the microphone plane, opposite to the speaker array. The −40 degree shows the maximum at the microphones at the top but also displays the first side lobe of the speaker array at the bottom of the array plane.  Figure 9 shows the simulated case for comparison. The source distribution also shows some minor errors along the rotation axis due to the VRA, but source location and directivity are at the expected locations. Compared to test case 2, the ideal conditions in the simulation result in a sparser source distribution and directivity. This is mainly caused by the lack of noise in the simulation and the usage of ideal monopoles without casing instead of real loudspeakers.

Conclusions
This contribution demonstrates a methodology for the calculation of sound source distribution and discrete directivities for rotating sources using a planar microphone array. The motion compensation of the rotating sources is achieved using the virtual rotating array method for a spiral array geometry. A combination of two inverse methods is able to find the correct source positions for standing and moving sources. For the calculation of the source positions, CMF with regularization is used to minimize the influence of noise for the measurement data. The experimental results show that for measured data, as well as synthetic data, the method is able to calculate the correct source positions. The directivities of the loudspeaker array for the four different radiation angles are clearly distinguishable. The simulated data produce a sparse directivity result due to the lack of noise. The usage of a smoothing factor for the directivities might be helpful in these cases. For the measurement data, the directivity patterns are very plausible even without smoothing.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

AC
Alternating current CMF Covariance matrix fitting CSM Cross-spectral matrix RMS Root mean square ROSI Rotating source identifier SEM Spectral estimation method SODIX Source directivity modeling in the cross-spectral matrix VRA Virtual rotating array