From 00323af0f0e5d841a534273d62c1292e1abef671 Mon Sep 17 00:00:00 2001 From: awe Date: Thu, 26 Feb 2026 14:00:56 +0300 Subject: [PATCH] arccos to apply --- calib_envelope.npy | Bin 3164 -> 3164 bytes rfg_adc_plotter/constants.py | 7 + rfg_adc_plotter/gui/matplotlib_backend.py | 55 +++-- rfg_adc_plotter/gui/pyqtgraph_backend.py | 54 +++-- rfg_adc_plotter/processing/fourier.py | 254 +++++++++++++++++++-- rfg_adc_plotter/state/app_state.py | 4 +- rfg_adc_plotter/state/ring_buffer.py | 57 ++++- tests/test_fourier_phase_reconstruction.py | 75 ++++++ tests/test_ring_buffer_fft_axis.py | 40 ++++ 9 files changed, 472 insertions(+), 74 deletions(-) create mode 100644 tests/test_fourier_phase_reconstruction.py create mode 100644 tests/test_ring_buffer_fft_axis.py diff --git a/calib_envelope.npy b/calib_envelope.npy index 27a2e6426b406aa112708ff03da94dadd7ad5958..1628af64c25c6d2093545574ae784bd1e490841f 100644 GIT binary patch literal 3164 zcmbW2X>3$g6vuB3f?dFhfKUn<2+)Fs0MoKZUgteXS*k^N0^x%L(g~YTS|F$(mmP6o z6NHLT4UswnSV5E!mj_k@wM8hDeW?gant~OnR3+-~ym=qg4}Os3m$TgeIrrZ4?vxFW z+u1qeW7ItLLdxXaf@ktm2KZ72h0;^fd?}$>`E&9o&U|`S{^VTA51lxpAeVf>l!hF(gl{?a<&x-&m=d)%Pl0HVcAs^ETu+Rs%R8zQ=J87~NFt6X-I?+tkm+pH6O4G&ZuR&$8@|S>!_jdJO)~qpC)WYG zLx0WzY5*EPaT-0*I>~6;aUQe|6s3~yjc=hH%bwEvQ7n1Y8JfW@Gt#ZjwC&rb(sxXy z?qYv6l@p7$6}=O%G1SDNeVCq2=-FVZqIlxAAbDyXYBjY^4LvT??}`}-SL0uW(@@6x zN$e#!2UVsr$m=uIh~^e?J@Gxb4#&uShTV;Q54!_fjE$gs3x5svW#X0CCB)C;&%)-K ze*Y9R6q;@pZkcITl+PwN7Z$RPPOpH~rt+<2-;JgUzRB9##9QdIj`$_k7n7r(qShc* zg=k;FE`m+0?;_rYkKQOD{ss=gKJumHxtDW-xB`2bxDLxb_2=mK6LSak0elQ!!xh#* zF)uCixqz+`|0tY=AE>`bo}8P7Z46`7V3- zE|RBCe~3NIw?t!XH&ykk#LtnJ9?aX2zOPfu-Su+fRm9BPVIR#*6tmMHv0|p`MVLZv z4ruh!51p2?F>5X76Wvm*%$a>;U$uw$ePU)}c;DM$Mw+}K_bXyu&ct`o&k~P=r(g_( zU_Nj@`5x}Qc_)gy8QfD->+a2d?zfbEmyzRqZe1RhGpis-UlV%3aGSB^F&8bcEKB019pdl-dUipaibMP8bKpa1GAGarg=xI0(N% zJL-SI?uGB6De*1rVfY!EbH=8eF_!%O%q@}gC3E*w;&g}u>30)E*9`rGYd-*5(zXCFh> z=;Q7s{e$G?`%8Uvk4QK`Pk9fVO;e{bu-&m;K;DGRv@5wzSRb}6mbobA%{^V6(Vjdr zbC`!}Lcci9Dl=x@2J_bVj?7!WOLItqj%Zq-ZGwM??{p9D!*2}e4A9_6O+4CowDPv( zdp3qY(ccazE$^WoEAOQtVB5b3BH=&qW2l#&cUk)hTMOvXD8;%augL2f>^16cQYSqW zd&nDGK_2_IjMf zKMhrI3Qj^LYbU6cUcxabhcaLf_jgi?=2J8u5tp#v0qh~{QEVmo=o$GtOMD6BH+BYO z4>@bJuWcViTS_dy6X~@J-hs{Z-_F{*_r!2lncZfG;SxTl+w zdj1cGR?r5L(Iuf1zqPxE)J4xDHEECvve)s?>P0I$(MSD%o=MKln?@tE+KncRUyMCSEPdtda&D>n_kJ?7 zBOrTzp;if=(p`yv$rSRtgYG3SZ^9er)gQ&VN9}Nc#4eWfbM`?Uv_XiPUDQ-B-@<7)0bjy# zI0}kWe>-74EP;=q3lyu|TCaLie}#}m-h9xxvX6ms#u~;?8P!7wF2HSKF1qz~XAbgQ zg}oS5Px-2$9$zV}MK2^~nX3!AZ+YCKH@Q>uO|N$r{umeyeSX#e&F~|P#1;)7$p_t_ zKDY2 zNVJ#->tHJs;j{2fMPDITH?)KLp8;<{2{G>xvkW~CJsqu@e-NiWwO-X%?G(s?B6rkt zER1B{c(=|AXg^x(%LMg38*g58Li$D z8Fk)@(^}NK?sXz_bWgNLOBoxX0qUU+YKf^JrW9Sox_7b%dsurt`*09jJLnzQ5B2zV zVc*WUmRPM%J#K(%I0Aj-I!Vq$=>6zQbRk;3Yz3`bH5>7@ppT$i(Yn{_N%yxHo9>Z% zDuZ&@uDw?s^{q9uLML<)f1LQ?_usOP@f^Z_mazx>cj*3yu&l4}soy3z4L69piS0Ie z2>laU-{wrwtZ@~??hV2r*AJF~iekX4D z8{UUa{hqa33&cP)omwu;|t0j2Zo!5)}wC;Kiq(cg^NyI)+Y%H;p zu!Tb;cwidDxH_6K0v?5tpzqco=)LUZyLF!5?z`wPYK?$(o|kZOulfx89yy`m*t7;>#vx zB|6XX+4g!TuD541Hs#NOG*G`j$OPH%*Om6M)`iT=aPn&X!!e5WfyW*7902v70J4XJ z0eR3)KYQtChx@Ep?XhZALxVe(6Sv1%n{2jS=JagaOId@Lm~i}0<4;6CypQ}Zh|j{V z+z(?}+uho#r&{gM>}u7Qy~JM;v(2AcvChAGz(!~IkF8qZALw`)Ujn_n$X+DTM>0L7 zVoUpnmf(D9>yG%G-iGgk=EO25nmOZ{GX}o_gWQuo=!6EJWqs{y3Eo17F*k}l)7?C# zn==Kw@-M=k=IZ~}@_L`5-jhz9+L}MfcZJ;N$$b_bLLWs3otj>6rO&pDeVnI{v&H6v P1n@vO7!aDnKF;|YX&PSZ diff --git a/rfg_adc_plotter/constants.py b/rfg_adc_plotter/constants.py index 6ba9cba..e1e8855 100644 --- a/rfg_adc_plotter/constants.py +++ b/rfg_adc_plotter/constants.py @@ -5,6 +5,13 @@ LOG_EXP = 2.0 # основание экспоненты для опции --l # считаем, что сигнал «меньше нуля» и домножаем свип на -1 DATA_INVERSION_THRESHOLD = 10.0 +# Частотная сетка рабочего свипа (положительная часть), ГГц +FREQ_MIN_GHZ = 3.323 +FREQ_MAX_GHZ = 14.323 + +# Скорость света для перевода времени пролёта в one-way depth +SPEED_OF_LIGHT_M_S = 299_792_458.0 + # Параметры IFFT-спектра (временной профиль из спектра 3.2..14.3 ГГц) # Двусторонний спектр формируется как: [нули -14.3..-3.2 | нули -3.2..+3.2 | данные +3.2..+14.3] ZEROS_LOW = 758 # нули от -14.3 до -3.2 ГГц diff --git a/rfg_adc_plotter/gui/matplotlib_backend.py b/rfg_adc_plotter/gui/matplotlib_backend.py index 110c7e6..0c94c5a 100644 --- a/rfg_adc_plotter/gui/matplotlib_backend.py +++ b/rfg_adc_plotter/gui/matplotlib_backend.py @@ -7,9 +7,7 @@ from typing import Optional, Tuple import numpy as np -from rfg_adc_plotter.constants import FFT_LEN, FREQ_SPAN_GHZ, IFFT_LEN - -_IFFT_T_MAX_NS = float((IFFT_LEN - 1) / (FREQ_SPAN_GHZ * 1e9) * 1e9) +from rfg_adc_plotter.constants import FFT_LEN, FREQ_MAX_GHZ, FREQ_MIN_GHZ, IFFT_LEN from rfg_adc_plotter.io.sweep_reader import SweepReader from rfg_adc_plotter.processing.normalizer import build_calib_envelopes from rfg_adc_plotter.state.app_state import AppState, format_status @@ -146,8 +144,8 @@ def run_matplotlib(args): # График спектра fft_line_obj, = ax_fft.plot([], [], lw=1) ax_fft.set_title("FFT", pad=1) - ax_fft.set_xlabel("Время, нс") - ax_fft.set_ylabel("Мощность, дБ") + ax_fft.set_xlabel("Глубина, м") + ax_fft.set_ylabel("Амплитуда") # Водопад сырых данных img_obj = ax_img.imshow( @@ -166,8 +164,8 @@ def run_matplotlib(args): np.zeros((1, 1), dtype=np.float32), aspect="auto", interpolation="nearest", origin="lower", cmap=args.cmap, ) - ax_spec.set_title("B-scan (дБ)", pad=12) - ax_spec.set_ylabel("Время, нс") + ax_spec.set_title("B-scan", pad=12) + ax_spec.set_ylabel("Глубина, м") try: ax_spec.tick_params(axis="x", labelbottom=False) except Exception: @@ -176,7 +174,7 @@ def run_matplotlib(args): # Слайдеры и чекбокс contrast_slider = None try: - fft_bins = ring.fft_bins + fft_bins = ring.fft_bins if ring.fft_bins > 0 else IFFT_LEN ax_smin = fig.add_axes([0.92, 0.55, 0.02, 0.35]) ax_smax = fig.add_axes([0.95, 0.55, 0.02, 0.35]) ax_sctr = fig.add_axes([0.98, 0.55, 0.02, 0.35]) @@ -440,23 +438,36 @@ def run_matplotlib(args): calib_cb = None line_mode_state = {"value": "raw"} - FREQ_MIN = 3.323 - FREQ_MAX = 14.323 + FREQ_MIN = float(FREQ_MIN_GHZ) + FREQ_MAX = float(FREQ_MAX_GHZ) + + def _fft_depth_max() -> float: + axis = ring.fft_depth_axis_m + if axis is None or axis.size == 0: + return 1.0 + try: + vmax = float(axis[-1]) + except Exception: + vmax = float(np.nanmax(axis)) + if not np.isfinite(vmax) or vmax <= 0.0: + return 1.0 + return vmax # --- Инициализация imshow при первом свипе --- def _init_imshow_extents(): w = ring.width ms = ring.max_sweeps - fb = ring.fft_bins + fb = max(1, int(ring.fft_bins)) + depth_max = _fft_depth_max() img_obj.set_data(np.zeros((w, ms), dtype=np.float32)) img_obj.set_extent((0, ms - 1, FREQ_MIN, FREQ_MAX)) ax_img.set_xlim(0, ms - 1) ax_img.set_ylim(FREQ_MIN, FREQ_MAX) img_fft_obj.set_data(np.zeros((fb, ms), dtype=np.float32)) - img_fft_obj.set_extent((0, ms - 1, 0.0, _IFFT_T_MAX_NS)) + img_fft_obj.set_extent((0, ms - 1, 0.0, depth_max)) ax_spec.set_xlim(0, ms - 1) - ax_spec.set_ylim(0.0, _IFFT_T_MAX_NS) - ax_fft.set_xlim(0.0, _IFFT_T_MAX_NS) + ax_spec.set_ylim(0.0, depth_max) + ax_fft.set_xlim(0.0, depth_max) _imshow_initialized = [False] @@ -544,13 +555,16 @@ def run_matplotlib(args): ax_line.autoscale_view(scalex=False, scaley=True) ax_line.set_ylabel("Y") - # Спектр — используем уже вычисленный в ring IFFT (временной профиль) - if ring.last_fft_vals is not None and ring.fft_time_axis is not None: + # Профиль по глубине — используем уже вычисленный в ring IFFT. + if ring.last_fft_vals is not None and ring.fft_depth_axis_m is not None: fft_vals = ring.last_fft_vals - xs_fft = ring.fft_time_axis + xs_fft = ring.fft_depth_axis_m n = min(fft_vals.size, xs_fft.size) - fft_line_obj.set_data(xs_fft[:n], fft_vals[:n]) - if np.isfinite(np.nanmin(fft_vals)) and np.isfinite(np.nanmax(fft_vals)): + if n > 0: + fft_line_obj.set_data(xs_fft[:n], fft_vals[:n]) + else: + fft_line_obj.set_data([], []) + if n > 0 and np.isfinite(np.nanmin(fft_vals)) and np.isfinite(np.nanmax(fft_vals)): ax_fft.set_xlim(0, float(xs_fft[n - 1])) ax_fft.set_ylim(float(np.nanmin(fft_vals)), float(np.nanmax(fft_vals))) @@ -572,6 +586,9 @@ def run_matplotlib(args): disp_fft = ring.get_display_ring_fft() disp_fft = ring.subtract_recent_mean_fft(disp_fft, spec_mean_sec) img_fft_obj.set_data(disp_fft) + depth_max = _fft_depth_max() + img_fft_obj.set_extent((0, ring.max_sweeps - 1, 0.0, depth_max)) + ax_spec.set_ylim(0.0, depth_max) levels = ring.compute_fft_levels(disp_fft, spec_clip) if levels is not None: try: diff --git a/rfg_adc_plotter/gui/pyqtgraph_backend.py b/rfg_adc_plotter/gui/pyqtgraph_backend.py index b2246ed..cff7d7f 100644 --- a/rfg_adc_plotter/gui/pyqtgraph_backend.py +++ b/rfg_adc_plotter/gui/pyqtgraph_backend.py @@ -8,16 +8,13 @@ from typing import Optional, Tuple import numpy as np -from rfg_adc_plotter.constants import FREQ_SPAN_GHZ, IFFT_LEN +from rfg_adc_plotter.constants import FREQ_MAX_GHZ, FREQ_MIN_GHZ from rfg_adc_plotter.io.sweep_reader import SweepReader from rfg_adc_plotter.processing.normalizer import build_calib_envelopes from rfg_adc_plotter.state.app_state import AppState, format_status from rfg_adc_plotter.state.ring_buffer import RingBuffer from rfg_adc_plotter.types import SweepPacket -# Максимальное значение временной оси IFFT в нс -_IFFT_T_MAX_NS = float((IFFT_LEN - 1) / (FREQ_SPAN_GHZ * 1e9) * 1e9) - def _parse_ylim(ylim_str: Optional[str]) -> Optional[Tuple[float, float]]: if not ylim_str: @@ -202,11 +199,11 @@ def run_pyqtgraph(args): 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_fft.setLabel("bottom", "Глубина, м") + p_fft.setLabel("left", "Амплитуда") # Водопад спектров (справа-снизу) - p_spec = win.addPlot(row=1, col=1, title="B-scan (дБ)") + p_spec = win.addPlot(row=1, col=1, title="B-scan") p_spec.invertY(True) p_spec.showGrid(x=False, y=False) p_spec.setLabel("bottom", "Время (новое справа)") @@ -214,7 +211,7 @@ def run_pyqtgraph(args): p_spec.getAxis("bottom").setStyle(showValues=False) except Exception: pass - p_spec.setLabel("left", "Время, нс") + p_spec.setLabel("left", "Глубина, м") img_fft = pg.ImageItem() p_spec.addItem(img_fft) @@ -470,20 +467,33 @@ def run_pyqtgraph(args): _imshow_initialized = [False] - FREQ_MIN = 3.323 - FREQ_MAX = 14.323 + FREQ_MIN = float(FREQ_MIN_GHZ) + FREQ_MAX = float(FREQ_MAX_GHZ) + + def _fft_depth_max() -> float: + axis = ring.fft_depth_axis_m + if axis is None or axis.size == 0: + return 1.0 + try: + vmax = float(axis[-1]) + except Exception: + vmax = float(np.nanmax(axis)) + if not np.isfinite(vmax) or vmax <= 0.0: + return 1.0 + return vmax def _init_imshow_extents(): ms = ring.max_sweeps - fb = ring.fft_bins img.setImage(ring.ring.T, autoLevels=False) img.setRect(pg.QtCore.QRectF(0.0, FREQ_MIN, float(ms), FREQ_MAX - FREQ_MIN)) p_img.setRange(xRange=(0, ms - 1), yRange=(FREQ_MIN, FREQ_MAX), padding=0) p_line.setXRange(FREQ_MIN, FREQ_MAX, padding=0) - img_fft.setImage(ring.ring_fft.T, autoLevels=False) - img_fft.setRect(pg.QtCore.QRectF(0.0, 0.0, float(ms), _IFFT_T_MAX_NS)) - p_spec.setRange(xRange=(0, ms - 1), yRange=(0.0, _IFFT_T_MAX_NS), padding=0) - p_fft.setXRange(0.0, _IFFT_T_MAX_NS, padding=0) + disp_fft = ring.get_display_ring_fft() + img_fft.setImage(disp_fft, autoLevels=False) + depth_max = _fft_depth_max() + img_fft.setRect(pg.QtCore.QRectF(0.0, 0.0, float(ms), depth_max)) + p_spec.setRange(xRange=(0, ms - 1), yRange=(0.0, depth_max), padding=0) + p_fft.setXRange(0.0, depth_max, padding=0) def _img_rect(ms: int) -> "pg.QtCore.QRectF": return pg.QtCore.QRectF(0.0, FREQ_MIN, float(ms), FREQ_MAX - FREQ_MIN) @@ -573,13 +583,15 @@ def run_pyqtgraph(args): p_line.enableAutoRange(axis="y", enable=True) p_line.setLabel("left", "Y") - # Спектр — используем уже вычисленный в ring IFFT (временной профиль) - if ring.last_fft_vals is not None and ring.fft_time_axis is not None: + # Профиль по глубине — используем уже вычисленный в ring IFFT. + if ring.last_fft_vals is not None and ring.fft_depth_axis_m is not None: fft_vals = ring.last_fft_vals - xs_fft = ring.fft_time_axis + xs_fft = ring.fft_depth_axis_m n = min(fft_vals.size, xs_fft.size) - curve_fft.setData(xs_fft[:n], fft_vals[:n]) - p_fft.setYRange(float(np.nanmin(fft_vals)), float(np.nanmax(fft_vals)), padding=0) + if n > 0: + curve_fft.setData(xs_fft[:n], fft_vals[:n]) + p_fft.setXRange(0.0, float(xs_fft[n - 1]), padding=0) + p_fft.setYRange(float(np.nanmin(fft_vals)), float(np.nanmax(fft_vals)), padding=0) # Позиция подписи канала try: @@ -619,7 +631,7 @@ def run_pyqtgraph(args): img_fft.setImage(disp_fft, autoLevels=False, levels=levels) else: img_fft.setImage(disp_fft, autoLevels=False) - img_fft.setRect(pg.QtCore.QRectF(0.0, 0.0, float(ring.max_sweeps), _IFFT_T_MAX_NS)) + img_fft.setRect(pg.QtCore.QRectF(0.0, 0.0, float(ring.max_sweeps), _fft_depth_max())) timer = pg.QtCore.QTimer() timer.timeout.connect(update) diff --git a/rfg_adc_plotter/processing/fourier.py b/rfg_adc_plotter/processing/fourier.py index 54a7f8f..65cd733 100644 --- a/rfg_adc_plotter/processing/fourier.py +++ b/rfg_adc_plotter/processing/fourier.py @@ -1,43 +1,251 @@ -"""Преобразование свипа в IFFT-временной профиль (дБ).""" +"""Преобразование свипа в IFFT-профиль по глубине (м). +Новый pipeline перед IFFT: +1) нормировка по max(abs(sweep)) +2) clip в [-1, 1] +3) phi = arccos(x) +4) непрерывная развёртка фазы (nearest-continuous) +5) s_complex = exp(1j * phi) +6) IFFT с учётом смещения частотной сетки +""" + +from __future__ import annotations + +import logging from typing import Optional import numpy as np -from rfg_adc_plotter.constants import FREQ_SPAN_GHZ, IFFT_LEN, SWEEP_LEN, ZEROS_LOW, ZEROS_MID +from rfg_adc_plotter.constants import ( + FREQ_MAX_GHZ, + FREQ_MIN_GHZ, + FREQ_SPAN_GHZ, + IFFT_LEN, + SPEED_OF_LIGHT_M_S, +) + + +logger = logging.getLogger(__name__) + +_EPS = 1e-12 +_TWO_PI = float(2.0 * np.pi) + + +def _fallback_depth_response( + size: int, + values: Optional[np.ndarray] = None, +) -> tuple[np.ndarray, np.ndarray]: + """Безопасный fallback для GUI/ring: всегда возвращает ненулевую длину.""" + n = max(1, int(size)) + depth = np.linspace(0.0, 1.0, n, dtype=np.float32) + if values is None: + return depth, np.zeros((n,), dtype=np.float32) + + arr = np.asarray(values) + if arr.size == 0: + return depth, np.zeros((n,), dtype=np.float32) + + if np.iscomplexobj(arr): + src = np.abs(arr) + else: + src = np.abs(np.nan_to_num(arr, nan=0.0, posinf=0.0, neginf=0.0)) + src = np.asarray(src, dtype=np.float32).ravel() + + out = np.zeros((n,), dtype=np.float32) + take = min(n, src.size) + if take > 0: + out[:take] = src[:take] + return depth, out def build_ifft_time_axis_ns() -> np.ndarray: - """Временная ось IFFT в наносекундах.""" + """Legacy helper: старая временная ось IFFT в наносекундах (фиксированная длина).""" return ( np.arange(IFFT_LEN, dtype=np.float64) / (FREQ_SPAN_GHZ * 1e9) * 1e9 ).astype(np.float32) -def compute_ifft_db_profile(sweep: Optional[np.ndarray]) -> np.ndarray: - """Построить IFFT-профиль свипа в дБ. +def build_frequency_axis_hz(sweep_width: int) -> np.ndarray: + """Построить частотную сетку (Гц) для текущей длины свипа.""" + n = int(sweep_width) + if n <= 0: + return np.zeros((0,), dtype=np.float64) + if n == 1: + return np.array([FREQ_MIN_GHZ * 1e9], dtype=np.float64) + return np.linspace(FREQ_MIN_GHZ * 1e9, FREQ_MAX_GHZ * 1e9, n, dtype=np.float64) - Цепочка: - raw/processed sweep -> двусторонний спектр (заполнение нулями) -> - ifftshift -> ifft -> |x| -> 20log10. + +def normalize_sweep_for_phase(sweep: np.ndarray) -> np.ndarray: + """Нормировать свип на max(abs(.)) и вернуть float64.""" + x = np.asarray(sweep, dtype=np.float64).ravel() + if x.size == 0: + return x + x = np.nan_to_num(x, nan=0.0, posinf=0.0, neginf=0.0) + amax = float(np.max(np.abs(x))) + if (not np.isfinite(amax)) or amax <= _EPS: + return np.zeros_like(x, dtype=np.float64) + return x / amax + + +def unwrap_arccos_phase_continuous(x_norm: np.ndarray) -> np.ndarray: + """Непрерывно развернуть фазу, восстановленную через arccos. + + Для каждой точки рассматриваются ветви ±phi + 2πk и выбирается кандидат, + ближайший к предыдущей фазе (nearest continuous). """ - bins = IFFT_LEN + x = np.asarray(x_norm, dtype=np.float64).ravel() + if x.size == 0: + return np.zeros((0,), dtype=np.float64) + x = np.nan_to_num(x, nan=0.0, posinf=1.0, neginf=-1.0) + x = np.clip(x, -1.0, 1.0) + phi0 = np.arccos(x) + + out = np.empty_like(phi0, dtype=np.float64) + out[0] = float(phi0[0]) + for i in range(1, phi0.size): + base_phi = float(phi0[i]) + prev = float(out[i - 1]) + + best_cand: Optional[float] = None + best_key: Optional[tuple[float, float]] = None + + for sign in (1.0, -1.0): + base = sign * base_phi + # Ищем ближайший сдвиг 2πk относительно prev именно для этой ветви. + k_center = int(np.round((prev - base) / _TWO_PI)) + for k in (k_center - 1, k_center, k_center + 1): + cand = base + _TWO_PI * float(k) + step = abs(cand - prev) + # Tie-break: при равенстве шага предпочесть больший кандидат. + key = (step, -cand) + if best_key is None or key < best_key: + best_key = key + best_cand = cand + + out[i] = prev if best_cand is None else float(best_cand) + + return out + + +def reconstruct_complex_spectrum_from_real_trace(sweep: np.ndarray) -> np.ndarray: + """Восстановить комплексный спектр из вещественного следа через arccos+Euler.""" + x_norm = normalize_sweep_for_phase(sweep) + if x_norm.size == 0: + return np.zeros((0,), dtype=np.complex128) + x_norm = np.clip(x_norm, -1.0, 1.0) + phi = unwrap_arccos_phase_continuous(x_norm) + return np.exp(1j * phi).astype(np.complex128, copy=False) + + +def perform_ifft_depth_response( + s_array: np.ndarray, + frequencies_hz: np.ndarray, + *, + axis: str = "abs", + start_hz: float | None = None, + stop_hz: float | None = None, +) -> tuple[np.ndarray, np.ndarray]: + """Frequency-to-depth conversion with zero-padding and frequency offset handling.""" + try: + s_in = np.asarray(s_array, dtype=np.complex128).ravel() + f_in = np.asarray(frequencies_hz, dtype=np.float64).ravel() + m = min(s_in.size, f_in.size) + if m < 2: + raise ValueError("Not enough points") + + s = s_in[:m] + f = f_in[:m] + + lo = float(FREQ_MIN_GHZ * 1e9 if start_hz is None else start_hz) + hi = float(FREQ_MAX_GHZ * 1e9 if stop_hz is None else stop_hz) + if hi < lo: + lo, hi = hi, lo + + mask = ( + np.isfinite(f) + & np.isfinite(np.real(s)) + & np.isfinite(np.imag(s)) + & (f >= lo) + & (f <= hi) + ) + f = f[mask] + s = s[mask] + + n = int(f.size) + if n < 2: + raise ValueError("Not enough frequency points after filtering") + + if np.any(np.diff(f) <= 0.0): + raise ValueError("Non-increasing frequency grid") + + df = float((f[-1] - f[0]) / (n - 1)) + if not np.isfinite(df) or df <= 0.0: + raise ValueError("Invalid frequency step") + + k0 = int(np.round(float(f[0]) / df)) + if k0 < 0: + raise ValueError("Negative frequency offset index") + + min_len = int(2 * (k0 + n - 1)) + if min_len <= 0: + raise ValueError("Invalid FFT length") + n_fft = 1 << int(np.ceil(np.log2(float(min_len)))) + + dt = 1.0 / (n_fft * df) + t_sec = np.arange(n_fft, dtype=np.float64) * dt + + h = np.zeros((n_fft,), dtype=np.complex128) + end = k0 + n + if end > n_fft: + raise ValueError("Spectrum placement exceeds FFT buffer") + h[k0:end] = s + + y = np.fft.ifft(h) + depth_m = t_sec * SPEED_OF_LIGHT_M_S + + axis_name = str(axis).strip().lower() + if axis_name == "abs": + y_fin = np.abs(y) + elif axis_name == "real": + y_fin = np.real(y) + elif axis_name == "imag": + y_fin = np.imag(y) + elif axis_name == "phase": + y_fin = np.angle(y) + else: + raise ValueError(f"Invalid axis parameter: {axis!r}") + + return depth_m.astype(np.float32, copy=False), np.asarray(y_fin, dtype=np.float32) + + except Exception as exc: # noqa: BLE001 + logger.error("IFFT depth response failed: %r", exc) + return _fallback_depth_response(np.asarray(s_array).size, np.asarray(s_array)) + + +def compute_ifft_profile_from_sweep(sweep: Optional[np.ndarray]) -> tuple[np.ndarray, np.ndarray]: + """Высокоуровневый pipeline: sweep -> complex spectrum -> IFFT(abs) depth profile.""" if sweep is None: - return np.full((bins,), np.nan, dtype=np.float32) + return _fallback_depth_response(1, None) - s = np.asarray(sweep) - if s.size == 0: - return np.full((bins,), np.nan, dtype=np.float32) + try: + s = np.asarray(sweep, dtype=np.float64).ravel() + if s.size == 0: + return _fallback_depth_response(1, None) - sig = np.zeros(SWEEP_LEN, dtype=np.float32) - take = min(int(s.size), SWEEP_LEN) - seg = np.nan_to_num(s[:take], nan=0.0).astype(np.float32, copy=False) - sig[:take] = seg + freqs_hz = build_frequency_axis_hz(s.size) + s_complex = reconstruct_complex_spectrum_from_real_trace(s) + depth_m, y = perform_ifft_depth_response(s_complex, freqs_hz, axis="abs") + n = min(depth_m.size, y.size) + if n <= 0: + return _fallback_depth_response(s.size, s) + return depth_m[:n].astype(np.float32, copy=False), y[:n].astype(np.float32, copy=False) + except Exception as exc: # noqa: BLE001 + logger.error("compute_ifft_profile_from_sweep failed: %r", exc) + return _fallback_depth_response(np.asarray(sweep).size if sweep is not None else 1, sweep) - data = np.zeros(IFFT_LEN, dtype=np.complex64) - data[ZEROS_LOW + ZEROS_MID :] = sig - spec = np.fft.ifftshift(data) - result = np.fft.ifft(spec) - mag = np.abs(result).astype(np.float32) - return (mag + 1e-9).astype(np.float32) +def compute_ifft_db_profile(sweep: Optional[np.ndarray]) -> np.ndarray: + """Legacy wrapper (deprecated name): возвращает линейный |IFFT| профиль.""" + _depth_m, y = compute_ifft_profile_from_sweep(sweep) + return y + diff --git a/rfg_adc_plotter/state/app_state.py b/rfg_adc_plotter/state/app_state.py index 4adab5e..03913a3 100644 --- a/rfg_adc_plotter/state/app_state.py +++ b/rfg_adc_plotter/state/app_state.py @@ -227,7 +227,7 @@ class AppState: bg_stage = "bg[sub]" if self.background is not None else "bg[missing]" return ( f"pipeline ch{ch_txt}: parsed -> {reader_stage} -> raw -> " - f"live-calib-ref -> {calib_stage} -> {bg_stage} -> ring -> IFFT(dB)" + f"live-calib-ref -> {calib_stage} -> {bg_stage} -> ring -> IFFT(abs, depth_m)" ) calib_stage = self.format_calib_source_status() @@ -235,7 +235,7 @@ class AppState: return ( f"pipeline ch{ch_txt}: parsed -> {reader_stage} -> raw -> " - f"{calib_stage} -> {bg_stage} -> ring -> IFFT(dB)" + f"{calib_stage} -> {bg_stage} -> ring -> IFFT(abs, depth_m)" ) def _format_summary(self, summary) -> str: diff --git a/rfg_adc_plotter/state/ring_buffer.py b/rfg_adc_plotter/state/ring_buffer.py index 81e8355..d7dc09b 100644 --- a/rfg_adc_plotter/state/ring_buffer.py +++ b/rfg_adc_plotter/state/ring_buffer.py @@ -6,10 +6,11 @@ from typing import Optional, Tuple import numpy as np from rfg_adc_plotter.constants import ( - IFFT_LEN, + FREQ_MAX_GHZ, + FREQ_MIN_GHZ, WF_WIDTH, ) -from rfg_adc_plotter.processing.fourier import build_ifft_time_axis_ns, compute_ifft_db_profile +from rfg_adc_plotter.processing.fourier import compute_ifft_profile_from_sweep class RingBuffer: @@ -21,7 +22,8 @@ class RingBuffer: def __init__(self, max_sweeps: int): self.max_sweeps = max_sweeps - self.fft_bins = IFFT_LEN # = 1953 (полная длина IFFT-результата) + # Размер IFFT-профиля теперь динамический и определяется по первому успешному свипу. + self.fft_bins = 0 # Инициализируются при первом свипе (ensure_init) self.ring: Optional[np.ndarray] = None # (max_sweeps, WF_WIDTH) @@ -30,7 +32,7 @@ class RingBuffer: self.head: int = 0 self.width: Optional[int] = None self.x_shared: Optional[np.ndarray] = None - self.fft_time_axis: Optional[np.ndarray] = None # временная ось IFFT в нс + self.fft_depth_axis_m: Optional[np.ndarray] = None # ось глубины IFFT в метрах self.y_min_fft: Optional[float] = None self.y_max_fft: Optional[float] = None # FFT последнего свипа (для отображения без повторного вычисления) @@ -40,19 +42,21 @@ class RingBuffer: def is_ready(self) -> bool: return self.ring is not None + @property + def fft_time_axis(self) -> Optional[np.ndarray]: + """Legacy alias: старое имя поля (раньше было время в нс, теперь глубина в м).""" + return self.fft_depth_axis_m + def ensure_init(self, sweep_width: int): """Инициализировать буферы при первом свипе. Повторные вызовы — no-op (кроме x_shared).""" if self.ring is None: self.width = WF_WIDTH self.ring = np.full((self.max_sweeps, self.width), np.nan, dtype=np.float32) self.ring_time = np.full((self.max_sweeps,), np.nan, dtype=np.float64) - self.ring_fft = np.full((self.max_sweeps, self.fft_bins), np.nan, dtype=np.float32) - # Временная ось IFFT вынесена в processing.fourier для явного pipeline. - self.fft_time_axis = build_ifft_time_axis_ns() self.head = 0 # Обновляем x_shared если пришёл свип большего размера if self.x_shared is None or sweep_width > self.x_shared.size: - self.x_shared = np.linspace(3.323, 14.323, sweep_width, dtype=np.float32) + self.x_shared = np.linspace(FREQ_MIN_GHZ, FREQ_MAX_GHZ, sweep_width, dtype=np.float32) def push(self, s: np.ndarray): """Добавить строку свипа в кольцевой буфер, вычислить FFT-строку.""" @@ -69,8 +73,43 @@ class RingBuffer: self._push_fft(s) def _push_fft(self, s: np.ndarray): - fft_row = compute_ifft_db_profile(s) + depth_axis_m, fft_row = compute_ifft_profile_from_sweep(s) + fft_row = np.asarray(fft_row, dtype=np.float32).ravel() + depth_axis_m = np.asarray(depth_axis_m, dtype=np.float32).ravel() + n = min(int(fft_row.size), int(depth_axis_m.size)) + if n <= 0: + return + if n != fft_row.size: + fft_row = fft_row[:n] + if n != depth_axis_m.size: + depth_axis_m = depth_axis_m[:n] + + needs_reset = ( + self.ring_fft is None + or self.fft_depth_axis_m is None + or self.fft_bins != n + or self.ring_fft.shape != (self.max_sweeps, n) + or self.fft_depth_axis_m.size != n + ) + if (not needs_reset) and n > 0: + prev_axis = self.fft_depth_axis_m + assert prev_axis is not None + if prev_axis.size != n: + needs_reset = True + else: + # Если ось изменилась (например, изменилась длина/частотная сетка), сбрасываем FFT-водопад. + if not np.allclose(prev_axis[[0, -1]], depth_axis_m[[0, -1]], rtol=1e-6, atol=1e-9): + needs_reset = True + + if needs_reset: + self.fft_bins = n + self.ring_fft = np.full((self.max_sweeps, n), np.nan, dtype=np.float32) + self.fft_depth_axis_m = depth_axis_m.astype(np.float32, copy=True) + self.y_min_fft = None + self.y_max_fft = None + + assert self.ring_fft is not None prev_head = (self.head - 1) % self.ring_fft.shape[0] self.ring_fft[prev_head, :] = fft_row self.last_fft_vals = fft_row diff --git a/tests/test_fourier_phase_reconstruction.py b/tests/test_fourier_phase_reconstruction.py new file mode 100644 index 0000000..3cb909b --- /dev/null +++ b/tests/test_fourier_phase_reconstruction.py @@ -0,0 +1,75 @@ +import numpy as np + +from rfg_adc_plotter.processing.fourier import ( + build_frequency_axis_hz, + compute_ifft_profile_from_sweep, + normalize_sweep_for_phase, + perform_ifft_depth_response, + reconstruct_complex_spectrum_from_real_trace, + unwrap_arccos_phase_continuous, +) + + +def test_normalize_sweep_for_phase_max_abs_and_finite(): + sweep = np.array([np.nan, -10.0, 5.0, 20.0, -40.0, np.inf, -np.inf], dtype=np.float32) + x = normalize_sweep_for_phase(sweep) + assert x.dtype == np.float64 + assert np.all(np.isfinite(x)) + assert np.max(np.abs(x)) <= 1.0 + 1e-12 + + +def test_arccos_unwrap_continuous_recovers_complex_phase_without_large_jumps(): + phi_true = np.linspace(0.0, 4.0 * np.pi, 1000, dtype=np.float64) + x = np.cos(phi_true) + + phi_rec = unwrap_arccos_phase_continuous(x) + assert phi_rec.shape == phi_true.shape + assert np.max(np.abs(np.diff(phi_rec))) < 0.2 + + z_true = np.exp(1j * phi_true) + z_rec = np.exp(1j * phi_rec) + assert np.allclose(z_rec, z_true, atol=2e-2, rtol=0.0) + + +def test_reconstruct_complex_spectrum_from_real_trace_output_complex128(): + sweep = np.linspace(-1.0, 1.0, 64, dtype=np.float32) + z = reconstruct_complex_spectrum_from_real_trace(sweep) + assert z.dtype == np.complex128 + assert z.shape == sweep.shape + assert np.all(np.isfinite(np.real(z))) + assert np.all(np.isfinite(np.imag(z))) + + +def test_perform_ifft_depth_response_basic_abs(): + n = 128 + freqs = build_frequency_axis_hz(n) + s = np.exp(1j * np.linspace(0.0, 2.0 * np.pi, n, dtype=np.float64)) + + depth_m, y = perform_ifft_depth_response(s, freqs, axis="abs") + + assert depth_m.dtype == np.float32 + assert y.dtype == np.float32 + assert depth_m.ndim == 1 and y.ndim == 1 + assert depth_m.size == y.size + assert depth_m.size >= n + assert np.all(np.diff(depth_m) >= 0.0) + assert np.all(y >= 0.0) + + +def test_perform_ifft_depth_response_bad_grid_returns_fallback_not_exception(): + s = np.ones(16, dtype=np.complex128) + freqs_desc = np.linspace(10.0, 1.0, 16, dtype=np.float64) + depth_m, y = perform_ifft_depth_response(s, freqs_desc, axis="abs") + assert depth_m.size == y.size + assert depth_m.size == s.size + assert np.all(np.isfinite(depth_m)) + + +def test_compute_ifft_profile_from_sweep_returns_depth_and_linear_abs(): + sweep = np.linspace(-5.0, 7.0, 257, dtype=np.float32) + depth_m, y = compute_ifft_profile_from_sweep(sweep) + assert depth_m.dtype == np.float32 + assert y.dtype == np.float32 + assert depth_m.size == y.size + assert depth_m.size > 0 + assert np.all(np.diff(depth_m) >= 0.0) diff --git a/tests/test_ring_buffer_fft_axis.py b/tests/test_ring_buffer_fft_axis.py new file mode 100644 index 0000000..d145868 --- /dev/null +++ b/tests/test_ring_buffer_fft_axis.py @@ -0,0 +1,40 @@ +import numpy as np + +from rfg_adc_plotter.state.ring_buffer import RingBuffer + + +def test_ring_buffer_allocates_fft_buffers_from_first_push(): + ring = RingBuffer(max_sweeps=4) + ring.ensure_init(64) + + sweep = np.linspace(-1.0, 1.0, 64, dtype=np.float32) + ring.push(sweep) + + assert ring.ring_fft is not None + assert ring.fft_depth_axis_m is not None + assert ring.last_fft_vals is not None + assert ring.fft_bins == ring.ring_fft.shape[1] + assert ring.fft_bins == ring.fft_depth_axis_m.size + assert ring.fft_bins == ring.last_fft_vals.size + # Legacy alias kept for compatibility with existing GUI code paths. + assert ring.fft_time_axis is ring.fft_depth_axis_m + + +def test_ring_buffer_reallocates_fft_buffers_when_ifft_length_changes(): + ring = RingBuffer(max_sweeps=4) + ring.ensure_init(512) + + ring.push(np.linspace(-1.0, 1.0, 64, dtype=np.float32)) + first_bins = ring.fft_bins + first_shape = None if ring.ring_fft is None else ring.ring_fft.shape + + ring.push(np.linspace(-1.0, 1.0, 512, dtype=np.float32)) + second_bins = ring.fft_bins + second_shape = None if ring.ring_fft is None else ring.ring_fft.shape + + assert ring.ring is not None # raw ring сохраняется + assert first_shape is not None and second_shape is not None + assert first_bins != second_bins + assert second_shape == (ring.max_sweeps, second_bins) + assert ring.fft_depth_axis_m is not None + assert ring.fft_depth_axis_m.size == second_bins