diff --git a/backend/festival4d/audio_sync.py b/backend/festival4d/audio_sync.py index 0a0a8a6..5797ee0 100644 --- a/backend/festival4d/audio_sync.py +++ b/backend/festival4d/audio_sync.py @@ -68,29 +68,15 @@ def _peak_ratio(env: np.ndarray, peak: int, sr: int, guard_s: float = 0.004) -> return main / second -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 - SIG = np.fft.rfft(sig, n) - REF = np.fft.rfft(ref, n) +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 - cc = np.fft.irfft(R, n) # lag 0 at index 0; negative lags wrap to the tail + return np.fft.irfft(R, n) # lag 0 at index 0; negative lags wrap to the tail - 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)) + +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) @@ -107,6 +93,27 @@ def gcc_phat(sig: np.ndarray, ref: np.ndarray, sr: int, 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 # --------------------------------------------------------------------------- @@ -116,13 +123,30 @@ def pairwise_offsets(signals: dict[int, np.ndarray], sr: int) -> list[dict]: ``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] - delay_s, conf = gcc_phat(signals[a], signals[b], sr) + 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 diff --git a/backend/tests/test_audio_sync.py b/backend/tests/test_audio_sync.py index d1fbcd9..5ee496e 100644 --- a/backend/tests/test_audio_sync.py +++ b/backend/tests/test_audio_sync.py @@ -88,6 +88,28 @@ def test_pairwise_and_solve_recover_synthetic_offsets(): assert abs(solved[i + 1] - off) <= 10.0, (off, solved[i + 1]) +def test_pairwise_matches_direct_gcc_phat_with_unequal_lengths(): + """The cached-FFT pairwise path agrees with a direct per-pair gcc_phat call. + + Signals get different lengths on purpose: the shared padded size (2·max_len) then + differs from the per-pair size (la+lb), which is exactly where the two paths could + diverge if the lag windowing were wrong. + """ + sr = 16_000 + master = _broadband(sr * 8, seed=7) + signals = { + 1: master[: sr * 6], + 2: _delay_samples(master, 900)[: sr * 4], + 3: _delay_samples(master, -1500)[: sr * 5], + } + edges = audio_sync.pairwise_offsets(signals, sr) + assert len(edges) == 3 + for e in edges: + direct_delay, direct_conf = audio_sync.gcc_phat(signals[e["a"]], signals[e["b"]], sr) + assert abs(e["offset_s"] - -direct_delay) < 1.0 / sr, e + assert e["confidence"] > 5.0 and direct_conf > 5.0 + + def test_pairwise_edge_sign_is_relative_start_offset(): """Edge ``offset_s`` = offset_a − offset_b (seconds).""" sr = config.AUDIO_SAMPLE_RATE