"""Один составной опыт сравнения цифрового счётчика около 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())