From 465b78c535f20be0b6bb92938cffdc67307a75ad Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Thu, 13 Apr 2023 17:09:30 +0200 Subject: [PATCH] ZH: Determine SNR for Airshower vs Noise --- .../aa_generate_beacon.py | 10 +- .../ac_show_signal_to_noise.py | 92 +++++++++++++++++-- airshower_beacon_simulation/lib/snr.py | 25 ++++- 3 files changed, 113 insertions(+), 14 deletions(-) diff --git a/airshower_beacon_simulation/aa_generate_beacon.py b/airshower_beacon_simulation/aa_generate_beacon.py index d1ae538..0422661 100755 --- a/airshower_beacon_simulation/aa_generate_beacon.py +++ b/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) diff --git a/airshower_beacon_simulation/ac_show_signal_to_noise.py b/airshower_beacon_simulation/ac_show_signal_to_noise.py index ed98c35..6ca8a70 100755 --- a/airshower_beacon_simulation/ac_show_signal_to_noise.py +++ b/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") diff --git a/airshower_beacon_simulation/lib/snr.py b/airshower_beacon_simulation/lib/snr.py index 4044e2c..002c6a1 100644 --- a/airshower_beacon_simulation/lib/snr.py +++ b/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