diff --git a/docs/modules.rst b/docs/modules.rst
index 52a15ec..9da889f 100644
--- a/docs/modules.rst
+++ b/docs/modules.rst
@@ -7,4 +7,4 @@ according to their modules.
 .. toctree::
    :maxdepth: 1
 
-
+   modules/imkar.scattering.coefficient
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..e176307 100644
--- a/imkar/__init__.py
+++ b/imkar/__init__.py
@@ -3,3 +3,6 @@
 __author__ = """The pyfar developers"""
 __email__ = 'info@pyfar.org'
 __version__ = '0.1.0'
+
+
+from . import scattering  # noqa
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..6dbc11f
--- /dev/null
+++ b/imkar/scattering/coefficient.py
@@ -0,0 +1,127 @@
+import numpy as np
+import pyfar as pf
+
+
+def freefield(sample_pressure, reference_pressure, weights_microphones):
+    """
+    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 : pf.FrequencyData
+        Reflected sound pressure or directivity of the test sample. Its cshape
+        need to be (..., #microphones).
+    reference_pressure : pf.FrequencyData
+        Reflected sound pressure or directivity of the test
+        reference sample. It has the same shape as `sample_pressure`.
+    weights_microphones : np.ndarray
+        An array object with all weights for the microphone positions.
+        Its cshape need to be (#microphones). Microphone positions need to be
+        same for `sample_pressure` and `reference_pressure`.
+
+    Returns
+    -------
+    scattering_coefficients : FrequencyData
+        The scattering coefficient for each incident 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, June 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.weights)
+    >>> 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(weights_microphones, np.ndarray):
+        raise ValueError("weights_microphones have to be a numpy.array")
+    if not sample_pressure.cshape == reference_pressure.cshape:
+        raise ValueError(
+            "sample_pressure and reference_pressure have to have the "
+            "same cshape.")
+    if not weights_microphones.shape[0] == sample_pressure.cshape[-1]:
+        raise ValueError(
+            "the last dimension of sample_pressure need be same as the "
+            "weights_microphones.shape.")
+    if not np.allclose(
+            sample_pressure.frequencies, reference_pressure.frequencies):
+        raise ValueError(
+            "sample_pressure and reference_pressure have to have the "
+            "same frequencies.")
+
+    # calculate according to mommertz correlation method Equation (5)
+    p_sample = np.moveaxis(sample_pressure.freq, -1, 0)
+    p_reference = np.moveaxis(reference_pressure.freq, -1, 0)
+    p_sample_abs = np.abs(p_sample)
+    p_reference_abs = np.abs(p_reference)
+    p_sample_sq = p_sample_abs*p_sample_abs
+    p_reference_sq = p_reference_abs*p_reference_abs
+    p_cross = p_sample * np.conj(p_reference)
+
+    p_sample_sum = np.sum(p_sample_sq * weights_microphones, axis=-1)
+    p_ref_sum = np.sum(p_reference_sq * weights_microphones, axis=-1)
+    p_cross_sum = np.sum(p_cross * weights_microphones, axis=-1)
+
+    data_scattering_coefficient \
+        = 1 - ((np.abs(p_cross_sum)**2)/(p_sample_sum*p_ref_sum))
+
+    scattering_coefficients = pf.FrequencyData(
+        np.moveaxis(data_scattering_coefficient, 0, -1),
+        sample_pressure.frequencies)
+    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 : pf.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")
+
+    theta = incident_positions.get_sph().T[1]
+    weight = np.cos(theta) * incident_positions.weights
+    norm = np.sum(weight)
+    random_scattering = scattering_coefficient*weight/norm
+    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..dc61984
--- /dev/null
+++ b/imkar/testing/__init__.py
@@ -0,0 +1,7 @@
+from .stub_utils import (
+    frequency_data_from_shape,
+)
+
+__all__ = [
+    'frequency_data_from_shape',
+    ]
diff --git a/imkar/testing/stub_utils.py b/imkar/testing/stub_utils.py
new file mode 100644
index 0000000..3375398
--- /dev/null
+++ b/imkar/testing/stub_utils.py
@@ -0,0 +1,22 @@
+"""
+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 frequency_data_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/tests/test_imkar.py b/tests/test_imkar.py
index f801097..6b6c622 100644
--- a/tests/test_imkar.py
+++ b/tests/test_imkar.py
@@ -5,7 +5,7 @@
 import pytest
 
 
-from imkar import imkar
+from imkar import imkar  # noqa
 
 
 @pytest.fixture
@@ -18,7 +18,7 @@ def response():
     # return requests.get('https://github.com/mberz/cookiecutter-pypackage')
 
 
-def test_content(response):
+def test_content(response):  # noqa
     """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_scattering_coefficient.py b/tests/test_scattering_coefficient.py
new file mode 100644
index 0000000..22c9717
--- /dev/null
+++ b/tests/test_scattering_coefficient.py
@@ -0,0 +1,198 @@
+import pytest
+import numpy as np
+import pyfar as pf
+
+from imkar import scattering
+from imkar.testing import stub_utils
+
+
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_freefield_1(frequencies):
+    mics = pf.samplings.sph_gaussian(42)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    p_sample = stub_utils.frequency_data_from_shape(
+        mics.cshape, 1, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        mics.cshape, 0, frequencies)
+    p_sample.freq[5, :] = 0
+    p_reference.freq[5, :] = np.sum(p_sample.freq.flatten())/2
+    s = scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+    np.testing.assert_allclose(s.freq, 1)
+
+
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_freefield_wrong_input(frequencies):
+    mics = pf.samplings.sph_gaussian(42)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    p_sample = stub_utils.frequency_data_from_shape(
+        mics.cshape, 1, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        mics.cshape, 0, frequencies)
+    p_sample.freq[5, :] = 0
+    p_reference.freq[5, :] = np.sum(p_sample.freq.flatten())/2
+    with pytest.raises(ValueError, match='sample_pressure'):
+        scattering.coefficient.freefield(1, p_reference, mics.weights)
+    with pytest.raises(ValueError, match='reference_pressure'):
+        scattering.coefficient.freefield(p_sample, 1, mics.weights)
+    with pytest.raises(ValueError, match='weights_microphones'):
+        scattering.coefficient.freefield(p_sample, p_reference, 1)
+    with pytest.raises(ValueError, match='cshape'):
+        scattering.coefficient.freefield(
+            p_sample[:-2, ...], p_reference, mics.weights)
+    with pytest.raises(ValueError, match='weights_microphones'):
+        scattering.coefficient.freefield(
+            p_sample, p_reference, mics.weights[:10])
+    with pytest.raises(ValueError, match='same frequencies'):
+        p_sample.frequencies[0] = 1
+        scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+
+
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_freefield_05(frequencies):
+    mics = pf.samplings.sph_gaussian(42)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    p_sample = stub_utils.frequency_data_from_shape(
+        mics.cshape, 0, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        mics.cshape, 0, frequencies)
+    p_sample.freq[7, :] = 1
+    p_sample.freq[28, :] = 1
+    p_reference.freq[7, :] = 1
+    s = scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+    np.testing.assert_allclose(s.freq, 0.5)
+
+
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_freefield_0(frequencies):
+    mics = pf.samplings.sph_gaussian(42)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    p_sample = stub_utils.frequency_data_from_shape(
+        mics.cshape, 0, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        mics.cshape, 0, frequencies)
+    p_reference.freq[5, :] = 1
+    p_sample.freq[5, :] = 1
+    s = scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+    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_incident(frequencies):
+    mics = pf.samplings.sph_equal_area(42)
+    mics.weights = pf.samplings.calculate_sph_voronoi_weights(mics)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    incident_directions = pf.samplings.sph_gaussian(10)
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[1] <= np.pi/2]
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[0] <= np.pi/2]
+    data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten()
+    p_sample = stub_utils.frequency_data_from_shape(
+        data_shape, 0, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        data_shape, 0, frequencies)
+    p_reference.freq[:, 2, :] = 1
+    p_sample.freq[:, 2, :] = 1
+    s = scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+    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_incidence(
+        frequencies):
+    mics = pf.samplings.sph_equal_angle(10)
+    mics.weights = pf.samplings.calculate_sph_voronoi_weights(mics)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    incident_directions = pf.samplings.sph_gaussian(10)
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[1] <= np.pi/2]
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[0] <= np.pi/2]
+    data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten()
+    p_sample = stub_utils.frequency_data_from_shape(
+        data_shape, 0, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        data_shape, 0, frequencies)
+    p_reference.freq[:, 2, :] = 1
+    p_sample.freq[:, 3, :] = 1
+    s = scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+    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])])
+@pytest.mark.parametrize("mics",  [
+    (pf.samplings.sph_equal_angle(10)),
+    ])
+def test_freefield_05_with_inci(frequencies, mics):
+    mics.weights = pf.samplings.calculate_sph_voronoi_weights(mics)
+    mics = mics[mics.get_sph().T[1] <= np.pi/2]  # delete lower part of sphere
+    incident_directions = pf.samplings.sph_gaussian(10)
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[1] <= np.pi/2]
+    incident_directions = incident_directions[
+        incident_directions.get_sph().T[0] <= np.pi/2]
+    data_shape = np.array((incident_directions.cshape, mics.cshape)).flatten()
+    p_sample = stub_utils.frequency_data_from_shape(
+        data_shape, 0, frequencies)
+    p_reference = stub_utils.frequency_data_from_shape(
+        data_shape, 0, frequencies)
+    p_sample.freq[:, 7, :] = 1
+    p_sample.freq[:, 5, :] = 1
+    p_reference.freq[:, 5, :] = 1
+    s = scattering.coefficient.freefield(p_sample, p_reference, mics.weights)
+    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 = pf.samplings.sph_gaussian(10)
+    incident_directions = coords[coords.get_sph().T[1] <= np.pi/2]
+    s = stub_utils.frequency_data_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):
+    data = pf.samplings.sph_gaussian(10)
+    incident_directions = data[data.get_sph().T[1] <= np.pi/2]
+    incident_cshape = incident_directions.cshape
+    s_value = np.arange(
+        incident_cshape[0]).reshape(incident_cshape) / incident_cshape[0]
+    theta = incident_directions.get_sph().T[1]
+    actual_weight = np.cos(theta) * incident_directions.weights
+    actual_weight /= np.sum(actual_weight)
+    s = stub_utils.frequency_data_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..99a9003
--- /dev/null
+++ b/tests/test_sub_utils.py
@@ -0,0 +1,37 @@
+import pytest
+import numpy as np
+
+from imkar.testing import stub_utils
+
+
+@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_frequency_data_from_shape(shapes, data_in, frequency):
+    data = stub_utils.frequency_data_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)