A Study of the Complexity and Accuracy of Direction of Arrival Estimation Methods Based on GCC-PHAT for a Pair of Close Microphones

11/28/2018 ∙ by Francois Grondin, et al. ∙ MIT 0

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.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 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.

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 interpolation-based 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.

2 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 be the sample rate (in sample/sec), the distance between two microphones, and the speed of sound . The TDOA (in sample) is given by (1), where the angle in radian represents the DOA.

(1)

The range is discretized into Q points, which leads to discrete angles , where :

(2)

The Short-Time Fourier Transform (STFT) provides a useful representation of the signals needed to compute the GCC-PHAT. In this paper, the expression stands for the STFT frame size in samples and for the hop size in samples between successive frames. We define as the STFT coefficients of microphone for the frequency bins at frame . Without loss of generality, we omit the frame index from now on for clarity. The expression then stands for the phase transformed cross-correlation spectrum at frequency bin between microphones and :

(3)

where stands for the complex conjugate and for the absolute value.

Equation 4 computes the GCC-PHAT result, denoted as , where and extracts the real part only.

(4)

where the scalar defined in (5) normalizes the transform.

(5)

The estimated DOA then corresponds to the discrete angle which maximizes the GCC-PHAT result, as shown in Eq. 6.

(6)

It is also relevant to compute the maximum GCC-PHAT magnitude, denoted here as , which provides insightful information regarding true DOA and false detections.

2.1 Matrix multiplication

The most straightforward method consists in implementing Eq. 4 via matrix multiplication, as shown in Eq. 7, where , and .

(7)

Each element of is computed offline, and given by Eq. 8. Note that the scaling factor ensures that , where the operator sets all non-diagonal elements to zero, stands for the Hermitian operator and is a identity matrix.

(8)

Although appealing in its simplicity and exactitude, this method involves a significant amount of computations, as the complexity order reaches .

2.2 Fast Fourier Transform

As an alternative to matrix multiplication, the inverse fast Fourier transform (IFFT) offers an approximation to Eq. 4 where the TDOA is rounded to the closest integer (). To reduce the discretization error, it is possible to increase the spectrum size from to frequency bins, where the bins in the range are set to zero and is usually a power of

. This zero padding in the frequency domain reduces the discretization error as it results in an interpolation operation after performing the IFFT, denoted by

. The value that corresponds to the closest TDOA is then mapped to each point , as shown in Eq’s. 9-11.

(9)
(10)
(11)

Note that the modulo operator in Eq. 11 ensures that negative values are mapped within the range . Most general purpose processors (GPDs) and digital signal processors (DSPs) compute efficiently radiix-2 IFFTs, and therefore this approach leads a complexity order of .

2.3 Quadratic interpolation

A complementary approach to the IFFT consists of performing quadratic interpolation on signal obtained in Eq. 9. The rounded TDOA corresponds to , and the rounding error is defined as , where . For each point , a parabolic curve is fit to the samples associated to and its neighbors. The expression in Eq. 12 then estimates the parabolic curve parameters , , and .

(12)

Finally, the parabolic expression provides an estimation of the fractional TDOA from the rounding error :

(13)

The complexity of computations is similar to the complexity of the IFFT-based method, with an additional parabolic fitting step for all points, which results in a total complexity of .

2.4 Singular Value Decomposition

Another alternative to matrix multiplication consists of decomposing the matrix 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 is first represented by the sum of two real matrices and , such that . This step is optional but allows using existing SVD algorithms designed for real matrices only. Given , the decomposition of with the largest singular values leads to the approximation in Eq. 14.

(14)

where , and stands for the matrix transpose operator.

To reduce the amount of computations, it is desirable to make the rank as small as possible, yet it needs to be large enough to ensure an accurate reconstruction of . This can be achieved by selecting the smallest that satisfies the condition in Eq. 15, where the operator provides the trace of the matrix, and the parameter is a small positive value that models the tolerable reconstruction error.

(15)

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.

(16)

The matrices , , and are computed offline using SVD, and then the expression in Eq. 16 is evaluated online. This leads to a complexity of , which represents a significant reduction compared to matrix multiplication when .

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 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 lookups during the mapping procedure described in (11) for the FFT-based approach.

Matrix multiplication
Fast Fourier Transform
Quadratic Interpolation
Singular Value Decomp.
Table 1: Complexity Orders of DOA Estimation Methods

3 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 arc is discretized by 181 points, which provides a resolution of one degree. The STFT is computed with frames of msecs, spaced by intervals of msecs. Microphones are spaced by cm, which matches typical distances on smart speakers. The speed of sound is chosen to match normal indoor conditions, and the reconstruction error is set to a small value to make the SVD method accurate.

Table 2: GCC-PHAT Parameters

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

dB to dB. 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 , which provides insight on the relevance of each DOA.

(17)

where stands for the Euclidean norm.

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(a)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(b)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(c)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(d)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(e)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(f)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(g)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(h)

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

RMSE
(i)
Figure 1: RMSE for the investigated methods with different reflection coefficients () and SNR – smaller is better.

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 online111http://github.com/FrancoisGrondin/gccphat), 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.

MM

FFT01

FFT02

FFT04

FFT08

FFT16

FFT32

FFT01-QI

FFT02-QI

FFT04-QI

FFT08-QI

FFT16-QI

FFT32-QI

SVD

Execution time per frame (sec)
Figure 2: Execution time per frame on a PC – less is better.

4 Conclusion

In this paper we investigate some GCC-PHAT interpolated-based 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.

References

  • [1] Albert S Bregman, Auditory scene analysis: The perceptual organization of sound, MIT press, 1994.
  • [2] E. Habets, J. Benesty, I. Cohen, S. Gannot, and J. Dmochowski, “New insights into the mvdr beamformer in room acoustics,” IEEE Trans. Audio, Speech, Language Process., vol. 18, no. 1, pp. 158, 2010.
  • [3] D.H. Johnson and D.E. Dudgeon, Array signal processing: concepts and techniques, PTR Prentice Hall Englewood Cliffs, 1993.
  • [4] L.C. Parra and C.V. Alvino, “Geometric source separation: Merging convolutive source separation with geometric beamforming,” IEEE Trans. Speech Audio Process., vol. 10, no. 6, pp. 352–362, 2002.
  • [5] J.-M. Valin, J. Rouat, and F. Michaud, “Enhanced robot audition based on microphone array source separation with post-filter,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots & Systems, 2004, vol. 3, pp. 2123–2128.
  • [6] S. Yamamoto, J.-M. Valin, K. Nakadai, J. Rouat, F. Michaud, T. Ogata, and H.G. Okuno, “Enhanced robot speech recognition based on microphone array source separation and missing feature theory,” in Proc. IEEE Int. Conf. Robotics & Automation, 2005, pp. 1477–1482.
  • [7] M.S. Brandstein and H.F. Silverman, “A robust method for speech signal time-delay estimation in reverberant rooms,” in Proc. IEEE Int. Conf. Acoustics, Speech & Signals Process., 1997, vol. 1, pp. 375–378.
  • [8] M. Cobos, A. Marti, and J.J. Lopez, “A modified SRP-PHAT functional for robust real-time sound source localization with scalable spatial sampling,” IEEE Signal Process. Letters, vol. 18, no. 1, pp. 71–74, 2011.
  • [9] B. Kwon, Y. Park, and Y.-S. Park, “Analysis of the GCC-PHAT technique for multiple sources,” in Proc. IEEE Int. Conf. Control Autom. & Systems, 2010, pp. 2070–2073.
  • [10] F. Grondin, D. Létourneau, F. Ferland, V. Rousseau, and F. Michaud, “The ManyEars open framework,” Auton. Robots, vol. 34, no. 3, pp. 217–232, 2013.
  • [11] J.-M. Valin, F. Michaud, and J. Rouat, “Robust localization and tracking of simultaneous moving sound sources using beamforming and particle filtering,” Rob. Auton. Syst., vol. 55, no. 3, pp. 216–228, 2007.
  • [12] A. Agarwal, M. Jain, P. Kumar, and S. Patel, “Opportunistic sensing with mic arrays on smart speakers for distal interaction and exercise tracking,” in Proc. IEEE Int. Conf. Acoustics, Speech & Signals Process., 2018, pp. 6403–6407.
  • [13] G. Jacovitti and G. Scarano, “Discrete time techniques for time delay estimation,” IEEE Trans. Signal Process., vol. 41, no. 2, pp. 525–533, 1993.
  • [14] M. McCormick and T. Varghese, “An approach to unbiased subsample interpolation for motion tracking,” Ultrasonic imaging, vol. 35, no. 2, pp. 76–89, 2013.
  • [15] F. Viola and W.F. Walker, “A spline-based algorithm for continuous time-delay estimation using sampled data,” IEEE Trans Ultrason Ferroelectr Freq Control, vol. 52, no. 1, pp. 80–93, 2005.
  • [16] B. Van Den Broeck, A. Bertrand, P. Karsmakers, B. Vanrumste, and M. Moonen, “Time-domain generalized cross correlation phase transform sound source localization for small microphone arrays,” in Proc. European DSP Educ. & Research Conf., 2012, pp. 76–80.
  • [17] D.L. Maskell and G.S. Woods, “The estimation of subsample time delay of arrival in the discrete-time measurement of phase delay,” IEEE Trans. Instrum. Meas, vol. 48, no. 6, pp. 1227–1231, 1999.
  • [18] K.K. Sharma and S.D. Joshi, “Time delay estimation using fractional fourier transform,” EURASIP J. Audio Speech, vol. 87, no. 5, pp. 853–865, 2007.
  • [19] J.B. Allen and D.A. Berkley, “Image method for efficiently simulating small-room acoustics,” J. Acoust. Soc. Am., vol. 65, no. 4, pp. 943–950, 1979.
  • [20] V. Zue, S. Seneff, and J. Glass, “Speech database development at MIT: TIMIT and beyond,” Speech Comm., vol. 9, no. 4, pp. 351–356, 1990.
  • [21] Matteo Frigo and Steven G Johnson, “Fftw: An adaptive software architecture for the fft,” in Proc. IEEE Int. Conf. Acoustics, Speech & Signals Process., 1998, vol. 3, pp. 1381–1384.
  • [22] S.K. Raman, V. Pentkovski, and J. Keshava, “Implementing streaming simd extensions on the pentium iii processor,” IEEE micro, vol. 20, no. 4, pp. 47–57, 2000.