#!/usr/bin/env python3
"""Quantify timestamp-gap distortion in the public GL-MT6000 CSI capture.

The capture container is:
    <float64 timestamp> <uint16 payload-length> <payload>

The 24-byte CSI payload header is:
    <I timestamp_us, 6s transmitter, b RSSI, B SNR, B bandwidth,
     B channel, B mode, B epoch, H subcarriers, H tx, H rx, H chain_info>

This script deliberately depends only on NumPy and Matplotlib. It compares:
1. linear interpolation onto a uniform 10 Hz grid followed by Welch PSD;
2. a least-squares sinusoid fit on the original uneven timestamps, the
   Lomb-Scargle family of periodograms;
3. a known 0.25 Hz synthetic control placed on the real sampling times.

The capture has no human-event ground truth. The actual-data comparison shows
that the two estimators disagree; it does not identify a physiological cause.
"""

from __future__ import annotations

import argparse
import collections
import csv
import hashlib
import html
import json
import math
import struct
from pathlib import Path

import numpy as np


RECORD = struct.Struct("<dH")
HEADER = struct.Struct("<I6sbBBBBBHHHH")
GAP_THRESHOLD_S = 0.5
UNIFORM_RATE_HZ = 10.0
LS_STEP_HZ = 0.001
SYNTHETIC_FREQUENCIES_HZ = [0.10, 0.15, 0.20, 0.25, 0.30, 0.40, 0.50]


def rounded(value: float, digits: int = 9) -> float:
    return round(float(value), digits)


def read_groups(path: Path):
    data = path.read_bytes()
    offset = 0
    pending = []
    groups = []
    outer_records = 0
    malformed_records = 0

    while offset + RECORD.size <= len(data):
        timestamp, length = RECORD.unpack_from(data, offset)
        payload_start = offset + RECORD.size
        payload_end = payload_start + length
        if payload_end > len(data) or length < HEADER.size:
            malformed_records += 1
            break

        header = HEADER.unpack_from(data, payload_start)
        _, _, rssi, snr, bandwidth, channel, mode, epoch, count, tx, rx, chain = header
        iq_offset = payload_start + HEADER.size
        available_pairs = max(0, (payload_end - iq_offset) // 4)
        pairs = min(count, available_pairs)
        iq = np.frombuffer(
            data,
            dtype="<i2",
            count=pairs * 2,
            offset=iq_offset,
        )
        pending.append(
            {
                "outer_ts": timestamp,
                "rssi": rssi,
                "snr": snr,
                "bandwidth": bandwidth,
                "channel": channel,
                "mode": mode,
                "epoch": epoch,
                "subcarriers": count,
                "tx": tx,
                "rx": rx,
                "chain": chain,
                "iq": iq,
                "nonzero_iq": bool(iq.size and np.any(iq)),
            }
        )
        outer_records += 1
        if chain & 0x8000:
            groups.append(pending)
            pending = []
        offset = payload_end

    return {
        "data": data,
        "groups": groups,
        "outer_records": outer_records,
        "incomplete_tail_records": len(pending),
        "malformed_records": malformed_records,
    }


def amplitude(iq: np.ndarray) -> np.ndarray:
    return np.hypot(iq[0::2].astype(np.float64), iq[1::2].astype(np.float64))


def build_series(groups):
    times = []
    rows = []
    chain_counts = []
    selected_chains = collections.Counter()

    for group in groups:
        valid = [
            record
            for record in group
            if record["subcarriers"] == 256 and record["nonzero_iq"]
        ]
        if not valid:
            continue

        chosen = next(
            (
                record
                for record in valid
                if (record["tx"], record["rx"]) == (0, 0)
            ),
            None,
        )
        if chosen is None:
            chosen = max(valid, key=lambda record: amplitude(record["iq"]).mean())

        times.append(chosen["outer_ts"])
        rows.append(amplitude(chosen["iq"]))
        chain_counts.append(len(group))
        selected_chains[(int(chosen["tx"]), int(chosen["rx"]))] += 1

    if len(times) < 3:
        raise SystemExit("not enough usable measurement groups")

    times = np.asarray(times, dtype=np.float64)
    times -= times[0]
    matrix = np.vstack(rows)
    return {
        "times_s": times,
        "feature": matrix.mean(axis=1),
        "matrix_shape": list(matrix.shape),
        "chain_counts": np.asarray(chain_counts),
        "selected_chains": {
            f"tx{tx}/rx{rx}": count
            for (tx, rx), count in sorted(selected_chains.items())
        },
    }


def segment_durations(times: np.ndarray, intervals: np.ndarray, threshold: float):
    edges = np.r_[0, np.where(intervals > threshold)[0] + 1, len(times)]
    durations = times[edges[1:] - 1] - times[edges[:-1]]
    return durations, edges


def welch_psd(values: np.ndarray, rate_hz: float, nperseg: int = 1024):
    nperseg = max(16, min(nperseg, len(values)))
    step = max(1, nperseg // 2)
    window = np.hanning(nperseg)
    scale = rate_hz * np.sum(window**2)
    spectra = []

    for start in range(0, len(values) - nperseg + 1, step):
        segment = values[start : start + nperseg].astype(np.float64, copy=True)
        segment -= segment.mean()
        segment *= window
        spectrum = np.fft.rfft(segment)
        power = np.abs(spectrum) ** 2 / scale
        if power.size > 2:
            power[1:-1] *= 2.0
        spectra.append(power)

    if not spectra:
        raise ValueError("signal is shorter than one Welch window")
    return np.fft.rfftfreq(nperseg, 1.0 / rate_hz), np.asarray(spectra).mean(axis=0)


def irregular_least_squares_periodogram(
    times: np.ndarray,
    values: np.ndarray,
    frequencies: np.ndarray,
):
    centered = values.astype(np.float64, copy=True)
    centered -= centered.mean()
    denominator = float(np.dot(centered, centered))
    if denominator <= 0:
        raise ValueError("signal has zero variance")

    power = np.empty_like(frequencies)
    for index, frequency in enumerate(frequencies):
        phase = 2.0 * np.pi * frequency * times
        sine = np.sin(phase)
        cosine = np.cos(phase)
        ss = float(np.dot(sine, sine))
        cc = float(np.dot(cosine, cosine))
        sc = float(np.dot(sine, cosine))
        ys = float(np.dot(centered, sine))
        yc = float(np.dot(centered, cosine))
        determinant = ss * cc - sc * sc
        if abs(determinant) < 1e-18:
            power[index] = 0.0
            continue
        sine_coefficient = (ys * cc - yc * sc) / determinant
        cosine_coefficient = (yc * ss - ys * sc) / determinant
        explained = (
            sine_coefficient * ys + cosine_coefficient * yc
        ) / denominator
        power[index] = max(0.0, explained)
    return power


def band_power(frequencies, power, low, high):
    mask = (frequencies >= low) & (frequencies < high)
    if not np.any(mask):
        return 0.0
    integrate = getattr(np, "trapezoid", np.trapz)
    return float(integrate(power[mask], frequencies[mask]))


def peak_in_band(frequencies, power, low=0.05, high=0.60):
    mask = (frequencies >= low) & (frequencies <= high)
    index = int(np.argmax(power[mask]))
    return float(frequencies[mask][index]), float(power[mask][index])


def interpolation_mask(
    times: np.ndarray,
    intervals: np.ndarray,
    uniform_times: np.ndarray,
    threshold: float,
):
    left = np.searchsorted(times, uniform_times, side="right") - 1
    valid = (left >= 0) & (left < len(times) - 1)
    mask = np.zeros(len(uniform_times), dtype=bool)
    mask[valid] = intervals[left[valid]] > threshold
    return mask


def synthetic_metrics(
    times: np.ndarray,
    uniform_times: np.ndarray,
    frequency_hz: float,
):
    samples = np.sin(2.0 * np.pi * frequency_hz * times)
    interpolated = np.interp(uniform_times, times, samples)
    truth = np.sin(2.0 * np.pi * frequency_hz * uniform_times)
    error = interpolated - truth

    exact_f, exact_power = welch_psd(truth, UNIFORM_RATE_HZ)
    interp_f, interp_power = welch_psd(interpolated, UNIFORM_RATE_HZ)
    ls_frequencies = np.arange(0.05, 0.600 + LS_STEP_HZ / 2, LS_STEP_HZ)
    ls_power = irregular_least_squares_periodogram(
        times, samples, ls_frequencies
    )

    exact_peak, _ = peak_in_band(exact_f, exact_power)
    interp_peak, _ = peak_in_band(interp_f, interp_power)
    ls_peak, _ = peak_in_band(ls_frequencies, ls_power)
    signal_rms = float(np.sqrt(np.mean(truth**2)))
    error_rms = float(np.sqrt(np.mean(error**2)))

    exact_band = band_power(exact_f, exact_power, 0.10, 0.50)
    interp_band = band_power(interp_f, interp_power, 0.10, 0.50)
    exact_total = band_power(exact_f, exact_power, 0.05, 1.0)
    interp_total = band_power(interp_f, interp_power, 0.05, 1.0)

    return {
        "frequency_hz": frequency_hz,
        "samples": samples,
        "interpolated": interpolated,
        "truth": truth,
        "exact_frequencies_hz": exact_f,
        "exact_power": exact_power,
        "interpolated_frequencies_hz": interp_f,
        "interpolated_power": interp_power,
        "ls_frequencies_hz": ls_frequencies,
        "ls_power": ls_power,
        "metrics": {
            "true_frequency_hz": frequency_hz,
            "exact_uniform_peak_hz": rounded(exact_peak, 6),
            "interpolated_welch_peak_hz": rounded(interp_peak, 6),
            "ls_peak_hz": rounded(ls_peak, 6),
            "interpolation_error_rms": rounded(error_rms),
            "interpolation_error_max_abs": rounded(np.max(np.abs(error))),
            "signal_rms": rounded(signal_rms),
            "relative_rms_error": rounded(error_rms / signal_rms, 6),
            "respiratory_band_power_ratio_interpolated_to_exact": rounded(
                interp_band / exact_band if exact_band else 0.0, 6
            ),
            "in_band_power_fraction_exact": rounded(
                exact_band / exact_total if exact_total else 0.0, 6
            ),
            "in_band_power_fraction_interpolated": rounded(
                interp_band / interp_total if interp_total else 0.0, 6
            ),
        },
    }


def normalized(values):
    maximum = float(np.max(values))
    return values / maximum if maximum > 0 else values


def downsample(x, y, limit=1800):
    if len(x) <= limit:
        return x, y
    index = np.linspace(0, len(x) - 1, limit).astype(int)
    return x[index], y[index]


def svg_line(x1, y1, x2, y2, color="#334155", width=1, dash=None, opacity=1):
    dash_attr = f' stroke-dasharray="{dash}"' if dash else ""
    return (
        f'<line x1="{x1:.2f}" y1="{y1:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" '
        f'stroke="{color}" stroke-width="{width}" opacity="{opacity}"{dash_attr}/>'
    )


def svg_text(x, y, value, size=12, color="#e5e7eb", anchor="start", weight=400):
    return (
        f'<text x="{x:.2f}" y="{y:.2f}" fill="{color}" font-size="{size}" '
        f'font-family="ui-sans-serif,system-ui,sans-serif" '
        f'text-anchor="{anchor}" font-weight="{weight}">{html.escape(str(value))}</text>'
    )


def svg_polyline(
    x_values,
    y_values,
    box,
    x_range,
    y_range,
    color,
    width=1.5,
    opacity=1,
):
    left, top, width_px, height_px = box
    x_min, x_max = x_range
    y_min, y_max = y_range
    x_values, y_values = downsample(
        np.asarray(x_values, dtype=np.float64),
        np.asarray(y_values, dtype=np.float64),
    )
    x_span = x_max - x_min or 1.0
    y_span = y_max - y_min or 1.0
    points = []
    for x_value, y_value in zip(x_values, y_values):
        if not np.isfinite(y_value):
            continue
        px = left + (float(x_value) - x_min) / x_span * width_px
        py = top + height_px - (float(y_value) - y_min) / y_span * height_px
        points.append(f"{px:.2f},{py:.2f}")
    return (
        f'<polyline points="{" ".join(points)}" fill="none" stroke="{color}" '
        f'stroke-width="{width}" stroke-linejoin="round" '
        f'stroke-linecap="round" opacity="{opacity}"/>'
    )


def panel_frame(title, box, x_label, y_label):
    left, top, width, height = box
    elements = [
        f'<rect x="{left}" y="{top}" width="{width}" height="{height}" rx="8" fill="#111827" stroke="#334155"/>',
        svg_text(left + width / 2, top - 14, title, 15, anchor="middle", weight=600),
        svg_text(left + width / 2, top + height + 42, x_label, 12, "#cbd5e1", "middle"),
        f'<g transform="translate({left - 14},{top + height / 2}) rotate(-90)">'
        f'{svg_text(0, 0, y_label, 12, "#cbd5e1", "middle")}</g>',
    ]
    return elements


def legend(items, x, y):
    width = max(len(label) for _, label in items) * 6.6 + 55
    height = len(items) * 20 + 10
    elements = [
        f'<rect x="{x - 10:.2f}" y="{y - 18:.2f}" width="{width:.2f}" height="{height:.2f}" '
        'rx="5" fill="#111827" fill-opacity="0.9" stroke="#334155"/>'
    ]
    for index, (color, label) in enumerate(items):
        row_y = y + index * 20
        elements.append(svg_line(x, row_y - 4, x + 26, row_y - 4, color, 3))
        elements.append(svg_text(x + 34, row_y, label, 11, "#cbd5e1"))
    return elements


def tick_grid(
    box,
    x_ticks,
    y_ticks,
    x_range,
    y_range,
    x_formatter=lambda value: f"{value:g}",
    y_formatter=lambda value: f"{value:g}",
):
    left, top, width, height = box
    x_min, x_max = x_range
    y_min, y_max = y_range
    x_span = x_max - x_min or 1.0
    y_span = y_max - y_min or 1.0
    elements = []
    for value in x_ticks:
        x = left + (value - x_min) / x_span * width
        elements.append(svg_line(x, top, x, top + height, "#1f2937", 1))
        elements.append(svg_text(x, top + height + 17, x_formatter(value), 10, "#94a3b8", "middle"))
    for value in y_ticks:
        y = top + height - (value - y_min) / y_span * height
        elements.append(svg_line(left, y, left + width, y, "#1f2937", 1))
        elements.append(svg_text(left - 8, y + 3, y_formatter(value), 10, "#94a3b8", "end"))
    return elements


def plot_results(
    times,
    intervals,
    synthetic,
    actual,
    output: Path,
):
    width, height = 1400, 920
    panels = [
        (90, 115, 580, 285),
        (790, 115, 540, 285),
        (90, 505, 580, 285),
        (790, 505, 540, 285),
    ]
    elements = [
        f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {width} {height}" role="img" aria-labelledby="title desc">',
        '<title id="title">GL-MT6000 CSI timestamp gaps and spectral-method sensitivity</title>',
        '<desc id="desc">Four panels compare real CSI timing gaps, a known 0.25 Hz synthetic control, its spectra, and the two estimators on the real CSI feature.</desc>',
        f'<rect width="{width}" height="{height}" fill="#0b0f17"/>',
        svg_text(
            width / 2,
            38,
            "GL-MT6000 CSI timestamp gaps and spectral-method sensitivity",
            24,
            anchor="middle",
            weight=700,
        ),
        svg_text(
            width / 2,
            63,
            "Real timestamps · synthetic 0.25 Hz control · interpolated Welch vs irregular least squares",
            13,
            "#9aa7b8",
            "middle",
        ),
    ]

    elements.extend(panel_frame("Real CSI group intervals", panels[0], "Elapsed time (s)", "Interval (s, log10)"))
    left, top, panel_width, panel_height = panels[0]
    elements.extend(
        tick_grid(
            panels[0],
            [0, 150, 300, 450, 600],
            [-3, -2, -1, 0],
            (0, float(times[-1])),
            (-3.5, 0.5),
            y_formatter=lambda value: f"{10 ** int(value):g}",
        )
    )
    interval_x = times[1:]
    interval_y = np.log10(np.maximum(intervals, 1e-6))
    for interval, color, dash, label in [
        (0.5, "#fb923c", "8 5", "0.5 s"),
        (1.0, "#f87171", "3 5", "1.0 s"),
    ]:
        y = top + panel_height - (math.log10(interval) - (-3.5)) / 4.0 * panel_height
        elements.append(svg_line(left, y, left + panel_width, y, color, 1.2, dash))
        elements.append(svg_text(left + panel_width - 8, y - 6, label, 11, color, "end"))
    elements.append(
        svg_polyline(
            interval_x,
            interval_y,
            panels[0],
            (0, float(times[-1])),
            (-3.5, 0.5),
            "#38bdf8",
            1.2,
            0.9,
        )
    )

    largest_gap_index = int(np.argmax(intervals))
    gap_center = float((times[largest_gap_index] + times[largest_gap_index + 1]) / 2.0)
    plot_start = max(0.0, gap_center - 8.0)
    plot_end = min(float(times[-1]), gap_center + 8.0)
    elements.extend(
        panel_frame(
            f"Synthetic control around the largest {intervals[largest_gap_index]:.2f} s gap",
            panels[1],
            "Elapsed time (s)",
            "Amplitude",
        )
    )
    mask = (
        (synthetic["truth_times_s"] >= plot_start)
        & (synthetic["truth_times_s"] <= plot_end)
    )
    left, top, panel_width, panel_height = panels[1]
    elements.extend(
        tick_grid(
            panels[1],
            np.arange(math.ceil(plot_start / 4) * 4, plot_end + 0.1, 4),
            [-1, 0, 1],
            (plot_start, plot_end),
            (-1.5, 1.5),
            x_formatter=lambda value: f"{value:g}",
        )
    )
    elements.append(
        svg_polyline(
            synthetic["truth_times_s"][mask],
            synthetic["truth"][mask],
            panels[1],
            (plot_start, plot_end),
            (-1.5, 1.5),
            "#a3e635",
            1.8,
        )
    )
    elements.append(
        svg_polyline(
            synthetic["truth_times_s"][mask],
            synthetic["interpolated"][mask],
            panels[1],
            (plot_start, plot_end),
            (-1.5, 1.5),
            "#fb923c",
            1.3,
            0.9,
        )
    )
    elements.extend(
        legend(
            [
                ("#a3e635", "Known 0.25 Hz"),
                ("#fb923c", "Linear interpolation"),
            ],
            left + panel_width - 190,
            top + 24,
        )
    )

    elements.extend(panel_frame("Known 0.25 Hz control spectrum", panels[2], "Frequency (Hz)", "Normalized power"))
    left, top, panel_width, panel_height = panels[2]
    elements.extend(
        tick_grid(
            panels[2],
            [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
            [0, 0.5, 1],
            (0, 0.60),
            (0, 1.05),
            x_formatter=lambda value: f"{value:.1f}",
            y_formatter=lambda value: f"{value:.1f}",
        )
    )
    spectrum_mask = synthetic["exact_frequencies_hz"] <= 0.60
    elements.append(
        svg_polyline(
            synthetic["exact_frequencies_hz"][spectrum_mask],
            normalized(synthetic["exact_power"])[spectrum_mask],
            panels[2],
            (0, 0.60),
            (0, 1.05),
            "#a3e635",
            2.0,
        )
    )
    elements.append(
        svg_polyline(
            synthetic["interpolated_frequencies_hz"][spectrum_mask],
            normalized(synthetic["interpolated_power"])[spectrum_mask],
            panels[2],
            (0, 0.60),
            (0, 1.05),
            "#fb923c",
            1.5,
            0.9,
        )
    )
    ls_mask = synthetic["ls_frequencies_hz"] <= 0.60
    elements.append(
        svg_polyline(
            synthetic["ls_frequencies_hz"][ls_mask],
            normalized(synthetic["ls_power"])[ls_mask],
            panels[2],
            (0, 0.60),
            (0, 1.05),
            "#c084fc",
            1.5,
        )
    )
    elements.extend(
        legend(
            [
                ("#a3e635", "Known uniform"),
                ("#fb923c", "Interpolated Welch"),
                ("#c084fc", "Irregular least squares"),
            ],
            left + panel_width - 220,
            top + 24,
        )
    )

    elements.extend(panel_frame("Real CSI feature: estimator comparison", panels[3], "Frequency (Hz)", "Normalized power"))
    left, top, panel_width, panel_height = panels[3]
    elements.extend(
        tick_grid(
            panels[3],
            [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
            [0, 0.5, 1],
            (0, 0.60),
            (0, 1.05),
            x_formatter=lambda value: f"{value:.1f}",
            y_formatter=lambda value: f"{value:.1f}",
        )
    )
    actual_naive_mask = actual["naive_frequencies_hz"] <= 0.60
    elements.append(
        svg_polyline(
            actual["naive_frequencies_hz"][actual_naive_mask],
            normalized(actual["naive_power"])[actual_naive_mask],
            panels[3],
            (0, 0.60),
            (0, 1.05),
            "#fb923c",
            1.6,
        )
    )
    actual_ls_mask = actual["ls_frequencies_hz"] <= 0.60
    elements.append(
        svg_polyline(
            actual["ls_frequencies_hz"][actual_ls_mask],
            normalized(actual["ls_power"])[actual_ls_mask],
            panels[3],
            (0, 0.60),
            (0, 1.05),
            "#38bdf8",
            1.6,
        )
    )
    elements.extend(
        legend(
            [
                ("#fb923c", "Interpolated Welch"),
                ("#38bdf8", "Irregular least squares"),
            ],
            left + panel_width - 230,
            top + 24,
        )
    )

    elements.append(
        svg_text(
            width / 2,
            height - 18,
            "Source: 10-minute public capture, SHA-256 527466ce…; no human-event ground truth.",
            11,
            "#9aa7b8",
            "middle",
        )
    )
    elements.append("</svg>")
    output.write_text("\n".join(elements) + "\n", encoding="utf-8")


def main() -> int:
    parser = argparse.ArgumentParser()
    parser.add_argument("capture", type=Path)
    parser.add_argument("--summary", type=Path, required=True)
    parser.add_argument("--img", type=Path, required=True)
    parser.add_argument("--csv", type=Path, required=True)
    args = parser.parse_args()

    parsed = read_groups(args.capture)
    series = build_series(parsed["groups"])
    times = series["times_s"]
    feature = series["feature"]
    intervals = np.diff(times)
    duration = float(times[-1])

    uniform_times = np.arange(0.0, duration + 1e-12, 1.0 / UNIFORM_RATE_HZ)
    mask = interpolation_mask(
        times, intervals, uniform_times, GAP_THRESHOLD_S
    )

    normalized_feature = (feature - feature.mean()) / feature.std()
    naive_feature = np.interp(uniform_times, times, normalized_feature)
    naive_frequencies, naive_power = welch_psd(
        naive_feature, UNIFORM_RATE_HZ
    )
    ls_frequencies = np.arange(0.05, 0.600 + LS_STEP_HZ / 2, LS_STEP_HZ)
    ls_power = irregular_least_squares_periodogram(
        times, normalized_feature, ls_frequencies
    )
    naive_peak, _ = peak_in_band(naive_frequencies, naive_power)
    ls_peak, _ = peak_in_band(ls_frequencies, ls_power)

    synthetic = synthetic_metrics(times, uniform_times, 0.25)
    synthetic["truth_times_s"] = uniform_times

    sweep_rows = []
    for frequency in SYNTHETIC_FREQUENCIES_HZ:
        result = synthetic_metrics(times, uniform_times, frequency)
        sweep_rows.append(result["metrics"])

    durations_half, _ = segment_durations(times, intervals, 0.5)
    durations_one, _ = segment_durations(times, intervals, 1.0)

    summary = {
        "source": {
            "file": args.capture.name,
            "size_bytes": len(parsed["data"]),
            "sha256": hashlib.sha256(parsed["data"]).hexdigest(),
        },
        "capture": {
            "outer_records": parsed["outer_records"],
            "closed_groups": len(parsed["groups"]),
            "usable_groups": len(times),
            "incomplete_tail_records": parsed["incomplete_tail_records"],
            "malformed_records": parsed["malformed_records"],
            "duration_s": rounded(duration),
            "chain_records_per_second": rounded(
                parsed["outer_records"] / duration
            ),
            "usable_groups_per_second": rounded(len(times) / duration),
            "matrix_shape": series["matrix_shape"],
            "selected_chains": series["selected_chains"],
        },
        "timing": {
            "median_interval_s": rounded(np.median(intervals)),
            "p95_interval_s": rounded(np.percentile(intervals, 95)),
            "p99_interval_s": rounded(np.percentile(intervals, 99)),
            "max_interval_s": rounded(np.max(intervals)),
            "gap_threshold_s": GAP_THRESHOLD_S,
            "long_gaps_over_0_5_s": int(np.sum(intervals > 0.5)),
            "long_gaps_over_1_s": int(np.sum(intervals > 1.0)),
            "segments_ge_20_s_at_0_5_threshold": int(
                np.sum(durations_half >= 20.0)
            ),
            "coverage_s_ge_20_s_at_0_5_threshold": rounded(
                np.sum(durations_half[durations_half >= 20.0])
            ),
            "segments_ge_20_s_at_1_threshold": int(
                np.sum(durations_one >= 20.0)
            ),
            "coverage_s_ge_20_s_at_1_threshold": rounded(
                np.sum(durations_one[durations_one >= 20.0])
            ),
        },
        "resampling": {
            "uniform_rate_hz": UNIFORM_RATE_HZ,
            "uniform_samples": int(len(uniform_times)),
            "samples_inside_over_0_5_s_intervals": int(mask.sum()),
            "fraction_samples_inside_over_0_5_s_intervals": rounded(
                mask.mean(), 6
            ),
        },
        "synthetic": synthetic["metrics"],
        "actual": {
            "feature": "mean amplitude across the selected chain's 256 subcarriers",
            "naive_welch_peak_0_05_to_0_60_hz": rounded(naive_peak, 6),
            "irregular_ls_peak_0_05_to_0_60_hz": rounded(ls_peak, 6),
            "peak_difference_hz": rounded(ls_peak - naive_peak, 6),
            "naive_respiratory_band_power_0_10_to_0_50": rounded(
                band_power(naive_frequencies, naive_power, 0.10, 0.50)
            ),
            "naive_mid_band_power_0_50_to_1_00": rounded(
                band_power(naive_frequencies, naive_power, 0.50, 1.0)
            ),
            "naive_motion_band_power_1_00_to_4_00": rounded(
                band_power(naive_frequencies, naive_power, 1.0, 4.0)
            ),
        },
        "method": {
            "record_layout": "float64 outer timestamp + uint16 length + 24-byte CSI header + int16 I/Q",
            "group_boundary": "chain_info bit 15",
            "chain_selection": "prefer tx0/rx0; otherwise the valid chain with highest mean amplitude",
            "uniform_grid": "linear interpolation at 10 Hz; Hann Welch windows of 1024 samples with 50% overlap",
            "irregular_periodogram": "least-squares sinusoid fit on original timestamps over 0.05-0.60 Hz at 0.001 Hz spacing",
            "synthetic_control": "unit-amplitude sinusoids evaluated at the real timestamps, then compared on the real 10 Hz grid",
        },
        "limitations": [
            "The capture has no human-event or respiration ground truth.",
            "The synthetic control tests estimator behaviour; it is not a breathing recording.",
            "A single indoor capture cannot establish general behaviour across rooms, clients or traffic.",
            "The 10 Hz grid is suitable for the 0.05-4 Hz analysis here, but it is not a claim that the original stream was uniformly sampled.",
            "Agreement or disagreement between estimators does not by itself identify a person, intent or medical state.",
        ],
    }

    args.summary.parent.mkdir(parents=True, exist_ok=True)
    args.summary.write_text(
        json.dumps(summary, indent=2, ensure_ascii=True) + "\n",
        encoding="utf-8",
    )

    args.csv.parent.mkdir(parents=True, exist_ok=True)
    with args.csv.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(sweep_rows[0]))
        writer.writeheader()
        writer.writerows(sweep_rows)

    actual = {
        "naive_frequencies_hz": naive_frequencies,
        "naive_power": naive_power,
        "ls_frequencies_hz": ls_frequencies,
        "ls_power": ls_power,
    }
    args.img.parent.mkdir(parents=True, exist_ok=True)
    plot_results(
        times,
        intervals,
        synthetic,
        actual,
        args.img,
    )
    print(
        json.dumps(
            {
                "usable_groups": len(times),
                "duration_s": rounded(duration, 3),
                "synthetic_interpolation_error_rms": summary["synthetic"][
                    "interpolation_error_rms"
                ],
                "actual_naive_peak_hz": summary["actual"][
                    "naive_welch_peak_0_05_to_0_60_hz"
                ],
                "actual_ls_peak_hz": summary["actual"][
                    "irregular_ls_peak_0_05_to_0_60_hz"
                ],
            },
            indent=2,
        )
    )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())