[FEATURE] add pytta.response with quality factor, cutoffs and step metrics

pytta/__init__.py

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

pytta/response.py

@@ -0,0 +1,173 @@
+"""
+PyTTa system response measures:
+------------------------------
+Measures of the response of a linear system, from its impulse response or
+its step response: resonance frequency and quality factor (from the -3 dB
+bandwidth or from the decay), -3 dB cutoff frequencies, passband ripple,
+the rise time, overshoot and settling time of a step response, and the
+onset and hysteresis of a self-sustained oscillation. They use the first
+channel of the SignalObj.
+
+    >>> f0, Q, bandwidth = pytta.response.resonance(ir)
+
+"""
+import numpy as np
+import scipy.signal as ss
+
+
+def _magnitude(signalObj, padding=8):
+    """Frequencies and magnitude spectrum of the first channel, zero padded
+    `padding` times for finer interpolation."""
+    x = np.asarray(signalObj.timeSignal[:, 0], dtype=float)
+    n = padding * x.size
+    return np.fft.rfftfreq(n, 1 / signalObj.samplingRate), \
+        np.abs(np.fft.rfft(x, n))
+
+
+def _crossing(f, y, level, start, step):
+    """Frequency where y falls to `level` walking from `start` by `step`
+    (+1 or -1), linearly interpolated; nan when it never does."""
+    i = start
+    while 0 <= i + step < y.size:
+        if y[i + step] <= level:
+            j = i + step
+            return f[i] + (f[j] - f[i]) * (y[i] - level) / (y[i] - y[j])
+        i += step
+    return np.nan
+
+
+def resonance(signalObj, fmin=0., fmax=None):
+    """
+    Resonance frequency, quality factor and -3 dB bandwidth of the highest
+    peak of the magnitude response between fmin and fmax, in Hz:
+    (f0, Q = f0 / bandwidth, bandwidth). The half-power points are
+    interpolated on the squared magnitude.
+    """
+    f, H = _magnitude(signalObj)
+    band = (f >= fmin) & (f <= (f[-1] if fmax is None else fmax))
+    peak = np.flatnonzero(band)[np.argmax(H[band])]
+    power = H ** 2
+    low = _crossing(f, power, power[peak] / 2, peak, -1)
+    high = _crossing(f, power, power[peak] / 2, peak, +1)
+    return f[peak], f[peak] / (high - low), high - low
+
+
+def q_from_decay(signalObj, f0, top=-5., bottom=-35.):
+    """
+    Quality factor from the decay of a single resonance at f0, in Hz:
+    Q = pi f0 tau, with tau the time constant of the amplitude envelope
+    (Hilbert), fitted on its decay between `top` and `bottom` dB below the
+    maximum.
+    """
+    x = np.asarray(signalObj.timeSignal[:, 0], dtype=float)
+    with np.errstate(divide='ignore'):
+        level = 20 * np.log10(np.abs(ss.hilbert(x)))
+    level -= level.max()
+    start = np.argmax(level)
+    fit = np.arange(x.size) >= start
+    fit &= (level <= top) & (level >= bottom)
+    fit &= np.arange(x.size) < start + np.argmax(level[start:] < bottom)
+    t = np.arange(x.size) / signalObj.samplingRate
+    slope = np.polyfit(t[fit], level[fit], 1)[0]           # dB per second
+    tau = -20 * np.log10(np.e) / slope
+    return np.pi * f0 * tau
+
+
+def cutoff_frequencies(signalObj, reference=None):
+    """
+    Lower and upper -3 dB frequencies, in Hz, around the maximum of the
+    magnitude response (re `reference` in dB, the maximum by default); nan
+    for a side where the response does not fall by 3 dB.
+    """
+    f, H = _magnitude(signalObj)
+    with np.errstate(divide='ignore'):
+        L = 20 * np.log10(H)
+    peak = np.argmax(L)
+    level = (L[peak] if reference is None else reference) - 10 * np.log10(2)
+    return _crossing(f, L, level, peak, -1), _crossing(f, L, level, peak, +1)
+
+
+def passband_ripple(signalObj, fmin, fmax):
+    """Peak-to-valley variation of the magnitude response between fmin
+    and fmax, in dB."""
+    f, H = _magnitude(signalObj)
+    band = (f >= fmin) & (f <= fmax)
+    L = 20 * np.log10(H[band])
+    return L.max() - L.min()
+
+
+def step_metrics(signalObj, final=None, tolerance=0.02):
+    """
+    Rise time (10 % to 90 % of the final value) in s, overshoot as a
+    fraction of the final value, and settling time in s (from the start of
+    the signal until it stays within `tolerance` of the final value) of a
+    step response. The final value is the mean of the last 10 % of the
+    signal unless given.
+    """
+    y = np.asarray(signalObj.timeSignal[:, 0], dtype=float)
+    fs = signalObj.samplingRate
+    if final is None:
+        final = np.mean(y[-max(1, y.size // 10):])
+
+    def first_time(fraction):
+        reached = np.flatnonzero(y >= fraction * final)
+        if reached.size == 0 or reached[0] == 0:
+            return np.nan
+        i = reached[0]
+        return (i - 1 + (fraction * final - y[i - 1]) / (y[i] - y[i - 1])) / fs
+
+    rise = first_time(0.9) - first_time(0.1)
+    overshoot = max(0.0, (y.max() - final) / final)
+    error = np.abs(y - final) - tolerance * abs(final)
+    outside = np.flatnonzero(error > 0)
+    if outside.size == 0:
+        return rise, overshoot, 0.0
+    i = outside[-1]
+    if i == y.size - 1:                       # never settles
+        return rise, overshoot, np.nan
+    settling = (i + error[i] / (error[i] - error[i + 1])) / fs
+    return rise, overshoot, settling
+
+
+def bifurcation(params, recordings, freq=None, rise=20.):
+    """
+    Onset, extinction and hysteresis of a self-sustained oscillation along
+    a sweep of a control parameter (e.g. gain or blowing pressure), first
+    up and then down.
+
+    Parameters (default), (type):
+    -----------------------------
+        * params (), (array-like): value of the control parameter at each
+          recording, increasing up to its maximum and then decreasing;
+        * recordings (), (list of SignalObj): one recording per value;
+        * freq (None), (float): oscillation frequency in Hz; its amplitude
+          is the magnitude of the spectral peak within 5 % of freq. With
+          None, the rms of each recording is used;
+        * rise (20), (float): how many dB above the quietest recording an
+          amplitude has to be to count as oscillating.
+
+    Return: dict with 'amplitude' (one value per recording), 'onset' (the
+    first value going up that oscillates), 'extinction' (the last value
+    going down that still oscillates) and 'hysteresis' (onset minus
+    extinction). A missing onset or extinction is nan.
+    """
+    params = np.asarray(params, dtype=float)
+    amplitude = []
+    for rec in recordings:
+        x = np.asarray(rec.timeSignal[:, 0], dtype=float)
+        if freq is None:
+            amplitude.append(np.sqrt(np.mean(x ** 2)))
+            continue
+        f = np.fft.rfftfreq(x.size, 1 / rec.samplingRate)
+        X = 2 * np.abs(np.fft.rfft(x * np.hanning(x.size))) \
+            / np.hanning(x.size).sum()
+        amplitude.append(X[np.abs(f - freq) <= 0.05 * freq].max())
+    amplitude = np.asarray(amplitude)
+    oscillating = 20 * np.log10(amplitude / amplitude.min()) >= rise
+    top = int(np.argmax(params))
+    up = np.flatnonzero(oscillating[:top + 1])
+    down = np.flatnonzero(oscillating[top:]) + top
+    onset = params[up[0]] if up.size else np.nan
+    extinction = params[down[-1]] if down.size else np.nan
+    return {'amplitude': amplitude, 'onset': onset,
+            'extinction': extinction, 'hysteresis': onset - extinction}

tests/test_response.py

@@ -0,0 +1,107 @@
+"""
+Tests for pytta.response, measures of the response of a linear system:
+resonance frequency and quality factor (from the -3 dB bandwidth and from
+the decay), -3 dB cutoff frequencies, passband ripple, and the rise time,
+overshoot and settling time of a step response. Each test pins one rule
+on a system whose values are known exactly.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+import scipy.signal as ss
+
+import pytta
+from pytta import response
+
+FS = 48000
+
+
+def resonator(f0, Q, seconds=2.0):
+    """Impulse response of a second-order band-pass resonator."""
+    b, a = ss.iirpeak(f0, Q, fs=FS)
+    return ss.lfilter(b, a, np.eye(1, int(seconds * FS))[0])
+
+
+class TestResponse(unittest.TestCase):
+
+    def test_quality_factor_from_bandwidth_and_from_decay(self):
+        # Both definitions give the Q of a second-order resonance.
+        for f0, Q in ((1000.0, 10.0), (250.0, 30.0)):
+            ir = pytta.SignalObj(resonator(f0, Q), 'time', FS)
+            fr, Qb, bw = response.resonance(ir)
+            Qd = response.q_from_decay(ir, f0)
+            print(f'\n  f0 {f0}: found {fr:.2f} Hz, Q {Qb:.3f} (bandwidth), '
+                  f'{Qd:.3f} (decay), expected {Q}')
+            self.assertAlmostEqual(fr / f0, 1, delta=1e-3)
+            self.assertAlmostEqual(Qb / Q, 1, delta=0.01)
+            self.assertAlmostEqual(Qd / Q, 1, delta=0.02)
+
+    def test_cutoffs_and_ripple_of_designed_filters(self):
+        # A Butterworth band-pass is -3.01 dB at its design edges; a
+        # Chebyshev type I low-pass ripples by rp in its passband.
+        sos = ss.butter(4, [500, 2000], 'bandpass', fs=FS, output='sos')
+        ir = pytta.SignalObj(ss.sosfilt(sos, np.eye(1, FS)[0]), 'time', FS)
+        low, high = response.cutoff_frequencies(ir)
+        self.assertAlmostEqual(low / 500, 1, delta=2e-3)
+        self.assertAlmostEqual(high / 2000, 1, delta=2e-3)
+        sos = ss.cheby1(6, 1.0, 4000, fs=FS, output='sos')
+        ir = pytta.SignalObj(ss.sosfilt(sos, np.eye(1, FS)[0]), 'time', FS)
+        self.assertAlmostEqual(response.passband_ripple(ir, 20, 3500), 1.0,
+                               delta=0.02)
+
+    def test_step_response_of_a_first_order_system(self):
+        # y = 1 - exp(-t / tau): rise time (10 to 90 %) tau ln 9, no
+        # overshoot, settling time within 2 % tau ln 50.
+        tau = 0.01
+        t = np.arange(FS // 2) / FS
+        step = pytta.SignalObj(1 - np.exp(-t / tau), 'time', FS)
+        rise, overshoot, settling = response.step_metrics(step)
+        self.assertAlmostEqual(rise, tau * np.log(9), delta=2 / FS)
+        self.assertAlmostEqual(overshoot, 0, delta=1e-9)
+        self.assertAlmostEqual(settling, tau * np.log(50), delta=2 / FS)
+
+    def test_overshoot_of_a_second_order_system(self):
+        # Underdamped step response with damping ratio zeta: overshoot
+        # exp(-pi zeta / sqrt(1 - zeta^2)).
+        for zeta in (0.2, 0.5):
+            wn = 2 * np.pi * 50
+            b, a = ss.bilinear([wn ** 2], [1, 2 * zeta * wn, wn ** 2], FS)
+            y = ss.lfilter(b, a, np.ones(FS))
+            _, overshoot, _ = response.step_metrics(
+                pytta.SignalObj(y, 'time', FS))
+            expected = np.exp(-np.pi * zeta / np.sqrt(1 - zeta ** 2))
+            print(f'\n  zeta {zeta}: overshoot {overshoot:.4f} (expected '
+                  f'{expected:.4f})')
+            self.assertAlmostEqual(overshoot, expected, delta=2e-3)
+
+
+if __name__ == '__main__':
+    unittest.main()
+
+
+class TestBifurcation(unittest.TestCase):
+
+    def test_onset_extinction_and_hysteresis(self):
+        # Subcritical behaviour: going up, the oscillation starts at
+        # mu = 1.0; coming down, it persists until mu = 0.6. Below that only
+        # noise 40 dB down. Hysteresis 0.4.
+        up = np.round(np.arange(0.2, 1.41, 0.1), 2)
+        down = up[::-1][1:]
+        params = np.r_[up, down]
+        rng = np.random.default_rng(0)
+        t = np.arange(FS // 4) / FS
+        recordings = []
+        for i, mu in enumerate(params):
+            on = mu >= 1.0 if i < up.size else mu >= 0.6
+            amp = (0.5 + 0.2 * (mu - 0.6)) if on else 0.0
+            x = amp * np.sin(2 * np.pi * 220 * t) \
+                + 0.005 * rng.standard_normal(t.size)
+            recordings.append(pytta.SignalObj(x, 'time', FS))
+        result = response.bifurcation(params, recordings, freq=220)
+        self.assertAlmostEqual(result['onset'], 1.0)
+        self.assertAlmostEqual(result['extinction'], 0.6)
+        self.assertAlmostEqual(result['hysteresis'], 0.4)
+        self.assertEqual(len(result['amplitude']), params.size)