From 6174d2ee5841a20edaf14b9b5981d57cbe01d1eb Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:25:01 +0100 Subject: [PATCH] [FEATURE] add pytta.descriptors with crest factor, kurtosis and spectral shape Crest factor (ratio or dB), Pearson kurtosis of the samples, and the centroid, spread, flatness and rolloff of the power spectrum, one value per channel. Tests use signals with exact values: crest factor of a sine and a square wave, kurtosis of a sine and of Gaussian noise, the spectral shape of a sine, and the flatness, rolloff and centroid of white noise. Co-Authored-By: Claude Fable 5.1 --- pytta/__init__.py | 2 ++ pytta/descriptors.py | 68 +++++++++++++++++++++++++++++++++++++++ tests/test_descriptors.py | 64 ++++++++++++++++++++++++++++++++++++ 3 files changed, 134 insertions(+) create mode 100644 pytta/descriptors.py create mode 100644 tests/test_descriptors.py diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..6a40e5a 100644 --- a/pytta/__init__.py +++ b/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 diff --git a/pytta/descriptors.py b/pytta/descriptors.py new file mode 100644 index 0000000..04a9d6b --- /dev/null +++ b/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)] diff --git a/tests/test_descriptors.py b/tests/test_descriptors.py new file mode 100644 index 0000000..6793d43 --- /dev/null +++ b/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()