diff --git a/echopype/calibrate/__init__.py b/echopype/calibrate/__init__.py index 94613c73b..8fb4e732e 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_Sp, compute_Sv, compute_Sv_spectrum, compute_TS, compute_TS_spectrum -__all__ = ["compute_Sv", "compute_TS"] +__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 d12b27b7c..00d33d26e 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 @@ -21,7 +23,7 @@ def _compute_cal( - cal_type, + cal_type: str, echodata: EchoData, env_params=None, cal_params=None, @@ -30,9 +32,20 @@ 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 + # 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 @@ -80,13 +93,27 @@ 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 = { + "Sp": "compute_Sp", + "TS": "compute_TS", + "Sv": "compute_Sv", + # add Sp_spectrum?? + "TS_spectrum": "compute_TS_spectrum", + } - return cal_ds + try: + 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 if echodata.sonar_model not in ["EK80", "ES80", "EA640"]: @@ -161,7 +188,10 @@ 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") + + cal_ds_list.append(cal_ds_iteration) # # Alternative? # for channel in echodata[ed_beam_group]["channel"].values: @@ -201,12 +231,27 @@ 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 "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)", + "Sp": "Point scattering strength (Sp re 1 m^2)", "TS": "Target strength (TS re 1 m^2)", + "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", } @@ -345,6 +390,29 @@ def compute_Sv(echodata: EchoData, **kwargs) -> xr.Dataset: return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs) +def compute_Sv_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute frequency-dependent volume backscattering strength Sv(f) + from broadband EK80 complex data. + + Notes + ----- + This functionality is not yet implemented. + """ + raise NotImplementedError("compute_Sv_spectrum is not yet implemented.") + + +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. @@ -447,3 +515,48 @@ 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_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + 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 + 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``, + 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 beam-compensated 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 b58be2ca7..2b469a27c 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -2,17 +2,27 @@ import numpy as np import xarray as xr +from scipy.signal import get_window 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_ts_spectrum, + _compute_ts_spectrum_calibrated, + _compute_ts_spectrum_power, + _extract_target_from_range_gate, + _get_autocorrelation, + _get_average_signal, + _get_pulse_compressed_signal, + _get_splitbeam_angles, get_filter_coeff, - get_norm_fac, get_tau_effective, get_transmit_signal, ) @@ -171,17 +181,38 @@ 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"] + + 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() @@ -264,6 +295,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 +514,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" @@ -535,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 ------- @@ -626,16 +656,47 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: out.name = "Sv" # out = out.rename_vars({list(out.data_vars.keys())[0]: "Sv"}) + 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 = ( + sp = ( 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" + + 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=signal_for_angles, + gamma_alongship=self.cal_params["angle_sensitivity_alongship"], + gamma_athwartship=self.cal_params["angle_sensitivity_athwartship"], + ) + + beam_correction_db = self._get_beam_correction( + theta=angle_alongship, + phi=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"], + ) + + 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) @@ -658,53 +719,570 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: return out - def _compute_cal(self, cal_type) -> xr.Dataset: + def _select_param( + self, + 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", + ) + + def _get_beam_correction( + self, + theta, + phi, + angle_offset_alongship, + angle_offset_athwartship, + beamwidth_alongship, + beamwidth_athwartship, + ): + """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) + + 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)``. + + # 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. + """ - Private method to compute Sv or TS from EK80 data, called by compute_Sv or compute_TS. + 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, + ) + + 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) + + 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 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)) + + 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, + ) + + 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) + ) + else: + target_range_min = None + target_range_max = None + + 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 + + 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) + else: + 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, + ) + + 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") + + def _cal_complex_samples_f( + self, + cal_type: str, + frequency_resolution: float | None = None, + 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 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 == "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 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 + band-averaged gridded TS 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 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. Parameters ---------- cal_type : str - 'Sv' for calculating volume backscattering strength, or - 'TS' for calculating target strength + Calibration type. Supported values include ``"Sv"``, ``"Sp"``, + ``"TS"``, 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 flag_complex: - # Complex samples can be BB or CW + 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) + elif flag_complex: + # Complex-sample calibration: BB/FM data are pulse-compressed and + # 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: - # 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) from raw data. + + 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 + ------- + 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). + """ + 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 + strength spectra, use ``compute_TS_spectrum()``. Returns ------- - TS : xr.DataSet + TS : xr.Dataset A DataSet containing target strength (``TS``) and the corresponding range (``echo_range``) in units meter. """ return self._compute_cal(cal_type="TS") + + 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 spectrum TS(f). + + Returns + ------- + TS_f : xr.Dataset + A Dataset containing frequency-dependent target strength. + """ + 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..5a5bcff44 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -9,6 +9,272 @@ 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, + 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 _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, + 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 _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 _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..2ae7adaa1 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 @@ -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/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..1acc185f8 --- /dev/null +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -0,0 +1,286 @@ +import numpy as np +import pytest +import xarray as xr + +import echopype as ep +from echopype.calibrate.ek80_complex import _extract_target_from_range_gate + +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"] + + +@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): + """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( + 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_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.""" + ds = ep.calibrate.compute_Sp( + ts_echodata, + waveform_mode="FM", + 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): + """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 + 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_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( + 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, + ) + + 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.35, 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.""" + 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): + """Check that window=None is equivalent to a boxcar window.""" + 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"]) + + +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, + 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): + """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, + 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, + ) + + +# 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 diff --git a/echopype/tests/calibrate/test_ek80_complex.py b/echopype/tests/calibrate/test_ek80_complex.py index 595800496..e5bf54853 100644 --- a/echopype/tests/calibrate/test_ek80_complex.py +++ b/echopype/tests/calibrate/test_ek80_complex.py @@ -2,7 +2,15 @@ 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_ts_spectrum_power, + _align_autocorrelation, + _compute_ts_spectrum, + _compute_ts_spectrum_calibrated, +) pytestmark = pytest.mark.unit @@ -73,3 +81,103 @@ 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_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 734d9836c..f1219b5e8 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: