Add harp.data.synchronization for aligning non-Harp devices to the Harp clock

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

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-->

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"])

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

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

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",
+]