diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..aea93e3 100644 --- a/pytta/__init__.py +++ b/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 diff --git a/pytta/response.py b/pytta/response.py new file mode 100644 index 0000000..dbf69de --- /dev/null +++ b/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} diff --git a/tests/test_response.py b/tests/test_response.py new file mode 100644 index 0000000..863d7f0 --- /dev/null +++ b/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)