diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..d61eb01 100644 --- a/docs/pyttadoc/functions.rst +++ b/docs/pyttadoc/functions.rst @@ -26,4 +26,12 @@ Functions .. autofunction:: pytta.plot_bars .. autofunction:: pytta.plot_spectrogram .. autofunction:: pytta.plot_waterfall - +.. autofunction:: pytta.envelope +.. autofunction:: pytta.start_ir +.. autofunction:: pytta.time_shift +.. autofunction:: pytta.append +.. autofunction:: pytta.normalize +.. autofunction:: pytta.amplify_to +.. autofunction:: pytta.zero_phase +.. autofunction:: pytta.randomize_phase +.. autofunction:: pytta.quantize diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..d2927b8 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -113,7 +113,9 @@ 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,\ + envelope, start_ir, time_shift, append, normalize, amplify_to,\ + zero_phase, randomize_phase, quantize # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +148,15 @@ 'resample', 'corr_coef', 'peak_time', + 'envelope', + 'start_ir', + 'time_shift', + 'append', + 'normalize', + 'amplify_to', + 'zero_phase', + 'randomize_phase', + 'quantize', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..f318da7 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -23,6 +23,15 @@ >>> pytta.plot_bars(signalObj1, signalObj2, ..., signalObjN) >>> pytta.save(fileName, obj1, ..., objN) >>> pytta.load(fileName) + >>> pytta.envelope(signalObj) + >>> pytta.start_ir(signalObj) + >>> pytta.time_shift(signalObj, shift) + >>> pytta.append(signalObj1, signalObj2) + >>> pytta.normalize(signalObj) + >>> pytta.amplify_to(signalObj, level) + >>> pytta.zero_phase(signalObj) + >>> pytta.randomize_phase(signalObj) + >>> pytta.quantize(signalObj, bits) For further information, check the function specific documentation. @@ -274,6 +283,190 @@ def peak_time(signal): else: return peaks_time[0] +def _new_signal(signal, data, domain='time'): + """SignalObj like signal (sampling rate, channels) holding data.""" + newSignal = SignalObj(data, domain, samplingRate=signal.samplingRate, + signalType=signal.signalType) + newSignal.channels = cp.deepcopy(signal.channels) + return newSignal + + +def envelope(signal): + """ + Envelope of each channel, the magnitude of the analytic signal + (Hilbert transform), computed on twice the signal length to keep the + circular wrap out of it, as ita_envelope. + + >>> env = pytta.envelope(signalObj) + """ + x = np.asarray(signal.timeSignal, dtype=float) + n = x.shape[0] + return _new_signal(signal, np.abs(ss.hilbert(x, N=2 * n, axis=0)[:n])) + + +def start_ir(signal, threshold=20.): + """ + Start of each channel of an impulse response, ISO 3382-1:2009, A.3: the + last sample before the peak that lies more than threshold dB below the + peak of the squared response, as ita_start_IR. + + If a sample 6 dB above that level comes earlier, the start moves to the + last sample below the threshold before it (an oscillating onset). A + channel whose peak is less than 20 dB above its noise (the mean of the + last 10 %) or lies in that last 10 % gets 0. + + Returns + ------- + numpy.ndarray + 0-based start sample of each channel. + + >>> pytta.start_ir(ir) + """ + x2 = np.asarray(signal.timeSignal, dtype=float) ** 2 + n = x2.shape[0] + noise = np.mean(x2[int(round(0.9 * n)) - 1:], axis=0) + peak, at = np.max(x2, axis=0), np.argmax(x2, axis=0) + thr = -abs(threshold) + start = np.zeros(x2.shape[1], dtype=int) + for c in range(x2.shape[1]): + if peak[c] < 100 * noise[c] or at[c] + 1 > round(0.9 * n) or at[c] == 0: + continue + with np.errstate(divide='ignore'): + level = 10 * np.log10(x2[:at[c] + 1, c]) - 10 * np.log10(peak[c]) + below = np.flatnonzero(level < thr) + start[c] = below[-1] if below.size else 0 + above = np.flatnonzero(level[:start[c] + 1] > thr + 6) + if above.size: + before = np.flatnonzero(level[:above[0] + 1] < thr) + start[c] = before[-1] if before.size else 0 + return start + + +def time_shift(signal, shift=None, unit='time'): + """ + Cyclic shift of each channel, as ita_time_shift. + + Parameters + ---------- + signal : SignalObj + shift : float, array_like or None + Shift in seconds (unit='time') or samples (unit='samples'), one for + all channels or one per channel; a fraction of a sample is applied + as a linear phase. None moves the earliest start of the impulse + response (start_ir) to the first sample. + unit : str + 'time' (default) or 'samples'. + + >>> shifted = pytta.time_shift(ir, 1e-3) + """ + x = np.asarray(signal.timeSignal, dtype=float) + fs = signal.samplingRate + if shift is None: + samples = np.full(x.shape[1], -np.min(start_ir(signal)), dtype=float) + elif unit == 'samples': + samples = np.broadcast_to(np.asarray(shift, dtype=float), (x.shape[1],)) + elif unit == 'time': + samples = np.broadcast_to(np.round(np.asarray(shift, dtype=float) * fs), (x.shape[1],)) + else: + raise ValueError("unit must be 'time' or 'samples'.") + out = np.empty_like(x) + n = x.shape[0] + for c in range(x.shape[1]): + k = samples[c] + if k == round(k): + out[:, c] = np.roll(x[:, c], int(k)) + else: + f = np.fft.rfftfreq(n) + out[:, c] = np.fft.irfft(np.fft.rfft(x[:, c]) * np.exp(-2j * np.pi * f * k), n) + return _new_signal(signal, out) + + +def append(signal1, signal2): + """ + signal2 after signal1 in time, as ita_append; both need the same + sampling rate and number of channels. + + >>> longer = pytta.append(signalObj1, signalObj2) + """ + if signal1.samplingRate != signal2.samplingRate: + raise ValueError('the sampling rates differ.') + if signal1.numChannels != signal2.numChannels: + raise ValueError('the numbers of channels differ.') + return _new_signal(signal1, np.vstack([np.asarray(signal1.timeSignal, dtype=float), + np.asarray(signal2.timeSignal, dtype=float)])) + + +def normalize(signal, domain='time', per_channel=False): + """ + Signal scaled so that its largest magnitude, in time or in the + spectrum, is 1, as ita_normalize_dat and ita_normalize_spk. + + Returns + ------- + (SignalObj, numpy.ndarray) + The scaled signal and the gain removed (one value, or one per + channel with per_channel=True). + + >>> normed, peak = pytta.normalize(signalObj) + """ + x = np.asarray(signal.timeSignal, dtype=float) + if domain == 'time': + mag = np.abs(x) + elif domain == 'freq': + mag = np.abs(np.fft.rfft(x, axis=0)) + else: + raise ValueError("domain must be 'time' or 'freq'.") + gain = np.max(mag, axis=0) if per_channel else np.max(mag) + return _new_signal(signal, x / gain), np.atleast_1d(gain) + + +def amplify_to(signal, level): + """ + Signal scaled so that the mean RMS of its channels is level dB, re + 20 uPa for channels in Pa and re 1 otherwise, as ita_amplify_to. + + >>> at_70dB = pytta.amplify_to(signalObj, 70) + """ + x = np.asarray(signal.timeSignal, dtype=float) + ref = 2e-5 if signal.channels[signal.channels.mapping[0]].unit == 'Pa' else 1. + rms = np.mean(np.sqrt(np.mean(x ** 2, axis=0))) + return _new_signal(signal, x * 10 ** ((level - 20 * np.log10(rms / ref)) / 20)) + + +def zero_phase(signal): + """Signal with the magnitude spectrum of signal and zero phase (ita_zerophase).""" + x = np.asarray(signal.timeSignal, dtype=float) + return _new_signal(signal, np.fft.irfft(np.abs(np.fft.rfft(x, axis=0)), x.shape[0], axis=0)) + + +def randomize_phase(signal, seed=None): + """ + Signal with the magnitude spectrum of signal and a uniformly random + phase (ita_randomize_phase); DC and Nyquist keep their sign so that the + result is real. + """ + x = np.asarray(signal.timeSignal, dtype=float) + X = np.fft.rfft(x, axis=0) + phase = np.random.default_rng(seed).uniform(0, 2 * np.pi, X.shape) + phase[0] = np.angle(X[0]) + if x.shape[0] % 2 == 0: + phase[-1] = np.angle(X[-1]) + return _new_signal(signal, np.fft.irfft(np.abs(X) * np.exp(1j * phase), x.shape[0], axis=0)) + + +def quantize(signal, bits=24): + """ + Signal rounded to 2^(bits - 1) + 1 intervals across its peak, halves + away from zero, as ita_quantize. + """ + x = np.asarray(signal.timeSignal, dtype=float) + peak = np.max(np.abs(x)) + steps = 2 ** (bits - 1) + 1 + v = x / peak * steps / 2 + v = np.sign(v) * np.floor(np.abs(v) + 0.5) # halves away from zero + return _new_signal(signal, v / steps * 2 * peak) + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/golden/dsp_edit_ita.mat b/tests/golden/dsp_edit_ita.mat new file mode 100644 index 0000000..329c698 Binary files /dev/null and b/tests/golden/dsp_edit_ita.mat differ diff --git a/tests/golden/make_dsp_edit_ita.m b/tests/golden/make_dsp_edit_ita.m new file mode 100644 index 0000000..340df19 --- /dev/null +++ b/tests/golden/make_dsp_edit_ita.m @@ -0,0 +1,22 @@ +% Golden data for tests/test_dsp_edit.py: ita_start_IR, ita_envelope and +% ita_quantize of the ITA-Toolbox on synthetic signals. Indices are 1-based. +% Run from tests/golden with the ITA-Toolbox on the path. +fs = 44100; n = 2048; +rng(5, 'twister'); +ir = zeros(n, 3); +delays = [37 300 900]; +for c = 1:3 + t = (0:n-delays(c)-1)'; + ir(delays(c)+1:end, c) = exp(-t/200) .* (rand(numel(t),1) - 0.5); + ir(delays(c)-4:delays(c), c) = [0.001 -0.002 0.003 -0.004 0.005]'; % weak lead-in +end +ir = double(single(ir + 1e-5 * (rand(n, 3) - 0.5))); +a = itaAudio(ir, fs, 'time'); +start = ita_start_IR(a); +t = (0:n-1)'/fs; +am = (1 + 0.5*cos(2*pi*50*t)) .* cos(2*pi*2000*t); +env = ita_envelope(itaAudio(am, fs, 'time')).timeData; +q = ita_quantize(itaAudio(ir(:,1), fs, 'time'), 'bits', 8).timeData; +ir = single(ir); env = single(env); q = single(q); % am is rebuilt by the test +save('dsp_edit_ita.mat', 'fs', 'ir', 'start', 'env', 'q'); +fprintf('start %s\n', mat2str(start)); diff --git a/tests/test_dsp_edit.py b/tests/test_dsp_edit.py new file mode 100644 index 0000000..4da43a5 --- /dev/null +++ b/tests/test_dsp_edit.py @@ -0,0 +1,135 @@ +""" +Tests for the editing utilities of pytta.functions that follow the ITA-Toolbox +(ita_envelope, ita_start_IR, ita_time_shift, ita_append, ita_normalize_dat, +ita_normalize_spk, ita_amplify_to, ita_zerophase, ita_randomize_phase, +ita_quantize). Each test pins one rule. + +The golden data come from golden/make_dsp_edit_ita.m (MATLAB R2026b with the +ITA-Toolbox); its indices are 1-based. + +SignalObj inspects the caller's frame (classes/_instanceinfo.py) and fails +when the caller's first local variable is an array, so each test starts +with a scalar local. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import pytta + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden') + + +def golden(): + return sio.loadmat(os.path.join(GOLDEN, 'dsp_edit_ita.mat'), squeeze_me=True) + + +class TestEdit(unittest.TestCase): + + def test_start_ir_matches_ita(self): + # ISO 3382-1, A.3: the last sample 20 dB below the peak before it, + # here after a weak lead-in; ITA indices are 1-based. + label = 'start' + d = golden() + ir = pytta.SignalObj(d['ir'].astype(float), 'time', int(d['fs'])) + start = pytta.start_ir(ir) + print(f'\n start {start}, ITA {d["start"] - 1}') + np.testing.assert_array_equal(start, d['start'] - 1) + + def test_envelope_of_an_am_tone(self): + # The envelope of (1 + m cos(2 pi fm t)) cos(2 pi fc t) is + # 1 + m cos(2 pi fm t) away from the ends, and the ITA envelope of + # the same samples agrees. + label = 'envelope' + d = golden() + fs = int(d['fs']) + t = np.arange(2048) / fs + am = (1 + 0.5 * np.cos(2 * np.pi * 50 * t)) * np.cos(2 * np.pi * 2000 * t) + env = pytta.envelope(pytta.SignalObj(am, 'time', fs)).timeSignal[:, 0] + mid = slice(200, -200) + e_rule = np.abs(env[mid] - (1 + 0.5 * np.cos(2 * np.pi * 50 * t[mid]))).max() + e_ita = np.abs(env - d['env']).max() + print(f'\n max |env - (1 + m cos)| {e_rule:.1e} (away from the ends), against ITA {e_ita:.1e}') + self.assertLess(e_rule, 1e-2) + self.assertLess(e_ita, 1e-6) + + def test_time_shift(self): + # An integer shift is a cyclic roll; half a sample twice is one + # sample; None brings the start of the response to sample 0. + label = 'shift' + fs = 48000 + x = np.zeros(256) + x[40:60] = np.hanning(20) + sig = pytta.SignalObj(x, 'time', fs) + rolled = pytta.time_shift(sig, 5, unit='samples').timeSignal[:, 0] + half = pytta.time_shift(pytta.time_shift(sig, 0.5, 'samples'), 0.5, 'samples').timeSignal[:, 0] + auto = pytta.time_shift(sig) + e_half = np.abs(half - np.roll(x, 1)).max() + print(f'\n roll {np.abs(rolled - np.roll(x, 5)).max():.1e}, two half samples {e_half:.1e}, ' + f'start after auto {pytta.start_ir(auto)[0]}') + np.testing.assert_allclose(rolled, np.roll(x, 5), atol=1e-7) + self.assertLess(e_half, 1e-6) + self.assertEqual(pytta.start_ir(auto)[0], 0) + + def test_levels(self): + # normalize brings the peak (time or spectrum) to 1 and returns it; + # amplify_to sets the mean RMS level re 20 uPa for channels in Pa. + label = 'levels' + fs = 48000 + t = np.arange(fs) / fs + x = np.column_stack([0.2 * np.sin(2 * np.pi * 1000 * t), 0.05 * np.sin(2 * np.pi * 500 * t)]) + sig = pytta.SignalObj(x, 'time', fs) + for num in sig.channels.mapping: + sig.channels[num].unit = 'Pa' + normed, gain = pytta.normalize(sig) + per, gains = pytta.normalize(sig, per_channel=True) + spec, _ = pytta.normalize(sig, domain='freq') + at70 = pytta.amplify_to(sig, 70).timeSignal + level = 20 * np.log10(np.mean(np.sqrt(np.mean(at70 ** 2, axis=0))) / 2e-5) + print(f'\n peaks {np.abs(normed.timeSignal).max():.6f}, {np.abs(per.timeSignal).max(axis=0)}, ' + f'gain {gain[0]:.4f}; level {level:.6f} dB') + self.assertAlmostEqual(np.abs(normed.timeSignal).max(), 1, places=6) + np.testing.assert_allclose(np.abs(per.timeSignal).max(axis=0), 1, atol=1e-6) + np.testing.assert_allclose(gains, [0.2, 0.05], rtol=1e-6) + self.assertAlmostEqual(np.abs(np.fft.rfft(spec.timeSignal, axis=0)).max(), 1, places=5) + self.assertAlmostEqual(level, 70, delta=1e-3) # SignalObj stores float32 + + def test_phase_edits_keep_the_magnitude(self): + # zero_phase and randomize_phase keep the magnitude spectrum; the + # zero-phase signal is even (x[n] = x[N - n]). + label = 'phase' + fs = 48000 + x = np.random.default_rng(3).standard_normal(1024) + sig = pytta.SignalObj(x, 'time', fs) + mag = np.abs(np.fft.rfft(x)) + z = pytta.zero_phase(sig).timeSignal[:, 0] + r = pytta.randomize_phase(sig, seed=1).timeSignal[:, 0] + e_z = np.abs(np.abs(np.fft.rfft(z)) - mag).max() + e_r = np.abs(np.abs(np.fft.rfft(r)) - mag).max() + print(f'\n magnitude kept: zero phase {e_z:.1e}, random phase {e_r:.1e}') + self.assertLess(max(e_z, e_r) / mag.max(), 1e-5) + np.testing.assert_allclose(z[1:], z[1:][::-1], atol=1e-6) + + def test_append_and_quantize(self): + # append concatenates in time; quantize rounds to 2^(b-1) + 1 steps + # across the peak, as ITA does with the same samples. + label = 'append' + d = golden() + fs = int(d['fs']) + a = pytta.SignalObj(np.ones(100), 'time', fs) + b = pytta.SignalObj(np.zeros(50), 'time', fs) + self.assertEqual(pytta.append(a, b).numSamples, 150) + q = pytta.quantize(pytta.SignalObj(d['ir'][:, 0].astype(float), 'time', fs), bits=8) + e = np.abs(q.timeSignal[:, 0] - d['q']).max() + print(f'\n quantize against ITA {e:.1e}') + self.assertLess(e, 1e-7) + with self.assertRaises(ValueError): + pytta.append(a, pytta.SignalObj(np.zeros(50), 'time', 48000)) + + +if __name__ == '__main__': + unittest.main()