From 22810198535e1006395898ed23f2e11a311a04bd Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Tue, 14 Nov 2023 16:46:35 +0100 Subject: [PATCH] Pulse: as generating Thesis plots --- simulations/11_pulsed_timing.py | 347 ++++++++++++++++++++++++++------ 1 file changed, 289 insertions(+), 58 deletions(-) diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index 191ba46..8c2a6d9 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -1,4 +1,10 @@ #!/usr/bin/env python3 +# vim: fdm=marker fmr=<<<,>>> + +# TODO: compare with Peak Hilbert Envelope Timing + +# Remove non-cross Points in SNR plot +# extrapolate exponential to lower snr values from lib import util @@ -79,7 +85,7 @@ def antenna_bp(trace, low_bp, high_bp, dt, order=3): return bandpassed -def my_correlation(in1, template, lags=None): +def my_correlation(in1, template, lags=None, normalise=True): template_length = len(template) in1_length = len(in1) @@ -118,6 +124,9 @@ def my_correlation(in1, template, lags=None): corrs[i] = np.dot(in1_slice, template_slice) + if normalise: + corrs /= np.amax(corrs) + return corrs, (in1, template, lags) def trace_upsampler(trace, template_t, trace_t): @@ -149,7 +158,7 @@ def trace_downsampler(trace, template_t, trace_t, offset): pass -def hilbert_envelope_max_amplitude_time(trace, trace_t, do_plot=False, fname_distinguish=None): +def hilbert_envelope_max_amplitude_time(trace, trace_t, do_plot=False, fname_distinguish='', zoom_wx=50, inset_zoom_extent=(0.03, 0.4, 0.53, 0.57)): analytic_signal = signal.hilbert(trace) envelope = abs(analytic_signal) @@ -173,9 +182,43 @@ def hilbert_envelope_max_amplitude_time(trace, trace_t, do_plot=False, fname_dis if True: ax.legend() + ax.grid() fig.tight_layout() fig.savefig(f'figures/11_hilbert_timing{fname_distinguish}.pdf') + if zoom_wx: + xlims = ax.get_xlim() + + if not hasattr(zoom_wx, '__len__'): + zoom_wx = (zoom_wx, zoom_wx) + + if inset_zoom_extent: # do inset axes + orig_ax = ax + + axins = orig_ax.inset_axes(inset_zoom_extent) + axins.patch.set_alpha(0.9) + axins.set_yticklabels([]) + axins.set_xlim(t_max - zoom_wx[0], t_max + zoom_wx[-1]) + axins.grid() + + # replot data + axins.plot(trace_t, trace, label='Waveform') + axins.plot(trace_t, envelope, ls='dashed', label='Envelope') + + # indicate maximum and t_max + axins.axhline(envelope[max_idx], ls='dotted', color='g') + axins.axvline(t_max, ls='dotted', color='g') + + # increase margins and indicate inset zoom + orig_ax.margins(y=0.09) + orig_ax.indicate_inset_zoom(axins) + else: + ax.set_xlim(t_max - zoom_wx[0], t_max + zoom_wx[-1]) + + fig.tight_layout() + fig.savefig(f'figures/11_hilbert_timing{fname_distinguish}_zoom.pdf') + + ax.set_xlim(*xlims) plt.close(fig) @@ -239,7 +282,7 @@ def get_time_residuals_for_template( snr_sigma_factor=10,bp_freq=(0,np.inf), normalise_noise=False, h5_cache_fname=None, read_cache=True, write_cache=None, rng=rng, tqdm=tqdm, - peak_window=[0.2, 0.8], + peak_window=[0.6, 0.65], ): # Read in cached time residuals if read_cache: @@ -248,11 +291,14 @@ def get_time_residuals_for_template( else: cached_time_residuals, cached_snrs, cached_hilbert_time_residuals = np.array([]), np.array([]), np.array([]) + print(cached_hilbert_time_residuals.shape) + print(cached_time_residuals.shape) + # # Find difference between true and templated times # - hilbert_interp_t_max, _ = hilbert_envelope_max_amplitude_time(interp_template.signal, interp_template.t) + hilbert_interp_t_max, _ = hilbert_envelope_max_amplitude_time(interp_template.signal, interp_template.t, zoom_wx=None) time_residuals = np.zeros(max(0, (N_residuals - len(cached_time_residuals)))) snrs = np.zeros_like(time_residuals) @@ -305,7 +351,7 @@ def get_time_residuals_for_template( # Show signals if do_plots: - fig, axs = plt.subplots(2, sharex=True) + fig, axs = plt.subplots(1, sharex=True) if not hasattr(axs, '__len__'): axs = [axs] @@ -319,18 +365,16 @@ def get_time_residuals_for_template( if True: # indicate signal and noise levels level_kwargs = dict(ls='dashed', alpha=0.4) - axs[0].axhline(antenna.signal_level, color=l2[0].get_color(), **level_kwargs, label='Signal Level') - axs[0].axhline(antenna.noise_level, color=l3[0].get_color(), **level_kwargs, label='Noise Level') + axs[0].axhline(antenna.signal_level, color=l2[0].get_color(), **level_kwargs)#, label='Signal Level') + axs[0].axhline(antenna.noise_level, color=l3[0].get_color(), **level_kwargs)#, label='Noise Level') - axs[0].legend(title=f'SNR = {antenna.signal_to_noise:.1g}') + axs[0].legend(title=f'SNR = {antenna.signal_to_noise:.2g}', loc='lower right') axs[0].grid() if len(axs) > 1: - axs[1].set_title("Template") axs[1].set_ylabel("Amplitude") - axs[1].plot(template.t, template.signal, label='orig') - axs[1].plot(template.t + true_time_offset, template.signal, label='true moved orig') + axs[1].plot(template.t + true_time_offset, template.signal, label='Template') axs[1].legend() axs[1].grid() @@ -338,11 +382,53 @@ def get_time_residuals_for_template( fig.savefig(f'figures/11_antenna_signals_tdt{template.dt:.1g}.pdf') if True: # zoom - wx = 100 - x0 = true_time_offset + wx = 50 + x0 = true_time_offset + wx/2 old_xlims = axs[0].get_xlim() - axs[0].set_xlim( x0-wx, x0+wx) + if True: # do inset axes + extent = [0.03, 0.4, 0.53, 0.57] + orig_ax = axs[0] + + axins = orig_ax.inset_axes(extent) + axins.patch.set_alpha(0.9) + axins.set_yticklabels([]) + axins.set_xlim(x0-wx, x0+wx) + axins.grid() + + # replot data + l1 = axins.plot(antenna.t, antenna.signal, label='Filtered w/ noise', alpha=0.7) + l2 = axins.plot(antenna.t, antenna.signal - filtered_noise, label='Filtered w/o noise', alpha=0.7) + l3 = axins.plot(antenna.t, filtered_noise, label='Noise', alpha=0.7) + + if True: # indicate signal and noise levels + level_kwargs = dict(ls='dashed', alpha=0.4) + + axins.axhline(antenna.signal_level, color=l2[0].get_color(), **level_kwargs)#, label='Signal Level') + axins.axhline(antenna.noise_level, color=l3[0].get_color(), **level_kwargs)#, label='Noise Level') + + # increase margins and indicate inset zoom + orig_ax.margins(y=0.09) + orig_ax.indicate_inset_zoom(axins) + + if len(axs) > 1: + orig_ax = axs[1] + axins2 = orig_ax.inset_axes(extent) + axins2.patch.set_alpha(axins.patch.get_alpha()) + axins2.set_yticklabels([]) + axins2.set_xlim(x0-wx, x0+wx) + axins2.grid() + + # replot data + axins2.plot(template.t + true_time_offset, template.signal) + + # increase margins and indicate inset zoom + orig_ax.margins(y=0.1) + orig_ax.indicate_inset_zoom(axins2) + else: + axs[0].set_xlim( x0-wx, x0+wx) + + fig.tight_layout() fig.savefig(f'figures/11_antenna_signals_tdt{template.dt:.1g}_zoom.pdf') # restore @@ -361,15 +447,39 @@ def get_time_residuals_for_template( axs2[-1].set_xlabel("Time [ns]") axs2[0].set_ylabel("Amplitude") - axs2[0].plot(antenna.t, antenna.signal, marker='o', label='orig') + axs2[0].plot(antenna.t, antenna.signal, marker='o', label='waveform') axs2[0].plot(upsampled_t, upsampled_trace, label='upsampled') axs2[0].legend(loc='upper right') + axs2[0].grid() + + fig2.tight_layout() fig2.savefig(f'figures/11_upsampled_tdt{template.dt:.1g}.pdf') - wx = 1e2 - x0 = upsampled_t[0] + wx - 5 - axs2[0].set_xlim(x0-wx, x0+wx) + wx = 0.25e2 + x0 = upsampled_t[np.argmax(upsampled_trace)] - 5 + + if True: # do inset axes + extent = [0.03, 0.4, 0.47, 0.57] + orig_ax = axs2[0] + + axins = orig_ax.inset_axes(extent) + axins.patch.set_alpha(0.9) + axins.set_yticklabels([]) + axins.set_xlim(x0-wx, x0+wx) + axins.grid() + + # replot data + axins.plot(antenna.t, antenna.signal, marker='o') + axins.plot(upsampled_t, upsampled_trace) + + # increase margins and indicate inset zoom + orig_ax.margins(y=0.1) + orig_ax.indicate_inset_zoom(axins) + else: + axs2[0].set_xlim(x0-wx, x0+wx) + + fig2.tight_layout() fig2.savefig(f'figures/11_upsampled_tdt{template.dt:.1g}_zoom.pdf') if True: @@ -385,7 +495,7 @@ def get_time_residuals_for_template( best_time_lag = best_sample_lag * lag_dt # Find Hilbert Envelope t0 - hilbert_best_time_lag, _ = hilbert_envelope_max_amplitude_time(upsampled_trace, upsampled_t, do_plot=do_plots) + hilbert_best_time_lag, _ = hilbert_envelope_max_amplitude_time(upsampled_trace, upsampled_t, do_plot=do_plots, zoom_wx=(6,12)) else: # downsampled template raise NotImplementedError @@ -411,13 +521,14 @@ def get_time_residuals_for_template( template_amp_scaler = max(abs(template.signal)) / max(abs(antenna.signal)) # start the figure - fig, axs = plt.subplots(2, sharex=True) + fig, axs = plt.subplots(2, sharex=True, figsize=(9,6)) ylabel_kwargs = dict( #rotation=0, - ha='right', + #ha='right', va='center' ) axs[-1].set_xlabel("Time [ns]") + axs[-1].set_xlabel("Time [ns]") offset_list = [ [best_time_lag, dict(label=template.name, color='orange')], @@ -437,10 +548,12 @@ def get_time_residuals_for_template( l = axs[i].plot(offset + template.t, template_amp_scaler * template.signal, **this_kwargs) axs[i].legend() + axs[i].grid() # Correlation i=1 axs[i].set_ylabel("Correlation", **ylabel_kwargs) + axs[i].grid() axs[i].plot(lags * lag_dt, corrs) # Lines across both axes @@ -454,20 +567,83 @@ def get_time_residuals_for_template( axs[0].axvline(offset + len(template.signal) * (template.t[1] - template.t[0]), color=this_kwargs['color'], alpha=0.7) + # separated axes + for i, myax in enumerate(axs): + [ axes.set_visible(False) for axes in axs] + myax.set_visible(True) + + fig.tight_layout() + fig.savefig(f'figures/11_corrs_tdt{template.dt:.1g}_axes{i}.pdf') + + # re enable all axes + [ axes.set_visible(True) for axes in axs] + fig.tight_layout() + fig.savefig(f'figures/11_corrs_tdt{template.dt:.1g}.pdf') + if True: # zoom wx = len(template.signal) * (min(1,template.dt))/4 t0 = true_time_offset old_xlims = axs[0].get_xlim() - axs[i].set_xlim( x0-wx, x0+3*wx) + if True: # do inset axes + extent = [0.03, 0.4, 0.47, 0.57] + + axins = [] + for i in [0,1]: + orig_ax = axs[i] + axins.append(orig_ax.inset_axes(extent)) + axins[i].patch.set_alpha(0.9) + axins[i].set_yticklabels([]) + axins[i].set_xlim(x0-wx, x0+wx) + axins[i].grid() + + # replot data + if i == 0: + axins[i].plot(antenna.t, antenna.signal, label=antenna.name) + + # Plot the template + for offset_args in offset_list: + this_kwargs = offset_args[1] + offset = offset_args[0] + + l = axins[i].plot(offset + template.t, template_amp_scaler * template.signal, **this_kwargs) + + elif i == 1: # correlation + axins[i].plot(lags*lag_dt, corrs) + + # Lines across both axes + for offset_args in offset_list: + this_kwargs = offset_args[1] + offset = offset_args[0] + + for j in [0,1]: + axins[j].axvline(offset, ls='--', color=this_kwargs['color'], alpha=0.7) + + axins[i].axvline(offset + len(template.signal) * (template.t[1] - template.t[0]), color=this_kwargs['color'], alpha=0.7) + + # increase margins and indicate inset zoom + orig_ax.margins(y=0.1) + orig_ax.indicate_inset_zoom(axins[i]) + + else: + axs[i].set_xlim( t0-wx, t0+2*wx) + + # separated axes + for i, myax in enumerate(axs): + [ axes.set_visible(False) for axes in axs] + myax.set_visible(True) + + fig.tight_layout() + fig.savefig(f'figures/11_corrs_tdt{template.dt:.1g}_axes{i}_zoom.pdf') + + # re enable all axes + [ axes.set_visible(True) for axes in axs] + fig.tight_layout() fig.savefig(f'figures/11_corrs_tdt{template.dt:.1g}_zoom.pdf') # restore axs[i].set_xlim(*old_xlims) - fig.tight_layout() - fig.savefig(f'figures/11_corrs_tdt{template.dt:.1g}.pdf') - if True: plt.close(fig) @@ -496,27 +672,45 @@ if __name__ == "__main__": if os.name == 'posix' and "DISPLAY" not in os.environ: matplotlib.use('Agg') + figsize = (8,6) + fontsize = 12 + if True: + from matplotlib import rcParams + #rcParams["text.usetex"] = True + rcParams["font.family"] = "serif" + rcParams["font.size"] = fontsize + if not True:# small + figsize = (6, 4) + rcParams["font.size"] = "15" # 15 at 6,4 looks fine + elif True: # large + figsize = (9, 6) + rcParams["font.size"] = "16" # 15 at 9,6 looks fine + rcParams["grid.linestyle"] = 'dotted' + rcParams["figure.figsize"] = figsize + fontsize = rcParams['font.size'] - if False: - plt.rc('font', size=25) - - figsize = (12,12) bp_freq = (30e-3, 80e-3) # GHz interp_template_dt = 5e-5 # ns template_length = 200 # ns antenna_dt = 2 # ns - antenna_timelength = 1024 # ns + antenna_timelength = 2048 # ns N_residuals = 50*3 if len(sys.argv) < 2 else int(sys.argv[1]) - template_dts = np.array([antenna_dt, 5e-1, 5e-2]) # ns + if True: + template_dts = np.array([ 5e-1, 1e-1, 1e-2]) # ns + elif True: + template_dts = np.array([1e-2]) # ns + else: + template_dts = np.array([antenna_dt, 5e-1]) # ns 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] + #[5, 50] ), axis=None, dtype=float) @@ -536,16 +730,18 @@ if __name__ == "__main__": # to create an 'analog' sampled antenna interp_template, _deltapeak = create_template(dt=interp_template_dt, timelength=template_length, bp_freq=bp_freq, name='Interpolation Template', normalise=True) - interp_template.interpolate = interpolate.interp1d(interp_template.t, interp_template.signal, assume_sorted=True, fill_value=0, bounds_error=False, copy=False) + interp_template.interpolate = interpolate.interp1d(interp_template.t, interp_template.signal, assume_sorted=True, fill_value=0, bounds_error=False, copy=False)#, kind='nearest') - if True: # show interpolation template + if not True: # show interpolation template fig, ax = plt.subplots() ax.set_title("Filter Response") ax.set_xlabel("Time [ns]") ax.set_ylabel("Amplitude") ax.plot(interp_template.t, max(interp_template.signal)*_deltapeak[0], label='Impulse') ax.plot(interp_template.t, interp_template.signal, label='Filtered Signal') - ax.legend() + ax.legend(loc='upper right') + ax.grid() + fig.tight_layout() fig.savefig('figures/11_filter_response.pdf') if True: # show filtering equivalence samplerates @@ -556,7 +752,8 @@ if __name__ == "__main__": ax.plot(_time, max(_bandpassed)*_deltapeak[0], label='Impulse Antenna') ax.plot(_time, _bandpassed, label='Filtered Antenna') - ax.legend() + ax.legend(loc='upper right') + fig.tight_layout() fig.savefig('figures/11_interpolation_deltapeak+antenna.pdf') if True: @@ -615,18 +812,19 @@ if __name__ == "__main__": hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step') fig, ax = plt.subplots() - ax.set_title( - "Template Correlation Lag finding" - + f"\n template dt: {template_dt: .1e}ns" - + f"; antenna dt: {antenna_dt: .1e}ns" - + ";" if not mask_count else "\n" - + f"snr_factor: {snr_sigma_factor: .1e}" - + "" if not mask_count else f"; N_masked: {mask_count}" - ) + #ax.set_title( + # "Template Correlation Lag finding" + # + f"\n template dt: {template_dt: .1e}ns" + # + f"; antenna dt: {antenna_dt: .1e}ns" + # + ";" if not mask_count else "\n" + # + f"snr_factor: {snr_sigma_factor: .1e}" + # + "" if not mask_count else f"; N_masked: {mask_count}" + # ) ax.set_xlabel("Time Residual [ns]") ax.set_ylabel("#") + ax.grid() - if True: + if not True: # indicate boxcar accuracy limits for sign in [-1, 1]: ax.axvline( sign*template_dt/np.sqrt(12), ls='--', alpha=0.5, color='green') @@ -671,8 +869,17 @@ if __name__ == "__main__": chisq_strs ) - ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes) + ax.text( *(0.02, 0.95), text_str, ha='left', va='top', transform=ax.transAxes) + if True: + ax.legend(title=f"$\\langle SNR \\rangle$ = {snr_sigma_factor:.2g}", loc='upper right') + + if True: + this_lim = 55 + if ax.get_ylim()[1] <= this_lim: + ax.set_ylim([None, this_lim]) + + fig.tight_layout() if mask_count: fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{snr_sigma_factor:.1e}_masked.pdf") else: @@ -687,20 +894,27 @@ if __name__ == "__main__": # if True: enable_threshold_markers = [False, False, True, True] - threshold_markers = ['^', 'v', '8', 'X'] # make sure to have filled markers here + threshold_markers = ['^', 'v', '8', 'o'] # make sure to have filled markers here mask_thresholds = np.array([np.inf, N_residuals*0.5, N_residuals*0.1, 1, 0]) fig, ax = plt.subplots(figsize=figsize) ax.set_title(f"Template matching SNR vs time accuracy") - ax.set_xlabel("Signal to Noise Factor") + ax.set_xlabel("Signal to Noise") ax.set_ylabel("Time Accuracy [ns]") ax.grid() - ax.legend(title="\n".join([ + ax.legend(title=", ".join([ f"N={N_residuals}", #f"template_dt={template_dt:0.1e}ns", - f"antenna_dt={antenna_dt:0.1e}ns", - ])) + f"$1/f_s$ ={antenna_dt}ns", + ]), loc='lower left') + + if not True: + ax.set_title(f"Template matching, $N={N_residuals}$, $dt={antenna_dt}\\mathrm{{ns}}$") + + if False: + pass + # add wrong_peak_condition_multiple into plot # plot the values per template_dt slice template_dt_colors = [None]*len(template_dts) @@ -733,11 +947,11 @@ if __name__ == "__main__": snr_sigma_factor *= 2 # plot all invalid datapoints - if True: + if False: ax.plot(snrs[~valid_mask], y_values[~valid_mask], color='grey', **scatter_kwargs) # plot valid datapoints - if True: + if False: if template_dt_colors[a] is not None: scatter_kwargs['color'] = template_dt_colors[a] @@ -747,20 +961,26 @@ if __name__ == "__main__": masked_count = np.count_nonzero(~valid_mask) + threshold_index = np.argmin(masked_count <= mask_thresholds) -1 + + if not enable_threshold_markers[threshold_index]: + continue + # plot accuracy indicating masking counts kwargs = dict( ls='none', color= None if template_dt_colors[a] is None else template_dt_colors[a], - marker=threshold_markers[np.argmin( masked_count <= mask_thresholds)-1], + marker=threshold_markers[threshold_index], ms=10, markeredgecolor='white', markeredgewidth=1, + alpha=0.8 ) #l = ax.plot(snr_sigma_factor, np.sqrt(np.mean(y_values[valid_mask])**2), **{**kwargs, **dict(ms=50)}) if False: - l = ax.errorbar(snr_sigma_factor, time_accuracy, yerr=time_accuracy_std, xerr=snr_std, **kwargs) + l = ax.errorbar(snr_sigma_factor, time_accuracy, yerr=time_accuracy_std, xerr=snr_std, **kwargs, capsize=5) else: l = ax.plot(snr_sigma_factor, time_accuracy, **kwargs) @@ -771,14 +991,19 @@ if __name__ == "__main__": # indicate boxcar threshold if True: - ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color=template_dt_colors[a], label=f'Template dt:{template_dt:0.1e}ns') + ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color=template_dt_colors[a], label=f'{template_dt}ns') + text_coord = (0.03, template_dt/np.sqrt(12)) + ax.text( *text_coord, f'${template_dt}\mathrm{{\,ns}} / \sqrt{{12}}$', va='bottom', ha='left', color=template_dt_colors[a], fontsize=fontsize-1, transform=ax.get_yaxis_transform()) # Set horizontal line at 1 ns - if True: + if not True: ax.axhline(1, ls='--', alpha=0.8, color='g') - ax.legend() + if not True: + ax.legend(title="Template dt", loc='lower left') + elif True: + ax.legend().remove() fig.tight_layout() fig.savefig(f"figures/11_time_res_vs_snr_full_linear.pdf") @@ -793,6 +1018,11 @@ if __name__ == "__main__": this_lim = 1e1 if ax.get_ylim()[1] >= this_lim: ax.set_ylim([None, this_lim]) + # but keep it above 1 + if True: + this_lim = 1e0 + if ax.get_ylim()[1] <= this_lim: + ax.set_ylim([None, this_lim]) # require y-axis lower limit to be at least 1e-1 if True: @@ -807,12 +1037,13 @@ if __name__ == "__main__": low_ylims = ax.get_ylim()[0] if low_ylims <= this_lim: ax.set_ylim([this_lim, None]) - + # require x-axis lower limit to be under 1e0 if True: this_lim = 1e0 if ax.get_xlim()[0] >= this_lim: ax.set_xlim([this_lim, None]) + fig.tight_layout() if len(template_dts) == 1: