3

I am having some trouble to find or implement an algorithm to find a signal source. The objective of my work is to find the sound emitter position.

To accomplish this I am using three vibration sensors. The technique that I am using is multilateration that is based on the time difference of arrival.

The time difference of arrival between each sensor are found using Cross Correlation of the received signals.

I already implemented the algorithm to find the time difference of arrival, but my problem is more on how multilateration works, it's unclear for me based on my reference, and I couldn't find any other good reference for this that are free/open.

I saw this post Trilateration using TDOA But I can't figure out how to solve the set of equations(7) of the wikipedia page of multilateration as i have only the three TDOA.

Any help on this would be much appreciated

Community
  • 1
  • 1
Youssef Kamoun
  • 43
  • 1
  • 8

2 Answers2

5

You have three sensor coordinates A,B,C, unknown coordinate of signal source P, unknown time of signal start t0, and three times of signal registration ta, tb, tc.

Example: Let's sensor A caught a signal in 12:00:05, sensor B - in 12:00:00, sensor C - 12:00:07. So assign time differences: ta=5, tb=0, tc=7

Squared distances from sensors to source correspond to times of signal walk with velocity v (speed of sound in air or another environment)

(Ax-Px)^2 + (Ay-Py)^2 = (v*(ta-t0))^2   {1}
(Bx-Px)^2 + (By-Py)^2 = (v*(tb-t0))^2   {2}
(Cx-Px)^2 + (Cy-Py)^2 = (v*(tc-t0))^2   {3}

Open the brackets, then subtract equations ({2}-{1}, {3}-{2},{1}-{3}) to discard squares of unknown terms.

 Ax^2-2*Ax*Px + Px^2 + Ay^2-2*Ay*Py + Py^2 = v^2*(ta^2 - 2*ta*t0 + t0^2)
 Bx^2-2*Bx*Px + Px^2 + By^2-2*By*Py + Py^2 = v^2*(tb^2 - 2*tb*t0 + t0^2)
 Cx^2-2*Cx*Px + Px^2 + Cy^2-2*Cy*Py + Py^2 = v^2*(tc^2 - 2*tc*t0 + t0^2)

 Bx^2-Ax^2 -2*(Bx-Ax)*Px + By^2-Ay^2 -2*(By-Ay)*Py = v^2*(tb^2-ta^2 -2*(tb-ta)*t0)    {1'}
 Cx^2-Bx^2 -2*(Cx-Bx)*Px + Cy^2-By^2 -2*(Cy-By)*Py = v^2*(tc^2-tb^2 -2*(tc-tb)*t0)    {2'}
 Ax^2-Cx^2 -2*(Ax-Cx)*Px + Ay^2-Cy^2 -2*(Ay-Cy)*Py = v^2*(ta^2-tc^2 -2*(ta-tc)*t0)    {3'}

Now you have system of three linear equations with three unknowns. It might be solved with some widespread algorithms - Gauss elimination, LU decomposition etc.

Note that solution precision strongly depends on small errors in coordinates and time measurements (this method is not very robust).

MBo
  • 77,366
  • 5
  • 53
  • 86
  • Is velocity v = peak value ? ,and for t0 how do I get it? and Ta= TDOA of AB and Tb of BC and TC of AC ? – Youssef Kamoun Mar 23 '16 at 11:41
  • v is sound speed. t0 is unknown value, you will found it as aside effect. ta,tb,tc are moments of signal registration. You can assign zero to the earliest time, and differences to this moment to others. Example added. – MBo Mar 23 '16 at 11:55
  • Okay then for me i have Tbc = 88ms, Tab =44ms and Tac=44ms It is right? i found these by cross correlation between a,b,c signals ,so t0 and x,y of P ,i have to find three of them or t0 will be deleted after substracting different equations? then I have vibration sensor,does it have also velocity like speed of sound ? and If it is indeed ,in normal environment means in the air ,is it like something uniform? Thanks and sorry for troubling you :) – Youssef Kamoun Mar 23 '16 at 12:15
  • If both vibration sensor and sound source are mounted at the hard ground (like earthshake detector), you have to consider sound speed in this ground. And it is not clear for me - what sensor was the first? – MBo Mar 23 '16 at 12:25
  • it is like both of them on the table and vibration source is like tapping on the table . the first was the B so should it be tb =0 ,ta = 44 and tc =44? – Youssef Kamoun Mar 23 '16 at 12:29
  • Yes, here you have to use sound speed in solid table plate material. I wonder - if `tb =0 ,ta = 44,tc =44` - why Tbc = 88? May be `tb=0, ta=44, tc=88`? – MBo Mar 23 '16 at 12:36
  • Ah Yes as you said tb=0,ta=44 and tc=88 , then for t0 ,what do you mean ,i will find it as asdie effect?and the sound speed in solid table plate material is universal or would i have to calculate it? – Youssef Kamoun Mar 23 '16 at 12:39
  • I found here http://www.engineeringtoolbox.com/sound-speed-solids-d_713.html that the speed of sound Wood 3300 - 3600 (m/s) ,as the table is in wood,is it the v that i should use? – Youssef Kamoun Mar 23 '16 at 12:58
  • t0 could be calculated during system solution, but you probably don't need it. Yes, you can use this speed – MBo Mar 23 '16 at 14:09
  • Thanks,gonna try it and tell you how it went ;) – Youssef Kamoun Mar 23 '16 at 14:15
  • For Tac=308ms,Tab=44,Tbc=132 ,knowing that it arrived first at B ,i tried like Tb=0,Ta=44,Tc=132 and i find t0=0 and let's not talk about x or y that i found ridiculous values,is it correct For the T ? – Youssef Kamoun Mar 24 '16 at 09:40
  • Strange - for Ta=44,Tc=132 - Tac should be 88? – MBo Mar 24 '16 at 10:49
  • You are right ,now i have made many test , I have Ta=0,Tb=129 and Tc=86 ,The three equations up after some modification will be like this : -2*(BX-AX),-2*(BY-AY), 2*V*V*(Tb-Ta)= V*V*((Tb*Tb)-(Ta*Ta))-(BX*BX)+(AX*AX)-(BY*BY)+(AY*AY) -2*(CX-BX),-2*(CY-BY),2*V*V*(Tc-Tb)=V*V*((Tc*Tc)-(Ta*Ta))-(CX*CX)+(BX*BX)-(CY*CY)+(BY*BY) -2*(AX-CX),-2*(AY-CY),2*V*V*(Ta-Tc)=V*V*((Ta*Ta)-(Tc*Tc))-(AX*AX)+(CX*CX)-(AY*AY)+(CY*CY) But when i do computations like this i get a singular matrix where det=0 , the equations are good,right? – Youssef Kamoun Mar 30 '16 at 15:32
  • Is it because of the Time difference that that matrix is singular or there is an error with the equations? – Youssef Kamoun Mar 31 '16 at 13:31
  • Any help would be appreciated. – Youssef Kamoun Apr 04 '16 at 11:16
  • Sorry, I elaborated all the math for this problem from scratch. I believe that this is real help. But I have not a lot of time to do all work for you. I recommend to make simulation - get some position, calculate theoretic times for it, use these times in equations - to see possible mistakes. – MBo Apr 04 '16 at 11:26
  • I tried everything but still gives no solution. i want to get at least 10 or 20cm near the point ,it gives me really bad solutions like 20km near the point ^^ – Youssef Kamoun Apr 04 '16 at 11:45
  • Short example with real table sizes and approx. sound speed: A=(0.0),B=(0.5,0),C=(0,0.4) (in meters). P=(0.3,0). distPB=0.2m, tb=0.2/2000=10^-4s=100 microseconds (not milli!), ta=150mks, tb-ta=50mks. Are you sure that your delays are in milliseconds? – MBo Apr 04 '16 at 11:58
  • Yes from a python script the function millis which return ms – Youssef Kamoun Apr 04 '16 at 13:56
  • normally we search for P ,right? and with just X and Y so P (X,Y) ,what is the 2000 to find Tb? i don't understand your example.. – Youssef Kamoun Apr 04 '16 at 13:58
  • Time to signal walk from P to B is t=d/v = 0.2 meters / 2000 meters per second = 100 **micro**seconds – MBo Apr 04 '16 at 14:25
  • But i don't have distPB,i have Ta,Tb,Tc so Tp = d/v ,which d?it is getting confused in here – Youssef Kamoun Apr 04 '16 at 15:16
  • this is simple model simulation. d is distance – MBo Apr 04 '16 at 15:49
  • What would this look like in three dimensions? – 10GeV Nov 15 '20 at 01:27
  • @KeithMadison I think - similar. One more sensor, z-coordinate added, 4 equations. – MBo Nov 15 '20 at 03:15
1

Geometrically, a hyperbola represents the cloud of points with a constant difference in distance between two points. You have 3 points, but taken pairwise, the time differences between the 3 possible pairs will allow you to draw 3 hyperbolas. Look for a spot at or between where the hyperbolas intersect on a plot. Or solve the equivalent algebra (least squares).

hotpaw2
  • 70,107
  • 14
  • 90
  • 153