From 2616155c30b00d9779c4158207d1b193e4570d61 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Wed, 21 Sep 2022 14:40:03 +0200 Subject: [PATCH 01/39] add spherical integration over theta and phi --- imkar/__init__.py | 3 ++ imkar/imkar.py | 1 - imkar/integrate/__init__.py | 5 +++ imkar/integrate/integrate.py | 42 ++++++++++++++++++++++ tests/test_imkar.py | 24 ------------- tests/test_integrate.py | 68 ++++++++++++++++++++++++++++++++++++ 6 files changed, 118 insertions(+), 25 deletions(-) delete mode 100644 imkar/imkar.py create mode 100644 imkar/integrate/__init__.py create mode 100644 imkar/integrate/integrate.py delete mode 100644 tests/test_imkar.py create mode 100644 tests/test_integrate.py diff --git a/imkar/__init__.py b/imkar/__init__.py index 2d1e55d..a079d97 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -3,3 +3,6 @@ __author__ = """The pyfar developers""" __email__ = 'info@pyfar.org' __version__ = '0.1.0' + + +from . import integrate diff --git a/imkar/imkar.py b/imkar/imkar.py deleted file mode 100644 index dd0b80e..0000000 --- a/imkar/imkar.py +++ /dev/null @@ -1 +0,0 @@ -"""Main module.""" diff --git a/imkar/integrate/__init__.py b/imkar/integrate/__init__.py new file mode 100644 index 0000000..6c49894 --- /dev/null +++ b/imkar/integrate/__init__.py @@ -0,0 +1,5 @@ +from .integrate import (spherical) + + +__all__ = [ + 'spherical'] diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py new file mode 100644 index 0000000..f44bd00 --- /dev/null +++ b/imkar/integrate/integrate.py @@ -0,0 +1,42 @@ +import warnings +from scipy.integrate import trapezoid +import numpy as np + + +def spherical(data, coords): + """Intergrated over a sphere in spherical coordinates. + + .. math:: + + S &= \\int(\\int(data(\\phi, \\theta)d\\phi)sin(\\theta)\\dtheta), + + Parameters + ---------- + data : FrequencyData + Input data to integrate. Its shape is (..., #phi, #theta) + coords : Coordinates + A Coordinate object with the sample points corresponding to the data + values. It has the shape of (#phi, #theta). + + Returns + ------- + FrequencyData + the integrated object. Its shape is (...) and 2 dimentions shorter + then before. + """ + ptr = coords.get_sph(convention='top_colat') + phi = ptr[1, :, 0] + theta = ptr[:, 0, 1] + r_mean = np.mean(ptr[:, :, 2]) + if not np.allclose(ptr[:, :, 2], r_mean): + warnings.warn(r'all radii should be almost same') + last_close_to_0 = np.isclose(phi[-1], 0, atol=1e-16) + pi2_consistant = np.isclose(phi[-2]+np.diff(phi)[0], 2*np.pi, atol=1e-16) + if last_close_to_0 and pi2_consistant: + phi[-1] = 2*np.pi + weights = np.reshape(np.sin(theta), (len(theta), 1)) + result_raw = trapezoid(data.freq, x=phi, axis=-2) + result_raw1 = trapezoid(result_raw*weights, x=theta, axis=-2) + result = data.copy() + result.freq = result_raw1 + return result diff --git a/tests/test_imkar.py b/tests/test_imkar.py deleted file mode 100644 index f801097..0000000 --- a/tests/test_imkar.py +++ /dev/null @@ -1,24 +0,0 @@ -#!/usr/bin/env python - -"""Tests for `imkar` package.""" - -import pytest - - -from imkar import imkar - - -@pytest.fixture -def response(): - """Sample pytest fixture. - - See more at: http://doc.pytest.org/en/latest/fixture.html - """ - # import requests - # return requests.get('https://github.com/mberz/cookiecutter-pypackage') - - -def test_content(response): - """Sample pytest test function with the pytest fixture as an argument.""" - # from bs4 import BeautifulSoup - # assert 'GitHub' in BeautifulSoup(response.content).title.string diff --git a/tests/test_integrate.py b/tests/test_integrate.py new file mode 100644 index 0000000..013bc2f --- /dev/null +++ b/tests/test_integrate.py @@ -0,0 +1,68 @@ +import pytest +import numpy as np + +from imkar import integrate +from pyfar import FrequencyData, Coordinates + + +@pytest.mark.parametrize("radius", [1, 2, 3, 4, 5]) +@pytest.mark.parametrize("delta,rtol", [(1, 1e-4), (10, 5e-3)]) +def test_spherical_radius(radius, delta, rtol): + data, coords = _create_test_data( + np.arange(0, 360+delta, delta), + np.arange(0, 180+delta, delta)) + result = integrate.spherical(data, coords) + actual = np.real(result.freq[0, 0]) + desired = 4*np.pi + np.testing.assert_allclose(actual, desired, rtol=rtol) + + +@pytest.mark.parametrize("radius", [1, 2, 3, 4, 5]) +@pytest.mark.parametrize("delta,rtol", [(1, 1e-4), (10, 5e-3)]) +def test_spherical_limits(radius, delta, rtol): + data, coords = _create_test_data( + np.arange(0, 180+delta, delta), + np.arange(0, 90+delta, delta)) + result = integrate.spherical(data, coords) + desired = np.pi + actual = np.real(result.freq[0, 0]) + np.testing.assert_allclose(actual, desired, rtol=rtol) + + +@pytest.mark.parametrize("radius", [1, 2, 3, 4, 5]) +@pytest.mark.parametrize("delta,rtol", [(1, 1e-4), (10, 5e-3)]) +def test_spherical_data_dimensions(radius, delta, rtol): + data_raw = np.arange(1, 7).reshape(2, 3) + data, coords = _create_test_data( + np.arange(0, 360+delta, delta), + np.arange(0, 180+delta, delta), data_raw=data_raw) + result = integrate.spherical(data, coords) + desired = data_raw*4*np.pi + actual = np.squeeze(np.real(result.freq)) + np.testing.assert_allclose(actual, desired, rtol=rtol) + + +def test_spherical_warning_wrong_radius(): + data, coords = _create_test_data( + np.arange(0, 360+10, 10), + np.arange(0, 180+10, 10)) + sph = coords.get_sph() + # manipualte radius + sph[:, 1, 2] = 2 + coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) + with pytest.warns(Warning, match='radi'): + integrate.spherical(data, coords) + + +def _create_test_data(phi_deg, theta_deg, data_raw=None, n_bins=1, radius=1): + phi, theta = np.meshgrid(phi_deg, theta_deg) + coords = Coordinates( + phi, theta, np.ones(phi.shape)*radius, 'sph', 'top_colat', 'deg') + if data_raw is None: + data_raw = np.ones((1)) + for dim in phi.shape: + data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) + data_raw = np.repeat(data_raw[..., np.newaxis], n_bins, axis=-1) + freq_data = np.arange(1, n_bins+1)*100 + data = FrequencyData(data_raw, freq_data) + return data, coords From 5381d3bd8401ca8e8294f7dd8ed616aec260fe05 Mon Sep 17 00:00:00 2001 From: ahms5 <64446926+ahms5@users.noreply.github.com> Date: Thu, 22 Sep 2022 11:21:33 +0200 Subject: [PATCH 02/39] Apply suggestions from code review Co-authored-by: Marco Berzborn --- imkar/integrate/integrate.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index f44bd00..62bf41e 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -3,8 +3,8 @@ import numpy as np -def spherical(data, coords): - """Intergrated over a sphere in spherical coordinates. +def surface_sphere(data, coords): + """Integrate over a set of points sampled on a spherical surface. .. math:: @@ -13,16 +13,14 @@ def spherical(data, coords): Parameters ---------- data : FrequencyData - Input data to integrate. Its shape is (..., #phi, #theta) + Input data to integrate. Its cshape need to be (..., #phi, #theta) coords : Coordinates - A Coordinate object with the sample points corresponding to the data - values. It has the shape of (#phi, #theta). + Coordinate points at which the data is sampled. Its cshape needs to be (#phi, #theta), matching the cshape of `data`. Returns ------- FrequencyData - the integrated object. Its shape is (...) and 2 dimentions shorter - then before. + The integration result. Its dimension is reduced by the last two, which are consumed by the integration. """ ptr = coords.get_sph(convention='top_colat') phi = ptr[1, :, 0] From c0ddd18e44921f523dcdf26c619c8b93f2189021 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 14:55:58 +0200 Subject: [PATCH 03/39] apply suggestions from @mberz and renew tests --- imkar/integrate/__init__.py | 4 +- imkar/integrate/integrate.py | 23 +++++--- tests/test_integrate.py | 106 +++++++++++++++++++++++++---------- 3 files changed, 91 insertions(+), 42 deletions(-) diff --git a/imkar/integrate/__init__.py b/imkar/integrate/__init__.py index 6c49894..0663a81 100644 --- a/imkar/integrate/__init__.py +++ b/imkar/integrate/__init__.py @@ -1,5 +1,5 @@ -from .integrate import (spherical) +from .integrate import (surface_sphere) __all__ = [ - 'spherical'] + 'surface_sphere'] diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 62bf41e..228b27b 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -13,26 +13,31 @@ def surface_sphere(data, coords): Parameters ---------- data : FrequencyData - Input data to integrate. Its cshape need to be (..., #phi, #theta) + Input data to integrate. Its cshape need to be (..., #theta, #phi) coords : Coordinates - Coordinate points at which the data is sampled. Its cshape needs to be (#phi, #theta), matching the cshape of `data`. + Coordinate points at which the data is sampled. Its cshape needs to be + (#theta, #phi), matching the cshape of `data`. Returns ------- FrequencyData - The integration result. Its dimension is reduced by the last two, which are consumed by the integration. + The integration result. Its dimension is reduced by the last two, + which are consumed by the integration. """ - ptr = coords.get_sph(convention='top_colat') - phi = ptr[1, :, 0] - theta = ptr[:, 0, 1] - r_mean = np.mean(ptr[:, :, 2]) - if not np.allclose(ptr[:, :, 2], r_mean): + if coords.cshape != data.cshape[-2:]: + raise ValueError(r'Coordinates.cshape should be same as data.cshape[-2:]') + coords_spherical = coords.get_sph(convention='top_colat') + phi = coords_spherical[1, :, 0] + theta = coords_spherical[:, 0, 1] + radius = coords_spherical[:, :, 2] + r_mean = np.mean(radius) + if not np.allclose(radius, r_mean): warnings.warn(r'all radii should be almost same') last_close_to_0 = np.isclose(phi[-1], 0, atol=1e-16) pi2_consistant = np.isclose(phi[-2]+np.diff(phi)[0], 2*np.pi, atol=1e-16) if last_close_to_0 and pi2_consistant: phi[-1] = 2*np.pi - weights = np.reshape(np.sin(theta), (len(theta), 1)) + weights = np.transpose(np.atleast_2d(np.sin(theta))) result_raw = trapezoid(data.freq, x=phi, axis=-2) result_raw1 = trapezoid(result_raw*weights, x=theta, axis=-2) result = data.copy() diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 013bc2f..0085cb2 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -5,59 +5,103 @@ from pyfar import FrequencyData, Coordinates -@pytest.mark.parametrize("radius", [1, 2, 3, 4, 5]) -@pytest.mark.parametrize("delta,rtol", [(1, 1e-4), (10, 5e-3)]) -def test_spherical_radius(radius, delta, rtol): +@pytest.mark.parametrize( + " phi_0, phi_1, theta_0, theta_1, desired", [ + (0, 2*np.pi, 0, np.pi, 4*np.pi), + (0, np.pi, 0, np.pi, 2*np.pi), + (0, np.pi, 0, np.pi/2, np.pi), + (np.pi, 2*np.pi, 0, np.pi/2, np.pi), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2, np.pi), + (np.pi, 2*np.pi, 0, np.pi, 2*np.pi), + ]) +def test_surface_sphere_uniform_data_different_limits( + phi_0, phi_1, theta_0, theta_1, desired): + delta = np.deg2rad(10) data, coords = _create_test_data( - np.arange(0, 360+delta, delta), - np.arange(0, 180+delta, delta)) - result = integrate.spherical(data, coords) + np.arange(phi_0, phi_1+delta, delta), + np.arange(theta_0, theta_1+delta, delta)) + result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) - desired = 4*np.pi - np.testing.assert_allclose(actual, desired, rtol=rtol) + np.testing.assert_allclose(actual, desired, rtol=5e-3) -@pytest.mark.parametrize("radius", [1, 2, 3, 4, 5]) -@pytest.mark.parametrize("delta,rtol", [(1, 1e-4), (10, 5e-3)]) -def test_spherical_limits(radius, delta, rtol): +@pytest.mark.parametrize( + " phi_0, phi_1, theta_0, theta_1, desired", [ + (0, 2*np.pi, 0, np.pi, 4*np.pi), + (0, np.pi, 0, np.pi, 2*np.pi), + (0, np.pi, 0, np.pi/2, np.pi), + (np.pi, 2*np.pi, 0, np.pi/2, np.pi), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2, np.pi), + (np.pi, 2*np.pi, 0, np.pi, 2*np.pi), + ]) +def test_surface_sphere_data_preserve_shape_with_different_limits( + phi_0, phi_1, theta_0, theta_1, desired): + delta = np.deg2rad(10) + data_raw = np.arange(1, 7).reshape(2, 3) data, coords = _create_test_data( - np.arange(0, 180+delta, delta), - np.arange(0, 90+delta, delta)) - result = integrate.spherical(data, coords) - desired = np.pi - actual = np.real(result.freq[0, 0]) - np.testing.assert_allclose(actual, desired, rtol=rtol) + np.arange(phi_0, phi_1+delta, delta), + np.arange(theta_0, theta_1+delta, delta), + data_raw=data_raw) + result = integrate.surface_sphere(data, coords) + actual = np.squeeze(np.real(result.freq)) + np.testing.assert_allclose(actual, data_raw*desired, rtol=5e-3) -@pytest.mark.parametrize("radius", [1, 2, 3, 4, 5]) -@pytest.mark.parametrize("delta,rtol", [(1, 1e-4), (10, 5e-3)]) -def test_spherical_data_dimensions(radius, delta, rtol): - data_raw = np.arange(1, 7).reshape(2, 3) +@pytest.mark.parametrize( + " phi_0, phi_1, theta_0, theta_1", [ + (0, 2*np.pi, 0, np.pi), + (0, np.pi, 0, np.pi), + (0, np.pi, 0, np.pi/2), + (np.pi, 2*np.pi, 0, np.pi/2), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2), + (0, 2*np.pi, np.pi/2, np.pi), + (np.pi-.1, 2*np.pi-.1, 0, np.pi), + (np.pi, 2*np.pi, 0, np.pi), + ]) +def test_surface_sphere_nonuniform_data_different_limits( + phi_0, phi_1, theta_0, theta_1): + delta = np.deg2rad(10) data, coords = _create_test_data( - np.arange(0, 360+delta, delta), - np.arange(0, 180+delta, delta), data_raw=data_raw) - result = integrate.spherical(data, coords) - desired = data_raw*4*np.pi - actual = np.squeeze(np.real(result.freq)) - np.testing.assert_allclose(actual, desired, rtol=rtol) + np.arange(phi_0, phi_1+delta, delta), + np.arange(theta_0, theta_1+delta, delta)) + theta = coords.get_sph()[..., 1] + data.freq = data.freq[0, ...] + data.freq[..., 0] = np.cos(theta) + result = integrate.surface_sphere(data, coords) + actual = np.real(result.freq[0, 0]) + desired = ((-np.cos(2*theta_1)/2) - (-np.cos(2*theta_0)/2)) * (phi_1-phi_0) + np.testing.assert_allclose(actual, desired/2, rtol=5e-3, atol=0.02) def test_spherical_warning_wrong_radius(): + delta = np.deg2rad(10) data, coords = _create_test_data( - np.arange(0, 360+10, 10), - np.arange(0, 180+10, 10)) + np.arange(0, 2*np.pi+delta, delta), + np.arange(0, np.pi+delta, delta)) sph = coords.get_sph() # manipualte radius sph[:, 1, 2] = 2 coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.warns(Warning, match='radi'): - integrate.spherical(data, coords) + integrate.surface_sphere(data, coords) + + +def test_surface_sphere_error_invalid_coordinates_shape(): + delta = np.deg2rad(10) + data, coords = _create_test_data( + np.arange(0, 2*np.pi+delta, delta), + np.arange(0, np.pi+delta, delta)) + sph = coords.get_sph() + sph = sph[1:, :, :] + coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) + with pytest.raises(ValueError, match='Coordinates.cshape'): + integrate.surface_sphere(data, coords) def _create_test_data(phi_deg, theta_deg, data_raw=None, n_bins=1, radius=1): phi, theta = np.meshgrid(phi_deg, theta_deg) coords = Coordinates( - phi, theta, np.ones(phi.shape)*radius, 'sph', 'top_colat', 'deg') + phi, theta, np.ones(phi.shape)*radius, 'sph') if data_raw is None: data_raw = np.ones((1)) for dim in phi.shape: From 9663447204eb810baf9e8fe26d50592f2784efad Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 15:28:20 +0200 Subject: [PATCH 04/39] flake8 fixes and test improvements --- imkar/integrate/integrate.py | 3 ++- tests/test_integrate.py | 8 +++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 228b27b..1962f84 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -25,7 +25,8 @@ def surface_sphere(data, coords): which are consumed by the integration. """ if coords.cshape != data.cshape[-2:]: - raise ValueError(r'Coordinates.cshape should be same as data.cshape[-2:]') + raise ValueError( + r'Coordinates.cshape should be same as data.cshape[-2:]') coords_spherical = coords.get_sph(convention='top_colat') phi = coords_spherical[1, :, 0] theta = coords_spherical[:, 0, 1] diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 0085cb2..4b9b149 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -60,7 +60,7 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( ]) def test_surface_sphere_nonuniform_data_different_limits( phi_0, phi_1, theta_0, theta_1): - delta = np.deg2rad(10) + delta = np.pi/18 data, coords = _create_test_data( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) @@ -69,8 +69,10 @@ def test_surface_sphere_nonuniform_data_different_limits( data.freq[..., 0] = np.cos(theta) result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) - desired = ((-np.cos(2*theta_1)/2) - (-np.cos(2*theta_0)/2)) * (phi_1-phi_0) - np.testing.assert_allclose(actual, desired/2, rtol=5e-3, atol=0.02) + theta_upper = -np.cos(2*theta_1)/4 + theta_lower = -np.cos(2*theta_0)/4 + desired = (theta_upper - theta_lower) * (phi_1 - phi_0) + np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) def test_spherical_warning_wrong_radius(): From 647e2bb55d310546566fa164a3c0abbe0e27af89 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 15:44:35 +0200 Subject: [PATCH 05/39] add example and value name tests --- imkar/integrate/integrate.py | 16 ++++++++++++++++ tests/test_integrate.py | 4 ++-- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 1962f84..8195847 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -23,6 +23,22 @@ def surface_sphere(data, coords): FrequencyData The integration result. Its dimension is reduced by the last two, which are consumed by the integration. + + Examples + -------- + >>> from imkar import integrate + >>> from pyfar import FrequencyData, Coordinates + >>> import numpy as np + >>> phi, theta = np.meshgrid( + >>> np.linspace(0, 2*np.pi, num=361), + >>> np.linspace(0, np.pi, num=181)) + >>> coords = Coordinates( + >>> phi, theta, np.ones(phi.shape), 'sph') + >>> data_raw = np.ones(phi.shape) + >>> data = FrequencyData(data_raw[..., None], 100) + >>> result = integrate.surface_sphere(data, coords) + >>> result.freq + array([[12.56605162+0.j]]) """ if coords.cshape != data.cshape[-2:]: raise ValueError( diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 4b9b149..fb7db42 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -100,8 +100,8 @@ def test_surface_sphere_error_invalid_coordinates_shape(): integrate.surface_sphere(data, coords) -def _create_test_data(phi_deg, theta_deg, data_raw=None, n_bins=1, radius=1): - phi, theta = np.meshgrid(phi_deg, theta_deg) +def _create_test_data(phi_rad, theta_rad, data_raw=None, n_bins=1, radius=1): + phi, theta = np.meshgrid(phi_rad, theta_rad) coords = Coordinates( phi, theta, np.ones(phi.shape)*radius, 'sph') if data_raw is None: From 3c2584d4205bee876d0756d35d97dbf97a3bcd01 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 16:16:30 +0200 Subject: [PATCH 06/39] debug doc --- docs/conf.py | 15 ++++++++++++++- docs/modules.rst | 10 ++++++++++ docs/modules/imkar.integrate.rst | 7 +++++++ imkar/integrate/integrate.py | 4 ++-- requirements_dev.txt | 2 ++ 5 files changed, 35 insertions(+), 3 deletions(-) create mode 100644 docs/modules.rst create mode 100644 docs/modules/imkar.integrate.rst diff --git a/docs/conf.py b/docs/conf.py index 12123b8..4f9c2b2 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -31,7 +31,20 @@ # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.viewcode'] +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.viewcode', + 'sphinx.ext.napoleon', + 'sphinx.ext.autosummary', + 'matplotlib.sphinxext.plot_directive', + 'sphinx.ext.imgmath', + 'autodocsumm'] + +imgmath_latex_preamble = r'\usepackage{array}' + +# show tocs for classes and functions of modules using the autodocsumm +# package +autodoc_default_options = {'autosummary': True} # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] diff --git a/docs/modules.rst b/docs/modules.rst new file mode 100644 index 0000000..f06f880 --- /dev/null +++ b/docs/modules.rst @@ -0,0 +1,10 @@ +Modules +======= + +The following gives detailed information about all pyfar functions sorted +according to their modules. + +.. toctree:: + :maxdepth: 1 + + modules/imkar.integrate diff --git a/docs/modules/imkar.integrate.rst b/docs/modules/imkar.integrate.rst new file mode 100644 index 0000000..3357ef2 --- /dev/null +++ b/docs/modules/imkar.integrate.rst @@ -0,0 +1,7 @@ +imkar.integrate +============= + +.. automodule:: imkar.integrate + :members: + :undoc-members: + :show-inheritance: diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 8195847..8aa3d68 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -4,11 +4,11 @@ def surface_sphere(data, coords): - """Integrate over a set of points sampled on a spherical surface. + r"""Integrate over a set of points sampled on a spherical surface. .. math:: - S &= \\int(\\int(data(\\phi, \\theta)d\\phi)sin(\\theta)\\dtheta), + S = \int \int data(\phi, \theta) sin(\theta) d\phi d\theta Parameters ---------- diff --git a/requirements_dev.txt b/requirements_dev.txt index e63930f..6222c36 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -11,3 +11,5 @@ pytest numpy scipy pyfar +insipid-sphinx-theme +autodocsumm From e0e4aa1656357ac06f41e220fe8dac7ceac9c4f9 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 16:17:18 +0200 Subject: [PATCH 07/39] flake8 --- docs/conf.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 4f9c2b2..b6932e3 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -170,6 +170,3 @@ 'One line description of project.', 'Miscellaneous'), ] - - - From e83d947d7882c0ed148a586a0c032982409f2a78 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 16:18:39 +0200 Subject: [PATCH 08/39] draft scattering freefield --- imkar/__init__.py | 1 + imkar/scattering/__init__.py | 4 ++ imkar/scattering/coefficient.py | 81 +++++++++++++++++++++ tests/test_scattering_coefficient.py | 103 +++++++++++++++++++++++++++ 4 files changed, 189 insertions(+) create mode 100644 imkar/scattering/__init__.py create mode 100644 imkar/scattering/coefficient.py create mode 100644 tests/test_scattering_coefficient.py diff --git a/imkar/__init__.py b/imkar/__init__.py index a079d97..100419d 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -6,3 +6,4 @@ from . import integrate +from . import scattering diff --git a/imkar/scattering/__init__.py b/imkar/scattering/__init__.py new file mode 100644 index 0000000..d8a87cd --- /dev/null +++ b/imkar/scattering/__init__.py @@ -0,0 +1,4 @@ +from . import coefficient + +__all__ = [ + 'freefield'] diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py new file mode 100644 index 0000000..479245b --- /dev/null +++ b/imkar/scattering/coefficient.py @@ -0,0 +1,81 @@ +import numpy as np +import pyfar as pf +from imkar import integrate + + +def freefield(p_sample, p_reference, mics, incident_directions=None): + """ + This function calculates the free-field scattering coefficient using the + Mommertz correlation method [#]_. The input objects are the reflection + directivity of the test sample as well as a reference surface of equal + dimensions. The third input argument is the microphone array positions. + + If an optional fourth argument is given, the direction of the different + plane waves are expected. They will be used to determine the + random-incidence value using Paris' formula. + + Parameters + ---------- + p_sample : FrequencyData + Reflection directivity of the test sample. + p_reference : FrequencyData + Reflection directivity of the test reference sample. + mics : _type_ + _description_ + incident_directions : _type_ + _description_ + + Returns + ------- + s : FrequencyData + The output is the scattering coefficient for each plane wave direction. + If a second output argument is requested, the random-incidence value is + also returned. + s_rand : FrequencyData, optional + If a second output argument is requested, the random-incidence value is + also returned. + + References + ---------- + .. [#] E. Mommertz, „Determination of scattering coefficients from the + reflection directivity of architectural surfaces“, Applied + Acoustics, Bd. 60, Nr. 2, S. 201–203, Juni 2000, + doi: 10.1016/S0003-682X(99)00057-2. + + """ + + if not isinstance(p_sample, pf.FrequencyData): + raise ValueError("p_sample has to be FrequencyData") + if not isinstance(p_reference, pf.FrequencyData): + raise ValueError("p_reference has to be FrequencyData") + if not isinstance(mics, pf.Coordinates): + raise ValueError("mics have to be Coordinates") + if (incident_directions is not None) & \ + ~isinstance(incident_directions, pf.Coordinates): + raise ValueError("incident_directions have to be None or Coordinates") + + # calculate according to mommertz correlation method (inlcuding weights) + p_sample_sq = p_sample*p_sample + p_reference_sq = p_reference*p_reference + p_reference_conj = p_reference.copy() + p_reference_conj.freq = np.conj(p_reference_conj.freq) + p_cross = p_sample*p_reference_conj + + p_sample_sum = integrate.spherical_2d(p_sample_sq, mics) + p_ref_sum = integrate.spherical_2d(p_reference_sq, mics) + p_cross_sum = integrate.spherical_2d(p_cross, mics) + + s = (1 - ((p_cross_sum*p_cross_sum)/(p_sample_sum*p_ref_sum))) + s.comment = 'scattering coefficient' + # calculate random-incidence scattering coefficient + if incident_directions is not None: + sph = incident_directions.get_sph() + theta = sph[..., 1, None] + weight = np.sin(2*theta) # sin(2*theta) + norm = np.sum(weight) + s_rand = s*weight/norm + s_rand.freq = np.sum(s_rand.freq, axis=-2) + s_rand.freq = np.sum(s_rand.freq, axis=-2) + s_rand.comment = 'random-incidence scattering coefficient' + return s, s_rand + return s diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py new file mode 100644 index 0000000..dd6aff2 --- /dev/null +++ b/tests/test_scattering_coefficient.py @@ -0,0 +1,103 @@ +import pytest +import numpy as np + +from imkar import scattering +from pyfar import FrequencyData, Coordinates + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_1(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + p_sample = _create_frequencydata(mics.cshape, 1, frequencies) + p_reference = _create_frequencydata(mics.cshape, 0, frequencies) + p_sample.freq[5, 0, :] = 0 + p_reference.freq[5, 0, :] = np.sum(p_sample.freq.flatten())/2 + s = scattering.coefficient.freefield(p_sample, p_reference, mics) + np.testing.assert_allclose(s.freq, 1) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_05(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + p_sample = _create_frequencydata(mics.cshape, 0, frequencies) + p_reference = _create_frequencydata(mics.cshape, 0, frequencies) + p_sample.freq[5, 7, :] = 1 + p_sample.freq[5, 5, :] = 1 + p_reference.freq[5, 5, :] = 1 + s = scattering.coefficient.freefield(p_sample, p_reference, mics) + np.testing.assert_allclose(s.freq, 0.5) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_0(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + p_sample = _create_frequencydata(mics.cshape, 0, frequencies) + p_reference = _create_frequencydata(mics.cshape, 0, frequencies) + p_reference.freq[5, 0, :] = 1 + p_sample.freq[5, 0, :] = 1 + s = scattering.coefficient.freefield(p_sample, p_reference, mics) + np.testing.assert_allclose(s.freq, 0) + assert s.freq.shape[-1] == len(frequencies) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_0_with_inci(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + incident_directions = _create_coordinates( + np.arange(0, 360+30, 30), np.arange(0, 90+30, 30)) + data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() + p_sample = _create_frequencydata(data_shape, 0, frequencies) + p_reference = _create_frequencydata(data_shape, 0, frequencies) + p_reference.freq[:, :, 1, 2, :] = 1 + p_sample.freq[:, :, 1, 2, :] = 1 + s, s_rand = scattering.coefficient.freefield( + p_sample, p_reference, mics, incident_directions=incident_directions) + np.testing.assert_allclose(s.freq, 0) + np.testing.assert_allclose(s_rand.freq, 0) + assert s.freq.shape[-1] == len(frequencies) + assert s.cshape == incident_directions.cshape + assert s_rand.freq.shape[-1] == len(frequencies) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_1_with_inci(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + incident_directions = _create_coordinates( + np.arange(0, 360+30, 30), np.arange(0, 90+30, 30)) + data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() + p_sample = _create_frequencydata(data_shape, 0, frequencies) + p_reference = _create_frequencydata(data_shape, 0, frequencies) + p_reference.freq[:, :, 1, 2, :] = 1 + p_sample.freq[:, :, 2, 3, :] = 1 + s, s_rand = scattering.coefficient.freefield( + p_sample, p_reference, mics, incident_directions=incident_directions) + np.testing.assert_allclose(s.freq, 1) + np.testing.assert_allclose(s_rand.freq, 1) + assert s.freq.shape[-1] == len(frequencies) + assert s.cshape == incident_directions.cshape + assert s_rand.freq.shape[-1] == len(frequencies) + + +def _create_coordinates(phi_deg, theta_deg): + phi, theta = np.meshgrid(phi_deg, theta_deg) + coords = Coordinates( + phi, theta, np.ones(phi.shape), 'sph', 'top_colat', 'deg') + return coords + + +def _create_frequencydata(shape, value, frequencies): + frequencies = np.atleast_1d(frequencies) + shape_new = np.append(shape, frequencies.shape) + data = np.zeros(shape_new) + value + p_reference = FrequencyData(data, frequencies) + return p_reference From 3bb225260a9ae3b1111e8b85a141ea7e95947583 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 16:39:27 +0200 Subject: [PATCH 09/39] fix docs and bugs in scattering.freefield --- docs/modules.rst | 1 + docs/modules/imkar.scattering.coefficient.rst | 7 +++ imkar/scattering/coefficient.py | 45 +++++++++---------- tests/test_scattering_coefficient.py | 42 +++++++++++++++++ 4 files changed, 72 insertions(+), 23 deletions(-) create mode 100644 docs/modules/imkar.scattering.coefficient.rst diff --git a/docs/modules.rst b/docs/modules.rst index f06f880..3d9cc24 100644 --- a/docs/modules.rst +++ b/docs/modules.rst @@ -8,3 +8,4 @@ according to their modules. :maxdepth: 1 modules/imkar.integrate + modules/imkar.scattering.coefficient diff --git a/docs/modules/imkar.scattering.coefficient.rst b/docs/modules/imkar.scattering.coefficient.rst new file mode 100644 index 0000000..15fadb8 --- /dev/null +++ b/docs/modules/imkar.scattering.coefficient.rst @@ -0,0 +1,7 @@ +imkar.scattering.coefficient +============= + +.. automodule:: imkar.scattering.coefficient + :members: + :undoc-members: + :show-inheritance: diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index 479245b..43cccc3 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -6,34 +6,32 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): """ This function calculates the free-field scattering coefficient using the - Mommertz correlation method [#]_. The input objects are the reflection - directivity of the test sample as well as a reference surface of equal - dimensions. The third input argument is the microphone array positions. - - If an optional fourth argument is given, the direction of the different - plane waves are expected. They will be used to determine the - random-incidence value using Paris' formula. + Mommertz correlation method [#]_. Parameters ---------- p_sample : FrequencyData - Reflection directivity of the test sample. + Reflection sound pressure or directivity of the test sample. Its cshape + need to be (..., #theta_incident_directions, #phi_incident_directions, + #theta_mics, #phi_mics) p_reference : FrequencyData - Reflection directivity of the test reference sample. - mics : _type_ - _description_ - incident_directions : _type_ - _description_ + Reflection Reflection sound pressure or directivity of the test + reference sample. It has the same shape as p_sample. + mics : Coordinates + A Coordinate object with all microphone directions. Its cshape need to + be (#theta_mics, #phi_mics) + incident_directions : Coordinates, optinal + A Coordinate object with all incident directions. Its cshape need to + be (#theta_incident, #phi_incident). This input is optinal, if it is + the random-incidence scattering coefficient can be calculated. Returns ------- s : FrequencyData - The output is the scattering coefficient for each plane wave direction. - If a second output argument is requested, the random-incidence value is - also returned. + The scattering coefficient for each plane wave direction. s_rand : FrequencyData, optional - If a second output argument is requested, the random-incidence value is - also returned. + The random-incidence scattering coefficient. Is just retruned if the + incident_directions are given. References ---------- @@ -43,7 +41,7 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): doi: 10.1016/S0003-682X(99)00057-2. """ - + # check inputs if not isinstance(p_sample, pf.FrequencyData): raise ValueError("p_sample has to be FrequencyData") if not isinstance(p_reference, pf.FrequencyData): @@ -54,19 +52,20 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): ~isinstance(incident_directions, pf.Coordinates): raise ValueError("incident_directions have to be None or Coordinates") - # calculate according to mommertz correlation method (inlcuding weights) + # calculate according to mommertz correlation method p_sample_sq = p_sample*p_sample p_reference_sq = p_reference*p_reference p_reference_conj = p_reference.copy() p_reference_conj.freq = np.conj(p_reference_conj.freq) p_cross = p_sample*p_reference_conj - p_sample_sum = integrate.spherical_2d(p_sample_sq, mics) - p_ref_sum = integrate.spherical_2d(p_reference_sq, mics) - p_cross_sum = integrate.spherical_2d(p_cross, mics) + p_sample_sum = integrate.surface_sphere(p_sample_sq, mics) + p_ref_sum = integrate.surface_sphere(p_reference_sq, mics) + p_cross_sum = integrate.surface_sphere(p_cross, mics) s = (1 - ((p_cross_sum*p_cross_sum)/(p_sample_sum*p_ref_sum))) s.comment = 'scattering coefficient' + # calculate random-incidence scattering coefficient if incident_directions is not None: sph = incident_directions.get_sph() diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index dd6aff2..f325ac6 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -18,6 +18,25 @@ def test_freefield_1(frequencies): np.testing.assert_allclose(s.freq, 1) +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_wrong_input(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + p_sample = _create_frequencydata(mics.cshape, 1, frequencies) + p_reference = _create_frequencydata(mics.cshape, 0, frequencies) + p_sample.freq[5, 0, :] = 0 + p_reference.freq[5, 0, :] = np.sum(p_sample.freq.flatten())/2 + with pytest.raises(ValueError, match='p_sample'): + scattering.coefficient.freefield(1, p_reference, mics) + with pytest.raises(ValueError, match='p_reference'): + scattering.coefficient.freefield(p_sample, 1, mics, mics) + with pytest.raises(ValueError, match='mics'): + scattering.coefficient.freefield(p_sample, p_reference, 1, mics) + with pytest.raises(ValueError, match='incident_directions'): + scattering.coefficient.freefield(p_sample, p_reference, mics, 1) + + @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) def test_freefield_05(frequencies): @@ -88,6 +107,28 @@ def test_freefield_1_with_inci(frequencies): assert s_rand.freq.shape[-1] == len(frequencies) +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_05_with_inci(frequencies): + mics = _create_coordinates( + np.arange(0, 370, 10), np.arange(0, 100, 10)) + incident_directions = _create_coordinates( + np.arange(0, 360+30, 30), np.arange(0, 90+30, 30)) + data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() + p_sample = _create_frequencydata(data_shape, 0, frequencies) + p_reference = _create_frequencydata(data_shape, 0, frequencies) + p_sample.freq[:, :, 5, 7, :] = 1 + p_sample.freq[:, :, 5, 5, :] = 1 + p_reference.freq[:, :, 5, 5, :] = 1 + s, s_rand = scattering.coefficient.freefield( + p_sample, p_reference, mics, incident_directions=incident_directions) + np.testing.assert_allclose(s.freq, 0.5) + np.testing.assert_allclose(s_rand.freq, 0.5) + assert s.freq.shape[-1] == len(frequencies) + assert s.cshape == incident_directions.cshape + assert s_rand.freq.shape[-1] == len(frequencies) + + def _create_coordinates(phi_deg, theta_deg): phi, theta = np.meshgrid(phi_deg, theta_deg) coords = Coordinates( @@ -101,3 +142,4 @@ def _create_frequencydata(shape, value, frequencies): data = np.zeros(shape_new) + value p_reference = FrequencyData(data, frequencies) return p_reference + From 8cb5b1a7b95a996a7e2031e48ddbdce7935ba61f Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Thu, 22 Sep 2022 16:39:55 +0200 Subject: [PATCH 10/39] flake8 --- tests/test_scattering_coefficient.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index f325ac6..6f2f63d 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -142,4 +142,3 @@ def _create_frequencydata(shape, value, frequencies): data = np.zeros(shape_new) + value p_reference = FrequencyData(data, frequencies) return p_reference - From 12344e9fbc2d2868d18d823c0c86e71e6f1bb3c0 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Fri, 23 Sep 2022 09:39:55 +0200 Subject: [PATCH 11/39] separate calc random_incidence from freefield --- imkar/scattering/coefficient.py | 39 ++++++++++++---- tests/test_scattering_coefficient.py | 69 ++++++++++++++++++++++------ 2 files changed, 85 insertions(+), 23 deletions(-) diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index 43cccc3..6db322e 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -68,13 +68,36 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): # calculate random-incidence scattering coefficient if incident_directions is not None: - sph = incident_directions.get_sph() - theta = sph[..., 1, None] - weight = np.sin(2*theta) # sin(2*theta) - norm = np.sum(weight) - s_rand = s*weight/norm - s_rand.freq = np.sum(s_rand.freq, axis=-2) - s_rand.freq = np.sum(s_rand.freq, axis=-2) - s_rand.comment = 'random-incidence scattering coefficient' + s_rand = random_incidence(s, incident_directions) return s, s_rand return s + + +def random_incidence(s, incident_directions): + """This function claculates the random-incidence scattering coefficient + according to Paris formula. + + Parameters + ---------- + s : FrequencyData + The scattering coefficient for each plane wave direction. Its cshape + need to be (..., #theta, #phi) + incident_directions : Coordinates + A Coordinate object with all incident directions. Its cshape need to + be (#theta, #phi). + + + Returns + ------- + s_rand : FrequencyData + The random-incidence scattering coefficient. + """ + sph = incident_directions.get_sph() + theta = sph[..., 1, None] + weight = np.sin(2*theta) # sin(2*theta) + norm = np.sum(weight) + s_rand = s*weight/norm + s_rand.freq = np.sum(s_rand.freq, axis=-2) + s_rand.freq = np.sum(s_rand.freq, axis=-2) + s_rand.comment = 'random-incidence scattering coefficient' + return s_rand diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index 6f2f63d..46be8f8 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -9,7 +9,8 @@ ([100, 200]), ([100])]) def test_freefield_1(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=10)) p_sample = _create_frequencydata(mics.cshape, 1, frequencies) p_reference = _create_frequencydata(mics.cshape, 0, frequencies) p_sample.freq[5, 0, :] = 0 @@ -22,7 +23,8 @@ def test_freefield_1(frequencies): ([100, 200]), ([100])]) def test_freefield_wrong_input(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=10)) p_sample = _create_frequencydata(mics.cshape, 1, frequencies) p_reference = _create_frequencydata(mics.cshape, 0, frequencies) p_sample.freq[5, 0, :] = 0 @@ -41,7 +43,8 @@ def test_freefield_wrong_input(frequencies): ([100, 200]), ([100])]) def test_freefield_05(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=10)) p_sample = _create_frequencydata(mics.cshape, 0, frequencies) p_reference = _create_frequencydata(mics.cshape, 0, frequencies) p_sample.freq[5, 7, :] = 1 @@ -55,7 +58,8 @@ def test_freefield_05(frequencies): ([100, 200]), ([100])]) def test_freefield_0(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=10)) p_sample = _create_frequencydata(mics.cshape, 0, frequencies) p_reference = _create_frequencydata(mics.cshape, 0, frequencies) p_reference.freq[5, 0, :] = 1 @@ -69,9 +73,10 @@ def test_freefield_0(frequencies): ([100, 200]), ([100])]) def test_freefield_0_with_inci(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=10)) incident_directions = _create_coordinates( - np.arange(0, 360+30, 30), np.arange(0, 90+30, 30)) + np.linspace(0, 2*np.pi, num=8), np.linspace(0, np.pi/2, num=4)) data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() p_sample = _create_frequencydata(data_shape, 0, frequencies) p_reference = _create_frequencydata(data_shape, 0, frequencies) @@ -90,9 +95,10 @@ def test_freefield_0_with_inci(frequencies): ([100, 200]), ([100])]) def test_freefield_1_with_inci(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=19)) incident_directions = _create_coordinates( - np.arange(0, 360+30, 30), np.arange(0, 90+30, 30)) + np.linspace(0, 2*np.pi, num=8), np.linspace(0, np.pi/2, num=4)) data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() p_sample = _create_frequencydata(data_shape, 0, frequencies) p_reference = _create_frequencydata(data_shape, 0, frequencies) @@ -111,9 +117,10 @@ def test_freefield_1_with_inci(frequencies): ([100, 200]), ([100])]) def test_freefield_05_with_inci(frequencies): mics = _create_coordinates( - np.arange(0, 370, 10), np.arange(0, 100, 10)) + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi, num=10)) incident_directions = _create_coordinates( - np.arange(0, 360+30, 30), np.arange(0, 90+30, 30)) + np.linspace(0, 2*np.pi, num=8), np.linspace(0, np.pi/2, num=4)) data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() p_sample = _create_frequencydata(data_shape, 0, frequencies) p_reference = _create_frequencydata(data_shape, 0, frequencies) @@ -129,16 +136,48 @@ def test_freefield_05_with_inci(frequencies): assert s_rand.freq.shape[-1] == len(frequencies) -def _create_coordinates(phi_deg, theta_deg): - phi, theta = np.meshgrid(phi_deg, theta_deg) +@pytest.mark.parametrize("s_value", [ + (0), (0.2), (0.5), (0.8), (1)]) +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_random_incidence_constant_s(s_value, frequencies): + incident_directions = _create_coordinates( + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi/2, num=10)) + s = _create_frequencydata(incident_directions.cshape, s_value, frequencies) + s_rand = scattering.coefficient.random_incidence(s, incident_directions) + np.testing.assert_allclose(s_rand.freq, s_value) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_random_incidence_non_constant_s(frequencies): + incident_directions = _create_coordinates( + np.linspace(0, 2*np.pi, num=37), + np.linspace(0, np.pi/2, num=10)) + s_value = np.arange(37*10).reshape(incident_directions.cshape) / 370 + sph = incident_directions.get_sph() + actual_weight = np.sin(2*sph[..., 1]) # sin(2*theta) + actual_weight /= np.sum(actual_weight) + s = _create_frequencydata(incident_directions.cshape, s_value, frequencies) + s_rand = scattering.coefficient.random_incidence(s, incident_directions) + desired = np.sum(s_value*actual_weight) + np.testing.assert_allclose(s_rand.freq, desired) + + +def _create_coordinates(phi_rad, theta_rad): + phi, theta = np.meshgrid(phi_rad, theta_rad) coords = Coordinates( - phi, theta, np.ones(phi.shape), 'sph', 'top_colat', 'deg') + phi, theta, np.ones(phi.shape), 'sph') return coords -def _create_frequencydata(shape, value, frequencies): +def _create_frequencydata(shape, data_raw, frequencies): frequencies = np.atleast_1d(frequencies) shape_new = np.append(shape, frequencies.shape) - data = np.zeros(shape_new) + value + if hasattr(data_raw, "__len__"): + data_raw = np.repeat( + data_raw[..., np.newaxis], len(frequencies), axis=-1) + data = np.zeros(shape_new) + data_raw p_reference = FrequencyData(data, frequencies) return p_reference From 92fd3610adec9b23f62065a384abbf6f58b59da9 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Fri, 23 Sep 2022 11:24:47 +0200 Subject: [PATCH 12/39] clean up tests --- imkar/testing/__init__.py | 5 + imkar/testing/stub_utils.py | 37 ++++++++ tests/conftest.py | 52 +++++++++++ tests/test_integrate.py | 48 ++++------ tests/test_scattering_coefficient.py | 134 ++++++++++++--------------- 5 files changed, 172 insertions(+), 104 deletions(-) create mode 100644 imkar/testing/__init__.py create mode 100644 imkar/testing/stub_utils.py create mode 100644 tests/conftest.py diff --git a/imkar/testing/__init__.py b/imkar/testing/__init__.py new file mode 100644 index 0000000..0504466 --- /dev/null +++ b/imkar/testing/__init__.py @@ -0,0 +1,5 @@ +"""Testing sub-package for pyfar.""" + +from . import (stub_utils) + +__all__ = ['stub_utils'] diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py new file mode 100644 index 0000000..ed600da --- /dev/null +++ b/imkar/testing/stub_utils.py @@ -0,0 +1,37 @@ +""" +Contains tools to easily generate stubs for the most common pyfar Classes. + +Stubs are used instead of pyfar objects for testing functions that have pyfar +objects as input arguments. This makes testing such functions independent from +the pyfar objects themselves and helps to find bugs. +""" +import numpy as np +import pyfar as pf + + +def create_coordinates_sph(phi_rad, theta_rad, r=1): + phi, theta = np.meshgrid(phi_rad, theta_rad) + coords = pf.Coordinates( + phi, theta, np.ones(phi.shape)*r, 'sph') + return coords + + +def create_const_frequencydata_from_shape(shape, value, frequencies): + frequencies = np.atleast_1d(frequencies) + shape_new = np.append(shape, frequencies.shape) + data = np.zeros(shape_new) + value + p_reference = pf.FrequencyData(data, frequencies) + return p_reference + + +def create_frequencydata_from_shape(shape, data_raw, frequencies): + frequencies = np.atleast_1d(frequencies) + + if len(shape) > 0: + shape_new = np.append(shape, frequencies.shape) + for dim in shape: + data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) + data_raw = np.repeat(data_raw[..., np.newaxis], len(frequencies), axis=-1) + # data = np.zeros(shape_new) + data_raw + p_reference = pf.FrequencyData(data_raw, frequencies) + return p_reference diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..31b5ded --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,52 @@ +import pytest +import numpy as np +import pyfar as pf +from imkar.testing import stub_utils + + +@pytest.fixture +def coords_half_sphere_1_deg(): + delta = 1 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_5_deg(): + delta = 5 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_10_deg(): + delta = 10 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_30_deg(): + delta = 30 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_45_deg(): + delta = 45 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_sphere_10_deg(): + delta = 10 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi, num=int(180/delta)+1)) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index fb7db42..9582d05 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -3,6 +3,7 @@ from imkar import integrate from pyfar import FrequencyData, Coordinates +from imkar.testing import stub_utils @pytest.mark.parametrize( @@ -17,9 +18,11 @@ def test_surface_sphere_uniform_data_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.deg2rad(10) - data, coords = _create_test_data( + coords = stub_utils.create_coordinates_sph( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) + data = stub_utils.create_const_frequencydata_from_shape( + coords.cshape, 1, 100) result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) np.testing.assert_allclose(actual, desired, rtol=5e-3) @@ -38,10 +41,11 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.deg2rad(10) data_raw = np.arange(1, 7).reshape(2, 3) - data, coords = _create_test_data( + coords = stub_utils.create_coordinates_sph( np.arange(phi_0, phi_1+delta, delta), - np.arange(theta_0, theta_1+delta, delta), - data_raw=data_raw) + np.arange(theta_0, theta_1+delta, delta)) + data = stub_utils.create_frequencydata_from_shape( + coords.cshape, data_raw, 100) result = integrate.surface_sphere(data, coords) actual = np.squeeze(np.real(result.freq)) np.testing.assert_allclose(actual, data_raw*desired, rtol=5e-3) @@ -61,11 +65,12 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( def test_surface_sphere_nonuniform_data_different_limits( phi_0, phi_1, theta_0, theta_1): delta = np.pi/18 - data, coords = _create_test_data( + coords = stub_utils.create_coordinates_sph( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) + data = stub_utils.create_const_frequencydata_from_shape( + coords.cshape, 1, 100) theta = coords.get_sph()[..., 1] - data.freq = data.freq[0, ...] data.freq[..., 0] = np.cos(theta) result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) @@ -75,11 +80,10 @@ def test_surface_sphere_nonuniform_data_different_limits( np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) -def test_spherical_warning_wrong_radius(): - delta = np.deg2rad(10) - data, coords = _create_test_data( - np.arange(0, 2*np.pi+delta, delta), - np.arange(0, np.pi+delta, delta)) +def test_spherical_warning_wrong_radius(coords_sphere_10_deg): + coords = coords_sphere_10_deg + data = stub_utils.create_const_frequencydata_from_shape( + coords.cshape, 1, 100) sph = coords.get_sph() # manipualte radius sph[:, 1, 2] = 2 @@ -88,27 +92,13 @@ def test_spherical_warning_wrong_radius(): integrate.surface_sphere(data, coords) -def test_surface_sphere_error_invalid_coordinates_shape(): - delta = np.deg2rad(10) - data, coords = _create_test_data( - np.arange(0, 2*np.pi+delta, delta), - np.arange(0, np.pi+delta, delta)) +def test_surface_sphere_error_invalid_coordinates_shape(coords_sphere_10_deg): + coords = coords_sphere_10_deg + data = stub_utils.create_const_frequencydata_from_shape( + coords.cshape, 1, 100) sph = coords.get_sph() sph = sph[1:, :, :] coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.raises(ValueError, match='Coordinates.cshape'): integrate.surface_sphere(data, coords) - -def _create_test_data(phi_rad, theta_rad, data_raw=None, n_bins=1, radius=1): - phi, theta = np.meshgrid(phi_rad, theta_rad) - coords = Coordinates( - phi, theta, np.ones(phi.shape)*radius, 'sph') - if data_raw is None: - data_raw = np.ones((1)) - for dim in phi.shape: - data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) - data_raw = np.repeat(data_raw[..., np.newaxis], n_bins, axis=-1) - freq_data = np.arange(1, n_bins+1)*100 - data = FrequencyData(data_raw, freq_data) - return data, coords diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index 46be8f8..6103715 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -3,16 +3,17 @@ from imkar import scattering from pyfar import FrequencyData, Coordinates +from imkar.testing import stub_utils @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_1(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=10)) - p_sample = _create_frequencydata(mics.cshape, 1, frequencies) - p_reference = _create_frequencydata(mics.cshape, 0, frequencies) +def test_freefield_1(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 1, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 0, frequencies) p_sample.freq[5, 0, :] = 0 p_reference.freq[5, 0, :] = np.sum(p_sample.freq.flatten())/2 s = scattering.coefficient.freefield(p_sample, p_reference, mics) @@ -21,12 +22,12 @@ def test_freefield_1(frequencies): @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_wrong_input(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=10)) - p_sample = _create_frequencydata(mics.cshape, 1, frequencies) - p_reference = _create_frequencydata(mics.cshape, 0, frequencies) +def test_freefield_wrong_input(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 1, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 0, frequencies) p_sample.freq[5, 0, :] = 0 p_reference.freq[5, 0, :] = np.sum(p_sample.freq.flatten())/2 with pytest.raises(ValueError, match='p_sample'): @@ -41,12 +42,12 @@ def test_freefield_wrong_input(frequencies): @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_05(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=10)) - p_sample = _create_frequencydata(mics.cshape, 0, frequencies) - p_reference = _create_frequencydata(mics.cshape, 0, frequencies) +def test_freefield_05(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 0, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 0, frequencies) p_sample.freq[5, 7, :] = 1 p_sample.freq[5, 5, :] = 1 p_reference.freq[5, 5, :] = 1 @@ -56,12 +57,12 @@ def test_freefield_05(frequencies): @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_0(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=10)) - p_sample = _create_frequencydata(mics.cshape, 0, frequencies) - p_reference = _create_frequencydata(mics.cshape, 0, frequencies) +def test_freefield_0(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 0, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + mics.cshape, 0, frequencies) p_reference.freq[5, 0, :] = 1 p_sample.freq[5, 0, :] = 1 s = scattering.coefficient.freefield(p_sample, p_reference, mics) @@ -71,15 +72,15 @@ def test_freefield_0(frequencies): @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_0_with_inci(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=10)) - incident_directions = _create_coordinates( - np.linspace(0, 2*np.pi, num=8), np.linspace(0, np.pi/2, num=4)) +def test_freefield_0_with_inci( + frequencies, coords_half_sphere_10_deg, coords_half_sphere_30_deg): + mics = coords_half_sphere_10_deg + incident_directions = coords_half_sphere_30_deg data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() - p_sample = _create_frequencydata(data_shape, 0, frequencies) - p_reference = _create_frequencydata(data_shape, 0, frequencies) + p_sample = stub_utils.create_const_frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + data_shape, 0, frequencies) p_reference.freq[:, :, 1, 2, :] = 1 p_sample.freq[:, :, 1, 2, :] = 1 s, s_rand = scattering.coefficient.freefield( @@ -93,15 +94,15 @@ def test_freefield_0_with_inci(frequencies): @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_1_with_inci(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=19)) - incident_directions = _create_coordinates( - np.linspace(0, 2*np.pi, num=8), np.linspace(0, np.pi/2, num=4)) +def test_freefield_1_with_inci( + frequencies, coords_half_sphere_10_deg, coords_half_sphere_30_deg): + mics = coords_half_sphere_10_deg + incident_directions = coords_half_sphere_30_deg data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() - p_sample = _create_frequencydata(data_shape, 0, frequencies) - p_reference = _create_frequencydata(data_shape, 0, frequencies) + p_sample = stub_utils.create_const_frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + data_shape, 0, frequencies) p_reference.freq[:, :, 1, 2, :] = 1 p_sample.freq[:, :, 2, 3, :] = 1 s, s_rand = scattering.coefficient.freefield( @@ -115,15 +116,15 @@ def test_freefield_1_with_inci(frequencies): @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_freefield_05_with_inci(frequencies): - mics = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi, num=10)) - incident_directions = _create_coordinates( - np.linspace(0, 2*np.pi, num=8), np.linspace(0, np.pi/2, num=4)) +def test_freefield_05_with_inci( + frequencies, coords_half_sphere_10_deg, coords_half_sphere_30_deg): + mics = coords_half_sphere_10_deg + incident_directions = coords_half_sphere_30_deg data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten() - p_sample = _create_frequencydata(data_shape, 0, frequencies) - p_reference = _create_frequencydata(data_shape, 0, frequencies) + p_sample = stub_utils.create_const_frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference = stub_utils.create_const_frequencydata_from_shape( + data_shape, 0, frequencies) p_sample.freq[:, :, 5, 7, :] = 1 p_sample.freq[:, :, 5, 5, :] = 1 p_reference.freq[:, :, 5, 5, :] = 1 @@ -140,44 +141,27 @@ def test_freefield_05_with_inci(frequencies): (0), (0.2), (0.5), (0.8), (1)]) @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_random_incidence_constant_s(s_value, frequencies): - incident_directions = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi/2, num=10)) - s = _create_frequencydata(incident_directions.cshape, s_value, frequencies) +def test_random_incidence_constant_s( + s_value, frequencies, coords_half_sphere_10_deg): + incident_directions = coords_half_sphere_10_deg + s = stub_utils.create_const_frequencydata_from_shape( + incident_directions.cshape, s_value, frequencies) s_rand = scattering.coefficient.random_incidence(s, incident_directions) np.testing.assert_allclose(s_rand.freq, s_value) @pytest.mark.parametrize("frequencies", [ ([100, 200]), ([100])]) -def test_random_incidence_non_constant_s(frequencies): - incident_directions = _create_coordinates( - np.linspace(0, 2*np.pi, num=37), - np.linspace(0, np.pi/2, num=10)) +def test_random_incidence_non_constant_s( + frequencies, coords_half_sphere_10_deg): + incident_directions = coords_half_sphere_10_deg s_value = np.arange(37*10).reshape(incident_directions.cshape) / 370 sph = incident_directions.get_sph() actual_weight = np.sin(2*sph[..., 1]) # sin(2*theta) actual_weight /= np.sum(actual_weight) - s = _create_frequencydata(incident_directions.cshape, s_value, frequencies) + s = stub_utils.create_frequencydata_from_shape( + [], s_value, frequencies) s_rand = scattering.coefficient.random_incidence(s, incident_directions) desired = np.sum(s_value*actual_weight) np.testing.assert_allclose(s_rand.freq, desired) - -def _create_coordinates(phi_rad, theta_rad): - phi, theta = np.meshgrid(phi_rad, theta_rad) - coords = Coordinates( - phi, theta, np.ones(phi.shape), 'sph') - return coords - - -def _create_frequencydata(shape, data_raw, frequencies): - frequencies = np.atleast_1d(frequencies) - shape_new = np.append(shape, frequencies.shape) - if hasattr(data_raw, "__len__"): - data_raw = np.repeat( - data_raw[..., np.newaxis], len(frequencies), axis=-1) - data = np.zeros(shape_new) + data_raw - p_reference = FrequencyData(data, frequencies) - return p_reference From bb27c7ab5002f36961faf87f52171484aaf6641a Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Fri, 23 Sep 2022 12:17:59 +0200 Subject: [PATCH 13/39] cosmetics flake8 --- imkar/testing/stub_utils.py | 3 --- tests/test_integrate.py | 1 - 2 files changed, 4 deletions(-) diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py index ed600da..f6f8dff 100644 --- a/imkar/testing/stub_utils.py +++ b/imkar/testing/stub_utils.py @@ -26,12 +26,9 @@ def create_const_frequencydata_from_shape(shape, value, frequencies): def create_frequencydata_from_shape(shape, data_raw, frequencies): frequencies = np.atleast_1d(frequencies) - if len(shape) > 0: - shape_new = np.append(shape, frequencies.shape) for dim in shape: data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) data_raw = np.repeat(data_raw[..., np.newaxis], len(frequencies), axis=-1) - # data = np.zeros(shape_new) + data_raw p_reference = pf.FrequencyData(data_raw, frequencies) return p_reference diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 9582d05..1c2bd42 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -2,7 +2,6 @@ import numpy as np from imkar import integrate -from pyfar import FrequencyData, Coordinates from imkar.testing import stub_utils From 8886b4e9fd552660488293c021931bb67bc3e676 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Fri, 23 Sep 2022 12:20:02 +0200 Subject: [PATCH 14/39] flake8 --- imkar/scattering/__init__.py | 3 ++- tests/test_integrate.py | 1 - tests/test_scattering_coefficient.py | 2 -- 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/imkar/scattering/__init__.py b/imkar/scattering/__init__.py index d8a87cd..8cfb124 100644 --- a/imkar/scattering/__init__.py +++ b/imkar/scattering/__init__.py @@ -1,4 +1,5 @@ from . import coefficient __all__ = [ - 'freefield'] + 'freefield', + 'random_incidence'] diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 1c2bd42..fcab5b3 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -100,4 +100,3 @@ def test_surface_sphere_error_invalid_coordinates_shape(coords_sphere_10_deg): coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.raises(ValueError, match='Coordinates.cshape'): integrate.surface_sphere(data, coords) - diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index 6103715..6e7e9fb 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -2,7 +2,6 @@ import numpy as np from imkar import scattering -from pyfar import FrequencyData, Coordinates from imkar.testing import stub_utils @@ -164,4 +163,3 @@ def test_random_incidence_non_constant_s( s_rand = scattering.coefficient.random_incidence(s, incident_directions) desired = np.sum(s_value*actual_weight) np.testing.assert_allclose(s_rand.freq, desired) - From 575de91843890c4f0b79f3635ec91940a463e1b1 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Tue, 29 Nov 2022 14:14:05 +0100 Subject: [PATCH 15/39] bugfix for different arragements of angles in coords --- imkar/integrate/integrate.py | 9 +++++++-- tests/test_integrate.py | 25 +++++++++---------------- 2 files changed, 16 insertions(+), 18 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 8aa3d68..46d68ca 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -45,7 +45,12 @@ def surface_sphere(data, coords): r'Coordinates.cshape should be same as data.cshape[-2:]') coords_spherical = coords.get_sph(convention='top_colat') phi = coords_spherical[1, :, 0] - theta = coords_spherical[:, 0, 1] + theta = coords_spherical[:, 1, 1] + axis_index = -2 + if np.sum(np.diff(phi[1:-2])) < 1e-3: + phi = coords_spherical[:, 1, 0] + theta = coords_spherical[1, :, 1] + axis_index = -3 radius = coords_spherical[:, :, 2] r_mean = np.mean(radius) if not np.allclose(radius, r_mean): @@ -55,7 +60,7 @@ def surface_sphere(data, coords): if last_close_to_0 and pi2_consistant: phi[-1] = 2*np.pi weights = np.transpose(np.atleast_2d(np.sin(theta))) - result_raw = trapezoid(data.freq, x=phi, axis=-2) + result_raw = trapezoid(np.abs(data.freq), x=phi, axis=axis_index) result_raw1 = trapezoid(result_raw*weights, x=theta, axis=-2) result = data.copy() result.freq = result_raw1 diff --git a/tests/test_integrate.py b/tests/test_integrate.py index fb7db42..eac825e 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -48,30 +48,23 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( @pytest.mark.parametrize( - " phi_0, phi_1, theta_0, theta_1", [ - (0, 2*np.pi, 0, np.pi), - (0, np.pi, 0, np.pi), - (0, np.pi, 0, np.pi/2), - (np.pi, 2*np.pi, 0, np.pi/2), - (np.pi-.1, 2*np.pi-.1, 0, np.pi/2), - (0, 2*np.pi, np.pi/2, np.pi), - (np.pi-.1, 2*np.pi-.1, 0, np.pi), - (np.pi, 2*np.pi, 0, np.pi), + " phi_0, phi_1, theta_0, theta_1, desired", [ + (0, 2*np.pi, 0, np.pi, 4*np.pi), + (0, np.pi, 0, np.pi, 4*np.pi/2), + (0, np.pi, 0, np.pi/2, 4*np.pi/4), + (np.pi, 2*np.pi, 0, np.pi/2, 4*np.pi/4), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2, 4*np.pi/4), + (0, 2*np.pi, np.pi/2, np.pi, 4*np.pi/2), + (np.pi-.1, 2*np.pi-.1, 0, np.pi, 4*np.pi/2), ]) def test_surface_sphere_nonuniform_data_different_limits( - phi_0, phi_1, theta_0, theta_1): + phi_0, phi_1, theta_0, theta_1, desired): delta = np.pi/18 data, coords = _create_test_data( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) - theta = coords.get_sph()[..., 1] - data.freq = data.freq[0, ...] - data.freq[..., 0] = np.cos(theta) result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) - theta_upper = -np.cos(2*theta_1)/4 - theta_lower = -np.cos(2*theta_0)/4 - desired = (theta_upper - theta_lower) * (phi_1 - phi_0) np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) From 2880a59ae0ff64c8f974fe4125ae95352eb4bc81 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 07:49:19 +0100 Subject: [PATCH 16/39] include review --- imkar/integrate/integrate.py | 28 +++++++++++++++++++--------- tests/test_integrate.py | 6 +++--- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 46d68ca..4dbe35c 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -1,6 +1,7 @@ import warnings from scipy.integrate import trapezoid import numpy as np +import pyfar as pf def surface_sphere(data, coords): @@ -13,10 +14,11 @@ def surface_sphere(data, coords): Parameters ---------- data : FrequencyData - Input data to integrate. Its cshape need to be (..., #theta, #phi) + Input data to integrate. Its `cshape` needs to be (..., #theta, #phi) + or (..., #phi, #theta). coords : Coordinates Coordinate points at which the data is sampled. Its cshape needs to be - (#theta, #phi), matching the cshape of `data`. + (#theta, #phi) or (#phi, #theta), matching the cshape of `data`. Returns ------- @@ -40,9 +42,12 @@ def surface_sphere(data, coords): >>> result.freq array([[12.56605162+0.j]]) """ + if coords.cshape != data.cshape[-2:]: raise ValueError( - r'Coordinates.cshape should be same as data.cshape[-2:]') + f'Coordinates.cshape should be same as {data.cshape[-2:]}') + + # parse angles coords_spherical = coords.get_sph(convention='top_colat') phi = coords_spherical[1, :, 0] theta = coords_spherical[:, 1, 1] @@ -55,13 +60,18 @@ def surface_sphere(data, coords): r_mean = np.mean(radius) if not np.allclose(radius, r_mean): warnings.warn(r'all radii should be almost same') - last_close_to_0 = np.isclose(phi[-1], 0, atol=1e-16) - pi2_consistant = np.isclose(phi[-2]+np.diff(phi)[0], 2*np.pi, atol=1e-16) - if last_close_to_0 and pi2_consistant: + + # pf.Coordinates turns phi = 2*pi to 0, due to cylindircal + # for the integration the upper limit cannot be zero but 2pi + last_phi_is_zero = np.isclose(phi[-1], 0, atol=1e-16) + increasing_phi_to_2pi = np.isclose( + phi[-2]+np.diff(phi)[0], 2*np.pi, atol=1e-16) + if last_phi_is_zero and increasing_phi_to_2pi: phi[-1] = 2*np.pi + + # integrate over angles with wight for sperical integration weights = np.transpose(np.atleast_2d(np.sin(theta))) result_raw = trapezoid(np.abs(data.freq), x=phi, axis=axis_index) result_raw1 = trapezoid(result_raw*weights, x=theta, axis=-2) - result = data.copy() - result.freq = result_raw1 - return result + + return pf.FrequencyData(result_raw1, data.frequencies) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index eac825e..8e87b01 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -2,7 +2,7 @@ import numpy as np from imkar import integrate -from pyfar import FrequencyData, Coordinates +import pyfar as pf @pytest.mark.parametrize( @@ -95,7 +95,7 @@ def test_surface_sphere_error_invalid_coordinates_shape(): def _create_test_data(phi_rad, theta_rad, data_raw=None, n_bins=1, radius=1): phi, theta = np.meshgrid(phi_rad, theta_rad) - coords = Coordinates( + coords = pf.Coordinates( phi, theta, np.ones(phi.shape)*radius, 'sph') if data_raw is None: data_raw = np.ones((1)) @@ -103,5 +103,5 @@ def _create_test_data(phi_rad, theta_rad, data_raw=None, n_bins=1, radius=1): data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) data_raw = np.repeat(data_raw[..., np.newaxis], n_bins, axis=-1) freq_data = np.arange(1, n_bins+1)*100 - data = FrequencyData(data_raw, freq_data) + data = pf.FrequencyData(data_raw, freq_data) return data, coords From d351258b1764e6c3b2ed0814f98d97b98d6c30c0 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 07:58:24 +0100 Subject: [PATCH 17/39] include review --- imkar/integrate/integrate.py | 2 +- tests/test_integrate.py | 13 ++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index 4dbe35c..eb6462f 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -13,7 +13,7 @@ def surface_sphere(data, coords): Parameters ---------- - data : FrequencyData + data : FrequencyData, Signal Input data to integrate. Its `cshape` needs to be (..., #theta, #phi) or (..., #phi, #theta). coords : Coordinates diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 8e87b01..06ffb11 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -1,7 +1,6 @@ import pytest import numpy as np - -from imkar import integrate +import imkar as ik import pyfar as pf @@ -20,7 +19,7 @@ def test_surface_sphere_uniform_data_different_limits( data, coords = _create_test_data( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) - result = integrate.surface_sphere(data, coords) + result = ik.integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) np.testing.assert_allclose(actual, desired, rtol=5e-3) @@ -42,7 +41,7 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta), data_raw=data_raw) - result = integrate.surface_sphere(data, coords) + result = ik.integrate.surface_sphere(data, coords) actual = np.squeeze(np.real(result.freq)) np.testing.assert_allclose(actual, data_raw*desired, rtol=5e-3) @@ -63,7 +62,7 @@ def test_surface_sphere_nonuniform_data_different_limits( data, coords = _create_test_data( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) - result = integrate.surface_sphere(data, coords) + result = ik.integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) @@ -78,7 +77,7 @@ def test_spherical_warning_wrong_radius(): sph[:, 1, 2] = 2 coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.warns(Warning, match='radi'): - integrate.surface_sphere(data, coords) + ik.integrate.surface_sphere(data, coords) def test_surface_sphere_error_invalid_coordinates_shape(): @@ -90,7 +89,7 @@ def test_surface_sphere_error_invalid_coordinates_shape(): sph = sph[1:, :, :] coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.raises(ValueError, match='Coordinates.cshape'): - integrate.surface_sphere(data, coords) + ik.integrate.surface_sphere(data, coords) def _create_test_data(phi_rad, theta_rad, data_raw=None, n_bins=1, radius=1): From d8812d4b12830ecd6728c58afb72d86b8a74b617 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 11:45:16 +0100 Subject: [PATCH 18/39] merge scattering-freefield into integrate --- imkar/testing/__init__.py | 5 ++++ imkar/testing/stub_utils.py | 29 +++++++++++++++++++ tests/conftest.py | 51 ++++++++++++++++++++++++++++++++ tests/test_integrate.py | 58 +++++++++++++++---------------------- 4 files changed, 108 insertions(+), 35 deletions(-) create mode 100644 imkar/testing/__init__.py create mode 100644 imkar/testing/stub_utils.py create mode 100644 tests/conftest.py diff --git a/imkar/testing/__init__.py b/imkar/testing/__init__.py new file mode 100644 index 0000000..0504466 --- /dev/null +++ b/imkar/testing/__init__.py @@ -0,0 +1,5 @@ +"""Testing sub-package for pyfar.""" + +from . import (stub_utils) + +__all__ = ['stub_utils'] diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py new file mode 100644 index 0000000..9999676 --- /dev/null +++ b/imkar/testing/stub_utils.py @@ -0,0 +1,29 @@ +""" +Contains tools to easily generate stubs for the most common pyfar Classes. +Stubs are used instead of pyfar objects for testing functions that have pyfar +objects as input arguments. This makes testing such functions independent from +the pyfar objects themselves and helps to find bugs. +""" +import numpy as np +import pyfar as pf + + +def create_coordinates_sph(phi_rad, theta_rad, r=1): + phi, theta = np.meshgrid(phi_rad, theta_rad) + coords = pf.Coordinates( + phi, theta, np.ones(phi.shape)*r, 'sph') + return coords + + +def frequencydata_from_shape(shape, data_raw, frequencies): + frequencies = np.atleast_1d(frequencies) + shape_new = np.append(shape, frequencies.shape) + if hasattr(data_raw, "__len__"): # is array + if len(shape) > 0: + for dim in shape: + data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) + data = np.repeat(data_raw[..., np.newaxis], len(frequencies), axis=-1) + else: + data = np.zeros(shape_new) + data_raw + p_reference = pf.FrequencyData(data, frequencies) + return p_reference diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..54314d8 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,51 @@ +import pytest +import numpy as np +from imkar.testing import stub_utils + + +@pytest.fixture +def coords_half_sphere_1_deg(): + delta = 1 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_5_deg(): + delta = 5 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_10_deg(): + delta = 10 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_30_deg(): + delta = 30 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_half_sphere_45_deg(): + delta = 45 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi/2, num=int(90/delta)+1)) + + +@pytest.fixture +def coords_sphere_10_deg(): + delta = 10 + return stub_utils.create_coordinates_sph( + np.linspace(0, 2*np.pi, num=int(360/delta)+1), + np.linspace(0, np.pi, num=int(180/delta)+1)) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 06ffb11..f94ecfe 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -1,6 +1,8 @@ import pytest import numpy as np -import imkar as ik + +from imkar import integrate +from imkar.testing import stub_utils import pyfar as pf @@ -16,10 +18,12 @@ def test_surface_sphere_uniform_data_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.deg2rad(10) - data, coords = _create_test_data( + coords = stub_utils.create_coordinates_sph( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) - result = ik.integrate.surface_sphere(data, coords) + data = stub_utils.frequencydata_from_shape( + coords.cshape, 1, 100) + result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) np.testing.assert_allclose(actual, desired, rtol=5e-3) @@ -37,11 +41,11 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.deg2rad(10) data_raw = np.arange(1, 7).reshape(2, 3) - data, coords = _create_test_data( + coords = stub_utils.create_coordinates_sph( np.arange(phi_0, phi_1+delta, delta), - np.arange(theta_0, theta_1+delta, delta), - data_raw=data_raw) - result = ik.integrate.surface_sphere(data, coords) + np.arange(theta_0, theta_1+delta, delta)) + data = stub_utils.frequencydata_from_shape(coords.cshape, data_raw, 100) + result = integrate.surface_sphere(data, coords) actual = np.squeeze(np.real(result.freq)) np.testing.assert_allclose(actual, data_raw*desired, rtol=5e-3) @@ -59,48 +63,32 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( def test_surface_sphere_nonuniform_data_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.pi/18 - data, coords = _create_test_data( + coords = stub_utils.create_coordinates_sph( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) - result = ik.integrate.surface_sphere(data, coords) + data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) + result = integrate.surface_sphere(data, coords) actual = np.real(result.freq[0, 0]) np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) -def test_spherical_warning_wrong_radius(): - delta = np.deg2rad(10) - data, coords = _create_test_data( - np.arange(0, 2*np.pi+delta, delta), - np.arange(0, np.pi+delta, delta)) +def test_spherical_warning_wrong_radius(coords_sphere_10_deg): + coords = coords_sphere_10_deg + data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) sph = coords.get_sph() # manipualte radius sph[:, 1, 2] = 2 coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.warns(Warning, match='radi'): - ik.integrate.surface_sphere(data, coords) + integrate.surface_sphere(data, coords) -def test_surface_sphere_error_invalid_coordinates_shape(): - delta = np.deg2rad(10) - data, coords = _create_test_data( - np.arange(0, 2*np.pi+delta, delta), - np.arange(0, np.pi+delta, delta)) +def test_surface_sphere_error_invalid_coordinates_shape(coords_sphere_10_deg): + coords = coords_sphere_10_deg + data = stub_utils.frequencydata_from_shape( + coords.cshape, 1, 100) sph = coords.get_sph() sph = sph[1:, :, :] coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.raises(ValueError, match='Coordinates.cshape'): - ik.integrate.surface_sphere(data, coords) - - -def _create_test_data(phi_rad, theta_rad, data_raw=None, n_bins=1, radius=1): - phi, theta = np.meshgrid(phi_rad, theta_rad) - coords = pf.Coordinates( - phi, theta, np.ones(phi.shape)*radius, 'sph') - if data_raw is None: - data_raw = np.ones((1)) - for dim in phi.shape: - data_raw = np.repeat(data_raw[..., np.newaxis], dim, axis=-1) - data_raw = np.repeat(data_raw[..., np.newaxis], n_bins, axis=-1) - freq_data = np.arange(1, n_bins+1)*100 - data = pf.FrequencyData(data_raw, freq_data) - return data, coords + integrate.surface_sphere(data, coords) From 6ab317186e68a47b1a4405b92f5c954df524df80 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 11:46:59 +0100 Subject: [PATCH 19/39] flake 8 --- imkar/testing/stub_utils.py | 1 - tests/conftest.py | 1 - 2 files changed, 2 deletions(-) diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py index e880e53..fbd8944 100644 --- a/imkar/testing/stub_utils.py +++ b/imkar/testing/stub_utils.py @@ -5,7 +5,6 @@ objects as input arguments. This makes testing such functions independent from the pyfar objects themselves and helps to find bugs. """ -import collections.abc import numpy as np import pyfar as pf diff --git a/tests/conftest.py b/tests/conftest.py index 31b5ded..54314d8 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,6 +1,5 @@ import pytest import numpy as np -import pyfar as pf from imkar.testing import stub_utils From 1546bc2c3a8b428de29cbc1911d62e318c4755ed Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 11:48:11 +0100 Subject: [PATCH 20/39] flake8 --- imkar/testing/stub_utils.py | 1 - tests/test_integrate.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py index fbd8944..9999676 100644 --- a/imkar/testing/stub_utils.py +++ b/imkar/testing/stub_utils.py @@ -1,6 +1,5 @@ """ Contains tools to easily generate stubs for the most common pyfar Classes. - Stubs are used instead of pyfar objects for testing functions that have pyfar objects as input arguments. This makes testing such functions independent from the pyfar objects themselves and helps to find bugs. diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 74edeef..008cf01 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -3,7 +3,7 @@ from imkar import integrate from imkar.testing import stub_utils - +import pyfar as pf @pytest.mark.parametrize( " phi_0, phi_1, theta_0, theta_1, desired", [ From 37dab6421eedfc2e933b1712dbea2450daaf036c Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 11:53:19 +0100 Subject: [PATCH 21/39] flake 8 --- tests/test_integrate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index f94ecfe..74edeef 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -3,7 +3,6 @@ from imkar import integrate from imkar.testing import stub_utils -import pyfar as pf @pytest.mark.parametrize( From 56a79702dd4999179b3f5d773913f3b08d4e79d1 Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Wed, 30 Nov 2022 12:01:28 +0100 Subject: [PATCH 22/39] Apply suggestions from code review Co-authored-by: sikersten <70263411+sikersten@users.noreply.github.com> --- docs/modules/imkar.scattering.coefficient.rst | 2 +- imkar/scattering/coefficient.py | 18 +++++++++--------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/modules/imkar.scattering.coefficient.rst b/docs/modules/imkar.scattering.coefficient.rst index 15fadb8..f384e68 100644 --- a/docs/modules/imkar.scattering.coefficient.rst +++ b/docs/modules/imkar.scattering.coefficient.rst @@ -1,5 +1,5 @@ imkar.scattering.coefficient -============= +============================ .. automodule:: imkar.scattering.coefficient :members: diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index 9693b04..e66e507 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -5,19 +5,19 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): """ - This function calculates the free-field scattering coefficient using the + Calculate the free-field scattering coefficient using the Mommertz correlation method [#]_. Parameters ---------- - p_sample : FrequencyData + sample_pressure : FrequencyData Reflection sound pressure or directivity of the test sample. Its cshape need to be (..., #theta_incident_directions, #phi_incident_directions, #theta_mics, #phi_mics) - p_reference : FrequencyData + reference_pressure : FrequencyData Reflection Reflection sound pressure or directivity of the test - reference sample. It has the same shape as p_sample. - mics : Coordinates + reference sample. It has the same shape as `p_sample`. + mic_positions : Coordinates A Coordinate object with all microphone directions. Its cshape need to be (#theta_mics, #phi_mics) incident_directions : Coordinates, optinal @@ -27,7 +27,7 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): Returns ------- - s : FrequencyData + scattering_coefficients : FrequencyData The scattering coefficient for each plane wave direction. s_rand : FrequencyData, optional The random-incidence scattering coefficient. Is just retruned if the @@ -63,7 +63,7 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): p_ref_sum = integrate.surface_sphere(p_reference_sq, mics) p_cross_sum = integrate.surface_sphere(p_cross, mics) - s = (1 - ((p_cross_sum*p_cross_sum)/(p_sample_sum*p_ref_sum))) + s = 1 - ((p_cross_sum*p_cross_sum)/(p_sample_sum*p_ref_sum)) s.comment = 'scattering coefficient' # calculate random-incidence scattering coefficient @@ -74,12 +74,12 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): def random_incidence(s, incident_directions): - """This function claculates the random-incidence scattering coefficient + """Calculate the random-incidence scattering coefficient according to Paris formula. Parameters ---------- - s : FrequencyData + scattering_coefficient : FrequencyData The scattering coefficient for each plane wave direction. Its cshape need to be (..., #theta, #phi) incident_directions : Coordinates From 6e020a898e7c9fba0385bf2e4dcc48b0a4309f1d Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 12:58:02 +0100 Subject: [PATCH 23/39] include review for scattering freefield --- imkar/scattering/coefficient.py | 125 +++++++++++++++------------ tests/test_scattering_coefficient.py | 24 +++-- 2 files changed, 81 insertions(+), 68 deletions(-) diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index e66e507..844ec9d 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -3,35 +3,33 @@ from imkar import integrate -def freefield(p_sample, p_reference, mics, incident_directions=None): +def freefield(sample_pressure, reference_pressure, mic_positions): """ - Calculate the free-field scattering coefficient using the - Mommertz correlation method [#]_. + Calculate the free-field scattering coefficient for each incident direction + using the Mommertz correlation method [#]_. See :py:func:`random_incidence` + to calculate the random incidence scattering coefficient. + Parameters ---------- sample_pressure : FrequencyData - Reflection sound pressure or directivity of the test sample. Its cshape - need to be (..., #theta_incident_directions, #phi_incident_directions, - #theta_mics, #phi_mics) + Reflected sound pressure or directivity of the test sample. Its cshape + need to be (..., #theta_incident_positions, #phi_incident_positions, + #angle1, #angle2) reference_pressure : FrequencyData Reflection Reflection sound pressure or directivity of the test - reference sample. It has the same shape as `p_sample`. + reference sample. It has the same shape as `sample_pressure`. mic_positions : Coordinates - A Coordinate object with all microphone directions. Its cshape need to - be (#theta_mics, #phi_mics) - incident_directions : Coordinates, optinal - A Coordinate object with all incident directions. Its cshape need to - be (#theta_incident, #phi_incident). This input is optinal, if it is - the random-incidence scattering coefficient can be calculated. + A Coordinate object with all microphone positions. Its cshape need to + be (#angle1, #angle2). In sperical coordinates the radii need to + be constant. It need to be same for `sample_pressure` and + `reference_pressure`. Returns ------- scattering_coefficients : FrequencyData The scattering coefficient for each plane wave direction. - s_rand : FrequencyData, optional - The random-incidence scattering coefficient. Is just retruned if the - incident_directions are given. + References ---------- @@ -40,51 +38,62 @@ def freefield(p_sample, p_reference, mics, incident_directions=None): Acoustics, Bd. 60, Nr. 2, S. 201–203, Juni 2000, doi: 10.1016/S0003-682X(99)00057-2. + Examples + -------- + Read signal and orientations objects stored in a .far file. + + >>> import imkar + >>> scattering_coefficients = imkar.scattering.coefficient.freefield( + >>> sample_pressure, reference_pressure, mic_positions) + >>> random_s = imkar.scattering.coefficient.random_incidence( + >>> scattering_coefficients, incident_positions) """ # check inputs - if not isinstance(p_sample, pf.FrequencyData): - raise ValueError("p_sample has to be FrequencyData") - if not isinstance(p_reference, pf.FrequencyData): - raise ValueError("p_reference has to be FrequencyData") - if not isinstance(mics, pf.Coordinates): - raise ValueError("mics have to be Coordinates") - if (incident_directions is not None) & \ - ~isinstance(incident_directions, pf.Coordinates): - raise ValueError("incident_directions have to be None or Coordinates") - - # calculate according to mommertz correlation method - p_sample_sq = p_sample*p_sample - p_reference_sq = p_reference*p_reference - p_reference_conj = p_reference.copy() - p_reference_conj.freq = np.conj(p_reference_conj.freq) - p_cross = p_sample*p_reference_conj - - p_sample_sum = integrate.surface_sphere(p_sample_sq, mics) - p_ref_sum = integrate.surface_sphere(p_reference_sq, mics) - p_cross_sum = integrate.surface_sphere(p_cross, mics) - - s = 1 - ((p_cross_sum*p_cross_sum)/(p_sample_sum*p_ref_sum)) - s.comment = 'scattering coefficient' - - # calculate random-incidence scattering coefficient - if incident_directions is not None: - s_rand = random_incidence(s, incident_directions) - return s, s_rand - return s - - -def random_incidence(s, incident_directions): + if not isinstance(sample_pressure, pf.FrequencyData): + raise ValueError("sample_pressure has to be FrequencyData") + if not isinstance(reference_pressure, pf.FrequencyData): + raise ValueError("reference_pressure has to be FrequencyData") + if not isinstance(mic_positions, pf.Coordinates): + raise ValueError("microphone positions have to be Coordinates") + + # calculate according to mommertz correlation method Equation (5) + p_sample_abs = pf.FrequencyData( + np.abs(sample_pressure.freq), sample_pressure.frequencies) + p_reference_abs = pf.FrequencyData( + np.abs(reference_pressure.freq), reference_pressure.frequencies) + p_sample_sq = p_sample_abs*p_sample_abs + p_reference_sq = p_reference_abs*p_reference_abs + p_reference_conj = pf.FrequencyData( + np.conj(reference_pressure.freq), reference_pressure.frequencies) + p_cross = sample_pressure*p_reference_conj + + p_sample_sum = integrate.surface_sphere(p_sample_sq, mic_positions) + p_ref_sum = integrate.surface_sphere(p_reference_sq, mic_positions) + p_cross_sum = integrate.surface_sphere(p_cross, mic_positions) + p_cross_sum_abs = pf.FrequencyData( + np.abs(p_cross_sum.freq), p_cross_sum.frequencies) + + scattering_coefficients \ + = 1 - ((p_cross_sum_abs**2)/(p_sample_sum*p_ref_sum)) + scattering_coefficients.comment = 'scattering coefficient' + + return scattering_coefficients + + +def random_incidence(scattering_coefficient, incident_positions): """Calculate the random-incidence scattering coefficient - according to Paris formula. + according to Paris formula. Note that the incident directions should be + equally distributed to get a valid result. Parameters ---------- scattering_coefficient : FrequencyData The scattering coefficient for each plane wave direction. Its cshape - need to be (..., #theta, #phi) - incident_directions : Coordinates - A Coordinate object with all incident directions. Its cshape need to - be (#theta, #phi). + need to be (..., #angle1, #angle2) + incident_positions : Coordinates + Defines the incidence directions of each `scattering_coefficient` in a + Coordinates object. Its cshape need to be (#angle1, #angle2). In + sperical coordinates the radii need to be constant. Returns @@ -92,11 +101,17 @@ def random_incidence(s, incident_directions): s_rand : FrequencyData The random-incidence scattering coefficient. """ - sph = incident_directions.get_sph() + if not isinstance(scattering_coefficient, pf.FrequencyData): + raise ValueError("scattering_coefficient has to be FrequencyData") + if (incident_positions is not None) & \ + ~isinstance(incident_positions, pf.Coordinates): + raise ValueError("incident_positions have to be None or Coordinates") + + sph = incident_positions.get_sph() theta = sph[..., 1] weight = np.sin(2*theta) # sin(2*theta) norm = np.sum(weight) - s_rand = s*weight/norm + s_rand = scattering_coefficient*weight/norm s_rand.freq = np.sum(s_rand.freq, axis=-2) s_rand.freq = np.sum(s_rand.freq, axis=-2) s_rand.comment = 'random-incidence scattering coefficient' diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index cac8c7a..8357544 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -29,14 +29,12 @@ def test_freefield_wrong_input(frequencies, coords_half_sphere_10_deg): mics.cshape, 0, frequencies) p_sample.freq[5, 0, :] = 0 p_reference.freq[5, 0, :] = np.sum(p_sample.freq.flatten())/2 - with pytest.raises(ValueError, match='p_sample'): + with pytest.raises(ValueError, match='sample_pressure'): scattering.coefficient.freefield(1, p_reference, mics) - with pytest.raises(ValueError, match='p_reference'): - scattering.coefficient.freefield(p_sample, 1, mics, mics) - with pytest.raises(ValueError, match='mics'): - scattering.coefficient.freefield(p_sample, p_reference, 1, mics) - with pytest.raises(ValueError, match='incident_directions'): - scattering.coefficient.freefield(p_sample, p_reference, mics, 1) + with pytest.raises(ValueError, match='reference_pressure'): + scattering.coefficient.freefield(p_sample, 1, mics) + with pytest.raises(ValueError, match='microphone positions'): + scattering.coefficient.freefield(p_sample, p_reference, 1) @pytest.mark.parametrize("frequencies", [ @@ -82,8 +80,8 @@ def test_freefield_0_with_inci( data_shape, 0, frequencies) p_reference.freq[:, :, 1, 2, :] = 1 p_sample.freq[:, :, 1, 2, :] = 1 - s, s_rand = scattering.coefficient.freefield( - p_sample, p_reference, mics, incident_directions=incident_directions) + s = scattering.coefficient.freefield(p_sample, p_reference, mics) + s_rand = scattering.coefficient.random_incidence(s, incident_directions) np.testing.assert_allclose(s.freq, 0) np.testing.assert_allclose(s_rand.freq, 0) assert s.freq.shape[-1] == len(frequencies) @@ -104,8 +102,8 @@ def test_freefield_1_with_inci( data_shape, 0, frequencies) p_reference.freq[:, :, 1, 2, :] = 1 p_sample.freq[:, :, 2, 3, :] = 1 - s, s_rand = scattering.coefficient.freefield( - p_sample, p_reference, mics, incident_directions=incident_directions) + s = scattering.coefficient.freefield(p_sample, p_reference, mics) + s_rand = scattering.coefficient.random_incidence(s, incident_directions) np.testing.assert_allclose(s.freq, 1) np.testing.assert_allclose(s_rand.freq, 1) assert s.freq.shape[-1] == len(frequencies) @@ -127,8 +125,8 @@ def test_freefield_05_with_inci( p_sample.freq[:, :, 5, 7, :] = 1 p_sample.freq[:, :, 5, 5, :] = 1 p_reference.freq[:, :, 5, 5, :] = 1 - s, s_rand = scattering.coefficient.freefield( - p_sample, p_reference, mics, incident_directions=incident_directions) + s = scattering.coefficient.freefield(p_sample, p_reference, mics) + s_rand = scattering.coefficient.random_incidence(s, incident_directions) np.testing.assert_allclose(s.freq, 0.5) np.testing.assert_allclose(s_rand.freq, 0.5) assert s.freq.shape[-1] == len(frequencies) From a1fc6dbfba697a7220889916a9801b0394970bdc Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 13:00:47 +0100 Subject: [PATCH 24/39] fix random scattering --- imkar/scattering/coefficient.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index 844ec9d..eead6d2 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -40,7 +40,8 @@ def freefield(sample_pressure, reference_pressure, mic_positions): Examples -------- - Read signal and orientations objects stored in a .far file. + Calculate freefield scattering coefficients and then the random incidence + scattering coefficient. >>> import imkar >>> scattering_coefficients = imkar.scattering.coefficient.freefield( @@ -72,7 +73,7 @@ def freefield(sample_pressure, reference_pressure, mic_positions): p_cross_sum = integrate.surface_sphere(p_cross, mic_positions) p_cross_sum_abs = pf.FrequencyData( np.abs(p_cross_sum.freq), p_cross_sum.frequencies) - + scattering_coefficients \ = 1 - ((p_cross_sum_abs**2)/(p_sample_sum*p_ref_sum)) scattering_coefficients.comment = 'scattering coefficient' @@ -95,10 +96,9 @@ def random_incidence(scattering_coefficient, incident_positions): Coordinates object. Its cshape need to be (#angle1, #angle2). In sperical coordinates the radii need to be constant. - Returns ------- - s_rand : FrequencyData + random_scattering : FrequencyData The random-incidence scattering coefficient. """ if not isinstance(scattering_coefficient, pf.FrequencyData): @@ -111,8 +111,8 @@ def random_incidence(scattering_coefficient, incident_positions): theta = sph[..., 1] weight = np.sin(2*theta) # sin(2*theta) norm = np.sum(weight) - s_rand = scattering_coefficient*weight/norm - s_rand.freq = np.sum(s_rand.freq, axis=-2) - s_rand.freq = np.sum(s_rand.freq, axis=-2) - s_rand.comment = 'random-incidence scattering coefficient' - return s_rand + random_scattering = scattering_coefficient*weight/norm + random_scattering.freq = np.sum(random_scattering.freq, axis=-2) + random_scattering.freq = np.sum(random_scattering.freq, axis=-2) + random_scattering.comment = 'random-incidence scattering coefficient' + return random_scattering From 00bd4e8fa12a21dcb84acc6c5035f298bb51ae98 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 13:09:29 +0100 Subject: [PATCH 25/39] fix init --- imkar/scattering/__init__.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/imkar/scattering/__init__.py b/imkar/scattering/__init__.py index 8cfb124..5e3fc4d 100644 --- a/imkar/scattering/__init__.py +++ b/imkar/scattering/__init__.py @@ -1,5 +1,9 @@ -from . import coefficient +from .coefficient import ( + random_incidence, + freefield +) __all__ = [ 'freefield', - 'random_incidence'] + 'random_incidence' + ] From 9874c810c2e29b6448aa7da3c087c89de95d9ba1 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 13:13:03 +0100 Subject: [PATCH 26/39] fix flake8 --- imkar/__init__.py | 2 +- imkar/integrate/__init__.py | 4 +++- tests/test_import.py | 5 +++++ 3 files changed, 9 insertions(+), 2 deletions(-) create mode 100644 tests/test_import.py diff --git a/imkar/__init__.py b/imkar/__init__.py index a079d97..1691ca3 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -5,4 +5,4 @@ __version__ = '0.1.0' -from . import integrate +from . import integrate # noqa diff --git a/imkar/integrate/__init__.py b/imkar/integrate/__init__.py index 0663a81..f8ec396 100644 --- a/imkar/integrate/__init__.py +++ b/imkar/integrate/__init__.py @@ -1,4 +1,6 @@ -from .integrate import (surface_sphere) +from .integrate import ( + surface_sphere + ) __all__ = [ diff --git a/tests/test_import.py b/tests/test_import.py new file mode 100644 index 0000000..e6d07f7 --- /dev/null +++ b/tests/test_import.py @@ -0,0 +1,5 @@ +def test_import_imkar(): + try: + import imkar # noqa + except ImportError: + assert False From 73825ae8b25779887fa3c57757787e22c688e2d4 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 13:14:35 +0100 Subject: [PATCH 27/39] fix flake8 --- imkar/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imkar/__init__.py b/imkar/__init__.py index d33a575..b65e907 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -6,4 +6,4 @@ from . import integrate # noqa -from . import scattering +from . import scattering # noqa From d4ff85ebcdb8b59d62b3024dcd52b02ef414f57b Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 14:01:46 +0100 Subject: [PATCH 28/39] added test and rename sub_utils --- imkar/testing/stub_utils.py | 2 +- tests/conftest.py | 12 ++++---- tests/test_integrate.py | 6 ++-- tests/test_sub_utils.py | 58 +++++++++++++++++++++++++++++++++++++ 4 files changed, 68 insertions(+), 10 deletions(-) create mode 100644 tests/test_sub_utils.py diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py index 9999676..dc54762 100644 --- a/imkar/testing/stub_utils.py +++ b/imkar/testing/stub_utils.py @@ -8,7 +8,7 @@ import pyfar as pf -def create_coordinates_sph(phi_rad, theta_rad, r=1): +def spherical_coordinates(phi_rad, theta_rad, r=1): phi, theta = np.meshgrid(phi_rad, theta_rad) coords = pf.Coordinates( phi, theta, np.ones(phi.shape)*r, 'sph') diff --git a/tests/conftest.py b/tests/conftest.py index 54314d8..9a348d0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -6,7 +6,7 @@ @pytest.fixture def coords_half_sphere_1_deg(): delta = 1 - return stub_utils.create_coordinates_sph( + return stub_utils.spherical_coordinates( np.linspace(0, 2*np.pi, num=int(360/delta)+1), np.linspace(0, np.pi/2, num=int(90/delta)+1)) @@ -14,7 +14,7 @@ def coords_half_sphere_1_deg(): @pytest.fixture def coords_half_sphere_5_deg(): delta = 5 - return stub_utils.create_coordinates_sph( + return stub_utils.spherical_coordinates( np.linspace(0, 2*np.pi, num=int(360/delta)+1), np.linspace(0, np.pi/2, num=int(90/delta)+1)) @@ -22,7 +22,7 @@ def coords_half_sphere_5_deg(): @pytest.fixture def coords_half_sphere_10_deg(): delta = 10 - return stub_utils.create_coordinates_sph( + return stub_utils.spherical_coordinates( np.linspace(0, 2*np.pi, num=int(360/delta)+1), np.linspace(0, np.pi/2, num=int(90/delta)+1)) @@ -30,7 +30,7 @@ def coords_half_sphere_10_deg(): @pytest.fixture def coords_half_sphere_30_deg(): delta = 30 - return stub_utils.create_coordinates_sph( + return stub_utils.spherical_coordinates( np.linspace(0, 2*np.pi, num=int(360/delta)+1), np.linspace(0, np.pi/2, num=int(90/delta)+1)) @@ -38,7 +38,7 @@ def coords_half_sphere_30_deg(): @pytest.fixture def coords_half_sphere_45_deg(): delta = 45 - return stub_utils.create_coordinates_sph( + return stub_utils.spherical_coordinates( np.linspace(0, 2*np.pi, num=int(360/delta)+1), np.linspace(0, np.pi/2, num=int(90/delta)+1)) @@ -46,6 +46,6 @@ def coords_half_sphere_45_deg(): @pytest.fixture def coords_sphere_10_deg(): delta = 10 - return stub_utils.create_coordinates_sph( + return stub_utils.spherical_coordinates( np.linspace(0, 2*np.pi, num=int(360/delta)+1), np.linspace(0, np.pi, num=int(180/delta)+1)) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 74edeef..abf376f 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -17,7 +17,7 @@ def test_surface_sphere_uniform_data_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.deg2rad(10) - coords = stub_utils.create_coordinates_sph( + coords = stub_utils.spherical_coordinates( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) data = stub_utils.frequencydata_from_shape( @@ -40,7 +40,7 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.deg2rad(10) data_raw = np.arange(1, 7).reshape(2, 3) - coords = stub_utils.create_coordinates_sph( + coords = stub_utils.spherical_coordinates( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) data = stub_utils.frequencydata_from_shape(coords.cshape, data_raw, 100) @@ -62,7 +62,7 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( def test_surface_sphere_nonuniform_data_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.pi/18 - coords = stub_utils.create_coordinates_sph( + coords = stub_utils.spherical_coordinates( np.arange(phi_0, phi_1+delta, delta), np.arange(theta_0, theta_1+delta, delta)) data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) diff --git a/tests/test_sub_utils.py b/tests/test_sub_utils.py new file mode 100644 index 0000000..9a04a8b --- /dev/null +++ b/tests/test_sub_utils.py @@ -0,0 +1,58 @@ +import pytest +import numpy as np +import numpy.testing as npt + +from imkar.testing import stub_utils + + +@pytest.mark.parametrize( + " phi_0, phi_1, theta_0, theta_1, radius", [ + (0, 2*np.pi, 0, np.pi, 1), + (0, np.pi, 0, np.pi, 2), + (0, np.pi, 0, np.pi/2, 0.1), + (np.pi, 2*np.pi, 0, np.pi/2, np.pi), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2, np.pi), + (np.pi, 2*np.pi, 0, np.pi, 2), + ]) +def test_spherical_coordinates(phi_0, phi_1, theta_0, theta_1, radius): + delta = np.pi/18 + phi = np.arange(phi_0, phi_1+delta, delta) + theta = np.arange(theta_0, theta_1+delta, delta) + coords = stub_utils.spherical_coordinates(phi, theta, radius) + npt.assert_allclose(coords.get_sph()[1, :, 0], phi) + npt.assert_allclose(coords.get_sph()[:, 1, 1], theta) + npt.assert_allclose( + coords.get_sph()[:, :, 2], np.zeros(coords.cshape)+radius) + + +@pytest.mark.parametrize( + "shapes", [ + (3, 2), + (5, 2), + (3, 2, 7), + ]) +@pytest.mark.parametrize( + "data_in", [ + 0.1, + 0, + np.array([0.1, 1]), + np.arange(4*5).reshape(4, 5), + ]) +@pytest.mark.parametrize( + "frequency", [ + [100], + [100, 200], + ]) +def test_frequencydata_from_shape_data(shapes, data_in, frequency): + data = stub_utils.frequencydata_from_shape(shapes, data_in, frequency) + # npt.assert_allclose(data.freq, data_in) + if hasattr(data_in, '__len__'): + for idx in range(len(data_in.shape)): + assert data.cshape[idx] == data_in.shape[idx] + for idx in range(len(shapes)): + assert data.cshape[idx+len(data_in.shape)] == shapes[idx] + + else: + for idx in range(len(shapes)): + assert data.cshape[idx] == shapes[idx] + assert data.n_bins == len(frequency) From 42193aa15f767f0d1587f1c09a0d69777f28972d Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 14:17:18 +0100 Subject: [PATCH 29/39] add test for non uniform surface --- imkar/integrate/integrate.py | 2 +- tests/test_integrate.py | 31 ++++++++++++++++++++++++++++++- 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index eb6462f..d04572f 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -71,7 +71,7 @@ def surface_sphere(data, coords): # integrate over angles with wight for sperical integration weights = np.transpose(np.atleast_2d(np.sin(theta))) - result_raw = trapezoid(np.abs(data.freq), x=phi, axis=axis_index) + result_raw = trapezoid(data.freq, x=phi, axis=axis_index) result_raw1 = trapezoid(result_raw*weights, x=theta, axis=-2) return pf.FrequencyData(result_raw1, data.frequencies) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index abf376f..7c55453 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -59,7 +59,7 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( (0, 2*np.pi, np.pi/2, np.pi, 4*np.pi/2), (np.pi-.1, 2*np.pi-.1, 0, np.pi, 4*np.pi/2), ]) -def test_surface_sphere_nonuniform_data_different_limits( +def test_surface_sphere_uniform_data_different_limits( phi_0, phi_1, theta_0, theta_1, desired): delta = np.pi/18 coords = stub_utils.spherical_coordinates( @@ -91,3 +91,32 @@ def test_surface_sphere_error_invalid_coordinates_shape(coords_sphere_10_deg): coords.set_sph(sph[..., 0], sph[..., 1], sph[..., 2]) with pytest.raises(ValueError, match='Coordinates.cshape'): integrate.surface_sphere(data, coords) + + +@pytest.mark.parametrize( + " phi_0, phi_1, theta_0, theta_1", [ + (0, 2*np.pi, 0, np.pi), + (0, np.pi, 0, np.pi), + (0, np.pi, 0, np.pi/2), + (np.pi, 2*np.pi, 0, np.pi/2), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2), + (0, 2*np.pi, np.pi/2, np.pi), + (np.pi-.1, 2*np.pi-.1, 0, np.pi), + (np.pi, 2*np.pi, 0, np.pi), + ]) +def test_surface_sphere_nonuniform_data_different_limits( + phi_0, phi_1, theta_0, theta_1): + delta = np.pi/18 + coords = stub_utils.spherical_coordinates( + np.arange(phi_0, phi_1+delta, delta), + np.arange(theta_0, theta_1+delta, delta)) + data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) + theta = coords.get_sph()[..., 1] + # data.freq = data.freq[0, ...] + data.freq[..., 0] = np.cos(theta) + result = integrate.surface_sphere(data, coords) + actual = np.real(result.freq[0, 0]) + theta_upper = -np.cos(2*theta_1)/4 + theta_lower = -np.cos(2*theta_0)/4 + desired = (theta_upper - theta_lower) * (phi_1 - phi_0) + np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) From 11d019a2307e961a3ea142367b7af69587decb3e Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 14:25:31 +0100 Subject: [PATCH 30/39] remove duplicated test --- tests/test_integrate.py | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 7c55453..5dbe6e8 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -49,28 +49,6 @@ def test_surface_sphere_data_preserve_shape_with_different_limits( np.testing.assert_allclose(actual, data_raw*desired, rtol=5e-3) -@pytest.mark.parametrize( - " phi_0, phi_1, theta_0, theta_1, desired", [ - (0, 2*np.pi, 0, np.pi, 4*np.pi), - (0, np.pi, 0, np.pi, 4*np.pi/2), - (0, np.pi, 0, np.pi/2, 4*np.pi/4), - (np.pi, 2*np.pi, 0, np.pi/2, 4*np.pi/4), - (np.pi-.1, 2*np.pi-.1, 0, np.pi/2, 4*np.pi/4), - (0, 2*np.pi, np.pi/2, np.pi, 4*np.pi/2), - (np.pi-.1, 2*np.pi-.1, 0, np.pi, 4*np.pi/2), - ]) -def test_surface_sphere_uniform_data_different_limits( - phi_0, phi_1, theta_0, theta_1, desired): - delta = np.pi/18 - coords = stub_utils.spherical_coordinates( - np.arange(phi_0, phi_1+delta, delta), - np.arange(theta_0, theta_1+delta, delta)) - data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) - result = integrate.surface_sphere(data, coords) - actual = np.real(result.freq[0, 0]) - np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) - - def test_spherical_warning_wrong_radius(coords_sphere_10_deg): coords = coords_sphere_10_deg data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) From b5d4a9a0fc5f8638d2f05c8481ab7e7058fc3778 Mon Sep 17 00:00:00 2001 From: ahms5 Date: Wed, 30 Nov 2022 14:26:17 +0100 Subject: [PATCH 31/39] fix flake8 --- tests/test_integrate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 80d1683..5dbe6e8 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -4,6 +4,7 @@ from imkar import integrate from imkar.testing import stub_utils + @pytest.mark.parametrize( " phi_0, phi_1, theta_0, theta_1, desired", [ (0, 2*np.pi, 0, np.pi, 4*np.pi), From fe1337e6105df47d3931cd2bdc94203815daad73 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Tue, 24 Jan 2023 12:52:56 +0100 Subject: [PATCH 32/39] add shapes and f of pressure need to be same --- imkar/scattering/coefficient.py | 8 ++++++++ tests/test_scattering_coefficient.py | 5 +++++ 2 files changed, 13 insertions(+) diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index eead6d2..8789496 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -56,6 +56,14 @@ def freefield(sample_pressure, reference_pressure, mic_positions): raise ValueError("reference_pressure has to be FrequencyData") if not isinstance(mic_positions, pf.Coordinates): raise ValueError("microphone positions have to be Coordinates") + if not sample_pressure.cshape == reference_pressure.cshape: + raise ValueError( + "sample_presure and reference_pressure have to have the " + "same cshape.") + if any(sample_pressure.frequencies != reference_pressure.frequencies): + raise ValueError( + "sample_presure and reference_pressure have to have the " + "same frequencies.") # calculate according to mommertz correlation method Equation (5) p_sample_abs = pf.FrequencyData( diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py index 8357544..8b28dcd 100644 --- a/tests/test_scattering_coefficient.py +++ b/tests/test_scattering_coefficient.py @@ -35,6 +35,11 @@ def test_freefield_wrong_input(frequencies, coords_half_sphere_10_deg): scattering.coefficient.freefield(p_sample, 1, mics) with pytest.raises(ValueError, match='microphone positions'): scattering.coefficient.freefield(p_sample, p_reference, 1) + with pytest.raises(ValueError, match='cshape'): + scattering.coefficient.freefield(p_sample[:-2, ...], p_reference, mics) + with pytest.raises(ValueError, match='same frequencies'): + p_sample.frequencies[0] = 1 + scattering.coefficient.freefield(p_sample, p_reference, mics) @pytest.mark.parametrize("frequencies", [ From b3fe3bd6f97333e4f910a60e6b23dffd5f37c4cc Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Tue, 24 Jan 2023 13:25:35 +0100 Subject: [PATCH 33/39] fix documentation --- docs/modules/imkar.integrate.rst | 2 +- imkar/scattering/coefficient.py | 9 ++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/modules/imkar.integrate.rst b/docs/modules/imkar.integrate.rst index 3357ef2..8a6b0c6 100644 --- a/docs/modules/imkar.integrate.rst +++ b/docs/modules/imkar.integrate.rst @@ -1,5 +1,5 @@ imkar.integrate -============= +=============== .. automodule:: imkar.integrate :members: diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index 8789496..aa1dc5f 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -14,16 +14,15 @@ def freefield(sample_pressure, reference_pressure, mic_positions): ---------- sample_pressure : FrequencyData Reflected sound pressure or directivity of the test sample. Its cshape - need to be (..., #theta_incident_positions, #phi_incident_positions, - #angle1, #angle2) + need to be (..., #angle1, #angle2), see `mic_positions` for more detail. reference_pressure : FrequencyData - Reflection Reflection sound pressure or directivity of the test + Reflected sound pressure or directivity of the test reference sample. It has the same shape as `sample_pressure`. mic_positions : Coordinates A Coordinate object with all microphone positions. Its cshape need to be (#angle1, #angle2). In sperical coordinates the radii need to - be constant. It need to be same for `sample_pressure` and - `reference_pressure`. + be constant. Microphone positions need to be same for + `sample_pressure` and `reference_pressure`. Returns ------- From 36c28bfba49d725e77a5bc6e03e51e952e129123 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Tue, 24 Jan 2023 13:28:58 +0100 Subject: [PATCH 34/39] fix flake8 --- imkar/scattering/coefficient.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index aa1dc5f..bfc0bac 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -14,7 +14,8 @@ def freefield(sample_pressure, reference_pressure, mic_positions): ---------- sample_pressure : FrequencyData Reflected sound pressure or directivity of the test sample. Its cshape - need to be (..., #angle1, #angle2), see `mic_positions` for more detail. + need to be (..., #angle1, #angle2), see `mic_positions` for more + detail. reference_pressure : FrequencyData Reflected sound pressure or directivity of the test reference sample. It has the same shape as `sample_pressure`. From 0366f97b4bb9c52de9dd6ec3ca3c7ce30a344c65 Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Tue, 7 Feb 2023 09:05:50 +0100 Subject: [PATCH 35/39] init --- imkar/__init__.py | 1 + imkar/diffusion/__init__.py | 7 +++ imkar/diffusion/coefficient.py | 104 +++++++++++++++++++++++++++++++++ 3 files changed, 112 insertions(+) create mode 100644 imkar/diffusion/__init__.py create mode 100644 imkar/diffusion/coefficient.py diff --git a/imkar/__init__.py b/imkar/__init__.py index b65e907..8756e6c 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -7,3 +7,4 @@ from . import integrate # noqa from . import scattering # noqa +from . import diffusion # noqa diff --git a/imkar/diffusion/__init__.py b/imkar/diffusion/__init__.py new file mode 100644 index 0000000..7a2b700 --- /dev/null +++ b/imkar/diffusion/__init__.py @@ -0,0 +1,7 @@ +from .coefficient import ( + freefield + ) + + +__all__ = [ + 'freefield'] diff --git a/imkar/diffusion/coefficient.py b/imkar/diffusion/coefficient.py new file mode 100644 index 0000000..d089953 --- /dev/null +++ b/imkar/diffusion/coefficient.py @@ -0,0 +1,104 @@ +import numpy as np +import pyfar as pf +from imkar import integrate + + +def freefield(sample_pressure, mic_positions): + """ + Calculate the free-field diffusion coefficient for each incident direction + using the Mommertz correlation method [1]_. See :py:func:`random_incidence` + to calculate the random incidence diffusion coefficient. + + + Parameters + ---------- + sample_pressure : FrequencyData + Reflected sound pressure or directivity of the test sample. Its cshape + need to be (..., #theta_incident_positions, #phi_incident_positions, + #angle1, #angle2) + mic_positions : Coordinates + A Coordinate object with all microphone positions. Its cshape need to + be (#angle1, #angle2). In sperical coordinates the radii need to + be constant. It need to be same for `sample_pressure`. + + Returns + ------- + diffusion_coefficients : FrequencyData + The diffusion coefficient for each plane wave direction. + + + References + ---------- + .. [1] ISO 17497-2:2012, Sound-scattering properties of surfaces. + Part 2: Measurement of the directional diffusion coefficient in a + free field. Geneva, Switzerland: International Organization for + Standards, 2012. + + + Examples + -------- + Calculate freefield diffusion coefficients and then the random incidence + diffusion coefficient. + + >>> import imkar as ik + >>> diffusion_coefficients = ik.diffusion.coefficient.freefield( + >>> sample_pressure, mic_positions) + >>> random_d = ik.diffusion.coefficient.random_incidence( + >>> diffusion_coefficients, incident_positions) + """ + # check inputs + if not isinstance(sample_pressure, pf.FrequencyData): + raise ValueError("sample_pressure has to be FrequencyData") + if not isinstance(mic_positions, pf.Coordinates): + raise ValueError("microphone positions have to be Coordinates") + + # calculate according to mommertz correlation method Equation (6) + p_sample_sq = pf.FrequencyData( + np.abs(sample_pressure.freq)**2, sample_pressure.frequencies) + + p_sample_sum_sq = integrate.surface_sphere( + p_sample_sq*p_sample_sq, mic_positions) + p_sample_sq_sum = integrate.surface_sphere(p_sample_sq, mic_positions)**2 + n = 1 + diffusion_coefficients \ + = (p_sample_sum_sq - p_sample_sq_sum) / ((n-1) * p_sample_sq_sum) + diffusion_coefficients.comment = 'diffusion coefficients' + + return diffusion_coefficients + + +def random_incidence(scattering_coefficient, incident_positions): + """Calculate the random-incidence scattering coefficient + according to Paris formula. Note that the incident directions should be + equally distributed to get a valid result. + + Parameters + ---------- + scattering_coefficient : FrequencyData + The scattering coefficient for each plane wave direction. Its cshape + need to be (..., #angle1, #angle2) + incident_positions : Coordinates + Defines the incidence directions of each `scattering_coefficient` in a + Coordinates object. Its cshape need to be (#angle1, #angle2). In + sperical coordinates the radii need to be constant. + + Returns + ------- + random_scattering : FrequencyData + The random-incidence scattering coefficient. + """ + if not isinstance(scattering_coefficient, pf.FrequencyData): + raise ValueError("scattering_coefficient has to be FrequencyData") + if (incident_positions is not None) & \ + ~isinstance(incident_positions, pf.Coordinates): + raise ValueError("incident_positions have to be None or Coordinates") + + sph = incident_positions.get_sph() + theta = sph[..., 1] + weight = np.sin(2*theta) # sin(2*theta) + norm = np.sum(weight) + random_scattering = scattering_coefficient*weight/norm + random_scattering.freq = np.sum(random_scattering.freq, axis=-2) + random_scattering.freq = np.sum(random_scattering.freq, axis=-2) + random_scattering.comment = 'random-incidence scattering coefficient' + return random_scattering From 698f1af3f5744b537a920a6fab10df324760681c Mon Sep 17 00:00:00 2001 From: Anne Heimes Date: Tue, 7 Feb 2023 09:08:56 +0100 Subject: [PATCH 36/39] fix doc --- imkar/scattering/coefficient.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py index bfc0bac..6638027 100644 --- a/imkar/scattering/coefficient.py +++ b/imkar/scattering/coefficient.py @@ -6,7 +6,7 @@ def freefield(sample_pressure, reference_pressure, mic_positions): """ Calculate the free-field scattering coefficient for each incident direction - using the Mommertz correlation method [#]_. See :py:func:`random_incidence` + using the Mommertz correlation method [1]_. See :py:func:`random_incidence` to calculate the random incidence scattering coefficient. @@ -33,7 +33,7 @@ def freefield(sample_pressure, reference_pressure, mic_positions): References ---------- - .. [#] E. Mommertz, „Determination of scattering coefficients from the + .. [1] E. Mommertz, „Determination of scattering coefficients from the reflection directivity of architectural surfaces“, Applied Acoustics, Bd. 60, Nr. 2, S. 201–203, Juni 2000, doi: 10.1016/S0003-682X(99)00057-2. From c53e36f90f0af4a89984cd87497d40155630ffe2 Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Thu, 23 Feb 2023 16:51:27 +0100 Subject: [PATCH 37/39] remove random incidence --- imkar/diffusion/coefficient.py | 41 ++-------------------------------- 1 file changed, 2 insertions(+), 39 deletions(-) diff --git a/imkar/diffusion/coefficient.py b/imkar/diffusion/coefficient.py index d089953..751a83d 100644 --- a/imkar/diffusion/coefficient.py +++ b/imkar/diffusion/coefficient.py @@ -37,7 +37,7 @@ def freefield(sample_pressure, mic_positions): Examples -------- - Calculate freefield diffusion coefficients and then the random incidence + Calculate free-field diffusion coefficients and then the random incidence diffusion coefficient. >>> import imkar as ik @@ -52,7 +52,7 @@ def freefield(sample_pressure, mic_positions): if not isinstance(mic_positions, pf.Coordinates): raise ValueError("microphone positions have to be Coordinates") - # calculate according to mommertz correlation method Equation (6) + # calculate according to Mommertz correlation method Equation (6) p_sample_sq = pf.FrequencyData( np.abs(sample_pressure.freq)**2, sample_pressure.frequencies) @@ -65,40 +65,3 @@ def freefield(sample_pressure, mic_positions): diffusion_coefficients.comment = 'diffusion coefficients' return diffusion_coefficients - - -def random_incidence(scattering_coefficient, incident_positions): - """Calculate the random-incidence scattering coefficient - according to Paris formula. Note that the incident directions should be - equally distributed to get a valid result. - - Parameters - ---------- - scattering_coefficient : FrequencyData - The scattering coefficient for each plane wave direction. Its cshape - need to be (..., #angle1, #angle2) - incident_positions : Coordinates - Defines the incidence directions of each `scattering_coefficient` in a - Coordinates object. Its cshape need to be (#angle1, #angle2). In - sperical coordinates the radii need to be constant. - - Returns - ------- - random_scattering : FrequencyData - The random-incidence scattering coefficient. - """ - if not isinstance(scattering_coefficient, pf.FrequencyData): - raise ValueError("scattering_coefficient has to be FrequencyData") - if (incident_positions is not None) & \ - ~isinstance(incident_positions, pf.Coordinates): - raise ValueError("incident_positions have to be None or Coordinates") - - sph = incident_positions.get_sph() - theta = sph[..., 1] - weight = np.sin(2*theta) # sin(2*theta) - norm = np.sum(weight) - random_scattering = scattering_coefficient*weight/norm - random_scattering.freq = np.sum(random_scattering.freq, axis=-2) - random_scattering.freq = np.sum(random_scattering.freq, axis=-2) - random_scattering.comment = 'random-incidence scattering coefficient' - return random_scattering From 93901db8c801c49a84e157b372053842b822c6fc Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Thu, 23 Feb 2023 17:17:48 +0100 Subject: [PATCH 38/39] weights can be manipulated now --- imkar/integrate/integrate.py | 30 ++++++++++++++++++------------ tests/test_integrate.py | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 12 deletions(-) diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py index d04572f..710a9b8 100644 --- a/imkar/integrate/integrate.py +++ b/imkar/integrate/integrate.py @@ -4,7 +4,7 @@ import pyfar as pf -def surface_sphere(data, coords): +def surface_sphere(data, coords, weights=None): r"""Integrate over a set of points sampled on a spherical surface. .. math:: @@ -17,8 +17,12 @@ def surface_sphere(data, coords): Input data to integrate. Its `cshape` needs to be (..., #theta, #phi) or (..., #phi, #theta). coords : Coordinates - Coordinate points at which the data is sampled. Its cshape needs to be - (#theta, #phi) or (#phi, #theta), matching the cshape of `data`. + Coordinate points at which the data is sampled. Its `cshape` needs to + be (#theta, #phi) or (#phi, #theta), matching the `cshape` of `data`. + weights : np.array + weights for the integration over the coordinates. Its `cshape` + needs to be same as for `coords`. + By default the its sin(theta). Returns ------- @@ -51,27 +55,29 @@ def surface_sphere(data, coords): coords_spherical = coords.get_sph(convention='top_colat') phi = coords_spherical[1, :, 0] theta = coords_spherical[:, 1, 1] - axis_index = -2 + axis_index = -1 if np.sum(np.diff(phi[1:-2])) < 1e-3: phi = coords_spherical[:, 1, 0] theta = coords_spherical[1, :, 1] - axis_index = -3 + axis_index = -2 radius = coords_spherical[:, :, 2] r_mean = np.mean(radius) if not np.allclose(radius, r_mean): warnings.warn(r'all radii should be almost same') - # pf.Coordinates turns phi = 2*pi to 0, due to cylindircal + # pf.Coordinates turns phi = 2*pi to 0, due to circular # for the integration the upper limit cannot be zero but 2pi last_phi_is_zero = np.isclose(phi[-1], 0, atol=1e-16) increasing_phi_to_2pi = np.isclose( phi[-2]+np.diff(phi)[0], 2*np.pi, atol=1e-16) if last_phi_is_zero and increasing_phi_to_2pi: - phi[-1] = 2*np.pi + phi[-1] = 2 * np.pi - # integrate over angles with wight for sperical integration - weights = np.transpose(np.atleast_2d(np.sin(theta))) - result_raw = trapezoid(data.freq, x=phi, axis=axis_index) - result_raw1 = trapezoid(result_raw*weights, x=theta, axis=-2) + # integrate over angles with wight for spherical integration + data_int = np.moveaxis(data.freq, -1, 0) + if weights is None: + weights = np.sin(coords_spherical[..., 1]) + result_raw = trapezoid(data_int*weights, x=phi, axis=axis_index) + result_raw1 = trapezoid(result_raw, x=theta, axis=-1) - return pf.FrequencyData(result_raw1, data.frequencies) + return pf.FrequencyData(np.moveaxis(result_raw1, 0, -1), data.frequencies) diff --git a/tests/test_integrate.py b/tests/test_integrate.py index 5dbe6e8..3abc702 100644 --- a/tests/test_integrate.py +++ b/tests/test_integrate.py @@ -98,3 +98,35 @@ def test_surface_sphere_nonuniform_data_different_limits( theta_lower = -np.cos(2*theta_0)/4 desired = (theta_upper - theta_lower) * (phi_1 - phi_0) np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) + + +@pytest.mark.parametrize( + " phi_0, phi_1, theta_0, theta_1", [ + (0, 2*np.pi, 0, np.pi), + (0, np.pi, 0, np.pi), + (0, np.pi, 0, np.pi/2), + (np.pi, 2*np.pi, 0, np.pi/2), + (np.pi-.1, 2*np.pi-.1, 0, np.pi/2), + (0, 2*np.pi, np.pi/2, np.pi), + (np.pi-.1, 2*np.pi-.1, 0, np.pi), + (np.pi, 2*np.pi, 0, np.pi), + ]) +def test_surface_sphere_nonuniform_data_different_limits_over_weights( + phi_0, phi_1, theta_0, theta_1): + # same as test_surface_sphere_nonuniform_data_different_limits, but here + # we set the cos(theta) over the weights instead of the data + delta = np.pi/18 + coords = stub_utils.spherical_coordinates( + np.arange(phi_0, phi_1+delta, delta), + np.arange(theta_0, theta_1+delta, delta)) + data = stub_utils.frequencydata_from_shape(coords.cshape, 1, 100) + theta = coords.get_sph()[..., 1] + # data.freq = data.freq[0, ...] + data.freq[..., 0] = 1 + weights = np.cos(theta)*np.sin(theta) + result = integrate.surface_sphere(data, coords, weights=weights) + actual = np.real(result.freq[0, 0]) + theta_upper = -np.cos(2*theta_1)/4 + theta_lower = -np.cos(2*theta_0)/4 + desired = (theta_upper - theta_lower) * (phi_1 - phi_0) + np.testing.assert_allclose(actual, desired, rtol=5e-3, atol=0.04) From 9770ad090196b73f3202a187784470d3f9f9e995 Mon Sep 17 00:00:00 2001 From: Anne Heimes <64446926+ahms5@users.noreply.github.com> Date: Thu, 23 Feb 2023 18:48:11 +0100 Subject: [PATCH 39/39] implement and test diffusion calculation --- imkar/diffusion/coefficient.py | 41 +++++++++++++++++------ tests/test_diffusion_coefficient.py | 52 +++++++++++++++++++++++++++++ 2 files changed, 82 insertions(+), 11 deletions(-) create mode 100644 tests/test_diffusion_coefficient.py diff --git a/imkar/diffusion/coefficient.py b/imkar/diffusion/coefficient.py index 751a83d..37a232d 100644 --- a/imkar/diffusion/coefficient.py +++ b/imkar/diffusion/coefficient.py @@ -1,12 +1,11 @@ import numpy as np import pyfar as pf -from imkar import integrate def freefield(sample_pressure, mic_positions): """ Calculate the free-field diffusion coefficient for each incident direction - using the Mommertz correlation method [1]_. See :py:func:`random_incidence` + after ISO 17497-2:2012 [1]_. See :py:func:`random_incidence` to calculate the random incidence diffusion coefficient. @@ -43,7 +42,7 @@ def freefield(sample_pressure, mic_positions): >>> import imkar as ik >>> diffusion_coefficients = ik.diffusion.coefficient.freefield( >>> sample_pressure, mic_positions) - >>> random_d = ik.diffusion.coefficient.random_incidence( + >>> random_d = ik.scattering.coefficient.random_incidence( >>> diffusion_coefficients, incident_positions) """ # check inputs @@ -52,16 +51,36 @@ def freefield(sample_pressure, mic_positions): if not isinstance(mic_positions, pf.Coordinates): raise ValueError("microphone positions have to be Coordinates") + # parse angles + coords_spherical = mic_positions.get_sph(convention='top_colat') + phi = coords_spherical[1, :, 0] + theta = coords_spherical[:, 1, 1] + if np.sum(np.diff(phi[1:-2])) < 1e-3: + phi = coords_spherical[:, 1, 0] + theta = coords_spherical[1, :, 1] + d_theta = np.mean(np.diff(theta)) + d_phi = np.mean(np.diff(phi)) + + # calculate weights + thetas = coords_spherical[..., 1] + A_i = 2 * np.sin(thetas) * np.sin(d_theta/2) + A_i[thetas == 0] = 2 * np.sin(d_theta/2) + A_i[thetas == np.pi/2] = 4 * np.pi / d_phi * (np.sin(d_theta/4)**2) + N_i = A_i / np.min(A_i) + # calculate according to Mommertz correlation method Equation (6) - p_sample_sq = pf.FrequencyData( - np.abs(sample_pressure.freq)**2, sample_pressure.frequencies) + p_sample_abs_sq = np.moveaxis(np.abs(sample_pressure.freq)**2, -1, 0) - p_sample_sum_sq = integrate.surface_sphere( - p_sample_sq*p_sample_sq, mic_positions) - p_sample_sq_sum = integrate.surface_sphere(p_sample_sq, mic_positions)**2 - n = 1 - diffusion_coefficients \ - = (p_sample_sum_sq - p_sample_sq_sum) / ((n-1) * p_sample_sq_sum) + p_sample_sum_sq = np.sum(np.sum( + p_sample_abs_sq**2 * N_i, axis=-1), axis=-1) + p_sample_sq_sum = np.sum(np.sum( + p_sample_abs_sq * N_i, axis=-1), axis=-1)**2 + n = np.sum(N_i) + diffusion_array \ + = (p_sample_sq_sum - p_sample_sum_sq) / ((n-1) * p_sample_sum_sq) + diffusion_coefficients = pf.FrequencyData( + np.moveaxis(diffusion_array, 0, -1), + sample_pressure.frequencies) diffusion_coefficients.comment = 'diffusion coefficients' return diffusion_coefficients diff --git a/tests/test_diffusion_coefficient.py b/tests/test_diffusion_coefficient.py new file mode 100644 index 0000000..f9d78f8 --- /dev/null +++ b/tests/test_diffusion_coefficient.py @@ -0,0 +1,52 @@ +import pytest +import numpy as np + +from imkar import diffusion +from imkar.testing import stub_utils + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield(frequencies): + + mics = stub_utils.spherical_coordinates( + np.arange(0, 360, 10)/180*np.pi, + np.arange(10, 90, 10)/180*np.pi, + 1) + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, 1, frequencies) + d = diffusion.coefficient.freefield(p_sample, mics) + np.testing.assert_allclose(d.freq, 1) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +@pytest.mark.parametrize("radius", [ + (1), (10)]) +@pytest.mark.parametrize("magnitude", [ + (1), (10), (np.ones((5, 5)))]) +def test_freefield_with_theta_0(frequencies, radius, magnitude): + mics = stub_utils.spherical_coordinates( + np.arange(0, 360, 10)/180*np.pi, + np.arange(0, 100, 10)/180*np.pi, + radius) + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, magnitude, frequencies) + d = diffusion.coefficient.freefield(p_sample, mics) + np.testing.assert_allclose(d.freq, 1) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +@pytest.mark.parametrize("radius", [ + (1), (10)]) +def test_freefield_not_one(frequencies, radius): + mics = stub_utils.spherical_coordinates( + np.arange(0, 360, 10)/180*np.pi, + np.arange(0, 100, 10)/180*np.pi, + radius) + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, 1, frequencies) + p_sample.freq[1, 0, :] = 2 + d = diffusion.coefficient.freefield(p_sample, mics) + np.testing.assert_allclose(d.freq, 0.9918, atol=0.003)