XIE et al:SYNCHRONIZE INERTIAL READINGS FROM MULTIPLE MOBILE DEVICES IN SPATIAL DIMENSION 2153 Positive Time Interval 80 20 Axs:g' Ag8:日 Fig.11.Identify the positive time interval according to the PCA results. Fig.12.Estimate the rotation matrix At.t+6t. those acceleration measurements only from the positive time intervals to extract the accumulated positive magnitudes of According to Theorem 2,our solution of estimating the the consistent acceleration fe(t).According to the formulation direction of fe(t)is as follows:First,we select a sufficiently in Eq.(1),for the accumulated acceleration measurements, large time interval to obtain the corresponding PCA results, the accumulated positive magnitudes of f(t)dominates the and identify the positive time intervals from the PCA results result than the accumulated magnitudes of the inconsistent where the accelerations are greater than 0.Then,for each axis acceleration f(t),which is very close to 0.In this way, of the reference local coordinate,we compute the expected we are able to further estimate the direction of fc(t)in value of the acceleration measurements within the above posi- its local coordinate,according to its positive magnitudes in tive time intervals.After that,we can compute the angles oi,Bi different axes.Therefore,we can estimate the direction of()andaccording to Eq.(),and estimate the directionoff(t). according to Theorem 2. As the value of E[fe(t)]is the same in all three formulations Theorem 2:Let i,Bi and Yi be the constant angles in Eq.(8),it is actually unnecessary to compute the value between the consistent acceleration fe(t)and the z,y and z of E[fe(t)].As a matter of fact,according to the previous axes in the reference local coordinate,respectively.Let fe(t) analysis,we have vE[fi(t)]2 +E[fi.(t)]2 +E[fi.2(t)2 be the magnitude of fe(t),and fi(t),fi.(t)and fi.=(t)be E[fe(t)]. the acceleration measurements in the z,y,and z-axes from the accelerometer,respectively.The direction of f(t)can be VII.GYROSCOPE BASED ORIENTATION TRACKING represented by a unit vector (cos ai,cos Bi,cos i),where the A.Rotation Based Modeling value of cos ai,cos Bi,and cosYi can be estimated as follows: Our solution is based on the observation that the extracted E[fi.(t)] cos ai= direction of gravity is stable in the global coordinate.This ELf(t切 observation can be further leveraged to calibrate the rotation cos Bi= E[fi,y(t】 (8) matrix A.t+t derived from gyroscope based tracking.For E[fe(t)] device Di,suppose that the gravity is denoted as a vector E[fi.z(t)] g=(gz(t),gu(t),g=(t))in the local coordinate system Li(t) cos= E[fe(t)] at time t.Then,after a period of tracking,the gravity could be denoted as a different vector g'=(gz(t +6t),gu(t+ Proof:Without loss of generality,we prove the equation EIfi.z(t)] 6t),g2(t +6t))in the local coordinate system Li(t +6t)at cosai According to Eg.(1).for the z-axis, DLfe(t打 time t+6t.The difference between g and g'is caused by fi.(t)=fe(t)cosi+f(t).After we obtain a sequence of the rotation of the body frame from time t to time t+ot. fi.(t)from the positive time intervals in a sufficiently large we denote the corresponding rotation matrix as A.t+.Thus time window Wp,the sum of fi.(t)in these positive time we have g'=At.t+6g.According to the above relationship, intervals,ie..tew(t).can be depicted as follows: the rotation matrix At.t+6t can be computed as follows.As ∑fi()=∑[f():cos剑+∑fi(t). shown in Fig.12,note that due to the rotation of the body frame,there may exist a non-zero angle between the vectors tEWp tEWp g'and g.it can be calculated as:-arccos().where According to Theorem 1.tewf(t)0.Besides, the.operation refers to the inner product of the two vectors. since fe(t)is extracted from the positive time intervals,then During the process of rotation,to align the vector g to the ∑tew,f.()·cos]≥∑tew,fi.().Thus, vector g',it is essential to rotate axis u by an angle 0.The rotation axis u is orthogonal to the plane where the vectors ∑f.()≈∑[f.()cos创. g and g'lies on.Then,the rotation axis u can be obtained tEW by computing the cross product of g and g':u EXg' g×gT Moreover,since cosa;is constant,then, We call the above rotation as the vertical rotation.Therefore, ∑tew,fiz( according to the Rodrigues'rotation formula [28].we compute Cos≈ E[fi.=(t)] the corresponding rotation matrix R as: ∑tew,f.(O E[fe(t)] R=I+(sin0)U+(1-cos0)U2 (9) Similarly,we can prove that cos Elfi(t)] and Efe(t) Ef.=(t)】 where I is a 3 x 3 unit matrix,and U denote the cross-product c0s7= E fe(t) matrix for vector u.XIE et al.: SYNCHRONIZE INERTIAL READINGS FROM MULTIPLE MOBILE DEVICES IN SPATIAL DIMENSION 2153 Fig. 11. Identify the positive time interval according to the PCA results. those acceleration measurements only from the positive time intervals to extract the accumulated positive magnitudes of the consistent acceleration fc(t). According to the formulation in Eq.(1), for the accumulated acceleration measurements, the accumulated positive magnitudes of fc(t) dominates the result than the accumulated magnitudes of the inconsistent acceleration f i(t), which is very close to 0. In this way, we are able to further estimate the direction of fc(t) in its local coordinate, according to its positive magnitudes in different axes. Therefore, we can estimate the direction of fc(t) according to Theorem 2. Theorem 2: Let αi, βi and γi be the constant angles between the consistent acceleration fc(t) and the x, y and z axes in the reference local coordinate, respectively. Let fc(t) be the magnitude of fc(t), and fi,x(t), fi,y(t) and fi,z(t) be the acceleration measurements in the x, y, and z-axes from the accelerometer, respectively. The direction of fc(t) can be represented by a unit vector cos αi, cos βi, cos γi, where the value of cos αi, cos βi, and cos γi can be estimated as follows: ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ cos αi = E[fi,x(t)] E[fc(t)] cos βi = E[fi,y(t)] E[fc(t)] cos γi = E[fi,z(t)] E[fc(t)] . (8) Proof: Without loss of generality, we prove the equation cos αi = E[fi,x(t)] E[fc(t)] . According to Eq.(1), for the x-axis, fi,x(t) = fc(t) cos αi +f i,x(t). After we obtain a sequence of fi,x(t) from the positive time intervals in a sufficiently large time window Wp, the sum of fi,x(t) in these positive time intervals, i.e., t∈Wp fi,x(t), can be depicted as follows:  t∈Wp fi,x(t) =  t∈Wp [fc(t) · cos αi] +  t∈Wp f i,x(t). According to Theorem 1, t∈Wp f i,x(t) ≈ 0. Besides, since fc(t) is extracted from the positive time intervals, then t∈Wp [fc(t) · cos αi]  t∈Wp f i,x(t). Thus,  t∈Wp fi,x(t) ≈  t∈Wp [fc(t) · cos αi]. Moreover, since cos αi is constant, then, cos αi ≈ t∈Wp fi,x(t) t∈Wp fc(t) ≈ E[fi,x(t)] E[fc(t)] . Similarly, we can prove that cos βi = E[fi,y(t)] E[fc(t)] and cos γi = E[fi,z(t)] E[fc(t)] . Fig. 12. Estimate the rotation matrix At,t+δt. According to Theorem 2, our solution of estimating the direction of fc(t) is as follows: First, we select a sufficiently large time interval to obtain the corresponding PCA results, and identify the positive time intervals from the PCA results where the accelerations are greater than 0. Then, for each axis of the reference local coordinate, we compute the expected value of the acceleration measurements within the above posi￾tive time intervals. After that, we can compute the angles αi, βi and γi according to Eq. (8), and estimate the direction of fc(t). As the value of E[fc(t)] is the same in all three formulations in Eq. (8), it is actually unnecessary to compute the value of E[fc(t)]. As a matter of fact, according to the previous analysis, we have E[fi,x(t)]2 + E[fi,y(t)]2 + E[fi,z(t)]2 ≈ E[fc(t)]. VII. GYROSCOPE BASED ORIENTATION TRACKING A. Rotation Based Modeling Our solution is based on the observation that the extracted direction of gravity is stable in the global coordinate. This observation can be further leveraged to calibrate the rotation matrix A t,t+δt derived from gyroscope based tracking. For device Di, suppose that the gravity is denoted as a vector g = gx(t), gy(t), gz(t) in the local coordinate system Li(t) at time t. Then, after a period of tracking, the gravity could be denoted as a different vector g’ = gx(t + δt), gy(t + δt), gz(t + δt) in the local coordinate system Li(t + δt) at time t + δt. The difference between g and g’ is caused by the rotation of the body frame from time t to time t + δt, we denote the corresponding rotation matrix as At,t+δt. Thus we have g’ = At,t+δtg. According to the above relationship, the rotation matrix At,t+δt can be computed as follows. As shown in Fig. 12, note that due to the rotation of the body frame, there may exist a non-zero angle θ between the vectors g’ and g, it can be calculated as: θ = arccos( g·g’ g|g’ ), where the · operation refers to the inner product of the two vectors. During the process of rotation, to align the vector g to the vector g’, it is essential to rotate axis u by an angle θ. The rotation axis u is orthogonal to the plane where the vectors g and g’ lies on. Then, the rotation axis u can be obtained by computing the cross product of g and g’: u = g×g’ g×g’ . We call the above rotation as the vertical rotation. Therefore, according to the Rodrigues’ rotation formula [28], we compute the corresponding rotation matrix R as: R = I + (sin θ)U + (1 − cos θ)U2 (9) where I is a 3×3 unit matrix, and U denote the cross-product matrix for vector u.