Repository: OSOceanAcoustics/echopype (language: python)

Issue:
Implement Noise and Attenuation Functions from Ryan et. al. 2015
From [Reducing bias due to noise and attenuation in open-ocean echo integration data](https://academic.oup.com/icesjms/article/72/8/2482/2459093), implement the following noise filters:

- [x] IN filter
- [x] AS filter
- [x] TN filter

Implementation of some (if not all) of these filters can be found in [echopy](https://github.com/open-ocean-sounding/echopy) written in pure numpy. We want to implement these filters using Xarray, and if Dask Delaying is not too involved, implement that as well. If not, we can write a separate issue for that.


Discussion / hints:
Related to #1211, since the existing noise functions will have to be renamed to something more specific.

Fix patch:
diff --git a/echopype/clean/__init__.py b/echopype/clean/__init__.py
index bfcef8186..9ce5d2908 100644
--- a/echopype/clean/__init__.py
+++ b/echopype/clean/__init__.py
@@ -1,6 +1,15 @@
-from .api import estimate_noise, remove_noise
+from .api import (
+    estimate_background_noise,
+    mask_attenuated_signal,
+    mask_impulse_noise,
+    mask_transient_noise,
+    remove_background_noise,
+)
 
 __all__ = [
-    "estimate_noise",
-    "remove_noise",
+    "estimate_background_noise",
+    "mask_attenuated_signal",
+    "mask_impulse_noise",
+    "mask_transient_noise",
+    "remove_background_noise",
 ]
diff --git a/echopype/clean/api.py b/echopype/clean/api.py
index 7ef198e8b..b771cf476 100644
--- a/echopype/clean/api.py
+++ b/echopype/clean/api.py
@@ -2,74 +2,380 @@
 Functions for reducing variabilities in backscatter data.
 """
 
+from functools import partial
+from typing import Union
+
+import numpy as np
+import xarray as xr
+
+from ..commongrid.utils import _setup_and_validate
+from ..utils.compute import _lin2log, _log2lin
+from ..utils.log import _init_logger
 from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level
-from .noise_est import NoiseEst
+from .utils import (
+    add_remove_background_noise_attrs,
+    calc_transient_noise_pooled_Sv,
+    downsample_upsample_along_depth,
+    echopy_attenuated_signal_mask,
+    echopy_impulse_noise_mask,
+)
+
+logger = _init_logger(__name__)
 
 
-def estimate_noise(ds_Sv, ping_num, range_sample_num, noise_max=None):
+def mask_transient_noise(
+    ds_Sv: xr.Dataset,
+    func: str = "nanmean",
+    depth_bin: str = "20m",
+    num_side_pings: int = 50,
+    exclude_above: Union[int, float] = 250.0,
+    transient_noise_threshold: Union[int, float] = 12.0,
+) -> xr.DataArray:
     """
-    Estimate background noise by computing mean calibrated power of a collection of pings.
+    Locate and create a mask for transient noise using a pooling comparison.
+
+    Parameters
+    ----------
+    ds_Sv : xarray.Dataset
+        Calibrated Sv data with depth data variable.
+    func: str, default `nanmean`
+        Pooling function used in the pooled Sv aggregation.
+    depth_bin : str, default `20`m
+        Pooling radius along ``depth`` in meters.
+    num_side_pings : int, default `50`
+        Number of side pings to look at for the pooling.
+    exclude_above : Union[int, float], default `250`m
+        Exclude all depth above (closer to the surface than) this value.
+    transient_noise_threshold : Union[int, float], default `10.0`dB
+        Transient noise threshold value (dB) for the pooling comparison.
+
+    Returns
+    -------
+    xr.Dataset
+        Xarray boolean array transient noise mask.
+
+    References
+    ----------
+    This function's implementation is based on the following text reference:
+
+        Ryan et al. (2015) Reducing bias due to noise and attenuation in
+        open-ocean echo integration data, ICES Journal of Marine Science, 72: 2482–2493.
+
+    Additionally, this code was derived from echopy's numpy single-channel implementation of
+    transient noise masking and translated into xarray code:
+    https://github.com/open-ocean-sounding/echopy/blob/master/echopy/processing/mask_transient.py # noqa
+    """
+    if "depth" not in ds_Sv.data_vars:
+        raise ValueError(
+            "Masking attenuated signal requires depth data variable in `ds_Sv`. "
+            "Consider adding depth with `ds_Sv = ep.consolidate.add_depth(ds_Sv)`."
+        )
+
+    # Copy `ds_Sv`
+    ds_Sv = ds_Sv.copy()
+
+    # Check for appropriate function passed in
+    if func != "nanmean" and func != "nanmedian":
+        raise ValueError("Input `func` is `nanmode`. `func` must be `nanmean` or `nanmedian`.")
+
+    # Warn when `func=nanmedian` since the sorting overhead makes it incredibly slow compared to
+    # `nanmean`.
+    logger.warning(
+        "Consider using `func=nanmean`. `func=nanmedian` is an incredibly slow operation due to "
+        "the overhead sorting."
+    )
+
+    # Setup depth bin
+    ds_Sv, depth_bin = _setup_and_validate(ds_Sv, "depth", depth_bin)
+
+    # Compute pooled Sv
+    pooled_Sv = calc_transient_noise_pooled_Sv(
+        ds_Sv, func, depth_bin, num_side_pings, exclude_above
+    )
+
+    # Compute transient noise mask
+    transient_noise_mask = ds_Sv["Sv"] - pooled_Sv > transient_noise_threshold
+
+    return transient_noise_mask
+
+
+def mask_impulse_noise(
+    ds_Sv: xr.Dataset,
+    depth_bin: str = "5m",
+    num_side_pings: int = 2,
+    impulse_noise_threshold: Union[int, float] = 10.0,
+) -> xr.DataArray:
+    """
+    Locate and create a mask for impulse noise using a ping-wise two-sided comparison.
+
+    Parameters
+    ----------
+    ds_Sv : xarray.Dataset
+        Calibrated Sv data with depth data variable.
+    depth_bin : str, default `5`m
+        Donwsampling bin size along ``depth`` in meters.
+    num_side_pings : int, default `2`
+        Number of side pings to look at for the two-side comparison.
+    impulse_noise_threshold : Union[int, float], default `10.0`dB
+        Impulse noise threshold value (dB) for the two-side comparison.
+
+    Returns
+    -------
+    xr.Dataset
+        Xarray boolean array impulse noise mask.
+
+    References
+    ----------
+    This function's implementation is based on the following text reference:
+
+        Ryan et al. (2015) Reducing bias due to noise and attenuation in
+        open-ocean echo integration data, ICES Journal of Marine Science, 72: 2482–2493.
+
+    Additionally, code was derived from echopy's numpy single-channel implementation of
+    impulse noise masking and translated into xarray code:
+    https://github.com/open-ocean-sounding/echopy/blob/master/echopy/processing/mask_impulse.py # noqa
+    """
+    if "depth" not in ds_Sv.data_vars:
+        raise ValueError(
+            "Masking attenuated signal requires depth data variable in `ds_Sv`. "
+            "Consider adding depth with `ds_Sv = ep.consolidate.add_depth(ds_Sv)`."
+        )
+
+    # Copy `ds_Sv`
+    ds_Sv = ds_Sv.copy()
+
+    # Setup and validate Sv and depth bin
+    ds_Sv, depth_bin = _setup_and_validate(ds_Sv, "depth", depth_bin)
+
+    # Downsample and Upsample Sv along depth
+    _, upsampled_Sv = downsample_upsample_along_depth(ds_Sv, depth_bin)
 
-    See ``remove_noise`` for reference.
+    # Create partial of `echopy_impulse_noise_mask`
+    partial_echopy_impulse_noise_mask = partial(
+        echopy_impulse_noise_mask,
+        num_side_pings=num_side_pings,
+        impulse_noise_threshold=impulse_noise_threshold,
+    )
+
+    # Create noise mask
+    impulse_noise_mask = xr.apply_ufunc(
+        partial_echopy_impulse_noise_mask,
+        upsampled_Sv,
+        input_core_dims=[["range_sample", "ping_time"]],
+        output_core_dims=[["range_sample", "ping_time"]],
+        dask="parallelized",
+        vectorize=True,
+        output_dtypes=[np.float64],
+    )
+
+    return impulse_noise_mask
+
+
+def mask_attenuated_signal(
+    ds_Sv: xr.Dataset,
+    upper_limit_sl: Union[int, float] = 400.0,
+    lower_limit_sl: Union[int, float] = 500.0,
+    num_side_pings: int = 15,
+    attenuation_signal_threshold: Union[int, float] = 8.0,
+) -> xr.DataArray:
+    """
+    Locate attenuated signals and create an attenuated signal mask.
+
+    Parameters
+    ----------
+    ds_Sv : xarray.Dataset
+        Calibrated Sv data with depth data variable.
+    upper_limit_sl : Union[int, float], default `400`m
+        Upper limit of deep scattering layer line (m).
+    lower_limit_sl : Union[int, float], default `500`m
+        Lower limit of deep scattering layer line (m).
+    num_side_pings : int, default `15`
+        Number of preceding & subsequent pings defining the block.
+    attenuation_signal_threshold : Union[int, float], default `8.0`dB
+        Attenuation signal threshold value (dB) for the ping-block comparison.
+
+    Returns
+    -------
+    xr.Dataset
+        Xarray boolean array attenuated signal mask.
+
+    References
+    ----------
+    This function's implementation is based on the following text reference:
+
+        Ryan et al. (2015) Reducing bias due to noise and attenuation in
+        open-ocean echo integration data, ICES Journal of Marine Science, 72: 2482–2493.
+
+    Additionally, code was derived from echopy's numpy single-channel implementation of
+    attenuation signal masking and translated into xarray code:
+    https://github.com/open-ocean-sounding/echopy/blob/master/echopy/processing/mask_attenuated.py # noqa
+    """
+    if "depth" not in ds_Sv.data_vars:
+        raise ValueError(
+            "Masking attenuated signal requires depth data variable in `ds_Sv`. "
+            "Consider adding depth with `ds_Sv = ep.consolidate.add_depth(ds_Sv)`."
+        )
+
+    # Check range values
+    if upper_limit_sl > lower_limit_sl:
+        raise ValueError("Minimum range has to be shorter than maximum range")
+
+    # Copy `ds_Sv`
+    ds_Sv = ds_Sv.copy()
+
+    # Return empty masks if searching range is outside the echosounder range
+    if (upper_limit_sl > ds_Sv["depth"].max()) or (lower_limit_sl < ds_Sv["depth"].min()):
+        attenuated_mask = xr.zeros_like(ds_Sv["Sv"], dtype=bool)
+        return attenuated_mask
+
+    # Create partial of echopy attenuation mask computation
+    partial_echopy_attenuation_mask = partial(
+        echopy_attenuated_signal_mask,
+        upper_limit_sl=upper_limit_sl,
+        lower_limit_sl=lower_limit_sl,
+        num_side_pings=num_side_pings,
+        attenuation_signal_threshold=attenuation_signal_threshold,
+    )
+
+    # Compute attenuated signal mask
+    attenuated_mask = xr.apply_ufunc(
+        partial_echopy_attenuation_mask,
+        ds_Sv["Sv"],
+        ds_Sv["depth"],
+        input_core_dims=[["ping_time", "range_sample"], ["ping_time", "range_sample"]],
+        output_core_dims=[["ping_time", "range_sample"]],
+        dask="parallelized",
+        vectorize=True,
+        output_dtypes=[bool],
+    )
+
+    return attenuated_mask
+
+
+def estimate_background_noise(
+    ds_Sv: xr.Dataset, ping_num: int, range_sample_num: int, noise_max: float = None
+) -> xr.DataArray:
+    """
+    Estimate background noise by computing mean calibrated power of a collection of pings.
 
     Parameters
     ----------
     ds_Sv : xr.Dataset
-        dataset containing ``Sv`` and ``echo_range`` [m]
+        Dataset containing ``Sv`` and ``echo_range`` [m].
     ping_num : int
-        number of pings to obtain noise estimates
+        Number of pings to obtain noise estimates
     range_sample_num : int
-        number of samples along the ``range_sample`` dimension to obtain noise estimates
-    noise_max : float
-        the upper limit for background noise expected under the operating conditions
+        Number of samples along the ``range_sample`` dimension to obtain noise estimates.
+    noise_max : float, default `None`
+        The upper limit for background noise expected under the operating conditions.
 
     Returns
     -------
     A DataArray containing noise estimated from the input ``ds_Sv``
+
+    Notes
+    -----
+    This function's implementation is based on the following text reference:
+
+        De Robertis & Higginbottom. 2007.
+        A post-processing technique to estimate the signal-to-noise ratio
+        and remove echosounder background noise.
+        ICES Journal of Marine Sciences 64(6): 1282–1291.
     """
-    noise_obj = NoiseEst(ds_Sv=ds_Sv.copy(), ping_num=ping_num, range_sample_num=range_sample_num)
-    noise_obj.estimate_noise(noise_max=noise_max)
-    return noise_obj.Sv_noise
+    # Compute transmission loss
+    spreading_loss = 20 * np.log10(ds_Sv["echo_range"].where(ds_Sv["echo_range"] >= 1, other=1))
+    absorption_loss = 2 * ds_Sv["sound_absorption"] * ds_Sv["echo_range"]
+
+    # Compute power binned averages
+    power_cal = _log2lin(ds_Sv["Sv"] - spreading_loss - absorption_loss)
+    power_cal_binned_avg = 10 * np.log10(
+        power_cal.coarsen(
+            ping_time=ping_num,
+            range_sample=range_sample_num,
+            boundary="pad",
+        ).mean()
+    )
+
+    # Compute noise
+    noise = power_cal_binned_avg.min(dim="range_sample", skipna=True)
+
+    # Align noise `ping_time` to the first index of each coarsened `ping_time` bin
+    noise = noise.assign_coords(ping_time=ping_num * np.arange(len(noise["ping_time"])))
+
+    # Limit max noise level
+    noise = noise.where(noise < noise_max, noise_max) if noise_max is not None else noise
+
+    # Upsample noise to original ping time dimension
+    Sv_noise = (
+        noise.reindex({"ping_time": power_cal["ping_time"]}, method="ffill")
+        + spreading_loss
+        + absorption_loss
+    )
+
+    return Sv_noise
 
 
 @add_processing_level("L*B")
-def remove_noise(ds_Sv, ping_num, range_sample_num, noise_max=None, SNR_threshold=3):
+def remove_background_noise(
+    ds_Sv: xr.Dataset,
+    ping_num: int,
+    range_sample_num: int,
+    noise_max: float = None,
+    SNR_threshold: float = 3.0,
+) -> xr.Dataset:
     """
     Remove noise by using estimates of background noise
     from mean calibrated power of a collection of pings.
 
-    Reference: De Robertis & Higginbottom. 2007.
-    A post-processing technique to estimate the signal-to-noise ratio
-    and remove echosounder background noise.
-    ICES Journal of Marine Sciences 64(6): 1282–1291.
-
     Parameters
     ----------
     ds_Sv : xr.Dataset
         dataset containing ``Sv`` and ``echo_range`` [m]
     ping_num : int
-        number of pings to obtain noise estimates
+        Number of pings to obtain noise estimates.
     range_sample_num : int
-        number of samples along the ``range_sample`` dimension to obtain noise estimates
-    noise_max : float
-        the upper limit for background noise expected under the operating conditions
-    SNR_threshold : float
-        acceptable signal-to-noise ratio, default to 3 dB
+        Number of samples along the ``range_sample`` dimension to obtain noise estimates.
+    noise_max : float, default `None`
+        The upper limit for background noise expected under the operating conditions.
+    SNR_threshold : float, default `3.0`
+        Acceptable signal-to-noise ratio, default to 3 dB.
 
     Returns
     -------
     The input dataset with additional variables, including
     the corrected Sv (``Sv_corrected``) and the noise estimates (``Sv_noise``)
+
+    Notes
+    -----
+    This function's implementation is based on the following text reference:
+
+        De Robertis & Higginbottom. 2007.
+        A post-processing technique to estimate the signal-to-noise ratio
+        and remove echosounder background noise.
+        ICES Journal of Marine Sciences 64(6): 1282–1291.
     """
-    noise_obj = NoiseEst(ds_Sv=ds_Sv.copy(), ping_num=ping_num, range_sample_num=range_sample_num)
-    noise_obj.remove_noise(noise_max=noise_max, SNR_threshold=SNR_threshold)
-    ds_Sv = noise_obj.ds_Sv
+    # Compute Sv_noise
+    Sv_noise = estimate_background_noise(ds_Sv, ping_num, range_sample_num, noise_max=noise_max)
+
+    # Correct Sv for noise
+    linear_corrected_Sv = _log2lin(ds_Sv["Sv"]) - _log2lin(Sv_noise)
+    corrected_Sv = _lin2log(linear_corrected_Sv.where(linear_corrected_Sv > 0, other=np.nan))
+    corrected_Sv = corrected_Sv.where(corrected_Sv - Sv_noise > SNR_threshold, other=np.nan)
 
+    # Assemble output dataset
+    ds_Sv["Sv_noise"] = Sv_noise
+    ds_Sv["Sv_noise"] = add_remove_background_noise_attrs(
+        ds_Sv["Sv_noise"], "noise", ping_num, range_sample_num, SNR_threshold, noise_max
+    )
+    ds_Sv["Sv_corrected"] = corrected_Sv
+    ds_Sv["Sv_corrected"] = add_remove_background_noise_attrs(
+        ds_Sv["Sv_corrected"], "corrected", ping_num, range_sample_num, SNR_threshold, noise_max
+    )
     prov_dict = echopype_prov_attrs(process_type="processing")
-    prov_dict["processing_function"] = "clean.remove_noise"
+    prov_dict["processing_function"] = "clean.remove_background_noise"
     ds_Sv = ds_Sv.assign_attrs(prov_dict)
 
-    # The output ds_Sv is built as a copy of the input ds_Sv, so the step below is
+    # The output `ds_Sv` is built as a copy of the input `ds_Sv`, so the step below is
     # not needed, strictly speaking. But doing makes the decorator function more generic
     ds_Sv = insert_input_processing_level(ds_Sv, input_ds=ds_Sv)
 
diff --git a/echopype/clean/noise_est.py b/echopype/clean/noise_est.py
deleted file mode 100644
index b9fbb2355..000000000
--- a/echopype/clean/noise_est.py
+++ /dev/null
@@ -1,134 +0,0 @@
-import numpy as np
-
-from ..utils import uwa
-
-
-class NoiseEst:
-    """
-    Attributes
-    ----------
-    ds_Sv : xr.Dataset
-        dataset containing ``Sv`` and ``echo_range`` [m]
-    ping_num : int
-        number of pings to obtain noise estimates
-    range_sample_num : int
-        number of samples along ``echo_range`` to obtain noise estimates
-    """
-
-    def __init__(self, ds_Sv, ping_num, range_sample_num):
-        self.ds_Sv = ds_Sv
-        self.ping_num = ping_num
-        self.range_sample_num = range_sample_num
-        self.spreading_loss = None
-        self.absorption_loss = None
-        self.Sv_noise = None
-
-        self._compute_transmission_loss()
-        self._compute_power_cal()
-
-    def _compute_transmission_loss(self):
-        """Compute transmission loss"""
-        if "sound_absorption" not in self.ds_Sv:
-            sound_absorption = uwa.calc_absorption(
-                frequency=self.ds_Sv.frequency_nominal,
-                temperature=self.ds_Sv["temperature"],
-                salinity=self.ds_Sv["salinity"],
-                pressure=self.ds_Sv["pressure"],
-            )
-        else:
-            sound_absorption = self.ds_Sv["sound_absorption"]
-
-        # Transmission loss
-        self.spreading_loss = 20 * np.log10(
-            self.ds_Sv["echo_range"].where(self.ds_Sv["echo_range"] >= 1, other=1)
-        )
-        self.absorption_loss = 2 * sound_absorption * self.ds_Sv["echo_range"]
-
-    def _compute_power_cal(self):
-        """Compute calibrated power without TVG, linear domain"""
-        self.power_cal = 10 ** (
-            (self.ds_Sv["Sv"] - self.spreading_loss - self.absorption_loss) / 10
-        )
-
-    def estimate_noise(self, noise_max=None):
-        """Estimate noise from a collected of pings
-
-        Parameters
-        ----------
-        noise_max : Union[int, float]
-            the upper limit for background noise expected under the operating conditions
-        """
-        power_cal_binned_avg = 10 * np.log10(  # binned averages of calibrated power
-            self.power_cal.coarsen(
-                ping_time=self.ping_num,
-                range_sample=self.range_sample_num,
-                boundary="pad",
-            ).mean()
-        )
-        noise = power_cal_binned_avg.min(dim="range_sample", skipna=True)
-
-        # align ping_time to first of each ping collection
-        noise["ping_time"] = self.power_cal["ping_time"][:: self.ping_num]
-
-        if noise_max is not None:
-            noise = noise.where(noise < noise_max, noise_max)  # limit max noise level
-        self.Sv_noise = (
-            noise.reindex(
-                {"ping_time": self.power_cal["ping_time"]}, method="ffill"
-            )  # forward fill empty index
-            + self.spreading_loss
-            + self.absorption_loss
-        )
-
-    def remove_noise(self, noise_max=None, SNR_threshold=3):
-        """
-        Remove noise by using estimates of background noise
-        from mean calibrated power of a collection of pings.
-
-        This method adds two data variables to the input ``ds_Sv``:
-        - corrected Sv (``Sv_corrected``)
-        - noise estimates (``Sv_noise``)
-
-        Reference: De Robertis & Higginbottom. 2007.
-        A post-processing technique to estimate the signal-to-noise ratio
-        and remove echosounder background noise.
-        ICES Journal of Marine Sciences 64(6): 1282–1291.
-
-        Parameters
-        ----------
-        noise_max : float
-            the upper limit for background noise expected under the operating conditions
-        SNR_threshold : float
-            acceptable signal-to-noise ratio, default to 3 dB
-        """
-        # Compute Sv_noise
-        self.estimate_noise(noise_max=noise_max)
-
-        # Sv corrected for noise
-        # linear domain
-        fac = 10 ** (self.ds_Sv["Sv"] / 10) - 10 ** (self.Sv_noise / 10)
-        Sv_corr = 10 * np.log10(fac.where(fac > 0, other=np.nan))
-        Sv_corr = Sv_corr.where(
-            Sv_corr - self.Sv_noise > SNR_threshold, other=np.nan
-        )  # other=-999 (from paper)
-
-        # Assemble output dataset
-        def add_attrs(sv_type, da):
-            da.attrs = {
-                "long_name": f"Volume backscattering strength, {sv_type} (Sv re 1 m-1)",
-                "units": "dB",
-                "actual_range": [
-                    round(float(da.min().values), 2),
-                    round(float(da.max().values), 2),
-                ],
-                "noise_ping_num": self.ping_num,
-                "noise_range_sample_num": self.range_sample_num,
-                "SNR_threshold": SNR_threshold,
-                "noise_max": noise_max,
-            }
-
-        self.ds_Sv["Sv_noise"] = self.Sv_noise
-        add_attrs("noise", self.ds_Sv["Sv_noise"])
-
-        self.ds_Sv["Sv_corrected"] = Sv_corr
-        add_attrs("corrected", self.ds_Sv["Sv_corrected"])
diff --git a/echopype/clean/utils.py b/echopype/clean/utils.py
new file mode 100644
index 000000000..e33910a84
--- /dev/null
+++ b/echopype/clean/utils.py
@@ -0,0 +1,238 @@
+import flox.xarray
+import numpy as np
+import xarray as xr
+
+from ..commongrid.utils import _convert_bins_to_interval_index
+from ..utils.compute import _lin2log, _log2lin
+
+
+def calc_transient_noise_pooled_Sv(
+    ds_Sv: xr.Dataset, func: str, depth_bin: float, num_side_pings: int, exclude_above: float
+) -> xr.DataArray:
+    """
+    Compute pooled Sv array for transient noise masking.
+    """
+    # Create ping time indices array
+    ping_time_indices = xr.DataArray(
+        np.arange(len(ds_Sv["ping_time"]), dtype=int),
+        dims=["ping_time"],
+        coords=[ds_Sv["ping_time"]],
+        name="ping_time_indices",
+    )
+
+    # Create NaN pooled Sv array
+    pooled_Sv = xr.full_like(ds_Sv["Sv"], np.nan)
+
+    # Set min max values
+    depth_values_min = ds_Sv["depth"].min()
+    depth_values_max = ds_Sv["depth"].max()
+    ping_time_index_min = 0
+    ping_time_index_max = len(ds_Sv["ping_time"])
+
+    # Iterate through the channel dimension
+    for channel_index in range(len(ds_Sv["channel"])):
+
+        # Set channel arrays
+        chan_Sv = ds_Sv["Sv"].isel(channel=channel_index)
+        chan_depth = ds_Sv["depth"].isel(channel=channel_index)
+
+        # Iterate through the range sample dimension
+        for range_sample_index in range(len(ds_Sv["range_sample"])):
+
+            # Iterate through the ping time dimension
+            for ping_time_index in range(len(ds_Sv["ping_time"])):
+
+                # Grab current depth
+                current_depth = ds_Sv["depth"].isel(
+                    channel=channel_index,
+                    range_sample=range_sample_index,
+                    ping_time=ping_time_index,
+                )
+
+                # Check if current value is within a valid window
+                if (
+                    (current_depth - depth_bin >= depth_values_min)
+                    & (current_depth + depth_bin <= depth_values_max)
+                    & (current_depth - depth_bin >= exclude_above)
+                    & (ping_time_index - num_side_pings >= ping_time_index_min)
+                    & (ping_time_index + num_side_pings <= ping_time_index_max)
+                ):
+
+                    # Compute aggregate window Sv value
+                    window_mask = (
+                        (current_depth - depth_bin <= chan_depth)
+                        & (chan_depth <= current_depth + depth_bin)
+                        & (ping_time_index - num_side_pings <= ping_time_indices)
+                        & (ping_time_indices <= ping_time_index + num_side_pings)
+                    )
+                    window_Sv = chan_Sv.where(window_mask, other=np.nan).pipe(_log2lin)
+                    aggregate_window_Sv_value = window_Sv.pipe(
+                        np.nanmean if func == "nanmean" else np.nanmedian
+                    )
+                    aggregate_window_Sv_value = _lin2log(aggregate_window_Sv_value)
+
+                    # Put aggregate value in pooled Sv array
+                    pooled_Sv[
+                        dict(
+                            channel=channel_index,
+                            range_sample=range_sample_index,
+                            ping_time=ping_time_index,
+                        )
+                    ] = aggregate_window_Sv_value
+
+    return pooled_Sv
+
+
+def downsample_upsample_along_depth(ds_Sv: xr.Dataset, depth_bin: float) -> xr.DataArray:
+    """
+    Downsample and upsample Sv to mimic what was done in echopy impulse
+    noise masking.
+    """
+    # Validate and compute range interval
+    depth_min = ds_Sv["depth"].min()
+    depth_max = ds_Sv["depth"].max()
+    range_interval = np.arange(depth_min, depth_max + depth_bin, depth_bin)
+    range_interval = _convert_bins_to_interval_index(range_interval)
+
+    # Downsample Sv along range sample
+    downsampled_Sv = flox.xarray.xarray_reduce(
+        ds_Sv["Sv"].pipe(_log2lin),
+        ds_Sv["channel"],
+        ds_Sv["ping_time"],
+        ds_Sv["depth"],
+        expected_groups=(None, None, range_interval),
+        isbin=[False, False, True],
+        method="map-reduce",
+        func="nanmean",
+        skipna=True,
+    ).pipe(_lin2log)
+
+    # Assign a depth bin index to each Sv depth value
+    depth_bin_assignment = xr.DataArray(
+        np.digitize(
+            ds_Sv["depth"], [interval.left for interval in downsampled_Sv["depth_bins"].data]
+        ),
+        dims=["channel", "ping_time", "range_sample"],
+    )
+
+    # Initialize upsampled Sv
+    upsampled_Sv = ds_Sv["Sv"].copy()
+
+    # Iterate through all channels
+    for channel_index in range(len(depth_bin_assignment["channel"])):
+        # Iterate through all ping times
+        for ping_time_index in range(len(depth_bin_assignment["ping_time"])):
+            # Get unique range sample values along a single ping from the digitized depth array:
+            # NOTE: The unique index corresponds to the first unique value's position which in
+            # turn corresponds to the first range sample value contained in each depth bin.
+            _, unique_range_sample_indices = np.unique(
+                depth_bin_assignment.isel(channel=channel_index, ping_time=ping_time_index).data,
+                return_index=True,
+            )
+
+            # Select a single ping downsampled Sv vector
+            subset_downsampled_Sv = downsampled_Sv.isel(
+                channel=channel_index, ping_time=ping_time_index
+            )
+
+            # Substitute depth bin coordinate in the downsampled Sv to be the range sample value
+            # corresponding to the first element (lowest depth value) of each depth bin, and rename
+            # `depth_bin` coordinate to `range_sample`.
+            subset_downsampled_Sv = subset_downsampled_Sv.assign_coords(
+                {"depth_bins": unique_range_sample_indices}
+            ).rename({"depth_bins": "range_sample"})
+
+            # Upsample via `reindex` `ffill`
+            upsampled_Sv[dict(channel=channel_index, ping_time=ping_time_index)] = (
+                subset_downsampled_Sv.reindex(
+                    {"range_sample": ds_Sv["range_sample"]}, method="ffill"
+                )
+            )
+
+    return downsampled_Sv, upsampled_Sv
+
+
+def echopy_impulse_noise_mask(
+    Sv: np.ndarray, num_side_pings: int, impulse_noise_threshold: float
+) -> np.ndarray:
+    """Single-channel impulse noise mask computation from echopy."""
+    # Construct the two ping side-by-side comparison arrays
+    dummy = np.zeros((Sv.shape[0], num_side_pings)) * np.nan
+    comparison_forward = Sv - np.c_[Sv[:, num_side_pings:], dummy]
+    comparison_backward = Sv - np.c_[dummy, Sv[:, 0:-num_side_pings]]
+    comparison_forward[np.isnan(comparison_forward)] = np.inf
+    comparison_backward[np.isnan(comparison_backward)] = np.inf
+
+    # Create mask by checking if comparison arrays are above `impulse_noise_threshold`
+    maskf = comparison_forward > impulse_noise_threshold
+    maskb = comparison_backward > impulse_noise_threshold
+    mask = maskf & maskb
+
+    return mask
+
+
+def echopy_attenuated_signal_mask(
+    Sv: np.ndarray,
+    depth: np.ndarray,
+    upper_limit_sl: float,
+    lower_limit_sl: float,
+    num_side_pings: int,
+    attenuation_signal_threshold: float,
+) -> np.ndarray:
+    """Single-channel attenuated signal mask computation from echopy."""
+    # Initialize mask
+    attenuated_mask = np.zeros(Sv.shape, dtype=bool)
+
+    for ping_time_idx in range(Sv.shape[0]):
+
+        # Find indices for upper and lower SL limits
+        up = np.argmin(abs(depth[ping_time_idx, :] - upper_limit_sl))
+        lw = np.argmin(abs(depth[ping_time_idx, :] - lower_limit_sl))
+
+        # Mask when attenuation masking is feasible
+        if not (
+            (ping_time_idx - num_side_pings < 0)
+            | (ping_time_idx + num_side_pings > Sv.shape[0] - 1)
+            | np.all(np.isnan(Sv[ping_time_idx, up:lw]))
+        ):
+            # Compare ping and block medians, and mask ping if difference greater than
+            # threshold.
+            pingmedian = _lin2log(np.nanmedian(_log2lin(Sv[ping_time_idx, up:lw])))
+            blockmedian = _lin2log(
+                np.nanmedian(
+                    _log2lin(
+                        Sv[
+                            (ping_time_idx - num_side_pings) : (ping_time_idx + num_side_pings),
+                            up:lw,
+                        ]
+                    )
+                )
+            )
+            if (pingmedian - blockmedian) < attenuation_signal_threshold:
+                attenuated_mask[ping_time_idx, :] = True
+
+    return attenuated_mask
+
+
+def add_remove_background_noise_attrs(
+    da: xr.DataArray,
+    sv_type: str,
+    ping_num: int,
+    range_sample_num: int,
+    SNR_threshold: float,
+    noise_max: float,
+) -> xr.DataArray:
+    """Add attributes to a `remove_background_noise` data array."""
+    da.attrs = {
+        "long_name": f"Volume backscattering strength, {sv_type} (Sv re 1 m-1)",
+        "units": "dB",
+        "actual_range": [
+            round(float(da.min().values), 2),
+            round(float(da.max().values), 2),
+        ],
+        "noise_ping_num": ping_num,
+        "noise_range_sample_num": range_sample_num,
+        "SNR_threshold": SNR_threshold,
+        "noise_max": noise_max,
+    }
+    return da


Test patch:
diff --git a/echopype/tests/clean/test_noise.py b/echopype/tests/clean/test_noise.py
index 92423cced..9364973fd 100644
--- a/echopype/tests/clean/test_noise.py
+++ b/echopype/tests/clean/test_noise.py
@@ -2,11 +2,509 @@
 import xarray as xr
 import echopype as ep
 import pytest
-from numpy.random import default_rng
 
+from echopype.clean.utils import calc_transient_noise_pooled_Sv, downsample_upsample_along_depth
+from echopype.utils.compute import _lin2log, _log2lin
 
-def test_remove_noise():
-    """Test remove_noise on toy data"""
+
+@pytest.mark.integration
+def test_mask_functions_with_no_depth():
+    """Test mask functions when the depth variable is not within `ds_Sv`."""
+    # Open raw and calibrate
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+
+    # `depth` is not contained in `ds_Sv`. Ensure that `ValueError` is raised
+    # for all masking functions.
+    with pytest.raises(ValueError):
+        ep.clean.mask_attenuated_signal(ds_Sv)
+    with pytest.raises(ValueError):
+        ep.clean.mask_impulse_noise(ds_Sv)
+    with pytest.raises(ValueError):
+        ep.clean.mask_transient_noise(ds_Sv)
+
+
+@pytest.mark.integration
+def test_mask_functions_dimensions():
+    """Test mask functions' output dimensions."""
+    # Open raw, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Select Sv subset
+    ds_Sv = ds_Sv.isel(ping_time=slice(294, 300), range_sample=slice(1794,1800))
+
+    # Compute masks and check that dimensions match `ds_Sv`
+    attenuated_signal_mask = ep.clean.mask_attenuated_signal(ds_Sv)
+    impulse_noise_mask = ep.clean.mask_impulse_noise(ds_Sv)
+    transient_noise_mask = ep.clean.mask_transient_noise(ds_Sv)
+    for mask in [attenuated_signal_mask, impulse_noise_mask, transient_noise_mask]:
+        assert ds_Sv["channel"].equals(mask["channel"])
+        assert np.allclose(ds_Sv["range_sample"].data, mask["range_sample"].data)
+        assert ds_Sv["ping_time"].equals(mask["ping_time"])
+
+
+@pytest.mark.integration
+def test_transient_mask_noise_func_error_and_warnings(caplog):
+    """Check if appropriate warnings and errors are raised for transient noise mask func input."""
+    # Open raw, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Select Sv subset
+    ds_Sv = ds_Sv.isel(ping_time=slice(290, 300), range_sample=slice(1794,1800))
+
+    # Set window args
+    depth_bin="0.2m"
+    num_side_pings=2
+    exclude_above = 250
+
+    ### Check for `nanmedian` warning:
+
+    # Turn on logger verbosity
+    ep.utils.log.verbose(override=False)
+
+    # Compute transient noise mask
+    ep.clean.mask_transient_noise(
+        ds_Sv,
+        depth_bin=depth_bin,
+        num_side_pings=num_side_pings,
+        exclude_above=exclude_above,
+        func="nanmedian",
+    )
+
+    # Check for expected warning
+    expected_warning = (
+        "Consider using `func=nanmean`. `func=nanmedian` is an incredibly slow operation due to "
+        "the overhead sorting."
+    )
+    assert any(expected_warning in record.message for record in caplog.records)
+
+    # Turn off logger verbosity
+    ep.utils.log.verbose(override=True)
+
+    # Check for func value error:
+    with pytest.raises(ValueError, match="Input `func` is `nanmode`. `func` must be `nanmean` or `nanmedian`."):
+        ep.clean.mask_transient_noise(
+            ds_Sv,
+            depth_bin=depth_bin,
+            num_side_pings=num_side_pings,
+            exclude_above=exclude_above,
+            func="nanmode",
+        )
+
+
+@pytest.mark.integration
+@pytest.mark.parametrize(
+    ("chunk", "func"),
+    [
+        (False, "nanmean"),
+        (True, "nanmean"),
+        (False, "nanmedian"),
+        (True, "nanmedian"),
+    ],
+)
+def test_calc_transient_noise_pooled_Sv_values(chunk, func):
+    """
+    Manually check if the pooled Sv for transient noise masking contains 
+    the correct nan boundary and the correct bin aggregate values.
+    """
+    # Open raw, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Select Sv subset
+    ds_Sv = ds_Sv.isel(ping_time=slice(294, 300), range_sample=slice(1794,1800))
+
+    if chunk:
+        # Chunk calibrated Sv
+        ds_Sv = ds_Sv.chunk("auto")
+
+    # Set window args
+    depth_bin = 0.2 # depth values ~0.2m apart per range sample
+    num_side_pings=2
+    exclude_above = 250
+
+    # Compute pooled Sv
+    pooled_Sv = calc_transient_noise_pooled_Sv(
+        ds_Sv, func, depth_bin, num_side_pings, exclude_above
+    ).compute()
+
+    # Compute min and max values
+    depth_values_min = ds_Sv["depth"].min()
+    depth_values_max = ds_Sv["depth"].max()
+    ping_time_index_min = 0
+    ping_time_index_max = len(ds_Sv["ping_time"])
+
+    # Create ping time indices array
+    ping_time_indices = xr.DataArray(
+        np.arange(len(ds_Sv["ping_time"]), dtype=int),
+        dims=["ping_time"],
+        coords=[ds_Sv["ping_time"]],
+        name="ping_time_indices",
+    )
+
+    # Check appropriate NaN boundaries
+    within_mask = (
+        (ds_Sv["depth"] - depth_bin >= depth_values_min) &
+        (ds_Sv["depth"] + depth_bin <= depth_values_max) &
+        (ds_Sv["depth"] - depth_bin >= exclude_above) &
+        (ping_time_indices - num_side_pings >= ping_time_index_min) &
+        (ping_time_indices + num_side_pings <= ping_time_index_max)
+    )
+    unique_pool_boundary_values = np.unique(
+        pooled_Sv.where(~within_mask)
+    )
+
+    # Check that NaN is the only pool boundary unique value
+    assert np.isclose(unique_pool_boundary_values, np.array([np.nan]), equal_nan=True)
+
+    # Compute Sv
+    ds_Sv["Sv"] = ds_Sv["Sv"].compute()
+
+    # Check correct binning and aggregation values
+    for channel_index in range(len(ds_Sv["channel"])):
+        for ping_time_index in range(len(ds_Sv["ping_time"])):
+            for range_sample_index in range(len(ds_Sv["range_sample"])):
+                # Grab pooled value
+                pooled_value = pooled_Sv.isel(
+                    channel=channel_index,
+                    ping_time=ping_time_index,
+                    range_sample=range_sample_index
+                ).data
+                if not np.isnan(pooled_value):
+                    assert np.isclose(
+                        pooled_value,
+                        # Compute pooled value using the fact that the depth bin covers
+                        # 3 depth values (1 below, middle, 1 above) and that the ping bin
+                        # covers 5 ping values (2 below, middle, 2 above).
+                        _lin2log(
+                            _log2lin(
+                                ds_Sv["Sv"].isel(
+                                    channel=channel_index,
+                                    ping_time=slice(ping_time_index-2, ping_time_index+3),
+                                    range_sample=slice(range_sample_index-1, range_sample_index+2)
+                                )
+                            ).pipe(
+                                np.nanmean if func == "nanmean" else np.nanmedian
+                            )
+                        ),
+                        rtol=1e-10,
+                        atol=1e-10,
+                    )
+
+
+@pytest.mark.integration
+@pytest.mark.parametrize(
+    ("chunk", "func"),
+    [
+        (False, "nanmean"),
+        (True, "nanmean"),
+        (False, "nanmedian"),
+        (True, "nanmedian"),
+    ],
+)
+def test_transient_noise_mask_values(chunk, func):
+    """Manually check if impulse noise mask removes transient noise values."""
+    # Open raw, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Select Sv subset
+    ds_Sv = ds_Sv.isel(ping_time=slice(294, 300), range_sample=slice(1794,1800))
+
+    if chunk:
+        # Chunk calibrated Sv
+        ds_Sv = ds_Sv.chunk("auto")
+
+    # Set window args
+    depth_bin = "0.2m" # depth values ~0.2m apart per range sample
+    num_side_pings=2
+    exclude_above = 250
+    transient_noise_threshold = 12
+
+    # Compute transient noise mask
+    transient_noise_mask = ep.clean.mask_transient_noise(
+        ds_Sv, func, depth_bin, num_side_pings, exclude_above, transient_noise_threshold
+    ).compute()
+
+    # Remove transient noise from Sv
+    ds_Sv["Sv_cleaned_of_transient_noise"] = xr.where(
+        transient_noise_mask,
+        np.nan,
+        ds_Sv["Sv"]
+    ).compute()
+
+    # Compute Sv
+    ds_Sv["Sv"] = ds_Sv["Sv"].compute()
+
+    # Check if transient noise values have been removed:
+    for channel_index in range(len(ds_Sv["channel"])):
+        for ping_time_index in range(len(ds_Sv["ping_time"])):
+            for range_sample_index in range(len(ds_Sv["range_sample"])):
+                # Grab cleaned value
+                Sv_value = ds_Sv["Sv_cleaned_of_transient_noise"].isel(
+                    channel=channel_index,
+                    ping_time=ping_time_index,
+                    range_sample=range_sample_index
+                ).data
+                # Compute pooled value using the fact that the depth bin covers
+                # 3 depth values (1 below, middle, 1 above) and that the ping bin
+                # covers 5 ping values (2 below, middle, 2 above).
+                pooled_value = _lin2log(
+                    _log2lin(
+                        ds_Sv["Sv"].isel(
+                            channel=channel_index,
+                            ping_time=slice(ping_time_index-2, ping_time_index+3),
+                            range_sample=slice(range_sample_index-1, range_sample_index+2)
+                        )
+                    ).pipe(
+                        np.nanmean if func == "nanmean" else np.nanmedian
+                    )
+                )
+                # Check negation of transient noise condition only when both values are not NaN:
+                if not np.isnan(Sv_value) and not np.isnan(pooled_value):
+                    assert Sv_value - pooled_value <= transient_noise_threshold
+
+
+@pytest.mark.integration
+@pytest.mark.parametrize(
+    ("chunk"),
+    [
+        (False),
+        (True),
+    ],
+)
+def test_downsample_upsample_along_depth(chunk):
+    """Test downsample bins and upsample repeating values"""
+    # Open raw, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Subset `ds_Sv`
+    ds_Sv = ds_Sv.isel(ping_time=slice(100, 120), range_sample=slice(0, 100))
+
+    if chunk:
+        # Chunk calibrated Sv
+        ds_Sv = ds_Sv.chunk("auto")
+
+    # Run downsampling and upsampling
+    downsampled_Sv, upsampled_Sv = downsample_upsample_along_depth(ds_Sv, 2)
+
+    # Compute DataArrays
+    downsampled_Sv = downsampled_Sv.compute()
+    upsampled_Sv = upsampled_Sv.compute()
+    original_resolution_depth = ds_Sv["depth"].compute()
+
+    # Check for appropriate binning behavior
+    # Test every channel
+    for channel_index in range(len(downsampled_Sv["channel"])):
+        # Test every depth bin
+        for depth_bin_index in range(len(downsampled_Sv["depth_bins"])):
+            # Test every ping time
+            for ping_time_index in range(len(downsampled_Sv["ping_time"])):
+                # Check that manual and flox downsampled bin Sv are equal
+                flox_downsampled_bin_Sv = downsampled_Sv.isel(
+                    channel=channel_index, ping_time=ping_time_index, depth_bins=depth_bin_index
+                ).data
+                flox_downsampled_bin_Sv_indices = np.where(
+                    upsampled_Sv.isel(channel=channel_index, ping_time=ping_time_index).data == flox_downsampled_bin_Sv
+                )[0]
+                manual_downsampled_bin_Sv = _lin2log(
+                    np.nanmean(
+                        ds_Sv["Sv"].compute().isel(
+                        channel=channel_index, ping_time=ping_time_index, range_sample=flox_downsampled_bin_Sv_indices
+                        ).pipe(_log2lin)
+                    )
+                )
+                assert np.isclose(manual_downsampled_bin_Sv, flox_downsampled_bin_Sv, atol=1e-10, rtol=1e-10)
+  
+                # Check that depth bins encapsulated the correct original resolution depth values
+                manual_depth_array = original_resolution_depth.isel(
+                    channel=channel_index, ping_time=ping_time_index, range_sample=flox_downsampled_bin_Sv_indices
+                ).data
+                flox_depth_bin = downsampled_Sv["depth_bins"].data[depth_bin_index]
+                for manual_depth in manual_depth_array:
+                    if not np.isnan(manual_depth):
+                        assert flox_depth_bin.left <= manual_depth < flox_depth_bin.right
+
+
+@pytest.mark.integration
+@pytest.mark.parametrize(
+    ("chunk"),
+    [
+        (False),
+        (True),
+    ],
+)
+def test_impulse_noise_mask_values(chunk):
+    """Manually check if impulse noise mask removes impulse noise values."""
+    # Open raw, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR230-D20091215-T121917.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    if chunk:
+        # Chunk calibrated Sv
+        ds_Sv = ds_Sv.chunk("auto")
+
+    # Subset `ds_Sv`
+    ds_Sv = ds_Sv.isel(ping_time=slice(100, 120), range_sample=slice(0, 100))
+
+    # Create impulse noise mask
+    impulse_noise_mask = ep.clean.mask_impulse_noise(ds_Sv, "2m")
+
+    # Compute upsampled data
+    _, upsampled_Sv = downsample_upsample_along_depth(ds_Sv, 2)
+    upsampled_Sv = upsampled_Sv.compute()
+
+    # Remove impulse noise from Sv
+    ds_Sv["upsampled_Sv_cleaned_of_impulse_noise"] = xr.where(
+        impulse_noise_mask,
+        np.nan,
+        upsampled_Sv
+    ).compute()
+
+    # Iterate through every channel
+    for channel_index in range(len(ds_Sv["channel"])):
+        # Iterate through every range sample
+        for range_sample_index in range(len(ds_Sv["range_sample"])):
+            # Iterate through every valid ping time
+            for ping_time_index in range(2, len(ds_Sv["ping_time"]) - 2):
+                # Grab range sample row array
+                row_array = ds_Sv["upsampled_Sv_cleaned_of_impulse_noise"].isel(
+                    channel=channel_index,
+                    ping_time=slice(ping_time_index - 2, ping_time_index + 3),
+                    range_sample=range_sample_index
+                ).data
+                # Compute left and right subtraction values
+                left_subtracted_value = row_array[2] - row_array[0]
+                right_subtracted_value = row_array[2] - row_array[4]
+                # Check negation of impulse condition if middle array value and subtraction values are not NaN
+                if not (
+                    np.isnan(row_array[2]) or np.isnan(left_subtracted_value) or np.isnan(right_subtracted_value)
+                ):
+                    assert (left_subtracted_value <= 10.0 or right_subtracted_value <= 10.0)
+
+
+@pytest.mark.integration
+def test_mask_attenuated_signal_limit_error():
+    """Test `mask_attenuated_signal` limit error."""
+    # Parse, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Attempt to create mask with `upper_limit_sl > lower_limit_sl`
+    with pytest.raises(ValueError):
+        ep.clean.mask_attenuated_signal(
+            ds_Sv,
+            upper_limit_sl=180,
+            lower_limit_sl=170,
+        )
+
+
+@pytest.mark.integration
+def test_mask_attenuated_signal_outside_searching_range():
+    """Test `mask_attenuated_signal` values errors."""
+    # Parse, calibrate, and add_depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    # Create mask
+    upper_limit_sl, lower_limit_sl, num_side_pings, attenuation_signal_threshold = 1800, 2800, 15, -6 # units: (m, m, pings, dB)
+    attenuated_mask = ep.clean.mask_attenuated_signal(
+        ds_Sv,
+        upper_limit_sl,
+        lower_limit_sl,
+        num_side_pings,
+        attenuation_signal_threshold
+    )
+    
+    # Check outputs
+    assert np.allclose(attenuated_mask, xr.zeros_like(ds_Sv["Sv"], dtype=bool))
+
+
+@pytest.mark.integration
+@pytest.mark.parametrize(
+    ("chunk"),
+    [
+        (False),
+        (True),
+    ],
+)
+def test_mask_attenuated_signal_against_echopy(chunk):
+    """Test `attenuated_signal` to see if Echopype output matches echopy output mask."""
+    # Parse, calibrate, and add depth
+    ed = ep.open_raw(
+        "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645.raw",
+        sonar_model="EK60"
+    )
+    ds_Sv = ep.calibrate.compute_Sv(ed)
+    ds_Sv = ep.consolidate.add_depth(ds_Sv)
+
+    if chunk:
+        # Chunk dataset
+        ds_Sv = ds_Sv.chunk("auto")
+
+    # Create mask
+    upper_limit_sl, lower_limit_sl, num_side_pings, attenuation_signal_threshold = 180, 280, 30, -6 # units: (m, m, pings, dB)
+    attenuated_mask = ep.clean.mask_attenuated_signal(
+        ds_Sv,
+        upper_limit_sl,
+        lower_limit_sl,
+        num_side_pings,
+        attenuation_signal_threshold
+    )
+
+    # Grab echopy attenuated signal mask
+    echopy_attenuated_mask = xr.open_dataset(
+        "echopype/test_data/ek60/from_echopy/JR161-D20061118-T010645_echopy_attenuated_masks.zarr",
+        engine="zarr"
+    )
+
+    # Check that Echopype 38kHz mask matches echopy mask
+    assert np.allclose(
+        echopy_attenuated_mask["attenuated_mask"],
+        attenuated_mask.isel(channel=0).transpose("range_sample", "ping_time")
+    )
+
+
+def test_remove_background_noise():
+    """Test remove_background_noise on toy data"""
 
     # Parameters for fake data
     nchan, npings, nrange_samples = 1, 10, 100
@@ -40,7 +538,7 @@ def test_remove_noise():
     )
     ds_Sv = ds_Sv.assign(sound_absorption=0.001)
     # Run noise removal
-    ds_Sv = ep.clean.remove_noise(
+    ds_Sv = ep.clean.remove_background_noise(
         ds_Sv, ping_num=2, range_sample_num=5, SNR_threshold=0
     )
 
@@ -57,7 +555,7 @@ def test_remove_noise():
     data = np.random.normal(
         loc=-100, scale=2, size=(nchan, npings, nrange_samples)
     )
-    # Make Dataset to pass into remove_noise
+    # Make Dataset to pass into remove_background_noise
     Sv = xr.DataArray(
         data,
         coords=[
@@ -77,7 +575,7 @@ def test_remove_noise():
     )
     ds_Sv = ds_Sv.assign(sound_absorption=0.001)
     # Run noise removal
-    ds_Sv = ep.clean.remove_noise(
+    ds_Sv = ep.clean.remove_background_noise(
         ds_Sv, ping_num=2, range_sample_num=5, SNR_threshold=0
     )
     null = ds_Sv.Sv_corrected.isnull()
@@ -86,13 +584,3 @@ def test_remove_noise():
         np.count_nonzero(null.isel(channel=0, range_sample=slice(None, 50)))
         == 6
     )
-
-
-def test_remove_noise_no_sound_absorption():
-    """
-    Tests remove_noise on toy data that does
-    not have sound absorption as a variable.
-    """
-
-    pytest.xfail(f"Tests for remove_noise have not been implemented" +
-                 " when no sound absorption is provided!")