125 lines
4.0 KiB
Python
Executable File
125 lines
4.0 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
# vim: fdm=indent ts=4
|
|
|
|
"""
|
|
Report best time offset per frequency for each antenna
|
|
"""
|
|
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from os import path
|
|
|
|
import aa_generate_beacon as beacon
|
|
from lib import figlib
|
|
|
|
if __name__ == "__main__":
|
|
import sys
|
|
import os
|
|
import matplotlib
|
|
if os.name == 'posix' and "DISPLAY" not in os.environ:
|
|
matplotlib.use('Agg')
|
|
|
|
from scriptlib import MyArgumentParser
|
|
parser = MyArgumentParser()
|
|
args = parser.parse_args()
|
|
|
|
figsize = (12,8)
|
|
|
|
fig_dir = args.fig_dir # set None to disable saving
|
|
show_plots = args.show_plots
|
|
|
|
####
|
|
fname_dir = args.data_dir
|
|
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
|
|
time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname
|
|
beacon_snr_fname = path.join(fname_dir, beacon.beacon_snr_fname)
|
|
|
|
# create fig_dir
|
|
if fig_dir:
|
|
os.makedirs(fig_dir, exist_ok=True)
|
|
|
|
# Read in antennas from file
|
|
_, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
|
|
|
|
# Read in snr info
|
|
beacon_snrs = beacon.read_snr_file(beacon_snr_fname)
|
|
snr_str = f"$\\langle SNR \\rangle$ = {beacon_snrs['mean']: .1g}"
|
|
|
|
# For now only implement using one freq_name
|
|
freq_names = antennas[0].beacon_info.keys()
|
|
if len(freq_names) > 1:
|
|
raise NotImplementedError
|
|
|
|
freq_name = next(iter(freq_names))
|
|
f_beacon = antennas[0].beacon_info[freq_name]['freq']
|
|
|
|
|
|
# TODO: redo matrix sweeping for new timing??
|
|
measured_antenna_time_shifts = {}
|
|
for i, ant in enumerate(antennas):
|
|
clock_phase_time = ant.beacon_info[freq_name]['clock_phase_mean']/(2*np.pi*f_beacon)
|
|
best_k_time = ant.beacon_info[freq_name]['best_k_time']
|
|
|
|
total_clock_time = best_k_time + clock_phase_time
|
|
measured_antenna_time_shifts[ant.name] = -1*total_clock_time
|
|
|
|
###
|
|
# Compare actual vs measured time shifts
|
|
###
|
|
actual_antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) }
|
|
|
|
N_ant = len(antennas)
|
|
|
|
if True:
|
|
# keep dataset in the same ordering
|
|
antenna_names = [int(k)-1 for k,v in actual_antenna_time_shifts.items()]
|
|
actual_time_shifts = np.array([ v for k,v in actual_antenna_time_shifts.items()])
|
|
measured_time_shifts = np.array([ measured_antenna_time_shifts[k] for k,v in actual_antenna_time_shifts.items() ])
|
|
|
|
# remove global shift
|
|
global_shift = actual_time_shifts[0] - measured_time_shifts[0]
|
|
actual_time_shifts -= global_shift
|
|
|
|
for i in range(2):
|
|
plot_residuals = i == 1
|
|
|
|
|
|
true_phases = actual_time_shifts
|
|
measured_phases = measured_time_shifts
|
|
|
|
hist_kwargs = {}
|
|
if plot_residuals:
|
|
measured_phases = measured_phases - true_phases
|
|
hist_kwargs['histtype'] = 'stepfilled'
|
|
|
|
fig = figlib.phase_comparison_figure(
|
|
measured_phases,
|
|
true_phases,
|
|
plot_residuals=plot_residuals,
|
|
f_beacon=f_beacon,
|
|
figsize=figsize,
|
|
hist_kwargs=hist_kwargs,
|
|
secondary_axis='phase',
|
|
fit_gaussian=True,
|
|
)
|
|
|
|
axs = fig.get_axes()
|
|
|
|
axs[0].legend(title=snr_str)
|
|
|
|
if plot_residuals:
|
|
fig.suptitle("Difference between Measured and Actual clock offsets")
|
|
axs[-1].set_xlabel("Antenna Time Offset Residual $\\Delta_t$ [ns]")
|
|
else:
|
|
fig.suptitle("Comparison Measured and Actual clock offset")
|
|
axs[-1].set_xlabel("Antenna Time Offset $t_c = \\left(\\frac{\\Delta\\varphi}{2\\pi} + k\\right) / f_{beac}$ [ns]")
|
|
|
|
if fig_dir:
|
|
extra_name = "comparison"
|
|
if plot_residuals:
|
|
extra_name = "residuals"
|
|
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".time.{extra_name}.pdf"))
|
|
|
|
if show_plots:
|
|
plt.show()
|