diff --git a/gsrp_smart_util.py b/gsrp_smart_util.py
index ba9c0a63900f8ee7cb60636472beae5c95d2384a..a5463296616224a3cc9486d43e9b3c0200673f2b 100644
--- a/gsrp_smart_util.py
+++ b/gsrp_smart_util.py
@@ -90,9 +90,11 @@ def mul(mem1, mem2, cc, t_max, id1, id2, n_ind):
# assume len(id2) == 1
idx1, idx2 = np.meshgrid(np.arange(len(mem1[0])), np.arange(len(mem2[0])))
idx1, idx2 = idx1.ravel(), idx2.ravel()
- out_tij = np.concatenate((mem1[1][:, idx1], mem2[1][:, idx2]), axis=0)
- mask = (np.abs(out_tij[:-1] - out_tij[-1:]) <= t_max).all(0)
- out_tij = out_tij[:, mask]
+ out_tij = np.empty((len(id1) + 1, len(idx1)), mem1[1].dtype)
+ np.take(mem1[1], idx1, axis=1, out=out_tij[:-1])
+ np.take(mem2[1], idx2, axis=1, out=out_tij[-1:])
+ mask = reduce_all(np.abs(out_tij[:-1] - out_tij[-1:]) <= t_max) # Faster than numpy for large array
+ out_tij = np.compress(mask, out_tij, axis=1)
idx1, idx2 = idx1[mask], idx2[mask]
out_val = mem1[0][idx1] * mem2[0][idx2]
tij_dep = out_tij[:-1] - out_tij[-1:]
@@ -104,12 +106,19 @@ def mul(mem1, mem2, cc, t_max, id1, id2, n_ind):
def mask_val(mem, val):
mask = mem[0] >= val
- return mem[0][mask], mem[1][:, mask]
+ return mem[0][mask], np.compress(mask, mem[1], axis=1)
+
+
+def reduce_all(arr):
+ out = arr[0].copy()
+ for a in arr[1:]:
+ out &= a
+ return out
def mask_lim(mem, tij_min, tij_max, t_max):
- mask = ((tij_max - t_max <= mem[1]) & (mem[1] <= tij_min + t_max)).all(0)
- return mem[0][mask], mem[1][:, mask]
+ mask = reduce_all((tij_max - t_max <= mem[1]) & (mem[1] <= tij_min + t_max)) # Faster than numpy for large array
+ return mem[0][mask], np.compress(mask, mem[1], axis=1)
def constrained_argmax(mem, cc, tij_ind, curr_tij, used_tij, t_max, n_ind):
@@ -134,7 +143,7 @@ def smart_gsrp(cc, n_ind, n_tot, t_max, tree, program, clean_list):
if op == 'mem':
if i == 0:
memory[(0, j)] = mask_val((np.concatenate((cc[op.right, -t_max:], cc[op.right, :t_max + 1])),
- np.arange(-t_max, t_max + 1)[np.newaxis]), val)
+ np.arange(-t_max, t_max + 1, dtype=np.int32)[np.newaxis]), val)
else:
memory[(i, j)] = mask_val(memory[(op.left, op.right)], val)
else: # op == 'mul'
@@ -153,7 +162,7 @@ def smart_gsrp(cc, n_ind, n_tot, t_max, tree, program, clean_list):
except ValueError: # search of potential maxima ended outside of possible values
tij_min, tij_max = memory[(i, j)][1].min(), memory[(i, j)][1].max()
for k in range(j):
- memory[(i, k)] = mask_lim(memory[(i, k)], tij_min, tij_max, t_max)
+ memory[(i, k)] = mask_lim(memory[(i, k)], tij_min, tij_max, t_max)
# print('tdoa:', tij, 'val:', val, 'mem size:', (lambda x: f'{x} ({100 * x / (n_ind//(i+1)) / (2 * t_max + 1) ** (i+1)}%)')(sum(len(o[0].T) for o in memory.values())))
@@ -161,10 +170,4 @@ def smart_gsrp(cc, n_ind, n_tot, t_max, tree, program, clean_list):
for p in clean_list[i]:
del memory[p]
- return val, tij
-
-
-
-
-
-
+ return np.log10(val), tij
diff --git a/gsrp_tdoa_hyperres.py b/gsrp_tdoa_hyperres.py
index 334868d830f540ed06fc7e2c0d4fa16fde2ec923..ff59caf669b567d50091b03ab69acd340e093555 100755
--- a/gsrp_tdoa_hyperres.py
+++ b/gsrp_tdoa_hyperres.py
@@ -78,10 +78,10 @@ def corr(data, pos, w_size, max_tdoa, decimate=1, mode='prepare', hyper=True):
else:
raise ValueError(f'Unknown mode {mode}')
- tdoas = np.zeros((len(pos), num_channel_pairs + 1), np.float32)
+ tdoas = np.zeros((len(pos), num_channel_pairs + 2), np.float32)
if hyper: # prepare hyper res
- tdoas2 = np.zeros((len(pos), num_channel_pairs + 1), np.float32)
+ tdoas2 = np.zeros((len(pos), num_channel_pairs + 2), np.float32)
poly = PolynomialFeatures(2)
lin = LinearRegression()
pipe = Pipeline([('poly', poly), ('lin', lin)])
@@ -98,8 +98,8 @@ def corr(data, pos, w_size, max_tdoa, decimate=1, mode='prepare', hyper=True):
).named_steps['lin'].coef_
der = np.zeros((num_channels - 1, num_channels - 1))
der[ind] = coef[4:]
- poly_min = np.linalg.lstsq(der + der.T, -coef[1:4], rcond=None)[0] + mean
- return pipe.predict(poly_min[np.newaxis]).item(), mat @ poly_min
+ poly_min = np.linalg.lstsq(der + der.T, -coef[1:4], rcond=None)[0]
+ return np.log10(pipe.predict(poly_min[np.newaxis]).item()), mat @ (poly_min + mean)
cc = np.empty((num_channel_pairs, dw_size), np.float32)
for i in trange(len(pos)):
@@ -109,25 +109,27 @@ def corr(data, pos, w_size, max_tdoa, decimate=1, mode='prepare', hyper=True):
fft = np.asarray(fft, dtype=np.complex64)
cc[:] = irfft(fft[:, v2] * np.conj(fft[:, v1]), axis=0).T
cc -= cc.min(-1, keepdims=True)
+ maxs = cc.max(1, keepdims=True)
+ cc /= maxs
+ maxs = np.log10(maxs.prod())
if mode == 'prepare':
- tdoas[i, 0], index = c_corr.c_corr_at(cc, tausf)
- tdoas[i, 1:] = ((tausf[:, index] + dw_size // 2) % dw_size) - dw_size // 2
+ tdoas[i, :2], index = c_corr.c_corr_at(cc, tausf)
+ tdoas[i, 2:] = ((tausf[:, index] + dw_size // 2) % dw_size) - dw_size // 2
elif mode == 'on-the-fly':
- tdoas[i, 0], tdoas[i, 1:] = c_corr.c_corr_all(cc, cc_size//2, num_channels - 1)
+ tdoas[i, :2], tdoas[i, 2:] = c_corr.c_corr_all(cc, cc_size//2, num_channels - 1)
elif mode == 'smart':
- maxs = cc.max(1, keepdims=True)
- cc /= maxs
- val, tdoas[i, 1:] = smart_gsrp(cc, num_channels - 1, num_channel_pairs, cc_size // 2,
+ tdoas[i, :2], tdoas[i, 2:] = smart_gsrp(cc, num_channels - 1, num_channel_pairs, cc_size // 2,
tree, program, clean_list)
- tdoas[i, 0] = np.log10(val * maxs.prod())
- cc *= maxs
else:
raise ValueError(f'Unknown mode {mode}')
+ tdoas[i, 1] += maxs
if hyper:
- tdoas2[i, 0], tdoas2[i, 1:] = _hyperres(tdoas[i, 1:], cc)
-
+ tdoas2[i, :2], tdoas2[i, 2:] = _hyperres(tdoas[i, 2:], cc)
+ tdoas2[i, 1] += maxs
+ tdoas[:, :2] *= 20
if hyper:
+ tdoas2[:, :2] *= 20
return np.hstack((np.expand_dims(pos, -1), tdoas)), np.hstack((np.expand_dims(pos, -1), tdoas2))
else:
return np.hstack((np.expand_dims(pos, -1), tdoas))
@@ -219,8 +221,9 @@ if __name__ == "__main__":
parser.add_argument('infile', type=str, help='The sound file to process.')
parser.add_argument('outfile', type=str, help='The text or npy file to write results to. Each row gives the '
- 'position (in samples), cross-correlation product, the independent '
- 'TDOAs (in samples), and TDOAs derived from the independent ones. ')
+ 'position (in samples), cross-correlation product in decibel '
+ '(normalized and unormalized), the independent TDOAs (in samples), '
+ 'and TDOAs derived from the independent ones.')
group1 = parser.add_argument_group('Channels')
group1.add_argument('-c', '--channels', type=intlist, default=None,
@@ -257,7 +260,7 @@ if __name__ == "__main__":
group4.add_argument('-e', '--erase', action='store_false', help='Erase existing outfile. If outfile exist and '
'--erase is not provide, the script will exit.')
group4.add_argument('-n', '--no-hyperres', action='store_true', help='Disable the hyper resolution evalutation of '
- 'the TDOA')
+ 'the TDOA')
group4.add_argument('-M', '--mode', choices={'prepare', 'on-the-fly', 'smart', 'auto'}, default='smart',
help=f'R|How to explore the TDOA space (default: %(default)s).\n'
f'{BColors.BOLD}prepare{BColors.ENDC} precomputes all the possible TDOA pairs and then '