diff --git a/src/__init__.py b/src/__init__.py
index d5e0b9b9ec70e0b4714ad2ad85b930e48c3b677b..787e158ab7c6cf893696e6cfd0adace350a27e35 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -77,7 +77,8 @@ from .sugar import *
 
 from .plot import Plot
 
-from .library.special_distributions import InverseGammaOperator, UniformOperator
+from .library.special_distributions import (InverseGammaOperator, UniformOperator,
+                                            LaplaceOperator)
 from .library.los_response import LOSResponse
 from .library.dynamic_operator import (dynamic_operator,
                                        dynamic_lightcone_operator)
diff --git a/src/library/special_distributions.py b/src/library/special_distributions.py
index a46f17c044f6abde8fd692dabca49cc08d894aaa..4913ce0f36bafd1a1217ea657c5aef3fc3823b63 100644
--- a/src/library/special_distributions.py
+++ b/src/library/special_distributions.py
@@ -17,7 +17,7 @@
 
 import numpy as np
 from scipy.interpolate import CubicSpline
-from scipy.stats import invgamma, norm
+from scipy.stats import invgamma, norm, laplace
 
 from .. import random
 from ..domain_tuple import DomainTuple
@@ -158,3 +158,29 @@ class UniformOperator(Operator):
     def inverse(self, field):
         res = norm._ppf(field.val/self._scale - self._loc)
         return Field(field.domain, res)
+
+
+def LaplaceOperator(domain, alpha, delta=1e-2):
+    """Transforms a Gaussian with unit covariance and zero mean into a
+    Laplacian distribution.
+
+    The pdf of the laplace distribution is defined as follows:
+
+    .. math::
+        \\frac{\\alpha}{2} \\exp \\left(-\\alpha |x| \\right)
+
+    This transformation is implemented as a spline interpolation which maps a
+    Gaussian onto a Laplace distribution.
+
+    Parameters
+    ----------
+    domain : Domain, tuple of Domain or DomainTuple
+        The domain on which the field shall be defined. This is at the same
+        time the domain and the target of the operator.
+    alpha : float
+        The alpha-parameter of the Laplace distribution.
+    delta : float
+        Distance between sampling points for spline interpolation.
+    """
+    return _InterpolationOperator(domain, lambda x: laplace.ppf(norm._cdf(x), float(alpha)),
+                                  -8.2, 8.2, delta)