Node Localization in Unsynchronized Time of Arrival Sensor Networks
Node Localization in Unsynchronized Time of Arrival Sensor Networks ˚ om Simon Burgess, Yubin Kuang and Kalle Astr¨ Centre for Mathematical Sciences, Lund University, Lund, Sweden {simonb, yubin, kalle}@maths.lth.se
Abstract We present a method for solving the previously unstudied problem of localizing a set of receivers and directions from transmitters placed far from the receivers, measuring unsynchronized time of arrival data. The same problem is present in node localization of microphone and antenna arrays. The solution algorithm using 5 receivers and 9 transmitters is extended to the overdetermined case in a straightforward manner. Degenerate cases are shown to be when i) the measurement matrix has rank 4 or less or ii) the directions from the transmitters to the receivers lie on an intersection between the unit sphere and another quadric surface. In simulated experiments we explore how sensitive the solution is with respect to different degrees of far field approximations of the transmitters and with respect to noise in the data. Using real data we get a reconstruction of the receivers with a relative error of 14%.
1. Introduction In this paper we study the problem of node localization using unsynchronized time of arrival data. The problem arises naturally in microphone arrays for audio sensing. Using multiple unsynchronized microphones or recorded sound files, and matching sounds emanating from unknown locations and unknown time, is it possible to both calculate the microphone positions as well as the timings and directions of the sound sources? An example application could be when several cell phones record distinct environmental sounds, and the microphone positions and sound directions is to be recovered without the need of synchronizing the phones. The identical mathematical problem occurs as a part of dealing with uncalibrated arrays of antennas. By measuring phase response for each antenna from a transmitter and repeating the measurement with different transmitters, how can the constellation and phase
offset of the antennas as well as offsets and directions to the transmitters be obtained? The phase difference for each signal can be measured up to an unknown offset for each antenna, and an offset for each transmitter. If the array is smaller than half a wave length in width, which is often the case for radio applications, the phase differences can readily be translated into distance measurements with the phase offsets now as distance offsets, and using the known speed of the signals we have unsynchronized time of arrival measurements. Although time of arrival (TOA) and time difference of arrival (TDOA) problems have been studied extensively in the literature in the form of localization of e.g. a sound source using a calibrated array, array [3, 4, 5, 6], the problem of calibration of a sensor array using only measurement, i.e. the node localization problem, has received much less attention. One solution is to manually measure the inter-distance between pairs of receivers and use multi-dimensional scaling to compute node locations [1]. Another option is to use GPS [9] or to use additional transmitters close to each receiver [11, 12]. Sensor network calibration is treated in [2]. Initialization of TDOA network node localization is studied in [10], where solutions were given to two non-minimal cases of ten transmitters and five receivers, whereas the minimal solution for far field approximation in [8] are six transmitters and four receivers. In all of these cases it has been assumed that there is synchronization at either the receiver or transmitter end. In this paper we study the problem of sensor network initialization of both unsynchronized transmitters and unsynchronized receivers in a far field setting, giving unsynchronized time of arrival (UTOA) measurements. This problem have to the authors’ best knowledge not been studied before. We show that for at least 5 receivers and at least 9 transmitter the problem is solvable. Numerically stable algorithms for parameter estimation have been developed both for initialization and for nonlinear least squares minimization. The methods is tested on both simulated and real data.
1.1. Far field UTOA-measurements
2.1. Gauge freedom
In the following treatment, we make no difference between real and virtual transmitters. Assume that the transmitters are stationary at position bj ∈ R3 , j = 1, . . . , k and that the receivers are at positions ri ∈ R3 , i = 1, . . . , m. By measuring how long time the signals take to reach the receiver and knowing the speed of the signals, distances dij = |ri − bj | can be measured. These are TOA measurements. When neither receivers or transmitters are synchronised, for instance external sound sources recorded on different computers, the measurements will be of the form dij = |ri − bj | + fi + gˆj where fi , gˆj are unknown offsets for receivers and transmitters respectively. We denote measurements of this kind Unsynchronized Time Of Arrival (UTOA) measurements. Furthermore, if the transmitters are so far from the receivers that a transmitter can be considered to have a common direction to the receivers, the measurements can be approximated by dij = |ri − bj | + fi + gˆj ≈ |r1 − bj | + (ri − r1 )T nj + fi + gˆj = riT nj + g¯j + fi + gˆj where g¯j = |r1 − bj | − r1T nj and nj is the direction of unit length from transmitter j to the receivers. By setting gj = g¯j + gˆj we get the far field approximation
The unknown parameters (r, n, f, g) have certain gauge freedom. There is a translation t, rotation matrix R and offset change K which can be applied to the solution according to
dij ≈ riT nj + fi + gj . When the approximation is good, we will call these Far Field UTOA (FFUTOA) measurements.
2. The FFUTOA localization problem Problem 2.1 Given mk FFUTOA measurements dij ∈ R, i = 1, . . . , m, j = 1, . . . , k, taken from m receivers and k transmitters, estimate receiver positions r = (r1 , . . . , rm ), ri ∈ R3 , directions n = (n1 , . . . , nk ), nj ∈ R3 from transmitter j to receivers such that d˜ij = riT nj + fi + gj (1) ||nj ||2 = 1 where fi ∈ R, gj ∈ R are constant distance offsets for receiver i and transmitter j respectively. Note that the problem is symmetric in receivers and transmitters, i.e. if each receiver instead could be viewed as having a common direction to all transmitters, the same problem can be solved for transmitter positions and receiver directions. We denote f = [f1 , . . . , fm ]T , g = [g1 , . . . , gk ] and n = [n1 , . . . , nk ].
r¯i = ri + t, g¯j = gj − tT nj r¯i = Rri , n ¯ j = Rnj f¯i = fi + K, g¯j = gj − K without changing the measurements. We lock the gauge freedom by setting r1 = 0, r2 r3 r4 = U, f1 = 0 where U is upper triangular with Uii > 0.
2.2. A method for solving problem 2.1 Without loss of generality we assume that the solution is normalized for gauge freedom as in section 2.1. Using the FFUTOA measurements d˜ij , collected in ˜ = [d˜ij ]m×k we immediately obtain the the matrix D unknowns gj since d˜1j = r1T nj + f1 + gj = gj , since r1 = 0 and f1 = 0. Using these values we can subtract gj from all measurements and obtain a new matrix that fulfill n (2) D = rT f ~ 1 where ~1 is a vector of ones. D is a product of two matrices of rank ≤ 4, given that the number of receivers m and transmitters k are ≥ 3. This gives the first result. ˜ has Lemma 2.1 The measurement matrix D ˜ ≤ 5 and the matrix D has Rank(D) ≤ 4. Rank(D) To solve for the unknowns we use singular value decomposition (svd) D = U ΣV˜ T = U V where U and V are m × 4 and 4 × k matrices respectively. This gives the solution up to an unknown 4 × 4 invertible matrix H, i.e. D = U H −1 HV . where U H −1 and HV are the two matrices in (2). We write the 4 × 4 matrix H as H1 H= h4 where H1 is a 3 × 4 matrix and h4 a 1 × 4 vector. Noting that the last row of HV should be ~1 we obtain the linear constraints h4 Vj = 1. The quadratic constraints nTj nj = 1 gives constraints on the upper 3 × 4 part H1 of H according to VjT H1T H1 Vj = 1, which can be phrased as a linear
problem in C = H1T H1 , i.e. find the symmetric rank 3 matrix C such that VjT CVj = 1.
(3)
This is a linear problem in the ten unknowns in C, called c. Ignoring the rank constraint, the linear problem can be formulated as Ac = 1. (4) C = hT4 h4 is one solution, and there exists at least one other solution of rank 3, C = H1T H1 . Thus A must have Rank(A) ≤ 9. With 9 or more directions nj , A in general has rank 9. Section 2.4 characterizes degenerate cases. The solution of (3) thus gives a one parameter family of solutions C (t) = C0 + tC1 . Enforcing the constraint det(C(t)) = 0
(5)
gives a fourth order polynomial in t. The polynomial can be calculated in a efficient and stable manner by interpolation using [7]. The following lemma explains what we can expect of the solutions to (5). Lemma 2.2 The fourth order polynomial equation (5) has a real triple root t1 encoding C(t1 ) = hT4 h4 and a single real root t2 encoding C(t2 ) = H1T H1 Proof: As the solution space {C(t)} to (3) contains C(t1 ) = hT4 h4 , we know that there exists a rank 1 matrix C(t1 ), thus having three singular values equal to zero. As C(t1 ) is hermitian, and hermitian matrices singular values are absolute values of the eigenvalues, C(t1 ) has 3 eigenvalues equal to zero and thus t1 is a triple root to (5). t1 must be real as C(t1 ) = hT4 h4 , C0 and C1 are real. There is one more real root t2 to (5) as it is a fourth order polynomial with real coefficients. As we know that C = H1T H1 exists in the solution space and it is not of full rank, C(t2 ) = H1T H1 . The matrix H1 can then be calculated from a svd ˆ SˆU ˆ T which gives H1 as the three first of C(t2 )p= U ˆT rows of SˆU .T Using these calculated values we obT tain H = H1 hT4 and then change coordinates to obtain T n r f = U H −1 , ~ = HV. 1 Accordingly, we have the following algorithm: Algorithm 2.1 ˜ of size 5 × 9. Input: FFUTOA measurement matrix D Output: Receiver positions r, transmitter directions n, receiver offsets f and transmitter offsets g ˜ 1j and Dij := D ˜ ij − gj 1. Set gj := D
2. Remove the first row of D ˜ S V˜ T and set U to first 4 3. Calculate the svd D = U ˜ columns of U and V to first four columns of S V˜ T 4. Solve for h4 , h4 V = ~1 5. Find the solution space {C(t)} for the ten unknowns in the symmetric matrix C using the nine linear constraints VjT CVj = 1 6. Solve for t the fourth order polynomial equation det(C(t)) = 0 and set t2 to the single root ˆ SˆU ˆ T and set H1 to 7. Calculate the svd p C(t2 ) = U ˆT SˆU first three rows of T T T 8. Set H := H1 h4 , rT f := U H −1 and set n to three first rows of HV Note that the algorithm is not a true minimal case as when the data are not true FFUTOA measurements, step 4 will not have a solution as there are four unknowns in h4 but 9 linear constraints. Section 2.3 addresses this. The algorithm is however close to minimal, as when using m = 5 receivers and k = 9 transmitters, the number of measurements mk = 45 are close to the number of unknowns, 4m + 3k − 7 = 40 when taking gauge freedom into account.
2.3. Extension to overdetermined cases When data are corrupted by noise, the far field approximation is not good or when having more than 5 receivers and 9 transmitters, the extension of Algorithm 2.1 can be made as follows: (i) The best rank 4 approximation in Frobenius norm can still be found in step 3 by svd. (ii) In step 4, the linear least squares solution of h4 is obtained. (iii) As VjT hT4 h4 Vj 6= 1, this is no longer a true solution to (3), and in general A in (4) will have full rank. A is thus approximated by its closest rank 9 matrix by svd. (iv) In step 6, the roots are no longer necessarily one triple root and one single root, but a perturbation of this constellation. Thus, t2 is selected amongst the four roots of (5) as the real root such that C(t2 ) differs the most from hT4 h4 in Frobenius norm. This extended algorithm has no guarantee of optimality, but can serve as an initialization for nonlinear optimization schemes as seen in section 3.1. From here on the extended algorithm is the one that will be used. Note that when data are true FFUTOA measurements, the extended algorithm is equivalent to Algorithm 2.1.
2.4. Degenerate cases It is enlightening to know the failure modes of the algorithm. This is captured by the following theorem.
Theorem 2.1 Algorithm 2.1 fails if the measurement ˜ has Rank(D) ˜ ≤ 4. It also fails if the dimatrix D rections nj from the transmitters lie on the intersection of a quadric surface and the unit sphere.
120
100
where Cs is the symmetric matrix giving the unit sphere, and C2 6= λCs , λ ∈ R is another symmetric matrix. This is equivalent to VjT H T Cs HVj = 0, VjT H T C2 HVj = 0,
Frequency
80
˜ has rank ≤ 4, Proof: If the measurement matrix D the matrix D will have rank ≤ 3 and step 3 in the algorithm will not be meaningful. That nj also lies on another quadric surface than the unit sphere means that nj fulfills T T T nj 1 C2 nTj 1 = 0, nTj 1 Cs nTj 1 = 0
60
40
20
0 −15
−14.5
−14
−13.5 −13 −12.5 −12 log10 of error in position
−11.5
−11
−10.5
Figure 1. Histogram over errors of receiver position for Algorithm 2.1
(6)
as nTj 1 = VjT H. This implies that VjT CVj = 0 has two linearly independent solutions and thus (4) has a null space of at least dimension two. Then step 5 in the algorithm will fail as it needs the null space to be of dimension one. ˜ ≤ 4 is when Note that a special case of Rank(D) either the receivers ri or transmitters directions nj lie in a plane. If both lie in the same plane, a similar algorithm for 2D can be constructed following the same ideas as for the 3D case.
To evaluate the assumption that each transmitter has a common direction to receivers, data was generated using 3D positions for both transmitters and receivers at different distances in-between receivers and transmitters. In each experiment, transmitters positions are placed randomly over the sphere of radius d. A graph showing the mean error as a function of the radius is displayed in Figure 2. 0
10
r5 − t9 r5 − t10 r6 − t9 r6 − t10 r7 − t11
−1
10
3. Experimental validation Error in Position
For all experiments, errors are defined as the l2 -norm of difference between reconstructed receiver positions r and ground truth.
−2
10
−3
10
−4
10
−5
10
−6
10
−7
3.1. Simulations
10
−8
10
2
4
10
10
6
10
Approximate distance from tranmitters to receivers
For all simulations the following setup was used. For the minimal case, four receivers are placed in the corners of a tetrahedron inscribed in a cube of unit volume centered around the origin, and one receiver is placed in the origin. For the overdetermined cases, additional receivers are drawn from an i.i.d. uniform distribution over the same cube. Offsets fi and gj are drawn from i.i.d. uniform distributions over [0, 10]. The numerical performance of the algorithm was evaluated by generating problems where the measure˜ are FFUTOA. Nine transmitter direcment matrix D tions were simulated over the unit sphere. Measurements were then constructed according to (1). The errors over 1000 such experiments can be seen in Figure 1
Figure 2. Effect of quality of far filed approximation on number of receivers (r) and transmitters (t)
To evaluate the sensitivity of noise, white gaussian noise was added to UTOA measurements, simulated with transmitters 107 from the origin. 10 transmitters and 6 receivers were used. A graph showing the mean error of 100 experiments plotted against the standard deviation of the noise is shown in Figure 3. The reconstruction was also optimized under the function 2 P˜ min r,n,f,g dij − rT nj − fi − gj using a Gaussi
Newton like approach. The result is also shown in Figure 3.
−2
”Wearable Visual Systems” and the Swedish Research Council project ”Polynomial Equations and Geometry in Computer Vision”.
−3
References
−1
10
initialization opitmized
Error in Position
10
10
−4
10
−5
10
−5
10
−4
10
−3
10
−2
10
Standerd deviation
Figure 3. Effect of additive white gaussian noise.
3.2. Real Data An experimental setup using 8 microphones as receivers and random distinct manually made sounds as transmitters was used. The 19 sound sources were ≈30 m away from the receivers. Microphones were set in the corners of a cuboid of roughly 100 × 105 × 60 cm3 . The beginning of each sound were matched by a heuristic algorithm to create TDOA measurements. As microphones were recorded through a sound card, an artificial offset had to be added to each microphone to create UTOA measurements. Using the same real measurements, but adding new receiver offsets for each experiment, the algorithm was run on 100 experiments. A mean error of 13 cm was obtained. 85% of the errors are in the floorto-roof direction. This can be explained by the sounds all being made close to ground level and thus the transmitter directions will be close to being in a plane, giving poor resolution in floor-to-roof direction.
4. Conclusions We have presented a novel method for solving the, to the authors’ best knowledge, previously unstudied problem of initializing a sensor network using fully unsynchronized time of arrival data in a far field setting. Degenerate constellations are also characterized. By using real and simulated data we show the applicability of the method. Experimental validation also shows that the method is fairly robust to both noise and settings where the approximation that the transmitters are far from the receivers is not good.
5
Acknowledgements
The research leading to these results has received funding from the swedish strategic research project ELLIIT, the Swedish Strategic Foundation funded project
[1] S. T. Birchfield and A. Subramanya. Microphone array position calibration by basis-point classical multidimensional scaling. IEEE transactions on Speech and Audio Processing, 13(5), 2005. [2] R. Biswas and S. Thrun. A passive approach to sensor network localization. In IROS, 2004. [3] M. Brandstein, J. Adcock, and H. Silverman. A closedform location estimator for use with room environment microphone arrays. Speech and Audio Processing, IEEE Transactions on, 5(1):45 –50, jan 1997. [4] A. Cirillo, R. Parisi, and A. Uncini. Sound mapping in reverberant rooms by a robust direct method. In Acoustics, Speech and Signal Processing, 2008. ICASSP 2008. IEEE International Conference on, pages 285 – 288, 31 2008-april 4 2008. [5] M. Cobos, A. Marti, and J. Lopez. A modified srp-phat functional for robust real-time sound source localization with scalable spatial sampling. Signal Processing Letters, IEEE, 18(1):71 –74, jan. 2011. [6] H. Do, H. Silverman, and Y. Yu. A real-time srpphat source location implementation using stochastic region contraction(src) on a large-aperture microphone array. In Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference on, volume 1, pages I–121 –I–124, april 2007. [7] M. Hromcik and M. Sebek. Numerical and symbolic computation of polynomial matrix determinant. In Decision and Control, 1999. Proceedings of the 38th IEEE Conference on, volume 2, pages 1887 –1888, 1999. ˚ om. Under[8] Y. Kuang, E. Ask, S. Burgess, and K. Astr¨ standing toa and tdoa network calibration using far field approximation as initial estimate. In ICPRAM, 2012. [9] D. Niculescu and B. Nath. Ad hoc positioning system (aps). In GLOBECOM-01, 2001. [10] M. Pollefeys and D. Nister. Direct computation of sound and microphone locations from time-differenceof-arrival data. In ICASSP, 2008. [11] V. C. Raykar, I. V. Kozintsev, and R. Lienhart. Position calibration of microphones and loudspeakers in distributed computing platforms. IEEE transactions on Speech and Audio Processing, 13(1), 2005. [12] J. Sallai, G. Balogh, M. Maroti, and A. Ledeczi. Acoustic ranging in resource-constrained sensor networks. In eCOTS-04, 2004.
Abstract We present a method for solving the previously unstudied problem of localizing a set of receivers and directions from transmitters placed far from the receivers, measuring unsynchronized time of arrival data. The same problem is present in node localization of microphone and antenna arrays. The solution algorithm using 5 receivers and 9 transmitters is extended to the overdetermined case in a straightforward manner. Degenerate cases are shown to be when i) the measurement matrix has rank 4 or less or ii) the directions from the transmitters to the receivers lie on an intersection between the unit sphere and another quadric surface. In simulated experiments we explore how sensitive the solution is with respect to different degrees of far field approximations of the transmitters and with respect to noise in the data. Using real data we get a reconstruction of the receivers with a relative error of 14%.
1. Introduction In this paper we study the problem of node localization using unsynchronized time of arrival data. The problem arises naturally in microphone arrays for audio sensing. Using multiple unsynchronized microphones or recorded sound files, and matching sounds emanating from unknown locations and unknown time, is it possible to both calculate the microphone positions as well as the timings and directions of the sound sources? An example application could be when several cell phones record distinct environmental sounds, and the microphone positions and sound directions is to be recovered without the need of synchronizing the phones. The identical mathematical problem occurs as a part of dealing with uncalibrated arrays of antennas. By measuring phase response for each antenna from a transmitter and repeating the measurement with different transmitters, how can the constellation and phase
offset of the antennas as well as offsets and directions to the transmitters be obtained? The phase difference for each signal can be measured up to an unknown offset for each antenna, and an offset for each transmitter. If the array is smaller than half a wave length in width, which is often the case for radio applications, the phase differences can readily be translated into distance measurements with the phase offsets now as distance offsets, and using the known speed of the signals we have unsynchronized time of arrival measurements. Although time of arrival (TOA) and time difference of arrival (TDOA) problems have been studied extensively in the literature in the form of localization of e.g. a sound source using a calibrated array, array [3, 4, 5, 6], the problem of calibration of a sensor array using only measurement, i.e. the node localization problem, has received much less attention. One solution is to manually measure the inter-distance between pairs of receivers and use multi-dimensional scaling to compute node locations [1]. Another option is to use GPS [9] or to use additional transmitters close to each receiver [11, 12]. Sensor network calibration is treated in [2]. Initialization of TDOA network node localization is studied in [10], where solutions were given to two non-minimal cases of ten transmitters and five receivers, whereas the minimal solution for far field approximation in [8] are six transmitters and four receivers. In all of these cases it has been assumed that there is synchronization at either the receiver or transmitter end. In this paper we study the problem of sensor network initialization of both unsynchronized transmitters and unsynchronized receivers in a far field setting, giving unsynchronized time of arrival (UTOA) measurements. This problem have to the authors’ best knowledge not been studied before. We show that for at least 5 receivers and at least 9 transmitter the problem is solvable. Numerically stable algorithms for parameter estimation have been developed both for initialization and for nonlinear least squares minimization. The methods is tested on both simulated and real data.
1.1. Far field UTOA-measurements
2.1. Gauge freedom
In the following treatment, we make no difference between real and virtual transmitters. Assume that the transmitters are stationary at position bj ∈ R3 , j = 1, . . . , k and that the receivers are at positions ri ∈ R3 , i = 1, . . . , m. By measuring how long time the signals take to reach the receiver and knowing the speed of the signals, distances dij = |ri − bj | can be measured. These are TOA measurements. When neither receivers or transmitters are synchronised, for instance external sound sources recorded on different computers, the measurements will be of the form dij = |ri − bj | + fi + gˆj where fi , gˆj are unknown offsets for receivers and transmitters respectively. We denote measurements of this kind Unsynchronized Time Of Arrival (UTOA) measurements. Furthermore, if the transmitters are so far from the receivers that a transmitter can be considered to have a common direction to the receivers, the measurements can be approximated by dij = |ri − bj | + fi + gˆj ≈ |r1 − bj | + (ri − r1 )T nj + fi + gˆj = riT nj + g¯j + fi + gˆj where g¯j = |r1 − bj | − r1T nj and nj is the direction of unit length from transmitter j to the receivers. By setting gj = g¯j + gˆj we get the far field approximation
The unknown parameters (r, n, f, g) have certain gauge freedom. There is a translation t, rotation matrix R and offset change K which can be applied to the solution according to
dij ≈ riT nj + fi + gj . When the approximation is good, we will call these Far Field UTOA (FFUTOA) measurements.
2. The FFUTOA localization problem Problem 2.1 Given mk FFUTOA measurements dij ∈ R, i = 1, . . . , m, j = 1, . . . , k, taken from m receivers and k transmitters, estimate receiver positions r = (r1 , . . . , rm ), ri ∈ R3 , directions n = (n1 , . . . , nk ), nj ∈ R3 from transmitter j to receivers such that d˜ij = riT nj + fi + gj (1) ||nj ||2 = 1 where fi ∈ R, gj ∈ R are constant distance offsets for receiver i and transmitter j respectively. Note that the problem is symmetric in receivers and transmitters, i.e. if each receiver instead could be viewed as having a common direction to all transmitters, the same problem can be solved for transmitter positions and receiver directions. We denote f = [f1 , . . . , fm ]T , g = [g1 , . . . , gk ] and n = [n1 , . . . , nk ].
r¯i = ri + t, g¯j = gj − tT nj r¯i = Rri , n ¯ j = Rnj f¯i = fi + K, g¯j = gj − K without changing the measurements. We lock the gauge freedom by setting r1 = 0, r2 r3 r4 = U, f1 = 0 where U is upper triangular with Uii > 0.
2.2. A method for solving problem 2.1 Without loss of generality we assume that the solution is normalized for gauge freedom as in section 2.1. Using the FFUTOA measurements d˜ij , collected in ˜ = [d˜ij ]m×k we immediately obtain the the matrix D unknowns gj since d˜1j = r1T nj + f1 + gj = gj , since r1 = 0 and f1 = 0. Using these values we can subtract gj from all measurements and obtain a new matrix that fulfill n (2) D = rT f ~ 1 where ~1 is a vector of ones. D is a product of two matrices of rank ≤ 4, given that the number of receivers m and transmitters k are ≥ 3. This gives the first result. ˜ has Lemma 2.1 The measurement matrix D ˜ ≤ 5 and the matrix D has Rank(D) ≤ 4. Rank(D) To solve for the unknowns we use singular value decomposition (svd) D = U ΣV˜ T = U V where U and V are m × 4 and 4 × k matrices respectively. This gives the solution up to an unknown 4 × 4 invertible matrix H, i.e. D = U H −1 HV . where U H −1 and HV are the two matrices in (2). We write the 4 × 4 matrix H as H1 H= h4 where H1 is a 3 × 4 matrix and h4 a 1 × 4 vector. Noting that the last row of HV should be ~1 we obtain the linear constraints h4 Vj = 1. The quadratic constraints nTj nj = 1 gives constraints on the upper 3 × 4 part H1 of H according to VjT H1T H1 Vj = 1, which can be phrased as a linear
problem in C = H1T H1 , i.e. find the symmetric rank 3 matrix C such that VjT CVj = 1.
(3)
This is a linear problem in the ten unknowns in C, called c. Ignoring the rank constraint, the linear problem can be formulated as Ac = 1. (4) C = hT4 h4 is one solution, and there exists at least one other solution of rank 3, C = H1T H1 . Thus A must have Rank(A) ≤ 9. With 9 or more directions nj , A in general has rank 9. Section 2.4 characterizes degenerate cases. The solution of (3) thus gives a one parameter family of solutions C (t) = C0 + tC1 . Enforcing the constraint det(C(t)) = 0
(5)
gives a fourth order polynomial in t. The polynomial can be calculated in a efficient and stable manner by interpolation using [7]. The following lemma explains what we can expect of the solutions to (5). Lemma 2.2 The fourth order polynomial equation (5) has a real triple root t1 encoding C(t1 ) = hT4 h4 and a single real root t2 encoding C(t2 ) = H1T H1 Proof: As the solution space {C(t)} to (3) contains C(t1 ) = hT4 h4 , we know that there exists a rank 1 matrix C(t1 ), thus having three singular values equal to zero. As C(t1 ) is hermitian, and hermitian matrices singular values are absolute values of the eigenvalues, C(t1 ) has 3 eigenvalues equal to zero and thus t1 is a triple root to (5). t1 must be real as C(t1 ) = hT4 h4 , C0 and C1 are real. There is one more real root t2 to (5) as it is a fourth order polynomial with real coefficients. As we know that C = H1T H1 exists in the solution space and it is not of full rank, C(t2 ) = H1T H1 . The matrix H1 can then be calculated from a svd ˆ SˆU ˆ T which gives H1 as the three first of C(t2 )p= U ˆT rows of SˆU .T Using these calculated values we obT tain H = H1 hT4 and then change coordinates to obtain T n r f = U H −1 , ~ = HV. 1 Accordingly, we have the following algorithm: Algorithm 2.1 ˜ of size 5 × 9. Input: FFUTOA measurement matrix D Output: Receiver positions r, transmitter directions n, receiver offsets f and transmitter offsets g ˜ 1j and Dij := D ˜ ij − gj 1. Set gj := D
2. Remove the first row of D ˜ S V˜ T and set U to first 4 3. Calculate the svd D = U ˜ columns of U and V to first four columns of S V˜ T 4. Solve for h4 , h4 V = ~1 5. Find the solution space {C(t)} for the ten unknowns in the symmetric matrix C using the nine linear constraints VjT CVj = 1 6. Solve for t the fourth order polynomial equation det(C(t)) = 0 and set t2 to the single root ˆ SˆU ˆ T and set H1 to 7. Calculate the svd p C(t2 ) = U ˆT SˆU first three rows of T T T 8. Set H := H1 h4 , rT f := U H −1 and set n to three first rows of HV Note that the algorithm is not a true minimal case as when the data are not true FFUTOA measurements, step 4 will not have a solution as there are four unknowns in h4 but 9 linear constraints. Section 2.3 addresses this. The algorithm is however close to minimal, as when using m = 5 receivers and k = 9 transmitters, the number of measurements mk = 45 are close to the number of unknowns, 4m + 3k − 7 = 40 when taking gauge freedom into account.
2.3. Extension to overdetermined cases When data are corrupted by noise, the far field approximation is not good or when having more than 5 receivers and 9 transmitters, the extension of Algorithm 2.1 can be made as follows: (i) The best rank 4 approximation in Frobenius norm can still be found in step 3 by svd. (ii) In step 4, the linear least squares solution of h4 is obtained. (iii) As VjT hT4 h4 Vj 6= 1, this is no longer a true solution to (3), and in general A in (4) will have full rank. A is thus approximated by its closest rank 9 matrix by svd. (iv) In step 6, the roots are no longer necessarily one triple root and one single root, but a perturbation of this constellation. Thus, t2 is selected amongst the four roots of (5) as the real root such that C(t2 ) differs the most from hT4 h4 in Frobenius norm. This extended algorithm has no guarantee of optimality, but can serve as an initialization for nonlinear optimization schemes as seen in section 3.1. From here on the extended algorithm is the one that will be used. Note that when data are true FFUTOA measurements, the extended algorithm is equivalent to Algorithm 2.1.
2.4. Degenerate cases It is enlightening to know the failure modes of the algorithm. This is captured by the following theorem.
Theorem 2.1 Algorithm 2.1 fails if the measurement ˜ has Rank(D) ˜ ≤ 4. It also fails if the dimatrix D rections nj from the transmitters lie on the intersection of a quadric surface and the unit sphere.
120
100
where Cs is the symmetric matrix giving the unit sphere, and C2 6= λCs , λ ∈ R is another symmetric matrix. This is equivalent to VjT H T Cs HVj = 0, VjT H T C2 HVj = 0,
Frequency
80
˜ has rank ≤ 4, Proof: If the measurement matrix D the matrix D will have rank ≤ 3 and step 3 in the algorithm will not be meaningful. That nj also lies on another quadric surface than the unit sphere means that nj fulfills T T T nj 1 C2 nTj 1 = 0, nTj 1 Cs nTj 1 = 0
60
40
20
0 −15
−14.5
−14
−13.5 −13 −12.5 −12 log10 of error in position
−11.5
−11
−10.5
Figure 1. Histogram over errors of receiver position for Algorithm 2.1
(6)
as nTj 1 = VjT H. This implies that VjT CVj = 0 has two linearly independent solutions and thus (4) has a null space of at least dimension two. Then step 5 in the algorithm will fail as it needs the null space to be of dimension one. ˜ ≤ 4 is when Note that a special case of Rank(D) either the receivers ri or transmitters directions nj lie in a plane. If both lie in the same plane, a similar algorithm for 2D can be constructed following the same ideas as for the 3D case.
To evaluate the assumption that each transmitter has a common direction to receivers, data was generated using 3D positions for both transmitters and receivers at different distances in-between receivers and transmitters. In each experiment, transmitters positions are placed randomly over the sphere of radius d. A graph showing the mean error as a function of the radius is displayed in Figure 2. 0
10
r5 − t9 r5 − t10 r6 − t9 r6 − t10 r7 − t11
−1
10
3. Experimental validation Error in Position
For all experiments, errors are defined as the l2 -norm of difference between reconstructed receiver positions r and ground truth.
−2
10
−3
10
−4
10
−5
10
−6
10
−7
3.1. Simulations
10
−8
10
2
4
10
10
6
10
Approximate distance from tranmitters to receivers
For all simulations the following setup was used. For the minimal case, four receivers are placed in the corners of a tetrahedron inscribed in a cube of unit volume centered around the origin, and one receiver is placed in the origin. For the overdetermined cases, additional receivers are drawn from an i.i.d. uniform distribution over the same cube. Offsets fi and gj are drawn from i.i.d. uniform distributions over [0, 10]. The numerical performance of the algorithm was evaluated by generating problems where the measure˜ are FFUTOA. Nine transmitter direcment matrix D tions were simulated over the unit sphere. Measurements were then constructed according to (1). The errors over 1000 such experiments can be seen in Figure 1
Figure 2. Effect of quality of far filed approximation on number of receivers (r) and transmitters (t)
To evaluate the sensitivity of noise, white gaussian noise was added to UTOA measurements, simulated with transmitters 107 from the origin. 10 transmitters and 6 receivers were used. A graph showing the mean error of 100 experiments plotted against the standard deviation of the noise is shown in Figure 3. The reconstruction was also optimized under the function 2 P˜ min r,n,f,g dij − rT nj − fi − gj using a Gaussi
Newton like approach. The result is also shown in Figure 3.
−2
”Wearable Visual Systems” and the Swedish Research Council project ”Polynomial Equations and Geometry in Computer Vision”.
−3
References
−1
10
initialization opitmized
Error in Position
10
10
−4
10
−5
10
−5
10
−4
10
−3
10
−2
10
Standerd deviation
Figure 3. Effect of additive white gaussian noise.
3.2. Real Data An experimental setup using 8 microphones as receivers and random distinct manually made sounds as transmitters was used. The 19 sound sources were ≈30 m away from the receivers. Microphones were set in the corners of a cuboid of roughly 100 × 105 × 60 cm3 . The beginning of each sound were matched by a heuristic algorithm to create TDOA measurements. As microphones were recorded through a sound card, an artificial offset had to be added to each microphone to create UTOA measurements. Using the same real measurements, but adding new receiver offsets for each experiment, the algorithm was run on 100 experiments. A mean error of 13 cm was obtained. 85% of the errors are in the floorto-roof direction. This can be explained by the sounds all being made close to ground level and thus the transmitter directions will be close to being in a plane, giving poor resolution in floor-to-roof direction.
4. Conclusions We have presented a novel method for solving the, to the authors’ best knowledge, previously unstudied problem of initializing a sensor network using fully unsynchronized time of arrival data in a far field setting. Degenerate constellations are also characterized. By using real and simulated data we show the applicability of the method. Experimental validation also shows that the method is fairly robust to both noise and settings where the approximation that the transmitters are far from the receivers is not good.
5
Acknowledgements
The research leading to these results has received funding from the swedish strategic research project ELLIIT, the Swedish Strategic Foundation funded project
[1] S. T. Birchfield and A. Subramanya. Microphone array position calibration by basis-point classical multidimensional scaling. IEEE transactions on Speech and Audio Processing, 13(5), 2005. [2] R. Biswas and S. Thrun. A passive approach to sensor network localization. In IROS, 2004. [3] M. Brandstein, J. Adcock, and H. Silverman. A closedform location estimator for use with room environment microphone arrays. Speech and Audio Processing, IEEE Transactions on, 5(1):45 –50, jan 1997. [4] A. Cirillo, R. Parisi, and A. Uncini. Sound mapping in reverberant rooms by a robust direct method. In Acoustics, Speech and Signal Processing, 2008. ICASSP 2008. IEEE International Conference on, pages 285 – 288, 31 2008-april 4 2008. [5] M. Cobos, A. Marti, and J. Lopez. A modified srp-phat functional for robust real-time sound source localization with scalable spatial sampling. Signal Processing Letters, IEEE, 18(1):71 –74, jan. 2011. [6] H. Do, H. Silverman, and Y. Yu. A real-time srpphat source location implementation using stochastic region contraction(src) on a large-aperture microphone array. In Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference on, volume 1, pages I–121 –I–124, april 2007. [7] M. Hromcik and M. Sebek. Numerical and symbolic computation of polynomial matrix determinant. In Decision and Control, 1999. Proceedings of the 38th IEEE Conference on, volume 2, pages 1887 –1888, 1999. ˚ om. Under[8] Y. Kuang, E. Ask, S. Burgess, and K. Astr¨ standing toa and tdoa network calibration using far field approximation as initial estimate. In ICPRAM, 2012. [9] D. Niculescu and B. Nath. Ad hoc positioning system (aps). In GLOBECOM-01, 2001. [10] M. Pollefeys and D. Nister. Direct computation of sound and microphone locations from time-differenceof-arrival data. In ICASSP, 2008. [11] V. C. Raykar, I. V. Kozintsev, and R. Lienhart. Position calibration of microphones and loudspeakers in distributed computing platforms. IEEE transactions on Speech and Audio Processing, 13(1), 2005. [12] J. Sallai, G. Balogh, M. Maroti, and A. Ledeczi. Acoustic ranging in resource-constrained sensor networks. In eCOTS-04, 2004.
