diff --git a/pytta/generate.py b/pytta/generate.py
index 1979f54..5b440b5 100644
--- a/pytta/generate.py
+++ b/pytta/generate.py
@@ -695,6 +695,39 @@ def measurement(kind='playrec',
         return frfObj
 
 
+def stochastic_ir(reverbTime, volume, distance=1, c=343, samplingRate=None,
+                  timeLength=None, seed=None):
+    """
+    Synthetic room impulse response of a diffuse room, after the statistical
+    model of Kuttruff, Room Acoustics, 6th ed. (2016), sec. 4.2, eq. 4.6 and
+    4.7, as ita_stochastic_room_impulse_response without an HRTF: the direct sound
+    1 / distance at distance / c, and reflections arriving as a Poisson
+    process of density 4 pi c^3 t^2 / volume, each of random sign and of
+    amplitude 10^(-3 t / reverbTime) / (c t), so that the energy decays by
+    60 dB in reverbTime. timeLength defaults to 1.2 reverbTime.
+
+    >>> ir = pytta.generate.stochastic_ir(0.8, 200, distance=3, seed=1)
+    """
+    if samplingRate is None:
+        samplingRate = default.samplingRate
+    if timeLength is None:
+        timeLength = 1.2 * reverbTime
+    rng = np.random.default_rng(seed)
+    n = int(round(timeLength * samplingRate))
+    t0 = distance / c
+    # Unit rate Poisson points mapped through the inverse of the integrated
+    # density, Lambda(t) = 4 pi c^3 (t^3 - t0^3) / (3 V).
+    k = 4 * np.pi * c**3 / (3 * volume)
+    total = k * (timeLength**3 - t0**3)
+    u = np.sort(rng.uniform(0, total, rng.poisson(total)))
+    t = np.cbrt(t0**3 + u / k)
+    x = np.zeros(n)
+    np.add.at(x, np.minimum((t * samplingRate).round().astype(int), n - 1),
+              rng.choice([-1, 1], t.size) * 10**(-3 * t / reverbTime) / (c * t))
+    x[int(round(t0 * samplingRate))] += 1 / distance
+    return SignalObj(x, 'time', samplingRate, signalType='energy')
+
+
 def stream(IO='IO',
            device=None,
            integration=None,
diff --git a/tests/test_stochastic_ir.py b/tests/test_stochastic_ir.py
new file mode 100644
index 0000000..b7f68a5
--- /dev/null
+++ b/tests/test_stochastic_ir.py
@@ -0,0 +1,47 @@
+"""
+Tests for pytta.generate.stochastic_ir, which follows the statistical model
+of Kuttruff, Room Acoustics, 6th ed. (2016), sec. 4.2, as
+ita_stochastic_room_impulse_response of the ITA-Toolbox does. Each test pins
+one property of the model.
+"""
+import unittest
+
+import numpy as np
+
+import pytta
+
+FS, V, T, R, C = 48000, 200, 0.8, 3, 343
+
+
+class TestStochasticIR(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        cls.x = pytta.generate.stochastic_ir(T, V, R, C, FS, seed=4).timeSignal[:, 0].astype(float)
+
+    def test_direct_sound_comes_first(self):
+        i = round(R / C * FS)
+        self.assertFalse(self.x[:i].any())
+        self.assertAlmostEqual(self.x[i], 1 / R, places=6)
+
+    def test_reflection_density(self):
+        # Number of reflections between 20 and 40 ms against the integral of
+        # 4 pi c^3 t^2 / V (eq. 4.6), where they rarely share a sample.
+        t1, t2 = 0.02, 0.04
+        count = np.count_nonzero(self.x[round(t1 * FS):round(t2 * FS)])
+        expected = 4 * np.pi * C**3 * (t2**3 - t1**3) / (3 * V)
+        print(f'\n  reflections: {count}, expected {expected:.0f}', end='')
+        self.assertLess(abs(count - expected), 4 * np.sqrt(expected))
+
+    def test_decay_gives_the_reverberation_time(self):
+        # T30 from the backward integrated energy (-5 to -35 dB).
+        e = np.cumsum(self.x[::-1] ** 2)[::-1]
+        L = 10 * np.log10(e / e[0])
+        i5, i35 = np.argmax(L < -5), np.argmax(L < -35)
+        slope = np.polyfit(np.arange(i5, i35) / FS, L[i5:i35], 1)[0]
+        print(f'\n  T30 {-60 / slope:.3f} s for {T} s', end='')
+        self.assertLess(abs(-60 / slope / T - 1), 0.02)
+
+
+if __name__ == '__main__':
+    unittest.main()