diff --git a/AGENTS.md b/AGENTS.md index 43bb42e9..1cd6f62e 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -16,6 +16,8 @@ Marc uses AI to help create the JavaScript mini-apps and solve issues like when ## How to build locally +If the prompt does not ask to build it, then don't build it. + - Activate the project virtual environment which should be in the root of this repo under .venv - Build the site using: @@ -50,8 +52,11 @@ PySDR's prose is intentionally instructional, conversational, and example-driven - Keep technical terminology precise, but avoid sounding overly academic or formal. - When a section already has figures or code, make the surrounding prose point the reader to them and explain what they should notice. - Avoid hype, filler, and motivational fluff. - -When in doubt, read the surrounding chapter aloud mentally and aim for the same cadence and readability rather than a generic documentation voice. +- Lead with a concrete scenario before any equation. Pose a small "what if" with real numbers (e.g. "the emitter is 100 m closer to one sensor"), then generalize. Introduce the named concept (hyperbola, foci) in plain words before showing the formula. +- Read the equation back in plain English. Right after a .. math:: block, add a sentence translating it ("which reads: distance to one sensor minus distance to the other equals..."). +- Explain why, not just what. For each fact, give the intuition behind it (why a range difference can't exceed the baseline) rather than stating it as a rule. +- Define jargon inline the moment it appears ("the baseline," "ill-conditioned, meaning small errors move the estimate a lot") instead of assuming the reader knows it. +- Use second person and rhetorical questions to walk the reader through the reasoning as if thinking aloud. ## Notes diff --git a/_images/tdoa_cramer_rao.svg b/_images/tdoa_cramer_rao.svg new file mode 100644 index 00000000..90e2e616 --- /dev/null +++ b/_images/tdoa_cramer_rao.svg @@ -0,0 +1,1468 @@ + + + + + + + + 2026-06-25T12:54:18.794956 + image/svg+xml + + + Matplotlib v3.10.9, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/_images/tdoa_gdop.svg b/_images/tdoa_gdop.svg new file mode 100644 index 00000000..eb79e378 --- /dev/null +++ b/_images/tdoa_gdop.svg @@ -0,0 +1,1634 @@ + + + + + + + + 2026-06-23T03:53:09.404485 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + −8 + + + + + + + + + + −6 + + + + + + + + + + −4 + + + + + + + + + + −2 + + + + + + + + + + 0 + + + + + + + + + + 2 + + + + + + + + + + 4 + + + + + + + + + + 6 + + + + + + + + + + 8 + + + + x + + + + + + + + + + + + + + −8 + + + + + + + + + + −6 + + + + + + + + + + −4 + + + + + + + + + + −2 + + + + + + + + + + 0 + + + + + + + + + + 2 + + + + + + + + + + 4 + + + + + + + + + + 6 + + + + + + + + + + 8 + + + + y + + + + + + + + + + + + + + + + + + + + + S1 + + + S2 + + + S3 + + + triangular array + + + + 2 + + + + + 3 + + + + + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + −8 + + + + + + + + + + −6 + + + + + + + + + + −4 + + + + + + + + + + −2 + + + + + + + + + + 0 + + + + + + + + + + 2 + + + + + + + + + + 4 + + + + + + + + + + 6 + + + + + + + + + + 8 + + + + x + + + + + + + + + + + −8 + + + + + + + + + + −6 + + + + + + + + + + −4 + + + + + + + + + + −2 + + + + + + + + + + 0 + + + + + + + + + + 2 + + + + + + + + + + 4 + + + + + + + + + + 6 + + + + + + + + + + 8 + + + + y + + + + + + + + + + + + + + + + + + + + + S1 + + + S2 + + + S3 + + + nearly collinear array + + + + 2 + + + + + 2 + + + + + 3 + + + + + 3 + + + + + 5 + + + + + 5 + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.2 + + + + + + + + + + 2 + + + + + + + + + + 3 + + + + + + + + + + 5 + + + + + + + + + + ≥8 + + + + + + + + + + + + + + + 4 + + + + + + + + + + + + + 6 + + + + + + + + 7 + + + + + + + + + + + + GDOP (position error / range-difference error) + + + + + + + + + + + + + + + + diff --git a/_images/tdoa_hyperbola.svg b/_images/tdoa_hyperbola.svg new file mode 100644 index 00000000..282affa7 --- /dev/null +++ b/_images/tdoa_hyperbola.svg @@ -0,0 +1,768 @@ + + + + + + + + 2026-06-23T03:53:09.170870 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + −3 + + + + + + + + + + + + + −2 + + + + + + + + + + + + + −1 + + + + + + + + + + + + + 0 + + + + + + + + + + + + + 1 + + + + + + + + + + + + + 2 + + + + + + + + + + + + + 3 + + + + x + + + + + + + + + + + + + + + + + −3 + + + + + + + + + + + + + −2 + + + + + + + + + + + + + −1 + + + + + + + + + + + + + 0 + + + + + + + + + + + + + 1 + + + + + + + + + + + + + 2 + + + + + + + + + + + + + 3 + + + + y + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + S₁ + + + S₂ + + + baseline + + + Δr > 0 + (source nearer S₁) + + + Δr < 0 + (source nearer S₂) + + + Δr = 0 + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/_images/tdoa_principle.svg b/_images/tdoa_principle.svg new file mode 100644 index 00000000..a9bdac04 --- /dev/null +++ b/_images/tdoa_principle.svg @@ -0,0 +1,619 @@ + + + + + + + + 2026-06-23T03:58:57.663473 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + wavefront + expands at + speed c + + + Source + (unknown position, + unknown transmit time t₀) + + + + + + r₁ + + + S₁ + + + + + + r₂ + + + S₂ + + + + + + r₃ + + + S₃ + + + The TDOA principle: different path lengths give different arrival times + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + time + + + + + + + t₀ (unknown + transmit time) + + + S₁ + + + t₁ + + + S₂ + + + t₂ + + + S₃ + + + t₃ + + + + + + + + τ₂₁ = t₂ − t₁ + + + + + + + + τ₃₁ = t₃ − t₁ + + + + TDOA measures only the gaps τ between arrivals — the unknown transmit time t₀ cancels in every difference. + + + + + + + + + + + diff --git a/_images/tdoa_python_heatmap.svg b/_images/tdoa_python_heatmap.svg new file mode 100644 index 00000000..e8909cb8 --- /dev/null +++ b/_images/tdoa_python_heatmap.svg @@ -0,0 +1,1452 @@ + + + + + + + + 2026-06-25T13:28:25.517802 + image/svg+xml + + + Matplotlib v3.10.9, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/_images/tdoa_python_integer.svg b/_images/tdoa_python_integer.svg new file mode 100644 index 00000000..abdd2cf7 --- /dev/null +++ b/_images/tdoa_python_integer.svg @@ -0,0 +1,2078 @@ + + + + + + + + 2026-06-24T23:58:47.498877 + image/svg+xml + + + Matplotlib v3.10.9, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/_images/tdoa_python_subsample.svg b/_images/tdoa_python_subsample.svg new file mode 100644 index 00000000..ecd3ceed --- /dev/null +++ b/_images/tdoa_python_subsample.svg @@ -0,0 +1,2494 @@ + + + + + + + + 2026-06-24T23:58:47.790598 + image/svg+xml + + + Matplotlib v3.10.9, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/_static/js/tdoa.js b/_static/js/tdoa.js new file mode 100644 index 00000000..dcf15875 --- /dev/null +++ b/_static/js/tdoa.js @@ -0,0 +1,649 @@ +function tdoa_app(containerId) { + // ----- configuration ------------------------------------------------------- + const c = 3e8; // propagation speed [m/s] (free space / RF) + const W = 600; // canvas width [px] + const H = 480; // canvas height [px] + const worldSpan = 1000; // world width represented across the canvas [m] + const nodeRadius = 9; // hit/draw radius for draggable handles [px] + const edgeMargin = 12; // keep handles at least this far inside the canvas [px] + const maxSensors = 10; // upper bound on how many sensors the user can add + const palette = ["#e6550d", "#3182bd", "#31a354"]; // first few sensor-pair colors + + // heatmap: give each hyperbola some width and add them together so their + // overlap lights up where the emitter actually is + let heatHalfWidth = 150; // band half-width around each hyperbola [m] + const heatRes = 2; // pixel block size of the heatmap grid (quality vs speed) + const heatMaxAlpha = 0.55; // overlay opacity where every band coincides + + // Color for the k-th sensor pair: use the fixed palette first, then spread the + // remaining hues around the color wheel so every pair stays distinguishable. + function pairColor(k) { + if (k < palette.length) return palette[k]; + return `hsl(${(k * 47) % 360}, 65%, 45%)`; + } + + // ----- DOM setup ----------------------------------------------------------- + const container = document.getElementById(containerId || "tdoaApp") || document.body; + + // canvas sits in a flex row next to the noise slider on its right + const row = document.createElement("div"); + row.style.display = "flex"; + row.style.alignItems = "stretch"; + row.style.gap = "10px"; + container.appendChild(row); + + // wrapper lets us overlay controls (the Add-sensor button) on top of the canvas + const canvasWrap = document.createElement("div"); + canvasWrap.style.position = "relative"; + canvasWrap.style.lineHeight = "0"; // avoid extra space under the canvas + row.appendChild(canvasWrap); + + const canvas = document.createElement("canvas"); + canvas.width = W; + canvas.height = H; + canvas.style.border = "1px solid #888"; + canvas.style.touchAction = "none"; // let us handle touch-drag ourselves + canvas.style.cursor = "grab"; + canvasWrap.appendChild(canvas); + + // vertical noise slider: standard deviation of the Gaussian noise [m] + const sliderBox = document.createElement("div"); + sliderBox.style.display = "flex"; + sliderBox.style.flexDirection = "column"; + sliderBox.style.alignItems = "center"; + sliderBox.style.fontFamily = "sans-serif"; + sliderBox.style.fontSize = "12px"; + row.appendChild(sliderBox); + + const sliderLabel = document.createElement("div"); + sliderLabel.style.textAlign = "center"; + sliderLabel.style.marginBottom = "6px"; + sliderBox.appendChild(sliderLabel); + + let noiseStd = 0; // std dev of Gaussian noise added to each range diff [m] + function updateSliderLabel() { + sliderLabel.innerHTML = `Noise
${noiseStd.toFixed(0)} m`; + } + updateSliderLabel(); + + const slider = document.createElement("input"); + slider.type = "range"; + slider.min = "0"; + slider.max = "100"; + slider.step = "1"; + slider.value = "0"; + // make the range input vertical (with a fallback for older browsers) + slider.setAttribute("orient", "vertical"); + slider.style.writingMode = "vertical-lr"; + slider.style.direction = "rtl"; // 0 at the bottom, max at the top + slider.style.height = H / 2 + "px"; // half the canvas height + slider.style.width = "24px"; + sliderBox.appendChild(slider); + + slider.addEventListener("input", () => { + noiseStd = parseFloat(slider.value); + updateSliderLabel(); + render(); + }); + + // vertical dynamic-range slider: a gamma exponent applied to the heatmap + // intensity. Higher values darken the weak single bands and let only the + // strong overlaps near the emitter stand out; lower values flatten it out. + let heatGamma = 1.5; // exponent applied to the normalized heatmap intensity + const drBox = document.createElement("div"); + drBox.style.display = "flex"; + drBox.style.flexDirection = "column"; + drBox.style.alignItems = "center"; + drBox.style.fontFamily = "sans-serif"; + drBox.style.fontSize = "12px"; + row.appendChild(drBox); + + const drLabel = document.createElement("div"); + drLabel.style.textAlign = "center"; + drLabel.style.marginBottom = "6px"; + drBox.appendChild(drLabel); + + function updateDrLabel() { + drLabel.innerHTML = `Dyn.
range
${heatGamma.toFixed(1)}`; + } + updateDrLabel(); + + const drSlider = document.createElement("input"); + drSlider.type = "range"; + drSlider.min = "0.5"; + drSlider.max = "5"; + drSlider.step = "0.1"; + drSlider.value = String(heatGamma); + drSlider.setAttribute("orient", "vertical"); + drSlider.style.writingMode = "vertical-lr"; + drSlider.style.direction = "rtl"; // low at the bottom, high at the top + drSlider.style.height = H / 2 + "px"; + drSlider.style.width = "24px"; + drBox.appendChild(drSlider); + + drSlider.addEventListener("input", () => { + heatGamma = parseFloat(drSlider.value); + updateDrLabel(); + render(false); // visual-only: keep the existing noise samples + }); + + // vertical width slider: half-width of the band drawn around each hyperbola. + // Wider bands overlap more readily (good with lots of noise); narrow bands + // pin the emitter down tightly. + const widthBox = document.createElement("div"); + widthBox.style.display = "flex"; + widthBox.style.flexDirection = "column"; + widthBox.style.alignItems = "center"; + widthBox.style.fontFamily = "sans-serif"; + widthBox.style.fontSize = "12px"; + row.appendChild(widthBox); + + const widthLabel = document.createElement("div"); + widthLabel.style.textAlign = "center"; + widthLabel.style.marginBottom = "6px"; + widthBox.appendChild(widthLabel); + + function updateWidthLabel() { + widthLabel.innerHTML = `Width
${heatHalfWidth.toFixed(0)} m`; + } + updateWidthLabel(); + + const widthSlider = document.createElement("input"); + widthSlider.type = "range"; + widthSlider.min = "10"; + widthSlider.max = "200"; + widthSlider.step = "5"; + widthSlider.value = String(heatHalfWidth); + widthSlider.setAttribute("orient", "vertical"); + widthSlider.style.writingMode = "vertical-lr"; + widthSlider.style.direction = "rtl"; // narrow at the bottom, wide at the top + widthSlider.style.height = H / 2 + "px"; + widthSlider.style.width = "24px"; + widthBox.appendChild(widthSlider); + + widthSlider.addEventListener("input", () => { + heatHalfWidth = parseFloat(widthSlider.value); + updateWidthLabel(); + render(false); // visual-only: keep the existing noise samples + }); + + // buttons to add/remove a sensor, overlaid on the top-left corner of the canvas + const addBtn = document.createElement("button"); + addBtn.style.position = "absolute"; + addBtn.style.top = "8px"; + addBtn.style.left = "8px"; + addBtn.style.fontFamily = "sans-serif"; + addBtn.style.fontSize = "13px"; + addBtn.style.lineHeight = "normal"; + addBtn.style.padding = "4px 10px"; + addBtn.style.cursor = "pointer"; + canvasWrap.appendChild(addBtn); + + const removeBtn = document.createElement("button"); + removeBtn.style.position = "absolute"; + removeBtn.style.top = "40px"; + removeBtn.style.left = "8px"; + removeBtn.style.fontFamily = "sans-serif"; + removeBtn.style.fontSize = "13px"; + removeBtn.style.lineHeight = "normal"; + removeBtn.style.padding = "4px 10px"; + removeBtn.style.cursor = "pointer"; + canvasWrap.appendChild(removeBtn); + + // checkbox to toggle the heatmap overlay, overlaid below the buttons + let showHeatmap = true; // heatmap is on by default + const heatToggle = document.createElement("label"); + heatToggle.style.position = "absolute"; + heatToggle.style.top = "72px"; + heatToggle.style.left = "8px"; + heatToggle.style.display = "flex"; + heatToggle.style.alignItems = "center"; + heatToggle.style.gap = "4px"; + heatToggle.style.fontFamily = "sans-serif"; + heatToggle.style.fontSize = "13px"; + heatToggle.style.color = "#222"; + heatToggle.style.cursor = "pointer"; + heatToggle.style.userSelect = "none"; + + const heatCheckbox = document.createElement("input"); + heatCheckbox.type = "checkbox"; + heatCheckbox.checked = showHeatmap; + heatCheckbox.style.cursor = "pointer"; + heatCheckbox.style.margin = "0"; + heatToggle.appendChild(heatCheckbox); + heatToggle.appendChild(document.createTextNode("Heatmap")); + canvasWrap.appendChild(heatToggle); + + heatCheckbox.addEventListener("change", () => { + showHeatmap = heatCheckbox.checked; + render(false); // visual-only: keep the existing noise samples + }); + + // collapsible details panel: a thick bar you click to reveal the text readout, + // collapsed by default to keep the figure compact + const details = document.createElement("details"); + details.style.marginTop = "8px"; + details.style.width = W + "px"; + details.style.border = "1px solid #ccc"; + details.style.borderRadius = "4px"; + details.style.overflow = "hidden"; + container.appendChild(details); + + const summary = document.createElement("summary"); + summary.textContent = "Show Debug Info"; + summary.style.cursor = "pointer"; + summary.style.userSelect = "none"; + summary.style.fontFamily = "sans-serif"; + summary.style.fontSize = "13px"; + summary.style.fontWeight = "bold"; + summary.style.padding = "10px 12px"; + summary.style.background = "#f0f0f0"; + summary.style.color = "#333"; + details.appendChild(summary); + + // swap the label between collapsed/expanded states + details.addEventListener("toggle", () => { + summary.textContent = details.open ? "Hide Debug Info" : "Show Debug Info"; + }); + + const readout = document.createElement("div"); + readout.style.fontFamily = "monospace"; + readout.style.fontSize = "13px"; + readout.style.padding = "8px 12px"; + details.appendChild(readout); + + const ctx = canvas.getContext("2d"); + + // ----- scene state (world coordinates, meters, origin at center) ----------- + // y points up in world coordinates (flipped when drawing to the canvas). + const emitter = { x: -155, y: 226, label: "Emitter" }; + const sensors = [ + { x: -350, y: -200, label: "Sensor 0" }, + { x: 350, y: -200, label: "Sensor 1" }, + { x: 0, y: 300, label: "Sensor 2" } + ]; + + // ----- coordinate transforms ---------------------------------------------- + const scale = W / worldSpan; // px per meter + function worldToPx(p) { + return { x: W / 2 + p.x * scale, y: H / 2 - p.y * scale }; + } + function pxToWorld(px) { + return { x: (px.x - W / 2) / scale, y: (H / 2 - px.y) / scale }; + } + + function dist(a, b) { + return Math.hypot(a.x - b.x, a.y - b.y); + } + + // ----- TDOA simulation ----------------------------------------------------- + // Standard normal sample via the Box-Muller transform. + function randn() { + let u = 0; + let v = 0; + while (u === 0) u = Math.random(); + while (v === 0) v = Math.random(); + return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v); + } + + // Every unique sensor pair (i < j); order matches the original [0,1],[0,2],[1,2]. + function sensorPairs() { + const pairs = []; + for (let i = 0; i < sensors.length; i++) { + for (let j = i + 1; j < sensors.length; j++) pairs.push([i, j]); + } + return pairs; + } + + // Returns the per-sensor TOA and, for each pair, the TDOA and range diff. + function simulate() { + const toa = sensors.map((s) => dist(emitter, s) / c); // seconds + const measurements = sensorPairs().map(([i, j]) => { + // ideal TDOA, then add Gaussian noise (slider sets its std dev in meters, + // so we convert that range error into the equivalent time error) + const tdoa = toa[i] - toa[j] + (noiseStd * randn()) / c; // seconds + return { i, j, tdoa, dr: c * tdoa }; // dr = r_i - r_j [m] + }); + return { toa, measurements }; + } + + // ----- hyperbola drawing --------------------------------------------------- + // Locus of points u with |u - s_i| - |u - s_j| = dr is a hyperbola with foci + // at the two sensors. We parametrize it in a frame centered on the midpoint + // of the foci, with the transverse axis along the baseline. + function drawHyperbola(si, sj, dr, color) { + const mid = { x: (si.x + sj.x) / 2, y: (si.y + sj.y) / 2 }; + const baseline = dist(si, sj); + const cFoci = baseline / 2; // half the focal separation + const a = dr / 2; // signed semi-transverse axis; sign picks the branch + if (Math.abs(a) >= cFoci) return; // |dr| can't exceed the baseline + const b = Math.sqrt(cFoci * cFoci - a * a); + + // Unit vector u from s_i toward s_j (axis), and perpendicular v. + const ux = (sj.x - si.x) / baseline; + const uy = (sj.y - si.y) / baseline; + const vx = -uy; + const vy = ux; + + // Sweep the parameter t; x = a*cosh(t) keeps us on the correct branch + // because a carries the sign of dr. + ctx.beginPath(); + let first = true; + for (let t = -3; t <= 3.0001; t += 0.05) { + const xl = a * Math.cosh(t); + const yl = b * Math.sinh(t); + const wx = mid.x + xl * ux + yl * vx; + const wy = mid.y + xl * uy + yl * vy; + const p = worldToPx({ x: wx, y: wy }); + if (first) { + ctx.moveTo(p.x, p.y); + first = false; + } else { + ctx.lineTo(p.x, p.y); + } + } + ctx.strokeStyle = color; + ctx.lineWidth = 2; + ctx.stroke(); + } + + // ----- heatmap ------------------------------------------------------------- + // Reuse the same hyperbolas, but instead of an infinitely thin curve give + // each one a band of width 2*heatHalfWidth whose intensity follows a raised + // sine (raised cosine): 1 right on the hyperbola, tapering to 0 at the band + // edges. Summing every band makes their common crossing — the emitter — the + // brightest spot. We compute on a coarse grid and let drawImage smooth it. + const heatCanvas = document.createElement("canvas"); + const heatCtx = heatCanvas.getContext("2d"); + + // warm colormap: yellow at low intensity ramping to red at high intensity + function heatColor(t) { + return [255, Math.round(220 * (1 - t)), Math.round(40 * (1 - t))]; + } + + function drawHeatmap(measurements) { + if (measurements.length === 0) return; + const gw = Math.ceil(W / heatRes); + const gh = Math.ceil(H / heatRes); + heatCanvas.width = gw; + heatCanvas.height = gh; + const img = heatCtx.createImageData(gw, gh); + const data = img.data; + const nPairs = measurements.length; + + for (let gy = 0; gy < gh; gy++) { + for (let gx = 0; gx < gw; gx++) { + // world coordinates at the center of this grid block + const w = pxToWorld({ + x: gx * heatRes + heatRes / 2, + y: gy * heatRes + heatRes / 2 + }); + + let sum = 0; + for (let k = 0; k < nPairs; k++) { + const m = measurements[k]; + const si = sensors[m.i]; + const sj = sensors[m.j]; + const ri = Math.hypot(w.x - si.x, w.y - si.y); + const rj = Math.hypot(w.x - sj.x, w.y - sj.y); + // how far this point's range difference is from the measured one; + // 0 means the point sits exactly on pair k's hyperbola + const residual = ri - rj - m.dr; + if (Math.abs(residual) < heatHalfWidth) { + sum += 0.5 * (1 + Math.cos((Math.PI * residual) / heatHalfWidth)); + } + } + + // normalized intensity (1 where all bands overlap, the emitter), then + // a gamma to set the dynamic range: >1 suppresses the weak bands + const norm = Math.pow(sum / nPairs, heatGamma); + const idx = (gy * gw + gx) * 4; + const rgb = heatColor(norm); + data[idx] = rgb[0]; + data[idx + 1] = rgb[1]; + data[idx + 2] = rgb[2]; + data[idx + 3] = Math.round(norm * heatMaxAlpha * 255); + } + } + heatCtx.putImageData(img, 0, 0); + + // scale the coarse buffer up to full size; bilinear smoothing hides the grid + ctx.imageSmoothingEnabled = true; + ctx.drawImage(heatCanvas, 0, 0, gw, gh, 0, 0, W, H); + } + + // ----- rendering ----------------------------------------------------------- + function drawGrid() { + ctx.clearRect(0, 0, W, H); + ctx.strokeStyle = "#eee"; + ctx.lineWidth = 1; + const step = 100 * scale; // grid every 100 m + for (let x = (W / 2) % step; x < W; x += step) { + ctx.beginPath(); + ctx.moveTo(x, 0); + ctx.lineTo(x, H); + ctx.stroke(); + } + for (let y = (H / 2) % step; y < H; y += step) { + ctx.beginPath(); + ctx.moveTo(0, y); + ctx.lineTo(W, y); + ctx.stroke(); + } + // axes + ctx.strokeStyle = "#ccc"; + ctx.beginPath(); + ctx.moveTo(W / 2, 0); + ctx.lineTo(W / 2, H); + ctx.moveTo(0, H / 2); + ctx.lineTo(W, H / 2); + ctx.stroke(); + + // axis labels + ctx.fillStyle = "#999"; + ctx.font = "13px sans-serif"; + ctx.textBaseline = "alphabetic"; + ctx.textAlign = "right"; + ctx.fillText("X", W - 6, H / 2 - 6); + ctx.textAlign = "left"; + ctx.fillText("Y", W / 2 + 6, 14); + ctx.textAlign = "left"; + } + + function drawSensor(s) { + const p = worldToPx(s); + ctx.fillStyle = "#222"; + ctx.beginPath(); + ctx.moveTo(p.x, p.y - nodeRadius); + ctx.lineTo(p.x + nodeRadius, p.y + nodeRadius); + ctx.lineTo(p.x - nodeRadius, p.y + nodeRadius); + ctx.closePath(); + ctx.fill(); + ctx.fillStyle = "#222"; + ctx.font = "bold 16px sans-serif"; + ctx.fillText(s.label, p.x + nodeRadius + 2, p.y - 2); + } + + // small label near a node showing its world coordinates in meters + function drawCoordTooltip(node) { + const p = worldToPx(node); + const text = `(${node.x.toFixed(0)}, ${node.y.toFixed(0)}) m`; + ctx.font = "12px monospace"; + ctx.textAlign = "left"; + ctx.textBaseline = "alphabetic"; + const padX = 5; + const padY = 3; + const tw = ctx.measureText(text).width; + const boxW = tw + padX * 2; + const boxH = 16 + padY * 2; + // sit the box just below the node, nudged on-screen if it would clip + let bx = p.x + nodeRadius + 2; + let by = p.y + nodeRadius + 2; + if (bx + boxW > W) bx = W - boxW; + if (by + boxH > H) by = p.y - nodeRadius - boxH - 2; + ctx.fillStyle = "rgba(0, 0, 0, 0.8)"; + ctx.fillRect(bx, by, boxW, boxH); + ctx.fillStyle = "#fff"; + ctx.fillText(text, bx + padX, by + boxH - padY - 3); + } + + function drawEmitter() { + const p = worldToPx(emitter); + // black outline, no fill, so the heatmap underneath stays visible + ctx.beginPath(); + ctx.arc(p.x, p.y, nodeRadius, 0, 2 * Math.PI); + ctx.strokeStyle = "#000"; + ctx.lineWidth = 2; + ctx.stroke(); + ctx.fillStyle = "#000"; + ctx.font = "bold 16px sans-serif"; + ctx.fillText(emitter.label, p.x + nodeRadius + 2, p.y - 2); + } + + // cache the last simulation so purely-visual changes (heatmap width, dynamic + // range, hover) can repaint without drawing fresh noise samples + let lastSim = null; + function render(resimulate = true) { + if (resimulate || !lastSim) lastSim = simulate(); + const { toa, measurements } = lastSim; + + drawGrid(); + if (showHeatmap) drawHeatmap(measurements); + measurements.forEach((m, k) => { + drawHyperbola(sensors[m.i], sensors[m.j], m.dr, pairColor(k)); + }); + sensors.forEach(drawSensor); + drawEmitter(); + + // coordinate tooltip for the node under the cursor (or being dragged) + const activeNode = dragTarget || hoverNode; + if (activeNode) drawCoordTooltip(activeNode); + + // text readout of the simulated quantities + let html = ""; + toa.forEach((t, i) => { + html += `TOA(${sensors[i].label}) = ${(t * 1e9).toFixed(1)} ns   (range ${dist(emitter, sensors[i]).toFixed(0)} m)
`; + }); + measurements.forEach((m, k) => { + html += `TDOA(${sensors[m.i].label},${sensors[m.j].label}) = ${(m.tdoa * 1e9).toFixed(1)} ns   Δr = ${m.dr.toFixed(0)} m
`; + }); + readout.innerHTML = html; + } + + // ----- dragging ------------------------------------------------------------ + let dragTarget = null; + let hoverNode = null; // node currently under the cursor (for the tooltip) + + function eventPx(e) { + const rect = canvas.getBoundingClientRect(); + const src = e.touches ? e.touches[0] : e; + return { x: src.clientX - rect.left, y: src.clientY - rect.top }; + } + + function pickNode(px) { + const all = [emitter, ...sensors]; + for (const node of all) { + if (dist(worldToPx(node), px) <= nodeRadius + 4) return node; + } + return null; + } + + function onDown(e) { + const px = eventPx(e); + dragTarget = pickNode(px); + if (dragTarget) { + canvas.style.cursor = "grabbing"; + e.preventDefault(); + } + } + + function onMove(e) { + const px = eventPx(e); + if (!dragTarget) { + const hit = pickNode(px); + canvas.style.cursor = hit ? "grab" : "default"; + // only repaint when the hovered node actually changes, so the tooltip + // appears/disappears without re-noising the scene on every mouse move + if (hit !== hoverNode) { + hoverNode = hit; + render(false); // visual-only: keep the existing noise samples + } + return; + } + // keep the handle inside the canvas (with a small margin) so it can't be + // dragged off-screen and lost + px.x = Math.max(edgeMargin, Math.min(W - edgeMargin, px.x)); + px.y = Math.max(edgeMargin, Math.min(H - edgeMargin, px.y)); + const w = pxToWorld(px); + dragTarget.x = w.x; + dragTarget.y = w.y; + render(); + e.preventDefault(); + } + + function onUp() { + dragTarget = null; + canvas.style.cursor = "grab"; + } + + canvas.addEventListener("mouseleave", () => { + if (hoverNode) { + hoverNode = null; + render(false); // visual-only: keep the existing noise samples + } + }); + + canvas.addEventListener("mousedown", onDown); + window.addEventListener("mousemove", onMove); + window.addEventListener("mouseup", onUp); + canvas.addEventListener("touchstart", onDown, { passive: false }); + canvas.addEventListener("touchmove", onMove, { passive: false }); + canvas.addEventListener("touchend", onUp); + + // ----- adding / removing sensors ------------------------------------------- + const minSensors = 2; // at least 2 sensors to form one TDOA pair + + function updateSensorButtons() { + const atMax = sensors.length >= maxSensors; + addBtn.disabled = atMax; + addBtn.textContent = atMax + ? `Max sensors reached (${maxSensors})` + : `Add sensor`; + + const atMin = sensors.length <= minSensors; + removeBtn.disabled = atMin; + removeBtn.textContent = atMin + ? `Min sensors reached (${minSensors})` + : "Remove sensor"; + } + + function addSensor() { + if (sensors.length >= maxSensors) return; + const idx = sensors.length; + // place the new sensor on a ring, stepping by the golden angle so successive + // sensors spread out rather than landing on top of each other + const ang = idx * 2.399963229728653; + const r = 320; + sensors.push({ + x: r * Math.cos(ang), + y: r * Math.sin(ang), + label: "Sensor " + idx + }); + updateSensorButtons(); + render(); + } + + function removeSensor() { + if (sensors.length <= minSensors) return; + sensors.pop(); // drop the most recently added sensor + updateSensorButtons(); + render(); + } + + addBtn.addEventListener("click", addSensor); + removeBtn.addEventListener("click", removeSensor); + updateSensorButtons(); + + // ----- go ------------------------------------------------------------------ + render(); +} diff --git a/conf.py b/conf.py index 08c8b1b1..5750e653 100644 --- a/conf.py +++ b/conf.py @@ -201,7 +201,8 @@ def setup(app): 'js/beamforming_slider_app.js', 'js/FFT.js', 'js/cyclostationary_app.js', - 'js/homepage_app.js' + 'js/homepage_app.js', + 'js/tdoa.js' # we also include the index.js file from the PhasedArrayVisualizer directory in setup() above ] diff --git a/content/tdoa.rst b/content/tdoa.rst new file mode 100644 index 00000000..9145667d --- /dev/null +++ b/content/tdoa.rst @@ -0,0 +1,753 @@ +.. _tdoa-chapter: + +#### +TDOA +#### + +Time Difference of Arrival (TDOA) is a technique that can find the position of a transmitter (a.k.a. emitter) using multiple synchronized receivers (a.k.a. sensors), by comparing differences in signal arrival time. This chapter covers the full TDOA pipeline: geometry, GCC-PHAT time-delay estimation, closed-form and maximum-likelihood localization, accuracy bounds (CRLB and GDOP), and challenges like synchronization and multipath. TDOA is commonly used in both RF and acoustic/sonar applications. + +Before diving in, try playing with the interactive demo below to get a quick feel for how TDOA works, which involves the intersection of hyperbolas. + +.. raw:: html + +
+ + +************ +Introduction +************ + +A common problem across RF and acoustics/sonar is the desire to find the position of an emitter, also known as the process of geolocation. The emitter may be cooperative (a cell phone trying to be found) or non-cooperative (a radar emitter that would rather not be), stationary or moving, and the medium may be air, water, or free space. TDOA-based localization appears in cellular emergency-caller location, acoustics with microphone arrays (e.g., gunshot-detection systems mounted on city streetlights), passive sonar, passive (non-emitting) radar, electronic warfare and signals intelligence, and even wildlife tracking. In each case the engineering details differ, but the mathematical skeleton is the same. + +The key behind TDOA is that when the same wavefront hits two sensors, the difference in arrival times depends only on geometry, not on when the emitter transmitted. To see why, consider the propagation time from an emitter to sensor :math:`i`: :math:`t_i = t_0 + r_i / c`, where :math:`t_0` is the (unknown) start of transmission, :math:`r_i` is the emitter-to-sensor distance, and :math:`c` is the propagation speed. If we subtract the arrival times at two sensors, + +.. math:: + + \tau_{ij} = t_i - t_j = \frac{r_i - r_j}{c}, + +the unknown :math:`t_0` vanishes, which is good because we will likely never know :math:`t_0`. The TDOA depends only on the *difference* of ranges, which depends only on emitter and sensor geometry. This single fact is why TDOA dominates for non-cooperative emitters; we never need to know when the emitter transmitted, only that the same wavefront reached our synchronized receivers at measurable relative delays. You still need to isolate the signal so that you're only observing one emitter, so signal detection and classification may be necessary, plus filtering. + +Each pair of sensors yields one TDOA, and each TDOA traces out one hyperbola, so the number of hyperbolas we can draw is just the number of sensor pairs. With :math:`N` sensors that is + +.. math:: + + \binom{N}{2} = \frac{N(N-1)}{2}, + +i.e. 3 sensors give 3 hyperbolas, 4 give 6, 5 give 10, and so on. Not all of these are independent, as we will see below, only :math:`N-1` carry new geometric information, but the full set is still useful for averaging out noise. + +The price we pay is that the receivers must share a precise common time reference, a requirement that, as discussed later, is itself a demanding engineering problem because a timing error of one nanosecond corresponds to about 1 foot or 0.3 meters of range error, at least in RF applications. + +************* +TDOA Geometry +************* + +From Time Difference to Range Difference +=============================================== + +Multiplying a measured time difference by the propagation speed converts it into a *range difference*: + +.. math:: + + \Delta r_{ij} = c\,\tau_{ij} = r_i - r_j . + +For acoustic problems :math:`c \approx 343` m/s in air; for radio problems :math:`c \approx 2.998\times10^8` m/s. Note immediately the consequence for accuracy: in air, a :math:`0.1` ms timing error is only :math:`\sim`\3 cm, whereas in free space the same timing error is 30 km. Radio TDOA therefore demands extraordinarily precise timing, a theme we return to repeatedly. + +The diagram below shows an example of an emitter and three sensors, with a time domain plot of the signal being received by each sensor at different times. + +.. image:: ../_images/tdoa_principle.svg + :align: center + :target: ../_images/tdoa_principle.svg + :alt: An emitter and three sensors, with a time domain plot of the signal being received by each sensor at different times. + +The Hyperbola +=================== + +Now think about what a single range difference actually tells us. Suppose we have two sensors, and we have measured that the emitter is, say, 100 meters closer to one sensor than the other. Where could the emitter be? Not at a single spot, it turns out, but anywhere along a curved line. As you slide along that line, the two distances to the sensors both change, but their *difference* stays fixed at 100 meters the whole way. + +That curve has a name: it is a **hyperbola**, with the two sensors sitting at its two focal points. (In 3D the same idea sweeps out a curved surface called a hyperboloid, but the 2D picture is easier to reason about and everything carries over). If we write the emitter position as :math:`\mathbf{u}` and the two sensor positions as :math:`\mathbf{s}_i` and :math:`\mathbf{s}_j`, the hyperbola is just the set of points obeying + +.. math:: + + |\mathbf{u}-\mathbf{s}_i| - |\mathbf{u}-\mathbf{s}_j| = \Delta r_{ij} = \text{constant}, + +which reads "distance to one sensor minus distance to the other equals our measured range difference". A few practical consequences fall right out of this: + +* **A range difference can't exceed the spacing between the two sensors.** Intuitively, the difference of two distances is largest when the emitter lies directly out beyond one sensor along the line connecting them, and even then it can only equal that sensor-to-sensor spacing (often called the *baseline*). So if you ever measure a range difference bigger than the baseline, something is wrong, most likely noise, multipath, or a timing/synchronization error. +* **The sign tells you which side you're on.** A hyperbola actually has two mirror-image branches, one curving toward each sensor. Whether your range difference came out positive or negative picks the branch nearer the closer sensor, so you don't get confused between the two halves. +* **The shape depends on the measurement.** When the range difference is close to the full baseline, the hyperbola hugs the line between the sensors. When the range difference is near zero (the emitter is roughly equidistant), the curve straightens out into the line that perpendicularly bisects the baseline. Near both of these extremes the geometry becomes "ill-conditioned," meaning small measurement errors push the estimated position around a lot, so positioning there is less reliable. + +A single TDOA thus constrains the emitter's position to a curve, not a point. To fix a position we intersect several such curves. Below we plot two sensors, and several hyperbola branches drawn for :math:`\Delta r < 0`, :math:`\Delta r = 0` (the perpendicular bisector), and :math:`\Delta r > 0`. On each hyperbola, the TDOA between the two sensors is constant. If you calculated the TDOA with just two sensors, you would know it is somewhere on that line, but you would need a third sensor to find where on that line it is (performing geolocation). + +.. image:: ../_images/tdoa_hyperbola.svg + :align: center + :target: ../_images/tdoa_hyperbola.svg + :alt: Two sensors, and several hyperbola branches drawn + +Multilateration +===================== + +With :math:`N` sensors we can form pairs and intersect their hyperbolas; the emitter lies at (or near) their common intersection. This process is **hyperbolic multilateration**. Counting degrees of freedom tells us how many sensors we need: + +* In **2D** the emitter position has 2 unknowns :math:`(x,y)`. Each independent TDOA gives one equation, so we need at least 2 independent TDOAs, which requires **3 sensors**. For example, if we know the emitter is on land and we're not having to take into account curvature of the Earth, this would work. +* In **3D** the emitter position has 3 unknowns :math:`(x,y,z)`, requiring 3 independent TDOAs and therefore **4 sensors**. + +In the noiseless case the hyperbolas meet at a single point (with an occasional geometric ambiguity resolved by branch signs or an extra sensor). With more sensors than the minimum the system is *overdetermined*: noisy hyperbolas no longer share an exact common point, and we must solve a least-squares or maximum-likelihood problem, as described below. + +Reference Sensor and Independent Pairs +============================================= + +From :math:`N` sensors one can form :math:`\binom{N}{2}` pairwise TDOAs, but they are not all independent. Choosing one sensor as a **reference** (say sensor 1) and forming :math:`\tau_{i1}` for :math:`i = 2,\dots,N` yields :math:`N-1` TDOAs from which every other pairwise difference can be reconstructed, since :math:`\tau_{ij} = \tau_{i1} - \tau_{j1}`. These :math:`N-1` are the *independent* measurements that carry all the geometric information. + +The redundant pairs are not worthless, however. Because each measured TDOA carries independent *noise*, using all :math:`\binom{N}{2}` pairs (with a correctly modeled, correlated noise covariance, where the reference sensor's noise is common to every :math:`\tau_{i1}`) can improve the estimate. + +Example: A Three-Sensor 2D Fix +============================== + +Place three sensors at + +.. math:: + + \mathbf{s}_1=(0,0),\quad \mathbf{s}_2=(100,0),\quad \mathbf{s}_3=(0,100)\ \text{(meters)}, + +and suppose the true emitter is at :math:`\mathbf{u}=(40,30)`. The emitter-to-sensor distances are + +.. math:: + + r_1=\sqrt{40^2+30^2}=50,\quad + r_2=\sqrt{60^2+30^2}=\sqrt{4500}\approx 67.08,\quad + r_3=\sqrt{40^2+70^2}=\sqrt{6500}\approx 80.62 . + +Taking sensor 1 as reference, the range-difference measurements are + +.. math:: + + \Delta r_{21}=r_2-r_1\approx 17.08\ \text{m},\qquad + \Delta r_{31}=r_3-r_1\approx 30.62\ \text{m}. + +Each defines a hyperbola with foci :math:`\{\mathbf{s}_2,\mathbf{s}_1\}` and :math:`\{\mathbf{s}_3,\mathbf{s}_1\}` respectively; their intersection is the emitter's position. Solving the two hyperbola equations by hand is awkward, which is precisely the motivation for the algebraic linearization developed below, where we will recover :math:`(40,30)` from exactly these numbers in closed form. + +************************************* +The Signal and Measurement Model +************************************* + +Received-Signal Model +============================ + +Each sensor receives a time-delayed, scaled, noisy copy of whatever the emitter is transmitting. Specifically, let :math:`s(t)` be the transmit waveform. Sensor :math:`i` receives: + +.. math:: + + x_i(t) = a_i \, s(t - t_i) + n_i(t), \qquad i = 1,\dots,N, + +where :math:`a_i` is a real (or complex, for passband signals) gain capturing propagation loss and antenna response, :math:`t_i = t_0 + r_i/c` is the absolute arrival time, and :math:`n_i(t)` is additive noise. This model assumes a single dominant line-of-sight path; multipath and non-line-of-sight effects are deferred to a later section. + +Defining the TDOA +======================== + +The pairwise TDOA is the difference of arrival times, + +.. math:: + + \tau_{ij} = t_i - t_j = \frac{r_i - r_j}{c} = \frac{|\mathbf{u}-\mathbf{s}_i| - |\mathbf{u}-\mathbf{s}_j|}{c}. + +The right-hand side makes explicit that the TDOA is a nonlinear function of the emitter coordinates :math:`\mathbf{u}`. The *measurement* problem is to estimate :math:`\tau_{ij}` from the waveforms :math:`x_i, x_j`; the *localization* problem is to invert the nonlinear map from :math:`\mathbf{u}` to the collection of TDOAs. + +Noise Assumptions +======================== + +We assume each :math:`n_i(t)` is zero-mean, wide-sense stationary, Gaussian, and independent of the transmitted signal and of the noise at other sensors. The per-sensor signal-to-noise ratio is + +.. math:: + + \mathrm{SNR}_i = \frac{a_i^2 \sigma_s^2}{\sigma_{n_i}^2}, + +with :math:`\sigma_s^2` and :math:`\sigma_{n_i}^2` the signal and noise powers. These are idealizations; real noise is often colored and partially correlated across sensors, but they lead to estimators and bounds that perform well in practice, and the framework extends to a general noise covariance when needed. + +The Nonlinear Measurement Equations +========================================== + +Collecting the :math:`N-1` reference-based range differences into a vector :math:`\mathbf{m}` with entries :math:`m_i = c\,\tau_{i1} = r_i - r_1`, the noiseless model is + +.. math:: + + \mathbf{m} = \mathbf{h}(\mathbf{u}), \qquad + h_i(\mathbf{u}) = |\mathbf{u}-\mathbf{s}_i| - |\mathbf{u}-\mathbf{s}_1|, + +and the noisy measurement is :math:`\tilde{\mathbf{m}} = \mathbf{h}(\mathbf{u}) + \boldsymbol{\varepsilon}`, where :math:`\boldsymbol{\varepsilon}` is the range-difference error induced by time-delay estimation errors. The function :math:`\mathbf{h}` is nonlinear because of the Euclidean norms, and this nonlinearity is the source of every algorithmic complication that follows. Two broad strategies address it: algebraically *linearize* by introducing an auxiliary variable (described in the next section), or *iteratively* linearize about a current estimate (described further below). + +************************************************* +Time-Delay Estimation (the Measurement Front End) +************************************************* + +Before any geometry can be exploited we must extract the delays :math:`\tau_{ij}` from the raw waveforms. This is the *time-delay estimation* (TDE) problem, and its accuracy ultimately caps the accuracy of the entire system. + +Cross-Correlation +======================== + +The natural estimator exploits the fact that :math:`x_i` and :math:`x_j` are noisy, shifted copies of the same waveform. Their cross-correlation, + +.. math:: + + R_{x_i x_j}(\tau) = \mathbb{E}\!\left[ x_i(t)\, x_j(t+\tau) \right], + +is maximized when the shift :math:`\tau` aligns the two copies, i.e. at :math:`\tau = \tau_{ij}`. The estimator is therefore + +.. math:: + + \hat{\tau}_{ij} = \arg\max_{\tau} \, \hat{R}_{x_i x_j}(\tau). + +In practice the correlation is usually computed efficiently in the frequency domain via the FFT (just like how large convolutions typically use an FFT), using the cross-power spectral density :math:`G_{x_i x_j}(f) = \mathcal{F}\{R_{x_i x_j}(\tau)\}` and an inverse transform. + +Python Simulation +================== + +Enough math, let's see how this all looks with a simple Python example. First we'll set up the simulation with some high level parameters such as emitter and sensor position, and the simulated sample rate, which is essentially how much spectrum the receivers will "see". + +.. code-block:: python + + import numpy as np + import matplotlib.pyplot as plt + from matplotlib.lines import Line2D + from itertools import combinations + from scipy.signal import firwin, lfilter + + sample_rate = 50e6 + c = 3e8 # speed of light [m/s] + snr_db = 10 # SNR of the received signal at each receiver [dB] + tx_len_samples = 1000 # samples to transmit + rx_positions = np.array([ + [65, 229], # Rx0 + [676, 123], # Rx1 + [153, 543], # Rx2 + ]) + num_rx = rx_positions.shape[0] + tx_position = np.array([153, 355]) + pairs = list(combinations(range(num_rx), 2)) # For 3 receivers it's (Rx0,Rx1), (Rx0,Rx2), (Rx1,Rx2) -> 3 pairs + +TDOA is not very dependent on the specific signal the transmitter emits, although the bandwidth of the signal does matter, so to keep this simple we'll have it transmit random noise that is band-limited to a specified bandwidth. If we were to use something like QPSK of the same bandwidth instead, nothing would really change. + +.. code-block:: python + + bandwidth = 20e6 + taps = firwin(numtaps=129, cutoff=bandwidth / 2, fs=sample_rate) + tx_signal = lfilter(taps, 1.0, np.random.randn(tx_len_samples) + 1j * np.random.randn(tx_len_samples)) + +Next we will simulate the receivers receiving the signal at a delay based on their position. We will use a fractional delay filter like we learned about in the :ref:`sync-chapter` Chapter. The rest of the code should look relatively straightforward. We make sure to apply unique AWGN per receiver. + +.. code-block:: python + + # Simulate what each receiver records + true_distances = np.linalg.norm(rx_positions - tx_position, axis=1) + true_delays = true_distances / c + unknown_tx_time = 1.234e-5 # seconds. arbitrary, unknown to receivers and we won't use it in any TDOA calcs + + # Calc the actual TDOAs to act as ground truth + for k, (a, b) in enumerate(pairs): + true_rd = true_distances[b] - true_distances[a] + + # Figure out how many samples we have to simulate + total_delay_samples = (unknown_tx_time + true_delays.max()) * sample_rate + buffer_len = tx_len_samples + int(np.ceil(total_delay_samples)) + 10 + + # Taken from Synchronization chapter + def frac_delay_filter(delay): # delay is in samples, but it can (and will be) not an integer + N = 21 # number of taps, keep this odd + n = np.arange(-(N-1)//2, N//2+1) # -10,-9,...,0,...,9,10 + h = np.sinc(n - delay) # calc filter taps + h *= np.hamming(N) # window the filter to make sure it decays to 0 on both sides + h /= np.sum(h) # normalize to get unity gain, we don't want to change the amplitude/power + return h + + # Simulate the delayed signal being received by each sensor + rx_signals = np.zeros((num_rx, buffer_len), dtype=complex) + for i in range(num_rx): + tau = unknown_tx_time + true_delays[i] # absolute delay at this Rx, in seconds + tau_samples = tau * sample_rate + tau_integer_samps = int(np.round(tau_samples)) + tau_frac_samps = tau_samples - tau_integer_samps + rx = np.zeros(buffer_len, dtype=complex) + rx[tau_integer_samps:tau_integer_samps+tx_len_samples] = tx_signal + frac_delay_i = frac_delay_filter(tau_frac_samps) + rx = np.convolve(rx, frac_delay_i, "same") + + # Each receiver adds its own thermal noise, scaled to hit the SNR set at the top + signal_power = np.mean(np.abs(tx_signal)**2) + noise_power = signal_power / 10**(snr_db / 10) + noise = np.sqrt(noise_power / 2) * (np.random.randn(buffer_len) + 1j * np.random.randn(buffer_len)) + rx_signals[i] = rx + noise + +Everything so far was purely for simulation, the rest represents what you would actually do to calculate the TDOA, typically at a central location or one of the sensors, but it needs access to the samples received at all three sensors. It's not a lot of code, we simply loop through each pair of sensors, calculate the cross-correlation between their received samples, and pull out the peak. Later we will see how to do the subsample version of this for more granularity. + +.. code-block:: python + + # Estimate the TDOAs using a normal cross-correlation + range_diff = np.zeros(len(pairs)) # meters + for k, (a, b) in enumerate(pairs): + xcorr = np.correlate(rx_signals[b], rx_signals[a], mode='full') + peak_lag = np.argmax(np.abs(xcorr)) - (buffer_len - 1) # 'full' puts zero lag at index buffer_len-1 + range_diff[k] = (peak_lag / sample_rate) * c # meters + +Not much to it! This gives us the following results: + +.. image:: ../_images/tdoa_python_integer.svg + :align: center + :target: ../_images/tdoa_python_integer.svg + :alt: Python simulation output when doing integer correlation + +Note that this code doesn't fully "solve" the problem, even though it might seem like it does at first glance because the lines intersect exactly at the position of the emitter, but it's really your brain doing the final "solving" of the position, by looking at the intersection of the hyperbolas. Also, if there was more noise, the hyperbolas would not all intersect at one point. We will dive into automated solutions later in this chapter. + +The full Python code (including the plotting portion) can be found `here `_. + +Resolution and Sub-Sample Estimation +========================================== + +With sampling rate :math:`f_s`, the correlation is computed on a lag grid spaced :math:`1/f_s` apart, so the naive peak resolution is one sample, i.e. :math:`c/f_s` in range. This is usually far too coarse, especially if the sensors (and emitter) are close together, e.g., less than 100 meters. There are two options for sub-sample refinement, one way is to interpolate the signals as part of the cross-correlation, and another is to fit a model to the samples around the discrete peak. For the latter, parabolic interpolation through the peak and its two neighbors is the simplest, while sinc-based interpolation is more accurate because the true correlation of a band-limited signal is a sinc-like function. Good interpolation routinely yields delay estimates one to two orders of magnitude finer than the sample period. Below we show an example of doing the interpolated cross-correlation. + +.. code-block:: python + + U = 16 # correlation upsampling factor + half = (buffer_len + 1) // 2 # number of DC + positive-frequency bins + range_diff = np.zeros(len(pairs)) # meters + for k, (a, b) in enumerate(pairs): + # Cross-correlation in the frequency domain + X = np.conj(np.fft.fft(rx_signals[a])) * np.fft.fft(rx_signals[b]) + + # Insert zeros in the high-frequency MIDDLE: DC + positive freqs at the front, negative freqs at the back, so it stays a valid FFT layout. + X_padded = np.zeros(U * buffer_len, dtype=complex) + X_padded[:half] = X[:half] + X_padded[U * buffer_len - (buffer_len - half):] = X[half:] + + # Now IFFT to finish the crosscorrelation + xcorr = np.abs(np.fft.ifft(X_padded)) * U + + # Peak index -> signed lag; indices past the midpoint are negative lags + peak_idx = np.argmax(xcorr) + if peak_idx > U * buffer_len // 2: + peak_idx -= U * buffer_len + peak_lag = peak_idx / U # sub-sample lag, +ve => Rx_b farther + range_diff[k] = (peak_lag / sample_rate) * c # meters + +When doing the same Python simulation as before, but with subsampling, we get the following results. You would have to zoom in on the left-hand plot to see the accuracy difference. + +.. image:: ../_images/tdoa_python_subsample.svg + :align: center + :target: ../_images/tdoa_python_subsample.svg + :alt: Python simulation output when doing subsample correlation + +Looking at the right-hand plot, we can see how the original integer-only method was off by a decent margin. + +The Generalized Cross-Correlation Framework +================================================== + +Plain cross-correlation is fragile: if the transmitted signal is narrow in bandwidth or the channel contains multipath, the correlation peak is broad and easily shifted by noise. Knapp and Carter's *Generalized Cross-Correlation* (GCC) addresses this by inserting a frequency weighting :math:`\Psi(f)` before transforming back to the lag domain: + +.. math:: + + R^{\mathrm{GCC}}_{x_i x_j}(\tau) = \int_{-\infty}^{\infty} \Psi(f)\, G_{x_i x_j}(f)\, e^{j 2\pi f \tau}\, df . + +The weighting reshapes the spectrum to sharpen and stabilize the peak. Different choices of :math:`\Psi(f)` correspond to different classical estimators, and selecting it well is the heart of robust TDE. Common weightings include: + +* **Cross-correlation** (:math:`\Psi = 1`): the maximum-likelihood choice only in the high-SNR, broadband-flat limit; otherwise suboptimal. +* **Roth** (:math:`\Psi = 1/G_{x_i x_i}(f)`): suppresses frequencies where one sensor is noisy. +* **SCOT** (Smoothed Coherence Transform, :math:`\Psi = 1/\sqrt{G_{x_i x_i}G_{x_j x_j}}`): symmetric whitening of both channels. +* **PHAT** (Phase Transform, :math:`\Psi = 1/|G_{x_i x_j}(f)|`): the most widely used choice in acoustics. + +Diving deeper into the **GCC-PHAT** estimator- by dividing out the magnitude of the cross-spectrum it retains *only the phase*: + +.. math:: + + R^{\mathrm{PHAT}}_{x_i x_j}(\tau) = \int \frac{G_{x_i x_j}(f)}{\bigl|G_{x_i x_j}(f)\bigr|} e^{j2\pi f \tau} df . + +Because the delay between two copies of a signal is encoded entirely in the *linear phase* term :math:`e^{-j2\pi f \tau_{ij}}`, while the magnitude carries the (often unhelpful) spectral shape and reverberant coloring, whitening to unit magnitude weights every frequency equally and produces a sharp, near-impulsive peak at the true delay. This makes PHAT strikingly robust to multipath. Its weakness is that it also whitens noise-dominated frequencies, so at low SNR the equal weighting amplifies noise; SNR-aware variants reintroduce a coherence-based weighting to compensate. In real systems, because most signals are not always on, you typically have to determine the time-frequency bounding box of the target signal before performing TDOA, so unless there is a lot of interference you can estimate the SNR pretty easily. + +Practical Considerations +================================ + +Several effects govern the accuracy of the TDOA results: + +* The **integration window** :math:`T` trades estimator variance (longer is better, since variance falls roughly as :math:`1/T`) against the stationarity assumption and, for moving emitters, against blurring of the delay over the window. Many times this value is determined by the signal itself, e.g. you might perform TDOA on a per-packet basis. +* **Coherence bandwidth** limits which frequencies actually carry usable phase. +* **Signal bandwidth** is decisive: as the Cramér-Rao analysis below shows, delay variance falls as the *square* of the bandwidth, so wideband signals localize far better than narrowband ones, unlike DOA where we didn't care about bandwidth (in fact, many of the DOA concepts used a narrowband assumption). That being said, you don't need to capture the entire signal (from a frequency domain perspective) to perform TDOA, if you are limited by your SDR's maximum sample rate, and you can only receive a portion of the signal bandwidth, you can still do TDOA! + +From a compute perspective, the TDOA computation is dominated by FFTs and is :math:`O(M\log M)` per sensor pair for records of :math:`M` samples, which is what makes large sensor networks tractable. + +GCC-PHAT Python Example +============================================ + +The nice thing about PHAT is that we can fold it into the simulation we already built with almost no new code. Recall that the sub-sample estimator above already worked in the frequency domain: it formed the cross-spectrum :math:`X_a^*(f)\,X_b(f)`, zero-padded it to interpolate, and inverse-FFT to recover the lag-domain correlation. PHAT is just one extra line: before transforming back to the lag domain, we divide the cross-spectrum by its own magnitude, so that every frequency bin contributes with unit weight and only the phase, which is where the delay lives, survives. + +.. code-block:: python + + U = 16 # correlation upsampling factor + half = (buffer_len + 1) // 2 # number of DC + positive-frequency bins + range_diff = np.zeros(len(pairs)) # meters + for k, (a, b) in enumerate(pairs): + # Cross-spectrum, same as the sub-sample example + X = np.conj(np.fft.fft(rx_signals[a])) * np.fft.fft(rx_signals[b]) + + # PHAT weighting: divide out the magnitude so only the phase remains + X = X / (np.abs(X) + 1e-12) # small epsilon avoids divide-by-zero + + # Zero-pad in the high-frequency middle to interpolate, then IFFT + X_padded = np.zeros(U * buffer_len, dtype=complex) + X_padded[:half] = X[:half] + X_padded[U * buffer_len - (buffer_len - half):] = X[half:] + xcorr = np.abs(np.fft.ifft(X_padded)) * U + + # Peak index -> signed lag; indices past the midpoint are negative lags + peak_idx = np.argmax(xcorr) + if peak_idx > U * buffer_len // 2: + peak_idx -= U * buffer_len + peak_lag = peak_idx / U # sub-sample lag, +ve => Rx_b farther + range_diff[k] = (peak_lag / sample_rate) * c # meters + +The only change from the sub-sample code is the single ``X = X / (np.abs(X) + 1e-12)`` line; the small epsilon in the denominator keeps frequency bins that hold almost no energy from blowing up when we divide. With our wideband, high-SNR simulated signal the result barely differs from plain cross-correlation, because PHAT and cross-correlation coincide in exactly that broadband, high-SNR limit. The payoff shows up in harder conditions: when the spectrum is colored or multipath smears the peak, whitening to unit magnitude collapses the correlation back to a sharp, near-impulsive spike at the true delay, which is why PHAT is the default in acoustics. + +************************************* +Closed-Form Localization Algorithms +************************************* + +The measurement equations above are nonlinear and, taken directly, require iterative solution with a good starting point. *Closed-form* (non-iterative) estimators sidestep this by an algebraic trick: introduce an auxiliary variable that absorbs the nonlinearity and renders the system linear. They are fast, need no initial guess, and cannot get stuck in local minima, making them invaluable both on their own and as initializers for the iterative methods described below. + +The Linearization Strategy +================================== + +The trick is to square the range equations and subtract pairs, which cancels the nonlinear :math:`x^2+y^2` term and introduces :math:`r_1`, the range to the reference sensor, as a single auxiliary unknown. Starting with the squared range from the emitter :math:`\mathbf{u}=(x,y)` to sensor :math:`i` at :math:`\mathbf{s}_i=(x_i,y_i)`: + +.. math:: + + r_i^2 = (x-x_i)^2 + (y-y_i)^2 = K_i - 2x_i x - 2y_i y + (x^2+y^2), + \qquad K_i \equiv x_i^2 + y_i^2 . + +The troublesome term is :math:`x^2+y^2`, common to every sensor. Take sensor 1 as reference and subtract its equation from sensor :math:`i`'s: + +.. math:: + + r_i^2 - r_1^2 = (K_i - K_1) - 2(x_i-x_1)x - 2(y_i-y_1)y . + +Now use the measured range difference :math:`r_{i1}\equiv r_i - r_1 = c\,\tau_{i1}`. Since :math:`r_i = r_{i1}+r_1`, we have :math:`r_i^2 = r_{i1}^2 + 2r_{i1}r_1 + r_1^2`, so :math:`r_i^2 - r_1^2 = r_{i1}^2 + 2 r_{i1} r_1`. Substituting and rearranging, + +.. math:: + + \boxed{2(x_i-x_1)\,x + 2(y_i-y_1)\,y + 2 r_{i1} r_1 = K_i - K_1 - r_{i1}^2} + +This equation is **linear** in the unknowns :math:`(x, y, r_1)`, where the range to the reference :math:`r_1` is treated as an auxiliary variable. Stacking it for :math:`i=2,\dots,N` gives a linear system :math:`\mathbf{A}\boldsymbol{\theta} = \mathbf{b}` with :math:`\boldsymbol{\theta}=[x,y,r_1]^\top`, solvable by ordinary or weighted least squares. The nonlinearity has been quarantined into the single extra unknown :math:`r_1`. + +Spherical Interpolation and Spherical Intersection +========================================================= + +The earliest closed-form estimators, Spherical Interpolation (SI) and Spherical Intersection (SX) of Schau and Robinson, exploit exactly this structure. They first solve the linear system for :math:`(x,y)` as a function of :math:`r_1`, then impose the constraint that ties them together, namely :math:`r_1^2 = (x-x_1)^2+(y-y_1)^2`, to pin down :math:`r_1`. SI obtains :math:`r_1` by a least-squares projection; SX substitutes the linear solution into the quadratic constraint and solves the resulting scalar quadratic. They are simple and fast but treat the auxiliary variable somewhat crudely, leaving accuracy on the table at higher noise. + +Fang's Method +==================== + +Fang's algorithm provides an exact algebraic solution for the *minimum* configuration (3 sensors in 2D, 4 in 3D), giving a determined system rather than an overdetermined one. It is elegant and computationally trivial but does not use redundant sensors, so it cannot average down measurement noise and is sensitive to geometry. It is best viewed as the exact-determined special case that the least-squares methods generalize. + +To see how it works, look again at the boxed linear equation above. With three sensors there are exactly two such equations (:math:`i=2,3`) but three unknowns :math:`(x,y,r_1)`. That looks underdetermined, but :math:`r_1` is not free: it is glued to the position by :math:`r_1^2=(x-x_1)^2+(y-y_1)^2`. Fang's trick is to *defer* that constraint, treat :math:`r_1` as a known constant, and solve the two linear equations for :math:`x` and :math:`y`. Because :math:`r_1` enters linearly, inverting the :math:`2\times2` matrix (well-conditioned as long as the sensors are not collinear) gives both coordinates as straight-line functions of the still-unknown range, + +.. math:: + + x = g_x + h_x\,r_1, \qquad y = g_y + h_y\,r_1 , + +with constants :math:`g_x,h_x,g_y,h_y` from the matrix inverse. Now cash in the deferred constraint: substituting these into :math:`r_1^2=(x-x_1)^2+(y-y_1)^2` collapses everything to a single scalar quadratic :math:`a\,r_1^2 + b\,r_1 + c = 0`. Solve it, keep the physical root (a range must be positive; the other root typically lands on the wrong hyperbola branch), and back-substitute to read off :math:`(x,y)`. That is the whole method: one :math:`2\times2` solve and one quadratic, no iteration and no initial guess. The worked example in the next subsection is precisely this procedure carried out with numbers. + +Chan's Method (Two-Step Weighted Least Squares) +====================================================== + +The estimator that became the practical standard is Chan and Ho's two-step weighted least squares (WLS). It is built on the linear system above but treats the statistics correctly and refines the auxiliary variable, achieving accuracy close to the Cramér-Rao bound at small-to-moderate noise. + +**First step.** Treat :math:`\boldsymbol{\theta}=[x,y,r_1]^\top` as if its three components were independent and solve the linear system by weighted least squares, + +.. math:: + + \hat{\boldsymbol{\theta}} = (\mathbf{A}^\top \mathbf{W}\mathbf{A})^{-1}\mathbf{A}^\top \mathbf{W}\,\mathbf{b}, + +with the weight :math:`\mathbf{W}` chosen as the inverse covariance of the equation errors. Because that covariance itself depends on the unknown ranges, in practice one first solves with :math:`\mathbf{W}=\mathbf{I}` (or the raw TDOA noise covariance), then recomputes :math:`\mathbf{W}` from the resulting range estimates and re-solves, a one- or two-pass refinement. + +**Second step.** The first step ignored the known relationship :math:`r_1^2 = (x-x_1)^2+(y-y_1)^2` that couples the auxiliary variable to the position. The second step restores it: form a new small least-squares problem in the squared quantities :math:`[(x-x_1)^2,(y-y_1)^2,r_1^2]`, using the first-step covariance to weight it, and solve for a corrected position. This second WLS removes much of the bias of the naive linear solution and is what brings Chan's estimator close to optimal. + +The method returns a position directly, with computational cost dominated by inverting small :math:`3\times3` matrices, negligible compared with the FFTs of the front end. Its limitations appear at high noise or unfavorable geometry, where the squared-range manipulation amplifies errors and the second step can pick the wrong root; there, the iterative refinement described below, seeded by Chan's output, is the standard remedy. + +Example, Continued: Solving the Three-Sensor Fix in Closed Form +================================================================ + +Let's pick up right where the Python simulation left off. At that point we had the ``range_diff`` array holding one measured range difference per sensor pair, and earlier we let our brain do the final step by eyeballing where the hyperbolas crossed. Now we'll replace that eyeballing with the closed-form algebra developed above, recovering the emitter position directly from ``range_diff`` and ``rx_positions``. Because we have exactly three sensors in 2D, this is Fang's minimum-configuration case: two boxed linear equations and one quadratic, no iteration and no initial guess. + +We take ``Rx0`` as the reference sensor. The pairs were built as ``(0,1)``, ``(0,2)``, ``(1,2)``, and recall that ``range_diff[k]`` for pair ``(a,b)`` is :math:`r_b - r_a`, so the two pairs that include the reference, ``(0,1)`` and ``(0,2)``, hand us exactly the reference-based range differences :math:`r_{i0}=r_i-r_0` that the boxed equation needs. + +.. code-block:: python + + # Solve for the emitter position in closed form (Fang's method, 3 sensors in 2D) + ref = 0 # use Rx0 as the reference sensor + s = rx_positions.astype(float) + K = np.sum(s**2, axis=1) # K_i = x_i^2 + y_i^2 for each sensor + + # Reference-based range differences r_i0 = r_i - r_ref for the two non-reference sensors + others = [i for i in range(num_rx) if i != ref] + r_i0 = np.array([range_diff[pairs.index((ref, i))] for i in others]) # pair (ref,i) holds r_i - r_ref + + # Build the 2x2 linear system that gives (x, y) as a function of the unknown range r_ref + M = 2 * (s[others] - s[ref]) # rows: [2(x_i - x_ref), 2(y_i - y_ref)] + d = K[others] - K[ref] - r_i0**2 # right-hand side constants + Minv = np.linalg.inv(M) # well-conditioned as long as the sensors aren't collinear + g = Minv @ d # part of (x, y) that doesn't depend on r_ref + h = -2 * (Minv @ r_i0) # how (x, y) slide with r_ref: [x, y] = g + h * r_ref + + # Cash in the deferred constraint r_ref^2 = (x - x_ref)^2 + (y - y_ref)^2 -> scalar quadratic in r_ref + p = g - s[ref] # constant part of (x - x_ref, y - y_ref) + a_q = h[0]**2 + h[1]**2 - 1 + b_q = 2 * (p[0]*h[0] + p[1]*h[1]) + c_q = p[0]**2 + p[1]**2 + roots = np.roots([a_q, b_q, c_q]) + + # Keep the physical root (a range must be positive and real), then back-substitute + r_ref = roots[(roots.real > 0) & (np.abs(roots.imag) < 1e-6)].real.max() + emitter_est = g + h * r_ref + + print("Estimated emitter position:", emitter_est) # ~[153, 355] + print("True emitter position: ", tx_position) + +The structure mirrors the math exactly: ``M`` and ``d`` are the two boxed linear equations, ``g`` and ``h`` express :math:`x` and :math:`y` as straight-line functions of the still-unknown reference range :math:`r_1` (called ``r_ref`` here), and substituting those into :math:`r_1^2=(x-x_1)^2+(y-y_1)^2` collapses everything to the scalar quadratic that ``np.roots`` solves. We discard the non-physical (negative or complex) root, keep the positive real one, and back-substitute to read off the position. With our high-SNR, wideband simulation the estimate lands right on top of the true emitter at :math:`(153, 355)`, with no human in the loop reading off a hyperbola intersection. + +With noisier measurements the two linear equations would no longer be perfectly consistent, the quadratic root would be perturbed, and, because three sensors give us no redundancy to average over, the error would pass straight through. That is exactly where the redundant pairs and the weighting and second step of Chan's method earn their keep, governing how gracefully the estimate degrades. + +***************************************** +Iterative and Statistical Estimation +***************************************** + +Closed-form methods are fast but make algebraic approximations that cost accuracy at high noise or poor geometry. When the best possible estimate is required, we solve the nonlinear estimation problem directly, typically initialized by a closed-form result. + +Nonlinear Least Squares +============================== + +Define the residual between measured and predicted range differences and minimize its weighted squared norm: + +.. math:: + + \hat{\mathbf{u}} = \arg\min_{\mathbf{u}} \bigl[\tilde{\mathbf{m}} - \mathbf{h}(\mathbf{u})\bigr]^\top \mathbf{C}^{-1} \bigl[\tilde{\mathbf{m}} - \mathbf{h}(\mathbf{u})\bigr], + +where :math:`\mathbf{C}` is the covariance of the range-difference errors. This cost has no closed-form minimizer because :math:`\mathbf{h}` is nonlinear, so we descend it iteratively. + +Taylor-Series (Gauss-Newton) Method +========================================== + +Foy's classical approach linearizes :math:`\mathbf{h}` about the current estimate :math:`\mathbf{u}^{(k)}` using its Jacobian :math:`\mathbf{J}`, whose row :math:`i` is the gradient of :math:`h_i`: + +.. math:: + + \frac{\partial h_i}{\partial \mathbf{u}} = \frac{\mathbf{u}-\mathbf{s}_i}{|\mathbf{u}-\mathbf{s}_i|} - \frac{\mathbf{u}-\mathbf{s}_1}{|\mathbf{u}-\mathbf{s}_1|} + = \hat{\mathbf{e}}_i - \hat{\mathbf{e}}_1, + +a difference of *unit vectors* pointing from the candidate emitter toward sensor :math:`i` and the reference. The Gauss-Newton update is + +.. math:: + + \mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} + (\mathbf{J}^\top \mathbf{C}^{-1}\mathbf{J})^{-1}\mathbf{J}^\top \mathbf{C}^{-1}\bigl[\tilde{\mathbf{m}}-\mathbf{h}(\mathbf{u}^{(k)})\bigr], + +iterated to convergence. Each step solves a small linear system. The method converges quickly *when started near the solution*, which is exactly why Chan's closed-form estimate is the preferred seed: it places the iteration in the basin of the global minimum and avoids the spurious local minima that plague hyperbolic cost surfaces, especially in poor geometry. + +Maximum-Likelihood Estimation +===================================== + +Under Gaussian noise, the negative log-likelihood is, up to constants, exactly the weighted squared residual above. So **the maximum-likelihood estimator coincides with weighted nonlinear least squares**, the Gauss-Newton iteration is not a heuristic, it is the statistically optimal estimator under the assumed model. This is also the estimator whose covariance the Cramér-Rao bound below predicts. + +Let's put that to work by continuing the Python example one more time. We already have a position from the closed-form solver, ``emitter_est``, and the theory tells us two things: the maximum-likelihood estimate is just the Gauss-Newton iteration above, and the closed-form fix is the ideal seed for it because it drops us right inside the basin of the true minimum. So we'll start at ``emitter_est`` and take a few Gauss-Newton steps, each one re-linearizing the range-difference model at the current guess and solving a tiny least-squares problem for the correction. Unlike Fang's solver, which used only the two pairs touching the reference sensor, this one uses *all three* pairs in ``range_diff``, so the extra pair acts as redundancy that the iteration averages over. + +.. code-block:: python + + # Refine the closed-form fix with Gauss-Newton (= maximum likelihood under Gaussian noise) + u = emitter_est.copy() # seed the iteration with the closed-form estimate + for _ in range(10): + h = np.zeros(len(pairs)) # predicted range differences at the current guess + J = np.zeros((len(pairs), 2)) # Jacobian, one row per pair + for k, (a, b) in enumerate(pairs): + e_a = (u - s[a]) / np.linalg.norm(u - s[a]) # unit vector from Rx_a toward the guess + e_b = (u - s[b]) / np.linalg.norm(u - s[b]) # unit vector from Rx_b toward the guess + h[k] = np.linalg.norm(u - s[b]) - np.linalg.norm(u - s[a]) # predicted r_b - r_a + J[k] = e_b - e_a # row of the Jacobian is a difference of unit bearing vectors + + residual = range_diff - h # measured minus predicted range differences + delta, *_ = np.linalg.lstsq(J, residual, rcond=None) # Gauss-Newton step (J^T J)^-1 J^T residual + u = u + delta + if np.linalg.norm(delta) < 1e-9: # stop once the update stops moving the estimate + break + + emitter_ml = u + print("ML (Gauss-Newton) estimate:", emitter_ml) # ~[153, 355] + print("True emitter position: ", tx_position) + +A couple of details worth pointing out. Because we assumed the range-difference errors are independent with equal variance, the weight :math:`\mathbf{C}^{-1}=\sigma^{-2}\mathbf{I}` is a scalar that cancels out of the update, which is why a plain ``np.linalg.lstsq`` (no weight matrix) computes the step exactly; if the pairs had unequal quality we would fold their inverse variances in here. The Jacobian rows are literally the ``e_b - e_a`` differences of unit bearing vectors from the math above, so you can watch the geometry enter the estimator directly. Starting from the already-good closed-form seed, the iteration converges in just a handful of steps and lands on the true emitter at :math:`(153, 355)`. In our high-SNR simulation it barely moves off the closed-form answer, but with noisier measurements this is where the extra pair and the iterative refinement pay off, and it is this same :math:`\mathbf{J}^\top\mathbf{C}^{-1}\mathbf{J}` that reappears in the Cramér-Rao bound below as the estimator's covariance. + +Robust, Recursive, and Bayesian Extensions +================================================== + +Real measurements contain outliers, a multipath-corrupted TDOA can be wildly wrong while the rest are fine. Plain least squares, which squares residuals, is badly distorted by such outliers. *Robust* estimators replace the squared loss with one that grows more slowly (e.g. Huber's), or explicitly detect and discard inconsistent TDOAs via residual tests or RANSAC-style consensus. + +When the emitter *moves*, we want to fuse measurements over time rather than localize each instant independently. State-space filtering does this by modeling the emitter's position (and velocity) as an evolving state. The **Kalman filter** is optimal for linear-Gaussian dynamics, but the TDOA measurement is nonlinear, so practitioners use the **Extended Kalman Filter** (which linearizes the measurement with the same Jacobian as above), the **Unscented Kalman Filter** (which propagates a deterministic set of sigma points through the nonlinearity, avoiding explicit Jacobians and handling stronger nonlinearity better), or, for multimodal or heavily non-Gaussian problems, the **particle filter** (which represents the posterior by a weighted sample cloud). These trackers also naturally enforce motion continuity, which suppresses the per-snapshot ambiguities of static localization. + +**************************** +Brute-Force Heatmap Approach +**************************** + +Every method so far has been algebraic or iterative: we manipulated equations or descended a gradient. But there is a refreshingly simple alternative that needs neither. Lay a grid over the search area, and at every candidate position ask a single question: *if the emitter were here, what range differences would the sensors see, and how far off are those from what we actually measured?* Squaring and summing those mismatches gives a cost at each grid point, and the emitter is wherever that cost is smallest. The result is a heatmap of the same cost surface the Gauss-Newton iteration was quietly walking down, except now we can see all of it at once. + +Back to the Python example, we can do this with the variables we already have, ``range_diff``, ``rx_positions``, and ``pairs``: + +.. code-block:: python + + # Evaluate the TDOA cost on a grid of candidate emitter positions + gx = np.linspace(0, 700, 400) + gy = np.linspace(0, 700, 400) + GX, GY = np.meshgrid(gx, gy) + + cost = np.zeros_like(GX) + for k, (a, b) in enumerate(pairs): + r_a = np.hypot(GX - rx_positions[a, 0], GY - rx_positions[a, 1]) # range to Rx_a + r_b = np.hypot(GX - rx_positions[b, 0], GY - rx_positions[b, 1]) # range to Rx_b + cost += ((r_b - r_a) - range_diff[k])**2 # squared mismatch for this pair, summed over pairs + + # The best estimate is simply the grid cell with the lowest cost + iy, ix = np.unravel_index(np.argmin(cost), cost.shape) + emitter_grid = np.array([gx[ix], gy[iy]]) + print("Grid estimate:", emitter_grid) # ~[153, 355] + + # Invert the cost into a likelihood-style surface so higher = more likely emitter location + likelihood = -np.log10(cost + 1e-9) + +Plotting ``likelihood`` as an image reveals the geometry directly: the bright ridges trace out the hyperbolas from earlier, and they all funnel into one bright peak at the true emitter. We take the negative log of the cost so that the most likely location is the maximum rather than a minimum, which is easier to read off visually. The trade-offs are exactly what you'd expect. The method is dead simple, needs no initial guess, and cannot diverge or land on the wrong root, so it is a great sanity check and a robust way to *seed* the iterative refiner. It also handles multimodal cost surfaces gracefully, since it sees every minimum, not just the nearest one. The price is resolution and speed: accuracy is limited by the grid spacing, and cost grows with the number of grid cells, so for a fine answer over a large area you would localize coarsely first and then refine, either by zooming the grid or by handing the result to Gauss-Newton. Below shows the heatmap approach applied to our Python example. + +.. image:: ../_images/tdoa_python_heatmap.svg + :align: center + :target: ../_images/tdoa_python_heatmap.svg + :alt: Adding heatmap to the tdoa plot shown earlier + +One nice part about the heatmap approach is if there is a lot of error, or sensors with low SNR without realizing it, there may be multiple hot spots on the heatmap, which your brain can notice. The heatmap can even be overlaid on top of a satellite view of the area! + +This brute-force approach is not very computationally efficient, and it's not really an option at all for 3D TDOA. One alternative to calculating every grid point but still brute-forcing it is to draw all of the hyperbolas in 2D but with "width" applied to each one, e.g. by applying a lobe shaped function along the hyperbola so it tapers off. + +*********************************************** +Performance Analysis and Fundamental Bounds +*********************************************** + +Having estimators in hand, we ask: how accurate *can* a TDOA system be, and what governs that accuracy? Two ideas answer this: the Cramér-Rao bound, which sets a noise floor from the signals, and geometric dilution of precision, which describes how sensor-emitter geometry amplifies that floor. + +Error Propagation +======================== + +System accuracy is a two-stage cascade. First, finite SNR and bandwidth limit how precisely each delay can be measured (TDE error). Second, the geometry maps those range-difference errors into a position error. Writing :math:`\delta\mathbf{u}` for the position error and :math:`\boldsymbol{\varepsilon}` for the range-difference errors, the linearized relation near the solution is :math:`\boldsymbol{\varepsilon}\approx \mathbf{J}\,\delta\mathbf{u}`, so the position-error covariance is + +.. math:: + + \mathrm{Cov}(\hat{\mathbf{u}}) \approx (\mathbf{J}^\top \mathbf{C}^{-1}\mathbf{J})^{-1}. + +This single expression contains both stages: :math:`\mathbf{C}` is the measurement quality (from TDE) and :math:`\mathbf{J}` is the geometry. + +The Time-Delay Estimation Bound +======================================= + +We can bound how well any estimator can measure a single delay. For a signal of RMS bandwidth :math:`\beta` observed over time :math:`T`, the variance of any unbiased delay estimate obeys + +.. math:: + + \mathrm{var}(\hat\tau_{ij}) \gtrsim \frac{1}{8\pi^2 \beta^2 T \gamma}, + +where :math:`\beta` is the *RMS (Gabor) bandwidth* of the signal and :math:`\gamma` is an effective SNR factor combining the two sensors' SNRs. Three design lessons fall straight out: variance improves with **integration time** :math:`T`, with **effective SNR** :math:`\gamma`, and, most strikingly, with the **square of bandwidth** :math:`\beta^2`. Doubling the bandwidth quarters the delay variance. This is why wideband and spread-spectrum waveforms are so prized for ranging, and why narrowband emitters are intrinsically hard to localize by TDOA alone. + +The Localization Cramér-Rao Lower Bound +============================================== + +Combining measurement quality and geometry, the Fisher information matrix for the emitter position is + +.. math:: + + \mathbf{F} = \mathbf{J}^\top \mathbf{C}^{-1} \mathbf{J}, + +and the Cramér-Rao Lower Bound states that *any* unbiased estimator has covariance no smaller than its inverse: + +.. math:: + + \mathrm{Cov}(\hat{\mathbf{u}}) \succeq \mathbf{F}^{-1} = (\mathbf{J}^\top \mathbf{C}^{-1}\mathbf{J})^{-1}. + +The bound is the benchmark against which estimators are judged: a method that attains it is *efficient*. The maximum-likelihood estimator above attains it asymptotically (large :math:`T`, high SNR), and Chan's closed-form method attains it at small noise, which is exactly why both are used. The CRLB also cleanly separates the two influences on accuracy: :math:`\mathbf{C}` (signal-and-noise quality, improvable by more bandwidth, power, or integration) and :math:`\mathbf{J}` (geometry, improvable by sensor placement), studied next. The plot below shows a few example bandwidths and the lower bound over SNR, to give you a feel for how much error you should expect, or at least the floor. The y-axis is the 1-:math:`\mathrm{\sigma}` value (one standard deviation). + +.. image:: ../_images/tdoa_cramer_rao.svg + :align: center + :target: ../_images/tdoa_cramer_rao.svg + :alt: Plot of cramer rao lower bound + + +Geometric Dilution of Precision +======================================= + +Suppose your sensors can measure range differences to about 1 m of accuracy, a respectable number for a well-synchronized radio system. You might expect to then pin down the emitter to roughly 1 m as well. But where is the emitter? Picture it sitting comfortably inside a triangle of three sensors: the hyperbolas from each sensor pair slice across one another at steep, nearly right angles, and where they cross is pinned down tightly, so your 1 m of ranging error turns into maybe 1.5 m of position error. Now slide that same emitter far off to one side, well outside the cluster. The hyperbolas now graze each other at a shallow angle, like two gently curving lines that nearly overlap, and the crossing point smears out along the direction they share. The very same 1 m of ranging error can now balloon into tens of meters of position error. Nothing about your hardware changed, only the geometry did. + +That blow-up factor has a name: **Geometric Dilution of Precision** (GDOP). It captures how much the sensor-emitter layout magnifies measurement error into position error. If the range-difference errors are independent and each has the same standard deviation :math:`\sigma`, so the covariance is :math:`\mathbf{C}=\sigma^2\mathbf{I}` (a diagonal matrix with :math:`\sigma^2` on the diagonal), then + +.. math:: + + \mathrm{GDOP} = \sqrt{\mathrm{tr}\bigl[(\mathbf{J}^\top\mathbf{J})^{-1}\bigr]}, \qquad + \sigma_{\text{position}} = \mathrm{GDOP}\cdot \sigma . + +Your position error is just your ranging error multiplied by GDOP, so GDOP is a unitless number, always :math:`\ge 1`, telling you the factor by which ranging error gets magnified at a given emitter location. + +Where does the magnification come from? It is baked into the Jacobian :math:`\mathbf{J}`, whose rows are differences of unit bearing vectors :math:`\hat{\mathbf{e}}_i - \hat{\mathbf{e}}_1` (the direction to one sensor minus the direction to another). When those directions point all over the place, :math:`\mathbf{J}^\top\mathbf{J}` is *well-conditioned* (far from singular, so its inverse stays small) and GDOP is small. When they nearly line up, :math:`\mathbf{J}^\top\mathbf{J}` becomes nearly singular and GDOP blows up. So an emitter surrounded by the sensors, with bearing vectors well-spread and hyperbolas crossing at large angles, gets a small GDOP (good), while an emitter far outside the cluster, or sensors nearly collinear (almost in a straight line), leaves the bearing vectors nearly parallel and the hyperbolas grazing at shallow angles, giving a huge GDOP (bad). + +This is the same effect we saw when hyperbolas degenerate near the ends of the baseline. The takeaway: a TDOA system can be limited far more by *where its sensors sit* than by *how well it measures time*, all the nanosecond synchronization and wide bandwidth in the world won't save a fix in a high-GDOP region of the map. + +The figure below shows GDOP heat maps over a plane for (left) three sensors at the vertices of an equilateral triangle and (right) three nearly collinear sensors, showing a broad low-GDOP region inside the triangle versus a narrow usable corridor for the collinear array, with GDOP rising sharply outside the convex hull in both cases. + +.. image:: ../_images/tdoa_gdop.svg + :align: center + :target: ../_images/tdoa_gdop.svg + :alt: GDOP heat maps over a plane for (left) three sensors at the vertices of an equilateral triangle and (right) three nearly collinear sensors, showing a broad low-GDOP region inside the triangle versus a narrow usable corridor for the collinear array, with GDOP rising sharply outside the convex hull in both cases. + +Sensor-Placement Optimization +===================================== + +Because geometry is often a *design* variable, we can place sensors to minimize error. Common objectives minimize a scalar derived from :math:`\mathbf{F}^{-1}`, such as its trace (equivalent to GDOP), its determinant (the confidence-ellipse volume), or its largest eigenvalue (worst-case error). The qualitative results are intuitive: spread the sensors widely so long baselines sharpen angular resolution, surround the region of interest so emitters fall inside the convex hull, avoid collinear or coplanar layouts that create ill-conditioned directions, and add sensors where redundancy both lowers variance and guards against outliers. For a moving target or large coverage area, placement is optimized over the whole region, minimizing average or worst-case GDOP, usually by numerical search. + +***************************************** +Practical Challenges in Real Systems +***************************************** + +The model used in this chapter so far omits a few effects that usually dominate the error budget in actual TDOA deployments. Three deserve detailed treatment: + +Receiver Synchronization +================================ + +TDOA's defining advantage, that it needs no synchronized transmitter, comes paired with its defining burden: the *receivers* must share a common time reference, and any error in that reference enters the measurement directly. If sensor :math:`i`'s clock is offset from truth by :math:`\delta t_i`, the measured TDOA is corrupted by :math:`\delta t_i - \delta t_j`, an error multiplied by :math:`c` in range. The scale is unforgiving for radio systems: + +.. math:: + + c \times 1\ \text{ns} = (3\times10^8\,\text{m/s})(10^{-9}\,\text{s}) = 0.30\ \text{m}. + +So a 1 ns synchronization error already costs :math:`\sim`\0.3 m, and 100 ns costs 30 m. Achieving and holding nanosecond-level synchronization across distributed sensors is therefore central to system design. Common mechanisms include GPS-disciplined oscillators (each sensor recovers a :math:`\sim`\10-100 ns timing reference from satellites), the Precision Time Protocol (IEEE 1588, distributing time over a network to sub-microsecond or, with hardware timestamping, sub-100 ns accuracy), and for the most demanding installations White Rabbit, which reaches sub-nanosecond synchronization over fiber. Two further subtleties matter: clocks not only have a static *offset* but *drift* over time, requiring continual discipline; and in acoustic systems, where :math:`c` is a million times smaller, the same absolute timing error is a million times less harmful, which is why microphone-array TDOA is comparatively forgiving while radio TDOA lives or dies by its clocks. + +Off-the-shelf SDRs that can be easily synchronized include any of the Ettus Research USRPs that can take a GPS disciplined oscillator (GPSDO), such as a `B200 `_ with a `TCXO `_ (a type of GPSDO). If the sensors are close enough, they can also be synchronized by sharing a PPS signal over a cable, e.g., generated using an `OctoClock `_, which also generates a 10 MHz signal used to frequency synchronize them. Nearly all of the USRPs have a PPS and 10 MHz input, and most have room for a GPSDO, or come with one. + +Multipath and Non-Line-of-Sight Propagation +=================================================== + +Everything so far has assumed a single line-of-sight path between the emitter and receivers. Real environments add reflections (multipath) and can block the direct path entirely. Multipath superimposes delayed copies of the signal, which distort or split the correlation peak and bias the delay estimate; this is exactly the failure GCC-PHAT was designed to resist, since whitening sharpens the direct-path peak relative to the smeared reflections. When the direct path is entirely obstructed, the *earliest* arriving energy travels an excess distance, so the measured TDOA is biased *long* in a way no amount of averaging removes, because the error is systematic rather than random. Mitigation strategies include identifying these non-line-of-sight links statistically (non-line-of-sight measurements often show larger variance or violate geometric consistency among redundant sensors), down-weighting or discarding them (requires having way more sensors than three), and exploiting redundancy so that a few corrupted links among many can be detected and rejected by the robust estimators described earlier. In dense indoor multipath, which is an extremely difficult environment for TDOA, model-based delay estimation and machine-learning approaches increasingly outperform classical correlation. + +Sensor-Position Uncertainty and Calibration +=================================================== + +The geometry assumed exact knowledge of the sensor coordinates :math:`\mathbf{s}_i`. Errors in those coordinates propagate into the position estimate just as measurement errors do, and for distant emitters can be amplified by the same poor geometry that inflates GDOP. Careful survey of fixed installations, GPS positioning of mobile sensors, and *self-calibration*, jointly estimating sensor positions and emitter locations from emitters of opportunity at known or constrained locations, are the standard responses. A full error budget must include sensor-position uncertainty alongside timing and TDE error; in well-synchronized systems it is often the next-largest term. + +******************* +Advanced Topics +******************* + +Joint TDOA/FDOA Estimation +================================== + +When the emitter, the sensors, or both are *moving*, the relative motion imparts a Doppler shift that differs between sensors, a **Frequency Difference of Arrival** (FDOA). FDOA carries information about the emitter's *velocity* and, crucially, adds an independent geometric constraint that improves position observability, especially for the difficult far-field and few-sensor cases where TDOA alone is poorly conditioned. TDOA and FDOA are estimated jointly by maximizing the **Complex Ambiguity Function** (CAF) over both delay and frequency offset: + +.. math:: + + A(\tau,\nu) = \int_0^T x_i(t)\, x_j^{*}(t-\tau)\, e^{-j2\pi \nu t}\, dt, + +whose two-dimensional peak gives :math:`(\hat\tau_{ij},\hat\nu_{ij})` simultaneously. The CAF generalizes the cross-correlation technique above by adding a frequency-search dimension, at correspondingly higher computational cost. Joint TDOA/FDOA processing is the backbone of satellite and airborne geolocation of radio emitters, where a single pair of moving platforms can localize a stationary emitter from the combined delay and Doppler constraints. diff --git a/figure-generating-scripts/tdoa.py b/figure-generating-scripts/tdoa.py new file mode 100644 index 00000000..47012146 --- /dev/null +++ b/figure-generating-scripts/tdoa.py @@ -0,0 +1,231 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.lines import Line2D +from itertools import combinations +from scipy.signal import firwin, lfilter + +sample_rate = 50e6 +c = 3e8 # speed of light [m/s] +snr_db = 10 # SNR of the received signal at each receiver [dB] +tx_len_samples = 1000 # samples to transmit +rx_positions = np.array([ + [65, 229], # Rx0 + [676, 123], # Rx1 + [153, 543], # Rx2 +]) +num_rx = rx_positions.shape[0] +tx_position = np.array([153, 355]) +pairs = list(combinations(range(num_rx), 2)) # For 3 receivers it's (Rx0,Rx1), (Rx0,Rx2), (Rx1,Rx2) -> 3 pairs + +# For the tx signal itself it's arbitrary, although bandwidth matters, we'll transmit band-limited noise +bandwidth = 20e6 +taps = firwin(numtaps=129, cutoff=bandwidth / 2, fs=sample_rate) +tx_signal = lfilter(taps, 1.0, np.random.randn(tx_len_samples) + 1j * np.random.randn(tx_len_samples)) + +# Simulate what each receiver records +true_distances = np.linalg.norm(rx_positions - tx_position, axis=1) +true_delays = true_distances / c +unknown_tx_time = 1.234e-5 # seconds. arbitray, unknown to receivers and we wont use it in any TDOA calcs + +# Calc the actual TDOAs to act as ground truth +for k, (a, b) in enumerate(pairs): + true_rd = true_distances[b] - true_distances[a] + +# Figure out how many samples we have to simulate +total_delay_samples = (unknown_tx_time + true_delays.max()) * sample_rate +buffer_len = tx_len_samples + int(np.ceil(total_delay_samples)) + 10 + +# Taken from Synchronization chapter +def frac_delay_filter(delay): # delay is in samples, but it can (and will be) not an integer + N = 21 # number of taps, keep this odd + n = np.arange(-(N-1)//2, N//2+1) # -10,-9,...,0,...,9,10 + h = np.sinc(n - delay) # calc filter taps + h *= np.hamming(N) # window the filter to make sure it decays to 0 on both sides + h /= np.sum(h) # normalize to get unity gain, we don't want to change the amplitude/power + return h + +# Simulate the delayed signal being received by each sensor +rx_signals = np.zeros((num_rx, buffer_len), dtype=complex) +for i in range(num_rx): + tau = unknown_tx_time + true_delays[i] # absolute delay at this Rx, in seconds + tau_samples = tau * sample_rate + tau_integer_samps = int(np.round(tau_samples)) + tau_frac_samps = tau_samples - tau_integer_samps + rx = np.zeros(buffer_len, dtype=complex) + rx[tau_integer_samps:tau_integer_samps+tx_len_samples] = tx_signal + frac_delay_i = frac_delay_filter(tau_frac_samps) + rx = np.convolve(rx, frac_delay_i, "same") + + # Each receiver adds its own thermal noise, scaled to hit the SNR set at the top + signal_power = np.mean(np.abs(tx_signal)**2) + noise_power = signal_power / 10**(snr_db / 10) + noise = np.sqrt(noise_power / 2) * (np.random.randn(buffer_len) + 1j * np.random.randn(buffer_len)) + rx_signals[i] = rx + noise + +# Estimate the TDOAs using a normal cross-correlation +range_diff = np.zeros(len(pairs)) # meters +for k, (a, b) in enumerate(pairs): + xcorr = np.correlate(rx_signals[b], rx_signals[a], mode='full') + peak_lag = np.argmax(np.abs(xcorr)) - (buffer_len - 1) # 'full' puts zero lag at index buffer_len-1 + range_diff[k] = (peak_lag / sample_rate) * c # meters + +# FIGURE 1: the integer-only result. +# Precompute the distance from each receiver to every grid point, this will get used in the contour plot +grid_x = np.linspace(-200, 800, 400) +grid_y = np.linspace(-200, 800, 400) +GX, GY = np.meshgrid(grid_x, grid_y) +rx_dist = [] +for i in range(num_rx): + rx_dist.append(np.sqrt((GX - rx_positions[i, 0])**2 + (GY - rx_positions[i, 1])**2)) +fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6)) +# Left: the hyperbolas, which all cross at the transmitter. +hyperbola_handles = [] +pair_colors = ['tab:blue', 'tab:orange', 'tab:green'] +for k, (a, b) in enumerate(pairs): + # the next line is what calculates the hyperbola, note levels=[0] means we're making a contour map but only one level, specifically the level where the difference is zero + ax1.contour(GX, GY, (rx_dist[b] - rx_dist[a]) - range_diff[k], levels=[0], colors=pair_colors[k], linestyles='--') + hyperbola_handles.append(Line2D([0], [0], color=pair_colors[k], linestyle='--', label=f'Rx{a}-Rx{b}')) +ax1.scatter(rx_positions[:, 0], rx_positions[:, 1], c='tab:blue', marker='^', s=120, edgecolors='k', label='Receivers', zorder=5) +for i in range(num_rx): + ax1.annotate(f'Rx{i}', rx_positions[i], textcoords='offset points', xytext=(8, 8), fontweight='bold', zorder=6) +ax1.scatter(*tx_position, c='red', marker='*', s=300, edgecolors='k', label='True Tx', zorder=5) +ax1.set_xlim(grid_x[0], grid_x[-1]); ax1.set_ylim(grid_y[0], grid_y[-1]) +ax1.set_xlabel('x [m]'); ax1.set_ylabel('y [m]') +ax1.set_title('TDOA hyperbolas') +ax1.legend(handles=ax1.get_legend_handles_labels()[0] + hyperbola_handles, loc='upper right') +ax1.set_aspect('equal') +# Cross-correlation of one pair, integer only +a, b = pairs[1] # the (Rx0, Rx2) pair +xcorr = np.abs(np.correlate(rx_signals[b], rx_signals[a], mode='full')) +lags = np.arange(xcorr.size) - (buffer_len - 1) +peak = lags[np.argmax(xcorr)] +ax2.plot(lags, xcorr, 'o-', markersize=5, label='correlation samples') +ax2.axvline(peak, color='red', linestyle='--', label=f'integer peak = {peak} samples') +ax2.set_xlim(peak - 6, peak + 6) +ax2.set_xlabel('lag [samples]'); ax2.set_ylabel('|cross-correlation|') +ax2.set_title(f'Cross-correlation of Rx{b} vs Rx{a}') +ax2.legend() +ax2.grid() +#fig1.savefig('../_images/tdoa_python_integer.svg', bbox_inches='tight') +fig1.tight_layout() + +# Subsample TDOA calc using a freq domain cross-correlation that was padded as a way to interpolate +U = 16 # correlation upsampling factor +half = (buffer_len + 1) // 2 # number of DC + positive-frequency bins +range_diff = np.zeros(len(pairs)) # meters +for k, (a, b) in enumerate(pairs): + # Cross-correlation in the frequency domain + X = np.conj(np.fft.fft(rx_signals[a])) * np.fft.fft(rx_signals[b]) + + # Insert zeros in the high-frequency MIDDLE: DC + positive freqs at the front, negative freqs at the back, so it stays a valid FFT layout. + X_padded = np.zeros(U * buffer_len, dtype=complex) + X_padded[:half] = X[:half] + X_padded[U * buffer_len - (buffer_len - half):] = X[half:] + + # Now IFFT to finish the crosscorrelation + xcorr = np.abs(np.fft.ifft(X_padded)) * U + + # Peak index -> signed lag; indices past the midpoint are negative lags + peak_idx = np.argmax(xcorr) + if peak_idx > U * buffer_len // 2: + peak_idx -= U * buffer_len + peak_lag = peak_idx / U # sub-sample lag, +ve => Rx_b farther + range_diff[k] = (peak_lag / sample_rate) * c # meters + +print("METHOD 2 (sub-sample, zero-padded FFT)") +print(" Pair | true range diff [m] | measured range diff [m]") +for k, (a, b) in enumerate(pairs): + true_rd = true_distances[b] - true_distances[a] + print(f"Rx{b}-Rx{a} | {true_rd:9.1f} | {range_diff[k]:9.1f}") + +# 8. FIGURE 2: the sub-sample result, same layout as Figure 1. +fig2, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6)) + +# Left: the hyperbolas from the refined range differences. +hyperbola_handles = [] +for k, (a, b) in enumerate(pairs): + ax1.contour(GX, GY, (rx_dist[b] - rx_dist[a]) - range_diff[k], levels=[0], + colors=pair_colors[k], linewidths=1.5, linestyles='--') + hyperbola_handles.append(Line2D([0], [0], color=pair_colors[k], + linestyle='--', label=f'Rx{a}-Rx{b}')) +ax1.scatter(rx_positions[:, 0], rx_positions[:, 1], c='tab:blue', marker='^', + s=120, edgecolors='k', label='Receivers', zorder=5) +for i in range(num_rx): + ax1.annotate(f'Rx{i}', rx_positions[i], textcoords='offset points', + xytext=(8, 8), fontweight='bold', zorder=6) +ax1.scatter(*tx_position, c='red', marker='*', s=300, edgecolors='k', + label='True Tx', zorder=5) +ax1.set_xlim(grid_x[0], grid_x[-1]); ax1.set_ylim(grid_y[0], grid_y[-1]) +ax1.set_xlabel('x [m]'); ax1.set_ylabel('y [m]') +ax1.set_title('TDOA hyperbolas') +ax1.legend(handles=ax1.get_legend_handles_labels()[0] + hyperbola_handles, loc='upper right') +ax1.set_aspect('equal') + +# Right: coarse (1 sample/lag) correlation as dots vs the U-times upsampled +# correlation as a smooth curve, so the sub-sample shift is visible. +a, b = pairs[1] # the (Rx0, Rx2) pair +X = np.conj(np.fft.fft(rx_signals[a])) * np.fft.fft(rx_signals[b]) +cc_coarse = np.abs(np.fft.ifft(X)) +X_padded = np.zeros(U * buffer_len, dtype=complex) +X_padded[:half] = X[:half] +X_padded[U * buffer_len - (buffer_len - half):] = X[half:] +cc_fine = np.abs(np.fft.ifft(X_padded)) * U +lags_coarse = np.where(np.arange(buffer_len) <= buffer_len // 2, np.arange(buffer_len), np.arange(buffer_len) - buffer_len) +lags_fine = np.arange(U * buffer_len) / U +lags_fine = np.where(lags_fine <= buffer_len / 2, lags_fine, lags_fine - buffer_len) +peak = lags_coarse[np.argmax(cc_coarse)] +subsample_peak = lags_fine[np.argmax(cc_fine)] +# The lag axes wrap from + back to - partway through, so sort before plotting, +# otherwise the connecting line jumps across the figure. +order_c = np.argsort(lags_coarse) +order_f = np.argsort(lags_fine) +ax2.plot(lags_coarse[order_c], cc_coarse[order_c], 'o', markersize=6, label='coarse (1 sample/lag)') +ax2.plot(lags_fine[order_f], cc_fine[order_f], '.-', label=f'{U}x interpolation') +ax2.axvline(subsample_peak, color='red', linestyle='--', + label=f'sub-sample peak = {subsample_peak:.3f}') +ax2.set_xlim(peak - 6, peak + 6) +ax2.set_xlabel('lag [samples]'); ax2.set_ylabel('|cross-correlation|') +ax2.set_title(f'Cross-correlation of Rx{b} vs Rx{a}') +ax2.legend() +ax2.grid() +fig2.tight_layout() +#fig2.savefig('../_images/tdoa_python_subsample.svg', bbox_inches='tight') + + +# Heatmap portion: evaluate the TDOA cost on the grid and display it under the ax1 map +cost = np.zeros_like(GX) +for k, (a, b) in enumerate(pairs): + cost += ((rx_dist[b] - rx_dist[a]) - range_diff[k])**2 # squared mismatch for this pair, summed over pairs + +# The best estimate is simply the grid cell with the lowest cost +iy, ix = np.unravel_index(np.argmin(cost), cost.shape) +emitter_grid = np.array([grid_x[ix], grid_y[iy]]) +print("Grid estimate:", emitter_grid) # ~[153, 355] + +fig3, ax1 = plt.subplots(1, 1, figsize=(7, 6)) +# Invert the cost into a likelihood-style surface so higher (brighter) = more likely emitter location +likelihood = -np.log10(cost + 1e-9) +# Cap the colormap a bit below the peak so the bright region around the emitter spreads out and is easier to see +vmax = likelihood.min() + 0.5 * (likelihood.max() - likelihood.min()) # 0.5 was adjusted manually to look good +im = ax1.imshow(likelihood, origin='lower', cmap='viridis', vmax=vmax, + extent=[grid_x[0], grid_x[-1], grid_y[0], grid_y[-1]]) +fig3.colorbar(im, ax=ax1, label='likelihood (higher = more likely)') +# Overlay the hyperbolas, receivers, true Tx, and the grid-search estimate +hyperbola_handles = [] +for k, (a, b) in enumerate(pairs): + ax1.contour(GX, GY, (rx_dist[b] - rx_dist[a]) - range_diff[k], levels=[0], + colors=pair_colors[k], linewidths=1.5, linestyles='--') + hyperbola_handles.append(Line2D([0], [0], color=pair_colors[k], linestyle='--', label=f'Rx{a}-Rx{b}')) +ax1.scatter(rx_positions[:, 0], rx_positions[:, 1], c='tab:cyan', marker='^', s=120, edgecolors='k', label='Receivers', zorder=5) +for i in range(num_rx): + ax1.annotate(f'Rx{i}', rx_positions[i], textcoords='offset points', xytext=(8, 8), color='w', fontweight='bold', zorder=6) +ax1.scatter(*tx_position, c='red', marker='*', s=300, edgecolors='k', label='True Tx', zorder=5) +#ax1.scatter(*emitter_grid, c='white', marker='x', s=120, linewidths=2, label='Grid estimate', zorder=6) +ax1.set_xlim(grid_x[0], grid_x[-1]); ax1.set_ylim(grid_y[0], grid_y[-1]) +ax1.set_xlabel('x [m]'); ax1.set_ylabel('y [m]') +ax1.legend(handles=ax1.get_legend_handles_labels()[0] + hyperbola_handles, loc='upper right') +ax1.set_aspect('equal') +fig3.tight_layout() +fig3.savefig('../_images/tdoa_python_heatmap.svg', bbox_inches='tight') + +plt.show() diff --git a/figure-generating-scripts/tdoa_cramer_rao.py b/figure-generating-scripts/tdoa_cramer_rao.py new file mode 100644 index 00000000..3f2390a6 --- /dev/null +++ b/figure-generating-scripts/tdoa_cramer_rao.py @@ -0,0 +1,38 @@ +import numpy as np +import matplotlib.pyplot as plt + +# This script generates a figure showing the Cramer-Rao Lower Bound (CRLB) on +# time-delay estimation accuracy as a function of SNR, for three different +# signal bandwidths. The bound on the variance of any unbiased delay estimate is +# +# var(tau_hat) >= 1 / (8 * pi^2 * beta^2 * T * gamma) +# +# where beta is the RMS (Gabor) bandwidth of the signal, T is the integration +# time, and gamma is an effective SNR factor. We convert the resulting delay +# standard deviation into a ranging error in meters (multiplying by the speed of +# light) since that is the more intuitive quantity for localization. + +c = 3e8 # speed of light [m/s] +T = 100e-6 # integration time [s] + +snr_db = np.linspace(0, 30, 200) # SNR sweep [dB] +gamma = 10 ** (snr_db / 10) # effective SNR factor (linear) + +# For band-limited noise that is flat from -B/2 to +B/2, the RMS bandwidth is +# beta = B / sqrt(12). We show three signal bandwidths. +bandwidths = [1e6, 10e6, 50e6] # signal bandwidths [Hz] + +fig, ax = plt.subplots(figsize=(8, 5)) +for B in bandwidths: + beta = B / np.sqrt(12) # RMS bandwidth [Hz] + var_tau = 1.0 / (8 * np.pi**2 * beta**2 * T * gamma) # delay variance [s^2] + range_std = c * np.sqrt(var_tau) # ranging error [m] + ax.semilogy(snr_db, range_std, label=f'{B/1e6:.0f} MHz bandwidth') + +ax.set_xlabel('SNR [dB]') +ax.set_ylabel('Ranging error (CRLB) [m]') +ax.grid(True, which='both', alpha=0.4) +ax.legend() +ax.set_xlim(snr_db[0], snr_db[-1]) +fig.savefig('../_images/tdoa_cramer_rao.svg', bbox_inches='tight') +plt.show() diff --git a/index.rst b/index.rst index d897abd9..ebe6c0f9 100644 --- a/index.rst +++ b/index.rst @@ -35,6 +35,7 @@ content/pyqt content/detection content/fpv_video + content/tdoa content/about_author .. raw:: html diff --git a/spelling_wordlist.txt b/spelling_wordlist.txt index 424acd36..74b4023c 100644 --- a/spelling_wordlist.txt +++ b/spelling_wordlist.txt @@ -331,3 +331,36 @@ Vandermonde eigendecomposition ULA geolocation +observability +timestamping +coplanar +collinear +GDOP +Cramér +Rao +nonlinearity +Jacobians +linearizes +Foy +iteratively +minimizer +Jacobian +linearized +multimodal +overdetermined +Schau +Linearization +initializers +stationarity +reverberant +suboptimal +linearize +foci +Multilateration +subsample +subsampling +linearization +hyperboloid +underdetermined +unitless +linearizing