[WIP] Source/microphone directivity for ray tracing

examples/plot_directivity_2D.py

@@ -0,0 +1,32 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pyroomacoustics import dB
+from pyroomacoustics.directivities import DirectivityPattern, \
+    DirectionVector, Directivity
+
+ORIENTATION = DirectionVector(azimuth=0, colatitude=90, degrees=True)
+LOWER_GAIN = -20
+
+# plot each directivity
+angles = np.linspace(start=0, stop=360, num=361, endpoint=True)
+angles = np.radians(angles)
+
+# plot each pattern
+fig = plt.figure()
+ax = plt.subplot(111, projection="polar")
+for _dir in DirectivityPattern.values():
+
+    dir_obj = Directivity(orientation=ORIENTATION, pattern=_dir)
+    gains = []
+    for a in angles:
+        direction = DirectionVector(azimuth=a, degrees=False)
+        gains.append(dir_obj.get_response(direction))
+    gains_db = dB(np.array(gains))
+    ax.plot(angles, gains_db, label=_dir)
+
+plt.legend(bbox_to_anchor=(1, 1))
+plt.ylim([LOWER_GAIN, 0])
+ax.yaxis.set_ticks(np.arange(start=LOWER_GAIN, stop=5, step=5))
+plt.tight_layout()
+plt.show()

examples/plot_directivity_3D.py

@@ -0,0 +1,19 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pyroomacoustics.directivities import DirectivityPattern, \
+    DirectionVector, Directivity
+
+
+PATTERN = DirectivityPattern.HYPERCARDIOID
+ORIENTATION = DirectionVector(azimuth=0, colatitude=45, degrees=True)
+
+# create cardioid object
+dir_obj = Directivity(orientation=ORIENTATION, pattern=PATTERN)
+
+# plot
+azimuth = np.linspace(start=0, stop=360, num=361, endpoint=True)
+colatitude = np.linspace(start=0, stop=180, num=180, endpoint=True)
+# colatitude = None   # for 2D plot
+dir_obj.plot_response(azimuth=azimuth, colatitude=colatitude, degrees=True)
+plt.show()

pyroomacoustics/__init__.py

@@ -130,6 +130,7 @@
 from .multirate import *
 from .acoustics import *
 from .recognition import *
+from .directivities import *
 
 from . import doa
 from . import adaptive
@@ -140,6 +141,7 @@
 from . import bss
 from . import denoise
 from . import phase
+from . import random
 
 import warnings
 warnings.warn(

pyroomacoustics/directivities.py

@@ -0,0 +1,183 @@
+import numpy as np
+from pyroomacoustics.utilities import Options
+
+
+class DirectivityPattern(Options):
+    FIGURE_EIGHT = "figure_eight"
+    HYPERCARDIOID = "hypercardioid"
+    CARDIOID = "cardioid"
+    SUBCARDIOID = "subcardioid"
+    OMNI = "omni"
+
+
+PATTERN_TO_CONSTANT = {
+    DirectivityPattern.FIGURE_EIGHT: 0,
+    DirectivityPattern.HYPERCARDIOID: 0.25,
+    DirectivityPattern.CARDIOID: 0.5,
+    DirectivityPattern.SUBCARDIOID: 0.75,
+    DirectivityPattern.OMNI: 1.0
+}
+
+
+class DirectionVector(object):
+    """
+    Parameters
+    ----------
+    azimuth : float
+    colatitude : float, optional
+        Default to PI / 2, only XY plane.
+    degrees : bool
+        Whether provided values are in degrees (True) or radians (False).
+    """
+    def __init__(self, azimuth, colatitude=None, degrees=True):
+        if degrees is True:
+            azimuth = np.radians(azimuth)
+            if colatitude is not None:
+                colatitude = np.radians(colatitude)
+        self._azimuth = azimuth
+        if colatitude is None:
+            colatitude = np.pi / 2
+        assert colatitude <= np.pi and colatitude >= 0
+        self._colatitude = colatitude
+
+    def get_azimuth(self, degrees=False):
+        if degrees:
+            return np.degrees(self._azimuth)
+        else:
+            return self._azimuth
+
+    def get_colatitude(self, degrees=False):
+        if degrees:
+            return np.degrees(self._colatitude)
+        else:
+            return self._colatitude
+
+
+class Directivity(object):
+    """
+
+    Parameters
+    ----------
+    orientation : DirectionVector
+        Indicates direction of the pattern.
+    pattern : str
+        One of directivities in `DirectivityPattern`.
+        TODO : support arbitrary directivities.
+
+    """
+    def __init__(self, orientation, pattern):
+        assert pattern in DirectivityPattern.values()
+        assert isinstance(orientation, DirectionVector)
+        self._orientation = orientation
+        self._pattern = pattern
+        self._p = PATTERN_TO_CONSTANT[pattern]
+
+    def get_directivity_pattern(self):
+        return self._pattern
+
+    def get_azimuth(self, degrees=True):
+        return self._orientation.get_azimuth(degrees)
+
+    def get_colatitude(self, degrees=True):
+        return self._orientation.get_colatitude(degrees)
+
+    def get_response(self, direction):
+        """
+        Get response for a single direction.
+
+        TODO : vectorize
+
+        Parameters
+        ----------
+        direction : DirectionVector
+            Direction for which to compute gain
+
+        """
+        assert isinstance(direction, DirectionVector)
+        if self._pattern == DirectivityPattern.OMNI:
+            return 1
+        else:
+            # inspiration from Habets: https://github.com/ehabets/RIR-Generator/blob/5eb70f066b74ff18c2be61c97e8e666f8492c149/rir_generator.cpp#L111
+            # we use colatitude instead of elevation
+            gain = np.sin(self._orientation.get_colatitude()) * \
+                   np.sin(direction.get_colatitude()) * \
+                   np.cos(
+                       self._orientation.get_azimuth() -
+                       direction.get_azimuth()
+                   ) + \
+                   np.cos(self._orientation.get_colatitude()) * \
+                   np.cos(direction.get_colatitude())
+            return np.abs(self._p + (1 - self._p) * gain)
+
+    def plot_response(self, azimuth, colatitude=None, degrees=True):
+        """
+        Parameters
+        ----------
+        azimuth : array_like
+            Azimuth values for plotting.
+        colatitude : array_like, optional
+            Colatitude values for plotting. If not provided, 2D plot.
+        degrees : bool
+            Whether provided values are in degrees (True) or radians (False).
+        """
+
+        try:
+            import matplotlib.pyplot as plt
+        except ImportError:
+            import warnings
+            warnings.warn('Matplotlib is required for plotting')
+            return
+
+        fig = plt.figure()
+        if colatitude is not None:
+
+            # compute response
+            gains = np.zeros(shape=(len(azimuth), len(colatitude)))
+            for i, a in enumerate(azimuth):
+                for j, c in enumerate(colatitude):
+                    direction = DirectionVector(
+                        azimuth=a,
+                        colatitude=c,
+                        degrees=degrees
+                    )
+                    gains[i, j] = self.get_response(direction)
+
+            # create surface plot, need cartesian coordinates
+            if degrees:
+                azimuth = np.radians(azimuth)
+                colatitude = np.radians(colatitude)
+            AZI, COL = np.meshgrid(azimuth, colatitude)
+
+            X = gains.T * np.sin(COL) * np.cos(AZI)
+            Y = gains.T * np.sin(COL) * np.sin(AZI)
+            Z = gains.T * np.cos(COL)
+
+            ax = fig.add_subplot(1, 1, 1, projection='3d')
+            ax.plot_surface(X, Y, Z)
+            ax.set_title("{}, azimuth={}, colatitude={}".format(
+                self.get_directivity_pattern(),
+                self.get_azimuth(),
+                self.get_colatitude()
+            ))
+            ax.set_xlabel("x")
+            ax.set_ylabel("y")
+            ax.set_zlabel("z")
+            ax.set_xlim([-1, 1])
+            ax.set_ylim([-1, 1])
+            ax.set_zlim([-1, 1])
+
+        else:
+
+            # compute response
+            gains = np.zeros(len(azimuth))
+            for i, a in enumerate(azimuth):
+                direction = DirectionVector(azimuth=a, degrees=degrees)
+                gains[i] = self.get_response(direction)
+
+            # plot
+            ax = plt.subplot(111, projection="polar")
+            if degrees:
+                angles = np.radians(azimuth)
+            else:
+                angles = azimuth
+            ax.plot(angles, gains)

pyroomacoustics/random/__init__.py

@@ -0,0 +1,2 @@
+from . import distributions, sampler
+from .spherical import power_spherical, uniform_spherical

pyroomacoustics/random/distributions.py

@@ -0,0 +1,116 @@
+import abc
+
+import numpy as np
+import scipy
+
+from .spherical import power_spherical, uniform_spherical
+
+
+class Distribution(abc.ABC):
+    """
+    Abstract base class for distributions. Derived classes
+    should implement a ``pdf`` method that provides the pdf value for samples
+    and a ``sample`` method that allows to sample from the distribution
+    """
+
+    def __init__(self, dim):
+        self._dim = dim
+
+    @property
+    def dim(self):
+        return self._dim
+
+    @abc.abstractmethod
+    def pdf(self):
+        pass
+
+    @abc.abstractmethod
+    def sample(self):
+        pass
+
+    def __call__(self, *args, **kwargs):
+        return self.sample(*args, **kwargs)
+
+
+class AdHoc(Distribution):
+    """
+    A helper class to construct distribution from separate pdf/sample
+    functions
+    """
+
+    def __init__(self, pdf, sampler, dim):
+        super().__init__(dim)
+
+        self._pdf = pdf
+        self._sampler = sampler
+
+    def pdf(self, *args, **kwargs):
+        return self._pdf(*args, **kwargs)
+
+    def sample(self, *args, **kwargs):
+        return self._sampler(*args, **kwargs)
+
+
+class UniformSpherical(Distribution):
+    """
+    Uniform distribution on the n-sphere
+    """
+
+    def __init__(self, dim=3):
+        super().__init__(dim)
+        self._area = (
+            self._dim
+            * np.pi ** (self._dim / 2)
+            / scipy.special.gamma(self._dim / 2 + 1)
+        )
+
+    @property
+    def _n_sphere_area(self):
+        return self._area
+
+    def pdf(self, x):
+        return 1.0 / self._n_sphere_area
+
+    def sample(self, size=None):
+        return uniform_spherical(dim=self.dim, size=size)
+
+
+class PowerSpherical(Distribution):
+    def __init__(self, dim=3, loc=None, scale=None):
+        super().__init__(dim)
+
+        # first canonical basis vector
+        if loc is None:
+            self._loc = np.zeros(self.dim)
+            self._loc[0] = 1.0
+        else:
+            self._loc = loc
+
+        if not scale:
+            self._scale = 1.0
+        else:
+            self._scale = scale
+
+        # computation the normalizing constant
+        a = (self.dim - 1) / 2 + self.scale
+        b = (self.dim - 1) / 2
+        self._Z_inv = 1.0 / (
+            2 ** (a + b)
+            * np.pi ** b
+            * scipy.special.gamma(a)
+            / scipy.special.gamma(a + b)
+        )
+
+    @property
+    def loc(self):
+        return self._loc
+
+    @property
+    def scale(self):
+        return self._scale
+
+    def pdf(self, x):
+        assert (
+            x.shape[-1] == self.dim
+        ), "Input dimension does not match distribution dimension"
+        return self._Z_inv * (1.0 + np.matmul(x, self.loc)) ** self.scale