diff --git a/.github/workflows/pythonpackage.yml b/.github/workflows/pythonpackage.yml index b8bbc7f7..3572911a 100644 --- a/.github/workflows/pythonpackage.yml +++ b/.github/workflows/pythonpackage.yml @@ -35,9 +35,17 @@ jobs: uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} - - name: Install dependencies + - name: Upgrade pip run: | python -m pip install --upgrade pip + - name: Workaround for windows and python 3.8 + if: matrix.os == 'windows-2019' && matrix.python-version == 3.8 + run: | + pip install netCDF4<=1.5.8 + # There is no binary package of netCF4>=1.6.0 for windows and python 3.7 + # the largest supported version is 1.5.8 + - name: Install dependencies + run: | pip install -r requirements.txt - name: Build package run: | diff --git a/.gitignore b/.gitignore index 9ad4abe1..bd81b87b 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,4 @@ local_examples/ pyroomacoustics.egg-info/ pyroomacoustics/build_rir.c pyroomacoustics/**/*.so +pyroomacoustics/directivities/tests/data/*.pdf diff --git a/CHANGELOG.rst b/CHANGELOG.rst index dde500a8..d33499ce 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -11,7 +11,55 @@ adheres to `Semantic Versioning `_. `Unreleased`_ ------------- -Nothing yet +This new version introduces some major changes. In particular, it introduces +the use of measured microphone and source directivities. This new feature +allows to simulate more accurately real recording equipments. See the +`documentation +`_ +for more details. + +Added +~~~~~ + +- New global parameters to control the octave bands used for simulation. + + - ``octave_bands_base_freq``: the base frequency used for the octave bands (default ``125``), + note that together with the sampling frequency this will determine the number of sub-bands + used in simulation. + - ``octave_bands_n_fft``: lengths of the octave band filters (current is default ``512`` + but will be changed to ``128`` in the next release) + - ``octave_bands_keep_dc``: extends the lowest band to include the DC offset, + the current default is ``False`` to match past behavior, but will be changed to + ``True`` in the next release because the filters have less oscillations this way. + +- New parameter ``min_phase`` of class ``Room``. If set to ``True``, the band-pass + filters are designed as minimum phase filters. + +- Simulation with measured directivity responses in SOFA format (limited file types) is + possible with the image source model. + + - A flexible and extensible SOFA file reader with optional spherical interpolation. + - A small database of files is bundled, including the `DIRPAT database + `_, collected by + Manuel Brandner, Matthias Frank, and Daniel Rudrich University of Music and + Performing Arts Graz, Graz. The file also include two HRTF of the + `MIT KEMAR dummy head `_. + +- Introduces a new class ``pyroomacoustics.directivities.Rotation3D`` to specify + the orientation of sources and receivers in 3D. + +- Adds `soxr `_ as a dependency since resampling + can often be necessary when working with SOFA files. + +Changed +~~~~~~~ + +- In ray tracing, the histograms are now linearly interpolated between + the bins to obtain smoother RIR +- Changed the API of ``CardioidFamily`` to take a float parameter. + New class wrappers for ``Cardioid``, ``Hypercardioid``, ``Supercardioid``, + ``Bidirectional`` and ``Omnidirectional`` are added in the ``directivity`` + sub-module. The enum of type ``DirectivityPattern`` has been removed. `0.7.7`_ - 2024-09-09 --------------------- diff --git a/MANIFEST.in b/MANIFEST.in index b5c0eb8d..0b82c592 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,5 +1,10 @@ # Acoustic data for the simulation include pyroomacoustics/data/materials.json +include pyroomacoustics/data/sofa_files.json +include pyroomacoustics/data/sofa/AKG_c480_c414_CUBE.sofa +include pyroomacoustics/data/sofa/EM32_Directivity.sofa +include pyroomacoustics/data/sofa/mit_kemar_large_pinna.sofa +include pyroomacoustics/data/sofa/mit_kemar_normal_pinna.sofa # The header files for the compiled extension include pyroomacoustics/libroom_src/*.h @@ -19,9 +24,12 @@ graft pyroomacoustics/tests graft pyroomacoustics/adaptive/tests graft pyroomacoustics/bss/tests graft pyroomacoustics/datasets/tests +graft pyroomacoustics/denoise/tests +graft pyroomacoustics/directivities/tests graft pyroomacoustics/doa/tests graft pyroomacoustics/experimental/tests graft pyroomacoustics/libroom_src/tests +graft pyroomacoustics/phase/tests graft pyroomacoustics/transform/tests global-exclude *.py[co] diff --git a/docs/conf.py b/docs/conf.py index 81ddbe94..3d84eb91 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -52,6 +52,7 @@ "mpl_toolkits.mplot3d", "joblib", "scipy.io", + "soxr", ] try: diff --git a/docs/index.rst b/docs/index.rst index b90bd82a..74d99231 100755 --- a/docs/index.rst +++ b/docs/index.rst @@ -14,6 +14,7 @@ Table of contents contributing changelog pyroomacoustics.room + pyroomacoustics.directivities pyroomacoustics.materials.database pyroomacoustics.transform pyroomacoustics.datasets diff --git a/docs/modules.rst b/docs/modules.rst deleted file mode 100644 index ece765c7..00000000 --- a/docs/modules.rst +++ /dev/null @@ -1,7 +0,0 @@ -pyroomacoustics -=============== - -.. toctree:: - :maxdepth: 4 - - pyroomacoustics diff --git a/docs/pyroomacoustics.datasets.rst b/docs/pyroomacoustics.datasets.rst index 6c463630..6ffb410c 100644 --- a/docs/pyroomacoustics.datasets.rst +++ b/docs/pyroomacoustics.datasets.rst @@ -16,6 +16,7 @@ Datasets Available pyroomacoustics.datasets.cmu_arctic pyroomacoustics.datasets.google_speech_commands + pyroomacoustics.datasets.sofa pyroomacoustics.datasets.timit Tools and Helpers diff --git a/docs/pyroomacoustics.datasets.sofa.rst b/docs/pyroomacoustics.datasets.sofa.rst new file mode 100644 index 00000000..33cc2f13 --- /dev/null +++ b/docs/pyroomacoustics.datasets.sofa.rst @@ -0,0 +1,5 @@ +SOFA Database +============= + +.. automodule:: pyroomacoustics.datasets.sofa + :members: diff --git a/docs/pyroomacoustics.directivities.rst b/docs/pyroomacoustics.directivities.rst index 37dcdead..a64f40a2 100644 --- a/docs/pyroomacoustics.directivities.rst +++ b/docs/pyroomacoustics.directivities.rst @@ -1,7 +1,75 @@ -pyroomacoustics.directivities module -==================================== +Directional Sources and Microphones +=================================== .. automodule:: pyroomacoustics.directivities :members: :undoc-members: :show-inheritance: + +Analytic Directional Responses +------------------------------ + +.. automodule:: pyroomacoustics.directivities.analytic + :members: + :undoc-members: + :show-inheritance: + +Measured Directivities +---------------------- + +.. automodule:: pyroomacoustics.directivities.measured + :members: + :undoc-members: + :show-inheritance: + +Built-in SOFA Files Database +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. automodule:: pyroomacoustics.datasets.sofa + :members: SOFADatabase, get_sofa_db + :noindex: + + +Reading Other or Custom File Types +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +It is possible to read other file types by providing a custom reader function to +:py:class:`~pyroomacoustics.directivities.measured.MeasuredDirectivityFile` with the +argument ``file_reader_callback``. +The function should have the same signature as :py:func:`~pyroomacoustics.directivities.sofa.open_sofa_file`. + + +SOFA File Readers +................. + +.. automodule:: pyroomacoustics.directivities.sofa + :members: + :show-inheritance: + +Direction of the Patterns +------------------------- + +.. automodule:: pyroomacoustics.directivities.direction + :members: + :show-inheritance: + + +Creating New Types of Directivities +----------------------------------- + +.. automodule:: pyroomacoustics.directivities.base + :members: + :undoc-members: + :show-inheritance: + +Spherical Interpolation +----------------------- + +.. automodule:: pyroomacoustics.directivities.interp + :members: spherical_interpolation + +Numerical Spherical Integral +---------------------------- + +.. automodule:: pyroomacoustics.directivities.integration + :members: spherical_integral diff --git a/docs/pyroomacoustics.libroom.rst b/docs/pyroomacoustics.libroom.rst deleted file mode 100644 index 4591d24c..00000000 --- a/docs/pyroomacoustics.libroom.rst +++ /dev/null @@ -1,7 +0,0 @@ -pyroomacoustics.libroom module -============================== - -.. automodule:: pyroomacoustics.libroom - :members: - :undoc-members: - :show-inheritance: diff --git a/docs/pyroomacoustics.rst b/docs/pyroomacoustics.rst index 29deac40..8a61fa34 100644 --- a/docs/pyroomacoustics.rst +++ b/docs/pyroomacoustics.rst @@ -29,6 +29,7 @@ Submodules pyroomacoustics.parameters pyroomacoustics.recognition pyroomacoustics.room + pyroomacoustics.simulation pyroomacoustics.soundsource pyroomacoustics.stft pyroomacoustics.sync diff --git a/docs/pyroomacoustics.simulation.ism.rst b/docs/pyroomacoustics.simulation.ism.rst new file mode 100644 index 00000000..93672000 --- /dev/null +++ b/docs/pyroomacoustics.simulation.ism.rst @@ -0,0 +1,7 @@ +pyroomacoustics.simulation.ism module +===================================== + +.. automodule:: pyroomacoustics.simulation.ism + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/pyroomacoustics.simulation.rst b/docs/pyroomacoustics.simulation.rst new file mode 100644 index 00000000..2bb1cf95 --- /dev/null +++ b/docs/pyroomacoustics.simulation.rst @@ -0,0 +1,21 @@ +Simulation Routines Sub-package +=============================== + +This sub-package contains some internal routines to simulate room acoustics. + +Submodules +---------- + +.. toctree:: + :maxdepth: 4 + + pyroomacoustics.simulation.ism + pyroomacoustics.simulation.rt + +Module contents +--------------- + +.. automodule:: pyroomacoustics.simulation + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/pyroomacoustics.simulation.rt.rst b/docs/pyroomacoustics.simulation.rt.rst new file mode 100644 index 00000000..d7b55f7a --- /dev/null +++ b/docs/pyroomacoustics.simulation.rt.rst @@ -0,0 +1,7 @@ +pyroomacoustics.simulation.rt module +==================================== + +.. automodule:: pyroomacoustics.simulation.rt + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/pyroomacoustics.version.rst b/docs/pyroomacoustics.version.rst deleted file mode 100644 index e71268d2..00000000 --- a/docs/pyroomacoustics.version.rst +++ /dev/null @@ -1,7 +0,0 @@ -pyroomacoustics\.version module -=============================== - -.. automodule:: pyroomacoustics.version - :members: - :undoc-members: - :show-inheritance: diff --git a/examples/directivities/cardioid_function.py b/examples/directivities/cardioid_function.py index 948acbf2..408ee2be 100644 --- a/examples/directivities/cardioid_function.py +++ b/examples/directivities/cardioid_function.py @@ -15,7 +15,7 @@ direction = spher2cart(azimuth=225, degrees=True) # compute response -resp = cardioid_func(x=cart, direction=direction, coef=0.5, magnitude=True) +resp = cardioid_func(x=cart, direction=direction, p=0.5, magnitude=True) resp_db = dB(np.array(resp)) # plot @@ -33,7 +33,7 @@ direction = spher2cart(azimuth=0, colatitude=45, degrees=True) # compute response -resp = cardioid_func(x=cart, direction=direction, coef=0.25, magnitude=True) +resp = cardioid_func(x=cart, direction=direction, p=0.25, magnitude=True) # plot (surface plot) fig = plt.figure() diff --git a/examples/directivities/circular_mic_array.py b/examples/directivities/circular_mic_array.py index d9611048..aac53637 100644 --- a/examples/directivities/circular_mic_array.py +++ b/examples/directivities/circular_mic_array.py @@ -1,11 +1,7 @@ import matplotlib.pyplot as plt import pyroomacoustics as pra -from pyroomacoustics.directivities import ( - CardioidFamily, - DirectionVector, - DirectivityPattern, -) +from pyroomacoustics.directivities import Cardioid, DirectionVector three_dim = True # 2D or 3D @@ -28,9 +24,8 @@ room.add_source(source_loc) # add circular microphone array -pattern = DirectivityPattern.CARDIOID orientation = DirectionVector(azimuth=mic_rotation, colatitude=colatitude, degrees=True) -directivity = CardioidFamily(orientation=orientation, pattern_enum=pattern) +directivity = Cardioid(orientation=orientation) mic_array = pra.beamforming.circular_microphone_array_xyplane( center=center, M=7, diff --git a/examples/directivities/mic_array.py b/examples/directivities/mic_array.py index e1c8dd20..bfbc782b 100644 --- a/examples/directivities/mic_array.py +++ b/examples/directivities/mic_array.py @@ -2,13 +2,8 @@ import numpy as np import pyroomacoustics as pra -from pyroomacoustics.directivities import ( - CardioidFamily, - DirectionVector, - DirectivityPattern, -) +from pyroomacoustics.directivities import DirectionVector, HyperCardioid -pattern = DirectivityPattern.HYPERCARDIOID orientation = DirectionVector(azimuth=0, colatitude=0, degrees=True) # create room @@ -26,7 +21,7 @@ M = 3 R = pra.linear_2D_array(center=[5, 5], M=M, phi=0, d=0.7) R = np.concatenate((R, np.ones((1, M)))) -directivity = CardioidFamily(orientation=orientation, pattern_enum=pattern) +directivity = HyperCardioid(orientation=orientation) room.add_microphone_array(R, directivity=directivity) # plot room diff --git a/examples/directivities/mic_array_diff_directivities.py b/examples/directivities/mic_array_diff_directivities.py index adf5d057..953a64a4 100644 --- a/examples/directivities/mic_array_diff_directivities.py +++ b/examples/directivities/mic_array_diff_directivities.py @@ -2,24 +2,17 @@ import numpy as np import pyroomacoustics as pra -from pyroomacoustics.directivities import ( - CardioidFamily, - DirectionVector, - DirectivityPattern, -) +from pyroomacoustics.directivities import Cardioid, DirectionVector, HyperCardioid -dir_1 = CardioidFamily( +dir_1 = Cardioid( orientation=DirectionVector(azimuth=180, colatitude=30, degrees=True), - pattern_enum=DirectivityPattern.HYPERCARDIOID, ) -dir_2 = CardioidFamily( +dir_2 = HyperCardioid( orientation=DirectionVector(azimuth=0, colatitude=30, degrees=True), - pattern_enum=DirectivityPattern.HYPERCARDIOID, ) # source_dir = None -source_dir = CardioidFamily( +source_dir = HyperCardioid( orientation=DirectionVector(azimuth=45, colatitude=90, degrees=True), - pattern_enum=DirectivityPattern.CARDIOID, ) # create room diff --git a/examples/directivities/plot_directivity_2D.py b/examples/directivities/plot_directivity_2D.py index 5019ff24..9b5d39ec 100644 --- a/examples/directivities/plot_directivity_2D.py +++ b/examples/directivities/plot_directivity_2D.py @@ -3,9 +3,11 @@ from pyroomacoustics import dB from pyroomacoustics.directivities import ( - CardioidFamily, + Cardioid, DirectionVector, - DirectivityPattern, + FigureEight, + HyperCardioid, + SubCardioid, ) orientation = DirectionVector(azimuth=0, colatitude=90, degrees=True) @@ -18,11 +20,11 @@ # plot each pattern fig = plt.figure() ax = plt.subplot(111, projection="polar") -for pattern in DirectivityPattern: - dir_obj = CardioidFamily(orientation=orientation, pattern_enum=pattern) +for obj in [SubCardioid, Cardioid, HyperCardioid, FigureEight]: + dir_obj = obj(orientation=orientation) resp = dir_obj.get_response(azimuth=angles, magnitude=True, degrees=False) resp_db = dB(np.array(resp)) - ax.plot(angles, resp_db, label=pattern.name) + ax.plot(angles, resp_db, label=dir_obj.directivity_pattern) plt.legend(bbox_to_anchor=(1, 1)) plt.ylim([lower_gain, 0]) diff --git a/examples/directivities/plot_directivity_3D.py b/examples/directivities/plot_directivity_3D.py index dc47e7b3..bafa3764 100644 --- a/examples/directivities/plot_directivity_3D.py +++ b/examples/directivities/plot_directivity_3D.py @@ -1,21 +1,16 @@ import matplotlib.pyplot as plt import numpy as np -from pyroomacoustics.directivities import ( - CardioidFamily, - DirectionVector, - DirectivityPattern, -) +from pyroomacoustics.directivities import DirectionVector, HyperCardioid -pattern = DirectivityPattern.HYPERCARDIOID orientation = DirectionVector(azimuth=0, colatitude=45, degrees=True) # create cardioid object -dir_obj = CardioidFamily(orientation=orientation, pattern_enum=pattern) +dir_obj = HyperCardioid(orientation=orientation) # plot azimuth = np.linspace(start=0, stop=360, num=361, endpoint=True) -colatitude = np.linspace(start=0, stop=180, num=180, endpoint=True) +colatitude = np.linspace(start=0, stop=180, num=181, endpoint=True) # colatitude = None # for 2D plot dir_obj.plot_response(azimuth=azimuth, colatitude=colatitude, degrees=True) plt.show() diff --git a/examples/directivities/simulate_binaural_recording.py b/examples/directivities/simulate_binaural_recording.py new file mode 100644 index 00000000..9b65e1cc --- /dev/null +++ b/examples/directivities/simulate_binaural_recording.py @@ -0,0 +1,164 @@ +# 2023 Robin Scheibler +""" +This scripts simulates a binaural room impulse response + +The directivity patterns are loaded from SOFA files. +The SOFA format is described at https://www.sofaconventions.org + +Get more SOFA files at https://www.sofaconventions.org/mediawiki/index.php/Files +""" +import argparse +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from scipy import signal +from scipy.fft import fft, fftfreq +from scipy.io import wavfile +from scipy.signal import fftconvolve + +import pyroomacoustics as pra +from pyroomacoustics.directivities import MeasuredDirectivityFile, Rotation3D + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Simulate a binaural room impulse response" + ) + parser.add_argument( + "--hrtf", + type=Path, + default="mit_kemar_normal_pinna.sofa", + help="Path to HRTF file", + ) + parser.add_argument( + "--source", + type=Path, + default=Path(__file__).parents[1] / "input_samples/cmu_arctic_us_axb_a0004.wav", + help="Path to speech or other source audio file", + ) + parser.add_argument( + "--output", + type=Path, + default=Path(__file__).parents[1] + / "output_samples/simulate_binaural_recording.wav", + help="Path to output file", + ) + parser.add_argument( + "--interp-order", + type=int, + default=12, + help="Maximum order to use in the spherical harmonics interpolation. " + "Setting to -1 will disable interpolation", + ) + parser.add_argument( + "--interp-n-points", + type=int, + default=1000, + help="Maximum order to use in the spherical harmonics interpolation", + ) + args = parser.parse_args() + + interp_order = 12 + interp_n_points = 1000 + + fs, speech = wavfile.read(args.source) + speech = speech * (0.95 / abs(speech).max()) + azimuth_deg = -45.0 + colatitude_deg = 90.0 + + hrtf = MeasuredDirectivityFile( + path=args.hrtf, + fs=fs, + interp_order=args.interp_order, + interp_n_points=args.interp_n_points, + ) + orientation = Rotation3D([colatitude_deg, azimuth_deg], "yz", degrees=True) + dir_left = hrtf.get_mic_directivity("left", orientation=orientation) + dir_right = hrtf.get_mic_directivity("right", orientation=orientation) + + room_dim = [6, 6, 2.4] + + all_materials = { + "east": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "west": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "north": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "south": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "ceiling": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "floor": pra.Material( + energy_absorption={ + "coeffs": [0.11, 0.14, 0.37, 0.43, 0.27, 0.25], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + } + + room = pra.ShoeBox( + room_dim, + fs=fs, + max_order=40, + # materials=pra.Material(0.5), + materials=all_materials, + air_absorption=True, + ray_tracing=False, + min_phase=False, + use_rand_ism=True, + max_rand_disp=0.05, + ) + + room.add_source([1.5, 3.01, 1.044], signal=speech) + + room.add_microphone([1.1, 3.01, 2.2], directivity=dir_left) + room.add_microphone([1.1, 3.01, 2.2], directivity=dir_right) + + room.simulate() + + signals = room.mic_array.signals + signals *= 0.95 / abs(signals).max() + signals = (signals * 2**15).astype(np.int16) + wavfile.write(args.output, fs, signals.T) + + room.plot_rir(FD=True) + room.plot_rir(FD=False) + + fig = plt.figure() + for idx, fb in enumerate(range(44)): + if idx >= 5 * 10: + break + ax = fig.add_subplot(5, 10, idx + 1, projection="3d") + dir_left.plot(freq_bin=fb, ax=ax, depth=True) + # ax.set_xticks([]) + # ax.set_yticks([]) + # ax.set_zticks([]) + ax.set_title(idx) + plt.show() diff --git a/examples/directivities/simulation_with_measured_directivity.py b/examples/directivities/simulation_with_measured_directivity.py new file mode 100644 index 00000000..f4291666 --- /dev/null +++ b/examples/directivities/simulation_with_measured_directivity.py @@ -0,0 +1,169 @@ +# 2022 (c) Prerak SRIVASTAVA +# 2024/11/27 Modified by Robin Scheibler +""" +Simulating RIRs with measured directivity patterns from DIRPAT +============================================================== + +In this example, we show how we can use measured directivity patterns +from the DIRPAT dataset in a simulation. + +The procedure to use the directivity patterns is as follows. + +1. Read the files potentially containing multiple measurements. +2. Get a directivity object from the file object with desired orientation. + The directivities can be accessed by index or label (if existing). + The same pattern can be used multiple times with different orientations. +3. Provide the directivity pattern object to the microphone object. + +The DIRPAT database has three different files. + +The ``AKG_c480_c414_CUBE.sofa`` DIRPAT file include mic patterns for CARDIOID, +FIGURE_EIGHT, HYPERCARDIOID, OMNI, SUBCARDIOID. + +a)AKG_c480 +b)AKG_c414K +c)AKG_c414N +d)AKG_c414S +e)AKG_c414A + +The Eigenmic directivity pattern file ``EM32_Directivity.sofa``, specify mic +name at the end to retrive directivity pattern for that particular mic from the +eigenmike. This file contains 32 patterns of the form ``EM_32_*``, where ``*`` +is one of 0, 1, ..., 31. For example, ``EM_32_9`` will retrive pattern of mic +number "10" from the eigenmic. + +The ``LSPs_HATS_GuitarCabinets_Akustikmessplatz.sofa`` DIRPAT file includes +some source patterns. + +a) Genelec_8020 +b) Lambda_labs_CX-1A +c) HATS_4128C +d) Tannoy_System_1200 +e) Neumann_KH120A +f) Yamaha_DXR8 +g) BM_1x12inch_driver_closed_cabinet +h) BM_1x12inch_driver_open_cabinet +i) BM_open_stacked_on_closed_withCrossoverNetwork +j) BM_open_stacked_on_closed_fullrange +k) Palmer_1x12inch +l) Vibrolux_2x10inch +""" + +import matplotlib.pyplot as plt + +import pyroomacoustics as pra +from pyroomacoustics.directivities import ( + Cardioid, + DirectionVector, + FigureEight, + MeasuredDirectivityFile, + Rotation3D, +) + +# Reads the file containing the Eigenmike's directivity measurements +eigenmike = MeasuredDirectivityFile("EM32_Directivity", fs=16000) +# Reads the file containing the directivity measurements of another microphones +akg = MeasuredDirectivityFile("AKG_c480_c414_CUBE", fs=16000) + +# Create a rotation object to orient the microphones. +rot_54_73 = Rotation3D([73, 54], "yz", degrees=True) + +# Get the directivity objects from the two files +dir_obj_Dmic = akg.get_mic_directivity("AKG_c414K", orientation=rot_54_73) +dir_obj_Emic = eigenmike.get_mic_directivity("EM_32_9", orientation=rot_54_73) + +# Create two analytical directivities for comparison +dir_obj_Cmic = FigureEight( + orientation=DirectionVector(azimuth=90, colatitude=123, degrees=True), +) +dir_obj_Csrc = Cardioid( + orientation=DirectionVector(azimuth=56, colatitude=123, degrees=True), +) + + +room_dim = [6, 6, 2.4] + + +all_materials = { + "east": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "west": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "north": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "south": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "ceiling": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "floor": pra.Material( + energy_absorption={ + "coeffs": [0.11, 0.14, 0.37, 0.43, 0.27, 0.25], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), +} + + +room = pra.ShoeBox( + room_dim, + fs=16000, + max_order=20, + materials=pra.Material(0.5), + air_absorption=True, + ray_tracing=False, + min_phase=False, +) + +dir_mic = dir_obj_Emic + +room.add_source([1.52, 0.883, 1.044], directivity=dir_obj_Csrc) + + +room.add_microphone([2.31, 1.65, 1.163], directivity=dir_mic) + +dir_mic.set_orientation(Rotation3D([73, 54], rot_order="yz")) + + +room.compute_rir() +room.plot_rir(FD=True) + +# print(dir_mic.obj_open_sofa_inter.freq_angles_fft.shape) +# dir_mic.obj_open_sofa_inter.interpolate = False + +fig = plt.figure() +for idx, fb in enumerate(range(44)): + if idx >= 5 * 10: + break + ax = fig.add_subplot(5, 10, idx + 1, projection="3d") + dir_mic.plot(freq_bin=fb, ax=ax, depth=True) + ax.set_title(idx) +plt.show() + + +rir_1_0 = room.rir[0][0] diff --git a/examples/room_from_rt60.py b/examples/room_from_rt60.py index d33141e0..9ee7ebfc 100644 --- a/examples/room_from_rt60.py +++ b/examples/room_from_rt60.py @@ -39,10 +39,17 @@ # We invert Sabine's formula to obtain the parameters for the ISM simulator e_absorption, max_order = pra.inverse_sabine(rt60_tgt, room_dim) + pra.constants.set("octave_bands_keep_dc", True) + # Create the room if args.method == "ism": room = pra.ShoeBox( - room_dim, fs=fs, materials=pra.Material(e_absorption), max_order=max_order + room_dim, + fs=fs, + materials=pra.Material(e_absorption), + max_order=max_order, + use_rand_ism=True, + air_absorption=True, ) elif args.method == "hybrid": room = pra.ShoeBox( diff --git a/pyroomacoustics/acoustics.py b/pyroomacoustics/acoustics.py index 7d6f1864..26bc176f 100644 --- a/pyroomacoustics/acoustics.py +++ b/pyroomacoustics/acoustics.py @@ -31,7 +31,7 @@ import numpy as np from scipy.fftpack import dct from scipy.interpolate import interp1d -from scipy.signal import butter, fftconvolve, sosfiltfilt +from scipy.signal import butter, fftconvolve from .parameters import constants from .transform import stft @@ -159,10 +159,11 @@ class OctaveBandsFactory(object): Use third octave bands if True (default: False) """ - def __init__(self, base_frequency=125.0, fs=16000, n_fft=512): + def __init__(self, base_frequency=125.0, fs=16000, n_fft=512, keep_dc=False): self.base_freq = base_frequency self.fs = fs self.n_fft = n_fft + self.keep_dc = keep_dc # compute the number of bands self.n_bands = math.floor(np.log2(fs / base_frequency)) @@ -258,27 +259,7 @@ def __call__(self, coeffs=0.0, center_freqs=None, interp_kind="linear", **kwargs def _make_filters(self): """ - Create the band-pass filters for the octave bands - - Parameters - ---------- - order: int, optional - The order of the IIR filters (default: 8) - output: {'ba', 'zpk', 'sos'} - Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or - second-order sections ('sos'). Default is 'ba'. - - Returns - ------- - A list of callables that will each apply one of the band-pass filters - """ - - """ - filter_bank = bandpass_filterbank( - self.bands, fs=self.fs, order=order, output=output - ) - - return [lambda sig: sosfiltfilt(bpf, sig) for bpf in filter_bank] + Creates the band-pass filters for the octave bands """ # This seems to work only for Octave bands out of the box @@ -286,16 +267,27 @@ def _make_filters(self): n = len(self.centers) new_bands = [[centers[0] / 2, centers[1]]] + for i in range(1, n - 1): new_bands.append([centers[i - 1], centers[i + 1]]) new_bands.append([centers[-2], self.fs / 2]) n_freq = self.n_fft // 2 + 1 freq_resp = np.zeros((n_freq, n)) - freq = np.arange(n_freq) / self.n_fft * self.fs + + freq = ( + np.arange(n_freq) / self.n_fft * self.fs + ) # This only contains positive newfrequencies for b, (band, center) in enumerate(zip(new_bands, centers)): + if b == 0 and self.keep_dc: + # Converting Octave bands so that the minimum phase filters do not + # have ripples. + make_one = freq < center + freq_resp[make_one, b] = 1.0 + lo = np.logical_and(band[0] <= freq, freq < center) + freq_resp[lo, b] = 0.5 * (1 + np.cos(2 * np.pi * freq[lo] / center)) if b != n - 1: @@ -309,10 +301,12 @@ def _make_filters(self): np.fft.irfft(freq_resp, n=self.n_fft, axis=0), axes=[0], ) - # remove the first sample to make them odd-length symmetric filters self.filters = filters[1:, :] + # Octave band filters in frequency domain + self.filters_freq_domain = np.fft.fft(filters, axis=0, n=self.n_fft) + def critical_bands(): """ diff --git a/pyroomacoustics/beamforming.py b/pyroomacoustics/beamforming.py index f385d205..42cf2b19 100644 --- a/pyroomacoustics/beamforming.py +++ b/pyroomacoustics/beamforming.py @@ -24,20 +24,13 @@ from __future__ import division -import copy - import numpy as np import scipy.linalg as la from . import transform from . import utilities as u from . import windows -from .directivities import ( - CardioidFamily, - DirectionVector, - Directivity, - DirectivityPattern, -) +from .directivities import DirectionVector, Directivity, Omnidirectional from .parameters import constants from .soundsource import build_rir_matrix @@ -308,10 +301,7 @@ def circular_microphone_array_xyplane( if directivity is not None: assert isinstance(directivity, Directivity) else: - orientation = DirectionVector(azimuth=0, colatitude=colatitude, degrees=True) - directivity = CardioidFamily( - orientation=orientation, pattern_enum=DirectivityPattern.OMNI - ) + directivity = Omnidirectional() # for plotting azimuth_plot = None @@ -345,6 +335,7 @@ def circular_microphone_array_xyplane( # ========================================================================= # Classes (microphone array and beamformer related) # ========================================================================= +import copy class MicrophoneArray(object): @@ -372,15 +363,16 @@ def __init__(self, R, fs, directivity=None): self.R = R # array geometry self.fs = fs # sampling frequency of microphones - self.directivity = None - - if directivity is not None: - self.set_directivity(directivity) + self.set_directivity(directivity) self.signals = None self.center = np.mean(R, axis=1, keepdims=True) + @property + def is_directive(self): + return any([d is not None for d in self.directivity]) + def set_directivity(self, directivities): """ This functions sets self.directivity as a list of directivities with `n_mics` entries, @@ -399,7 +391,7 @@ def set_directivity(self, directivities): self.directivity = directivities else: # only 1 directivity specified - assert isinstance(directivities, Directivity) + assert directivities is None or isinstance(directivities, Directivity) self.directivity = [directivities] * self.nmic def record(self, signals, fs): @@ -462,8 +454,8 @@ def to_wav(self, filename, mono=False, norm=False, bitdepth=float): if true, normalize the signal to fit in the dynamic range (default `False`) bitdepth: int, optional - the format of output samples [np.int8/16/32/64 or np.float - (default)] + the format of output samples [np.int8/16/32/64 or + float (default)/np.float32/np.float64] """ from scipy.io import wavfile @@ -472,7 +464,7 @@ def to_wav(self, filename, mono=False, norm=False, bitdepth=float): else: signal = self.signals.T # each column is a channel - float_types = [float, np.float32, np.float64] + float_types = [float, float, np.float32, np.float64] if bitdepth in float_types: bits = None @@ -509,6 +501,7 @@ def append(self, locs): """ if isinstance(locs, MicrophoneArray): self.R = np.concatenate((self.R, locs.R), axis=1) + self.directivity += locs.directivity else: self.R = np.concatenate((self.R, locs), axis=1) diff --git a/pyroomacoustics/build_rir.pyx b/pyroomacoustics/build_rir.pyx index 5c0941d7..8cb119e8 100644 --- a/pyroomacoustics/build_rir.pyx +++ b/pyroomacoustics/build_rir.pyx @@ -57,29 +57,110 @@ def fast_rir_builder( # create a look-up table of the sinc function and # then use linear interpolation cdef float delta = 1. / lut_gran + + #Total number of points in 81 fractional delay cdef int lut_size = (fdl + 1) * lut_gran + 1 + + #equal space between -41 to +41 for 1641 length as each point between -40 to +40 represents 20 samples in the sinc n = np.linspace(-fdl2-1, fdl2 + 1, lut_size) - cdef double [:] sinc_lut = np.sinc(n) - cdef double [:] hann = np.hanning(fdl) + cdef double [:] sinc_lut = np.sinc(n) #Sinc over linspace n + cdef double [:] hann = np.hanning(fdl) #Hanning window of size 81 cdef int lut_pos, i, f, time_ip cdef float x_off, x_off_frac, sample_frac + g =[] + print_filter=0 + pf=[] + #Loop through each image source for i in range(n_times): if visibility[i] == 1: # decompose integer and fractional delay - sample_frac = fs * time[i] - time_ip = int(floor(sample_frac)) - time_fp = sample_frac - time_ip + sample_frac = fs * time[i] #Samples in fraction eg 250.567 sample , actual time of arrival of the image source to the microphone + time_ip = int(floor(sample_frac)) #Get the integer value of the sample eg 250th sample as int(250.567) = 250 + time_fp = sample_frac - time_ip #Get the fractional sample 250.567-250= 0.567 samples # do the linear interpolation - x_off_frac = (1. - time_fp) * lut_gran - lut_gran_off = int(floor(x_off_frac)) - x_off = (x_off_frac - lut_gran_off) - lut_pos = lut_gran_off + x_off_frac = (1. - time_fp) * lut_gran #(1-0.567) *20 = 8.66 , as each point represents 20 samples in the sinc table + lut_gran_off = int(floor(x_off_frac)) #int(8.66) = 8 sample in the sinc table + x_off = (x_off_frac - lut_gran_off) #fractional in the sinc table 8.66-8=0.66 + lut_pos = lut_gran_off #lut_pos=8 k = 0 + + #Loop through -40 to 41 , which accounts for 81 , as it is the amount of fractional delay every dirac goes through, the sinc table helps spread the energy to -40 to +40 samples in the RIR. for f in range(-fdl2, fdl2+1): - rir[time_ip + f] += alpha[i] * hann[k] * (sinc_lut[lut_pos] + rir[time_ip + f] += alpha[i] * hann[k] * (sinc_lut[lut_pos] + x_off * (sinc_lut[lut_pos+1] - sinc_lut[lut_pos])) + pf.append(rir[time_ip+f]) lut_pos += lut_gran k += 1 + +def fast_window_sinc_interpolator(double [:] vectorized_time_fp, int window_length, double [:,:] vectorized_interpolated_sinc): #Takes fractional part of the delay of IS k + + cdef double [:] hann_wd = np.hanning(window_length) + cdef int fdl2 = (window_length - 1) // 2 + cdef int lut_gran=20 + cdef int lut_size = (window_length + 1) * lut_gran + 1 + n_ = np.linspace(-fdl2 - 1, fdl2 + 1, lut_size) + cdef double [:] sinc_lut=np.sinc(n_) + cdef int img_src = 0 + + for time_fp in vectorized_time_fp: + x_off_frac = (1 - time_fp) * lut_gran + lut_gran_off = int(np.floor(x_off_frac)) + x_off = x_off_frac - lut_gran_off + lut_pos = lut_gran_off + filter_sample = 0 + + for f in range(-fdl2, fdl2 + 1): + vectorized_interpolated_sinc[img_src,filter_sample] = hann_wd[filter_sample]*(sinc_lut[lut_pos] + x_off * (sinc_lut[lut_pos + 1] - sinc_lut[lut_pos])) + lut_pos += lut_gran + filter_sample += 1 + + img_src+=1 + + return vectorized_interpolated_sinc + +cdef int val_i +import multiprocessing + +nthread = multiprocessing.cpu_count() + +def fast_convolution_4 ( + double complex [:] a, + double complex [:] b, + double complex [:] c, + double complex [:] d, + int final_fir_IS_len,): + + cdef double complex [:] out = np.zeros(final_fir_IS_len,dtype=np.complex_) + + a=np.fft.fft(a) + b=np.fft.fft(b) + c=np.fft.fft(c) + d=np.fft.fft(d) + + for val_i in range(final_fir_IS_len): + out[val_i]=a[val_i]*b[val_i]*c[val_i]*d[val_i] + + out=np.fft.ifft(out) + return out + + +def fast_convolution_3 ( + double complex [:] a, + double complex [:] b, + double complex [:] c, + int final_fir_IS_len,): + + cdef double complex [:] out = np.zeros(final_fir_IS_len,dtype=np.complex_) + + a=np.fft.fft(a) + b=np.fft.fft(b) + c=np.fft.fft(c) + + for val_i in range(final_fir_IS_len): + out[val_i]=a[val_i]*b[val_i]*c[val_i] + + out=np.fft.ifft(out) + return out diff --git a/pyroomacoustics/data/sofa/AKG_c480_c414_CUBE.sofa b/pyroomacoustics/data/sofa/AKG_c480_c414_CUBE.sofa new file mode 100644 index 00000000..bc111fff Binary files /dev/null and b/pyroomacoustics/data/sofa/AKG_c480_c414_CUBE.sofa differ diff --git a/pyroomacoustics/data/sofa/EM32_Directivity.sofa b/pyroomacoustics/data/sofa/EM32_Directivity.sofa new file mode 100644 index 00000000..7c03ba2b Binary files /dev/null and b/pyroomacoustics/data/sofa/EM32_Directivity.sofa differ diff --git a/pyroomacoustics/data/sofa/mit_kemar_large_pinna.sofa b/pyroomacoustics/data/sofa/mit_kemar_large_pinna.sofa new file mode 100644 index 00000000..16241403 Binary files /dev/null and b/pyroomacoustics/data/sofa/mit_kemar_large_pinna.sofa differ diff --git a/pyroomacoustics/data/sofa/mit_kemar_normal_pinna.sofa b/pyroomacoustics/data/sofa/mit_kemar_normal_pinna.sofa new file mode 100644 index 00000000..deea7917 Binary files /dev/null and b/pyroomacoustics/data/sofa/mit_kemar_normal_pinna.sofa differ diff --git a/pyroomacoustics/data/sofa_files.json b/pyroomacoustics/data/sofa_files.json new file mode 100644 index 00000000..26ad6e08 --- /dev/null +++ b/pyroomacoustics/data/sofa_files.json @@ -0,0 +1,120 @@ +{ + "Oktava_MK4012_CUBE": { + "supported": false, + "type": "microphones", + "url": "https://fedora.kug.ac.at/fedora/objects/o:68229/methods/bdef:Container/get?param0=Oktava_MK4012_CUBE.sofa", + "homepage": "https://phaidra.kug.ac.at/detail/o:68229", + "license": "CC0", + "contains": [] + }, + + "LSPs_HATS_GuitarCabinets_Akustikmessplatz": { + "supported": true, + "type": "sources", + "url": "https://fedora.kug.ac.at/fedora/objects/o:68229/methods/bdef:Container/get?param0=LSPs_HATS_GuitarCabinets_Akustikmessplatz.sofa", + "homepage": "https://phaidra.kug.ac.at/detail/o:68229", + "license": "CC0", + "contains": [ + "Genelec_8020", + "Lambda_labs_CX-1A", + "HATS_4128C", + "Tannoy_System_1200", + "Neumann_KH120A", + "Yamaha_DXR8", + "BM_1x12inch_driver_closed_cabinet", + "BM_1x12inch_driver_open_cabinet", + "BM_open_stacked_on_closed_withCrossoverNetwork", + "BM_open_stacked_on_closed_fullrange", + "Palmer_1x12inch", + "Vibrolux_2x10inch" + ] + }, + + "Soundfield_ST450_CUBE": { + "supported": false, + "type": "microphones", + "url": "https://fedora.kug.ac.at/fedora/objects/o:68229/methods/bdef:Container/get?param0=Soundfield_ST450_CUBE.sofa", + "homepage": "https://phaidra.kug.ac.at/detail/o:68229", + "license": "CC0", + "contains": [] + }, + + "AKG_c480_c414_CUBE": { + "supported": true, + "type": "microphones", + "url": "https://fedora.kug.ac.at/fedora/objects/o:68229/methods/bdef:Container/get?param0=AKG_c480_c414_CUBE.sofa", + "homepage": "https://phaidra.kug.ac.at/detail/o:68229", + "license": "CC0", + "contains": [ + "AKG_c480", + "AKG_c414K", + "AKG_c414N", + "AKG_c414S", + "AKG_c414A" + ] + }, + + "EM32_Directivity": { + "supported": true, + "type": "microphones", + "url": "https://phaidra.kug.ac.at/download/o:69292", + "homepage": "https://phaidra.kug.ac.at/detail/o:69292", + "license": "CC0", + "contains": [ + "EM_32_0", + "EM_32_1", + "EM_32_2", + "EM_32_3", + "EM_32_4", + "EM_32_5", + "EM_32_6", + "EM_32_7", + "EM_32_8", + "EM_32_9", + "EM_32_10", + "EM_32_11", + "EM_32_12", + "EM_32_13", + "EM_32_14", + "EM_32_15", + "EM_32_16", + "EM_32_17", + "EM_32_18", + "EM_32_19", + "EM_32_20", + "EM_32_21", + "EM_32_22", + "EM_32_23", + "EM_32_24", + "EM_32_25", + "EM_32_26", + "EM_32_27", + "EM_32_28", + "EM_32_29", + "EM_32_30", + "EM_32_31" + ] + }, + "mit_kemar_large_pinna": { + "supported": true, + "type": "microphones", + "url": "http://sofacoustics.org/data/database/mit/mit_kemar_large_pinna.sofa", + "homepage": "https://sound.media.mit.edu/resources/KEMAR/", + "license": "This data is Copyright 1994 by the MIT Media Laboratory. It is provided free with no restrictions on use, provided the authors are cited when the data is used in any research or commercial application.", + "contains": [ + "left", + "right" + ] + }, + "mit_kemar_normal_pinna": { + "supported": true, + "type": "microphones", + "url": "http://sofacoustics.org/data/database/mit/mit_kemar_normal_pinna.sofa", + "homepage": "https://sound.media.mit.edu/resources/KEMAR/", + "license": "This data is Copyright 1994 by the MIT Media Laboratory. It is provided free with no restrictions on use, provided the authors are cited when the data is used in any research or commercial application.", + "contains": [ + "left", + "right" + ] + } +} diff --git a/pyroomacoustics/datasets/__init__.py b/pyroomacoustics/datasets/__init__.py index a8cea886..5673058c 100644 --- a/pyroomacoustics/datasets/__init__.py +++ b/pyroomacoustics/datasets/__init__.py @@ -140,4 +140,5 @@ from .base import AudioSample, Dataset, Meta, Sample from .cmu_arctic import CMUArcticCorpus, CMUArcticSentence, cmu_arctic_speakers from .google_speech_commands import GoogleSample, GoogleSpeechCommands +from .sofa import SOFADatabase, download_sofa_files from .timit import Sentence, TimitCorpus, Word diff --git a/pyroomacoustics/datasets/sofa.py b/pyroomacoustics/datasets/sofa.py new file mode 100644 index 00000000..fd07fbfc --- /dev/null +++ b/pyroomacoustics/datasets/sofa.py @@ -0,0 +1,272 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2024 Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Pyroomacoustics contains a small database of SOFA files that have been tested +and can be used for simultions. +The database can be loaded using the function +:py:class:`~pyroomacoustics.datasets.sofa.SOFADatabase`. + +.. code-block:: python + + # load and display the list of available SOFA files and their content + from pyroomacoustics.datasets import SOFADatabase + + db = SOFADatabase() + db.list() + +The database contains the following files. + +- Three files from the `DIRPAT database + `_ + collected by Manuel Brandner, Matthias Frank, and Daniel Rudrich University + of Music and Performing Arts, Graz. + + - ``AKG_c480_c414_CUBE.sofa`` containing the directive responses of a + microphone with 4 different patterns. + - ``EM32_Directivity.sofa`` that contains the directional response of the + `Eigenmike em32 `_ microphone array. + - ``LSPs_HATS_GuitarCabinets_Akustikmessplatz.sofa`` that contains 12 source + directivities. This file is dynamically downloaded upon its first use. + - The files are public domain + (`CC0 `_), + but if you use them in your research, please cite the following + `paper `_. + + :: + + M. Brandner, M. Frank, and D. Rudrich, "DirPat—Database and + Viewer of 2D/3D Directivity Patterns of Sound Sources and Receivers," + in Audio Engineering Society Convention 144, Paper 425, 2018. + +- Two head-related transfer functions of the MIT KEMAR dummy head + with normal and large pinna. The data was collected by Bill Gardner + and Keith Martin from MIT and is free to use provided the authors are + cited. See the + `full license `_ + for more details. + + +""" + + +import json +import typing as tp +from dataclasses import dataclass +from pathlib import Path + +from .utils import AttrDict, download_multiple + +_pra_data_folder = Path(__file__).parents[1] / "data" +DEFAULT_SOFA_PATH = _pra_data_folder / "sofa" +SOFA_INFO = _pra_data_folder / "sofa_files.json" + +_DIRPAT_FILES = [ + "Soundfield_ST450_CUBE", + "AKG_c480_c414_CUBE", + "Oktava_MK4012_CUBE", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", +] + + +def is_dirpat(name): + if isinstance(name, Path): + name = name.stem + return name in _DIRPAT_FILES + + +def get_sofa_db(): + """ + A helper function to quickly load the SOFA database + """ + # we want to avoid loading the database multiple times + global sofa_db + try: + return sofa_db + except NameError: + sofa_db = SOFADatabase() + return sofa_db + + +def resolve_sofa_path(path): + path = Path(path) + + if path.exists(): + return path + + sofa_db = get_sofa_db() + if path.stem in sofa_db: + return Path(sofa_db[path.stem].path) + + raise ValueError(f"SOFA file {path} could not be found") + + +def get_sofa_db_info(): + with open(SOFA_INFO, "r") as f: + sofa_info = json.load(f) + return sofa_info + + +def download_sofa_files(path=None, overwrite=False, verbose=False, no_fail=False): + """ + Download the SOFA files containing source/receiver impulse responses + + Parameters + ---------- + path: str or Path, optional + A path to a directory where the files will be downloaded + overwrite: bool, optional + If set to `True`, forces the download even if the files already exist + verbose: bool + Print some information about the download status (default `False`) + + Returns + ------- + files: list of Path + The list of path to the files downloaded + """ + if path is None: + path = DEFAULT_SOFA_PATH + + path = Path(path) + + path.mkdir(exist_ok=True, parents=True) + + sofa_info = get_sofa_db_info() + + files = { + path / f"{name}.sofa": info["url"] + for name, info in sofa_info.items() + if info["supported"] + } + download_multiple(files, overwrite=overwrite, verbose=verbose, no_fail=no_fail) + + return list(files.keys()) + + +@dataclass +class SOFAFileInfo: + """ + A class to store information about a SOFA file + + Parameters + ---------- + path: Path + The path to the SOFA file + supported: bool + Whether the SOFA file is supported by Pyroom Acoustics + type: str + The type of device (e.g., 'sources' or 'microphones') + url: str + The URL where the SOFA file can be downloaded + homepage: str + The URL of the SOFA file homepage + license: str + The license of the SOFA file + contains: List[str] + The labels of the sources/microphones contained in the SOFA file, + or``None`` if the information is not available + """ + + path: Path + supported: bool = True + type: str = "unknown" + url: str = "unknown" + homepage: str = "unknown" + license: str = "unknown" + contains: tp.List[str] = None + + +class SOFADatabase(dict): + """ + A small database of SOFA files containing source/microphone directional + impulse responses + + The database object is a dictionary-like object where the keys are the + names of the SOFA files and the values are objects with the following + attributes: + + .. code-block:: python + + db = SOFADatabase() + + # type of device: 'sources' or 'microphones' + db["Soundfield_ST450_CUBE"].type + + # list of the labels of the sources/microphones + db["Soundfield_ST450_CUBE"].contains + + + Parameters + ---------- + download: bool, optional + If set to `True`, the SOFA files are downloaded if they are not already + present in the default folder + """ + + def __init__(self, download=True): + super().__init__() + + if download: + # specify "no_fail" to avoid errors if internet is not available + download_sofa_files(path=self.root, no_fail=True) + + self._db = {} + for name, info in get_sofa_db_info().items(): + path = self.root / f"{name}.sofa" + if path.exists(): + dict.__setitem__(self, name, SOFAFileInfo(path=path, **info)) + + for path in DEFAULT_SOFA_PATH.glob("*.sofa"): + name = path.stem + if name not in self: + dict.__setitem__( + self, + name, + SOFAFileInfo(path=path), + ) + + def list(self): + """ + Print a list of the available SOFA files and the labels of the + different devices they contain + """ + for name, info in self.items(): + print(f"- {name} ({info.type})") + if info.contains is not None: + for channel in info.contains: + print(f" - {channel}") + + @property + def root(self): + """The path to the folder containing the SOFA files""" + return DEFAULT_SOFA_PATH + + @property + def db_info_path(self): + """The path to the JSON file containing the SOFA files information""" + return SOFA_INFO + + def __setitem__(self, key, val): + # disallow writing elements + raise RuntimeError(f"{self.__class__} is not writable") diff --git a/pyroomacoustics/datasets/utils.py b/pyroomacoustics/datasets/utils.py index d106cb3c..44448fec 100644 --- a/pyroomacoustics/datasets/utils.py +++ b/pyroomacoustics/datasets/utils.py @@ -25,11 +25,40 @@ import bz2 import os import tarfile +from pathlib import Path try: - from urllib.request import urlopen + from urllib.request import urlopen, urlretrieve except ImportError: - from urllib import urlopen + # support for python 2.7, should be able to remove by now + from urllib import urlopen, urlretrieve + + +class AttrDict(object): + """Convert a dictionary into an object""" + + def __init__(self, dictionary): + for key, val in dictionary.items(): + if isinstance(val, dict): + setattr(self, key, AttrDict(val)) + elif isinstance(val, list): + setattr( + self, key, [AttrDict(v) if isinstance(v, dict) else v for v in val] + ) + else: + setattr(self, key, val) + + def __getitem__(self, key): + return getattr(self, key) + + def __setitem__(self, key, val): + return setattr(self, key, val) + + def __str__(self): + return str(self.__dict__) + + def __repr__(self): + return repr(self.__dict__) def download_uncompress(url, path=".", compression=None, context=None): @@ -69,3 +98,41 @@ def download_uncompress(url, path=".", compression=None, context=None): stream = urlopen(url) tf = tarfile.open(fileobj=stream, mode=mode) tf.extractall(path) + + +def download_multiple(files_dict, overwrite=False, verbose=False, no_fail=False): + """ + A utility to download multiple files + + Parameters + ---------- + files_dict: dict + A dictionary of files to download with key=local_path and value=url + overwrite: bool + If `True` if the local file exists, it will be overwritten. If `False` + (default), existing files are skipped. + """ + skip_ok = not overwrite + + for path, url in files_dict.items(): + path = Path(path) + if path.exists() and skip_ok: + if verbose: + print( + f"{path} exists: skip. Use `overwrite` option to download anyway." + ) + continue + + if verbose: + print(f"Download {url} -> {path}...", end="") + + try: + urlretrieve(url, path) + except URLError: + if no_fail: + continue + else: + raise URLError(f"Failed to download {url}") + + if verbose: + print(" done.") diff --git a/pyroomacoustics/directivities.py b/pyroomacoustics/directivities.py deleted file mode 100644 index d900bd47..00000000 --- a/pyroomacoustics/directivities.py +++ /dev/null @@ -1,370 +0,0 @@ -import abc -from enum import Enum - -import numpy as np - -from pyroomacoustics.doa import spher2cart -from pyroomacoustics.utilities import all_combinations, requires_matplotlib - - -class DirectivityPattern(Enum): - """ - Common cardioid patterns and their corresponding coefficient for the expression: - - .. math:: - - r = a (1 + \\cos \\theta), - - where :math:`a` is the coefficient that determines the cardioid pattern and :math:`r` yields - the gain at angle :math:`\\theta`. - - """ - - FIGURE_EIGHT = 0 - HYPERCARDIOID = 0.25 - CARDIOID = 0.5 - SUBCARDIOID = 0.75 - OMNI = 1.0 - - -class DirectionVector(object): - """ - Object for representing direction vectors in 3D, parameterized by an azimuth and colatitude - angle. - - 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 - - self._unit_v = np.array( - [ - np.cos(self._azimuth) * np.sin(self._colatitude), - np.sin(self._azimuth) * np.sin(self._colatitude), - np.cos(self._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 - - @property - def unit_vector(self): - """Direction vector in cartesian coordinates.""" - return self._unit_v - - -class Directivity(abc.ABC): - """ - Abstract class for directivity patterns. - - """ - - def __init__(self, orientation): - assert isinstance(orientation, DirectionVector) - self._orientation = orientation - - 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 set_orientation(self, orientation): - """ - Set orientation of directivity pattern. - - Parameters - ---------- - orientation : DirectionVector - New direction for the directivity pattern. - """ - assert isinstance(orientation, DirectionVector) - self._orientation = orientation - - @abc.abstractmethod - def get_response( - self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True - ): - """ - Get response for provided angles and frequency. - """ - return - - -class CardioidFamily(Directivity): - """ - Object for directivities coming from the - `cardioid family `_. - - Parameters - ---------- - orientation : DirectionVector - Indicates direction of the pattern. - pattern_enum : DirectivityPattern - Desired pattern for the cardioid. - """ - - def __init__(self, orientation, pattern_enum, gain=1.0): - Directivity.__init__(self, orientation) - self._p = pattern_enum.value - self._gain = gain - self._pattern_name = pattern_enum.name - - @property - def directivity_pattern(self): - """Name of cardioid directivity pattern.""" - return self._pattern_name - - def get_response( - self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True - ): - """ - Get response for provided angles. - - Parameters - ---------- - azimuth : array_like - Azimuth in degrees - colatitude : array_like, optional - Colatitude in degrees. Default is to be on XY plane. - magnitude : bool, optional - Whether to return magnitude of response. - frequency : float, optional - For which frequency to compute the response. Cardioid are frequency-independent so this - value has no effect. - degrees : bool, optional - Whether provided angles are in degrees. - - - Returns - ------- - resp : :py:class:`~numpy.ndarray` - Response at provided angles. - """ - - if colatitude is not None: - assert len(azimuth) == len(colatitude) - if self._p == DirectivityPattern.OMNI: - return np.ones(len(azimuth)) - else: - coord = spher2cart(azimuth=azimuth, colatitude=colatitude, degrees=degrees) - resp = self._gain * self._p + (1 - self._p) * np.matmul( - self._orientation.unit_vector, coord - ) - if magnitude: - return np.abs(resp) - else: - return resp - - @requires_matplotlib - def plot_response( - self, azimuth, colatitude=None, degrees=True, ax=None, offset=None - ): - """ - Plot directivity response at specified angles. - - 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). - ax : axes object - offset : list - 3-D coordinates of the point where the response needs to be plotted. - - Return - ------ - ax : :py:class:`~matplotlib.axes.Axes` - """ - import matplotlib.pyplot as plt - - if offset is not None: - x_offset = offset[0] - y_offset = offset[1] - else: - x_offset = 0 - y_offset = 0 - - if degrees: - azimuth = np.radians(azimuth) - - if colatitude is not None: - if degrees: - colatitude = np.radians(colatitude) - - if ax is None: - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - - if offset is not None: - z_offset = offset[2] - else: - z_offset = 0 - - spher_coord = all_combinations(azimuth, colatitude) - azi_flat = spher_coord[:, 0] - col_flat = spher_coord[:, 1] - - # compute response - resp = self.get_response( - azimuth=azi_flat, colatitude=col_flat, magnitude=True, degrees=False - ) - RESP = resp.reshape(len(azimuth), len(colatitude)) - - # create surface plot, need cartesian coordinates - AZI, COL = np.meshgrid(azimuth, colatitude) - X = RESP.T * np.sin(COL) * np.cos(AZI) + x_offset - Y = RESP.T * np.sin(COL) * np.sin(AZI) + y_offset - Z = RESP.T * np.cos(COL) + z_offset - - ax.plot_surface(X, Y, Z) - - if ax is None: - ax.set_title( - "{}, azimuth={}, colatitude={}".format( - self.directivity_pattern, - self.get_azimuth(), - self.get_colatitude(), - ) - ) - else: - ax.set_title("Directivity Plot") - - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - - else: - if ax is None: - fig = plt.figure() - ax = plt.subplot(111) - - # compute response - resp = self.get_response(azimuth=azimuth, magnitude=True, degrees=False) - RESP = resp - - # create surface plot, need cartesian coordinates - X = RESP.T * np.cos(azimuth) + x_offset - Y = RESP.T * np.sin(azimuth) + y_offset - ax.plot(X, Y) - - return ax - - -def cardioid_func(x, direction, coef, gain=1.0, normalize=True, magnitude=False): - """ - One-shot function for computing cardioid response. - - Parameters - ----------- - x: array_like, shape (..., n_dim) - Cartesian coordinates - direction: array_like, shape (n_dim) - Direction vector, should be normalized. - coef: float - Parameter for the cardioid function. - gain: float - The gain. - normalize : bool - Whether to normalize coordinates and direction vector. - magnitude : bool - Whether to return magnitude, default is False. - - Returns - ------- - resp : :py:class:`~numpy.ndarray` - Response at provided angles for the speficied cardioid function. - """ - assert coef >= 0.0 - assert coef <= 1.0 - - # normalize positions - if normalize: - x /= np.linalg.norm(x, axis=0) - direction /= np.linalg.norm(direction) - - # compute response - resp = gain * (coef + (1 - coef) * np.matmul(direction, x)) - if magnitude: - return np.abs(resp) - else: - return resp - - -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 diff --git a/pyroomacoustics/directivities/__init__.py b/pyroomacoustics/directivities/__init__.py new file mode 100644 index 00000000..671b7a25 --- /dev/null +++ b/pyroomacoustics/directivities/__init__.py @@ -0,0 +1,95 @@ +# Directivity module that provides routines to use analytic and mesured directional +# responses for sources and microphones. +# Copyright (C) 2024 Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Real-world microphones and sound sources usually exhibit directional responses. +That is, the impulse response (or frequency response) depends on the emission +or reception angle (for sources and microphones, respectively). +A concrete example is the human ear attached to the head. The left ear is +typically more sensitive to sounds coming from the left side than from the right. + +This sub-module provides an interface to add such directional responses to +microphones and sources in the room impulse response simulation. + +.. warning:: + The directional responses are currently only supported for the + image source method based simulation. + +.. warning:: + Directional responses are only supported for 3D rooms. + + +The directivities are described by an object of a class derived from :py:class:`~pyroomacoustics.directivities.base.Directivity`. + +Let's dive right in with an example. +Here, we simulate a shoebox room with a cardioid source and a dummy head +receiver with two ears (i.e., microphones). This simulates a binaural response. + +.. code-block:: python + + import pyroomacoustics as pra + + room = pra.ShoeBox( + p=[5, 3, 3], + materials=pra.Material(energy_absorption), + fs=16000, + max_order=40, + ) + + # add a cardioid source + dir = pra.directivities.Cardioid(DirectionVector(azimuth=-65, colatitude=90) , gain=1.0) + room.add_source([3.75, 2.13, 1.41], directivity=dir) + + # add a dummy head receiver from the MIT KEMAR database + hrtf = MeasuredDirectivityFile( + path="mit_kemar_normal_pinna.sofa", # SOFA file is in the database + fs=room.fs, + interp_order=12, # interpolation order + interp_n_points=1000, # number of points in the interpolation grid + ) + + # provide the head rotation + orientation = Rotation3D([90.0, 30.0], "yz", degrees=True) + + # choose and interpolate the directivities + dir_left = hrtf.get_mic_directivity("left", orientation=orientation) + dir_right = hrtf.get_mic_directivity("right", orientation=orientation) + + # for a head-related transfer function, the microphone should be co-located + mic_pos = [1.05, 1.74, 1.81] + room.add_microphone(mic_pos, directivity=dir) + room.add_microphone(mic_pos, directivity=dir) +""" +from .analytic import ( + Cardioid, + CardioidFamily, + FigureEight, + HyperCardioid, + Omnidirectional, + SubCardioid, + cardioid_func, +) +from .base import Directivity +from .direction import DirectionVector, Rotation3D +from .measured import MeasuredDirectivity, MeasuredDirectivityFile diff --git a/pyroomacoustics/directivities/analytic.py b/pyroomacoustics/directivities/analytic.py new file mode 100644 index 00000000..59a3b28d --- /dev/null +++ b/pyroomacoustics/directivities/analytic.py @@ -0,0 +1,445 @@ +# Directivity module that provides routines to use analytic and mesured directional +# responses for sources and microphones. +# Copyright (C) 2020-2024 Robin Scheibler, Satvik Dixit, Eric Bezzam +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +A class of directional responses can be defined analytically. +Such respones include in particular the cardioid family of patterns +that describes cardioid, super-cardioid, and figure-of-eight microphones, (see +`cardioid family `_, under Polar patterns, with cardioid, hypercardioid, cardioid, subcardioid, figure-eight, and omnidirectional). +In three dimensions, for an orientation given by unit vector :math:`\boldsymbol{u}`, a parameter :math:`p \in [0, 1]`, +and a gain :math:`G`, the response to direction :math:`\boldsymbol{r}` (also a unit vector) is given by the following equation. + +.. math:: + f(\boldsymbol{r}\,;\,\boldsymbol{d}, p, G) = G (p + (1 - p) \boldsymbol{d}^\top \boldsymbol{r}), + +Note that :math:`\boldsymbol{d}^\top \boldsymbol{r}` is the inner product of two unit +vectors, that is, the cosine of the angle between them. + +Different values of :math:`p` correspond to different patterns: 0 for +figure-eight, 0.25 for hyper-cardioid, 0.5 for cardioid, 0.75 for +sub-cardioid, and 1 for omni. + +Specialized objects +:py:class:`~pyroomacoustics.directivities.analytic.Cardioid`, +:py:class:`~pyroomacoustics.directivities.analytic.FigureEight`, +:py:class:`~pyroomacoustics.directivities.analytic.SubCardioid`, +:py:class:`~pyroomacoustics.directivities.analytic.HyperCardioid`, +and :py:class:`~pyroomacoustics.directivities.analytic.Omnidirectional` are provided +for the different patterns. +The class :py:class:`~pyroomacoustics.directivities.analytic.CardioidFamily` can be used to make +a pattern with arbitrary parameter :math:`p`. + +.. code-block:: python + + # a cardioid pointing toward the ``z`` direction + from pyroomacoustics.directivities import CardioidFamily + + dir = Cardioid([0, 0, 1], gain=1.0) +""" +import numpy as np + +from ..doa import spher2cart +from ..utilities import all_combinations, requires_matplotlib +from .base import Directivity +from .direction import DirectionVector + +_FIGURE_EIGHT = 0 +_HYPERCARDIOID = 0.25 +_CARDIOID = 0.5 +_SUBCARDIOID = 0.75 +_OMNI = 1.0 + + +class CardioidFamily(Directivity): + r""" + Object for directivities coming from the + `cardioid family`_. + In three dimensions, for an orientation given by unit vector :math:`\\boldsymbol{u}`, a parameter :math:`p \in [0, 1]`, + and a gain :math:`G`, the pattern is given by the following equation. + + .. math:: + f(\boldsymbol{r}\,;\,\boldsymbol{d}, p, G) = G (p + (1 - p) \boldsymbol{d}^\top \boldsymbol{r}), + + Different values of :math:`p` correspond to different patterns: 0 for + figure-eight, 0.25 for hyper-cardioid, 0.5 for cardioid, 0.75 for + sub-cardioid, and 1 for omni. + + Note that all the patterns are cylindrically symmetric around the + orientation vector. + + Parameters + ---------- + orientation: DirectionVector or numpy.ndarray + Indicates direction of the pattern. + 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 + gain: float + The linear gain of the directivity pattern (default is 1.0) + """ + + def __init__(self, orientation, p, gain=1.0): + if isinstance(orientation, list) or hasattr(orientation, "__array__"): + orientation = np.array(orientation) + # check if valid direction vector, normalize, make object + if orientation.shape != (3,): + raise ValueError("Orientation must be a 3D vector.") + orientation = orientation / np.linalg.norm(orientation) + azimuth, colatitude = spher2cart(orientation) # returns radians + self._orientation = DirectionVector(azimuth, colatitude, degrees=False) + elif isinstance(orientation, DirectionVector): + self._orientation = orientation + else: + raise ValueError("Orientation must be a DirectionVector or a 3D vector.") + + self._p = p + if not 0 <= self._p <= 1: + raise ValueError("The parameter p must be between 0 and 1.") + + self._gain = gain + self._pattern_name = f"cardioid family, p={self._p}" + + @property + def is_impulse_response(self): + # this is not an impulse response, do not make docstring to avoid clutter in the + # documentation + return False + + @property + def filter_len_ir(self): + # no impulse response means length 1 + return 1 + + @property + def directivity_pattern(self): + """Name of cardioid directivity pattern.""" + return self._pattern_name + + def get_azimuth(self, degrees=True): + return self._orientation.get_azimuth(degrees, degrees=degrees) + + def get_colatitude(self, degrees=True): + return self._orientation.get_colatitude(degrees, degrees=degrees) + + def set_orientation(self, orientation): + """ + Set orientation of directivity pattern. + + Parameters + ---------- + orientation : DirectionVector + New direction for the directivity pattern. + """ + assert isinstance(orientation, DirectionVector) + self._orientation = orientation + + def get_response( + self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True + ): + """ + Get response for provided angles. + + Parameters + ---------- + azimuth : array_like + Azimuth + colatitude : array_like, optional + Colatitude. Default is to be on XY plane. + magnitude : bool, optional + Whether to return magnitude of response. + frequency : float, optional + For which frequency to compute the response. Cardioid are frequency-independent so this + value has no effect. + degrees : bool, optional + If ``True``, ``azimuth`` and ``colatitude`` are in degrees. + Otherwise, they are in radians. + + Returns + ------- + resp : :py:class:`~numpy.ndarray` + Response at provided angles. + """ + + if colatitude is not None: + assert len(azimuth) == len(colatitude) + if self._p == 1.0: + return self._gain * np.ones(len(azimuth)) + else: + coord = spher2cart(azimuth=azimuth, colatitude=colatitude, degrees=degrees) + + resp = self._gain * ( + self._p + + (1 - self._p) * np.matmul(self._orientation.unit_vector, coord) + ) + + if magnitude: + return np.abs(resp) + else: + return resp + + @requires_matplotlib + def plot_response( + self, azimuth, colatitude=None, degrees=True, ax=None, offset=None + ): + """ + Plot directivity response at specified angles. + + 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). + ax : axes object + offset : list + 3-D coordinates of the point where the response needs to be plotted. + + Return + ------ + ax : :py:class:`~matplotlib.axes.Axes` + """ + import matplotlib.pyplot as plt + + if offset is not None: + x_offset = offset[0] + y_offset = offset[1] + else: + x_offset = 0 + y_offset = 0 + + if degrees: + azimuth = np.radians(azimuth) + + if colatitude is not None: + if degrees: + colatitude = np.radians(colatitude) + + if ax is None: + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1, projection="3d") + + if offset is not None: + z_offset = offset[2] + else: + z_offset = 0 + + spher_coord = all_combinations(azimuth, colatitude) + azi_flat = spher_coord[:, 0] + col_flat = spher_coord[:, 1] + + # compute response + resp = self.get_response( + azimuth=azi_flat, colatitude=col_flat, magnitude=True, degrees=False + ) + RESP = resp.reshape(len(azimuth), len(colatitude)) + + # create surface plot, need cartesian coordinates + AZI, COL = np.meshgrid(azimuth, colatitude) + X = RESP.T * np.sin(COL) * np.cos(AZI) + x_offset + Y = RESP.T * np.sin(COL) * np.sin(AZI) + y_offset + Z = RESP.T * np.cos(COL) + z_offset + + ax.plot_surface(X, Y, Z) + + if ax is None: + ax.set_title( + "{}, azimuth={}, colatitude={}".format( + self.directivity_pattern, + self.get_azimuth(), + self.get_colatitude(), + ) + ) + else: + ax.set_title("Directivity Plot") + + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + + else: + if ax is None: + fig = plt.figure() + ax = plt.subplot(111) + + # compute response + resp = self.get_response(azimuth=azimuth, magnitude=True, degrees=False) + RESP = resp + + # create surface plot, need cartesian coordinates + X = RESP.T * np.cos(azimuth) + x_offset + Y = RESP.T * np.sin(azimuth) + y_offset + ax.plot(X, Y) + + return ax + + +class Cardioid(CardioidFamily): + """ + Cardioid directivity pattern. + + Parameters + ---------- + orientation: DirectionVector + Indicates direction of the pattern. + gain: float + The linear gain of the directivity pattern (default is 1.0) + """ + + def __init__(self, orientation, gain=1.0): + super().__init__(orientation, p=_CARDIOID, gain=gain) + self._pattern_name = "cardioid" + + +class FigureEight(CardioidFamily): + """ + Figure-of-eight directivity pattern. + + Parameters + ---------- + orientation: DirectionVector + Indicates direction of the pattern. + gain: float + The linear gain of the directivity pattern (default is 1.0) + """ + + def __init__(self, orientation, gain=1.0): + super().__init__(orientation, p=_FIGURE_EIGHT, gain=gain) + self._pattern_name = "figure-eight" + + +class SubCardioid(CardioidFamily): + """ + Sub-cardioid directivity pattern. + + Parameters + ---------- + orientation: DirectionVector + Indicates direction of the pattern. + gain: float + The linear gain of the directivity pattern (default is 1.0) + """ + + def __init__(self, orientation, gain=1.0): + super().__init__(orientation, p=_SUBCARDIOID, gain=gain) + self._pattern_name = "sub-cardioid" + + +class HyperCardioid(CardioidFamily): + """ + Hyper-cardioid directivity pattern. + + Parameters + ---------- + orientation: DirectionVector + Indicates direction of the pattern. + gain: float + The linear gain of the directivity pattern (default is 1.0) + """ + + def __init__(self, orientation, gain=1.0): + CardioidFamily.__init__(self, orientation, p=_HYPERCARDIOID, gain=gain) + self._pattern_name = "hyper-cardioid" + + +class Omnidirectional(CardioidFamily): + """ + Hyper-cardioid directivity pattern. + + Parameters + ---------- + orientation: DirectionVector + Indicates direction of the pattern. + gain: float + The linear gain of the directivity pattern (default is 1.0) + """ + + def __init__(self, gain=1.0): + CardioidFamily.__init__(self, DirectionVector(0.0, 0.0), p=_OMNI, gain=gain) + self._pattern_name = "omni" + + +def cardioid_func(x, direction, p, gain=1.0, normalize=True, magnitude=False): + """ + One-shot function for computing cardioid response. + + Parameters + ----------- + x: array_like, shape (n_dim, ...) + Cartesian coordinates + direction: array_like, shape (n_dim) + Direction vector, should be normalized + p: float + Parameter for the cardioid function (between 0 and 1) + gain: float + The gain + normalize : bool + Whether to normalize coordinates and direction vector + magnitude : bool + Whether to return magnitude, default is False + + Returns + ------- + resp : :py:class:`~numpy.ndarray` + Response at provided angles for the speficied cardioid function. + """ + if not 0.0 <= p <= 1.0: + raise ValueError("The parameter p must be between 0 and 1.") + + # normalize positions + if normalize: + x /= np.linalg.norm(x, axis=0) + direction /= np.linalg.norm(direction) + + # compute response + resp = gain * (p + (1 - p) * np.matmul(direction, x)) + if magnitude: + return np.abs(resp) + else: + return resp + + +def cardioid_energy(p, gain=1.0): + r""" + This function gives the exact value of the surface integral of the cardioid + (family) function on the unit sphere + + .. math:: + + E(p, G) = \iint_{\mathbb{S}^2} G^2 \left( p + (1 - p) \boldsymbol{d}^\top \boldsymbol{r} \right)^2 d\boldsymbol{r} + = \frac{4 \pi}{3} G^2 \left( 4 p^2 - 2 p + 1 \right). + + This can be used to normalize the energy sent/received. + + Parameters + --------- + p: float + The parameter of the cardioid function (between 0 and 1) + gain: float + The gain of the cardioid function + """ + return gain**2 * (4.0 * np.pi / 3.0) * (4 * p**2 - 2 * p + 1) diff --git a/pyroomacoustics/directivities/base.py b/pyroomacoustics/directivities/base.py new file mode 100644 index 00000000..349f25f9 --- /dev/null +++ b/pyroomacoustics/directivities/base.py @@ -0,0 +1,94 @@ +# Directivity module that provides routines to use analytic and mesured directional +# responses for sources and microphones. +# Copyright (C) 2020-2024 Robin Scheibler, Satvik Dixit, Eric Bezzam +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +The `Directivity` class is an abstract class that can be subclassed to create new +types of directivities. The class provides a common interface to access the +directivity patterns. + +The class should implement the following methods: + +- ``get_response`` to get the response for a given angle and frequency +- ``is_impulse_response`` to indicate whether the directivity is an impulse response or + just band coefficients +- ``filter_len_ir`` to return the length of the impulse response. This should return 1 + if the directivity is not an impulse response. +""" +import abc + + +class Directivity(abc.ABC): + """ + Abstract class for directivity patterns. + """ + + @property + @abc.abstractmethod + def is_impulse_response(self): + """ + Indicates whether the array contains coefficients for octave bands + (returns ``False``) or is a full-size impulse response (returns + ``True``). + """ + raise NotImplementedError + + @property + @abc.abstractmethod + def filter_len_ir(self): + """ + When ``is_impulse_response`` returns ``True``, this property returns the + lengths of the impulse responses returned. + All impulse responses are assumed to have the same length. + """ + raise NotImplementedError + + @abc.abstractmethod + def get_response( + self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True + ): + """ + Get response for provided angles. + + Parameters + ---------- + azimuth : array_like + Azimuth + colatitude : array_like, optional + Colatitude in degrees. Default is to be on XY plane. + magnitude : bool, optional + Whether to return magnitude of response. + frequency : float, optional + For which frequency to compute the response. + If the response is frequency independent, this parameter is ignored. + degrees : bool, optional + If ``True``, ``azimuth`` and ``colatitude`` are in degrees. + Otherwise, they are in radians. + + + Returns + ------- + resp : :py:class:`~numpy.ndarray` + Response at provided angles. + """ + raise NotImplementedError diff --git a/pyroomacoustics/directivities/direction.py b/pyroomacoustics/directivities/direction.py new file mode 100644 index 00000000..c0463ea1 --- /dev/null +++ b/pyroomacoustics/directivities/direction.py @@ -0,0 +1,230 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2020-2024 Robin Scheibler, Satvik Dixit, Eric Bezzam +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Using directivities makes sources and microphones having a different response +depending on the location of other objects. This means that their orientation +in 3D space matters. + +Some types of directivities such as ``CardioidFamily`` and derived classes +are defined only by a vector (i.e., direction). The response is then symmetric +around the axis defined by this vector. + +However, in general, not all directivities are symmetric in this way. +For the general case, the orientation can be defined by `Euler angles `_. +This is implemented in the class :py:class:`pyroomacoustics.direction.Rotation3D`. +""" +import numpy as np + + +class DirectionVector(object): + """ + Object for representing direction vectors in 3D, parameterized by an azimuth and colatitude + angle. + + 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 + + self._unit_v = np.array( + [ + np.cos(self._azimuth) * np.sin(self._colatitude), + np.sin(self._azimuth) * np.sin(self._colatitude), + np.cos(self._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 + + @property + def unit_vector(self): + """Direction vector in cartesian coordinates.""" + return self._unit_v + + +def _make_rot_matrix(a, axis): + """ + Compute the rotation matrix for a single rotation around the z-axis. + + .. reference: + https://en.wikipedia.org/wiki/Rotation_matrix + """ + if not 0 <= axis <= 2: + raise ValueError("Axis must be 0, 1, or 2.") + if axis == 2: + return np.array( + [ + [np.cos(a), -np.sin(a), 0], + [np.sin(a), np.cos(a), 0], + [0, 0, 1], + ] + ) + elif axis == 1: + return np.array( + [ + [np.cos(a), 0, np.sin(a)], + [0, 1, 0], + [-np.sin(a), 0, np.cos(a)], + ] + ) + else: + return np.array( + [ + [1, 0, 0], + [0, np.cos(a), -np.sin(a)], + [0, np.sin(a), np.cos(a)], + ] + ) + + +_axis_map = {"x": 0, "y": 1, "z": 2, 0: 0, 1: 1, 2: 2} + + +class Rotation3D: + """ + An object representing 3D rotations by their + `Euler angles `_. + + A rotation in 3D space can be fully described by 3 angles (i.e., the Euler angles). Each rotation is applied + around one of the three axes and there are 12 possible ways of pickinig the order or the rotations. + + This class can apply full or partial rotations to sets of points. + + The angles are provided as an array ``angles``. The axes of rotation for + the angles are provided in ``rot_order`` as a string of three characters + out of ``["x", "y", "z"]`` or a list of three integers out of ``[0, 1, + 2]``. Each axis can be repeated. To obtain full rotations, the same + axis should not be used twice in a row. + + By default, the angles are specified in degrees. This can be changed by setting + ``degrees=False``. In that case, the angles are assumed to be in radians. + + Parameters + ---------- + angles : array_like + An array containing between 0 and 3 angles. + rot_order : str of List[int] + The order of the rotations. The default is "zyx". + The order indicates around which axis the rotation is performed. + For example, "zyx" means that the rotation is first around the z-axis, + then the y-axis, and finally the x-axis. + degrees : bool + Whether the angles are in degrees (True) or radians (False). + """ + + def __init__(self, angles, rot_order="zyx", degrees=True): + angles = np.array(angles) + + if not 0 <= len(angles) <= 3: + raise ValueError("The number of angles must be between 0 and 3.") + + if degrees: + angles = np.radians(angles) + + self._angles = angles + + self._rot_order = [] + for ax in rot_order: + if ax not in _axis_map: + raise ValueError( + "Axis must be 'x', 'y', or 'z'. Alternatively, 0, 1, 2 can be used instead." + ) + self._rot_order.append(_axis_map[ax]) + + if len(angles) != len(self._rot_order): + raise ValueError( + f"The number of rotation angles ({len(angles)}) must the " + f"same as that of rotations axes ({len(self._rot_order)})." + ) + + self._rot_matrix = self._compute_matrix() + + def _compute_matrix(self): + """ + Compute the rotation matrix. + """ + mat = np.eye(3) + + for angle, axis in zip(self._angles, self._rot_order): + mat = np.dot(_make_rot_matrix(angle, axis), mat) + + return mat + + def rotate(self, points): + """ + Rotate a set of points. + + Parameters + ---------- + points : array_like (3, ...) + The points to rotate. The first dimension must be 3. + + Returns + ------- + rotated_points : np.ndarray + The rotated points. + """ + return np.dot(self._rot_matrix, points) + + def rotate_transpose(self, points): + """ + Transposed rotations of a set of points. + + Parameters + ---------- + points : array_like (3, ...) + The points to rotate. The last dimension must be 3. + + Returns + ------- + rotated_points : np.ndarray + The rotated points. + """ + return np.dot(self._rot_matrix.T, points) diff --git a/pyroomacoustics/directivities/integration.py b/pyroomacoustics/directivities/integration.py new file mode 100644 index 00000000..8e488f80 --- /dev/null +++ b/pyroomacoustics/directivities/integration.py @@ -0,0 +1,64 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2024 Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Provides a function to numerically integrate a function over the 3D unit-sphere (:math:`\mathbb{S}^2`). + +.. math:: + + \iint_{\mathbb{S}^2} f(\mathbf{x})\, d\mathbf{x} + +""" +import numpy as np +from scipy.spatial import SphericalVoronoi + +from pyroomacoustics.doa import fibonacci_spherical_sampling + + +def spherical_integral(func, n_points): + """ + Numerically integrate a function over the sphere. + + Parameters + ----------- + func: callable + The function to integrate. It should take an array of shape (3, n_points) + and return an array of shape (n_points,) + n_points: int + The number of points to use for integration + + Returns + ------- + value: np.ndarray + The value of the integral + """ + + points = fibonacci_spherical_sampling(n_points).T # shape (n_points, 3) + + # The weights are the areas of the voronoi cells + sv = SphericalVoronoi(points) + w_ = sv.calculate_areas() + + f = func(points.T) + + return np.sum(w_ * f) diff --git a/pyroomacoustics/directivities/interp.py b/pyroomacoustics/directivities/interp.py new file mode 100644 index 00000000..529b40a0 --- /dev/null +++ b/pyroomacoustics/directivities/interp.py @@ -0,0 +1,237 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2022-2024 Prerak Srivastava, Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +This module provides functions to interpolate impulse responses on a sphere. +The interpolation is done in the spherical harmonics domain. +""" +import numpy as np +import scipy +from scipy.spatial import SphericalVoronoi + +from ..doa import detect_regular_grid + + +def cal_sph_basis(azimuth, colatitude, degree): # theta_target,phi_target + """ + Calculate a spherical basis matrix + + Parameters + ----------- + azimuth: array_like + Azimuth of the spherical coordinates of the grid points + phi: array_like + Colatitude spherical coordinates of the grid points + degree:(int) + spherical harmonic degree + + Return + ------ + Ysh (np.array) shape: (no_of_nodes, (degree + 1)**2) + Spherical harmonics basis matrix + + """ + # build linear array of indices + # 0 + # -1 0 1 + # -2 -1 0 1 2 + # ... + ms = [] + ns = [] + for i in range(degree + 1): + for order in range(-i, i + 1): + ms.append(order) + ns.append(i) + m, n = np.array([ms]), np.array([ns]) + + # compute all the spherical harmonics at once + Ysh = scipy.special.sph_harm(m, n, azimuth[:, None], colatitude[:, None]) + + return Ysh + + +def _weighted_pinv(weights, Y, rcond=1e-2): + return np.linalg.pinv( + weights[:, None] * Y, rcond=rcond + ) # rcond is inverse of the condition number + + +def calculation_pinv_voronoi_cells(Ysh, colatitude, colatitude_grid, len_azimuth_grid): + """ + Weighted least square solution "Analysis and Synthesis of Sound-Radiation with Spherical Arrays: Franz Zotter Page 76" + + Calculation of pseudo inverse and voronoi cells for regular sampling in spherical coordinates + + Parameters + ----------- + Ysh: (np.ndarray) + Spherical harmonic basis matrix + colatitude: (np.ndarray) + The colatitudes of the measurements + colatitude_grid: (int) + the colatitudes of the grid lines + len_azimuth_grid: + The number of distinct azimuth values in the grid + + Returns: + ------------------------------- + Ysh_tilda_inv : (np.ndarray) + Weighted psuedo inverse of spherical harmonic basis matrix Ysh + w_ : (np.ndarray) + Weight on the original grid + """ + # compute the areas of the voronoi regions of the grid analytically + # assuming that the grid is regular in azimuth/colatitude + res = (colatitude_grid[:-1] + colatitude_grid[1:]) / 2 + res = np.insert(res, len(res), np.pi) + res = np.insert(res, 0, 0) + w = -np.diff(np.cos(res)) * (2 * np.pi / len_azimuth_grid) + w_dict = {t: ww for t, ww in zip(colatitude_grid, w)} + + # select the weights + w_ = np.array([w_dict[col] for col in colatitude]) + w_ /= 4 * np.pi # normalizing by unit sphere area + + return _weighted_pinv(w_, Ysh), w_ + + +def calculation_pinv_voronoi_cells_general(Ysh, points): + """ + Weighted least square solution "Analysis and Synthesis of Sound-Radiation with + Spherical Arrays: Franz Zotter Page 76" + + Calculation of pseudo inverse and voronoi cells for arbitrary sampling of the sphere. + + Parameters + ----------- + Ysh: (np.ndarray) + Spherical harmonic basis matrix + points: numpy.ndarray, (n_points, 3) + The sampling points on the sphere + + Returns: + ------------------------------- + Ysh_tilda_inv : (np.ndarray) + Weighted pseudo inverse of spherical harmonic basis matrix Ysh + w_ : (np.ndarray) + Weight on the original grid + """ + + # The weights are the areas of the voronoi cells + sv = SphericalVoronoi(points) + w_ = sv.calculate_areas() + w_ /= 4 * np.pi # normalizing by unit sphere area + + Ysh_tilda_inv = np.linalg.pinv( + w_[:, None] * Ysh, rcond=1e-2 + ) # rcond is inverse of the condition number + + return _weighted_pinv(w_, Ysh), w_ + + +def spherical_interpolation( + grid, + impulse_responses, + new_grid, + spherical_harmonics_order=12, + axis=-2, + nfft=None, +): + """ + Parameters + ---------- + grid: pyroomacoustics.doa.GridSphere + The grid of the measurements + impulse_responses: numpy.ndarray, (..., n_measurements, ..., n_samples) + The impulse responses to interpolate, the last axis is time and one other + axis should have dimension matching the length of the grid. By default, + it is assumed to be second from the end, but can be specified with the + `axis` argument. + new_grid: pyroomacoustics.doa.GridSphere + Grid of points at which to interpolate + spherical_harmonics_order: int + The order of spherical harmonics to use for interpolation + axis: int + The axis of the grid in the impulse responses array + nfft: int + The length of the FFT to use for the interpolation (default ``n_samples``) + """ + ir = np.swapaxes(impulse_responses, axis, -2) + if nfft is None: + nfft = ir.shape[-1] + + if len(grid) != ir.shape[-2]: + raise ValueError( + "The length of the grid should be the same as the number of impulse" + f"responses provide (grid={len(grid)}, impulse response={ir.shape[-2]})" + ) + + # Calculate spherical basis for the original grid + Ysh = cal_sph_basis(grid.azimuth, grid.colatitude, spherical_harmonics_order) + + # Calculate spherical basis for the target grid (fibonacci grid) + Ysh_fibo = cal_sph_basis( + new_grid.azimuth, + new_grid.colatitude, + spherical_harmonics_order, + ) + + # this will check if the points are on a regular grid. + # If they are, then the azimuths and colatitudes of the grid + # are returned + regular_grid = detect_regular_grid(grid.azimuth, grid.colatitude) + + # calculate pinv and voronoi cells for least square solution for the whole grid + if regular_grid is not None: + Ysh_tilda_inv, w_ = calculation_pinv_voronoi_cells( + Ysh, + grid.colatitude, + regular_grid.colatitude, + len(regular_grid.azimuth), + ) + else: + Ysh_tilda_inv, w_ = calculation_pinv_voronoi_cells_general( + Ysh, grid.cartesian.T + ) + + # Do the interpolation in the frequency domain + + # shape: (..., n_measurements, n_samples // 2 + 1) + tf = np.fft.rfft(ir, axis=-1, n=nfft) + + g_tilda = w_[:, None] * tf + + gamma_full_scale = np.matmul(Ysh_tilda_inv, g_tilda) + + interpolated_original_grid = np.fft.irfft( + np.matmul(Ysh, gamma_full_scale), n=nfft, axis=-1 + ) + interpolated_target_grid = np.fft.irfft( + np.matmul(Ysh_fibo, gamma_full_scale), n=nfft, axis=-1 + ) + + # restore the order of the axes + interpolated_target_grid = np.swapaxes(interpolated_target_grid, -2, axis) + interpolated_original_grid = np.swapaxes(interpolated_original_grid, -2, axis) + + return interpolated_target_grid, interpolated_original_grid diff --git a/pyroomacoustics/directivities/measured.py b/pyroomacoustics/directivities/measured.py new file mode 100644 index 00000000..2f95f04a --- /dev/null +++ b/pyroomacoustics/directivities/measured.py @@ -0,0 +1,458 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2022-2024 Prerak Srivastava, Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Source and microphone directivities can be measured in an anechoic chamber. +Such measurements result in a collection of impulse responses or transfer functions +each associated with a specific source and receiver (i.e., microphone) location. +The `SOFA file format `_ has been proposed +as a standard for the storage of such measurements. + +This sub-module offers a way to read such measurements from (SOFA) files and +use the measurement to obtain a more faithful simulation. + +The workhorse of this module is the class :py:class:`MeasuredDirectivityFile` +which reads the content of a file and standardize the data for futher use. +A single SOFA file can contain multiple measurements (for example corresponding +to different devices). The class provies a method to retrieve measurements +from individual sources and turn them into a py:class:`MeasuredDirectivity` object +that can be used to create a py:class:`pyroomacoustics.MicrophoneArray` object +with this directivity. + +Such measurements do not provide impulse responses for every possible impinging +direction. Instead, during simulation the impulse response closest to the +desired direction is used instead. To avoid sharp transitions, the py:class:`MeasuredDirectivityFile` +provides an interpolation method in the spherical harmonics domain. +This can be activated by providing an order for the interpolation, e.g, `interp_order=12`. + +Here is an example of loading a head-related transfer function and load +the directivities for left and right ears of a dummy head HRTF. + +.. code-block:: python + + from pyroomacoustics.directivities import MeasuredDirectivityFile, Rotation3D + + # the file reader object reads the file and optionally performs interpolation + # if the file contains multiple directivities, they are all read + hrtf = MeasuredDirectivityFile( + path="mit_kemar_normal_pinna.sofa", # SOFA file is in the database + fs=fs, + interp_order=12, + interp_n_points=1000, + ) + + # orientations can be provided as rotation matrices + orientation = Rotation3D([colatitude_deg, azimuth_deg], "yz", degrees=True) + + # we can then choose which directivities we want from the file + dir_left = hrtf.get_mic_directivity("left", orientation=orientation) + dir_right = hrtf.get_mic_directivity("right", orientation=orientation) + +""" +from pathlib import Path + +import numpy as np +from scipy.interpolate import griddata +from scipy.spatial import cKDTree + +from ..datasets import SOFADatabase +from ..doa import Grid, GridSphere, cart2spher, fibonacci_spherical_sampling, spher2cart +from ..utilities import requires_matplotlib +from .base import Directivity +from .direction import Rotation3D +from .interp import spherical_interpolation +from .sofa import open_sofa_file + + +class MeasuredDirectivity(Directivity): + """ + A class to store directivity patterns obtained by measurements. + + Parameters + ---------- + orientation: Rotation3D + A rotation to apply to the pattern + grid: doa.Grid + The grid of unit vectors where the measurements were taken + impulse_responses: np.ndarray, (n_grid, n_ir) + The impulse responses corresponding to the grid points + fs: int + The sampling frequency of the impulse responses + """ + + def __init__(self, orientation, grid, impulse_responses, fs): + if not isinstance(orientation, Rotation3D): + raise ValueError("Orientation must be a Rotation3D object") + + if not isinstance(grid, Grid): + raise ValueError("Grid must be a Grid object") + + self._original_grid = grid + self._irs = impulse_responses + self.fs = fs + + # set the initial orientation + self.set_orientation(orientation) + + @property + def is_impulse_response(self): + return True + + @property + def filter_len_ir(self): + """Length of the impulse response in samples""" + return self._irs.shape[-1] + + 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 ValueError("Orientation must be a Rotation3D object") + + self._orientation = orientation + + # rotate the grid and re-build the KD-tree + self._grid = GridSphere( + cartesian_points=self._orientation.rotate(self._original_grid.cartesian) + ) + # create the kd-tree + self._kdtree = cKDTree(self._grid.cartesian.T) + + def get_response( + self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True + ): + """ + Get response for provided angles and frequency. + + Parameters + ---------- + azimuth: np.ndarray, (n_points,) + The azimuth of the desired responses + colatitude: np.ndarray, (n_points,) + The colatitude of the desired responses + magnitude: bool + Ignored + frequency: np.ndarray, (n_freq,) + Ignored + degrees: bool + If ``True``, indicates that azimuth and colatitude are provided in degrees + """ + if degrees: + azimuth = np.radians(azimuth) + colatitude = np.radians(colatitude) + + cart = spher2cart(azimuth, colatitude) + + _, index = self._kdtree.query(cart.T) + return self._irs[index, :] + + @requires_matplotlib + def plot(self, freq_bin=0, n_grid=100, ax=None, depth=False, offset=None): + """ + Plot the directivity pattern at a given frequency. + + Parameters + ---------- + freq_bin: int + The frequency bin to plot + n_grid: int + The number of points to use for the interpolation grid + ax: matplotlib.axes.Axes, optional + The axes to plot on. If not provided, a new figure is created + depth: bool + If ``True``, directive response is both depicted by color and depth + of the surface. If ``False``, then only the color map denotes the + intensity. (default ``False``) + offset: float + An offset to apply to the directivity pattern + + Returns + ------- + ax: matplotlib.axes.Axes + The axes on which the directivity is plotted + """ + import matplotlib.pyplot as plt + from matplotlib import cm + + cart = self._grid.cartesian.T + length = np.abs(np.fft.rfft(self._irs, axis=-1)[:, freq_bin]) + + # regrid the data on a 2D grid + g = np.linspace(-1, 1, n_grid) + AZ, COL = np.meshgrid( + np.linspace(0, 2 * np.pi, n_grid), np.linspace(0, np.pi, n_grid // 2) + ) + # multiply by 0.99 to make sure the interpolation grid is inside the convex hull + # of the original points, otherwise griddata returns NaN + shrink_factor = 0.99 + while True: + X = np.cos(AZ) * np.sin(COL) * shrink_factor + Y = np.sin(AZ) * np.sin(COL) * shrink_factor + Z = np.cos(COL) * shrink_factor + grid = np.column_stack((X.flatten(), Y.flatten(), Z.flatten())) + interp_len = griddata(cart, length, grid, method="linear") + V = interp_len.reshape((n_grid // 2, n_grid)) + + # there may be some nan + if np.any(np.isnan(V)): + shrink_factor *= 0.99 + else: + break + + if ax is None: + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1, projection="3d") + + # Colour the plotted surface according to the sign of Y. + cmap = plt.cm.ScalarMappable(cmap=plt.get_cmap("coolwarm")) + cmap.set_clim(0, V.max()) + + if depth: + X *= V + Y *= V + Z *= V + + surf = ax.plot_surface( + X, Y, Z, facecolors=cmap.to_rgba(V), linewidth=0, antialiased=False + ) + + return ax + + +class MeasuredDirectivityFile: + """ + This class reads measured directivities from a + `SOFA`_ format file. + Optionally, it can perform interpolation of the impulse responses onto a finer grid. + The interpolation is done in the spherical harmonics domain. + + + Parameters + -------------- + path : (string) + Path towards the specific DIRPAT file + fs: (int) + The desired sampling frequency. If the impulse responses were stored at + a different sampling frequency, they are resampled at ``fs``. + interp_order: (int) + The order of spherical harmonics to use for interpolation. + If ``None`` interpolation is not used. + interp_n_points: (int) + Number of points for the interpolation grid. The interpolation grid is a + Fibonnaci pseudo-uniform sampling of the sphere. + file_reader_callback: (callable) + A callback function that reads the SOFA file and returns the impulse responses + The signature should be the same as the function `open_sofa_file` + mic_labels: (list of strings) + List of labels for the microphones. If not provided, the labels are simply the + indices of the microphones in the array + source_labels: (list of strings) + List of labels for the sources. If not provided, the labels are simply the + indices of the measurements in the array + """ + + def __init__( + self, + path, + fs=None, + interp_order=None, + interp_n_points=1000, + file_reader_callback=None, + mic_labels=None, + source_labels=None, + ): + self.path = Path(path) + + if file_reader_callback is None: + # default reader is for SOFA files + file_reader_callback = open_sofa_file + + ( + self.impulse_responses, # (n_sources, n_mics, taps) + self.fs, + self.source_locs, # (3, n_sources), spherical coordinates + self.mic_locs, # (3, n_mics), cartesian coordinates + src_labels_file, + mic_labels_file, + ) = file_reader_callback( + path=self.path, + fs=fs, + ) + + if mic_labels is None: + self.mic_labels = mic_labels_file + else: + if len(mic_labels) != self.mic_locs.shape[1]: + breakpoint() + raise ValueError( + f"Number of labels provided ({len(mic_labels)}) does not match the " + f"number of microphones ({self.mic_locs.shape[1]})" + ) + self.mic_labels = mic_labels + + if source_labels is None: + self.source_labels = src_labels_file + else: + if len(source_labels) != self.source_locs.shape[1]: + raise ValueError( + f"Number of labels provided ({len(source_labels)}) does not match " + f"the number of sources ({self.source_locs.shape[1]})" + ) + self.source_labels = source_labels + + self.interp_order = interp_order + self.interp_n_points = interp_n_points + + self._ir_interp_cache = {} + if interp_order is not None: + self.interp_grid = GridSphere( + cartesian_points=fibonacci_spherical_sampling(n_points=interp_n_points) + ) + else: + self.interp_grid = None + + def _interpolate(self, type, mid, grid, impulse_responses): + if self.interp_order is None: + return grid, impulse_responses + + label = f"{type}_{mid}" + + if label not in self._ir_interp_cache: + self._ir_interp_cache[label], _ = spherical_interpolation( + grid, + impulse_responses, + self.interp_grid, + spherical_harmonics_order=self.interp_order, + axis=-2, + ) + + return self.interp_grid, self._ir_interp_cache[label] + + def _get_measurement_index(self, meas_id, labels): + if isinstance(meas_id, int): + return meas_id + elif labels is not None: + idx = labels.index(meas_id) + if idx >= 0: + return idx + else: + raise KeyError(f"Measurement id {meas_id} not found") + + raise ValueError(f"Measurement id {meas_id} not found") + + def get_mic_position(self, measurement_id): + """ + Get the position of source with id `measurement_id` + + Parameters + ---------- + measurement_id: int or str + The id of the source + """ + mid = self._get_measurement_index(measurement_id, self.mic_labels) + + if not (0 <= mid < self.mic_locs.shape[1]): + raise ValueError(f"Microphone id {mid} not found") + + return self.mic_locs[:, mid] + + def get_source_position(self, measurement_id): + """ + Get the position of source with id `measurement_id` + + Parameters + ---------- + measurement_id: int or str + The id of the source + """ + mid = self._get_measurement_index(measurement_id, self.source_labels) + + if not (0 <= mid < self.source_locs.shape[1]): + raise ValueError(f"Source id {mid} not found") + + # convert to cartesian since the sources are stored by + # default in spherical coordinates + pos = spher2cart(*self.source_locs[:, mid]) + + return pos + + def get_mic_directivity(self, measurement_id, orientation): + """ + Get a directivity for a microphone + + Parameters + ---------- + measurement_id: int or str + The id of the microphone + orientation: Rotation3D + The orientation of the directivity pattern + """ + mid = self._get_measurement_index(measurement_id, self.mic_labels) + + if not (0 <= mid < self.mic_locs.shape[1]): + raise ValueError(f"Microphone id {mid} not found") + + # select the measurements corresponding to the mic id + ir = self.impulse_responses[:, mid, :] + src_grid = GridSphere(spherical_points=self.source_locs[:2]) + + # interpolate the IR + grid, ir = self._interpolate("mic", mid, src_grid, ir) + + dir_obj = MeasuredDirectivity(orientation, grid, ir, self.fs) + return dir_obj + + def get_source_directivity(self, measurement_id, orientation): + """ + Get a directivity for a source + + Parameters + ---------- + measurement_id: int or str + The id of the source + orientation: Rotation3D + The orientation of the directivity pattern + """ + mid = self._get_measurement_index(measurement_id, self.source_labels) + + if not (0 <= mid < self.source_locs.shape[1]): + raise ValueError(f"Source id {mid} not found") + + # select the measurements corresponding to the mic id + ir = self.impulse_responses[mid, :, :] + + # here we need to swap the coordinate types + mic_pos = np.array(cart2spher(self.mic_locs)) + mic_grid = GridSphere(spherical_points=mic_pos[:2]) + + # interpolate the IR + grid, ir = self._interpolate("source", mid, mic_grid, ir) + + dir_obj = MeasuredDirectivity(orientation, grid, ir, self.fs) + return dir_obj diff --git a/pyroomacoustics/directivities/sofa.py b/pyroomacoustics/directivities/sofa.py new file mode 100644 index 00000000..4ecff49e --- /dev/null +++ b/pyroomacoustics/directivities/sofa.py @@ -0,0 +1,332 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2022-2024 Prerak Srivastava, Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +`SOFA `_ is a very flexible file format for storing +direction impulse responses. This module provides a function to read SOFA files and +extract the impulse responses in a format that can be used for simulation. +""" +from pathlib import Path + +import numpy as np + +try: + import sofa + + has_sofa = True +except ImportError: + has_sofa = False + + +from ..datasets.sofa import get_sofa_db, is_dirpat, resolve_sofa_path +from ..doa import cart2spher, spher2cart +from ..utilities import resample + + +def open_sofa_file(path, fs=16000): + """ + Open a SOFA file and read the impulse responses + + Parameters + ---------- + path: str or Path + Path to the SOFA file + fs: int, optional + The desired sampling frequency. If the impulse responses were stored at + a different sampling frequency, they are resampled at ``fs``. + + Returns + ------- + ir: np.ndarray (n_sources, n_mics, taps) + The impulse responses + fs: int + The sampling frequency of the impulse responses + source_dir: np.ndarray (3, n_sources) + The direction of the sources in spherical coordinates + rec_loc: np.ndarray (3, n_mics) + The location of the receivers in cartesian coordinates + source_labels: List[str] + If available, a list of human readable labels for the sources is + returned. Otherwise, ``None`` is returned + mic_labels: List[str] + If available, a list of human readable labels for the microphones is + returned. Otherwise, ``None`` is returned + """ + # Memo for notation of SOFA dimensions + # From: https://www.sofaconventions.org/mediawiki/index.php/SOFA_conventions#AnchorDimensions + # M number of measurements integer >0 + # R number of receivers or harmonic coefficients describing receivers integer >0 + # E number of emitters or harmonic coefficients describing emitters integer >0 + # N number of data samples describing one measurement integer >0 + # S number of characters in a string integer ≥0 + # I singleton dimension, constant always 1 + # C coordinate triplet, constant always 3 + + # Open DirPat database + if not has_sofa: + raise ValueError( + "The package 'python-sofa' needs to be installed to call this function. Install by doing `pip install python-sofa`" + ) + + sofa_db = get_sofa_db() + + path = resolve_sofa_path(path) + + file_sofa = sofa.Database.open(path) + + # we have a special case for DIRPAT files because they need surgery + if is_dirpat(path): + # special case for DIRPAT files + # this is needed because the DIRPAT files are not following the SOFA conventions + # and require some custom processing + ir, fs, src_loc, mic_loc = _read_dirpat(file_sofa, path.stem, fs) + + else: + # readers for SOFA files following the conventions + + conv_name = file_sofa.convention.name + + if conv_name == "SimpleFreeFieldHRIR": + ir, fs, src_loc, mic_loc = _read_simple_free_field_hrir(file_sofa, fs) + + elif conv_name == "GeneralFIR": + ir, fs, src_loc, mic_loc = _read_general_fir(file_sofa, fs) + + else: + raise NotImplementedError(f"SOFA convention {conv_name} not implemented") + + # assign the source labels if known + source_labels = None + mic_labels = None + + file_id = path.stem + if file_id in sofa_db: + info = sofa_db[file_id] + if info.type == "sources": + source_labels = info.contains + mic_labels = None + if source_labels is not None and len(source_labels) != ir.shape[0]: + raise ValueError( + f"Number of labels ({len(source_labels)}) does not match the " + f"number of sources ({ir.shape[0]})" + ) + elif info.type == "microphones": + source_labels = None + mic_labels = info.contains + if mic_labels is not None and len(mic_labels) != ir.shape[1]: + raise ValueError( + f"Number of labels ({len(mic_labels)}) does not match the " + f"number of microphones ({ir.shape[1]})" + ) + + return ir, fs, src_loc, mic_loc, source_labels, mic_labels + + +def _parse_locations(sofa_pos, target_format): + """ + Reads and normalize a position stored in a SOFA file + + Parameters + ---------- + sofa_pos: + SOFA position object + target_format: + One of 'spherical' or 'cartesian'. For 'spherical', the + angles are always in radians + + Returns + ------- + A numpy array in the correct format + """ + + if target_format not in ("spherical", "cartesian"): + raise ValueError("Target format should be 'spherical' or 'cartesian'") + + # SOFA dimensions + dim = sofa_pos.dimensions() + + # source positions + pos = sofa_pos.get_values() + + if len(dim) == 3 and dim[-1] == "I": + pos = pos[..., 0] + dim = dim[:-1] + + # get units + pos_units = sofa_pos.Units + if "," in pos_units: + pos_units = pos_units.split(",") + pos_units = [p.strip() for p in pos_units] + else: + pos_units = [pos_units] * pos.shape[1] + + pos_type = sofa_pos.Type + + if pos_type == "cartesian": + if any([p != "metre" for p in pos_units]): + raise ValueError(f"Found unit '{pos_units}' in SOFA file") + + pos = pos.T + + if target_format == "spherical": + return np.array(cart2spher(pos)) + else: + return pos + + elif pos_type == "spherical": + azimuth = pos[:, 0] if pos_units[0] != "degree" else np.deg2rad(pos[:, 0]) + colatitude = pos[:, 1] if pos_units[0] != "degree" else np.deg2rad(pos[:, 1]) + distance = pos[:, 2] + + if np.any(colatitude < 0.0): + # it looks like the data is using elevation format + colatitude = np.pi / 2.0 - colatitude + + if target_format == "cartesian": + return spher2cart(azimuth, colatitude, distance) + else: + return np.array([azimuth, colatitude, distance]) + + else: + raise NotImplementedError(f"{pos_type} not implemented") + + +def _read_simple_free_field_hrir(file_sofa, fs): + """ + Reads the HRIRs stored in a SOFA file with the SimpleFreeFieldHRIR convention + + Parameters + ---------- + file_sofa: SOFA object + Path to the SOFA file + fs: int + The desired sampling frequency. If the impulse responses were stored at + a different sampling frequency, they are resampled at ``fs`` + + Returns + ------- + ir: np.ndarray + The impulse responses in format ``(n_sources, n_mics, taps)`` + source_dir: np.ndarray + The direction of the sources in spherical coordinates + ``(3, n_sources)`` where the first row is azimuth and the second is colatitude + and the third is distance + rec_loc: np.ndarray + The location of the receivers in cartesian coordinates with respect to the + origin of the SOFA file + fs: int + The sampling frequency of the impulse responses + """ + # read the mesurements (source direction, receiver location, taps) + msr = file_sofa.Data.IR.get_values() + + # downsample the fir filter. + fs_file = file_sofa.Data.SamplingRate.get_values()[0] + if fs is None: + fs = fs_file + else: + msr = resample(msr, fs_file, fs) + + # Source positions + source_loc = _parse_locations(file_sofa.Source.Position, target_format="spherical") + + # Receivers locations (i.e., "ears" for HRIR) + rec_loc = _parse_locations(file_sofa.Receiver.Position, target_format="cartesian") + + return msr, fs, source_loc, rec_loc + + +def _read_general_fir(file_sofa, fs): + # read the mesurements (source direction, receiver location, taps) + msr = file_sofa.Data.IR.get_values() + + # downsample the fir filter. + fs_file = file_sofa.Data.SamplingRate.get_values()[0] + if fs is None: + fs = fs_file + else: + msr = resample(msr, fs_file, fs) + + # Source positions: (azimuth, colatitude, distance) + source_loc = _parse_locations(file_sofa.Source.Position, target_format="spherical") + + # Receivers locations (i.e., "ears" for HRIR) + rec_loc = _parse_locations(file_sofa.Receiver.Position, target_format="cartesian") + + return msr, fs, source_loc, rec_loc + + +def _read_dirpat(file_sofa, filename, fs=None): + sofa_db = get_sofa_db() + + if filename not in sofa_db: + raise ValueError(f"DIRPAT file {filename} not found in the SOFA database") + + # read the mesurements + msr = file_sofa.Data.IR.get_values() # (n_sources, n_mics, taps) + + # downsample the fir filter. + fs_file = file_sofa.Data.SamplingRate.get_values()[0] + if fs is None: + fs = fs_file + else: + msr = resample(msr, fs_file, fs) + + # Receiver positions + mic_pos = file_sofa.Receiver.Position.get_values() # (3, n_mics) + mic_pos_units = file_sofa.Receiver.Position.Units.split(",") + + # Source positions + src_pos = file_sofa.Source.Position.get_values() + src_pos_units = file_sofa.Source.Position.Units.split(",") + + # There is a bug in the DIRPAT measurement files where the array of + # measurement locations were not flattened correctly + src_pos_units[0:1] = "radian" + if filename == "LSPs_HATS_GuitarCabinets_Akustikmessplatz": + # this is a source file + mic_pos_RS = np.reshape(mic_pos, [36, -1, 3]) + mic_pos = np.swapaxes(mic_pos_RS, 0, 1).reshape([mic_pos.shape[0], -1]) + + if np.any(mic_pos[:, 1] < 0.0): + # it looks like the data is using elevation format + mic_pos[:, 1] = np.pi / 2.0 - mic_pos[:, 1] + + # by convention, we keep the microphone locations in cartesian coordinates + mic_pos = spher2cart(*mic_pos.T).T + + # create source locations, they are all at the center + src_pos = np.zeros((msr.shape[0], 3)) + + else: + src_pos_RS = np.reshape(src_pos, [30, -1, 3]) + src_pos = np.swapaxes(src_pos_RS, 0, 1).reshape([src_pos.shape[0], -1]) + + if np.any(src_pos[:, 1] < 0.0): + # it looks like the data is using elevation format + src_pos[:, 1] = np.pi / 2.0 - src_pos[:, 1] + + # create fake microphone locations, they are all at the center + mic_pos = np.zeros((msr.shape[1], 3)) + + return msr, fs, src_pos.T, mic_pos.T diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_False-cardioid.npy new file mode 100644 index 00000000..e4932be8 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_False-oneside.npy new file mode 100644 index 00000000..f05f2b99 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_True-cardioid.npy new file mode 100644 index 00000000..3608f08d Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_True-oneside.npy new file mode 100644 index 00000000..7bb7ee17 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414A-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_False-cardioid.npy new file mode 100644 index 00000000..60cb159c Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_False-oneside.npy new file mode 100644 index 00000000..c77b16de Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_True-cardioid.npy new file mode 100644 index 00000000..0c20cd81 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_True-oneside.npy new file mode 100644 index 00000000..5ab0a3d1 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414K-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_False-cardioid.npy new file mode 100644 index 00000000..ac79fcbb Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_False-oneside.npy new file mode 100644 index 00000000..93cf99cc Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_True-cardioid.npy new file mode 100644 index 00000000..79b9fcd1 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_True-oneside.npy new file mode 100644 index 00000000..ec96d3fe Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414N-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_False-cardioid.npy new file mode 100644 index 00000000..4947e611 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_False-oneside.npy new file mode 100644 index 00000000..0d0273e2 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_True-cardioid.npy new file mode 100644 index 00000000..179840a1 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_True-oneside.npy new file mode 100644 index 00000000..4ca6c646 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c414S-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_False-cardioid.npy new file mode 100644 index 00000000..0dca3675 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_False-oneside.npy new file mode 100644 index 00000000..030d7a7d Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_True-cardioid.npy new file mode 100644 index 00000000..f64a5ccc Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_True-oneside.npy new file mode 100644 index 00000000..37f7010c Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/AKG_c480_c414_CUBE-AKG_c480-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_False-cardioid.npy new file mode 100644 index 00000000..556a99f5 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_False-oneside.npy new file mode 100644 index 00000000..1750dc60 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_True-cardioid.npy new file mode 100644 index 00000000..4574617f Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_True-oneside.npy new file mode 100644 index 00000000..97e63320 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_0-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_False-cardioid.npy new file mode 100644 index 00000000..d87fceee Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_False-oneside.npy new file mode 100644 index 00000000..8c518b6e Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_True-cardioid.npy new file mode 100644 index 00000000..410ea359 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_True-oneside.npy new file mode 100644 index 00000000..27b3038c Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/EM32_Directivity-EM_32_31-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-AKG_c480_c414_CUBE-AKG_c480-minphase_False-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-AKG_c480_c414_CUBE-AKG_c480-minphase_False-twosides.npy new file mode 100644 index 00000000..64f76fe8 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-AKG_c480_c414_CUBE-AKG_c480-minphase_False-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-AKG_c480_c414_CUBE-AKG_c480-minphase_True-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-AKG_c480_c414_CUBE-AKG_c480-minphase_True-twosides.npy new file mode 100644 index 00000000..0cffb067 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-AKG_c480_c414_CUBE-AKG_c480-minphase_True-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-EM32_Directivity-EM_32_31-minphase_False-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-EM32_Directivity-EM_32_31-minphase_False-twosides.npy new file mode 100644 index 00000000..30de291d Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-EM32_Directivity-EM_32_31-minphase_False-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-EM32_Directivity-EM_32_31-minphase_True-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-EM32_Directivity-EM_32_31-minphase_True-twosides.npy new file mode 100644 index 00000000..224395ff Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-EM32_Directivity-EM_32_31-minphase_True-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_False-cardioid.npy new file mode 100644 index 00000000..eb82026e Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_False-oneside.npy new file mode 100644 index 00000000..37a934aa Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_True-cardioid.npy new file mode 100644 index 00000000..81fb085e Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_True-oneside.npy new file mode 100644 index 00000000..a1f537ff Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Genelec_8020-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-AKG_c480_c414_CUBE-AKG_c414K-minphase_False-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-AKG_c480_c414_CUBE-AKG_c414K-minphase_False-twosides.npy new file mode 100644 index 00000000..ba957bfe Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-AKG_c480_c414_CUBE-AKG_c414K-minphase_False-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-AKG_c480_c414_CUBE-AKG_c414K-minphase_True-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-AKG_c480_c414_CUBE-AKG_c414K-minphase_True-twosides.npy new file mode 100644 index 00000000..864c67cb Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-AKG_c480_c414_CUBE-AKG_c414K-minphase_True-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-EM32_Directivity-EM_32_0-minphase_False-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-EM32_Directivity-EM_32_0-minphase_False-twosides.npy new file mode 100644 index 00000000..2a606f95 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-EM32_Directivity-EM_32_0-minphase_False-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-EM32_Directivity-EM_32_0-minphase_True-twosides.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-EM32_Directivity-EM_32_0-minphase_True-twosides.npy new file mode 100644 index 00000000..d22bfec2 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-EM32_Directivity-EM_32_0-minphase_True-twosides.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_False-cardioid.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_False-cardioid.npy new file mode 100644 index 00000000..5d868346 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_False-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_False-oneside.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_False-oneside.npy new file mode 100644 index 00000000..9a0e038f Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_False-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_True-cardioid.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_True-cardioid.npy new file mode 100644 index 00000000..1a0d1f10 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_True-cardioid.npy differ diff --git a/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_True-oneside.npy b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_True-oneside.npy new file mode 100644 index 00000000..d3ffd239 Binary files /dev/null and b/pyroomacoustics/directivities/tests/data/LSPs_HATS_GuitarCabinets_Akustikmessplatz-Vibrolux_2x10inch-minphase_True-oneside.npy differ diff --git a/pyroomacoustics/directivities/tests/test_cardioid_energy.py b/pyroomacoustics/directivities/tests/test_cardioid_energy.py new file mode 100644 index 00000000..4a91bca8 --- /dev/null +++ b/pyroomacoustics/directivities/tests/test_cardioid_energy.py @@ -0,0 +1,50 @@ +""" +This test compares the theoretical value of the energy integral of the cardioid function +with the numerical value obtained by integrating the squared cardioid function over the +unit sphere. The test is performed for different values of the parameter p and the gain. +""" + +import numpy as np +import pytest + +from pyroomacoustics.directivities.analytic import ( + CardioidFamily, + cardioid_energy, + cardioid_func, +) +from pyroomacoustics.directivities.direction import DirectionVector +from pyroomacoustics.directivities.integration import spherical_integral +from pyroomacoustics.doa import cart2spher + +PARAMETERS = [(p, G) for p in [0.0, 0.25, 0.5, 0.75, 1.0] for G in [1.0, 0.5, 2.0]] + + +@pytest.mark.parametrize("p,gain", PARAMETERS) +def test_cardioid_func_energy(p, gain): + n_points = 10000 + direction = np.array([0.0, 0.0, 1.0]) + + def func(points): + return cardioid_func(points, direction, p, gain=gain) ** 2 + + num = spherical_integral(func, n_points) + thy = cardioid_energy(p, gain=gain) + assert abs(num - thy) < 1e-4 + + +@pytest.mark.parametrize("p,gain", PARAMETERS) +def test_cardioid_family_energy(p, gain): + n_points = 10000 + direction = np.array([0.0, 0.0, 1.0]) + + e3 = DirectionVector(0.0, 0.0) # vector pointing up + + card_obj = CardioidFamily(e3, p, gain=gain) + + def func(points): + az, co, _ = cart2spher(points) + return card_obj.get_response(az, co, degrees=False) ** 2 + + num = spherical_integral(func, n_points) + thy = cardioid_energy(p, gain=gain) + assert abs(num - thy) < 1e-4 diff --git a/pyroomacoustics/directivities/tests/test_sofa_directivities.py b/pyroomacoustics/directivities/tests/test_sofa_directivities.py new file mode 100644 index 00000000..be92561c --- /dev/null +++ b/pyroomacoustics/directivities/tests/test_sofa_directivities.py @@ -0,0 +1,694 @@ +""" +This is a regression test for the SOFA source and receiver measured directivity patterns + +The tests compare the output of the simulation with some pre-generated samples. + +To generate the samples run this file: `python ./test_sofa_directivities.py` +""" + +import argparse +import os +import warnings +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pytest + +import pyroomacoustics as pra +from pyroomacoustics.datasets.sofa import ( + DEFAULT_SOFA_PATH, + download_sofa_files, + get_sofa_db_info, +) +from pyroomacoustics.directivities import ( + CardioidFamily, + DirectionVector, + FigureEight, + MeasuredDirectivityFile, + Rotation3D, +) +from pyroomacoustics.directivities.interp import ( + calculation_pinv_voronoi_cells, + calculation_pinv_voronoi_cells_general, +) +from pyroomacoustics.doa import GridSphere + +sofa_info = get_sofa_db_info() +supported_sofa = [name for name, info in sofa_info.items() if info["supported"] == True] +save_plot = False +interp_order = 12 + +TEST_DATA = Path(__file__).parent / "data" + +# the processing delay due to the band-pass filters was removed in +# after the test files were created +# we need to subtract this delay from the reference signal +ref_delay = 0 + +# tolerances for the regression tests +atol = 5e-3 +rtol = 3e-2 + +room_dim = [6, 6, 2.4] + +all_materials = { + "east": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "west": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "north": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "south": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "ceiling": pra.Material( + energy_absorption={ + "coeffs": [0.02, 0.02, 0.03, 0.03, 0.04, 0.05], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), + "floor": pra.Material( + energy_absorption={ + "coeffs": [0.11, 0.14, 0.37, 0.43, 0.27, 0.25], + "center_freqs": [125, 250, 500, 1000, 2000, 4000], + }, + scattering=0.54, + ), +} + +""" +all_materials = pra.Material(0.1) +all_materials = pra.Material( + energy_absorption={ + "coeffs": [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1], + "center_freqs": [125, 250, 500, 1000, 2000, 4000, 8000], + }, +) +""" + + +def test_dirpat_download(): + files = download_sofa_files(verbose=True, no_fail=False) + for file in files: + assert file.exists() + + +SOFA_ONE_SIDE_PARAMETERS = [ + ("AKG_c480", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414K", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414N", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414S", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414A", "AKG_c480_c414_CUBE", False, False, save_plot), + ("EM_32_0", "EM32_Directivity", False, False, save_plot), + ("EM_32_31", "EM32_Directivity", False, False, save_plot), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + False, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + False, + False, + save_plot, + ), + ("AKG_c480", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414K", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414N", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414S", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414A", "AKG_c480_c414_CUBE", True, False, save_plot), + ("EM_32_0", "EM32_Directivity", True, False, save_plot), + ("EM_32_31", "EM32_Directivity", True, False, save_plot), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + True, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + True, + False, + save_plot, + ), +] + + +@pytest.mark.parametrize( + "pattern_id,sofa_file_name,min_phase,save_flag,plot_flag", SOFA_ONE_SIDE_PARAMETERS +) +def test_sofa_one_side(pattern_id, sofa_file_name, min_phase, save_flag, plot_flag): + """ + Tests with only microphone *or* source from a SOFA file + """ + + if min_phase: + pra.constants.set("octave_bands_n_fft", 128) + pra.constants.set("octave_bands_keep_dc", True) + else: + pra.constants.set("octave_bands_n_fft", 512) + pra.constants.set("octave_bands_keep_dc", False) + + # make sure all the files are present + # in case tests are run out of order + download_sofa_files(overwrite=False, verbose=False) + + # create room + room = pra.ShoeBox( + room_dim, + fs=16000, + max_order=2, + materials=all_materials, + air_absorption=True, + ray_tracing=False, + min_phase=min_phase, + ) + + # define source with figure_eight directivity + dir_factory = MeasuredDirectivityFile( + path=sofa_file_name, + fs=room.fs, + interp_order=interp_order, + interp_n_points=1000, + ) + + if sofa_info[sofa_file_name]["type"] == "sources": + directivity = dir_factory.get_source_directivity( + pattern_id, + # orientation=DirectionVector(azimuth=0, colatitude=0, degrees=False), + orientation=Rotation3D([0.0, 0.0], rot_order="yz"), + ) + else: + directivity = dir_factory.get_mic_directivity( + pattern_id, + # orientation=DirectionVector(azimuth=0, colatitude=0, degrees=False), + orientation=Rotation3D([0.0, 0.0], rot_order="yz"), + ) + + # add source with figure_eight directivity + if sofa_info[sofa_file_name]["type"] == "microphones": + directivity_SRC = None + directivity_MIC = dir_factory.get_mic_directivity( + pattern_id, + # orientation=DirectionVector(azimuth=0, colatitude=0, degrees=False), + orientation=Rotation3D([0.0, 0.0], rot_order="yz"), + ) + directivity = directivity_MIC + elif sofa_info[sofa_file_name]["type"] == "sources": + directivity_SRC = dir_factory.get_source_directivity( + pattern_id, + # orientation=DirectionVector(azimuth=0, colatitude=0, degrees=False), + orientation=Rotation3D([0.0, 0.0], rot_order="yz"), + ) + directivity_MIC = None + directivity = directivity_SRC + else: + raise ValueError("unknown pattern type") + + room.add_source([1.52, 0.883, 1.044], directivity=directivity_SRC) + + # add microphone in its null + room.add_microphone([2.31, 1.65, 1.163], directivity=directivity_MIC) + + # Check set different orientation after intailization of the DIRPATRir class + # directivity.set_orientation(DirectionVector(0.0, 0.0)) + directivity.set_orientation(Rotation3D([0.0, 0.0], rot_order="yz")) + + room.compute_rir() + + rir_1_0 = room.rir[0][0] + + filename = ( + "-".join([sofa_file_name.split(".")[0], pattern_id]) + + f"-minphase_{min_phase}-oneside.npy" + ) + test_file_path = TEST_DATA / filename + if save_flag: + TEST_DATA.mkdir(exist_ok=True, parents=True) + np.save(test_file_path, rir_1_0) + elif test_file_path.exists(): + reference_data = np.load(test_file_path) + reference_data = reference_data[ref_delay : ref_delay + rir_1_0.shape[0]] + rir_1_0 = rir_1_0[: reference_data.shape[0]] + print("Max diff.:", abs(reference_data - rir_1_0).max()) + print( + "Rel diff.:", + abs(reference_data - rir_1_0).max() / abs(reference_data).max(), + ) + if plot_flag: + fig, ax = plt.subplots(1, 1) + ax.plot(rir_1_0, label="test") + ax.plot(reference_data, label="ref") + ax.legend() + fig.savefig(test_file_path.with_suffix(".pdf")) + plt.close(fig) + else: + assert np.allclose(reference_data, rir_1_0, atol=atol, rtol=rtol) + else: + warnings.warn("Did not find the reference data. Output was not checked.") + + +SOFA_TWO_SIDES_PARAMETERS = [ + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "AKG_c480", + "AKG_c480_c414_CUBE", + False, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "AKG_c414K", + "AKG_c480_c414_CUBE", + False, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "EM_32_0", + "EM32_Directivity", + False, + False, + save_plot, + ), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "EM_32_31", + "EM32_Directivity", + False, + False, + save_plot, + ), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "AKG_c480", + "AKG_c480_c414_CUBE", + True, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "AKG_c414K", + "AKG_c480_c414_CUBE", + True, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "EM_32_0", + "EM32_Directivity", + True, + False, + save_plot, + ), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + "EM_32_31", + "EM32_Directivity", + True, + False, + save_plot, + ), +] + + +@pytest.mark.parametrize( + "src_pattern_id, src_sofa_file_name, mic_pattern_id, " + "mic_sofa_file_name, min_phase, save_flag, plot_flag", + SOFA_TWO_SIDES_PARAMETERS, +) +def test_sofa_two_sides( + src_pattern_id, + src_sofa_file_name, + mic_pattern_id, + mic_sofa_file_name, + min_phase, + save_flag, + plot_flag, +): + """ + Tests with only microphone *or* source from a SOFA file + """ + + if min_phase: + pra.constants.set("octave_bands_n_fft", 128) + pra.constants.set("octave_bands_keep_dc", True) + else: + pra.constants.set("octave_bands_n_fft", 512) + pra.constants.set("octave_bands_keep_dc", False) + + # make sure all the files are present + # in case tests are run out of order + download_sofa_files(overwrite=False, verbose=False) + + # create room + room = pra.ShoeBox( + room_dim, + fs=16000, + max_order=2, + materials=all_materials, + air_absorption=True, + ray_tracing=False, + min_phase=min_phase, + ) + + src_factory = MeasuredDirectivityFile( + path=src_sofa_file_name, + fs=room.fs, + interp_order=interp_order, + ) + src_directivity = src_factory.get_source_directivity( + src_pattern_id, + Rotation3D( + [0, 0], rot_order="yz" + ), # DirectionVector(azimuth=0, colatitude=0, degrees=True) + ) + + mic_factory = MeasuredDirectivityFile( + path=mic_sofa_file_name, + fs=room.fs, + interp_order=interp_order, + ) + mic_directivity = mic_factory.get_mic_directivity( + mic_pattern_id, + Rotation3D( + [0, 0], rot_order="yz" + ), # DirectionVector(azimuth=0, colatitude=0, degrees=True) + ) + + room.add_source([1.52, 0.883, 1.044], directivity=src_directivity) + + # add microphone in its null + room.add_microphone([2.31, 1.65, 1.163], directivity=mic_directivity) + + # Check set different orientation after intailization of the DIRPATRir class + # mic_directivity.set_orientation(DirectionVector(np.radians(np.pi), 0, degrees=True)) + mic_directivity.set_orientation( + Rotation3D([0.0, np.radians(np.pi)], rot_order="yz", degrees=True) + ) + # src_directivity.set_orientation( + # DirectionVector(0, np.radians(np.pi / 2.0), degrees=True) + # ) + src_directivity.set_orientation( + Rotation3D([np.radians(np.pi / 2.0), 0.0], rot_order="yz", degrees=True) + ) + + room.compute_rir() + + rir_1_0 = room.rir[0][0] + + filename = ( + "-".join( + [ + src_sofa_file_name.split(".")[0], + src_pattern_id, + mic_sofa_file_name.split(".")[0], + mic_pattern_id, + ] + ) + + f"-minphase_{min_phase}-twosides.npy" + ) + test_file_path = TEST_DATA / filename + if save_flag: + TEST_DATA.mkdir(exist_ok=True, parents=True) + np.save(test_file_path, rir_1_0) + elif test_file_path.exists(): + reference_data = np.load(test_file_path) + reference_data = reference_data[ref_delay : ref_delay + rir_1_0.shape[0]] + + if plot_flag: + fig, ax = plt.subplots(1, 1) + ax.plot(rir_1_0, label="test") + ax.plot(reference_data, label="ref") + ax.legend() + fig.savefig(test_file_path.with_suffix(".pdf")) + plt.close(fig) + + print("Max diff.:", abs(reference_data - rir_1_0).max()) + print( + "Rel diff.:", + abs(reference_data - rir_1_0).max() / abs(reference_data).max(), + ) + + if not plot_flag: + assert np.allclose(reference_data, rir_1_0, atol=atol, rtol=rtol) + else: + warnings.warn("Did not find the reference data. Output was not checked.") + + +SOFA_CARDIOID_PARAMETERS = [ + ("AKG_c480", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414K", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414N", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414S", "AKG_c480_c414_CUBE", False, False, save_plot), + ("AKG_c414A", "AKG_c480_c414_CUBE", False, False, save_plot), + ("EM_32_0", "EM32_Directivity", False, False, save_plot), + ("EM_32_31", "EM32_Directivity", False, False, save_plot), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + False, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + False, + False, + save_plot, + ), + ("AKG_c480", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414K", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414N", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414S", "AKG_c480_c414_CUBE", True, False, save_plot), + ("AKG_c414A", "AKG_c480_c414_CUBE", True, False, save_plot), + ("EM_32_0", "EM32_Directivity", True, False, save_plot), + ("EM_32_31", "EM32_Directivity", True, False, save_plot), + ( + "Genelec_8020", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + True, + False, + save_plot, + ), + ( + "Vibrolux_2x10inch", + "LSPs_HATS_GuitarCabinets_Akustikmessplatz", + True, + False, + save_plot, + ), +] + + +@pytest.mark.parametrize( + "pattern_id,sofa_file_name,min_phase,save_flag,plot_flag", SOFA_CARDIOID_PARAMETERS +) +def test_sofa_and_cardioid(pattern_id, sofa_file_name, min_phase, save_flag, plot_flag): + """ + Tests with only microphone *or* source from a SOFA file + """ + + if min_phase: + pra.constants.set("octave_bands_n_fft", 128) + pra.constants.set("octave_bands_keep_dc", True) + else: + pra.constants.set("octave_bands_n_fft", 512) + pra.constants.set("octave_bands_keep_dc", False) + + # make sure all the files are present + # in case tests are run out of order + download_sofa_files(overwrite=False, verbose=False) + + # create room + room = pra.ShoeBox( + room_dim, + fs=16000, + max_order=2, + materials=all_materials, + air_absorption=True, + ray_tracing=False, + min_phase=False, + ) + + # define source with figure_eight directivity + dir_factory = MeasuredDirectivityFile( + path=sofa_file_name, + fs=16000, + interp_order=interp_order, + ) + + # add source with figure_eight directivity + if sofa_info[sofa_file_name]["type"] == "microphones": + directivity_SRC = FigureEight( + orientation=DirectionVector(azimuth=90, colatitude=90, degrees=True), + ) + directivity_MIC = dir_factory.get_mic_directivity( + pattern_id, # DirectionVector(azimuth=0, colatitude=0, degrees=True) + Rotation3D([0, 0], rot_order="yz"), + ) + directivity = directivity_MIC + elif sofa_info[sofa_file_name]["type"] == "sources": + directivity_SRC = dir_factory.get_source_directivity( + pattern_id, # DirectionVector(azimuth=0, colatitude=0, degrees=True) + Rotation3D([0, 0], rot_order="yz", degrees=True), + ) + directivity_MIC = FigureEight( + orientation=DirectionVector(azimuth=90, colatitude=90, degrees=True), + ) + directivity = directivity_SRC + else: + raise ValueError("unknown pattern type") + + room.add_source([1.52, 0.883, 1.044], directivity=directivity_SRC) + + # add microphone in its null + room.add_microphone([2.31, 1.65, 1.163], directivity=directivity_MIC) + + # Check set different orientation after intailization of the DIRPATRir class + directivity.set_orientation(Rotation3D([0.0, 0.0], "yz")) + + room.compute_rir() + + rir_1_0 = room.rir[0][0] + + filename = ( + "-".join([sofa_file_name.split(".")[0], pattern_id]) + + f"-minphase_{min_phase}-cardioid.npy" + ) + test_file_path = TEST_DATA / filename + if save_flag: + TEST_DATA.mkdir(exist_ok=True, parents=True) + np.save(test_file_path, rir_1_0) + elif test_file_path.exists(): + reference_data = np.load(test_file_path) + reference_data = reference_data[ref_delay : ref_delay + rir_1_0.shape[0]] + + if plot_flag: + fig, ax = plt.subplots(1, 1) + ax.plot(rir_1_0, label="test") + ax.plot(reference_data, label="ref") + ax.legend() + fig.savefig(test_file_path.with_suffix(".pdf")) + plt.close(fig) + + print("Max diff.:", abs(reference_data - rir_1_0).max()) + print( + "Rel diff.:", + abs(reference_data - rir_1_0).max() / abs(reference_data).max(), + ) + + if not plot_flag: + assert np.allclose(reference_data, rir_1_0, atol=atol, rtol=rtol) + else: + warnings.warn("Did not find the reference data. Output was not checked.") + + +PINV_PARAMETERS = [ + (30, 16, np.pi / 32, np.pi - np.pi / 32, "F", 5e-5, 1e-8), + (30, 16, np.pi / 32, np.pi - np.pi / 32, "C", 5e-5, 1e-8), + (30, 16, 1e-5, np.pi - np.pi / 32, "C", 5e-5, 1e-8), + (30, 16, 1e-5, np.pi - np.pi / 32, "F", 5e-5, 1e-8), + (30, 16, np.pi / 32, np.pi - 1e-5, "C", 5e-5, 1e-8), + (30, 16, np.pi / 32, np.pi - 1e-5, "F", 5e-5, 1e-8), + (35, 17, np.pi / 32, np.pi - np.pi / 32, "F", 5e-5, 1e-8), + (24, 24, np.pi / 32, np.pi - np.pi / 32, "C", 5e-5, 1e-8), + (40, 20, np.pi / 32, np.pi - 1e-5, "C", 5e-5, 1e-8), + (40, 20, np.pi / 32, np.pi - 1e-5, "F", 5e-5, 1e-8), +] + + +@pytest.mark.parametrize( + "n_azimuth, n_col, col_start, col_end, order, atol, rtol", + PINV_PARAMETERS, +) +def test_weighted_pinv(n_azimuth, n_col, col_start, col_end, order, atol, rtol): + azimuth = np.linspace(0, 2 * np.pi, n_azimuth, endpoint=False) + colatitude = np.linspace(col_start, col_end, n_col) + A, C = np.meshgrid(azimuth, colatitude) + alin = A.flatten(order=order) + clin = C.flatten(order=order) + + Ysh = np.random.randn(alin.shape[0]) + + points = np.array( + [ + np.cos(alin) * np.sin(clin), + np.sin(alin) * np.sin(clin), + np.cos(clin), + ] + ).T + + Y_reg, w_reg = calculation_pinv_voronoi_cells(Ysh, clin, colatitude, len(azimuth)) + Y_gen, w_gen = calculation_pinv_voronoi_cells_general(Ysh, points) + + assert np.allclose(w_reg, w_gen, rtol=rtol, atol=atol) + assert np.allclose(Y_reg, Y_gen, rtol=rtol, atol=atol) + + +if __name__ == "__main__": + # generate the test files for regression testing + parser = argparse.ArgumentParser() + parser.add_argument( + "--save", action="store_true", help="save the signal as a reference" + ) + parser.add_argument( + "--plot", action="store_true", help="plot the generated signals" + ) + args = parser.parse_args() + + download_sofa_files(verbose=True) + + for params in SOFA_ONE_SIDE_PARAMETERS: + new_params = params[:-2] + (args.save, args.plot) + test_sofa_one_side(*new_params) + + for params in SOFA_TWO_SIDES_PARAMETERS: + new_params = params[:-2] + (args.save, args.plot) + test_sofa_two_sides(*new_params) + + for params in SOFA_CARDIOID_PARAMETERS: + new_params = params[:-2] + (args.save, args.plot) + test_sofa_and_cardioid(*new_params) + + for params in PINV_PARAMETERS: + test_weighted_pinv(*params) diff --git a/pyroomacoustics/tests/test_source_directivities.py b/pyroomacoustics/directivities/tests/test_source_directivities.py similarity index 80% rename from pyroomacoustics/tests/test_source_directivities.py rename to pyroomacoustics/directivities/tests/test_source_directivities.py index f7176e24..73dd10ba 100644 --- a/pyroomacoustics/tests/test_source_directivities.py +++ b/pyroomacoustics/directivities/tests/test_source_directivities.py @@ -3,11 +3,7 @@ import numpy as np import pyroomacoustics as pra -from pyroomacoustics.directivities import ( - CardioidFamily, - DirectionVector, - DirectivityPattern, -) +from pyroomacoustics.directivities import DirectionVector, FigureEight # create room room = pra.ShoeBox( @@ -18,9 +14,8 @@ ) # define source with figure_eight directivity -PATTERN = DirectivityPattern.FIGURE_EIGHT ORIENTATION = DirectionVector(azimuth=90, colatitude=90, degrees=True) -directivity = CardioidFamily(orientation=ORIENTATION, pattern_enum=PATTERN) +directivity = FigureEight(orientation=ORIENTATION) # add source with figure_eight directivity room.add_source([2.5, 2.5, 2.5], directivity=directivity) diff --git a/pyroomacoustics/doa/grid.py b/pyroomacoustics/doa/grid.py index 09d1e22b..c0c4ed85 100644 --- a/pyroomacoustics/doa/grid.py +++ b/pyroomacoustics/doa/grid.py @@ -10,7 +10,7 @@ import scipy.spatial as sp # import ConvexHull, SphericalVoronoi from .detect_peaks import detect_peaks -from .utils import great_circ_dist +from .utils import cart2spher, fibonacci_spherical_sampling, great_circ_dist, spher2cart class Grid: @@ -38,6 +38,9 @@ def __init__(self, n_points): self.dim = 0 self.values = None + def __len__(self): + return self.cartesian.shape[1] + @abstractmethod def apply(self, func, spherical=False): return NotImplemented @@ -143,8 +146,10 @@ def plot(self, mark_peaks=0): class GridSphere(Grid): """ - This function computes nearly equidistant points on the sphere - using the fibonacci method + This object represents a grid of points of the sphere. + + If the points are not provided, pseudo-uniform points computed + according to the Fibonnaci method are used. Parameters ---------- @@ -152,7 +157,16 @@ class GridSphere(Grid): The number of points to sample spherical_points: ndarray, optional A 2 x n_points array of spherical coordinates with azimuth in - the top row and colatitude in the second row. Overrides n_points. + the top row and colatitude in the second row. Overrides ``n_points`` + and ``cartesian_points``. + cartesian_points: ndarray, optional + A 3 x n_points array of Cartesian coordinates with x, y, z coordinates + in the rows. The vectors are normalized to unit-norm in the constructor. + Overrides ``n_points``. + precompute_neighbors: bool, optional + If `True`, the convex hull algorithm is used to find all + the neighbors of the grid points. This is used for the peak finding + algorithm. References ---------- @@ -160,12 +174,21 @@ class GridSphere(Grid): http://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere """ - def __init__(self, n_points=1000, spherical_points=None): + def __init__( + self, + n_points=1000, + spherical_points=None, + cartesian_points=None, + precompute_neighbors=False, + ): if spherical_points is not None: if spherical_points.ndim != 2 or spherical_points.shape[0] != 2: raise ValueError("spherical_points must be a 2D array with two rows.") - n_points = spherical_points.shape[1] + elif cartesian_points is not None: + if cartesian_points.ndim != 2 or cartesian_points.shape[0] != 3: + raise ValueError("cartesian_points must be a 3D array with two rows.") + n_points = cartesian_points.shape[1] # Parent constructor Grid.__init__(self, n_points) @@ -173,36 +196,37 @@ def __init__(self, n_points=1000, spherical_points=None): self.dim = 3 if spherical_points is not None: - # If a list of points was provided, use it - self.spherical[:, :] = spherical_points + self.cartesian[:] = spher2cart(self.azimuth, self.colatitude) - # transform to cartesian coordinates - self.x[:] = np.cos(self.azimuth) * np.sin(self.colatitude) - self.y[:] = np.sin(self.azimuth) * np.sin(self.colatitude) - self.z[:] = np.cos(self.colatitude) + elif cartesian_points is not None: + # normalize all + norms = np.linalg.norm(cartesian_points, axis=0, keepdims=True) + self.cartesian[:] = cartesian_points / norms + self.azimuth[:], self.colatitude[:], _ = cart2spher(self.cartesian) else: # If no list was provided, samples points on the sphere # as uniformly as possible - # Fibonnaci sampling - offset = 2.0 / n_points - increment = np.pi * (3.0 - np.sqrt(5.0)) - - self.z[:] = (np.arange(n_points) * offset - 1) + offset / 2 - rho = np.sqrt(1.0 - self.z**2) - - phi = np.arange(n_points) * increment - - self.x[:] = np.cos(phi) * rho - self.y[:] = np.sin(phi) * rho + self.x, self.y, self.z = fibonacci_spherical_sampling(n_points) # Create convenient arrays # to access both in cartesian and spherical coordinates self.azimuth[:] = np.arctan2(self.y, self.x) self.colatitude[:] = np.arctan2(np.sqrt(self.x**2 + self.y**2), self.z) + self._neighbors = None + if precompute_neighbors: + self._compute_neighbors() + + @property + def neighbors(self): + if self._neighbors is None: + self._compute_neighbors() + return self._neighbors + + def _compute_neighbors(self): # To perform the peak detection in 2D on a non-squared grid it is # necessary to know the neighboring points of each grid point. The # Convex Hull of points on the sphere is equivalent to the Delauney @@ -226,7 +250,7 @@ def __init__(self, n_points=1000, spherical_points=None): adjacency[tri[2]].add(tri[1]) # convert to list of lists - self.neighbors = [list(x) for x in adjacency] + self._neighbors = [list(x) for x in adjacency] def apply(self, func, spherical=False): """ diff --git a/pyroomacoustics/doa/tests/test_grid.py b/pyroomacoustics/doa/tests/test_grid.py new file mode 100644 index 00000000..bf3c26d6 --- /dev/null +++ b/pyroomacoustics/doa/tests/test_grid.py @@ -0,0 +1,45 @@ +import numpy as np +import pytest +from scipy.spatial import SphericalVoronoi + +from pyroomacoustics.doa import fibonacci_spherical_sampling + + +@pytest.mark.parametrize("n", [20, 100, 200, 500, 1000, 2000, 5000, 10000]) +@pytest.mark.parametrize("tol", [0.06]) +def test_voronoi_area(n, tol): + """ + check the implementation of Fibonacci spherical sampling + + The idea is that the area of the voronoi regions generated + by the sampling should be very close to each other, i.e., + for ``n_points`` they should be approx. ``4.0 * np.pi / n_points``. + + We observed empirically that the relative max error is + around 6% so we set that as the threshold for the test + """ + points = fibonacci_spherical_sampling(n_points=n) + sphere_area = 4.0 * np.pi + area_one_pt = sphere_area / n + + sv = SphericalVoronoi(points=points.T) + areas_voronoi = sv.calculate_areas() + + err = abs(areas_voronoi - area_one_pt) / area_one_pt + max_err = err.max() + avg_err = err.mean() + min_err = err.min() + + print(f"{n=} {max_err=} {avg_err=} {min_err=}") + assert max_err < tol + + +if __name__ == "__main__": + test_voronoi_area(20, 0.01) + test_voronoi_area(100, 0.01) + test_voronoi_area(200, 0.01) + test_voronoi_area(500, 0.01) + test_voronoi_area(1000, 0.001) + test_voronoi_area(2000, 0.001) + test_voronoi_area(5000, 0.001) + test_voronoi_area(10000, 0.001) diff --git a/pyroomacoustics/doa/utils.py b/pyroomacoustics/doa/utils.py index 3bb4af9a..0ff9d149 100644 --- a/pyroomacoustics/doa/utils.py +++ b/pyroomacoustics/doa/utils.py @@ -5,8 +5,12 @@ from __future__ import division +import collections + import numpy as np +RegularGrid = collections.namedtuple("RegularGrid", ["azimuth", "colatitude"]) + def circ_dist(azimuth1, azimuth2, r=1.0): """ @@ -95,7 +99,7 @@ def spher2cart(azimuth, colatitude=None, r=1, degrees=False): r: radius degrees: - Returns values in degrees instead of radians if set to ``True`` + If True, indicates that the input angles are in degree (instead of radian) Returns ------- @@ -167,3 +171,91 @@ def polar_distance(x1, x2): else: index = np.array([index, 1]) return d, index + + +def fibonacci_spherical_sampling(n_points): + """ + This function computes nearly equidistant points on the sphere + using the fibonacci method + + Parameters + ---------- + n_points: int + The number of points to sample + + Returns + ------- + points: numpy.ndarray, (3, n_points) + The cartesian coordinates of the points + + References + ---------- + http://lgdv.cs.fau.de/uploads/publications/spherical_fibonacci_mapping.pdf + http://stackoverflow.com/questions/9600801/evenly-distributing-n-points-on-a-sphere + """ + + points = np.zeros((3, n_points)) + + # Fibonnaci sampling + offset = 2.0 / n_points + increment = np.pi * (3.0 - np.sqrt(5.0)) + + points[2, :] = (np.arange(n_points) * offset - 1) + offset / 2 + rho = np.sqrt(1.0 - points[2] ** 2) + + phi = np.arange(n_points) * increment + + points[0, :] = np.cos(phi) * rho + points[1, :] = np.sin(phi) * rho + + return points + + +def detect_regular_grid(azimuth, colatitude): + """ + This function checks that the linearized azimuth/colatitude where sampled + from a regular grid. + + It also checks that the azimuth are uniformly spread in [0, 2 * np.pi). + The colatitudes can have arbitrary positions. + + Parameters + ---------- + azimuth: numpy.ndarray (npoints,) + The azimuth values in radian + colatitude: numpy.ndarray (npoints,) + The colatitude values in radian + + Returns + ------- + regular_grid: dict["azimuth", "colatitude"] or None + A dictionary with entries for the sorted distinct azimuth an colatitude values + of the grid, if the points form a grid. + Returns `None` if the points do not form a grid. + """ + if len(azimuth) != len(colatitude): + return None + + azimuth_unique = np.unique(azimuth) + colatitude_unique = np.unique(colatitude) + regular_grid = None + if len(azimuth_unique) * len(colatitude_unique) == len(azimuth): + # check that the azimuth are uniformly spread + az_loop = np.insert( + azimuth_unique, len(azimuth_unique), azimuth_unique[0] + 2 * np.pi + ) + delta_az = np.diff(az_loop) + if np.allclose(delta_az, 2 * np.pi / len(azimuth_unique)): + # remake the grid from the unique points and check + # that it matches the original + A, C = np.meshgrid(azimuth_unique, colatitude_unique) + regrid = np.column_stack([A.flatten(), C.flatten()]) + regrid = regrid[np.lexsort(regrid.T), :] + ogrid = np.column_stack([azimuth, colatitude]) + ogrid = ogrid[np.lexsort(ogrid.T), :] + if np.allclose(regrid, ogrid): + regular_grid = RegularGrid( + azimuth=azimuth_unique, colatitude=colatitude_unique + ) + + return regular_grid diff --git a/pyroomacoustics/libroom_src/rir_builder.cpp b/pyroomacoustics/libroom_src/rir_builder.cpp index a2300263..215431bb 100644 --- a/pyroomacoustics/libroom_src/rir_builder.cpp +++ b/pyroomacoustics/libroom_src/rir_builder.cpp @@ -46,15 +46,12 @@ void threaded_rir_builder_impl( py::array_t rir, const py::array_t &time, const py::array_t &alpha, - const py::array_t - &visibility, int fs, size_t fdl, size_t lut_gran, size_t num_threads) { auto pi = get_pi(); // accessors for the arrays auto rir_acc = rir.mutable_unchecked(); - auto vis_acc = visibility.unchecked(); auto tim_acc = time.unchecked(); auto alp_acc = alpha.unchecked(); @@ -65,9 +62,6 @@ void threaded_rir_builder_impl( // error handling if (n_times != size_t(alpha.size())) throw std::runtime_error("time and alpha arrays should have the same size"); - if (n_times != size_t(visibility.size())) - throw std::runtime_error( - "time and visibility arrays should have the same size"); if (fdl % 2 != 1) throw std::runtime_error("the fractional filter length should be odd"); @@ -123,25 +117,23 @@ void threaded_rir_builder_impl( results.emplace_back(pool.enqueue( [&](size_t t_start, size_t t_end, size_t offset) { for (size_t idx = t_start; idx < t_end; idx++) { - if (vis_acc(idx)) { - // decompose integral/fractional parts - T sample_frac = fs * tim_acc(idx); - T time_ip = std::floor(sample_frac); - T time_fp = sample_frac - time_ip; - - // LUT interpolation - T x_off_frac = (1. - time_fp) * lut_gran_f; - T lut_gran_off = std::floor(x_off_frac); - T x_off = (x_off_frac - lut_gran_off); - - int lut_pos = int(lut_gran_off); - int f = int(time_ip) - fdl2; - for (size_t k = 0; k < fdl; lut_pos += lut_gran, f++, k++) - rir_out[offset + f] += - alp_acc(idx) * hann[k] * - (sinc_lut[lut_pos] + - x_off * (sinc_lut[lut_pos + 1] - sinc_lut[lut_pos])); - } + // decompose integral/fractional parts + T sample_frac = fs * tim_acc(idx); + T time_ip = std::floor(sample_frac); + T time_fp = sample_frac - time_ip; + + // LUT interpolation + T x_off_frac = (1. - time_fp) * lut_gran_f; + T lut_gran_off = std::floor(x_off_frac); + T x_off = (x_off_frac - lut_gran_off); + + int lut_pos = int(lut_gran_off); + int f = int(time_ip) - fdl2; + for (size_t k = 0; k < fdl; lut_pos += lut_gran, f++, k++) + rir_out[offset + f] += + alp_acc(idx) * hann[k] * + (sinc_lut[lut_pos] + + x_off * (sinc_lut[lut_pos + 1] - sinc_lut[lut_pos])); } }, t_start, t_end, offset)); @@ -169,16 +161,14 @@ void threaded_rir_builder_impl( } void rir_builder(py::buffer rir, const py::buffer &time, const py::buffer &alpha, - const py::buffer &visibility, int fs, size_t fdl, - size_t lut_gran, size_t num_threads) { - + int fs, size_t fdl, size_t lut_gran, size_t num_threads) { // dispatch to correct implementation depending on input type auto buf = pybind11::array::ensure(rir); if (py::isinstance>(buf)) { - threaded_rir_builder_impl(rir, time, alpha, visibility, fs, fdl, + threaded_rir_builder_impl(rir, time, alpha, fs, fdl, lut_gran, num_threads); } else if (py::isinstance>(buf)) { - threaded_rir_builder_impl(rir, time, alpha, visibility, fs, fdl, + threaded_rir_builder_impl(rir, time, alpha, fs, fdl, lut_gran, num_threads); } else { std::runtime_error("wrong type array for rir builder"); diff --git a/pyroomacoustics/libroom_src/rir_builder.hpp b/pyroomacoustics/libroom_src/rir_builder.hpp index fe7f3557..52f24ceb 100644 --- a/pyroomacoustics/libroom_src/rir_builder.hpp +++ b/pyroomacoustics/libroom_src/rir_builder.hpp @@ -8,8 +8,7 @@ namespace py = pybind11; void rir_builder(py::buffer rir, const py::buffer &time, const py::buffer &alpha, - const py::buffer &visibility, int fs, size_t fdl, - size_t lut_gran, size_t num_threads); + int fs, size_t fdl, size_t lut_gran, size_t num_threads); void delay_sum(const py::buffer irs, const py::buffer delays, py::buffer output, size_t num_threads); diff --git a/pyroomacoustics/parameters.py b/pyroomacoustics/parameters.py index 4c7178fd..978de47f 100644 --- a/pyroomacoustics/parameters.py +++ b/pyroomacoustics/parameters.py @@ -79,6 +79,9 @@ def get_num_threads(): "room_isinside_max_iter": 20, # Max iterations for checking if point is inside room "sinc_lut_granularity": 20, # num. points in integer interval in the sinc interp. LUT "num_threads": get_num_threads(), # num. of threads to use + "octave_bands_n_fft": 512, + "octave_bands_base_freq": 125.0, + "octave_bands_keep_dc": False, } diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py index 3a7c9f63..329842d2 100644 --- a/pyroomacoustics/room.py +++ b/pyroomacoustics/room.py @@ -98,19 +98,20 @@ Until recently, rooms would take an ``absorption`` parameter that was actually **not** the energy absorption we use now. The ``absorption`` parameter is now deprecated and will be removed in the future. - - + + + Randomized Image Method ~~~~~~~~~~~~~~~~~~~~~~~~~ -In highly symmetric shoebox rooms, the regularity of image sources’ positions + +In highly symmetric shoebox rooms, the regularity of image sources’ positions leads to a monotonic convergence in the time arrival of far-field image pairs. -This causes sweeping echoes. The randomized image method adds a small random -displacement to the image source positions, so that they are no longer +This causes sweeping echoes. The randomized image method adds a small random +displacement to the image source positions, so that they are no longer time-aligned, thus reducing sweeping echoes considerably. - To use the randomized method, set the flag ``use_rand_ism`` to True while creating -a room. Additionally, the maximum displacement of the image sources can be +a room. Additionally, the maximum displacement of the image sources can be chosen by setting the parameter ``max_rand_disp``. The default value is 8cm. For a full example see examples/randomized_image_method.py @@ -697,20 +698,12 @@ def callback_mix(premix, snr=0, sir=0, ref_mic=0, n_src=None, n_tgt=None): from . import libroom from .acoustics import OctaveBandsFactory, rt60_eyring, rt60_sabine from .beamforming import MicrophoneArray -from .directivities import CardioidFamily, source_angle_shoebox -from .doa import GridCircle, GridSphere +from .directivities import CardioidFamily, MeasuredDirectivity from .experimental import measure_rt60 from .libroom import Wall, Wall2D -from .parameters import ( - Material, - Physics, - constants, - eps, - get_num_threads, - make_materials, -) +from .parameters import Material, Physics, constants, eps, make_materials +from .simulation import compute_ism_rir, compute_rt_rir from .soundsource import SoundSource -from .utilities import angle_function, fractional_delay def wall_factory(corners, absorption, scattering, name=""): @@ -723,32 +716,6 @@ def wall_factory(corners, absorption, scattering, name=""): raise ValueError("Rooms can only be 2D or 3D") -def sequence_generation(volume, duration, c, fs, max_rate=10000): - # repeated constant - fpcv = 4 * np.pi * c**3 / volume - - # initial time - t0 = ((2 * np.log(2)) / fpcv) ** (1.0 / 3.0) - times = [t0] - - while times[-1] < t0 + duration: - # uniform random variable - z = np.random.rand() - # rate of the point process at this time - mu = np.minimum(fpcv * (t0 + times[-1]) ** 2, max_rate) - # time interval to next point - dt = np.log(1 / z) / mu - - times.append(times[-1] + dt) - - # convert from continuous to discrete time - indices = (np.array(times) * fs).astype(np.int64) - seq = np.zeros(indices[-1] + 1) - seq[indices] = np.random.choice([1, -1], size=len(indices)) - - return seq - - def find_non_convex_walls(walls): """ Finds the walls that are not in the convex hull @@ -869,6 +836,9 @@ class Room(object): max_rand_disp: float, optional; If using randomized image source method, what is the maximum displacement of the image sources? + min_phase: bool, optional + If set to True, generated RIRs will have a minimum phase response. + Cannot be used with ray tracing model. """ def __init__( @@ -886,6 +856,7 @@ def __init__( ray_tracing=False, use_rand_ism=False, max_rand_disp=0.08, + min_phase=False, ): self.walls = walls @@ -907,6 +878,7 @@ def __init__( ray_tracing, use_rand_ism, max_rand_disp, + min_phase, ) # initialize the C++ room engine @@ -936,6 +908,7 @@ def _var_init( ray_tracing, use_rand_ism, max_rand_disp, + min_phase, ): self.fs = fs @@ -948,7 +921,12 @@ def _var_init( self.max_order = max_order self.sigma2_awgn = sigma2_awgn - self.octave_bands = OctaveBandsFactory(fs=self.fs) + self.octave_bands = OctaveBandsFactory( + fs=self.fs, + n_fft=constants.get("octave_bands_n_fft"), + keep_dc=constants.get("octave_bands_keep_dc"), + base_frequency=constants.get("octave_bands_base_freq"), + ) self.max_rand_disp = max_rand_disp # Keep track of the state of the simulator @@ -986,6 +964,8 @@ def _var_init( # initialize the attribute for the impulse responses self.rir = None + self.min_phase = min_phase + def _init_room_engine(self, *args): args = list(args) @@ -1090,7 +1070,7 @@ def _set_ray_tracing_options( if use_ray_tracing: if hasattr(self, "mic_array") and self.mic_array is not None: - if self.mic_array.directivity is not None: + if self.mic_array.is_directive: raise NotImplementedError( "Directivity not supported with ray tracing." ) @@ -1328,7 +1308,12 @@ def from_corners( materials = [Material(0.0, 0.0)] * n_walls # Resample material properties at octave bands - octave_bands = OctaveBandsFactory(fs=fs) + octave_bands = OctaveBandsFactory( + fs=fs, + n_fft=constants.get("octave_bands_n_fft"), + keep_dc=constants.get("octave_bands_keep_dc"), + base_frequency=constants.get("octave_bands_base_freq"), + ) if not Material.all_flat(materials): for mat in materials: mat.resample(octave_bands) @@ -1570,7 +1555,7 @@ def plot( c="k", ) - if plot_directivity and self.mic_array.directivity is not None: + if plot_directivity and self.mic_array.directivity[i] is not None: azimuth_plot = np.linspace( start=0, stop=360, num=361, endpoint=True ) @@ -1783,7 +1768,7 @@ def plot( c="k", ) - if plot_directivity and self.mic_array.directivity is not None: + if plot_directivity and self.mic_array.directivity[i] is not None: azimuth_plot = np.linspace( start=0, stop=360, num=361, endpoint=True ) @@ -2032,6 +2017,9 @@ def add_microphone(self, loc, fs=None, directivity=None): if self.simulator_state["rt_needed"] and directivity is not None: raise NotImplementedError("Directivity not supported with ray tracing.") + if self.dim != 3 and directivity is not None: + raise NotImplementedError("Directivity is only supported for 3D rooms.") + # make sure this is a loc = np.array(loc) @@ -2068,14 +2056,17 @@ def add_microphone_array(self, mic_array, directivity=None): if self.simulator_state["rt_needed"] and directivity is not None: raise NotImplementedError("Directivity not supported with ray tracing.") + if self.dim != 3 and directivity is not None: + raise NotImplementedError("Directivity is only supported for 3D rooms.") + if not isinstance(mic_array, MicrophoneArray): # if the type is not a microphone array, try to parse a numpy array mic_array = MicrophoneArray(mic_array, self.fs, directivity) else: # if the type is microphone array - if directivity is not None: - mic_array.set_directivity(directivity) - if self.simulator_state["rt_needed"] and mic_array.directivity is not None: + mic_array.set_directivity(directivity) + + if self.simulator_state["rt_needed"] and mic_array.is_directive: raise NotImplementedError("Directivity not supported with ray tracing.") return self.add(mic_array) @@ -2104,6 +2095,9 @@ def add_source(self, position, signal=None, delay=0, directivity=None): if self.simulator_state["rt_needed"] and directivity is not None: raise NotImplementedError("Directivity not supported with ray tracing.") + 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 @@ -2114,18 +2108,26 @@ def add_source(self, position, signal=None, delay=0, directivity=None): if isinstance(position, SoundSource): if directivity is not None: - assert isinstance(directivity, CardioidFamily) - return self.add(SoundSource(position, directivity=directivity)) + if isinstance(directivity, CardioidFamily) or isinstance( + directivity, MeasuredDirectivity + ): + return self.add(SoundSource(position, directivity=directivity)) else: return self.add(position) else: if directivity is not None: - assert isinstance(directivity, CardioidFamily) - return self.add( - SoundSource( - position, signal=signal, delay=delay, directivity=directivity + if isinstance(directivity, CardioidFamily) or isinstance( + directivity, MeasuredDirectivity + ): + return self.add( + SoundSource( + position, + signal=signal, + delay=delay, + directivity=directivity, + ) ) - ) + else: return self.add(SoundSource(position, signal=signal, delay=delay)) @@ -2153,15 +2155,27 @@ def image_source_model(self): if n_visible_sources > 0: # Copy to python managed memory - source.images = self.room_engine.sources.copy() - source.orders = self.room_engine.orders.copy() - source.orders_xyz = self.room_engine.orders_xyz.copy() - source.walls = self.room_engine.gen_walls.copy() - source.damping = self.room_engine.attenuations.copy() + + source.images = ( + self.room_engine.sources.copy() + ) # Positions of the image source (3,n) n: n_sources + source.orders = ( + self.room_engine.orders.copy() + ) # Reflection order for each image source shape n:n_sources + source.orders_xyz = ( + self.room_engine.orders_xyz.copy() + ) # Reflection order for each image source for each coordinate shape (3,n) n:n_sources + source.walls = ( + self.room_engine.gen_walls.copy() + ) # Something that i don't get [-1,-1,-1,-1,-1...] shape n:n_sources + source.damping = ( + self.room_engine.attenuations.copy() + ) # Octave band damping's shape (no_of_octave_bands*n_sources) damping value for each image source for each octave bands source.generators = -np.ones(source.walls.shape) # if randomized image method is selected, add a small random # displacement to the image sources + if self.simulator_state["random_ism_needed"]: n_images = np.shape(source.images)[1] @@ -2184,7 +2198,7 @@ def image_source_model(self): else: # if we are here, this means even the direct path is not visible # we set the visibility of the direct path as 0 - self.visibility.append(np.zeros((self.mic_array.M, 1))) + self.visibility.append(np.zeros((self.mic_array.M, 1), dtype=np.int32)) # Update the state self.simulator_state["ism_done"] = True @@ -2208,6 +2222,7 @@ def ray_tracing(self): # reset all the receivers' histograms self.room_engine.reset_mics() + # Basically, histograms for 2 mics corresponding to each source , the histograms are in each octave bands hence (7,2500) 2500 histogram length # update the state self.simulator_state["rt_done"] = True @@ -2226,6 +2241,8 @@ def compute_rir(self): volume_room = self.get_volume() + # Loop over ever microphone present in the room and then for each + # microphone and source pair present in the room for m, mic in enumerate(self.mic_array.R.T): self.rir.append([]) for s, src in enumerate(self.sources): @@ -2236,162 +2253,283 @@ def compute_rir(self): """ # fractional delay length fdl = constants.get("frac_delay_length") - fdl2 = fdl // 2 - # default, just in case both ism and rt are disabled (should never happen) - N = fdl + rir_parts = [] if self.simulator_state["ism_needed"]: - # compute azimuth and colatitude angles for receiver - if self.mic_array.directivity is not None: - angle_function_array = angle_function(src.images, mic) - azimuth = angle_function_array[0] - colatitude = angle_function_array[1] - - # compute azimuth and colatitude angles for source - if self.sources[s].directivity is not None: - azimuth_s, colatitude_s = source_angle_shoebox( - image_source_loc=src.images, - wall_flips=abs(src.orders_xyz), - mic_loc=mic, - ) + ir_ism = compute_ism_rir( + src, + mic, + self.mic_array.directivity[m], + self.visibility[s][m, :], + fdl, + self.c, + self.fs, + self.octave_bands, + air_abs_coeffs=self.air_absorption, + min_phase=self.min_phase, + ) + rir_parts.append(ir_ism) - # compute the distance from image sources - dist = np.sqrt(np.sum((src.images - mic[:, None]) ** 2, axis=0)) - # the RIR building routine works in float32, so we cast here - time = (dist / self.c).astype(np.float32) - t_max = time.max() - N = int(math.ceil(t_max * self.fs)) + if self.simulator_state["rt_needed"]: + ir_rt = compute_rt_rir( + self.rt_histograms[m][s], + self.rt_args["hist_bin_size"], + self.rt_args["hist_bin_size_samples"], + volume_room, + fdl, + self.c, + self.fs, + self.octave_bands, + air_abs_coeffs=self.air_absorption, + ) + rir_parts.append(ir_rt) + if len(rir_parts) == 0: + raise ValueError("Both ISM and RT are disabled") + elif len(rir_parts) == 1: + rir = rir_parts[0] else: - t_max = 0.0 + max_len = max([r.shape[0] for r in rir_parts]) + rir = np.zeros(max_len) + for r in rir_parts: + rir[: r.shape[0]] += r - if self.simulator_state["rt_needed"]: - # get the maximum length from the histograms - nz_bins_loc = np.nonzero(self.rt_histograms[m][s][0].sum(axis=0))[0] - if len(nz_bins_loc) == 0: - n_bins = 0 - else: - n_bins = nz_bins_loc[-1] + 1 - - t_max = np.maximum(t_max, n_bins * self.rt_args["hist_bin_size"]) - - # the number of samples needed - # round up to multiple of the histogram bin size - # add the lengths of the fractional delay filter - hbss = int(self.rt_args["hist_bin_size_samples"]) - N = int(math.ceil(t_max * self.fs / hbss) * hbss) - - # this is where we will compose the RIR - - # here we create an array of the right length to - # receiver the full RIR - # the +1 is due to some rare cases where numerical - # errors push the last sample over the end of the - # array - ir = np.zeros(N + fdl + 1) - - # This is the distance travelled wrt time - distance_rir = np.arange(N) / self.fs * self.c - - # this is the random sequence for the tail generation - seq = sequence_generation(volume_room, N / self.fs, self.c, self.fs) - seq = seq[:N] - - # Do band-wise RIR construction - is_multi_band = self.is_multi_band - bws = self.octave_bands.get_bw() if is_multi_band else [self.fs / 2] - rir_bands = [] - - for b, bw in enumerate(bws): - ir_loc = np.zeros_like(ir, dtype=np.float32) - - # IS method - if self.simulator_state["ism_needed"]: - alpha = src.damping[b, :] / dist - - if self.mic_array.directivity is not None: - alpha *= self.mic_array.directivity[m].get_response( - azimuth=azimuth, - colatitude=colatitude, - frequency=bw, - degrees=False, - ) - - if self.sources[s].directivity is not None: - alpha *= self.sources[s].directivity.get_response( - azimuth=azimuth_s, - colatitude=colatitude_s, - frequency=bw, - degrees=False, - ) - - vis = self.visibility[s][m, :].astype(np.int32) - # we add the delay due to the factional delay filter to - # the arrival times to avoid problems when propagation - # is shorter than the delay to to the filter - # hence: time + fdl2 - time_adjust = time + fdl2 / self.fs - libroom.rir_builder( - ir_loc, - time_adjust.astype(np.float32), - alpha.astype(np.float32), - vis, - self.fs, - fdl, - constants.get("sinc_lut_granularity"), - constants.get("num_threads"), - ) + self.rir[m].append(rir) - if is_multi_band: - ir_loc = self.octave_bands.analysis(ir_loc, band=b) + self.simulator_state["rir_done"] = True - ir += ir_loc + def dft_scale_rir_calc( + self, + attenuations, + dist, + time, + bws, + N, + azi_m, + col_m, + azi_s, + col_s, + src_pos=0, + mic_pos=0, + ): + """ + Full DFT scale RIR construction. - # Ray Tracing - if self.simulator_state["rt_needed"]: - if is_multi_band: - seq_bp = self.octave_bands.analysis(seq, band=b) - else: - seq_bp = seq.copy() + This function also takes into account the FIR's of the source and receiver retrieved from the SOFA file. - # interpolate the histogram and multiply the sequence - seq_bp_rot = seq_bp.reshape((-1, hbss)) - new_n_bins = seq_bp_rot.shape[0] - hist = self.rt_histograms[m][s][0][b, :new_n_bins] - normalization = np.linalg.norm(seq_bp_rot, axis=1) - indices = normalization > 0.0 - seq_bp_rot[indices, :] /= normalization[indices, None] - seq_bp_rot *= np.sqrt(hist[:, None]) + Parameters + ---------- + attenuations: arr + Dampings for all the image sources Shape : ( No_of_octave_band x no_img_src) + dist : arr + distance of all the image source present in the room from this particular mic Shape : (no_img_src) + time : arr + Time of arrival of all the image source Shape : (no_img_src) + bws : + bandwidth of all the octave bands + N : + azi_m : arr + Azimuth angle of arrival of this particular mic for all image sources Shape : (no_img_src) + col_m : arr + Colatitude angle of arrival of this particular mic for all image sources Shape : (no_img_src) + azi_s : arr + Azimuth angle of departure of this particular source for all image sources Shape : (no_img_src) + col_s : arr + Colatitude angle of departure of this particular source for all image sources Shape : (no_img_src) + src_pos : int + The particular source we are calculating RIR + mic_pos : int + The particular mic we are calculating RIR - # Normalize the band power - # The bands should normally sum up to fs / 2 - seq_bp *= np.sqrt(bw / self.fs * 2.0) + Returns + ------- + rir : :py:class:`~numpy.ndarray` + Constructed RIR for this particlar src mic pair . - ir_loc[fdl2 : fdl2 + N] += seq_bp + The constructed RIR still lacks air absorption and distance absorption because in the old pyroom these calculation happens on the octave band level. - # keep for further processing - rir_bands.append(ir_loc) - # Do Air absorption - if self.simulator_state["air_abs_needed"]: - # In case this was not multi-band, do the band pass filtering - if len(rir_bands) == 1: - rir_bands = self.octave_bands.analysis(rir_bands[0]).T + """ - # Now apply air absorption - for band, air_abs in zip(rir_bands, self.air_absorption): - air_decay = np.exp(-0.5 * air_abs * distance_rir) - band[fdl2 : N + fdl2] *= air_decay + attenuations = attenuations / dist + alp = [] + window_length = 81 - # Sum up all the bands - np.sum(rir_bands, axis=0, out=ir) + no_imag_src = attenuations.shape[1] - self.rir[-1].append(ir) + fp_im = N + fir_length_octave_band = self.octave_bands.n_fft - self.simulator_state["rir_done"] = True + from .build_rir import ( + fast_convolution_3, + fast_convolution_4, + fast_window_sinc_interpolator, + ) + + rec_presence = True if (len(azi_m) > 0 and len(col_m) > 0) else False + source_presence = True if (len(azi_s) > 0 and len(col_s) > 0) else False + + final_fir_IS_len = ( + (self.mic_array.directivity[mic_pos].filter_len_ir if (rec_presence) else 1) + + ( + self.sources[src_pos].directivity.filter_len_ir + if (source_presence) + else 1 + ) + + window_length + + fir_length_octave_band + ) - 3 + + if rec_presence and source_presence: + resp_mic = self.mic_array.directivity[mic_pos].get_response( + azimuth=azi_m, colatitude=col_m, degrees=False + ) # Return response as an array of number of (img_sources * length of filters) + resp_src = self.sources[src_pos].directivity.get_response( + azimuth=azi_s, colatitude=col_s, degrees=False + ) + + if self.mic_array.directivity[mic_pos].filter_len_ir == 1: + resp_mic = np.array(resp_mic).reshape(-1, 1) + + else: + assert ( + self.fs == self.mic_array.directivity[mic_pos].fs + ), "Mic directivity: frequency of simulation should be same as frequency of interpolation" + + if self.sources[src_pos].directivity.filter_len_ir == 1: + resp_src = np.array(resp_src).reshape(-1, 1) + else: + assert ( + self.fs == self.sources[src_pos].directivity.fs + ), "Source directivity: frequency of simulation should be same as frequency of interpolation" + + else: + if source_presence: + assert ( + self.fs == self.sources[src_pos].directivity.fs + ), "Directivity source frequency of simulation should be same as frequency of interpolation" + + resp_src = self.sources[src_pos].directivity.get_response( + azimuth=azi_s, + colatitude=col_s, + degrees=False, + ) + + elif rec_presence: + assert ( + self.fs == self.mic_array.directivity[mic_pos].fs + ), "Directivity mic frequency of simulation should be same as frequency of interpolation" + + resp_mic = self.mic_array.directivity[mic_pos].get_response( + azimuth=azi_m, + colatitude=col_m, + degrees=False, + ) + + # else: + # txt = "No" + # final_fir_IS_len = (fir_length_octave_band + window_length) - 1 + + time_arrival_is = time # For min phase + + # Calculating fraction delay sinc filter + sample_frac = time_arrival_is * self.fs # Find the fractional sample number + + ir_diff = np.zeros(N + (final_fir_IS_len)) # 2050 #600 + + # Create arrays for fractional delay low pass filter, sum of {damping coeffiecients * octave band filter}, source response, receiver response. + + cpy_ir_len_1 = np.zeros((no_imag_src, final_fir_IS_len), dtype=np.complex_) + cpy_ir_len_2 = np.zeros((no_imag_src, final_fir_IS_len), dtype=np.complex_) + cpy_ir_len_3 = np.zeros((no_imag_src, final_fir_IS_len), dtype=np.complex_) + cpy_ir_len_4 = np.zeros((no_imag_src, final_fir_IS_len), dtype=np.complex_) + att_in_dft_scale = np.zeros( + (no_imag_src, fir_length_octave_band), dtype=np.complex_ + ) + + # Vectorized sinc filters + + vectorized_interpolated_sinc = np.zeros( + (no_imag_src, window_length), dtype=np.double + ) + vectorized_time_ip = np.array( + [int(math.floor(sample_frac[img_src])) for img_src in range(no_imag_src)] + ) + vectorized_time_fp = [ + sample_frac[img_src] - int(math.floor(sample_frac[img_src])) + for img_src in range(no_imag_src) + ] + vectorized_time_fp = np.array(vectorized_time_fp, dtype=np.double) + vectorized_interpolated_sinc = fast_window_sinc_interpolator( + vectorized_time_fp, window_length, vectorized_interpolated_sinc + ) + + for i in range(no_imag_src): # Loop through Image source + att_in_octave_band = attenuations[:, i] + att_in_dft_scale_ = att_in_dft_scale[i, :] + + # Interpolating attenuations given in the single octave band to a DFT scale. + + att_in_dft_scale_ = self.octave_bands.octave_band_dft_interpolation( + att_in_octave_band, + self.air_absorption, + dist[i], + att_in_dft_scale_, + bws, + self.min_phase, + ) + + # time_ip = int(math.floor(sample_frac[i])) # Calculating the integer sample + + # time_fp = sample_frac[i] - time_ip # Calculating the fractional sample + + # windowed_sinc_filter = fast_window_sinc_interpolater(time_fp) + + cpy_ir_len_1[i, : att_in_dft_scale_.shape[0]] = np.fft.ifft( + att_in_dft_scale_ + ) + cpy_ir_len_2[i, :window_length] = vectorized_interpolated_sinc[i, :] + + if source_presence and rec_presence: + cpy_ir_len_3[i, : resp_src[i, :].shape[0]] = resp_src[i, :] + + cpy_ir_len_4[i, : resp_mic[i, :].shape[0]] = resp_mic[i, :] + + out = fast_convolution_4( + cpy_ir_len_1[i, :], + cpy_ir_len_2[i, :], + cpy_ir_len_3[i, :], + cpy_ir_len_4[i, :], + final_fir_IS_len, + ) + + ir_diff[ + vectorized_time_ip[i] : (vectorized_time_ip[i] + final_fir_IS_len) + ] += np.real(out) + + else: + if source_presence: + resp = resp_src[i, :] + elif rec_presence: + resp = resp_mic[i, :] + + cpy_ir_len_3[i, : resp.shape[0]] = resp + + out = fast_convolution_3( + cpy_ir_len_1[i, :], + cpy_ir_len_2[i, :], + cpy_ir_len_3[i, :], + final_fir_IS_len, + ) + + ir_diff[ + vectorized_time_ip[i] : (vectorized_time_ip[i] + final_fir_IS_len) + ] += np.real(out) + + return ir_diff def simulate( self, @@ -2873,6 +3011,7 @@ def __init__( ray_tracing=False, use_rand_ism=False, max_rand_disp=0.08, + min_phase=False, ): p = np.array(p, dtype=np.float32) @@ -2896,6 +3035,7 @@ def __init__( ray_tracing, use_rand_ism, max_rand_disp, + min_phase, ) # Keep the correctly ordered naming of walls diff --git a/pyroomacoustics/simulation/__init__.py b/pyroomacoustics/simulation/__init__.py new file mode 100644 index 00000000..ede5538c --- /dev/null +++ b/pyroomacoustics/simulation/__init__.py @@ -0,0 +1,2 @@ +from .ism import compute_ism_rir +from .rt import compute_rt_rir diff --git a/pyroomacoustics/simulation/ism.py b/pyroomacoustics/simulation/ism.py new file mode 100644 index 00000000..02479179 --- /dev/null +++ b/pyroomacoustics/simulation/ism.py @@ -0,0 +1,334 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2024 Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Internal routines used for simulation using the image source method. In +particular, how to transform the image sources, attenuations, etc, obtained +from the core simulation engine into impulse responses. +""" +import math + +import numpy as np +from scipy.signal import fftconvolve, hilbert + +from .. import libroom +from ..parameters import constants +from ..utilities import angle_function + + +def multi_convolve(*signals): + max_len = signals[0].shape[-1] + shape = signals[0].shape[:-1] + for s in signals[1:]: + if shape != s.shape[:-1]: + raise ValueError("All signals to convolve should have same batch shape") + + max_len = max_len + s.shape[1] - 1 + + pow2_len = 2 ** math.ceil(np.log2(max_len)) + + fd = np.fft.rfft(signals[0], axis=-1, n=pow2_len) + for s in signals[1:]: + fd *= np.fft.rfft(s, axis=-1, n=pow2_len) + + conv = np.fft.irfft(fd, axis=-1, n=pow2_len) + conv = conv[:, :max_len] + + return conv + + +def apply_air_aborption( + oct_band_amplitude, + air_abs_coeffs, + distance, +): + air_abs_factor = np.exp(-0.5 * air_abs_coeffs[:, None] * distance) + return oct_band_amplitude * air_abs_factor + + +def interpolate_octave_bands( + octave_bands, + att_in_octave_bands, + min_phase=True, +): + """ + Convert octave band dampings to dft scale, interpolates octave band values to full dft scale values. + + Parameters + ---------- + octave_bands: OctaveBands + The octave bands object that contains the filters + att_in_octave_band : np.ndarray + Dampings in octave band Shape : (no_of_octave_band) + air_abs_band : np.ndarray + air absorption in octave band Shape : (no_of_octave_band) + min_phase : Boolean + decides if the final filter is minimum phase (causal) or (non-causal) linear phase sinc filter + + Returns: + ------------- + att_in_dft_scale : np.ndarray + Dampings in octave bands interpolated to full scale frequency domain. + + """ + n_bands = octave_bands.n_bands + + att_in_dft_scale = np.einsum( + "bi,fb->if", att_in_octave_bands, octave_bands.filters_freq_domain + ) + + if min_phase: + att_in_dft_scale += ( + 1e-07 # To avoid divide by zero error when performing hilbert transform. + ) + m_p = np.imag(-hilbert(np.log(np.abs(att_in_dft_scale)), axis=-1)) + att_in_dft_scale = np.abs(att_in_dft_scale) * np.exp(1j * m_p) + + ir = np.fft.ifft(att_in_dft_scale, n=octave_bands.n_fft, axis=-1).real + + 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, + is_visible, + fdl, + c, + fs, + octave_bands, + min_phase=True, + air_abs_coeffs=None, +): + fdl2 = fdl // 2 + + images = src.images[:, is_visible] + att = src.damping[:, is_visible] + + dist = np.sqrt( + np.sum((images - mic[:, None]) ** 2, axis=0) + ) # Calculate distance between image sources and for each microphone + + # dist shape (n) : n0 of image sources + time = dist / c # Calculate time of arrival for each image source + + # we add the delay due to the factional delay filter to + # the arrival times to avoid problems when propagation + # is shorter than the delay to to the filter + # hence: time + fdl2 + delay = fdl2 / fs + time += delay + + t_max = ( + time.max() + ) # The image source which takes the most time to arrive to this particular microphone + + # Here we create an array of the right length to + # receiver the full RIR + # The +1 is due to some rare cases where numerical + # errors push the last sample over the end of the + # array + N = int(math.ceil(t_max * fs + fdl2 + 1)) + 1 + + oct_band_amplitude = att / dist + full_band_imp_resp = [] + + if air_abs_coeffs is not None: + oct_band_amplitude = apply_air_aborption( + oct_band_amplitude, air_abs_coeffs, dist + ) + + if mic_dir is not None: + angle_function_array = angle_function(images, mic) + azimuth_m = angle_function_array[0] + colatitude_m = angle_function_array[1] + mic_gain = mic_dir.get_response( + azimuth=azimuth_m, + colatitude=colatitude_m, + degrees=False, + ) + + if mic_dir.is_impulse_response: + full_band_imp_resp.append(mic_gain) + else: + 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, + ) + src_gain = src.directivity.get_response( + azimuth=azimuth_s, + colatitude=colatitude_s, + degrees=False, + ) + + if src.directivity.is_impulse_response: + full_band_imp_resp.append(src_gain) + else: + oct_band_amplitude *= src_gain + + # there should be 3 possibile shapes for the gains + # 1) (n_images,) for freq flat sources + # 2) (n_images, n_octave_bands) for directivites defined by octave bands + # 3) (n_images, n_taps) for directivites defined as impulse responses + # Cases 2-3 are ambiguous, although we will usually have no_of_octave_bands == 7 + # and n_taps > 7 + # Proposed solution: add a type for IR type of impulse responses + # (MeasuredDirectivity only for now) + # + # Then, we can have damping coefficients either + # 1) (1, n_images) flat + # 2) (n_octave_bands, n_images) per octave bands + # + # flat/octave bands + # 1) mic_gain * src_gain * damping + # 2) run rir_builder + # + # with impulse response + # 1) compute flat/oct. bands gains + # mic_gain * damping or src_gain * damping or damping + # 2) compute impulse response part + # mic_gain or src_gain or convolve(mic_gain, src_gain) + # 3) run the dft_scale_rir_calc routine + + n_bands = oct_band_amplitude.shape[0] + + if len(full_band_imp_resp) > 0: + # Case 3) Full band RIR construction + sample_frac = time * fs + time_ip = np.floor(sample_frac).astype(np.int32) + time_fp = (sample_frac - time_ip).astype(np.float32) + + # create fractional delay filters + frac_delays = np.zeros((time_fp.shape[0], fdl), dtype=np.float32) + libroom.fractional_delay( + frac_delays, + time_fp, + constants.get("sinc_lut_granularity"), + constants.get("num_threads"), + ) + full_band_imp_resp.append(frac_delays) + + # convolve all the impulse responses + if n_bands == 1: + irs = multi_convolve(*full_band_imp_resp) + irs *= oct_band_amplitude.T + + else: + ir_att = interpolate_octave_bands( + octave_bands, oct_band_amplitude, min_phase=min_phase + ) + full_band_imp_resp.append(ir_att) + irs = multi_convolve(*full_band_imp_resp) + + # now overlap-add all the short impulse responses + n_max = int(time_ip.max() + irs.shape[1]) + rir = np.zeros(n_max, dtype=np.float32) + libroom.delay_sum( + irs.astype(np.float32), time_ip, rir, constants.get("num_threads") + ) + + if n_bands > 1 and not min_phase: + # we want to trim the extra samples introduced by the octave + # band filters + s = (constants.get("octave_bands_n_fft")) // 2 + rir = rir[s:] + + else: + # Case 1) or 2) + # Single- or Octave-bands RIR construction + n_bands = oct_band_amplitude.shape[0] + rir = np.zeros(N, dtype=np.float32) # ir for every band + for b in range(n_bands): # Loop through every band + ir_loc = np.zeros(N, dtype=np.float32) # ir for every band + libroom.rir_builder( + ir_loc, + time.astype(np.float32), + oct_band_amplitude[b].astype(np.float32), + fs, + fdl, + constants.get("sinc_lut_granularity"), + constants.get("num_threads"), + ) + + if n_bands > 1: + rir += octave_bands.analysis(ir_loc, band=b) + else: + rir += ir_loc + + return rir diff --git a/pyroomacoustics/simulation/rt.py b/pyroomacoustics/simulation/rt.py new file mode 100644 index 00000000..8ee6f3c2 --- /dev/null +++ b/pyroomacoustics/simulation/rt.py @@ -0,0 +1,175 @@ +# Some classes to apply rotate objects or indicate directions in 3D space. +# Copyright (C) 2024 Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +Internal routines used for simulation using the ray tracing method. +In particular, how to transform the histograms obtained from the core +simulation engine into impulse responses. +""" +import math + +import numpy as np +from scipy.interpolate import interp1d + + +def sequence_generation(volume, duration, c, fs, max_rate=10000): + # repeated constant + fpcv = 4 * np.pi * c**3 / volume + + # initial time + t0 = ((2 * np.log(2)) / fpcv) ** (1.0 / 3.0) + times = [t0] + + while times[-1] < t0 + duration: + # uniform random variable + z = np.random.rand() + # rate of the point process at this time + mu = np.minimum(fpcv * (t0 + times[-1]) ** 2, max_rate) + # time interval to next point + dt = np.log(1 / z) / mu + + times.append(times[-1] + dt) + + # convert from continuous to discrete time + + indices = (np.array(times) * fs).astype(np.int64) + seq = np.zeros(indices[-1] + 1) + seq[indices] = np.random.choice([1, -1], size=len(indices)) + + return seq + + +def interp_hist(hist, N): + """ + interpolate the histogram on N samples + + we use the bin centers as anchors + since we can't interpolate outside the anchors, we just repeat + the first and last value to fill the array + """ + hbss = N // hist.shape[0] + pad = hbss // 2 + t = np.linspace(pad, N - 1 - pad, hist.shape[0]) + f = interp1d(t, hist) + out = np.zeros(N) + out[pad:-pad] = f(np.arange(pad, N - pad)) + out[:pad] = out[pad] + out[-pad:] = out[-pad] + return out + + +def seq_bin_power(seq, hbss): + seq_rot = seq.reshape((-1, hbss)) # shape 72,64 + power = np.sum(seq_rot**2, axis=1) + power[power <= 0.0] = 1.0 + return power + + +def compute_rt_rir( + histograms, + hist_bin_size, + hist_bin_size_samples, + volume_room, + fdl, + c, + fs, + octave_bands, + air_abs_coeffs=None, +): + # get the maximum length from the histograms + # Sum vertically across octave band for each value in + # histogram (7,2500) -> (2500) -> np .nonzero( + nz_bins_loc = np.nonzero(histograms[0].sum(axis=0))[0] + + if len(nz_bins_loc) == 0: + # the histogram is all zeros, there is no RIR to build + # we return only an RIR that contains the default delay + return np.zeros(fdl // 2) + else: + n_bins = nz_bins_loc[-1] + 1 + + t_max = n_bins * hist_bin_size + + # N changes here , the length of RIR changes if we apply RT method. + # the number of samples needed + # round up to multiple of the histogram bin size + # add the lengths of the fractional delay filter + + hbss = int(hist_bin_size_samples) + + fdl2 = fdl // 2 # delay due to fractional delay filter + # make the length a multiple of the bin size + n_bins = math.ceil(t_max * fs / hbss) + N = n_bins * hbss + + # this is the random sequence for the tail generation + seq = sequence_generation(volume_room, N / fs, c, fs) + seq = seq[:N] # take values according to N as seq is larger + + n_bands = histograms[0].shape[0] + bws = octave_bands.get_bw() if n_bands > 1 else [fs / 2] + + rir_bands = np.zeros((n_bands, N)) + for b, bw in enumerate(bws): # Loop through every band + if n_bands > 1: + seq_bp = octave_bands.analysis(seq, band=b) + else: + seq_bp = seq.copy() + + # Take only those bins which have some non-zero values for that specific octave bands. + hist = histograms[0][b, :n_bins] + + # we normalize the histogram by the sequence power in that bin + seq_power = seq_bin_power(seq_bp, hbss) + + # we linarly interpolate the histogram to smoothly cover all the samples + hist_lin_interp = interp_hist(hist / seq_power, N) + + # Normalize the band power + # the (bw / fs * 2.0) is necessary to normalize the band power + # this is the contribution of the octave band to total energy + seq_bp *= np.sqrt(bw / fs * 2.0 * hist_lin_interp) + + # Impulse response for every octave band for each microphone + rir_bands[b] = seq_bp + + if air_abs_coeffs is not None: + if rir_bands.shape[0] == 1: + # do the octave band analysis if it was not done in the first step + rir_bands = np.array( + [ + octave_bands.analysis(rir_bands[0], band=b) + for b in range(octave_bands.n_bands) + ] + ) + + air_abs = np.exp( + -0.5 * air_abs_coeffs[:, None] * np.arange(N)[None, :] / fs * c + ) + + rir_bands *= air_abs + + rir = np.zeros(fdl2 + N) + rir[fdl2:] = np.sum(rir_bands, axis=0) + + return rir diff --git a/pyroomacoustics/tests/test_build_rir.py b/pyroomacoustics/tests/test_build_rir.py index 81feb025..85503b51 100644 --- a/pyroomacoustics/tests/test_build_rir.py +++ b/pyroomacoustics/tests/test_build_rir.py @@ -77,7 +77,7 @@ ) -def build_rir_wrap(time, alpha, visibility, fs, fdl): +def build_rir_wrap(time, alpha, fs, fdl): # fractional delay length fdl = pra.constants.get("frac_delay_length") fdl2 = (fdl - 1) // 2 @@ -91,12 +91,11 @@ def build_rir_wrap(time, alpha, visibility, fs, fdl): ir_cpp_f = ir_cython.astype(np.float32) # Try to use the Cython extension - # build_rir.fast_rir_builder(ir_cython, time, alpha, visibility, fs, fdl) + # build_rir.fast_rir_builder(ir_cython, time, alpha, fs, fdl) libroom.rir_builder( ir_cpp_f, time.astype(np.float32), alpha.astype(np.float32), - visibility.astype(np.int32), fs, fdl, 20, @@ -105,12 +104,11 @@ def build_rir_wrap(time, alpha, visibility, fs, fdl): # fallback to pure Python implemenation for i in range(time.shape[0]): - if visibility[i] == 1: - time_ip = int(np.round(fs * time[i])) - time_fp = (fs * time[i]) - time_ip - ir_ref[time_ip - fdl2 : time_ip + fdl2 + 1] += alpha[ - i - ] * pra.fractional_delay(time_fp) + time_ip = int(np.round(fs * time[i])) + time_fp = (fs * time[i]) - time_ip + ir_ref[time_ip - fdl2 : time_ip + fdl2 + 1] += alpha[i] * pra.fractional_delay( + time_fp + ) return ir_ref, ir_cpp_f @@ -120,9 +118,7 @@ def test_build_rir(): return for t, a, v in zip(times, alphas, visibilities): - ir_ref, ir_cython = build_rir_wrap( - times[0], alphas[0], visibilities[0], fs, fdl - ) + ir_ref, ir_cython = build_rir_wrap(times[0], alphas[0], fs, fdl) assert np.max(np.abs(ir_ref - ir_cython)) < tol @@ -139,14 +135,12 @@ def test_short(): time = np.array([0.0], dtype=np.float32) alpha = np.array([1.0], dtype=np.float32) - visibility = np.array([1], dtype=np.int32) with pytest.raises(RuntimeError): libroom.rir_builder( rir, time, alpha, - visibility, fs, fdl, 20, @@ -167,10 +161,9 @@ def test_long(): time = np.array([(N - 1) / fs], dtype=np.float32) alpha = np.array([1.0], dtype=np.float32) - visibility = np.array([1], dtype=np.int32) with pytest.raises(RuntimeError): - libroom.rir_builder(rir, time, alpha, visibility, fs, fdl, 20, 2) + libroom.rir_builder(rir, time, alpha, fs, fdl, 20, 2) def test_errors(): @@ -186,16 +179,12 @@ def test_errors(): time = np.array([100 / fs, 200 / fs], dtype=np.float32) alpha = np.array([1.0, 1.0], dtype=np.float32) - visibility = np.array([1, 1], dtype=np.int32) - - with pytest.raises(RuntimeError): - libroom.rir_builder(rir, time, alpha[:1], visibility, fs, fdl, 20, 2) with pytest.raises(RuntimeError): - libroom.rir_builder(rir, time, alpha, visibility[:1], fs, fdl, 20, 2) + libroom.rir_builder(rir, time, alpha[:1], fs, fdl, 20, 2) with pytest.raises(RuntimeError): - libroom.rir_builder(rir, time, alpha, visibility, fs, 80, 20, 2) + libroom.rir_builder(rir, time, alpha, fs, 80, 20, 2) @pytest.mark.parametrize("dtype,tol", [(np.float32, 1e-5), (np.float64, 1e-7)]) @@ -260,22 +249,17 @@ def measure_runtime(dtype=np.float32, num_threads=4): alpha = np.random.randn(n_img).astype(dtype) rir_len = int(np.ceil(time_arr.max() * fs) + fdl) rir = np.zeros(rir_len, dtype=dtype) - visibility = np.random.rand(n_img) > 0.5 tick = time.perf_counter() rir[:] = 0.0 for i in range(n_repeat): - libroom.rir_builder( - rir, time_arr, alpha, visibility.astype(np.int32), fs, fdl, 20, 1 - ) + libroom.rir_builder(rir, time_arr, alpha, fs, fdl, 20, 1) tock_1 = (time.perf_counter() - tick) / n_repeat tick = time.perf_counter() rir[:] = 0.0 for i in range(n_repeat): - libroom.rir_builder( - rir, time_arr, alpha, visibility.astype(np.int32), fs, fdl, 20, num_threads - ) + libroom.rir_builder(rir, time_arr, alpha, fs, fdl, 20, num_threads) tock_8 = (time.perf_counter() - tick) / n_repeat tick = time.perf_counter() @@ -294,9 +278,7 @@ def measure_runtime(dtype=np.float32, num_threads=4): tick = time.perf_counter() rir[:] = 0.0 for i in range(n_repeat): - build_rir.fast_rir_builder( - rir, time_arr, alpha, visibility.astype(np.int32), fs, fdl - ) + build_rir.fast_rir_builder(rir, time_arr, alpha, fs, fdl) tock_old = (time.perf_counter() - tick) / n_repeat print("runtime:") diff --git a/pyroomacoustics/tests/test_source_directivity_flipping.py b/pyroomacoustics/tests/test_source_directivity_flipping.py index 46bcfc76..0495e617 100644 --- a/pyroomacoustics/tests/test_source_directivity_flipping.py +++ b/pyroomacoustics/tests/test_source_directivity_flipping.py @@ -3,7 +3,7 @@ import numpy as np import pyroomacoustics as pra -import pyroomacoustics.directivities +from pyroomacoustics.simulation.ism import source_angle_shoebox class TestSourceDirectivityFlipping(TestCase): @@ -24,7 +24,7 @@ def test_x(self): room.image_source_model() # compute azimuth_s and colatitude_s pair for images along x-axis - source_angle_array = pyroomacoustics.directivities.source_angle_shoebox( + source_angle_array = source_angle_shoebox( image_source_loc=room.sources[0].images, wall_flips=abs(room.sources[0].orders_xyz), mic_loc=mic_loc, @@ -61,7 +61,7 @@ def test_y(self): room.image_source_model() # compute azimuth_s and colatitude_s pair for images along x-axis - source_angle_array = pyroomacoustics.directivities.source_angle_shoebox( + source_angle_array = source_angle_shoebox( image_source_loc=room.sources[0].images, wall_flips=abs(room.sources[0].orders_xyz), mic_loc=mic_loc, @@ -95,7 +95,7 @@ def test_z(self): room.image_source_model() # compute azimuth_s and colatitude_s pair for images along z-axis - source_angle_array = pyroomacoustics.directivities.source_angle_shoebox( + source_angle_array = source_angle_shoebox( image_source_loc=room.sources[0].images, wall_flips=abs(room.sources[0].orders_xyz), mic_loc=mic_loc, diff --git a/pyroomacoustics/utilities.py b/pyroomacoustics/utilities.py index 7155b49b..cd5cc3d9 100644 --- a/pyroomacoustics/utilities.py +++ b/pyroomacoustics/utilities.py @@ -23,9 +23,11 @@ # not, see . from __future__ import division +import functools import itertools import numpy as np +import soxr from scipy import signal from scipy.io import wavfile @@ -35,6 +37,7 @@ def requires_matplotlib(func): + @functools.wraps(func) # preserves name, docstrings, and signature of function def function_wrapper(*args, **kwargs): try: import matplotlib.pyplot as plt @@ -231,9 +234,9 @@ def highpass(signal, Fs, fc=None, plot=False): wc = 2.0 * fc / Fs # design the filter - from scipy.signal import freqz, iirfilter, lfilter + from scipy.signal import iirfilter, sosfiltfilt, sosfreqz - b, a = iirfilter(n, Wn=wc, rp=rp, rs=rs, btype="highpass", ftype=type) + sos = iirfilter(n, Wn=wc, rp=rp, rs=rs, btype="highpass", ftype=type, output="sos") # plot frequency response of filter if requested if plot: @@ -245,7 +248,7 @@ def highpass(signal, Fs, fc=None, plot=False): warnings.warn("Matplotlib is required for plotting") return - w, h = freqz(b, a) + w, h = sosfreqz(sos) plt.figure() plt.title("Digital filter frequency response") @@ -256,7 +259,7 @@ def highpass(signal, Fs, fc=None, plot=False): plt.grid() # apply the filter - signal = lfilter(b, a, signal.copy()) + signal = sosfiltfilt(sos, signal.copy()) return signal @@ -816,3 +819,43 @@ def angle_function(s1, v2): co[:] = np.pi / 2 return np.vstack((az, co)) + + +def resample(data, old_fs, new_fs): + """ + Resample an ndarray from ``old_fs`` to ``new_fs`` along the last axis. + + Parameters + ---------- + data : numpy array + Input data to be resampled. + old_fs : int + Original sampling rate. + new_fs : int + New sampling rate. + + Returns + ------- + numpy array + """ + ndim = data.ndim + + if ndim == 1: + data = data[:, None] + elif ndim == 2: + data = data.T + else: + shape = data.shape + data = data.reshape(-1, data.shape[-1]).T + + resampled_data = soxr.resample(data, old_fs, new_fs) + + if ndim == 1: + resampled_data = resampled_data[:, 0] + elif ndim == 2: + resampled_data = resampled_data.T + else: + new_shape = shape[:-1] + resampled_data.shape[:1] + resampled_data = resampled_data.T.reshape(new_shape) + + return resampled_data diff --git a/requirements.txt b/requirements.txt index cbfe8fcf..bf8d42a7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,5 @@ matplotlib sounddevice samplerate mir_eval +python-sofa +soxr diff --git a/setup.py b/setup.py index d62fd050..c90d7ce2 100644 --- a/setup.py +++ b/setup.py @@ -180,6 +180,7 @@ def build_extensions(self): "pyroomacoustics.bss", "pyroomacoustics.denoise", "pyroomacoustics.phase", + "pyroomacoustics.simulation", ], # Libroom C extension ext_modules=ext_modules, @@ -190,6 +191,7 @@ def build_extensions(self): "numpy>=1.13.0", "scipy>=0.18.0", "pybind11>=2.2", + "soxr", ], cmdclass={"build_ext": BuildExt}, # taken from pybind11 example zip_safe=False,