User:Tedburke/Sandbox

= Draft Solution of the Tapping Localisation Problem =

Introduction
This is a draft solution for what we are currently referring to as the tapping localisation problem. The problem relates to an experimental human-machine interface which is under development. A rectangular board sits on a horizontal surface. At each corner of the board is an accelerometer which senses vibrations. When the board is "tapped", vibrations propagate outwards in all directions from the point of impact. Because the distance from the point of impact to each corner is different, the vibrations are detected by each accelerometer at a different time. Knowing only the dimensions of the rectangle and the time at which each accelerometer sensed the impact, the objective is to calculate at what point on the board the tap occurred. This is the tapping localisation problem.

It is easy to think of many other problems which resemble this one - e.g. finding the epicentre of an earthquake using readings from multiple measurement stations or identifying the position of an audio source using a microphone array - so it must certainly have been solved before. However, since the problem seemed like a straightforward geometrical one, we tried to solve it ourselves from first principles. Unfortunately, although we have a working solution (as described below), it is not an entirely satisfactory one. Part of the solution is carried out using an iterative numerical method - something which my instinct tells me should not be necessary. Nevertheless, a purely analytical solution has so far eluded us, so what follows will have to do for the moment.

Assumptions and Notation
The tapping area is assumed to be a rectangle of width $$w$$ and height $$h$$. As shown in the diagram below, the four corners of the rectangle are as follows:


 * A is the bottom left corner, at coordinates (0,0)
 * B is the bottom right corner, at coordinates (w,0)
 * C is the top right corner, at coordinates (w,h)
 * D is the top left corner, at coordinates (0,h)



At time $$t_0$$, a "tap" occurs at a point $$p=(x,y)$$ somewhere inside the rectangle. There is one accelerometer at each corner of the rectangle. Each accelerometer detects the tap only when the resulting vibrations have propagated out from $$p$$ as far as its location. The speed of propagation, $$v$$, is assumed to be constant throughout the medium.

Initially, $$v$$ is not known. Also, when a tap occurs, the time of impact is not known initially. All that we have to work with are the times at which the tap's vibrations reach each of the accelerometers. These four times are $$ t_a, t_b, t_c$$ and $$t_d$$.

Obviously, $$ t_a, t_b, t_c$$ and $$t_d$$ are greater than $$t_0$$ (i.e. no accelerometer can detect the tap before it occurs). Furthermore, the propagation speed, $$v$$, is constant and positive.

Finally, $$ a, b, c$$ and $$d$$ are the distances from $$p$$ to corners $$ A, B, C$$ and $$D$$ respectively.

Applying Pythagoras' Theorem
Using Pythagoras' theorem we can the following expressions for $$a^2, b^2, c^2$$ and $$d^2$$ in terms of $$x$$ and $$y$$.


 * $$ a^2 = x^2 + y^2 $$


 * $$ b^2 = (w - x)^2 + y^2 $$


 * $$ c^2 = (w-x)^2 + (h-y)^2 $$


 * $$ d^2 = x^2 + (h-y)^2 $$

Combining these equations in two pairs yields the following two equations.


 * $$ a^2 + c^2 = 2x^2 + 2y^2 - 2wx - 2hy + w^2 + h^2 $$


 * $$ b^2 + d^2 = 2x^2 + 2y^2 - 2wx - 2hy + w^2 + h^2 $$

It is clear from these that


 * $$ a^2 -b^2 + c^2 - d^2 = 0 $$

Propagation Delays
We can also express a, b, c and d in terms of the detection times $$t_a, t_b, t_c, t_d$$ and the unknowns $$t_0$$ (the time at which the tap occurred) and $$v$$ (the propagation speed).


 * $$ a = v(t_a - t_0) $$


 * $$ b = v(t_b - t_0) $$


 * $$ c = v(t_c - t_0) $$


 * $$ d = v(t_d - t_0) $$

As we did above with the corresponding pythagorean expressions, we now write expressions for $$a^2, b^2, c^2$$ and $$d^2$$.


 * $$ a^2 = v^2(t_a - t_0)^2 $$


 * $$ b^2 = v^2(t_b - t_0)^2 $$


 * $$ c^2 = v^2(t_c - t_0)^2 $$


 * $$ d^2 = v^2(t_d - t_0)^2 $$

Now, since we know that $$a^2 - b^2 + c^2 - d^2 = 0$$, we can write


 * $$v^2\left((t_a-t_0)^2 - (t_b-t_0)^2 + (t_c-t_0)^2 - (t_d-t_0)^2\right) = 0$$

We know that v, the propagation velocity cannot be zero. Therefore, we can say that


 * $$(t_a-t_0)^2 - (t_b-t_0)^2 + (t_c-t_0)^2 - (t_d-t_0)^2 = 0$$

If we multiply out each square term above, then rearrange the equation, it yields the following expression for $$t_0$$, the first of our unknowns.


 * $$t_0 = \frac{1}{2} \left( \frac{t_a^2 - t_b^2 + t_c^2 - t_d^2}{t_a-t_b+t_c-t_d} \right) $$

Since $$t_a, t_b, t_c$$ and $$t_d$$ are all known (i.e. they were detected directly from the accelerometer readings), we can now determine exactly how long before the first readings were detected the actual tap took place.

We define four new time variables, $$\tau_a, \tau_b, \tau_c$$ and $$\tau_d$$, which record the actual propagation time taken for the tap vibrations to travel from $$p$$ to each of the four corners.


 * $$\tau_a = t_a - t_0$$


 * $$\tau_b = t_b - t_0$$


 * $$\tau_c = t_c - t_0$$


 * $$\tau_d = t_d - t_0$$

Scale Replica Rectangle
We now construct a scale replica of the original rectangle using $$\tau_a, \tau_b, \tau_c, \tau_d$$ in place of the lengths $$a, b, c, d$$. The width and height of our scale replica ($$w_s$$ and $$h_s$$ respectively) are not known initially, but by specifying the lengths of the four lines joining $$p_s$$ to each of the corners $$A_s, B_s, C_s, D_s$$, they will be uniquely determined.

The four corners of our scale replica rectangle are


 * $$A_s$$ is the bottom left corner, at coordinates $$(0,0)$$
 * $$B_s$$ is the bottom right corner, at coordinates $$(w_s,0)$$
 * $$C_s$$ is the top right corner, at coordinates $$(w_s,h_s)$$
 * $$D_s$$ is the top left corner, at coordinates $$(0,h_s)$$



The relative proportions of $$\tau_a, \tau_b, \tau_c$$ and $$\tau_d$$ are exactly the same as those of $$a,b,c$$ and $$d$$. The point $$p_s = (x_s, y_s)$$ within our replica rectangle corresponds to $$p$$ within the original rectangle.

Applying Pythagoras' theorem again,


 * $$\tau_a^2 = x_s^2 + y_s^2$$

Therefore,


 * $$y_s = \sqrt{\tau_a^2 - x_s^2}$$

From the rectangle, we can also write these two equations:


 * $$w_s = x_s + \sqrt{\tau_b^2 - y_s^2} $$


 * $$h_s = y_s + \sqrt{\tau_d^2 - x_s^2} $$

Combining these two equations,


 * $$\frac{x_s}{w_s} - \frac{y_s}{h_s} + \frac{\sqrt{\tau_b^2 - y_s^2}}{w_s} - \frac{\sqrt{\tau_d^2 - x_s^2}}{h_s} = 0$$

Multiplying all terms by $$h_s$$ gives


 * $$\frac{h_s}{w_s}x_s - y_s + \frac{h_s}{w_s}\sqrt{\tau_b^2 - y_s^2} - \sqrt{\tau_d^2 - x_s^2} = 0$$

Since our replica rectangle has the same proportions as the original one, $$\frac{h_s}{w_s} = \frac{h}{w}$$, we can write


 * $$\frac{h}{w}x_s - y_s + \frac{h}{w}\sqrt{\tau_b^2 - y_s^2} - \sqrt{\tau_d^2 - x_s^2} = 0$$

Now, substitute in the previous expression for $$y_s$$ in terms of $$x_s$$.


 * $$\frac{h}{w}x_s - \sqrt{\tau_a^2 - x_s^2} + \frac{h}{w}\sqrt{\tau_b^2 - \tau_a^2 + x_s^2} - \sqrt{\tau_d^2 - x_s^2} = 0$$

We should now be able to solve for $$x_s$$, the only unknown in the previous equation. However, I've spent quite a few hours trying to hammer out an expression for $$x_s$$ in terms of $$\tau_a, \tau_b, \tau_d, w$$ and $$h$$ and I'm embarrassed to say that my efforts have been fruitless. Having been frustrated in my search for an elegant analytical solution, I've had to shelve my pride and settle for a numerical solution for the time being.

Numerical solution for $$x_s$$
Since we know that $$x_s$$ is real and positive (as are $$\tau_a, \tau_b, \tau_c, \tau_d$$), we can immediately identify some bounds within which it must lie.


 * $$x_s < \tau_a \,$$


 * $$x_s < \tau_d \,$$


 * $$x_s^2 > \tau_a^2 - \tau_b^2 \quad \mbox{for} \quad \tau_a > \tau_b$$


 * $$x_s^2 > 0 \quad \mbox{for} \quad \tau_a \leq \tau_b$$

The first three of these conditions arise from the fact that if they are violated, square root terms in our equation will become imaginary (and we know $$x_s$$ must be real).

The value of $$x_s$$ can be estimated with arbitrary precision using the following iterative algorithm. First define a function $$f(x)$$.


 * $$f(x_s) = \frac{h}{w}x_s - \sqrt{\tau_a^2 - x_s^2} + \frac{h}{w}\sqrt{\tau_b^2 - \tau_a^2 + x_s^2} - \sqrt{\tau_d^2 - x_s^2}$$

The desired value of x_s is that for which $$f(x_s) = 0$$. We define variables $$x_{min}$$ and $$x_{max}$$ and assign them the following initial values:


 * $$x_{min} = \left\{

\begin{array}{ll} \sqrt{\tau_a^2 - \tau_b^2} & ,\tau_a > \tau_b \\ \\ 0 &, \tau_a \leq \tau_b \end{array} \right.$$


 * $$x_{max} = \min{\left(\tau_a, \tau_d\right)} $$

I'm making the assumption here that $$f(x_s)$$ is monotonic on the interval $$x_{min} < x_s < x_{max}$$ (and hence that there is only one solution for $$f(x_s)=0$$ in that interval). Based on my visual inspection of the geometry in this problem, the solution is unique. However, it would probably be advisable to revisit this to confirm that the assumption is valid.

Now, perform the following steps as many times as necessary.


 * 1) Set $$x_s = \frac{x_{min}+x_{max}}{2}$$
 * 2) If $$f(x_s) \times f(x_{min}) < 0 $$ then set $$x_{min} = x_s \,$$. Otherwise, set $$x_{max} = x_s \,$$
 * 3) If $$x_{max} - x_{min} < \epsilon \,$$ for some predefined tolerance $$\epsilon \,$$ then cease iterating and take the current value of $$x_s \,$$ as the final value. Otherwise, repeat from step 1 above.

What is happening here is that $$x_{min}$$ and $$x_{max}$$ begin on opposite sides of the zero-crossing point of $$f(x_s)$$. Each iteration of the algorithm ratchets $$x_{min}$$ and $$x_{max}$$ closer to each other, but keeps the zero-crossing point in between them. At each iteration, the distance between $$x_{min}$$ and $$x_{max}$$ is reduced by half. After several iterations, there will only be a very short distance between $$x_{min}$$ and $$x_{max}$$ and the zero-crossing point is known to be somewhere in that very small range.

Finally, calculate $$p=(x,y)$$
Now, we have the value of $$x_s$$, the x-coordinate of the point we seek in the scale replica square. To find $$y_s$$, we simply use the following equation from above.


 * $$y_s = \sqrt{\tau_a^2 - x_s^2}$$

We also need to find out the width and height of the replica rectangle. For this, we can use the following equations from above.


 * $$w_s = x_s + \sqrt{\tau_b^2 - y_s^2} $$


 * $$h_s = y_s + \sqrt{\tau_d^2 - x_s^2} $$

Now that we know the values $$x_s, y_s$$ and $$w_s, h_s$$, we can easily find the actual point, $$p$$, at which the tap occurred in the original (unscaled) rectangle. Simply scale $$(x_s, y_s)$$ using the ratio of the dimensions of the replica rectangle to those of the original rectangle.


 * $$x = \frac{w}{w_s} x_s$$


 * $$y = \frac{w}{w_s} y_s$$

These are the coordinates of $$p$$, the point at which the tap occurred.

Summary of calculation
To recap, here is the sequence of calculations required to find $$p=(x,p)$$ when $$t_a, t_b, t_c, t_d, w$$ and $$h$$ are given.

First, find $$t_0$$:


 * $$t_0 = \frac{1}{2} \left( \frac{t_a^2 - t_b^2 + t_c^2 - t_d^2}{t_a-t_b+t_c-t_d} \right) $$

Next calculate $$\tau_a, \tau_b, \tau_c$$ and $$\tau_d$$:


 * $$\tau_a = t_a - t_0$$


 * $$\tau_b = t_b - t_0$$


 * $$\tau_c = t_c - t_0$$


 * $$\tau_d = t_d - t_0$$

Define the function $$f(x_s)$$:


 * $$f(x_s) = \frac{h}{w}x_s - \sqrt{\tau_a^2 - x_s^2} + \frac{h}{w}\sqrt{\tau_b^2 - \tau_a^2 + x_s^2} - \sqrt{\tau_d^2 - x_s^2}$$

Calculate initial values for $$x_{min}$$ and $$x_{max}$$:


 * $$x_{min} = \left\{

\begin{array}{ll} \sqrt{\tau_a^2 - \tau_b^2} & ,\tau_a > \tau_b \\ \\ 0 &, \tau_a \leq \tau_b \end{array} \right.$$


 * $$x_{max} = \min{\left(\tau_a, \tau_d\right)} $$

Perform the following steps as many times as necessary.


 * 1) Set $$x_s = \frac{x_{min}+x_{max}}{2}$$
 * 2) If $$f(x_s) \times f(x_{min}) < 0 $$ then set $$x_{min} = x_s \,$$. Otherwise, set $$x_{max} = x_s \,$$
 * 3) If $$x_{max} - x_{min} < \epsilon \,$$ for some predefined tolerance $$\epsilon \,$$ then cease iterating and take the current value of $$x_s \,$$ as the final value. Otherwise, repeat from step 1 above.

Calculate $$y_s$$:


 * $$y_s = \sqrt{\tau_a^2 - x_s^2}$$

Calculate $$w_s$$:


 * $$w_s = x_s + \sqrt{\tau_b^2 - y_s^2} $$

Finally, calculate $$p=(x,y)$$:


 * $$x = \frac{w}{w_s} x_s$$


 * $$y = \frac{w}{w_s} y_s$$