diff --git a/docs/api_reference.rst b/docs/api_reference.rst
index 8f0383c..5556bd2 100644
--- a/docs/api_reference.rst
+++ b/docs/api_reference.rst
@@ -16,7 +16,7 @@ Modules
 .. toctree::
    :maxdepth: 1
 
-   modules/imkar
+   modules/imkar.utils
 
 
 .. _examples gallery: https://pyfar-gallery.readthedocs.io/en/latest/examples_gallery.html
diff --git a/docs/modules/imkar.rst b/docs/modules/imkar.utils.rst
similarity index 51%
rename from docs/modules/imkar.rst
rename to docs/modules/imkar.utils.rst
index 0ea5287..89d6335 100644
--- a/docs/modules/imkar.rst
+++ b/docs/modules/imkar.utils.rst
@@ -1,7 +1,7 @@
-imkar
-=====
+imkar.utils
+===========
 
-.. automodule:: imkar
+.. automodule:: imkar.utils
    :members:
    :undoc-members:
    :show-inheritance:
diff --git a/imkar/__init__.py b/imkar/__init__.py
index 0dcfbd5..7d314cf 100644
--- a/imkar/__init__.py
+++ b/imkar/__init__.py
@@ -5,3 +5,10 @@
 __author__ = """The pyfar developers"""
 __email__ = ''
 __version__ = '0.1.0'
+
+
+from . import utils
+
+__all__ = [
+    'utils',
+]
diff --git a/imkar/utils.py b/imkar/utils.py
new file mode 100644
index 0000000..0a7d526
--- /dev/null
+++ b/imkar/utils.py
@@ -0,0 +1,62 @@
+"""Utilities for imkar."""
+import numpy as np
+import pyfar as pf
+
+
+def paris_formula(coefficients, colatitude_rad, area_weights):
+    r"""
+    Calculate the random-incidence coefficient
+    according to the Paris formula.
+
+    The Paris formula computes the random-incidence coefficient based on
+    directional coefficients. It is valid for scattering and absorption
+    coefficients and requires equally distributed incident directions
+    to get a valid result.
+    
+    The implementation follows the Equation 2.53 from [#]_ and is
+    discretized as:
+
+    .. math::
+        c_{rand} = \sum_{\Omega_S} c(\Omega_S) \cdot \cos(\theta) \cdot w
+
+    with the `coefficients` :math:`c`, and the
+    area weights :math:`w` from the `incident_directions`.
+    :math:`\theta` represents the angle between the
+    surface normal and the incident direction.
+
+    Parameters
+    ----------
+    coefficients : pyfar.FrequencyData
+        coefficients for different incident directions. Its cshape
+        needs to be (..., n_incident_directions).
+    colatitude_rad : pyfar.Coordinates
+        Defines the angle between the surface normal and the sound
+        incidence in radiant. Its cshape
+        needs to be (n_incident_directions).
+    area_weights : np.ndarray
+        Area weights for each incident direction. Its cshape
+        needs to be (n_incident_directions).
+
+    Returns
+    -------
+    random_coefficient : pyfar.FrequencyData
+        The random-incidence coefficient.
+
+    References
+    ----------
+    .. [#]  H. Kuttruff, Room acoustics, Sixth edition. Boca Raton:
+            CRC Press/Taylor & Francis Group, 2017.
+    """
+    if not isinstance(coefficients, pf.FrequencyData):
+        raise ValueError("coefficients has to be FrequencyData")
+    if colatitude_rad.shape != area_weights.shape:
+        raise ValueError(
+            "colatitude_rad and area_weights need to have the same shape.")
+
+    theta = colatitude_rad
+    weight = np.cos(theta) * area_weights
+    weight = weight[..., np.newaxis]
+    norm = np.sum(weight)
+    random_coefficient = coefficients*weight/norm
+    random_coefficient.freq = np.sum(random_coefficient.freq, axis=-2)
+    return random_coefficient
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..02c72ac
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,60 @@
+import pytest
+import numpy as np
+import pyfar as pf
+
+from imkar import utils
+
+
+@pytest.mark.parametrize("c_value",  [
+    (0), (0.2), (0.5), (0.8), (1)])
+@pytest.mark.parametrize("frequencies",  [
+    ([100, 200]), ([100])])
+def test_random_constant_coefficient(
+        c_value, frequencies):
+    incident_directions = pf.Coordinates.from_spherical_colatitude(
+        0, np.linspace(0, np.pi/2, 10), 1)
+    incident_directions.weights = np.sin(incident_directions.colatitude)
+    shape = np.append(incident_directions.cshape, len(frequencies))
+    coefficient = pf.FrequencyData(np.zeros(shape)+c_value, frequencies)
+    c_rand = utils.paris_formula(coefficient, incident_directions)
+    np.testing.assert_allclose(c_rand.freq, c_value)
+    assert c_rand.comment == 'random-incidence coefficient'
+
+
+def test_random_non_constant_coefficient():
+    incident_directions = pf.samplings.sph_gaussian(10)
+    incident_directions = incident_directions[incident_directions.z >= 0]
+    incident_cshape = incident_directions.cshape
+    s_value = np.arange(
+        incident_cshape[0]).reshape(incident_cshape) / incident_cshape[0]
+    theta = incident_directions.colatitude
+    actual_weight = np.cos(theta) * incident_directions.weights
+    actual_weight /= np.sum(actual_weight)
+    coefficient = pf.FrequencyData(s_value.reshape((50, 1)), [100])
+    c_rand = utils.paris_formula(coefficient, incident_directions)
+    desired = np.sum(s_value*actual_weight)
+    np.testing.assert_allclose(c_rand.freq, desired)
+
+
+def test_paris_formula_coefficients_type_error():
+    incident_directions = pf.samplings.sph_gaussian(10)
+    with pytest.raises(
+            ValueError, match="coefficients has to be FrequencyData"):
+        utils.paris_formula(np.zeros(10), incident_directions)
+
+
+def test_paris_formula_incident_directions_type_error():
+    coefficients = pf.FrequencyData(np.zeros((10, 1)), [100])
+    with pytest.raises(
+            ValueError,
+            match="incident_directions have to be None or Coordinates"):
+        utils.paris_formula(coefficients, None)
+
+
+def test_paris_formula_shape_mismatch_error():
+    incident_directions = pf.samplings.sph_gaussian(10)
+    coefficients = pf.FrequencyData(np.zeros((5, 1)), [100])
+    with pytest.raises(
+            ValueError,
+            match="the last dimension of coefficients needs be same as"):
+        utils.paris_formula(coefficients, incident_directions)