diffusion coefficient

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,8 @@
 __author__ = """The pyfar developers"""
 __email__ = '[email protected]'
 __version__ = '0.1.0'
+
+
+from . import integrate  # noqa
+from . import scattering  # noqa
+from . import diffusion  # noqa

imkar/diffusion/__init__.py

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

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

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