From f0fa7e8a46061d015b4c0c5ef9aa0ef7082276fd Mon Sep 17 00:00:00 2001 From: LittleSam129 <1qwe432@gmail.com> Date: Wed, 19 Aug 2026 18:14:13 +0300 Subject: [PATCH] Lab043: validate pilot-aided short BPSK link --- docs/lab043_pluto_to_rtlsdr_spec.md | 521 +++++++ experiments/lab043_pluto_to_rtlsdr.py | 1832 +++++++++++++++++++++++++ experiments/lab043_prbs11_hardware.py | 485 +++++++ protocol/bpsk_radio.py | 531 +++++++ tests/test_bpsk_radio.py | 156 +++ tests/test_lab043_calibration.py | 891 ++++++++++++ 6 files changed, 4416 insertions(+) create mode 100644 docs/lab043_pluto_to_rtlsdr_spec.md create mode 100644 experiments/lab043_pluto_to_rtlsdr.py create mode 100644 experiments/lab043_prbs11_hardware.py create mode 100644 tests/test_lab043_calibration.py diff --git a/docs/lab043_pluto_to_rtlsdr_spec.md b/docs/lab043_pluto_to_rtlsdr_spec.md new file mode 100644 index 0000000..a975e61 --- /dev/null +++ b/docs/lab043_pluto_to_rtlsdr_spec.md @@ -0,0 +1,521 @@ +# Lab043. Передача JPEG с Pluto+ на независимый RTL-SDR + +## 1. Цель лабораторной + +Передать тот же подготовленный JPEG, который использован в Lab042, с +Pluto+ на отдельный RTL-SDR и восстановить его побайтово. + +В Lab042 передатчик и приёмник находились в одном Pluto+ и использовали +общий опорный генератор. Lab043 должна убрать это упрощение и раздельно +измерить три эффекта: + +1. рассогласование несущих частот передатчика и приёмника; +2. относительное рассогласование частот дискретизации двух устройств; +3. начальную фазу выбора отсчёта символа внутри символьного интервала. + +Лабораторная считается успешной только при выполнении радиокритериев, +прикладных критериев и программных проверок, определённых в разделе 11. +Запуск программы сам по себе успехом не считается. + +## 2. Границы работы + +Lab043 остаётся кабельной. Эфир добавил бы многолучёвость, внешние сигналы +и неконтролируемые потери, хотя предмет этой работы состоит в +синхронизации независимых устройств. + +В Lab043 не входят: + +- передача через антенны; +- выбор рабочего диапазона 200-250 МГц; +- ARQ и повторная передача; +- исправление стираний и новая FEC; +- потоковое видео; +- скачки по частоте; +- изменение логики управления и безопасности ровера. + +Эфирный тракт остаётся предметом следующей лабораторной. + +## 3. Подтверждённое оборудование и окружение + +Состояние оборудования проверено непосредственным опросом. Перед +реализацией и перед каждым аппаратным прогоном его необходимо проверить +повторно. + +### 3.1. Передатчик + +- Pluto+; +- устройство должно открываться успешным `iio.Context`, а не проверяться + только командой `ping`; +- чип AD9361; +- используется первый передающий канал. + +Конкретный рабочий URI IIO определяется повторным опросом перед опытом и +сохраняется в метаданных запуска. Исторически работали сетевой URI +`ip:192.168.2.1` и прямой USB IIO, но ни один из них нельзя считать +доступным без текущей проверки. + +### 3.2. Приёмник + +- USB ID `VID_0BDA&PID_2838`; +- имя библиотеки: `Generic RTL2832U OEM`; +- тюнер, сообщённый библиотекой: Fitipower FC0013; +- основной интерфейс `MI_00` использует WinUSB от libwdi; +- устройство открыто из Python через `pyrtlsdr` и `librtlsdr`; +- установка 435 МГц подтверждена чтением значения обратно; +- частота дискретизации 2 400 000 отсчётов/с работает; +- получены комплексные отсчёты; +- доступные аппаратные усиления находятся в диапазоне от -9,9 до + 19,7 дБ и задаются дискретными ступенями устройства. + +Ошибка второго USB-интерфейса не считается отказом SDR. Основной +интерфейс `MI_00`, через который работает `librtlsdr`, исправен. + +### 3.3. Зафиксированное программное окружение + +- Python 3.13.1 x64; +- `pyrtlsdr` 0.5.0; +- `pyrtlsdrlib` 0.0.5; +- поставляемая пакетом нативная `librtlsdr` v0.9.0; +- `pyadi-iio` и `pylibiio` для Pluto+. + +Принятый способ доступа к RTL-SDR: API `pyrtlsdr` с поставляемой +`pyrtlsdrlib` нативной библиотекой. Малый собственный DLL-адаптер не +нужен. Реализация должна проверять версии и выдавать понятную ошибку, если +библиотека или устройство недоступны. Пользовательский путь к DLL +зашивать в код запрещено. + +## 4. Физическая схема и безопасность + +Схема: + +`Pluto+ TX -> SMA-кабель -> AT30S 30 дБ -> переходник/кабель -> RTL-SDR RX` + +Обязательные условия перед каждым включением TX: + +1. Антенны Pluto+ и RTL-SDR сняты. +2. AT30S включён последовательно в тракт, а не подключён к свободному + порту. +3. Все разъёмы затянуты, кабель не отсоединяется при включённом TX. +4. Пользователь отдельно подтверждает эту физическую схему. +5. Первый заранее определённый рабочий режим: TX -30 дБ, RTL-SDR + 19,7 дБ, несущая 435 МГц, номинально 2,4 Мвыб/с на обоих устройствах. +6. Автоматический перебор усилений, частот, порогов и способов + синхронизации запрещён. +7. Если первый режим не проходит, программа сохраняет исходные данные и + диагностику и останавливается. Следующий прогон разрешён только после + разбора причины и с изменением ровно одного заранее названного + параметра. + +При TX -30 дБ расчётный уровень после AT30S составляет около -53 дБм. +Это значение следует из бюджета Lab042: около +7 дБм при TX 0 дБ, +минус 30 дБ настройки TX и минус 30 дБ аттенюатора. Расчёт не заменяет +проверку наличия аттенюатора и правильности схемы. + +## 5. Номинальные параметры сигнала и обязательное чтение обратно + +| Параметр | Запрошенное значение | Происхождение | +|---|---:|---| +| Несущая TX и RX | 435 МГц | Кабельная рабочая точка Lab042 | +| Символьная скорость | 20 000 симв/с | Lab018-Lab023 и Lab042 | +| Частота дискретизации TX | 2 400 000 отсчётов/с | Поддерживается Pluto+ и проверена на RTL-SDR | +| Частота дискретизации RX | 2 400 000 отсчётов/с | Проверена непосредственным опросом RTL-SDR | +| Номинальное число отсчётов на символ | 120 | 2 400 000 / 20 000 | +| Полоса TX | 200 кГц | Минимальная полоса AD9361, использованная в Lab042 | +| Размер данных фрагмента | 512 байт | Lab042 | +| Модуляция | BPSK | Lab018-Lab023 и Lab042 | + +Одинаковое номинальное значение 2,4 Мвыб/с выбрано специально. Lab043 +измеряет физическую ошибку генераторов, а не искусственно созданную +разницу номинальных частот дискретизации. + +Перед аппаратным опытом программа обязана установить и прочитать обратно: + +- частоту дискретизации TX Pluto+; +- частоту дискретизации RX RTL-SDR; +- центральную частоту TX Pluto+; +- центральную частоту RX RTL-SDR; +- усиление TX Pluto+; +- фактически выбранную ступень усиления RTL-SDR. + +Для каждого параметра в отчёте и метаданных хранятся отдельные поля +`requested` и `actual`. Если API не предоставляет независимого чтения +обратно, поле `actual` имеет значение `N/A`, а не копию запрошенного +значения. Передача не начинается, если частота дискретизации или несущая +не прочитаны обратно либо отличаются от запрошенных больше, чем допускает +API устройства и заранее установленная проверка. + +## 6. Двухтоновая калибровка + +Перед BPSK Pluto+ передаёт два комплексных тона, симметричных относительно +несущей: + +`-Fcal` и `+Fcal`, где `Fcal = 50 кГц`. + +50 кГц выбраны до измерения: оба тона находятся внутри полосы 200 кГц, +не сливаются около нуля и дают разнос 100 кГц для оценки масштаба частот. +RTL-SDR измеряет положения двух пиков относительно своей настроенной +несущей: `f_low` и `f_high`. + +### 6.1. Рассогласование несущей + +`carrier_offset_hz = (f_low + f_high) / 2` + +`carrier_offset_ppm = carrier_offset_hz / 435000000 * 1e6` + +Положительный результат означает, что принятый спектр сдвинут вверх по +частоте. Для компенсации сырые комплексные отсчёты умножаются на: + +`exp(-j * 2*pi*carrier_offset_hz*n/Fs_rx_actual)` + +Знак поправки проверяется независимо: после компенсации среднее положение +двух тонов должно стать ближе к нулю. Рассогласование несущей не должно +подменяться оценкой ошибки частоты дискретизации. + +### 6.2. Относительная ошибка частот дискретизации + +`clock_scale = (f_high - f_low) / (2*Fcal)` + +`sample_clock_error_ppm = (clock_scale - 1) * 1e6` + +Здесь положительная ошибка означает, что масштаб передающего такта больше +масштаба приёмного. Для перехода принятых отсчётов на временную сетку +передатчика ожидаемая длина после передискретизации равна: + +`len_corrected = round(len_received * clock_scale)` + +Направление этой операции не принимается на веру. Оно подтверждается +машинными тестами из раздела 8. После передискретизации накопленная ошибка +символьного времени должна уменьшаться. + +### 6.3. Начальная фаза выбора отсчёта символа + +После грубой CFO-коррекции и передискретизации приёмник отдельно ищет +начальную фазу выбора отсчёта в диапазоне от 0 до 119 номинальных +отсчётов на символ. Это не CFO и не ошибка такта. В результатах отдельно +хранятся: + +- `carrier_offset_hz` и `carrier_offset_ppm`; +- `clock_scale` и `sample_clock_error_ppm`; +- `symbol_sample_phase` и остаточная ошибка символьного времени. + +### 6.4. Формальное определение шумового фона + +Шумовой фон и превышение пиков вычисляются одинаковым детерминированным +алгоритмом: + +1. Из комплексных отсчётов вычитается среднее. +2. Оценивается спектральная плотность мощности методом Welch с окном Hann, + `nfft = nperseg = 65536` и перекрытием 50 процентов. +3. Разрешение одного БПФ вычисляется по прочитанной обратно частоте как + `delta_f = Fs_rx_actual / 65536`. При 2,4 Мвыб/с оно равно + `36,62109375 Гц`. +4. Анализируется полоса от -100 до +100 кГц относительно настроенной + несущей. +5. Нижний пик ищется в диапазоне от -80 до -20 кГц, верхний от +20 до + +80 кГц. Эти непересекающиеся окна заданы до измерения. +6. После нахождения кандидатов ожидаемых тонов из оценки шума исключаются + защитные зоны шириной +/-2 кГц вокруг измеренного положения каждого + кандидата. Каждая полузона занимает около 55 разрешающих элементов БПФ + при 2,4 Мвыб/с, поэтому утечка основного лепестка и ближайших боковых + лепестков тона не попадает в медиану шума. +7. Шумовая мощность `P_noise` равна медиане линейных значений мощности + всех оставшихся спектральных элементов анализируемой полосы. +8. Для каждого тона вычисляется + `peak_excess_db = 10*log10(P_peak/P_noise)`. + +Оба пика должны превышать шумовой фон не менее чем на 10 дБ. Порог 10 дБ, +полоса анализа, окна поиска и защитные зоны зафиксированы до аппаратного +измерения. Если нет двух конечных пиков, недостаточно элементов шума или +хотя бы один пик не проходит порог, оценка этого захвата недействительна. + +### 6.5. Повторы и представление недействительных оценок + +Без изменения схемы выполняются три отдельных калибровочных захвата. Для +каждого захвата отдельно сохраняются: + +- `f_low_hz` и `f_high_hz`; +- `low_peak_excess_db` и `high_peak_excess_db`; +- `carrier_offset_hz` и `carrier_offset_ppm`; +- `clock_scale` и `sample_clock_error_ppm`; +- остаточная CFO после коррекции; +- остаточная ошибка разноса тонов в герцах и ppm; +- признак достоверности и причина отказа. + +По трём захватам вычисляются среднее, минимум, максимум и стандартное +отклонение каждой конечной оценки, а также число недействительных +захватов. Стандартное отклонение считается по совокупности всех конечных +оценок с `ddof=0`. Недействительное значение записывается как `NaN` в CSV +и как `N/A` в текстовом отчёте. Оно никогда не заменяется нулём. Если +конечных значений нет, все агрегаты этой величины равны `N/A`. + +## 7. Диагностическая последовательность передачи + +Аппаратная часть выполняется по ступеням. Переход к следующей ступени +разрешён только после сохранения результатов предыдущей: + +1. двухтоновые калибровочные захваты; +2. известная короткая BPSK-последовательность; +3. символьная синхронизация и измерение BER известной последовательности; +4. один пакет с CRC; +5. десять пакетированных фрагментов JPEG; +6. сборка JPEG; +7. проверка размера и SHA-256. + +Известная последовательность задаётся до измерения как PRBS11 длиной +2047 бит с полиномом `x^11 + x^9 + 1` и фиксированным ненулевым начальным +состоянием, записанным в метаданные. BER вычисляется прямым сравнением +каждого принятого бита с ожидаемым до упаковки в пакет и до проверки CRC. +CRC не может использоваться вместо измерения BER. + +### 7.1. Контрольный опыт PRBS11 после проверки непрерывного async-приёма + +Короткая запись S и длинная запись L являются двумя отдельными аппаратными +захватами. Каждый захват выполняется одним непрерывным +`rtlsdr_read_async()` и содержит в этой последовательности: входной запас, +двухтоновую калибровку, фиксированный защитный участок 0,10 с, один +начальный маркер, непрерывную PRBS11 и выходной запас. Полезная PRBS11 +длится 0,25 с в S и не менее 2,0 с в L. Внутри полезной PRBS11 нет +периодической повторной синхронизации. + +Оценки CFO и масштаба частоты дискретизации для обработки PRBS11 берутся +только из двух тонов в той же самой IQ-записи. Оценку из предыдущего +двухтонового опыта или из другой записи применять запрещено. Если тоны в +текущем захвате недостоверны, PRBS11 этой записи не интерпретируется. + +Для обоих захватов заранее фиксируются 435 МГц, номинальные 2 400 000 +отсчётов/с TX и RX, 20 000 символов/с, TX -30 дБ, запрос усиления RTL-SDR +19,7 дБ и `Fcal = 50 кГц`. Асинхронный приём использует callback 262144 +байта, 15 буферов и отбрасывает первый callback как разогревочный без +перезапуска сеанса. Автоматический подбор параметров запрещён. + +Одна и та же IQ-запись обрабатывается режимами A, B и C. В B и C грубая +CFO берётся из тонов этой записи; в C дополнительно применяется +`clock_scale` этих же тонов. Разрешённая защищённая тонкая CFO-коррекция не +заменяет двухтоновое измерение. Для длинной записи SRO дополнительно +оценивается независимо по накопленному дрейфу символьной фазы. + +## 8. Обязательные программные проверки до аппаратного опыта + +На синтетических комплексных данных проверяются четыре ошибки частоты +дискретизации: +20, -20, +100 и -100 ppm. Для каждого случая тест обязан +проверить: + +- знак оценки `sample_clock_error_ppm`; +- величину оценки в заранее заданном допуске; +- направление изменения длины при передискретизации; +- уменьшение остаточной накопленной ошибки символьного времени после + поправки. + +Для положительной ошибки длина должна изменяться в направлении, +предписанном формулой `round(N * clock_scale)`, для отрицательной в +противоположном. Тест должен упасть, если вместо `clock_scale` применена +обратная величина. Допуск оценки и критерий уменьшения остаточной ошибки +задаются в тесте до аппаратного измерения и печатаются в его результате. + +Дополнительно синтетические тесты проверяют: + +- знак грубой CFO-коррекции для положительной и отрицательной CFO; +- уменьшение остаточной CFO после поправки; +- отдельное восстановление начальной фазы выбора отсчёта символа; +- явный `N/A`, если два тона не найдены достоверно; +- формулу шумового фона и порог 10 дБ; +- BER известной последовательности до CRC; +- отсутствие ложного успеха при неполном пакете или изображении; +- формирование кода возврата только по полному набору критериев. + +До прохождения этих проверок передатчик включать нельзя. + +## 9. Порядок аппаратного эксперимента + +1. Открыть Pluto+ успешным `iio.Context` и записать URI и идентификаторы. +2. Открыть RTL-SDR через `pyrtlsdr` и записать имя устройства и тюнера. +3. Установить и прочитать обратно параметры из раздела 5. +4. Получить отдельное подтверждение пользователя, что кабельная схема из + раздела 4 собрана и AT30S находится в разрыве. +5. Запустить один непрерывный приём RTL-SDR за 0,5 с до передачи Pluto+. +6. Выполнить три отдельных двухтоновых калибровочных захвата. +7. Сохранить оценки каждого захвата, агрегаты и не менее одного + эталонного сырого IQ-захвата. +8. Передать известную BPSK-последовательность и измерить BER. +9. Принять один пакет и проверить CRC. +10. Передать десять фрагментов JPEG и проверить критерии раздела 11. + +Если любая ступень не проходит, переход к следующей запрещён. Исходные +данные и отрицательный результат сохраняются без подгонки. Изменение +одного параметра для следующего запуска оформляется как отдельное заранее +объяснённое действие. + +Весь опыт принимается RTL-SDR как одна непрерывная временная +последовательность. Продолжительность приёма вычисляется из фактически +сформированной передаваемой последовательности: + +`T_rx = 0,5 с + N_tx / Fs_tx_actual + 0,5 с`. + +Здесь `N_tx` является фактическим числом комплексных отсчётов всей +сформированной TX-последовательности. Общая длительность опыта не задаётся +произвольной константой. Поток разрешено читать блоками, например около +262144 отсчётов, но блоки объединяются без пропусков и перестановок в одну +последовательность. Размер блока чтения является параметром реализации, а +не физическим защитным интервалом. Десять отдельных запусков приёмника для +десяти кадров запрещены, потому что они уничтожили бы наблюдаемое +накопление ошибки тактов. + +Между BPSK-кадрами не добавляются новые интервалы. Используется структура +Lab042: `GUARD_SYMBOL_COUNT` импортируется из её текущей конфигурации без +копии числовой константы в Lab043, а формирование кадра повторяет ту же +последовательность переднего интервала, символов, заднего интервала и RRC. +Готовая функция Lab042 жёстко связана с 128 отсчётами на символ, поэтому +при 120 отсчётах Lab043 переиспользуются её параметры, а не скрытая копия +её глобальной частоты дискретизации. Если в дальнейшем отдельный интервал +исчезнет из Lab042, Lab043 не создаёт новый без отдельного диагностического +решения. + +## 10. Три режима обработки одной записи + +Одна и та же сохранённая BPSK/IQ-запись обрабатывается тремя режимами: + +- режим A: без грубой CFO-коррекции и без коррекции частоты + дискретизации; +- режим B: только грубая CFO-коррекция из двухтоновой калибровки; +- режим C: грубая CFO-коррекция, коррекция частоты дискретизации и + защищённая тонкая CFO-коррекция. + +Во всех режимах начальная фаза выбора отсчёта символа оценивается и +записывается отдельно. Для каждого режима сохраняются BER известной +последовательности, остаточная CFO, остаточная временная ошибка, +корреляция маркера, число найденных кадров, число разобранных заголовков и +CRC пакетов. + +В режиме C защищённая функция `should_apply_cfo_correction` из +`protocol/bpsk_radio.py` остаётся тонкой коррекцией после грубого +измерения. Она не используется как единственный способ поиска сдвига в +десятки килогерц, поскольку символьная оценка неоднозначна за пределами +половины символьной скорости. + +Провал режимов A или B не считается провалом лабораторной. Они являются +контрольными измерениями. Критерии передачи применяются к режиму C. + +## 11. Критерии приёмки + +### 11.1. Радиокритерии режима C + +- найдены все 10 кадров; +- разобраны все 10 пакетов; +- CRC32 сошлась у всех 10 пакетов. + +### 11.2. Прикладные критерии режима C + +- восстановлены все 10 уникальных фрагментов; +- JPEG собран без заполнения отсутствующих данных; +- размер принятого JPEG совпадает с размером переданного; +- SHA-256 принятого JPEG совпадает с SHA-256 переданного. + +### 11.3. Общий результат + +Код возврата 0 разрешён только при одновременном выполнении: + +- всех радиокритериев; +- всех прикладных критериев; +- всех обязательных программных проверок раздела 8. + +Для подтверждения воспроизводимости заранее выбранная рабочая точка +запускается три раза подряд без изменения схемы и параметров. Каждый +повтор сохраняется отдельно. Три повтора являются проверкой +воспроизводимости, а не статистической оценкой вероятности отказа. + +Отрицательный результат допустим и сохраняется полностью. Запрещено +автоматически поднимать усиление, менять частоту, пороги или алгоритм +синхронизации после неудачного прогона. + +## 12. Измерения и артефакты + +### 12.1. Измерения передачи + +Для каждого запуска и каждого режима обработки сохраняются: + +- запрошенные и фактические усиления, несущие и частоты дискретизации; +- число принятых отсчётов; +- отдельные оценки CFO, ошибки частоты дискретизации и начальной + символьной фазы; +- BER известной последовательности до CRC; +- число найденных кадров и разобранных заголовков; +- CRC каждого пакета; +- решение защищённой тонкой CFO-коррекции и его диагностические величины; +- корреляция маркера и корректно определённая EVM; +- число и номера восстановленных фрагментов; +- размер и SHA-256 переданного и принятого JPEG; +- длительность сигнала и полезная скорость. + +Каждое число вычисляется из текущего запуска. Перенесённые константы +снабжаются ссылкой на лабораторную-источник. Невычислимая величина +записывается как `NaN` или `N/A`, но не как правдоподобный ноль. + +### 12.2. План программных артефактов + +Предполагаемые файлы после отдельного согласования реализации: + +- `experiments/lab043_pluto_to_rtlsdr.py`; +- `tests/test_lab043_calibration.py`; +- `data/processed/lab043/lab043_calibration.csv`; +- `data/processed/lab043/lab043_frames.csv`; +- `data/processed/lab043/lab043_summary.csv`; +- `data/processed/lab043/lab043_report.txt`; +- графики калибровочного спектра, остаточной ошибки и доставки + фрагментов. + +### 12.3. Эталонный сырой IQ-захват + +Не менее одного полного эталонного захвата до CFO- и SRO-коррекций +сохраняется в `data/raw/lab043/` и не добавляется в Git. + +Формат отсчётов: NumPy `.npy`, одномерный массив `complex64` в порядке +приёма. Нормирование исходных 8-битных I/Q отсчётов и порядок I/Q должны +быть однозначно описаны в метаданных. Рядом сохраняется UTF-8 JSON с тем +же базовым именем и полями: + +- локальная метка времени и UTC; +- идентификаторы устройств и модель тюнера; +- версии Python, `pyrtlsdr`, `pyrtlsdrlib`, `librtlsdr`, `pyadi-iio` и + `pylibiio`; +- запрошенные и фактические частоты дискретизации, несущие и усиления; +- число отсчётов, длительность, dtype и порядок байтов; +- идентификатор и SHA-256 переданной формы сигнала; +- результаты калибровки, применённые к этой записи; +- SHA-256 файла `.npy`. + +IQ-файл содержит всю единую последовательность приёма. Ожидаемое число +отсчётов вычисляется как `ceil(T_rx * Fs_rx_actual)`. Чтение выполняется +последовательными блоками до достижения этого числа, последний блок +обрезается только после получения требуемого количества отсчётов. + +### 12.4. Разделение общего и лабораторного кода + +Общие математические примитивы обнаружения двух известных тонов, оценки +грубой CFO и относительной ошибки частот дискретизации, грубой частотной +коррекции и передискретизации размещаются в `protocol/bpsk_radio.py` либо +в другом подходящем существующем модуле `protocol/`. Каждый такой перенос +сопровождается быстрым синтетическим тестом. + +В `experiments/lab043_pluto_to_rtlsdr.py` остаются сценарий лабораторной, +конкретное `Fcal`, три калибровочных захвата, методика оценки шума Lab043, +режимы A/B/C, PRBS11, контрольный пакет, JPEG, отчётные артефакты, +критерии приёмки и аппаратная последовательность. Общий код не зависит от +конкретного JPEG и номера лабораторной. + +## 13. Точка остановки перед аппаратным опытом + +Программная реализация и синтетические проверки разрешены. Перед +аппаратным опытом необходимо отдельно: + +1. пройти все обязательные программные проверки раздела 8; +2. повторно проверить доступность устройств и чтение параметров обратно; +3. получить подтверждение пользователя о схеме и отдельное разрешение на + включение TX. + +Способ Python-доступа к RTL-SDR уже выбран и не является вопросом точки +остановки. + +До отдельного разрешения передатчик не включается, усиление Pluto+ не +изменяется, аппаратная калибровка и передача тонов, BPSK или JPEG не +выполняются. Операции `git add`, `commit`, `push`, `pull` и `fetch` также +не выполняются. diff --git a/experiments/lab043_pluto_to_rtlsdr.py b/experiments/lab043_pluto_to_rtlsdr.py new file mode 100644 index 0000000..6e7d0af --- /dev/null +++ b/experiments/lab043_pluto_to_rtlsdr.py @@ -0,0 +1,1832 @@ +"""Lab043: программная часть тракта Pluto+ -> независимый RTL-SDR. + +Текущий этап содержит только детерминированные синтетические проверки и +чистые функции подготовки обработки. Модуль не импортирует аппаратные +библиотеки, не открывает SDR и не включает передатчик. + +Аппаратная постановка описана в docs/lab043_pluto_to_rtlsdr_spec.md. +""" + +from __future__ import annotations + +import csv +import hashlib +import json +import math +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Iterable + +import numpy as np +from matplotlib.figure import Figure +from scipy.signal import fftconvolve, welch + +from experiments.lab042_pluto_image_loopback import ( + GUARD_SYMBOL_COUNT as LAB042_GUARD_SYMBOL_COUNT, + RRC_ROLLOFF as LAB042_RRC_ROLLOFF, + RRC_SPAN_SYMBOLS as LAB042_RRC_SPAN_SYMBOLS, + TX_AMPLITUDE as LAB042_TX_AMPLITUDE, +) +from protocol import bpsk_radio as radio +from protocol.image_fragments import MissingFragmentsError, reassemble_image +from protocol.packet import MESSAGE_TYPE_TEXT, build_packet, parse_packet + + +SYMBOL_RATE = 20_000 +SAMPLE_RATE_HZ = 2_400_000 +SAMPLES_PER_SYMBOL = SAMPLE_RATE_HZ // SYMBOL_RATE +CARRIER_HZ = 435_000_000 +F_CAL_HZ = 50_000.0 + +CALIBRATION_NFFT = 65_536 +CALIBRATION_ANALYSIS_HALF_BAND_HZ = 100_000.0 +CALIBRATION_SEARCH_HALF_WIDTH_HZ = 30_000.0 +CALIBRATION_GUARD_HALF_WIDTH_HZ = 2_000.0 +MINIMUM_TONE_EXCESS_DB = 10.0 +MAXIMUM_ABS_SAMPLE_CLOCK_ERROR_PPM = 1_000.0 + +RX_LEADING_MARGIN_SECONDS = 0.5 +RX_TRAILING_MARGIN_SECONDS = 0.5 +DEFAULT_READ_BLOCK_SAMPLES = 262_144 + +RTL_ASYNC_BUFFER_BYTES = 262_144 +RTL_ASYNC_BUFFER_COUNT = 15 +RTL_ASYNC_WARMUP_CALLBACK_COUNT = 1 +RTL_BYTES_PER_COMPLEX_SAMPLE = 2 + +PRBS11_LENGTH = 2_047 +PRBS11_INITIAL_STATE = 0x7FF +PRBS11_POLYNOMIAL = "x^11 + x^9 + 1" +PRBS_SHORT_DURATION_SECONDS = 0.25 +PRBS_LONG_DURATION_SECONDS = 2.0 +PILOT_SYMBOL_COUNT = 1_280 +PILOT_INITIAL_STATE = 0x7F +PILOT_POLYNOMIAL = "x^7 + x^6 + 1" +CALIBRATION_TONE_DURATION_SECONDS = 0.25 +CALIBRATION_TO_BPSK_GUARD_SECONDS = 0.10 +PRBS_MARKER_MINIMUM_CORRELATION = 0.65 +EXPECTED_FRAME_COUNT = 10 + + +@dataclass(frozen=True) +class CalibrationResult: + """Единый явный контракт полной двухтоновой калибровки.""" + + valid: bool + invalid_reason: str + f_low_hz: float + f_high_hz: float + carrier_offset_hz: float + carrier_offset_ppm: float + clock_scale: float + sample_clock_error_ppm: float + residual_cfo_hz: float + phase_fit_rmse_rad: float + peak_margin_low_db: float + peak_margin_high_db: float + noise_power: float + + @property + def low_peak_excess_db(self) -> float: + """Совместимое имя для старых отчётных потребителей.""" + + return self.peak_margin_low_db + + @property + def high_peak_excess_db(self) -> float: + """Совместимое имя для старых отчётных потребителей.""" + + return self.peak_margin_high_db + + @property + def failure_reason(self) -> str: + """Совместимое имя причины отказа.""" + + return self.invalid_reason + + +# Старое имя оставлено только как импортная совместимость. Объект и контракт +# один: новые функции и отчёты используют CalibrationResult. +CalibrationEstimate = CalibrationResult + + +@dataclass(frozen=True) +class ModeResult: + """Результат обработки одной записи одним из режимов A/B/C.""" + + mode: str + bit_error_rate: float + bit_errors: int + bit_count: int + marker_correlation: float + symbol_sample_phase: int + fine_cfo_applied: bool + estimated_fine_cfo_hz: float + + +@dataclass(frozen=True) +class PrbsTransmissionPlan: + """Заранее фиксированный состав одной калибровочно-PRBS11 посылки.""" + + label: str + useful_duration_seconds: float + payload_bits: np.ndarray + marker_bits: np.ndarray + pilot_bits: np.ndarray + tx_samples: np.ndarray + calibration_sample_count: int + fixed_guard_sample_count: int + bpsk_start_sample: int + + +@dataclass(frozen=True) +class PrbsModeMetrics: + """Метрики одного режима обработки той же непрерывной IQ-записи.""" + + mode: str + detected: bool + failure_reason: str + transmitted_bit_count: int + matched_bit_count: int + bit_errors: int + bit_error_rate: float + quarter_bit_error_rates: tuple[float, float, float, float] + evm_percent: float + marker_correlation: float + symbol_phase_start_samples: float + symbol_phase_end_samples: float + accumulated_timing_drift_samples: float + accumulated_timing_drift_symbols: float + residual_cfo_hz: float + maximum_correct_run_bits: int + fine_cfo_applied: bool + estimated_fine_cfo_hz: float + timing_sro_ppm: float + coarse_cfo_applied: bool + pilot_cfo_applied: bool + pilot_estimate_valid: bool + estimated_pilot_cfo_hz: float + residual_pilot_cfo_hz: float + pilot_phase_fit_rmse_rad: float + pilot_mean_block_coherence: float + legacy_prbs_adjacent_phase_hz: float + bit_pattern_diagnostics: dict + + +@dataclass(frozen=True) +class AcceptanceResult: + """Раздельные критерии радиотракта и прикладного уровня.""" + + frames_found: int + packets_parsed: int + packets_crc_valid: int + fragments_recovered: int + image_reassembled: bool + image_size_matches: bool + image_sha256_matches: bool + software_tests_passed: bool + + @property + def radio_passed(self) -> bool: + return ( + self.frames_found == EXPECTED_FRAME_COUNT + and self.packets_parsed == EXPECTED_FRAME_COUNT + and self.packets_crc_valid == EXPECTED_FRAME_COUNT + ) + + @property + def application_passed(self) -> bool: + return ( + self.fragments_recovered == EXPECTED_FRAME_COUNT + and self.image_reassembled + and self.image_size_matches + and self.image_sha256_matches + ) + + @property + def passed(self) -> bool: + return self.radio_passed and self.application_passed and self.software_tests_passed + + +@dataclass(frozen=True) +class FunctionalTestResult: + """Результат одной встроенной программной проверки.""" + + name: str + passed: bool + detail: str + + +@dataclass(frozen=True) +class AsyncCallbackBoundary: + """Граница одного callback непрерывного RTL-SDR-захвата.""" + + callback_index: int + source_sample_count: int + accepted_start_sample: int + accepted_end_sample: int + accepted_sample_count: int + warmup_discarded: bool + trimmed_at_target: bool + + +class ContinuousAsyncIqCollector: + """Собрать одну запись из последовательно пронумерованных callback-блоков. + + Коллектор ничего не знает об аппаратуре и проверяет программную часть + непрерывного захвата: повтор, пропуск или перестановка номера callback + приводят к ошибке. Первый заранее заданный callback можно отбросить как + разогревочный, не перезапуская асинхронный сеанс. + """ + + def __init__( + self, + expected_sample_count: int, + warmup_callback_count: int = RTL_ASYNC_WARMUP_CALLBACK_COUNT, + ) -> None: + if expected_sample_count <= 0: + raise ValueError("Ожидаемое число отсчётов должно быть положительным") + if warmup_callback_count < 0: + raise ValueError("Число разогревочных callback не может быть отрицательным") + self.expected_sample_count = int(expected_sample_count) + self.warmup_callback_count = int(warmup_callback_count) + self._expected_callback_index = 0 + self._accepted_sample_count = 0 + self._blocks: list[np.ndarray] = [] + self._boundaries: list[AsyncCallbackBoundary] = [] + + @property + def complete(self) -> bool: + return self._accepted_sample_count == self.expected_sample_count + + @property + def accepted_sample_count(self) -> int: + return self._accepted_sample_count + + @property + def callback_boundaries(self) -> tuple[AsyncCallbackBoundary, ...]: + return tuple(self._boundaries) + + def add_callback(self, callback_index: int, samples: np.ndarray) -> bool: + """Добавить ровно следующий callback и вернуть признак завершения.""" + + if self.complete: + raise RuntimeError("Захват уже набрал заданное число отсчётов") + if callback_index != self._expected_callback_index: + if callback_index < self._expected_callback_index: + raise ValueError("Callback продублирован или переставлен назад") + raise ValueError("В последовательности callback обнаружен пропуск") + + block = np.asarray(samples, dtype=np.complex64) + if block.ndim != 1 or block.size == 0: + raise ValueError("Callback должен содержать непустой одномерный IQ-блок") + + warmup = callback_index < self.warmup_callback_count + start = self._accepted_sample_count + if warmup: + accepted = 0 + else: + accepted = min(block.size, self.expected_sample_count - start) + self._blocks.append(block[:accepted].copy()) + self._accepted_sample_count += accepted + end = self._accepted_sample_count + self._boundaries.append( + AsyncCallbackBoundary( + callback_index=callback_index, + source_sample_count=int(block.size), + accepted_start_sample=int(start), + accepted_end_sample=int(end), + accepted_sample_count=int(accepted), + warmup_discarded=warmup, + trimmed_at_target=bool(not warmup and accepted < block.size), + ) + ) + self._expected_callback_index += 1 + return self.complete + + def finalize(self) -> np.ndarray: + """Вернуть запись только после получения точного целевого размера.""" + + if not self.complete: + raise RuntimeError("Асинхронный захват короче заданной длины") + combined = np.concatenate(self._blocks).astype(np.complex64, copy=False) + if len(combined) != self.expected_sample_count: + raise RuntimeError("Сумма принятых callback не совпала с целевой длиной") + return combined + + +def capture_continuous_rtlsdr_async( + sdr, + expected_sample_count: int, + capture_ready_event=None, + buffer_bytes: int = RTL_ASYNC_BUFFER_BYTES, + buffer_count: int = RTL_ASYNC_BUFFER_COUNT, + warmup_callback_count: int = RTL_ASYNC_WARMUP_CALLBACK_COUNT, +) -> tuple[np.ndarray, tuple[AsyncCallbackBoundary, ...], dict]: + """Принять RTL-SDR IQ одним непрерывным ``rtlsdr_read_async``. + + Используется существующая ctypes-привязка pyrtlsdr. На установленной + Windows-сборке высокоуровневый ``read_bytes_async`` после длительной + штатной отмены пытается повторно закрыть уже недоступный USB-дескриптор. + Поэтому здесь напрямую вызывается уже объявленная функция librtlsdr, без + собственной DLL-обёртки. Вызывающий код не должен повторно использовать + объект SDR, если ``read_async_result`` отрицателен. + """ + + if buffer_bytes <= 0 or buffer_bytes % 512 != 0: + raise ValueError("Размер async-буфера должен быть положительным и кратным 512 байтам") + if buffer_bytes % 16_384 != 0: + raise ValueError("Размер async-буфера должен быть кратным 16384 байтам") + if buffer_bytes % RTL_BYTES_PER_COMPLEX_SAMPLE != 0: + raise ValueError("Async-буфер должен содержать целое число комплексных отсчётов") + if buffer_count <= 0: + raise ValueError("Число async-буферов должно быть положительным") + + from ctypes import py_object + import threading + import time + + from rtlsdr import librtlsdr + from rtlsdr.librtlsdr import rtlsdr_read_async_cb_t + + collector = ContinuousAsyncIqCollector( + expected_sample_count, + warmup_callback_count=warmup_callback_count, + ) + done = threading.Event() + cancel_errors: list[str] = [] + callback_index = 0 + started = time.perf_counter() + + def receive_bytes(raw_bytes, _context) -> None: + nonlocal callback_index + iq = sdr.packed_bytes_to_iq(raw_bytes) + complete = collector.add_callback(callback_index, iq) + callback_index += 1 + if callback_index == warmup_callback_count and capture_ready_event is not None: + capture_ready_event.set() + if complete: + done.set() + + def cancel_when_complete() -> None: + done.wait() + try: + sdr.cancel_read_async() + except Exception as error: # pragma: no cover - зависит от USB-библиотеки + cancel_errors.append(f"{type(error).__name__}: {error}") + + sdr.DEFAULT_ASYNC_BUF_NUMBER = int(buffer_count) + sdr._callback_bytes = receive_bytes + sdr.read_async_canceling = False + callback = rtlsdr_read_async_cb_t(sdr._bytes_converter_callback) + cancel_thread = threading.Thread(target=cancel_when_complete, daemon=True) + cancel_thread.start() + result = librtlsdr.rtlsdr_read_async( + sdr.dev_p, + callback, + py_object(sdr), + int(buffer_count), + int(buffer_bytes), + ) + elapsed = time.perf_counter() - started + cancel_thread.join(timeout=5.0) + + if not collector.complete: + raise RuntimeError( + f"Асинхронное чтение завершилось до целевого размера, код {result}" + ) + if result not in (0, -5): + raise RuntimeError(f"Асинхронное чтение завершилось с кодом {result}") + samples = collector.finalize() + diagnostics = { + "read_async_result": int(result), + "elapsed_seconds": float(elapsed), + "callback_count": int(callback_index), + "callback_boundary_count": max(0, int(callback_index) - 1), + "buffer_bytes": int(buffer_bytes), + "buffer_count": int(buffer_count), + "warmup_callback_count": int(warmup_callback_count), + "cancel_errors": tuple(cancel_errors), + } + return samples, collector.callback_boundaries, diagnostics + + +def receive_duration_seconds(tx_sample_count: int, actual_tx_sample_rate_hz: float) -> float: + """Рассчитать длительность единого приёма из фактического TX-сигнала.""" + + if tx_sample_count < 0: + raise ValueError("Число TX-отсчётов не может быть отрицательным") + if not np.isfinite(actual_tx_sample_rate_hz) or actual_tx_sample_rate_hz <= 0.0: + raise ValueError("Фактическая частота TX должна быть положительной") + return ( + RX_LEADING_MARGIN_SECONDS + + tx_sample_count / actual_tx_sample_rate_hz + + RX_TRAILING_MARGIN_SECONDS + ) + + +def receive_sample_count( + tx_sample_count: int, + actual_tx_sample_rate_hz: float, + actual_rx_sample_rate_hz: float, +) -> int: + """Вернуть число RX-отсчётов для единого непрерывного захвата.""" + + if not np.isfinite(actual_rx_sample_rate_hz) or actual_rx_sample_rate_hz <= 0.0: + raise ValueError("Фактическая частота RX должна быть положительной") + duration = receive_duration_seconds(tx_sample_count, actual_tx_sample_rate_hz) + return math.ceil(duration * actual_rx_sample_rate_hz) + + +def combine_continuous_blocks( + blocks: Iterable[np.ndarray], + expected_sample_count: int, +) -> np.ndarray: + """Объединить последовательные блоки чтения в одну временную запись.""" + + if expected_sample_count < 0: + raise ValueError("Ожидаемое число отсчётов не может быть отрицательным") + arrays = [np.asarray(block, dtype=np.complex64) for block in blocks] + if any(array.ndim != 1 for array in arrays): + raise ValueError("Каждый блок должен быть одномерным") + combined = np.concatenate(arrays) if arrays else np.empty(0, dtype=np.complex64) + if len(combined) < expected_sample_count: + raise ValueError("Непрерывный захват короче рассчитанной длительности") + return combined[:expected_sample_count] + + +def shape_frame_like_lab042( + symbols: np.ndarray, + samples_per_symbol: int = SAMPLES_PER_SYMBOL, +) -> np.ndarray: + """Сформировать кадр с существующими защитными интервалами Lab042.""" + + symbols = np.asarray(symbols, dtype=np.complex128) + if symbols.ndim != 1 or symbols.size == 0: + raise ValueError("Нужна непустая одномерная последовательность символов") + if samples_per_symbol <= 0: + raise ValueError("Число отсчётов на символ должно быть положительным") + + taps = radio.root_raised_cosine_taps( + LAB042_RRC_ROLLOFF, + samples_per_symbol, + LAB042_RRC_SPAN_SYMBOLS, + ) + upsampled = np.zeros(len(symbols) * samples_per_symbol, dtype=np.complex128) + upsampled[::samples_per_symbol] = symbols + guard = np.zeros( + LAB042_GUARD_SYMBOL_COUNT * samples_per_symbol, + dtype=np.complex128, + ) + return ( + fftconvolve(np.concatenate((guard, upsampled, guard)), taps, mode="full") + * LAB042_TX_AMPLITUDE + ) + + +def concatenate_frames_without_new_gaps(frame_waveforms: Iterable[np.ndarray]) -> np.ndarray: + """Соединить готовые кадры без добавочного интервала Lab043.""" + + frames = [np.asarray(frame, dtype=np.complex128) for frame in frame_waveforms] + if any(frame.ndim != 1 for frame in frames): + raise ValueError("Каждый кадр должен быть одномерным") + return np.concatenate(frames) if frames else np.empty(0, dtype=np.complex128) + + +def prbs11( + length: int = PRBS11_LENGTH, + initial_state: int = PRBS11_INITIAL_STATE, +) -> np.ndarray: + """Сформировать детерминированную PRBS11 с полиномом x^11+x^9+1.""" + + if length <= 0: + raise ValueError("Длина PRBS11 должна быть положительной") + if not 1 <= initial_state <= 0x7FF: + raise ValueError("Начальное состояние PRBS11 должно быть ненулевым 11-битным") + + state = initial_state + bits = np.empty(length, dtype=np.uint8) + for index in range(length): + bits[index] = (state >> 10) & 1 + feedback = ((state >> 10) ^ (state >> 8)) & 1 + state = ((state << 1) & 0x7FF) | feedback + return bits + + +def pilot_prbs7( + length: int = PILOT_SYMBOL_COUNT, + initial_state: int = PILOT_INITIAL_STATE, +) -> np.ndarray: + """Сформировать отдельный детерминированный пилот x^7+x^6+1. + + Последовательность и начальное состояние фиксированы до аппаратного + опыта и не зависят от полезной PRBS11 либо сохранённого IQ. + """ + + if length <= 0: + raise ValueError("Длина пилота должна быть положительной") + if not 1 <= initial_state <= 0x7F: + raise ValueError("Начальное состояние пилота должно быть ненулевым 7-битным") + + state = initial_state + bits = np.empty(length, dtype=np.uint8) + for index in range(length): + bits[index] = (state >> 6) & 1 + feedback = ((state >> 6) ^ (state >> 5)) & 1 + state = ((state << 1) & 0x7F) | feedback + return bits + + +def build_prbs11_transmission_plan( + label: str, + useful_duration_seconds: float, + sample_rate_hz: int = SAMPLE_RATE_HZ, + samples_per_symbol: int = SAMPLES_PER_SYMBOL, +) -> PrbsTransmissionPlan: + """Собрать один непрерывный TX: два тона, фиксированный guard и PRBS11. + + Перед PRBS используется ровно один известный маркер для начальной + синхронизации. Внутри полезной PRBS периодических маркеров и повторных + запусков синхронизации нет. + """ + + if label not in {"S", "L"}: + raise ValueError("Метка опыта должна быть S или L") + if useful_duration_seconds <= 0.0: + raise ValueError("Полезная длительность PRBS11 должна быть положительной") + if sample_rate_hz != samples_per_symbol * SYMBOL_RATE: + raise ValueError("Fs должна быть целым числом отсчётов на символ") + + payload_bit_count = round(useful_duration_seconds * SYMBOL_RATE) + if not math.isclose( + payload_bit_count / SYMBOL_RATE, + useful_duration_seconds, + rel_tol=0.0, + abs_tol=1e-12, + ): + raise ValueError("Полезная длительность должна содержать целое число символов") + + payload_bits = prbs11(payload_bit_count, PRBS11_INITIAL_STATE) + marker_bits, marker_symbols = radio.build_frame_marker() + pilot_bits = ( + pilot_prbs7(PILOT_SYMBOL_COUNT, PILOT_INITIAL_STATE) + if label == "S" + else np.empty(0, dtype=np.uint8) + ) + bpsk_symbols = np.concatenate( + ( + marker_symbols, + radio.bpsk_modulate(pilot_bits), + radio.bpsk_modulate(payload_bits), + ) + ) + bpsk_samples = shape_frame_like_lab042(bpsk_symbols, samples_per_symbol) + + calibration_sample_count = round(CALIBRATION_TONE_DURATION_SECONDS * sample_rate_hz) + indexes = np.arange(calibration_sample_count, dtype=np.float64) + calibration = 0.35 * ( + np.exp(-1j * 2.0 * np.pi * F_CAL_HZ * indexes / sample_rate_hz) + + np.exp(1j * 2.0 * np.pi * F_CAL_HZ * indexes / sample_rate_hz) + ) + fixed_guard_sample_count = round(CALIBRATION_TO_BPSK_GUARD_SECONDS * sample_rate_hz) + fixed_guard = np.zeros(fixed_guard_sample_count, dtype=np.complex128) + tx_samples = np.concatenate((calibration, fixed_guard, bpsk_samples)) + peak = float(np.max(np.abs(tx_samples))) + if peak > 1.0: + raise ValueError(f"Пик TX-последовательности {peak:.6f} превышает единицу") + + return PrbsTransmissionPlan( + label=label, + useful_duration_seconds=float(useful_duration_seconds), + payload_bits=payload_bits, + marker_bits=marker_bits, + pilot_bits=pilot_bits, + tx_samples=tx_samples, + calibration_sample_count=calibration_sample_count, + fixed_guard_sample_count=fixed_guard_sample_count, + bpsk_start_sample=calibration_sample_count + fixed_guard_sample_count, + ) + + +def to_pluto_tx_samples(samples: np.ndarray) -> np.ndarray: + """Преобразовать нормированные комплексные отсчёты в формат Pluto+.""" + + waveform = np.asarray(samples, dtype=np.complex128) + if waveform.ndim != 1 or waveform.size == 0: + raise ValueError("TX-последовательность должна быть непустой и одномерной") + peak = float(np.max(np.abs(waveform))) + if peak > 1.0: + raise ValueError("Амплитуда TX-последовательности превышает единицу") + return (waveform * (2**14)).astype(np.complex64) + + +def locate_calibration_tone_start( + samples: np.ndarray, + sample_rate_hz: float, + block_samples: int = 8_192, + hop_samples: int = 4_096, +) -> int: + """Найти первый сильный участок заранее первой двухтоновой посылки.""" + + received = np.asarray(samples, dtype=np.complex128) + if received.ndim != 1 or len(received) < 4 * block_samples: + raise ValueError("Захват слишком короток для поиска калибровки") + if sample_rate_hz <= 0.0 or block_samples <= 0 or hop_samples <= 0: + raise ValueError("Параметры поиска должны быть положительными") + + centered = received - np.mean(received[: min(len(received), round(0.3 * sample_rate_hz))]) + starts = np.arange(0, len(centered) - block_samples + 1, hop_samples, dtype=np.int64) + powers = np.asarray( + [np.mean(np.abs(centered[start : start + block_samples]) ** 2) for start in starts], + dtype=np.float64, + ) + baseline_limit = max(3, int(0.30 * sample_rate_hz / hop_samples)) + baseline = float(np.median(powers[:baseline_limit])) + threshold = max(baseline * 10.0, np.finfo(np.float64).tiny) + active = powers >= threshold + for index in range(len(active) - 2): + if bool(np.all(active[index : index + 3])): + return int(starts[index]) + raise RuntimeError("Начало двухтоновой калибровки не найдено по фиксированному порогу 10 дБ") + + +def split_calibration_and_bpsk( + samples: np.ndarray, + plan: PrbsTransmissionPlan, + actual_rx_sample_rate_hz: float, + actual_tx_sample_rate_hz: float, +) -> tuple[np.ndarray, np.ndarray, dict]: + """Выделить тоны и BPSK по одной временной шкале одного захвата.""" + + received = np.asarray(samples, dtype=np.complex128) + tone_start = locate_calibration_tone_start(received, actual_rx_sample_rate_hz) + tone_duration_rx = round( + plan.calibration_sample_count / actual_tx_sample_rate_hz * actual_rx_sample_rate_hz + ) + edge = round(0.03 * actual_rx_sample_rate_hz) + calibration_start = tone_start + edge + calibration_end = tone_start + tone_duration_rx - edge + if calibration_end <= calibration_start: + raise RuntimeError("После удаления переходных краёв калибровочный участок пуст") + + bpsk_nominal_start = tone_start + round( + plan.bpsk_start_sample / actual_tx_sample_rate_hz * actual_rx_sample_rate_hz + ) + search_margin = round(0.025 * actual_rx_sample_rate_hz) + bpsk_tx_count = len(plan.tx_samples) - plan.bpsk_start_sample + bpsk_rx_count = round(bpsk_tx_count / actual_tx_sample_rate_hz * actual_rx_sample_rate_hz) + bpsk_start = max(0, bpsk_nominal_start - search_margin) + bpsk_end = min(len(received), bpsk_nominal_start + bpsk_rx_count + search_margin) + if bpsk_end <= bpsk_start: + raise RuntimeError("BPSK-участок не помещается в захват") + return ( + received[calibration_start:calibration_end], + received[bpsk_start:bpsk_end], + { + "tone_start_sample": int(tone_start), + "calibration_start_sample": int(calibration_start), + "calibration_end_sample": int(calibration_end), + "bpsk_nominal_start_sample": int(bpsk_nominal_start), + "bpsk_slice_start_sample": int(bpsk_start), + "bpsk_slice_end_sample": int(bpsk_end), + }, + ) + + +def estimate_refined_calibration( + calibration_samples: np.ndarray, + actual_sample_rate_hz: float, +) -> CalibrationResult: + """Оценить тоны спектром и уточнить обе частоты по фазовому наклону.""" + + coarse = estimate_calibration(calibration_samples, actual_sample_rate_hz) + if not coarse.valid: + return coarse + try: + refined = radio.refine_two_tone_frequencies_from_phase( + calibration_samples, + actual_sample_rate_hz, + coarse.f_low_hz, + coarse.f_high_hz, + block_samples=max(256, round(actual_sample_rate_hz / 10_000.0)), + hop_samples=max(128, round(actual_sample_rate_hz / 20_000.0)), + ) + except (ValueError, RuntimeError) as error: + return _invalid_calibration( + f"фазовое уточнение не выполнено: {error}", + peak_margin_low_db=coarse.peak_margin_low_db, + peak_margin_high_db=coarse.peak_margin_high_db, + noise_power=coarse.noise_power, + ) + offsets = radio.estimate_two_tone_offsets( + refined.low_frequency_hz, + refined.high_frequency_hz, + F_CAL_HZ, + ) + valid = refined.valid and all(np.isfinite(value) for value in offsets.values()) + if not valid: + return _invalid_calibration( + refined.invalid_reason or "фазовое уточнение двух тонов недостоверно", + peak_margin_low_db=coarse.peak_margin_low_db, + peak_margin_high_db=coarse.peak_margin_high_db, + noise_power=coarse.noise_power, + ) + carrier_offset_hz = float(offsets["carrier_offset_hz"]) + return CalibrationResult( + valid=True, + invalid_reason="", + f_low_hz=refined.low_frequency_hz, + f_high_hz=refined.high_frequency_hz, + carrier_offset_hz=carrier_offset_hz, + carrier_offset_ppm=carrier_offset_hz / CARRIER_HZ * 1e6, + clock_scale=float(offsets["clock_scale"]), + sample_clock_error_ppm=float(offsets["sample_clock_error_ppm"]), + residual_cfo_hz=refined.residual_cfo_hz, + phase_fit_rmse_rad=refined.phase_fit_rmse_rad, + peak_margin_low_db=coarse.peak_margin_low_db, + peak_margin_high_db=coarse.peak_margin_high_db, + noise_power=coarse.noise_power, + ) + + +def _invalid_calibration( + reason: str, + *, + peak_margin_low_db: float = float("nan"), + peak_margin_high_db: float = float("nan"), + noise_power: float = float("nan"), +) -> CalibrationResult: + """Вернуть полный недействительный контракт, не опуская поля.""" + + nan = float("nan") + return CalibrationResult( + valid=False, + invalid_reason=str(reason), + f_low_hz=nan, + f_high_hz=nan, + carrier_offset_hz=nan, + carrier_offset_ppm=nan, + clock_scale=nan, + sample_clock_error_ppm=nan, + residual_cfo_hz=nan, + phase_fit_rmse_rad=nan, + peak_margin_low_db=float(peak_margin_low_db), + peak_margin_high_db=float(peak_margin_high_db), + noise_power=float(noise_power), + ) + + +def bit_error_rate(expected_bits: np.ndarray, received_bits: np.ndarray) -> tuple[int, float]: + """Измерить ошибки битов до пакетного разбора и CRC.""" + + expected = np.asarray(expected_bits, dtype=np.uint8) + received = np.asarray(received_bits, dtype=np.uint8) + if expected.ndim != 1 or received.ndim != 1: + raise ValueError("Битовые последовательности должны быть одномерными") + if expected.size == 0 or len(expected) != len(received): + return 0, float("nan") + if not np.all((expected <= 1) & (received <= 1)): + raise ValueError("Последовательности должны содержать только нули и единицы") + errors = int(np.count_nonzero(expected != received)) + return errors, errors / len(expected) + + +def _parabolic_peak_frequency( + frequencies_hz: np.ndarray, + powers: np.ndarray, + peak_index: int, +) -> float: + """Уточнить частоту максимума параболой по логарифму мощности.""" + + if peak_index <= 0 or peak_index >= len(powers) - 1: + return float(frequencies_hz[peak_index]) + local = np.maximum(powers[peak_index - 1 : peak_index + 2], np.finfo(float).tiny) + left, center, right = np.log(local) + denominator = left - 2.0 * center + right + if not np.isfinite(denominator) or abs(denominator) < 1e-15: + return float(frequencies_hz[peak_index]) + offset_bins = 0.5 * (left - right) / denominator + offset_bins = float(np.clip(offset_bins, -1.0, 1.0)) + bin_width_hz = float(frequencies_hz[peak_index + 1] - frequencies_hz[peak_index]) + return float(frequencies_hz[peak_index] + offset_bins * bin_width_hz) + + +def estimate_calibration_from_spectrum( + frequency_axis_hz: np.ndarray, + power_spectrum: np.ndarray, +) -> CalibrationResult: + """Применить зафиксированную методику Lab043 к готовому спектру мощности.""" + + frequencies = np.asarray(frequency_axis_hz, dtype=np.float64) + powers = np.asarray(power_spectrum, dtype=np.float64) + if frequencies.ndim != 1 or powers.ndim != 1 or len(frequencies) != len(powers): + raise ValueError("Ось частот и спектр мощности должны быть одномерными и равными") + + analysis = ( + (np.abs(frequencies) <= CALIBRATION_ANALYSIS_HALF_BAND_HZ) + & np.isfinite(frequencies) + & np.isfinite(powers) + & (powers >= 0.0) + ) + candidates = radio.find_known_tone_peaks( + frequencies[analysis], + powers[analysis], + F_CAL_HZ, + CALIBRATION_SEARCH_HALF_WIDTH_HZ, + ) + low_candidate = float(candidates["low_frequency_hz"]) + high_candidate = float(candidates["high_frequency_hz"]) + + if np.isfinite(low_candidate): + low_index = int(np.nanargmin(np.abs(frequencies - low_candidate))) + low_candidate = _parabolic_peak_frequency(frequencies, powers, low_index) + if np.isfinite(high_candidate): + high_index = int(np.nanargmin(np.abs(frequencies - high_candidate))) + high_candidate = _parabolic_peak_frequency(frequencies, powers, high_index) + + noise_mask = analysis.copy() + if np.isfinite(low_candidate): + noise_mask &= np.abs(frequencies - low_candidate) > CALIBRATION_GUARD_HALF_WIDTH_HZ + if np.isfinite(high_candidate): + noise_mask &= np.abs(frequencies - high_candidate) > CALIBRATION_GUARD_HALF_WIDTH_HZ + noise_values = powers[noise_mask] + noise_power = float(np.median(noise_values)) if noise_values.size else float("nan") + + def excess_db(peak_power: float) -> float: + if not np.isfinite(peak_power) or not np.isfinite(noise_power) or noise_power <= 0.0: + return float("nan") + return float(10.0 * np.log10(peak_power / noise_power)) + + low_excess = excess_db(float(candidates["low_power"])) + high_excess = excess_db(float(candidates["high_power"])) + peaks_valid = ( + np.isfinite(low_candidate) + and np.isfinite(high_candidate) + and np.isfinite(low_excess) + and np.isfinite(high_excess) + and low_excess >= MINIMUM_TONE_EXCESS_DB + and high_excess >= MINIMUM_TONE_EXCESS_DB + ) + + if not peaks_valid: + reason = "два тона не превышают медианный фон на 10 дБ" + return _invalid_calibration( + reason, + peak_margin_low_db=low_excess, + peak_margin_high_db=high_excess, + noise_power=noise_power, + ) + + offsets = radio.estimate_two_tone_offsets(low_candidate, high_candidate, F_CAL_HZ) + valid = all(np.isfinite(offsets[name]) for name in offsets) + if valid and abs(float(offsets["sample_clock_error_ppm"])) > MAXIMUM_ABS_SAMPLE_CLOCK_ERROR_PPM: + return _invalid_calibration( + "разнос тонов означает ошибку такта более 1000 ppm", + peak_margin_low_db=low_excess, + peak_margin_high_db=high_excess, + noise_power=noise_power, + ) + carrier_offset_hz = float(offsets["carrier_offset_hz"]) + if not valid: + return _invalid_calibration( + "невычислимая геометрия двух тонов", + peak_margin_low_db=low_excess, + peak_margin_high_db=high_excess, + noise_power=noise_power, + ) + return CalibrationResult( + valid=True, + invalid_reason="", + f_low_hz=low_candidate, + f_high_hz=high_candidate, + carrier_offset_hz=carrier_offset_hz, + carrier_offset_ppm=carrier_offset_hz / CARRIER_HZ * 1e6, + clock_scale=float(offsets["clock_scale"]), + sample_clock_error_ppm=float(offsets["sample_clock_error_ppm"]), + residual_cfo_hz=float("nan"), + phase_fit_rmse_rad=float("nan"), + peak_margin_low_db=low_excess, + peak_margin_high_db=high_excess, + noise_power=noise_power, + ) + + +def estimate_calibration(samples: np.ndarray, actual_sample_rate_hz: float) -> CalibrationResult: + """Оценить калибровку Lab043 по комплексным отсчётам.""" + + received = np.asarray(samples, dtype=np.complex128) + if received.ndim != 1 or len(received) < CALIBRATION_NFFT: + raise ValueError(f"Нужно не менее {CALIBRATION_NFFT} комплексных отсчётов") + if not np.isfinite(actual_sample_rate_hz) or actual_sample_rate_hz <= 0.0: + raise ValueError("Фактическая частота дискретизации должна быть положительной") + + centered = received - np.mean(received) + frequencies, powers = welch( + centered, + fs=actual_sample_rate_hz, + window="hann", + nperseg=CALIBRATION_NFFT, + noverlap=CALIBRATION_NFFT // 2, + nfft=CALIBRATION_NFFT, + return_onesided=False, + scaling="density", + ) + order = np.argsort(frequencies) + return estimate_calibration_from_spectrum(frequencies[order], powers[order]) + + +def summarize_calibrations(estimates: Iterable[CalibrationResult]) -> dict: + """Свести отдельные оценки без замены невычислимых значений нулём.""" + + rows = tuple(estimates) + if len(rows) != 3: + raise ValueError("Для Lab043 нужны ровно три калибровочных захвата") + fields = ( + "f_low_hz", + "f_high_hz", + "peak_margin_low_db", + "peak_margin_high_db", + "carrier_offset_hz", + "carrier_offset_ppm", + "clock_scale", + "sample_clock_error_ppm", + ) + summary: dict[str, object] = { + "capture_count": len(rows), + "failure_count": sum(not row.valid for row in rows), + } + for field_name in fields: + values = np.asarray([getattr(row, field_name) for row in rows], dtype=np.float64) + finite = values[np.isfinite(values)] + summary[field_name] = { + "mean": float(np.mean(finite)) if finite.size else float("nan"), + "minimum": float(np.min(finite)) if finite.size else float("nan"), + "maximum": float(np.max(finite)) if finite.size else float("nan"), + "standard_deviation": float(np.std(finite, ddof=0)) if finite.size else float("nan"), + } + return summary + + +def synthesize_known_bpsk_capture( + carrier_offset_hz: float, + clock_scale: float, + samples_per_symbol: int = 16, + initial_sample_phase: int = 5, +) -> tuple[np.ndarray, np.ndarray]: + """Синтетически сформировать одну запись маркера и полной PRBS11.""" + + _marker_bits, marker_symbols = radio.build_frame_marker() + payload_bits = prbs11() + symbols = np.concatenate((marker_symbols, radio.bpsk_modulate(payload_bits))) + prefix = np.zeros(3 * samples_per_symbol + initial_sample_phase, dtype=np.complex128) + suffix = np.zeros(3 * samples_per_symbol, dtype=np.complex128) + transmitter_samples = np.concatenate((prefix, np.repeat(symbols, samples_per_symbol), suffix)) + + received_length = max(1, math.floor(len(transmitter_samples) / clock_scale)) + receive_indexes = np.arange(received_length, dtype=np.float64) + transmitter_positions = np.minimum( + receive_indexes * clock_scale, + len(transmitter_samples) - 1.0, + ) + source_indexes = np.arange(len(transmitter_samples), dtype=np.float64) + received = np.interp(transmitter_positions, source_indexes, transmitter_samples.real).astype( + np.complex128 + ) + sample_rate_hz = SYMBOL_RATE * samples_per_symbol + received *= np.exp( + 1j * 2.0 * np.pi * carrier_offset_hz * receive_indexes / sample_rate_hz + ) + return received, payload_bits + + +def process_bpsk_modes( + received_samples: np.ndarray, + expected_payload_bits: np.ndarray, + coarse_carrier_offset_hz: float, + clock_scale: float, + samples_per_symbol: int, +) -> dict[str, ModeResult]: + """Обработать одну запись режимами A, B и C без изменения исходника.""" + + source = np.asarray(received_samples, dtype=np.complex128) + marker_bits, marker_symbols = radio.build_frame_marker() + expected_payload = np.asarray(expected_payload_bits, dtype=np.uint8) + required_symbol_count = len(marker_bits) + len(expected_payload) + sample_rate_hz = SYMBOL_RATE * samples_per_symbol + results: dict[str, ModeResult] = {} + + for mode in ("A", "B", "C"): + processed = source.copy() + if mode in ("B", "C"): + processed = radio.apply_coarse_frequency_correction( + processed, + coarse_carrier_offset_hz, + sample_rate_hz, + ) + if mode == "C": + processed = radio.resample_for_clock_scale(processed, clock_scale) + + try: + found = radio.find_radio_frame(processed, marker_symbols, samples_per_symbol) + start = int(found["start_symbol_index"]) + symbols = found["symbol_samples"][start : start + required_symbol_count] + if len(symbols) != required_symbol_count: + raise RuntimeError("Известная последовательность обрезана") + estimate = radio.estimate_carrier_parameters(symbols[: len(marker_bits)], marker_symbols) + fine_cfo_hz = float(estimate["estimated_frequency_hz"]) + fine_applied = mode == "C" and radio.should_apply_cfo_correction( + fine_cfo_hz, + float(estimate["phase_consistency"]), + float(estimate["coherence_gain"]), + ) + if fine_applied: + corrected_symbols = radio.correct_phase_and_frequency( + symbols, + float(estimate["initial_phase_after_cfo"]), + float(estimate["phase_increment"]), + ) + else: + corrected_symbols = radio.correct_phase_and_frequency( + symbols, + float(estimate["constant_phase"]), + 0.0, + ) + received_payload = radio.bpsk_demodulate(corrected_symbols[len(marker_bits) :]) + errors, ber = bit_error_rate(expected_payload, received_payload) + results[mode] = ModeResult( + mode, + ber, + errors, + len(expected_payload), + float(found["score"]), + int(found["sample_phase"]), + fine_applied, + fine_cfo_hz, + ) + except (RuntimeError, ValueError): + results[mode] = ModeResult( + mode, + float("nan"), + 0, + 0, + 0.0, + 0, + False, + float("nan"), + ) + return results + + +def _maximum_true_run(values: np.ndarray) -> int: + best = current = 0 + for value in np.asarray(values, dtype=bool): + current = current + 1 if value else 0 + best = max(best, current) + return int(best) + + +def _best_local_timing_offset( + matched_samples: np.ndarray, + nominal_start_sample: int, + expected_symbols: np.ndarray, + samples_per_symbol: int, +) -> float: + """Оценить локальное смещение символьной сетки без повторной синхронизации.""" + + expected = np.asarray(expected_symbols, dtype=np.complex128) + if expected.size == 0: + return float("nan") + count = min(len(expected), 2_048) + expected = expected[:count] + best_offset = 0 + best_score = -1.0 + for offset in range(-samples_per_symbol // 2, samples_per_symbol // 2 + 1): + indexes = nominal_start_sample + offset + np.arange(count) * samples_per_symbol + if indexes[0] < 0 or indexes[-1] >= len(matched_samples): + continue + observed = matched_samples[indexes] + denominator = math.sqrt( + float(np.sum(np.abs(observed) ** 2) * np.sum(np.abs(expected) ** 2)) + ) + 1e-12 + score = float(abs(np.vdot(expected, observed)) / denominator) + if score > best_score: + best_score = score + best_offset = offset + if best_score < 0.0: + return float("nan") + return float(best_offset) + + +def _failed_prbs_metrics(mode: str, transmitted_bit_count: int, reason: str) -> PrbsModeMetrics: + nan = float("nan") + return PrbsModeMetrics( + mode=mode, + detected=False, + failure_reason=reason, + transmitted_bit_count=int(transmitted_bit_count), + matched_bit_count=0, + bit_errors=0, + bit_error_rate=nan, + quarter_bit_error_rates=(nan, nan, nan, nan), + evm_percent=nan, + marker_correlation=0.0, + symbol_phase_start_samples=nan, + symbol_phase_end_samples=nan, + accumulated_timing_drift_samples=nan, + accumulated_timing_drift_symbols=nan, + residual_cfo_hz=nan, + maximum_correct_run_bits=0, + fine_cfo_applied=False, + estimated_fine_cfo_hz=nan, + timing_sro_ppm=nan, + coarse_cfo_applied=False, + pilot_cfo_applied=False, + pilot_estimate_valid=False, + estimated_pilot_cfo_hz=nan, + residual_pilot_cfo_hz=nan, + pilot_phase_fit_rmse_rad=nan, + pilot_mean_block_coherence=nan, + legacy_prbs_adjacent_phase_hz=nan, + bit_pattern_diagnostics={}, + ) + + +def _prbs_bit_pattern_diagnostics( + payload_symbols: np.ndarray, + expected_bits: np.ndarray, +) -> dict: + """Сгруппировать диагностические ошибки по тройкам соседних битов.""" + + symbols = np.asarray(payload_symbols, dtype=np.complex128) + bits = np.asarray(expected_bits, dtype=np.uint8) + expected_symbols = radio.bpsk_modulate(bits) + if len(symbols) != len(bits) or len(bits) < 3: + return {} + gain = np.vdot(expected_symbols, symbols) / ( + np.vdot(expected_symbols, expected_symbols) + 1e-12 + ) + if abs(gain) <= 1e-12: + return {} + normalized_despread = symbols / gain * expected_symbols + received_bits = radio.bpsk_demodulate(symbols) + error_mask = received_bits != bits + groups: dict[str, dict] = {} + indexes = np.arange(1, len(bits) - 1) + for previous in (0, 1): + for current in (0, 1): + for following in (0, 1): + selected = indexes[ + (bits[indexes - 1] == previous) + & (bits[indexes] == current) + & (bits[indexes + 1] == following) + ] + values = normalized_despread[selected] + errors = int(np.count_nonzero(error_mask[selected])) + groups[f"{previous}{current}{following}"] = { + "count": int(len(selected)), + "bit_errors": errors, + "bit_error_rate": errors / len(selected) if len(selected) else float("nan"), + "mean_real": float(np.mean(values.real)) if len(values) else float("nan"), + "mean_imag": float(np.mean(values.imag)) if len(values) else float("nan"), + "mean_magnitude": float(np.mean(np.abs(values))) if len(values) else float("nan"), + } + return { + "gain_real": float(gain.real), + "gain_imag": float(gain.imag), + "groups": groups, + } + + +def analyze_prbs_modes( + bpsk_samples: np.ndarray, + expected_payload_bits: np.ndarray, + coarse_carrier_offset_hz: float, + clock_scale: float, + actual_sample_rate_hz: float, + samples_per_symbol: int = SAMPLES_PER_SYMBOL, + known_pilot_bits: np.ndarray | None = None, +) -> dict[str, PrbsModeMetrics]: + """Измерить A/B/C/D на одной записи без знания полезной PRBS при коррекции. + + A: без грубой CFO и SRO; B: только грубая CFO; C: грубая CFO и + отдельный пилот; D: то же, что C, плюс компенсация SRO. Известная PRBS + используется только после рабочей коррекции для расчёта метрик. + """ + + source = np.asarray(bpsk_samples, dtype=np.complex128) + expected_bits = np.asarray(expected_payload_bits, dtype=np.uint8) + pilot_bits = ( + pilot_prbs7() + if known_pilot_bits is None + else np.asarray(known_pilot_bits, dtype=np.uint8) + ) + if pilot_bits.ndim != 1 or len(pilot_bits) == 0: + raise ValueError("Для короткого S нужен отдельный непустой известный пилот") + marker_bits, marker_symbols = radio.build_frame_marker() + pilot_symbols = radio.bpsk_modulate(pilot_bits) + expected_payload_symbols = radio.bpsk_modulate(expected_bits) + pilot_start_symbol = len(marker_bits) + payload_start_symbol = pilot_start_symbol + len(pilot_bits) + required_symbol_count = payload_start_symbol + len(expected_bits) + taps = radio.root_raised_cosine_taps( + LAB042_RRC_ROLLOFF, + samples_per_symbol, + LAB042_RRC_SPAN_SYMBOLS, + ) + results: dict[str, PrbsModeMetrics] = {} + + for mode in ("A", "B", "C", "D"): + try: + processed = source.copy() + coarse_applied = mode in ("B", "C", "D") + pilot_applied = mode in ("C", "D") + if coarse_applied: + processed = radio.apply_coarse_frequency_correction( + processed, + coarse_carrier_offset_hz, + actual_sample_rate_hz, + ) + if mode == "D": + processed = radio.resample_for_clock_scale(processed, clock_scale) + + matched = fftconvolve(processed, taps, mode="full") + found = radio.find_radio_frame(matched, marker_symbols, samples_per_symbol) + marker_correlation = float(found["score"]) + if marker_correlation < PRBS_MARKER_MINIMUM_CORRELATION: + raise RuntimeError( + f"корреляция маркера {marker_correlation:.4f} ниже фиксированного порога " + f"{PRBS_MARKER_MINIMUM_CORRELATION:.2f}" + ) + + start_symbol = int(found["start_symbol_index"]) + symbols = found["symbol_samples"][start_symbol : start_symbol + required_symbol_count] + if len(symbols) != required_symbol_count: + raise RuntimeError("маркер, пилот или полезная PRBS11 обрезаны") + + marker_carrier = radio.estimate_carrier_parameters( + symbols[: len(marker_bits)], + marker_symbols, + ) + constant_corrected = radio.correct_phase_and_frequency( + symbols, + float(marker_carrier["constant_phase"]), + 0.0, + ) + + pilot_estimate = None + if coarse_applied: + pilot_estimate = radio.estimate_known_pilot_carrier( + constant_corrected[pilot_start_symbol:payload_start_symbol], + pilot_symbols, + ) + if pilot_applied and (pilot_estimate is None or not pilot_estimate.valid): + reason = ( + "пилотная оценка отсутствует" + if pilot_estimate is None + else pilot_estimate.invalid_reason + ) + raise RuntimeError(f"пилотная оценка недостоверна: {reason}") + + if pilot_applied and pilot_estimate is not None: + corrected_tail = radio.correct_phase_and_frequency( + constant_corrected[pilot_start_symbol:], + pilot_estimate.initial_phase_rad, + pilot_estimate.phase_increment_rad_per_symbol, + ) + corrected_pilot = corrected_tail[: len(pilot_bits)] + payload_symbols = corrected_tail[len(pilot_bits) :] + residual_estimate = radio.estimate_known_pilot_carrier( + corrected_pilot, + pilot_symbols, + ) + else: + corrected_pilot = constant_corrected[ + pilot_start_symbol:payload_start_symbol + ] + payload_symbols = constant_corrected[payload_start_symbol:] + residual_estimate = pilot_estimate + + pilot_estimate_valid = bool( + pilot_estimate is not None and pilot_estimate.valid + ) + estimated_pilot_cfo_hz = ( + float(pilot_estimate.frequency_hz) + if pilot_estimate_valid and pilot_estimate is not None + else float("nan") + ) + residual_pilot_cfo_hz = ( + float(residual_estimate.frequency_hz) + if residual_estimate is not None and residual_estimate.valid + else float("nan") + ) + pilot_phase_fit_rmse_rad = ( + float(pilot_estimate.phase_fit_rmse_rad) + if pilot_estimate is not None + else float("nan") + ) + pilot_mean_block_coherence = ( + float(pilot_estimate.mean_block_coherence) + if pilot_estimate is not None + else float("nan") + ) + + if len(payload_symbols) != len(expected_bits): + raise RuntimeError("полезная PRBS11 обрезана после коррекции") + received_bits = radio.bpsk_demodulate(payload_symbols) + errors, ber = bit_error_rate(expected_bits, received_bits) + + quarter_bers: list[float] = [] + for quarter in np.array_split(np.arange(len(expected_bits)), 4): + quarter_errors = int( + np.count_nonzero(expected_bits[quarter] != received_bits[quarter]) + ) + quarter_bers.append(quarter_errors / len(quarter)) + + gain = np.vdot(expected_payload_symbols, payload_symbols) / ( + np.vdot(expected_payload_symbols, expected_payload_symbols) + 1e-12 + ) + reference = gain * expected_payload_symbols + evm_percent = float( + 100.0 + * math.sqrt( + float(np.mean(np.abs(payload_symbols - reference) ** 2)) + / (float(np.mean(np.abs(reference) ** 2)) + 1e-12) + ) + ) + + # Историческая величина оставлена только с явным названием. Она + # не используется как CFO и не участвует ни в одной коррекции. + diagnostic_despread = payload_symbols * expected_payload_symbols + diagnostic_adjacent = ( + diagnostic_despread[1:] * np.conj(diagnostic_despread[:-1]) + ) + legacy_prbs_adjacent_phase_hz = float( + np.angle(np.sum(diagnostic_adjacent)) * SYMBOL_RATE / (2.0 * np.pi) + ) + + marker_start_sample = int(found["sample_phase"]) + start_symbol * samples_per_symbol + payload_start_sample = marker_start_sample + payload_start_symbol * samples_per_symbol + timing_window = min(2_048, max(1, len(expected_bits) // 4)) + start_offset = _best_local_timing_offset( + matched, + payload_start_sample, + expected_payload_symbols[:timing_window], + samples_per_symbol, + ) + end_bit = len(expected_bits) - timing_window + end_offset = _best_local_timing_offset( + matched, + payload_start_sample + end_bit * samples_per_symbol, + expected_payload_symbols[end_bit:], + samples_per_symbol, + ) + timing_drift = end_offset - start_offset + timing_denominator = max(1, end_bit) * samples_per_symbol + timing_sro_ppm = -timing_drift / timing_denominator * 1e6 + correct = received_bits == expected_bits + + results[mode] = PrbsModeMetrics( + mode=mode, + detected=True, + failure_reason="", + transmitted_bit_count=len(expected_bits), + matched_bit_count=len(received_bits), + bit_errors=errors, + bit_error_rate=ber, + quarter_bit_error_rates=tuple(float(value) for value in quarter_bers), + evm_percent=evm_percent, + marker_correlation=marker_correlation, + symbol_phase_start_samples=float(found["sample_phase"]) + start_offset, + symbol_phase_end_samples=float(found["sample_phase"]) + end_offset, + accumulated_timing_drift_samples=timing_drift, + accumulated_timing_drift_symbols=timing_drift / samples_per_symbol, + residual_cfo_hz=residual_pilot_cfo_hz, + maximum_correct_run_bits=_maximum_true_run(correct), + fine_cfo_applied=bool(pilot_applied), + estimated_fine_cfo_hz=estimated_pilot_cfo_hz, + timing_sro_ppm=timing_sro_ppm, + coarse_cfo_applied=coarse_applied, + pilot_cfo_applied=pilot_applied, + pilot_estimate_valid=pilot_estimate_valid, + estimated_pilot_cfo_hz=estimated_pilot_cfo_hz, + residual_pilot_cfo_hz=residual_pilot_cfo_hz, + pilot_phase_fit_rmse_rad=pilot_phase_fit_rmse_rad, + pilot_mean_block_coherence=pilot_mean_block_coherence, + legacy_prbs_adjacent_phase_hz=legacy_prbs_adjacent_phase_hz, + bit_pattern_diagnostics=_prbs_bit_pattern_diagnostics( + payload_symbols, + expected_bits, + ), + ) + except (RuntimeError, ValueError) as error: + results[mode] = _failed_prbs_metrics(mode, len(expected_bits), str(error)) + return results + + +def control_packet_crc_roundtrip(payload: bytes = b"Lab043 control packet") -> bool: + """Проверить один настоящий пакет с CRC внутри радиокадра.""" + + packet = build_packet(payload, MESSAGE_TYPE_TEXT, sequence_number=43) + bits, _, _ = radio.build_radio_frame(packet) + recovered_packet = radio.parse_radio_frame(bits) + return recovered_packet is not None and parse_packet(recovered_packet).payload == payload + + +def try_reassemble_complete_image(fragments: Iterable) -> bytes | None: + """Собрать изображение либо явно вернуть отсутствие полного результата.""" + + try: + return reassemble_image(fragments) + except MissingFragmentsError: + return None + + +def exit_code_for_acceptance(result: AcceptanceResult) -> int: + """Вернуть ноль только по полному критерию Lab043.""" + + return 0 if result.passed else 1 + + +def save_reference_iq_capture( + base_path: Path, + samples: np.ndarray, + metadata: dict, +) -> tuple[Path, Path]: + """Сохранить эталонный IQ и воспроизводимые метаданные. + + Функция только подготавливает поддержку формата. Аппаратный файл в + текущем программном этапе не создаётся. + """ + + iq = np.asarray(samples, dtype=np.complex64) + if iq.ndim != 1: + raise ValueError("IQ-захват должен быть одномерным") + required = { + "actual_sample_rate_hz", + "center_frequency_hz", + "rx_gain_db", + "capture_started_utc", + "capture_order", + "tx_parameters", + "calibration", + } + missing = sorted(required - set(metadata)) + if missing: + raise ValueError("Не хватает метаданных: " + ", ".join(missing)) + + base_path = Path(base_path) + base_path.parent.mkdir(parents=True, exist_ok=True) + iq_path = base_path.with_suffix(".npy") + metadata_path = base_path.with_suffix(".json") + np.save(iq_path, iq, allow_pickle=False) + file_sha256 = hashlib.sha256(iq_path.read_bytes()).hexdigest() + document = dict(metadata) + document.update( + { + "sample_count": len(iq), + "dtype": "complex64", + "iq_file": iq_path.name, + "iq_sha256": file_sha256, + } + ) + metadata_path.write_text( + json.dumps(document, ensure_ascii=False, indent=2, sort_keys=True) + "\n", + encoding="utf-8", + ) + return iq_path, metadata_path + + +def _write_json_atomically(path: Path, document: dict) -> None: + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text( + json.dumps(document, ensure_ascii=False, indent=2, sort_keys=True) + "\n", + encoding="utf-8", + ) + temporary.replace(path) + + +def save_raw_iq_before_processing( + base_path: Path, + samples: np.ndarray, + metadata: dict, +) -> tuple[Path, Path, dict]: + """Сохранить и повторно проверить raw до запуска любого анализа.""" + + iq = np.asarray(samples, dtype=np.complex64) + if iq.ndim != 1 or iq.size == 0: + raise ValueError("Raw IQ должен быть непустым одномерным массивом") + if not np.all(np.isfinite(iq)): + raise ValueError("Raw IQ содержит нечисловые значения") + required = { + "capture_status", + "processing_status", + "processing_error", + "capture_started_utc", + "tx_parameters", + "rx_parameters", + "waveform_sha256", + "callback_count", + "callback_boundaries", + "tx_waveform_duration_seconds", + "expected_intervals", + } + missing = sorted(required - set(metadata)) + if missing: + raise ValueError("Не хватает raw-метаданных: " + ", ".join(missing)) + if metadata["capture_status"] != "captured": + raise ValueError("До обработки capture_status должен быть captured") + if metadata["processing_status"] != "pending": + raise ValueError("До обработки processing_status должен быть pending") + + base_path = Path(base_path) + base_path.parent.mkdir(parents=True, exist_ok=True) + iq_path = base_path.with_suffix(".npy") + metadata_path = base_path.with_suffix(".json") + + # Порядок обязателен: NPY -> первичный JSON -> SHA-256 -> проверка NPY. + np.save(iq_path, iq, allow_pickle=False) + document = dict(metadata) + document.update( + { + "sample_count": int(len(iq)), + "dtype": "complex64", + "iq_file": iq_path.name, + "iq_sha256": None, + "raw_verified": False, + } + ) + _write_json_atomically(metadata_path, document) + iq_sha256 = hashlib.sha256(iq_path.read_bytes()).hexdigest() + restored = np.load(iq_path, allow_pickle=False) + verified = ( + restored.dtype == np.dtype(np.complex64) + and restored.shape == iq.shape + and np.array_equal(restored, iq) + and hashlib.sha256(iq_path.read_bytes()).hexdigest() == iq_sha256 + ) + if not verified: + raise RuntimeError("Повторная проверка сохранённого raw IQ не прошла") + document["iq_sha256"] = iq_sha256 + document["raw_verified"] = True + _write_json_atomically(metadata_path, document) + return iq_path, metadata_path, document + + +def update_processing_metadata( + metadata_path: Path, + *, + processing_status: str, + processing_error: str | None, + analysis: dict | None = None, +) -> dict: + """Безопасно обновить только состояние обработки уже сохранённого raw.""" + + if processing_status not in {"success", "processing_failed"}: + raise ValueError("Неизвестный processing_status") + path = Path(metadata_path) + document = json.loads(path.read_text(encoding="utf-8")) + document["processing_status"] = processing_status + document["processing_error"] = processing_error + if analysis is not None: + document["analysis"] = analysis + _write_json_atomically(path, document) + return document + + +def save_then_process_capture( + base_path: Path, + samples: np.ndarray, + metadata: dict, + processor, +) -> tuple[Path, Path, object | None, str | None]: + """Зафиксировать raw, затем обработать; ошибку анализа оставить в JSON.""" + + iq_path, metadata_path, _document = save_raw_iq_before_processing( + base_path, + samples, + metadata, + ) + try: + result = processor(np.asarray(samples, dtype=np.complex64)) + except Exception as error: + message = f"{type(error).__name__}: {error}" + update_processing_metadata( + metadata_path, + processing_status="processing_failed", + processing_error=message, + ) + return iq_path, metadata_path, None, message + update_processing_metadata( + metadata_path, + processing_status="success", + processing_error=None, + analysis=asdict(result) if hasattr(result, "__dataclass_fields__") else result, + ) + return iq_path, metadata_path, result, None + + +def save_diagnostic_artifacts( + output_directory: Path, + calibrations: Iterable[CalibrationResult], + mode_results: dict[str, ModeResult], + acceptance: AcceptanceResult, +) -> tuple[Path, ...]: + """Сохранить CSV, TXT и PNG диагностики без обращения к аппаратуре.""" + + output_directory = Path(output_directory) + output_directory.mkdir(parents=True, exist_ok=True) + calibration_rows = tuple(calibrations) + + calibration_csv = output_directory / "lab043_calibration.csv" + with calibration_csv.open("w", encoding="utf-8", newline="") as stream: + field_names = list(CalibrationResult.__dataclass_fields__) + writer = csv.DictWriter(stream, fieldnames=field_names) + writer.writeheader() + for row in calibration_rows: + writer.writerow(asdict(row)) + + modes_csv = output_directory / "lab043_modes.csv" + with modes_csv.open("w", encoding="utf-8", newline="") as stream: + field_names = list(ModeResult.__dataclass_fields__) + writer = csv.DictWriter(stream, fieldnames=field_names) + writer.writeheader() + for mode in ("A", "B", "C"): + if mode in mode_results: + writer.writerow(asdict(mode_results[mode])) + + summary_csv = output_directory / "lab043_summary.csv" + summary_row = asdict(acceptance) | { + "radio_passed": acceptance.radio_passed, + "application_passed": acceptance.application_passed, + "passed": acceptance.passed, + "exit_code": exit_code_for_acceptance(acceptance), + } + with summary_csv.open("w", encoding="utf-8", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=list(summary_row)) + writer.writeheader() + writer.writerow(summary_row) + + report_path = output_directory / "lab043_report.txt" + report_lines = [ + "Lab043: программная диагностика", + f"Калибровочных захватов: {len(calibration_rows)}.", + f"Радиокритерий: {'пройден' if acceptance.radio_passed else 'не пройден'}.", + f"Прикладной критерий: {'пройден' if acceptance.application_passed else 'не пройден'}.", + f"Программные проверки: {'пройдены' if acceptance.software_tests_passed else 'не пройдены'}.", + f"Код возврата: {exit_code_for_acceptance(acceptance)}.", + "Аппаратный тракт этим отчётом не подтверждается.", + ] + report_path.write_text("\n".join(report_lines) + "\n", encoding="utf-8") + + plot_path = output_directory / "lab043_calibration.png" + figure = Figure(figsize=(7, 4)) + axis = figure.subplots() + capture_indexes = np.arange(1, len(calibration_rows) + 1) + low = [row.f_low_hz for row in calibration_rows] + high = [row.f_high_hz for row in calibration_rows] + axis.plot(capture_indexes, low, "o-", label="Нижний тон") + axis.plot(capture_indexes, high, "o-", label="Верхний тон") + axis.set_xlabel("Номер калибровочного захвата") + axis.set_ylabel("Частота относительно несущей, Гц") + axis.grid(True, alpha=0.3) + axis.legend() + figure.tight_layout() + figure.savefig(plot_path, dpi=140) + + return calibration_csv, modes_csv, summary_csv, report_path, plot_path + + +def _run_check(name: str, function) -> FunctionalTestResult: + try: + detail = function() + return FunctionalTestResult(name, True, str(detail)) + except Exception as error: + return FunctionalTestResult(name, False, f"{type(error).__name__}: {error}") + + +def run_functional_tests() -> tuple[FunctionalTestResult, ...]: + """Выполнить быстрые синтетические проверки без аппаратуры и файлов.""" + + def duration_is_derived() -> str: + duration = receive_duration_seconds(2_400_000, 2_400_000.0) + assert duration == 2.0 + assert receive_sample_count(2_400_000, 2_400_000.0, 2_400_000.0) == 4_800_000 + return "длительность получена из числа TX-отсчётов" + + def blocks_are_continuous() -> str: + blocks = (np.arange(4), np.arange(4, 9)) + combined = combine_continuous_blocks(blocks, 8) + assert np.array_equal(combined.real, np.arange(8)) + return "последовательные блоки образуют одну запись" + + def prbs_has_full_period() -> str: + sequence = prbs11() + assert len(sequence) == PRBS11_LENGTH + assert not np.array_equal(sequence, np.roll(sequence, 1)) + return "PRBS11 содержит 2047 детерминированных битов" + + def separate_pilot_is_fixed() -> str: + pilot = pilot_prbs7() + assert len(pilot) == PILOT_SYMBOL_COUNT + assert np.array_equal(pilot[:127], pilot[127:254]) + assert int(np.count_nonzero(pilot[:127])) == 64 + assert len(pilot) / SYMBOL_RATE == 0.064 + return "отдельный PRBS7-пилот содержит 1280 символов и длится 64 мс" + + def cfo_signs_are_correct() -> str: + sample_rate = 10_000.0 + indexes = np.arange(10_000) + for offset in (-700.0, 700.0): + impaired = np.exp(1j * 2.0 * np.pi * offset * indexes / sample_rate) + corrected = radio.apply_coarse_frequency_correction(impaired, offset, sample_rate) + assert np.max(np.abs(corrected - 1.0)) < 1e-9 + return "грубая CFO обоих знаков компенсируется правильным знаком" + + def clock_scale_cases_are_correct() -> str: + for ppm in (20.0, -20.0, 100.0, -100.0): + scale = 1.0 + ppm * 1e-6 + low = -F_CAL_HZ * scale + high = F_CAL_HZ * scale + estimate = radio.estimate_two_tone_offsets(low, high, F_CAL_HZ) + assert math.isclose(estimate["sample_clock_error_ppm"], ppm, abs_tol=1e-6) + source = np.arange(100_000, dtype=np.complex128) + corrected = radio.resample_for_clock_scale(source, scale) + assert len(corrected) == round(len(source) * scale) + return "знак, величина и направление проверены для +/-20 и +/-100 ppm" + + def weak_tones_are_rejected() -> str: + frequencies = np.linspace(-100_000.0, 100_000.0, 4001) + powers = np.ones_like(frequencies) + powers[np.argmin(np.abs(frequencies + F_CAL_HZ))] = 5.0 + powers[np.argmin(np.abs(frequencies - F_CAL_HZ))] = 5.0 + estimate = estimate_calibration_from_spectrum(frequencies, powers) + assert not estimate.valid and math.isnan(estimate.carrier_offset_hz) + return "слабые тоны дают явную невычислимую оценку" + + def packet_crc_is_valid() -> str: + assert control_packet_crc_roundtrip() + return "контрольный пакет прошёл CRC" + + def modes_use_one_capture() -> str: + samples_per_symbol = 16 + sample_rate = SYMBOL_RATE * samples_per_symbol + plan = build_prbs11_transmission_plan( + "S", + 0.025, + sample_rate_hz=sample_rate, + samples_per_symbol=samples_per_symbol, + ) + capture = plan.tx_samples[plan.bpsk_start_sample :] + indexes = np.arange(len(capture), dtype=np.float64) + capture = capture * np.exp( + 1j * 2.0 * np.pi * 701.2 * indexes / sample_rate + ) + results = analyze_prbs_modes( + capture, + plan.payload_bits, + 700.0, + 1.0, + sample_rate, + samples_per_symbol=samples_per_symbol, + known_pilot_bits=plan.pilot_bits, + ) + assert set(results) == {"A", "B", "C", "D"} + assert results["C"].pilot_estimate_valid + assert results["C"].bit_error_rate == 0.0 + return "режим C восстановил PRBS11 по маркеру и отдельному пилоту" + + def incomplete_acceptance_fails() -> str: + complete = AcceptanceResult(10, 10, 10, 10, True, True, True, True) + incomplete = AcceptanceResult(10, 10, 10, 9, False, False, False, True) + assert exit_code_for_acceptance(complete) == 0 + assert exit_code_for_acceptance(incomplete) == 1 + return "код 0 выдаётся только по полному критерию" + + checks = ( + ("01. Расчёт непрерывного захвата", duration_is_derived), + ("02. Склейка блоков без разрывов", blocks_are_continuous), + ("03. Известная PRBS11", prbs_has_full_period), + ("04. Отдельный пилот", separate_pilot_is_fixed), + ("05. Знак грубой CFO", cfo_signs_are_correct), + ("06. Направление clock_scale", clock_scale_cases_are_correct), + ("07. Порог двух тонов", weak_tones_are_rejected), + ("08. Один пакет с CRC", packet_crc_is_valid), + ("09. Режимы A/B/C/D", modes_use_one_capture), + ("10. Полный критерий успеха", incomplete_acceptance_fails), + ) + return tuple(_run_check(name, function) for name, function in checks) + + +def main() -> int: + """Запустить только программные проверки Lab043.""" + + results = run_functional_tests() + for result in results: + print(f"{'PASS' if result.passed else 'FAIL'} | {result.name} | {result.detail}") + passed = all(result.passed for result in results) + print("Аппаратный тракт не запускался.") + return 0 if passed else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/experiments/lab043_prbs11_hardware.py b/experiments/lab043_prbs11_hardware.py new file mode 100644 index 0000000..eceec2d --- /dev/null +++ b/experiments/lab043_prbs11_hardware.py @@ -0,0 +1,485 @@ +"""Один разрешённый аппаратный захват PRBS11 для Lab043. + +Скрипт выполняет ровно один из заранее определённых опытов S или L. Он не +подбирает параметры, не повторяет передачу при ошибке и не переходит к +пакетам, CRC либо JPEG. +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import math +import threading +import time +from dataclasses import asdict, dataclass +from datetime import datetime, timezone +from pathlib import Path + +import numpy as np + +from experiments import lab043_pluto_to_rtlsdr as lab043 + + +PLUTO_URI = "ip:192.168.2.1" +TX_GAIN_DB = -30.0 +RTL_GAIN_DB = 19.7 +RF_BANDWIDTH_HZ = 200_000 +CONTINUITY_MAX_PHASE_JUMP_RAD = 0.25 + + +@dataclass(frozen=True) +class CapturedPrbsAnalysis: + """Полный результат обработки одного уже принятого S/L.""" + + processing_success: bool + invalid_reason: str + sections: dict + calibration: lab043.CalibrationResult + continuity: dict + modes: dict[str, lab043.PrbsModeMetrics] + + +def preflight_pluto() -> dict: + """Подтвердить IIO-контекст, PHY и первый TX-канал без передачи.""" + + import iio + + context = iio.Context(PLUTO_URI) + phy = context.find_device("ad9361-phy") + tx_device = context.find_device("cf-ad9361-dds-core-lpc") + tx_channel = tx_device.find_channel("voltage0", True) if tx_device is not None else None + if phy is None: + raise RuntimeError("IIO-контекст открыт, но ad9361-phy не найден") + if tx_device is None or tx_channel is None: + raise RuntimeError("IIO-контекст открыт, но первый TX-канал недоступен") + return { + "uri": PLUTO_URI, + "context_opened": True, + "phy_name": str(phy.name), + "tx_device_name": str(tx_device.name), + "first_tx_channel": str(tx_channel.id), + } + + +def _configure_rtlsdr(sdr) -> dict: + sdr.sample_rate = float(lab043.SAMPLE_RATE_HZ) + sdr.center_freq = float(lab043.CARRIER_HZ) + sdr.gain = float(RTL_GAIN_DB) + actual_sample_rate = float(sdr.sample_rate) + actual_center = float(sdr.center_freq) + raw_gain = float(sdr.gain) + actual_gain = raw_gain if math.isclose(raw_gain, RTL_GAIN_DB, abs_tol=0.2) else None + return { + "requested_center_frequency_hz": float(lab043.CARRIER_HZ), + "actual_center_frequency_hz": actual_center, + "requested_sample_rate_hz": float(lab043.SAMPLE_RATE_HZ), + "actual_sample_rate_hz": actual_sample_rate, + "requested_gain_db": RTL_GAIN_DB, + "raw_gain_readback_db": raw_gain, + "actual_gain_db": actual_gain, + } + + +def preflight_rtlsdr() -> dict: + """Открыть RTL-SDR, прочитать короткий блок и закрыть до TX.""" + + from rtlsdr import RtlSdr + + sdr = RtlSdr(device_index=0) + try: + parameters = _configure_rtlsdr(sdr) + probe = np.asarray(sdr.read_samples(16_384), dtype=np.complex64) + if len(probe) != 16_384 or not np.all(np.isfinite(probe)): + raise RuntimeError("RTL-SDR не вернул полный конечный пробный блок") + parameters["probe_sample_count"] = int(len(probe)) + parameters["probe_rms"] = float(np.sqrt(np.mean(np.abs(probe) ** 2))) + return parameters + finally: + sdr.close() + + +def configure_pluto() -> tuple[object, dict]: + import adi + + device = adi.Pluto(uri=PLUTO_URI) + device.sample_rate = int(lab043.SAMPLE_RATE_HZ) + device.tx_lo = int(lab043.CARRIER_HZ) + device.tx_rf_bandwidth = int(RF_BANDWIDTH_HZ) + device.tx_hardwaregain_chan0 = float(TX_GAIN_DB) + device.tx_cyclic_buffer = False + return device, { + "requested_center_frequency_hz": int(lab043.CARRIER_HZ), + "actual_center_frequency_hz": int(device.tx_lo), + "requested_sample_rate_hz": int(lab043.SAMPLE_RATE_HZ), + "actual_sample_rate_hz": int(device.sample_rate), + "requested_gain_db": TX_GAIN_DB, + "actual_gain_db": float(device.tx_hardwaregain_chan0), + "requested_rf_bandwidth_hz": RF_BANDWIDTH_HZ, + "actual_rf_bandwidth_hz": int(device.tx_rf_bandwidth), + } + + +def transmit_noncyclic_buffer_once( + device, + device_samples: np.ndarray, + actual_sample_rate_hz: float, + sleep_function=time.sleep, +) -> float: + """Передать один нециклический буфер и не уничтожать его раньше времени. + + В libiio v1 постановка блока в поток может завершиться до того, как DMA + физически выведет все отсчёты. Поэтому буфер остаётся жив не меньше его + расчётной длительности. Возвращаемое значение сохраняется в диагностике. + """ + + samples = np.asarray(device_samples, dtype=np.complex64) + if samples.ndim != 1 or samples.size == 0: + raise ValueError("TX-буфер должен быть непустым и одномерным") + if not np.isfinite(actual_sample_rate_hz) or actual_sample_rate_hz <= 0.0: + raise ValueError("Фактическая частота TX должна быть положительной") + hold_seconds = len(samples) / actual_sample_rate_hz + device.tx_destroy_buffer() + device.tx(samples) + sleep_function(hold_seconds) + device.tx_destroy_buffer() + return float(hold_seconds) + + +def calibration_boundary_continuity( + samples: np.ndarray, + callback_boundaries, + sections: dict, + calibration: lab043.CalibrationResult, + sample_rate_hz: float, +) -> dict: + """Проверить фазу обоих тонов на callback-границах внутри калибровки.""" + + relevant = [ + int(row.accepted_end_sample) + for row in callback_boundaries + if sections["calibration_start_sample"] + 4_096 + <= row.accepted_end_sample + <= sections["calibration_end_sample"] - 4_096 + ] + jumps: list[dict] = [] + for boundary in relevant: + row = {"boundary_sample": boundary} + for name, frequency_hz in ( + ("low", calibration.f_low_hz), + ("high", calibration.f_high_hz), + ): + before_indexes = np.arange(boundary - 4_096, boundary, dtype=np.float64) + after_indexes = np.arange(boundary, boundary + 4_096, dtype=np.float64) + before = np.mean( + samples[boundary - 4_096 : boundary] + * np.exp(-1j * 2.0 * np.pi * frequency_hz * before_indexes / sample_rate_hz) + ) + after = np.mean( + samples[boundary : boundary + 4_096] + * np.exp(-1j * 2.0 * np.pi * frequency_hz * after_indexes / sample_rate_hz) + ) + row[f"{name}_phase_jump_rad"] = float(np.angle(after * np.conj(before))) + jumps.append(row) + maximum = max( + ( + abs(value) + for row in jumps + for key, value in row.items() + if key.endswith("phase_jump_rad") + ), + default=0.0, + ) + return { + "checked_boundary_count": len(relevant), + "maximum_absolute_phase_jump_rad": float(maximum), + "fixed_limit_rad": CONTINUITY_MAX_PHASE_JUMP_RAD, + "confirmed_discontinuity_count": int( + sum( + max(abs(row["low_phase_jump_rad"]), abs(row["high_phase_jump_rad"])) + > CONTINUITY_MAX_PHASE_JUMP_RAD + for row in jumps + ) + ), + "boundaries": jumps, + } + + +def analyze_captured_prbs( + samples: np.ndarray, + plan: lab043.PrbsTransmissionPlan, + actual_rx_rate: float, + actual_tx_rate: float, + callback_boundaries=(), +) -> CapturedPrbsAnalysis: + """Пройти тем же полным путём, что аппаратный S, но без обращения к SDR.""" + + calibration_samples, bpsk_samples, sections = lab043.split_calibration_and_bpsk( + samples, + plan, + actual_rx_rate, + actual_tx_rate, + ) + calibration = lab043.estimate_refined_calibration(calibration_samples, actual_rx_rate) + if not calibration.valid: + reason = f"калибровка недостоверна: {calibration.invalid_reason}" + modes = { + mode: lab043._failed_prbs_metrics(mode, len(plan.payload_bits), reason) + for mode in ("A", "B", "C", "D") + } + return CapturedPrbsAnalysis( + processing_success=False, + invalid_reason=reason, + sections=sections, + calibration=calibration, + continuity={ + "checked_boundary_count": 0, + "confirmed_discontinuity_count": 0, + "reason": "тоны недостоверны", + }, + modes=modes, + ) + + continuity = calibration_boundary_continuity( + np.asarray(samples, dtype=np.complex64), + callback_boundaries, + sections, + calibration, + actual_rx_rate, + ) + if continuity["confirmed_discontinuity_count"]: + reason = "в той же записи подтверждён фазовый разрыв" + modes = { + mode: lab043._failed_prbs_metrics(mode, len(plan.payload_bits), reason) + for mode in ("A", "B", "C", "D") + } + return CapturedPrbsAnalysis( + processing_success=False, + invalid_reason=reason, + sections=sections, + calibration=calibration, + continuity=continuity, + modes=modes, + ) + + modes = lab043.analyze_prbs_modes( + bpsk_samples, + plan.payload_bits, + calibration.carrier_offset_hz, + calibration.clock_scale, + actual_rx_rate, + known_pilot_bits=plan.pilot_bits, + ) + mode_d = modes["D"] + success = bool(mode_d.detected and mode_d.matched_bit_count == len(plan.payload_bits)) + reason = "" if success else (mode_d.failure_reason or "режим D не восстановил полную PRBS11") + return CapturedPrbsAnalysis( + processing_success=success, + invalid_reason=reason, + sections=sections, + calibration=calibration, + continuity=continuity, + modes=modes, + ) + + +def run_one_capture(label: str, output_directory: Path) -> dict: + if label != "S": + raise RuntimeError("На текущем этапе разрешён только один короткий опыт S") + duration = ( + lab043.PRBS_SHORT_DURATION_SECONDS + if label == "S" + else lab043.PRBS_LONG_DURATION_SECONDS + ) + plan = lab043.build_prbs11_transmission_plan(label, duration) + capture_started = datetime.now(timezone.utc) + + pluto_preflight = preflight_pluto() + rtl_preflight = preflight_rtlsdr() + time.sleep(0.25) + + pluto, tx_parameters = configure_pluto() + from rtlsdr import RtlSdr + + sdr = RtlSdr(device_index=0) + rx_parameters = _configure_rtlsdr(sdr) + actual_tx_rate = float(tx_parameters["actual_sample_rate_hz"]) + actual_rx_rate = float(rx_parameters["actual_sample_rate_hz"]) + expected_rx_samples = lab043.receive_sample_count( + len(plan.tx_samples), + actual_tx_rate, + actual_rx_rate, + ) + + ready = threading.Event() + tx_errors: list[str] = [] + tx_buffer_hold_seconds: list[float] = [] + + def transmit_once() -> None: + if not ready.wait(timeout=10.0): + tx_errors.append("RTL-SDR async-приём не подтвердил готовность") + return + time.sleep(lab043.RX_LEADING_MARGIN_SECONDS) + try: + hold_seconds = transmit_noncyclic_buffer_once( + pluto, + lab043.to_pluto_tx_samples(plan.tx_samples), + actual_tx_rate, + ) + tx_buffer_hold_seconds.append(hold_seconds) + except Exception as error: # pragma: no cover - аппаратный путь + tx_errors.append(f"{type(error).__name__}: {error}") + + tx_thread = threading.Thread(target=transmit_once, daemon=True) + tx_thread.start() + samples, boundaries, async_diagnostics = lab043.capture_continuous_rtlsdr_async( + sdr, + expected_rx_samples, + capture_ready_event=ready, + buffer_bytes=lab043.RTL_ASYNC_BUFFER_BYTES, + buffer_count=lab043.RTL_ASYNC_BUFFER_COUNT, + warmup_callback_count=lab043.RTL_ASYNC_WARMUP_CALLBACK_COUNT, + ) + tx_thread.join(timeout=15.0) + if tx_thread.is_alive(): + raise RuntimeError("Поток единственной передачи не завершился") + if tx_errors: + raise RuntimeError("; ".join(tx_errors)) + + received = np.asarray(samples, dtype=np.complex64) + if received.ndim != 1 or len(received) != expected_rx_samples: + raise RuntimeError("RX завершён, но размер массива не совпал с ожидаемым") + if not np.all(np.isfinite(received)): + raise RuntimeError("RX завершён, но массив содержит нечисловые значения") + boundary_rows = [asdict(row) for row in boundaries] + waveform_sha256 = hashlib.sha256( + np.asarray(plan.tx_samples, dtype=np.complex64).tobytes() + ).hexdigest() + tx_waveform_duration_seconds = len(plan.tx_samples) / actual_tx_rate + calibration_end_seconds = plan.calibration_sample_count / actual_tx_rate + guard_end_seconds = plan.bpsk_start_sample / actual_tx_rate + prbs_end_seconds = len(plan.tx_samples) / actual_tx_rate + overload = { + "peak_magnitude": float(np.max(np.abs(received))), + "clipped_component_fraction": float( + np.mean((np.abs(received.real) >= 0.999) | (np.abs(received.imag) >= 0.999)) + ), + } + metadata = { + "experiment": "Lab043 PRBS11 continuous async", + "capture_label": label, + "capture_started_utc": capture_started.isoformat().replace("+00:00", "Z"), + "capture_order": 1 if label == "S" else 2, + "capture_status": "captured", + "processing_status": "pending", + "processing_error": None, + "actual_sample_rate_hz": actual_rx_rate, + "center_frequency_hz": rx_parameters["actual_center_frequency_hz"], + "rx_gain_db": rx_parameters["actual_gain_db"], + "pluto_preflight": pluto_preflight, + "rtl_preflight": rtl_preflight, + "tx_parameters": tx_parameters, + "rx_parameters": rx_parameters, + "physical_scheme": "Pluto+ TX -> SMA -> AT30S 30 dB -> RTL-SDR RX; antennas removed", + "waveform_sha256": waveform_sha256, + "tx_waveform_duration_seconds": tx_waveform_duration_seconds, + "callback_count": int(async_diagnostics["callback_count"]), + "callback_boundaries": boundary_rows, + "expected_intervals": { + "relative_to_tx_start_seconds": { + "calibration": [0.0, calibration_end_seconds], + "guard": [calibration_end_seconds, guard_end_seconds], + "bpsk_marker_pilot_prbs": [guard_end_seconds, prbs_end_seconds], + }, + "rx_leading_margin_seconds": lab043.RX_LEADING_MARGIN_SECONDS, + "rx_trailing_margin_seconds": lab043.RX_TRAILING_MARGIN_SECONDS, + }, + "waveform": { + "calibration_tones_hz": [-lab043.F_CAL_HZ, lab043.F_CAL_HZ], + "calibration_duration_seconds": lab043.CALIBRATION_TONE_DURATION_SECONDS, + "fixed_guard_seconds": lab043.CALIBRATION_TO_BPSK_GUARD_SECONDS, + "single_initial_marker_bit_count": len(plan.marker_bits), + "pilot_symbol_count": len(plan.pilot_bits), + "pilot_duration_seconds": len(plan.pilot_bits) / lab043.SYMBOL_RATE, + "pilot_polynomial": lab043.PILOT_POLYNOMIAL, + "pilot_initial_state_hex": f"0x{lab043.PILOT_INITIAL_STATE:02X}", + "pilot_sha256": hashlib.sha256(plan.pilot_bits.tobytes()).hexdigest(), + "periodic_resynchronization": False, + "prbs_polynomial": lab043.PRBS11_POLYNOMIAL, + "prbs_initial_state_hex": f"0x{lab043.PRBS11_INITIAL_STATE:03X}", + "prbs_useful_duration_seconds": plan.useful_duration_seconds, + "transmitted_prbs_bit_count": len(plan.payload_bits), + "transmitted_prbs_sha256": hashlib.sha256(plan.payload_bits.tobytes()).hexdigest(), + "tx_sample_count": len(plan.tx_samples), + "tx_samples_sha256": waveform_sha256, + }, + "async_capture": async_diagnostics | { + "expected_sample_count": expected_rx_samples, + "callback_boundaries": boundary_rows, + }, + "tx_buffer_hold_seconds": ( + tx_buffer_hold_seconds[0] if tx_buffer_hold_seconds else None + ), + "overload": overload, + } + stamp = capture_started.strftime("%Y%m%d_%H%M%S") + base = output_directory / f"lab043_prbs11_{label.lower()}_{stamp}" + iq_path, json_path, analysis, processing_error = lab043.save_then_process_capture( + base, + received, + metadata, + lambda stored: analyze_captured_prbs( + stored, + plan, + actual_rx_rate, + actual_tx_rate, + boundaries, + ), + ) + if analysis is None: + return { + "iq_path": str(iq_path), + "metadata_path": str(json_path), + "capture_status": "captured", + "processing_status": "processing_failed", + "processing_error": processing_error, + "tx_parameters": tx_parameters, + "rx_parameters": rx_parameters, + "async_capture": async_diagnostics, + "overload": overload, + } + result = { + "iq_path": str(iq_path), + "metadata_path": str(json_path), + "capture_status": "captured", + "processing_status": "success", + "processing_error": None, + "processing_success": analysis.processing_success, + "invalid_reason": analysis.invalid_reason, + "calibration": asdict(analysis.calibration), + "continuity": analysis.continuity, + "modes": {mode: asdict(value) for mode, value in analysis.modes.items()}, + "tx_parameters": tx_parameters, + "rx_parameters": rx_parameters, + "async_capture": async_diagnostics, + "overload": overload, + } + return result + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--label", choices=("S",), required=True) + parser.add_argument("--output-directory", type=Path, default=Path("data/raw/lab043")) + arguments = parser.parse_args() + result = run_one_capture(arguments.label, arguments.output_directory) + print(json.dumps(result, ensure_ascii=False, indent=2)) + if result["processing_status"] != "success": + return 2 + mode_d = result["modes"]["D"] + return 0 if result["calibration"]["valid"] and mode_d["detected"] else 2 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/protocol/bpsk_radio.py b/protocol/bpsk_radio.py index 9f517c2..89a3692 100644 --- a/protocol/bpsk_radio.py +++ b/protocol/bpsk_radio.py @@ -17,6 +17,7 @@ Lab019 и Lab023 при импорте создают каталоги и сох from __future__ import annotations import struct +from dataclasses import dataclass import numpy as np @@ -38,8 +39,48 @@ CFO_REFINEMENT_STEP_HZ = 1.0 CFO_DEAD_ZONE_HZ = 15.0 MINIMUM_PHASE_CONSISTENCY = 0.55 +# Критерии рабочей оценки малой остаточной CFO по отдельному известному +# пилоту Lab043. Они зафиксированы до нового аппаратного опыта. +KNOWN_PILOT_BLOCK_SYMBOL_COUNT = 20 +MINIMUM_KNOWN_PILOT_BLOCK_COUNT = 8 +MINIMUM_KNOWN_PILOT_COHERENCE = 0.50 +MAXIMUM_KNOWN_PILOT_PHASE_FIT_RMSE_RAD = 0.25 + MAXIMUM_PROTOCOL_PACKET_BYTES = 4096 + +@dataclass(frozen=True) +class TwoTonePhaseRefinement: + """Явный результат фазового уточнения двух известных тонов.""" + + valid: bool + invalid_reason: str + low_frequency_hz: float + high_frequency_hz: float + carrier_offset_hz: float + tone_spacing_hz: float + residual_cfo_hz: float + phase_fit_rmse_rad: float + low_phase_fit_rms_rad: float + high_phase_fit_rms_rad: float + paired_phase_fit_rms_rad: float + block_count: int + + +@dataclass(frozen=True) +class KnownPilotCarrierEstimate: + """Оценка малой остаточной CFO по отдельному известному BPSK-пилоту.""" + + valid: bool + invalid_reason: str + frequency_hz: float + phase_increment_rad_per_symbol: float + initial_phase_rad: float + phase_fit_rmse_rad: float + mean_block_coherence: float + block_count: int + known_symbol_count: int + # Перенесено дословно из Lab023 (tests/lab023_guarded_cfo_correction.py). def bytes_to_bits( data: bytes, @@ -136,6 +177,496 @@ def bpsk_demodulate( ).astype(np.uint8) +def find_known_tone_peaks( + frequency_axis_hz: np.ndarray, + power_spectrum: np.ndarray, + expected_offset_hz: float, + search_half_width_hz: float, +) -> dict: + """Найти максимумы мощности около двух известных симметричных тонов. + + Функция не задаёт способ оценки шумового фона и порог достоверности: + эти решения принадлежат конкретному эксперименту. Если в одном из + поисковых окон нет конечных значений мощности, соответствующие частота + и мощность возвращаются как ``NaN``. + """ + + frequencies = np.asarray(frequency_axis_hz, dtype=np.float64) + powers = np.asarray(power_spectrum, dtype=np.float64) + + if frequencies.ndim != 1 or powers.ndim != 1: + raise ValueError("Ось частот и спектр мощности должны быть одномерными") + if len(frequencies) != len(powers): + raise ValueError("Ось частот и спектр мощности должны иметь одинаковую длину") + if not np.isfinite(expected_offset_hz) or expected_offset_hz <= 0.0: + raise ValueError("Ожидаемый отступ тона должен быть положительным") + if not np.isfinite(search_half_width_hz) or search_half_width_hz <= 0.0: + raise ValueError("Полуширина окна поиска должна быть положительной") + + def peak_in_window(center_hz: float) -> tuple[float, float]: + in_window = ( + (frequencies >= center_hz - search_half_width_hz) + & (frequencies <= center_hz + search_half_width_hz) + & np.isfinite(frequencies) + & np.isfinite(powers) + & (powers >= 0.0) + ) + indexes = np.flatnonzero(in_window) + if indexes.size == 0: + return float("nan"), float("nan") + peak_index = int(indexes[np.argmax(powers[indexes])]) + return float(frequencies[peak_index]), float(powers[peak_index]) + + low_frequency_hz, low_power = peak_in_window(-expected_offset_hz) + high_frequency_hz, high_power = peak_in_window(expected_offset_hz) + return { + "low_frequency_hz": low_frequency_hz, + "high_frequency_hz": high_frequency_hz, + "low_power": low_power, + "high_power": high_power, + } + + +def estimate_two_tone_offsets( + low_frequency_hz: float, + high_frequency_hz: float, + known_tone_offset_hz: float, +) -> dict: + """Оценить грубую CFO и отношение тактов по двум известным тонам. + + ``clock_scale`` определён как отношение масштаба TX к масштабу RX. + Поэтому для перехода принятой последовательности на временную сетку + передатчика её ожидаемая длина равна ``round(N * clock_scale)``. + Невычислимая оценка представляется только значениями ``NaN``. + """ + + if not np.isfinite(known_tone_offset_hz) or known_tone_offset_hz <= 0.0: + raise ValueError("Известный отступ тона должен быть положительным") + + invalid = { + "carrier_offset_hz": float("nan"), + "clock_scale": float("nan"), + "sample_clock_error_ppm": float("nan"), + } + if not np.isfinite(low_frequency_hz) or not np.isfinite(high_frequency_hz): + return invalid + if high_frequency_hz <= low_frequency_hz: + return invalid + + carrier_offset_hz = (low_frequency_hz + high_frequency_hz) / 2.0 + clock_scale = ( + (high_frequency_hz - low_frequency_hz) + / (2.0 * known_tone_offset_hz) + ) + if not np.isfinite(clock_scale) or clock_scale <= 0.0: + return invalid + + return { + "carrier_offset_hz": float(carrier_offset_hz), + "clock_scale": float(clock_scale), + "sample_clock_error_ppm": float((clock_scale - 1.0) * 1e6), + } + + +def _tone_block_projections( + samples: np.ndarray, + sample_rate_hz: float, + coarse_frequency_hz: float, + block_samples: int, + hop_samples: int, +) -> tuple[np.ndarray, np.ndarray]: + """Спроецировать фазово-непрерывный сигнал на грубую частоту тона.""" + + received = np.asarray(samples, dtype=np.complex128) + if received.ndim != 1: + raise ValueError("Комплексные отсчёты должны быть одномерными") + if not np.isfinite(sample_rate_hz) or sample_rate_hz <= 0.0: + raise ValueError("Частота дискретизации должна быть положительной") + if not np.isfinite(coarse_frequency_hz): + raise ValueError("Грубая частота тона должна быть конечной") + if not isinstance(block_samples, (int, np.integer)) or block_samples < 3: + raise ValueError("Размер блока должен быть целым числом не меньше трёх") + if not isinstance(hop_samples, (int, np.integer)) or hop_samples <= 0: + raise ValueError("Шаг блоков должен быть положительным целым числом") + if len(received) < block_samples + 2 * hop_samples: + raise ValueError("Для фазовой оценки нужны не менее трёх блоков") + + starts = np.arange( + 0, + len(received) - block_samples + 1, + hop_samples, + dtype=np.int64, + ) + sample_indexes = np.arange(len(received), dtype=np.float64) + mixed = received * np.exp( + -1j * 2.0 * np.pi * coarse_frequency_hz * sample_indexes / sample_rate_hz + ) + window = np.hanning(block_samples) + projections = np.asarray( + [ + np.sum(mixed[start : start + block_samples] * window) + for start in starts + ], + dtype=np.complex128, + ) + center_times_seconds = ( + starts.astype(np.float64) + (block_samples - 1) / 2.0 + ) / sample_rate_hz + return center_times_seconds, projections + + +def _weighted_phase_frequency( + center_times_seconds: np.ndarray, + projections: np.ndarray, +) -> tuple[float, float]: + """Оценить наклон развёрнутой фазы и СКО остатка линейной модели.""" + + times = np.asarray(center_times_seconds, dtype=np.float64) + values = np.asarray(projections, dtype=np.complex128) + weights = np.abs(values) ** 2 + weight_sum = float(np.sum(weights)) + if ( + times.ndim != 1 + or values.ndim != 1 + or len(times) != len(values) + or len(times) < 3 + or not np.all(np.isfinite(times)) + or not np.all(np.isfinite(values)) + or not np.isfinite(weight_sum) + or weight_sum <= 0.0 + ): + raise ValueError("Фазовые проекции не позволяют оценить частоту") + + phases = np.unwrap(np.angle(values)) + mean_time = float(np.sum(weights * times) / weight_sum) + mean_phase = float(np.sum(weights * phases) / weight_sum) + centered_times = times - mean_time + denominator = float(np.sum(weights * centered_times**2)) + if not np.isfinite(denominator) or denominator <= 0.0: + raise ValueError("Временные точки не позволяют оценить наклон фазы") + + slope_rad_per_second = float( + np.sum(weights * centered_times * (phases - mean_phase)) / denominator + ) + intercept = mean_phase - slope_rad_per_second * mean_time + residuals = phases - (intercept + slope_rad_per_second * times) + phase_fit_rms_rad = float( + np.sqrt(np.sum(weights * residuals**2) / weight_sum) + ) + return slope_rad_per_second / (2.0 * np.pi), phase_fit_rms_rad + + +def refine_tone_frequency_from_phase( + samples: np.ndarray, + sample_rate_hz: float, + coarse_frequency_hz: float, + block_samples: int, + hop_samples: int | None = None, +) -> dict: + """Уточнить частоту тона по наклону фазы когерентных проекций. + + Входной участок должен быть фазово-непрерывным. Функцию нельзя применять + через границы отдельных аппаратных чтений, если непрерывность их фазы не + доказана. Грубая частота должна быть достаточно точной, чтобы остаточная + фаза между соседними блоками разворачивалась без неоднозначности. + """ + + hop = block_samples // 2 if hop_samples is None else hop_samples + times, projections = _tone_block_projections( + samples, + sample_rate_hz, + coarse_frequency_hz, + block_samples, + hop, + ) + residual_frequency_hz, phase_fit_rms_rad = _weighted_phase_frequency( + times, + projections, + ) + return { + "frequency_hz": float(coarse_frequency_hz + residual_frequency_hz), + "residual_frequency_hz": float(residual_frequency_hz), + "phase_fit_rms_rad": float(phase_fit_rms_rad), + "block_count": int(len(projections)), + } + + +def refine_two_tone_frequencies_from_phase( + samples: np.ndarray, + sample_rate_hz: float, + low_coarse_frequency_hz: float, + high_coarse_frequency_hz: float, + block_samples: int, + hop_samples: int | None = None, +) -> TwoTonePhaseRefinement: + """Уточнить несущую и расстояние двух тонов на непрерывном участке. + + Несущая получается из двух индивидуальных наклонов фазы. Расстояние тонов + оценивается по фазе произведения верхней проекции на сопряжённую нижнюю: + общая фазовая ошибка приёмника при этом сокращается. Итоговые частоты + строятся из общей оценки центра и парной оценки расстояния. + """ + + if high_coarse_frequency_hz <= low_coarse_frequency_hz: + raise ValueError("Частота верхнего тона должна быть выше частоты нижнего") + hop = block_samples // 2 if hop_samples is None else hop_samples + times, low_projections = _tone_block_projections( + samples, + sample_rate_hz, + low_coarse_frequency_hz, + block_samples, + hop, + ) + high_times, high_projections = _tone_block_projections( + samples, + sample_rate_hz, + high_coarse_frequency_hz, + block_samples, + hop, + ) + if not np.array_equal(times, high_times): + raise RuntimeError("Временные сетки двух тонов не совпали") + + low_residual_hz, low_rms_rad = _weighted_phase_frequency(times, low_projections) + high_residual_hz, high_rms_rad = _weighted_phase_frequency(times, high_projections) + spacing_residual_hz, paired_rms_rad = _weighted_phase_frequency( + times, + high_projections * np.conj(low_projections), + ) + + individual_low_hz = low_coarse_frequency_hz + low_residual_hz + individual_high_hz = high_coarse_frequency_hz + high_residual_hz + center_hz = (individual_low_hz + individual_high_hz) / 2.0 + spacing_hz = ( + high_coarse_frequency_hz + - low_coarse_frequency_hz + + spacing_residual_hz + ) + low_frequency_hz = float(center_hz - spacing_hz / 2.0) + high_frequency_hz = float(center_hz + spacing_hz / 2.0) + residual_cfo_hz = float( + center_hz + - (low_coarse_frequency_hz + high_coarse_frequency_hz) / 2.0 + ) + phase_fit_rmse_rad = float( + np.sqrt(np.mean(np.square((low_rms_rad, high_rms_rad, paired_rms_rad)))) + ) + numeric = ( + low_frequency_hz, + high_frequency_hz, + center_hz, + spacing_hz, + residual_cfo_hz, + phase_fit_rmse_rad, + low_rms_rad, + high_rms_rad, + paired_rms_rad, + ) + valid = bool(all(np.isfinite(value) for value in numeric)) + if not valid: + nan = float("nan") + return TwoTonePhaseRefinement( + valid=False, + invalid_reason="фазовое уточнение вернуло невычислимое значение", + low_frequency_hz=nan, + high_frequency_hz=nan, + carrier_offset_hz=nan, + tone_spacing_hz=nan, + residual_cfo_hz=nan, + phase_fit_rmse_rad=nan, + low_phase_fit_rms_rad=nan, + high_phase_fit_rms_rad=nan, + paired_phase_fit_rms_rad=nan, + block_count=int(len(low_projections)), + ) + return TwoTonePhaseRefinement( + valid=True, + invalid_reason="", + low_frequency_hz=low_frequency_hz, + high_frequency_hz=high_frequency_hz, + carrier_offset_hz=float(center_hz), + tone_spacing_hz=float(spacing_hz), + residual_cfo_hz=residual_cfo_hz, + phase_fit_rmse_rad=phase_fit_rmse_rad, + low_phase_fit_rms_rad=float(low_rms_rad), + high_phase_fit_rms_rad=float(high_rms_rad), + paired_phase_fit_rms_rad=float(paired_rms_rad), + block_count=int(len(low_projections)), + ) + + +def apply_coarse_frequency_correction( + samples: np.ndarray, + carrier_offset_hz: float, + sample_rate_hz: float, +) -> np.ndarray: + """Убрать измеренный сдвиг несущей из комплексных отсчётов.""" + + received = np.asarray(samples, dtype=np.complex128) + if received.ndim != 1: + raise ValueError("Комплексные отсчёты должны быть одномерными") + if not np.isfinite(carrier_offset_hz): + raise ValueError("Сдвиг несущей должен быть конечным") + if not np.isfinite(sample_rate_hz) or sample_rate_hz <= 0.0: + raise ValueError("Частота дискретизации должна быть положительной") + + sample_indexes = np.arange(len(received), dtype=np.float64) + correction = np.exp( + -1j * 2.0 * np.pi * carrier_offset_hz * sample_indexes / sample_rate_hz + ) + return received * correction + + +def resample_for_clock_scale( + samples: np.ndarray, + clock_scale: float, +) -> np.ndarray: + """Передискретизировать отсчёты на временную сетку передатчика. + + Используется комплексная линейная интерполяция. Она детерминирована, + не требует рационального приближения ppm-отношения и подходит для + сильно передискретизированного BPSK-тракта. Направление преобразования + закрепляется синтетическими тестами для ошибок обоих знаков. + """ + + received = np.asarray(samples, dtype=np.complex128) + if received.ndim != 1: + raise ValueError("Комплексные отсчёты должны быть одномерными") + if not np.isfinite(clock_scale) or clock_scale <= 0.0: + raise ValueError("Масштаб такта должен быть положительным") + if received.size == 0: + return received.copy() + if received.size == 1: + return np.repeat(received, max(1, round(clock_scale))) + + output_length = max(1, round(len(received) * clock_scale)) + output_indexes = np.arange(output_length, dtype=np.float64) + source_positions = np.minimum(output_indexes / clock_scale, len(received) - 1.0) + source_indexes = np.arange(len(received), dtype=np.float64) + real = np.interp(source_positions, source_indexes, received.real) + imaginary = np.interp(source_positions, source_indexes, received.imag) + return real + 1j * imaginary + + +def estimate_known_pilot_carrier( + received_symbols: np.ndarray, + known_symbols: np.ndarray, + symbol_rate: float = SYMBOL_RATE, + block_symbol_count: int = KNOWN_PILOT_BLOCK_SYMBOL_COUNT, + minimum_block_coherence: float = MINIMUM_KNOWN_PILOT_COHERENCE, + maximum_phase_fit_rmse_rad: float = MAXIMUM_KNOWN_PILOT_PHASE_FIT_RMSE_RAD, +) -> KnownPilotCarrierEstimate: + """Оценить малую остаточную CFO по отдельному известному BPSK-пилоту. + + Известные знаки удаляются, затем символы объединяются в короткие + когерентные блоки. Развёрнутая фаза блоков аппроксимируется устойчивой + линейной моделью Huber. Полезная нагрузка для оценки не используется. + """ + + received = np.asarray(received_symbols, dtype=np.complex128) + known = np.asarray(known_symbols, dtype=np.complex128) + if received.ndim != 1 or known.ndim != 1: + raise ValueError("Принятый и известный пилоты должны быть одномерными") + if len(received) != len(known): + raise ValueError("Принятый и известный пилоты должны иметь одинаковую длину") + if not np.isfinite(symbol_rate) or symbol_rate <= 0.0: + raise ValueError("Символьная скорость должна быть положительной") + if not isinstance(block_symbol_count, (int, np.integer)) or block_symbol_count < 2: + raise ValueError("Размер когерентного блока должен быть целым и не меньше двух") + if not 0.0 < minimum_block_coherence <= 1.0: + raise ValueError("Порог когерентности должен находиться в интервале (0, 1]") + if not np.isfinite(maximum_phase_fit_rmse_rad) or maximum_phase_fit_rmse_rad <= 0.0: + raise ValueError("Порог ошибки фазовой модели должен быть положительным") + if not np.all(np.isfinite(received)) or not np.all(np.isfinite(known)): + raise ValueError("Пилот не должен содержать нечисловые значения") + if np.any(np.abs(known) <= 0.0): + raise ValueError("Все известные символы пилота должны быть ненулевыми") + + block_count = len(received) // block_symbol_count + if block_count < MINIMUM_KNOWN_PILOT_BLOCK_COUNT: + raise ValueError( + f"Для оценки нужны не менее {MINIMUM_KNOWN_PILOT_BLOCK_COUNT} когерентных блоков" + ) + usable_count = block_count * block_symbol_count + despread = received[:usable_count] * np.conj(known[:usable_count]) + blocks = despread.reshape(block_count, block_symbol_count) + block_sums = np.sum(blocks, axis=1) + block_magnitude_sums = np.sum(np.abs(blocks), axis=1) + 1e-12 + block_coherences = np.abs(block_sums) / block_magnitude_sums + mean_block_coherence = float(np.mean(block_coherences)) + + centers_symbols = ( + np.arange(block_count, dtype=np.float64) * block_symbol_count + + (block_symbol_count - 1) / 2.0 + ) + center_times_seconds = centers_symbols / symbol_rate + phases = np.unwrap(np.angle(block_sums)) + base_weights = np.maximum(np.abs(block_sums), 1e-12) + design = np.column_stack((np.ones(block_count), center_times_seconds)) + weights = base_weights.copy() + coefficients = np.zeros(2, dtype=np.float64) + + for _iteration in range(12): + root_weights = np.sqrt(weights) + new_coefficients, *_unused = np.linalg.lstsq( + design * root_weights[:, np.newaxis], + phases * root_weights, + rcond=None, + ) + residuals = phases - design @ new_coefficients + residual_median = float(np.median(residuals)) + robust_scale = ( + 1.4826 * float(np.median(np.abs(residuals - residual_median))) + + 1e-9 + ) + huber_limit = 1.345 * robust_scale + robust_weights = np.ones_like(residuals) + outliers = np.abs(residuals) > huber_limit + robust_weights[outliers] = huber_limit / np.abs(residuals[outliers]) + weights = base_weights * robust_weights + if np.allclose(coefficients, new_coefficients, rtol=0.0, atol=1e-12): + coefficients = new_coefficients + break + coefficients = new_coefficients + + residuals = phases - design @ coefficients + phase_fit_rmse_rad = float( + np.sqrt(np.sum(base_weights * residuals**2) / np.sum(base_weights)) + ) + frequency_hz = float(coefficients[1] / (2.0 * np.pi)) + initial_phase_rad = float(coefficients[0]) + phase_increment = float(2.0 * np.pi * frequency_hz / symbol_rate) + + invalid_reasons: list[str] = [] + if not all( + np.isfinite(value) + for value in (frequency_hz, initial_phase_rad, phase_fit_rmse_rad, mean_block_coherence) + ): + invalid_reasons.append("оценка содержит нечисловое значение") + if mean_block_coherence < minimum_block_coherence: + invalid_reasons.append( + f"когерентность {mean_block_coherence:.4f} ниже порога {minimum_block_coherence:.2f}" + ) + if phase_fit_rmse_rad > maximum_phase_fit_rmse_rad: + invalid_reasons.append( + f"RMSE фазы {phase_fit_rmse_rad:.4f} выше порога {maximum_phase_fit_rmse_rad:.2f} рад" + ) + + valid = not invalid_reasons + nan = float("nan") + return KnownPilotCarrierEstimate( + valid=valid, + invalid_reason="; ".join(invalid_reasons), + frequency_hz=frequency_hz if valid else nan, + phase_increment_rad_per_symbol=phase_increment if valid else nan, + initial_phase_rad=initial_phase_rad if valid else nan, + phase_fit_rmse_rad=phase_fit_rmse_rad, + mean_block_coherence=mean_block_coherence, + block_count=block_count, + known_symbol_count=len(known), + ) + + # Перенесено дословно из Lab023 (tests/lab023_guarded_cfo_correction.py). def estimate_carrier_parameters( received_symbols: np.ndarray, diff --git a/tests/test_bpsk_radio.py b/tests/test_bpsk_radio.py index 96c777b..ab74248 100644 --- a/tests/test_bpsk_radio.py +++ b/tests/test_bpsk_radio.py @@ -9,6 +9,7 @@ from __future__ import annotations import ast +import math import pathlib import struct @@ -210,3 +211,158 @@ def test_frame_search_finds_the_start_in_a_shaped_signal() -> None: recovered = found["symbol_samples"][start : start + len(bits)] assert radio.parse_radio_frame(radio.bpsk_demodulate(recovered)) == packet + + +# ------------------------------------------------------ Lab043: общие примитивы + + +def test_known_tone_peak_search_finds_both_windows() -> None: + frequencies = np.linspace(-100_000.0, 100_000.0, 4001) + powers = np.ones_like(frequencies) + powers[np.argmin(np.abs(frequencies + 51_200.0))] = 100.0 + powers[np.argmin(np.abs(frequencies - 49_300.0))] = 80.0 + + peaks = radio.find_known_tone_peaks(frequencies, powers, 50_000.0, 30_000.0) + + assert math.isclose(peaks["low_frequency_hz"], -51_200.0, abs_tol=1.0) + assert math.isclose(peaks["high_frequency_hz"], 49_300.0, abs_tol=1.0) + assert peaks["low_power"] == 100.0 + assert peaks["high_power"] == 80.0 + + +def test_known_tone_peak_search_returns_nan_without_bins() -> None: + frequencies = np.linspace(-1_000.0, 1_000.0, 101) + powers = np.ones_like(frequencies) + peaks = radio.find_known_tone_peaks(frequencies, powers, 50_000.0, 1_000.0) + + assert math.isnan(peaks["low_frequency_hz"]) + assert math.isnan(peaks["high_frequency_hz"]) + + +@pytest.mark.parametrize("ppm", [20.0, -20.0, 100.0, -100.0]) +def test_two_tone_clock_estimate_has_correct_sign_and_magnitude(ppm: float) -> None: + scale = 1.0 + ppm * 1e-6 + carrier_offset_hz = -1_250.0 + low = carrier_offset_hz - 50_000.0 * scale + high = carrier_offset_hz + 50_000.0 * scale + + estimate = radio.estimate_two_tone_offsets(low, high, 50_000.0) + + assert math.isclose(estimate["carrier_offset_hz"], carrier_offset_hz, abs_tol=1e-9) + assert math.isclose(estimate["clock_scale"], scale, abs_tol=1e-12) + assert math.isclose(estimate["sample_clock_error_ppm"], ppm, abs_tol=1e-6) + + +def test_two_tone_invalid_estimate_is_nan_not_zero() -> None: + estimate = radio.estimate_two_tone_offsets(float("nan"), 50_000.0, 50_000.0) + assert all(math.isnan(value) for value in estimate.values()) + + +@pytest.mark.parametrize("carrier_offset_hz", [730.0, -730.0]) +def test_coarse_frequency_correction_handles_both_signs(carrier_offset_hz: float) -> None: + sample_rate_hz = 20_000.0 + indexes = np.arange(20_000, dtype=np.float64) + impaired = np.exp(1j * 2.0 * np.pi * carrier_offset_hz * indexes / sample_rate_hz) + + corrected = radio.apply_coarse_frequency_correction( + impaired, + carrier_offset_hz, + sample_rate_hz, + ) + + assert np.max(np.abs(corrected - 1.0)) < 1e-9 + + +@pytest.mark.parametrize( + "carrier_offset_hz", + [-10.0, -5.0, -2.0, -1.2, -1.0, -0.5, 0.5, 1.0, 1.2, 2.0, 5.0, 10.0], +) +def test_known_pilot_estimator_resolves_sub_hertz_cfo_with_hardware_like_phase_noise( + carrier_offset_hz: float, +) -> None: + symbol_count = 1_280 + indexes = np.arange(symbol_count, dtype=np.float64) + known = np.where((indexes.astype(np.int64) * 73 + 19) % 127 < 64, 1.0, -1.0).astype( + np.complex128 + ) + estimates: list[float] = [] + valid_flags: list[bool] = [] + for repetition in range(96): + random_generator = np.random.default_rng( + 43_000_000 + + int(round((carrier_offset_hz + 20.0) * 1_000.0)) + + repetition + ) + phase_noise = random_generator.normal(0.0, 0.4691, symbol_count) + received = known * np.exp( + 1j + * ( + 2.0 * np.pi * carrier_offset_hz * indexes / radio.SYMBOL_RATE + + phase_noise + ) + ) + estimate = radio.estimate_known_pilot_carrier(received, known) + valid_flags.append(estimate.valid) + estimates.append(estimate.frequency_hz) + + errors = np.asarray(estimates) - carrier_offset_hz + assert all(valid_flags) + assert np.all(np.sign(estimates) == np.sign(carrier_offset_hz)) + assert abs(float(np.mean(errors))) < 0.08 + assert float(np.std(errors, ddof=1)) < 0.18 + assert float(np.percentile(np.abs(errors), 95.0)) < 0.35 + + +def test_known_pilot_estimator_rejects_noise_instead_of_reporting_false_cfo() -> None: + random_generator = np.random.default_rng(43_043) + known = np.resize(np.asarray([-1.0, 1.0], dtype=np.complex128), 1_280) + noise = ( + random_generator.normal(0.0, 1.0, len(known)) + + 1j * random_generator.normal(0.0, 1.0, len(known)) + ) + + estimate = radio.estimate_known_pilot_carrier(noise, known) + + assert not estimate.valid + assert estimate.invalid_reason + assert math.isnan(estimate.frequency_hz) + assert math.isnan(estimate.phase_increment_rad_per_symbol) + + +def _sample_at_positions(signal: np.ndarray, positions: np.ndarray) -> np.ndarray: + indexes = np.arange(len(signal), dtype=np.float64) + real = np.interp(positions, indexes, signal.real) + imaginary = np.interp(positions, indexes, signal.imag) + return real + 1j * imaginary + + +@pytest.mark.parametrize("ppm", [20.0, -20.0, 100.0, -100.0]) +def test_clock_resampling_direction_reduces_timing_error(ppm: float) -> None: + scale = 1.0 + ppm * 1e-6 + sample_count = 200_000 + indexes = np.arange(sample_count, dtype=np.float64) + reference = np.exp(1j * 2.0 * np.pi * 0.071 * indexes) + + received_length = math.floor(sample_count / scale) + received_positions = np.arange(received_length, dtype=np.float64) * scale + received = _sample_at_positions(reference, received_positions) + + corrected = radio.resample_for_clock_scale(received, scale) + wrong_direction = radio.resample_for_clock_scale(received, 1.0 / scale) + + uncorrected_count = min(len(received), len(reference)) + corrected_count = min(len(corrected), len(reference)) - 2 + wrong_count = min(len(wrong_direction), len(reference)) - 2 + uncorrected_error = float( + np.mean(np.abs(received[:uncorrected_count] - reference[:uncorrected_count]) ** 2) + ) + corrected_error = float( + np.mean(np.abs(corrected[:corrected_count] - reference[:corrected_count]) ** 2) + ) + wrong_error = float( + np.mean(np.abs(wrong_direction[:wrong_count] - reference[:wrong_count]) ** 2) + ) + + assert len(corrected) == round(len(received) * scale) + assert corrected_error < uncorrected_error + assert corrected_error < wrong_error diff --git a/tests/test_lab043_calibration.py b/tests/test_lab043_calibration.py new file mode 100644 index 0000000..7d738f3 --- /dev/null +++ b/tests/test_lab043_calibration.py @@ -0,0 +1,891 @@ +"""Быстрые синтетические проверки программной части Lab043.""" + +from __future__ import annotations + +import hashlib +import json +import math + +import numpy as np +import pytest +from scipy.signal import fftconvolve + +from experiments import lab043_pluto_to_rtlsdr as lab043 +from experiments import lab043_prbs11_hardware as lab043_hardware +from protocol import bpsk_radio as radio +from protocol.image_fragments import split_image_bytes +from protocol.packet import CRCError, MESSAGE_TYPE_TEXT, build_packet, parse_packet + + +def test_receive_duration_is_derived_from_actual_transmit_length() -> None: + assert lab043.receive_duration_seconds(4_800_000, 2_400_000.0) == 3.0 + assert lab043.receive_sample_count(4_800_000, 2_400_000.0, 2_399_900.0) == math.ceil( + 3.0 * 2_399_900.0 + ) + + +def test_continuous_blocks_are_joined_once_and_trimmed_at_end() -> None: + blocks = ( + np.asarray([0, 1, 2], dtype=np.complex64), + np.asarray([3, 4, 5], dtype=np.complex64), + ) + combined = lab043.combine_continuous_blocks(blocks, expected_sample_count=5) + assert np.array_equal(combined.real, np.arange(5)) + + +def test_short_continuous_capture_is_rejected() -> None: + with pytest.raises(ValueError, match="короче"): + lab043.combine_continuous_blocks((np.zeros(9),), expected_sample_count=10) + + +def test_async_collector_preserves_order_trims_target_and_records_boundaries() -> None: + collector = lab043.ContinuousAsyncIqCollector( + expected_sample_count=7, + warmup_callback_count=1, + ) + + assert not collector.add_callback(0, np.asarray([100, 101], dtype=complex)) + assert not collector.add_callback(1, np.asarray([0, 1, 2], dtype=complex)) + assert not collector.add_callback(2, np.asarray([3, 4, 5], dtype=complex)) + assert collector.add_callback(3, np.asarray([6, 7, 8], dtype=complex)) + + combined = collector.finalize() + assert np.array_equal(combined.real, np.arange(7)) + assert len(combined) == 7 + + boundaries = collector.callback_boundaries + assert [row.callback_index for row in boundaries] == [0, 1, 2, 3] + assert [row.accepted_start_sample for row in boundaries] == [0, 0, 3, 6] + assert [row.accepted_end_sample for row in boundaries] == [0, 3, 6, 7] + assert boundaries[0].warmup_discarded + assert boundaries[-1].trimmed_at_target + assert boundaries[-1].source_sample_count == 3 + assert boundaries[-1].accepted_sample_count == 1 + + +def test_async_collector_rejects_duplicate_callback() -> None: + collector = lab043.ContinuousAsyncIqCollector(4, warmup_callback_count=0) + collector.add_callback(0, np.asarray([0, 1], dtype=complex)) + with pytest.raises(ValueError, match="продублирован"): + collector.add_callback(0, np.asarray([2, 3], dtype=complex)) + + +def test_async_collector_rejects_missing_callback() -> None: + collector = lab043.ContinuousAsyncIqCollector(4, warmup_callback_count=0) + collector.add_callback(0, np.asarray([0, 1], dtype=complex)) + with pytest.raises(ValueError, match="пропуск"): + collector.add_callback(2, np.asarray([2, 3], dtype=complex)) + + +def test_async_collector_rejects_early_finalize_and_extra_callback() -> None: + collector = lab043.ContinuousAsyncIqCollector(2, warmup_callback_count=0) + collector.add_callback(0, np.asarray([0], dtype=complex)) + with pytest.raises(RuntimeError, match="короче"): + collector.finalize() + assert collector.add_callback(1, np.asarray([1], dtype=complex)) + with pytest.raises(RuntimeError, match="уже набрал"): + collector.add_callback(2, np.asarray([2], dtype=complex)) + + +def test_frame_concatenation_adds_no_lab043_gap() -> None: + first = lab043.shape_frame_like_lab042(np.ones(16, dtype=complex), samples_per_symbol=8) + second = lab043.shape_frame_like_lab042(-np.ones(12, dtype=complex), samples_per_symbol=8) + continuous = lab043.concatenate_frames_without_new_gaps((first, second)) + assert len(continuous) == len(first) + len(second) + assert np.array_equal(continuous[: len(first)], first) + assert np.array_equal(continuous[len(first) :], second) + + +def test_spectrum_method_finds_two_strong_tones_and_cfo() -> None: + sample_rate_hz = 2_400_000.0 + sample_count = 2 * lab043.CALIBRATION_NFFT + indexes = np.arange(sample_count, dtype=np.float64) + carrier_offset_hz = 1_250.0 + clock_scale = 1.0 + 100.0e-6 + low_hz = carrier_offset_hz - lab043.F_CAL_HZ * clock_scale + high_hz = carrier_offset_hz + lab043.F_CAL_HZ * clock_scale + random_generator = np.random.default_rng(43) + samples = ( + np.exp(1j * 2.0 * np.pi * low_hz * indexes / sample_rate_hz) + + 0.8 * np.exp(1j * 2.0 * np.pi * high_hz * indexes / sample_rate_hz) + + 0.002 + * ( + random_generator.standard_normal(sample_count) + + 1j * random_generator.standard_normal(sample_count) + ) + ) + + estimate = lab043.estimate_calibration(samples, sample_rate_hz) + + assert estimate.valid + assert estimate.low_peak_excess_db > 10.0 + assert estimate.high_peak_excess_db > 10.0 + assert math.isclose(estimate.carrier_offset_hz, carrier_offset_hz, abs_tol=5.0) + assert math.isclose(estimate.sample_clock_error_ppm, 100.0, abs_tol=60.0) + + +def _synthetic_two_tones( + low_hz: float, + high_hz: float, + seed: int, +) -> tuple[np.ndarray, float]: + sample_rate_hz = 240_000.0 + sample_count = 120_000 + indexes = np.arange(sample_count, dtype=np.float64) + random_generator = np.random.default_rng(seed) + samples = ( + np.exp(1j * 2.0 * np.pi * low_hz * indexes / sample_rate_hz) + + 0.8 * np.exp(1j * 2.0 * np.pi * high_hz * indexes / sample_rate_hz) + + 0.003 + * ( + random_generator.standard_normal(sample_count) + + 1j * random_generator.standard_normal(sample_count) + ) + ) + return samples, sample_rate_hz + + +def _phase_refinement_from_spectrum( + samples: np.ndarray, + sample_rate_hz: float, +) -> tuple[lab043.CalibrationResult, radio.TwoTonePhaseRefinement]: + spectrum_estimate = lab043.estimate_calibration(samples, sample_rate_hz) + assert spectrum_estimate.valid + phase_estimate = radio.refine_two_tone_frequencies_from_phase( + samples, + sample_rate_hz, + spectrum_estimate.f_low_hz, + spectrum_estimate.f_high_hz, + block_samples=240, + hop_samples=120, + ) + return spectrum_estimate, phase_estimate + + +@pytest.mark.parametrize( + "tone_shift_hz", + [-1.0, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1.0], +) +def test_phase_refinement_improves_known_sub_hertz_common_shift( + tone_shift_hz: float, +) -> None: + true_low_hz = -lab043.F_CAL_HZ + tone_shift_hz + true_high_hz = lab043.F_CAL_HZ + tone_shift_hz + samples, sample_rate_hz = _synthetic_two_tones( + true_low_hz, + true_high_hz, + seed=43_000 + int((tone_shift_hz + 2.0) * 100), + ) + + old, refined = _phase_refinement_from_spectrum(samples, sample_rate_hz) + old_max_error_hz = max( + abs(old.f_low_hz - true_low_hz), + abs(old.f_high_hz - true_high_hz), + ) + refined_max_error_hz = max( + abs(refined.low_frequency_hz - true_low_hz), + abs(refined.high_frequency_hz - true_high_hz), + ) + + assert refined_max_error_hz < 0.001 + assert refined_max_error_hz < old_max_error_hz + + +@pytest.mark.parametrize( + "half_spacing_shift_hz", + [-1.0, -0.5, -0.25, -0.1, 0.1, 0.25, 0.5, 1.0], +) +def test_phase_refinement_improves_known_sub_hertz_spacing_shift( + half_spacing_shift_hz: float, +) -> None: + true_low_hz = -lab043.F_CAL_HZ - half_spacing_shift_hz + true_high_hz = lab043.F_CAL_HZ + half_spacing_shift_hz + true_offsets = radio.estimate_two_tone_offsets( + true_low_hz, + true_high_hz, + lab043.F_CAL_HZ, + ) + samples, sample_rate_hz = _synthetic_two_tones( + true_low_hz, + true_high_hz, + seed=44_000 + int((half_spacing_shift_hz + 2.0) * 100), + ) + + old, refined = _phase_refinement_from_spectrum(samples, sample_rate_hz) + refined_offsets = radio.estimate_two_tone_offsets( + refined.low_frequency_hz, + refined.high_frequency_hz, + lab043.F_CAL_HZ, + ) + old_error_ppm = abs( + old.sample_clock_error_ppm - true_offsets["sample_clock_error_ppm"] + ) + refined_error_ppm = abs( + refined_offsets["sample_clock_error_ppm"] + - true_offsets["sample_clock_error_ppm"] + ) + + assert refined_error_ppm < 0.001 + assert refined_error_ppm < old_error_ppm + + +@pytest.mark.parametrize("sample_clock_error_ppm", [100.0, -100.0]) +def test_public_refined_calibration_returns_complete_contract( + sample_clock_error_ppm: float, +) -> None: + carrier_offset_hz = 1_234.5 + scale = 1.0 + sample_clock_error_ppm * 1e-6 + samples, sample_rate_hz = _synthetic_two_tones( + carrier_offset_hz - lab043.F_CAL_HZ * scale, + carrier_offset_hz + lab043.F_CAL_HZ * scale, + seed=45_000 + int(sample_clock_error_ppm), + ) + + result = lab043.estimate_refined_calibration(samples, sample_rate_hz) + + assert isinstance(result, lab043.CalibrationResult) + assert result.valid + assert result.invalid_reason == "" + assert math.isclose(result.carrier_offset_hz, carrier_offset_hz, abs_tol=0.01) + assert math.isclose( + result.sample_clock_error_ppm, + sample_clock_error_ppm, + abs_tol=0.01, + ) + assert math.isfinite(result.residual_cfo_hz) + assert math.isfinite(result.phase_fit_rmse_rad) + assert result.peak_margin_low_db >= lab043.MINIMUM_TONE_EXCESS_DB + assert result.peak_margin_high_db >= lab043.MINIMUM_TONE_EXCESS_DB + + +def test_public_refined_calibration_rejects_weak_or_missing_tone_without_exception() -> None: + sample_rate_hz = 240_000.0 + sample_count = 120_000 + indexes = np.arange(sample_count, dtype=np.float64) + random_generator = np.random.default_rng(46_000) + noise = 0.003 * ( + random_generator.standard_normal(sample_count) + + 1j * random_generator.standard_normal(sample_count) + ) + one_tone = np.exp(1j * 2.0 * np.pi * lab043.F_CAL_HZ * indexes / sample_rate_hz) + + weak = lab043.estimate_refined_calibration(noise, sample_rate_hz) + missing = lab043.estimate_refined_calibration(one_tone + noise, sample_rate_hz) + + for result in (weak, missing): + assert isinstance(result, lab043.CalibrationResult) + assert not result.valid + assert result.invalid_reason + assert math.isnan(result.f_low_hz) + assert math.isnan(result.f_high_hz) + assert math.isnan(result.carrier_offset_hz) + assert math.isnan(result.sample_clock_error_ppm) + + +def test_tones_below_ten_decibels_are_rejected_with_nan() -> None: + frequencies = np.linspace(-100_000.0, 100_000.0, 4001) + powers = np.ones_like(frequencies) + powers[np.argmin(np.abs(frequencies + lab043.F_CAL_HZ))] = 9.0 + powers[np.argmin(np.abs(frequencies - lab043.F_CAL_HZ))] = 9.0 + + estimate = lab043.estimate_calibration_from_spectrum(frequencies, powers) + + assert not estimate.valid + assert math.isnan(estimate.f_low_hz) + assert math.isnan(estimate.f_high_hz) + assert math.isnan(estimate.carrier_offset_hz) + assert estimate.failure_reason + + +def _estimate(carrier_offset_hz: float, valid: bool = True) -> lab043.CalibrationResult: + if not valid: + return lab043._invalid_calibration( + "нет тонов", + noise_power=1.0, + ) + return lab043.CalibrationResult( + valid=True, + invalid_reason="", + f_low_hz=-50_000.0 + carrier_offset_hz, + f_high_hz=50_000.0 + carrier_offset_hz, + carrier_offset_hz=carrier_offset_hz, + carrier_offset_ppm=carrier_offset_hz / lab043.CARRIER_HZ * 1e6, + clock_scale=1.0, + sample_clock_error_ppm=0.0, + residual_cfo_hz=0.0, + phase_fit_rmse_rad=0.0, + peak_margin_low_db=20.0, + peak_margin_high_db=21.0, + noise_power=1.0, + ) + + +def test_three_calibrations_keep_individual_failures_and_finite_statistics() -> None: + summary = lab043.summarize_calibrations((_estimate(10.0), _estimate(14.0), _estimate(0, False))) + + assert summary["capture_count"] == 3 + assert summary["failure_count"] == 1 + carrier = summary["carrier_offset_hz"] + assert carrier["mean"] == 12.0 + assert carrier["minimum"] == 10.0 + assert carrier["maximum"] == 14.0 + assert carrier["standard_deviation"] == 2.0 + + +def test_all_invalid_calibrations_keep_nan_aggregates() -> None: + summary = lab043.summarize_calibrations((_estimate(0, False),) * 3) + assert summary["failure_count"] == 3 + assert math.isnan(summary["carrier_offset_hz"]["mean"]) + + +def test_calibration_summary_requires_exactly_three_captures() -> None: + with pytest.raises(ValueError, match="ровно три"): + lab043.summarize_calibrations((_estimate(10.0), _estimate(12.0))) + + +def test_prbs11_has_full_2047_bit_period() -> None: + sequence = lab043.prbs11(2 * lab043.PRBS11_LENGTH) + first = sequence[: lab043.PRBS11_LENGTH] + second = sequence[lab043.PRBS11_LENGTH :] + + assert np.array_equal(first, second) + assert not np.array_equal(first, np.roll(first, 23)) + assert not np.array_equal(first, np.roll(first, 89)) + + +def test_separate_pilot_is_fixed_balanced_and_independent_from_prbs11() -> None: + pilot = lab043.pilot_prbs7() + first_period = pilot[:127] + second_period = pilot[127:254] + + assert len(pilot) == 1_280 + assert np.array_equal(first_period, second_period) + assert int(np.count_nonzero(first_period)) == 64 + assert not np.array_equal(pilot, lab043.prbs11(len(pilot))) + longest_run = max( + len(group) + for group in np.split(pilot, np.flatnonzero(np.diff(pilot)) + 1) + ) + assert longest_run <= 7 + + +@pytest.mark.parametrize( + "carrier_offset_hz", + [-10.0, -5.0, -2.0, -1.2, -1.0, -0.5, 0.5, 1.0, 1.2, 2.0, 5.0, 10.0], +) +def test_fixed_pilot_estimator_statistics_match_hardware_s_phase_dispersion( + carrier_offset_hz: float, +) -> None: + known = radio.bpsk_modulate(lab043.pilot_prbs7()) + indexes = np.arange(len(known), dtype=np.float64) + estimates: list[float] = [] + for repetition in range(96): + random_generator = np.random.default_rng( + 43_043_000 + + int(round((carrier_offset_hz + 20.0) * 1_000.0)) + + repetition + ) + phase_noise = random_generator.normal(0.0, 0.4691, len(known)) + received = known * np.exp( + 1j + * ( + 2.0 * np.pi * carrier_offset_hz * indexes / lab043.SYMBOL_RATE + + phase_noise + ) + ) + estimate = radio.estimate_known_pilot_carrier(received, known) + assert estimate.valid + estimates.append(estimate.frequency_hz) + + errors = np.asarray(estimates) - carrier_offset_hz + assert np.all(np.sign(estimates) == np.sign(carrier_offset_hz)) + assert abs(float(np.mean(errors))) < 0.08 + assert float(np.std(errors, ddof=1)) < 0.18 + assert float(np.percentile(np.abs(errors), 95.0)) < 0.35 + + +def test_fixed_pilot_estimator_rejects_noise_only_input() -> None: + random_generator = np.random.default_rng(43_043) + known = radio.bpsk_modulate(lab043.pilot_prbs7()) + noise = ( + random_generator.normal(0.0, 1.0, len(known)) + + 1j * random_generator.normal(0.0, 1.0, len(known)) + ) + + estimate = radio.estimate_known_pilot_carrier(noise, known) + + assert not estimate.valid + assert math.isnan(estimate.frequency_hz) + + +@pytest.mark.parametrize( + ("label", "duration_seconds", "expected_bits"), + [("S", 0.25, 5_000), ("L", 2.0, 40_000)], +) +def test_prbs_transmission_plan_is_fixed_and_has_one_initial_marker( + label: str, + duration_seconds: float, + expected_bits: int, +) -> None: + plan = lab043.build_prbs11_transmission_plan(label, duration_seconds) + + assert len(plan.payload_bits) == expected_bits + assert np.array_equal( + plan.payload_bits, + lab043.prbs11(expected_bits, lab043.PRBS11_INITIAL_STATE), + ) + assert plan.bpsk_start_sample == ( + plan.calibration_sample_count + plan.fixed_guard_sample_count + ) + assert plan.calibration_sample_count == 600_000 + assert plan.fixed_guard_sample_count == 240_000 + assert np.max(np.abs(plan.tx_samples)) <= 1.0 + if label == "S": + assert len(plan.pilot_bits) == lab043.PILOT_SYMBOL_COUNT + assert len(plan.pilot_bits) / lab043.SYMBOL_RATE == 0.064 + else: + assert len(plan.pilot_bits) == 0 + + +def test_noncyclic_tx_buffer_is_held_for_full_waveform_duration() -> None: + events: list[object] = [] + + class FakePluto: + def tx_destroy_buffer(self) -> None: + events.append("destroy") + + def tx(self, samples: np.ndarray) -> None: + events.append(("tx", samples.copy())) + + samples = np.arange(12, dtype=np.float32).astype(np.complex64) + duration = lab043_hardware.transmit_noncyclic_buffer_once( + FakePluto(), + samples, + actual_sample_rate_hz=48.0, + sleep_function=lambda seconds: events.append(("sleep", seconds)), + ) + + assert duration == 0.25 + assert events[0] == "destroy" + assert events[1][0] == "tx" + assert np.array_equal(events[1][1], samples) + assert events[2] == ("sleep", 0.25) + assert events[3] == "destroy" + + +def _synthetic_prbs_channel( + transmitted: np.ndarray, + sample_rate_hz: float, + carrier_offset_hz: float, + clock_scale: float, +) -> np.ndarray: + received_length = max(1, math.floor(len(transmitted) / clock_scale)) + receive_indexes = np.arange(received_length, dtype=np.float64) + source_positions = np.minimum(receive_indexes * clock_scale, len(transmitted) - 1.0) + source_indexes = np.arange(len(transmitted), dtype=np.float64) + received = np.interp(source_positions, source_indexes, transmitted.real) + 1j * np.interp( + source_positions, + source_indexes, + transmitted.imag, + ) + return received * np.exp( + 1j * 2.0 * np.pi * carrier_offset_hz * receive_indexes / sample_rate_hz + ) + + +def test_prbs_modes_measure_drift_and_resampling_reduces_it() -> None: + samples_per_symbol = 16 + sample_rate_hz = lab043.SYMBOL_RATE * samples_per_symbol + clock_scale = 1.0 - 100.0e-6 + carrier_offset_hz = 730.0 + plan = lab043.build_prbs11_transmission_plan( + "S", + 0.25, + sample_rate_hz=sample_rate_hz, + samples_per_symbol=samples_per_symbol, + ) + transmitted_bpsk = plan.tx_samples[plan.bpsk_start_sample :] + received = _synthetic_prbs_channel( + transmitted_bpsk, + sample_rate_hz, + carrier_offset_hz, + clock_scale, + ) + + modes = lab043.analyze_prbs_modes( + received, + plan.payload_bits, + carrier_offset_hz, + clock_scale, + sample_rate_hz, + samples_per_symbol=samples_per_symbol, + ) + + assert modes["B"].detected + assert modes["C"].detected + assert modes["D"].detected + assert modes["B"].matched_bit_count == len(plan.payload_bits) + assert modes["D"].bit_error_rate == 0.0 + assert modes["B"].timing_sro_ppm < 0.0 + assert math.isclose(modes["B"].timing_sro_ppm, -100.0, abs_tol=40.0) + assert abs(modes["D"].accumulated_timing_drift_samples) < abs( + modes["B"].accumulated_timing_drift_samples + ) + + +def test_marker_pilot_correction_and_unknown_prbs_follow_one_software_path() -> None: + samples_per_symbol = 16 + sample_rate_hz = lab043.SYMBOL_RATE * samples_per_symbol + coarse_cfo_hz = 700.0 + residual_cfo_hz = 1.2 + plan = lab043.build_prbs11_transmission_plan( + "S", + 0.25, + sample_rate_hz=sample_rate_hz, + samples_per_symbol=samples_per_symbol, + ) + transmitted_bpsk = plan.tx_samples[plan.bpsk_start_sample :] + received = _synthetic_prbs_channel( + transmitted_bpsk, + sample_rate_hz, + coarse_cfo_hz + residual_cfo_hz, + 1.0, + ) + + modes = lab043.analyze_prbs_modes( + received, + plan.payload_bits, + coarse_cfo_hz, + 1.0, + sample_rate_hz, + samples_per_symbol=samples_per_symbol, + known_pilot_bits=plan.pilot_bits, + ) + wrong_reference_modes = lab043.analyze_prbs_modes( + received, + np.zeros_like(plan.payload_bits), + coarse_cfo_hz, + 1.0, + sample_rate_hz, + samples_per_symbol=samples_per_symbol, + known_pilot_bits=plan.pilot_bits, + ) + + assert set(modes) == {"A", "B", "C", "D"} + assert modes["C"].pilot_estimate_valid + assert modes["C"].pilot_cfo_applied + assert math.isclose( + modes["C"].estimated_pilot_cfo_hz, + residual_cfo_hz, + abs_tol=0.25, + ) + assert abs(modes["C"].residual_pilot_cfo_hz) < 0.05 + assert modes["C"].bit_error_rate < modes["B"].bit_error_rate + assert modes["C"].bit_error_rate == 0.0 + assert math.isclose( + wrong_reference_modes["C"].estimated_pilot_cfo_hz, + modes["C"].estimated_pilot_cfo_hz, + abs_tol=1e-12, + ) + assert wrong_reference_modes["C"].bit_error_rate > 0.4 + + +def test_calibration_and_bpsk_are_split_from_the_same_capture() -> None: + samples_per_symbol = 16 + sample_rate_hz = lab043.SYMBOL_RATE * samples_per_symbol + plan = lab043.build_prbs11_transmission_plan( + "S", + 0.025, + sample_rate_hz=sample_rate_hz, + samples_per_symbol=samples_per_symbol, + ) + leading = np.zeros(round(0.5 * sample_rate_hz), dtype=np.complex128) + trailing = np.zeros(round(0.5 * sample_rate_hz), dtype=np.complex128) + capture = np.concatenate((leading, plan.tx_samples, trailing)) + + calibration, bpsk, sections = lab043.split_calibration_and_bpsk( + capture, + plan, + sample_rate_hz, + sample_rate_hz, + ) + + assert len(calibration) > 0.15 * sample_rate_hz + assert len(bpsk) > len(plan.tx_samples) - plan.bpsk_start_sample + assert sections["calibration_end_sample"] < sections["bpsk_nominal_start_sample"] + + +def _continuous_synthetic_capture( + plan: lab043.PrbsTransmissionPlan, + tx_samples: np.ndarray | None = None, +) -> np.ndarray: + leading = np.zeros( + round(lab043.RX_LEADING_MARGIN_SECONDS * lab043.SAMPLE_RATE_HZ), + dtype=np.complex64, + ) + trailing = np.zeros( + round(lab043.RX_TRAILING_MARGIN_SECONDS * lab043.SAMPLE_RATE_HZ), + dtype=np.complex64, + ) + waveform = plan.tx_samples if tx_samples is None else tx_samples + return np.concatenate((leading, waveform.astype(np.complex64), trailing)) + + +def test_full_short_s_processing_path_returns_calibration_and_modes() -> None: + plan = lab043.build_prbs11_transmission_plan("S", 0.025) + capture = _continuous_synthetic_capture(plan) + + result = lab043_hardware.analyze_captured_prbs( + capture, + plan, + lab043.SAMPLE_RATE_HZ, + lab043.SAMPLE_RATE_HZ, + ) + + assert result.processing_success + assert result.calibration.valid + assert set(result.modes) == {"A", "B", "C", "D"} + assert result.modes["D"].detected + assert result.modes["D"].matched_bit_count == len(plan.payload_bits) + assert result.modes["D"].bit_error_rate == 0.0 + assert math.isfinite(result.modes["D"].evm_percent) + + +def test_full_short_s_missing_prbs_is_not_reported_as_zero_ber() -> None: + plan = lab043.build_prbs11_transmission_plan("S", 0.025) + without_prbs = plan.tx_samples.copy() + without_prbs[plan.bpsk_start_sample :] = 0.0 + capture = _continuous_synthetic_capture(plan, without_prbs) + + result = lab043_hardware.analyze_captured_prbs( + capture, + plan, + lab043.SAMPLE_RATE_HZ, + lab043.SAMPLE_RATE_HZ, + ) + + assert result.calibration.valid + assert not result.processing_success + assert not result.modes["D"].detected + assert math.isnan(result.modes["D"].bit_error_rate) + + +def test_full_short_s_invalid_calibration_makes_all_modes_na() -> None: + plan = lab043.build_prbs11_transmission_plan("S", 0.025) + one_tone_waveform = plan.tx_samples.copy() + indexes = np.arange(plan.calibration_sample_count, dtype=np.float64) + one_tone_waveform[: plan.calibration_sample_count] = 0.35 * np.exp( + 1j * 2.0 * np.pi * lab043.F_CAL_HZ * indexes / lab043.SAMPLE_RATE_HZ + ) + capture = _continuous_synthetic_capture(plan, one_tone_waveform) + + result = lab043_hardware.analyze_captured_prbs( + capture, + plan, + lab043.SAMPLE_RATE_HZ, + lab043.SAMPLE_RATE_HZ, + ) + + assert not result.calibration.valid + assert result.calibration.invalid_reason + assert not result.processing_success + for mode in result.modes.values(): + assert not mode.detected + assert math.isnan(mode.bit_error_rate) + + +def test_raw_survives_processing_exception_and_json_records_failure(tmp_path) -> None: + samples = np.asarray([0.25 + 0.5j, -0.5 + 0.25j], dtype=np.complex64) + metadata = { + "capture_status": "captured", + "processing_status": "pending", + "processing_error": None, + "capture_started_utc": "2026-08-19T12:00:00Z", + "tx_parameters": {"gain_db": -30.0}, + "rx_parameters": {"gain_db": None}, + "waveform_sha256": "a" * 64, + "callback_count": 2, + "callback_boundaries": [{"callback_index": 0}, {"callback_index": 1}], + "tx_waveform_duration_seconds": 0.25, + "expected_intervals": {"calibration": [0.0, 0.1]}, + } + + iq_path, json_path, result, error = lab043.save_then_process_capture( + tmp_path / "failed_after_rx", + samples, + metadata, + lambda _samples: (_ for _ in ()).throw(RuntimeError("synthetic failure")), + ) + + document = json.loads(json_path.read_text(encoding="utf-8")) + assert iq_path.exists() + assert json_path.exists() + assert result is None + assert error == "RuntimeError: synthetic failure" + assert document["capture_status"] == "captured" + assert document["processing_status"] == "processing_failed" + assert document["processing_error"] == error + assert document["raw_verified"] + assert document["iq_sha256"] == hashlib.sha256(iq_path.read_bytes()).hexdigest() + + +def test_ber_is_measured_before_crc() -> None: + expected = lab043.prbs11() + received = expected.copy() + received[[0, 100, 1000]] ^= 1 + + errors, bit_error_rate = lab043.bit_error_rate(expected, received) + + assert errors == 3 + assert bit_error_rate == 3 / lab043.PRBS11_LENGTH + + +def test_ber_is_nan_for_truncated_known_sequence() -> None: + expected = lab043.prbs11() + errors, bit_error_rate = lab043.bit_error_rate(expected, expected[:-1]) + assert errors == 0 + assert math.isnan(bit_error_rate) + + +@pytest.mark.parametrize("initial_sample_phase", [0, 3, 7, 15]) +def test_symbol_sample_phase_is_selected_separately(initial_sample_phase: int) -> None: + samples_per_symbol = 16 + _, marker_symbols = radio.build_frame_marker() + taps = radio.root_raised_cosine_taps(0.35, samples_per_symbol, 10) + upsampled = np.zeros(len(marker_symbols) * samples_per_symbol, dtype=complex) + upsampled[::samples_per_symbol] = marker_symbols + transmitted = fftconvolve(upsampled, taps, mode="full") + shifted = np.concatenate((np.zeros(initial_sample_phase), transmitted)) + matched = fftconvolve(shifted, taps, mode="full") + + found = radio.find_radio_frame(matched, marker_symbols, samples_per_symbol) + + assert found["sample_phase"] == initial_sample_phase + assert found["score"] > 0.99 + + +def test_one_control_packet_passes_crc() -> None: + assert lab043.control_packet_crc_roundtrip(b"known payload") + + +def test_corrupted_control_packet_fails_crc() -> None: + packet = bytearray(build_packet(b"known payload", MESSAGE_TYPE_TEXT, 43)) + packet[-1] ^= 1 + bits, _, _ = radio.build_radio_frame(bytes(packet)) + recovered = radio.parse_radio_frame(bits) + assert recovered is not None + with pytest.raises(CRCError): + parse_packet(recovered) + + +def test_modes_a_b_c_process_the_same_capture_and_c_recovers_prbs() -> None: + clock_scale = 1.0 + 100.0e-6 + capture, expected = lab043.synthesize_known_bpsk_capture( + carrier_offset_hz=730.0, + clock_scale=clock_scale, + samples_per_symbol=16, + ) + + results = lab043.process_bpsk_modes( + capture, + expected, + coarse_carrier_offset_hz=700.0, + clock_scale=clock_scale, + samples_per_symbol=16, + ) + + assert set(results) == {"A", "B", "C"} + assert results["C"].bit_error_rate == 0.0 + assert results["C"].fine_cfo_applied + assert not math.isfinite(results["A"].bit_error_rate) or ( + results["A"].bit_error_rate > results["C"].bit_error_rate + ) + assert results["B"].bit_error_rate > results["C"].bit_error_rate + + +def test_missing_jpeg_fragment_never_produces_an_image() -> None: + fragments = split_image_bytes(bytes(range(250)) * 5, image_id=43, fragment_data_size=128) + assert len(fragments) == 10 + assert lab043.try_reassemble_complete_image(fragments[:-1]) is None + + +def _complete_acceptance() -> lab043.AcceptanceResult: + return lab043.AcceptanceResult(10, 10, 10, 10, True, True, True, True) + + +@pytest.mark.parametrize( + "result", + [ + lab043.AcceptanceResult(9, 10, 10, 10, True, True, True, True), + lab043.AcceptanceResult(10, 10, 9, 10, True, True, True, True), + lab043.AcceptanceResult(10, 10, 10, 9, False, False, False, True), + lab043.AcceptanceResult(10, 10, 10, 10, True, True, True, False), + ], +) +def test_exit_zero_requires_every_acceptance_layer(result: lab043.AcceptanceResult) -> None: + assert lab043.exit_code_for_acceptance(result) == 1 + + +def test_exit_zero_is_allowed_for_complete_acceptance() -> None: + result = _complete_acceptance() + assert result.radio_passed + assert result.application_passed + assert lab043.exit_code_for_acceptance(result) == 0 + + +def test_reference_iq_format_contains_reprocessing_metadata(tmp_path) -> None: + metadata = { + "actual_sample_rate_hz": 2_400_000.0, + "center_frequency_hz": 435_000_000.0, + "rx_gain_db": 19.7, + "capture_started_utc": "2026-08-19T12:00:00Z", + "capture_order": 1, + "tx_parameters": {"gain_db": -10.0}, + "calibration": {"carrier_offset_hz": 123.0, "clock_scale": 1.00002}, + } + samples = np.asarray([1 + 2j, 3 + 4j], dtype=np.complex64) + + iq_path, metadata_path = lab043.save_reference_iq_capture( + tmp_path / "reference", + samples, + metadata, + ) + + restored = np.load(iq_path, allow_pickle=False) + document = json.loads(metadata_path.read_text(encoding="utf-8")) + assert np.array_equal(restored, samples) + assert document["sample_count"] == 2 + assert document["dtype"] == "complex64" + assert len(document["iq_sha256"]) == 64 + assert document["calibration"]["clock_scale"] == 1.00002 + + +def test_reference_iq_rejects_incomplete_metadata(tmp_path) -> None: + with pytest.raises(ValueError, match="Не хватает метаданных"): + lab043.save_reference_iq_capture(tmp_path / "reference", np.zeros(4), {}) + + +def test_diagnostic_csv_txt_and_png_are_created_from_supplied_results(tmp_path) -> None: + calibrations = (_estimate(10.0), _estimate(12.0), _estimate(14.0)) + capture, expected = lab043.synthesize_known_bpsk_capture(30.0, 1.0, 16) + modes = lab043.process_bpsk_modes(capture, expected, 0.0, 1.0, 16) + + paths = lab043.save_diagnostic_artifacts( + tmp_path, + calibrations, + modes, + _complete_acceptance(), + ) + + assert len(paths) == 5 + assert all(path.exists() and path.stat().st_size > 0 for path in paths) + summary = (tmp_path / "lab043_summary.csv").read_text(encoding="utf-8") + assert "radio_passed" in summary + assert summary.rstrip().endswith(",0") + assert (tmp_path / "lab043_calibration.png").stat().st_size > 1000 + + +def test_all_embedded_lab043_checks_pass_without_hardware() -> None: + results = lab043.run_functional_tests() + assert len(results) == 10 + assert all(result.passed for result in results), results