diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index 8aaaec2..98ba4b8 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -127,20 +127,20 @@ def trace_upsampler(trace, template_t, trace_t): def trace_downsampler(trace, template_t, trace_t, offset): pass -def read_time_residuals_cache(cache_fname, template_dt, antenna_dt, noise_sigma_factor, N=None): +def read_time_residuals_cache(cache_fname, template_dt, antenna_dt, snr_sigma_factor, N=None): try: with h5py.File(cache_fname, 'r') as fp: pgroup = fp['time_residuals'] pgroup2 = pgroup[f'{template_dt}_{antenna_dt}'] - ds_name = str(noise_sigma_factor) + ds_name = str(snr_sigma_factor) ds = pgroup2[ds_name] if N is None: return deepcopy(ds[:]) else: return deepcopy(ds[:min(N, len(ds))]) - except KeyError: + except (KeyError, FileNotFoundError): return np.array([]) def write_time_residuals_cache(cache_fname, time_residuals, template_dt, antenna_dt, noise_sigma_factor): @@ -164,13 +164,23 @@ if __name__ == "__main__": bp_freq = (30e-3, 80e-3) # GHz template_dt = 5e-2 # ns - template_length = 500 # ns + interp_template_dt = 5e-5 # ns + template_length = 200 # ns N_residuals = 50*3 if len(sys.argv) < 2 else int(sys.argv[1]) - noise_factors = [1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2, 1e-1, 3e-1, 5e-1, 7e-1] # amplitude factor + snr_factors = np.concatenate( # 1/noise_amplitude factor + ( + [0.25, 0.5, 0.75], + [1, 1.5, 2, 2.5, 3, 4, 5, 7], + [10, 20, 30, 50], + [100, 200, 300, 500] + ), + axis=None) antenna_dt = 2 # ns - antenna_timelength = 2048 # ns + antenna_timelength = 1024 # ns + + cut_wrong_peak_matches = True # # Create the template @@ -180,7 +190,16 @@ if __name__ == "__main__": template.signal = antenna_bp(_deltapeak[0], *bp_freq, template_dt) template.peak_sample = _deltapeak[1] template.peak_time = template.dt * template.peak_sample - interp1d_template = interpolate.interp1d(template.t, template.signal, assume_sorted=True, fill_value=0, bounds_error=False) + + # Interpolation Template + # to create an 'analog' sampled antenna + interp_template = Waveform(None, dt=interp_template_dt, name='Interpolation Template') + _interp_deltapeak = util.deltapeak(timelength=template_length, samplerate=1/interp_template.dt, offset=0) + interp_template.signal = antenna_bp(_interp_deltapeak[0], *bp_freq, interp_template.dt) + interp_template.peak_sample = _interp_deltapeak[1] + interp_template.peak_time = interp_template.dt * interp_template.peak_sample + interp_template.signal *= max(template.signal)/max(interp_template.signal) + interp_template.interpolate = interpolate.interp1d(interp_template.t, interp_template.signal, assume_sorted=True, fill_value=0, bounds_error=False) if True: # show template fig, ax = plt.subplots() @@ -189,6 +208,8 @@ if __name__ == "__main__": ax.set_ylabel("Amplitude") ax.plot(template.t, max(template.signal)*_deltapeak[0], label='Impulse Template') ax.plot(template.t, template.signal, label='Filtered Template') + ax.plot(interp_template.t, interp_template.signal, label='Filtered Interpolation Template', ls='dashed') + ax.legend() fig.savefig('figures/11_template_deltapeak.pdf') if True: # show filtering equivalence samplerates @@ -210,12 +231,13 @@ if __name__ == "__main__": # h5_cache_fname = f'11_pulsed_timing.hdf5' - time_accuracies = np.zeros(len(noise_factors)) - for k, noise_sigma_factor in tqdm(enumerate(noise_factors)): - print() #separating tqdm - + time_accuracies = np.zeros(len(snr_factors)) + for k, snr_sigma_factor in tqdm(enumerate(snr_factors)): # Read in cached time residuals - cached_time_residuals = read_time_residuals_cache(h5_cache_fname, template.dt, antenna_dt, noise_sigma_factor) + if True: + cached_time_residuals = read_time_residuals_cache(h5_cache_fname, template.dt, antenna_dt, snr_sigma_factor) + else: + cached_time_residuals = np.array([]) # # Find difference between true and templated times @@ -243,7 +265,8 @@ if __name__ == "__main__": antenna.t = util.sampled_time(1/antenna.dt, start=0, end=antenna_timelength) - antenna.signal = interp1d_template(antenna.t - antenna.peak_time) + # Sample the interpolation template + antenna.signal = interp_template.interpolate(antenna.t - antenna.peak_time) antenna.peak_sample = antenna.peak_time/antenna.dt antenna_true_signal = antenna.signal @@ -254,7 +277,7 @@ if __name__ == "__main__": antenna.signal *= -1 ## Add noise - noise_amplitude = max(template.signal) * noise_sigma_factor + noise_amplitude = max(template.signal) * 1/snr_sigma_factor noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal)) filtered_noise = antenna_bp(noise_realisation, *bp_freq, antenna.dt) @@ -323,16 +346,18 @@ if __name__ == "__main__": lag_dt = upsampled_t[1] - upsampled_t[0] corrs, (out1_signal, out2_template, lags) = my_correlation(upsampled_trace, template.signal) + # Determine best correlation time + idx = np.argmax(abs(corrs)) + best_sample_lag = lags[idx] + best_time_lag = best_sample_lag * lag_dt + else: # downsampled template raise NotImplementedError - corrs, (out1_signal, out2_template, lags) = my_downsampling_correlation(template.signal, antenna.signal, template.t, antenna.t) + corrs, (_, _, lags) = my_downsampling_correlation(antenna.signal, antenna.t, template.signal, template.t) lag_dt = upsampled_t[1] - upsampled_t[0] - # Determine best correlation time - idx = np.argmax(abs(corrs)) - best_sample_lag = lags[idx] - best_time_lag = best_sample_lag * lag_dt + # Calculate the time residual time_residuals[j] = best_time_lag - true_time_offset if not do_plots: @@ -408,13 +433,16 @@ if __name__ == "__main__": if True: plt.close(fig) + print()# separating tqdm print()# separating tqdm # Were new time residuals calculated? # Add them to the cache file if len(time_residuals) > 1: # merge cached and calculated time residuals time_residuals = np.concatenate((cached_time_residuals, time_residuals), axis=None) - write_time_residuals_cache(h5_cache_fname, time_residuals, template_dt, antenna_dt, noise_sigma_factor) + + if True: # write the cache + write_time_residuals_cache(h5_cache_fname, time_residuals, template_dt, antenna_dt, snr_sigma_factor) else: time_residuals = cached_time_residuals @@ -498,13 +526,17 @@ if __name__ == "__main__": ax.set_yscale('log') # plot the values - ax.plot(1/np.asarray(noise_factors), time_accuracies, ls='none', marker='o') + ax.plot(np.asarray(snr_factors), time_accuracies, ls='none', marker='o') + + if True: # limit y-axis to 1e0 + ax.set_ylim([None, 1e1]) # Set horizontal line at 1 ns if True: ax.axhline(1, ls='--', alpha=0.8, color='g') ax.grid() + ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color='b') fig.tight_layout() fig.savefig(f"figures/11_time_res_vs_snr_tdt{template_dt:0.1e}.pdf")