ssvep.classification — the decode engine

Part of the API reference, generated from the docstrings in src/ssvep.

Decoding

ssvep.classification.decoding

SSVEP frequency-recognition decoders (domain-generic, array-based).

All functions take plain NumPy arrays and a sampling rate, so the same code runs offline (batch) and online (per-window in an LSL loop). No MNE / BIDS / study dependencies.

• CCA, FBCCA: calibration-free (sinusoidal reference templates; no training).
• TRCA: template-trained; gate it with trca_feasible() — it needs several phase-consistent trials per class and degrades to chance otherwise.

Trial array convention: (n_channels, n_times); batches: (n_trials, n_channels, n_times). Reference frequencies are in Hz.

UndecodableWindow

Bases: ValueError

No candidate can be scored: every one sits on a mains line at this window's resolution.

Raised rather than returning an argmax over all-missing scores, which would silently name the first candidate (review of #302).

Source code in src/ssvep/classification/decoding.py

class UndecodableWindow(ValueError):
    """No candidate can be scored: every one sits on a mains line at this window's resolution.

    Raised rather than returning an argmax over all-missing scores, which would silently name the
    first candidate (review of #302).
    """

TRCA

Ensemble Task-Related Component Analysis (Nakanishi et al. 2018).

Calibration-based: one spatial filter + template per class. Requires several phase-consistent trials per class — check trca_feasible(y) first.

Epoch-length rules (every epoch must be time-locked to stimulus onset): - Within a class, training epochs must be equal length — the cross-trial covariance is sample-aligned, so ragged within-class input raises. - Across classes, templates may differ in length: each class's filter and template are learned independently, with no cross-frequency coupling. So a whole-number-of-cycles (per-frequency) window is fine, and there is no need for a common length across classes. - predict compares a test window to each class over their shared onset-locked interval (both cropped to the shorter), so the test window need not match the training length either. - Caveat: a correlation over more samples trends higher, so scoring classes at different lengths biases the argmax toward the longer templates. Prefer a common decision-window length unless you specifically want per-frequency lengths.

Source code in src/ssvep/classification/decoding.py

class TRCA:
    """Ensemble Task-Related Component Analysis (Nakanishi et al. 2018).

    Calibration-based: one spatial filter + template per class. Requires several
    phase-consistent trials per class — check `trca_feasible(y)` first.

    Epoch-length rules (every epoch must be time-locked to stimulus onset):
    - Within a class, training epochs must be equal length — the cross-trial
      covariance is sample-aligned, so ragged within-class input raises.
    - Across classes, templates may differ in length: each class's filter and
      template are learned independently, with no cross-frequency coupling. So a
      whole-number-of-cycles (per-frequency) window is fine, and there is no need
      for a common length across classes.
    - `predict` compares a test window to each class over their shared
      onset-locked interval (both cropped to the shorter), so the test window need
      not match the training length either.
    - Caveat: a correlation over more samples trends higher, so scoring classes at
      different lengths biases the argmax toward the longer templates. Prefer a
      common decision-window length unless you specifically want per-frequency
      lengths.
    """

    def __init__(self, n_harmonics: int = 3):
        self.n_harmonics = n_harmonics
        self.classes_ = self.templates_ = self.W_ = None

    @staticmethod
    def _trca_filter(trials):
        """Fit one class's TRCA spatial filter, restricted to the rank-revealing subspace of the
        (per-trial, time-centred) channel data (:func:`_rank_revealing_basis`, :data:`RANK_RTOL` —
        the same rank fix #367 gives CCA).

        After a common-average reference over every decoded channel, the channels sum to zero at
        every timepoint by construction, so ``Q`` is exactly singular along the all-ones direction
        — but ``np.linalg.solve`` only raises on an *exact* zero pivot, and floating-point rounding
        never delivers one (measured ``sigma_min/sigma_max(Q) ~ 1e-17``), so the ``pinv`` fallback
        this used to have was dead code: ``solve`` always "succeeded", returning a filter with a
        large, rounding-determined component along the null direction — the fitted filter was
        mostly the CAR null direction itself (``|w . 1| / |w|`` measured 0.77-1.0) (#368).

        Fix: project each trial onto the retained basis *before* forming ``S``/``Q`` (as the issue
        suggests), so the null direction never enters the eigenproblem, then map the generalised
        eigenvector found there back to channel space with the same basis. Unlike the old code,
        this does not depend on a solver noticing a near-singular matrix: the direction is removed
        outright.

        ``S`` is symmetric (every ordered pair ``i != j`` is summed alongside its transpose-partner
        ``j, i``) and ``Q`` is, by construction, symmetric positive definite on the retained
        subspace (see the comment at its use below) — so the generalised eigenproblem is solved
        with :func:`scipy.linalg.eigh`, which is exact for a symmetric-definite pencil: real
        eigenvalues in ascending order, no ``.real`` truncation of an assumed-real result from a
        general (non-symmetric) solver.

        Normalisation (found in review of #368): a bare generalised eigenvector is unique only up
        to an arbitrary, solver-chosen scale *and* sign — harmless for a single class scored on its
        own (:meth:`predict`'s per-class correlation is scale/sign-free), but ``predict`` actually
        flattens *every* class's filtered output into one vector before correlating (the ensemble
        score Nakanishi et al. describe), and a flattened multi-row vector's correlation is not
        invariant to an independent sign/scale change in only some of its rows. Two representations
        of the identical subspace (e.g. this basis vs. literally dropping a redundant channel) can
        then disagree on scale (each solver normalises its eigenvector's own norm, which depends on
        the space it was solved in) and sign (arbitrary per solve) even though the decision either
        model reaches alone would be the same. Both are pinned here so that two equivalent
        parametrisations produce byte-identical ensemble filters, not just equivalent ones:

        - **Scale**: rescale ``w`` to unit standard deviation of its own filtered output over this
          class's training trials — every class's filter then contributes equal weight to the
          flattened ensemble vector, regardless of which subspace representation produced it.
        - **Sign**: flip ``w`` so the filtered template's largest-magnitude sample is non-negative
          (not, as an earlier draft used, the *mean* over time of the filtered template): on
          per-trial time-centred input — ordinary preprocessing, and not distinguishable from the
          uncentred case by looking at this function alone, since it time-centres internally either
          way — the mean of an already near-zero-mean template is itself near zero, so a mean-based
          rule compares noise against noise and flips under a perturbation many orders of magnitude
          smaller than any real acquisition's noise floor (measured: a coin toss under a
          1e-15-relative perturbation). The largest-magnitude sample carries the residual DC the
          mean rule was chasing (band-pass crop leakage, or the trial-averaged noise mean on
          uncentred data) at the same absolute scale, so it is stable under the same perturbation on
          every input shape tested — uncentred, per-trial time-centred, band-passed, each with and
          without CAR (see ``test_filter_is_invariant_to_a_tiny_perturbation`` in
          ``tests/test_trca_rank.py``, including its time-centred case, which fails under the old
          mean-based rule). The sign is not claimed to be "physically meaningful" beyond that
          stability; it is whichever polarity the real, if small, cross-trial signal/noise structure
          at the data's own scale happens to produce, just measured at the sample that is least
          sensitive to rounding.
        """
        n_tr, n_ch, _ = trials.shape
        centred = [t - t.mean(1, keepdims=True) for t in trials]
        cat = np.concatenate(centred, axis=1)
        basis = _rank_revealing_basis(cat)          # (n_ch, r): channel directions the data spans
        if basis.shape[1] == 0:
            raise ValueError("TRCA: trial data has no rank (all channels constant or zero) — "
                             "cannot fit a spatial filter")
        proj = [basis.T @ x for x in centred]        # (r, n_times) per trial
        r = basis.shape[1]
        S = np.zeros((r, r))
        for i in range(n_tr):
            for j in range(n_tr):
                if i != j:
                    S += proj[i] @ proj[j].T
        catp = np.concatenate(proj, axis=1)
        Q = catp @ catp.T
        # Q is, by construction, exactly `basis`'s own retained singular values squared (on the
        # diagonal, up to the rounding a matrix product introduces) — every one of them is above
        # RANK_RTOL of the largest by the basis's own selection, so Q is provably symmetric positive
        # definite on this subspace. No pinv fallback is needed or reachable here; a LinAlgError
        # would mean the trial data itself is pathological in a way this rank fix does not address.
        vals, vecs = eigh(S, Q)                      # ascending: largest eigenvalue is last
        w = basis @ vecs[:, -1]

        out = np.concatenate([w @ x for x in centred])
        w = w / out.std()                            # unit output variance (see docstring: scale)
        y = w @ trials.mean(0)                        # this class's filtered, trial-averaged template
        if y[np.argmax(np.abs(y))] < 0:
            w = -w                                    # largest-magnitude sample made non-negative
        return w

    @staticmethod
    def _class_trials(X, y, c):
        """Trials for class ``c`` as ``(n_tr, n_ch, n_times)``.

        Accepts an ``(n_trials, n_ch, n_times)`` array or a sequence of
        ``(n_ch, n_times)`` trials (lengths may vary across classes). Unequal
        lengths *within* the class are an error, not a silent truncation.
        """
        trials = [np.asarray(X[i]) for i in np.where(y == c)[0]]
        try:
            return np.stack(trials)
        except ValueError as e:
            raise ValueError(
                f"TRCA: class {c!r} has unequal-length epochs; within a class, "
                "training epochs must be equal length and onset-aligned"
            ) from e

    def fit(self, X, y):
        y = np.asarray(y)
        self.classes_ = np.unique(y)
        filters, templates = [], {}
        for c in self.classes_:
            trials = self._class_trials(X, y, c)
            filters.append(self._trca_filter(trials))
            templates[c] = trials.mean(0)
        self.W_ = np.column_stack(filters)
        self.templates_ = templates
        return self

    def predict(self, X):
        preds = []
        for trial in X:
            trial = np.asarray(trial)
            scores = []
            for c in self.classes_:
                tmpl = self.templates_[c]
                n = min(trial.shape[-1], tmpl.shape[-1])   # shared onset-locked span
                r = np.corrcoef((self.W_.T @ trial[:, :n]).ravel(),
                                (self.W_.T @ tmpl[:, :n]).ravel())[0, 1]
                scores.append(0.0 if np.isnan(r) else r)
            preds.append(self.classes_[int(np.argmax(scores))])
        return np.array(preds)

Look

Bases: NamedTuple

One evaluated look of an :class:EarlyStopRule: the window (seconds, and the sample count it was scored on), the candidate scores, :func:decision_margin on them, and confident, which is :func:is_confident and "this is not the maximum window" (the last look is the decision whether or not it is confident, so it is never an early stop).

Source code in src/ssvep/classification/decoding.py

class Look(NamedTuple):
    """One evaluated look of an :class:`EarlyStopRule`: the window (seconds, and the sample count it
    was scored on), the candidate ``scores``, :func:`decision_margin` on them, and ``confident``,
    which is :func:`is_confident` **and** "this is not the maximum window" (the last look is the
    decision whether or not it is confident, so it is never an *early* stop)."""
    window_s: float
    n_samples: int
    scores: np.ndarray
    margin: float
    confident: bool

EarlyStopRule

The growing-window early-stopping rule (#23), defined once for the offline replay (:func:early_stop_decode) and the live decoder (:class:ssvep.runtime.online.LiveDecoder).

The configuration: candidate freqs (the frequencies the references are built at), the scoring (method, n_harmonics, band, line_hz, filter-bank keywords, as :func:classify takes them), and the schedule (min_window_s, max_window_s, step_s, confidence_margin; a None margin is the method's calibrated default). max_window_s is required here: "the whole trial" is something only a caller holding the whole trial can resolve. A schedule :func:check_decision_schedule refuses raises.

:attr:looks is the schedule as (window_s, n_samples) pairs, n_samples = round(w * sfreq) counted from the first usable sample. :meth:evaluate scores one look on a prefix of a trial, and :meth:stepper returns the incremental state machine both harnesses drive: the offline one hands it the whole trial at once, the live one hands it the growing buffer as samples arrive. The decision cannot depend on which, because each look sees only its own first n_samples.

Source code in src/ssvep/classification/decoding.py

class EarlyStopRule:
    """The growing-window early-stopping rule (#23), defined once for the offline replay
    (:func:`early_stop_decode`) and the live decoder (:class:`ssvep.runtime.online.LiveDecoder`).

    The configuration: candidate ``freqs`` (the frequencies the references are built at), the
    scoring (``method``, ``n_harmonics``, ``band``, ``line_hz``, filter-bank keywords, as
    :func:`classify` takes them), and the schedule (``min_window_s``, ``max_window_s``, ``step_s``,
    ``confidence_margin``; a ``None`` margin is the method's calibrated default). ``max_window_s``
    is required here: "the whole trial" is something only a caller holding the whole trial can
    resolve. A schedule :func:`check_decision_schedule` refuses raises.

    :attr:`looks` is the schedule as ``(window_s, n_samples)`` pairs, ``n_samples = round(w *
    sfreq)`` counted from the first usable sample. :meth:`evaluate` scores one look on a prefix of a
    trial, and :meth:`stepper` returns the incremental state machine both harnesses drive: the
    offline one hands it the whole trial at once, the live one hands it the growing buffer as
    samples arrive. The decision cannot depend on which, because each look sees only its own first
    ``n_samples``.
    """

    def __init__(self, freqs, sfreq, method="fbcca", n_harmonics=DEFAULT_N_HARMONICS, band=None,
                 line_hz=None, min_window_s=1.0, max_window_s=None, step_s=0.5,
                 confidence_margin=None, **fb_kw):
        if confidence_margin is None:
            confidence_margin = default_confidence_margin(method)
        check_decision_schedule(min_window_s, max_window_s, step_s, confidence_margin)
        if max_window_s is None:
            raise ValueError("EarlyStopRule needs a max_window_s; only a caller holding the whole "
                             "trial can read None as 'the whole trial' (early_stop_decode does)")
        self.freqs = list(freqs)
        self.sfreq = float(sfreq)
        self.method = method
        self.n_harmonics = n_harmonics
        self.band = band
        self.line_hz = line_hz
        self.fb_kw = dict(fb_kw)
        self.min_window_s = min_window_s
        self.max_window_s = float(max_window_s)
        self.step_s = step_s
        self.confidence_margin = confidence_margin
        self.looks = tuple((w, int(round(w * self.sfreq)))
                           for w in _decision_windows(min_window_s, max_window_s, step_s))
        #: Samples of usable data the last (maximum) look needs.
        self.n_max = int(round(self.max_window_s * self.sfreq))

    def summary(self) -> dict:
        """This rule as :func:`decision_rule_summary` records it."""
        return decision_rule_summary("early_stop", self.min_window_s, self.max_window_s,
                                     self.step_s, self.confidence_margin)

    def evaluate(self, prefix, i: int):
        """Score look ``i`` of :attr:`looks` on ``prefix`` (``(n_ch, >= n_samples)``; only its first
        ``n_samples`` are used). Returns a :class:`Look`, or ``None`` for a look that cannot be
        scored: under two samples, or every candidate undecodable (:class:`UndecodableWindow`).
        A skipped look is neither confident nor a decision."""
        w, n = self.looks[i]
        if n < 2:
            return None
        try:
            s = _scores(prefix[:, :n], self.freqs, self.sfreq, self.method, self.n_harmonics,
                        self.band, self.line_hz, **self.fb_kw)
        except UndecodableWindow:
            return None
        confident = (w < self.max_window_s - 1e-9
                     and is_confident(s, self.method, self.confidence_margin))
        return Look(w, n, s, decision_margin(s, self.method), confident)

    def stepper(self) -> "EarlyStopStepper":
        """A fresh per-trial state machine over this rule."""
        return EarlyStopStepper(self)

summary

summary()

This rule as :func:decision_rule_summary records it.

Source code in src/ssvep/classification/decoding.py

def summary(self) -> dict:
    """This rule as :func:`decision_rule_summary` records it."""
    return decision_rule_summary("early_stop", self.min_window_s, self.max_window_s,
                                 self.step_s, self.confidence_margin)

evaluate

evaluate(prefix, i)

Score look i of :attr:looks on prefix ((n_ch, >= n_samples); only its first n_samples are used). Returns a :class:Look, or None for a look that cannot be scored: under two samples, or every candidate undecodable (:class:UndecodableWindow). A skipped look is neither confident nor a decision.

Source code in src/ssvep/classification/decoding.py

def evaluate(self, prefix, i: int):
    """Score look ``i`` of :attr:`looks` on ``prefix`` (``(n_ch, >= n_samples)``; only its first
    ``n_samples`` are used). Returns a :class:`Look`, or ``None`` for a look that cannot be
    scored: under two samples, or every candidate undecodable (:class:`UndecodableWindow`).
    A skipped look is neither confident nor a decision."""
    w, n = self.looks[i]
    if n < 2:
        return None
    try:
        s = _scores(prefix[:, :n], self.freqs, self.sfreq, self.method, self.n_harmonics,
                    self.band, self.line_hz, **self.fb_kw)
    except UndecodableWindow:
        return None
    confident = (w < self.max_window_s - 1e-9
                 and is_confident(s, self.method, self.confidence_margin))
    return Look(w, n, s, decision_margin(s, self.method), confident)

stepper

stepper()

A fresh per-trial state machine over this rule.

Source code in src/ssvep/classification/decoding.py

def stepper(self) -> "EarlyStopStepper":
    """A fresh per-trial state machine over this rule."""
    return EarlyStopStepper(self)

EarlyStopStepper

One trial's progress through an :class:EarlyStopRule: looks are evaluated strictly in schedule order, never skipped (skipping one would change the rule's false-stop behaviour), and the first confident one ends the trial.

:attr:next_samples is how many usable samples the next pending look needs (None once :attr:done). :meth:step evaluates exactly that look on a prefix holding at least that many; :meth:feed evaluates every pending look a prefix covers. :meth:decision is the outcome so far, in :func:early_stop_decode's return shape.

Source code in src/ssvep/classification/decoding.py

class EarlyStopStepper:
    """One trial's progress through an :class:`EarlyStopRule`: looks are evaluated strictly in
    schedule order, never skipped (skipping one would change the rule's false-stop behaviour), and
    the first confident one ends the trial.

    :attr:`next_samples` is how many usable samples the next pending look needs (``None`` once
    :attr:`done`). :meth:`step` evaluates exactly that look on a prefix holding at least that many;
    :meth:`feed` evaluates every pending look a prefix covers. :meth:`decision` is the outcome so
    far, in :func:`early_stop_decode`'s return shape.
    """

    def __init__(self, rule: EarlyStopRule):
        self.rule = rule
        self._i = 0
        self.history: list[dict] = []
        self.last_scored: Look | None = None
        self.early: Look | None = None

    @property
    def done(self) -> bool:
        return self.early is not None or self._i >= len(self.rule.looks)

    @property
    def next_samples(self):
        return None if self.done else self.rule.looks[self._i][1]

    def step(self, prefix):
        """Evaluate the next pending look on ``prefix`` and return it (``None`` when it was skipped
        as unscoreable). Raises when there is no pending look, or ``prefix`` is too short for it:
        a look is never scored on less data than its window."""
        if self.done:
            raise RuntimeError("this early-stopping trial is already decided")
        n = self.rule.looks[self._i][1]
        if prefix.shape[-1] < n:
            raise ValueError(f"look {self._i} needs {n} samples; got {prefix.shape[-1]}")
        look = self.rule.evaluate(prefix, self._i)
        self._i += 1
        if look is not None:
            self.history.append({"window_s": look.window_s, "scores": look.scores.tolist(),
                                 "margin": look.margin})
            self.last_scored = look
            if look.confident:
                self.early = look
        return look

    def feed(self, prefix) -> list:
        """Evaluate every pending look ``prefix`` holds enough samples for, in order; return them
        (a skipped look as ``None``)."""
        out = []
        while not self.done and prefix.shape[-1] >= self.rule.looks[self._i][1]:
            out.append(self.step(prefix))
        return out

    @property
    def deciding_look(self):
        """The look the decision is made at so far: the early stop, else the last look scored."""
        return self.early if self.early is not None else self.last_scored

    def decision(self) -> dict:
        """``predicted_index``, ``decided_window_s``, ``decided_early`` and ``history``, as
        :func:`early_stop_decode` returns them. Provisional until :attr:`done`."""
        look = self.deciding_look
        if look is None:
            return {"predicted_index": None, "decided_window_s": None, "decided_early": False,
                    "history": list(self.history)}
        return {"predicted_index": pick(look.scores), "decided_window_s": look.window_s,
                "decided_early": self.early is not None, "history": list(self.history)}

deciding_look property

deciding_look

The look the decision is made at so far: the early stop, else the last look scored.

step

step(prefix)

Evaluate the next pending look on prefix and return it (None when it was skipped as unscoreable). Raises when there is no pending look, or prefix is too short for it: a look is never scored on less data than its window.

Source code in src/ssvep/classification/decoding.py

def step(self, prefix):
    """Evaluate the next pending look on ``prefix`` and return it (``None`` when it was skipped
    as unscoreable). Raises when there is no pending look, or ``prefix`` is too short for it:
    a look is never scored on less data than its window."""
    if self.done:
        raise RuntimeError("this early-stopping trial is already decided")
    n = self.rule.looks[self._i][1]
    if prefix.shape[-1] < n:
        raise ValueError(f"look {self._i} needs {n} samples; got {prefix.shape[-1]}")
    look = self.rule.evaluate(prefix, self._i)
    self._i += 1
    if look is not None:
        self.history.append({"window_s": look.window_s, "scores": look.scores.tolist(),
                             "margin": look.margin})
        self.last_scored = look
        if look.confident:
            self.early = look
    return look

feed

feed(prefix)

Evaluate every pending look prefix holds enough samples for, in order; return them (a skipped look as None).

Source code in src/ssvep/classification/decoding.py

def feed(self, prefix) -> list:
    """Evaluate every pending look ``prefix`` holds enough samples for, in order; return them
    (a skipped look as ``None``)."""
    out = []
    while not self.done and prefix.shape[-1] >= self.rule.looks[self._i][1]:
        out.append(self.step(prefix))
    return out

decision

decision()

predicted_index, decided_window_s, decided_early and history, as :func:early_stop_decode returns them. Provisional until :attr:done.

Source code in src/ssvep/classification/decoding.py

def decision(self) -> dict:
    """``predicted_index``, ``decided_window_s``, ``decided_early`` and ``history``, as
    :func:`early_stop_decode` returns them. Provisional until :attr:`done`."""
    look = self.deciding_look
    if look is None:
        return {"predicted_index": None, "decided_window_s": None, "decided_early": False,
                "history": list(self.history)}
    return {"predicted_index": pick(look.scores), "decided_window_s": look.window_s,
            "decided_early": self.early is not None, "history": list(self.history)}

expand_subepochs

expand_subepochs(X, y, win_samples, step_samples)

Split each trial into shorter sub-epochs; returns (Xs, ys, groups).

groups is the originating trial index so cross-validation can keep a trial's sub-epochs together (no leakage). NOTE: sub-epochs cut at arbitrary offsets have different SSVEP phase — fine for CCA, but they break TRCA's time-domain template averaging unless aligned to the stimulus cycle.

Source code in src/ssvep/classification/decoding.py

def expand_subepochs(X: np.ndarray, y: np.ndarray, win_samples: int,
                     step_samples: int):
    """Split each trial into shorter sub-epochs; returns (Xs, ys, groups).

    ``groups`` is the originating trial index so cross-validation can keep a
    trial's sub-epochs together (no leakage). NOTE: sub-epochs cut at arbitrary
    offsets have different SSVEP phase — fine for CCA, but they break TRCA's
    time-domain template averaging unless aligned to the stimulus cycle.
    """
    n_times = X.shape[-1]
    starts = range(0, n_times - win_samples + 1, step_samples)
    Xs, ys, groups = [], [], []
    for i, (trial, label) in enumerate(zip(X, y)):
        for s in starts:
            Xs.append(trial[:, s:s + win_samples]); ys.append(label); groups.append(i)
    return np.array(Xs), np.array(ys), np.array(groups)

as_mains_tuple

as_mains_tuple(line_hz)

Normalise a line_hz argument: None/0/empty -> None; a number -> (hz,); a sequence -> a sorted tuple of its positive values. Every mains-guard function accepts any form.

Source code in src/ssvep/classification/decoding.py

def as_mains_tuple(line_hz):
    """Normalise a ``line_hz`` argument: ``None``/``0``/empty -> ``None``; a number -> ``(hz,)``; a
    sequence -> a sorted tuple of its positive values. Every mains-guard function accepts any form."""
    if line_hz is None:
        return None
    vals = [line_hz] if np.isscalar(line_hz) else list(line_hz)
    out = tuple(sorted({float(v) for v in vals if v and float(v) > 0}))
    return out or None

mains_line_source

mains_line_source(line_hz)

How line_hz came to be what it is, for the record: "configured" (one frequency, as a site profile declares), "unknown" (exactly :data:UNKNOWN_MAINS_HZ, the unconfigured default), "given" (any other set a caller passed), or "off".

Source code in src/ssvep/classification/decoding.py

def mains_line_source(line_hz) -> str:
    """How ``line_hz`` came to be what it is, for the record: ``"configured"`` (one frequency, as a
    site profile declares), ``"unknown"`` (exactly :data:`UNKNOWN_MAINS_HZ`, the unconfigured
    default), ``"given"`` (any other set a caller passed), or ``"off"``."""
    lines = as_mains_tuple(line_hz)
    if lines is None:
        return "off"
    if lines == UNKNOWN_MAINS_HZ:
        return "unknown"
    return "configured" if len(lines) == 1 else "given"

single_line_hz

single_line_hz(line_hz)

One mains frequency for a display or quality check that needs exactly one number, plus whether it was assumed rather than configured: (hz, assumed).

A mains guard (:func:mains_safe_harmonics) happily takes several lines at once — guarding an extra one costs at most a harmonic order — but a PSD marker line, a mains-ratio quality check (:func:ssvep.runtime.signal_quality.line_noise_ratio) or the report's mains-ratio map each need one number to point at. line_hz's own single configured line when there is exactly one; otherwise (unconfigured, several lines, or the guard off) :data:DEFAULT_LINE_FREQUENCY_HZ with assumed=True — there is no single frequency to report a ratio against in those cases, and the long-standing 60 Hz default is the honest fallback, said out loud as assumed rather than silently presented as known.

Shared by the offline report's quality block (:func:ssvep.analysis.pipeline.run_analysis) and the live Run tab's readout (:class:ssvep.ui.run._LiveEegPlot, #362), so the mains-ratio number an operator watches live and the one the report writes afterward describe the same frequency, whichever one drew it.

Source code in src/ssvep/classification/decoding.py

def single_line_hz(line_hz) -> tuple[float, bool]:
    """One mains frequency for a display or quality check that needs exactly one number, plus
    whether it was *assumed* rather than configured: ``(hz, assumed)``.

    A mains **guard** (:func:`mains_safe_harmonics`) happily takes several lines at once — guarding
    an extra one costs at most a harmonic order — but a PSD marker line, a mains-ratio quality check
    (:func:`ssvep.runtime.signal_quality.line_noise_ratio`) or the report's mains-ratio map each need
    one number to point at. ``line_hz``'s own single configured line when there is exactly one;
    otherwise (unconfigured, several lines, or the guard off) :data:`DEFAULT_LINE_FREQUENCY_HZ` with
    ``assumed=True`` — there is no single frequency to report a ratio *against* in those cases, and
    the long-standing 60 Hz default is the honest fallback, said out loud as assumed rather than
    silently presented as known.

    Shared by the offline report's quality block (:func:`ssvep.analysis.pipeline.run_analysis`) and
    the live Run tab's readout (:class:`ssvep.ui.run._LiveEegPlot`, #362), so the mains-ratio number
    an operator watches live and the one the report writes afterward describe the same frequency,
    whichever one drew it.
    """
    lines = as_mains_tuple(line_hz)
    if lines and len(lines) == 1:
        return lines[0], False
    return DEFAULT_LINE_FREQUENCY_HZ, True

mains_lines

mains_lines(line_hz, sfreq)

Each mains frequency in line_hz and its harmonics below Nyquist ([] when off).

Power-line interference is not one line: non-linear loads on the circuit (rectifiers, switching supplies) put energy at multiples of it, and any multiple below Nyquist can reach the recording. Multiples above Nyquist are not included: the amplifier's anti-alias filter attenuates them before sampling — unlike a reference harmonic above Nyquist, which is generated digitally and aliases exactly (:func:folded_hz).

Source code in src/ssvep/classification/decoding.py

def mains_lines(line_hz, sfreq: float) -> list[float]:
    """Each mains frequency in ``line_hz`` and its harmonics below Nyquist (``[]`` when off).

    Power-line interference is not one line: non-linear loads on the circuit (rectifiers, switching
    supplies) put energy at multiples of it, and any multiple below Nyquist can reach the recording.
    Multiples above Nyquist are not included: the amplifier's anti-alias filter attenuates them
    before sampling — unlike a *reference* harmonic above Nyquist, which is generated digitally and
    aliases exactly (:func:`folded_hz`).
    """
    nyq = float(sfreq) / 2.0
    out = set()
    for line in as_mains_tuple(line_hz) or ():
        out.update(k * line for k in range(1, int(nyq // line) + 1) if k * line < nyq)
    return sorted(out)

folded_hz

folded_hz(f, sfreq)

The frequency a sinusoid at f is once sampled at sfreq: folded into [0, sfreq/2].

A reference harmonic above Nyquist is not absent — sampled, it is its alias. 3 x 43.5 Hz = 130.5 Hz at 250 Hz sampling is a 119.5 Hz sinusoid, on the 120 Hz mains harmonic (review of

302). So collisions are tested on the folded frequency.

Source code in src/ssvep/classification/decoding.py

def folded_hz(f: float, sfreq: float) -> float:
    """The frequency a sinusoid at ``f`` *is* once sampled at ``sfreq``: folded into ``[0, sfreq/2]``.

    A reference harmonic above Nyquist is not absent — sampled, it is its alias. 3 x 43.5 Hz =
    130.5 Hz at 250 Hz sampling is a 119.5 Hz sinusoid, on the 120 Hz mains harmonic (review of
    #302). So collisions are tested on the folded frequency.
    """
    r = float(f) % float(sfreq)
    return min(r, float(sfreq) - r)

mains_safe_harmonics

mains_safe_harmonics(freqs, n_harmonics, n_times, sfreq, line_hz, *, min_tol_hz=0.0)

Which reference harmonics to use when the recording carries mains at line_hz (#280).

Returns (orders, undecodable): the harmonic orders (1-based) every candidate's reference bank uses, and the indices of candidates left with none. line_hz is one frequency, several (:data:UNKNOWN_MAINS_HZ when the site has not said), or None to disable the guard.

Why. Mains is a large, narrow, stationary line in unfiltered EEG — far larger than any SSVEP in the band it occupies. A candidate whose reference bank includes a sinusoid at that frequency matches the mains line on every trial, whatever the participant was looking at, and wins. That is the single-class collapse #280 found: the teaching example's 10 Hz reversing target is decoded at 2f = 20 Hz, whose third harmonic is exactly 60 Hz. Plain flicker hits it the same way (20 Hz × 3, 30 Hz × 2); only which targets collide changes.

Which harmonics count as "on" a line. Over a window of T seconds a reference sinusoid cannot be told apart from a line closer than the window's own frequency resolution, 1/T — the main lobe of its spectral kernel. So a harmonic whose sampled frequency (:func:folded_hz: a harmonic above Nyquist is its alias) lies within max(1/T, MAINS_MIN_GUARD_HZ) of a mains line or one of its harmonics below Nyquist (:func:mains_lines) collides. The tolerance is the window's: 2 Hz for a 0.5 s decision, 0.25 Hz for a 4 s one. Because it is 1/T, a harmonic exactly at the boundary can flip when the window gains or loses a sample, which is why the orders used are recorded per window (per live trial, and per offline sweep window), not once.

What is dropped — the order, from every candidate. Dropping a colliding harmonic from only the candidate that has it leaves that candidate with fewer reference terms than the rest, and a canonical correlation over fewer terms runs lower on the same noise: that candidate is then under-picked, the mirror image of the bias being removed (measured, and pinned by a test). So a colliding order is removed from all candidates, and every candidate keeps the same number of reference terms.

Except the fundamental. A candidate whose own fundamental is on a mains line cannot be decoded at all — every harmonic of it is on a mains harmonic too. It gets no references (its score is NaN and it is never picked, which is the honest result for a target sitting on the power line), and it does not strip the fundamental or harmonics from anyone else. If every candidate is undecodable the scorers refuse (:func:pick raises :class:UndecodableWindow) rather than name one.

A window shorter than :data:MIN_GUARD_SAMPLES raises :class:UndecodableWindow too: its tolerance is so wide that the answer would be noise.

min_tol_hz widens the tolerance for a caller whose estimator reaches further than the window's resolution. The decoder never passes it; the target-SNR DV does (:func:ssvep.classification.spectral.snr_safe_harmonics, #309), because a line anywhere in an SNR bin's noise neighbourhood corrupts that bin's SNR too.

Source code in src/ssvep/classification/decoding.py

def mains_safe_harmonics(freqs, n_harmonics: int, n_times: int, sfreq: float, line_hz, *,
                         min_tol_hz: float = 0.0):
    """Which reference harmonics to use when the recording carries mains at ``line_hz`` (#280).

    Returns ``(orders, undecodable)``: the harmonic orders (1-based) every candidate's reference
    bank uses, and the indices of candidates left with none. ``line_hz`` is one frequency, several
    (:data:`UNKNOWN_MAINS_HZ` when the site has not said), or ``None`` to disable the guard.

    **Why.** Mains is a large, narrow, stationary line in unfiltered EEG — far larger than any SSVEP
    in the band it occupies. A candidate whose reference bank includes a sinusoid at that frequency
    matches the mains line on every trial, whatever the participant was looking at, and wins. That
    is the single-class collapse #280 found: the teaching example's 10 Hz reversing target is
    decoded at 2f = 20 Hz, whose third harmonic is exactly 60 Hz. Plain flicker hits it the same way
    (20 Hz × 3, 30 Hz × 2); only which targets collide changes.

    **Which harmonics count as "on" a line.** Over a window of ``T`` seconds a reference sinusoid
    cannot be told apart from a line closer than the window's own frequency resolution, ``1/T`` —
    the main lobe of its spectral kernel. So a harmonic whose **sampled** frequency
    (:func:`folded_hz`: a harmonic above Nyquist is its alias) lies within
    ``max(1/T, MAINS_MIN_GUARD_HZ)`` of a mains line or one of its harmonics below Nyquist
    (:func:`mains_lines`) collides. The tolerance is the window's: 2 Hz for a 0.5 s decision,
    0.25 Hz for a 4 s one. Because it is ``1/T``, a harmonic exactly at the boundary can flip when
    the window gains or loses a sample, which is why the orders used are recorded per window (per
    live trial, and per offline sweep window), not once.

    **What is dropped — the order, from every candidate.** Dropping a colliding harmonic from only
    the candidate that has it leaves that candidate with fewer reference terms than the rest, and a
    canonical correlation over fewer terms runs lower on the same noise: that candidate is then
    under-picked, the mirror image of the bias being removed (measured, and pinned by a test). So
    a colliding order is removed from **all** candidates, and every candidate keeps the same number
    of reference terms.

    **Except the fundamental.** A candidate whose own fundamental is on a mains line cannot be
    decoded at all — every harmonic of it is on a mains harmonic too. It gets no references (its
    score is NaN and it is never picked, which is the honest result for a target sitting on the
    power line), and it does not strip the fundamental or harmonics from anyone else. If *every*
    candidate is undecodable the scorers refuse (:func:`pick` raises :class:`UndecodableWindow`)
    rather than name one.

    A window shorter than :data:`MIN_GUARD_SAMPLES` raises :class:`UndecodableWindow` too: its
    tolerance is so wide that the answer would be noise.

    ``min_tol_hz`` widens the tolerance for a caller whose estimator reaches further than the
    window's resolution. The decoder never passes it; the target-SNR DV does
    (:func:`ssvep.classification.spectral.snr_safe_harmonics`, #309), because a line anywhere in an
    SNR bin's noise neighbourhood corrupts that bin's SNR too.
    """
    orders = list(range(1, int(n_harmonics) + 1))
    freqs = [float(f) for f in freqs]
    lines = mains_lines(line_hz, sfreq)
    if not lines or not freqs:
        return tuple(orders), ()
    if n_times < MIN_GUARD_SAMPLES:
        raise UndecodableWindow(f"mains guard: a {n_times}-sample window is too short to resolve "
                                f"which harmonics collide with mains (minimum {MIN_GUARD_SAMPLES})")
    tol = max(float(sfreq) / float(n_times), MAINS_MIN_GUARD_HZ, float(min_tol_hz or 0.0))

    def on_mains(f, h):
        fa = folded_hz(h * f, sfreq)
        return any(abs(fa - line) < tol for line in lines)

    undecodable = tuple(i for i, f in enumerate(freqs) if on_mains(f, 1))
    rest = [f for i, f in enumerate(freqs) if i not in undecodable]
    dropped = {h for h in orders if h > 1 for f in rest if on_mains(f, h)}
    keep = tuple(h for h in orders if h not in dropped)
    return keep, undecodable

mains_guard_summary

mains_guard_summary(freqs, n_harmonics, n_times, sfreq, line_hz)

What :func:mains_safe_harmonics did at one window, as a JSON-able record.

Source code in src/ssvep/classification/decoding.py

def mains_guard_summary(freqs, n_harmonics: int, n_times: int, sfreq: float, line_hz) -> dict:
    """What :func:`mains_safe_harmonics` did at one window, as a JSON-able record."""
    orders, undecodable = mains_safe_harmonics(freqs, n_harmonics, n_times, sfreq, line_hz)
    lines = as_mains_tuple(line_hz)
    return {
        "line_hz": list(lines) if lines else None,
        "line_hz_source": mains_line_source(line_hz),
        "window_s": round(float(n_times) / float(sfreq), 4),
        "harmonics_used": list(orders),
        "harmonics_dropped": [h for h in range(1, int(n_harmonics) + 1) if h not in orders],
        "undecodable_hz": [float(freqs[i]) for i in undecodable],
    }

pick

pick(scores)

Index of the best candidate score, skipping undecodable (NaN) candidates.

Raises :class:UndecodableWindow when no candidate has a score, instead of letting argmax return 0 for an all-NaN vector — a silent, confident, wrong answer.

Source code in src/ssvep/classification/decoding.py

def pick(scores) -> int:
    """Index of the best candidate score, skipping undecodable (NaN) candidates.

    Raises :class:`UndecodableWindow` when no candidate has a score, instead of letting ``argmax``
    return 0 for an all-NaN vector — a silent, confident, wrong answer.
    """
    s = np.asarray(scores, dtype=float)
    if not np.isfinite(s).any():
        raise UndecodableWindow("no candidate can be decoded at this window: every one sits on a "
                                "mains line at the window's frequency resolution")
    return int(np.nanargmax(np.where(np.isfinite(s), s, np.nan)))

alias_safe_orders

alias_safe_orders(freq, orders, sfreq)

Which of orders (1-based harmonic orders) are real references for freq, not aliases.

See :func:_is_alias for the boundary. One candidate's own answer — used where only one candidate's harmonics are in question; :func:nyquist_safe_harmonics is the shared-bank answer for a whole set of candidates decoded together.

Source code in src/ssvep/classification/decoding.py

def alias_safe_orders(freq: float, orders, sfreq: float) -> tuple:
    """Which of ``orders`` (1-based harmonic orders) are real references for ``freq``, not aliases.

    See :func:`_is_alias` for the boundary. One candidate's own answer — used where only one
    candidate's harmonics are in question; :func:`nyquist_safe_harmonics` is the shared-bank answer
    for a whole set of candidates decoded together.
    """
    return tuple(h for h in orders if not _is_alias(freq, h, sfreq))

nyquist_safe_harmonics

nyquist_safe_harmonics(freqs, orders, sfreq)

Shared harmonic orders after dropping any that alias at or above Nyquist for any decodable candidate (#310; #346 review finding 1).

Returns (per_freq_orders, undecodable): a tuple aligned to freqs (each entry the same shared subset of orders for every decodable candidate, or () for an undecodable one), and the indices left with none.

Why shared, not per candidate — the mains guard's own rule, applied here too. An earlier version of this function computed each candidate's surviving orders independently: a 44 Hz candidate's 3rd harmonic (132 Hz at 250 Hz sampling) aliases past the 125 Hz Nyquist edge while a 10 Hz candidate's does not, so the 44 Hz candidate kept fewer reference terms than the rest of the bank. That is exactly the bias :func:mains_safe_harmonics exists to prevent (see its own docstring, and DESIGN_PRINCIPLES #12): a candidate with fewer reference columns runs a lower canonical correlation on the same noise, so it is under-picked on null data and mis-decoded at low SNR — measured on the shipped high-band preset (review of #346) at a cost bigger than the alias bias #310 set out to remove. So an order that aliases for any decodable candidate is dropped from every decodable candidate, mirroring :func:mains_safe_harmonics's own rule for a colliding order, and every decodable candidate keeps the same reference dimensionality. A candidate whose own fundamental already aliases (h=1) has no real reference at any order — it was never going to be decodable, with or without this rule — and is named in undecodable rather than silently given an empty bank.

Source code in src/ssvep/classification/decoding.py

def nyquist_safe_harmonics(freqs, orders, sfreq: float):
    """Shared harmonic orders after dropping any that alias at or above Nyquist for *any* decodable
    candidate (#310; #346 review finding 1).

    Returns ``(per_freq_orders, undecodable)``: a tuple aligned to ``freqs`` (each entry the same
    shared subset of ``orders`` for every decodable candidate, or ``()`` for an undecodable one),
    and the indices left with none.

    **Why shared, not per candidate — the mains guard's own rule, applied here too.** An earlier
    version of this function computed each candidate's surviving orders independently: a 44 Hz
    candidate's 3rd harmonic (132 Hz at 250 Hz sampling) aliases past the 125 Hz Nyquist edge while
    a 10 Hz candidate's does not, so the 44 Hz candidate kept fewer reference terms than the rest of
    the bank. That is exactly the bias :func:`mains_safe_harmonics` exists to prevent (see its own
    docstring, and DESIGN_PRINCIPLES #12): a candidate with fewer reference columns runs a lower
    canonical correlation on the same noise, so it is under-picked on null data and mis-decoded at
    low SNR — measured on the shipped high-band preset (review of #346) at a cost bigger than the
    alias bias #310 set out to remove. So an order that aliases for *any* decodable candidate is
    dropped from *every* decodable candidate, mirroring :func:`mains_safe_harmonics`'s own rule for
    a colliding order, and every decodable candidate keeps the same reference dimensionality. A
    candidate whose own fundamental already aliases (``h=1``) has no real reference at any order —
    it was never going to be decodable, with or without this rule — and is named in
    ``undecodable`` rather than silently given an empty bank.
    """
    orders = tuple(orders)
    freqs_f = [float(f) for f in freqs]
    undecodable = tuple(i for i, f in enumerate(freqs_f) if _is_alias(f, 1, sfreq))
    rest = [f for i, f in enumerate(freqs_f) if i not in undecodable]
    dropped = {h for h in orders if h > 1 for f in rest if _is_alias(f, h, sfreq)}
    keep = tuple(h for h in orders if h not in dropped)
    per_freq = tuple(keep if i not in undecodable else () for i in range(len(freqs_f)))
    return per_freq, undecodable

usable_harmonics

usable_harmonics(freqs, n_harmonics, n_times, sfreq, line_hz)

Shared harmonic orders after both guards, for one window (#310; #346 review).

Returns (per_freq_orders, undecodable) in :func:nyquist_safe_harmonics's shape: every decodable candidate gets the same tuple of surviving orders, and an undecodable one gets () directly (not a leftover alias-filtered remnant of the mains-only set — the record and the decode must agree on this without a second rewrite, review of #346 finding 2).

Not a simple two-step compose. Chaining :func:mains_safe_harmonics then :func:nyquist_safe_harmonics tests a candidate's aliased harmonic for a mains collision on its folded frequency, even though that harmonic was never going to be built for that candidate at all — which strips a real harmonic from every other candidate over a collision that belongs to a reference nobody uses (review of #346 finding 4: a 43.5 Hz candidate's aliased 3rd harmonic folds onto 120 Hz mains and, chained naively, takes the real 108/114 Hz 3rd harmonics of 36/38 Hz candidates with it). So a colliding-with-mains test here only ever asks it of an order that is alias-safe for that candidate to begin with; an order that aliases is dropped for aliasing, not re-litigated for mains.

Both guards drop globally across the decodable candidates (mains: one external interferer shared by the whole bank; aliasing: symmetric reference dimensionality, see :func:nyquist_safe_harmonics), so this reduces to one shared order set, computed together rather than as two independent shared sets intersected afterwards.

This is what :func:cca_scores actually builds each candidate's reference bank from. A record that only carries the mains-only harmonics_used (:func:mains_guard_summary, computed without knowledge of aliasing) is not describing what was decoded whenever some candidate's own frequency puts a harmonic at or above Nyquist — use this (or :func:usable_harmonics_summary for the JSON form) wherever the record needs to match the decode exactly.

Source code in src/ssvep/classification/decoding.py

def usable_harmonics(freqs, n_harmonics: int, n_times: int, sfreq: float, line_hz):
    """Shared harmonic orders after **both** guards, for one window (#310; #346 review).

    Returns ``(per_freq_orders, undecodable)`` in :func:`nyquist_safe_harmonics`'s shape: every
    decodable candidate gets the *same* tuple of surviving orders, and an undecodable one gets
    ``()`` directly (not a leftover alias-filtered remnant of the mains-only set — the record and
    the decode must agree on this without a second rewrite, review of #346 finding 2).

    **Not a simple two-step compose.** Chaining :func:`mains_safe_harmonics` then
    :func:`nyquist_safe_harmonics` tests a candidate's *aliased* harmonic for a mains collision on
    its folded frequency, even though that harmonic was never going to be built for that candidate
    at all — which strips a *real* harmonic from every other candidate over a collision that
    belongs to a reference nobody uses (review of #346 finding 4: a 43.5 Hz candidate's aliased 3rd
    harmonic folds onto 120 Hz mains and, chained naively, takes the real 108/114 Hz 3rd harmonics
    of 36/38 Hz candidates with it). So a colliding-with-mains test here only ever asks it of an
    order that is alias-safe for that candidate to begin with; an order that aliases is dropped for
    aliasing, not re-litigated for mains.

    Both guards drop **globally** across the decodable candidates (mains: one external interferer
    shared by the whole bank; aliasing: symmetric reference dimensionality, see
    :func:`nyquist_safe_harmonics`), so this reduces to one shared order set, computed together
    rather than as two independent shared sets intersected afterwards.

    This is what :func:`cca_scores` actually builds each candidate's reference bank from. A record
    that only carries the mains-only ``harmonics_used`` (:func:`mains_guard_summary`, computed
    without knowledge of aliasing) is not describing what was decoded whenever some candidate's own
    frequency puts a harmonic at or above Nyquist — use this (or :func:`usable_harmonics_summary`
    for the JSON form) wherever the record needs to match the decode exactly.
    """
    orders = tuple(range(1, int(n_harmonics) + 1))
    freqs_f = [float(f) for f in freqs]
    lines = mains_lines(line_hz, sfreq)
    if lines and freqs_f and n_times < MIN_GUARD_SAMPLES:
        raise UndecodableWindow(f"mains guard: a {n_times}-sample window is too short to resolve "
                                f"which harmonics collide with mains (minimum {MIN_GUARD_SAMPLES})")
    tol = max(float(sfreq) / float(n_times), MAINS_MIN_GUARD_HZ) if lines else None

    def on_mains(f, h):
        # Only a real reference can collide with the power line: an order that already aliases for
        # this candidate was never built, so its folded frequency says nothing about what the
        # decoder actually references (review of #346 finding 4).
        if not lines or _is_alias(f, h, sfreq):
            return False
        fa = folded_hz(h * f, sfreq)
        return any(abs(fa - line) < tol for line in lines)

    undecodable = tuple(i for i, f in enumerate(freqs_f)
                        if _is_alias(f, 1, sfreq) or on_mains(f, 1))
    rest = [f for i, f in enumerate(freqs_f) if i not in undecodable]
    dropped = {h for h in orders if h > 1
              for f in rest if _is_alias(f, h, sfreq) or on_mains(f, h)}
    keep = tuple(h for h in orders if h not in dropped)
    per_freq = tuple(keep if i not in undecodable else () for i in range(len(freqs_f)))
    return per_freq, undecodable

usable_harmonics_summary

usable_harmonics_summary(freqs, n_harmonics, n_times, sfreq, line_hz)

What :func:usable_harmonics did at one window, as a JSON-able record (#310).

harmonics_used is a list per candidate, aligned to freqs — unlike the mains guard's own harmonics_used (:func:mains_guard_summary), which is computed without knowledge of aliasing, this one is folded together with the Nyquist guard. Every decodable candidate shares the same orders (:func:usable_harmonics); an undecodable one is [] directly, from usable_harmonics itself — nothing here rewrites it a second time (review of #346 finding 2). harmonics_reduced_hz names every candidate left with fewer orders than the mains guard's own (alias-unaware) shared set allowed, so a run that lost harmonics to aliasing says so in the record even for a candidate whose own frequency never aliased (its bank was still cut down to match the rest of the decodable set). undecodable_hz names any candidate left with no usable harmonic at all — refused (NaN score, never picked by :func:pick; see :func:cca_scores) rather than silently decoded on no reference at all.

Source code in src/ssvep/classification/decoding.py

def usable_harmonics_summary(freqs, n_harmonics: int, n_times: int, sfreq: float, line_hz) -> dict:
    """What :func:`usable_harmonics` did at one window, as a JSON-able record (#310).

    ``harmonics_used`` is a list **per candidate**, aligned to ``freqs`` — unlike the mains guard's
    own ``harmonics_used`` (:func:`mains_guard_summary`), which is computed without knowledge of
    aliasing, this one is folded together with the Nyquist guard. Every decodable candidate shares
    the same orders (:func:`usable_harmonics`); an undecodable one is ``[]`` directly, from
    ``usable_harmonics`` itself — nothing here rewrites it a second time (review of #346 finding 2).
    ``harmonics_reduced_hz`` names every candidate left with fewer orders than the mains guard's own
    (alias-unaware) shared set allowed, so a run that lost harmonics to aliasing says so in the
    record even for a candidate whose own frequency never aliased (its bank was still cut down to
    match the rest of the decodable set). ``undecodable_hz`` names any candidate left with *no*
    usable harmonic at all — refused (NaN score, never picked by :func:`pick`; see
    :func:`cca_scores`) rather than silently decoded on no reference at all.
    """
    freqs_f = [float(f) for f in freqs]
    mains_orders, _ = mains_safe_harmonics(freqs, n_harmonics, n_times, sfreq, line_hz)
    per_freq, undecodable = usable_harmonics(freqs, n_harmonics, n_times, sfreq, line_hz)
    return {
        "harmonics_used": [list(o) for o in per_freq],
        "harmonics_reduced_hz": [freqs_f[i] for i, o in enumerate(per_freq)
                                 if i not in undecodable and len(o) < len(mains_orders)],
        "undecodable_hz": [freqs_f[i] for i in undecodable],
    }

reference_signals

reference_signals(freq, n_harmonics, n_times, sfreq, harmonics=None)

Sin/cos reference matrix (n_times, 2*len(harmonics)) for one frequency.

harmonics is an explicit list of 1-based orders; None means 1..n_harmonics. An empty harmonics (()) raises rather than building a matrix with no columns: a candidate with no usable harmonic is undecodable (:func:usable_harmonics) and :func:cca_scores never reaches this function for it — an empty bank reaching here is a caller bug, and np.column_stack([]) raising ValueError: need at least one array to concatenate deep inside numpy is not a useful way to learn that (review of #346, finding 8).

Source code in src/ssvep/classification/decoding.py

def reference_signals(freq: float, n_harmonics: int, n_times: int,
                      sfreq: float, harmonics=None) -> np.ndarray:
    """Sin/cos reference matrix ``(n_times, 2*len(harmonics))`` for one frequency.

    ``harmonics`` is an explicit list of 1-based orders; ``None`` means ``1..n_harmonics``. An
    empty ``harmonics`` (``()``) raises rather than building a matrix with no columns: a candidate
    with no usable harmonic is undecodable (:func:`usable_harmonics`) and :func:`cca_scores` never
    reaches this function for it — an empty bank reaching here is a caller bug, and
    ``np.column_stack([])`` raising ``ValueError: need at least one array to concatenate`` deep
    inside numpy is not a useful way to learn that (review of #346, finding 8).
    """
    hs = tuple(harmonics) if harmonics is not None else tuple(range(1, n_harmonics + 1))
    if not hs:
        raise ValueError("reference_signals: harmonics is empty — there is no reference to build; "
                         "an undecodable candidate (usable_harmonics) should never reach here")
    t = np.arange(n_times) / sfreq
    cols = []
    for h in hs:
        cols.append(np.sin(2 * np.pi * h * freq * t))
        cols.append(np.cos(2 * np.pi * h * freq * t))
    return np.column_stack(cols)

cca_scores

cca_scores(trial, freqs, sfreq, n_harmonics=DEFAULT_N_HARMONICS, line_hz=None)

Max canonical correlation of one trial (n_ch, n_times) to each freq.

line_hz (the mains frequency or frequencies, e.g. 60.0) leaves out reference harmonics that would land on the power line — see :func:mains_safe_harmonics. Independently, any harmonic that aliases at or above Nyquist for a candidate's own frequency is left out of that candidate's bank alone — see :func:usable_harmonics (#310), which composes the two guards. A candidate left with no reference at all scores NaN; choose with :func:pick, not argmax. None (the default) builds every harmonic, as this function always has.

Source code in src/ssvep/classification/decoding.py

def cca_scores(trial: np.ndarray, freqs, sfreq: float,
               n_harmonics: int = DEFAULT_N_HARMONICS, line_hz=None) -> np.ndarray:
    """Max canonical correlation of one trial ``(n_ch, n_times)`` to each freq.

    ``line_hz`` (the mains frequency or frequencies, e.g. 60.0) leaves out reference harmonics that
    would land on the power line — see :func:`mains_safe_harmonics`. Independently, any harmonic
    that aliases at or above Nyquist for a *candidate's own* frequency is left out of that
    candidate's bank alone — see :func:`usable_harmonics` (#310), which composes the two guards. A
    candidate left with no reference at all scores NaN; choose with :func:`pick`, not ``argmax``.
    ``None`` (the default) builds every harmonic, as this function always has.
    """
    X = trial.T
    n_times = X.shape[0]
    per_freq, undecodable = usable_harmonics(freqs, n_harmonics, n_times, sfreq, line_hz)
    return np.array([
        np.nan if i in undecodable else
        _canon_corr_max(X, reference_signals(f, n_harmonics, n_times, sfreq, harmonics=per_freq[i]))
        for i, f in enumerate(freqs)])

fbcca_filter_bank

fbcca_filter_bank(band_low, band_high, sfreq, n_subbands=DEFAULT_N_SUBBANDS, a=DEFAULT_FB_A, b=DEFAULT_FB_B)

The filter bank FBCCA decodes through: [(passband, weight), ...], one per sub-band.

In order: sub-band k asks for k*band_low - 2 to k*band_high + 2 Hz, with the upper edge clipped to Nyquist - 1. It is dropped when that leaves it 2 Hz wide or less, measured from the unclamped lower edge, so a high candidate band or a low sampling rate gets fewer than n_subbands. If none survives, the single band band_low - 2 to band_high + 2 (upper edge clipped the same way) is used instead. Only then is the lower edge raised to 0.5 Hz (:func:_passband): passband is the (low, high) the band-pass is built with, or None when that leaves nothing, and the sub-band's data is used unfiltered. :func:fbcca_scores iterates exactly this list and :func:filter_bank_summary records it, so the decode and its record cannot disagree.

Source code in src/ssvep/classification/decoding.py

def fbcca_filter_bank(band_low, band_high, sfreq, n_subbands=DEFAULT_N_SUBBANDS, a=DEFAULT_FB_A,
                      b=DEFAULT_FB_B) -> list:
    """The filter bank FBCCA decodes through: ``[(passband, weight), ...]``, one per sub-band.

    In order: sub-band *k* asks for ``k*band_low - 2`` to ``k*band_high + 2`` Hz, with the upper
    edge clipped to Nyquist - 1. It is dropped when that leaves it 2 Hz wide or less, measured from
    the **unclamped** lower edge, so a high candidate band or a low sampling rate gets fewer than
    ``n_subbands``. If none survives, the single band ``band_low - 2`` to ``band_high + 2`` (upper
    edge clipped the same way) is used instead. Only then is the lower edge raised to 0.5 Hz
    (:func:`_passband`): ``passband`` is the ``(low, high)`` the band-pass is built with, or
    ``None`` when that leaves nothing, and the sub-band's data is used unfiltered.
    :func:`fbcca_scores` iterates exactly this list and :func:`filter_bank_summary` records it,
    so the decode and its record cannot disagree.
    """
    return [(_passband(lo, hi, sfreq), k ** (-a) + b)
            for k, (lo, hi) in enumerate(_fb_subbands(band_low, band_high, n_subbands, sfreq), 1)]

filter_bank_summary

filter_bank_summary(method, band, sfreq, n_subbands=DEFAULT_N_SUBBANDS, a=DEFAULT_FB_A, b=DEFAULT_FB_B)

The sub-band configuration a decode ran with, as a JSON-able record (filter_bank, #22).

Built from :func:fbcca_filter_bank, the list :func:fbcca_scores filters through, with the same keyword defaults, so pass it exactly the method, band and filter-bank keywords the decode was given. n_subbands in the record is how many sub-bands were actually filtered; n_subbands_requested is the setting. bands_hz are the passbands really built (None for one used unfiltered), weights their weights, and band_hz the candidate range they scale from. Plain CCA filters nothing, so for it every field but method is None.

The record is window-independent: the edges depend only on the band and the sampling rate, never on the decision window, so one record describes every window a harness decodes with that band (the pipeline's whole sweep, its operating window and its channel search; every trial of a live run).

Scope: this is FBCCA's (and plain CCA's) record only. TRCA is band-passed separately and is not described here at all — see :func:ssvep.analysis.pipeline.run_trca's own band_hz (the edges actually built, post-clamp since #325) and band_hz_requested.

The setting (n_subbands_requested, a, b) is what says two decodes were configured alike. The edges follow from it, from the candidate frequencies and from the sampling rate, so a contrast-reversing run decoded at 2f legitimately has different edges from a solid one decoded at f.

Source code in src/ssvep/classification/decoding.py

def filter_bank_summary(method, band, sfreq, n_subbands=DEFAULT_N_SUBBANDS, a=DEFAULT_FB_A,
                        b=DEFAULT_FB_B) -> dict:
    """The sub-band configuration a decode ran with, as a JSON-able record (``filter_bank``, #22).

    Built from :func:`fbcca_filter_bank`, the list :func:`fbcca_scores` filters through, with the
    same keyword defaults, so pass it exactly the ``method``, ``band`` and filter-bank keywords the
    decode was given. ``n_subbands`` in the record is how many sub-bands were actually filtered;
    ``n_subbands_requested`` is the setting. ``bands_hz`` are the passbands really built
    (``None`` for one used unfiltered), ``weights`` their weights, and ``band_hz`` the candidate
    range they scale from. Plain CCA filters nothing, so for it every field but ``method`` is
    ``None``.

    The record is **window-independent**: the edges depend only on the band and the sampling
    rate, never on the decision window, so one record describes every window a harness decodes
    with that band (the pipeline's whole sweep, its operating window and its channel search; every
    trial of a live run).

    **Scope: this is FBCCA's (and plain CCA's) record only.** TRCA is band-passed separately and
    is not described here at all — see :func:`ssvep.analysis.pipeline.run_trca`'s own ``band_hz``
    (the edges actually built, post-clamp since #325) and ``band_hz_requested``.

    The setting (``n_subbands_requested``, ``a``, ``b``) is what says two decodes were configured
    alike. The edges follow from it, from the candidate frequencies and from the sampling rate, so
    a contrast-reversing run decoded at 2f legitimately has different edges from a solid one
    decoded at f.
    """
    none = {"method": method, "n_subbands": None, "n_subbands_requested": None, "bands_hz": None,
            "weights": None, "a": None, "b": None, "band_hz": None}
    if method != "fbcca":
        return none
    rec = dict(none, n_subbands_requested=int(n_subbands), a=float(a), b=float(b))
    if band is None or not sfreq:
        return rec
    bank = fbcca_filter_bank(band[0], band[1], float(sfreq), n_subbands, a, b)
    rec.update(n_subbands=len(bank),
               bands_hz=[None if e is None else [float(e[0]), float(e[1])] for e, _ in bank],
               weights=[float(w) for _, w in bank],
               band_hz=[float(band[0]), float(band[1])])
    return rec

fbcca_scores

fbcca_scores(trial, freqs, sfreq, band_low, band_high, n_harmonics=DEFAULT_N_HARMONICS, n_subbands=DEFAULT_N_SUBBANDS, a=DEFAULT_FB_A, b=DEFAULT_FB_B, line_hz=None)

Filter-bank CCA combined scores (Chen et al. 2015), one per candidate freq.

The sub-bands and their weights are :func:fbcca_filter_bank's; a harness that records the configuration does it with :func:filter_bank_summary, which reads the same list.

line_hz is passed to :func:cca_scores for every sub-band — see :func:mains_safe_harmonics. The sub-band that contains the mains line is where a colliding reference does its damage, so the guard has to reach every one of them. The Nyquist-alias guard (:func:usable_harmonics, #310) reaches every sub-band the same way, for the same reason: each call builds the same candidate's reference bank at the same sfreq, so which harmonics alias does not depend on which sub-band is currently filtered.

Source code in src/ssvep/classification/decoding.py

def fbcca_scores(trial, freqs, sfreq, band_low, band_high, n_harmonics=DEFAULT_N_HARMONICS,
                 n_subbands=DEFAULT_N_SUBBANDS, a=DEFAULT_FB_A, b=DEFAULT_FB_B,
                 line_hz=None) -> np.ndarray:
    """Filter-bank CCA combined scores (Chen et al. 2015), one per candidate freq.

    The sub-bands and their weights are :func:`fbcca_filter_bank`'s; a harness that records the
    configuration does it with :func:`filter_bank_summary`, which reads the same list.

    ``line_hz`` is passed to :func:`cca_scores` for every sub-band — see
    :func:`mains_safe_harmonics`. The sub-band that contains the mains line is where a colliding
    reference does its damage, so the guard has to reach every one of them. The Nyquist-alias guard
    (:func:`usable_harmonics`, #310) reaches every sub-band the same way, for the same reason: each
    call builds the same candidate's reference bank at the same ``sfreq``, so which harmonics alias
    does not depend on which sub-band is currently filtered.
    """
    total = np.zeros(len(freqs))
    for edges, w in fbcca_filter_bank(band_low, band_high, sfreq, n_subbands, a, b):
        x = trial if edges is None else _bandpass(trial, edges[0], edges[1], sfreq)
        total += w * cca_scores(x, freqs, sfreq, n_harmonics, line_hz=line_hz) ** 2
    return total

trca_feasible

trca_feasible(y, min_trials_per_class=4)

True if every class has >= min_trials_per_class trials (TRCA pre-check).

Source code in src/ssvep/classification/decoding.py

def trca_feasible(y, min_trials_per_class: int = 4) -> bool:
    """True if every class has >= min_trials_per_class trials (TRCA pre-check)."""
    _, counts = np.unique(y, return_counts=True)
    return len(counts) >= 2 and counts.min() >= min_trials_per_class

classify

classify(trial, freqs, sfreq, method='fbcca', n_harmonics=DEFAULT_N_HARMONICS, band=None, line_hz=None, **fb_kw)

Predict the attended-frequency index for one trial. Calibration-free.

method='cca' or 'fbcca' (fbcca needs band=(low, high)). line_hz guards the references against the mains line (:func:mains_safe_harmonics); None disables it. Reference harmonics that alias at or above Nyquist for any decodable candidate are always left out of every decodable candidate (:func:usable_harmonics, #310, #346) — that guard has no on/off switch, because a harmonic that samples to a frequency the design never declared is never a legitimate reference. Raises :class:UndecodableWindow when no candidate can be scored.

Source code in src/ssvep/classification/decoding.py

def classify(trial, freqs, sfreq, method="fbcca", n_harmonics=DEFAULT_N_HARMONICS,
             band=None, line_hz=None, **fb_kw) -> int:
    """Predict the attended-frequency index for one trial. Calibration-free.

    method='cca' or 'fbcca' (fbcca needs band=(low, high)). ``line_hz`` guards the references
    against the mains line (:func:`mains_safe_harmonics`); ``None`` disables it. Reference harmonics
    that alias at or above Nyquist for any decodable candidate are always left out of every
    decodable candidate (:func:`usable_harmonics`, #310, #346) — that guard has no on/off switch,
    because a harmonic that samples to a frequency the design never declared is never a legitimate
    reference. Raises :class:`UndecodableWindow` when no candidate can be scored.
    """
    return pick(_scores(trial, freqs, sfreq, method, n_harmonics, band, line_hz, **fb_kw))

default_confidence_margin

default_confidence_margin(method)

:data:DEFAULT_CONFIDENCE_MARGINS for method; raises for a method with none.

Source code in src/ssvep/classification/decoding.py

def default_confidence_margin(method: str) -> float:
    """:data:`DEFAULT_CONFIDENCE_MARGINS` for ``method``; raises for a method with none."""
    try:
        return DEFAULT_CONFIDENCE_MARGINS[method]
    except KeyError:
        raise ValueError(f"unknown calibration-free method: {method}") from None

decision_margin

decision_margin(scores, method)

:data:DECISION_STATISTIC for one window's candidate scores from method: the gap between the best and second-best decodable (finite) candidate on the correlation scale, relative to the best, in [0, 1]. Scale-free: a ratio of two correlations, so it does not grow or shrink with how large chance correlations are at a given window length or channel count.

Fewer than two decodable candidates read as 0.0, no evidence: a forced choice (the mains guard left one candidate, or the design has one target) was not chosen over anything, so it is never "confident". Never raises on the scores: an early-stopping check that could raise on the very data it is deciding about would be worse than the fixed window it replaces.

Source code in src/ssvep/classification/decoding.py

def decision_margin(scores, method: str) -> float:
    """:data:`DECISION_STATISTIC` for one window's candidate ``scores`` from ``method``: the gap
    between the best and second-best *decodable* (finite) candidate on the correlation scale,
    relative to the best, in ``[0, 1]``. Scale-free: a ratio of two correlations, so it does not
    grow or shrink with how large chance correlations are at a given window length or channel
    count.

    Fewer than two decodable candidates read as ``0.0``, no evidence: a forced choice (the mains
    guard left one candidate, or the design has one target) was not chosen over anything, so it is
    never "confident". Never raises on the scores: an early-stopping check that could raise on the
    very data it is deciding about would be worse than the fixed window it replaces.
    """
    r = _correlation_scale(scores, method)
    finite = np.sort(r[np.isfinite(r)])
    if finite.size < 2:
        return 0.0
    top, second = finite[-1], finite[-2]
    if top <= 0:
        return 0.0
    return float((top - second) / top)

is_confident

is_confident(scores, method, margin=None)

Whether scores clears :func:decision_margin at margin (None: the method's calibrated :func:default_confidence_margin): the one predicate the offline replay (:func:early_stop_decode) calls to decide "stop now" versus "wait for more data", written so the live half of #23 can call the same one and the rule is defined once.

Source code in src/ssvep/classification/decoding.py

def is_confident(scores, method: str, margin: float | None = None) -> bool:
    """Whether ``scores`` clears :func:`decision_margin` at ``margin`` (``None``: the method's
    calibrated :func:`default_confidence_margin`): the one predicate the
    offline replay (:func:`early_stop_decode`) calls to decide "stop now" versus "wait for more
    data", written so the live half of #23 can call the same one and the rule is defined once."""
    if margin is None:
        margin = default_confidence_margin(method)
    return decision_margin(scores, method) >= float(margin)

check_decision_schedule

check_decision_schedule(min_window_s, max_window_s, step_s, confidence_margin)

Refuse an early-stopping schedule that cannot mean anything, naming why: the windows and the step must be finite and positive (a zero step never reaches the maximum), the minimum no longer than the maximum, and the margin strictly between 0 and 1 (0 stops every trial at the first look; 1 or more is never met, the fixed rule under another name). max_window_s may be None, meaning the whole trial.

Source code in src/ssvep/classification/decoding.py

def check_decision_schedule(min_window_s, max_window_s, step_s, confidence_margin) -> None:
    """Refuse an early-stopping schedule that cannot mean anything, naming why: the windows and the
    step must be finite and positive (a zero step never reaches the maximum), the minimum no longer
    than the maximum, and the margin strictly between 0 and 1 (0 stops every trial at the first
    look; 1 or more is never met, the fixed rule under another name). ``max_window_s`` may be
    ``None``, meaning the whole trial."""
    for name, v in (("min_window_s", min_window_s), ("step_s", step_s)):
        if not _finite_number(v) or v <= 0:
            raise ValueError(f"early stopping needs a positive, finite {name}; got {v!r}")
    if max_window_s is not None:
        if not _finite_number(max_window_s) or max_window_s <= 0:
            raise ValueError(f"early stopping needs a positive, finite max_window_s; "
                             f"got {max_window_s!r}")
        if min_window_s > max_window_s:
            raise ValueError(f"early stopping's min_window_s ({min_window_s:g} s) is longer than "
                             f"its max_window_s ({max_window_s:g} s)")
    if not _finite_number(confidence_margin) or not 0 < confidence_margin < 1:
        raise ValueError(f"early stopping's confidence_margin must be strictly between 0 and 1; "
                         f"got {confidence_margin!r}")

early_stop_decode

early_stop_decode(trial, freqs, sfreq, method='fbcca', n_harmonics=DEFAULT_N_HARMONICS, band=None, line_hz=None, min_window_s=1.0, max_window_s=None, step_s=0.5, confidence_margin=None, replay=False, **fb_kw)

Decode trial ((n_ch, n_times)) at growing windows, stopping at the first one whose scores are confident (:func:is_confident) rather than always using max_window_s: the dynamic decision window #23 asks for. Reuses :func:classify's own scoring (:func:_scores/:func:pick); only this stopping schedule is new.

Windows run min_window_s, min_window_s + step_s, ... up to max_window_s (the whole trial when not given), always including max_window_s itself as the last look even if it falls between steps. Each window is the trial's first w seconds scored, and for FBCCA band-pass filtered, on its own, as a live decoder holding only that much of the trial would have to: the filter's edge effects at a short window are the ones the live half will meet, not borrowed from data after the window ends.

A trial shorter than max_window_s is not decoded at all (predicted_index None), the rule :func:ssvep.analysis.pipeline.decode_at_window applies at a fixed window. So every trial that contributes was eligible for every look, the dynamic and fixed numbers at the same maximum describe the same trials, and a short trial cannot contribute only when it happened to stop early (which would keep the easy trials and drop the hard ones).

Returns a dict: predicted_index (or None), decided_window_s (the window the decision was made at), decided_early (True when that is shorter than max_window_s) and history (every window scored, with its scores and margin; diagnostic, not part of the decision). A window where every candidate is undecodable (:class:UndecodableWindow) is skipped, not treated as confident or as a decision. confidence_margin None is the method's calibrated default (:data:DEFAULT_CONFIDENCE_MARGINS). Raises :class:ValueError for a schedule :func:check_decision_schedule refuses.

replay=True is the replay of a trial the live decoder stopped (#23): its flicker ended at the stop, so it is shorter than max_window_s by design, and the short-trial rule above would drop exactly the trials the rule decided. With it (and an explicit max_window_s), a short trial is fed the looks its data covers, and decides only if one of them stops it early; one that runs out of data first stays undecoded (predicted_index None, with its history). The caller says a trial was stopped from its stim_off marker (stopped: "early"), never from its length. A trial at least max_window_s long decodes exactly as without it.

Source code in src/ssvep/classification/decoding.py

def early_stop_decode(trial, freqs, sfreq, method="fbcca", n_harmonics=DEFAULT_N_HARMONICS,
                      band=None, line_hz=None, min_window_s=1.0, max_window_s=None,
                      step_s=0.5, confidence_margin=None, replay=False, **fb_kw) -> dict:
    """Decode ``trial`` (``(n_ch, n_times)``) at growing windows, stopping at the first one whose
    scores are confident (:func:`is_confident`) rather than always using ``max_window_s``: the
    dynamic decision window #23 asks for. Reuses :func:`classify`'s own scoring
    (:func:`_scores`/:func:`pick`); only this stopping schedule is new.

    Windows run ``min_window_s, min_window_s + step_s, ...`` up to ``max_window_s`` (the whole
    trial when not given), always including ``max_window_s`` itself as the last look even if it
    falls between steps. Each window is the trial's first ``w`` seconds **scored, and for FBCCA
    band-pass filtered, on its own**, as a live decoder holding only that much of the trial would
    have to: the filter's edge effects at a short window are the ones the live half will meet, not
    borrowed from data after the window ends.

    **A trial shorter than ``max_window_s`` is not decoded at all** (``predicted_index`` ``None``),
    the rule :func:`ssvep.analysis.pipeline.decode_at_window` applies at a fixed window. So every
    trial that contributes was eligible for every look, the dynamic and fixed numbers at the same
    maximum describe the same trials, and a short trial cannot contribute only when it happened to
    stop early (which would keep the easy trials and drop the hard ones).

    Returns a dict: ``predicted_index`` (or ``None``), ``decided_window_s`` (the window the decision
    was made at), ``decided_early`` (``True`` when that is shorter than ``max_window_s``) and
    ``history`` (every window scored, with its scores and margin; diagnostic, not part of the
    decision). A window where every candidate is undecodable (:class:`UndecodableWindow`) is
    skipped, not treated as confident or as a decision. ``confidence_margin`` ``None`` is the
    method's calibrated default (:data:`DEFAULT_CONFIDENCE_MARGINS`). Raises :class:`ValueError`
    for a schedule :func:`check_decision_schedule` refuses.

    ``replay=True`` is the **replay** of a trial the live decoder stopped (#23): its flicker ended
    at the stop, so it is shorter than ``max_window_s`` by design, and the short-trial rule above
    would drop exactly the trials the rule decided. With it (and an explicit ``max_window_s``), a
    short trial is fed the looks its data covers, and decides **only if one of them stops it
    early**; one that runs out of data first stays undecoded (``predicted_index`` ``None``, with its
    ``history``). The caller says a trial was stopped from its ``stim_off`` marker (``stopped:
    "early"``), never from its length. A trial at least ``max_window_s`` long decodes exactly as
    without it.
    """
    if confidence_margin is None:
        confidence_margin = default_confidence_margin(method)
    check_decision_schedule(min_window_s, max_window_s, step_s, confidence_margin)
    n_times = trial.shape[-1]
    max_w = float(max_window_s) if max_window_s is not None else n_times / float(sfreq)
    short = int(round(max_w * sfreq)) > n_times
    if (short and not (replay and max_window_s is not None)) or min_window_s > max_w:
        return {"predicted_index": None, "decided_window_s": None, "decided_early": False,
                "history": []}
    rule = EarlyStopRule(freqs, sfreq, method=method, n_harmonics=n_harmonics, band=band,
                         line_hz=line_hz, min_window_s=min_window_s, max_window_s=max_w,
                         step_s=step_s, confidence_margin=confidence_margin, **fb_kw)
    stepper = rule.stepper()
    stepper.feed(trial)
    if short and stepper.early is None:
        # Replay of a short trial that ran out of data before any look stopped it: no decision.
        return {"predicted_index": None, "decided_window_s": None, "decided_early": False,
                "history": list(stepper.history)}
    return stepper.decision()

early_stop_span_s

early_stop_span_s(max_window_s, onset_offset_s)

How long a trial that may stop early must flicker when it runs to its maximum (s), so that its maximum look is always decidable, live and in the replay (#23).

Both harnesses score the max look on round(onset_offset_s * sfreq) lead-in samples plus round(max_window_s * sfreq) window samples (:attr:EarlyStopRule.n_max), taken from the samples stamped inside [stim_on, stim_off]. That is at most (onset_offset_s + max_window_s) * sfreq + 1 samples, and a presentation of D seconds always holds at least floor(D * sfreq) of them, so D >= onset_offset_s + max_window_s + 1 / sfreq is enough on an exact clock. :data:EARLY_STOP_SPAN_GUARD_S (25 ms) covers that one sample at any rate of 100 Hz or more (10 ms) plus 15 ms of marker-to-EEG stamp jitter. Presenting exactly onset_offset_s + max_window_s rounded to whole frames is not enough: at 60 Hz 4.14 s rounds to 248 frames = 4.133 s, and every trial that ran to its maximum went undecided.

The schedule rounds this UP to whole frames (:func:ssvep.runtime.presentation.stim_frames) and the live decoder sizes its buffer from it; the max look itself stays exactly max_window_s (design.timing.stim_s), and what was actually presented is recorded on each stim_off (presented_s).

Source code in src/ssvep/classification/decoding.py

def early_stop_span_s(max_window_s, onset_offset_s) -> float:
    """How long a trial that may stop early must flicker when it runs to its maximum (s), so that
    its maximum look is always decidable, live and in the replay (#23).

    Both harnesses score the max look on ``round(onset_offset_s * sfreq)`` lead-in samples plus
    ``round(max_window_s * sfreq)`` window samples (:attr:`EarlyStopRule.n_max`), taken from the
    samples stamped inside ``[stim_on, stim_off]``. That is at most ``(onset_offset_s +
    max_window_s) * sfreq + 1`` samples, and a presentation of ``D`` seconds always holds at least
    ``floor(D * sfreq)`` of them, so ``D >= onset_offset_s + max_window_s + 1 / sfreq`` is enough on
    an exact clock. :data:`EARLY_STOP_SPAN_GUARD_S` (25 ms) covers that one sample at any rate of
    100 Hz or more (10 ms) plus 15 ms of marker-to-EEG stamp jitter. Presenting exactly
    ``onset_offset_s + max_window_s`` rounded to whole frames is not enough: at 60 Hz 4.14 s rounds
    to 248 frames = 4.133 s, and every trial that ran to its maximum went undecided.

    The schedule rounds this UP to whole frames (:func:`ssvep.runtime.presentation.stim_frames`)
    and the live decoder sizes its buffer from it; the max look itself stays exactly
    ``max_window_s`` (``design.timing.stim_s``), and what was actually presented is recorded on
    each ``stim_off`` (``presented_s``).
    """
    return float(onset_offset_s) + float(max_window_s) + EARLY_STOP_SPAN_GUARD_S

decision_rule_summary

decision_rule_summary(kind='fixed', min_window_s=None, max_window_s=None, step_s=None, confidence_margin=None)

The decision-window rule a decode ran with, as a JSON-able record (decision_rule, #23).

"fixed" is ordinary single-operating-window decoding (every recording analysed before this existed, and any run that opts out) and carries nothing but its own name. "early_stop" records the growing-window schedule :func:early_stop_decode used and the stopping statistic by name (:data:DECISION_STATISTIC), so a later change of statistic or of its threshold reads as a different rule. The setting, not which window any particular trial stopped at: that is a result, the distinction :func:filter_bank_summary and :func:mains_guard_summary draw between a decode's configuration and its outcome. Any other kind raises: there is no third rule yet (:data:DECISION_RULE_KINDS).

Source code in src/ssvep/classification/decoding.py

def decision_rule_summary(kind="fixed", min_window_s=None, max_window_s=None, step_s=None,
                          confidence_margin=None) -> dict:
    """The decision-window rule a decode ran with, as a JSON-able record (``decision_rule``, #23).

    ``"fixed"`` is ordinary single-operating-window decoding (every recording analysed before this
    existed, and any run that opts out) and carries nothing but its own name. ``"early_stop"``
    records the growing-window schedule :func:`early_stop_decode` used and the stopping statistic
    by name (:data:`DECISION_STATISTIC`), so a later change of statistic or of its threshold reads
    as a different rule. The setting, not which window any particular trial stopped at: that is a
    *result*, the distinction :func:`filter_bank_summary` and :func:`mains_guard_summary` draw
    between a decode's configuration and its outcome. Any other ``kind`` raises: there is no third
    rule yet (:data:`DECISION_RULE_KINDS`).
    """
    if kind not in DECISION_RULE_KINDS:
        raise ValueError(f"unknown decision rule {kind!r}; known: {', '.join(DECISION_RULE_KINDS)}")
    if kind == "fixed":
        return {"kind": "fixed", "statistic": None, "min_window_s": None, "max_window_s": None,
                "step_s": None, "confidence_margin": None}
    check_decision_schedule(min_window_s, max_window_s, step_s, confidence_margin)
    return {"kind": "early_stop", "statistic": DECISION_STATISTIC,
            "min_window_s": float(min_window_s),
            "max_window_s": float(max_window_s) if max_window_s is not None else None,
            "step_s": float(step_s), "confidence_margin": float(confidence_margin)}

Spectral analysis

ssvep.classification.spectral

Spectral analysis and SSVEP SNR (array-based, scipy only — no MNE).

SNR follows the standard SSVEP definition (Norcia et al. 2015): power at the stimulus frequency divided by the mean power of neighbouring bins, excluding a small guard band, reported in dB. A harmonic-summed variant aggregates f,2f,3f.

Data convention: last axis is time. (n_times,), (n_ch, n_times) or (n_trials, n_ch, n_times) all work; PSD/SNR keep the leading axes. The one exception is :func:snr_spectrum_db, which always averages over a channel axis — see its own docstring.

psd_nperseg

psd_nperseg(sfreq, n_times, nperseg_s=PSD_NPERSEG_S)

The Welch segment length :func:compute_psd uses for an n_times-sample input — the window, capped at nperseg_s seconds. With sfreq it fixes the bin width, sfreq / nperseg, exactly (:func:snr_kernel_bins).

Source code in src/ssvep/classification/spectral.py

def psd_nperseg(sfreq, n_times, nperseg_s=PSD_NPERSEG_S) -> int:
    """The Welch segment length :func:`compute_psd` uses for an ``n_times``-sample input — the
    window, capped at ``nperseg_s`` seconds. With ``sfreq`` it fixes the bin width,
    ``sfreq / nperseg``, exactly (:func:`snr_kernel_bins`)."""
    return min(int(nperseg_s * sfreq), int(n_times))

compute_psd

compute_psd(data, sfreq, fmin=1.0, fmax=None, nperseg_s=PSD_NPERSEG_S)

Welch PSD. Returns (freqs, psd) with frequency on the last axis.

Source code in src/ssvep/classification/spectral.py

def compute_psd(data, sfreq, fmin=1.0, fmax=None, nperseg_s=PSD_NPERSEG_S):
    """Welch PSD. Returns (freqs, psd) with frequency on the last axis."""
    nperseg = psd_nperseg(sfreq, data.shape[-1], nperseg_s)
    freqs, psd = welch(data, fs=sfreq, nperseg=nperseg,
                       noverlap=nperseg // 2, axis=-1)
    keep = _psd_keep(freqs, sfreq, fmin, fmax)
    return freqs[keep], psd[..., keep]

snr_kernel_bins

snr_kernel_bins(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

(n_neigh, n_guard): the SNR kernel's half-widths in bins for an n_times-sample window — the one place the kernel size is chosen (#392). :func:snr_spectrum (through every caller: the target-SNR DV, the live readout, the plotted and the live SNR spectrum), the mains guard's reach (:func:snr_reach_hz) and the edge rule (:func:snr_edge_status) all take it from here.

n_neigh = max(round(neighbor_hz / bin_w), 1) and n_guard = max(round(guard_hz / bin_w), 0), with the bin width bin_w = sfreq / nperseg (:func:psd_nperseg) exact: the quotients are computed as rationals (:func:_decimal), so the answer depends on sfreq, the window and the kernel's own parameters, and on nothing else. It used to be read off the spacing of the first two bins of whatever frequency axis a caller passed, and two slicings of one bin grid can differ in that spacing by an ulp — 0.6666666666666667 on a sub-band, 0.6666666666666666 on the full grid at 250 Hz / 1.5 s — which at a rounding tie rounded to different kernels: the same data got a different SNR depending only on how the axis was cut.

Tie rule: round half to even, Python's own round on the exact quotient. A tie is a window where a quotient is exactly a half — 1.5, 2.5 and 3.5 s at any integer rate the window is a whole number of samples at (neighbor_hz 1.0 Hz: 1.5, 2.5, 3.5 bins; guard_hz 0.2 Hz: 0.5 bins at 2.5 s); both neighbours of a tie miss the nominal width by the same half bin, so neither is more correct, and the rule is chosen to change nothing that is not broken. Half to even is what round was already doing to the target-SNR DV's quotient (computed on the full 0 Hz-to-Nyquist grid since #387), the live readout's, and the mains guard's reach (sfreq / nperseg) at every window 0.25–4 s at 125, 250, 500 and 1000 Hz — so none of them changes value here. Round half up would have moved the DV at every 2.5 s window (2 neighbour and 0 guard bins to 3 and 1) for no gain in correctness. Over every integer rate 100–1024 Hz (windows from 8 samples to 4 s) half to even differs from the old full-grid float rounding at 73 of 1 852 tie windows, all at unusual rates (196 Hz at 1.5 s, say); half up would at 519, 250, 500 and 1000 Hz among them. What changes at the four common rates is the plotted and live SNR spectrum at a rounding-tie window (1.5, 2.5 or 3.5 s, or another rate's own ties), depending on where the band's cut fell, since their sub-band axes rounded the other way.

Source code in src/ssvep/classification/spectral.py

def snr_kernel_bins(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ,
                    nperseg_s=PSD_NPERSEG_S) -> tuple:
    """``(n_neigh, n_guard)``: the SNR kernel's half-widths in bins for an ``n_times``-sample
    window — the **one** place the kernel size is chosen (#392). :func:`snr_spectrum` (through every
    caller: the target-SNR DV, the live readout, the plotted and the live SNR spectrum), the mains
    guard's reach (:func:`snr_reach_hz`) and the edge rule (:func:`snr_edge_status`) all take it
    from here.

    ``n_neigh = max(round(neighbor_hz / bin_w), 1)`` and ``n_guard = max(round(guard_hz / bin_w),
    0)``, with the bin width ``bin_w = sfreq / nperseg`` (:func:`psd_nperseg`) **exact**: the
    quotients are computed as rationals (:func:`_decimal`), so the answer depends on ``sfreq``,
    the window and the kernel's own parameters, and on nothing else. It used to be read off the
    spacing of the first two bins of whatever frequency axis a caller passed, and two slicings of
    one bin grid can differ in that spacing by an ulp — ``0.6666666666666667`` on a sub-band,
    ``0.6666666666666666`` on the full grid at 250 Hz / 1.5 s — which at a rounding tie rounded to
    different kernels: the same data got a different SNR depending only on how the axis was cut.

    **Tie rule: round half to even**, Python's own ``round`` on the exact quotient. A tie is a
    window where a quotient is exactly a half — 1.5, 2.5 and 3.5 s at any integer rate the window is
    a whole number of samples at (``neighbor_hz`` 1.0 Hz: 1.5, 2.5, 3.5 bins; ``guard_hz`` 0.2 Hz:
    0.5 bins at 2.5 s); both neighbours of a tie miss the nominal width by the same half bin, so
    neither is more correct, and the rule is chosen to change nothing that is not broken. Half to
    even is what ``round`` was already doing to the target-SNR DV's quotient (computed on the full
    0 Hz-to-Nyquist grid since #387), the live readout's, and the mains guard's reach
    (``sfreq / nperseg``) at **every** window 0.25–4 s at 125, 250, 500 and 1000 Hz — so none of
    them changes value here. Round half *up* would have moved the DV at every 2.5 s window (2
    neighbour and 0 guard bins to 3 and 1) for no gain in correctness. Over every integer rate
    100–1024 Hz (windows from 8 samples to 4 s) half to even differs from the old full-grid float
    rounding at 73 of 1 852 tie windows, all at unusual rates (196 Hz at 1.5 s, say); half up would
    at 519, 250, 500 and 1000 Hz among them. What changes at the four common rates is the plotted
    and live SNR spectrum at a rounding-tie window (1.5, 2.5 or 3.5 s, or another rate's own
    ties), depending on where the band's cut fell, since their sub-band axes rounded the other way."""
    nperseg = psd_nperseg(sfreq, n_times, nperseg_s)
    if nperseg < 1:
        raise ValueError(f"snr_kernel_bins: no Welch segment for a {n_times}-sample window")
    per_bin = Fraction(nperseg) / _decimal(sfreq)            # 1 / bin_w, exactly
    return (max(round(_decimal(neighbor_hz) * per_bin), 1),
            max(round(_decimal(guard_hz) * per_bin), 0))

snr_spectrum

snr_spectrum(psd, freqs, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, kernel=None)

SNR at every bin: power / mean(neighbouring bins, guard excluded). Same shape.

kernel is (n_neigh, n_guard) from :func:snr_kernel_bins — every SNR this toolbox computes passes it, so the kernel is fixed by sfreq and the window, never by how freqs was sliced (#392). Without it (a caller that has only an axis), the kernel is sized the old way, from freqs[1] - freqs[0], which at a rounding-tie window can depend on the slicing.

The kernel needs at least one real neighbour on each side. At a very low sfreq (or a very narrow band), the PSD's frequency axis can be shorter than that: compute_psd's bin width is sfreq / nperseg, so a low rate coarsens the bins and leaves fewer of them in a given band at the same time (found during #325, on a synthetic decimation-testing config). Unshrunk, the guard band + neighbours kernel would be built wider than the data, and np.convolve(..., mode="same") then returns an array as long as the kernel, not the data — psd / noise then fails with a shape-mismatch, far from where the real cause is. So the guard is dropped first, then the neighbour half-width is shrunk to whatever fits; if even one bin each side doesn't fit, this refuses with a clear message instead of failing inside a broadcast.

Source code in src/ssvep/classification/spectral.py

def snr_spectrum(psd, freqs, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, kernel=None):
    """SNR at every bin: power / mean(neighbouring bins, guard excluded). Same shape.

    ``kernel`` is ``(n_neigh, n_guard)`` from :func:`snr_kernel_bins` — every SNR this toolbox
    computes passes it, so the kernel is fixed by ``sfreq`` and the window, never by how ``freqs``
    was sliced (#392). Without it (a caller that has only an axis), the kernel is sized the old way,
    from ``freqs[1] - freqs[0]``, which at a rounding-tie window can depend on the slicing.

    The kernel needs at least one real neighbour on each side. At a very low ``sfreq`` (or a very
    narrow band), the PSD's frequency axis can be shorter than that: ``compute_psd``'s bin width is
    ``sfreq / nperseg``, so a low rate coarsens the bins *and* leaves fewer of them in a given band
    at the same time (found during #325, on a synthetic decimation-testing config). Unshrunk, the
    guard band + neighbours kernel would be built wider than the data, and
    ``np.convolve(..., mode="same")`` then returns an array as long as the *kernel*, not the data —
    ``psd / noise`` then fails with a shape-mismatch, far from where the real cause is. So the guard
    is dropped first, then the neighbour half-width is shrunk to whatever fits; if even one bin each
    side doesn't fit, this refuses with a clear message instead of failing inside a broadcast.
    """
    n_bins = len(freqs)
    if n_bins < 2:
        raise ValueError(
            f"snr_spectrum: only {n_bins} frequency bin(s) in the PSD — SNR needs at least one "
            "neighbouring bin to compare against. This happens at a very low sfreq or over a very "
            "narrow band, where compute_psd's own bins are too few to resolve more. Raise sfreq or "
            "widen the analysed band.")
    if kernel is None:
        n_neigh, n_guard = _float_kernel_bins(float(freqs[1] - freqs[0]), neighbor_hz, guard_hz)
    else:
        n_neigh, n_guard = (int(v) for v in kernel)
    kernel_len = 2 * n_neigh + 2 * n_guard + 1
    if kernel_len > n_bins:
        max_neigh = (n_bins - 1) // 2
        if max_neigh < 1:
            raise ValueError(
                f"snr_spectrum: the smoothing kernel ({kernel_len} bins: {n_neigh} neighbour + "
                f"{n_guard} guard each side of the target bin) is wider than the {n_bins}-bin PSD "
                "— there is no room for even one real neighbouring bin. This happens at a very low "
                "sfreq or over a very narrow band, where compute_psd's own bins are too coarse and "
                "too few at once. Raise sfreq, widen the analysed band, or shorten "
                "neighbor_hz/guard_hz.")
        n_guard, n_neigh = 0, max_neigh   # drop the guard band first, then shrink to what fits
    kernel = np.concatenate([np.ones(n_neigh), np.zeros(n_guard * 2 + 1), np.ones(n_neigh)])
    kernel /= kernel.sum()
    noise = np.apply_along_axis(lambda m: np.convolve(m, kernel, mode="same"), -1, psd)
    with np.errstate(divide="ignore", invalid="ignore"):
        return psd / noise

snr_at

snr_at(freqs, snr_spec, f)

Index an SNR spectrum at the bin nearest frequency f.

Source code in src/ssvep/classification/spectral.py

def snr_at(freqs, snr_spec, f):
    """Index an SNR spectrum at the bin nearest frequency ``f``."""
    return snr_spec[..., int(np.argmin(np.abs(freqs - f)))]

target_snr_db

target_snr_db(psd, freqs, f, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, n_harmonics=1, harmonics=None, kernel=None)

SNR (dB) at target f, optionally harmonic-summed over f,2f,3f...

kernel is passed to :func:snr_spectrum (:func:snr_kernel_bins, #392).

harmonics is an explicit tuple/list of 1-based harmonic orders to sum (e.g. (1,) when a mains guard has dropped the 2nd harmonic); when given, it overrides n_harmonics (which otherwise implies range(1, n_harmonics + 1)) — see :func:ssvep.analysis.pipeline.mean_target_snr_db (#309). An empty harmonics (every order guarded away — the target's own fundamental sits on the mains line) returns NaN: there is nothing left to sum, not zero.

Source code in src/ssvep/classification/spectral.py

def target_snr_db(psd, freqs, f, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, n_harmonics=1,
                  harmonics=None, kernel=None):
    """SNR (dB) at target ``f``, optionally harmonic-summed over f,2f,3f...

    ``kernel`` is passed to :func:`snr_spectrum` (:func:`snr_kernel_bins`, #392).

    ``harmonics`` is an explicit tuple/list of 1-based harmonic orders to sum (e.g. ``(1,)`` when a
    mains guard has dropped the 2nd harmonic); when given, it overrides ``n_harmonics`` (which
    otherwise implies ``range(1, n_harmonics + 1)``) — see
    :func:`ssvep.analysis.pipeline.mean_target_snr_db` (#309). An empty ``harmonics`` (every order
    guarded away — the target's own fundamental sits on the mains line) returns NaN: there is
    nothing left to sum, not zero.
    """
    orders = list(harmonics) if harmonics is not None else list(range(1, n_harmonics + 1))
    hz = [f * k for k in orders if f * k <= freqs[-1]]
    if not hz:
        return np.full(psd.shape[:-1], np.nan, dtype=float)
    spec = snr_spectrum(psd, freqs, neighbor_hz, guard_hz, kernel=kernel)
    vals = np.stack([snr_at(freqs, spec, h) for h in hz], axis=0)
    return 10.0 * np.log10(vals.mean(axis=0))

target_snr_from_data

target_snr_from_data(data, sfreq, f, n_harmonics=3, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, average_trials=True, harmonics=None)

Convenience: per-channel SNR (dB) at f straight from a data array.

data (n_trials, n_ch, n_times) or (n_ch, n_times). With trials and average_trials, the PSD is averaged across trials first. harmonics, when given, is an explicit set of orders, overriding n_harmonics (which otherwise implies range(1, n_harmonics + 1)).

The summed value of :func:target_snr_orders_from_data — one computation, so the live readout (ssvep.ui.run) and the offline DV cannot disagree about which bins a target's SNR may use. An order without a full real noise neighbourhood (too close to 0 Hz or Nyquist for the window,

387) is left out of the sum; NaN when no order is left.

Source code in src/ssvep/classification/spectral.py

def target_snr_from_data(data, sfreq, f, n_harmonics=3, neighbor_hz=SNR_NEIGHBOR_HZ,
                         guard_hz=SNR_GUARD_HZ, average_trials=True, harmonics=None):
    """Convenience: per-channel SNR (dB) at ``f`` straight from a data array.

    data ``(n_trials, n_ch, n_times)`` or ``(n_ch, n_times)``. With trials and
    average_trials, the PSD is averaged across trials first. ``harmonics``, when given, is an
    explicit set of orders, overriding ``n_harmonics`` (which otherwise implies
    ``range(1, n_harmonics + 1)``).

    The summed value of :func:`target_snr_orders_from_data` — one computation, so the live readout
    (``ssvep.ui.run``) and the offline DV cannot disagree about which bins a target's SNR may use.
    An order without a full real noise neighbourhood (too close to 0 Hz or Nyquist for the window,
    #387) is left out of the sum; NaN when no order is left.
    """
    orders = list(harmonics) if harmonics is not None else list(range(1, n_harmonics + 1))
    return target_snr_orders_from_data(data, sfreq, f, orders, neighbor_hz=neighbor_hz,
                                       guard_hz=guard_hz, average_trials=average_trials)[0]

psd_bin_grid

psd_bin_grid(sfreq, n_times, nperseg_s=PSD_NPERSEG_S)

Exactly the frequency axis compute_psd(x, sfreq, fmin=0.0) returns for an n_times-sample input — the axis the target-SNR DV is computed on (:func:target_snr_orders_from_data): Welch's bins (:func:_welch_freqs) under :func:compute_psd's own mask (:func:_psd_keep), which since #392 keeps every one of them up to and including the Nyquist bin. Lets the edge rule be evaluated without any data (a record's D7 key derives it from sfreq and a window alone).

Source code in src/ssvep/classification/spectral.py

def psd_bin_grid(sfreq, n_times, nperseg_s=PSD_NPERSEG_S):
    """Exactly the frequency axis ``compute_psd(x, sfreq, fmin=0.0)`` returns for an
    ``n_times``-sample input — the axis the target-SNR DV is computed on
    (:func:`target_snr_orders_from_data`): Welch's bins (:func:`_welch_freqs`) under
    :func:`compute_psd`'s own mask (:func:`_psd_keep`), which since #392 keeps every one of them up
    to and including the Nyquist bin. Lets the edge rule be evaluated without any data (a record's
    D7 key derives it from ``sfreq`` and a window alone)."""
    raw = _welch_freqs(sfreq, n_times, nperseg_s)
    return raw[_psd_keep(raw, sfreq, 0.0, None)]

snr_edge_status

snr_edge_status(f, harmonics, sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

{order: EDGE_OK | EDGE_REMEASURED | EDGE_UNAVAILABLE} for each of target f's harmonics that has a bin at all, at an n_times-sample window: f * k at or below the top bin of the fixed axis (:func:psd_bin_grid) or of the pre-#387 one. An order past both was never summed, before or after #387, and is absent here; one past the fixed axis's top but not the old axis's (the old fetch kept a Nyquist bin one ulp above sfreq / 2 that :func:compute_psd's default mask drops) was summed before and is not now, so it is EDGE_UNAVAILABLE.

A function of f, the orders, sfreq and the window alone — no data — so the record (:func:ssvep.analysis.pipeline.target_snr_harmonics_summary) and the D7 key's reading of a record written before the fix (:func:ssvep.analysis.contrast.snr_edge_key) apply one rule.

The rule (#387). Let bin_w = sfreq / min(n_times, nperseg_s * sfreq) and margin = max(round(neighbor_hz / bin_w), 1) + round(guard_hz / bin_w) bins (the SNR kernel's half-width). An order at h = k * f, read at its nearest bin j:

• is available iff margin <= j <= J - margin, J the Nyquist bin: roughly margin * bin_w <= h <= Nyquist - margin * bin_w (1.25 Hz at a 4 s window);
• is unchanged (bit-identical to the pre-#387 value) iff the old fetch, [max(1, f - 5), f * (max(harmonics) + 1) + 5], gave it that same full kernel. That is two conditions. (a) A full kernel of real bins inside the old axis: roughly h - margin * bin_w >= max(1, f - 5) in bins, so h above about 1 + (margin + 0.5) * bin_w for a low target — everything from 2.15 Hz at a 4 s window, and from about 3 Hz at any window from 0.5 s (the bin quantisation makes the exact edge move with the window, not monotonically) — and the mirror-image condition at the top. (b) The same kernel size. :func:snr_spectrum sizes its kernel from round(neighbor_hz / bin_w) and round(guard_hz / bin_w), with bin_w read off its axis's first two bins, and at a window where either quotient is exactly a half (1.5, 2.5 or 3.5 s, with sfreq * window a whole number of samples) the old masked axis's spacing could differ from the full grid's by one ulp and round the other way — a 1-bin rather than 2-bin neighbourhood at 1.5 s, systematically about 2.5 dB apart on a real response. Whether it did depended on the target, through where the mask started. Such an order is EDGE_REMEASURED however far it is from both edges. The fixed path's kernel comes from sfreq and the window alone (:func:snr_kernel_bins, #392), so it never depends on the target.

"Now" is the current toolbox: since #392 its axis always has the Nyquist bin (:func:_psd_keep) and its kernel is exact, so this compares a pre-#387 record with what the toolbox computes today, not with what #387 alone computed.

Source code in src/ssvep/classification/spectral.py

def snr_edge_status(f, harmonics, sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ,
                    guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S) -> dict:
    """``{order: EDGE_OK | EDGE_REMEASURED | EDGE_UNAVAILABLE}`` for each of target ``f``'s
    ``harmonics`` that has a bin at all, at an ``n_times``-sample window: ``f * k`` at or below the
    top bin of the fixed axis (:func:`psd_bin_grid`) or of the pre-#387 one. An order past both
    was never summed, before or after #387, and is absent here; one past the fixed axis's top but
    not the old axis's (the old fetch kept a Nyquist bin one ulp above ``sfreq / 2`` that
    :func:`compute_psd`'s default mask drops) was summed before and is not now, so it is
    ``EDGE_UNAVAILABLE``.

    A function of ``f``, the orders, ``sfreq`` and the window alone — no data — so the record
    (:func:`ssvep.analysis.pipeline.target_snr_harmonics_summary`) and the D7 key's reading of a
    record written before the fix (:func:`ssvep.analysis.contrast.snr_edge_key`) apply one rule.

    **The rule (#387).** Let ``bin_w = sfreq / min(n_times, nperseg_s * sfreq)`` and ``margin =
    max(round(neighbor_hz / bin_w), 1) + round(guard_hz / bin_w)`` bins (the SNR kernel's
    half-width). An order at ``h = k * f``, read at its nearest bin ``j``:

    * is **available** iff ``margin <= j <= J - margin``, ``J`` the Nyquist bin: roughly
      ``margin * bin_w <= h <= Nyquist - margin * bin_w`` (1.25 Hz at a 4 s window);
    * is **unchanged** (bit-identical to the pre-#387 value) iff the old fetch,
      ``[max(1, f - 5), f * (max(harmonics) + 1) + 5]``, gave it that same full kernel. That is
      two conditions. (a) A full kernel of real bins inside the old axis: roughly
      ``h - margin * bin_w >= max(1, f - 5)`` in bins, so ``h`` above about
      ``1 + (margin + 0.5) * bin_w`` for a low target — everything from 2.15 Hz at a 4 s window,
      and from about 3 Hz at any window from 0.5 s (the bin quantisation makes the exact edge move
      with the window, not monotonically) — and the mirror-image condition at the top. (b) The
      same kernel *size*. :func:`snr_spectrum` sizes its kernel from ``round(neighbor_hz / bin_w)``
      and ``round(guard_hz / bin_w)``, with ``bin_w`` read off its axis's first two bins, and at a
      window where either quotient is exactly a half (1.5, 2.5 or 3.5 s, with ``sfreq * window`` a
      whole number of samples) the old masked axis's spacing could differ from the full grid's by
      one ulp and round the other way — a 1-bin rather than 2-bin neighbourhood at 1.5 s,
      systematically about 2.5 dB apart on a real response. Whether it did depended on the target,
      through where the mask started. Such an order is ``EDGE_REMEASURED`` however far it is from
      both edges. The fixed path's kernel comes from ``sfreq`` and the window alone
      (:func:`snr_kernel_bins`, #392), so it never depends on the target.

    "Now" is the current toolbox: since #392 its axis always has the Nyquist bin
    (:func:`_psd_keep`) and its kernel is exact, so this compares a pre-#387 record with what the
    toolbox computes today, not with what #387 alone computed.
    """
    harmonics = [int(k) for k in harmonics]
    raw = _welch_freqs(sfreq, n_times, nperseg_s)                # the pre-fix axis's source
    grid = psd_bin_grid(sfreq, n_times, nperseg_s)              # the fixed path's own axis
    if not harmonics or len(grid) == 0:
        return {}
    kernel = snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
    old = _pre_fix_axis(raw, sfreq, f, harmonics)
    old_kernel = _kernel_of(old, neighbor_hz, guard_hz)
    out = {}
    for k in harmonics:
        h = f * k
        if h > grid[-1] and not (len(old) and h <= old[-1]):
            continue            # past the top bin both before and after the fix: never summed
        now_ok = _full_kernel_at(grid, h, kernel)
        was_ok = (old_kernel == kernel and h <= old[-1]
                  and _full_kernel_at(old, h, old_kernel))
        out[k] = (EDGE_OK if now_ok and was_ok else EDGE_REMEASURED if now_ok
                  else EDGE_UNAVAILABLE)
    return out

snr_kernel_tie

snr_kernel_tie(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

Whether an n_times-sample window is a rounding tie for :func:snr_spectrum's kernel size: neighbor_hz / bin_w or guard_hz / bin_w (bin_w = sfreq / nperseg) within :data:_KERNEL_TIE_TOL of a half-integer, where two float spacings of the same bin grid one ulp apart can round to different bin counts (#387 review; :func:snr_edge_status). A superset of the windows where the pre-#387 axis actually sized a target's kernel differently from the full grid (a half that both roundings clamp to the same count, such as a 0.5-bin neighbourhood, which is at least 1 either way, is not a tie) — whether it did depends on the target — so a caller that cannot afford :func:snr_edge_status at every window (the D7 key's pre-fix derivation) may screen with it and only ever over-check. At an integer sfreq with the default kernel: windows of 1.5, 2.5 and 3.5 s whose sample count is a whole number.

Source code in src/ssvep/classification/spectral.py

def snr_kernel_tie(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ,
                   nperseg_s=PSD_NPERSEG_S) -> bool:
    """Whether an ``n_times``-sample window is a rounding tie for :func:`snr_spectrum`'s kernel
    size: ``neighbor_hz / bin_w`` or ``guard_hz / bin_w`` (``bin_w = sfreq / nperseg``) within
    :data:`_KERNEL_TIE_TOL` of a half-integer, where two float spacings of the same bin grid one ulp
    apart can round to different bin counts (#387 review; :func:`snr_edge_status`). A superset of
    the windows where the pre-#387 axis actually sized a target's kernel differently from the full
    grid (a half that both roundings clamp to the same count, such as a 0.5-bin neighbourhood, which
    is at least 1 either way, is not a tie) — whether it did depends on the target — so a caller that cannot afford
    :func:`snr_edge_status` at every window (the D7 key's pre-fix derivation) may screen with it and
    only ever over-check. At an integer ``sfreq`` with the default kernel: windows of 1.5, 2.5 and
    3.5 s whose sample count is a whole number."""
    nperseg = psd_nperseg(sfreq, n_times, nperseg_s)
    if nperseg < 1:
        return False
    for hz, least in ((neighbor_hz, 1), (guard_hz, 0)):     # snr_spectrum's own clamps
        q = float(hz) * nperseg / float(sfreq)
        if abs(q % 1.0 - 0.5) < _KERNEL_TIE_TOL and max(int(q), least) != max(int(q) + 1, least):
            return True
    return False

snr_exact_bins_changed

snr_exact_bins_changed(f, harmonics, sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

The orders of target f's harmonics whose target SNR at an n_times-sample window differs between a record written after #387 but before #392 and the current toolbox — what the D7 key derives for such a record (:func:ssvep.analysis.contrast.snr_edge_key).

That record's DV ran on :func:_pre_exact_grid (no Nyquist bin where rfftfreq put it one ulp high) with the kernel :func:_kernel_of read off that grid; today's runs on :func:psd_bin_grid with :func:snr_kernel_bins. An order differs when one path sums it and the other does not (in practice: an order at the top of the band that the extra Nyquist bin gives a full kernel), or both sum it with differently sized kernels (never, at 125, 250, 500 or 1000 Hz, where round-half-to-even matches the full grid's float rounding at every window; it does at 73 tie windows at other integer rates 100-1024, e.g. 196 Hz / 1.5 s). Otherwise its value is bit-identical: same bins, same kernel, and the extra Nyquist bin lies outside its kernel. Which orders the mains guard kept is the record's own answer, taken as given here; the D7 key separately counts every order of a guarded target at a window where the guard's reach itself moved (:func:reach_kernel_moved).

Source code in src/ssvep/classification/spectral.py

def snr_exact_bins_changed(f, harmonics, sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ,
                           guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S) -> tuple:
    """The orders of target ``f``'s ``harmonics`` whose target SNR at an ``n_times``-sample window
    differs between a record written after #387 but before #392 and the current toolbox — what the
    D7 key derives for such a record (:func:`ssvep.analysis.contrast.snr_edge_key`).

    That record's DV ran on :func:`_pre_exact_grid` (no Nyquist bin where ``rfftfreq`` put it one
    ulp high) with the kernel :func:`_kernel_of` read off that grid; today's runs on
    :func:`psd_bin_grid` with :func:`snr_kernel_bins`. An order differs when one path sums it and
    the other does not (in practice: an order at the top of the band that the extra Nyquist bin
    gives a full kernel), or both sum it with differently sized kernels (never, at 125, 250, 500 or
    1000 Hz, where round-half-to-even matches the full grid's float rounding at every window; it
    does at 73 tie windows at other integer rates 100-1024, e.g. 196 Hz / 1.5 s). Otherwise its
    value is bit-identical: same bins, same kernel, and the extra Nyquist bin lies outside its
    kernel. Which orders the mains guard kept is the record's own answer, taken as given here; the
    D7 key separately counts every order of a guarded target at a window where the guard's reach
    itself moved (:func:`reach_kernel_moved`)."""
    harmonics = [int(k) for k in harmonics]
    then_grid = _pre_exact_grid(sfreq, n_times, nperseg_s)
    now_grid = psd_bin_grid(sfreq, n_times, nperseg_s)
    if not harmonics or len(now_grid) < 2:
        return ()
    then_k = _kernel_of(then_grid, neighbor_hz, guard_hz)
    now_k = snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
    out = []
    for k in harmonics:
        h = f * k
        then = bool(len(then_grid)) and h <= then_grid[-1] and _full_kernel_at(then_grid, h, then_k)
        now = h <= now_grid[-1] and _full_kernel_at(now_grid, h, now_k)
        if then != now or (then and then_k != now_k):
            out.append(k)
    return tuple(out)

exact_bins_window_differs

exact_bins_window_differs(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

Whether :func:snr_exact_bins_changed can name anything at this window: the pre-#392 grid lacks the Nyquist bin, or either pre-#392 kernel (the grid's, or the mains guard reach's) differs from :func:snr_kernel_bins. A cheap screen for the D7 key's candidate windows.

Source code in src/ssvep/classification/spectral.py

def exact_bins_window_differs(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ,
                              nperseg_s=PSD_NPERSEG_S) -> bool:
    """Whether :func:`snr_exact_bins_changed` can name anything at this window: the pre-#392 grid
    lacks the Nyquist bin, or either pre-#392 kernel (the grid's, or the mains guard reach's)
    differs from :func:`snr_kernel_bins`. A cheap screen for the D7 key's candidate windows."""
    then_grid = _pre_exact_grid(sfreq, n_times, nperseg_s)
    now_grid = psd_bin_grid(sfreq, n_times, nperseg_s)
    if len(now_grid) < 2:
        return False
    if len(then_grid) != len(now_grid):
        return True
    now_k = snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
    return (_kernel_of(then_grid, neighbor_hz, guard_hz) != now_k
            or reach_kernel_moved(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s))

reach_kernel_moved

reach_kernel_moved(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

Whether the mains guard's reach (:func:snr_reach_hz) at this window differs between the float arithmetic that wrote every record before #392 (:func:float_snr_guard_tolerance_hz) and the exact kernel it uses now. Measured never — at any window from 8 samples to 4 s at every integer rate 100–1024 Hz and every quarter-integer rate 100–525 Hz, round of the float quotient landed where round half to even on the exact one does — so this is insurance: where it would fire, an older record's guard may have kept different orders, so the D7 key counts every order of a guarded target there.

Source code in src/ssvep/classification/spectral.py

def reach_kernel_moved(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ,
                       nperseg_s=PSD_NPERSEG_S) -> bool:
    """Whether the mains guard's reach (:func:`snr_reach_hz`) at this window differs between the
    float arithmetic that wrote every record before #392 (:func:`float_snr_guard_tolerance_hz`) and
    the exact kernel it uses now. Measured never — at any window from 8 samples to 4 s at every
    integer rate 100–1024 Hz and every quarter-integer rate 100–525 Hz, ``round`` of the float
    quotient landed where round half to even on the exact one does — so this is insurance: where it
    would fire, an older record's guard may have kept different orders, so the D7 key counts every
    order of a guarded target there."""
    nperseg = psd_nperseg(sfreq, n_times, nperseg_s)
    if nperseg < 1:
        return False
    return (_float_kernel_bins(float(sfreq) / float(nperseg), neighbor_hz, guard_hz)
            != snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s))

snr_edge_report

snr_edge_report(f, harmonics, sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

The record entries for target f's orders whose SNR the #387 edge fix touched — every order :func:snr_edge_status does not call EDGE_OK — each with its order, harmonic, window, status and a checkable reason; [] for a target whose every order is far enough from both spectrum edges — every target at or above about 2.25 Hz at a 4 s window and about 3 Hz at any window from 0.5 s, so every shipped design — except at a rounding-tie window (1.5, 2.5 or 3.5 s; :func:snr_edge_status), where an order far from both edges can be listed too because the old path sized its kernel differently. The same shape for the real pipeline and the synthetic record.

Source code in src/ssvep/classification/spectral.py

def snr_edge_report(f, harmonics, sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ,
                    guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S) -> list:
    """The record entries for target ``f``'s orders whose SNR the #387 edge fix touched — every
    order :func:`snr_edge_status` does not call ``EDGE_OK`` — each with its order, harmonic, window,
    status and a checkable reason; ``[]`` for a target whose every order is far enough from both
    spectrum edges — every target at or above about 2.25 Hz at a 4 s window and about 3 Hz at any
    window from 0.5 s, so every shipped design — except at a rounding-tie window (1.5, 2.5 or
    3.5 s; :func:`snr_edge_status`), where an order far from both edges can be listed too because
    the old path sized its kernel differently. The same shape for the real pipeline and the
    synthetic record."""
    harmonics = [int(k) for k in harmonics]
    status = snr_edge_status(f, harmonics, sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
    raw = _welch_freqs(sfreq, n_times, nperseg_s)                # the pre-fix axis's source
    grid = psd_bin_grid(sfreq, n_times, nperseg_s)              # the fixed path's own axis
    out = []
    for k, st in status.items():
        if st == EDGE_OK:
            continue
        n_neigh, n_guard = snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
        old_kernel = _kernel_of(_pre_fix_axis(raw, sfreq, f, harmonics), neighbor_hz, guard_hz)
        reach = (n_neigh + n_guard) * float(grid[1] - grid[0])
        h = float(f * k)
        edge = "0 Hz" if h < float(grid[-1]) / 2.0 else "Nyquist"
        if st == EDGE_UNAVAILABLE:
            why = (f"{h:g} Hz is within the SNR kernel's {reach:g} Hz half-width of {edge} at this "
                   f"{int(n_times)}-sample window, so it has no full set of real neighbouring bins "
                   "to compare against; left out of the target SNR rather than reported inflated")
        elif old_kernel is not None and old_kernel != (n_neigh, n_guard):
            why = (f"{h:g} Hz is measured with {n_neigh} neighbour and {n_guard} guard bins each "
                   "side; before #387 the narrower PSD fetch's bin spacing differed from the full "
                   f"spectrum's by a rounding error, which at this {int(n_times)}-sample window (a "
                   f"rounding tie) sized the kernel at {old_kernel[0]} neighbour and "
                   f"{old_kernel[1]} guard bins instead, so its value changes")
        else:
            why = (f"{h:g} Hz now has a full set of real neighbouring bins; before #387 the PSD "
                   f"fetch stopped short of its {reach:g} Hz neighbourhood (on the {edge} side) and "
                   "the missing bins were zero-padded, which read high")
        out.append({"response_hz": float(f), "order": int(k), "harmonic_hz": h,
                    "n_samples": int(n_times), "status": st, "reason": why})
    return out

target_snr_orders_from_data

target_snr_orders_from_data(data, sfreq, f, harmonics, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, average_trials=True)

Target f's SNR summed over the explicit harmonics, plus each order's own SNR.

Returns (summed_db, by_order_db): summed_db is exactly what target_snr_from_data(data, sfreq, f, harmonics=harmonics, ...) returns (that function is this one's summed value), and by_order_db maps each order that was actually summed to its own per-channel SNR in dB, from that same spectrum.

Why both (#309 review): a target that lost its 2nd harmonic to the mains guard sums fewer orders than one that did not, so comparing the two sums compares different estimators. Recording each order lets a contrast re-sum both arms on the orders they have in common (:func:ssvep.analysis.contrast.snr_db_on_orders).

Real neighbours only (#387). The PSD is taken over every bin from 0 Hz to Nyquist, and an order is summed only when its bin gets :func:snr_spectrum's full kernel of real neighbouring bins (:func:_full_kernel_at). The fetch used to start at max(1.0, f - 5), so a target below about 2.25 Hz (at a 4 s window; higher at shorter ones) had part of its noise neighbourhood zero-padded and its SNR read high — +3 to +7 dB at 1.5 Hz. Fetching from 0 Hz gives such a bin its real neighbours; a bin within the kernel's half-width of 0 Hz or Nyquist has none to get, so it is left out (absent from by_order_db), never reported inflated. Every order whose neighbourhood the old fetch already covered, with the same kernel size, reads bit-for-bit what it did; at a rounding-tie window (1.5, 2.5 or 3.5 s) the old, narrower axis could size the kernel differently even far from both edges (:func:snr_edge_status is the rule, from the window alone, and calls such an order EDGE_REMEASURED).

Source code in src/ssvep/classification/spectral.py

def target_snr_orders_from_data(data, sfreq, f, harmonics, neighbor_hz=SNR_NEIGHBOR_HZ,
                                guard_hz=SNR_GUARD_HZ, average_trials=True):
    """Target ``f``'s SNR summed over the explicit ``harmonics``, plus each order's own SNR.

    Returns ``(summed_db, by_order_db)``: ``summed_db`` is exactly what
    ``target_snr_from_data(data, sfreq, f, harmonics=harmonics, ...)`` returns (that function is
    this one's summed value), and ``by_order_db`` maps each order that was actually summed to its
    own per-channel SNR in dB, from that same spectrum.

    Why both (#309 review): a target that lost its 2nd harmonic to the mains guard sums fewer orders
    than one that did not, so comparing the two sums compares different estimators. Recording each
    order lets a contrast re-sum both arms on the orders they have in common
    (:func:`ssvep.analysis.contrast.snr_db_on_orders`).

    **Real neighbours only (#387).** The PSD is taken over every bin from 0 Hz to Nyquist, and an
    order is summed only when its bin gets :func:`snr_spectrum`'s full kernel of *real* neighbouring
    bins (:func:`_full_kernel_at`). The fetch used to start at ``max(1.0, f - 5)``, so a target below
    about 2.25 Hz (at a 4 s window; higher at shorter ones) had part of its noise neighbourhood
    zero-padded and its SNR read high — +3 to +7 dB at 1.5 Hz. Fetching from 0 Hz gives such a bin
    its real neighbours; a bin within the kernel's half-width of 0 Hz or Nyquist has none to get, so
    it is left out (absent from ``by_order_db``), never reported inflated. Every order whose
    neighbourhood the old fetch already covered, with the same kernel size, reads bit-for-bit what
    it did; at a rounding-tie window (1.5, 2.5 or 3.5 s) the old, narrower axis could size the
    kernel differently even far from both edges (:func:`snr_edge_status` is the rule, from the
    window alone, and calls such an order ``EDGE_REMEASURED``).
    """
    harmonics = [int(k) for k in harmonics]
    arr = np.asarray(data)
    lead = arr.shape[1:-1] if (arr.ndim == 3 and average_trials) else arr.shape[:-1]
    if not harmonics:
        return np.full(lead, np.nan, dtype=float), {}
    freqs, psd = compute_psd(arr, sfreq, fmin=0.0)           # every bin, DC to Nyquist (#387)
    if psd.ndim == 3 and average_trials:
        psd = psd.mean(axis=0)
    kernel = (snr_kernel_bins(sfreq, arr.shape[-1], neighbor_hz, guard_hz)
              if psd_nperseg(sfreq, arr.shape[-1]) >= 1 else None)
    kept = [k for k in harmonics
            if len(freqs) and f * k <= freqs[-1]
            and _full_kernel_at(freqs, f * k, kernel)]
    if not kept:
        return np.full(psd.shape[:-1], np.nan, dtype=float), {}
    spec = snr_spectrum(psd, freqs, neighbor_hz, guard_hz, kernel=kernel)
    lin = {k: snr_at(freqs, spec, f * k) for k in kept}
    summed = 10.0 * np.log10(np.stack([lin[k] for k in kept], axis=0).mean(axis=0))
    with np.errstate(divide="ignore", invalid="ignore"):
        return summed, {k: 10.0 * np.log10(v) for k, v in lin.items()}

snr_reach_hz

snr_reach_hz(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

How far (Hz) from a target bin a narrowband line can still change that bin's SNR (#309).

:func:snr_spectrum divides a bin's power by the mean of the bins n_guard + 1 to n_guard + n_neigh either side of it, so a line anywhere in that neighbourhood inflates the noise estimate and pulls the SNR down, exactly as a line in the signal bin pushes it up. And a line is not one bin wide: through the Welch segment's Hann window it spreads over :data:WELCH_MAIN_LOBE_HALF_BINS bins either side of its own, so a line just outside the neighbourhood still leaks into its outermost bin. The reach is therefore (n_neigh + n_guard + WELCH_MAIN_LOBE_HALF_BINS) * bin_w, with n_neigh, n_guard and the bin width sfreq / nperseg computed exactly as :func:compute_psd and :func:snr_spectrum compute them for an n_times-sample input: derived from the estimator's own parameters, so changing them moves the guard with them. At 250 Hz with 4 s trials that is (4 + 1 + 2) * 0.25 = 1.75 Hz; the decoder's own guard tolerance, 1/T, is 0.25 Hz there, which is why the SNR guard is wider than the decoder's.

It covers the main lobe, not the side lobes: an off-bin line tens of dB above the noise floor still leaks a smaller, decaying error beyond this, which no finite guard removes.

Not monotonic in n_times (#362 review): n_neigh/n_guard round to whole bins, so a shorter window can occasionally come out narrower than a longer one (e.g. at 250 Hz: 1.75 Hz at 1000 samples, 1.667 Hz at 600, 2.0 Hz at 500) — a shorter window does not always widen the tolerance. Harmless for the target-SNR DV: each target's reach is evaluated at the same n_times its own SNR is summed over (:func:ssvep.analysis.pipeline._snr_guard_by_target), so the two never disagree about which window they describe.

Source code in src/ssvep/classification/spectral.py

def snr_reach_hz(sfreq, n_times, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ,
                 nperseg_s=PSD_NPERSEG_S) -> float:
    """How far (Hz) from a target bin a narrowband line can still change that bin's SNR (#309).

    :func:`snr_spectrum` divides a bin's power by the mean of the bins ``n_guard + 1`` to
    ``n_guard + n_neigh`` either side of it, so a line anywhere in that neighbourhood inflates the
    noise estimate and pulls the SNR *down*, exactly as a line in the signal bin pushes it up. And a
    line is not one bin wide: through the Welch segment's Hann window it spreads over
    :data:`WELCH_MAIN_LOBE_HALF_BINS` bins either side of its own, so a line just outside the
    neighbourhood still leaks into its outermost bin. The reach is therefore
    ``(n_neigh + n_guard + WELCH_MAIN_LOBE_HALF_BINS) * bin_w``, with ``n_neigh``, ``n_guard`` and
    the bin width ``sfreq / nperseg`` computed exactly as :func:`compute_psd` and
    :func:`snr_spectrum` compute them for an ``n_times``-sample input: derived from the estimator's
    own parameters, so changing them moves the guard with them. At 250 Hz with 4 s trials that is
    (4 + 1 + 2) * 0.25 = 1.75 Hz; the decoder's own guard tolerance, ``1/T``, is 0.25 Hz there,
    which is why the SNR guard is wider than the decoder's.

    It covers the main lobe, not the side lobes: an off-bin line tens of dB above the noise floor
    still leaks a smaller, decaying error beyond this, which no finite guard removes.

    **Not monotonic in ``n_times``** (#362 review): ``n_neigh``/``n_guard`` round to whole bins, so a
    shorter window can occasionally come out *narrower* than a longer one (e.g. at 250 Hz: 1.75 Hz at
    1000 samples, 1.667 Hz at 600, 2.0 Hz at 500) — a shorter window does not always widen the
    tolerance. Harmless for the target-SNR DV: each target's reach is evaluated at the same ``n_times``
    its own SNR is summed over (:func:`ssvep.analysis.pipeline._snr_guard_by_target`), so the two never
    disagree about which window they describe.
    """
    nperseg = psd_nperseg(sfreq, n_times, nperseg_s)
    if nperseg < 1:
        return float("inf")
    bin_w = float(sfreq) / float(nperseg)
    n_neigh, n_guard = snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
    return (n_neigh + n_guard + WELCH_MAIN_LOBE_HALF_BINS) * bin_w

float_snr_guard_tolerance_hz

float_snr_guard_tolerance_hz(sfreq, n_times)

:func:snr_guard_tolerance_hz with the reach's kernel sized by the float arithmetic every record before #392 was written with (:func:_float_kernel_bins of sfreq / nperseg). The D7 key pins an older record's window by reproducing its tolerance_hz exactly (:func:ssvep.analysis.contrast.snr_edge_key), so it must use the arithmetic that wrote it — the same number as :func:snr_guard_tolerance_hz wherever :func:reach_kernel_moved is false, which is every window measured (:func:reach_kernel_moved).

Source code in src/ssvep/classification/spectral.py

def float_snr_guard_tolerance_hz(sfreq, n_times) -> float:
    """:func:`snr_guard_tolerance_hz` with the reach's kernel sized by the float arithmetic every
    record before #392 was written with (:func:`_float_kernel_bins` of ``sfreq / nperseg``). The D7
    key pins an older record's window by reproducing its ``tolerance_hz`` exactly
    (:func:`ssvep.analysis.contrast.snr_edge_key`), so it must use the arithmetic that wrote it —
    the same number as :func:`snr_guard_tolerance_hz` wherever :func:`reach_kernel_moved` is false,
    which is every window measured (:func:`reach_kernel_moved`)."""
    from . import decoding
    nperseg = psd_nperseg(sfreq, n_times)
    if nperseg < 1:
        return float("inf")
    bin_w = float(sfreq) / float(nperseg)
    n_neigh, n_guard = _float_kernel_bins(bin_w)
    reach = (n_neigh + n_guard + WELCH_MAIN_LOBE_HALF_BINS) * bin_w
    return float(max(sfreq / n_times, decoding.MAINS_MIN_GUARD_HZ, reach))

snr_guard_tolerance_hz

snr_guard_tolerance_hz(sfreq, n_times)

How close to a mains line a target's harmonic counts as colliding, for the target-SNR DV's guard at an n_times-sample window (#309, #362): the widest of the decoder's own resolution sfreq / n_times, its floor :data:~ssvep.classification.decoding.MAINS_MIN_GUARD_HZ, and the SNR estimator's reach (:func:snr_reach_hz). What target_snr_harmonics.tolerance_hz records per target, named so a record's reader (:func:ssvep.analysis.contrast.snr_edge_key,

387) can recompute it with the identical arithmetic.

Source code in src/ssvep/classification/spectral.py

def snr_guard_tolerance_hz(sfreq, n_times) -> float:
    """How close to a mains line a target's harmonic counts as colliding, for the target-SNR DV's
    guard at an ``n_times``-sample window (#309, #362): the widest of the decoder's own resolution
    ``sfreq / n_times``, its floor :data:`~ssvep.classification.decoding.MAINS_MIN_GUARD_HZ`, and the
    SNR estimator's reach (:func:`snr_reach_hz`). What ``target_snr_harmonics.tolerance_hz``
    records per target, named so a record's reader (:func:`ssvep.analysis.contrast.snr_edge_key`,
    #387) can recompute it with the identical arithmetic."""
    from . import decoding
    return float(max(sfreq / n_times, decoding.MAINS_MIN_GUARD_HZ, snr_reach_hz(sfreq, n_times)))

snr_spectrum_db

snr_spectrum_db(data, sfreq, lo_hz, hi_hz, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, require_full_neighbours=False)

Channel-averaged SNR spectrum (dB) over [lo_hz, hi_hz] — one curve vs frequency, from the same :func:snr_spectrum kernel every other SNR number in this module comes from (#21, #53).

data is (n_ch, n_times) — one window (the live Run-tab pane's own averaging tail) — or (n_trials, n_ch, n_times), averaged over the leading trial axis after the PSD (power, not amplitude, averages correctly across trials of unrelated phase) exactly as :func:ssvep.analysis.pipeline.mean_snr_spectrum always has. Unlike every other function in this module, a bare (n_times,) array does not work here: np.nanmean(snr, axis=0) always treats the leading axis as channels to average over, so a 1-D input silently collapses frequency to a scalar instead (and then raises IndexError indexing it). Pass a single channel as (1, n_times) instead.

The returned bins are independent of which band is requested: :func:compute_psd's bin grid is fixed by sfreq and the window length alone, and [lo_hz, hi_hz] only selects which already-computed bins to keep. So a live pane showing a wide band (to fit the mains line and every harmonic marker alongside the PSD view it sits next to) and an offline report showing a narrow one (just the targets ± a margin) agree bin-for-bin wherever the two overlap, at the same window length (the bin grid depends on it too) — this is the mechanism behind "live and offline SNR agree by construction" (#21), not a coincidence of matching parameters.

The kernel is :func:snr_kernel_bins's, from sfreq and the window (#392) — not from the fetched sub-band's own first two bins, whose float spacing depends on where fmin cut the grid and, at a rounding-tie window (1.5, 2.5 or 3.5 s), used to size the live call's kernel (fetch floor 1 Hz) differently from the report's (fetch from 0 Hz or lo_hz - pad). So the two agree bin for bin with each other and with the target-SNR DV's own full-grid spectrum at every window.

The analysed PSD is padded by the neighbour/guard kernel's half-width before :func:snr_spectrum, then trimmed back to [lo_hz, hi_hz] after — near a spectrum edge the neighbour kernel would otherwise run off the array and the estimate would be invalid (spuriously high). require_full_neighbours=True (the live Run-tab view's own call, #380) keeps its existing 1.0 Hz fetch floor (:func:compute_psd's own default fmin; the target-SNR path, :func:target_snr_orders_from_data, used the same floor until #387) and drops any bin that does not end up with a full kernel of real neighbours on both sides — at the low-frequency edge the floor creates, and symmetrically at the top edge when hi_hz + pad runs into Nyquist — instead of returning it. Its own numbers are unchanged by anything below (#381 is scoped to the offline path only; a prior review already fixed the live one).

require_full_neighbours=False (the default, the offline report's own call) instead pads from true 0 Hz (fmin = max(0.0, lo_hz - pad)) rather than flooring at 1.0 Hz: :func:compute_psd's bins are real all the way down to DC, so the 1.0 Hz floor was never a data limit, just a convention this kernel used to share with the target-SNR path (fixed there too, #387) — and sharing it here silently zero-padded the low edge instead of giving it the real data that was available all along (#381; measured on pink noise with no response: a 1.25 Hz bin read 14.3 dB zero-padded, 3.7 dB with its real neighbours). That recovers a full kernel for every bin except the few nearest the two edges no amount of padding can supply real data for: below 0 Hz, and above Nyquist. Those — margin = round(neighbor_hz / bin_w) + round(guard_hz / bin_w) bins at each edge (checking whether the fetch itself reached 0 Hz/Nyquist before dropping is wrong below about a 0.6 s window — see the code comment — so this is applied the same way regardless of window length) — are dropped the same way require_full_neighbours drops them, rather than returned zero-padded and unmarked. The dropped bins are only the margin bins nearest 0 Hz/Nyquist — about 1.25 Hz at a 4 s window, so no shipped target (minimum 4 Hz) loses a bin. A wider band of bins changes value without being dropped: any bin whose own kernel neighbourhood used to reach the old 1.0 Hz floor now gets real data instead, which is every bin below about 1.0 + margin — roughly 2.25 Hz at a 4 s window, further at shorter windows where margin is wider in Hz. The shipped 4 Hz sweep's own target bins sit above that and are bit-identical (this display curve is not the target-SNR path, whose own rounding-tie exception :func:snr_edge_status describes); what changes for it is only the unlabelled display curve near 1–2 Hz, where no target sits, which no longer shows an inflated artefact.

Source code in src/ssvep/classification/spectral.py

def snr_spectrum_db(data, sfreq, lo_hz, hi_hz, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ,
                    require_full_neighbours=False):
    """Channel-averaged SNR spectrum (dB) over ``[lo_hz, hi_hz]`` — one curve vs frequency, from the
    same :func:`snr_spectrum` kernel every other SNR number in this module comes from (#21, #53).

    ``data`` is ``(n_ch, n_times)`` — one window (the live Run-tab pane's own averaging tail) — or
    ``(n_trials, n_ch, n_times)``, averaged over the leading trial axis *after* the PSD (power, not
    amplitude, averages correctly across trials of unrelated phase) exactly as
    :func:`ssvep.analysis.pipeline.mean_snr_spectrum` always has. Unlike every other function in
    this module, a bare ``(n_times,)`` array does not work here: ``np.nanmean(snr, axis=0)`` always
    treats the leading axis as channels to average over, so a 1-D input silently collapses frequency
    to a scalar instead (and then raises ``IndexError`` indexing it). Pass a single channel as
    ``(1, n_times)`` instead.

    The returned bins are **independent of which band is requested**: :func:`compute_psd`'s bin grid
    is fixed by ``sfreq`` and the window length alone, and ``[lo_hz, hi_hz]`` only selects which
    already-computed bins to keep. So a live pane showing a wide band (to fit the mains line and
    every harmonic marker alongside the PSD view it sits next to) and an offline report showing a
    narrow one (just the targets ± a margin) agree bin-for-bin wherever the two overlap, **at the
    same window length** (the bin grid depends on it too) — this is the mechanism behind "live and
    offline SNR agree by construction" (#21), not a coincidence of matching parameters.

    The kernel is :func:`snr_kernel_bins`'s, from ``sfreq`` and the window (#392) — not from the
    fetched sub-band's own first two bins, whose float spacing depends on where ``fmin`` cut the grid
    and, at a rounding-tie window (1.5, 2.5 or 3.5 s), used to size the live call's kernel (fetch floor
    1 Hz) differently from the report's (fetch from 0 Hz or ``lo_hz - pad``). So the two agree bin
    for bin with each other and with the target-SNR DV's own full-grid spectrum at every window.

    The analysed PSD is padded by the neighbour/guard kernel's half-width before
    :func:`snr_spectrum`, then trimmed back to ``[lo_hz, hi_hz]`` after — near a spectrum edge the
    neighbour kernel would otherwise run off the array and the estimate would be invalid (spuriously
    high). ``require_full_neighbours=True`` (the live Run-tab view's own call, #380) keeps its
    existing 1.0 Hz fetch floor (:func:`compute_psd`'s own default ``fmin``; the target-SNR path,
    :func:`target_snr_orders_from_data`, used the same floor until #387) and drops any bin that
    does not end up with a full kernel of real neighbours on both sides — at the low-frequency edge
    the floor creates, and symmetrically at the top edge when ``hi_hz + pad`` runs into Nyquist —
    instead of returning it. Its own numbers are unchanged by anything below (#381 is scoped to the offline
    path only; a prior review already fixed the live one).

    ``require_full_neighbours=False`` (the default, the offline report's own call) instead pads
    from **true 0 Hz** (``fmin = max(0.0, lo_hz - pad)``) rather than flooring at 1.0 Hz:
    :func:`compute_psd`'s bins are real all the way down to DC, so the 1.0 Hz floor was never a
    data limit, just a convention this kernel used to share with the target-SNR path (fixed there
    too, #387) — and sharing it here silently zero-padded the low edge instead of giving it the
    real data that was available all along (#381; measured on pink noise with no response: a 1.25 Hz bin read
    14.3 dB zero-padded, 3.7 dB with its real neighbours). That recovers a full kernel for every bin
    except the few nearest the two edges no amount of padding can supply real data for: below 0 Hz,
    and above Nyquist. Those — ``margin = round(neighbor_hz / bin_w) + round(guard_hz / bin_w)``
    bins at *each* edge (checking whether the fetch itself reached 0 Hz/Nyquist before dropping is
    wrong below about a 0.6 s window — see the code comment — so this is applied the same way
    regardless of window length) — are dropped the same way ``require_full_neighbours`` drops them,
    rather than returned zero-padded and unmarked. The *dropped* bins are only the ``margin`` bins
    nearest 0 Hz/Nyquist — about 1.25 Hz at a 4 s window, so no shipped target (minimum 4 Hz) loses
    a bin. A wider band of bins *changes value* without being dropped: any bin whose own kernel
    neighbourhood used to reach the old 1.0 Hz floor now gets real data instead, which is every bin
    below about 1.0 + margin — roughly 2.25 Hz at a 4 s window, further at shorter windows where
    ``margin`` is wider in Hz. The shipped 4 Hz sweep's own target bins sit above that and are
    bit-identical (this display curve is not the target-SNR path, whose own rounding-tie
    exception :func:`snr_edge_status` describes); what changes for it is only the unlabelled *display* curve near 1–2 Hz, where no
    target sits, which no longer shows an inflated artefact.
    """
    pad = neighbor_hz + guard_hz + 0.5                           # kernel half-width + a bin
    # The live call keeps its own 1.0 Hz fetch floor exactly as before (#380); the offline default
    # fetches from true 0 Hz instead, per the docstring above.
    fetch_floor = 1.0 if require_full_neighbours else 0.0
    fmin = max(fetch_floor, lo_hz - pad)
    fx, psd = compute_psd(data, sfreq, fmin=fmin, fmax=hi_hz + pad)
    if psd.ndim == 3:
        psd = psd.mean(axis=0)                                   # avg trials -> (n_ch, n_freq)
    # The kernel from sfreq and the window, exactly (#392) — not from this sub-band's own first
    # two bins, whose float spacing depends on where `fmin` cut the grid.
    n_times = np.shape(data)[-1]
    kernel = snr_kernel_bins(sfreq, n_times, neighbor_hz, guard_hz) if n_times >= 1 else None
    snr = snr_spectrum(psd, fx, neighbor_hz=neighbor_hz, guard_hz=guard_hz, kernel=kernel)
    with warnings.catch_warnings(), np.errstate(invalid="ignore"):
        # An idle/all-zero board (or any bin with no finite channel) makes a whole column of `snr`
        # all-NaN; np.nanmean's own "Mean of empty slice" plus log10's "invalid value" would
        # otherwise spam the console at every recompute (#21 review) — the result is correctly NaN
        # either way, just quietly.
        warnings.simplefilter("ignore", category=RuntimeWarning)
        snr_db = 10.0 * np.log10(np.nanmean(snr, axis=0))         # avg channels -> (n_freq,)
    keep = (fx >= lo_hz) & (fx <= hi_hz)
    if len(fx) >= 2:
        n_neigh, n_guard = kernel
        margin = n_neigh + n_guard
        n_bins = len(fx)
        idx = np.arange(n_bins)
        # Applied unconditionally at both edges, for both callers (#381 review finding 3): a bin
        # still needs margin real neighbours on each side regardless of whether the fetch happened
        # to land exactly on 0 Hz/Nyquist. Checking "did the fetch hit the edge" instead (as a
        # previous version of this function did for the offline path) is wrong whenever the kernel
        # margin in Hz exceeds `pad` — any window under ~0.6 s — because then the fetch never
        # reaches the true edge at all, the check never fires, and the lowest/highest kept bin is
        # still short real neighbours. For every window length where the fetch does reach (or pass)
        # the edge, this is identical to the old conditional check, so ordinary ≥1 s windows are
        # unaffected.
        if 2 * margin < n_bins:
            keep = keep & (idx >= margin) & (idx <= n_bins - 1 - margin)
        else:                                          # no bin anywhere has room for a full kernel
            keep = keep & np.zeros(n_bins, dtype=bool)
    # Provenance marker, not a D7 pooling key: every call from here on sets it so a ``metrics.json``
    # can be told apart from one written before #381, but no pooled DV ever reads this block (see
    # #381's PR discussion) — ``report._fig_snr_spectrum`` is its only reader.
    return {"freqs_hz": [float(v) for v in fx[keep]], "snr_db": [float(v) for v in snr_db[keep]],
            "edge_fix": True}

snr_safe_harmonics

snr_safe_harmonics(f, n_harmonics, n_times, sfreq, line_hz, neighbor_hz=SNR_NEIGHBOR_HZ, guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S)

The harmonic orders (1-based) target f's own SNR sum can trust when the recording carries mains at line_hz (#309); () when not even its fundamental can be trusted.

The decoder's collision rule (:func:ssvep.classification.decoding.mains_safe_harmonics: a harmonic, as sampled, near a mains line or one of its harmonics below Nyquist), with two differences:

• Per target. Called on [f] alone, so only f's own collisions affect f. The decoder drops a colliding order from every candidate because its argmax compares candidates; one target's SNR is computed on its own. (A contrast that compares SNR between arms does need one order set for both, and builds it from the per-order record: :func:ssvep.analysis.contrast.common_snr_orders.)
• The estimator's reach, not the decoder's resolution. The tolerance is at least :func:snr_reach_hz (a line in the harmonic's noise neighbourhood corrupts its SNR as much as one in its signal bin), as well as the decoder's own max(1/T, MAINS_MIN_GUARD_HZ).

Shared by the offline pipeline and the live readout, so both drop the same orders.

Source code in src/ssvep/classification/spectral.py

def snr_safe_harmonics(f, n_harmonics, n_times, sfreq, line_hz, neighbor_hz=SNR_NEIGHBOR_HZ,
                       guard_hz=SNR_GUARD_HZ, nperseg_s=PSD_NPERSEG_S) -> tuple:
    """The harmonic orders (1-based) target ``f``'s own SNR sum can trust when the recording
    carries mains at ``line_hz`` (#309); ``()`` when not even its fundamental can be trusted.

    The decoder's collision rule (:func:`ssvep.classification.decoding.mains_safe_harmonics`: a
    harmonic, as sampled, near a mains line or one of its harmonics below Nyquist), with two
    differences:

    * **Per target.** Called on ``[f]`` alone, so only ``f``'s own collisions affect ``f``. The
      decoder drops a colliding order from every candidate because its argmax compares candidates;
      one target's SNR is computed on its own. (A *contrast* that compares SNR between arms does
      need one order set for both, and builds it from the per-order record:
      :func:`ssvep.analysis.contrast.common_snr_orders`.)
    * **The estimator's reach, not the decoder's resolution.** The tolerance is at least
      :func:`snr_reach_hz` (a line in the harmonic's noise neighbourhood corrupts its SNR as much
      as one in its signal bin), as well as the decoder's own ``max(1/T, MAINS_MIN_GUARD_HZ)``.

    Shared by the offline pipeline and the live readout, so both drop the same orders.
    """
    from . import decoding
    reach = snr_reach_hz(sfreq, n_times, neighbor_hz, guard_hz, nperseg_s)
    orders, undecodable = decoding.mains_safe_harmonics([f], n_harmonics, n_times, sfreq, line_hz,
                                                        min_tol_hz=reach)
    return () if undecodable else tuple(orders)

Metrics

ssvep.classification.metrics

Classification metrics for SSVEP BCI: accuracy, macro precision, Wolpaw ITR.

Pure NumPy + scikit-learn. Rigor helper: always pass chance = 1/n_classes alongside any accuracy you report.

wolpaw_itr

wolpaw_itr(accuracy, n_classes, selection_time_s)

Information transfer rate in bits/min (Wolpaw et al. 1998).

selection_time_s = decision window + inter-selection overhead. B = log2(N) + Plog2(P) + (1-P)log2((1-P)/(N-1)); ITR = B*60/selection_time.

Source code in src/ssvep/classification/metrics.py

def wolpaw_itr(accuracy: float, n_classes: int, selection_time_s: float) -> float:
    """Information transfer rate in bits/min (Wolpaw et al. 1998).

    selection_time_s = decision window + inter-selection overhead.
    B = log2(N) + P*log2(P) + (1-P)*log2((1-P)/(N-1)); ITR = B*60/selection_time.
    """
    p = float(np.clip(accuracy, 0.0, 1.0))
    n = int(n_classes)
    if n <= 1 or selection_time_s <= 0:
        return 0.0
    if p >= 1.0:
        bits = np.log2(n)
    elif p <= 0.0:
        bits = 0.0
    else:
        bits = np.log2(n) + p * np.log2(p) + (1 - p) * np.log2((1 - p) / (n - 1))
    return float(max(bits, 0.0) * 60.0 / selection_time_s)

classification_metrics

classification_metrics(y_true, y_pred, n_classes, selection_time_s, labels=None)

Accuracy, macro precision, ITR, and n for one set of predictions.

Source code in src/ssvep/classification/metrics.py

def classification_metrics(y_true, y_pred, n_classes, selection_time_s, labels=None):
    """Accuracy, macro precision, ITR, and n for one set of predictions."""
    acc = accuracy_score(y_true, y_pred)
    return {
        "accuracy": float(acc),
        "precision_macro": float(precision_score(y_true, y_pred, labels=labels,
                                                 average="macro", zero_division=0)),
        "itr_bits_per_min": wolpaw_itr(acc, n_classes, selection_time_s),
        "n_trials": int(len(y_true)),
        "chance": 1.0 / n_classes,
    }

Visualisation

ssvep.classification.viz

Domain-generic SSVEP plotting primitives (matplotlib only).

• plot_channel_array: discrete per-channel map (NO interpolation) for sparse / single-region montages where interpolated topomaps mislead. Bichromatic diverging for signed values (SNR), monochromatic sequential for non-negative values (impedance).
• plot_window_sweep: mean ± SD vs decision-window length, with a chance line.
• plot_confusion: normalised confusion matrix.

plot_channel_array

plot_channel_array(positions, ch_names, values, ax=None, title=None, vlim=None, cbar=True, label='value', diverging=True, cmap=None, show_values=None, norm=None, cbar_ticks=None)

One coloured disc per channel at its 2-D montage position; no interpolation.

positions: dict ch -> (x, y) (e.g. cm from a reference landmark). Discs are drawn in data coordinates with a radius just under the nearest-neighbour spacing, so they never overlap regardless of channel count (a 64-ch cap reads as cleanly as an 8-ch montage); the figure grows with the montage. show_values prints the numeric value inside each disc — defaults to on for sparse montages (≤20 ch) and off for dense ones (colour + colourbar carry the value there); the label font shrinks with channel count so it fits the disc on a dense cap. diverging=True -> bichromatic, centred at 0 (default RdBu_r) for signed values. diverging=False -> sequential from 0 (default Greens) for non-negative values. norm (+ a cmap that may be a Colormap object) overrides the auto scale — e.g. a BoundaryNorm + ListedColormap for discrete pass/marginal/fail zones; cbar_ticks then labels the boundaries.

Source code in src/ssvep/classification/viz.py

def plot_channel_array(positions, ch_names, values, ax=None, title=None, vlim=None,
                       cbar=True, label="value", diverging=True, cmap=None, show_values=None,
                       norm=None, cbar_ticks=None):
    """One coloured disc per channel at its 2-D montage position; no interpolation.

    positions: dict ``ch -> (x, y)`` (e.g. cm from a reference landmark). Discs are drawn in **data
    coordinates** with a radius just under the nearest-neighbour spacing, so they never overlap
    regardless of channel count (a 64-ch cap reads as cleanly as an 8-ch montage); the figure grows
    with the montage. ``show_values`` prints the numeric value inside each disc — defaults to on for
    sparse montages (≤20 ch) and off for dense ones (colour + colourbar carry the value there); the
    label font shrinks with channel count so it fits the disc on a dense cap.
    diverging=True  -> bichromatic, centred at 0 (default RdBu_r) for signed values.
    diverging=False -> sequential from 0 (default Greens) for non-negative values.
    ``norm`` (+ a ``cmap`` that may be a Colormap object) overrides the auto scale — e.g. a
    ``BoundaryNorm`` + ``ListedColormap`` for discrete pass/marginal/fail zones; ``cbar_ticks`` then
    labels the boundaries.
    """
    vals = np.asarray(values, float)
    finite = vals[np.isfinite(vals)]
    if norm is not None:
        cmap = cmap or "Greens"                                 # caller supplied the scale; use as-is
    elif diverging:
        cmap = cmap or "RdBu_r"
        m = (np.nanmax(np.abs(finite)) if finite.size else 1.0) if vlim is None else max(abs(vlim[0]), abs(vlim[1]))
        norm = TwoSlopeNorm(vcenter=0.0, vmin=-m or -1e-6, vmax=m or 1e-6)
    else:
        cmap = cmap or "Greens"
        vlim = vlim or (0.0, np.nanmax(finite) if finite.size else 1.0)
        norm = Normalize(vmin=vlim[0], vmax=max(vlim[1], vlim[0] + 1e-6))
    cmap_obj = cmap if isinstance(cmap, Colormap) else plt.get_cmap(cmap)
    xs = np.array([positions[c][0] for c in ch_names], float)
    ys = np.array([positions[c][1] for c in ch_names], float)
    n = len(ch_names)
    if show_values is None:
        show_values = n <= 20

    # Declutter for display: some montages place electrodes almost on top of each other (the actiCAP
    # schematic puts TP9/TP7 ~0.6 cm apart), and one near-coincident pair would force a tiny disc radius
    # for the WHOLE map. Nudge only the too-close pairs apart to a floor of 0.8×(typical spacing); a
    # well-spaced montage (8-ch) is untouched. Positions here are display-only — the montage is unchanged.
    xs, ys = _declutter(xs, ys, 0.8 * _typical_spacing(xs, ys))
    radius = 0.46 * _nearest_neighbour_dist(xs, ys)            # data units → discs (near-)touch at most
    if ax is None:
        span_x = (xs.max() - xs.min()) + 3 * radius
        span_y = (ys.max() - ys.min()) + 3 * radius
        base = float(np.clip(4.6 + 0.05 * n, 4.6, 9.0))         # bigger canvas for denser montages
        big = max(span_x, span_y)
        _, ax = plt.subplots(figsize=(base * span_x / big + (1.2 if cbar else 0.3),
                                      base * span_y / big + 0.6))

    faces = [cmap_obj(norm(v)) if np.isfinite(v) else (0.85, 0.85, 0.85, 1.0) for v in vals]
    ax.add_collection(PatchCollection([Circle((x, y), radius) for x, y in zip(xs, ys)],
                                      facecolors=faces, edgecolors="k", linewidths=0.8, zorder=2))

    # Shrink the label font as the montage densifies so it fits inside the (smaller) discs. The old
    # 200/√n constant saturated at the 9 pt cap for every montage ≤493 ch, so a 64-ch cap got the same
    # 9 pt as an 8-ch one and the 10-10 names overflowed. 40/√n gives ~9 pt at 8-16 ch, ~7 at 32, ~5 at 64.
    fontsize = float(np.clip(40.0 / np.sqrt(max(n, 1)), 4.5, 9.0))
    for c, x, y, v, face in zip(ch_names, xs, ys, vals, faces):
        r, g, b = face[:3]
        tcol = "white" if np.isfinite(v) and (0.299 * r + 0.587 * g + 0.114 * b) < 0.5 else "black"
        if show_values:
            txt = f"{c}\n{v:.1f}" if np.isfinite(v) else f"{c}\nn/a"
        else:
            txt = c
        ax.text(x, y, txt, ha="center", va="center", fontsize=fontsize, color=tcol, zorder=3)

    pad = 1.5 * radius
    ax.set_xlim(xs.min() - pad, xs.max() + pad)
    ax.set_ylim(ys.min() - pad, ys.max() + pad)
    ax.set_aspect("equal")
    ax.set_xticks([]); ax.set_yticks([])
    if title:
        ax.set_title(title, fontsize=10)
    if cbar:
        sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap_obj)
        sm.set_array([])
        cb = ax.figure.colorbar(sm, ax=ax, fraction=0.046, label=label,   # figure-safe (embedded canvases)
                                ticks=cbar_ticks)
        if cbar_ticks is not None:
            cb.ax.tick_params(labelsize=8)
    return ax.figure

impedance_zones

impedance_zones(board, electrode_type=None)

The traffic-light zones for one board + electrode prep, or None for the fixed gradient.

actiCHamp has its own; every Cyton variant (8-ch, Daisy, WiFi) gets the Cyton zones whatever the electrode prep — calibrated for dry, borrowed for everything else (see :data:CYTON_IMPEDANCE_ZONES). A board with neither (Unicorn, synthetic, an unknown one) falls back to the fixed gradient. Pure, so the choice is testable without a widget.

Source code in src/ssvep/classification/viz.py

def impedance_zones(board: str | None, electrode_type: str | None = None):
    """The traffic-light zones for one board + electrode prep, or None for the fixed gradient.

    actiCHamp has its own; every Cyton variant (8-ch, Daisy, WiFi) gets the Cyton zones whatever the
    electrode prep — calibrated for dry, borrowed for everything else (see
    :data:`CYTON_IMPEDANCE_ZONES`). A board with neither (Unicorn, synthetic, an unknown one) falls
    back to the fixed gradient. Pure, so the choice is testable without a widget.
    """
    board = (board or "").lower()
    if board == "actichamp":
        return ACTICHAMP_IMPEDANCE_ZONES
    if board.startswith("cyton"):
        # Kept as two names on purpose, though they are the same bounds today: giving gold-cup its
        # own calibration is a change to CYTON_IMPEDANCE_ZONES alone, and dry must not move with it.
        return (CYTON_DRY_IMPEDANCE_ZONES if "dry" in (electrode_type or "").lower()
                else CYTON_IMPEDANCE_ZONES)
    return None

zone_bounds

zone_bounds(zones)

Colour-bar bounds for zones: [0, t1, …, tn, top], a fixed function of the zones alone.

The last band is open-ended (inf), so the bar needs a finite top to draw it; it is the last threshold plus the width of the band below it, so the red band is drawn as wide as the yellow one. It never depends on the data: a top that followed the worst channel slid the threshold ticks along the bar on every update (#564).

Source code in src/ssvep/classification/viz.py

def zone_bounds(zones) -> list[float]:
    """Colour-bar bounds for ``zones``: ``[0, t1, …, tn, top]``, a fixed function of the zones alone.

    The last band is open-ended (``inf``), so the bar needs a finite top to draw it; it is the last
    threshold plus the width of the band below it, so the red band is drawn as wide as the yellow one.
    It never depends on the data: a top that followed the worst channel slid the threshold ticks
    along the bar on every update (#564).
    """
    ups = [float(u) for u, _ in zones][:-1]           # the finite thresholds, e.g. [30, 60]
    if not ups:
        return [0.0, IMPEDANCE_VMAX_KOHM]
    below = ups[-2] if len(ups) >= 2 else 0.0
    return [0.0] + ups + [ups[-1] + (ups[-1] - below)]

plot_impedance

plot_impedance(positions, ch_names, kohm, ax=None, vmax=None, title='Electrode impedance', cmap=IMPEDANCE_CMAP, zones=None)

Discrete channel-array map of impedance (kΩ).

Thin wrapper over :func:plot_channel_array (sequential, non-negative). kohm aligns to ch_names; NaN renders as an 'n/a' (unmeasured) marker.

zones (e.g. :func:impedance_zones) gives discrete traffic-light bands: a (upper_kohm, colour) list, worst bound last. A band includes its lower threshold and excludes its upper one (exactly 30 kΩ is yellow on actiCHamp), and the colourbar shows the bands with ticks at the thresholds. vmax is ignored when zones is given.

Without zones it is a green→red gradient from 0 to vmax (default :data:IMPEDANCE_VMAX_KOHM), clipped above. The scale is never derived from the data, either way (#564): the same colour must mean the same kΩ for the whole check.

Source code in src/ssvep/classification/viz.py

def plot_impedance(positions, ch_names, kohm, ax=None, vmax=None, title="Electrode impedance",
                   cmap=IMPEDANCE_CMAP, zones=None):
    """Discrete channel-array map of impedance (kΩ).

    Thin wrapper over :func:`plot_channel_array` (sequential, non-negative). ``kohm`` aligns to
    ``ch_names``; NaN renders as an 'n/a' (unmeasured) marker.

    ``zones`` (e.g. :func:`impedance_zones`) gives **discrete traffic-light bands**: a
    ``(upper_kohm, colour)`` list, worst bound last. A band includes its lower threshold and excludes
    its upper one (exactly 30 kΩ is yellow on actiCHamp), and the colourbar shows the bands with ticks
    at the thresholds. ``vmax`` is ignored when ``zones`` is given.

    Without zones it is a green→red gradient from 0 to ``vmax`` (default
    :data:`IMPEDANCE_VMAX_KOHM`), clipped above. **The scale is never derived from the data**,
    either way (#564): the same colour must mean the same kΩ for the whole check.
    """
    if zones:
        listed = ListedColormap([c for _, c in zones])
        bounds = zone_bounds(zones)
        norm = BoundaryNorm(bounds, listed.N, clip=True)        # >top → clipped into the last (red) band
        return plot_channel_array(positions, ch_names, kohm, ax=ax, title=title,
                                  label="impedance (kΩ)", diverging=False, cmap=listed, norm=norm,
                                  cbar_ticks=bounds[1:-1])       # tick the real thresholds (30, 60)
    vmax = IMPEDANCE_VMAX_KOHM if vmax is None else float(vmax)
    norm = Normalize(vmin=0.0, vmax=vmax, clip=True)
    return plot_channel_array(positions, ch_names, kohm, ax=ax, title=title,
                              label="impedance (kΩ)", diverging=False, cmap=cmap, norm=norm)

plot_window_sweep

plot_window_sweep(windows, series, ylabel, chance=None, ax=None, ypct=False, title=None, ylim=None)

series: dict label -> (mean_array, sd_array) over windows.

ylim fixes the y-axis to (lo, hi) so plots are comparable across reports; None keeps matplotlib's autoscale. An auto-scaled axis makes every run's peak sit at the same height even when the runs differ by tens of bits/min, so the reader can't eyeball which run was better — the report layer passes a fixed range for exactly that reason.

Source code in src/ssvep/classification/viz.py

def plot_window_sweep(windows, series, ylabel, chance=None, ax=None, ypct=False, title=None,
                      ylim=None):
    """series: dict ``label -> (mean_array, sd_array)`` over ``windows``.

    ``ylim`` fixes the y-axis to ``(lo, hi)`` so plots are comparable across reports; ``None`` keeps
    matplotlib's autoscale. An auto-scaled axis makes every run's peak sit at the same height even
    when the runs differ by tens of bits/min, so the reader can't eyeball which run was better — the
    report layer passes a fixed range for exactly that reason.
    """
    if ax is None:
        _, ax = plt.subplots(figsize=(6.5, 4.3))
    for lab, (mean, sd) in series.items():
        mean, sd = np.asarray(mean, float), np.asarray(sd, float)
        if ypct:
            mean, sd = mean * 100, sd * 100
        line, = ax.plot(windows, mean, "-o", label=lab, zorder=3)
        ax.fill_between(windows, mean - sd, mean + sd, color=line.get_color(), alpha=0.2)
    if chance is not None:
        ax.axhline(chance * 100 if ypct else chance, color="grey", ls="--", lw=1,
                   label=f"chance ({chance*100:.1f}%)")
    ax.set_xlabel("Stimulation-window length (s)")
    ax.set_ylabel(ylabel + (" (%)" if ypct else ""))
    if ylim is not None:
        ax.set_ylim(*ylim)
    if title:
        ax.set_title(title, fontsize=10)
    ax.legend(fontsize=8); ax.grid(alpha=0.3)
    return ax.figure