Add conversion between power series and harmonic series of a nonlinearity

docs/pyttadoc/functions.rst

@@ -26,4 +26,5 @@ Functions
 .. autofunction:: pytta.plot_bars
 .. autofunction:: pytta.plot_spectrogram
 .. autofunction:: pytta.plot_waterfall
-
+.. autofunction:: pytta.power_to_harmonic_series
+.. autofunction:: pytta.harmonic_to_power_series

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

pytta/functions.py

@@ -274,6 +274,46 @@ def peak_time(signal):
     else:
         return peaks_time[0]
 
+def _harmonic_matrix(order, gain, phase):
+    # Column k - 1: amplitudes of the harmonics 1 to order that x**k gives for
+    # x = gain cos(wt), divided by gain as the sweep deconvolution does, with
+    # harmonic m turned by (m - 1) phase (Dietrich 2013, sec. 3.3).
+    T = np.zeros((order, order), dtype=complex)
+    for k in range(1, order + 1):
+        T[:k, k - 1] = np.polynomial.chebyshev.poly2cheb(np.eye(k + 1)[k])[1:] * gain ** (k - 1)
+    return T * np.exp(1j * phase * np.arange(order))[:, None]
+
+
+def power_to_harmonic_series(coefficients, gain=1, phase=0):
+    """
+    Harmonic series of a power series, as
+    ita_nonlinear_power_to_harmonic_series: the amplitudes of the harmonics
+    1 to N that y = sum_k a_k x**k gives for a cosine of amplitude gain,
+    divided by gain, from the coefficients a_1 to a_N along the first axis
+    (scalars, or arrays converted entry by entry). Harmonic m is turned by
+    (m - 1) phase, the phase offset of an exponential sweep. The DC term is
+    left out.
+
+    >>> pytta.power_to_harmonic_series([1, 1, 1, 1], gain=0.5)  # 1.1875, 0.3125, 0.0625, 0.015625
+    """
+    a = np.asarray(coefficients)
+    T = _harmonic_matrix(a.shape[0], gain, phase)
+    h = np.tensordot(T, a, axes=1)
+    return h.real if np.isrealobj(a) and phase == 0 else h
+
+
+def harmonic_to_power_series(harmonics, gain=1, phase=0):
+    """
+    Power series coefficients a_1 to a_N of the harmonic series 1 to N, the
+    inverse of power_to_harmonic_series, as
+    ita_nonlinear_harmonic_to_power_series.
+    """
+    h = np.asarray(harmonics)
+    T = _harmonic_matrix(h.shape[0], gain, phase)
+    a = np.tensordot(np.linalg.inv(T), h, axes=1)
+    return a.real if np.isrealobj(h) and phase == 0 else a
+
+
 def fft_degree(*args,**kwargs):
     """
     DEPRECATED

tests/test_nonlinear_series.py

@@ -0,0 +1,46 @@
+"""
+Tests for power_to_harmonic_series and harmonic_to_power_series of
+pytta.functions, which follow ita_nonlinear_power_to_harmonic_series and
+ita_nonlinear_harmonic_to_power_series of the ITA-Toolbox and Dietrich,
+Uncertainties in acoustical transfer functions (2013), sec. 3.3.
+"""
+import unittest
+
+import numpy as np
+
+import pytta
+
+
+class TestNonlinearSeries(unittest.TestCase):
+
+    def test_example_of_dietrich(self):
+        # y = x + x^2 + x^3 + x^4 driven with G = 1/2 (Dietrich 2013, p. 51).
+        h = pytta.power_to_harmonic_series([1, 1, 1, 1], gain=0.5)
+        print(f'\n  harmonics {h}', end='')
+        np.testing.assert_allclose(h, [1.1875, 0.3125, 0.0625, 0.0156], atol=5e-5)
+
+    def test_harmonics_of_a_driven_polynomial(self):
+        # The DFT of the polynomial of gain cos(wt), over one period and
+        # divided by gain, has the predicted harmonics.
+        a = np.random.default_rng(1).standard_normal(6)
+        g, n = 0.8, 64
+        x = g * np.cos(2 * np.pi * np.arange(n) / n)
+        y = sum(a[k] * x ** (k + 1) for k in range(6))
+        dft = 2 * np.fft.rfft(y)[1:7].real / n / g
+        e = np.abs(pytta.power_to_harmonic_series(a, gain=g) - dft).max()
+        print(f'\n  against the DFT: {e:.1e}', end='')
+        self.assertLess(e, 1e-13)
+
+    def test_phase_and_inverse(self):
+        # Harmonic m turns by (m - 1) phase, as in the ITA-Toolbox, and the
+        # inverse gives back coefficients per frequency.
+        h = pytta.power_to_harmonic_series([1, 1, 1, 1], phase=0.3)
+        np.testing.assert_allclose(np.angle(h), [0, 0.3, 0.6, 0.9])
+        a = np.random.default_rng(2).standard_normal((5, 100)) + 1j
+        back = pytta.harmonic_to_power_series(pytta.power_to_harmonic_series(a, 0.3, 1.1), 0.3, 1.1)
+        print(f'\n  inverse: {np.abs(back - a).max():.1e}', end='')
+        self.assertLess(np.abs(back - a).max(), 1e-12)
+
+
+if __name__ == '__main__':
+    unittest.main()