From 935ab12cc61b172a5b88212e158ed34165bd381e Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Tue, 29 Jul 2025 12:00:47 +0200 Subject: [PATCH 1/4] add paris fomula --- docs/api_reference.rst | 2 +- docs/modules/{imkar.rst => imkar.utils.rst} | 6 +- imkar/__init__.py | 7 +++ imkar/utils.py | 65 +++++++++++++++++++++ tests/test_utils.py | 60 +++++++++++++++++++ 5 files changed, 136 insertions(+), 4 deletions(-) rename docs/modules/{imkar.rst => imkar.utils.rst} (51%) create mode 100644 imkar/utils.py create mode 100644 tests/test_utils.py diff --git a/docs/api_reference.rst b/docs/api_reference.rst index 8f0383c..5556bd2 100644 --- a/docs/api_reference.rst +++ b/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 diff --git a/docs/modules/imkar.rst b/docs/modules/imkar.utils.rst similarity index 51% rename from docs/modules/imkar.rst rename to docs/modules/imkar.utils.rst index 0ea5287..89d6335 100644 --- a/docs/modules/imkar.rst +++ b/docs/modules/imkar.utils.rst @@ -1,7 +1,7 @@ -imkar -===== +imkar.utils +=========== -.. automodule:: imkar +.. automodule:: imkar.utils :members: :undoc-members: :show-inheritance: diff --git a/imkar/__init__.py b/imkar/__init__.py index 0dcfbd5..7d314cf 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -5,3 +5,10 @@ __author__ = """The pyfar developers""" __email__ = '' __version__ = '0.1.0' + + +from . import utils + +__all__ = [ + 'utils', +] diff --git a/imkar/utils.py b/imkar/utils.py new file mode 100644 index 0000000..c20d45f --- /dev/null +++ b/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 diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 0000000..02c72ac --- /dev/null +++ b/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) From 56f680f4630be955edcd5d50fb055b63d4489f9a Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Tue, 29 Jul 2025 12:50:23 +0200 Subject: [PATCH 2/4] Apply suggestions from code review Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- imkar/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imkar/utils.py b/imkar/utils.py index c20d45f..3c5add7 100644 --- a/imkar/utils.py +++ b/imkar/utils.py @@ -47,7 +47,7 @@ def paris_formula(coefficients, incident_directions): 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") + raise ValueError("incident_directions have to be Coordinates") if incident_directions.cshape[0] != coefficients.cshape[-1]: raise ValueError( "the last dimension of coefficients needs be same as " From 4240fdff444c764d281df97f9d7bc0c19bc272d9 Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Wed, 3 Dec 2025 14:00:05 +0100 Subject: [PATCH 3/4] Apply suggestions from code review Co-authored-by: Fabian Brinkmann --- imkar/utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/imkar/utils.py b/imkar/utils.py index 3c5add7..c3fc306 100644 --- a/imkar/utils.py +++ b/imkar/utils.py @@ -6,7 +6,7 @@ def paris_formula(coefficients, incident_directions): r""" Calculate the random-incidence coefficient - according to Paris formula. + according to the Paris formula. The implementation follows the Equation 2.53 from [#]_ and is discretized as: @@ -14,8 +14,8 @@ def paris_formula(coefficients, incident_directions): .. 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``. + 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. @@ -31,8 +31,8 @@ def paris_formula(coefficients, 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. + sperical coordinates the radii needs to be constant. The area weights must be + stored in ``ìncedent_directions.weights``. Returns ------- From a2b54a2c4b2d3c2a9def2631b9f34c2a50ce912e Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Wed, 3 Dec 2025 14:17:03 +0100 Subject: [PATCH 4/4] start working on review --- imkar/utils.py | 51 ++++++++++++++++++++++++-------------------------- 1 file changed, 24 insertions(+), 27 deletions(-) diff --git a/imkar/utils.py b/imkar/utils.py index c3fc306..0a7d526 100644 --- a/imkar/utils.py +++ b/imkar/utils.py @@ -3,41 +3,44 @@ import pyfar as pf -def paris_formula(coefficients, incident_directions): +def paris_formula(coefficients, colatitude_rad, area_weights): r""" Calculate the random-incidence coefficient according to the Paris formula. + The Paris formula computes the random-incidence coefficient based on + directional coefficients. It is valid for scattering and absorption + coefficients and requires equally distributed incident directions + to get a valid result. + 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 + c_{rand} = \sum_{\Omega_S} c(\Omega_S) \cdot \cos(\theta) \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 + :math:`\theta` represents 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 area weights must be - stored in ``ìncedent_directions.weights``. + needs to be (..., n_incident_directions). + colatitude_rad : pyfar.Coordinates + Defines the angle between the surface normal and the sound + incidence in radiant. Its cshape + needs to be (n_incident_directions). + area_weights : np.ndarray + Area weights for each incident direction. Its cshape + needs to be (n_incident_directions). Returns ------- random_coefficient : pyfar.FrequencyData - The random-incidence scattering coefficient. + The random-incidence coefficient. References ---------- @@ -46,20 +49,14 @@ def paris_formula(coefficients, incident_directions): """ 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 Coordinates") - if incident_directions.cshape[0] != coefficients.cshape[-1]: + if colatitude_rad.shape != area_weights.shape: raise ValueError( - "the last dimension of coefficients needs be same as " - "the incident_directions.cshape.") + "colatitude_rad and area_weights need to have the same shape.") - theta = incident_directions.colatitude - weight = np.cos(theta) * incident_directions.weights + theta = colatitude_rad + weight = np.cos(theta) * area_weights + weight = weight[..., np.newaxis] 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', - ) + random_coefficient = coefficients*weight/norm + random_coefficient.freq = np.sum(random_coefficient.freq, axis=-2) return random_coefficient