diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..aa332ae 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -113,7 +113,7 @@ 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, time_window # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +146,7 @@ 'resample', 'corr_coef', 'peak_time', + 'time_window', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..579dd64 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -17,6 +17,7 @@ >>> pytta.corr_coef(signalObj1, signalObj2) >>> pytta.resample(signalObj, newSamplingRate) >>> pytta.peak_time(signalObj1, signalObj2, ..., signalObjN) + >>> pytta.time_window(signalObj, interval) >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN) @@ -274,6 +275,87 @@ def peak_time(signal): else: return peaks_time[0] +def time_window(signalObj, interval, unit='time', window='hann', + symmetric=False): + """ + Apply a time window to a SignalObj and return a new SignalObj. + + The window is built from pairs of bounds. Each pair [start, end] with + start <= end is a falling slope: ones before start, the falling half of + the window between start and end (both included), zeros after end. A pair + with start > end is the mirror image, a rising slope: zeros before end, + rising half between end and start, ones after start. Several pairs + multiply. The usual impulse response window is therefore + [startFade, 0, endFadeStart, endFadeEnd]: a rise at the beginning and a + fall at the end. + + Parameters (default), (type): + ----------------------------- + * signalObj (), (SignalObj): + signal to be windowed. It is not modified; + * interval (), (list): + an even number of bounds, in seconds or samples; + * unit ('time'), (str): + 'time' for seconds (converted with round(t * fs), halves up) or + 'samples' for integer, zero based sample indices; + * window ('hann'), (str): + 'hann', 'rect' or any window name scipy.signal.get_window + accepts without parameters. The slope has abs(end - start) + 1 + samples and is one half of a symmetric window twice that long, + so a 'hann' slope ends at exactly zero; + * symmetric (False), (bool): + mirror the window about the centre of the signal, for impulse + responses that are slightly non causal. + + The behaviour follows ita_time_window of the ITA-Toolbox (RWTH Aachen), + except that the bounds are not moved to even sample numbers. + + >>> ir = pytta.time_window(ir, [0.02, 0, 0.8, 1.0]) + """ + if not isinstance(signalObj, SignalObj): + raise TypeError('signalObj must be a SignalObj.') + bounds = np.asarray(interval, dtype=float).ravel() + if bounds.size == 0 or bounds.size % 2: + raise ValueError('interval must hold an even number of bounds.') + if unit == 'time': + bounds = np.floor(bounds * signalObj.samplingRate + 0.5) + elif unit != 'samples': + raise ValueError("unit must be 'time' or 'samples'.") + numSamples = signalObj.numSamples + if np.any(bounds < 0) or np.any(bounds >= numSamples): + raise ValueError('interval bounds must lie within the signal.') + bounds = bounds.astype(int) + name = {'rect': 'boxcar', 'rectwin': 'boxcar'}.get(window, window) + try: + ss.get_window(name, 4, fftbins=False) + except ValueError: + raise ValueError(f"unknown window '{window}'.") + fullWindow = np.ones(numSamples) + for start, end in bounds.reshape(-1, 2): + length = abs(end - start) + 1 + halves = ss.get_window(name, 2 * length, fftbins=False) + slope = np.ones(numSamples) + if start <= end: + slope[start:end + 1] = halves[length:] + slope[end + 1:] = 0 + else: + slope[:end] = 0 + slope[end:start + 1] = halves[:length] + fullWindow *= slope + if symmetric: + fullWindow = fullWindow + fullWindow[::-1] + fullWindow -= fullWindow.min() + newSignal = SignalObj(signalObj.timeSignal * fullWindow[:, np.newaxis], + 'time', + samplingRate=signalObj.samplingRate, + signalType=signalObj.signalType, + freqMin=signalObj.freqMin, + freqMax=signalObj.freqMax, + comment=signalObj.comment) + newSignal.channels = cp.deepcopy(signalObj.channels) + return newSignal + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/golden/make_time_window_ita.m b/tests/golden/make_time_window_ita.m new file mode 100644 index 0000000..61a8395 --- /dev/null +++ b/tests/golden/make_time_window_ita.m @@ -0,0 +1,49 @@ +function make_time_window_ita(itaRoot) +% MAKE_TIME_WINDOW_ITA Golden data for tests/test_time_window.py. +% +% make_time_window_ita(itaRoot) runs ita_time_window from the ITA-Toolbox +% (RWTH Aachen, https://git.rwth-aachen.de/ita/toolbox) on a fixed random +% signal and writes time_window_ita.mat next to this file. +% +% Intervals are given in SAMPLES, 1-based as MATLAB indexes them, so the +% file pins the window shape itself and not the seconds-to-samples +% rounding of either implementation. The Python test subtracts 1. +% +% The signal is stored in the file: the test never regenerates it. +% +% Usage: matlab -batch "make_time_window_ita('/path/to/ITA-Toolbox')" + +addpath(genpath(itaRoot)); + +fs = 1000; +rng(0, 'twister'); +x = randn(1000, 2); +a = itaAudio(x, fs, 'time'); + +% name, interval (1-based samples), window, symmetric +spec = { + 'falling_hann', [201 400], @hann, false + 'rising_hann', [600 450], @hann, false + 'both_ends_hann', [100 1 800 1000], @hann, false + 'both_ends_rect', [50 1 900 1000], @rectwin, false + 'symmetric_hann', [301 500], @hann, true +}; + +cases = struct('name', {}, 'interval', {}, 'window', {}, 'symmetric', {}, 'y', {}); +for k = 1:size(spec, 1) + if spec{k, 4} + r = ita_time_window(a, spec{k, 2}, spec{k, 3}, 'samples', 'symmetric'); + else + r = ita_time_window(a, spec{k, 2}, spec{k, 3}, 'samples'); + end + cases(k).name = spec{k, 1}; + cases(k).interval = spec{k, 2}; + cases(k).window = func2str(spec{k, 3}); + cases(k).symmetric = spec{k, 4}; + cases(k).y = r.timeData; +end + +out = fullfile(fileparts(mfilename('fullpath')), 'time_window_ita.mat'); +save(out, 'x', 'fs', 'cases', '-v7'); +fprintf('wrote %s (%d cases)\n', out, numel(cases)); +end diff --git a/tests/golden/time_window_ita.mat b/tests/golden/time_window_ita.mat new file mode 100644 index 0000000..9c61a1a Binary files /dev/null and b/tests/golden/time_window_ita.mat differ diff --git a/tests/test_time_window.py b/tests/test_time_window.py new file mode 100644 index 0000000..4246b0a --- /dev/null +++ b/tests/test_time_window.py @@ -0,0 +1,144 @@ +""" +Tests for pytta.time_window. + +Each test pins one rule of the function. The shape rules are checked on a +signal of ones, so the output is the window itself. The last test compares +the whole thing with the ITA-Toolbox (RWTH Aachen) on data generated by +golden/make_time_window_ita.m. + +Index convention: an interval is a list of [start, end] pairs. start <= end +is a falling slope (ones before start, zeros after end); start > end is a +rising slope (zeros before end, ones after start). Both bounds belong to the +slope, which therefore has abs(end - start) + 1 samples. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +import pytta + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden', 'time_window_ita.mat') +FS = 1000 +N = 1000 +TOL32 = 1e-6 # rounding of the float32 SignalObj storage + + +def ones(numChannels=1): + return pytta.SignalObj(np.ones((N, numChannels)), 'time', FS) + + +def window_of(interval, **kwargs): + """Apply time_window to ones and return the resulting window, 1-D.""" + out = pytta.time_window(ones(), interval, unit='samples', **kwargs) + return out.timeSignal[:, 0].astype(float) + + +class TestTimeWindow(unittest.TestCase): + + def test_falling_slope_is_hann_half_between_bounds(self): + # [200, 399]: 200 ones, then the falling half of hann(400), i.e. + # 0.5 * (1 - cos(2 pi k / 399)) for k = 200..399, then zeros. + # The slope ends at exactly 0 (symmetric hann, not hanning). + s, e = 200, 399 + w = window_of([s, e]) + L = e - s + 1 + k = np.arange(L, 2 * L) + expected = 0.5 * (1 - np.cos(2 * np.pi * k / (2 * L - 1))) + worst = np.max(np.abs(w[s:e + 1] - expected)) + print(f'\n falling slope vs analytic hann, worst={worst:.2e} (tol {TOL32})') + self.assertTrue(np.all(w[:s] == 1)) + self.assertLess(worst, TOL32) + self.assertEqual(w[e], 0) + self.assertTrue(np.all(w[e + 1:] == 0)) + + def test_rising_slope_mirrors_falling(self): + # start > end reverses the slope: zeros, rise, ones. It must be the + # sample-by-sample mirror of the falling window on mirrored indices. + rising = window_of([600, 450]) + falling = window_of([N - 1 - 600, N - 1 - 450]) + self.assertTrue(np.array_equal(rising, falling[::-1])) + + def test_pairs_compose_by_product(self): + # Two pairs are two windows multiplied, independent of order. + both = window_of([120, 0, 800, 999]) + rise = window_of([120, 0]) + fall = window_of([800, 999]) + self.assertTrue(np.allclose(both, rise * fall, atol=TOL32)) + self.assertTrue(np.array_equal(both, window_of([800, 999, 120, 0]))) + + def test_rect_window_is_hard_cut(self): + # A rectangular window keeps the slope at 1: a cut right after end. + w = window_of([300, 499], window='rect') + self.assertTrue(np.all(w[:500] == 1)) + self.assertTrue(np.all(w[500:] == 0)) + + def test_time_unit_rounds_to_nearest_sample(self): + # Seconds are converted with round(t * fs), halves up. 0.2004 s at + # 1 kHz is sample 200, 0.3995 s is sample 400 (not 399). No further + # adjustment of the bounds (the ITA-Toolbox moves them to even + # numbers; this implementation does not). + inSeconds = pytta.time_window(ones(), [0.2004, 0.3995], unit='time') + inSamples = pytta.time_window(ones(), [200, 400], unit='samples') + self.assertTrue(np.array_equal(inSeconds.timeSignal, inSamples.timeSignal)) + + def test_symmetric_mirrors_window_about_centre(self): + # symmetric=True adds the window to its reverse and removes the + # offset: for a falling slope in the first half it yields + # ones, fall, zeros, rise, ones, symmetric about the centre. + s, e = 300, 499 + w = window_of([s, e], symmetric=True) + self.assertTrue(np.array_equal(w, w[::-1])) + self.assertTrue(np.all(w[:s] == 1)) + self.assertEqual(w[e], 0) + self.assertEqual(w[N // 2], 0) + + def test_returns_new_object_and_keeps_metadata(self): + # The input is not modified; the output keeps sampling rate and + # channel count and windows every channel the same way. + sig = pytta.SignalObj(np.random.default_rng(0).standard_normal((N, 3)), + 'time', FS) + before = sig.timeSignal.copy() + out = pytta.time_window(sig, [0.1, 0.4]) + self.assertTrue(np.array_equal(sig.timeSignal, before)) + self.assertIsNot(out, sig) + self.assertEqual(out.samplingRate, FS) + self.assertEqual(out.numChannels, 3) + w = window_of([100, 400]) + for ch in range(3): + self.assertTrue(np.allclose(out.timeSignal[:, ch], + before[:, ch] * w, atol=TOL32)) + + def test_rejects_invalid_interval(self): + with self.assertRaises(ValueError): + pytta.time_window(ones(), [100, 200, 300], unit='samples') # odd + with self.assertRaises(ValueError): + pytta.time_window(ones(), [100, N], unit='samples') # past end + with self.assertRaises(ValueError): + pytta.time_window(ones(), [0.1, 0.2], unit='seconds') # unknown + with self.assertRaises(ValueError): + pytta.time_window(ones(), [100, 200], unit='samples', window='no') + + def test_matches_ita_toolbox(self): + # Sample-by-sample against ita_time_window on the stored signal, all + # five cases of golden/make_time_window_ita.m. Intervals there are + # 1-based MATLAB samples; subtract 1. Tolerance is the float32 + # storage of SignalObj relative to the signal peak. + self.assertTrue(os.path.exists(GOLDEN), 'run golden/make_time_window_ita.m') + g = sio.loadmat(GOLDEN, squeeze_me=True, struct_as_record=False) + x, fs = g['x'], int(g['fs']) + sig = pytta.SignalObj(x, 'time', fs) + peak = np.max(np.abs(x)) + for case in np.atleast_1d(g['cases']): + interval = (np.atleast_1d(case.interval) - 1).tolist() + window = {'hann': 'hann', 'rectwin': 'rect'}[case.window] + out = pytta.time_window(sig, interval, unit='samples', + window=window, symmetric=bool(case.symmetric)) + diff = np.max(np.abs(out.timeSignal - case.y)) / peak + print(f'\n {case.name:16s} rel peak diff vs ITA = {diff:.2e} (tol {TOL32})') + self.assertLess(diff, TOL32, case.name) + + +if __name__ == '__main__': + unittest.main()