95b71788a3
pairwise_offsets now computes each signal's rFFT once at a shared padded length (2*max_len) and reuses it for every pair: O(K) forward FFTs + cheap per-pair spectrum products instead of O(K^2) full FFTs. Each pair's lag search stays capped at +/-(la+lb)/2 — the same window the per-pair transform used — and gcc_phat() itself is unchanged for callers (drift windows, tests). Measured 1.8x on 6 cams x 10 min; the advantage grows with camera count. New test locks the cached path against direct gcc_phat, including unequal-length signals where the padded sizes genuinely differ. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
439 lines
18 KiB
Python
439 lines
18 KiB
Python
"""GCC-PHAT pairwise audio offsets + global graph solve + drift (spec M1, lane A).
|
||
|
||
Algorithm (spec M1):
|
||
|
||
- **GCC-PHAT**, not raw cross-correlation: whiten the cross-spectrum
|
||
``R = X·conj(Y); R /= |R| + eps; r = irfft(R)`` — robust to loud/clipped/reverberant
|
||
concert audio. numpy FFTs; librosa only for loading (in :func:`run_sync`).
|
||
- All **pairwise** offsets, each scored by peak-to-second-peak ratio; then a global
|
||
least-squares solve over the offset graph, weighted by confidence, rejecting pairs that
|
||
violate cycle consistency (``|off_ab + off_bc − off_ac| > 50 ms``, generalized here as a
|
||
residual test against the fitted global solution).
|
||
- **Drift**: re-run GCC-PHAT on 10 s windows every 30 s along the overlap; fit a line to
|
||
``delay(t)``; store the slope (dimensionless rate) as ``drift_ppm`` (store 0 if < 5 ppm).
|
||
- Disconnected sync-graph components → ``offset_ms = None`` (never a guess; spec pitfall #5).
|
||
- Persist via ``db.update_video_sync`` and export ``config.SYNC_JSON``.
|
||
|
||
Sign conventions (kept consistent with ``config``'s timebase):
|
||
|
||
A video *v*'s audio sample at local time ``τ`` equals the master signal at global time
|
||
``offset_s + τ`` (drift 0). So for two videos ``a``/``b``:
|
||
``offset_a − offset_b = −gcc_phat(a, b)`` seconds, where :func:`gcc_phat` returns the
|
||
delay of its first argument relative to its second. Under drift ``d`` (= ``drift_ppm·1e−6``),
|
||
the windowed delay of ``v`` vs the reference is ``d·t − offset_s`` — a line whose slope is
|
||
the drift, recovered by :func:`estimate_drift`.
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import json
|
||
import logging
|
||
|
||
import numpy as np
|
||
|
||
from festival4d import config, db
|
||
from festival4d.ingest import audio_wav_path
|
||
|
||
log = logging.getLogger("festival4d.audio_sync")
|
||
|
||
_EPS = 1e-10
|
||
# Edges weaker than this peak-to-second-peak ratio carry no reliable alignment
|
||
# (flat/ambiguous correlation) and are dropped before the global solve.
|
||
_MIN_EDGE_CONFIDENCE = 1.15
|
||
# A pairwise offset whose residual against the fitted global solution exceeds this is a
|
||
# cycle-consistency violation (spec M1); it is rejected and the component is re-solved.
|
||
_CYCLE_TOLERANCE_MS = 50.0
|
||
# Drift smaller than this magnitude is stored as exactly 0 (spec M1).
|
||
_DRIFT_DEADBAND_PPM = 5.0
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Core correlation
|
||
# ---------------------------------------------------------------------------
|
||
def _peak_ratio(env: np.ndarray, peak: int, sr: int, guard_s: float = 0.004) -> float:
|
||
"""Peak-to-second-peak ratio of a correlation envelope (the spec's edge score).
|
||
|
||
The second peak is the largest value outside a small guard band around the main peak,
|
||
so the main lobe's own shoulders don't count as competition.
|
||
"""
|
||
main = float(env[peak])
|
||
if main <= _EPS:
|
||
return 0.0
|
||
guard = max(1, int(round(guard_s * sr)))
|
||
masked = env.copy()
|
||
masked[max(0, peak - guard):min(len(env), peak + guard + 1)] = 0.0
|
||
second = float(masked.max()) if masked.size else 0.0
|
||
if second <= _EPS:
|
||
return 1.0e6
|
||
return main / second
|
||
|
||
|
||
def _whitened_cc(SIG: np.ndarray, REF: np.ndarray, n: int) -> np.ndarray:
|
||
"""PHAT-whitened cross-correlation from precomputed rFFTs (lag 0 at index 0)."""
|
||
R = SIG * np.conj(REF)
|
||
R /= np.abs(R) + _EPS # PHAT weighting: flatten magnitude, keep phase
|
||
return np.fft.irfft(R, n) # lag 0 at index 0; negative lags wrap to the tail
|
||
|
||
|
||
def _peak_offset(cc: np.ndarray, sr: int, max_lag: int) -> tuple[float, float]:
|
||
"""Peak of a whitened correlation → ``(offset_s, confidence)`` over ±``max_lag`` samples."""
|
||
# Re-center to contiguous lags [−max_lag … +max_lag].
|
||
cc = np.concatenate((cc[-max_lag:], cc[:max_lag + 1]))
|
||
env = np.abs(cc)
|
||
peak = int(np.argmax(env))
|
||
# Parabolic sub-sample refinement of the peak location (3-point fit around the max),
|
||
# so precision isn't capped at one sample (62.5 µs at 16 kHz) and drift slopes stay clean.
|
||
delta = 0.0
|
||
if 0 < peak < len(env) - 1:
|
||
a, b, c = env[peak - 1], env[peak], env[peak + 1]
|
||
denom = a - 2.0 * b + c
|
||
if abs(denom) > _EPS:
|
||
delta = float(np.clip(0.5 * (a - c) / denom, -0.5, 0.5))
|
||
offset_s = ((peak - max_lag) + delta) / sr
|
||
return offset_s, _peak_ratio(env, peak, sr)
|
||
|
||
|
||
def gcc_phat(sig: np.ndarray, ref: np.ndarray, sr: int,
|
||
max_lag_s: float | None = None) -> tuple[float, float]:
|
||
"""GCC-PHAT time delay of ``sig`` relative to ``ref`` (spec M1, lane A).
|
||
|
||
Returns ``(offset_s, confidence)`` where ``offset_s`` is positive when ``sig`` *lags*
|
||
``ref`` (i.e. ``sig[n] ≈ ref[n − offset_s·sr]``), and confidence is the
|
||
peak-to-second-peak ratio.
|
||
"""
|
||
sig = np.asarray(sig, dtype=np.float64)
|
||
ref = np.asarray(ref, dtype=np.float64)
|
||
# Remove DC so the whitened correlation keys on structure, not a bias term.
|
||
sig = sig - sig.mean()
|
||
ref = ref - ref.mean()
|
||
|
||
n = len(sig) + len(ref) # linear (zero-padded) correlation → no circular wrap
|
||
cc = _whitened_cc(np.fft.rfft(sig, n), np.fft.rfft(ref, n), n)
|
||
max_lag = n // 2 if max_lag_s is None else int(round(max_lag_s * sr))
|
||
max_lag = max(1, min(max_lag, n // 2))
|
||
return _peak_offset(cc, sr, max_lag)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Pairwise graph
|
||
# ---------------------------------------------------------------------------
|
||
def pairwise_offsets(signals: dict[int, np.ndarray], sr: int) -> list[dict]:
|
||
"""All pairwise GCC-PHAT offsets with confidences (spec M1, lane A).
|
||
|
||
``signals`` maps ``video_id -> mono samples``. Returns a list of
|
||
``{a, b, offset_s, confidence}`` edges, where ``offset_s`` is the *relative start
|
||
offset* ``offset_a − offset_b`` in seconds (see the module's sign conventions).
|
||
|
||
Each signal's rFFT is computed **once** at a shared padded length and reused across
|
||
every pair — O(K) FFTs + O(K²) cheap spectrum products, instead of the O(K²) full
|
||
FFTs a per-pair :func:`gcc_phat` would cost. The shared length (2·max_len ≥ any
|
||
pair's la+lb) keeps the correlation linear (no circular wrap), and each pair's lag
|
||
search stays capped at ±(la+lb)/2 — the same window the per-pair transform used.
|
||
"""
|
||
ids = sorted(signals)
|
||
if len(ids) < 2:
|
||
return []
|
||
lengths = {i: len(signals[i]) for i in ids}
|
||
n = 2 * max(lengths.values())
|
||
spectra: dict[int, np.ndarray] = {}
|
||
for i in ids:
|
||
s = np.asarray(signals[i], dtype=np.float64)
|
||
s = s - s.mean() # same DC removal as gcc_phat
|
||
spectra[i] = np.fft.rfft(s, n)
|
||
edges: list[dict] = []
|
||
for i in range(len(ids)):
|
||
for j in range(i + 1, len(ids)):
|
||
a, b = ids[i], ids[j]
|
||
cc = _whitened_cc(spectra[a], spectra[b], n)
|
||
max_lag = max(1, min(n // 2, (lengths[a] + lengths[b]) // 2))
|
||
delay_s, conf = _peak_offset(cc, sr, max_lag)
|
||
edges.append({"a": a, "b": b, "offset_s": -delay_s, "confidence": conf})
|
||
return edges
|
||
|
||
|
||
def _components(ids: list[int], edges: list[dict]) -> dict[int, list[int]]:
|
||
"""Connected components of the offset graph (union-find over usable edges)."""
|
||
parent = {i: i for i in ids}
|
||
|
||
def find(x: int) -> int:
|
||
while parent[x] != x:
|
||
parent[x] = parent[parent[x]]
|
||
x = parent[x]
|
||
return x
|
||
|
||
for e in edges:
|
||
parent[find(e["a"])] = find(e["b"])
|
||
comps: dict[int, list[int]] = {}
|
||
for i in ids:
|
||
comps.setdefault(find(i), []).append(i)
|
||
return comps
|
||
|
||
|
||
def _connected(nodes: list[int], root: int, edges: list[dict]) -> bool:
|
||
"""Whether every node is reachable from ``root`` over ``edges``."""
|
||
adj: dict[int, list[int]] = {n: [] for n in nodes}
|
||
for e in edges:
|
||
adj[e["a"]].append(e["b"])
|
||
adj[e["b"]].append(e["a"])
|
||
seen = {root}
|
||
stack = [root]
|
||
while stack:
|
||
for y in adj[stack.pop()]:
|
||
if y not in seen:
|
||
seen.add(y)
|
||
stack.append(y)
|
||
return len(seen) == len(nodes)
|
||
|
||
|
||
def _weighted_lstsq(nodes: list[int], ref: int, edges: list[dict]) -> dict[int, float]:
|
||
"""Confidence-weighted least-squares offsets (seconds) with ``ref`` pinned to 0.
|
||
|
||
Each edge contributes the constraint ``offset_a − offset_b = offset_s``, weighted by
|
||
``sqrt(confidence)``. Solved for all nodes except the pinned reference.
|
||
"""
|
||
unknown = [n for n in nodes if n != ref]
|
||
col = {n: k for k, n in enumerate(unknown)}
|
||
A = np.zeros((len(edges), len(unknown)))
|
||
b = np.zeros(len(edges))
|
||
w = np.zeros(len(edges))
|
||
for r, e in enumerate(edges):
|
||
if e["a"] in col:
|
||
A[r, col[e["a"]]] += 1.0
|
||
if e["b"] in col:
|
||
A[r, col[e["b"]]] -= 1.0
|
||
b[r] = e["offset_s"] # ref terms are 0, so they drop out of A
|
||
w[r] = np.sqrt(max(e["confidence"], _EPS))
|
||
x, *_ = np.linalg.lstsq(A * w[:, None], b * w, rcond=None)
|
||
sol = {ref: 0.0}
|
||
for n in unknown:
|
||
sol[n] = float(x[col[n]])
|
||
return sol
|
||
|
||
|
||
def _solve_component(nodes: list[int], edges: list[dict]) -> dict[int, float]:
|
||
"""Solve one connected component, rejecting cycle-inconsistent edges iteratively.
|
||
|
||
The reference (smallest ``video_id``) is pinned to 0. While the graph is over-determined
|
||
(more edges than a spanning tree needs) and the worst edge residual exceeds
|
||
``_CYCLE_TOLERANCE_MS``, drop that edge (unless doing so would disconnect a node) and
|
||
re-solve. Returns offsets in seconds.
|
||
"""
|
||
nodes = sorted(nodes)
|
||
ref = min(nodes)
|
||
edges = list(edges)
|
||
sol = _weighted_lstsq(nodes, ref, edges)
|
||
while len(edges) > len(nodes) - 1:
|
||
worst, worst_res = None, 0.0
|
||
for e in edges:
|
||
res = abs((sol[e["a"]] - sol[e["b"]]) - e["offset_s"])
|
||
if res > worst_res:
|
||
worst, worst_res = e, res
|
||
if worst is None or worst_res * 1000.0 <= _CYCLE_TOLERANCE_MS:
|
||
break
|
||
trial = [e for e in edges if e is not worst]
|
||
if not _connected(nodes, ref, trial):
|
||
break # can't drop it without orphaning a node
|
||
log.info("audio_sync: rejecting inconsistent edge %d-%d (residual %.1f ms)",
|
||
worst["a"], worst["b"], worst_res * 1000.0)
|
||
edges = trial
|
||
sol = _weighted_lstsq(nodes, ref, edges)
|
||
return sol
|
||
|
||
|
||
def solve_global_offsets(edges: list[dict], video_ids: list[int]) -> dict[int, float | None]:
|
||
"""Least-squares global solve over the pairwise-offset graph (spec M1, lane A).
|
||
|
||
The reference component (largest; ties broken by smallest ``video_id``) is solved with
|
||
its smallest-id node pinned to 0. Every video outside it — including any with no usable
|
||
edges — gets ``None`` (spec pitfall #5). Returns ``video_id -> offset_ms`` (or ``None``).
|
||
"""
|
||
ids = sorted(set(video_ids))
|
||
usable = [
|
||
e for e in edges
|
||
if e["a"] in ids and e["b"] in ids
|
||
and e["confidence"] >= _MIN_EDGE_CONFIDENCE
|
||
and np.isfinite(e["offset_s"])
|
||
]
|
||
result: dict[int, float | None] = {i: None for i in ids}
|
||
if not ids:
|
||
return result
|
||
|
||
comps = _components(ids, usable)
|
||
# Reference component: the one that aligns the most cameras. A lone reference (no
|
||
# overlap with anyone) is still "aligned" at 0; everyone else stays None.
|
||
ref_comp = max(comps.values(), key=lambda c: (len(c), -min(c)))
|
||
if len(ref_comp) == 1:
|
||
result[ref_comp[0]] = 0.0
|
||
return result
|
||
|
||
ref_set = set(ref_comp)
|
||
comp_edges = [e for e in usable if e["a"] in ref_set and e["b"] in ref_set]
|
||
for i, off_s in _solve_component(ref_comp, comp_edges).items():
|
||
result[i] = off_s * 1000.0
|
||
return result
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Drift
|
||
# ---------------------------------------------------------------------------
|
||
def estimate_drift(sig: np.ndarray, ref: np.ndarray, sr: int,
|
||
window_s: float = 10.0, hop_s: float = 30.0) -> float:
|
||
"""Estimate clock drift in ppm by fitting window-offsets over the overlap (spec M1).
|
||
|
||
Slides a ``window_s`` window every ``hop_s`` over the common span, measures each
|
||
window's GCC-PHAT delay, and fits a line ``delay ≈ slope·t + b``. ``slope`` is the
|
||
dimensionless drift rate; ``drift_ppm = slope·1e6``. Returns 0 when there are too few
|
||
windows to fit or the estimate is within the ±5 ppm deadband.
|
||
"""
|
||
sig = np.asarray(sig, dtype=np.float64)
|
||
ref = np.asarray(ref, dtype=np.float64)
|
||
n = min(len(sig), len(ref))
|
||
win = int(round(window_s * sr))
|
||
hop = int(round(hop_s * sr))
|
||
if win < sr or n < win or hop < 1: # need a window of at least ~1 s that fits
|
||
return 0.0
|
||
|
||
centers: list[float] = []
|
||
delays: list[float] = []
|
||
start = 0
|
||
while start + win <= n:
|
||
delay_s, conf = gcc_phat(sig[start:start + win], ref[start:start + win], sr,
|
||
max_lag_s=window_s / 2.0)
|
||
if conf >= _MIN_EDGE_CONFIDENCE:
|
||
centers.append((start + win / 2.0) / sr)
|
||
delays.append(delay_s)
|
||
start += hop
|
||
|
||
if len(centers) < 2:
|
||
return 0.0
|
||
slope = float(np.polyfit(np.array(centers), np.array(delays), 1)[0])
|
||
drift_ppm = slope * 1.0e6
|
||
return 0.0 if abs(drift_ppm) < _DRIFT_DEADBAND_PPM else drift_ppm
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Orchestration
|
||
# ---------------------------------------------------------------------------
|
||
def _reference_id(offsets_ms: dict[int, float | None]) -> int | None:
|
||
"""The reference video: smallest id among those with a solved (non-None) offset."""
|
||
solved = [i for i, o in offsets_ms.items() if o is not None]
|
||
return min(solved) if solved else None
|
||
|
||
|
||
def _confidence_by_video(edges: list[dict], offsets_ms: dict[int, float | None],
|
||
ref_id: int | None) -> dict[int, float | None]:
|
||
"""Per-video sync confidence in [0, 1): best incident edge, mapped ``1 − 1/ratio``.
|
||
|
||
The reference is 1.0; unsolved videos are None. A raw peak ratio of 1 (ambiguous) maps
|
||
to 0, 2 to 0.5, large ratios approach 1 — a bounded, honest confidence.
|
||
"""
|
||
best: dict[int, float] = {}
|
||
for e in edges:
|
||
if offsets_ms.get(e["a"]) is None or offsets_ms.get(e["b"]) is None:
|
||
continue
|
||
r = e["confidence"]
|
||
best[e["a"]] = max(best.get(e["a"], 0.0), r)
|
||
best[e["b"]] = max(best.get(e["b"], 0.0), r)
|
||
out: dict[int, float | None] = {}
|
||
for vid, off in offsets_ms.items():
|
||
if off is None:
|
||
out[vid] = None
|
||
elif vid == ref_id:
|
||
out[vid] = 1.0
|
||
else:
|
||
ratio = best.get(vid)
|
||
out[vid] = None if ratio is None else float(max(0.0, 1.0 - 1.0 / ratio))
|
||
return out
|
||
|
||
|
||
def run_sync() -> dict:
|
||
"""Sync every ingested video: pairwise GCC-PHAT, global solve, drift, persist.
|
||
|
||
Entrypoint for ``python -m festival4d sync``. Writes ``offset_ms``/``drift_ppm``/
|
||
``sync_confidence`` via ``db`` and exports ``config.SYNC_JSON``. Returns the solution.
|
||
"""
|
||
import librosa # spec: librosa only for loading; import lazily (heavy module)
|
||
|
||
db.init_engine()
|
||
db.init_db()
|
||
videos = db.get_videos()
|
||
if not videos:
|
||
raise RuntimeError("no videos in DB — run `python -m festival4d ingest` first")
|
||
|
||
sr = config.AUDIO_SAMPLE_RATE
|
||
signals: dict[int, np.ndarray] = {}
|
||
missing: list[str] = []
|
||
for v in videos:
|
||
wav = audio_wav_path(v.id)
|
||
if not wav.exists():
|
||
missing.append(v.filename)
|
||
continue
|
||
y, _ = librosa.load(str(wav), sr=sr, mono=True)
|
||
signals[v.id] = y.astype(np.float64)
|
||
if missing:
|
||
log.warning("audio_sync: no extracted audio for %s — run ingest; "
|
||
"they will be left unsynced (offset=None)", ", ".join(missing))
|
||
|
||
edges = pairwise_offsets(signals, sr) if len(signals) >= 2 else []
|
||
offsets_ms = solve_global_offsets(edges, list(signals.keys()))
|
||
# Videos with no audio at all are unsynced too.
|
||
for v in videos:
|
||
offsets_ms.setdefault(v.id, None)
|
||
|
||
ref_id = _reference_id(offsets_ms)
|
||
ref_sig = signals.get(ref_id) if ref_id is not None else None
|
||
confidence = _confidence_by_video(edges, offsets_ms, ref_id)
|
||
|
||
solution = {
|
||
"reference_video_id": ref_id,
|
||
"reference_filename": next((v.filename for v in videos if v.id == ref_id), None),
|
||
"sample_rate": sr,
|
||
"videos": [],
|
||
"edges": [
|
||
{"a": e["a"], "b": e["b"],
|
||
"offset_ms": round(e["offset_s"] * 1000.0, 3),
|
||
"confidence": round(float(e["confidence"]), 4)}
|
||
for e in edges
|
||
],
|
||
"components": sorted(
|
||
(sorted(c) for c in _components(sorted(signals), [
|
||
e for e in edges if e["confidence"] >= _MIN_EDGE_CONFIDENCE
|
||
]).values()),
|
||
key=lambda c: (-len(c), c),
|
||
),
|
||
"generated_by": "festival4d.audio_sync",
|
||
}
|
||
|
||
for v in videos:
|
||
off = offsets_ms.get(v.id)
|
||
if off is None:
|
||
drift = None
|
||
elif v.id == ref_id:
|
||
drift = 0.0
|
||
elif ref_sig is not None and v.id in signals:
|
||
drift = estimate_drift(signals[v.id], ref_sig, sr)
|
||
else:
|
||
drift = 0.0
|
||
conf = confidence.get(v.id)
|
||
db.update_video_sync(v.id, off, drift, conf)
|
||
solution["videos"].append({
|
||
"video_id": v.id, "filename": v.filename,
|
||
"offset_ms": None if off is None else round(off, 3),
|
||
"drift_ppm": None if drift is None else round(drift, 4),
|
||
"sync_confidence": None if conf is None else round(conf, 4),
|
||
"connected": off is not None,
|
||
})
|
||
log.info("audio_sync: %s offset=%s drift=%s conf=%s",
|
||
v.filename,
|
||
"None" if off is None else f"{off:+.1f}ms",
|
||
"None" if drift is None else f"{drift:+.2f}ppm",
|
||
"None" if conf is None else f"{conf:.2f}")
|
||
|
||
config.SYNC_JSON.parent.mkdir(parents=True, exist_ok=True)
|
||
config.SYNC_JSON.write_text(json.dumps(solution, indent=2))
|
||
log.info("audio_sync: wrote %s (%d synced, %d unsynced)",
|
||
config.SYNC_JSON,
|
||
sum(1 for vv in solution["videos"] if vv["connected"]),
|
||
sum(1 for vv in solution["videos"] if not vv["connected"]))
|
||
return solution
|