"""Tests for tde.py. The load-bearing checks, per the arc's standing rules: - the delay CRLB closed form is verified against a numerically built Fisher information (finite differences of the actual delayed mean vector) AND against an independent analytic special case (the Gaussian pulse), because a bound quoted from memory in the wrong parameterisation is this arc's recurring trap; - the two-sensor pairwise bound is verified against a numerically built 2x2 Fisher matrix, including the matrix's predicted structure [[A+B, B], [B, B]]; - the estimator's variance is Monte-Carlo-compared to the bound, and the documented failure modes (naive unsigned-frequency bound, PHAT's low-SNR fragility) are demonstrated, not just asserted away. """ import numpy as np import pytest from scipy.signal import lfilter from tde import (cross_correlation, centroid_delay, crlb_delay, crlb_delay_pair, delay_fisher_information_numeric, estimate_delay, fractional_delay, gcc, parabolic_interpolation, rms_bandwidth, signal_centroid) from tde import _bandlimited_derivative def bandlimited_pulse(n=512, f0=0.05, f1=0.20, centre=0.5, width=0.06): """A Gaussian-windowed chirp, well contained in the record.""" t = np.arange(n, dtype=float) tc = centre * n rate = (f1 - f0) / (2 * 0.4 * n) phase = 2 * np.pi * (f0 * (t - tc) + rate * (t - tc) ** 2) return np.exp(-0.5 * ((t - tc) / (width * n)) ** 2) * np.cos(phase) class TestFractionalDelay: def test_integer_delay_is_a_roll(self): rng = np.random.default_rng(0) x = rng.standard_normal(64) assert np.allclose(fractional_delay(x, 5.0), np.roll(x, 5), atol=1e-12) def test_subsample_delay_matches_analytic_shift(self): # A band-limited waveform sampled at shifted instants is the # exact answer the FFT method must reproduce. t = np.arange(512, dtype=float) def sig(tv): return (np.exp(-0.5 * ((tv - 256) / 30.0) ** 2) * np.cos(2 * np.pi * 0.1 * (tv - 256))) got = fractional_delay(sig(t), 3.37) assert np.max(np.abs(got - sig(t - 3.37))) < 1e-10 def test_negative_delay_inverts(self): x = bandlimited_pulse() back = fractional_delay(fractional_delay(x, 2.6), -2.6) assert np.max(np.abs(back - x)) < 1e-10 def test_rejects_bad_input(self): with pytest.raises(ValueError): fractional_delay(np.zeros((4, 4)), 1.0) with pytest.raises(ValueError): fractional_delay(np.zeros(8), 1.0, fs=0.0) class TestCrossCorrelation: def test_peak_sits_at_positive_lag_for_delayed_copy(self): # Orientation is the easiest thing to get backwards: y lagging # x by D must put the peak at +D, matching R_xy(tau). rng = np.random.default_rng(1) x = rng.standard_normal(256) y = np.roll(x, 9) lags, r = cross_correlation(x, y) assert lags[np.argmax(r)] == 9 def test_max_lag_restricts_the_window(self): rng = np.random.default_rng(2) x = rng.standard_normal(128) lags, r = cross_correlation(x, x, max_lag=10) assert lags.min() == -10 and lags.max() == 10 assert lags[np.argmax(r)] == 0 def test_matches_direct_sum(self): x = np.array([1.0, 2.0, -1.0, 0.5]) y = np.array([0.0, 1.0, 2.0, -1.0]) lags, r = cross_correlation(x, y) for lag, val in zip(lags, r): direct = sum(x[n] * y[n + lag] for n in range(len(x)) if 0 <= n + lag < len(y)) assert abs(val - direct) < 1e-12 def test_rejects_empty(self): with pytest.raises(ValueError): cross_correlation(np.array([]), np.array([1.0])) class TestGcc: def test_direct_weighting_matches_cross_correlation(self): rng = np.random.default_rng(3) x = rng.standard_normal(200) y = rng.standard_normal(200) lags_a, r_a = cross_correlation(x, y, max_lag=50) lags_b, r_b = gcc(x, y, weighting="direct", max_lag=50) assert np.array_equal(lags_a, lags_b) assert np.allclose(r_a, r_b, atol=1e-9) def test_phat_sharpens_the_peak_under_multipath(self): # Colored source plus a strong echo: the direct correlation # peak is a wide blob (the source's own autocorrelation # convolved with the channel); PHAT whitens and resolves it. rng = np.random.default_rng(4) src = lfilter([1.0], [1.0, -0.95], rng.standard_normal(4096)) y = fractional_delay(src, 23.0) + 0.8 * fractional_delay(src, 61.0) widths = {} for wgt in ("direct", "phat"): lags, r = gcc(src, y, weighting=wgt, max_lag=100) widths[wgt] = int(np.sum(r / r.max() > 0.5)) assert lags[np.argmax(r)] == 23 assert widths["phat"] * 10 <= widths["direct"] def test_phat_degrades_faster_at_low_snr(self): # The documented price of whitening: noise-only bins vote too. rng = np.random.default_rng(5) src = lfilter([1.0], [1.0, -0.95], rng.standard_normal(4096)) yd = fractional_delay(src, 23.0) sig = np.sqrt(np.var(src) / 10 ** (-10 / 10)) # -10 dB SNR err_d, err_p = [], [] for _ in range(20): xx = src + sig * rng.standard_normal(len(src)) yy = yd + sig * rng.standard_normal(len(src)) err_d.append(abs(estimate_delay(xx, yy, max_lag=100) - 23.0)) err_p.append(abs(estimate_delay(xx, yy, weighting="phat", max_lag=100) - 23.0)) assert np.median(err_p) > np.median(err_d) def test_rejects_mismatched_lengths(self): with pytest.raises(ValueError): gcc(np.zeros(8), np.zeros(9)) class TestParabolicInterpolation: def test_recovers_exact_quadratic_vertex(self): delta_true = 0.31 i = np.arange(7, dtype=float) y = -(i - (3 + delta_true)) ** 2 assert abs(parabolic_interpolation(y, 3) - delta_true) < 1e-12 def test_refuses_non_maximum(self): assert parabolic_interpolation(np.array([0.0, 0.0, 0.0]), 1) == 0.0 assert parabolic_interpolation(np.array([2.0, 1.0, 2.0]), 1) == 0.0 def test_rejects_edge_index(self): with pytest.raises(ValueError): parabolic_interpolation(np.array([1.0, 2.0, 1.0]), 0) class TestEstimateDelay: def test_noise_free_subsample_accuracy(self): s = bandlimited_pulse() got = estimate_delay(s, fractional_delay(s, 7.37), max_lag=50) assert abs(got - 7.37) < 0.02 def test_seconds_scaling(self): s = bandlimited_pulse() fs = 48000.0 got = estimate_delay(s, fractional_delay(s, 12.0), fs=fs) assert abs(got - 12.0 / fs) < 1e-6 def test_negative_delay(self): s = bandlimited_pulse() got = estimate_delay(s, fractional_delay(s, -4.6), max_lag=50) assert abs(got - (-4.6)) < 0.02 class TestCrlbDelay: """The standing rule: never ship a bound checked only against itself.""" def test_closed_form_matches_numeric_fisher_information(self): s = bandlimited_pulse() sigma2 = 0.09 j_closed = 1.0 / crlb_delay(s, sigma2) j_numeric = delay_fisher_information_numeric(s, sigma2) assert abs(j_closed / j_numeric - 1.0) < 1e-6 def test_gaussian_pulse_analytic_special_case(self): # Independent closed form: for s(t) = A exp(-t^2 / 2w^2), # integral of s'(t)^2 dt = sqrt(pi) A^2 / (4 w) * ... worked: # s' = -(t/w^2) s, int t^2 exp(-t^2/w^2) dt = sqrt(pi) w^3 / 2, # so int s'^2 = A^2 sqrt(pi) / (2 w). With Ts = 1 the sampled # sum approximates the integral directly. amp, w = 1.3, 8.0 t = np.arange(512, dtype=float) - 256 g = amp * np.exp(-t ** 2 / (2 * w ** 2)) analytic = np.sqrt(np.pi) * amp ** 2 / (2 * w) sigma2 = 0.25 assert abs(1.0 / crlb_delay(g, sigma2) - analytic / sigma2) < 1e-6 * analytic / sigma2 def test_equals_inverse_enr_times_beta_squared(self): s = bandlimited_pulse() sigma2 = 0.04 beta = rms_bandwidth(s) enr = np.sum(s ** 2) / sigma2 assert abs(crlb_delay(s, sigma2) * beta ** 2 * enr - 1.0) < 1e-9 def test_unsigned_frequency_reading_is_badly_wrong(self): # The trap this module's docstrings warn about: summing # (2 pi i / N)^2 over DFT indices 0..N-1 treats the negative # frequencies as ever-higher positive ones and wildly inflates # the Fisher information for a real lowpass signal. s = bandlimited_pulse() sigma2 = 0.09 n = len(s) spec = np.abs(np.fft.fft(s)) ** 2 w_naive = 2 * np.pi * np.arange(n) / n j_naive = np.sum(w_naive ** 2 * spec) / n / sigma2 j_true = 1.0 / crlb_delay(s, sigma2) assert j_naive > 50 * j_true def test_scales_with_noise_and_sample_rate(self): s = bandlimited_pulse() assert np.isclose(crlb_delay(s, 0.2) / crlb_delay(s, 0.1), 2.0) # Same waveform played fs times faster: derivative energy up # by fs^2 per sample, so the variance drops by fs^2. assert np.isclose(crlb_delay(s, 0.1, fs=100.0) * 100.0 ** 2, crlb_delay(s, 0.1, fs=1.0)) def test_monte_carlo_estimator_sits_on_the_bound(self): # The correlation-peak estimator is the MLE; at comfortable # ENR its measured variance must sit on (not under, not far # above) the bound. rng = np.random.default_rng(6) s = bandlimited_pulse() enr = 10 ** (25 / 10) sigma2 = np.sum(s ** 2) / enr delayed = fractional_delay(s, 7.37) est = [estimate_delay(s, delayed + rng.standard_normal(len(s)) * np.sqrt(sigma2), max_lag=50) for _ in range(400)] ratio = np.var(est) / crlb_delay(s, sigma2) assert 0.8 < ratio < 1.5 def test_rejects_bad_input(self): s = bandlimited_pulse() with pytest.raises(ValueError): crlb_delay(s, 0.0) with pytest.raises(ValueError): crlb_delay(np.zeros(64), 1.0) # zero energy -> no information class TestCrlbDelayPair: def test_numeric_two_parameter_fisher_matrix(self): # Build J for theta = (tau1, dtau) by finite differences of # the actual two-channel mean, invert it, and compare # [J^-1]_22 to the closed form 1/A + 1/B. This is the check # that resolves the archived derivation's "which is strange" # ending: 1/J_22 alone is NOT the bound. s = bandlimited_pulse(n=600) sig1, sig2, gain = 0.15, 0.4, 0.6 h = 1e-3 def mean(tau1, dtau): return (fractional_delay(s, tau1), gain * fractional_delay(s, tau1 + dtau)) def dmean(k): args = [(h, 5.0), (-h, 5.0)] if k == 0 else [(0.0, 5.0 + h), (0.0, 5.0 - h)] (a1, b1), (a2, b2) = mean(*args[0]), mean(*args[1]) return (a1 - a2) / (2 * h), (b1 - b2) / (2 * h) d = [dmean(0), dmean(1)] jm = np.array([[np.sum(d[i][0] * d[j][0]) / sig1 ** 2 + np.sum(d[i][1] * d[j][1]) / sig2 ** 2 for j in range(2)] for i in range(2)]) closed = crlb_delay_pair(s, sig1 ** 2, sig2 ** 2, gain=gain) assert abs(np.linalg.inv(jm)[1, 1] / closed - 1.0) < 1e-4 # The matrix structure the archived note derived: with # A and B the single-channel informations, # J = [[A + B, B], [B, B]]. a = 1.0 / crlb_delay(s, sig1 ** 2) b = 1.0 / crlb_delay(gain * s, sig2 ** 2) assert abs(jm[0, 1] / jm[1, 1] - 1.0) < 1e-6 assert abs(jm[1, 1] / b - 1.0) < 1e-4 assert abs((jm[0, 0] - jm[0, 1]) / a - 1.0) < 1e-4 def test_equal_channels_double_the_single_channel_bound(self): # The factor-of-two trap: per-sensor bound != pairwise bound. s = bandlimited_pulse() single = crlb_delay(s, 0.09) pair = crlb_delay_pair(s, 0.09, 0.09) assert np.isclose(pair, 2.0 * single) def test_clean_reference_recovers_single_channel_bound(self): s = bandlimited_pulse() near_clean = crlb_delay_pair(s, 1e-12, 0.09) assert abs(near_clean / crlb_delay(s, 0.09) - 1.0) < 1e-6 class TestCentroid: def test_shift_equivariance(self): # The property that makes centroid differences measure delays. pulse = bandlimited_pulse(n=2048, centre=0.25) got = centroid_delay(pulse, fractional_delay(pulse, 13.71)) assert abs(got - 13.71) < 1e-9 def test_centroid_of_symmetric_pulse_is_its_centre(self): t = np.arange(512, dtype=float) g = np.exp(-0.5 * ((t - 200.0) / 12.0) ** 2) assert abs(signal_centroid(g) - 200.0) < 1e-9 def test_equals_energy_weighted_mean_group_delay(self): # The Parseval-for-moments identity from the archived centroid # derivation: send the pulse through an all-pass channel with # known analytic group delay tau0 + 3 b w^2 and check the # centroid moves by exactly the |S|^2-weighted mean of it. n = 2048 t = np.arange(n, dtype=float) pulse = (np.exp(-0.5 * ((t - 500) / 40.0) ** 2) * np.cos(2 * np.pi * 0.08 * (t - 500))) w = 2 * np.pi * np.fft.fftfreq(n) tau0, b = 200.0, 800.0 channel = np.exp(-1j * (tau0 * w + b * w ** 3)) out = np.fft.ifft(np.fft.fft(pulse) * channel).real spec = np.abs(np.fft.fft(pulse)) ** 2 predicted = np.sum(spec * (tau0 + 3 * b * w ** 2)) / np.sum(spec) assert abs(centroid_delay(pulse, out) - predicted) < 1e-6 def test_seconds_scaling(self): t = np.arange(512, dtype=float) g = np.exp(-0.5 * ((t - 200.0) / 12.0) ** 2) assert np.isclose(signal_centroid(g, fs=1000.0), 0.2) def test_rejects_zero_energy(self): with pytest.raises(ValueError): signal_centroid(np.zeros(32)) class TestRmsBandwidth: def test_gaussian_pulse_analytic_value(self): # For a Gaussian envelope exp(-t^2/2w^2), beta^2 = # int s'^2 / int s^2 = 1/(2 w^2) (radians, Ts = 1). w = 8.0 t = np.arange(512, dtype=float) - 256 g = np.exp(-t ** 2 / (2 * w ** 2)) assert abs(rms_bandwidth(g) - 1.0 / (w * np.sqrt(2))) < 1e-6 def test_carrier_dominates_for_narrowband(self): # A slow envelope on a fast carrier: beta is set by the # carrier, not the envelope (the accuracy/ambiguity story). t = np.arange(4096, dtype=float) env = np.exp(-0.5 * ((t - 2048) / 400.0) ** 2) f0 = 0.2 beta = rms_bandwidth(env * np.cos(2 * np.pi * f0 * (t - 2048))) assert abs(beta - 2 * np.pi * f0) < 0.02 * 2 * np.pi * f0 def test_derivative_identity(self): s = bandlimited_pulse() ds = _bandlimited_derivative(s, 1.0) assert np.isclose(rms_bandwidth(s) ** 2, np.sum(ds ** 2) / np.sum(s ** 2))