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

Pulse finding for multiple SNR

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-04-24 18:37:13 +02:00
parent 81b3502de5
commit a011aee28e
1 changed files with 241 additions and 211 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

452
simulations/11_pulsed_timing.py
View file

@ -136,9 +136,9 @@ if __name__ == "__main__":
bp_freq = (30e-3, 80e-3) # GHz bp_freq = (30e-3, 80e-3) # GHz
template_dt = 5e-2 # ns template_dt = 5e-2 # ns
template_length = 500 # ns template_length = 500 # ns
noise_sigma_factor = 1e-1 # amplitude factor
N_residuals = 50*3 if len(sys.argv) < 2 else int(sys.argv[1]) 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
antenna_dt = 2 # ns antenna_dt = 2 # ns
antenna_timelength = 2048 # ns antenna_timelength = 2048 # ns
@ -176,254 +176,284 @@ if __name__ == "__main__":
if True: if True:
plt.close(fig) plt.close(fig)
# time_accuracies = np.zeros(len(noise_factors))
# Find difference between true and templated times for k, noise_sigma_factor in tqdm(enumerate(noise_factors)):
# print() #separating tqdm
time_residuals = np.zeros(N_residuals) #
for j in tqdm(range(N_residuals)): # Find difference between true and templated times
do_plots = j==0 #
time_residuals = np.zeros(N_residuals)
for j in tqdm(range(N_residuals)):
do_plots = j==0
# receive at antenna # receive at antenna
## place the deltapeak signal at a random location ## place the deltapeak signal at a random location
antenna = Waveform(None, dt=antenna_dt, name='Signal') antenna = Waveform(None, dt=antenna_dt, name='Signal')
if not True: # Create antenna trace without template if not True: # Create antenna trace without template
antenna_true_signal, antenna_peak_sample = util.deltapeak(timelength=antenna_timelength, samplerate=1/antenna.dt, offset=[0.2, 0.8], rng=rng) antenna_true_signal, antenna_peak_sample = util.deltapeak(timelength=antenna_timelength, samplerate=1/antenna.dt, offset=[0.2, 0.8], rng=rng)
antenna.peak_sample = antenna_peak_sample antenna.peak_sample = antenna_peak_sample
antenna.peak_time = antenna.dt * antenna.peak_sample antenna.peak_time = antenna.dt * antenna.peak_sample
antenna.signal = antenna_bp(antenna.signal, *bp_freq, antenna.dt) antenna.signal = antenna_bp(antenna.signal, *bp_freq, antenna.dt)
print(f"Antenna Peak Time: {antenna.peak_time}")
print(f"Antenna Peak Sample: {antenna.peak_sample}")
else: # Sample the template at some offset else: # Sample the template at some offset
antenna.peak_time = antenna_timelength * ((0.8 - 0.2) *rng.random(1) + 0.2) antenna.peak_time = antenna_timelength * ((0.8 - 0.2) *rng.random(1) + 0.2)
sampling_offset = rng.random(1)*antenna.dt sampling_offset = rng.random(1)*antenna.dt
antenna.t = util.sampled_time(1/antenna.dt, start=0, end=antenna_timelength) antenna.t = util.sampled_time(1/antenna.dt, start=0, end=antenna_timelength)
antenna.signal = interp1d_template(antenna.t - antenna.peak_time) antenna.signal = interp1d_template(antenna.t - antenna.peak_time)
antenna.peak_sample = antenna.peak_time/antenna.dt antenna.peak_sample = antenna.peak_time/antenna.dt
antenna_true_signal = antenna.signal antenna_true_signal = antenna.signal
true_time_offset = antenna.peak_time - template.peak_time
if do_plots: true_time_offset = antenna.peak_time - template.peak_time
print(f"Antenna Peak Time: {antenna.peak_time}")
print(f"Antenna Peak Sample: {antenna.peak_sample}")
if False: # flip polarisation if False: # flip polarisation
antenna.signal *= -1 antenna.signal *= -1
## Add noise ## Add noise
noise_amplitude = max(template.signal) * noise_sigma_factor noise_amplitude = max(template.signal) * noise_sigma_factor
noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal)) noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal))
filtered_noise = antenna_bp(noise_realisation, *bp_freq, antenna.dt) filtered_noise = antenna_bp(noise_realisation, *bp_freq, antenna.dt)
antenna.signal += filtered_noise antenna.signal += filtered_noise
if do_plots: # show signals if do_plots: # show signals
fig, axs = plt.subplots(2, sharex=True) fig, axs = plt.subplots(2, sharex=True)
axs[0].set_title("Antenna Waveform") axs[0].set_title("Antenna Waveform")
axs[-1].set_xlabel("Time [ns]") axs[-1].set_xlabel("Time [ns]")
axs[0].set_ylabel("Amplitude") axs[0].set_ylabel("Amplitude")
axs[0].plot(antenna.t, antenna.signal, label='bandpassed w/ noise', alpha=0.9) axs[0].plot(antenna.t, antenna.signal, label='bandpassed w/ noise', alpha=0.9)
axs[0].plot(antenna.t, antenna.signal - filtered_noise, label='bandpassed w/o noise', alpha=0.9) axs[0].plot(antenna.t, antenna.signal - filtered_noise, label='bandpassed w/o noise', alpha=0.9)
axs[0].legend() axs[0].legend()
axs[1].set_title("Template") axs[1].set_title("Template")
axs[1].set_ylabel("Amplitude") axs[1].set_ylabel("Amplitude")
axs[1].plot(template.t, template.signal, label='orig') 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='true moved orig')
axs[1].legend() axs[1].legend()
axs[0].grid() axs[0].grid()
axs[1].grid() axs[1].grid()
fig.savefig('figures/11_antenna_signals.pdf') fig.savefig('figures/11_antenna_signals.pdf')
if True: # zoom if True: # zoom
wx = 100 wx = 100
x0 = true_time_offset x0 = true_time_offset
old_xlims = axs[0].get_xlim() old_xlims = axs[0].get_xlim()
axs[0].set_xlim( x0-wx, x0+wx) axs[0].set_xlim( x0-wx, x0+wx)
fig.savefig('figures/11_antenna_signals_zoom.pdf') fig.savefig('figures/11_antenna_signals_zoom.pdf')
# restore # restore
axs[0].set_xlim(*old_xlims) axs[0].set_xlim(*old_xlims)
if False:
plt.close(fig)
axs2 = None
if True: # upsampled trace
upsampled_trace, upsampled_t = trace_upsampler(antenna.signal, template.t, antenna.t)
if do_plots: # Show upsampled traces
fig2, axs2 = plt.subplots(1, sharex=True)
if not hasattr(axs2, '__len__'):
axs2 = [axs2]
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(upsampled_t, upsampled_trace, label='upsampled')
axs2[0].legend(loc='upper right')
fig2.savefig('figures/11_upsampled.pdf')
wx = 1e2
x0 = upsampled_t[0] + wx - 5
axs2[0].set_xlim(x0-wx, x0+wx)
fig2.savefig('figures/11_upsampled_zoom.pdf')
if True: if True:
plt.close(fig2) plt.close(fig)
# determine correlations with arguments axs2 = None
lag_dt = upsampled_t[1] - upsampled_t[0] if True: # upsampled trace
corrs, (out1_signal, out2_template, lags) = my_correlation(upsampled_trace, template.signal) upsampled_trace, upsampled_t = trace_upsampler(antenna.signal, template.t, antenna.t)
else: # downsampled template if do_plots: # Show upsampled traces
raise NotImplementedError fig2, axs2 = plt.subplots(1, sharex=True)
if not hasattr(axs2, '__len__'):
axs2 = [axs2]
corrs, (out1_signal, out2_template, lags) = my_downsampling_correlation(template.signal, antenna.signal, template.t, antenna.t) axs2[-1].set_xlabel("Time [ns]")
lag_dt = upsampled_t[1] - upsampled_t[0] axs2[0].set_ylabel("Amplitude")
axs2[0].plot(antenna.t, antenna.signal, marker='o', label='orig')
axs2[0].plot(upsampled_t, upsampled_trace, label='upsampled')
axs2[0].legend(loc='upper right')
# Determine best correlation time fig2.savefig('figures/11_upsampled.pdf')
idx = np.argmax(abs(corrs))
best_sample_lag = lags[idx]
best_time_lag = best_sample_lag * lag_dt
time_residuals[j] = best_time_lag - true_time_offset
if not do_plots: wx = 1e2
continue x0 = upsampled_t[0] + wx - 5
axs2[0].set_xlim(x0-wx, x0+wx)
fig2.savefig('figures/11_upsampled_zoom.pdf')
if do_plots and axs2: if True:
axs2[-1].axvline(best_time_lag, color='r', alpha=0.5, linewidth=2) plt.close(fig2)
axs2[-1].axvline(true_time_offset, color='g', alpha=0.5, linewidth=2)
# Show the final signals correlated # determine correlations with arguments
if do_plots: lag_dt = upsampled_t[1] - upsampled_t[0]
# amplitude scaling required for single axis plotting corrs, (out1_signal, out2_template, lags) = my_correlation(upsampled_trace, template.signal)
template_amp_scaler = max(abs(template.signal)) / max(abs(antenna.signal))
# start the figure else: # downsampled template
fig, axs = plt.subplots(2, sharex=True) raise NotImplementedError
ylabel_kwargs = dict(
#rotation=0,
ha='right',
va='center'
)
axs[-1].set_xlabel("Time [ns]")
offset_list = [ corrs, (out1_signal, out2_template, lags) = my_downsampling_correlation(template.signal, antenna.signal, template.t, antenna.t)
[best_time_lag, dict(label=template.name, color='orange')], lag_dt = upsampled_t[1] - upsampled_t[0]
[true_time_offset, dict(label='True offset', color='green')],
]
# Signal # Determine best correlation time
i=0 idx = np.argmax(abs(corrs))
axs[i].set_ylabel("Amplitude", **ylabel_kwargs) best_sample_lag = lags[idx]
axs[i].plot(antenna.t, antenna.signal, label=antenna.name) best_time_lag = best_sample_lag * lag_dt
time_residuals[j] = best_time_lag - true_time_offset
# Plot the template if not do_plots:
for offset_args in offset_list: continue
this_kwargs = offset_args[1]
offset = offset_args[0]
l = axs[i].plot(offset + template.t, template_amp_scaler * template.signal, **this_kwargs) if do_plots and axs2:
axs2[-1].axvline(best_time_lag, color='r', alpha=0.5, linewidth=2)
axs2[-1].axvline(true_time_offset, color='g', alpha=0.5, linewidth=2)
axs[i].legend() # Show the final signals correlated
if do_plots:
# amplitude scaling required for single axis plotting
template_amp_scaler = max(abs(template.signal)) / max(abs(antenna.signal))
# Correlation # start the figure
i=1 fig, axs = plt.subplots(2, sharex=True)
axs[i].set_ylabel("Correlation", **ylabel_kwargs) ylabel_kwargs = dict(
axs[i].plot(lags * lag_dt, corrs) #rotation=0,
ha='right',
va='center'
)
axs[-1].set_xlabel("Time [ns]")
# Lines across both axes offset_list = [
for offset_args in offset_list: [best_time_lag, dict(label=template.name, color='orange')],
this_kwargs = offset_args[1] [true_time_offset, dict(label='True offset', color='green')],
offset = offset_args[0]
for i in [0,1]:
axs[i].axvline(offset, ls='--', color=this_kwargs['color'], alpha=0.7)
axs[0].axvline(offset + len(template.signal) * (template.t[1] - template.t[0]), color=this_kwargs['color'], alpha=0.7)
if True: # zoom
wx = len(template.signal) * (template.dt)/2
t0 = best_time_lag
old_xlims = axs[0].get_xlim()
axs[i].set_xlim( x0-wx, x0+3*wx)
fig.savefig('figures/11_corrs_zoom.pdf')
# restore
axs[i].set_xlim(*old_xlims)
fig.tight_layout()
fig.savefig('figures/11_corrs.pdf')
if False:
plt.close(fig)
# Make a plot of the time residuals
if len(time_residuals) > 1:
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*1e3: .1e}ps"
+ f"; antenna dt: {antenna.dt: .1e}ns"
+ f"; noise_factor: {noise_sigma_factor: .1e}"
)
ax.set_xlabel("Time Residual [ns]")
ax.set_ylabel("#")
counts, bins, _patches = ax.hist(time_residuals, **hist_kwargs)
if True: # fit gaussian to histogram
min_x = min(time_residuals)
max_x = max(time_residuals)
dx = bins[1] - bins[0]
scale = len(time_residuals) * dx
xs = np.linspace(min_x, max_x)
# do the fit
name = "Norm"
param_names = [ "$\\mu$", "$\\sigma$" ]
distr_func = stats.norm
label = name +"(" + ','.join(param_names) + ')'
# plot
fit_params = distr_func.fit(time_residuals)
fit_ys = scale * distr_func.pdf(xs, *fit_params)
ax.plot(xs, fit_ys, label=label)
# chisq
ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts)
if True:
ct *= np.sum(counts)/np.sum(ct)
c2t = stats.chisquare(counts, ct, ddof=len(fit_params))
chisq_strs = [
f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}"
] ]
# text on plot # Signal
text_str = "\n".join( i=0
[label] axs[i].set_ylabel("Amplitude", **ylabel_kwargs)
+ axs[i].plot(antenna.t, antenna.signal, label=antenna.name)
[ f"{param} = {value: .2e}" for param, value in zip_longest(param_names, fit_params, fillvalue='?') ]
+
chisq_strs
)
ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes) # Plot the template
for offset_args in offset_list:
this_kwargs = offset_args[1]
offset = offset_args[0]
fig.savefig("figures/11_time_residual_hist.pdf") l = axs[i].plot(offset + template.t, template_amp_scaler * template.signal, **this_kwargs)
axs[i].legend()
# Correlation
i=1
axs[i].set_ylabel("Correlation", **ylabel_kwargs)
axs[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 i in [0,1]:
axs[i].axvline(offset, ls='--', color=this_kwargs['color'], alpha=0.7)
axs[0].axvline(offset + len(template.signal) * (template.t[1] - template.t[0]), color=this_kwargs['color'], alpha=0.7)
if True: # zoom
wx = len(template.signal) * (template.dt)/2
t0 = best_time_lag
old_xlims = axs[0].get_xlim()
axs[i].set_xlim( x0-wx, x0+3*wx)
fig.savefig('figures/11_corrs_zoom.pdf')
# restore
axs[i].set_xlim(*old_xlims)
fig.tight_layout()
fig.savefig('figures/11_corrs.pdf')
if True:
plt.close(fig)
print()# separating tqdm
# Make a plot of the time residuals
if len(time_residuals) > 1:
time_accuracies[k] = np.std(time_residuals)
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*1e3: .1e}ps"
+ f"; antenna dt: {antenna.dt: .1e}ns"
+ f"; noise_factor: {noise_sigma_factor: .1e}"
)
ax.set_xlabel("Time Residual [ns]")
ax.set_ylabel("#")
counts, bins, _patches = ax.hist(time_residuals, **hist_kwargs)
if True: # fit gaussian to histogram
min_x = min(time_residuals)
max_x = max(time_residuals)
dx = bins[1] - bins[0]
scale = len(time_residuals) * dx
xs = np.linspace(min_x, max_x)
# do the fit
name = "Norm"
param_names = [ "$\\mu$", "$\\sigma$" ]
distr_func = stats.norm
label = name +"(" + ','.join(param_names) + ')'
# plot
fit_params = distr_func.fit(time_residuals)
fit_ys = scale * distr_func.pdf(xs, *fit_params)
ax.plot(xs, fit_ys, label=label)
# chisq
ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts)
if True:
ct *= np.sum(counts)/np.sum(ct)
c2t = stats.chisquare(counts, ct, ddof=len(fit_params))
chisq_strs = [
f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}"
]
# text on plot
text_str = "\n".join(
[label]
+
[ f"{param} = {value: .2e}" for param, value in zip_longest(param_names, fit_params, fillvalue='?') ]
+
chisq_strs
)
ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes)
fig.savefig("figures/11_time_residual_hist_{noise_sigma_factor: .1e}.pdf")
if True:
plt.close(fig)
# SNR time accuracy plot
if True:
fig, ax = plt.subplots()
ax.set_title("Template matching SNR vs time accuracy")
ax.set_xlabel("Signal to Noise Factor")
ax.set_ylabel("Time Accuracy [ns]")
if True:
ax.set_xscale('log')
ax.set_yscale('log')
# plot the values
ax.plot(1/np.asarray(noise_factors), time_accuracies, ls='none', marker='o')
# Set horizontal line at 1 ns
if True:
ax.axhline(1, ls='--', alpha=0.8, color='g')
ax.grid()
fig.tight_layout()
fig.savefig("figures/11_time_res_vs_snr.pdf")
plt.show() plt.show()