"""Time-delay estimation: correlation peaks, subsample interpolation, bounds. Every "where" question a signal can answer is a "when" question in disguise: a sonar range is an echo delay, a direction of arrival is an inter-sensor delay, a cable-fault location is a difference in arrival times. This module estimates that delay from sampled data: the cross-correlation peak, its subsample refinement, the generalised cross-correlation (GCC-PHAT) variant for multipath, the Cramer-Rao bound that says how well any of them can do, and the signal-centroid delay for dispersive channels where "the" delay stops being a single number. Conventions, stated up front per the arc's standing rule -------------------------------------------------------- - ``estimate_delay(x, y)`` returns the delay of *y* relative to *x*: positive when y is a delayed copy of x. This matches the cross-correlation convention R_xy(tau) = E[X(t) Y(t+tau)], which peaks at tau = +D when y lags x by D. - The delay CRLB is quoted against the **mean-square radian bandwidth** beta^2 = sum_i (2 pi f_i)^2 |S_i|^2 / sum_i |S_i|^2 taken over **signed** frequencies (numpy.fft.fftfreq), as the second moment about zero, carrier included. Evaluating the same sum over unsigned DFT indices 0..N-1 silently inflates the Fisher information for a real signal (the upper half of the spectrum is the *negative* frequencies, not ever-higher positive ones); the tests demonstrate the failure. And the signal-strength axis is the energy-to-noise ratio Es/sigma^2 of the detection-theory topic, NOT a per-sample SNR: mislabeling that axis moves every curve by 10 log10(N) dB. - All delays are in seconds for a given ``fs``; with the default ``fs = 1.0`` they are in samples. """ import numpy as np from scipy.signal import correlate, correlation_lags def fractional_delay(x: np.ndarray, delay: float, fs: float = 1.0) -> np.ndarray: """Delay a signal by an arbitrary, possibly subsample, amount. Multiplies the spectrum by exp(-j 2 pi f * delay) and inverse transforms: the exact delay operator for a band-limited, periodic signal. The operation is **circular**, so use it on records with enough zero padding that nothing meaningful wraps around the ends; energy pushed past the last sample reappears at the first. For real input with even length the Nyquist bin has no conjugate partner, so a non-integer delay cannot keep it real; taking the real part of the result is exact whenever the signal has no energy at Nyquist, which any properly band-limited record satisfies. Parameters ---------- x : ndarray Input record (real or complex). delay : float Delay in seconds (samples if fs = 1); may be negative. fs : float Sample rate in Hz. Returns ------- ndarray The delayed record, same length and dtype family as x. """ x = np.asarray(x) if x.ndim != 1: raise ValueError("x must be one-dimensional") if fs <= 0: raise ValueError("fs must be strictly positive") f = np.fft.fftfreq(len(x), d=1.0 / fs) shifted = np.fft.ifft(np.fft.fft(x) * np.exp(-2j * np.pi * f * delay)) return shifted.real if np.isrealobj(x) else shifted def cross_correlation(x: np.ndarray, y: np.ndarray, max_lag: int | None = None ) -> tuple[np.ndarray, np.ndarray]: """Linear cross-correlation r[m] = sum_n x[n] y[n+m]. The sample counterpart of R_xy(tau) = E[X(t) Y(t+tau)]: when y is a delayed copy of x, the peak sits at the **positive** lag equal to the delay. Computed via FFT (zero-padded, so linear, not circular). Parameters ---------- x, y : ndarray The two records; lengths may differ. max_lag : int, optional Keep only lags in [-max_lag, +max_lag]. Restricting the search range to the physically possible delays (e.g. +-d/c samples for a microphone pair) is not an optimisation but part of the estimator: it rejects correlation peaks that no true delay could produce. Returns ------- lags : ndarray of int Lag of y relative to x, in samples. r : ndarray Correlation value at each lag. """ x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) if x.ndim != 1 or y.ndim != 1: raise ValueError("x and y must be one-dimensional") if len(x) == 0 or len(y) == 0: raise ValueError("x and y must be non-empty") # scipy's correlate(y, x) computes sum_n y[n + m] x[n] at lag m. r = correlate(y, x, mode="full", method="fft") lags = correlation_lags(len(y), len(x), mode="full") if max_lag is not None: if max_lag < 1: raise ValueError("max_lag must be at least 1") keep = np.abs(lags) <= max_lag lags, r = lags[keep], r[keep] return lags, r def gcc(x: np.ndarray, y: np.ndarray, weighting: str = "direct", max_lag: int | None = None) -> tuple[np.ndarray, np.ndarray]: """Generalised cross-correlation of two equal-length records. Computes the cross-spectrum G(f) = X*(f) Y(f) on a zero-padded grid, applies a frequency weighting, and inverse transforms: - ``'direct'``: no weighting; identical (to rounding) to :func:`cross_correlation`. - ``'phat'``: the phase transform, G(f)/|G(f)|. Every frequency bin keeps only its phase, so the result is the correlation of spectrally *whitened* signals: the peak sharpens toward an impulse and multipath copies separate into distinct peaks instead of one smeared blob. The price is that bins containing only noise vote with the same weight as bins full of signal, which is why PHAT degrades faster than direct correlation as SNR falls. Bins with |G| below 1e-12 of the maximum are zeroed rather than divided. Parameters ---------- x, y : ndarray Equal-length records. weighting : str 'direct' or 'phat'. max_lag : int, optional As in :func:`cross_correlation`. Returns ------- lags : ndarray of int Lag of y relative to x, in samples. r : ndarray (Weighted) correlation at each lag. """ x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) if x.shape != y.shape or x.ndim != 1: raise ValueError("x and y must be one-dimensional and equal length") if weighting not in ("direct", "phat"): raise ValueError("weighting must be 'direct' or 'phat'") n = len(x) if n == 0: raise ValueError("x and y must be non-empty") nfft = 2 * n # zero padding: linear correlation, lags -n+1 .. n-1 g = np.conj(np.fft.fft(x, nfft)) * np.fft.fft(y, nfft) if weighting == "phat": mag = np.abs(g) floor = 1e-12 * mag.max() g = np.where(mag > floor, g / np.maximum(mag, floor), 0.0) r_circ = np.fft.ifft(g).real lag_span = n - 1 if max_lag is None else min(max_lag, n - 1) if max_lag is not None and max_lag < 1: raise ValueError("max_lag must be at least 1") lags = np.arange(-lag_span, lag_span + 1) r = r_circ[lags % nfft] return lags, r def parabolic_interpolation(r: np.ndarray, i: int) -> float: """Subsample offset of a peak from the parabola through three points. Fits a parabola through (i-1, r[i-1]), (i, r[i]), (i+1, r[i+1]) and returns the vertex offset delta = (r[i-1] - r[i+1]) / (2 (r[i-1] - 2 r[i] + r[i+1])), in (-0.5, 0.5) when r[i] is a strict local maximum. This is the same move as the FFT-peak interpolators of the estimating-a-sinusoid topic, applied on the lag axis instead of the frequency axis, and it inherits the same caveat: the parabola is a model of the peak's shape, and it is only as good as the match between that model and the true correlation peak (good for a smooth, oversampled peak; biased when the peak is cusp-like or barely sampled). Parameters ---------- r : ndarray Correlation (or any) sequence. i : int Index of the sample peak; must have both neighbours. Returns ------- float Offset delta such that the interpolated peak is at i + delta. Returns 0.0 when the three points do not form a maximum (non-negative curvature), rather than extrapolating nonsense. """ r = np.asarray(r, dtype=float) if not 1 <= i <= len(r) - 2: raise ValueError("peak index must have a neighbour on each side") denom = r[i - 1] - 2.0 * r[i] + r[i + 1] if denom >= 0.0: return 0.0 return float(0.5 * (r[i - 1] - r[i + 1]) / denom) def estimate_delay(x: np.ndarray, y: np.ndarray, fs: float = 1.0, max_lag: int | None = None, weighting: str = "direct", interpolate: bool = True) -> float: """Delay of y relative to x from the (generalised) correlation peak. Finds the largest correlation lag, optionally refines it with :func:`parabolic_interpolation`, and converts to seconds. This is the maximum-likelihood estimator for a known waveform in white Gaussian noise (the matched-filter peak), and the workhorse estimator for the two-noisy-sensor case. Parameters ---------- x, y : ndarray The two records (equal length required for 'phat'). fs : float Sample rate in Hz; the returned delay is in seconds. max_lag : int, optional Restrict the peak search to |lag| <= max_lag samples. weighting : str 'direct' or 'phat' (see :func:`gcc`). interpolate : bool Apply parabolic subsample interpolation around the peak. Returns ------- float Estimated delay in seconds (samples if fs = 1). """ if fs <= 0: raise ValueError("fs must be strictly positive") if weighting == "direct" and len(x) != len(y): lags, r = cross_correlation(x, y, max_lag=max_lag) else: lags, r = gcc(x, y, weighting=weighting, max_lag=max_lag) i = int(np.argmax(r)) delta = 0.0 if interpolate and 1 <= i <= len(r) - 2: delta = parabolic_interpolation(r, i) return float((lags[i] + delta) / fs) def rms_bandwidth(s: np.ndarray, fs: float = 1.0) -> float: """Root-mean-square (Gabor) bandwidth, in radians per second. beta = sqrt( sum (2 pi f_i)^2 |S_i|^2 / sum |S_i|^2 ) over signed frequencies: the square root of the second moment of the energy spectrum about **zero** frequency, carrier included. This is the exact quantity the delay CRLB contains: a bandpass signal at a high carrier has a large beta (delay information rides on carrier phase) even when its envelope is slow, which is precisely the accuracy-versus-ambiguity trade of narrowband time-delay estimation. Parameters ---------- s : ndarray The (real) signal record. fs : float Sample rate in Hz. Returns ------- float RMS bandwidth in rad/s. """ s = np.asarray(s, dtype=float) if s.ndim != 1 or len(s) < 2: raise ValueError("s must be a one-dimensional record") if fs <= 0: raise ValueError("fs must be strictly positive") spec = np.abs(np.fft.fft(s)) ** 2 e = spec.sum() if e == 0: raise ValueError("signal has zero energy") w = 2.0 * np.pi * np.fft.fftfreq(len(s), d=1.0 / fs) return float(np.sqrt(np.sum(w ** 2 * spec) / e)) def _bandlimited_derivative(s: np.ndarray, fs: float) -> np.ndarray: """ds/dt of the band-limited interpolation of s, via the FFT.""" w = 2.0 * np.pi * np.fft.fftfreq(len(s), d=1.0 / fs) return np.fft.ifft(1j * w * np.fft.fft(s)).real def crlb_delay(s: np.ndarray, sigma2: float, fs: float = 1.0) -> float: """Cramer-Rao bound for the delay of a known waveform, one channel. Model: x[k] = s(k Ts - tau) + w[k], with s known, w white Gaussian of per-sample variance sigma2. The Fisher information is J(tau) = (1/sigma2) sum_k s'(k Ts - tau)^2 = (1/sigma2) (1/N) sum_i (2 pi f_i)^2 |S_i|^2 (Parseval, signed frequencies f_i), so var(tau_hat) >= sigma2 / sum_k s'(k)^2 = 1 / (beta^2 * Es / sigma2) with beta the :func:`rms_bandwidth` and Es = sum s^2: the bound falls with the energy-to-noise ratio like any CRLB, and with the *bandwidth* squared, because delay information lives in how fast the waveform changes. Sharp edges locate; slow ripples do not. Parameterisation warnings, per the arc's standing rule: Es/sigma2 is a total-energy ratio, not a per-sample SNR (factor N); the continuous-time literature writes the same bound as 1/(beta^2 * 2E/N0) where the factor 2 is absorbed by the two-sided noise density N0/2 (sigma2 = (N0/2) fs for ideally band-limited white noise); and this is the ONE-channel, reference-known bound: see :func:`crlb_delay_pair` before applying it to a sensor pair. The tests check this closed form against a numerically built Fisher information and against an independent analytic special case. Parameters ---------- s : ndarray The known waveform, sampled; assumed well contained in the record (negligible energy at the ends and at Nyquist). sigma2 : float Noise variance per sample. fs : float Sample rate in Hz. Returns ------- float Lower bound on the delay variance, in seconds squared. """ s = np.asarray(s, dtype=float) if s.ndim != 1 or len(s) < 2: raise ValueError("s must be a one-dimensional record") if sigma2 <= 0: raise ValueError("sigma2 must be strictly positive") if fs <= 0: raise ValueError("fs must be strictly positive") ds = _bandlimited_derivative(s, fs) j = np.sum(ds ** 2) / sigma2 if j == 0: raise ValueError("signal carries no delay information (zero bandwidth)") return float(1.0 / j) def delay_fisher_information_numeric(s: np.ndarray, sigma2: float, fs: float = 1.0, step: float = 1e-3) -> float: """Fisher information for the delay, built numerically. The check the arc's standing rule demands: instead of trusting the closed form, differentiate the model itself. The mean vector is mu(tau) = s(k Ts - tau); its derivative is approximated by a central difference of two fractionally delayed copies, dmu/dtau ~= (mu(tau + h) - mu(tau - h)) / (2 h), and J = sum (dmu/dtau)^2 / sigma2. No Parseval, no analytic derivative: if this and :func:`crlb_delay` disagree, the closed form is wrong (this catch is what found the sinusoid CRLB's factor-of-two slip in piece 2 of the arc, and the naive unsigned-frequency reading of the delay bound here). Parameters ---------- s : ndarray The known waveform, sampled, well contained in the record. sigma2 : float Noise variance per sample. fs : float Sample rate in Hz. step : float Half-width h of the central difference, in samples. Returns ------- float Fisher information J, in 1/seconds^2; the CRLB is 1/J. """ if step <= 0: raise ValueError("step must be strictly positive") if sigma2 <= 0: raise ValueError("sigma2 must be strictly positive") h = step / fs dmu = (fractional_delay(s, +h, fs) - fractional_delay(s, -h, fs)) / (2 * h) return float(np.sum(dmu ** 2) / sigma2) def crlb_delay_pair(s: np.ndarray, sigma2_ref: float, sigma2_del: float, fs: float = 1.0, gain: float = 1.0) -> float: """CRLB for the delay *difference* between two noisy sensors. Model: the reference sensor sees x1[k] = s(k Ts - tau1) + w1[k] and the second sensor sees x2[k] = gain * s(k Ts - tau1 - dtau) + w2[k], with s known but the reference arrival time tau1 **unknown** and estimated jointly with dtau. The 2x2 Fisher matrix is [[A + B, B], [B, B]] with A and B the single-channel informations of the two sensors, and inverting it gives var(dtau_hat) >= [J^-1]_22 = 1/A + 1/B : the individual channel bounds **add**. With equal channels the pairwise bound is twice the single-channel bound: quoting a per-sensor bound for a sensor-pair measurement is a silent factor-of-two error, the delay estimation edition of this arc's recurring parameterisation trap. Note 1/B alone (what you would get by pretending tau1 were known) is NOT the bound; the origin of that distinction is worked through on the topic page, where it resolves the open puzzle in the archived derivation this function descends from. Parameters ---------- s : ndarray The waveform at the reference sensor. sigma2_ref, sigma2_del : float Noise variance per sample at the reference and delayed sensor. fs : float Sample rate in Hz. gain : float Amplitude gain of the second sensor's copy (attenuation < 1). Returns ------- float Lower bound on the delay-difference variance, seconds squared. """ if gain == 0: raise ValueError("gain must be nonzero (no signal, no bound)") a_inv = crlb_delay(s, sigma2_ref, fs) b_inv = crlb_delay(gain * np.asarray(s, dtype=float), sigma2_del, fs) return float(a_inv + b_inv) def signal_centroid(x: np.ndarray, fs: float = 1.0) -> float: """Gabor centroid of a signal's energy: tau_c = m1 / m0. tau_c = sum_k t_k |x[k]|^2 / sum_k |x[k]|^2, t_k = k / fs. The energy-weighted mean arrival time, measured from the first sample. Two properties make it a good reference point for timing a pulse: it is shift-equivariant (delaying the signal by d moves the centroid by exactly d, so centroid differences measure delays), and in the frequency domain it equals the energy-weighted mean **group delay**, which stays well defined even when a dispersive channel has smeared the pulse so badly that "the peak" and "the onset" have become arbitrary choices. That identity, and the partial-discharge location method built on it, is derived on the topic page. Parameters ---------- x : ndarray Signal record (real or complex); energy must be nonzero and should be well contained in the record. fs : float Sample rate in Hz. Returns ------- float Centroid time in seconds from the start of the record. """ x = np.asarray(x) if x.ndim != 1 or len(x) < 1: raise ValueError("x must be a one-dimensional record") if fs <= 0: raise ValueError("fs must be strictly positive") e = np.abs(x) ** 2 m0 = e.sum() if m0 == 0: raise ValueError("signal has zero energy") t = np.arange(len(x)) / fs return float(np.sum(t * e) / m0) def centroid_delay(x: np.ndarray, y: np.ndarray, fs: float = 1.0) -> float: """Delay of y relative to x as the difference of energy centroids. For a non-dispersive channel this equals the true delay (both records aligned to the same clock); for a dispersive channel it equals the energy-weighted mean group delay of the channel, the quantity the centroid method of the topic page localises partial-discharge sources with. Unlike the correlation peak it needs no waveform model at all, but it inherits the centroid's sensitivity to noise far from the pulse: window the records to the pulse's neighbourhood first. Parameters ---------- x, y : ndarray The two records, sampled on a common clock. fs : float Sample rate in Hz. Returns ------- float Centroid delay in seconds. """ return signal_centroid(y, fs) - signal_centroid(x, fs)