Real-time separation of non-stationary sound fields with pressure and particle acceleration measurements Chuan-Xing Bi,a) Lin Geng, and Xiao-Zheng Zhang Institute of Sound and Vibration Research, Hefei University of Technology, 193 Tunxi Road, Hefei 230009, People’s Republic of China
(Received 24 September 2013; revised 28 February 2014; accepted 28 April 2014) To extract the desired non-stationary sound field generated by a target source in the presence of disturbing sources, a real-time sound field separation method with pressure and particle acceleration measurements is proposed. In this method, the pressure and particle acceleration signals at a time instant are first measured on one measurement plane, where the particle acceleration is obtained by the finite difference approximation with the aid of an auxiliary measurement plane; then, the desired pressure signal generated by the target source at the same time instant can be extracted in a timely manner, by a simple superposition of the measured pressure and the convolution between the measured particle acceleration and the derived impulse response function. Thereby, the proposed method possesses a significant feature of real-time separation of non-stationary sound fields, which provides the potential to in situ analyze the radiation characteristics of a non-stationary source. The proposed method was examined through numerical simulation and experiment. Results demonstrated that the proposed method can not only extract the desired time-evolving pressure signal generated by the target source at any space point, but can also obtain the desired spatial distribution of the pressure field generated by the target source at any time instant. C 2014 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4875576] V PACS number(s): 43.60.Gk, 43.20.Px, 43.60.Pt [MS]
I. INTRODUCTION
The sound field separation method (FSM) has been proven to be an effective technique for extracting the spatial acoustic characteristics of a target source in the presence of disturbing sources. Currently, several different methods, such as the spatial Fourier transform method,1–6 the spherical wave expansion method,7–9 the statistically optimized nearfield acoustic holography,10–13 the boundary element method,14–16 and the equivalent source method17–19 have been developed to realize this separation for distinguishing the sources on both sides of the measurement surface and suppressing the influence of disturbing sources. These FSMs require measurements of pressures or particle velocities on two measurement surfaces, or measurements of pressure and particle velocity on one measurement surface, and they all have the ability to extract the sound field generated by the target source in the frequency domain. However, for a nonstationary source whose statistical properties fluctuate in time, it is more reasonable to separate the sound field at a time instant for further analyzing the acoustic characteristics related to time. Recently, on the basis of the impulse response approach,20,21 a method for separating the non-stationary sound fields by measuring the pressures on two closely spaced measurement planes was developed by Zhang et al.22 In that method, the pressure signal generated by the target source at each time step is solved via a deconvolution process. The deconvolution process implies a complicated a)
Author to whom correspondence should be addressed. Electronic mail: [email protected]
3474
J. Acoust. Soc. Am. 135 (6), June 2014
Pages: 3474–3482
iterative calculation, i.e., the pressure signal generated by the target source at the present time step has to be calculated in the need of the pressure signals solved at all the previous time steps. Furthermore, the deconvolution process is an illposed inverse problem. To obtain an appropriate solution, the singular value decomposition (SVD) and the Tikhonov regularization23 have to be used, which leads to a large computational workload. In the present paper, an alternative method, namely, real-time sound field separation method (RT-FSM), is proposed to separate the non-stationary sound fields as well. In RT-FSM, the desired pressure signal generated by the target source at a specific time instant can be extracted by a simple superposition of the measured pressure at the same time instant and the convolution between the measured particle acceleration and the corresponding impulse response function. Here, the particle acceleration is obtained by the finite difference approximation with the pressures measured on two closely spaced measurement planes. Compared to the deconvolution process involved in Zhang’s method, the method proposed in the present paper only contains a simple forward solving process, i.e., the pressure signal generated by the target source at the present time step can be extracted with no need of any pressure signal solved at the previous time steps. And more importantly, the forward solving process in RT-FSM avoids the use of the SVD and the regularization, which reduces the computational workload significantly. Therefore, the RT-FSM provides a better feature of real-time separation than Zhang’s method, which gives the potential to in situ analyze the radiation characteristics of a non-stationary source. The paper is organized as follows. Section II presents a theoretical description of RT-FSM. In Sec. III, to examine
0001-4966/2014/135(6)/3474/9/$30.00
C 2014 Acoustical Society of America V
the performance of RT-FSM in both time and space domains, a numerical simulation with three monopole sources driven by different non-stationary signals is investigated. An experiment involving three speakers driven by a Morlet wavelet signal is carried out in Sec. IV. Finally, conclusions are summarized in Sec. V. II. REAL-TIME SOUND FIELD SEPARATION METHOD
Figure 1 illustrates the geometry of RT-FSM in the Cartesian coordinate system oðx; y; zÞ. Because of the joint contributions from the target source Mo and the disturbing source Md , the pressure and particle acceleration on the measurement plane H are both mixed. The purpose of the present paper is to remove the influence of the disturbing source Md and extract the time-evolving pressure signal generated by the target source Mo from the mixed one. According to the superposition principle of waves, the mixed pressure signal pðx; y; zH ; tÞ and the mixed particle acceleration signal aðx; y; zH ; tÞ at the measurement point ðx; y; zH Þ at the time t can be expressed, respectively, as pðx; y; zH ; tÞ ¼ po ðx; y; zH ; tÞ þ pd ðx; y; zH ; tÞ;
(1)
aðx; y; zH ; tÞ ¼ ao ðx; y; zH ; tÞ ad ðx; y; zH ; tÞ;
(2)
where the subscripts “o” and “d” denote the sound fields generated by the target source and the disturbing source, respectively. By applying a two-dimensional Fourier transform with respect to x and y to Eqs. (1) and (2), it yields
Pd ðkx ; ky ; zH ; tÞ ¼ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ;
(6)
where the asterisk denotes the convolution of two time functions. The expression of the impulse response function hðkx ; ky ; 0; tÞ is derived in the Appendix and given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðkx ; ky ; 0; tÞ ¼ q0 cJ0 (7) c2 ðkx2 þ ky2 Þt HðtÞ; where q0 is the medium density, c is the sound velocity, J0 is the Bessel function of the first kind of order zero, and HðtÞ is the Heaviside function. By replacing the term Pd ðkx ; ky ; zH ; tÞ in Eq. (3) with the right part of Eq. (6), Eq. (3) can be rewritten as Pðkx ; ky ; zH ; tÞ ¼ Po ðkx ; ky ; zH ; tÞ þ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ:
(8)
By convolving with the impulse response function hðkx ; ky ; 0; tÞ on both sides of Eq. (4), it yields Aðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ ¼ Ao ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ :
(9)
In accordance with Eq. (5), Eq. (9) can be further expressed as Aðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ ¼ Po ðkx ; ky ; zH ; tÞ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ:
Pðkx ; ky ; zH ; tÞ ¼ Po ðkx ; ky ; zH ; tÞ þ Pd ðkx ; ky ; zH ; tÞ;
(3)
(10)
Aðkx ; ky ; zH ; tÞ ¼ Ao ðkx ; ky ; zH ; tÞ Ad ðkx ; ky ; zH ; tÞ;
(4)
By combining Eq. (8) with Eq. (10), the time–wave number spectrum of the pressure Po ðkx ; ky ; zH ; tÞ can be finally obtained by
where kx and ky represent the wave number components in x and y directions, respectively. Referring to the derivation in Appendix, the time–wave number spectrum of the pressure can be related to that of the particle acceleration produced by the same source with the following formulas: Po ðkx ; ky ; zH ; tÞ ¼ Ao ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ;
(5)
Po ðkx ;ky ;zH ;tÞ ¼ 0:5½Pðkx ;ky ;zH ;tÞþ Aðkx ;ky ;zH ;tÞ hðkx ;ky ;0;tÞ:
(11)
Equation (11) indicates that the pressure spectrum Po ðkx ; ky ; zH ; tÞ generated by the target source on the measurement plane H can be extracted from the mixed pressure spectrum Pðkx ; ky ; zH ; tÞ with the aid of the mixed particle acceleration spectrum Aðkx ; ky ; zH ; tÞ. In Eq. (11), the convolution is continuous in the time domain. To meet the requirement of numerical computation, the time t should be discretized as tn ¼ ðn 1ÞDt;
(12)
where n ¼ 1; 2; …; N, N is the total number of time steps, and Dt is the time step. Accordingly, Eq. (11) can be expressed in a discrete form as follows: Po ðkx ;ky ;zH ;tn Þ ¼0:5 Pðkx ;ky ;zH ;tn Þ n X Aðkx ;ky ;zH ;ti Þhðkx ;ky ;0;tniþ1 Þ : þ FIG. 1. (Color online) Geometry of RT-FSM in the Cartesian coordinate system oðx; y; zÞ. J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
i¼1
(13) Bi et al.: Real-time sound field separation
3475
From Eq. (13), it is noted that once the mixed pressure spectrum at the time instant tn and the mixed particle acceleration spectrums from the time instant t1 to the time instant tn are known, the pressure spectrum Po ðkx ; ky ; zH ; tn Þ at the time instant tn can be obtained. In fact, the mixed particle acceleration spectrums from the time instant t1 to the time instant tn1 have been achieved before performing the measurement at the time instant tn , and therefore, the pressure spectrum Po ðkx ; ky ; zH ; tn Þ can be extracted in real time only by measuring the mixed pressure and particle acceleration at the time instant tn . Since the pressure field in the space domain is desired, the two-dimensional inverse Fourier transform with respect to kx and ky is then applied to the pressure spectrum Po ðkx ; ky ; zH ; tn Þ to obtain the pressure po ðx; y; zH ; tn Þ. The procedure of RT-FSM can be described clearly with a flow chart as shown in Fig. 2. In RT-FSM, the particle acceleration should be measured as one of the inputs to start the separation. Unfortunately, it is difficult to measure the particle acceleration directly. In the following, an alternative way to acquire the particle acceleration is provided, which is based on the Euler’s equation given by aðx; y; zH ; tÞ ¼
1 @pðx; y; zH ; tÞ : q0 @z
(14)
Equation (14) indicates that the particle acceleration is related to the partial derivative of the pressure. Since the partial derivative is often approximated by the finite difference, an auxiliary measurement plane H1 is placed near the measurement plane H and the partial derivative in Eq. (14) is approximated by the finite difference of the pressures on the planes H and H1 as @pðx; y; zH ; tÞ pðx; y; zH1 ; tÞ pðx; y; zH ; tÞ ; @z zH1 zH
(15)
where pðx; y; zH1 ; tÞ is the mixed pressure measured on the auxiliary measurement plane H1.
Substituting Eq. (15) into Eq. (14), it finally yields the approximate expression of the particle acceleration as aðx; y; zH ; tÞ
pðx; y; zH1 ; tÞ pðx; y; zH ; tÞ : q0 ðzH1 zH Þ
(16)
III. NUMERICAL SIMULATION
A numerical simulation with three monopole sources was investigated to verify the feasibility of the proposed method. In the simulation, the monopole sources S1 and S3 were fixed on one side of the measurement plane as the target sources and the monopole source S2 was put on the other side as the disturbing source. The geometry of sources and measurement planes is shown in Fig. 3. The source S1 located at (0.1 m, 0.45 m, 0.145 m) generated a non-stationary signal with a linear frequency modulation in the [200, 1800] Hz band and a Gaussian amplitude modulation. The source S2 located at (0.4 m, 0.45 m, 0.14 m) generated a Morlet wavelet signal defined by sðtÞ ¼ cosð2pf0 tÞet
2
=2
;
(17)
where f0 ¼ 344 Hz. The source S3 located at (0.55 m, 0.35 m, 0.14 m) generated a non-stationary signal with a linear frequency modulation in the [150, 2400] Hz band and a Gaussian amplitude modulation. The measurement plane H was located at z ¼ 0 m in the Cartesian coordinate system oðx; y; zÞ, and the auxiliary measurement plane H1 located 0.02 m away from the measurement plane H was used to acquire the particle acceleration. These two planes both provided 1515 measurement points and gave the grid spacing of 0.05 m in both x and y directions. The theoretical pressure at each measurement point was calculated by the following equation: pðR; tÞ ¼
sðt R=cÞ ; 4pR
(18)
FIG. 2. (Color online) Flow chart of RT-FSM. Fxy and F1 xy denote twodimensional Fourier transform and the two-dimensional inverse Fourier transform with respect to x and y, respectively.
3476
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
FIG. 3. (Color online) Geometry of the measurement plane H, the auxiliary measurement plane H1, and the locations of three monopole sources in the Cartesian coordinate system oðx; y; zÞ. Three space points R1 , R2 , and R3 are chosen and marked with the symbol “þ.” The points R1 and R3 are facing the target sources S1 and S3, respectively, and the point R2 is facing the disturbing source S2.
where sðtÞ denotes the input signal, and R denotes the distance between the measurement point and the monopole source. The signals were sampled at a frequency of fe ¼25.6 kHz, providing 512 sampling points. A Gaussian white noise with a signal-to-noise ratio of 20 dB was added to the measured signals to simulate the practical measurement. To lessen the wrap-around errors due to the use of discrete twodimensional spatial Fourier transforms associated with the finite size of measurement planes, the 1515 pressure and particle acceleration matrices on the measurement plane H were extended to 2929 matrices by zero-padding. To assess the relevance of the proposed method in the time domain, three space points on the measurement plane H, identified in Fig. 3, were chosen and their positions were R1 (0.1 m, 0.45 m, 0 m), R2 (0.4 m, 0.45 m, 0 m), and R3 (0.55 m, 0.35 m, 0 m). The points R1 and R3 were facing the target sources S1 and S3, respectively, and the point R2 was facing the disturbing source S2. The time-domain waveform
comparisons among the theoretical pressures calculated after physically removing the disturbing source (line with plus sign), the mixed pressures generated by the target sources and the disturbing source together (dotted line), and the extracted pressures of target sources (solid line) at these three space points are shown in Fig. 4. The comparison results between the lines with a plus sign and the dotted lines indicate that the theoretical pressure signals radiated by the target sources are affected by the disturbing source in different degrees at these three points. To remove the influence of the disturbing source and extract the pressure generated by the target sources, the RT-FSM was performed. The comparison results between the lines with a plus sign and the solid lines indicate that the extracted pressures using the proposed method match well with the theoretical pressures radiated by the target sources. To estimate the consistency between the theoretical pressure calculated after physically removing the disturbing
FIG. 4. (Color online) Time-domain waveform comparisons among the theoretical pressures calculated after physically removing the disturbing source (line with plus sign), the mixed pressures generated by the target sources and the disturbing source together (dotted line), and the extracted pressures of target sources (solid line) at the three space points (a) R1 , (b) R2 , and (c) R3 .
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
3477
FIG. 5. (Color online) Spatial maps of (a) Ep with a contour line at the value of 0.9 and (b) Ea with a contour line at the value of 0.1. The space points R1 , R2 , and R3 are marked with the symbol “þ.”
source and the extracted pressure of target sources quantitatively, two indicators related to the phase and the amplitude were defined, respectively, by24 hpm ðx; y; z; tÞpe ðx; y; z; tÞi ; Epðx; y; zÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hp2m ðx; y; z; tÞihp2e ðx; y; z; tÞi
(19)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hp2 ðx; y; z; tÞi hp2 ðx; y; z; tÞi m e pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Eaðx; y; zÞ ¼ ; hp2m ðx; y; z; tÞi
(20)
where pm ðx; y; z; tÞ is the theoretical pressure calculated after physically removing the disturbing source, pe ðx; y; z; tÞ is the
FIG. 6. (Color online) Spatial distributions of the theoretical pressure fields calculated after physically removing the disturbing source at (a) t1 ¼ 9.18 ms and (b) t2 ¼ 12.86 ms, the mixed pressure fields generated by the target sources and the disturbing source together at (c) t1 ¼ 9.18 ms and (d) t2 ¼ 12.86 ms, and the extracted pressure fields of target sources at (e) t1 ¼ 9.18 ms and (f) t2 ¼ 12.86 ms on the measurement plane H. The three monopole sources are marked with the symbol “þ.”
3478
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
FIG. 7. (Color online) Experimental setup: three speakers, a microphone array, the recording system and the trigger. The speakers S1 and S3 are used as the target sources, and the speaker S2 is used as the disturbing source.
extracted pressure of target sources, and h i denotes the time averaged value. Here, it is noted that when Ep is close to unity, there exists a high phase similarity, and when Ea is near zero, the amplitude difference is small. To quantify the accuracy of the extracted results at all points on the measurement plane H, the spatial maps of the indicators Ep and Ea calculated by Eqs. (19) and (20), respectively, are shown in Fig. 5. The values of Ep and Ea at a majority of space points are more than 0.9 and below 0.1, respectively, which indicates that the extracted phases and amplitudes agree with the theoretical ones very well. Here, the values of Ep at the points R1 , R2 , and R3 are 0.9701, 0.9254, and 0.9713, respectively, and the values of Ea are 0.0513, 0.1399, and 0.0376, respectively.
Similarly, for showing the extracted results in the space domain, two time instants, t1 ¼9.18 ms and t2 ¼12.86 ms, were chosen. Figures 6(a), 6(c), and 6(e) depict the theoretical pressure field calculated after physically removing the disturbing source, the mixed pressure field generated by the target sources and the disturbing source together, and the extracted pressure field of target sources at t1 ¼9.18 ms, respectively. Figures 6(b), 6(d), and 6(f) depict the same pressure fields but at t2 ¼12.86 ms. The comparisons between the theoretical pressure fields in Figs. 6(a) and 6(b) and the mixed pressure fields in Figs. 6(c) and 6(d) at these two time instants illustrate that the spatial distributions of theoretical pressure fields have been changed by the disturbing source. After performing the RT-FSM, it is obvious that the extracted pressure fields in Figs. 6(e) and 6(f) are almost the same as the theoretical ones in Figs. 6(a) and 6(b), which indicates that the proposed method seems effective to suppress the sound radiated by the disturbing source in the space domain. IV. EXPERIMENTAL STUDY
An experiment was carried out in a semi-anechoic chamber to further verify the proposed method. In the experiment, three speakers S1, S2, and S3 with the same sizes were used as the sources, among which the speakers S1 and S3 were the target sources and the speaker S2 was the disturbing source, as shown in Fig. 7. All speakers were connected with a laptop and all driven by the same Morlet wavelet signal as that used in the simulation case. In the experiment, an array with 15 microphones was moved 15 times on each measurement surface to measure the pressure signals at 1515 measurement points. To keep
FIG. 8. (Color online) Time-domain waveform comparisons among the pressures measured after physically removing the disturbing source (line with plus sign), the mixed pressures generated by the target sources and the disturbing source together (dotted line), and the extracted pressures of target sources (solid line) at the three space points (a) R1 , (b) R2 , and (c) R3 .
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
3479
FIG. 9. (Color online) Spatial maps of (a) Ep with a contour line at the value of 0.9 and (b) Ea with a contour line at the value of 0.1. The space points R1 , R2 , and R3 are marked with the symbol “þ.”
the measured signals the same at each measurement, the output of the source signal was set as the trigger to activate the acquisition system recording the data. In the experiment, the centers of the three speaker cones were positioned at the same locations of the three sources in the simulation case, and the locations of measurement planes and the grid spacing of measurement points were all the same
as those in the simulation case. The signals were sampled at a frequency of fe ¼25.6 kHz, providing 512 sampling points. The 1515 pressure and particle acceleration matrices on the measurement plane H were also extended to 2929 matrices by zero-padding to lessen the wrap-around errors. Just as the simulation case, the same three space points R1 , R2 , and R3 were selected to show the time-domain
FIG. 10. (Color online) Spatial distributions of the pressure fields measured after physically removing the disturbing source at (a) t1 ¼ 9.38 ms and (b) t2 ¼ 10.78 ms, the mixed pressure fields generated by the target sources and the disturbing source together at (c) t1 ¼ 9.38 ms and (d) t2 ¼ 10.78 ms, and the extracted pressure fields of target sources at (e) t1 ¼ 9.38 ms and (f) t2 ¼ 10.78 ms on the measurement plane H. The centers of the three speaker cones are marked with the symbol “þ.”
3480
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
waveform comparisons among the pressures measured after physically removing the disturbing source, the mixed pressures generated by the target sources and the disturbing source together, and the extracted pressures of target sources, as depicted in Fig. 8. The comparisons between the pressures measured after physically removing the disturbing source and the mixed pressures illustrate that the disturbing source S2 imposes a serious disturbance on the pressure signal generated by the target sources S1 and S3. To remove the disturbance due to the disturbing source and extract the desired signal related to the target sources from the mixed one, the RT-FSM was performed. As shown in Fig. 8, the extracted pressures of target sources at these three space points are in good agreement with their measured values. According to Eqs. (19) and (20), the values of Ep at the points R1 , R2 , and R3 are 0.9843, 0.9856, and 0.9275, respectively, and the values of Ea are 0.1116, 0.0275, and 0.0475, respectively. Figure 9 shows the values of Ep and Ea on the measurement plane H, which illustrates that the pressures measured after physically removing the disturbing source and the extracted pressures of target sources achieve the high phase similarity at all space points where the values of Ep are above 0.9 and the little amplitude difference at most space points where the values of Ea are below 0.1. Because the values at all the measurement points at any time instant were contaminated by the disturbing speaker S2, the spatial distribution of the pressure field at any time instant would be changed. Compared to Figs. 10(a) and 10(b) which show the pressure fields measured after physically removing the disturbing source at the times 9.38 ms and 10.78 ms, respectively, Figs. 10(c) and 10(d) with the mixed pressure fields at the same two time instants clearly present the difference of the spatial distribution. In Figs. 10(a) and 10(b), the locations of the two target speakers S1 and S3 are clearly revealed on the measurement plane, while the location of the disturbing speaker S2 dominates in Figs. 10(c) and 10(d). The extracted pressure fields of target sources at the same two time instants by RT-FSM are shown in Figs. 10(e) and 10(f). It is seen that the influence of the disturbing speaker S2 is removed effectively and the extracted pressure fields of target sources offer almost the same spatial distributions as those in Figs. 10(a) and 10(b), which indicates that the RT-FSM is effective in extracting the desired sound field generated by the target sources at any time instant from the measured pressure and particle acceleration fields in presence of the disturbing source. V. CONCLUSIONS
A method was developed to separate the time-evolving pressure field generated by the target source from the mixed one in real time by measuring the pressure and particle acceleration on one measurement plane, where the particle acceleration was obtained by the finite difference approximation with the aid of an auxiliary measurement plane. The separation formulation was deduced in the time-wave number domain on the basis of the superposition principle of waves and the impulse response function relating the J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
pressure to the particle acceleration. The proposed method has the ability of real-time separation of non-stationary sound fields, which provides the potential to in situ analyze the non-stationary radiation characteristics of the target source in presence of disturbing sources. A numerical simulation with two monopole sources as the target sources and one as the disturbing source has been carried out to verify the validity of RT-FSM. The results demonstrated that the proposed method can remove the influence of the disturbing source and extract the desired pressure signals in both time and space domains well. An experiment with three speakers driven by a Morlet wavelet signal also evidenced the feasibility of the proposed method. Furthermore, the extracted pressure signals could be used as the input of nearfield acoustic holography25–27 for further studying the vibration and sound radiation of target sources. ACKNOWLEDGMENT
This work was supported by the National Natural Science Foundation of China (Grant Nos. 51322505 and 11274087). APPENDIX: THE DERIVATION OF THE IMPULSE RESPONSE FUNCTION
Referring to the Euler’s equation, the relationship between the time-independent pressure and particle acceleration can be expressed as aðx; y; z; tÞ ¼
1 @pðx; y; z; tÞ ; q0 @z
(A1)
where aðx; y; z; tÞ is the particle acceleration in the z direction, pðx; y; z; tÞ is the pressure, and q0 is the medium density. By applying a two-dimensional Fourier transforms with respect to x and y and a Laplace transform with respect to t to Eq. (A1), it yields Aðkx ; ky ; z; sÞ ¼
1 @Pðkx ; ky ; z; sÞ ; q0 @z
(A2)
where kx and ky represent the wave-number components in x and y directions, respectively, and s represents the time domain Laplace variable. Considering a measurement plane located at z ¼ 0, the propagation of the sound field to the plane z > 0 can be expressed as28 pffiffiffiffiffiffiffiffiffi 2 Pðkx ; ky ; z; sÞ ¼ Pðkx ; ky ; 0; sÞez as þb ; (A3) where a ¼ 1=c2 and b ¼ kx2 þ ky2 . After substituting Eq. (A3) into Eq. (A2) and taking the derivation, it can be solved as q0 ffi Aðkx ; ky ; z; sÞ: Pðkx ; ky ; z; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi as2 þ b
(A4)
By setting q0 ffi; hðkx ; ky ; 0; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi as2 þ b Bi et al.: Real-time sound field separation
(A5)
3481
and taking the inverse Laplace transform to Eq. (A4), finally yields Pðkx ; ky ; zh ; tÞ ¼ Aðkx ; ky ; zs ; tÞ hðkx ; ky ; 0; tÞ;
(A6)
where the asterisk denotes a convolution of two time functions, and hðkx ; ky ; 0; tÞ is the impulse response function relating the pressure to the particle acceleration, whose analytical expression is achieved by solving the inverse Laplace transform of hðkx ; ky ; 0; sÞ in Eq. (A5) as29 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðkx ; ky ; 0; tÞ ¼ q0 cJ0 (A7) c2 ðkx2 þ ky2 Þt HðtÞ: In Eq. (A7), J0 is the Bessel function of the first kind of order zero, and HðtÞ is the Heaviside function. 1
G. V. Frisk, A. V. Oppenheim, and D. R. Martinez, “A technique for measuring the plane-wave reflection coefficient of the ocean bottom,” J. Acoust. Soc. Am. 68, 602–612 (1980). 2 D. L. Hallman and J. S. Bolton, “A technique for performing source identification in reflective environments by using nearfield acoustical holography,” in Proceedings of Noise-Con 93, Williamsburg, VA (May 2–5, 1993), pp. 479–484. 3 M. Tamura, “Spatial Fourier transform method of measuring reflection coefficients at oblique incidence. I: Theory and numerical examples,” J. Acoust. Soc. Am. 88, 2259–2264 (1990). 4 M. Tamura, J. F. Allard, and D. Lafarge, “Spatial Fourier-transform method for measuring reflection coefficients at oblique incidence. II. Experimental results,” J. Acoust. Soc. Am. 97, 2255–2262 (1995). 5 M. T. Cheng, J. A. Mann III, and A. Pate, “Wave-number domain separation of the incident and scattered sound field in Cartesian and cylindrical coordinates,” J. Acoust. Soc. Am. 97, 2293–2303 (1995). 6 E. G. Williams, “The nearfield acoustical holography (NAH) experimental method applied to vibration and radiation in light and heavy fluids,” Comput. Struct. 65, 323–335 (1997). 7 G. Weinreich and E. B. Arnold, “Method for measuring acoustic radiation fields,” J. Acoust. Soc. Am. 68, 404–411 (1980). 8 X. Zhao and S. F. Wu, “Reconstruction of vibroacoustic fields in halfspace by using hybrid near-field acoustical holography,” J. Acoust. Soc. Am. 117, 555–565 (2005). 9 Bavu, and A. Garcia, “Evaluation Y. Braikia, M. Melon, C. Langrenne, E. of a separation method for source identification in small spaces,” J. Acoust. Soc. Am. 134, 323–331 (2013). 10 J. Hald, “Patch holography in cabin environments using a two-layer handheld array with an extended SONAH algorithm,” in Proceedings of Euronoise 2006, Tampere, Finland (May 30–June 1, 2006). 11 F. Jacobsen and V. Jaud, “Statistically optimized near field acoustic holography using an array of pressure-velocity probes,” J. Acoust. Soc. Am. 121, 1550–1558 (2007).
3482
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
12
E. Fernandez-Grande and F. Jacobsen, “Sound field separation with a double layer velocity transducer array (L),” J. Acoust. Soc. Am. 130, 5–8 (2011). 13 Y. B. Zhang, X. Z. Chen, and F. Jacobsen, “A sound filed separation technique based on measurements with pressure-velocity probes (L),” J. Acoust. Soc. Am. 125, 3518–3521 (2009). 14 C. Langrenne, M. Melon, and A. Garcia, “Boundary element method for the acoustic characterization of a machine in bounded noisy environment,” J. Acoust. Soc. Am. 121, 2750–2757 (2007). 15 N. P. Valdivia, E. G. Williams, and P. C. Herdic, “Approximations of inverse boundary element methods with partial measurements of the pressure field,” J. Acoust. Soc. Am. 123, 109–120 (2008). 16 C. Langrenne, M. Melon, and A. Garcia, “Measurement of confined acoustic sources using near-field acoustic holography,” J. Acoust. Soc. Am. 126, 1250–1256 (2009). 17 C. X. Bi, X. Z. Chen, and J. Chen, “Sound field separation technique based on equivalent source method and its application in nearfield acoustic holography,” J. Acoust. Soc. Am. 123, 1472–1478 (2008). 18 C. X. Bi and J. S. Bolton, “An equivalent source technique for recovering the free sound field in a noisy environment,” J. Acoust. Soc. Am. 131, 1260–1270 (2012). 19 E. Fernandez-Grande, F. Jacobsen, and Q. Lecle`re, “Sound field separation with sound pressure and particle velocity measurements,” J. Acoust. Soc. Am. 132, 3818–3825 (2012). 20 M. Forbes, S. Letcher, and P. Stepanishen, “A wave vector, time-domain method of forward projecting time-dependent pressure fields,” J. Acoust. Soc. Am. 90, 2782–2793 (1991). 21 J. Sun and P. R. Stepanishen, “Modal impulse response and wave vector-time domain methods to evaluate transient radiation from axisymmetric planar radiators,” J. Acoust. Soc. Am. 102, 979–988 (1997). 22 X. Z. Zhang, J. H. Thomas, C. X. Bi, and J. C. Pascal, “Separation of nonstationary sound fields in the time-wave number domain,” J. Acoust. Soc. Am. 131, 2180–2189 (2012). 23 E. G. Williams, “Regularization methods for near-field acoustical holography,” J. Acoust. Soc. Am. 110, 1976–1988 (2001). 24 S. Paillasseur, J. H. Thomas, and J. C. Pascal, “Regularization for improving the deconvolution in real-time near-field acoustic holography,” J. Acoust. Soc. Am. 129, 3777–3787 (2011). 25 J. D. Maynard, E. G. Williams, and Y. Lee, “Nearfield acoustic holography: I. Theory of generalized holography and the development of NAH,” J. Acoust. Soc. Am. 78, 1395–1413 (1985). 26 W. A. Veronesi and J. D. Maynard, “Nearfield acoustic holography (NAH) II. Holographic reconstruction algorithms and computer implementation,” J. Acoust. Soc. Am. 81, 1307–1322 (1987). 27 E. G. Williams, Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography (Academic Press, San Diego, CA, 1999), pp. 1–306. 28 V. Grulier, S. Paillasseur, J. H. Thomas, J. C. Pascal, and J. C. Le Roux, “Forward propagation of time evolving acoustic pressure: Formulation and investigation of the impulse response in time-wave number domain,” J. Acoust. Soc. Am. 126, 2367–2378 (2009). 29 J. L. Schiff, The Laplace Transform: Theory and Applications (Springer, New York, 1999), pp. 1–233.
Bi et al.: Real-time sound field separation
Copyright of Journal of the Acoustical Society of America is the property of American Institute of Physics and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
(Received 24 September 2013; revised 28 February 2014; accepted 28 April 2014) To extract the desired non-stationary sound field generated by a target source in the presence of disturbing sources, a real-time sound field separation method with pressure and particle acceleration measurements is proposed. In this method, the pressure and particle acceleration signals at a time instant are first measured on one measurement plane, where the particle acceleration is obtained by the finite difference approximation with the aid of an auxiliary measurement plane; then, the desired pressure signal generated by the target source at the same time instant can be extracted in a timely manner, by a simple superposition of the measured pressure and the convolution between the measured particle acceleration and the derived impulse response function. Thereby, the proposed method possesses a significant feature of real-time separation of non-stationary sound fields, which provides the potential to in situ analyze the radiation characteristics of a non-stationary source. The proposed method was examined through numerical simulation and experiment. Results demonstrated that the proposed method can not only extract the desired time-evolving pressure signal generated by the target source at any space point, but can also obtain the desired spatial distribution of the pressure field generated by the target source at any time instant. C 2014 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4875576] V PACS number(s): 43.60.Gk, 43.20.Px, 43.60.Pt [MS]
I. INTRODUCTION
The sound field separation method (FSM) has been proven to be an effective technique for extracting the spatial acoustic characteristics of a target source in the presence of disturbing sources. Currently, several different methods, such as the spatial Fourier transform method,1–6 the spherical wave expansion method,7–9 the statistically optimized nearfield acoustic holography,10–13 the boundary element method,14–16 and the equivalent source method17–19 have been developed to realize this separation for distinguishing the sources on both sides of the measurement surface and suppressing the influence of disturbing sources. These FSMs require measurements of pressures or particle velocities on two measurement surfaces, or measurements of pressure and particle velocity on one measurement surface, and they all have the ability to extract the sound field generated by the target source in the frequency domain. However, for a nonstationary source whose statistical properties fluctuate in time, it is more reasonable to separate the sound field at a time instant for further analyzing the acoustic characteristics related to time. Recently, on the basis of the impulse response approach,20,21 a method for separating the non-stationary sound fields by measuring the pressures on two closely spaced measurement planes was developed by Zhang et al.22 In that method, the pressure signal generated by the target source at each time step is solved via a deconvolution process. The deconvolution process implies a complicated a)
Author to whom correspondence should be addressed. Electronic mail: [email protected]
3474
J. Acoust. Soc. Am. 135 (6), June 2014
Pages: 3474–3482
iterative calculation, i.e., the pressure signal generated by the target source at the present time step has to be calculated in the need of the pressure signals solved at all the previous time steps. Furthermore, the deconvolution process is an illposed inverse problem. To obtain an appropriate solution, the singular value decomposition (SVD) and the Tikhonov regularization23 have to be used, which leads to a large computational workload. In the present paper, an alternative method, namely, real-time sound field separation method (RT-FSM), is proposed to separate the non-stationary sound fields as well. In RT-FSM, the desired pressure signal generated by the target source at a specific time instant can be extracted by a simple superposition of the measured pressure at the same time instant and the convolution between the measured particle acceleration and the corresponding impulse response function. Here, the particle acceleration is obtained by the finite difference approximation with the pressures measured on two closely spaced measurement planes. Compared to the deconvolution process involved in Zhang’s method, the method proposed in the present paper only contains a simple forward solving process, i.e., the pressure signal generated by the target source at the present time step can be extracted with no need of any pressure signal solved at the previous time steps. And more importantly, the forward solving process in RT-FSM avoids the use of the SVD and the regularization, which reduces the computational workload significantly. Therefore, the RT-FSM provides a better feature of real-time separation than Zhang’s method, which gives the potential to in situ analyze the radiation characteristics of a non-stationary source. The paper is organized as follows. Section II presents a theoretical description of RT-FSM. In Sec. III, to examine
0001-4966/2014/135(6)/3474/9/$30.00
C 2014 Acoustical Society of America V
the performance of RT-FSM in both time and space domains, a numerical simulation with three monopole sources driven by different non-stationary signals is investigated. An experiment involving three speakers driven by a Morlet wavelet signal is carried out in Sec. IV. Finally, conclusions are summarized in Sec. V. II. REAL-TIME SOUND FIELD SEPARATION METHOD
Figure 1 illustrates the geometry of RT-FSM in the Cartesian coordinate system oðx; y; zÞ. Because of the joint contributions from the target source Mo and the disturbing source Md , the pressure and particle acceleration on the measurement plane H are both mixed. The purpose of the present paper is to remove the influence of the disturbing source Md and extract the time-evolving pressure signal generated by the target source Mo from the mixed one. According to the superposition principle of waves, the mixed pressure signal pðx; y; zH ; tÞ and the mixed particle acceleration signal aðx; y; zH ; tÞ at the measurement point ðx; y; zH Þ at the time t can be expressed, respectively, as pðx; y; zH ; tÞ ¼ po ðx; y; zH ; tÞ þ pd ðx; y; zH ; tÞ;
(1)
aðx; y; zH ; tÞ ¼ ao ðx; y; zH ; tÞ ad ðx; y; zH ; tÞ;
(2)
where the subscripts “o” and “d” denote the sound fields generated by the target source and the disturbing source, respectively. By applying a two-dimensional Fourier transform with respect to x and y to Eqs. (1) and (2), it yields
Pd ðkx ; ky ; zH ; tÞ ¼ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ;
(6)
where the asterisk denotes the convolution of two time functions. The expression of the impulse response function hðkx ; ky ; 0; tÞ is derived in the Appendix and given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðkx ; ky ; 0; tÞ ¼ q0 cJ0 (7) c2 ðkx2 þ ky2 Þt HðtÞ; where q0 is the medium density, c is the sound velocity, J0 is the Bessel function of the first kind of order zero, and HðtÞ is the Heaviside function. By replacing the term Pd ðkx ; ky ; zH ; tÞ in Eq. (3) with the right part of Eq. (6), Eq. (3) can be rewritten as Pðkx ; ky ; zH ; tÞ ¼ Po ðkx ; ky ; zH ; tÞ þ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ:
(8)
By convolving with the impulse response function hðkx ; ky ; 0; tÞ on both sides of Eq. (4), it yields Aðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ ¼ Ao ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ :
(9)
In accordance with Eq. (5), Eq. (9) can be further expressed as Aðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ ¼ Po ðkx ; ky ; zH ; tÞ Ad ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ:
Pðkx ; ky ; zH ; tÞ ¼ Po ðkx ; ky ; zH ; tÞ þ Pd ðkx ; ky ; zH ; tÞ;
(3)
(10)
Aðkx ; ky ; zH ; tÞ ¼ Ao ðkx ; ky ; zH ; tÞ Ad ðkx ; ky ; zH ; tÞ;
(4)
By combining Eq. (8) with Eq. (10), the time–wave number spectrum of the pressure Po ðkx ; ky ; zH ; tÞ can be finally obtained by
where kx and ky represent the wave number components in x and y directions, respectively. Referring to the derivation in Appendix, the time–wave number spectrum of the pressure can be related to that of the particle acceleration produced by the same source with the following formulas: Po ðkx ; ky ; zH ; tÞ ¼ Ao ðkx ; ky ; zH ; tÞ hðkx ; ky ; 0; tÞ;
(5)
Po ðkx ;ky ;zH ;tÞ ¼ 0:5½Pðkx ;ky ;zH ;tÞþ Aðkx ;ky ;zH ;tÞ hðkx ;ky ;0;tÞ:
(11)
Equation (11) indicates that the pressure spectrum Po ðkx ; ky ; zH ; tÞ generated by the target source on the measurement plane H can be extracted from the mixed pressure spectrum Pðkx ; ky ; zH ; tÞ with the aid of the mixed particle acceleration spectrum Aðkx ; ky ; zH ; tÞ. In Eq. (11), the convolution is continuous in the time domain. To meet the requirement of numerical computation, the time t should be discretized as tn ¼ ðn 1ÞDt;
(12)
where n ¼ 1; 2; …; N, N is the total number of time steps, and Dt is the time step. Accordingly, Eq. (11) can be expressed in a discrete form as follows: Po ðkx ;ky ;zH ;tn Þ ¼0:5 Pðkx ;ky ;zH ;tn Þ n X Aðkx ;ky ;zH ;ti Þhðkx ;ky ;0;tniþ1 Þ : þ FIG. 1. (Color online) Geometry of RT-FSM in the Cartesian coordinate system oðx; y; zÞ. J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
i¼1
(13) Bi et al.: Real-time sound field separation
3475
From Eq. (13), it is noted that once the mixed pressure spectrum at the time instant tn and the mixed particle acceleration spectrums from the time instant t1 to the time instant tn are known, the pressure spectrum Po ðkx ; ky ; zH ; tn Þ at the time instant tn can be obtained. In fact, the mixed particle acceleration spectrums from the time instant t1 to the time instant tn1 have been achieved before performing the measurement at the time instant tn , and therefore, the pressure spectrum Po ðkx ; ky ; zH ; tn Þ can be extracted in real time only by measuring the mixed pressure and particle acceleration at the time instant tn . Since the pressure field in the space domain is desired, the two-dimensional inverse Fourier transform with respect to kx and ky is then applied to the pressure spectrum Po ðkx ; ky ; zH ; tn Þ to obtain the pressure po ðx; y; zH ; tn Þ. The procedure of RT-FSM can be described clearly with a flow chart as shown in Fig. 2. In RT-FSM, the particle acceleration should be measured as one of the inputs to start the separation. Unfortunately, it is difficult to measure the particle acceleration directly. In the following, an alternative way to acquire the particle acceleration is provided, which is based on the Euler’s equation given by aðx; y; zH ; tÞ ¼
1 @pðx; y; zH ; tÞ : q0 @z
(14)
Equation (14) indicates that the particle acceleration is related to the partial derivative of the pressure. Since the partial derivative is often approximated by the finite difference, an auxiliary measurement plane H1 is placed near the measurement plane H and the partial derivative in Eq. (14) is approximated by the finite difference of the pressures on the planes H and H1 as @pðx; y; zH ; tÞ pðx; y; zH1 ; tÞ pðx; y; zH ; tÞ ; @z zH1 zH
(15)
where pðx; y; zH1 ; tÞ is the mixed pressure measured on the auxiliary measurement plane H1.
Substituting Eq. (15) into Eq. (14), it finally yields the approximate expression of the particle acceleration as aðx; y; zH ; tÞ
pðx; y; zH1 ; tÞ pðx; y; zH ; tÞ : q0 ðzH1 zH Þ
(16)
III. NUMERICAL SIMULATION
A numerical simulation with three monopole sources was investigated to verify the feasibility of the proposed method. In the simulation, the monopole sources S1 and S3 were fixed on one side of the measurement plane as the target sources and the monopole source S2 was put on the other side as the disturbing source. The geometry of sources and measurement planes is shown in Fig. 3. The source S1 located at (0.1 m, 0.45 m, 0.145 m) generated a non-stationary signal with a linear frequency modulation in the [200, 1800] Hz band and a Gaussian amplitude modulation. The source S2 located at (0.4 m, 0.45 m, 0.14 m) generated a Morlet wavelet signal defined by sðtÞ ¼ cosð2pf0 tÞet
2
=2
;
(17)
where f0 ¼ 344 Hz. The source S3 located at (0.55 m, 0.35 m, 0.14 m) generated a non-stationary signal with a linear frequency modulation in the [150, 2400] Hz band and a Gaussian amplitude modulation. The measurement plane H was located at z ¼ 0 m in the Cartesian coordinate system oðx; y; zÞ, and the auxiliary measurement plane H1 located 0.02 m away from the measurement plane H was used to acquire the particle acceleration. These two planes both provided 1515 measurement points and gave the grid spacing of 0.05 m in both x and y directions. The theoretical pressure at each measurement point was calculated by the following equation: pðR; tÞ ¼
sðt R=cÞ ; 4pR
(18)
FIG. 2. (Color online) Flow chart of RT-FSM. Fxy and F1 xy denote twodimensional Fourier transform and the two-dimensional inverse Fourier transform with respect to x and y, respectively.
3476
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
FIG. 3. (Color online) Geometry of the measurement plane H, the auxiliary measurement plane H1, and the locations of three monopole sources in the Cartesian coordinate system oðx; y; zÞ. Three space points R1 , R2 , and R3 are chosen and marked with the symbol “þ.” The points R1 and R3 are facing the target sources S1 and S3, respectively, and the point R2 is facing the disturbing source S2.
where sðtÞ denotes the input signal, and R denotes the distance between the measurement point and the monopole source. The signals were sampled at a frequency of fe ¼25.6 kHz, providing 512 sampling points. A Gaussian white noise with a signal-to-noise ratio of 20 dB was added to the measured signals to simulate the practical measurement. To lessen the wrap-around errors due to the use of discrete twodimensional spatial Fourier transforms associated with the finite size of measurement planes, the 1515 pressure and particle acceleration matrices on the measurement plane H were extended to 2929 matrices by zero-padding. To assess the relevance of the proposed method in the time domain, three space points on the measurement plane H, identified in Fig. 3, were chosen and their positions were R1 (0.1 m, 0.45 m, 0 m), R2 (0.4 m, 0.45 m, 0 m), and R3 (0.55 m, 0.35 m, 0 m). The points R1 and R3 were facing the target sources S1 and S3, respectively, and the point R2 was facing the disturbing source S2. The time-domain waveform
comparisons among the theoretical pressures calculated after physically removing the disturbing source (line with plus sign), the mixed pressures generated by the target sources and the disturbing source together (dotted line), and the extracted pressures of target sources (solid line) at these three space points are shown in Fig. 4. The comparison results between the lines with a plus sign and the dotted lines indicate that the theoretical pressure signals radiated by the target sources are affected by the disturbing source in different degrees at these three points. To remove the influence of the disturbing source and extract the pressure generated by the target sources, the RT-FSM was performed. The comparison results between the lines with a plus sign and the solid lines indicate that the extracted pressures using the proposed method match well with the theoretical pressures radiated by the target sources. To estimate the consistency between the theoretical pressure calculated after physically removing the disturbing
FIG. 4. (Color online) Time-domain waveform comparisons among the theoretical pressures calculated after physically removing the disturbing source (line with plus sign), the mixed pressures generated by the target sources and the disturbing source together (dotted line), and the extracted pressures of target sources (solid line) at the three space points (a) R1 , (b) R2 , and (c) R3 .
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
3477
FIG. 5. (Color online) Spatial maps of (a) Ep with a contour line at the value of 0.9 and (b) Ea with a contour line at the value of 0.1. The space points R1 , R2 , and R3 are marked with the symbol “þ.”
source and the extracted pressure of target sources quantitatively, two indicators related to the phase and the amplitude were defined, respectively, by24 hpm ðx; y; z; tÞpe ðx; y; z; tÞi ; Epðx; y; zÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hp2m ðx; y; z; tÞihp2e ðx; y; z; tÞi
(19)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hp2 ðx; y; z; tÞi hp2 ðx; y; z; tÞi m e pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Eaðx; y; zÞ ¼ ; hp2m ðx; y; z; tÞi
(20)
where pm ðx; y; z; tÞ is the theoretical pressure calculated after physically removing the disturbing source, pe ðx; y; z; tÞ is the
FIG. 6. (Color online) Spatial distributions of the theoretical pressure fields calculated after physically removing the disturbing source at (a) t1 ¼ 9.18 ms and (b) t2 ¼ 12.86 ms, the mixed pressure fields generated by the target sources and the disturbing source together at (c) t1 ¼ 9.18 ms and (d) t2 ¼ 12.86 ms, and the extracted pressure fields of target sources at (e) t1 ¼ 9.18 ms and (f) t2 ¼ 12.86 ms on the measurement plane H. The three monopole sources are marked with the symbol “þ.”
3478
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
FIG. 7. (Color online) Experimental setup: three speakers, a microphone array, the recording system and the trigger. The speakers S1 and S3 are used as the target sources, and the speaker S2 is used as the disturbing source.
extracted pressure of target sources, and h i denotes the time averaged value. Here, it is noted that when Ep is close to unity, there exists a high phase similarity, and when Ea is near zero, the amplitude difference is small. To quantify the accuracy of the extracted results at all points on the measurement plane H, the spatial maps of the indicators Ep and Ea calculated by Eqs. (19) and (20), respectively, are shown in Fig. 5. The values of Ep and Ea at a majority of space points are more than 0.9 and below 0.1, respectively, which indicates that the extracted phases and amplitudes agree with the theoretical ones very well. Here, the values of Ep at the points R1 , R2 , and R3 are 0.9701, 0.9254, and 0.9713, respectively, and the values of Ea are 0.0513, 0.1399, and 0.0376, respectively.
Similarly, for showing the extracted results in the space domain, two time instants, t1 ¼9.18 ms and t2 ¼12.86 ms, were chosen. Figures 6(a), 6(c), and 6(e) depict the theoretical pressure field calculated after physically removing the disturbing source, the mixed pressure field generated by the target sources and the disturbing source together, and the extracted pressure field of target sources at t1 ¼9.18 ms, respectively. Figures 6(b), 6(d), and 6(f) depict the same pressure fields but at t2 ¼12.86 ms. The comparisons between the theoretical pressure fields in Figs. 6(a) and 6(b) and the mixed pressure fields in Figs. 6(c) and 6(d) at these two time instants illustrate that the spatial distributions of theoretical pressure fields have been changed by the disturbing source. After performing the RT-FSM, it is obvious that the extracted pressure fields in Figs. 6(e) and 6(f) are almost the same as the theoretical ones in Figs. 6(a) and 6(b), which indicates that the proposed method seems effective to suppress the sound radiated by the disturbing source in the space domain. IV. EXPERIMENTAL STUDY
An experiment was carried out in a semi-anechoic chamber to further verify the proposed method. In the experiment, three speakers S1, S2, and S3 with the same sizes were used as the sources, among which the speakers S1 and S3 were the target sources and the speaker S2 was the disturbing source, as shown in Fig. 7. All speakers were connected with a laptop and all driven by the same Morlet wavelet signal as that used in the simulation case. In the experiment, an array with 15 microphones was moved 15 times on each measurement surface to measure the pressure signals at 1515 measurement points. To keep
FIG. 8. (Color online) Time-domain waveform comparisons among the pressures measured after physically removing the disturbing source (line with plus sign), the mixed pressures generated by the target sources and the disturbing source together (dotted line), and the extracted pressures of target sources (solid line) at the three space points (a) R1 , (b) R2 , and (c) R3 .
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
3479
FIG. 9. (Color online) Spatial maps of (a) Ep with a contour line at the value of 0.9 and (b) Ea with a contour line at the value of 0.1. The space points R1 , R2 , and R3 are marked with the symbol “þ.”
the measured signals the same at each measurement, the output of the source signal was set as the trigger to activate the acquisition system recording the data. In the experiment, the centers of the three speaker cones were positioned at the same locations of the three sources in the simulation case, and the locations of measurement planes and the grid spacing of measurement points were all the same
as those in the simulation case. The signals were sampled at a frequency of fe ¼25.6 kHz, providing 512 sampling points. The 1515 pressure and particle acceleration matrices on the measurement plane H were also extended to 2929 matrices by zero-padding to lessen the wrap-around errors. Just as the simulation case, the same three space points R1 , R2 , and R3 were selected to show the time-domain
FIG. 10. (Color online) Spatial distributions of the pressure fields measured after physically removing the disturbing source at (a) t1 ¼ 9.38 ms and (b) t2 ¼ 10.78 ms, the mixed pressure fields generated by the target sources and the disturbing source together at (c) t1 ¼ 9.38 ms and (d) t2 ¼ 10.78 ms, and the extracted pressure fields of target sources at (e) t1 ¼ 9.38 ms and (f) t2 ¼ 10.78 ms on the measurement plane H. The centers of the three speaker cones are marked with the symbol “þ.”
3480
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
Bi et al.: Real-time sound field separation
waveform comparisons among the pressures measured after physically removing the disturbing source, the mixed pressures generated by the target sources and the disturbing source together, and the extracted pressures of target sources, as depicted in Fig. 8. The comparisons between the pressures measured after physically removing the disturbing source and the mixed pressures illustrate that the disturbing source S2 imposes a serious disturbance on the pressure signal generated by the target sources S1 and S3. To remove the disturbance due to the disturbing source and extract the desired signal related to the target sources from the mixed one, the RT-FSM was performed. As shown in Fig. 8, the extracted pressures of target sources at these three space points are in good agreement with their measured values. According to Eqs. (19) and (20), the values of Ep at the points R1 , R2 , and R3 are 0.9843, 0.9856, and 0.9275, respectively, and the values of Ea are 0.1116, 0.0275, and 0.0475, respectively. Figure 9 shows the values of Ep and Ea on the measurement plane H, which illustrates that the pressures measured after physically removing the disturbing source and the extracted pressures of target sources achieve the high phase similarity at all space points where the values of Ep are above 0.9 and the little amplitude difference at most space points where the values of Ea are below 0.1. Because the values at all the measurement points at any time instant were contaminated by the disturbing speaker S2, the spatial distribution of the pressure field at any time instant would be changed. Compared to Figs. 10(a) and 10(b) which show the pressure fields measured after physically removing the disturbing source at the times 9.38 ms and 10.78 ms, respectively, Figs. 10(c) and 10(d) with the mixed pressure fields at the same two time instants clearly present the difference of the spatial distribution. In Figs. 10(a) and 10(b), the locations of the two target speakers S1 and S3 are clearly revealed on the measurement plane, while the location of the disturbing speaker S2 dominates in Figs. 10(c) and 10(d). The extracted pressure fields of target sources at the same two time instants by RT-FSM are shown in Figs. 10(e) and 10(f). It is seen that the influence of the disturbing speaker S2 is removed effectively and the extracted pressure fields of target sources offer almost the same spatial distributions as those in Figs. 10(a) and 10(b), which indicates that the RT-FSM is effective in extracting the desired sound field generated by the target sources at any time instant from the measured pressure and particle acceleration fields in presence of the disturbing source. V. CONCLUSIONS
A method was developed to separate the time-evolving pressure field generated by the target source from the mixed one in real time by measuring the pressure and particle acceleration on one measurement plane, where the particle acceleration was obtained by the finite difference approximation with the aid of an auxiliary measurement plane. The separation formulation was deduced in the time-wave number domain on the basis of the superposition principle of waves and the impulse response function relating the J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
pressure to the particle acceleration. The proposed method has the ability of real-time separation of non-stationary sound fields, which provides the potential to in situ analyze the non-stationary radiation characteristics of the target source in presence of disturbing sources. A numerical simulation with two monopole sources as the target sources and one as the disturbing source has been carried out to verify the validity of RT-FSM. The results demonstrated that the proposed method can remove the influence of the disturbing source and extract the desired pressure signals in both time and space domains well. An experiment with three speakers driven by a Morlet wavelet signal also evidenced the feasibility of the proposed method. Furthermore, the extracted pressure signals could be used as the input of nearfield acoustic holography25–27 for further studying the vibration and sound radiation of target sources. ACKNOWLEDGMENT
This work was supported by the National Natural Science Foundation of China (Grant Nos. 51322505 and 11274087). APPENDIX: THE DERIVATION OF THE IMPULSE RESPONSE FUNCTION
Referring to the Euler’s equation, the relationship between the time-independent pressure and particle acceleration can be expressed as aðx; y; z; tÞ ¼
1 @pðx; y; z; tÞ ; q0 @z
(A1)
where aðx; y; z; tÞ is the particle acceleration in the z direction, pðx; y; z; tÞ is the pressure, and q0 is the medium density. By applying a two-dimensional Fourier transforms with respect to x and y and a Laplace transform with respect to t to Eq. (A1), it yields Aðkx ; ky ; z; sÞ ¼
1 @Pðkx ; ky ; z; sÞ ; q0 @z
(A2)
where kx and ky represent the wave-number components in x and y directions, respectively, and s represents the time domain Laplace variable. Considering a measurement plane located at z ¼ 0, the propagation of the sound field to the plane z > 0 can be expressed as28 pffiffiffiffiffiffiffiffiffi 2 Pðkx ; ky ; z; sÞ ¼ Pðkx ; ky ; 0; sÞez as þb ; (A3) where a ¼ 1=c2 and b ¼ kx2 þ ky2 . After substituting Eq. (A3) into Eq. (A2) and taking the derivation, it can be solved as q0 ffi Aðkx ; ky ; z; sÞ: Pðkx ; ky ; z; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi as2 þ b
(A4)
By setting q0 ffi; hðkx ; ky ; 0; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi as2 þ b Bi et al.: Real-time sound field separation
(A5)
3481
and taking the inverse Laplace transform to Eq. (A4), finally yields Pðkx ; ky ; zh ; tÞ ¼ Aðkx ; ky ; zs ; tÞ hðkx ; ky ; 0; tÞ;
(A6)
where the asterisk denotes a convolution of two time functions, and hðkx ; ky ; 0; tÞ is the impulse response function relating the pressure to the particle acceleration, whose analytical expression is achieved by solving the inverse Laplace transform of hðkx ; ky ; 0; sÞ in Eq. (A5) as29 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðkx ; ky ; 0; tÞ ¼ q0 cJ0 (A7) c2 ðkx2 þ ky2 Þt HðtÞ: In Eq. (A7), J0 is the Bessel function of the first kind of order zero, and HðtÞ is the Heaviside function. 1
G. V. Frisk, A. V. Oppenheim, and D. R. Martinez, “A technique for measuring the plane-wave reflection coefficient of the ocean bottom,” J. Acoust. Soc. Am. 68, 602–612 (1980). 2 D. L. Hallman and J. S. Bolton, “A technique for performing source identification in reflective environments by using nearfield acoustical holography,” in Proceedings of Noise-Con 93, Williamsburg, VA (May 2–5, 1993), pp. 479–484. 3 M. Tamura, “Spatial Fourier transform method of measuring reflection coefficients at oblique incidence. I: Theory and numerical examples,” J. Acoust. Soc. Am. 88, 2259–2264 (1990). 4 M. Tamura, J. F. Allard, and D. Lafarge, “Spatial Fourier-transform method for measuring reflection coefficients at oblique incidence. II. Experimental results,” J. Acoust. Soc. Am. 97, 2255–2262 (1995). 5 M. T. Cheng, J. A. Mann III, and A. Pate, “Wave-number domain separation of the incident and scattered sound field in Cartesian and cylindrical coordinates,” J. Acoust. Soc. Am. 97, 2293–2303 (1995). 6 E. G. Williams, “The nearfield acoustical holography (NAH) experimental method applied to vibration and radiation in light and heavy fluids,” Comput. Struct. 65, 323–335 (1997). 7 G. Weinreich and E. B. Arnold, “Method for measuring acoustic radiation fields,” J. Acoust. Soc. Am. 68, 404–411 (1980). 8 X. Zhao and S. F. Wu, “Reconstruction of vibroacoustic fields in halfspace by using hybrid near-field acoustical holography,” J. Acoust. Soc. Am. 117, 555–565 (2005). 9 Bavu, and A. Garcia, “Evaluation Y. Braikia, M. Melon, C. Langrenne, E. of a separation method for source identification in small spaces,” J. Acoust. Soc. Am. 134, 323–331 (2013). 10 J. Hald, “Patch holography in cabin environments using a two-layer handheld array with an extended SONAH algorithm,” in Proceedings of Euronoise 2006, Tampere, Finland (May 30–June 1, 2006). 11 F. Jacobsen and V. Jaud, “Statistically optimized near field acoustic holography using an array of pressure-velocity probes,” J. Acoust. Soc. Am. 121, 1550–1558 (2007).
3482
J. Acoust. Soc. Am., Vol. 135, No. 6, June 2014
12
E. Fernandez-Grande and F. Jacobsen, “Sound field separation with a double layer velocity transducer array (L),” J. Acoust. Soc. Am. 130, 5–8 (2011). 13 Y. B. Zhang, X. Z. Chen, and F. Jacobsen, “A sound filed separation technique based on measurements with pressure-velocity probes (L),” J. Acoust. Soc. Am. 125, 3518–3521 (2009). 14 C. Langrenne, M. Melon, and A. Garcia, “Boundary element method for the acoustic characterization of a machine in bounded noisy environment,” J. Acoust. Soc. Am. 121, 2750–2757 (2007). 15 N. P. Valdivia, E. G. Williams, and P. C. Herdic, “Approximations of inverse boundary element methods with partial measurements of the pressure field,” J. Acoust. Soc. Am. 123, 109–120 (2008). 16 C. Langrenne, M. Melon, and A. Garcia, “Measurement of confined acoustic sources using near-field acoustic holography,” J. Acoust. Soc. Am. 126, 1250–1256 (2009). 17 C. X. Bi, X. Z. Chen, and J. Chen, “Sound field separation technique based on equivalent source method and its application in nearfield acoustic holography,” J. Acoust. Soc. Am. 123, 1472–1478 (2008). 18 C. X. Bi and J. S. Bolton, “An equivalent source technique for recovering the free sound field in a noisy environment,” J. Acoust. Soc. Am. 131, 1260–1270 (2012). 19 E. Fernandez-Grande, F. Jacobsen, and Q. Lecle`re, “Sound field separation with sound pressure and particle velocity measurements,” J. Acoust. Soc. Am. 132, 3818–3825 (2012). 20 M. Forbes, S. Letcher, and P. Stepanishen, “A wave vector, time-domain method of forward projecting time-dependent pressure fields,” J. Acoust. Soc. Am. 90, 2782–2793 (1991). 21 J. Sun and P. R. Stepanishen, “Modal impulse response and wave vector-time domain methods to evaluate transient radiation from axisymmetric planar radiators,” J. Acoust. Soc. Am. 102, 979–988 (1997). 22 X. Z. Zhang, J. H. Thomas, C. X. Bi, and J. C. Pascal, “Separation of nonstationary sound fields in the time-wave number domain,” J. Acoust. Soc. Am. 131, 2180–2189 (2012). 23 E. G. Williams, “Regularization methods for near-field acoustical holography,” J. Acoust. Soc. Am. 110, 1976–1988 (2001). 24 S. Paillasseur, J. H. Thomas, and J. C. Pascal, “Regularization for improving the deconvolution in real-time near-field acoustic holography,” J. Acoust. Soc. Am. 129, 3777–3787 (2011). 25 J. D. Maynard, E. G. Williams, and Y. Lee, “Nearfield acoustic holography: I. Theory of generalized holography and the development of NAH,” J. Acoust. Soc. Am. 78, 1395–1413 (1985). 26 W. A. Veronesi and J. D. Maynard, “Nearfield acoustic holography (NAH) II. Holographic reconstruction algorithms and computer implementation,” J. Acoust. Soc. Am. 81, 1307–1322 (1987). 27 E. G. Williams, Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography (Academic Press, San Diego, CA, 1999), pp. 1–306. 28 V. Grulier, S. Paillasseur, J. H. Thomas, J. C. Pascal, and J. C. Le Roux, “Forward propagation of time evolving acoustic pressure: Formulation and investigation of the impulse response in time-wave number domain,” J. Acoust. Soc. Am. 126, 2367–2378 (2009). 29 J. L. Schiff, The Laplace Transform: Theory and Applications (Springer, New York, 1999), pp. 1–233.
Bi et al.: Real-time sound field separation
Copyright of Journal of the Acoustical Society of America is the property of American Institute of Physics and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
