Acoustics SPL bugfix (#379)

Fixes a bug where calculating SPLs in lower frequency bands was failing because the bandwidth resolution was on par with the frequency resolution, i.e. numpy.trapz was failing to calculate the integral because the input was only 1D over frequency. Also expanded testing to cover this particular case.

Author: jmcvey3

mhkit/acoustics/analysis.py

@@ -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/tests/acoustics/test_analysis.py

@@ -57,7 +57,7 @@ def test_averaging(self):
 
         # Frequency average into # octave bands
         octave = 3
-        td_spsdl_mean = acoustics.band_aggregate(td_spsdl, octave, fmin=50)
+        td_spsdl_mean = acoustics.band_aggregate(td_spsdl, octave, fmin=10, fmax=100000)
 
         # Time average into 30 s bins
         lbin = 30
@@ -81,29 +81,29 @@ def test_averaging(self):
         )
         cd_spsdl_50 = np.array(
             [
-                [73.71803613, 70.97557445, 69.79906778, 69.04934313, 67.56449352],
-                [73.72245955, 71.53327285, 70.55206775, 68.69638127, 67.75243522],
-                [73.64022645, 72.24548986, 70.09995522, 69.00394292, 68.22919418],
-                [73.1301846, 71.99940268, 70.56372046, 69.01366589, 67.19515351],
-                [74.67880072, 71.27235403, 70.23024477, 67.4915765, 66.73024553],
+                [63.45507, 64.753525, 65.04905, 67.15576, 73.47938],
+                [62.77437, 64.58199, 65.18464, 66.37395, 72.30796],
+                [64.76277, 64.950264, 65.80557, 67.88482, 73.24013],
+                [63.654488, 62.31394, 65.598816, 67.370674, 71.52472],
+                [62.45623, 62.461388, 62.111694, 66.06419, 72.324936],
             ]
         )
         cd_spsdl_25 = np.array(
             [
-                [72.42136105, 70.37422873, 68.60783404, 67.56108417, 66.4751517],
-                [71.95173902, 71.03281659, 69.59019407, 67.79615712, 66.73980611],
-                [71.12756436, 70.68228634, 69.53891917, 68.126758, 67.48463198],
-                [71.71909635, 70.1849931, 69.22647784, 68.14102709, 66.18740693],
-                [72.25521793, 70.18087912, 68.97354823, 66.71295946, 65.35302077],
+                [59.33189297, 62.89503765, 61.60455799, 64.80938911, 70.59576607],
+                [60.37440872, 60.69928551, 61.9694643, 64.91986465, 70.00148964],
+                [61.1297617, 63.02504444, 64.41207123, 66.37802315, 71.38513947],
+                [59.52737236, 59.45869541, 62.48176765, 66.0959053, 70.06054497],
+                [58.55439758, 59.88098335, 59.66310596, 63.86431885, 70.20335197],
             ]
         )
         cd_spsdl_75 = np.array(
             [
-                [75.29614796, 71.86901413, 71.08418954, 69.6835928, 68.26993291],
-                [74.51608597, 72.82376854, 71.31219865, 70.38580566, 69.01731822],
-                [75.17013043, 73.45962974, 71.30593827, 71.50687178, 69.49805535],
-                [74.38176106, 73.13456376, 72.13861655, 70.45825381, 67.93458589],
-                [75.52387419, 72.99604074, 71.26831962, 68.90629303, 67.79114848],
+                [66.33672714, 67.13593102, 67.34234238, 68.7525959, 75.30982399],
+                [64.58539009, 66.84792709, 67.11526108, 69.7322197, 74.50746346],
+                [66.56425095, 67.85562325, 69.30602646, 69.83069992, 74.79984283],
+                [67.34252357, 65.65701294, 67.48604202, 70.948246, 73.59340286],
+                [66.26214409, 65.43437958, 64.36196518, 67.67719078, 74.33639717],
             ]
         )
 
@@ -118,13 +118,17 @@ def test_freq_loss(self):
         self.assertEqual(fmin, 39.84375)
 
     def test_spl(self):
-        td_spl = acoustics.sound_pressure_level(self.spsd, fmin=50)
+        td_spl = acoustics.sound_pressure_level(self.spsd, fmin=10, fmax=100000)
 
         # Decidecade octave sound pressure level
-        td_spl10 = acoustics.decidecade_sound_pressure_level(self.spsd, fmin=50)
+        td_spl10 = acoustics.decidecade_sound_pressure_level(
+            self.spsd, fmin=10, fmax=100000
+        )
 
         # Median third octave sound pressure level
-        td_spl3 = acoustics.third_octave_sound_pressure_level(self.spsd, fmin=50)
+        td_spl3 = acoustics.third_octave_sound_pressure_level(
+            self.spsd, fmin=10, fmax=100000
+        )
 
         cc = np.array(
             [
@@ -136,31 +140,74 @@ def test_spl(self):
             ],
             dtype="datetime64[ns]",
         )
-        cd_spl = np.array(
-            [97.48727775, 98.21888437, 96.99586637, 97.43571891, 96.60915502]
+        cd_spl_head = np.array([98.12284, 98.639824, 97.62718, 97.85709, 96.98539])
+        cd_spl_tail = np.array([98.420975, 98.10879, 97.430115, 97.99395, 97.95798])
+
+        cd_spl10_freq_head = np.array(
+            [10.0, 10.717735, 11.486984, 12.311444, 13.195079]
+        )
+        cd_spl10_head = np.array(
+            [
+                [64.40389, 64.78221, 63.64469, 67.782845, 73.05421],
+                [56.934277, 62.80422, 66.329056, 67.336395, 65.79995],
+                [67.16593, 69.089096, 69.5835, 67.69844, 63.657196],
+                [66.010025, 67.567635, 68.16686, 66.72678, 64.52246],
+                [63.887203, 68.73698, 72.71424, 72.98322, 69.08172],
+            ]
         )
-        cd_spl10 = np.array(
+        cd_spl10_freq_tail = np.array(
+            [38217.031333, 40960.0, 43899.841025, 47050.684621, 50427.675171]
+        )
+        cd_spl10_tail = np.array(
             [
-                [82.06503071, 78.20349846, 79.78088446, 75.31281183, 82.1194826],
-                [82.66175023, 79.77804574, 82.86005403, 77.57078269, 76.7598224],
-                [77.48975416, 82.72580274, 83.88251531, 74.71242694, 74.01377947],
-                [79.11312683, 76.56114947, 82.18953494, 75.40888015, 74.80285354],
-                [81.26751434, 82.29074565, 80.08831394, 75.75364773, 73.52176641],
+                [77.38338, 73.43317, 72.7755, 72.53339, np.nan],
+                [77.72596, 73.57787, 73.16561, 72.637436, np.nan],
+                [77.61121, 73.59171, 73.20265, 72.57601, np.nan],
+                [77.3753, 73.35718, 72.89339, 72.386765, np.nan],
+                [77.31649, 73.806496, 73.296, 72.73348, np.nan],
             ]
         )
-        cd_spl3 = np.array(
+        cd_spl3_freq_head = np.array([10.0, 12.59921, 15.874011, 20.0, 25.198421])
+        cd_spl3_head = np.array(
             [
-                [86.5847236, 84.98068691, 85.61056131, 83.55067796, 84.41810962],
-                [87.5449842, 84.48841036, 84.09406069, 85.81895309, 86.71437852],
-                [86.37334939, 84.08914125, 86.01614536, 83.36059983, 84.54635288],
-                [84.21413445, 84.63996392, 82.52906024, 84.54731095, 83.45652422],
-                [86.90033232, 84.8217658, 83.85297355, 82.92231618, 81.39163217],
+                [68.88561, 75.65294, 68.29522, 75.80323, 82.53724],
+                [62.806908, 69.76993, 62.64113, 73.26091, 83.27883],
+                [71.73166, 68.541534, 68.056076, 75.438034, 84.268715],
+                [70.84345, 68.65471, 63.4681, 72.818085, 77.38771],
+                [69.23148, 74.04387, 64.49707, 74.146164, 79.52727],
+            ]
+        )
+        cd_spl3_freq_tail = np.array(
+            [20480.0, 25803.183102, 32509.973544, 40960.0, 51606.366204]
+        )
+        cd_spl3_tail = np.array(
+            [
+                [80.37833, 81.21788, 83.5725, 80.37073, 72.06452],
+                [81.848434, 81.772064, 83.928505, 80.70311, 72.164345],
+                [81.13474, 81.67803, 83.96902, 80.6636, 72.07929],
+                [81.532005, 81.694954, 83.796875, 80.38368, 71.94872],
+                [80.70353, 81.6905, 83.76083, 80.53248, 72.248276],
             ]
         )
 
-        np.testing.assert_allclose(td_spl.head().values, cd_spl, atol=1e-6)
-        np.testing.assert_allclose(td_spl10.head().values, cd_spl10, atol=1e-6)
-        np.testing.assert_allclose(td_spl3.head().values, cd_spl3, atol=1e-6)
+        np.testing.assert_allclose(td_spl.head().values, cd_spl_head, atol=1e-6)
+        np.testing.assert_allclose(td_spl.tail().values, cd_spl_tail, atol=1e-6)
+        np.testing.assert_allclose(
+            td_spl10["freq_bins"].head().values, cd_spl10_freq_head, atol=1e-6
+        )
+        np.testing.assert_allclose(td_spl10.head().values, cd_spl10_head, atol=1e-6)
+        np.testing.assert_allclose(
+            td_spl10["freq_bins"].tail().values, cd_spl10_freq_tail, atol=1e-6
+        )
+        np.testing.assert_allclose(td_spl10.tail().values, cd_spl10_tail, atol=1e-6)
+        np.testing.assert_allclose(
+            td_spl3["freq_bins"].head().values, cd_spl3_freq_head, atol=1e-6
+        )
+        np.testing.assert_allclose(td_spl3.head().values, cd_spl3_head, atol=1e-6)
+        np.testing.assert_allclose(
+            td_spl3["freq_bins"].tail().values, cd_spl3_freq_tail, atol=1e-6
+        )
+        np.testing.assert_allclose(td_spl3.tail().values, cd_spl3_tail, atol=1e-6)
         np.testing.assert_equal(td_spl["time"].head().values, cc)
 
     def test_sound_pressure_spectral_density(self):