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: switched filter + sample template for fine time offset

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-04-24 16:59:54 +02:00
parent 354ec93a79
commit 8f42db2a99
1 changed files with 64 additions and 34 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

98
simulations/11_pulsed_timing.py
View file

@ -54,14 +54,10 @@ def white_noise_realisation(N_samples, noise_sigma=1, rng=rng):
return rng.normal(0, noise_sigma or 0, size=N_samples) return rng.normal(0, noise_sigma or 0, size=N_samples)
def antenna_bp(trace, low_bp, high_bp, dt, order=3): def antenna_bp(trace, low_bp, high_bp, dt, order=3):
# order < 6 for stability
fs = 1/dt fs = 1/dt
nyq = 1* fs
low_bp = low_bp / nyq
high_bp = high_bp / nyq
bp_filter = signal.butter(order, [low_bp, high_bp], 'band', fs=fs) bp_filter = signal.butter(order, [low_bp, high_bp], 'band', fs=fs, output='sos')
bandpassed = signal.lfilter(*bp_filter, trace) bandpassed = signal.sosfilt(bp_filter, trace)
return bandpassed return bandpassed
@ -139,7 +135,7 @@ 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 = 1e0 # keep between 10 and 0.1 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])
@ -150,21 +146,32 @@ if __name__ == "__main__":
# Create the template # Create the template
# #
template = Waveform(None, dt=template_dt, name='Template') template = Waveform(None, dt=template_dt, name='Template')
_deltapeak = util.deltapeak(timelength=template_length, samplerate=1/template.dt, offset=10/template.dt) _deltapeak = util.deltapeak(timelength=template_length, samplerate=1/template.dt, offset=0)
template.signal = antenna_bp(_deltapeak[0], *bp_freq, (np.sqrt(antenna_dt/template_dt))*template_dt) # TODO: fix sqrt constant template.signal = antenna_bp(_deltapeak[0], *bp_freq, template_dt)
template.peak_sample = _deltapeak[1] template.peak_sample = _deltapeak[1]
template.peak_time = template.dt * template.peak_sample 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)
if True: # show template if True: # show template
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.set_title("Deltapeak and Bandpassed Template") ax.set_title("Deltapeak and Bandpassed Template")
ax.set_xlabel("Time [ns]") ax.set_xlabel("Time [ns]")
ax.set_ylabel("Amplitude") ax.set_ylabel("Amplitude")
ax.plot(template.t, max(template.signal)*_deltapeak[0]) ax.plot(template.t, max(template.signal)*_deltapeak[0], label='Impulse Template')
ax.plot(template.t, template.signal) ax.plot(template.t, template.signal, label='Filtered Template')
fig.savefig('figures/11_template_deltapeak.pdf') fig.savefig('figures/11_template_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_template_deltapeak+antenna.pdf')
if True: if True:
plt.close(fig) plt.close(fig)
@ -178,29 +185,41 @@ 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')
antenna_true_signal, antenna_peak_sample = util.deltapeak(timelength=antenna_timelength, samplerate=1/antenna.dt, offset=[0.2, 0.8], rng=rng)
antenna.signal = antenna_true_signal if not True: # Create antenna trace without template
antenna.peak_sample = antenna_peak_sample antenna_true_signal, antenna_peak_sample = util.deltapeak(timelength=antenna_timelength, samplerate=1/antenna.dt, offset=[0.2, 0.8], rng=rng)
antenna.peak_time = antenna.dt * antenna.peak_sample
antenna.peak_sample = antenna_peak_sample
antenna.peak_time = antenna.dt * antenna.peak_sample
antenna.signal = antenna_bp(antenna.signal, *bp_freq, antenna.dt)
else: # Sample the template at some offset
antenna.peak_time = antenna_timelength * ((0.8 - 0.2) *rng.random(1) + 0.2)
sampling_offset = rng.random(1)*antenna.dt
antenna.t = util.sampled_time(1/antenna.dt, start=0, end=antenna_timelength)
antenna.signal = interp1d_template(antenna.t - antenna.peak_time)
antenna.peak_sample = antenna.peak_time/antenna.dt
antenna_true_signal = antenna.signal
true_time_offset = antenna.peak_time - template.peak_time
if do_plots: if do_plots:
print(f"Antenna Peak Time: {antenna.peak_time}") print(f"Antenna Peak Time: {antenna.peak_time}")
print(f"Antenna Peak Sample: {antenna.peak_sample}") print(f"Antenna Peak Sample: {antenna.peak_sample}")
if False: # bandpass when emitting the signal
antenna.signal = antenna_bp(antenna.signal, *bp_freq, antenna.dt)
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)
antenna_unfiltered_signal = antenna.signal + noise_realisation antenna.signal += filtered_noise
antenna.signal = antenna_bp(antenna_unfiltered_signal, *bp_freq, antenna.dt)
true_time_offset = antenna.peak_time - template.peak_time
if do_plots: # show signals if do_plots: # show signals
fig, axs = plt.subplots(2, sharex=True) fig, axs = plt.subplots(2, sharex=True)
@ -208,8 +227,7 @@ if __name__ == "__main__":
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_unfiltered_signal, label='true signal 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_true_signal, label='true signal w/o noise', alpha=0.9)
axs[0].legend() axs[0].legend()
axs[1].set_title("Template") axs[1].set_title("Template")
@ -218,6 +236,9 @@ if __name__ == "__main__":
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[1].grid()
fig.savefig('figures/11_antenna_signals.pdf') fig.savefig('figures/11_antenna_signals.pdf')
if True: # zoom if True: # zoom
@ -277,9 +298,14 @@ if __name__ == "__main__":
if do_plots and axs2: if do_plots and axs2:
axs2[-1].axvline(best_time_lag, color='r', alpha=0.5, linewidth=2) 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)
# Show the final signals correlated # Show the final signals correlated
if do_plots: if do_plots:
# amplitude scaling required for single axis plotting
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)
ylabel_kwargs = dict( ylabel_kwargs = dict(
#rotation=0, #rotation=0,
@ -298,25 +324,30 @@ if __name__ == "__main__":
axs[i].set_ylabel("Amplitude", **ylabel_kwargs) axs[i].set_ylabel("Amplitude", **ylabel_kwargs)
axs[i].plot(antenna.t, antenna.signal, label=antenna.name) axs[i].plot(antenna.t, antenna.signal, label=antenna.name)
# Put template on an twinned axis (magnitudes are different) # Plot the template
ax2 = axs[i].twinx() for offset_args in offset_list:
for i, offset_args in enumerate(offset_list):
this_kwargs = offset_args[1] this_kwargs = offset_args[1]
offset = offset_args[0] offset = offset_args[0]
ax2.axvline(offset + len(template.signal) * (template.t[1] - template.t[0]), color=this_kwargs['color'], alpha=0.7) l = axs[i].plot(offset + template.t, template_amp_scaler * template.signal, **this_kwargs)
ax2.axvline(offset, color=this_kwargs['color'], alpha=0.7)
ax2.plot(offset + template.t, template.signal, **this_kwargs) axs[i].legend()
# Correlation # Correlation
i=1 i=1
axs[i].set_ylabel("Correlation", **ylabel_kwargs) axs[i].set_ylabel("Correlation", **ylabel_kwargs)
axs[i].plot(lags * lag_dt, corrs) axs[i].plot(lags * lag_dt, corrs)
for i, offset_args in enumerate(offset_list):
# Lines across both axes
for offset_args in offset_list:
this_kwargs = offset_args[1] this_kwargs = offset_args[1]
offset = offset_args[0] offset = offset_args[0]
axs[i].axvline(offset, ls='--', **this_kwargs) 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 if True: # zoom
wx = len(template.signal) * (template.dt)/2 wx = len(template.signal) * (template.dt)/2
@ -330,7 +361,6 @@ if __name__ == "__main__":
axs[i].set_xlim(*old_xlims) axs[i].set_xlim(*old_xlims)
fig.tight_layout() fig.tight_layout()
fig.legend()
fig.savefig('figures/11_corrs.pdf') fig.savefig('figures/11_corrs.pdf')
if False: if False: