Merge pull request #72 from Eomys/minor

Minor corrections

Author: wantysal

README.md

@@ -42,18 +42,6 @@ MOSQITO is available on [pip](https://pypi.org/project/pip/). Simply type in a s
 
 This command line should download and install MOSQITO on your computer, along with the dependencies needed to compute SQ metrics.
 
-If you need to import .uff or .unv files, you will need the pyuff package dependency. Note that 'pyuff' is released under the GPL license which prevents MOSQITO from being used in other software that must be under a more permissive license. To include the 'pyuff' dependancy anyway, type the following command:
-
-    pip install mosqito[uff]
-
-If you want to use MOSQITO coupled with SciDataTool, you will need SDT package dependency. To install it along with MOSQITO, use the following command:
-
-    pip install mosqito[SciDataTool]
-
-Note that all the depencies can be installed at once using:
-
-    pip install mosqito[all]
-
 ## Contact
 
 You can contact us on Github by opening an issue (to request a feature, ask a question or report a bug). 

full_requirements.txt

@@ -2,5 +2,6 @@ numpy
 scipy
 matplotlib
 pyuff
-pytest
-scidatatool
+pytest>=5.4.1
+pandas
+openpyxl

mosqito/sq_metrics/loudness/loudness_zwst/loudness_zwst.py

@@ -9,16 +9,8 @@
 from mosqito.sq_metrics.loudness.loudness_zwst._calc_slopes import _calc_slopes
 from mosqito.utils.conversion import amp2db
 
-# Optional package import
-try:
-    from SciDataTool import DataTime, DataLinspace, DataFreq
-except ImportError:
-    DataTime = None
-    DataLinspace = None
-    DataFreq = None
 
-
-def loudness_zwst(signal, fs=None, field_type="free", is_sdt_output=False):
+def loudness_zwst(signal, fs, field_type="free"):
     """Zwicker-loudness calculation for stationary signals
 
     Calculates the acoustic loudness according to Zwicker method for
@@ -36,33 +28,32 @@ def loudness_zwst(signal, fs=None, field_type="free", is_sdt_output=False):
 
     Parameters
     ----------
-    signal : numpy.array or DataTime object
+    signal : numpy.array
         Signal time values [Pa]
-    fs : float, optional
-        Sampling frequency, can be omitted if the input is a DataTime
-        object. Default to None
+    fs : float
+        Sampling frequency [Hz]
     field_type : str
         Type of soundfield corresponding to spec_third ("free" by
         default or "diffuse").
-    is_sdt_output : Bool, optional
-        If True, the outputs are returned as SciDataTool objects.
-        Default to False
 
     Outputs
     -------
     N : float or numpy.array
         The overall loudness array [sones], size (Ntime,).
-    N_specific : numpy.ndarray or DataFreq object
+    N_specific : numpy.ndarray
         The specific loudness array [sones/bark], size (Nbark, Ntime).
     bark_axis: numpy.array
         The Bark axis array, size (Nbark,).
     """
 
-    # Manage SciDataTool input type
-    if DataTime is not None and isinstance(signal, DataTime):
-        time = signal.get_along("time")["time"]
-        fs = 1 / (time[1] - time[0])
-        signal = signal.get_along("time")[signal.symbol]
+    if fs < 48000:
+        print(
+            "[Warning] Signal resampled to 48 kHz to allow calculation. To fulfill the standard requirements fs should be >=48 kHz."
+        )
+        from scipy.signal import resample
+
+        signal = resample(signal, int(48000 * len(signal) / fs))
+        fs = 48000
 
     # Compute third octave band spectrum
     spec_third, _ = noct_spectrum(signal, fs, fmin=24, fmax=12600)
@@ -80,27 +71,4 @@ def loudness_zwst(signal, fs=None, field_type="free", is_sdt_output=False):
     # Define Bark axis
     bark_axis = np.linspace(0.1, 24, int(24 / 0.1))
 
-    # Manage SciDataTool output type
-    if is_sdt_output:
-        if DataLinspace is None:
-            raise RuntimeError(
-                "In order to handle Data objects you need the 'SciDataTool' package."
-            )
-        else:
-            bark_data = DataLinspace(
-                name="Critical band rate",
-                unit="Bark",
-                initial=0,
-                final=24,
-                number=int(24 / 0.1),
-                include_endpoint=True,
-            )
-            N_specific = DataFreq(
-                name="Specific loudness (Zwicker method for stationnary signal)",
-                symbol="N'_{zwst}",
-                axes=[bark_data],
-                values=N_specific,
-                unit="sone/Bark",
-            )
-
     return N, N_specific, bark_axis

mosqito/sq_metrics/loudness/loudness_zwst/loudness_zwst_freq.py

@@ -1,7 +1,8 @@
 # -*- coding: utf-8 -*-
 
 # Third party imports
-import numpy as np
+from numpy import linspace, arange, tile, diff
+from scipy.interpolate import interp1d
 
 # Local application imports
 from mosqito.sound_level_meter.noct_spectrum.noct_synthesis import noct_synthesis
@@ -45,8 +46,45 @@ def loudness_zwst_freq(spectrum, freqs, field_type="free"):
     bark_axis : numpy.array
         Frequency axis in bark, size (Nbark,).
     """
-    if len(spectrum) != len(freqs):
-        raise ValueError('Input spectrum and frequency axis must have the same shape')
+
+    # 1D spectrum
+    if len(spectrum.shape) == 1:
+        if len(spectrum) != len(freqs):
+            raise ValueError(
+                'Input spectrum and frequency axis do not have the same length')
+        if (freqs.max() < 24000) or (freqs.min() > 24):
+            print("[WARNING] freqs argument is not wide enough to cover the full audio range. Missing frequency bands will be filled with 0. To fulfill the standard requirements, the frequency axis should go from 24Hz up to 24 kHz."
+            )
+            df = freqs[1] - freqs[0]
+            spectrum = interp1d(freqs, spectrum, axis=0, bounds_error=False, fill_value=0)(linspace(0, 24000, int(24000//df)))
+            freqs = linspace(0, 24000, int(24000//df))
+    # 2D spectrum
+    elif len(spectrum.shape) > 1:
+        nseg = spectrum.shape[1]
+        # one frequency axis per segment
+        if len(freqs.shape) > 1:
+            if spectrum.shape != freqs.shape:
+                raise ValueError(
+                    'Input spectrum and frequency axis do not have the same shape.')
+            if (freqs.max() < 24000) or (freqs.min() > 24):
+                print("[WARNING] freqs argument is not wide enough to cover the full audio range. Missing frequency bands will be filled with 0. To fulfill the standard requirements, the frequency axis should go from 24Hz up to 24 kHz."
+                )
+                df = diff(freqs, axis=0).min()
+                nperseg = int(24000//df)
+                spectrum = interp1d(freqs.ravel(), spectrum.ravel(), bounds_error=False, fill_value=0)(tile(linspace(0, 24000, nperseg), nseg)).reshape((nseg, nperseg)).T
+                freqs = tile(linspace(0, 24000, nperseg), (nseg, 1)).T   
+        # one frequency axis for all the segments
+        elif len(freqs.shape) == 1:
+            if spectrum.shape[0] != len(freqs):
+                raise ValueError(
+                    'Input spectra and frequency axis do not have the same length. "freqs" must have dim(nperseg) and "spectra" must have dim(nperseg,nseg).')
+            if (freqs.max() < 24000) or (freqs.min() > 24):
+                print("[WARNING] freqs argument is not wide enough to cover the full audio range. Missing frequency bands will be filled with 0. To fulfill the standard requirements, the frequency axis should go from 24Hz up to 24 kHz."
+                )
+                df = freqs[1] - freqs[0]
+                nperseg = int(24000//df)
+                spectrum = interp1d(freqs, spectrum, axis=0, bounds_error=False, fill_value=0)(linspace(0, 24000, nperseg))
+                freqs = tile(linspace(0, 24000, nperseg), (nseg, 1)).T
 
     # Compute third octave band spectrum
     spec_third, _ = noct_synthesis(spectrum, freqs, fmin=24, fmax=12600)
@@ -62,6 +100,6 @@ def loudness_zwst_freq(spectrum, freqs, field_type="free"):
     N, N_specific = _calc_slopes(Nm)
 
     # Define Bark axis
-    bark_axis = np.linspace(0.1, 24, int(24 / 0.1))
+    bark_axis = linspace(0.1, 24, int(24 / 0.1))
 
     return N, N_specific, bark_axis

mosqito/sq_metrics/loudness/loudness_zwst/loudness_zwst_perseg.py

@@ -4,18 +4,8 @@
 from mosqito.utils import time_segmentation
 from mosqito.sq_metrics.loudness.loudness_zwst.loudness_zwst import loudness_zwst
 
-# Optional package import
-try:
-    from SciDataTool import DataTime, DataLinspace, DataFreq, Norm_func
-except ImportError:
-    DataTime = None
-    DataLinspace = None
-    DataFreq = None
 
-
-def loudness_zwst_perseg(
-    signal, fs=None, nperseg=4096, noverlap=None, field_type="free", is_sdt_output=False
-):
+def loudness_zwst_perseg(signal, fs, nperseg=4096, noverlap=None, field_type="free"):
     """Zwicker-loudness calculation for stationary signals
 
     Calculates the acoustic loudness according to Zwicker method for
@@ -33,10 +23,10 @@ def loudness_zwst_perseg(
 
     Parameters
     ----------
-    signal : numpy.array or DataTime object
+    signal : numpy.array
         Time signal values [Pa].
-    fs : float, optional
-        Sampling frequency, can be omitted if the input is a DataTime
+    fs : float
+        Sampling frequency [Hz]
         object. Default to None
     nperseg: int, optional
         Length of each segment. Defaults to 4096.
@@ -46,74 +36,32 @@ def loudness_zwst_perseg(
     field_type : str
         Type of soundfield corresponding to spec_third ("free" by
         default or "diffuse").
-    is_sdt_output : Bool, optional
-        If True, the outputs are returned as SciDataTool objects.
-        Default to False
 
     Outputs
     -------
-    N : numpy.array or DataTime object
+    N : numpy.array
         The overall loudness array [sones], size (Ntime,).
-    N_specific : numpy.ndarray or DataFreq object
+    N_specific : numpy.ndarray
         The specific loudness array [sones/bark], size (Nbark, Ntime).
     bark_axis: numpy.array
         The Bark axis array, size (Nbark,).
     time_axis: numpy.array
-        The time axis array, size (Ntime,) or None.
+        The time axis array, size (Ntime,).
 
     """
+    if fs < 48000:
+        print(
+            "[Warning] Signal resampled to 48 kHz to allow calculation. To fulfill the standard requirements fs should be >=48 kHz."
+        )
+        from scipy.signal import resample
 
-    # Manage input type
-    if DataTime is not None and isinstance(signal, DataTime):
-        time = signal.get_along("time")["time"]
-        fs = 1 / (time[1] - time[0])
-        signal = signal.get_along("time")[signal.symbol]
+        signal = resample(signal, int(48000 * len(signal) / fs))
+        fs = 48000
 
     # Time signal segmentation
     signal, time_axis = time_segmentation(signal, fs, nperseg, noverlap)
 
     # Compute loudness
     N, N_specific, bark_axis = loudness_zwst(signal, fs, field_type=field_type)
 
-    # Manage SciDataTool output type
-    if is_sdt_output:
-        if DataLinspace is None:
-            raise RuntimeError(
-                "In order to handle Data objects you need the 'SciDataTool' package."
-            )
-        else:
-            bark_data = DataLinspace(
-                name="Critical band rate",
-                unit="Bark",
-                initial=0,
-                final=24,
-                number=int(24 / 0.1),
-                include_endpoint=True,
-                normalizations={
-                    "Hz": Norm_func(function=lambda x: 1960 * (x + 0.53) / (26.28 - x))
-                },
-            )
-            time = DataLinspace(
-                name="time",
-                unit="s",
-                initial=time_axis[0],
-                final=time_axis[-1],
-                number=len(time_axis),
-                include_endpoint=True,
-            )
-            N_specific = DataFreq(
-                name="Specific loudness (Zwicker method for stationnary signal)",
-                symbol="N'_{zwst}",
-                axes=[bark_data, time],
-                values=N_specific,
-                unit="sone/Bark",
-            )
-            N = DataTime(
-                name="Loudness (Zwicker method for stationnary signal)",
-                symbol="N_{zwst}",
-                axes=[time],
-                values=N,
-                unit="sone",
-            )
-
     return N, N_specific, bark_axis, time_axis

mosqito/sq_metrics/loudness/loudness_zwtv/loudness_zwtv.py

@@ -45,6 +45,14 @@ def loudness_zwtv(signal, fs, field_type='free'):
 
     """        
 
+    if fs < 48000:
+        print("[Warning] Signal resampled to 48 kHz to allow calculation. To fulfill the standard requirements fs should be >=48 kHz."
+             )
+        from scipy.signal import resample
+        signal = resample(signal, int(48000 * len(signal) / fs))
+        fs = 48000
+
+
     # Compute third octave band spectrum vs. time
     spec_third, time_axis, _ = _third_octave_levels(signal, fs)