Impulse Response Calculation

[vinz-uts]

Hi, I'm not sure that the following is a bug or it is correct, I report what I notice:
I was looking the arrivals_to_impulse_response function in uwapm.py file and its computation. I notice that the amplitude of the arrival, row.arrival_amplitude, is setted in the right ndx position for each arrival, but, because ndx is computed as an approximation of the time of arrival (ndx = int(_np.round((row.time_of_arrival.real-t0)*fs))) can be happen that 2 different arrival have the same ndx index. In this case is correct to replace the ir[ndx] value with the last arrival amplitude, or should be added the new amplitude to the previous (ir[ndx] += row.arrival_amplitude)? In the first case aren't we losing one (or more) arrival?
I tried to generate arrivals with the following script:

import arlpy.uwapm as pm

fs = 10000  
env = pm.create_env2d()
rays = pm.compute_eigenrays(env, debug=True)
arrivals = pm.compute_arrivals(env)
ir = pm.arrivals_to_impulse_response(arrivals, fs, True)

Then I try to do some elaboration to investigate the results, I build two dictionary, A_arr, A_ir, which contains the pairs: time of arrival (approximated as in arrivals_to_impulse_response function) and a list of the amplitudes of each arrival received at this time. In lost_arr list is contained the tuple of "lost arrivals".

A_arr = {} # Amplitude of arrivals
A_ir  = {} # Amplitude of impulse response

for i in range(len(ir)):
    if ir[i] != 0:
        ir_T = (i/fs)
        A_ir[ir_T] = ir[i]

for _, arr in arrivals.iterrows():
    T_ = round(arr.time_of_arrival, 4)
    if T_ in A_arr.keys():
        A_arr[T_.real].append(arr.arrival_amplitude)
    else:
        A_arr[T_.real] = [arr.arrival_amplitude]

lost_arr = []
for (T, A) in A_arr.items():
    for A_ in A:
        if A_ != A_ir[T]:
            lost_arr.append((T, A_))

print(f'Lost arrivals: {lost_arr}')

I also try to plot the sum of the amplitude of the impulse response and the sum of the amplitudes of the arrival (which I expected to be the same), but they differ in the quantities present inlost_arr.

>>> print(f'Sum of arrival amplitudes: {sum(list(arrivals.arrival_amplitude))}')
>>> print(f'Sum of impulse response amplitudes: {sum(ir)}')
>>> print(f'Sum of impulse response amplitudes + lost arrivals: {sum(ir) + sum([e[1] for e in lost_arr])}')
Sum of arrival amplitudes: (0.00023585833449997007+0.001292865475452424j)
Sum of impulse response amplitudes: (-7.266767664939084e-05+0.0012367161871826558j)
Sum of impulse response amplitudes + lost arrivals: (0.00023585833449997015+0.0012928654754524244j)

I'm not sure if there is nothing wrong in my reasoning, if so I apologize, otherwise can you check too? In case it is correct to fix it replacing with ir[ndx] += row.arrival_amplitude the similar line in uwapm.arrivals_to_impulse_response function?

Thanks!

[mchitre]

Yes, it would make sense to add the amplitudes (as complex values). The current implementation implicitly assumes that the sampling rate is high enough to avoid overlapping arrivals. But even in that case, geometry can contrive to have two arrivals overlap, and we should handle it properly.

PR welcome, if you feel like taking a stab at fixing it @vinz-uts