From 488a4446f20c3a30f1d7df29630652a11dc8da73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jakub=20K=C3=A1kona?= Date: Sun, 6 Sep 2026 11:46:16 +0200 Subject: [PATCH] Detect firmware event-buffer truncation and plot count rate, not raw sums MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The AIRDOS03 firmware caps above-threshold events at MAX_EVENTS per interval; anything past that is silently dropped while the true count (evCount in $STOP) keeps growing. Nothing surfaced this to the user. Separately, since the firmware can now flush early when that buffer fills (see the corresponding AIRDOS03 firmware change), interval length is no longer fixed, making raw per-record totals hard to compare. - Warn (file parser and live UART thread) whenever a $STOP's evCount exceeds the number of $E lines actually received, i.e. the firmware's per-interval event buffer was exceeded and the spectrum above THRESHOLD is undercounted for that interval. - Parse captStart/stopSystime from $START/$STOP and derive each record's real duration: prefer the delta between consecutive GNSS-synced timestamps (or host wall-clock in the live case), falling back to the raw TCNT1 tick pair (128 µs/tick) when no absolute time is available yet. Persisted through save_as()/NpzLogParser so reopened sessions keep it. - Evolution plot now shows counts/sec ("Count rate", cps) whenever a duration is known, since a shortened (early-flushed) interval and a full one are no longer directly comparable as raw sums. Falls back to the previous "Total count per exposition" raw-count display for sources with no duration info (older logs/NPZ archives). Verified against the existing test suite plus manual smoke tests (synthetic multi-run logs exercising GNSS-diff, tick-fallback and NPZ round-trip paths); no display-side (PyQt) smoke test since this environment has no display, only offscreen import/construction checks. --- dosview/__init__.py | 111 ++++++++++++++++++++++++++++++++++++++++---- dosview/parsers.py | 78 ++++++++++++++++++++++++++++++- 2 files changed, 178 insertions(+), 11 deletions(-) diff --git a/dosview/__init__.py b/dosview/__init__.py index a2ccf02..f622f27 100644 --- a/dosview/__init__.py +++ b/dosview/__init__.py @@ -44,6 +44,7 @@ OldLogParser, get_parser_for_file, parse_file, + TICK_SECONDS, ) from .eeprom_widget import EepromManagerWidget from .rtc_widget import RTCManagerWidget @@ -144,6 +145,23 @@ def __init__(self, parent=None, file_path=None): self._updating = False self._evo_region = None self._spec_region = None + self._durations = None + self._cur_display = np.array([]) + + def _extract_durations(self, data): + """Per-record interval length (seconds), when the source provides one. + + Present for AIRDOS v2 logs/live sessions (fw 1.3+, see parsers.py); + absent for older logs, in which case the evolution plot falls back + to raw counts per record. + """ + dur = data[6] if len(data) > 6 else None + if dur is None: + return None + dur = np.asarray(dur, dtype=float) + if dur.shape[0] != self._sums.shape[0]: + return None + return dur def plot(self, data): start_time = time.time() @@ -163,6 +181,7 @@ def plot(self, data): self._sums = np.asarray(data[1], dtype=float) self._hist = np.asarray(data[2], dtype=float) self._channels = np.arange(len(self._hist), dtype=float) + self._durations = self._extract_durations(data) sm = data[5] if len(data) > 5 else None if sm is not None and hasattr(sm, "ndim") and sm.ndim == 2 and sm.shape[0] >= 1 and sm.shape[1] > 0: self._spectral_matrix = np.asarray(sm, dtype=float) @@ -176,7 +195,15 @@ def plot(self, data): 'bottom': SymLogAxisItem(orientation='bottom')}) self.plot_evolution.showGrid(x=True, y=True) - self.plot_evolution.setLabel("left", "Total count per exposition", units="#") + # Per-record intervals are no longer all the same length once a high-rate + # burst triggers an early firmware flush (see main.cpp flushDataOut()), so + # raw counts per record stop being comparable across records — plot a + # rate instead whenever a duration is known. Falls back to raw counts + # for sources with no duration info (older logs/firmware). + if self._durations is not None: + self.plot_evolution.setLabel("left", "Count rate", units="cps") + else: + self.plot_evolution.setLabel("left", "Total count per exposition", units="#") self.plot_evolution.setLabel("bottom", "Time", units="min") self._curve_evolution = self.plot_evolution.plot([], [], @@ -240,8 +267,13 @@ def _on_mouse_moved(self, evt): idx = int(np.argmin(np.abs(self._time_axis_min - x))) t = float(self._time_axis_min[idx]) count = float(self._cur_sums[idx]) - px, py = t, float(symlog(count)) - text = f"t = {t:.2f} min\ncount = {count:.0f}" + value = float(self._cur_display[idx]) + px, py = t, float(symlog(value)) + if self._durations is not None and idx < len(self._durations): + duration = float(self._durations[idx]) + text = f"t = {t:.2f} min\nrate = {value:.2f} cps\ncount = {count:.0f} in {duration:.2f} s" + else: + text = f"t = {t:.2f} min\ncount = {count:.0f}" else: if not len(self._channels): continue @@ -257,13 +289,25 @@ def _on_mouse_moved(self, evt): item.setVisible(True) def _set_evolution_curve(self, sums): - """Draw the evolution curve (and rolling average) in symlog-y.""" + """Draw the evolution curve (and rolling average) in symlog-y. + + Plotted as counts/sec when per-record durations are known, since + records can now be shorter than the nominal interval (early flush on + a full event buffer — see flushDataOut() in firmware) and raw sums + of different-length records aren't directly comparable. Falls back + to raw counts per record otherwise. + """ sums = np.asarray(sums, dtype=float) - self._cur_sums = sums # currently displayed values (for the hover tooltip) - self._curve_evolution.setData(self._time_axis_min, symlog(sums)) + self._cur_sums = sums # raw per-record counts (for the hover tooltip) + if self._durations is not None and self._durations.shape[0] == sums.shape[0]: + display = sums / self._durations + else: + display = sums + self._cur_display = display # currently displayed values (for the hover tooltip) + self._curve_evolution.setData(self._time_axis_min, symlog(display)) window_size = self.WINDOW_SIZE - if len(sums) >= window_size: - rolling_avg = np.convolve(sums, np.ones(window_size) / window_size, mode='valid') + if len(display) >= window_size: + rolling_avg = np.convolve(display, np.ones(window_size) / window_size, mode='valid') self._curve_rolling_avg.setData(self._time_axis_min[window_size - 1:], symlog(rolling_avg)) else: self._curve_rolling_avg.setData([], []) @@ -393,6 +437,7 @@ def update_data(self, data): self._sums = np.asarray(data[1], dtype=float) self._hist = np.asarray(data[2], dtype=float) self._channels = np.arange(len(self._hist), dtype=float) + self._durations = self._extract_durations(data) sm = data[5] if len(data) > 5 else None if sm is not None and hasattr(sm, "ndim") and sm.ndim == 2 and sm.shape[0] >= 1 and sm.shape[1] > 0: self._spectral_matrix = np.asarray(sm, dtype=float) @@ -778,6 +823,7 @@ def run(self): hist = np.zeros(65536, dtype=int) time_axis = _GrowableArray(dtype=float) sums = _GrowableArray(dtype=int) + durations = _GrowableArray(dtype=float) spectral_records = _GrowableMatrix(width=len(hist), dtype=int) metadata = {"log_runs_count": 0, "log_device_info": {}} fmt = None # 'old' or 'v2' @@ -788,6 +834,8 @@ def run(self): # v2-specific inter-record state current_hist = None + current_cap_start = 0 + prev_abs_time = None current_counts = 0 def _resolve_timestamp(parts): @@ -950,6 +998,10 @@ def _build_telemetry(env_recs): elif msg == "$START": current_hist = np.zeros_like(hist) current_counts = 0 + try: + current_cap_start = int(parts[2]) if len(parts) > 2 else 0 + except ValueError: + current_cap_start = 0 elif msg == "$E" and current_hist is not None and len(parts) >= 3: try: ch = int(parts[2]) @@ -961,6 +1013,20 @@ def _build_telemetry(env_recs): logger.warning("Dropped malformed $E record: %r", line) elif msg == "$STOP" and current_hist is not None: try: + if len(parts) > 4: + try: + ev_count = int(parts[4]) + except ValueError: + ev_count = None + if ev_count is not None and ev_count > current_counts: + logger.warning( + "Device counted %d above-threshold events " + "but only %d were transmitted — the " + "firmware's per-interval event buffer was " + "exceeded, so the spectrum is undercounted " + "above THRESHOLD for this interval.", + ev_count, current_counts, + ) for idx, val in enumerate(parts[5:]): try: current_hist[idx] += int(val) @@ -973,13 +1039,38 @@ def _build_telemetry(env_recs): # fall back to a non-time field such as the run # index, which would silently compress the x-axis. t = _resolve_timestamp(parts) + + # Interval duration: _resolve_timestamp already gives a + # monotonically increasing clock (device time when + # GNSS-synced, host wall-clock otherwise), so a plain + # delta against the previous record works for the + # normal case. Only fall back to the free-running + # TCNT1 ticks in captStart/stopSystime — and finally to + # the old fixed cadence — for the first record, or if + # a GNSS sync/unsync transition made time run backwards. + if prev_abs_time is not None and t > prev_abs_time: + duration = t - prev_abs_time + else: + try: + stop_ticks = int(parts[3]) if len(parts) > 3 else None + except ValueError: + stop_ticks = None + if stop_ticks is not None: + ticks_elapsed = (stop_ticks - current_cap_start) & 0xFFFF + duration = max(ticks_elapsed, 1) * TICK_SECONDS + else: + duration = 10.0 + prev_abs_time = t + durations.append(duration) + spectral_records.append(current_hist) hist += current_hist time_axis.append(t) sums.append(int(current_hist.sum())) self.dataUpdated.emit( [time_axis.view(), sums.view(), hist.copy(), metadata, - _build_telemetry(env_records), spectral_records.view()] + _build_telemetry(env_records), spectral_records.view(), + durations.view()] ) except (ValueError, IndexError): dropped_records += 1 @@ -1840,6 +1931,8 @@ def save_as(self): arrays[f"telemetry_value_{key}"] = v if len(data) > 5 and data[5] is not None and hasattr(data[5], "shape") and data[5].ndim == 2: arrays["spectral_matrix"] = data[5] + if len(data) > 6 and data[6] is not None: + arrays["durations"] = data[6] np.savez_compressed(path, **arrays) diff --git a/dosview/parsers.py b/dosview/parsers.py index 99dc864..7891c75 100644 --- a/dosview/parsers.py +++ b/dosview/parsers.py @@ -8,11 +8,40 @@ from typing import List, Sequence, Tuple +import logging import time from pathlib import Path import numpy as np +logger = logging.getLogger(__name__) + +# Seconds per TCNT1 tick on the AIRDOS03 AVR firmware: Timer1 runs off the +# 8 MHz clock with a /1024 prescaler (see fw/*/src/main.cpp, TCCR1B). Used to +# turn the raw captStart/stopSystime tick counts in $START/$STOP into a +# duration when the device has no GNSS fix (so tm.tm_s100 is always 0). +TICK_SECONDS = 128e-6 + +# TCNT1 is 16-bit, so tick-based durations are only unambiguous up to one +# full wrap of the counter. +MAX_TICK_DURATION_S = 65536 * TICK_SECONDS + + +def _fill_invalid_durations(durations: np.ndarray) -> np.ndarray: + """Replace non-positive/non-finite duration entries with a sane fallback. + + Individual $STOP records can fail to yield a duration (malformed line). + Rather than let a stray zero/NaN divide-by-zero corrupt a counts/sec + plot, backfill with the median of the valid durations in this log (or + 1 s, i.e. "treat as raw counts", if none are valid at all). + """ + valid = np.isfinite(durations) & (durations > 0) + if valid.all(): + return durations + durations = durations.copy() + durations[~valid] = np.median(durations[valid]) if valid.any() else 1.0 + return durations + class BaseLogParser: """Base parser class.""" @@ -60,10 +89,13 @@ def parse(self): total_counts = 0 sums: List[int] = [] time_axis: List[float] = [] + durations: List[float] = [] spectral_records: List[np.ndarray] = [] inside_run = False current_hist = None current_counts = 0 + current_cap_start = 0 + prev_abs_time: float | None = None device_type = "unknown" env_records: List[Tuple[float, ...]] = [] batt_records: List[Tuple[float, ...]] = [] @@ -102,6 +134,10 @@ def parse(self): inside_run = True current_hist = np.zeros_like(hist) current_counts = 0 + try: + current_cap_start = int(parts[2]) if len(parts) > 2 else 0 + except ValueError: + current_cap_start = 0 case "$E": if inside_run and len(parts) >= 3: channel = int(parts[2]) @@ -110,6 +146,20 @@ def parse(self): current_counts += 1 case "$STOP": if inside_run: + if len(parts) > 4: + try: + ev_count = int(parts[4]) + except ValueError: + ev_count = None + if ev_count is not None and ev_count > current_counts: + logger.warning( + "Run %d: device counted %d above-threshold " + "events but only %d were transmitted — the " + "firmware's per-interval event buffer was " + "exceeded, so the spectrum is undercounted " + "above THRESHOLD for this interval.", + metadata["log_runs_count"], ev_count, current_counts, + ) if len(parts) > 5: for idx, val in enumerate(parts[5:]): try: @@ -120,7 +170,27 @@ def parse(self): hist += current_hist total_counts += current_counts sums.append(current_counts) - time_axis.append(float(parts[2])) + cur_abs_time = float(parts[2]) + time_axis.append(cur_abs_time) + + # Interval duration: prefer the GNSS-synced wall-clock + # delta between consecutive $STOP timestamps (works for + # any length); fall back to the free-running TCNT1 tick + # count in captStart/stopSystime (works without GNSS, + # but only unambiguous up to one 16-bit timer wrap). + if prev_abs_time and cur_abs_time > prev_abs_time: + durations.append(cur_abs_time - prev_abs_time) + else: + try: + stop_ticks = int(parts[3]) if len(parts) > 3 else None + except ValueError: + stop_ticks = None + if stop_ticks is not None: + ticks_elapsed = (stop_ticks - current_cap_start) & 0xFFFF + durations.append(max(ticks_elapsed, 1) * TICK_SECONDS) + else: + durations.append(float("nan")) + prev_abs_time = cur_abs_time if cur_abs_time else None inside_run = False current_hist = None case "$ENV": @@ -179,8 +249,9 @@ def parse(self): telemetry["temperature"] = (ba[:, 0], ba[:, 5]) spectral_matrix = np.array(spectral_records) if spectral_records else np.zeros((0, hist.shape[0]), dtype=int) + durations_arr = _fill_invalid_durations(np.array(durations, dtype=float)) if durations else np.zeros(0) print("Parsed AIRDOS v2 format in", time.time() - start_time, "s") - return [np.array(time_axis), np.array(sums), hist, metadata, telemetry, spectral_matrix] + return [np.array(time_axis), np.array(sums), hist, metadata, telemetry, spectral_matrix, durations_arr] # Backwards-compatible alias @@ -319,6 +390,8 @@ def parse(self): spectral_matrix = npz["spectral_matrix"] else: spectral_matrix = np.zeros((0, hist.shape[0]), dtype=int) + if "durations" in npz.files: + return [time_axis, sums, hist, metadata, telemetry, spectral_matrix, npz["durations"]] return [time_axis, sums, hist, metadata, telemetry, spectral_matrix] @@ -345,4 +418,5 @@ def parse_file(file_path: str | Path): "NpzLogParser", "get_parser_for_file", "parse_file", + "TICK_SECONDS", ]