Separation of the Sound Power Spectrum of Multiple Sources by Three-Dimensional Sound Intensity Decomposition

The identification and separation of sources are the prerequisite of industrial noise control. Industrial machinery usually contains multiple noise sources sharing same-frequency components. There are usually multiple noise sources in mechanical equipment, and there are few effective methods available to separate the spectrum intensity of each sound source. This study tries to solve the problem by the radiation relationship between three-dimensional sound intensity vectors and the power of the sources. When the positions of the probe and the sound source are determined, the sound power of the sound source at each frequency can be solved by the particle swarm optimization algorithm. The solution results at each frequency are combined to obtain the sound power spectrum of each sound source. The proposed method is first verified by a simulation on two point sources. The experiment is carried out on a fault simulation test bed in an ordinary laboratory; we used three three-dimensional sound intensity probes to form a line array and conducted spectrum separation of the nine main noise sources. The sound intensity on the main frequency band of each sound source was close to the result of the near-field measurement of the one-dimensional sound intensity probe. The proposed spectral separation method of the sound power of multiple sound sources provides a new method for accurate noise identification in industrial environments.


Introduction
The operation of mechanical equipment is often accompanied by strong noise, which contains a wealth of information about the status of equipment parts. If each sound source inside the machine can be separated, it is of great significance for monitoring the equipment state through noise [1]. In rotating machinery, the sound sources are mainly the rotating parts. Due to the transmission relationship, the noise signal usually contains same-frequency components. Therefore, the noise separation method must have the ability of accurate intensity separation for sounds of the same frequency.
In the field of signal processing, there are many ways to separate signals. The amplitude of each signal source obtained by blind signal separation is uncertain-that is, it is difficult to obtain the accurate strength of each source. The signal sources separated by this method need to be independent of each other, and sources with same-frequency components cannot be separated [2]. Empirical mode decomposition is based on the time-scale characteristics of the data [3], and methods such as wavelet decomposition are based on the frequency characteristics of the signal [4]. The components obtained by those methods are not directly related to the actual signal sources. Thus, these methods are difficult to apply to the separation of multiple sound sources in rotating machinery.
Among the sound source identification methods, near-field acoustic holography uses a microphone array to accurately reconstruct the sound pressure, particle velocity, and sound

Basic Principles of the Three-Dimensional Sound Intensity Measurement
The four microphones of the three-dimensional sound intensity probe were placed at the four vertices of a cube [18], as shown in Figure 1. First, we calculated the sound intensity components from the center of the cube to the vertices where the four microphones were located according to the dual-microphone cross-spectrum method. Then, we performed vector calculations according to the geometric characteristics of the cube and obtained the sound intensity components in the X, Y, and Z directions of the center point of the cube [19,20]. obtained the sound intensity components in the X, Y, and Z directions of the cente of the cube [19,20]. In Figure 1, the diagonal length of the cube surface is 1 d , and the distance fro vertex to the center O is a , .
Among the 4 microphones, the conne any 3 microphones can form an equilateral triangle. According to the sound intensity measurement of the cross-power spectrum m for a one-dimensional sound intensity probe composed of two microphones, the intensity at the center is:

 
In the formula, ( ) are the Fourier transform of the sound p ( ) P t and the particle velocity ( ) U t   at the center point, respectively; Re is the re of the complex number, and Im is the imaginary part. The superscript "*" indica conjugate of the complex number. 12 G is the cross-power spectrum on the soun sure of the two microphones, and r Δ is the distance between the microphones. As shown in Figure 1, the vector synthesis of the four particle velocities from t microphones to the center O corresponds to the total particle vibration velocity at t ter O. By decomposing the particle vibration velocity into 3 coordinate axes, we tain: In Figure 1, the diagonal length of the cube surface is d 1 , and the distance from each vertex to the center O is a, a = ( Among the 4 microphones, the connection of any 3 microphones can form an equilateral triangle. → U x (t), → U y (t), and → U Z (t) are the particle velocities in the X, Y, and Z axes of the three-dimensional coordinate system, respectively. Supposing the sound pressure measured by the four microphones is p 1 (t), p 2 (t), p 3 (t), and p 4 (t) and their Fourier transform is P 1 (ω), P 2 (ω), P 3 (ω), and P 4 (ω), abbreviated as P 1 , P 2 , P 3 , and P 4 . Then, the sound pressure P O at the center point O is approximately equal to the average of the sound pressure measured by the four microphones, namely: According to the sound intensity measurement of the cross-power spectrum method, for a one-dimensional sound intensity probe composed of two microphones, the sound intensity at the center is: In the formula, P(ω) and → U(ω) are the Fourier transform of the sound pressure P(t) and the particle velocity → U(t) at the center point, respectively; Re is the real part of the complex number, and Im is the imaginary part. The superscript "*" indicates the conjugate of the complex number. G 12 is the cross-power spectrum on the sound pressure of the two microphones, and ∆r is the distance between the microphones.
As shown in Figure 1, the vector synthesis of the four particle velocities from the four microphones to the center O corresponds to the total particle vibration velocity at the center O. By decomposing the particle vibration velocity into 3 coordinate axes, we can obtain: Taking Equation (2) as a reference, we can obtain the sound intensity vector in the X direction at the center point O as: In Equation (6), → U x (ω) is the particle velocity in the X direction and G ij is the crosspower spectrum of the i-th microphone and the j-th microphone signals.
Similarly, the sound intensity vector in the Y direction at the center point O is: The sound intensity vector in the Z direction at the center point O is: According to Equations (6)-(8), the sound intensity component at the center point O can be calculated from the sound pressure on the 4 microphones. The total sound intensity vector is composed of components in the X, Y, and Z directions. If there are multiple sound sources, the measured sound intensity vectors are the synthesis of multiple sound sources, and the intensity of each sound source cannot be obtained. However, if the number of three-dimensional sound intensity probes is increased and the sound field is sampled at multiple points, it is possible to completely determine the sound power spectrum of the sound source. Lu Yi used this method to identify the position and intensity of two sound sources by three three-dimensional sound intensity probes and showed that this method has high accuracy.

The Separation of the Sound Power Spectrum of Multiple Sound Sources
For N point sources, the coordinates in the spherical coordinate system are → P Sn = (ρ Sn , α Sn , θ Sn ), n = 1, 2, . . . , N, and the corresponding sound power is W n , n = 1, 2, . . . , N. M three-dimensional sound intensity probes are used for the measurement, and the probe coordinates are → P Tm = (ρ Tm , α Tm , θ Tm ), m= 1, 2, . . . , M, respectively. Then, the path vector from the sound source n to the probe m is → P Tm − → P Sn . According to the propagation law of the monopole sound source in the free field, without considering the actual environmental factors such as scattering, the sound intensity generated by the sound source n on the probe m is: The total sound intensity measured on each probe is the superposition of the sound intensities produced by each of the N sound sources, which is: When the positions of sound sources and probes are known, Equation (10) constitutes a linear equation system for measuring sound intensity and sound power, since the sound intensity is vectoral and the m vectors can be decomposed into a spatial rectangular coordinate system to obtain 3M equations. In theory, as long as the number of probes is no less than 1/3 of the number of sound sources, the equation set can be solved directly to obtain the sound power of each sound source on each frequency (i.e., W n (ω)). However, in practice, due to the fact that the sound sources are not ideal monopole sources as well as directionality of intensity radiation, sound field environment, and noise interference, it is more feasible to solve the optimal solution of the equations by numerical methods.
Particle swarm optimization (PSO) is used to solve the 3M equations decomposed from Equation (10). The authors of the particle swarm optimization algorithm are Kennedy and Eberhart [21]. The characteristics of this algorithm are random search and parallel optimization. The algorithm can be well adapted to engineering applications [22]. The basic concept of PSO comes from the study of the foraging behavior of bird flocks. It uses particles to imitate the social behavior of individual birds when searching for food and searches for the global optimum in parallel in a multi-dimensional space [23]. In multi-dimensional space, each particle can be regarded as a search individual, the current position of the particle is the candidate solution of the corresponding optimization problem, and the flight process of the particle is the search process of the search individual. The optimal individual extreme value in the particle swarm is taken as the current global optimal value. In the iterative process, the particles update their velocity and position through individual optimal values and global optimal values. The iterative formula of particle velocity and position in particle swarm algorithm is: where is a solution of the optimization problem corresponding to the position of the i-th particle in the particle swarm; is the optimal position at all distances experienced by the i-th particle; P g = (p 1 g , p 2 g , . . . , p D g ) is the optimal position at all distances experienced by all particles; is the search speed of particle i; ω is the non-negative inertia factor, c 1 and c 2 are nonnegative acceleration constants; r 1 and r 2 are random numbers in the range of 0 to 1; α is the constraint factor to control the weight of speed.
The process flow of the PSO algorithm is as follows: (1) Initialize the particle swarm, including the population size (set to 100), inertia weight ω (set to 0.1), acceleration constant c 1 (set to 2.2), c 2 (set to 2.1), and the random position and the search speed of each particle (set the maximum speed to 3 m/s). (2) Calculate the objective function that is the fitness function.
(3) Update the individual optimal value of the particle and the global optimal value of the group according to the fitness function value and update the optimal position and optimal speed of the particle according to Equations (11) and (12). (4) Determine whether the maximum number of iterations (set to 1000) or the global optimal position is reached and whether the minimum limit is satisfied. If satisfied, end the iteration; otherwise, repeat steps 2-4.
The same parameters were used both in the experiment and in the simulation. Since both the sound intensity I Tm and the sound power W n are the function of the frequency, the particle swarm algorithm can be used to solve Equation (10) at each frequency point of the spectrum to obtain the sound power at each frequency point. It is equivalent to the power spectral density of sound power, and the frequency spectrum separation of each sound source can be realized.

Simulation Analysis
A simulation experiment of two sound sources measured with one three-dimensional sound intensity probe was established. The spherical coordinates of the three-dimensional intensity probe were at 500 Hz was 0.00605 W, the sound power at 1000 Hz was 0.02 W, and the sound power at 2000 Hz was 0.0128 W. The sound power of the three frequencies in the sound source 2 was 0.08, 0.0072, and 0.03125 W, respectively.
We substituted the signals on the four microphones into Equations (6)- (8) to calculate the sound intensity values in the X, Y, and Z directions of the probe, as shown in Table 1. We then substituted the above position and sound intensity into Equation (10) to solve the sound power of the sound source. According to the calculation results, the sound power spectra of the two sound sources were separated, as shown in Figure 2. It can be seen that the separated sound power spectrum is consistent with the theoretical sound power spectrum of the simulated sound source, and there is no obvious error. . Both sound sources contained three frequencies of 500, 1000, and 2000 Hz. In sound source 1, the sound power at 500 Hz was 0.00605 W, the sound power at 1000 Hz was 0.02 W, and the sound power at 2000 Hz was 0.0128 W. The sound power of the three frequencies in the sound source 2 was 0.08, 0.0072, and 0.03125 W, respectively.
We substituted the signals on the four microphones into Equations (6)- (8) to calculate the sound intensity values in the X, Y, and Z directions of the probe, as shown in Table 1. We then substituted the above position and sound intensity into Equation (10) to solve the sound power of the sound source. According to the calculation results, the sound power spectra of the two sound sources were separated, as shown in Figure 2. It can be seen that the separated sound power spectrum is consistent with the theoretical sound power spectrum of the simulated sound source, and there is no obvious error.

Experimental Verification
In order to verify the feasibility of the three-dimensional sound intensity vector decomposition method for spectrum separation in the actual environment, experimental research was carried out on the QPZZ-II fault simulation test bench in an ordinary laboratory environment. The experimental bench was located at the corner of the laboratory, and the distances to the two walls were 940 and 610 mm. The nearby walls were treated with simple sound absorption. The speed of the motor was 888.75 r/min and the shaft frequency was 888.75/60 = 14.8125 Hz.
As shown in Figure 3, the sound sources on the test bench included motors, pulleys, couplings, and various support bearings. The vertical distance between the probe and the transmission shaft plane was 500 mm. Since the distance between the center of the nine sound sources and the test point was more than two times the maximum geometric size of the sound sources, the nine sound sources were all regarded as point sound sources. The spherical coordinate system was set up with the axial position of bearing seat I (S1) as the origin. The specific information of the sound sources is shown in Table 2.

Experimental Verification
In order to verify the feasibility of the three-dimensional sound intensity vector decomposition method for spectrum separation in the actual environment, experimental research was carried out on the QPZZ-II fault simulation test bench in an ordinary laboratory environment. The experimental bench was located at the corner of the laboratory, and the distances to the two walls were 940 and 610 mm. The nearby walls were treated with simple sound absorption. The speed of the motor was 888.75 r/min and the shaft frequency was 888.75/60 = 14.8125 Hz.
As shown in Figure 3, the sound sources on the test bench included motors, pulleys, couplings, and various support bearings. The vertical distance between the probe and the transmission shaft plane was 500 mm. Since the distance between the center of the nine sound sources and the test point was more than two times the maximum geometric size of the sound sources, the nine sound sources were all regarded as point sound sources. The spherical coordinate system was set up with the axial position of bearing seat I (S1) as the origin. The specific information of the sound sources is shown in Table 2.    In order to obtain the true strength of each sound source, the one-dimension intensity probe was used to collect the intensity of each source at a distance of 30 m the centers of each sound source.
Since the one-dimensional sound intensity probe was close to the sound so pointed to the sound source, the influence of other sound sources can be ignored ing to the characteristics of sound intensity. Based on the point sound source hy the sound power spectrum of each sound source can be obtained referring to Equ This sound power spectrum is taken as the ground truth to verify the proposed m Since there were nine sound sources, three three-dimensional sound intensit were used for testing. Four MPA416 array microphones were installed on each dynamic signal instrument was used to collect the microphone signals, and the s frequency was 25.6 kHz. The vertical distance between the probe and the tran shaft plane was 500 mm. Probe 1 was located above the loading device (sound so and the spherical coordinates were 1

S1
The axis of bearing seat I The axis of the loading device The axis of bearing seat II The axis of the rigid coupling The axis of the midpoint of the bearing support The axis of the pulley I The axis of the pulleyII The rotor in the middle of the motor The fan at the rear of the motor In order to obtain the true strength of each sound source, the one-dimensional sound intensity probe was used to collect the intensity of each source at a distance of 30 mm from the centers of each sound source.
Since the one-dimensional sound intensity probe was close to the sound source and pointed to the sound source, the influence of other sound sources can be ignored according to the characteristics of sound intensity. Based on the point sound source hypothesis, the sound power spectrum of each sound source can be obtained referring to Equation (9). This sound power spectrum is taken as the ground truth to verify the proposed method.
Since there were nine sound sources, three three-dimensional sound intensity probes were used for testing. Four MPA416 array microphones were installed on each probe. A dynamic signal instrument was used to collect the microphone signals, and the sampling frequency was 25.6 kHz. The vertical distance between the probe and the transmission shaft plane was 500 mm. Probe 1 was located above the loading device (sound source S6), and the spherical coordinates were → P T1 = (93, 90 • , 57 • ). Probe 2 was located above bearing seat II (sound source S3), and the coordinates were → P T2 = (67, 90 • , 41 • ). Probe 3 was located above pulley I (sound source S6), with the coordinates → P T3 = (51, 90 • , 9.0 • ). The sound intensity spectra of the main frequency bands measured by three threedimensional sound intensity probes are shown in Figure 4. We selected the frequency point with a sound power value greater than 1 × 10 −13 W for sound power spectrum separation. The reasons are that firstly, the sound power value at other frequency points is too small, and the contribution value to the sound power spectrum is not large; secondly, the amount of calculations can be reduced. small, and the contribution value to the sound power spectrum is not large; secondly, the amount of calculations can be reduced.
(b) Probe 2.  The positions of the sound sources, the positions of the probes, and the sound intensity were substituted into Equation (10) as known quantities, and the PSO was used to solve the sound power of the nine sound sources at each frequency point. The separated spectra were compared with the true power spectra measured by the one-dimensional sound intensity probe, as shown in Figure 5. (a) Sound power spectra of sound source 1 The positions of the sound sources, the positions of the probes, and the sound intensity were substituted into Equation (10) as known quantities, and the PSO was used to solve the sound power of the nine sound sources at each frequency point. The separated spectra were compared with the true power spectra measured by the one-dimensional sound intensity probe, as shown in Figure 5.  The positions of the sound sources, the positions of the probes, and the sound intensity were substituted into Equation (10) as known quantities, and the PSO was used to solve the sound power of the nine sound sources at each frequency point. The separated spectra were compared with the true power spectra measured by the one-dimensional sound intensity probe, as shown in Figure 5.   (g) Sound power spectra of sound source 7.  The sound power spectra of the nine sound sources calculated by the sound intensity vector decomposition method are very close to the true values, and the errors are relatively small at the frequencies with higher intensity.
The dominating frequency in the sound sources S1-S6 is 472 Hz, and the recognition errors are all small. The specifications of the two pulleys are the same, and the number of teeth is 32. The meshing frequency of the toothed belt is equal to the product of the shaft frequency and the number of teeth, and the meshing frequency of the ruler belt is 474 Hz. Here, 472 Hz is close to the meshing frequency of the toothed belt. The six sound sources The sound power spectra of the nine sound sources calculated by the sound intensity vector decomposition method are very close to the true values, and the errors are relatively small at the frequencies with higher intensity.
The dominating frequency in the sound sources S1-S6 is 472 Hz, and the recognition errors are all small. The specifications of the two pulleys are the same, and the number of teeth is 32. The meshing frequency of the toothed belt is equal to the product of the shaft frequency and the number of teeth, and the meshing frequency of the ruler belt is 474 Hz. Here, 472 Hz is close to the meshing frequency of the toothed belt. The six sound sources are distributed on two rigidly connected slender shafts, thus sharing the same peak frequency. The power on S1, S4, and S6 is greater than on the other three sources. S1 is the test bearing of the test bench, which is often disassembled and assembled, and the bearing clearance is large. S4 is the coupling position, and there may be a matching problem. S6 is pulley I, which is directly excited by the meshing of the belt and the gear teeth.
For the sound source S7, the dominating frequency is 475 Hz. S7 is the pulley on the driving side. Compared with the pulley S6 on the load side, the peak frequency shifted slightly.
The dominating frequency of S8 is 447 Hz. S8 is at the center of the motor, and the sound intensity is the highest among all sources. The frequency spectrum structure is more complicated. The sound source S9 is a motor fan with five blades and its peak frequency is 519 Hz, which is seven times the blade-passing frequency of 74 Hz.
It can be seen that the three-dimensional sound intensity vector decomposition method can recover the intensity of the main frequency range of each sound source accurately, and the result is equivalent to that of one-dimensional sound intensity measured at close distance. The recovery ratio Rr is used to quantify the accuracy. Rr is the proportion of the separated sound power W s (ω) and the true sound power W t (ω), as: To omit the background noise, W s (ω) and W t (ω) were selected from the sound power spectra with values greater than 1 × 10 −12 W.
where W s (ω) is the sound power spectrum separated by the proposed method, and W t (ω) is the true sound power spectrum calculated from the near-field sound intensity measurement. The recovery ratio Rr values of the nine sources are shown in Figure 6. The closer this ratio is to 1, the closer the separation sound power value is to the true value and the smaller the separation error. The highest value is 0.9 for source 5, while the lowest is 0.7 for source 2. are distributed on two rigidly connected slender shafts, thus sharing the same peak frequency. The power on S1, S4, and S6 is greater than on the other three sources. S1 is the test bearing of the test bench, which is often disassembled and assembled, and the bearing clearance is large. S4 is the coupling position, and there may be a matching problem. S6 is pulley I, which is directly excited by the meshing of the belt and the gear teeth. For the sound source S7, the dominating frequency is 475 Hz. S7 is the pulley on the driving side. Compared with the pulley S6 on the load side, the peak frequency shifted slightly.
The dominating frequency of S8 is 447 Hz. S8 is at the center of the motor, and the sound intensity is the highest among all sources. The frequency spectrum structure is more complicated. The sound source S9 is a motor fan with five blades and its peak frequency is 519 Hz, which is seven times the blade-passing frequency of 74 Hz.
It can be seen that the three-dimensional sound intensity vector decomposition method can recover the intensity of the main frequency range of each sound source accurately, and the result is equivalent to that of one-dimensional sound intensity measured at close distance. The recovery ratio Rr is used to quantify the accuracy. Rr is the proportion of the separated sound power To omit the background noise, where ( ) s W ω is the sound power spectrum separated by the proposed method, and ( ) t W ω is the true sound power spectrum calculated from the near-field sound intensity measurement. The recovery ratio Rr values of the nine sources are shown in Figure 6. The closer this ratio is to 1, the closer the separation sound power value is to the true value and the smaller the separation error. The highest value is 0.9 for source 5, while the lowest is 0.7 for source 2.  Figure 6. The recovery ratios of the nine sound sources in the experiment.

Conclusions
Based on the three-dimensional sound intensity measurement method, this paper proposes the sound power spectrum separation method for multiple sound sources. Given the positions of the sound sources, by solving the equations of the three-dimensional sound intensity on the sound power of the sources, the sound power of each sound source at each frequency is calculated. Thus, the sound power spectrum of each sound source is obtained. Taking the traditional near-field one-dimensional sound intensity measurement as the reference, the main frequency bands of multiple sound sources on the test bench were separated in an ordinary laboratory, and the accuracy of the three-dimensional sound intensity vector decomposition spectrum separation method was verified.
Compared with other sound source identification methods, this method can separate the sound power spectra of multiple noise sources over a longer distance, providing a new method for noise monitoring and diagnosis of rotating machinery. However, there are some shortcomings. First, it is assumed that the sound sources are ideal point sound sources, but in fact, the sound source may have a certain volume or have multiple types, which will also affect the coordinates of the sound sources to a certain extent. Furthermore, there are interfering sound sources in an ordinary laboratory, but we ignored the noise interference from the external environment. In addition, when there are too many sound sources or too many frequency components, too much calculation may reduce the calculation accuracy. In the future, further research will be conducted from the above-mentioned aspects.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.