Coffee Space

Two Microphone Single Observation

Imagine that you have two microphones at the left and right positions of a centre point of reference, $\{m_l, m_r\}$, separated by a distance $d$. What we want to know is the time difference $\Delta t$ for which the source sound arrives at each microphone, where $\Delta t = t_l - t_r$ (difference of left and right microphones). Wikipedia says that speed of sound, $V$, is $343 m/s$. A function for this would look as follows:

0001 import math
0002 
0003 # Constants
0004 V = 343 # m/s
0005 
0006 def calc_delta_t(d = 0.2, src = (0, 0)) :
0007   m_l_xy = { "x": -d / 2.0, "y": 0, "r": 0 }
0008   m_r_xy = { "x":  d / 2.0, "y": 0, "r": 0 }
0009   for m in [m_l_xy, m_r_xy] :
0010     a = src[1] - m["y"]
0011     b = src[0] - m["x"]
0012     m["r"] = math.sqrt((a ** 2) + (b ** 2))
0013   return (m_l_xy["r"] - m_r_xy["r"]) / V

Now for an environment of width and height $(w, h)$ in meters, with a stepping resolution $s$ in meters, we simulate a source sound for different locations about the centre $(x = 0, y = 0)$:

0014 def map_grid(d = 0.2, width = 10, height = 10, w_step = 0.5, h_step = 0.5) :
0015   grid = []
0016   y = -height / 2.0
0017   while y <= height / 2.0 :
0018     grid += [ [] ]
0019     x = -width / 2.0
0020     while x <= width / 2.0 :
0021       grid[-1] += [ [ calc_delta_t(d = d, src = (x, y)) ] ]
0022       x += w_step
0023     y += h_step
0024   return grid

Next we want to render this:

0025 import sys
0026 
0027 def display_grid(grid = [ [] ], fig_w = 800, fig_h = 800) :
0028   # Analyse
0029   g_max = -sys.float_info.max
0030   g_min = sys.float_info.max
0031   g_dmax = -sys.float_info.max
0032   for y in range(len(grid)) :
0033     for x in range(len(grid[y])) :
0034       if len(grid[y][x]) >= 1 :
0035         if grid[y][x][0] > g_max : g_max = grid[y][x][0]
0036         if grid[y][x][0] < g_min : g_min = grid[y][x][0]
0037       if len(grid[y][x]) >= 2 :
0038         d = abs(grid[y][x][0] - grid[y][x][1])
0039         if d > g_dmax : g_dmax = d
0040   assert(g_max != 0.0)
0041   assert(g_min != 0.0)
0042   assert(g_dmax!= 0.0)
0043   # Display diagram
0044   g = Svgg(width = fig_w, height = fig_h, bg = "#FFF")
0045   for y in range(len(grid)) :
0046     for x in range(len(grid[y])) :
0047       p_w = fig_w / len(grid[y])
0048       p_h = fig_h / len(grid)
0049       _r = 255
0050       _g = 255
0051       _b = 255
0052       for z in range(len(grid[y][x])) :
0053         if z == 0 :
0054           if grid[y][x][z] < 0 :
0055             _b = (255.0 / g_min) * (g_min - grid[y][x][z])
0056           elif grid[y][x][z] > 0 :
0057             _g = (255.0 / g_max) * (g_max - grid[y][x][z])
0058         if z == 1 :
0059           d = abs(grid[y][x][z - 1] - grid[y][x][z])
0060           _r = (255.0 / g_dmax) * (g_dmax - d)
0061       bg = "#" + hex(int(_r) + 256)[-2:] + hex(int(_b) + 256)[-2:] + hex(int(_g) + 256)[-2:]
0062       g.add("rect", {"x": x * p_w, "y": y * p_h, "width": p_w, "height": p_h, "fill": bg, "stroke": "#000"})
0063       for z in range(len(grid[y][x])) :
0064         data = str(round(grid[y][x][z] * 1000, 4))
0065         if z == 1 : data = str(round(grid[y][x][z - 1] - grid[y][x][z] * 1000, 4))
0066         if z >= 2 : data = str(grid[y][x][z])
0067         g.add("text", {"x": (x + 0.5) * p_w, "y": (y + ((1.0 / (len(grid[y][x]) + 1)) * (z + 1))) * p_h, "dominant-baseline": "middle", "text-anchor": "middle", "fill": "#000", "value": data})
0068   print(g.to_uri())
0069   return
0070 
0071 grid = map_grid(d = 0.2, width = 5, height = 5)
0072 display_grid(grid)

The arrival time delta is given in milliseconds, with each number representing the received difference if the audio source originated in that location. Negative numbers show that the sound was first detected in the left microphone $m_l$, and positive numbers show a bias for the right microphone $m_r$.

Next we want to plot the frequency in which each arrival time is observed:

0073 def sort_key(e) :
0074   return float(e[0])
0075 
0076 def plot_grid(grid, labels) :
0077   p = Plott(
0078     title = "Dual Mic Receive Delta Distribution for " + str(len(grid[0])) + "x" + str(len(grid)) + " Grid",
0079     x_title = "Time (ms)",
0080     y_title = "Frequency",
0081   )
0082   for z in range(len(grid[0][0])) :
0083     dist = {}
0084     for j in range(len(grid)) :
0085       for i in range(len(grid[j])) :
0086         v = grid[j][i][z]
0087         if z >= 1 : v = grid[j][i][0] - grid[j][i][z]
0088         v = str(round(v * 1000.0, 4))
0089         if not v in dist : dist[v] = 0
0090         dist[v] += 1
0091     dist = [(k, v) for k, v in dist.items()]
0092     dist.sort(key = sort_key, reverse = False)
0093     x = []
0094     y = []
0095     for d in range(len(dist)) :
0096       x += [ float(dist[d][0]) ]
0097       y += [ int(dist[d][1]) ]
0098     p.plot(x, y, labels[z])
0099   print(p.to_uri())
0100   return
0101 
0102 plot_grid(grid, ["Recieve Delta"])

As you can see, every observation has symmetry ($y \ge 2$). For an observing agent, you would not be able to differentiate between a source arriving in front or behind.

Turn Improvement

To break the symmetry, we have the observing agent turn, such that their microphones are rotated about $(0, 0)$. This is functionally equivalent of rotating the source of the sound:

0103 def rotate_source(src = (0, 0), angle = math.pi / 2) :
0104   x = src[0] * math.cos(angle) - src[1] * math.sin(angle)
0105   y = src[0] * math.sin(angle) + src[1] * math.cos(angle)
0106   return (x, y)

Now we compute both the initial observation delta $\Delta t_0 = t_l - t_r$, rotate the observing agent by an angle $\alpha$, and note the new observation delta $\Delta t_{\alpha} = t_l - t_r$:

0107 def map_delta(d = 0.2, a = 0, width = 10, height = 10, w_step = 0.5, h_step = 0.5) :
0108   grid = []
0109   y = -height / 2.0
0110   while y <= height / 2.0 :
0111     grid += [ [] ]
0112     x = -width / 2.0
0113     while x <= width / 2.0 :
0114       grid[-1] += [ [
0115         calc_delta_t(d = d, src = (x, y)),
0116         calc_delta_t(d = d, src = rotate_source((x, y), a)),
0117       ] ]
0118       x += w_step
0119     y += h_step
0120   return grid
0121 
0122 grid = map_delta(d = 0.2, a = math.radians(10), width = 5, height = 5)
0123 display_grid(grid)

The top number in each cell is as before, but the second value in each cell is $\Delta \Delta \mathbf{t} = \Delta t_0 - \Delta t_{\alpha}$, which is the difference in delta observations given a turn angle of $\alpha$.

Observe that the symmetry is now broken, with top left observations $(-,-)$, top right $(+,-)$, bottom left $(-,+)$ and bottom right $(+,+)$.

Now we plot this again:

0124 plot_grid(grid, ["Recieve Delta", "Recieve Turn α=10° Difference Delta"])

What’s not great is that for the chosen turn angle, $\Delta \Delta \mathbf{t}$ is a much smaller time. We should see what the optimal $\alpha$ is…

Optimal Turn Angle

Lastly, is there any preference for a given second observation angle? We assume that $0 \ge \alpha \ge 360$, and for each $alpha$ we look at the range $r$ (difference in time delta, $\Delta \Delta \mathbf{t}$), standard deviation of the data, the cost $c$ of turning to an angle with a linear time assumption (no acceleration), and finally the cost-benefit as $c \times r$:

0125 import statistics
0126 
0127 def grid_analyse(grid, z = 1) :
0128   vals = []
0129   g_max = -sys.float_info.max
0130   g_min = sys.float_info.max
0131   for y in range(len(grid)) :
0132     for x in range(len(grid[y])) :
0133       v = grid[y][x][z]
0134       if z >= 1 : v = grid[y][x][0] - grid[y][x][z]
0135       vals += [ v ]
0136       if v > g_max : g_max = v
0137       if v < g_min : g_min = v
0138   stddev = statistics.stdev(vals)
0139   return (g_min, g_max, stddev)
0140 
0141 x = []
0142 y = []
0143 z = []
0144 c = []
0145 cb = []
0146 for a in range(0, 360 - 1) :
0147   a += 1
0148   grid = map_delta(d = 0.2, a = math.radians(a), width = 5, height = 5)
0149   g_min, g_max, stddev = grid_analyse(grid)
0150   x += [ a ]
0151   y += [ round(abs(g_min - g_max) * 1000.0, 4) ]
0152   z += [ stddev ]
0153   c += [ y[-1] / a ]
0154   cb += [ y[-1] * c[-1] ]
0155   if y[-1] > 0.0 : y[-1] = math.log(y[-1])
0156   if z[-1] > 0.0 : z[-1] = math.log(z[-1])
0157   if c[-1] > 0.0 : c[-1] = math.log(c[-1])
0158   if cb[-1] > 0.0 : cb[-1] = math.log(cb[-1])
0159 
0160 p = Plott(
0161   title = "Dual Mic Receive Angle α Optimisation",
0162   x_title = "Angle α (degrees °)",
0163   y_title = "Min-Max Delta Range (ms)",
0164 )
0165 p.plot(x, y, "Range of time delta")
0166 p.plot(x, z, "Standard deviation of time delta")
0167 p.plot(x, c, "Cost with linear assumption")
0168 p.plot(x, cb, "Cost-benefit with linear assumption")
0169 print(p.to_uri())
0170 
0171 best_a = x[cb.index(max(cb))]
0172 print("Best turn angle α = " + str(best_a) + "°")

Best turn angle α = 135.0°

The more the observer turns between observations, the greater the delta observation, up until 180 degrees where the observation is mirrored. The standard deviation of the points also increases in a similar fashion. An approximate linear time cost for turning for the second observation shows that there is a bias to turning less. Cost benefit indicates that there is an optimum turn angle of $\alpha = 135$, but realistically the observer should turn as near to this angle as they can.

Let’s see the best second observation angle re-plotted:

0173 grid = map_delta(d = 0.2, a = math.radians(best_a), width = 5, height = 5)
0174 plot_grid(grid, ["Recieve Delta", "Recieve Turn α=" + str(best_a) + "° Difference Delta"])

We appear to have a nice spread, with only a few observations less 0.1ms. Next we display the grid again:

0175 display_grid(grid)

Notice how the symmetry is broken! There appears to be an observation bias! We can show that by rotating the other way:

0176 grid = map_delta(d = 0.2, a = math.radians(-best_a), width = 5, height = 5)
0177 display_grid(grid)

Turn Angle Optimisation

It appears as though the direction, and potentially the optimal turn angle, may be better optimised using the initial observation! Let’s find it!

The question we want to ask is slightly subtle: Given an initial sound event observation, and knowing the real source angle $\theta$, and knowing that a difficult to differentiate signal would come from $180 - \theta$, which is the best angle $\alpha$ to turn the observing agent?

0178 def map_delta_opt(d = 0.2, max_a = 180, width = 10, height = 10, w_step = 0.5, h_step = 0.5) :
0179   grid = []
0180   y = -height / 2.0
0181   while y <= height / 2.0 :
0182     grid += [ [] ]
0183     x = -width / 2.0
0184     while x <= width / 2.0 :
0185       best_a = 0
0186       best_o = -sys.float_info.max
0187       for a in range(-max_a, max_a + 1, 1) :
0188         if a == 0 : continue # Avoid divide by zero case
0189         _o = calc_delta_t(d = d, src = (x, y))
0190         o_ = calc_delta_t(d = d, src = rotate_source((x, y), math.radians(a)))
0191         o__ = calc_delta_t(d = d, src = rotate_source((x, y), math.radians(-a)))
0192         benefit = abs(o_ - o__)
0193         if benefit > best_o :
0194           best_a = a
0195           best_o = benefit
0196       if best_o <= 0.0 : best_a = 0 # Don't turn if nothing gained
0197       grid[-1] += [ [
0198         calc_delta_t(d = d, src = (x, y)),
0199         calc_delta_t(d = d, src = rotate_source((x, y), math.radians(best_a))),
0200         best_a,
0201       ] ]
0202       x += w_step
0203     y += h_step
0204   return grid
0205 
0206 grid = map_delta_opt(d = 0.2, max_a = 180, width = 5, height = 5)
0207 display_grid(grid)

Interesting! It appears the optimal angle for the observing agent to turn is 90 degrees, with no preference for direction!