Source directivity and various extensions to gpuRIR

README.md

@@ -14,6 +14,7 @@
   * [`beta_SabineEstimation`](#beta_sabineestimation)
   * [`att2t_SabineEstimator`](#att2t_sabineestimator)
   * [`t2n`](#t2n)
+  * [`extensions`](#t2n)
 - [References](#references)
 
 
@@ -177,6 +178,13 @@ Estimation of the number of images needed for a correct RIR simulation.
 *3 elements list of integers.*
         The number of images sources to compute in each dimension.
 
+### `extensions`
+
+
+[@corrooli](https://github.com/corrooli) and [@fuerbringer](https://github.com/fuerbringer) made an awesome job implementing some interesting extra features for gpuRIR in python, such as simulating air absorption or frequency dependant wall absorption coefficients. 
+
+I would like to better integrate these features with the rest of the library, but sadly I have no time to do this now. You can find more information about how to use these functions in the [documentation of the subpackage](gpuRIR/extensions), in their [pull request](https://github.com/DavidDiazGuerra/gpuRIR/pull/27), or their [examples](examples/extensions). I am including this as they developed it and I cannot offer support for it.
+
 
 ## References
 

examples/extensions/mono_filters.py

@@ -0,0 +1,102 @@
+import numpy as np
+
+import gpuRIR.extensions.room_parameters as rp
+import gpuRIR.extensions.generate_RIR as gRIR
+import gpuRIR.extensions.filtering as filtering
+
+from gpuRIR.extensions.wall_absorption.materials import Materials as mat
+import gpuRIR.extensions.wall_absorption.freq_dep_abs_coeff as fdac
+
+from gpuRIR.extensions.filters.air_absorption_bandpass import AirAbsBandpass
+from gpuRIR.extensions.filters.air_absorption_stft import AirAbsSTFT
+from gpuRIR.extensions.filters.characteristic_filter import CharacteristicFilter
+from gpuRIR.extensions.filters import characteristic_models as model
+from gpuRIR.extensions.filters.linear_filter import LinearFilter
+
+''' Example script for applying various gpuRIR extensions on an RIR generated by gpuRIR.
+Optionally writes a mono impulse response (IR) wave file for convolving signals yourself, or convolves a sound file you provide.
+
+Usable extensions in this context are:
+
+*   Frequency dependent absorption coefficients (virtual room materials) -> freq_dep_abs_coeff and wall_materials
+*   Air absorption (via Bandpassing or STFT) -> AirAbsSTFT() or AirAbsBandpass()
+*   Source and receiver characteristics (virtual mic / speaker models -> CharacteristicFilter() and LinearFilter()
+*   Benchmarking -> verbose
+
+See comments in code for usage examples
+'''
+# Write impulse response (IR) wave file for convolving signals yourself
+write_IR_file = False
+
+# Enables adaptive gain. Raises volume without clipping. If disabled, the IR wave file amplitude may be too small to be audible.
+enable_adaptive_gain = False
+
+# Visualizes waveform and spectrogram of each generated IR file. Depending on filter, additional graphs are drawn.
+visualize = False
+
+# Prints calculation times and parameter/processing info onto the terminal if True. Needed for benchmarking, debugging and further info.
+verbose = False
+
+# If True, apply frequency dependent wall absorption coefficients to simulate realistic wall/ceiling/floor materials.
+# Caution: Needs more resources!
+freq_dep_abs_coeff = False
+
+# Wall, floor and ceiling materials the room is consisting of
+# Structure: Array of six materials (use 'mat.xxx') corresponding to:
+# Left wall | Right wall | Front wall | Back wall | Floor | Ceiling
+wall_materials = 4 * [mat.wallpaper_on_lime_cement_plaster] + \
+    [mat.parquet_glued] + [mat.concrete]
+
+# Define room parameters
+params = rp.RoomParameters(
+    room_sz=[5, 4, 3],  # Size of the room [m]
+    pos_src=[[3, 3,  1.8]],  # Positions of the sources [m]
+    pos_rcv=[[1.5, 1.5, 1.6]],  # Positions of the receivers [m]
+    orV_src=[0, -1, 0],  # Steering vector of source(s)
+    orV_rcv=[0.1, 1, 0],  # Steering vector of receiver(s)
+    spkr_pattern="omni",  # Source polar pattern
+    mic_pattern="card",  # Receiver polar pattern
+    T60=1.0,  # Time for the RIR to reach 60dB of attenuation [s]
+    # Attenuation when start using the diffuse reverberation model [dB]
+    att_diff=15.0,
+    att_max=60.0,  # Attenuation at the end of the simulation [dB]
+    fs=44100,  # Sampling frequency [Hz]
+    # Bit depth of WAV file. Either np.int8 for 8 bit, np.int16 for 16 bit or np.int32 for 32 bit
+    bit_depth=np.int32,
+    # Absorption coefficient of walls, ceiling and floor.
+    wall_materials=wall_materials
+)
+
+# Generate room impulse response (RIR) with given parameters
+if freq_dep_abs_coeff:
+    receiver_channels = fdac.generate_RIR_freq_dep_walls(
+        params, LR=False, order=10, band_width=100, factor=1.1, use_bandpass_on_borders=False, visualize=visualize, verbose=verbose)
+    filename_appendix = "freqdep"
+else:
+    receiver_channels = gRIR.generate_RIR(params)
+    filename_appendix = ""
+
+for i in range(len(params.pos_rcv)):
+    # All listed filters wil be applied in that order.
+    # Leave filters array empty if no filters should be applied.
+
+    filters = [
+        # Speaker simulation.
+        # Comment either one out
+        # CharacteristicFilter(model.tiny_speaker, params.fs, visualize=visualize),
+        # LinearFilter(101, (0, 100, 150, 7000, 7001, params.fs/2), (0, 0, 1, 1, 0, 0), params.fs),
+
+        # Air absorption simulation.
+        # Comment either one out
+        # AirAbsBandpass(divisions=50, max_frequency=params.fs/2, LR=True, order=10, use_bandpass_on_borders=True, verbose=verbose, visualize=visualize),
+        # AirAbsSTFT(),
+
+        # Mic simulation.
+        # Comment either one out
+        # CharacteristicFilter(model.sm57, params.fs, visualize=visualize),
+        # LinearFilter(101, (0, 100, 150, 7000, 7001, params.fs/2), (0, 0, 1, 1, 0, 0), params.fs, visualize=visualize)
+    ]
+
+    # The filters above are applied on the RIR generated by gpuRIR and returned into filtered_IR
+    filtered_IR = filtering.filter_mono_IR(
+        receiver_channels[i], filters, params.bit_depth, params.fs, filename_appendix=filename_appendix, write_wave=write_IR_file, enable_adaptive_gain=enable_adaptive_gain, visualize=visualize, verbose=verbose)

examples/extensions/stereo_filters.py

@@ -0,0 +1,170 @@
+import numpy as np
+
+import gpuRIR.extensions.room_parameters as rp
+import gpuRIR.extensions.generate_RIR as gRIR
+import gpuRIR.extensions.filtering as filtering
+
+from gpuRIR.extensions.wall_absorption.materials import Materials as mat
+import gpuRIR.extensions.wall_absorption.freq_dep_abs_coeff as fdac
+
+from gpuRIR.extensions.filters.air_absorption_bandpass import AirAbsBandpass
+from gpuRIR.extensions.filters.air_absorption_stft import AirAbsSTFT
+from gpuRIR.extensions.filters.characteristic_filter import CharacteristicFilter
+from gpuRIR.extensions.filters import characteristic_models as model
+from gpuRIR.extensions.filters.linear_filter import LinearFilter
+from gpuRIR.extensions.filters.hrtf_filter import HRTF_Filter
+
+from gpuRIR.extensions.hrtf.hrtf_binaural_receiver import BinauralReceiver
+
+''' Example script for applying various gpuRIR extensions on an RIR generated by gpuRIR.
+Optionally writes a stereo impulse response (IR) wave file for convolving signals yourself, or convolves a sound file you provide.
+
+Please note:
+*   Number of receivers is restricted to 2 for each ear due to HRTF.
+
+Extensions are:
+*   Head-related transfer function (HRTF, simulated human hearing)
+*   Frequency dependent absorption coefficients (virtual room materials) -> freq_dep_abs_coeff and wall_materials
+*   Air absorption (via Bandpassing or STFT) -> AirAbsSTFT() or AirAbsBandpass()
+*   Source and receiver characteristics (virtual mic / speaker models -> CharacteristicFilter()
+*   Binaural receiver (two receivers in room) -> BinauralReceiver()
+*   Head-related transfer function (simulated human hearing) -> use_hrtf
+*   Benchmarking -> verbose
+
+See comments in code for usage examples.
+'''
+
+# Write impulse response (IR) wave file for convolving signals yourself
+write_IR_file = False
+
+# Enables adaptive gain. Raises volume without clipping. If disabled, the IR wave file amplitude may be too small to be audible.
+enable_adaptive_gain = False
+
+# Visualizes waveform and spectrogram of each generated IR file. Depending on filter, additional graphs are drawn.
+visualize = False
+
+# Prints calculation times and parameter/processing info onto the terminal if True. Needed for benchmarking, debugging and further info.
+verbose = False
+
+# If True, apply frequency dependent wall absorption coefficients (virtual room materials) to simulate realistic wall/ceiling/floor materials.
+# Caution: Needs more resources!
+freq_dep_abs_coeff = False
+
+# Uses binaural (HRTF) processing to enable realistic sound in 3D spaces. Use of headphones recommended.
+use_hrtf = False
+
+# Wall, floor and ceiling materials the room is consisting of
+# Structure: Array of six materials (use 'mat.xxx') corresponding to:
+# Left wall | Right wall | Front wall | Back wall | Floor | Ceiling
+wall_materials = 4 * [mat.wallpaper_on_lime_cement_plaster] + \
+    [mat.parquet_glued] + [mat.concrete]
+
+# Setup of binaural receiver with head position [m] and head direction [m].
+head = BinauralReceiver(
+    head_position=[1.5, 1.5, 1.6], head_direction=[0, -1, 0], verbose=verbose)
+
+
+# Common gpuRIR parameters (applied to both channels)
+room_sz = [5, 4, 3]  # Size of the room [m]
+pos_src = [[1.5, 1.8, 1.8]]  # Positions of the sources [m]
+orV_src = [0, -1, 0]  # Steering vector of source(s)
+spkr_pattern = "omni"  # Source polar pattern
+mic_pattern = "homni"  # Receiver polar patterny
+T60 = 1.0  # Time for the RIR to reach 60dB of attenuation [s]
+# Attenuation when start using the diffuse reverberation model [dB]
+att_diff = 15.0
+att_max = 60.0  # Attenuation at the end of the simulation [dB]
+fs = 44100  # Sampling frequency [Hz]
+# Bit depth of WAV file. Either np.int8 for 8 bit, np.int16 for 16 bit or np.int32 for 32 bit
+bit_depth = np.int32
+
+if visualize:
+    head.visualize(room_sz, pos_src, orV_src)
+
+# Define room parameters
+params_left = rp.RoomParameters(
+    room_sz=room_sz,
+    pos_src=pos_src,
+    orV_src=orV_src,
+    spkr_pattern=spkr_pattern,
+    mic_pattern=mic_pattern,
+    T60=T60,
+    att_diff=att_diff,
+    att_max=att_max,
+    fs=fs,
+    bit_depth=bit_depth,
+    wall_materials=wall_materials,
+
+    # Positions of the receivers [m]
+    pos_rcv=[head.ear_position_l],  # Position of left ear
+    orV_rcv=head.ear_direction_l,  # Steering vector of left ear
+    head_direction=head.direction,
+    head_position=head.position
+)
+
+params_right = rp.RoomParameters(
+    room_sz=room_sz,
+    pos_src=pos_src,
+    orV_src=orV_src,
+    spkr_pattern=spkr_pattern,
+    mic_pattern=mic_pattern,
+    T60=T60,
+    att_diff=att_diff,
+    att_max=att_max,
+    fs=fs,
+    bit_depth=bit_depth,
+    wall_materials=wall_materials,
+
+    # Positions of the receivers [m]
+    pos_rcv=[head.ear_position_r],  # Position of right ear
+    orV_rcv=head.ear_direction_r,  # Steering vector of right ear
+    head_direction=head.direction,
+    head_position=head.position
+)
+
+# Generate two room impulse responses (RIR) with given parameters for each ear
+if freq_dep_abs_coeff:
+    receiver_channel_r = fdac.generate_RIR_freq_dep_walls(
+        params_right, LR=False, order=10, band_width=100, factor=1.1, use_bandpass_on_borders=False, visualize=visualize, verbose=verbose)
+    receiver_channel_l = fdac.generate_RIR_freq_dep_walls(
+        params_left, LR=False, order=10, band_width=100, factor=1.1, use_bandpass_on_borders=False, visualize=visualize, verbose=verbose
+    )
+    filename_appendix = "freqdep"
+
+else:
+    receiver_channel_r = gRIR.generate_RIR(params_right)
+    receiver_channel_l = gRIR.generate_RIR(params_left)
+    filename_appendix = ""
+
+# Common filters, applied to both channels.
+# All listed filters wil be applied in that order.
+# Leave filters array empty if no filters should be applied.
+filters_both = [
+    # Speaker simulation.
+    # Comment either one out
+    # CharacteristicFilter(model.tiny_speaker, fs, visualize=visualize),
+    # LinearFilter(101, (0, 100, 150, 7000, 7001, fs/2), (0, 0, 1, 1, 0, 0), fs),
+
+    # Air absorption simulation.
+    # Comment either one out
+    # AirAbsBandpass(divisions=50, max_frequency=fs/2, LR=True, order=10, use_bandpass_on_borders=False, verbose=verbose, visualize=visualize),
+    # AirAbsSTFT(),
+
+    # Mic simulation.
+    # Comment either one out
+    # CharacteristicFilter(model.sm57, fs, visualize=visualize),
+    # LinearFilter(101, (0, 100, 150, 7000, 7001, fs/2), (0, 0, 1, 1, 0, 0), fs, visualize=visualize)
+]
+
+if use_hrtf:
+    filters_r = filters_both + \
+        [HRTF_Filter('r', params_right, verbose=verbose)]
+    filters_l = filters_both + [HRTF_Filter('l', params_left, verbose=verbose)]
+else:
+    filters_r = filters_both
+    filters_l = filters_both
+
+# The filters above are applied on the RIR generated by gpuRIR and returned into filtered_IR_right and filtered_IR_left for both stereo channels
+filtered_IR_right, filtered_IR_left = filtering.filter_stereo_IR(receiver_channel_r[0], receiver_channel_l[0], filters_r, filters_l,
+                                                                 bit_depth, fs, filename_appendix=filename_appendix, write_wave=write_IR_file,
+                                                                 enable_adaptive_gain=enable_adaptive_gain, verbose=verbose, visualize=visualize)

examples/polar_plots.py

@@ -9,16 +9,20 @@
 import plotly.graph_objects as go
 from plotly.subplots import make_subplots
 
+import gpuRIR.extensions.filters.air_absorption_calculation as aa
+import gpuRIR.extensions.room_parameters as rp
+import gpuRIR.extensions.generate_RIR as gRIR
+
 """ 
 Visualizes polar patterns of the source by rotating source by 360, repeatedly calling gpuRIR.
-Warning: Takes up to 2-3 minutes depending on your hardware!
+Warning: Takes up to 2-3 minutes depending on your hardware!
 """
 
 # Feel free to change these parameters
 # Resolution of polar plot (amount to divide 360 degrees into)
 PARTITIONS = 360
 
-# gpuRIR parameters
+# gpuRIR parameters
 gpuRIR.activateMixedPrecision(False)
 gpuRIR.activateLUT(False)
 
@@ -31,7 +35,7 @@
 def normalize_amps(amps):
     '''
     Normalizing amplitude. Change all amplitude such that loudest sample has the amplitude of 1.
-    :param amps Array of samples.
+    :param amps Array of samples.
     '''
     amps_normalized = np.copy(amps)
     max_value = np.max(amps_normalized)
@@ -48,7 +52,7 @@ def create_polar_plot(fig, i, amps, name):
     '''
     Creates single polar plot.
 
-    :param fig Plot figure.
+    :param fig Plot figure.
     :param i Index of plot
     :param name Name of polar pattern
     '''
@@ -67,7 +71,7 @@ def create_polar_plots(title=""):
     '''
     Creates compilation of polar plots using plotly.
 
-    :param title Title of plot.
+    :param title Title of plot.
     '''
     fig_polar = make_subplots(rows=2, cols=3, specs=[[{'type': 'polar'}]*3]*2)
     for i in range(6):
@@ -120,29 +124,30 @@ def create_polar_plots(title=""):
             rad = degree*np.pi/180
 
             # RIR parameters
-            room_sz=[16, 8, 3]  # Size of the room [m]
-            pos_src=np.array([[4, 4, 1.7]])  # Positions of the sources [m]
-            pos_rcv=np.array([[10, 4, 2]])  # Positions of the receivers [m]
-            # Steering vector of source(s)
-            orV_src=np.matlib.repmat(np.array([np.cos(rad), np.sin(rad), 0]), 1, 1)
-            # Steering vector of receiver(s)
-            orV_rcv=np.matlib.repmat(np.array([1, 0, 0]), 1, 1)
-            spkr_pattern=POLAR_PATTERNS[p]  # Source polar pattern
-            mic_pattern="omni"  # Receiver polar pattern
-            T60=0.21  # Time for the RIR to reach 60dB of attenuation [s]
-            # Attenuation when start using the diffuse reverberation model [dB]
-            att_diff=15.0
-            att_max=60.0  # Attenuation at the end of the simulation [dB]
-            fs=44100  # Sampling frequency [Hz]
-            # Bit depth of WAV file. Either np.int8 for 8 bit, np.int16 for 16 bit or np.int32 for 32 bit
-            bit_depth=np.int32
-            beta=6*[0.1] # Reflection coefficients
-            Tdiff= gpuRIR.att2t_SabineEstimator(att_diff, T60) # Time to start the diffuse reverberation model [s]
-            Tmax = gpuRIR.att2t_SabineEstimator(att_max, T60)	 # Time to stop the simulation [s]
-            nb_img = gpuRIR.t2n( Tdiff, room_sz )	# Number of image sources in each dimension
+            params = rp.RoomParameters(
+                room_sz=[16, 8, 3],  # Size of the room [m]
+                pos_src=[[4, 4, 1.7]],  # Positions of the sources [m]
+                pos_rcv=[[10, 4, 2]],  # Positions of the receivers [m]
+                # Steering vector of source(s)
+                orV_src=np.matlib.repmat(
+                    np.array([np.cos(rad), np.sin(rad), 0]), 1, 1),
+                # Steering vector of receiver(s)
+                orV_rcv=np.matlib.repmat(np.array([1, 0, 0]), 1, 1),
+                spkr_pattern=POLAR_PATTERNS[p],  # Source polar pattern
+                mic_pattern="omni",  # Receiver polar pattern
+                T60=0.21,  # Time for the RIR to reach 60dB of attenuation [s]
+                # Attenuation when start using the diffuse reverberation model [dB]
+                att_diff=15.0,
+                att_max=60.0,  # Attenuation at the end of the simulation [dB]
+                fs=44100,  # Sampling frequency [Hz]
+                # Bit depth of WAV file. Either np.int8 for 8 bit, np.int16 for 16 bit or np.int32 for 32 bit
+                bit_depth=np.int32,
+                beta=6*[0.1]
+            )
 
             # Prepare sound data arrays.
-            receiver_channels = gpuRIR.simulateRIR(room_sz, beta, pos_src, pos_rcv, nb_img, Tmax, fs, Tdiff=Tdiff, orV_rcv=orV_rcv, mic_pattern=mic_pattern, orV_src=orV_src, spkr_pattern=spkr_pattern)
+            receiver_channels = gRIR.generate_RIR(params)
+
 
             # Stack array vertically
             impulseResponseArray = np.vstack(receiver_channels[0])