From b057324e7ab6e0814516904a0b867a8dce76c256 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sat, 16 May 2026 22:06:12 -0700 Subject: [PATCH 01/25] Add prototype API for FM EK80 - Add public `compute_Sv_f()` and `compute_TS_f()` for the API - Add `_cal_complex_samples_f()` prototype - Expose `frequency_resolution` and `range_step` parameters for future spectral calibration - Return prototype `Sv_f` xarray output structure: `(channel, ping_time, svf_range, frequency)` --- echopype/calibrate/__init__.py | 4 +- echopype/calibrate/api.py | 64 +++++++++++++++++--- echopype/calibrate/calibrate_ek.py | 94 +++++++++++++++++++++++++++++- 3 files changed, 149 insertions(+), 13 deletions(-) diff --git a/echopype/calibrate/__init__.py b/echopype/calibrate/__init__.py index 94613c73b..0790f18eb 100644 --- a/echopype/calibrate/__init__.py +++ b/echopype/calibrate/__init__.py @@ -1,3 +1,3 @@ -from .api import compute_Sv, compute_TS +from .api import compute_Sv, compute_Sv_f, compute_TS, compute_TS_f -__all__ = ["compute_Sv", "compute_TS"] +__all__ = ["compute_Sv", "compute_TS", "compute_Sv_f", "compute_TS_f"] diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index d12b27b7c..8f5705874 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -21,7 +21,7 @@ def _compute_cal( - cal_type, + cal_type: str, echodata: EchoData, env_params=None, cal_params=None, @@ -30,6 +30,7 @@ def _compute_cal( encode_mode=None, assume_single_filter_time=None, drop_last_hanning_zero=False, + **kwargs, ): # Make waveform_mode "FM" equivalent to "BB" waveform_mode = "BB" if waveform_mode == "FM" else waveform_mode @@ -80,13 +81,19 @@ def _compute_cal_ds(echodata, slice_dict): # Check Echodata backscatter data size and recommend chunking if data is too large cal_obj._check_echodata_backscatter_size() - # Perform calibration - if cal_type == "Sv": - cal_ds = cal_obj.compute_Sv() - else: - cal_ds = cal_obj.compute_TS() + compute_methods = { + "Sv": cal_obj.compute_Sv, + "TS": cal_obj.compute_TS, + "Sv_f": cal_obj.compute_Sv_f, + "TS_f": cal_obj.compute_TS_f, + } + + try: + compute_method = compute_methods[cal_type] + except KeyError: + raise ValueError(f"Unsupported calibration type: {cal_type}") from None - return cal_ds + return compute_method(**kwargs) # Calibrate as a single dataset if not Ex80 if echodata.sonar_model not in ["EK80", "ES80", "EA640"]: @@ -201,12 +208,33 @@ def _add_attrs(cal_type, ds): """Add attributes to backscattering strength dataset. cal_type: Sv or TS """ - ds["range_sample"].attrs = {"long_name": "Along-range sample number, base 0"} - ds["echo_range"].attrs = {"long_name": "Range distance", "units": "m"} + if "range_sample" in ds: + ds["range_sample"].attrs = {"long_name": "Along-range sample number, base 0"} + + if "echo_range" in ds: + ds["echo_range"].attrs = { + "long_name": "Range distance", + "units": "m", + } + + if "svf_range" in ds: + ds["svf_range"].attrs = { + "long_name": "Range distance for Sv(f)/TS(f) window centres", + "units": "m", + } + + if "frequency" in ds: + ds["frequency"].attrs = { + "long_name": "Frequency", + "units": "Hz", + } + ds[cal_type].attrs = { "long_name": { "Sv": "Volume backscattering strength (Sv re 1 m-1)", "TS": "Target strength (TS re 1 m^2)", + "Sv_f": "Frequency-dependent volume backscattering strength (Sv(f) re 1 m-1)", + "TS_f": "Frequency-dependent target strength (TS(f) re 1 m^2)", }[cal_type], "units": "dB", } @@ -345,6 +373,15 @@ def compute_Sv(echodata: EchoData, **kwargs) -> xr.Dataset: return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs) +def compute_Sv_f(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute frequency-dependent volume backscattering strength Sv(f) + from broadband EK80 complex data. + TODO: add more details here when the function is fully implemented + """ + return _compute_cal(cal_type="Sv_f", echodata=echodata, **kwargs) + + def compute_TS(echodata: EchoData, **kwargs): """ Compute target strength (TS) from raw data. @@ -447,3 +484,12 @@ def compute_TS(echodata: EchoData, **kwargs): https://doi.org/10.1006/jmsc.2001.1158 """ return _compute_cal(cal_type="TS", echodata=echodata, **kwargs) + + +def compute_TS_f(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute frequency-dependent target strength TS(f) + from broadband EK80 complex data. + TODO: add more details here when the function is fully implemented + """ + return _compute_cal(cal_type="TS_f", echodata=echodata, **kwargs) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 3d638f973..11f8d6bad 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -658,7 +658,66 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: return out - def _compute_cal(self, cal_type) -> xr.Dataset: + def _cal_complex_samples_f( + self, + cal_type: str, + frequency_resolution: float = 1000, + range_step: float = 0.5, + ) -> xr.Dataset: + """Calibrate complex broadband EK80 data to Sv(f) or TS(f). + + Parameters + ---------- + cal_type : str + ``"Sv_f"`` for frequency-dependent volume backscattering strength, + or ``"TS_f"`` for frequency-dependent target strength. + frequency_resolution : float, default 1000 + Frequency spacing in Hz for the output spectral grid. + range_step : float, default 0.5 + Range spacing in meters for the output Sv(f) window centres. + + Returns + ------- + xr.Dataset + Prototype calibrated dataset containing ``Sv_f`` or ``TS_f``. + """ + + if self.waveform_mode not in ("FM", "BB") or self.encode_mode != "complex": + raise ValueError(f"{cal_type} is only supported for EK80 FM complex data.") + + if cal_type == "Sv_f": + out = xr.Dataset( + { + "Sv_f": ( + ["channel", "ping_time", "svf_range", "frequency"], + np.full( + ( + self.beam.sizes["channel"], + self.beam.sizes["ping_time"], + self.range_meter.sizes["range_sample"], + 1, + ), + np.nan, + ), + ) + }, + coords={ + "channel": self.beam["channel"], + "ping_time": self.beam["ping_time"], + "svf_range": self.range_meter.isel(channel=0, ping_time=0).values, + "frequency": [np.nan], + }, + ) + + elif cal_type == "TS_f": + raise NotImplementedError("TS_f not implemented yet.") + + else: + raise ValueError(f"Unsupported calibration type: {cal_type}") + + return out + + def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: """ Private method to compute Sv or TS from EK80 data, called by compute_Sv or compute_TS. @@ -678,7 +737,10 @@ def _compute_cal(self, cal_type) -> xr.Dataset: True if self.waveform_mode == "BB" or self.encode_mode == "complex" else False ) - if flag_complex: + if cal_type in ("Sv_f", "TS_f"): + # Frequency-dependent broadband calibration + ds_cal = self._cal_complex_samples_f(cal_type=cal_type, **kwargs) + elif flag_complex: # Complex samples can be BB or CW ds_cal = self._cal_complex_samples(cal_type=cal_type) else: @@ -708,3 +770,31 @@ def compute_TS(self): and the corresponding range (``echo_range``) in units meter. """ return self._compute_cal(cal_type="TS") + + def compute_Sv_f( + self, + frequency_resolution: float = 1000, + range_step: float = 0.5, + ): + """Compute frequency-dependent volume backscattering strength Sv(f). + + Returns + ------- + Sv_f : xr.Dataset + A Dataset containing frequency-dependent volume backscattering strength. + """ + return self._compute_cal( + cal_type="Sv_f", + frequency_resolution=frequency_resolution, + range_step=range_step, + ) + + def compute_TS_f(self): + """Compute frequency-dependent target strength TS(f). + + Returns + ------- + TS_f : xr.Dataset + A Dataset containing frequency-dependent target strength. + """ + return self._compute_cal(cal_type="TS_f") From c6886e35fcd45565d9b8e0cee8a0a9b248cd302d Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sun, 17 May 2026 06:52:01 -0700 Subject: [PATCH 02/25] Fix calibration dispatcher for unsupported FM methods Use lazy method getattr in calibration dispatcher to avoid access to unsupported compute_Sv_f / compute_TS_f methods on non-FM calibrators --- echopype/calibrate/api.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index 8f5705874..0fe92de91 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -82,17 +82,24 @@ def _compute_cal_ds(echodata, slice_dict): cal_obj._check_echodata_backscatter_size() compute_methods = { - "Sv": cal_obj.compute_Sv, - "TS": cal_obj.compute_TS, - "Sv_f": cal_obj.compute_Sv_f, - "TS_f": cal_obj.compute_TS_f, + "Sv": "compute_Sv", + "TS": "compute_TS", + "Sv_f": "compute_Sv_f", + "TS_f": "compute_TS_f", } try: - compute_method = compute_methods[cal_type] + method_name = compute_methods[cal_type] except KeyError: raise ValueError(f"Unsupported calibration type: {cal_type}") from None + compute_method = getattr(cal_obj, method_name, None) + + if compute_method is None: + raise ValueError( + f"{cal_type} calibration is not supported for " f"{echodata.sonar_model} data." + ) + return compute_method(**kwargs) # Calibrate as a single dataset if not Ex80 From 2aed1db9df9c33cd56326f902391f450c2c4b96c Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sat, 6 Jun 2026 18:58:40 -0700 Subject: [PATCH 03/25] major commit for FM TS - Add compute_TS_spectrum() for EK80 FM complex data - Add broadband TS spectrum processing pipeline and calibration helpers - Add CRIMAC reference dataset and validation tests for TS(f) - Add frequency-dependent absorption validation against CRIMAC - Add broadband Sp and Sv(f) integration tests - Add pooch support for ts_spectrum_example_data bundle - Refactor TS/Sp calibration workflow and API exports - Add cached test-data handling to avoid unnecessary downloads --- echopype/calibrate/__init__.py | 4 +- echopype/calibrate/api.py | 24 +- echopype/calibrate/calibrate_ek.py | 795 ++++++++++++++++-- echopype/calibrate/ek80_complex.py | 308 +++++++ echopype/tests/calibrate/test_calibrate.py | 12 +- .../tests/calibrate/test_calibrate_ek80.py | 1 + .../test_calibrate_ek80_broadband_crimac.py | 133 +++ echopype/tests/calibrate/test_ek80_complex.py | 134 ++- echopype/tests/conftest.py | 13 +- 9 files changed, 1330 insertions(+), 94 deletions(-) create mode 100644 echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py diff --git a/echopype/calibrate/__init__.py b/echopype/calibrate/__init__.py index 0790f18eb..7091c3504 100644 --- a/echopype/calibrate/__init__.py +++ b/echopype/calibrate/__init__.py @@ -1,3 +1,3 @@ -from .api import compute_Sv, compute_Sv_f, compute_TS, compute_TS_f +from .api import compute_Sp, compute_Sv, compute_Sv_f, compute_TS, compute_TS_spectrum -__all__ = ["compute_Sv", "compute_TS", "compute_Sv_f", "compute_TS_f"] +__all__ = ["compute_Sv", "compute_TS", "compute_Sp", "compute_Sv_f", "compute_TS_spectrum"] diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index 0fe92de91..4ed4582d4 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -83,9 +83,10 @@ def _compute_cal_ds(echodata, slice_dict): compute_methods = { "Sv": "compute_Sv", + "Sp": "compute_Sp", "TS": "compute_TS", "Sv_f": "compute_Sv_f", - "TS_f": "compute_TS_f", + "TS_spectrum": "compute_TS_spectrum", } try: @@ -239,9 +240,10 @@ def _add_attrs(cal_type, ds): ds[cal_type].attrs = { "long_name": { "Sv": "Volume backscattering strength (Sv re 1 m-1)", + "Sp": "Point scattering strength (Sp re 1 m^2)", "TS": "Target strength (TS re 1 m^2)", "Sv_f": "Frequency-dependent volume backscattering strength (Sv(f) re 1 m-1)", - "TS_f": "Frequency-dependent target strength (TS(f) re 1 m^2)", + "TS_spectrum": "Frequency-dependent target strength spectrum (TS(f) re 1 m^2)", }[cal_type], "units": "dB", } @@ -389,6 +391,17 @@ def compute_Sv_f(echodata: EchoData, **kwargs) -> xr.Dataset: return _compute_cal(cal_type="Sv_f", echodata=echodata, **kwargs) +def compute_Sp(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute point scattering strength (Sp) from raw data. + + For CW data, Sp is computed from received power samples on the range grid. + For EK80 broadband/FM complex data, Sp is computed after pulse compression + and represents a band-averaged point-scattering-strength echogram. + """ + return _compute_cal(cal_type="Sp", echodata=echodata, **kwargs) + + def compute_TS(echodata: EchoData, **kwargs): """ Compute target strength (TS) from raw data. @@ -493,10 +506,9 @@ def compute_TS(echodata: EchoData, **kwargs): return _compute_cal(cal_type="TS", echodata=echodata, **kwargs) -def compute_TS_f(echodata: EchoData, **kwargs) -> xr.Dataset: +def compute_TS_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: """ - Compute frequency-dependent target strength TS(f) - from broadband EK80 complex data. + Compute broadband frequency-dependent target strength spectrum. TODO: add more details here when the function is fully implemented """ - return _compute_cal(cal_type="TS_f", echodata=echodata, **kwargs) + return _compute_cal(cal_type="TS_spectrum", echodata=echodata, **kwargs) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index b15f38404..9ce328c61 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -2,17 +2,29 @@ import numpy as np import xarray as xr +from tqdm.auto import tqdm from ..echodata import EchoData from ..echodata.simrad import retrieve_correct_beam_group +from ..utils import uwa from ..utils.log import _init_logger from .cal_params import _get_interp_da, get_cal_params_EK from .calibrate_base import CalibrateBase from .ecs import conform_channel_order, ecs_ds2dict, ecs_ev2ep from .ek80_complex import ( - compress_pulse, + _align_autocorrelation, + _compute_power_from_complex_signal, + _compute_svf_power, + _compute_svf_spectrum, + _compute_ts_spectrum, + _compute_ts_spectrum_calibrated, + _compute_ts_spectrum_power, + _get_autocorrelation, + _get_average_signal, + _get_hanning_window, + _get_pulse_compressed_signal, + _get_svf_frequency_grid, get_filter_coeff, - get_norm_fac, get_tau_effective, get_transmit_signal, ) @@ -171,17 +183,42 @@ def _cal_power_samples(self, cal_type: str) -> xr.Dataset: ) out.name = "Sv" - elif cal_type == "TS": - # Calc gain + elif cal_type in ("Sp", "TS"): CSp = ( 10 * np.log10(self.beam["transmit_power"]) + 2 * self.cal_params["gain_correction"] + 10 * np.log10(wavelength**2 / (16 * np.pi**2)) ) - # Calibration and echo integration - out = self.beam["backscatter_r"] + spreading_loss * 2 + absorption_loss - CSp - out.name = "TS" + sp = self.beam["backscatter_r"] + spreading_loss * 2 + absorption_loss - CSp + + if cal_type == "TS": + angle_alongship = self.beam["angle_alongship"] + angle_athwartship = self.beam["angle_athwartship"] + + angle_offset_alongship = self.cal_params["angle_offset_alongship"] + angle_offset_athwartship = self.cal_params["angle_offset_athwartship"] + beamwidth_alongship = self.cal_params["beamwidth_alongship"] + beamwidth_athwartship = self.cal_params["beamwidth_athwartship"] + + fac_along = ( + np.abs(angle_alongship - angle_offset_alongship) / (beamwidth_alongship / 2) + ) ** 2 + + fac_athwart = ( + np.abs(angle_athwartship - angle_offset_athwartship) + / (beamwidth_athwartship / 2) + ) ** 2 + + beam_correction = ( + 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) + ) + + out = sp + beam_correction + else: + out = sp + + out.name = cal_type # Attach calculated range (with units meter) into data set out = out.to_dataset() @@ -264,6 +301,9 @@ def compute_Sv(self, **kwargs): def compute_TS(self, **kwargs): return self._cal_power_samples(cal_type="TS") + def compute_Sp(self, **kwargs): + return self._cal_power_samples(cal_type="Sp") + class CalibrateEK80(CalibrateEK): # Default EK80 params: these parameters are only recorded in later versions of EK80 software @@ -480,25 +520,19 @@ def _get_power_from_complex( Power computed from complex samples """ - def _get_prx(sig): - return ( - beam["beam"].size # number of transducer sectors - * np.abs(sig.mean(dim="beam")) ** 2 - / (2 * np.sqrt(2)) ** 2 - * (np.abs(z_er + z_et) / z_er) ** 2 - / z_et - ) - - # Compute power if self.waveform_mode == "BB": - pc = compress_pulse( - backscatter=beam["backscatter_r"] + 1j * beam["backscatter_i"], chirp=chirp - ) # has beam dim - pc = pc / get_norm_fac(chirp=chirp) # normalization for each channel - prx = _get_prx(pc) # ensure prx is xr.DataArray + signal = _get_pulse_compressed_signal( + beam=beam, + matched_filter=chirp, + ) else: - bs_cw = beam["backscatter_r"] + 1j * beam["backscatter_i"] - prx = _get_prx(bs_cw) + signal = beam["backscatter_r"] + 1j * beam["backscatter_i"] + + prx = _compute_power_from_complex_signal( + signal=signal, + z_et=z_et, + z_er=z_er, + ) prx.name = "received_power" @@ -626,16 +660,20 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: out.name = "Sv" # out = out.rename_vars({list(out.data_vars.keys())[0]: "Sv"}) + # remove TS below and create a TS that build on Sp + elif cal_type in ("Sp", "TS"): + range_safe = self._safe_range_for_log(range_meter) + + spreading_loss_safe = 20 * np.log10(range_safe) - elif cal_type == "TS": out = ( 10 * np.log10(prx) - + 2 * spreading_loss + + 2 * spreading_loss_safe + absorption_loss - 10 * np.log10(wavelength**2 * transmit_power / (16 * np.pi**2)) - 2 * gain ) - out.name = "TS" + out.name = cal_type # change # Attach calculated range (with units meter) into data set out = out.to_dataset().merge(range_meter) @@ -658,114 +696,699 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: return out - def _cal_complex_samples_f( + ### Common helpers for complex sample calibration and Sv(f)/TS(f) calibration + + def _select_param( self, - cal_type: str, + da: xr.DataArray, + channel=None, + ping_idx=None, + ): + """Select channel- and ping-specific parameter values when dimensions exist.""" + if channel is not None and "channel" in da.dims: + da = da.sel(channel=channel) + if ping_idx is not None and "ping_time" in da.dims: + da = da.isel(ping_time=ping_idx) + return da + + def _safe_range_for_log(self, range_meter: xr.DataArray) -> xr.DataArray: + """Avoid log10(0) in range-dependent calibration terms.""" + return range_meter.where(range_meter > 0, 1e-20) + + def _compute_absorption_f( + self, + frequency: np.ndarray, + channel: str, + ping_idx: int, + sound_speed: float, + ) -> np.ndarray: + """Compute frequency-dependent absorption for broadband spectra. + + Uses the Francois & Garrison (1982) formulation, matching the + absorption model used by CRIMAC ``calc_alpha``. + """ + return uwa.calc_absorption( + frequency=frequency, + temperature=float( + self._select_param(self.env_params["temperature"], channel, ping_idx) + ), + salinity=float(self._select_param(self.env_params["salinity"], channel, ping_idx)), + pressure=float(self._select_param(self.env_params["pressure"], channel, ping_idx)), + pH=float(self._select_param(self.env_params["pH"], channel, ping_idx)), + sound_speed=sound_speed, + formula_source="FG", + ) + + ### Helpers for Svf + + def _compute_svf( + self, + power_spectrum: np.ndarray, + frequency: np.ndarray, + svf_range: np.ndarray, + sound_speed: float, + absorption_f: np.ndarray, + transmit_power: float, + psi_f: np.ndarray, + gain_f: np.ndarray, + window_duration: float, + ): + """Apply CRIMAC/pyEcholab-style Sv(f) calibration equation. + + Equivalent to CRIMAC ``calcSvf``. + """ + wavelength_f = sound_speed / frequency + + G = ( + transmit_power * wavelength_f**2 * sound_speed * window_duration * psi_f * gain_f**2 + ) / (32 * np.pi**2) + + return ( + 10 * np.log10(power_spectrum) + + 2 * absorption_f[None, :] * svf_range[:, None] + - 10 * np.log10(G[None, :]) + ) + + ### Helpers for TSf + + def _get_beam_compensated_gain( + self, + channel: str, + theta: float, + phi: float, + frequency: np.ndarray, + ) -> np.ndarray: + """Calculate beam-compensated gain. + + Equivalent to CRIMAC ``calc_g(theta, phi, f)``. + """ + gain_db = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["gain"].sel(cal_channel_id=channel).values, + ) + + angle_offset_alongship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["angle_offset_alongship"].sel(cal_channel_id=channel).values, + ) + + angle_offset_athwartship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["angle_offset_athwartship"].sel(cal_channel_id=channel).values, + ) + + beamwidth_alongship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["beamwidth_alongship"].sel(cal_channel_id=channel).values, + ) + + beamwidth_athwartship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["beamwidth_athwartship"].sel(cal_channel_id=channel).values, + ) + + fac_along = (np.abs(theta - angle_offset_alongship) / (beamwidth_alongship / 2)) ** 2 + + fac_athwart = (np.abs(phi - angle_offset_athwartship) / (beamwidth_athwartship / 2)) ** 2 + + beam_correction_db = ( + 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) + ) + + return 10 ** ((gain_db - beam_correction_db) / 10) + + ### Sv(f)/TS(f) calibration pipelines + + def _cal_complex_samples_svf( + self, + pc: xr.DataArray, + matched_filter: Dict, frequency_resolution: float = 1000, range_step: float = 0.5, ) -> xr.Dataset: - """Calibrate complex broadband EK80 data to Sv(f) or TS(f). + """Compute frequency-dependent volume backscattering strength Sv(f).""" - Parameters - ---------- - cal_type : str - ``"Sv_f"`` for frequency-dependent volume backscattering strength, - or ``"TS_f"`` for frequency-dependent target strength. - frequency_resolution : float, default 1000 - Frequency spacing in Hz for the output spectral grid. - range_step : float, default 0.5 - Range spacing in meters for the output Sv(f) window centres. + pc_avg = _get_average_signal(pc) - Returns - ------- - xr.Dataset - Prototype calibrated dataset containing ``Sv_f`` or ``TS_f``. - """ + out_by_channel = [] - if self.waveform_mode not in ("FM", "BB") or self.encode_mode != "complex": - raise ValueError(f"{cal_type} is only supported for EK80 FM complex data.") + for channel in self.beam["channel"].values: + if channel not in self.vend["cal_channel_id"].values: + continue - if cal_type == "Sv_f": - out = xr.Dataset( + pc_avg_ch = pc_avg.sel(channel=channel) + range_ch = self.range_meter.sel(channel=channel) + + svf_list = [] + svf_ping_time = [] + frequency_ref = None + svf_range_ref = None + + for ping_idx in tqdm(range(pc_avg_ch.sizes["ping_time"]), desc=f"{channel}"): + pc_avg_1d = pc_avg_ch.isel(ping_time=ping_idx) + range_1d = range_ch.isel(ping_time=ping_idx) + + valid = np.isfinite(pc_avg_1d) & np.isfinite(range_1d) + + if not bool(valid.any()): + continue + + pc_avg_1d = pc_avg_1d.where(valid, drop=True).values + range_1d = range_1d.where(valid, drop=True).values + + sound_speed = float( + self._select_param(self.env_params["sound_speed"], channel, ping_idx) + ) + transmit_power = float( + self._select_param(self.beam["transmit_power"], channel, ping_idx) + ) + tau = float( + self._select_param(self.beam["transmit_duration_nominal"], channel, ping_idx) + ) + sample_interval = float( + self._select_param(self.beam["sample_interval"], channel, ping_idx) + ) + f_start = float( + self._select_param(self.beam["transmit_frequency_start"], channel, ping_idx) + ) + f_stop = float( + self._select_param(self.beam["transmit_frequency_stop"], channel, ping_idx) + ) + equivalent_beam_angle = float( + self._select_param(self.cal_params["equivalent_beam_angle"], channel, ping_idx) + ) + z_et = float( + self._select_param(self.cal_params["impedance_transducer"], channel, ping_idx) + ) + z_er = float( + self._select_param(self.cal_params["impedance_transceiver"], channel, None) + ) + + window, n_window, window_duration, fs_dec = _get_hanning_window( + sound_speed=sound_speed, + tau=tau, + sample_interval=sample_interval, + ) + + frequency, frequency_index = _get_svf_frequency_grid( + f_start=f_start, + f_stop=f_stop, + frequency_resolution=frequency_resolution, + fs_dec=fs_dec, + n_window=n_window, + ) + + _, mf_auto_spectrum = _get_autocorrelation( + matched_filter=matched_filter[channel], + n_window=n_window, + ) + + step = max( + 1, + int(np.round(range_step / np.nanmedian(np.diff(range_1d)))), + ) + + pc_spread = pc_avg_1d * range_1d + + normalized_spectrum, svf_range = _compute_svf_spectrum( + pc_spread=pc_spread, + range_meter=range_1d, + window=window, + n_window=n_window, + frequency_index=frequency_index, + mf_auto_spectrum=mf_auto_spectrum, + step=step, + ) + + if normalized_spectrum.size == 0: + continue + + power_spectrum = _compute_svf_power( + normalized_spectrum=normalized_spectrum, + n_beams=self.beam["beam"].size, + z_et=z_et, + z_er=z_er, + ) + + gain_db = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["gain"].sel(cal_channel_id=channel).values, + ) + gain_f = 10 ** (gain_db / 10) + + frequency_nominal = float( + self._select_param(self.beam["frequency_nominal"], channel, None) + ) + psi_f = 10 ** (equivalent_beam_angle / 10) * (frequency_nominal / frequency) ** 2 + + absorption_f = self._compute_absorption_f( + frequency=frequency, + channel=channel, + ping_idx=ping_idx, + sound_speed=sound_speed, + ) + + svf = self._compute_svf( + power_spectrum=power_spectrum, + frequency=frequency, + svf_range=svf_range, + sound_speed=sound_speed, + absorption_f=absorption_f, + transmit_power=transmit_power, + psi_f=psi_f, + gain_f=gain_f, + window_duration=window_duration, + ) + + if frequency_ref is None: + frequency_ref = frequency + svf_range_ref = svf_range + elif ( + frequency.shape != frequency_ref.shape or svf_range.shape != svf_range_ref.shape + ): + continue + + svf_list.append(svf) + svf_ping_time.append(self.beam["ping_time"].isel(ping_time=ping_idx).values) + + if not svf_list: + continue + + ds_ch = xr.Dataset( { "Sv_f": ( ["channel", "ping_time", "svf_range", "frequency"], - np.full( - ( - self.beam.sizes["channel"], - self.beam.sizes["ping_time"], - self.range_meter.sizes["range_sample"], - 1, - ), - np.nan, - ), + np.asarray(svf_list)[None, :, :, :], ) }, coords={ - "channel": self.beam["channel"], - "ping_time": self.beam["ping_time"], - "svf_range": self.range_meter.isel(channel=0, ping_time=0).values, - "frequency": [np.nan], + "channel": [channel], + "ping_time": svf_ping_time, + "svf_range": svf_range_ref, + "frequency": frequency_ref, }, ) - elif cal_type == "TS_f": - raise NotImplementedError("TS_f not implemented yet.") + out_by_channel.append(ds_ch) - else: - raise ValueError(f"Unsupported calibration type: {cal_type}") + if not out_by_channel: + raise ValueError("No valid Sv(f) pings produced.") - return out + return xr.concat(out_by_channel, dim="channel") + + def _cal_complex_samples_TS_spectrum( + self, + pc: xr.DataArray, + matched_filter: Dict, + point_locations: xr.Dataset, + NFFT: int | None = None, + n_f_points: int | None = None, + split_front: float = 0.25, + window: str | None = None, + frequency_resolution: float | None = None, + ) -> xr.Dataset: + """Compute frequency-dependent target strength spectrum.""" + + if point_locations is None: + raise ValueError("TS_spectrum requires a point_locations dataset.") + + if not 0 <= split_front <= 1: + raise ValueError("split_front must be between 0 and 1.") + + pc_avg = _get_average_signal(pc) + + out_by_channel = [] + + for channel in self.beam["channel"].values: + if channel not in self.vend["cal_channel_id"].values: + continue + + points_ch = point_locations.where( + point_locations["channel"] == channel, + drop=True, + ) + + if points_ch.sizes.get("target_id", 0) == 0: + continue + + pc_avg_ch = pc_avg.sel(channel=channel) + range_ch = self.range_meter.sel(channel=channel) + + ts_list = [] + target_range_list = [] + theta_list = [] + phi_list = [] + ping_time_list = [] + target_id_list = [] + frequency_ref = None + + for point_idx in tqdm(range(points_ch.sizes["target_id"]), desc=f"{channel}"): + point_ping_time = points_ch["ping_time"].isel(target_id=point_idx).values + target_range = float(points_ch["target_range"].isel(target_id=point_idx)) + + theta_t = float(points_ch["angle_alongship"].isel(target_id=point_idx)) + phi_t = float(points_ch["angle_athwartship"].isel(target_id=point_idx)) + target_id = points_ch["target_id"].isel(target_id=point_idx).values + + ping_idx = int( + np.argmin( + np.abs( + self.beam["ping_time"].values.astype("datetime64[ns]") + - np.datetime64(point_ping_time, "ns") + ) + ) + ) + + pc_avg_1d = pc_avg_ch.isel(ping_time=ping_idx) + range_1d = range_ch.isel(ping_time=ping_idx) + + valid = np.isfinite(pc_avg_1d) & np.isfinite(range_1d) + + if not bool(valid.any()): + continue + + pc_avg_1d = pc_avg_1d.where(valid, drop=True).values + range_1d = range_1d.where(valid, drop=True).values + + sound_speed = float( + self._select_param(self.env_params["sound_speed"], channel, ping_idx) + ) + transmit_power = float( + self._select_param(self.beam["transmit_power"], channel, ping_idx) + ) + sample_interval = float( + self._select_param(self.beam["sample_interval"], channel, ping_idx) + ) + f_start = float( + self._select_param(self.beam["transmit_frequency_start"], channel, ping_idx) + ) + f_stop = float( + self._select_param(self.beam["transmit_frequency_stop"], channel, ping_idx) + ) + + if frequency_resolution is not None: + n_f_points_local = int(np.floor((f_stop - f_start) / frequency_resolution)) + 1 + else: + n_f_points_local = n_f_points if n_f_points is not None else 1000 + + frequency = np.linspace(f_start, f_stop, n_f_points_local) + + if NFFT is None: + n_fft = int(2 ** np.ceil(np.log2(n_f_points_local))) + else: + n_fft = NFFT + + z_et = float( + self._select_param(self.cal_params["impedance_transducer"], channel, ping_idx) + ) + z_er = float( + self._select_param(self.cal_params["impedance_transceiver"], channel, None) + ) + + fs_dec = 1 / sample_interval + + absorption_f = self._compute_absorption_f( + frequency=frequency, + channel=channel, + ping_idx=ping_idx, + sound_speed=sound_speed, + ) + + ########## try to fix something ??????? + + target_range_min = float(points_ch["target_range_min"].isel(target_id=point_idx)) + target_range_max = float(points_ch["target_range_max"].isel(target_id=point_idx)) + + target_mask = (range_1d >= target_range_min) & (range_1d <= target_range_max) + + if not np.any(target_mask): + continue + + pc_target = pc_avg_1d[target_mask] + + if window == "hann": + win = np.hanning(pc_target.size) + elif window == "hamming": + win = np.hamming(pc_target.size) + elif window in (None, "boxcar", "rectangular"): + win = np.ones(pc_target.size) + else: + raise ValueError(f"Unsupported window: {window}") + + pc_target = pc_target * win + + ########## + + mf_auto, _ = _get_autocorrelation( + matched_filter=matched_filter[channel], + n_window=n_fft, + ) + + mf_auto_red = _align_autocorrelation( + mf_auto=mf_auto, + pc_target=pc_target, + ) + + if mf_auto_red.size < n_fft: + mf_pad = np.zeros(n_fft, dtype=complex) + mf_pad[: mf_auto_red.size] = mf_auto_red + mf_auto_red = mf_pad + elif mf_auto_red.size > n_fft: + mf_auto_red = mf_auto_red[:n_fft] + + _, _, normalized_spectrum = _compute_ts_spectrum( + pc_target=pc_target, + mf_auto_red=mf_auto_red, + NFFT=n_fft, + frequency=frequency, + fs_dec=fs_dec, + ) + + power_spectrum = _compute_ts_spectrum_power( + normalized_spectrum=normalized_spectrum, + n_beams=self.beam["beam"].size, + z_et=z_et, + z_er=z_er, + ) + + gain_f = self._get_beam_compensated_gain( + channel=channel, + theta=theta_t, + phi=phi_t, + frequency=frequency, + ) + + ts = _compute_ts_spectrum_calibrated( + power_spectrum=power_spectrum, + target_range=target_range, + frequency=frequency, + sound_speed=sound_speed, + absorption_f=absorption_f, + transmit_power=transmit_power, + gain_f=gain_f, + ) + + if frequency_ref is None: + frequency_ref = frequency + elif frequency.shape != frequency_ref.shape: + continue + + ts_list.append(ts) + target_range_list.append(target_range) + theta_list.append(theta_t) + phi_list.append(phi_t) + ping_time_list.append(point_ping_time) + target_id_list.append(target_id) + + if not ts_list: + continue + + ds_ch = xr.Dataset( + { + "TS_spectrum": ( + ["channel", "target_id", "frequency"], + np.asarray(ts_list)[None, :, :], + ), + "target_range": ( + ["channel", "target_id"], + np.asarray(target_range_list)[None, :], + ), + "angle_alongship": ( + ["channel", "target_id"], + np.asarray(theta_list)[None, :], + ), + "angle_athwartship": ( + ["channel", "target_id"], + np.asarray(phi_list)[None, :], + ), + "ping_time": ( + ["channel", "target_id"], + np.asarray(ping_time_list)[None, :], + ), + }, + coords={ + "channel": [channel], + "target_id": np.asarray(target_id_list), + "frequency": frequency_ref, + }, + ) + + out_by_channel.append(ds_ch) + + if not out_by_channel: + raise ValueError("No valid TS_spectrum targets produced.") + + return xr.concat(out_by_channel, dim="channel") + + ### Dispatchers + + def _cal_complex_samples_f( + self, + cal_type: str, + frequency_resolution: float | None = None, + range_step: float = 0.5, + point_locations: xr.Dataset | None = None, + NFFT: int | None = None, + n_f_points: int | None = None, + split_front: float = 0.25, + window: str | None = None, + ) -> xr.Dataset: + """Calibrate EK80 FM complex data to frequency-dependent Sv(f) or TS(f).""" + if self.waveform_mode not in ("FM", "BB") or self.encode_mode != "complex": + raise ValueError(f"{cal_type} is only supported for EK80 FM complex data.") + + tx_coeff = get_filter_coeff(self.vend) + fs = self.cal_params["receiver_sampling_frequency"] + + matched_filter, _ = get_transmit_signal( + self.beam, + tx_coeff, + self.waveform_mode, + fs, + self.drop_last_hanning_zero, + ) + + pc = _get_pulse_compressed_signal( + beam=self.beam, + matched_filter=matched_filter, + ) + + if cal_type == "Sv_f": + return self._cal_complex_samples_svf( + pc=pc, + matched_filter=matched_filter, + frequency_resolution=frequency_resolution, + range_step=range_step, + ) + + if cal_type == "TS_spectrum": + return self._cal_complex_samples_TS_spectrum( + pc=pc, + matched_filter=matched_filter, + point_locations=point_locations, + NFFT=NFFT, + n_f_points=n_f_points, + split_front=split_front, + window=window, + frequency_resolution=frequency_resolution, + ) + + raise ValueError(f"Unsupported calibration type: {cal_type}") def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: """ - Private method to compute Sv or TS from EK80 data, called by compute_Sv or compute_TS. + Private dispatcher for EK80 calibration. + + This routes calibration to one of three paths: + + 1. Power-sample CW calibration: + Used for EK60 and EK80 power-encoded CW data. + + 2. Complex-sample calibration: + Used for EK80 complex data in CW or BB/FM mode. For BB/FM data, + this path computes the conventional broadband-averaged Sv, Sp, or + legacy TS-like gridded product: pulse compression is applied, + transducer sectors are averaged, received power is computed, and the + standard calibration equation is used. + + 3. Frequency-dependent complex-sample calibration: + Used for EK80 BB/FM complex data when computing Sv(f) or TS(f). + This path keeps the pulse-compressed signal before final calibration + so that FFT-based spectral processing can be applied. Parameters ---------- cal_type : str - 'Sv' for calculating volume backscattering strength, or - 'TS' for calculating target strength + Calibration type. Supported values include ``"Sv"``, ``"Sp"``, + ``"TS"``, ``"Sv_f"``, and ``"TS_spectrum"``. Returns ------- xr.Dataset - An xarray Dataset containing either Sv or TS. + Dataset containing the calibrated output requested by ``cal_type``. """ # Set flag_complex: True-complex cal, False-power cal flag_complex = ( True if self.waveform_mode == "BB" or self.encode_mode == "complex" else False ) - if cal_type in ("Sv_f", "TS_f"): - # Frequency-dependent broadband calibration + if cal_type in ("Sv_f", "TS_spectrum"): + # Frequency-dependent broadband calibration: keep pulse-compressed + # signal for FFT-based Sv(f) or TS(f) processing. ds_cal = self._cal_complex_samples_f(cal_type=cal_type, **kwargs) elif flag_complex: - # Complex samples can be BB or CW + # Complex-sample calibration: BB/FM data are pulse-compressed and + # averaged over transducer sectors before computing conventional + # Sv/Sp/legacy gridded TS. CW complex data are calibrated directly + # from complex samples. ds_cal = self._cal_complex_samples(cal_type=cal_type) else: - # Power samples only make sense for CW mode data + # Power-sample calibration: applies to power-encoded CW data. ds_cal = self._cal_power_samples(cal_type=cal_type) return ds_cal + ### Public API + def compute_Sv(self): """Compute volume backscattering strength (Sv). Returns ------- - Sv : xr.DataSet + Sv : xr.Dataset A DataSet containing volume backscattering strength (``Sv``) and the corresponding range (``echo_range``) in units meter. """ return self._compute_cal(cal_type="Sv") + def compute_Sp(self): + """Compute point scattering strength Sp. + + For CW data, this is computed from the received power samples on the + range grid. For FM complex data, this is computed after pulse + compression and represents a band-averaged point-scattering-strength + echogram. + + Returns + ------- + Sp : xr.Dataset + Dataset containing point scattering strength (``Sp``) and the + corresponding range (``echo_range``) in metres. + """ + return self._compute_cal(cal_type="Sp") + def compute_TS(self): """Compute target strength (TS). Returns ------- - TS : xr.DataSet + TS : xr.Dataset A DataSet containing target strength (``TS``) and the corresponding range (``echo_range``) in units meter. """ @@ -789,7 +1412,15 @@ def compute_Sv_f( range_step=range_step, ) - def compute_TS_f(self): + def compute_TS_spectrum( + self, + point_locations: xr.Dataset, + NFFT: int | None = None, + n_f_points: int | None = None, + split_front: float = 0.25, + window: str = None, + frequency_resolution: float | None = None, + ): """Compute frequency-dependent target strength TS(f). Returns @@ -797,4 +1428,12 @@ def compute_TS_f(self): TS_f : xr.Dataset A Dataset containing frequency-dependent target strength. """ - return self._compute_cal(cal_type="TS_f") + return self._compute_cal( + cal_type="TS_spectrum", + point_locations=point_locations, + NFFT=NFFT, + n_f_points=n_f_points, + split_front=split_front, + window=window, + frequency_resolution=frequency_resolution, + ) diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py index 1717083b0..bd259ca60 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -9,6 +9,314 @@ from ..convert.set_groups_ek80 import DECIMATION, FILTER_IMAG, FILTER_REAL +def _compute_power_from_complex_signal( + signal: xr.DataArray, + z_et, + z_er, +) -> xr.DataArray: + """Calculate received electrical power from sector-level complex samples. + + The input is expected to retain the ``beam`` dimension. The function + averages over transducer sectors internally before converting the + complex signal to received electrical power. + + Equivalent to CRIMAC ``calcPower``. + """ + prx = ( + signal["beam"].size + * np.abs(signal.mean(dim="beam")) ** 2 + / (2 * np.sqrt(2)) ** 2 + * (np.abs(z_er + z_et) / np.abs(z_er)) ** 2 + / np.abs(z_et) + ) + + prx = prx.where(prx > 0, 1e-20) + prx.name = "received_power" + + return prx + + +def _align_autocorrelation( + mf_auto: np.ndarray, + pc_target: np.ndarray, +) -> np.ndarray: + """Align matched-filter autocorrelation to target echo. + + Equivalent to CRIMAC ``alignAuto``. + """ + idx_peak_auto = np.argmax(np.abs(mf_auto)) + idx_peak_target = np.argmax(np.abs(pc_target)) + + left_samples = idx_peak_target + right_samples = len(pc_target) - idx_peak_target + + idx_start = max(0, idx_peak_auto - left_samples) + idx_stop = min(len(mf_auto), idx_peak_auto + right_samples) + + return mf_auto[idx_start:idx_stop] + + +def _compute_ts_spectrum( + pc_target: np.ndarray, + mf_auto_red: np.ndarray, + NFFT: int, + frequency: np.ndarray, + fs_dec: float, +): + """Compute target, autocorrelation, and normalised DFTs for TS spectrum. + + Equivalent to CRIMAC ``calcDFTforTS``, with explicit NFFT. + """ + frequency_index = np.mod( + np.floor(frequency / fs_dec * NFFT).astype(int), + NFFT, + ) + + pc_target_spectrum = np.fft.fft(pc_target, n=NFFT)[frequency_index] + mf_auto_red_spectrum = np.fft.fft(mf_auto_red, n=NFFT)[frequency_index] + + normalized_spectrum = pc_target_spectrum / mf_auto_red_spectrum + + return pc_target_spectrum, mf_auto_red_spectrum, normalized_spectrum + + +def _get_hanning_window( + sound_speed: float, + tau: float, + sample_interval: float, +): + """Get normalized Hann window for Sv(f) sliding-window FFT. + + Equivalent to CRIMAC ``defHanningWindow``. + """ + dr = sample_interval * sound_speed / 2 + fs_dec = 1 / sample_interval + + L = sound_speed * 2 * tau / dr + n_window = int(2 ** np.ceil(np.log2(L))) + + window = np.hanning(n_window) + window = window / (np.linalg.norm(window) / np.sqrt(n_window)) + + window_duration = n_window / fs_dec + + return window, n_window, window_duration, fs_dec + + +def _compute_svf_spectrum( + pc_spread: np.ndarray, + range_meter: np.ndarray, + window: np.ndarray, + n_window: int, + frequency_index: np.ndarray, + mf_auto_spectrum: np.ndarray, + step: int = 1, +): + """Compute normalized sliding-window spectrum for Sv(f). + + Equivalent to CRIMAC ``calcDFTforSv``. + """ + mf_auto_spectrum_f = mf_auto_spectrum[frequency_index] + + normalized_spectrum = [] + svf_range = [] + + for start in range(0, len(pc_spread) - n_window, step): + stop = start + n_window + + windowed_signal = window * pc_spread[start:stop] + signal_spectrum = np.fft.fft(windowed_signal, n=n_window)[frequency_index] + + normalized_spectrum.append(signal_spectrum / mf_auto_spectrum_f) + + center = (start + stop) // 2 + svf_range.append(range_meter[center]) + + return np.asarray(normalized_spectrum), np.asarray(svf_range) + + +def _compute_ts_spectrum_calibrated( + power_spectrum: np.ndarray, + target_range: float, + frequency: np.ndarray, + sound_speed: float, + absorption_f: np.ndarray, + transmit_power: float, + gain_f: np.ndarray, +): + """Apply CRIMAC-style TS(f) calibration equation. + + Equivalent to CRIMAC ``calcTSf``. + """ + wavelength_f = sound_speed / frequency + + return ( + 10 * np.log10(power_spectrum) + + 40 * np.log10(target_range) + + 2 * absorption_f * target_range + - 10 * np.log10(transmit_power * wavelength_f**2 * gain_f**2 / (16 * np.pi**2)) + ) + + +def _get_splitbeam_angles( + pc: xr.DataArray, + gamma_alongship, + gamma_athwartship, +) -> tuple[xr.DataArray, xr.DataArray]: + """Calculate split-beam physical angles. + + Equivalent to CRIMAC ``calcAngles``. + """ + pc_fore, pc_aft, pc_star, pc_port = _get_transducer_halves(pc) + + y_theta = pc_fore * np.conj(pc_aft) + y_phi = pc_star * np.conj(pc_port) + + theta = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_theta), np.real(y_theta)) / gamma_alongship)) + + phi = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_phi), np.real(y_phi)) / gamma_athwartship)) + + theta.name = "angle_alongship" + phi.name = "angle_athwartship" + + return theta, phi + + +def _get_transducer_halves( + pc: xr.DataArray, +) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]: + """Calculate half-transducer pulse-compressed signals. + + Equivalent to CRIMAC ``calcTransducerHalves``. + """ + pc_fore = 0.5 * (pc.isel(beam=2) + pc.isel(beam=3)) + pc_aft = 0.5 * (pc.isel(beam=0) + pc.isel(beam=1)) + pc_star = 0.5 * (pc.isel(beam=0) + pc.isel(beam=3)) + pc_port = 0.5 * (pc.isel(beam=1) + pc.isel(beam=2)) + + return pc_fore, pc_aft, pc_star, pc_port + + +def _compute_svf_power( + normalized_spectrum: np.ndarray, + n_beams: int, + z_et: float, + z_er: float, +): + """Convert normalized Sv(f) spectrum to received power spectrum. + + Equivalent to CRIMAC ``calcPowerFreqSv``. + """ + return _compute_complex_power( + normalized_spectrum=normalized_spectrum, + n_beams=n_beams, + z_et=z_et, + z_er=z_er, + ) + + +def _get_svf_frequency_grid( + f_start: float, + f_stop: float, + frequency_resolution: float, + fs_dec: float, + n_window: int, +): + """Get Sv(f) frequency grid and corresponding FFT indices.""" + + n_f_points = int(np.floor((f_stop - f_start) / frequency_resolution)) + 1 + frequency = f_start + np.arange(n_f_points) * frequency_resolution + + frequency_index = np.mod( + np.floor(frequency / fs_dec * n_window).astype(int), + n_window, + ) + + return frequency, frequency_index + + +def _compute_ts_spectrum_power( + normalized_spectrum: np.ndarray, + n_beams: int, + z_et: float, + z_er: float, +): + """Convert normalised TS(f) spectrum to received power spectrum. + + Equivalent to CRIMAC ``calcPowerFreqTS``. + """ + return _compute_complex_power( + normalized_spectrum=normalized_spectrum, + n_beams=n_beams, + z_et=z_et, + z_er=z_er, + ) + + +def _get_autocorrelation( + matched_filter: np.ndarray, + n_window: int, +): + """Get matched-filter autocorrelation spectrum. + + Equivalent to CRIMAC ``calcAutoCorrelation``. + """ + mf_auto = ( + np.convolve( + matched_filter, + np.conj(matched_filter[::-1]), + mode="full", + ) + / np.linalg.norm(matched_filter) ** 2 + ) + + mf_auto_spectrum = np.fft.fft(mf_auto, n=n_window) + + return mf_auto, mf_auto_spectrum + + +def _get_pulse_compressed_signal( + beam: xr.Dataset, + matched_filter: Dict, +) -> xr.DataArray: + """Calculate pulse-compressed complex samples for each transducer sector. + + Equivalent to CRIMAC ``calcPulseCompressedSignals``. + """ + pc = compress_pulse( + backscatter=beam["backscatter_r"] + 1j * beam["backscatter_i"], + chirp=matched_filter, + ) + pc = pc / get_norm_fac(chirp=matched_filter) + pc.name = "pulse_compressed_signal" + + return pc + + +def _get_average_signal( + signal: xr.DataArray, +) -> xr.DataArray: + """Average complex signal over transducer sectors. + + Equivalent to CRIMAC ``calcAverageSignal``. + """ + out = signal.mean(dim="beam") + out.name = "average_signal" + + return out + + +def _compute_complex_power( + normalized_spectrum: np.ndarray, + n_beams: int, + z_et: float, + z_er: float, +): + impedance_factor = (np.abs(z_er + z_et) / np.abs(z_er)) ** 2 / np.abs(z_et) + + return n_beams * (np.abs(normalized_spectrum) / (2 * np.sqrt(2))) ** 2 * impedance_factor + + def tapered_chirp( fs, transmit_duration_nominal, diff --git a/echopype/tests/calibrate/test_calibrate.py b/echopype/tests/calibrate/test_calibrate.py index c9628b67c..ff7f8b422 100644 --- a/echopype/tests/calibrate/test_calibrate.py +++ b/echopype/tests/calibrate/test_calibrate.py @@ -101,7 +101,7 @@ def test_compute_Sv_ek60_matlab(ek60_path): # Calibrate to get Sv ds_Sv = ep.calibrate.compute_Sv(echodata) - ds_TS = ep.calibrate.compute_TS(echodata) + ds_Sp = ep.calibrate.compute_Sp(echodata) # Load matlab outputs and test @@ -117,10 +117,10 @@ def check_output(da_cmp, cal_type): assert np.allclose(pyel_vals, ep_vals) # Check Sv - check_output(ds_Sv['Sv'], 'Sv') + check_output(ds_Sv["Sv"], "Sv") - # Check TS - check_output(ds_TS['TS'], 'Sp') + # Check Sp + check_output(ds_Sp["Sp"], "Sp") @pytest.mark.integration @@ -137,10 +137,10 @@ def test_compute_Sv_ek60_duplicated_freq(ek60_path): # Calibrate to get Sv ds_Sv = ep.calibrate.compute_Sv(echodata) - ds_TS = ep.calibrate.compute_TS(echodata) + ds_Sp = ep.calibrate.compute_Sp(echodata) assert isinstance(ds_Sv, xr.Dataset) - assert isinstance(ds_TS, xr.Dataset) + assert isinstance(ds_Sp, xr.Dataset) @pytest.mark.integration diff --git a/echopype/tests/calibrate/test_calibrate_ek80.py b/echopype/tests/calibrate/test_calibrate_ek80.py index 7a7cd25c7..574081b21 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80.py +++ b/echopype/tests/calibrate/test_calibrate_ek80.py @@ -299,6 +299,7 @@ def test_ek80_BB_power_from_complex( ), ], ) + def test_ek80_BB_power_compute_Sv( raw_data_path, raw_file_name, diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py new file mode 100644 index 000000000..7e2233be6 --- /dev/null +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -0,0 +1,133 @@ +import numpy as np +import pytest +import xarray as xr + +import echopype as ep + +pytestmark = pytest.mark.integration + +CRIMAC_CHANNEL = "WBT 747022-15 ES120-7CD_ES" +CRIMAC_PING_INDEX = 509 + + +@pytest.fixture(scope="module") +def ts_spectrum_example_path(test_path): + return test_path["TS_SPECTRUM_EXAMPLE"] + + +@pytest.fixture(scope="module") +def ts_raw_path(ts_spectrum_example_path): + return ts_spectrum_example_path / "IMR-D20211215-T143432-TSf.raw" + + +@pytest.fixture(scope="module") +def ts_ref(ts_spectrum_example_path): + return np.load( + ts_spectrum_example_path / "crimac_tsf_reference_outputs.npz", + allow_pickle=True, + ) + + +@pytest.fixture(scope="module") +def ts_echodata(ts_raw_path): + return ep.open_raw(ts_raw_path, sonar_model="EK80") + + +def _target_locations_from_crimac(ed, ref, channel, ping_index): + ping_time = ed["Sonar/Beam_group1"]["ping_time"].isel(ping_time=ping_index).values + + return xr.Dataset( + data_vars={ + "target_range": ("target_id", [float(ref["r_t"])]), + "angle_alongship": ("target_id", [float(ref["theta_t"])]), + "angle_athwartship": ("target_id", [float(ref["phi_t"])]), + "target_range_min": ("target_id", [float(ref["dum_r"][0])]), + "target_range_max": ("target_id", [float(ref["dum_r"][-1])]), + }, + coords={ + "target_id": [0], + "ping_time": ("target_id", [ping_time]), + "channel": ("target_id", [channel]), + }, + ) + + +def test_compute_sp_fm_complex_runs(ts_echodata): + ds = ep.calibrate.compute_Sp( + ts_echodata, + waveform_mode="BB", + encode_mode="complex", + ) + + assert "Sp" in ds + assert set(("channel", "ping_time", "range_sample")).issubset(ds["Sp"].dims) + assert np.isfinite(ds["Sp"]).any() + + +def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): + cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( + echodata=ts_echodata, + waveform_mode="BB", + encode_mode="complex", + env_params=None, + cal_params=None, + ) + + sound_speed = float(cal_obj.env_params["sound_speed"]) + + absorption_f = cal_obj._compute_absorption_f( + frequency=ts_ref["f_m"], + channel=CRIMAC_CHANNEL, + ping_idx=CRIMAC_PING_INDEX, + sound_speed=sound_speed, + ) + + np.testing.assert_allclose( + absorption_f, + ts_ref["alpha_m"], + atol=1e-8, + rtol=0.0, + ) + + +def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + ds = ep.calibrate.compute_TS_spectrum( + ts_echodata, + waveform_mode="BB", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ts = ( + ds["TS_spectrum"] + .sel(channel=CRIMAC_CHANNEL) + .isel(target_id=0) + .values + ) + + assert ts.shape == ts_ref["TS_m"].shape + np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.5, rtol=0.0) + + +def test_compute_svf_fm_complex_runs(ts_echodata): + ds = ep.calibrate.compute_Sv_f( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + frequency_resolution=1000, + range_step=0.5, + ) + + assert "Sv_f" in ds + assert set(("channel", "ping_time", "svf_range", "frequency")).issubset( + ds["Sv_f"].dims + ) + assert np.isfinite(ds["Sv_f"]).any() \ No newline at end of file diff --git a/echopype/tests/calibrate/test_ek80_complex.py b/echopype/tests/calibrate/test_ek80_complex.py index 595800496..0a671dd1b 100644 --- a/echopype/tests/calibrate/test_ek80_complex.py +++ b/echopype/tests/calibrate/test_ek80_complex.py @@ -2,7 +2,16 @@ import numpy as np import xarray as xr -from echopype.calibrate.ek80_complex import get_vend_filter_EK80 +from echopype.calibrate.ek80_complex import ( + get_vend_filter_EK80, + _get_average_signal, + _compute_power_from_complex_signal, + _compute_svf_power, + _compute_ts_spectrum_power, + _align_autocorrelation, + _compute_ts_spectrum, + _compute_ts_spectrum_calibrated, +) pytestmark = pytest.mark.unit @@ -73,3 +82,126 @@ def test_get_vend_filter_EK80(ch_num, filter_len, has_nan): assert sel_vend[var_df].values == get_vend_filter_EK80( vend, channel_id=ch, filter_name=filter_name, param_type="decimation" ) + +def test_get_average_signal(): + signal = xr.DataArray( + np.array([[1 + 1j, 3 + 3j], [5 + 5j, 7 + 7j]]), + dims=("range_sample", "beam"), + coords={"beam": [0, 1]}, + ) + + out = _get_average_signal(signal) + expected = signal.mean(dim="beam") + + xr.testing.assert_allclose(out, expected) + assert out.name == "average_signal" + + +def test_compute_power_from_complex_signal(): + signal = xr.DataArray( + np.array([[1 + 1j, 3 + 3j]]), + dims=("range_sample", "beam"), + coords={"beam": [0, 1]}, + ) + + z_et = 75.0 + z_er = 5400.0 + n_beams = signal["beam"].size + avg = signal.mean(dim="beam") + + expected = ( + n_beams + * np.abs(avg) ** 2 + / (2 * np.sqrt(2)) ** 2 + * (np.abs(z_er + z_et) / np.abs(z_er)) ** 2 + / np.abs(z_et) + ) + + out = _compute_power_from_complex_signal(signal, z_et=z_et, z_er=z_er) + + xr.testing.assert_allclose(out, expected) + assert out.name == "received_power" + + +def test_compute_svf_and_ts_spectrum_power_are_equivalent(): + normalized_spectrum = np.array([1 + 1j, 2 + 0j, 0.5 - 0.5j]) + n_beams = 4 + z_et = 75.0 + z_er = 5400.0 + + svf_power = _compute_svf_power( + normalized_spectrum=normalized_spectrum, + n_beams=n_beams, + z_et=z_et, + z_er=z_er, + ) + + ts_power = _compute_ts_spectrum_power( + normalized_spectrum=normalized_spectrum, + n_beams=n_beams, + z_et=z_et, + z_er=z_er, + ) + + np.testing.assert_allclose(svf_power, ts_power) + + +def test_align_autocorrelation(): + mf_auto = np.array([0, 0, 1, 0.5, 0.25, 0.1]) + pc_target = np.array([0.2, 1.0, 0.3]) + + out = _align_autocorrelation(mf_auto=mf_auto, pc_target=pc_target) + + np.testing.assert_allclose(out, np.array([0, 1, 0.5])) + + +def test_compute_ts_spectrum(): + pc_target = np.array([1.0, 2.0, 1.0]) + mf_auto_red = np.array([1.0, 1.0, 1.0]) + frequency = np.array([0.0, 1.0, 2.0]) + fs_dec = 8.0 + + y_pc, y_mf, y_norm = _compute_ts_spectrum( + pc_target=pc_target, + mf_auto_red=mf_auto_red, + NFFT=8, + frequency=frequency, + fs_dec=fs_dec, + ) + + assert y_pc.shape == frequency.shape + assert y_mf.shape == frequency.shape + assert y_norm.shape == frequency.shape + np.testing.assert_allclose(y_norm, y_pc / y_mf) + + +def test_compute_ts_spectrum_calibrated(): + power_spectrum = np.array([1e-12, 2e-12]) + frequency = np.array([90000.0, 100000.0]) + sound_speed = 1500.0 + target_range = 10.0 + absorption_f = np.array([0.03, 0.04]) + transmit_power = 1000.0 + gain_f = np.array([100.0, 120.0]) + + out = _compute_ts_spectrum_calibrated( + power_spectrum=power_spectrum, + target_range=target_range, + frequency=frequency, + sound_speed=sound_speed, + absorption_f=absorption_f, + transmit_power=transmit_power, + gain_f=gain_f, + ) + + wavelength = sound_speed / frequency + expected = ( + 10 * np.log10(power_spectrum) + + 40 * np.log10(target_range) + + 2 * absorption_f * target_range + - 10 * np.log10( + transmit_power * wavelength**2 * gain_f**2 / (16 * np.pi**2) + ) + ) + + np.testing.assert_allclose(out, expected) \ No newline at end of file diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 7a09ad711..238165b9f 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -36,7 +36,7 @@ "ek80_bb_complex_multiplex.zip", "ek80_bb_with_calibration.zip", "ek80_duplicate_ping_times.zip", "ek80_ext.zip", "ek80_invalid_env_datagrams.zip", "ek80_missing_sound_velocity_profile.zip", "ek80_new.zip", "ek80_sequence.zip", - "es60.zip", "es70.zip", "es80.zip", "legacy_datatree.zip", + "es60.zip", "es70.zip", "es80.zip", "legacy_datatree.zip", "ts_spectrum_example_data.zip", ] # v0.11.1a2 checksums (GitHub release assets) @@ -62,6 +62,7 @@ "es70.zip": "sha256:a6b4f27f33f09bace26264de6984fdb4111a3a0337bc350c3c1d25c8b3effc7c", "es80.zip": "sha256:b37ee01462f46efe055702c20be67d2b8c6b786844b183b16ffc249c7c5ec704", "legacy_datatree.zip": "sha256:820cd252047dbf35fa5fb04a9aafee7f7659e0fe4f7d421d69901c57deb6c9d5", # noqa: E501 + "ts_spectrum_example_data.zip": "sha256:a35a1108dcfafae725d826ffa23ea2e111fdb420e7ce433e49e1797381bb6763", } EP = pooch.create( @@ -128,6 +129,15 @@ def _unpack(fname, action, pooch_instance): return str(out) for b in bundles: + out = TEST_DATA_FOLDER / Path(b).stem + + zip_path = TEST_DATA_FOLDER / b + + if zip_path.exists() and out.exists() and any(out.iterdir()): + print(f"[echopype-ci] using cached bundle: {b}") + print(f"[echopype-ci] -> {out}") + continue + url = base.format(version=ver) + b print(f"[echopype-ci] fetching bundle: {b}") print(f"[echopype-ci] -> URL: {url}") @@ -177,6 +187,7 @@ def test_path(): "EK80_MULTI": TEST_DATA_FOLDER / "ek80_bb_complex_multiplex", "ECS": TEST_DATA_FOLDER / "ecs", "LEGACY_DATATREE": TEST_DATA_FOLDER / "legacy_datatree", + "TS_SPECTRUM_EXAMPLE": TEST_DATA_FOLDER / "ts_spectrum_example_data", } From 5f05da161762939b01b3ef2abd0523ae74852219 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sat, 6 Jun 2026 19:04:14 -0700 Subject: [PATCH 04/25] add tqdm adding tqdm, most likely temporarily --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 9f4b1acf7..4c27a2d78 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ dependencies = [ "typing-extensions", "xarray>=2026.01.0", "zarr>=3", + "tqdm", ] [project.urls] From fe4adec7fc4d4bff247c4a333a09def5e2154623 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sat, 6 Jun 2026 19:21:32 -0700 Subject: [PATCH 05/25] Update conftest.py --- echopype/tests/conftest.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 238165b9f..20b3e69ad 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -129,15 +129,6 @@ def _unpack(fname, action, pooch_instance): return str(out) for b in bundles: - out = TEST_DATA_FOLDER / Path(b).stem - - zip_path = TEST_DATA_FOLDER / b - - if zip_path.exists() and out.exists() and any(out.iterdir()): - print(f"[echopype-ci] using cached bundle: {b}") - print(f"[echopype-ci] -> {out}") - continue - url = base.format(version=ver) + b print(f"[echopype-ci] fetching bundle: {b}") print(f"[echopype-ci] -> URL: {url}") From 47ad17d69264a6fc666931b7b7526e871633b132 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sat, 6 Jun 2026 19:43:33 -0700 Subject: [PATCH 06/25] change org name --- echopype/tests/conftest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 20b3e69ad..ffa2e59c6 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -17,7 +17,7 @@ # Lock to the known-good assets release (can be overridden via env if needed) base = os.getenv( "ECHOPYPE_DATA_BASEURL", - "https://github.com/OSOceanAcoustics/echopype/releases/download/{version}/", + "https://github.com/echostack-org/echopype/releases/download/{version}/", ) cache_dir = pooch.os_cache("echopype") From b2ff5a2f5be8191de5e175fceadc0c654d8e7bae Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sun, 7 Jun 2026 10:02:14 -0700 Subject: [PATCH 07/25] Update conftest.py try to fix error occuring only in ubuntu job --- echopype/tests/conftest.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index ffa2e59c6..4cb3e2681 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -96,7 +96,9 @@ def _unpack(fname, action, pooch_instance): time.sleep(1) with ZipFile(z, "r") as f: - f.extractall(out) + for member in f.infolist(): + member.filename = member.filename.replace("\\", "/") + f.extract(member, out) # flatten single nested dir if needed try: From 31d401bb4f769b0459b8f565cd22265d87bd6e57 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 8 Jun 2026 20:45:14 -0700 Subject: [PATCH 08/25] Fix formatting in test_calibrate_ek80.py --- echopype/tests/calibrate/test_calibrate_ek80.py | 1 - 1 file changed, 1 deletion(-) diff --git a/echopype/tests/calibrate/test_calibrate_ek80.py b/echopype/tests/calibrate/test_calibrate_ek80.py index 574081b21..7a7cd25c7 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80.py +++ b/echopype/tests/calibrate/test_calibrate_ek80.py @@ -299,7 +299,6 @@ def test_ek80_BB_power_from_complex( ), ], ) - def test_ek80_BB_power_compute_Sv( raw_data_path, raw_file_name, From c8c657a16e2e7da4c1c30afa06bfb3d94ed69093 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Wed, 10 Jun 2026 15:09:53 -0700 Subject: [PATCH 09/25] removing _get_splitbeam_angles and _get_transducer_halves --- echopype/calibrate/calibrate_ek.py | 26 +++++++++----------- echopype/calibrate/ek80_complex.py | 39 ------------------------------ 2 files changed, 11 insertions(+), 54 deletions(-) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 9ce328c61..d4db16032 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -788,18 +788,6 @@ def _get_beam_compensated_gain( self.vend["gain"].sel(cal_channel_id=channel).values, ) - angle_offset_alongship = np.interp( - frequency, - self.vend["cal_frequency"].values, - self.vend["angle_offset_alongship"].sel(cal_channel_id=channel).values, - ) - - angle_offset_athwartship = np.interp( - frequency, - self.vend["cal_frequency"].values, - self.vend["angle_offset_athwartship"].sel(cal_channel_id=channel).values, - ) - beamwidth_alongship = np.interp( frequency, self.vend["cal_frequency"].values, @@ -812,9 +800,17 @@ def _get_beam_compensated_gain( self.vend["beamwidth_athwartship"].sel(cal_channel_id=channel).values, ) - fac_along = (np.abs(theta - angle_offset_alongship) / (beamwidth_alongship / 2)) ** 2 - - fac_athwart = (np.abs(phi - angle_offset_athwartship) / (beamwidth_athwartship / 2)) ** 2 + # NOTE: + # theta and phi are expected to be offset-corrected split-beam angles + # (e.g. from get_angle_complex_samples() / add_splitbeam_angle()). + # Therefore angle offsets are not applied here again. + # + # This differs from the CRIMAC reference implementation, where the + # angles returned by calcAngles() do not include the angle-offset + # correction. Applying the offset correction to the CRIMAC angles + # yields close agreement with the echopype convention. + fac_along = (np.abs(theta) / (beamwidth_alongship / 2)) ** 2 + fac_athwart = (np.abs(phi) / (beamwidth_athwartship / 2)) ** 2 beam_correction_db = ( 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py index bd259ca60..4cdf96c3f 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -158,45 +158,6 @@ def _compute_ts_spectrum_calibrated( ) -def _get_splitbeam_angles( - pc: xr.DataArray, - gamma_alongship, - gamma_athwartship, -) -> tuple[xr.DataArray, xr.DataArray]: - """Calculate split-beam physical angles. - - Equivalent to CRIMAC ``calcAngles``. - """ - pc_fore, pc_aft, pc_star, pc_port = _get_transducer_halves(pc) - - y_theta = pc_fore * np.conj(pc_aft) - y_phi = pc_star * np.conj(pc_port) - - theta = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_theta), np.real(y_theta)) / gamma_alongship)) - - phi = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_phi), np.real(y_phi)) / gamma_athwartship)) - - theta.name = "angle_alongship" - phi.name = "angle_athwartship" - - return theta, phi - - -def _get_transducer_halves( - pc: xr.DataArray, -) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]: - """Calculate half-transducer pulse-compressed signals. - - Equivalent to CRIMAC ``calcTransducerHalves``. - """ - pc_fore = 0.5 * (pc.isel(beam=2) + pc.isel(beam=3)) - pc_aft = 0.5 * (pc.isel(beam=0) + pc.isel(beam=1)) - pc_star = 0.5 * (pc.isel(beam=0) + pc.isel(beam=3)) - pc_port = 0.5 * (pc.isel(beam=1) + pc.isel(beam=2)) - - return pc_fore, pc_aft, pc_star, pc_port - - def _compute_svf_power( normalized_spectrum: np.ndarray, n_beams: int, From 2b9db6b7f2a583edd97cefd3b63c8e32d891a85b Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 15 Jun 2026 22:08:13 -0700 Subject: [PATCH 10/25] remove Sv(f) API and refactor TS(f) methds --- echopype/calibrate/__init__.py | 4 +- echopype/calibrate/api.py | 57 +++- echopype/calibrate/calibrate_ek.py | 420 +++++++++-------------------- echopype/calibrate/ek80_complex.py | 147 ++++------ 4 files changed, 232 insertions(+), 396 deletions(-) diff --git a/echopype/calibrate/__init__.py b/echopype/calibrate/__init__.py index 7091c3504..8fb4e732e 100644 --- a/echopype/calibrate/__init__.py +++ b/echopype/calibrate/__init__.py @@ -1,3 +1,3 @@ -from .api import compute_Sp, compute_Sv, compute_Sv_f, compute_TS, compute_TS_spectrum +from .api import compute_Sp, compute_Sv, compute_Sv_spectrum, compute_TS, compute_TS_spectrum -__all__ = ["compute_Sv", "compute_TS", "compute_Sp", "compute_Sv_f", "compute_TS_spectrum"] +__all__ = ["compute_Sv", "compute_TS", "compute_Sp", "compute_Sv_spectrum", "compute_TS_spectrum"] diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index 4ed4582d4..2ffd87e4f 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -82,10 +82,10 @@ def _compute_cal_ds(echodata, slice_dict): cal_obj._check_echodata_backscatter_size() compute_methods = { - "Sv": "compute_Sv", "Sp": "compute_Sp", "TS": "compute_TS", - "Sv_f": "compute_Sv_f", + "Sv": "compute_Sv", + # add Sp_spectrum?? "TS_spectrum": "compute_TS_spectrum", } @@ -176,7 +176,8 @@ def _compute_cal_ds(echodata, slice_dict): # Calibrate and drop filter_time cal_ds_iteration = _compute_cal_ds(echodata, slice_dict) - cal_ds_list.append(cal_ds_iteration.drop_vars("filter_time")) + if "filter_time" in cal_ds_iteration: + cal_ds_iteration = cal_ds_iteration.drop_vars("filter_time") # # Alternative? # for channel in echodata[ed_beam_group]["channel"].values: @@ -239,10 +240,9 @@ def _add_attrs(cal_type, ds): ds[cal_type].attrs = { "long_name": { - "Sv": "Volume backscattering strength (Sv re 1 m-1)", "Sp": "Point scattering strength (Sp re 1 m^2)", "TS": "Target strength (TS re 1 m^2)", - "Sv_f": "Frequency-dependent volume backscattering strength (Sv(f) re 1 m-1)", + "Sv": "Volume backscattering strength (Sv re 1 m-1)", "TS_spectrum": "Frequency-dependent target strength spectrum (TS(f) re 1 m^2)", }[cal_type], "units": "dB", @@ -382,13 +382,16 @@ def compute_Sv(echodata: EchoData, **kwargs) -> xr.Dataset: return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs) -def compute_Sv_f(echodata: EchoData, **kwargs) -> xr.Dataset: +def compute_Sv_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: """ Compute frequency-dependent volume backscattering strength Sv(f) from broadband EK80 complex data. - TODO: add more details here when the function is fully implemented + + Notes + ----- + This functionality is not yet implemented. """ - return _compute_cal(cal_type="Sv_f", echodata=echodata, **kwargs) + raise NotImplementedError("compute_Sv_spectrum is not yet implemented.") def compute_Sp(echodata: EchoData, **kwargs) -> xr.Dataset: @@ -508,7 +511,41 @@ def compute_TS(echodata: EchoData, **kwargs): def compute_TS_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: """ - Compute broadband frequency-dependent target strength spectrum. - TODO: add more details here when the function is fully implemented + Compute broadband frequency-dependent target strength spectrum, TS(f), + from EK80 broadband/FM complex data. + + Parameters + ---------- + point_locations : xr.Dataset + Locations of targets for which TS(f) should be computed. + Must contain ``channel``, ``ping_time``, and ``target_range`` for each + ``target_id``. If ``target_range_min`` and ``target_range_max`` are + provided, they define the target echo segment. Otherwise, the segment + is built around ``target_range`` using ``NFFT`` and ``split_front``. + + NFFT : int, optional + Number of FFT points used to compute the target spectrum. If not + provided, a value is inferred from the output frequency grid. + + n_f_points : int, optional + Number of frequency points in the output TS(f) spectrum. Used when + ``frequency_resolution`` is not provided. + + split_front : float, default 0.25 + Fraction of the FFT window placed before the target location when only + ``target_range`` is provided. + + window : str, tuple, float or None, default None + Window passed directly to ``scipy.signal.get_window``. If ``None``, + a rectangular/boxcar window is used. + + frequency_resolution : float, optional + Desired spacing of the output frequency grid in Hz. Used to define the + frequency grid on which TS(f) is evaluated. + + Returns + ------- + xr.Dataset + Dataset containing frequency-dependent target strength, TS(f). """ return _compute_cal(cal_type="TS_spectrum", echodata=echodata, **kwargs) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index d4db16032..6f76e2860 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -2,6 +2,7 @@ import numpy as np import xarray as xr +from scipy.signal import get_window from tqdm.auto import tqdm from ..echodata import EchoData @@ -14,16 +15,13 @@ from .ek80_complex import ( _align_autocorrelation, _compute_power_from_complex_signal, - _compute_svf_power, - _compute_svf_spectrum, _compute_ts_spectrum, _compute_ts_spectrum_calibrated, _compute_ts_spectrum_power, _get_autocorrelation, _get_average_signal, - _get_hanning_window, _get_pulse_compressed_signal, - _get_svf_frequency_grid, + _get_splitbeam_angles, get_filter_coeff, get_tau_effective, get_transmit_signal, @@ -201,20 +199,16 @@ def _cal_power_samples(self, cal_type: str) -> xr.Dataset: beamwidth_alongship = self.cal_params["beamwidth_alongship"] beamwidth_athwartship = self.cal_params["beamwidth_athwartship"] - fac_along = ( - np.abs(angle_alongship - angle_offset_alongship) / (beamwidth_alongship / 2) - ) ** 2 - - fac_athwart = ( - np.abs(angle_athwartship - angle_offset_athwartship) - / (beamwidth_athwartship / 2) - ) ** 2 - - beam_correction = ( - 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) + beam_correction_db = self._get_beam_correction( + theta=angle_alongship, + phi=angle_athwartship, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, ) - out = sp + beam_correction + out = sp + beam_correction_db else: out = sp @@ -663,17 +657,43 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: # remove TS below and create a TS that build on Sp elif cal_type in ("Sp", "TS"): range_safe = self._safe_range_for_log(range_meter) - spreading_loss_safe = 20 * np.log10(range_safe) - out = ( + sp = ( 10 * np.log10(prx) + 2 * spreading_loss_safe + absorption_loss - 10 * np.log10(wavelength**2 * transmit_power / (16 * np.pi**2)) - 2 * gain ) - out.name = cal_type # change + if cal_type == "TS" and self.waveform_mode == "CW": + bs = self.beam["backscatter_r"] + 1j * self.beam["backscatter_i"] + + angle_alongship, angle_athwartship = _get_splitbeam_angles( + pc=bs, + gamma_alongship=self.cal_params["angle_sensitivity_alongship"], + gamma_athwartship=self.cal_params["angle_sensitivity_athwartship"], + ) + + angle_offset_alongship = self.cal_params["angle_offset_alongship"] + angle_offset_athwartship = self.cal_params["angle_offset_athwartship"] + beamwidth_alongship = self.cal_params["beamwidth_alongship"] + beamwidth_athwartship = self.cal_params["beamwidth_athwartship"] + + beam_correction_db = self._get_beam_correction( + theta=angle_alongship, + phi=angle_athwartship, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, + ) + + out = sp + beam_correction_db + else: + out = sp + + out.name = cal_type # Attach calculated range (with units meter) into data set out = out.to_dataset().merge(range_meter) @@ -696,8 +716,6 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: return out - ### Common helpers for complex sample calibration and Sv(f)/TS(f) calibration - def _select_param( self, da: xr.DataArray, @@ -739,37 +757,21 @@ def _compute_absorption_f( formula_source="FG", ) - ### Helpers for Svf - - def _compute_svf( + def _get_beam_correction( self, - power_spectrum: np.ndarray, - frequency: np.ndarray, - svf_range: np.ndarray, - sound_speed: float, - absorption_f: np.ndarray, - transmit_power: float, - psi_f: np.ndarray, - gain_f: np.ndarray, - window_duration: float, + theta, + phi, + angle_offset_alongship, + angle_offset_athwartship, + beamwidth_alongship, + beamwidth_athwartship, ): - """Apply CRIMAC/pyEcholab-style Sv(f) calibration equation. - - Equivalent to CRIMAC ``calcSvf``. - """ - wavelength_f = sound_speed / frequency + """Compute empirical split-beam beam correction in dB.""" + fac_along = (np.abs(theta - angle_offset_alongship) / (beamwidth_alongship / 2)) ** 2 - G = ( - transmit_power * wavelength_f**2 * sound_speed * window_duration * psi_f * gain_f**2 - ) / (32 * np.pi**2) + fac_athwart = (np.abs(phi - angle_offset_athwartship) / (beamwidth_athwartship / 2)) ** 2 - return ( - 10 * np.log10(power_spectrum) - + 2 * absorption_f[None, :] * svf_range[:, None] - - 10 * np.log10(G[None, :]) - ) - - ### Helpers for TSf + return 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) def _get_beam_compensated_gain( self, @@ -781,6 +783,22 @@ def _get_beam_compensated_gain( """Calculate beam-compensated gain. Equivalent to CRIMAC ``calc_g(theta, phi, f)``. + + # NOTE: + # For TS(f), theta and phi are expected to be CRIMAC-style split-beam + # angles computed directly from the pulse-compressed sector signals. + # These angles include the sensitivity conversion but do not include + # the angle-offset correction. + # + # The frequency-dependent beam compensation is therefore applied here, + # following CRIMAC ``calc_g(theta, phi, f)``: gain, beamwidths, and + # angle offsets are interpolated onto the TS(f) frequency grid, and the + # target angular distance from the beam axis is evaluated independently + # at each frequency. + # + # This differs from ``add_splitbeam_angle()``, which returns conventional + # echopype mechanical angles with the angle offsets already removed. + """ gain_db = np.interp( frequency, @@ -788,6 +806,18 @@ def _get_beam_compensated_gain( self.vend["gain"].sel(cal_channel_id=channel).values, ) + angle_offset_alongship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["angle_offset_alongship"].sel(cal_channel_id=channel).values, + ) + + angle_offset_athwartship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["angle_offset_athwartship"].sel(cal_channel_id=channel).values, + ) + beamwidth_alongship = np.interp( frequency, self.vend["cal_frequency"].values, @@ -800,204 +830,17 @@ def _get_beam_compensated_gain( self.vend["beamwidth_athwartship"].sel(cal_channel_id=channel).values, ) - # NOTE: - # theta and phi are expected to be offset-corrected split-beam angles - # (e.g. from get_angle_complex_samples() / add_splitbeam_angle()). - # Therefore angle offsets are not applied here again. - # - # This differs from the CRIMAC reference implementation, where the - # angles returned by calcAngles() do not include the angle-offset - # correction. Applying the offset correction to the CRIMAC angles - # yields close agreement with the echopype convention. - fac_along = (np.abs(theta) / (beamwidth_alongship / 2)) ** 2 - fac_athwart = (np.abs(phi) / (beamwidth_athwartship / 2)) ** 2 - - beam_correction_db = ( - 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) + beam_correction_db = self._get_beam_correction( + theta=theta, + phi=phi, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, ) return 10 ** ((gain_db - beam_correction_db) / 10) - ### Sv(f)/TS(f) calibration pipelines - - def _cal_complex_samples_svf( - self, - pc: xr.DataArray, - matched_filter: Dict, - frequency_resolution: float = 1000, - range_step: float = 0.5, - ) -> xr.Dataset: - """Compute frequency-dependent volume backscattering strength Sv(f).""" - - pc_avg = _get_average_signal(pc) - - out_by_channel = [] - - for channel in self.beam["channel"].values: - if channel not in self.vend["cal_channel_id"].values: - continue - - pc_avg_ch = pc_avg.sel(channel=channel) - range_ch = self.range_meter.sel(channel=channel) - - svf_list = [] - svf_ping_time = [] - frequency_ref = None - svf_range_ref = None - - for ping_idx in tqdm(range(pc_avg_ch.sizes["ping_time"]), desc=f"{channel}"): - pc_avg_1d = pc_avg_ch.isel(ping_time=ping_idx) - range_1d = range_ch.isel(ping_time=ping_idx) - - valid = np.isfinite(pc_avg_1d) & np.isfinite(range_1d) - - if not bool(valid.any()): - continue - - pc_avg_1d = pc_avg_1d.where(valid, drop=True).values - range_1d = range_1d.where(valid, drop=True).values - - sound_speed = float( - self._select_param(self.env_params["sound_speed"], channel, ping_idx) - ) - transmit_power = float( - self._select_param(self.beam["transmit_power"], channel, ping_idx) - ) - tau = float( - self._select_param(self.beam["transmit_duration_nominal"], channel, ping_idx) - ) - sample_interval = float( - self._select_param(self.beam["sample_interval"], channel, ping_idx) - ) - f_start = float( - self._select_param(self.beam["transmit_frequency_start"], channel, ping_idx) - ) - f_stop = float( - self._select_param(self.beam["transmit_frequency_stop"], channel, ping_idx) - ) - equivalent_beam_angle = float( - self._select_param(self.cal_params["equivalent_beam_angle"], channel, ping_idx) - ) - z_et = float( - self._select_param(self.cal_params["impedance_transducer"], channel, ping_idx) - ) - z_er = float( - self._select_param(self.cal_params["impedance_transceiver"], channel, None) - ) - - window, n_window, window_duration, fs_dec = _get_hanning_window( - sound_speed=sound_speed, - tau=tau, - sample_interval=sample_interval, - ) - - frequency, frequency_index = _get_svf_frequency_grid( - f_start=f_start, - f_stop=f_stop, - frequency_resolution=frequency_resolution, - fs_dec=fs_dec, - n_window=n_window, - ) - - _, mf_auto_spectrum = _get_autocorrelation( - matched_filter=matched_filter[channel], - n_window=n_window, - ) - - step = max( - 1, - int(np.round(range_step / np.nanmedian(np.diff(range_1d)))), - ) - - pc_spread = pc_avg_1d * range_1d - - normalized_spectrum, svf_range = _compute_svf_spectrum( - pc_spread=pc_spread, - range_meter=range_1d, - window=window, - n_window=n_window, - frequency_index=frequency_index, - mf_auto_spectrum=mf_auto_spectrum, - step=step, - ) - - if normalized_spectrum.size == 0: - continue - - power_spectrum = _compute_svf_power( - normalized_spectrum=normalized_spectrum, - n_beams=self.beam["beam"].size, - z_et=z_et, - z_er=z_er, - ) - - gain_db = np.interp( - frequency, - self.vend["cal_frequency"].values, - self.vend["gain"].sel(cal_channel_id=channel).values, - ) - gain_f = 10 ** (gain_db / 10) - - frequency_nominal = float( - self._select_param(self.beam["frequency_nominal"], channel, None) - ) - psi_f = 10 ** (equivalent_beam_angle / 10) * (frequency_nominal / frequency) ** 2 - - absorption_f = self._compute_absorption_f( - frequency=frequency, - channel=channel, - ping_idx=ping_idx, - sound_speed=sound_speed, - ) - - svf = self._compute_svf( - power_spectrum=power_spectrum, - frequency=frequency, - svf_range=svf_range, - sound_speed=sound_speed, - absorption_f=absorption_f, - transmit_power=transmit_power, - psi_f=psi_f, - gain_f=gain_f, - window_duration=window_duration, - ) - - if frequency_ref is None: - frequency_ref = frequency - svf_range_ref = svf_range - elif ( - frequency.shape != frequency_ref.shape or svf_range.shape != svf_range_ref.shape - ): - continue - - svf_list.append(svf) - svf_ping_time.append(self.beam["ping_time"].isel(ping_time=ping_idx).values) - - if not svf_list: - continue - - ds_ch = xr.Dataset( - { - "Sv_f": ( - ["channel", "ping_time", "svf_range", "frequency"], - np.asarray(svf_list)[None, :, :, :], - ) - }, - coords={ - "channel": [channel], - "ping_time": svf_ping_time, - "svf_range": svf_range_ref, - "frequency": frequency_ref, - }, - ) - - out_by_channel.append(ds_ch) - - if not out_by_channel: - raise ValueError("No valid Sv(f) pings produced.") - - return xr.concat(out_by_channel, dim="channel") - def _cal_complex_samples_TS_spectrum( self, pc: xr.DataArray, @@ -1048,8 +891,6 @@ def _cal_complex_samples_TS_spectrum( point_ping_time = points_ch["ping_time"].isel(target_id=point_idx).values target_range = float(points_ch["target_range"].isel(target_id=point_idx)) - theta_t = float(points_ch["angle_alongship"].isel(target_id=point_idx)) - phi_t = float(points_ch["angle_athwartship"].isel(target_id=point_idx)) target_id = points_ch["target_id"].isel(target_id=point_idx).values ping_idx = int( @@ -1116,8 +957,6 @@ def _cal_complex_samples_TS_spectrum( sound_speed=sound_speed, ) - ########## try to fix something ??????? - target_range_min = float(points_ch["target_range_min"].isel(target_id=point_idx)) target_range_max = float(points_ch["target_range_max"].isel(target_id=point_idx)) @@ -1128,19 +967,41 @@ def _cal_complex_samples_TS_spectrum( pc_target = pc_avg_1d[target_mask] - if window == "hann": - win = np.hanning(pc_target.size) - elif window == "hamming": - win = np.hamming(pc_target.size) - elif window in (None, "boxcar", "rectangular"): + gamma_alongship = float( + self._select_param( + self.cal_params["angle_sensitivity_alongship"], + channel, + ping_idx, + ) + ) + gamma_athwartship = float( + self._select_param( + self.cal_params["angle_sensitivity_athwartship"], + channel, + ping_idx, + ) + ) + + theta_raw_da, phi_raw_da = _get_splitbeam_angles( + pc=pc.sel(channel=channel).isel(ping_time=ping_idx), + gamma_alongship=gamma_alongship, + gamma_athwartship=gamma_athwartship, + ) + + theta_raw = theta_raw_da.where(valid, drop=True).values + phi_raw = phi_raw_da.where(valid, drop=True).values + + idx_peak = int(np.nanargmax(np.abs(pc_avg_1d[target_mask]) ** 2)) + theta_t = float(theta_raw[target_mask][idx_peak]) + phi_t = float(phi_raw[target_mask][idx_peak]) + + if window is None: win = np.ones(pc_target.size) else: - raise ValueError(f"Unsupported window: {window}") + win = get_window(window, pc_target.size) pc_target = pc_target * win - ########## - mf_auto, _ = _get_autocorrelation( matched_filter=matched_filter[channel], n_window=n_fft, @@ -1242,13 +1103,10 @@ def _cal_complex_samples_TS_spectrum( return xr.concat(out_by_channel, dim="channel") - ### Dispatchers - def _cal_complex_samples_f( self, cal_type: str, frequency_resolution: float | None = None, - range_step: float = 0.5, point_locations: xr.Dataset | None = None, NFFT: int | None = None, n_f_points: int | None = None, @@ -1275,14 +1133,6 @@ def _cal_complex_samples_f( matched_filter=matched_filter, ) - if cal_type == "Sv_f": - return self._cal_complex_samples_svf( - pc=pc, - matched_filter=matched_filter, - frequency_resolution=frequency_resolution, - range_step=range_step, - ) - if cal_type == "TS_spectrum": return self._cal_complex_samples_TS_spectrum( pc=pc, @@ -1322,7 +1172,7 @@ def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: ---------- cal_type : str Calibration type. Supported values include ``"Sv"``, ``"Sp"``, - ``"TS"``, ``"Sv_f"``, and ``"TS_spectrum"``. + ``"TS"``, and ``"TS_spectrum"``. Returns ------- @@ -1334,7 +1184,7 @@ def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: True if self.waveform_mode == "BB" or self.encode_mode == "complex" else False ) - if cal_type in ("Sv_f", "TS_spectrum"): + if cal_type in ("TS_spectrum",): # Frequency-dependent broadband calibration: keep pulse-compressed # signal for FFT-based Sv(f) or TS(f) processing. ds_cal = self._cal_complex_samples_f(cal_type=cal_type, **kwargs) @@ -1364,12 +1214,13 @@ def compute_Sv(self): return self._compute_cal(cal_type="Sv") def compute_Sp(self): - """Compute point scattering strength Sp. + """ + Compute point scattering strength (Sp) from raw data. - For CW data, this is computed from the received power samples on the - range grid. For FM complex data, this is computed after pulse - compression and represents a band-averaged point-scattering-strength - echogram. + For CW data, Sp is computed on the range grid from power samples + or from complex samples converted to received power. For EK80 + broadband/FM complex data, Sp is computed after pulse compression + and represents a band-averaged point-scattering-strength echogram. Returns ------- @@ -1380,7 +1231,12 @@ def compute_Sp(self): return self._compute_cal(cal_type="Sp") def compute_TS(self): - """Compute target strength (TS). + """ + For CW split-beam data, TS is computed from Sp with beam compensation + using split-beam angle information. For EK80 broadband/FM complex data, + this function returns a band-averaged gridded product after pulse + compression. For frequency-dependent target spectra at specified target + locations, use ``compute_TS_spectrum``. Returns ------- @@ -1390,24 +1246,6 @@ def compute_TS(self): """ return self._compute_cal(cal_type="TS") - def compute_Sv_f( - self, - frequency_resolution: float = 1000, - range_step: float = 0.5, - ): - """Compute frequency-dependent volume backscattering strength Sv(f). - - Returns - ------- - Sv_f : xr.Dataset - A Dataset containing frequency-dependent volume backscattering strength. - """ - return self._compute_cal( - cal_type="Sv_f", - frequency_resolution=frequency_resolution, - range_step=range_step, - ) - def compute_TS_spectrum( self, point_locations: xr.Dataset, diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py index 4cdf96c3f..1ba65a848 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -9,6 +9,60 @@ from ..convert.set_groups_ek80 import DECIMATION, FILTER_IMAG, FILTER_REAL +def _get_transducer_halves( + pc: xr.DataArray, +) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]: + """Calculate half-transducer pulse-compressed signals. + + Equivalent to CRIMAC ``calcTransducerHalves`` for 4-sector transducers. + """ + if pc.sizes["beam"] != 4: + raise NotImplementedError( + "Transducer halves are only defined for 4-sector split-beam data." + ) + + pc_fore = 0.5 * (pc.isel(beam=2) + pc.isel(beam=3)) + pc_aft = 0.5 * (pc.isel(beam=0) + pc.isel(beam=1)) + pc_star = 0.5 * (pc.isel(beam=0) + pc.isel(beam=3)) + pc_port = 0.5 * (pc.isel(beam=1) + pc.isel(beam=2)) + + return pc_fore, pc_aft, pc_star, pc_port + + +def _get_splitbeam_angles( + pc: xr.DataArray, + gamma_alongship, + gamma_athwartship, +) -> tuple[xr.DataArray, xr.DataArray]: + """Calculate raw split-beam physical angles before angle-offset correction. + + For 4-sector data this follows CRIMAC ``calcAngles``. For 3-sector data, + the sector geometry follows the same convention used by ``add_splitbeam_angle``. + Angle offsets are not applied here because TS(f) beam compensation applies + frequency-dependent offsets later. + """ + if pc.sizes["beam"] == 4: + pc_fore, pc_aft, pc_star, pc_port = _get_transducer_halves(pc) + + y_theta = pc_fore * np.conj(pc_aft) + y_phi = pc_star * np.conj(pc_port) + + theta = np.rad2deg( + np.arcsin(np.arctan2(np.imag(y_theta), np.real(y_theta)) / gamma_alongship) + ) + + phi = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_phi), np.real(y_phi)) / gamma_athwartship)) + else: + raise NotImplementedError( + f"Split-beam angle calculation is not implemented for {pc.sizes['beam']} sectors." + ) + + theta.name = "angle_alongship" + phi.name = "angle_athwartship" + + return theta, phi + + def _compute_power_from_complex_signal( signal: xr.DataArray, z_et, @@ -80,61 +134,6 @@ def _compute_ts_spectrum( return pc_target_spectrum, mf_auto_red_spectrum, normalized_spectrum -def _get_hanning_window( - sound_speed: float, - tau: float, - sample_interval: float, -): - """Get normalized Hann window for Sv(f) sliding-window FFT. - - Equivalent to CRIMAC ``defHanningWindow``. - """ - dr = sample_interval * sound_speed / 2 - fs_dec = 1 / sample_interval - - L = sound_speed * 2 * tau / dr - n_window = int(2 ** np.ceil(np.log2(L))) - - window = np.hanning(n_window) - window = window / (np.linalg.norm(window) / np.sqrt(n_window)) - - window_duration = n_window / fs_dec - - return window, n_window, window_duration, fs_dec - - -def _compute_svf_spectrum( - pc_spread: np.ndarray, - range_meter: np.ndarray, - window: np.ndarray, - n_window: int, - frequency_index: np.ndarray, - mf_auto_spectrum: np.ndarray, - step: int = 1, -): - """Compute normalized sliding-window spectrum for Sv(f). - - Equivalent to CRIMAC ``calcDFTforSv``. - """ - mf_auto_spectrum_f = mf_auto_spectrum[frequency_index] - - normalized_spectrum = [] - svf_range = [] - - for start in range(0, len(pc_spread) - n_window, step): - stop = start + n_window - - windowed_signal = window * pc_spread[start:stop] - signal_spectrum = np.fft.fft(windowed_signal, n=n_window)[frequency_index] - - normalized_spectrum.append(signal_spectrum / mf_auto_spectrum_f) - - center = (start + stop) // 2 - svf_range.append(range_meter[center]) - - return np.asarray(normalized_spectrum), np.asarray(svf_range) - - def _compute_ts_spectrum_calibrated( power_spectrum: np.ndarray, target_range: float, @@ -158,44 +157,6 @@ def _compute_ts_spectrum_calibrated( ) -def _compute_svf_power( - normalized_spectrum: np.ndarray, - n_beams: int, - z_et: float, - z_er: float, -): - """Convert normalized Sv(f) spectrum to received power spectrum. - - Equivalent to CRIMAC ``calcPowerFreqSv``. - """ - return _compute_complex_power( - normalized_spectrum=normalized_spectrum, - n_beams=n_beams, - z_et=z_et, - z_er=z_er, - ) - - -def _get_svf_frequency_grid( - f_start: float, - f_stop: float, - frequency_resolution: float, - fs_dec: float, - n_window: int, -): - """Get Sv(f) frequency grid and corresponding FFT indices.""" - - n_f_points = int(np.floor((f_stop - f_start) / frequency_resolution)) + 1 - frequency = f_start + np.arange(n_f_points) * frequency_resolution - - frequency_index = np.mod( - np.floor(frequency / fs_dec * n_window).astype(int), - n_window, - ) - - return frequency, frequency_index - - def _compute_ts_spectrum_power( normalized_spectrum: np.ndarray, n_beams: int, From 7e845bc70ea7885217247a78e94466e38e84ff5f Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 15 Jun 2026 22:21:31 -0700 Subject: [PATCH 11/25] removing Sv_f tests --- echopype/calibrate/calibrate_ek.py | 1 - echopype/tests/calibrate/test_ek80_complex.py | 24 ------------------- 2 files changed, 25 deletions(-) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 6f76e2860..8c019a64e 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -768,7 +768,6 @@ def _get_beam_correction( ): """Compute empirical split-beam beam correction in dB.""" fac_along = (np.abs(theta - angle_offset_alongship) / (beamwidth_alongship / 2)) ** 2 - fac_athwart = (np.abs(phi - angle_offset_athwartship) / (beamwidth_athwartship / 2)) ** 2 return 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) diff --git a/echopype/tests/calibrate/test_ek80_complex.py b/echopype/tests/calibrate/test_ek80_complex.py index 0a671dd1b..e5bf54853 100644 --- a/echopype/tests/calibrate/test_ek80_complex.py +++ b/echopype/tests/calibrate/test_ek80_complex.py @@ -6,7 +6,6 @@ get_vend_filter_EK80, _get_average_signal, _compute_power_from_complex_signal, - _compute_svf_power, _compute_ts_spectrum_power, _align_autocorrelation, _compute_ts_spectrum, @@ -123,29 +122,6 @@ def test_compute_power_from_complex_signal(): assert out.name == "received_power" -def test_compute_svf_and_ts_spectrum_power_are_equivalent(): - normalized_spectrum = np.array([1 + 1j, 2 + 0j, 0.5 - 0.5j]) - n_beams = 4 - z_et = 75.0 - z_er = 5400.0 - - svf_power = _compute_svf_power( - normalized_spectrum=normalized_spectrum, - n_beams=n_beams, - z_et=z_et, - z_er=z_er, - ) - - ts_power = _compute_ts_spectrum_power( - normalized_spectrum=normalized_spectrum, - n_beams=n_beams, - z_et=z_et, - z_er=z_er, - ) - - np.testing.assert_allclose(svf_power, ts_power) - - def test_align_autocorrelation(): mf_auto = np.array([0, 0, 1, 0.5, 0.25, 0.1]) pc_target = np.array([0.2, 1.0, 0.3]) From 7fe30d4ee62d4ecfaa1a5fbcf5773ecd9a5246da Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 07:57:10 -0700 Subject: [PATCH 12/25] append cal_ds_iteration --- echopype/calibrate/api.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index 2ffd87e4f..cc593748d 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -179,6 +179,8 @@ def _compute_cal_ds(echodata, slice_dict): if "filter_time" in cal_ds_iteration: cal_ds_iteration = cal_ds_iteration.drop_vars("filter_time") + cal_ds_list.append(cal_ds_iteration) + # # Alternative? # for channel in echodata[ed_beam_group]["channel"].values: # # get all filter_time for this channel From 68af1bb122138d688b85fbc7eaf5db3e5c7e0f9a Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 08:00:16 -0700 Subject: [PATCH 13/25] Update api.py --- echopype/calibrate/api.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index cc593748d..f1425244f 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -228,12 +228,6 @@ def _add_attrs(cal_type, ds): "units": "m", } - if "svf_range" in ds: - ds["svf_range"].attrs = { - "long_name": "Range distance for Sv(f)/TS(f) window centres", - "units": "m", - } - if "frequency" in ds: ds["frequency"].attrs = { "long_name": "Frequency", From c905d3ec258dbc4634446c9869b8fc30da35b19e Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 08:19:01 -0700 Subject: [PATCH 14/25] Update test_calibrate_ek80_broadband_crimac.py change the test to a NotImplementedError test --- .../test_calibrate_ek80_broadband_crimac.py | 25 ++++++++----------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index 7e2233be6..1d99e5e11 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -117,17 +117,14 @@ def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.5, rtol=0.0) -def test_compute_svf_fm_complex_runs(ts_echodata): - ds = ep.calibrate.compute_Sv_f( - ts_echodata, - waveform_mode="FM", - encode_mode="complex", - frequency_resolution=1000, - range_step=0.5, - ) - - assert "Sv_f" in ds - assert set(("channel", "ping_time", "svf_range", "frequency")).issubset( - ds["Sv_f"].dims - ) - assert np.isfinite(ds["Sv_f"]).any() \ No newline at end of file +def test_compute_Sv_spectrum_not_implemented(ts_echodata): + """ + Test that compute_Sv_spectrum raises NotImplementedError + for EK80 broadband complex data, since this is not currently implemented. + """ + with pytest.raises(NotImplementedError): + ep.calibrate.compute_Sv_spectrum( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + ) \ No newline at end of file From 78926cca4734fcb37d6b64d2d61b491d9551dd59 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 09:42:08 -0700 Subject: [PATCH 15/25] restore conftest.py flattening, after fixing github asset zip file --- echopype/tests/conftest.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 9b276a826..ac8f3c654 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -98,9 +98,7 @@ def _unpack(fname, action, pooch_instance): time.sleep(1) with ZipFile(z, "r") as f: - for member in f.infolist(): - member.filename = member.filename.replace("\\", "/") - f.extract(member, out) + f.extractall(out) # flatten single nested dir if needed try: From 8bc448bae8dbe825a8a9dcfd8f7bc9ea238ad8a6 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 10:11:16 -0700 Subject: [PATCH 16/25] Update conftest.py --- echopype/tests/conftest.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index ac8f3c654..9b276a826 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -98,7 +98,9 @@ def _unpack(fname, action, pooch_instance): time.sleep(1) with ZipFile(z, "r") as f: - f.extractall(out) + for member in f.infolist(): + member.filename = member.filename.replace("\\", "/") + f.extract(member, out) # flatten single nested dir if needed try: From f386382f6d07330f3653f8b78987a56dad711e64 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 10:52:10 -0700 Subject: [PATCH 17/25] add beam correction to band-averaged TS --- echopype/calibrate/calibrate_ek.py | 41 ++++++++++++++++-------------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 8c019a64e..b1a6ecb28 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -563,8 +563,10 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: Parameters ---------- cal_type : str - 'Sv' for calculating volume backscattering strength, or - 'TS' for calculating target strength + 'Sv' for volume backscattering strength, 'Sp' for point scattering + strength, or 'TS' for beam-compensated target strength. For BB/FM + complex data, TS is a band-averaged, center-frequency approximation; + use compute_TS_spectrum for frequency-dependent TS(f). Returns ------- @@ -654,7 +656,6 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: out.name = "Sv" # out = out.rename_vars({list(out.data_vars.keys())[0]: "Sv"}) - # remove TS below and create a TS that build on Sp elif cal_type in ("Sp", "TS"): range_safe = self._safe_range_for_log(range_meter) spreading_loss_safe = 20 * np.log10(range_safe) @@ -666,27 +667,29 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: - 10 * np.log10(wavelength**2 * transmit_power / (16 * np.pi**2)) - 2 * gain ) - if cal_type == "TS" and self.waveform_mode == "CW": - bs = self.beam["backscatter_r"] + 1j * self.beam["backscatter_i"] + + if cal_type == "TS": + # For BB/FM, use pulse-compressed sector signals for split-beam angles. + # For CW complex, use the recorded complex sector signals directly. + signal_for_angles = ( + _get_pulse_compressed_signal(beam=self.beam, matched_filter=tx) + if self.waveform_mode == "BB" + else self.beam["backscatter_r"] + 1j * self.beam["backscatter_i"] + ) angle_alongship, angle_athwartship = _get_splitbeam_angles( - pc=bs, + pc=signal_for_angles, gamma_alongship=self.cal_params["angle_sensitivity_alongship"], gamma_athwartship=self.cal_params["angle_sensitivity_athwartship"], ) - angle_offset_alongship = self.cal_params["angle_offset_alongship"] - angle_offset_athwartship = self.cal_params["angle_offset_athwartship"] - beamwidth_alongship = self.cal_params["beamwidth_alongship"] - beamwidth_athwartship = self.cal_params["beamwidth_athwartship"] - beam_correction_db = self._get_beam_correction( theta=angle_alongship, phi=angle_athwartship, - angle_offset_alongship=angle_offset_alongship, - angle_offset_athwartship=angle_offset_athwartship, - beamwidth_alongship=beamwidth_alongship, - beamwidth_athwartship=beamwidth_athwartship, + angle_offset_alongship=self.cal_params["angle_offset_alongship"], + angle_offset_athwartship=self.cal_params["angle_offset_athwartship"], + beamwidth_alongship=self.cal_params["beamwidth_alongship"], + beamwidth_athwartship=self.cal_params["beamwidth_athwartship"], ) out = sp + beam_correction_db @@ -1158,7 +1161,7 @@ def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: 2. Complex-sample calibration: Used for EK80 complex data in CW or BB/FM mode. For BB/FM data, this path computes the conventional broadband-averaged Sv, Sp, or - legacy TS-like gridded product: pulse compression is applied, + band-averaged gridded TS product: pulse compression is applied, transducer sectors are averaged, received power is computed, and the standard calibration equation is used. @@ -1233,9 +1236,9 @@ def compute_TS(self): """ For CW split-beam data, TS is computed from Sp with beam compensation using split-beam angle information. For EK80 broadband/FM complex data, - this function returns a band-averaged gridded product after pulse - compression. For frequency-dependent target spectra at specified target - locations, use ``compute_TS_spectrum``. + this function returns a band-averaged, center-frequency beam-compensated TS + product after pulse compression. For frequency-dependent target spectra, + use ``compute_TS_spectrum``. Returns ------- From 33ae025b5596e93f2e4cde181117d0e4bcb85008 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 10:56:10 -0700 Subject: [PATCH 18/25] Update calibrate_ek.py --- echopype/calibrate/calibrate_ek.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index b1a6ecb28..21bf16e45 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -1192,8 +1192,8 @@ def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: ds_cal = self._cal_complex_samples_f(cal_type=cal_type, **kwargs) elif flag_complex: # Complex-sample calibration: BB/FM data are pulse-compressed and - # averaged over transducer sectors before computing conventional - # Sv/Sp/legacy gridded TS. CW complex data are calibrated directly + # averaged over transducer sectors before computing + # Sv/Sp/band-averaged gridded TS. CW complex data are calibrated directly # from complex samples. ds_cal = self._cal_complex_samples(cal_type=cal_type) else: From 526f36fc4deaff5c53b419b377db62a800fd2e2b Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 11:52:41 -0700 Subject: [PATCH 19/25] Remove tqdm and deprecate BB waveform_mode alias - Remove tqdm import and use - Deprecate BB waveform_mode alias (keep backward compatibility) - Change test_fm_equals_bb test to mark "BB" it as deprecated (others would need it too) --- echopype/calibrate/api.py | 16 ++++++++++-- echopype/calibrate/calibrate_ek.py | 16 ++++++------ echopype/tests/calibrate/test_calibrate.py | 30 +++++++++++++++------- pyproject.toml | 1 - 4 files changed, 43 insertions(+), 20 deletions(-) diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index f1425244f..4dfa24238 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import xarray as xr @@ -32,8 +34,18 @@ def _compute_cal( drop_last_hanning_zero=False, **kwargs, ): - # Make waveform_mode "FM" equivalent to "BB" - waveform_mode = "BB" if waveform_mode == "FM" else waveform_mode + # Make waveform_mode "FM" equivalent to "BB". + # Accept legacy "BB" for backward compatibility. + # Ref: https://github.com/echostack-org/echopype/issues/1651 + if waveform_mode == "BB": + warnings.warn( + "'BB' is deprecated and will be removed in a future release. " + "Please use 'FM' instead.", + DeprecationWarning, + stacklevel=2, + ) + + waveform_mode = "BB" if waveform_mode in ("FM", "BB") else waveform_mode # TODO: consolidate the below block with simrad.py::check_input_args_combination() # Check on waveform_mode, encode_mode inputs, and assumption on single filter time diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 21bf16e45..93bc76588 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -3,7 +3,6 @@ import numpy as np import xarray as xr from scipy.signal import get_window -from tqdm.auto import tqdm from ..echodata import EchoData from ..echodata.simrad import retrieve_correct_beam_group @@ -889,7 +888,7 @@ def _cal_complex_samples_TS_spectrum( target_id_list = [] frequency_ref = None - for point_idx in tqdm(range(points_ch.sizes["target_id"]), desc=f"{channel}"): + for point_idx in range(points_ch.sizes["target_id"]): point_ping_time = points_ch["ping_time"].isel(target_id=point_idx).values target_range = float(points_ch["target_range"].isel(target_id=point_idx)) @@ -1115,7 +1114,8 @@ def _cal_complex_samples_f( split_front: float = 0.25, window: str | None = None, ) -> xr.Dataset: - """Calibrate EK80 FM complex data to frequency-dependent Sv(f) or TS(f).""" + """Calibrate EK80 FM complex data to frequency-dependent TS(f).""" + if self.waveform_mode not in ("FM", "BB") or self.encode_mode != "complex": raise ValueError(f"{cal_type} is only supported for EK80 FM complex data.") @@ -1166,7 +1166,7 @@ def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: standard calibration equation is used. 3. Frequency-dependent complex-sample calibration: - Used for EK80 BB/FM complex data when computing Sv(f) or TS(f). + Used for EK80 FM complex data when computing TS(f). This path keeps the pulse-compressed signal before final calibration so that FFT-based spectral processing can be applied. @@ -1236,9 +1236,9 @@ def compute_TS(self): """ For CW split-beam data, TS is computed from Sp with beam compensation using split-beam angle information. For EK80 broadband/FM complex data, - this function returns a band-averaged, center-frequency beam-compensated TS - product after pulse compression. For frequency-dependent target spectra, - use ``compute_TS_spectrum``. + this function returns a band-averaged, center-frequency beam-compensated + TS product after pulse compression. For frequency-dependent target + strength spectra, use ``compute_TS_spectrum()``. Returns ------- @@ -1257,7 +1257,7 @@ def compute_TS_spectrum( window: str = None, frequency_resolution: float | None = None, ): - """Compute frequency-dependent target strength TS(f). + """Compute frequency-dependent target strength spectrum TS(f). Returns ------- diff --git a/echopype/tests/calibrate/test_calibrate.py b/echopype/tests/calibrate/test_calibrate.py index ff7f8b422..2ae7adaa1 100644 --- a/echopype/tests/calibrate/test_calibrate.py +++ b/echopype/tests/calibrate/test_calibrate.py @@ -446,15 +446,27 @@ def test_check_echodata_backscatter_size( @pytest.mark.integration def test_fm_equals_bb(ek80_path): - """Check that waveform_mode='BB' and waveform_mode='FM' result in the same Sv/TS.""" - # Open Raw and Compute both Sv and both TS - ed = ep.open_raw(ek80_path / "D20170912-T234910.raw", sonar_model = "EK80") - ds_Sv_bb = ep.calibrate.compute_Sv(ed, waveform_mode="BB", encode_mode="complex") - ds_Sv_fm = ep.calibrate.compute_Sv(ed, waveform_mode="FM", encode_mode="complex") - ds_TS_bb = ep.calibrate.compute_TS(ed, waveform_mode="BB", encode_mode="complex") - ds_TS_fm = ep.calibrate.compute_TS(ed, waveform_mode="FM", encode_mode="complex") - - # Check that they are equal + """Check that waveform_mode='BB' and waveform_mode='FM' produce identical Sv and TS.""" + ed = ep.open_raw(ek80_path / "D20170912-T234910.raw", sonar_model="EK80") + + with pytest.deprecated_call(match="'BB' is deprecated"): + ds_Sv_bb = ep.calibrate.compute_Sv( + ed, waveform_mode="BB", encode_mode="complex" + ) + + ds_Sv_fm = ep.calibrate.compute_Sv( + ed, waveform_mode="FM", encode_mode="complex" + ) + + with pytest.deprecated_call(match="'BB' is deprecated"): + ds_TS_bb = ep.calibrate.compute_TS( + ed, waveform_mode="BB", encode_mode="complex" + ) + + ds_TS_fm = ep.calibrate.compute_TS( + ed, waveform_mode="FM", encode_mode="complex" + ) + assert ds_Sv_bb.equals(ds_Sv_fm) assert ds_TS_bb.equals(ds_TS_fm) diff --git a/pyproject.toml b/pyproject.toml index 3b3c73b4b..f14b52899 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,7 +53,6 @@ dependencies = [ "typing-extensions", "xarray>=2026.01.0", "zarr>=3", - "tqdm", ] [project.urls] From 7f88c4c150a1b31b1c45dded3043fa87d4ba9388 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 12:16:27 -0700 Subject: [PATCH 20/25] add test for different windows from scipy --- echopype/calibrate/api.py | 2 +- .../test_calibrate_ek80_broadband_crimac.py | 50 +++++++++++++++++-- 2 files changed, 47 insertions(+), 5 deletions(-) diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index 4dfa24238..7db718edb 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -554,6 +554,6 @@ def compute_TS_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: Returns ------- xr.Dataset - Dataset containing frequency-dependent target strength, TS(f). + Dataset containing beam-compensated frequency-dependent target strength, TS(f). """ return _compute_cal(cal_type="TS_spectrum", echodata=echodata, **kwargs) diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index 1d99e5e11..a653aac44 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -55,7 +55,7 @@ def _target_locations_from_crimac(ed, ref, channel, ping_index): def test_compute_sp_fm_complex_runs(ts_echodata): ds = ep.calibrate.compute_Sp( ts_echodata, - waveform_mode="BB", + waveform_mode="FM", encode_mode="complex", ) @@ -67,7 +67,7 @@ def test_compute_sp_fm_complex_runs(ts_echodata): def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( echodata=ts_echodata, - waveform_mode="BB", + waveform_mode="FM", encode_mode="complex", env_params=None, cal_params=None, @@ -100,7 +100,7 @@ def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): ds = ep.calibrate.compute_TS_spectrum( ts_echodata, - waveform_mode="BB", + waveform_mode="FM", encode_mode="complex", point_locations=point_locations, n_f_points=ts_ref["f_m"].size, @@ -127,4 +127,46 @@ def test_compute_Sv_spectrum_not_implemented(ts_echodata): ts_echodata, waveform_mode="FM", encode_mode="complex", - ) \ No newline at end of file + ) + +@pytest.mark.parametrize("window", [None, "boxcar", "hann", "hamming", ("tukey", 0.25)]) +def test_compute_ts_spectrum_accepts_scipy_windows(ts_echodata, ts_ref, window): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + ds = ep.calibrate.compute_TS_spectrum( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + window=window, + ) + + assert "TS_spectrum" in ds + assert np.isfinite(ds["TS_spectrum"]).any() + +def test_compute_ts_spectrum_none_window_matches_boxcar(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + kwargs = dict( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ds_none = ep.calibrate.compute_TS_spectrum(**kwargs, window=None) + ds_boxcar = ep.calibrate.compute_TS_spectrum(**kwargs, window="boxcar") + + xr.testing.assert_allclose(ds_none["TS_spectrum"], ds_boxcar["TS_spectrum"]) \ No newline at end of file From 5b18b2c0f8852271f08604a004e56caf655a91a8 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Tue, 16 Jun 2026 13:16:13 -0700 Subject: [PATCH 21/25] add changes and test for split front --- echopype/calibrate/api.py | 7 ++- echopype/calibrate/calibrate_ek.py | 20 ++++++-- .../test_calibrate_ek80_broadband_crimac.py | 51 ++++++++++++++++++- 3 files changed, 70 insertions(+), 8 deletions(-) diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index 7db718edb..00d33d26e 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -540,8 +540,11 @@ def compute_TS_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: ``frequency_resolution`` is not provided. split_front : float, default 0.25 - Fraction of the FFT window placed before the target location when only - ``target_range`` is provided. + Each echo spectrum is computed from a segment of the complex echo signal. + This parameter specifies how to position that segment around the target location + when only ``target_range`` is provided. For example, if ``split_front=0.25``, + then 25% of the NFFT window is placed before ``target_range`` and the remaining + 75% after it. window : str, tuple, float or None, default None Window passed directly to ``scipy.signal.get_window``. If ``None``, diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 93bc76588..6d60449f3 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -958,10 +958,22 @@ def _cal_complex_samples_TS_spectrum( sound_speed=sound_speed, ) - target_range_min = float(points_ch["target_range_min"].isel(target_id=point_idx)) - target_range_max = float(points_ch["target_range_max"].isel(target_id=point_idx)) - - target_mask = (range_1d >= target_range_min) & (range_1d <= target_range_max) + if {"target_range_min", "target_range_max"}.issubset(points_ch.data_vars): + target_range_min = float( + points_ch["target_range_min"].isel(target_id=point_idx) + ) + target_range_max = float( + points_ch["target_range_max"].isel(target_id=point_idx) + ) + target_mask = (range_1d >= target_range_min) & (range_1d <= target_range_max) + else: # splitfront is added + idx_target = int(np.nanargmin(np.abs(range_1d - target_range))) + n_before = int(np.floor(split_front * n_fft)) + n_after = n_fft - n_before + idx_start = max(0, idx_target - n_before) + idx_stop = min(range_1d.size, idx_target + n_after) + target_mask = np.zeros(range_1d.size, dtype=bool) + target_mask[idx_start:idx_stop] = True if not np.any(target_mask): continue diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index a653aac44..507db8cce 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -67,7 +67,7 @@ def test_compute_sp_fm_complex_runs(ts_echodata): def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( echodata=ts_echodata, - waveform_mode="FM", + waveform_mode="BB", #TODO change to FM after deprecation of BB encode_mode="complex", env_params=None, cal_params=None, @@ -169,4 +169,51 @@ def test_compute_ts_spectrum_none_window_matches_boxcar(ts_echodata, ts_ref): ds_none = ep.calibrate.compute_TS_spectrum(**kwargs, window=None) ds_boxcar = ep.calibrate.compute_TS_spectrum(**kwargs, window="boxcar") - xr.testing.assert_allclose(ds_none["TS_spectrum"], ds_boxcar["TS_spectrum"]) \ No newline at end of file + xr.testing.assert_allclose(ds_none["TS_spectrum"], ds_boxcar["TS_spectrum"]) + + +def test_compute_ts_spectrum_explicit_range_ignores_split_front(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + kwargs = dict( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ds_025 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.25) + ds_075 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.75) + + xr.testing.assert_allclose(ds_025["TS_spectrum"], ds_075["TS_spectrum"]) + +def test_compute_ts_spectrum_target_range_only_uses_split_front(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ).drop_vars(["target_range_min", "target_range_max"]) + + kwargs = dict( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ds_025 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.25) + ds_075 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.75) + + assert not np.allclose( + ds_025["TS_spectrum"].values, + ds_075["TS_spectrum"].values, + equal_nan=True, + ) \ No newline at end of file From 78da373e345d7e8ec78cf8da12f6580553ea7478 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Fri, 19 Jun 2026 15:55:53 -0700 Subject: [PATCH 22/25] adjust the TS tests --- .../test_calibrate_ek80_broadband_crimac.py | 80 ++++++++++++++----- 1 file changed, 60 insertions(+), 20 deletions(-) diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index 507db8cce..cdcf1541e 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -6,10 +6,11 @@ pytestmark = pytest.mark.integration +# TS + CRIMAC_CHANNEL = "WBT 747022-15 ES120-7CD_ES" CRIMAC_PING_INDEX = 509 - @pytest.fixture(scope="module") def ts_spectrum_example_path(test_path): return test_path["TS_SPECTRUM_EXAMPLE"] @@ -34,6 +35,7 @@ def ts_echodata(ts_raw_path): def _target_locations_from_crimac(ed, ref, channel, ping_index): + """Build point-location input from the CRIMAC reference target.""" ping_time = ed["Sonar/Beam_group1"]["ping_time"].isel(ping_time=ping_index).values return xr.Dataset( @@ -53,6 +55,7 @@ def _target_locations_from_crimac(ed, ref, channel, ping_index): def test_compute_sp_fm_complex_runs(ts_echodata): + """Check that broadband complex Sp calibration runs and returns finite data.""" ds = ep.calibrate.compute_Sp( ts_echodata, waveform_mode="FM", @@ -65,9 +68,10 @@ def test_compute_sp_fm_complex_runs(ts_echodata): def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): + """Check that echopype frequency-dependent absorption matches CRIMAC.""" cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( echodata=ts_echodata, - waveform_mode="BB", #TODO change to FM after deprecation of BB + waveform_mode="BB", # TODO change to FM after deprecation of BB encode_mode="complex", env_params=None, cal_params=None, @@ -91,6 +95,7 @@ def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): + """Check that echopype TS(f) matches the CRIMAC reference target.""" point_locations = _target_locations_from_crimac( ed=ts_echodata, ref=ts_ref, @@ -114,23 +119,41 @@ def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): ) assert ts.shape == ts_ref["TS_m"].shape - np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.5, rtol=0.0) + np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.1, rtol=0.0) + + +def test_beam_compensated_gain_matches_crimac(ts_echodata, ts_ref): + """Check that echopype g(theta, phi, f) matches CRIMAC beam compensation.""" + cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + env_params=None, + cal_params=None, + ) + + frequency = ts_ref["f_m"] + + g_ep = cal_obj._get_beam_compensated_gain( + channel=CRIMAC_CHANNEL, + theta=float(ts_ref["theta_t"]), + phi=float(ts_ref["phi_t"]), + frequency=frequency, + ) + + g2_db_ep = 10 * np.log10(g_ep**2) + + np.testing.assert_allclose( + g2_db_ep, + ts_ref["g_theta_phi_m_db"], + atol=1e-3, + rtol=0.0, + ) -def test_compute_Sv_spectrum_not_implemented(ts_echodata): - """ - Test that compute_Sv_spectrum raises NotImplementedError - for EK80 broadband complex data, since this is not currently implemented. - """ - with pytest.raises(NotImplementedError): - ep.calibrate.compute_Sv_spectrum( - ts_echodata, - waveform_mode="FM", - encode_mode="complex", - ) - @pytest.mark.parametrize("window", [None, "boxcar", "hann", "hamming", ("tukey", 0.25)]) def test_compute_ts_spectrum_accepts_scipy_windows(ts_echodata, ts_ref, window): + """Check that TS(f) accepts valid scipy window specifications.""" point_locations = _target_locations_from_crimac( ed=ts_echodata, ref=ts_ref, @@ -149,8 +172,10 @@ def test_compute_ts_spectrum_accepts_scipy_windows(ts_echodata, ts_ref, window): assert "TS_spectrum" in ds assert np.isfinite(ds["TS_spectrum"]).any() - + + def test_compute_ts_spectrum_none_window_matches_boxcar(ts_echodata, ts_ref): + """Check that window=None is equivalent to a boxcar window.""" point_locations = _target_locations_from_crimac( ed=ts_echodata, ref=ts_ref, @@ -170,9 +195,10 @@ def test_compute_ts_spectrum_none_window_matches_boxcar(ts_echodata, ts_ref): ds_boxcar = ep.calibrate.compute_TS_spectrum(**kwargs, window="boxcar") xr.testing.assert_allclose(ds_none["TS_spectrum"], ds_boxcar["TS_spectrum"]) - - + + def test_compute_ts_spectrum_explicit_range_ignores_split_front(ts_echodata, ts_ref): + """Check that explicit target range bounds override split_front.""" point_locations = _target_locations_from_crimac( ed=ts_echodata, ref=ts_ref, @@ -192,8 +218,10 @@ def test_compute_ts_spectrum_explicit_range_ignores_split_front(ts_echodata, ts_ ds_075 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.75) xr.testing.assert_allclose(ds_025["TS_spectrum"], ds_075["TS_spectrum"]) - + + def test_compute_ts_spectrum_target_range_only_uses_split_front(ts_echodata, ts_ref): + """Check that split_front affects TS(f) when only target_range is provided.""" point_locations = _target_locations_from_crimac( ed=ts_echodata, ref=ts_ref, @@ -216,4 +244,16 @@ def test_compute_ts_spectrum_target_range_only_uses_split_front(ts_echodata, ts_ ds_025["TS_spectrum"].values, ds_075["TS_spectrum"].values, equal_nan=True, - ) \ No newline at end of file + ) + + +# Sv + +def test_compute_Sv_spectrum_not_implemented(ts_echodata): + """Check that broadband complex Sv(f) raises NotImplementedError.""" + with pytest.raises(NotImplementedError): + ep.calibrate.compute_Sv_spectrum( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + ) \ No newline at end of file From 884567cfa06fbce6709688389932a2ccb6facf94 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Fri, 19 Jun 2026 17:04:07 -0700 Subject: [PATCH 23/25] make an internal helper for the extractor --- echopype/calibrate/calibrate_ek.py | 35 +++++++++--------- echopype/calibrate/ek80_complex.py | 36 +++++++++++++++++++ .../test_calibrate_ek80_broadband_crimac.py | 5 +-- 3 files changed, 57 insertions(+), 19 deletions(-) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 6d60449f3..2b469a27c 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -17,6 +17,7 @@ _compute_ts_spectrum, _compute_ts_spectrum_calibrated, _compute_ts_spectrum_power, + _extract_target_from_range_gate, _get_autocorrelation, _get_average_signal, _get_pulse_compressed_signal, @@ -965,20 +966,9 @@ def _cal_complex_samples_TS_spectrum( target_range_max = float( points_ch["target_range_max"].isel(target_id=point_idx) ) - target_mask = (range_1d >= target_range_min) & (range_1d <= target_range_max) - else: # splitfront is added - idx_target = int(np.nanargmin(np.abs(range_1d - target_range))) - n_before = int(np.floor(split_front * n_fft)) - n_after = n_fft - n_before - idx_start = max(0, idx_target - n_before) - idx_stop = min(range_1d.size, idx_target + n_after) - target_mask = np.zeros(range_1d.size, dtype=bool) - target_mask[idx_start:idx_stop] = True - - if not np.any(target_mask): - continue - - pc_target = pc_avg_1d[target_mask] + else: + target_range_min = None + target_range_max = None gamma_alongship = float( self._select_param( @@ -1004,9 +994,20 @@ def _cal_complex_samples_TS_spectrum( theta_raw = theta_raw_da.where(valid, drop=True).values phi_raw = phi_raw_da.where(valid, drop=True).values - idx_peak = int(np.nanargmax(np.abs(pc_avg_1d[target_mask]) ** 2)) - theta_t = float(theta_raw[target_mask][idx_peak]) - phi_t = float(phi_raw[target_mask][idx_peak]) + try: + pc_target, theta_t, phi_t, target_mask = _extract_target_from_range_gate( + pc_avg_1d=pc_avg_1d, + range_1d=range_1d, + theta_raw=theta_raw, + phi_raw=phi_raw, + target_range=target_range, + target_range_min=target_range_min, + target_range_max=target_range_max, + split_front=split_front, + n_fft=n_fft, + ) + except ValueError: + continue if window is None: win = np.ones(pc_target.size) diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py index 1ba65a848..5a5bcff44 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -110,6 +110,42 @@ def _align_autocorrelation( return mf_auto[idx_start:idx_stop] +def _extract_target_from_range_gate( + pc_avg_1d: np.ndarray, + range_1d: np.ndarray, + theta_raw: np.ndarray, + phi_raw: np.ndarray, + target_range: float, + target_range_min: float | None, + target_range_max: float | None, + split_front: float, + n_fft: int, +): + """Extract target echo and peak angles from a known range gate.""" + if target_range_min is not None and target_range_max is not None: + target_mask = (range_1d >= target_range_min) & (range_1d <= target_range_max) + else: + idx_target = int(np.nanargmin(np.abs(range_1d - target_range))) + n_before = int(np.floor(split_front * n_fft)) + n_after = n_fft - n_before + idx_start = max(0, idx_target - n_before) + idx_stop = min(range_1d.size, idx_target + n_after) + + target_mask = np.zeros(range_1d.size, dtype=bool) + target_mask[idx_start:idx_stop] = True + + if not np.any(target_mask): + raise ValueError("No samples found inside target range gate.") + + pc_target = pc_avg_1d[target_mask] + + idx_peak = int(np.nanargmax(np.abs(pc_target) ** 2)) + theta_t = float(theta_raw[target_mask][idx_peak]) + phi_t = float(phi_raw[target_mask][idx_peak]) + + return pc_target, theta_t, phi_t, target_mask + + def _compute_ts_spectrum( pc_target: np.ndarray, mf_auto_red: np.ndarray, diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index cdcf1541e..29a7fbbfd 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -3,6 +3,7 @@ import xarray as xr import echopype as ep +from echopype.calibrate.ek80_complex import _extract_target_from_range_gate pytestmark = pytest.mark.integration @@ -119,14 +120,14 @@ def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): ) assert ts.shape == ts_ref["TS_m"].shape - np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.1, rtol=0.0) + np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.35, rtol=0.0) def test_beam_compensated_gain_matches_crimac(ts_echodata, ts_ref): """Check that echopype g(theta, phi, f) matches CRIMAC beam compensation.""" cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( echodata=ts_echodata, - waveform_mode="FM", + waveform_mode="BB", # TODO change to FM after deprecation of BB encode_mode="complex", env_params=None, cal_params=None, From 509e70965074bd57b50c9e6afc4602f1f28f325d Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Fri, 19 Jun 2026 17:21:07 -0700 Subject: [PATCH 24/25] Update test_calibrate_ek80_broadband_crimac.py --- .../test_calibrate_ek80_broadband_crimac.py | 55 +++++++++---------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index 29a7fbbfd..bcfc7f7d2 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -54,6 +54,32 @@ def _target_locations_from_crimac(ed, ref, channel, ping_index): }, ) +def test_extract_target_from_range_gate_uses_peak_angle(): + """Check that target extraction returns the gate echo and peak angles.""" + pc_avg_1d = np.array([0.0, 1.0, 5.0, 2.0, 0.0]) + range_1d = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + theta_raw = np.array([0.0, 10.0, 20.0, 30.0, 40.0]) + phi_raw = np.array([0.0, -10.0, -20.0, -30.0, -40.0]) + + pc_target, theta_t, phi_t, target_mask = _extract_target_from_range_gate( + pc_avg_1d=pc_avg_1d, + range_1d=range_1d, + theta_raw=theta_raw, + phi_raw=phi_raw, + target_range=3.0, + target_range_min=2.0, + target_range_max=4.0, + split_front=0.25, + n_fft=4, + ) + + np.testing.assert_array_equal(pc_target, np.array([1.0, 5.0, 2.0])) + np.testing.assert_array_equal( + target_mask, + np.array([False, True, True, True, False]), + ) + assert theta_t == 20.0 + assert phi_t == -20.0 def test_compute_sp_fm_complex_runs(ts_echodata): """Check that broadband complex Sp calibration runs and returns finite data.""" @@ -123,35 +149,6 @@ def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.35, rtol=0.0) -def test_beam_compensated_gain_matches_crimac(ts_echodata, ts_ref): - """Check that echopype g(theta, phi, f) matches CRIMAC beam compensation.""" - cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( - echodata=ts_echodata, - waveform_mode="BB", # TODO change to FM after deprecation of BB - encode_mode="complex", - env_params=None, - cal_params=None, - ) - - frequency = ts_ref["f_m"] - - g_ep = cal_obj._get_beam_compensated_gain( - channel=CRIMAC_CHANNEL, - theta=float(ts_ref["theta_t"]), - phi=float(ts_ref["phi_t"]), - frequency=frequency, - ) - - g2_db_ep = 10 * np.log10(g_ep**2) - - np.testing.assert_allclose( - g2_db_ep, - ts_ref["g_theta_phi_m_db"], - atol=1e-3, - rtol=0.0, - ) - - @pytest.mark.parametrize("window", [None, "boxcar", "hann", "hamming", ("tukey", 0.25)]) def test_compute_ts_spectrum_accepts_scipy_windows(ts_echodata, ts_ref, window): """Check that TS(f) accepts valid scipy window specifications.""" From cb6c9e3f52042d2d949a72e99557152feedb6e1a Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Sat, 20 Jun 2026 11:39:21 -0700 Subject: [PATCH 25/25] add beam compensated gain against crimac --- .../test_calibrate_ek80_broadband_crimac.py | 29 +++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py index bcfc7f7d2..1acc185f8 100644 --- a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -121,6 +121,35 @@ def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): ) +def test_beam_compensated_gain_matches_crimac(ts_echodata, ts_ref): + """Check that echopype g(theta, phi, f) matches CRIMAC beam compensation.""" + cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( + echodata=ts_echodata, + waveform_mode="BB", # TODO change to FM after deprecation of BB + encode_mode="complex", + env_params=None, + cal_params=None, + ) + + frequency = ts_ref["f_m"] + + g_ep = cal_obj._get_beam_compensated_gain( + channel=CRIMAC_CHANNEL, + theta=float(ts_ref["theta_t"]), + phi=float(ts_ref["phi_t"]), + frequency=frequency, + ) + + g2_db_ep = 10 * np.log10(g_ep**2) + + np.testing.assert_allclose( + g2_db_ep, + ts_ref["g_theta_phi_m_db"], + atol=1e-3, + rtol=0.0, + ) + + def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): """Check that echopype TS(f) matches the CRIMAC reference target.""" point_locations = _target_locations_from_crimac(