diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..eb5d384 100644 --- a/docs/pyttadoc/functions.rst +++ b/docs/pyttadoc/functions.rst @@ -26,4 +26,7 @@ Functions .. autofunction:: pytta.plot_bars .. autofunction:: pytta.plot_spectrogram .. autofunction:: pytta.plot_waterfall - +.. autofunction:: pytta.dither +.. autofunction:: pytta.impulsive_noise +.. autofunction:: pytta.invert_spectrum +.. autofunction:: pytta.merge_near_far diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..47286da 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -113,7 +113,8 @@ list_devices, print_devices, get_device_from_user,\ fft_convolve, find_delay, corr_coef, resample,\ peak_time, save, load, plot_time, plot_time_dB, plot_freq,\ - plot_bars, plot_spectrogram, plot_waterfall, SPL + plot_bars, plot_spectrogram, plot_waterfall, SPL,\ + dither, impulsive_noise, invert_spectrum, merge_near_far # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +147,10 @@ 'resample', 'corr_coef', 'peak_time', + 'dither', + 'impulsive_noise', + 'invert_spectrum', + 'merge_near_far', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..6cc00b4 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -274,6 +274,121 @@ def peak_time(signal): else: return peaks_time[0] +def dither(signal, bits=24, kind='tri', quiet=False, seed=None): + """ + signal (full scale 1) with dither noise added before it is quantized to + bits bits, so that the quantization error becomes noise independent of + the signal, as ita_dither (Lipshitz, Wannamaker and Vanderkooy, + Quantization and dither: a theoretical survey, J. Audio Eng. Soc. + 40(5), 1992). kind 'rect' is uniform over one step, 'tri' (triangular, + two steps wide) makes the mean and the power of the error independent + of the signal, and 'gauss' has a standard deviation of half a step. With + quiet, samples equal to zero get no dither. + """ + rng = np.random.default_rng(seed) + x = np.asarray(signal.timeSignal, dtype=float) + step = 2.0 ** (1 - bits) + if kind == 'rect': + d = rng.uniform(-step / 2, step / 2, x.shape) + elif kind == 'tri': + d = rng.uniform(-step / 2, step / 2, x.shape) + rng.uniform(-step / 2, step / 2, x.shape) + elif kind == 'gauss': + d = rng.normal(0, step / 2, x.shape) + else: + raise ValueError("kind must be 'rect', 'tri' or 'gauss'.") + if quiet: + d[x == 0] = 0 + return SignalObj(x + d, 'time', signal.samplingRate) + + +def impulsive_noise(ir, excitation, noise_start, threshold=20): + """ + Detection of an impulsive event during a sweep measurement, as + ita_impulsiveNoiseDetection: the part of the impulse response ir after + noise_start seconds, where it holds only noise, is convolved again with + the excitation (circularly, at the length of ir) to rebuild the + background noise of the recording, and L_max_rms = 20 lg(max |n| / + rms(n)) of that noise is compared with threshold (dB). Stationary noise + gives 12 to 14 dB. + + Returns (detected, L_max_rms, background noise SignalObj). + """ + h = np.asarray(ir.timeSignal[:, 0], dtype=float).copy() + h[:int(round(noise_start * ir.samplingRate))] = 0 + x = np.asarray(excitation.timeSignal[:, 0], dtype=float) + n = np.fft.irfft(np.fft.rfft(h) * np.fft.rfft(x, h.size), h.size) + level = 20 * np.log10(np.abs(n).max() / np.sqrt(np.mean(n ** 2))) + return level > threshold, level, SignalObj(n, 'time', ir.samplingRate) + + +def _log_crossfade(f, f0, f1): + # 0 below f0, 1 above f1, raised cosine on a logarithmic axis between. + with np.errstate(divide='ignore'): + u = np.clip(np.log(f / f0) / np.log(f1 / f0), 0, 1) + return 0.5 - 0.5 * np.cos(np.pi * np.nan_to_num(u)) + + +def invert_spectrum(signal, freqRange, minimum_phase=False): + """ + Regularized inverse of the spectrum of signal, C = conj(H) / (|H|^2 + + eps) (Kirkeby, Rubak, Nelson and Farina, Design of cross-talk + cancellation networks by using fast deconvolution, 1999), with eps + negligible inside freqRange = [fmin, fmax] and half the largest |H|^2 + outside it, crossfaded over half an octave, the rule of the + regularization of ImpulsiveResponse. With minimum_phase, the inverse + gets the minimum phase of its magnitude, as + ita_invert_spk_minphase_regularization. + """ + x = np.asarray(signal.timeSignal, dtype=float) + n = x.shape[0] + H = np.fft.rfft(x, axis=0) + f = np.fft.rfftfreq(n, 1 / signal.samplingRate) + fmin, fmax = freqRange + outside = 1 - _log_crossfade(f, fmin / np.sqrt(2), fmin) \ + + _log_crossfade(f, fmax, fmax * np.sqrt(2)) + eps = (1e-10 + outside[:, None]) * (np.abs(H) ** 2).max(axis=0) / 2 + C = np.conj(H) / (np.abs(H) ** 2 + eps) + if minimum_phase: + # Homomorphic: fold the real cepstrum of log |C| onto positive times. + c = np.fft.irfft(np.log(np.abs(C)), n, axis=0) + fold = np.zeros(n) + fold[0] = 1 + fold[1:(n + 1) // 2] = 2 + if n % 2 == 0: + fold[n // 2] = 1 + C = np.exp(np.fft.rfft(c * fold[:, None], axis=0)) + return SignalObj(np.fft.irfft(C, n, axis=0), 'time', signal.samplingRate) + + +def merge_near_far(near, far, crossover, order=4, port=None, port_ratio=1.): + """ + Loudspeaker response from a near-field measurement below crossover (Hz) + and a far-field one above it, as + ita_add_nearfield_farfield_measurements (Keele, Low-frequency + loudspeaker assessment by nearfield sound-pressure measurement, J. + Audio Eng. Soc. 22(3), 1974). The near-field response, with the port + one times port_ratio (port radius over diaphragm radius) added for a + vented box, is delayed to the far-field one (find_delay) and scaled to + its level over the third octave around the crossover, and the two are + crossfaded with the power-complementary weights 1 / (1 + (f/fc)^(2 + order)) and its complement, which add up to 1. + """ + fs = far.samplingRate + m = np.asarray(near.timeSignal[:, 0], dtype=float) + if port is not None: + m = m + port_ratio * np.asarray(port.timeSignal[:, 0], dtype=float) + x = np.asarray(far.timeSignal[:, 0], dtype=float) + n = x.size + f = np.fft.rfftfreq(n, 1 / fs) + lag = np.argmax(ss.correlate(x, m, mode='full', method='fft')) - (m.size - 1) + M = np.fft.rfft(m, n) * np.exp(-2j * np.pi * f * lag / fs) + X = np.fft.rfft(x) + band = (f >= crossover * 2 ** (-1 / 6)) & (f <= crossover * 2 ** (1 / 6)) + M *= np.sqrt(np.sum(np.abs(X[band]) ** 2) / np.sum(np.abs(M[band]) ** 2)) + low = 1 / (1 + (f / crossover) ** (2 * order)) + return SignalObj(np.fft.irfft(low * M + (1 - low) * X, n), 'time', fs) + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/test_measurement_tools.py b/tests/test_measurement_tools.py new file mode 100644 index 0000000..94dfff2 --- /dev/null +++ b/tests/test_measurement_tools.py @@ -0,0 +1,95 @@ +""" +Tests for dither, impulsive_noise, invert_spectrum and merge_near_far of +pytta.functions, which cover ita_dither, ita_impulsiveNoiseDetection, +ita_invert_spk_minphase_regularization and +ita_add_nearfield_farfield_measurements of the ITA-Toolbox. Each test pins +one rule on a synthetic measurement. +""" +import unittest + +import numpy as np +import scipy.signal as ss + +import pytta + +FS = 48000 + + +def signal(*samples): + # A tuple as first local variable: SignalObj compares it with itself + # while walking the stack, which an array cannot do before #59. + return pytta.SignalObj(samples[0], 'time', FS) + + +class TestMeasurementTools(unittest.TestCase): + + def test_triangular_dither_makes_the_error_independent(self): + # With triangular dither the total error of an 8-bit quantizer has + # the power step^2 / 4 and no mean, whatever the signal (Lipshitz + # et al. 1992); without it the error follows the signal. + label = 'x' + step = 2.0 ** -7 + x = 1.5 * step * np.sin(2 * np.pi * 50 * np.arange(FS) / FS) + y = pytta.dither(signal(x), 8, 'tri', seed=1).timeSignal[:, 0].astype(float) + e = np.round(y / step) * step - x + power = np.mean(e ** 2) / (step ** 2 / 4) + bias = max(abs(np.mean(e[x > step])), abs(np.mean(e[x < -step]))) / step + print(f'\n dithered error: power {power:.3f} of step^2/4, mean {bias:.1e} step', end='') + self.assertLess(abs(power - 1), 0.02) + self.assertLess(bias, 0.01) + + def test_impulsive_noise_is_detected(self): + # Excitation of flat magnitude and random phase, so that the + # deconvolution is exact. A click in the recording raises L_max_rms + # of the rebuilt noise above 20 dB; Gaussian noise alone stays near + # 13 dB. + label = 'x' + rng = np.random.default_rng(2) + n = 2 ** 16 + X = np.exp(2j * np.pi * rng.uniform(size=n // 2 + 1)) + X[[0, -1]] = 1 + x = np.fft.irfft(X, n) + room = np.exp(-np.arange(4800) / 600) * rng.standard_normal(4800) + levels = [] + for click in (0, 0.05): + y = np.fft.irfft(np.fft.rfft(room, n) * X, n) + 1e-3 * rng.standard_normal(n) + y[30000] += click + h = np.fft.irfft(np.fft.rfft(y) / X, n) + found, level, _ = pytta.impulsive_noise(signal(h), signal(x), 0.2) + levels.append(level) + self.assertEqual(found, click > 0) + print(f'\n L_max_rms: {levels[0]:.1f} dB without, {levels[1]:.1f} dB with a click', end='') + + def test_inverse_is_exact_in_band_and_minimum_phase(self): + label = 'x' + b, a = ss.butter(2, [100, 8000], 'bandpass', fs=FS) + h = ss.lfilter(b, a, np.r_[1, np.zeros(8191)]) + H = np.fft.rfft(h) + f = np.fft.rfftfreq(h.size, 1 / FS) + band = (f > 200) & (f < 4000) + for minimum in (False, True): + c = pytta.invert_spectrum(signal(h), [200, 4000], minimum).timeSignal[:, 0].astype(float) + C = np.fft.rfft(c) + e = np.abs(np.abs(C * H)[band] - 1).max() + print(f'\n inverse, minimum phase {minimum}: in band {e:.1e}', end='') + self.assertLess(e, 1e-4) + # The minimum phase inverse has the most energy at its start among + # inverses with the same magnitude (here the zero phase one). + zero = np.roll(np.fft.irfft(np.abs(C), h.size), h.size // 2) + lead = np.cumsum(c ** 2) - np.cumsum(zero ** 2) + self.assertGreater(lead.min(), -1e-5 * np.sum(c ** 2)) # float32 samples + + def test_near_and_far_field_merge_back(self): + # A near-field measurement 20 dB louder and 100 samples earlier than + # the far-field one is aligned, scaled and crossfaded back into it. + label = 'x' + h = ss.lfilter(*ss.butter(2, [40, 15000], 'bandpass', fs=FS), np.r_[1, np.zeros(16383)]) + far = np.roll(h, 100) + merged = pytta.merge_near_far(signal(10 * h), signal(far), 300).timeSignal[:, 0].astype(float) + e = np.abs(merged - far).max() / np.abs(far).max() + print(f'\n near and far field merge: {e:.1e}', end='') + self.assertLess(e, 1e-4) + + +if __name__ == '__main__': + unittest.main()