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: move timeresidual matching to function

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-04-26 15:45:42 +02:00
parent 1f00a3fe76
commit 168b0a60bc
1 changed files with 227 additions and 211 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

184
simulations/11_pulsed_timing.py
View file

@ -172,82 +172,17 @@ def create_template(dt=1, timelength=1, bp_freq=(0, np.inf), name=None, normalis
return template, _deltapeak return template, _deltapeak
if __name__ == "__main__": def get_time_residuals_for_template(
import os N_residuals, template, interpolation_template=None,
import matplotlib antenna_dt=1, antenna_timelength=100,
import sys snr_sigma_factor=10,bp_freq=(0,np.inf),
if os.name == 'posix' and "DISPLAY" not in os.environ: normalise_noise=False, h5_cache_fname=None, read_cache=True, write_cache=None,
matplotlib.use('Agg') rng=rng, tqdm=tqdm,
):
bp_freq = (30e-3, 80e-3) # GHz
template_dt = 5e-2 # 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])
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 = 1024 # ns
cut_wrong_peak_matches = True
normalise_noise = False
h5_cache_fname = f'11_pulsed_timing.hdf5'
#
# Interpolation Template
# 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)
if True: # show interpolation template
fig, ax = plt.subplots()
ax.set_title("Deltapeak and Bandpassed Template")
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Amplitude")
ax.plot(interp_template.t, max(interp_template.signal)*_deltapeak[0], label='Impulse Template')
ax.plot(interp_template.t, interp_template.signal, label='Filtered Template')
ax.legend()
fig.savefig('figures/11_interpolation_deltapeak.pdf')
if True: # show filtering equivalence samplerates
_deltapeak = util.deltapeak(timelength=template_length, samplerate=1/antenna_dt, offset=0)
_time = util.sampled_time(end=template_length, sample_rate=1/antenna_dt)
_bandpassed = antenna_bp(_deltapeak[0], *bp_freq, antenna_dt)
ax.plot(_time, max(_bandpassed)*_deltapeak[0], label='Impulse Antenna')
ax.plot(_time, _bandpassed, label='Filtered Antenna')
ax.legend()
fig.savefig('figures/11_interpolation_deltapeak+antenna.pdf')
if True:
plt.close(fig)
#
# Create the template
# This is sampled at a lower samplerate than the interpolation template
#
template, _ = create_template(dt=template_dt, timelength=template_length, bp_freq=bp_freq, name='Template')
#
# Find time accuracies as a function of signal strength
#
time_accuracies = np.zeros(len(snr_factors))
mask_counts = np.zeros(len(snr_factors))
for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
# Read in cached time residuals # Read in cached time residuals
if True: if read_cache:
cached_time_residuals = read_time_residuals_cache(h5_cache_fname, template.dt, antenna_dt, snr_sigma_factor) cached_time_residuals = read_time_residuals_cache(h5_cache_fname, template.dt, antenna_dt, snr_sigma_factor)
else: else:
cached_time_residuals = np.array([]) cached_time_residuals = np.array([])
@ -261,8 +196,7 @@ if __name__ == "__main__":
# 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 interpolation_template is None: # Create antenna trace without interpolation template
if False: # Create antenna trace without interpolation 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
@ -278,7 +212,7 @@ if __name__ == "__main__":
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)
# Sample the interpolation template # Sample the interpolation template
antenna.signal = interp_template.interpolate(antenna.t - antenna.peak_time) antenna.signal = interpolation_template.interpolate(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
@ -295,7 +229,8 @@ if __name__ == "__main__":
antenna.signal += filtered_noise antenna.signal += filtered_noise
if do_plots: # show signals # Show signals
if do_plots:
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]")
@ -332,7 +267,6 @@ if __name__ == "__main__":
axs2 = None axs2 = None
if True: # upsampled trace if True: # upsampled trace
upsampled_trace, upsampled_t = trace_upsampler(antenna.signal, template.t, antenna.t) upsampled_trace, upsampled_t = trace_upsampler(antenna.signal, template.t, antenna.t)
if do_plots: # Show upsampled traces if do_plots: # Show upsampled traces
fig2, axs2 = plt.subplots(1, sharex=True) fig2, axs2 = plt.subplots(1, sharex=True)
if not hasattr(axs2, '__len__'): if not hasattr(axs2, '__len__'):
@ -445,23 +379,105 @@ if __name__ == "__main__":
if True: if True:
plt.close(fig) plt.close(fig)
print()# separating tqdm
print()# separating tqdm
# Were new time residuals calculated? # Were new time residuals calculated?
# Add them to the cache file # Add them to the cache file
if len(time_residuals) > 1: if len(time_residuals) > 1:
# merge cached and calculated time residuals # merge cached and calculated time residuals
time_residuals = np.concatenate((cached_time_residuals, time_residuals), axis=None) time_residuals = np.concatenate((cached_time_residuals, time_residuals), axis=None)
if True: # write the cache if write_cache or read_cache and write_cache is None: # write the cache
write_time_residuals_cache(h5_cache_fname, time_residuals, template_dt, antenna_dt, snr_sigma_factor) write_time_residuals_cache(h5_cache_fname, time_residuals, template_dt, antenna_dt, snr_sigma_factor)
else: else:
time_residuals = cached_time_residuals time_residuals = cached_time_residuals
# Only return N_residuals (even if more have been cached)
return time_residuals[:N_residuals]
if __name__ == "__main__":
import os
import matplotlib
import sys
if os.name == 'posix' and "DISPLAY" not in os.environ:
matplotlib.use('Agg')
bp_freq = (30e-3, 80e-3) # GHz
template_dt = 5e-2 # 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])
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, dtype=float)
antenna_dt = 2 # ns
antenna_timelength = 1024 # ns
cut_wrong_peak_matches = True
normalise_noise = False
h5_cache_fname = f'11_pulsed_timing.hdf5'
#
# Interpolation Template
# 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)
if True: # show interpolation template
fig, ax = plt.subplots()
ax.set_title("Deltapeak and Bandpassed Template")
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Amplitude")
ax.plot(interp_template.t, max(interp_template.signal)*_deltapeak[0], label='Impulse Template')
ax.plot(interp_template.t, interp_template.signal, label='Filtered Template')
ax.legend()
fig.savefig('figures/11_interpolation_deltapeak.pdf')
if True: # show filtering equivalence samplerates
_deltapeak = util.deltapeak(timelength=template_length, samplerate=1/antenna_dt, offset=0)
_time = util.sampled_time(end=template_length, sample_rate=1/antenna_dt)
_bandpassed = antenna_bp(_deltapeak[0], *bp_freq, antenna_dt)
ax.plot(_time, max(_bandpassed)*_deltapeak[0], label='Impulse Antenna')
ax.plot(_time, _bandpassed, label='Filtered Antenna')
ax.legend()
fig.savefig('figures/11_interpolation_deltapeak+antenna.pdf')
if True:
plt.close(fig)
#
# Create the template
# This is sampled at a lower samplerate than the interpolation template
#
template, _ = create_template(dt=template_dt, timelength=template_length, bp_freq=bp_freq, name='Template')
#
# Find time accuracies as a function of signal strength
#
time_accuracies = np.zeros(len(snr_factors))
mask_counts = np.zeros(len(snr_factors))
for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
time_residuals = get_time_residuals_for_template(
N_residuals, template, interpolation_template=interp_template,
antenna_dt=antenna_dt, antenna_timelength=antenna_timelength,
snr_sigma_factor=snr_sigma_factor, bp_freq=bp_freq, normalise_noise=normalise_noise,
h5_cache_fname=h5_cache_fname, rng=rng, tqdm=tqdm)
print()# separating tqdm
print()# separating tqdm
# Make a plot of the time residuals # Make a plot of the time residuals
if N_residuals > 1: if N_residuals > 1:
time_residuals = time_residuals[:N_residuals]
for i in range(1 + cut_wrong_peak_matches): for i in range(1 + cut_wrong_peak_matches):
mask_count = 0 mask_count = 0