udpscope/src/measure.hpp

// SPDX-License-Identifier: MIT
// SPDX-FileCopyrightText: 2026 Martino Ferrari
#pragma once

#include <cmath>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <string>
#include <vector>

namespace udpscope {

/// What a measurement computes over the visible window of one signal.
enum class MeasureType {
  Average,        ///< the mean of every visible sample
  Min,            ///< the smallest visible value
  Max,            ///< the largest visible value
  /// The largest swing in the window: `max - min`, the two extremes the eye would measure
  /// with a ruler across the waveform -- ripple and overshoot included, which is what makes it
  /// the *extremes* and not a level.
  PeakToPeak,
  /// The delta between the two **flat** levels, `top - bottom` (a bench scope's `Vamplitude`,
  /// which is `Vtop - Vbase`). Deliberately *not* the peak-to-peak swing: a band-limited edge
  /// overshoots and rings, so the extremes of a nominal 2 V square read about 2.36 V, while the
  /// levels the signal actually dwells at are ±1 and their delta exactly 2. Where the window
  /// has no density estimate to separate a plateau from its ripple (shorter than one bucket, or
  /// an HF-rejected walk) it degenerates to the window's own peak to peak, which is the same
  /// number on a waveform with nothing to exclude. Never `max(|value|)`: that is a *distance
  /// from zero*, which for a signal riding on an offset is a different and less useful number.
  Amplitude,
  /// The signal's high level: the densest band of the window's upper side, whose samples' exact
  /// mean is the value -- the *flat top*, so a plateau's ripple and an edge's overshoot neither
  /// define it nor enter it.
  Top,
  /// The signal's low level, the same way from the window's lower side.
  Bottom,
  PulseDuration,  ///< how long the signal spends on one side of the threshold
  Period,         ///< the time between two successive crossings of the threshold
  Frequency,      ///< 1 / period
};

/// Type-specific behaviour. Every field is meaningful only for the types named below; the
/// others ignore it (an option that does nothing for a type is left out of the dialog
/// rather than stored and ignored).
struct MeasureOptions {
  /// `Top`/`Bottom`/`PulseDuration`/`Period`/`Frequency`: where the split is, as a share of
  /// the window's own min..max. The middle by default, which is what a symmetric waveform
  /// wants; a 3.3 V logic signal with ringing wants it nearer the top.
  double threshold_pct = 50.0;
  /// `PulseDuration`: which side of the threshold is "the pulse". `Period`/`Frequency`
  /// use it as the edge direction, so one cycle is counted once.
  bool above = true;
  /// Draw the value on the plot: a ruler at a level (top/bottom), both rulers at the
  /// extremes (min/max), the pulse's two edges as vertical markers (pulse duration). The
  /// markers are readouts -- they follow the window, and the user cannot drag them.
  bool show_on_plot = false;
  /// `PulseDuration`/`Period`/`Frequency`: the half-width of the hysteresis band about the threshold, as
  /// a percentage of the window's own peak-to-peak. A crossing counts only once the signal has passed the
  /// band, so an edge count is a count of the *waveform's* cycles rather than of its noise.
  ///
  /// The default is measured rather than chosen: 5 % reads the mock's 0.1 interferer correctly and fails
  /// at 0.2 (`period=0.075` for a true 0.2), while 10 % holds to an interferer amplitude of 0.3 -- a third
  /// of the waveform -- and to white noise at sigma=0.05. Above that the band starts eating a real
  /// waveform's own extremes on anything less square than a sine. 0 reduces the walk exactly to a single
  /// level.
  double noise_pct = 10.0;
  /// Every type: the corner of a one-pole low pass applied to each sample before *either* walk reads it,
  /// so the extremes, the mean and the crossings all describe the filtered signal. 0 is off. In Hz.
  ///
  /// It rejects what a band cannot: an impulse large enough to cross the band on its own. It is not free
  /// -- a filtered window cannot use the LOD pyramid, because a bucket aggregates unfiltered samples --
  /// so a non-zero corner costs a full raw walk.
  double hf_fc = 0.0;
};

/// One entry in the measurements pane, and everything needed to recompute it: the signal by
/// *name* (a CONFIG renumbers channels, so an index would start measuring another signal)
/// and the type with its options.
struct Measure {
  std::string signal;
  MeasureType type = MeasureType::Average;
  MeasureOptions options;
};

/// One thing to draw on the plot for a measurement: a *level* (a horizontal ruler) or a
/// *time* (a vertical marker). A result carries up to two, and they need not be the same kind
/// -- `min` is a ruler at the value and a marker where the value occurred.
struct MeasureMarkSpec {
  bool level = false;
  double value = 0.0;
};

/// What a measurement computed. `valid` is what the pane shows as `--` when a type has
/// nothing to measure in the window (a period on a flat signal, a pulse that never crosses).
struct MeasureResult {
  double primary = 0.0;    ///< the value the pane prints
  double secondary = 0.0;  ///< pulse duration also reports the duty; unused otherwise
  bool has_secondary = false;
  std::vector<MeasureMarkSpec> marks;  ///< what to draw, empty when the type draws nothing
  bool valid = false;
};

/// The density estimate behind `Top`/`Bottom`/`Amplitude`: one bin per band of the window's
/// value range, over the whole LOD buckets pass one walked.
///
/// A flat level is a *density* question, not an extreme one: on a band-limited square the
/// plateau holds most of the samples while the edge's ringing holds a few percent, so the
/// densest band is the level the signal dwells at and the ringing is what a mean over a side
/// reports instead. One bucket is one observation at its midpoint -- a bucket is 64 samples
/// held as one aggregate, so counting its samples would weight a flat region by nothing at all
/// and a ringing one by 64 -- and each bin also carries the exact sum and count of the samples
/// its buckets hold, which is what makes the reported level an exact mean rather than a bin
/// centre. The bucket list is the one pass one already fetches; no pass reads raw samples for
/// this, which is the whole reason it is built from buckets and not from a sample histogram
/// (that would be a second O(window) walk).
struct MeasureFlatHist {
  /// Bins across the buckets' own range. The band the level is read from is picked at this
  /// resolution and is a share of the *window's* peak-to-peak, so the bin count only sets how
  /// finely the densest band can be located.
  static constexpr std::size_t kBins = 64;
  double lo = 0.0;  ///< the bin range: the whole buckets' own smallest `lo`..
  double hi = 0.0;  ///< ..and largest `hi`. Empty when nothing was folded.
  double sum[kBins] = {};
  std::size_t count[kBins] = {};
  std::size_t buckets[kBins] = {};
  /// Whole buckets folded in: 0 means this window has no histogram at all (a window shorter
  /// than a bucket, or an HF-rejected walk, which cannot use the pyramid), and the flat-level
  /// types then fall back to their split-relative readings.
  std::size_t folded = 0;
};

/// The window's own extremes, which every threshold-relative type scales from and which
/// several types report directly. Public so a caller holding *several* measurements on one
/// signal can pay for the scan once: it is O(visible samples) and a live view at 1 MS/s makes
/// that the difference between a frame budget and a stall.
struct MeasureStats {
  double lo = 0.0;
  double hi = 0.0;
  /// The sum and count of the finite samples: the mean, and the mean over a level, are both
  /// derived from them, so the window is read once for every type that is not an edge count.
  double sum = 0.0;
  std::size_t count = 0;
  /// The *window index* of the first sample to reach each extreme -- "min" is a level and a
  /// place, and finding the place is why this is not just the two numbers. An index rather than
  /// a stamp because a stamp costs a ring read and the walk moves an extreme many times, while
  /// the caller needs the time exactly twice.
  std::size_t lo_i = 0;
  std::size_t hi_i = 0;
  /// How many samples from `lo_i`/`hi_i` still have to be searched to turn it into a *sample*,
  /// or 0 when it already is one. A bucket pass knows only which bucket holds the extreme, so
  /// the caller finishes with one raw `measure_scan()` over that bucket alone -- two buckets at
  /// most, whatever the window's size.
  std::size_t lo_search = 0;
  std::size_t hi_search = 0;
  /// The time of each extreme. Filled by the caller from `lo_i`/`hi_i` once the walk is done;
  /// `compute_measure()` reads these, so a caller that skips the resolve reports t = 0.
  double at_lo = 0.0;
  double at_hi = 0.0;
  /// The flat levels' density estimate, folded from the same whole buckets pass one walked.
  /// The *values* are resolved by `measure_flat_levels()`, which needs the split and the
  /// window's span and so cannot be done until the walk is over.
  MeasureFlatHist flat;
  bool valid = false;
};

/// The HF reject, as a walk carries it: one pole, one state word.
///
/// `a == 1` is a pass-through and is the default, so an unfiltered walk pays one multiply and no branch,
/// and every caller can step it unconditionally. A non-finite sample passes through unchanged and leaves
/// the state alone -- a dropout is a real reading, and letting it poison `y` would make every sample
/// after it non-finite too.
struct MeasureFilter {
  double a = 1.0;
  double y = 0.0;
  bool have = false;

  double step(double v) noexcept {
    if (!std::isfinite(v)) return v;
    if (!have) {
      y = v;
      have = true;
      return y;
    }
    y += a * (v - y);
    return y;
  }
  void reset() noexcept {
    y = 0.0;
    have = false;
  }
};

/// `a = 1 - exp(-2*pi*fc*dt)`, clamped just under Nyquist so `a` cannot exceed 1. `fc == 0` or `dt == 0`
/// leaves it at 1, which is the pass-through. Resets the state.
void measure_filter_set(MeasureFilter& f, double fc, double dt);

/// One chunk of a window's values, in time order, exactly as the rings hold them.
///
/// The engine is *fed* the window this way rather than handed a materialized copy. A visible
/// window is millions of samples and the rings hold `f32`, so widening the whole thing into two
/// `std::vector<double>` cost a 16-byte-per-sample write and then a read of it for every pass --
/// 60 % of the pane's time before a single statistic had been computed, and 12 ns a sample over a
/// 1.7 M-sample window. Only the *values* are streamed: a stamp is needed at the extremes and at
/// the level crossings, which are a handful in a window of millions, so it is fetched by index.
struct MeasureChunk {
  const float* value = nullptr;
  std::size_t count = 0;
  /// The window offset of `value[0]`, so `time_of()` below can be asked for one absolute sample.
  std::size_t offset = 0;
};

/// One whole LOD bucket of a window, exactly as `plot::data::MinMaxPyramid` already holds it.
///
/// The pyramid *is* "evaluate only on the new samples": the producer folds every sample into it
/// as it arrives, at amortised O(1), so a bucket is 64 samples the measurement walks never have
/// to read. A 1.7 M-sample window is 26 k buckets, and only the handful that *straddle* a level
/// are opened.
///
/// A bucket aggregates **unfiltered** samples, so a walk with an HF reject must not use these at all --
/// that is what makes a non-zero `hf_fc` cost a full raw walk.
struct MeasureBucket {
  double lo = 0.0;
  double hi = 0.0;
  double sum = 0.0;
  std::size_t count = 0;   ///< finite samples folded in; 0 = a bucket with nothing usable
  std::size_t offset = 0;  ///< window offset of the bucket's first sample
  std::size_t span = 0;    ///< samples the bucket covers, finite or not
};

/// Pass one, chunk by chunk: both extremes with the *index* of the first sample to reach each,
/// and the sum and count of the finite samples. Feed the chunks in window order into a zeroed
/// `stats` -- the first sample to reach an extreme is the one the mark names, so a later chunk
/// that merely ties it must not move it.
void measure_scan(const MeasureChunk& chunk, MeasureStats& stats, MeasureFilter& filter);

/// Pass one from whole buckets: the same figures, exactly, for 1/64 of the reads. The extreme's
/// position comes out as a *bucket* (`lo_search`/`hi_search` non-zero), which one raw
/// `measure_scan()` of that bucket then pins. Interleave with `measure_scan()` in window order.
///
/// `flat_hist` folds the same list into `stats.flat` for the flat-level types, in a second tight
/// loop over the array (the bin range is the buckets' own, which the first loop is what
/// discovers). That loop is O(buckets) and reads no sample, so the `measure cost:` line's `raw`
/// is unchanged by it -- the cost is buckets, which is the budget the pane is on. It is the
/// caller's because only the caller knows whether any measurement reads a flat level.
void measure_scan_buckets(const MeasureBucket* buckets, std::size_t count, MeasureStats& stats,
                          bool flat_hist);

/// Everything the threshold-relative types need from a second walk, at one level: the sum and
/// count of the samples on either side of it, and its crossings reduced to running figures -- the
/// first and last *enter* edge of each side, how many there were, and the summed width of the
/// intervals that closed. All of it is O(1) state, so `Top`, `Bottom`, `PulseDuration`, `Period`
/// and `Frequency` are all answered from one pass: adding a *type* must not add a pass, and nor
/// must adding a second measurement of the same type at the same level.
struct MeasureLevel {
  double level = 0.0;
  /// Half-width of the hysteresis band about `level`: a crossing is counted only once the signal
  /// has passed `level + hysteresis` going up or `level - hysteresis` going down, and a sample
  /// inside the band decides nothing. It is a Schmitt trigger, and it is what makes an edge count
  /// a count of the *waveform's* cycles rather than of its noise: a single level is recrossed
  /// several times per true crossing by any interferer faster than the fundamental, because near
  /// the mid level the fundamental's own slope is at its smallest. Measured on a 5 Hz sine with a
  /// 500 Hz interferer at a tenth of the amplitude: 9 true crossings counted as 79, reading
  /// `period=0.025622 freq=39.029` for a true 0.2 s / 5 Hz. With the band, `period=0.200188`.
  ///
  /// It moves the *edge decision* only. The `sum`/`count` split below stays at the bare `level`,
  /// so a band can never change it -- the split is the split, and a sample's side of it is a
  /// fact rather than a threshold decision. `measure_set_level()` is what sizes it; zero reduces
  /// the walk exactly to a single level.
  double hysteresis = 0.0;
  /// The samples either side of `level`, exactly: `Top`/`Bottom`'s reading on a window with no
  /// bucket histogram to read a flat level from (shorter than one bucket, or an HF-rejected
  /// walk). They are folded on every window -- pass two walks it for the crossing types anyway
  /// -- but a window that has a histogram is read from that instead, because a mean over a side
  /// reports the ripple and the overshoot rather than the level the signal dwells at.
  /// Index 0 is the side below the level (a sample exactly *at* it counts here), index 1 the
  /// side above it.
  double sum[2] = {0.0, 0.0};
  std::size_t count[2] = {0, 0};
  double first_enter[2] = {0.0, 0.0};
  double last_enter[2] = {0.0, 0.0};
  std::size_t enters[2] = {0, 0};
  double width_sum[2] = {0.0, 0.0};
  std::size_t widths[2] = {0, 0};
  /// Walk state carried between chunks: the side the signal was last *decided* to be on (a sample
  /// equal to the level decides nothing, so the walk holds the previous side), and the interval
  /// per side that has been entered but not yet left.
  bool have_side = false;
  bool side = false;
  double open_at[2] = {0.0, 0.0};
  bool open[2] = {false, false};
};

/// Fill a level accumulator's `level` and `hysteresis` from the window's extremes and the share
/// `measure_level_pct()` gave. The one place the band is sized, so every caller gets the same rule:
/// a tenth of the window's peak-to-peak, *clamped* to half the room between the level and the
/// nearer extreme. The clamp is what keeps a threshold near an extreme usable -- `PulseDuration`
/// lets the user put one at 90 %, and a fixed band there puts the upper threshold above every
/// sample in the window, which reads as nothing measured at all rather than as a noisy reading.
/// `noise_frac` is `MeasureOptions::noise_pct / 100`. It is the caller's because the band is a
/// per-measurement setting and this is called once per *level*, which several measurements may share; the
/// clamp to half the room to the nearer extreme stays here, so every caller gets it.
void measure_set_level(MeasureLevel& level, const MeasureStats& stats, double pct, double noise_frac);

/// Pass two: one chunk into every level at once. The usual case is a single level -- every
/// crossing type reads the mid level and `Top`/`Bottom` default to it too -- so a pane full of
/// measurements on one signal costs two walks, not two per measurement.
void measure_accumulate(const MeasureChunk& chunk,
                        const std::function<double(std::size_t)>& time_of, MeasureLevel* levels,
                        std::size_t level_count, MeasureFilter& filter);

/// Pass two from one whole bucket, into one level. True when the bucket answered it: a bucket
/// every one of whose samples is on the same side of the level holds no crossing, so its sums
/// fold in whole and the walk's side is decided by the bucket rather than by 64 samples. An edge
/// *into* such a bucket lands on its first sample, exactly, which is why that is the stamp taken.
///
/// False means the caller must feed the bucket's raw samples to `measure_accumulate()` for this
/// level instead. That is the only place a measurement still reads raw samples in bulk, and on a
/// clean waveform it is a bucket or two per crossing. Three things make it happen: the bucket
/// straddles the level (so the sums cannot be binned); it holds a non-finite sample, which would
/// make the edge's position a guess; or it is where the signal passes `level +/- hysteresis` and the
/// *position* of that crossing is inside the bucket. The band is deliberately not folded into the
/// first test -- see the comment in the implementation for why refusing every bucket inside the band
/// costs a share of the window rather than a bucket.
bool measure_accumulate_bucket(const MeasureBucket& bucket,
                               const std::function<double(std::size_t)>& time_of,
                               MeasureLevel& level);

/// The level a type reads, as a share of the window's own min..max, or false when it reads none
/// (the pure statistics, answered from `MeasureStats` alone). `Top`/`Bottom` read the
/// measurement's own threshold; the crossing types read the mid level, deliberately -- see
/// `measure_uses_threshold()`. `Amplitude` reads no level either: its value is the delta between
/// the two flat levels, which the histogram answers, and the split's two means are only a
/// fallback it takes from the window's own extremes -- see `compute_measure()`.
bool measure_level_pct(MeasureType type, const MeasureOptions& options, double& pct);

/// The window's two flat levels, from the bucket histogram in `stats.flat`. `top` is the densest
/// band of the bins *above* `lo + pct * (hi - lo)`, `bottom` the same below it -- for a
/// crossing type's mid level that is the two plateaus of a pulse train, and for `Top`/`Bottom`'s
/// own threshold it is "which side is the level I want". Each side is false when it holds no
/// bucket, and both are false when the window has no histogram at all.
///
/// Two things this does not do. It does not read a sample: every figure comes from buckets pass
/// one already fetched. And it does not read the window raw, which is why the answer is
/// independent of how the window happens to be aligned -- the band is a *value* band, so the
/// level it reports does not move when the window does.
struct MeasureFlatLevels {
  double top = 0.0;
  bool top_ok = false;
  double bottom = 0.0;
  bool bottom_ok = false;
};
MeasureFlatLevels measure_flat_levels(const MeasureStats& stats, double pct);

/// The engine: one measurement's answer, O(1), from the two passes' results. `level` is the
/// accumulator for `measure_level_pct(type, options)` and may be null for the types that read
/// none. The function reads nothing else, which is what keeps the answer independent of the plot.
MeasureResult compute_measure(const MeasureStats& stats, const MeasureLevel* level,
                              MeasureType type, const MeasureOptions& options);

/// A measurement type's name, as the pane's row and the dialog's dropdown spell it.
const char* measure_type_name(MeasureType type);
/// The number of entries the dialog's dropdown holds; the index is `MeasureType`.
int measure_type_count();
MeasureType measure_type_at(int index);
/// Whether a type uses the threshold (see `MeasureOptions`), so the dialog can show the
/// field only then. Deliberately not the case for `Period`/`Frequency`: a period is the
/// spacing of the crossings of the *mid* level, one rising edge per cycle, and offering a
/// threshold there only invites a setting that makes the count wrong.
bool measure_uses_threshold(MeasureType type);
/// Whether a type uses the `above` side (only `PulseDuration`, where it says which half of
/// the waveform is the pulse).
bool measure_uses_side(MeasureType type);
/// Whether a type reads the hysteresis band, so the dialog can show the field only then. The band moves
/// the *edge* decision and nothing else, so it is exactly the crossing types: `Top`/`Bottom` read a flat
/// level out of the bucket histogram, which no band touches.
bool measure_uses_noise(MeasureType type);
/// Whether a type draws its result on the plot when asked.
bool measure_can_show_on_plot(MeasureType type);

/// Running statistics over a measurement's successive results.
///
/// One observation per *recompute* of the measurements pane, never one per sample: a
/// measurement is a figure derived from a whole window, and counting its samples would
/// describe the window rather than the measurement's stability. `n` is therefore a count of
/// passes, which is what makes the readout comparable with a bench scope's.
///
/// The mean and spread are Welford's recurrence rather than running sums of `v` and `v*v`:
/// the results are absolute signal values, so a channel sitting at a large offset with a
/// small ripple -- the case the readout is *for* -- cancels to nothing in a naive
/// sum-of-squares long before the ripple is visible.
struct MeasureStatsAccum {
  std::uint64_t n = 0;
  double mean = 0.0;  ///< Welford's running mean.
  double m2 = 0.0;    ///< Welford's sum of squared deviations from the running mean.
  double lo = 0.0;
  double hi = 0.0;

  void add(double v) noexcept {
    if (n == 0) {
      lo = v;
      hi = v;
    } else {
      if (v < lo) lo = v;
      if (v > hi) hi = v;
    }
    ++n;
    const double delta = v - mean;
    mean += delta / static_cast<double>(n);
    m2 += delta * (v - mean);
  }

  void reset() noexcept { *this = MeasureStatsAccum{}; }

  /// The sample standard deviation (n-1), and 0 rather than NaN for fewer than two
  /// observations: the readout is a figure in a label, and a single pass has no spread.
  double sd() const noexcept {
    return n < 2 ? 0.0 : std::sqrt(m2 / static_cast<double>(n - 1));
  }
};

}  // namespace udpscope