Lab043: validate pilot-aided short BPSK link

2026-08-19 18:14:13 +03:00
parent fa2e473c6d
commit f0fa7e8a46
6 changed files with 4416 additions and 0 deletions
--- a/protocol/bpsk_radio.py
+++ b/protocol/bpsk_radio.py
@@ -17,6 +17,7 @@ Lab019 и Lab023 при импорте создают каталоги и сох
 from __future__ import annotations

 import struct
+from dataclasses import dataclass

 import numpy as np

@@ -38,8 +39,48 @@ CFO_REFINEMENT_STEP_HZ = 1.0
 CFO_DEAD_ZONE_HZ = 15.0
 MINIMUM_PHASE_CONSISTENCY = 0.55

+# Критерии рабочей оценки малой остаточной CFO по отдельному известному
+# пилоту Lab043. Они зафиксированы до нового аппаратного опыта.
+KNOWN_PILOT_BLOCK_SYMBOL_COUNT = 20
+MINIMUM_KNOWN_PILOT_BLOCK_COUNT = 8
+MINIMUM_KNOWN_PILOT_COHERENCE = 0.50
+MAXIMUM_KNOWN_PILOT_PHASE_FIT_RMSE_RAD = 0.25
+
 MAXIMUM_PROTOCOL_PACKET_BYTES = 4096

+
+@dataclass(frozen=True)
+class TwoTonePhaseRefinement:
+    """Явный результат фазового уточнения двух известных тонов."""
+
+    valid: bool
+    invalid_reason: str
+    low_frequency_hz: float
+    high_frequency_hz: float
+    carrier_offset_hz: float
+    tone_spacing_hz: float
+    residual_cfo_hz: float
+    phase_fit_rmse_rad: float
+    low_phase_fit_rms_rad: float
+    high_phase_fit_rms_rad: float
+    paired_phase_fit_rms_rad: float
+    block_count: int
+
+
+@dataclass(frozen=True)
+class KnownPilotCarrierEstimate:
+    """Оценка малой остаточной CFO по отдельному известному BPSK-пилоту."""
+
+    valid: bool
+    invalid_reason: str
+    frequency_hz: float
+    phase_increment_rad_per_symbol: float
+    initial_phase_rad: float
+    phase_fit_rmse_rad: float
+    mean_block_coherence: float
+    block_count: int
+    known_symbol_count: int
+
 # Перенесено дословно из Lab023 (tests/lab023_guarded_cfo_correction.py).
 def bytes_to_bits(
    data: bytes,
@@ -136,6 +177,496 @@ def bpsk_demodulate(
    ).astype(np.uint8)


+def find_known_tone_peaks(
+    frequency_axis_hz: np.ndarray,
+    power_spectrum: np.ndarray,
+    expected_offset_hz: float,
+    search_half_width_hz: float,
+) -> dict:
+    """Найти максимумы мощности около двух известных симметричных тонов.
+
+    Функция не задаёт способ оценки шумового фона и порог достоверности:
+    эти решения принадлежат конкретному эксперименту. Если в одном из
+    поисковых окон нет конечных значений мощности, соответствующие частота
+    и мощность возвращаются как ``NaN``.
+    """
+
+    frequencies = np.asarray(frequency_axis_hz, dtype=np.float64)
+    powers = np.asarray(power_spectrum, dtype=np.float64)
+
+    if frequencies.ndim != 1 or powers.ndim != 1:
+        raise ValueError("Ось частот и спектр мощности должны быть одномерными")
+    if len(frequencies) != len(powers):
+        raise ValueError("Ось частот и спектр мощности должны иметь одинаковую длину")
+    if not np.isfinite(expected_offset_hz) or expected_offset_hz <= 0.0:
+        raise ValueError("Ожидаемый отступ тона должен быть положительным")
+    if not np.isfinite(search_half_width_hz) or search_half_width_hz <= 0.0:
+        raise ValueError("Полуширина окна поиска должна быть положительной")
+
+    def peak_in_window(center_hz: float) -> tuple[float, float]:
+        in_window = (
+            (frequencies >= center_hz - search_half_width_hz)
+            & (frequencies <= center_hz + search_half_width_hz)
+            & np.isfinite(frequencies)
+            & np.isfinite(powers)
+            & (powers >= 0.0)
+        )
+        indexes = np.flatnonzero(in_window)
+        if indexes.size == 0:
+            return float("nan"), float("nan")
+        peak_index = int(indexes[np.argmax(powers[indexes])])
+        return float(frequencies[peak_index]), float(powers[peak_index])
+
+    low_frequency_hz, low_power = peak_in_window(-expected_offset_hz)
+    high_frequency_hz, high_power = peak_in_window(expected_offset_hz)
+    return {
+        "low_frequency_hz": low_frequency_hz,
+        "high_frequency_hz": high_frequency_hz,
+        "low_power": low_power,
+        "high_power": high_power,
+    }
+
+
+def estimate_two_tone_offsets(
+    low_frequency_hz: float,
+    high_frequency_hz: float,
+    known_tone_offset_hz: float,
+) -> dict:
+    """Оценить грубую CFO и отношение тактов по двум известным тонам.
+
+    ``clock_scale`` определён как отношение масштаба TX к масштабу RX.
+    Поэтому для перехода принятой последовательности на временную сетку
+    передатчика её ожидаемая длина равна ``round(N * clock_scale)``.
+    Невычислимая оценка представляется только значениями ``NaN``.
+    """
+
+    if not np.isfinite(known_tone_offset_hz) or known_tone_offset_hz <= 0.0:
+        raise ValueError("Известный отступ тона должен быть положительным")
+
+    invalid = {
+        "carrier_offset_hz": float("nan"),
+        "clock_scale": float("nan"),
+        "sample_clock_error_ppm": float("nan"),
+    }
+    if not np.isfinite(low_frequency_hz) or not np.isfinite(high_frequency_hz):
+        return invalid
+    if high_frequency_hz <= low_frequency_hz:
+        return invalid
+
+    carrier_offset_hz = (low_frequency_hz + high_frequency_hz) / 2.0
+    clock_scale = (
+        (high_frequency_hz - low_frequency_hz)
+        / (2.0 * known_tone_offset_hz)
+    )
+    if not np.isfinite(clock_scale) or clock_scale <= 0.0:
+        return invalid
+
+    return {
+        "carrier_offset_hz": float(carrier_offset_hz),
+        "clock_scale": float(clock_scale),
+        "sample_clock_error_ppm": float((clock_scale - 1.0) * 1e6),
+    }
+
+
+def _tone_block_projections(
+    samples: np.ndarray,
+    sample_rate_hz: float,
+    coarse_frequency_hz: float,
+    block_samples: int,
+    hop_samples: int,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Спроецировать фазово-непрерывный сигнал на грубую частоту тона."""
+
+    received = np.asarray(samples, dtype=np.complex128)
+    if received.ndim != 1:
+        raise ValueError("Комплексные отсчёты должны быть одномерными")
+    if not np.isfinite(sample_rate_hz) or sample_rate_hz <= 0.0:
+        raise ValueError("Частота дискретизации должна быть положительной")
+    if not np.isfinite(coarse_frequency_hz):
+        raise ValueError("Грубая частота тона должна быть конечной")
+    if not isinstance(block_samples, (int, np.integer)) or block_samples < 3:
+        raise ValueError("Размер блока должен быть целым числом не меньше трёх")
+    if not isinstance(hop_samples, (int, np.integer)) or hop_samples <= 0:
+        raise ValueError("Шаг блоков должен быть положительным целым числом")
+    if len(received) < block_samples + 2 * hop_samples:
+        raise ValueError("Для фазовой оценки нужны не менее трёх блоков")
+
+    starts = np.arange(
+        0,
+        len(received) - block_samples + 1,
+        hop_samples,
+        dtype=np.int64,
+    )
+    sample_indexes = np.arange(len(received), dtype=np.float64)
+    mixed = received * np.exp(
+        -1j * 2.0 * np.pi * coarse_frequency_hz * sample_indexes / sample_rate_hz
+    )
+    window = np.hanning(block_samples)
+    projections = np.asarray(
+        [
+            np.sum(mixed[start : start + block_samples] * window)
+            for start in starts
+        ],
+        dtype=np.complex128,
+    )
+    center_times_seconds = (
+        starts.astype(np.float64) + (block_samples - 1) / 2.0
+    ) / sample_rate_hz
+    return center_times_seconds, projections
+
+
+def _weighted_phase_frequency(
+    center_times_seconds: np.ndarray,
+    projections: np.ndarray,
+) -> tuple[float, float]:
+    """Оценить наклон развёрнутой фазы и СКО остатка линейной модели."""
+
+    times = np.asarray(center_times_seconds, dtype=np.float64)
+    values = np.asarray(projections, dtype=np.complex128)
+    weights = np.abs(values) ** 2
+    weight_sum = float(np.sum(weights))
+    if (
+        times.ndim != 1
+        or values.ndim != 1
+        or len(times) != len(values)
+        or len(times) < 3
+        or not np.all(np.isfinite(times))
+        or not np.all(np.isfinite(values))
+        or not np.isfinite(weight_sum)
+        or weight_sum <= 0.0
+    ):
+        raise ValueError("Фазовые проекции не позволяют оценить частоту")
+
+    phases = np.unwrap(np.angle(values))
+    mean_time = float(np.sum(weights * times) / weight_sum)
+    mean_phase = float(np.sum(weights * phases) / weight_sum)
+    centered_times = times - mean_time
+    denominator = float(np.sum(weights * centered_times**2))
+    if not np.isfinite(denominator) or denominator <= 0.0:
+        raise ValueError("Временные точки не позволяют оценить наклон фазы")
+
+    slope_rad_per_second = float(
+        np.sum(weights * centered_times * (phases - mean_phase)) / denominator
+    )
+    intercept = mean_phase - slope_rad_per_second * mean_time
+    residuals = phases - (intercept + slope_rad_per_second * times)
+    phase_fit_rms_rad = float(
+        np.sqrt(np.sum(weights * residuals**2) / weight_sum)
+    )
+    return slope_rad_per_second / (2.0 * np.pi), phase_fit_rms_rad
+
+
+def refine_tone_frequency_from_phase(
+    samples: np.ndarray,
+    sample_rate_hz: float,
+    coarse_frequency_hz: float,
+    block_samples: int,
+    hop_samples: int | None = None,
+) -> dict:
+    """Уточнить частоту тона по наклону фазы когерентных проекций.
+
+    Входной участок должен быть фазово-непрерывным. Функцию нельзя применять
+    через границы отдельных аппаратных чтений, если непрерывность их фазы не
+    доказана. Грубая частота должна быть достаточно точной, чтобы остаточная
+    фаза между соседними блоками разворачивалась без неоднозначности.
+    """
+
+    hop = block_samples // 2 if hop_samples is None else hop_samples
+    times, projections = _tone_block_projections(
+        samples,
+        sample_rate_hz,
+        coarse_frequency_hz,
+        block_samples,
+        hop,
+    )
+    residual_frequency_hz, phase_fit_rms_rad = _weighted_phase_frequency(
+        times,
+        projections,
+    )
+    return {
+        "frequency_hz": float(coarse_frequency_hz + residual_frequency_hz),
+        "residual_frequency_hz": float(residual_frequency_hz),
+        "phase_fit_rms_rad": float(phase_fit_rms_rad),
+        "block_count": int(len(projections)),
+    }
+
+
+def refine_two_tone_frequencies_from_phase(
+    samples: np.ndarray,
+    sample_rate_hz: float,
+    low_coarse_frequency_hz: float,
+    high_coarse_frequency_hz: float,
+    block_samples: int,
+    hop_samples: int | None = None,
+) -> TwoTonePhaseRefinement:
+    """Уточнить несущую и расстояние двух тонов на непрерывном участке.
+
+    Несущая получается из двух индивидуальных наклонов фазы. Расстояние тонов
+    оценивается по фазе произведения верхней проекции на сопряжённую нижнюю:
+    общая фазовая ошибка приёмника при этом сокращается. Итоговые частоты
+    строятся из общей оценки центра и парной оценки расстояния.
+    """
+
+    if high_coarse_frequency_hz <= low_coarse_frequency_hz:
+        raise ValueError("Частота верхнего тона должна быть выше частоты нижнего")
+    hop = block_samples // 2 if hop_samples is None else hop_samples
+    times, low_projections = _tone_block_projections(
+        samples,
+        sample_rate_hz,
+        low_coarse_frequency_hz,
+        block_samples,
+        hop,
+    )
+    high_times, high_projections = _tone_block_projections(
+        samples,
+        sample_rate_hz,
+        high_coarse_frequency_hz,
+        block_samples,
+        hop,
+    )
+    if not np.array_equal(times, high_times):
+        raise RuntimeError("Временные сетки двух тонов не совпали")
+
+    low_residual_hz, low_rms_rad = _weighted_phase_frequency(times, low_projections)
+    high_residual_hz, high_rms_rad = _weighted_phase_frequency(times, high_projections)
+    spacing_residual_hz, paired_rms_rad = _weighted_phase_frequency(
+        times,
+        high_projections * np.conj(low_projections),
+    )
+
+    individual_low_hz = low_coarse_frequency_hz + low_residual_hz
+    individual_high_hz = high_coarse_frequency_hz + high_residual_hz
+    center_hz = (individual_low_hz + individual_high_hz) / 2.0
+    spacing_hz = (
+        high_coarse_frequency_hz
+        - low_coarse_frequency_hz
+        + spacing_residual_hz
+    )
+    low_frequency_hz = float(center_hz - spacing_hz / 2.0)
+    high_frequency_hz = float(center_hz + spacing_hz / 2.0)
+    residual_cfo_hz = float(
+        center_hz
+        - (low_coarse_frequency_hz + high_coarse_frequency_hz) / 2.0
+    )
+    phase_fit_rmse_rad = float(
+        np.sqrt(np.mean(np.square((low_rms_rad, high_rms_rad, paired_rms_rad))))
+    )
+    numeric = (
+        low_frequency_hz,
+        high_frequency_hz,
+        center_hz,
+        spacing_hz,
+        residual_cfo_hz,
+        phase_fit_rmse_rad,
+        low_rms_rad,
+        high_rms_rad,
+        paired_rms_rad,
+    )
+    valid = bool(all(np.isfinite(value) for value in numeric))
+    if not valid:
+        nan = float("nan")
+        return TwoTonePhaseRefinement(
+            valid=False,
+            invalid_reason="фазовое уточнение вернуло невычислимое значение",
+            low_frequency_hz=nan,
+            high_frequency_hz=nan,
+            carrier_offset_hz=nan,
+            tone_spacing_hz=nan,
+            residual_cfo_hz=nan,
+            phase_fit_rmse_rad=nan,
+            low_phase_fit_rms_rad=nan,
+            high_phase_fit_rms_rad=nan,
+            paired_phase_fit_rms_rad=nan,
+            block_count=int(len(low_projections)),
+        )
+    return TwoTonePhaseRefinement(
+        valid=True,
+        invalid_reason="",
+        low_frequency_hz=low_frequency_hz,
+        high_frequency_hz=high_frequency_hz,
+        carrier_offset_hz=float(center_hz),
+        tone_spacing_hz=float(spacing_hz),
+        residual_cfo_hz=residual_cfo_hz,
+        phase_fit_rmse_rad=phase_fit_rmse_rad,
+        low_phase_fit_rms_rad=float(low_rms_rad),
+        high_phase_fit_rms_rad=float(high_rms_rad),
+        paired_phase_fit_rms_rad=float(paired_rms_rad),
+        block_count=int(len(low_projections)),
+    )
+
+
+def apply_coarse_frequency_correction(
+    samples: np.ndarray,
+    carrier_offset_hz: float,
+    sample_rate_hz: float,
+) -> np.ndarray:
+    """Убрать измеренный сдвиг несущей из комплексных отсчётов."""
+
+    received = np.asarray(samples, dtype=np.complex128)
+    if received.ndim != 1:
+        raise ValueError("Комплексные отсчёты должны быть одномерными")
+    if not np.isfinite(carrier_offset_hz):
+        raise ValueError("Сдвиг несущей должен быть конечным")
+    if not np.isfinite(sample_rate_hz) or sample_rate_hz <= 0.0:
+        raise ValueError("Частота дискретизации должна быть положительной")
+
+    sample_indexes = np.arange(len(received), dtype=np.float64)
+    correction = np.exp(
+        -1j * 2.0 * np.pi * carrier_offset_hz * sample_indexes / sample_rate_hz
+    )
+    return received * correction
+
+
+def resample_for_clock_scale(
+    samples: np.ndarray,
+    clock_scale: float,
+) -> np.ndarray:
+    """Передискретизировать отсчёты на временную сетку передатчика.
+
+    Используется комплексная линейная интерполяция. Она детерминирована,
+    не требует рационального приближения ppm-отношения и подходит для
+    сильно передискретизированного BPSK-тракта. Направление преобразования
+    закрепляется синтетическими тестами для ошибок обоих знаков.
+    """
+
+    received = np.asarray(samples, dtype=np.complex128)
+    if received.ndim != 1:
+        raise ValueError("Комплексные отсчёты должны быть одномерными")
+    if not np.isfinite(clock_scale) or clock_scale <= 0.0:
+        raise ValueError("Масштаб такта должен быть положительным")
+    if received.size == 0:
+        return received.copy()
+    if received.size == 1:
+        return np.repeat(received, max(1, round(clock_scale)))
+
+    output_length = max(1, round(len(received) * clock_scale))
+    output_indexes = np.arange(output_length, dtype=np.float64)
+    source_positions = np.minimum(output_indexes / clock_scale, len(received) - 1.0)
+    source_indexes = np.arange(len(received), dtype=np.float64)
+    real = np.interp(source_positions, source_indexes, received.real)
+    imaginary = np.interp(source_positions, source_indexes, received.imag)
+    return real + 1j * imaginary
+
+
+def estimate_known_pilot_carrier(
+    received_symbols: np.ndarray,
+    known_symbols: np.ndarray,
+    symbol_rate: float = SYMBOL_RATE,
+    block_symbol_count: int = KNOWN_PILOT_BLOCK_SYMBOL_COUNT,
+    minimum_block_coherence: float = MINIMUM_KNOWN_PILOT_COHERENCE,
+    maximum_phase_fit_rmse_rad: float = MAXIMUM_KNOWN_PILOT_PHASE_FIT_RMSE_RAD,
+) -> KnownPilotCarrierEstimate:
+    """Оценить малую остаточную CFO по отдельному известному BPSK-пилоту.
+
+    Известные знаки удаляются, затем символы объединяются в короткие
+    когерентные блоки. Развёрнутая фаза блоков аппроксимируется устойчивой
+    линейной моделью Huber. Полезная нагрузка для оценки не используется.
+    """
+
+    received = np.asarray(received_symbols, dtype=np.complex128)
+    known = np.asarray(known_symbols, dtype=np.complex128)
+    if received.ndim != 1 or known.ndim != 1:
+        raise ValueError("Принятый и известный пилоты должны быть одномерными")
+    if len(received) != len(known):
+        raise ValueError("Принятый и известный пилоты должны иметь одинаковую длину")
+    if not np.isfinite(symbol_rate) or symbol_rate <= 0.0:
+        raise ValueError("Символьная скорость должна быть положительной")
+    if not isinstance(block_symbol_count, (int, np.integer)) or block_symbol_count < 2:
+        raise ValueError("Размер когерентного блока должен быть целым и не меньше двух")
+    if not 0.0 < minimum_block_coherence <= 1.0:
+        raise ValueError("Порог когерентности должен находиться в интервале (0, 1]")
+    if not np.isfinite(maximum_phase_fit_rmse_rad) or maximum_phase_fit_rmse_rad <= 0.0:
+        raise ValueError("Порог ошибки фазовой модели должен быть положительным")
+    if not np.all(np.isfinite(received)) or not np.all(np.isfinite(known)):
+        raise ValueError("Пилот не должен содержать нечисловые значения")
+    if np.any(np.abs(known) <= 0.0):
+        raise ValueError("Все известные символы пилота должны быть ненулевыми")
+
+    block_count = len(received) // block_symbol_count
+    if block_count < MINIMUM_KNOWN_PILOT_BLOCK_COUNT:
+        raise ValueError(
+            f"Для оценки нужны не менее {MINIMUM_KNOWN_PILOT_BLOCK_COUNT} когерентных блоков"
+        )
+    usable_count = block_count * block_symbol_count
+    despread = received[:usable_count] * np.conj(known[:usable_count])
+    blocks = despread.reshape(block_count, block_symbol_count)
+    block_sums = np.sum(blocks, axis=1)
+    block_magnitude_sums = np.sum(np.abs(blocks), axis=1) + 1e-12
+    block_coherences = np.abs(block_sums) / block_magnitude_sums
+    mean_block_coherence = float(np.mean(block_coherences))
+
+    centers_symbols = (
+        np.arange(block_count, dtype=np.float64) * block_symbol_count
+        + (block_symbol_count - 1) / 2.0
+    )
+    center_times_seconds = centers_symbols / symbol_rate
+    phases = np.unwrap(np.angle(block_sums))
+    base_weights = np.maximum(np.abs(block_sums), 1e-12)
+    design = np.column_stack((np.ones(block_count), center_times_seconds))
+    weights = base_weights.copy()
+    coefficients = np.zeros(2, dtype=np.float64)
+
+    for _iteration in range(12):
+        root_weights = np.sqrt(weights)
+        new_coefficients, *_unused = np.linalg.lstsq(
+            design * root_weights[:, np.newaxis],
+            phases * root_weights,
+            rcond=None,
+        )
+        residuals = phases - design @ new_coefficients
+        residual_median = float(np.median(residuals))
+        robust_scale = (
+            1.4826 * float(np.median(np.abs(residuals - residual_median)))
+            + 1e-9
+        )
+        huber_limit = 1.345 * robust_scale
+        robust_weights = np.ones_like(residuals)
+        outliers = np.abs(residuals) > huber_limit
+        robust_weights[outliers] = huber_limit / np.abs(residuals[outliers])
+        weights = base_weights * robust_weights
+        if np.allclose(coefficients, new_coefficients, rtol=0.0, atol=1e-12):
+            coefficients = new_coefficients
+            break
+        coefficients = new_coefficients
+
+    residuals = phases - design @ coefficients
+    phase_fit_rmse_rad = float(
+        np.sqrt(np.sum(base_weights * residuals**2) / np.sum(base_weights))
+    )
+    frequency_hz = float(coefficients[1] / (2.0 * np.pi))
+    initial_phase_rad = float(coefficients[0])
+    phase_increment = float(2.0 * np.pi * frequency_hz / symbol_rate)
+
+    invalid_reasons: list[str] = []
+    if not all(
+        np.isfinite(value)
+        for value in (frequency_hz, initial_phase_rad, phase_fit_rmse_rad, mean_block_coherence)
+    ):
+        invalid_reasons.append("оценка содержит нечисловое значение")
+    if mean_block_coherence < minimum_block_coherence:
+        invalid_reasons.append(
+            f"когерентность {mean_block_coherence:.4f} ниже порога {minimum_block_coherence:.2f}"
+        )
+    if phase_fit_rmse_rad > maximum_phase_fit_rmse_rad:
+        invalid_reasons.append(
+            f"RMSE фазы {phase_fit_rmse_rad:.4f} выше порога {maximum_phase_fit_rmse_rad:.2f} рад"
+        )
+
+    valid = not invalid_reasons
+    nan = float("nan")
+    return KnownPilotCarrierEstimate(
+        valid=valid,
+        invalid_reason="; ".join(invalid_reasons),
+        frequency_hz=frequency_hz if valid else nan,
+        phase_increment_rad_per_symbol=phase_increment if valid else nan,
+        initial_phase_rad=initial_phase_rad if valid else nan,
+        phase_fit_rmse_rad=phase_fit_rmse_rad,
+        mean_block_coherence=mean_block_coherence,
+        block_count=block_count,
+        known_symbol_count=len(known),
+    )
+
+
 # Перенесено дословно из Lab023 (tests/lab023_guarded_cfo_correction.py).
 def estimate_carrier_parameters(
    received_symbols: np.ndarray,