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,110 @@
+"""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 scattering coefficient from free-field reflection pattern.
+
+    This function uses the Mommertz correlation method [#]_ to compute the
+    scattering coefficient by averaging over the microphone positions:
+
+    .. math::
+        s = 1 -
+            \frac{|\sum_{\Omega_r} \underline{p}_{\text{sample}}(\Omega_r)
+            \cdot \underline{p}_{\text{reference}}^*(\Omega_r)
+            \cdot w(\Omega_r)|^2}
+            {\sum_{\Omega_r} |\underline{p}_{\text{sample}}(\Omega_r)|^2
+            \cdot w(\Omega_r) \cdot \sum_{\Omega_r}
+            |\underline{p}_{\text{reference}}(\Omega_r)|^2
+            \cdot w(\Omega_r) }
+
+    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:`\Omega_r` represents the solid angle of the microphone
+          positions and :math:`w(\Omega_r)` represents its area weight.
+
+    Parameters
+    ----------
+    sample_pressure : pyfar.FrequencyData, pyfar.Signal
+        Reflected sound pressure or directivity of the test sample. Its
+        `cshape` must be ``(..., microphone_weights.size)`` and broadcastable
+        to the `cshape` of ``reference_pressure``. The frequency vectors
+        of both ``sample_pressure`` and ``reference_pressure`` must match.
+    reference_pressure : pyfar.FrequencyData, pyfar.Signal
+        Reflected sound pressure or directivity of the reference sample. Its
+        `cshape` must be (..., microphone_weights.size) and broadcastable
+        to the `cshape` of ``sample_pressure``. The frequency vectors of both
+        ``sample_pressure`` and ``reference_pressure`` must match.
+    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 : pyfar.FrequencyData
+        The scattering coefficient of the broadcasted `cshape` of
+        ``sample_pressure`` and ``reference_pressure``, excluding the
+        last dimension.
+
+    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, pf.Signal)):
+        raise TypeError(
+            "sample_pressure must be of type pyfar.FrequencyData or "
+            "pyfar.Signal")
+    if not isinstance(reference_pressure, (pf.FrequencyData, pf.Signal)):
+        raise TypeError(
+            "reference_pressure must be of type pyfar.FrequencyData or "
+            "pyfar.Signal")
+    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")
+    # check frequency vectors
+    if not np.allclose(
+            reference_pressure.frequencies, sample_pressure.frequencies,
+            atol=1e-15):
+        raise ValueError("The frequencies do not match.")
+
+    # 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,168 @@
+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):
+    """Generate a plane wave for a given direction and sampling.
+
+    This function is just used for testing purposes, so the frequency is
+    fixed to 5000 Hz.
+
+    Parameters
+    ----------
+    amplitude : float
+        The amplitude of the plane wave.
+    direction : pf.Coordinates
+        The direction of the plane wave.
+    sampling : pf.Sampling
+        The sampling grid for the plane wave.
+
+    Returns
+    -------
+    pf.FrequencyData
+        The generated plane wave in the frequency domain.
+    """
+
+    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_zero_scattering():
+    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_fractional_scattering(s_scatter, Phi_scatter_deg):
+    """
+    Test analytic scattering coefficient calculation for non 0 or 1.
+
+    For the sample pressure, two plane waves are used. The scattered and the
+    specular wave. The reflection factor of both waves is calculated based on
+    the scattering coefficient and the scattering angle.
+    """
+    # calculate the specular and scattered reflection factors
+    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)))
+
+    # define the sampling grid and generate the plane waves
+    # for the sample pressure and the reference pressure
+    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)
+
+    # Calculate the scattering coefficient with for the specular wave
+    # in other words, the reference pressure is the specular wave (R=1)
+    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,
+    )
+    npt.assert_almost_equal(sd_spec.freq, s_spec, 1)
+
+    # Calculate the scattering coefficient for the scattered wave
+    # in other words, the reference pressure is the scattered wave (R=1)
+    # so the result will be 1-s_scatter
+    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_scatter.freq, s_scatter, 1)
+
+    # Calculate the scattering coefficient for 5° more than the specular
+    # wave, then the scattering coefficient should be 0
+    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_one_scattering():
+    """Two plane waves in the different directions should result in s=1."""
+    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, 1, 0), 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, 1, decimal=3)
+
+
+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]]).T, [100])
+    reference_pressure = pf.FrequencyData(
+        np.array([[1, 2]]).T, [100])
+    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]]).T, [100])
+    reference_pressure = pf.FrequencyData(
+        np.array([[1, 2, 3]]).T, [100])
+    microphone_weights = np.array([0.5, 0.5])
+    with pytest.raises(
+            ValueError, match="The last dimension of reference_pressure must "
+            "match the size of microphone_weights"):
+        sff.correlation_method(
+            sample_pressure, reference_pressure, microphone_weights)