Pulse: analog_template interpol + rename snr_factor

simulations/11_pulsed_timing.py

@@ -127,20 +127,20 @@ def trace_upsampler(trace, template_t, trace_t):
 def trace_downsampler(trace, template_t, trace_t, offset):
     pass
 
-def read_time_residuals_cache(cache_fname, template_dt, antenna_dt, noise_sigma_factor, N=None):
+def read_time_residuals_cache(cache_fname, template_dt, antenna_dt, snr_sigma_factor, N=None):
     try:
         with h5py.File(cache_fname, 'r') as fp:
             pgroup = fp['time_residuals']
             pgroup2 = pgroup[f'{template_dt}_{antenna_dt}']
 
-            ds_name = str(noise_sigma_factor)
+            ds_name = str(snr_sigma_factor)
 
             ds = pgroup2[ds_name]
             if N is None:
                 return deepcopy(ds[:])
             else:
                 return deepcopy(ds[:min(N, len(ds))])
-    except KeyError:
+    except (KeyError, FileNotFoundError):
         return np.array([])
 
 def write_time_residuals_cache(cache_fname, time_residuals, template_dt, antenna_dt, noise_sigma_factor):
@@ -164,13 +164,23 @@ if __name__ == "__main__":
 
     bp_freq = (30e-3, 80e-3) # GHz
     template_dt = 5e-2 # ns
-    template_length = 500 # 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])
-    noise_factors = [1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2, 1e-1, 3e-1, 5e-1, 7e-1] # amplitude factor
+    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 = 2048 # ns
+    antenna_timelength = 1024 # ns
+
+    cut_wrong_peak_matches = True
 
     #
     # Create the template
@@ -180,7 +190,16 @@ if __name__ == "__main__":
     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)
+
+    # Interpolation Template
+    # to create an 'analog' sampled antenna
+    interp_template = Waveform(None, dt=interp_template_dt, name='Interpolation Template')
+    _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
         fig, ax = plt.subplots()
@@ -189,6 +208,8 @@ if __name__ == "__main__":
         ax.set_ylabel("Amplitude")
         ax.plot(template.t, max(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 Interpolation Template', ls='dashed')
+        ax.legend()
         fig.savefig('figures/11_template_deltapeak.pdf')
 
         if True: # show filtering equivalence samplerates
@@ -210,12 +231,13 @@ if __name__ == "__main__":
     #
     h5_cache_fname = f'11_pulsed_timing.hdf5'
 
-    time_accuracies = np.zeros(len(noise_factors))
+    time_accuracies = np.zeros(len(snr_factors))
-    for k, noise_sigma_factor in tqdm(enumerate(noise_factors)):
+    for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
-        print() #separating tqdm
-
         # Read in cached time residuals
-        cached_time_residuals = read_time_residuals_cache(h5_cache_fname, template.dt, antenna_dt, noise_sigma_factor)
+        if True:
+            cached_time_residuals = read_time_residuals_cache(h5_cache_fname, template.dt, antenna_dt, snr_sigma_factor)
+        else:
+            cached_time_residuals = np.array([])
 
         #
         # Find difference between true and templated times
@@ -243,7 +265,8 @@ if __name__ == "__main__":
 
                 antenna.t = util.sampled_time(1/antenna.dt, start=0, end=antenna_timelength)
 
-                antenna.signal = interp1d_template(antenna.t - antenna.peak_time)
+                # Sample the interpolation template
+                antenna.signal = interp_template.interpolate(antenna.t - antenna.peak_time)
 
                 antenna.peak_sample = antenna.peak_time/antenna.dt
                 antenna_true_signal = antenna.signal
@@ -254,7 +277,7 @@ if __name__ == "__main__":
                 antenna.signal *= -1
 
             ## Add noise
-            noise_amplitude = max(template.signal) * noise_sigma_factor
+            noise_amplitude = max(template.signal) * 1/snr_sigma_factor
             noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal))
             filtered_noise = antenna_bp(noise_realisation, *bp_freq, antenna.dt)
 
@@ -323,16 +346,18 @@ if __name__ == "__main__":
                 lag_dt = upsampled_t[1] - upsampled_t[0]
                 corrs, (out1_signal, out2_template, lags) = my_correlation(upsampled_trace, template.signal)
 
-            else: # downsampled template
-                raise NotImplementedError
-
-                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]
-
                 # Determine best correlation time
                 idx = np.argmax(abs(corrs))
                 best_sample_lag = lags[idx]
                 best_time_lag = best_sample_lag * lag_dt
+
+            else: # downsampled template
+                raise NotImplementedError
+
+                corrs, (_, _, lags) = my_downsampling_correlation(antenna.signal, antenna.t, template.signal, template.t)
+                lag_dt = upsampled_t[1] - upsampled_t[0]
+
+            # Calculate the time residual
             time_residuals[j] = best_time_lag - true_time_offset
 
             if not do_plots:
@@ -408,13 +433,16 @@ if __name__ == "__main__":
                 if True:
                     plt.close(fig)
 
+        print()# separating tqdm
         print()# separating tqdm
         # Were new time residuals calculated?
         # Add them to the cache file
         if len(time_residuals) > 1:
             # merge cached and calculated time residuals
             time_residuals = np.concatenate((cached_time_residuals, time_residuals), axis=None)
-            write_time_residuals_cache(h5_cache_fname, time_residuals, template_dt, antenna_dt, noise_sigma_factor)
+
+            if True: # write the cache
+                write_time_residuals_cache(h5_cache_fname, time_residuals, template_dt, antenna_dt, snr_sigma_factor)
         else:
             time_residuals = cached_time_residuals
 
@@ -498,13 +526,17 @@ if __name__ == "__main__":
             ax.set_yscale('log')
 
         # plot the values
-        ax.plot(1/np.asarray(noise_factors), time_accuracies, ls='none', marker='o')
+        ax.plot(np.asarray(snr_factors), time_accuracies, ls='none', marker='o')
+
+        if True: # limit y-axis to 1e0
+            ax.set_ylim([None, 1e1])
 
 
         # Set horizontal line at 1 ns
         if True:
             ax.axhline(1, ls='--', alpha=0.8, color='g')
             ax.grid()
+            ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color='b')
 
         fig.tight_layout()
         fig.savefig(f"figures/11_time_res_vs_snr_tdt{template_dt:0.1e}.pdf")