Add functions for broadband TS spectrum

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"]

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)