turov_av/SDR-Rover
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
SDR-Rover/experiments/lab022_noise_and_frequency_offset.py
LittleSam129 c486039053 Split experiments from tests 
The tests/ directory held 50 laboratory programs and no tests. They model
channels, run hundreds of repetitions and write CSV, PNG and reports;
calling that a test suite blocked introducing a real one, because any
pytest run would have collected the labs and re-executed every
experiment.

- move all 50 lab programs to experiments/ with git mv, preserving history
- rewrite the 38 cross-imports between labs from tests.labNNN to
  experiments.labNNN
- leave tests/ empty for actual fast checks of protocol/
- point quick_gate and the hook at the new layout and add experiments/ to
  the syntax sweep
- update the paths quoted in the Lab042 specification and the verifier
  agent definition

This also defuses the import-time work finding without touching 41 files:
the labs still create directories and write files on import, but nothing
imports them now except the gate, which does so deliberately.

Gate passes: syntax clean, protocol imports, 15 lab modules import, 2
functional suites run.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
2026-08-10 14:34:58 +03:00

1716 lines
37 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
"""
Lab022. Совместное действие AWGN и частотного рассогласования.
Программа многократно передаёт BPSK-радиокадр из Lab018.
В каждой попытке случайно выбираются:
- задержка начала кадра;
- начальный фазовый поворот;
- частотное рассогласование;
- реализация AWGN-шума.
Сравниваются два варианта приёма:
1. Только компенсация постоянной фазы.
2. Компенсация постоянной фазы и CFO
(Carrier Frequency Offset).
Для каждой комбинации Eb/N0 и максимального CFO
рассчитываются:
- FER без компенсации CFO;
- FER после компенсации CFO;
- число необнаруженных кадров;
- средняя корреляционная оценка;
- средняя и среднеквадратичная ошибка оценки CFO.
"""
from csv import DictWriter
from pathlib import Path
import struct
import matplotlib.pyplot as plt
import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
from scipy.signal import fftconvolve
from protocol.packet import (
CRCError,
MESSAGE_TYPE_TEXT,
PacketError,
parse_packet,
)
# ============================================================
# Параметры радиокадра
# ============================================================
RADIO_SYNC_WORD = 0xD391
PREAMBLE_BIT_COUNT = 64
RADIO_FRAME_BIT_COUNT = 320
SAMPLES_PER_SYMBOL = 32
SYMBOL_RATE = 20_000
SAMPLE_RATE = (
SYMBOL_RATE
* SAMPLES_PER_SYMBOL
)
RRC_ROLLOFF = 0.35
RRC_SPAN_SYMBOLS = 10
# ============================================================
# Контрольные данные внутреннего пакета
# ============================================================
EXPECTED_MESSAGE = "ПРИВЕТ SDR"
EXPECTED_SEQUENCE_NUMBER = 18
# ============================================================
# Настройки эксперимента
# ============================================================
EB_N0_VALUES_DB = [
6.0,
7.0,
8.0,
10.0,
]
# В каждой попытке реальный CFO выбирается случайно
# из диапазона:
#
# -MAX_ABS_CFO ... +MAX_ABS_CFO
MAX_ABS_CFO_VALUES_HZ = [
0.0,
250.0,
1000.0,
]
TRIALS_PER_POINT = 30
RANDOM_SEED = 2026
MAX_RANDOM_DELAY_SAMPLES = (
2
* SAMPLES_PER_SYMBOL
)
# При меньшей оценке маркер считаем необнаруженным.
DETECTION_THRESHOLD = 0.45
# После грубой оценки CFO выполняется уточняющий поиск.
CFO_REFINEMENT_HALF_WIDTH_HZ = 300.0
CFO_REFINEMENT_STEP_HZ = 1.0
MAX_PROTOCOL_PACKET_SIZE = 4096
# ============================================================
# Пути
# ============================================================
INPUT_IQ_PATH = Path(
"data/processed/lab018/"
"lab018_bpsk_tx_iq.npy"
)
OUTPUT_DIRECTORY = Path(
"data/processed/lab022"
)
OUTPUT_DIRECTORY.mkdir(
parents=True,
exist_ok=True,
)
CSV_PATH = (
OUTPUT_DIRECTORY
/ "lab022_noise_cfo_results.csv"
)
GRAPH_PATH = (
OUTPUT_DIRECTORY
/ "lab022_noise_cfo.png"
)
REPORT_PATH = (
OUTPUT_DIRECTORY
/ "lab022_noise_cfo_report.txt"
)
# ============================================================
# Статусы приёма
# ============================================================
STATUS_SUCCESS = "SUCCESS"
STATUS_DETECTION_FAIL = "DETECTION FAIL"
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)
# ============================================================
# Root Raised Cosine
# ============================================================
def root_raised_cosine_taps(
rolloff: float,
samples_per_symbol: int,
span_symbols: int,
) -> np.ndarray:
"""
Рассчитать коэффициенты RRC-фильтра.
"""
if not 0.0 < rolloff <= 1.0:
raise ValueError(
"rolloff должен быть в диапазоне 0...1"
)
if samples_per_symbol <= 0:
raise ValueError(
"samples_per_symbol должен быть положительным"
)
if span_symbols <= 0:
raise ValueError(
"span_symbols должен быть положительным"
)
if span_symbols % 2 != 0:
raise ValueError(
"span_symbols должен быть чётным"
)
half_sample_count = (
span_symbols
* samples_per_symbol
// 2
)
sample_indexes = np.arange(
-half_sample_count,
half_sample_count + 1,
dtype=np.float64,
)
time_values = (
sample_indexes
/ samples_per_symbol
)
taps = np.zeros_like(
time_values
)
beta = rolloff
for index, time_value in enumerate(
time_values
):
if np.isclose(
time_value,
0.0,
):
taps[index] = (
1.0
- beta
+ 4.0 * beta / np.pi
)
continue
if np.isclose(
abs(time_value),
1.0 / (4.0 * beta),
):
taps[index] = (
beta
/ np.sqrt(2.0)
* (
(
1.0
+ 2.0 / np.pi
)
* np.sin(
np.pi
/ (4.0 * beta)
)
+ (
1.0
- 2.0 / np.pi
)
* np.cos(
np.pi
/ (4.0 * beta)
)
)
)
continue
numerator = (
np.sin(
np.pi
* time_value
* (1.0 - beta)
)
+ (
4.0
* beta
* time_value
* np.cos(
np.pi
* time_value
* (1.0 + beta)
)
)
)
denominator = (
np.pi
* time_value
* (
1.0
- (
4.0
* beta
* time_value
) ** 2
)
)
taps[index] = (
numerator
/ denominator
)
taps /= np.sqrt(
np.sum(
taps ** 2
)
)
return taps
# ============================================================
# Маркер PREAMBLE + RADIO SYNC
# ============================================================
def build_frame_marker(
) -> tuple[np.ndarray, np.ndarray]:
"""
Сформировать известный маркер радиокадра.
"""
preamble_bits = np.tile(
np.array(
[1, 0],
dtype=np.uint8,
),
PREAMBLE_BIT_COUNT // 2,
)
sync_bits = bytes_to_bits(
struct.pack(
">H",
RADIO_SYNC_WORD,
)
)
marker_bits = np.concatenate(
[
preamble_bits,
sync_bits,
]
)
marker_symbols = bpsk_modulate(
marker_bits
)
return marker_bits, marker_symbols
# ============================================================
# Искажения канала
# ============================================================
def apply_frequency_offset(
iq_samples: np.ndarray,
frequency_offset_hz: float,
phase_offset_radians: float,
) -> np.ndarray:
"""
Добавить начальную фазу и частотное рассогласование.
"""
sample_indexes = np.arange(
len(iq_samples),
dtype=np.float64,
)
phase_values = (
phase_offset_radians
+ 2.0
* np.pi
* frequency_offset_hz
* sample_indexes
/ SAMPLE_RATE
)
return (
iq_samples
* np.exp(
1j * phase_values
)
)
def add_awgn_for_eb_n0(
clean_iq: np.ndarray,
signal_energy_per_bit: float,
eb_n0_db: float,
random_generator: np.random.Generator,
) -> np.ndarray:
"""
Добавить комплексный AWGN для заданного Eb/N0.
"""
eb_n0_linear = 10.0 ** (
eb_n0_db / 10.0
)
complex_noise_variance = (
signal_energy_per_bit
/ eb_n0_linear
)
component_sigma = np.sqrt(
complex_noise_variance / 2.0
)
noise = component_sigma * (
random_generator.standard_normal(
len(clean_iq)
)
+ 1j
* random_generator.standard_normal(
len(clean_iq)
)
)
return clean_iq + noise
# ============================================================
# Поиск кадра и оценка CFO
# ============================================================
def find_frame_and_frequency_offset(
matched_iq: np.ndarray,
marker_symbols: np.ndarray,
) -> dict:
"""
Найти кадр и оценить частотное рассогласование.
Алгоритм состоит из двух этапов.
Этап 1:
грубая оценка CFO по межсимвольному
изменению фазы известного маркера.
Этап 2:
уточняющий частотный поиск около
полученной грубой оценки.
"""
marker_length = len(
marker_symbols
)
marker_indexes = np.arange(
marker_length,
dtype=np.float64,
)
marker_energy = float(
np.sum(
np.abs(marker_symbols) ** 2
)
)
best_result = None
# --------------------------------------------------------
# Грубый поиск по всем фазам дискретизации
# --------------------------------------------------------
for sample_phase in range(
SAMPLES_PER_SYMBOL
):
symbol_samples = matched_iq[
sample_phase::SAMPLES_PER_SYMBOL
]
if len(symbol_samples) < marker_length:
continue
windows = sliding_window_view(
symbol_samples,
marker_length,
)
# Удаляем известные знаки BPSK-маркера.
despread_windows = (
windows
* marker_symbols[
np.newaxis,
:
]
)
adjacent_products = (
despread_windows[:, 1:]
* np.conj(
despread_windows[:, :-1]
)
)
coarse_phase_increments = np.angle(
np.sum(
adjacent_products,
axis=1,
)
)
coarse_compensation = np.exp(
-1j
* coarse_phase_increments[
:, np.newaxis
]
* marker_indexes[
np.newaxis, :
]
)
coherent_sums = np.sum(
despread_windows
* coarse_compensation,
axis=1,
)
window_energy = np.sum(
np.abs(windows) ** 2,
axis=1,
)
scores = (
np.abs(coherent_sums)
/ (
np.sqrt(
window_energy
* marker_energy
)
+ 1e-12
)
)
start_index = int(
np.argmax(
scores
)
)
score = float(
scores[start_index]
)
if (
best_result is None
or score > best_result["score"]
):
best_result = {
"score": score,
"sample_phase": sample_phase,
"start_symbol_index": start_index,
"symbol_samples": symbol_samples,
"despread_marker": (
despread_windows[
start_index
].copy()
),
"coarse_phase_increment": float(
coarse_phase_increments[
start_index
]
),
}
if best_result is None:
raise RuntimeError(
"Не удалось выполнить поиск радиокадра"
)
# --------------------------------------------------------
# Уточнение CFO частотным поиском
# --------------------------------------------------------
coarse_frequency_hz = (
best_result[
"coarse_phase_increment"
]
* SYMBOL_RATE
/ (2.0 * np.pi)
)
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, :
]
)
refined_coherent_sums = np.sum(
best_result[
"despread_marker"
][
np.newaxis, :
]
* candidate_compensation,
axis=1,
)
best_frequency_index = int(
np.argmax(
np.abs(
refined_coherent_sums
)
)
)
estimated_frequency_hz = float(
frequency_candidates_hz[
best_frequency_index
]
)
estimated_phase_increment = float(
phase_increment_candidates[
best_frequency_index
]
)
estimated_initial_phase = float(
np.angle(
refined_coherent_sums[
best_frequency_index
]
)
)
selected_marker_energy = float(
np.sum(
np.abs(
best_result[
"despread_marker"
]
) ** 2
)
)
refined_score = float(
np.abs(
refined_coherent_sums[
best_frequency_index
]
)
/ (
np.sqrt(
selected_marker_energy
* marker_length
)
+ 1e-12
)
)
best_result.update(
{
"score": refined_score,
"frequency_offset_hz": (
estimated_frequency_hz
),
"phase_increment": (
estimated_phase_increment
),
"initial_phase": (
estimated_initial_phase
),
}
)
return best_result
# ============================================================
# Разбор восстановленного радиокадра
# ============================================================
def decode_radio_frame(
corrected_symbols: np.ndarray,
frame_start_symbol: int,
) -> str:
"""
Демодулировать радиокадр и вернуть статус.
"""
received_bits = bpsk_demodulate(
corrected_symbols
)
available_bits = received_bits[
frame_start_symbol:
]
radio_header_start = (
PREAMBLE_BIT_COUNT
)
radio_header_end = (
radio_header_start
+ 32
)
if len(available_bits) < radio_header_end:
return STATUS_HEADER_ERROR
try:
radio_header = bits_to_bytes(
available_bits[
radio_header_start:
radio_header_end
]
)
(
received_radio_sync,
protocol_packet_length,
) = struct.unpack(
">HH",
radio_header,
)
except (ValueError, struct.error):
return STATUS_HEADER_ERROR
if received_radio_sync != RADIO_SYNC_WORD:
return STATUS_HEADER_ERROR
if not (
1
<= protocol_packet_length
<= MAX_PROTOCOL_PACKET_SIZE
):
return STATUS_HEADER_ERROR
protocol_packet_start = (
radio_header_end
)
protocol_packet_end = (
protocol_packet_start
+ protocol_packet_length * 8
)
if len(available_bits) < protocol_packet_end:
return STATUS_HEADER_ERROR
try:
protocol_packet = bits_to_bytes(
available_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
# ============================================================
# Загрузка IQ-сигнала
# ============================================================
if not INPUT_IQ_PATH.exists():
raise FileNotFoundError(
f"Не найден IQ-файл: {INPUT_IQ_PATH}. "
"Сначала необходимо выполнить Lab018."
)
transmitted_iq = np.load(
INPUT_IQ_PATH
).astype(
np.complex128
)
if transmitted_iq.ndim != 1:
raise ValueError(
"IQ-массив должен быть одномерным"
)
if not np.iscomplexobj(
transmitted_iq
):
raise ValueError(
"IQ-массив должен быть комплексным"
)
# ============================================================
# Подготовка приёмника
# ============================================================
rrc_taps = root_raised_cosine_taps(
rolloff=RRC_ROLLOFF,
samples_per_symbol=SAMPLES_PER_SYMBOL,
span_symbols=RRC_SPAN_SYMBOLS,
)
marker_bits, marker_symbols = (
build_frame_marker()
)
signal_energy_per_bit = (
np.sum(
np.abs(transmitted_iq) ** 2
)
/ RADIO_FRAME_BIT_COUNT
)
# ============================================================
# Основной эксперимент
# ============================================================
results = []
for eb_n0_db in EB_N0_VALUES_DB:
for max_abs_cfo_hz in (
MAX_ABS_CFO_VALUES_HZ
):
random_generator = (
np.random.default_rng(
RANDOM_SEED
+ int(
eb_n0_db * 100
)
+ int(
max_abs_cfo_hz
)
)
)
success_without_cfo = 0
success_with_cfo = 0
detection_fail_count = 0
correlation_scores = []
frequency_errors_hz = []
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,
)
)
phase_offset_radians = float(
random_generator.uniform(
-np.pi,
np.pi,
)
)
sample_delay = int(
random_generator.integers(
0,
MAX_RANDOM_DELAY_SAMPLES + 1,
)
)
impaired_iq = apply_frequency_offset(
iq_samples=transmitted_iq,
frequency_offset_hz=(
true_frequency_hz
),
phase_offset_radians=(
phase_offset_radians
),
)
clean_received_iq = np.concatenate(
[
np.zeros(
sample_delay,
dtype=np.complex128,
),
impaired_iq,
np.zeros(
2 * SAMPLES_PER_SYMBOL,
dtype=np.complex128,
),
]
)
noisy_received_iq = add_awgn_for_eb_n0(
clean_iq=clean_received_iq,
signal_energy_per_bit=(
signal_energy_per_bit
),
eb_n0_db=eb_n0_db,
random_generator=(
random_generator
),
)
matched_iq = fftconvolve(
noisy_received_iq,
rrc_taps,
mode="full",
)
search_result = (
find_frame_and_frequency_offset(
matched_iq=matched_iq,
marker_symbols=marker_symbols,
)
)
correlation_score = (
search_result["score"]
)
correlation_scores.append(
correlation_score
)
if (
correlation_score
< DETECTION_THRESHOLD
):
detection_fail_count += 1
continue
symbol_samples = search_result[
"symbol_samples"
]
frame_start_symbol = search_result[
"start_symbol_index"
]
estimated_frequency_hz = (
search_result[
"frequency_offset_hz"
]
)
frequency_errors_hz.append(
estimated_frequency_hz
- true_frequency_hz
)
estimated_initial_phase = (
search_result[
"initial_phase"
]
)
estimated_phase_increment = (
search_result[
"phase_increment"
]
)
symbol_indexes = np.arange(
len(symbol_samples),
dtype=np.float64,
)
relative_symbol_indexes = (
symbol_indexes
- frame_start_symbol
)
# -----------------------------------------------
# Только постоянная фазовая коррекция
# -----------------------------------------------
constant_phase_corrected = (
symbol_samples
* np.exp(
-1j
* estimated_initial_phase
)
)
status_without = decode_radio_frame(
corrected_symbols=(
constant_phase_corrected
),
frame_start_symbol=(
frame_start_symbol
),
)
if status_without == STATUS_SUCCESS:
success_without_cfo += 1
# -----------------------------------------------
# Фазовая и частотная коррекция
# -----------------------------------------------
complete_phase_model = (
estimated_initial_phase
+ estimated_phase_increment
* relative_symbol_indexes
)
frequency_corrected = (
symbol_samples
* np.exp(
-1j
* complete_phase_model
)
)
status_with = decode_radio_frame(
corrected_symbols=(
frequency_corrected
),
frame_start_symbol=(
frame_start_symbol
),
)
if status_with == STATUS_SUCCESS:
success_with_cfo += 1
fer_without_cfo = (
1.0
- success_without_cfo
/ TRIALS_PER_POINT
)
fer_with_cfo = (
1.0
- success_with_cfo
/ TRIALS_PER_POINT
)
mean_correlation = float(
np.mean(
correlation_scores
)
)
if frequency_errors_hz:
frequency_errors_array = np.array(
frequency_errors_hz,
dtype=np.float64,
)
mean_abs_frequency_error_hz = float(
np.mean(
np.abs(
frequency_errors_array
)
)
)
frequency_rmse_hz = float(
np.sqrt(
np.mean(
frequency_errors_array ** 2
)
)
)
else:
mean_abs_frequency_error_hz = float(
"nan"
)
frequency_rmse_hz = float(
"nan"
)
results.append(
{
"eb_n0_db": eb_n0_db,
"max_abs_cfo_hz": (
max_abs_cfo_hz
),
"trials": TRIALS_PER_POINT,
"success_without_cfo": (
success_without_cfo
),
"success_with_cfo": (
success_with_cfo
),
"detection_fail_count": (
detection_fail_count
),
"fer_without_cfo": (
fer_without_cfo
),
"fer_with_cfo": (
fer_with_cfo
),
"mean_correlation": (
mean_correlation
),
"mean_abs_frequency_error_hz": (
mean_abs_frequency_error_hz
),
"frequency_rmse_hz": (
frequency_rmse_hz
),
}
)
# ============================================================
# Вывод таблицы
# ============================================================
print(
"=== Lab022. Шум и частотное рассогласование ==="
)
print("\nПопыток на каждую точку:")
print(TRIALS_PER_POINT)
print("\nРезультаты:")
print(
f"{'Eb/N0':>9}"
f"{'CFO ±':>10}"
f"{'Успех без':>12}"
f"{'Успех с':>11}"
f"{'FER без':>10}"
f"{'FER с':>10}"
f"{'Нет кадра':>11}"
f"{'CFO RMSE':>12}"
f"{'Коррел.':>10}"
)
print("-" * 95)
for result in results:
print(
f"{result['eb_n0_db']:>6.1f} дБ"
f"{result['max_abs_cfo_hz']:>7.0f} Гц"
f"{result['success_without_cfo']:>12}"
f"{result['success_with_cfo']:>11}"
f"{result['fer_without_cfo']:>10.3f}"
f"{result['fer_with_cfo']:>10.3f}"
f"{result['detection_fail_count']:>11}"
f"{result['frequency_rmse_hz']:>9.2f} Гц"
f"{result['mean_correlation']:>10.3f}"
)
# ============================================================
# Сохранение 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)
# ============================================================
# Графики
# ============================================================
figure, axes = plt.subplots(
2,
2,
figsize=(13, 10),
)
# ------------------------------------------------------------
# 1. FER без CFO-компенсации
# ------------------------------------------------------------
measurement_floor = (
0.5
/ TRIALS_PER_POINT
)
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
)
]
eb_values = [
result["eb_n0_db"]
for result in selected_results
]
fer_values = [
max(
result["fer_without_cfo"],
measurement_floor,
)
for result in selected_results
]
axes[0, 0].semilogy(
eb_values,
fer_values,
marker="o",
label=(
f"CFO ±{max_abs_cfo_hz:.0f} Гц"
),
)
axes[0, 0].set_xlabel(
"Eb/N0, дБ"
)
axes[0, 0].set_ylabel(
"FER"
)
axes[0, 0].set_title(
"FER без компенсации CFO"
)
axes[0, 0].grid(
True,
which="both",
)
axes[0, 0].legend()
# ------------------------------------------------------------
# 2. FER после 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
)
]
eb_values = [
result["eb_n0_db"]
for result in selected_results
]
fer_values = [
max(
result["fer_with_cfo"],
measurement_floor,
)
for result in selected_results
]
axes[0, 1].semilogy(
eb_values,
fer_values,
marker="s",
label=(
f"CFO ±{max_abs_cfo_hz:.0f} Гц"
),
)
axes[0, 1].set_xlabel(
"Eb/N0, дБ"
)
axes[0, 1].set_ylabel(
"FER"
)
axes[0, 1].set_title(
"FER после компенсации CFO"
)
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["frequency_rmse_hz"]
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 RMSE, Гц"
)
axes[1, 0].set_title(
"Среднеквадратичная ошибка оценки CFO"
)
axes[1, 0].grid(
True
)
axes[1, 0].legend()
# ------------------------------------------------------------
# 4. Корреляционная оценка
# ------------------------------------------------------------
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["mean_correlation"]
for result in selected_results
],
marker="o",
label=(
f"CFO ±{max_abs_cfo_hz:.0f} Гц"
),
)
axes[1, 1].axhline(
DETECTION_THRESHOLD,
linestyle=":",
label="Порог обнаружения",
)
axes[1, 1].set_xlabel(
"Eb/N0, дБ"
)
axes[1, 1].set_ylabel(
"Средняя корреляция"
)
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 = [
"Lab022. AWGN and carrier frequency offset",
"",
(
"Trials per point: "
f"{TRIALS_PER_POINT}"
),
(
"Detection threshold: "
f"{DETECTION_THRESHOLD:.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"
),
(
"Success without CFO correction: "
f"{result['success_without_cfo']}"
),
(
"Success with CFO correction: "
f"{result['success_with_cfo']}"
),
(
"FER without correction: "
f"{result['fer_without_cfo']:.6f}"
),
(
"FER with correction: "
f"{result['fer_with_cfo']:.6f}"
),
(
"CFO RMSE: "
f"{result['frequency_rmse_hz']:.3f} Hz"
),
"",
]
)
REPORT_PATH.write_text(
"\n".join(
report_lines
),
encoding="utf-8",
)
# ============================================================
# Автоматические проверки
# ============================================================
result_10_db_no_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
)
)
total_success_without_large_cfo = sum(
result["success_without_cfo"]
for result in results
if result["max_abs_cfo_hz"] == 1000.0
)
total_success_with_large_cfo = sum(
result["success_with_cfo"]
for result in results
if result["max_abs_cfo_hz"] == 1000.0
)
assert len(results) == (
len(EB_N0_VALUES_DB)
* len(MAX_ABS_CFO_VALUES_HZ)
)
assert (
result_10_db_no_cfo["success_with_cfo"]
>= 25
)
assert (
result_10_db_large_cfo["success_with_cfo"]
>= 24
)
assert (
total_success_with_large_cfo
> total_success_without_large_cfo
)
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 исследовано."
)
Reference in New Issue View Git Blame Copy Permalink