add paris fomula

Author: ahms5

docs/api_reference.rst

@@ -16,7 +16,7 @@ Modules
 .. toctree::
    :maxdepth: 1
 
-   modules/imkar
+   modules/imkar.utils
 
 
 .. _examples gallery: https://pyfar-gallery.readthedocs.io/en/latest/examples_gallery.html

docs/modules/imkar.rst docs/modules/imkar.utils.rstdocs/modules/imkar.rst renamed to docs/modules/imkar.utils.rst

@@ -1,7 +1,7 @@
-imkar
-=====
+imkar.utils
+===========
 
-.. automodule:: imkar
+.. automodule:: imkar.utils
    :members:
    :undoc-members:
    :show-inheritance:

imkar/__init__.py

@@ -5,3 +5,10 @@
 __author__ = """The pyfar developers"""
 __email__ = ''
 __version__ = '0.1.0'
+
+
+from . import utils
+
+__all__ = [
+    'utils',
+]

imkar/utils.py

@@ -0,0 +1,65 @@
+"""Utilities for imkar."""
+import numpy as np
+import pyfar as pf
+
+
+def paris_formula(coefficients, incident_directions):
+    r"""
+    Calculate the random-incidence coefficient
+    according to Paris formula.
+
+    The implementation follows the Equation 2.53 from [#]_ and is
+    discretized as:
+
+    .. math::
+        c_{rand} = \sum_{\Omega_S} c(\Omega_S) \cdot |\Omega_S \cdot n| \cdot w
+
+    with the ``coefficients`` :math:`c`, and the
+    area weights :math:`w` from the ``incident_directions``.
+    :math:`|\Omega_S \cdot n|` represent the cosine of the angle between the
+    surface normal and the incident direction.
+
+    .. note::
+        The incident directions should be
+        equally distributed to get a valid result.
+
+    Parameters
+    ----------
+    coefficients : pyfar.FrequencyData
+        coefficients for different incident directions. Its cshape
+        needs to be (..., n_incident_directions)
+    incident_directions : pyfar.Coordinates
+        Defines the incidence directions of each `coefficients` in a
+        Coordinates object. Its cshape needs to be (n_incident_directions). In
+        sperical coordinates the radii needs to be constant. The weights need
+        to reflect the area weights.
+
+    Returns
+    -------
+    random_coefficient : pyfar.FrequencyData
+        The random-incidence scattering coefficient.
+
+    References
+    ----------
+    .. [#]  H. Kuttruff, Room acoustics, Sixth edition. Boca Raton:
+            CRC Press/Taylor & Francis Group, 2017.
+    """
+    if not isinstance(coefficients, pf.FrequencyData):
+        raise ValueError("coefficients has to be FrequencyData")
+    if not isinstance(incident_directions, pf.Coordinates):
+        raise ValueError("incident_directions have to be None or Coordinates")
+    if incident_directions.cshape[0] != coefficients.cshape[-1]:
+        raise ValueError(
+            "the last dimension of coefficients needs be same as "
+            "the incident_directions.cshape.")
+
+    theta = incident_directions.colatitude
+    weight = np.cos(theta) * incident_directions.weights
+    norm = np.sum(weight)
+    coefficients_freq = np.swapaxes(coefficients.freq, -1, -2)
+    random_coefficient = pf.FrequencyData(
+        np.sum(coefficients_freq*weight/norm, axis=-1),
+        coefficients.frequencies,
+        comment='random-incidence coefficient',
+    )
+    return random_coefficient

tests/test_utils.py

@@ -0,0 +1,60 @@
+import pytest
+import numpy as np
+import pyfar as pf
+
+from imkar import utils
+
+
+@pytest.mark.parametrize("c_value",  [
+    (0), (0.2), (0.5), (0.8), (1)])
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_random_constant_coefficient(
+        c_value, frequencies):
+    incident_directions = pf.Coordinates.from_spherical_colatitude(
+        0, np.linspace(0, np.pi/2, 10), 1)
+    incident_directions.weights = np.sin(incident_directions.colatitude)
+    shape = np.append(incident_directions.cshape, len(frequencies))
+    coefficient = pf.FrequencyData(np.zeros(shape)+c_value, frequencies)
+    c_rand = utils.paris_formula(coefficient, incident_directions)
+    np.testing.assert_allclose(c_rand.freq, c_value)
+    assert c_rand.comment == 'random-incidence coefficient'
+
+
+def test_random_non_constant_coefficient():
+    incident_directions = pf.samplings.sph_gaussian(10)
+    incident_directions = incident_directions[incident_directions.z >= 0]
+    incident_cshape = incident_directions.cshape
+    s_value = np.arange(
+        incident_cshape[0]).reshape(incident_cshape) / incident_cshape[0]
+    theta = incident_directions.colatitude
+    actual_weight = np.cos(theta) * incident_directions.weights
+    actual_weight /= np.sum(actual_weight)
+    coefficient = pf.FrequencyData(s_value.reshape((50, 1)), [100])
+    c_rand = utils.paris_formula(coefficient, incident_directions)
+    desired = np.sum(s_value*actual_weight)
+    np.testing.assert_allclose(c_rand.freq, desired)
+
+
+def test_paris_formula_coefficients_type_error():
+    incident_directions = pf.samplings.sph_gaussian(10)
+    with pytest.raises(
+            ValueError, match="coefficients has to be FrequencyData"):
+        utils.paris_formula(np.zeros(10), incident_directions)
+
+
+def test_paris_formula_incident_directions_type_error():
+    coefficients = pf.FrequencyData(np.zeros((10, 1)), [100])
+    with pytest.raises(
+            ValueError,
+            match="incident_directions have to be None or Coordinates"):
+        utils.paris_formula(coefficients, None)
+
+
+def test_paris_formula_shape_mismatch_error():
+    incident_directions = pf.samplings.sph_gaussian(10)
+    coefficients = pf.FrequencyData(np.zeros((5, 1)), [100])
+    with pytest.raises(
+            ValueError,
+            match="the last dimension of coefficients needs be same as"):
+        utils.paris_formula(coefficients, incident_directions)