Files
SDR-Rover/experiments/lab043_hackrf_zero_if_compare.py

1142 lines
49 KiB
Python
Raw 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.
"""Один составной опыт сравнения цифрового счётчика около 0 и +50 кГц.
Скрипт передаёт ровно один нециклический буфер Pluto+ и принимает его одним
непрерывным RX-сеансом HackRF. Исходные signed int8 I/Q-байты сохраняются до
преобразования вместе с полученным из них массивом complex64.
"""
from __future__ import annotations
import argparse
import csv
import gc
import hashlib
import json
import math
import sys
import threading
import time
from dataclasses import asdict
from datetime import datetime, timezone
from pathlib import Path
import numpy as np
from scipy.signal import fftconvolve, welch
if __package__ in {None, ""}:
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
from experiments import lab043_hackrf_counter as counter
TONE_OFFSET_HZ = 50_000.0
ZERO_SECONDS = 0.5
TONE_SECONDS = 0.5
GUARD_SECONDS = 0.25
RX_LEADING_SECONDS = 0.75
RX_TRAILING_SECONDS = 0.75
TX_UPLOAD_MARGIN_SECONDS = 10.0
RX_TIMEOUT_SECONDS = 35.0
TONE_SEARCH_HALF_WIDTH_HZ = 5_000.0
TONE_LOCAL_NOISE_GUARD_HZ = 2_000.0
COARSE_SEARCH_SECONDS = 0.25
MATCH_REFERENCE_START_SYMBOL = 10_000
MATCH_REFERENCE_SYMBOL_COUNT = 30_000
CFO_BLOCK_SYMBOLS = 128
RAW_DIRECTORY = Path("data/raw/lab043_hackrf")
PROCESSED_DIRECTORY = Path("data/processed/lab043_hackrf_zero_if_compare")
def _section(
name: str,
start_sample: int,
sample_count: int,
sample_rate_hz: float,
purpose: str,
) -> dict:
end_sample = start_sample + sample_count
return {
"name": name,
"purpose": purpose,
"start_sample": int(start_sample),
"end_sample_exclusive": int(end_sample),
"sample_count": int(sample_count),
"start_seconds": float(start_sample / sample_rate_hz),
"end_seconds": float(end_sample / sample_rate_hz),
"duration_seconds": float(sample_count / sample_rate_hz),
}
def build_composite_waveform(actual_sample_rate_hz: float) -> tuple[np.ndarray, dict]:
"""Построить семь заранее фиксированных участков одного TX-буфера."""
counter_waveform, counter_parameters = counter.build_counter_waveform(
actual_sample_rate_hz
)
zero_samples = int(round(ZERO_SECONDS * actual_sample_rate_hz))
tone_samples = int(round(TONE_SECONDS * actual_sample_rate_hz))
guard_samples = int(round(GUARD_SECONDS * actual_sample_rate_hz))
counts = (
zero_samples,
tone_samples,
guard_samples,
len(counter_waveform),
guard_samples,
len(counter_waveform),
zero_samples,
)
names = (
"leading_zero",
"tone_plus_50khz",
"counter_plus_50khz",
"counter_plus_50khz",
"guard_after_shifted_counter",
"counter_zero_if",
"trailing_zero",
)
purposes = (
"известный участок без полезного сигнала",
"одиночный контрольный тон +50 кГц",
"защитный интервал",
"цифровой счётчик +50 кГц",
"защитный интервал",
"тот же цифровой счётчик около 0 Гц",
"известный участок после сигнала",
)
# Исправить имя третьего участка отдельно, сохраняя явный порядок семи секций.
names = (
names[0],
names[1],
"guard_before_shifted_counter",
names[3],
names[4],
names[5],
names[6],
)
total_samples = int(sum(counts))
waveform = np.zeros(total_samples, dtype=np.complex64)
sections = []
start = 0
for name, count, purpose in zip(names, counts, purposes, strict=True):
sections.append(_section(name, start, count, actual_sample_rate_hz, purpose))
start += count
by_name = {row["name"]: row for row in sections}
tone = by_name["tone_plus_50khz"]
tone_indexes = np.arange(tone["sample_count"], dtype=np.float64)
waveform[tone["start_sample"] : tone["end_sample_exclusive"]] = (
counter.TX_AMPLITUDE
* np.exp(2j * np.pi * TONE_OFFSET_HZ * tone_indexes / actual_sample_rate_hz)
).astype(np.complex64)
shifted = by_name["counter_plus_50khz"]
shifted_indexes = np.arange(shifted["sample_count"], dtype=np.float64)
shifted_carrier = np.exp(
2j * np.pi * TONE_OFFSET_HZ * shifted_indexes / actual_sample_rate_hz
)
waveform[shifted["start_sample"] : shifted["end_sample_exclusive"]] = (
counter_waveform * shifted_carrier
).astype(np.complex64)
centered = by_name["counter_zero_if"]
waveform[centered["start_sample"] : centered["end_sample_exclusive"]] = (
counter_waveform
)
symbol_ids = counter.counter_symbol_ids(0, counter.TRANSMITTED_SYMBOL_COUNT)
symbol_sha256 = hashlib.sha256(symbol_ids.tobytes()).hexdigest()
metadata = {
"sample_rate_hz": float(actual_sample_rate_hz),
"total_sample_count": total_samples,
"total_duration_seconds": float(total_samples / actual_sample_rate_hz),
"tone_offset_hz": TONE_OFFSET_HZ,
"tone_amplitude": counter.TX_AMPLITUDE,
"counter_amplitude": counter.TX_AMPLITUDE,
"counter_symbol_ids_sha256_shifted": symbol_sha256,
"counter_symbol_ids_sha256_centered": symbol_sha256,
"digital_sequences_bit_identical": True,
"sections": sections,
"counter": counter_parameters
| {
"lfsr_polynomial": counter.LFSR_POLYNOMIAL,
"lfsr_initial_state_hex": f"0x{counter.LFSR_INITIAL_STATE:05X}",
"pulse_shape": "rectangular",
},
"peak_amplitude": float(np.max(np.abs(waveform))),
}
return waveform, metadata
def scale_pluto_samples_in_place(waveform: np.ndarray) -> np.ndarray:
"""Преобразовать нормированный complex64-буфер в шкалу Pluto без копии."""
samples = np.asarray(waveform)
if samples.dtype != np.dtype(np.complex64) or samples.ndim != 1 or not len(samples):
raise ValueError("TX-буфер должен быть непустым одномерным complex64")
if not np.all(np.isfinite(samples)) or float(np.max(np.abs(samples))) > 1.0:
raise ValueError("Нормированный TX-буфер недействителен")
np.multiply(samples, np.float32(2**14), out=samples)
return samples
def validate_pluto_working_point(parameters: dict) -> None:
"""Проверить все заданные параметры Pluto+ до единственного TX."""
counter.validate_pluto_working_point(parameters)
if parameters.get("actual_rf_bandwidth_hz") != 200_000:
raise RuntimeError(
"Рабочая точка Pluto+ не подтверждена: фактическая RF bandwidth не 200 кГц"
)
def _hash_file(path: Path) -> str:
digest = hashlib.sha256()
with Path(path).open("rb") as stream:
while block := stream.read(4 * 1024 * 1024):
digest.update(block)
return digest.hexdigest()
def _write_json_atomically(path: Path, document: dict) -> None:
temporary = path.with_suffix(path.suffix + ".tmp")
temporary.write_text(
json.dumps(document, ensure_ascii=False, indent=2, sort_keys=True) + "\n",
encoding="utf-8",
)
temporary.replace(path)
def verify_int8_to_complex64(
int8_path: Path,
iq_path: Path,
chunk_samples: int = 1_000_000,
) -> dict:
"""Проверить signed int8, порядок I/Q и точное преобразование по частям."""
raw = np.memmap(int8_path, dtype=np.int8, mode="r")
iq = np.load(iq_path, mmap_mode="r", allow_pickle=False)
if len(raw) != 2 * len(iq) or len(iq) == 0 or iq.dtype != np.dtype(np.complex64):
raise RuntimeError("Размеры int8 и complex64 raw не согласованы")
for start in range(0, len(iq), chunk_samples):
end = min(len(iq), start + chunk_samples)
codes = np.asarray(raw[2 * start : 2 * end]).reshape(-1, 2).astype(np.float32)
restored = ((codes[:, 0] + 1j * codes[:, 1]) / 128.0).astype(np.complex64)
if not np.array_equal(restored, np.asarray(iq[start:end])):
raise RuntimeError(f"int8 → complex64 не совпал начиная с отсчёта {start}")
return {
"verified": True,
"int8_dtype": "signed int8",
"byte_order": "interleaved I then Q",
"scale_divisor": 128.0,
"complex_dtype": "complex64",
"sample_count": int(len(iq)),
}
def save_dual_raw_before_processing(
base_path: Path,
raw_bytes: bytes,
samples: np.ndarray,
metadata: dict,
) -> tuple[Path, Path, Path, dict]:
"""Сохранить int8 и complex64, перечитать и проверить их до анализа."""
iq = np.asarray(samples)
if not isinstance(raw_bytes, bytes) or len(raw_bytes) != 2 * len(iq):
raise ValueError("Исходный int8 raw не соответствует числу complex64-отсчётов")
if iq.dtype != np.dtype(np.complex64) or iq.ndim != 1 or not np.all(np.isfinite(iq)):
raise ValueError("Complex64 raw недействителен")
if metadata.get("capture_status") != "captured":
raise ValueError("До обработки capture_status должен быть captured")
if metadata.get("processing_status") != "pending":
raise ValueError("До обработки processing_status должен быть pending")
base = Path(base_path)
base.parent.mkdir(parents=True, exist_ok=True)
int8_path = base.with_suffix(".int8")
iq_path = base.with_suffix(".npy")
metadata_path = base.with_suffix(".json")
int8_path.write_bytes(raw_bytes)
np.save(iq_path, iq, allow_pickle=False)
document = dict(metadata)
document.update(
{
"int8_file": int8_path.name,
"int8_byte_count": len(raw_bytes),
"int8_sha256": None,
"iq_file": iq_path.name,
"sample_count": int(len(iq)),
"dtype": "complex64",
"iq_sha256": None,
"raw_verified": False,
"int8_to_complex64": None,
}
)
_write_json_atomically(metadata_path, document)
initial_int8_sha = hashlib.sha256(raw_bytes).hexdigest()
int8_sha = _hash_file(int8_path)
iq_sha = _hash_file(iq_path)
if int8_sha != initial_int8_sha or _hash_file(int8_path) != int8_sha:
raise RuntimeError("Повторная SHA-256 проверка int8 raw не прошла")
if _hash_file(iq_path) != iq_sha:
raise RuntimeError("Повторная SHA-256 проверка complex64 raw не прошла")
conversion = verify_int8_to_complex64(int8_path, iq_path)
document.update(
{
"int8_sha256": int8_sha,
"iq_sha256": iq_sha,
"raw_verified": True,
"int8_to_complex64": conversion,
}
)
_write_json_atomically(metadata_path, document)
return int8_path, iq_path, metadata_path, document
def _code_distribution(codes: np.ndarray) -> dict[str, int]:
values, counts = np.unique(np.asarray(codes, dtype=np.int8), return_counts=True)
return {str(int(value)): int(count) for value, count in zip(values, counts, strict=True)}
def background_metrics(samples: np.ndarray, codes: np.ndarray, sample_rate_hz: float) -> dict:
"""Измерить фон без вычитания постоянного комплексного смещения."""
values = np.asarray(samples, dtype=np.complex64)
pairs = np.asarray(codes, dtype=np.int8)
if values.ndim != 1 or pairs.shape != (len(values), 2) or not len(values):
raise ValueError("Фоновый int8/complex64 интервал недействителен")
frequencies, density = welch(
values,
fs=sample_rate_hz,
window="hann",
nperseg=min(65_536, len(values)),
noverlap=0,
return_onesided=False,
scaling="density",
)
finite_density = density[np.isfinite(density) & (density > 0.0)]
mean_value = complex(np.mean(values, dtype=np.complex128))
mean_power = float(np.mean(np.abs(values.astype(np.complex128)) ** 2))
return {
"sample_count": int(len(values)),
"i_code_distribution": _code_distribution(pairs[:, 0]),
"q_code_distribution": _code_distribution(pairs[:, 1]),
"mean_i": float(mean_value.real),
"mean_q": float(mean_value.imag),
"complex_dc": {"i": float(mean_value.real), "q": float(mean_value.imag)},
"rms": float(math.sqrt(mean_power)),
"mean_power": mean_power,
"median_noise_psd_dbfs_per_hz": float(
10.0 * math.log10(max(float(np.median(finite_density)), 1e-300))
),
"welch_bin_hz": float(sample_rate_hz / min(65_536, len(values))),
"dc_correction_applied": False,
}
def _parabolic_frequency(frequencies: np.ndarray, powers: np.ndarray, peak: int) -> float:
if peak <= 0 or peak >= len(powers) - 1:
return float(frequencies[peak])
y0, y1, y2 = np.log(np.maximum(powers[peak - 1 : peak + 2], 1e-300))
denominator = y0 - 2.0 * y1 + y2
fraction = 0.0 if abs(denominator) < 1e-18 else 0.5 * (y0 - y2) / denominator
return float(frequencies[peak] + fraction * (frequencies[1] - frequencies[0]))
def locate_and_measure_tone(
samples: np.ndarray,
sample_rate_hz: float,
expected_start_sample: int,
expected_count: int,
) -> dict:
"""Найти фактический тон около +50 кГц в априорном временном окне."""
search_margin = int(round(0.35 * sample_rate_hz))
search_start = max(0, expected_start_sample - search_margin)
search_end = min(len(samples), expected_start_sample + expected_count + search_margin)
search = np.asarray(samples[search_start:search_end], dtype=np.complex64)
fft_samples = min(524_288, 1 << (len(search).bit_length() - 1))
if fft_samples < 65_536:
raise RuntimeError("Недостаточно raw для поиска контрольного тона")
hop = max(fft_samples // 8, 1)
window = np.hanning(fft_samples).astype(np.float32)
frequencies = np.fft.fftfreq(fft_samples, d=1.0 / sample_rate_hz)
tone_band = np.flatnonzero(
np.abs(frequencies - TONE_OFFSET_HZ) <= TONE_SEARCH_HALF_WIDTH_HZ
)
best_candidate: tuple[float, int, int, np.ndarray] | None = None
for start in range(0, len(search) - fft_samples + 1, hop):
power = np.abs(np.fft.fft(search[start : start + fft_samples] * window)) ** 2
peak = int(tone_band[np.argmax(power[tone_band])])
candidate = (float(power[peak]), start, peak, power)
if best_candidate is None or candidate[0] > best_candidate[0]:
best_candidate = candidate
if best_candidate is None:
raise RuntimeError("В априорном окне нет полного FFT-блока контрольного тона")
_, coarse_start, coarse_peak, coarse_power = best_candidate
frequency_hz = _parabolic_frequency(frequencies, coarse_power, coarse_peak)
detector_samples = 8_192
detector_hop = 2_048
detector_scores = []
detector_starts = []
detector_indexes = np.arange(detector_samples, dtype=np.float64)
detector_oscillator = np.exp(
-2j * np.pi * frequency_hz * detector_indexes / sample_rate_hz
)
for start in range(0, len(search) - detector_samples + 1, detector_hop):
gain = np.mean(
search[start : start + detector_samples].astype(np.complex128)
* detector_oscillator
* np.exp(-2j * np.pi * frequency_hz * start / sample_rate_hz)
)
detector_starts.append(start)
detector_scores.append(float(abs(gain) ** 2))
scores = np.asarray(detector_scores, dtype=np.float64)
baseline = float(np.percentile(scores, 20.0))
threshold = baseline + 0.35 * (float(np.max(scores)) - baseline)
active = scores >= threshold
best_run: tuple[int, int] | None = None
run_start = None
for index, is_active in enumerate(np.append(active, False)):
if is_active and run_start is None:
run_start = index
elif not is_active and run_start is not None:
if best_run is None or index - run_start > best_run[1] - best_run[0]:
best_run = (run_start, index)
run_start = None
if best_run is None:
raise RuntimeError("Не удалось выделить временной интервал контрольного тона")
detected_start = search_start + detector_starts[best_run[0]]
detected_end = (
search_start
+ detector_starts[best_run[1] - 1]
+ detector_samples
)
center_sample = (detected_start + detected_end) // 2
tone_start = int(round(center_sample - expected_count / 2.0))
tone_start = int(np.clip(tone_start, 0, len(samples) - expected_count))
tone_end = tone_start + expected_count
tone = np.asarray(samples[tone_start:tone_end], dtype=np.complex64)
indexes = np.arange(expected_count, dtype=np.float64)
oscillator = np.exp(-2j * np.pi * frequency_hz * indexes / sample_rate_hz)
gain = complex(np.mean(tone.astype(np.complex128) * oscillator))
predicted = gain * np.conj(oscillator)
residual_power = float(np.mean(np.abs(tone.astype(np.complex128) - predicted) ** 2))
spectrum = np.abs(np.fft.fft(tone[:fft_samples] * window)) ** 2
local = (
np.abs(frequencies - frequency_hz) <= TONE_SEARCH_HALF_WIDTH_HZ
) & (np.abs(frequencies - frequency_hz) >= TONE_LOCAL_NOISE_GUARD_HZ)
local_noise = float(np.median(spectrum[local]))
peak_power = float(spectrum[int(np.argmin(np.abs(frequencies - frequency_hz)))])
harmonics = {}
for order in (2, 3):
harmonic_frequency = order * frequency_hz
if abs(harmonic_frequency) >= sample_rate_hz / 2.0:
continue
harmonic = complex(
np.mean(
tone.astype(np.complex128)
* np.exp(-2j * np.pi * harmonic_frequency * indexes / sample_rate_hz)
)
)
harmonics[str(order)] = {
"frequency_hz": float(harmonic_frequency),
"amplitude": float(abs(harmonic)),
"relative_to_tone_db": float(
20.0 * math.log10(max(abs(harmonic), 1e-300) / max(abs(gain), 1e-300))
),
}
component = np.maximum(np.abs(tone.real), np.abs(tone.imag))
dc = complex(np.mean(tone, dtype=np.complex128))
return {
"start_sample": tone_start,
"end_sample_exclusive": tone_end,
"frequency_hz": frequency_hz,
"frequency_offset_from_requested_hz": float(frequency_hz - TONE_OFFSET_HZ),
"amplitude": float(abs(gain)),
"power": float(abs(gain) ** 2),
"residual_power_without_dc_correction": residual_power,
"tone_to_residual_db": float(
10.0 * math.log10(max(abs(gain) ** 2, 1e-300) / max(residual_power, 1e-300))
),
"peak_above_local_noise_db": float(
10.0 * math.log10(max(peak_power, 1e-300) / max(local_noise, 1e-300))
),
"complex_dc": {"i": float(dc.real), "q": float(dc.imag)},
"harmonics": harmonics,
"maximum_component": float(np.max(component)),
"clipped_component_fraction": float(np.mean(component >= 127.0 / 128.0)),
"dc_correction_applied": False,
}
def _mix_samples(
samples: np.ndarray,
absolute_start_sample: int,
frequency_hz: float,
sample_rate_hz: float,
) -> np.ndarray:
values = np.asarray(samples, dtype=np.complex64)
result = np.empty(len(values), dtype=np.complex64)
chunk = 1_000_000
for start in range(0, len(values), chunk):
end = min(len(values), start + chunk)
indexes = absolute_start_sample + np.arange(start, end, dtype=np.float64)
result[start:end] = values[start:end] * np.exp(
-2j * np.pi * frequency_hz * indexes / sample_rate_hz
)
return result
def _integrate_symbols(
samples: np.ndarray,
start_sample: int,
symbol_count: int,
carrier_hz: float,
sample_rate_hz: float,
) -> np.ndarray:
sample_count = symbol_count * counter.NOMINAL_SAMPLES_PER_SYMBOL
mixed = _mix_samples(
samples[start_sample : start_sample + sample_count],
start_sample,
carrier_hz,
sample_rate_hz,
)
symbols = mixed.reshape(symbol_count, counter.NOMINAL_SAMPLES_PER_SYMBOL).mean(axis=1)
return symbols.astype(np.complex64)
def _find_counter_start(
samples: np.ndarray,
expected_start_sample: int,
carrier_hz: float,
sample_rate_hz: float,
) -> int:
margin = int(round(COARSE_SEARCH_SECONDS * sample_rate_hz))
region_start = max(0, expected_start_sample - margin)
region_end = min(
len(samples),
expected_start_sample
+ counter.TRANSMITTED_SYMBOL_COUNT * counter.NOMINAL_SAMPLES_PER_SYMBOL
+ margin,
)
region_count = (region_end - region_start) // counter.NOMINAL_SAMPLES_PER_SYMBOL
mixed = _mix_samples(
samples[
region_start : region_start
+ region_count * counter.NOMINAL_SAMPLES_PER_SYMBOL
],
region_start,
carrier_hz,
sample_rate_hz,
)
coarse_symbols = mixed.reshape(
region_count, counter.NOMINAL_SAMPLES_PER_SYMBOL
).mean(axis=1)
del mixed
reference = counter.counter_qpsk_symbols(
MATCH_REFERENCE_START_SYMBOL,
MATCH_REFERENCE_SYMBOL_COUNT,
).astype(np.complex128)
correlation = fftconvolve(
coarse_symbols.astype(np.complex128),
np.conj(reference[::-1]),
mode="valid",
)
match_index = int(np.argmax(np.abs(correlation)))
coarse_start = (
region_start
+ (match_index - MATCH_REFERENCE_START_SYMBOL)
* counter.NOMINAL_SAMPLES_PER_SYMBOL
)
candidates = np.arange(
coarse_start - counter.NOMINAL_SAMPLES_PER_SYMBOL + 1,
coarse_start + counter.NOMINAL_SAMPLES_PER_SYMBOL,
dtype=np.int64,
)
subset_starts = (
candidates
+ MATCH_REFERENCE_START_SYMBOL * counter.NOMINAL_SAMPLES_PER_SYMBOL
)
valid = (subset_starts >= 0) & (
subset_starts
+ MATCH_REFERENCE_SYMBOL_COUNT * counter.NOMINAL_SAMPLES_PER_SYMBOL
<= len(samples)
)
candidates = candidates[valid]
subset_starts = subset_starts[valid]
if not len(candidates):
raise RuntimeError("Не удалось уточнить начало цифрового участка")
mixed_start = int(np.min(subset_starts))
mixed_end = int(
np.max(subset_starts)
+ MATCH_REFERENCE_SYMBOL_COUNT * counter.NOMINAL_SAMPLES_PER_SYMBOL
)
mixed = _mix_samples(
samples[mixed_start:mixed_end],
mixed_start,
carrier_hz,
sample_rate_hz,
)
cumulative = np.concatenate(
(np.zeros(1, dtype=np.complex128), np.cumsum(mixed, dtype=np.complex128))
)
symbol_offsets = (
np.arange(MATCH_REFERENCE_SYMBOL_COUNT, dtype=np.int64)
* counter.NOMINAL_SAMPLES_PER_SYMBOL
)
best: tuple[float, int] | None = None
for candidate, subset_start in zip(candidates, subset_starts, strict=True):
starts = int(subset_start - mixed_start) + symbol_offsets
symbols = (
cumulative[starts + counter.NOMINAL_SAMPLES_PER_SYMBOL]
- cumulative[starts]
) / counter.NOMINAL_SAMPLES_PER_SYMBOL
gain = float(abs(np.mean(symbols * np.conj(reference))))
row = (gain, -int(candidate))
if best is None or row > best:
best = row
if best is None:
raise RuntimeError("Не удалось уточнить начало цифрового участка")
return -best[1]
def _constellation_metrics(symbols: np.ndarray, expected_ids: np.ndarray) -> dict:
centroids = []
residual_sum = 0.0
for state in range(4):
group = symbols[expected_ids == state]
centroid = complex(np.mean(group, dtype=np.complex128))
centroids.append(centroid)
residual_sum += float(np.sum(np.abs(group - centroid) ** 2))
distances = [
abs(centroids[left] - centroids[right])
for left in range(4)
for right in range(left + 1, 4)
]
within_rms = math.sqrt(residual_sum / len(symbols))
minimum_distance = float(min(distances))
return {
"centroids": [
{"i": float(value.real), "q": float(value.imag)} for value in centroids
],
"minimum_centroid_distance": minimum_distance,
"pooled_within_state_rms": float(within_rms),
"centroid_separation_to_within_rms": float(
minimum_distance / max(within_rms, 1e-300)
),
}
def analyze_counter_section(
samples: np.ndarray,
expected_start_sample: int,
prior_carrier_hz: float,
sample_rate_hz: float,
background_power: float,
) -> dict:
"""Одинаково обработать один известный цифровой участок без DC-коррекции."""
start = _find_counter_start(
samples,
expected_start_sample,
prior_carrier_hz,
sample_rate_hz,
)
end = start + counter.TRANSMITTED_SYMBOL_COUNT * counter.NOMINAL_SAMPLES_PER_SYMBOL
if start < 0 or end > len(samples):
raise RuntimeError("Полный цифровой участок не помещается в raw")
symbols = _integrate_symbols(
samples,
start,
counter.TRANSMITTED_SYMBOL_COUNT,
prior_carrier_hz,
sample_rate_hz,
).astype(np.complex128)
expected_ids = counter.counter_symbol_ids(0, counter.TRANSMITTED_SYMBOL_COUNT)
expected = counter.QPSK_STATES[expected_ids].astype(np.complex128)
matched = symbols * np.conj(expected)
symbol_rate = sample_rate_hz / counter.NOMINAL_SAMPLES_PER_SYMBOL
time_axis = (np.arange(len(symbols), dtype=np.float64) + 0.5) / symbol_rate
power = np.abs(np.fft.fft(matched * np.hanning(len(matched)))) ** 2
frequencies = np.fft.fftfreq(len(matched), d=1.0 / symbol_rate)
peak = int(np.argmax(power))
residual_initial = _parabolic_frequency(frequencies, power, peak)
block_count = len(matched) // CFO_BLOCK_SYMBOLS
initial = matched * np.exp(-2j * np.pi * residual_initial * time_axis)
block_products = initial[: block_count * CFO_BLOCK_SYMBOLS].reshape(
block_count, CFO_BLOCK_SYMBOLS
).mean(axis=1)
block_times = time_axis[: block_count * CFO_BLOCK_SYMBOLS].reshape(
block_count, CFO_BLOCK_SYMBOLS
).mean(axis=1)
slope, _ = np.polyfit(block_times, np.unwrap(np.angle(block_products)), 1)
residual_cfo = float(residual_initial + slope / (2.0 * np.pi))
derotation = np.exp(-2j * np.pi * residual_cfo * time_axis)
gain = complex(np.mean(matched * derotation))
predicted = gain * expected * np.exp(2j * np.pi * residual_cfo * time_axis)
residual_power = float(np.mean(np.abs(symbols - predicted) ** 2))
signal_power = float(abs(gain) ** 2)
normalized = symbols * derotation / gain
received_ids = np.argmin(
np.abs(normalized[:, None] - counter.QPSK_STATES.astype(np.complex128)[None, :]),
axis=1,
).astype(np.uint8)
symbol_errors = int(np.count_nonzero(received_ids != expected_ids))
position = counter.analyze_counter_symbols(received_ids)
section = np.asarray(samples[start:end], dtype=np.complex64)
dc = complex(np.mean(section, dtype=np.complex128))
total_power = float(np.mean(np.abs(section.astype(np.complex128)) ** 2))
excess_power = total_power - background_power
corrected_symbols = normalized.astype(np.complex64)
return {
"start_sample": int(start),
"end_sample_exclusive": int(end),
"prior_carrier_hz": float(prior_carrier_hz),
"measured_carrier_hz": float(prior_carrier_hz + residual_cfo),
"residual_cfo_hz": residual_cfo,
"complex_gain": {"i": float(gain.real), "q": float(gain.imag)},
"useful_signal_power_symbol": signal_power,
"noise_and_distortion_power_symbol": residual_power,
"symbol_snr_db": float(
10.0 * math.log10(max(signal_power, 1e-300) / max(residual_power, 1e-300))
),
"section_total_power": total_power,
"background_power": float(background_power),
"sample_excess_power": float(excess_power) if excess_power > 0.0 else None,
"sample_power_snr_db": (
float(10.0 * math.log10(excess_power / background_power))
if excess_power > 0.0
else None
),
"sample_power_snr_status": "computed" if excess_power > 0.0 else "not_computed",
"complex_dc": {"i": float(dc.real), "q": float(dc.imag)},
"constellation": _constellation_metrics(corrected_symbols, expected_ids),
"symbol_error_count": symbol_errors,
"symbol_error_rate": float(symbol_errors / len(expected_ids)),
"recovered_position_count": int(position["recovered_symbol_count"]),
"anchor_count": int(position["anchor_count"]),
"position_map_available": bool(position["recovered_symbol_count"] > 0),
"position_analysis": position,
"dc_correction_applied": False,
"new_filter_applied": False,
"equalizer_applied": False,
}
def _mapped_section_bounds(section: dict, raw_tx_start: int, tx_fs: float, rx_fs: float) -> tuple[int, int]:
start = raw_tx_start + int(round(section["start_sample"] / tx_fs * rx_fs))
end = raw_tx_start + int(round(section["end_sample_exclusive"] / tx_fs * rx_fs))
return start, end
def analyze_saved_capture(
samples: np.ndarray,
raw_codes: np.ndarray,
document: dict,
) -> dict:
"""Выполнить весь анализ только после проверки обоих raw-файлов."""
rx_fs = float(document["rx_parameters"]["set_sample_rate_hz"])
tx_fs = float(document["tx_parameters"]["actual_sample_rate_hz"])
sections = {row["name"]: row for row in document["waveform"]["sections"]}
boundaries = document["callback_boundaries"]
stream_start = float(document["continuous_rx"]["stream_started_monotonic"])
estimated_tx_return = (
float(document["tx_timing"]["tx_completed_monotonic"])
- float(document["tx_timing"]["tx_buffer_hold_seconds"])
)
relative = estimated_tx_return - stream_start
times = np.asarray([row["monotonic_seconds"] for row in boundaries], dtype=np.float64)
ends = np.asarray([row["accepted_end_sample"] for row in boundaries], dtype=np.float64)
expected_tx_start = int(round(np.interp(relative, times, ends)))
tone_section = sections["tone_plus_50khz"]
expected_tone_start, expected_tone_end = _mapped_section_bounds(
tone_section, expected_tx_start, tx_fs, rx_fs
)
tone = locate_and_measure_tone(
samples,
rx_fs,
expected_tone_start,
expected_tone_end - expected_tone_start,
)
raw_tx_start = tone["start_sample"] - int(
round(tone_section["start_sample"] / tx_fs * rx_fs)
)
raw_tx_end = raw_tx_start + int(
round(document["waveform"]["total_sample_count"] / tx_fs * rx_fs)
)
if raw_tx_start < 0 or raw_tx_end > len(samples):
raise RuntimeError("Полный составной TX-буфер не помещается в raw")
quiet_margin = int(round(0.05 * rx_fs))
quiet_count = int(round(0.35 * rx_fs))
pre_start = max(0, raw_tx_start - quiet_margin - quiet_count)
pre_end = max(pre_start, raw_tx_start - quiet_margin)
post_start = min(len(samples), raw_tx_end + quiet_margin)
post_end = min(len(samples), post_start + quiet_count)
if pre_end <= pre_start or post_end <= post_start:
raise RuntimeError("Не сохранились оба фоновых участка до и после TX")
background_samples = np.concatenate(
(np.asarray(samples[pre_start:pre_end]), np.asarray(samples[post_start:post_end]))
)
pairs = np.asarray(raw_codes).reshape(-1, 2)
background_codes = np.concatenate(
(np.asarray(pairs[pre_start:pre_end]), np.asarray(pairs[post_start:post_end]))
)
background = background_metrics(background_samples, background_codes, rx_fs)
shifted_expected, _ = _mapped_section_bounds(
sections["counter_plus_50khz"], raw_tx_start, tx_fs, rx_fs
)
centered_expected, _ = _mapped_section_bounds(
sections["counter_zero_if"], raw_tx_start, tx_fs, rx_fs
)
shifted = analyze_counter_section(
samples,
shifted_expected,
tone["frequency_hz"],
rx_fs,
background["mean_power"],
)
centered = analyze_counter_section(
samples,
centered_expected,
tone["frequency_hz"] - TONE_OFFSET_HZ,
rx_fs,
background["mean_power"],
)
snr_improvement = shifted["symbol_snr_db"] - centered["symbol_snr_db"]
ser_improvement = centered["symbol_error_rate"] - shifted["symbol_error_rate"]
if snr_improvement >= 3.0 and ser_improvement >= 0.05:
interpretation = "A"
conclusion = "Перенос на +50 кГц значительно улучшил приём относительно центра"
elif abs(snr_improvement) < 1.5 and max(
shifted["symbol_error_rate"], centered["symbol_error_rate"]
) >= 0.25:
interpretation = "B"
conclusion = "Оба цифровых участка одинаково плохи; основной предел остаётся уровневым"
elif max(shifted["symbol_error_rate"], centered["symbol_error_rate"]) <= 0.02:
interpretation = "C"
conclusion = "Оба цифровых участка достоверно декодируются"
else:
interpretation = "D"
conclusion = "Результат прямого сравнения неоднозначен"
return {
"raw_tx_start_sample": int(raw_tx_start),
"raw_tx_end_sample_exclusive": int(raw_tx_end),
"raw_tx_start_seconds": float(raw_tx_start / rx_fs),
"raw_tx_end_seconds": float(raw_tx_end / rx_fs),
"background_intervals": [
{"start_sample": pre_start, "end_sample_exclusive": pre_end},
{"start_sample": post_start, "end_sample_exclusive": post_end},
],
"background": background,
"tone_plus_50khz": tone,
"counter_plus_50khz": shifted,
"counter_zero_if": centered,
"comparison": {
"symbol_snr_improvement_plus_50khz_db": float(snr_improvement),
"symbol_error_rate_improvement_plus_50khz": float(ser_improvement),
"interpretation": interpretation,
"conclusion": conclusion,
"zero_if_degradation_proven": interpretation == "A",
"insufficient_level_remains_primary": interpretation == "B",
},
"usb_loss_analysis": "not_performed",
"dc_correction_applied": False,
"new_filter_applied": False,
"equalizer_applied": False,
}
def save_processed_artifacts(analysis: dict, output_directory: Path) -> tuple[Path, ...]:
output = Path(output_directory)
output.mkdir(parents=True, exist_ok=True)
_write_json_atomically(output / "analysis.json", analysis)
shifted = analysis["counter_plus_50khz"]
centered = analysis["counter_zero_if"]
rows = [
("мощность полезного сигнала на символ", shifted["useful_signal_power_symbol"], centered["useful_signal_power_symbol"]),
("мощность шума и искажений на символ", shifted["noise_and_distortion_power_symbol"], centered["noise_and_distortion_power_symbol"]),
("SNR на символ, дБ", shifted["symbol_snr_db"], centered["symbol_snr_db"]),
("разделимость состояний", shifted["constellation"]["centroid_separation_to_within_rms"], centered["constellation"]["centroid_separation_to_within_rms"]),
("доля ошибочных символов", shifted["symbol_error_rate"], centered["symbol_error_rate"]),
("восстановленные позиции", shifted["recovered_position_count"], centered["recovered_position_count"]),
("постоянное смещение I", shifted["complex_dc"]["i"], centered["complex_dc"]["i"]),
("постоянное смещение Q", shifted["complex_dc"]["q"], centered["complex_dc"]["q"]),
]
with (output / "comparison.csv").open("w", newline="", encoding="utf-8-sig") as stream:
writer = csv.writer(stream)
writer.writerow(("показатель", "+50 кГц", "около 0 Гц"))
writer.writerows(rows)
report = [
"Lab043: прямое сравнение цифрового счётчика +50 кГц и около 0 Гц",
"",
f"Интерпретация: {analysis['comparison']['interpretation']}.",
analysis["comparison"]["conclusion"] + ".",
f"SNR на символ +50 кГц: {shifted['symbol_snr_db']:.3f} дБ.",
f"SNR на символ около 0 Гц: {centered['symbol_snr_db']:.3f} дБ.",
f"Доля ошибочных символов +50 кГц: {shifted['symbol_error_rate']:.6f}.",
f"Доля ошибочных символов около 0 Гц: {centered['symbol_error_rate']:.6f}.",
f"Восстановленные позиции +50 кГц: {shifted['recovered_position_count']}.",
f"Восстановленные позиции около 0 Гц: {centered['recovered_position_count']}.",
"",
"SNR означает отношение мощности сигнала к мощности шума и искажений.",
"I/Q означает синфазные и квадратурные комплексные отсчёты.",
"Постоянное смещение не вычиталось; эквалайзер и новый фильтр не применялись.",
]
(output / "report.txt").write_text("\n".join(report) + "\n", encoding="utf-8")
return tuple(sorted(output.iterdir()))
def update_processing_metadata(
metadata_path: Path,
processing_status: str,
processing_error: str | None,
analysis: dict | None = None,
) -> dict:
if processing_status not in {"success", "processing_failed"}:
raise ValueError("Неизвестный processing_status")
document = json.loads(Path(metadata_path).read_text(encoding="utf-8"))
document["processing_status"] = processing_status
document["processing_error"] = processing_error
if analysis is not None:
document["analysis"] = analysis
_write_json_atomically(Path(metadata_path), document)
return document
def preflight_hardware() -> dict:
"""Открыть оба SDR без TX и проверить единственную разрешённую точку."""
from experiments import lab043_hackrf_hardware as hackrf
from experiments import lab043_prbs11_hardware as pluto_hardware
pluto_uri = pluto_hardware.discover_pluto_usb_uri()
return {
"pluto": pluto_hardware.preflight_pluto(pluto_uri),
"hackrf": counter.preflight_counter_receiver(hackrf),
}
def run_one_zero_if_comparison(
raw_directory: Path = RAW_DIRECTORY,
processed_directory: Path = PROCESSED_DIRECTORY,
) -> dict:
"""Выполнить ровно один разрешённый составной аппаратный опыт."""
from experiments import lab043_hackrf_hardware as hackrf
from experiments import lab043_prbs11_hardware as pluto_hardware
pluto_uri = pluto_hardware.discover_pluto_usb_uri()
pluto_preflight = pluto_hardware.preflight_pluto(pluto_uri)
hackrf_preflight = counter.preflight_counter_receiver(hackrf)
pluto, tx_parameters = pluto_hardware.configure_pluto(pluto_uri)
validate_pluto_working_point(tx_parameters)
actual_tx_rate = float(tx_parameters["actual_sample_rate_hz"])
waveform, waveform_parameters = build_composite_waveform(actual_tx_rate)
waveform_sha256 = hashlib.sha256(waveform.tobytes()).hexdigest()
scale_pluto_samples_in_place(waveform)
target_samples = math.ceil(
(
RX_LEADING_SECONDS
+ waveform_parameters["total_duration_seconds"]
+ RX_TRAILING_SECONDS
+ TX_UPLOAD_MARGIN_SECONDS
)
* hackrf.HACKRF_SAMPLE_RATE_HZ
) + hackrf.HACKRF_EXPECTED_CALLBACK_BYTES // 2
capture_started = datetime.now(timezone.utc)
capture_holder = []
errors: list[str] = []
tx_timing: dict[str, float] = {}
with hackrf.HackrfRxDevice() as receiver:
rx_parameters = counter.configure_counter_receiver(receiver)
def receive() -> None:
try:
capture_holder.append(
receiver.capture(target_samples, timeout_seconds=RX_TIMEOUT_SECONDS)
)
except Exception as error: # pragma: no cover - аппаратный путь
errors.append(f"RX {type(error).__name__}: {error}")
rx_thread = threading.Thread(target=receive, daemon=True)
rx_thread.start()
deadline = time.perf_counter() + 5.0
stream_confirmed = False
while time.perf_counter() < deadline:
if errors:
raise RuntimeError(
"HackRF RX завершился ошибкой до TX; Pluto TX не включён: "
+ "; ".join(errors)
)
if receiver.library.hackrf_is_streaming(receiver.device) == 1:
stream_confirmed = True
break
time.sleep(0.01)
if not stream_confirmed:
raise RuntimeError("HackRF RX не подтвердил поток; Pluto TX не включён")
time.sleep(RX_LEADING_SECONDS)
if errors:
raise RuntimeError(
"HackRF RX завершился ошибкой до TX; Pluto TX не включён: "
+ "; ".join(errors)
)
if receiver.library.hackrf_is_streaming(receiver.device) != 1:
raise RuntimeError("HackRF RX остановился до TX; Pluto TX не включён")
tx_timing["tx_call_started_monotonic"] = time.perf_counter()
tx_timing["tx_buffer_hold_seconds"] = pluto_hardware.transmit_noncyclic_buffer_once(
pluto,
waveform,
actual_tx_rate,
)
tx_timing["tx_completed_monotonic"] = time.perf_counter()
del waveform
gc.collect()
rx_thread.join(timeout=RX_TIMEOUT_SECONDS + 10.0)
if rx_thread.is_alive():
raise RuntimeError("Единственный непрерывный HackRF RX-сеанс не завершился")
if errors:
raise RuntimeError("; ".join(errors))
if len(capture_holder) != 1:
raise RuntimeError("HackRF не вернул единственный захват")
capture = capture_holder[0]
if (
len(capture.samples) != target_samples
or len(capture.raw_bytes) != 2 * target_samples
or not np.all(np.isfinite(capture.samples))
):
raise RuntimeError("Размер или содержимое raw RX недействительны")
boundary_rows = [asdict(row) for row in capture.boundaries]
metadata = {
"experiment": "Lab043 HackRF zero-IF direct comparison",
"capture_status": "captured",
"processing_status": "pending",
"processing_error": None,
"capture_started_utc": capture_started.isoformat().replace("+00:00", "Z"),
"physical_scheme": "Pluto+ TX -> SMA -> AT30S 30 dB -> HackRF ANTENNA",
"tx_parameters": tx_parameters,
"rx_parameters": rx_parameters
| {
"serial_number": hackrf_preflight["serial_number"],
"firmware_version": hackrf_preflight["firmware_version"],
"libhackrf_version": hackrf_preflight["libhackrf_version"],
},
"pluto_preflight": pluto_preflight,
"hackrf_preflight": hackrf_preflight,
"waveform_sha256_normalized_complex64": waveform_sha256,
"waveform": waveform_parameters,
"tx_waveform_duration_seconds": waveform_parameters["total_duration_seconds"],
"tx_timing": tx_timing,
"callback_count": capture.diagnostics["callback_count_total"],
"callback_boundaries": boundary_rows,
"continuous_rx": capture.diagnostics | {"boundaries": boundary_rows},
"expected_intervals": {
"rx_leading_seconds": RX_LEADING_SECONDS,
"rx_trailing_seconds": RX_TRAILING_SECONDS,
"tx_upload_margin_seconds": TX_UPLOAD_MARGIN_SECONDS,
},
"forbidden_modes_executed": [],
"automatic_gain_change": False,
"dc_correction_applied": False,
}
stamp = capture_started.strftime("%Y%m%d_%H%M%S")
base = Path(raw_directory) / f"lab043_hackrf_zero_if_compare_{stamp}"
int8_path, iq_path, metadata_path, document = save_dual_raw_before_processing(
base,
capture.raw_bytes,
capture.samples,
metadata,
)
del capture_holder, capture
gc.collect()
samples = np.load(iq_path, mmap_mode="r", allow_pickle=False)
raw_codes = np.memmap(int8_path, dtype=np.int8, mode="r")
if _hash_file(int8_path) != document["int8_sha256"]:
raise RuntimeError("SHA-256 int8 raw изменился до анализа")
if _hash_file(iq_path) != document["iq_sha256"]:
raise RuntimeError("SHA-256 complex64 raw изменился до анализа")
try:
analysis = analyze_saved_capture(samples, raw_codes, document)
save_processed_artifacts(analysis, processed_directory)
except Exception as error:
message = f"{type(error).__name__}: {error}"
update_processing_metadata(metadata_path, "processing_failed", message)
return {
"int8_path": str(int8_path),
"iq_path": str(iq_path),
"metadata_path": str(metadata_path),
"int8_sha256": document["int8_sha256"],
"iq_sha256": document["iq_sha256"],
"processing_status": "processing_failed",
"processing_error": message,
}
update_processing_metadata(metadata_path, "success", None, analysis)
return {
"int8_path": str(int8_path),
"iq_path": str(iq_path),
"metadata_path": str(metadata_path),
"int8_sha256": document["int8_sha256"],
"iq_sha256": document["iq_sha256"],
"processing_status": "success",
"tx_parameters": tx_parameters,
"rx_parameters": metadata["rx_parameters"],
"analysis": analysis,
}
def reprocess_saved_capture(
metadata_path: Path,
processed_directory: Path = PROCESSED_DIRECTORY,
) -> dict:
"""Повторить только анализ сохранённых raw без открытия SDR и без TX."""
path = Path(metadata_path)
document = json.loads(path.read_text(encoding="utf-8"))
int8_path = path.parent / document["int8_file"]
iq_path = path.parent / document["iq_file"]
if not document.get("raw_verified"):
raise RuntimeError("Raw не отмечены как проверенные")
if _hash_file(int8_path) != document.get("int8_sha256"):
raise RuntimeError("SHA-256 int8 raw не совпал")
if _hash_file(iq_path) != document.get("iq_sha256"):
raise RuntimeError("SHA-256 complex64 raw не совпал")
verify_int8_to_complex64(int8_path, iq_path)
samples = np.load(iq_path, mmap_mode="r", allow_pickle=False)
raw_codes = np.memmap(int8_path, dtype=np.int8, mode="r")
try:
analysis = analyze_saved_capture(samples, raw_codes, document)
save_processed_artifacts(analysis, processed_directory)
except Exception as error:
message = f"{type(error).__name__}: {error}"
update_processing_metadata(path, "processing_failed", message)
return {"processing_status": "processing_failed", "processing_error": message}
update_processing_metadata(path, "success", None, analysis)
return {"processing_status": "success", "analysis": analysis}
def main() -> int:
parser = argparse.ArgumentParser()
group = parser.add_mutually_exclusive_group(required=True)
group.add_argument("--preflight-only", action="store_true")
group.add_argument("--run-one-zero-if-compare", action="store_true")
group.add_argument("--reprocess-metadata", type=Path)
arguments = parser.parse_args()
if arguments.preflight_only:
print(json.dumps(preflight_hardware(), ensure_ascii=False, indent=2))
return 0
result = (
reprocess_saved_capture(arguments.reprocess_metadata)
if arguments.reprocess_metadata is not None
else run_one_zero_if_comparison()
)
print(json.dumps(result, ensure_ascii=False, indent=2))
return 0 if result["processing_status"] == "success" else 2
if __name__ == "__main__":
raise SystemExit(main())