Add multitone, spectrum crossfade, impulse train cutting and velocity from a p-p probe

docs/pyttadoc/functions.rst

@@ -26,4 +26,6 @@ Functions
 .. autofunction:: pytta.plot_bars
 .. autofunction:: pytta.plot_spectrogram
 .. autofunction:: pytta.plot_waterfall
-
+.. autofunction:: pytta.crossfade_spectra
+.. autofunction:: pytta.cut_impulse_train
+.. autofunction:: pytta.velocity_from_pressure_gradient

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,\
+    crossfade_spectra, cut_impulse_train, velocity_from_pressure_gradient
 
 # Must go on v0.1.0 (DEPRECATED)
 from .functions import fft_degree
@@ -146,6 +147,9 @@
            'resample',
            'corr_coef',
            'peak_time',
+           'crossfade_spectra',
+           'cut_impulse_train',
+           'velocity_from_pressure_gradient',
            'save',
            'load',
            'weighting',

pytta/functions.py

@@ -274,6 +274,92 @@ def peak_time(signal):
     else:
         return peaks_time[0]
 
+def crossfade_spectra(signal_low, signal_high, freqs):
+    """
+    Spectrum of signal_low below freqs[0] Hz and of signal_high above
+    freqs[1] Hz, with a raised cosine crossfade between them, as
+    ita_xfade_spk; a single frequency takes half of each at that bin.
+    Both signals need the same sampling rate and number of samples.
+
+    >>> combined = pytta.crossfade_spectra(ir_simulated, ir_measured, [200, 400])
+    """
+    fs = signal_low.samplingRate
+    if fs != signal_high.samplingRate:
+        raise ValueError('The signals must have the same sampling rate.')
+    f0, f1 = np.broadcast_to(np.asarray(freqs, dtype=float), (2,))
+    n = signal_low.timeSignal.shape[0]
+    A = np.fft.rfft(signal_low.timeSignal, axis=0)
+    if f0 >= fs / 2:
+        return SignalObj(np.fft.irfft(A, n, axis=0), 'time', fs)
+    B = np.fft.rfft(signal_high.timeSignal, axis=0)
+    b0, b1 = (int(round(f * n / fs)) for f in (f0, min(f1, fs / 2)))
+    if b1 > b0:
+        rise = 0.5 - 0.5 * np.cos(np.pi * np.arange(b1 - b0 + 1) / (b1 - b0))
+    else:
+        rise = np.array([0.5])
+    gain_low = np.zeros(A.shape[0])
+    gain_low[:b0] = 1
+    gain_low[b0:b1 + 1] = 1 - rise
+    gain_high = np.zeros(A.shape[0])
+    gain_high[b1 + 1:] = 1
+    gain_high[b0:b1 + 1] = rise
+    X = A * gain_low[:, None] + B * gain_high[:, None]
+    return SignalObj(np.fft.irfft(X, n, axis=0), 'time', fs)
+
+
+def cut_impulse_train(signal, threshold=0.7, length=1, start_before=0.004, ref_channel=0):
+    """
+    Impulses of an impulse train cut out of signal, as ita_cutImpulseTrain.
+    A sample of ref_channel above the threshold, threshold * L_max +
+    (1 - threshold) * L_rms dB, begins an impulse when it comes length
+    seconds or more after the previous sample above it; each impulse is
+    length seconds from start_before seconds before that sample, cut at the
+    same samples in every channel.
+
+    Returns a list with one SignalObj per channel, whose channels are the
+    impulses, and the 0-based start samples.
+
+    >>> impulses, starts = pytta.cut_impulse_train(take, length=0.5)
+    """
+    fs = signal.samplingRate
+    x = np.asarray(signal.timeSignal)
+    n = x.shape[0]
+    nlen, nbefore = round(length * fs), round(start_before * fs)
+    ref = x[:, ref_channel]
+    with np.errstate(divide='ignore'):
+        level = 20 * np.log10(np.abs(ref))
+    thr = threshold * level.max() + (1 - threshold) * 10 * np.log10(np.mean(ref ** 2))
+    idx = np.flatnonzero(level > thr)
+    idx = idx[idx <= n - nlen + nbefore]
+    idx = idx[np.r_[True, np.diff(idx) >= nlen]] if idx.size else idx
+    starts = idx[idx >= nbefore] - nbefore
+    take = starts[:, None] + np.arange(nlen)
+    return [SignalObj(x[take, c].T, 'time', fs) for c in range(x.shape[1])], starts
+
+
+def velocity_from_pressure_gradient(signal, distance, rho=1.2, c=343, normalized=False):
+    """
+    Particle velocity (m/s) midway between two microphones distance m apart
+    from the pressures of channels 0 and 1 (Pa), by Euler's equation,
+    v = (p0 - p1) / (j omega rho distance), as
+    ita_velocity_from_pressure_gradient; times rho c if normalized. The DC
+    bin is set to zero.
+
+    >>> v = pytta.velocity_from_pressure_gradient(pp_probe, 0.012)
+    """
+    x = np.asarray(signal.timeSignal, dtype=float)
+    if x.shape[1] != 2:
+        raise ValueError('The signal must have two channels.')
+    n = x.shape[0]
+    P = np.fft.rfft(x[:, 0] - x[:, 1])
+    omega = 2 * np.pi * np.fft.rfftfreq(n, 1 / signal.samplingRate)
+    V = np.zeros_like(P)
+    V[1:] = P[1:] / (1j * omega[1:] * rho * distance)
+    if normalized:
+        V *= rho * c
+    return SignalObj(np.fft.irfft(V, n), 'time', signal.samplingRate)
+
+
 def fft_degree(*args,**kwargs):
     """
     DEPRECATED

pytta/generate.py

@@ -504,6 +504,29 @@ def impulse(samplingRate=None,
     return impulseSignal
 
 
+def multitone(freqs, amplitude=1, phase=0, samplingRate=None, fftDegree=None):
+    """
+    Sum of sines of peak amplitude amplitude and phase phase (radians, one
+    for all tones or one per tone) at the frequencies freqs rounded to the
+    nearest DFT bin of 2**fftDegree samples, so that every tone has a whole
+    number of periods, as ita_multitone. Logarithmic or fixed factor spacing
+    is freqs = numpy.geomspace(...).
+
+    >>> sig = pytta.generate.multitone(np.geomspace(100, 10000, 5))
+    """
+    if samplingRate is None:
+        samplingRate = default.samplingRate
+    if fftDegree is None:
+        fftDegree = default.fftDegree
+    n = 2**fftDegree
+    f = np.round(np.asarray(freqs, dtype=float) * n / samplingRate) * samplingRate / n
+    if f.max() > samplingRate / 2:
+        raise ValueError('The frequencies must not exceed samplingRate/2.')
+    t = np.arange(n)[:, None] / samplingRate
+    x = (amplitude * np.sin(2 * np.pi * f * t + np.broadcast_to(phase, f.shape))).sum(axis=1)
+    return SignalObj(x, 'time', samplingRate, signalType='power')
+
+
 def measurement(kind='playrec',
                 samplingRate=None,
                 freqMin=None,

tests/golden/make_cut_impulse_train_ita.m

@@ -0,0 +1,20 @@
+% Golden data for tests/test_dsp_misc.py: ita_cutImpulseTrain of the
+% ITA-Toolbox on a synthetic train of four decaying bursts in two channels,
+% some with an earlier weak arrival. Indices are 1-based.
+% Run from tests/golden with the ITA-Toolbox on the path.
+fs = 8000; n = 12000;
+rng(7, 'twister');
+x = 1e-3 * (rand(n, 2) - 0.5);
+onsets = [800 3700 6100 9400];
+for k = 1:numel(onsets)
+    t = (0:1499)';
+    burst = exp(-t/150) .* (rand(1500, 1) - 0.5);
+    x(onsets(k):onsets(k)+1499, 1) = x(onsets(k):onsets(k)+1499, 1) + burst;
+    x(onsets(k)+20:onsets(k)+1519, 2) = x(onsets(k)+20:onsets(k)+1519, 2) + 0.5*burst;
+end
+x = double(single(x));
+[cuts, ~, edges] = ita_cutImpulseTrain(itaAudio(x, fs, 'time'), 'length', 0.25, 'startBeforeThreshold', 0.002);
+start = edges(:, 1);
+c1 = single(cuts(1).timeData); c2 = single(cuts(2).timeData); x = single(x);
+save('cut_impulse_train_ita.mat', 'fs', 'x', 'start', 'c1', 'c2');
+fprintf('start %s, size %s\n', mat2str(start'), mat2str(size(c1)));

tests/test_dsp_misc.py

@@ -0,0 +1,87 @@
+"""
+Tests for pytta.generate.multitone and the functions crossfade_spectra,
+cut_impulse_train and velocity_from_pressure_gradient of pytta.functions,
+which follow the ITA-Toolbox (ita_multitone, ita_xfade_spk,
+ita_cutImpulseTrain, ita_velocity_from_pressure_gradient). Each test pins
+one rule against an analytic result or golden data of the ITA-Toolbox from
+tests/golden/make_cut_impulse_train_ita.m (MATLAB R2026b).
+"""
+import os
+import unittest
+
+import numpy as np
+from scipy.io import loadmat
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import pytta
+
+GOLDEN = os.path.join(os.path.dirname(__file__), 'golden')
+
+
+def tone(f, fs, n, phase=0):
+    return np.cos(2 * np.pi * f * np.arange(n) / fs + phase)
+
+
+class TestDSPMisc(unittest.TestCase):
+
+    def test_multitone_tones_fall_on_bins(self):
+        # Every tone has a whole number of periods, so the DFT is A N / 2 at
+        # the rounded bins and zero elsewhere.
+        label = 'x'
+        fs, deg = 48000, 14
+        n = 2 ** deg
+        sig = pytta.generate.multitone([100, 316.2, 1000, 3162, 10000], 0.3, samplingRate=fs, fftDegree=deg)
+        X = np.abs(np.fft.rfft(sig.timeSignal[:, 0].astype(float)))
+        bins = np.round(np.array([100, 316.2, 1000, 3162, 10000]) * n / fs).astype(int)
+        e_on = np.abs(X[bins] - 0.3 * n / 2).max() / (0.3 * n / 2)
+        e_off = np.delete(X, bins).max() / (0.3 * n / 2)
+        print(f'\n  multitone: on bins {e_on:.1e}, off bins {e_off:.1e}', end='')
+        self.assertLess(e_on, 1e-6)
+        self.assertLess(e_off, 1e-6)
+
+    def test_crossfade_keeps_each_band_and_sums_to_one(self):
+        # Below f0 the low signal, above f1 the high one, and the raised
+        # cosine weights add up to 1, so a signal crossfaded with itself
+        # comes back.
+        label = 'x'
+        fs, n = 8000, 4096
+        rng = np.random.default_rng(3)
+        a = pytta.SignalObj(rng.standard_normal(n), 'time', fs)
+        b = pytta.SignalObj(rng.standard_normal(n), 'time', fs)
+        A, B = (np.fft.rfft(s.timeSignal[:, 0].astype(float)) for s in (a, b))
+        X = np.fft.rfft(pytta.crossfade_spectra(a, b, [500, 1000]).timeSignal[:, 0].astype(float))
+        b0, b1 = 500 * n // fs, 1000 * n // fs
+        e_low = np.abs(X[:b0] - A[:b0]).max() / np.abs(A).max()
+        e_high = np.abs(X[b1 + 1:] - B[b1 + 1:]).max() / np.abs(B).max()
+        same = pytta.crossfade_spectra(a, a, [500, 1000]).timeSignal - a.timeSignal
+        print(f'\n  crossfade: low band {e_low:.1e}, high band {e_high:.1e}, '
+              f'with itself {np.abs(same).max():.1e}', end='')
+        self.assertLess(max(e_low, e_high), 1e-6)
+        self.assertLess(np.abs(same).max(), 1e-6)
+
+    def test_cut_impulse_train_matches_ita(self):
+        g = loadmat(os.path.join(GOLDEN, 'cut_impulse_train_ita.mat'))
+        sig = pytta.SignalObj(g['x'].astype(float), 'time', g['fs'].item())
+        cuts, starts = pytta.cut_impulse_train(sig, length=0.25, start_before=0.002)
+        print(f'\n  cut_impulse_train: starts {starts}, ITA {g["start"].ravel() - 1}', end='')
+        np.testing.assert_array_equal(starts, g['start'].ravel() - 1)
+        np.testing.assert_array_equal(cuts[0].timeSignal, g['c1'])
+        np.testing.assert_array_equal(cuts[1].timeSignal, g['c2'])
+
+    def test_pressure_gradient_of_a_plane_wave(self):
+        # For a plane wave the finite difference gives the velocity at the
+        # midpoint, p / (rho c), times sin(kd/2) / (kd/2).
+        label = 'x'
+        fs, n, f, d, rho, c = 48000, 4800, 1000, 0.05, 1.2, 343
+        k = 2 * np.pi * f / c
+        p = pytta.SignalObj(np.column_stack([tone(f, fs, n), tone(f, fs, n, -k * d)]), 'time', fs)
+        v = pytta.velocity_from_pressure_gradient(p, d, rho, c).timeSignal[:, 0]
+        expected = np.sinc(k * d / 2 / np.pi) * tone(f, fs, n, -k * d / 2) / (rho * c)
+        e = np.abs(v - expected).max() * rho * c
+        print(f'\n  pressure gradient: {e:.1e} of p / (rho c)', end='')
+        self.assertLess(e, 1e-6)
+
+
+if __name__ == '__main__':
+    unittest.main()