Acoustics SPL bugfix

mhkit/acoustics/__init__.py

@@ -1,17 +1,17 @@
 """
 The passive acoustics module provides a set of functions
-for analyzing and visualizing passive acoustic monitoring 
+for analyzing and visualizing passive acoustic monitoring
 data deployed in water bodies. This package reads in raw
-*.wav* files and conducts basic acoustics analysis and 
+*.wav* files and conducts basic acoustics analysis and
 visualization.
 
 To start using the module, import it directly from MHKiT:
 ``from mhkit import acoustics``. The analysis functions
-are available directly from the main import, while the 
-I/O and graphics submodules are available from 
+are available directly from the main import, while the
+I/O and graphics submodules are available from
 ``acoustics.io`` and  ``acoustics.graphics``, respectively.
-The base functions are intended to be used on top of the I/O submodule, and 
-include functionality to calibrate data, create spectral densities, sound 
+The base functions are intended to be used on top of the I/O submodule, and
+include functionality to calibrate data, create spectral densities, sound
 pressure levels, and time or band aggregate spectral data.
 """
 

mhkit/acoustics/analysis.py

@@ -337,7 +337,7 @@ def sound_pressure_spectral_density_level(spsd: xr.DataArray) -> xr.DataArray:
 
 
 def _validate_method(
-    method: Union[str, Dict[str, Union[float, int]]]
+    method: Union[str, Dict[str, Union[float, int]]],
 ) -> Tuple[str, Optional[Union[float, int]]]:
     """
     Validates the 'method' parameter and returns the method name and its argument (if any)
@@ -473,6 +473,44 @@ def _validate_method(
     return method_name, method_arg
 
 
+def _create_frequency_bands(octave, fmin, fmax):
+    """
+    Calculates frequency bands based on the specified octave, minimum and
+    maximum frequency limits.
+
+    Parameters
+    ----------
+    octave: int
+        Octave to subdivide spectral density level by.
+    fmin : int, optional
+        Lower frequency band limit (lower limit of the hydrophone). Default is 10 Hz.
+    fmax : int, optional
+        Upper frequency band limit (Nyquist frequency). Default is 100,000 Hz.
+
+    Returns
+    -------
+    octave_bins: numpy.array
+        Array of octave bin edges
+    band: dict(str, numpy.array)
+        Dictionary containing the frequency band edges and center frequency
+    """
+
+    bandwidth = 2 ** (1 / octave)
+    half_bandwidth = 2 ** (1 / (octave * 2))
+
+    band = {}
+    band["center_freq"] = 10 ** np.arange(
+        np.log10(fmin),
+        np.log10(fmax * bandwidth),
+        step=np.log10(bandwidth),
+    )
+    band["lower_limit"] = band["center_freq"] / half_bandwidth
+    band["upper_limit"] = band["center_freq"] * half_bandwidth
+    octave_bins = np.append(band["lower_limit"], band["upper_limit"][-1])
+
+    return octave_bins, band
+
+
 def band_aggregate(
     spsdl: xr.DataArray,
     octave: int = 3,
@@ -531,18 +569,7 @@ def band_aggregate(
     fn = spsdl["freq"].max().values
     fmax = _fmax_warning(fn, fmax)
 
-    bandwidth = 2 ** (1 / octave)
-    half_bandwidth = 2 ** (1 / (octave * 2))
-
-    band = {}
-    band["center_freq"] = 10 ** np.arange(
-        np.log10(fmin),
-        np.log10(fmax * bandwidth),
-        step=np.log10(bandwidth),
-    )
-    band["lower_limit"] = band["center_freq"] / half_bandwidth
-    band["upper_limit"] = band["center_freq"] * half_bandwidth
-    octave_bins = np.append(band["lower_limit"], band["upper_limit"][-1])
+    octave_bins, band = _create_frequency_bands(octave, fmin, fmax)
 
     # Use xarray binning methods
     spsdl_group = spsdl.groupby_bins("freq", octave_bins, labels=band["center_freq"])
@@ -732,8 +759,7 @@ def sound_pressure_level(
 
 def _band_sound_pressure_level(
     spsd: xr.DataArray,
-    bandwidth: int,
-    half_bandwidth: int,
+    octave: int,
     fmin: int = 10,
     fmax: int = 100000,
 ) -> xr.DataArray:
@@ -744,10 +770,8 @@ def _band_sound_pressure_level(
     ----------
     spsd: xarray.DataArray (time, freq)
         Mean square sound pressure spectral density.
-    bandwidth : int or float
-        Bandwidth to average over.
-    half_bandwidth : int or float
-        Half-bandwidth, used to set upper and lower bandwidth limits.
+    octave: int
+        Octave to subdivide spectral density level by.
     fmin : int, optional
         Lower frequency band limit (lower limit of the hydrophone). Default is 10 Hz.
     fmax : int, optional
@@ -763,10 +787,8 @@ def _band_sound_pressure_level(
     # Type checks
     if not isinstance(spsd, xr.DataArray):
         raise TypeError("'spsd' must be an xarray.DataArray.")
-    if not isinstance(bandwidth, (int, float)):
-        raise TypeError("'bandwidth' must be a numeric type (int or float).")
-    if not isinstance(half_bandwidth, (int, float)):
-        raise TypeError("'half_bandwidth' must be a numeric type (int or float).")
+    if not isinstance(octave, int) or (octave <= 0):
+        raise TypeError("'octave' must be a positive integer.")
     if not isinstance(fmin, int):
         raise TypeError("'fmin' must be an integer.")
     if not isinstance(fmax, int):
@@ -795,28 +817,35 @@ def _band_sound_pressure_level(
     # Reference value of sound pressure
     reference = 1e-12  # Pa^2, = 1 uPa^2
 
-    band = {}
-    band["center_freq"] = 10 ** np.arange(
-        np.log10(fmin),
-        np.log10(fmax * bandwidth),
-        step=np.log10(bandwidth),
-    )
-    band["lower_limit"] = band["center_freq"] / half_bandwidth
-    band["upper_limit"] = band["center_freq"] * half_bandwidth
-    octave_bins = np.append(band["lower_limit"], band["upper_limit"][-1])
+    _, band = _create_frequency_bands(octave, fmin, fmax)
 
     # Manual trapezoidal rule to get Pa^2
     pressure_squared = xr.DataArray(
         coords={"time": spsd["time"], "freq_bins": band["center_freq"]},
         dims=["time", "freq_bins"],
     )
     for i, key in enumerate(band["center_freq"]):
-        band_min = octave_bins[i]
-        band_max = octave_bins[i + 1]
-        pressure_squared.loc[{"freq_bins": key}] = np.trapz(
-            spsd.sel(freq=slice(band_min, band_max)),
-            spsd["freq"].sel(freq=slice(band_min, band_max)),
-        )
+        # Min and max band limits
+        band_range = [band["lower_limit"][i], band["upper_limit"][i]]
+
+        # Integrate spectral density by frequency
+        x = spsd["freq"].sel(freq=slice(*band_range))
+        if len(x) < 2:
+            # Interpolate between band frequencies if width is narrow
+            bandwidth = band_range[1] / band_range[0]
+            # Use smaller set of dataset to speed up interpolation
+            spsd_slc = spsd.sel(
+                freq=slice(
+                    None,  # Only happens at low frequency
+                    band_range[1] * bandwidth * 2,
+                )
+            )
+            spsd_slc = spsd_slc.interp(freq=band_range)
+            x = band_range
+        else:
+            spsd_slc = spsd.sel(freq=slice(*band_range))
+
+        pressure_squared.loc[{"freq_bins": key}] = np.trapz(spsd_slc, x)
 
     # Mean square sound pressure level in dB rel 1 uPa
     mspl = 10 * np.log10(pressure_squared / reference)
@@ -845,12 +874,8 @@ def third_octave_sound_pressure_level(
     mspl: xarray.DataArray (time, freq_bins)
         Sound pressure level [dB re 1 uPa] indexed by time and third octave bands
     """
-
-    # Third octave bin frequencies
-    bandwidth = 2 ** (1 / 3)
-    half_bandwidth = 2 ** (1 / 6)
-
-    mspl = _band_sound_pressure_level(spsd, bandwidth, half_bandwidth, fmin, fmax)
+    octave = 3
+    mspl = _band_sound_pressure_level(spsd, octave, fmin, fmax)
     mspl.attrs = {
         "units": "dB re 1 uPa",
         "long_name": "Third Octave Sound Pressure Level",
@@ -881,11 +906,8 @@ def decidecade_sound_pressure_level(
         Sound pressure level [dB re 1 uPa] indexed by time and decidecade bands
     """
 
-    # Decidecade bin frequencies
-    bandwidth = 2 ** (1 / 10)
-    half_bandwidth = 2 ** (1 / 20)
-
-    mspl = _band_sound_pressure_level(spsd, bandwidth, half_bandwidth, fmin, fmax)
+    octave = 10
+    mspl = _band_sound_pressure_level(spsd, octave, fmin, fmax)
     mspl.attrs = {
         "units": "dB re 1 uPa",
         "long_name": "Decidecade Sound Pressure Level",

mhkit/acoustics/graphics.py

@@ -1,6 +1,6 @@
 """
-This submodule provides essential plotting functions for visualizing passive acoustics 
-data. The functions allow for customizable plotting of sound pressure spectral density 
+This submodule provides essential plotting functions for visualizing passive acoustics
+data. The functions allow for customizable plotting of sound pressure spectral density
 levels across time and frequency dimensions.
 
 Each plotting function leverages the flexibility of Matplotlib, allowing for passthrough
@@ -11,12 +11,12 @@
 -------------
 1. **plot_spectrogram**:
 
-   - Generates a spectrogram plot from sound pressure spectral density level data, 
+   - Generates a spectrogram plot from sound pressure spectral density level data,
      with a logarithmic frequency scale by default for improved readability of acoustic data.
 
 2. **plot_spectra**:
 
-   - Produces a spectral density plot with a log-transformed x-axis, allowing for clear 
+   - Produces a spectral density plot with a log-transformed x-axis, allowing for clear
      visualization of spectral density across frequency bands.
 """
 

mhkit/acoustics/io.py

@@ -1,7 +1,7 @@
 """
 This submodule provides input/output functions for passive acoustics data,
 focusing on hydrophone recordings stored in WAV files. The main functionality
-includes reading and processing hydrophone data from various manufacturers 
+includes reading and processing hydrophone data from various manufacturers
 and exporting audio files for easy playback and analysis.
 
 Supported Hydrophone Models
@@ -14,28 +14,28 @@
 
 1. **Data Reading**:
 
-   - `read_hydrophone`: Main function to read a WAV file from a hydrophone and 
-     convert it to either a voltage or pressure time series, depending on the 
+   - `read_hydrophone`: Main function to read a WAV file from a hydrophone and
+     convert it to either a voltage or pressure time series, depending on the
      availability of sensitivity data.
 
-   - `read_soundtrap`: Wrapper for reading Ocean Instruments SoundTrap hydrophone 
+   - `read_soundtrap`: Wrapper for reading Ocean Instruments SoundTrap hydrophone
      files, automatically using appropriate metadata.
 
-   - `read_iclisten`: Wrapper for reading Ocean Sonics icListen hydrophone files, 
-     including metadata processing to apply hydrophone sensitivity for direct 
+   - `read_iclisten`: Wrapper for reading Ocean Sonics icListen hydrophone files,
+     including metadata processing to apply hydrophone sensitivity for direct
      sound pressure calculation.
 
 2. **Audio Export**:
 
-   - `export_audio`: Converts processed sound pressure data back into a WAV file 
+   - `export_audio`: Converts processed sound pressure data back into a WAV file
      format, with optional gain adjustment to improve playback quality.
 
 3. **Data Extraction**:
 
-   - `_read_wav_metadata`: Extracts metadata from a WAV file, including bit depth 
+   - `_read_wav_metadata`: Extracts metadata from a WAV file, including bit depth
      and other header information.
 
-   - `_calculate_voltage_and_time`: Converts raw WAV data into voltage values and 
+   - `_calculate_voltage_and_time`: Converts raw WAV data into voltage values and
      generates a time index based on the sampling frequency.
 """
 

mhkit/dolfyn/adv/clean.py

@@ -1,5 +1,4 @@
-"""Module containing functions to clean data
-"""
+"""Module containing functions to clean data"""
 
 import warnings
 import numpy as np

mhkit/loads/__init__.py

@@ -1,7 +1,7 @@
 """
 The `loads` package of the MHKiT (Marine and Hydrokinetic Toolkit) library
 provides tools and functionalities for analyzing and visualizing loads data
-from marine and hydrokinetic (MHK) devices. This package is designed to 
+from marine and hydrokinetic (MHK) devices. This package is designed to
 assist engineers, researchers, and analysts in understanding the forces and
 stresses applied to MHK devices under various operational and environmental
 conditions.

mhkit/loads/extreme/__init__.py

@@ -3,7 +3,7 @@
 and wave data statistics.
 
 It includes methods for calculating peaks over threshold, estimating
-short-term extreme distributions,and performing wave amplitude 
+short-term extreme distributions,and performing wave amplitude
 normalization for most likely extreme response analysis.
 """
 

mhkit/loads/extreme/extremes.py

@@ -1,29 +1,29 @@
 """
 This module provides functionality for estimating the short-term and
-long-term extreme distributions of responses in a time series. It 
-includes methods for analyzing peaks, block maxima, and applying 
-statistical distributions to model extreme events. The module supports 
-various methods for short-term extreme estimation, including peaks 
-fitting with Weibull, tail fitting, peaks over threshold, and block 
-maxima methods with GEV (Generalized Extreme Value) and Gumbel 
-distributions. Additionally, it offers functionality to approximate 
-the long-term extreme distribution by weighting short-term extremes 
+long-term extreme distributions of responses in a time series. It
+includes methods for analyzing peaks, block maxima, and applying
+statistical distributions to model extreme events. The module supports
+various methods for short-term extreme estimation, including peaks
+fitting with Weibull, tail fitting, peaks over threshold, and block
+maxima methods with GEV (Generalized Extreme Value) and Gumbel
+distributions. Additionally, it offers functionality to approximate
+the long-term extreme distribution by weighting short-term extremes
 across different sea states.
 
 Functions:
-- ste_peaks: Estimates the short-term extreme distribution from peaks 
+- ste_peaks: Estimates the short-term extreme distribution from peaks
     distribution using specified statistical methods.
 - block_maxima: Finds the block maxima in a time-series data to be used
     in block maxima methods.
-- ste_block_maxima_gev: Approximates the short-term extreme distribution 
+- ste_block_maxima_gev: Approximates the short-term extreme distribution
     using the block maxima method with the GEV distribution.
-- ste_block_maxima_gumbel: Approximates the short-term extreme 
+- ste_block_maxima_gumbel: Approximates the short-term extreme
     distribution using the block maxima method with the Gumbel distribution.
-- ste: Alias for `short_term_extreme`, facilitating easier access to the 
+- ste: Alias for `short_term_extreme`, facilitating easier access to the
     primary functionality of estimating short-term extremes.
-- short_term_extreme: Core function to approximate the short-term extreme 
+- short_term_extreme: Core function to approximate the short-term extreme
     distribution from a time series using chosen methods.
-- full_seastate_long_term_extreme: Combines short-term extreme 
+- full_seastate_long_term_extreme: Combines short-term extreme
     distributions using weights to estimate the long-term extreme distribution.
 """
 

mhkit/loads/extreme/mler.py

@@ -1,16 +1,16 @@
 """
-This module provides functionalities to calculate and analyze Most 
+This module provides functionalities to calculate and analyze Most
 Likely Extreme Response (MLER) coefficients for wave energy converter
 design and risk assessment. It includes functions to:
 
   - Calculate MLER coefficients (`mler_coefficients`) from a sea state
     spectrum and a response Amplitude Response Operator (ARO).
   - Define and manipulate simulation parameters (`mler_simulation`) used
     across various MLER analyses.
-  - Renormalize the incoming amplitude of the MLER wave 
+  - Renormalize the incoming amplitude of the MLER wave
     (`mler_wave_amp_normalize`) to match the desired peak height for more
     accurate modeling and analysis.
-  - Export the wave amplitude time series (`mler_export_time_series`) 
+  - Export the wave amplitude time series (`mler_export_time_series`)
     based on the calculated MLER coefficients for further analysis or
     visualization.
 """