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