From 1db0deaf33d06a44715ef9d53b1f6c11fd42f57d Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 01:49:38 +0100 Subject: [PATCH] [FEATURE] add warp, warp_zpk and bark_warping for frequency warping Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/functions.rst | 4 ++- pytta/__init__.py | 6 ++++- pytta/functions.py | 51 +++++++++++++++++++++++++++++++++++++ tests/test_warping.py | 46 +++++++++++++++++++++++++++++++++ 4 files changed, 105 insertions(+), 2 deletions(-) create mode 100644 tests/test_warping.py diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..df53962 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.bark_warping +.. autofunction:: pytta.warp +.. autofunction:: pytta.warp_zpk diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..d0b6e61 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,\ + bark_warping, warp, warp_zpk # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +147,9 @@ 'resample', 'corr_coef', 'peak_time', + 'bark_warping', + 'warp', + 'warp_zpk', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..868f998 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -274,6 +274,57 @@ def peak_time(signal): else: return peaks_time[0] +def bark_warping(samplingRate): + """ + Warping coefficient that makes the frequency axis approach the Bark + scale, 0.8517 sqrt(arctan(0.06583 fs / 1000)) - 0.1916 (Smith and Abel, + Bark and ERB bilinear transforms, IEEE Trans. Speech Audio Process. 7(6), + 1999), as the barkwarp of the ITA-Toolbox. + """ + return 0.8517 * np.sqrt(np.arctan(0.06583 * samplingRate / 1000)) - 0.1916 + + +def warp(signal, lam, numSamples=None): + """ + Frequency warped signal, as ita_audio_warp: every delay z^-1 of the + signal, taken as an FIR filter, is replaced by the allpass + (z^-1 - lam) / (1 - lam z^-1), so that the new spectrum at omega is the + old one at omega + 2 arctan(lam sin(omega) / (1 - lam cos(omega))). + Positive lam stretches the low frequencies; warp with -lam undoes it. + The result has numSamples samples (those of signal if None). + + >>> w = pytta.warp(ir, pytta.bark_warping(ir.samplingRate)) + """ + x = np.asarray(signal.timeSignal, dtype=float) + n = numSamples or x.shape[0] + y = np.zeros((n, x.shape[1])) + g = np.zeros(n) + g[0] = 1 + for k in range(x.shape[0]): + y += g[:, None] * x[k] + g = ss.lfilter([-lam, 1], [1, -lam], g) + return SignalObj(y, 'time', signal.samplingRate) + + +def warp_zpk(z, p, k, lam): + """ + Zeros, poles and gain of the filter H(D(z)) with z^-1 replaced by the + allpass D(z) = (z^-1 - lam) / (1 - lam z^-1), as ita_zpk_warp: each root + r moves to (r + lam) / (1 + lam r), and the difference between the + numbers of zeros and poles leaves roots at lam. warp_zpk with -lam undoes + it. + """ + z, p = np.atleast_1d(z).astype(complex), np.atleast_1d(p).astype(complex) + k = k * np.prod(1 + lam * z) / np.prod(1 + lam * p) + extra = np.full(abs(z.size - p.size), lam, dtype=complex) + zw, pw = (z + lam) / (1 + lam * z), (p + lam) / (1 + lam * p) + if z.size > p.size: + pw = np.r_[pw, extra] + else: + zw = np.r_[zw, extra] + return zw, pw, k.real if abs(k.imag) <= 1e-12 * abs(k) else k + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/test_warping.py b/tests/test_warping.py new file mode 100644 index 0000000..e3c3427 --- /dev/null +++ b/tests/test_warping.py @@ -0,0 +1,46 @@ +""" +Tests for warp and warp_zpk of pytta.functions, which follow ita_audio_warp +and ita_zpk_warp of the ITA-Toolbox. Each test compares with the analytic +frequency mapping of the allpass D(z) = (z^-1 - lam) / (1 - lam z^-1). +""" +import unittest + +import numpy as np +import scipy.signal as ss + +import pytta + +LAM = 0.6 + + +def theta(w, lam=LAM): + # Phase of the allpass: D(e^jw) = e^(-j theta(w)). + return w + 2 * np.arctan2(lam * np.sin(w), 1 - lam * np.cos(w)) + + +class TestWarping(unittest.TestCase): + + def test_warped_signal_has_the_mapped_spectrum(self): + label = 'x' + h = np.random.default_rng(5).standard_normal(64) * 0.9 ** np.arange(64) + y = pytta.warp(pytta.SignalObj(h, 'time', 48000), LAM, 4096).timeSignal[:, 0].astype(float) + w = np.linspace(0.05, 3, 25) + Y = np.exp(-1j * np.outer(w, np.arange(y.size))) @ y + H = np.exp(-1j * np.outer(theta(w), np.arange(h.size))) @ h + e = np.abs(Y - H).max() / np.abs(H).max() + print(f'\n warped spectrum: {e:.1e}', end='') + self.assertLess(e, 1e-6) + + def test_warped_filter_has_the_mapped_response_and_unwarps(self): + z, p, k = ss.ellip(4, 1, 40, 0.2, output='zpk') + zw, pw, kw = pytta.warp_zpk(z, p, k, LAM) + w = np.linspace(0.05, 3, 25) + e = np.abs(ss.freqz_zpk(zw, pw, kw, w)[1] - ss.freqz_zpk(z, p, k, theta(w))[1]).max() + zb, pb, kb = pytta.warp_zpk(zw, pw, kw, -LAM) + back = np.abs(ss.freqz_zpk(zb, pb, kb, w)[1] - ss.freqz_zpk(z, p, k, w)[1]).max() + print(f'\n warped filter: {e:.1e}, unwarped {back:.1e}', end='') + self.assertLess(max(e, back), 1e-10) + + +if __name__ == '__main__': + unittest.main()