diff --git a/docs/conf.py b/docs/conf.py index 12123b8..b6932e3 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'] @@ -157,6 +170,3 @@ 'One line description of project.', 'Miscellaneous'), ] - - - diff --git a/docs/modules.rst b/docs/modules.rst new file mode 100644 index 0000000..3d9cc24 --- /dev/null +++ b/docs/modules.rst @@ -0,0 +1,11 @@ +Modules +======= + +The following gives detailed information about all pyfar functions sorted +according to their modules. + +.. toctree:: + :maxdepth: 1 + + modules/imkar.integrate + modules/imkar.scattering.coefficient diff --git a/docs/modules/imkar.integrate.rst b/docs/modules/imkar.integrate.rst new file mode 100644 index 0000000..8a6b0c6 --- /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/docs/modules/imkar.scattering.coefficient.rst b/docs/modules/imkar.scattering.coefficient.rst new file mode 100644 index 0000000..f384e68 --- /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/__init__.py b/imkar/__init__.py index 2d1e55d..8756e6c 100644 --- a/imkar/__init__.py +++ b/imkar/__init__.py @@ -3,3 +3,8 @@ __author__ = """The pyfar developers""" __email__ = 'info@pyfar.org' __version__ = '0.1.0' + + +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..37a232d --- /dev/null +++ b/imkar/diffusion/coefficient.py @@ -0,0 +1,86 @@ +import numpy as np +import pyfar as pf + + +def freefield(sample_pressure, mic_positions): + """ + Calculate the free-field diffusion coefficient for each incident direction + after ISO 17497-2:2012 [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 free-field 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.scattering.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") + + # 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_abs_sq = np.moveaxis(np.abs(sample_pressure.freq)**2, -1, 0) + + 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/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..f8ec396 --- /dev/null +++ b/imkar/integrate/__init__.py @@ -0,0 +1,7 @@ +from .integrate import ( + surface_sphere + ) + + +__all__ = [ + 'surface_sphere'] diff --git a/imkar/integrate/integrate.py b/imkar/integrate/integrate.py new file mode 100644 index 0000000..710a9b8 --- /dev/null +++ b/imkar/integrate/integrate.py @@ -0,0 +1,83 @@ +import warnings +from scipy.integrate import trapezoid +import numpy as np +import pyfar as pf + + +def surface_sphere(data, coords, weights=None): + r"""Integrate over a set of points sampled on a spherical surface. + + .. math:: + + S = \int \int data(\phi, \theta) sin(\theta) d\phi d\theta + + Parameters + ---------- + data : FrequencyData, Signal + 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`. + 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 + ------- + 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( + 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] + 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 = -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 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 + + # 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(np.moveaxis(result_raw1, 0, -1), data.frequencies) diff --git a/imkar/scattering/__init__.py b/imkar/scattering/__init__.py new file mode 100644 index 0000000..5e3fc4d --- /dev/null +++ b/imkar/scattering/__init__.py @@ -0,0 +1,9 @@ +from .coefficient import ( + random_incidence, + freefield +) + +__all__ = [ + 'freefield', + 'random_incidence' + ] diff --git a/imkar/scattering/coefficient.py b/imkar/scattering/coefficient.py new file mode 100644 index 0000000..6638027 --- /dev/null +++ b/imkar/scattering/coefficient.py @@ -0,0 +1,126 @@ +import numpy as np +import pyfar as pf +from imkar import integrate + + +def freefield(sample_pressure, reference_pressure, mic_positions): + """ + Calculate the free-field scattering coefficient for each incident direction + using the Mommertz correlation method [1]_. See :py:func:`random_incidence` + to calculate the random incidence scattering coefficient. + + + Parameters + ---------- + 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. + reference_pressure : FrequencyData + 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. Microphone positions need to be same for + `sample_pressure` and `reference_pressure`. + + Returns + ------- + scattering_coefficients : FrequencyData + The scattering coefficient for each plane wave direction. + + + References + ---------- + .. [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. + + Examples + -------- + Calculate freefield scattering coefficients and then the random incidence + scattering coefficient. + + >>> 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(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") + 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( + 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. 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 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..dc54762 --- /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 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') + 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/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 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..9a348d0 --- /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.spherical_coordinates( + 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.spherical_coordinates( + 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.spherical_coordinates( + 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.spherical_coordinates( + 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.spherical_coordinates( + 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.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_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) 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_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 diff --git a/tests/test_integrate.py b/tests/test_integrate.py new file mode 100644 index 0000000..3abc702 --- /dev/null +++ b/tests/test_integrate.py @@ -0,0 +1,132 @@ +import pytest +import numpy as np + +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), + (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) + 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) + + +@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) + 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) + result = integrate.surface_sphere(data, coords) + actual = np.squeeze(np.real(result.freq)) + np.testing.assert_allclose(actual, data_raw*desired, rtol=5e-3) + + +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'): + integrate.surface_sphere(data, coords) + + +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'): + 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) + + +@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) diff --git a/tests/test_scattering_coefficient.py b/tests/test_scattering_coefficient.py new file mode 100644 index 0000000..8b28dcd --- /dev/null +++ b/tests/test_scattering_coefficient.py @@ -0,0 +1,168 @@ +import pytest +import numpy as np + +from imkar import scattering +from imkar.testing import stub_utils + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_1(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, 1, frequencies) + p_reference = stub_utils.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) + np.testing.assert_allclose(s.freq, 1) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +def test_freefield_wrong_input(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, 1, frequencies) + p_reference = stub_utils.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='sample_pressure'): + scattering.coefficient.freefield(1, p_reference, mics) + 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) + 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", [ + ([100, 200]), ([100])]) +def test_freefield_05(frequencies, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, 0, frequencies) + p_reference = stub_utils.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 + 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, coords_half_sphere_10_deg): + mics = coords_half_sphere_10_deg + p_sample = stub_utils.frequencydata_from_shape( + mics.cshape, 0, frequencies) + p_reference = stub_utils.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) + 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, 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 = stub_utils.frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference = stub_utils.frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference.freq[:, :, 1, 2, :] = 1 + p_sample.freq[:, :, 1, 2, :] = 1 + 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) + 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, 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 = stub_utils.frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference = stub_utils.frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference.freq[:, :, 1, 2, :] = 1 + p_sample.freq[:, :, 2, 3, :] = 1 + 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) + assert s.cshape == incident_directions.cshape + assert s_rand.freq.shape[-1] == len(frequencies) + + +@pytest.mark.parametrize("frequencies", [ + ([100, 200]), ([100])]) +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 = stub_utils.frequencydata_from_shape( + data_shape, 0, frequencies) + p_reference = stub_utils.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 + 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) + assert s.cshape == incident_directions.cshape + assert s_rand.freq.shape[-1] == len(frequencies) + + +@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, coords_half_sphere_10_deg): + incident_directions = coords_half_sphere_10_deg + s = stub_utils.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, 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 = stub_utils.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) 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)