Implement source directivity and polar pattern (#28)

* Add our custom helper scripts.

Co-Authored-By: fuerbringer <[email protected]>

* Create air absorption script.

Co-Authored-By: fuerbringer <[email protected]>

* Create script for generation IR file with air
absorption. Not yet working

Co-Authored-By: fuerbringer <[email protected]>

* Continue on the bandpass route. Doesn't work too
well right now

Co-Authored-By: fuerbringer <[email protected]>

* Fixing a small mistake

Co-Authored-By: fuerbringer <[email protected]>

* Successfully implemented air absorption through
splitting frequency ranges into
different bands.

Co-Authored-By: fuerbringer <[email protected]>

* Change default filter order settings.

Co-Authored-By: fuerbringer <[email protected]>

* Expanded scripts for more versatile utilization
(importing into other scripts, terminal arguments)

Co-Authored-By: fuerbringer <[email protected]>

* Implement STFT air absorption

Co-Authored-By: fuerbringer <[email protected]>

* Not working yet.

Co-Authored-By: fuerbringer <[email protected]>

* Not working yet

Co-Authored-By: fuerbringer <[email protected]>

* Implemented strategy pattern for filters.

Co-Authored-By: fuerbringer <[email protected]>

* Slight improvements, everything works

Co-Authored-By: fuerbringer <[email protected]>

* Calculate difference

* Finalize work on air absorption

Co-Authored-By: fuerbringer <[email protected]>

* Adjust standard parameters to ease repeated
testing

Co-Authored-By: fuerbringer <[email protected]>

* Update gitignore.

Co-Authored-By: fuerbringer <[email protected]>

* Attempt to multithread bandpass filter.

* Improve bandpass performance drastically.

Co-Authored-By: fuerbringer <[email protected]>
Co-Authored-By: TheBlueFireFox <[email protected]>

* Further optimizations

Co-Authored-By: fuerbsev <[email protected]>

* Make preparations. Not working yet.

Co-Authored-By: fuerbringer <[email protected]>

* Small cleanup, add possibility for forwards
and backwards filtering

Co-Authored-By: fuerbringer <[email protected]>

* Conducted experiments with linear filters
Ready to implement into code

Co-Authored-By: fuerbringer <[email protected]>

* Implement linear filter using strategy pattern.

* Cleaned up code, improved project structure

* Further small improvements and cleanup

* Implement characteristic filtering via STFT.

Add SM57 and iPhone X frequency responce model.

Co-Authored-By: fuerbringer <[email protected]>

* Change file name, small cleanup

Co-Authored-By: fuerbringer <[email protected]>

* Change file name

Co-Authored-By: fuerbringer <[email protected]>

* Implement source characteristic.
We brought greatness upon the world.

Co-Authored-By: fuerbringer <[email protected]>

* Finished frequency dependent sender and receiver
characteristics.
Next up: Refactoring and consolidating generate_IR script to simplify
the application of filters/frequency responses.

Co-Authored-By: fuerbringer <[email protected]>

* Improve spectrogram.
Further discussion needed.

Co-Authored-By: fuerbringer <[email protected]>

* Add colorbar for dB Legend to spectro.

* Consolidate generate_IR script.
Implement linear filter.
Big cleanup.

Co-Authored-By: fuerbringer <[email protected]>

* Implement speaker directivity.
Just hardcoded values yet, yet to analyze results and expand parameters.

Co-Authored-By: fuerbringer <[email protected]>

* Did the heckin' removerino ;)

Co-Authored-By: fuerbringer <[email protected]>

* Fix db colorbar scale and add db cap.

* Fixing logical errors in source directivity.

Co-Authored-By: fuerbringer <[email protected]>

* Improved spectrogram, working perfectly now

Co-Authored-By: fuerbringer <[email protected]>

* Make font size bigger, relabel Hz to KHz

Co-Authored-By: fuerbringer <[email protected]>

* Squashed a bug in source directivity.

Co-Authored-By: fuerbringer <[email protected]>

* Added parameters, updated python binding
Minor refactoring

Co-Authored-By: fuerbringer <[email protected]>

* Create automated directivity test script.

Co-Authored-By: fuerbringer <[email protected]>

* Fixed comments.

Co-Authored-By: fuerbringer <[email protected]>

* Create polar pattern plot script.

Co-Authored-By: fuerbringer <[email protected]>

* Small code cleanup

Co-Authored-By: fuerbringer <[email protected]>

* Set mic/speaker interpolation to interp1d

Co-Authored-By: fuerbringer <[email protected]>

* Corrected errors in polar plot

Co-Authored-By: fuerbringer <[email protected]>

* New polar pattern using plotly library

Co-Authored-By: fuerbringer <[email protected]>

* Consolidated everythign around directivity testing
and plotting in one script.

Co-Authored-By: fuerbringer <[email protected]>

* Delete old polar plot scripts.

Co-Authored-By: fuerbringer <[email protected]>

* Add material list

Co-Authored-By: fuerbringer <[email protected]>

* Work in progress

Co-Authored-By: fuerbringer <[email protected]>

* Make frequency dependant wall absorption
coefficents work

Co-Authored-By: fuerbringer <[email protected]>

* Consolidate frequency dependent
absoption coefficient except multiple receiver channels

Co-Authored-By: fuerbringer <[email protected]>

* Make multiple receivers work.

Co-Authored-By: fuerbringer <[email protected]>

* Consolidate further, frequency dependent wall
materials broken at the moment

Co-Authored-By: fuerbringer <[email protected]>

* Finish implementation, everything works

Co-Authored-By: fuerbringer <[email protected]>

* Corrected a wrong value in concrete_coarse
and small formatting

Co-Authored-By: fuerbringer <[email protected]>

* Change interpolation
Instead of extrapolating values, keep first/last value in interp1d

* Expand comment

Co-Authored-By: fuerbringer <[email protected]>

* Concentrated band division on actual frequency
response of materials

Co-Authored-By: fuerbringer <[email protected]>

* Implement logarithmic distribution of bands.

Co-Authored-By: fuerbringer <[email protected]>

* Fixed verbose printouts.

Co-Authored-By: fuerbringer <[email protected]>

* Change materials to accomodate a more reputable source.

Co-Authored-By: fuerbringer <[email protected]>

* Update README.md

Add spkr_pattern and orV_src.

* Create README.md

Not yet done!

* Made a sketch for LR filtering.
Still having the low frequency band anomaly

Co-Authored-By: fuerbringer <[email protected]>

* Improve bandpassing
by adding high and low pass on frequency boundaries

Co-Authored-By: fuerbsev <[email protected]>

* Improve bandpassing in air absorption
With low- and highpass filters

Co-Authored-By: fuerbringer <[email protected]>

* Cleanup

Co-Authored-By: fuerbringer <[email protected]>

* Implement HRTF.
Not tested yet.

Co-Authored-By: fuerbringer <[email protected]>

* Made HRTF work.

Co-Authored-By: fuerbringer <[email protected]>

* Add 3D plot for visualizing hrtf

Co-Authored-By: fuerbringer <[email protected]>

* Fixed an error concering rotation of head

Co-Authored-By: fuerbringer <[email protected]>

* Fix some issues with hrft, verify results
Expand spectrogram method
Expand adaptive gain to support stereo
Various small improvements and bug fixes

Co-Authored-By: fuerbringer <[email protected]>

* Corrected ear offset

Co-Authored-By: fuerbringer <[email protected]>

* Update README.md

* Update README.md

* Update README.md

* Cleanup

Co-Authored-By: fuerbringer <[email protected]>

* Fixed Azimuth

Co-Authored-By: fuerbringer <[email protected]>

* Fix elevation, not done yet

Fixed bug in azimuth

Co-Authored-By: fuerbringer <[email protected]>

* Write HRTF calculation tests to verify results.
Elevation still needs some work.

Co-Authored-By: fuerbringer <[email protected]>

* Begin refactoring

Co-Authored-By: fuerbringer <[email protected]>

* Massive refactoring.
Consolidated all scripts and improved UX massively.
All functionality is now in gpuRIR module.

Co-Authored-By: fuerbringer <[email protected]>

* Changed butterworth filters to sos.

Co-Authored-By: fuerbringer <[email protected]>

* Make butterworth helper class
implement sos in freq_dep_abs_coeff

Co-Authored-By: fuerbringer <[email protected]>

* Fix elevation, complete test cases

Co-Authored-By: fuerbringer <[email protected]>

* Write tests for binaural receiver
Fix pinna offset errors

Co-Authored-By: fuerbringer <[email protected]>

* Comments and cleanup

Co-Authored-By: fuerbringer <[email protected]>

* Refactoring polar_plot.py
Extending room parameters with beta

Co-Authored-By: fuerbringer <[email protected]>

* Update readme

Co-Authored-By: fuerbringer <[email protected]>

* Additions to documentation / comments

Co-Authored-By: fuerbringer <[email protected]>

* More comments, renamed variables for consistency

Co-Authored-By: fuerbringer <[email protected]>

* Fixed bug with spectrogram generation
Further refactoring

Co-Authored-By: fuerbringer <[email protected]>

* Fixed bug with spectrogram generation
Further refactoring
More comments

Co-Authored-By: fuerbringer <[email protected]>

* Added further comments

Co-Authored-By: fuerbringer <[email protected]>

* Added bandpass/hi-low-pass switch

Co-Authored-By: fuerbringer <[email protected]>

* Update readme.

* Small corrections

* Reverted stereo_filters.py to standard settings.

* Prepare source directivity for upstream pull request.

Co-Authored-By: fuerbringer <[email protected]>

* Remove gpuRIR.extensions imports not yet in this branch.

Co-Authored-By: fuerbringer <[email protected]>

Co-authored-by: fuerbringer <[email protected]>
Co-authored-by: corrooli <[email protected]>
Co-authored-by: TheBlueFireFox <[email protected]>
Co-authored-by: fuerbsev <[email protected]>

.gitignore

@@ -37,6 +37,17 @@
 ehthumbs.db
 Thumbs.db
 
+# Outputs and graphs #
+######################
+*.png
+*.wav
+
+# Generated by VSCode #
+#######################
+config.py
+objectdb
+*.pyc
+
 # Folders #
 ###########
 build/

README.md

@@ -61,16 +61,20 @@ Room Impulse Responses (RIRs) simulation using the Image Source Method (ISM). Fo
         RIRs sampling frequency (in Hertz).
 * **Tdiff** : *float, optional*
         Time (in seconds) when the ISM is replaced by a diffuse reverberation model. Default is Tmax (full ISM simulation).
+* **spkr_pattern** : *{"omni", "homni", "card", "hypcard", "subcard", "bidir"}, optional.*
+        Polar pattern of the sources (the same for all of them).
 * **mic_pattern** : *{"omni", "homni", "card", "hypcard", "subcard", "bidir"}, optional.*
-        Polar pattern of the receivers (the same for all of them).
+	Polar pattern of the receivers (the same for all of them).
 	* *"omni"* : Omnidireccional (default).
-	* *"homni"*: Half omnidireccional, 1 in front of the microphone, 0 backwards.
+	* *"homni"*: Half omnidirectional, 1 in front of the microphone, 0 backwards.
 	* *"card"*: Cardioid.
 	* *"hypcard"*: Hypercardioid.
 	* *"subcard"*: Subcardioid.
 	* *"bidir"*: Bidirectional, a.k.a. figure 8.
+* **orV_src** : *ndarray with 2 dimensions and 3 columns or None, optional.*
+        Orientation of the sources as vectors pointing in the same direction. Applies to each source. 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 omnidireccional patterns.
+        Orientation of the receivers as vectors pointing in the same direction. Applies to each receiver. None (default) is only valid for omnidirectional patterns.
 * **c** : *float, optional.*
         Speed of sound (in m/s). The default is 343.0.
 

examples/polar_plots.py

@@ -0,0 +1,156 @@
+import numpy as np
+import numpy.matlib
+import matplotlib.pyplot as plt
+from math import ceil
+import time
+import gpuRIR
+from scipy.io import wavfile
+from scipy import signal
+import plotly.graph_objects as go
+from plotly.subplots import make_subplots
+
+""" 
+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!
+"""
+
+# Feel free to change these parameters
+# Resolution of polar plot (amount to divide 360 degrees into)
+PARTITIONS = 360
+
+# gpuRIR parameters
+gpuRIR.activateMixedPrecision(False)
+gpuRIR.activateLUT(False)
+
+# Don't change these things
+POLAR_PATTERNS = np.array(
+    ["omni", "homni", "card", "hypcard", "subcard", "bidir"])
+MAX_VALUES = np.zeros((POLAR_PATTERNS.shape[0], PARTITIONS))
+
+
+def normalize_amps(amps):
+    '''
+    Normalizing amplitude. Change all amplitude such that loudest sample has the amplitude of 1.
+    :param amps Array of samples.
+    '''
+    amps_normalized = np.copy(amps)
+    max_value = np.max(amps_normalized)
+    for i in range(0, len(amps_normalized)):
+        amps_normalized[i] /= max_value
+    return np.abs(amps_normalized)
+
+
+COLORS = ['mediumseagreen', 'darkorange',
+          'mediumpurple', 'magenta', 'limegreen', 'royalblue']
+
+
+def create_polar_plot(fig, i, amps, name):
+    '''
+    Creates single polar plot.
+
+    :param fig Plot figure.
+    :param i Index of plot
+    :param name Name of polar pattern
+    '''
+    fig.add_trace(go.Scatterpolar(
+        r=normalize_amps(amps),
+        theta=np.linspace(0, 360, len(amps)),
+        mode='lines',
+        name=name,
+        line_color=COLORS[i],
+        line_width=3,
+        subplot=f"polar{i+1}"
+    ), 1 if i <= 2 else 2, (i % 3) + 1)
+
+
+def create_polar_plots(title=""):
+    '''
+    Creates compilation of polar plots using plotly.
+
+    :param title Title of plot.
+    '''
+    fig_polar = make_subplots(rows=2, cols=3, specs=[[{'type': 'polar'}]*3]*2)
+    for i in range(6):
+        create_polar_plot(fig_polar, i, MAX_VALUES[i], POLAR_PATTERNS[i])
+    fig_polar.update_layout(
+        font_size=22,
+        showlegend=True,
+
+        polar1=dict(
+            radialaxis=dict(type="log", tickangle=45),
+            radialaxis_range=[-1, 0.1],
+            radialaxis_tickfont_size=18
+        ),
+        polar2=dict(
+            radialaxis=dict(type="log", tickangle=45),
+            radialaxis_range=[-1, 0.1],
+            radialaxis_tickfont_size=18
+        ),
+        polar3=dict(
+            radialaxis=dict(type="log", tickangle=45),
+            radialaxis_range=[-1, 0.1],
+            radialaxis_tickfont_size=18
+        ),
+        polar4=dict(
+            radialaxis=dict(type="log", tickangle=45),
+            radialaxis_range=[-1, 0.1],
+            radialaxis_tickfont_size=18
+        ),
+        polar5=dict(
+            radialaxis=dict(type="log", tickangle=45),
+            radialaxis_range=[-1, 0.1],
+            radialaxis_tickfont_size=18
+        ),
+        polar6=dict(
+            radialaxis=dict(type="log", tickangle=45),
+            radialaxis_range=[-1, 0.1],
+            radialaxis_tickfont_size=18
+        ),
+    )
+    fig_polar.show()
+
+
+if __name__ == "__main__":
+    fig_wf = plt.figure(1)
+    fig_sp = plt.figure(2)
+
+    for p in range(len(POLAR_PATTERNS)):
+        for i in range(0, PARTITIONS):
+            degree = i * (360 / PARTITIONS)
+            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
+
+            # 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)
+
+            # Stack array vertically
+            impulseResponseArray = np.vstack(receiver_channels[0])
+
+            # Extract biggest peak for polar pattern plotting
+            MAX_VALUES[p][i] = np.max(impulseResponseArray)
+            negative_peak = np.abs(np.min(impulseResponseArray))
+            if MAX_VALUES[p][i] < negative_peak:
+                MAX_VALUES[p][i] = negative_peak
+
+    create_polar_plots()

examples/simulate_trajectory.py

@@ -30,7 +30,6 @@
 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
 RIRs = gpuRIR.simulateRIR(room_sz, beta, pos_traj, pos_rcv, nb_img, Tmax, fs, Tdiff=Tdiff, orV_rcv=orV_rcv, mic_pattern=mic_pattern)
-
 filtered_signal = gpuRIR.simulateTrajectory(source_signal, RIRs)
 wavfile.write('filtered_signal.wav', fs, filtered_signal)
 plt.plot(filtered_signal)

gpuRIR/__init__.py

@@ -8,9 +8,9 @@
 
 from gpuRIR_bind import gpuRIR_bind
 
-__all__ = ["mic_patterns", "beta_SabineEstimation", "att2t_SabineEstimator", "t2n", "simulateRIR", "simulateTrajectory", "activate_mixed_precision", "activate_lut"]
+__all__ = ["polar_patterns", "beta_SabineEstimation", "att2t_SabineEstimator", "t2n", "simulateRIR", "simulateTrajectory", "activate_mixed_precision", "activate_lut"]
 
-mic_patterns =	{
+polar_patterns =	{
   "omni": 0,
   "homni": 1,
   "card": 2,
@@ -92,7 +92,7 @@ def t2n(T, rooms_sz, c=343.0):
 	nb_img = 2 * T / (np.array(rooms_sz) / c)
 	return [ int(n) for n in np.ceil(nb_img) ]
 
-def simulateRIR(room_sz, beta, pos_src, pos_rcv, nb_img, Tmax, fs, Tdiff=None, mic_pattern="omni", orV_rcv=None, c=343.0):
+def simulateRIR(room_sz, beta, pos_src, pos_rcv, nb_img, Tmax, fs, Tdiff=None, spkr_pattern="omni", mic_pattern="omni", orV_src=None, orV_rcv=None, c=343.0):
 	''' Room Impulse Responses (RIRs) simulation using the Image Source Method (ISM).
 
 	Parameters
@@ -114,17 +114,28 @@ def simulateRIR(room_sz, beta, pos_src, pos_rcv, nb_img, Tmax, fs, Tdiff=None, m
 	Tdiff : float, optional
 		Time (in seconds) when the ISM is replaced by a diffuse reverberation model.
 		Default is Tmax (full ISM simulation).
+	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 omnidireccional, 1 in front of the microphone, 0 backwards.
+			"homni": Half omnidirectional, 1 in front of the microphone, 0 backwards.
 			"card": Cardioid.
 			"hypcard": Hypercardioid.
 			"subcard": Subcardioid.
 			"bidir": Bidirectional, a.k.a. figure 8.
+	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 omnidireccional patterns.
+		None (default) is only valid for omnidirectional patterns.
 	c : float, optional
 		Speed of sound [m/s] (the default is 343.0).
 
@@ -139,21 +150,29 @@ def simulateRIR(room_sz, beta, pos_src, pos_rcv, nb_img, Tmax, fs, Tdiff=None, m
 	the GPU memory and crash the kernel.
 
 	'''
+
 	assert not ((pos_src >= room_sz).any() or (pos_src <= 0).any()), "The sources must be inside the room"
 	assert not ((pos_rcv >= room_sz).any() or (pos_rcv <= 0).any()), "The receivers must be inside the room"
 	assert Tdiff is None or Tdiff <= Tmax, "Tmax must be equal or greater than Tdiff"
-	assert mic_pattern in mic_patterns, "mic_pattern must be omni, homni, card, hypcard, subcard or bidir"
+	assert mic_pattern in polar_patterns, "mic_pattern must be omni, homni, card, hypcard, subcard or bidir"
 	assert mic_pattern == "omni" or orV_rcv is not None, "the mics are not omni but their orientation is undefined"
+	assert spkr_pattern in spkr_pattern, "spkr_pattern must be omni, homni, card, hypcard, subcard or bidir"
+	assert spkr_pattern == "omni" or orV_src is not None, "the sources are not omni but their orientation is undefined"
+
 
 	pos_src = pos_src.astype('float32', order='C', copy=False)
 	pos_rcv = pos_rcv.astype('float32', order='C', copy=False)
 
 	if Tdiff is None: Tdiff = Tmax
 	if mic_pattern is None: mic_pattern = "omni"
+	if spkr_pattern is None: spkr_pattern = "omni"
 	if orV_rcv is None: orV_rcv = np.zeros_like(pos_rcv)
 	else: orV_rcv = orV_rcv.astype('float32', order='C', copy=False)
+	if orV_src is None: orV_src = np.zeros_like(pos_src)
+	else: orV_src = orV_src.astype('float32', order='C', copy=False)
+
 
-	return gpuRIR_bind_simulator.simulateRIR_bind(room_sz, beta, pos_src, pos_rcv, orV_rcv, mic_patterns[mic_pattern], nb_img, Tdiff, Tmax, fs, c)
+	return gpuRIR_bind_simulator.simulateRIR_bind(room_sz, beta, pos_src, pos_rcv, orV_src, orV_rcv, polar_patterns[spkr_pattern], polar_patterns[mic_pattern], nb_img, Tdiff, Tmax, fs, c)
 
 def simulateTrajectory(source_signal, RIRs, timestamps=None, fs=None):
 	''' Filter an audio signal by the RIRs of a motion trajectory recorded with a microphone array.