From 93b52e5367b4578e977388a83e144d565e65c1fd Mon Sep 17 00:00:00 2001 From: Oliver Corrodi <61276147+corrooli@users.noreply.github.com> Date: Thu, 16 Dec 2021 14:08:18 +0100 Subject: [PATCH] Implement source directivity and polar pattern (#28) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add our custom helper scripts. Co-Authored-By: fuerbringer * Create air absorption script. Co-Authored-By: fuerbringer * Create script for generation IR file with air absorption. Not yet working Co-Authored-By: fuerbringer * Continue on the bandpass route. Doesn't work too well right now Co-Authored-By: fuerbringer * Fixing a small mistake Co-Authored-By: fuerbringer * Successfully implemented air absorption through splitting frequency ranges into different bands. Co-Authored-By: fuerbringer * Change default filter order settings. Co-Authored-By: fuerbringer * Expanded scripts for more versatile utilization (importing into other scripts, terminal arguments) Co-Authored-By: fuerbringer * Implement STFT air absorption Co-Authored-By: fuerbringer * Not working yet. Co-Authored-By: fuerbringer * Not working yet Co-Authored-By: fuerbringer * Implemented strategy pattern for filters. Co-Authored-By: fuerbringer * Slight improvements, everything works Co-Authored-By: fuerbringer * Calculate difference * Finalize work on air absorption Co-Authored-By: fuerbringer * Adjust standard parameters to ease repeated testing Co-Authored-By: fuerbringer * Update gitignore. Co-Authored-By: fuerbringer * Attempt to multithread bandpass filter. * Improve bandpass performance drastically. Co-Authored-By: fuerbringer Co-Authored-By: TheBlueFireFox * Further optimizations Co-Authored-By: fuerbsev * Make preparations. Not working yet. Co-Authored-By: fuerbringer * Small cleanup, add possibility for forwards and backwards filtering Co-Authored-By: fuerbringer * Conducted experiments with linear filters Ready to implement into code Co-Authored-By: fuerbringer * 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 * Change file name, small cleanup Co-Authored-By: fuerbringer * Change file name Co-Authored-By: fuerbringer * Implement source characteristic. We brought greatness upon the world. Co-Authored-By: fuerbringer * 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 * Improve spectrogram. Further discussion needed. Co-Authored-By: fuerbringer * Add colorbar for dB Legend to spectro. * Consolidate generate_IR script. Implement linear filter. Big cleanup. Co-Authored-By: fuerbringer * Implement speaker directivity. Just hardcoded values yet, yet to analyze results and expand parameters. Co-Authored-By: fuerbringer * Did the heckin' removerino ;) Co-Authored-By: fuerbringer * Fix db colorbar scale and add db cap. * Fixing logical errors in source directivity. Co-Authored-By: fuerbringer * Improved spectrogram, working perfectly now Co-Authored-By: fuerbringer * Make font size bigger, relabel Hz to KHz Co-Authored-By: fuerbringer * Squashed a bug in source directivity. Co-Authored-By: fuerbringer * Added parameters, updated python binding Minor refactoring Co-Authored-By: fuerbringer * Create automated directivity test script. Co-Authored-By: fuerbringer * Fixed comments. Co-Authored-By: fuerbringer * Create polar pattern plot script. Co-Authored-By: fuerbringer * Small code cleanup Co-Authored-By: fuerbringer * Set mic/speaker interpolation to interp1d Co-Authored-By: fuerbringer * Corrected errors in polar plot Co-Authored-By: fuerbringer * New polar pattern using plotly library Co-Authored-By: fuerbringer * Consolidated everythign around directivity testing and plotting in one script. Co-Authored-By: fuerbringer * Delete old polar plot scripts. Co-Authored-By: fuerbringer * Add material list Co-Authored-By: fuerbringer * Work in progress Co-Authored-By: fuerbringer * Make frequency dependant wall absorption coefficents work Co-Authored-By: fuerbringer * Consolidate frequency dependent absoption coefficient except multiple receiver channels Co-Authored-By: fuerbringer * Make multiple receivers work. Co-Authored-By: fuerbringer * Consolidate further, frequency dependent wall materials broken at the moment Co-Authored-By: fuerbringer * Finish implementation, everything works Co-Authored-By: fuerbringer * Corrected a wrong value in concrete_coarse and small formatting Co-Authored-By: fuerbringer * Change interpolation Instead of extrapolating values, keep first/last value in interp1d * Expand comment Co-Authored-By: fuerbringer * Concentrated band division on actual frequency response of materials Co-Authored-By: fuerbringer * Implement logarithmic distribution of bands. Co-Authored-By: fuerbringer * Fixed verbose printouts. Co-Authored-By: fuerbringer * Change materials to accomodate a more reputable source. Co-Authored-By: fuerbringer * 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 * Improve bandpassing by adding high and low pass on frequency boundaries Co-Authored-By: fuerbsev * Improve bandpassing in air absorption With low- and highpass filters Co-Authored-By: fuerbringer * Cleanup Co-Authored-By: fuerbringer * Implement HRTF. Not tested yet. Co-Authored-By: fuerbringer * Made HRTF work. Co-Authored-By: fuerbringer * Add 3D plot for visualizing hrtf Co-Authored-By: fuerbringer * Fixed an error concering rotation of head Co-Authored-By: fuerbringer * 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 * Corrected ear offset Co-Authored-By: fuerbringer * Update README.md * Update README.md * Update README.md * Cleanup Co-Authored-By: fuerbringer * Fixed Azimuth Co-Authored-By: fuerbringer * Fix elevation, not done yet Fixed bug in azimuth Co-Authored-By: fuerbringer * Write HRTF calculation tests to verify results. Elevation still needs some work. Co-Authored-By: fuerbringer * Begin refactoring Co-Authored-By: fuerbringer * Massive refactoring. Consolidated all scripts and improved UX massively. All functionality is now in gpuRIR module. Co-Authored-By: fuerbringer * Changed butterworth filters to sos. Co-Authored-By: fuerbringer * Make butterworth helper class implement sos in freq_dep_abs_coeff Co-Authored-By: fuerbringer * Fix elevation, complete test cases Co-Authored-By: fuerbringer * Write tests for binaural receiver Fix pinna offset errors Co-Authored-By: fuerbringer * Comments and cleanup Co-Authored-By: fuerbringer * Refactoring polar_plot.py Extending room parameters with beta Co-Authored-By: fuerbringer * Update readme Co-Authored-By: fuerbringer * Additions to documentation / comments Co-Authored-By: fuerbringer * More comments, renamed variables for consistency Co-Authored-By: fuerbringer * Fixed bug with spectrogram generation Further refactoring Co-Authored-By: fuerbringer * Fixed bug with spectrogram generation Further refactoring More comments Co-Authored-By: fuerbringer * Added further comments Co-Authored-By: fuerbringer * Added bandpass/hi-low-pass switch Co-Authored-By: fuerbringer * Update readme. * Small corrections * Reverted stereo_filters.py to standard settings. * Prepare source directivity for upstream pull request. Co-Authored-By: fuerbringer * Remove gpuRIR.extensions imports not yet in this branch. Co-Authored-By: fuerbringer Co-authored-by: fuerbringer Co-authored-by: corrooli Co-authored-by: TheBlueFireFox Co-authored-by: fuerbsev --- .gitignore | 11 +++ README.md | 10 +- examples/polar_plots.py | 156 ++++++++++++++++++++++++++++++++ examples/simulate_trajectory.py | 1 - gpuRIR/__init__.py | 33 +++++-- src/gpuRIR_cuda.cu | 62 +++++++++---- src/gpuRIR_cuda.h | 6 +- src/python_bind.cpp | 21 ++++- 8 files changed, 264 insertions(+), 36 deletions(-) create mode 100644 examples/polar_plots.py diff --git a/.gitignore b/.gitignore index 567565c..7a52c51 100644 --- a/.gitignore +++ b/.gitignore @@ -37,6 +37,17 @@ ehthumbs.db Thumbs.db +# Outputs and graphs # +###################### +*.png +*.wav + +# Generated by VSCode # +####################### +config.py +objectdb +*.pyc + # Folders # ########### build/ diff --git a/README.md b/README.md index 6b8d802..1b5e80e 100644 --- a/README.md +++ b/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. diff --git a/examples/polar_plots.py b/examples/polar_plots.py new file mode 100644 index 0000000..28e29b7 --- /dev/null +++ b/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() diff --git a/examples/simulate_trajectory.py b/examples/simulate_trajectory.py index 7977d4e..fff93e7 100644 --- a/examples/simulate_trajectory.py +++ b/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) diff --git a/gpuRIR/__init__.py b/gpuRIR/__init__.py index 73554c7..9d21ef0 100644 --- a/gpuRIR/__init__.py +++ b/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. diff --git a/src/gpuRIR_cuda.cu b/src/gpuRIR_cuda.cu index 779dc8c..22328d3 100644 --- a/src/gpuRIR_cuda.cu +++ b/src/gpuRIR_cuda.cu @@ -129,7 +129,7 @@ __device__ __forceinline__ float SabineT60( float room_sz_x, float room_sz_y, fl return 0.161f * V / Sa; } -__device__ __forceinline__ float mic_directivity(float doaVec[3], float orVec[3], micPattern pattern) { +__device__ __forceinline__ float directivity(float doaVec[3], float orVec[3], polarPattern pattern) { if (pattern == DIR_OMNI) return 1.0f; float cosTheta = doaVec[0]*orVec[0] + doaVec[1]*orVec[1] + doaVec[2]*orVec[2]; @@ -142,7 +142,7 @@ __device__ __forceinline__ float mic_directivity(float doaVec[3], float orVec[3] case DIR_HYPCARD: return 0.25f + 0.75f*cosTheta; case DIR_SUBCARD: return 0.75f + 0.25f*cosTheta; case DIR_BIDIR: return cosTheta; - default: printf("Invalid microphone pattern.\n"); return 0.0f; + default: printf("Invalid polar pattern.\n"); return 0.0f; } } @@ -245,14 +245,14 @@ __device__ __forceinline__ float image_sample_lut(float amp, float tau, float t, __global__ void calcAmpTau_kernel(float* g_amp /*[M_src]M_rcv][nb_img_x][nb_img_y][nb_img_z]*/, float* g_tau /*[M_src]M_rcv][nb_img_x][nb_img_y][nb_img_z]*/, float* g_tau_dp /*[M_src]M_rcv]*/, - float* g_pos_src/*[M_src][3]*/, float* g_pos_rcv/*[M_rcv][3]*/, float* g_orV_rcv/*[M_rcv][3]*/, - micPattern mic_pattern, float room_sz_x, float room_sz_y, float room_sz_z, + float* g_pos_src/*[M_src][3]*/, float* g_pos_rcv/*[M_rcv][3]*/, float* g_orV_src/*[M_src][3]*/, float* g_orV_rcv/*[M_rcv][3]*/, + polarPattern spkr_pattern, polarPattern mic_pattern, float room_sz_x, float room_sz_y, float room_sz_z, float beta_x1, float beta_x2, float beta_y1, float beta_y2, float beta_z1, float beta_z2, int nb_img_x, int nb_img_y, int nb_img_z, int M_src, int M_rcv, float c, float Fs) { extern __shared__ float sdata[]; - + int n[3]; n[0] = blockIdx.x * blockDim.x + threadIdx.x; n[1] = blockIdx.y * blockDim.y + threadIdx.y; @@ -308,6 +308,16 @@ __global__ void calcAmpTau_kernel(float* g_amp /*[M_src]M_rcv][nb_img_x][nb_img_ sh_orV_rcv[m*3+2] = g_orV_rcv[m*3+2]; } } + + // Copy g_orV_src to shared memory + float* sh_orV_src = &sh_orV_rcv[M_src*3]; + if (threadIdx.x==0 && threadIdx.y==0) { + for (int m=threadIdx.z; m -1 * orV_src. If rflx_id = 0 => 1 * orV_src + im_src_to_rcv[d] = -vec[d]; // change vector direction (mirror through origin) } dist = sqrtf(dist); float amp = rflx_att / (4*PI*dist); - amp *= mic_directivity(vec, &sh_orV_rcv[m_rcv], mic_pattern); + amp *= directivity(vec, &sh_orV_rcv[m_rcv], mic_pattern); // apply microphone directivity dampening + amp *= directivity(im_src_to_rcv, orV_src, spkr_pattern); // apply speaker directivity dampening + g_amp[m_src*M_rcv*prodN + m_rcv*prodN + n_idx] = amp; g_tau[m_src*M_rcv*prodN + m_rcv*prodN + n_idx] = dist / c * Fs; - if (direct_path) g_tau_dp[m_src*M_rcv + m_rcv] = dist / c * Fs; + if (direct_path){ + g_tau_dp[m_src*M_rcv + m_rcv] = dist / c * Fs; + // printf("%f,",amp); // For creating polar pattern diagrams + } } } } @@ -739,19 +761,23 @@ int gpuRIR_cuda::PadData(float *signal, float **padded_signal, int segment_len, /***********************/ float* gpuRIR_cuda::cuda_simulateRIR(float room_sz[3], float beta[6], float* h_pos_src, int M_src, - float* h_pos_rcv, float* h_orV_rcv, micPattern mic_pattern, int M_rcv, int nb_img[3], + float* h_pos_rcv, float* h_orV_src, float* h_orV_rcv, + polarPattern spkr_pattern, polarPattern mic_pattern, int M_rcv, int nb_img[3], float Tdiff, float Tmax, float Fs, float c) { // function float* cuda_simulateRIR(float room_sz[3], float beta[6], float* h_pos_src, int M_src, - // float* h_pos_rcv, float* h_orV_rcv, micPattern mic_pattern, int M_rcv, int nb_img[3], - // float Tdiff, float Tmax, float Fs, float c); + // float* h_pos_rcv, float* h_orV_src, float* h_orV_rcv, + // polarPattern spkr_pattern, polarPattern mic_pattern, int M_rcv, int nb_img[3], + // float Tdiff, float Tmax, float Fs, float c); // Input parameters: // float room_sz[3] : Size of the room [m] // float beta[6] : Reflection coefficients [beta_x1 beta_x2 beta_y1 beta_y2 beta_z1 beta_z2] // float* h_pos_src : M_src x 3 matrix with the positions of the sources [m] // int M_src : Number of sources // float* h_pos_rcv : M_rcv x 3 matrix with the positions of the receivers [m] + // float* h_orV_src : M_rcv x 3 matrix with vectors pointing in the same direction than the source // float* h_orV_rcv : M_rcv x 3 matrix with vectors pointing in the same direction than the receivers - // micPattern mic_pattern : Polar pattern of the receivers (see gpuRIR_cuda.h) + // polarPattern spkr_pattern : Polar pattern of the sources (see gpuRIR_cuda.h) + // polarPattern mic_pattern : Polar pattern of the receivers (see gpuRIR_cuda.h) // int M_rcv : Number of receivers // int nb_img[3] : Number of sources in each dimension // float Tdiff : Time when the ISM is replaced by a diffusse reverberation model [s] @@ -760,12 +786,14 @@ float* gpuRIR_cuda::cuda_simulateRIR(float room_sz[3], float beta[6], float* h_p // float c : Speed of sound [m/s] // Copy host memory to GPU - float *pos_src, *pos_rcv, *orV_rcv; + float *pos_src, *pos_rcv, *orV_src, *orV_rcv; gpuErrchk( cudaMalloc(&pos_src, M_src*3*sizeof(float)) ); gpuErrchk( cudaMalloc(&pos_rcv, M_rcv*3*sizeof(float)) ); + gpuErrchk( cudaMalloc(&orV_src, M_src*3*sizeof(float)) ); gpuErrchk( cudaMalloc(&orV_rcv, M_rcv*3*sizeof(float)) ); gpuErrchk( cudaMemcpy(pos_src, h_pos_src, M_src*3*sizeof(float), cudaMemcpyHostToDevice ) ); gpuErrchk( cudaMemcpy(pos_rcv, h_pos_rcv, M_rcv*3*sizeof(float), cudaMemcpyHostToDevice ) ); + gpuErrchk( cudaMemcpy(orV_src, h_orV_src, M_src*3*sizeof(float), cudaMemcpyHostToDevice ) ); gpuErrchk( cudaMemcpy(orV_rcv, h_orV_rcv, M_rcv*3*sizeof(float), cudaMemcpyHostToDevice ) ); @@ -785,7 +813,7 @@ float* gpuRIR_cuda::cuda_simulateRIR(float room_sz[3], float beta[6], float* h_p calcAmpTau_kernel<<>> ( amp, tau, tau_dp, - pos_src, pos_rcv, orV_rcv, mic_pattern, + pos_src, pos_rcv, orV_src, orV_rcv, spkr_pattern, mic_pattern, room_sz[0], room_sz[1], room_sz[2], beta[0], beta[1], beta[2], beta[3], beta[4], beta[5], nb_img[0], nb_img[1], nb_img[2], diff --git a/src/gpuRIR_cuda.h b/src/gpuRIR_cuda.h index 21b7f6c..abe0e48 100644 --- a/src/gpuRIR_cuda.h +++ b/src/gpuRIR_cuda.h @@ -1,6 +1,6 @@ -// Accepted polar patterns for the receivers: -typedef int micPattern; +// Accepted polar patterns for the receivers and sources: +typedef int polarPattern; #define DIR_OMNI 0 #define DIR_HOMNI 1 #define DIR_CARD 2 @@ -17,7 +17,7 @@ class gpuRIR_cuda { public: gpuRIR_cuda(bool, bool); - float* cuda_simulateRIR(float[3], float[6], float*, int, float*, float*, micPattern, int, int[3], float, float, float, float); + float* cuda_simulateRIR(float[3], float[6], float*, int, float*, float*, float*, polarPattern, polarPattern, int, int[3], float, float, float, float); float* cuda_convolutions(float*, int, int, float*, int, int); bool activate_mixed_precision(bool); bool activate_lut(bool); diff --git a/src/python_bind.cpp b/src/python_bind.cpp index e1c0682..09bd25d 100644 --- a/src/python_bind.cpp +++ b/src/python_bind.cpp @@ -14,7 +14,7 @@ class gpuRIR_bind { public: gpuRIR_bind(bool mPrecision=false, bool lut=true) : mixed_precision(mPrecision), lookup_table(lut), gpuRIR_cuda_simulator(mPrecision, lut) {}; - py::array simulateRIR_bind(std::vector, std::vector, py::array_t, py::array_t, py::array_t, micPattern, std::vector ,float, float, float, float); + py::array simulateRIR_bind(std::vector, std::vector, py::array_t, py::array_t, py::array_t, py::array_t, polarPattern, polarPattern, std::vector ,float, float, float, float); py::array gpu_conv(py::array_t, py::array_t); bool activate_mixed_precision_bind(bool); bool activate_lut_bind(bool); @@ -30,8 +30,10 @@ py::array gpuRIR_bind::simulateRIR_bind(std::vector room_sz, // Size of t std::vector beta, // Reflection coefficients py::array_t pos_src, // positions of the sources [m] py::array_t pos_rcv, // positions of the receivers [m] + py::array_t orV_src, // orientation of the sources py::array_t orV_rcv, // orientation of the receivers - micPattern mic_pattern, // Polar pattern of the receivers (see gpuRIR_cuda.h) + polarPattern spkr_pattern, // Polar pattern of the sources (see gpuRIR_cuda.h) + polarPattern mic_pattern, // Polar pattern of the receivers (see gpuRIR_cuda.h) std::vector nb_img, // Number of sources in each dimension float Tdiff, // Time when the ISM is replaced by a diffusse reverberation model [s] float Tmax, // RIRs length [s] @@ -41,6 +43,7 @@ py::array gpuRIR_bind::simulateRIR_bind(std::vector room_sz, // Size of t { py::buffer_info info_pos_src = pos_src.request(); py::buffer_info info_pos_rcv = pos_rcv.request(); + py::buffer_info info_orV_src = orV_src.request(); py::buffer_info info_orV_rcv = orV_rcv.request(); // Check input sizes @@ -49,17 +52,25 @@ py::array gpuRIR_bind::simulateRIR_bind(std::vector room_sz, // Size of t assert(nb_img.size() == 3); assert(pos_src.ndim == 2); assert(pos_rcv.ndim == 2); + assert(orV_src.ndim == 2); assert(orV_rcv.ndim == 2); assert(info_pos_src.shape[1] == 3); assert(info_pos_rcv.shape[1] == 3); + assert(info_orV_src.shape[1] == 3); assert(info_orV_rcv.shape[1] == 3); assert(info_pos_rcv.shape[0] == info_orV_rcv.shape[0]); int M_src = info_pos_src.shape[0]; int M_rcv = info_pos_rcv.shape[0]; float* rir = gpuRIR_cuda_simulator.cuda_simulateRIR(&room_sz[0], &beta[0], - (float*) info_pos_src.ptr, M_src, - (float*) info_pos_rcv.ptr, (float*) info_orV_rcv.ptr, mic_pattern, M_rcv, + (float*) info_pos_src.ptr, + M_src, + (float*) info_pos_rcv.ptr, + (float*) info_orV_src.ptr, + (float*) info_orV_rcv.ptr, + spkr_pattern, + mic_pattern, + M_rcv, &nb_img[0], Tdiff, Tmax, Fs, c); py::capsule free_when_done(rir, [](void *f) { @@ -122,7 +133,7 @@ PYBIND11_MODULE(gpuRIR_bind,m) py::class_(m, "gpuRIR_bind") .def(py::init(), py::arg("mixed_precision")=false) .def("simulateRIR_bind", &gpuRIR_bind::simulateRIR_bind, "RIR simulation", py::arg("room_sz"), py::arg("beta"), py::arg("pos_src"), - py::arg("pos_rcv"), py::arg("orV_rcv"), py::arg("mic_pattern"), py::arg("nb_img"), py::arg("Tdiff"), py::arg("Tmax"), + py::arg("pos_rcv"), py::arg("orV_src"), py::arg("orV_rcv"), py::arg("spkr_pattern"), py::arg("mic_pattern"), py::arg("nb_img"), py::arg("Tdiff"), py::arg("Tmax"), py::arg("Fs"), py::arg("c")=343.0f ) .def("gpu_conv", &gpuRIR_bind::gpu_conv, "Batched convolution using FFTs in GPU", py::arg("source_segments"), py::arg("RIR")) .def("activate_mixed_precision_bind", &gpuRIR_bind::activate_mixed_precision_bind, "Activate the mixed precision mode, only for Pascal GPU architecture or superior",