Skip to content

Commit

Permalink
Switches to using source directions computed as part of the image sou…
Browse files Browse the repository at this point in the history
…rce model for source directivities. This enables the use of source directivities with non-shoebox rooms.
  • Loading branch information
fakufaku committed Dec 31, 2024
1 parent 3e177a7 commit cf35908
Show file tree
Hide file tree
Showing 5 changed files with 142 additions and 87 deletions.
24 changes: 7 additions & 17 deletions examples/directivities/source_and_microphone.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,21 +2,15 @@
import numpy as np

import pyroomacoustics as pra
from pyroomacoustics.directivities import (
CardioidFamily,
DirectionVector,
DirectivityPattern,
)
from pyroomacoustics.directivities import DirectionVector, FigureEight, HyperCardioid

three_dim = True # 2D or 3D
shoebox = True # source directivity not supported for non-shoebox!
energy_absorption = 0.4
source_pos = [2, 1.8]
# source_dir = None # to disable
source_dir = DirectivityPattern.FIGURE_EIGHT
mic_pos = [3.5, 1.8]
# mic_dir = None # to disable
mic_dir = DirectivityPattern.HYPERCARDIOID


# make 2-D room
Expand All @@ -42,19 +36,15 @@
colatitude = None

# add source with directivity
if source_dir is not None:
source_dir = CardioidFamily(
orientation=DirectionVector(azimuth=90, colatitude=colatitude, degrees=True),
pattern_enum=source_dir,
)
source_dir = FigureEight(
orientation=DirectionVector(azimuth=90, colatitude=colatitude, degrees=True),
)
room.add_source(position=source_pos, directivity=source_dir)

# add microphone with directivity
if mic_dir is not None:
mic_dir = CardioidFamily(
orientation=DirectionVector(azimuth=0, colatitude=colatitude, degrees=True),
pattern_enum=mic_dir,
)
mic_dir = HyperCardioid(
orientation=DirectionVector(azimuth=0, colatitude=colatitude, degrees=True),
)
room.add_microphone(loc=mic_pos, directivity=mic_dir)

# plot room
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
import numpy as np
import pytest

import pyroomacoustics as pra

"""
+ - - - - - - - - - - - +
| |
| |
| |
| + - - - - - - - +
| |
| |
| |
+ - - - +
"""

corners = np.array([[0, 0], [4, 0], [4, 4], [12, 4], [12, 8], [0, 8]])
source_location = [6, 6, 1.5]
mic_location = [1, 6, 1.5]

images_expected = np.array(
[
source_location, # direct
[6, 10, 1.5], # north
[-6, 6, 1.5], # west
[18, 6, 1.5], # east
[6, 6, 4.5], # ceilling
[6, 6, -1.5], # floor
]
)
directions_expected = np.array(
[
[-1, 0, 0],
np.array([-2.5, 2, 0]) / np.sqrt(2.5**2 + 2**2),
[-1, 0, 0],
[1, 0, 0],
np.array([-2.5, 0, 1.5]) / np.sqrt(2.5**2 + 1.5**2),
np.array([-2.5, 0, -1.5]) / np.sqrt(2.5**2 + 1.5**2),
]
)


room = pra.Room.from_corners(corners.T, materials=pra.Material(0.1), max_order=1)
room.extrude(3.0, materials=pra.Material(0.1))
room.add_source(
source_location,
directivity=pra.directivities.Cardioid(
orientation=pra.directivities.DirectionVector(
azimuth=0, colatitude=np.pi / 2, degrees=False
)
),
).add_microphone(mic_location).add_microphone([3, 6, 1.5])
# We add one microphone with other images visible to make sure we have some
# image sources not visible.


def test_source_directions_nonshoebox():
room.image_source_model()
visible = room.visibility[0][0, :]
directions_obtained = room.sources[0].directions[0, :, visible]
images_obtained = room.sources[0].images[:, visible]

order_obtained = np.lexsort(images_obtained)
images_obtained = images_obtained[:, order_obtained].T
directions_obtained = directions_obtained[order_obtained, :]

order_expected = np.lexsort(images_expected.T)
images_expected_reordered = images_expected[order_expected, :]
directions_expected_reordered = directions_expected[order_expected, :]

np.testing.assert_almost_equal(images_obtained, images_expected_reordered)
np.testing.assert_almost_equal(directions_obtained, directions_expected_reordered)
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,57 @@
import numpy as np

import pyroomacoustics as pra
from pyroomacoustics.simulation.ism import source_angle_shoebox


def source_angle_shoebox(image_source_loc, wall_flips, mic_loc):
"""
Determine outgoing angle for each image source for a ShoeBox configuration.
Implementation of the method described in the paper:
https://www2.ak.tu-berlin.de/~akgroup/ak_pub/2018/000458.pdf
Parameters
-----------
image_source_loc : array_like
Locations of image sources.
wall_flips: array_like
Number of x, y, z flips for each image source.
mic_loc: array_like
Microphone location.
Returns
-------
azimuth : :py:class:`~numpy.ndarray`
Azimith for each image source, in radians
colatitude : :py:class:`~numpy.ndarray`
Colatitude for each image source, in radians.
"""

image_source_loc = np.array(image_source_loc)
wall_flips = np.array(wall_flips)
mic_loc = np.array(mic_loc)

dim, n_sources = image_source_loc.shape
assert wall_flips.shape[0] == dim
assert mic_loc.shape[0] == dim

p_vector_array = image_source_loc - np.array(mic_loc)[:, np.newaxis]
d_array = np.linalg.norm(p_vector_array, axis=0)

# Using (12) from the paper
power_array = np.ones_like(image_source_loc) * -1
power_array = np.power(power_array, (wall_flips + np.ones_like(image_source_loc)))
p_dash_array = p_vector_array * power_array

# Using (13) from the paper
azimuth = np.arctan2(p_dash_array[1], p_dash_array[0])
if dim == 2:
colatitude = np.ones(n_sources) * np.pi / 2
else:
colatitude = np.pi / 2 - np.arcsin(p_dash_array[2] / d_array)

return azimuth, colatitude


class TestSourceDirectivityFlipping(TestCase):
Expand Down
20 changes: 8 additions & 12 deletions pyroomacoustics/room.py
Original file line number Diff line number Diff line change
Expand Up @@ -1409,8 +1409,8 @@ def extrude(self, height, v_vec=None, absorption=None, materials=None):
if len(self.walls[i].absorption) == 1:
# Single band
wall_materials[name] = Material(
energy_absorption=float(self.walls[i].absorption),
scattering=float(self.walls[i].scatter),
energy_absorption=float(self.walls[i].absorption[0]),
scattering=float(self.walls[i].scatter[0]),
)
elif len(self.walls[i].absorption) == self.octave_bands.n_bands:
# Multi-band
Expand Down Expand Up @@ -2139,14 +2139,6 @@ def add_source(self, position, signal=None, delay=0, directivity=None):
if self.dim != 3 and directivity is not None:
raise NotImplementedError("Directivity is only supported for 3D rooms.")

if directivity is not None:
from pyroomacoustics import ShoeBox

if not isinstance(self, ShoeBox):
raise NotImplementedError(
"Source directivity only supported for ShoeBox room."
)

if isinstance(position, SoundSource):
if directivity is not None:
if isinstance(directivity, CardioidFamily) or isinstance(
Expand Down Expand Up @@ -2241,8 +2233,11 @@ def image_source_model(self):
self.visibility[-1][m, :] = 0
else:
# if we are here, this means even the direct path is not visible
# we set the visibility of the direct path as 0
# we set the visibility of the direct path as 0.
self.visibility.append(np.zeros((self.mic_array.M, 1), dtype=np.int32))
# We also need a fake array of directions as this is expected later in
# the code.
source.directions = np.zeros((self.mic_array.M, self.dim, 1), dtype=np.float32)

# Update the state
self.simulator_state["ism_done"] = True
Expand Down Expand Up @@ -2305,13 +2300,14 @@ def compute_rir(self):
src,
mic,
self.mic_array.directivity[m],
src.directions[m, :, :],
self.visibility[s][m, :],
fdl,
self.c,
self.fs,
self.octave_bands,
air_abs_coeffs=self.air_absorption,
min_phase=self.min_phase,
air_abs_coeffs=self.air_absorption,
)
rir_parts.append(ir_ism)

Expand Down
60 changes: 3 additions & 57 deletions pyroomacoustics/simulation/ism.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
import numpy as np
from scipy.signal import fftconvolve, hilbert

from .. import libroom
from .. import doa, libroom
from ..parameters import constants
from ..utilities import angle_function

Expand Down Expand Up @@ -109,61 +109,11 @@ def interpolate_octave_bands(
return ir


def source_angle_shoebox(image_source_loc, wall_flips, mic_loc):
"""
Determine outgoing angle for each image source for a ShoeBox configuration.
Implementation of the method described in the paper:
https://www2.ak.tu-berlin.de/~akgroup/ak_pub/2018/000458.pdf
Parameters
-----------
image_source_loc : array_like
Locations of image sources.
wall_flips: array_like
Number of x, y, z flips for each image source.
mic_loc: array_like
Microphone location.
Returns
-------
azimuth : :py:class:`~numpy.ndarray`
Azimith for each image source, in radians
colatitude : :py:class:`~numpy.ndarray`
Colatitude for each image source, in radians.
"""

image_source_loc = np.array(image_source_loc)
wall_flips = np.array(wall_flips)
mic_loc = np.array(mic_loc)

dim, n_sources = image_source_loc.shape
assert wall_flips.shape[0] == dim
assert mic_loc.shape[0] == dim

p_vector_array = image_source_loc - np.array(mic_loc)[:, np.newaxis]
d_array = np.linalg.norm(p_vector_array, axis=0)

# Using (12) from the paper
power_array = np.ones_like(image_source_loc) * -1
power_array = np.power(power_array, (wall_flips + np.ones_like(image_source_loc)))
p_dash_array = p_vector_array * power_array

# Using (13) from the paper
azimuth = np.arctan2(p_dash_array[1], p_dash_array[0])
if dim == 2:
colatitude = np.ones(n_sources) * np.pi / 2
else:
colatitude = np.pi / 2 - np.arcsin(p_dash_array[2] / d_array)

return azimuth, colatitude


def compute_ism_rir(
src,
mic,
mic_dir,
src_directions,
is_visible,
fdl,
c,
Expand Down Expand Up @@ -226,11 +176,7 @@ def compute_ism_rir(
oct_band_amplitude *= mic_gain

if src.directivity is not None:
azimuth_s, colatitude_s = source_angle_shoebox(
image_source_loc=images,
wall_flips=abs(src.orders_xyz[:, is_visible]),
mic_loc=mic,
)
azimuth_s, colatitude_s, _ = doa.cart2spher(src_directions[:, is_visible])
src_gain = src.directivity.get_response(
azimuth=azimuth_s,
colatitude=colatitude_s,
Expand Down

0 comments on commit cf35908

Please sign in to comment.