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..8f6d0da --- /dev/null +++ b/pytta/response.py @@ -0,0 +1,247 @@ +""" +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, the slope of a +decaying envelope, the inharmonicity of a tone and phase plane plots. 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 envelope_slope(signalObj, top=-5., bottom=-35.): + """ + Slope, in dB per second, of the amplitude envelope (Hilbert) of a + decaying transient, fitted on its decay from the maximum between `top` + and `bottom` dB below it. + """ + 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 + return np.polyfit(t[fit], level[fit], 1)[0] + + +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, + from envelope_slope between `top` and `bottom` dB below the maximum. + """ + tau = -20 * np.log10(np.e) / envelope_slope(signalObj, top, bottom) + 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} + + +def inharmonicity(signalObj, f0, numPartials=10): + """ + Inharmonicity of a tone with fundamental near f0, in Hz, after the + stiff string model f_n = n f0 sqrt(1 + B n^2) (Fletcher, Normal + vibration frequencies of a stiff piano string, J. Acoust. Soc. Am. 36, + 203, 1964). Each partial is the spectral peak (Hann window, zero padded + 8 times, parabolic interpolation) within f0 / 4 of where the partials + found so far predict it. + + Return: dict with 'f0' and 'B' fitted by least squares on + f_n^2 / n^2 = f0^2 + f0^2 B n^2, 'partials' (Hz) and 'cents', the + deviation of each partial from n f0. + """ + x = np.asarray(signalObj.timeSignal[:, 0], dtype=float) + n_fft = 8 * x.size + f = np.fft.rfftfreq(n_fft, 1 / signalObj.samplingRate) + with np.errstate(divide='ignore'): + L = 20 * np.log10(np.abs(np.fft.rfft(x * np.hanning(x.size), n_fft))) + fit_f0, B, partials = f0, 0., [] + for n in range(1, numPartials + 1): + guess = n * fit_f0 * np.sqrt(1 + B * n ** 2) + band = np.flatnonzero(np.abs(f - guess) <= f0 / 4) + k = band[np.argmax(L[band])] + d = 0.5 * (L[k - 1] - L[k + 1]) / (L[k - 1] - 2 * L[k] + L[k + 1]) + partials.append(f[k] + d * (f[1] - f[0])) + if n >= 2: + m = np.arange(1, n + 1) + c1, c0 = np.polyfit(m ** 2., np.asarray(partials) ** 2 / m ** 2, 1) + fit_f0, B = np.sqrt(c0), c1 / c0 + partials = np.asarray(partials) + cents = 1200 * np.log2(partials / (fit_f0 * np.arange(1, numPartials + 1))) + return {'f0': fit_f0, 'B': B, 'partials': partials, 'cents': cents} + + +def plot_phase_plane(signalObj, quantity='position'): + """ + Phase plane plots of a vibration: position against velocity, velocity + against acceleration and acceleration against position. The signal is + the `quantity` ('position', 'velocity' or 'acceleration'); the others + are derived from it in the frequency domain (multiplying or dividing by + j 2 pi f, DC set to zero), so it is taken as periodic. Returns the + Figure. + """ + import matplotlib.pyplot as plt + order = ['position', 'velocity', 'acceleration'] + x = np.asarray(signalObj.timeSignal[:, 0], dtype=float) + X = np.fft.rfft(x) + jw = 2j * np.pi * np.fft.rfftfreq(x.size, 1 / signalObj.samplingRate) + q = {} + for i, name in enumerate(order): + gain = np.zeros_like(jw) + gain[1:] = jw[1:] ** (i - order.index(quantity)) + q[name] = x if name == quantity else np.fft.irfft(X * gain, x.size) + units = {'position': 'Position', 'velocity': 'Velocity', + 'acceleration': 'Acceleration'} + fig, axes = plt.subplots(1, 3, figsize=(12, 4)) + for ax, (a, b) in zip(axes, [('position', 'velocity'), + ('velocity', 'acceleration'), + ('acceleration', 'position')]): + ax.plot(q[a], q[b]) + ax.set(xlabel=units[a], ylabel=units[b]) + ax.grid(True) + fig.tight_layout() + return fig diff --git a/tests/test_response.py b/tests/test_response.py new file mode 100644 index 0000000..a4050c3 --- /dev/null +++ b/tests/test_response.py @@ -0,0 +1,147 @@ +""" +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) + + +class TestResponseMore(unittest.TestCase): + + def test_envelope_slope_of_an_exponential_decay(self): + # exp(-t / tau) falls by 20 log10(e) / tau dB per second. + label = 'x' + tau, t = 0.05, np.arange(FS // 2) / FS + x = np.exp(-t / tau) * np.sin(2 * np.pi * 1000 * t) + slope = response.envelope_slope(pytta.SignalObj(x, 'time', FS)) + expected = -20 * np.log10(np.e) / tau + print(f'\n envelope slope {slope:.2f} dB/s for {expected:.2f}', end='') + self.assertLess(abs(slope / expected - 1), 1e-3) + + def test_inharmonicity_of_a_stiff_string(self): + # Twelve partials at n f0 sqrt(1 + B n^2) with B = 4e-4. + label = 'x' + f0, B, t = 110., 4e-4, np.arange(2 * FS) / FS + n = np.arange(1, 13) + fn = n * f0 * np.sqrt(1 + B * n ** 2) + x = (np.sin(2 * np.pi * np.outer(t, fn)) / n).sum(axis=1) + result = response.inharmonicity(pytta.SignalObj(x, 'time', FS), 110, 12) + print(f"\n inharmonicity: B {result['B']:.4e}, f0 {result['f0']:.4f}, " + f"partials {np.abs(result['partials'] - fn).max():.1e} Hz", end='') + self.assertLess(abs(result['B'] / B - 1), 0.01) + self.assertLess(abs(result['f0'] - f0), 0.01) + np.testing.assert_allclose(result['cents'], 1200 * np.log2(np.sqrt(1 + B * n ** 2)), atol=0.1) + + def test_phase_plane_of_a_harmonic_oscillator(self): + # x = cos(wt) gives the ellipse x^2 + (v / w)^2 = 1, from the + # position or from the acceleration (float32 samples of SignalObj, + # whose rounding the integration raises at low frequencies). + label = 'x' + w, t = 2 * np.pi * 50, np.arange(FS // 10) / FS + for quantity, x in (('position', np.cos(w * t)), ('acceleration', -w ** 2 * np.cos(w * t))): + ax = response.plot_phase_plane(pytta.SignalObj(x, 'time', FS), quantity).axes[0] + p, v = ax.lines[0].get_xdata(), ax.lines[0].get_ydata() + e = np.abs(p ** 2 + (v / w) ** 2 - 1).max() + print(f'\n phase plane from {quantity}: {e:.1e}', end='') + self.assertLess(e, 1e-4)