From 8f42db2a9969cb01e23dfe728eff8adcf498aab6 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 24 Apr 2023 16:59:54 +0200 Subject: [PATCH] Pulse: switched filter + sample template for fine time offset --- simulations/11_pulsed_timing.py | 98 +++++++++++++++++++++------------ 1 file changed, 64 insertions(+), 34 deletions(-) diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index 77c9840..e67c8ba 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -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) def antenna_bp(trace, low_bp, high_bp, dt, order=3): - # order < 6 for stability 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) - bandpassed = signal.lfilter(*bp_filter, trace) + bp_filter = signal.butter(order, [low_bp, high_bp], 'band', fs=fs, output='sos') + bandpassed = signal.sosfilt(bp_filter, trace) return bandpassed @@ -139,7 +135,7 @@ if __name__ == "__main__": bp_freq = (30e-3, 80e-3) # GHz template_dt = 5e-2 # 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]) @@ -150,21 +146,32 @@ if __name__ == "__main__": # Create the template # template = Waveform(None, dt=template_dt, name='Template') - _deltapeak = util.deltapeak(timelength=template_length, samplerate=1/template.dt, offset=10/template.dt) - template.signal = antenna_bp(_deltapeak[0], *bp_freq, (np.sqrt(antenna_dt/template_dt))*template_dt) # TODO: fix sqrt constant + _deltapeak = util.deltapeak(timelength=template_length, samplerate=1/template.dt, offset=0) + template.signal = antenna_bp(_deltapeak[0], *bp_freq, template_dt) template.peak_sample = _deltapeak[1] 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 fig, ax = plt.subplots() ax.set_title("Deltapeak and Bandpassed Template") ax.set_xlabel("Time [ns]") ax.set_ylabel("Amplitude") - ax.plot(template.t, max(template.signal)*_deltapeak[0]) - ax.plot(template.t, template.signal) + ax.plot(template.t, max(template.signal)*_deltapeak[0], label='Impulse Template') + ax.plot(template.t, template.signal, label='Filtered Template') 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: plt.close(fig) @@ -178,29 +185,41 @@ if __name__ == "__main__": # receive at antenna ## place the deltapeak signal at a random location 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 - antenna.peak_sample = antenna_peak_sample - antenna.peak_time = antenna.dt * antenna.peak_sample + + 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.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: print(f"Antenna Peak Time: {antenna.peak_time}") 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 antenna.signal *= -1 ## Add noise noise_amplitude = max(template.signal) * noise_sigma_factor 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 = antenna_bp(antenna_unfiltered_signal, *bp_freq, antenna.dt) - - true_time_offset = antenna.peak_time - template.peak_time + antenna.signal += filtered_noise if do_plots: # show signals fig, axs = plt.subplots(2, sharex=True) @@ -208,8 +227,7 @@ if __name__ == "__main__": axs[-1].set_xlabel("Time [ns]") 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_unfiltered_signal, label='true signal w/ noise', alpha=0.9) - axs[0].plot(antenna.t, antenna_true_signal, label='true signal 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[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].legend() + axs[0].grid() + axs[1].grid() + fig.savefig('figures/11_antenna_signals.pdf') if True: # zoom @@ -277,9 +298,14 @@ if __name__ == "__main__": 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) # 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)) + + # start the figure fig, axs = plt.subplots(2, sharex=True) ylabel_kwargs = dict( #rotation=0, @@ -298,25 +324,30 @@ if __name__ == "__main__": axs[i].set_ylabel("Amplitude", **ylabel_kwargs) axs[i].plot(antenna.t, antenna.signal, label=antenna.name) - # Put template on an twinned axis (magnitudes are different) - ax2 = axs[i].twinx() - for i, offset_args in enumerate(offset_list): + # Plot the template + for offset_args in offset_list: this_kwargs = offset_args[1] offset = offset_args[0] - ax2.axvline(offset + len(template.signal) * (template.t[1] - template.t[0]), color=this_kwargs['color'], alpha=0.7) - ax2.axvline(offset, color=this_kwargs['color'], alpha=0.7) - ax2.plot(offset + template.t, template.signal, **this_kwargs) + 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) - for i, offset_args in enumerate(offset_list): + + # Lines across both axes + for offset_args in offset_list: this_kwargs = offset_args[1] 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 wx = len(template.signal) * (template.dt)/2 @@ -330,7 +361,6 @@ if __name__ == "__main__": axs[i].set_xlim(*old_xlims) fig.tight_layout() - fig.legend() fig.savefig('figures/11_corrs.pdf') if False: