seismo-relay/minimateplus/histogram_codec.py

"""
histogram_codec.py — decoder for MiniMate Plus histogram-mode event bodies.

FULLY DECODED 2026-05-20.  Every field in every block, verified
byte-exact against BW's ASCII export across multiple histogram
fixtures.

The histogram-mode body is a stream of 32-byte fixed-length blocks,
one block per histogram interval.  Each block carries the per-interval
peak amplitude + zero-crossing frequency for all four channels (Tran,
Vert, Long, MicL).

────────────────────────────────────────────────────────────────────────────
Body layout (CONFIRMED 2026-05-20)
────────────────────────────────────────────────────────────────────────────

    [stream of 32-byte blocks]

Body length is approximately ``n_intervals * 32`` bytes plus a small
trailing remnant (1-9 bytes typically) at the very end.  Walker should
iterate 32-stride and stop before the tail.

────────────────────────────────────────────────────────────────────────────
32-byte block layout
────────────────────────────────────────────────────────────────────────────

    [0]    0x00                      always-zero tag
    [1]    segment_id  (uint8)       0x00..0x03 - 256 blocks per segment
    [2:4]  block_ctr  (uint16 LE)    resets each segment (0x0100, 0x0101, ...)
    [4]    0x0a (uint8)              constant marker (= 10)
    [5:7]  T_peak_count   uint16 BE  Tran peak (count x 0.005 -> in/s)
    [7:9]  T_halfperiod   uint16 BE  Tran half-period in samples (freq = 512 / halfp)
    [9:11] V_peak_count   uint16 BE
    [11:13] V_halfperiod  uint16 BE
    [13:15] L_peak_count  uint16 BE
    [15:17] L_halfperiod  uint16 BE
    [17:19] M_peak_count  uint16 BE  MicL peak (count -> dB via mic_count_to_db)
    [19:21] M_halfperiod  uint16 BE  MicL half-period in samples
    [21:23] 0x00 0x00                constant on standard blocks
    [24:28] 4-byte variable          purpose unknown (possibly CRC or timestamp delta)
    [28:32] block-end signature      see "Two block tails" below

**Every per-channel field is uint16 BIG-endian** (confirmed 2026-08-25).
Only ``block_ctr`` at [2:4] is little-endian.

HISTORY - two earlier readings of this block were wrong in ways that
cancelled out on quiet data:

  1. *peak as uint16 LE at [6:8]* - produced 268 in/s peaks on any
     interval whose next byte was non-zero.
  2. *peak as uint8 at [6] with an "annotation" byte at [7]* - correct
     for every peak below 256 counts (1.275 in/s), but it silently
     **clipped larger peaks**: the final interval of
     BE18193/T193LQ9K.OE0H reads 8.270 in/s in BW's export
     (1654 counts = 0x0676) and decoded as 0x76 = 118 = 0.590 in/s.
     The "annotation" byte was never an annotation - it is the high
     byte of the big-endian half-period, which is why it was non-zero
     exactly on the sub-Hz intervals BW renders as "<1.0".

Both readings also forced ``block[5] == 0`` via a bogus ``uint16 LE``
marker check at [4:6], which is what capped the peak at one byte.
The marker is ``block[4]`` alone.

Verified 2026-08-25 against 1211 production histograms paired with
their Blastware ASCII exports: **1211/1211 decode exactly** (interval
count plus every per-interval peak), and 842,442 per-interval
frequency comparisons match with **zero** mismatches.

Two block tails
---------------
Standard blocks end with ``1e 0a 00 00``.  The **final block of the
stream** ends with ``9c 06 00 42`` instead, and carries arbitrary bytes
at [21:23].  Rejecting it dropped the last interval of nearly every
histogram - and the last interval is frequently the one holding the
event peak, so the file's reported PPV came out low.  Observed in 1206
of 1211 production histograms, always positioned after every
standard-tail block.

Block-identification anchor: ``block[0] == 0x00`` AND
``block[4] == 0x0A`` AND the tail is one of the two signatures above;
standard-tail blocks additionally require ``block[22] == 0x00``.

────────────────────────────────────────────────────────────────────────────
Per-channel encoding
────────────────────────────────────────────────────────────────────────────

Geophone channels (Tran, Vert, Long):
  - peak_count × 0.005 = peak amplitude in in/s at Normal range
  - half-period in samples → freq_Hz = 512 / half-period

Microphone channel (MicL):
  - peak_count → dB via the same formula used by the waveform codec:
        dB = sign(c) × (81.94 + 20·log10(|c|))    for |c| ≥ 1
        dB = 0                                    for c == 0
  - half-period → freq_Hz = 512 / half-period (same as geo)

Frequency `>100 Hz` sentinel: the device emits half-period ≤ 5 when the
measured zero-crossing rate exceeds the geophone's measurement range
(since 512/5 = 102 Hz; the BW display rounds anything > 100 to ">100").

────────────────────────────────────────────────────────────────────────────
Output shape
────────────────────────────────────────────────────────────────────────────

``decode_histogram_body`` returns a per-channel dict matching the
waveform codec's shape so the rest of the pipeline (.h5 writer,
sidecar, viewer) consumes it without special-casing:

    {"Tran": [peak_count_i for each interval i],
     "Vert": [peak_count_i ...],
     "Long": [peak_count_i ...],
     "MicL": [peak_count_i ...]}

Values are in **16-count units for geo** (LSB = 0.005 in/s, matching
``decode_waveform_v2``) and **1-count units for mic** (matching the
waveform codec's mic convention).  Run through
``waveform_codec.decoded_to_adc_counts`` to scale geo to 1-count ADC.

Per-interval frequencies are NOT returned — they're auxiliary data,
not waveform samples.  Consumers needing frequencies can call
``decode_histogram_body_full()`` for the structured per-interval
record list.
"""

from __future__ import annotations

import struct
from typing import List, Optional, Tuple

# Block-end signature: constant `1e 0a 00 00` in bytes [28:32] of every
# real data block.  More distinctive than the byte-22 `00 00` (which
# matches many false positives), so we anchor on this.
_BLOCK_TAIL = b"\x1e\x0a\x00\x00"

# The final block of a histogram stream ends with this instead.  It is a
# real data block - same layout - and holds the last interval.  See the
# module docstring, "Two block tails".
_BLOCK_TAIL_TERMINAL = b"\x9c\x06\x00\x42"

_BLOCK_SIZE = 32

# Marker byte at block[4:6] of every histogram data block.  Used as
# additional validation that we're looking at a real block.
_BLOCK_MARKER = 10

# Geo peak scaling: stored as "count × 0.005 in/s" where 1 count = one
# 0.005 in/s display quantum.  Equivalent to the waveform codec's
# 16-count-unit output (1 unit = 0.005 in/s = 16 ADC counts).
_GEO_LSB_INS = 0.005

# Frequency formula: freq_Hz = _FREQ_NUMERATOR / half_period_samples.
# Empirically determined to be 512 (= sample_rate / 2, where sample rate
# is 1024 sps for the standard MiniMate Plus configuration).
_FREQ_NUMERATOR = 512


def _is_data_block(block: bytes) -> bool:
    """Tight identification of a histogram data block.

    Accepts both tail signatures.  ``block[4]`` alone is the marker -
    ``block[5]`` is the high byte of the Tran peak and is non-zero on any
    interval above 1.275 in/s, so it must not be part of the marker test.
    The ``block[22] == 0`` constraint is what keeps trailer content out,
    but it applies only to standard-tail blocks: terminal blocks carry
    arbitrary bytes there.
    """
    if len(block) < _BLOCK_SIZE:
        return False
    if block[0] != 0x00:
        return False
    if block[4] != _BLOCK_MARKER:
        return False
    tail = block[28:32]
    if tail == _BLOCK_TAIL:
        return block[22] == 0x00
    return tail == _BLOCK_TAIL_TERMINAL


def _decode_block(block: bytes) -> Optional[dict]:
    """Decode one 32-byte histogram block.  Caller must have validated
    with ``_is_data_block`` first.

    Returns a record with per-channel peak counts (uint8) and
    half-periods (uint16 LE).
    """
    # Every per-channel field is uint16 BIG-endian; only block_ctr is LE.
    # See the module docstring for the two superseded readings and why
    # each looked correct on quiet data.
    def _be16(i: int) -> int:
        return (block[i] << 8) | block[i + 1]

    t_peak  = _be16(5)
    t_halfp = _be16(7)
    v_peak  = _be16(9)
    v_halfp = _be16(11)
    l_peak  = _be16(13)
    l_halfp = _be16(15)
    m_peak  = _be16(17)
    m_halfp = _be16(19)
    segment_id = block[1]
    block_ctr  = block[2] | (block[3] << 8)
    var_meta   = bytes(block[24:28])
    return {
        "segment_id":  segment_id,
        "block_ctr":   block_ctr,
        "t_peak":      t_peak,
        "t_halfp":     t_halfp,
        "v_peak":      v_peak,
        "v_halfp":     v_halfp,
        "l_peak":      l_peak,
        "l_halfp":     l_halfp,
        "m_peak":      m_peak,
        "m_halfp":     m_halfp,
        "meta_var":    var_meta,
        "is_terminal": block[28:32] == _BLOCK_TAIL_TERMINAL,
    }


def walk_body(body: bytes) -> List[dict]:
    """Walk the body and return one dict per histogram interval.

    Iterates 32-byte strides from offset 0.  Yields a decoded record
    for every block that passes ``_is_data_block`` validation.  Stops
    when the remaining bytes are too short to form a complete block.

    In Histogram+Continuous mode the body interleaves data blocks with
    other 32-byte content (likely continuous-mode waveform blocks) that
    fail the data-block validation; the walker naturally skips them
    without losing 32-byte alignment.  Use ``block_ctr`` from each
    returned record to map back to the original interval index — the
    record list is sparse when other block types are interleaved.
    """
    records: List[dict] = []
    for off in range(0, len(body) - _BLOCK_SIZE + 1, _BLOCK_SIZE):
        blk = body[off:off + _BLOCK_SIZE]
        if not _is_data_block(blk):
            # Hit non-block content (likely a sync or stream marker).
            # Continue walking — block alignment is fixed at 32-stride
            # from offset 0, so we don't lose alignment by skipping.
            continue
        decoded = _decode_block(blk)
        if decoded is None:
            # Block validated as a histogram block but had peak fields
            # outside the plausible range — undocumented extension.
            # Skip rather than propagating bogus PVS contributions.
            continue
        records.append(decoded)
    return records


def decode_histogram_body(body: bytes) -> Optional[dict]:
    """Decode a histogram-mode body into per-channel peak-sample arrays.

    Returns ``{"Tran": [...], "Vert": [...], "Long": [...], "MicL": [...]}``
    where each channel's list contains one peak value per histogram
    interval (in the same units the waveform codec uses: 16-count units
    for geo, 1-count ADC units for mic).  Returns ``None`` if the body
    doesn't contain any valid histogram blocks.

    To convert to physical units:
      - Geo channels: ``count * 0.005`` = peak in in/s at Normal range
        (or run through ``waveform_codec.decoded_to_adc_counts`` first
         to get 1-count ADC values, then ``count / 32767 * 10.0`` for in/s)
      - Mic channel:  use ``waveform_codec.mic_count_to_db(count)``
    """
    records = walk_body(body)
    if not records:
        return None
    return {
        "Tran": [r["t_peak"] for r in records],
        "Vert": [r["v_peak"] for r in records],
        "Long": [r["l_peak"] for r in records],
        "MicL": [r["m_peak"] for r in records],
    }


def decode_histogram_body_full(body: bytes) -> Optional[List[dict]]:
    """Decode a histogram-mode body into the full per-interval record list.

    Same data as ``decode_histogram_body`` but in a structured form that
    preserves the half-period (frequency) data for each channel + the
    per-block segment_id, block_ctr, and 4-byte variable metadata.
    Useful for diagnostic tools, sidecar enrichment, and future-codec
    work.

    Returns ``None`` if the body has no valid blocks.
    """
    records = walk_body(body)
    return records if records else None


def half_period_to_hz(halfp: int) -> Optional[float]:
    """Convert a half-period in samples to frequency in Hz.

    Returns ``None`` for half-period ≤ 5 — the device emits values in
    that range when the measured zero-crossing rate exceeds 100 Hz
    (the BW display reports `>100 Hz` for such cases).  Callers can
    treat ``None`` as the `>100 Hz` sentinel.
    """
    if halfp <= 5:
        return None
    return _FREQ_NUMERATOR / halfp


def geo_count_to_ins(count: int) -> float:
    """Convert a histogram geo peak count to in/s at Normal range."""
    return count * _GEO_LSB_INS