Add pytta.fourpole with two-port matrices of lumped elements

docs/pyttadoc/fourpole.rst

@@ -0,0 +1,7 @@
+.. currentmodule:: pytta
+
+Four-pole matrices
+==================
+
+.. automodule:: pytta.fourpole
+   :members:

docs/pyttadoc/index.rst

@@ -64,4 +64,5 @@ Modules:
    functions
    generate
    rooms
+   fourpole
    properties

pytta/__init__.py

@@ -72,6 +72,9 @@
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.fourpole:
+        four-pole matrices of lumped elements;
+
 
 You can find usage examples in the toolbox's examples folder.
 
@@ -108,6 +111,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import fourpole
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +134,7 @@
            'generate',
            'utils',
            'iso3741',
+           'fourpole',
            'rooms',
 
            # Functions

pytta/fourpole.py

@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+"""
+Four-pole (two-port) matrices of lumped mechanical and acoustic elements,
+as itaFourpole and ita_make_fourpole of the ITA-Toolbox.
+
+The transmission matrix A of an element relates the force (or pressure) and
+the velocity (or volume velocity) at its input to those at its output,
+
+    [F1, v1] = A @ [F2, v2],
+
+with v2 leaving the element at the output. The matrices are arrays of shape
+(nFreqs, 2, 2), and a chain of elements is the product of their matrices.
+
+>>> f = np.linspace(1, 1000, 1000)
+>>> A = pytta.fourpole.chain(pytta.fourpole.element('mass', 0.1, f),
+...                          pytta.fourpole.element('spring', 1e-4, f))
+>>> Zin = pytta.fourpole.input_impedance(A)   # output blocked
+
+"""
+
+import numpy as np
+
+__all__ = ['element', 'chain', 'input_impedance', 'a_to_z', 'z_to_a']
+
+
+def element(kind, value, freqs):
+    """
+    Transmission matrix of one element at the frequencies freqs (Hz):
+
+        * 'mass': a mass of value kg, moving with both ports, A12 = j omega m;
+        * 'spring': a spring of compliance value m/N, A21 = j omega n;
+        * 'damper': a dashpot of resistance value N s/m, A21 = 1 / r;
+        * 'series': an impedance value (scalar or one per frequency) in
+          series, A12 = Z;
+        * 'shunt': an impedance value to ground, A21 = 1 / Z;
+        * 'transformer': A = diag(value, 1 / value);
+        * 'gyrator': A = [[0, value], [1 / value, 0]].
+    """
+    jw = 2j * np.pi * np.asarray(freqs, dtype=float)
+    value = np.broadcast_to(np.asarray(value, dtype=complex), jw.shape)
+    one, zero = np.ones_like(value), np.zeros_like(value)
+    a = {'mass': (one, jw * value, zero, one),
+         'spring': (one, zero, jw * value, one),
+         'damper': (one, zero, 1 / value, one),
+         'series': (one, value, zero, one),
+         'shunt': (one, zero, 1 / value, one),
+         'transformer': (value, zero, zero, 1 / value),
+         'gyrator': (zero, value, 1 / value, zero)}
+    if kind not in a:
+        raise ValueError(f'kind must be one of {", ".join(a)}.')
+    return np.stack(a[kind], axis=-1).reshape(jw.shape + (2, 2))
+
+
+def chain(*matrices):
+    """Transmission matrix of elements connected one after the other."""
+    result = matrices[0]
+    for m in matrices[1:]:
+        result = result @ m
+    return result
+
+
+def input_impedance(A, load=None):
+    """
+    F1 / v1 of the chain A with the impedance load at its output; None is a
+    blocked output (v2 = 0) and 0 a free one (F2 = 0).
+    """
+    if load is None:
+        return A[..., 0, 0] / A[..., 1, 0]
+    return (A[..., 0, 0] * load + A[..., 0, 1]) / (A[..., 1, 0] * load + A[..., 1, 1])
+
+
+def a_to_z(A):
+    """
+    Impedance matrix, [F1, F2] = Z @ [v1, -v2], of the transmission matrix
+    A, that is with both velocities entering the element, as
+    itaFourpole.A2Z; Z12 = Z21 for a reciprocal element (det A = 1).
+    """
+    a11, a12, a21, a22 = A[..., 0, 0], A[..., 0, 1], A[..., 1, 0], A[..., 1, 1]
+    det = a11 * a22 - a12 * a21
+    return np.stack([a11, det, np.ones_like(a11), a22], axis=-1).reshape(A.shape) / a21[..., None, None]
+
+
+def z_to_a(Z):
+    """Transmission matrix of the impedance matrix Z, the inverse of a_to_z."""
+    z11, z12, z21, z22 = Z[..., 0, 0], Z[..., 0, 1], Z[..., 1, 0], Z[..., 1, 1]
+    det = z11 * z22 - z12 * z21
+    return np.stack([z11, det, np.ones_like(z11), z22], axis=-1).reshape(Z.shape) / z21[..., None, None]

tests/test_fourpole.py

@@ -0,0 +1,47 @@
+"""
+Tests for pytta.fourpole, which follows itaFourpole and ita_make_fourpole of
+the ITA-Toolbox. Each test pins one rule against an analytic result.
+"""
+import unittest
+
+import numpy as np
+
+from pytta import fourpole
+
+F = np.linspace(10, 1000, 200)
+
+
+class TestFourpole(unittest.TestCase):
+
+    def test_mass_on_a_spring(self):
+        # A mass on a grounded spring has Z = j omega m + 1 / (j omega n),
+        # and a free output leaves only the mass.
+        m, n = 0.1, 1e-4
+        A = fourpole.chain(fourpole.element('mass', m, F), fourpole.element('spring', n, F))
+        jw = 2j * np.pi * F
+        e = np.abs(fourpole.input_impedance(A) / (jw * m + 1 / (jw * n)) - 1).max()
+        e_free = np.abs(fourpole.input_impedance(A, 0) / (jw * m) - 1).max()
+        print(f'\n  mass on spring: blocked {e:.1e}, free {e_free:.1e}', end='')
+        self.assertLess(max(e, e_free), 1e-12)
+
+    def test_impedance_matrix_is_reciprocal(self):
+        # A chain of passive elements has det A = 1, so Z12 = Z21, and the
+        # conversion goes back.
+        A = fourpole.chain(fourpole.element('mass', 0.2, F), fourpole.element('damper', 50, F),
+                           fourpole.element('spring', 3e-4, F), fourpole.element('shunt', 20 + 5j, F))
+        Z = fourpole.a_to_z(A)
+        e = np.abs(Z[:, 0, 1] - Z[:, 1, 0]).max() / np.abs(Z[:, 1, 0]).max()
+        back = np.abs(fourpole.z_to_a(Z) - A).max() / np.abs(A).max()
+        print(f'\n  reciprocity {e:.1e}, A -> Z -> A {back:.1e}', end='')
+        self.assertLess(max(e, back), 1e-12)
+
+    def test_gyrator_twice_and_transformer(self):
+        G = fourpole.element('gyrator', 7.5, F)
+        T = fourpole.element('transformer', 3, F)
+        np.testing.assert_allclose(fourpole.chain(G, G), np.broadcast_to(np.eye(2), G.shape), atol=1e-15)
+        # A load Z behind a transformer of ratio t is seen as t^2 Z.
+        np.testing.assert_allclose(fourpole.input_impedance(T, 4 + 1j), 9 * (4 + 1j))
+
+
+if __name__ == '__main__':
+    unittest.main()