Refactor DOLfYN PSD Scripts to use Scipy.Signal.Welch

mhkit/acoustics/analysis.py

@@ -201,21 +201,13 @@ def sound_pressure_spectral_density(
     fs: Union[int, float],
     bin_length: Union[int, float] = 1,
     fft_length: Optional[Union[int, float]] = None,
-    rms: bool = True,
+    pct_overlap: float = 0.5,
 ) -> xr.DataArray:
     """
     Calculates the sound pressure spectral density (SPSD) from audio
     samples split into FFTs with a specified bin length in seconds,
-    using Hanning windowing with 50% overlap.
-
-    By default (`rms=True`), this function returns the mean-squared SPSD,
-    which found by scaling the total spectral power (frequency domain) with
-    the time-domain averaged mean-squared power, in accordance with
-    Parseval's theorem.
-
-    Setting `rms=False` disables this scaling and returns the
-    power spectral density of the sound pressure signal.
-    Both forms have units of [Pa^2/Hz] or [V^2/Hz].
+    using Hanning windowing with 50% overlap. Uses Welch's method to
+    average overlapping FFT windows within each bin.
 
     Parameters
     ----------
@@ -227,14 +219,13 @@ def sound_pressure_spectral_density(
         Length of time in seconds to create FFTs. Default: 1.
     fft_length: int or float, optional
         Length of FFT to use. If None, uses bin_length * fs. Default: None.
-    rms: bool
-        If True, calculates the mean-squared SPSD. Set to False to
-        calculate standard SPSD. Default: True.
+    pct_overlap: float
+        Percentage of overlap between FFT segments. Default: 0.5 (50%).
 
     Returns
     -------
-    spsd: xarray.DataArray (time, freq)
-        Spectral density [Pa^2/Hz] indexed by time and frequency
+    spsd: xarray.DataArray (time_psd, freq)
+        Spectral density [Pa^2/Hz] or [V^2/Hz] indexed by time and frequency
     """
 
     # Type checks
@@ -256,36 +247,21 @@ def sound_pressure_spectral_density(
         nfft = nbin
 
     # Use dolfyn PSD functionality
-    binner = VelBinner(n_bin=nbin, fs=fs, n_fft=nfft)
-    # Always 50% overlap if numbers reshape perfectly
-    # Mean square sound pressure
-    psd = binner.power_spectral_density(pressure, freq_units="Hz")
-    if rms:
-        # Scale PSD by mean square of original signal
-        samples = (
-            binner.reshape(pressure.values) - binner.mean(pressure.values)[:, None]
-        )
-        # mean squared pressure ("power") in time domain
-        t_power = np.sum(samples**2, axis=1) / nbin
-        # pressure ("power") in frequency domain
-        f_power = psd.sum("freq") * (fs / nbin)
-        # Adjust the amplitude of the PSD to return the mean-squared PSD
-        # based on Parseval's theorem: total energy computed in the time
-        # domain must equal the total energy computed in the frequency domain
-        psd = psd * t_power[:, None] / f_power
-        long_name = "Mean Square Sound Pressure Spectral Density"
-    else:
-        long_name = "Sound Pressure Spectral Density"
+    binner = VelBinner(fs=fs, n_bin=nbin, n_fft=nfft)
+    # Sound pressure spectral densities with 50% overlap between FFT windows
+    psd = binner.power_spectral_density(
+        pressure, freq_units="Hz", window="hann", pct_overlap=pct_overlap
+    )
 
     out = xr.DataArray(
         psd,
-        coords={"time": psd["time"], "freq": psd["freq"]},
+        coords={"time_psd": psd["time_psd"], "freq": psd["freq"]},
         attrs={
             "units": pressure.units + "^2/Hz",
-            "long_name": long_name,
+            "long_name": "Mean Square Sound Pressure Spectral Density",
             "fs": fs,
             "bin_length": bin_length,
-            "overlap": "50%",
+            "overlap": f"{pct_overlap*100}%",
             "n_fft": nfft,
         },
     )
@@ -304,7 +280,7 @@ def apply_calibration(
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral density in V^2/Hz.
     sensitivity_curve: xarray.DataArray (freq)
         Calibrated sensitivity curve in units of dB rel 1 V^2/uPa^2.
@@ -737,10 +713,11 @@ def _convert_to_band_spectral_density(
         partial_pts=partial_pts,
         weights=weights,
     )
+    time_dim = spsd.dims[0]
     out = xr.DataArray(
         band_spsd,
-        coords={"time": spsd["time"], "freq": bands[:, 1]},
-        dims=["time", "freq"],
+        coords={time_dim: spsd[time_dim], "freq": bands[:, 1]},
+        dims=[time_dim, "freq"],
         attrs=spsd.attrs,
     )
 

mhkit/acoustics/sel.py

@@ -91,7 +91,7 @@ def sound_exposure_level(
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Sound pressure spectral density in [Pa^2/Hz] with a bin length
         equal to the time over which sound exposure should be computed.
     group: str
@@ -136,9 +136,10 @@ def sound_exposure_level(
     # Sound exposure level (L_{E,p}) = (L_{p,rms} + 10log10(t))
     sel = 10 * np.log10(exposure / reference) + 10 * np.log10(spsd.attrs["bin_length"])
 
+    time_dim = spsd.dims[0]
     out = xr.DataArray(
         sel.astype(np.float32),
-        coords={"time": spsd["time"]},
+        coords={time_dim: spsd[time_dim]},
         attrs={
             "units": "dB re 1 uPa^2 s",
             "long_name": long_name,

mhkit/acoustics/spl.py

@@ -54,7 +54,7 @@ def _argument_check(spsd, fmin, fmax):
     _check_numeric(fmax, "fmax")
 
     # Ensure 'freq' and 'time' dimensions are present
-    if ("freq" not in spsd.dims) or ("time" not in spsd.dims):
+    if ("freq" not in spsd.dims) or ("time" not in spsd.dims[0]):
         raise ValueError("'spsd' must have 'time' and 'freq' as dimensions.")
 
     # Check that 'fs' (sampling frequency) is available in attributes
@@ -89,7 +89,7 @@ def sound_pressure_level(
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral density in [Pa^2/Hz]
     fmin: int
         Lower frequency band limit (lower limit of the hydrophone). Default: 10 Hz
@@ -116,9 +116,10 @@ def sound_pressure_level(
     # Mean square sound pressure level
     mspl = 10 * np.log10(pressure_squared / reference)
 
+    time_dim = spsd.dims[0]
     out = xr.DataArray(
         mspl.astype(np.float32),
-        coords={"time": spsd["time"]},
+        coords={time_dim: spsd[time_dim]},
         attrs={
             "units": "dB re 1 uPa",
             "long_name": "Sound Pressure Level",
@@ -138,7 +139,7 @@ def _band_sound_pressure_level(spsd: xr.DataArray, octave: int, base: int):
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral density in [Pa^2/Hz]
     octave: int
         Octave subdivision (1 = full octave, 3 = third-octave, etc.)
@@ -191,13 +192,14 @@ def _band_sound_pressure_level(spsd: xr.DataArray, octave: int, base: int):
             )
 
     # Mean square sound pressure level in dB rel 1 uPa
+    time_dim = spsd.dims[0]
     out_spl = xr.DataArray(
         10 * np.log10(out_sp / reference),
         coords={
-            "time": spsd["time"],
+            time_dim: spsd[time_dim],
             "freq_bins": bands[:, 1],
         },
-        dims=["time", "freq_bins"],
+        dims=[time_dim, "freq_bins"],
     )
     return out_spl
 
@@ -211,7 +213,7 @@ def third_octave_sound_pressure_level(
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral in [Pa^2/Hz].
     fmin: int
         Lower frequency band limit (lower limit of the hydrophone).
@@ -253,7 +255,7 @@ def decidecade_sound_pressure_level(
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral density in [Pa^2/Hz].
     fmin: int
         Lower frequency band limit (lower limit of the hydrophone).

mhkit/acoustics/spsdl.py

@@ -29,7 +29,7 @@ def sound_pressure_spectral_density_level(spsd: xr.DataArray) -> xr.DataArray:
 
     Parameters
     ----------
-    spsd: xarray.DataArray (time, freq)
+    spsd: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral density in Pa^2/Hz
 
     Returns
@@ -47,7 +47,7 @@ def sound_pressure_spectral_density_level(spsd: xr.DataArray) -> xr.DataArray:
 
     spsdl = xr.DataArray(
         lpf.astype(np.float32),
-        coords={"time": spsd["time"], "freq": spsd["freq"]},
+        coords=spsd.coords,
         attrs={
             "units": "dB re 1 uPa^2/Hz",
             "long_name": "Sound Pressure Spectral Density Level",
@@ -260,7 +260,7 @@ def band_aggregate(
         raise ValueError("'fmax' must be greater than 'fmin'.")
 
     # Value checks
-    if ("freq" not in spsdl.dims) or ("time" not in spsdl.dims):
+    if ("freq" not in spsdl.dims) or ("time" not in spsdl.dims[0]):
         raise ValueError("'spsdl' must have 'time' and 'freq' as dimensions.")
 
     # Validate method and get method_name and method_arg
@@ -307,7 +307,7 @@ def time_aggregate(
 
     Parameters
     ----------
-    spsdl: xarray.DataArray (time, freq)
+    spsdl: xarray.DataArray (time_psd, freq)
         Mean square sound pressure spectral density level in dB rel 1 uPa^2/Hz
     window: int
         Time in seconds to subdivide spectral density level into. Default: 60 s.
@@ -322,43 +322,31 @@ def time_aggregate(
         Time-averaged sound pressure spectral density level [dB re 1 uPa^2/Hz]
         indexed by time and frequency
     """
-    warnings.warn(
-        "The 'time_aggregate' function is deprecated and will be removed in a future release. "
-        "Please use one of the following alternatives instead to convert the SPSD to the "
-        "appropriate time-aggregated form before calculating the SPSDL using "
-        "'sound_pressure_spectral_density_level':\n"
-        "- For time-averaged SPSDLs, use the 'mhkit.acoustics.time_average' function.\n"
-        "- For time-summed SPSDLs, use the 'mhkit.acoustics.time_summation' function.\n"
-        "If you are using this function for a different purpose, please reach out to the MHKiT "
-        "developers to discuss how we can support your use case with a more specific function.",
-        DeprecationWarning,
-        stacklevel=2,
-    )
-
     # Type checks
     if not isinstance(spsdl, xr.DataArray):
         raise TypeError("'spsdl' must be an xarray.DataArray.")
     if not isinstance(window, int):
         raise TypeError("'window' must be an integer.")
     if not isinstance(method, (str, dict)):
         raise TypeError("'method' must be a string or dictionary.")
-    if "time" not in spsdl.dims:
+    if "time" not in spsdl.dims[0]:
         raise ValueError("'spsdl' must have 'time' dimension.")
 
     # Value checks
     if window <= 0:
         raise ValueError("'window' must be a positive integer.")
 
     # Ensure 'time' coordinate is of datetime64 dtype
-    if not np.issubdtype(spsdl["time"].dtype, np.datetime64):
-        raise TypeError("'spsdl['time']' must be of dtype 'datetime64'.")
+    time_dim = spsdl.dims[0]
+    if not np.issubdtype(spsdl[time_dim].dtype, np.datetime64):
+        raise TypeError("spsdl 'time' must be of dtype 'datetime64'.")
 
     # Validate method and get method_name and method_arg
     method_name, method_arg = _validate_method(method)
 
     window = np.timedelta64(window, "s")
     time_bins_lower = np.arange(
-        spsdl["time"][0].values, spsdl["time"][-1].values, window
+        spsdl[time_dim][0].values, spsdl[time_dim][-1].values, window
     )
     time_bins_upper = time_bins_lower + window
     time_bins = np.append(time_bins_lower, time_bins_upper[-1])
@@ -367,7 +355,7 @@ def time_aggregate(
     )
 
     # Use xarray binning methods
-    spsdl_group = spsdl.groupby_bins("time", time_bins, labels=center_time)
+    spsdl_group = spsdl.groupby_bins(time_dim, time_bins, labels=center_time)
 
     # Handle method being a string or a dict
     if isinstance(method, str):
@@ -388,20 +376,23 @@ def time_aggregate(
     if method == "quantile":
         out = out.drop_vars("quantile")
 
+    out = out.rename({"time_psd_bins": "time_bins"})
     return out
 
 
 def time_average(spsdl, window):
     """
-    Reorganizes spectral density level frequency tensor into time windows and takes the
-    spectral average of all of the inputs along the time dimension. This is effectively
-    time-averaging the original SPSDs.
-    Note: 'window' must be larger than the original 'bin_length' of the SPSD
+    Reorganizes spectral density level frequency tensor into time windows and computes
+    the energy-averaged SPSDL for each window. Values are converted from dB to linear
+    power, averaged across the window, then converted back to dB. This is equivalent
+    to Welch's method: it produces the same result as recomputing the SPSD from the
+    original time series using 'bin_length=window'.
+    Note: 'window' must be larger than the original 'bin_length' of the SPSD.
 
     Parameters
     ----------
     spsdl: xarray.DataArray
-        Sound pressure spectral density level with dimensions (time, freq)
+        Sound pressure spectral density level with dimensions (time_psd, freq)
     window: int
         Time in seconds to group spectral density level into.
 
@@ -413,12 +404,14 @@ def time_average(spsdl, window):
     """
 
     def spectral_average(x):
+        # time dimension name
+        time_dim = x.dims[0]
         # Convert value in decibels to absolute magnitude, still relevant to original units
         magnitude = 10 ** (x / 10)
         # Sum energy in each time bin
-        summed_magnitude = magnitude.sum("time")
+        summed_magnitude = magnitude.sum(time_dim)
         # Take average
-        average_magnitude = summed_magnitude / magnitude["time"].size
+        average_magnitude = summed_magnitude / magnitude[time_dim].size
         # Convert back to decibels
         result = 10 * np.log10(average_magnitude)
 
@@ -429,16 +422,19 @@ def spectral_average(x):
 
 def time_summation(spsdl, window):
     """
-    Reorganizes spectral density level frequency tensor into time windows and takes the
-    spectral sum of each window. This is the equivalent of recalculating the SPSD using
-    `mhkit.acoustics.sound_presssure_spectral_density` with 'bin_length=window' instead
-    of the original 'bin_length'.
-    Note: 'window' must be larger than the original 'bin_length' of the SPSD
+    Reorganizes spectral density level frequency tensor into time windows and computes
+    the spectral sum for each window. Values are converted from dB to linear power,
+    summed across the window, then converted back to dB. This represents the total
+    accumulated spectral energy within each window and is proportional to N times the
+    window-averaged SPSDL (where N is the number of input bins per window). It is NOT
+    equivalent to recomputing the SPSD with a longer 'bin_length'; use 'time_average'
+    for that purpose.
+    Note: 'window' must be larger than the original 'bin_length' of the SPSD.
 
     Parameters
     ----------
     spsdl: xarray.DataArray
-        Sound pressure spectral density level with dimensions (time, freq)
+        Sound pressure spectral density level with dimensions (time_psd, freq)
     window: int
         Time in seconds to group spectral density level into.
 
@@ -450,10 +446,12 @@ def time_summation(spsdl, window):
     """
 
     def spectral_sum(x):
+        # time dimension name
+        time_dim = x.dims[0]
         # Convert value in decibels to absolute magnitude, still relevant to original units
         magnitude = 10 ** (x / 10)
         # Sum energy in each time bin
-        summed_magnitude = magnitude.sum("time")
+        summed_magnitude = magnitude.sum(time_dim)
         # Convert back to decibels
         result = 10 * np.log10(summed_magnitude)
 

mhkit/dolfyn/adp/clean.py

@@ -5,7 +5,7 @@
 import numpy as np
 import xarray as xr
 from scipy.signal import medfilt
-from ..tools.misc import medfiltnan
+from ..tools import medfiltnan
 from ..rotate.api import rotate2
 from ..rotate.base import quaternion2orient