Разработка цифрового радиолюбительского протокола на базе OFDM. PHY-уровень

Статья посвящена разработке радиолюбительского протокола связи на базе OFDM (Orthogonal Frequency-Division Multiplexing). В ней рассматривается процесс создания физического уровня, сложности, с которыми пришлось столкнуться, а также его испытания в эфире.

Несмотря на то, что статья в первую очередь ориентирована на радиолюбителей, она охватывает сразу несколько смежных тем, таких как беспроводные сети, методы коррекции ошибок, мультиплексирование, цифровая обработка сигналов. Статья может быть интересна тем, кто интересуется устройством и принципами работы беспроводных сетей, какие математические принципы в них используются.

Предыстория

Чуть больше года назад, 28.07.2025г. на хабр была опубликована статья о радиолюбительском протоколе VarAC, который по описанию выглядел как что-то "невероятное", позволяющее превратить радиолюбительский SSB-трансивер в сетевую карту и передавать любые двоичные данные (а также играть P2P в незамысловатые игры, такие как реверси или шашки с шахматами), при этом иметь относительно высокую скорость и быть устойчивым к помехам. По описанию, VarAC основан на применении OFDM — методе модуляции с мультиплексированием на ортогональных частотах, методе, широко используемом в цифровом телевидении, Wi-Fi и LTE/5G сетях, тем самым вызвав интерес, как удалось адаптировать это для любительской радиосвязи.

Так как VarAC является разработкой с закрытым исходным кодом, а также требует платных подписок и вообще, вся из себя полувоенная, было решено самостоятельно разобраться с OFDM и попробовать сделать свою реализацию любительского протокола.

TL;DR

Далее в статье идет лонгрид с большим количеством математических формул и кода на Python. Тем, кому это не интересно, можно сразу перемотать на параграфы "Тестирование в симуляторе" и "Тестирование в эфире".

Под спойлером в формате вопрос-ответ приведено краткое повествование статьи:

Hidden text

• В: как работает?
  О: используется OFDM с PSK модуляцией, сформированный сигнал передается и декодируется звуком.

• В: как передается сигнал?
  О: используется любительский SSB-трансивер.

• В: какие PSK используются?
  О: BPSK и QPSK (QAM в перспективе).

• В: устойчив ли к многолучевости?
  О: да, применяется циклический префикс длиной 25% от символа.

• В: работает ли в эфире?
  О: в лабораторных условиях подтверждена работоспособность в реальном эфире.

• В: может ли работать под шумами?
  О: да, работает под шумами, максимальный предел чувствительности -9 дБ, рабочий предел -7 дБ.

• В: какая канальная скорость передачи данных?
  О: от 353.9823 до 707.9646 бит/сек.

• В: какая скорость передачи данных?
  О: BPSK от 117 до 265 бит/сек; QPSK от 235 до 530 бит/сек.

• В: сколько несущих и пилотов в сигнале?
  О: рассматривается модель 23 поднесущих, из которых 7 преамбул.

• В: какой диапазон частот между поднесущими?
  О: 93.75 Гц

• В: можно ли плотнее?
  О: да, например модель шагом 46.87 Гц, но она не рассматривается в статье.

• В: устойчив ли протокол к сдвигу частоты?
  О: да, устойчив, в пределах 300Гц.

• В: методы коррекции ошибок?
  О: применяется сверточный код с базовой скоростью ⅓.

• В: можно ли быстрее?
  О: да, она трансформируется в скорости ½, ⅔ и ¾.

• В: как декодируется сверточный код?
  О: алгоритм витерби на мягких решениях бит.

• В: как устроен процесс синхронизации?
  О: используется функция с идеальной автокорреляционной функцией (последовательность Задова-Чу).

• В: какая полоса частот используется?
  О: полоса от 300 до 2700 Гц.

• В: на каких диапазонах планируется передача?
  О: тест был на УКВ 2 м (144 МГц), формальных ограничений нет.

• В: чем этот протокол лучше чем FT8 или VarAC?
  О: ничем, FT8 использует принципиально другой тип модуляции и скорость; VarAC — закрытый код.

Немного теории

На тему OFDM уже написано много тем, проливающих свет на принцип его работы, в данном разделе будут поверхностно рассмотрены только основные принципы и особенности.

Для тех, кто знаком с устройством OFDM этот раздел можно пропустить. Более детально особенности рассматриваются в следующем разделе.

Кратко о устройстве OFDM

Orthogonal Frequency Division Multiplexing — это метод цифровой модуляции, разделяющий один высокоскоростной поток данных на множество медленных, передающий их одновременно на разных частотах.

Как видно из названия, особенностью этой модуляции является ортогональность частот, это означает, что когда сигнал одной из несущей имеет максимальное значение амплитуды, значения амплитуд остальных несущих равны нулю.

Ортогональность частот достигается за счет применения обратного преобразования Фурье.

Формирование сигнала осуществляется в частотной области, в которой задаются комплексными значениями задаются амплитуды и фазы частот, затем подаваемые в ОБПФ (IFFT), который, в свою очередь, формирует комплексный IQ-сигнал во временной области.

Рисунок 1. Формирование OFDM-символа через IFFT.

Рисунок 2. Формирование OFDM-символа из гармонических частот.

Так как в частотной области формирование осуществляется путем задавания амплитуды и фазы, символьное кодирование (маппинг) может осуществляться фазовой манипуляцией, такими как BPSK, QPSK и QAM разных порядков, позволяя передавать на одну несущей от одного до бит информации; при этом, чем больше информации переносится несущей, тем сильнее она подвержена влиянию помех и требует большего отношения сигнал/шум (SNR).

На принимающей стороне осуществляется обратное действие — через [прямое] преобразование фурье (FFT) из сигнала извлекаются комплексные значения амплитуд и фаз каждой из несущих.

Таким образом, базовая реализация OFDM заключается в использовании обратного преобразования Фурье для кодирования сигнала и прямого преобразования для его декодирования и извлечения из него информации.

Циклический префикс

Одна из особенностей OFDM-систем заключается в применении циклического префикса, задача которого минимизировать межсимвольную интерференцию (ISI — Inter Symbol Interferenct). Так как в основном OFDM применяется в высокочастотных системах, где радиоволны начинают активно проявлять свойства, близкие оптическим, что проявляется в виде эха и интерференций при многократном отражении и многолучевом распространении сигнала (рисунок 3). Так как любая интерференция разрушает OFDM-символ на принимающей стороне, нужно разнести символы во времени так, чтобы при приеме они не наползали друг на друга. Для этой цели можно было бы добавить защитный интервал, когда передатчик просто не передает, но такой подход приводит к усложнению частотного выравнивания на приемнике, разрушению ортогональности частот и возникновению интерференций между несущими (ICI — Inter Carrier Interference).

Циклический префикс позволяет минимизировать подобные проблемы. Циклическим он назван потому, что является копией хвоста передаваемого OFDM-символа и, как следствие, сохраняет когерентность, фазу и ортогональность несущих.

Рисунок 3. Многолучевое распространение сигнала.

Размер циклического префикса подбирается так, чтобы его продолжительность в эфире составляла максимальное расчетное время, в течении которого присутствует многолучевость. При этом она не должна превышать длительность самого символа.

Рисунок 4. Формирование циклического префикса.

Рисунок 5. Осциллограмма OFDM-символа с циклическим префиксом.

Помимо борьбы с ISI/ICI циклический префикс также позволяет определять и компенсировать частотный дрейф (CFO — Carrier Frequency Offset), а также позволяет применить круговую свертку при обработке сигнала, вместо линейной.

Рисунок 6. Межсимвольная интерференция в зоне циклического префикса.

Таким образом, все возможные искажения и интерференции, вызванные многолучевостью, остаются в циклическом префиксе, который на приемнике исключается из процесса демодуляции. Это позволяет получить менее искаженный OFDM-символ, но, для его корректной работы требуется очень точная синхронизация во времени.

Пилот-сигналы

OFDM-системы как правило использую достаточно широкий диапазон частот для передачи данных, при этом радиосигнал в канале может подвергаться быстрым, медленных и частотно-селективных затуханиям, а также интерферировать с другими сигналами, что приводит к сильным искажениями фазовых созвездий, затрудняя (или делая невозможным) декодирование сигнала. Для того, чтобы связь была возможной, помимо поднесущих, передающих данные, в OFDM выделяется группа поднесущих, называемых пилотами. Их задача заключается в передаваче эталонного сигнала постоянной амплитуды. Так как и приемник и передатчик заранее знают на каких несущих, с какой амплитудой и фазой передаются пилоты, приемник эквализирует (например, используя линейную интерполяцию) поднесущие принимаемых OFDM-символов, тем самым компенсируя искажения фазовых созвездий, что значительно повышает вероятность правильного декодирования сигнала.

Рисунок 7. Процесс эквализации принятых символов по пилот-сигналам.

Помимо эквализации пилот-сигналы используются для оценки CFO и синхронизации. В OFDM-системах пилоты могут передаваться как в составе одного символа вместе с данными, так и согласно сетке. При использовании сетки, пилоты также выступают маркерами синхронизации.

Концепция протокола

Из статей про VarAC было сделано заключение, что если та программа способна работать в диапазоне звуковых частот, поверх SSB, то можно попробовать пойти тем же путем.

Примечание автора:

В начале разработки казалось, что всего-то делов — BPSK, FFT/IFFT и срез массива для получения префикса и квази-физический уровень протокола готов. Так-то оно так, да не так, эта простота на физическом уровне позже проявится в сложностях уровнями выше.

VarAC работает с сигналами в акустической полосе частот от 300 Гц до 2700 Гц, располагая несущие с шагом примерно в 50 Гц.

Так как большинство радиолюбительских SSB трансиверов для передачи голоса (т.н. телефон) использую полосу частот 2.7-3 КГц, для формирования сигнала была выбрана частота дискретизации 12 КГц, что дает запас частот до 6 КГц.

DEFAULT_SAMPLE_RATE = 12000
sample_rate=DEFAULT_SAMPLE_RATE,

Шаг между несущими было решено сделать больше, чем в VarAC, так как любительские трансиверы не могут гарантировать высокую точность частоты генератора опорного сигнала и при передаче возможно расхождение частот в десятки герц. Размер FFT был задан как 128 полос, что при заданной частоте дискретизации составляет 12000/128=93.75 Гц, при этом допустимый CFO будет в пределах +-46.875 Гц. При меньших значениях вероятность, что несущие будут располагаться не на своих местах становится выше, что приводит к неправильной демодуляции сигнала.

fft_bins=128,

Примечание автора:

На начальных этапах разработки эта проблема была критической. Символ синхронизирующей последовательности при смещении частоты давал пик корреляции со смещением во времени, что приводило не только неправильную оценку CFO но и начала передачи. В последствии эту проблему удалось обойти, применив двухступенчатую преамбулу и в симуляции CFO удавалось компенсировать более чем на 300 Гц.

Полоса частот ограничена диапазоном 300-2400 Гц.

freq_lo=300,
freq_hi=2400,

При такой конфигурации получается 23 несущих, которым соответствуют индексы: 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25.

Центральная полоса №14 соответствующая частоте 1312.5 Гц.

self._sample_rate = sample_rate
self._fft_bins = fft_bins
self._freq_step = self._sample_rate / self._fft_bins
self._freq_lo = freq_lo
self._freq_hi = freq_hi
self._bin_lo = int(self._freq_lo / self._freq_step)
self._bin_hi = int(self._freq_hi / self._freq_step)
self._bin_center = (self._bin_lo + self._bin_hi) // 2

self._channel_indices = np.arange(self._bin_lo, self._bin_hi + 1)
self._subcarriers = len(self._channel_indices)

Примечание:

Индексы 125, 124, 123, 122, 121, 120, 119, 118, 117, 116, 115, 114, 113, 112, 111, 110, 109, 108, 107, 106, 105, 104, 103 — область отрицательных частот, комплексно сопряженных с вышеприведенными.

Остальные полосы частот всегда заполняются нулями и не переносят информацию.

При этих параметрах длительность одиночного ODFM-символа без циклического префикса составляет 128/12000=0.01 сек.

Из 23 несущих выделено 7 для передачи пилот-сигналов.

self._pilots_count = pilots_count

self._pilot_local_indices = np.linspace(0, len(self._channel_indices) - 1, self._pilots_count, dtype=int)
self._data_local_indices = np.delete(np.arange(len(self._channel_indices)), self._pilot_local_indices)

self._pilot_carriers = self._channel_indices[self._pilot_local_indices]
self._pilot_carriers_len = len(self._pilot_carriers)

self._data_carriers = self._channel_indices[self._data_local_indices]
self._data_carriers_len = len(self._data_carriers)

Размещение пилотов формируется так, чтобы минимум два пилота находились на границах спектра; остальные поднесущие распределены равномерно. Таким образом, пилоты используют полосы 3, 6, 10, 14, 17, 21, 25 (комплексно-сопряженные с ними 125, 122, 118, 114, 111, 107, 103).

Рисунок 8. Размещение несущих в спектре.

Кодирование информации осуществляется фазовой манипуляцией. BPSK (Binary Phase Shift Keying) является самой помехоустойчивой, но при этом является и самым медленным, так как переносит 1 бит информации. Этот вид манипуляции используется для передачи заголовков пакетов и при плохом качестве канала.

Рисунок 9. Созвездие для BPSK.

На рисунке 9 приведено созвездие для BPSK. Схема кодирования задана в соответствии с таблицей:

_MAPPING_TABLE = {
    (0,): -1 + 0j,
    (1,): 1 + 0j,
}

Помимо BPSK было реализовано более быстрое кодирование — QPSK, созвездие которого приведено на рисунке 10. При этом кодировании, объем передаваемой информации становится вдвое больше, соответственно и увеличивается скорость передачи, однако такой метод кодирования более чувствителен к фазовым искажениям и требует большего значения отношения сигнал/шум (SNR).

Рисунок 10. Созвездие для QPSK.

Таблица кодирования для QPSK:

_NORM = np.sqrt(2.0)
_MAPPING_TABLE = {
    (0, 0): (-1.0 - 1j) / _NORM,
    (0, 1): (-1.0 + 1j) / _NORM,
    (1, 1): (1.0 + 1j) / _NORM,
    (1, 0): (1.0 - 1j) / _NORM,
}

В QPSK биты кодируются парами, по этой причине ключи _MAPPING_TABLE заданы кортежами.

Для повышения помехоустойчивости протокола и возможности осуществления связи при отрицательных значениях SNR было принято решение принести в жертву скорость передачи данных, за счет повторения символов во временной области.

sym_tile=4,

Параметр sym_tile определяет, сколько раз дублируется каждый кодируемый OFDM-символ. На принимающей стороне приемник выполняет когерентное сложение сигнала, тем самым получая прирост SNR на = 6.02 dB. При этом скорость передачи данных снижается в 4 раза.

Примечание автора:

Значение sym_tile поределено эмпирическим путем.

Так как формируемый OFDM-символ является когерентным, операция тайлинга осуществляется путем простой конкатенации массива.

Расчет предела Шеннона

Расчет предела Шеннона позволяет определить максимальную физическую пропускную способность канала с аддитивным белым Гауссовским шумом (AWGN — Additive White Gaussian Noise).

Расчет осуществляется по формуле Шеннона-Хартли:

Где:

• — пропускная способность канала (предел Шеннона) в битах в секунду (бит/с)

• — ширина полосы пропускания канала в Герцах (Гц).

• — линейное отношение мощностей сигнал/шум

Примечание:

Так как значения SNR записываются в децибелах, формула перевода SNR в линейный вид выглядит так:

Подставив в формулу Шеннона-Хартли ширину полосы пропускания 2400-300=2100 Гц и значения SNR можно получить максимальную пропускную способность канала в условиях шума.

В таблице 1 приведены рассчитанные значения предела Шеннона в условиях положительных и отрицательных значений SNR.

SNR (dB)	SNR (линейный)	Предел Шеннона (бит/сек)
5	3.1623	4320.48 (max: 4200.00)
4	2.5119	3805.72
3	1.9953	3323.63
2	1.5849	2877.22
1	1.2589	2468.84
0	1	2100.00
-1	0.7943	1771.23
-2	0.631	1482.01
-3	0.5012	1230.82
-4	0.3981	1015.30
-5	0.3162	832.46
-6	0.2512	678.93
-7	0.1995	551.18
-8	0.1585	445.71
-9	0.1259	359.25

Верхняя граница пропускной способности ограничена частотой Найквиста и не может превышать = 2*2100 = 4200 бит/сек, поэтому значения SNR свыше 5 дБ не имеют смысла.

На основе рассчитанных значений пределов Шеннона можно произвести оценку эффективности передачи данных. Оценка рассчитывается как отношение скорости передаваемой информации к пределу Шеннона.

Исходя из концептуальных параметров, продолжительность одного OFDM-символа с учетом тайлинга и циклического префикса составляет 4*0.0106+0.0026 = 0.0452 сек. Вместимость одного OFDM-символа составляет 16 бит для BPSK и 32 для QPSK, скорости передачи данных при этом составляют 353.9823 бит/сек и 707.9646 соответственно.

В таблице 2 приведены расчеты оценки эффективности протокола:

SNR (dB)	Предел Шеннона (бит/сек)	Модуляция	Скорость протокола (бит/сек)	Эффективность (%)	Оценка
5	4200.00	QPSK	707.96	16.85	Низкая эффективность, можно применить QAM
4	3805.72	QPSK	707.96	18.60	Низкая эффективность, можно применить QAM
3	3323.63	QPSK	707.96	21.30	Эффективен с применением FEC
2	2877.22	QPSK	707.96	24.60	Эффективен с применением FEC
1	2468.84	QPSK	707.96	28.67	Эффективен с применением FEC
0	2100.00	QPSK	707.96	33.71	Эффективен с применением FEC
-1	1771.23	BPSK	353.98	19.98	Низкая эффективность, можно применить QPSK
-2	1482.01	BPSK	353.98	23.88	Эффективен с применением FEC
-3	1230.82	BPSK	353.98	28.75	Эффективен с применением FEC
-4	1015.30	BPSK	353.98	34.86	Эффективен с применением FEC
-5	832.46	BPSK	353.98	42.52	Эффективен с применением FEC
-6	678.93	BPSK	353.98	52.13	Обязательно применение FEC
-7	551.18	BPSK	353.98	64.22	Обязательно применение FEC с понижением скорости
-8	445.71	BPSK	353.98	79.41	Обязательно применение FEC с понижением скорости
-9	359.25	BPSK	353.98	98.53	Предел

Как видно из таблицы 2, теоретический предел работы протокола составляет порядка -9 dB SNR, при этом видно, что возможно применение QAM. Также видно, что QPSK вполне может работать до уровня SNR в -2dB при наличии алгоритмов коррекции ошибок (FEC — Forward Error Correction).

Архитектура протокола

Физический уровень протокола базируется на ортогональном частотном разделении каналов. На рисунке 11 приведены упрощенные структурные схемы передатчика и приемника.

Рисунок 11. Упрощенная схема передающего и приемного трактов.

В передающем тракте данные разбиваются на пакеты и проходят процедуру помехозащищенного кодирования. Закодированные данные подвергаются интерливингу и скремблингу. На основе полученной последовательности формируются PSK-символы, для переноса на поднесущие в IDFT. Обратное преобразование Фурье из полученных данных формирует дискретный сигнал во временной области. К полученному сигналу добавляется преамбула для синхронизации, формируемая также с помощью IDFT.

Схема приемного тракта зеркально противоположна передающему. После успешной синхронизации, сигнал переводится в частотную область преобразованием Фурье (DFT). PSK-демодулятор вычисляет вероятностные (мягкие) значения бит данных. Перед декодированием алгоритмом Витерби, последовательность мягких бит проходит процедуру дескремблирования и деинтерливинга. Декодированные биты проверяются на целостность с помощью контрольных сумм (CRC), и в случае их совпадения данные передаются на прикладной уровень.

Рисунок 11. Структурная схема передающего тракта.

Рисунок 12. Структурная схема приемного тракта.

На рисунках 11 и 12 продемонстрированы основные функциональные блоки и логика их взаимодействия. В правой части рисунков детализированы процессы, входящие в состав ключевых узлов, а также операции, выполняемые на границах их интерфейсов.

Структура пакета

В протоколе реализована пакетная передача данных. Пакет представляет собой модульную структуру (рисунок 13), включающую в себя преамбулу синхронизации и заголовок (Header). Заголовок имеет фиксированную длину, передается строго в BPSK и скоростью кодирования ⅓. Содержащаяся в заголовке служебная информация определяет объем, тип модуляции и скорость кодирования последующего блока данных (Data).

Рисунок 13. Структурная схема пакета данных.

На рисунке 14 приведена структура заголовка пакета.

Рисунок 14. Структура заголовка пакета.

Описание заголовка:

class PacketType(Enum):
    BEACON = 0
    DATA = 4


class ModType(Enum):
    BPSK = 0
    QPSK = 1
    ...


class CCSpeed(Enum):
    R13 = 0
    R12 = 1
    R23 = 2
    R34 = 3


@dataclass
class Header(PacketCRC8):
    PACKET_SIZE: ClassVar[Final[int]] = 17 + PacketCRC8.CRC_SIZE

    ver: int
    typ: PacketType
    mod: ModType
    spd: CCSpeed
    len: int

    def encode(self) -> npt.NDArray[np.uint8]:
        ver = self._to_bits(self.ver, 2)
        typ = self._to_bits(self.typ.value, 3)
        mod = self._to_bits(self.mod.value, 2)
        spd = self._to_bits(self.spd.value, 2)
        len = self._to_bits(self.len, 8)

        data = np.concatenate([ver, typ, mod, spd, len])

        packet = self._append_crc(data)
        return packet

    @classmethod
    def decode(cls, bits: npt.NDArray[np.uint8], check_crc: bool = True) -> 'Header':
        assert bits.shape == (cls.PACKET_SIZE,)

        data = cls._extract_data(bits, check_crc)

        ver = cls._from_bits(data[0:2])
        typ = cls._from_bits(data[2:5])
        mod = cls._from_bits(data[5:7])
        spd = cls._from_bits(data[7:9])
        len = cls._from_bits(data[9:17])

        return cls(ver, PacketType(typ), ModType(mod), CCSpeed(spd), len)

Поле ver (2 бита) определяет версию протокола (в данном случае всегда равна 1). Поле typ (3 бита) определяет тип передаваемых данных. Поля mod и spd (по 2 бита каждый) соответствуют типу модуляции и скорости кодирования данных, что позволяет реализовать адаптивный канал связи. Размер идущего после заголовка блока данных содержится в поле len (8 бит), ограничивая максимальный объем до 255 бит.

Таким образом, объем заголовка составляет 25 бит (17 бит данных и 8 бит CRC).

Пример кодирования:

Header(ver=1, typ= < PacketType.BEACON: 0 >, mod = < ModType.BPSK: 0 >, spd = < CCSpeed.R13: 0 >, len = 90)

В двоичном представлении: 0b01_000_00_00_01011010_10011110.

Контроль целостности передачи служебной информации осуществляется сверкой контрольной суммы CRC-8. Принятый заголовок, в котором не совпадает расчетный и переданный CRC признается ошибочным. Что приводит к прерыванию процесса декодирования и пропуску пакета.

В рамках статьи рассматриваются два типа прикладных сообщений: маяк, посылающий в эфир свой позывной с QTH-локатором (рисунок 15) и пакет произвольных данных (рисунок 16).

Рисунок 15. Структура пакета данных маяка.

Описание пакета данных маяка:

class BeaconMode(Enum):
    ...
    BEACON = 15


@dataclass
class Beacon(PacketCRC8):
    PACKET_SIZE: ClassVar[Final[int]] = 82 + PacketCRC8.CRC_SIZE

    callsign: str
    qth: str
    mode: BeaconMode = BeaconMode.BEACON

    def encode(self) -> npt.NDArray[np.uint8]:
        call_val = callsign_to_int(self.callsign)
        qth_val = qth6_to_int(self.qth)

        call_bits = val_to_bits(call_val, 53)
        qth_bits = val_to_bits(qth_val, 25)
        mode_bits = val_to_bits(self.mode.value, 4)

        data = np.concatenate([call_bits, qth_bits, mode_bits])

        packet = self._append_crc(data)
        return packet

    @classmethod
    def decode(cls, bits: npt.NDArray, check_crc: bool = True) -> 'Beacon':
        assert bits.shape == (cls.PACKET_SIZE,)

        data = cls._extract_data(bits, check_crc)

        call_bits = data[0:53]
        qth_bits = data[53:78]
        flag_bits = data[78:82]

        call_val = bits_to_val(call_bits)
        qth_val = bits_to_val(qth_bits)
        mode = BeaconMode(bits_to_val(flag_bits))

        args = int_to_callsign(call_val), int_to_qth6(qth_val), mode
        return cls(*args)

Данные в полях callsign и qth занимают 53 и 25 бит. Поле mode (4 бита) определяет тип работы маяка, значение BEACON информирует, что передатчик работает только в режиме передающего маяка.

Сигнал маяка рассчитан на прием в условиях сильных шумов, в связи с чем его передача осуществляется с применением BPSK и с самой медленной скоростью кодирования.

Позывной кодируется алгоритмом Base38:

ALPHABET = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ/"
BASE_CALL = len(ALPHABET)
MAX_CALL_LEN = 10


def callsign_to_int(callsign: str) -> int:
    s = callsign.strip().upper()[:MAX_CALL_LEN].ljust(MAX_CALL_LEN, " ")
    val = 0
    for char in s:
        val = val * BASE_CALL + ALPHABET.index(char)
    return val

def int_to_callsign(val: int) -> str:
    chars = []
    for _ in range(MAX_CALL_LEN):
        chars.append(ALPHABET[val % BASE_CALL])
        val //= BASE_CALL
    return "".join(reversed(chars)).strip()

Примечание автора:

В данном случае можно было добиться еще большей плотности кодирования позывного воспользовавшись алгоритмами из FT8, но фактическая экономия в 4 бита слишком несущественна в рамках прототипа.

Пример кодирования позывного R9FEU: 4671407026402144 (0b10000100110001001111010110100011010001111001101100000)

Алгоритм кодирования QTH-локатора заимствован из FT8, но расширены до 6 значного представления:

def qth6_to_int(qth: str) -> int:
    q = qth.upper()
    assert len(q) == 6
  
    off = ord('A')    
    lon_field = ord(q[0]) - off
    lat_field = ord(q[1]) - off
    lon_sq = int(q[2])
    lat_sq = int(q[3])
    lon_sub = ord(q[4]) - off
    lat_sub = ord(q[5]) - off

    val = lon_field
    val = val * 18 + lat_field
    val = val * 10 + lon_sq
    val = val * 10 + lat_sq
    val = val * 24 + lon_sub
    val = val * 24 + lat_sub
    return val

def int_to_qth6(val: int) -> str:
    lat_sub = val % 24
    val //= 24
    lon_sub = val % 24
    val //= 24
    lat_sq = val % 10
    val //= 10
    lon_sq = val % 10
    val //= 10
    lat_field = val % 18
    val //= 18
    lon_field = val % 18

    off = ord('A')
    return f"{chr(lon_field + off)}{chr(lat_field + off)}{lon_sq}{lat_sq}{chr(lon_sub + off)}{chr(lat_sub + off)}"

Пример кодирования 6 значного локатора LO88ca: 12261936 (0b101110110001101000110000).

Ниже приведены вспомогательные функции конвертации данных в массивы бит:

@staticmethod
def _to_bits(val: int, count: int) -> npt.NDArray[np.uint8]:
    val = val & ((1 << count) - 1)
    bytes_count = (count + 7) // 8
    val_bytes = val.to_bytes(bytes_count, byteorder='big')
    val_array = np.frombuffer(val_bytes, dtype=np.uint8)
    return np.unpackbits(val_array)[-count:]

@staticmethod
def _from_bits(bits: npt.NDArray[np.uint8]) -> int:
    pad_width = (8 - (len(bits) % 8)) % 8
    padded_bits = np.pad(bits, (pad_width, 0), mode="constant", constant_values=0)
    bytes_array = np.packbits(padded_bits)
    return int.from_bytes(bytes_array.tobytes(), byteorder='big')

Рисунок 16. Структура пакета произвольных данных.

Описание пакета произвольных данных:

@dataclass
class Data(PacketCRC16):
    reserved: int
    payload: bytes

    def encode(self) -> npt.NDArray[np.uint8]:
        reserved = self._to_bits(self.reserved, 20)
        payload = np.unpackbits(np.frombuffer(self.payload, dtype=np.uint8))

        data = np.concatenate([reserved, payload])

        packet = self._append_crc(data)
        return packet

    @classmethod
    def decode(cls, bits: npt.NDArray[np.uint8], check_crc: bool = True) -> 'Data':
        assert bits.shape >= (44,)

        data = cls._extract_data(bits, check_crc)

        reserved = cls._from_bits(data[0:20])
        payload = np.packbits(data[20:]).tobytes()

        return cls(reserved, payload)

Поле reserved (20 бит) предназначено для передачи служебной информации о передаваемых в пакете данных. Поле payload произвольной длины в диапазоне от 1 до 27 байт.

Так как размер пакета составляет десятки байт, для контроля целостности используется 16 битный CRC.

Аналогично заголовку, несовпадение контрольных сумм в принятом пакете приводит к его игнорированию.

CRC-8/CRC-16

Контроль целостности данных в принимаемых пакетах осуществляется алгоритмами  CRC-8 для пакетов малой длины и CRC-16 для пакетов больших размеров.

def crc(bits: npt.NDArray[np.uint8], polynomial: int, seed: int, crc_bits: int) -> int:
    mask = (1 << crc_bits) - 1

    crc = seed & mask
    for bit in bits:
        crc_msb = (crc >> (crc_bits - 1)) & 1

        crc = (crc << 1) & mask

        if crc_msb ^ (bit & 1):
            crc ^= polynomial

    return crc

Функция crc принимает на вход массив информационных бит bits, значение полинома polynomial. Аргумент seed определяет начальное состояние регистров для расчета CRC; аргумент crc_bits определяет битность результата.

Для расчета контрольных сумм взяты проверенные полиномы из LTE и Wi-Fi:

POLY_ITU_LTE = 0x07
POLY_CCITT_WIFI = 0x1021

Функции crc8_lte и crc16_ccitt реализуют обертки, вычисляющие значения соответствующих CRC для передаваемых данных:

def crc8(bits: npt.NDArray[np.uint8], polynomial: int) -> int:
    return crc(bits, polynomial, 0xff, 8)


def crc16(bits: npt.NDArray[np.uint8], polynomial: int) -> int:
    return crc(bits, polynomial, 0xffff, 16)


def crc8_lte(bits: npt.NDArray[np.uint8]) -> int:
    return crc8(bits, POLY_ITU_LTE)


def crc16_ccitt(bits: npt.NDArray[np.uint8]) -> int:
    return crc16(bits, POLY_CCITT_WIFI)

Для массива бит 0b1000010011000100111101011010001101000111100110110000001011101100011000000001110000 CRC-8 соответствует значению 0b10101101.

Для данных 0b00000011000001010111001000000010000000100000001000000101010001101000011011110111010101100111011010000010000001110100011010000110100101110011001000000110001001100101001000000110110101100001011001000110111001100101011100110111001100101100 CRC-16 равен 0b1001001100011111.

Механизм проверки CRC в пакетах данных:

@classmethod
def _check_crc(cls, data: npt.NDArray[np.uint8], crc_data: npt.NDArray[np.uint8],
                    raise_exception: bool = True) -> bool:
    crc_our = cls._calc_crc(data)
    crc_their = cls._from_bits(crc_data)

    ok = crc_our == crc_their

    if raise_exception and not ok:
        raise PacketCRCMissmatch("Packet CRC missmatch")

    return ok

Преамбула пакета

На стороне приемника происходит непрерывная потоковая обработка данных, принимаемых из эфира. Чтобы приемник мог точно определить присутствие в принимаемом сигнале пакетов, передатчик перед отправкой OFDM-символов формирует преамбулу пакета.

Преамбула пакета позволяет решить сразу несколько задач:

1. определить в сигнале присутствие пакета данных;

2. определить CFO;

3. определить начало OFDM-символов.

Пункт номер 1 по сути является следствием пунктов 2 и 3.

Так как OFDM требует коррекции частоты и очень точной синхронизации по времени, автокорреляционная функция (АКФ) должна иметь одиночный очень четкий пик корреляции, при этом иметь как можно сильнее подавленные боковые лепестки (побочные максимумы). Такими свойствами АКФ обладает последовательность Задова-Чу (Zadoff-Chu, далее ZC-последовательность), широко применяемая в системах связи.

Последовательность Задова-Чу относится к группе CAZAC-последовательностей, то есть Constant Amplitude Zero AutoCorrelation. Ее амплитуда всегда постоянна, а АКФ представляет собой дельта-функцию и дает один единственный максимальный пик корреляции (рисунок 19).

Функция генерации ZC-последовательности:

@staticmethod
def _gen_zc_seq(root: int, length: int) -> npt.NDArray[np.complex64]:
    n = np.arange(length)
    seq = np.exp(-1j * np.pi * root * n * (n + length % 2) / length)
    return seq

ZC-последовательность обладает уникальным свойством: ее циклическая автокорреляционная функция равна дискретной дельта-функции (рисунок 19), при условии, что корневое значение и длина последовательности являются взаимно простыми числами.

В реализации протокола были выбраны значения:

zc_preamble_root=17,
zc_preamble_len=23,

Длина последовательности была выбрана равной 23, так как это число является простым. Значение корня задано как 17, поскольку оно является взаимно простым для числа 23 (при этом еще и само является простым). Значение корня должно быть строго меньше длины последовательности, так как из-за ее цикличности любые другие корни будут тождественны остатку от их деления на длину последовательности.

Примечание автора:

Как можно заметить, количество несущих и длина ZC-последовательности равны 23, такая конфигурация была подобрана специально, чтобы полностью вместить преамбулу в частотный диапазон и не оставлять пустых поднесущих и не ломать логику генерации последовательности.

Рисунок 17. Созвездие ZC-последовательности.

Рисунок 18. Осциллограмма сигнала для ZC последовательности.

На рисунках 17 и 18 представлены фазовое созвездие и комплексная осциллограмма сигнала для последовательности Задова-Чу на приведенных параметрах.

Рисунок 19. Графики периодический и апериодической АКФ для ZC-последовательности.

Сильное подавление боковых лепестков АКФ (рисунок 19) возникает за счет деструктивной интерференции фазовых сдвигов.

Преамбула с ZC-последовательностью формируется как обычный OFDM символ. ZC-последовательность обладает высокой помехоустойчивостью и может успешно обнаруживаться при очень низком SNR. Для дополнительного повышения помехоустойчивости (примечание автора: в рамках экспериментов), преамбула подвергается операции тайлинга, при которой преамбула повторяется несколько раз.

zc_preamble_count=4,

Для защиты от ISI и обеспечения круговой свертки в частотной области к сформированной последовательности добавляется циклический префикс.

def _gen_preamble_symbol(self):
    symbol_freq = np.zeros(self._fft_bins, dtype=np.complex64)

    zc_seq = self._gen_zc_seq(self._zc_preamble_root, self._zc_preamble_len)

    start_bin = self._bin_center - (self._zc_preamble_len // 2)
    symbol_freq[start_bin: start_bin + self._zc_preamble_len] = zc_seq

    symbol = self._idft(symbol_freq)
    return symbol
...
def gen_zc_preamble(self) -> npt.NDArray[np.complex64]:
    preamble_zc = self._gen_preamble_symbol()
    preamble_base = np.tile(preamble_zc, self._zc_preamble_count)
    preamble_block = self._add_cp(preamble_base)

    return preamble_block

ZC-последовательность позволяет достичь очень высокой точности при обнаружении сигнала во временной области, но дает низкую оценку CFO, при этом, обладает свойством, что если CFO будет кратно количеству несущих, пик корреляции сместится во времени, что приведет к неправильному декодированию OFDM, при этом оценка CFO окажется неправильной. Решение этой задачи требует ресурсоемкого анализа сигнала, вычисления грубой оценки CFO согласно сетке гипотез и последующей точной оценке фаз.

Помимо ZC-преамбулы, в начало сигнала добавляется тональная преамбула, представляющая собой последовательность из 4 тонов, сменяющих друг друга по времени.

Примечание автора:

Главная причина добавления еще одной преамбулы заключается в том, что ZC-последовательность достаточно чувствительна к искажениям и при первом же тестировании на реальном эфире преамбула разрушалась из-за того, что на принимающей стороне был включен АРУ (AGC), приводящий к клиппингу и перегрузке принимаемого сигнала. С помощью тональной преамбулы планировалось нивелировать этот эффект и использовать ее для стабилизации АРУ. Но, забегая вперед, от разрушительного воздействия АРУ это помогло, к сожалению, слабо, но позволило реализовать менее ресурсоемкую детекцию сигнала с предварительной компенсацией CFO, устраняющей эффект сдвига ZC-преамбулы по времени.

Как и OFDM-символы, тональная преамбула формируется в частотной области. Относительно центральной частоты располагаются 4 тона с промежутком в 4 полосы друг от друга.

self._newman_preamble_bins = np.arange(self._bin_center - 6, self._bin_center + 6 + 1, 4)

Индексы полос: 8, 12, 16, 20.

Одиночная четверка тонов позволяет определить вероятность присутствия преамбулы в сигнале, но не дает оценку во времени, и, как следствие, приведет к ложно-положительным срабатываниям детектора. Для этого формируется еще одна четверка тонов со сдвигом в 2 полосы относительно предыдущей, что снижает вероятность ложноположительных срабатываний, а также сужает область обнаружения полезного сигнала.

self._newman_preamble_shifts = [0, 2]

Индексы полос: 10, 14, 18, 22.

При синтезе многочастотных сигналов в частотной области через IDFT совпадение фаз приводит к конструктивной интерференции, в результате которой формируется пик высокой интенсивности, эквивалентный дельта-функции; говоря простым языком, образуется громкий щелчок.

В частотной области синтез тонов рассчитывается по формуле Ньюмана, предназначенной для вычислением значения фаз тонов, тем самым избегая увеличения пик-фактора (в отличие от CAZAC-последовательностей, не обеспечивает константной амплитуды).

@staticmethod
def _gen_newman_phases(M: int) -> npt.NDArray[np.float64]:
    m = np.arange(M)
    newman_phases = (np.pi * (m ** 2)) / M
    return newman_phases

Функция gen_newman_preamble_tiles формирует заполняет спектр в частотной области.

def gen_newman_preamble_tiles(self) -> npt.NDArray[np.complex64]:
    phases = self._gen_newman_phases(M=len(self._newman_preamble_bins))
    tone_values = np.exp(1j * phases).astype(np.complex64)

    sym_mat = np.zeros((len(self._newman_preamble_shifts), self._fft_bins), dtype=np.complex64)

    for i, shift in enumerate(self._newman_preamble_shifts):
        current_bins = self._newman_preamble_bins + shift

        sym_mat[i, current_bins] = tone_values
        sym_mat[i, self._fft_bins - current_bins] = np.conj(tone_values)

    return sym_mat

Функция gen_newman_preamble переводит ее во временную область. Как упоминалось выше, тональная преамбула изначально проектировалась для стабилизации АРУ, по этой причине управление ее продолжительностью осуществляется значением _newmanpreamble_tile.

self._newman_preamble_tile = 10
...
def gen_newman_preamble(self) -> typing.Iterator[npt.NDArray[np.complex64]]:
    sym_mat = self.gen_newman_preamble_tiles()
    for i, symbol in enumerate(sym_mat):
        signal = self._idft(symbol)
        signal_tile = np.tile(signal, self._newman_preamble_tile * (2 if i == 0 else 1))

        yield signal_tile

Для первой четверки тонов продолжительность в два раза больше, чем у последующей, так как планировалось это время выделить для участка нарастания сигнала в АРУ приемника.

def gen_preamble(self) -> typing.Iterator[npt.NDArray[np.complex64]]:
    yield self.gen_zc_preamble()
...
def gen_preamble(self) -> typing.Iterator[npt.NDArray[np.complex64]]:
     yield from self.gen_newman_preamble()
     yield from super().gen_preamble()

Таким образом преамбула сигнала формируется из конкатенации тональной и ZC-последовательностей.

На рисунке 20 приведен спектр составной преамбулы.

Рисунок 20. Осциллограмма и спектр тональной и ZC-преамбул.

Поиск преамбулы

Анализ сигнала на наличие полезного сигнала осуществляется в два этапа: поиск тональной преамбулы и поиск ZC-преамбулы.

def detect_preamble(self, signal: npt.NDArray[np.complex64]) -> typing.Optional[typing.Tuple[int, float]]:
   return self.detect_zc_preamble(signal)
...
def detect_preamble(self, signal: npt.NDArray[np.complex64]) -> typing.Optional[typing.Tuple[int, float]]:
    if not (coarse := self.detect_newman_preamble(signal)):
        return None

    coarse_sample, coarse_cfo = coarse
    coarse_shifted = freq_shift(self._sample_rate, signal[coarse_sample:], coarse_cfo)

    if not (fine := super().detect_preamble(coarse_shifted)):
        return None

    fine_sample, fine_cfo = fine

    cfo = coarse_cfo + fine_cfo
    return (fine_sample + coarse_sample, cfo)

Фазы тональной преамбулы не несут информации, задача поиска сводится к анализу спектра и сопоставлению амплитуд эталонного и анализируемого сигналов. Этот метод поиска проще когерентной детекции при наличии неизвестного CFO.

def detect_newman_preamble(self, signal: npt.NDArray[np.complex64]) -> typing.Optional[typing.Tuple[int, float]]:
    N = self._fft_bins
    T = self._newman_preamble_tile
    bin_width = self._sample_rate / N

    sym_mat = self.gen_newman_preamble_tiles()
    K = sym_mat.shape[0]
    masks = (np.abs(sym_mat) > 0).astype(np.float32)

    cfo_shifts = np.arange(-self._max_cfo, self._max_cfo + 1)
    num_shifts = len(cfo_shifts)
...

Для грубой оценки CFO формируется матрица эталонных сигналов с заранее заданным сдвигом частоты.

...
    shifted_masks_mat = np.zeros((K, num_shifts, N), dtype=np.float32)
    for s_idx, shift in enumerate(cfo_shifts):
        if shift >= 0:
            shifted_masks_mat[:, s_idx, shift:] = masks[:, :-shift] if shift > 0 else masks
        else:
            shifted_masks_mat[:, s_idx, :shift] = masks[:, -shift:]

    num_blocks = len(signal) // N
    total_preamble_blocks = 2 * T + (K - 1) * T
    if num_blocks < total_preamble_blocks:
        return None
...

Для каждого блока, равной длине преамбулы, формируется спектр энергии частот.

...
    spectrogram = np.zeros((num_blocks, N), dtype=np.float32)
    for b in range(num_blocks):
        spectrogram[b, :] = np.abs(np.fft.fft(signal[b * N: (b + 1) * N])) ** 2

    best_metric = -1.0
    best_block_idx = 0
    best_cfo_idx = 0
...

Накопление энергий для каждого из проанализированных блоков с последующей нормировкой для рассчета метрик.

...
    for b in range(num_blocks - total_preamble_blocks + 1):
        accum_symbols = np.zeros((K, N), dtype=np.float32)
        current_block_offset = b
        for k in range(K):
            blocks_to_sum = 2 * T if k == 0 else T
            accum_symbols[k, :] = np.sum(
                spectrogram[current_block_offset: current_block_offset + blocks_to_sum, :],
                axis=0
            )
            current_block_offset += blocks_to_sum

        signal_powers = np.sum(accum_symbols[:, np.newaxis, :] * shifted_masks_mat, axis=2)
        total_powers = np.sum(accum_symbols, axis=1)[:, np.newaxis] + 1e-10
        norm_powers = signal_powers / total_powers
...

Накопление и анализ метрик для всей сетки эталонных сигналов с учетом CFO.

...
        metrics_for_all_shifts = np.prod(norm_powers, axis=0)
        max_shift_idx = np.argmax(metrics_for_all_shifts)
        current_best_metric = metrics_for_all_shifts[max_shift_idx]

        if current_best_metric > best_metric:
            best_metric = current_best_metric
            best_block_idx = b
            best_cfo_idx = max_shift_idx
...

Анализ метрики с пороговым значением. Если метрика ниже порога, проанализированный сигнал считается шумом.

...
    threshold = 0.05
    if best_metric < threshold:
        return None
...

Грубая и точная оценки CFO. Грубая оценка осуществляется относительно индекса в спектре. Для точной оценки сигнал сдвигается на значение coarse_cfo_hz, после чего вычисляется автокорреляция преамбулы. Аргумент комплексного числа показывает разность фаз, линейно связанных с остатком CFO.

...
    sample_index = best_block_idx * N
    coarse_cfo_bin = cfo_shifts[best_cfo_idx]
    coarse_cfo_hz = coarse_cfo_bin * bin_width

    preamble_part = signal[sample_index: sample_index + 2 * T * N]

    t = np.arange(len(preamble_part))
    preamble_corrected = preamble_part * np.exp(-1j * 2 * np.pi * coarse_cfo_hz * t / self._sample_rate)

    r_sum = np.sum(np.conj(preamble_corrected[:-N]) * preamble_corrected[N:])
    phase_diff = np.angle(r_sum)

    cfo_error = phase_diff * self._sample_rate / (2 * np.pi * N)
    cfo = coarse_cfo_hz + cfo_error

    return sample_index, cfo

Детектор преамбулы Ньюмана обеспечивает точность оценивания CFO до нескольких герц, но, так как анализ сигнала осуществляется в частотной области, алгоритм теряет фазовую информацию, что приводит к росту вероятности ложноположительных срабатываний в условиях шума.

Ложноположительные срабатывания детектора, финальная коррекция частоты и определение начала блока данных осуществляется детектором ZC-преамбулы, в основе которого применен согласованный фильтр.

def detect_zc_preamble(self, signal: npt.NDArray[np.complex64]) -> typing.Optional[typing.Tuple[int, float]]:
    N = self._fft_bins
    L = self._zc_preamble_count
    CP = self._cyclic_prefix
    preamble_len = L * N

    assert len(signal) >= preamble_len
...

Аналогично предыдущему детектору, формируется матрица гипотез, в которых эталонный сигнал имеет заранее заданный CFO в диапазоне от -max_cfo до max_cfo.

...
    max_cfo = self._max_cfo
    hypotheses = np.arange(-max_cfo, max_cfo + 1)

    global_metric_max = -1.0
    best_m = 0
    best_start_zc_coarse = 0

    preamble_block = self.gen_zc_preamble()
    one_zc_base = preamble_block[-N:]

    sig_energy_single = np.abs(signal) ** 2
    energy_cumsum = np.cumsum(np.insert(sig_energy_single, 0, 0.0))
...

Для каждой гипотезы на основе эталонного сигнала формируется ядро согласованного фильтра (комплексно сопряженный и развернутый сигнал). fftconvolve реализует согласованный фильтр.

...
    for m in hypotheses:
        t = np.arange(N)
        cfo_shift_vector = np.exp(2j * np.pi * m * t / N)
        modulated_zc = one_zc_base * cfo_shift_vector

        one_zc_kernel = np.conj(modulated_zc[::-1])
        corr_single = np.abs(fftconvolve(signal, one_zc_kernel, mode="valid"))

        coarse_corr = np.zeros(len(corr_single) - (L - 1) * N)
        for i in range(L):
            shift = i * N
            coarse_corr += corr_single[shift: shift + len(coarse_corr)]

        local_max_idx = np.argmax(coarse_corr)
        local_max_val = coarse_corr[local_max_idx]
...

Нормализация адаптивным пороговым значением позволяет снизить вероятность ложных срабатываний при получении пика корреляции.

...
        idx_start = local_max_idx
        idx_end = local_max_idx + preamble_len

        if idx_end <= len(signal):
            window_energy = energy_cumsum[idx_end] - energy_cumsum[idx_start]

            epsilon = 1e-8
            if window_energy > epsilon:
                local_metric = local_max_val / (window_energy * preamble_len) ** 0.5
            else:
                local_metric = 0.0

            mean_corr_floor = np.mean(coarse_corr)
            peak_to_floor_ratio = local_max_val / (mean_corr_floor + epsilon)
            base_threshold = self._detection_threshold
            if peak_to_floor_ratio < 3.5:
                dynamic_threshold = base_threshold + 0.1 * (3.5 - peak_to_floor_ratio)
                dynamic_threshold = min(dynamic_threshold, 0.45)
            else:
                dynamic_threshold = base_threshold
        else:
            local_metric = 0.0
            dynamic_threshold = self._detection_threshold
...

Запоминается гипотеза с наилучшими характеристиками.

...
        if local_metric > global_metric_max and local_metric >= dynamic_threshold:
            global_metric_max = local_metric
            best_start_zc_coarse = local_max_idx
            best_m = m

    if global_metric_max <= 0.0:
        return None

    start_zc_coarse = best_start_zc_coarse
    true_start_zc = start_zc_coarse
    est_time = true_start_zc - CP

    if true_start_zc + preamble_len > len(signal) or true_start_zc < 0:
        return None
...

На основе наилучшей гипотезы происходит деротация фазы сигнала, тем самым компенсируя грубое значение CFO. Полученная разность фаз между скорректированным и эталонным сигналом дает точную оценку CFO.

...
    preamble_signal = signal[true_start_zc: true_start_zc + preamble_len]

    t_full = np.arange(preamble_len)
    de_rotated_preamble = preamble_signal * np.exp(-2j * np.pi * best_m * t_full / N)

    sig_current = de_rotated_preamble[: (L - 1) * N]
    sig_delayed = de_rotated_preamble[N: preamble_len]

    r_vector = np.sum(np.conj(sig_current) * sig_delayed)
    phase_diff = np.angle(r_vector)

    fractional_cfo = phase_diff / (2 * np.pi * (N / self._sample_rate))

    subcarrier_spacing = self._sample_rate / N
    final_cfo = (best_m * subcarrier_spacing) + fractional_cfo
    final_time = est_time + len(preamble_block)

    return final_time, final_cfo

Так как ZC-последовательность дает точность определения сигнала с точностью до семпла, в est_time записывается номер сэмпла, где была обнаружена преамбула, с учетом присутствующего в сигнале циклического префикса. Значение final_time указывает на сэмпл, с которого начинаются OFDM-символы.

Пилот-сигналы

Основное назначение пилот-сигналов — передача эталонного сигнала с заранее известными параметрами для последующей оценки канала и коррекции фазовых искажений при демодуляции.

Как было указано выше, в каждом OFDM-символе под пилот-сигналы выделены поднесущие. Использование в качестве эталона простейшего значения 1+0j означает, что на данной поднесущей передается чистый гармонический сигнал с постоянной амплитудой и нулевой фазой. Такой подход приводит к значительному росту пик-фактора сигнала. Для минимизации этого эффекта пилот-сигналы формируются на основе последовательности Задова-Чу.

def _get_pilot_values(self) -> npt.NDArray[np.complex64]:
    return self._gen_zc_seq(self._zc_pilots_root, self._pilots_count)

Для ZC-последовательности заданы следующие параметры:

zc_pilots_root=3,
zc_pilots_count=7,

Примечание автора:

Как и в случае с преамбулой, количество выделенных поднесущих для пилот-сигналов также задано простым числом. При этом корреляционные свойства ZC-последовательности здесь не используются — она применяется только для снижения пик-фактора сигнала.

PSK маппинг и формирование OFDM-символа

Использование нескольких поднесущих в OFDM позволяет реализовать параллельную передачу данных. При этом битовая емкость одного OFDM-символа складывается из битовой емкости поднесущих. Для этого входящий поток бит группируется в блоки, объем которых равен суммарной вместительности одного OFDM-символа.

def _serial_to_parallel(self, bits: npt.NDArray[np.complex64]):
    return bits.reshape((self.data_carriers_len, self._mapper.MU))

Функция _serialto_parallel преобразует исходный последовательный поток бит в матрицу, количество строк в которой соответствует числу используемых несущих, а количество столбцов определяется разрядностью модуляции — то есть количеством бит, способных поместиться в один PSK-символ (параметр MU). Так, для BPSK это значение составляет 1 бит, а для QPSK — 2 бита.

Маппинг битовых групп в комплексные PSK-символы осуществляется методом map класса PSKMapper:

class PSKMapper(ABC):
    ...
    @classmethod
    def map(cls, bits: npt.NDArray) -> npt.NDArray:
        m_table = cls._mapping_table()
        return np.array([m_table[tuple(b)] for b in bits])

class BPSKMapper(PSKMapper):
    MU = 1

    _MAPPING_TABLE = {
        (0,): -1 + 0j,
        (1,): 1 + 0j,
    }

class QPSKMapper(PSKMapper):
    MU = 2

    _NORM = np.sqrt(2.0)
    _MAPPING_TABLE = {
        (0, 0): (-1.0 - 1j) / _NORM,
        (0, 1): (-1.0 + 1j) / _NORM,
        (1, 1): (1.0 + 1j) / _NORM,
        (1, 0): (1.0 - 1j) / _NORM,
    }

На вход метода map подается матрица, сформированная функцией _serialto_parallel. Значение MU () определяет количество бит на символ.

def _map_subsymbols_bits(self, bits: npt.NDArray) -> npt.NDArray[np.complex64]:
    bits_sp = self._serial_to_parallel(bits)
    subsymbols = self._mapper.map(bits_sp)
    return subsymbols

Метод _mapsubsymbols_bits преобразует полученную матрицу бит в массив PSK символов.

Массив PSK-символов преобразуется в OFDM-символ в частотной области:

def gen_symbol(self, bits: npt.NDArray) -> npt.NDArray[np.complex64]:
    subsymbols = self._map_subsymbols_bits(bits)

    symbol = np.zeros(self._fft_bins, dtype=np.complex64)
    symbol[self._channel_indices] = subsymbols

    mirror_bins = self._fft_bins - np.array(self._channel_indices)
    symbol[mirror_bins] = np.conj(subsymbols)

    return symbol

Завершающий этап — перенос OFDM-символа из частотной во временную область с помощью обратного преобразования Фурье:

def _add_cp(self, signal: npt.NDArray[np.complex64]) -> npt.NDArray[np.complex64]:
    cp = signal[-self.cyclic_prefix:]
    return np.hstack([cp, signal])

def _modulate_signal(self, symbol: npt.NDArray[np.complex64]) -> npt.NDArray[np.complex64]:
    signal = self._idft(symbol)
    return signal

def modulate_symbol_cp(self, bits: npt.NDArray) -> npt.NDArray[np.complex64]:
    symbol = self.gen_symbol(bits)
    signal = self._modulate_signal(symbol)
    signal_cp = self._add_cp(signal)
    return signal_cp

Демодуляция и эквализация сигнала

Процесс демодуляции представляет собой обратное преобразование принятого сигнала. На первом шаге удаляется циклический префикс OFDM-символа:

def demodulate_symbol_cp(self, symbol: npt.NDArray[np.complex64]) -> npt.NDArray[np.complex64]:
    assert len(symbol) == self.symbol_len

    signal_no_cp = self._remove_cp(symbol)
    demod_all = self._demodulate_signal(signal_no_cp)
    return demod_all

Так как при формировании сигнала символы дублировались, прием осуществляется когерентным накоплением принятых символов.

Принимаемый сигнал разделяется на одиночные символы и собирается в матрицу.

def _demodulate_signal(self, signal: npt.NDArray[np.complex64]) -> npt.NDArray[np.complex64]:
    repeats = signal.reshape(self.sym_tile, -1)
    tiles = self.sym_tile
...

Поиск относительного сдвига фаз между символами. Применяется комплексная автокорреляция во временной области.

...
    x_data = [0]
    y_data = [0.0]

    step = min(3, tiles - 1)
    phase_diffs = {}

    for step in range(1, step + 1):
        for i in range(tiles - step):
            j = i + step
            R_ij = np.sum(repeats[i] * np.conj(repeats[j]))
            phase_diffs[(i, j)] = np.angle(R_ij)
...

Расчет фаз для каждого отдельного символа относительно самого первого.

...
    for target in range(1, tiles):
        if target <= step:
            x_data.append(target)
            y_data.append(phase_diffs[(0, target)])

        chain_phase = 0.0
        for i in range(target):
            chain_phase += phase_diffs[(i, i + 1)]

        x_data.append(target)
        y_data.append(chain_phase)

        if target >= 2:
            step2_phase = 0.0
            curr = target

            while curr >= 2:
                step2_phase += phase_diffs[(curr - 2, curr)]
                curr -= 2

            if curr == 1:
                step2_phase += phase_diffs[(0, 1)]

            x_data.append(target)
            y_data.append(step2_phase)
...

Полиномиальное сглаживание шума. Сформированные на предыдущем шаге точки аппроксимируются полиномом 2-й степени, что позволяет описать и сгладить дрейф фазы, вызванный CFO. Фаза первого символа принудительно обнуляется.

...
    poly_coeffs = np.polyfit(x_data, y_data, deg=2)

    idx = np.arange(tiles)
    smooth_phases = np.polyval(poly_coeffs, idx)

    smooth_phases -= smooth_phases[0]
...

Восстановление фаз для каждого из символов и прямое преобразование Фурье.

...
    phase_corrs = np.exp(1j * smooth_phases)[:, np.newaxis]
    corrected = repeats * phase_corrs

    acc = np.mean(corrected, axis=0)
    demod = super()._demodulate_signal(acc)
    return demod

Далее OFDM-символ переносится из временной области в частотную путем применения к нему прямого преобразования Фурье:

def _demodulate_signal(self, signal: npt.NDArray[np.complex64]) -> npt.NDArray[np.complex64]:
    demod = self._dft(signal)
    return demod

def demodulate_symbol_cp_soft(self, symbol: npt.NDArray[np.complex64]) -> typing.Tuple[
    npt.NDArray[np.float64], float, npt.NDArray[np.float64]]:

    demod_all = self.demodulate_symbol_cp(symbol)
    demod_channel = demod_all[self._channel_indices]
    ...

Этап частотной коррекции. Производится грубая оценка канала на пилот-сигналах с последующей интерполяцией на спектр символа. На этом этапе выполняется первичная грубая эквализация:

    ...
    pilots_ref = self._get_pilot_values()
    pilots = demod_all[self._pilot_carriers]

    H_est_init = self._channel_estimate_zf(pilots_ref, pilots, self._pilot_carriers)
    eq_init = demod_channel / (H_est_init + 1e-6)
    data_init = eq_init[self._data_local_indices]
    ...

Грубая коррекция фаз принятого сигнала, путем подтягивания к идеальным точкам на созвездиях BPSK или QPSK:

    ...
    mu = self._mapper.MU

    if mu == 1:
        ideal_init = np.sign(np.real(data_init))
        ideal_init[ideal_init == 0] = 1.0
    elif mu == 2:
        ideal_init_real = np.sign(np.real(data_init))
        ideal_init_real[ideal_init_real == 0] = 1.0
        ideal_init_imag = np.sign(np.imag(data_init))
        ideal_init_imag[ideal_init_imag == 0] = 1.0
        ideal_init = (ideal_init_real + 1j * ideal_init_imag) / np.sqrt(2.0)
    ...

Оценка шума путем вычисления среднеквадратичной ошибки (MSE) между принятыми и идеальными значениями созвездия:

    ...
    noise_var_eq_init = np.mean(np.abs(data_init - ideal_init) ** 2)
    mean_H_power_init = np.mean(np.abs(H_est_init[self._data_local_indices]) ** 2)
    curr_noise_var = max(noise_var_eq_init * mean_H_power_init, 1e-6)
    curr_noise_var = min(curr_noise_var, 60.0)
    ...

На основе полученной начальной оценки уровня шума осуществляется точная оценка канала с применением фильтра Винера. После этого к сигналу применяется MMSE-эквалайзер (Minimum Mean Square Error):

    ...
    H_est = self._channel_estimate_wiener(pilots_ref, pilots, self._pilot_carriers, curr_noise_var)


    equalized_all = (demod_channel * np.conj(H_est)) / (np.abs(H_est) ** 2 + curr_noise_var)
    equalized_data_symbols = equalized_all[self._data_local_indices]
    data_symbols = equalized_data_symbols
    ...

Для отфильтрованного сигнала производится повторная, более точная коррекция фаз относительно BPSK и QPSK созвездий:

    ...
    if mu == 1:
        bpsk_dec = np.sign(real_parts)
        bpsk_dec[bpsk_dec == 0] = 1.0

        amplitude = np.mean(real_parts * bpsk_dec)

        noise_var_eq_real = np.mean((real_parts - amplitude * bpsk_dec) ** 2)
        noise_var_eq_imag = np.mean(imag_parts ** 2)
        noise_var_eq = noise_var_eq_real + noise_var_eq_imag
    elif mu == 2:
        qpsk_dec_real = np.sign(real_parts)
        qpsk_dec_real[qpsk_dec_real == 0] = 1.0
        qpsk_dec_imag = np.sign(imag_parts)
        qpsk_dec_imag[qpsk_dec_imag == 0] = 1.0

        amplitude_real = np.mean(real_parts * qpsk_dec_real)
        amplitude_imag = np.mean(imag_parts * qpsk_dec_imag)

        noise_var_eq_real = np.mean((real_parts - amplitude_real * qpsk_dec_real) ** 2)
        noise_var_eq_imag = np.mean((imag_parts - amplitude_imag * qpsk_dec_imag) ** 2)
        noise_var_eq = noise_var_eq_real + noise_var_eq_imag
    ...

Сглаживание шумов в канале:

    ...
    H_data = H_est[self._data_local_indices]
    H_power_sq = np.abs(H_data) ** 2
    mean_H_power = np.mean(H_power_sq)

    updated_noise_var = noise_var_eq * mean_H_power

    alpha = 0.1
    noise_var = alpha * updated_noise_var + (1 - alpha) * curr_noise_var
    noise_var = max(noise_var, 1e-6)
    ...

Рисунок 21. Сравнение BPSK созвездий до и после эквализации (при -9 и -3 dB SNR).

Рисунок 22. Сравнение QPSK созвездий до и после эквализации (при -9 и -3 dB SNR).

Финальный этап демодуляции — расчет значений энергии сигнала на бит к мощности шума (Es/N0) для канала и формирование мягких решений (LLR):

    ...
    es_n0_linear_per_carrier = H_power_sq / noise_var
    es_n0_linear = mean_H_power / noise_var

    if mu == 1:
        llr_outputs = 4.0 * real_parts * es_n0_linear_per_carrier
    elif mu == 2:
        llr_i = 2.0 * real_parts * es_n0_linear_per_carrier
        llr_q = 2.0 * imag_parts * es_n0_linear_per_carrier

        llr_outputs = np.empty(2 * len(data_symbols), dtype=np.float64)
        llr_outputs[0::2] = llr_i
        llr_outputs[1::2] = llr_q

    llr_outputs = np.clip(llr_outputs, -20.0, 20.0)

    return es_n0_linear_per_carrier, es_n0_linear, llr_outputs

Логарифмические значения вероятностей бит (LLR) принудительно ограничены диапазоном от -20 до 20. Ограничение связано с тем, что при высоких значениях SNR диапазон может стать слишком широким, что может привести к невозможности декодирования.

Полученный в методе demodulate_symbol_cp_soft массив мягких бит далее подается на вход декодера Витерби для восстановления информации.

Коррекция ошибок

Передаваемый сигнал подвергается разного рода помехам и замираниями, что приводит к фазовым искажениям OFDM и, как следствие, невозможности правильно декодировать принятую информацию. Коды коррекции ошибок с добавлением избыточности решают задачу восстановления правильности информации на принимающей стороне, приближая канал связи к пределу Шеннона.

Передаваемый объем информации в разрабатываемом протоколе составляет 25 бит для заголовка пакета и до 255 бит для блока данных. При таких объемах данных, одним из эффективных алгоритмов коррекции ошибок является сверточный код.

Примечание автора:

Сверточный код — это не единственное эффективное решение на малых объемах данных, но при этом достаточно простое и быстро декодируется алгоритмом Витерби при использовании коротких полиномов.

Решение остановиться в выборе на сверточном коде также было принято для избежания состояния бесконечной разработки; то есть если на сверточном коде разработка покажет свою работоспособность, то в перспективе можно будет его заменить на более эффективные алгоритмы.

В основе реализации используется не систематический сверточный код со скоростью ⅓ и LTE полиномы: 0o133, 0o171, 0o165 (из которых 0o133, 0o171 — полиномы NASA, зарекомендовавшие себя в космической программе).

class ConvCodec:
    @classmethod
    def _encode(cls, bits: npt.NDArray[np.uint8]) -> npt.NDArray[np.uint8]:
        k = 0
        state = 0
        symbols = np.zeros(len(bits) * len(cls.POLY), dtype=np.uint8)

        for bit in bits:
            state = ((state << 1) | int(bit)) & cls.MASK

            for poly in cls.POLY:
                n = state & poly
                symbols[k] = n.bit_count() % 2
                k += 1

        return symbols

    @classmethod
    def calc_cc_elements(cls, bits_count: int) -> int:
        return (bits_count + cls.PAD_LEN) * cls.SPEED

    @classmethod
    def encode(cls, bits: npt.NDArray[np.uint8]) -> npt.NDArray[np.uint8]:
        data_tailed = np.pad(bits, pad_width=(0, cls.PAD_LEN), mode="constant")
        cc_data = cls._encode(data_tailed)
        return cc_data

Класс ConvCodec реализует базовый алгоритм сверточного кодирования. Метод _encode — алгоритм сверточного кодирования. Метод encode реализует интерфейс кодировщика. Функция calccc_elements предназначена для расчета размера выходных данных при заданной скорости кодирования.

Примечание автора:

В алгоритме не предусмотрено выталкивание последних бит из регистров. Вместо этого перед кодированием к входящему потоку добавляется массивом нулей, длина которого определяется значением PAD_LEN (см. функцию encode).

Согласно границам Шеннона, при высоких значениях SNR базовая скорость передачи ⅓ слишком избыточная и может быть увеличена для повышения эффективности.

Скорости ½, ⅔, и ¾ достигаются путем пунктирования (выкалывания) бит из последовательности, изначально полученной на базовой скорости ⅓.

class ConvCodecPunctured(ConvCodec):
    @classmethod
    def calc_cc_elements(cls, bits_count: int) -> int:
        k = cls.PUNCTURE_MASK.shape[1]
        n = int(np.sum(cls.PUNCTURE_MASK == 1))

        num_bits = bits_count + cls.PAD_LEN
        blocks = (num_bits + k - 1) // k
        return blocks * n

    @classmethod
    def puncture(cls, bits: npt.NDArray[np.uint8]) -> npt.NDArray[np.uint8]:
        mask = cls.PUNCTURE_MASK

        period = mask.shape[1]

        stream = np.array(bits).reshape(-1, cls.SPEED).T
        num_elements = stream.shape[1]

        reps = (num_elements + period - 1) // period
        full_mask = np.tile(mask, (1, reps))[:, :num_elements]

        return stream.T[full_mask.T.astype(bool)]

    @classmethod
    def encode(cls, bits: npt.NDArray[np.uint8]) -> npt.NDArray[np.uint8]:
        cc_data = super().encode(bits)
        cc_data_punct = cls.puncture(cc_data)
        return cc_data_punct

Метод puncture прореживает закодированную последовательность бит согласно матрице PUNCTURE_MASK для увеличения скорости. Значение 0 в матрице исключает бит из выходной последовательности.

class CCLTEBPSK(ConvCodecPunctured, metaclass=MetaConvCodec):
    PUNCTURE_MASK = np.array([[1], [1], [1]])


    POLY = [0o133, 0o171, 0o165]
    K = 7
    SOFT_BIT_SCALE = 1.0

CCLTEBPSK_13 = CCLTEBPSK


class CCLTEBPSK_12(CCLTEBPSK, metaclass=MetaConvCodec):
    PUNCTURE_MASK = np.array([[1, 1], [1, 0], [0, 1]])


class CCLTEBPSK_23(CCLTEBPSK, metaclass=MetaConvCodec):
    PUNCTURE_MASK = np.array([[1, 1], [1, 0], [0, 0]])


class CCLTEBPSK_34(CCLTEBPSK, metaclass=MetaConvCodec):
    PUNCTURE_MASK = np.array([[1, 1, 1], [1, 0, 0], [0, 0, 0]])

Классы CCLTEBPSK_13, CCLTEBPSK_12, CCLTEBPSK_23, CCLTEBPSK_34 реализуют сверточные кодеры с разными скоростями.

Отклик кодировщика со скоростью ⅓ на дельта-функцию (единичный бит): 111100001110111011111.

Пример кодирования на примере данных заголовка:

• заголовок: 0100000000101101010011110;

• скорость ⅓: 000111100001110111011111000000111100110101010011010010110001010011001011100011101011100111000;

• скорость ½: 00111001111101110000111011110101010011010101000110011001101100;

• скорость ⅔: 00110011101100011111101001011001000010010010100;

• скорость ¾: 001100110110001111100100110001001001011100.

Декодирование осуществляется алгоритмом Витерби. На вход подается последовательность мягких бит (LLR — Log-Likelihood Ratio), отражающих апостериорную вероятность принятого бита. Знак определяет наиболее вероятное логическое значение: отрицательные величины соответствуют единице, положительные — нулю. Модуль LLR задает степень уверенности в принятом решении, а значение близкое к нулю соответствует максимальной неопределенности.

Скорость декодера составляет ⅓. Поскольку для повышения скорости применяется пунктирование массива бит, перед декодированием исходная избыточность восстанавливается депунктором. Депунктор вставляет биты с нулевой LLR на позиции удаленных при кодировании бит.

Реализация депунктора:

class ConvCodecPunctured(ConvCodec):
    ...
    @classmethod
    def depuncture(cls, soft_bits: npt.NDArray[np.float128], bits_count: int) -> npt.NDArray[np.float128]:
        mask = cls.PUNCTURE_MASK
        base_streams = mask.shape[0]
        period = mask.shape[1]

        num_elements = bits_count + cls.PAD_LEN
        if (mod := num_elements % period) != 0:
            num_elements += period - mod

        reps = num_elements // period
        full_mask = np.tile(mask, (1, reps))

        depunct_T = np.zeros((num_elements, base_streams), dtype=soft_bits.dtype)

        bool_mask_T = full_mask.T.astype(bool)

        slots = np.sum(bool_mask_T)
        actual_bits = min(soft_bits.size, slots)

        if slots > actual_bits:
            indices = np.where(bool_mask_T)
            bool_mask_T[indices[0][actual_bits:], indices[1][actual_bits:]] = False

        depunct_T[bool_mask_T] = soft_bits[:actual_bits]
        return depunct_T.ravel()

     @classmethod
     def decode(cls, soft_bits: npt.NDArray[np.float128], bits_count: int) -> npt.NDArray[np.uint8]:
        soft_bits_depuct = cls.depuncture(soft_bits, bits_count)
        bits = super().decode(soft_bits_depuct, bits_count)
        return bits

Реализация декодера Витерби:

class ConvCodec:
    ...
    @classmethod
    def _decode(cls, soft_bits: npt.NDArray[np.float128], total_steps: int) -> npt.NDArray[np.uint8]:
        path_metrics = np.full(cls.NUM_STATES, -np.inf)
        path_metrics[0] = 0.0

        tb_bits = np.zeros((total_steps, cls.NUM_STATES), dtype=np.uint8)
        tb_prev_states = np.zeros((total_steps, cls.NUM_STATES), dtype=np.uint8)

        num_polys = cls.EXPECTED_SYMS.shape[2]

        for step in range(total_steps):
            rx_chunk = soft_bits[step * num_polys: (step + 1) * num_polys]
            next_path_metrics = np.full(cls.NUM_STATES, -np.inf)

            for curr_state in range(cls.NUM_STATES):
                if path_metrics[curr_state] == -np.inf:
                    continue

                for bit in (0, 1):
                    next_state = cls.NEXT_STATES[curr_state, bit]
                    exp_syms = cls.EXPECTED_SYMS[curr_state, bit]

                    branch_metric = np.sum(rx_chunk * exp_syms)
                    candidate_metric = path_metrics[curr_state] + branch_metric

                    if candidate_metric > next_path_metrics[next_state]:
                        next_path_metrics[next_state] = candidate_metric
                        tb_bits[step, next_state] = bit
                        tb_prev_states[step, next_state] = curr_state

            path_metrics = next_path_metrics

        decoded_bits = np.zeros(total_steps, dtype=int)
        best_state = np.argmax(path_metrics)

        for step in range(total_steps - 1, -1, -1):
            decoded_bits[step] = tb_bits[step, best_state]
            best_state = tb_prev_states[step, best_state]

        return decoded_bits

    @classmethod
    def decode(cls, soft_bits: npt.NDArray[np.float128], bits_count: int) -> npt.NDArray[np.uint8]:
        soft_bits_cropped = soft_bits[:(bits_count + cls.PAD_LEN) * cls.SPEED]
        bits = cls._decode(soft_bits_cropped, bits_count)
        bits_cropped = bits[:bits_count]
        return bits_cropped

Функция _decode реализует классический алгоритм декодирования методом Витерби. Метод decode реализует интерфейс декодера, параметр bitscount определяет количество выходных бит.

Для оптимизации алгоритма используются предварительно рассчитанные карты переходов алгоритма NEXT_STATES. EXPECTED_SYMS содержит эталонные значения софт-бит формируемые для каждого типа модуляции.

Интерливер

Особенностью сверточных кодов является высокая корректирующая способность при воздействии одиночных ошибок, при этом они критически уязвимы к групповым искажениям. Для снижения вероятности пакетных ошибок в протоколе применен механизм интерливинга. Задача которого заключается в перераспределении информационных битов таким образом, чтобы длинные серии ошибок перераспределялись в независимые одиночные.

def interleave(bits: npt.NDArray, num_carriers: int) -> npt.NDArray:
    data_len = len(bits)

    num_rows = int(np.ceil(data_len / num_carriers))
    total_cells = num_rows * num_carriers

    indices = np.arange(total_cells)
    interleaved_indices = indices.reshape(num_rows, num_carriers).T.flatten()
    pruned_indices = interleaved_indices[interleaved_indices < data_len]

    return bits[pruned_indices]

Функция interleave построчно заполняет матрицу исходными битами, затем читает матрицу по столбцам, то есть выдает транспонированную матрицу.

def deinterleave(symbols: npt.NDArray, num_carriers: int) -> npt.NDArray:
    data_len = len(symbols)

    num_rows = int(np.ceil(data_len / num_carriers))
    total_cells = num_rows * num_carriers

    indices = np.arange(total_cells)
    interleaved_indices = indices.reshape(num_rows, num_carriers).T.flatten()
    pruned_indices = interleaved_indices[interleaved_indices < data_len]

    deinterleaved = np.zeros(data_len, dtype=symbols.dtype)
    deinterleaved[pruned_indices] = symbols

    return deinterleaved

Функция deinterleave решает обратную задачу.

Пример соответствия индексов исходной последовательности из 32 бит после интерливинга на 16 поднесущих:

0, 16, 1, 17, 2, 18, 3, 19, 4, 20, 5, 21, 6, 22, 7, 23, 8, 24, 9, 25, 10, 26, 11, 27, 12, 28, 13, 29, 14, 30, 15, 31

PAPR и Скремблирование

Одной из особенностей OFDM-сигналов является высокий уровень пик-фактора (PAPR), возникающий из-за конструктивной интерференцией фаз отдельных поднесущих. В процессе формирования OFDM-символа их амплитуды могут складываться синфазно, что приводит к резким выбросам, выражающихся в виде щелчков, которые вызывают перегрузку усилителя и клиппингу сигнала. При клиппинге в сигнале возникают внеполосные излучения и высокочастотные помехи, препятствующие декодированию OFDM.

Причина конструктивной интерференции заключается в случайном совпадении фаз поднесущих при определенных комбинациях передаваемых информационных бит (PSK-символов).

Для снижения вероятности возникновения нежелательных комбинаций данных в OFDM-системах применяется скремблирование, придающее входной битовой последовательности свойства псевдослучайного шума, обеспечивая равновероятное распределение нулей и единиц в информационном потоке.

DEFAULT_SEED = 0x5A


def scrambler(data: npt.NDArray, seed: int, operator: typing.Callable[[int, typing.Any], typing.Any]) -> npt.NDArray:
    scrambled = np.zeros_like(data)

    lfsr = seed & 0x7FFF
    if lfsr == 0:
        lfsr = 1

    for i in range(len(data)):
        fb = ((lfsr >> 6) ^ (lfsr >> 3)) & 1
        prbs_bit = lfsr & 1
        scrambled[i] = operator(prbs_bit, data[i])
        lfsr = ((lfsr << 1) | fb) & 0x7FFF

    return scrambled


def scramble(bits: npt.NDArray[np.uint8], seed: int = DEFAULT_SEED) -> npt.NDArray[np.uint8]:
    return scrambler(bits, seed, operator=lambda prbs_bit, value: value ^ prbs_bit)

Функция scrambler реализует скремблирование на основе сдвигового регистра с обратной связью. Так как начальное состояние не может быть нулевым (иначе всегда будет выдаваться 0), выбрано значение 0x5A, поскольку оно обеспечивает высокую частоту переходов между логическими 0 и 1.

Аргумент operator определяет операцию xor над битами.

Ниже приведен пример распределения бит до и после прохождения скремблера:

• исходная последовательность: 011011011110000001110001111111111000110110011111000011001001000000101100100010100110011010111000

• скремблированная последовательность: 011001110101111011100101100100011011001000011000011101011111010010101101100110010111000111100011

Можно заметить существенное снижение длин последовательностей из одинаковых значений бит.

Использование операции XOR позволяет алгоритму скремблирования быть симметричным: повторное применение скремблера восстанавливает исходные данные.

Особенностью дескремблирования является то, что на вход подается последовательность из мягких бит. При работе с мягкими битами логическая операция XOR заменяется на арифметическую операцию смены знака:

def descramble(llr_bits: npt.NDArray[np.float64], seed: int = DEFAULT_SEED) -> npt.NDArray[np.float64]:
    return scrambler(llr_bits, seed, operator=lambda prbs_bit, value: value * (-1.0 if prbs_bit == 1 else 1.0))

Операция Clip and Filter

Независимо от используемых методов снижения PAPR в OFDM, в формируемом сигнале так или иначе остаются пиковые выбросы. Для их дополнительного подавления применяется схема клиппирования и фильтрации (Clip and Filter, CAF), которая ограничивает пиковую амплитуду сигнала и фильтрует возникающие внеполосные спектральные составляющие.

Примечание автора:

Следует отметить, что в системах связи для снижения PAPR применяют DFT-распределение спектра, вводящее в OFDM схему дополнительную пару DFT-IDFT. Во время разработки этот подход требовал больших архитектурных изменений протокола и приводил к снижению помехоустойчивости. Графики BER/PER смещались по оси SNR в область больших значений. Впоследствии на замену этой идее был использован CAF.

def clip_and_filter(
       data: npt.NDArray[np.float64], papr_cutoff_db: float, cutoff_freq_hz: float, fs_hz: float, filter_order: int = 8
) -> npt.NDArray[np.float64]:
    rms = np.sqrt(np.mean(data ** 2))

    clip_threshold = rms * (10 ** (papr_cutoff_db / 20.0))

    clipped_data = np.clip(data, -clip_threshold, clip_threshold)

    nyquist = 0.5 * fs_hz
    normal_cutoff = cutoff_freq_hz / nyquist
    b, a = signal.butter(filter_order, normal_cutoff, btype='low', analog=False)

    filtered_data = signal.filtfilt(b, a, clipped_data)
    return filtered_data

Функция clip_and_filter работает с дискретным вещественным сигналом. Аргумент papr_cutoff_db определяет порог ограничения PAPR в децибелах, относительно среднеквадратичного значения амплитуд сигнала.

Клиппированный сигнал обрабатывается Low-Pass фильтром Баттерворта 8 порядка, подавляющим в сигнале все частоты выше cutoff_freq_hz.

Финальный клиппинг осуществляется с ограничением в 6 дБ и верхней частотой фильтра в 3.0 КГц.

signal_caf = clip_and_filter(
    signal_full,
    papr_cutoff_db=6.0, cutoff_freq_hz=3000,
    fs_hz=self.DEFAULT_SAMPLE_RATE)
...
signal = (signal_caf * 32768).astype(np.int16)

Полученный сигнал с дискретностью 16 бит и частотой дискретизации 12.0 КГц может быть подан на вход звукового тракта трансивера.

Тестирование в симуляторе

Оценка работоспособности протокола и тестирование гипотез производились в симуляторе эфира. На нем моделировались условия AWGN, многолучевое распространение, частотный дрейф, а также BEC и BSC — модели каналов с замираниями.

Моделирование осуществлялось со следующими параметрами:

SNRdb = -6
channelResponse = np.array([1.0, 0, 0.4, 0, 0, 0.2])
bit_flip_prob = 0.001
erasure_prob = 0.02

SNRdb — отношение сигнал/шум в децибелах.
channelResponse — вектор отклика канала на многолучевое распространение (первый элемент задает амплитуду основного сигнала, последующие — отраженных).
bit_flip_prob — вероятность инверсии бит в сигнале.
erasure_prob — вероятность глубокого замирания сигнала.

Код симулятора:

def simulate_channel(
        signal, time_shift: int, freq_shift_hz: float, sample_rate: int,
        frame_size: int = 160
):
    rx_time = np.zeros(len(signal) + time_shift, dtype=float)
    rx_time[time_shift:] = signal

    analytic_signal = hilbert(rx_time)
...

Симуляция частотного дрейфа. Вещественный сигнал переводится в аналитическую форму через преобразование Гильберта. После смещения спектра путем умножения на комплексную экспоненту сигнал возвращается в вещественную форму.

...
    t = np.arange(len(analytic_signal)) / sample_rate

    dynamic_cfo = freq_shift_hz + 5.0 * (t / t[-1]) ** 2
    phase_drift = 2 * np.pi * np.cumsum(dynamic_cfo) / sample_rate
    cfo_vector = np.exp(1j * phase_drift)
    shifted_complex = analytic_signal * cfo_vector

    shifted_real = shifted_complex.real
...

Симуляция многолучевого распространения. Искажение сигнала моделируется сверткой сигнала с вектором отклика. Полученный сигнал складывается с белым шумом, мощность которого рассчитывается на основе значения SNRdb переведенного из децибел в линейный масштаб.

...
    convolved = np.convolve(shifted_real, channelResponse, mode='full')

    signal_power = np.mean(convolved ** 2)
    if signal_power > 0:
        sigma2 = signal_power * 10 ** (-SNRdb / 10)
        noise = np.sqrt(sigma2) * np.random.randn(*convolved.shape)
        result = convolved + noise
    else:
        result = np.random.randn(*convolved.shape) * 10 ** (-SNRdb / 20)

    num_samples = len(result)
    num_frames = int(np.ceil(num_samples / frame_size))
...

Симуляция условий BSC. Для моделирования инверсии бит блоки сигнала умножается на -1. Это действие изменяет фазу на π радиан, что при демодуляции приводит к зеркальной смене бит.

...
    if bit_flip_prob > 0.0 and frame_size > 0:
         flip_decisions = np.random.choice([-1, 1], size=num_frames, p=[bit_flip_prob, 1.0 - bit_flip_prob])
        flip_mask = np.repeat(flip_decisions, frame_size)[:num_samples]
        result *= flip_mask
...

Симуляция условий BEC. Моделирование глубоких замираний в сигнале осуществляется умножением блока сигнала на 0. Образующаяся тишина имитирует полную потерю данных пакета.

...
    if erasure_prob > 0.0 and frame_size > 0:
        erasure_decisions = np.random.choice([0, 1], size=num_frames, p=[erasure_prob, 1.0 - erasure_prob])
        erasure_mask = np.repeat(erasure_decisions, frame_size)[:num_samples]
        result *= erasure_mask
...

Нормализация выходного дискретного сигнала.

Результаты моделирования зависимостей BER (Bit Error Rate) и PER (Packet Error Rate) от уровня SNR для всех типов модуляций приведены на рисунках 23 и 24.

Рисунок 23. График зависимости BER от уровня SNR.

Рисунок 24. График зависимости BER от уровня SNR.

Тестирование производилось с применением алгоритма контроля целостности данных. Метрика BER определяет частоту появления ошибочных бит в декодированных данных. Характеристика PER отражает долю некорректно декодированных пакетов. Согласно стандарту телекоммуникаций, по которому метрика PER не должна превышать 10%, чувствительность протокола находится в окрестности SNR = -7 дБ (более точно -7.5 дБ).

Тестирование в эфире

После проверки работоспособности протокола в симуляторе и определения границ устойчивости к шумам, замираниям и искажениям была проведена серия натурных испытаний.

В качестве передатчика сигнала был использован трансивер ICOM IC-9700 (рисунок 25), подключенный к компьютеру и настроенный в системе как звуковая карта. На первом этапе прием осуществлялся на аналогичный трансивер, удаленный на 14 км и подключенный к направленной антенне с высоким коэффициентом усиления. Принимаемый сигнал записывался в формате WAV встроенными средствами трансивера. Передача велась на любительском диапазоне 2 м, на частоте 144.300 МГц в однополосной модуляции (USB).

Рисунок 25. Любительский трансивер IC-9700.

Первый эксперимент завершился неуспешно, так как на приемнике была включена автоматическая регулировка усиления (AGC), при этом первая реализация протокола не содержала в себе защитной преамбулы для стабилизации АРУ. Анализ результатов показал, что включенная АРУ исказила преамбулу синхронизации. При этом детектор выдавал ложноположительные срабатывания с неестественно большим значением CFO.

По итогам испытаний протокол был доработан. В частности, была переделана модель симулятора, в нее была добавлена симуляция воздействия АРУ и клиппинга сигнала.  Поскольку на начальном этапе под подозрением также были тип модуляции и эквалайзер звука в приемнике, их модели также были добавлены в симулятор. В результате был модифицирован согласованный фильтр ZC-преамбулы, добавлена тональная преамбула для противодействия влиянию АРУ, а также предприняты попытки снижения пик-фактора (PAPR) сигнала.

Последующие испытания осуществлялись с использованием SDR-приемника RX-888 (рисунок 27) под управлением программы SDR Console (рисунок 28). Демодулированный сигнал записывался в WAV. Запись велась как со включенной АРУ, так и без нее — с фиксированным коэффициентом усиления (рисунки 29 и 30).

Тестовая трансляция содержала в себе три сигнала BEACON (R9FEU LO88CA) и две пары пакетов DATA, содержащих строки "    Though this be madness," и " yet there is method in't. " (обе по 27 байт). При этом, для второго пакета применялась QPSK-модуляция.

На рисунке 26 изображены осциллограмма и спектр передаваемого сигнала. В правой половине графика отчетливо заметна разница в длительности BPSK- и QPSK-пакетов.

Рисунок 26. Спектр и осциллограмма передаваемого сигнала.

Рисунок 27. SDR-приемник RX-888.

Рисунок 28. Прием сигнала в SDR Console.

Рисунок 29. Осциллограмма и спектр сигнала с включенной АРУ.

Рисунок 30. Осциллограмма и спектр сигнала с отключенной АРУ.

Из записи с включенной АРУ была успешно декодирована часть пакетов. Детектирование оставшихся пакетов срывалось из-за искажений, вносимых АРУ.

В записи, сделанной с фиксированным усилением, только в двух пакетах были ошибки детектирования. Анализ результатов показал необходимость предварительного усиления сигнала в приемном тракте самого декодера. После внедрения алгоритма предварительного позволило успешно принять все пакеты. Помимо этого, усилитель сделал частично возможным декодирование сигнала с включенным АРУ.

На основе чего было выдвинуто требование об обязательном отключении АРУ приемника.

Ниже приведен лог декодирования записи сигнала.

18:30:28 [   INFO] Packet: Beacon(callsign='R9FEU', qth='LO88CA', mode=<BeaconMode.BEACON: 15>)
18:30:28 [   INFO] Es/N0 carriers: -5.073dB, -2.487dB, 0.038dB, 0.183dB, 0.105dB, -2.126dB, -3.169dB, -2.747dB, -2.081dB, -1.789dB, -1.675dB, -2.011dB, -1.138dB, 1.060dB, 1.333dB, 1.815dBdB
18:30:28 [   INFO] Es/N0: -0.871dB
18:30:28 [   INFO] SNR: -6.892dB
18:30:28 [   INFO] BER: 3.41%
18:30:28 [   INFO] Packet: Beacon(callsign='R9FEU', qth='LO88CA', mode=<BeaconMode.BEACON: 15>)
18:30:28 [   INFO] Es/N0 carriers: -3.940dB, -1.089dB, 2.176dB, 2.472dB, 2.449dB, 0.260dB, -0.733dB, -0.290dB, 0.368dB, 0.664dB, 0.748dB, 0.335dB, 1.150dB, 3.362dB, 3.568dB, 3.919dBdB
18:30:28 [   INFO] Es/N0: 1.370dB
18:30:28 [   INFO] SNR: -4.651dB
18:30:28 [   INFO] BER: 5.77%
18:30:28 [WARNING] Demod failed: head
18:30:28 [   INFO] Packet: Beacon(callsign='R9FEU', qth='LO88CA', mode=<BeaconMode.BEACON: 15>)
18:30:28 [   INFO] Es/N0 carriers: -3.324dB, -0.758dB, 2.651dB, 2.997dB, 2.991dB, 1.071dB, 0.267dB, 0.639dB, 1.120dB, 1.375dB, 1.525dB, 1.193dB, 1.917dB, 3.939dB, 4.112dB, 4.359dBdB
18:30:28 [   INFO] Es/N0: 2.003dB
18:30:28 [   INFO] SNR: -4.018dB
18:30:28 [   INFO] BER: 6.30%
18:30:29 [   INFO] Packet: Data(reserved=123, payload=b'    Though this be madness,')
18:30:29 [   INFO] Es/N0 carriers: -2.827dB, -0.081dB, 2.605dB, 2.760dB, 2.682dB, 0.442dB, -0.332dB, 0.382dB, 0.743dB, 0.779dB, 0.949dB, 0.517dB, 1.295dB, 3.547dB, 3.751dB, 4.138dBdB
18:30:29 [   INFO] Es/N0: 1.674dB
18:30:29 [   INFO] SNR: -4.347dB
18:30:29 [   INFO] BER: 3.69%
18:30:29 [   INFO] Packet: Data(reserved=123, payload=b" yet there is method in't. ")
18:30:29 [   INFO] Es/N0 carriers: -0.259dB, 2.071dB, 5.400dB, 5.646dB, 5.550dB, 3.653dB, 3.127dB, 3.720dB, 3.628dB, 3.335dB, 3.336dB, 2.969dB, 3.806dB, 6.304dB, 6.589dB, 6.845dBdB
18:30:29 [   INFO] Es/N0: 4.459dB
18:30:29 [   INFO] SNR: -1.561dB
18:30:29 [   INFO] BER: 5.42%
18:30:29 [WARNING] Demod failed: head
18:30:29 [   INFO] Packet: Data(reserved=123, payload=b'    Though this be madness,')
18:30:29 [   INFO] Es/N0 carriers: -2.237dB, 0.281dB, 3.202dB, 3.427dB, 3.371dB, 1.420dB, 0.742dB, 1.286dB, 1.563dB, 1.589dB, 1.777dB, 1.441dB, 2.110dB, 4.120dB, 4.299dB, 4.589dBdB
18:30:29 [   INFO] Es/N0: 2.367dB
18:30:29 [   INFO] SNR: -3.653dB
18:30:29 [   INFO] BER: 5.07%
18:30:29 [   INFO] Packet: Data(reserved=123, payload=b" yet there is method in't. ")
18:30:29 [   INFO] Es/N0 carriers: -1.564dB, 2.677dB, 6.098dB, 6.236dB, 6.185dB, 3.370dB, 2.325dB, 3.490dB, 3.658dB, 3.410dB, 3.323dB, 2.048dB, 3.123dB, 6.820dB, 7.059dB, 7.412dBdB
18:30:29 [   INFO] Es/N0: 4.662dB
18:30:29 [   INFO] SNR: -1.358dB
18:30:29 [   INFO] BER: 6.00%

Декодер полностью извлек данные всех пакетов из сигнала, включая пакеты с QPSK-модуляцией. Поскольку все пакеты были приняты, записи "Demod failed: head" свидетельствуют о ложноположительных срабатываниях детектора. Ненулевое значение BER указывает на наличие остаточных искажений в сигнале.

Заключение

Подводя итоги разработки и тестов в реальном эфире, можно утверждать, что реализация любительского протокола связи с использованием OFDM имеет право на существование и доказывает работоспособность даже при отрицательных значениях SNR. Также можно отметить, что OFDM достаточно капризен и его применение на КВ диапазонах может быть крайне затруднительным, даже если увеличить продолжительность символов для когерентного накопления сигнала.

Кроме того, передача OFDM-сигнала поверх SSB — не самое рациональное решение, не позволяющее эффективно использовать выделенную полосу частот. Гораздо оптимальнее формировать дискретный IQ-сигнал с последующим переносом его на несущую частоту. Такой подход накладывает дополнительные аппаратные требования и сильно ограничивает выбор любительских трансиверов.

В основу архитектуры разрабатываемого протокола легли технические решения и принципы, заимствованные из стандартов связи LTE и Wi-Fi. По этой причине, хоть статья и нацелена на радиолюбительскую аудиторию, фактически она выходит далеко за рамки радиолюбительства, затрагивая темы цифровой обработки связей, беспроводных телекоммуникаций и методов коррекции ошибок.

В заключение остается лишь процитировать Шекспира: Though this be madness, yet there is method in't.