1716 lines
37 KiB
Python
1716 lines
37 KiB
Python
"""
|
||
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 исследовано."
|
||
) |