ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-12-21 19:13:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

Pulse: as generating Thesis plots

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-11-14 16:46:35 +01:00
parent b9314a2800
commit 2281019853
1 changed files with 289 additions and 58 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

341
simulations/11_pulsed_timing.py
View file

@ -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()
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
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:
@ -814,6 +1044,7 @@ if __name__ == "__main__":
if ax.get_xlim()[0] >= this_lim:
ax.set_xlim([this_lim, None])
fig.tight_layout()
if len(template_dts) == 1:
fig.savefig(f"figures/11_time_res_vs_snr_tdt{template_dt:0.1e}.pdf")