From 9cf7c0345a0441e3a8f59ee6e84c924e12bd518d Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Fri, 8 Nov 2024 23:53:48 +0900 Subject: [PATCH] Adds the ray sampler for measured directivities. --- pyroomacoustics/directivities/analytic.py | 5 ++- pyroomacoustics/directivities/measured.py | 37 ++++++++++++++++++- pyroomacoustics/doa/histogram.py | 4 -- .../random/tests/test_rejection.py | 8 +++- 4 files changed, 46 insertions(+), 8 deletions(-) diff --git a/pyroomacoustics/directivities/analytic.py b/pyroomacoustics/directivities/analytic.py index 5fd5ba19..c7102089 100644 --- a/pyroomacoustics/directivities/analytic.py +++ b/pyroomacoustics/directivities/analytic.py @@ -413,7 +413,7 @@ def __init__(self, loc=None, p=0.5): self._coeff = p def _pattern(self, x): - return cardioid_func( + response = cardioid_func( x.T, direction=self._loc, p=self._coeff, @@ -421,6 +421,9 @@ def _pattern(self, x): normalize=False, magnitude=True, ) + # The number of rays needs to be proportional to the + # response energy. + return response**2 def cardioid_func(x, direction, p, gain=1.0, normalize=True, magnitude=False): diff --git a/pyroomacoustics/directivities/measured.py b/pyroomacoustics/directivities/measured.py index 9c63e097..3ecc28cc 100644 --- a/pyroomacoustics/directivities/measured.py +++ b/pyroomacoustics/directivities/measured.py @@ -82,6 +82,35 @@ from .direction import Rotation3D from .interp import spherical_interpolation from .sofa import open_sofa_file +from .. import random + + +class MeasuredDirectivitySampler(random.sampler.DirectionalSampler): + """ + This object draws samples from the distribution defined by the energy + of the measured directional response object. + + Parameters + ---------- + loc: array_like + The unit vector pointing in the main direction of the cardioid + p: float + Parameter of the cardioid pattern. A value of 0 corresponds to a + figure-eight pattern, 0.5 to a cardioid pattern, and 1 to an omni + pattern + The parameter must be between 0 and 1 + """ + + def __init__(self, kdtree, energy): + super().__init__() + self._kdtree = kdtree + # Normalize to maximum energy 1 because that is also the + # maximum value of the proposal unnormalized uniform distribution. + self._energy = energy / energy.max() + + def _pattern(self, x): + _, index = self._kdtree.query(x) + return self._energy[index] class MeasuredDirectivity(Directivity): @@ -144,6 +173,10 @@ def set_orientation(self, orientation): # create the kd-tree self._kdtree = cKDTree(self._grid.cartesian.T) + # create the ray sampler + ir_energy = np.square(self._irs).mean(axis=-1) + self._ray_sampler = MeasuredDirectivitySampler(self._kdtree, ir_energy) + def get_response( self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True ): @@ -173,7 +206,7 @@ def get_response( return self._irs[index, :] def sample_rays(self, n_rays): - raise NotImplementedError() + return self._ray_sampler(n_rays).T @requires_matplotlib def plot(self, freq_bin=0, n_grid=100, ax=None, depth=False, offset=None): @@ -241,7 +274,7 @@ def plot(self, freq_bin=0, n_grid=100, ax=None, depth=False, offset=None): Y *= V Z *= V - surf = ax.plot_surface( + _ = ax.plot_surface( X, Y, Z, facecolors=cmap.to_rgba(V), linewidth=0, antialiased=False ) diff --git a/pyroomacoustics/doa/histogram.py b/pyroomacoustics/doa/histogram.py index 84e97750..5c9c4cf4 100644 --- a/pyroomacoustics/doa/histogram.py +++ b/pyroomacoustics/doa/histogram.py @@ -43,10 +43,6 @@ def __init__(self, n_bins, dim=3, enable_peak_finding=False): else: raise NotImplementedError("Only 3D histogram has been implemented") - import pdb - - pdb.set_trace() - # we need to know the area of each bin self._voronoi = SphericalVoronoi(self._grid.cartesian.T) self._areas = self._voronoi.calculate_areas() diff --git a/pyroomacoustics/random/tests/test_rejection.py b/pyroomacoustics/random/tests/test_rejection.py index 6166f4b4..ee989dbe 100644 --- a/pyroomacoustics/random/tests/test_rejection.py +++ b/pyroomacoustics/random/tests/test_rejection.py @@ -23,7 +23,13 @@ # Figure-of-eight sampler = CardioidFamilySampler(loc=loc, p=0) - points = sampler(size=10000).T # shape (n_dim, n_points) + # Measured eigenmike response + eigenmike = pra.MeasuredDirectivityFile("EM32_Directivity", fs=16000) + rot_54_73 = pra.Rotation3D([73, 54], "yz", degrees=True) + dir_obj_Emic = eigenmike.get_mic_directivity("EM_32_9", orientation=rot_54_73) + sampler = dir_obj_Emic._ray_sampler + + points = sampler(size=100000).T # shape (n_dim, n_points) # Create a spherical histogram hist = pra.doa.SphericalHistogram(n_bins=500)