Lab052: suppress false preamble candidates

experiments/lab052_candidate_detection.py

@@ -0,0 +1,287 @@
+"""Lab052: локальные группы пиков и дешёвая проверка частей PRBS7.
+
+Кандидат никогда не становится пакетом без прежних PHY/FEC/CRC. Параметры
+заданы длинами PHY и моделью шума, без ожидаемых стартов и содержимого TX.
+"""
+
+from __future__ import annotations
+
+from collections import Counter
+from contextlib import contextmanager
+from dataclasses import dataclass
+import math
+from pathlib import Path
+import sys
+import time
+
+if __package__ in (None, ""):
+    sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+
+import numpy as np
+from scipy.signal import fftconvolve
+
+from experiments import lab051_streaming_receiver as previous
+
+phy = previous.phy
+offline = previous.offline
+SAMPLES_PER_SYMBOL = 16
+PREAMBLE_SAMPLES = phy.lab045.PREAMBLE_BIT_COUNT * SAMPLES_PER_SYMBOL
+COARSE_ERROR_SAMPLES = 2000 + 8  # исходное уточнение ±2000 плюс сетка децимации
+MIN_FRAME_SAMPLES = (phy.lab045.PREAMBLE_BIT_COUNT + phy.lab045.PHY_HEADER_BIT_COUNT
+                     + phy.lab045.conv.encoded_bit_count(phy.lab045.MIN_PACKET_BYTES * 8)
+                     + 2 * phy.lab045.GUARD_SYMBOL_COUNT + phy.lab045.RRC_SPAN_SYMBOLS) * SAMPLES_PER_SYMBOL
+CLUSTER_WIDTH_SAMPLES = PREAMBLE_SAMPLES
+PART_FALSE_ALARM = .01
+
+
+@dataclass(frozen=True)
+class DetectorRule:
+    """parts=0 оставляет только группировку; 1/2/3 нужны для сравнения правил."""
+
+    parts: int = 2
+    cluster_width_samples: int = CLUSTER_WIDTH_SAMPLES
+
+    def __post_init__(self):
+        if type(self.parts) is not int or not 0 <= self.parts <= 3:
+            raise ValueError("Нужно от нуля до трёх частей")
+        if (type(self.cluster_width_samples) is not int or self.cluster_width_samples < 0
+                or self.cluster_width_samples % 8
+                or self.cluster_width_samples >= MIN_FRAME_SAMPLES - 2 * COARSE_ERROR_SAMPLES):
+            raise ValueError("Окно должно различать два минимальных физических кадра")
+
+
+def peak_groups(indexes: np.ndarray, width: int):
+    """Полная ширина от первого пика, без продолжения по цепочке соседей.
+
+    На границе IQ-блока группа обрывается: возможна лишняя проверка, но ни
+    один следующий пик не теряется из-за представителя из другого блока.
+    """
+    cursor = 0
+    while cursor < len(indexes):
+        end = int(np.searchsorted(indexes, indexes[cursor] + width, side="right"))
+        yield indexes[cursor:end]
+        cursor = end
+
+
+def cluster_peak_indexes(scores: np.ndarray, width: int, threshold: float) -> np.ndarray:
+    """Первые представители ограниченных групп; равенство → правый пик."""
+    peaks = offline._all_local_peak_indexes(scores)
+    peaks = peaks[scores[peaks] >= threshold]
+    return np.asarray([group[len(group) - 1 - int(np.argmax(scores[group[::-1]]))]
+                       for group in peak_groups(peaks, width)], dtype=np.int64)
+
+
+def part_correlations(values: np.ndarray, reference: np.ndarray, count: int) -> np.ndarray:
+    """Несколько vdot на уже децимированном согласованном входе, без timing.
+
+    Модуль сохраняет неизвестную фазу/полярность и не вводит ограничения CFO
+    через согласование фаз между половинами. Каждая часть нормируется своей
+    энергией, поэтому усиление и смена мощности фона не требуют TX reference.
+    """
+    if values.size != reference.size or count < 1:
+        return np.zeros(max(count, 0))
+    boundaries = np.linspace(0, reference.size, count + 1, dtype=int)
+    output = np.zeros(count)
+    for number, (low, high) in enumerate(zip(boundaries[:-1], boundaries[1:])):
+        part = values[low:high]
+        expected = reference[low:high]
+        energy = float(np.vdot(part, part).real) * float(np.vdot(expected, expected).real)
+        output[number] = abs(np.vdot(expected, part)) / math.sqrt(max(energy, 1e-30))
+    return output
+
+
+def part_threshold(count: int) -> float:
+    """Номинальная H0-модель: независимые комплексные отсчёты по символам.
+
+    |rho|²~Beta(1,N-1), alpha=0.01 на часть. Это исходная инженерная настройка,
+    не обещание точной вероятности после RRC/выбора пиков: фон проверяется отдельно.
+    Она фиксируется до положительных опытов и не подгоняется под их offsets.
+    """
+    symbols = phy.lab045.PREAMBLE_BIT_COUNT // count
+    return math.sqrt(1 - PART_FALSE_ALARM ** (1 / (symbols - 1)))
+
+
+class CandidateFrontEnd(previous.BlockFrontEnd):
+    """Общая корреляция — согласованный фильтр преамбулы на двух выборках/символ.
+
+    Промежуточному этапу достаточно его исходных данных; полный общий RRC-
+    выход создаётся прежним строгим трактом только при прошедшем кандидате.
+    """
+
+    def coarse_scores(self, samples, reference):
+        raw = self.decimator.get(samples, self.base)
+        mean = np.mean(raw, dtype=np.complex128).astype(np.complex64)
+        needed = math.ceil(self.owned_end / 8) + reference.size + 1
+        self.detection_values = raw[:needed].copy()
+        self.detection_values -= mean
+        self.detection_values *= previous.detection_rotation(self.detection_values.size, 8)
+        if self.detection_values.size < reference.size:
+            return np.empty(0, np.float32)
+        scores = previous.lab043._normalized_sliding_correlation(self.detection_values, reference).astype(np.float32)
+        scores[~np.isfinite(scores)] = 0
+        np.clip(scores, 0, 1, out=scores)
+        return scores
+
+    def split_scores(self, reference, count, score_count):
+        """Общие корреляции частей для всех кандидатов, без per-candidate FFT.
+
+        Проверка выполняется ДО выбора представителя: сильный случайный пик
+        одной половины не может скрыть согласованный слабый кандидат рядом.
+        """
+        output = np.empty((score_count, count), np.float32)
+        boundaries = np.linspace(0, reference.size, count + 1, dtype=int)
+        energy = np.concatenate(([0.], np.cumsum(np.abs(self.detection_values) ** 2, dtype=np.float64)))
+        for number, (low, high) in enumerate(zip(boundaries[:-1], boundaries[1:])):
+            expected = reference[low:high]
+            correlation = fftconvolve(self.detection_values[low:], np.conj(expected[::-1]), mode="valid")[:score_count]
+            local_energy = energy[high:high + score_count] - energy[low:low + score_count]
+            denominator = np.sqrt(np.maximum(local_energy * float(np.vdot(expected, expected).real), 1e-30))
+            output[:, number] = np.abs(correlation) / denominator
+        return output
+
+
+class StreamingReceiver(previous.StreamingReceiver):
+    """Ограниченный IQ и кандидаты только одного блока, API получает только IQ."""
+
+    def __init__(self, symbol_rate: int = 500000, *, block_samples: int = 65536,
+                 rule: DetectorRule = DetectorRule(), observer=None):
+        super().__init__(symbol_rate, block_samples=block_samples)
+        self.rule, self.observer = rule, observer
+        self.detection_statistics = Counter()
+        self.cheap_seconds = self.cluster_seconds = 0.0
+        self.true_strict_seconds = self.false_strict_seconds = 0.0
+        self.peak_cluster_count = self.peak_eligible_count = self.candidate_state_count = 0
+        self.working_detection_arrays_peak_bytes = 0
+        with previous.previous.streaming_mode(self.mode):
+            self.raw_reference = phy.lab045.lab044._preamble_reference_samples().copy()
+            self.matched_taps = previous.lab043._matched_filter_taps().copy()
+            self.preamble_bits = phy.lab045.lab044.preamble_bits().copy()
+        for values in (self.raw_reference, self.matched_taps, self.preamble_bits):
+            values.flags.writeable = False
+
+    @contextmanager
+    def constant_templates(self):
+        """Старый адаптер очищает lru на каждом push; эти три константы локальны.
+
+        Возвращаем именно функции до выхода streaming_mode, чтобы его
+        cache_clear и соседний старый режим работали по прежним правилам.
+        """
+        targets = ((phy.lab045.lab044, "_preamble_reference_samples", lambda: self.raw_reference),
+                   (previous.lab043, "_matched_filter_taps", lambda: self.matched_taps),
+                   (phy.lab045.lab044, "preamble_bits", lambda: self.preamble_bits))
+        originals = [(owner, name, getattr(owner, name)) for owner, name, function in targets]
+        try:
+            for owner, name, function in targets:
+                setattr(owner, name, function)
+            yield
+        finally:
+            for owner, name, function in reversed(originals):
+                setattr(owner, name, function)
+
+    def _scan(self, owned_end):
+        with self.constant_templates():
+            return self._scan_with_templates(owned_end)
+
+    def _scan_with_templates(self, owned_end):
+        width = self.rule.cluster_width_samples // self.decimation
+        front = CandidateFrontEnd(self.buffer, owned_end - self.base_sample,
+                                  self.front_statistics, decimator=self.decimator, base=self.base_sample)
+        started = time.perf_counter()
+        scores = front.coarse_scores(self.buffer, self.reference)
+        self.coarse_seconds += time.perf_counter() - started
+        threshold = np.float32(phy.lab047.FIXED_HARDWARE_CANDIDATE_THRESHOLD)
+        peaks = offline._all_local_peak_indexes(scores)
+        offsets = self.base_sample + peaks * self.decimation
+        owned = (offsets >= self.scan_sample) & (offsets < owned_end)
+        owned_peaks = peaks[owned]
+        first_stage = owned_peaks[scores[owned_peaks] >= threshold]
+        self.detection_statistics["total_local_peaks"] += len(owned_peaks)
+        self.detection_statistics["first_stage_candidates"] += len(first_stage)
+        self.peak_candidate_count = max(self.peak_candidate_count, len(first_stage))
+        if self.observer is not None:
+            self.observer.scores(scores, self.base_sample, self.scan_sample, owned_end, self.decimation)
+        started = time.perf_counter()
+        parts = (front.split_scores(self.reference, self.rule.parts, len(scores))
+                 if self.rule.parts and len(scores) else np.empty((len(scores), 0)))
+        passed = (np.all(parts[first_stage] >= part_threshold(self.rule.parts), axis=1)
+                  if self.rule.parts else np.ones(len(first_stage), dtype=bool))
+        self.cheap_seconds += time.perf_counter() - started
+        self.detection_statistics["cheap_stage_checked"] += len(first_stage)
+        self.detection_statistics["cheap_stage_rejected"] += int(np.sum(~passed))
+        indexes = first_stage[passed]
+        if self.observer is not None and hasattr(self.observer, "cheap"):
+            for index, keep in zip(first_stage, passed):
+                self.observer.cheap(self.base_sample + int(index) * self.decimation,
+                                    float(scores[index]), parts[index], bool(keep))
+        started = time.perf_counter()
+        groups = list(peak_groups(indexes, width))
+        self.cluster_seconds += time.perf_counter() - started
+        self.detection_statistics["clusters"] += len(groups)
+        self.peak_cluster_count = max(self.peak_cluster_count, len(groups))
+        self.candidate_state_count = len(indexes)
+        self.peak_eligible_count = max(self.peak_eligible_count, len(indexes))
+        self.working_detection_arrays_peak_bytes = max(self.working_detection_arrays_peak_bytes,
+            front.detection_values.nbytes + scores.nbytes * 4 + parts.nbytes + (len(peaks) + len(indexes)) * 24)
+        emitted = []
+        try:
+            with previous.optimized_context(front):
+                # Внутри группы — самый сильный первым. Отказ не даёт права
+                # подавить альтернативы: пробуем следующий. Только CRC PASS
+                # разрешает прежний пропуск уже принятого физического кадра.
+                ordered = [group[np.argsort(scores[group], kind="stable")[::-1]] for group in groups]
+                attempted_groups = set()
+                for group_number, index in ((number, item) for number, group in enumerate(ordered) for item in group):
+                    absolute = self.base_sample + int(index) * self.decimation
+                    # Сетка грубого поиска может округлить следующий стык на
+                    # одну выборку назад. Не терять соседний физический кадр.
+                    if absolute < self.accepted_until - self.decimation:
+                        self.detection_statistics["inside_accepted_frame"] += 1
+                        continue
+                    self.detection_statistics["alternative_checks" if group_number in attempted_groups
+                                              else "representative_checks"] += 1
+                    attempted_groups.add(group_number)
+                    self.candidate_count += 1
+                    candidate = offline.CoarseCandidate(int(index) * self.decimation, float(scores[index]))
+                    started = time.perf_counter()
+                    outcome = offline.decode_candidate_from_memmap(
+                        self.buffer.view(np.float32), self.buffer.size, candidate, self.candidate_count,
+                        read_iq_window=lambda raw, low, high: raw[2 * low:2 * high].view(np.complex64))
+                    elapsed = time.perf_counter() - started
+                    self.strict_seconds += elapsed
+                    if self.observer is not None:
+                        self.observer.outcome(self.base_sample, candidate, outcome)
+                    if outcome.accepted and outcome.crc == phy.lab045.CRC_PASS and outcome.frame is not None:
+                        self.true_strict_seconds += elapsed
+                        self.accepted_until = self.base_sample + outcome.next_sample
+                        self.packet_count += 1
+                        emitted.append(previous.previous.ReceivedPacket(
+                            self.base_sample + outcome.refined_start_sample, self.accepted_until,
+                            self.total_samples, bytes(outcome.frame.reception.payload_bytes)))
+                    else:
+                        self.false_strict_seconds += elapsed
+                        self.rejected_reasons[outcome.reason] += 1
+        finally:
+            self.candidate_state_count = 0
+        self.blocks_processed += 1
+        self.scan_sample = owned_end
+        retain_start = max(0, self.scan_sample - self.margin)
+        trim = retain_start - self.base_sample
+        if trim > 0:
+            self.buffer = self.buffer[trim:].copy()
+            self.base_sample = retain_start
+        return emitted
+
+
+def run_functional_tests() -> None:
+    assert CLUSTER_WIDTH_SAMPLES < MIN_FRAME_SAMPLES - 2 * COARSE_ERROR_SAMPLES
+    assert .08 < part_threshold(2) < .09
+    assert StreamingReceiver().finish() == []
+
+
+def main() -> None:
+    run_functional_tests()
+    print("Lab052 detector functional checks PASS; RF TX=0")
+
+
+if __name__ == "__main__":
+    main()