From a011aee28e96534f48bdbbbcdc0461ccde53c18a Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 24 Apr 2023 18:37:13 +0200 Subject: [PATCH] Pulse finding for multiple SNR --- simulations/11_pulsed_timing.py | 452 +++++++++++++++++--------------- 1 file changed, 241 insertions(+), 211 deletions(-) diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index ae4eb38..b242bfd 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -136,9 +136,9 @@ if __name__ == "__main__": bp_freq = (30e-3, 80e-3) # GHz template_dt = 5e-2 # 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]) + 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_timelength = 2048 # ns @@ -176,254 +176,284 @@ if __name__ == "__main__": if True: plt.close(fig) - # - # Find difference between true and templated times - # - time_residuals = np.zeros(N_residuals) - for j in tqdm(range(N_residuals)): - do_plots = j==0 + time_accuracies = np.zeros(len(noise_factors)) + for k, noise_sigma_factor in tqdm(enumerate(noise_factors)): + print() #separating tqdm + # + # Find difference between true and templated times + # + time_residuals = np.zeros(N_residuals) + for j in tqdm(range(N_residuals)): + do_plots = j==0 - # receive at antenna - ## place the deltapeak signal at a random location - antenna = Waveform(None, dt=antenna_dt, name='Signal') + # receive at antenna + ## place the deltapeak signal at a random location + antenna = Waveform(None, dt=antenna_dt, name='Signal') - 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) + 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) + antenna.peak_sample = antenna_peak_sample + antenna.peak_time = antenna.dt * antenna.peak_sample + 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 - antenna.peak_time = antenna_timelength * ((0.8 - 0.2) *rng.random(1) + 0.2) - sampling_offset = rng.random(1)*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.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_true_signal = antenna.signal - true_time_offset = antenna.peak_time - template.peak_time + antenna.peak_sample = antenna.peak_time/antenna.dt + antenna_true_signal = antenna.signal - if do_plots: - print(f"Antenna Peak Time: {antenna.peak_time}") - print(f"Antenna Peak Sample: {antenna.peak_sample}") + true_time_offset = antenna.peak_time - template.peak_time - if False: # flip polarisation - antenna.signal *= -1 + 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) + ## 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.signal += filtered_noise + antenna.signal += filtered_noise - if do_plots: # show signals - fig, axs = plt.subplots(2, sharex=True) - axs[0].set_title("Antenna Waveform") - 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.signal - filtered_noise, label='bandpassed w/o noise', alpha=0.9) - axs[0].legend() + if do_plots: # show signals + fig, axs = plt.subplots(2, sharex=True) + axs[0].set_title("Antenna Waveform") + 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.signal - filtered_noise, label='bandpassed w/o noise', alpha=0.9) + axs[0].legend() - 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].legend() + 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].legend() - axs[0].grid() - axs[1].grid() + axs[0].grid() + axs[1].grid() - fig.savefig('figures/11_antenna_signals.pdf') + fig.savefig('figures/11_antenna_signals.pdf') - if True: # zoom - wx = 100 - x0 = true_time_offset + if True: # zoom + wx = 100 + x0 = true_time_offset - old_xlims = axs[0].get_xlim() - axs[0].set_xlim( x0-wx, x0+wx) - fig.savefig('figures/11_antenna_signals_zoom.pdf') + old_xlims = axs[0].get_xlim() + axs[0].set_xlim( x0-wx, x0+wx) + fig.savefig('figures/11_antenna_signals_zoom.pdf') - # restore - 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') + # restore + axs[0].set_xlim(*old_xlims) if True: - plt.close(fig2) + plt.close(fig) - # determine correlations with arguments - lag_dt = upsampled_t[1] - upsampled_t[0] - corrs, (out1_signal, out2_template, lags) = my_correlation(upsampled_trace, template.signal) + axs2 = None + if True: # upsampled trace + upsampled_trace, upsampled_t = trace_upsampler(antenna.signal, template.t, antenna.t) - else: # downsampled template - raise NotImplementedError + if do_plots: # Show upsampled traces + 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) - lag_dt = upsampled_t[1] - upsampled_t[0] + 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') - # Determine best correlation time - 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 + fig2.savefig('figures/11_upsampled.pdf') - if not do_plots: - continue + wx = 1e2 + 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: - 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) + if True: + plt.close(fig2) - # 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)) + # determine correlations with arguments + lag_dt = upsampled_t[1] - upsampled_t[0] + corrs, (out1_signal, out2_template, lags) = my_correlation(upsampled_trace, template.signal) - # start the figure - fig, axs = plt.subplots(2, sharex=True) - ylabel_kwargs = dict( - #rotation=0, - ha='right', - va='center' - ) - axs[-1].set_xlabel("Time [ns]") + else: # downsampled template + raise NotImplementedError - offset_list = [ - [best_time_lag, dict(label=template.name, color='orange')], - [true_time_offset, dict(label='True offset', color='green')], - ] + corrs, (out1_signal, out2_template, lags) = my_downsampling_correlation(template.signal, antenna.signal, template.t, antenna.t) + lag_dt = upsampled_t[1] - upsampled_t[0] - # Signal - i=0 - axs[i].set_ylabel("Amplitude", **ylabel_kwargs) - axs[i].plot(antenna.t, antenna.signal, label=antenna.name) + # Determine best correlation time + 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 - # Plot the template - for offset_args in offset_list: - this_kwargs = offset_args[1] - offset = offset_args[0] + if not do_plots: + continue - 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 - i=1 - axs[i].set_ylabel("Correlation", **ylabel_kwargs) - axs[i].plot(lags * lag_dt, corrs) + # start the figure + fig, axs = plt.subplots(2, sharex=True) + ylabel_kwargs = dict( + #rotation=0, + ha='right', + va='center' + ) + axs[-1].set_xlabel("Time [ns]") - # 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 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)}" + offset_list = [ + [best_time_lag, dict(label=template.name, color='orange')], + [true_time_offset, dict(label='True offset', color='green')], ] - # 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 - ) + # Signal + i=0 + axs[i].set_ylabel("Amplitude", **ylabel_kwargs) + axs[i].plot(antenna.t, antenna.signal, label=antenna.name) - 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()