Merge pull request #25 from mglesser/master

Documentation updates + bug fix in Loudness computation

Author: mglesser

CITATION.cff

@@ -0,0 +1,18 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+authors:
+  - family-names: Green Forge Coop
+    given-names:
+title: MOSQITO
+version: 0.3.3
+doi: 10.5281/zenodo.5419590
+date-released: 2021-09-03
+keywords:
+  - audio
+  - python
+  - signal-processing
+  - audio-analysis
+  - acoustic
+  - sound-quality
+  - sq-metrics
+license: Apache-2.0

README.md

@@ -35,11 +35,34 @@ It is written in Python, one of the most popular free programming language in th
 
 Tutorials are available in the [tutorials](./tutorials/) folder. Documentation and validation of the MOSQITO functions are available in the [documentation](./documentation/) folder.
 
+## Getting MOSQITO
+MOSQITO is available on [pip](https://pypi.org/project/pip/). Simply type in a shell the following command:
+
+    pip install mosqito
+
+This command line should download and install MOSQITO on your computer, along with all the needed dependencies.
+
+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]
+
 ## Contact
 
-You can contact us on Github by opening an issue (to request a feature, ask a question or report a bug).
+You can contact us on Github by opening an issue (to request a feature, ask a question or report a bug). 
+
+## Citing MOSQITO
+
+If you are using MOSQITO in your research activities, please help our scientific visibility by citing our work! You can use the following citation in APA format:
+
+Green Forge Coop. MOSQITO [Computer software]. https://doi.org/10.5281/zenodo.5284054
+
+If you need to cite the current release of MOSQITO, please use the "Cite this repository" feature in the "About" section of this Github repository.
+
+
+## Publications citing MOSQITO
+
+Glesser, M., Ni, S., Degrendele, K., Wanty, S., & Le Besnerais, J. (2021, October 13-14). *Sound quality analysis of Electric Drive Units under different switching control strategies* [Paper presentation]. Automotive NVH Comfort Congress 2021, Le Mans , France.
+
+San Millán-Castillo, R., Latorre-Iglesias, E., Glesser, M., Wanty, S., Jiménez-Caminero, D., & Álvarez-Jimeno, J.M. (2021). MOSQITO: an open-source and free toolbox for sound quality metrics in the industry and education. *INTER-NOISE and NOISE-CON Congress and Conference Proceedings*, 12, 1164-1175. https://doi.org/10.3397/IN-2021-1767
 
-## How to cite MOSQITO
 
-If you use MOSQITO for your research activities and need to cite the software in a publication, please use the following citation:
-TODO

mosqito/functions/loudness_zwicker/loudness_zwicker_nonlinear_decay.py

@@ -26,7 +26,7 @@ def calc_nl_loudness(core_loudness):
     """
     # Initialization
     sample_rate = 2000
-    nl_loudness = core_loudness
+    nl_loudness = core_loudness.copy()
     # Factor for virtual upsampling/inner iterations
     nl_iter = 24
     # Time constants for non_linear temporal decay
@@ -52,17 +52,17 @@ def calc_nl_loudness(core_loudness):
     ]
     nl_lp = {"B": B}
 
-    for i_cl in np.arange(np.shape(core_loudness)[0]):
+    for i_cl in np.arange(np.shape(nl_loudness)[0]):
         # At beginning capacitors C1 and C2 are discharged
         nl_lp["uo_last"] = 0
         nl_lp["u2_last"] = 0
 
-        for i_time in np.arange(np.shape(core_loudness)[1] - 1):
+        for i_time in np.arange(np.shape(nl_loudness)[1] - 1):
             # interpolation steps between current and next sample
             delta = (
-                core_loudness[i_cl, i_time + 1] - core_loudness[i_cl, i_time]
+                nl_loudness[i_cl, i_time + 1] - nl_loudness[i_cl, i_time]
             ) / nl_iter
-            ui = core_loudness[i_cl, i_time]
+            ui = nl_loudness[i_cl, i_time]
             nl_lp = calc_nl_lp(ui, nl_lp)
             nl_loudness[i_cl, i_time] = nl_lp["uo_last"]
             ui += delta
@@ -71,9 +71,9 @@ def calc_nl_loudness(core_loudness):
             for i_in in np.arange(1, nl_iter):
                 nl_lp = calc_nl_lp(ui, nl_lp)
                 ui += delta
-        nl_lp = calc_nl_lp(core_loudness[i_cl, i_time + 1], nl_lp)
+        nl_lp = calc_nl_lp(nl_loudness[i_cl, i_time + 1], nl_lp)
         nl_loudness[i_cl, i_time + 1] = nl_lp["uo_last"]
-    return core_loudness
+    return nl_loudness
 
 
 def calc_nl_lp(ui, nl_lp):
@@ -104,7 +104,7 @@ def calc_nl_lp(ui, nl_lp):
                 uo = ui  # lower than ui
             u2 = uo
     else:
-        if abs(ui - nl_lp["uo_last"] < 1e-5):  # case 2 (charge)
+        if abs(ui - nl_lp["uo_last"]) < 1e-5:  # case 2 (charge)
             uo = ui
             if uo > nl_lp["u2_last"]:  # case 2.1
                 u2 = (nl_lp["u2_last"] - ui) * nl_lp["B"][5] + ui

mosqito/functions/loudness_zwicker/loudness_zwicker_shared.py

@@ -169,6 +169,219 @@ def calc_main_loudness(spec_third, field_type):
     return nm
 
 
+def calc_main_loudness_ea(spec_third, field_type):
+    """Calculate core loudness
+
+    The code is based on BASIC program published in "Program for
+    calculating loudness according to DIN 45631 (ISO 532B)", E.Zwicker
+    and H.Fastl, J.A.S.J (E) 12, 1 (1991).
+    It also corresponds to the following functions of the C program
+    published with ISO 532-1:2017:
+    - corr_third_octave_intensities
+    - f_calc_lcbs
+    - f_calc_core_loudness
+    - f_corr_loudness
+
+    Parameters
+    ----------
+    spec_third : numpy.ndarray
+        A third octave band spectrum [dB ref. 2e-5 Pa]
+    field_type : str
+        Type of soundfield correspondin to spec_third ("free" or
+        "diffuse")
+
+    Outputs
+    -------
+    nm :  numpy.ndarray
+        Core loudness
+    """
+    #
+    # Date tables definition (variable names and description according to
+    # Zwicker:1991)
+    # Ranges of 1/3 octave band levels for correction at low frequencies
+    # according to equal loudness contours
+    rap = np.array([45, 55, 65, 71, 80, 90, 100, 120])
+    # Reduction of 1/3 octave band levels at low frequencies according to
+    # equal loudness contours within the eight ranges defined by RAP
+    dll = np.array(
+        [
+            (-32, -24, -16, -10, -5, 0, -7, -3, 0, -2, 0),
+            (-29, -22, -15, -10, -4, 0, -7, -2, 0, -2, 0),
+            (-27, -19, -14, -9, -4, 0, -6, -2, 0, -2, 0),
+            (-25, -17, -12, -9, -3, 0, -5, -2, 0, -2, 0),
+            (-23, -16, -11, -7, -3, 0, -4, -1, 0, -1, 0),
+            (-20, -14, -10, -6, -3, 0, -4, -1, 0, -1, 0),
+            (-18, -12, -9, -6, -2, 0, -3, -1, 0, -1, 0),
+            (-15, -10, -8, -4, -2, 0, -3, -1, 0, -1, 0),
+        ]
+    )
+    # Critical band level at absolute threshold without taking into
+    # account the transmission characteristics of the ear
+    ltq = np.array([30, 18, 12, 8, 7, 6, 5, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3])
+    # Correction of levels according to the transmission characteristics
+    # of the ear
+    a0 = np.array(
+        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -1.6, -3.2, -5.4, -5.6, -4, -1.5, 2, 5, 12]
+    )
+    # Level difference between free and diffuse sound fields
+    ddf = np.array(
+        [
+            0,
+            0,
+            0.5,
+            0.9,
+            1.2,
+            1.6,
+            2.3,
+            2.8,
+            3,
+            2,
+            0,
+            -1.4,
+            -2,
+            -1.9,
+            -1,
+            0.5,
+            3,
+            4,
+            4.3,
+            4,
+        ]
+    )
+    # Adaptation of 1/3 oct. band levels to the corresponding critical
+    # band level
+    dcb = np.array(
+        [
+            -0.25,
+            -0.6,
+            -0.8,
+            -0.8,
+            -0.5,
+            0,
+            0.5,
+            1.1,
+            1.5,
+            1.7,
+            1.8,
+            1.8,
+            1.7,
+            1.6,
+            1.4,
+            1.2,
+            0.8,
+            0.5,
+            0,
+            -0.5,
+        ]
+    )
+    #
+    # Correction of 1/3 oct. band levels according to equal loudness
+    # contours 'xp' and calculation of the intensities for 1/3 oct.
+    # bands up to 315 Hz
+
+    # Prepare al arrays to work with
+    spec_third_aux = spec_third[: dll.shape[1], :]
+    spec_third_aux[:, -1] = 0
+
+    # Convert rap, dll in 3 dimensional array
+    # 1. generate the array shape
+    base_mat = np.ones(dll.shape[0] * dll.shape[1] * spec_third_aux.shape[1]).reshape(
+        dll.shape[0], dll.shape[1], spec_third_aux.shape[1]
+    )
+    # 2. start saving data in rap and DLL in a 3D array
+    rap_mat = np.array(
+        [base_mat[:, i, :].T * rap for i in np.arange(dll.shape[1])]
+    ).transpose(2, 0, 1)
+    dll_mat = np.array(
+        [np.multiply(base_mat[:, :, i], dll) for i in np.arange(spec_third.shape[1])]
+    ).transpose(1, 2, 0)
+    spec_third_aux_mat = np.array(
+        [
+            np.multiply(base_mat[i, :, :], spec_third_aux)
+            for i in np.arange(dll.shape[0])
+        ]
+    )
+    # create the the array rap-dll
+    rap_dll_mat = rap_mat - dll_mat
+
+    # This part substitutes the while loop.
+    # create the mask to operate
+    logic_mat = spec_third_aux_mat > rap_dll_mat
+    dll_result = dll_mat[0, :, :]
+    dll_result[logic_mat[0, :, :]] = 0
+    for i in np.arange(1, dll_mat.shape[0] - 1):
+        mask = np.logical_xor(logic_mat[i - 1, :, :], logic_mat[i, :, :])
+        dll_result[mask] = dll_mat[i, mask]
+
+    xp = dll_result + spec_third_aux
+    ti = np.power(10, (xp / 10))
+
+    # Determination of levels LCB(1), LCB(2) and LCB(3) within the
+    # first three critical bands
+
+    gi = np.zeros([3, dll_result.shape[1]])
+    gi[0, :] = ti[0:6, :].sum(axis=0)
+    gi[1, :] = ti[6:9, :].sum(axis=0)
+    gi[2, :] = ti[9:11, :].sum(axis=0)
+
+    logic_gi = gi > 0
+    lcb = np.zeros([3, dll_result.shape[1]])
+    lcb = 10 * np.log10(gi[logic_gi])
+    lcb = lcb.reshape(3, dll_result.shape[1])
+
+    # Calculation of main loudness
+    s = 0.25
+    nm = np.zeros([20, spec_third_aux.shape[1]])
+    le = np.copy(spec_third[8:, :])
+    # le = le.reshape((20))
+    le[0:3, :] = lcb
+    a0_mat = np.ones(a0.shape[0] * spec_third.shape[1]).reshape(
+        a0.shape[0], spec_third.shape[1]
+    )
+    a0_mat = (a0_mat.T * a0).T
+    le = le - a0_mat
+
+    if field_type == "diffuse":
+        ddf_mat = np.ones(ddf.shape[0] * spec_third.shape[1]).reshape(
+            ddf.shape[0], spec_third.shape[1]
+        )
+        ddf_mat = (ddf_mat.T * ddf).T
+        le += ddf_mat
+
+    ltq_mat = np.ones(ltq.shape[0] * spec_third.shape[1]).reshape(
+        ltq.shape[0], spec_third.shape[1]
+    )
+    ltq_mat = (ltq_mat.T * ltq).T
+    i = le > ltq_mat
+    dcb_mat = np.ones(dcb.shape[0] * spec_third.shape[1]).reshape(
+        dcb.shape[0], spec_third.shape[1]
+    )
+    dcb_mat = (dcb_mat.T * dcb).T
+    le[i] -= dcb_mat[i]
+
+    mp1 = 0.0635 * np.power(10, 0.025 * ltq_mat)
+    mp1[i == False] = 0
+    mp2 = np.power(1 - s + s * np.power(10, 0.1 * (le - ltq_mat)), 0.25) - 1
+    mp2[i == False] = 0
+
+    nm = np.multiply(mp1, mp2)
+
+    nm[i == False] = 0
+    nm[nm < 0] = 0
+    current_shape = nm.shape
+    nm = np.append(nm, np.zeros(current_shape[1])).reshape(
+        current_shape[0] + 1, current_shape[1]
+    )
+    #
+    # Correction of specific loudness in the lowest critical band
+    # taking into account the dependance of absolute threshold
+    # within this critical band
+    korry = 0.4 + 0.32 * nm[0] ** 0.2
+    nm[0, korry <= 1] *= korry
+    nm[:, -1] = 0
+    return nm
+
+
 def calc_slopes(nm):
     """"""
     # Upper limits of approximated critical bands in terms of critical

mosqito/functions/loudness_zwicker/loudness_zwicker_time.py

@@ -10,6 +10,7 @@
 # Local applications imports
 from mosqito.functions.loudness_zwicker.loudness_zwicker_shared import (
     calc_main_loudness,
+    calc_main_loudness_ea,
 )
 from mosqito.functions.loudness_zwicker.loudness_zwicker_nonlinear_decay import (
     calc_nl_loudness,
@@ -54,10 +55,14 @@ def loudness_zwicker_time(third_octave_levels, field_type):
     """
 
     # Calculate core loudness
-    num_sample_level = np.shape(third_octave_levels)[1]
-    core_loudness = np.zeros((21, num_sample_level))
-    for i in np.arange(num_sample_level - 1):
-        core_loudness[:, i] = calc_main_loudness(third_octave_levels[:, i], field_type)
+    # num_sample_level = np.shape(third_octave_levels)[1]
+    # core_loudness = np.zeros((21, num_sample_level))
+    # for i in np.arange(num_sample_level - 1):
+    #     core_loudness[:, i] = calc_main_loudness(third_octave_levels[:, i], field_type)
+
+    # Calculate core loudness (vectorized version, need debugging)
+    core_loudness = calc_main_loudness_ea(third_octave_levels, field_type)
+
     #
     # Nonlinearity
     core_loudness = calc_nl_loudness(core_loudness)