diff --git a/CHANGELOG.rst b/CHANGELOG.rst index ea7a0efe..5abbbc25 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -11,6 +11,11 @@ adheres to `Semantic Versioning `_. `Unreleased`_ ------------- +Added +~~~~~~~ + +- Added ``orientation`` parameter to ``RealSphericalHarmonicsDirectivity`` object. + Fixed ~~~~~ diff --git a/pyroomacoustics/directivities/harmonics.py b/pyroomacoustics/directivities/harmonics.py index bd5f6be5..fec01721 100644 --- a/pyroomacoustics/directivities/harmonics.py +++ b/pyroomacoustics/directivities/harmonics.py @@ -52,8 +52,9 @@ import numpy as np import scipy.special -from ..doa import cart2spher +from ..doa import cart2spher, spher2cart from .base import Directivity +from .direction import Rotation3D def get_mn_in_acn_order(degree): @@ -140,16 +141,44 @@ class RealSphericalHarmonicsDirectivity(Directivity): Order of the spherical harmonic. n: int Degree of the spherical harmonic. + orientation: Rotation3D, optional + A rotation to apply to the pattern. If not provided, the + identity rotation is used, i.e., the harmonic is evaluated in + its canonical frame with no reorientation. 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): + def __init__(self, m, n, orientation=None, condon_shortley_phase: bool = False): self.m = m self.n = n - self.condon_shortley_phase = condon_shortley_phase + if orientation is None: + orientation = Rotation3D( + [0.0, 0.0, 0.0] + ) # identity rotation, no-op default + elif not isinstance(orientation, Rotation3D): + raise TypeError( + f"orientation must be a Rotation3D object, got {type(orientation).__name__}" + ) + self._orientation = orientation + + def set_orientation(self, orientation): + """ + Set orientation of directivity pattern. + + Parameters + ---------- + orientation: Rotation3D + New direction for the directivity pattern. + """ + if not isinstance(orientation, Rotation3D): + raise TypeError( + f"orientation must be a Rotation3D object, got {type(orientation).__name__}" + ) + self._orientation = orientation + @property def is_impulse_response(self): return True @@ -159,18 +188,10 @@ def filter_len_ir(self): return 1 def get_response_cartesian(self, directions, magnitude=False, frequency=None): - az, co, _ = cart2spher(directions.T) - return self.get_response( - az, co, magnitude=magnitude, frequency=frequency, degrees=False - ) + local_dirs = self._orientation.rotate_transpose(directions.T) - def get_response( - self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True - ): + azimuth, colatitude, _ = cart2spher(local_dirs) - if degrees: - azimuth = np.radians(azimuth) - colatitude = np.radians(colatitude) return real_sph_harm( self.n, self.m, @@ -179,6 +200,15 @@ def get_response( condon_shortley_phase=self.condon_shortley_phase, )[:, np.newaxis] + def get_response( + self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True + ): + + # world-frame cartesian unit vectors + directions = spher2cart(azimuth, colatitude, degrees=degrees).T + + return self.get_response_cartesian(directions, magnitude=magnitude) + def sample_rays(self, n_rays, rng=None): """Not yet implemented.""" raise NotImplementedError diff --git a/tests/directivities/test_harmonics_directivities.py b/tests/directivities/test_harmonics_directivities.py index 30ea9cc7..b82b134a 100644 --- a/tests/directivities/test_harmonics_directivities.py +++ b/tests/directivities/test_harmonics_directivities.py @@ -6,6 +6,25 @@ real_sph_harm, ) +FS = 16000 +MIC_POS = [1.0, 1.2, 1.5] +SRC_POS = [3.5, 2.8, 1.7] +ROOM_DIM = [6.0, 5.0, 3.0] + + +def _simulate_rir(directivity, max_order=0): + """Run a room simulation with a single directional microphone.""" + if max_order == 0: + room = pra.AnechoicRoom(fs=FS, dim=3) + else: + mat = pra.Material(energy_absorption=0.3) + room = pra.ShoeBox(ROOM_DIM, fs=FS, materials=mat, max_order=max_order) + + room.add_source(SRC_POS) + room.add_microphone(MIC_POS, directivity=directivity) + room.compute_rir() + return room.rir[0][0] + def test_harmonics_directivity(): """ @@ -38,3 +57,69 @@ def test_harmonics_directivity(): ) np.allclose(rir_sum, sh, atol=1e-2) + + +def test_harmonics_rotated_z(): + """ + Check that the `orientation` attribute of `RealSphericalHarmonicsDirectivity` + correctly rotates the pattern about the z-axis: an (m=1, n=1) harmonic + rotated by 90 degrees about z should produce the same RIR as an + unrotated (m=-1, n=1) harmonic at the same position. + """ + max_order = 2 + + rot_90_z = pra.directivities.Rotation3D([90.0], "z", degrees=True) + + dir_rotated = pra.directivities.RealSphericalHarmonicsDirectivity( + m=1, n=1, orientation=rot_90_z + ) + dir_unrotated = pra.directivities.RealSphericalHarmonicsDirectivity(m=-1, n=1) + + rir_rotated = _simulate_rir(dir_rotated, max_order=max_order) + rir_reference = _simulate_rir(dir_unrotated, max_order=max_order) + + np.testing.assert_allclose(rir_rotated, rir_reference, atol=1e-9, rtol=0) + + +def test_harmonics_rotated_y(): + """ + Check that the `orientation` attribute of `RealSphericalHarmonicsDirectivity` + correctly rotates the pattern about the y-axis: an (m=0, n=1) harmonic + rotated by 90 degrees about y should produce the same RIR as an + unrotated (m=1, n=1) harmonic at the same position. + """ + max_order = 2 + + rot_90_y = pra.directivities.Rotation3D([90.0], "y", degrees=True) + + dir_rotated = pra.directivities.RealSphericalHarmonicsDirectivity( + m=0, n=1, orientation=rot_90_y + ) + dir_unrotated = pra.directivities.RealSphericalHarmonicsDirectivity(m=1, n=1) + + rir_rotated = _simulate_rir(dir_rotated, max_order=max_order) + rir_reference = _simulate_rir(dir_unrotated, max_order=max_order) + + np.testing.assert_allclose(rir_rotated, rir_reference, atol=1e-9, rtol=0) + + +def test_harmonics_rotated_x(): + """ + Check that the `orientation` attribute of `RealSphericalHarmonicsDirectivity` + correctly rotates the pattern about the x-axis: an (m=0, n=1) harmonic + rotated by -90 degrees about x should produce the same RIR as an + unrotated (m=-1, n=1) harmonic at the same position. + """ + max_order = 2 + + rot_m90_x = pra.directivities.Rotation3D([-90.0], "x", degrees=True) + + dir_rotated = pra.directivities.RealSphericalHarmonicsDirectivity( + m=0, n=1, orientation=rot_m90_x + ) + dir_unrotated = pra.directivities.RealSphericalHarmonicsDirectivity(m=-1, n=1) + + rir_rotated = _simulate_rir(dir_rotated, max_order=max_order) + rir_reference = _simulate_rir(dir_unrotated, max_order=max_order) + + np.testing.assert_allclose(rir_rotated, rir_reference, atol=1e-9, rtol=0)