Add frequency warping of signals and filters

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

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

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

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