Repository: OSOceanAcoustics/echopype (language: python)

Issue:
Overhaul calibration code
There are a number of issues with the current calibration implementation -- listed below and will be updated as we progress. Now that we've gotten the data format settled down, it's time to carefully review and overhaul this part of the codebase.

### Issues tracked:
- [x] #743 
- [ ] #702 
- [x] #737 
- [x] #458
  - I think this is already addressed, but let's review and make sure of it
- [x] #308 
- [x] #887
- [ ] #925
Allow feeding in frequency-dependent gain factor for broadband calibration
While working on #446 it became obvious that there is a need for users to be able to put in frequency-dependent gain factors externally and not relying on what is or isn't in the data. This is especially true since a lot of files don't already have the frequency-dependent gain factors stored (i.e. the `Vendor` group does not have a data variable called `gain` with dimension `cal_channel_id` and `cal_frequency`).

This seems can be handled in 2 ways:
1. Similar to how we update the platform group to add missing GPS data using `echodata.update_platform()`
2. Feeding in a `gain` xarray DataArray as part of the input argument `cal_params` when calling `compute_Sv`

Pinging @emiliom here since the route we choose is related to what we consider should be an integral part of the data file. 
Perhaps both should be supported? Users may want to apply different gain factors from different calibration sessions for comparison, so route 2 is useful. 
Implement `add_splitbeam_angle` function
This PR addresses issue #890 by creating the structure of the `add_splitbeam_angle` function and implementing several of the core functions of the routine.

Note:

This PR is currently a draft and the following items need to be addressed/discussed: 

- [x] We currently have the function `check_waveform_encode_mode` that ensures the user provided input is appropriate. I have delayed creating a unit test for this function until we agree on its implementation. This unit test should be completed once the function is finalized. 
- [x] Calculating the split-beam angle for the options `waveform_mode="BB"`, `encode_mode="complex"`, and `pulse_compression=False` has been implemented. However, the output is currently not matching the output of the echoview file we are comparing it against in `test_add_splitbeam_angle` for test ID `"ek80_BB_complex_no_pulse"`. This should be investigated.
- [x] Currently, there is no implementation for `waveform_mode="BB"`, `encode_mode="complex"`, and `pulse_compression=True`. We need to decide if this option will be added during this PR. 
- [x] How/should we account for the case where the Sv data sit on disk as a zarr store?
    - Originally when considering #890 we had thought that we could add the angle data to the Sv zarr store on disk. However, I am finding it difficult to locate the zarr store path that our `ds_Sv` is linked to. I think this may be do to xarray not including (or hiding) this information when `open_mfdataset` is used. If we were using just `open_dataset`, we could just look at the encoding (e.g. `ds.encoding['source']`) and it would tell us where the zarr store was. This should be investigated further, if we want to directly write to the zarr store that Sv is linked to. 


Discussion / hints:

Are there specific, convention defined variables for gain factors? Can you point me to them (name(s) and group)
It's a little involved...

There is a `gain_correction` data variable under the `Beam` group, which I decided to use to store the gain factor for CW mode data. There is another `gain_correction` data variable in the `Vendor` group that I decided to store, because it contains info that is more comprehensive than that stored under the `Beam` group (basically this is a look up table with coordinate dimension `frequency` and `pulse_length` -- the `gain_correction` in the `Beam` group is the result of indexing using the actual `pulse_length`)

For broadband data there could be additional calibration results (frequency-dependent gain factors) saved in the config XML, which are currently stored under data variable `gain` in the `Vendor` group with coordinate dimensions `cal_channel_id` and `cal_frequency`. The data variable name `gain` here should probably be changed to something like `gain_correction_broadband` or something similar to be consistent but distinguishable from the `gain_correction` defined in the convention.

Note for BB data files the narrowband gain (the ones used for CW data) still exist, probably because the echosounder can function in both modes so it is more convenient for them to just be included regardless.



Fix patch:
diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py
index f74ac0400..588bf27a4 100644
--- a/echopype/calibrate/api.py
+++ b/echopype/calibrate/api.py
@@ -26,6 +26,9 @@ def _compute_cal(
     waveform_mode=None,
     encode_mode=None,
 ):
+    # TODO: change below to use echodata.simrad::_check_input_args_combination
+    #       for checking Simrad sonar models
+
     # Check on waveform_mode and encode_mode inputs
     if echodata.sonar_model == "EK80":
         if waveform_mode is None or encode_mode is None:
@@ -59,6 +62,15 @@ def _compute_cal(
         encode_mode=encode_mode,
     )
 
+    # Perform calibration
+    if cal_type == "Sv":
+        cal_ds = cal_obj.compute_Sv()
+    elif cal_type == "TS":
+        cal_ds = cal_obj.compute_TS()
+    else:
+        raise ValueError("cal_type must be Sv or TS")
+
+    # Add attributes
     def add_attrs(cal_type, ds):
         """Add attributes to backscattering strength dataset.
         cal_type: Sv or TS
@@ -84,13 +96,9 @@ def add_attrs(cal_type, ds):
                 }
             )
 
-    # Perform calibration
-    if cal_type == "Sv":
-        cal_ds = cal_obj.compute_Sv(waveform_mode=waveform_mode, encode_mode=encode_mode)
-    else:
-        cal_ds = cal_obj.compute_TS(waveform_mode=waveform_mode, encode_mode=encode_mode)
-
     add_attrs(cal_type, cal_ds)
+
+    # Add provinance
     # Provenance source files may originate from raw files (echodata.source_files)
     # or converted files (echodata.converted_raw_path)
     if echodata.source_file is not None:
@@ -109,8 +117,8 @@ def add_attrs(cal_type, ds):
         .assign_attrs(prov_dict)
     )
 
+    # Add water_level to the created xr.Dataset
     if "water_level" in echodata["Platform"].data_vars.keys():
-        # add water_level to the created xr.Dataset
         cal_ds["water_level"] = echodata["Platform"].water_level
 
     return cal_ds
diff --git a/echopype/calibrate/cal_params.py b/echopype/calibrate/cal_params.py
new file mode 100644
index 000000000..a9974a5b7
--- /dev/null
+++ b/echopype/calibrate/cal_params.py
@@ -0,0 +1,224 @@
+from typing import Dict, Union
+
+import numpy as np
+import xarray as xr
+
+from ..echodata import EchoData
+
+CAL_PARAMS = {
+    "EK": ("sa_correction", "gain_correction", "equivalent_beam_angle"),
+    "AZFP": ("EL", "DS", "TVR", "VTX", "equivalent_beam_angle", "Sv_offset"),
+}
+
+
+# TODO: need a function (something like "_check_param_freq_dep")
+# to check user input cal_params and env_params
+
+
+def get_cal_params_AZFP(echodata: EchoData, user_cal_dict: dict) -> dict:
+    """
+    Get cal params using user inputs or values from data file.
+
+    Parameters
+    ----------
+    echodata : EchoData
+        An EchoData object containing data to be calibrated
+    user_cal_dict : dict
+        A dictionary containing user-defined calibration parameters.
+        The user-defined calibration parameters will overwrite values in the data file.
+
+    Returns
+    -------
+    A dict containing the calibration parameters for the AZFP echosounder
+    """
+    out_dict = dict.fromkeys(CAL_PARAMS["AZFP"])
+    if user_cal_dict is None:
+        user_cal_dict = {}
+
+    # Get params from Beam_group1
+    out_dict["equivalent_beam_angle"] = (
+        user_cal_dict["equivalent_beam_angle"]
+        if "equivalent_beam_angle" in user_cal_dict
+        else echodata["Sonar/Beam_group1"]["equivalent_beam_angle"]
+    )
+
+    # Get params from the Vendor_specific group
+    for p in ["EL", "DS", "TVR", "VTX", "Sv_offset"]:
+        # substitute if None in user input
+        out_dict[p] = user_cal_dict[p] if p in user_cal_dict else echodata["Vendor_specific"][p]
+
+    return out_dict
+
+
+def get_cal_params_EK(
+    beam: xr.Dataset, vend: xr.Dataset, user_cal_dict: Dict[str, xr.DataArray]
+) -> Dict:
+    """
+    Get cal params using user inputs or values from data file.
+
+    Parameters
+    ----------
+    beam : xr.Dataset
+        A subset of Sonar/Beam_groupX that contains only the channels specified for calibration
+    vend : xr.Dataset
+        A subset of Vendor_specific that contains only the channels specified for calibration
+    user_cal_dict : dict
+        A dictionary containing user-defined calibration parameters.
+        The user-defined calibration parameters will overwrite values in the data file.
+
+    Returns
+    -------
+    A dictionary containing the calibration parameters for the EK echosounders
+    """
+    out_dict = dict.fromkeys(CAL_PARAMS["EK"])
+    if user_cal_dict is None:
+        user_cal_dict = {}
+
+    # Params from the Vendor_specific group
+    params_from_vend = ["sa_correction", "gain_correction"]
+    for p in params_from_vend:
+        # substitute if p not in user input
+        out_dict[p] = (
+            user_cal_dict[p]
+            if p in user_cal_dict
+            else get_vend_cal_params_power(beam=beam, vend=vend, param=p)
+        )
+
+    # Other params
+    out_dict["equivalent_beam_angle"] = (
+        user_cal_dict["equivalent_beam_angle"]
+        if "equivalent_beam_angle" in user_cal_dict
+        else beam["equivalent_beam_angle"]
+    )
+
+    return out_dict
+
+
+def get_vend_filter_EK80(
+    vend: xr.Dataset, channel_id: str, filter_name: str, param_type: str
+) -> Union[np.ndarray, int]:
+    """
+    Get filter coefficients stored in the Vendor_specific group attributes.
+
+    Parameters
+    ----------
+    vend: xr.Dataset
+        EchoData["Vendor_specific"]
+    channel_id : str
+        channel id for which the param to be retrieved
+    filter_name : str
+        name of filter coefficients to retrieve
+    param_type : str
+        'coeff' or 'decimation'
+
+    Returns
+    -------
+    np.ndarray or int
+        The filter coefficient or the decimation factor
+    """
+    if param_type == "coeff":
+        v = vend.attrs["%s %s filter_r" % (channel_id, filter_name)] + 1j * np.array(
+            vend.attrs["%s %s filter_i" % (channel_id, filter_name)]
+        )
+        if v.size == 1:
+            v = np.expand_dims(v, axis=0)  # expand dims for convolution
+        return v
+    else:
+        return vend.attrs["%s %s decimation" % (channel_id, filter_name)]
+
+
+def get_vend_cal_params_power(beam: xr.Dataset, vend: xr.Dataset, param: str) -> xr.DataArray:
+    """
+    Get cal parameters stored in the Vendor_specific group
+    by matching the transmit_duration_nominal with allowable pulse_length.
+
+    Parameters
+    ----------
+    beam : xr.Dataset
+        A subset of Sonar/Beam_groupX that contains only the channels specified for calibration
+    vend : xr.Dataset
+        A subset of Vendor_specific that contains only the channels specified for calibration
+    param : str {"sa_correction", "gain_correction"}
+        name of parameter to retrieve
+
+    Returns
+    -------
+    An xr.DataArray containing the matched ``param``
+    """
+
+    # Check parameter is among those allowed
+    if param not in ["sa_correction", "gain_correction"]:
+        raise ValueError(f"Unknown parameter {param}")
+
+    # Check parameter exists
+    if param not in vend:
+        raise ValueError(f"{param} does not exist in the Vendor_specific group!")
+
+    # Find idx to select the corresponding param value
+    # by matching beam["transmit_duration_nominal"] with ds_vend["pulse_length"]
+    transmit_isnull = beam["transmit_duration_nominal"].isnull()
+    idxmin = np.abs(beam["transmit_duration_nominal"] - vend["pulse_length"]).idxmin(
+        dim="pulse_length_bin"
+    )
+
+    # fill nan position with 0 (will remove before return)
+    # and convert to int for indexing
+    idxmin = idxmin.where(~transmit_isnull, 0).astype(int)
+
+    # Get param dataarray into correct shape
+    da_param = (
+        vend[param]
+        .expand_dims(dim={"ping_time": idxmin["ping_time"]})  # expand dims for direct indexing
+        .sortby(idxmin.channel)  # sortby in case channel sequence differs in vend and beam
+    )
+
+    # Select corresponding index and clean up the original nan elements
+    da_param = da_param.sel(pulse_length_bin=idxmin, drop=True)
+    return da_param.where(~transmit_isnull, np.nan)  # set the nan elements back to nan
+
+
+def get_gain_BB(
+    vend: xr.Dataset, freq_center: xr.DataArray, cal_params_CW: Dict[str, xr.DataArray]
+) -> xr.DataArray:
+    """
+    Get broadband gain factor for calibrating complex samples.
+
+    Interpolate ``gain`` in the Vendor_specific group to the center frequency of each ping
+    for BB mode samples if nominal frequency is within the calibrated frequency range
+
+    Parameters
+    ----------
+    waveform_mode : str
+        ``CW`` for CW-mode samples, either recorded as complex or power samples
+        ``BB`` for BB-mode samples, recorded as complex samples
+    chan_sel : xr.DataArray
+        Nominal channel for CW mode samples
+        and a xr.DataArray of selected channels for BB mode samples
+
+    Returns
+    -------
+    An xr.DataArray containing the interpolated gain factors
+    """
+    gain = []
+    for ch_id in freq_center["channel"]:
+        # if frequency-dependent gain exists in Vendor group, interpolate at center frequency
+        if ch_id in vend["cal_channel_id"]:
+            gain_temp = (
+                vend["gain"]
+                .sel(cal_channel_id=ch_id)
+                .interp(cal_frequency=freq_center.sel(channel=ch_id))
+                .drop(["cal_channel_id", "cal_frequency"])
+                .expand_dims("channel")
+            )
+        # if no frequency-dependent gain exists, use CW gain
+        else:
+            gain_temp = (
+                cal_params_CW["gain_correction"].sel(channel=ch_id)
+                # .reindex_like(echodata["Sonar/Beam_group1"]["backscatter_r"], method="nearest")
+                .expand_dims("channel")
+            )
+        gain_temp.name = "gain"
+        gain.append(gain_temp)
+    gain = xr.merge(gain).gain  # select the single data variable
+
+    return gain
diff --git a/echopype/calibrate/calibrate_azfp.py b/echopype/calibrate/calibrate_azfp.py
index b61c79c7f..10ba8b793 100644
--- a/echopype/calibrate/calibrate_azfp.py
+++ b/echopype/calibrate/calibrate_azfp.py
@@ -1,86 +1,24 @@
 import numpy as np
 
 from ..echodata import EchoData
-from ..utils import uwa
-from .calibrate_base import CAL_PARAMS
+from .cal_params import get_cal_params_AZFP
 from .calibrate_ek import CalibrateBase
+from .env_params import get_env_params_AZFP
+from .range import compute_range_AZFP
 
 
 class CalibrateAZFP(CalibrateBase):
-    def __init__(self, echodata: EchoData, env_params, cal_params, **kwargs):
+    def __init__(self, echodata: EchoData, env_params=None, cal_params=None, **kwargs):
         super().__init__(echodata, env_params)
 
-        # initialize cal params
-        self.cal_params = dict.fromkeys(CAL_PARAMS["AZFP"])
-
         # load env and cal parameters
-        self.get_env_params()
-        if cal_params is None:
-            cal_params = {}
-        self.get_cal_params(cal_params)
+        self.env_params = get_env_params_AZFP(echodata=echodata, user_env_dict=env_params)
+        self.cal_params = get_cal_params_AZFP(echodata=echodata, user_cal_dict=cal_params)
 
         # self.range_meter computed under self._cal_power()
         # because the implementation is different for Sv and TS
 
-    def get_cal_params(self, cal_params):
-        """Get cal params using user inputs or values from data file.
-
-        Parameters
-        ----------
-        cal_params : dict
-        """
-
-        # Get params from Beam_group1
-        self.cal_params["equivalent_beam_angle"] = (
-            cal_params["equivalent_beam_angle"]
-            if "equivalent_beam_angle" in cal_params
-            else self.echodata["Sonar/Beam_group1"]["equivalent_beam_angle"]
-        )
-
-        # Get params from the Vendor_specific group
-        for p in ["EL", "DS", "TVR", "VTX", "Sv_offset"]:
-            # substitute if None in user input
-            self.cal_params[p] = (
-                cal_params[p] if p in cal_params else self.echodata["Vendor_specific"][p]
-            )
-
-    def get_env_params(self):
-        """Get env params using user inputs or values from data file.
-
-        Parameters
-        ----------
-        env_params : dict
-        """
-        # Temperature comes from either user input or data file
-        self.env_params["temperature"] = (
-            self.env_params["temperature"]
-            if "temperature" in self.env_params
-            else self.echodata["Environment"]["temperature"]
-        )
-
-        # Salinity and pressure always come from user input
-        if ("salinity" not in self.env_params) or ("pressure" not in self.env_params):
-            raise ReferenceError("Please supply both salinity and pressure in env_params.")
-        else:
-            self.env_params["salinity"] = self.env_params["salinity"]
-            self.env_params["pressure"] = self.env_params["pressure"]
-
-        # Always calculate sound speed and absorption
-        self.env_params["sound_speed"] = uwa.calc_sound_speed(
-            temperature=self.env_params["temperature"],
-            salinity=self.env_params["salinity"],
-            pressure=self.env_params["pressure"],
-            formula_source="AZFP",
-        )
-        self.env_params["sound_absorption"] = uwa.calc_absorption(
-            frequency=self.echodata["Sonar/Beam_group1"]["frequency_nominal"],
-            temperature=self.env_params["temperature"],
-            salinity=self.env_params["salinity"],
-            pressure=self.env_params["pressure"],
-            formula_source="AZFP",
-        )
-
-    def compute_range_meter(self, cal_type):
+    def compute_echo_range(self, cal_type):
         """Calculate range (``echo_range``) in meter using AZFP formula.
 
         Note the range calculation differs for Sv and TS per AZFP matlab code.
@@ -91,9 +29,11 @@ def compute_range_meter(self, cal_type):
             'Sv' for calculating volume backscattering strength, or
             'TS' for calculating target strength
         """
-        self.range_meter = self.echodata.compute_range(self.env_params, azfp_cal_type=cal_type)
+        self.range_meter = compute_range_AZFP(
+            echodata=self.echodata, env_params=self.env_params, cal_type=cal_type
+        )
 
-    def _cal_power(self, cal_type, **kwargs):
+    def _cal_power_samples(self, cal_type, **kwargs):
         """Calibrate to get volume backscattering strength (Sv) from AZFP power data.
 
         The calibration formulae used here is based on Appendix G in
@@ -103,9 +43,8 @@ def _cal_power(self, cal_type, **kwargs):
         See calc_Sv_offset() in convert/azfp.py
         """
         # Compute range in meters
-        self.compute_range_meter(
-            cal_type=cal_type
-        )  # range computation different for Sv and TS per AZFP matlab code
+        # range computation different for Sv and TS per AZFP matlab code
+        self.compute_echo_range(cal_type=cal_type)
 
         # Compute derived params
 
@@ -172,7 +111,7 @@ def _cal_power(self, cal_type, **kwargs):
         return out.squeeze("beam", drop=True)
 
     def compute_Sv(self, **kwargs):
-        return self._cal_power(cal_type="Sv")
+        return self._cal_power_samples(cal_type="Sv")
 
     def compute_TS(self, **kwargs):
-        return self._cal_power(cal_type="TS")
+        return self._cal_power_samples(cal_type="TS")
diff --git a/echopype/calibrate/calibrate_base.py b/echopype/calibrate/calibrate_base.py
index ae47eccc8..cbd911bf8 100644
--- a/echopype/calibrate/calibrate_base.py
+++ b/echopype/calibrate/calibrate_base.py
@@ -1,7 +1,7 @@
 import abc
 
 from ..echodata import EchoData
-from .env_params import EnvParams
+from .env_params_old import EnvParams
 
 CAL_PARAMS = {
     "EK": ("sa_correction", "gain_correction", "equivalent_beam_angle"),
@@ -12,8 +12,10 @@
 class CalibrateBase(abc.ABC):
     """Class to handle calibration for all sonar models."""
 
-    def __init__(self, echodata: EchoData, env_params=None):
+    def __init__(self, echodata: EchoData, env_params=None, cal_params=None):
         self.echodata = echodata
+
+        # TODO: clean up init of env_params here: can just leave it as None
         if isinstance(env_params, EnvParams):
             env_params = env_params._apply(echodata)
         elif env_params is None:
@@ -22,22 +24,14 @@ def __init__(self, echodata: EchoData, env_params=None):
             raise ValueError(
                 "invalid env_params type; provide an EnvParams instance, a dict, or None"
             )
-        self.env_params = env_params  # env_params are set in child class
-        self.cal_params = None  # cal_params are set in child class
+        self.env_params = env_params
+        self.cal_params = cal_params
 
         # range_meter is computed in compute_Sv/TS in child class
         self.range_meter = None
 
     @abc.abstractmethod
-    def get_env_params(self, **kwargs):
-        pass
-
-    @abc.abstractmethod
-    def get_cal_params(self, **kwargs):
-        pass
-
-    @abc.abstractmethod
-    def compute_range_meter(self, **kwargs):
+    def compute_echo_range(self, **kwargs):
         """Calculate range (``echo_range``) in units meter.
 
         Returns
@@ -48,7 +42,7 @@ def compute_range_meter(self, **kwargs):
         pass
 
     @abc.abstractmethod
-    def _cal_power(self, cal_type, **kwargs):
+    def _cal_power_samples(self, cal_type, **kwargs):
         """Calibrate power data for EK60, EK80, and AZFP.
 
         Parameters
diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py
index 11b03d7e2..986a40714 100644
--- a/echopype/calibrate/calibrate_ek.py
+++ b/echopype/calibrate/calibrate_ek.py
@@ -1,141 +1,42 @@
+from collections import defaultdict
+from typing import Dict
+
 import numpy as np
 import xarray as xr
-from scipy import signal
 
 from ..echodata import EchoData
-from ..echodata.simrad import retrieve_correct_beam_group
-from ..utils import uwa
+from ..echodata.simrad import check_input_args_combination, retrieve_correct_beam_group
 from ..utils.log import _init_logger
-from .calibrate_base import CAL_PARAMS, CalibrateBase
+from .cal_params import get_cal_params_EK, get_gain_BB, get_vend_filter_EK80
+from .calibrate_base import CalibrateBase
+from .ek80_complex import compress_pulse, get_tau_effective, get_transmit_signal
+from .env_params import get_env_params_EK60, get_env_params_EK80
+from .range import compute_range_EK
 
 logger = _init_logger(__name__)
 
 
 class CalibrateEK(CalibrateBase):
-    def __init__(self, echodata: EchoData, env_params):
-        super().__init__(echodata, env_params)
-
-        # cal params specific to EK echosounders
-        self.cal_params = dict.fromkeys(CAL_PARAMS["EK"])
+    def __init__(self, echodata: EchoData, env_params, cal_params):
+        super().__init__(echodata, env_params, cal_params)
 
-    def compute_range_meter(self, waveform_mode, encode_mode):
+    def compute_echo_range(self):
         """
-        Parameters
-        ----------
-        waveform_mode : {"CW", "BB"}
-            Type of transmit waveform.
-            Required only for data from the EK80 echosounder.
-
-            - `"CW"` for narrowband transmission,
-              returned echoes recorded either as complex or power/angle samples
-            - `"BB"` for broadband transmission,
-              returned echoes recorded as complex samples
-
-        encode_mode : {"complex", "power"}
-            Type of encoded return echo data.
-            Required only for data from the EK80 echosounder.
-
-            - `"complex"` for complex samples
-            - `"power"` for power/angle samples, only allowed when
-              the echosounder is configured for narrowband transmission
+        Compute echo range for EK echosounders.
 
         Returns
         -------
         range_meter : xr.DataArray
             range in units meter
         """
-        self.range_meter = self.echodata.compute_range(
-            self.env_params,
-            ek_waveform_mode=waveform_mode,
-            ek_encode_mode=encode_mode,
-        )
-
-    def _get_vend_cal_params_power(self, param, waveform_mode):
-        """Get cal parameters stored in the Vendor_specific group.
-
-        Parameters
-        ----------
-        param : str {"sa_correction", "gain_correction"}
-            name of parameter to retrieve
-        """
-        ds_vend = self.echodata["Vendor_specific"]
-
-        if ds_vend is None or param not in ds_vend:
-            return None
-
-        if param not in ["sa_correction", "gain_correction"]:
-            raise ValueError(f"Unknown parameter {param}")
-
-        if waveform_mode == "CW" and self.echodata["Sonar/Beam_group2"] is not None:
-            beam = self.echodata["Sonar/Beam_group2"]
-        else:
-            beam = self.echodata["Sonar/Beam_group1"]
-
-        # indexes of frequencies that are for power, not complex
-        relevant_indexes = np.where(
-            np.isin(ds_vend["frequency_nominal"], beam["frequency_nominal"])
-        )[0]
-
-        # Find idx to select the corresponding param value
-        # by matching beam["transmit_duration_nominal"] with ds_vend["pulse_length"]
-        transmit_isnull = beam["transmit_duration_nominal"].isnull()
-        idxmin = np.abs(
-            beam["transmit_duration_nominal"] - ds_vend["pulse_length"][relevant_indexes]
-        ).idxmin(dim="pulse_length_bin")
-
-        # fill nan position with 0 (witll remove before return)
-        # and convert to int for indexing
-        idxmin = idxmin.where(~transmit_isnull, 0).astype(int)
-
-        # Get param dataarray into correct shape
-        da_param = (
-            ds_vend[param][relevant_indexes]
-            .expand_dims(dim={"ping_time": idxmin["ping_time"]})  # expand dims for direct indexing
-            .sortby(idxmin.channel)  # sortby in case channel sequence different in vendor and beam
-        )
-
-        # Select corresponding index and clean up the original nan elements
-        da_param = da_param.sel(pulse_length_bin=idxmin, drop=True)
-        return da_param.where(~transmit_isnull, np.nan)  # set the nan elements back to nan
-
-    def get_cal_params(self, cal_params, waveform_mode, encode_mode):
-        """Get cal params using user inputs or values from data file.
-
-        Parameters
-        ----------
-        cal_params : dict
-        """
-
-        if (
-            encode_mode == "power"
-            and waveform_mode == "CW"
-            and self.echodata["Sonar/Beam_group2"] is not None
-        ):
-            beam = self.echodata["Sonar/Beam_group2"]
-        else:
-            beam = self.echodata["Sonar/Beam_group1"]
-
-        # Params from the Vendor_specific group
-
-        # only execute this if cw and power
-        if waveform_mode == "CW" and beam is not None:
-            params_from_vend = ["sa_correction", "gain_correction"]
-            for p in params_from_vend:
-                # substitute if None in user input
-                self.cal_params[p] = (
-                    cal_params[p]
-                    if p in cal_params
-                    else self._get_vend_cal_params_power(p, waveform_mode=waveform_mode)
-                )
-
-        # Other params
-        self.cal_params["equivalent_beam_angle"] = (
-            cal_params["equivalent_beam_angle"]
-            if "equivalent_beam_angle" in cal_params
-            else beam["equivalent_beam_angle"]
+        self.range_meter = compute_range_EK(
+            echodata=self.echodata,
+            env_params=self.env_params,
+            waveform_mode=self.waveform_mode,
+            encode_mode=self.encode_mode,
         )
 
-    def _cal_power(self, cal_type: str, power_ed_group: str = None) -> xr.Dataset:
+    def _cal_power_samples(self, cal_type: str, power_ed_group: str = None) -> xr.Dataset:
         """Calibrate power data from EK60 and EK80.
 
         Parameters
@@ -226,450 +127,182 @@ def _cal_power(self, cal_type: str, power_ed_group: str = None) -> xr.Dataset:
 
 class CalibrateEK60(CalibrateEK):
     def __init__(self, echodata, env_params, cal_params, **kwargs):
-        super().__init__(echodata, env_params)
+        super().__init__(echodata, env_params, cal_params)
 
-        # load env and cal parameters
-        self.get_env_params()
-        if cal_params is None:
-            cal_params = {}
-        self.get_cal_params(cal_params, waveform_mode="CW", encode_mode="power")
+        # Get env_params
+        self.env_params = get_env_params_EK60(echodata=echodata, user_env_dict=env_params)
+        self.waveform_mode = "CW"
+        self.encode_mode = "power"
 
-        # default to CW mode recorded as power samples
-        self.compute_range_meter(waveform_mode="CW", encode_mode="power")
+        # Compute range
+        self.compute_echo_range()
 
-    def get_env_params(self, **kwargs):
-        """Get env params using user inputs or values from data file.
+        # Get the right ed_group for CW power samples
+        self.ed_group = retrieve_correct_beam_group(
+            echodata=self.echodata, waveform_mode=self.waveform_mode, encode_mode=self.encode_mode
+        )
 
-        EK60 file by default contains only sound speed and absorption.
-        In cases when temperature, salinity, and pressure values are supplied
-        by the user simultaneously, the sound speed and absorption are re-calculated.
+        # Set the channels to calibrate
+        # For EK60 this is all channels
+        self.chan_sel = self.echodata[self.ed_group]["channel"]
 
-        Parameters
-        ----------
-        env_params : dict
-        """
-        # Re-calculate environment parameters if user supply all env variables
-        if (
-            ("temperature" in self.env_params)
-            and ("salinity" in self.env_params)
-            and ("pressure" in self.env_params)
-        ):
-            self.env_params["sound_speed"] = uwa.calc_sound_speed(
-                temperature=self.env_params["temperature"],
-                salinity=self.env_params["salinity"],
-                pressure=self.env_params["pressure"],
-            )
-            self.env_params["sound_absorption"] = uwa.calc_absorption(
-                frequency=self.echodata["Sonar/Beam_group1"]["frequency_nominal"],
-                temperature=self.env_params["temperature"],
-                salinity=self.env_params["salinity"],
-                pressure=self.env_params["pressure"],
-            )
-        # Otherwise get sound speed and absorption from user inputs or raw data file
-        else:
-            self.env_params["sound_speed"] = (
-                self.env_params["sound_speed"]
-                if "sound_speed" in self.env_params
-                else self.echodata["Environment"]["sound_speed_indicative"]
-            )
-            self.env_params["sound_absorption"] = (
-                self.env_params["sound_absorption"]
-                if "sound_absorption" in self.env_params
-                else self.echodata["Environment"]["absorption_indicative"]
-            )
+        # Get cal_params
+        self.cal_params = get_cal_params_EK(
+            beam=echodata[self.ed_group].sel(channel=self.chan_sel),
+            vend=echodata["Vendor_specific"].sel(channel=self.chan_sel),
+            user_cal_dict=cal_params,
+        )
 
     def compute_Sv(self, **kwargs):
-        power_ed_group = retrieve_correct_beam_group(
-            echodata=self.echodata, waveform_mode="CW", encode_mode="power", pulse_compression=False
-        )
-        return self._cal_power(cal_type="Sv", power_ed_group=power_ed_group)
+        return self._cal_power_samples(cal_type="Sv", power_ed_group=self.ed_group)
 
     def compute_TS(self, **kwargs):
-        power_ed_group = retrieve_correct_beam_group(
-            echodata=self.echodata, waveform_mode="CW", encode_mode="power", pulse_compression=False
-        )
-        return self._cal_power(cal_type="TS", power_ed_group=power_ed_group)
+        return self._cal_power_samples(cal_type="TS", power_ed_group=self.ed_group)
 
 
 class CalibrateEK80(CalibrateEK):
+
+    # TODO: add option to get these from data file
     fs = 1.5e6  # default full sampling frequency [Hz]
     z_et = 75
     z_er = 1000
 
     def __init__(self, echodata, env_params, cal_params, waveform_mode, encode_mode):
-        super().__init__(echodata, env_params)
-
-        # initialize cal params
-        # cal params are those used by both complex and power data calibration
-        # TODO: add complex data-specific params, like the freq-dependent gain factor
-        self.cal_params = dict.fromkeys(CAL_PARAMS["EK"])
-        # TODO: make waveform_mode and encode_mode class attributes
+        super().__init__(echodata, env_params, cal_params)
 
-        # load env and cal parameters
-        self.get_env_params(waveform_mode=waveform_mode, encode_mode=encode_mode)
-        if cal_params is None:
-            cal_params = {}
-        self.get_cal_params(cal_params, waveform_mode=waveform_mode, encode_mode=encode_mode)
+        # Check the combination of waveform and encode mode makes sense
+        check_input_args_combination(waveform_mode, encode_mode)
+        self.waveform_mode = waveform_mode
+        self.encode_mode = encode_mode
+        self.echodata = echodata
 
-        # self.range_meter computed under self._compute_cal()
-        # because the implementation is different depending on waveform_mode and encode_mode
+        # Get the right ed_group given waveform and encode mode
+        self.ed_group = retrieve_correct_beam_group(
+            echodata=self.echodata, waveform_mode=self.waveform_mode, encode_mode=self.encode_mode
+        )
 
-    def get_env_params(self, waveform_mode=None, encode_mode="complex"):
-        """Get env params using user inputs or values from data file.
+        # Select the channels to calibrate
+        if self.encode_mode == "power":
+            # Power sample only possible under CW mode,
+            # and all power samples will live in the same group
+            self.chan_sel = self.echodata[self.ed_group]["channel"]
+        else:
+            # Complex samples can be CW or BB, so select based on waveform mode
+            chan_dict = self._get_chan_dict(self.echodata[self.ed_group])
+            self.chan_sel = chan_dict[self.waveform_mode]
 
-        EK80 file by default contains sound speed, temperature, depth, salinity, and acidity,
-        therefore absorption is always calculated unless it is supplied by the user.
-        In cases when temperature, salinity, and pressure values are supplied
-        by the user simultaneously, both the sound speed and absorption are re-calculated.
+        # Subset of the right Sonar/Beam_groupX group given the selected channels
+        beam = self.echodata[self.ed_group].sel(channel=self.chan_sel)
 
-        Parameters
-        ----------
-        env_params : dict
+        # Use center frequency if in BB mode, else use nominal channel frequency
+        if self.waveform_mode == "BB":
+            # use true center frequency to interpolate for gain factor
+            self.freq_center = (beam["frequency_start"] + beam["frequency_end"]).sel(
+                channel=self.chan_sel
+            ).isel(beam=0).drop("beam") / 2
+        else:
+            # use nominal channel frequency for CW pulse
+            self.freq_center = beam["frequency_nominal"].sel(channel=self.chan_sel)
 
-        waveform_mode : {"CW", "BB"}
-            Type of transmit waveform.
+        # Get env_params: depends on waveform mode
+        self.env_params = get_env_params_EK80(
+            echodata=echodata, freq=self.freq_center, user_env_dict=env_params
+        )
 
-            - `"CW"` for narrowband transmission,
-              returned echoes recorded either as complex or power/angle samples
-            - (default) `"BB"` for broadband transmission,
-              returned echoes recorded as complex samples
+        # Get cal_params: depends on waveform and encode mode
+        self.cal_params = get_cal_params_EK(
+            beam=beam,
+            vend=self.echodata["Vendor_specific"].sel(channel=self.chan_sel),
+            user_cal_dict=cal_params,
+        )
 
-        encode_mode : {"complex", "power"}
-            Type of encoded return echo data.
+        # Compute echo range in meters
+        self.compute_echo_range()
 
-            - (default) `"complex"` for complex samples
-            - `"power"` for power/angle samples, only allowed when
-              the echosounder is configured for narrowband transmission
+    @staticmethod
+    def _get_chan_dict(beam: xr.Dataset) -> Dict:
         """
+        Build dict to select BB and CW channels from complex samples where data
+        from both waveform modes may co-exist.
+        """
+        # Use center frequency for each ping to select BB or CW channels
+        # when all samples are encoded as complex samples
+        if "frequency_start" in beam and "frequency_end" in beam:
+            # At least some channels are BB
+            # frequency_start and frequency_end are NaN for CW channels
+            freq_center = (beam["frequency_start"] + beam["frequency_end"]) / 2  # has beam dim
+
+            return {
+                # For BB: drop channels containing CW samples (nan in freq start/end)
+                "BB": freq_center.dropna(dim="channel").channel,
+                # For CW: drop channels containing BB samples (not nan in freq start/end)
+                "CW": freq_center.where(np.isnan(freq_center), drop=True).channel,
+            }
 
-        if (
-            encode_mode == "power"
-            and waveform_mode == "CW"
-            and self.echodata["Sonar/Beam_group2"] is not None
-        ):
-            beam = self.echodata["Sonar/Beam_group2"]
         else:
-            beam = self.echodata["Sonar/Beam_group1"]
+            # All channels are CW
+            return {"BB": None, "CW": beam.channel}
 
-        # Use center frequency if in BB mode, else use nominal channel frequency
-        if waveform_mode == "BB":
-            freq = (beam["frequency_start"] + beam["frequency_end"]) / 2
-        else:
-            freq = beam["frequency_nominal"]
-
-        # Re-calculate environment parameters if user supply all env variables
-        if (
-            ("temperature" in self.env_params)
-            and ("salinity" in self.env_params)
-            and ("pressure" in self.env_params)
-        ):
-            self.env_params["sound_speed"] = uwa.calc_sound_speed(
-                temperature=self.env_params["temperature"],
-                salinity=self.env_params["salinity"],
-                pressure=self.env_params["pressure"],
-            )
-            self.env_params["sound_absorption"] = uwa.calc_absorption(
-                frequency=freq,
-                temperature=self.env_params["temperature"],
-                salinity=self.env_params["salinity"],
-                pressure=self.env_params["pressure"],
-            )
-        # Otherwise
-        #  get temperature, salinity, and pressure from raw data file
-        #  get sound speed from user inputs or raw data file
-        #  get absorption from user inputs or computing from env params stored in raw data file
-        else:
-            # pressure is encoded as "depth" in EK80  # TODO: change depth to pressure in EK80 file?
-            for p1, p2 in zip(
-                ["temperature", "salinity", "pressure"],
-                ["temperature", "salinity", "depth"],
-            ):
-                self.env_params[p1] = (
-                    self.env_params[p1]
-                    if p1 in self.env_params
-                    else self.echodata["Environment"][p2]
-                )
-            self.env_params["sound_speed"] = (
-                self.env_params["sound_speed"]
-                if "sound_speed" in self.env_params
-                else self.echodata["Environment"]["sound_speed_indicative"]
-            )
-            self.env_params["sound_absorption"] = (
-                self.env_params["sound_absorption"]
-                if "sound_absorption" in self.env_params
-                else uwa.calc_absorption(
-                    frequency=freq,
-                    temperature=self.env_params["temperature"],
-                    salinity=self.env_params["salinity"],
-                    pressure=self.env_params["pressure"],
-                )
-            )
-
-    def _get_vend_cal_params_complex(self, channel_id, filter_name, param_type):
-        """Get filter coefficients stored in the Vendor_specific group attributes.
-
-        Parameters
-        ----------
-        channel_id : str
-            channel id for which the param to be retrieved
-        filter_name : str
-            name of filter coefficients to retrieve
-        param_type : str
-            'coeff' or 'decimation'
+    def _get_filter_coeff(self, channel: xr.DataArray) -> Dict:
         """
-        if param_type == "coeff":
-            v = self.echodata["Vendor_specific"].attrs[
-                "%s %s filter_r" % (channel_id, filter_name)
-            ] + 1j * np.array(
-                self.echodata["Vendor_specific"].attrs["%s %s filter_i" % (channel_id, filter_name)]
-            )
-            if v.size == 1:
-                v = np.expand_dims(v, axis=0)  # expand dims for convolution
-            return v
-        else:
-            return self.echodata["Vendor_specific"].attrs[
-                "%s %s decimation" % (channel_id, filter_name)
-            ]
-
-    def _tapered_chirp(
-        self,
-        transmit_duration_nominal,
-        slope,
-        transmit_power,
-        frequency_nominal=None,
-        frequency_start=None,
-        frequency_end=None,
-    ):
-        """Create a baseline chirp template."""
-        if frequency_start is None and frequency_end is None:  # CW waveform
-            frequency_start = frequency_nominal
-            frequency_end = frequency_nominal
-
-        t = np.arange(0, transmit_duration_nominal, 1 / self.fs)
-        nwtx = int(2 * np.floor(slope * t.size))  # length of tapering window
-        wtx_tmp = np.hanning(nwtx)  # hanning window
-        nwtxh = int(np.round(nwtx / 2))  # half length of the hanning window
-        wtx = np.concatenate(
-            [wtx_tmp[0:nwtxh], np.ones((t.size - nwtx)), wtx_tmp[nwtxh:]]
-        )  # assemble full tapering window
-        y_tmp = (
-            np.sqrt((transmit_power / 4) * (2 * self.z_et))  # amplitude
-            * signal.chirp(t, frequency_start, t[-1], frequency_end)
-            * wtx
-        )  # taper and scale linear chirp
-        return y_tmp / np.max(np.abs(y_tmp)), t  # amp has no actual effect
-
-    def _filter_decimate_chirp(self, y, ch_id):
-        """Filter and decimate the chirp template.
+        Get WBT and PC filter coefficients for constructing the transmit replica.
 
-        Parameters
-        ----------
-        y : np.array
-            chirp from _tapered_chirp
-        ch_id : str
-            channel_id to select the right coefficients and factors
+        Returns
+        -------
+        A dictionary indexed by ``channel`` and values being dictionaries containing
+        filter coefficients and decimation factors for constructing the transmit replica.
         """
-        # filter coefficients and decimation factor
-        wbt_fil = self._get_vend_cal_params_complex(ch_id, "WBT", "coeff")
-        pc_fil = self._get_vend_cal_params_complex(ch_id, "PC", "coeff")
-        wbt_decifac = self._get_vend_cal_params_complex(ch_id, "WBT", "decimation")
-        pc_decifac = self._get_vend_cal_params_complex(ch_id, "PC", "decimation")
-
-        # WBT filter and decimation
-        ytx_wbt = signal.convolve(y, wbt_fil)
-        ytx_wbt_deci = ytx_wbt[0::wbt_decifac]
-
-        # PC filter and decimation
-        if len(pc_fil.squeeze().shape) == 0:  # in case it is a single element
-            pc_fil = [pc_fil.squeeze()]
-        ytx_pc = signal.convolve(ytx_wbt_deci, pc_fil)
-        ytx_pc_deci = ytx_pc[0::pc_decifac]
-        ytx_pc_deci_time = np.arange(ytx_pc_deci.size) * 1 / self.fs * wbt_decifac * pc_decifac
-
-        return ytx_pc_deci, ytx_pc_deci_time
-
-    @staticmethod
-    def _get_tau_effective(ytx, fs_deci, waveform_mode):
-        """Compute effective pulse length.
-
-        Parameters
-        ----------
-        ytx :array
-            transmit signal
-        fs_deci : float
-            sampling frequency of the decimated (recorded) signal
-        waveform_mode : str
-            ``CW`` for CW-mode samples, either recorded as complex or power samples
-            ``BB`` for BB-mode samples, recorded as complex samples
+        vend = self.echodata["Vendor_specific"]
+        coeff = defaultdict(dict)
+        for ch_id in channel.values:
+            # filter coefficients and decimation factor
+            coeff[ch_id]["wbt_fil"] = get_vend_filter_EK80(vend, ch_id, "WBT", "coeff")
+            coeff[ch_id]["pc_fil"] = get_vend_filter_EK80(vend, ch_id, "PC", "coeff")
+            coeff[ch_id]["wbt_decifac"] = get_vend_filter_EK80(vend, ch_id, "WBT", "decimation")
+            coeff[ch_id]["pc_decifac"] = get_vend_filter_EK80(vend, ch_id, "PC", "decimation")
+
+        return coeff
+
+    def _get_power_from_complex(
+        self, beam: xr.Dataset, chan_sel: xr.DataArray, chirp: Dict
+    ) -> xr.DataArray:
         """
-        if waveform_mode == "BB":
-            ytxa = signal.convolve(ytx, np.flip(np.conj(ytx))) / np.linalg.norm(ytx) ** 2
-            ptxa = abs(ytxa) ** 2
-        elif waveform_mode == "CW":
-            ptxa = np.abs(ytx) ** 2  # energy of transmit signal
-        return ptxa.sum() / (
-            ptxa.max() * fs_deci
-        )  # TODO: verify fs_deci = 1.5e6 in spheroid data sets
-
-    def get_transmit_chirp(self, waveform_mode):
-        """Reconstruct transmit signal and compute effective pulse length.
+        Get power from complex samples.
 
         Parameters
         ----------
-        waveform_mode : str
-            ``CW`` for CW-mode samples, either recorded as complex or power samples
-            ``BB`` for BB-mode samples, recorded as complex samples
-        """
-        # Make sure it is BB mode data
-        if waveform_mode == "BB" and (
-            ("frequency_start" not in self.echodata["Sonar/Beam_group1"])
-            or ("frequency_end" not in self.echodata["Sonar/Beam_group1"])
-        ):
-            raise TypeError("File does not contain BB mode complex samples!")
-
-        y_all = {}
-        y_time_all = {}
-        tau_effective = {}
-        for chan in self.echodata["Sonar/Beam_group1"].channel.values:
-            # TODO: currently only deal with the case with
-            # a fixed tx key param values within a channel
-            if waveform_mode == "BB":
-                tx_param_names = [
-                    "transmit_duration_nominal",
-                    "slope",
-                    "transmit_power",
-                    "frequency_start",
-                    "frequency_end",
-                ]
-            else:
-                tx_param_names = [
-                    "transmit_duration_nominal",
-                    "slope",
-                    "transmit_power",
-                    "frequency_nominal",
-                ]
-            tx_params = {}
-            for p in tx_param_names:
-                tx_params[p] = np.unique(self.echodata["Sonar/Beam_group1"][p].sel(channel=chan))
-                if tx_params[p].size != 1:
-                    raise TypeError("File contains changing %s!" % p)
-            y_tmp, _ = self._tapered_chirp(**tx_params)
-
-            # Filter and decimate chirp template
-            fs_deci = (
-                1 / self.echodata["Sonar/Beam_group1"].sel(channel=chan)["sample_interval"].values
-            )
-            y_tmp, y_tmp_time = self._filter_decimate_chirp(y_tmp, chan)
-
-            # Compute effective pulse length
-            tau_effective_tmp = self._get_tau_effective(y_tmp, fs_deci, waveform_mode=waveform_mode)
-
-            y_all[chan] = y_tmp
-            y_time_all[chan] = y_tmp_time
-            tau_effective[chan] = tau_effective_tmp
-
-        return y_all, y_time_all, tau_effective
+        beam : xr.Dataset
+            EchoData["Sonar/Beam_group1"]
 
-    def compress_pulse(self, chirp, chan_BB=None):
-        """Perform pulse compression on the backscatter data.
-
-        Parameters
-        ----------
-        chirp : dict
-            transmit chirp replica indexed by channel_id
-        chan_BB : str
-            channels that transmit in BB mode
-            (since CW mode can be in mixed in complex samples too)
+        Returns
+        -------
+        prx : xr.DataArray
+            Power computed from complex samples
         """
-        backscatter = self.echodata["Sonar/Beam_group1"]["backscatter_r"].sel(
-            channel=chan_BB
-        ) + 1j * self.echodata["Sonar/Beam_group1"]["backscatter_i"].sel(channel=chan_BB)
-
-        pc_all = []
-        for chan in chan_BB:
-            backscatter_chan = (
-                backscatter.sel(channel=chan)
-                .dropna(dim="range_sample", how="all")
-                .dropna(dim="beam", how="all")
-                .dropna(dim="ping_time")
-            )
-            replica = xr.DataArray(np.conj(chirp[str(chan.values)]), dims="window")
-            # Pulse compression via rolling
-            pc = (
-                backscatter_chan.rolling(range_sample=replica.size).construct("window").dot(replica)
-                / np.linalg.norm(chirp[str(chan.values)]) ** 2
-            )
-            # Expand dimension and add name to allow merge
-            pc = pc.expand_dims(dim="channel")
-            pc.name = "pulse_compressed_output"
-            pc_all.append(pc)
-
-        pc_merge = xr.merge(pc_all)
 
-        return pc_merge
+        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(self.z_er + self.z_et) / self.z_er) ** 2
+                / self.z_et
+            )
 
-    def _get_gain_for_complex(self, waveform_mode, chan_sel) -> xr.DataArray:
-        """Get gain factor for calibrating complex samples.
+        if self.waveform_mode == "BB":
+            pc = compress_pulse(beam=beam, chirp=chirp, chan_BB=chan_sel)  # has beam dim
+            prx = _get_prx(pc["pulse_compressed_output"])  # ensure prx is xr.DataArray
+        else:
+            bs_cw = beam["backscatter_r"].sel(channel=chan_sel) + 1j * beam["backscatter_i"].sel(
+                channel=chan_sel
+            )
+            prx = _get_prx(bs_cw)
 
-        Use values from ``gain_correction`` in the Vendor_specific group for CW mode samples,
-        or interpolate ``gain`` to the center frequency of each ping for BB mode samples
-        if nominal frequency is within the calibrated frequencies range
+        prx.name = "received_power"
 
-        Parameters
-        ----------
-        waveform_mode : str
-            ``CW`` for CW-mode samples, either recorded as complex or power samples
-            ``BB`` for BB-mode samples, recorded as complex samples
-        chan_sel : xr.DataArray
-            Nominal channel for CW mode samples
-            and a xr.DataArray of selected channels for BB mode samples
+        return prx
 
-        Returns
-        -------
-        An xr.DataArray
-        """
-        if waveform_mode == "BB":
-            gain_single = self._get_vend_cal_params_power(
-                "gain_correction", waveform_mode=waveform_mode
-            )
-            gain = []
-            if "gain" in self.echodata["Vendor_specific"].data_vars:
-                # index using channel_id as order of frequency across channel can be arbitrary
-                # reference to freq_center in case some channels are CW complex samples
-                # (already dropped when computing freq_center in the calling function)
-                for ch_id in chan_sel:
-                    # if channel gain exists in data
-                    if ch_id in self.echodata["Vendor_specific"].cal_channel_id:
-                        gain_vec = self.echodata["Vendor_specific"].gain.sel(cal_channel_id=ch_id)
-                        gain_temp = (
-                            gain_vec.interp(
-                                cal_frequency=self.echodata[
-                                    "Vendor_specific"
-                                ].frequency_nominal.sel(channel=ch_id)
-                            ).drop(["cal_channel_id", "cal_frequency"])
-                        ).expand_dims("channel")
-                    # if no freq-dependent gain use CW gain
-                    else:
-                        gain_temp = (
-                            gain_single.sel(channel=ch_id)
-                            .reindex_like(
-                                self.echodata["Sonar/Beam_group1"].backscatter_r, method="nearest"
-                            )
-                            .expand_dims("channel")
-                        )
-                    gain_temp.name = "gain"
-                    gain.append(gain_temp)
-                gain = xr.merge(gain).gain  # select the single data variable
-            else:
-                gain = gain_single
-        elif waveform_mode == "CW":
-            gain = self._get_vend_cal_params_power(
-                "gain_correction", waveform_mode=waveform_mode
-            ).sel(channel=chan_sel)
-
-        return gain
-
-    def _cal_complex(self, cal_type, waveform_mode) -> xr.Dataset:
+    def _cal_complex_samples(self, cal_type: str, complex_ed_group: str) -> xr.Dataset:
         """Calibrate complex data from EK80.
 
         Parameters
@@ -677,142 +310,83 @@ def _cal_complex(self, cal_type, waveform_mode) -> xr.Dataset:
         cal_type : str
             'Sv' for calculating volume backscattering strength, or
             'TS' for calculating target strength
-        waveform_mode : {"CW", "BB"}
-            Type of transmit waveform.
-
-            - `"CW"` for narrowband transmission,
-              returned echoes recorded either as complex or power/angle samples
-            - `"BB"` for broadband transmission,
-              returned echoes recorded as complex samples
+        complex_ed_group : str
+            The ``EchoData`` beam group path containing complex data
 
         Returns
         -------
         xr.Dataset
             The calibrated dataset containing Sv or TS
         """
-        # Transmit replica and effective pulse length
-        chirp, _, tau_effective = self.get_transmit_chirp(waveform_mode=waveform_mode)
-
-        # use center frequency for each ping to select BB or CW channels
-        # when all samples are encoded as complex samples
-        if (
-            "frequency_start" in self.echodata["Sonar/Beam_group1"]
-            and "frequency_end" in self.echodata["Sonar/Beam_group1"]
-        ):
-            freq_center = (
-                self.echodata["Sonar/Beam_group1"]["frequency_start"]
-                + self.echodata["Sonar/Beam_group1"]["frequency_end"]
-            ) / 2  # has beam dim
-        else:
-            freq_center = None
-
-        if waveform_mode == "BB":
-            if freq_center is None:
-                raise ValueError(
-                    "frequency_start and frequency_end should exist in BB mode data, "
-                    "double check the EchoData object!"
-                )
-            # if CW and BB complex samples co-exist
-            # drop those that contain CW samples (nan in freq start/end)
-            chan_sel = freq_center.dropna(dim="channel").channel
-
-            # backscatter data
-            pc = self.compress_pulse(chirp, chan_BB=chan_sel)  # has beam dim
-            prx = (
-                self.echodata["Sonar/Beam_group1"].beam.size
-                * np.abs(pc.mean(dim="beam")) ** 2
-                / (2 * np.sqrt(2)) ** 2
-                * (np.abs(self.z_er + self.z_et) / self.z_er) ** 2
-                / self.z_et
-            )
-        else:
-            if freq_center is None:
-                # when only have CW complex samples
-                chan_sel = self.echodata["Sonar/Beam_group1"].channel
-            else:
-                # if BB and CW complex samples co-exist
-                # drop those that contain BB samples (not nan in freq start/end)
-                chan_sel = freq_center.where(np.isnan(freq_center), drop=True).channel
-
-            # backscatter data
-            backscatter_cw = (
-                self.echodata["Sonar/Beam_group1"]["backscatter_r"]
-                + 1j * self.echodata["Sonar/Beam_group1"]["backscatter_i"]
-            )
-            prx = (
-                self.echodata["Sonar/Beam_group1"].beam.size
-                * np.abs(backscatter_cw.mean(dim="beam")) ** 2
-                / (2 * np.sqrt(2)) ** 2
-                * (np.abs(self.z_er + self.z_et) / self.z_er) ** 2
-                / self.z_et
-            )
-            prx.name = "received_power"
-            prx = prx.to_dataset()
+        # Select source of backscatter data
+        beam = self.echodata[complex_ed_group]
+
+        # Get transmit signal
+        tx_coeff = self._get_filter_coeff(channel=self.chan_sel)
+        tx, tx_time = get_transmit_signal(
+            beam=beam,
+            coeff=tx_coeff,
+            waveform_mode=self.waveform_mode,
+            channel=self.chan_sel,
+            fs=self.fs,
+            z_et=self.z_et,
+        )
 
-        # Compute derived params
+        # Get power from complex samples
+        prx = self._get_power_from_complex(beam=beam, chan_sel=self.chan_sel, chirp=tx)
 
         # Harmonize time coordinate between Beam_groupX data and env_params
         # Use self.echodata["Sonar/Beam_group1"] because complex sample is always in Beam_group1
         for p in self.env_params.keys():
             if "channel" in self.env_params[p].coords:
-                self.env_params[p] = self.env_params[p].sel(channel=chan_sel)
+                self.env_params[p] = self.env_params[p].sel(channel=self.chan_sel)
             self.env_params[p] = self.echodata._harmonize_env_param_time(
-                self.env_params[p], ping_time=self.echodata["Sonar/Beam_group1"].ping_time
+                self.env_params[p], ping_time=beam["ping_time"]
             )
 
+        # Common params
         sound_speed = self.env_params["sound_speed"]
-        absorption = self.env_params["sound_absorption"].sel(channel=chan_sel)
-        range_meter = self.range_meter.sel(channel=chan_sel)
-        if waveform_mode == "BB":
-            # use true center frequency for BB pulse
-            wavelength = sound_speed / self.echodata["Sonar/Beam_group1"].frequency_nominal.sel(
-                channel=chan_sel
+        absorption = self.env_params["sound_absorption"].sel(channel=self.chan_sel)
+        range_meter = self.range_meter.sel(channel=self.chan_sel)
+        wavelength = sound_speed / self.freq_center
+
+        if self.waveform_mode == "BB" and "gain" in self.echodata["Vendor_specific"]:
+            # If frequency-dependent gain exists, interpolate at true center frequency
+            gain = get_gain_BB(
+                vend=self.echodata["Vendor_specific"],
+                freq_center=self.freq_center,
+                cal_params_CW=self.cal_params,
             )
-
-            # use true center frequency to interpolate for gain factor
-            gain = self._get_gain_for_complex(waveform_mode=waveform_mode, chan_sel=chan_sel)
-
         else:
-            # use nominal channel frequency for CW pulse
-            wavelength = sound_speed / self.echodata["Sonar/Beam_group1"].frequency_nominal.sel(
-                channel=chan_sel
-            )
+            # use gain already retrieved in init
+            gain = self.cal_params["gain_correction"]
 
-            # use nominal channel frequency to select gain factor
-            gain = self._get_gain_for_complex(waveform_mode=waveform_mode, chan_sel=chan_sel)
-
-        # Transmission loss
         spreading_loss = 20 * np.log10(range_meter.where(range_meter >= 1, other=1))
         absorption_loss = 2 * absorption * range_meter
+        transmit_power = beam["transmit_power"].sel(channel=self.chan_sel)
 
-        # TODO: both Sv and TS are off by ~<0.5 dB from matlab outputs.
-        #  Is this due to the use of 'single' in matlab code?
+        # Compute based on cal_type
         if cal_type == "Sv":
-            # effective pulse length
-            tau_effective = xr.DataArray(
-                data=list(tau_effective.values()),
-                coords=[
-                    self.echodata["Sonar/Beam_group1"].channel,
-                    self.echodata["Sonar/Beam_group1"].ping_time,
-                ],
-                dims=["channel", "ping_time"],
-            ).sel(channel=chan_sel)
-
-            # other params
-            transmit_power = self.echodata["Sonar/Beam_group1"]["transmit_power"].sel(
-                channel=chan_sel
+            # Compute effective pulse length
+            tau_effective = get_tau_effective(
+                ytx_dict=tx,
+                fs_deci_dict={k: np.diff(v[:2]) for (k, v) in tx_time.items()},
+                waveform_mode=self.waveform_mode,
+                channel=self.chan_sel,
+                ping_time=beam["ping_time"],
             )
-            # equivalent_beam_angle has beam dim
-            if waveform_mode == "BB":
-                psifc = self.echodata["Sonar/Beam_group1"]["equivalent_beam_angle"].sel(
-                    channel=chan_sel
-                ) + 10 * np.log10(
-                    self.echodata["Vendor_specific"].frequency_nominal.sel(channel=chan_sel)
-                    / freq_center
-                )
-            elif waveform_mode == "CW":
-                psifc = self.echodata["Sonar/Beam_group1"]["equivalent_beam_angle"].sel(
-                    channel=chan_sel
+
+            # equivalent_beam_angle:
+            # identical within each channel regardless of ping_time/beam
+            # but drop only the beam dimension here
+            # to allow scaling for potential center frequency changes
+            psifc = (
+                beam["equivalent_beam_angle"].sel(channel=self.chan_sel).isel(beam=0).drop("beam")
+            )
+            if self.waveform_mode == "BB":
+                # if BB scale according to true center frequency
+                psifc += 20 * np.log10(  # TODO: BUGS! should be 20 * log10 [WJ resolved 2022/12/27]
+                    beam["frequency_nominal"].sel(channel=self.chan_sel) / self.freq_center
                 )
 
             out = (
@@ -824,27 +398,24 @@ def _cal_complex(self, cal_type, waveform_mode) -> xr.Dataset:
                 - 10 * np.log10(tau_effective)
                 - psifc
             )
-            out = out.rename_vars({list(out.data_vars.keys())[0]: "Sv"})
+            out.name = "Sv"
+            # out = out.rename_vars({list(out.data_vars.keys())[0]: "Sv"})
 
         elif cal_type == "TS":
-            transmit_power = self.echodata["Sonar/Beam_group1"]["transmit_power"].sel(
-                channel=chan_sel
-            )
-
             out = (
                 10 * np.log10(prx)
                 + 2 * spreading_loss
-                + absorption_loss  # has beam dim
+                + absorption_loss
                 - 10 * np.log10(wavelength**2 * transmit_power / (16 * np.pi**2))
                 - 2 * gain
             )
-            out = out.rename_vars({list(out.data_vars.keys())[0]: "TS"})
+            out.name = "TS"
 
         # Attach calculated range (with units meter) into data set
-        out = out.merge(range_meter)
+        out = out.to_dataset().merge(range_meter)
 
         # Add frequency_nominal to data set
-        out["frequency_nominal"] = self.echodata["Sonar/Beam_group1"]["frequency_nominal"]
+        out["frequency_nominal"] = beam["frequency_nominal"]
 
         # Add env and cal parameters
         out = self._add_params_to_output(out)
@@ -852,9 +423,12 @@ def _cal_complex(self, cal_type, waveform_mode) -> xr.Dataset:
         # Squeeze out the beam dim
         # out has beam dim, which came from absorption and absorption_loss
         # self.cal_params["equivalent_beam_angle"] also has beam dim
+
+        # TODO: out should not have beam dimension at this stage
+        # once that dimension is removed from equivalent_beam_angle
         return out.isel(beam=0).drop("beam")
 
-    def _compute_cal(self, cal_type, waveform_mode, encode_mode) -> xr.Dataset:
+    def _compute_cal(self, cal_type) -> xr.Dataset:
         """
         Private method to compute Sv or TS from EK80 data, called by compute_Sv or compute_TS.
 
@@ -864,111 +438,45 @@ def _compute_cal(self, cal_type, waveform_mode, encode_mode) -> xr.Dataset:
             'Sv' for calculating volume backscattering strength, or
             'TS' for calculating target strength
 
-        waveform_mode : {"CW", "BB"}
-            Type of transmit waveform.
-
-            - `"CW"` for narrowband transmission,
-              returned echoes recorded either as complex or power/angle samples
-            - `"BB"` for broadband transmission,
-              returned echoes recorded as complex samples
-
-        encode_mode : {"complex", "power"}
-            Type of encoded return echo data.
-
-            - `"complex"` for complex samples
-            - `"power"` for power/angle samples, only allowed when
-              the echosounder is configured for narrowband transmission
-
         Returns
         -------
         xr.Dataset
             An xarray Dataset containing either Sv or TS.
         """
-
-        power_ed_group = retrieve_correct_beam_group(
-            echodata=self.echodata,
-            waveform_mode=waveform_mode,
-            encode_mode=encode_mode,
-            pulse_compression=False,
+        # Set flag_complex: True-complex cal, False-power cal
+        flag_complex = (
+            True if self.waveform_mode == "BB" or self.encode_mode == "complex" else False
         )
 
-        # Set flag_complex
-        #  - True: complex cal
-        #  - False: power cal
-        flag_complex = False
-        if waveform_mode == "BB":
-            flag_complex = True
-        elif encode_mode == "complex":
-            flag_complex = True
-
-        # Compute Sv
         if flag_complex:
             # Complex samples can be BB or CW
-            self.compute_range_meter(waveform_mode=waveform_mode, encode_mode=encode_mode)
-            ds_cal = self._cal_complex(cal_type=cal_type, waveform_mode=waveform_mode)
+            self.compute_echo_range()
+            ds_cal = self._cal_complex_samples(cal_type=cal_type, complex_ed_group=self.ed_group)
         else:
             # Power samples only make sense for CW mode data
-            self.compute_range_meter(waveform_mode="CW", encode_mode=encode_mode)
-            ds_cal = self._cal_power(cal_type=cal_type, power_ed_group=power_ed_group)
+            self.compute_echo_range()
+            ds_cal = self._cal_power_samples(cal_type=cal_type, power_ed_group=self.ed_group)
 
         return ds_cal
 
-    def compute_Sv(self, waveform_mode="BB", encode_mode="complex"):
+    def compute_Sv(self):
         """Compute volume backscattering strength (Sv).
 
-        Parameters
-        ----------
-        waveform_mode : {"CW", "BB"}
-            Type of transmit waveform.
-
-            - `"CW"` for narrowband transmission,
-              returned echoes recorded either as complex or power/angle samples
-            - (default) `"BB"` for broadband transmission,
-              returned echoes recorded as complex samples
-
-        encode_mode : {"complex", "power"}
-            Type of encoded return echo data.
-
-            - (default) `"complex"` for complex samples
-            - `"power"` for power/angle samples, only allowed when
-              the echosounder is configured for narrowband transmission
-
         Returns
         -------
         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", waveform_mode=waveform_mode, encode_mode=encode_mode
-        )
+        return self._compute_cal(cal_type="Sv")
 
-    def compute_TS(self, waveform_mode="BB", encode_mode="complex"):
+    def compute_TS(self):
         """Compute target strength (TS).
 
-        Parameters
-        ----------
-        waveform_mode : {"CW", "BB"}
-            Type of transmit waveform.
-
-            - `"CW"` for narrowband transmission,
-              returned echoes recorded either as complex or power/angle samples
-            - (default) `"BB"` for broadband transmission,
-              returned echoes recorded as complex samples
-
-        encode_mode : {"complex", "power"}
-            Type of encoded return echo data.
-
-            - (default) `"complex"` for complex samples
-            - `"power"` for power/angle samples, only allowed when
-              the echosounder is configured for narrowband transmission
-
         Returns
         -------
         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", waveform_mode=waveform_mode, encode_mode=encode_mode
-        )
+        return self._compute_cal(cal_type="TS")
diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py
new file mode 100644
index 000000000..57584c641
--- /dev/null
+++ b/echopype/calibrate/ek80_complex.py
@@ -0,0 +1,247 @@
+from typing import Dict
+
+import numpy as np
+import xarray as xr
+from scipy import signal
+
+
+def tapered_chirp(
+    fs,
+    z_et,
+    transmit_duration_nominal,
+    slope,
+    transmit_power,
+    frequency_nominal=None,
+    frequency_start=None,
+    frequency_end=None,
+):
+    """Create a baseline chirp template."""
+    if frequency_start is None and frequency_end is None:  # CW waveform
+        frequency_start = frequency_nominal
+        frequency_end = frequency_nominal
+
+    t = np.arange(0, transmit_duration_nominal, 1 / fs)
+    nwtx = int(2 * np.floor(slope * t.size))  # length of tapering window
+    wtx_tmp = np.hanning(nwtx)  # hanning window
+    nwtxh = int(np.round(nwtx / 2))  # half length of the hanning window
+    wtx = np.concatenate(
+        [wtx_tmp[0:nwtxh], np.ones((t.size - nwtx)), wtx_tmp[nwtxh:]]
+    )  # assemble full tapering window
+    chirp_fac = (
+        (frequency_end - frequency_start) / transmit_duration_nominal
+    ) * t / 2 + frequency_start
+    y_tmp = (
+        np.sqrt((transmit_power / 4) * (2 * z_et))  # amplitude
+        * np.cos(2 * np.pi * chirp_fac * t)  # chirp
+        * wtx  # tapering
+    )  # taper and scale linear chirp
+    return y_tmp / np.max(np.abs(y_tmp)), t  # amp has no actual effect
+
+
+def filter_decimate_chirp(coeff_ch: Dict, y_ch: np.array, fs: float):
+    """Filter and decimate the transmit replica for one channel.
+
+    Parameters
+    ----------
+    coeff_ch : dict
+        a dictionary containing filter coefficients and decimation factors for ``ch_id``
+    y_ch : np.array
+        chirp from _tapered_chirp
+    fs : float
+        system sampling frequency [Hz]
+    """
+    # Get values
+
+    # WBT filter and decimation
+    ytx_wbt = signal.convolve(y_ch, coeff_ch["wbt_fil"])
+    ytx_wbt_deci = ytx_wbt[0 :: coeff_ch["wbt_decifac"]]
+
+    # PC filter and decimation
+    if len(coeff_ch["pc_fil"].squeeze().shape) == 0:  # in case it is a single element
+        coeff_ch["pc_fil"] = [coeff_ch["pc_fil"].squeeze()]
+    ytx_pc = signal.convolve(ytx_wbt_deci, coeff_ch["pc_fil"])
+    ytx_pc_deci = ytx_pc[0 :: coeff_ch["pc_decifac"]]
+    ytx_pc_deci_time = (
+        np.arange(ytx_pc_deci.size) * 1 / fs * coeff_ch["wbt_decifac"] * coeff_ch["pc_decifac"]
+    )
+
+    return ytx_pc_deci, ytx_pc_deci_time
+
+
+def get_tau_effective(
+    ytx_dict: Dict[str, np.array],
+    fs_deci_dict: Dict[str, float],
+    waveform_mode: str,
+    channel: xr.DataArray,
+    ping_time: xr.DataArray,
+):
+    """Compute effective pulse length.
+
+    Parameters
+    ----------
+    ytx_dict : dict
+        A dict of transmit signals, with keys being the ``channel`` and
+        values being either a vector when the transmit signals are identical across all pings
+        or a 2D array when the transmit signals vary across ping
+    fs_deci_dict : dict
+        A dict of sampling frequency of the decimated (recorded) signal,
+        with keys being the ``channel``
+    waveform_mode : str
+        ``CW`` for CW-mode samples, either recorded as complex or power samples
+        ``BB`` for BB-mode samples, recorded as complex samples
+    """
+    tau_effective = {}
+    for ch, ytx in ytx_dict.items():
+        if waveform_mode == "BB":
+            ytxa = signal.convolve(ytx, np.flip(np.conj(ytx))) / np.linalg.norm(ytx) ** 2
+            ptxa = np.abs(ytxa) ** 2
+        elif waveform_mode == "CW":
+            ptxa = np.abs(ytx) ** 2  # energy of transmit signal
+        tau_effective[ch] = ptxa.sum() / (ptxa.max() * fs_deci_dict[ch])
+
+    # set up coordinates
+    if len(ytx.shape) == 1:  # ytx is a vector (transmit signals are identical across pings)
+        coords = {"channel": channel}
+    elif len(ytx.shape) == 2:  # ytx is a matrix (transmit signals vary across pings)
+        coords = {"channel": channel, "ping_time": ping_time}
+
+    vals = np.array(list(tau_effective.values())).squeeze()
+    if vals.size == 1:
+        vals = np.expand_dims(vals, axis=0)
+
+    tau_effective = xr.DataArray(
+        data=vals,
+        coords=coords,
+    )
+
+    return tau_effective
+
+
+def get_transmit_signal(
+    beam: xr.Dataset, coeff: Dict, waveform_mode: str, channel: xr.DataArray, fs: float, z_et: float
+):
+    """Reconstruct transmit signal and compute effective pulse length.
+
+    Parameters
+    ----------
+    beam : xr.Dataset
+        EchoData["Sonar/Beam_group1"]
+    coeff : dict
+        a dictionary indexed by ``channel`` and values being dictionaries containing
+        filter coefficients and decimation factors for constructing the transmit replica.
+    waveform_mode : str
+        ``CW`` for CW-mode samples, either recorded as complex or power samples
+        ``BB`` for BB-mode samples, recorded as complex samples
+    channel : list or xr.DataArray
+        channel names (channel id), either as a list or an xr.DataArray
+
+    Return
+    ------
+    y_all
+        Transmit replica (BB: broadband chirp, CW: constant frequency sinusoid)
+    y_time_all
+        Timestamp for the transmit replica
+    """
+    # Make sure it is BB mode data
+    # This is already checked in calibrate_ek
+    # but keeping this here for use as standalone function
+    if waveform_mode == "BB" and (("frequency_start" not in beam) or ("frequency_end" not in beam)):
+        raise TypeError("File does not contain BB mode complex samples!")
+
+    # Build channel list
+    if not isinstance(channel, (list, xr.DataArray)):
+        raise ValueError("channel must be a list or an xr.DataArray!")
+    else:
+        if isinstance(channel, xr.DataArray):
+            ch_list = channel.values
+
+    # Generate all transmit replica
+    y_all = {}
+    y_time_all = {}
+    for ch in ch_list:
+        # TODO: expand to deal with the case with varying tx param across ping_time
+        if waveform_mode == "BB":
+            tx_param_names = [
+                "transmit_duration_nominal",
+                "slope",
+                "transmit_power",
+                "frequency_start",
+                "frequency_end",
+            ]
+        else:
+            tx_param_names = [
+                "transmit_duration_nominal",
+                "slope",
+                "transmit_power",
+                "frequency_nominal",
+            ]
+        tx_params = {}
+        for p in tx_param_names:
+            tx_params[p] = np.unique(beam[p].sel(channel=ch))
+            if tx_params[p].size != 1:
+                raise TypeError("File contains changing %s!" % p)
+        tx_params["fs"] = fs
+        tx_params["z_et"] = z_et
+        y_ch, _ = tapered_chirp(**tx_params)
+
+        # Filter and decimate chirp template
+        y_ch, y_tmp_time = filter_decimate_chirp(coeff_ch=coeff[ch], y_ch=y_ch, fs=fs)
+
+        # Fill into output dict
+        y_all[ch] = y_ch
+        y_time_all[ch] = y_tmp_time
+
+    return y_all, y_time_all
+
+
+def compress_pulse(beam: xr.Dataset, chirp: Dict, chan_BB=None):
+    """Perform pulse compression on the backscatter data.
+
+    Parameters
+    ----------
+    beam : xr.Dataset
+        EchoData["Sonar/Beam_group1"]
+    chirp : dict
+        transmit chirp replica indexed by ``channel``
+    chan_BB : str
+        channels that transmit in BB mode
+        (since CW mode can be in mixed in complex samples too)
+    """
+    backscatter = beam["backscatter_r"].sel(channel=chan_BB) + 1j * beam["backscatter_i"].sel(
+        channel=chan_BB
+    )
+
+    pc_all = []
+    for chan in chan_BB:
+        backscatter_chan = (
+            backscatter.sel(channel=chan)
+            # .dropna(dim="range_sample", how="all")
+            .dropna(dim="beam", how="all")
+            # .dropna(dim="ping_time")
+        )
+
+        tx = chirp[str(chan.values)]
+        replica = np.flipud(np.conj(tx))
+        pc = xr.apply_ufunc(
+            lambda m: np.apply_along_axis(
+                lambda m: (
+                    signal.convolve(m, replica, mode="full")[tx.size - 1 :]
+                    / np.linalg.norm(tx) ** 2
+                ),
+                axis=2,
+                arr=m,
+            ),
+            backscatter_chan,
+            input_core_dims=[["range_sample"]],
+            output_core_dims=[["range_sample"]],
+            # exclude_dims={"range_sample"},
+        )
+
+        # Expand dimension and add name to allow merge
+        pc = pc.expand_dims(dim="channel")
+        pc.name = "pulse_compressed_output"
+        pc_all.append(pc)
+
+    pc_merge = xr.merge(pc_all)
+
+    return pc_merge
diff --git a/echopype/calibrate/env_params.py b/echopype/calibrate/env_params.py
index ba6f2d588..77b223600 100644
--- a/echopype/calibrate/env_params.py
+++ b/echopype/calibrate/env_params.py
@@ -1,203 +1,204 @@
-from typing import Dict, List
+from typing import Dict, Optional
 
 import numpy as np
-import scipy.interpolate
 import xarray as xr
-from typing_extensions import Literal
-
-DataKind = Literal["stationary", "mobile", "organized"]
-InterpMethod = Literal["linear", "nearest", "zero", "slinear", "quadratic", "cubic"]
-ExtrapMethod = Literal["linear", "nearest"]
-VALID_INTERP_METHODS: Dict[DataKind, List[InterpMethod]] = {
-    "stationary": ["linear", "nearest", "zero", "slinear", "quadratic", "cubic"],
-    "mobile": ["linear", "nearest", "cubic"],
-    "organized": ["linear", "nearest"],
-}
-
-
-class EnvParams:
-    def __init__(
-        self,
-        env_params: xr.Dataset,
-        data_kind: DataKind,
-        interp_method: InterpMethod = "linear",
-        extrap_method: ExtrapMethod = "linear",
-    ):
-        """
-        Class to hold and interpolate external environmental data for calibration purposes.
-
-        This class can be used as the `env_params` parameter in `echopype.calibrate.compute_Sv`
-        or `echopype.calibrate.compute_TS`. It is intended to be used with environmental parameters
-        indexed by time. Environmental parameters will be interpolated onto dimensions within
-        the Platform group of the `EchoData` object being used for calibration.
-
-        Parameters
-        ----------
-        env_params : xr.Dataset
-            The environmental parameters to use for calibration. This data will be interpolated with
-            a provided `EchoData` object.
-
-            When `data_kind` is `"stationary"`, env_params must have a coordinate `"time3"`.
-            When `data_kind` is `"mobile"`, env_params must have coordinates `"latitude"`
-            and `"longitude"`.
-            When `data_kind` is `"organized"`, env_params must have coordinates `"time"`,
-            `"latitude"`, and `"longitude"`. This `data_kind` is not currently supported.
-        data_kind : {"stationary", "mobile", "organized"}
-            The type of the environmental parameters.
-
-            `"stationary"`: environmental parameters from a fixed location
-            (for example, a single CTD).
-            `"mobile"` environmental parameters from a moving location (for example, a ship).
-            `"organized"`: environmental parameters from many fixed locations
-            (for example, multiple CTDs).
-        interp_method: {"linear", "nearest", "zero", "slinear", "quadratic", "cubic"}
-            Method for interpolation of environmental parameters with the data from the
-            provided `EchoData` object.
-
-            When `data_kind` is `"stationary"`, valid `interp_method`s are `"linear"`, `"nearest"`,
-            `"zero"`, `"slinear"`, `"quadratic"`, and `"cubic"`
-            (see <https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.interpolate.interp1d.html>).
-            When `data_kind` is `"mobile"`, valid `interp_method`s are `"linear"`, `"nearest"`, and `"cubic"`
-            (see <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html>).
-            When `data_kind` is `"organized"`, valid `interp_method`s are `"linear"` and `"nearest"`
-            (see <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interpn.html>).
-        extrap_method: {"linear", "nearest"}
-            Method for extrapolation of environmental parameters with the data from the
-            provided `EchoData` object. Currently only supported when `data_kind` is `"stationary"`.
-
-        Notes
-        -----
-        Currently cases where `data_kind` is `"organized"` are not supported; support will be added
-        in a future version.
-
-        Examples
-        --------
-        >>> env_params = xr.open_dataset("env_params.nc")
-        >>> EnvParams(env_params, data_kind="mobile", interp_method="linear")
-        >>> echopype.calibrate.compute_Sv(echodata, env_params=env_params)
-        """  # noqa
-        if interp_method not in VALID_INTERP_METHODS[data_kind]:
-            raise ValueError(f"invalid interp_method {interp_method} for data_kind {data_kind}")
-
-        self.env_params = env_params
-        self.data_kind = data_kind
-        self.interp_method = interp_method
-        self.extrap_method = extrap_method
-
-    def _apply(self, echodata) -> Dict[str, xr.DataArray]:
-        if self.data_kind == "stationary":
-            dims = ["time3"]
-        elif self.data_kind == "mobile":
-            dims = ["latitude", "longitude"]
-        elif self.data_kind == "organized":
-            dims = ["time", "latitude", "longitude"]
-        else:
-            raise ValueError("invalid data_kind")
-
-        for dim in dims:
-            if dim not in echodata["Platform"]:
-                raise ValueError(
-                    f"could not interpolate env_params; EchoData is missing dimension {dim}"
-                )
-
-        env_params = self.env_params
-
-        if self.data_kind == "mobile":
-            if np.isnan(echodata["Platform"]["time1"]).all():
-                raise ValueError("cannot perform mobile interpolation without time1")
-            # compute_range needs indexing by ping_time
-            interp_plat = echodata["Platform"].interp(
-                {"time1": echodata["Sonar/Beam_group1"]["ping_time"]}
-            )
 
-            result = {}
-            for var, values in env_params.data_vars.items():
-                points = np.column_stack(
-                    (env_params["latitude"].data, env_params["longitude"].data)
-                )
-                values = values.data
-                xi = np.column_stack(
-                    (
-                        interp_plat["latitude"].data,
-                        interp_plat["longitude"].data,
-                    )
-                )
-                interp = scipy.interpolate.griddata(points, values, xi, method=self.interp_method)
-                result[var] = ("time1", interp)
-            env_params = xr.Dataset(
-                # we expect env_params to have coordinate time1
-                data_vars=result,
-                coords={"time1": interp_plat["ping_time"].data},
-            )
-        else:
-            # TODO: organized case
-            min_max = {
-                dim: {"min": env_params[dim].min(), "max": env_params[dim].max()} for dim in dims
-            }
-
-            extrap = env_params.interp(
-                {dim: echodata["Platform"][dim].data for dim in dims},
-                method=self.extrap_method,
-                # scipy interp uses "extrapolate" but scipy interpn uses None
-                kwargs={"fill_value": "extrapolate" if len(dims) == 1 else None},
-            )
-            # only keep unique indexes; xarray requires that indexes be unique
-            extrap_unique_idx = {dim: np.unique(extrap[dim], return_index=True)[1] for dim in dims}
-            extrap = extrap.isel(**extrap_unique_idx)
-            interp = env_params.interp(
-                {dim: echodata["Platform"][dim].data for dim in dims},
-                method=self.interp_method,
+from ..echodata import EchoData
+from ..utils import uwa
+
+# TODO: create default dict with empty values but specific keys for out_dict
+# like in cal_params, right now it is initiated as {}
+
+
+def get_env_params_AZFP(echodata: EchoData, user_env_dict: Optional[dict] = None):
+    """Get env params using user inputs or values from data file.
+
+    Parameters
+    ----------
+    user_env_dict : dict
+
+    Returns
+    -------
+    A dictionary containing the calibration parameters.
+    """
+    out_dict = {}
+
+    # Temperature comes from either user input or data file
+    out_dict["temperature"] = (
+        user_env_dict["temperature"]
+        if "temperature" in user_env_dict
+        else echodata["Environment"]["temperature"]
+    )
+
+    # Salinity and pressure always come from user input
+    if ("salinity" not in user_env_dict) or ("pressure" not in user_env_dict):
+        raise ReferenceError("Please supply both salinity and pressure in env_params.")
+    else:
+        out_dict["salinity"] = user_env_dict["salinity"]
+        out_dict["pressure"] = user_env_dict["pressure"]
+
+    # Always calculate sound speed and absorption
+    out_dict["sound_speed"] = uwa.calc_sound_speed(
+        temperature=user_env_dict["temperature"],
+        salinity=user_env_dict["salinity"],
+        pressure=user_env_dict["pressure"],
+        formula_source="AZFP",
+    )
+    out_dict["sound_absorption"] = uwa.calc_absorption(
+        frequency=echodata["Sonar/Beam_group1"]["frequency_nominal"],
+        temperature=user_env_dict["temperature"],
+        salinity=user_env_dict["salinity"],
+        pressure=user_env_dict["pressure"],
+        formula_source="AZFP",
+    )
+
+    return out_dict
+
+
+def get_env_params_EK60(echodata: EchoData, user_env_dict: Optional[Dict] = None) -> Dict:
+    """Get env params using user inputs or values from data file.
+
+    EK60 file by default contains only sound speed and absorption.
+    In cases when temperature, salinity, and pressure values are supplied
+    by the user simultaneously, the sound speed and absorption are re-calculated.
+
+    Parameters
+    ----------
+    user_env_dict : dict
+
+    Returns
+    -------
+    A dictionary containing the calibration parameters.
+    """
+    # Initiate input/output dict
+    out_dict = {}
+    if user_env_dict is None:
+        user_env_dict = {}
+
+    # Re-calculate environment parameters if user supply all env variables
+    tsp_all_exist = np.all([p in user_env_dict for p in ["temperature", "salinity", "pressure"]])
+
+    if tsp_all_exist:
+        out_dict["sound_speed"] = uwa.calc_sound_speed(
+            temperature=user_env_dict["temperature"],
+            salinity=user_env_dict["salinity"],
+            pressure=user_env_dict["pressure"],
+        )
+        out_dict["sound_absorption"] = uwa.calc_absorption(
+            frequency=echodata["Sonar/Beam_group1"]["frequency_nominal"],
+            temperature=user_env_dict["temperature"],
+            salinity=user_env_dict["salinity"],
+            pressure=user_env_dict["pressure"],
+        )
+    # Otherwise get sound speed and absorption from user inputs or raw data file
+    else:
+        out_dict["sound_speed"] = (
+            user_env_dict["sound_speed"]
+            if "sound_speed" in user_env_dict
+            else echodata["Environment"]["sound_speed_indicative"]
+        )
+        out_dict["sound_absorption"] = (
+            user_env_dict["sound_absorption"]
+            if "sound_absorption" in user_env_dict
+            else echodata["Environment"]["absorption_indicative"]
+        )
+
+    return out_dict
+
+
+def get_env_params_EK80(
+    echodata: EchoData,
+    freq: xr.DataArray,
+    user_env_dict: Optional[Dict] = None,
+) -> Dict:
+    """Get env params using user inputs or values from data file.
+
+    EK80 file by default contains sound speed, temperature, depth, salinity, and acidity,
+    therefore absorption is always calculated unless it is supplied by the user.
+    In cases when temperature, salinity, and pressure values are supplied
+    by the user simultaneously, both the sound speed and absorption are re-calculated.
+
+    Parameters
+    ----------
+    user_env_dict : dict
+
+    waveform_mode : {"CW", "BB"}
+        Type of transmit waveform.
+
+        - `"CW"` for narrowband transmission,
+            returned echoes recorded either as complex or power/angle samples
+        - (default) `"BB"` for broadband transmission,
+            returned echoes recorded as complex samples
+
+    encode_mode : {"complex", "power"}
+        Type of encoded return echo data.
+
+        - (default) `"complex"` for complex samples
+        - `"power"` for power/angle samples, only allowed when
+            the echosounder is configured for narrowband transmission
+
+    Returns
+    -------
+    A dictionary containing the calibration parameters.
+    """
+
+    # Initiate input/output dict
+    out_dict = {}
+    if user_env_dict is None:
+        user_env_dict = {}
+
+    # Re-calculate environment parameters if user supply all env variables
+    tsp_all_exist = np.all([p in user_env_dict for p in ["temperature", "salinity", "pressure"]])
+
+    if tsp_all_exist:
+        out_dict["sound_speed"] = uwa.calc_sound_speed(
+            temperature=user_env_dict["temperature"],
+            salinity=user_env_dict["salinity"],
+            pressure=user_env_dict["pressure"],
+        )
+        out_dict["sound_absorption"] = uwa.calc_absorption(
+            frequency=freq,
+            temperature=user_env_dict["temperature"],
+            salinity=user_env_dict["salinity"],
+            pressure=user_env_dict["pressure"],
+        )
+    # Otherwise
+    #  get temperature, salinity, and pressure from raw data file
+    #  get sound speed from user inputs or raw data file
+    #  get absorption from user inputs or computing from env params stored in raw data file
+    else:
+        # pressure is encoded as "depth" in EK80  # TODO: change depth to pressure in EK80 file?
+        for p1, p2 in zip(
+            ["temperature", "salinity", "pressure"],
+            ["temperature", "salinity", "depth"],
+        ):
+            out_dict[p1] = user_env_dict[p1] if p1 in user_env_dict else echodata["Environment"][p2]
+        out_dict["sound_speed"] = (
+            user_env_dict["sound_speed"]
+            if "sound_speed" in user_env_dict
+            else echodata["Environment"]["sound_speed_indicative"]
+        )
+        out_dict["sound_absorption"] = (
+            user_env_dict["sound_absorption"]
+            if "sound_absorption" in user_env_dict
+            else uwa.calc_absorption(
+                frequency=freq,
+                temperature=out_dict["temperature"],
+                salinity=out_dict["salinity"],
+                pressure=out_dict["pressure"],
             )
-            interp_unique_idx = {dim: np.unique(interp[dim], return_index=True)[1] for dim in dims}
-            interp = interp.isel(**interp_unique_idx)
-
-            if self.extrap_method is not None:
-                less = extrap.sel(
-                    {dim: extrap[dim][extrap[dim] < min_max[dim]["min"]] for dim in dims}
-                )
-                middle = interp.sel(
-                    {
-                        dim: interp[dim][
-                            np.logical_and(
-                                interp[dim] >= min_max[dim]["min"],
-                                interp[dim] <= min_max[dim]["max"],
-                            )
-                        ]
-                        for dim in dims
-                    }
-                )
-                greater = extrap.sel(
-                    {dim: extrap[dim][extrap[dim] > min_max[dim]["max"]] for dim in dims}
-                )
-
-                # remove empty datasets (xarray does not allow any dims from any datasets
-                # to be length 0 in combine_by_coords)
-                non_zero_dims = [
-                    ds
-                    for ds in (less, middle, greater)
-                    if all(dim_len > 0 for dim_len in ds.dims.values())
-                ]
-                env_params = xr.combine_by_coords(non_zero_dims)
-
-        # if self.data_kind == "organized":
-        #     # get platform latitude and longitude indexed by ping_time
-        #     interp_plat = echodata["Platform"].interp(
-        #         {"time": echodata["Platform"]["ping_time"]}
-        #     )
-        #     # get env_params latitude and longitude indexed by ping_time
-        #     env_params = env_params.interp(
-        #         {
-        #             "latitude": interp_plat["latitude"],
-        #             "longitude": interp_plat["longitude"],
-        #         }
-        #     )
-
-        if self.data_kind == "stationary":
-            # renaming time3 (from Platform group) to time1 because we expect
-            # environmental parameters to have dimension time1
-            return {
-                var: env_params[var].rename({"time3": "time1"})
-                for var in ("temperature", "salinity", "pressure")
-            }
-        else:
-            return {var: env_params[var] for var in ("temperature", "salinity", "pressure")}
+        )
+
+    return out_dict
+
+
+# TODO: this function is currently unused, consider removing
+def get_env_params(
+    sonar_model: str, echodata: EchoData, env_params: Optional[Dict] = None, **kwarg
+):
+    if sonar_model == "AZFP":
+        return get_env_params_AZFP(echodata=echodata, user_env_dict=env_params)
+    elif sonar_model in ["EK60", "ES70"]:
+        return get_env_params_EK60(echodata=echodata, user_env_dict=env_params)
+    elif sonar_model in ["EK80", "ES80", "EA640"]:
+        return get_env_params_EK80(echodata=echodata, user_env_dict=env_params, **kwarg)
diff --git a/echopype/calibrate/env_params_old.py b/echopype/calibrate/env_params_old.py
new file mode 100644
index 000000000..ba6f2d588
--- /dev/null
+++ b/echopype/calibrate/env_params_old.py
@@ -0,0 +1,203 @@
+from typing import Dict, List
+
+import numpy as np
+import scipy.interpolate
+import xarray as xr
+from typing_extensions import Literal
+
+DataKind = Literal["stationary", "mobile", "organized"]
+InterpMethod = Literal["linear", "nearest", "zero", "slinear", "quadratic", "cubic"]
+ExtrapMethod = Literal["linear", "nearest"]
+VALID_INTERP_METHODS: Dict[DataKind, List[InterpMethod]] = {
+    "stationary": ["linear", "nearest", "zero", "slinear", "quadratic", "cubic"],
+    "mobile": ["linear", "nearest", "cubic"],
+    "organized": ["linear", "nearest"],
+}
+
+
+class EnvParams:
+    def __init__(
+        self,
+        env_params: xr.Dataset,
+        data_kind: DataKind,
+        interp_method: InterpMethod = "linear",
+        extrap_method: ExtrapMethod = "linear",
+    ):
+        """
+        Class to hold and interpolate external environmental data for calibration purposes.
+
+        This class can be used as the `env_params` parameter in `echopype.calibrate.compute_Sv`
+        or `echopype.calibrate.compute_TS`. It is intended to be used with environmental parameters
+        indexed by time. Environmental parameters will be interpolated onto dimensions within
+        the Platform group of the `EchoData` object being used for calibration.
+
+        Parameters
+        ----------
+        env_params : xr.Dataset
+            The environmental parameters to use for calibration. This data will be interpolated with
+            a provided `EchoData` object.
+
+            When `data_kind` is `"stationary"`, env_params must have a coordinate `"time3"`.
+            When `data_kind` is `"mobile"`, env_params must have coordinates `"latitude"`
+            and `"longitude"`.
+            When `data_kind` is `"organized"`, env_params must have coordinates `"time"`,
+            `"latitude"`, and `"longitude"`. This `data_kind` is not currently supported.
+        data_kind : {"stationary", "mobile", "organized"}
+            The type of the environmental parameters.
+
+            `"stationary"`: environmental parameters from a fixed location
+            (for example, a single CTD).
+            `"mobile"` environmental parameters from a moving location (for example, a ship).
+            `"organized"`: environmental parameters from many fixed locations
+            (for example, multiple CTDs).
+        interp_method: {"linear", "nearest", "zero", "slinear", "quadratic", "cubic"}
+            Method for interpolation of environmental parameters with the data from the
+            provided `EchoData` object.
+
+            When `data_kind` is `"stationary"`, valid `interp_method`s are `"linear"`, `"nearest"`,
+            `"zero"`, `"slinear"`, `"quadratic"`, and `"cubic"`
+            (see <https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.interpolate.interp1d.html>).
+            When `data_kind` is `"mobile"`, valid `interp_method`s are `"linear"`, `"nearest"`, and `"cubic"`
+            (see <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html>).
+            When `data_kind` is `"organized"`, valid `interp_method`s are `"linear"` and `"nearest"`
+            (see <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interpn.html>).
+        extrap_method: {"linear", "nearest"}
+            Method for extrapolation of environmental parameters with the data from the
+            provided `EchoData` object. Currently only supported when `data_kind` is `"stationary"`.
+
+        Notes
+        -----
+        Currently cases where `data_kind` is `"organized"` are not supported; support will be added
+        in a future version.
+
+        Examples
+        --------
+        >>> env_params = xr.open_dataset("env_params.nc")
+        >>> EnvParams(env_params, data_kind="mobile", interp_method="linear")
+        >>> echopype.calibrate.compute_Sv(echodata, env_params=env_params)
+        """  # noqa
+        if interp_method not in VALID_INTERP_METHODS[data_kind]:
+            raise ValueError(f"invalid interp_method {interp_method} for data_kind {data_kind}")
+
+        self.env_params = env_params
+        self.data_kind = data_kind
+        self.interp_method = interp_method
+        self.extrap_method = extrap_method
+
+    def _apply(self, echodata) -> Dict[str, xr.DataArray]:
+        if self.data_kind == "stationary":
+            dims = ["time3"]
+        elif self.data_kind == "mobile":
+            dims = ["latitude", "longitude"]
+        elif self.data_kind == "organized":
+            dims = ["time", "latitude", "longitude"]
+        else:
+            raise ValueError("invalid data_kind")
+
+        for dim in dims:
+            if dim not in echodata["Platform"]:
+                raise ValueError(
+                    f"could not interpolate env_params; EchoData is missing dimension {dim}"
+                )
+
+        env_params = self.env_params
+
+        if self.data_kind == "mobile":
+            if np.isnan(echodata["Platform"]["time1"]).all():
+                raise ValueError("cannot perform mobile interpolation without time1")
+            # compute_range needs indexing by ping_time
+            interp_plat = echodata["Platform"].interp(
+                {"time1": echodata["Sonar/Beam_group1"]["ping_time"]}
+            )
+
+            result = {}
+            for var, values in env_params.data_vars.items():
+                points = np.column_stack(
+                    (env_params["latitude"].data, env_params["longitude"].data)
+                )
+                values = values.data
+                xi = np.column_stack(
+                    (
+                        interp_plat["latitude"].data,
+                        interp_plat["longitude"].data,
+                    )
+                )
+                interp = scipy.interpolate.griddata(points, values, xi, method=self.interp_method)
+                result[var] = ("time1", interp)
+            env_params = xr.Dataset(
+                # we expect env_params to have coordinate time1
+                data_vars=result,
+                coords={"time1": interp_plat["ping_time"].data},
+            )
+        else:
+            # TODO: organized case
+            min_max = {
+                dim: {"min": env_params[dim].min(), "max": env_params[dim].max()} for dim in dims
+            }
+
+            extrap = env_params.interp(
+                {dim: echodata["Platform"][dim].data for dim in dims},
+                method=self.extrap_method,
+                # scipy interp uses "extrapolate" but scipy interpn uses None
+                kwargs={"fill_value": "extrapolate" if len(dims) == 1 else None},
+            )
+            # only keep unique indexes; xarray requires that indexes be unique
+            extrap_unique_idx = {dim: np.unique(extrap[dim], return_index=True)[1] for dim in dims}
+            extrap = extrap.isel(**extrap_unique_idx)
+            interp = env_params.interp(
+                {dim: echodata["Platform"][dim].data for dim in dims},
+                method=self.interp_method,
+            )
+            interp_unique_idx = {dim: np.unique(interp[dim], return_index=True)[1] for dim in dims}
+            interp = interp.isel(**interp_unique_idx)
+
+            if self.extrap_method is not None:
+                less = extrap.sel(
+                    {dim: extrap[dim][extrap[dim] < min_max[dim]["min"]] for dim in dims}
+                )
+                middle = interp.sel(
+                    {
+                        dim: interp[dim][
+                            np.logical_and(
+                                interp[dim] >= min_max[dim]["min"],
+                                interp[dim] <= min_max[dim]["max"],
+                            )
+                        ]
+                        for dim in dims
+                    }
+                )
+                greater = extrap.sel(
+                    {dim: extrap[dim][extrap[dim] > min_max[dim]["max"]] for dim in dims}
+                )
+
+                # remove empty datasets (xarray does not allow any dims from any datasets
+                # to be length 0 in combine_by_coords)
+                non_zero_dims = [
+                    ds
+                    for ds in (less, middle, greater)
+                    if all(dim_len > 0 for dim_len in ds.dims.values())
+                ]
+                env_params = xr.combine_by_coords(non_zero_dims)
+
+        # if self.data_kind == "organized":
+        #     # get platform latitude and longitude indexed by ping_time
+        #     interp_plat = echodata["Platform"].interp(
+        #         {"time": echodata["Platform"]["ping_time"]}
+        #     )
+        #     # get env_params latitude and longitude indexed by ping_time
+        #     env_params = env_params.interp(
+        #         {
+        #             "latitude": interp_plat["latitude"],
+        #             "longitude": interp_plat["longitude"],
+        #         }
+        #     )
+
+        if self.data_kind == "stationary":
+            # renaming time3 (from Platform group) to time1 because we expect
+            # environmental parameters to have dimension time1
+            return {
+                var: env_params[var].rename({"time3": "time1"})
+                for var in ("temperature", "salinity", "pressure")
+            }
+        else:
+            return {var: env_params[var] for var in ("temperature", "salinity", "pressure")}
diff --git a/echopype/calibrate/range.py b/echopype/calibrate/range.py
new file mode 100644
index 000000000..e18d48f1d
--- /dev/null
+++ b/echopype/calibrate/range.py
@@ -0,0 +1,270 @@
+import datetime
+from typing import Dict, Optional, Union
+
+import numpy as np
+import xarray as xr
+
+from ..echodata import EchoData
+from ..echodata.simrad import retrieve_correct_beam_group
+
+TVG_CORRECTION_FACTOR = {
+    "EK60": 2,
+    "ES70": 2,
+    "EK80": 0,
+    "ES80": 0,
+    "EA640": 0,
+}
+
+
+def _harmonize_env_param_time(
+    p: Union[int, float, xr.DataArray],
+    ping_time: Optional[Union[xr.DataArray, datetime.datetime]] = None,
+):
+    """
+    Harmonize time coordinate between Beam_groupX data and env_params to make sure
+    the timestamps are broacast correctly in calibration and range calculations.
+
+    Regardless of the source, if `p` is an xr.DataArray, the time coordinate name
+    needs to be `time1` to be consistent with the time coordinate in EchoData["Environment"].
+    If `time1` is of length=1, the dimension `time1` is dropped.
+    Otherwise, `p` is interpolated to `ping_time`.
+    If `p` is not an xr.DataArray it is returned directly.
+
+    Parameters
+    ----------
+    p
+        The environment parameter for timestamp check/correction
+    ping_time
+        Beam_groupX ping_time to interpolate env_params timestamps to.
+        Only used if p.time1 has length >1
+
+    Returns
+    -------
+    Environment parameter with correctly broadcasted timestamps
+    """
+    if isinstance(p, xr.DataArray):
+        if "time1" not in p.coords:
+            return p
+        else:
+            # If there's only 1 time1 value,
+            # or if after dropping NaN there's only 1 time1 value
+            if p["time1"].size == 1 or p.dropna(dim="time1").size == 1:
+                return p.dropna(dim="time1").squeeze(dim="time1").drop("time1")
+
+            # Direct assignment if all timestamps are identical (EK60 data)
+            elif np.all(p["time1"].values == ping_time.values):
+                return p.rename({"time1": "ping_time"})
+
+            elif ping_time is None:
+                raise ValueError(f"ping_time needs to be provided for interpolating {p.name}")
+
+            else:
+                return p.dropna(dim="time1").interp(time1=ping_time)
+    else:
+        return p
+
+
+def compute_range_AZFP(echodata: EchoData, env_params: Dict, cal_type: str) -> xr.DataArray:
+    """
+    Computes the range (``echo_range``) of AZFP backscatter data in meters.
+
+    Parameters
+    ----------
+    echodata : EchoData
+        An EchoData object holding data from an AZFP echosounder
+    env_params : dict
+        A dictionary holding environmental parameters needed for computing range
+        See echopype.calibrate.env_params.get_env_params_AZFP()
+    cal_type : {"Sv", "TS"}
+
+        - `"Sv"` for calculating volume backscattering strength
+        - `"TS"` for calculating target strength.
+
+        This parameter needs to be specified for data from the AZFP echosounder
+        due to a difference in the range computation given by the manufacturer
+
+    Returns
+    -------
+    xr.DataArray
+        The range (``echo_range``) of the data in meters.
+
+    Notes
+    -----
+    For AZFP echosounder, the returned ``echo_range`` is duplicated along ``ping_time``
+    to conform with outputs from other echosounders, even though within each data
+    file the range is held constant.
+    """
+    # sound_speed should exist already
+    if "sound_speed" not in env_params:
+        raise RuntimeError(
+            "sounds_speed not included in env_params, "
+            "use echopype.calibrate.env_params.get_env_params_AZFP() to compute env_params "
+            "by supplying temperature, salinity, and pressure."
+        )
+    else:
+        sound_speed = env_params["sound_speed"]
+
+    # Check cal_type
+    if cal_type is None:
+        raise ValueError('cal_type must be "Sv" or "TS"')
+
+    # Groups to use
+    vend = echodata["Vendor_specific"]
+    beam = echodata["Sonar/Beam_group1"]
+
+    # Notation below follows p.86 of user manual
+    N = vend["number_of_samples_per_average_bin"]  # samples per bin
+    f = vend["digitization_rate"]  # digitization rate
+    L = vend["lockout_index"]  # number of lockout samples
+
+    # keep this in ref of AZFP matlab code,
+    # set to 1 since we want to calculate from raw data
+    bins_to_avg = 1
+
+    # Harmonize sound_speed time1 and Beam_group1 ping_time
+    sound_speed = _harmonize_env_param_time(
+        p=sound_speed,
+        ping_time=beam.ping_time,
+    )
+
+    # Calculate range using parameters for each freq
+    # This is "the range to the centre of the sampling volume
+    # for bin m" from p.86 of user manual
+    if cal_type == "Sv":
+        range_offset = 0
+    else:
+        range_offset = sound_speed * beam["transmit_duration_nominal"] / 4  # from matlab code
+    range_meter = (
+        sound_speed * L / (2 * f)
+        + (sound_speed / 4)
+        * (
+            ((2 * (beam["range_sample"] + 1) - 1) * N * bins_to_avg - 1) / f
+            + beam["transmit_duration_nominal"]
+        )
+        - range_offset
+    )
+
+    range_meter.name = "echo_range"  # add name to facilitate xr.merge
+
+    return range_meter
+
+
+def compute_range_EK(
+    echodata: EchoData, env_params: Dict, waveform_mode: str = "CW", encode_mode: str = "power"
+):
+    """
+    Computes the range (``echo_range``) of EK backscatter data in meters.
+
+    Parameters
+    ----------
+    echodata : EchoData
+        An EchoData object holding data from an AZFP echosounder
+    env_params : dict
+        A dictionary holding environmental parameters needed for computing range
+        See echopype.calibrate.env_params.get_env_params_EK60() or .get_env_params_EK80()
+    waveform_mode : {"CW", "BB"}
+        Type of transmit waveform.
+        Required only for data from the EK80 echosounder.
+
+        - `"CW"` for narrowband transmission,
+            returned echoes recorded either as complex or power/angle samples
+        - `"BB"` for broadband transmission,
+            returned echoes recorded as complex samples
+
+    encode_mode : {"complex", "power"}
+        Type of encoded data format.
+        Required only for data from the EK80 echosounder.
+
+        - `"complex"` for complex samples
+        - `"power"` for power/angle samples, only allowed when
+          the echosounder is configured for narrowband transmission
+
+    Returns
+    -------
+    xr.DataArray
+        The range (``echo_range``) of the data in meters.
+
+    Notes
+    -----
+    The EK80 echosounder can be configured to transmit
+    either broadband (``waveform_mode="BB"``) or narrowband (``waveform_mode="CW"``) signals.
+    When transmitting in broadband mode, the returned echoes must be
+    encoded as complex samples (``encode_mode="complex"``).
+    When transmitting in narrowband mode, the returned echoes can be encoded
+    either as complex samples (``encode_mode="complex"``)
+    or as power/angle combinations (``encode_mode="power"``) in a format
+    similar to those recorded by EK60 echosounders (the "power/angle" format).
+    """
+    # sound_speed should exist already
+    if echodata.sonar_model in ("EK60", "ES70"):
+        ek_str = "EK60"
+    elif echodata.sonar_model in ("EK80", "ES80", "EA640"):
+        ek_str = "EK80"
+    else:
+        raise ValueError("The specified sonar_model is not supported!")
+
+    if "sound_speed" not in env_params:
+        raise RuntimeError(
+            "sounds_speed not included in env_params, "
+            f"use echopype.calibrate.env_params.get_env_params_{ek_str}() to compute env_params "
+        )
+    else:
+        sound_speed = env_params["sound_speed"]
+
+    # Get the right Sonar/Beam_groupX group according to encode_mode
+    ed_group = retrieve_correct_beam_group(echodata, waveform_mode, encode_mode)
+    beam = echodata[ed_group]
+
+    # TVG correction factor changes depending when the echo recording starts
+    # wrt when the transmit signal is sent out.
+    # This implementation is different for EK60 and EK80.
+    tvg_correction_factor = TVG_CORRECTION_FACTOR[echodata.sonar_model]
+
+    if waveform_mode == "BB":
+
+        shift = beam["transmit_duration_nominal"]  # based on Lar Anderson's Matlab code
+
+        # Harmonize sound_speed time1 and Beam_group1 ping_time
+        sound_speed = _harmonize_env_param_time(
+            p=sound_speed,
+            ping_time=beam.ping_time,
+        )
+
+        # TODO: once we allow putting in arbitrary sound_speed,
+        # change below to use linearly-interpolated values
+        range_meter = (beam["range_sample"] * beam["sample_interval"] - shift) * sound_speed / 2
+
+    else:  # CW
+        # Harmonize sound_speed time1 and Beam_groupX ping_time
+        sound_speed = _harmonize_env_param_time(
+            p=sound_speed,
+            ping_time=beam.ping_time,
+        )
+
+        sample_thickness = beam["sample_interval"] * sound_speed / 2
+
+        # TODO: Check with the AFSC about the half sample difference
+        range_meter = (
+            beam.range_sample - tvg_correction_factor
+        ) * sample_thickness  # [frequency x range_sample]
+
+    # make order of dims conform with the order of backscatter data
+    range_meter = range_meter.transpose("channel", "ping_time", "range_sample")
+    range_meter = range_meter.where(range_meter > 0, 0)  # set negative ranges to 0
+
+    # set entries with NaN backscatter data to NaN
+    if "beam" in beam["backscatter_r"].dims:
+        # Drop beam because echo_range does not have beam dimension
+        valid_idx = ~beam["backscatter_r"].isel(beam=0).drop("beam").isnull()
+    else:
+        valid_idx = ~beam["backscatter_r"].isnull()
+    range_meter = range_meter.where(valid_idx)
+
+    # remove time1 if exists as a coordinate
+    if "time1" in range_meter.coords:
+        range_meter = range_meter.drop("time1")
+
+    # add name to facilitate xr.merge
+    range_meter.name = "echo_range"
+
+    return range_meter
diff --git a/echopype/echodata/echodata.py b/echopype/echodata/echodata.py
index cfa68e492..f98e39af8 100644
--- a/echopype/echodata/echodata.py
+++ b/echopype/echodata/echodata.py
@@ -14,7 +14,7 @@
 if TYPE_CHECKING:
     from ..core import EngineHint, FileFormatHint, PathHint, SonarModelsHint
 
-from ..calibrate.env_params import EnvParams
+from ..calibrate.env_params_old import EnvParams
 from ..utils.coding import set_time_encodings
 from ..utils.io import check_file_existence, sanitize_file_path
 from ..utils.log import _init_logger
@@ -519,7 +519,6 @@ def compute_range(
                 ) * sample_thickness  # [frequency x range_sample]
             elif waveform_mode == "BB":
                 beam = self["Sonar/Beam_group1"]  # always use the Beam group
-                # TODO: bug: right now only first ping_time has non-nan range
                 shift = beam["transmit_duration_nominal"]  # based on Lar Anderson's Matlab code
 
                 # Harmonize sound_speed time1 and Beam_group1 ping_time
@@ -530,6 +529,8 @@ def compute_range(
 
                 # TODO: once we allow putting in arbitrary sound_speed,
                 # change below to use linearly-interpolated values
+                # TODO: change to not shift range since now we remove
+                # the equivalent length of samples from pulse compression output
                 range_meter = (
                     (beam.range_sample * beam["sample_interval"] - shift) * sound_speed / 2
                 )
diff --git a/echopype/echodata/simrad.py b/echopype/echodata/simrad.py
index 195f68e1f..f5f9facf1 100644
--- a/echopype/echodata/simrad.py
+++ b/echopype/echodata/simrad.py
@@ -6,8 +6,8 @@
 from .echodata import EchoData
 
 
-def _check_input_args_combination(
-    waveform_mode: str, encode_mode: str, pulse_compression: bool
+def check_input_args_combination(
+    waveform_mode: str, encode_mode: str, pulse_compression: bool = None
 ) -> None:
     """
     Checks that the ``waveform_mode`` and ``encode_mode`` have
@@ -35,11 +35,12 @@ def _check_input_args_combination(
         raise ValueError("encode_mode='power' not allowed when waveform_mode='BB'!")
 
     # make sure that we have BB and complex inputs, if pulse compression is selected
-    if pulse_compression and ((waveform_mode != "BB") or (encode_mode != "complex")):
-        raise ValueError(
-            "Pulse compression can only be used with "
-            "waveform_mode='BB' and encode_mode='complex'"
-        )
+    if pulse_compression is not None:
+        if pulse_compression and ((waveform_mode != "BB") or (encode_mode != "complex")):
+            raise RuntimeError(
+                "Pulse compression can only be used with "
+                "waveform_mode='BB' and encode_mode='complex'"
+            )
 
 
 def _retrieve_correct_beam_group_EK60(
@@ -174,9 +175,7 @@ def _retrieve_correct_beam_group_EK80(
     return power_ed_group, complex_ed_group
 
 
-def retrieve_correct_beam_group(
-    echodata: EchoData, waveform_mode: str, encode_mode: str, pulse_compression: bool
-) -> str:
+def retrieve_correct_beam_group(echodata: EchoData, waveform_mode: str, encode_mode: str) -> str:
     """
     A function to make sure that the user has provided the correct
     ``waveform_mode`` and ``encode_mode`` inputs based off of the
@@ -201,9 +200,6 @@ def retrieve_correct_beam_group(
         The ``EchoData`` beam group path corresponding to the ``encode_mode`` input
     """
 
-    # checks input and logic of modes without referencing data
-    _check_input_args_combination(waveform_mode, encode_mode, pulse_compression)
-
     if echodata.sonar_model in ["EK60", "ES70"]:
 
         # initialize complex_data_location (needed only for EK60)


Test patch:
diff --git a/echopype/tests/calibrate/test_calibrate.py b/echopype/tests/calibrate/test_calibrate.py
index 1bbf4d1fa..df49c3738 100644
--- a/echopype/tests/calibrate/test_calibrate.py
+++ b/echopype/tests/calibrate/test_calibrate.py
@@ -4,7 +4,8 @@
 from scipy.io import loadmat
 import echopype as ep
 from echopype.calibrate.calibrate_ek import CalibrateEK80
-from echopype.calibrate.env_params import EnvParams
+from echopype.calibrate.env_params_old import EnvParams
+from echopype.calibrate.ek80_complex import get_transmit_signal, compress_pulse
 import xarray as xr
 
 
@@ -114,6 +115,9 @@ def check_output(da_cmp, cal_type):
 
 
 def test_compute_Sv_ek60_duplicated_freq(ek60_path):
+
+    # TODO: add comparison of actual values in this test
+
     ek60_raw_path = str(
         ek60_path.joinpath('DY1002_EK60-D20100318-T023008_rep_freq.raw')
     )
@@ -231,16 +235,16 @@ def test_compute_Sv_ek80_pc_echoview(ek80_path):
         waveform_mode="BB",
         encode_mode="complex",
     )
-    cal_obj.compute_range_meter(
-        waveform_mode="BB", encode_mode="complex"
-    )  # compute range [m]
-    chirp, _, tau_effective = cal_obj.get_transmit_chirp(waveform_mode="BB")
-    freq_center = (
-        echodata["Sonar/Beam_group1"]["frequency_start"] + echodata["Sonar/Beam_group1"]["frequency_end"]
-    ).dropna(
-        dim="channel"
-    ) / 2  # drop those that contain CW samples (nan in freq start/end)
-    pc = cal_obj.compress_pulse(chirp, chan_BB=freq_center.channel)
+    cal_obj.compute_echo_range()  # compute range [m]
+    beam = echodata["Sonar/Beam_group1"]
+    chan_sel = beam["channel"]  # only BB data exist
+
+    coeff = cal_obj._get_filter_coeff(channel=chan_sel)
+    chirp, _ = get_transmit_signal(
+        beam=beam, coeff=coeff, channel=chan_sel, waveform_mode="BB",
+        fs=cal_obj.fs, z_et=cal_obj.z_et)
+
+    pc = compress_pulse(beam=beam, chirp=chirp, chan_BB=chan_sel)
     pc_mean = (
         pc.pulse_compressed_output.isel(channel=1)
         .mean(dim='beam')
diff --git a/echopype/tests/consolidate/test_consolidate.py b/echopype/tests/consolidate/test_consolidate.py
index 41070a574..6501b2354 100644
--- a/echopype/tests/consolidate/test_consolidate.py
+++ b/echopype/tests/consolidate/test_consolidate.py
@@ -12,12 +12,12 @@
 import os
 
 """
-For future reference: 
+For future reference:
 
 For ``test_add_splitbeam_angle`` the test data is in the following locations:
-- the EK60 raw file is in `test_data/ek60/DY1801_EK60-D20180211-T164025.raw` and the 
+- the EK60 raw file is in `test_data/ek60/DY1801_EK60-D20180211-T164025.raw` and the
 associated echoview split-beam data is in `test_data/ek60/splitbeam`.
-- the EK80 raw file is in `test_data/ek80_bb_with_calibration/2018115-D20181213-T094600.raw` and 
+- the EK80 raw file is in `test_data/ek80_bb_with_calibration/2018115-D20181213-T094600.raw` and
 the associated echoview split-beam data is in `test_data/ek80_bb_with_calibration/splitbeam`
 """
 
@@ -127,7 +127,7 @@ def test_swap_dims_channel_frequency(test_data_samples):
 def _build_ds_Sv(channel, range_sample, ping_time, sample_interval):
     return xr.Dataset(
         data_vars={
-            "Sv": ( 
+            "Sv": (
                 ("channel", "range_sample", "ping_time"),
                 np.random.random((len(channel), range_sample.size, ping_time.size)),
             ),
@@ -336,4 +336,4 @@ def test_add_splitbeam_angle(sonar_model, test_path_key, raw_file_name, test_pat
 
 
 # TODO: need a test for power/angle data, with mock EchoData object
-# containing some channels with single-beam data and some channels with split-beam data
\ No newline at end of file
+# containing some channels with single-beam data and some channels with split-beam data
diff --git a/echopype/tests/echodata/test_echodata.py b/echopype/tests/echodata/test_echodata.py
index 4037fca74..e2863aa7a 100644
--- a/echopype/tests/echodata/test_echodata.py
+++ b/echopype/tests/echodata/test_echodata.py
@@ -9,7 +9,7 @@
 from zarr.errors import GroupNotFoundError
 
 import echopype
-from echopype.calibrate.env_params import EnvParams
+from echopype.calibrate.env_params_old import EnvParams
 from echopype.echodata import EchoData
 from echopype import open_converted
 
diff --git a/echopype/tests/echodata/test_echodata_simrad.py b/echopype/tests/echodata/test_echodata_simrad.py
index e8b318472..6be3b535c 100644
--- a/echopype/tests/echodata/test_echodata_simrad.py
+++ b/echopype/tests/echodata/test_echodata_simrad.py
@@ -2,7 +2,7 @@
 Tests functions contained within echodata/simrad.py
 """
 import pytest
-from echopype.echodata.simrad import retrieve_correct_beam_group, _check_input_args_combination
+from echopype.echodata.simrad import retrieve_correct_beam_group, check_input_args_combination
 
 
 @pytest.mark.parametrize(
@@ -58,7 +58,7 @@ def test_check_input_args_combination(waveform_mode: str, encode_mode: str,
         States whether pulse compression should be used
     """
 
-    _check_input_args_combination(waveform_mode, encode_mode, pulse_compression)
+    check_input_args_combination(waveform_mode, encode_mode, pulse_compression)
 
 
 def test_retrieve_correct_beam_group():
@@ -66,4 +66,3 @@ def test_retrieve_correct_beam_group():
     # TODO: create this test once we are happy with the form of retrieve_correct_beam_group
 
     pytest.skip("We need to add tests for retrieve_correct_beam_group!")
-