Lab061: integrate T2 into the strict receiver
This commit is contained in:
LittleSam129
230
protocol/bpsk_carrier_tracking.py
Normal file
230
protocol/bpsk_carrier_tracking.py
Normal file
@@ -0,0 +1,230 @@
|
||||
"""Причинное сопровождение несущей BPSK после известной преамбулы.
|
||||
|
||||
Модуль не разбирает заголовок, полезную нагрузку, FEC или CRC. Начальное
|
||||
состояние P3 оценивается только по переданным в функцию символам преамбулы.
|
||||
После неё T2 получает текущий комплексный символ, собственное решение и своё
|
||||
предыдущее состояние.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
import math
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
PREAMBLE_SYMBOLS = 1280
|
||||
P3_BLOCK_SYMBOLS = 20
|
||||
NORMALIZED_LOOP_BANDWIDTH = 0.0005
|
||||
DAMPING_FACTOR = 1.0 / math.sqrt(2.0)
|
||||
DETECTOR_CLIP_RAD = math.pi / 3.0
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class P3CarrierEstimate:
|
||||
"""Начальное состояние в физических единицах при времени символа 0."""
|
||||
|
||||
frequency_hz: float
|
||||
gain: complex
|
||||
phase_fit_rmse_rad: float
|
||||
squared_coherence: float
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class T2Parameters:
|
||||
"""Коэффициенты петли на один символ."""
|
||||
|
||||
normalized_bandwidth: float
|
||||
damping_factor: float
|
||||
phase_coefficient: float
|
||||
frequency_coefficient: float
|
||||
detector_clip_rad: float
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class T2TrackingResult:
|
||||
"""Применённое состояние и решения для одного последовательного прохода."""
|
||||
|
||||
corrected: np.ndarray
|
||||
hard_decisions: np.ndarray
|
||||
phase_rad: np.ndarray
|
||||
frequency_hz: np.ndarray
|
||||
detector_error_rad: np.ndarray
|
||||
initial_estimate: P3CarrierEstimate
|
||||
parameters: T2Parameters
|
||||
|
||||
|
||||
def t2_parameters(
|
||||
normalized_bandwidth: float = NORMALIZED_LOOP_BANDWIDTH,
|
||||
damping_factor: float = DAMPING_FACTOR,
|
||||
) -> T2Parameters:
|
||||
"""Перевести нормированную полосу в коэффициенты Lab060 на символ."""
|
||||
|
||||
if not np.isfinite(normalized_bandwidth) or normalized_bandwidth <= 0.0:
|
||||
raise ValueError("Нормированная полоса должна быть положительной")
|
||||
if not np.isfinite(damping_factor) or damping_factor <= 0.0:
|
||||
raise ValueError("Коэффициент затухания должен быть положительным")
|
||||
theta = normalized_bandwidth / (
|
||||
damping_factor + 0.25 / damping_factor
|
||||
)
|
||||
denominator = 1.0 + 2.0 * damping_factor * theta + theta * theta
|
||||
return T2Parameters(
|
||||
normalized_bandwidth=float(normalized_bandwidth),
|
||||
damping_factor=float(damping_factor),
|
||||
phase_coefficient=float(4.0 * damping_factor * theta / denominator),
|
||||
frequency_coefficient=float(4.0 * theta * theta / denominator),
|
||||
detector_clip_rad=DETECTOR_CLIP_RAD,
|
||||
)
|
||||
|
||||
|
||||
def _huber_phase_line(
|
||||
block_values: np.ndarray,
|
||||
block_times_seconds: np.ndarray,
|
||||
) -> tuple[float, float, float, float]:
|
||||
"""Повторить зафиксированную в P3 устойчивую прямую фазы."""
|
||||
|
||||
values = np.asarray(block_values, dtype=np.complex128)
|
||||
times = np.asarray(block_times_seconds, dtype=np.float64)
|
||||
if values.ndim != 1 or times.ndim != 1 or len(values) != len(times) or len(values) < 3:
|
||||
raise ValueError("Недостаточно фазовых блоков")
|
||||
phases = np.unwrap(np.angle(values))
|
||||
base_weights = np.maximum(np.abs(values), 1e-12)
|
||||
center_time = float(np.sum(base_weights * times) / np.sum(base_weights))
|
||||
centered = times - center_time
|
||||
design = np.column_stack((np.ones(len(times)), centered))
|
||||
weights = base_weights.copy()
|
||||
coefficients = np.zeros(2, dtype=np.float64)
|
||||
for _iteration in range(12):
|
||||
root = np.sqrt(weights / max(float(np.max(weights)), 1e-30))
|
||||
updated, *_unused = np.linalg.lstsq(
|
||||
design * root[:, None], phases * root, rcond=None
|
||||
)
|
||||
residual = phases - design @ updated
|
||||
median = float(np.median(residual))
|
||||
scale = 1.4826 * float(np.median(np.abs(residual - median))) + 1e-9
|
||||
limit = 1.345 * scale
|
||||
robust = np.ones_like(residual)
|
||||
outside = np.abs(residual) > limit
|
||||
robust[outside] = limit / np.abs(residual[outside])
|
||||
weights = base_weights * robust
|
||||
if np.allclose(updated, coefficients, rtol=0.0, atol=1e-12):
|
||||
coefficients = updated
|
||||
break
|
||||
coefficients = updated
|
||||
residual = phases - design @ coefficients
|
||||
rmse = float(np.sqrt(np.sum(base_weights * residual**2) / np.sum(base_weights)))
|
||||
phase_at_zero = float(coefficients[0] - coefficients[1] * center_time)
|
||||
coherence = float(abs(np.sum(values)) / max(float(np.sum(abs(values))), 1e-30))
|
||||
return phase_at_zero, float(coefficients[1]), rmse, coherence
|
||||
|
||||
|
||||
def estimate_p3(
|
||||
preamble_symbols: np.ndarray,
|
||||
known_preamble_symbols: np.ndarray,
|
||||
symbol_rate: float,
|
||||
) -> P3CarrierEstimate:
|
||||
"""Оценить P3 только по полной 1280-символьной преамбуле."""
|
||||
|
||||
received = np.asarray(preamble_symbols, dtype=np.complex128)
|
||||
known = np.asarray(known_preamble_symbols, dtype=np.complex128)
|
||||
if received.shape != (PREAMBLE_SYMBOLS,) or known.shape != received.shape:
|
||||
raise ValueError("P3 требуется ровно 1280 принятых и известных символов")
|
||||
if not np.isfinite(symbol_rate) or symbol_rate <= 0.0:
|
||||
raise ValueError("Скорость символов должна быть положительной")
|
||||
if not np.all(np.isfinite(received)) or not np.all(np.isfinite(known)):
|
||||
raise ValueError("P3 не принимает NaN или Inf")
|
||||
|
||||
indexes = np.arange(PREAMBLE_SYMBOLS, dtype=np.int64)
|
||||
block_count = PREAMBLE_SYMBOLS // P3_BLOCK_SYMBOLS
|
||||
block_values = np.mean(
|
||||
(received**2).reshape(block_count, P3_BLOCK_SYMBOLS), axis=1
|
||||
)
|
||||
block_times = np.mean(
|
||||
((indexes + 0.5) / symbol_rate).reshape(block_count, P3_BLOCK_SYMBOLS),
|
||||
axis=1,
|
||||
)
|
||||
_phase0, squared_slope, rmse, coherence = _huber_phase_line(
|
||||
block_values, block_times
|
||||
)
|
||||
frequency_hz = float(squared_slope / (4.0 * math.pi))
|
||||
times = (indexes + 0.5) / symbol_rate
|
||||
gain = complex(np.mean(
|
||||
received
|
||||
* np.conj(known)
|
||||
* np.exp(-2j * math.pi * frequency_hz * times)
|
||||
))
|
||||
if not np.isfinite(gain.real) or not np.isfinite(gain.imag) or abs(gain) == 0.0:
|
||||
raise ValueError("P3 получил нулевой или неконечный коэффициент")
|
||||
return P3CarrierEstimate(frequency_hz, gain, rmse, coherence)
|
||||
|
||||
|
||||
def normalized_phase_detector(
|
||||
corrected_symbol: complex,
|
||||
detector_clip_rad: float = DETECTOR_CLIP_RAD,
|
||||
) -> tuple[int, float]:
|
||||
"""Вернуть собственное решение BPSK и ограниченную фазовую ошибку."""
|
||||
|
||||
value = complex(corrected_symbol)
|
||||
if not np.isfinite(value.real) or not np.isfinite(value.imag):
|
||||
raise ValueError("Фазовый детектор не принимает NaN или Inf")
|
||||
hard = 1 if value.real >= 0.0 else -1
|
||||
error = math.atan2(value.imag * hard, value.real * hard)
|
||||
return hard, float(np.clip(error, -detector_clip_rad, detector_clip_rad))
|
||||
|
||||
|
||||
def track_t2(
|
||||
symbols: np.ndarray,
|
||||
symbol_rate: float,
|
||||
initial_estimate: P3CarrierEstimate,
|
||||
*,
|
||||
start_index: int = PREAMBLE_SYMBOLS,
|
||||
parameters: T2Parameters | None = None,
|
||||
) -> T2TrackingResult:
|
||||
"""Последовательно применить T2; будущие символы не используются."""
|
||||
|
||||
values = np.asarray(symbols, dtype=np.complex128)
|
||||
if values.ndim != 1 or not len(values) or not np.all(np.isfinite(values)):
|
||||
raise ValueError("T2 требуется конечная одномерная последовательность")
|
||||
if not np.isfinite(symbol_rate) or symbol_rate <= 0.0:
|
||||
raise ValueError("Скорость символов должна быть положительной")
|
||||
if not 0 <= start_index <= len(values):
|
||||
raise ValueError("Некорректная граница начала сопровождения")
|
||||
if abs(initial_estimate.gain) == 0.0:
|
||||
raise ValueError("Начальный коэффициент P3 не может быть нулевым")
|
||||
selected = t2_parameters() if parameters is None else parameters
|
||||
|
||||
corrected = np.empty(len(values), dtype=np.complex128)
|
||||
decisions = np.empty(len(values), dtype=np.int8)
|
||||
phases = np.empty(len(values), dtype=np.float64)
|
||||
frequencies = np.empty(len(values), dtype=np.float64)
|
||||
detector = np.empty(len(values), dtype=np.float64)
|
||||
omega = 2.0 * math.pi * initial_estimate.frequency_hz / symbol_rate
|
||||
phase0 = float(np.angle(initial_estimate.gain))
|
||||
phase = phase0 + omega * (start_index + 0.5)
|
||||
amplitude = abs(initial_estimate.gain)
|
||||
|
||||
for index, symbol in enumerate(values):
|
||||
current_phase = phase0 + omega * (index + 0.5) if index < start_index else phase
|
||||
value = symbol * np.exp(-1j * current_phase) / amplitude
|
||||
hard, error = normalized_phase_detector(value, selected.detector_clip_rad)
|
||||
corrected[index] = value
|
||||
decisions[index] = hard
|
||||
phases[index] = current_phase
|
||||
frequencies[index] = omega * symbol_rate / (2.0 * math.pi)
|
||||
detector[index] = error
|
||||
if index >= start_index:
|
||||
phase = phase + omega + selected.phase_coefficient * error
|
||||
omega = omega + selected.frequency_coefficient * error
|
||||
|
||||
for array in (corrected, decisions, phases, frequencies, detector):
|
||||
array.flags.writeable = False
|
||||
return T2TrackingResult(
|
||||
corrected=corrected,
|
||||
hard_decisions=decisions,
|
||||
phase_rad=phases,
|
||||
frequency_hz=frequencies,
|
||||
detector_error_rad=detector,
|
||||
initial_estimate=initial_estimate,
|
||||
parameters=selected,
|
||||
)
|
||||
Reference in New Issue
Block a user