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", ]