Multiple Sound Source Localization with SVD-PHAT
Francois Grondin James Glass



Differentially private sub-Gaussian location estimators
    Marco Avella-Medina Columbia University, Department of Statistics, marco.avella@columbia.edu Victor-Emmanuel Brunel ENSAE ParisTech, Department of Statistics, victor.emmanuel.brunel@ensae.fr
    December 30, 2023
==================================================================================================================================================================================================




This paper introduces a modification of phase transform on singular value decomposition (SVD-PHAT) to localize multiple sound sources.
This work aims to improve localization accuracy and keeps the algorithm complexity low for real-time applications.
This method relies on multiple scans of the search space, with projection of each  low-dimensional observation onto orthogonal subspaces.
We show that this method localizes multiple sound sources more accurately than discrete SRP-PHAT, with a reduction in the Root Mean Square Error up to 0.0395 radians.


Index Terms: multiple sound source socalization, srp-phat, svd-phat, direction of arrival




§ INTRODUCTION



The cocktail party effect consists of the ability to focus on a specific conversation in a noisy environment.
While humans can usually perform this task efficiently, distant speech processing remains challenging for automatic speech recognition (ASR) systems <cit.>.
To improve ASR performances, it is common to use a beamformer with multiple microphones as a preprocessing step to enhance the corrupted speech signal <cit.>.
Some beamforming methods, such as the delay and sum and the minimum variance distortionless response (MVDR) <cit.>, require the target source direction of arrival (DOA).
On the other hand, methods like geometric sound separation require both the target and interference sources direction of arrival <cit.>.
It is therefore desirable to estimate the direction of arrival of multiple sound sources.

High resolution methods such as Multiple Signal Classification (MUSIC) <cit.> and the Estimation of Signal Parameters via Rotational Invariance Technique (ESPRIT) <cit.> were initially designed for narrowband signals, and subsequently adapted to broadband signals such as speech <cit.>.
However, MUSIC-based methods involve online computations of eigenvectors, which makes real-time implementation challenging on low-cost embedded hardware.
On the other hand, ESPRIT-based techniques require significantly less computations, but need twice as many sensors as MUSIC to perform with similar performance, which is problematic for microphone arrays with few sensors.

Alternatively, the Steered-Response Power Phase Transform (SRP-PHAT) robustly estimates the direction of arrival of speech sources and can be computed with low-cost embedded hardware <cit.>.
SRP-PHAT relies on the Generalized Cross-Correlation with Phase Transform (GCC-PHAT) between each pair of microphones.
The Fast Fourier Transform is often used to speed up the computation of GCC-PHAT, but this also reduces localization accuracy.
This discrete SRP-PHAT approach can localize many sound sources by scanning the search space multiple times, and nulling the GCC-PHAT region related to each found DOA <cit.>.
Hierarchical search also reduces the number of lookups in memory <cit.>.
The discrete SRP-PHAT approach however relies on rounded TDOA values, which may reduce the localization accuracy.

Alternatively, Cai et al. propose using multiple subbands to individually localize one sound source per band <cit.>.
Similarly, it is possible to localize multiple speech sources based on their distinct pitch values <cit.>.
Pavlidi et al. introduce a technique to identify single-source zones in the time-frequency range and generate a histogram to count and localize multiple sound sources <cit.>.
However, these methods rely on narrow bands to localize sound sources, which makes localization more sensitive to reverberation.
On the other hand, localization can also exploit interesting properties of microphone arrays with symmetrical geometries.
For instance, wavefield decomposition enables localizing multiple sound source with spherical arrays <cit.>.
Similarly, low-complexity multiple sources localization is possible in 2-D with circular arrays <cit.>.
These methods offer interesting performance, but rely on a specific microphone array geometry, which restricts their scalability.

We recently proposed a new method called SVD-PHAT that relies on singular value decomposition to map the observations to a small subspace, and then uses a nearest neighbor search algorithm like a k-d tree to find the DOA <cit.>.
This single source localization method is appealing as it preserves exact SRP-PHAT accuracy while greatly reducing the computational complexity, and can adapt to microphone array with arbitrary shapes.
In this paper, we extend SVD-PHAT to localize multiple sound sources.



§ SRP-PHAT



We first introduce SRP-PHAT with rounded TDOA that allows efficient localization of multiple sound sources with arbitrary array shapes.
Let X^l_m[k] ∈ℂ be the Short Time Fourier Transform (STFT) coefficients, where N ∈ℕ and Δ N ∈ℕ stand for the frame and hop sizes in samples, respectively, and k ∈{0,1,…,N/2}, m ∈ℳ = {1,2,…,M} and l ∈ℕ stand for the frequency bin, microphone and frame indexes, respectively.
The cross-correlation X^l_i,j[k] for each microphone pair (i,j) ∈𝒫 = {(x,y) ∈ℳ^2: x < y} is obtained with the following recursive estimation with α∈ [0,1]:


    X^l_i,j[k] = (1-α) X^l-1_i,j[k] + α X^l_i[k] (X^l_j[k])^*


where {…}^* stands for the complex conjugate.
For clarity, the frame index l is omitted in this paper without loss of generality.
The phase transformed spectrum X̂_i,j[k] ∈ℂ for each microphone pair is then obtained in (<ref>), where |…| stands for the absolute value.


    X̂_i,j[k] = X_i,j[k]/|X_i,j[k]|


The generalized cross correlation with phase transform (GCC-PHAT) for each pair of microphone and TDOA τ∈ℝ is given in (<ref>), where W[k,τ] = exp(2π√(-1)kτ/N).


    x_i,j[τ] = ∑_k=0^N/2X̂_i,j[k]W[k,τ]


The TDOA τ_i,j,q∈ℝ (in samples) corresponds to the difference between the direction of arrival (DOA) from a source 𝐬_q ∈{𝐱∈ℝ^3: ‖𝐱‖_2 = 1} to microphone i at position 𝐫_i∈ℝ^3, and the DOA between the same source and another microphone j at position 𝐫_j∈ℝ^3, scaled with the speed of sound in air c ∈ℝ^+ (in m/sec) and the sample rate (f_S ∈ℕ):


    τ_i,j,q = f_S/c(𝐫_j - 𝐫_i) ·𝐬_q


where {·} stands for the dot product.

It is common to discretize τ_i,j,q by rounding to the closest integer (denoted as ⌊τ_i,j,q⌉∈ℤ), and compute the GCC-PHAT in (<ref>) using an Inverse Fourier Transform (IFFT) for all τ = n ∈𝒩 = {0, 1, …, N-1}.
The expression Y_q∈ℝ is then obtained as follows for every possible DOA 𝐝_q, where q ∈𝒬 = {1, 2, …, Q}::


    Y_q = ∑_(i,j)∈𝒫x_i,j[⌊τ_i,j,q⌉ N]


The index of the most likely DOA then corresponds to:


    q^* = _q ∈𝒬{Y_q}


Once the DOA at index q^* is found, we can remove its contribution from the current observations and perform a new scan to detect other active sound sources.
A naive approach consists in nulling the expression Y_q^* and some of its closest neighbors, and then scan for a new maximum value.
However, this approach ignores the possible side lobes generated by the found source, and these may lead to false positives in the next scan iteration.
To address this issue, a popular solution consists in nulling some regions in the GCC-PHAT results instead, recompute Y_q ∀ q ∈𝒬 with (<ref>), and then find a new maximum as in (<ref>).
For each DOA index q, we define a subset of DOA indexes 𝒬_q = {x ∈{1, 2, …, Q}: arccos(𝐬_x ·𝐬_q) ≤Δθ} that gathers DOAs close in space to the DOA 𝐬_q, where Δθ is a user-defined parameter that stands for the maximum angle difference.
We then define the set 𝒯_i,j,q = {τ_i,j,q^min, τ_i,j,q^min+1, …, τ_i,j,q^max-1, τ_i,j,q^max} that contains all TDOAs related to the DOA 𝐬_q and its closest neighbors, for the microphone pair (i,j), where:


    τ_i,j,q^min = ⌊min_p ∈𝒬_q{τ_i,j,p}⌋ and τ_i,j,q^max = ⌈max_p ∈𝒬_q{τ_i,j,p}⌉


and the GCC-PHAT values in range of 𝒯_i,j,q^* are then set to zero for all pairs.

Algorithm <ref> summarizes how SRP-PHAT can be adapted to localize R multiple sources.
At each scan r ∈ℛ = {1, 2, …, R}, the GCC-PHAT values are updated, and the following scans are thus performed without the contribution of the source recently found.
The expressions 𝐝_r and e_r stand for the DOA and energy level found at scan r.



Although appealing as it relies on an efficient implementation of GCC-PHAT with IFFTs, this approach relies on discrete cross-correlation results, which reduces the accuracy.
We therefore propose to adapt SVD-PHAT to estimate the direction of arrival (DOA) of multiple sound sources with more accuracy.



§ SVD-PHAT



To define the SVD-PHAT method, it is convenient to start from SRP-PHAT in matrix form.
Let us define the vector 𝐗_i,j∈ℂ^(N/2+1)× 1 for the microphone pair (i,j) ∈𝒫 that holds the phase normalized cross-correlation coefficients for all bins k ∈{0, 1, …, N/2}:


    𝐗_i,j = [
            [   X̂_i,j[0]   X̂_i,j[1]           ⋯ X̂_i,j[N/2];             ]]^T


where {…}^T stands for the transpose operator.

A single vector 𝐗∈ℂ^P(N/2+1)× 1 then holds all these vectors:


    𝐗 = [
            [   (𝐗_1,2)^T   (𝐗_1,3)^T           ⋯ (𝐗_M-1,M)^T;             ]]^T


For each pair of microphones (i,j), all coefficients W[k,τ] ∈ℂ are concatenated in a matrix 𝐖_i,j∈ℂ^Q× (N/2+1):


    𝐖_i,j = [
        [   W[0,τ_1,i,j]              ⋯ W[N/2,τ_1,i,j];   W[0,τ_2,i,j]              ⋯ W[N/2,τ_2,i,j];              ⋮              ⋱              ⋮;   W[0,τ_Q,i,j]              ⋯ W[N/2,τ_Q,i,j];                ]]


The supermatrix 𝐖∈ℂ^Q× P(N/2+1) then holds all the matrices 𝐖_i,j ∀ (i,j) ∈𝒫:


    𝐖 = [
    [   𝐖_1,2   𝐖_1,3       ⋯ 𝐖_M-1,M ]]


Finally, the vector 𝐘∈ℝ^Q× 1 holds the results Y_q ∀ q ∈{1, 2, …, Q}, where {…} returns the real part:


    𝐘 = [
            [ Y_1   … Y_Q ]]^T = {𝐖𝐗}


The supermatrix 𝐖 can be estimated with SVD of rank K ∈{1, 2, …, K_max}, with K_max = max{Q,P(N/2+1)} and where 𝐔∈ℂ^Q× K, 𝐒∈ℂ^K× K and
𝐕∈ℂ^P(N/2+1)× K:


    𝐖≈𝐔𝐒𝐕^H


The rank K corresponds to the minimum value for which the following condition holds, where δ∈ (0,1) is a user-defined parameter that stands for the reconstruction tolerable error:


    {𝐒𝐒^T}≥ (1 - δ){𝐖𝐖^H}


where {…} stands for the trace of the matrix.

The vector 𝐙∈ℂ^K × 1 then results from the projection of the observations 𝐗 in the K-dimensions subpace:


    𝐙 = 𝐕^H𝐗


Similarly, we define the dictionary 𝐃∈ℂ^Q × K, made of the vectors 𝐃_q ∈ℂ^1 × K ∀ q ∈{1, 2, …, Q}:


    𝐃 = 𝐔𝐒 = [
        [ (𝐃_1)^T (𝐃_2)^T       … (𝐃_Q)^T;         ]]^T


As explained in <cit.>, the DOA index then corresponds to q^*, obtained as follows:


    q^* = _q ∈𝒬{{𝐃_q ·𝐙^H}}


which can be converted into the following nearest neighbor problem with an algorithm such as k-d tree:


    q^* = _q ∈𝒬{‖𝐃̂_q-𝐙̂^H‖_2^2}


where 𝐃̂_q = 𝐃_q / ‖𝐃_q ‖_2 and 𝐙̂ = 𝐙 / ‖𝐙‖_2.

Intuitively, we would like to remove the component in 𝐙 that spans the space spanned by (𝐃_q^*)^*, and then perform a new scan to find another source.
We thus define the vector 𝐯_r ∈ℂ^1 × K as follows:


    𝐯_r = (𝐃_q^*)^*


The Gram-Schmidt process then makes the current vector 𝐯_r at scan r orthogonal to all the vectors previously found (𝐮̂_n ∀ n ∈{1,2,…,r-1}), and generates 𝐮_r:


    𝐮_r = 𝐯_r - ∑_n=1^r-1(𝐮̂_n ·𝐯_r)𝐮̂_n


which is then normalize to have a unit norm:


    𝐮̂_r =𝐮_r/‖𝐮_r ‖_2


Finally, the current observation 𝐙 is projected in the subspace orthogonal to 𝐮̂_r to remove the current contribution of the source previously found:


    𝐙' = 𝐙 - (𝐮̂_r ·𝐙)𝐮̂_r


Algorithm <ref> summarizes these steps for SVD-PHAT.
This approach is appealing as it involves R k-d tree search instead of computing R times Y_q ∀ q ∈𝒬 as in (<ref>), which reduces the algorithm complexity.





§ RESULTS



We investigate three different microphone array geometries: a 1-D linear array, a 2-D planar array and a 3-D array.
The microphones xyz-positions with respect to the center of the array are given in cm in Table <ref>.



Simulations are conducted to measure the accuracy of the proposed method and compare it to the SRP-PHAT approach discretized with IFFTs.
The microphone array is positioned and rotated randomly in a 10m× 10m× 3m rectangular room, with a minimum distance of 0.5m from the walls, ceiling and floor.
The target sources are also positioned randomly in the room, and the random setup ensures a minimum angle difference of 30^∘ between each source, a distance of at least 0.5m between each source and the center of the microphone array, and a distance of at least 0.5m between each source and the walls, ceiling and floor.
For each configuration, the room reverberation is modeled with Room Impulse Responses (RIRs) generated with the image method <cit.>, where the reverberation time (RT60) is sampled randomly in the uniform interval between 200 and 500 msecs, which corresponds to the levels previously used in <cit.>.
Sound segments selected randomly from the TIMIT dataset <cit.> are normalized to have the same energy levels, and are convolved with the generated RIRs.
For each type of array (1-D, 2-D and 3-D) and number of active sources (1, 2 and 3), we perform 1000 simulations.

Table <ref> introduces the parameters used with SRP-PHAT and SVD-PHAT.
The sample rate f_S captures all the frequency content of speech (including wideband fricatives that contain relevant localization information), and the speed of sound c matches typical indoor conditions at room temperature.
The frame size N analyze speech segments of 32 msecs, and the hop size provides a 50% overlap.
The DOAs are scanned on a unit sphere generated recursively from a tetrahedron, for a total of 2562 points, as in <cit.>.
Moreover, the cross-correlation adaptation rate α estimates the sound statistic over the past 400 msecs.
In the case of SRP-PHAT, the maximum angle difference Δθ corresponds to 0.1745 radians to null the current source within a region of 10^∘.
In the specific case of SVD-PHAT, the user-defined parameter is set to δ = 10^-5 as in <cit.>.



For a 1-D array with all microphones on the x-axis, the spatial resolution is limited to an arc that goes from 0^∘ to 180^∘ in the xy-plane.
The 3-D DOA is therefore projected to this subspace as follows:


    f_1(𝐱) = [cos(g(𝐱)), sin(g(𝐱)), 0]


where:


    g(𝐱) = {(𝐱)_x,√((𝐱)_y^2+(𝐱)_z^2)}


Similarly, a 2-D array that spans the xy-plane allows DOA estimation on a half hemisphere oriented in the z-axis, and therefore all DOAs are projected to the positive z-axis:


    f_2(𝐱) = [(𝐱)_x, (𝐱)_y, |(𝐱)_z|]


Finally, for a 3-D array, the DOA can span the full space:


    f_3(𝐱) = 𝐱


For each speech source t ∀ {1,2,…,T}, the minimum angle difference (in radians) between the theoretical DOA and all estimated DOAs at frame l ∈{1, 2, …, L} is given as follows:


    ϕ_t^l = min_r ∈ℛ{arccos(f_β(𝐝^l_r) · f_β(𝐜_t))}


where β∈{1,2,3} matches the array geometry.
The goal is therefore to have at least one DOA estimation that matches each speech source true DOA.

In the proposed experiments, the number of sources varies with T = {1,2,3}, and the number of scans R matches this number.
The root mean square error (RMSE) in rad for a simulation therefore corresponds to:


    RMSE = √(1/LT∑_l=1^L∑_t=1^T(ϕ^l_t)^2)


Figure <ref> shows the estimated DOAs obtained with SRP-PHAT and SVD-PHAT for a 1-D array with three speech sources located at -1.2192 rad, -0.4335 rad and 0.4015 rad, and a reverberation time (RT60) of 238 msecs.
In this example, the SRP-PHAT method fails to detect the source at -0.4335 rad at different times, whereas SVD-PHAT detects this source most of the time.
The RMSEs of SRP-PHAT and SVD-PHAT correspond to 0.3009 rad and 0.2027 rad, respectively, which indicates that SVD-PHAT outperforms SRP-PHAT in this specific example.



The RMSE gap between both SRP-PHAT and SVD-PHAT is however usually smaller than the one shown in Figure <ref>.
To better compare both methods, Table <ref> shows the mean of all RMSEs for the 1000 simulations with each configuration.
In all cases, the proposed multiple source SVD-PHAT reduces the RMSE compared to the discrete SRP-PHAT, but with a smaller gap that oscillates between 0.0244 rad and 0.0395 rad.
It is interesting to note that for multiple sound sources (T>1), the RMSE increases rapidly.
This is expected as multiple active sources partially overlap each other in the time-frequency domain, which makes localization more challenging.
The best improvement for multiple sources occurs for the 3-D array when two sources are active, with a reduction in the RMSE of 0.0395 rad.





§ CONCLUSION



This paper extends SVD-PHAT for multiple sound source localization.
This technique outperforms the discrete SRP-PHAT approach in terms of accuracy, while preserving the low complexity of the original SVD-PHAT.
On average, the reduction in the RMSE varies between 0.0244 and 0.0395 radians, and the best improvement is observed for an array that spans 3-D space with two simultaneous speech sources.

In future work, we will investigate alternatives to k-d tree search to address the curse of dimensionality during the nearest neighbor search <cit.>.
The method could also be extended to deal with speed of sound mismatch and the near-field effect.
Microphone directivity could also be combined with SVD-PHAT to make the propagation model more realistic <cit.>.
The sound source tracking method proposed in <cit.> could also be combined to SVD-PHAT to estimate the number of sound sources and track their positions over time.
Finally, it would be interesting to implement SVD-PHAT in C code for easy deployment on real-time embedded systems.







IEEEtran


