1637 lines
36 KiB
Python
1637 lines
36 KiB
Python
"""
|
||
Lab023. Защищённая компенсация частотной ошибки.
|
||
|
||
Сравниваются три стратегии:
|
||
|
||
1. PHASE_ONLY
|
||
Компенсируется только постоянный фазовый поворот.
|
||
|
||
2. ALWAYS_CFO
|
||
Оценённое частотное рассогласование применяется всегда.
|
||
|
||
3. GUARDED_CFO
|
||
Частотная коррекция применяется только тогда, когда:
|
||
|
||
- модуль оценённого CFO превышает мёртвую зону;
|
||
- межсимвольная фазовая оценка достаточно достоверна.
|
||
|
||
Эксперимент выполняется на полном 320-битном радиокадре:
|
||
|
||
PREAMBLE
|
||
+ RADIO SYNC
|
||
+ LENGTH
|
||
+ SDR Rover Link packet
|
||
|
||
В каждой попытке добавляются:
|
||
|
||
- случайная начальная фаза;
|
||
- случайный CFO;
|
||
- AWGN-шум.
|
||
|
||
Символьная синхронизация в этой лабораторной считается
|
||
уже выполненной. Это позволяет исследовать именно логику
|
||
управления CFO-компенсацией.
|
||
"""
|
||
|
||
from csv import DictWriter
|
||
from pathlib import Path
|
||
import struct
|
||
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
|
||
from protocol.packet import (
|
||
CRCError,
|
||
MESSAGE_TYPE_TEXT,
|
||
PacketError,
|
||
build_packet,
|
||
parse_packet,
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Параметры радиокадра
|
||
# ============================================================
|
||
|
||
RADIO_SYNC_WORD = 0xD391
|
||
|
||
PREAMBLE_BIT_COUNT = 64
|
||
|
||
SYMBOL_RATE = 20_000
|
||
|
||
EXPECTED_MESSAGE = "ПРИВЕТ SDR"
|
||
EXPECTED_SEQUENCE_NUMBER = 18
|
||
|
||
|
||
# ============================================================
|
||
# Параметры защищённой CFO-коррекции
|
||
# ============================================================
|
||
|
||
# Если модуль оценки меньше этого значения,
|
||
# приёмник считает CFO практически нулевым.
|
||
CFO_DEAD_ZONE_HZ = 15.0
|
||
|
||
# Достоверность фазового приращения:
|
||
#
|
||
# 0.0 — полностью случайная фаза;
|
||
# 1.0 — идеально постоянное межсимвольное вращение.
|
||
MINIMUM_PHASE_CONSISTENCY = 0.55
|
||
# Уточняющий поиск около грубой оценки CFO.
|
||
#
|
||
# Грубая оценка даёт приблизительный центр,
|
||
# после чего перебираются частоты вокруг него.
|
||
CFO_REFINEMENT_HALF_WIDTH_HZ = 200.0
|
||
|
||
CFO_REFINEMENT_STEP_HZ = 1.0
|
||
|
||
# ============================================================
|
||
# Параметры эксперимента
|
||
# ============================================================
|
||
|
||
EB_N0_VALUES_DB = [
|
||
6.0,
|
||
7.0,
|
||
8.0,
|
||
10.0,
|
||
]
|
||
|
||
# В каждой попытке CFO выбирается случайно:
|
||
#
|
||
# -MAX_CFO ... +MAX_CFO
|
||
MAX_ABS_CFO_VALUES_HZ = [
|
||
0.0,
|
||
250.0,
|
||
1000.0,
|
||
]
|
||
|
||
TRIALS_PER_POINT = 200
|
||
|
||
RANDOM_SEED = 2026
|
||
|
||
|
||
# ============================================================
|
||
# Выходные файлы
|
||
# ============================================================
|
||
|
||
OUTPUT_DIRECTORY = Path(
|
||
"data/processed/lab023"
|
||
)
|
||
|
||
OUTPUT_DIRECTORY.mkdir(
|
||
parents=True,
|
||
exist_ok=True,
|
||
)
|
||
|
||
CSV_PATH = (
|
||
OUTPUT_DIRECTORY
|
||
/ "lab023_guarded_cfo_results.csv"
|
||
)
|
||
|
||
GRAPH_PATH = (
|
||
OUTPUT_DIRECTORY
|
||
/ "lab023_guarded_cfo.png"
|
||
)
|
||
|
||
REPORT_PATH = (
|
||
OUTPUT_DIRECTORY
|
||
/ "lab023_guarded_cfo_report.txt"
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Статусы декодирования
|
||
# ============================================================
|
||
|
||
STATUS_SUCCESS = "SUCCESS"
|
||
STATUS_HEADER_ERROR = "HEADER ERROR"
|
||
STATUS_CRC_ERROR = "CRC ERROR"
|
||
STATUS_PACKET_ERROR = "PACKET ERROR"
|
||
|
||
|
||
# ============================================================
|
||
# Преобразование bytes и bits
|
||
# ============================================================
|
||
|
||
def bytes_to_bits(
|
||
data: bytes,
|
||
) -> np.ndarray:
|
||
"""
|
||
Преобразовать bytes в одномерный массив битов.
|
||
"""
|
||
|
||
if not isinstance(data, bytes):
|
||
raise TypeError(
|
||
"data должен иметь тип bytes"
|
||
)
|
||
|
||
return np.unpackbits(
|
||
np.frombuffer(
|
||
data,
|
||
dtype=np.uint8,
|
||
)
|
||
)
|
||
|
||
|
||
def bits_to_bytes(
|
||
bits: np.ndarray,
|
||
) -> bytes:
|
||
"""
|
||
Упаковать массив битов обратно в bytes.
|
||
"""
|
||
|
||
bits = np.asarray(
|
||
bits,
|
||
dtype=np.uint8,
|
||
)
|
||
|
||
if bits.ndim != 1:
|
||
raise ValueError(
|
||
"bits должен быть одномерным массивом"
|
||
)
|
||
|
||
if len(bits) % 8 != 0:
|
||
raise ValueError(
|
||
"Количество битов должно быть кратно восьми"
|
||
)
|
||
|
||
if not np.all(
|
||
(bits == 0) | (bits == 1)
|
||
):
|
||
raise ValueError(
|
||
"bits должен содержать только 0 и 1"
|
||
)
|
||
|
||
return np.packbits(
|
||
bits
|
||
).tobytes()
|
||
|
||
|
||
# ============================================================
|
||
# BPSK
|
||
# ============================================================
|
||
|
||
def bpsk_modulate(
|
||
bits: np.ndarray,
|
||
) -> np.ndarray:
|
||
"""
|
||
Преобразовать биты в BPSK-символы.
|
||
|
||
0 → -1
|
||
1 → +1
|
||
"""
|
||
|
||
bits = np.asarray(
|
||
bits,
|
||
dtype=np.uint8,
|
||
)
|
||
|
||
symbols = (
|
||
2.0
|
||
* bits.astype(np.float64)
|
||
- 1.0
|
||
)
|
||
|
||
return symbols.astype(
|
||
np.complex128
|
||
)
|
||
|
||
|
||
def bpsk_demodulate(
|
||
symbols: np.ndarray,
|
||
) -> np.ndarray:
|
||
"""
|
||
Демодулировать BPSK по знаку компоненты I.
|
||
"""
|
||
|
||
return (
|
||
symbols.real >= 0.0
|
||
).astype(np.uint8)
|
||
|
||
|
||
# ============================================================
|
||
# Формирование полного радиокадра
|
||
# ============================================================
|
||
|
||
def build_radio_frame(
|
||
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
||
"""
|
||
Сформировать:
|
||
|
||
- биты полного радиокадра;
|
||
- BPSK-символы полного радиокадра;
|
||
- известные символы PREAMBLE + RADIO SYNC.
|
||
"""
|
||
|
||
payload = EXPECTED_MESSAGE.encode(
|
||
"utf-8"
|
||
)
|
||
|
||
protocol_packet = build_packet(
|
||
payload=payload,
|
||
message_type=MESSAGE_TYPE_TEXT,
|
||
sequence_number=EXPECTED_SEQUENCE_NUMBER,
|
||
)
|
||
|
||
if len(protocol_packet) > 65535:
|
||
raise ValueError(
|
||
"Внутренний пакет слишком велик"
|
||
)
|
||
|
||
radio_header = struct.pack(
|
||
">HH",
|
||
RADIO_SYNC_WORD,
|
||
len(protocol_packet),
|
||
)
|
||
|
||
preamble_bits = np.tile(
|
||
np.array(
|
||
[1, 0],
|
||
dtype=np.uint8,
|
||
),
|
||
PREAMBLE_BIT_COUNT // 2,
|
||
)
|
||
|
||
radio_header_bits = bytes_to_bits(
|
||
radio_header
|
||
)
|
||
|
||
protocol_packet_bits = bytes_to_bits(
|
||
protocol_packet
|
||
)
|
||
|
||
frame_bits = np.concatenate(
|
||
[
|
||
preamble_bits,
|
||
radio_header_bits,
|
||
protocol_packet_bits,
|
||
]
|
||
)
|
||
|
||
frame_symbols = bpsk_modulate(
|
||
frame_bits
|
||
)
|
||
|
||
marker_bit_count = (
|
||
PREAMBLE_BIT_COUNT
|
||
+ 16
|
||
)
|
||
|
||
marker_symbols = frame_symbols[
|
||
:marker_bit_count
|
||
]
|
||
|
||
return (
|
||
frame_bits,
|
||
frame_symbols,
|
||
marker_symbols,
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Искажения канала
|
||
# ============================================================
|
||
|
||
def apply_carrier_impairments(
|
||
symbols: np.ndarray,
|
||
frequency_offset_hz: float,
|
||
initial_phase_radians: float,
|
||
) -> np.ndarray:
|
||
"""
|
||
Добавить начальный фазовый поворот и CFO.
|
||
|
||
Между соседними символами фаза изменяется на:
|
||
|
||
2π · CFO / SYMBOL_RATE
|
||
"""
|
||
|
||
symbol_indexes = np.arange(
|
||
len(symbols),
|
||
dtype=np.float64,
|
||
)
|
||
|
||
phase_increment = (
|
||
2.0
|
||
* np.pi
|
||
* frequency_offset_hz
|
||
/ SYMBOL_RATE
|
||
)
|
||
|
||
phase_values = (
|
||
initial_phase_radians
|
||
+ phase_increment
|
||
* symbol_indexes
|
||
)
|
||
|
||
return (
|
||
symbols
|
||
* np.exp(
|
||
1j * phase_values
|
||
)
|
||
)
|
||
|
||
|
||
def add_awgn(
|
||
symbols: np.ndarray,
|
||
eb_n0_db: float,
|
||
random_generator: np.random.Generator,
|
||
) -> np.ndarray:
|
||
"""
|
||
Добавить комплексный AWGN.
|
||
|
||
Энергия одного BPSK-символа равна единице.
|
||
Один символ переносит один бит.
|
||
"""
|
||
|
||
eb_n0_linear = 10.0 ** (
|
||
eb_n0_db / 10.0
|
||
)
|
||
|
||
component_sigma = np.sqrt(
|
||
1.0
|
||
/ (
|
||
2.0
|
||
* eb_n0_linear
|
||
)
|
||
)
|
||
|
||
noise = component_sigma * (
|
||
random_generator.standard_normal(
|
||
len(symbols)
|
||
)
|
||
+ 1j
|
||
* random_generator.standard_normal(
|
||
len(symbols)
|
||
)
|
||
)
|
||
|
||
return symbols + noise
|
||
|
||
|
||
# ============================================================
|
||
# Оценка фазы и CFO по известному маркеру
|
||
# ============================================================
|
||
|
||
def estimate_carrier_parameters(
|
||
received_symbols: np.ndarray,
|
||
marker_symbols: np.ndarray,
|
||
) -> dict:
|
||
"""
|
||
Оценить фазу и частотное рассогласование по маркеру.
|
||
|
||
Алгоритм:
|
||
|
||
1. Удалить известные BPSK-знаки маркера.
|
||
2. Получить грубую оценку CFO по соседним символам.
|
||
3. Выполнить уточняющий частотный поиск.
|
||
4. Оценить начальную фазу после компенсации CFO.
|
||
5. Рассчитать достоверность оценки.
|
||
"""
|
||
|
||
marker_length = len(
|
||
marker_symbols
|
||
)
|
||
|
||
received_marker = received_symbols[
|
||
:marker_length
|
||
]
|
||
|
||
if len(received_marker) != marker_length:
|
||
raise ValueError(
|
||
"Недостаточно символов маркера"
|
||
)
|
||
|
||
# Известные BPSK-знаки равны -1 или +1.
|
||
# Умножение удаляет переданную манипуляцию,
|
||
# оставляя фазу канала и шум.
|
||
despread_marker = (
|
||
received_marker
|
||
* marker_symbols
|
||
)
|
||
|
||
marker_indexes = np.arange(
|
||
marker_length,
|
||
dtype=np.float64,
|
||
)
|
||
|
||
marker_magnitude_sum = float(
|
||
np.sum(
|
||
np.abs(
|
||
despread_marker
|
||
)
|
||
)
|
||
)
|
||
|
||
# ========================================================
|
||
# Постоянная фаза без CFO-компенсации
|
||
# ========================================================
|
||
|
||
constant_coherent_sum = np.sum(
|
||
despread_marker
|
||
)
|
||
|
||
constant_phase = float(
|
||
np.angle(
|
||
constant_coherent_sum
|
||
)
|
||
)
|
||
|
||
constant_coherence = float(
|
||
np.abs(
|
||
constant_coherent_sum
|
||
)
|
||
/ (
|
||
marker_magnitude_sum
|
||
+ 1e-12
|
||
)
|
||
)
|
||
|
||
# ========================================================
|
||
# Грубая оценка CFO
|
||
# ========================================================
|
||
|
||
adjacent_products = (
|
||
despread_marker[1:]
|
||
* np.conj(
|
||
despread_marker[:-1]
|
||
)
|
||
)
|
||
|
||
adjacent_sum = np.sum(
|
||
adjacent_products
|
||
)
|
||
|
||
coarse_phase_increment = float(
|
||
np.angle(
|
||
adjacent_sum
|
||
)
|
||
)
|
||
|
||
coarse_frequency_hz = (
|
||
coarse_phase_increment
|
||
* SYMBOL_RATE
|
||
/ (2.0 * np.pi)
|
||
)
|
||
|
||
phase_consistency = float(
|
||
np.abs(
|
||
adjacent_sum
|
||
)
|
||
/ (
|
||
np.sum(
|
||
np.abs(
|
||
adjacent_products
|
||
)
|
||
)
|
||
+ 1e-12
|
||
)
|
||
)
|
||
|
||
# ========================================================
|
||
# Уточняющий поиск CFO
|
||
# ========================================================
|
||
|
||
frequency_candidates_hz = np.arange(
|
||
(
|
||
coarse_frequency_hz
|
||
- CFO_REFINEMENT_HALF_WIDTH_HZ
|
||
),
|
||
(
|
||
coarse_frequency_hz
|
||
+ CFO_REFINEMENT_HALF_WIDTH_HZ
|
||
+ CFO_REFINEMENT_STEP_HZ / 2.0
|
||
),
|
||
CFO_REFINEMENT_STEP_HZ,
|
||
dtype=np.float64,
|
||
)
|
||
|
||
phase_increment_candidates = (
|
||
2.0
|
||
* np.pi
|
||
* frequency_candidates_hz
|
||
/ SYMBOL_RATE
|
||
)
|
||
|
||
candidate_compensation = np.exp(
|
||
-1j
|
||
* phase_increment_candidates[
|
||
:, np.newaxis
|
||
]
|
||
* marker_indexes[
|
||
np.newaxis, :
|
||
]
|
||
)
|
||
|
||
coherent_sums = np.sum(
|
||
despread_marker[
|
||
np.newaxis, :
|
||
]
|
||
* candidate_compensation,
|
||
axis=1,
|
||
)
|
||
|
||
best_candidate_index = int(
|
||
np.argmax(
|
||
np.abs(
|
||
coherent_sums
|
||
)
|
||
)
|
||
)
|
||
|
||
estimated_frequency_hz = float(
|
||
frequency_candidates_hz[
|
||
best_candidate_index
|
||
]
|
||
)
|
||
|
||
phase_increment = float(
|
||
phase_increment_candidates[
|
||
best_candidate_index
|
||
]
|
||
)
|
||
|
||
best_coherent_sum = (
|
||
coherent_sums[
|
||
best_candidate_index
|
||
]
|
||
)
|
||
|
||
initial_phase_after_cfo = float(
|
||
np.angle(
|
||
best_coherent_sum
|
||
)
|
||
)
|
||
|
||
cfo_coherence = float(
|
||
np.abs(
|
||
best_coherent_sum
|
||
)
|
||
/ (
|
||
marker_magnitude_sum
|
||
+ 1e-12
|
||
)
|
||
)
|
||
|
||
coherence_gain = (
|
||
cfo_coherence
|
||
- constant_coherence
|
||
)
|
||
|
||
return {
|
||
"constant_phase": constant_phase,
|
||
"constant_coherence": constant_coherence,
|
||
"phase_increment": phase_increment,
|
||
"estimated_frequency_hz": (
|
||
estimated_frequency_hz
|
||
),
|
||
"initial_phase_after_cfo": (
|
||
initial_phase_after_cfo
|
||
),
|
||
"phase_consistency": phase_consistency,
|
||
"cfo_coherence": cfo_coherence,
|
||
"coherence_gain": coherence_gain,
|
||
}
|
||
|
||
# ============================================================
|
||
# Коррекция несущей
|
||
# ============================================================
|
||
|
||
def correct_constant_phase(
|
||
received_symbols: np.ndarray,
|
||
phase_radians: float,
|
||
) -> np.ndarray:
|
||
"""
|
||
Компенсировать только постоянную фазу.
|
||
"""
|
||
|
||
return (
|
||
received_symbols
|
||
* np.exp(
|
||
-1j * phase_radians
|
||
)
|
||
)
|
||
|
||
|
||
def correct_phase_and_frequency(
|
||
received_symbols: np.ndarray,
|
||
initial_phase_radians: float,
|
||
phase_increment: float,
|
||
) -> np.ndarray:
|
||
"""
|
||
Компенсировать постоянную фазу и CFO.
|
||
"""
|
||
|
||
symbol_indexes = np.arange(
|
||
len(received_symbols),
|
||
dtype=np.float64,
|
||
)
|
||
|
||
phase_model = (
|
||
initial_phase_radians
|
||
+ phase_increment
|
||
* symbol_indexes
|
||
)
|
||
|
||
return (
|
||
received_symbols
|
||
* np.exp(
|
||
-1j * phase_model
|
||
)
|
||
)
|
||
|
||
|
||
def should_apply_cfo_correction(
|
||
estimated_frequency_hz: float,
|
||
phase_consistency: float,
|
||
coherence_gain: float,
|
||
) -> bool:
|
||
"""
|
||
Разрешить CFO-коррекцию только при наличии
|
||
достаточных оснований.
|
||
|
||
Требования:
|
||
|
||
1. Оценка находится вне мёртвой зоны.
|
||
2. Межсимвольное вращение достаточно согласованно.
|
||
3. Компенсация CFO действительно повышает
|
||
когерентность известного маркера.
|
||
"""
|
||
|
||
MINIMUM_COHERENCE_GAIN = 0.02
|
||
|
||
return (
|
||
abs(
|
||
estimated_frequency_hz
|
||
)
|
||
>= CFO_DEAD_ZONE_HZ
|
||
and phase_consistency
|
||
>= MINIMUM_PHASE_CONSISTENCY
|
||
and coherence_gain
|
||
>= MINIMUM_COHERENCE_GAIN
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Разбор радиокадра
|
||
# ============================================================
|
||
|
||
def decode_radio_frame(
|
||
corrected_symbols: np.ndarray,
|
||
) -> str:
|
||
"""
|
||
Демодулировать полный радиокадр.
|
||
|
||
Возвращает статус приёма.
|
||
"""
|
||
|
||
received_bits = bpsk_demodulate(
|
||
corrected_symbols
|
||
)
|
||
|
||
radio_header_start = (
|
||
PREAMBLE_BIT_COUNT
|
||
)
|
||
|
||
radio_header_end = (
|
||
radio_header_start
|
||
+ 32
|
||
)
|
||
|
||
if len(received_bits) < radio_header_end:
|
||
return STATUS_HEADER_ERROR
|
||
|
||
try:
|
||
radio_header = bits_to_bytes(
|
||
received_bits[
|
||
radio_header_start:
|
||
radio_header_end
|
||
]
|
||
)
|
||
|
||
(
|
||
received_sync,
|
||
protocol_packet_length,
|
||
) = struct.unpack(
|
||
">HH",
|
||
radio_header,
|
||
)
|
||
|
||
except (ValueError, struct.error):
|
||
return STATUS_HEADER_ERROR
|
||
|
||
if received_sync != RADIO_SYNC_WORD:
|
||
return STATUS_HEADER_ERROR
|
||
|
||
if not 1 <= protocol_packet_length <= 4096:
|
||
return STATUS_HEADER_ERROR
|
||
|
||
protocol_packet_start = (
|
||
radio_header_end
|
||
)
|
||
|
||
protocol_packet_end = (
|
||
protocol_packet_start
|
||
+ protocol_packet_length * 8
|
||
)
|
||
|
||
if len(received_bits) < protocol_packet_end:
|
||
return STATUS_HEADER_ERROR
|
||
|
||
try:
|
||
protocol_packet = bits_to_bytes(
|
||
received_bits[
|
||
protocol_packet_start:
|
||
protocol_packet_end
|
||
]
|
||
)
|
||
|
||
parsed_packet = parse_packet(
|
||
protocol_packet
|
||
)
|
||
|
||
except CRCError:
|
||
return STATUS_CRC_ERROR
|
||
|
||
except (PacketError, ValueError):
|
||
return STATUS_PACKET_ERROR
|
||
|
||
try:
|
||
restored_message = (
|
||
parsed_packet.payload.decode(
|
||
"utf-8"
|
||
)
|
||
)
|
||
|
||
except UnicodeDecodeError:
|
||
return STATUS_PACKET_ERROR
|
||
|
||
if (
|
||
parsed_packet.message_type
|
||
!= MESSAGE_TYPE_TEXT
|
||
or parsed_packet.sequence_number
|
||
!= EXPECTED_SEQUENCE_NUMBER
|
||
or restored_message
|
||
!= EXPECTED_MESSAGE
|
||
):
|
||
return STATUS_PACKET_ERROR
|
||
|
||
return STATUS_SUCCESS
|
||
|
||
|
||
# ============================================================
|
||
# Подготовка исходного кадра
|
||
# ============================================================
|
||
|
||
(
|
||
transmitted_frame_bits,
|
||
transmitted_symbols,
|
||
marker_symbols,
|
||
) = build_radio_frame()
|
||
|
||
print(
|
||
"=== Lab023. Защищённая CFO-компенсация ==="
|
||
)
|
||
|
||
print("\nРазмер полного радиокадра:")
|
||
|
||
print(
|
||
len(transmitted_frame_bits),
|
||
"битов"
|
||
)
|
||
|
||
print("\nМёртвая зона CFO:")
|
||
|
||
print(
|
||
f"±{CFO_DEAD_ZONE_HZ:.1f}",
|
||
"Гц"
|
||
)
|
||
|
||
print("\nМинимальная достоверность фазовой оценки:")
|
||
|
||
print(
|
||
f"{MINIMUM_PHASE_CONSISTENCY:.2f}"
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Основной эксперимент
|
||
# ============================================================
|
||
|
||
results = []
|
||
|
||
for eb_n0_db in EB_N0_VALUES_DB:
|
||
|
||
for max_abs_cfo_hz in (
|
||
MAX_ABS_CFO_VALUES_HZ
|
||
):
|
||
|
||
seed = (
|
||
RANDOM_SEED
|
||
+ int(
|
||
eb_n0_db * 1000
|
||
)
|
||
+ int(
|
||
max_abs_cfo_hz
|
||
)
|
||
)
|
||
|
||
random_generator = np.random.default_rng(
|
||
seed
|
||
)
|
||
|
||
success_phase_only = 0
|
||
success_always_cfo = 0
|
||
success_guarded_cfo = 0
|
||
|
||
guarded_apply_count = 0
|
||
|
||
frequency_errors = []
|
||
phase_consistency_values = []
|
||
|
||
for trial_index in range(
|
||
TRIALS_PER_POINT
|
||
):
|
||
|
||
if max_abs_cfo_hz == 0.0:
|
||
true_frequency_hz = 0.0
|
||
|
||
else:
|
||
true_frequency_hz = float(
|
||
random_generator.uniform(
|
||
-max_abs_cfo_hz,
|
||
max_abs_cfo_hz,
|
||
)
|
||
)
|
||
|
||
initial_phase = float(
|
||
random_generator.uniform(
|
||
-np.pi,
|
||
np.pi,
|
||
)
|
||
)
|
||
|
||
impaired_symbols = (
|
||
apply_carrier_impairments(
|
||
symbols=transmitted_symbols,
|
||
frequency_offset_hz=(
|
||
true_frequency_hz
|
||
),
|
||
initial_phase_radians=(
|
||
initial_phase
|
||
),
|
||
)
|
||
)
|
||
|
||
received_symbols = add_awgn(
|
||
symbols=impaired_symbols,
|
||
eb_n0_db=eb_n0_db,
|
||
random_generator=random_generator,
|
||
)
|
||
|
||
estimate = estimate_carrier_parameters(
|
||
received_symbols=received_symbols,
|
||
marker_symbols=marker_symbols,
|
||
)
|
||
|
||
estimated_frequency_hz = (
|
||
estimate[
|
||
"estimated_frequency_hz"
|
||
]
|
||
)
|
||
|
||
phase_consistency = (
|
||
estimate[
|
||
"phase_consistency"
|
||
]
|
||
)
|
||
|
||
frequency_errors.append(
|
||
estimated_frequency_hz
|
||
- true_frequency_hz
|
||
)
|
||
|
||
phase_consistency_values.append(
|
||
phase_consistency
|
||
)
|
||
|
||
# =================================================
|
||
# Стратегия 1. Только постоянная фаза
|
||
# =================================================
|
||
|
||
phase_only_symbols = (
|
||
correct_constant_phase(
|
||
received_symbols=(
|
||
received_symbols
|
||
),
|
||
phase_radians=(
|
||
estimate[
|
||
"constant_phase"
|
||
]
|
||
),
|
||
)
|
||
)
|
||
|
||
status_phase_only = (
|
||
decode_radio_frame(
|
||
phase_only_symbols
|
||
)
|
||
)
|
||
|
||
if (
|
||
status_phase_only
|
||
== STATUS_SUCCESS
|
||
):
|
||
success_phase_only += 1
|
||
|
||
# =================================================
|
||
# Стратегия 2. CFO применяется всегда
|
||
# =================================================
|
||
|
||
always_cfo_symbols = (
|
||
correct_phase_and_frequency(
|
||
received_symbols=(
|
||
received_symbols
|
||
),
|
||
initial_phase_radians=(
|
||
estimate[
|
||
"initial_phase_after_cfo"
|
||
]
|
||
),
|
||
phase_increment=(
|
||
estimate[
|
||
"phase_increment"
|
||
]
|
||
),
|
||
)
|
||
)
|
||
|
||
status_always_cfo = (
|
||
decode_radio_frame(
|
||
always_cfo_symbols
|
||
)
|
||
)
|
||
|
||
if (
|
||
status_always_cfo
|
||
== STATUS_SUCCESS
|
||
):
|
||
success_always_cfo += 1
|
||
|
||
# =================================================
|
||
# Стратегия 3. Защищённая CFO-коррекция
|
||
# =================================================
|
||
|
||
apply_cfo = (
|
||
should_apply_cfo_correction(
|
||
estimated_frequency_hz=(
|
||
estimated_frequency_hz
|
||
),
|
||
phase_consistency=(
|
||
phase_consistency
|
||
),
|
||
coherence_gain=(
|
||
estimate[
|
||
"coherence_gain"
|
||
]
|
||
),
|
||
)
|
||
)
|
||
|
||
if apply_cfo:
|
||
guarded_apply_count += 1
|
||
|
||
guarded_symbols = (
|
||
always_cfo_symbols
|
||
)
|
||
|
||
else:
|
||
guarded_symbols = (
|
||
phase_only_symbols
|
||
)
|
||
|
||
status_guarded_cfo = (
|
||
decode_radio_frame(
|
||
guarded_symbols
|
||
)
|
||
)
|
||
|
||
if (
|
||
status_guarded_cfo
|
||
== STATUS_SUCCESS
|
||
):
|
||
success_guarded_cfo += 1
|
||
|
||
# =====================================================
|
||
# Итог одной экспериментальной точки
|
||
# =====================================================
|
||
|
||
frequency_errors_array = np.asarray(
|
||
frequency_errors,
|
||
dtype=np.float64,
|
||
)
|
||
|
||
frequency_rmse_hz = float(
|
||
np.sqrt(
|
||
np.mean(
|
||
frequency_errors_array ** 2
|
||
)
|
||
)
|
||
)
|
||
|
||
mean_phase_consistency = float(
|
||
np.mean(
|
||
phase_consistency_values
|
||
)
|
||
)
|
||
|
||
fer_phase_only = (
|
||
1.0
|
||
- success_phase_only
|
||
/ TRIALS_PER_POINT
|
||
)
|
||
|
||
fer_always_cfo = (
|
||
1.0
|
||
- success_always_cfo
|
||
/ TRIALS_PER_POINT
|
||
)
|
||
|
||
fer_guarded_cfo = (
|
||
1.0
|
||
- success_guarded_cfo
|
||
/ TRIALS_PER_POINT
|
||
)
|
||
|
||
guarded_apply_rate = (
|
||
guarded_apply_count
|
||
/ TRIALS_PER_POINT
|
||
)
|
||
|
||
results.append(
|
||
{
|
||
"eb_n0_db": eb_n0_db,
|
||
"max_abs_cfo_hz": (
|
||
max_abs_cfo_hz
|
||
),
|
||
"trials": TRIALS_PER_POINT,
|
||
"success_phase_only": (
|
||
success_phase_only
|
||
),
|
||
"success_always_cfo": (
|
||
success_always_cfo
|
||
),
|
||
"success_guarded_cfo": (
|
||
success_guarded_cfo
|
||
),
|
||
"fer_phase_only": (
|
||
fer_phase_only
|
||
),
|
||
"fer_always_cfo": (
|
||
fer_always_cfo
|
||
),
|
||
"fer_guarded_cfo": (
|
||
fer_guarded_cfo
|
||
),
|
||
"guarded_apply_count": (
|
||
guarded_apply_count
|
||
),
|
||
"guarded_apply_rate": (
|
||
guarded_apply_rate
|
||
),
|
||
"frequency_rmse_hz": (
|
||
frequency_rmse_hz
|
||
),
|
||
"mean_phase_consistency": (
|
||
mean_phase_consistency
|
||
),
|
||
}
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Вывод результатов
|
||
# ============================================================
|
||
|
||
print("\nПопыток на каждую точку:")
|
||
|
||
print(TRIALS_PER_POINT)
|
||
|
||
print("\nРезультаты:")
|
||
|
||
print(
|
||
f"{'Eb/N0':>9}"
|
||
f"{'CFO ±':>10}"
|
||
f"{'Усп. фаза':>12}"
|
||
f"{'Усп. всегда':>14}"
|
||
f"{'Усп. guard':>13}"
|
||
f"{'FER фаза':>11}"
|
||
f"{'FER всегда':>12}"
|
||
f"{'FER guard':>11}"
|
||
f"{'Guard on':>10}"
|
||
f"{'CFO RMSE':>11}"
|
||
)
|
||
|
||
print("-" * 114)
|
||
|
||
for result in results:
|
||
|
||
print(
|
||
f"{result['eb_n0_db']:>6.1f} дБ"
|
||
f"{result['max_abs_cfo_hz']:>7.0f} Гц"
|
||
f"{result['success_phase_only']:>12}"
|
||
f"{result['success_always_cfo']:>14}"
|
||
f"{result['success_guarded_cfo']:>13}"
|
||
f"{result['fer_phase_only']:>11.3f}"
|
||
f"{result['fer_always_cfo']:>12.3f}"
|
||
f"{result['fer_guarded_cfo']:>11.3f}"
|
||
f"{result['guarded_apply_rate'] * 100:>8.1f} %"
|
||
f"{result['frequency_rmse_hz']:>8.2f} Гц"
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Сохранение CSV
|
||
# ============================================================
|
||
|
||
with CSV_PATH.open(
|
||
"w",
|
||
newline="",
|
||
encoding="utf-8-sig",
|
||
) as csv_file:
|
||
|
||
writer = DictWriter(
|
||
csv_file,
|
||
fieldnames=list(
|
||
results[0].keys()
|
||
),
|
||
)
|
||
|
||
writer.writeheader()
|
||
writer.writerows(results)
|
||
|
||
|
||
# ============================================================
|
||
# Подготовка графиков
|
||
# ============================================================
|
||
|
||
measurement_floor = (
|
||
0.5
|
||
/ TRIALS_PER_POINT
|
||
)
|
||
|
||
figure, axes = plt.subplots(
|
||
2,
|
||
2,
|
||
figsize=(13, 10),
|
||
)
|
||
|
||
|
||
# ------------------------------------------------------------
|
||
# 1. FER при истинном CFO = 0
|
||
# ------------------------------------------------------------
|
||
|
||
zero_cfo_results = [
|
||
result
|
||
for result in results
|
||
if result["max_abs_cfo_hz"] == 0.0
|
||
]
|
||
|
||
zero_eb_values = [
|
||
result["eb_n0_db"]
|
||
for result in zero_cfo_results
|
||
]
|
||
|
||
axes[0, 0].semilogy(
|
||
zero_eb_values,
|
||
[
|
||
max(
|
||
result["fer_phase_only"],
|
||
measurement_floor,
|
||
)
|
||
for result in zero_cfo_results
|
||
],
|
||
marker="o",
|
||
label="Только постоянная фаза",
|
||
)
|
||
|
||
axes[0, 0].semilogy(
|
||
zero_eb_values,
|
||
[
|
||
max(
|
||
result["fer_always_cfo"],
|
||
measurement_floor,
|
||
)
|
||
for result in zero_cfo_results
|
||
],
|
||
marker="s",
|
||
label="CFO применяется всегда",
|
||
)
|
||
|
||
axes[0, 0].semilogy(
|
||
zero_eb_values,
|
||
[
|
||
max(
|
||
result["fer_guarded_cfo"],
|
||
measurement_floor,
|
||
)
|
||
for result in zero_cfo_results
|
||
],
|
||
marker="^",
|
||
label="Защищённая CFO-коррекция",
|
||
)
|
||
|
||
axes[0, 0].set_xlabel(
|
||
"Eb/N0, дБ"
|
||
)
|
||
|
||
axes[0, 0].set_ylabel(
|
||
"FER"
|
||
)
|
||
|
||
axes[0, 0].set_title(
|
||
"Истинный CFO = 0 Гц"
|
||
)
|
||
|
||
axes[0, 0].grid(
|
||
True,
|
||
which="both",
|
||
)
|
||
|
||
axes[0, 0].legend()
|
||
|
||
|
||
# ------------------------------------------------------------
|
||
# 2. FER при CFO ±1000 Гц
|
||
# ------------------------------------------------------------
|
||
|
||
large_cfo_results = [
|
||
result
|
||
for result in results
|
||
if result["max_abs_cfo_hz"] == 1000.0
|
||
]
|
||
|
||
large_eb_values = [
|
||
result["eb_n0_db"]
|
||
for result in large_cfo_results
|
||
]
|
||
|
||
axes[0, 1].semilogy(
|
||
large_eb_values,
|
||
[
|
||
max(
|
||
result["fer_phase_only"],
|
||
measurement_floor,
|
||
)
|
||
for result in large_cfo_results
|
||
],
|
||
marker="o",
|
||
label="Только постоянная фаза",
|
||
)
|
||
|
||
axes[0, 1].semilogy(
|
||
large_eb_values,
|
||
[
|
||
max(
|
||
result["fer_always_cfo"],
|
||
measurement_floor,
|
||
)
|
||
for result in large_cfo_results
|
||
],
|
||
marker="s",
|
||
label="CFO применяется всегда",
|
||
)
|
||
|
||
axes[0, 1].semilogy(
|
||
large_eb_values,
|
||
[
|
||
max(
|
||
result["fer_guarded_cfo"],
|
||
measurement_floor,
|
||
)
|
||
for result in large_cfo_results
|
||
],
|
||
marker="^",
|
||
label="Защищённая CFO-коррекция",
|
||
)
|
||
|
||
axes[0, 1].set_xlabel(
|
||
"Eb/N0, дБ"
|
||
)
|
||
|
||
axes[0, 1].set_ylabel(
|
||
"FER"
|
||
)
|
||
|
||
axes[0, 1].set_title(
|
||
"Случайный CFO в диапазоне ±1000 Гц"
|
||
)
|
||
|
||
axes[0, 1].grid(
|
||
True,
|
||
which="both",
|
||
)
|
||
|
||
axes[0, 1].legend()
|
||
|
||
|
||
# ------------------------------------------------------------
|
||
# 3. Доля включения CFO-коррекции
|
||
# ------------------------------------------------------------
|
||
|
||
for max_abs_cfo_hz in (
|
||
MAX_ABS_CFO_VALUES_HZ
|
||
):
|
||
|
||
selected_results = [
|
||
result
|
||
for result in results
|
||
if (
|
||
result["max_abs_cfo_hz"]
|
||
== max_abs_cfo_hz
|
||
)
|
||
]
|
||
|
||
axes[1, 0].plot(
|
||
[
|
||
result["eb_n0_db"]
|
||
for result in selected_results
|
||
],
|
||
[
|
||
result["guarded_apply_rate"]
|
||
for result in selected_results
|
||
],
|
||
marker="o",
|
||
label=(
|
||
f"CFO ±{max_abs_cfo_hz:.0f} Гц"
|
||
),
|
||
)
|
||
|
||
axes[1, 0].set_xlabel(
|
||
"Eb/N0, дБ"
|
||
)
|
||
|
||
axes[1, 0].set_ylabel(
|
||
"Доля включения CFO-коррекции"
|
||
)
|
||
|
||
axes[1, 0].set_ylim(
|
||
-0.03,
|
||
1.03,
|
||
)
|
||
|
||
axes[1, 0].set_title(
|
||
"Решения защищённого корректора"
|
||
)
|
||
|
||
axes[1, 0].grid(
|
||
True
|
||
)
|
||
|
||
axes[1, 0].legend()
|
||
|
||
|
||
# ------------------------------------------------------------
|
||
# 4. Ошибка оценки CFO
|
||
# ------------------------------------------------------------
|
||
|
||
for max_abs_cfo_hz in (
|
||
MAX_ABS_CFO_VALUES_HZ
|
||
):
|
||
|
||
selected_results = [
|
||
result
|
||
for result in results
|
||
if (
|
||
result["max_abs_cfo_hz"]
|
||
== max_abs_cfo_hz
|
||
)
|
||
]
|
||
|
||
axes[1, 1].plot(
|
||
[
|
||
result["eb_n0_db"]
|
||
for result in selected_results
|
||
],
|
||
[
|
||
result["frequency_rmse_hz"]
|
||
for result in selected_results
|
||
],
|
||
marker="o",
|
||
label=(
|
||
f"CFO ±{max_abs_cfo_hz:.0f} Гц"
|
||
),
|
||
)
|
||
|
||
axes[1, 1].set_xlabel(
|
||
"Eb/N0, дБ"
|
||
)
|
||
|
||
axes[1, 1].set_ylabel(
|
||
"CFO RMSE, Гц"
|
||
)
|
||
|
||
axes[1, 1].set_title(
|
||
"Точность оценки частотной ошибки"
|
||
)
|
||
|
||
axes[1, 1].grid(
|
||
True
|
||
)
|
||
|
||
axes[1, 1].legend()
|
||
|
||
|
||
figure.tight_layout()
|
||
|
||
figure.savefig(
|
||
GRAPH_PATH,
|
||
dpi=160,
|
||
)
|
||
|
||
plt.close(
|
||
figure
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Текстовый отчёт
|
||
# ============================================================
|
||
|
||
report_lines = [
|
||
"Lab023. Guarded CFO correction",
|
||
"",
|
||
(
|
||
"Trials per point: "
|
||
f"{TRIALS_PER_POINT}"
|
||
),
|
||
(
|
||
"CFO dead zone: "
|
||
f"{CFO_DEAD_ZONE_HZ:.2f} Hz"
|
||
),
|
||
(
|
||
"Minimum phase consistency: "
|
||
f"{MINIMUM_PHASE_CONSISTENCY:.3f}"
|
||
),
|
||
"",
|
||
]
|
||
|
||
for result in results:
|
||
|
||
report_lines.extend(
|
||
[
|
||
(
|
||
"Eb/N0: "
|
||
f"{result['eb_n0_db']:.1f} dB"
|
||
),
|
||
(
|
||
"Maximum absolute CFO: "
|
||
f"{result['max_abs_cfo_hz']:.1f} Hz"
|
||
),
|
||
(
|
||
"FER phase only: "
|
||
f"{result['fer_phase_only']:.6f}"
|
||
),
|
||
(
|
||
"FER always CFO: "
|
||
f"{result['fer_always_cfo']:.6f}"
|
||
),
|
||
(
|
||
"FER guarded CFO: "
|
||
f"{result['fer_guarded_cfo']:.6f}"
|
||
),
|
||
(
|
||
"Guard apply rate: "
|
||
f"{result['guarded_apply_rate']:.6f}"
|
||
),
|
||
(
|
||
"CFO RMSE: "
|
||
f"{result['frequency_rmse_hz']:.3f} Hz"
|
||
),
|
||
"",
|
||
]
|
||
)
|
||
|
||
REPORT_PATH.write_text(
|
||
"\n".join(
|
||
report_lines
|
||
),
|
||
encoding="utf-8",
|
||
)
|
||
|
||
|
||
# ============================================================
|
||
# Автоматические проверки
|
||
# ============================================================
|
||
|
||
result_10_db_zero_cfo = next(
|
||
result
|
||
for result in results
|
||
if (
|
||
result["eb_n0_db"] == 10.0
|
||
and result["max_abs_cfo_hz"] == 0.0
|
||
)
|
||
)
|
||
|
||
result_10_db_large_cfo = next(
|
||
result
|
||
for result in results
|
||
if (
|
||
result["eb_n0_db"] == 10.0
|
||
and result["max_abs_cfo_hz"] == 1000.0
|
||
)
|
||
)
|
||
|
||
assert len(transmitted_frame_bits) == 320
|
||
|
||
# При реальном CFO = 0 защищённый режим
|
||
# не должен постоянно включать частотную коррекцию.
|
||
assert (
|
||
result_10_db_zero_cfo[
|
||
"guarded_apply_rate"
|
||
]
|
||
< 0.25
|
||
)
|
||
|
||
# При большом CFO защищённый корректор должен
|
||
# быть намного лучше режима только с постоянной фазой.
|
||
assert (
|
||
result_10_db_large_cfo[
|
||
"success_guarded_cfo"
|
||
]
|
||
> result_10_db_large_cfo[
|
||
"success_phase_only"
|
||
]
|
||
)
|
||
|
||
# При хорошем SNR защищённый режим должен сохранять
|
||
# практически все кадры.
|
||
assert (
|
||
result_10_db_large_cfo[
|
||
"success_guarded_cfo"
|
||
]
|
||
>= int(
|
||
TRIALS_PER_POINT * 0.90
|
||
)
|
||
)
|
||
|
||
assert CSV_PATH.exists()
|
||
assert GRAPH_PATH.exists()
|
||
assert REPORT_PATH.exists()
|
||
|
||
|
||
print("\nCSV:")
|
||
|
||
print(CSV_PATH)
|
||
|
||
print("\nГрафик:")
|
||
|
||
print(GRAPH_PATH)
|
||
|
||
print("\nОтчёт:")
|
||
|
||
print(REPORT_PATH)
|
||
|
||
print(
|
||
"\nПроверка пройдена: "
|
||
"мёртвая зона исключает лишнюю CFO-коррекцию, "
|
||
"но сохраняет компенсацию значительного рассогласования."
|
||
) |