Pulse: move templating to function

simulations/11_pulsed_timing.py

@@ -160,6 +160,18 @@ def write_time_residuals_cache(cache_fname, time_residuals, template_dt, antenna
 
         ds = pgroup2.create_dataset(ds_name, (len(time_residuals)), dtype='f', data=time_residuals, maxshape=(None))
 
+def create_template(dt=1, timelength=1, bp_freq=(0, np.inf), name=None, normalise=False):
+    template = Waveform(None, dt=dt, name=name)
+    _deltapeak = util.deltapeak(timelength=timelength, samplerate=1/dt, offset=0)
+    template.signal = antenna_bp(_deltapeak[0], *bp_freq, dt)
+    template.peak_sample = _deltapeak[1]
+    template.peak_time = template.dt * template.peak_sample
+
+    if normalise:
+        template.signal /= max(template.signal)
+
+    return template, _deltapeak
+
 if __name__ == "__main__":
     import os
     import matplotlib
@@ -186,36 +198,26 @@ if __name__ == "__main__":
     antenna_timelength = 1024 # ns
 
     cut_wrong_peak_matches = True
+    normalise_noise = False
+
+    h5_cache_fname = f'11_pulsed_timing.hdf5'
 
     #
-    # Create the template
-    #
-    template = Waveform(None, dt=template_dt, name='Template')
-    _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
-
     # Interpolation Template
     # to create an 'analog' sampled antenna
-    interp_template = Waveform(None, dt=interp_template_dt, name='Interpolation Template')
+    interp_template, _deltapeak = create_template(dt=interp_template_dt, timelength=template_length, bp_freq=bp_freq, name='Interpolation Template', normalise=True)
-    _interp_deltapeak = util.deltapeak(timelength=template_length, samplerate=1/interp_template.dt, offset=0)
-    interp_template.signal = antenna_bp(_interp_deltapeak[0], *bp_freq, interp_template.dt)
-    interp_template.peak_sample = _interp_deltapeak[1]
-    interp_template.peak_time = interp_template.dt * interp_template.peak_sample
-    interp_template.signal *= max(template.signal)/max(interp_template.signal)
-    interp_template.interpolate = interpolate.interp1d(interp_template.t, interp_template.signal, assume_sorted=True, fill_value=0, bounds_error=False)
 
-    if True: # show template
+    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(template.t, max(template.signal)*_deltapeak[0], label='Impulse Template')
+        ax.plot(interp_template.t, max(interp_template.signal)*_deltapeak[0], label='Impulse Template')
-        ax.plot(template.t, template.signal, label='Filtered Template')
+        ax.plot(interp_template.t, interp_template.signal, label='Filtered Template')
-        ax.plot(interp_template.t, interp_template.signal, label='Filtered Interpolation Template', ls='dashed')
         ax.legend()
-        fig.savefig('figures/11_template_deltapeak.pdf')
+        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)
@@ -226,16 +228,20 @@ if __name__ == "__main__":
             ax.plot(_time, _bandpassed, label='Filtered Antenna')
 
             ax.legend()
-            fig.savefig('figures/11_template_deltapeak+antenna.pdf')
+            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
     #
-    h5_cache_fname = f'11_pulsed_timing.hdf5'
-
     time_accuracies = np.zeros(len(snr_factors))
     mask_counts = np.zeros(len(snr_factors))
     for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
@@ -256,7 +262,7 @@ if __name__ == "__main__":
             ## place the deltapeak signal at a random location
             antenna = Waveform(None, dt=antenna_dt, name='Signal')
 
-            if False: # Create antenna trace without 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.peak_sample = antenna_peak_sample
@@ -265,7 +271,7 @@ if __name__ == "__main__":
                 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 interpolation 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
 
@@ -284,7 +290,7 @@ if __name__ == "__main__":
 
             ## Add noise
             noise_amplitude = max(template.signal) * 1/snr_sigma_factor
-            noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal), normalise=False)
+            noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal), normalise=normalise_noise)
             filtered_noise = antenna_bp(noise_realisation, *bp_freq, antenna.dt)
 
             antenna.signal += filtered_noise
@@ -569,7 +575,7 @@ if __name__ == "__main__":
 
             l = ax.plot(snr_factors[mask], time_accuracies[mask], **kwargs)
 
-        if True: # limit y-axis to 1e0
+        if True: # limit y-axis to 1e1
             ax.set_ylim([None, 1e1])