Adds the ray sampler for measured directivities.

pyroomacoustics/directivities/analytic.py

@@ -413,14 +413,17 @@ 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,
             gain=1.0,
             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):

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
         )
 

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()

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)