diff --git a/utils/python/hybrid_observables_plot_sample.py b/utils/python/hybrid_observables_plot_sample.py index a2c617233..a0baac188 100755 --- a/utils/python/hybrid_observables_plot_sample.py +++ b/utils/python/hybrid_observables_plot_sample.py @@ -27,6 +27,7 @@ from lib.gnss_sdr_conf import ( ConfigError, add_conf_argument, load_gnss_sdr_conf, + signal_pretty_name, ) from lib.plot_format import add_output_format_argument, apply_publication_style from lib.read_hybrid_observables_dump import read_hybrid_observables_dump @@ -104,6 +105,20 @@ def apply_conf_defaults(args): else: args.channels = DEFAULT_CHANNELS + args.channel_signals = channel_signal_labels(conf, args.channels) + + +def channel_signal_labels(conf, channels): + labels = [""] * channels + if conf is None: + return labels + + for signal in conf.enabled_signals: + for channel in signal.channels: + if 0 <= channel < channels: + labels[channel] = signal.signal + return labels + def first_valid_observable(gnss_observables, channels): min_tow_idx = None @@ -122,6 +137,100 @@ def first_valid_observable(gnss_observables, channels): return min_tow_idx, obs_idx +def valid_channel_data(gnss_observables, channel, observable): + valid_mask = np.array(gnss_observables["valid"][channel]) > 0 + rx_time = np.array(gnss_observables["RX_time"][channel]) + values = np.array(gnss_observables[observable][channel]) + return rx_time[valid_mask], values[valid_mask] + + +def valid_channel_samples(gnss_observables, channel, observable): + valid_mask = np.array(gnss_observables["valid"][channel]) > 0 + rx_time = np.array(gnss_observables["RX_time"][channel]) + values = np.array(gnss_observables[observable][channel]) + prns = np.array(gnss_observables["PRN"][channel]) + return rx_time[valid_mask], values[valid_mask], prns[valid_mask] + + +def valid_time_bounds(gnss_observables, channels): + valid_times = [] + for channel in range(channels): + rx_time, _ = valid_channel_data(gnss_observables, channel, "valid") + if len(rx_time) > 0: + valid_times.append(rx_time) + + if not valid_times: + raise ValueError("No valid observables found in the dump.") + + valid_times = np.concatenate(valid_times) + return valid_times.min() - 100, valid_times.max() + 100 + + +def series_by_tracked_prn(gnss_observables, channels, observable, channel_signals): + parts = {} + for channel in range(channels): + rx_time, values, prns = valid_channel_samples( + gnss_observables, channel, observable + ) + signal = channel_signals[channel] if channel < len(channel_signals) else "" + for prn in np.unique(prns): + if int(prn) == 0: + continue + selected = prns == prn + parts.setdefault((signal, int(prn)), []).append( + (rx_time[selected], values[selected]) + ) + + series = {} + for series_key, chunks in parts.items(): + rx_time = np.concatenate([chunk[0] for chunk in chunks]) + values = np.concatenate([chunk[1] for chunk in chunks]) + order = np.argsort(rx_time, kind="stable") + series[series_key] = (rx_time[order], values[order]) + return series + + +def tracked_prn_label(series_key): + signal, prn = series_key + if signal: + return f"{signal_pretty_name(signal)} PRN {prn}" + return f"PRN {prn}" + + +def plot_observable_by_tracked_prn( + gnss_observables, + args, + *, + observable, + title, + ylabel, + window_title, + output_name, + time_start, + time_end, +): + plt.figure() + plt.title(title) + series = series_by_tracked_prn( + gnss_observables, args.channels, observable, args.channel_signals + ) + for series_key in sorted(series): + rx_time, values = series[series_key] + plt.scatter( + rx_time, + values, + s=1, + label=tracked_prn_label(series_key), + ) + plt.xlim(time_start, time_end) + plt.xlabel("TOW [s]") + plt.ylabel(ylabel) + plt.grid(True) + plt.legend() + plt.gcf().canvas.manager.set_window_title(window_title) + save_figure(args.fig_path, output_name, args.show, args.output_format) + + def save_figure(fig_path, name, show, output_format): plt.tight_layout() plt.savefig(fig_path / f"{name}.{output_format}") @@ -140,74 +249,56 @@ def main(): directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) observables_file = directory / f"{base}.dat" gnss_observables = read_hybrid_observables_dump(args.channels, observables_file) - min_tow_idx, obs_idx = first_valid_observable(gnss_observables, args.channels) - time_start = gnss_observables["RX_time"][obs_idx][min_tow_idx] - 100 - time_end = gnss_observables["RX_time"][obs_idx][-1] + 100 + first_valid_observable(gnss_observables, args.channels) + time_start, time_end = valid_time_bounds(gnss_observables, args.channels) - plt.figure() - plt.title("Pseudorange") - for channel in range(args.channels): - plt.scatter( - gnss_observables["RX_time"][channel][min_tow_idx:], - gnss_observables["Pseudorange_m"][channel][min_tow_idx:], - s=1, - label=f"Channel {channel}", - ) - plt.xlim(time_start, time_end) - plt.grid(True) - plt.xlabel("TOW [s]") - plt.ylabel("Pseudorange [m]") - plt.legend() - plt.gcf().canvas.manager.set_window_title("Pseudorange.png") - save_figure(args.fig_path, "Pseudorange", args.show, args.output_format) + plot_observable_by_tracked_prn( + gnss_observables, + args, + observable="Pseudorange_m", + title="Pseudorange", + ylabel="Pseudorange [m]", + window_title="Pseudorange.png", + output_name="Pseudorange", + time_start=time_start, + time_end=time_end, + ) - plt.figure() - plt.title("Carrier Phase") - for channel in range(args.channels): - plt.scatter( - gnss_observables["RX_time"][channel][min_tow_idx:], - gnss_observables["Carrier_phase_hz"][channel][min_tow_idx:], - s=1, - label=f"Channel {channel}", - ) - plt.xlim(time_start, time_end) - plt.xlabel("TOW [s]") - plt.ylabel("Accumulated Carrier Phase [cycles]") - plt.grid(True) - plt.legend() - plt.gcf().canvas.manager.set_window_title("AccumulatedCarrierPhase.png") - save_figure(args.fig_path, "AccumulatedCarrierPhase", args.show, args.output_format) + plot_observable_by_tracked_prn( + gnss_observables, + args, + observable="Carrier_phase_hz", + title="Carrier Phase", + ylabel="Accumulated Carrier Phase [cycles]", + window_title="AccumulatedCarrierPhase.png", + output_name="AccumulatedCarrierPhase", + time_start=time_start, + time_end=time_end, + ) - plt.figure() - plt.title("Doppler Effect") - for channel in range(args.channels): - plt.scatter( - gnss_observables["RX_time"][channel][min_tow_idx:], - gnss_observables["Carrier_Doppler_hz"][channel][min_tow_idx:], - s=1, - label=f"Channel {channel}", - ) - plt.xlim(time_start, time_end) - plt.xlabel("TOW [s]") - plt.ylabel("Doppler Frequency [Hz]") - plt.grid(True) - plt.legend() - plt.gcf().canvas.manager.set_window_title("DopplerFrequency.png") - save_figure(args.fig_path, "DopplerFrequency", args.show, args.output_format) + plot_observable_by_tracked_prn( + gnss_observables, + args, + observable="Carrier_Doppler_hz", + title="Doppler Effect", + ylabel="Doppler Frequency [Hz]", + window_title="DopplerFrequency.png", + output_name="DopplerFrequency", + time_start=time_start, + time_end=time_end, + ) plt.figure() plt.title("GNSS Channels captured") for channel in range(args.channels): - label = "unknown" - for prn in gnss_observables["PRN"][channel][min_tow_idx:]: - if int(prn) != 0: - label = str(int(prn)) - break + rx_time, prns = valid_channel_data(gnss_observables, channel, "PRN") + if len(rx_time) == 0: + continue plt.scatter( - gnss_observables["RX_time"][channel][min_tow_idx:], - gnss_observables["PRN"][channel][min_tow_idx:], + rx_time, + prns, s=1, - label=f"PRN {channel} = {label}", + label=f"Channel {channel}", ) plt.xlim(time_start, time_end) plt.xlabel("TOW [s]") diff --git a/utils/python/lib/read_hybrid_observables_dump.py b/utils/python/lib/read_hybrid_observables_dump.py index 9b37df886..6b0d81be6 100644 --- a/utils/python/lib/read_hybrid_observables_dump.py +++ b/utils/python/lib/read_hybrid_observables_dump.py @@ -29,13 +29,13 @@ def read_hybrid_observables_dump(channels, filename): double_size_bytes = 8 bytes_shift = 0 - RX_time = [[] for _ in range(channels+1)] - d_TOW_at_current_symbol = [[] for _ in range(channels+1)] - Carrier_Doppler_hz = [[] for _ in range(channels+1)] - Carrier_phase_hz = [[] for _ in range(channels+1)] - Pseudorange_m = [[] for _ in range(channels+1)] - PRN = [[] for _ in range(channels+1)] - valid = [[] for _ in range(channels+1)] + RX_time = [[] for _ in range(channels)] + d_TOW_at_current_symbol = [[] for _ in range(channels)] + Carrier_Doppler_hz = [[] for _ in range(channels)] + Carrier_phase_hz = [[] for _ in range(channels)] + Pseudorange_m = [[] for _ in range(channels)] + PRN = [[] for _ in range(channels)] + valid = [[] for _ in range(channels)] f = open(filename, 'rb') if f is None: @@ -43,8 +43,7 @@ def read_hybrid_observables_dump(channels, filename): else: while True: try: - # There is an empty channel at the end (Channel-6) - for N in range(0, channels+1): + for N in range(0, channels): f.seek(bytes_shift, 0) RX_time[N].append(struct.unpack( @@ -86,19 +85,6 @@ def read_hybrid_observables_dump(channels, filename): # Reached a partial record at end of file: stop reading. break - # Delete the trailing empty channel written after the configured ones. - RX_time = [row for i, row in enumerate(RX_time) if i != channels] - d_TOW_at_current_symbol = [row for i, row in enumerate( - d_TOW_at_current_symbol)if i != channels] - Carrier_Doppler_hz = [row for i, row in enumerate( - Carrier_Doppler_hz) if i != channels] - Carrier_phase_hz = [row for i, row in enumerate( - Carrier_phase_hz) if i != channels] - Pseudorange_m = [ - row for i, row in enumerate(Pseudorange_m) if i != channels] - PRN = [row for i, row in enumerate(PRN) if i != channels] - valid = [row for i, row in enumerate(valid) if i != channels] - observables = { 'RX_time': RX_time, 'd_TOW_at_current_symbol': d_TOW_at_current_symbol,