{'fraction_non_alphanumeric': 0.04385042799219102, 'fraction_numerical': 0.01794563748310557, 'mean_word_length': 4.361916264090177, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Study of the Complexity and Accuracy of Direction of Arrival Estimation Methods Based on GCC-PHAT for a Pair of Close Microphones

This paper investigates the accuracy of various Generalized Cross-Correlation with Phase Transform (GCC-PHAT) methods for a close pair of microphones. We investigate interpolation-based methods and also propose another approach based on Singular Value Decomposition (SVD). All investigated methods are implemented in C code, and the execution time is measured to determine which approach is the most appealing for real-time applications on low-cost embedded hardware.


INTRODUCTION
Sound Source Localization (SSL) consists of estimating the Direction of Arrival (DOA) of a target source in space with respect to a microphone array. DOA information is useful for audio scene analysis [1] and is often used by common beamforming methods [2,3,4,5,6] to enhance the target source and reduce interferences. One of the most popular SSL method relies on the Steered-Response Power Phase Transform (SRP-PHAT) [7,8]. SRP-PHAT consists of computing the Generalized Cross-Correlation with Phase Transform (GCC-PHAT) on each pair of microphones [9]. The Fast Fourier Transform (FFT) provides an efficient way to compute GCC-PHAT [10,11], which is central for real-time systems that rely on SRP-PHAT as the number of GCC-PHATs needed is proportional to the square of the number of microphones. Using FFTs reduces the amount of computation, but also decreases the accuracy as the Time Difference of Arrival (TDOA) is rounded to the closest sample integer value. Recent smart speakers such as the Amazon Echo, Google Home and Apple HomePod are equipped with multiple microphones spaced by a few centimeters [12]. With these devices, DOA estimation can be particulary challenging as microphones are close to each other in space and the sample rate is low.
This work was supported in part by the Toyota Research Institute and by the Fonds de recherche du Québec -Nature et technologies.
To cope with the discretization artifacts, interpolation can be performed on the GCC-PHAT results [13,14,15,16]. Another approach consists of using fractional delay estimation [17]. Finally, the fractional Fourier transform [18] attempts to overcome the FFT discretization drawback.
In this paper, we compare the accuracy of some interpolationbased methods with the exact GCC-PHAT computation that involves a significant amount of computation. We propose another approach to estimate the DOA with two microphones based on Singular Value Decomposition (SVD). All investigated methods are implemented in C code, and the executation time is also measured to determine which approach is the most appealing for real-time applications on low-cost embedded hardware.

DOA ESTIMATION METHODS
The direction of arrival (DOA) of a sound source can be derived directly from the Time Difference of Arrival (TDOA) between two microphones. Let f S ∈ N be the sample rate (in sample/sec), d ∈ R >0 the distance between two microphones, and c ∈ R >0 the speed of sound . The TDOA τ ∈ R (in sample) is given by (1), where the angle in radian θ ∈ R ∩ [−π/2, π/2] represents the DOA.
The range [−π/2, π/2] is discretized into Q points, which leads to discrete angles θ q ∈ R ∩ [−π/2, π/2], where q ∈ Z ∩ [0, Q − 1]: The Short-Time Fourier Transform (STFT) provides a useful representation of the signals needed to compute the GCC-PHAT. In this paper, the expression N ∈ N stands for the STFT frame size in samples and ∆N ∈ N for the hop size in samples between successive frames. We define X l m [k] as the STFT coefficients of microphone m ∈ N ∩ [1, 2] for the frequency bins k ∈ Z ∩ [0, N/2] at frame l ∈ N. Without loss of generality, we omit the frame index l from now on for clarity. The expression X 12 [k] ∈ C then stands for the phase transformed cross-correlation spectrum at frequency bin k between microphones 1 and 2: where {. . . } * stands for the complex conjugate and | . . . | for the absolute value. Equation 4 computes the GCC-PHAT result, denoted as x 12 [q] ∈ R, where j = √ −1 and R{. . . } extracts the real part only.
where the scalar g k ∈ R >0 defined in (5) normalizes the transform.
The estimated DOA θ est then corresponds to the discrete angle θ qmax which q maximizes the GCC-PHAT result, as shown in Eq. 6.
It is also relevant to compute the maximum GCC-PHAT magnitude, denoted here as E est = x 12 [q max ], which provides insightful information regarding true DOA and false detections.
Each element w q,k ∈ C of W is computed offline, and given by Eq. 8. Note that the scaling factor ensures that diag{WW H } = I Q , where the operator diag{. . . } sets all non-diagonal elements to zero, {. . . } H stands for the Hermitian operator and I Q is a Q × Q identity matrix.
Although appealing in its simplicity and exactitude, this method involves a significant amount of computations, as the complexity order reaches O(QN ).

Fast Fourier Transform
As an alternative to matrix multiplication, the inverse fast Fourier transform (IFFT) offers an approximation to Eq. 4 where the TDOA τ q is rounded to the closest integer (⌊τ q ⌉ ∈ Z). To reduce the discretization error, it is possible to increase the spectrum size from (N/2 + 1) to (N i/2 + 1) frequency bins, where the bins in the range [N/2 + 2, N i/2 + 1] are set to zero and i is usually a power of 2. This zero padding in the frequency domain reduces the discretization error as it results in an interpolation operation after performing the IFFT, denoted by y 12 ∈ R N ×1 . The value that corresponds to the closest TDOA is then mapped to each point q, as shown in Eq's. 9-11. y 12 = IFFT{X 12 } (9) map : Note that the modulo operator in Eq. 11 ensures that negative values are mapped within the range [0, iN − 1]. Most general purpose processors (GPDs) and digital signal processors (DSPs) compute efficiently radiix-2 IFFTs, and therefore this approach leads a complexity order of O(iN log iN ).

Quadratic interpolation
A complementary approach to the IFFT consists of performing quadratic interpolation on signal y 12 obtained in Eq. 9. The rounded TDOA corresponds toτ q = ⌊iτ q ⌉, and the rounding error is defined as ∆τ q ∈ R, where ∆τ q = iτ q − ⌊iτ q ⌉. For each point q, a parabolic curve is fit to the y 12 samples associated toτ q and its neighbors. The expression in Eq. 12 then estimates the parabolic curve parameters a q ∈ R, b q ∈ R, and c q ∈ R.
Finally, the parabolic expression provides an estimation of the fractional TDOA from the rounding error ∆τ q : The complexity of computations is similar to the complexity of the IFFT-based method, with an additional parabolic fitting step for all Q points, which results in a total complexity of O(iN log iN + Q).

Singular Value Decomposition
Another alternative to matrix multiplication consists of decomposing the matrix W into a different set of bases. As the microphone distance gets smaller, the bases that compose this matrix become more linearly dependent. This implies that a smaller set of orthogonal bases obtained via Singular Value Decomposition (SVD) can span the same space. The matrix W is first represented by the sum of two real matrices W R ∈ R Q×(N/2+1) and W I ∈ R Q×(N/2+1) , such that W = W R + jW I . This step is optional but allows using existing SVD algorithms designed for real matrices only. Given α ∈ {R, I}, the decomposition of W α with the K α largest singular values leads to the approximation in Eq. 14.
To reduce the amount of computations, it is desirable to make the rank K α as small as possible, yet it needs to be large enough to ensure an accurate reconstruction of W α . This can be achieved by selecting the smallest K α that satisfies the condition in Eq. 15, where the operator Tr{. . . } provides the trace of the matrix, and the parameter δ is a small positive value that models the tolerable reconstruction error.
The decomposition obtained in Eq. 14 is then substituted in the initial matrix multiplication in Eq. 7, which leads to Eq. 16, where ℜ{. . . } and ℑ{. . . } capture respectively the real and imaginary parts of the corresponding expression.
The matrices U R , T R , U R and U I are computed offline using SVD, and then the expression in Eq. 16 is evaluated online. This leads to a complexity of O((K R +K I )(N +Q)), which represents a significant reduction compared to matrix multiplication when (K R + K I ) ≪ (N + Q). Table 1 presents a summary of all methods and their corresponding complexity order. It is interesting to note that the number of potential source positions Q affects the complexity order of Matrix Multiplication, Quadratic Interpolation and Singular Value Decomposition, as opposed to the Fast Fourier Transform method. Although no arithmetic is involved, there are Q lookups during the mapping procedure described in (11) for the FFT-based approach.

RESULTS
We perform experiments to analyze the accuracy of the previously described GCC-PHAT methods: Matrix Multiplication (MM), Fast Fourier Transform with interpolation rate from 1 to 32 (FFT01, FFT02, FFT04, FFT08, FFT16 and FFT32), Fast Fourier Transform with interpolation rate from 1 to 32, followed by Quadratic Interpolation (FFT01-QI, FFT02-QI, FFT04-QI, FFT08-QI, FFT16-QI, and FFT32-QI), and Singular Value Decomposition (SVD). Table 2 lists the parameters used with the investigated methods. The search space modeled by an 180 • arc is discretized by 181 points, which provides a resolution of one degree. The STFT is computed with frames of 32 msecs, spaced by intervals of 10 msecs. Microphones are spaced by d = 5 cm, which matches typical distances on smart speakers. The speed of sound c is chosen to match normal indoor conditions, and the reconstruction error δ is set to a small value to make the SVD method accurate. Simulations are conducted to measure the accuracy of the proposed method. The two microphones are positioned randomly in rooms of various sizes. Three categories of rooms are investigated: small (between 5m x 5m x 3m and 10m x 10m x 5m), medium (between 10m x 10m x 3m and 20m x 20m x 5m) and large (between 20m x 20m x 5m and 20m x 20m x 10m). For each configuration, the room category is chosen randomly, and then the dimensions are generated randomly within the corresponding size range. The reflection coefficients are set to 0, 0.3 and 0.6 to investigate different reverberation levels. Reverberation is modeled with Room Impulse Responses (RIRs) generated with the image method [19]. Sound segments from the TIMIT dataset are then convolved with the generated RIRs [20]. Diffuse white noise is added to each channel, for a signal-to-noise ratio (SNR) that varies from 10dB to 40dB. A total of 1000 different configurations are generated.
The Root Mean Square Error (RMSE) is calculated in Eq. 17 by summing all DOAs for each room and speech signal, with the estimated DOA weighted with the energy E est , which provides insight on the relevance of each DOA.
where . . . 2 stands for the Euclidean norm. Figure 1 shows the RMSE with multiple reverberation levels and SNR values. Note how the accuracy decreases as the reverberation increases (when the reflection coefficient β increases) and as the SNR decreases. For all cases, the MM method provides the best performance. Moreover, the FFT{02,04,08,16,32}-QI matches the MM performance, followed closely by the SVD approach. It is interesting to note the reduction in accuracy when the quadratic interpolation alone on non-interpolated FFT results (FFT01-QI) is used. Note how the FFT interpolation also reduces the RMSE, but requires a significant factor (FFT32) to achieve a performance level similar to MM, whereas other interpolation rates increase the RMSE (FFT{01,02,04,08,16}). All methods are implemented in C to minimize execution time on an Intel Xeon E5-1620 3.70GHz (source code is available online 1 ), and use the FFTW library to optimize the computations of FFTs [21]. Figure 2 shows the execution time per frame for all methods. Among methods that minimize RMSE, FFT02-QI and SVD also minimize the execution time. As expected, the MM approach is the most expensive in terms of computational requirements, and the FFT interpolation increases rapidly the execution time. Note that, as opposed to the FFT from the FFTW library, the matrix multiplication in the SVD method is performed elementwise, and not optimized for special sets of instruction such as Streaming SIMD Extensions (SSE) [22]. Using these special instructions should reduce the SVD method execution time.

CONCLUSION
In this paper we investigate some GCC-PHAT interpolatedbased methods and a SVD-based method for accurate DOA  estimation for physically close microphones. Results demonstrate that interpolation by increasing the size of the FFT and quadratic interpolation must be combined together to perform accurate DOA estimation. It also demonstrates that the proposed SVD transform also offers good accuracy. Finally, doubling the size of the FFT with quadratic interpolation and the SVD method both offer the best trade-off between accuracy and computational load. As future work, we could evaluate the performance for various distances between a pair of microphones, and then choose the optimal GCC-PHAT building blocks with microphone arrays of arbitrary shapes.