diff --git a/content-ukraine/detection.rst b/content-ukraine/detection.rst
new file mode 100644
index 00000000..4a3b8c9b
--- /dev/null
+++ b/content-ukraine/detection.rst
@@ -0,0 +1,1128 @@
+.. _detection-chapter:
+
+#####################################################
+Виявлення за допомогою кореляції
+#####################################################
+
+.. raw:: html
+
+ У співавторстві із Семом Брауном
+
+У цьому розділі ми навчимося виявляти наявність сигналів і відновлювати їхню часову прив'язку, обчислюючи взаємну кореляцію прийнятих відліків із заздалегідь відомою частиною сигналу, наприклад преамбулою пакета. Звідси природно випливає проста форма класифікації за допомогою банку кореляторів. Ми введемо основні поняття виявлення сигналів і зосередимося на тому, як вирішити, чи присутній заданий сигнал у шумному середовищі. Водночас розглянемо теорію та практичні методи ухвалення надійних рішень в умовах невизначеності.
+
+****************************************************
+Основи виявлення сигналів і кореляторів
+****************************************************
+
+Виявлення сигналу — це задача визначення, чи є спостережуваний сплеск енергії корисним сигналом, чи лише фоновим шумом.
+
+Складність полягає в тому, що в радіолокаційних і гідроакустичних системах шум присутній усюди. Надто чутливий детектор створює хибні тривоги, а недостатньо чутливий пропускає справжню ціль.
+
+Розв'язання починається з детектора Неймана—Пірсона, який математично знаходить оптимальний компроміс: максимізує ймовірність виявлення сигналу, водночас утримуючи хибні тривоги нижче заданої межі. Детектори CFAR розвивають цю ідею, пристосовуючись до змін рівня шуму. Вони особливо корисні, коли статистика шуму нестаціонарна, тобто шумовий поріг і розподіл змінюються через завади або мінливі умови каналу. Мета полягає в автоматичному регулюванні порога виявлення разом із фоновим шумом зі збереженням заданої частоти хибних тривог. Для цього потрібно оцінювати шумовий поріг у часі.
+
+Коли система визначила, що сигнал присутній, їй ще потрібно знайти точний початок даних. Цифрові пакети LTE, 5G і WiFi починаються з преамбули — відомого повторюваного цифрового шаблону. Корелятор преамбули працює подібно до замка з ключем: ключем є відома приймачу послідовність символів, унікальна для відновлюваного сигналу. Приймач переміщує копію преамбули вздовж вхідного сигналу та для кожної затримки обчислює скалярний добуток, вимірюючи подібність шаблону до прийнятих відліків у кожному положенні. Коли вони точно суміщаються, виникає гострий максимум, який указує, звідки починати читати дані. Складніші варіанти також ураховують частотні зсуви через невелику різницю налаштування телефона й базової станції або через ефект Доплера.
+
+Коли відомий сигнал або преамбула передається каналом, спотвореним лише адитивним білим гаусовим шумом (AWGN), потрібно визначити, чи присутній сигнал. Це найпростіша й фундаментальна задача виявлення.
+
+Функція взаємної кореляції
+##########################
+
+У найпростішому вигляді корелятор обчислює взаємну кореляцію між прийнятим сигналом і шаблоном. Взаємна кореляція — це скалярний добуток двох векторів, коли один із них ковзає вздовж іншого. Якщо ви вже вивчали згортку, операція майже така сама, але другий вектор не перевертається, тому вона трохи простіша. Для комплексних сигналів, із якими ми працюватимемо, один із входів також потрібно комплексно спрягати. У Python це можна реалізувати так:
+
+.. code-block:: python
+
+ def correlate(a, v):
+ n = len(a)
+ m = len(v)
+ result = []
+ for i in range(n - m + 1):
+ s = 0
+ for j in range(m):
+ s += a[i + j] * v[j].conjugate()
+ result.append(s)
+ return result
+
+ # Example usage:
+ a = [1+2j, 2+1j, 3+0j, 4-1j, 5-2j]
+ v = [0+1j, 1+0j, 0.5-0.5j]
+ correlate(a, v)
+
+Зверніть увагу, як ми переміщуємо :code:`a`, комплексно спряжуючи :code:`v`, і що цикл із :code:`j` та :code:`s` насправді лише обчислює скалярний добуток векторів. На щастя, реалізовувати взаємну кореляцію з нуля не потрібно: у Python можна скористатися функцією NumPy :code:`correlate`. У SciPy також є власна версія для експериментів.
+
+Приклад взаємної кореляції в Python
+###################################
+
+Для простого прикладу корелятора спочатку потрібний сигнал із відомою преамбулою, зануреною в шум. Як преамбулу використаємо послідовність Задова—Чу завдяки її чудовим автокореляційним властивостям і поширеності в системах зв'язку. Іншу частину сигналу моделювати не будемо, хоча в більшості систем після відомої преамбули йдуть невідомі дані. Послідовність Задова—Чу можна згенерувати так:
+
+.. code-block:: python
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ N = 839 # Length of Zadoff-Chu sequence
+ u = 25 # Root of ZC sequence
+ t = np.arange(N)
+ zadoff_chu = np.exp(-1j * np.pi * u * t * (t + 1) / N)
+
+Отримана послідовність сама є сигналом. IQ-відліки в :code:`zadoff_chu` представляють комплексний сигнал основної смуги, подібний до багатьох уже розглянутих у підручнику, але не представляють біти. Реалістичний сценарій можна змоделювати, додавши сигнал Задова—Чу у випадковій позиції до довшого потоку AWGN:
+
+.. code-block:: python
+
+ signal_length = 10 * N # overall simulated signal length
+ offset = np.random.randint(N, signal_length - N)
+ print(f"True offset: {offset}")
+ snr_db = -15
+ noise_power = 1 / (2 * (10**(snr_db / 10)))
+ signal = np.sqrt(noise_power/2) * (np.random.randn(signal_length) + 1j * np.random.randn(signal_length))
+ signal[offset:offset+N] += zadoff_chu # place our ZC signal at the random offset
+
+Зауважте, що використовується дуже низьке SNR. Воно настільки мале, що в часовій області послідовність Задова—Чу взагалі не видно. Вона має 839 відліків серед приблизно 8000 змодельованих і настільки глибоко занурена в шум, що неможливо помітити навіть невелике збільшення модуля сигналу.
+
+.. image:: ../_images/detection_basic_1.svg
+ :align: center
+ :target: ../_images/detection_basic_1.svg
+ :alt: Сигнал у часовій області з послідовністю Задова—Чу
+
+Тепер реалізуємо корелятор, обчисливши взаємну кореляцію прийнятого сигналу з відомою послідовністю Задова—Чу за допомогою :code:`np.correlate()`. Припускається, що приймач знає точну використану преамбулу. У попередньому коді :code:`zadoff_chu` спочатку створено для моделювання сигналу, а тепер та сама змінна є шаблоном преамбули на приймачі. Корелятор реалізується одним рядком:
+
+.. code-block:: python
+
+ correlation = np.correlate(signal, zadoff_chu, mode='valid')
+
+Режим :code:`valid` буде пояснено нижче. Також нормалізуємо результат на довжину послідовності та беремо квадрат модуля, щоб отримати потужність, хоча можна було б узяти лише модуль. Головною операцією є саме :code:`np.correlate()`.
+
+.. code-block:: python
+
+ correlation = np.abs(correlation / N)**2 # normalize by N, and take magnitude squared
+
+Побудуємо графік квадрата модуля й позначимо справжній початок послідовності, щоб перевірити роботу корелятора:
+
+.. image:: ../_images/detection_basic_2.svg
+ :align: center
+ :target: ../_images/detection_basic_2.svg
+ :alt: Вихід корелятора
+
+Попри дуже низьке SNR, на виході корелятора видно чіткий максимум саме там, де розміщено послідовність Задова—Чу. Він позначає початок послідовності, тож наступні 839 відліків містять преамбулу. У цьому й полягає сила кореляційного виявлення в поєднанні з довгою преамбулою. Поріг, який визначає, чи є максимум цільовим сигналом або лише шумом, ми ще не задали — поки що лише розглядаємо результат візуально. Решта розділу присвячена автоматизації цього рішення, особливо за змінного шумового порога та фонових завад.
+
+Режими valid, same і full
+#########################
+
+Можливо, ви помітили, що :code:`np.correlate()` і :code:`np.convolve()` підтримують три режими: :code:`valid`, :code:`same` та :code:`full`. Вони визначають довжину вихідного масиву відносно вхідних. Ми використали :code:`valid`, тому вихід містить лише точки, у яких два вхідні масиви повністю перекриваються. Довжина результату становить :code:`len(signal) - len(zadoff_chu) + 1`. У режимі :code:`same` вихід мав би довжину довшого вхідного сигналу. Режим :code:`full` дає повну дискретну лінійну згортку — дещо довший масив довжиною :code:`max(M, N) - min(M, N) + 1`, де :code:`M` і :code:`N` — довжини входів. У радіочастотній обробці згортку часто застосовують як FIR-фільтр, де однакова довжина входу й виходу зручна, тому поширений режим :code:`same`. Для кореляційного виявлення зазвичай потрібний :code:`valid`, адже нас цікавлять лише точки повного перекриття преамбули з прийнятим сигналом, особливо якщо передбачається, що сигнал почався вже після початку приймання.
+
+Детектор Неймана—Пірсона
+########################
+
+Еталонним способом вибору порога для виходу корелятора є детектор Неймана—Пірсона. Теорія дає оптимальне рішення за конкретного обмеження: знаходить поріг, який максимізує ймовірність виявлення :math:`P_{D}` за фіксованої допустимої ймовірності хибної тривоги :math:`P_{FA}`. Простіше кажучи, ви визначаєте прийнятну кількість хибних виявлень, наприклад одну хибну тривогу на годину, а детектор Неймана—Пірсона задає найкращий поріг, щоб виявити якомога більше справжніх сигналів. Для відомої преамбули в AWGN підхід простий: обчислюється кореляція між прийнятим сигналом і відомим шаблоном. Якщо значення перевищує наперед заданий поріг :math:`\tau`, сигнал вважається присутнім; інакше вважається, що є лише шум.
+
+Характеристики детектора, вимірювані :math:`P_{D}` і :math:`P_{FA}`, залежать від порога :math:`\tau`, SNR та довжини преамбули :math:`L`. Імовірність хибної тривоги залежить від порога й дисперсії шуму :math:`\sigma_n^2`:
+
+:math:`P_{FA} = Q\left(\frac{\tau}{\sigma_n}\right)`
+
+Імовірність виявлення є функцією порога, дисперсії шуму та енергії преамбули (:math:`E_s = L \cdot S`, де :math:`S` — середня потужність символу):
+
+:math:`P_{D} = Q\left(\frac{\tau - \sqrt{E_s}}{\sigma_n}\right) = Q\left(\frac{\tau - \sqrt{L \cdot S}}{\sigma_n}\right)`
+
+Тут :math:`Q(x)` — Q-функція, тобто хвостова ймовірність стандартного нормального розподілу: імовірність того, що стандартна нормальна випадкова величина перевищить :math:`x`.
+
+Аналіз характеристик: ROC і залежність Pd від SNR
+##################################################
+
+Для кількісної оцінки кореляційного детектора за наявності шуму використовують дві основні візуалізації: робочу характеристику приймача (Receiver Operating Characteristic, ROC) і графік імовірності виявлення :math:`P_d` залежно від SNR.
+
+ROC показує ймовірність виявлення :math:`P_D` залежно від імовірності хибної тривоги :math:`P_{FA}` за фіксованого SNR. Змінюючи поріг на виході корелятора, ми вибираємо точку на цій кривій, тобто фундаментальний компроміс. Нижчий поріг збільшує :math:`P_D`, бо знаходить більше справжніх сигналів, але водночас збільшує :math:`P_{FA}`, частіше спрацьовуючи на шум. Чим сильніше крива вигинається до верхнього лівого кута, тим кращий детектор. Ідеальний детектор досягає цього кута зі 100% :math:`P_D` і 0% :math:`P_{FA}`, тоді як діагональ відповідає випадковому вгадуванню.
+
+.. image:: ../_images/detection_pd_vs_snr.svg
+ :align: center
+ :target: ../_images/detection_pd_vs_snr.svg
+ :alt: Залежність Pd від SNR і ROC
+
+Рівняння та інтуїція разом показують, що довжина преамбули :math:`L` є критичним параметром проєктування, бо безпосередньо визначає виграш обробки, а отже, якість виявлення. За фіксованих порога та SNR імовірність :math:`P_D` зростає з :math:`L`. Довша преамбула дає змогу накопичити більше енергії сигналу й легше відрізнити його від фонового шуму. Це поліпшення називають виграшем обробки й зазвичай вимірюють у децибелах як :math:`10\log_{10}(L)`. Воно критично важливе для виявлення слабких сигналів, які інакше було б пропущено. Інтегруючи енергію більшої кількості відліків, можна витягнути сигнал із шуму, навіть коли він лежить нижче шумового порога. Хороший реальний приклад — GPS, приймач якого має відновлювати надзвичайно слабкі сигнали з відомою кодовою структурою.
+
+****************************************************
+Приклад: виявлення GPS нижче шумового порога
+****************************************************
+
+Короткий вступ до сигналів GPS
+##############################
+
+Станом на березень 2026 року американське сузір'я GPS містить 31 діючий супутник на середній навколоземній орбіті (MEO), кожен із яких двічі на добу облітає Землю. Усі супутники передають на одній несній частоті 1575,42 МГц, яка називається L1. На поверхні Землі сигнал надзвичайно слабкий і лежить далеко нижче шумового порога. Ортогональність супутників забезпечується унікальним для кожного 1023-чиповим псевдовипадковим шумовим кодом (PRN), що називається C/A-кодом; тому сигнал також позначають L1 C/A. C/A-коди є кодами Голда й спеціально побудовані так, щоб будь-яка їхня пара була майже ортогональною: взаємна кореляція кодів двох супутників майже нульова. C/A-код передається зі швидкістю 1,023 мільйона чипів за секунду й має лише 1023 чипи, тому повторюється кожну 1 мс. Поверх повторюваного коду кожен супутник повільно модулює навігаційні дані, зокрема параметри орбіти й поправки годинника, зі швидкістю лише 50 біт/с. Один біт даних охоплює 20 повних повторень коду. Використання окремого коду для кожного передавача називається кодовим розділенням каналів (Code Division Multiple Access, CDMA); ту саму ідею застосовували в мобільному зв'язку 3G.
+
+На боці приймача пошук одного з 31 супутника означає генерування локальної копії його PRN-послідовності й використання корелятора для знаходження початку кодового періоду. Хоча GPS передає безперервно, цей початок можна розглядати як початок пакета або кадру. Точне положення кореляційного максимуму також використовують для оцінювання пройденої сигналом відстані. Маючи такі оцінки для чотирьох або більше супутників, приймач визначає своє положення на Землі трилатерацією. Супутники рухаються відносно вас приблизно зі швидкістю 4 км/с, тому приймач також має перебрати сітку можливих частотних зсувів і знайти найбільший кореляційний максимум. Це двовимірний пошук. Максимальний доплерівський зсув становить приблизно +/-20 кГц (:code:`4e3 / 3e8 * 1.575e9`). Процес повторюється кожну 1 мс, але приймач відстежує затримку та доплерівський зсув, тому не виконує повний пошук щоразу. Початковий пошук супутника називають захопленням, а подальше супроводження сигналу — стеженням. Захоплення обчислювально складніше й може тривати кілька хвилин, якщо приймач починає без попередньої інформації про видимі супутники, доплерівські зсуви або власне положення.
+
+Кореляційний підхід
+###################
+
+Обчислюємо взаємну кореляцію вхідного сигналу — у цьому випадку запису L1 — із локально згенерованою копією коду кожного супутника. Великий максимум означає, що супутник видно, і задає початок кодового періоду 1 мс. Для одночасного пошуку за частотою використовуємо кореляцію на основі ШПФ у частотній області. Вона дає змогу ефективно перевіряти багато частотних зсувів, переміщуючи відліки ШПФ локальної копії коду. Нарешті, накопичуємо квадрат модуля кореляції для кількох блоків по 1 мс, щоб підвищити SNR. Це називають некогерентним інтегруванням; воно допомагає виявляти GPS-сигнали, прийняті нижче шумового порога. Для нормалізації порівнюємо з порогом відношення виходу корелятора до середньої потужності кореляції за всіма затримками.
+
+Приклад запису
+##############
+
+Скористаємося прикладом запису GPS від Даніеля Естевеса, який можна `завантажити тут `_. Це файл комплексних float32 із частотою дискретизації 4 МГц і центральною частотою 1575,42 МГц.
+
+Нижче показано спектрограму запису. На ній майже нічого не видно, а вертикальна лінія не є справжнім GPS-сигналом — імовірно, це вузькосмугова завада. Сигнали GPS L1 використовують швидкість 1,023 Мчип/с із дуже повільними даними поверх коду, тому займають приблизно 2 МГц. Проте на спектрограмі цієї смуги просто не видно. Це наочний приклад приймання GPS нижче шумового порога й причини, чому для його пошуку потрібне кореляційне виявлення.
+
+.. image:: ../_images/detection_gps_spectrogram.svg
+ :align: center
+ :target: ../_images/detection_gps_spectrogram.svg
+ :alt: Спектрограма запису GPS L1
+
+Для зацікавлених: цей запис є невеликою частиною значно більшого файла на `IQEngine `_ у каталозі :code:`estevez/GPS and other GNSS`; знайдіть запис :code:`GPS-L1-2022-03-27`. На IQEngine це файл int16 у форматі SigMF.
+
+Приклад Python
+##############
+
+Змініть :code:`filename` відповідно до розташування завантаженого IQ-файла. Параметр :code:`num_integrations` визначає обсяг запису, який читається й обробляється. Загальна тривалість дорівнює цьому числу, помноженому на 1 мс; для скороченого запису максимальне значення становить 10.
+
+.. code-block:: python
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+
+ filename = "GPS_L1_recording_10ms_4MHz_cf32.iq"
+ sample_rate = 4e6
+ chip_rate = 1023000 # chips / sec (part of the GPS spec)
+ num_chips = 1023 # chips per C/A code period
+ samples_per_code = int(round(sample_rate / chip_rate * num_chips)) # Exact number of samples in one 1 ms code period at 4 MHz
+ doppler_min_hz = -5e3 # GPS Doppler ≈ ±4 kHz for stationary receiver
+ doppler_max_hz = 5e3
+ doppler_step_hz = 500 # good enough for a coarse search
+ num_integrations = 10 # non-coherent power integrations (so 10 ms total), determines how much of the IQ recording we read in and process!
+ detection_thresh_dB = 14.0 # Peak-to-mean ratio (PMR) threshold in dB to declare a detection, GPS C/A signals are typically 14–20 dB PMR above threshold with 10ms of integration
+ gps_svs = list(range(1, 33)) # 1–32
+
+ ##### C/A Code Generation #####
+ # The GPS C/A code is a Gold code formed by XOR-ing two 10-stage maximal-length
+ # shift registers (G1 and G2). G2 is effectively delayed by a satellite-
+ # specific number of chips before the XOR
+ # Reference: IS-GPS-200, Table 3-Ia
+ G2_DELAY = [ # G2 phase delay (chips) for gps_svs 1–32
+ 5, 6, 7, 8, 17, 18, 139, 140, # 1– 8
+ 141, 251, 252, 254, 255, 256, 257, 258, # 9–16
+ 469, 470, 471, 472, 473, 474, 509, 512, # 17–24
+ 513, 514, 515, 516, 859, 860, 861, 862, # 25–32
+ ]
+
+ """G1 LFSR: polynomial x^10 + x^3 + 1, all-ones init, output at stage 10."""
+ reg = np.ones(10, dtype=np.int8)
+ G1 = np.empty(num_chips, dtype=np.int8)
+ for i in range(num_chips):
+ G1[i] = reg[9]
+ fb = reg[2] ^ reg[9] # stages 3 and 10 (0-indexed: 2 and 9)
+ reg = np.roll(reg, 1)
+ reg[0] = fb
+
+ """G2 LFSR: polynomial x^10+x^9+x^8+x^6+x^3+x^2+1, all-ones init."""
+ reg = np.ones(10, dtype=np.int8)
+ G2 = np.empty(num_chips, dtype=np.int8)
+ for i in range(num_chips):
+ G2[i] = reg[9]
+ fb = reg[1]^reg[2]^reg[5]^reg[7]^reg[8]^reg[9] # taps 2,3,6,8,9,10
+ reg = np.roll(reg, 1)
+ reg[0] = fb
+
+ # 1023-chip C/A PRN code for SV sv (1-32) as float32, 1's and -1's, so BPSK
+ def make_prn(sv: int) -> np.ndarray:
+ g2_delayed = np.roll(G2, G2_DELAY[sv - 1])
+ bits = G1 ^ g2_delayed # {0, 1}
+ return (1 - 2 * bits).astype(np.float32) # BPSK: {+1, −1}
+
+ def upsample_prn(sv: int) -> np.ndarray:
+ """Nearest-neighbour upsample 1023-chip C/A code → samples_per_code samples."""
+ code = make_prn(sv)
+ idx = (np.arange(samples_per_code) * num_chips / samples_per_code).astype(int)
+ return code[idx]
+
+ # Pre-compute template signals - conjugate FFTs of all upsampled PRN codes
+ template_signals = {sv: np.conj(np.fft.fft(upsample_prn(sv))) for sv in gps_svs}
+
+ # Read in IQ file
+ n_needed = samples_per_code * num_integrations
+ iq = np.fromfile(filename, dtype=np.complex64, count=n_needed)
+ # For the full version from IQEngine use the following instead
+ #iq = np.fromfile(filename, dtype=np.int16, count=n_needed * 2)
+ #iq = (iq[0::2] + 1j * iq[1::2]).astype(np.complex64)
+
+ # Search each satellite across Doppler and code phase
+ results = []
+ detected = []
+ print(f" {'SV':>3} {'Doppler (Hz)':>13} {'Phase (chips)':>14}"
+ f" {'Phase (samp)':>13} {'Delay (µs)':>11} {'PMR (dB)':>9}")
+ doppler_bins = np.arange(doppler_min_hz, doppler_max_hz + doppler_step_hz, doppler_step_hz)
+ for sv in gps_svs:
+ corr_map = np.zeros((len(doppler_bins), samples_per_code))
+ n_total = samples_per_code * num_integrations
+ for di, f_d in enumerate(doppler_bins):
+ t = np.arange(n_total) / sample_rate # Time vector
+ mixed = iq[:n_total] * np.exp(-2j*np.pi*float(f_d)*t) # Apply the frequency shift
+
+ # Accumulate squared correlation magnitude non-coherently
+ for k in range(num_integrations):
+ blk = mixed[k * samples_per_code:(k + 1) * samples_per_code]
+ sig_fft = np.fft.fft(blk)
+ corr = np.fft.ifft(sig_fft * template_signals[sv]) # Frequency-domain correlation
+ corr_map[di] += np.abs(corr)**2
+
+ # Normalize by the mean and convert to dB
+ peak_val = float(np.max(corr_map))
+ mean_val = float(np.mean(corr_map))
+ pmr_db = 10.0 * np.log10(peak_val / mean_val)
+
+ peak_idx = np.unravel_index(np.argmax(corr_map), corr_map.shape)
+ best_doppler_hz = float(doppler_bins[peak_idx[0]])
+ best_phase_samp = int(peak_idx[1])
+ best_phase_chips = best_phase_samp * num_chips / samples_per_code
+
+ r = {
+ "sv": sv,
+ "detected": pmr_db >= detection_thresh_dB,
+ "doppler_hz": best_doppler_hz,
+ "code_phase_samp": best_phase_samp, # sample offset = "start of packet"
+ "code_phase_chip": best_phase_chips,
+ "pmr_db": pmr_db,
+ "corr_map": corr_map,
+ "doppler_bins": doppler_bins,
+ }
+ results.append(r)
+
+ # Print the result row
+ delay_us = r['code_phase_samp'] / sample_rate * 1e6
+ flag = " ← DETECTED" if r['detected'] else ""
+ print(f" {sv:>3} {r['doppler_hz']:>+13.0f} {r['code_phase_chip']:>14.2f}"
+ f" {r['code_phase_samp']:>13d} {delay_us:>11.3f} {r['pmr_db']:>9.1f}{flag}")
+
+Цей код має дати такий результат:
+
+.. code-block::
+
+ SV Doppler (Hz) Phase (chips) Phase (samp) Delay (µs) PMR (dB)
+ 1 -3000 757.79 2963 740.750 5.6
+ 2 +1500 264.19 1033 258.250 9.1
+ 3 -2000 316.62 1238 309.500 5.8
+ 4 +5000 577.48 2258 564.500 5.0
+ 5 +1000 64.96 254 63.500 5.3
+ 6 +1500 511.76 2001 500.250 5.0
+ 7 -4000 763.41 2985 746.250 5.0
+ 8 +3500 961.62 3760 940.000 5.4
+ 9 +3500 118.67 464 116.000 4.9
+ 10 +0 890.52 3482 870.500 5.4
+ 11 +2500 837.33 3274 818.500 14.6 ← DETECTED
+ 12 -500 871.60 3408 852.000 16.4 ← DETECTED
+ 13 +1000 137.85 539 134.750 5.9
+ 14 +2500 287.72 1125 281.250 5.0
+ 15 -5000 908.68 3553 888.250 5.3
+ 16 +1500 292.58 1144 286.000 5.9
+ 17 +500 994.61 3889 972.250 5.3
+ 18 +4500 1005.61 3932 983.000 5.4
+ 19 +5000 588.48 2301 575.250 5.0
+ 20 +0 768.53 3005 751.250 5.4
+ 21 -3000 749.60 2931 732.750 5.0
+ 22 +2500 558.05 2182 545.500 14.4 ← DETECTED
+ 23 -5000 390.02 1525 381.250 5.3
+ 24 +2500 955.48 3736 934.000 5.9
+ 25 +1500 597.94 2338 584.500 15.5 ← DETECTED
+ 26 -1500 239.89 938 234.500 6.2
+ 27 -2500 488.74 1911 477.750 4.7
+ 28 +3000 858.81 3358 839.500 5.2
+ 29 -4000 998.70 3905 976.250 5.2
+ 30 -2000 937.58 3666 916.500 5.2
+ 31 +5000 463.42 1812 453.000 15.9 ← DETECTED
+ 32 +1000 342.45 1339 334.750 16.2 ← DETECTED
+
+Як бачимо, виявлено шість супутників. Хоча поріг дорівнював 14,0, зі списку легко визначити, що більшість інших супутників не перебувала в зоні видимості. Винятком є SV-2, який, імовірно, було видно, але він трохи не досяг порога. Для охочих перевірити: запис зроблено 2022-03-27 о 11:32:04 десь в Іспанії.
+
+Побудова графіків
+#################
+
+Побудуймо результати для супутника 11 — першого з виявлених. Перший графік є двовимірною картою кореляції за доплерівським зсувом і часом/затримкою. Другий — зріз карти в найкращому доплерівському відліку, який показує потужність кореляції в часі, як у попередньому розділі.
+
+.. code-block:: python
+
+ # Plotting
+ sv = 11 # we detected 11, 12, 22, 25, 31, 32 although try looking at one we didnt find as well!
+ r = results[sv - 1] # print the dict of results for this SV to see what we got
+ cmap = r['corr_map'] # 2-D array of correlation power vs Doppler and code phase
+ d_bins = r['doppler_bins'] # Doppler bins corresponding
+ chips_axis = np.arange(samples_per_code) * num_chips / samples_per_code
+
+ # 2-D Doppler × code-phase map
+ plt.figure(0, figsize=(10, 6))
+ im = plt.pcolormesh(chips_axis, d_bins, cmap, shading='auto', cmap='viridis')
+ plt.xlabel("Code Phase (chips)")
+ plt.ylabel("Doppler (Hz)")
+ plt.title(f"SV {sv} — 2-D Acquisition Map (PMR = {r['pmr_db']:.1f} dB)")
+ plt.legend(fontsize=8, loc='upper right')
+ plt.colorbar(im, label="Correlation Power")
+
+ # Code-phase slice at the best Doppler
+ best_di = int(np.argmin(np.abs(d_bins - r['doppler_hz'])))
+ plt.figure(1, figsize=(10, 6))
+ plt.plot(chips_axis, cmap[best_di], lw=1, color='steelblue')
+ plt.xlabel("Code Phase (chips)")
+ plt.ylabel("Correlation Power")
+ plt.title(f"SV {sv} — Code-Phase Slice (Doppler = {r['doppler_hz']:+.0f} Hz)")
+ plt.legend(fontsize=8)
+ plt.grid(True, alpha=0.3)
+
+ plt.show()
+
+.. image:: ../_images/detection_gps_2d_map.png
+ :align: center
+ :width: 700px
+ :alt: Двовимірна карта захоплення
+
+.. image:: ../_images/detection_gps_code_phase_slice.svg
+ :align: center
+ :target: ../_images/detection_gps_code_phase_slice.svg
+ :alt: Зріз за кодовою фазою
+
+Трилатерацію тут не розглядатимемо, але саме точне положення цього максимуму зрештою дає GPS-приймачу змогу визначити відстань до супутника. Поєднавши таку інформацію від чотирьох або більше супутників, він визначає своє положення на Землі.
+
+****************************************************
+Детектори CFAR у мінливому середовищі
+****************************************************
+
+Детектор Неймана—Пірсона оптимальний за сталого рівня шуму, але реальні умови рідко бувають настільки стабільними. У динамічному середовищі — наприклад, коли радар супроводжує літак крізь дощ або бездротовий приймач працює в переповненому місті — рівні фонового шуму й завад постійно змінюються. Саме тут потрібний детектор зі сталою частотою хибних тривог (Constant False Alarm Rate, CFAR).
+
+CFAR є основним інструментом систем, де непередбачуваний фон не дає підтримувати фіксований поріг:
+
+- Радіолокація та гідроакустика виявляють цілі, як-от літаки чи підводні човни, на тлі «місцевих предметів» — відбиттів від хвиль, дощу або землі, які змінюються під час руху сенсора.
+- Системи бездротового зв'язку, зокрема когнітивне радіо та LTE/5G, застосовують CFAR для пошуку вільного спектра або вхідних пакетів за імпульсних і непередбачуваних завад від інших пристроїв.
+- У медичній візуалізації CFAR допомагає автоматичному аналізу УЗД або МРТ відрізняти справжні особливості тканин від змінного електронного шуму.
+
+Літера «C» у CFAR означає Constant («сталий»), бо мета полягає в утриманні ймовірності хибної тривоги :math:`P_{FA}` на сталому передбачуваному рівні.
+
+Для встановлення порога потрібно прийняти статистичну модель шуму, тобто його розподіл. У простому AWGN шум має гаусів розподіл. Радіолокаційні відбиття від місцевих предметів можуть натомість мати розподіл Релея або Вейбулла. Якщо модель неправильна, :math:`P_{FA}` «дрейфуватиме», через що система або перестане бачити цілі, або буде перевантажена хибними спрацьовуваннями.
+
+Замість жорстко заданого значення CFAR оцінює потужність шуму в локальному «околі» сигналу й множить її на коефіцієнт масштабування :math:`T`, отриманий із бажаного :math:`P_{FA}`. Тому разом зі зростанням шумового порога зростає й поріг виявлення.
+
+Хибні тривоги для однієї затримки та всієї системи
+##################################################
+
+Це важлива відмінність, яку початківці часто пропускають. Під час пошуку преамбули зазвичай виконується ковзна кореляція, а поріг щосекунди перевіряється в тисячах часових зсувів, або затримок.
+
+:math:`P_{FA}` **для однієї затримки** — імовірність того, що одна конкретна перевірка кореляції дасть хибну тривогу. Якщо розрахунок задає :math:`P_{FA}=0.001`, кожна окрема затримка має шанс 1 до 1000 стати «примарним» сигналом.
+
+:math:`P_{FA}` **рівня системи (глобальна)** — імовірність того, що система створить хоча б одну хибну тривогу протягом усього вікна пошуку, наприклад серед 2048 затримок.
+
+Якщо :math:`P_{FA}` для однієї затримки дорівнює :math:`p`, то ймовірність принаймні однієї хибної тривоги серед :math:`N` затримок приблизно дорівнює :math:`1-(1-p)^{N}`.
+
+Отже, для 1000 затримок і :math:`P_{FA}=0.001` на одну затримку система насправді видаватиме хибну тривогу майже в 63% пошуків. Щоб глобальна частота хибних тривог залишалася низькою, :math:`P_{FA}` для однієї затримки має бути надзвичайно малим.
+
+Приклад Python
+##############
+
+Щоб поекспериментувати з власним CFAR, спочатку змоделюємо повторювані QPSK-пакети з відомою преамбулою, які проходять каналом зі змінним у часі шумовим порогом. Потім реалізуємо простий CFAR з усередненням комірок (Cell-Averaging CFAR, CA-CFAR), щоб знаходити преамбули в прийнятому сигналі. Наступний код генерує прийнятий сигнал:
+
+.. code-block:: python
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ from scipy.signal import correlate
+
+ def generate_qpsk_packets(num_packets, sps, preamble):
+ """Generates repeating QPSK packets with gaps and varying noise."""
+ qpsk_map = np.array([1+1j, -1+1j, -1-1j, 1-1j]) / np.sqrt(2)
+ data_len = 200
+ gap_len = 100
+ full_signal = []
+
+ # Precompute the upsampled preamble for correlation
+ upsampled_preamble = np.repeat(preamble, sps)
+
+ for _ in range(num_packets):
+ data = qpsk_map[np.random.randint(0, 4, data_len)]
+ packet = np.concatenate([preamble, data])
+ full_signal.extend(np.repeat(packet, sps))
+ full_signal.extend(np.zeros(gap_len * sps))
+
+ return np.array(full_signal), upsampled_preamble
+
+ # Simulation parameters
+ sps = 4
+ preamble_syms = np.array([1+1j, 1+1j, -1-1j, -1-1j, 1-1j, -1+1j]) / np.sqrt(2)
+ tx_signal, ref_preamble = generate_qpsk_packets(5, sps, preamble_syms)
+
+ # Time-varying noise floor
+ t = np.arange(len(tx_signal))
+ noise_env = 0.05 + 0.3 * np.sin(2 * np.pi * 0.0003 * t)**2
+ noise = (np.random.randn(len(tx_signal)) + 1j*np.random.randn(len(tx_signal))) * noise_env
+ rx_signal = tx_signal + noise
+
+Перший крок — один раз обчислити кореляцію прийнятого сигналу з відомою преамбулою. На практиці це зазвичай роблять блоками відліків, але поки що обробимо все одним блоком:
+
+.. code-block:: python
+
+ # Correlation spike appears when the reference matches the received segment
+ corr_out = correlate(rx_signal, ref_preamble, mode='same')
+ corr_power = np.abs(corr_out)**2
+
+Тепер реалізуємо CFAR, застосуємо його до виходу корелятора й покажемо результат:
+
+.. code-block:: python
+
+ # CFAR detection on the correlator output
+ def ca_cfar_adaptive(data, num_train, num_guard, pfa):
+ num_cells = len(data)
+ thresholds = np.zeros(num_cells)
+ alpha = num_train * (pfa**(-1/num_train) - 1) # Scaling factor
+ half_window = (num_train + num_guard) // 2
+ guard_half = num_guard // 2
+ for i in range(half_window, num_cells - half_window):
+ # Build the training set around the cell under test (CUT)
+ lagging_win = data[i - half_window : i - guard_half]
+ leading_win = data[i + guard_half + 1 : i + half_window + 1]
+ noise_floor_est = np.mean(np.concatenate([lagging_win, leading_win]))
+ thresholds[i] = alpha * noise_floor_est
+ return thresholds
+
+ # Detect peaks in correlator power
+ cfar_thresholds = ca_cfar_adaptive(corr_power, num_train=60, num_guard=20, pfa=1e-5)
+ detections = np.where(corr_power > cfar_thresholds)[0]
+ # Remove edge detections where the threshold is undefined
+ detections = detections[cfar_thresholds[detections] > 0]
+
+ # Subplot 1: received signal and raw power
+ plt.figure(figsize=(14, 8))
+ plt.subplot(2, 1, 1)
+ plt.plot(np.abs(rx_signal)**2, color='gray', alpha=0.4, label='Rx Signal Power ($|r(t)|^2$)')
+ plt.title("Time-Domain Received Signal")
+ plt.ylabel("Power")
+ plt.legend()
+ plt.grid(True, alpha=0.3)
+
+ # Subplot 2: correlator output vs adaptive threshold
+ plt.subplot(2, 1, 2)
+ plt.plot(corr_power, label='Correlator Output $|r(t) * p^*(-t)|^2$', color='blue')
+ plt.plot(cfar_thresholds, label='CFAR Adaptive Threshold', color='red', linestyle='--', linewidth=1.5)
+ if len(detections) > 0: # Overlay the detections
+ plt.scatter(detections, corr_power[detections], color='lime', edgecolors='black', label='Detections (Preamble Found)', zorder=5)
+ plt.title("Preamble Correlator Output with Adaptive CFAR Threshold")
+ plt.xlabel("Sample Index")
+ plt.ylabel("Correlation Power")
+ plt.legend()
+ plt.grid(True, alpha=0.3)
+ plt.show()
+
+.. image:: ../_images/detection_cfar.svg
+ :align: center
+ :target: ../_images/detection_cfar.svg
+ :alt: Приклад виходу детектора CFAR
+
+Корелятори преамбули, стійкі до частотного зсуву
+#################################################
+
+Коли центральна частота невідома, виявлення преамбули перетворюється на багатовимірний пошук. В ідеально синхронізованій системі когерентний корелятор працює як узгоджений фільтр і максимізує SNR. Однак частотний зсув створює змінне в часі обертання фази, яке декорелює сигнал із локальним шаблоном і значно знижує чутливість виявлення.
+
+Вплив частотного зсуву :math:`\Delta f` залежить від його величини відносно тривалості преамбули :math:`T_p`.
+
+Невеликий зсув, наприклад через ефект Доплера або дрейф годинника, зазвичай спричинений неточністю локального генератора в ppm або повільним рухом. У цьому випадку :math:`\Delta f \cdot T_p \ll 1`. Кореляційний максимум трохи послаблюється, але часову прив'язку все ще можна відновити.
+
+Якщо частотний зсув зовсім невідомий, наприклад під час холодного захоплення супутника або у високодинамічному каналі БПЛА, когерентна сума може стати нульовою, коли фаза за час преамбули повернеться більш ніж на :math:`180^{\circ}` (:math:`\Delta f > 1/(2T_p)`). Тоді виявлення стає неможливим незалежно від SNR.
+
+Втрату модуля кореляції через частотний зсув описує ядро Діріхле, або періодична sinc-функція. Зі збільшенням зсуву когерентна сума повернутих векторів спадає за sinc-подібним законом.
+
+Втрату в децибелах можна наближено обчислити так:
+
+:math:`L_{dB}(\Delta f) = 20 \log_{10} \left| \frac{\sin(\pi \Delta f N T_{s})}{N \sin(\pi \Delta f T_{s})} \right|`
+
+де:
+
+ - :math:`N` — кількість символів преамбули;
+ - :math:`T_s` — період символу;
+ - :math:`\Delta f` — частотний зсув у герцах.
+
+Зі збільшенням :math:`\Delta f` чисельник коливається, а знаменник зростає, утворюючи нулі чутливості детектора. Для звичайного корелятора перший нуль виникає при :math:`\Delta f = 1/(N T_s)`. Якщо зсув дорівнює половині ширини частотного відліку, втрата становить приблизно 3,9 дБ, що істотно погіршує ефективні SNR і :math:`P_d`.
+
+Методи стійкості до частотного зсуву
+####################################
+
+А. Когерентний сегментований корелятор
+
+Преамбула довжиною :math:`N` ділиться на :math:`M` сегментів довжиною :math:`L=N/M`. Кожен сегмент корелюється когерентно, а результати поєднуються з компенсацією дрейфу фази між сегментами:
+
+:math:`Y_{coh} = \sum_{m=0}^{M-1} \left( \sum_{k=0}^{L-1} r[k+mL] \cdot p^{*}[k] \right) e^{-j \hat{\phi}_m}`
+
+Тут :math:`\hat{\phi}_m` — оцінка повороту фази відповідного сегмента. Метод зберігає виграш SNR повної преамбули, але для суміщення фаз потребує точної оцінки частоти.
+
+Б. Некогерентний сегментований корелятор
+
+Сегменти корелюються когерентно, але додаються квадрати модулів, тож інформація про фазу відкидається:
+
+:math:`Y_{non-coh} = \sum_{m=0}^{M-1} \left| \sum_{k=0}^{L-1} r[k+mL] \cdot p^{*}[k] \right|^{2}`
+
+Цей підхід надзвичайно стійкий до частотних зсувів аж до :math:`1/(L T_s)`. Проте він має втрати некогерентного інтегрування. Додавання модулів замість комплексних значень дає шуму змогу накопичуватися швидше за сигнал, фактично зменшуючи SNR після детектування.
+
+В. Повний перебір частоти
+
+Приймач запускає кілька паралельних кореляторів, кожен із яких зсунутий на окрему частоту :math:`\Delta f_i`.
+
+Цей метод дає найкращу якість за SNR, тобто повний когерентний виграш, але потребує найбільше обчислень. Крок частотних відліків потрібно вибрати за формулою Діріхле досить малим, щоб найгірша втрата між ними залишалася прийнятною, наприклад меншою за 1 дБ.
+
+Під час реалізації в часовій області відліки згортаються з фіксованим набором ваг. Для частотного пошуку потрібний окремий банк FIR для кожного частотного відліку. Для коротких преамбул це ефективно реалізується у FPGA за допомогою блоків Xilinx DSP48. У частотній області для пошуку виконують ШПФ вхідного сигналу й преамбули, адже множення у спектрі відповідає кореляції. Існує корисний прийом частотного зсуву: для перевірки різних зсувів не потрібні окремі ШПФ. Достатньо циклічно переміщувати відліки ШПФ преамбули відносно сигналу перед поелементним множенням і оберненим ШПФ. Для безперервних потоків застосовують обробку блоками Overlap-Save або Overlap-Add, щоб не втрачати кореляційні максимуми на межах вікон ШПФ.
+
+Стійкість до частотного зсуву є компромісом між виграшем обробки та обчислювальною складністю. Некогерентна сегментована кореляція найстійкіша за великої невизначеності, але потребує більшого енергетичного запасу каналу. Когерентні сегментовані методи та повний пошук на основі ШПФ мають вищу чутливість, але споживають значно більше апаратних ресурсів. Розуміння втрат за законом Діріхле критично важливе для вибору щільності частотних відліків у будь-якому приймачі з частотним пошуком.
+
+.. image:: ../_images/detection_freq_offset.svg
+ :align: center
+ :target: ../_images/detection_freq_offset.svg
+ :alt: Вплив частотного зсуву на кореляцію
+
+*****************************************************************
+Виявлення сигналів DSSS
+*****************************************************************
+
+У системі прямого розширення спектра (Direct Sequence Spread Spectrum, DSSS) кореляційний детектор витягує змістовний сигнал із того, що спочатку нагадує випадковий шум. Високошвидкісна чипова послідовність, або код розширення, розподіляє енергію сигналу в значно ширшій смузі, ніж потрібно початковим даним. Загальна потужність залишається сталою, тому її розподіл у ширшому діапазоні знижує спектральну густину потужності (PSD). Таке спектральне розрідження може опустити сигнал нижче теплового шумового порога, зробивши його майже невидимим для звичайних вузькосмугових приймачів. Цільовий приймач натомість застосовує ту саму чипову послідовність для зворотного стискання спектра: енергія концентрується в початковій вузькій смузі, а вузькосмугові завади, навпаки, розширюються. Саме це забезпечує надійне виявлення навіть у дуже шумному середовищі. Далі розглянемо часову сторону цієї задачі.
+
+Роль автокореляційних властивостей
+##################################
+
+Правильний вибір послідовності критично важливий для синхронізації та придушення багатопроменевості. В ідеалі вона повинна мати досконалу автокореляцію: високий максимум за точного суміщення й майже нульові значення за будь-якого іншого часового зсуву. Гострий автокореляційний максимум дає приймачу змогу синхронізуватися із субчиповою точністю. Якщо сигнал відбивається від будівлі й приходить пізніше, хороша автокореляція дає змогу сприйняти затриману копію як некорельований шум, а не руйнівну заваду, послаблюючи багатопроменевість.
+
+Поширені послідовності розширення
+#################################
+
+Різним застосуванням потрібні різні математичні властивості послідовностей. Ось кілька прикладів:
+
+- Коди Баркера мають найкращі можливі автокореляційні властивості для коротких довжин до 13 і відомі застосуванням у Wi-Fi 802.11b.
+- M-послідовності максимальної довжини генеруються регістрами зсуву з лінійним зворотним зв'язком (LFSR) і мають чудову випадковість та автокореляцію протягом дуже довгих періодів.
+- Коди Голда, утворені з пар m-послідовностей, дають велику множину кодів із контрольованою взаємною кореляцією. Тому вони є стандартом GPS і CDMA, де одночасно співіснує багато сигналів.
+- Послідовності Задова—Чу (ZC) — комплексні послідовності зі сталою амплітудою та нульовою автокореляцією для всіх ненульових зсувів; нині вони широко застосовуються для синхронізації LTE і 5G.
+- Коди Касамі подібні до кодів Голда, але мають ще нижчу взаємну кореляцію за тієї самої довжини, що корисно в середовищах із великою щільністю сигналів.
+
+Синхронізація чипів у DSSS
+##########################
+
+Здатність DSSS-приймача відновлювати дані повністю залежить від синхронізації з вхідною чиповою послідовністю. Чипи значно коротші за біти даних, тому навіть мала дробова похибка часу, коли приймач бере відлік між чипами, істотно зменшує кореляційний максимум. Дослідимо вплив дробового часового зсуву, змоделювавши просту DSSS-систему й побудувавши кореляційний вихід для зсувів від 0 до 1 чипа. Тут не виконується повна кореляція: ми лише беремо скалярний добуток за нульової затримки, бо вже знаємо, що саме там буде максимум.
+
+.. code-block:: python
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+
+ # Barker 11 sequence
+ barker11 = np.array([1, -1, 1, 1, -1, 1, 1, 1, -1, -1, -1])
+ samples_per_chip = 100
+
+ # Upsample the sequence to simulate continuous time
+ sig = np.repeat(barker11, samples_per_chip)
+
+ offsets = np.linspace(-1.5, 1.5, 500) # Fractional chip offsets
+ peaks = []
+
+ for offset in offsets:
+ # Shift the signal by a fractional number of chips, converted to samples
+ shift_samples = int(offset * samples_per_chip)
+ if shift_samples > 0:
+ shifted_sig = np.pad(sig, (shift_samples, 0))[:len(sig)]
+ elif shift_samples < 0:
+ shifted_sig = np.pad(sig, (0, abs(shift_samples)))[abs(shift_samples):]
+ else:
+ shifted_sig = sig
+
+ # Compute normalized correlation at zero lag for this offset
+ correlation = np.vdot(sig, shifted_sig) / np.vdot(sig, sig)
+ peaks.append(np.abs(correlation))
+
+ plt.figure(figsize=(10, 5))
+ plt.plot(offsets, peaks, label='Normalized Correlation', color='blue', linewidth=2)
+ plt.axvline(0, color='red', linestyle='--', alpha=0.5, label='Perfect Alignment')
+ plt.title('DSSS Correlation Peak vs. Fractional Chip Timing Offset')
+ plt.xlabel('Offset (Fraction of a Chip)')
+ plt.ylabel('Normalized Correlation Peak Magnitude')
+ plt.grid(True, which='both', linestyle='--', alpha=0.6)
+ plt.legend()
+ plt.savefig('../_images/detection_dsss.svg', bbox_inches='tight')
+ plt.show()
+
+.. image:: ../_images/detection_dsss.svg
+ :align: center
+ :target: ../_images/detection_dsss.svg
+ :alt: DSSS
+
+Як і очікувалося, максимум виникає за нульового зсуву й лінійно спадає, досягаючи половини значення при зсуві на пів чипа. Після зсуву більш ніж на один чип кореляція може начебто знову зростати, але справжній максимум залишається низьким, бо сигнал уже не суміщений із послідовністю.
+
+****************************************************
+Виявлення пакетів у безперервному IQ-потоці
+****************************************************
+
+Досі ми розглядали теоретичні основи виявлення сигналів: корелятори, CFAR і системи з розширеним спектром. Тепер поєднаємо їх для розв'язання поширеної практичної задачі — **виявлення переривчастих пакетів у безперервному потоці IQ-відліків від SDR**. Уявімо, що модем або пристрій IoT передає пакет даних раз на секунду чи з нерегулярними інтервалами. SDR безперервно приймає, наприклад, мільйон відліків за секунду. Пакети надходять у непередбачувані моменти й занурені в шум та завади. Потрібно:
+
+1. Виявити момент надходження пакета.
+2. Визначити точний індекс його першого відліку.
+3. Виділити пакет для подальшої обробки — демодуляції, декодування тощо.
+4. Робити це в реальному часі, не пропускаючи пакетів.
+
+Це принципово відрізняється від обробки заздалегідь записаного IQ-файла, де весь сигнал доступний одразу. Тут відліки надходять безперервно, а рішення потрібно ухвалювати в реальному часі за обмежених обчислювальних ресурсів. Поєднаємо кілька розглянутих методів:
+
+1. **Взаємну кореляцію** для пошуку відомої преамбули.
+2. **CFAR** для адаптивного встановлення порога за змінного шуму.
+3. **Керування буферами** для роботи з безперервним потоком.
+4. **Пошук максимумів** для точного визначення часу пакета.
+
+Для роботи в реальному часі накопичуватимемо відліки в **буферах**, наприклад по 100 000 відліків, запускатимемо детектор для кожного буфера й зберігатимемо стан між ними, щоб не втрачати пакети, які перетинають межу двох буферів.
+
+Реалізація
+##########
+
+Детектор працює за такою схемою:
+
+.. mermaid::
+
+ flowchart TD
+ A("Безперервний IQ-потік від SDR
(частота дискретизації 1 МГц)")
+ B("Накопичення буфера
(100 тис. відліків = 0,1 с)")
+ C("Взаємна кореляція з відомою преамбулою")
+ D("Обчислення порога CFAR")
+ E("Пошук максимумів
(кореляція > поріг)")
+ F("Виділення та перевірка пакета")
+ A --> B --> C --> D --> E --> F
+
+Щоб не пропустити пакет на межі буферів, застосуємо **overlap-save**: кожен буфер міститиме останні ``N_preamble`` відліків попереднього. Тоді пакет, який починається наприкінці буфера ``i``, повністю потрапить до буфера ``i+1``. Обсяг додаткових обчислень невеликий і значно менш шкідливий, ніж пропущені пакети.
+
+Побудуймо повний детектор пакетів у Python крок за кроком. Використаємо коротшу, ніж раніше, преамбулу Задова—Чу й адаптивний CFAR.
+
+Крок 1. Преамбула та параметри
+******************************
+
+.. code-block:: python
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ from scipy.signal import correlate
+
+ # Preamble: Zadoff-Chu sequence (excellent correlation properties)
+ N_zc = 63 # ZC sequence length (typically prime or power of 2 - 1)
+ u = 5 # ZC root
+ t = np.arange(N_zc)
+ preamble = np.exp(-1j * np.pi * u * t * (t + 1) / N_zc)
+
+ # System parameters
+ sample_rate = 1e6
+ buffer_size = 100000
+ overlap_size = len(preamble) # Overlap to catch boundary packets
+
+ # CFAR parameters
+ cfar_guard = 10
+ cfar_train = 50
+ pfa_target = 1e-6
+
+ # Packet parameters (for simulation)
+ packet_length = 500 # Total packet length in samples (preamble + data)
+ snr_db = -5
+
+Крок 2. Функція детектора CFAR
+******************************
+
+Скористаємося CA-CFAR із попереднього прикладу, трохи його оптимізувавши:
+
+.. code-block:: python
+
+ def ca_cfar_1d(signal, num_train, num_guard, pfa):
+ """
+ 1D Cell-Averaging CFAR detector.
+
+ Args:
+ signal: Input signal (typically correlation magnitude)
+ num_train: Number of training cells (on each side)
+ num_guard: Number of guard cells (on each side)
+ pfa: Target probability of false alarm
+
+ Returns:
+ threshold: Adaptive threshold array
+ """
+ n = len(signal)
+ threshold = np.zeros(n)
+ alpha = num_train * (pfa**(-1/num_train) - 1)
+
+ for i in range(n):
+ # Define training window indices
+ train_start_left = max(0, i - num_guard - num_train)
+ train_end_left = max(0, i - num_guard)
+ train_start_right = min(n, i + num_guard + 1)
+ train_end_right = min(n, i + num_guard + num_train + 1)
+
+ # Collect training cells (avoid guard cells and CUT)
+ train_cells = np.concatenate([
+ signal[train_start_left:train_end_left],
+ signal[train_start_right:train_end_right]
+ ])
+
+ if len(train_cells) > 0:
+ noise_est = np.mean(train_cells)
+ threshold[i] = alpha * noise_est
+
+ return threshold
+
+Крок 3. Функція виявлення пакетів
+*********************************
+
+.. code-block:: python
+
+ def detect_packets(buffer, preamble, cfar_guard, cfar_train, pfa,
+ min_spacing=None):
+ """
+ Detect packets in a buffer of IQ samples.
+
+ Args:
+ buffer: Complex IQ samples
+ preamble: Known preamble sequence
+ cfar_guard: CFAR guard cells
+ cfar_train: CFAR training cells
+ pfa: Target false alarm probability
+ min_spacing: Minimum samples between detections (prevents duplicates)
+
+ Returns:
+ detections: List of sample indices where packets start
+ """
+ # Correlate buffer with preamble
+ corr = correlate(buffer, preamble, mode='same')
+ corr_power = np.abs(corr)**2
+
+ # Compute adaptive threshold
+ threshold = ca_cfar_1d(corr_power, cfar_train, cfar_guard, pfa)
+
+ # Find peaks above threshold
+ detections_raw = np.where(corr_power > threshold)[0]
+
+ # Compensate for correlation offset (peak occurs at len(preamble)//2 after true start)
+ half_preamble = len(preamble) // 2
+ detections_raw = detections_raw - half_preamble
+
+ # Remove edge detections (unreliable)
+ half_preamble = len(preamble) // 2
+ detections_raw = detections_raw[
+ (detections_raw > half_preamble) &
+ (detections_raw < len(buffer) - half_preamble)
+ ]
+
+ # Remove duplicate detections (peaks close together)
+ if min_spacing is None:
+ min_spacing = len(preamble)
+
+ detections = []
+ if len(detections_raw) > 0:
+ detections.append(detections_raw[0])
+ for det in detections_raw[1:]:
+ if det - detections[-1] > min_spacing:
+ detections.append(det)
+
+ return detections, corr_power, threshold
+
+Крок 4. Моделювання тестового сигналу
+*************************************
+
+.. code-block:: python
+
+ def generate_packet_stream(preamble, packet_length, num_packets,
+ sample_rate, snr_db):
+ """
+ Generate a simulated IQ stream with intermittent packets.
+
+ Returns:
+ signal: Complex IQ samples
+ true_starts: Ground truth packet start indices
+ """
+ # Calculate noise power from SNR
+ signal_power = 1.0 # Normalized preamble power
+ noise_power = signal_power / (10**(snr_db/10))
+ noise_std = np.sqrt(noise_power / 2) # Complex noise
+
+ # Generate QPSK data (random payload after preamble)
+ qpsk_map = np.array([1+1j, -1+1j, -1-1j, 1-1j]) / np.sqrt(2)
+
+ # Time between packets (1 second +/- 20% jitter)
+ packets_per_sec = 1
+ avg_gap = int(sample_rate / packets_per_sec)
+
+ signal = []
+ true_starts = []
+
+ for i in range(num_packets):
+ # Add gap (noise only)
+ if i == 0:
+ gap_length = np.random.randint(avg_gap//2, avg_gap)
+ else:
+ gap_length = np.random.randint(int(avg_gap*0.8), int(avg_gap*1.2))
+
+ noise = noise_std * (np.random.randn(gap_length) +
+ 1j*np.random.randn(gap_length))
+ signal.extend(noise)
+
+ # Record true packet start
+ true_starts.append(len(signal))
+
+ # Add packet (preamble + data)
+ data_length = packet_length - len(preamble)
+ data = qpsk_map[np.random.randint(0, 4, data_length)]
+ packet = np.concatenate([preamble, data])
+
+ # Add noise to packet
+ packet_noisy = packet + noise_std * (np.random.randn(len(packet)) +
+ 1j*np.random.randn(len(packet)))
+ signal.extend(packet_noisy)
+
+ # Add final gap
+ gap_length = np.random.randint(avg_gap//2, avg_gap)
+ noise = noise_std * (np.random.randn(gap_length) +
+ 1j*np.random.randn(gap_length))
+ signal.extend(noise)
+
+ return np.array(signal), true_starts
+
+ # Generate 5 seconds of signal with ~5 packets
+ signal, true_starts = generate_packet_stream(
+ preamble, packet_length, num_packets=5,
+ sample_rate=sample_rate, snr_db=snr_db
+ )
+
+ print(f"Generated {len(signal)} samples ({len(signal)/sample_rate:.1f} sec)")
+ print(f"True packet starts: {true_starts}")
+
+Крок 5. Виявлення в потоковому режимі
+*************************************
+
+Тепер обробимо сигнал блоками, імітуючи потік у реальному часі:
+
+.. code-block:: python
+
+ def process_stream(signal, preamble, buffer_size, overlap_size,
+ cfar_guard, cfar_train, pfa):
+ """
+ Process continuous IQ stream in buffers (simulates real-time).
+
+ Returns:
+ all_detections: List of detected packet starts (global indices)
+ """
+ all_detections = []
+ n_samples = len(signal)
+ current_pos = 0
+
+ while current_pos < n_samples:
+ # Define buffer with overlap
+ buffer_start = max(0, current_pos - overlap_size)
+ buffer_end = min(n_samples, current_pos + buffer_size)
+ buffer = signal[buffer_start:buffer_end]
+
+ # Detect packets in this buffer
+ detections, corr_power, threshold = detect_packets(
+ buffer, preamble, cfar_guard, cfar_train, pfa
+ )
+
+ # Convert buffer-relative indices to global indices
+ for det in detections:
+ global_idx = buffer_start + det
+
+ # Avoid duplicate detections from overlap region
+ if len(all_detections) == 0 or \
+ global_idx - all_detections[-1] > len(preamble):
+ all_detections.append(global_idx)
+
+ current_pos += buffer_size
+
+ return all_detections
+
+
+ detected_starts = process_stream(
+ signal, preamble, buffer_size, overlap_size,
+ cfar_guard, cfar_train, pfa_target
+ )
+
+ print(f"\nDetection Results:")
+ print(f"True packets: {len(true_starts)}")
+ print(f"Detected packets: {len(detected_starts)}")
+ print(f"Detected starts: {detected_starts}")
+
+Крок 6. Оцінювання характеристик
+********************************
+
+.. code-block:: python
+
+ # Calculate detection statistics
+ tolerance = len(preamble)
+
+ matched_detections = []
+ false_alarms = []
+
+ for det in detected_starts:
+ # Check if detection matches any true packet
+ matched = False
+ for true_start in true_starts:
+ if abs(det - true_start) <= tolerance:
+ matched_detections.append(det)
+ matched = True
+ break
+ if not matched:
+ false_alarms.append(det)
+
+ missed_packets = len(true_starts) - len(matched_detections)
+
+ print(f"\nPerformance Metrics:")
+ print(f" Correct detections: {len(matched_detections)}/{len(true_starts)}")
+ print(f" Missed packets: {missed_packets}")
+ print(f" False alarms: {len(false_alarms)}")
+
+ # Calculate timing errors
+ timing_errors = []
+ for det in matched_detections:
+ errors = [abs(det - ts) for ts in true_starts]
+ timing_errors.append(min(errors))
+
+ if len(timing_errors) > 0:
+ print(f" Timing error (avg): {np.mean(timing_errors):.1f} samples")
+ print(f" Timing error (max): {np.max(timing_errors):.1f} samples")
+
+Крок 7. Візуалізація результатів
+********************************
+
+.. code-block:: python
+
+ # Process one buffer for detailed visualization
+ buffer_start = max(0, true_starts[0] - 5000)
+ buffer_end = min(len(signal), true_starts[0] + 20000)
+ viz_buffer = signal[buffer_start:buffer_end]
+
+ detections_viz, corr_viz, thresh_viz = detect_packets(
+ viz_buffer, preamble, cfar_guard, cfar_train, pfa_target
+ )
+
+ # Convert to global indices for plotting
+ detections_viz_global = [d + buffer_start for d in detections_viz]
+
+ # Create visualization
+ fig, axes = plt.subplots(3, 1, figsize=(14, 10))
+ time_axis = (np.arange(len(viz_buffer)) + buffer_start) / sample_rate * 1000 # ms
+
+ # Subplot 1: Received signal power
+ axes[0].plot(time_axis, np.abs(viz_buffer)**2, 'gray', alpha=0.6, linewidth=0.5)
+ axes[0].set_ylabel('Power')
+ axes[0].set_title('Received IQ Signal Power')
+ axes[0].grid(True, alpha=0.3)
+
+ # Mark true packet locations
+ for ts in true_starts:
+ if buffer_start <= ts <= buffer_end:
+ t_ms = ts / sample_rate * 1000
+ axes[0].axvline(t_ms, color='green', linestyle='--', alpha=0.7,
+ label='True Packet' if ts == true_starts[0] else '')
+ axes[0].legend()
+
+ # Subplot 2: Correlation output
+ axes[1].plot(time_axis, corr_viz, 'blue', linewidth=1, label='Correlation')
+ axes[1].plot(time_axis, thresh_viz, 'red', linestyle='--', linewidth=1.5,
+ label='CFAR Threshold')
+ axes[1].set_ylabel('Correlation Power')
+ axes[1].set_title('Preamble Correlation with Adaptive CFAR Threshold')
+ axes[1].grid(True, alpha=0.3)
+ axes[1].legend()
+
+ # Subplot 3: Detections
+ detection_mask = np.zeros(len(viz_buffer))
+ for det in detections_viz:
+ detection_mask[det] = corr_viz[det]
+
+ axes[2].plot(time_axis, corr_viz, 'blue', alpha=0.4, linewidth=0.8)
+ axes[2].scatter(time_axis[detection_mask > 0], detection_mask[detection_mask > 0],
+ color='lime', edgecolors='black', s=100, zorder=5,
+ label='Detected Packets')
+ axes[2].set_xlabel('Time (ms)')
+ axes[2].set_ylabel('Correlation Power')
+ axes[2].set_title('Detected Packet Locations')
+ axes[2].grid(True, alpha=0.3)
+ axes[2].legend()
+
+ plt.tight_layout()
+ plt.show()
+
+Візуалізація має показати:
+
+1. **Верхній графік** — необроблену потужність сигналу з позначеними справжніми положеннями пакетів.
+2. **Середній графік** — вихід корелятора й адаптивний поріг CFAR, який слідує за шумовим порогом.
+3. **Нижній графік** — виявлені пакети, виділені як максимуми над порогом.
+
+.. image:: ../_images/detection_realtime.png
+ :align: center
+ :scale: 50 %
+ :alt: Результати виявлення пакетів у реальному часі
+
+Практичні міркування та налаштування
+####################################
+
+Компроміси розміру буфера
+*************************
+
+**Більші буфери**, наприклад 1 млн відліків:
+
+- ✅ Краще оцінювання шуму CFAR завдяки більшій кількості навчальних комірок.
+- ✅ Менші накладні обчислювальні витрати через рідші виклики обробки.
+- ❌ Вища затримка, бо потрібно чекати заповнення буфера.
+- ❌ Більше споживання пам'яті.
+
+**Менші буфери**, наприклад 10 тис. відліків:
+
+- ✅ Нижча затримка й швидша реакція.
+- ✅ Менше споживання пам'яті.
+- ❌ Гірша робота CFAR через меншу кількість навчальних комірок.
+- ❌ Вище завантаження процесора через частішу обробку.
+
+**Рекомендація:** почніть із буфера, довжина якого у 10–100 разів більша за преамбулу. Для преамбули з 63 відліків і частоти 1 Мвідл/с спробуйте від 10 до 100 тис. відліків.
+
+Налаштування параметрів CFAR
+****************************
+
+Поведінку детектора визначають три параметри:
+
+**num_guard** — захисні комірки:
+
+- Не дають енергії сигналу потрапити в оцінку шуму.
+- Надто мале значення: сигнал потрапляє в навчальну область, підвищує поріг і спричиняє пропуски.
+- Надто велике: залишається менше навчальних комірок, тому оцінка шуму погіршується.
+- Практичне правило: приблизно 0,5–1,0 довжини преамбули.
+
+**num_train** — навчальні комірки:
+
+- Оцінюють локальний шумовий поріг.
+- Надто мале значення: шумний поріг і більше хибних тривог або пропусків.
+- Надто велике: поріг недостатньо швидко реагує на зміни шуму.
+- Практичне правило: приблизно 3–5 довжин преамбули.
+
+**pfa** — імовірність хибної тривоги:
+
+- Керує чутливістю виявлення.
+- Надто високе значення, наприклад 1e-2: багато хибних тривог.
+- Надто низьке, наприклад 1e-10: пропуски слабких пакетів.
+- Практичне правило: почніть із 1e-5 для однієї затримки, а потім налаштуйте за глобальною частотою хибних тривог системи.
+
+Пам'ятайте наведений раніше зв'язок між імовірністю хибної тривоги для однієї затримки й усієї системи.
+
+****************
+Rake-приймачі
+****************
+
+Уявіть, що ви стоїте надворі з телефоном, а сигнал базової станції доходить до вас кількома шляхами. Частина енергії надходить прямо, інша копія відбивається від будівлі за 150 м і приходить трохи пізніше, а третя, можливо, відбивається від схилу й затримується ще більше. Кожен шлях має різну довжину, тому копії потрапляють до приймача в різні моменти, з різною потужністю та фазою. Це багатопроменеве поширення — звичайний стан будь-якого реального бездротового середовища, а не рідкісний виняток.
+
+Як багатопроменевість виглядає на виході корелятора? Коли преамбула або чипова послідовність ковзає вздовж прийнятих відліків, виникає не один чистий максимум, а кілька: один для прямого шляху й менші для кожного відлуння. Найпростіше було б вибрати найвищий максимум, а решту відкинути. Проте менші максимуми — не шум, а ті самі дані, що прийшли іншим маршрутом. Відкидаючи їх, ми втрачаємо корисну енергію сигналу.
+
+Rake-приймач відповідає на просте запитання: навіщо відкидати відлуння, якщо їх можна зібрати й додати? Назва походить від садових грабель, кожен зубець яких збирає окрему смугу. Тут кожен «палець» є окремим корелятором, налаштованим на одну з багатопроменевих затримок. Перший палець синхронізується з прямим шляхом, другий — із відбиттям від будівлі, третій — із відбиттям від схилу тощо. Кожен виконує ту саму взаємну кореляцію, що використовувалася протягом розділу, але його шаблон зсунутий на іншу затримку й відстежує конкретну копію сигналу.
+
+Пошук пальців
+#############
+
+Перш ніж rake-приймач зможе щось об'єднати, він має визначити положення копій. Для цього використовується вже знайомий вихід корелятора. Набір багатопроменевих затримок і їхніх відносних потужностей називають *профілем потужності затримок* каналу. Це просто графік модулів кореляційних максимумів залежно від затримки. Для налаштування rake-приймача профіль переглядають, вибирають кілька найсильніших максимумів і призначають кожному окремий палець. Приймач може мати три або чотири пальці, адже виділяти палець для крихітного максимуму, зануреного в шум, коштує більше, ніж дає користі.
+
+Канал змінюється, коли ви рухаєтеся, повз проїжджають автомобілі або змінюється середовище, тому пальці не можна налаштувати один раз і забути. Приймач постійно повторно сканує профіль затримок і перепризначає пальці, коли максимуми зростають, завмирають або зміщуються в часі. Якщо сильне відбиття, яке відстежував палець, згасає, приймач перепризначає його новому шляху, що став потужнішим.
+
+Об'єднання пальців
+##################
+
+Після синхронізації кожного пальця зі своєю копією приймач має кілька незалежних вимірювань тих самих переданих символів. Як перетворити їх на одне рішення? Просте додавання було б помилкою, адже сильна чиста копія та слабка зашумлена не заслуговують однакової ваги. Слабкий палець переважно містить шум, тому рівне зважування погіршило б спільний результат.
+
+Стандартне рішення — *об'єднання з максимальним відношенням* (Maximal Ratio Combining, MRC), яке перед сумуванням зважує кожен палець відповідно до його потужності. Максимум із високим SNR отримує велику вагу, а палець ледве над шумом — малу. Позначимо комплексний вихід пальця :math:`k` для заданого символу як :math:`r_k`, а коефіцієнт каналу цього пальця як :math:`h_k`. Тоді об'єднаний вихід дорівнює
+
+.. math::
+
+ y = \sum_{k} h_k^{*} \, r_k
+
+Інакше кажучи, вихід кожного пальця множиться на комплексно спряжений коефіцієнт його каналу, після чого всі результати додаються. Спряження одночасно виконує дві дії: модуль :math:`|h_k|` масштабує пальці за їхньою потужністю, тому сильні мають більший вплив, а фаза повертає внесок кожного пальця так, щоб усі копії суміщалися й додавалися конструктивно, а не частково компенсували одна одну.
+
+Виграш має дві складові. Перша — просто більша енергія сигналу: збираючи відлуння замість їх відкидання, rake-приймач повертає потужність, яку детектор одного максимуму втратив би. Друга, часто важливіша перевага — *рознесення*. Глибокі завмирання, що руйнують бездротовий канал, виникають, коли шляхи в певній точці тимчасово взаємно компенсуються. Малоймовірно, що кілька шляхів із різними затримками завмруть одночасно. Тому коли прямий шлях потрапляє в завмирання, одне з відбиттів, імовірно, залишається сильним, а rake-приймач спирається на ті пальці, які в цей момент мають кращу якість.
+
+Зауважте, що rake-приймачі можливі лише для широкосмугових сигналів із гострою автокореляцією, коли багатопроменеві копії утворюють чітко розділені максимуми в профілі затримок. Вузькосмугові сигнали розмивають копії так, що їх неможливо розрізнити, тому rake-приймач для них не працює.
diff --git a/content-ukraine/fpv_video.rst b/content-ukraine/fpv_video.rst
new file mode 100644
index 00000000..fc8f3e35
--- /dev/null
+++ b/content-ukraine/fpv_video.rst
@@ -0,0 +1,152 @@
+.. _fpv-chapter:
+
+############################
+Аналогові відеосигнали FPV
+############################
+
+У цьому розділі ми розглянемо аналогові відеосигнали, які використовують у більшості аматорських і саморобних FPV-дронів. Вони являють собою ЧМ-модульовані сигнали NTSC або PAL. Ми проаналізуємо ці сигнали та покажемо, як їх демодулювати, щоб відновити відеозображення.
+
+****************
+Вступ
+****************
+
+Аналоговий FPV (First-Person View, «вид від першої особи») — це традиційний спосіб передавання відео в реальному часі з радіокерованої моделі, найчастіше дрона або квадрокоптера, до пілота. Замість цифрового кодування відео аналогова система передає сигнал з камери у форматі NTSC або PAL, модулює ним частоту несної та випромінює його, зазвичай у діапазоні 5,8 ГГц (хоча також використовують діапазони 1,2–1,3 ГГц і 2,4 ГГц). Це не цифровий сигнал: тут немає ані стиснення, ані шифрування. Визначальна перевага аналогового FPV — надзвичайно мала затримка. Оскільки відео не стискається й майже не обробляється, пілот бачить зображення з камери практично миттєво, що критично важливо для швидкого й точного керування. Такі системи зазвичай недорогі: аналоговий FPV-відеопередавач можна придбати приблизно за 10 доларів, а моноблок із камерою та антеною — приблизно за 20 доларів. Більшість аналогових FPV-відеопередавачів також можуть передавати звуковий сигнал, який додається до відеосигналу. Для керування моделлю аналогове відео зазвичай поєднують з окремою радіосистемою, наприклад FrSky, FlySky, Spektrum або ELRS.
+
+****************
+Будова сигналу
+****************
+
+З погляду приймача після частотної демодуляції сигнал складається з таких компонентів:
+
+.. image:: ../_images/fpv_baseband_spectrum_after_demod.svg
+ :align: center
+ :target: ../_images/fpv_baseband_spectrum_after_demod.svg
+ :alt: Спектр сигналу в основній смузі після частотної демодуляції
+
+Одна з приємних переваг ЧМ полягає в тому, що приймач не обов'язково має бути ідеально налаштований на центр сигналу. Поки весь сигнал перебуває в смузі огляду приймача, частотний демодулятор працюватиме належним чином. Причина в тому, що частотна демодуляція спирається на зміни частоти, а не на її абсолютне значення. Отже, достатньо, щоб сигнал мав належну потужність і вміщувався в смугу пропускання приймача. Проте бажано все ж розташувати його приблизно по центру, щоб перед частотною демодуляцією відфільтрувати зайвий шум.
+
+Розгляньмо приклад. IQ-запис сигналу NTSC, використаний у коді цього розділу, можна завантажити `тут `_. Зауважте, що запис містить лише кілька кадрів.
+
+На графіку спектральної густини потужності необробленого радіочастотного сигналу видно ЧМ-сигнал із центром на 0 Гц. Це відповідає частоті 5,925 ГГц, на яку було налаштовано SDR. Вона є центральною частотою одного зі стандартних каналів FPV. Ширина смуги сигналу становить приблизно 6 МГц.
+
+.. image:: ../_images/fpv_psd_raw_rf.svg
+ :align: center
+ :target: ../_images/fpv_psd_raw_rf.svg
+ :alt: Спектральна густина потужності необробленого радіочастотного сигналу
+
+Після частотної демодуляції, яку можна виконати одним рядком Python — :code:`np.angle(x[1:] * np.conj(x[:-1]))`, — отримаємо таку спектральну густину потужності:
+
+.. image:: ../_images/fpv_psd_after_fm_demod.svg
+ :align: center
+ :target: ../_images/fpv_psd_after_fm_demod.svg
+ :alt: Спектральна густина потужності сигналу після частотної демодуляції
+
+У цьому прикладі звукового сигналу немає. Натомість добре видно складову кольоровості. Якщо збільшити низькочастотну ділянку, побачимо гармоніки на частотах, кратних 15,734 кГц (для PAL — 15,625 кГц). Вони відповідають сигналу горизонтальної синхронізації, який з'являється один раз у кожному рядку відео.
+
+.. image:: ../_images/fpv_psd_after_fm_demod_harmomics.svg
+ :align: center
+ :target: ../_images/fpv_psd_after_fm_demod_harmomics.svg
+ :alt: Низькочастотні гармоніки в спектрі сигналу після частотної демодуляції
+
+Щоб краще зрозуміти сигнал, погляньмо на нього в часовій області. На наступному рисунку показано один рядок відеосигналу, знову ж таки після частотної демодуляції. На початку та наприкінці бачимо імпульс горизонтальної синхронізації, а невеликі коливання одразу після нього — це пакет колірної синхронізації (color burst). Він слугує опорним сигналом, за яким приймач декодує інформацію про колір. Решта сигналу містить відеодані — як яскравість чорно-білого зображення, так і інформацію про колір.
+
+.. image:: ../_images/fpv_time_domain_one_line.svg
+ :align: center
+ :target: ../_images/fpv_time_domain_one_line.svg
+ :alt: Один рядок сигналу в часовій області
+
+Якщо збільшити часовий масштаб, можна побачити особливу послідовність синхронізації, яка з'являється один раз на кадр і називається імпульсом вертикальної синхронізації. Це незмінна послідовність імпульсів, яка повідомляє приймачу про початок нового кадру. На рисунку нижче її показано в першій половині графіка. Далі наведено ще 13 рядків сигналу, подібних до розглянутого вище, але в меншому масштабі.
+
+.. image:: ../_images/fpv_time_domain.svg
+ :align: center
+ :target: ../_images/fpv_time_domain.svg
+ :alt: Сигнал у часовій області із синхронізацією кадру та кількома рядками
+
+********************************
+Демодуляція відео
+********************************
+
+Щоб демодулювати відеосигнал і відновити зображення, виконаємо такі кроки:
+
+#. Відфільтруємо звуковий сигнал.
+#. Передискретизуємо сигнали яскравості та кольоровості рівно до 508 відліків на рядок, щоб кожен відлік відповідав одному пікселю.
+#. Перетворимо одновимірний масив відліків на двовимірне зображення.
+#. Масштабуємо значення до діапазону 0–255 і покажемо результат як зображення у відтінках сірого.
+
+Зауважте, що цей процес відновлює лише чорно-білу складову. Інформація про колір закодована інакше, тому відновити її складніше.
+
+Нижче наведено повний робочий приклад для запису, який можна завантажити `тут `_.
+
+.. code-block:: python
+
+ import numpy as np
+ import matplotlib.pyplot as plt
+ import scipy.signal as signal
+
+ filename = 'ntsc_remy_10MHz_5925Hz_500ksamples_cf32.iq'
+ x = np.fromfile(filename, dtype=np.complex64)
+ sample_rate = 10e6
+ color_subcarrier_freq = 3.579545e6 # NTSC. higher than luma carrier, not relative to center freq
+ # color_subcarrier_freq = 4.43361875e6 # PAL and SECAM
+ relative_audio_subcarrier_freq = 3.5e6 # the audio might show up at 5.5, 6.0, or 6.5 MHz
+
+ # NTSC constants
+ samples_per_line = 508
+ lines_per_frame = 525
+ refresh_Hz = 30.0/1.001 # almost exactly 29.97 # not exactly 30 Hz!! makes difference
+
+ # PAL constants
+ #samples_per_line = 512
+ #lines_per_frame = 625 # (576 visible lines)
+ #refresh_Hz = 25
+
+ samples_per_frame = samples_per_line * lines_per_frame // 2 # NTSC's vertical sync repeats every field (half-frame), not every full frame
+ line_Hz = refresh_Hz * lines_per_frame
+
+ # FM demodulation
+ x_demod = np.angle(x[1:] * np.conj(x[:-1]))
+
+ # Filter out audio from demodded signal
+ h = signal.firwin(301, 3e6, fs=sample_rate) # for the 10 Mhz recording
+ x_demod = np.convolve(x_demod, h, 'same')
+
+ # Resample luma and chroma to exactly L samples per line
+ resampling_rate = samples_per_line / (sample_rate / line_Hz)
+ resampling_rate *= 1.00003 # fixes the drift, not 100% sure where it comes from, perhaps sample clock offset
+ x_demod = signal.resample(x_demod, int(len(x_demod)*resampling_rate))
+ print("Resampling rate:", resampling_rate)
+
+ # crop to 1 frames worth of samples (essentially a manual sync)
+ if False:
+ manually_tuned_offset = 122250 # for both frame sync and horizontal sync
+ x_demod = x_demod[manually_tuned_offset:manually_tuned_offset+samples_per_frame]
+
+ # reshape into 2D
+ x_demod = x_demod[:len(x_demod) - (len(x_demod) % samples_per_line)] # trim to multiple of samples_per_line
+ frame = x_demod.reshape(-1, samples_per_line) # type: ignore
+
+ # Normalize to 0-255 and convert to uint8
+ frame_norm = frame - np.min(frame)
+ frame_norm = frame_norm / np.max(frame_norm)
+ frame_uint8 = (frame_norm * 255).astype(np.uint8)
+
+ # Display as single image with fixed scaling
+ plt.imshow(frame_uint8, cmap='gray', aspect='auto', vmin=0, vmax=255)
+ plt.axis('off')
+ plt.show()
+
+Якщо запустити цей код без змін, обробка почнеться від початку запису, тобто з випадкового моменту часу. Без синхронізації за імпульсом горизонтальної розгортки кожен рядок просто зсувається на однакову величину. Тому зображення все одно залишається зрозумілим, але зміщується по горизонталі та вертикалі. Крім того, ми одночасно бачимо відліки з кількох кадрів.
+
+.. image:: ../_images/fpv_image_no_sync.svg
+ :align: center
+ :target: ../_images/fpv_image_no_sync.svg
+ :alt: Демодульоване зображення без виділення одного кадру
+
+Обмежити дані одним кадром нескладно, адже значення :code:`samples_per_frame` ми вже обчислили. Проте спочатку потрібно синхронізуватися з початком кадру. Це можна зробити різними способами. Наприклад, можна знайти послідовність вертикальної синхронізації за допомогою кореляції — відтворивши її самостійно або використавши її запис із високим відношенням сигнал/шум. Також можна побудувати сигнал у часовій області й знайти послідовність візуально. Нижче показано зображення після синхронізації. У наведеному коді її виконано вручну: відомо, що новий кадр починається з відліку 122250.
+
+.. image:: ../_images/fpv_image_one_frame.svg
+ :align: center
+ :target: ../_images/fpv_image_one_frame.svg
+ :alt: Демодульоване зображення після виділення одного кадру
+
+Якщо хтось захоче запропонувати надійний демодулятор кольору мовою Python, будь ласка, зв'яжіться з автором. Водночас потрібно показати, що він працює з різними записами — наприклад, як із синтетичними, так і з реальними сигналами — без ручного налаштування параметрів.
diff --git a/content-ukraine/tdoa.rst b/content-ukraine/tdoa.rst
new file mode 100644
index 00000000..66f3a53a
--- /dev/null
+++ b/content-ukraine/tdoa.rst
@@ -0,0 +1,752 @@
+.. _tdoa-chapter:
+
+####
+TDOA
+####
+
+Різниця часу прибуття (Time Difference of Arrival, TDOA) — це метод визначення положення передавача (також званого випромінювачем) за допомогою кількох синхронізованих приймачів (також званих сенсорами) шляхом порівняння часу прибуття сигналу. У цьому розділі розглянуто весь ланцюжок TDOA: геометрію, оцінювання часової затримки методом GCC-PHAT, локалізацію в замкненій формі та методом максимальної правдоподібності, межі точності (CRLB і GDOP), а також такі проблеми, як синхронізація й багатопроменеве поширення. TDOA широко застосовують як у радіочастотних, так і в акустичних та гідроакустичних системах.
+
+Перш ніж заглиблюватися в теорію, спробуйте інтерактивну демонстрацію нижче. Вона допоможе швидко відчути принцип роботи TDOA, заснований на перетині гіпербол.
+
+.. raw:: html
+
+
+
+
+************
+Вступ
+************
+
+У радіотехніці, акустиці та гідроакустиці часто виникає потреба визначити положення випромінювача — виконати його геолокацію. Випромінювач може бути кооперативним, як мобільний телефон, який потрібно знайти, або некооперативним, як радіолокатор, оператор якого не бажає розкривати своє положення. Він може бути нерухомим чи рухомим, а середовищем поширення можуть бути повітря, вода або вільний простір. Локалізацію за TDOA використовують для визначення положення абонента під час екстреного виклику, в акустичних мікрофонних решітках (наприклад, у системах виявлення пострілів на міських ліхтарях), пасивній гідроакустиці, пасивній радіолокації, радіоелектронній боротьбі, радіотехнічній розвідці та навіть для стеження за дикими тваринами. Інженерні деталі в кожному випадку різні, але математична основа однакова.
+
+Головна ідея TDOA полягає в тому, що різниця часу, коли один і той самий хвильовий фронт досягає двох сенсорів, залежить лише від геометрії, а не від моменту випромінювання. Розгляньмо час поширення від випромінювача до сенсора :math:`i`: :math:`t_i = t_0 + r_i / c`, де :math:`t_0` — невідомий момент початку передавання, :math:`r_i` — відстань від випромінювача до сенсора, а :math:`c` — швидкість поширення. Віднімемо час прибуття на два сенсори:
+
+.. math::
+
+ \tau_{ij} = t_i - t_j = \frac{r_i - r_j}{c},
+
+Невідоме :math:`t_0` зникає, і це добре, адже ми, найімовірніше, ніколи його не дізнаємося. TDOA залежить лише від *різниці* дальностей, яка визначається тільки взаємним розташуванням випромінювача й сенсорів. Саме тому TDOA особливо корисний для некооперативних випромінювачів: нам не потрібно знати момент передавання, достатньо виміряти відносні затримки, з якими той самий хвильовий фронт досяг синхронізованих приймачів. Водночас сигнал потрібно виділити так, щоб спостерігався лише один випромінювач, тому можуть знадобитися виявлення, класифікація та фільтрація сигналу.
+
+Кожна пара сенсорів дає одне значення TDOA, а кожному значенню TDOA відповідає одна гіпербола. Отже, кількість гіпербол дорівнює кількості пар сенсорів. Для :math:`N` сенсорів вона становить
+
+.. math::
+
+ \binom{N}{2} = \frac{N(N-1)}{2},
+
+тобто три сенсори дають три гіперболи, чотири — шість, п'ять — десять тощо. Не всі ці вимірювання незалежні: як побачимо далі, нову геометричну інформацію несуть лише :math:`N-1` із них. Однак повний набір усе одно корисний для усереднення шуму.
+
+Ціна такої можливості — вимога до приймачів мати спільну надзвичайно точну опору часу. Як буде показано далі, це саме по собі складне інженерне завдання: принаймні в радіочастотних системах похибка часу в одну наносекунду відповідає похибці дальності приблизно 0,3 метра.
+
+****************
+Геометрія TDOA
+****************
+
+Від різниці часу до різниці дальностей
+======================================
+
+Множення виміряної різниці часу на швидкість поширення перетворює її на *різницю дальностей*:
+
+.. math::
+
+ \Delta r_{ij} = c\,\tau_{ij} = r_i - r_j .
+
+Для акустичних задач у повітрі :math:`c \approx 343` м/с, а для радіосигналів :math:`c \approx 2.998\times10^8` м/с. Звідси одразу видно наслідок для точності: у повітрі похибка часу :math:`0.1` мс відповідає лише приблизно :math:`\sim`\3 см, тоді як у вільному просторі та сама похибка дає 30 км. Тому радіочастотний TDOA потребує винятково точного вимірювання часу — до цієї теми ми ще неодноразово повернемося.
+
+На схемі нижче показано випромінювач і три сенсори, а також часові графіки сигналу, який досягає кожного сенсора в різні моменти.
+
+.. image:: ../_images/tdoa_principle.svg
+ :align: center
+ :target: ../_images/tdoa_principle.svg
+ :alt: Випромінювач і три сенсори з часовими графіками сигналів, прийнятих у різні моменти
+
+Гіпербола
+=========
+
+Що насправді повідомляє одна різниця дальностей? Припустімо, є два сенсори, а вимірювання показало, що випромінювач на 100 метрів ближчий до одного з них. Де він може перебувати? Не в одній конкретній точці, а будь-де вздовж кривої лінії. Під час руху вздовж цієї лінії обидві відстані до сенсорів змінюються, але їхня *різниця* весь час залишається рівною 100 метрам.
+
+Ця крива називається **гіперболою**, а два сенсори розташовані в її фокусах. У тривимірному просторі та сама ідея утворює криву поверхню — гіперболоїд, проте двовимірний випадок простіше уявити, а всі міркування переносяться й у 3D. Позначимо положення випромінювача як :math:`\mathbf{u}`, а положення сенсорів як :math:`\mathbf{s}_i` і :math:`\mathbf{s}_j`. Тоді гіпербола є множиною точок, які задовольняють рівняння
+
+.. math::
+
+ |\mathbf{u}-\mathbf{s}_i| - |\mathbf{u}-\mathbf{s}_j| = \Delta r_{ij} = \text{constant},
+
+яке читається так: відстань до одного сенсора мінус відстань до другого дорівнює виміряній різниці дальностей. Звідси безпосередньо випливають кілька практичних наслідків:
+
+* **Різниця дальностей не може перевищувати відстань між двома сенсорами.** Інтуїтивно різниця двох відстаней найбільша, коли випромінювач лежить за одним із сенсорів на прямій, що їх сполучає. Навіть тоді вона може лише дорівнювати відстані між сенсорами, яку часто називають *базою*. Якщо виміряна різниця більша за базу, щось негаразд: найімовірнішими причинами є шум, багатопроменеве поширення або похибка часу чи синхронізації.
+* **Знак указує, з якого боку розташований випромінювач.** Гіпербола має дві дзеркальні гілки, кожна з яких вигинається до одного із сенсорів. Додатний або від'ємний знак різниці дальностей вибирає гілку біля ближчого сенсора, усуваючи неоднозначність між половинами.
+* **Форма залежить від виміряного значення.** Коли різниця дальностей близька до повної довжини бази, гіпербола притискається до прямої між сенсорами. Коли різниця близька до нуля, тобто випромінювач майже рівновіддалений від них, крива випрямляється в серединний перпендикуляр до бази. Поблизу обох крайніх випадків геометрія стає погано обумовленою: малі похибки вимірювань сильно зміщують оцінку положення, тому локалізація менш надійна.
+
+Отже, одне значення TDOA обмежує положення випромінювача кривою, а не точкою. Щоб визначити точку, потрібно перетнути кілька таких кривих. Нижче показано два сенсори та кілька гілок гіпербол для :math:`\Delta r < 0`, :math:`\Delta r = 0` (серединний перпендикуляр) і :math:`\Delta r > 0`. На кожній гіперболі TDOA між двома сенсорами постійна. Маючи лише два сенсори, ми знали б, що випромінювач лежить десь на відповідній лінії; третій сенсор потрібний, щоб знайти конкретну точку на ній, тобто виконати геолокацію.
+
+.. image:: ../_images/tdoa_hyperbola.svg
+ :align: center
+ :target: ../_images/tdoa_hyperbola.svg
+ :alt: Два сенсори та кілька гілок гіпербол
+
+Мультилатерація
+===============
+
+Маючи :math:`N` сенсорів, можна утворити їхні пари й перетнути відповідні гіперболи; випромінювач лежить у точці або поблизу точки їхнього спільного перетину. Цей процес називають **гіперболічною мультилатерацією**. Підрахунок ступенів вільності показує потрібну кількість сенсорів:
+
+* У **2D** положення випромінювача має дві невідомі координати :math:`(x,y)`. Кожне незалежне TDOA дає одне рівняння, тому потрібні щонайменше два незалежні TDOA, а отже, **три сенсори**. Наприклад, такий підхід придатний, якщо відомо, що випромінювач розташований на поверхні землі й кривизною Землі можна знехтувати.
+* У **3D** положення має три невідомі координати :math:`(x,y,z)`, тому потрібні три незалежні TDOA і, відповідно, **чотири сенсори**.
+
+За відсутності шуму гіперболи перетинаються в одній точці. Іноді виникає геометрична неоднозначність, яку усувають знаки гілок або додатковий сенсор. Якщо сенсорів більше за мінімальну кількість, система стає *перевизначеною*: через шум гіперболи вже не мають точної спільної точки, тож потрібно розв'язувати задачу найменших квадратів або максимальної правдоподібності, як описано далі.
+
+Опорний сенсор і незалежні пари
+===============================
+
+Із :math:`N` сенсорів можна утворити :math:`\binom{N}{2}` попарних TDOA, але не всі вони незалежні. Виберімо один сенсор опорним, наприклад сенсор 1, і сформуймо :math:`\tau_{i1}` для :math:`i = 2,\dots,N`. Отримаємо :math:`N-1` TDOA, з яких можна відновити будь-яку іншу попарну різницю, оскільки :math:`\tau_{ij} = \tau_{i1} - \tau_{j1}`. Саме ці :math:`N-1` *незалежних* вимірювань несуть усю геометричну інформацію.
+
+Однак надлишкові пари не є марними. Кожне виміряне TDOA має власний *шум*, тому використання всіх :math:`\binom{N}{2}` пар може поліпшити оцінку. Для цього потрібно правильно врахувати корельовану коваріацію шуму: шум опорного сенсора є спільним для всіх :math:`\tau_{i1}`.
+
+Приклад: двовимірна локалізація за трьома сенсорами
+===================================================
+
+Розташуємо три сенсори в точках
+
+.. math::
+
+ \mathbf{s}_1=(0,0),\quad \mathbf{s}_2=(100,0),\quad \mathbf{s}_3=(0,100)\ \text{(meters)},
+
+а справжнє положення випромінювача нехай дорівнює :math:`\mathbf{u}=(40,30)`. Відстані від випромінювача до сенсорів становлять
+
+.. 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 .
+
+Вибравши сенсор 1 опорним, отримаємо такі різниці дальностей:
+
+.. 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}.
+
+Кожне рівняння задає гіперболу з фокусами :math:`\{\mathbf{s}_2,\mathbf{s}_1\}` і :math:`\{\mathbf{s}_3,\mathbf{s}_1\}` відповідно; точка їхнього перетину є положенням випромінювача. Розв'язувати два рівняння гіпербол вручну незручно. Саме тому далі буде розроблено алгебричну лінеаризацію, за допомогою якої ми відновимо :math:`(40,30)` із цих чисел у замкненій формі.
+
+********************************
+Модель сигналу та вимірювань
+********************************
+
+Модель прийнятого сигналу
+=========================
+
+Кожен сенсор приймає затриману в часі, масштабовану й зашумлену копію сигналу випромінювача. Нехай :math:`s(t)` — форма переданого сигналу. Тоді сенсор :math:`i` приймає
+
+.. math::
+
+ x_i(t) = a_i \, s(t - t_i) + n_i(t), \qquad i = 1,\dots,N,
+
+де :math:`a_i` — дійсний (або комплексний для смугових сигналів) коефіцієнт підсилення, який враховує втрати поширення й характеристику антени; :math:`t_i = t_0 + r_i/c` — абсолютний час прибуття; :math:`n_i(t)` — адитивний шум. Ця модель передбачає один домінантний шлях прямої видимості. Багатопроменеве поширення та відсутність прямої видимості розглянуто в наступному розділі.
+
+Визначення TDOA
+===============
+
+Попарна TDOA — це різниця часу прибуття:
+
+.. 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}.
+
+Права частина явно показує, що TDOA є нелінійною функцією координат випромінювача :math:`\mathbf{u}`. Задача *вимірювання* полягає в оцінюванні :math:`\tau_{ij}` за сигналами :math:`x_i, x_j`, а задача *локалізації* — в оберненні нелінійного відображення від :math:`\mathbf{u}` до набору TDOA.
+
+Припущення щодо шуму
+====================
+
+Припустімо, що кожен :math:`n_i(t)` має нульове середнє, є стаціонарним у широкому сенсі, гаусовим і незалежним від переданого сигналу та шуму інших сенсорів. Відношення сигнал/шум для окремого сенсора дорівнює
+
+.. math::
+
+ \mathrm{SNR}_i = \frac{a_i^2 \sigma_s^2}{\sigma_{n_i}^2},
+
+де :math:`\sigma_s^2` і :math:`\sigma_{n_i}^2` — потужності сигналу та шуму. Це ідеалізовані припущення: реальний шум часто забарвлений і частково корельований між сенсорами. Проте вони приводять до оцінювачів і меж, які добре працюють на практиці, а за потреби цей підхід можна узагальнити на довільну коваріацію шуму.
+
+Нелінійні рівняння вимірювань
+=============================
+
+Зберемо :math:`N-1` різниць дальностей відносно опорного сенсора у вектор :math:`\mathbf{m}` з елементами :math:`m_i = c\,\tau_{i1} = r_i - r_1`. Тоді модель без шуму має вигляд
+
+.. math::
+
+ \mathbf{m} = \mathbf{h}(\mathbf{u}), \qquad
+ h_i(\mathbf{u}) = |\mathbf{u}-\mathbf{s}_i| - |\mathbf{u}-\mathbf{s}_1|,
+
+а зашумлене вимірювання дорівнює :math:`\tilde{\mathbf{m}} = \mathbf{h}(\mathbf{u}) + \boldsymbol{\varepsilon}`, де :math:`\boldsymbol{\varepsilon}` — похибка різниці дальностей, спричинена похибками оцінювання часової затримки. Функція :math:`\mathbf{h}` нелінійна через евклідові норми, і ця нелінійність є джерелом усіх подальших алгоритмічних ускладнень. Є дві основні стратегії: алгебрично *лінеаризувати* задачу введенням допоміжної змінної, як описано в наступному розділі, або *ітеративно* лінеаризувати її біля поточної оцінки, як описано нижче.
+
+****************************************************
+Оцінювання часової затримки: вимірювальний тракт
+****************************************************
+
+Перш ніж використовувати геометрію, необхідно виділити затримки :math:`\tau_{ij}` із необроблених сигналів. Це задача *оцінювання часової затримки* (Time-Delay Estimation, TDE), і її точність зрештою обмежує точність усієї системи.
+
+Взаємна кореляція
+=================
+
+Природний оцінювач використовує те, що :math:`x_i` і :math:`x_j` є зашумленими та зсунутими в часі копіями одного сигналу. Їхня взаємна кореляція
+
+.. math::
+
+ R_{x_i x_j}(\tau) = \mathbb{E}\!\left[ x_i(t)\, x_j(t+\tau) \right],
+
+набуває максимуму, коли зсув :math:`\tau` суміщає дві копії, тобто при :math:`\tau = \tau_{ij}`. Отже, оцінка має вигляд
+
+.. math::
+
+ \hat{\tau}_{ij} = \arg\max_{\tau} \, \hat{R}_{x_i x_j}(\tau).
+
+На практиці кореляцію зазвичай ефективно обчислюють у частотній області за допомогою ШПФ, подібно до великих згорток. Для цього використовують взаємну спектральну густину потужності :math:`G_{x_i x_j}(f) = \mathcal{F}\{R_{x_i x_j}(\tau)\}` і обернене перетворення.
+
+Моделювання в Python
+====================
+
+Досить формул — погляньмо, як усе це працює в простому прикладі Python. Спочатку задаємо основні параметри моделювання: положення випромінювача й сенсорів та частоту дискретизації, яка фактично визначає ширину спектра, що його «бачать» приймачі.
+
+.. 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 мало залежить від конкретного сигналу передавача, хоча його ширина смуги має значення. Щоб спростити приклад, передаватимемо випадковий шум, обмежений заданою смугою. Якби замість нього використали, наприклад, QPSK із такою самою шириною смуги, принципово нічого не змінилося б.
+
+.. 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))
+
+Далі змоделюємо приймання сигналу із затримкою, що залежить від положення кожного приймача. Скористаємося фільтром дробової затримки, з яким ознайомилися в розділі :ref:`sync-chapter`. Решта коду доволі проста. Для кожного приймача обов'язково додаємо окрему реалізацію AWGN.
+
+.. 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
+
+Увесь попередній код стосувався лише моделювання. Решта відтворює дії, які справді виконують для обчислення TDOA — зазвичай у центральному вузлі або на одному із сенсорів, що має доступ до відліків усіх трьох приймачів. Коду небагато: перебираємо всі пари сенсорів, обчислюємо взаємну кореляцію прийнятих сигналів і знаходимо її максимум. Далі розглянемо субдискретний варіант із вищою роздільною здатністю.
+
+.. 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
+
+Ось і все. Отримаємо такий результат:
+
+.. image:: ../_images/tdoa_python_integer.svg
+ :align: center
+ :target: ../_images/tdoa_python_integer.svg
+ :alt: Результат моделювання в Python із цілочисловою кореляцією
+
+Зауважте, що код ще не «розв'язує» задачу повністю. На перший погляд може здаватися інакше, адже лінії перетинаються точно в положенні випромінювача. Проте останній крок розв'язання насправді виконує ваш мозок, коли знаходить перетин гіпербол на графіку. За більшого шуму всі гіперболи також не перетиналися б в одній точці. Автоматизовані способи розв'язання розглянемо далі в цьому розділі.
+
+Повний код Python, включно з побудовою графіків, можна знайти `тут `_.
+
+Роздільна здатність і субдискретне оцінювання
+=============================================
+
+За частоти дискретизації :math:`f_s` кореляцію обчислюють на сітці затримок із кроком :math:`1/f_s`, тому найпростіша роздільна здатність максимуму становить один відлік, або :math:`c/f_s` за дальністю. Зазвичай це надто грубо, особливо коли сенсори й випромінювач розміщені близько, наприклад на відстанях менше 100 метрів. Є два способи субдискретного уточнення. Перший — інтерполювати сигнали під час обчислення взаємної кореляції. Другий — підігнати модель до відліків навколо дискретного максимуму. У другому випадку найпростішою є параболічна інтерполяція за максимумом і двома сусідніми відліками. Точнішу оцінку дає sinc-інтерполяція, адже справжня кореляція сигналу з обмеженою смугою має форму, подібну до sinc. Якісна інтерполяція зазвичай дає оцінки затримки в 10–100 разів точніші за період дискретизації. Нижче показано приклад інтерпольованої взаємної кореляції.
+
+.. 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
+
+Повторивши попереднє моделювання із субдискретним оцінюванням, отримаємо такий результат. Щоб побачити різницю точності, лівий графік доведеться збільшити.
+
+.. image:: ../_images/tdoa_python_subsample.svg
+ :align: center
+ :target: ../_images/tdoa_python_subsample.svg
+ :alt: Результат моделювання в Python із субдискретною кореляцією
+
+На правому графіку добре видно, наскільки значною була похибка початкового цілочислового методу.
+
+Узагальнена взаємна кореляція
+=============================
+
+Звичайна взаємна кореляція нестійка: якщо переданий сигнал вузькосмуговий або канал має багатопроменевість, кореляційний максимум стає широким і легко зміщується під дією шуму. *Узагальнена взаємна кореляція* (Generalized Cross-Correlation, GCC) Кнаппа й Картера розв'язує цю проблему введенням частотної ваги :math:`\Psi(f)` перед зворотним перетворенням в область затримок:
+
+.. 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 .
+
+Вагова функція змінює форму спектра, загострюючи й стабілізуючи максимум. Різні варіанти :math:`\Psi(f)` відповідають класичним оцінювачам, а правильний вибір ваги є основою стійкого TDE. Поширені варіанти:
+
+* **Взаємна кореляція** (:math:`\Psi = 1`): вибір максимальної правдоподібності лише для широкосмугового плаского спектра за високого SNR; в інших умовах неоптимальний.
+* **Roth** (:math:`\Psi = 1/G_{x_i x_i}(f)`): пригнічує частоти, на яких один із сенсорів має багато шуму.
+* **SCOT** (Smoothed Coherence Transform, :math:`\Psi = 1/\sqrt{G_{x_i x_i}G_{x_j x_j}}`): симетричне відбілювання обох каналів.
+* **PHAT** (Phase Transform, :math:`\Psi = 1/|G_{x_i x_j}(f)|`): найпоширеніший в акустиці варіант.
+
+Розгляньмо оцінювач **GCC-PHAT** докладніше. Ділення на модуль взаємного спектра залишає *лише фазу*:
+
+.. 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 .
+
+Затримку між двома копіями сигналу повністю закодовано в члені *лінійної фази* :math:`e^{-j2\pi f \tau_{ij}}`, тоді як модуль містить форму спектра та ревербераційне забарвлення, які часто лише заважають. Відбілювання до одиничного модуля надає всім частотам однакову вагу й утворює гострий, майже імпульсний максимум на справжній затримці. Саме тому PHAT напрочуд стійкий до багатопроменевого поширення. Його слабкість полягає в тому, що він також відбілює частоти, де переважає шум. За низького SNR рівне зважування підсилює шум, тому варіанти з урахуванням SNR повертають вагу на основі когерентності. У реальних системах більшість сигналів передаються не безперервно, тому перед TDOA зазвичай визначають часово-частотні межі цільового сигналу. Якщо завад небагато, це також дає змогу досить легко оцінити SNR.
+
+Практичні міркування
+====================
+
+На точність TDOA впливають кілька чинників:
+
+* **Вікно інтегрування** :math:`T` задає компроміс між дисперсією оцінки (довше вікно краще, оскільки дисперсія спадає приблизно як :math:`1/T`) і припущенням стаціонарності. Для рухомого випромінювача довге вікно також розмиває затримку. Часто його тривалість визначає сам сигнал, наприклад TDOA можна обчислювати окремо для кожного пакета.
+* **Смуга когерентності** обмежує частоти, фаза яких справді придатна для вимірювання.
+* **Ширина смуги сигналу** має вирішальне значення. Як покаже аналіз межі Крамера—Рао, дисперсія затримки зменшується пропорційно *квадрату* ширини смуги. Тому широкосмугові сигнали локалізуються набагато краще за вузькосмугові. Цим TDOA відрізняється від DOA, де ширина смуги не була важливою, а багато методів навіть спиралися на вузькосмугове припущення. Водночас для TDOA не обов'язково приймати весь спектр сигналу. Навіть якщо максимальна частота дискретизації SDR дає змогу захопити лише частину його смуги, TDOA усе одно можна виконати.
+
+З погляду обчислень у TDOA переважають ШПФ зі складністю :math:`O(M\log M)` для кожної пари сенсорів і записів довжиною :math:`M`. Саме це робить практичними великі сенсорні мережі.
+
+Приклад GCC-PHAT у Python
+=========================
+
+Перевага PHAT полягає в тому, що його можна додати до вже створеної моделі майже без нового коду. Наведений вище субдискретний оцінювач уже працював у частотній області: формував взаємний спектр :math:`X_a^*(f)\,X_b(f)`, доповнював його нулями для інтерполяції та виконував обернене ШПФ, щоб отримати кореляцію в області затримок. PHAT додає лише один рядок: перед зворотним перетворенням ділимо взаємний спектр на його модуль. Тоді кожен частотний відлік має одиничну вагу, а зберігається лише фаза, у якій міститься затримка.
+
+.. 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
+
+Єдина відмінність від субдискретного коду — рядок ``X = X / (np.abs(X) + 1e-12)``. Мале значення епсилон у знаменнику не дає частотним відлікам із майже нульовою енергією необмежено зростати під час ділення. Для змодельованого широкосмугового сигналу з високим SNR результат майже не відрізняється від звичайної взаємної кореляції, адже саме в таких умовах PHAT і взаємна кореляція збігаються. Перевага проявляється у складніших умовах: коли спектр забарвлений або багатопроменевість розмиває максимум, відбілювання до одиничного модуля знову стискає кореляцію до гострого, майже імпульсного максимуму на справжній затримці. Саме тому PHAT є типовим вибором в акустиці.
+
+***************************************
+Алгоритми локалізації в замкненій формі
+***************************************
+
+Наведені вище рівняння вимірювань нелінійні й у прямому вигляді потребують ітеративного розв'язання з якісним початковим наближенням. *Оцінювачі в замкненій формі*, тобто неітеративні, обходять цю проблему алгебричним прийомом: вводять допоміжну змінну, яка поглинає нелінійність і перетворює систему на лінійну. Такі методи швидкі, не потребують початкової оцінки й не можуть застрягнути в локальному мінімумі. Тому вони корисні як самостійно, так і для ініціалізації описаних далі ітеративних методів.
+
+Стратегія лінеаризації
+======================
+
+Потрібно піднести рівняння дальностей до квадрата й відняти їх попарно. Це усуває нелінійний член :math:`x^2+y^2` і вводить :math:`r_1`, дальність до опорного сенсора, як єдину допоміжну невідому. Почнімо з квадрата дальності від випромінювача :math:`\mathbf{u}=(x,y)` до сенсора :math:`i` у точці :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 .
+
+Проблемний член :math:`x^2+y^2` однаковий для всіх сенсорів. Виберімо сенсор 1 опорним і віднімемо його рівняння від рівняння сенсора :math:`i`:
+
+.. math::
+
+ r_i^2 - r_1^2 = (K_i - K_1) - 2(x_i-x_1)x - 2(y_i-y_1)y .
+
+Тепер використаємо виміряну різницю дальностей :math:`r_{i1}\equiv r_i - r_1 = c\,\tau_{i1}`. Оскільки :math:`r_i = r_{i1}+r_1`, маємо :math:`r_i^2 = r_{i1}^2 + 2r_{i1}r_1 + r_1^2`, а отже, :math:`r_i^2 - r_1^2 = r_{i1}^2 + 2 r_{i1} r_1`. Після підстановки й перегрупування:
+
+.. 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}
+
+Це рівняння **лінійне** відносно невідомих :math:`(x, y, r_1)`, де дальність до опорного сенсора :math:`r_1` розглядається як допоміжна змінна. Об'єднавши рівняння для :math:`i=2,\dots,N`, отримаємо лінійну систему :math:`\mathbf{A}\boldsymbol{\theta} = \mathbf{b}` з :math:`\boldsymbol{\theta}=[x,y,r_1]^\top`, яку можна розв'язати звичайним або зваженим методом найменших квадратів. Уся нелінійність тепер зосереджена в єдиній додатковій невідомій :math:`r_1`.
+
+Сферична інтерполяція та сферичний перетин
+==========================================
+
+Найперші оцінювачі в замкненій формі — сферична інтерполяція (Spherical Interpolation, SI) та сферичний перетин (Spherical Intersection, SX) Шау й Робінсона — використовують саме цю структуру. Спочатку вони розв'язують лінійну систему для :math:`(x,y)` як функції :math:`r_1`, а потім накладають зв'язувальне обмеження :math:`r_1^2 = (x-x_1)^2+(y-y_1)^2`, щоб визначити :math:`r_1`. Метод SI знаходить :math:`r_1` проєкцією за методом найменших квадратів. SX підставляє лінійний розв'язок у квадратичне обмеження й розв'язує отримане скалярне квадратне рівняння. Обидва методи прості та швидкі, але дещо грубо поводяться з допоміжною змінною, через що за вищого рівня шуму втрачають точність.
+
+Метод Фанга
+===========
+
+Алгоритм Фанга дає точний алгебричний розв'язок для *мінімальної* конфігурації: трьох сенсорів у 2D або чотирьох у 3D. У такому разі система визначена, а не перевизначена. Метод елегантний і майже не потребує обчислень, але не використовує надлишкові сенсори, тому не може усереднювати шум вимірювань і чутливий до геометрії. Його можна розглядати як точно визначений окремий випадок, який методи найменших квадратів узагальнюють.
+
+Щоб зрозуміти його роботу, знову погляньмо на рівняння в рамці. Для трьох сенсорів маємо рівно два такі рівняння (:math:`i=2,3`), але три невідомі :math:`(x,y,r_1)`. Система здається недовизначеною, проте :math:`r_1` не є вільною змінною: вона пов'язана з положенням рівнянням :math:`r_1^2=(x-x_1)^2+(y-y_1)^2`. Прийом Фанга полягає в тому, щоб *відкласти* це обмеження, тимчасово вважати :math:`r_1` відомою сталою й розв'язати два лінійні рівняння відносно :math:`x` і :math:`y`. Оскільки :math:`r_1` входить лінійно, обернення матриці :math:`2\times2`, яка добре обумовлена, якщо сенсори не лежать на одній прямій, дає координати як лінійні функції ще невідомої дальності:
+
+.. math::
+
+ x = g_x + h_x\,r_1, \qquad y = g_y + h_y\,r_1 ,
+
+де сталі :math:`g_x,h_x,g_y,h_y` отримано з оберненої матриці. Тепер скористаємося відкладеним обмеженням. Підстановка цих виразів у :math:`r_1^2=(x-x_1)^2+(y-y_1)^2` зводить усе до одного скалярного квадратного рівняння :math:`a\,r_1^2 + b\,r_1 + c = 0`. Розв'язуємо його, залишаємо фізично допустимий корінь — дальність має бути додатною, а інший корінь зазвичай відповідає неправильній гілці гіперболи — і зворотною підстановкою знаходимо :math:`(x,y)`. Ось і весь метод: одне розв'язання системи :math:`2\times2`, одне квадратне рівняння, жодних ітерацій або початкового наближення. У наступному підрозділі цю процедуру буде виконано з конкретними числами.
+
+Метод Чана: двоетапний зважений метод найменших квадратів
+=========================================================
+
+Практичним стандартом став двоетапний зважений метод найменших квадратів (Weighted Least Squares, WLS) Чана й Хо. Він спирається на наведену вище лінійну систему, але правильно враховує статистику похибок і уточнює допоміжну змінну. За малого та помірного шуму його точність наближається до межі Крамера—Рао.
+
+**Перший етап.** Тимчасово вважатимемо три компоненти :math:`\boldsymbol{\theta}=[x,y,r_1]^\top` незалежними й розв'яжемо лінійну систему зваженим методом найменших квадратів:
+
+.. math::
+
+ \hat{\boldsymbol{\theta}} = (\mathbf{A}^\top \mathbf{W}\mathbf{A})^{-1}\mathbf{A}^\top \mathbf{W}\,\mathbf{b},
+
+де вагу :math:`\mathbf{W}` вибирають як обернену коваріацію похибок рівнянь. Оскільки ця коваріація сама залежить від невідомих дальностей, на практиці спочатку розв'язують систему з :math:`\mathbf{W}=\mathbf{I}` або з початковою коваріацією шуму TDOA. Потім за отриманими дальностями повторно обчислюють :math:`\mathbf{W}` і розв'язують систему ще раз. Зазвичай достатньо одного-двох таких уточнень.
+
+**Другий етап.** На першому етапі було проігноровано відомий зв'язок :math:`r_1^2 = (x-x_1)^2+(y-y_1)^2` між допоміжною змінною та положенням. Другий етап відновлює його: складається невелика нова задача найменших квадратів для квадратів :math:`[(x-x_1)^2,(y-y_1)^2,r_1^2]`, зважена коваріацією першого етапу, і знаходиться виправлене положення. Саме другий WLS усуває значну частину зміщення наївного лінійного розв'язку й наближає оцінювач Чана до оптимального.
+
+Метод одразу повертає положення, а основні витрати зводяться до обернення малих матриць :math:`3\times3`, що є незначним порівняно з ШПФ у вимірювальному тракті. Обмеження проявляються за сильного шуму або несприятливої геометрії: операції з квадратами дальностей підсилюють похибки, а на другому етапі можна вибрати неправильний корінь. Типовий спосіб виправлення — використати результат Чана як початкове наближення для описаного нижче ітеративного уточнення.
+
+Продовження прикладу: розв'язання задачі з трьома сенсорами
+===========================================================
+
+Повернімося до місця, де завершилося моделювання в Python. Масив ``range_diff`` уже містить одну виміряну різницю дальностей для кожної пари сенсорів. Раніше останній крок виконувався візуально — ми дивилися, де перетинаються гіперболи. Тепер замінимо це алгебричним розв'язком у замкненій формі й безпосередньо відновимо положення випромінювача за ``range_diff`` і ``rx_positions``. Оскільки в 2D маємо рівно три сенсори, це мінімальна конфігурація Фанга: два лінійні рівняння в рамці й одне квадратне рівняння, без ітерацій та початкового наближення.
+
+Виберімо ``Rx0`` опорним сенсором. Пари сформовано як ``(0,1)``, ``(0,2)``, ``(1,2)``. Нагадаємо, що ``range_diff[k]`` для пари ``(a,b)`` дорівнює :math:`r_b-r_a`. Тому пари з опорним сенсором, ``(0,1)`` і ``(0,2)``, безпосередньо дають різниці :math:`r_{i0}=r_i-r_0`, потрібні для рівняння в рамці.
+
+.. 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)
+
+Структура коду точно повторює математику. ``M`` і ``d`` задають два лінійні рівняння в рамці; ``g`` і ``h`` виражають :math:`x` та :math:`y` як лінійні функції ще невідомої опорної дальності :math:`r_1`, яка тут називається ``r_ref``. Підстановка цих функцій у :math:`r_1^2=(x-x_1)^2+(y-y_1)^2` зводить задачу до скалярного квадратного рівняння, яке розв'язує ``np.roots``. Відкидаємо нефізичний, тобто від'ємний або комплексний, корінь, залишаємо додатний дійсний і зворотною підстановкою знаходимо положення. У нашій широкосмуговій моделі з високим SNR оцінка точно збігається зі справжнім положенням :math:`(153, 355)` без ручного пошуку перетину гіпербол.
+
+За зашумленіших вимірювань два лінійні рівняння вже не були б цілком узгодженими, а корінь квадратного рівняння змістився б. Три сенсори не дають надлишковості для усереднення, тому похибка безпосередньо перейшла б у результат. Саме тут корисними стають надлишкові пари, зважування та другий етап методу Чана: вони визначають, наскільки плавно погіршуватиметься оцінка.
+
+****************************************
+Ітеративне та статистичне оцінювання
+****************************************
+
+Методи в замкненій формі швидкі, але використовують алгебричні наближення, які знижують точність за сильного шуму або несприятливої геометрії. Коли потрібна найкраща можлива оцінка, нелінійну задачу розв'язують безпосередньо, зазвичай починаючи з результату методу в замкненій формі.
+
+Нелінійний метод найменших квадратів
+====================================
+
+Визначимо залишок між виміряними й передбаченими різницями дальностей і мінімізуємо його зважену квадратичну норму:
+
+.. 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],
+
+де :math:`\mathbf{C}` — коваріація похибок різниці дальностей. Через нелінійність :math:`\mathbf{h}` ця цільова функція не має розв'язку в замкненій формі, тому її мінімізують ітеративно.
+
+Метод ряду Тейлора (Гаусса—Ньютона)
+===================================
+
+Класичний підхід Фоя лінеаризує :math:`\mathbf{h}` біля поточної оцінки :math:`\mathbf{u}^{(k)}` за допомогою якобіана :math:`\mathbf{J}`. Його рядок :math:`i` є градієнтом :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,
+
+тобто різницею *одиничних векторів*, спрямованих від передбачуваного положення випромінювача до сенсора :math:`i` та опорного сенсора. Крок Гаусса—Ньютона має вигляд
+
+.. 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],
+
+і повторюється до збіжності. На кожному кроці розв'язується мала лінійна система. Метод швидко збігається, *якщо початкова точка близька до розв'язку*. Саме тому оцінка Чана є бажаним початковим наближенням: вона потрапляє в область притягання глобального мінімуму й допомагає уникнути хибних локальних мінімумів гіперболічної цільової поверхні, особливо за несприятливої геометрії.
+
+Оцінювання методом максимальної правдоподібності
+=================================================
+
+За гаусового шуму від'ємний логарифм правдоподібності з точністю до сталих збігається з наведеним вище зваженим квадратом залишку. Отже, **оцінювач максимальної правдоподібності збігається зі зваженим нелінійним методом найменших квадратів**. Ітерація Гаусса—Ньютона — не евристика, а статистично оптимальний оцінювач для прийнятої моделі. Саме його коваріацію передбачає наведена далі межа Крамера—Рао.
+
+Продовжимо приклад Python. Уже маємо положення ``emitter_est`` від розв'язувача в замкненій формі. Теорія підказує, що оцінка максимальної правдоподібності є наведеною вище ітерацією Гаусса—Ньютона, а розв'язок у замкненій формі є ідеальною початковою точкою, бо лежить в області справжнього мінімуму. Почнемо з ``emitter_est`` і виконаємо кілька кроків. На кожному з них повторно лінеаризуємо модель різниць дальностей у поточній точці й розв'язуємо малу задачу найменших квадратів для поправки. На відміну від методу Фанга, який використовував лише дві пари з опорним сенсором, цей алгоритм бере всі три елементи ``range_diff``. Додаткова пара створює надлишковість, яку ітерація використовує для усереднення.
+
+.. 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)
+
+Зверніть увагу на кілька деталей. Ми припустили, що похибки різниць дальностей незалежні й мають однакову дисперсію, тому вага :math:`\mathbf{C}^{-1}=\sigma^{-2}\mathbf{I}` є скаляром, який скорочується в рівнянні поправки. Саме тому звичайний ``np.linalg.lstsq`` без матриці ваг точно обчислює крок. Якби пари мали різну якість, сюди потрібно було б включити обернені дисперсії. Рядки якобіана буквально є різницями одиничних векторів напрямку ``e_b - e_a`` з наведеного вище рівняння, тож у коді безпосередньо видно вплив геометрії. Починаючи з якісної оцінки в замкненій формі, ітерація збігається за кілька кроків до справжнього положення :math:`(153, 355)`. У моделі з високим SNR вона майже не змінює початкову відповідь, але за сильнішого шуму саме додаткова пара та ітеративне уточнення дають перевагу. Та сама матриця :math:`\mathbf{J}^\top\mathbf{C}^{-1}\mathbf{J}` далі з'явиться в межі Крамера—Рао як коваріація оцінювача.
+
+Стійкі, рекурсивні та баєсівські розширення
+===========================================
+
+Реальні вимірювання містять викиди: TDOA, спотворене багатопроменевістю, може бути дуже хибним, тоді як решта значень правильні. Звичайний метод найменших квадратів підносить залишки до квадрата, тому викиди сильно викривляють результат. *Стійкі* оцінювачі замінюють квадратичну функцію втрат на повільніше зростаючу, наприклад функцію Г'юбера, або явно знаходять і відкидають неузгоджені TDOA за тестами залишків чи консенсусом на зразок RANSAC.
+
+Коли випромінювач *рухається*, доцільно поєднувати вимірювання в часі, а не локалізувати кожен момент незалежно. Фільтрація в просторі станів моделює положення та швидкість як стан, що змінюється. **Фільтр Калмана** оптимальний для лінійної гаусової динаміки, але вимірювання TDOA нелінійне. Тому застосовують **розширений фільтр Калмана**, який лінеаризує вимірювання тим самим якобіаном; **ансцентний фільтр Калмана**, який проводить детермінований набір сигма-точок через нелінійність, не потребує явного якобіана й краще працює із сильнішою нелінійністю; або, для мультимодальних чи істотно негаусових задач, **частинковий фільтр**, що представляє апостеріорний розподіл зваженою хмарою відліків. Такі трекери природно забезпечують неперервність руху й пригнічують неоднозначності окремих статичних оцінок.
+
+*********************************
+Повний перебір із тепловою картою
+*********************************
+
+Усі попередні методи були алгебричними або ітеративними: ми перетворювали рівняння чи рухалися вздовж градієнта. Проте є напрочуд проста альтернатива, яка не потребує ні того, ні іншого. Накладемо сітку на область пошуку й у кожній можливій точці поставимо одне запитання: *якби випромінювач був тут, які різниці дальностей побачили б сенсори та наскільки вони відрізнялися б від виміряних?* Сума квадратів невідповідностей задає вартість кожної точки, а випромінювач розташований там, де ця вартість найменша. Результатом є теплова карта тієї самої цільової поверхні, уздовж якої непомітно рухалася ітерація Гаусса—Ньютона, але тепер ми бачимо всю поверхню одразу.
+
+У прикладі Python для цього вже є всі потрібні змінні: ``range_diff``, ``rx_positions`` і ``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)
+
+Зображення ``likelihood`` безпосередньо показує геометрію: яскраві гребені повторюють розглянуті раніше гіперболи й сходяться в одному яскравому максимумі на справжньому положенні випромінювача. Беремо від'ємний логарифм вартості, щоб найімовірніше положення було максимумом, а не мінімумом, — так графік легше читати. Компроміси цілком очікувані. Метод надзвичайно простий, не потребує початкового наближення, не може розбігтися чи вибрати неправильний корінь, тому є чудовою перевіркою та надійним способом *ініціалізувати* ітеративний оцінювач. Він також природно працює з мультимодальними поверхнями, бо бачить усі мінімуми, а не лише найближчий. Ціна — роздільна здатність і швидкість: точність обмежена кроком сітки, а обсяг обчислень зростає з кількістю її комірок. Для точного пошуку у великій області спочатку виконують грубу локалізацію, а потім уточнюють результат, збільшуючи масштаб сітки або передаючи оцінку методу Гаусса—Ньютона. Нижче показано теплову карту для нашого прикладу.
+
+.. image:: ../_images/tdoa_python_heatmap.svg
+ :align: center
+ :target: ../_images/tdoa_python_heatmap.svg
+ :alt: Теплова карта на раніше побудованому графіку TDOA
+
+Перевага теплової карти полягає ще й у тому, що за великої похибки або непоміченого низького SNR деяких сенсорів на ній можуть виникнути кілька гарячих ділянок — людина легко це помітить. Карту навіть можна накласти на супутникове зображення місцевості.
+
+Такий повний перебір обчислювально неефективний, а для тривимірного TDOA майже непридатний. Альтернатива, яка також залишається перебором, але не обчислює кожну точку сітки, — побудувати всі гіперболи в 2D із певною «шириною», наприклад застосувати вздовж кожної гіперболи пелюсткоподібну функцію зі спаданням убік.
+
+******************************************
+Аналіз характеристик і фундаментальні межі
+******************************************
+
+Маючи оцінювачі, поставимо два запитання: якої точності *в принципі* може досягти система TDOA і що цю точність визначає? Відповідь дають дві ідеї: межа Крамера—Рао задає обумовлену сигналом нижню межу шуму, а геометричний фактор погіршення точності описує, як взаємне розташування сенсорів і випромінювача підсилює цю межу.
+
+Поширення похибок
+=================
+
+Точність системи формується у два етапи. Спочатку скінченні SNR і ширина смуги обмежують точність вимірювання кожної затримки, утворюючи похибку TDE. Потім геометрія перетворює похибки різниць дальностей на похибку положення. Позначимо похибку положення як :math:`\delta\mathbf{u}`, а похибки різниць дальностей як :math:`\boldsymbol{\varepsilon}`. Лінеаризований зв'язок поблизу розв'язку має вигляд :math:`\boldsymbol{\varepsilon}\approx \mathbf{J}\,\delta\mathbf{u}`, тому коваріація похибки положення дорівнює
+
+.. math::
+
+ \mathrm{Cov}(\hat{\mathbf{u}}) \approx (\mathbf{J}^\top \mathbf{C}^{-1}\mathbf{J})^{-1}.
+
+Цей один вираз містить обидва етапи: :math:`\mathbf{C}` характеризує якість вимірювання TDE, а :math:`\mathbf{J}` — геометрію.
+
+Межа оцінювання часової затримки
+================================
+
+Можна обмежити точність будь-якого оцінювача однієї затримки. Для сигналу із середньоквадратичною шириною смуги :math:`\beta`, який спостерігається протягом часу :math:`T`, дисперсія будь-якої незміщеної оцінки затримки задовольняє нерівність
+
+.. math::
+
+ \mathrm{var}(\hat\tau_{ij}) \gtrsim \frac{1}{8\pi^2 \beta^2 T \gamma},
+
+де :math:`\beta` — середньоквадратична, або габорівська, ширина смуги сигналу, а :math:`\gamma` — ефективний коефіцієнт SNR, який поєднує SNR двох сенсорів. Звідси випливають три правила проєктування: дисперсія зменшується зі збільшенням **часу інтегрування** :math:`T`, **ефективного SNR** :math:`\gamma` і, що особливо важливо, **квадрата ширини смуги** :math:`\beta^2`. Подвоєння ширини смуги зменшує дисперсію затримки вчетверо. Саме тому широкосмугові сигнали та сигнали з розширеним спектром цінні для вимірювання дальності, а вузькосмугові випромінювачі за своєю природою складно локалізувати лише методом TDOA.
+
+Нижня межа Крамера—Рао для локалізації
+======================================
+
+Поєднавши якість вимірювань і геометрію, отримаємо матрицю інформації Фішера для положення випромінювача:
+
+.. math::
+
+ \mathbf{F} = \mathbf{J}^\top \mathbf{C}^{-1} \mathbf{J}.
+
+Нижня межа Крамера—Рао стверджує, що коваріація *будь-якого* незміщеного оцінювача не може бути меншою за обернену матрицю:
+
+.. math::
+
+ \mathrm{Cov}(\hat{\mathbf{u}}) \succeq \mathbf{F}^{-1} = (\mathbf{J}^\top \mathbf{C}^{-1}\mathbf{J})^{-1}.
+
+Ця межа є еталоном для порівняння оцінювачів: метод, який її досягає, називають *ефективним*. Наведений вище оцінювач максимальної правдоподібності асимптотично досягає межі за великого :math:`T` і високого SNR, а метод Чана — за малого шуму. Саме тому обидва широко застосовують. CRLB також чітко розділяє два чинники точності: :math:`\mathbf{C}` описує якість сигналу й шуму, яку можна поліпшити ширшою смугою, більшою потужністю або довшим інтегруванням; :math:`\mathbf{J}` описує геометрію, яку поліпшують розташуванням сенсорів. Нижче для кількох значень ширини смуги показано нижню межу залежно від SNR. Це дає уявлення про очікувану похибку або принаймні її мінімум. Вісь ординат показує значення :math:`1\sigma`, тобто одне стандартне відхилення.
+
+.. image:: ../_images/tdoa_cramer_rao.svg
+ :align: center
+ :target: ../_images/tdoa_cramer_rao.svg
+ :alt: Графік нижньої межі Крамера—Рао
+
+Геометричний фактор погіршення точності
+=======================================
+
+Припустімо, сенсори вимірюють різниці дальностей із точністю приблизно 1 м — цілком пристойний результат для добре синхронізованої радіосистеми. Можна було б очікувати, що положення випромінювача також визначатиметься приблизно до 1 м. Але де саме він розташований? Уявіть випромінювач усередині трикутника з трьох сенсорів. Гіперболи для кожної пари перетинаються під великими, майже прямими кутами, тому точка перетину визначена чітко, а метрова похибка дальності перетворюється, наприклад, на 1,5 м похибки положення. Тепер перемістімо той самий випромінювач далеко вбік, за межі групи сенсорів. Гіперболи торкаються одна одної під малим кутом, немов дві криві, які майже збігаються, а точка перетину розмивається вздовж їхнього спільного напрямку. Та сама метрова похибка може перетворитися на десятки метрів похибки положення. Апаратне забезпечення не змінилося — змінилася лише геометрія.
+
+Цей коефіцієнт підсилення називають **геометричним фактором погіршення точності** (Geometric Dilution of Precision, GDOP). Він показує, наскільки розташування сенсорів і випромінювача збільшує похибку вимірювання під час перетворення на похибку положення. Якщо похибки різниць дальностей незалежні й мають однакове стандартне відхилення :math:`\sigma`, тобто :math:`\mathbf{C}=\sigma^2\mathbf{I}`, де на діагоналі стоять :math:`\sigma^2`, тоді
+
+.. math::
+
+ \mathrm{GDOP} = \sqrt{\mathrm{tr}\bigl[(\mathbf{J}^\top\mathbf{J})^{-1}\bigr]}, \qquad
+ \sigma_{\text{position}} = \mathrm{GDOP}\cdot \sigma .
+
+Похибка положення дорівнює похибці дальності, помноженій на GDOP. Отже, GDOP — безрозмірне число, завжди :math:`\ge 1`, яке показує, у скільки разів похибка дальності збільшується в заданому положенні випромінювача.
+
+Звідки береться це підсилення? Воно закладене в якобіані :math:`\mathbf{J}`, рядки якого є різницями одиничних векторів напрямку :math:`\hat{\mathbf{e}}_i - \hat{\mathbf{e}}_1`, тобто напрямку на один сенсор мінус напрямок на інший. Коли ці напрямки добре рознесені, :math:`\mathbf{J}^\top\mathbf{J}` добре обумовлена — далека від виродження, тож її обернена матриця залишається малою — і GDOP невеликий. Коли напрямки майже збігаються, матриця стає майже виродженою, а GDOP різко зростає. Тому для випромінювача, оточеного сенсорами, вектори напрямків рознесені, гіперболи перетинаються під великими кутами й GDOP малий. Якщо випромінювач далеко за межами групи або сенсори майже колінеарні, вектори майже паралельні, гіперболи перетинаються під малими кутами й GDOP дуже великий.
+
+Це той самий ефект, який спостерігався під час виродження гіпербол біля кінців бази. Головний висновок: система TDOA може значно сильніше обмежуватися тим, *де розміщені її сенсори*, ніж тим, *наскільки точно вона вимірює час*. Навіть наносекундна синхронізація та широка смуга не врятують локалізацію в області з великим GDOP.
+
+На рисунку нижче показано теплові карти GDOP на площині: ліворуч три сенсори розташовано у вершинах рівностороннього трикутника, праворуч — майже на одній прямій. Усередині трикутника видно широку область малого GDOP, тоді як колінеарна решітка має лише вузький придатний коридор. В обох випадках за межами опуклої оболонки GDOP швидко зростає.
+
+.. image:: ../_images/tdoa_gdop.svg
+ :align: center
+ :target: ../_images/tdoa_gdop.svg
+ :alt: Теплові карти GDOP для трикутного й майже колінеарного розташування трьох сенсорів
+
+Оптимізація розташування сенсорів
+=================================
+
+Оскільки геометрію часто можна вибирати під час проєктування, сенсори розміщують так, щоб мінімізувати похибку. Поширені критерії мінімізують скаляр, отриманий із :math:`\mathbf{F}^{-1}`: слід матриці, еквівалентний GDOP; визначник, що відповідає об'єму еліпса довіри; або найбільше власне значення, яке характеризує найгіршу похибку. Якісні висновки інтуїтивні: сенсори слід широко розносити, щоб довгі бази підвищували кутову роздільну здатність; оточувати ними область інтересу, щоб випромінювачі потрапляли всередину опуклої оболонки; уникати колінеарних і компланарних конфігурацій із погано обумовленими напрямками; додавати сенсори там, де надлишковість одночасно зменшує дисперсію та захищає від викидів. Для рухомої цілі або великої території розташування оптимізують по всій області, зазвичай чисельним пошуком мінімуму середнього чи найгіршого GDOP.
+
+****************************************
+Практичні проблеми реальних систем
+****************************************
+
+Використана досі модель не враховувала кількох ефектів, які зазвичай переважають у балансі похибок реальної системи TDOA. Три з них потребують докладного розгляду.
+
+Синхронізація приймачів
+=======================
+
+Головна перевага TDOA — відсутність потреби в синхронізованому передавачі — нерозривно пов'язана з головною складністю: *приймачі* мають спільно використовувати точну опору часу, а будь-яка її похибка безпосередньо входить у вимірювання. Якщо годинник сенсора :math:`i` має зміщення :math:`\delta t_i`, виміряне TDOA містить похибку :math:`\delta t_i - \delta t_j`, яка під час перетворення на дальність множиться на :math:`c`. Для радіосистем масштаб дуже жорсткий:
+
+.. math::
+
+ c \times 1\ \text{ns} = (3\times10^8\,\text{m/s})(10^{-9}\,\text{s}) = 0.30\ \text{m}.
+
+Отже, похибка синхронізації 1 нс уже дає приблизно :math:`\sim`\0,3 м, а 100 нс — 30 м. Досягнення й утримання наносекундної синхронізації розподілених сенсорів є центральним завданням проєктування. Поширені засоби включають генератори, дисципліновані GPS (кожен сенсор отримує від супутників опору часу з точністю приблизно 10–100 нс), протокол точного часу PTP (IEEE 1588), який розподіляє час мережею з точністю менше мікросекунди, а з апаратними часовими мітками — менше 100 нс, і White Rabbit для найвимогливіших систем, що забезпечує субнаносекундну синхронізацію оптоволокном. Важливі ще дві тонкощі. Годинник має не лише сталий *зсув*, а й *дрейф* у часі, тому його потрібно неперервно дисциплінувати. В акустичних системах швидкість :math:`c` приблизно в мільйон разів менша, тому та сама абсолютна похибка часу в мільйон разів менш шкідлива. Саме тому TDOA мікрофонних решіток порівняно невибагливий, тоді як успіх радіочастотного TDOA повністю залежить від годинників.
+
+Серед готових SDR, які легко синхронізувати, є пристрої Ettus Research USRP із підтримкою генератора, дисциплінованого GPS (GPSDO), наприклад `B200 `_ із модулем `TCXO `_, що є різновидом GPSDO. Якщо сенсори розташовані достатньо близько, їм можна подати спільний сигнал PPS кабелем, наприклад від `OctoClock `_. Він також формує сигнал 10 МГц для частотної синхронізації. Майже всі USRP мають входи PPS і 10 МГц, а більшість підтримують установлення GPSDO або вже комплектуються ним.
+
+Багатопроменеве поширення та відсутність прямої видимості
+=========================================================
+
+Досі ми припускали єдиний шлях прямої видимості між випромінювачем і приймачами. У реальному середовищі виникають відбиття, тобто багатопроменевість, а прямий шлях іноді повністю перекритий. Багатопроменеве поширення накладає затримані копії сигналу, які спотворюють або розщеплюють кореляційний максимум і зміщують оцінку затримки. Саме проти цієї проблеми створено GCC-PHAT: відбілювання загострює максимум прямого шляху відносно розмитих відбиттів. Якщо прямий шлях повністю перекритий, навіть *найраніша* енергія проходить зайву відстань. Виміряне TDOA отримує систематичне завищення, якого не усуне жодне усереднення. Методи протидії включають статистичне виявлення каналів без прямої видимості — такі вимірювання часто мають більшу дисперсію або порушують геометричну узгодженість надлишкових сенсорів — зменшення їхньої ваги чи відкидання, для чого потрібно значно більше трьох сенсорів, і використання надлишковості, щоб стійкі оцінювачі могли виявити та відкинути кілька пошкоджених зв'язків. У щільному приміщенні з багатопроменевістю, яке є надзвичайно складним для TDOA, оцінювання затримки на основі моделі та методи машинного навчання дедалі частіше перевершують класичну кореляцію.
+
+Невизначеність положення сенсорів і калібрування
+================================================
+
+У геометричній моделі координати сенсорів :math:`\mathbf{s}_i` вважалися точними. Їхні похибки переходять в оцінку положення так само, як похибки вимірювання, а для віддалених випромінювачів можуть підсилюватися тією самою несприятливою геометрією, що збільшує GDOP. Типовими рішеннями є точна геодезична прив'язка стаціонарних установок, GPS-позиціювання рухомих сенсорів і *самокалібрування* — спільне оцінювання положень сенсорів та випромінювачів за сторонніми випромінювачами з відомими або обмеженими положеннями. Повний баланс похибок має враховувати невизначеність координат сенсорів разом із похибками часу та TDE. У добре синхронізованих системах це часто наступна за величиною складова.
+
+*******************
+Додаткові теми
+*******************
+
+Спільне оцінювання TDOA/FDOA
+============================
+
+Коли випромінювач, сенсори або всі вони *рухаються*, відносний рух створює різний доплерівський зсув у різних сенсорах — **різницю частоти прибуття** (Frequency Difference of Arrival, FDOA). FDOA несе інформацію про *швидкість* випромінювача й, що особливо важливо, додає незалежне геометричне обмеження, яке покращує спостережуваність положення. Це корисно у складних випадках далекого поля та малої кількості сенсорів, де TDOA погано обумовлений. TDOA й FDOA спільно оцінюють максимізацією **комплексної функції невизначеності** (Complex Ambiguity Function, CAF) за затримкою та частотним зсувом:
+
+.. math::
+
+ A(\tau,\nu) = \int_0^T x_i(t)\, x_j^{*}(t-\tau)\, e^{-j2\pi \nu t}\, dt,
+
+двовимірний максимум якої одночасно дає :math:`(\hat\tau_{ij},\hat\nu_{ij})`. CAF узагальнює наведену вище взаємну кореляцію, додаючи вимір частотного пошуку ціною більших обчислювальних витрат. Зауважте, що це не та сама CAF, яку в розділі про циклостаціонарність було введено як циклічну автокореляційну функцію. Спільна обробка TDOA/FDOA лежить в основі супутникової та повітряної геолокації радіовипромінювачів: одна пара рухомих платформ може визначити положення нерухомого випромінювача за спільними обмеженнями затримки й доплерівського зсуву.
diff --git a/index-ukraine.rst b/index-ukraine.rst
index 4e24f6fc..17f166e9 100644
--- a/index-ukraine.rst
+++ b/index-ukraine.rst
@@ -31,8 +31,11 @@
content-ukraine/doa
content-ukraine/2d_beamforming
content-ukraine/phaser
- content-ukraine/cyclostationar
+ content-ukraine/cyclostationary
content-ukraine/pyqt
+ content-ukraine/detection
+ content-ukraine/fpv_video
+ content-ukraine/tdoa
content-ukraine/about_author
.. raw:: html