[FEATURE] add pytta.descriptors with crest factor, kurtosis and spectral shape

pytta/__init__.py

@@ -108,6 +108,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import descriptors
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +131,7 @@
            'generate',
            'utils',
            'iso3741',
+           'descriptors',
            'rooms',
 
            # Functions

pytta/descriptors.py

@@ -0,0 +1,68 @@
+"""
+PyTTa signal descriptors:
+------------------------
+Descriptors of a whole signal, one value per channel: crest factor and
+kurtosis of the samples, and the centroid, spread, flatness and rolloff of
+the power spectrum.
+
+    >>> pytta.descriptors.crest_factor(signalObj, dB=True)
+
+"""
+import numpy as np
+
+
+def _samples(signalObj):
+    return np.asarray(signalObj.timeSignal, dtype=float)
+
+
+def _power(signalObj):
+    """Frequencies and power spectrum of each channel (one-sided, DC
+    excluded)."""
+    x = _samples(signalObj)
+    f = np.fft.rfftfreq(x.shape[0], 1 / signalObj.samplingRate)
+    return f[1:], np.abs(np.fft.rfft(x, axis=0)[1:]) ** 2
+
+
+def crest_factor(signalObj, dB=False):
+    """Peak magnitude over rms; in dB when dB is True."""
+    x = _samples(signalObj)
+    cf = np.max(np.abs(x), axis=0) / np.sqrt(np.mean(x ** 2, axis=0))
+    return 20 * np.log10(cf) if dB else cf
+
+
+def kurtosis(signalObj):
+    """Pearson kurtosis of the samples: 3 for Gaussian noise, 1.5 for a
+    sine, larger for impulsive signals."""
+    x = _samples(signalObj)
+    x = x - x.mean(axis=0)
+    return np.mean(x ** 4, axis=0) / np.mean(x ** 2, axis=0) ** 2
+
+
+def spectral_centroid(signalObj):
+    """Power-weighted mean frequency, in Hz."""
+    f, P = _power(signalObj)
+    return f @ P / P.sum(axis=0)
+
+
+def spectral_spread(signalObj):
+    """Power-weighted standard deviation of frequency around the
+    centroid, in Hz."""
+    f, P = _power(signalObj)
+    c = f @ P / P.sum(axis=0)
+    return np.sqrt(np.sum((f[:, None] - c) ** 2 * P, axis=0)
+                   / P.sum(axis=0))
+
+
+def spectral_flatness(signalObj):
+    """Geometric over arithmetic mean of the power spectrum: near 0 for a
+    tone, exp(-gamma) = 0.56 for white noise."""
+    f, P = _power(signalObj)
+    with np.errstate(divide='ignore'):
+        return np.exp(np.mean(np.log(P), axis=0)) / np.mean(P, axis=0)
+
+
+def spectral_rolloff(signalObj, fraction=0.85):
+    """Frequency below which `fraction` of the power lies, in Hz."""
+    f, P = _power(signalObj)
+    cumulative = np.cumsum(P, axis=0) / P.sum(axis=0)
+    return f[np.argmax(cumulative >= fraction, axis=0)]

tests/test_descriptors.py

@@ -0,0 +1,64 @@
+"""
+Tests for pytta.descriptors, signal descriptors with exact values for
+known signals. Each test pins one rule.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+
+import pytta
+from pytta import descriptors
+
+FS = 48000
+T = np.arange(FS) / FS
+
+
+class TestDescriptors(unittest.TestCase):
+
+    def test_crest_factor(self):
+        # Peak over rms: sqrt(2) (3.01 dB) for a sine, 1 (0 dB) for a square
+        # wave, per channel.
+        x = np.column_stack((np.sin(2 * np.pi * 1000 * T),
+                             np.sign(np.sin(2 * np.pi * 1000 * T + 0.1))))
+        cf = descriptors.crest_factor(pytta.SignalObj(x, 'time', FS))
+        self.assertTrue(np.allclose(cf, [np.sqrt(2), 1], rtol=1e-4))
+        cf_dB = descriptors.crest_factor(pytta.SignalObj(x, 'time', FS),
+                                         dB=True)
+        self.assertTrue(np.allclose(cf_dB, [3.0103, 0], atol=1e-3))
+
+    def test_kurtosis(self):
+        # Pearson kurtosis: 1.5 for a sine, 3 for Gaussian noise.
+        noise = np.random.default_rng(0).standard_normal(10 * FS)
+        k = descriptors.kurtosis(pytta.SignalObj(noise, 'time', FS))
+        self.assertAlmostEqual(k[0], 3, delta=0.02)
+        k = descriptors.kurtosis(pytta.SignalObj(
+            np.sin(2 * np.pi * 1000 * T), 'time', FS))
+        self.assertAlmostEqual(k[0], 1.5, delta=1e-4)
+
+    def test_spectral_shape_of_a_sine(self):
+        # All the power at f: centroid f, spread 0, flatness about 0.
+        sig = pytta.SignalObj(np.sin(2 * np.pi * 1000 * T), 'time', FS)
+        self.assertAlmostEqual(descriptors.spectral_centroid(sig)[0], 1000,
+                               delta=0.5)
+        self.assertLess(descriptors.spectral_spread(sig)[0], 1)
+        self.assertLess(descriptors.spectral_flatness(sig)[0], 1e-6)
+
+    def test_spectral_shape_of_white_noise(self):
+        # The power of each bin of white Gaussian noise is exponentially
+        # distributed: geometric over arithmetic mean is exp(-gamma) =
+        # 0.5615. The power is flat, so 85 % of it lies below 0.85 fs/2.
+        noise = np.random.default_rng(1).standard_normal(10 * FS)
+        sig = pytta.SignalObj(noise, 'time', FS)
+        self.assertAlmostEqual(descriptors.spectral_flatness(sig)[0],
+                               np.exp(-np.euler_gamma), delta=0.01)
+        self.assertAlmostEqual(descriptors.spectral_rolloff(sig)[0],
+                               0.85 * FS / 2, delta=0.01 * FS / 2)
+        self.assertAlmostEqual(descriptors.spectral_centroid(sig)[0], FS / 4,
+                               delta=0.01 * FS / 2)
+
+
+if __name__ == '__main__':
+    unittest.main()