Repository: pyfar/pyfar
Pull request #636 (2024-06-25 08:29:57)

Feature/transmission matrix
## Feature: New class representing transmission matrices

A transmission matrix (T-matrix) is used to model a two-port or quadrupole, which can represent electrical circuits, mechanical systems (e.g. loudspeakers), acoustic transmission lines, and so on. Systems can be cascaded by multiplying matrices. Functions are provided to calculate the input and output impedance as well as transfer functions.

### Extension towards acoustic transmission lines and other applications
My personal goal is to create a class / module that allows to handle and cascade different type of transmission line segments, which are internally represented by T-matrices. Generally, other applications are possible (e.g. loudspeaker modeling), but would be beyone the scope of my implementation.
- Should such a representation be part of pyfar or a new package? -> not pyfar
- Could be a new pyfar module? -> for now personal module

### Numpy 1.24 (Python 3.8) and logical indicing
Numpy >= 1.25 does not support Python 3.8 anymore (see [here](https://endoflife.date/numpy)). I had some issues with numpy 1.24 (latest Version for Python 3.8) regarding logical indicing. I found a work around but generally logical indicing seems to a little more tricky in older versions.
For higher Python Versions, numpy 2.x was used and my original code worked. I did not check higher numpy 1.XX versions.

### Notes
- A change of the internal `_arithmetic()` method of the `_Audio` class was required



Patch:
diff --git a/docs/api_reference.rst b/docs/api_reference.rst
index cee206fe..b354b1ed 100644
--- a/docs/api_reference.rst
+++ b/docs/api_reference.rst
@@ -21,6 +21,7 @@ Classes
    classes/pyfar.filter
    classes/pyfar.coordinates
    classes/pyfar.orientations
+   classes/pyfar.transmission_matrix
 
 Modules
 -------
diff --git a/docs/classes/pyfar.transmission_matrix.rst b/docs/classes/pyfar.transmission_matrix.rst
new file mode 100644
index 00000000..0b0d8802
--- /dev/null
+++ b/docs/classes/pyfar.transmission_matrix.rst
@@ -0,0 +1,7 @@
+TransmissionMatrix
+------------------
+
+.. automodule:: pyfar.classes.transmission_matrix
+   :members:
+   :undoc-members:
+   :inherited-members: pyfar.FrequencyData
diff --git a/pyfar/__init__.py b/pyfar/__init__.py
index 081c7bb9..935cc01a 100644
--- a/pyfar/__init__.py
+++ b/pyfar/__init__.py
@@ -14,6 +14,7 @@ from .classes.coordinates import Coordinates
 from .classes.coordinates import (deg2rad, rad2deg)
 from .classes.orientations import Orientations
 from .classes.filter import FilterFIR, FilterIIR, FilterSOS
+from .classes.transmission_matrix import TransmissionMatrix
 
 from . import plot
 from . import samplings
@@ -40,6 +41,7 @@ __all__ = [
     'FilterFIR',
     'FilterIIR',
     'FilterSOS',
+    'TransmissionMatrix',
     'plot',
     'samplings',
     'io',
diff --git a/pyfar/classes/__init__.py b/pyfar/classes/__init__.py
index f9a9e3a2..1080dfa3 100644
--- a/pyfar/classes/__init__.py
+++ b/pyfar/classes/__init__.py
@@ -2,6 +2,7 @@ from . import audio
 from . import coordinates
 from . import orientations
 from . import filter
+from . import transmission_matrix
 from . import warnings
 
 __all__ = [
@@ -9,5 +10,6 @@ __all__ = [
     'coordinates',
     'orientations',
     'filter',
+    'transmission_matrix',
     'warnings',
 ]
diff --git a/pyfar/classes/audio.py b/pyfar/classes/audio.py
index 7342139c..e79e0a08 100644
--- a/pyfar/classes/audio.py
+++ b/pyfar/classes/audio.py
@@ -1370,6 +1370,8 @@ def matrix_multiplication(
 
 def _arithmetic(data: tuple, domain: str, operation: Callable, **kwargs):
     """Apply arithmetic operations."""
+    #NOTE: The import is done here to avoid a circular import
+    from pyfar.classes.transmission_matrix import TransmissionMatrix
 
     # check input and obtain meta data of new signal
     division = True if operation == _divide else False
@@ -1400,6 +1402,8 @@ def _arithmetic(data: tuple, domain: str, operation: Callable, **kwargs):
         result = TimeData(result, times)
     elif audio_type == FrequencyData:
         result = FrequencyData(result, frequencies)
+    elif audio_type == TransmissionMatrix:
+        result = TransmissionMatrix(result, frequencies)
 
     return result
 
@@ -1444,6 +1448,8 @@ def _assert_match_for_arithmetic(data: tuple, domain: str, division: bool,
         Otherwise empty tuple.
 
     """
+    #NOTE: The import is done here to avoid a circular import
+    from pyfar.classes.transmission_matrix import TransmissionMatrix
 
     # we need at least two signals
     if not isinstance(data, tuple):
@@ -1467,7 +1473,7 @@ def _assert_match_for_arithmetic(data: tuple, domain: str, division: bool,
     # check input types and meta data
     n_audio_objects = 0
     for d in data:
-        if isinstance(d, (Signal, TimeData, FrequencyData)):
+        if isinstance(d,(Signal, TimeData, FrequencyData, TransmissionMatrix)):
             n_audio_objects += 1
             # store meta data upon first appearance
             if n_audio_objects == 1:
@@ -1479,7 +1485,7 @@ def _assert_match_for_arithmetic(data: tuple, domain: str, division: bool,
                     if domain != "time":
                         raise ValueError("The domain must be 'time'.")
                     times = d.times
-                elif isinstance(d, FrequencyData):
+                elif isinstance(d, (FrequencyData, TransmissionMatrix)):
                     if domain != "freq":
                         raise ValueError("The domain must be 'freq'.")
                     frequencies = d.frequencies
@@ -1502,7 +1508,7 @@ def _assert_match_for_arithmetic(data: tuple, domain: str, division: bool,
                     if not np.allclose(times, d.times, atol=1e-15):
                         raise ValueError(
                             "The times does not match.")
-                elif isinstance(d, FrequencyData):
+                elif isinstance(d, (FrequencyData, TransmissionMatrix)):
                     if not np.allclose(
                             frequencies, d.frequencies, atol=1e-15):
                         raise ValueError(
diff --git a/pyfar/classes/transmission_matrix.py b/pyfar/classes/transmission_matrix.py
new file mode 100644
index 00000000..693e4612
--- /dev/null
+++ b/pyfar/classes/transmission_matrix.py
@@ -0,0 +1,796 @@
+"""
+The following introduces the :py:class:`~TransmissionMatrix` class available
+in pyfar.
+
+Transmission matrices (short T-matrices) are a convenient representation of
+`Two-ports` or `Quadrupoles`. These can represent systems from various
+fields, e.g. electrical circuits, mechanical vibration, (acoustic) transmission
+lines.
+
+System properties like input impedance or transfer functions can directly be
+derived from a T-matrix. Furthermore, systems can be cascaded my multiplying
+consecutive T-matrices simply using the ``@`` operator:
+
+>>> import numpy as np
+>>> import pyfar as pf
+>>> frequencies = (100,200,300)
+>>>
+>>> # T-matrix with arbitrary data
+>>> A = np.ones(len(frequencies))
+>>> (B,C,D) = (A+1,A+2,A+3)
+>>> tmat = pf.TransmissionMatrix.from_abcd(A,B,C,D,frequencies)
+>>>
+>>> # T-matrix of a bypass system
+>>> tmat_bypass = pf.TransmissionMatrix.create_identity(frequencies)
+>>>
+>>> # Cascade
+>>> tmat_out = tmat @ tmat_bypass
+>>> tmat_out.freq == tmat.freq
+
+"""
+from __future__ import annotations # required for Python <= 3.9
+import numpy as np
+import numpy.testing as npt
+from pyfar.classes.audio import FrequencyData
+
+
+class TransmissionMatrix(FrequencyData):
+    r"""Class representing a transmission matrix.
+
+    This implementation is based on a paper by Lampton [#]_ and uses the ABCD-
+    representation. A single T-matrix is a (2x2)-matrix of the form:
+
+    .. math::
+        T = \begin{bmatrix}
+                A & B \\
+                C & D
+            \end{bmatrix}
+
+    This class is derived from :py:class:`~pyfar.classes.audio.FrequencyData`,
+    representing a frequency-dependent matrix of a multi-dimensional form
+    In the easiest case, each matrix entry (A,B,C,D) is a single-channel
+    FrequencyData object (i.e. a vector depending on frequency). However,
+    additional dimensions can be used to represent additional variables (e.g.
+    multiple layouts of an electrical circuit).
+
+    Notes
+    -----
+    This class is derived from
+    :py:class:`~pyfar.classes.audio.FrequencyData` but has a special
+    constraint: The underlying frequency domain data must match a shape like
+    (..., 2, 2, N). The last axis refers to the frequency, and the two axes
+    before to the ABCD-Matrix (which is a 2x2 matrix), so the resulting cshape
+    is (...,2,2).
+
+    Another point is the numerical handling when deriving input/output
+    impedance or transfer functions (see :py:func:`~input_impedance`,
+    :py:func:`~output_impedance`, :py:func:`~transfer_function`):
+    In the respective equations, the load impedance Zl is usually multiplied.
+    However, Zl is applied in the form of 1/Zl (as admittance) for Zl = inf
+    to avoid numerical problems.
+    There are additional cases where certain entries of the
+    T-matrix being zero might still lead to the main denominator becoming zero.
+    For these cases, the denominator is set to eps.
+
+    References
+    ----------
+    .. [#] M. Lampton. “Transmission Matrices in Electroacoustics". Acustica.
+           Vol. 39, pp. 239-251. 1978.
+
+    """
+
+    def __init__(self, data, frequencies, comment = ""):
+        """Initialize TransmissionMatrix with data, and frequencies.
+
+        This should not be used directly. Instead use :py:func:`~from_abcd`.
+        or :py:func:`~from_tmatrix`.
+
+        """
+        shape = np.shape(data)
+        n_dim = len(shape)
+        if n_dim < 3 or shape[-3] != 2 or shape[-2] != 2:
+            raise ValueError("'data' must have a shape like "
+                             "(..., 2, 2, n_bins), e.g. (2, 2, 100).")
+
+        super().__init__(data, frequencies, comment)
+
+    @classmethod
+    def from_tmatrix(cls, data, frequencies, comment = ""):
+        """Create TransmissionMatrix using data in T-matrix shape
+        and frequencies.
+
+        For using individual objects for A, B, C, D matrix entries,
+        see :py:func:`~from_abcd`.
+
+        Parameters
+        ----------
+        data : array_like, double
+            Raw data in the frequency domain.
+            In contrast to the :py:class:`~pyfar.classes.audio.FrequencyData`
+            class, data must have a shape of the form (..., 2, 2, N), e.g.
+            (2, 2, 1024) or (3, 2, 2, 1024). In those examples 1024 refers to
+            the number of frequency bins and the two dimensions before that to
+            the (2x2) ABCD matrices. Supported data types are ``int``,
+            ``float`` or ``complex``, where ``int`` is converted to ``float``.
+        frequencies : array_like, double
+            Frequencies of the data in Hz. The number of frequencies must match
+            the size of the last dimension of data.
+        comment : str, optional
+            A comment related to the data. The default is ``""``, which
+            initializes an empty string.
+
+        Examples
+        --------
+        >>> import pyfar as pf
+        >>> frequencies = (100,200,300)
+        >>> data = np.ones( (2, 2, len(frequencies)) )
+        >>> tmat = pf.TransmissionMatrix.from_tmatrix(data, frequencies)
+
+        >>> data = np.ones( (3, 2, 2, len(frequencies)) )
+        >>> tmat = pf.TransmissionMatrix.from_tmatrix(data, frequencies)
+
+        """
+        return cls(data, frequencies, comment)
+
+    @classmethod
+    def from_abcd(cls, A, B, C, D, frequencies = None):
+        """Create a TransmissionMatrix object from A-, B-, C-, D-data, and
+        frequencies.
+
+        Parameters
+        ----------
+        A : FrequencyData, array_like, double
+            Raw data for the matrix entries A. The data need to match the
+            B, C and D entry or be broadcastable into one ``shape``.
+        B : FrequencyData, array_like, double
+            See A.
+        C : FrequencyData, array_like, double
+            See A.
+        D : FrequencyData, array_like, double
+            See A.
+        frequencies : array, double, None
+            Frequencies of the data in Hz. This is optional if using the
+            FrequencyData type for A, B, C, D.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> import pyfar as pf
+        >>> frequencies = (100,200,300)
+
+        >>> # From np.array
+        >>> A = np.ones(len(frequencies))
+        >>> (B,C,D) = (A+1, A+2, A+3)
+        >>> tmat = pf.TransmissionMatrix.from_abcd(A,B,C,D, frequencies)
+        >>> tmat
+
+        >>> # From FrequencyData objects
+        >>> A = pf.FrequencyData( A, frequencies )
+        >>> (B,C,D) = (A+1, A+2, A+3)
+        >>> tmat = pf.TransmissionMatrix.from_abcd(A,B,C,D)
+        >>> tmat
+
+        >>> # Data with higher abcd dimension
+        >>> A = np.ones( (3, 4, len(frequencies)) )
+        >>> (B,C,D) = (A+1, A+2, A+3)
+        >>> tmat = pf.TransmissionMatrix.from_abcd(A,B,C,D, frequencies)
+        >>> tmat
+
+        """
+        num_freqdata = 0
+        for obj in (A,B,C,D):
+            if isinstance(obj, FrequencyData):
+                num_freqdata = num_freqdata + 1
+
+        if num_freqdata == 4:
+            frequencies = A.frequencies
+            for obj in (B,C,D): #Frequency bins must match
+                npt.assert_allclose(frequencies, obj.frequencies, atol=1e-15)
+            (A,B,C,D) = (A.freq, B.freq, C.freq, D.freq)
+        elif num_freqdata != 0:
+            raise ValueError(
+                        "If using FrequencyData objects, all matrix entries "
+                        "A, B, C, D, must be FrequencyData objects.")
+
+        if frequencies is None:
+            raise ValueError("'frequencies' must be specified if not using "
+                             "'FrequencyData' objects as input.")
+        # broadcast shapes
+        shape = np.broadcast_shapes(
+            np.array(A).shape, np.array(B).shape,
+            np.array(C).shape, np.array(D).shape)
+        A = np.broadcast_to(A, shape)
+        B = np.broadcast_to(B, shape)
+        C = np.broadcast_to(C, shape)
+        D = np.broadcast_to(D, shape)
+
+        data = np.array([[A, B], [C, D]])
+        # Switch dimension order so that T matrices refer to
+        # third and second last dimension (axes -3 and -2)
+        order = np.array(range(data.ndim-1))
+        order = np.roll(order, -2)
+        order = np.append(order, data.ndim-1)
+        data = np.transpose(data, order)
+
+        return cls(data, frequencies)
+
+    @property
+    def abcd_caxes(self):
+        """The indices of the channel axes referring to the transmission
+        matrix, namely (-2, -1)"""
+        return (-2, -1)
+
+    @property
+    def abcd_cshape(self):
+        """The channel shape of the transmission matrix entries (A, B, C, D)
+
+        This is the same as 'cshape' without the last two elements.
+        As an exception, a matrix with cshape (2,2) will return
+        (1,) as abcd_cshape.
+
+        """
+        abcd_cshape = self.cshape[:-2]
+        if abcd_cshape == ():
+            abcd_cshape = (1,)
+        return abcd_cshape
+
+    @property
+    def A(self) -> FrequencyData:
+        """A entry of the transmission matrix."""
+        return FrequencyData(self.freq[..., 0, 0, :], self.frequencies,
+                             self.comment)
+
+    @property
+    def B(self) -> FrequencyData:
+        """B entry of the transmission matrix."""
+        return FrequencyData(self.freq[..., 0, 1, :], self.frequencies,
+                             self.comment)
+
+    @property
+    def C(self) -> FrequencyData:
+        """C entry of the transmission matrix."""
+        return FrequencyData(self.freq[..., 1, 0, :], self.frequencies,
+                             self.comment)
+
+    @property
+    def D(self) -> FrequencyData:
+        """D entry of the transmission matrix."""
+        return FrequencyData(self.freq[..., 1, 1, :], self.frequencies,
+                             self.comment)
+
+    def _check_for_inf(self, Zl: complex | FrequencyData):
+        """Check given load impedance for np.inf values
+
+        Returns
+        -------
+        idx_inf : logical array
+            An np.ndarray of logicals pointing to elements referring to
+            Zl = inf. The shape is broadcasted to self.A.freq so that it can be
+            applied to the A-,B-,C-,D-entries.
+        idx_inf_Zl : logical array
+            An np.ndarray of logicals pointing to elements referring to
+            Zl = inf. The shape refers to an indexable version of given Zl.
+        Zl_indexable : np.ndarray
+            An indexable version of Zl.
+
+        """
+        if isinstance(Zl, FrequencyData):
+            Zl = Zl.freq
+        elif np.isscalar(Zl):
+            Zl = np.atleast_1d(Zl)
+
+        target_shape = list(self.abcd_cshape) + [self.n_bins]
+        idx_inf_Zl = Zl == np.inf
+        if target_shape is not None:
+            idx_inf = np.broadcast_to(idx_inf_Zl, target_shape)
+
+        return idx_inf, idx_inf_Zl, Zl
+
+    def input_impedance(self, Zl: complex | FrequencyData) -> FrequencyData:
+        """Calculates the input impedance given the load impedance Zl at the
+        output.
+
+        Two-port representation::
+
+                    o---xxxxxxxxx---o
+                        x       x   |
+            Zin-->      x       x   Zl
+                        x       x   |
+                    o---xxxxxxxxx---o
+
+        See Equation (2-6) in Reference [1]_:
+        :math:`Z_\\mathrm{in} = \\frac{AZ_L + B}{CZ_L + D}`
+
+        Parameters
+        ----------
+        Zl : scalar | FrequencyData
+            The load impedance data as scalar or FrequencyData. In latter case,
+            the shape must match the entries of the T-matrix, i.e.
+            shape(tmat.A.freq) == shape(Zl.freq), or must be broadcastable.
+
+        Returns
+        -------
+        Zout : FrequencyData
+            A FrequencyData object with the resulting output impedance. The
+            cshape is identical to the entries of the T-matrix, i.e.
+            tmat.A.cshape == Zout.cshape.
+
+        Example
+        -------
+
+        >>> import numpy as np
+        >>> import pyfar as pf
+        >>>
+        >>> # Frequency-dependent load impedance
+        >>> frequencies = (100,200,300)
+        >>> load_impedance = pf.FrequencyData((0, 1, np.inf), frequencies)
+        >>>
+        >>> # T-Matrix for frequency-independent series impedance R = 1 Ohm
+        >>> R = pf.FrequencyData((1,1,1), frequencies)
+        >>> tmat = pf.TransmissionMatrix.create_series_impedance(R)
+        >>>
+        >>> # Expected result: (1+0, 1+1, 1+inf) = (1, 2, inf)
+        >>> # Note, that due to numerical limitations infinite load will
+        >>> # result in Zin > 1e15.
+        >>> Zin = tmat.input_impedance(load_impedance)
+        >>> Zin.freq
+
+        """
+        nominator = (self.A * Zl + self.B)
+        denominator = (self.C * Zl + self.D)
+
+        # Admittance form for Zl = inf
+        idx_inf, __, __ = self._check_for_inf(Zl)
+        nominator.freq[idx_inf] = (self.A + self.B / Zl).freq[idx_inf]
+        denominator.freq[idx_inf] = (self.C + self.D / Zl).freq[idx_inf]
+
+        # Avoid cases where denominator is zero, examples
+        # Zl = inf & C = 0; Zl = 0 & D = 0
+        denominator.freq[denominator.freq == 0] = np.finfo(float).eps
+        return nominator / denominator
+
+    def output_impedance(self, Zl: complex | FrequencyData) -> FrequencyData:
+        """Calculates the output impedance given the load impedance Zl at the
+        input.
+
+        Two-port representation::
+
+                    o---xxxxxxxxx---o
+                    |   x       x
+                    Zl  x       x      <--Zout
+                    |   x       x
+                    o---xxxxxxxxx---o
+
+        See Equation (2-6) in Reference [1]_:
+        :math:`Z_\\mathrm{out} = \\frac{DZ_L + B}{CZ_L + A}`
+
+        For a code example, see :py:func:`~input_impedance` and exchange
+        respective method call with `output_impedance`.
+
+        Parameters
+        ----------
+        Zl : scalar | FrequencyData
+            The load impedance data as scalar or FrequencyData. In latter case,
+            the shape must match the entries of the T-matrix, i.e.
+            shape(tmat.A.freq) == shape(Zl.freq), or must be broadcastable.
+
+        Returns
+        -------
+        Zout : FrequencyData
+            A FrequencyData object with the resulting output impedance. The
+            cshape is identical to the entries of the T-matrix, i.e.
+            tmat.A.cshape == Zout.cshape.
+
+        """
+        nominator = (self.D * Zl + self.B)
+        denominator = (self.C * Zl + self.A)
+
+        # Admittance form for Zl = inf
+        idx_inf, __, __ = self._check_for_inf(Zl)
+        nominator.freq[idx_inf] = (self.D + self.B / Zl).freq[idx_inf]
+        denominator.freq[idx_inf] = (self.C + self.A / Zl).freq[idx_inf]
+
+        # Avoid cases where denominator is zero, examples
+        # Zl = 0 & A = 0; Zl = inf & C = 0
+        denominator.freq[denominator.freq == 0] = np.finfo(float).eps
+        return nominator / denominator
+
+    def transfer_function(self, quantity_indices,
+                          Zl: complex | FrequencyData) -> FrequencyData:
+        """Returns the transfer function (output/input) for specified
+        quantities and a given load impedance.
+
+        The transfer function is the relation between an output and input
+        quantity of modelled two-port and depends on the load impedance at the
+        output. Since there are two quantities at input and output
+        respectively, four transfer functions exist in total. The first usually
+        refers to the "voltage-like" quantity (:math:`Q_1`) whereas the second
+        refers to the "current-like" quantity (:math:`Q_2`).
+
+        The transfer functions can be derived from Equation (2-1)
+        in Reference [1]_:
+
+        .. math::
+            Q_{1,\\mathrm{in}} = AQ_{1,\\mathrm{out}} + BQ_{2,\\mathrm{in}}
+
+            Q_{2,\\mathrm{in}} = CQ_{1,\\mathrm{out}} + DQ_{2,\\mathrm{in}}
+
+        The four transfer functions are defined as:
+
+        * :math:`Q_{1,\\mathrm{out}} / Q_{1,\\mathrm{in}}` using
+          :math:`Q_{2,\\mathrm{out}} = Q_{1,\\mathrm{out}}/Z_L`
+        * :math:`Q_{2,\\mathrm{out}} / Q_{1,\\mathrm{in}} =
+          Q_{1,\\mathrm{out}} / Q_{1,\\mathrm{in}} \\cdot \\frac{1}{Z_L}`
+        * :math:`Q_{2,\\mathrm{out}} / Q_{2,\\mathrm{in}}` using
+          :math:`Q_{1,\\mathrm{out}} = Q_{2,\\mathrm{out}}\\cdot Z_L`
+        * :math:`Q_{1,\\mathrm{out}} / Q_{2,\\mathrm{in}} =
+          Q_{2,\\mathrm{out}} / Q_{2,\\mathrm{in}} \\cdot Z_L`
+
+        Parameters
+        ----------
+        quantity_indices : array_like (int, int)
+            Array-like object with two integer elements referring to the
+            indices of the utilized quantity at the output (first integer)
+            and input (second integer). For example, (1,0) refers to the
+            transfer function :math:`Q_{2,\\mathrm{out}} / Q_{1,\\mathrm{in}}`.
+        Zl : scalar | FrequencyData
+            The load impedance data as scalar or FrequencyData. In latter case,
+            the shape must match the entries of the T-matrix, i.e.
+            shape(tmat.A.freq) == shape(Zl.freq), or must be broadcastable.
+
+        Returns
+        -------
+        transfer_function : FrequencyData
+            A FrequencyData object with the resulting transfer function. The
+            cshape is identical to the entries of the T-matrix, i.e.
+            tmat.A.cshape == transfer_function.cshape.
+
+        """
+        is_scalar = np.isscalar(quantity_indices)
+        quantity_indices = np.array(quantity_indices)
+        is_numeric = np.issubdtype(quantity_indices.dtype, np.number)
+        if is_scalar or not is_numeric or not len(quantity_indices) == 2:
+            raise ValueError("'quantity_indices' must be an array-like type "
+                             "with two numeric elements.")
+        if not all(np.logical_or(quantity_indices==0, quantity_indices==1)):
+            raise ValueError("'quantity_indices' must contain two integers "
+                             "between 0 and 1.")
+
+        if quantity_indices[0] == 0 and quantity_indices[1] == 0:
+            return self._transfer_function_q1q1(Zl)
+        if quantity_indices[0] == 1 and quantity_indices[1] == 1:
+            return self._transfer_function_q2q2(Zl)
+        if quantity_indices[0] == 1 and quantity_indices[1] == 0:
+            return self._transfer_function_q2q1(Zl)
+        if quantity_indices[0] == 0 and quantity_indices[1] == 1:
+            return self._transfer_function_q1q2(Zl)
+
+    def _transfer_function_q1q1(self,
+                                Zl: complex | FrequencyData) -> FrequencyData:
+        """Returns the first quantity's transfer function (Q1_out/Q1_in)."""
+        idx_inf, __, __ = self._check_for_inf(Zl)
+        denominator = (self.A * Zl + self.B)
+        nominator = Zl*FrequencyData(np.ones_like(denominator.freq),
+                                     self.frequencies)
+
+        # Admittance form for Zl = inf
+        nominator.freq[idx_inf] = 1
+        denominator.freq[idx_inf] = (self.A + self.B / Zl).freq[idx_inf]
+
+        # Avoid cases where denominator is zero, examples
+        # Zl = 0 & B = 0; Zl = inf & A = 0
+        denominator.freq[denominator.freq == 0] = np.finfo(float).eps
+        return nominator / denominator
+
+    def _transfer_function_q2q1(self,
+                                Zl: complex | FrequencyData) -> FrequencyData:
+        """Returns the transfer function Q2_out / Q1_in."""
+        denominator = (self.A * Zl + self.B)
+        nominator = FrequencyData(np.ones_like(denominator.freq),
+                                  self.frequencies)
+
+        # In cases where the denominator is zero, e.g. Zl = 0 & B = 0,
+        # are related to short-circuated outputs (undefined current) => NaN
+        denominator.freq[denominator.freq == 0] = np.nan
+        return nominator / denominator
+
+    def _transfer_function_q2q2(self,
+                                Zl: complex | FrequencyData) -> FrequencyData:
+        """Returns the second quantity's transfer function (Q2_out/Q2_in)."""
+        denominator = (self.C * Zl + self.D)
+        nominator = FrequencyData(np.ones_like(denominator.freq),
+                                  self.frequencies)
+
+        # Admittance form for Zl = inf
+        idx_inf, idx_inf_Zl, Zl_indexable = self._check_for_inf(Zl)
+        nominator.freq[idx_inf] = 1 / Zl_indexable[idx_inf_Zl]
+        denominator.freq[idx_inf] = (self.C + self.D / Zl).freq[idx_inf]
+
+        # Avoid cases where denominator is zero, examples
+        # Zl = 0 & D = 0; Zl = inf & C = 0
+        denominator.freq[denominator.freq == 0] = np.finfo(float).eps
+        return nominator / denominator
+
+    def _transfer_function_q1q2(self,
+                                Zl: complex | FrequencyData) -> FrequencyData:
+        """Returns the transfer function Q1_out / Q2_in."""
+        denominator = (self.C * Zl + self.D)
+        nominator = Zl*FrequencyData(np.ones_like(denominator.freq),
+                                     self.frequencies)
+
+        # Admittance form for Zl = inf
+        idx_inf, __, __  = self._check_for_inf(Zl)
+        nominator.freq[idx_inf] = 1
+        denominator.freq[idx_inf] = (self.C + self.D / Zl).freq[idx_inf]
+
+        # In cases where the denominator is zero, e.g. Zl = inf & C = 0,
+        # are related to non-physical cases (e.g. undefined current/voltage)
+        # => NaN
+        denominator.freq[denominator.freq == 0] = np.nan
+        return nominator / denominator
+
+    @staticmethod
+    def create_identity(frequencies = None):
+        r"""Creates an object with identity matrix entries (bypass).
+
+        See Equation (2-7) in Table I of Reference [1]_:
+
+        .. math::
+            T = \begin{bmatrix}
+                1 & 0 \\
+                0 & 1
+                \end{bmatrix}
+
+        Parameters
+        ----------
+        frequencies : None | array_like, optional
+            The frequency sampling points in Hz. The default is `None` which
+            will result in this function to return an np.ndarray instead of an
+            TransmissionMatrix object.
+
+        Returns
+        -------
+        tmat : np.ndarray | TransmissionMatrix
+            If frequencies are specified, a TransmissionMatrix object that
+            contains one 2x2 identity matrix per bin is returned. Otherwise,
+            this returns a 2x2 np.ndarray.
+
+        """
+        if frequencies is None or len(frequencies) == 0:
+            return np.eye(2)
+
+        return TransmissionMatrix.from_abcd(np.ones_like(frequencies),
+                                            0, 0, 1, frequencies)
+
+    @staticmethod
+    def create_series_impedance(impedance: complex | FrequencyData,
+                                ) -> np.ndarray | TransmissionMatrix:
+        r"""Creates a transmission matrix representing a series impedance.
+
+        This means the impedance is connected in series with a potential load
+        impedance. See Equation (2-8) in Table I of Reference [1]_:
+
+        .. math::
+            T = \begin{bmatrix}
+                1 & Z \\
+                0 & 1
+                \end{bmatrix}
+
+        Parameters
+        ----------
+        impedance : scalar | FrequencyData
+            The impedance data of the series impedance.
+
+        Returns
+        -------
+        tmat : np.ndarray | TransmissionMatrix
+            A transmission matrix representing the series connection
+            and can be cascaded with TransmissionMatrix objects.
+            If a scalar was used as input a frequency-independent
+            matrix is returned, namely an np.ndarray of shape (2,2).
+
+        """
+        if np.isscalar(impedance) and not isinstance(impedance, str):
+            tmat = np.identity(2)
+            tmat[0,1] = impedance
+            return tmat
+
+        if not isinstance(impedance, FrequencyData):
+            raise ValueError("'impedance' must be a "
+                             "numerical scalar or FrequencyData object.")
+        return TransmissionMatrix.from_abcd(
+            1, impedance.freq, 0, 1, impedance.frequencies)
+
+    @staticmethod
+    def create_shunt_admittance(admittance: complex | FrequencyData,
+                                ) -> np.ndarray | TransmissionMatrix:
+        r"""Creates a transmission matrix representing a shunt admittance
+        (parallel connection).
+
+        In this case, the impedance (= 1 / admittance) is connected in parallel
+        with a potential load impedance.
+        See Equation (2-9) in Table I of Reference [1]_:
+
+        .. math::
+            T = \begin{bmatrix}
+                1 & 0 \\
+                Y & 1
+                \end{bmatrix}
+
+        Parameters
+        ----------
+        admittance : scalar | FrequencyData
+            The admittance data of the element connected in parallel.
+
+        Returns
+        -------
+        tmat : np.ndarray | TransmissionMatrix
+            A transmission matrix representing a parallel connection
+            and can be cascaded with TransmissionMatrix objects.
+            If a scalar was used as input a frequency-independent
+            matrix is returned, namely an np.ndarray of shape (2,2).
+
+        """
+        if np.isscalar(admittance) and not isinstance(admittance, str):
+            tmat = np.identity(2)
+            tmat[1,0] = admittance
+            return tmat
+
+        if not isinstance(admittance, FrequencyData):
+            raise ValueError("'admittance' must be a "
+                             "numerical scalar or FrequencyData object.")
+        return TransmissionMatrix.from_abcd(
+            1, 0, admittance.freq, 1, admittance.frequencies)
+
+    @staticmethod
+    def create_transformer(transducer_constant :
+                complex | FrequencyData) -> np.ndarray | TransmissionMatrix:
+        r"""Creates a transmission matrix representing a transformer.
+
+        See Equation (2-12) in Table I of Reference [1]_:
+
+        .. math::
+            T = \begin{bmatrix}
+                N & 0 \\
+                0 & 1/N
+                \end{bmatrix}
+
+        Parameters
+        ----------
+        transducer_constant : scalar | FrequencyData
+            The transmission ratio with respect to voltage-like quantity,
+            i.e. :math:`N=U_\mathrm{out}/U_\mathrm{in}`. If a scalar is given,
+            i.e. a frequency-independent transformer matrix is requested, the
+            return value will be a 2x2 np.ndarray isntead.
+
+        Returns
+        -------
+        tmat : np.ndarray | TransmissionMatrix
+            A transmission matrix representing the transformer
+            and can be cascaded with TransmissionMatrix objects.
+            If a scalar was used as input a frequency-independent
+            matrix is returned, namely an np.ndarray of shape (2,2).
+
+        """
+        if np.isscalar(transducer_constant) and not isinstance(
+            transducer_constant, str):
+            tmat = np.identity(2)
+            tmat[0,0] = transducer_constant
+            tmat[1,1] = 1/transducer_constant
+            return tmat
+
+        if not isinstance(transducer_constant, FrequencyData):
+            raise ValueError("'transducer_constant' must be a "
+                             "numerical scalar or FrequencyData object.")
+        A = transducer_constant.freq
+        D = (1/transducer_constant).freq
+        frequencies = transducer_constant.frequencies
+        return TransmissionMatrix.from_abcd(A, 0, 0, D, frequencies)
+
+    @staticmethod
+    def create_gyrator(transducer_constant :
+                complex | FrequencyData) -> np.ndarray | TransmissionMatrix:
+        r"""Creates a transmission matrix representing a gyrator.
+
+        The T-matrix is defined by a transducer constant (:math:`M`),
+        see Equation (2-14) in Table I of Reference [1]_:
+
+        .. math::
+            T = \begin{bmatrix}
+                0 & M \\
+                1/M & 0
+                \end{bmatrix}
+
+        :math:`M` connects the first input and second output quantity (e.g.,
+        :math:`U_\mathrm{out} = I_\mathrm{in} \cdot M`). A respective system
+        with load :math:`Z_L` has the input impedance
+        :math:`Z_\mathrm{in} = M^2 / Z_L`.
+
+        Parameters
+        ----------
+        transducer_constant : scalar | FrequencyData
+            The transducer constant :math:`M`. If a scalar is given,
+            i.e. a frequency-independent transformer matrix is requested, the
+            return value will be a 2x2 np.ndarray instead.
+
+        Returns
+        -------
+        tmat : np.ndarray | TransmissionMatrix
+            A the transmission matrix representing the gyrator and can be
+            cascaded with TransmissionMatrix objects. If a scalar was used as
+            input a frequency-independent matrix is returned, namely an
+            np.ndarray of shape (2,2).
+
+        """
+        if np.isscalar(transducer_constant) and not isinstance(
+            transducer_constant, str):
+            tmat = np.zeros([2,2])
+            tmat[0,1] = transducer_constant
+            tmat[1,0] = 1/transducer_constant
+            return tmat
+
+        if not isinstance(transducer_constant, FrequencyData):
+            raise ValueError("'transducer_constant' must be a "
+                             "numerical scalar or FrequencyData object.")
+        B = transducer_constant.freq
+        C = (1/transducer_constant).freq
+        frequencies = transducer_constant.frequencies
+        return TransmissionMatrix.from_abcd(0, B, C, 0, frequencies)
+
+    def __repr__(self):
+        """String representation of TransmissionMatrix class."""
+        repr_string = (
+            f"TransmissionMatrix:\n"
+            f"{self.abcd_cshape} channels per matrix entry "
+            f"with {self.n_bins} frequencies\n"
+            f"Total number of {self.cshape} channels")
+
+        return repr_string
+
+    @classmethod
+    def _decode(cls, obj_dict):
+        """Decode object based on its respective `_encode` counterpart."""
+        obj = cls(
+            obj_dict['_data'],
+            obj_dict['_frequencies'],
+            obj_dict['_comment'])
+        obj.__dict__.update(obj_dict)
+        return obj
+
+    def is_indexable(self) -> bool:
+        """Returns true if ABCD-entries have more than one channel and are
+        therefore indexable.
+        """
+        return len(self.cshape) > 2
+
+    def __getitem__(self, key):
+        """Get copied slice of the TransmissionMatrix at key.
+
+            Note, that slicing ABCD or frequency dimensions is not possible.
+        """
+
+        if not self.is_indexable():
+            raise IndexError(
+                "Object is not indexable, since ABCD-entries " \
+                "only have a single channel")
+
+        # Add three empty slices at the end to always get all data contained
+        # in frequency and T-Matrix dimensions (last three dimensions).
+        if hasattr(key, '__iter__'):
+            key = (*key, slice(None), slice(None), slice(None))
+
+        # try indexing and raise verbose errors if it fails
+        try:
+            data = self._data[key]
+        except IndexError as Error:
+            if 'too many indices for array' in str(Error):
+                raise IndexError((
+                    f'Indexed dimensions must not exceed the ABCD '
+                    f'channel dimension (abcd_cdim), which is '
+                    f'{len(self.abcd_cshape)}'))
+            else:
+                raise Error
+
+        return TransmissionMatrix.from_tmatrix(
+            data, frequencies=self.frequencies,
+            comment=self.comment)
diff --git a/pyfar/io/_codec.py b/pyfar/io/_codec.py
index c78cfe55..59cacade 100644
--- a/pyfar/io/_codec.py
+++ b/pyfar/io/_codec.py
@@ -398,6 +398,7 @@ def _is_pyfar_type(obj):
         'SphericalVoronoi',
         'TimeData',
         'FrequencyData',
+        'TransmissionMatrix',
         'BuiltinsWrapper']
 
 


Test patch:
diff --git a/tests/test_imports.py b/tests/test_imports.py
index ab36e7bb..b814f371 100644
--- a/tests/test_imports.py
+++ b/tests/test_imports.py
@@ -18,6 +18,7 @@ def test_import_classes():
     from pyfar import Signal                 # noqa: F401
     from pyfar import TimeData               # noqa: F401
     from pyfar import FrequencyData          # noqa: F401
+    from pyfar import TransmissionMatrix     # noqa: F401
 
     from pyfar import Coordinates            # noqa: F401
     from pyfar import Orientations           # noqa: F401
diff --git a/tests/test_io.py b/tests/test_io.py
index 8be336a0..94da63c0 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -16,7 +16,7 @@ from pyfar import Signal
 from pyfar import Coordinates
 from pyfar.samplings import SphericalVoronoi
 import pyfar.classes.filter as fo
-from pyfar import FrequencyData, TimeData
+from pyfar import FrequencyData, TimeData, TransmissionMatrix
 
 
 @pytest.mark.parametrize('input_type', ('filename', 'path_object'))
@@ -261,6 +261,17 @@ def test_write_read_frequencydata(frequency_data, tmpdir):
     assert isinstance(actual, FrequencyData)
     assert actual == frequency_data
 
+def test_write_read_transmissionmatrix(tmpdir):
+    """ TransmissionMatrix
+    Make sure `read` understands the bits written by `write`
+    """
+    tmatrix = TransmissionMatrix.create_transformer(FrequencyData(42, 100))
+    filename = os.path.join(tmpdir, 'transmissionmatrix.far')
+    io.write(filename, transmissionmatrix=tmatrix)
+    actual = io.read(filename)['transmissionmatrix']
+    assert isinstance(actual, FrequencyData)
+    assert actual == tmatrix
+
 
 def test_write_read_sphericalvoronoi(sphericalvoronoi, tmpdir):
     """ SphericalVoronoi
diff --git a/tests/test_transmission_matrix.py b/tests/test_transmission_matrix.py
new file mode 100644
index 00000000..a9e23970
--- /dev/null
+++ b/tests/test_transmission_matrix.py
@@ -0,0 +1,300 @@
+import numpy as np
+import numpy.testing as npt
+import pytest
+import re
+import pyfar as pf
+from pyfar import TransmissionMatrix
+from pyfar import FrequencyData
+
+@pytest.fixture(scope="module")
+def frequencies():
+    return [100, 200, 300]
+@pytest.fixture(scope="module")
+def A_list():
+    """Test data for a matrix-entry (e.g. A) using a list type"""
+    return [1, 2, 3]
+@pytest.fixture(scope="module")
+def A_np(A_list):
+    """Test data for a matrix-entry (e.g. A) using an np.ndarray"""
+    return np.array(A_list)
+@pytest.fixture(scope="module")
+def A_FreqDat(A_np, frequencies):
+    """Test data for a matrix-entry (e.g. A) using a FrequencyData object"""
+    return FrequencyData(A_np, frequencies)
+
+def _expect_data_with_wrong_abcd_dims(data: np.ndarray, frequencies):
+    error_msg = re.escape("'data' must have a shape like "
+                          "(..., 2, 2, n_bins), e.g. (2, 2, 100).")
+    with pytest.raises(ValueError, match=error_msg):
+        TransmissionMatrix(data, frequencies)
+    with pytest.raises(ValueError, match=error_msg):
+        TransmissionMatrix.from_tmatrix(data, frequencies)
+    with pytest.raises(ValueError, match=error_msg):
+        TransmissionMatrix(np.ndarray.tolist(data), frequencies)
+    with pytest.raises(ValueError, match=error_msg):
+        TransmissionMatrix.from_tmatrix(np.ndarray.tolist(data), frequencies)
+
+def test_tmatrix_init(frequencies):
+    """Test constructor with data of valid and invalid shapes."""
+    num_bins = len(frequencies)
+    TransmissionMatrix(np.ones([2, 2, num_bins]), frequencies)
+    TransmissionMatrix(np.ones([4, 2, 2, num_bins]), frequencies)
+    TransmissionMatrix.from_tmatrix(np.ones([2, 2, num_bins]), frequencies)
+    TransmissionMatrix.from_tmatrix(np.ones([4, 2, 2, num_bins]), frequencies)
+
+    _expect_data_with_wrong_abcd_dims(
+        np.ones([2, num_bins]), frequencies)
+    _expect_data_with_wrong_abcd_dims(
+        np.ones([3, 2, num_bins]), frequencies)
+    _expect_data_with_wrong_abcd_dims(
+        np.ones([2, 5, num_bins]), frequencies)
+    _expect_data_with_wrong_abcd_dims(
+        np.ones([7, 4, 2, num_bins]), frequencies)
+    _expect_data_with_wrong_abcd_dims(
+        np.ones([7, 8, 4, 2, num_bins]), frequencies)
+
+def _expect_error_abcd_same_type(A, B, C, D):
+    with pytest.raises(
+        ValueError, match=
+                    "If using FrequencyData objects, all matrix entries "
+                    "A, B, C, D, must be FrequencyData objects."):
+        TransmissionMatrix.from_abcd(A, B, C, D, 1000)
+
+def test_tmatrix_from_abcd_input_types(frequencies, A_list, A_np, A_FreqDat):
+    """Test 'from_abcd' with valid and invalid data types."""
+    TransmissionMatrix.from_abcd(A_list, A_list,
+                                 A_list, A_list, frequencies)
+    TransmissionMatrix.from_abcd(A_np, A_np, A_list, A_np, frequencies)
+    TransmissionMatrix.from_abcd(A_np, A_np, A_np, A_np, frequencies)
+    TransmissionMatrix.from_abcd(A_FreqDat, A_FreqDat, A_FreqDat, A_FreqDat)
+
+    _expect_error_abcd_same_type(A_np, A_np, A_np, A_FreqDat)
+    _expect_error_abcd_same_type(A_np, A_np, A_FreqDat, A_np)
+    _expect_error_abcd_same_type(A_np, A_FreqDat, A_np, A_np)
+    _expect_error_abcd_same_type(A_FreqDat, A_np, A_np, A_np)
+
+def test_tmatrix_from_abcd_optional_frequencies(A_list, A_FreqDat):
+    """Test 'from_abcd' throws error if handing in arrays but no frequencies"""
+    TransmissionMatrix.from_abcd(A_FreqDat, A_FreqDat, A_FreqDat, A_FreqDat)
+    with pytest.raises(ValueError, match="'frequencies' must be specified if "
+                       "not using 'FrequencyData' objects as input"):
+        TransmissionMatrix.from_abcd(A_list, A_list, A_list, A_list)
+
+
+# -------------------------
+# TESTS FOR HIGHER DIM DATA
+# -------------------------
+@pytest.fixture(scope="module")
+def abcd_data_1x2():
+    """ABCD matrices with 2 frequency bins and one additional
+    dimension of size 3"""
+    frequencies = [100, 200]
+    A = FrequencyData([[1, 1]], frequencies)
+    B = FrequencyData([[2, 2]], frequencies)
+    C = FrequencyData([[3, 3]], frequencies)
+    D = FrequencyData([[4, 4]], frequencies)
+    tmat = TransmissionMatrix.from_abcd(A, B, C, D)
+    return tmat, A, B, C, D
+@pytest.fixture(scope="module")
+def abcd_data_3x2():
+    """ABCD matrices with 2 frequency bins and one additional
+    dimension of size 3"""
+    frequencies = [100, 200]
+    A = FrequencyData([[1, 1], [1, 1], [1, 1]], frequencies)
+    B = FrequencyData([[2, 2], [2, 2], [2, 2]], frequencies)
+    C = FrequencyData([[3, 3], [3, 3], [3, 3]], frequencies)
+    D = FrequencyData([[4, 4], [4, 4], [4, 4]], frequencies)
+    tmat = TransmissionMatrix.from_abcd(A, B, C, D)
+    return tmat, A, B, C, D
+@pytest.fixture(scope="module")
+def abcd_data_3x3x1():
+    """ABCD matrices with 1 frequency bin and two additional
+    dimensions of size 3"""
+    A = FrequencyData(
+        [[[1.1], [1.1], [1.1]], [[1.2], [1.2], [1.2]], [[1.3], [1.3], [1.3]]],
+        100)
+    B = A + 1
+    C = A + 2
+    D = A + 3
+    tmat = TransmissionMatrix.from_abcd(A, B, C, D)
+    return tmat, A, B, C, D
+
+def test_tmatrix_abcd_cshape(abcd_data_1x2, abcd_data_3x2, abcd_data_3x3x1):
+    """Test whether abcd_cshape matches cshape of A-property."""
+    tmat, A, __, __, __ = abcd_data_1x2
+    assert tmat.abcd_cshape == A.cshape
+    tmat, A, __, __, __ = abcd_data_3x2
+    assert tmat.abcd_cshape == A.cshape
+    tmat, A, __, __, __  = abcd_data_3x3x1
+    assert tmat.abcd_cshape == A.cshape
+
+def _compare_tmat_vs_abcd(tmat, A, B, C, D):
+    """Test whether ABCD entries of T-Matrix match given data sets"""
+    if isinstance(A, FrequencyData):
+        assert tmat.A == A
+        assert tmat.B == B
+        assert tmat.C == C
+        assert tmat.D == D
+    else:
+        assert np.all(tmat.A.freq == A)
+        assert np.all(tmat.B.freq == B)
+        assert np.all(tmat.C.freq == C)
+        assert np.all(tmat.D.freq == D)
+
+def test_tmatrix_abcd_entries(abcd_data_3x2, abcd_data_3x3x1):
+    """Test whether ABCD entries of T-Matrix match ABCD data used for
+     initialization"""
+    tmat, A, B, C, D = abcd_data_3x2
+    _compare_tmat_vs_abcd(tmat, A, B, C, D)
+
+    tmat, A, B, C, D = abcd_data_3x3x1
+    _compare_tmat_vs_abcd(tmat, A, B, C, D)
+
+
+# ------------------------
+# TESTS FOR CREATE METHODS
+# ------------------------
+def test_tmatrix_create_identity(frequencies):
+    """Test whether creation of identity matrix with frequencies"""
+    tmat_eye = TransmissionMatrix.create_identity(frequencies)
+    assert isinstance(tmat_eye, TransmissionMatrix)
+    assert tmat_eye.abcd_cshape == (1,)
+    _compare_tmat_vs_abcd(tmat_eye, 1, 0, 0, 1)
+
+@pytest.mark.parametrize("no_freqs", [None, ()])
+def test_tmatrix_create_identity_scalar_input(no_freqs):
+    """Test whether creation of identity matrix without frequencies"""
+    eye = TransmissionMatrix.create_identity(no_freqs)
+    assert isinstance(eye, np.ndarray)
+    npt.assert_allclose(eye, np.eye(2), atol=1e-15)
+
+@pytest.mark.parametrize("method_name",
+                         ["create_series_impedance", "create_shunt_admittance",
+                          "create_transformer", "create_gyrator"])
+def test_tmatrix_create_methods_wrong_input(method_name):
+    """Test create methods raise error on invalid input data."""
+    if method_name == "create_series_impedance":
+        input_name = "impedance"
+    elif method_name == "create_shunt_admittance":
+        input_name = "admittance"
+    else:
+        input_name = "transducer_constant"
+
+    func = getattr(TransmissionMatrix, method_name)
+    err_msg = "'" + input_name + "' "
+    "must be a numerical scalar or FrequencyData object."
+    with pytest.raises(ValueError, match=err_msg):
+        func("wrong_input")
+    with pytest.raises(ValueError, match=err_msg):
+        func('wrong_input')
+    with pytest.raises(ValueError, match=err_msg):
+        func((1,2,3))
+    with pytest.raises(ValueError, match=err_msg):
+        func([1,2,3])
+
+
+@pytest.mark.parametrize("abcd_cshape", [(1,), (4,5)])
+def test_tmatrix_create_series_impedance(A_FreqDat, abcd_cshape):
+    """Test `create_series_impedance` using FrequencyData input."""
+    Z = pf.utils.broadcast_cshape(A_FreqDat, abcd_cshape)
+    tmat = TransmissionMatrix.create_series_impedance(Z)
+    assert isinstance(tmat, TransmissionMatrix)
+    assert tmat.abcd_cshape == abcd_cshape
+    _compare_tmat_vs_abcd(tmat, 1, Z.freq, 0, 1)
+
+def test_tmatrix_create_series_impedance_scalar_input():
+    """Test create series impedance using scalar input."""
+    Z = 42
+    tmat = TransmissionMatrix.create_series_impedance(Z)
+    assert isinstance(tmat, np.ndarray)
+    assert tmat.shape == (2,2)
+    npt.assert_allclose(tmat, [[1, Z],[0, 1]], atol=1e-15)
+
+@pytest.mark.parametrize("abcd_cshape", [(1,), (4,5)])
+def test_tmatrix_create_shunt_admittance(A_FreqDat, abcd_cshape):
+    """Test `create_shunt_admittance` using FrequencyData input."""
+    Y = pf.utils.broadcast_cshape(A_FreqDat, abcd_cshape)
+    tmat = TransmissionMatrix.create_shunt_admittance(Y)
+    assert isinstance(tmat, TransmissionMatrix)
+    assert tmat.abcd_cshape == abcd_cshape
+    _compare_tmat_vs_abcd(tmat, 1, 0, Y.freq, 1)
+
+def test_tmatrix_create_shunt_admittance_scalar_input():
+    """Test `create_shunt_admittance` using scalar input."""
+    Y = 42
+    tmat = TransmissionMatrix.create_shunt_admittance(Y)
+    assert isinstance(tmat, np.ndarray)
+    assert tmat.shape == (2,2)
+    npt.assert_allclose(tmat, [[1, 0],[Y, 1]], atol=1e-15)
+
+@pytest.mark.parametrize("transducer_constant",
+                         [2.5, (2.5), FrequencyData([2.5, 5, 10], [1, 2, 3])])
+def test_tmatrix_create_transformer(transducer_constant, frequencies):
+    """Test `create_transformer` for FrequencyData and scalar input."""
+    tmat = TransmissionMatrix.create_transformer(transducer_constant)
+    if isinstance(transducer_constant, FrequencyData):
+        assert(isinstance(tmat, TransmissionMatrix))
+        N = transducer_constant.freq
+    else:
+        assert(isinstance(tmat, np.ndarray))
+        N = transducer_constant
+        # Convert to T-Matrix object
+        tmat = TransmissionMatrix.create_identity(frequencies) @ tmat
+
+    Zl = 100
+    Zin_expected = N*N * Zl
+    Zin = tmat.input_impedance(Zl)
+    npt.assert_allclose(Zin.freq, Zin_expected, atol = 1e-15)
+
+@pytest.mark.parametrize("transducer_constant",
+                         [2.5, (2.5), FrequencyData([2.5, 5, 10], [1, 2, 3])])
+def test_tmatrix_create_gyrator(transducer_constant, frequencies):
+    """Test `create_gyrator` for FrequencyData and scalar input."""
+    tmat = TransmissionMatrix.create_gyrator(transducer_constant)
+    if isinstance(transducer_constant, FrequencyData):
+        assert(isinstance(tmat, TransmissionMatrix))
+        N = transducer_constant.freq
+    else:
+        assert(isinstance(tmat, np.ndarray))
+        N = transducer_constant
+        # Convert to T-Matrix object
+        tmat = TransmissionMatrix.create_identity(frequencies) @ tmat
+
+    Zl = 100
+    Zin_expected = N*N / Zl
+    Zin = tmat.input_impedance(Zl)
+    npt.assert_allclose(Zin.freq, Zin_expected, atol = 1e-15)
+
+
+def test_tmatrix_slicing(frequencies):
+    """Test whether slicing a T-Matrix object return T-Matrix or raises correct
+    error for invalid keys."""
+    eye_2x2 = TransmissionMatrix.create_identity(frequencies)
+    eye_1x2x2 = pf.utils.broadcast_cshape(eye_2x2, (1, 2, 2))
+    eye_3x2x2 = pf.utils.broadcast_cshape(eye_2x2, (3, 2, 2))
+    eye_4x3x2x2 = pf.utils.broadcast_cshape(eye_2x2, (4, 3, 2, 2))
+
+    npt.assert_allclose(eye_1x2x2[0].freq, eye_2x2.freq, atol = 1e-15)
+    npt.assert_allclose(eye_3x2x2[0].freq, eye_2x2.freq, atol = 1e-15)
+    npt.assert_allclose(eye_3x2x2[1].freq, eye_2x2.freq, atol = 1e-15)
+    npt.assert_allclose(eye_3x2x2[2].freq, eye_2x2.freq, atol = 1e-15)
+    npt.assert_allclose(eye_4x3x2x2[0,0].freq, eye_2x2.freq, atol = 1e-15)
+    npt.assert_allclose(eye_4x3x2x2[1,2].freq, eye_2x2.freq, atol = 1e-15)
+
+    error_msg = "Object is not indexable, since ABCD-entries " \
+    "only have a single channel"
+    with pytest.raises(IndexError, match=error_msg):
+        eye_2x2[0]
+
+    error_msg = "Indexed dimensions must not exceed the ABCD " \
+    "channel dimension (abcd_cdim), which is "
+    with pytest.raises(IndexError, match=re.escape(error_msg + "1")):
+        eye_1x2x2[0,1]
+    with pytest.raises(IndexError, match=re.escape(error_msg + "1")):
+        eye_1x2x2[0,:]
+    with pytest.raises(IndexError, match=re.escape(error_msg + "2")):
+        eye_4x3x2x2[0,0,0]
+    with pytest.raises(IndexError, match=re.escape(error_msg + "2")):
+        eye_4x3x2x2[0,:,0]
+    with pytest.raises(IndexError, match=re.escape(error_msg + "2")):
+        eye_4x3x2x2[0,:,:]
diff --git a/tests/test_transmission_matrix_derived_parameters.py b/tests/test_transmission_matrix_derived_parameters.py
new file mode 100644
index 00000000..ea639243
--- /dev/null
+++ b/tests/test_transmission_matrix_derived_parameters.py
@@ -0,0 +1,303 @@
+import numpy as np
+import numpy.testing as npt
+import pytest
+import re
+from pyfar import TransmissionMatrix
+from pyfar import FrequencyData
+
+
+@pytest.fixture(scope="module")
+def frequencies():
+    return [100,200,300]
+
+#---------------
+#| INPUT TESTS |
+#---------------
+@pytest.fixture(scope="module", params=
+                ["input_impedance", "output_impedance",
+                 "TF_voltage", "TF_current"])
+def parameter_function(request):
+    """Parametrized fixture with names of methods for derived parameters."""
+    return request.param
+@pytest.fixture(scope="module", params=
+                ["abcd_cshape()", "abcd_cshape(4,)", "abcd_cshape(4,5)"])
+def abcd_cshape(request):
+    """Parametrized fixture for different abcd_cshapes."""
+    if request.param == "abcd_cshape()":
+        return ()
+    elif request.param == "abcd_cshape(4,)":
+        return (4,)
+    elif request.param == "abcd_cshape(4,5)":
+        return (4,5)
+@pytest.fixture(scope="module")
+def tmatrix_random_data(abcd_cshape, frequencies):
+    """Fixture creating T-Matrix with rng data for different abcd_cshapes."""
+    shape = (*abcd_cshape, 2,2, len(frequencies))
+    rng = np.random.default_rng()
+    data = rng.uniform(0.0001, 2000, shape)
+    return TransmissionMatrix(data, frequencies)
+
+@pytest.fixture(scope="module", params=
+                ["scalar", "FrequencyData_vector",
+                 "FrequencyData_abcd_cshape"])
+def load_impedance_with_correct_format(request, frequencies, abcd_cshape):
+    """Fixture creating load_impedance in different formats:
+    1) scalar,
+    2) array matching the frequency vector,
+    3) FrequencyData object match abcd_cshape."""
+    if request.param == "scalar":
+        return 2
+    elif request.param == "FrequencyData_vector":
+        return FrequencyData(np.ones_like(frequencies)*2, frequencies)
+    elif request.param == "FrequencyData_abcd_cshape":
+        shape = (*abcd_cshape, len(frequencies))
+        return FrequencyData(np.ones(shape), frequencies)
+
+@pytest.mark.parametrize("impedance_type", ["input", "output"])
+def test_valid_input_formats_impedance(impedance_type,
+                                         load_impedance_with_correct_format,
+                                         tmatrix_random_data):
+    """Test whether input/output impedance runs without errors for Zl with
+    valid input formats."""
+    Zl = load_impedance_with_correct_format
+    if impedance_type == "input":
+        tmatrix_random_data.input_impedance(Zl)
+    elif impedance_type == "output":
+        tmatrix_random_data.output_impedance(Zl)
+
+@pytest.fixture(scope="module")
+def simple_tmat():
+    return TransmissionMatrix(np.ones([2,2,1]), 100)
+@pytest.mark.parametrize("quantity_indices",
+                         [(0,0), (0,1), (1,0), (1,1), [0,0], np.array([1,0])])
+def test_TF_valid_quantity_input(quantity_indices, simple_tmat):
+    """Test whether TF method runs without errors for valid input types
+    for indices."""
+    simple_tmat.transfer_function(quantity_indices, np.inf)
+
+@pytest.mark.parametrize("quantity_indices",
+                         [1, "string", (1, "string"), [1,1,1]])
+def test_TF_quantity_input_wrong_numel(quantity_indices, simple_tmat):
+    """Test whether TF method raises error for invalid types or size
+    of indices."""
+    error_msg = re.escape("'quantity_indices' must be an array-like type "
+                          "with two numeric elements.")
+    with pytest.raises(ValueError, match=error_msg):
+        simple_tmat.transfer_function(quantity_indices, np.inf)
+
+@pytest.mark.parametrize("quantity_indices", [(-1,0), (0, 1.9)])
+def test_TF_quantity_input_wrong_ints(quantity_indices, simple_tmat):
+    """Test whether TF method raises error for invalid integer values
+    for indices."""
+    error_msg = re.escape("'quantity_indices' must contain two "
+                          "integers between 0 and 1.")
+    with pytest.raises(ValueError, match=error_msg):
+        simple_tmat.transfer_function(quantity_indices, np.inf)
+
+def test_TF_valid_load_input_format(load_impedance_with_correct_format,
+                                            tmatrix_random_data):
+    """Test whether TF method runs without errors for valid Zl input."""
+    Zl = load_impedance_with_correct_format
+    tmatrix_random_data.transfer_function((1,1), Zl)
+
+
+#---------------
+#| RESULT TESTS|
+#---------------
+@pytest.fixture(scope="module")
+def impedance_random(frequencies) -> FrequencyData:
+    """Fixture returning impedance as FrequencyData object with random data."""
+    rng = np.random.default_rng()
+    return FrequencyData(rng.random(len(frequencies)), frequencies)
+
+@pytest.fixture(scope="module",
+                params=["random_load", "inf_load", "zero_load", "mixed_load"])
+def load_impedance(request, frequencies) -> FrequencyData:
+    """Parametrized fixture for load impedance as FrequencyData object with
+    different types of data:
+    1) random, 2) infinite load, 3) zero load, 4) mix of all."""
+    if request.param == "random_load":
+        rng = np.random.default_rng()
+        return FrequencyData(rng.random(len(frequencies)), frequencies)
+    elif request.param == "inf_load":
+        return FrequencyData(np.ones_like(frequencies)*np.inf, frequencies)
+    elif request.param == "zero_load":
+        return FrequencyData(np.zeros_like(frequencies), frequencies)
+    elif request.param == "mixed_load":
+        return FrequencyData([0, 2, np.inf], frequencies)
+
+def _special_twoport_tmatrix(
+        twoport_type, Zl: FrequencyData, Z: FrequencyData):
+    """Returns a T-Matrix representing a twoport for special circuit types:
+    1) By-pass system
+    2) A series impedance != load_impedance
+    3) A parallel impedance != load_impedance
+    4) A series impedance == load_impedance
+    5) A parallel impedance == load_impedance"""
+    if twoport_type == "bypass":
+        return TransmissionMatrix.create_identity(Zl.frequencies)
+    if twoport_type == "series_impedance":
+        return TransmissionMatrix.create_series_impedance(Z)
+    if twoport_type == "parallel_impedance":
+        return TransmissionMatrix.create_shunt_admittance(1/Z)
+    if twoport_type == "series_load_impedance":
+        return TransmissionMatrix.create_series_impedance(Zl)
+    if twoport_type == "parallel_load_impedance":
+        return TransmissionMatrix.create_shunt_admittance(1/Zl)
+    else:
+        raise ValueError("Unexpected value for 'twoport_type'")
+
+def _twoport_type_list():
+    """A list including all valid twoport types used in the tests."""
+    return ["bypass", "series_impedance", "parallel_impedance",
+            "series_load_impedance", "parallel_load_impedance"]
+
+def _expected_impedance(
+        twoport_type, Zl: FrequencyData, Z: FrequencyData) -> FrequencyData:
+    """Returns the expected input/output impedance for testes twoport types."""
+    if twoport_type == "bypass":
+        Zexpected =  Zl
+    elif twoport_type == "series_impedance":
+        Zexpected =  Zl + Z
+    elif twoport_type == "parallel_impedance":
+        Zexpected =  1/(1/Zl + 1/Z)
+    elif twoport_type == "series_load_impedance":
+        Zexpected =  2 * Zl
+        Zexpected.freq[np.isinf(Zl.freq)] = np.nan #non-phyical case
+    elif twoport_type == "parallel_load_impedance":
+        Zexpected =  Zl / 2
+        Zexpected.freq[Zl.freq == 0] = np.nan #non-phyical case
+    else:
+        raise ValueError("Unexpected value for 'twoport_type'")
+    return Zexpected
+
+def _expected_voltage_tf(
+        twoport_type, Zl: FrequencyData, Z: FrequencyData) -> FrequencyData:
+    """Returns the expected TF Uout/Uin for testes twoport types."""
+    if twoport_type not in _twoport_type_list():
+        raise ValueError("Unexpected value for 'twoport_type'")
+
+    if twoport_type == "series_impedance":
+        TF = Zl / (Zl + Z)
+    else:
+        TF = FrequencyData(np.ones_like(Zl.freq), Zl.frequencies)
+
+    # Special TFs for zero of inf load
+    TF.freq[Zl.freq == 0] = 0
+    TF.freq[np.isinf(Zl.freq)] = 1
+
+    if twoport_type == "series_load_impedance":
+        TF = TF * (1/2)
+        TF.freq[np.isinf(Zl.freq)] = np.nan
+    return TF
+
+def _expected_current_to_voltage_tf(
+        twoport_type, Zl: FrequencyData, Z: FrequencyData) -> FrequencyData:
+    """Returns the expected TF Iout/Uin for testes twoport types."""
+    if twoport_type not in _twoport_type_list():
+        raise ValueError("Unexpected value for 'twoport_type'")
+
+    if twoport_type == "series_impedance":
+        TF = 1 / (Zl + Z)
+    elif twoport_type == "series_load_impedance":
+        TF = 1 / (2 * Zl)
+        TF.freq[Zl.freq == 0] = np.nan
+    else:
+        TF = 1 / Zl
+        TF.freq[Zl.freq == 0] = np.nan
+
+    return TF
+
+def _expected_current_tf(
+        twoport_type, Zl: FrequencyData, Z: FrequencyData) -> FrequencyData:
+    """Returns the expected TF Iout/Iin for testes twoport types."""
+    if twoport_type not in _twoport_type_list():
+        raise ValueError("Unexpected value for 'twoport_type'")
+
+    if twoport_type == "parallel_impedance":
+        TF = Z / (Zl + Z)
+    else:
+        TF = FrequencyData(np.ones_like(Zl.freq), Zl.frequencies)
+
+    # Special TFs for zero of inf load
+    TF.freq[Zl.freq == 0] = 1
+    TF.freq[np.isinf(Zl.freq)] = 0
+
+    if twoport_type == "parallel_load_impedance":
+        TF = TF * (1/2)
+        TF.freq[Zl.freq == 0] = np.nan
+    return TF
+
+def _expected_voltage_to_current_tf(
+        twoport_type, Zl: FrequencyData, Z: FrequencyData) -> FrequencyData:
+    """Returns the expected TF Uout/Iin for testes twoport types."""
+    if twoport_type not in _twoport_type_list():
+        raise ValueError("Unexpected value for 'twoport_type'")
+
+    idx_inf = np.isinf(Zl.freq)
+    if twoport_type == "parallel_impedance":
+        TF = Z * Zl / (Zl + Z)
+        TF.freq[idx_inf] = Z.freq[idx_inf]
+    elif twoport_type == "parallel_load_impedance":
+        TF = Zl / 2
+        TF.freq[Zl.freq == 0] = np.nan
+        TF.freq[idx_inf] = np.nan
+    else: #bypass, series_impedance + series_load_impedance
+        TF = Zl.copy()
+        TF.freq[idx_inf] = np.nan
+
+    return TF
+
+@pytest.mark.parametrize("impedance_type", ["input", "output"])
+@pytest.mark.parametrize("twoport_type", _twoport_type_list())
+def test_input_impedance(impedance_type : str, twoport_type : str,
+                         load_impedance : FrequencyData, impedance_random):
+    """Significantly parametrized test for result of input/output impedance
+    method."""
+    tmat = _special_twoport_tmatrix(
+        twoport_type, load_impedance, impedance_random)
+    if impedance_type == "input":
+        Zres = tmat.input_impedance(load_impedance)
+    else:
+        Zres = tmat.output_impedance(load_impedance)
+    Zexpected = _expected_impedance(
+        twoport_type, load_impedance, impedance_random)
+
+    idx_inf = Zexpected.freq == np.inf
+    idx_default = np.logical_not(idx_inf)
+    npt.assert_allclose(
+        Zres.freq[idx_default], Zexpected.freq[idx_default], atol=1e-15)
+    assert(np.all(np.abs(Zres.freq[idx_inf]) > 1e15))
+
+@pytest.mark.parametrize("tf_type", ["voltage", "current", "voltage/current",
+                                     "current/voltage"])
+@pytest.mark.parametrize("twoport_type", _twoport_type_list())
+def test_transfer_function(tf_type, twoport_type,
+                           load_impedance, impedance_random):
+    """Significantly parametrized test for result of TF method."""
+    tmat = _special_twoport_tmatrix(
+        twoport_type, load_impedance, impedance_random)
+
+    if tf_type == "voltage":
+        quantity_indices = (0,0)
+        TF_expected = _expected_voltage_tf(twoport_type,
+                                           load_impedance, impedance_random)
+    elif tf_type == "current":
+        quantity_indices = (1,1)
+        TF_expected = _expected_current_tf(twoport_type,
+                                           load_impedance, impedance_random)
+    elif tf_type == "voltage/current":
+        quantity_indices = (0,1)
+        TF_expected = _expected_voltage_to_current_tf(
+            twoport_type, load_impedance, impedance_random)
+    elif tf_type == "current/voltage":
+        quantity_indices = (1,0)
+        TF_expected = _expected_current_to_voltage_tf(
+            twoport_type, load_impedance, impedance_random)
+    else:
+        raise ValueError("Unexpected value for 'tf_type'.")
+
+    TF = tmat.transfer_function(quantity_indices, load_impedance)
+    if isinstance(TF_expected, FrequencyData):
+        TF_expected = TF_expected.freq
+    npt.assert_allclose(TF.freq, TF_expected, atol=1e-15)
diff --git a/tests/test_transmission_matrix_multiplication.py b/tests/test_transmission_matrix_multiplication.py
new file mode 100644
index 00000000..3d70a172
--- /dev/null
+++ b/tests/test_transmission_matrix_multiplication.py
@@ -0,0 +1,101 @@
+import numpy as np
+import numpy.testing as npt
+import pytest
+from pyfar import TransmissionMatrix
+
+
+def _abcd_matrix_stack(abcd, shape_extra_dims=None,
+                       frequencies=[100, 200, 300]):
+    """Create T-Matrix with given extra-dimensions for testing purposes."""
+    n_freq_bins = len(frequencies)
+    if shape_extra_dims is None:
+        shape = n_freq_bins
+    else:
+        shape = np.append(shape_extra_dims, n_freq_bins)
+
+    A = np.ones(shape, dtype=type(abcd[0])) * abcd[0]
+    B = np.ones(shape, dtype=type(abcd[1])) * abcd[1]
+    C = np.ones(shape, dtype=type(abcd[2])) * abcd[2]
+    D = np.ones(shape, dtype=type(abcd[3])) * abcd[3]
+    return TransmissionMatrix.from_abcd(A, B, C, D, frequencies)
+
+def _check_matrix_multiplication_result(
+    abcd_mat1,
+    abcd_mat2,
+    abcd_target,
+    shape_extra_dims=None,
+    frequencies=[100, 200, 300],
+):
+    """Assertion check if matrix multiplication result matches given target."""
+    tmat1 = _abcd_matrix_stack(abcd_mat1, shape_extra_dims, frequencies)
+    tmat2 = _abcd_matrix_stack(abcd_mat2, shape_extra_dims, frequencies)
+    tmat_out = tmat1 @ tmat2
+    npt.assert_allclose(tmat_out.A.freq, abcd_target[0], atol=1e-15)
+    npt.assert_allclose(tmat_out.B.freq, abcd_target[1], atol=1e-15)
+    npt.assert_allclose(tmat_out.C.freq, abcd_target[2], atol=1e-15)
+    npt.assert_allclose(tmat_out.D.freq, abcd_target[3], atol=1e-15)
+
+
+@pytest.mark.parametrize("shape_extra_dims", [None, 4, [4, 5]])
+def test_tmatrix_multiplication_identity(
+    shape_extra_dims, frequencies=[100, 200, 300]):
+    """Test result of T-Matrix multiplication with eye matrix."""
+    rng = np.random.default_rng()
+    abcd_identity = (1, 0, 0, 1)
+    abcd_rng = rng.random(4)
+    _check_matrix_multiplication_result(
+        abcd_rng, abcd_identity, abcd_rng, shape_extra_dims, frequencies)
+
+@pytest.mark.parametrize("shape_extra_dims", [None, 4, [4, 5]])
+def test_tmatrix_multiplication_negative_identity(
+    shape_extra_dims, frequencies=[100, 200, 300]):
+    """Test result of T-Matrix multiplication with negated eye matrix."""
+    rng = np.random.default_rng()
+    abcd_identity = (-1, 0, 0, -1)
+    abcd_rng = rng.random(4)
+    _check_matrix_multiplication_result(
+        abcd_rng, abcd_identity, -abcd_rng, shape_extra_dims, frequencies)
+
+@pytest.mark.parametrize("shape_extra_dims", [None, 4, [4, 5]])
+def test_tmatrix_multiplication_bottom_row_zero(
+    shape_extra_dims, frequencies=[100, 200, 300]):
+    """
+    Inputs:
+        [ 2  1 ]    [ a  b ]
+        [ 0  0 ]    [ c  d ]
+    Expected Output:
+        [2a+c  2b+d]
+        [  0     0 ]
+    """
+    rng = np.random.default_rng()
+    abcd_ApluB = (2, 1, 0, 0)
+    abcd_rng = rng.random(4)
+    abcd_target = (2 * abcd_rng[0] + abcd_rng[2],
+                   2 * abcd_rng[1] + abcd_rng[3], 0, 0)
+    _check_matrix_multiplication_result(
+        abcd_ApluB, abcd_rng, abcd_target, shape_extra_dims, frequencies)
+
+@pytest.mark.parametrize("shape_extra_dims", [None, 4, [4, 5]])
+def test_tmatrix_multiplication_random(
+    shape_extra_dims, frequencies=[100, 200, 300]):
+    """
+    Inputs:
+        [ a1  b1 ]    [ a2  b2 ]
+        [ c1  d1 ]    [ c2  d2 ]
+    Expected Output:
+        [a1*a2 + b1*c2   a1*b2 + b1*d2]
+        [c1*a2 + d1*c2   c1*b2 + d1*d2]
+    """
+    rng = np.random.default_rng()
+    abcd_rng1 = rng.random(4)
+    abcd_rng2 = rng.random(4)
+    a1, b1, c1, d1 = abcd_rng1[0], abcd_rng1[1], abcd_rng1[2], abcd_rng1[3]
+    a2, b2, c2, d2 = abcd_rng2[0], abcd_rng2[1], abcd_rng2[2], abcd_rng2[3]
+    abcd_target = (
+        a1 * a2 + b1 * c2,
+        a1 * b2 + b1 * d2,
+        c1 * a2 + d1 * c2,
+        c1 * b2 + d1 * d2,
+    )
+    _check_matrix_multiplication_result(
+        abcd_rng1, abcd_rng2, abcd_target, shape_extra_dims, frequencies)