[FEATURE] add time_window for SignalObj

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',

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

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

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()