ZH: Determine SNR for Airshower vs Noise

airshower_beacon_simulation/aa_generate_beacon.py

@@ -18,7 +18,8 @@ import lib
 # {{{ vim marker
 tx_fname = 'tx.json'
 antennas_fname = 'antennas.hdf5'
-beacon_snr_fname = 'snr.json'
+beacon_snr_fname = 'beacon_snr.json'
+airshower_snr_fname = 'airshower_snr.json'
 c_light = lib.c_light
 
 def read_antenna_clock_repair_offsets(antennas, mode='all', freq_name=None):
@@ -418,6 +419,12 @@ if __name__ == "__main__":
         # Collect all data to be saved (with the first 3 values the E fields)
         traces = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon, noise_realisation])
 
+        append_antenna_hdf5( antennas_fname, antenna, traces, name='original_traces', prepend_time=True)
+
+        E_AxB = [np.dot(ev.uAxB,[ex,ey,ez]) for ex,ey,ez in zip(traces[0], traces[1], traces[2])]
+        t_AxB = antenna.t
+        append_antenna_hdf5( antennas_fname, antenna, [t_AxB, E_AxB, t_AxB, t_AxB], name='original_E_AxB', prepend_time=False)# Note the 4 element list containing dummys at idx 2 and 3
+
         # add beacon and noise to relevant polarisations
         for j, amp in enumerate(beacon_amplitudes):
             traces[j] = traces[j] + amp*beacon + noise_realisation
@@ -433,7 +440,6 @@ if __name__ == "__main__":
 
         # Save filtered E field in E_AxB
         E_AxB = [np.dot(ev.uAxB,[ex,ey,ez]) for ex,ey,ez in zip(traces[0], traces[1], traces[2])]
-        t_AxB = antenna.t
 
         append_antenna_hdf5( antennas_fname, antenna, [t_AxB, E_AxB], name='E_AxB', prepend_time=False)
 

airshower_beacon_simulation/ac_show_signal_to_noise.py

@@ -43,6 +43,7 @@ if __name__ == "__main__":
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
     tx_fname = path.join(fname_dir, beacon.tx_fname)
     beacon_snr_fname = path.join(fname_dir, beacon.beacon_snr_fname)
+    airshower_snr_fname = path.join(fname_dir, beacon.airshower_snr_fname)
 
     # create fig_dir
     if fig_dir:
@@ -52,6 +53,17 @@ if __name__ == "__main__":
     f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname, traces_key='filtered_traces')
     _, __, txdata = beacon.read_tx_file(tx_fname)
 
+    # Read zeropadded traces
+    _, __, signal_antennas = beacon.read_beacon_hdf5(antennas_fname, traces_key='original_E_AxB', read_AxB=False )
+    # !!HACK!! Repack traces in signal_antennas to antennas
+    for i, ant in enumerate(signal_antennas):
+        if antennas[i].name != ant.name:
+            print("Error!")
+            import sys
+            sys.exit()
+
+        antennas[i].orig_E_AxB = ant.Ex
+
     # general properties
     dt = antennas[0].t[1] - antennas[0].t[0] # ns
     beacon_pb = lib.passband(f_beacon, None) # GHz
@@ -76,7 +88,7 @@ if __name__ == "__main__":
         ant = antennas[0]
 
         fig, ax = plt.subplots(figsize=figsize)
-        _debug_snrs = lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), myfilter(ant.noise), samplerate=1/dt, signal_band=beacon_pb, noise_band=noise_pb, debug_ax=ax)
+        _debug_snrs = lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), myfilter(ant.noise), samplerate=1/dt, signal_band=beacon_pb, noise_band=noise_pb, debug_ax=ax, mode='sine')
 
         ax.legend(title="$\\langle SNR \\rangle$ = {: .1e}".format(np.mean(_debug_snrs)))
 
@@ -87,13 +99,14 @@ if __name__ == "__main__":
         ax.set_xlim(low_x, high_x)
 
         if fig_dir:
-            fig.savefig(path.join(fig_dir, path.basename(__file__) + f".debug_plot.pdf"))
+            fig.savefig(path.join(fig_dir, path.basename(__file__) + f".beacon_vs_noise_snr.debug_plot.pdf"))
+
     ##
     ## Beacon vs Noise SNR
     ##
     if True:
         N_samples = len(antennas[0].beacon)
-        beacon_snrs = [ lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), myfilter(ant.noise), samplerate=1/dt, signal_band=beacon_pb, noise_band=noise_pb) for i, ant in enumerate(antennas) ]
+        beacon_snrs = [ lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), myfilter(ant.noise), samplerate=1/dt, signal_band=beacon_pb, noise_band=noise_pb, mode='sine') for i, ant in enumerate(antennas) ]
 
         # write mean and std to file
         beacon.write_snr_file(beacon_snr_fname, beacon_snrs)
@@ -111,7 +124,7 @@ if __name__ == "__main__":
     ## Beacon vs Total SNR
     ##
     if True:
-        beacon_snrs = [ lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), ant.E_AxB, samplerate=1/dt, signal_band=beacon_pb, noise_band=pb) for ant in antennas ]
+        beacon_snrs = [ lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), ant.E_AxB, samplerate=1/dt, signal_band=beacon_pb, noise_band=pb, mode='sine') for ant in antennas ]
 
         fig, ax = plt.subplots(figsize=figsize)
         ax.set_title("Maximum Beacon/Total SNR")
@@ -122,26 +135,89 @@ if __name__ == "__main__":
         if fig_dir:
             fig.savefig(path.join(fig_dir, path.basename(__file__) + f".beacon_vs_total_snr.pdf"))
 
+    ##
+    ## Debug plot of Signal vs Noise SNR
+    ##
+    if True:
+        ant = antennas[0]
+
+        fig, ax = plt.subplots(figsize=figsize)
+        _debug_snrs = lib.signal_to_noise(myfilter(ant.orig_E_AxB), myfilter(ant.noise), samplerate=1/dt, debug_ax=ax, mode='pulse')
+
+        ax.legend(title="$\\langle SNR \\rangle$ = {: .1e}".format(np.mean(_debug_snrs)))
+
+        ax.set_title("Signal (max amp) and Noise (rms)")
+        ax.set_xlabel("Samples")
+        ax.set_ylabel("Amplitude")
+
+        if fig_dir:
+            fig.savefig(path.join(fig_dir, path.basename(__file__) + f".airshower_vs_noise_snr.debug_plot.pdf"))
+
     ##
     ## Signal vs Noise SNR
     ##
     if True:
-        beacon_snrs = [ lib.signal_to_noise(myfilter(ant.E_AxB - beacon_amp*ant.beacon), myfilter(ant.noise), samplerate=1/dt, signal_band=beacon_pb, noise_band=pb) for ant in antennas ]
+        airshower_snrs = [ lib.signal_to_noise(myfilter(ant.orig_E_AxB), myfilter(ant.noise), samplerate=1/dt, mode='pulse') for ant in antennas ]
+
+        # write mean and std to file
+        beacon.write_snr_file(airshower_snr_fname, airshower_snrs)
 
         fig, ax = plt.subplots(figsize=figsize)
         ax.set_title("Maximum Airshower/Noise SNR")
         ax.set_xlabel("Antenna no.")
         ax.set_ylabel("SNR")
-        ax.plot([ int(ant.name) for ant in antennas], beacon_snrs, 'o', ls='none')
+        ax.plot([ int(ant.name) for ant in antennas], airshower_snrs, 'o', ls='none')
 
         if fig_dir:
             fig.savefig(path.join(fig_dir, path.basename(__file__) + f".airshower_vs_noise_snr.pdf"))
 
     ##
-    ## Airshower signal vs Noise SNR
+    ## Debug plot of Signal vs Beacon SNR
     ##
     if True:
-        shower_snrs = [ lib.signal_to_noise(ant.E_AxB, myfilter(ant.noise), samplerate=1/dt, signal_band=pb, noise_band=pb) for ant in antennas ]
+        ant = antennas[0]
+
+        fig, ax = plt.subplots(figsize=figsize)
+        if False: #indirect SNR max_amp(signal) vs max_amp(beacon)
+            _debug_snrs_E_AxB = lib.signal_to_noise(myfilter(ant.orig_E_AxB), myfilter(ant.noise), samplerate=1/dt, debug_ax=ax, mode='pulse')
+            _debug_snrs_sine = lib.signal_to_noise(myfilter(beacon_amp*ant.beacon), myfilter(ant.noise), samplerate=1/dt, debug_ax=ax, mode='pulse')
+
+            _debug_snrs = _debug_snrs_E_AxB / _debug_snrs_sine
+        else: # direct max_amp(signal) vs rms(beacon)
+            _debug_snrs = lib.signal_to_noise(myfilter(ant.orig_E_AxB), myfilter(beacon_amp*ant.beacon), samplerate=1/dt, debug_ax=ax, mode='pulse')
+
+        ax.legend(title="$\\langle SNR \\rangle$ = {: .1e}".format(np.mean(_debug_snrs)))
+
+        ax.set_title("Signal (max amp), Beacon (max amp) and Noise (rms)")
+        ax.set_xlabel("Samples")
+        ax.set_ylabel("Amplitude")
+
+        if fig_dir:
+            fig.savefig(path.join(fig_dir, path.basename(__file__) + f".airshower_vs_beacon_snr.debug_plot.pdf"))
+
+
+    ##
+    ## Signal vs Beacon SNR
+    ##
+    if True:
+        shower_beacon_snrs = [ lib.signal_to_noise(myfilter(ant.orig_E_AxB), myfilter(beacon_amp*ant.beacon), samplerate=1/dt, mode='pulse') for ant in antennas ]
+
+        fig, ax = plt.subplots(figsize=figsize)
+        ax.set_title("Maximum Airshower/Beacon RMS SNR")
+        ax.set_xlabel("Antenna no.")
+        ax.set_ylabel("SNR")
+        ax.plot([ int(ant.name) for ant in antennas], beacon_snrs, 'o', ls='none')
+
+        if fig_dir:
+            fig.savefig(path.join(fig_dir, path.basename(__file__) + f".airshower_vs_beacon_snr.pdf"))
+
+
+
+    ##
+    ## Total signal vs Noise SNR
+    ##
+    if True:
+        shower_snrs = [ lib.signal_to_noise(ant.E_AxB, myfilter(ant.noise), samplerate=1/dt, mode='pulse') for ant in antennas ]
 
         fig, ax = plt.subplots(figsize=figsize)
         ax.set_title("Total (Signal+Beacon+Noise)/Noise SNR")

airshower_beacon_simulation/lib/snr.py

@@ -80,7 +80,7 @@ def bandpower(samples, samplerate=1, band=passband(), normalise_bandsize=True, d
 
     return power
 
-def signal_to_noise(samples, noise, samplerate=1, signal_band=passband(), noise_band=None, debug_ax=False):
+def signal_to_noise(samples, noise, samplerate=1, signal_band=passband(), noise_band=None, debug_ax=False, mode='sine'):
     if noise_band is None:
         noise_band = signal_band
 
@@ -90,8 +90,25 @@ def signal_to_noise(samples, noise, samplerate=1, signal_band=passband(), noise_
     if debug_ax is True:
         debug_ax = plt.gca()
 
-    noise_power = bandpower(noise, samplerate, noise_band, debug_ax=debug_ax)
+    if mode == 'sine':
+        noise_power = bandpower(noise, samplerate, noise_band, debug_ax=debug_ax)
+        noise_amplitude = np.sqrt(noise_power)
 
-    signal_power = bandpower(samples, samplerate, signal_band, debug_ax=debug_ax)
+        signal_power = bandpower(samples, samplerate, signal_band, debug_ax=debug_ax)
+        signal_amplitude = np.sqrt(signal_power)
 
-    return (signal_power/noise_power)**0.5
+    elif mode == 'pulse':
+        noise_amplitude = np.sqrt(np.mean(noise**2))
+
+        signal_amplitude = max(np.abs(samples))
+
+        if debug_ax:
+            l1 = debug_ax.plot(noise, alpha=0.5)
+            debug_ax.axhline(noise_amplitude, alpha=0.9, color=l1[0].get_color())
+
+            l2 = debug_ax.plot(samples, alpha=0.5)
+            debug_ax.axhline(signal_amplitude, alpha=0.9, color=l2[0].get_color())
+    else:
+        raise NotImplementedError("mode not in ['sine', 'pulse']")
+
+    return signal_amplitude/noise_amplitude