try to modern fft

This commit is contained in:
awe
2026-01-30 12:38:17 +03:00
parent 508c835368
commit e84c155e25

View File

@ -92,8 +92,8 @@ def try_open_pyserial(path: str, baud: int, timeout: float):
return None
try:
ser = serial.Serial(path, baudrate=baud, timeout=timeout)
# Включаем hardware flow control для предотвращения потери данных
ser.rtscts = True
# ВРЕМЕННО ОТКЛЮЧЕН: hardware flow control для проверки
# ser.rtscts = True
# Увеличиваем буфер приема ядра до 64KB
try:
ser.set_buffer_size(rx_size=65536, tx_size=4096)
@ -310,6 +310,7 @@ class SweepReader(threading.Thread):
self._total_empty_lines: int = 0 # Пустых строк
self._max_buf_size: int = 0 # Максимальный размер буфера парсинга
self._read_errors: int = 0 # Ошибок чтения из порта
self._last_diag_time: float = 0.0 # Время последнего вывода диагностики
def _finalize_current(self, xs, ys):
if not xs:
@ -400,6 +401,18 @@ class SweepReader(threading.Thread):
"max_buf": self._max_buf_size,
}
# Периодический вывод детальной диагностики в stderr (каждые 10 секунд)
now = time.time()
if now - self._last_diag_time > 10.0:
self._last_diag_time = now
sys.stderr.write(
f"[DIAG] sweep={self._sweep_idx} n_valid={n_valid:.1f} "
f"lines={self._total_lines_received} parse_err={self._total_parse_errors} "
f"read_err={self._read_errors} max_buf={self._max_buf_size} "
f"dropped={self._dropped_sweeps}\n"
)
sys.stderr.flush()
# Кладём готовый свип (если очередь полна — выбрасываем самый старый)
try:
self._q.put_nowait((sweep, info))
@ -502,6 +515,50 @@ class SweepReader(threading.Thread):
pass
def apply_temporal_unwrap(
current_phase: np.ndarray,
prev_phase: Optional[np.ndarray],
phase_offset: Optional[np.ndarray],
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Применяет phase unwrapping по времени (между последовательными свипами).
Args:
current_phase: Текущая фаза (обёрнутая в [-π, π]) для всех бинов
prev_phase: Предыдущая фаза (обёрнутая), может быть None при первом вызове
phase_offset: Накопленные смещения для каждого бина, может быть None
Returns:
(unwrapped_phase, new_prev_phase, new_phase_offset)
unwrapped_phase - развёрнутая фаза
new_prev_phase - обновлённая предыдущая фаза (для следующего вызова)
new_phase_offset - обновлённые смещения (для следующего вызова)
"""
n_bins = current_phase.size
# Инициализация при первом вызове
if prev_phase is None:
prev_phase = np.zeros(n_bins, dtype=np.float32)
if phase_offset is None:
phase_offset = np.zeros(n_bins, dtype=np.float32)
# Вычисляем разницу между текущей и предыдущей фазой
delta = current_phase - prev_phase
# Обнаруживаем скачки больше π и корректируем offset
# Скачок вверх (от π к -π): delta < -π → добавляем +2π к offset
# Скачок вниз (от -π к π): delta > π → вычитаем -2π из offset
phase_offset = phase_offset - 2.0 * np.pi * np.round(delta / (2.0 * np.pi))
# Применяем накопленные смещения
unwrapped_phase = current_phase + phase_offset
# Сохраняем текущую фазу как предыдущую для следующего свипа
new_prev_phase = current_phase.copy()
new_phase_offset = phase_offset.copy()
return unwrapped_phase, new_prev_phase, new_phase_offset
def main():
parser = argparse.ArgumentParser(
description=(
@ -571,12 +628,12 @@ def main():
reader = SweepReader(args.port, args.baud, q, stop_event, fancy=bool(args.fancy))
reader.start()
# Графика
fig, axs = plt.subplots(2, 2, figsize=(12, 8))
(ax_line, ax_img), (ax_fft, ax_spec) = axs
# Графика (3 ряда x 2 колонки = 6 графиков)
fig, axs = plt.subplots(3, 2, figsize=(12, 12))
(ax_line, ax_img), (ax_fft, ax_spec), (ax_phase, ax_phase_wf) = axs
fig.canvas.manager.set_window_title(args.title) if hasattr(fig.canvas.manager, "set_window_title") else None
# Увеличим расстояния и оставим место справа под ползунки оси Y B-scan
fig.subplots_adjust(wspace=0.25, hspace=0.35, left=0.07, right=0.90, top=0.92, bottom=0.08)
fig.subplots_adjust(wspace=0.25, hspace=0.35, left=0.07, right=0.90, top=0.95, bottom=0.05)
# Состояние для отображения
current_sweep: Optional[np.ndarray] = None
@ -593,6 +650,11 @@ def main():
ring_fft = None # type: Optional[np.ndarray]
y_min_fft, y_max_fft = None, None
freq_shared: Optional[np.ndarray] = None
# Phase состояние
ring_phase = None # type: Optional[np.ndarray]
prev_phase_per_bin: Optional[np.ndarray] = None
phase_offset_per_bin: Optional[np.ndarray] = None
y_min_phase, y_max_phase = None, None
# Параметры контраста водопада спектров
spec_clip = _parse_spec_clip(getattr(args, "spec_clip", None))
# Ползунки управления Y для B-scan и контрастом
@ -668,6 +730,29 @@ def main():
ax_spec.tick_params(axis="x", labelbottom=False)
except Exception:
pass
# График фазы текущего свипа
phase_line_obj, = ax_phase.plot([], [], lw=1)
ax_phase.set_title("Фаза спектра", pad=1)
ax_phase.set_xlabel("Бин")
ax_phase.set_ylabel("Фаза, радианы")
# Водопад фазы
img_phase_obj = ax_phase_wf.imshow(
np.zeros((1, 1), dtype=np.float32),
aspect="auto",
interpolation="nearest",
origin="lower",
cmap=args.cmap,
)
ax_phase_wf.set_title("Водопад фазы", pad=12)
ax_phase_wf.set_xlabel("")
ax_phase_wf.set_ylabel("Бин")
# Не показываем численные значения по времени
try:
ax_phase_wf.tick_params(axis="x", labelbottom=False)
except Exception:
pass
# Слайдеры для управления осью Y B-scan (мин/макс) и контрастом
try:
ax_smin = fig.add_axes([0.92, 0.55, 0.02, 0.35])
@ -700,6 +785,7 @@ def main():
def ensure_buffer(_w: int):
nonlocal ring, width, head, x_shared, ring_fft, freq_shared, ring_time
nonlocal ring_phase, prev_phase_per_bin, phase_offset_per_bin
if ring is not None:
return
width = WF_WIDTH
@ -719,6 +805,14 @@ def main():
ax_spec.set_xlim(0, max_sweeps - 1)
ax_spec.set_ylim(0, max(1, fft_bins - 1))
freq_shared = np.arange(fft_bins, dtype=np.int32)
# Phase буферы: время по X, бин по Y
ring_phase = np.full((max_sweeps, fft_bins), np.nan, dtype=np.float32)
prev_phase_per_bin = np.zeros(fft_bins, dtype=np.float32)
phase_offset_per_bin = np.zeros(fft_bins, dtype=np.float32)
img_phase_obj.set_data(np.zeros((fft_bins, max_sweeps), dtype=np.float32))
img_phase_obj.set_extent((0, max_sweeps - 1, 0, fft_bins - 1))
ax_phase_wf.set_xlim(0, max_sweeps - 1)
ax_phase_wf.set_ylim(0, max(1, fft_bins - 1))
def _visible_levels_matplotlib(data: np.ndarray, axis) -> Optional[Tuple[float, float]]:
"""(vmin, vmax) по текущей видимой области imshow (без накопления по времени)."""
@ -754,6 +848,7 @@ def main():
def push_sweep(s: np.ndarray):
nonlocal ring, head, ring_fft, y_min_fft, y_max_fft, ring_time
nonlocal ring_phase, prev_phase_per_bin, phase_offset_per_bin, y_min_phase, y_max_phase
if s is None or s.size == 0 or ring is None:
return
# Нормализуем длину до фиксированной ширины
@ -765,13 +860,14 @@ def main():
if ring_time is not None:
ring_time[head] = time.time()
head = (head + 1) % ring.shape[0]
# FFT строка (дБ)
# FFT строка (дБ) и фаза
if ring_fft is not None:
bins = ring_fft.shape[1]
# Подготовка входа FFT_LEN, замена NaN на 0
take_fft = min(int(s.size), FFT_LEN)
if take_fft <= 0:
fft_row = np.full((bins,), np.nan, dtype=np.float32)
phase_row = np.full((bins,), np.nan, dtype=np.float32)
else:
fft_in = np.zeros((FFT_LEN,), dtype=np.float32)
seg = s[:take_fft]
@ -789,6 +885,19 @@ def main():
if fft_row.shape[0] != bins:
# rfft длиной FFT_LEN даёт bins == FFT_LEN//2+1
fft_row = fft_row[:bins]
# Расчет фазы
phase = np.angle(spec).astype(np.float32)
if phase.shape[0] > bins:
phase = phase[:bins]
# Unwrapping по частоте (внутри свипа)
phase_unwrapped_freq = np.unwrap(phase)
# Unwrapping по времени (между свипами)
phase_unwrapped_time, prev_phase_per_bin, phase_offset_per_bin = apply_temporal_unwrap(
phase_unwrapped_freq, prev_phase_per_bin, phase_offset_per_bin
)
phase_row = phase_unwrapped_time
ring_fft[(head - 1) % ring_fft.shape[0], :] = fft_row
# Экстремумы для цветовой шкалы
fr_min = np.nanmin(fft_row)
@ -799,6 +908,17 @@ def main():
if y_max_fft is None or (not np.isnan(fr_max) and fr_max > y_max_fft):
y_max_fft = float(fr_max)
# Сохраняем фазу в буфер
if ring_phase is not None:
ring_phase[(head - 1) % ring_phase.shape[0], :] = phase_row
# Экстремумы для цветовой шкалы фазы
ph_min = np.nanmin(phase_row)
ph_max = np.nanmax(phase_row)
if y_min_phase is None or (not np.isnan(ph_min) and ph_min < y_min_phase):
y_min_phase = float(ph_min)
if y_max_phase is None or (not np.isnan(ph_max) and ph_max > y_max_phase):
y_max_phase = float(ph_max)
def drain_queue():
nonlocal current_sweep, current_info
drained = 0
@ -833,6 +953,12 @@ def main():
base = ring_fft if head == 0 else np.roll(ring_fft, -head, axis=0)
return base.T # (bins, time)
def make_display_ring_phase():
if ring_phase is None:
return np.zeros((1, 1), dtype=np.float32)
base = ring_phase if head == 0 else np.roll(ring_phase, -head, axis=0)
return base.T # (bins, time)
def update(_frame):
nonlocal frames_since_ylim_update
changed = drain_queue() > 0
@ -861,7 +987,7 @@ def main():
y1 += pad
ax_line.set_ylim(y0, y1)
# Обновление спектра текущего свипа
# Обновление спектра и фазы текущего свипа
take_fft = min(int(current_sweep.size), FFT_LEN)
if take_fft > 0 and freq_shared is not None:
fft_in = np.zeros((FFT_LEN,), dtype=np.float32)
@ -880,6 +1006,18 @@ def main():
ax_fft.set_xlim(0, max(1, xs_fft.size - 1))
ax_fft.set_ylim(float(np.nanmin(fft_vals)), float(np.nanmax(fft_vals)))
# Расчет и отображение фазы текущего свипа
phase = np.angle(spec).astype(np.float32)
if phase.size > xs_fft.size:
phase = phase[: xs_fft.size]
# Unwrapping по частоте
phase_unwrapped = np.unwrap(phase)
phase_line_obj.set_data(xs_fft[: phase_unwrapped.size], phase_unwrapped)
# Авто-диапазон по Y для фазы
if np.isfinite(np.nanmin(phase_unwrapped)) and np.isfinite(np.nanmax(phase_unwrapped)):
ax_phase.set_xlim(0, max(1, xs_fft.size - 1))
ax_phase.set_ylim(float(np.nanmin(phase_unwrapped)), float(np.nanmax(phase_unwrapped)))
# Обновление водопада
if changed and ring is not None:
disp = make_display_ring()
@ -926,11 +1064,29 @@ def main():
vmax_eff = vmin_v + c * (vmax_v - vmin_v)
img_fft_obj.set_clim(vmin=vmin_v, vmax=vmax_eff)
# Обновление водопада фазы
if changed and ring_phase is not None:
disp_phase = make_display_ring_phase()
img_phase_obj.set_data(disp_phase)
# Автодиапазон для фазы
try:
mean_phase = np.nanmean(disp_phase, axis=1)
vmin_p = float(np.nanmin(mean_phase))
vmax_p = float(np.nanmax(mean_phase))
except Exception:
vmin_p = vmax_p = None
# Фолбэк к отслеживаемым минимум/максимумам
if (vmin_p is None or not np.isfinite(vmin_p)) or (vmax_p is None or not np.isfinite(vmax_p)) or vmin_p == vmax_p:
if y_min_phase is not None and y_max_phase is not None and np.isfinite(y_min_phase) and np.isfinite(y_max_phase) and y_min_phase != y_max_phase:
vmin_p, vmax_p = y_min_phase, y_max_phase
if vmin_p is not None and vmax_p is not None and vmin_p != vmax_p:
img_phase_obj.set_clim(vmin=vmin_p, vmax=vmax_p)
if changed and current_info:
status_text.set_text(_format_status_kv(current_info))
# Возвращаем обновлённые артисты
return (line_obj, img_obj, fft_line_obj, img_fft_obj, status_text)
return (line_obj, img_obj, fft_line_obj, img_fft_obj, phase_line_obj, img_phase_obj, status_text)
ani = FuncAnimation(fig, update, interval=interval_ms, blit=False)
@ -970,7 +1126,7 @@ def run_pyqtgraph(args):
pg.setConfigOptions(useOpenGL=True, antialias=False)
app = pg.mkQApp(args.title)
win = pg.GraphicsLayoutWidget(show=True, title=args.title)
win.resize(1200, 600)
win.resize(1200, 900)
# Плот последнего свипа (слева-сверху)
p_line = win.addPlot(row=0, col=0, title="Сырые данные")
@ -992,14 +1148,14 @@ def run_pyqtgraph(args):
img = pg.ImageItem()
p_img.addItem(img)
# FFT (слева-снизу)
# FFT (слева-средний ряд)
p_fft = win.addPlot(row=1, col=0, title="FFT")
p_fft.showGrid(x=True, y=True, alpha=0.3)
curve_fft = p_fft.plot(pen=pg.mkPen((255, 120, 80), width=1))
p_fft.setLabel("bottom", "Бин")
p_fft.setLabel("left", "Амплитуда, дБ")
# Водопад спектров (справа-снизу)
# Водопад спектров (справа-средний ряд)
p_spec = win.addPlot(row=1, col=1, title="B-scan (дБ)")
p_spec.invertY(True)
p_spec.showGrid(x=False, y=False)
@ -1012,9 +1168,29 @@ def run_pyqtgraph(args):
img_fft = pg.ImageItem()
p_spec.addItem(img_fft)
# График фазы (слева-снизу)
p_phase = win.addPlot(row=2, col=0, title="Фаза спектра")
p_phase.showGrid(x=True, y=True, alpha=0.3)
curve_phase = p_phase.plot(pen=pg.mkPen((120, 255, 80), width=1))
p_phase.setLabel("bottom", "Бин")
p_phase.setLabel("left", "Фаза, радианы")
# Водопад фазы (справа-снизу)
p_phase_wf = win.addPlot(row=2, col=1, title="Водопад фазы")
p_phase_wf.invertY(True)
p_phase_wf.showGrid(x=False, y=False)
p_phase_wf.setLabel("bottom", "Время, с (новое справа)")
try:
p_phase_wf.getAxis("bottom").setStyle(showValues=False)
except Exception:
pass
p_phase_wf.setLabel("left", "Бин (0 снизу)")
img_phase = pg.ImageItem()
p_phase_wf.addItem(img_phase)
# Статусная строка (внизу окна)
status = pg.LabelItem(justify="left")
win.addItem(status, row=2, col=0, colspan=2)
win.addItem(status, row=3, col=0, colspan=2)
# Состояние
ring: Optional[np.ndarray] = None
@ -1029,6 +1205,11 @@ def run_pyqtgraph(args):
ring_fft: Optional[np.ndarray] = None
freq_shared: Optional[np.ndarray] = None
y_min_fft, y_max_fft = None, None
# Phase состояние
ring_phase: Optional[np.ndarray] = None
prev_phase_per_bin: Optional[np.ndarray] = None
phase_offset_per_bin: Optional[np.ndarray] = None
y_min_phase, y_max_phase = None, None
# Параметры контраста водопада спектров (процентильная обрезка)
spec_clip = _parse_spec_clip(getattr(args, "spec_clip", None))
# Диапазон по Y: авто по умолчанию (поддерживает отрицательные значения)
@ -1044,6 +1225,7 @@ def run_pyqtgraph(args):
def ensure_buffer(_w: int):
nonlocal ring, head, width, x_shared, ring_fft, freq_shared
nonlocal ring_phase, prev_phase_per_bin, phase_offset_per_bin
if ring is not None:
return
width = WF_WIDTH
@ -1060,6 +1242,13 @@ def run_pyqtgraph(args):
p_spec.setRange(xRange=(0, max_sweeps - 1), yRange=(0, max(1, fft_bins - 1)), padding=0)
p_fft.setXRange(0, max(1, fft_bins - 1), padding=0)
freq_shared = np.arange(fft_bins, dtype=np.int32)
# Phase: время по оси X, бин по оси Y
ring_phase = np.full((max_sweeps, fft_bins), np.nan, dtype=np.float32)
prev_phase_per_bin = np.zeros(fft_bins, dtype=np.float32)
phase_offset_per_bin = np.zeros(fft_bins, dtype=np.float32)
img_phase.setImage(ring_phase.T, autoLevels=False)
p_phase_wf.setRange(xRange=(0, max_sweeps - 1), yRange=(0, max(1, fft_bins - 1)), padding=0)
p_phase.setXRange(0, max(1, fft_bins - 1), padding=0)
def _visible_levels_pyqtgraph(data: np.ndarray) -> Optional[Tuple[float, float]]:
"""(vmin, vmax) по текущей видимой области ImageItem (без накопления по времени)."""
@ -1094,6 +1283,7 @@ def run_pyqtgraph(args):
def push_sweep(s: np.ndarray):
nonlocal ring, head, ring_fft, y_min_fft, y_max_fft
nonlocal ring_phase, prev_phase_per_bin, phase_offset_per_bin, y_min_phase, y_max_phase
if s is None or s.size == 0 or ring is None:
return
w = ring.shape[1]
@ -1102,7 +1292,7 @@ def run_pyqtgraph(args):
row[:take] = s[:take]
ring[head, :] = row
head = (head + 1) % ring.shape[0]
# FFT строка (дБ)
# FFT строка (дБ) и фаза
if ring_fft is not None:
bins = ring_fft.shape[1]
take_fft = min(int(s.size), FFT_LEN)
@ -1116,8 +1306,22 @@ def run_pyqtgraph(args):
fft_row = 20.0 * np.log10(mag + 1e-9)
if fft_row.shape[0] != bins:
fft_row = fft_row[:bins]
# Расчет фазы
phase = np.angle(spec).astype(np.float32)
if phase.shape[0] > bins:
phase = phase[:bins]
# Unwrapping по частоте (внутри свипа)
phase_unwrapped_freq = np.unwrap(phase)
# Unwrapping по времени (между свипами)
phase_unwrapped_time, prev_phase_per_bin, phase_offset_per_bin = apply_temporal_unwrap(
phase_unwrapped_freq, prev_phase_per_bin, phase_offset_per_bin
)
phase_row = phase_unwrapped_time
else:
fft_row = np.full((bins,), np.nan, dtype=np.float32)
phase_row = np.full((bins,), np.nan, dtype=np.float32)
ring_fft[(head - 1) % ring_fft.shape[0], :] = fft_row
fr_min = np.nanmin(fft_row)
fr_max = np.nanmax(fft_row)
@ -1126,6 +1330,17 @@ def run_pyqtgraph(args):
if y_max_fft is None or (not np.isnan(fr_max) and fr_max > y_max_fft):
y_max_fft = float(fr_max)
# Сохраняем фазу в буфер
if ring_phase is not None:
ring_phase[(head - 1) % ring_phase.shape[0], :] = phase_row
# Экстремумы для цветовой шкалы фазы
ph_min = np.nanmin(phase_row)
ph_max = np.nanmax(phase_row)
if y_min_phase is None or (not np.isnan(ph_min) and ph_min < y_min_phase):
y_min_phase = float(ph_min)
if y_max_phase is None or (not np.isnan(ph_max) and ph_max > y_max_phase):
y_max_phase = float(ph_max)
def drain_queue():
nonlocal current_sweep, current_info
drained = 0
@ -1165,7 +1380,7 @@ def run_pyqtgraph(args):
margin = 0.05 * max(1.0, (y1 - y0))
p_line.setYRange(y0 - margin, y1 + margin, padding=0)
# Обновим спектр
# Обновим спектр и фазу
take_fft = min(int(current_sweep.size), FFT_LEN)
if take_fft > 0 and freq_shared is not None:
fft_in = np.zeros((FFT_LEN,), dtype=np.float32)
@ -1181,6 +1396,15 @@ def run_pyqtgraph(args):
curve_fft.setData(xs_fft[: fft_vals.size], fft_vals)
p_fft.setYRange(float(np.nanmin(fft_vals)), float(np.nanmax(fft_vals)), padding=0)
# Расчет и отображение фазы текущего свипа
phase = np.angle(spec).astype(np.float32)
if phase.size > xs_fft.size:
phase = phase[: xs_fft.size]
# Unwrapping по частоте
phase_unwrapped = np.unwrap(phase)
curve_phase.setData(xs_fft[: phase_unwrapped.size], phase_unwrapped)
p_phase.setYRange(float(np.nanmin(phase_unwrapped)), float(np.nanmax(phase_unwrapped)), padding=0)
if changed and ring is not None:
disp = ring if head == 0 else np.roll(ring, -head, axis=0)
disp = disp.T[:, ::-1] # (width, time with newest at left)
@ -1226,6 +1450,28 @@ def run_pyqtgraph(args):
else:
img_fft.setImage(disp_fft, autoLevels=False)
# Обновление водопада фазы
if changed and ring_phase is not None:
disp_phase = ring_phase if head == 0 else np.roll(ring_phase, -head, axis=0)
disp_phase = disp_phase.T[:, ::-1]
# Автодиапазон для фазы
levels_phase = None
try:
mean_phase = np.nanmean(disp_phase, axis=1)
vmin_p = float(np.nanmin(mean_phase))
vmax_p = float(np.nanmax(mean_phase))
if np.isfinite(vmin_p) and np.isfinite(vmax_p) and vmin_p != vmax_p:
levels_phase = (vmin_p, vmax_p)
except Exception:
levels_phase = None
# Фолбэк к отслеживаемым минимум/максимумам
if levels_phase is None and y_min_phase is not None and y_max_phase is not None and np.isfinite(y_min_phase) and np.isfinite(y_max_phase) and y_min_phase != y_max_phase:
levels_phase = (y_min_phase, y_max_phase)
if levels_phase is not None:
img_phase.setImage(disp_phase, autoLevels=False, levels=levels_phase)
else:
img_phase.setImage(disp_phase, autoLevels=False)
timer = pg.QtCore.QTimer()
timer.timeout.connect(update)
timer.start(interval_ms)