diff --git a/README.md b/README.md index 1b5e80e..3e70c4c 100644 --- a/README.md +++ b/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 diff --git a/examples/extensions/mono_filters.py b/examples/extensions/mono_filters.py new file mode 100644 index 0000000..51a9d1c --- /dev/null +++ b/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) diff --git a/examples/extensions/stereo_filters.py b/examples/extensions/stereo_filters.py new file mode 100644 index 0000000..8469a12 --- /dev/null +++ b/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) diff --git a/examples/polar_plots.py b/examples/polar_plots.py index 28e29b7..3a78018 100644 --- a/examples/polar_plots.py +++ b/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]) diff --git a/gpuRIR/__init__.py b/gpuRIR/__init__.py index 9d21ef0..119f5db 100644 --- a/gpuRIR/__init__.py +++ b/gpuRIR/__init__.py @@ -4,11 +4,10 @@ import numpy as np from scipy.optimize import minimize -from scipy.signal import convolve from gpuRIR_bind import gpuRIR_bind -__all__ = ["polar_patterns", "beta_SabineEstimation", "att2t_SabineEstimator", "t2n", "simulateRIR", "simulateTrajectory", "activate_mixed_precision", "activate_lut"] +__all__ = ["polar_patterns", "beta_SabineEstimation", "att2t_SabineEstimator", "t2n", "simulateRIR", "simulateTrajectory", "activateMixedPrecision", "activateLUT"] polar_patterns = { "omni": 0, diff --git a/gpuRIR/extensions/README.md b/gpuRIR/extensions/README.md new file mode 100644 index 0000000..8d9fa89 --- /dev/null +++ b/gpuRIR/extensions/README.md @@ -0,0 +1,148 @@ +# gpuRIR extensions +Expands the functionality of gpuRIR by multiple optional extensions, listed as follows: + +* Frequency dependent absorption coefficients (virtual room materials) +* Air absorption (via Bandpassing or STFT) +* Source and receiver characteristics (virtual mic / speaker models +* Binaural receiver (two receivers in room) +* Head-related transfer function (HRTF, simulated human hearing) +* Room parameter class for gpuRIR (consolidating and automatically calculating parameters) +* Easy impulse response (IR) file generation with stereo support + +We have example usages of these features in our example scripts: +`examples/extensions/mono_filters.py` and `examples/extensions/stereo_filters.py`. + +## Frequency dependent absorption coefficients + +Models the impact different wall materials have on the reflective properties of a room. Depending on the material, different frequency bands are absorbed more or less. More materials can be defined manually in `wall_absorption/materials.py`. + +We divide up the RIR into frequency bands using bandpass filters with high- and +lowpass filters on the edges. The frequency bands in the lower frequencies are +narrower in order to apply wall absorption data accurately in those lower +registers. This is because the source data is defined up to 4kHz in our case. + +### Usage +Set `freq_dep_abs_coeff` to `True` in order to use this feature. +Add six materials (`mat.name_of_material`) to `wall_materials` array corresponding to left wall, right wall, front wall, back wall, floor, ceiling. + +### Relevant scripts + +* `gpuRIR/extensions/wall_absorption/freq_dep_abs_coeff.py` +* `gpuRIR/extensions/wall_absorption/materials.py` + +### Parameters +* **params** gpuRIR parameter collection. +* **band_width** Starting frequency width of band. The smaller the value the less performant but more realistic the algorithm is. Must be greater than 1. +* **factor** By what factor the starting `band_width` is multiplied each round. The smaller the value, the less performant but more realistic the algorithm is. Must be greater than 1. +* **order** Order of LR or Butterworth filter. If LR is used, the order is internally halved and 2 bandpass filters are used. +* **LR** If true, uses Linkwitz-Riley (LR) filtering. If false, use Butterworth filtering. +* **use_bandpass_on_borders** Enables bandpass filters on the lowest and highest frequency band instead of using high- and lowpass filters. +* **visualize** Plot the band borders and spectrogram. +* **verbose** Print debugging information. + +# Filters +### Usage of filters +Instantiate new filters inside the `filters` array on the `mono_filters.py` example script. Some commented-out code are placed there as a guide. All filters have pre-defined parameters as standard values, but can be overridden with user values. + +The filters are applied in the order the user provided, topmost filters are applied first. + +The array can be left empty if no filters are to be applied. + +**Example** +Following is a simple example with a tiny simulated speaker as source, air absorption and a simulated microphone as receiver. +``` +filters = [ + CharacteristicFilter(model.tiny_speaker) + AirAbsBandpass(), + CharacteristicFilter(model.sm57_freq_response, params.fs) + ] +``` + +## Air absorption +Applying air absorption using bandpass or STFT (Short Time Fourier Transformation) methods. + +### Relevant scripts + +* `gpuRIR/extensions/filters/air_absorption_bandpass.py` +* `gpuRIR/extensions/filters/air_absorption_calculation.py` +* `gpuRIR/extensions/filters/air_absorption_stft.py` + +### Parameters +* **f** Frequency of pure tone. +* **T** Temperature(degrees Celsius) +* **hr** Relative humidity +* **ps** Atmospheric pressure ratio. + +### Bandpass specific parameters +* **max_frequency** Where the last band stops filtering (band- or highpass). +* **min_frequency** Where the first band starts filtering (band- or lowpass). +* **divisions** Into how many bands to separate the spectrum between min- and max frequency. +* **fs** Sampling rate (must match source file). +* **order** Butterworth order +* **LR** Enables Linkwitz-Riley filter. +* **use_bandpass_on_borders** Uses bandpass instead of high/lowpass filters on lowest and highest band. +* **verbose** Terminal output for debugging or further information. +* **visualize** Plots band divisions in a graph. + +### STFT specific parameters +* **fs** Sampling rate (must match source file) +* **nFFT** Length of fourier transformations used in number of samples +* **noverlap** How much overlap in samples the segments have +* **window** Determines how big the Hann window is. + +## CharacteristicFilter: Source and receiver characteristics +Allows the simulation of microphone and speaker frequency responses, for example using a small bluetooth speaker as a source and a phone microphone as the receiver. +Uses 1D interpolation to create models, and STFT (Short Time Fourier Transformation) to apply the modelled and interpolated characteristics (defined in `filters/characteristic_models.py`). More microphone and speaker models can be defined in `filters/characteristic_models.py` by translating a standard frequency response graph into an array. + +### Relevant scripts + +* `gpuRIR/extensions/filters/characteristic_filter.py` +* `gpuRIR/extensions/filters/characteristic_models.py` + +### Parameters +* **model** Model of speaker/microphone (examples defined in *filters/characteristic_models.py*) +* **fs** Sampling rate (must match source file) +* **nFFT** Length of fourier transformations used in number of samples +* **noverlap** How much overlap in samples the segments have +* **window** Determines how big the Hann window is. +* **visualize** Plots the interpolated frequency response to aid comparison to source data. + +## Linear Filter + +Allows the simulation of microphone and speaker frequency responses, but using discrete linear filters. + +### Relevant scripts + +* `gpuRIR/extensions/filters/linear_filter.py` + +### Parameters +* **numtabs** Length of the filter. +* **bands** Number of frequency bands. +* **desired** Desired behavior of filter. +* **fs** Sampling rate (must match source file). +* **visualize** Plots the filter. + +More info on those parameters can be found in the SciPy documentation on [firls](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firls.html). + +## Head-related transfer function (HRTF) + +We use the [CIPIC +database](https://www.ece.ucdavis.edu/cipic/spatial-sound/hrtf-data/) to apply +the effect of the HRTF to RIRs generated by gpuRIR. This is accomplished by +applying the HRIR corresponding to the azimuth and elevation of the direct path +to the source. The `HRTF_Filter` is best used with our `BinauralReceiver`, +which simulates human ears as two receivers placed in the room. + +Please see `examples/extensions/stereo_filters.py` for an example on how these two classes +are used. + +### Relevant scripts + +* `gpuRIR/extensions/filters/hrtf_filter.py` +* `gpuRIR/extensions/hrtf/hrtf_binaural_receiver.py` +* `gpuRIR/extensions/hrtf/hrtf_rir.py` + +### Parameters +* **channel** Which channel (side) to fetch from the database. Options are 'l' for left or 'r' for right. +* **params** An instance RoomParameters used for simulateRIR. Mostly used for head and source positions. +* **verbose** Enable verbose mode to print debugging info. diff --git a/gpuRIR/extensions/__init__.py b/gpuRIR/extensions/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gpuRIR/extensions/create_spectrogram.py b/gpuRIR/extensions/create_spectrogram.py new file mode 100644 index 0000000..34ceda6 --- /dev/null +++ b/gpuRIR/extensions/create_spectrogram.py @@ -0,0 +1,67 @@ +import sys +import numpy as np +import matplotlib.pyplot as plt +from scipy.io import wavfile +from scipy import signal + + +def create_spectrogram_from_file(inner_file_path, title="", channel_count=1): + """ Shows a spectrogram from a provided file. + + Parameters + ---------- + inner_file_path : str + Path to the sound file that we want to visualize (WAV format) + title : str, optional + Title of spectrogram + channel_count : int, optional + Number of channels + """ + fs, source = wavfile.read(inner_file_path) + for channel in range(0, channel_count): + channel_name = 'left' if channel == 0 else 'right' + x = source[:, channel] + plt.rcParams.update({'font.size': 18}) + f, t, Sxx = signal.spectrogram(x, fs, nfft=512) + plt.pcolormesh(t, f/1000, 10*np.log10(Sxx/Sxx.max()), + vmin=-80, vmax=0, cmap='inferno') + plt.ylabel('Frequency [kHz]') + plt.xlabel('Time [s]') + plt.title(f"{title} Channel: {channel_name}") + plt.colorbar(label='dB').ax.yaxis.set_label_position('left') + plt.show() + + +def create_spectrogram_from_data(source, fs, channel_name="", title=""): + """ Shows a spectrogram from a provided file. + + Parameters + ---------- + source : ndarray + Sound data that we want to visualize (stereo or mono) + fs : int + Sampling rate (Hertz) of sound data. + title : str, optional + Title of spectrogram + channel_count : int, optional + Number of channels + """ + plt.rcParams.update({'font.size': 18}) + f, t, Sxx = signal.spectrogram(source, fs, nfft=512) + plt.pcolormesh(t, f/1000, 10*np.log10(Sxx/Sxx.max()), + vmin=-80, vmax=0, cmap='inferno') + plt.ylabel('Frequency [kHz]') + plt.xlabel('Time [s]') + plt.title(f"{title} {channel_name}") + plt.colorbar(label='dB').ax.yaxis.set_label_position('left') + plt.show() + + +""" Can be run from a terminal. First argument is sound file path, second argument is diagram title. +""" +if len(sys.argv) > 1: + file_path = sys.argv[1] + if len(sys.argv) > 2: + create_spectrogram_from_file(file_path, sys.argv[2]) + else: + create_spectrogram_from_file(file_path, file_path) diff --git a/gpuRIR/extensions/filtering.py b/gpuRIR/extensions/filtering.py new file mode 100644 index 0000000..be6c687 --- /dev/null +++ b/gpuRIR/extensions/filtering.py @@ -0,0 +1,252 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.io import wavfile +import time + +from gpuRIR.extensions.filters.filter import Filter +from gpuRIR.extensions.create_spectrogram import create_spectrogram_from_data + + +""" Generates an impulse response WAV file (IR) with optional filters. +Example usage: Convolving (reverberating) an audio signal in an impulse response loader plug-in like Space Designer in Logic Pro X. +""" + + +def mono_adaptive_gain(source, bit_depth, ceiling): + ''' Increases amplitude (loudness) to defined ceiling. Operates in an adaptive manner to prevent clipping. + + Parameters: + ----------- + source : 2D ndarray + Sound data to process. + bit_depth : int + Bit depth of source sound data. + ceiling : int + Maximum loudness (relative dB, e.g. -1dB) the sound data should be amplified to. [dB] + This way, the loudest sample will result being the ceiling in dB (eg. -1dB) + + Returns + ------- + 2D ndarray + Amplified source sound data. + ''' + peak = np.max(source) + negative_peak = np.abs(np.min(source)) + + # Check if the negative or positive peak is of a higher magnitude + if peak < negative_peak: + peak = negative_peak + + max_gain = np.iinfo(bit_depth).max*10**(-ceiling/10) + factor = max_gain/peak + + return source * factor + + +def stereo_adaptive_gain(source_l, source_r, bit_depth, ceiling): + ''' Increases amplitude (loudness) to defined ceiling. + + Parameters: + ----------- + source_l : 2D ndarray + Left channel sound data to process. + source_r : 2D ndarray + Right channel sound data to process. + bit_depth : int + Bit depth of source sound data. + ceiling : float + Maximum loudness (relative dB, e.g. -1dB) the sound data should be amplified to. [dB] + This way, the loudest sample will result being the ceiling in dB (eg. -1dB) + + Returns + ------- + 2D ndarray + Amplified source sound data. + ''' + # Read out positive and negative peak values + peak_l = np.max(source_l) + peak_r = np.max(source_r) + negative_peak_l = np.abs(np.min(source_l)) + negative_peak_r = np.abs(np.min(source_r)) + + # Check if left or right negative or positive peak is of a higher magnitude + if peak_l < negative_peak_l: + peak_l = negative_peak_l + if peak_r < negative_peak_r: + peak_r = negative_peak_r + if peak_l < peak_r: + peak_l = peak_r + + # Calculate amplification factor based on bit depth and highest overall peak of both channels + max_gain = np.iinfo(bit_depth).max * 10 ** (-ceiling / 10) + factor = max_gain / peak_l + + return (source_l * factor, source_r * factor) + + +def filter_mono_IR(source, filters, bit_depth, fs, filename_appendix="", write_wave=False, enable_adaptive_gain=False, visualize=False, verbose=False): + ''' Filters a mono IR file out of given source sound data and an optional array of filters to be applied. + + Parameters: + ----------- + source : 2D ndarray + Sound data to be converted into an impulse response file. + filters : list + List of filters to be applied (in that order) + bit_depth : int + Bit depth of source sound data. + fs : int + Sampling rate [Hz] + filename_appendix : str, optional + Filename appendix to indicate which processing was done to the sound data before it reaches this method. + write_wave : bool, optional + Writes IR wave file for convolving signals yourself. + enable_adaptive_gain : bool, optional + Enables adaptive gain, amplifying the sound data to a defined ceiling value. + visualize : bool, optional + Enables waveform and spectrogram plots. + verbose : bool, optional + Terminal logging for benchmarking, debugging and further info. + + Returns + ------- + ndarray + Filtered IR sound data (array of samples) + ''' + # Prepare sound data arrays. + source_signal = np.copy(source) + + # Apply filters + for i in range(len(filters)): + start_time = time.time() + source_signal = Filter(filters[i]).apply(source_signal) + end_time = time.time() + + # Print processing time per filter (relevant for benchmarking) + if verbose: + print(f"{filters[i].NAME} time = {end_time-start_time} seconds") + filename_appendix = f"{filename_appendix}_{filters[i].NAME}" + + # Stack array vertically + IR_array = np.vstack(source_signal) + + # Increase Amplitude to usable levels + if enable_adaptive_gain: + IR_array = mono_adaptive_gain( + IR_array, bit_depth, 3) + + if write_wave: + # Create stereo file (dual mono) + IR_array_concatenated = np.concatenate((IR_array, IR_array), axis=1) + + # Write impulse response file + filename = f'IR_mono_{filename_appendix}_{time.time()}.wav' + wavfile.write(filename, fs, IR_array_concatenated.astype(bit_depth)) + + if visualize: + # Create spectrogram + create_spectrogram_from_data( + source_signal, fs, "Mono", filename_appendix) + + # Visualize waveform of IR + # plt.title(filename_appendix) + plt.plot(source_signal) + plt.show() + + return IR_array + + +def filter_stereo_IR(source_r, source_l, filters_r, filters_l, bit_depth, fs, filename_appendix="", write_wave=False, enable_adaptive_gain=False, verbose=False, visualize=False): + ''' Filters a stereo IR file out of given source sound data and an optional array of filters to be applied. + + Parameters: + ----------- + source_r : 2D ndarray + Right channel sound data to be converted into an impulse response file. + source_l : 2D ndarray + Left channel sound data to be converted into an impulse response file. + filters_r : list + List of right channel filters to be applied (in that order) + filters_l : list + List of left channel filters to be applied (in that order) + bit_depth : int + Bit depth of source sound data. + fs : int + Sampling rate [Hz] + filename_appendix : str, optional + Filename appendix to indicate which processing was done to the sound data before it reaches this method. + write_wave : bool, optional + Writes IR wave file for convolving signals yourself. + enable_adaptive_gain : bool, optional + Enables adaptive gain, amplifying the sound data to a defined ceiling value. + visualize : bool, optional + Enables waveform and spectrogram plots. + verbose : bool, optional + Terminal logging for benchmarking, debugging and further info. + + Returns + ------- + ndarray + Right channel filtered IR sound data (array of samples) + ndarray + Left channel filtered IR sound data (array of samples) + ''' + + # Prepare stereo sound data arrays. + source_signal_r = np.copy(source_r) + source_signal_l = np.copy(source_l) + + # Apply filters for both stereo channels + for i in range(len(filters_r)): + start_time = time.time() + source_signal_r = Filter(filters_r[i]).apply(source_signal_r) + end_time = time.time() + # Print processing time per filter (relevant for benchmarking) + if verbose: + print( + f"Right Channel {filters_r[i].NAME} time = {end_time-start_time} seconds") + filename_appendix = f"{filename_appendix}_{filters_r[i].NAME}" + + for i in range(len(filters_l)): + start_time = time.time() + source_signal_l = Filter(filters_l[i]).apply(source_signal_l) + end_time = time.time() + # Print processing time per filter (relevant for benchmarking) + if verbose: + print( + f"Left Channel {filters_l[i].NAME} time = {end_time-start_time} seconds") + filename_appendix = f"{filename_appendix}_{filters_l[i].NAME}" + + # Stack array vertically + IR_array_r = np.vstack(source_signal_r) + IR_array_l = np.vstack(source_signal_l) + + # Increase Amplitude to usable levels + if enable_adaptive_gain: + IR_array_l, IR_array_r = stereo_adaptive_gain( + IR_array_l, IR_array_r, bit_depth, 3) + + # Put both stereo channels together + IR_array_concatenated = np.concatenate((IR_array_l, IR_array_r), axis=1) + + if write_wave: + # Write impulse response file + filename = f'IR_stereo_{filename_appendix}_{time.time()}.wav' + wavfile.write(filename, fs, IR_array_concatenated.astype(bit_depth)) + + if visualize: + # Create spectrogram + create_spectrogram_from_data( + source_signal_l, fs, "Left", filename_appendix) + create_spectrogram_from_data( + source_signal_r, fs, "Right", filename_appendix) + + # Visualize waveform of IR + plt.plot(source_signal_l, label="Left channel") + plt.title("Left channel") + plt.show() + plt.plot(source_signal_r, label="Right channel") + plt.title("Right channel") + plt.show() + + return IR_array_r, IR_array_l diff --git a/gpuRIR/extensions/filters/__init__.py b/gpuRIR/extensions/filters/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gpuRIR/extensions/filters/air_absorption_bandpass.py b/gpuRIR/extensions/filters/air_absorption_bandpass.py new file mode 100644 index 0000000..7cee075 --- /dev/null +++ b/gpuRIR/extensions/filters/air_absorption_bandpass.py @@ -0,0 +1,193 @@ +import numpy as np +import matplotlib.pyplot as plt +import multiprocessing + +from gpuRIR.extensions.filters.filter import FilterStrategy +from gpuRIR.extensions.filters.air_absorption_calculation import air_absorption +from gpuRIR.extensions.filters.butterworth import Butterworth + + +class AirAbsBandpass(FilterStrategy): + ''' Applies air absorption attenuation using bandpassing. + Splits the frequency bands into defined amount of divisions, applies air absorption and combines the bands back into one RIR. + ''' + + def __init__(self, max_frequency=22050, min_frequency=1, divisions=50, fs=44100, order=4, LR=False, use_bandpass_on_borders=False, verbose=False, visualize=False): + ''' Instantiates bandpass-driven air absorption. + + Parameters + ---------- + max_frequency : int + Top end of frequency range (High pass cut). + min_frequency : int + Bottom end of frequency range (Low pass cut). + divisions : int + Amount of bands the RIR gets divided into. The higher the value the better the quality, but with performance penalties. + fs : int + Sample rate [Hz] + order : int + Order of Butterworth or Linkwitz-Riley filter. + LR : bool, optional + Enables Linkwitz-Riley filter. + use_bandpass_on_borders : bool, optional + Uses bandpass instead of high/lowpass filters on lowest and highest band. + verbose : bool, optional + Terminal output for debugging or further information. + visualize : bool, optional + Plots band divisions in a graph + ''' + assert(order >= 4), "Order must be greater than 4!" + + self.min_frequency = min_frequency + self.divisions = divisions + self.fs = fs + self.max_frequency = max_frequency + self.order = order + self.LR = LR + if LR: + order = order / 2 + self.use_bandpass_on_borders = use_bandpass_on_borders + self.NAME = "bandpass_air_abs" + self.visualize = visualize + self.verbose = verbose + + @staticmethod + def distance_travelled(sample_number, sampling_frequency, c): + ''' Calculates how much distance the sound has travelled in meters. + + Parameters + ---------- + sample_number : int + How many samples have passed. + sampling_frequency : int + Sample rate [Hz] + c : float + Speed of sound. + + Returns + ------- + Time passed since first sample [s] + ''' + seconds_passed = sample_number * (sampling_frequency ** (-1)) + return (seconds_passed * c) # [m] + + def apply_single_band(self, IR, band_num, frequency_range): + ''' Applies low/high/bandpass filtering and air absorption attenuation onto a single frequency band. + Decides which kind of pass filtering based on the band number. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array.. + band_num : int + number of frequency band. + frequency_range : float + Frequency range of band [Hz] + ''' + + band = frequency_range / self.divisions + + # Upper ceiling of each band + band_max = (band * band_num) + + # Lower ceiling of each band and handling of edge case + if band_num == 1: + band_min = self.min_frequency + else: + band_min = (band * (band_num - 1)) + + # Calculating mean frequency of band which determines the attenuation. + band_mean = (band_max + band_min) / 2 + + if self.verbose: + print( + f"Band {band_num}:\tMin:{band_min}\tMean:{band_mean}\tMax:{band_max}\tBand width:{band_max - band_mean}") + + # Prepare + apply bandpass filter + if band_num == 1 and not self.use_bandpass_on_borders: + filtered_signal = Butterworth.apply_pass_filter( + IR, band_max, self.fs, 'lowpass', self.order*4, self.visualize + ) + if self.verbose: print(f"Lowpass at {band_max}") + elif band_num == self.divisions and not self.use_bandpass_on_borders: + filtered_signal = Butterworth.apply_pass_filter( + IR, band_min, self.fs, 'highpass', self.order*4, self.visualize + ) + if self.verbose: print(f"Highpass at {band_min}") + + else: + filtered_signal = Butterworth.apply_bandpass_filter( + IR, band_min, band_max, self.fs, self.order, self.LR, self.visualize) + + # Calculate air absorption coefficients + alpha, _, c, _ = air_absorption(band_mean) + + # Apply attenuation + for k in range(0, len(filtered_signal)): + distance = self.distance_travelled(k, self.fs, c) + attenuation = distance*alpha # [dB] + filtered_signal[k] *= 10**(-attenuation / 10) + + return filtered_signal + + def air_absorption_bandpass(self, IR): + ''' Creates a multi processing pool and calls methods to apply bandpass based air absorption. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed IR array. + ''' + if self.visualize: + plt.xscale('log') + plt.title('Air Absorption: Butterworth filter frequency response') + plt.xlabel('Frequency [Hz]') + plt.ylabel('Amplitude [dB]') + plt.ylim(bottom=-40) + plt.margins(0, 0.1) + plt.grid(which='both', axis='both') + + pool = multiprocessing.Pool() + + frequency_range = self.max_frequency - self.min_frequency + + if not self.visualize: + if self.verbose: + print("Processed bands are out of order due to multi-processing.") + # Divide frequency range into defined frequency bands + filtered_signals = pool.starmap(self.apply_single_band, + ((IR, j, frequency_range) + for j in range(1, self.divisions + 1)) + ) + else: + filtered_signals = [] + for j in range(1, self.divisions + 1): + filtered_signals.append( + self.apply_single_band(IR, j, frequency_range)) + + if self.visualize: + plt.show() + + arr = np.array(filtered_signals) + return np.add(0, arr.sum(axis=0)) + + def apply(self, IR): + ''' + Calls method to apply air absorption to RIR using bandpass filtering. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed IR array. + ''' + return self.air_absorption_bandpass(IR) diff --git a/gpuRIR/extensions/filters/air_absorption_calculation.py b/gpuRIR/extensions/filters/air_absorption_calculation.py new file mode 100644 index 0000000..733081a --- /dev/null +++ b/gpuRIR/extensions/filters/air_absorption_calculation.py @@ -0,0 +1,177 @@ +import numpy as np + + +def air_absorption(f=100, T=25, hr=50, ps=1): + ''' + air_absorption: Calculates sound absorption (attenuation) in humid air + + ********************************************************************** + + Description + + Sound absorption (attenuation) in humid air depends on frequency, + temperature, relative humidity, and atmospheric pressure. + + Calculates sound absorption (attenuation) in humid air using the ISO + standard and the Bass forumla. The sound absorption is output in + dB/meter. + + See appropriate input and output variables sections below for + more details. + + ********************************************************************** + Input Variables + + f=100 # (Hz) frequency of pure tone + # default is f=100 + + T=25 # (degrees Celsius) temperature + # default is T=25 + + hr=50 # (Percent) Relative Humidity + # default is hr=50 + + ps=1 # atmospheric pressure ratio + # pa/pr ratio of ambient atmospoheric pressure to the + # standard atmosphere. + # default is ps=1 + + ********************************************************************** + Output Variables + + alpha sound absorption is (dB/meter) according to Bass + + alpha_iso sound absorption is (dB/meter) according to ISO standard + + c speed of sound in humid air (meters/second) according to Bass + + c_iso speed of sound in humid air (meters/second) according to ISO standard + + ********************************************************************** + + Example='1' + + + f=100 # frequency in Hz + + T=20 # 20 Degrees Celsius + + hr=80 # Relative humidity in percentage hr=80 means 80 percent humidity + + ps=1 # Is the barometric pressure ratio. Usually, ps=1 + + # Run the program + [alpha, alpha_iso, c, c_iso]=air_absorption(f, T, hr, ps) + + + # ********************************************************************** + # + # References + # + # ANSI/ASA S1.26-1995 (R2009) + # + # ISO 9613-1:1993 Acoustics -- Attenuation of sound during propagation + # outdoors -- Part 1: Calculation of the absorption of sound by + # the atmosphere + # + # + # Need both Articles below inclduing the (Further Developments) article + # from Bass to make the calculations correctly. + # + # Attenborough, K., S. Taherzadeh, H. E. Bass, X. Di, R. Raspet, + # G. R. Becker, A. Güdesen et al. "Benchmark cases for outdoor sound + # propagation models." The Journal of the Acoustical Society of America + # 97 (1995): 173. + # + # Bass, H. E., L. C. Sutherland, A. J. Zuckerwar, D. T. Blackstock, + # and D. M. Hester. "Atmospheric absorption of sound: Further + # developments." The Journal of the Acoustical Society of America + # 97, no. 1 (1995): 680-683. + # + # Pierce A.D., Acoustics an introduciton to its physical principles and + # aplications, Acoustical Society of America, Equation 1-9.5 pp.30 + # + # ********************************************************************** + # + # This program was Written by Edward L. Zechmann + # + # date September 2006 + # + # modified March 2010 Updated Comments chnged input units of + # T from Farhenheit to Celsius. + # + # modified October 24 2013 Updated formula for speed of sound by + # dividing, h the humidity in percent + # molecular concentration by 100 to + # convert it to a fraction. + # + # ported October 6 2021 Ported from Matlab to Python + # + # ********************************************************************** + # + # Please Feel Free to Modify This Program + # + # See Also: Atmosphere (its on the Matlab file exchange) + + Edward Zechmann (2021). Air Absorption (https://www.mathworks.com/matlabcentral/fileexchange/26835-air-absorption), + MATLAB Central File Exchange. Retrieved December 11, 2021. + ''' + + # convert T from Celsius to Kelvin + # T=273.15+5/9*(T-32) to convert from Fahrenheit + T = 273.15 + T # To convert from Celsius + + # listing of constants + + T01 = 273.16 # triple point in degrees Kelvin + T0 = 293.15 + # atmospheric pressure ratio is the ambient pressure/standard pressure + ps0 = 1 # ps0= standard pressure/standard pressure which is unity + + # Bass formula for saturation pressure ratio + psat_ps0 = 10 ** (10.79586*(1-T01/T) - 5.02808*np.log10(T/T01) + 1.50474*10**(-4) * + (1-10**(-8.29692*(T/T01-1))) - 4.2873*10**(-4)*(1-10**(-4.76955*(T01/T-1)))-2.2195983) + + # Iso formula for saturation pressure ratio + psat_ps0_iso = 10 ** (-6.8346 * (T01 / T) ** 1.261 + 4.6151) + + ps_ps0 = ps / ps0 + + h = hr*psat_ps0 / ps_ps0 # h is the humidity in percent molar concentration + # h is the humidity in percent molar concentration + h_iso = hr*psat_ps0_iso / ps_ps0 + + c0 = 331 # c0 is the reference sound speed + c = (1+0.16*h/100) * c0 * np.sqrt(T/T01) + c_iso = (1+0.16*h_iso/100)*c0*np.sqrt(T/T01) + + # # ********************************************************************** + # + # Bass formula + + F = f/ps + + Fr0 = 1/ps0*(24+4.04*10**4*h*(0.02+h)/(0.391+h)) + FrN = 1/ps0*(T0/T)**(1/2)*(9+280*h*np.exp(-4.17*((T0/T)**(1/3)-1))) + + # Calculate the air absorption in dB/meter using the Bass formula + alpha = 20*np.log10(np.exp(1))*ps*F**2*(1.84*10**(-11)*(T/T0)**(0.5)*ps0+((T/T0)**(-5/2)) + * (0.01275*np.exp(-2239.1/T)/(Fr0+F**2/Fr0)+0.1068*np.exp(-3352/T)/(FrN+F**2/FrN))) + + # # ********************************************************************** + # + # ISO formula + + taur = T/T0 + pr = ps/ps0 + + fr0 = pr*(24+40400*h_iso*(0.02+h_iso)/(0.391+h_iso)) + frN = pr*(taur)**(-1/2)*(9+280*h_iso*np.exp(-4.17*((taur)**(-1/3)-1))) + + b1 = 0.1068*np.exp(-3352/T)/(frN+f**2/frN) + b2 = 0.01275*np.exp(-2239.1/T)/(fr0+f**2/fr0) + + # Calculate the air absorption in dB/meter for the ISO standard + alpha_iso = 8.686*f**2*taur**(1/2)*(1.84*10**(-11)/pr+taur**(-3)*(b1+b2)) + + return [alpha, alpha_iso, c, c_iso] diff --git a/gpuRIR/extensions/filters/air_absorption_stft.py b/gpuRIR/extensions/filters/air_absorption_stft.py new file mode 100644 index 0000000..960ea2f --- /dev/null +++ b/gpuRIR/extensions/filters/air_absorption_stft.py @@ -0,0 +1,81 @@ +from gpuRIR.extensions.filters.filter import FilterStrategy +import gpuRIR.extensions.filters.air_absorption_calculation as aa + +import numpy as np +from scipy.signal import istft, stft + + +class AirAbsSTFT(FilterStrategy): + ''' + Applies air absorption attenuation using short-time fourier transformation (STFT). + Splits the frequency bands into defined amount of divisions, applies air absorption and combines the bands back into one RIR. + + Reference: + Hamilton Brian. Adding air attenuation to simulated room impulse responses: A modal approach, 2021. URL: https://arxiv.org/pdf/2107.11871.pdf + ''' + + def __init__(self, T=25, hr=50, ps=1, fs=44100, nFFT=256, noverlap=int(256*0.75), window='hanning'): + ''' Instatiates STFT-driven air absorption. + For further reference about parameters see scipy.signal.stft + ''' + self.T = T + self.hr = hr + self.ps = ps + self.fs = fs + self.nFFT = nFFT + self.noverlap = noverlap + self.window = window + self.NAME = "STFT_air_abs" + + def STFT_air_absorption(self, RIR): + ''' Applies air absorption to RIR using STFT. + + Parameters + ---------- + RIR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + ''' + # Calculate STFT of RIR (like spectrogram) + f, t, RIR_TF = stft(RIR, self.fs, self.window, nperseg=self.nFFT, noverlap=self.noverlap, + nfft=self.nFFT, boundary='even', padded=True, return_onesided=False) + + # Get air absorption coeffs for given configuration, coeffs are in dB/meter + alphas_t0, _, c, _ = aa.air_absorption(f, self.T, self.hr, self.ps) + + # Get air absorption coeffs over distance/time for each band + alphas = np.zeros((len(f), len(t))) + for ii in range(len(t)): + alphas[:, ii] = -alphas_t0*t[ii]*c + + # Set lower bound for attenuation + alphas[alphas < -100] = -100 + + # Get linear absorption coeffs and apply them to the STFT of the RIR + RIR_TF_processed = RIR_TF*np.power(10, (alphas/20)) + + # Transform processed STFT of RIR back to time domain + _, RIR_processed = istft(RIR_TF_processed, self.fs, self.window, + nperseg=self.nFFT, noverlap=self.noverlap, nfft=self.nFFT) + RIR_processed = RIR_processed[0:len(RIR)] + return RIR_processed + + def apply(self, IR): + ''' + Calls method to apply air absorption to RIR using STFT. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + ''' + return self.STFT_air_absorption(IR) diff --git a/gpuRIR/extensions/filters/butterworth.py b/gpuRIR/extensions/filters/butterworth.py new file mode 100644 index 0000000..e2547c8 --- /dev/null +++ b/gpuRIR/extensions/filters/butterworth.py @@ -0,0 +1,125 @@ +from scipy.signal import butter, sosfreqz, sosfilt +import matplotlib.pyplot as plt +import numpy as np + + +class Butterworth(): + ''' Encapsulates Butterworth/Linkwitz Riley functionality with visualization. + ''' + + @staticmethod + def create_bandpass_filter(lowcut, highcut, fs, order, visualize=False): + ''' Returns a butterworth bandpass filter. + + Parameters + ---------- + lowcut : int + Lower bound of bandpass. [Hz] + highcut : int + Upper bound of bandpass. [Hz] + fs : int + Sample rate. [Hz] + order : int + Order of Butterworth filter. + visualize : bool, optional + Plots band divisions in a graph. + + Returns + ------- + ndarray + Second order sections of IIR filter. + ''' + sos = butter(order, [lowcut, highcut], + btype='bandpass', fs=fs, output='sos') + if visualize: + w, h = sosfreqz(sos, fs=fs) + plt.plot(w, 20*np.log10(np.abs(h)+1e-7)) + return sos + + @staticmethod + def apply_bandpass_filter(data, lowcut, highcut, fs, order, LR=False, visualize=False): + ''' Applies a butterworth bandpass filter. + + Parameters + ---------- + lowcut : int + Lower bound of bandpass. [Hz] + highcut : int + Upper bound of bandpass. [Hz] + fs : int + Sample rate. [Hz] + order : int + Order of Butterworth filter. + LR : bool + Enables Linkwitz-Riley filter. + visualize : bool, optional + Plots band divisions in a graph. + + Returns + ------- + ndarray + y Filtered sound data. + ''' + sos = Butterworth.create_bandpass_filter( + lowcut, highcut, fs, order, visualize) + y = sosfilt(sos, data) + if LR: + # Filter once again for Linkwitz-Riley filtering + y = sosfilt(sos, data) + return y + + @staticmethod + def create_pass_filter(cut, fs, pass_type, order, visualize=False): + ''' Returns a butterworth filter. + + Parameters + ---------- + cut : int + Cut frequency [Hz] + fs : int + Sample rate. [Hz] + pass_type : str + Type of butterworth filter (e.g. 'lowpass' or 'highpass'). + order : int + Order of Butterworth filter. + visualize : bool, optional + Plots band divisions in a graph. + + Returns + ------- + ndarray + Second order section of butterworth filter. + ''' + sos = butter(order, cut, btype=pass_type, fs=fs, output='sos') + + if visualize: + w, h = sosfreqz(sos, fs=fs) + plt.plot(w, 20*np.log10(np.abs(h)+1e-7)) + return sos + + @staticmethod + def apply_pass_filter(data, cut, fs, pass_type, order, visualize=False): + ''' Applies a butterworth filter. + + Parameters + ---------- + cut : int + Cut frequency [Hz] + fs : int + Sample rate. [Hz] + pass_type : str + Type of butterworth filter (e.g. 'lowpass' or 'highpass'). + order : int + Order of Butterworth filter. + visualize : bool, optional + Plots band divisions in a graph. + + Returns + ------- + ndarray + Second order section of butterworth filter. + ''' + sos = Butterworth.create_pass_filter( + cut, fs, pass_type, order=order, visualize=visualize) + sos = sosfilt(sos, data) # Single Filter + return sos diff --git a/gpuRIR/extensions/filters/characteristic_filter.py b/gpuRIR/extensions/filters/characteristic_filter.py new file mode 100644 index 0000000..150e20e --- /dev/null +++ b/gpuRIR/extensions/filters/characteristic_filter.py @@ -0,0 +1,134 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +from scipy.signal import istft, stft + +from gpuRIR.extensions.filters.filter import FilterStrategy + + +class CharacteristicFilter(FilterStrategy): + ''' Applies deviations of a frequency response (e.g. a microphone frequency response) to an audio signal via short time fourier transformation. + ''' + + def __init__(self, freq_response, T=25, hr=50, ps=1, fs=44100, nFFT=256, noverlap=int(256*0.75), window='hanning', visualize=False): + ''' Instantiates frequency response filtering based on predefined models using short-time fourier transformation (STFT). + For further reference about STFT parameters see scipy.signal.stft + + Parameters + ---------- + freq_response : 2d ndarray + Frequency response array (mapping from frequency to attenuation/amplification) with examples defined in extensions/filters/characteristic_models.py (np array shape (n, 2), n = number of rows) + visualize : bool, optional + Enables plotting of frequency response. + ''' + self.T = T + self.hr = hr + self.ps = ps + self.fs = fs + self.nFFT = nFFT + self.noverlap = noverlap + self.window = window + self.freq_response = freq_response + self.visualize = visualize + self.NAME = "characteristic" + + @staticmethod + def interpolate_frequency_response(freq_response, plot=False): + ''' Interpolates frequency response array with 1D interpolated curves generated from given frequency response model. + + Parameters + ---------- + freq_response : 2D ndarray + See __init__ from CharacteristicFilter. + plot : bool, optional + If True, plot frequency response graph. + + Returns + ------- + function + Interpolation function f. + ''' + # y: relative response [dB] + # x: frequency [Hz] + y = freq_response[:, 1] + x = freq_response[:, 0] + f = interp1d(x, y, fill_value=-100, bounds_error=False) + x_interpolated = np.arange(1, 20000) + y_interpolated = f(x_interpolated) + if plot: + plt.title('Characteristics Filter: Interpolated Frequency Response') + plt.plot(x, y, 'o', x_interpolated, y_interpolated, "-") + plt.show + return f + + def apply_characteristic(self, RIR): + ''' Applies characteristics filtering on the source data. + + Parameters + ---------- + RIRs : 2D ndarray + Room impulse responses generated by gpuRIR. + + Returns + ------- + Processed RIRs. + ''' + characteristic = self.interpolate_frequency_response( + self.freq_response, self.visualize) + + # Calculate STFT of RIR (like spectrogram) + f, t, RIR_TF = stft( + RIR, + self.fs, + self.window, + nperseg=self.nFFT, + noverlap=self.noverlap, + nfft=self.nFFT, + boundary='even', + padded=True, + return_onesided=False + ) + + alphas = np.zeros((len(f), len(t))) + + # Create array for relative gains based on interpolated frequency response + gains = characteristic(f) + + # Apply gain deviations based on gains array. + for i in range(len(t)): + alphas[:, i] = gains + + # Set lower bound for attenuation + alphas[alphas < -100] = -100 + + # Get calculated coeffs and apply them to the STFT of the RIR + RIR_TF_processed = RIR_TF*np.power(10, (alphas/20)) + + # Transform processed STFT of RIR back to time domain + _, RIR_processed = istft( + RIR_TF_processed, + self.fs, + self.window, + nperseg=self.nFFT, + noverlap=self.noverlap, + nfft=self.nFFT + ) + + RIR_processed = RIR_processed[0:len(RIR)] + return RIR_processed + + def apply(self, IR): + ''' Calls method to apply characteristics filtering on the source data. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + + ''' + return self.apply_characteristic(IR) diff --git a/gpuRIR/extensions/filters/characteristic_models.py b/gpuRIR/extensions/filters/characteristic_models.py new file mode 100644 index 0000000..f38d908 --- /dev/null +++ b/gpuRIR/extensions/filters/characteristic_models.py @@ -0,0 +1,76 @@ +import numpy as np + +''' +Microphone models +''' +# Shure SM57 dynamic microphone. Standard mic for US presidential speeches +sm57 = np.array([ + # Frequency in Hz | Relative response in dB + [50, -10], + [100, -4], + [200, 0], + [400, -1], + [700, -0.1], + [1000, 0], + [1500, 0.05], + [2000, 0.1], + [3000, 2], + [4000, 3], + [5000, 5], + [6000, 6.5], + [7000, 3], + [8000, 2.5], + [9000, 4], + [10000, 3.5], + [12000, 1], + [14000, 2], + [15000, -5] +]) + +# iPhone X. +iphone_x_speaker = np.array([ + [30, -12], + [40, -3], + [50, 0], + [70, 1], + [100, 2], + [150, 3], + [170, 2], + [200, 1.5], + [250, 1.5], + [300, 3], + [400, 2.5], + [500, 2], + [600, 2], + [1000, 2.5], + [1500, 3], + [2000, 3], + [2500, 4], + [3000, 4], + [4000, 3], + [5000, 3.5], + [6000, 3.8], + [7000, 4], + [8000, 4.5], + [9000, 4.75], + [10000, 5], + [12500, 1], + [15000, -5], + [17000, -7], + [20000, -8] +]) + +''' +Speaker models +''' +# Tiny speaker, comparable to smartphone speaker +tiny_speaker = np.array([ + [200, -10], + [300, -2], + [1000, 0], + [5000, 3], + [7000, 5], + [10000, -2], + [12000, -5], + [14000, -10] +]) diff --git a/gpuRIR/extensions/filters/filter.py b/gpuRIR/extensions/filters/filter.py new file mode 100644 index 0000000..97e1718 --- /dev/null +++ b/gpuRIR/extensions/filters/filter.py @@ -0,0 +1,34 @@ +class Filter(object): + ''' Encapsulates FilterStrategy and provides interface for applying filters on audio. + ''' + def __init__(self, filter_strategy): + ''' Initializes a filter. + + Parameters + ---------- + filter_strategy : FilterStrategy + Instance of a filter. + ''' + self._filter_strategy = filter_strategy + + def apply(self, IR): + ''' Applies a filter to given impulse response (IR). + + Parameters + ---------- + IR : 2D ndarray + Raw impulse response (IR) data. + + Returns + ------- + 2D ndarray + Filtered impulse response (IR) data. + ''' + return self._filter_strategy.apply(IR) + +class FilterStrategy(object): + ''' Superclass of all filters. + ''' + def apply(self): + ''' Required method + ''' \ No newline at end of file diff --git a/gpuRIR/extensions/filters/hrtf_filter.py b/gpuRIR/extensions/filters/hrtf_filter.py new file mode 100644 index 0000000..c2dc91f --- /dev/null +++ b/gpuRIR/extensions/filters/hrtf_filter.py @@ -0,0 +1,234 @@ +import numpy as np + +from gpuRIR.extensions.filters.filter import FilterStrategy +from gpuRIR.extensions.hrtf.hrtf_rir import HRTF_RIR + + +class HRTF_Filter(FilterStrategy): + ''' Head related transfer function to simulate psychoacoustic effects of a human head and upper body. + + Reference: Algazi V.R., Duda R.O., Thompson D.M. and Avendano C. The CIPIC + HRTF database. In Proceedings of the 2001 IEEE Workshop on the Applications of + Signal Processing to Audio and Acoustics (Cat. No.01TH8575) (2001), pp. 99–102. + ''' + # 90 degree angle in radiants + ANGLE_90 = np.pi/2 + + # 180 degree angle in radiants + ANGLE_180 = np.pi + + def __init__(self, channel, params, verbose=False): + ''' Initialized HRTF filter. + + Parameters + ---------- + channel : str + Determines if left or right channel is being processed. Either 'l' or 'r'. + params : RoomParameters + Abstracted gpuRIR parameter class object. + verbose : bool, optional + Terminal output for debugging or further information + ''' + self.channel = channel + self.NAME = "HRTF" + self.params = params + self.hrtf_rir = HRTF_RIR() + self.verbose = verbose + + @staticmethod + def find_angle(u, v): + ''' Find angle between two vectors on a 2D plane. + + Parameters + ---------- + u : ndarray + Vector with two elements + v : ndarray + Vector with two elements + + Returns + ------- + float + Scalar angle in radiants. + + ''' + norm_product = (np.linalg.norm(u) * np.linalg.norm(v)) + + if norm_product != 0: + return np.arccos((u @ v) / norm_product) + + return 0 + + # Find elevation between head and source + + @staticmethod + def calculate_elevation(pos_src, pos_rcv, head_direction): + ''' Calculates elevation between head position / direction and signal source position. + + Parameters + ---------- + pos_src : 3D ndarray + Position of signal source. + pos_rcv : 3D ndarray + Position of signal receiver (center of head). + head_direction : 3D ndarray + Direction in which the head is pointing towards. + + Returns + ------- + float + Elevation angle between head position / direction and signal source. + ''' + # Height of source + opposite = np.abs(pos_src[2] - pos_rcv[2]) + + # Length of floor distance between head and source + adjacent = np.linalg.norm( + np.array([pos_src[0], pos_src[1]]) - np.array([pos_rcv[0], pos_rcv[1]])) + + # Find elevation between head and source positions + if adjacent != 0: + el_rcv_src = np.arctan(opposite / adjacent) + else: + el_rcv_src = np.arctan(np.inf) + + # Edge case if source is below head + if pos_rcv[2] > pos_src[2]: + el_rcv_src = -el_rcv_src + + # Height of receiver + opposite = np.abs(head_direction[2]) + + # Length of floor distance between head and head direction vector + adjacent = np.linalg.norm(np.array([head_direction[0], head_direction[1]])) + + # Calculate elevation between head and head direction + el_rcv_dir = np.arctan(opposite / adjacent) + + # Edge case if source is below head + if pos_rcv[2] > pos_src[2]: + elevation_angle = el_rcv_src + el_rcv_dir + else: + elevation_angle = el_rcv_src - el_rcv_dir + + # Edge case if source is behind head + angle, _, _ = HRTF_Filter.vector_between_points( + pos_src, pos_rcv, head_direction) + if angle > HRTF_Filter.ANGLE_90: + # Source is behind head + elevation_angle = HRTF_Filter.ANGLE_180 - elevation_angle + + # Subtract elevation between head and source and between head and head direction + return elevation_angle + + @staticmethod + def vector_between_points(pos_src, pos_rcv, head_direction): + ''' Calculates a vector between two points in a 2D plane. + + Parameters + ---------- + pos_src : 3D ndarray + Position of signal source. + pos_rcv : 3D ndarray + Position of signal receiver (center of head). + head_direction : 3D ndarray + Direction in which the head is pointing towards. + + Returns + ------- + float + Scalar angle in radiants. + 2D ndarray + 2D Vector between head and signal source. + 2D ndarray + 2D vector of head direction. + + ''' + # 3D vector from head position (origin) to source + head_to_src = pos_src - pos_rcv + + # Extract 2D array from 3D + head_to_src = np.array([head_to_src[0], head_to_src[1]]) + # Extract 2D array from 3D + headdir_xy = [head_direction[0], head_direction[1]] + + # Return angle using trigonometry + return HRTF_Filter.find_angle(headdir_xy, head_to_src), head_to_src, headdir_xy + + @staticmethod + def calculate_azimuth(pos_src, pos_rcv, head_direction): + ''' Calculates azimuth between head position / direction and signal source position. + + Parameters + ---------- + pos_src : 3D ndarray + Position of signal source. + pos_rcv : 3D ndarray + Position of signal receiver (center of head). + head_direction : 3D ndarray + Direction in which the head is pointing towards. + + Returns + ------- + float + Azimuth angle between head position / direction and signal source. + ''' + # Find angle using trigonometry + angle, head_to_src, headdir_xy = HRTF_Filter.vector_between_points( + pos_src, pos_rcv, head_direction) + + # Check if azimuth goes above 90° + if angle > HRTF_Filter.ANGLE_90: + angle = np.pi - angle + + # Check left/right. If positive direction is left, if negative direction is right. + side = np.sign(np.linalg.det([headdir_xy, head_to_src])) + + return angle * (-side) + + def hrtf_convolve(self, IR): + ''' + Convolves an impulse response (IR) array with a HRTF room impulse response (RIR) retrieved from the CIPIC database. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + ''' + elevation = self.calculate_elevation( + self.params.pos_src[0], self.params.head_position, self.params.head_direction) + + if self.verbose: + print(f"Elevation = {elevation * (180 / np.pi)}") + + azimuth = self.calculate_azimuth( + self.params.pos_src[0], self.params.head_position, self.params.head_direction) + + if self.verbose: + print(f"Azimuth = {azimuth * (180 / np.pi)}") + + hrir_channel = self.hrtf_rir.get_hrtf_rir( + elevation, azimuth, self.channel) + + return np.convolve(IR, hrir_channel, mode='same') + + def apply(self, IR): + ''' Calls method to apply HRTF filtering on the source data. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + + ''' + return self.hrtf_convolve(IR) diff --git a/gpuRIR/extensions/filters/linear_filter.py b/gpuRIR/extensions/filters/linear_filter.py new file mode 100644 index 0000000..5d13c55 --- /dev/null +++ b/gpuRIR/extensions/filters/linear_filter.py @@ -0,0 +1,64 @@ +import scipy.signal +import numpy as np +import matplotlib.pyplot as plt + +from gpuRIR.extensions.filters.filter import FilterStrategy + +class LinearFilter(FilterStrategy): + + def __init__(self, numtaps, bands, desired, fs, visualize=False): + ''' Instantiates linear filtering. + For parameters reference see scipy.signal.firls. + ''' + self.numtaps=numtaps + self.bands = bands + self.desired = desired + self.fs = fs + self.visualize=visualize + self.NAME = "linear" + + + def apply_linear_filter(self, data): + ''' Applies a linear filter on given sound data. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + ''' + b_siumulated_mic = scipy.signal.firls(self.numtaps, self.bands, self.desired, fs=self.fs) # design filter + + if self.visualize: + # Plot frequency resonse of filter + w, response = scipy.signal.freqz(b_siumulated_mic) + freq = w/np.pi*self.fs/2 + _, ax1 = plt.subplots() + ax1.set_title('Linear Filter: Simulated Frequency Response') + ax1.plot(freq, 20 * np.log10(abs(response)), 'b') + ax1.set_ylabel('Amplitude [dB]', color='b') + ax1.set_xlabel('Frequency [Hz]') + plt.show() + + # Apply filter to audio signal + return scipy.signal.lfilter(b_siumulated_mic, 1, data) + + def apply(self, IR): + ''' Calls method to apply linear filtering on the source data. + + Parameters + ---------- + IR : 2D ndarray + Room impulse response array. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + + ''' + return self.apply_linear_filter(IR) diff --git a/gpuRIR/extensions/generate_RIR.py b/gpuRIR/extensions/generate_RIR.py new file mode 100644 index 0000000..9d207f3 --- /dev/null +++ b/gpuRIR/extensions/generate_RIR.py @@ -0,0 +1,27 @@ +import gpuRIR + +def generate_RIR(param): + # TODO: Integrate this with the original simulateRIR function preserving backward compatibility + ''' + Generates RIRs from the gpuRIR library. + + Parameters + ---------- + params : RoomParameters + gpuRIR parameters + + Returns + ------- + 2D ndarray + Receiver channels from room impulse response (RIR) generated by gpuRIR. + ''' + gpuRIR.activateMixedPrecision(False) + gpuRIR.activateLUT(False) + + RIRs = gpuRIR.simulateRIR(param.room_sz, param.beta, param.pos_src, param.pos_rcv, param.nb_img, + param.Tmax, param.fs, Tdiff=param.Tdiff, + orV_src=param.orV_src, orV_rcv=param.orV_rcv, + spkr_pattern=param.spkr_pattern, mic_pattern=param.mic_pattern) + + # return receiver channels (mono), number of receivers, sampling frequency and bit depth from RIRs. + return RIRs[0] \ No newline at end of file diff --git a/gpuRIR/extensions/hrtf/__init__.py b/gpuRIR/extensions/hrtf/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gpuRIR/extensions/hrtf/cipic_hrir/cipic_copyright.txt b/gpuRIR/extensions/hrtf/cipic_hrir/cipic_copyright.txt new file mode 100755 index 0000000..5eaab82 --- /dev/null +++ b/gpuRIR/extensions/hrtf/cipic_hrir/cipic_copyright.txt @@ -0,0 +1,39 @@ +Copyright +--------- + +Copyright (c) 2001 The Regents of the University of California. All Rights Reserved + + +Disclaimer +--------- + +THE REGENTS OF THE UNIVERSITY OF CALIFORNIA MAKE NO REPRESENTATION OR +WARRANTIES WITH RESPECT TO THE CONTENTS HEREOF AND SPECIFICALLY DISCLAIM ANY +IMPLIED WARRANTIES OR MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE. + +Further, the Regents of the University of California reserve the right to +revise this software and/or documentation and to make changes from time to +time in the content hereof without obligation of the Regents of the University +of California to notify any person of such revision or change. + + +Use of Materials +---------------- + +The Regents of the University of California hereby grant users permission to +reproduce and/or use materials available therein for any purpose- educational, +research or commercial. However, each reproduction of any part of the +materials must include the copyright notice, if it is present. + +In addition, as a courtesy, if these materials are used in published research, +this use should be acknowledged in the publication. If these materials are +used in the development of commercial products, the Regents of the University +of California request that written acknowledgment of such use be sent to: + + CIPIC- Center for Image Processing and Integrated Computing + University of California + 1 Shields Avenue + Davis, CA 95616-8553 + + + diff --git a/gpuRIR/extensions/hrtf/cipic_hrir/hrir_final.mat b/gpuRIR/extensions/hrtf/cipic_hrir/hrir_final.mat new file mode 100644 index 0000000..0cf5d9f Binary files /dev/null and b/gpuRIR/extensions/hrtf/cipic_hrir/hrir_final.mat differ diff --git a/gpuRIR/extensions/hrtf/hrtf_binaural_receiver.py b/gpuRIR/extensions/hrtf/hrtf_binaural_receiver.py new file mode 100644 index 0000000..88223ba --- /dev/null +++ b/gpuRIR/extensions/hrtf/hrtf_binaural_receiver.py @@ -0,0 +1,174 @@ +import numpy as np +import matplotlib.pyplot as plt + +# Helper methods for HRTF. + + +class BinauralReceiver(): + ''' Binaural receiver in a 3D space (two sound receivers in a room, representing human hearing). + ''' + + HEAD_WIDTH = 0.1449 # [m] + PINNA_OFFSET_DOWN = 0.0303 # [m] + PINNA_OFFSET_BACK = 0.0046 # [m] + + def __init__(self, head_position, head_direction, verbose=False): + ''' Instantiates a binaural receivers in a 3D space. + + Parameters + ---------- + head_position : 3D ndarray + Position of head in 3D space. + head_direction : 3D ndarray + Direction the head is pointing towards in the 3D space. + verbose : bool, optional + Terminal output for debugging or further information. + ''' + assert(head_direction[0] != 0 or head_direction[1] != 0), \ + "Ear directions are undefined for head pointing straight up or down! Try tilting a bit" + self.position = head_position + self.direction = head_direction + self.update(head_position, head_direction, verbose) + + @staticmethod + def find_spine_vector(head_direction): + ''' Finding a vector going down 90° from head direction (direction of spine). + + Parameters + ---------- + head_direction : 3D ndarray + Direction the head is pointing towards in the 3D space. + + Returns + ------- + 3D ndarray + New vector pointing down 90° from head direction (direction of spine). + ''' + h = np.copy(head_direction) + if h[2] == 0: + h[0] = 0 + h[1] = 0 + h[2] = -1 + else: + h[2] = -(h[0]**2 + h[1]**2) / h[2] + return h / np.linalg.norm(h) + + @staticmethod + def rotate_z_plane(vec, angle): + ''' Rotates the Z plane of a vector by given angle. + + Parameters + ---------- + vec : 3D ndarray + Vector to turn. + angle : float + Angle the vector is being turned. + + Returns + ------- + 3D ndarray + Rotated vector. + ''' + vec_copy = np.copy(vec) + + z_rotation = np.array([ + [np.cos(angle), -np.sin(angle), 0], + [np.sin(angle), np.cos(angle), 0], + [0, 0, 1] + ]) + vec_copy[2] = 0 + + return vec_copy @ z_rotation + + def update(self, head_position, head_direction, verbose=False): + ''' Updates head position & direction in 3D space, calculates left and right ear position. + + Parameters + ---------- + head_position : 3D ndarray + Position of head in 3D space. + head_direction : 3D ndarray + Direction the head is pointing towards in the 3D space. + verbose : bool, optional + Terminal output for debugging or further information. + ''' + self.position = head_position + self.direction = head_direction + + self.ear_direction_r = np.round( + BinauralReceiver.rotate_z_plane(self.direction, np.pi/2)) + self.ear_direction_l = -self.ear_direction_r + + ear_offset_back = BinauralReceiver.PINNA_OFFSET_BACK * \ + (self.direction / np.linalg.norm(self.direction, 2)) + + ear_offset_down = BinauralReceiver.find_spine_vector( + head_direction) * BinauralReceiver.PINNA_OFFSET_DOWN + + self.ear_position_r = (self.position + self.ear_direction_r * (BinauralReceiver.HEAD_WIDTH / 2)) #+ ear_offset_down - ear_offset_back + + self.ear_position_l = (self.position + self.ear_direction_l * (BinauralReceiver.HEAD_WIDTH / 2)) #+ ear_offset_down - ear_offset_back + + if head_direction[2] < 0: + self.ear_position_r += -ear_offset_down - ear_offset_back + self.ear_position_l += -ear_offset_down - ear_offset_back + else: + self.ear_position_r += ear_offset_down - ear_offset_back + self.ear_position_l += ear_offset_down - ear_offset_back + + if verbose: + print(f"Head position: {self.position}") + print(f"Head direction: {self.direction}") + print(f"Ear position L: {self.ear_position_l}") + print(f"Ear position R: {self.ear_position_r}") + print(f"Ear direction L: {self.ear_direction_l}") + print(f"Ear direction R: {self.ear_direction_r}") + print(f"Ear offset down: {ear_offset_down}") + print(f"Ear offset back: {ear_offset_back}") + + pass + + def visualize(self, room_sz, pos_src, orV_src): + ''' Visualizes signal source and receiver on a 3D plot. + + Parameters + ---------- + room_sz : 3D ndarray + Size of 3D room. + pos_src : 3D ndarray + Position of sorce in 3D space. + orV_src : 3D ndarray + Steering vector of source, determines which direction the source is pointing towards. + ''' + ax = plt.figure().add_subplot(projection='3d') + np.meshgrid( + np.arange(0, room_sz[0], 0.2), + np.arange(0, room_sz[1], 0.2), + np.arange(0, room_sz[2], 0.2) + ) + + # Plot origin -> head position + ax.quiver(0, 0, 0, self.position[0], self.position[1], self.position[2], + arrow_length_ratio=0, color='gray', label="Head position") + + # Plot ear directions + ax.quiver(self.ear_position_l[0], self.ear_position_l[1], self.ear_position_l[2], self.ear_direction_l[0], + self.ear_direction_l[1], self.ear_direction_l[2], length=BinauralReceiver.HEAD_WIDTH/2, color='b', label="Left ear direction") + ax.quiver(self.ear_position_r[0], self.ear_position_r[1], self.ear_position_r[2], self.ear_direction_r[0], self.ear_direction_r[1], + self.ear_direction_r[2], length=BinauralReceiver.HEAD_WIDTH/2, color='r', label="Right ear direction") + + # Plot head direction + ax.quiver(self.position[0], self.position[1], self.position[2], self.direction[0], + self.direction[1], self.direction[2], length=0.2, color='orange', label="Head direction") + + # Plot head -> signal source + ax.quiver(self.position[0], self.position[1], self.position[2], pos_src[0][0] - self.position[0], pos_src[0][1] - + self.position[1], pos_src[0][2] - self.position[2], arrow_length_ratio=0.1, color='g', label="Signal source") + # Plot head -> signal source + ax.quiver(pos_src[0][0], pos_src[0][1], pos_src[0][2], orV_src[0], orV_src[1], + orV_src[2], arrow_length_ratio=0.1, color='black', label="Source steering") + + plt.legend() + plt.show() + + pass diff --git a/gpuRIR/extensions/hrtf/hrtf_rir.py b/gpuRIR/extensions/hrtf/hrtf_rir.py new file mode 100644 index 0000000..fc289c8 --- /dev/null +++ b/gpuRIR/extensions/hrtf/hrtf_rir.py @@ -0,0 +1,87 @@ +import numpy as np +from scipy.io import loadmat +from matplotlib import pyplot as plt + +class HRTF_RIR: + ''' Head related transfer function room impulse response (HRFT RIR) interface. + + Reference: Algazi V.R., Duda R.O., Thompson D.M. and Avendano C. The CIPIC + HRTF database. In Proceedings of the 2001 IEEE Workshop on the Applications of + Signal Processing to Audio and Acoustics (Cat. No.01TH8575) (2001), pp. 99–102. + ''' + def __init__(self): + ''' Instantates interface, creates pre-defined azimuth array. + ''' + self.azimuths = np.float32([ + -80, -65, -55, -45, -40, + -35, -30, -25, -20, -15, + -10, -5, 0, 5, 10, + 15, 20, 25, 30, 35, + 40, 45, 55, 65, 80 + ]) + + self.elev_step = (360/64) + self.elevations = np.arange(-45, 230.625 + self.elev_step, self.elev_step, dtype=np.float32) + self.hrir = loadmat('gpuRIR/extensions/hrtf/cipic_hrir/hrir_final.mat') + + def azimuth_to_idx(self, azimuth): + ''' Translates azimuth to idx in order to access CIPIC database. + + Parameters + ---------- + azimuth : float + Azimuth value to translate into idx. + + Returns + ------- + float + idx + ''' + return int(np.argmin(np.abs(self.azimuths - azimuth))) + + def elevation_to_idx(self, elevation): + ''' Translates elevation to idx in order to access CIPIC database. + + Parameters + ---------- + azimuth : float + Elevation value to translate into idx. + + Returns + ------- + float + idx + ''' + return int(np.argmin(np.abs(self.elevations - elevation))) + + + def get_hrtf_rir(self, elevation, azimuth, channel, visualize=False): + """ Retrieves Head related transfer function room impulse response (HRFT RIR) + + Parameters + ---------- + elevation : float + Elevation value between source and receiver. + azimuth : float + Azimuth value between source and receiver. + channel : str + Left or right channel, represented as 'l' or 'r'. + visualize : bool + Visualizes source and receiver in a 3D space. + + Returns + ------- + 2D ndarray + HRTF RIR + """ + hrir_channel = self.hrir['hrir_'+channel][:, self.elevation_to_idx(elevation), :] + + if visualize: + x = hrir_channel[self.azimuth_to_idx(azimuth)] + plt.title("HRTF Frequency Response") + plt.rcParams.update({'font.size': 18}) + plt.plot(x) + plt.xlabel("Sample [n]") + plt.ylabel("Amplitude") + plt.show() + return hrir_channel[self.azimuth_to_idx(azimuth)] \ No newline at end of file diff --git a/gpuRIR/extensions/room_parameters.py b/gpuRIR/extensions/room_parameters.py new file mode 100644 index 0000000..b8f8659 --- /dev/null +++ b/gpuRIR/extensions/room_parameters.py @@ -0,0 +1,121 @@ +import numpy as np +import numpy.matlib as ml +import gpuRIR + + +class RoomParameters: + ''' Room parameter class for gpuRIR. + ''' + def __init__( + self, + room_sz, + pos_src, + pos_rcv, + orV_src, + orV_rcv, + spkr_pattern, + mic_pattern, + T60, + att_diff, + att_max, + fs, + bit_depth, # TODO: THhis is only to write wav files, it migh be better not including it here + abs_weights = 5 * [0.9] + [0.5], + wall_materials = None, + head_position = None, # Only for HRTF + head_direction = None, # Only for HRTF + beta = None): + ''' Instantiates a consolidated room parameter object for gpuRIR. + + Parameters + ---------- + room_sz : 3D ndarray + Size of room [m] + pos_src : 3D ndarray + Position of signal source [m] + pos_rcv : 3D ndarray + Position of singal receiver [m] + orV_src : ndarray with 2 dimensions and 3 columns or None, optional + Orientation of the sources as vectors pointing in the same direction. + None (default) is only valid for omnidirectional patterns. + orV_rcv : ndarray with 2 dimensions and 3 columns or None, optional + Orientation of the receivers as vectors pointing in the same direction. + None (default) is only valid for omnidirectional patterns. + spkr_pattern : {"omni", "homni", "card", "hypcard", "subcard", "bidir"}, optional + Polar pattern of the sources (the same for all of them). + "omni" : Omnidireccional (default). + "homni": Half omnidirectional, 1 in front of the microphone, 0 backwards. + "card": Cardioid. + "hypcard": Hypercardioid. + "subcard": Subcardioid. + "bidir": Bidirectional, a.k.a. figure 8. + mic_pattern : {"omni", "homni", "card", "hypcard", "subcard", "bidir"}, optional + Polar pattern of the receivers (the same for all of them). + "omni" : Omnidireccional (default). + "homni": Half omnidirectional, 1 in front of the microphone, 0 backwards. + "card": Cardioid. + "hypcard": Hypercardioid. + "subcard": Subcardioid. + "bidir": Bidirectional, a.k.a. figure 8. + T60 : float + Reverberation time of the room (seconds to reach 60dB attenuation). + att_diff : float + Attenuation when start using the diffuse reverberation model [dB] + att_max : float + Attenuation at the end of the simulation [dB] + fs : float + RIRs sampling frequency (in Hertz). + bit_depth : int + Bit depth of source sound data. + abs_weights : array_like with 6 elements, optional + List of six float elements, determining absorption coefficient (0..1) + wall_materials : array_like with 6 elements, optional + List of six Materials objects, representing virtual room materials. + head_position : 3D ndarray, optional + Only for HRTF. Position the head is located at. + head_direction : 3D ndarray, optional + Only for HRTF. Steering vector of the head, determining the direction it is pointing towards. + beta : array_like with 6 elements, optional + Reflection coefficients of the walls as $[beta_{x0}, beta_{x1}, beta_{y0}, beta_{y1}, beta_{z0}, beta_{z1}]$, + where $beta_{x0}$ and $beta_{x1}$ are the reflection coefficents of the walls orthogonal to the x axis at + x=0 and x=room_sz[0], respectively. + ''' + + self.room_sz = room_sz + self.pos_src = np.array(pos_src) + self.pos_rcv = np.array(pos_rcv) + self.orV_src = ml.repmat(np.array(orV_src), len(pos_src), 1) + self.orV_rcv = ml.repmat(np.array(orV_rcv), len(pos_rcv), 1) + self.spkr_pattern = spkr_pattern + self.mic_pattern = mic_pattern + self.T60 = T60 + self.fs = fs + self.bit_depth = bit_depth + self.abs_weights = abs_weights + self.wall_materials = wall_materials + + + if head_position is None: + self.head_position = pos_rcv + else: # Only for HRTF + self.head_position = head_position + + if head_direction is None: + self.head_direction = orV_rcv + else: # Only for HRTF + self.head_direction = head_direction + + if beta is None: + # Switch between self-determined wall coefficients used for frequency dependent wall absorption coefficients + self.beta = gpuRIR.beta_SabineEstimation(room_sz, T60, abs_weights=abs_weights) # Reflection coefficients + else: + self.beta = beta + + # Time to start the diffuse reverberation model [s] + self.Tdiff = gpuRIR.att2t_SabineEstimator(att_diff, T60) + + # Time to stop the simulation [s] + self.Tmax = gpuRIR.att2t_SabineEstimator(att_max, T60) + + # Number of image sources in each dimension + self.nb_img = gpuRIR.t2n(self.Tdiff, room_sz) diff --git a/gpuRIR/extensions/test/__init__.py b/gpuRIR/extensions/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gpuRIR/extensions/test/test_binaural_receiver.py b/gpuRIR/extensions/test/test_binaural_receiver.py new file mode 100644 index 0000000..71a5cef --- /dev/null +++ b/gpuRIR/extensions/test/test_binaural_receiver.py @@ -0,0 +1,124 @@ +from pytest import approx +from gpuRIR.extensions.hrtf.hrtf_binaural_receiver import BinauralReceiver + + +def test_binaural_receiver(): + # Head pointing east + head_position = [0, 0, 0] + head_direction = [1, 0, 0] + head = BinauralReceiver(head_position, head_direction) + + assert(head.ear_position_l[0] == + head_position[0]-BinauralReceiver.PINNA_OFFSET_BACK) + assert(head.ear_position_r[0] == + head_position[0]-BinauralReceiver.PINNA_OFFSET_BACK) + + assert(head.ear_position_l[1] == + head_position[1]+BinauralReceiver.HEAD_WIDTH/2) + assert(head.ear_position_r[1] == + head_position[1]-BinauralReceiver.HEAD_WIDTH/2) + + assert(head.ear_position_l[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + assert(head.ear_position_r[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + + # Head pointing west + head_position = [0, 0, 0] + head_direction = [-1, 0, 0] + head = BinauralReceiver(head_position, head_direction) + + assert(head.ear_position_l[0] == + head_position[0]+BinauralReceiver.PINNA_OFFSET_BACK) + assert(head.ear_position_r[0] == + head_position[0]+BinauralReceiver.PINNA_OFFSET_BACK) + + assert(head.ear_position_l[1] == + head_position[1]-BinauralReceiver.HEAD_WIDTH/2) + assert(head.ear_position_r[1] == + head_position[1]+BinauralReceiver.HEAD_WIDTH/2) + + assert(head.ear_position_l[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + assert(head.ear_position_r[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + + # Head pointing north + head_position = [0, 0, 0] + head_direction = [0, 1, 0] + head = BinauralReceiver(head_position, head_direction) + + assert(head.ear_position_l[0] == + head_position[0]-BinauralReceiver.HEAD_WIDTH/2) + assert(head.ear_position_r[0] == + head_position[0]+BinauralReceiver.HEAD_WIDTH/2) + + assert(head.ear_position_l[1] == + head_position[1]-BinauralReceiver.PINNA_OFFSET_BACK) + assert(head.ear_position_r[1] == + head_position[1]-BinauralReceiver.PINNA_OFFSET_BACK) + + assert(head.ear_position_l[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + assert(head.ear_position_r[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + + # Head pointing south + head_position = [0, 0, 0] + head_direction = [0, -1, 0] + head = BinauralReceiver(head_position, head_direction) + + assert(head.ear_position_l[0] == + head_position[0]+BinauralReceiver.HEAD_WIDTH/2) + assert(head.ear_position_r[0] == + head_position[0]-BinauralReceiver.HEAD_WIDTH/2) + + assert(head.ear_position_l[1] == + head_position[1]+BinauralReceiver.PINNA_OFFSET_BACK) + assert(head.ear_position_r[1] == + head_position[1]+BinauralReceiver.PINNA_OFFSET_BACK) + + assert(head.ear_position_l[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + assert(head.ear_position_r[2] == + head_position[2]-BinauralReceiver.PINNA_OFFSET_DOWN) + + # Head pointing -45° + head_position = [0, 0, 0] + head_direction = [1, 1, -1] + head = BinauralReceiver(head_position, head_direction) + + assert(head.ear_position_l[0] == + approx(-0.08747573)) + assert(head.ear_position_r[0] == + approx(0.05742427)) + + assert(head.ear_position_l[1] == + approx(0.05742427)) + assert(head.ear_position_r[1] == + approx(-0.08747573)) + + assert(head.ear_position_l[2] == + approx(-0.02208404)) + assert(head.ear_position_r[2] == + approx(-0.02208404)) + + # Head pointing +45° + head_position = [0, 0, 0] + head_direction = [1, 1, 1] + head = BinauralReceiver(head_position, head_direction) + + assert(head.ear_position_l[0] == + approx(-0.06273589)) + assert(head.ear_position_r[0] == + approx(0.08216411)) + + assert(head.ear_position_l[1] == + approx(0.08216411)) + assert(head.ear_position_r[1] == + approx(-0.06273589)) + + assert(head.ear_position_l[2] == + approx(-0.02739566)) + assert(head.ear_position_r[2] == + approx(-0.02739566)) \ No newline at end of file diff --git a/gpuRIR/extensions/test/test_hrtf.py b/gpuRIR/extensions/test/test_hrtf.py new file mode 100644 index 0000000..5183ac1 --- /dev/null +++ b/gpuRIR/extensions/test/test_hrtf.py @@ -0,0 +1,150 @@ +from pytest import approx +import numpy as np +from gpuRIR.extensions.filters.hrtf_filter import HRTF_Filter + +ANGLE_180 = approx(np.pi) + +ANGLE_90 = approx(np.pi / 2) +ANGLE_NEG_90 = approx(-(np.pi / 2)) + +ANGLE_45 = approx(np.pi / 4) +ANGLE_NEG_45 = approx(-(np.pi / 4)) + +ANGLE_135 = approx(3 * (np.pi / 4)) + +ANGLE_225 = approx(np.pi + (np.pi / 4)) + + +def test_azimuth(): + # Front + pos_src = np.array([0, 0, 0]) + pos_rcv = np.array([1, 0, 0]) + head_direction = [-1, 0, 0] + assert HRTF_Filter.calculate_azimuth(pos_src, pos_rcv, head_direction) == 0 + + # Rear + pos_src = np.array([1, 0, 0]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [1, 0, 0] + assert HRTF_Filter.calculate_azimuth(pos_src, pos_rcv, head_direction) == 0 + + # Left (-90 degrees) + pos_src = np.array([0, 0, 0]) + pos_rcv = np.array([0, 1, 0]) + head_direction = [-1, 0, 0] + assert HRTF_Filter.calculate_azimuth( + pos_src, pos_rcv, head_direction) == ANGLE_NEG_90 + + # Right (+90 degrees) + pos_src = np.array([0, 0, 0]) + pos_rcv = np.array([0, 1, 0]) + head_direction = [1, 0, 0] + assert HRTF_Filter.calculate_azimuth( + pos_src, pos_rcv, head_direction) == ANGLE_90 + + # Front Left (-45 degrees) + pos_src = np.array([-1, 1, 0]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [0, 1, 0] + assert HRTF_Filter.calculate_azimuth( + pos_src, pos_rcv, head_direction) == ANGLE_NEG_45 + + # Front Right (+45 degrees) + pos_src = np.array([1, 1, 0]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [0, 1, 0] + assert HRTF_Filter.calculate_azimuth( + pos_src, pos_rcv, head_direction) == ANGLE_45 + + # Rear Left (-45 degrees) + pos_src = np.array([-1, -1, 0]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [0, 1, 0] + assert HRTF_Filter.calculate_azimuth( + pos_src, pos_rcv, head_direction) == ANGLE_NEG_45 + + # Rear Right (+45 degrees) + pos_src = np.array([-1, 1, 0]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [0, -1, 0] + assert HRTF_Filter.calculate_azimuth( + pos_src, pos_rcv, head_direction) == ANGLE_45 + + +def test_elevation(): + # In front of head (+45 degrees) + pos_src = np.array([0, 1, 1]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [0, 1, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_45 + + # In front of head (-45 degrees) + pos_src = np.array([0, 1, -1]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [0, 1, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_NEG_45 + + # Above head (90 degrees) + pos_src = np.array([0, 0, 1]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [1, 1, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_90 + + # Below head (-90 degrees) + pos_src = np.array([0, 0, 1]) + pos_rcv = np.array([0, 0, 2]) + head_direction = [1, 1, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_NEG_90 + + # In front of head (0 degrees) (same elevation) + pos_src = np.array([0, 0, 0]) + pos_rcv = np.array([1, 0, 0]) + head_direction = [-1, 0, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == 0 + + # Behind head (180 degrees) (same elevation) + pos_src = np.array([-1, 0, 0]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [1, 0, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_180 + + # Behind head (224 degrees) + pos_src = np.array([-1, 0, -1]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [1, 0, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_225 + + # Behind head (135 degrees) + pos_src = np.array([-1, 0, 1]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [1, 0, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_135 + + # Behind head (153 degrees) + pos_src = np.array([-1, 0, 0.5]) + pos_rcv = np.array([0, 0, 0]) + head_direction = [1, 0, 0] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == approx(np.pi * .8524163) + + # Above head, head tilted 45 degrees upward + pos_src = np.array([1, 1, 2]) + pos_rcv = np.array([1, 1, 1]) + head_direction = [0, 1, 1] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_45 + + # Above head, head tilted 45 degrees down + pos_src = np.array([1, 1, 0]) + pos_rcv = np.array([1, 1, 1]) + head_direction = [0, 1, -1] + assert HRTF_Filter.calculate_elevation( + pos_src, pos_rcv, head_direction) == ANGLE_NEG_45 diff --git a/gpuRIR/extensions/wall_absorption/__init__.py b/gpuRIR/extensions/wall_absorption/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gpuRIR/extensions/wall_absorption/freq_dep_abs_coeff.py b/gpuRIR/extensions/wall_absorption/freq_dep_abs_coeff.py new file mode 100644 index 0000000..de19c1d --- /dev/null +++ b/gpuRIR/extensions/wall_absorption/freq_dep_abs_coeff.py @@ -0,0 +1,183 @@ +import gpuRIR.extensions.generate_RIR as gRIR + +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +from gpuRIR.extensions.filters.butterworth import Butterworth + + +def interpolate_pair(abs_coeff, visualize=False): + """Interpolates frequency response array. Returns a function. + + Parameters + ---------- + abs_coeff : 2D ndarray + Per-frequency-range absorption coefficient of a virtual wall material, as (pre-) defined in extensions/wall_absorption/materials.py. + visualize : bool + Plots the interpolated frequency responses on a 2D graph. + + Returns + ------- + functions + Interpolation function of a single virtual room material as (pre-) defined in extensions/wall_absorption/materials.py. + """ + # y: absorption coefficient + # x: frequency [Hz] + y = abs_coeff[:, 1] + x = abs_coeff[:, 0] + f = interp1d(x, y, bounds_error=False, fill_value=(y[0], y[-1])) + x_interpolated = np.arange(1, 20000) # TODO: Should 20000 be adjustable? + y_interpolated = f(x_interpolated) + if visualize: + plt.title( + 'Frequency Dependent Absorption Coefficients: Interpolated Frequency Responses') + plt.plot(x, y, 'o') + plt.plot(x_interpolated, y_interpolated, "-") + return f + + +def generate_RIR_freq_dep_walls(params, band_width=100, factor=1.5, order=4, LR=False, use_bandpass_on_borders=False, visualize=False, verbose=False): + """ Generates a custom room impulse response (RIR) with virtual room materials applied using a bandpassing method. + + Parameters: + ----------- + params : RoomParameters + gpuRIR parameters + band_width : int, optional + Initial width of frequency band. Lower means higher quality but less performant. (recommended: 10) + factor : float, optional + Multiplication factor of frequency band. Lower means higher quality but less performant. (recommended: 1.1) + order : int, optional + Butterworth filter order. + LR : bool, optional + Enables Linkwitz-Riley filtering. LR filter order will get converted automatically. + use_bandpass_on_borders : bool, optional + Uses bandpass instead of high/lowpass filters on lowest and highest band. + plot : bool, optional + Plots the interpolated material frequency response curve. + verbose : bool, optional + Prints current band parameters. + + Returns + ------- + 2D ndarray + Processed Room impulse response array. + """ + + assert(band_width > 1), "Band width must be greater than 1!" + assert(factor > 1), "Factor must be greater than 1!" + assert(order >= 4), "Order must be greater than 4!" + + if LR: + order = order / 2 + + min_frequency = 1 + max_frequency = 20000 + + # Structure: min / mean / max + bands = [] + + # Find out minimum and maximum frequency value of all materials + material_frequencies = np.vstack(params.wall_materials)[:, 0] + min_mat_frequency = np.min(material_frequencies) + max_mat_frequency = np.max(material_frequencies) + + # Outside of the material bounds, absorption values are constant. (Not wasting bands) + bands.append( + [min_frequency, (min_frequency + min_mat_frequency) / 2, min_mat_frequency]) + + if verbose: + print( + f"Min:{min_frequency}\tMean:{(min_frequency + min_mat_frequency) / 2}\tMax:{min_mat_frequency}") + + reached_max = False + + current_min = min_mat_frequency + current_max = current_min + band_width + current_mean = (current_min + current_max) / 2 + + while not reached_max: + if current_max > max_mat_frequency: + reached_max = True + + bands.append([current_min, current_mean, current_max]) + + if verbose: + print( + f"Min:{current_min}\tMean:{current_mean}\tMax:{current_max}\tBand width:{band_width}") + + band_width *= factor + + current_min = current_max + current_max = current_min + band_width + current_mean = (current_min + current_max) / 2 + + bands.append( + [current_min, (current_min + max_frequency) / 2, max_frequency]) + if verbose: + print( + f"Min:{current_min}\tMean:{(current_min + max_frequency) / 2}\tMax:{max_frequency}") + + # We create 6 interpolating functions for each material: + wall_mat_interp = [interpolate_pair(mat, visualize) + for mat in params.wall_materials] + if visualize: + plt.xlim(right=20000) + plt.show() + + receiver_channels = np.zeros((len(params.pos_rcv), 1)) + + if visualize: + plt.xscale('log') + plt.title( + 'Frequency Dependent Absorption Coefficients: Butterworth filter frequency response') + plt.xlabel('Frequency [Hz]') + plt.ylabel('Amplitude [dB]') + plt.ylim(bottom=-40) + plt.margins(0, 0.1) + plt.grid(which='both', axis='both') + + for i in range(len(bands)): + band = bands[i] + abs_coeffs = np.zeros(len(wall_mat_interp)) + for j in range(len(wall_mat_interp)): + abs_coeffs[j] = wall_mat_interp[j](band[1]) + # Generate RIR + params.beta = 6 * [1.] - abs_coeffs + RIR = gRIR.generate_RIR(params) # TODO: Computing all the octaves in parallel in the GPU would be far more + # efficient, but this mean a lot of modifications on the CUDA code. + + # Apply band/lowpassing and re-compiling sound data + for rcv in range(len(params.pos_rcv)): + # Lowpass lowest frequency band + if i == 0 and not use_bandpass_on_borders: + processed = Butterworth.apply_pass_filter( + RIR[rcv], band[2], params.fs, 'lowpass', order * 4, visualize + ) + if verbose: + print(f"Lowpass frequency: {band[2]} Order: {order * 4}") + + # Highpass highest frequency band + elif i == (len(bands) - 1) and not use_bandpass_on_borders: + processed = Butterworth.apply_pass_filter( + RIR[rcv], band[0], params.fs, 'highpass', order * 4, visualize + ) + if verbose: + print(f"Highpass frequency: {band[0]} Order: {order * 4}") + + # All other bands are bandpassed + else: + processed = Butterworth.apply_bandpass_filter( + RIR[rcv], band[0], band[2], params.fs, order, LR, visualize) + if verbose: + print(f"Band {i} Bandpass frequency: {band[0]} Order: {order}") + + # Re-compiling sound data + receiver_channels.resize(len(params.pos_rcv), len(processed)) + receiver_channels[rcv] += processed + + if visualize: + plt.show() + + # Sum up all bandpassed RIR's per receiver + return receiver_channels diff --git a/gpuRIR/extensions/wall_absorption/materials.py b/gpuRIR/extensions/wall_absorption/materials.py new file mode 100644 index 0000000..9563582 --- /dev/null +++ b/gpuRIR/extensions/wall_absorption/materials.py @@ -0,0 +1,121 @@ +import numpy as np + +# Source of all material absorption coefficients: +# Wolfgang M. Willems, Kai Schild, Diana Stricker: Formeln und Tabellen Bauphysik. Fourth edition, p. 441. + +class Materials: + """ Virtual room material characteristics, represented as 2D array, built on frequency-attenuation value pairs. + """ + + + # Building materials + concrete = np.array([ + [125, 0.02], + [250, 0.02], + [500, 0.03], + [1000, 0.04], + [2000, 0.05], + [4000, 0.05] + ]) + + carpet_10mm = np.array([ + [125, 0.04], + [250, 0.07], + [500, 0.12], + [1000, 0.3], + [2000, 0.5], + [4000, 0.8] + ]) + + # Wood + parquet_glued = np.array([ + [125, 0.04], + [250, 0.04], + [500, 0.05], + [1000, 0.06], + [2000, 0.06], + [4000, 0.06] + ]) + + # Glass + window = np.array([ + [125, 0.28], + [250, 0.2], + [500, 0.1], + [1000, 0.06], + [2000, 0.03], + [4000, 0.02] + ]) + + mirror_on_wall = np.array([ + [125, 0.12], + [250, 0.1], + [500, 0.05], + [1000, 0.04], + [2000, 0.02], + [4000, 0.02] + ]) + + # Special + wallpaper_on_lime_cement_plaster = np.array([ + [125, 0.02], + [250, 0.03], + [500, 0.04], + [1000, 0.05], + [2000, 0.07], + [4000, 0.08] + ]) + + bookshelf = np.array([ + [125, 0.3], + [250, 0.4], + [500, 0.4], + [1000, 0.3], + [2000, 0.3], + [4000, 0.2] + ]) + + cinema_screen = np.array([ + [125, 0.1], + [250, 0.1], + [500, 0.2], + [1000, 0.3], + [2000, 0.5], + [4000, 0.6] + ]) + + total_absorption = np.array([ + [125, 0.9999999], + [250, 0.9999999], + [500, 0.9999999], + [1000, 0.9999999], + [2000, 0.9999999], + [4000, 0.9999999] + ]) + + total_reflection = np.array([ + [125, 0.000001], + [250, 0.000001], + [500, 0.000001], + [1000, 0.000001], + [2000, 0.000001], + [4000, 0.000001] + ]) + + absorb_low_freqs = np.array([ + [125, 0.9999999], + [250, 0.9999999], + [500, 0.000001], + [1000, 0.000001], + [2000, 0.000001], + [4000, 0.000001] + ]) + + absorb_high_freqs = np.array([ + [125, 0.000001], + [250, 0.000001], + [500, 0.000001], + [1000, 0.000001], + [2000, 0.000001], + [4000, 0.999999] + ])