diff --git a/retina_tracker/calibrate.py b/retina_tracker/calibrate.py
new file mode 100644
index 0000000..bb25051
--- /dev/null
+++ b/retina_tracker/calibrate.py
@@ -0,0 +1,333 @@
+"""Fit R and Q from recorded innovations, separating two constants that cancel.
+
+R and Q are both wrong today, in opposite directions. R is around 45x too large
+in delay, while the model error it is standing in for grows as the square of the
+prediction interval rather than the fifth power white jerk implies. The sum comes
+out plausible, which is why nothing has complained, and it is why neither can be
+fitted alone: a single scalar like NIS moves one way for an oversized R and the
+other for an undersized Q.
+
+What separates them is the prediction interval. The innovation covariance is
+HPH' + R, so scatter measured against how long the filter had been predicting is
+a constant plus a growing term:
+
+    var(nu | L) = R + a * L^p
+
+R is the intercept, and (a, p) are the shape and size of the model error. A grid
+over p makes the rest a two-parameter linear fit, which is all this needs.
+
+Deliberately stdlib-only, like live_score: it is meant to be copied onto a node
+and run under whatever python3 is there.
+
+    python -m retina_tracker.calibrate innovations.jsonl --blah2-config config.yml
+    python -m retina_tracker.calibrate jn1.jsonl ffb.jsonl --fc 177000000
+"""
+
+import argparse
+import json
+import math
+import statistics
+import sys
+from collections import defaultdict
+
+SPEED_OF_LIGHT = 299792458.0
+
+# Median of chi-squared with two degrees of freedom. A matched filter's NIS
+# should sit here; the plan measured 1.083 live, which is the signature of the
+# oversized R holding the whole thing down.
+CHI2_2DOF_MEDIAN = 2.0 * math.log(2.0)
+
+# 1.4826 * MAD estimates sigma for a Gaussian, and unlike a standard deviation
+# it does not let one mis-association set the answer. Calibration data comes
+# from a live node, so it has mis-associations in it.
+MAD_TO_SIGMA = 1.4826
+
+# Prediction intervals are quantised by the frame rate anyway; rounding groups
+# the ones that differ only by clock jitter.
+INTERVAL_QUANTUM_S = 0.05
+
+MIN_PER_INTERVAL = 30
+MIN_INTERVALS = 3
+
+# A track whose Doppler never leaves one bin is the fallback test for
+# interference, used only on records written before the map carried its verdict.
+FIXED_DOPPLER_BINS = 1.0
+MIN_RECORDS_TO_JUDGE_A_TRACK = 8
+
+
+def wavelength_km(fc_hz):
+    return SPEED_OF_LIGHT / fc_hz / 1000.0
+
+
+def delay_cell_km(fs_hz):
+    return SPEED_OF_LIGHT / fs_hz / 1000.0
+
+
+def load_blah2_capture(path):
+    """Read fc, fs and cpi out of a blah2 config without a YAML dependency."""
+    found = {}
+    section = []
+    with open(path) as handle:
+        for line in handle:
+            if not line.strip() or line.lstrip().startswith("#"):
+                continue
+            indent = len(line) - len(line.lstrip())
+            stripped = line.strip()
+            depth = indent // 2
+            section = section[:depth]
+            if stripped.endswith(":"):
+                section.append(stripped[:-1])
+                continue
+            key, _, value = stripped.partition(":")
+            value = value.strip()
+            if section[:1] == ["capture"] and key in ("fc", "fs"):
+                found[key] = float(value)
+            elif section[:2] == ["process", "data"] and key == "cpi":
+                found["cpi"] = float(value)
+    return found
+
+
+def load_records(path):
+    records = []
+    with open(path) as handle:
+        for line in handle:
+            line = line.strip()
+            if line:
+                records.append(json.loads(line))
+    return records
+
+
+def interfering_tracks(records):
+    """Tracks whose Doppler never left one bin, for records with no verdict.
+
+    Only a fallback. It costs the genuine constant-Doppler aircraft along with
+    the tone, which biases what is left toward manoeuvring targets, so it is
+    worth something for R and worth less for Q. Records written by a tracker
+    carrying the occupancy map say so per detection instead.
+    """
+    by_track = defaultdict(list)
+    for record in records:
+        doppler = record.get("doppler")
+        if doppler is not None:
+            by_track[record.get("birth")].append(doppler)
+    return {
+        birth
+        for birth, dopplers in by_track.items()
+        if len(dopplers) >= MIN_RECORDS_TO_JUDGE_A_TRACK and max(dopplers) - min(dopplers) < FIXED_DOPPLER_BINS
+    }
+
+
+def usable(records):
+    """Records that describe a real measurement of a real target.
+
+    Drops the interference, which is most of the record at an interfered node,
+    and the degenerate updates, whose innovation is NaN because a singular
+    covariance meant the measurement was skipped.
+    """
+    fallback = interfering_tracks(records) if not any("interfering" in r for r in records) else set()
+    kept = []
+    for record in records:
+        if record.get("interfering") or record.get("birth") in fallback:
+            continue
+        innovation = record.get("innovation") or []
+        if len(innovation) != 2 or any(x is None or not math.isfinite(x) for x in innovation):
+            continue
+        if not record.get("dt"):
+            continue
+        kept.append(record)
+    return kept
+
+
+def prediction_interval(record):
+    """How long the filter had been predicting when this innovation was taken.
+
+    dt is the interval of the last predict and n_missed the frames coasted
+    before this association landed, so this assumes a steady frame rate across
+    a coast. It is the quantity R and Q are separated along, and getting it
+    from the record is the whole reason the tracker emits n_missed.
+    """
+    return record["dt"] * (1 + record.get("n_missed", 0))
+
+
+def mad_sigma(values):
+    centre = statistics.median(values)
+    return MAD_TO_SIGMA * statistics.median([abs(v - centre) for v in values])
+
+
+def group_by_interval(records, component):
+    """Scatter of one innovation component against prediction interval."""
+    buckets = defaultdict(list)
+    for record in records:
+        interval = round(prediction_interval(record) / INTERVAL_QUANTUM_S) * INTERVAL_QUANTUM_S
+        buckets[round(interval, 3)].append(record["innovation"][component])
+    return sorted(
+        (interval, len(values), mad_sigma(values) ** 2)
+        for interval, values in buckets.items()
+        if interval > 0 and len(values) >= MIN_PER_INTERVAL
+    )
+
+
+def fit_noise_model(points, exponents=None):
+    """Least squares for var = R + a * L^p, weighted by sample count.
+
+    Linear in R and a once p is fixed, so a grid over p reduces the whole thing
+    to a 2x2 solve per candidate. p is what says whether Q has the right shape:
+    white jerk, which is what the filter implements, would be 5.
+    """
+    if len(points) < MIN_INTERVALS:
+        return None
+    if exponents is None:
+        exponents = [0.5 + 0.05 * i for i in range(91)]
+
+    best = None
+    for p in exponents:
+        s_w = s_x = s_xx = s_y = s_xy = 0.0
+        for interval, count, variance in points:
+            x = interval**p
+            w = float(count)
+            s_w += w
+            s_x += w * x
+            s_xx += w * x * x
+            s_y += w * variance
+            s_xy += w * x * variance
+        determinant = s_w * s_xx - s_x * s_x
+        if abs(determinant) < 1e-30:
+            continue
+        r = (s_y * s_xx - s_x * s_xy) / determinant
+        a = (s_w * s_xy - s_x * s_y) / determinant
+        if r < 0 or a < 0:
+            continue
+        cost = sum(count * (variance - (r + a * interval**p)) ** 2 for interval, count, variance in points)
+        if best is None or cost < best[0]:
+            best = (cost, r, a, p)
+
+    if best is None:
+        return None
+    _cost, r, a, p = best
+    return {"R": r, "a": a, "exponent": p}
+
+
+def normalised_scatter(records, component):
+    """Actual scatter over what the filter claimed, per axis.
+
+    NIS collapses both axes into one number, which is how an oversized R and an
+    undersized Q have been hiding in it. Above one means the filter is more
+    confident than it has earned on that axis; below one means the reverse.
+    """
+    ratios = [
+        record["innovation"][component] / math.sqrt(record["s_diag"][component])
+        for record in records
+        if record.get("s_diag") and record["s_diag"][component] > 0
+    ]
+    return mad_sigma(ratios) if len(ratios) >= MIN_PER_INTERVAL else None
+
+
+def calibrate(records, fc_hz, fs_hz, cpi_s):
+    kept = usable(records)
+    report = {
+        "n_records": len(records),
+        "n_used": len(kept),
+        "n_interfering": sum(1 for r in records if r.get("interfering")),
+        "n_tracks": len({r.get("birth") for r in kept}),
+        "delay_cell_km": delay_cell_km(fs_hz),
+        "doppler_bin_hz": 1.0 / cpi_s,
+        "wavelength_km": wavelength_km(fc_hz),
+        "nis_median": statistics.median([r["nis"] for r in kept if r.get("nis") is not None]) if kept else None,
+        "axes": {},
+    }
+    if not kept:
+        return report
+
+    for name, component in (("delay", 0), ("doppler", 1)):
+        fit = fit_noise_model(group_by_interval(kept, component))
+        axis = {"scatter_vs_claimed": normalised_scatter(kept, component)}
+        if fit:
+            sigma = math.sqrt(fit["R"])
+            if name == "delay":
+                axis["sigma"] = sigma
+                axis["unit"] = "km"
+                axis["cells"] = sigma / report["delay_cell_km"]
+            else:
+                # The innovation is range rate, which is what R and S are in.
+                sigma_hz = sigma / report["wavelength_km"]
+                axis["sigma"] = sigma_hz
+                axis["unit"] = "Hz"
+                axis["cells"] = sigma_hz / report["doppler_bin_hz"]
+            axis["exponent"] = fit["exponent"]
+            # At p = 2 the growing term is a velocity error held for the whole
+            # interval, which is the shape the live data showed.
+            axis["error_rate"] = math.sqrt(fit["a"])
+        report["axes"][name] = axis
+
+    return report
+
+
+def format_report(name, report):
+    lines = [f"{name}"]
+    lines.append(
+        f"  records {report['n_used']} used of {report['n_records']} "
+        f"({report['n_interfering']} interfering), {report['n_tracks']} tracks"
+    )
+    lines.append(f"  cell    {report['delay_cell_km']:.4f} km x {report['doppler_bin_hz']:.2f} Hz")
+    if report["nis_median"] is not None:
+        lines.append(f"  NIS     median {report['nis_median']:.3f}, matched {CHI2_2DOF_MEDIAN:.3f}")
+
+    for axis_name in ("delay", "doppler"):
+        axis = report["axes"].get(axis_name)
+        if not axis:
+            lines.append(f"  {axis_name:<8}not enough spread in prediction interval to separate R from Q")
+            continue
+        if "sigma" in axis:
+            lines.append(
+                f"  {axis_name:<8}R = {axis['sigma']:.4f} {axis['unit']} = {axis['cells']:.3f} cells, "
+                f"model error grows as L^{axis['exponent']:.2f} at {axis['error_rate']:.4g}"
+            )
+        if axis["scatter_vs_claimed"] is not None:
+            lines.append(f"  {'':<8}scatter is {axis['scatter_vs_claimed']:.2f}x what the filter claimed")
+    return "\n".join(lines)
+
+
+def main(argv=None):
+    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
+    parser.add_argument("files", nargs="+", help="innovations.jsonl written by --innovations")
+    parser.add_argument("--blah2-config", help="blah2 config.yml to read fc, fs and cpi from")
+    parser.add_argument("--fc", type=float, help="Centre frequency in Hz")
+    parser.add_argument("--fs", type=float, default=2000000.0, help="Sample rate in Hz (default: 2e6)")
+    parser.add_argument("--cpi", type=float, default=0.5, help="Coherent processing interval in s (default: 0.5)")
+    parser.add_argument("--json", action="store_true", help="Emit the report as JSON")
+    args = parser.parse_args(argv)
+
+    capture = load_blah2_capture(args.blah2_config) if args.blah2_config else {}
+    fc = args.fc or capture.get("fc")
+    if not fc:
+        parser.error("centre frequency unknown: pass --fc or --blah2-config")
+    fs = capture.get("fs", args.fs)
+    cpi = capture.get("cpi", args.cpi)
+
+    reports = {}
+    for path in args.files:
+        reports[path] = calibrate(load_records(path), fc, fs, cpi)
+
+    if args.json:
+        print(json.dumps(reports, indent=2))
+        return 0
+
+    for path, report in reports.items():
+        print(format_report(path, report))
+        print()
+
+    cells = [
+        report["axes"][axis]["cells"]
+        for report in reports.values()
+        for axis in ("delay", "doppler")
+        if report["axes"].get(axis, {}).get("cells") is not None
+    ]
+    if len(cells) > 1:
+        spread = max(cells) / min(cells)
+        print(f"k across {len(cells)} axis-site pairs: {min(cells):.3f} to {max(cells):.3f} cells, {spread:.2f}x apart")
+        print("One k is only shippable if that spread is small. Two axes agreeing is time-bandwidth duality holding.")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/retina_tracker/output.py b/retina_tracker/output.py
index 662ce33..1e0438c 100644
--- a/retina_tracker/output.py
+++ b/retina_tracker/output.py
@@ -172,6 +172,11 @@ def write_residual(self, track_id, timestamp, track, detection):
             # events file, which a track that never confirmed is absent from.
             "delay": detection.get("delay"),
             "doppler": detection.get("doppler"),
+            # The occupancy map's verdict on that measurement, so calibration
+            # can drop the interference exactly rather than inferring it from a
+            # track's Doppler holding still. Present whether or not the node is
+            # suppressing on it.
+            "interfering": bool(detection.get("interfering")),
             "n_missed": track.last_n_missed,
             "q_scale": track.last_q_scale,
             "innovation": [float(x) for x in residual.innovation],
diff --git a/tests/test_calibrate.py b/tests/test_calibrate.py
new file mode 100644
index 0000000..a672803
--- /dev/null
+++ b/tests/test_calibrate.py
@@ -0,0 +1,211 @@
+"""Fitting R and Q back out of recorded innovations.
+
+The two cancel: R is around 45x too large in delay while the model error it
+stands in for grows as the square of the prediction interval, not the fifth
+power white jerk implies, so the sum looks plausible and NIS alone cannot tell
+which is wrong. What separates them is that only one of the two grows with how
+long the filter had been predicting.
+
+These tests synthesise records from a known R and a known velocity error and
+assert the fit recovers both, which is the only way to know the tool is
+measuring rather than asserting before any of it is pointed at a node.
+"""
+
+import json
+import random
+
+import pytest
+
+from retina_tracker.calibrate import (
+    calibrate,
+    fit_noise_model,
+    group_by_interval,
+    load_blah2_capture,
+    prediction_interval,
+    usable,
+    wavelength_km,
+)
+from retina_tracker.output import InnovationWriter
+from retina_tracker.tracker import Tracker
+
+FC = 177000000.0
+FS = 2000000.0
+CPI = 0.5
+
+TRUE_DELAY_SIGMA_KM = 0.022
+TRUE_DOPPLER_SIGMA_HZ = 0.356
+TRUE_VELOCITY_ERROR = 0.015
+
+INTERVALS = (0.5, 1.0, 1.5, 2.0, 3.0)
+PER_INTERVAL = 1200
+
+
+def record(interval, delay_innovation, rate_innovation, doppler=-120.0, birth=0, **overrides):
+    record = {
+        "track_id": "250915-000001",
+        "birth": birth,
+        "timestamp": 1718747745000,
+        "dt": 0.5,
+        "n_missed": round(interval / 0.5) - 1,
+        "snr": 16.0,
+        "delay": 20.0,
+        "doppler": doppler,
+        "interfering": False,
+        "q_scale": 1.0,
+        "innovation": [delay_innovation, rate_innovation],
+        "s_diag": [1.0, 5e-6],
+        "nis": 1.0,
+    }
+    record.update(overrides)
+    return record
+
+
+def synthesise(seed=1, delay_sigma=TRUE_DELAY_SIGMA_KM, velocity_error=TRUE_VELOCITY_ERROR, exponent=2.0):
+    """Innovations from a filter whose R is `delay_sigma` and whose model error
+    is a velocity held for the whole prediction interval."""
+    rng = random.Random(seed)
+    rate_sigma = TRUE_DOPPLER_SIGMA_HZ * wavelength_km(FC)
+    records = []
+    for interval in INTERVALS:
+        growth = velocity_error * interval ** (exponent / 2.0)
+        delay_spread = (delay_sigma**2 + growth**2) ** 0.5
+        rate_spread = (rate_sigma**2 + (growth / 50.0) ** 2) ** 0.5
+        for i in range(PER_INTERVAL):
+            records.append(
+                record(
+                    interval,
+                    rng.gauss(0.0, delay_spread),
+                    rng.gauss(0.0, rate_spread),
+                    doppler=-120.0 + i % 40,
+                )
+            )
+    return records
+
+
+class TestTheFitRecoversWhatWasPutIn:
+    def test_r_comes_back_from_the_intercept(self):
+        report = calibrate(synthesise(), FC, FS, CPI)
+        assert report["axes"]["delay"]["sigma"] == pytest.approx(TRUE_DELAY_SIGMA_KM, rel=0.2)
+
+    def test_r_comes_back_in_resolution_cells(self):
+        """The number that has to agree across sites, and the one a shipped
+        constant is expressed in."""
+        report = calibrate(synthesise(), FC, FS, CPI)
+        assert report["axes"]["delay"]["cells"] == pytest.approx(0.147, rel=0.25)
+
+    def test_the_doppler_axis_comes_back_in_hertz(self):
+        """The innovation is a range rate. Reported in Hz it is comparable with
+        the configured value and with the other axis in cells."""
+        report = calibrate(synthesise(), FC, FS, CPI)
+        assert report["axes"]["doppler"]["sigma"] == pytest.approx(TRUE_DOPPLER_SIGMA_HZ, rel=0.25)
+
+    def test_the_shape_of_the_model_error_comes_back(self):
+        """Two, not five. An exponent near five would mean white jerk really is
+        what the filter is missing and range_jerk merely needs raising."""
+        report = calibrate(synthesise(), FC, FS, CPI)
+        assert report["axes"]["delay"]["exponent"] == pytest.approx(2.0, abs=0.4)
+
+    def test_the_size_of_the_model_error_comes_back(self):
+        report = calibrate(synthesise(), FC, FS, CPI)
+        assert report["axes"]["delay"]["error_rate"] == pytest.approx(TRUE_VELOCITY_ERROR, rel=0.3)
+
+    def test_a_larger_r_is_reported_larger(self):
+        """The oversized R the tracker ships with today, against the measured
+        one, on data whose model error is identical."""
+        loose = calibrate(synthesise(seed=2, delay_sigma=0.2), FC, FS, CPI)
+        tight = calibrate(synthesise(seed=2), FC, FS, CPI)
+        assert loose["axes"]["delay"]["sigma"] > 5 * tight["axes"]["delay"]["sigma"]
+
+    def test_white_jerk_is_told_apart_from_a_velocity_error(self):
+        report = calibrate(synthesise(seed=3, exponent=5.0), FC, FS, CPI)
+        assert report["axes"]["delay"]["exponent"] > 3.0
+
+
+class TestSeparatingTheTwo:
+    def test_a_constant_scatter_reads_as_all_r_and_no_q(self):
+        report = calibrate(synthesise(seed=4, velocity_error=0.0), FC, FS, CPI)
+        assert report["axes"]["delay"]["sigma"] == pytest.approx(TRUE_DELAY_SIGMA_KM, rel=0.15)
+        assert report["axes"]["delay"]["error_rate"] < TRUE_DELAY_SIGMA_KM
+
+    def test_the_interval_is_what_carries_the_separation(self):
+        """Without a spread of prediction intervals there is one equation and
+        two unknowns, and the honest answer is to decline."""
+        records = [r for r in synthesise(seed=5) if prediction_interval(r) == 0.5]
+        assert "sigma" not in calibrate(records, FC, FS, CPI)["axes"]["delay"]
+
+    def test_a_coasted_update_is_a_longer_prediction(self):
+        assert prediction_interval({"dt": 0.5, "n_missed": 3}) == pytest.approx(2.0)
+
+    def test_scatter_is_reported_against_what_the_filter_claimed(self):
+        """Per axis, because NIS sums the two and that is how an oversized R
+        and an undersized Q have been hiding inside it."""
+        report = calibrate(synthesise(seed=6), FC, FS, CPI)
+        # Far below one is the signature being hunted: the filter is carrying
+        # an R so large that its innovations never come close to filling it.
+        assert report["axes"]["delay"]["scatter_vs_claimed"] < 0.1
+
+
+class TestWhatIsExcluded:
+    def test_interference_is_dropped_on_the_records_own_verdict(self):
+        records = synthesise(seed=7)[:200]
+        for r in records[:120]:
+            r["interfering"] = True
+        assert len(usable(records)) == 80
+
+    def test_a_degenerate_update_is_not_a_measurement(self):
+        """A singular covariance skipped the measurement, so its innovation is
+        NaN and there is nothing in it to fit."""
+        records = [record(0.5, float("nan"), float("nan")), record(0.5, 0.01, 1e-5)]
+        assert len(usable(records)) == 1
+
+    def test_older_records_fall_back_to_a_tracks_doppler_holding_still(self):
+        """Written before the tracker carried the occupancy map's verdict."""
+        tone = [record(0.5, 0.01, 1e-5, doppler=27.9, birth=1) for _ in range(12)]
+        aircraft = [record(0.5, 0.01, 1e-5, doppler=-120.0 + 4 * i, birth=2) for i in range(12)]
+        for r in tone + aircraft:
+            del r["interfering"]
+        kept = usable(tone + aircraft)
+        assert {r["birth"] for r in kept} == {2}
+
+    def test_the_fallback_is_not_used_when_a_verdict_is_present(self):
+        """A real target can hold one Doppler bin. Where the map has spoken,
+        the heuristic that would convict it must not run."""
+        holding = [record(0.5, 0.01, 1e-5, doppler=27.9, birth=1) for _ in range(12)]
+        assert len(usable(holding)) == 12
+
+    def test_a_thin_interval_is_not_fitted(self):
+        assert fit_noise_model([(0.5, 40, 1e-4), (1.0, 40, 2e-4)]) is None
+
+    def test_a_sparse_interval_is_not_a_measurement(self):
+        records = [record(0.5, 0.01, 1e-5) for _ in range(5)]
+        assert group_by_interval(records, 0) == []
+
+
+class TestTheCaptureConfig:
+    def test_all_three_parameters_come_out_of_blah2s_own_file(self, tmp_path):
+        path = tmp_path / "config.yml"
+        path.write_text(
+            "capture:\n  fs: 2000000\n  fc: 177000000\n  device:\n    type: 'RspDuo'\n"
+            "process:\n  data:\n    cpi: 0.5\n    buffer: 1.5\n"
+        )
+        assert load_blah2_capture(str(path)) == {"fs": 2000000.0, "fc": 177000000.0, "cpi": 0.5}
+
+
+class TestItReadsWhatTheTrackerWrites:
+    """The contract between the two is a set of field names in a file. A test
+    that synthesises its own records would never notice one being renamed."""
+
+    def test_a_real_recording_calibrates(self, tmp_path):
+        path = tmp_path / "innovations.jsonl"
+        writer = InnovationWriter(str(path), max_bytes=0)
+        tracker = Tracker(innovation_writer=writer)
+        for i in range(60):
+            tracker.process_frame([{"delay": 20.0 - 0.05 * i, "doppler": -120.0, "snr": 16.0}], 1718747745000 + i * 500)
+        writer.close()
+
+        records = [json.loads(line) for line in path.read_text().splitlines() if line.strip()]
+        report = calibrate(records, FC, FS, CPI)
+
+        assert report["n_records"] == len(records)
+        assert report["n_used"] > 0
+        assert report["nis_median"] is not None
diff --git a/tests/test_innovations.py b/tests/test_innovations.py
index cd7780a..4efc228 100644
--- a/tests/test_innovations.py
+++ b/tests/test_innovations.py
@@ -93,6 +93,7 @@ def test_a_record_carries_what_calibration_needs(self, tmp_path):
             "snr",
             "delay",
             "doppler",
+            "interfering",
             "n_missed",
             "q_scale",
             "innovation",