Lab045: demonstrate real RF BPSK error correction

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

protocol/convolutional_fec.py

@@ -0,0 +1,216 @@
+"""Свёрточный код R=1/2, K=7 и soft-decision декодер Витерби.
+
+Настоящий помехоустойчивый код уровня бит для исправления неизвестных
+битовых ошибок реального BPSK-канала (Lab044 показал BER ≈0.5–0.6%, при
+котором незащищённый пакет практически не проходит CRC). В отличие от
+``protocol/packet_erasure_fec`` (GF(256), только стирания), этот код
+исправляет случайные битовые ошибки без знания их позиций.
+
+Параметры зафиксированы и не передаются в каждом кадре:
+
+- constraint length ``K = 7``; память ``m = 6``; число состояний ``2**6 = 64``;
+- скорость ``R = 1/2`` (два выходных бита на один входной);
+- порождающие многочлены ``G0 = 0o171``, ``G1 = 0o133`` (стандарт NASA/Voyager);
+- начальное состояние регистра ``0``;
+- терминирование: к информационным битам добавляются ``m = 6`` нулевых
+  хвостовых бита, так что финальное состояние тоже ``0``.
+
+Соглашение сдвигового регистра. Полный 7-битный регистр для вычисления
+выхода собирается как ``full = (input_bit << 6) | state``, где ``state`` —
+6 предыдущих входных битов, самый свежий в бите 5, самый старый в бите 0.
+Порождающий многочлен применяется как битовая маска этого 7-битного
+регистра, выходной бит = чётность числа единиц. Новое состояние =
+``full >> 1`` (самый старый бит вытесняется). Такое соглашение проверяется
+независимо посчитанными эталонными векторами в
+``tests/test_lab045_hackrf_bpsk_fec.py``.
+
+Соответствие BPSK (совпадает с ``protocol/bpsk_radio.bpsk_modulate``):
+бит ``0`` → символ ``-1``, бит ``1`` → символ ``+1``. Декодер работает с
+soft-значениями BPSK (проекция принятого символа на решающую ось,
+``normalized.real`` в Lab044) и использует евклидову ветвевую метрику к
+ожидаемым символам ``{-1, +1}``. Hard-decision — частный случай: знаки
+soft-значений.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+CONSTRAINT_LENGTH = 7
+MEMORY = CONSTRAINT_LENGTH - 1  # 6
+STATE_COUNT = 1 << MEMORY  # 64
+GENERATOR_POLYNOMIALS_OCTAL = (0o171, 0o133)  # (171, 133) octal
+OUTPUTS_PER_INPUT = len(GENERATOR_POLYNOMIALS_OCTAL)  # 2
+TAIL_BIT_COUNT = MEMORY  # 6
+CODE_RATE = (1, OUTPUTS_PER_INPUT)  # 1/2
+
+_STATE_MASK = STATE_COUNT - 1  # 0x3F
+_INPUT_BIT_POSITION = MEMORY  # входной бит занимает бит 6 полного регистра
+_LOWER_STATE_MASK = (1 << (MEMORY - 1)) - 1  # 0x1F, младшие 5 бит состояния
+
+
+def _parity(value: int) -> int:
+    """Чётность числа установленных битов (0 или 1)."""
+
+    return bin(value).count("1") & 1
+
+
+def _next_state(state: int, input_bit: int) -> int:
+    full = (input_bit << _INPUT_BIT_POSITION) | state
+    return full >> 1
+
+
+def _output_bits(state: int, input_bit: int) -> tuple[int, ...]:
+    full = (input_bit << _INPUT_BIT_POSITION) | state
+    return tuple(_parity(full & generator) for generator in GENERATOR_POLYNOMIALS_OCTAL)
+
+
+def _build_decode_tables() -> tuple[np.ndarray, np.ndarray]:
+    """Предвычислить предшественников и ожидаемые символы для решётки.
+
+    Для каждого следующего состояния ``next_state`` вход, приведший в него,
+    однозначен и равен старшему биту состояния (бит 5). У ``next_state``
+    ровно два предшественника, различающиеся вытесненным младшим битом.
+    Возвращает ``predecessors[next_state, branch]`` и
+    ``symbols[next_state, branch, output_index]`` со значениями в {-1, +1}.
+    """
+
+    predecessors = np.zeros((STATE_COUNT, 2), dtype=np.int64)
+    symbols = np.zeros((STATE_COUNT, 2, OUTPUTS_PER_INPUT), dtype=np.float64)
+    for next_state in range(STATE_COUNT):
+        input_bit = next_state >> (MEMORY - 1)  # старший бит состояния
+        base = (next_state & _LOWER_STATE_MASK) << 1
+        for branch in range(2):
+            previous_state = base | branch
+            assert _next_state(previous_state, input_bit) == next_state
+            predecessors[next_state, branch] = previous_state
+            output = _output_bits(previous_state, input_bit)
+            for output_index, bit in enumerate(output):
+                symbols[next_state, branch, output_index] = 2.0 * bit - 1.0
+    return predecessors, symbols
+
+
+_PREDECESSORS, _BRANCH_SYMBOLS = _build_decode_tables()
+_BRANCH_SYMBOL_0 = _BRANCH_SYMBOLS[:, :, 0]
+_BRANCH_SYMBOL_1 = _BRANCH_SYMBOLS[:, :, 1]
+
+
+def _validate_information_bits(information_bits: np.ndarray) -> np.ndarray:
+    raw_bits = np.asarray(information_bits)
+    if raw_bits.ndim != 1:
+        raise ValueError("information_bits должен быть одномерным")
+    if raw_bits.size == 0:
+        raise ValueError("information_bits не должен быть пустым")
+    if not np.all((raw_bits == 0) | (raw_bits == 1)):
+        raise ValueError("information_bits должен содержать только 0 и 1")
+    return raw_bits.astype(np.uint8, copy=False)
+
+
+def encoded_bit_count(information_bit_count: int) -> int:
+    """Число кодовых битов для заданного числа информационных битов."""
+
+    if not isinstance(information_bit_count, (int, np.integer)):
+        raise TypeError("information_bit_count должен быть целым")
+    if information_bit_count < 0:
+        raise ValueError("information_bit_count не может быть отрицательным")
+    return OUTPUTS_PER_INPUT * (information_bit_count + TAIL_BIT_COUNT)
+
+
+def encode(information_bits: np.ndarray) -> np.ndarray:
+    """Свёрточно закодировать биты, добавив внутреннее терминирование.
+
+    Вход — одномерные биты {0,1}. К ним добавляются ``TAIL_BIT_COUNT``
+    нулевых хвостовых бита, регистр стартует из нуля и заканчивает в нуле.
+    Выход — одномерные кодовые биты {0,1} длиной ``2*(len+6)``.
+    """
+
+    bits = _validate_information_bits(information_bits)
+    terminated = np.concatenate((bits, np.zeros(TAIL_BIT_COUNT, dtype=np.uint8)))
+    output = np.empty(OUTPUTS_PER_INPUT * len(terminated), dtype=np.uint8)
+    state = 0
+    for index, input_bit in enumerate(int(bit) for bit in terminated):
+        full = (input_bit << _INPUT_BIT_POSITION) | state
+        for output_index, generator in enumerate(GENERATOR_POLYNOMIALS_OCTAL):
+            output[OUTPUTS_PER_INPUT * index + output_index] = _parity(full & generator)
+        state = full >> 1
+    assert state == 0
+    return output
+
+
+def viterbi_decode_soft(
+    soft_values: np.ndarray,
+    information_bit_count: int,
+) -> np.ndarray:
+    """Декодировать soft-значения BPSK терминированным Витерби.
+
+    ``soft_values`` — одномерный массив вещественных решающих значений
+    (в идеале ±1, реально ``normalized.real``), длиной ровно
+    ``encoded_bit_count(information_bit_count)``. Метрика ветви — евклидово
+    расстояние до ожидаемых символов {-1, +1}. Решётка терминирована:
+    старт и финиш в состоянии 0, хвостовые биты отбрасываются. Возвращает
+    ``information_bit_count`` информационных битов. TX-данные не используются.
+    """
+
+    if not isinstance(information_bit_count, (int, np.integer)):
+        raise TypeError("information_bit_count должен быть целым")
+    if information_bit_count <= 0:
+        raise ValueError("information_bit_count должен быть положительным")
+    information_bit_count = int(information_bit_count)
+    values = np.asarray(soft_values, dtype=np.float64)
+    if values.ndim != 1:
+        raise ValueError("soft_values должен быть одномерным")
+    if not np.all(np.isfinite(values)):
+        raise ValueError("soft_values должен содержать только конечные значения")
+    expected_length = encoded_bit_count(information_bit_count)
+    if len(values) != expected_length:
+        raise ValueError(
+            f"soft_values имеет длину {len(values)}, ожидалось {expected_length}"
+        )
+
+    step_count = information_bit_count + TAIL_BIT_COUNT
+    received = values.reshape(step_count, OUTPUTS_PER_INPUT)
+
+    infinity = np.float64(1e18)
+    path_metric = np.full(STATE_COUNT, infinity, dtype=np.float64)
+    path_metric[0] = 0.0
+    branch_choice = np.empty((step_count, STATE_COUNT), dtype=np.uint8)
+
+    for step in range(step_count):
+        received_0 = received[step, 0]
+        received_1 = received[step, 1]
+        branch_metric = (
+            (received_0 - _BRANCH_SYMBOL_0) ** 2
+            + (received_1 - _BRANCH_SYMBOL_1) ** 2
+        )  # (STATE_COUNT, 2)
+        candidate = path_metric[_PREDECESSORS] + branch_metric  # (STATE_COUNT, 2)
+        chosen = np.argmin(candidate, axis=1).astype(np.uint8)
+        path_metric = candidate[np.arange(STATE_COUNT), chosen]
+        branch_choice[step] = chosen
+
+    # Терминированный трекбэк из состояния 0.
+    decoded = np.empty(step_count, dtype=np.uint8)
+    current_state = 0
+    for step in range(step_count - 1, -1, -1):
+        branch = branch_choice[step, current_state]
+        decoded[step] = current_state >> (MEMORY - 1)
+        current_state = int(_PREDECESSORS[current_state, branch])
+    assert current_state == 0
+    return decoded[:information_bit_count]
+
+
+def viterbi_decode_hard(
+    coded_bits: np.ndarray,
+    information_bit_count: int,
+) -> np.ndarray:
+    """Hard-decision вариант: биты {0,1} отображаются в ±1 и декодируются.
+
+    Служит базовой линией для честного сравнения soft и hard в тестах.
+    """
+
+    raw_bits = np.asarray(coded_bits)
+    if raw_bits.ndim != 1:
+        raise ValueError("coded_bits должен быть одномерным")
+    if not np.all((raw_bits == 0) | (raw_bits == 1)):
+        raise ValueError("coded_bits должен содержать только 0 и 1")
+    bits = raw_bits.astype(np.uint8, copy=False)
+    return viterbi_decode_soft(2.0 * bits.astype(np.float64) - 1.0, information_bit_count)