Lab048: accelerate BPSK receiver to real-time

experiments/lab043_hackrf_bpsk.py

@@ -25,6 +25,7 @@ BPSK (Binary Phase Shift Keying) — двоичная фазовая манип
 from __future__ import annotations
 
 import argparse
+from functools import lru_cache
 import hashlib
 import json
 import math
@@ -41,7 +42,7 @@ if __package__ in (None, ""):
     sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
 
 import numpy as np
-from scipy.signal import fftconvolve, resample_poly
+from scipy.signal import fftconvolve, oaconvolve, resample_poly
 
 from experiments.lab043_hackrf_hardware import (
     signed_iq_bytes_to_complex64,
@@ -362,14 +363,36 @@ def _refine_frame_start(samples: np.ndarray, tx_samples: np.ndarray, coarse_star
     return low + int(np.argmax(np.abs(correlation)))
 
 
+@lru_cache(maxsize=1)
+def _matched_filter_taps() -> np.ndarray:
+    """Вернуть неизменяемые RRC-коэффициенты доказанной рабочей точки."""
+
+    taps = radio.root_raised_cosine_taps(
+        0.35,
+        SAMPLES_PER_SYMBOL,
+        10,
+    ).astype(np.float32)
+    taps.flags.writeable = False
+    return taps
+
+
+@lru_cache(maxsize=4)
+def _carrier_downconversion_vector(start: int, length: int) -> np.ndarray:
+    """Один раз вычислить неизменяемый перенос несущей для размера окна."""
+
+    indexes = start + np.arange(length)
+    rotation = np.exp(
+        -2j * np.pi * CARRIER_OFFSET_HZ * indexes / SAMPLE_RATE_HZ
+    ).astype(np.complex64)
+    rotation.flags.writeable = False
+    return rotation
+
+
 def _matched_waveform(samples: np.ndarray, start: int, tx_length: int) -> tuple[np.ndarray, int]:
     segment = samples[start : start + tx_length] - np.mean(samples)
-    indexes = start + np.arange(len(segment))
-    baseband = segment * np.exp(
-        -2j * np.pi * CARRIER_OFFSET_HZ * indexes / SAMPLE_RATE_HZ
-    )
-    taps = radio.root_raised_cosine_taps(0.35, SAMPLES_PER_SYMBOL, 10)
-    matched = fftconvolve(baseband, taps, mode="full")
+    baseband = segment * _carrier_downconversion_vector(start, len(segment))
+    taps = _matched_filter_taps()
+    matched = oaconvolve(baseband, taps, mode="full")
     nominal = 16 * SAMPLES_PER_SYMBOL + 2 * ((len(taps) - 1) // 2)
     return matched, nominal
 
@@ -400,6 +423,34 @@ def _timing_offsets(matched: np.ndarray, nominal: int, expected_symbols: np.ndar
     return preamble_offset, oracle_offset
 
 
+def _preamble_timing_offset(
+    matched: np.ndarray,
+    nominal: int,
+    expected_symbols: np.ndarray,
+) -> int:
+    """Найти то же преамбульное смещение без второй одинаковой оценки."""
+
+    expected = np.asarray(expected_symbols)
+    offsets = np.arange(-800, 801)
+    scores = np.zeros(offsets.size, dtype=np.float64)
+    first_sample = nominal + int(offsets[0])
+    for phase in range(SAMPLES_PER_SYMBOL):
+        score_indexes = np.arange(phase, offsets.size, SAMPLES_PER_SYMBOL)
+        sequence_length = expected.size + score_indexes.size - 1
+        start = first_sample + phase
+        stop = start + sequence_length * SAMPLES_PER_SYMBOL
+        sequence = matched[start:stop:SAMPLES_PER_SYMBOL]
+        if start < 0 or sequence.size != sequence_length:
+            return _timing_offsets(matched, nominal, expected)[0]
+        windows = np.lib.stride_tricks.sliding_window_view(sequence, expected.size)
+        correlation = windows @ np.conj(expected)
+        energy = np.sum(np.abs(windows) ** 2, axis=1)
+        scores[score_indexes] = np.abs(correlation) / np.sqrt(
+            expected.size * energy + 1e-30
+        )
+    return int(offsets[int(np.argmax(scores))])
+
+
 def _estimate_cfo_gain(
     symbols: np.ndarray,
     expected_symbols: np.ndarray,