Feature: correlation method after mommertz

docs/api_reference.rst

@@ -16,7 +16,7 @@ Modules
 .. toctree::
    :maxdepth: 1
 
-   modules/imkar
+   modules/imkar.scattering.freefield
 
 
 .. _examples gallery: https://pyfar-gallery.readthedocs.io/en/latest/examples_gallery.html

docs/modules/imkar.scattering.freefield.rst

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

imkar/__init__.py

@@ -5,3 +5,10 @@
 __author__ = """The pyfar developers"""
 __email__ = ''
 __version__ = '0.1.0'
+
+
+from . import scattering
+
+__all__ = [
+    "scattering",
+]

imkar/scattering/__init__.py

@@ -0,0 +1,7 @@
+"""Imkar scattering module."""
+
+from . import freefield
+
+__all__ = [
+    "freefield",
+]

imkar/scattering/freefield.py

@@ -0,0 +1,111 @@
+"""Scattering calculation functions based on free-field data."""
+import numpy as np
+import pyfar as pf
+
+
+def correlation_method(
+        sample_pressure, reference_pressure, microphone_weights):
+    r"""
+    Calculate the incident-dependent free-field scattering coefficient.
+
+    This function uses the Mommertz correlation method [#]_ to compute the
+    scattering coefficient of the input data:
+
+    .. math::
+        s = 1 -
+            \frac{|\sum_w \underline{p}_{\text{sample}}(\vartheta,\varphi)
+            \cdot \underline{p}_{\text{reference}}^*(\vartheta,\varphi)
+            \cdot w(\vartheta,\varphi)|^2}
+            {\sum_w |\underline{p}_{\text{sample}}(\vartheta,\varphi)|^2
+            \cdot w(\vartheta,\varphi) \cdot \sum_w
+            |\underline{p}_{\text{reference}}(\vartheta,\varphi)|^2
+            \cdot w(\vartheta,\varphi) }
+
+    where:
+        - :math:`\underline{p}_{\text{sample}}` is the reflected sound
+          pressure of the sample under investigation.
+        - :math:`\underline{p}_{\text{reference}}` is the reflected sound
+          pressure from the reference sample.
+        - :math:`w` represents the area weights of the sampling, and
+          :math:`\vartheta` and :math:`\varphi` are the ``colatitude``
+          and ``azimuth`` angles from the
+          :py:class:`~pyfar.classes.coordinates.Coordinates` object.
+
+    The test sample is assumed to lie in the x-y-plane.
+
+    Parameters
+    ----------
+    sample_pressure : :py:class:`~pyfar.classes.audio.FrequencyData`
+        Reflected sound pressure or directivity of the test sample. Its cshape
+        must be (..., microphone_weights.size) and broadcastable to the
+        cshape of ``reference_pressure``.
+    reference_pressure : :py:class:`~pyfar.classes.audio.FrequencyData`
+        Reflected sound pressure or directivity of the reference sample. Its
+        cshape must be (..., microphone_weights.size) and broadcastable to the
+        cshape of ``sample_pressure``.
+    microphone_weights : array_like
+        1D array containing the area weights for the microphone positions.
+        No normalization is required. Its shape must match the last dimension
+        in the cshape of ``sample_pressure`` and ``reference_pressure``.
+
+    Returns
+    -------
+    scattering_coefficients : :py:class:`~pyfar.classes.audio.FrequencyData`
+        The scattering coefficient for each incident direction as a function
+        of frequency.
+
+    References
+    ----------
+    .. [#]  E. Mommertz, "Determination of scattering coefficients from the
+            reflection directivity of architectural surfaces," Applied
+            Acoustics, vol. 60, no. 2, pp. 201-203, June 2000,
+            doi: 10.1016/S0003-682X(99)00057-2.
+
+    """
+    # check input types
+    if not isinstance(sample_pressure, pf.FrequencyData):
+        raise TypeError("sample_pressure must be of type pyfar.FrequencyData")
+    if not isinstance(reference_pressure, pf.FrequencyData):
+        raise TypeError(
+            "reference_pressure must be of type pyfar.FrequencyData")
+    microphone_weights = np.atleast_1d(
+        np.asarray(microphone_weights, dtype=float))
+
+    # check input dimensions
+    if sample_pressure.cshape[-1] != microphone_weights.size:
+        raise ValueError(
+            "The last dimension of sample_pressure must match the size of "
+            "microphone_weights")
+    if reference_pressure.cshape[-1] != microphone_weights.size:
+        raise ValueError(
+            "The last dimension of reference_pressure must match the size of "
+            "microphone_weights")
+
+    if sample_pressure.cshape[:-1] != reference_pressure.cshape[:-1]:
+        raise ValueError(
+            "The cshape of sample_pressure and reference_pressure must be "
+            "broadcastable except for the last dimension")
+    # Test whether the objects are able to perform arithmetic operations.
+    # e.g. does the frequency vectors match
+    _ = sample_pressure + reference_pressure
+
+    # prepare data
+    microphone_weights = microphone_weights[:, np.newaxis]
+    p_sample = sample_pressure.freq
+    p_reference = reference_pressure.freq
+
+    # calculate according to mommertz correlation method Equation (5)
+    p_sample_sum = np.sum(microphone_weights * np.abs(p_sample)**2, axis=-2)
+    p_ref_sum = np.sum(microphone_weights * np.abs(p_reference)**2, axis=-2)
+    p_cross_sum = np.sum(
+            p_sample * np.conj(p_reference) * microphone_weights, axis=-2)
+
+    data_scattering_coefficient \
+        = 1 - ((np.abs(p_cross_sum)**2)/(p_sample_sum*p_ref_sum))
+
+    # create pyfar.FrequencyData object
+    scattering_coefficients = pf.FrequencyData(
+        data_scattering_coefficient,
+        sample_pressure.frequencies)
+
+    return scattering_coefficients

tests/test_scattering_freefield.py

@@ -0,0 +1,138 @@
+import pyfar as pf
+import numpy as np
+import numpy.testing as npt
+import pytest
+from imkar.scattering import freefield as sff
+
+
+def plane_wave(amplitude, direction, sampling):
+    f = 5000
+    c = 343
+    x = sampling
+    direction.cartesian = direction.cartesian/direction.radius
+    dot_product = direction.x*x.x+direction.y*x.y+direction.z*x.z
+    dot_product = dot_product[..., np.newaxis]
+    f = np.atleast_1d(f)
+    return pf.FrequencyData(
+        amplitude*np.exp(-1j*2*np.pi*f/c*dot_product),
+        frequencies=f,
+    )
+
+
+def test_correlation_method_0():
+    sampling = pf.samplings.sph_equal_area(5000)
+    sampling.weights = pf.samplings.calculate_sph_voronoi_weights(sampling)
+    sampling = sampling[sampling.z>0]
+    sample_pressure = plane_wave(1, pf.Coordinates(0, 0, 1), sampling)
+    reference_pressure = plane_wave(1, pf.Coordinates(0, 0, 1), sampling)
+    s = sff.correlation_method(
+        sample_pressure, reference_pressure, sampling.weights,
+    )
+    npt.assert_almost_equal(s.freq, 0)
+
+
+@pytest.mark.parametrize("s_scatter", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("Phi_scatter_deg", [30, 60, 90, 120, 150, 42])
+def test_correlation_method_with_plane_waves(s_scatter, Phi_scatter_deg):
+    s_spec = 1-s_scatter
+    Phi_spec = 45/180*np.pi
+    Phi_scatter = Phi_scatter_deg/180*np.pi
+    R_spec = np.sqrt(s_spec)
+    R_scatter = np.sqrt(np.abs(s_scatter*np.sin(Phi_spec)/np.sin(Phi_scatter)))
+    sampling = pf.samplings.sph_equal_area(10000)
+    sampling.weights = pf.samplings.calculate_sph_voronoi_weights(sampling)
+    sampling = sampling[sampling.z>0]
+    sample_pressure = plane_wave(
+        R_spec,
+        pf.Coordinates.from_spherical_front(np.pi/2, Phi_spec, 1), sampling)
+    sample_pressure += plane_wave(
+        R_scatter,
+        pf.Coordinates.from_spherical_front(np.pi/2, Phi_scatter, 1), sampling)
+    reference_pressure = plane_wave(
+        1, pf.Coordinates.from_spherical_front(np.pi/2, Phi_spec, 1), sampling)
+    sd_spec = 1-sff.correlation_method(
+        sample_pressure, reference_pressure, sampling.weights,
+    )
+    reference_pressure = plane_wave(
+        1, pf.Coordinates.from_spherical_front(
+            np.pi/2, Phi_scatter, 1), sampling)
+    sd_scatter = 1-sff.correlation_method(
+        sample_pressure, reference_pressure, sampling.weights,
+    )
+    npt.assert_almost_equal(sd_spec.freq, s_spec, 1)
+    npt.assert_almost_equal(sd_scatter.freq, s_scatter, 1)
+
+    reference_pressure = plane_wave(
+        1, pf.Coordinates.from_spherical_front(
+            np.pi/2, Phi_spec+5/180*np.pi, 1), sampling)
+    sd_scatter_0 = 1-sff.correlation_method(
+        sample_pressure, reference_pressure, sampling.weights,
+    )
+    npt.assert_almost_equal(sd_scatter_0.freq, 0, 1)
+
+
+def test_correlation_method():
+    sampling = pf.samplings.sph_equal_area(5000)
+    sampling.weights = pf.samplings.calculate_sph_voronoi_weights(sampling)
+    sampling = sampling[sampling.z>0]
+    sample_pressure = plane_wave(1, pf.Coordinates(0, 0, 1), sampling)
+    reference_pressure = plane_wave(1, pf.Coordinates(0, 0, 1), sampling)
+    s = sff.correlation_method(
+        sample_pressure, reference_pressure, sampling.weights,
+    )
+    npt.assert_almost_equal(s.freq, 0)
+
+
+def test_correlation_method_invalid_sample_pressure_type():
+    reference_pressure = pf.FrequencyData(np.array([1, 2, 3]), [100, 200, 300])
+    microphone_weights = np.array([0.5, 0.5, 0.5])
+    with pytest.raises(
+            TypeError, match="sample_pressure must be of type "
+            "pyfar.FrequencyData"):
+        sff.correlation_method(
+            "invalid_type", reference_pressure, microphone_weights)
+
+
+def test_correlation_method_invalid_reference_pressure_type():
+    sample_pressure = pf.FrequencyData(np.array([1, 2, 3]), [100, 200, 300])
+    microphone_weights = np.array([0.5, 0.5, 0.5])
+    with pytest.raises(
+            TypeError, match="reference_pressure must be of type "
+            "pyfar.FrequencyData"):
+        sff.correlation_method(
+            sample_pressure, "invalid_type", microphone_weights)
+
+
+def test_correlation_method_mismatched_sample_pressure_weights():
+    sample_pressure = pf.FrequencyData(np.array([[1, 2, 3]]), [100, 200, 300])
+    reference_pressure = pf.FrequencyData(
+        np.array([[1, 2, 3]]), [100, 200, 300])
+    microphone_weights = np.array([0.5, 0.5])
+    with pytest.raises(
+            ValueError, match="The last dimension of sample_pressure must "
+            "match the size of microphone_weights"):
+        sff.correlation_method(
+            sample_pressure, reference_pressure, microphone_weights)
+
+
+def test_correlation_method_mismatched_reference_pressure_weights():
+    sample_pressure = pf.FrequencyData(np.array([[1, 2, 3]]), [100, 200, 300])
+    reference_pressure = pf.FrequencyData(
+        np.array([[1, 2, 3]]), [100, 200, 300])
+    microphone_weights = np.array([0.5, 0.5])
+    with pytest.raises(
+            ValueError, match="The last dimension of sample_pressure must "
+            "match the size of microphone_weights"):
+        sff.correlation_method(
+            sample_pressure, reference_pressure, microphone_weights)
+
+
+def test_correlation_method_mismatched_sample_reference_shapes():
+    sample_pressure = pf.FrequencyData(np.array([[1, 2, 3]]), [100, 200, 300])
+    reference_pressure = pf.FrequencyData(np.array([[1, 2]]), [100, 200])
+    microphone_weights = np.array([0.5, 0.5, 0.5])
+    with pytest.raises(
+            ValueError, match="The last dimension of sample_pressure must "
+            "match the size of microphone_weights"):
+        sff.correlation_method(
+            sample_pressure, reference_pressure, microphone_weights)