diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..af4648e 100644 --- a/docs/pyttadoc/functions.rst +++ b/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 diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..b236377 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,\ + 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', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..5ab1c08 100644 --- a/pytta/functions.py +++ b/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 diff --git a/pytta/generate.py b/pytta/generate.py index 1979f54..b20a51a 100644 --- a/pytta/generate.py +++ b/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, diff --git a/tests/golden/cut_impulse_train_ita.mat b/tests/golden/cut_impulse_train_ita.mat new file mode 100644 index 0000000..d59e704 Binary files /dev/null and b/tests/golden/cut_impulse_train_ita.mat differ diff --git a/tests/golden/make_cut_impulse_train_ita.m b/tests/golden/make_cut_impulse_train_ita.m new file mode 100644 index 0000000..b51ddef --- /dev/null +++ b/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))); diff --git a/tests/test_dsp_misc.py b/tests/test_dsp_misc.py new file mode 100644 index 0000000..685d398 --- /dev/null +++ b/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()