Impact location, using four sound sensors.

Honesty, i have lost you at the first part already........ =(
But i will spend some time on this and see if i can get some
understanding of how this all fits together..

Thank you for the effort you put in it is appreciated..

robtillaart:
--update --
Sorry, my first reasoning was far too simplified ==> removed; a retry.

You need the distance between the microphones and the two delta T arrival times of the shockwave to do the math.

If you have three points A, B and C.
arrival times of soundwave (in order)
A - T0
B - T1
C - T2

The fact that the soundwave arrived at A first defines an area(1) of all points P: d.PA < d.PB and d.PA < d.PC (d.PA = distance PA)

The delta-time AB = T1 - T0 defines a distance d1 = (T1 - T0) /so (so = speed Sound)
The delta-time AC = T2 - T0 defines a distance d2 = (T2 - T0) /so

define the curve of all points Q: d.QA - d.QB = d1 (.QA = distance QA)
define the curve of all points R: d.RA - d.RC = d2

These two curves cross each other in area (1) => the point of impact ==> Q == R.

-- update --
if you knew the time of impact the d.QA, dQB and d.QC would be known, making the math faaaaar simpler as these curves are not trivial - quadratic asymptotic beasts with sqrts in it - :frowning:

Think it is easier to write an approximating algorithm that searches the point.
The fact that point A heard the soundwave first => d.QA < d.QB

-- update 2 --

from: - Formula and graph of a hyperbola. How to graph a hyperbola based on its formula

A hyperbola is a set of all points P such that the difference between the distances from P to the foci, F1 and F2, are a constant K

so my "quadratic asymptotic beasts with sqrts in it" can be rewritten as hyperbola with A and B (A & C) as foci.

using the drawing of the webpage above:
Assume A = (0,-c) and B = (0,c) and the constant K = d1 = (T1 - T0) /so. The point (0,-a) where the hyperbola crosses the Y-axis is (0, -d1/2)

--> Focus of a Hyperbola

To determine the foci one uses a^2 + b^2 = c^2 => b^2 = c^2 - a^2 = (d.AB/2)^2 - (d1/2)^2

The formula of the hyperbola becomes : y^2 / (-d1/2)^2 - x^2 / (d.AB/2)^2 - (-d1/2)^2 = 1

Same trick for the points A & C (hint: it is easier to use another reference framework to determine the formula and do a translation afterwards: X -> X-xdelta Y -> Y-ydelta)

TODO: determine intersection points of the two hyperbolas and then your close...

-- update 3 --

intersection points
-- http://www.analyzemath.com/HyperbolaProblems/hyperbola_intersection.html

Difference with the location problem is that the hyperbola defined by points AB and the one defined by points AC are 'orthogonal' - think of it as the red in the drawing rotated 90 degrees (make a drawing!!) There will be two intersection points and because the soundwave arrived first at A it becomes obvious which one to choose.

The code is left as an exercise ....

Another line of thought.
When the soundwaves arrive the soundwave has a definite form which will be slightly different at the three points ABC. The essence is that the peakvolume of the soundwave is proportional (linear, quadratic or otherwise) with the distance. The farther away the weaker the sound. As one could calibrate the micro's with a signalstrength/distance table the math would become simpler again.

Don't know if the differences are within the noise-level of the signals/ micro's/arduino ADC but a simple test could reveal this. Furthermore soundwaves can be deformed by obstacles in the open field etc. Still it has some potential worth investigating.

Question: How big is ABCD in square meters? smaller/bigger?

Some additional hyperbola math: - http://www.codecogs.com/reference/maths/geometry/coordinate/hyperbola.php

Hi Rob,

i have been working on the method you sent me most of the day, but something
just caught my attention..
Using the Hyperbola will give me a XY coordinate on the Hyperbola's Dimensions A-B, not on the the larger
grid..
I still have no Idea mathematically where on the target area A-B is located, so those coordinates
results in another set of unusable data..

Does this make any sense to you ?
Or am i missing something here ?

robtillaart:
--update --
Sorry, my first reasoning was far too simplified ==> removed; a retry.

You need the distance between the microphones and the two delta T arrival times of the shockwave to do the math.

If you have three points A, B and C.
arrival times of soundwave (in order)
A - T0
B - T1
C - T2

The fact that the soundwave arrived at A first defines an area(1) of all points P: d.PA < d.PB and d.PA < d.PC (d.PA = distance PA)

The delta-time AB = T1 - T0 defines a distance d1 = (T1 - T0) /so (so = speed Sound)
The delta-time AC = T2 - T0 defines a distance d2 = (T2 - T0) /so

define the curve of all points Q: d.QA - d.QB = d1 (.QA = distance QA)
define the curve of all points R: d.RA - d.RC = d2

These two curves cross each other in area (1) => the point of impact ==> Q == R.

-- update --
if you knew the time of impact the d.QA, dQB and d.QC would be known, making the math faaaaar simpler as these curves are not trivial - quadratic asymptotic beasts with sqrts in it - :frowning:

Think it is easier to write an approximating algorithm that searches the point.
The fact that point A heard the soundwave first => d.QA < d.QB

-- update 2 --

from: - Formula and graph of a hyperbola. How to graph a hyperbola based on its formula

A hyperbola is a set of all points P such that the difference between the distances from P to the foci, F1 and F2, are a constant K

so my "quadratic asymptotic beasts with sqrts in it" can be rewritten as hyperbola with A and B (A & C) as foci.

using the drawing of the webpage above:
Assume A = (0,-c) and B = (0,c) and the constant K = d1 = (T1 - T0) /so. The point (0,-a) where the hyperbola crosses the Y-axis is (0, -d1/2)

--> Focus of a Hyperbola

To determine the foci one uses a^2 + b^2 = c^2 => b^2 = c^2 - a^2 = (d.AB/2)^2 - (d1/2)^2

The formula of the hyperbola becomes : y^2 / (-d1/2)^2 - x^2 / (d.AB/2)^2 - (-d1/2)^2 = 1

Same trick for the points A & C (hint: it is easier to use another reference framework to determine the formula and do a translation afterwards: X -> X-xdelta Y -> Y-ydelta)

TODO: determine intersection points of the two hyperbolas and then your close...

-- update 3 --

intersection points
-- http://www.analyzemath.com/HyperbolaProblems/hyperbola_intersection.html

Difference with the location problem is that the hyperbola defined by points AB and the one defined by points AC are 'orthogonal' - think of it as the red in the drawing rotated 90 degrees (make a drawing!!) There will be two intersection points and because the soundwave arrived first at A it becomes obvious which one to choose.

The code is left as an exercise ....

What kind of distance are you expecting between the receivers? Will they all be connected to the uno via wires?

Yes they will be connected via wire to the Uno,

My current test is based on 1200mm x 800mm but later on i would like to scale it up
more to about 2.5m x 2m...

andywatson:
What kind of distance are you expecting between the receivers? Will they all be connected to the uno via wires?

My current test is based on 1200mm x 800mm but later on i would like to scale it up

Did you make considerations about the precision you'll obtain?
Using the standard max sample rate of the Arduino ADC (prescaled with 128) you will get these readings AFAIK:

4 sensors dividing an expected sample rate of 8.9 Khz:

340m/s speed of sound (22[sup]o[/sup]C)= 340000mm/s
340000/(8900/4) ~= 153 mm.

This means a medium sampling distance of sound vs distance of about 15cm for a single sensor and of about 4 cm for the single fronts delta between adjacent sensors (which makes the error considerations tricky).

Being new to Arduino, i have not really considered that, but i have read quite a few references
of using external timing methods to improve the accuracy..

Honestly at this moment i am less concerned about the accuracy than figuring out the
mathematics to solve this problem..

The hardware is something i will work on fine tuning later if needed, that is luckily the part where there is
way more information available to improve the solution....
The math seems to be a challenge any way i look at it.

I have read a post on this forum of a Guy using a similar setup who achieved something like
1.4mm accuracy, now if that is true it is more than i need..
unfortunately i need to find the location before i can determine the accuracy of the solution
and start working on that problem...

Is there anyone out there who has done something similar, who can help me PLEASE =(

scjurgen:

My current test is based on 1200mm x 800mm but later on i would like to scale it up

Did you make considerations about the precision you'll obtain?
Using the standard max sample rate of the Arduino ADC (prescaled with 128) you will get these readings AFAIK:

4 sensors dividing an expected sample rate of 8.9 Khz:

340m/s speed of sound (22[sup]o[/sup]C)= 340000mm/s

340000/(8900/4) ~= 153 mm.



This means a medium sampling distance of sound vs distance of about 15cm for a single sensor and of about 4 cm for the single fronts delta between adjacent sensors (which makes the error considerations tricky).

played with the math today and couldn't figure out a nice formula, so I tried another approach, simulation.

When an impact is made one knows three time differences (these are simulated in the code below). In the simulation I try (not) all points of the grid:

  • if this point is the point of impact would it (re)create the observation made?

Based upon a grid of 1000x1000 the first approach (brute force) worked well but too slow. The current version starts with the middle of the grid and searches for the points with a lower error around it. This repeats until no better point is found (and the error == 0).

The timing improved dramatically, brute force took several minutes while the current algorithm is in order of 100 msec on my 328 - at least as far as I tested.

Comments and improvements are as allways welcome.

  • update -
    patched a bug;
    timing now averaging around ~100 msec (sometimes one pixel off)
//
//    FILE: impactSearch.pde
//  AUTHOR: Rob Tillaart
//    DATE: 2011-02-21 
//
// PURPOSE: determine impact location based upon arrival time of signal at 4 sensors
//
 
#define SQR(x) ((x)*(x))

// THE GRID = 1000 x 1000; 
// p[] = point A, B, C, D
const long p[8] = { 0,0, 0,1000, 1000,0, 1000,1000 };

int best_x = 0;  // point to be found
int best_y = 0;  // idem
float best_SE = 0;   // Square Error

unsigned long start = 0;  // timer

void setup()
{
  Serial.begin(115200);
  Serial.println("Start...");
}

void loop()
{
  // simulation random impact near point A then B then c then D
  // as the grid is 1000x1000 the max location is 500-500 as otherwise 
  // it would be closer to B, C or D
  long px = random(499);
  long py = random(499);

  Serial.print(px);
  Serial.print(",");
  Serial.print(py);
  Serial.print("\t\t");

  // The distances to the 4 mikes are calculated
  // and the value of TA is subtracted as there the sound arrives first.
  // Note distance is time * speed of sound
  // the values here are calculated 
  // one could add random noise for the simulation
  float TA = sqrt(SQR(px) + SQR(py));
  float TB = sqrt(SQR(px) + SQR(1000L - py)) - TA;
  float TC = sqrt(SQR(1000L-px) + SQR(py)) - TA;
  float TD = sqrt(SQR(1000L-px) + SQR(1000L - py)) - TA;
  
  start = millis();
  impactSearch(TA, TB, TC, TD);    
  Serial.print(millis() - start);  // how fast
  Serial.print("\t");
  
  // print delta + error
  Serial.print(best_x - px); 
  Serial.print(",");
  Serial.print(best_y - py); 
  Serial.print("\t");
  Serial.println(best_SE);
}

// impactSearch detemines the point of impact by searching
// the point with a smart trial and error method
// in essence:
// It first takes the middle and tests if that was the point of 
// impact what would be the timing and compares that to the factual 
// times. Then it does the same for the points around it and 
// moves from point to point decreasing the error. 
// if there is no better point, the point of impact has been found
//
// This code uses an optimization by doing big steps in the beginning
// and decreasing the stepsize until 0 in the end
//
// Some optimizations are under investigation
// e.g. 3 mikes seems to be enough, but when noise is added
// a fourth mike reduces the overall error I think.
void impactSearch(float t1, float t2, float t3, float t4)
{
  // start in the middle 
  // in fact (250,250) is a better point to start as this
  // is the middle of the grid 
  best_x = 500;
  best_y = 500;
  best_SE = 100000000L;
  int step = 16;
  
  boolean found = false; // point not found yet
  while (false == found)
  {
    boolean decreaseStep = true;
    
    for (int x = best_x-step; x <= best_x+step; x+=step)
      for (int y = best_y-step; y <= best_y+step; y+=step)
      {
        // optimization, same point is never better
        if (y== best_y && x==best_x) continue;

        // determine what would be the arrival time
        // for point x,y
        // DA = distance to point A - etc
        float DA = sqrt( SQR(p[0]-x) + SQR(p[1]-y) );
        float DB = sqrt( SQR(p[2]-x) + SQR(p[3]-y) );
        float DC = sqrt( SQR(p[4]-x) + SQR(p[5]-y) );
        float DD = sqrt( SQR(p[6]-x) + SQR(p[7]-y) );
        
        // use square error when compared to the real arrival times.
        float se = SQR(DA + t2 - DB);
        se += SQR(DA + t3 - DC);
        se += SQR(DA + t4 - DD);
        
        // remember the one with the smallest error
        if (se < best_SE)
        {
          decreaseStep = false;
          best_SE = se;
          best_x  = x;
          best_y  = y;
        }
      }
    // if no better point found decrease the search area
    // by decreaing the step size
    if (decreaseStep) step--;  // was step = step/2;
    found = (step == 0);
  }
}

Nice, I like that approach, being a non-linear problem it seems to be a good idea.
Some considerations. The changes between to guessed positions are 'near' linear. So, a kind o 2-D triangled bsearch should be possible (I think a Newton approximation using the f'() is overkill).

Another point is the first guess. There are physical minimum and maximum values, some are trivial. If the 2 deltatimes (dt) are the same (impacts: 0,dt,dt) they must lie on the diagonal between the 2 other sensors. Therefore, the distance should be a trivial 500-dtsqrt(2)/2 (sound travels linear in time). On the other hand, if one delta is 0 (impacts: 0,0,dt), it is on a line that is orthogonal to the 2 sensors. And the position should be between 500 (dt near 0) and 0 for dt=5002/(1+sqrt(5)). Interpolating linear between these 2 should be a good first guess (I guess).

Ok, here a numerical solution for your problem. Couldn't find an elegant solution that resolves this complex system (but still, it looks like it is possible). 8)

Some facts: you need 4 sensors, as one does only determine the first time impact, so you actually have only 3 data values, and as someone maybe remembers, with 3 values you can determine a triangle, not with 2. The following algorithm proves that (generating a lot of solutions if you don't consider the 3rd value, and they are all ok :wink: ).

The algorithm is based on this idea.
Not knowing the time to first impact I assume that I know it anyway, I calculate for a guessed value what x y should be, knowing two delta times (easy application of law of cosines). The third delta time cross checks the initial guess of t. Varying the guess of t you find the minimum error.
This small program has also a routine for testing some points on a grid.

What you need to do:

  1. ordering of impacts and mirroring the result values for the 4 different quadrants.
  2. feed the sound data to the guessing routine with the sound correction constant (speed of sound + temperature.
  3. it is highly optimizable (i think in 20 steps you can obtain an error of 0.01%, not sure, if I have time I will try that, it doesn't look necessary)
#include <stdlib.h>
#include <math.h>
/* quick and dirty code, not written on the arduino
*/

double a=1000,b=700; // this is all ugly, better create a class for the whole thingy


/* equations system is:
x=(p^2-t^2-a^2)/(-2*a), // cosine law, p and q the t1+t and p=t2+t
y=(q^2-t^2-b^2)/(-2*b)
*/
#define SQR(x) ((x)*(x))

void getxy(double t, double p, double q, double *x, double *y)
{
	*x=(SQR(q)-SQR(t)-SQR(a))/(-2*a);
	*y=(SQR(p)-SQR(t)-SQR(b))/(-2*b);
}

void findXY(double t1, double t2, double t3, double eps, double *xres, double *yres)
{
	double t,x,y;
	double mineps=sqrt(a*a+b*b);

/* this can be highly optimized !!! 
  the minimum t depends on the maximum time on the diagonal of the opposite sensor (or something like that)
*/

	for (t=0; t < sqrt(a*a+b*b); t+=0.1)
	{
		getxy(t,t1+t,t2+t,&x,&y);
		double t3e=sqrt(SQR(a-x)+SQR(b-y)); // cross check with opposite point (latest arrival)
		// printf("t=%f x=%f y=%f t3=%f eps=%f\n ", t,x,y, t3e, fabs(t3+t-t3e));
		if (fabs(t3+t-t3e) < mineps)
		{
/* you can optimize this too, once we have a local minimum for 3 values (epsilon[0]>epsilon[1] < epsilon[2]), we can refine the search  (t+=0.01 starting with t-0.1 until t+0.1) or exit 
*/
			mineps=fabs(t3+t-t3e);
			*xres=x;
			*yres=y;
		}
	}
}

void test(double x, double y)
{
	double t1=sqrt(SQR(x)+SQR(y));
	double t2=sqrt(SQR(x)+SQR(b-y));
	double t3=sqrt(SQR(a-x)+SQR(y));
	double t4=sqrt(SQR(a-x)+SQR(b-y));

	printf("%.1f %.1f %.1f %.1f: ", t1,t2,t3,t4);

//	printf("%f x ", (SQR(t2)-SQR(t1)-SQR(b))/(-2*b));
//	printf("%f\n", (SQR(t3)-SQR(t1)-SQR(a))/(-2*a));

	t2-=t1;
	t3-=t1;
	t4-=t1;
	t1=0;
	printf("t1=%.1f t2=%.1f t3=%.1f ", t2,t3,t4);
	findXY(t2,t3,t4,0.01,&x,&y);
	printf("~ {%.2f;%.2f}\n",x,y);
}


void testgrid()
{	
	for (double y=0; y <= b/2; y+=50.0)
	{
		for (double x=0; x <= a/2; x+=50.0)
		{
			test(x,y);
		}
	}
}

int main(int argc, char *argv[])
{
	testgrid();
}

Thanks Juergen,

I will start playing with this one aswell,
Thank all you guys again, i appreciate all the effort you put in to try and help resolve this
problem..

Hopefully one day i can also contribute as mush to your or someone else s project..
Might take time but hope to get there...

Thanks Rob and Jeurgen.. i will test and post my results..

Regards
Marcel

Hi Gentlemen,

I have a new question... relating to the same topic...
I have been doing some reading on the various ways to get the audio data
into the Arduino.

Well i am at the point where i am not sure what is the correct direction to take..

Option one:
Using the pre-ams on the analog ports as i am doing now.

Option two:
Using a OpAmp as a comparator and putting the data into the Digital ports,
and using one of the Hardware timers to create an interrupt to signal the arrival
of a new input..

Would Digital be faster than Analog, or the other way around ?

Does anyone have any advice regarding these options, it seems the prediction that
i will loose accuracy due to the speed of the processor.
so any improvements i can make will help me get more reliable data.

currently i get very random data, still trying to figure it out why, that is why i am looking at way to improve that
data input...

here is what i currently get, five impacts at roughly a similar location on a surface 940mm x 640mm
the times are in millisecond....

------1---------
Time1 = 0.00
Time2 = 0.94
Time3 = 2.84
Time4 = 5.69

------2---------
Time1 = 0.00
Time2 = 2.34
Time3 = 4.69
Time4 = 13.19

------3---------
Time1 = 0.00
Time2 = 3.75
Time3 = 4.70
Time4 = 5.65

------4---------
Time1 = 0.00
Time2 = 0.47
Time3 = 1.88
Time4 = 3.78

------5---------
Time1 = 0.00
Time2 = 3.76
Time3 = 4.70
Time4 = 8.03

Any Ideas????

Thank you
Marcel

Hyperbolas are the way to go for this, because what you have is differences in arrival times, not absolute times. You know the arrival times at all four locations, t0, t1, t2, t3. That gives you six differences: (t3 - t2), (t3 - t1), (t3 - t0), (t2 - t1), (t2 - t0), (t1 - t0); each of these generates a hyperbola. The intersection of six hyperbolas gives you one and only one point.

Four is the minimum number of sensors this will work with, absent other constraints. You need a minimum of six hyperbolas to get a unique intersection point. If you have three sensors, you get three hyperbolas, which gives you two intersections points. If you can rule one out (i.e. you know the impact is in a given area), you're good, but four sensors gives you a unique solution as long as the placement of the sensors does not give degenerate geometry.

Marceldv:
I have been doing some reading on the various ways to get the audio data into the Arduino.
Well i am at the point where i am not sure what is the correct direction to take..

Option one:
Using the pre-ams on the analog ports as i am doing now.

Option two:
Using a OpAmp as a comparator and putting the data into the Digital ports,
and using one of the Hardware timers to create an interrupt to signal the arrival
of a new input..

Would Digital be faster than Analog, or the other way around ?

Digital is way faster, and the routines of the standard digital acquisition can be optimized.
I did an analysis (Arduino 2009) some time ago using different techniques (and then I switched the micro :wink: ).

digitalRead non optimized: 165672 samples/second
digitalRead optimized: 265111 samples/second

analogRead not optimized: 8927 samples/second
analogRead prescale=64: 16603 samples/second
analogRead prescale=32: 31250 samples/second
analogRead prescale=16: 52521 samples/second
analogRead prescale8: 82987 samples/second

These are crude times in a simple loop, meaning, not considering eventual processing. Analogread can be very much enhanced if you split the request for a value and the reading (request value, do other stuff, read value)
the routines for changing analog reads you find here (you will need to change the code I guess, it was for reading 8 channels circulary)
http://www.schwietering.com/jayduino/fasterAnalogRead/

Also consider in your code that the sampling of the data channels is not at the same moment and introduces a fixed delay between 2 channels (t=1/samplerate) that you need to compensate for on each channel when you detect a signal.

You also need to consider doing these actions in a timer interrupt, your timing must be impeccable and not disturbed by processing. So, in the interrupt you get the value, save it to a ring buffer and get quickly out of there. The processing you do in the main loop by reading the data in the ring buffers and doing what you need to do.

HTH

Digital is way faster, but it requires you to put all of the things you need to discriminate between an impact and other ambient sounds. There is a guy here in Seattle by the name of Tangent makes a nifty little OSHW pendant that flashes some LEDs in response to the beat of ambient music and he appears to have tackled most of the relevant issues. There are some kits in the vending machine at Metrix, but I can't seem to find a link to his website.

Would i be able to implement some form of an interrupt on analog sensors ?

You can set an interrupt on the acquisition of a sample. Then you can decide whether to fall through to whatever else you want to do.

As far as I know only digital signal can trigger IRQs so you need to make electronics that makes a square wave of the analog signal when a threshold (soundlevel) is reached.

The Arduino 328/UNO has only two direct IRQ lines but also has a collective IRQ on the other pins. So you need to determinew which pin was triggered. See - Arduino Playground - PinChangeInt

Another option could be using the a/d converter in free running mode. Than at least the samples are taken with a known interval.

Following this thread with interest,

Jeroen

if you have 4 known locations (x1,y1, x2,y2, x3,y3, and x4,y4),
3 times after the trigger event, and you want to determine x0,y0 that is the location
of the sound source then you have 7 equations and 6 variables that means that the task is
solvable.
the equations are:
(x_i-x0)^2 + (y_i-y0)^2 = R_i^2 for i=1..4

  • 3 equations of type:
    R_i-R_earliest = v_sound*(t_i-t_earliest)=v_sound*t_i since t_earliest=0 as it syncs everything.
    R_earliest is one of R_i, each time it's a different one but it doesn't change anything

rewrite this system to express x0 and y0 and you're done.