diff --git a/docs/api/data.md b/docs/api/data.md
index c66a43f..8d746f1 100644
--- a/docs/api/data.md
+++ b/docs/api/data.md
@@ -10,3 +10,7 @@
 ::: harp.data.to_file
 ::: harp.data.to_buffer
 ::: harp.data.REFERENCE_EPOCH
+::: harp.data.synchronization.decode_clock_from_samples
+::: harp.data.synchronization.decode_clock_from_transitions
+::: harp.data.synchronization.ClockAnchor
+::: harp.data.synchronization.DEFAULT_BAUD_RATE
diff --git a/docs/examples/align_to_harp_clock/align_to_harp_clock.md b/docs/examples/align_to_harp_clock/align_to_harp_clock.md
new file mode 100644
index 0000000..5eaf6ed
--- /dev/null
+++ b/docs/examples/align_to_harp_clock/align_to_harp_clock.md
@@ -0,0 +1,23 @@
+# Aligning Local Timestamps to the Harp Clock
+
+Devices that are not part of the Harp bus run on their own clock, so their
+timestamps drift with respect to Harp time. Some Harp clock emitters mirror the
+[Synchronization Clock](https://harp-tech.org/protocol/SynchronizationClock.html)
+on a digital output at a much lower baud rate — typically 1 kbps instead of
+100 kbps — precisely so that such devices can record it on a spare digital (or
+analog) input and be aligned afterwards.
+
+This example decodes that recording back into Harp seconds, which local timestamps
+can then be expressed against.
+
+!!! note
+    The decoded table is a set of anchors: local time → whole Harp second. How
+    timestamps are placed between them is up to you — interpolating between
+    neighbouring anchors absorbs the drift between the two clocks, whereas a global
+    fit trades that away for noise rejection.
+
+<!--codeinclude-->
+```python
+[](./align_to_harp_clock.py)
+```
+<!--/codeinclude-->
diff --git a/docs/examples/align_to_harp_clock/align_to_harp_clock.py b/docs/examples/align_to_harp_clock/align_to_harp_clock.py
new file mode 100644
index 0000000..a982285
--- /dev/null
+++ b/docs/examples/align_to_harp_clock/align_to_harp_clock.py
@@ -0,0 +1,28 @@
+import numpy as np
+from harp.data.synchronization import decode_clock_from_samples, decode_clock_from_transitions
+
+# A non-Harp acquisition system recorded the downsampled Harp clock on one of its
+# digital lines. `samples` is that line, sampled at the system's own rate.
+samples = np.load("sync_line.npy")  # digital states; analog input needs a `threshold`
+sample_rate = 30_000.0
+
+# Decode it: one row per whole Harp second, keyed on the sample the packet was
+# anchored on — the axis this system timestamps the rest of its data on too.
+clock = decode_clock_from_samples(samples, sample_rate, baud_rate=1000.0)
+print(clock.head())
+#              Time
+# Sample
+# 37200   3806874.0
+# 67203   3806875.0
+
+# Anchors, so any of the system's timestamps — spikes, video frames, stimulus onsets —
+# can be placed on the Harp axis. Interpolating between neighbouring anchors absorbs
+# the drift between the two clocks.
+spike_samples = np.load("spike_samples.npy")
+harp_times = np.interp(spike_samples, clock.index, clock["Time"])
+
+# Event-based systems report line transitions instead of a sampled waveform: a local
+# time and the level the line took. Anchors then carry local seconds.
+transitions = np.load("line_transitions.npy")
+clock = decode_clock_from_transitions(transitions[:, 0], transitions[:, 1])
+harp_times = np.interp(spike_samples / sample_rate, clock.index, clock["Time"])
diff --git a/mkdocs.yml b/mkdocs.yml
index 23d980e..9c3954b 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -79,6 +79,7 @@ nav:
     - Subscribing to Events: examples/subscribing_to_events/subscribing_to_events.md
     - Reading a Whole Dataset Folder: examples/read_dataset/read_dataset.md
     - Reading Data into a DataFrame: examples/read_data_to_dataframe/read_data_to_dataframe.md
+    - Aligning Local Timestamps to the Harp Clock: examples/align_to_harp_clock/align_to_harp_clock.md
   - API:
     - Protocol: api/protocol.md
     - Serial: api/serial.md
diff --git a/src/packages/harp-data/README.md b/src/packages/harp-data/README.md
index 89d0677..8b4d3ab 100644
--- a/src/packages/harp-data/README.md
+++ b/src/packages/harp-data/README.md
@@ -86,6 +86,36 @@ _data, timestamps, _msg, payload = AnalogData.parse_bulk(raw)
 df = payload_to_dataframe(payload)
 ```
 
+## Align a non-Harp device to the Harp clock
+
+Devices outside the Harp bus keep their own clock. Some clock emitters mirror the
+[Synchronization Clock](https://harp-tech.org/protocol/SynchronizationClock.html)
+on a digital output at a much lower baud rate (typically 1 kbps rather than
+100 kbps) so that those devices can record it and be aligned post-hoc.
+`harp.data.synchronization` turns such a recording back into a table of anchors,
+keyed on whichever axis the device timestamps its own data on — the sample the
+packet was received at, or its local time in seconds — against the Harp second
+(`"Time"`) that packet carries:
+
+```python
+import numpy as np
+from harp.data.synchronization import decode_clock_from_samples
+
+clock = decode_clock_from_samples(sync_line, sample_rate=30_000.0)  # Sample -> Time
+harp_times = np.interp(spike_samples, clock.index, clock["Time"])
+
+# event-based systems report transitions instead, so anchors carry local seconds
+clock = decode_clock_from_transitions(edge_times, edge_states)      # LocalTime -> Time
+harp_times = np.interp(spike_times, clock.index, clock["Time"])
+```
+
+Packets that fail their start/stop bit check, or whose seconds do not add up
+against the local clock, are dropped — a glitched packet costs one anchor, not the
+alignment around it. By default anchors sit on the last transmitted bit of each
+packet, mirroring the protocol's synchronization event; pass
+`anchor="first_edge"` for emitters that align the whole second to the start of the
+transmission instead.
+
 ## Write data back out
 
 `to_file` / `to_buffer` are the inverse of the readers — encode values as Harp
diff --git a/src/packages/harp-data/src/harp/data/synchronization/__init__.py b/src/packages/harp-data/src/harp/data/synchronization/__init__.py
new file mode 100644
index 0000000..d95715a
--- /dev/null
+++ b/src/packages/harp-data/src/harp/data/synchronization/__init__.py
@@ -0,0 +1,15 @@
+"""Aligning a non-Harp device's local timestamps to the Harp clock."""
+
+from ._clock import (
+    DEFAULT_BAUD_RATE,
+    ClockAnchor,
+    decode_clock_from_samples,
+    decode_clock_from_transitions,
+)
+
+__all__ = [
+    "decode_clock_from_samples",
+    "decode_clock_from_transitions",
+    "ClockAnchor",
+    "DEFAULT_BAUD_RATE",
+]
diff --git a/src/packages/harp-data/src/harp/data/synchronization/_clock.py b/src/packages/harp-data/src/harp/data/synchronization/_clock.py
new file mode 100644
index 0000000..d0f696e
--- /dev/null
+++ b/src/packages/harp-data/src/harp/data/synchronization/_clock.py
@@ -0,0 +1,316 @@
+"""Decode the Harp Synchronization Clock from a downsampled digital signal.
+
+Harp devices share time over a dedicated serial bus running at 100 kbps (see the
+[Synchronization Clock protocol](https://harp-tech.org/protocol/SynchronizationClock.html)).
+Some emitters also mirror that signal on a digital output at a much lower baud rate —
+typically 1 kbps — so that it can be recorded by non-Harp acquisition systems whose
+sample rates are far below 100 kHz. The functions here decode such a recording back into
+Harp seconds, paired with the local time each of them was received at.
+
+The wire format is plain RS-232 without parity: the line idles high and every byte is
+sent as one low start bit, eight data bits (least significant first) and one high stop
+bit. The payload is the current Harp time in whole seconds (``uint32``, little-endian),
+optionally preceded by the ``0xAA 0xAF`` header of the full protocol packet — both
+framings are recognized.
+
+Two ways in, depending on what the recording system gives you:
+
+- a uniformly sampled waveform (digital or analog) → :func:`decode_clock_from_samples`
+- a list of line transitions ``(time, 0/1)`` → :func:`decode_clock_from_transitions`
+
+Each returns a table of anchors keyed on the axis that system reports its own timestamps
+on — a sample number or local seconds — against the whole Harp second received at it. How
+those timestamps are then interpolated onto the Harp axis is left to the caller.
+"""
+
+import warnings
+from typing import Any, Literal, Union, get_args
+
+import numpy as np
+import pandas as pd
+from numpy.typing import ArrayLike, NDArray
+
+#: Baud rate used by existing downsampled clock emitters, in bits per second.
+DEFAULT_BAUD_RATE = 1000.0
+
+#: Bits per byte on the wire: one start bit, eight data bits, one stop bit.
+_BITS_PER_BYTE = 10
+
+#: Bytes carrying the whole-second payload of a clock packet.
+_PAYLOAD_BYTES = 4
+
+#: Leading bytes of a full Harp Synchronization Clock packet, when present.
+_HEADER = (0xAA, 0xAF)
+
+#: How far, as a fraction of a bit period, a start bit may sit from its predicted
+#: position and still be used to re-synchronize the frame that follows it.
+_RESYNC_WINDOW = 0.3
+
+#: Where in the packet the whole second is taken to lapse. ``"last_bit"`` anchors on
+#: the end of the last transmitted bit, mirroring the Harp Synchronization Clock where
+#: the last byte carries the synchronization event; ``"first_edge"`` anchors on the
+#: first falling edge, which is what current emitters align to the second boundary.
+ClockAnchor = Literal["last_bit", "first_edge"]
+
+_ANCHORS = get_args(ClockAnchor)
+
+
+def decode_clock_from_transitions(
+    timestamps: ArrayLike,
+    states: ArrayLike,
+    *,
+    baud_rate: float = DEFAULT_BAUD_RATE,
+    anchor: ClockAnchor = "last_bit",
+    max_drift: Union[float, None] = 0.1,
+) -> pd.DataFrame:
+    """Decode Harp seconds from a list of clock-line transitions.
+
+    ``timestamps`` are the local times, in seconds, of each transition of the clock
+    line and ``states`` the level the line took at each of them (0 low, 1 high); this
+    is the shape most event-based acquisition systems report. The returned DataFrame
+    maps the local time each decoded packet is anchored on (the index, named
+    ``"LocalTime"``, see ``anchor``) to the whole Harp second it carries (``"Time"``,
+    the name the readers also give the Harp time axis).
+
+    ``baud_rate`` must match the emitter; bits are read at the center of each bit
+    period and every frame re-synchronizes on its own start bit, so a mismatch of a
+    few percent is tolerated. Packets whose start/stop bits do not check out are
+    dropped silently — a partially recorded or glitched packet costs one anchor, not
+    the decoding of the ones around it. ``max_drift`` additionally drops packets whose
+    seconds do not advance consistently with the local clock, allowing for at most
+    that fractional rate mismatch between the two (pass ``None`` to keep everything).
+    """
+    times, levels = _as_transitions(timestamps, states)
+    local, harp = _decode(times, levels, baud_rate=baud_rate, anchor=anchor)
+    local, harp = _drop_inconsistent(local, harp, max_drift)
+    return _clock_frame(pd.Index(local, name="LocalTime", dtype=np.float64), harp)
+
+
+def decode_clock_from_samples(
+    samples: ArrayLike,
+    sample_rate: float,
+    *,
+    threshold: Union[float, None] = None,
+    baud_rate: float = DEFAULT_BAUD_RATE,
+    anchor: ClockAnchor = "last_bit",
+    max_drift: Union[float, None] = 0.1,
+) -> pd.DataFrame:
+    """Decode Harp seconds from a uniformly sampled recording of the clock line.
+
+    ``samples`` is the recorded waveform and ``sample_rate`` its sampling frequency in
+    Hz; sample ``i`` is taken to occur at local time ``i / sample_rate``. Boolean and
+    integer input is treated as the digital line state (nonzero is high); pass
+    ``threshold`` to binarize an analog recording as ``samples >= threshold``.
+
+    The returned DataFrame maps the sample each decoded packet is anchored on (the
+    index, named ``"Sample"``, rounded to the nearest sample) to the whole Harp second
+    it carries (``"Time"``). See :func:`decode_clock_from_transitions` for ``anchor``
+    and ``max_drift``. About five samples per bit are needed for reliable decoding, so
+    a 1 kbps signal wants a sample rate of at least ~5 kHz.
+    """
+    if sample_rate <= 0:
+        raise ValueError(f"sample_rate must be positive, got {sample_rate}.")
+    if baud_rate <= 0:
+        raise ValueError(f"baud_rate must be positive, got {baud_rate}.")
+    if sample_rate < 2 * baud_rate:
+        raise ValueError(
+            f"A {baud_rate} bps signal cannot be decoded from samples taken at "
+            f"{sample_rate} Hz; at least 2 samples per bit are required (5 recommended)."
+        )
+
+    level = _as_levels(samples, threshold)
+    edges = np.flatnonzero(np.diff(level)) + 1
+    times = edges / sample_rate
+    local, harp = _decode(times, level[edges], baud_rate=baud_rate, anchor=anchor)
+    local, harp = _drop_inconsistent(local, harp, max_drift)
+    sample = np.rint(local * sample_rate).astype(np.int64)
+    return _clock_frame(pd.Index(sample, name="Sample"), harp)
+
+
+def _as_transitions(
+    timestamps: ArrayLike, states: ArrayLike
+) -> tuple[NDArray[np.float64], NDArray[np.bool_]]:
+    """Validate a transition list into sorted times and boolean levels."""
+    times = np.asarray(timestamps, dtype=np.float64)
+    levels = np.asarray(states)
+    if times.ndim != 1 or levels.ndim != 1:
+        raise ValueError("timestamps and states must be one-dimensional.")
+    if times.size != levels.size:
+        raise ValueError(
+            f"timestamps and states must have the same length, got {times.size} and {levels.size}."
+        )
+    if np.any(np.diff(times) < 0):
+        raise ValueError("timestamps must be sorted in ascending order.")
+    return times, levels.astype(np.bool_)
+
+
+def _as_levels(samples: ArrayLike, threshold: Union[float, None]) -> NDArray[np.bool_]:
+    """Binarize a sampled waveform into digital line levels."""
+    values = np.asarray(samples)
+    if values.ndim != 1:
+        raise ValueError("samples must be one-dimensional.")
+    if threshold is not None:
+        return values >= threshold
+    if values.dtype == np.bool_ or np.issubdtype(values.dtype, np.integer):
+        return values != 0
+    raise ValueError(
+        f"Pass a threshold to binarize samples of dtype {values.dtype}; only boolean "
+        "and integer samples are taken as digital line states."
+    )
+
+
+def _level_at(
+    times: NDArray[np.float64], levels: NDArray[np.bool_], at: NDArray[np.floating[Any]]
+) -> NDArray[np.bool_]:
+    """The line level at each time in ``at``, from the transitions ``(times, levels)``."""
+    index = np.searchsorted(times, at, side="right") - 1
+    # transitions alternate, so before the first one the line held the opposite level
+    return np.where(index < 0, ~levels[0], levels[np.maximum(index, 0)])
+
+
+def _resync(falling: NDArray[np.float64], expected: float, bit: float) -> float:
+    """Snap a predicted start-bit time to a nearby falling edge, as a UART would."""
+    best, window = expected, _RESYNC_WINDOW * bit
+    index = int(np.searchsorted(falling, expected))
+    for candidate in (index - 1, index):
+        if 0 <= candidate < falling.size:
+            distance = abs(falling[candidate] - expected)
+            if distance <= window:
+                best, window = float(falling[candidate]), distance
+    return best
+
+
+def _decode_bytes(
+    times: NDArray[np.float64],
+    levels: NDArray[np.bool_],
+    falling: NDArray[np.float64],
+    origin: float,
+    bit: float,
+    count: int,
+) -> Union[tuple[NDArray[np.uint8], NDArray[np.float64]], None]:
+    """Read ``count`` back-to-back RS-232 frames starting at ``origin``.
+
+    Returns the decoded bytes and the start-bit time of each frame, or ``None`` if any
+    of the frames is not properly delimited by a low start bit and a high stop bit.
+    """
+    values = np.empty(count, dtype=np.uint8)
+    origins = np.empty(count, dtype=np.float64)
+    offsets = (np.arange(_BITS_PER_BYTE) + 0.5) * bit
+    start = origin
+    for frame in range(count):
+        if frame:
+            # predict from the previous start bit, not from the packet, so that a small
+            # baud rate error stays a fraction of a bit instead of accumulating
+            start = _resync(falling, start + _BITS_PER_BYTE * bit, bit)
+        bits = _level_at(times, levels, start + offsets)
+        if bits[0] or not bits[-1]:  # start bit must be low, stop bit high
+            return None
+        origins[frame] = start
+        values[frame] = np.packbits(bits[1:-1], bitorder="little")[0]
+    return values, origins
+
+
+def _decode_packet(
+    times: NDArray[np.float64],
+    levels: NDArray[np.bool_],
+    falling: NDArray[np.float64],
+    origin: float,
+    bit: float,
+) -> Union[tuple[int, NDArray[np.float64]], None]:
+    """Decode one clock packet starting at ``origin``, with or without header."""
+    headed = _decode_bytes(times, levels, falling, origin, bit, len(_HEADER) + _PAYLOAD_BYTES)
+    if headed is not None and tuple(headed[0][: len(_HEADER)]) == _HEADER:
+        payload, origins = headed[0][len(_HEADER) :], headed[1]
+    else:
+        bare = _decode_bytes(times, levels, falling, origin, bit, _PAYLOAD_BYTES)
+        if bare is None:
+            return None
+        payload, origins = bare
+    return int.from_bytes(payload.tobytes(), "little"), origins
+
+
+def _decode(
+    times: NDArray[np.float64],
+    levels: NDArray[np.bool_],
+    *,
+    baud_rate: float,
+    anchor: ClockAnchor,
+) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
+    """Decode every clock packet in a transition list into (local time, Harp time)."""
+    if anchor not in _ANCHORS:
+        raise ValueError(f"anchor must be one of {_ANCHORS}, got {anchor!r}.")
+    if baud_rate <= 0:
+        raise ValueError(f"baud_rate must be positive, got {baud_rate}.")
+
+    bit = 1.0 / baud_rate
+    falling = times[~levels]  # a packet can only start on a falling edge
+    local: list[float] = []
+    harp: list[float] = []
+    decoded_until = -np.inf
+    for origin in falling:
+        if origin < decoded_until:
+            continue  # inside a packet we already decoded
+        packet = _decode_packet(times, levels, falling, float(origin), bit)
+        if packet is None:
+            continue
+        seconds, origins = packet
+        end = origins[-1] + _BITS_PER_BYTE * bit
+        decoded_until = end
+        local.append(end if anchor == "last_bit" else float(origin))
+        harp.append(float(seconds))
+    return np.array(local, dtype=np.float64), np.array(harp, dtype=np.float64)
+
+
+def _drop_inconsistent(
+    local: NDArray[np.float64], harp: NDArray[np.float64], max_drift: Union[float, None]
+) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
+    """Drop anchors whose Harp seconds disagree with the elapsed local time."""
+    if max_drift is None or local.size < 2:
+        return local, harp
+    if max_drift < 0:
+        raise ValueError(f"max_drift must be non-negative, got {max_drift}.")
+
+    def agrees(earlier: int, later: int) -> bool:
+        """Whether two anchors are one plausible stretch of the same clock apart."""
+        elapsed_harp = harp[later] - harp[earlier]
+        elapsed_local = local[later] - local[earlier]
+        return bool(
+            elapsed_harp >= 1
+            and abs(elapsed_harp - elapsed_local) <= max_drift * abs(elapsed_local)
+        )
+
+    # grow the longest chain of mutually consistent anchors out of the first agreeing
+    # pair, so that one corrupt value costs one anchor instead of every anchor after it
+    keep = np.zeros(local.size, dtype=np.bool_)
+    seed = next((i for i in range(local.size - 1) if agrees(i, i + 1)), None)
+    if seed is not None:
+        keep[seed] = keep[seed + 1] = True
+        last = seed + 1
+        for candidate in range(seed + 2, local.size):
+            if agrees(last, candidate):
+                keep[candidate] = True
+                last = candidate
+        first = seed
+        for candidate in range(seed - 1, -1, -1):
+            if agrees(candidate, first):
+                keep[candidate] = True
+                first = candidate
+
+    dropped = int(np.count_nonzero(~keep))
+    if dropped:
+        warnings.warn(
+            f"Dropped {dropped} clock packet(s) whose Harp time is inconsistent with the "
+            "local clock. Pass max_drift=None to keep them.",
+            stacklevel=3,
+        )
+    return local[keep], harp[keep]
+
+
+def _clock_frame(anchors: pd.Index, harp: NDArray[np.float64]) -> pd.DataFrame:
+    """Assemble the decoded packets into an ``anchors`` → Harp time table.
+
+    ``anchors`` locate each packet on the local device's own axis — whichever one it
+    reports timestamps on. "Time" is the Harp time axis throughout ``harp.data``; here
+    it is the column, because the index is the axis being translated from.
+    """
+    return pd.DataFrame({"Time": harp}, index=anchors)
diff --git a/tests/data/test_clock.py b/tests/data/test_clock.py
new file mode 100644
index 0000000..f157634
--- /dev/null
+++ b/tests/data/test_clock.py
@@ -0,0 +1,405 @@
+from contextlib import nullcontext
+
+import numpy as np
+import pytest
+from harp.data.synchronization import decode_clock_from_samples, decode_clock_from_transitions
+
+BAUD_RATE = 1000.0
+HEADER = (0xAA, 0xAF)
+
+
+def _packet_bits(value, *, header=False):
+    """The RS-232 bit stream of one clock packet: start, 8 data bits (LSB first), stop."""
+    payload = int(value).to_bytes(4, "little")
+    data = bytes(HEADER) + payload if header else payload
+    bits = []
+    for byte in data:
+        bits.append(0)
+        bits.extend((byte >> position) & 1 for position in range(8))
+        bits.append(1)
+    return bits
+
+
+def transitions(
+    values,
+    *,
+    baud_rate=BAUD_RATE,
+    emitted_baud_rate=None,
+    period=1.0,
+    first=0.5,
+    header=False,
+    jitter=0.0,
+    seed=0,
+):
+    """An idle-high clock line carrying one packet per ``period``.
+
+    Returns the transition times, the level taken at each of them, and the local time
+    of the first falling edge of every packet. ``emitted_baud_rate`` sets the rate the
+    signal is actually generated at, which may differ from the nominal ``baud_rate``.
+    """
+    bit = 1.0 / (emitted_baud_rate or baud_rate)
+    rng = np.random.default_rng(seed)
+    times, states, starts = [], [], []
+    level = 1
+    for packet, value in enumerate(values):
+        start = first + packet * period
+        starts.append(start)
+        for position, expected in enumerate(_packet_bits(value, header=header)):
+            if expected != level:
+                offset = rng.uniform(-jitter, jitter) * bit if jitter else 0.0
+                times.append(start + position * bit + offset)
+                states.append(expected)
+                level = expected
+    return np.array(times), np.array(states), np.array(starts)
+
+
+def render(times, states, sample_rate, *, duration=None):
+    """Sample an idle-high transition list onto a uniform grid of digital states."""
+    duration = duration if duration is not None else times[-1] + 0.5
+    grid = np.arange(int(round(duration * sample_rate))) / sample_rate
+    index = np.searchsorted(times, grid, side="right") - 1
+    return np.where(index < 0, ~states[0].astype(bool), states[np.maximum(index, 0)]).astype(bool)
+
+
+def packet_duration(*, baud_rate=BAUD_RATE, header=False):
+    return len(_packet_bits(0, header=header)) / baud_rate
+
+
+# ---------------------------------------------------------------- decoding
+
+
+@pytest.mark.parametrize("anchor", ["last_bit", "first_edge"])
+def test_decodes_a_single_packet(anchor):
+    times, states, starts = transitions([3806874])
+    clock = decode_clock_from_transitions(times, states, anchor=anchor)
+
+    assert clock.index.name == "LocalTime"
+    assert clock.columns.tolist() == ["Time"]
+    assert clock["Time"].tolist() == [3806874]
+    expected = starts[0] + (packet_duration() if anchor == "last_bit" else 0.0)
+    assert clock.index.to_numpy() == pytest.approx([expected])
+
+
+def test_decodes_consecutive_seconds():
+    values = [3806874, 3806875, 3806876, 3806877]
+    times, states, starts = transitions(values)
+
+    clock = decode_clock_from_transitions(times, states)
+
+    assert clock["Time"].tolist() == values
+    assert clock.index.to_numpy() == pytest.approx(starts + packet_duration())
+
+
+def test_header_framed_packets_are_detected():
+    values = [1000, 1001, 1002]
+    times, states, starts = transitions(values, header=True)
+
+    clock = decode_clock_from_transitions(times, states)
+
+    assert clock["Time"].tolist() == values
+    assert clock.index.to_numpy() == pytest.approx(starts + packet_duration(header=True))
+
+
+@pytest.mark.parametrize("value", [0, 1, 0xFF, 0xAFAA, 0xFF00FF00, 0xFFFFFFFE, 3806874])
+def test_round_trips_every_byte_pattern(value):
+    """Bytes of all zeros or all ones leave long runs without a transition."""
+    times, states, _ = transitions([value, value + 1])
+
+    clock = decode_clock_from_transitions(times, states)
+
+    assert clock["Time"].tolist() == [value, value + 1]
+
+
+@pytest.mark.parametrize("sample_rate", [5_000.0, 30_000.0])
+def test_sampled_recording_matches_the_transition_list(sample_rate):
+    values = [3806874, 3806875, 3806876]
+    times, states, starts = transitions(values)
+    samples = render(times, states, sample_rate)
+
+    clock = decode_clock_from_samples(samples, sample_rate)
+
+    assert clock["Time"].tolist() == values
+    assert clock.index.name == "Sample"
+    assert clock.index.dtype == np.int64
+    assert clock.index.tolist() == [
+        round(anchor * sample_rate) for anchor in starts + packet_duration()
+    ]
+
+
+def test_sample_index_points_at_the_anchor_sample():
+    sample_rate = 10_000.0
+    times, states, starts = transitions([42, 43], first=0.5)
+    samples = render(times, states, sample_rate)
+
+    clock = decode_clock_from_samples(samples, sample_rate, anchor="first_edge")
+
+    assert clock.columns.tolist() == ["Time"]
+    # the anchor sample is the first sample of the packet's start bit
+    for sample, start in zip(clock.index, starts):
+        assert sample == round(start * sample_rate)
+        assert not samples[sample]
+        assert samples[sample - 1]
+
+
+def test_analog_recording_is_binarized_by_threshold():
+    sample_rate = 10_000.0
+    values = [7, 8]
+    times, states, _ = transitions(values)
+    analog = np.where(render(times, states, sample_rate), 3.3, 0.1) + 0.02
+
+    clock = decode_clock_from_samples(analog, sample_rate, threshold=1.5)
+
+    assert clock["Time"].tolist() == values
+
+
+def test_float_samples_without_threshold_raises():
+    with pytest.raises(ValueError, match="threshold"):
+        decode_clock_from_samples(np.zeros(100, dtype=np.float64), 10_000.0)
+
+
+def test_sample_rate_below_two_samples_per_bit_raises():
+    with pytest.raises(ValueError, match="cannot be decoded"):
+        decode_clock_from_samples(np.zeros(100, dtype=bool), 1500.0, baud_rate=BAUD_RATE)
+
+
+def test_recording_without_any_transition_yields_no_anchors():
+    clock = decode_clock_from_samples(np.ones(10_000, dtype=bool), 10_000.0)
+
+    assert clock.empty
+    assert clock.columns.tolist() == ["Time"]
+    assert clock.index.name == "Sample"
+
+
+def test_empty_transition_list_yields_no_anchors():
+    clock = decode_clock_from_transitions([], [])
+
+    assert clock.empty
+    assert clock.columns.tolist() == ["Time"]
+    assert clock["Time"].dtype == np.float64
+
+
+def test_wrong_baud_rate_yields_no_anchors():
+    times, states, _ = transitions([100, 101, 102])
+
+    clock = decode_clock_from_transitions(times, states, baud_rate=9600.0)
+
+    assert clock.empty
+    assert clock.index.name == "LocalTime"
+
+
+# ---------------------------------------------------------------- robustness
+
+
+def test_glitches_between_packets_are_ignored():
+    values = [500, 501]
+    times, states, starts = transitions(values)
+    # a 100 us spike in the idle line, well away from any packet
+    times = np.concatenate([times, [1.4, 1.4001]])
+    states = np.concatenate([states, [0, 1]])
+    order = np.argsort(times)
+
+    clock = decode_clock_from_transitions(times[order], states[order])
+
+    assert clock["Time"].tolist() == values
+    assert clock.index.to_numpy() == pytest.approx(starts + packet_duration())
+
+
+def test_packet_truncated_by_the_end_of_the_recording_is_dropped():
+    values = [500, 501, 502]
+    times, states, starts = transitions(values)
+    keep = times < starts[-1] + 0.5 * packet_duration()
+
+    clock = decode_clock_from_transitions(times[keep], states[keep])
+
+    assert clock["Time"].tolist() == values[:-1]
+
+
+def test_recording_starting_mid_packet_keeps_the_later_packets():
+    values = [500, 501, 502]
+    times, states, starts = transitions(values)
+    keep = times > starts[0] + 0.5 * packet_duration()
+
+    clock = decode_clock_from_transitions(times[keep], states[keep])
+
+    assert clock["Time"].tolist() == values[1:]
+
+
+def test_tolerates_jitter_on_every_edge():
+    values = [3806874, 3806875, 3806876]
+    times, states, starts = transitions(values, jitter=0.2)
+
+    clock = decode_clock_from_transitions(times, states)
+
+    assert clock["Time"].tolist() == values
+    assert clock.index.to_numpy() == pytest.approx(starts + packet_duration(), abs=0.3 / BAUD_RATE)
+
+
+def test_tolerates_a_baud_rate_error_that_would_slip_a_whole_bit():
+    """Per-frame re-synchronization keeps a 3% rate error from walking off the bits."""
+    values = [3806874, 3806875]
+    times, states, _ = transitions(values, emitted_baud_rate=BAUD_RATE * 1.03)
+
+    clock = decode_clock_from_transitions(times, states, baud_rate=BAUD_RATE)
+
+    assert clock["Time"].tolist() == values
+
+
+def test_corrupt_packet_is_dropped_with_a_warning():
+    values = [3806874, 3806875, 12345, 3806877]
+    times, states, starts = transitions(values)
+
+    with pytest.warns(UserWarning, match="inconsistent with the local clock"):
+        clock = decode_clock_from_transitions(times, states)
+
+    assert clock["Time"].tolist() == [3806874, 3806875, 3806877]
+    assert clock.index.to_numpy() == pytest.approx(np.delete(starts, 2) + packet_duration())
+
+
+def test_corrupt_first_packet_is_dropped_with_a_warning():
+    values = [12345, 3806875, 3806876, 3806877]
+    times, states, _ = transitions(values)
+
+    with pytest.warns(UserWarning):
+        clock = decode_clock_from_transitions(times, states)
+
+    assert clock["Time"].tolist() == values[1:]
+
+
+def test_max_drift_none_keeps_every_decoded_packet():
+    values = [3806874, 3806875, 12345, 3806877]
+    times, states, _ = transitions(values)
+
+    clock = decode_clock_from_transitions(times, states, max_drift=None)
+
+    assert clock["Time"].tolist() == values
+
+
+def test_missing_packets_do_not_break_the_chain():
+    values = [3806874, 3806875, 3806876, 3806877]
+    times, states, starts = transitions(values)
+    dropped = (times < starts[1]) | (times >= starts[2])
+
+    clock = decode_clock_from_transitions(times[dropped], states[dropped])
+
+    assert clock["Time"].tolist() == [3806874, 3806876, 3806877]
+
+
+def test_invalid_transition_lists_raise():
+    times, states, _ = transitions([1, 2])
+    with pytest.raises(ValueError, match="same length"):
+        decode_clock_from_transitions(times, states[:-1])
+    with pytest.raises(ValueError, match="sorted"):
+        decode_clock_from_transitions(times[::-1], states)
+    with pytest.raises(ValueError, match="one-dimensional"):
+        decode_clock_from_transitions(np.zeros((2, 2)), np.zeros((2, 2)))
+    with pytest.raises(ValueError, match="anchor"):
+        decode_clock_from_transitions(times, states, anchor="middle")  # type: ignore[arg-type]
+    with pytest.raises(ValueError, match="baud_rate"):
+        decode_clock_from_transitions(times, states, baud_rate=0.0)
+
+
+def test_a_long_drifting_sampled_recording_decodes_end_to_end():
+    """A minute of a 30 kHz recording: local clock drift, jitter and one lost packet."""
+    rate = 1.001  # Harp seconds per local second
+    sample_rate = 30_000.0
+    values = [3806874 + second for second in range(60)]
+    times, states, starts = transitions(
+        values, period=1.0 / rate, baud_rate=BAUD_RATE * rate, jitter=0.15
+    )
+    lost = (times >= starts[30]) & (times < starts[30] + 1.5 * packet_duration())
+    samples = render(times[~lost], states[~lost], sample_rate)
+
+    clock = decode_clock_from_samples(
+        samples, sample_rate, baud_rate=BAUD_RATE * rate, anchor="first_edge"
+    )
+
+    assert clock["Time"].tolist() == [value for value in values if value != values[30]]
+    # anchors land on the emitted start bits, within the jitter and one sample period
+    assert clock.index.to_numpy() / sample_rate == pytest.approx(np.delete(starts, 30), abs=5e-4)
+
+
+# ---------------------------------------------------------------- recorded signals
+
+# Transitions acquired by an Open Ephys system from a Behavior board emitting the
+# clock at 1 kbps, taken from https://github.com/harp-tech/harp-python/pull/38.
+# fmt: off
+RECORDED = [
+    {
+        # two valid packets, followed by one cut short by the end of the recording
+        "timestamps": np.array([
+            0.        , 0.96106667, 0.96306667, 0.96406667, 0.96506667,
+            0.96706667, 0.96906667, 0.97106667, 0.97306667, 0.97506667,
+            0.97606667, 0.97706667, 0.98006667, 0.98106667, 0.98306667,
+            0.98406667, 0.98506667, 0.98806667, 0.99006667, 0.99106667,
+            1.00006667, 1.96116667, 1.96216667, 1.96416667, 1.96516667,
+            1.96716667, 1.96916667, 1.97116667, 1.97316667, 1.97516667,
+            1.97616667, 1.97716667, 1.98016667, 1.98116667, 1.98316667,
+            1.98416667, 1.98516667, 1.98816667, 1.99016667, 1.99116667,
+            2.00016667, 2.96126667, 2.96426667, 2.96726667, 2.96926667,
+            2.97126667, 2.97326667, 2.97526667, 2.97626667, 2.97726667,
+            2.98026667, 2.98126667, 2.98326667, 2.98426667, 2.98526667,
+            2.98826667, 2.99026667, 2.99126667]),
+        "expected_start_times": np.array([0.96106667, 1.96116667]),
+        "expected_harp_times": [3806874, 3806875],
+        # the cut-short packet fails its framing check, so nothing reaches the drift check
+        "expected_warning": False,
+    },
+    {
+        # four valid packets and one corrupt one
+        "timestamps": np.array([
+            0.14036667, 1.10146667, 1.10246667, 1.10346667, 1.10446667,
+            1.11146667, 1.11246667, 1.11646667, 1.11746667, 1.11846667,
+            1.11946667, 1.12146667, 1.12246667, 1.12546667, 1.12746667,
+            1.12846667, 1.13046667, 1.13146667, 1.14046667, 2.10156667,
+            2.10356667, 2.11156667, 2.11256667, 2.11656667, 2.11756667,
+            2.11856667, 2.11956667, 2.12156667, 2.12256667, 2.12556667,
+            2.12756667, 2.12856667, 2.13056667, 2.13156667, 2.14056667,
+            3.10163333, 3.10263333, 3.11163333, 3.11263333, 3.11663333,
+            3.11763333, 3.11863333, 3.11963333, 3.12163333, 3.12263333,
+            3.12563333, 3.12763333, 3.12863333, 3.13063333, 3.13163333,
+            3.14063333, 4.10173333, 4.11073333, 4.11173333, 4.11673333,
+            4.11873333, 4.11973333, 4.12173333, 4.12273333, 4.12573333,
+            4.12773333, 4.12873333, 4.13073333, 4.13173333, 4.14073333,
+            5.1018    , 5.1028    , 5.1038    , 5.1108    , 5.1118    ,
+            5.1168    , 5.1188    , 5.1198    , 5.1218    , 5.1228    ,
+            5.1258    , 5.1278    , 5.1288    , 5.1308    , 5.1318    ,
+            5.1368    , 5.1388    , 5.1398    , 5.14183333, 5.1428    ,
+            5.14583333, 5.14783333, 5.14883333, 5.15083333, 5.15183333,
+            5.16083333, 6.1059    ]),
+        "expected_start_times": np.array([1.10146667, 2.10156667, 3.10163333, 4.10173333]),
+        "expected_harp_times": [2600957, 2600958, 2600959, 2600960],
+        # the corrupt packet decodes cleanly but its seconds do not add up
+        "expected_warning": True,
+    },
+]
+# fmt: on
+
+
+def _decode_recorded(recorded, **kwargs):
+    times = recorded["timestamps"]
+    states = np.resize([1, 0], times.size)  # transitions alternate, starting from high
+    expected = (
+        pytest.warns(UserWarning, match="inconsistent with the local clock")
+        if recorded["expected_warning"]
+        else nullcontext()
+    )
+    with expected:
+        return decode_clock_from_transitions(times, states, **kwargs)
+
+
+@pytest.mark.parametrize("recorded", RECORDED)
+def test_decodes_a_recorded_clock_signal(recorded):
+    clock = _decode_recorded(recorded, anchor="first_edge")
+
+    assert clock["Time"].tolist() == recorded["expected_harp_times"]
+    assert clock.index.to_numpy() == pytest.approx(recorded["expected_start_times"])
+
+
+@pytest.mark.parametrize("recorded", RECORDED)
+def test_anchors_a_recorded_signal_on_its_last_bit(recorded):
+    clock = _decode_recorded(recorded)
+
+    assert clock["Time"].tolist() == recorded["expected_harp_times"]
+    assert clock.index.to_numpy() == pytest.approx(
+        recorded["expected_start_times"] + packet_duration()
+    )