Adds Real Spherical Harmonics Directivity (#422)

* Add directivity pattern using real spherical harmonics.

* Adds a test for the real spherical harmonics directivity.

---------

Co-authored-by: Erik Fleischhauer <fleischhauer@iks.rwth-aachen.de>

Author: je-fleischhauer

CHANGELOG.rst

@@ -19,6 +19,8 @@ Changed
 - Improves ``pra.experimental.measure_rt60``: Compute the T60 using a fit. Default is
   log-domain. Adds option to fit in linear domain.
 
+- Add Directivity pattern for real spherical harmonics. 
+
 Bugfix
 ~~~~~~
 

pyroomacoustics/directivities/__init__.py

@@ -92,4 +92,5 @@
 )
 from .base import Directivity
 from .direction import DirectionVector, Rotation3D
+from .harmonics import RealSphericalHarmonicsDirectivity
 from .measured import MeasuredDirectivity, MeasuredDirectivityFile

pyroomacoustics/directivities/harmonics.py

@@ -0,0 +1,137 @@
+# Copyright (C) 2026  Erik Fleischhauer
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+#
+# You should have received a copy of the MIT License along with this program. If
+# not, see <https://opensource.org/licenses/MIT>.
+import numpy as np
+import scipy.special
+
+from .base import Directivity
+
+
+def get_mn_in_acn_order(order):
+    """
+    Calculates the (m,n) pairs in ACN order up to a given order.
+
+    Parameters:
+    ----------
+    order : int
+        Maximum degree of the spherical harmonics.
+
+    Returns:
+    ----------
+    all_m : ndarray
+        Array of orders m in ACN order.
+    all_n : ndarray
+        Array of degrees n in ACN order.
+    """
+    all_m = np.array([j - i for i in range(0, order + 1) for j in range(0, 2 * i + 1)])
+    all_n = np.array([i for i in range(0, order + 1) for _ in range(0, 2 * i + 1)])
+    return all_m, all_n
+
+
+def real_sph_harm(n, m, theta, phi, condon_shortley_phase=False):
+    """
+    Calculates the real spherical harmonics.
+
+    Parameters:
+    ----------
+    n : int
+        Degree of the spherical harmonic.
+    m : int
+        Order of the spherical harmonic.
+    theta : array_like
+        Polar (colatitudinal) coordinate in radians.
+    phi : array_like
+        Azimuthal coordinate in radians.
+    condon_shortley_phase : bool, optional
+        If True, includes the Condon-Shortley phase factor (-1)^m. Default is False.
+
+    Returns:
+    ----------
+    y_real : ndarray
+        Real spherical harmonics evaluated at the given angles.
+    """
+    m = np.atleast_1d(m)
+
+    try:
+        ysh_complx = scipy.special.sph_harm_y(n, m, theta, phi)
+    except AttributeError:
+        # Deprecated since scipy v1.15.0.
+        ysh_complx = scipy.special.sph_harm(m, n, phi, theta)
+
+    y_real = np.empty_like(ysh_complx, dtype=np.float64)
+    y_real[(m >= 0)] = np.real(ysh_complx[(m >= 0)])
+    y_real[(m < 0)] = np.imag(ysh_complx[(m < 0)])
+
+    if not condon_shortley_phase:
+        # Cancel Condon-Shortley Phase (term (-1) ** m) by multiplying with (-1) ** m
+        # In Rafaely's book, this step is not done
+        y_real[(m != 0)] *= np.sqrt(2) * np.array([-1.0]) ** m[(m != 0)]
+        # In the formular, for m<0, |m| is used. We consider this by the use of the constraint.
+        y_real[np.logical_and(m < 0, (m % 2) == 0)] = -y_real[
+            np.logical_and(m < 0, (m % 2) == 0)
+        ]
+
+    return y_real
+
+
+class RealSphericalHarmonicsDirectivity(Directivity):
+    """
+    A class for real spherical harmonic directivity patterns.
+
+    Parameters:
+    ----------
+    m: int
+        Order of the spherical harmonic.
+    n: int
+        Degree of the spherical harmonic.
+    condon_shortley_phase: bool, optional
+        If True, includes the Condon-Shortley phase factor (-1)^m. Default is False.
+    """
+
+    def __init__(self, m, n, condon_shortley_phase: bool = False):
+        self.m = m
+        self.n = n
+
+        self.condon_shortley_phase = condon_shortley_phase
+
+    @property
+    def is_impulse_response(self):
+        return True
+
+    @property
+    def filter_len_ir(self):
+        return 1
+
+    def get_response(
+        self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True
+    ):
+
+        if degrees:
+            azimuth = np.radians(azimuth)
+            colatitude = np.radians(colatitude)
+        return real_sph_harm(
+            self.m,
+            self.n,
+            colatitude,
+            azimuth,
+            condon_shortley_phase=self.condon_shortley_phase,
+        )[:, np.newaxis]

pyroomacoustics/directivities/tests/test_harmonics_directivities.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+import pyroomacoustics as pra
+from pyroomacoustics.directivities.harmonics import (
+    get_mn_in_acn_order,
+    real_sph_harm,
+)
+
+
+def test_harmonics_directivity():
+    """
+    Test the Room Impulse Response (RIR) generated by a set of microphones
+    with real spherical harmonics directivities against the analytical values
+    of the spherical harmonics.
+    """
+    order = 2
+    condon_shortley_phase = True
+    azimuth = np.deg2rad(50)
+
+    room = pra.AnechoicRoom(fs=16000)
+    room.add_source([np.cos(azimuth), np.sin(azimuth), 1.0])
+
+    for m, n in zip(*get_mn_in_acn_order(order)):
+        room.add_microphone(
+            [0.0, 0.0, 1.0],
+            directivity=pra.directivities.RealSphericalHarmonicsDirectivity(
+                m, n, condon_shortley_phase=condon_shortley_phase
+            ),
+        )
+
+    room.compute_rir()
+    rir_sum = np.asarray(room.rir).squeeze().sum(axis=1)
+
+    m, n = get_mn_in_acn_order(order)
+    sh = real_sph_harm(
+        n, m, np.pi / 2, azimuth, condon_shortley_phase=condon_shortley_phase
+    )
+
+    np.allclose(rir_sum, sh, atol=1e-2)