added s11 collections

This commit is contained in:
Ayzen
2026-03-26 15:48:16 +03:00
parent 9581730e41
commit 24f7ebb2fb
22 changed files with 184 additions and 859 deletions
-653
View File
@@ -1,653 +0,0 @@
"""
MIMO GPR — локализация через пересечение эллипсов
==================================================
Физика в двух словах:
Пик A-скана пары (Tx_i, Rx_j) на задержке τ означает:
|Tx → объект| + |объект → Rx| = v · τ
Это уравнение эллипса. Истинный отражатель лежит на
пересечении всех 16 эллипсов (по одному на пару).
Алгоритм:
1. S(f) → IFFT → 16 A-сканов
2. Поиск пиков: SNR = пик / медиана > порог
3. Для каждого пика → мягкий эллипс в аккумуляторе
(с компенсацией геометрического и углового затухания)
4. CLEAN: найти максимум → убрать его эллипсы → повторить
О параметре SHELL_SIGMA:
Аккумулятор — это «мягкое голосование». Каждый эллипс добавляет
не единицу, а гауссово-взвешенный вклад:
w = exp(−δ²/2σ²), где δ = |R_Tx + R_Rx v·τ|
SHELL_SIGMA — ширина этой гауссовой оболочки.
Слишком широко → ghost-цели не подавляются.
Слишком узко → вклад падает до нуля из-за дискретности сетки.
Оптимум: ~ 0.4 × δZ, где δZ = v/(2B) — разрешение по глубине.
О score:
score = количество пар (из 16), чей эллипс проходит
через данную точку с невязкой δ < 3σ.
Принимает целые значения от 0 до 16.
Максимальный score у истинного объекта = 16 (все пары согласны).
Ghost-цели имеют меньший score, т.к. согласуются только
с частью пар.
О компенсации затухания:
При генерации S(f) сигнал ослаблен:
geo(i,j) = 1/(R_Tx · R_Rx) — геометрическое ослабление
pat(i,j) = cos²(θ_Tx)·cos²(θ_Rx) — диаграмма направленности
Без компенсации глубокий/угловой отражатель будет недооценён.
Компенсация: делим вес каждого пика на ожидаемое затухание
в точке z_apparent, вычисленное для данной пары антенн.
Геометрия: плоскость XZ (X — вдоль антенн, Z — глубина).
"""
"""
MIMO GPR — локализация через пересечение эллипсов
==================================================
Версия для реальных данных
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.ndimage import gaussian_filter, label
from pathlib import Path
### Изменяемые параметры
INPUT_IDX = [0,1,2,3]
OUTPUT_IDX = [0,3]
MIN_DEPTH = 2.0 # [м] пропустить прямую волну
MAX_DEPTH = 14.0
COMP_POWER = 0.2 # степень компенсации затухания
F_START = 30*1e8 # Нижняя частота
F_STOP = 6*1e9 # Верхняя частота
# Вычитание среднего фона (background removal)
BG_SUBTRACT = True #True/False
BG_PATH = Path('/Users/ivan_root/Downloads/Telegram_dwnld/03-13-measure/2_cylinder_dif_side_and_mushrooms/preprocessed')
DATA_PATH = Path('/Users/ivan_root/Downloads/Telegram_dwnld/03-13-measure/2_cylinder_dif_side_and_mushrooms/preprocessed/0006_id1_ns1875848749103')
###Конфиги
MODE = 'point' # Один из 2х режимов `point` or `extended`
eps_r = 1.0 # <-- Диэлектрическая проницаемость среды
v = 3e8 / np.sqrt(eps_r) # скорость света в среде
# КООРДИНАТЫ АНТЕНН
# Координаты вдоль оси X на поверхности (z = 0), в метрах
# Физические координаты антенн по их реальным индексам
TX_POSITIONS = {0: 90.5*0.01, # Tx с индексом 0 → x = +0.9 м
3: -90.5*0.01} # Tx с индексом 3 → x = -0.9 м
RX_POSITIONS = {0: -18*0.01,
1: 48.5*0.01,
2: -49.0*0.01,
3: 18.5*0.01}
# Массивы позиций в порядке возрастания физических индексов
x_tx = np.array([TX_POSITIONS[i] for i in sorted(TX_POSITIONS)])
x_rx = np.array([RX_POSITIONS[j] for j in sorted(RX_POSITIONS)])
# Параметры алгоритма
SNR_THRESH = 3.0 # минимальный SNR пика
SNR_COMP_MAX = 20.0 # Верхний порог для компенсированного значения SNR
# Параметр: максимальное число объектов для поиска
MAX_OBJECTS = 15 # <-- настройте под вашу задачу
# Границы сетки аккумулятора
x_min, x_max = x_tx.min() - 2.0, x_tx.max() + 2.0 # [м]
z_min, z_max = 0.10, MAX_DEPTH # <-- глубина [м]
# ══════════════════════════════════════════════════════
# 1. ЗАГРУЗКА РЕАЛЬНЫХ ДАННЫХ
# ══════════════════════════════════════════════════════
def load_mimo_data(data_path, input_idx, output_idx):
data_path = Path(data_path)
s21_data = {}
freq_data = {}
for f in data_path.glob("i*_o*_s21.npy"):
name = f.stem
parts = name.split('_')
i_tx_phys = int(parts[1][1:])
i_rx_phys = int(parts[0][1:])
if i_tx_phys not in output_idx or i_rx_phys not in input_idx:
continue
# Переводим физический индекс → порядковый (0,1,2,...)
i_tx = sorted(output_idx).index(i_tx_phys)
i_rx = sorted(input_idx).index(i_rx_phys)
s21_data[(i_tx, i_rx)] = np.load(f)
freq_file = data_path / f"i{i_rx_phys}_o{i_tx_phys}_freq.npy"
freq_data[(i_tx, i_rx)] = np.load(freq_file)
tx_indices = sorted(set(k[0] for k in s21_data))
rx_indices = sorted(set(k[1] for k in s21_data))
n_tx = len(tx_indices)
n_rx = len(rx_indices)
return s21_data, freq_data, n_tx, n_rx
# Загрузка данных
s21_data, freq_data, N_tx, N_rx = load_mimo_data(DATA_PATH, INPUT_IDX, OUTPUT_IDX)
N_pairs = len(s21_data)
# ══════════════════════════════════════════════════════
# 1б. ВЫЧИСЛЕНИЕ СРЕДНЕГО ФОНА ПО ВСЕМ СНИМКАМ
# ══════════════════════════════════════════════════════
def compute_background(bg_path, input_idx, output_idx):
"""
Для каждой пары (i_tx, i_rx) усредняем S21 по всем снимкам в папке.
Возвращает:
bg : dict[(i_tx, i_rx)] → np.array (complex), усреднённый S21
"""
bg_path = Path(bg_path)
snapshots = sorted(bg_path.glob("*/")) # каждый подкаталог — один снимок
snapshots = [s for s in snapshots if s.is_dir()]
if len(snapshots) == 0:
print("⚠️ Снимков для фона не найдено, BG_SUBTRACT отключён.")
return None
# Накопитель: для каждой пары суммируем S21
bg_sum = {}
bg_count = {}
for snap_dir in snapshots:
for f in snap_dir.glob("i*_o*_s21.npy"):
name = f.stem
parts = name.split('_')
i_tx_phys = int(parts[1][1:])
i_rx_phys = int(parts[0][1:])
if i_tx_phys not in output_idx or i_rx_phys not in input_idx:
continue
i_tx = sorted(output_idx).index(i_tx_phys)
i_rx = sorted(input_idx).index(i_rx_phys)
key = (i_tx, i_rx)
s21 = np.load(f)
if key not in bg_sum:
bg_sum[key] = np.zeros_like(s21, dtype=complex)
bg_count[key] = 0
bg_sum[key] += s21
bg_count[key] += 1
bg = {key: bg_sum[key] / bg_count[key] for key in bg_sum}
# print(f"готово. Пар: {len(bg)}, снимков на пару: "
# f"{list(bg_count.values())[0] if bg_count else 0}")
return bg
if BG_SUBTRACT:
background = compute_background(BG_PATH, INPUT_IDX, OUTPUT_IDX)
if background is None:
BG_SUBTRACT = False # автоматически выключаем если нет данных
else:
background = None
# print("BG_SUBTRACT = False, вычитание фона отключено.")
# Проверка частот (берём из первой пары как референс)
first_key = list(freq_data.keys())[0]
freqs = freq_data[first_key]
# Проверим что частоты одинаковые для всех пар
for key, freq in freq_data.items():
if not np.allclose(freq, freqs):
print(f"⚠️ Частоты для пары {key} отличаются!")
mask_freq = (freqs >= F_START) & (freqs <= F_STOP)
freqs = freqs[mask_freq]
f_min, f_max = freqs[0], freqs[-1]
BW = f_max - f_min
N_f = len(freqs)
# ══════════════════════════════════════════════════════
# 2. ПАРАМЕТРЫ СИСТЕМЫ
# ══════════════════════════════════════════════════════
# Проверка соответствия координатов антенн
assert len(x_tx) == N_tx, f"x_tx должен содержать {N_tx} элементов"
assert len(x_rx) == N_rx, f"x_rx должен содержать {N_rx} элементов"
# SHELL_SIGMA — ширина гауссовой оболочки
SHELL_SIGMA = v / BW * 0.5 # [м]
# Сетка аккумулятора
x_grid = np.linspace(x_min, x_max, 300)
z_grid = np.linspace(z_min, z_max, 300)
XX, ZZ = np.meshgrid(x_grid, z_grid)
# Расстояния от сетки до каждой антенны
R_tx_grid = {i: np.sqrt((XX - x_tx[i])**2 + ZZ**2) for i in range(N_tx)}
R_rx_grid = {j: np.sqrt((XX - x_rx[j])**2 + ZZ**2) for j in range(N_rx)}
# ══════════════════════════════════════════════════════
# 3. ВЫЧИСЛЕНИЕ A-СКАНОВ ИЗ РЕАЛЬНЫХ S21
# ══════════════════════════════════════════════════════
def compute_ascan(s21, freq, f_start, f_stop, window=True):
"""
S21(f) → IFFT → A-скан с правильным частотным сдвигом.
Проблема наивного подхода (buf[:n] = s21):
IFFT считает, что спектр начинается с 0 Гц.
Реальные данные начинаются с f[0] > 0, поэтому
нулевая задержка смещается и в A-скане появляются биения.
Правильный подход — сдвиг спектра:
Шаг частотной сетки df вычисляется из данных.
Индекс первой частоты: k0 = round(f[0] / df).
Данные кладутся в H[k0 : k0+n], а не в H[0 : n].
Тогда IFFT корректно восстанавливает временной сигнал
с нулевой задержкой в t=0.
Размер FFT:
Минимум для покрытия всего диапазона [0, f[-1]]:
min_len = 2 * (k0 + n - 1)
Округляем вверх до степени двойки для скорости FFT.
"""
mask_freq_ = (freq >= f_start) & (freq <= f_stop)
freq = freq[mask_freq_]
s21 = s21[mask_freq_]
n = len(freq)
if n < 2:
raise ValueError("Слишком мало частотных точек")
# Шаг частотной сетки
df = (freq[-1] - freq[0]) / (n - 1)
if df <= 0:
raise ValueError("Частоты не возрастают")
# Индекс первой частоты в полной сетке от 0 до f[-1]
k0 = int(np.round(freq[0] / df))
# Минимальный размер FFT, округлённый до степени двойки
min_len = 2 * (k0 + n - 1)
n_fft = 1 << int(np.ceil(np.log2(min_len)))
# Временна́я ось — пересчитываем из нового n_fft
dt = 1.0 / (n_fft * df)
t_sec = np.arange(n_fft, dtype=float) * dt
# Оконная функция (подавление боковых лепестков IFFT)
s = s21 * np.hanning(n) if window else s21.copy()
# Спектр со сдвигом: данные на своём месте в частотной сетке
H = np.zeros(n_fft, dtype=np.complex128)
H[k0 : k0 + n] = s
y = np.abs(np.fft.ifft(H))
return t_sec[:y.size], y[:y.size]
A = {} # A[(i,j)] — амплитудный A-скан
T_h = {} # T_h[(i,j)] — временна́я ось для этой пары [с]
Z_h = {} # Z_h[(i,j)] — ось глубины [м]
for (i, j), s21 in s21_data.items():
s21_proc = s21.copy()
# Вычитание фона в частотной области
if BG_SUBTRACT and background is not None and (i, j) in background:
s21_proc = s21_proc - background[(i, j)]
# Примечание: вычитаем до обрезки по частоте и до окна —
# фон вычисляется из полных (необрезанных) данных,
# поэтому вычитание корректно в полном частотном диапазоне.
t_pair, a_pair = compute_ascan(s21_proc, freq_data[(i, j)],
f_start=F_START, f_stop=F_STOP)
T_h[(i, j)] = t_pair
Z_h[(i, j)] = t_pair * v / 2
A[(i, j)] = a_pair
bg_label = "с вычитанием фона" if BG_SUBTRACT else "без вычитания фона"
# Общая ось z для визуализации и поиска пиков
# (берём максимальный диапазон по всем парам)
z_h = Z_h[list(Z_h.keys())[0]] # все пары дают одинаковую ось, если freq совпадают
t_h = T_h[list(T_h.keys())[0]]
# 4. ВИЗУАЛИЗАЦИЯ A-СКАНОВ
def plot_ascans(A, z_h, n_tx, n_rx):
"""Отображение всех A-сканов"""
fig, axes = plt.subplots(n_tx, n_rx, figsize=(3*n_rx, 3*n_tx),
sharex=True, sharey=True)
if n_tx == 1:
axes = axes.reshape(1, -1)
if n_rx == 1:
axes = axes.reshape(-1, 1)
for i in range(n_tx):
for j in range(n_rx):
ax = axes[i, j]
if (i, j) in A:
ax.plot(z_h, A[(i, j)], 'b-', lw=0.8)
ax.set_title(f'Tx{i} → Rx{j}', fontsize=10)
ax.grid(True, alpha=0.3)
else:
ax.set_visible(False)
axes[-1, 0].set_xlabel('Глубина z [м]')
axes[0, 0].set_ylabel('Амплитуда')
fig.suptitle('A-сканы всех пар Tx-Rx', fontsize=12)
plt.tight_layout()
plt.show()
def plot_bscan(A, z_h, x_tx, x_rx):
"""B-скан: все A-сканы рядом, отсортированные по виртуальной позиции"""
pairs = sorted(A.keys(), key=lambda p: (x_tx[p[0]] + x_rx[p[1]]) / 2)
bscan = np.array([A[p] for p in pairs]).T
x_virt = [(x_tx[p[0]] + x_rx[p[1]]) / 2 for p in pairs]
plt.figure(figsize=(10, 6))
plt.imshow(bscan, aspect='auto', origin='lower',
extent=[min(x_virt), max(x_virt), z_h[0], z_h[-1]],
cmap='jet')
plt.colorbar(label='Амплитуда')
plt.xlabel('Виртуальная позиция X [м]')
plt.ylabel('Глубина Z [м]')
plt.title('B-скан (все пары)')
plt.show()
# ══════════════════════════════════════════════════════
# 5. ДЕТЕКТИРОВАНИЕ ПИКОВ
# ══════════════════════════════════════════════════════
def attenuation_at_depth(i_tx, i_rx, z_app):
"""
Ожидаемое ослабление geo·pattern для точки прямо под виртуальным
центром пары на глубине z_app.
Используется для компенсации: реальный SNR пика делится на это
значение, чтобы вес глубокого/углового объекта не занижался.
"""
xc = (x_tx[i_tx] + x_rx[i_rx]) / 2.0 # виртуальный центр
Rtx = np.sqrt((xc - x_tx[i_tx])**2 + z_app**2)
Rrx = np.sqrt((xc - x_rx[i_rx])**2 + z_app**2)
geo = 1.0 / (Rtx * Rrx + 1e-12)
pat = (z_app / (Rtx + 1e-12))**2 * (z_app / (Rrx + 1e-12))**2
return geo * pat + 1e-30 # +ε чтобы не делить на ноль
def find_peaks_snr(i_tx, i_rx, SNR_COMP_MAX = SNR_COMP_MAX):
"""
Поиск пиков A-скана.
Возвращает список dict:
z_app — кажущаяся глубина [м]
tau — задержка [с]
snr_raw — SNR без компенсации = пик / медиана
snr_comp— SNR с компенсацией ослабления (используется в аккумуляторе)
"""
ascan = A[(i_tx, i_rx)]
z_h_ij = Z_h[(i_tx, i_rx)]
t_h_ij = T_h[(i_tx, i_rx)]
noise = np.median(ascan)
i_min = np.searchsorted(z_h_ij, MIN_DEPTH)
i_max = np.searchsorted(z_h_ij, MAX_DEPTH)
min_dist = max(4, int(v / (2*BW) / (z_h_ij[1] - z_h_ij[0]) * 0.7))
idx, _ = find_peaks(ascan[i_min:i_max],
height=noise * SNR_THRESH,
distance=min_dist)
idx += i_min
result = []
for p in idx:
z_app = float(z_h_ij[p])
snr_raw = float(ascan[p] / noise)
atten = attenuation_at_depth(i_tx, i_rx, z_app)
atten_norm = atten / attenuation_at_depth(i_tx, i_rx, 2.0) # Референсная глубина - 2м
snr_comp = snr_raw / (atten_norm ** COMP_POWER + 1e-12)
snr_comp = min(snr_comp, SNR_COMP_MAX) # ← clipping
result.append({'z_app': z_app,
'tau': float(t_h_ij[p]),
'snr_raw': snr_raw,
'snr_comp': snr_comp})
return result
peaks = {(i, j): find_peaks_snr(i, j)
for i in range(N_tx) for j in range(N_rx)}
n_total = sum(len(v2) for v2 in peaks.values())
# ══════════════════════════════════════════════════════
# 5. АККУМУЛЯТОР
# ══════════════════════════════════════════════════════
def build_accumulator(exclude_z_ranges):
"""
Для каждого пика строим гауссову оболочку вокруг эллипса.
"""
acc = np.zeros_like(XX)
for i in range(N_tx):
for j in range(N_rx):
for pk in peaks[(i, j)]:
if any(lo <= pk['z_app'] <= hi for lo, hi in exclude_z_ranges):
continue
R_total = v * pk['tau']
residual = R_tx_grid[i] + R_rx_grid[j] - R_total
shell = np.exp(-0.5 * (residual / SHELL_SIGMA)**2)
acc += shell * pk['snr_comp'] # Это и есть метрика на карте накопления
return acc
def count_agreeing_ellipses(x_est, z_est, exclude_z_ranges):
"""
Score = количество пар, чей эллипс проходит
через точку (x_est, z_est) с невязкой δ < 3σ.
"""
count = 0
for i in range(N_tx):
for j in range(N_rx):
for pk in peaks[(i, j)]:
if any(lo <= pk['z_app'] <= hi for lo, hi in exclude_z_ranges):
continue
Rt = np.sqrt((x_est - x_tx[i])**2 + z_est**2)
Rr = np.sqrt((x_est - x_rx[j])**2 + z_est**2)
if abs(Rt + Rr - v*pk['tau']) < SHELL_SIGMA * 6:
count += 1
break # одна пара — один голос
return count
# ══════════════════════════════════════════════════════
# 6. ПОИСК ОБЪЕКТОВ
# ══════════════════════════════════════════════════════
def find_centroid(acc_s, iz, ix, rpz, rpx):
"""
Взвешенный центроид аккумулятора в окрестности (iz, ix).
"""
NZ, NX = acc_s.shape
iz0 = max(0, iz - rpz); iz1 = min(NZ, iz + rpz)
ix0 = max(0, ix - rpx); ix1 = min(NX, ix + rpx)
patch = acc_s[iz0:iz1, ix0:ix1].copy()
W = patch.sum()
if W <= 0:
return x_grid[ix], z_grid[iz]
rows = np.arange(iz0, iz1)[:, None] * np.ones(patch.shape)
cols = np.ones(patch.shape) * np.arange(ix0, ix1)[None, :]
iz_c = int(round(np.clip((rows * patch).sum() / W, 0, NZ-1)))
ix_c = int(round(np.clip((cols * patch).sum() / W, 0, NX-1)))
return x_grid[ix_c], z_grid[iz_c]
def clean_find(n_search=10, suppress_r_cm=7, thresh_frac=0.05):
"""
CLEAN-итерация для точечных объектов.
"""
dx = x_grid[1] - x_grid[0]
dz = z_grid[1] - z_grid[0]
rpx = int(suppress_r_cm / 100 / dx)
rpz = int(suppress_r_cm / 100 / dz)
excl_z = []
found = []
acc_initial = build_accumulator([])
for step in range(n_search):
acc = build_accumulator(excl_z)
acc_s = gaussian_filter(acc, sigma=3)
if acc_s.max() < thresh_frac * acc_initial.max():
break
iz, ix = np.unravel_index(acc_s.argmax(), acc_s.shape)
x_est, z_est = find_centroid(acc_s, iz, ix, rpz, rpx)
score = count_agreeing_ellipses(x_est, z_est, excl_z)
found.append({'x': x_est, 'z': z_est, 'score': score})
# Пики, соответствующие этому объекту
matched = [pk['z_app']
for i in range(N_tx) for j in range(N_rx)
for pk in peaks[(i, j)]
if not any(lo<=pk['z_app']<=hi for lo,hi in excl_z)
and abs(np.sqrt((x_est-x_tx[i])**2+z_est**2) +
np.sqrt((x_est-x_rx[j])**2+z_est**2) -
v*pk['tau']) < SHELL_SIGMA * 3] #Возможны изменения
if matched:
margin = SHELL_SIGMA * 1.0
excl_z.append((min(matched) - margin, max(matched) + margin))
return found, acc_initial
def extended_find(thresh_frac=0.75, min_area_cm2=2.0):
"""
Режим для протяжённых объектов (труба, плита и т.п.).
"""
acc = build_accumulator([])
acc_s = gaussian_filter(acc, sigma=3)
binary = acc_s > thresh_frac * acc_s.max()
dx = (x_grid[1]-x_grid[0])*100
dz = (z_grid[1]-z_grid[0])*100
min_pix = int(min_area_cm2 / (dx * dz))
labeled, n = label(binary)
regions = []
for k in range(1, n+1):
mask = labeled == k
if mask.sum() < min_pix:
continue
w = acc_s[mask]
xs = XX[mask]; zs = ZZ[mask]
xc = (xs * w).sum() / w.sum()
zc = (zs * w).sum() / w.sum()
score = count_agreeing_ellipses(xc, zc, [])
regions.append({'x': xc, 'z': zc, 'score': score,
'mask': mask, 'n_pix': mask.sum()})
return regions, acc
if MODE == 'point':
found, accum = clean_find(n_search=MAX_OBJECTS)
# print(f"найдено {len(found)} объектов.")
# for k, obj in enumerate(found):
# print(f" [{k+1}] x={obj['x']*100:+.1f} см, "
# f"z={obj['z']*100:.1f} см, "
# f"score={obj['score']}/{N_pairs}")
else:
regions, accum = extended_find()
found = regions
# print(f"найдено {len(regions)} регионов.")
# for k, r in enumerate(regions):
# print(f" [{k+1}] центр x={r['x']*100:+.1f} см, "
# f"z={r['z']*100:.1f} см, score={r['score']}/{N_pairs}")
# ══════════════════════════════════════════════════════
# 7. ГРАФИКИ
# ══════════════════════════════════════════════════════
# ─── График 2: карта накопления ────────────────────
fig, ax = plt.subplots(figsize=(12, 7))
acc_s = gaussian_filter(accum, sigma=3)
im = ax.imshow(acc_s,
extent=[x_grid[0]*100, x_grid[-1]*100,
z_grid[-1]*100, z_grid[0]*100],
aspect='auto', origin='upper', cmap='hot',
vmin=acc_s.max()*0.35, vmax=acc_s.max()*0.95)
plt.colorbar(im, ax=ax, label='Накопленный вес (SNR_comp × гауссова оболочка)')
# Позиции антенн
ax.plot(x_tx*100, np.zeros(N_tx), 'r^', ms=10, label='Tx', zorder=5)
ax.plot(x_rx*100, np.zeros(N_rx), 'bv', ms=10, label='Rx', zorder=5)
# Найденные объекты
for obj in found:
lbl = f"score={obj['score']}/{N_pairs}"
ax.plot(obj['x']*100, obj['z']*100, 'wD', ms=9, zorder=11,
markeredgecolor='black', mew=1.2)
ax.annotate(lbl, (obj['x']*100, obj['z']*100),
textcoords='offset points', xytext=(6, 4),
fontsize=8, color='white',
bbox=dict(boxstyle='round,pad=0.2', fc='black', alpha=0.5))
if MODE == 'extended':
for r in regions:
ax.contour(x_grid*100, z_grid*100, r['mask'].astype(float),
levels=[0.5], colors=['cyan'], linewidths=[1.2])
ax.plot([], [], 'wD', ms=9, markeredgecolor='k', mew=1.2, label='Найденные объекты')
ax.set_xlabel("X [см]"); ax.set_ylabel("Глубина Z [см]")
ax.set_title("Карта накопления эллипсов\n"
f"score = число пар из {N_pairs}, чей эллипс проходит через точку")
ax.set_xlim(x_grid[0]*100, x_grid[-1]*100)
ax.set_ylim(z_grid[-1]*100, z_grid[0]*100)
ax.legend(loc='lower right', fontsize=9)
ax.grid(alpha=0.25)
ax.invert_yaxis()
plt.tight_layout()
plt.show()
-1
View File
@@ -32,7 +32,6 @@ ORCH_SOURCES := \
data_acq_and_processing/sweep_orchestrator/device_drivers/radar/librevna/librevna_minimal_driver_protocol.cpp \ data_acq_and_processing/sweep_orchestrator/device_drivers/radar/librevna/librevna_minimal_driver_protocol.cpp \
data_acq_and_processing/sweep_orchestrator/device_drivers/switches/h7992_minimal_driver.cpp \ data_acq_and_processing/sweep_orchestrator/device_drivers/switches/h7992_minimal_driver.cpp \
data_acq_and_processing/sweep_orchestrator/device_drivers/switches/hmc349a_minimal_driver.cpp \ data_acq_and_processing/sweep_orchestrator/device_drivers/switches/hmc349a_minimal_driver.cpp \
data_acq_and_processing/sweep_orchestrator/src/sweep_plan.cpp \
data_acq_and_processing/sweep_orchestrator/src/sweep_orchestrator.cpp \ data_acq_and_processing/sweep_orchestrator/src/sweep_orchestrator.cpp \
data_acq_and_processing/sweep_orchestrator/src/main.cpp data_acq_and_processing/sweep_orchestrator/src/main.cpp
@@ -32,8 +32,10 @@ struct ComboKeyHash {
struct SweepTraceBlock { struct SweepTraceBlock {
ComboKey combo{}; ComboKey combo{};
// Frequency axis in Hz. Must have the same size as `s21`. // Frequency axis in Hz. Must have the same size as `s11` and `s21`.
std::vector<float> frequency_hz{}; std::vector<float> frequency_hz{};
// Complex S11 samples for matching frequency points.
std::vector<Complex32> s11{};
// Complex S21 samples for matching frequency points. // Complex S21 samples for matching frequency points.
std::vector<Complex32> s21{}; std::vector<Complex32> s21{};
}; };
@@ -12,8 +12,8 @@
namespace radar::ipc { namespace radar::ipc {
namespace { namespace {
constexpr std::uint32_t kRawCollectionMagic = 0x31574152U; // RAW1 constexpr std::uint32_t kRawCollectionMagic = 0x32574152U; // RAW2
constexpr std::uint32_t kPreprocessedCollectionMagic = 0x31525050U; // PRP1 constexpr std::uint32_t kPreprocessedCollectionMagic = 0x32525050U; // PRP2
constexpr std::uint32_t kResultCollectionMagic = 0x314C5352U; // RSL1 constexpr std::uint32_t kResultCollectionMagic = 0x314C5352U; // RSL1
template <typename T> template <typename T>
@@ -102,6 +102,7 @@ class BinaryReader {
void write_trace_block(BinaryWriter& writer, const SweepTraceBlock& trace) { void write_trace_block(BinaryWriter& writer, const SweepTraceBlock& trace) {
ensure_equal_sizes(trace.frequency_hz.size(), trace.s21.size(), "Sweep trace"); ensure_equal_sizes(trace.frequency_hz.size(), trace.s21.size(), "Sweep trace");
ensure_equal_sizes(trace.frequency_hz.size(), trace.s11.size(), "Sweep trace S11");
writer.write(trace.combo.input_pos); writer.write(trace.combo.input_pos);
writer.write(trace.combo.output_pos); writer.write(trace.combo.output_pos);
@@ -111,6 +112,11 @@ void write_trace_block(BinaryWriter& writer, const SweepTraceBlock& trace) {
writer.write(frequency_hz); writer.write(frequency_hz);
} }
for (const auto& point : trace.s11) {
writer.write(point.re);
writer.write(point.im);
}
for (const auto& point : trace.s21) { for (const auto& point : trace.s21) {
writer.write(point.re); writer.write(point.re);
writer.write(point.im); writer.write(point.im);
@@ -124,12 +130,20 @@ void write_trace_block(BinaryWriter& writer, const SweepTraceBlock& trace) {
const auto point_count = reader.read<std::uint32_t>(); const auto point_count = reader.read<std::uint32_t>();
trace.frequency_hz.reserve(point_count); trace.frequency_hz.reserve(point_count);
trace.s11.reserve(point_count);
trace.s21.reserve(point_count); trace.s21.reserve(point_count);
for (std::uint32_t index = 0; index < point_count; ++index) { for (std::uint32_t index = 0; index < point_count; ++index) {
trace.frequency_hz.push_back(reader.read<float>()); trace.frequency_hz.push_back(reader.read<float>());
} }
for (std::uint32_t index = 0; index < point_count; ++index) {
trace.s11.push_back(Complex32{
.re = reader.read<float>(),
.im = reader.read<float>(),
});
}
for (std::uint32_t index = 0; index < point_count; ++index) { for (std::uint32_t index = 0; index < point_count; ++index) {
trace.s21.push_back(Complex32{ trace.s21.push_back(Complex32{
.re = reader.read<float>(), .re = reader.read<float>(),
@@ -30,6 +30,11 @@ void validate_trace_layout(const ipc::SweepTraceBlock& trace, const std::string&
trace_label + " frequency/complex vector size mismatch for combo " + combo_to_string(trace.combo) trace_label + " frequency/complex vector size mismatch for combo " + combo_to_string(trace.combo)
); );
} }
if (trace.frequency_hz.size() != trace.s11.size()) {
throw std::runtime_error(
trace_label + " frequency/S11 vector size mismatch for combo " + combo_to_string(trace.combo)
);
}
} }
} // namespace } // namespace
@@ -93,6 +98,7 @@ auto CalibrationMaster::apply(const ipc::SweepTraceBlock& measured_trace) const
ipc::SweepTraceBlock output{}; ipc::SweepTraceBlock output{};
output.combo = measured_trace.combo; output.combo = measured_trace.combo;
output.frequency_hz = measured_trace.frequency_hz; output.frequency_hz = measured_trace.frequency_hz;
output.s11 = measured_trace.s11;
output.s21 = calibrator_impl_->apply(measured_trace.s21, standard.s21); output.s21 = calibrator_impl_->apply(measured_trace.s21, standard.s21);
return output; return output;
@@ -35,6 +35,11 @@ void validate_trace_layout(const ipc::SweepTraceBlock& trace, const std::string&
trace_label + " frequency/complex vector size mismatch for combo " + combo_to_string(trace.combo) trace_label + " frequency/complex vector size mismatch for combo " + combo_to_string(trace.combo)
); );
} }
if (trace.frequency_hz.size() != trace.s11.size()) {
throw std::runtime_error(
trace_label + " frequency/S11 vector size mismatch for combo " + combo_to_string(trace.combo)
);
}
} }
} // namespace } // namespace
@@ -115,6 +120,7 @@ auto ReferenceMaster::apply(const ipc::SweepTraceBlock& calibrated_trace) const
ipc::SweepTraceBlock output{}; ipc::SweepTraceBlock output{};
output.combo = calibrated_trace.combo; output.combo = calibrated_trace.combo;
output.frequency_hz = calibrated_trace.frequency_hz; output.frequency_hz = calibrated_trace.frequency_hz;
output.s11 = calibrated_trace.s11;
output.s21.resize(calibrated_trace.s21.size()); output.s21.resize(calibrated_trace.s21.size());
static_assert(sizeof(ipc::Complex32) == sizeof(float) * 2U, "Complex32 layout must be two contiguous floats"); static_assert(sizeof(ipc::Complex32) == sizeof(float) * 2U, "Complex32 layout must be two contiguous floats");
@@ -7,12 +7,13 @@
namespace radar::drivers { namespace radar::drivers {
/** /**
* @brief One S21 sweep acquired from the radar. * @brief One forward sweep acquired from the radar.
* *
* Both vectors must have equal size and aligned indices. * All vectors must have equal size and aligned indices.
*/ */
struct SweepTrace { struct SweepTrace {
std::vector<float> frequency_hz{}; std::vector<float> frequency_hz{};
std::vector<ipc::Complex32> s11{};
std::vector<ipc::Complex32> s21{}; std::vector<ipc::Complex32> s21{};
}; };
@@ -20,7 +21,8 @@ struct SweepTrace {
* @brief Minimal radar interface used by the sweep orchestrator. * @brief Minimal radar interface used by the sweep orchestrator.
* *
* Implementations are expected to be lightweight: configuration is handled by * Implementations are expected to be lightweight: configuration is handled by
* the Python layer, while this interface only opens/closes and acquires S21. * the Python layer, while this interface only opens/closes and acquires one
* forward sweep.
*/ */
class RadarDriver { class RadarDriver {
public: public:
@@ -30,8 +32,8 @@ class RadarDriver {
virtual void open() = 0; virtual void open() = 0;
/** @brief Release all allocated resources. */ /** @brief Release all allocated resources. */
virtual void close() = 0; virtual void close() = 0;
/** @brief Acquire one S21 sweep. */ /** @brief Acquire one sweep containing the forward traces exposed by the driver. */
[[nodiscard]] virtual auto acquire_s21_sweep() -> SweepTrace = 0; [[nodiscard]] virtual auto acquire_sweep() -> SweepTrace = 0;
}; };
} // namespace radar::drivers } // namespace radar::drivers
@@ -68,7 +68,7 @@ void LibreVnaMinimalDriver::close() {
is_open_ = false; is_open_ = false;
} }
auto LibreVnaMinimalDriver::acquire_s21_sweep() -> SweepTrace { auto LibreVnaMinimalDriver::acquire_sweep() -> SweepTrace {
if (!is_open_) { if (!is_open_) {
throw std::runtime_error("Radar driver is not open"); throw std::runtime_error("Radar driver is not open");
} }
@@ -111,6 +111,7 @@ auto LibreVnaMinimalDriver::acquire_s21_sweep() -> SweepTrace {
auto LibreVnaMinimalDriver::acquire_mock() -> SweepTrace { auto LibreVnaMinimalDriver::acquire_mock() -> SweepTrace {
SweepTrace trace{}; SweepTrace trace{};
trace.frequency_hz.reserve(settings_.sweep.points); trace.frequency_hz.reserve(settings_.sweep.points);
trace.s11.reserve(settings_.sweep.points);
trace.s21.reserve(settings_.sweep.points); trace.s21.reserve(settings_.sweep.points);
const auto span_hz = settings_.sweep.stop_hz - settings_.sweep.start_hz; const auto span_hz = settings_.sweep.stop_hz - settings_.sweep.start_hz;
@@ -127,7 +128,14 @@ auto LibreVnaMinimalDriver::acquire_mock() -> SweepTrace {
sample.re = envelope * std::cos(phase); sample.re = envelope * std::cos(phase);
sample.im = envelope * std::sin(phase); sample.im = envelope * std::sin(phase);
const auto reflection_phase = 0.7F * phase + 0.35F;
const auto reflection_envelope = 0.15F + 0.1F * std::cos(0.25F * phase);
ipc::Complex32 reflection{};
reflection.re = reflection_envelope * std::cos(reflection_phase);
reflection.im = reflection_envelope * std::sin(reflection_phase);
trace.frequency_hz.push_back(frequency_hz); trace.frequency_hz.push_back(frequency_hz);
trace.s11.push_back(reflection);
trace.s21.push_back(sample); trace.s21.push_back(sample);
} }
@@ -147,6 +155,7 @@ auto LibreVnaMinimalDriver::acquire_native() -> SweepTrace {
SweepTrace trace{}; SweepTrace trace{};
trace.frequency_hz.assign(settings_.sweep.points, 0.0F); trace.frequency_hz.assign(settings_.sweep.points, 0.0F);
trace.s11.assign(settings_.sweep.points, ipc::Complex32{});
trace.s21.assign(settings_.sweep.points, ipc::Complex32{}); trace.s21.assign(settings_.sweep.points, ipc::Complex32{});
std::vector<std::uint8_t> received(settings_.sweep.points, 0U); std::vector<std::uint8_t> received(settings_.sweep.points, 0U);
@@ -165,24 +174,22 @@ auto LibreVnaMinimalDriver::acquire_native() -> SweepTrace {
); );
} }
std::uint32_t point_number = 0; DecodedVnaDatapoint datapoint{};
float frequency_hz = 0.0F; if (!decode_vna_datapoint_traces(packet.payload, datapoint)) {
ipc::Complex32 s21{}; throw std::runtime_error("Failed to decode traces from VNADatapoint packet");
if (!decode_vna_datapoint_s21(packet.payload, point_number, frequency_hz, s21)) {
throw std::runtime_error("Failed to decode S21 from VNADatapoint packet");
} }
if (point_number >= settings_.sweep.points) { if (datapoint.point_number >= settings_.sweep.points) {
throw std::runtime_error("Received out-of-range VNADatapoint index"); throw std::runtime_error("Received out-of-range VNADatapoint index");
} }
if (received[point_number] == 0U) { if (received[datapoint.point_number] == 0U) {
received[point_number] = 1U; received[datapoint.point_number] = 1U;
++received_count; ++received_count;
} }
trace.frequency_hz[point_number] = frequency_hz; trace.frequency_hz[datapoint.point_number] = datapoint.frequency_hz;
trace.s21[point_number] = s21; trace.s11[datapoint.point_number] = datapoint.s11;
trace.s21[datapoint.point_number] = datapoint.s21;
} }
return trace; return trace;
@@ -35,11 +35,9 @@ auto LibreVnaMinimalDriver::encode_frame(
return frame; return frame;
} }
auto LibreVnaMinimalDriver::decode_vna_datapoint_s21( auto LibreVnaMinimalDriver::decode_vna_datapoint_traces(
std::span<const std::uint8_t> payload, std::span<const std::uint8_t> payload,
std::uint32_t& point_number_out, DecodedVnaDatapoint& datapoint_out
float& frequency_out,
ipc::Complex32& s21_out
) -> bool { ) -> bool {
// VNADatapoint payload layout: // VNADatapoint payload layout:
// [0..7]=freq_or_time, [8..9]=cdbm, [10..11]=point_number, // [0..7]=freq_or_time, [8..9]=cdbm, [10..11]=point_number,
@@ -58,17 +56,21 @@ auto LibreVnaMinimalDriver::decode_vna_datapoint_s21(
return false; return false;
} }
point_number_out = detail::read_u16_le(payload, 10); datapoint_out = DecodedVnaDatapoint{};
frequency_out = static_cast<float>(detail::read_u64_le(payload, 0)); datapoint_out.point_number = detail::read_u16_le(payload, 10);
datapoint_out.frequency_hz = static_cast<float>(detail::read_u64_le(payload, 0));
const auto real_offset = 12U; const auto real_offset = 12U;
const auto imag_offset = real_offset + (4U * num_values); const auto imag_offset = real_offset + (4U * num_values);
const auto flags_offset = imag_offset + (4U * num_values); const auto flags_offset = imag_offset + (4U * num_values);
std::array<std::complex<float>, 8> ref_by_stage{}; std::array<std::complex<float>, 8> ref_by_stage{};
std::array<std::complex<float>, 8> measured_by_stage{}; std::array<std::complex<float>, 8> s11_measured_by_stage{};
std::array<std::complex<float>, 8> s21_measured_by_stage{};
std::array<bool, 8> has_ref{}; std::array<bool, 8> has_ref{};
std::array<bool, 8> has_measured{}; std::array<bool, 8> has_s11_measured{};
std::array<bool, 8> has_s21_measured{};
constexpr float kReferenceMagnitudeSquaredEpsilon = 1e-12F;
for (std::size_t index = 0; index < num_values; ++index) { for (std::size_t index = 0; index < num_values; ++index) {
const auto flags = payload[flags_offset + index]; const auto flags = payload[flags_offset + index];
@@ -87,26 +89,43 @@ auto LibreVnaMinimalDriver::decode_vna_datapoint_s21(
ref_by_stage[stage] = value; ref_by_stage[stage] = value;
has_ref[stage] = true; has_ref[stage] = true;
} }
if ((flags & detail::kPort1Mask) != 0U && !is_reference) {
s11_measured_by_stage[stage] = value;
has_s11_measured[stage] = true;
}
if ((flags & detail::kPort2Mask) != 0U && !is_reference) { if ((flags & detail::kPort2Mask) != 0U && !is_reference) {
measured_by_stage[stage] = value; s21_measured_by_stage[stage] = value;
has_measured[stage] = true; has_s21_measured[stage] = true;
} }
} }
// We need one reference sample from port1 and one measured sample from port2. // For the forward sweep used in this project, port1 reference is the
// incident signal, port1 measured is reflection (S11 numerator), and
// port2 measured is transmission (S21 numerator).
for (std::size_t stage = 0; stage < ref_by_stage.size(); ++stage) { for (std::size_t stage = 0; stage < ref_by_stage.size(); ++stage) {
if (!has_ref[stage] || !has_measured[stage]) { if (!has_ref[stage]) {
continue; continue;
} }
if (std::norm(ref_by_stage[stage]) <= 0.0F) { if (std::norm(ref_by_stage[stage]) <= kReferenceMagnitudeSquaredEpsilon) {
continue; continue;
} }
const auto ratio = measured_by_stage[stage] / ref_by_stage[stage]; if (has_s11_measured[stage]) {
s21_out.re = ratio.real(); const auto s11_ratio = s11_measured_by_stage[stage] / ref_by_stage[stage];
s21_out.im = ratio.imag(); datapoint_out.s11.re = s11_ratio.real();
datapoint_out.s11.im = s11_ratio.imag();
datapoint_out.has_s11 = true;
}
if (has_s21_measured[stage]) {
const auto s21_ratio = s21_measured_by_stage[stage] / ref_by_stage[stage];
datapoint_out.s21.re = s21_ratio.real();
datapoint_out.s21.im = s21_ratio.imag();
datapoint_out.has_s21 = true;
}
if (datapoint_out.has_s21) {
return true; return true;
} }
}
return false; return false;
} }
@@ -127,4 +146,3 @@ auto LibreVnaMinimalDriver::crc32(std::span<const std::uint8_t> data) -> std::ui
} }
} // namespace radar::drivers } // namespace radar::drivers
@@ -30,7 +30,7 @@ struct LibreVnaMinimalDriverSettings {
}; };
/** /**
* @brief Minimal radar driver that acquires S21 sweeps from LibreVNA. * @brief Minimal radar driver that acquires forward sweeps from LibreVNA.
* *
* This class intentionally keeps scope narrow: open/close transport and * This class intentionally keeps scope narrow: open/close transport and
* acquire one sweep. Full device configuration is expected to be done by the * acquire one sweep. Full device configuration is expected to be done by the
@@ -42,7 +42,7 @@ class LibreVnaMinimalDriver final : public RadarDriver {
void open() override; void open() override;
void close() override; void close() override;
[[nodiscard]] auto acquire_s21_sweep() -> SweepTrace override; [[nodiscard]] auto acquire_sweep() -> SweepTrace override;
private: private:
/** /**
@@ -53,6 +53,18 @@ class LibreVnaMinimalDriver final : public RadarDriver {
std::vector<std::uint8_t> payload{}; std::vector<std::uint8_t> payload{};
}; };
/**
* @brief One decoded VNADatapoint containing complex values for one point.
*/
struct DecodedVnaDatapoint {
std::uint32_t point_number = 0;
float frequency_hz = 0.0F;
ipc::Complex32 s11{};
ipc::Complex32 s21{};
bool has_s11 = false;
bool has_s21 = false;
};
[[nodiscard]] auto acquire_mock() -> SweepTrace; [[nodiscard]] auto acquire_mock() -> SweepTrace;
[[nodiscard]] auto acquire_native() -> SweepTrace; [[nodiscard]] auto acquire_native() -> SweepTrace;
@@ -70,11 +82,9 @@ class LibreVnaMinimalDriver final : public RadarDriver {
[[nodiscard]] auto pop_packet(std::uint8_t packet_type) -> std::optional<NativePacket>; [[nodiscard]] auto pop_packet(std::uint8_t packet_type) -> std::optional<NativePacket>;
[[nodiscard]] static auto encode_frame(std::uint8_t packet_type, std::span<const std::uint8_t> payload) [[nodiscard]] static auto encode_frame(std::uint8_t packet_type, std::span<const std::uint8_t> payload)
-> std::vector<std::uint8_t>; -> std::vector<std::uint8_t>;
[[nodiscard]] static auto decode_vna_datapoint_s21( [[nodiscard]] static auto decode_vna_datapoint_traces(
std::span<const std::uint8_t> payload, std::span<const std::uint8_t> payload,
std::uint32_t& point_number_out, DecodedVnaDatapoint& datapoint_out
float& frequency_out,
ipc::Complex32& s21_out
) -> bool; ) -> bool;
[[nodiscard]] static auto crc32(std::span<const std::uint8_t> data) -> std::uint32_t; [[nodiscard]] static auto crc32(std::span<const std::uint8_t> data) -> std::uint32_t;
void validate_device_info_payload(std::span<const std::uint8_t> payload) const; void validate_device_info_payload(std::span<const std::uint8_t> payload) const;
@@ -6,7 +6,6 @@
#include "radar_driver.hpp" #include "radar_driver.hpp"
#include "run_config.hpp" #include "run_config.hpp"
#include "shm_ring.hpp" #include "shm_ring.hpp"
#include "sweep_plan.hpp"
#include "switch_driver.hpp" #include "switch_driver.hpp"
namespace radar::acq { namespace radar::acq {
@@ -45,7 +44,7 @@ class SweepOrchestrator {
private: private:
/** /**
* @brief Acquire one collection for all combinations in `plan_`. * @brief Acquire one collection for all configured combinations.
*/ */
[[nodiscard]] auto acquire_one_collection(std::uint64_t collection_id, const std::atomic<bool>& stop_requested) [[nodiscard]] auto acquire_one_collection(std::uint64_t collection_id, const std::atomic<bool>& stop_requested)
-> ipc::RawSweepCollection; -> ipc::RawSweepCollection;
@@ -56,7 +55,6 @@ class SweepOrchestrator {
drivers::SwitchDriver& output_switch_driver_; drivers::SwitchDriver& output_switch_driver_;
ipc::ShmRing& raw_ring_; ipc::ShmRing& raw_ring_;
ipc::ShmRing* raw_tap_ring_ = nullptr; ipc::ShmRing* raw_tap_ring_ = nullptr;
SweepPlan plan_{};
}; };
} // namespace radar::acq } // namespace radar::acq
@@ -1,22 +0,0 @@
#pragma once
#include <vector>
#include "run_config.hpp"
#include "shared_types.hpp"
namespace radar::acq {
/**
* @brief Pre-validated execution order of switch combinations for one collection.
*/
struct SweepPlan {
std::vector<ipc::ComboKey> ordered_combos{};
};
/**
* @brief Build and validate sweep execution plan from runtime config.
*/
[[nodiscard]] auto build_sweep_plan(const config::RunConfig& config) -> SweepPlan;
} // namespace radar::acq
@@ -21,7 +21,7 @@ void sleep_if_needed_ms(std::uint32_t delay_ms) {
} }
void validate_sweep(const drivers::SweepTrace& sweep) { void validate_sweep(const drivers::SweepTrace& sweep) {
if (sweep.frequency_hz.size() != sweep.s21.size()) { if (sweep.frequency_hz.size() != sweep.s11.size() || sweep.frequency_hz.size() != sweep.s21.size()) {
throw std::runtime_error("Radar driver returned inconsistent sweep vectors"); throw std::runtime_error("Radar driver returned inconsistent sweep vectors");
} }
} }
@@ -106,8 +106,7 @@ SweepOrchestrator::SweepOrchestrator(
input_switch_driver_(input_switch_driver), input_switch_driver_(input_switch_driver),
output_switch_driver_(output_switch_driver), output_switch_driver_(output_switch_driver),
raw_ring_(raw_ring), raw_ring_(raw_ring),
raw_tap_ring_(raw_tap_ring), raw_tap_ring_(raw_tap_ring) {}
plan_(build_sweep_plan(config)) {}
void SweepOrchestrator::run(const std::atomic<bool>& stop_requested) { void SweepOrchestrator::run(const std::atomic<bool>& stop_requested) {
DriverLifecycleGuard lifecycle_guard(radar_driver_, input_switch_driver_, output_switch_driver_); DriverLifecycleGuard lifecycle_guard(radar_driver_, input_switch_driver_, output_switch_driver_);
@@ -142,10 +141,10 @@ auto SweepOrchestrator::acquire_one_collection(
ipc::RawSweepCollection collection{}; ipc::RawSweepCollection collection{};
collection.collection_id = collection_id; collection.collection_id = collection_id;
collection.monotonic_ns = ipc::current_monotonic_ns(); collection.monotonic_ns = ipc::current_monotonic_ns();
collection.traces.reserve(plan_.ordered_combos.size()); collection.traces.reserve(config_.run_combos.size());
bool interrupted = false; bool interrupted = false;
for (const auto& combo : plan_.ordered_combos) { for (const auto& combo : config_.run_combos) {
if (should_stop(stop_requested)) { if (should_stop(stop_requested)) {
interrupted = true; interrupted = true;
break; break;
@@ -157,12 +156,13 @@ auto SweepOrchestrator::acquire_one_collection(
input_switch_driver_.switch_to(combo.input_pos); input_switch_driver_.switch_to(combo.input_pos);
sleep_if_needed_ms(config_.runtime.settling_ms); sleep_if_needed_ms(config_.runtime.settling_ms);
auto sweep = radar_driver_.acquire_s21_sweep(); auto sweep = radar_driver_.acquire_sweep();
validate_sweep(sweep); validate_sweep(sweep);
ipc::SweepTraceBlock trace{}; ipc::SweepTraceBlock trace{};
trace.combo = combo; trace.combo = combo;
trace.frequency_hz = std::move(sweep.frequency_hz); trace.frequency_hz = std::move(sweep.frequency_hz);
trace.s11 = std::move(sweep.s11);
trace.s21 = std::move(sweep.s21); trace.s21 = std::move(sweep.s21);
collection.traces.push_back(std::move(trace)); collection.traces.push_back(std::move(trace));
} }
@@ -1,19 +0,0 @@
#include "sweep_plan.hpp"
#include <stdexcept>
namespace radar::acq {
auto build_sweep_plan(const config::RunConfig& config) -> SweepPlan {
// RunConfig is already validated in common_cpp/config. Keep this function
// focused on plan construction only.
if (config.run_combos.empty()) {
throw std::runtime_error("run.combos must not be empty");
}
SweepPlan plan{};
plan.ordered_combos = config.run_combos;
return plan;
}
} // namespace radar::acq
@@ -2,7 +2,6 @@
from __future__ import annotations from __future__ import annotations
from python_app.gui.runtime.history import remove_last_aligned_histories
from python_app.hardware_full.librevna_service import LibreVnaService from python_app.hardware_full.librevna_service import LibreVnaService
from python_app.models.run_config_model import ComboModel, GprRxGeometryModel, GprTxGeometryModel, RunConfigModel from python_app.models.run_config_model import ComboModel, GprRxGeometryModel, GprTxGeometryModel, RunConfigModel
from python_app.models.run_config_validation import validate_gpr_model from python_app.models.run_config_validation import validate_gpr_model
@@ -200,6 +199,8 @@ class AppWindowConfigMixin:
self._clear_gpr_plot() self._clear_gpr_plot()
elif self._result_history: elif self._result_history:
self._draw_results(self._result_history[-1]) self._draw_results(self._result_history[-1])
else:
self._clear_trace_plots()
except Exception as exc: # noqa: BLE001 except Exception as exc: # noqa: BLE001
self._show_error(f"Failed to update live processing settings: {exc}") self._show_error(f"Failed to update live processing settings: {exc}")
@@ -218,72 +219,16 @@ class AppWindowConfigMixin:
self._processing_mode_pages.updateGeometry() self._processing_mode_pages.updateGeometry()
self._on_processing_live_settings_changed() self._on_processing_live_settings_changed()
def _on_bscan_clear_history_clicked(self) -> None:
"""Permanently clear all runtime histories, ring backlogs, and B-scan cache."""
self._apply_history_mode_deletion(mode_label="B-scan", remove_last_only=False)
def _on_bscan_remove_last_sweep_clicked(self) -> None:
"""Permanently delete the latest sweep from runtime histories and rings."""
self._apply_history_mode_deletion(mode_label="B-scan", remove_last_only=True)
def _on_gpr_clear_history_clicked(self) -> None:
"""Permanently clear all runtime histories, ring backlogs, and GPR cache."""
self._apply_history_mode_deletion(mode_label="GPR", remove_last_only=False)
def _on_gpr_remove_last_measurement_clicked(self) -> None:
"""Permanently delete the latest measurement from runtime histories and rings."""
self._apply_history_mode_deletion(mode_label="GPR", remove_last_only=True)
def _clear_history_mode_caches(self) -> None: def _clear_history_mode_caches(self) -> None:
"""Drop mode-specific cached render state.""" """Drop cached render state for pass-through, B-scan, and GPR views."""
self._bscan_history_floor_collection_id = 0
self._clear_bscan_plot_history() self._clear_bscan_plot_history()
if hasattr(self, "_bscan_plot"):
self._bscan_plot.clear()
self._clear_trace_plots()
if hasattr(self, "_gpr_plot"): if hasattr(self, "_gpr_plot"):
self._clear_gpr_plot() self._clear_gpr_plot()
def _apply_history_mode_deletion(self, *, mode_label: str, remove_last_only: bool) -> None:
"""Apply destructive history deletion via C++ processor history commands."""
resume_acquisition = self._supervisor.is_running()
history_command = "remove_last" if remove_last_only else "clear_all"
action = "last measurement removed" if remove_last_only else "history fully cleared"
dropped_results = 0
try:
if resume_acquisition:
self._stop_run()
if self._result_reader is not None:
dropped_results = self._result_reader.drop_all()
if remove_last_only:
retained_raw, retained_pre, retained_result = remove_last_aligned_histories(
list(self._raw_history),
list(self._pre_history),
list(self._result_history),
)
else:
retained_raw = []
retained_pre = []
retained_result = []
self._replace_runtime_history(
retained_raw=retained_raw,
retained_pre=retained_pre,
retained_result=retained_result,
)
self._clear_history_mode_caches()
self._write_live_processing_config(history_command=history_command, bump_history_seq=True)
if self._supervisor.is_processor_running():
self._drain_results_until_quiet(timeout_s=0.35, poll_s=0.01)
self._update_history_indicator()
self._redraw_after_history_deletion()
if resume_acquisition:
self._start_run()
self._log(f"{mode_label} {action}; dropped pending results={dropped_results}")
except Exception as exc: # noqa: BLE001
self._show_error(f"Failed to delete {mode_label} history: {exc}")
def _redraw_after_history_deletion(self) -> None: def _redraw_after_history_deletion(self) -> None:
"""Refresh plot immediately after destructive history deletion.""" """Refresh plot immediately after destructive history deletion."""
if self._processing_mode.currentText() == "bscan": if self._processing_mode.currentText() == "bscan":
@@ -536,7 +536,7 @@ class AppWindowPlotMixin:
plot.setLabel("left", "Depth", units="m") plot.setLabel("left", "Depth", units="m")
view_box = plot.getViewBox() view_box = plot.getViewBox()
view_box.invertY(True) view_box.invertY(False)
view_box.enableAutoRange(x=False, y=False) view_box.enableAutoRange(x=False, y=False)
if self._gpr_lookup_table is None: if self._gpr_lookup_table is None:
@@ -673,6 +673,7 @@ class AppWindowPlotMixin:
plot.setUpdatesEnabled(False) plot.setUpdatesEnabled(False)
try: try:
self._ensure_gpr_plot_items() self._ensure_gpr_plot_items()
plot.getViewBox().invertY(False)
self._clear_gpr_point_labels() self._clear_gpr_point_labels()
self._clear_gpr_region_labels() self._clear_gpr_region_labels()
self._clear_gpr_region_masks() self._clear_gpr_region_masks()
@@ -684,7 +685,7 @@ class AppWindowPlotMixin:
self._gpr_image_item.show() self._gpr_image_item.show()
plot.setXRange(x_min, x_max, padding=0.02) plot.setXRange(x_min, x_max, padding=0.02)
plot.setYRange(y_min, y_max, padding=0.02) plot.setYRange(min(0.0, y_min), y_max, padding=0.02)
x_tx, x_rx = self._selected_gpr_geometry() x_tx, x_rx = self._selected_gpr_geometry()
if x_tx.size > 0: if x_tx.size > 0:
@@ -6,7 +6,7 @@ from pathlib import Path
import time import time
from PyQt6.QtWidgets import QFileDialog from PyQt6.QtWidgets import QFileDialog
from python_app.gui.runtime.history import record_result_history from python_app.gui.runtime.history import record_result_history, remove_last_aligned_histories
class AppWindowSnapshotMixin: class AppWindowSnapshotMixin:
@@ -90,16 +90,27 @@ class AppWindowSnapshotMixin:
except Exception as exc: # noqa: BLE001 except Exception as exc: # noqa: BLE001
self._show_error(f"Failed to save VNA history JSON: {exc}") self._show_error(f"Failed to save VNA history JSON: {exc}")
def _remove_last_runtime_history(self) -> None:
"""Remove the newest runtime measurement from all stages and processor replay state."""
self._apply_runtime_history_deletion(remove_last_only=True)
def _clear_all_runtime_history(self) -> None: def _clear_all_runtime_history(self) -> None:
"""Clear all runtime histories, ring backlogs, and processor replay state.""" """Clear all runtime histories, ring backlogs, and processor replay state."""
self._apply_runtime_history_deletion(remove_last_only=False)
def _apply_runtime_history_deletion(self, *, remove_last_only: bool) -> None:
"""Apply destructive runtime-history deletion across readers, caches, and processor replay state."""
if self._capture_session is not None: if self._capture_session is not None:
self._show_error("Cannot clear runtime history during active capture sequence") self._show_error("Cannot modify runtime history during active capture sequence")
return return
resume_acquisition = self._supervisor.is_running() resume_acquisition = self._supervisor.is_running()
dropped_raw = 0 dropped_raw = 0
dropped_pre = 0 dropped_pre = 0
dropped_results = 0 dropped_results = 0
history_command = "remove_last" if remove_last_only else "clear_all"
action_label = "last runtime measurement removed" if remove_last_only else "runtime history fully cleared"
error_action = "remove last runtime measurement" if remove_last_only else "clear runtime history"
try: try:
if resume_acquisition: if resume_acquisition:
@@ -109,23 +120,36 @@ class AppWindowSnapshotMixin:
dropped_pre = self._pre_reader.drop_all() if self._pre_reader is not None else 0 dropped_pre = self._pre_reader.drop_all() if self._pre_reader is not None else 0
dropped_results = self._result_reader.drop_all() if self._result_reader is not None else 0 dropped_results = self._result_reader.drop_all() if self._result_reader is not None else 0
self._replace_runtime_history(retained_raw=[], retained_pre=[], retained_result=[]) if remove_last_only:
retained_raw, retained_pre, retained_result = remove_last_aligned_histories(
list(self._raw_history),
list(self._pre_history),
list(self._result_history),
)
else:
retained_raw = []
retained_pre = []
retained_result = []
self._replace_runtime_history(
retained_raw=retained_raw,
retained_pre=retained_pre,
retained_result=retained_result,
)
self._bscan_history_floor_collection_id = 0 self._bscan_history_floor_collection_id = 0
self._clear_history_mode_caches() self._clear_history_mode_caches()
# Clear processor-side replay cache so newly rendered B-scan starts clean. self._write_live_processing_config(history_command=history_command, bump_history_seq=True)
self._write_live_processing_config(history_command="clear_all", bump_history_seq=True)
if self._supervisor.is_processor_running(): if self._supervisor.is_processor_running():
self._drain_results_until_quiet(timeout_s=0.25, poll_s=0.01) self._drain_results_until_quiet(timeout_s=0.35, poll_s=0.01)
if self._result_reader is not None: if self._result_reader is not None:
dropped_results += self._result_reader.drop_all() dropped_results += self._result_reader.drop_all()
self._result_history.clear()
self._update_history_indicator() self._update_history_indicator()
self._redraw_after_history_deletion() self._redraw_after_history_deletion()
self._log( self._log(
"Runtime history fully cleared: " f"{action_label}: "
f"dropped raw={dropped_raw}, " f"dropped raw={dropped_raw}, "
f"preprocessed={dropped_pre}, " f"preprocessed={dropped_pre}, "
f"results={dropped_results}" f"results={dropped_results}"
@@ -134,7 +158,7 @@ class AppWindowSnapshotMixin:
if resume_acquisition: if resume_acquisition:
self._start_run() self._start_run()
except Exception as exc: # noqa: BLE001 except Exception as exc: # noqa: BLE001
self._show_error(f"Failed to clear runtime history: {exc}") self._show_error(f"Failed to {error_action}: {exc}")
def _browse_save_path(self) -> None: def _browse_save_path(self) -> None:
"""Open directory picker for snapshot output path.""" """Open directory picker for snapshot output path."""
@@ -18,6 +18,8 @@ def build_data_actions_group(owner) -> QGroupBox:
save_button.clicked.connect(owner._save_snapshot) save_button.clicked.connect(owner._save_snapshot)
save_vna_json_button = QPushButton("Save VNA History JSON") save_vna_json_button = QPushButton("Save VNA History JSON")
save_vna_json_button.clicked.connect(owner._save_vna_history_json) save_vna_json_button.clicked.connect(owner._save_vna_history_json)
remove_last_button = QPushButton("Remove Last Runtime Measurement")
remove_last_button.clicked.connect(owner._remove_last_runtime_history)
clear_history_button = QPushButton("Clear ALL Runtime History") clear_history_button = QPushButton("Clear ALL Runtime History")
clear_history_button.clicked.connect(owner._clear_all_runtime_history) clear_history_button.clicked.connect(owner._clear_all_runtime_history)
owner._save_count = QSpinBox() owner._save_count = QSpinBox()
@@ -35,6 +37,7 @@ def build_data_actions_group(owner) -> QGroupBox:
save_row.addWidget(save_button) save_row.addWidget(save_button)
save_row.addWidget(save_vna_json_button) save_row.addWidget(save_vna_json_button)
save_row.addWidget(remove_last_button)
save_row.addWidget(clear_history_button) save_row.addWidget(clear_history_button)
save_row.addWidget(QLabel("Last N")) save_row.addWidget(QLabel("Last N"))
save_row.addWidget(owner._save_count) save_row.addWidget(owner._save_count)
@@ -8,9 +8,7 @@ from PyQt6.QtWidgets import (
QDoubleSpinBox, QDoubleSpinBox,
QFormLayout, QFormLayout,
QGroupBox, QGroupBox,
QHBoxLayout,
QLineEdit, QLineEdit,
QPushButton,
QSizePolicy, QSizePolicy,
QSpinBox, QSpinBox,
QStackedWidget, QStackedWidget,
@@ -138,24 +136,12 @@ def build_processing_group(owner) -> QGroupBox:
owner._bscan_stop_freq_mhz.setSingleStep(10.0) owner._bscan_stop_freq_mhz.setSingleStep(10.0)
owner._bscan_stop_freq_mhz.setValue(8800.0) owner._bscan_stop_freq_mhz.setValue(8800.0)
owner._bscan_clear_history_button = QPushButton("Clear B-scan History")
owner._bscan_clear_history_button.clicked.connect(owner._on_bscan_clear_history_clicked)
owner._bscan_remove_last_button = QPushButton("Remove Last Sweep")
owner._bscan_remove_last_button.clicked.connect(owner._on_bscan_remove_last_sweep_clicked)
bscan_actions = QWidget(owner._processing_mode_pages)
bscan_actions_layout = QHBoxLayout(bscan_actions)
bscan_actions_layout.setContentsMargins(0, 0, 0, 0)
bscan_actions_layout.setSpacing(8)
bscan_actions_layout.addWidget(owner._bscan_remove_last_button)
bscan_actions_layout.addWidget(owner._bscan_clear_history_button)
bscan_form.addRow("Axis", owner._bscan_axis) bscan_form.addRow("Axis", owner._bscan_axis)
bscan_form.addRow("Cut m", owner._bscan_cut_m) bscan_form.addRow("Cut m", owner._bscan_cut_m)
bscan_form.addRow("Max depth m", owner._bscan_max_depth_m) bscan_form.addRow("Max depth m", owner._bscan_max_depth_m)
bscan_form.addRow("Gain", owner._bscan_gain) bscan_form.addRow("Gain", owner._bscan_gain)
bscan_form.addRow("Start MHz", owner._bscan_start_freq_mhz) bscan_form.addRow("Start MHz", owner._bscan_start_freq_mhz)
bscan_form.addRow("Stop MHz", owner._bscan_stop_freq_mhz) bscan_form.addRow("Stop MHz", owner._bscan_stop_freq_mhz)
bscan_form.addRow(bscan_actions)
owner._processing_mode_pages.addWidget(bscan_page) owner._processing_mode_pages.addWidget(bscan_page)
gpr_page = QWidget(owner._processing_mode_pages) gpr_page = QWidget(owner._processing_mode_pages)
@@ -206,17 +192,6 @@ def build_processing_group(owner) -> QGroupBox:
owner._gpr_background_mean_count.setRange(0, 10_000) owner._gpr_background_mean_count.setRange(0, 10_000)
owner._gpr_background_mean_count.setValue(10) owner._gpr_background_mean_count.setValue(10)
owner._gpr_clear_history_button = QPushButton("Clear GPR History")
owner._gpr_clear_history_button.clicked.connect(owner._on_gpr_clear_history_clicked)
owner._gpr_remove_last_button = QPushButton("Remove Last Measurement")
owner._gpr_remove_last_button.clicked.connect(owner._on_gpr_remove_last_measurement_clicked)
gpr_actions = QWidget(owner._processing_mode_pages)
gpr_actions_layout = QHBoxLayout(gpr_actions)
gpr_actions_layout.setContentsMargins(0, 0, 0, 0)
gpr_actions_layout.setSpacing(8)
gpr_actions_layout.addWidget(owner._gpr_remove_last_button)
gpr_actions_layout.addWidget(owner._gpr_clear_history_button)
gpr_form.addRow("Input positions", owner._gpr_input_positions_input) gpr_form.addRow("Input positions", owner._gpr_input_positions_input)
gpr_form.addRow("Output positions", owner._gpr_output_positions_input) gpr_form.addRow("Output positions", owner._gpr_output_positions_input)
gpr_form.addRow("Min depth m", owner._gpr_min_depth_m) gpr_form.addRow("Min depth m", owner._gpr_min_depth_m)
@@ -226,7 +201,6 @@ def build_processing_group(owner) -> QGroupBox:
gpr_form.addRow("Stop MHz", owner._gpr_stop_freq_mhz) gpr_form.addRow("Stop MHz", owner._gpr_stop_freq_mhz)
gpr_form.addRow(owner._gpr_background_subtract_enabled) gpr_form.addRow(owner._gpr_background_subtract_enabled)
gpr_form.addRow("Mean count", owner._gpr_background_mean_count) gpr_form.addRow("Mean count", owner._gpr_background_mean_count)
gpr_form.addRow(gpr_actions)
owner._processing_mode_pages.addWidget(gpr_page) owner._processing_mode_pages.addWidget(gpr_page)
owner._processing_mode.currentTextChanged.connect(owner._on_processing_mode_changed) owner._processing_mode.currentTextChanged.connect(owner._on_processing_mode_changed)
+5 -2
View File
@@ -14,8 +14,8 @@ from python_app.models.dataset_model import (
) )
from python_app.orchestration.shm.binary_cursor import ByteCursor from python_app.orchestration.shm.binary_cursor import ByteCursor
RAW_MAGIC = 0x31574152 RAW_MAGIC = 0x32574152
PREPROC_MAGIC = 0x31525050 PREPROC_MAGIC = 0x32525050
RESULT_MAGIC = 0x314C5352 RESULT_MAGIC = 0x314C5352
@@ -40,6 +40,9 @@ def decode_trace_collection(payload: bytes, expected_magic: int) -> SweepCollect
freq = np.frombuffer(cursor.read_bytes(freq_bytes), dtype="<f4").astype(np.float32, copy=False) freq = np.frombuffer(cursor.read_bytes(freq_bytes), dtype="<f4").astype(np.float32, copy=False)
interleaved_bytes = point_count * 8 interleaved_bytes = point_count * 8
# Runtime trace payloads now carry S11 before S21. The Python layer
# still works with S21 only for now, so consume and discard S11 here.
cursor.read_bytes(interleaved_bytes)
interleaved = np.frombuffer(cursor.read_bytes(interleaved_bytes), dtype="<f4") interleaved = np.frombuffer(cursor.read_bytes(interleaved_bytes), dtype="<f4")
s21 = (interleaved[0::2] + 1j * interleaved[1::2]).astype(np.complex64, copy=False) s21 = (interleaved[0::2] + 1j * interleaved[1::2]).astype(np.complex64, copy=False)
+17 -9
View File
@@ -8,31 +8,39 @@ import numpy as np
from python_app.models.dataset_model import ResultCollection, SweepCollection from python_app.models.dataset_model import ResultCollection, SweepCollection
RAW_MAGIC = 0x31574152 RAW_MAGIC = 0x32574152
PREPROC_MAGIC = 0x31525050 PREPROC_MAGIC = 0x32525050
RESULT_MAGIC = 0x314C5352 RESULT_MAGIC = 0x314C5352
def _write_interleaved_complex(buffer: bytearray, values: np.ndarray) -> None:
"""Append complex64 array as interleaved float32 real/imag pairs."""
interleaved = np.empty(values.size * 2, dtype="<f4")
interleaved[0::2] = values.real.astype("<f4", copy=False)
interleaved[1::2] = values.imag.astype("<f4", copy=False)
buffer.extend(interleaved.tobytes())
def serialize_trace_collection(collection: SweepCollection, magic: int) -> bytes: def serialize_trace_collection(collection: SweepCollection, magic: int) -> bytes:
"""Serialize one raw/preprocessed trace collection into ring-compatible binary format.""" """Serialize one raw/preprocessed trace collection into ring-compatible binary format."""
buffer = bytearray() buffer = bytearray()
buffer.extend( buffer.extend(struct.pack("<IQQI", magic, collection.collection_id, collection.monotonic_ns, len(collection.traces)))
struct.pack("<IQQI", magic, collection.collection_id, collection.monotonic_ns, len(collection.traces))
)
for trace in collection.traces: for trace in collection.traces:
freq = np.asarray(trace.frequency_hz, dtype=np.float32) freq = np.asarray(trace.frequency_hz, dtype=np.float32)
s21 = np.asarray(trace.s21, dtype=np.complex64) s21 = np.asarray(trace.s21, dtype=np.complex64)
if freq.size != s21.size: if freq.size != s21.size:
raise ValueError("Trace frequency and S21 sizes must match") raise ValueError("Trace frequency and S21 sizes must match")
# Python workflows still operate on S21 only. Emit a zero-filled S11
# channel so C++ trace bundles keep the same wire format as runtime
# rings while the Python layer remains unchanged.
s11 = np.zeros(freq.size, dtype=np.complex64)
buffer.extend(struct.pack("<III", trace.combo.input_pos, trace.combo.output_pos, int(freq.size))) buffer.extend(struct.pack("<III", trace.combo.input_pos, trace.combo.output_pos, int(freq.size)))
buffer.extend(freq.astype("<f4", copy=False).tobytes()) buffer.extend(freq.astype("<f4", copy=False).tobytes())
interleaved = np.empty(freq.size * 2, dtype="<f4") _write_interleaved_complex(buffer, s11)
interleaved[0::2] = s21.real.astype("<f4", copy=False) _write_interleaved_complex(buffer, s21)
interleaved[1::2] = s21.imag.astype("<f4", copy=False)
buffer.extend(interleaved.tobytes())
return bytes(buffer) return bytes(buffer)
-1
View File
@@ -46,7 +46,6 @@ class NpzStore(StoreApi):
suffix = f"i{trace.combo.input_pos}_o{trace.combo.output_pos}" suffix = f"i{trace.combo.input_pos}_o{trace.combo.output_pos}"
freq_key = f"freq_{suffix}" freq_key = f"freq_{suffix}"
s21_key = f"s21_{suffix}" s21_key = f"s21_{suffix}"
payload[freq_key] = np.asarray(trace.frequency_hz, dtype=np.float32) payload[freq_key] = np.asarray(trace.frequency_hz, dtype=np.float32)
payload[s21_key] = np.asarray(trace.s21, dtype=np.complex64) payload[s21_key] = np.asarray(trace.s21, dtype=np.complex64)
combo_records.append( combo_records.append(