{'fraction_non_alphanumeric': 0.04013719325466498, 'fraction_numerical': 0.013474337511739006, 'mean_word_length': 4.368697939500219, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-time separation of non-stationary sound fields on spheres.

The sound field separation methods can separate the target field from the interfering noises, facilitating the study of the acoustic characteristics of the target source, which is placed in a noisy environment. However, most of the existing sound field separation methods are derived in the frequency-domain, thus are best suited for separating stationary sound fields. In this paper, a time-domain sound field separation method is developed that can separate the non-stationary sound field generated by the target source over a sphere in real-time. A spherical array sets up a boundary between the target source and the interfering sources, such that the outgoing field on the array is only generated by the target source. The proposed method decomposes the pressure and the radial particle velocity measured by the array into spherical harmonics coefficients, and recoveries the target outgoing field based on the time-domain relationship between the decomposition coefficients and the theoretically derived spatial filter responses. Simulations show the proposed method can separate non-stationary sound fields both in free field and room environments, and over a longer duration with small errors. The proposed method could serve as a foundation for developing future time-domain spatial sound field manipulation algorithms.

The sound field separation methods can separate the target field from the interfering noises, facilitating the study of the acoustic characteristics of the target source, which is placed in a noisy environment. However, most of the existing sound field separation methods are derived in the frequency-domain, thus are best suited for separating stationary sound fields. In this paper, a time-domain sound field separation method is developed that can separate the non-stationary sound field generated by the target source over a sphere in real-time. A spherical array sets up a boundary between the target source and the interfering sources, such that the outgoing field on the array is only generated by the target source. The proposed method decomposes the pressure and the radial particle velocity measured by the array into spherical harmonics coefficients, and recoveries the target outgoing field based on the time-domain relationship between the decomposition coefficients and the theoretically derived spatial filter responses. Simulations show the proposed method can separate non-stationary sound fields both in free field and room environments, and over a longer duration with small errors. The proposed method could serve as a foundation for developing future time-domain spatial sound field manipulation algorithms. a) Author to whom correspondence should be addressed. Electronic mail: fei.ma@ieee.org

I. INTRODUCTION
To investigate the working condition of a machine or the acoustic characteristics of a musical instrument, we can use a sensor array to measure the target sound field generated by them and analyze the sensor array measurement. However, in practical acoustic environments, such as inside rooms, there are disturbing sources, which produce disturbing sound fields. The disturbing sound fields and the room reflected sound fields will interfere with the target sound field and contaminate the sensor array measurements, making it difficult to study the target source characteristics. Nonetheless, the target source does not necessarily co-locates with the disturbing sources. By exploiting the spatial characteristics of the sound sources, it is possible to separate the target sound field from the interfering sound fields using a sensor array. Thus we can study the target source through analyzing the separated target sound field.
Over the years, a number of sound field separation (SFS) methods have been reported in the literature. The spherical wave expansion based method uses the spherical harmonic modes as the basic modeling functions of a sound field, and are capable of separating the outgoing and incoming fields on a spherical sensor array. [1][2][3][4] The Spatial Fourier Transform based method uses the two-dimensional Fourier Transform to decompose the sound field into plane-wave components, and can separates the incident and reflected sound fields on a planar sensor array. 5 The statistically optimized near field holography method is also capable of performing SFS on a planar sensor array, and can mitigate the spatial window leakage problem of the Spatial Fourier Transform based method. 6 The recently developed boundary element based method 7,8 and the equivalent source based method 9-11 extend the application of SFS to arbitrary shape sensor arrays. Further, these two methods can separate the scattering from the target source surface and recover the free-field radiation of the target source.
The above mentioned SFS methods can separate the sound fields incoming from two sides of a sensor array (or two sensor arrays) apart. However, they are all derived in the frequencydomain. That is, they first accumulate and transfer a frame of the time-domain pressure (or particle velocity) measurements into the time-frequency domain using the Short Time Fourier Transform. They then conduct the SFS at each time-frequency bin. The Short Time Fourier Transform process inevitably introduces the spectrum leakage problem and frame-delays into the SFS methods. 12 The frame-delay is acceptable if the target sound field is stationary (or quasi-stationary). Nonetheless, when dealing with fast changing nonstationary sound field, the frame-delay will make the SFS methods unable to recover the target field in real-time. That further makes the separated target field unable to be used in time-critical applications, such as active noise control, 13 real-time beamforming, and machine anomaly diagnosis. 7 A time-domain SFS method, on the other head, can track the changes of the target sound field without introducing the frame latency issue.
Manipulations of the sound fields over spatial regions in the time-domain are difficult, and time-domain SFS methods are seldom developed. The only existing time-domain SFS methods are developed by Bi and his coauthors. 14-17 Based on the Spatial Fourier Transform, the method they developed can separate the sound field coming from two sides of a planar sensor array. 14,15,17 Based on the interpolated time-domain equivalent source method, 16 they further developed a method that can separate the sound generated by a particular source in a multiple non-stationary sources scene. However, these time-domain methods are best suited for SFS in free field (or semi-anechoic rooms). Because first, in room environments, due to wall reflections the interfering sounds and the target sound can arrive at a planar sensor array in the same directions, 18,19 such that the target sound and the interfering noise become indistinguishable. Second, the equivalent source based method requires perfect knowledge of the location and geometry information of all the sources. However, in a room environment, the image source effect exists, and it is difficult to characterize the image source especially for an arbitrary shape room with irregular boundaries. 18,19 In this paper, based on the spherical wave expansion, we develop a time-domain SFS method that can separate the non-stationary outgoing and incoming sound fields on a sphere in both free field and room environments. We first decompose the pressure and the radial particle velocity into corresponding spherical harmonics coefficients. The timedomain convolution between the spherical harmonic coefficients and the derived impulse response functions results in the outgoing/incoming field on the array. The performance of the proposed method is confirmed by simulations, and compared with the Spatial Fourier Transform based method. Possible applications of the proposed method include (1) reference signal generation for active noise control system, 13 (2) real-time beamforming, (3) de-reverberation for speech recognition and sound field reproduction system, (4) machine working condition monitoring in a noisy environment.
This paper is organized as follows. We introduce the problem formulation and the frequency-domain SFS method based on spherical harmonic expansion in Sec. II. In Sec. III, we provide the theoretical derivation of the time-domain SFS method. Sec. IV introduces the practical implementations of the proposed method, whose effectiveness is validated by simulations in Sec. V, and Sec. VI concludes this paper.

A. Problem formulation
Consider the systems shown in Fig. 1, where in Fig. 1 (a) the target sources (denoted as △) are placed in free field, and in Fig. 1 (b) the target sources (denoted as △) are placed inside a room. We use (x, y, z) and (r, θ, φ) to denote the Cartesian and spherical coordinates of a point with respect to the point O, respectively. The sound field p(t, R, θ, φ) on the sphere S 2 of radius R is the superposition of the outgoing field p o (t, R, θ, φ) and the incoming field p i (t, R, θ, φ) where t is the continuous time. The frequency-domain representation of (1) is where ω = 2πf is the angular frequency (f is the frequency).
As shown in Fig. 1, the incoming field p i (t, R, θ, φ) on the sphere S 2 is due to the external sources or room reflections. The outgoing field p o (t, R, θ, φ) on the sphere S 2 , on the other hand, is generated by the target sources, thus characterizes the target sources. In this paper, we aim to separate the outgoing field p o (t, R, θ, φ) on the sphere S 2 to facilitate the study of the target sources.
B. The frequency-domain sound field separation on a sphere In the frequency-domain, the sound field P (ω, R, θ, φ) on a sphere S 2 of radius R consists of the outgoing field P o (ω, R, θ, φ) and the incoming field P i (ω, R, θ, φ) 1 where c is the speed of sound, A µν (ω, R) are spherical harmonic coefficients correspond the sound field P (ω, R, θ, φ) measured on the sphere S 2 , A o µν (ω, R) and A i µν (ω, R) are the outgoing and incoming field coefficients, respectively. h µ (·) is the second kind spherical Hankel function of order u, j µ (·) is the first kind spherical Bessel function of order u, L µν (·) and K µν (·) are the (radial-independent) spherical Hankel and Bessel function coefficients, respectively, Y µν (·, ·) is the real spherical harmonic of order n and degree m, 20 i.e., P µ|ν| (·) is the associated Legendre function of order µ and degree |ν|. Hereafter, we abbreviate (θ, φ) as a single symbol Θ to simplify the notations.
The frequency-domain radial particle velocity on the sphere S 2 can be similarly expressed where B µν (ω, R) are spherical harmonic coefficients corresponding to the radial particle velocity V (ω, R, Θ) measured on the sphere S 2 , h ′ µ (·) and j ′ µ (·) are derivatives of h µ (·) and j µ (·) about the argument, respectively, i is the unit imaginary number, and ρ 0 is the density of air.
Based on (3) and (5), we obtain the outgoing and incoming field coefficients as 1 respectively. Substitution of the coefficients A o µν (ω, R) and A i µν (ω, R) back into (3)  duce the frame-delay, is more appropriate to separate the non-stationary outgoing field for applications which have strict low-latency requirements.

III. THE TIME-DOMAIN SOUND FIELD SEPARATION ON A SPHERE
In this section, we derive the time-domain SFS method on a sphere. We start with a theorem, which provides the expressions of the outgoing/incoming field on a sphere. The detailed derivation of the theorem is in the appendix.
where the outgoing field coefficients a o µν (t, R) and the incoming field coefficients a i µν (t, R) are given by respectively, τ R = R/c, p(t, R, Θ) and v(t, R, Θ) are the pressure and the radial particle velocity on the sphere at the point (R, Θ), respectively, and Θ dΘ ≡ 2π 0 π 0 sin θdθdφ. The We have the following remarks on the time-domain SFS method: 1. Equations (8), (9), (10), (11) show that the outgoing (incoming) field at a particular point (R, Θ) can be described by the pressure, the pressure gradient, and the radial partial velocity gradient observed over the sphere within a certain time frame.
2. The spatial filter functions g u 0 (t), g u 1 (t), g u 2 (t), g u 3 (t), and g u 4 (t) are specially designed such that they take finite values inside of the time interval [0, 2τ R ] and are zeros outside of the interval. This limits the integrals (10) and (11) to the interval [0, 2τ R ].

METHOD
The theoretically derived (8), (9), (10), and (11) of the proposed method involve continuous integrals, time derivatives of the sound pressure and the radial particle velocity, which are acoustic quantities that can not be measured directly using existing sensors. In this section, we further develop approaches to facilitate the implementation of the proposed method.

A. Discretization and sampling
We first approximate the continuous integral in (10) by a discrete form as follows where δ τ = δ t = 1/f s , and T n = ⌈2f s τ R ⌉ + 1 is the number of samples corresponding to 2τ R (f s is the sampling frequency, and ⌈·⌉ is the ceiling operator).
We further develop (17) by replacing the continuous integral over the sphere with a finite summation at (R, Θ q ) Q q=1 and by substituting the continuous time t and τ n using the discrete time n and n ′ . By denoting p(n, R, Θ q ) and v(n, R, Θ q ) as the pressure and the radial partial velocity at (R, Θ q ) at discrete time instant n, we have the outgoing field coefficient a o µν (n, R) at discrete time n as where {γ q } Q q=1 are the sampling weights. 28 Substitution of a o µν (t, R) by a o µν (n, R) into (8) results in an estimation of the outgoing field p o (n, R, Θ) at discrete time instance on the sphere S 2 . We can obtain the discrete form of (11) in a similar manner.
Due to the finite sampling over the sphere, we can only recover the outgoing field coefficients up to N, where N ≤ √ Q is the truncation order of the outgoing field. 23

B. Pressure and velocity approximation
We may use an array of acoustic vector sensors, which can measure both the pressure and the radial particle velocity, 25 to realize the SFS method. Alternatively, we can use two microphone arrays to realize the method.
Place Q microphones on a sphere of radius R − δ R at {R − δ R , Θ q } Q q=1 , and another Q microphones on a sphere of radius R + δ R at {R + δ R , Θ q } Q q=1 , where δ R ≪ R is a real positive number. We approximate the pressure on the middle sphere of radius R at positions Based on the Euler's equation, we approximate the radial particle velocity on the middle Substitutions of (19) and (20)

V. SIMULATION EXAMPLES
In this section, we conduct simulations to validate the effectiveness of the proposed timedomain SFS method.

A. Sound field separation in free field
The simulation environment is as shown in Fig. 1 (a). We place a point source at (0, 0, 0.3) m as the target source, and 100 plane waves, whose directions are determined according to the 100-point spherical packing, 27 as the incoming fields. We let the target source output a band-limited sound, which is simulated by passing a unit-variance Gauss white noise through a 64-tap Bandpass Butterworth filter of the frequency range [100, 600] Hz. The incoming plane waves have the same frequency components as the target source output. The strengths of the plane waves are normalized such that their contributions to the sound pressure on the sphere are approximately the same compared with the contribution from the target source.
We aim to estimate the target (outgoing) sound field on the sphere S 2 produced by the target source. We choose the radius of the sphere S 2 as R = 0.65 m, and the outgoing field truncation order to be N = 5. We use a sensor array arranged on the sphere according to the 6-th order Gauss sampling scheme to realize the time-domain SFS method. 28 We simulate the impulse responses, including the radial particle velocity response, between the source and the sensors based on the free-space Green function and the Euler's equation. 1 The impulse responses are truncated to 1024 taps long under the sampling frequency f s = 48 kHz. The speed of sound is c = 343 m/s, and the air density is ρ 0 = 1.225 kg/m 3 . The length of impulse response functions g µ 0 (t), g µ 1 (t), g µ 2 (t) is T n = ⌈2R/cf s ⌉ + 1 = 183 taps.
We add Gauss white noise to the sensor measurement such that the signal to noise ratio is 40 dB. The simulation results are from an average of 100 independent runs. The settings in this paragraph are used in all simulations unless otherwise stated.
We reconstruct the outgoing field on the sphere over a period of 10 ms according to (8).
We depict the amplitudes of the outgoing field p o (t, R, Θ 17 ), the total field p(t, R, Θ 17 ), the  As shown in Fig. 3, the total field differs from the outgoing field. The separated outgoing field approximates the outgoing field, and the field separation error is small over the whole We depict the amplitudes of the outgoing field p o (t, R, Θ q ), the total field p(t, R, Θ q ), the separated outgoing fieldp o (t, R, Θ q ), and the outgoing field separation error p e (t, R, , on the sphere S 2 at t = 10 ms in Fig. 4. Here (R, Θ q ) 180×360 q=1 are equal-angle sampling point positions on the sphere S 2 . 28 As shown in Fig. 4, the proposed method is able to recover the outgoing field accurately over the whole sphere, and the normalized separation error over the whole sphere is  where p o (ω, R, Θ q ) is the frequency-domain outgoing field,p o,N (ω) is the frequency-domain counterpart of the following equation We plot the normalized separation error ξ N (ω) in Fig. 5, which shows that ξ N (ω) increases along with the frequency. The higher the outgoing field truncation order N is the smaller the normalized separation error ξ N (ω) is. However, in the low frequency range, the truncation order N = 5 does not make the normalized separation error ξ N (ω) significantly smaller than that in the case of N = 4, i.e., ξ 5 (ω) ≈ ξ 4 (ω) for ω < 800 π/rad or f < 400 Hz. To conclude, we should choose an appropriate truncation order N for the separated target (outgoing) field depending on its frequency components such that the separation error is sufficiently small.

B. Comparisons with the Spatial Fourier Transform based method
In this section, we compare the performance of the proposed method with the Spatial Fourier Transform based method. 15 The setup of the simulation environment is shown in Fig. 6, which has a Cartesian coordinate system with the origin at the point O. On the z = 0 m plane H, we equally distribute 15×15 sensors on a square of size 0.7 m × 0.7 m, whose center is at the point O.
On the z = 0.01 m plane H 1 , we place another 15×15 sensors equally on a square of size 0.7 m × 0.7 m, whose center is the at the point (x, y, z) = (0.0, 0.0, 0.01) m. We use the sensors on these two planes for SFS using the Spatial Fourier Transform based method. We extend the array size to be 81×81 by zero-padding to reduce the aliasing error of the Spatial Fourier Transform. 15 There is a sphere S 2 of radius R = 0.35 m. On the sphere S 2 , we place sensors according to the 9-th order Gauss sampling scheme. 28 We use the sensors on this sphere S 2 to estimate the outgoing field up to 6-th order using the proposed method.
The sampling schemes used by these two methods are selected such that they can detect approximately the same number of acoustic quantities, and that the sizes of the two arrays are comparable.   As shown in Fig. 7, the proposed method can recover the target sound accurately over the whole 30 ms period. The Spatial Fourier Transform based method can separate the target field with small errors before t = 5 ms, but after that the SFS error starts to increase. That is because, in the Spatial Fourier Transform based method, the target sound at one time instant depends on all previous sounds. 15 The truncation of the infinite long Bessel function, i.e., (7) from Bi, 15 to a finite length introduces errors. We have conducted simulations using the Spatial Fourier Transform based method with a longer truncation length (8192 taps) for the Bessel functions. The sound field separation error are similar to the dotted line as shown in Fig. 7, thus are not shown for brevity. In the proposed method we have specially designed the spatial filter functions functions, g µ 0 (t), g µ 1 (t), g µ 2 (t), to be finite long, avoiding SFS error due to the impulse response function truncation.
The simulation results demonstrate that the Spatial Fourier Transform-based method, which does not require the sensor array to surround the target source, is best suited for separating transient plane wave sound fields. The proposed method, on the other hand, can accurately separate the target (outgoing) sound field on a spherical array over a longer duration.

C. Sound field separation in a room
In practical implementations of the proposed SFS method, the target sources may be placed on a surface. In case the surface is rigid or highly reflective, we can measure the sound field using a semi-spherical array around the target sources, and duplicate the measurements with respect to the surface to use the full spherical harmonics expansion. 3,4 In this section, we simulate such a case.

Sound field separation
We measure the sound pressure and the radial particle velocity on the semi-sphere S 2 at four sensors, whose locations are determined according to the first-order Gauss sampling scheme. 28 We simulate the impulse responses, including the radial particle velocity response, between each source and sensor using the image source method. 19 We take the first 2744 image sources into consideration, and the impulse responses are truncated to 8192 taps long under the sampling frequency of 48 kHz. The pressure and the radial partial velocity on the semi-sphere S 2 are duplicated with respect to the floor to use the full spherical harmonic expansions. 3,4 We estimate the outgoing field coefficients up to 0-th order, and reconstruct the outgoing field on the sphere over a period of 10 ms according to (8). We depict the amplitudes of the outgoing field p o (t, R, Θ 1 ), the total field p(t, R, Θ 1 ), the separated outgoing field p o (t, R, Θ 1 ), and the outgoing field separation error p e (t, R, at the first sensor in Fig. 9. The outgoing field is produced by convolving the target source outputs and corresponding the half-space Green function. 3,4 The total filed is the summation of the outgoing and incoming fields.
As shown in Fig. 9, in this case, the proposed method is still able to accurately estimate the outgoing field. The normalized separation error at the first sensor is The proposed method also recovers the outgoing field over the whole semi-sphere. The results are similar to Fig. 4, thus are not shown for brevity.  FIG. 9. Sound field separation at a point: The amplitudes of the outgoing field p o (t, R, Θ 1 ), the total field p(t, R, Θ 1 ), the separated outgoing fieldp o (t, R, Θ 1 ), and the outgoing field separation error p e (t, R, Θ 1 ), at the first sensor over 10 ms.

VI. CONCLUSION
In this paper, we developed a time-domain SFS method that can separate non-stationary sound fields over a sphere. We decompose the sound field and the radial particle velocity measure on the sphere into spherical harmonic coefficients, and recover the outgoing/incoming field based on the time-domain relationship between the coefficients and the derived sepa-ration filters. The simulations demonstrated that the proposed method can separate nonstationary sound fields in both free field and room environments. Further, the method is able to accurately recover the outgoing field over a long period of time. A future extension of the proposed method is to take the scattering from the target source surface into consideration. 7,8,10 The experimental validation of the proposed method will be our future work. Substitutions of τ R = R/c and 22 We define the time-domain outgoing field coefficients a o µν (t, R) as where F −1 denotes the inverse Fourier transform. 12 The terms in the second line of (A2) are defined based on (A1) as follows η µν (t) ≡ F −1 (iω)B µν (ω, R) , We derive expressions for these six terms in the following.
The expressions of the incoming field coefficients a i µν (t, R), g µ 3 (t), and g µ 4 (t) can be derived similarly, and are not shown for brevity.