Scattering-freefield

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'),
 ]
-
-
-

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

docs/modules/imkar.integrate.rst

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

docs/modules/imkar.scattering.coefficient.rst

@@ -0,0 +1,7 @@
+imkar.scattering.coefficient
+============================
+
+.. automodule:: imkar.scattering.coefficient
+   :members:
+   :undoc-members:
+   :show-inheritance:

imkar/__init__.py

@@ -3,3 +3,7 @@
 __author__ = """The pyfar developers"""
 __email__ = 'info@pyfar.org'
 __version__ = '0.1.0'
+
+
+from . import integrate  # noqa
+from . import scattering  # noqa

imkar/integrate/__init__.py

@@ -0,0 +1,7 @@
+from .integrate import (
+    surface_sphere
+    )
+
+
+__all__ = [
+    'surface_sphere']

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)

imkar/scattering/__init__.py

@@ -0,0 +1,9 @@
+from .coefficient import (
+    random_incidence,
+    freefield
+)
+
+__all__ = [
+    'freefield',
+    'random_incidence'
+    ]

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

imkar/testing/__init__.py

@@ -0,0 +1,5 @@
+"""Testing sub-package for pyfar."""
+
+from . import (stub_utils)
+
+__all__ = ['stub_utils']

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

requirements_dev.txt

@@ -11,3 +11,5 @@ pytest
 numpy
 scipy
 pyfar
+insipid-sphinx-theme
+autodocsumm

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))