TDOA Matrices: Algebraic Properties and their Application to Robust Denoising with Missing Data
Abstract
Measuring the Time delay of Arrival (TDOA) between a set of sensors is the basic setup for many applications, such as localization or signal beamforming. This paper presents the set of TDOA matrices, which are built from noisefree TDOA measurements, not requiring knowledge of the sensor array geometry. We prove that TDOA matrices are ranktwo and have a special SVD decomposition that leads to a compact linear parametric representation. Properties of TDOA matrices are applied in this paper to perform denoising, by finding the TDOA matrix closest to the matrix composed with noisy measurements. The paper shows that this problem admits a closedform solution for TDOA measurements contaminated with Gaussian noise which extends to the case of having missing data. The paper also proposes a novel robust denoising method resistant to outliers, missing data and inspired in recent advances in robust lowrank estimation. Experiments in synthetic and real datasets show significant improvements of the proposed denoising algorithms in TDOAbased localization, both in terms of TDOA accuracy estimation and localization error.
I Introduction
Time delay of arrival (TDOA) estimation is an essential preprocessing step for multiple applications in the context of sensor array processing, such as multichannel source localization [1], selfcalibration [2] and beamforming [3]. In all cases, performance is directly related to the accuracy of the estimated TDOAs [4].
Estimating TDOA in noisy environments has been subject of study during the last two decades [5, 6, 7], and is still an active area of research, benefiting from current advances in signal processing and optimization strategies [8, 9, 10, 11].
Typically, the TDOA between a single pair of sensors is obtained by measuring the peak of the generalized crosscorrelation (GCC) of the received signals on each sensor [12], which are assumed to be generated from a single source. Many factors, such as the spectral content of the signal, multipath propagation, and noise contribute to errors in the estimation of the TDOA.
Given a set of sensors, TDOA measurements can be obtained for every possible pair of sensors. This is commonly known as the full TDOA set or spherical set [13]. This paper studies how to reduce noise and errors from the full TDOA set. The intuition behind this denoising is to exploit redundancy of the full TDOA set. For sensors, the full set of measurements can be represented by values, which are referred to as the nonredundant set. This problem has been studied in the past, showing that one can optimally obtain the nonredundant set when TDOA measurements are contaminated with additive Gaussian noise. This is known as the GaussMarkov estimator [14]. However, in more realistic scenarios errors are not Gaussian and some of the TDOA measurements may contain outliers. In these cases the GaussMarkov estimator performs poorly.
This paper presents the TDOA matrix, which is created by the arrangement of the full TDOA set inside a skewsymmetric matrix, and studies the algebraic properties of this matrix, showing that it has rank and a SVD decomposition with degrees of freedom. Such matrices have been previously defined in the literature [15], but their properties and applications have not been studied until now.
These algebraic properties are used in this paper to perform denoising under different scenarios, which include the presence of missing TDOA measurements and outliers. These denoising algorithms are tested in the context of speaker localization with microphone arrays, using synthetic and publicly available real datasets. Our denoising algorithms are able to recover accurate TDOA values for high rates of missing data and outliers, significantly outperforming the GaussMarkov estimator in those cases. All the proposed methods don’t require knowledge of the sensor positions, so that they can also be used for calibration [2].
The main contributions of this work are threefold: i) Definition of the algebraic properties of TDOA matrices. ii) A closedform solution for TDOA denoising for Gaussian noise and the presence of missing data. iii) Novel robustdenoising methods for handling additive correlated noise, outliers and missing data.
Ia Notation
Real scalar values are represented by lowercase letters (e.g. ). Vectors are by default arranged columnwise and are represented by lowercase bold letters (e.g. ). Matrices are represented by uppercase bold letters (e.g. ). Lowercase letters are reserved to define vector and set sizes (e.g. vector is of size ), and denotes transpose of vector . Calligraphic fonts are reserved to represent generic sets (e.g. ) or functions applied to matrices (e.g. ). The norm will be written by default as for simplicity, and is the Frobenius norm, while is reserved to represent absolute values of scalars. The norm of a matrix, written , is defined as the number of nonzero elements of the matrix. is the Hadamard product between and , defined as the entrywise multiplication of the corresponding matrices. is the trace function.
We also define the normalized unitary vector as , and the null vector as , both of them having size . Finally, is a matrix with all elements equal to , is a diagonal matrix where its main diagonal is the vector , and is the identity matrix.
IB Paper Structure
The rest of the paper is distributed as follows. Section II describes the related work and Section III the problem statement. In section IV TDOA matrices are described along with a derivation of their properties. TDOA denoising in Gaussian noise case is addressed in section V, also providing a closedform solution. In sections VI and VII we propose novel algorithms for robust handling of noise and missing data, respectively. Section VIII combine the proposals of the previous two sections into an unified algorithm. We also provide an extensive experimentation to validate the proposed algorithms using both synthetic (Section IX) and real data (Section X). Finally, conclusions are drawn in section XI.
Ii Related Work
TDOA estimation is an essential first step for multiple applications related to 1) localization, 2) selfcalibration and 3) beamforming (among others):
1) Localization: widely used in radar, sonar and acoustics, since no synchronization between the source and sensor is needed. The TDOA information is combined with knowledge of the sensors’ positions to generate a Maximum Likelihood spatial estimator made from hyperboloids intersected in some optimal sense. A linear closedform solution of the former problem, valid when the TDOA estimation errors are small, is given in [16].
2) Selfcalibration: since knowing the position of sensors is mandatory for localization techniques, some strategies have been also proposed in order to calibrate them using only TDOA measurements. In [17, 2], the TDOA problem is converted in a Time of Arrival (TOA) problem estimating the departure time of signals. Then, selfcalibration techniques for TOA can be employed. The main drawback of this approach is that the conversion step from TDOA to TOA is very sensitive to outliers and correlated noise.
3) Beamforming: precise TDOA estimations is also critical for beamforming techniques and its applications. In [3], for example, additional steps are proposed for selecting the appropriate TDOA value among the correlation peaks, and also dealing with TDOA outliers. These steps include a Viterbi decoding based algorithm which maximizes the continuity of the TDOA estimations in several frames. However, the TDOA selection criteria is just based on their distance to surrounding TDOA values and their GCCPHAT values, thus not attempting to benefit from the actual redundancy of the TDOA measurements.
Hence, an accurate estimation of TDOA is essential for a good performance of any of the former applications based on these measurements.
Typically, when only two sensors are employed, the peak of the generalized crosscorrelation (GCC) function of the signals of two sensors is a good estimator for the TDOA, for reasonable noise and reverberation levels [12].
When more than two sensors are used (), there are different TDOA measurements from all possible pairs of sensors, forming the full TDOA set or spherical set [13]. However, all those TDOA measurements are redundant. In fact, usually one sensor is considered the reference sensor, and only the subset of TDOA measurements which involve that sensor are considered. That nonredundant set is the set of measurement used by the majority of TDOAbased positioning algorithms proposed in the literature [18, 16, 19, 20, 21, 22]. Nevertheless, an optimal (denoised) version of the nonredundant set can be estimated from the redundant set using a Bayesian Linear Unbiased Estimator (BLUE), also known as the GaussMarkov estimator [14].
A closedform solution for the BLUE estimator is provided in [23], also proving that it is equal to the standard least squares estimator, and that it reaches the CramerRao lower bound for positioning estimation. However, all the results in that work are based on the assumption of additive Gaussian noise, which is unrealistic in many practical applications [24], and doesn’t yield good results when correlated noise is present as a consequence, for instance, of multipath propagation. Additionally, the experimental results shown in their work are only applied to synthetic data, thus not allowing to assess the performance of their proposal in real scenarios (in section X we show the limitations of their method when evaluated on real data).
A leastsquares solution to TDOA denoising is proposed in [25]. It is based on projecting the nonredundant set of TDOA measurements into a set of “feasible” bivectors (rank 2, antisymmetric tensors) that show the same geometric properties of TDOA matrices. This denoising is also optimal for Gaussian noise but as in [24] the experimental analysis is based on simulated data and it does not cope with missing data or the presence of outliers in the TDOA measurements.
Periodicity in correlated signals, coherent noise and multipath due to reverberation are the major sources of nonGaussian error in TDOA estimation. Different approaches have been proposed to deal with them. A basic method consists in making the GCC function more robust, deemphasizing the frequencydependent weighting. The Phase Transform (PHAT) [26] is one example of this procedure which has received considerable attention as the basis of acoustic source localization systems due to its robustness in real world scenarios [27, 28]. Other approaches are based in blind estimation of multipath (room impulse response) [29] but they need a good initialization to perform well.
Some previous works have also proposed more complicated structures in order to represent TDOA redundancy, while not imposing strong assumptions on the noise distribution. In [30] a representation based on graphs allows to disambiguate if a peak in correlation was generated by the direct path or by reverberation applying an efficient search algorithm among all possible combinations. However, they do not explicitly attempt to provide improved TDOA estimations by exploiting their redundancy.
Also different matrix representations have been used in the bibliography regarding TDOA formulation. For example [31] uses a representation slightly different to the TDOA matrices we describe here, but such representation does not have the algebraic properties that TDOA matrices have, and their authors do not address an study in this sense.
So, to the best of our knowledge, there are no previous reported work dealing with improving TDOA estimations by exploiting their redundancy, while not imposing Gaussian noise restrictions, not requiring the sensor positions, and being able to deal with the presence of outliers and missing measurements (errors that will severely impact the performance of applications based on TDOA measurements). In this paper we show that TDOA matrices are a powerful tool that combined with recent advances in robust lowrank estimation, are able to generate novel solutions for these problems.
Iii Problem Statement
Hereafter, we assume only one source located at the position , and sensors synchronized between them and placed in different positions .
Given this setup, we will assume that the source is emitting an unknown signal . Then, the signal received by the sensor , , is without loss of generality, a delayed and attenuated version of (direct propagation) in addition to a signal which summarizes all the adverse effects, i.e. noise, interference, multipath, etc. Thus, , where is the time of arrival (TOA) of the signal at the sensor , being the propagation speed.
Assuming that TOA cannot be estimated directly, the time delay of arrival (TDOA) between the sensors and is estimated by correlating the received signals and (typically using the Generalized CrossCorrelation GCC [26]).
Iv TDOA matrices
In this section we define TDOA matrices, and develop their main properties. In a nutshell, given any TDOA matrix , we show that: i) is rank (Theorem 1), ii) can be decomposed as with (Lemma 1) and iii) the previous decomposition is bijective (Theorem 2).
These properties are the foundations of the denoising algorithms that we present in sections V and VI, and the missing data recovery proposal described in section VII, plus their combination described in section VIII.
Iva Definition of TDOA matrices
Definition 1.
A TDOA matrix , is a skewsymmetric matrix where the element is the time difference of arrival (TDOA) between the signals arriving at sensor and sensor :
(1) 
with
(2) 
where is the time of arrival of the signal at the sensor .
We will also express in terms of its columns as , being .
We denote as to the set of TDOA matrices of size .
Notice that there is a bijection between the full TDOA set and the corresponding TDOA matrix. Nevertheless expressing TDOA measurements as a matrix has important advantages, that we will discover throughout this article.
Note also that in the former definition, knowing the sensor array geometry is not required. For a given geometry, all the feasible TDOA matrices (those that are consistent with that particular geometry) are a subset of . Studying the properties of such subset is out of the scope of this paper and the interested reader can refer to [8, 9, 32] for further details.
On the other hand, given a particular TDOA matrix, there are infinite number of sensors geometries which match with it. Left side of Fig. 1 shows that, given a set of TOAs (,…, ) compatible with the set of TDOA measurements, the microphones can be situated in any place along the circumference (sphere in the 3D case) with center in the source (dotted lines), preserving its corespondent TOA (and therefore, its TDOA). Right side of Fig. 1 shows that there are a infinite number of TOA sets that comply with a given set of TDOA measurements.
IvB Rank of TDOA matrices
Theorem 1.
Let , then is rank 2 or 0 (trivial case).
Proof:
The matrix can be expressed as:
(3) 
where is a rank 1 matrix defined as:
(4) 
Applying the well known inequality:
(5) 
we can deduce that .
Moreover, since the rank of any skewsymmetric matrix must be even, rank 1 is not feasible. So we can conclude that, excepting the case that is the zero matrix (trivial case), the rank of is 2. This completes the proof. ∎
Note that the formerly referred as ’trivial case’ only occurs when the time of arrival is the same for all the sensors, i.e sensors are geometrically placed on a sphere with center in the source.
Rank deficiency of TDOA matrices means that their rows and columns are linearly dependent. That is consistent with the fact that, in the noisefree case, the full TDOA set can be generated from the nonredundant set using linear equations [23]. In fact, in a TDOA matrix, the column is the TDOA nonredundant set when the sensor is the reference for TDOA measurements. Hereafter, and without loss of generality, we will consider the first sensor as the reference for the non redundantset.
IvC Bijective mapping for TDOA matrices
We show next a particular representation of TDOA matrices with parameters that describes their spectral properties and also forms the base for our denoising algorithms.
Theorem 2.
Let be a dimensional hyperplane of such that , then there always exits an isomorphism between and of form:
Proof:
Theorem 2 states that a bijective linear map exists between an a dimensional subset of . Since the matrix in (4) can be rewritten as , where , we can define the following linear map:
that is clearly surjective but not injective. That is, any vector represents the same TDOA matrix. Indeed, the kernel of is the linear subspace of generated by .
Since we are looking for an isomorphism, we will restrict the domain of to ensure injectivity. It yields , where is a hyperplane of not containing . Note that is bijective as keeps the surjectivity of and, since , the function is also injective. ∎
Hereafter, we will only consider the particular case of for the hyperplane . It yields the following expression for any :
(6) 
Choosing perpendicular to is very convenient, since both it simplifies as shown in the corollary 2.1, and it also can be used for calculate the singular value decomposition (SVD) of any , as discussed in IVC1. This leads to a parametric representation of that has important properties that we will exploit later for TDOA denoising.
Corollary 2.1.
Given any , the corresponding vector , perpendicular to , can be calculated as:
(7) 
Corollary 2.2.
can also be expressed as:
(8) 
since and ,
IvC1 Singular Value Decomposition
Because is a skewsymmetric matrix of rank 2, it has the following singular value decomposition (SVD) [33, Supplementary material]:
(9) 
where and are orthonormal vectors and . Note that the SVD decomposition of is not unique. Given any orthogonal matrix , the vectors also represent a valid SVD decomposition:
(10) 
Among all possible SVD decompositions of note that theorem 2 guarantees that always exists one where . In such case SVD decomposition can be computed from (6) as stated in the following theorem:
Lemma 1.
Given , it admits the following SVD decomposition:
(11) 
Note that, according the previous lemma, not all rank 2 skewsymmetric matrices are TDOAMatrices. In fact, if and in (10) are not coplanar with (i.e. ), then the resulting matrix is rank 2 and skewsymmetric but not a TDOA matrix.
V TDOA Denoising
In this section we propose a denoising strategy to deal with Gaussian noise in the estimated TDOA measurements, deriving a closed form solution for the proposed optimization problem. This solution is also compared with the GaussMarkov Estimator.
Va Denoising Strategy
We assume now that each TDOA measurement is contaminated with uncorrelated Gaussian noise , such that . Therefore, the measured TDOA matrix is also a skewsymmetric matrix, sum of a noisefree and a skewsymmetric matrix containing noise :
(12) 
Because of the noise, and thus Theorem 1 is no longer satisfied. Consequently, the rank of may be higher than two. Nevertheless, we will show that we can take advantage of the structure of TDOA matrices in order to denoise the measured data.
For denoising, we propose finding the closest , to the measured matrix , in the sense of the Frobenius norm. This approach yields the following optimization problem:
(13) 
VB ClosedForm Solution
Theorem 3.
Problem (13) has the following closed form solution:
Proof:
From (6), the denoising problem (13) is equivalent to the following constrained convex optimization problem:
(14)  
Using the definition of Frobenius norm , and trace properties and we rewrite the cost as:
(15) 
To solve the constrained problem (14) we use the method of Lagrange multipliers, resulting in the following unconstrained equivalent:
(16) 
where is the Lagrange multiplier and
(17) 
We find extrema in (17) by taking first derivatives with respect to both and and solving the following system:
(18) 
Given that the objective function is strictly convex, the solution to the system is unique, and therefore because the proposed solution ( and ) satisfies the equations of a critical point, it is the global minimum
VC Equivalence with the GaussMarkov Estimator
By operating in (20), each element of the denoised matrix is obtained as follows:
(21) 
The closedform in (21) is identical to the one reported in [23, eq.(14)] as the GaussMarkov estimator of the TDOA measurements, so that all the properties there can be extrapolated to this work. This is not surprising as the leastsquares cost of (13) is optimal for Gaussian noise. The same denoising result was found in [25] by projecting TDOA measurements into the set of “feasible” bivectors in a leastsquares sense. Under the assumption of Gaussian noise, we can conclude that [23, 25] and our denoising result in (21) are completely equivalent.
Vi Robust TDOA Denoising
In some application scenarios, the assumption of uncorrelated white noise made in section V is fully unrealistic. In cases where the noise is correlated with the signal, measurements are prone to contain outliers in the TDOA measurements due to spurious peaks in the correlation. For such cases, a more complete model for the measured matrix is:
(22) 
where , is a skewsymmetric matrix containing Gaussian noise, much like in (12), and the new matrix models the addition of all the outliers. Since the number of outliers is usually small as compared with the number of measurements, we will assume to be sparse and unknown.
In order to denoise , we propose solving the following optimization problem, finding both matrices and :
(23)  
where is the maximum number of outliers supposed to be present in the TDOA measurements.
Robust denoising in (23) is a nonconvex optimization problem with constraints that are not even differentiable. This kind of optimization problems have been explored in Robust PCA (RPCA) [36] or robust lowrank factorizations such in GoDec [37]. Despite TDOA matrices are lowrank, these algorithms are not well suited here as they do not include all the algebraic constraints in TDOA matrices.
In order to solve (23), we propose an iterative algorithm, inspired in GoDec. It consists of an alternation method in which and are obtained in turns, with closeform solutions for these two steps (we use a subindex to denote the iteration count):
(24) 
The first subproblem of (24) is the same as our denoising problem in (13), therefore can be updated via (20). Then, is updated via entrywise hard thresholding of . Thus:
(25) 
where is an function which generates a matrix with the same size of , preserving the elements of with the largest absolute value, and making the rest of elements zero. Note that, since is skew symmetric in our application, the result provided by is also skew symmetric. The convergence to a local minimum of this algorithm is guaranteed in similar circumstances as GoDec [37], as the solutions to both subproblems in (25) are solved globally.
So, the proposed robust denoising algorithm is shown in Alg. 1.
From now on, we will refer to this algorithm as Robust DeN.
Vii Missing Data Recovery
Viia Recovery Strategy
In real scenarios, there may be situations where some of the elements of might not be available (for instance, due to communications failure) or even when they are available, there are reasons to avoid using them (for example, due to a priori knowledge of unreliable measurements, or when calculating the whole redundant set is computationally too demanding) [32]. In such cases, we want to be able to avoid some measurements, thus performing estimations when part of the values in are missing.
In this section, we address the matrix completion problem ([38, 39]) for TDOA matrices. We assume that in a measured TDOA matrix , some of its elements are unknown, and the rest are contaminated with additive Gaussian noise. We take advantage of the redundancy present in TDOA matrices to estimate a complete denoised TDOA matrix including the missing entries.
The matrix completion problem is stated as follows:
(26) 
where is a symmetric binary matrix whose element is if the TDOA between the sensor and is known, being otherwise. For convenience and without loss of generality, the elements on the main diagonal of will be set to .
Solving (26) is equivalent to finding the full TDOA matrix whose elements best fit the available elements of . Note that, , where is the result of setting the unknown elements of to zero.
ViiB ClosedForm Solution
Theorem 4.
The problem (26) has the following closed form solution: where is a diagonal matrix with as its main diagonal. is the number of missing measurements with the sensor .
Proof:
Using Corollary 2.2, problem (26) is rewritten as
(27)  
Since is the identity element of the hadamard product and is a diagonal matrix, we can rewrite (27) as:
(28)  
Operating in a similar manner to (15) we get:
(29) 
Using the identity we get:
(30) 
and finally:
(31) 
It is important to note that equations (15) and (31) are identical when there is no missing data in (i.e and ).
We use the method of Lagrange multipliers to express (27) as the following unconstrained optimization problem:
(32) 
By taking derivatives we obtain the following system:
(33) 
Since implies that , we substitute in (33) obtaining:
(34) 
where (logical not operator over all elements of ). Note that: is symmetric and, furthermore, is one of its eigenvectors.
(35) 
Therefore, if the two terms of (34) are multiplied on the right by we get:
(36) 
Then applying to (36) the fact that
(37) 
we can conclude that . Thus, the solution of (27) is:
(38) 
This completes the proof. ∎
From now on, we will refer to this algorithm as MC. It is noteworthy to comment that the matrix contains important information about the recoverability of missing data: if it is fullrank, then the solution of (26) is unique and if is rankdeficient, missing data is not recoverable uniquely without any further assumption.
Viii Robust TDOA Denoising with Missing Data
In this section we aim to combine the results of sections VI and VII, addressing the more general case in which both outliers and missing data are considered. Therefore, the problem is a combination of (23) and (26) defined as:
(40)  
In the same way as in section VI, (40) can be solved by alternatively solving the following two subproblems until convergence:
(41a)  
(41b) 
The subproblem (41a) is equivalent to the missing data problem solved in section VII but considering . Therefore, according to theorem 4, it has a closed form solution:
(42) 
Since (41b) is of the same form as the second subproblem in (24), it can also be solved by entrywise hard thresholding of .
The pseudocode shown in Alg. 2 summarizes the proposed algorithm for the general case.
Note that in line 2 the matrix can be precalculated in order to get an efficient implementation of the algorithm.
From now on, we will refer to this algorithm as Robust DeN+MC.
Ix Experiments with Synthetic Data
In this section computer simulations will be used to compare the proposed algorithms with some of the alternatives existing in the state of the art.
For evaluating the Robust DeN and Robust DeN+MC algorithms, two different metrics will be used:

The SignaltoNoiseRatio SNR [dB] of the nonredundant set referenced to the first sensor (). This is an application independent metric (where is the estimation of ), that will allow assessing the proposal improvements in the TDOA measurements per se.

The localization error, measured as the average distance between the source ground truth position and the position estimated using any given localization algorithm based on TDOA estimations (such as [16] in our case). This is an application dependent metric, that will allow assessing the actual benefits of the proposal in an example of a real task. Note that our proposal is not restricted to localization and can be used in other applications that could benefit from denoised TDOAs (such as selfcalibration or beamforming).
Ixa Experimental setup
For all the synthetic data experiments, a set of 10 sensors (which implies 45 different sensor pairs) and 1 source were randomly located. Therefore, 45 different TDOA measurements were generated per experiment, and independent Gaussian noise was added to them, using the same variance for all the measurements.
The sensor locations were uniformly distributed in a cube of 1 meter side, and the source positions were uniformly distributed in a 2 meter side cube. The propagation speed of the signal was set to 343.313 m/s. In all the experiments where it’s required, is set to . To increase the statistical significance of the results, they are provided as averages of 20 independent runs.
IxB Evaluation of Robust TDOA Denoising
In this first experiment, we evaluated the performance of the Robust DeN algorithm proposed in section VI, imposing that some TDOA values were outliers. To simulate this, we randomly chose some measurements (between 0 and 10) and replaced them with a zeromean Gaussian distributed noise, with a standard deviation of 0.1 ms. It is worth mentioning that the outlier values calculated that way are not related at all to the real TDOAs, thus being true outliers. The parameter of the proposed algorithm, which fixes the maximum number of identifiable outliers, was set to 8.
IxB1 SNR Improvements Evaluation
Fig. 2a shows the SNR values for the Robust DeN algorithm when modifying the noise standard deviation and the number of outliers, compared with that obtained by the GaussMarkov estimator (Fig. 2b), and also when only the non redundant set is used, i.e. not using the redundancy of TDOA measurements (Fig. 2c).
As predicted in section V, when no outliers are present, the performance of the Robust DeN algorithm is the same as GaussMarkov (see row 0 in Figs. 2a and 2b), hence it reaches the CramerRao Bound [23], while being much better than using no redundancy. Nevertheless, the proposed algorithm clearly outperforms the other two approaches when outliers are present in the measurements (rows 1 through 10 in the graphics of Fig. 2).
IxB2 Source Localization Improvements Evaluation
The optimized non redundant set provided by the algorithms applied in Section IXB1 were used in a localization algorithm using [16]. The average localization errors (in mm) are shown in Fig. 3. Again, the Robust DeN algorithm performs as GaussMarkov when there are no outliers, but is clearly superior when outliers are present.
It is also worth mentioning that the behaviour of the robust denoising keeps the improvements at roughly the same level for increasing number of outliers present, thus validating the ability of the algorithm to pinpoint and eliminate their presence.
IxC Evaluation of Missing Data Recovery
In this second experiment, we evaluated the capability of the MC algorithm proposed in section VII to recover missing values. For our purposes, the missing TDOA measurements were also chosen randomly but, in contrast to the previous experiment, the matrix positions of the missing measurements were known.
Fig. 4 and Fig. 5 show, respectively, the SNR values, and the localization error for the MC algorithm, when modifying the noise standard deviation and the percentage of missing TDOA values in the TDOA matrix, as compared with using only the nonredundant set (missing values were set to zero).
From the figures, it can be clearly seen that the proposed algorithm can take advantage of the knowledge about which measurements were missing, achieving even better results than when the positions of the outliers were unknown. For example, removing 50% of the full TDOA set of measurements implies 23 missing values of 45 measurements , much more than the maximum of 10 outliers evaluated in Fig. 2, while keeping good performance.
IxD Evaluation of Robust TDOA Denoising with Missing Data
In this third experiment, we evaluated the capability of the Robust DeN+MC algorithm proposed in section VIII to face both outliers and recover missing values.
To provide a wide range of evaluation scenarios, we defined: i) Two conditions related to noise, namely low and high. The former corresponds to a standard deviation of ms., and the latter to ms. ii) Two conditions related to the presence of outliers, imposing the existence of 2 or 6 outliers. iii) A variable number of missing TDOA measurements, defined as a percentage of missing TDOA values in the TDOA matrix.
In all cases, the number of outliers is fixed among the full TDOA set and then, some measurements are discarded, i.e. some of the discarded measurements may be outliers. Note that this is consistent with the real case, where it is not possible to anticipate where the outliers are.
Fig. 6 and Fig. 7 show, respectively, the SNR values, and the localization error for different algorithms, and for different evaluation scenarios.
As it can be seen in Fig. 5(a), 5(c), 6(a), and 6(c), when there are a low number of outliers (2 in this case), the best results are obtained for lower values. However, when the number of outliers increase, (Fig. 5(b), 5(d), 6(b) and 6(d), low values perform worse. So, we can conclude that must be a number as low as possible, but higher than the number of actual outliers.
Nevertheless, it is worth to observe the behaviour of Fig. 5(b), 5(d), 6(b) and 6(d) (with more outliers) when the percentage of missing data increases. It can be clearly seen that the lines corresponding to different values of are crossing among them. This seems to indicate that as the missing data percentage increases, the number of outliers that we are able to detect decreases.
Anyway, the results obtained by the Robust DeN+MC algorithm outperforms the GaussMarkov estimator, asymptotically approaching it when the noise is very high. Note also that for high values of noise, the noise and the outliers are practically indistinguishable.
X Experiments with Real Data
The aim of this section is to evaluate whether the improvements obtained in section IX using synthetic data are actually found in real environments. To do this, the proposed algorithms have been evaluated using audio recordings from the AV16.3 database [40], an audiovisual corpus recorded in the Smart Meeting Room of the IDIAP Research Institute, in Switzerland. The same metrics and localization algorithm of the previous section has been employed.
Additionally, the proposed Robust DeN algorithm is compared with other recent stateoftheart methods in section XD. In that case, we have employed the same framework used in [8], wherein many methods were already compared. In order to make the comparison as fair as possible, in this part the localization algorithm will be the same as in [8].
Xa Experimental Setup
The IDIAP Meeting Room (shown in Fig. 8) is a rectangular space containing a centrally located rectangular table, on top of which two circular microphone arrays of radius are located, each of them composed by 8 microphones. The centers of the two arrays are separated by and the origin of coordinates is located in the middle point between the two arrays. A detailed description of the meeting room can be found in [41].
The audio recordings are synchronously sampled at , and the complete database along with the corresponding annotation files containing the recordings ground truth (3D coordinates of the speaker’s mouth) is fully accessible online at [42]. It is composed by several sequences from which we are using sequence 01, with a single male speaker generating digit strings in 16 positions (which can be seen as small circles in Fig. 7(b)), distributed along the room. The sequence duration accounts for 208 seconds in total, with 823 ground truth frames.
The TDOA measurements , from which the measured TDOA matrix is built, where estimated using the highest peak of the GCCPHAT function[26].
As in a real scenario outliers are common but difficult to anticipate or enforce, the sweep over noise levels and the number of outliers that we performed with synthetic data are not feasible. Therefore, in our experiments with real data, we will only provide the SNR values and localization errors obtained after using each algorithm.
XB Evaluation of Robust TDOA Denoising
In this experiment, all the microphone pairs have been considered, hence 120 TDOA values have been computed for each frame.
In table I we show an example of the results for the Robust DeN algorithm with . As it also happened with synthetic data, in this case the proposed algorithm outperforms the GaussMarkov estimator, yielding great improvements in both SNR and localization precision. Fig. 9 shows that the selection of is not very critical in this dataset as improvements over GaussMarkov are obtained for a wide range of .
SNR (dB)  Average Localization error (mm)  

Robust DeN  27.46  354 
GaussMarkov  23.19  515 
Only nonredundant set  17.83  858 
.
These results are the baseline for the experiments with missing data described in the next subsection.
XC Evaluation of Robust TDOA Denoising with Missing Data
In the second experiment with real data, we randomly remove a set of TDOA measurements. Fig. 9 shows the obtained results. The dotted lines correspond to the performance (SNR and localization error) achieved by the GaussMarkov estimator when there are no missing measurements. The solid lines with circular marks are the results obtained by the MC algorithm described in section VII.
On the other hand, the solid lines with squared/triangular/diamond marks correspond with the results of the Robust DeN+MC algorithm presented in section VIII. The different colors/shapes indicate different values of the hyperparameter .
Fig. 9 highlights the relevance of the proposed robustdenoising algorithm (Robust DeN+MC) in realscenarios, with important improvements over its non robust version (MC): higher than 4dB absolute improvement in terms of SNR, and around 30% relative improvement (15 cm absolute) in terms of localization precision.
We again observe that as the percentage of missing data increases, the lines corresponding to different values of are crossing among them. This behaviour is very similar to that found in the synthetic experiments above (refer to Fig. 6 and Fig. 7) for a high number of outliers, what suggests that this is the case in the real experiment, also serving as a validation for our simulation conclusions.
It is also noteworthy that, in order to get the best result, the maximum number of outliers should be decreased when the number of missing measurements increases.
XD Comparison with Other Methods in a Localization Task
In this section we made use of the code and real data of the gtde MATLAB toolbox [43, 8]. In this toolbox, real recordings were performed under noisy conditions in a 4x4x4 (m) room, in which an array of 4 microphones forming a tetrahedron of 20 cm side has been placed. Sound waves coming from a loudspeaker placed at 189 different locations 1.7 m away from the array were recorded at 48 KHz.
Table II shows the results of the different algorithms on the localization task. The first 5 columns refer to the results of the Robust DeN algorithm for 5 different values of the parameter . The rest of the columns show the results of a selection of the algorithms implemented and evaluated in [8] (the nomenclature has been kept and the results are essentially the same).
From [8], we selected three multilateration methods (generically referred to as xmult) which are implementations of the algorithm described in [44]. These methods require to be initialized with the distance to the source^{1}^{1}1nmult, tmult and fmult, use values of , and respectively, following [8].. They were selected because they also make use of the redundancy between all the correlations, using the same input as the Robust DeN algorithm. We also compare our proposal with the branch & bound (b&b) method proposed by the authors of [8], as this is the one with the best performance in that work. The interested reader may refer to [8] for further details about the used methods, and the results of some other methods as well.
The Robust DeN algorithm is the only one of the evaluated algorithms that does not require information about the array geometry to perform denoising of the TDOA estimations. The localization algorithm (i.e converting TDOA to angles), implemented in [43], was used after all the methods in order to make the comparison as fair as possible.
The first row in table II shows the percentage of localization measurements with an angular error lower than 30º, the second row shows their mean angular error, and the third row their standard deviation.
Robust DeN  nmult  tmult  fmult  b&b  

% Measurements with angular error  23.78%  24.12%  25.20%  25.71%  23.52%  13.99%  17.38%  17.91%  27.60% 
Mean angular error  17.21  16.44  16.60  16.78  17.66  17.08  16.05  15.25  16.89 
Standard deviation of angular error  7.52  7.44  7.52  7.52  7.38  7.31  7.80  7.78  7.50 
To complement the data of Table II, we have also evaluated the execution timing details of the evaluated algorithms, with the results shown in Table III.
Table II shows that the Robust DeN algorithm performs better than the xmult algorithms. Note that this happens even when the input data is the same, and neither the geometry of the array, nor the distance to the source are used for denonising in our proposal. In what respect to computational demands, our proposal is over 50 times faster.
Comparing with the b&b algorithm, our proposal is close to its performance, which is also an important result provided that (again) we do not use the array geometry for denoising, and our execution time is over 110 times faster.
Robust DeN  xmult  b&b  

Mean time (s)  0.024  1.255  2.649 
Std (s)  0.0011  0.086  0.339 
Xi Conclusions
This paper has studied the properties of TDOA matrices, showing that they can be effectively used for solving TDOA denoising problems. In particular, the paper has investigated challenging scenarios where the TDOA matrix is contaminated with Gaussian noise, outliers and where a percentage of the measurements are missing. The paper shows that denoising in the presence of Gaussian noise and missing data can be solved in closedform. This result is important, as it is the basis of an iterative algorithm that can also cope with outliers. The paper has tested the proposed algorithms in the context of acoustic localization using microphone arrays. The experimental results, both on real and synthetic data have shown that our algorithms successfully perform denoising (up to 30% of improvement in localization accuracy) with a high rate of missing data (up to 50%) and outliers, without knowing the sensor positions. This is important as it opens its application to tasks where the sensors geometry is unknown. Interestingly, in real datasets our robust denoising algorithm is systematically better than the GaussMarkov estimator even when there is no missing data. This is also an important result as it proves that the assumption of Gaussian noise does not hold in real cases, while our robust model is capable of automatically discard erroneous measurements. The proposed robust denoising method has also been compared with other methods in literature on a localization task. Our results are very similar to the state of the art, even though we do not require knowing the array geometry in the denoising stage. Furthermore, our proposal is significantly less computationally demanding.
As for future work, we plan to further test our denoising algorithms in applications where the position of the sensors is unknown in advance, such in selflocalization and beamforming.
Acknowledgements
This work has been supported by the Spanish Ministry of Economy and Competitiveness under project SPACESUAH (TIN201347630C21R), and by the University of Alcalá under projects DETECTOR and ARMIS. Afsaneh Asaei is supported by funding from SNSF project PHASER, grant agreement number 200021153507. We specially thank the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions.
References
 [1] X. Sheng and Y.H. Hu, “Maximum likelihood multiplesource localization using acoustic energy measurements with wireless sensor networks,” IEEE Transactions on Signal Processing, vol. 53, no. 1, pp. 44 – 53, jan. 2005.
 [2] Y. Kuang and K. Astrom, “Stratified sensor network selfcalibration from TDOA measurements,” in Signal Processing Conference (EUSIPCO), 2013 Proceedings of the 21st European. IEEE, 2013, pp. 1–5.
 [3] X. Anguera, C. Wooters, and J. Hernando, “Acoustic beamforming for speaker diarization of meetings,” Audio, Speech, and Language Processing, IEEE Transactions on, vol. 15, no. 7, pp. 2011–2022, 2007.
 [4] M. Brandstein and D. Ward, Eds., Microphone Arrays : Signal Processing Techniques and Applications. Springer, 2001.
 [5] G. C. Carter et al., “Special issue on time delay estimation,” in IEEE Trans. Acoust., Speech, Sig. Process., G. C. Carter, Ed., Jun 1981, vol. ASSP29, no. 3, pp. 461 – 624.
 [6] J. Chen, J. Benesty, and Y. Huang, “Time delay estimation in room acoustic environments: an overview,” EURASIP Journal on applied signal processing, vol. 2006, pp. 170–170, 2006.
 [7] K. Ho, “Bias reduction for an explicit solution of source localization using TDOA,” Signal Processing, IEEE Transactions on, vol. 60, no. 5, pp. 2101–2114, 2012.
 [8] X. AlamedaPineda and R. Horaud, “A geometric approach to sound source localization from timedelay estimates,” Audio, Speech, and Language Processing, IEEE/ACM Transactions on, vol. 22, no. 6, pp. 1082–1095, 2014.
 [9] M. Compagnoni, R. Notari, F. Antonacci, and A. Sarti, “A comprehensive analysis of the geometry of TDOA maps in localization problems,” Inverse Problems, vol. 30, no. 3, p. 035004, 2014.
 [10] A. Nouvellet, M. Charbit, F. Roueff, and A. Le Pichon, “Slowness estimation from noisy time delays observed on nonplanar arrays,” Geophysical Journal International, vol. 198, no. 2, pp. 1199–1207, 2014.
 [11] B. Huang, L. Xie, and Z. Yang, “TDOAbased source localization with distancedependent noises,” Wireless Communications, IEEE Transactions on, vol. 14, no. 1, pp. 468–480, 2015.
 [12] J. DiBiase, H. Silverman, and M. Brandstein, “Robust localization in reverberant rooms,” in Microphone Arrays, ser. Digital Signal Processing, M. Brandstein and D. Ward, Eds. Springer Berlin Heidelberg, 2001, pp. 157–180. [Online]. Available: http://dx.doi.org/10.1007/9783662046197_8
 [13] B. Yang and J. Scheuing, “A theoretical analysis of 2D sensor arrays for TDOA based localization,” in Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on, vol. 4. IEEE, 2006, pp. IV–IV.
 [14] W. Hahn and S. Tretter, “Optimum processing for delayvector estimation in passive signal arrays,” Information Theory, IEEE Transactions on, vol. 19, no. 5, pp. 608–614, Sep 1973.
 [15] N. Zhu, “Locating and extracting acoustic and neural signals,” Wayne State University, Dissertation 422, 2011. [Online]. Available: http://digitalcommons.wayne.edu/oa_dissertations/422
 [16] Y. Chan and K. Ho, “A simple and efficient estimator for hyperbolic location,” Signal Processing, IEEE Transactions on, vol. 42, no. 8, pp. 1905–1915, 1994.
 [17] M. Pollefeys and D. Nister, “Direct computation of sound and microphone locations from timedifferenceofarrival data.” in ICASSP, 2008, pp. 2445–2448.
 [18] J. O. Smith and J. S. Abel, “Closedform leastsquares source location estimation from rangedifference measurements,” Acoustics, Speech and Signal Processing, IEEE Transactions on, vol. 35, no. 12, pp. 1661–1669, 1987.
 [19] M. D. Gillette and H. F. Silverman, “A linear closedform algorithm for source localization from timedifferences of arrival,” Signal Processing Letters, IEEE, vol. 15, pp. 1–4, 2008.
 [20] Y. Weng, W. Xiao, and L. Xie, “Total least squares method for robust source localization in sensor networks using TDOA measurements,” International Journal of Distributed Sensor Networks, vol. 2011, 2011.
 [21] L. Lin, H.C. So, F. K. Chan, Y. Chan, and K. Ho, “A new constrained weighted least squares algorithm for TDOAbased localization,” Signal Processing, vol. 93, no. 11, pp. 2872–2878, 2013.
 [22] H. JamaliRad and G. Leus, “Sparsityaware multisource TDOA localization,” Signal Processing, IEEE Transactions on, vol. 61, no. 19, pp. 4874–4887, 2013.
 [23] H. C. So, Y. T. Chan, and F. Chan, “Closedform formulae for timedifferenceofarrival estimation,” Signal Processing, IEEE Transactions on, vol. 56, no. 6, pp. 2614–2620, June 2008.
 [24] A. Renaux, P. Forster, E. Boyer, and P. Larzabal, “Unconditional maximum likelihood performance at finite number of samples and high signaltonoise ratio,” Signal Processing, IEEE Transactions on, vol. 55, no. 5, pp. 2358–2364, 2007.
 [25] R. Schmidt, “Least squares range difference location,” Aerospace and Electronic Systems, IEEE Transactions on, vol. 32, no. 1, pp. 234–242, Jan 1996.
 [26] C. Knapp and G. Carter, “The generalized correlation method for estimation of time delay,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 24, no. 4, pp. 320 – 327, aug 1976.
 [27] C. Zhang, D. Florencio, and Z. Zhang, “Why does PHAT work well in low noise, reverberative environments?” in Proceedings of ICASSP 2008, 31 2008april 4 2008, pp. 2565 –2568.
 [28] J. Velasco, C. J. MartínArguedas, J. MaciasGuarasa, D. Pizarro, and M. Mazo, “Proposal and validation of an analytical generative model of SRPPHAT power maps in reverberant scenarios,” Signal Processing, vol. 119, pp. 209 – 228, 2016. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0165168415002650
 [29] J. Benesty, “Adaptive eigenvalue decomposition algorithm for passive acoustic source localization,” The Journal of the Acoustical Society of America, vol. 107, no. 1, pp. 384–391, 2000.
 [30] J. Scheuing and B. Yang, “Disambiguation of TDOA estimation for multiple sources in reverberant environments,” Audio, Speech, and Language Processing, IEEE Transactions on, vol. 16, no. 8, pp. 1479–1489, 2008.
 [31] P. Annibale, J. Filos, P. Naylor, R. Rabenstein et al., “TDOAbased speed of sound estimation for air temperature and room geometry inference,” Audio, Speech, and Language Processing, IEEE Transactions on, vol. 21, no. 2, pp. 234–246, 2013.
 [32] M. Compagnoni, A. Canclini, P. Bestagini, F. Antonacci, A. Sarti, and S. Tubaro, “TDOA denoising for acoustic source localization,” CoRR, vol. abs/1509.02380, 2015. [Online]. Available: http://arxiv.org/abs/1509.02380
 [33] G. Ye, D. Liu, I.H. Jhuo, and S.F. Chang, “Robust late fusion with rank minimization,” in Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on. IEEE, 2012, pp. 3021–3028.
 [34] C. Eckart and G. Young, “The approximation of one matrix by another of lower rank,” Psychometrika, vol. 1, no. 3, pp. 211–218, 1936.
 [35] L. Mirsky, “Symmetric gauge functions and unitarily invariant norms,” The quarterly journal of mathematics, vol. 11, no. 1, pp. 50–59, 1960.
 [36] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis,” Journal of the ACM (JACM), vol. 58, no. 3, p. 11, 2011.
 [37] T. Zhou and D. Tao, “GoDec: Randomized lowrank & sparse matrix decomposition in noisy case,” in in International Conference on Machine Learning, 2011.
 [38] A. Goldberg, B. Recht, J. Xu, R. Nowak, and X. Zhu, “Transduction with matrix completion: Three birds with one stone,” in Advances in Neural Information Processing Systems 23, J. D. Lafferty, C. K. I. Williams, J. ShaweTaylor, R. S. Zemel, and A. Culotta, Eds. Curran Associates, Inc., 2010, pp. 757–765. [Online]. Available: http://papers.nips.cc/paper/3932transductionwithmatrixcompletionthreebirdswithonestone.pdf
 [39] E. J. Candès and T. Tao, “The power of convex relaxation: Nearoptimal matrix completion,” IEEE Trans. Inf. Theor., vol. 56, no. 5, pp. 2053–2080, May 2010. [Online]. Available: http://dx.doi.org/10.1109/TIT.2010.2044061
 [40] G. Lathoud, J.M. Odobez, and D. GaticaPerez, “AV16.3: An audiovisual corpus for speaker localization and tracking,” in Procceedings of the MLMI, ser. Lecture Notes in Computer Science, S. Bengio and H. Bourlard, Eds., vol. 3361. SpringerVerlag, 2004, pp. 182–195.
 [41] D. C. Moore, “The IDIAP smart meeting room,” IDIAP Research Institute, Switzerland, Tech. Rep., November 2004.
 [42] G. Lathoud, “AV16.3 dataset,” http://www.idiap.ch/dataset/av163/, [Last accessed in december 2015].
 [43] X. AlamedaPineda, “The gtde matlab toolbox.” [Online]. Available: https://team.inria.fr/perception/thegtdematlabtoolbox/
 [44] M. S. Brandstein and H. F. Silverman, “A practical methodology for speech source localization with microphone arrays,” Computer Speech & Language, vol. 11, no. 2, pp. 91–126, 1997.