[CC] avoid empty list in _spectrum_smoothing

Author: Martin Glesser

mosqito/sq_metrics/tonality/tone_to_noise_ecma/_spectrum_smoothing.py

@@ -37,18 +37,17 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
 
     """
     n = spec.shape[0]
-    if len(spec.shape)>1:
+    if len(spec.shape) > 1:
         m = spec.shape[1]
     else:
         m = spec.shape[0]
         n = 1
     spec = spec.ravel()
-    stop = np.arange(1,n+1)*m
+    stop = np.arange(1, n + 1) * m
     freqs_in = freqs_in.ravel()
 
     # n-th octave bands filter
-    filter_freqs = _getFrequencies(
-        low_freq, high_freq, noct, G=10, fr=1000)["f"]
+    filter_freqs = _getFrequencies(low_freq, high_freq, noct, G=10, fr=1000)["f"]
     filter_freqs[len(filter_freqs) - 1, 2] = high_freq
     filter_freqs[0, 0] = low_freq
 
@@ -63,7 +62,6 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
             (freqs_in >= filter_freqs[i, 0]) & (freqs_in <= filter_freqs[i, 2])
         )[0]
 
-
         # If the frequency bin is empty, it is deleted from the list
         if len(bin_index) == 0:
             smoothed_spectrum = np.delete(smoothed_spectrum, i, axis=1)
@@ -74,27 +72,43 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
             # The spectral components within the frequency bin are averaged on an energy basis
             spec_sum = np.zeros((n))
             for j in range(n):
-                for k in bin_index[(bin_index<stop[j]) & (bin_index>(stop[j]-m))]:
-                    spec_sum[j] += 10 ** (spec[k] / 10)
-            smoothed_spectrum[:,i] = 10 * np.log10(spec_sum / len(bin_index[(bin_index<stop[j]) & (bin_index>(stop[j]-m))]))
+                if (
+                    len(bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - m))])
+                    != 0
+                ):
+                    spec_sum[j] = np.mean(
+                        10
+                        ** (
+                            spec[
+                                bin_index[
+                                    (bin_index < stop[j]) & (bin_index > (stop[j] - m))
+                                ]
+                            ]
+                            / 10
+                        )
+                    )
+                else:
+                    spec_sum[j] = 1e-12
+            smoothed_spectrum[:, i] = 10 * np.log10(
+                spec_sum
+                # / len(bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - m))])
+            )
         nb_bands -= 1
         i += 1
     # Pose of the smoothed spectrum on the frequency-axis
     low = np.zeros((n, filter_freqs.shape[0]))
     high = np.zeros((n, filter_freqs.shape[0]))
-    
+
     # Index of the lower and higher limit of each frequency bin into the original spectrum
     for i in range(len(filter_freqs)):
-        low[:,i] = np.argmin(np.abs(freqs_out - filter_freqs[i, 0]))
-        high[:,i] = np.argmin(np.abs(freqs_out - filter_freqs[i, 2]))
+        low[:, i] = np.argmin(np.abs(freqs_out - filter_freqs[i, 0]))
+        high[:, i] = np.argmin(np.abs(freqs_out - filter_freqs[i, 2]))
     low = low.astype(int)
     high = high.astype(int)
-        
-    
-    smooth_spec = np.zeros((n,m))
+
+    smooth_spec = np.zeros((n, m))
     for i in range(n):
         for j in range(filter_freqs.shape[0]):
-            smooth_spec[i,low[i,j]: high[i,j]] = smoothed_spectrum[i,j] 
-
+            smooth_spec[i, low[i, j] : high[i, j]] = smoothed_spectrum[i, j]
 
     return smooth_spec

tests/sq_metrics/tonality/test_tnr_ecma_tv.py

@@ -3,10 +3,9 @@
 try:
     import pytest
 except ImportError:
-    raise RuntimeError(
-        "In order to perform the tests you need the 'pytest' package."
-    )
+    raise RuntimeError("In order to perform the tests you need the 'pytest' package.")
 
+import numpy as np
 
 # Local application imports
 from mosqito.utils import load
@@ -38,6 +37,9 @@ def test_tnr_ecma_tv():
 
     # Compute tone-to-noise ratio
     t_tnr, tnr, prom, freq, time = tnr_ecma_tv(audio, fs, prominence=True)
+    np.testing.assert_almost_equal(max(t_tnr), 34.995108238375025)
+    assert np.count_nonzero(prom == True) == 6
+
 
 if __name__ == "__main__":
     test_tnr_ecma_tv()