# A Multiple Sound Source Localization and Counting Method Based on the Kernel Density Estimator

Authors: Yuzhuo Fang; Xian Zang; Juan Yang; Hongcheng Zhou; Zhiyong Xu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

An unambiguous signal processing algorithm when using a wide intermicrophone distance is proposed in this paper for simultaneously locating and counting multiple active sound sources. Based on the kernel density estimator, a multistage structure in the time-frequency domain is used to suppress the influence of spatial aliasing, then the pooled angular spectrum is combined with a peak search method having an updated cut-off threshold and a source merging module. Complete source localization and counting is realized through the combination of these two steps. Simulation results show that the proposed method has a more robust performance than the classic counterpart, especially in adverse environments with spatial aliasing, reverberation, and interference between different sound sources.


---

## Body

## 1. Introduction

Multiple sound source localization is a key element in many applications ranging from national security to video conference [1, 2]. Especially in far-field speech interaction, localizing multiple active sound sources accurately is an essential prerequisite for high-quality speech enhancement and recognition [3, 4]. In these applications, there is often ambient noise, reverberation, and mutual interference between different sources, which makes the source localization performance seriously affected.In recent years, based on W-disjoint orthogonality (WDO) [5, 6] of observed signals between pairwise microphones, multiple sound source localization methods can be further divided into clustering, histogram, and angular spectrum.The clustering method, such as interphase difference (IPD) [7] and direction estimation of mixing matrix (DEMIX) [8], directly achieves localization results by iteratively clustering time-frequency (TF) bins associated with each sound source [9]. It is sensitive to the initial clustering parameter [10]. The histogram method computes the weighted histogram such as circular integrated cross spectrum (CICS) [11] and smooth histogram [12] to make it have high estimation accuracy. However, the used array topology is relatively difficult to implement and popularize. Most of the above two methods do not consider spatial aliasing [13] which may exist in multiple sound source localization.From the space sampling theorem, the intermicrophone distance should be smaller than half the minimum signal wavelength (e.g. 2 cm for a sampling rate of 16 kHz). When the distance exceeds half the wavelength, it can be regarded as a wide distance, and spatial aliasing generated from the high-frequency part will occur. The wider the distance, the more serious the aliasing. The excessively wide distance will have a significant impact on sound source localization. To obtain a high resolution in the time difference of arrival (TDOA) space, the arrangement of the microphones often uses a distance far exceeding half the wavelength, especially in the case of a small number of microphones, thus inevitably generating the influence of aliasing [14]. Hence, we focus on the angular spectrum method because of its applicability for a wide intermicrophone distance [10, 15]. It consists of two steps: (i) spectrum construction and (ii) localization and counting.In the first step, the angular spectrum is constructed by accumulating the local function related to all possible angles in each TF bin. Classical generalized crosscorrelation (GCC) [16] computes the crosspower spectrum phase [17] from the pairwise microphones. Because of the limitations of the ideal single-source propagation model [18], the GCC spectrum will be severely distorted in environments with the coexistence of reverberation and spatial aliasing. Generated from the generalized state coherence transform (GSCT) [19, 20], approximate kernel density estimator (KDE) using the nonlinear Gaussian kernel [20, 21] has a significantly improved angular spectrum than GCC under reverberant environment. The frequency-dependent weighting factor of the Gaussian kernel can suppress spatial aliasing to a certain extent, but the effect is limited [14].Existing sound source localization methods often regard the number of sound sources as a prior information [7, 22]. However, during the utterance of mixed speech signal, some sources may only last for a short period of time and the number of sources often changes dynamically [23], which makes it difficult to determine the number of sources in advance.Therefore, the second step requires the active sources to be located and counted simultaneously [8] from the constructed spectrum. Current source localization and counting methods mainly use iterative search methods, such as single-point peak amplitude in peak search (PS) [24, 25], inner product in matching pursuit (MP) [26, 27], and source contribution rate in iterative contribution removal (ICR) [15]. Each iteration selects the optimal value satisfying the corresponding conditions, removes the components corresponding to the current sound source from the spectrum, and restarts the next. MP and ICR are more accurate than PS when obvious distortion is not produced in the spectrum. However, in a high-resolution spectrum, the computation cost of the inner product and the contribution rate is much higher, which is not conductive to real-time tracking the number of sources at different times. When the spectrum is seriously distorted due to reverberation and spatial aliasing, the current incorrect estimation will deteriorate the reconstructed spectrum, which may exert considerable influence on the following iteration [28].In this paper, to effectively solve the spatial aliasing when using a wide intermicrophone distance, an unambiguous angular spectrum-based algorithm is proposed for simultaneously localizing and counting multiple active sound sources. In the spectrum construction step, the local KDE spectrum in each TF bin between the pairwise microphones is generated from a Gaussian kernel function; then a multi-stage structure [14] is used to divide the entire band into several sub-bands by the maximum unambiguous frequency. The sub-band spectra are pooled together across all the TF bins to suppress the influence of spatial aliasing. In the localization and counting step, PS is used considering the simplicity and real-time processing. An updated cut-off threshold according to the current peak amplitude is used to replace the fixed threshold in traditional PS, which strengthen the flexibility of counting. Then the preliminary estimation is passed through a source merging module [15] to eliminate the duplicate sound source.The remainder of this paper is organized as follows: in Section2, the KDE spectrum is constructed on the basis of the signal propagation model of multiple sound sources. In Section 3, the MS structure is used to construct the KDEMS spectrum. Then PS are combined to form the complete localization and counting process. Numerical comparisons between the proposed KDEMS-PS and the other classical ones are given in Section 4. Finally, Section 5 concludes the paper.


## 2. Signal Propagation Model of Multiple Sound Sources and KDE Spectrum Construction

Setsni|n=1,2,…,N as N sound source signals at the transmitting end, and xmi|m=1,2,…,M as the observed signals corresponding to M microphones at the receiving end. The sampling rate is fs. Then in the approximately far field, the discrete time signal propagation model of multiple sound sources can be expressed as follows:(1)xmi=∑n=1Nhm,nTsni+vmi,where hm,n=hm,n0,…,hm,nL−1T denotes the vector of length L corresponding to the impulse response between the n-th source and the m-th microphone, sni=sni,…,sni−L+1T denotes the vector corresponding to the n-th source, and vmi denotes the additive white Gaussian noise of the m-th microphone which is independent of source signals and impulse responses.ThroughNFFT point short-time Fourier transform (STFT), the expression in the discrete TF domain can be obtained as(2)Xr,k=∑n=1NHnkSnr,k+Vr,k,where r and k denote the frame and frequency indices, respectively, Xr,k=X1r,k,…,Xmr,k,…,XMr,kT, Snr,k and Xmr,k denote the STFT coefficients corresponding to sni and xmi, respectively, Vr,k∈ℂM×1 denotes the complex vector corresponding to the noise, and(3)Hnk=H1,nk,…,Hm,nk,…,HM,nkT,where Hm,nk is the transfer function between the n-th source and the m-th microphone, including the direct wave component Hm,ndirk and the reverberation component Hm,nrevk. Since the impulse response is time-invariant, the transfer function is only related to k.When the intensity of the reflected wave is uniformly distributed in all possible directions of propagation,Hm,nrevk can be regarded as spatially diffuse noise [29, 30]. Hm,nk can be decomposed as(4)Hm,nk=Hm,ndirk+Hm,nrevk,where Hm,ndirk=αm,nexp−j2πfkim,n/fs, fk represents the frequency at the k-th frequency bin, αm,n and im,n/fs denote the corresponding propagation attenuation and time, respectively. Under anechoic, noise-free, and WDO conditions [5], only Hm,ndirk exists, and at most one sound source energy dominates each TF bin. Then equation (2) can be simplified as follows:(5)Xr,k=α1,ηr,kkexp−j2πfki1,ηr,kfs⋮αM,ηr,kkexp−j2πfkiM,ηr,kfsSηr,kr,k,where ηr,k∈1,…,N is the index of the dominant sound source in the current TF bin.Considering one pair of microphonesma, mb in the array, the normalized cross-power spectrum (NCS) of the observed signal can be expressed as [21].(6)NCSr,k=Xar,kXb∗r,kXar,kXb∗r,k=exp−j2πfkτηr,k,where τηr,k=ia,ηr,k−ib,ηr,k/fs denotes the true TDOA of the dominant sound source between ma and mb. Without loss of generality, the subscript of τηr,k can be omitted and τ can be estimated by regarding it as a random variable that satisfies a given probability distribution. Then a Gaussian kernel function is introduced as [20].(7)geτ=12πfkexp−eτ/2πfk22hK2,where(8)eτ=exp−j2πfkτ−NCSr,k,is a function of τ to calculate the Euclidean distance between exp−j2πfkτ and NCSr,k,(9)hK=τmaxB=dmaxcB,is the bandwidth of the kernel function where τmax is the maximum possible TDOA, dmax is the maximum spacing between adjacent microphones across the given pairs in the array, c is the velocity of sound propagation, and B is a factor that affects the resolution in the TDOA space. If B is set too large, the kernel function may generate many burrs, which may subsequently lead to the distortion of the spectrum. If B is set too small, the main lobe is obese, which is not conducive to locate the source correctly. B is set as 20 empirically when fs=16kHz. Substitute equation (8) into equation (7), the local KDE spectrum in each TF bin can be expressed as(10)φr,k,τ=12πfkexp−exp−j2πfkτ−NCSr,k22σk2,where σk=2πfkhK.Pool the local spectra across all the TF bins [10, 21]; the total KDE spectrum can be constructed as(11)ΦKDEτ=∑kφkτ=∑r,kφr,k,τ,where φkτ denotes the narrow-band KDE spectrum corresponding to the k-th frequency bin after pooling the local spectra across all time frames. The estimated TDOAs of sound sources can be obtained from the peaks of ΦKDEτ.


## 3. Localization and Counting Based on Multistage Structure and Peak Search

There is a frequency-dependent weighting factor1/2πfk in φr,k,τ, which means that the higher the frequency is, the more the amplitude of the local spectrum is suppressed. Then spatial aliasing mainly caused by the high-frequency part is weakened to a certain degree, but it cannot be completely eliminated [14], especially in a strong reverberation environment where the spectrum is severely distorted. Hence, a multistage structure is used to eliminate the influence of spatial aliasing more efficiently in this section.


### 3.1. Spectrum Construction with Multistage Structure

When calculatingΦKDEτ in equation (11), the indices of the lowest frequency bin and the highest are set to kL and kH, respectively. Then the corresponding center frequencies can be set to fL=kL⋅fs/NFFT and fH=kH⋅fs/NFFT, respectively.According to the spatial Nyquist sampling theorem, the unambiguous condition can be expressed as(12)dmic≤λmin2,where dmic is the distance between ma and mb, and λmin represents the minimum signal wavelength corresponding to fH which can be expressed as(13)λmin=cfH.Substituting equation (13) into equation (12), the condition can be rewritten as(14)fH≤c2dmic=fUA,where fUA denotes the maximum unambiguous frequency. If equation (14) is fulfilled, no spatial aliasing exists. However, in order to produce a higher resolution angular spectrum to distinguish different sound sources, it is necessary to use a wider dmic, which makes fH exceed fUA inevitably.IffH>fUA, the entire frequency band can be divided into two sub-bands fL,fUA and fUA,fH where a single sub-band may contain one or more continuous frequency bins. Then if fH>2fUA, fUA,fH can be further divided into fUA,2fUA and 2fUA,fH. By analogy, if fH>pfUA where p>1, the sub-bands divided at the p-th stage can be deduced as p−1fUA,pfUA and pfUA,fH. The entire band can be finally divided into a total of P sub-bands fL,fUA,…,P−1fUA,fH. P can be obtained as(15)P=ceilfH−fLfUA,where ceil⋅ denotes the round-up operator.Then by poolingφkτ across all frequency bins contained in each sub-band, the sub-band KDE spectrum can be obtained as(16)ΦKDEpτ=∑k=kpLkpHφkτ,where p=1,…,P, kpL denotes the index of the lowest frequency bin contained in the p-th sub-band which can be expressed as(17)kpL=kL,p=1,p−1kUA+1,p=2,…,P,and kpH denotes the index of the highest frequency bin which can be expressed as(18)kpH=pkUA,p=1,…,P−1,kH,p=P.In equations (17) and (18), kUA is the index corresponding to fUA which can be obtained as(19)kUA=floorfUA⋅NFFTfs,where floor⋅ is the round-down operator. When p=1, the frequency will not exceed fUA, so there is no interference of spatial aliasing; when p>1, equation (14) is no longer fulfilled, and the p-th sub-band spectrum ΦKDEpτ will have at most one more false peak than the previous sub-band.Then the MS structure is used for sub-band processing, and the output of the first stage can be expressed as(20)ΦKDEMS1τ=ΦKDE1τ.By multiplyingΦKDEMS1τ with the second sub-band spectrum, the output of the second stage can be obtained as(21)ΦKDEMS2τ=ΦKDEMS1τΦKDE2τ=ΦKDE1τΦKDE2τ,where the false peaks in ΦKDE2τ can be suppressed by ΦKDEMS1τ. Then the weighted output ΦKDEMS2τ is substituted into the next stage. By analogy, the output of the P-th stage, that is, the KDEMS spectrum can be deduced as(22)ΦKDEMSτ=ΦKDEMSPτ=ΦKDEMSP−1τΦKDEPτ=∏p=1PΦKDEpτ.The spectrum of the lower sub-band has a wider main lobe and lower resolution but is less affected by spatial aliasing, while the spectrum of the higher sub-band has a narrower main lobe and higher resolution but is subject to spatial aliasing. Combining these two parts of the spectrum through the MS structure can effectively absorb the advantages of each part while making up for their shortcomings, resulting in the final output spectrum with obvious ambiguity suppression and high resolution.Comparing equations (11) and (22), it can be seen that the MS structure does not increase the computational complexity. The total number of operations for both KDE and KDEMS can be approximately expressed as(23)Ncom=Nφ×Nτ×kH−kL+1,where Nφ denotes the number of operations required for φkτ, and Nτ denotes the number of samples contained in the TDOA space.


### 3.2. Multiple Sound Source Localization and Counting by Peak Search

When the angular spectrum is constructed, an improved PS is used to realize multiple sound source localization and counting. Without loss of generality and to facilitate the subsequent processing, the normalized form of the spectrum is obtained as(24)Φτ=ΦALτ−minΦALτmaxΦALτ−minΦALτ,where AL represents KDE, KDEMS, or other angular spectrum-based methods (e.g. GCC) and min⋅ and max⋅ represent the operators to find the minimum and maximum values in the function, respectively. Then the sequence containing all the peaks in Φτ can be expressed as(25)Φ=Φτ1,…,ΦτQ,where τ1,…,τQ denote the TDOAs corresponding to the peaks in the TDOA space sampled with τgrid as the grid value and Q denotes the total number of peaks. The elements in equation (25) are sorted in the descending order of peak amplitude, so the condition Φτi≥Φτi+1 where i=1,…,Q−1 is fulfilled. When the number of sound sources N is known where N≤Q, only the first N elements need to be extracted from Φ, and the set of the estimated TDOAs can be expressed as τi|i=1,…,N. However, in most real cases, N is unknown. Hence, it is necessary to count the number of sound sources to obtain effective localization results.In traditional PS, a fixed cut-off threshold set asΓ is compared with the elements in Φ one by one iteratively according to the order of the sequence until the element is not greater than Γ. If Γ is set too large or too small, it will lead to excessive missing alarm rate or false alarm rate, respectively. So the threshold is set to vary as the peak amplitude changes to improve the counting flexibility. Set Γ1 as the first threshold. It is generally assumed that there is at least one active sound source, so Φτ1>Γ1 holds. Then the information of the first peak amplitude is introduced into the threshold setting and the threshold used in the next iteration can be updated as Γ2=maxΦτ1/2,Γ1 where max⋅ represents the operators to find the maximum in the set. If Φτ2>Γ2, the iteration continues. Set the index of the iteration to i. Then the threshold used in the i-th iteration can be updated as(26)Γi=maxΦτi−1/2,Γ1,where i=2,…,I. I is the total number of iterations. Give the following condition:(27)Φτi≤Γii>Imax,where Imax denotes the maximum number of iterations. When either of the two conditions in equation (27) is not fulfilled, the iteration stops. The set of the estimated TDOAs can be expressed as(28)Ω=τi|i=1,…,I−1.Due to the interference of reverberation and noise, there may be multiple peaks with similar amplitude near the peak corresponding to the true TDOA in the distorted angular spectrum, which will lead to repeated estimation of the same sound source, thus bringing undesirable counting results. To address this problem, it is necessary to merge the duplicate estimation. The source merging is implemented as follows: for any two estimated TDOAsτi and τi′ where i≠i′, the condition can be given as(29)τi−τi′<τmin,where τmin indicates the allowed minimum distance between the estimated TDOAs corresponding to two different sound sources. τmin can be obtained as [15](30)τmin=dmicsinAminc,where Amin is the minimum angular distance empirically set to 10° when the interval is 0,180°. When equation (29) is fulfilled, according to the peak amplitude, the estimated results can be reassigned as(31)τi′=τi,Φτi≥Φτi′,τi=τi′,otherwise.After the reassignment, all the estimated TDOAs with the same value are merged into only one. Then the final localization and counting results can be expressed as(32)Ω′=τ^n|n=1,…,N^,where Ω′⊆Ω, τ^n indicates the estimated TDOA of the n-th sound source when using PS with the source merging module, and N^=cardΩ′ indicates the estimated number of sound sources where card⋅ is the operator to find the number of elements in the set.A complete localization and counting method can be used by combining the constructed angular spectrum with PS, where the combination is marked with “-” in this paper. Hence, when the angular spectrum is generated by GCC, KDE, or KDEMS, the corresponding localization and counting method can be called GCC-PS, KDE-PS or KDEMS-PS, respectively. The block diagram of the proposed KDEMS-PS is shown in Figure1, and the steps can be summarized as follows:(i)
STFT: Transform the observed signals of one pair of microphonesxai, xbi into the corresponding coefficients Xar,k, Xbr,k in each TF bin by STFT.(ii)
Narrowband KDE spectrum calculation: CalculateNCSr,k from equation (6), and then pool the local KDE spectrum φr,k,τ in equation (10) across the k-th frequency bin to obtain the narrowband KDE spectrum φkτ.(iii)
Multi-stage processing: Determine the number of sub-bandsP and the index corresponding to the maximum unambiguous frequency kUA from equations (15) and (19), respectively, then obtain the sub-band KDE spectrum ΦKDEpτ where p=1,…,P from equation (16).(iv)
Pooling and normalization: Pool all the sub-band spectra to construct the KDEMS angular spectrumΦKDEMSτ from equation (22), then obtain the normalized form Φτ from equation (24).(v)
Peak search: Search the peaks inΦτ, sort them in descending order, and then compare them with the updated cut-off threshold Γi from updated equation (26) to obtain the estimated TDOAs in Ω.(vi)
Source merging: Merge the repeatedly estimated sound sources from equation (31) to obtain the final localization and counting results in Ω′.Figure 1
A block diagram of the multiple sound source localization and counting method combined KDEMS with PS.


## 3.1. Spectrum Construction with Multistage Structure

When calculatingΦKDEτ in equation (11), the indices of the lowest frequency bin and the highest are set to kL and kH, respectively. Then the corresponding center frequencies can be set to fL=kL⋅fs/NFFT and fH=kH⋅fs/NFFT, respectively.According to the spatial Nyquist sampling theorem, the unambiguous condition can be expressed as(12)dmic≤λmin2,where dmic is the distance between ma and mb, and λmin represents the minimum signal wavelength corresponding to fH which can be expressed as(13)λmin=cfH.Substituting equation (13) into equation (12), the condition can be rewritten as(14)fH≤c2dmic=fUA,where fUA denotes the maximum unambiguous frequency. If equation (14) is fulfilled, no spatial aliasing exists. However, in order to produce a higher resolution angular spectrum to distinguish different sound sources, it is necessary to use a wider dmic, which makes fH exceed fUA inevitably.IffH>fUA, the entire frequency band can be divided into two sub-bands fL,fUA and fUA,fH where a single sub-band may contain one or more continuous frequency bins. Then if fH>2fUA, fUA,fH can be further divided into fUA,2fUA and 2fUA,fH. By analogy, if fH>pfUA where p>1, the sub-bands divided at the p-th stage can be deduced as p−1fUA,pfUA and pfUA,fH. The entire band can be finally divided into a total of P sub-bands fL,fUA,…,P−1fUA,fH. P can be obtained as(15)P=ceilfH−fLfUA,where ceil⋅ denotes the round-up operator.Then by poolingφkτ across all frequency bins contained in each sub-band, the sub-band KDE spectrum can be obtained as(16)ΦKDEpτ=∑k=kpLkpHφkτ,where p=1,…,P, kpL denotes the index of the lowest frequency bin contained in the p-th sub-band which can be expressed as(17)kpL=kL,p=1,p−1kUA+1,p=2,…,P,and kpH denotes the index of the highest frequency bin which can be expressed as(18)kpH=pkUA,p=1,…,P−1,kH,p=P.In equations (17) and (18), kUA is the index corresponding to fUA which can be obtained as(19)kUA=floorfUA⋅NFFTfs,where floor⋅ is the round-down operator. When p=1, the frequency will not exceed fUA, so there is no interference of spatial aliasing; when p>1, equation (14) is no longer fulfilled, and the p-th sub-band spectrum ΦKDEpτ will have at most one more false peak than the previous sub-band.Then the MS structure is used for sub-band processing, and the output of the first stage can be expressed as(20)ΦKDEMS1τ=ΦKDE1τ.By multiplyingΦKDEMS1τ with the second sub-band spectrum, the output of the second stage can be obtained as(21)ΦKDEMS2τ=ΦKDEMS1τΦKDE2τ=ΦKDE1τΦKDE2τ,where the false peaks in ΦKDE2τ can be suppressed by ΦKDEMS1τ. Then the weighted output ΦKDEMS2τ is substituted into the next stage. By analogy, the output of the P-th stage, that is, the KDEMS spectrum can be deduced as(22)ΦKDEMSτ=ΦKDEMSPτ=ΦKDEMSP−1τΦKDEPτ=∏p=1PΦKDEpτ.The spectrum of the lower sub-band has a wider main lobe and lower resolution but is less affected by spatial aliasing, while the spectrum of the higher sub-band has a narrower main lobe and higher resolution but is subject to spatial aliasing. Combining these two parts of the spectrum through the MS structure can effectively absorb the advantages of each part while making up for their shortcomings, resulting in the final output spectrum with obvious ambiguity suppression and high resolution.Comparing equations (11) and (22), it can be seen that the MS structure does not increase the computational complexity. The total number of operations for both KDE and KDEMS can be approximately expressed as(23)Ncom=Nφ×Nτ×kH−kL+1,where Nφ denotes the number of operations required for φkτ, and Nτ denotes the number of samples contained in the TDOA space.


## 3.2. Multiple Sound Source Localization and Counting by Peak Search

When the angular spectrum is constructed, an improved PS is used to realize multiple sound source localization and counting. Without loss of generality and to facilitate the subsequent processing, the normalized form of the spectrum is obtained as(24)Φτ=ΦALτ−minΦALτmaxΦALτ−minΦALτ,where AL represents KDE, KDEMS, or other angular spectrum-based methods (e.g. GCC) and min⋅ and max⋅ represent the operators to find the minimum and maximum values in the function, respectively. Then the sequence containing all the peaks in Φτ can be expressed as(25)Φ=Φτ1,…,ΦτQ,where τ1,…,τQ denote the TDOAs corresponding to the peaks in the TDOA space sampled with τgrid as the grid value and Q denotes the total number of peaks. The elements in equation (25) are sorted in the descending order of peak amplitude, so the condition Φτi≥Φτi+1 where i=1,…,Q−1 is fulfilled. When the number of sound sources N is known where N≤Q, only the first N elements need to be extracted from Φ, and the set of the estimated TDOAs can be expressed as τi|i=1,…,N. However, in most real cases, N is unknown. Hence, it is necessary to count the number of sound sources to obtain effective localization results.In traditional PS, a fixed cut-off threshold set asΓ is compared with the elements in Φ one by one iteratively according to the order of the sequence until the element is not greater than Γ. If Γ is set too large or too small, it will lead to excessive missing alarm rate or false alarm rate, respectively. So the threshold is set to vary as the peak amplitude changes to improve the counting flexibility. Set Γ1 as the first threshold. It is generally assumed that there is at least one active sound source, so Φτ1>Γ1 holds. Then the information of the first peak amplitude is introduced into the threshold setting and the threshold used in the next iteration can be updated as Γ2=maxΦτ1/2,Γ1 where max⋅ represents the operators to find the maximum in the set. If Φτ2>Γ2, the iteration continues. Set the index of the iteration to i. Then the threshold used in the i-th iteration can be updated as(26)Γi=maxΦτi−1/2,Γ1,where i=2,…,I. I is the total number of iterations. Give the following condition:(27)Φτi≤Γii>Imax,where Imax denotes the maximum number of iterations. When either of the two conditions in equation (27) is not fulfilled, the iteration stops. The set of the estimated TDOAs can be expressed as(28)Ω=τi|i=1,…,I−1.Due to the interference of reverberation and noise, there may be multiple peaks with similar amplitude near the peak corresponding to the true TDOA in the distorted angular spectrum, which will lead to repeated estimation of the same sound source, thus bringing undesirable counting results. To address this problem, it is necessary to merge the duplicate estimation. The source merging is implemented as follows: for any two estimated TDOAsτi and τi′ where i≠i′, the condition can be given as(29)τi−τi′<τmin,where τmin indicates the allowed minimum distance between the estimated TDOAs corresponding to two different sound sources. τmin can be obtained as [15](30)τmin=dmicsinAminc,where Amin is the minimum angular distance empirically set to 10° when the interval is 0,180°. When equation (29) is fulfilled, according to the peak amplitude, the estimated results can be reassigned as(31)τi′=τi,Φτi≥Φτi′,τi=τi′,otherwise.After the reassignment, all the estimated TDOAs with the same value are merged into only one. Then the final localization and counting results can be expressed as(32)Ω′=τ^n|n=1,…,N^,where Ω′⊆Ω, τ^n indicates the estimated TDOA of the n-th sound source when using PS with the source merging module, and N^=cardΩ′ indicates the estimated number of sound sources where card⋅ is the operator to find the number of elements in the set.A complete localization and counting method can be used by combining the constructed angular spectrum with PS, where the combination is marked with “-” in this paper. Hence, when the angular spectrum is generated by GCC, KDE, or KDEMS, the corresponding localization and counting method can be called GCC-PS, KDE-PS or KDEMS-PS, respectively. The block diagram of the proposed KDEMS-PS is shown in Figure1, and the steps can be summarized as follows:(i)
STFT: Transform the observed signals of one pair of microphonesxai, xbi into the corresponding coefficients Xar,k, Xbr,k in each TF bin by STFT.(ii)
Narrowband KDE spectrum calculation: CalculateNCSr,k from equation (6), and then pool the local KDE spectrum φr,k,τ in equation (10) across the k-th frequency bin to obtain the narrowband KDE spectrum φkτ.(iii)
Multi-stage processing: Determine the number of sub-bandsP and the index corresponding to the maximum unambiguous frequency kUA from equations (15) and (19), respectively, then obtain the sub-band KDE spectrum ΦKDEpτ where p=1,…,P from equation (16).(iv)
Pooling and normalization: Pool all the sub-band spectra to construct the KDEMS angular spectrumΦKDEMSτ from equation (22), then obtain the normalized form Φτ from equation (24).(v)
Peak search: Search the peaks inΦτ, sort them in descending order, and then compare them with the updated cut-off threshold Γi from updated equation (26) to obtain the estimated TDOAs in Ω.(vi)
Source merging: Merge the repeatedly estimated sound sources from equation (31) to obtain the final localization and counting results in Ω′.Figure 1
A block diagram of the multiple sound source localization and counting method combined KDEMS with PS.


## 4. Numerical Analysis

In this section, compared with classic GCC and KDE, the KDEMS spectrum is analyzed. Then the localization and counting performance of the combined GCC-PS, KDE-PS, and KDEMS-PS is further discussed.The sound sources are taken from 16 pure speeches composed of 8 male and 8 female in the TIMIT database [31]. Each speech segment sampled with fs=16kHz lasts 2 seconds, and has the same average power through preprocessing. The room is 7.8m×8.1m×3m whose x-y plane diagram is shown in Figure 2. A pair of omni-directional microphones ma, mb with distance dmic parallel to the x-axis is located at the center of the room marked with o, where dmic is set to be an multiple of 0.5λmin with a positive integer u as the multiple factor. N sound sources are distributed on a semicircle with o as the centroid, where direction of arrival (DOA) of the n-th sound source θn∈0,180° is defined in an anti-clockwise manner with 90° being the direction perpendicular to the line connecting ma and mb. Then, the true TDOA can be obtained as τn=−dmiccosθn/c where c is the sound velocity set as 344m/s. dms is the microphone-source distance set as 3m. In Figure 2, z coordinates of the sound sources and the microphones are all set as 1.3m.Figure 2
Thex-y plane diagram of the room.When the reverberation time expressed asRT60 changes between 200ms and 500ms and dmic changes from 0.5λmin to 6λmin with 0.5λmin as the interval, 200 simulations are performed under each scenario. In each simulation, N sound sources are randomly selected from the 16 candidate segments, convolved with the room impulse response (RIR) generated by the image-source model [32, 33], and added white Gaussian noise with fixed signal-to-noise ratio (SNR) which can be obtained as(33)SNR=10log10Pa+Pb2Pv,where Pa and Pb indicate the average power of the noise-free observed signal corresponding to ma and mb, respectively, and Pv indicates the average power of the additive white Gaussian noise at the receiving end. In the simulation scenarios, multiple sound sources with equal angular intervals are used. When N is 2,4, or 6, the true DOA distribution is shown in Table 1.Table 1
Number of sourcesN versus DOA.
NDOA°260, 120430, 60, 120, 15060, 30, 60, 120, 150, 180Then the observed signal is passed through the bandpass filters between 0.2 kHz and 4 kHz [14] and transformed into the T-F coefficients with a 1024-point Hamming-weighted STFT, where the frame shift rate is 25%. Thirty consecutive frames are extracted to pool the spectrum. According to the current environment, the parameters are chosen empirically. In the Gaussian kernel function used by KDE and KDEMS, B is set to 20. τgrid in the TDOA space is uniformly set to 1×10−5s. In the localization and counting step, Γ1 is set to 0.35, and Imax is set to 10.WhenN=2 and RT60=200ms, with dmic set as 0.5λmin, 1.5λmin, 3.5λmin, or 5.5λmin, the normalized spectra generated by GCC, KDE, and KDEMS are shown in Figure 3. In each subfigure, three curves colored by red, green, and blue correspond to GCC, KDE, and KDEMS, respectively. For the intuitiveness and unity of graphic expression, a weighting factor 1/τmax is introduced to ensure the unified abscissa with different dmic. Then the true TDOAs of the sound sources are located at 0.5 and −0.5, which are marked with two vertical dashed lines.Figure 3
Normalized spectral functionΦτ of GCC, KDE, and KDEMS versus τ/τmax when RT60=200ms and (a) dmic=0.5λmin; (b) dmic=1.5λmin; (c) dmic=3.5λmin; (d) dmic=5.5λmin.
(a)(b)(c)(d)In Figure3(a), since equation (12) is fulfilled when dmic=0.5λmin, there is no spatial aliasing. Hence, the KDEMS spectrum coincides with KDE. However, the low-resolution spectra generated by the narrow dmic makes the peak too wide to produce good discrimination between different sound sources in the TDOA space. GCC has only one peak obviously deviated from the true TDOAs which make it difficult to give correct estimation.In Figure3(b), when dmic is widened to 1.5λmin, the resolution in the TDOA space improves and all the three spectra can form two peaks with the two highest amplitude near the true TDOAs. However, the main lobe corresponding to each sound source is still too wide and there is still a certain degree of deviation.In Figure3(c), when dmic is further widened to 3.5λmin, sharper peaks close to the true TDOAs make the deviation smaller than Figure 3(b). However, the simultaneously introduced spatial aliasing makes false peaks with a certain amplitude start to appear at −0.9, −0.2, 0.2, and 0.9 in GCC. Due to the suppression of high frequency components by the frequency-dependent 1/2πfk in Gaussian kernel function, the amplitudes of the false peaks in KDE are lower than GCC. Then by further MS processing, the false peaks are almost completely suppressed in KDEMS.In Figure3(d), when dmic=5.5λmin, the number of false peaks in GCC is obviously more than KDE and KDEMS, the overall amplitudes of the false peaks becomes higher than Figure 3(c). The amplitude at −0.3 in GCC is almost equivalent to the one at 0.5 corresponding to the true TDOA, which is easy to produce incorrect estimation. However, KDE and KDEMS maintain sharp peaks near the true TDOAs and still have a good suppression on the false peaks, thus have better recognition than GCC.The normalized spectra whenRT60 increases to 500ms are shown in Figure 4. The enhancement of reverberation makes the spectrum distortion more serious. In Figure 4(a), GCC, KDE, and KDEMS still face the influence of low resolution. In Figure 4(b), all the three spectra cannot estimate the TDOAs effectively because of the deviation. Compared with Figure 3(c), the amplitudes of the false peaks at −0.9, −0.2, 0.2, and 0.9 are higher in Figure 4(c), where the amplitude at −0.2, 0.2 is almost equivalent to the one at 0.5 in GCC. In Figure 4(d), the amplitudes of the false peaks at 0.3, 0.7 have exceeded the amplitude at 0.5 in GCC, the amplitude at 0.7 is close to 0.5 in KDE, while KDEMS has a relatively flat spectrum which shows that KDEMS is more robust than GCC and KDE under strong reverberation.Figure 4
Normalized spectral functionΦτ of GCC, KDE, and KDEMS versus τ/τmax when RT60=500ms and (a) dmic=0.5λmin; (b) dmic=1.5λmin; (c) dmic=3.5λmin; (d) dmic=5.5λmin.
(a)(b)(c)(d)After the spectrum analysis of a single simulation, the source localization and counting performance of GCC-PS, KDE-PS, and KDEMS-PS can be evaluated in terms of three measures defined as(34)Rrate=N^coN,Prate=N^coN^,Fscore=2Rrate⋅PrateRrate+Prate,where recall rate expressed as Rrate and precision rate expressed as Prate are used to measure miss-detections and false alarms, respectively, F-score expressed as Fscore is the combination of Rrate and Prate, and N^co denotes the number of correctly estimated sound sources. The condition for correct estimation can be expressed as(35)θ^−θ<5°,where θ^ and θ represent the DOAs corresponding to the estimated TDOA and the true TDOA, respectively.Setdmic as 3.5λmin, the average values of these three measures versus SNR when N=2 and N=6 are shown in Figures 5 and 6, respectively, where they all reach its best value at 1 and worst at 0. As SNR decreases from 30 to −10 dB, the measures of GCC-PS, KDE-PS and KDEMS-PS gradually decline. In Figure 5(a), since only two sound sources are used, the difference of Rrate between these three methods is not very obvious. But the false peaks in GCC and KDE spectra will produce more false alarms than KDEMS. With the decrease of SNR, the increase of the spectrum base further aggravates the influence of false peaks, so Prate of KDEMS-PS is significantly better than the other two methods. The situation of Fscore is similar to that of Prate.Figure 5
WhenN=2 and dmic=3.5λmin, the average measure of GCC, KDE and KDEMS versus SNR using (a) recall; (b) precision; (c) F-score.
(a)(b)(c)Figure 6
WhenN=6 and dmic=3.5λmin, the average measure of GCC, KDE and KDEMS versus SNR using (a) recall; (b) precision; (c) F-score.
(a)(b)(c)Due to the increase in the number of sound sources, the missed-detection situation has been enhanced. ThenRrate of KDEMS-PS is significantly better than the other two in Figure 6(a). False alarms still exist in Figure 6(b) but will be alleviated to a certain extent by the added sources. Prate of KDEMS-PS is still the best of the three methods. The situation of Fscore is similar to that of Rrate.It can be seen from Figures5 and 6 that when the number of sound sources is small, Prate can better distinguish the performance of different methods, and when the number of sound sources is large, Rrate is better, while Fscore is the harmonic average of Rrate and Prate to comprehensively evaluate the performance. Hence, without loss of generality, use Fscore to measure the performance of different methods. With SNR set as 20 dB, the average F-score versus dmic when RT60=200ms and RT60=500ms is shown in Figures 7 and 8, respectively.Figure 7
The average F-score of GCC-PS, KDE-PS, and KDEMS-PS versusdmic when RT60=200ms and (a) N=2; (b) N=4; (c) N=6.
(a)(b)(c)Figure 8
The average F-score of GCC-PS, KDE-PS, and KDEMS-PS versusdmic when RT60=500ms and (a) N=2; (b) N=4; (c) N=6.
(a)(b)(c)In Figure7(a), when dmic≤1.5λmin, the deviation caused by the low resolution makes N^co small, thus bringing low Fscore, then Fscore gradually rises as dmic widens. When dmic>1.5λmin, the growth of Fscore disappears, GCC-PS and KDE-PS begin to gradually decline; this is due to the increase in the amplitude of false peaks caused by spatial aliasing, which makes N^ increase. However, due to MS processing, KDEMS-PS is much higher than the other two and it can stabilize around 0.9 with no obvious downtrend when 1.5λmin<dmic≤3.5λmin, which demonstrates the capability to effectively suppress the influence of spatial aliasing.The overallFscore declines to a certain extent when N increases. When dmic=4λmin, KDEMS-PS can maintain around 0.8 and 0.6 in Figures 7(b) and 7(c), respectively, not decline as shown in Figure 7(a). This is because the increase of N requires higher resolution to distinguish different sound sources, so a higher Fscore can be brought about when dmic is appropriately increased.The rule in Figure8 is similar to that of Figure 7. The enhancement of reverberation makes N^co smaller on the one hand and N^ larger on the other hand, thus makes Fscore of different methods have a certain degree of decline. The overall performance is still KDEMS-PS >rbin KDE-PS >rbin GCC-PS. When 1.5λmin<dmic≤4λmin, KDEMS-PS can still maintain above 0.5 and 0.4 when N=4 and N=6, respectively, much higher than KDE-PS and GCC-PS, which shows the robustness of the method under strong reverberation and interference between different sound sources.


## 5. Conclusion

To correctly locate the sound sources in real scene, the number of them needs to be estimated simultaneously. Based on the kernel density estimator, a multiple sound source localization and counting method called KDEMS-PS is proposed in this paper. Divide the entire frequency band into several sub-bands and process them stage by stage, the KDEMS angular spectrum is constructed. Then PS with an updated threshold and a source merging module is combined to locate and count multiple sound sources. As shown in the computer simulation using spectrum analysis and comparison of F-score, KDEMS spectrum can effectively weaken the interference caused by false peaks and KDEMS-PS is a robust multiple sound source localization and counting method with good spatial aliasing suppression when using a wide intermicrophone distance. In the experiment, we found that the MS structure is mainly used to process sound sources with more low-frequency energy and more uniform spectrum coverage (e.g. speech). When the spectral components are concentrated in the higher frequency range (e.g. bird sound), the sub-band processing may face the problem that the unambiguous low-frequency spectrum has very weak energy, which may cause a decrease in localization performance due to insufficient peak energy at the sound source location. Whether to use some spectrum enhancement techniques or other unambiguous methods to extend the scope of application is still a question worthy of further study. Theoretically, the angular spectrum-based method based on pairwise microphones has no additional restrictions on the number of microphones or the topology. In this paper, the applied microphone array consists of only two microphones to focus on the improvement of the MS spectrum. More microphone pairs of different topologies (e.g. planar, spherical) can be added for some specific applications in the following research.



---

*Source: 6783184-2022-01-29.xml*