Room Impulse Response Simulation with Stochastic Ray Tracing
Room impulse response simulation aims to model the reverberant properties of a space without having to perform acoustic measurements. Many geometric and wave-based room acoustic simulation methods exist in the literature [1].
The image-source method is a popular geometric method (for an example, see Room Impulse Response Simulation with the Image-Source Method and HRTF Interpolation). One drawback of the image-source method is that it only models specular reflections between a transmitter and a receiver. There are other geometric methods that address this limitation by also taking sound diffusion into account. Stochastic ray tracing is one such method.
First, this example shows how to use acousticRoomResponse
to perform stochastic ray tracing in a room in one line of code. Then, the example showcases the stochastic ray tracing algorithm step by step for a simple "shoebox" (cuboid) room.
Use acousticRoomResponse to Perform Acoustic Ray Tracing
Synthesize Impulse Response of Shoebox Room
Start with performing ray tracing on a shoebox room with acousticRoomResponse
.
Define Room
Define the room dimensions, in meters (width, length, and height, respectively).
roomDimensions = [10 8 4];
Treat the transmitter and receiver as points within the space of the room.
tx = [2 2 2]; rx = [5 5 1.8];
Plot the room space along with the receiver (red circle) and transmitter (blue x).
h = figure; plotRoom(roomDimensions,rx,tx,h)
Define Reflection and Scattering Coefficients
A sound ray is reflected when it hits a surface. The reflection is a combination of a specular component and a diffused component. The relative strength of each component is determined by the reflection and scattering coefficients of the surfaces.
Define the absorption coefficients of the walls. The absorption coefficient is a measure of how much sound is absorbed (rather than reflected) when hitting a surface.
The frequency-dependent absorption coefficients are defined at the frequencies in the variable FVect
[4]. In this section, you use the same coefficients for all walls in the room.
FVect = [125 250 500 1000 2000 4000 8000]; A = [.02 .02 .03 .03 .04 .05 .05];
Define the frequency-dependent scattering coefficients [5]. The scattering coefficient is defined as one minus the ratio between the specular reflected acoustic energy and the total reflected acoustic energy. In this section, you use the same coefficients for all walls in the room.
D = [.13 .56 .95 .95 .95 .95 0.95];
Synthesize Impulse Response
Synthesize and plot the impulse response using acousticRoomResponse
. Specify the algorithm as "stochastic ray tracing". Use a sampling rate of 44.1 kHz.
fs = 44.1e3; ir = acousticRoomResponse(roomDimensions,tx,rx,... SampleRate=fs,... Algorithm="stochastic ray tracing",... BandCenterFrequencies=FVect,... MaterialAbsorption=A,... MaterialScattering=D); figure plot((1/fs)*(0:numel(ir)-1),ir) grid on xlabel("Time (s)") ylabel("Impulse Response")
Auralize
Apply the synthesized impulse response to an audio signal.
[audioIn,fs] = audioread("FunkyDrums-44p1-stereo-25secs.mp3");
audioIn = audioIn(:,1);
Simulate the received audio by filtering with the impulse response.
audioOut = filter(ir,1,audioIn); audioOut = audioOut/max(audioOut);
Listen to a few seconds of the original audio.
T = 10; sound(audioIn(1:T*fs),fs) pause(T)
Listen to a few seconds of the received audio.
sound(audioOut(1:T*fs),fs) pause(T)
Synthesize Impulse Response of Non-Shoebox Room
Synthesize the impulse response of a more complex, non-shoebox room.
3-D model information is often stored in stereolithography (STL) files.
Load the room setup from a sample STL file.
room = stlread("room.stl");
Visualize the room along with the source and receiver positions. Notice that the room includes a partition wall. In this example, there is no line-of-sight between the transmitter and the receiver.
tx = [1 1 1.8]; rx = [3.5 3 1.7]; figure visualizeGeneralRoom(room,tx,rx)
Synthesize the impulse response using acousticRoomResponse
.
ir = acousticRoomResponse(room,tx,rx,... SampleRate=fs,... Algorithm="stochastic ray tracing",... BandCenterFrequencies=FVect,... MaterialAbsorption=A,... MaterialScattering=D); figure plot((1/fs)*(0:numel(ir)-1),ir) grid on xlabel("Time (s)") ylabel("Impulse Response")
Simulate the received audio by filtering with the impulse response.
audioOut = filter(ir,1,audioIn); audioOut = audioOut/max(audioOut);
Listen to a few seconds of the received audio.
sound(audioOut(1:T*fs),fs) pause(T)
Construct Room with Triangulation Objects
In the previous section, you loaded a room from an STL file.
You can also construct a custom room with triangulation objects. You can assign different custom material characteristics to each wall or section in the room.
Start with the room exterior.
roomDimensions = [5 4 6]; X = roomDimensions(1); Y = roomDimensions(2); Z = roomDimensions(3); pts = [0 0 0; X 0 0; 0 Y 0; X Y 0; 0 0 Z; 0 Y Z; X 0 Z; X Y Z]; dt = delaunayTriangulation(pts(:,1), pts(:,2), pts(:,3)); [conn, verts] = freeBoundary(dt);
Define the absorption and scattering coefficients of the room. Use the same coefficients for all the walls in the room.
roomAbsorption = repmat([0.02 0.02 0.03 0.03 0.04 0.05 0.05],12,1); roomScattering = repmat([.01 .28 .7 .7 .74 .5 .5],12,1);
Next, design a partition wall in the room.
wallDimensions = [4 .2 4.5]; X = wallDimensions(1); Y = wallDimensions(2); Z = wallDimensions(3); pts = [0 0 0; X 0 0; 0 Y 0; X Y 0; 0 0 Z; 0 Y Z; X 0 Z; X Y Z] + [0 2 0]; dt2 = delaunayTriangulation(pts(:,1), pts(:,2), pts(:,3)); [conn2, verts2] = freeBoundary(dt2); conn2 = conn2 + size(verts,1);
Define the absorption and scattering coefficients of the partition wall.
wallAbsorption = repmat([0.02 0.03 0.04 0.05 0.04 0.03 0.02],12,1); wallScattering = repmat([.01 .28 .7 .7 .74 .5 .5],12,1);
Combine the room and the partition wall.
room = triangulation([conn;conn2],[verts;verts2])
room = triangulation with properties: Points: [16×3 double] ConnectivityList: [24×3 double]
Visualize the room using the helper function visualizeGeneralRoom
.
visualizeGeneralRoom(room,tx,rx)
Combine the scattering and absorption coefficients of the room and the partition wall.
roomAbsorption = [roomAbsorption;wallAbsorption]; roomScattering = [roomScattering;wallScattering];
Synthesize and plot the impulse response of the room using acousticRoomResponse
.
ir = acousticRoomResponse(room,tx,rx,... Algorithm="stochastic ray tracing",... BandCenterFrequencies=[125 250 500 1000 2000 4000 8000],... MaterialAbsorption=roomAbsorption,... MaterialScattering=roomScattering); t = (0:numel(ir)-1)/fs; figure plot(t,ir) grid on xlabel("Time (s)") ylabel("Impulse Response")
Simulate the received audio by filtering with the impulse response.
audioOut = filter(ir,1,audioIn); audioOut = audioOut/max(audioOut);
Listen to a few seconds of the received audio.
sound(audioOut(1:T*fs),fs)
Stochastic Ray Tracing Algorithm
In this section, you learn how acousticRoomResponse
operates by implementing the ray tracing algorithm step-by-step. For simplicity, you use a simple shoebox room.
Ray tracing assumes that sound energy travels around the room in rays. The rays start at the sound source, and are emitted in all directions, following a uniform random distribution. In this section, you follow (trace) rays as they bounce off obstacles (walls, floor and ceiling) in the room. You treat each ray wall incidence point as a secondary source which transmits some of its energy to the receiver. You use this energy to update a frequency-dependent histogram. You then compute the room impulse response by weighting a Poisson random process by the histogram values [2].
Define Room
Define the room dimensions, in meters (width, length, and height, respectively).
roomDimensions = [10 8 4];
Define the transmitter and receiver locations. Assume that the receiver is a sphere with radius 8.75 cm.
tx = [2 2 2]; rx = [5 5 1.8]; r = 0.0875;
Plot the room space along with the receiver (red circle) and transmitter (blue x).
h = figure; plotRoom(roomDimensions,rx,tx,h)
Define Reflection and Scattering Coefficients
Define the absorption coefficients of the walls. Each row represents the coefficient of an individual wall.
FVect = [125 250 500 1000 2000 4000 8000]; A = [.02, .02, .03, .03, .04, .05, .05; .02, .02, .03, .03, .03, .04, .04; .02, .02, .03, .03, .04, .07, .07; .02, .02, .03, .03, .04, .05, .05; .02, .02, .03, .03, .04, .05, .05; .02, .02, .03, .03, .04, .05, .05];
Derive the reflection coefficients of the six surfaces from the absorption coefficients.
R = 1-A;
Define the frequency-dependent scattering coefficients . Each row represents the coefficient of an individual wall.
D = [.13, .56, .95, .95, .95, .95, 0.95; .5, .9, .95, .95, .95, .95, 0.95; .5, .9, .95, .95, .95, .95, 0.95; .13, .56, .95, .95, .95, .95, 0.95; .13, .56, .95, .95, .95, .95, 0.95; .13, .56, .95, .95, .95, .95, 0.95];
Generate Random Rays
First, generate rays emanating from the source in random directions.
Set the number of rays.
N = 5000;
Generate the rays using the helper function RandSampleSphere
. rays
is a N-by-3 matrix. Each row of rays
holds the three-dimensional ray vector direction.
rng(0) rays = RandSampleSphere(N); size(rays)
ans = 1×2
5000 3
Initialize Energy Histogram
As you trace the rays, you update a two-dimensional histogram of the energy detected at the receiver. The histogram records values along time and frequency.
Set the histogram time resolution, in seconds. The time resolution is typically much larger than the inverse of the audio sample rate.
histTimeStep = 0.0040;
Compute the number of histogram time bins. In this example, limit the impulse response length to half a second.
impResTime = 0.5; nTBins = round(impResTime/histTimeStep);
The ray tracing algorithm is frequency-selective. In this example, focus on six frequency bands, centered around the frequencies in FVect
.
The number of frequency bins in the histogram is equal to the number of frequency bands.
nFBins = length(FVect);
Initialize the histogram.
TFHist = zeros(nTBins,nFBins);
Perform Ray Tracing
Compute the received energy histogram by tracing the rays over each frequency band. When a ray hits a surface, record the amount of diffused ray energy seen at the receiver based on the diffused rain model [2]. The new ray direction upon hitting a surface is a combination of a specular reflection and a random reflection. Continue tracing the ray until its travel time exceeds the impulse response duration.
for iBand = 1:nFBins % Perform ray tracing independently for each frequency band. for iRay = 1:size(rays,1) % Select ray direction ray = rays(iRay,:); % All rays start at the source/transmitter ray_xyz = tx; % Set initial ray direction. This direction changes as the ray is % reflected off surfaces. ray_dxyz = ray; % Initialize ray travel time. Ray tracing is terminated when the % travel time exceeds the impulse response length. ray_time = 0; % Initialize the ray energy. Energy decreases when the ray hits % a surface. ray_energy = 1/size(rays,1); while (ray_time <= impResTime) % Determine the surface that the ray encounters [surfaceofimpact,displacement] = getImpactWall(ray_xyz,... ray_dxyz,roomDimensions); % Determine the distance traveled by the ray distance = sqrt(sum(displacement.^2)); % Determine the coordinates of the impact point impactCoord = ray_xyz+displacement; % Update ray location/source ray_xyz = impactCoord; % Update cumulative ray travel time c = 343; % speed of light (m/s) ray_time = ray_time+distance/c; % Apply surface reflection to ray's energy % This is the amount of energy that is not lost through % absorption. ray_energy = ray_energy*R(surfaceofimpact,iBand); % Apply diffuse reflection to ray energy % This is the fraction of energy used to determine what is % detected at the receiver rayrecv_energy = ray_energy*D(surfaceofimpact,iBand); % Adjust the amount of reflected energy ray_energy = ray_energy * (1-D(surfaceofimpact,iBand)); % Determine impact point-to-receiver direction. rayrecvvector = rx-impactCoord; % Determine the ray's time of arrival at receiver. distance = sqrt(sum(rayrecvvector.*rayrecvvector)); recv_timeofarrival = ray_time+distance/c; if recv_timeofarrival>impResTime break end % Determine amount of diffuse energy that reaches the receiver. % See (5.20) in [2]. % Compute received energy N = getWallNormalVector(surfaceofimpact); cosTheta = sum(rayrecvvector.*N)/(sqrt(sum(rayrecvvector.^2))); cosAlpha = sqrt(sum(rayrecvvector.^2)-r^2)/sqrt(sum(rayrecvvector.^2)); E = (1-cosAlpha)*2*cosTheta*rayrecv_energy; % Scale the energy by the traveled distance E = E/(recv_timeofarrival*c)^2; % Update energy histogram tbin = floor(recv_timeofarrival/histTimeStep + 0.5); TFHist(tbin,iBand) = TFHist(tbin,iBand) + E; % Compute a new direction for the ray. % Pick a random direction that is in the hemisphere of the % normal to the impact surface. d = rand(1,3); d = d/norm(d); if sum(d.*N)<0 d = -d; end % Derive the specular reflection with respect to the incident % wall ref = ray_dxyz-2*(sum(ray_dxyz.*N))*N; % Combine the specular and random components d = d/norm(d); ref = ref/norm(ref); ray_dxyz = D(surfaceofimpact,iBand)*d+(1-D(surfaceofimpact,iBand))*ref; ray_dxyz = ray_dxyz/norm(ray_dxyz); end end end
View Energy Histogram
Plot the computed frequency-dependent histogram.
figure bar(histTimeStep*(0:size(TFHist,1)-1),TFHist) grid on xlabel("Time (s)") legend(["125 Hz","250 Hz","500 Hz","1000 Hz","2000 Hz","4000 Hz"])
Generate Room Impulse Response
The energy histogram represents the envelope of the room impulse response. Synthesize the impulse response using a Poisson-distributed noise process [2].
Generate Poisson Random Process
You model sound reflections as a Poisson random process.
Define the start time of the process.
V = prod(roomDimensions);
t0 = ((2*V*log(2))/(4*pi*c^3))^(1/3); % eq 5.45 in [2]
Initialize the random process vector and the vector containing the times at which events occur.
poissonProcess = []; timeValues = [];
Create the random process.
t = t0; while (t<impResTime) timeValues = [timeValues t]; %#ok % Determine polarity. if (round(t*fs)-t*fs) < 0 poissonProcess = [poissonProcess 1]; %#ok else poissonProcess = [poissonProcess -1];%#ok end % Determine the mean event occurrence (eq 5.44 in [2]) mu = min(1e4,4*pi*c^3*t^2/V); % Determine the interval size (eq. 5.44 in [2]) deltaTA = (1/mu)*log(1/rand); % eq. 5.43 in [2]) t = t+deltaTA; end
Create a random process sampled at the specified sample rate.
randSeq = zeros(ceil(impResTime*fs),1); for index=1:length(timeValues) randSeq(round(timeValues(index)*fs)) = poissonProcess(index); end
Pass Poisson Process Through Bandpass Filters
You create the impulse response by passing the Poisson process through bandpass filters centered at the frequencies in FVect
, and then weighting the filtered signals with the received energy envelope computed in the histogram.
Define the lower and upper cutoff frequencies of the bandpass filters.
flow = FVect/2; fhigh = FVect*2;
Set the FFT length.
NFFT = 8192;
Get the frequency vector corresponding to the FFT values.
sfft = dsp.STFT(FFTLength=NFFT,FrequencyRange="onesided");
F = sfft.getFrequencyVector(fs);
Create the bandpass filters (use equation 5.46 in [2]).
RCF = zeros(length(FVect),length(F)); for index0 = 1:length(FVect) for index=1:length(F) f = F(index); if f<FVect(index0) && f>=flow(index0) RCF(index0,index) = .5*(1+cos(2*pi*f/FVect(index0))); end if f<fhigh(index0) && f>=FVect(index0) RCF(index0,index) = .5*(1-cos(2*pi*f/(2*FVect(index0)))); end end end
Plot the bandpass filter frequency responses.
figure semilogx(F,RCF(1,:)) hold on semilogx(F,RCF(2,:)) semilogx(F,RCF(3,:)) semilogx(F,RCF(4,:)) semilogx(F,RCF(5,:)) semilogx(F,RCF(6,:)) semilogx(F,RCF(6,:)) xlabel("Frequency (Hz)") ylabel("Response") grid on
RCF = fftshift(ifft(RCF,NFFT,2,'symmetric'),2);
Filter the Poisson sequence through the six bandpass filters.
y = zeros(length(randSeq),6); for index=1:length(FVect) y(:,index) = conv(randSeq,RCF(index,:),'same'); end
Combine the Filtered Sequences
Construct the impulse response by weighting the filtered random sequences sample-wise using the envelope (histogram) values.
Compute the times corresponding to the impulse response samples.
impTimes = (1/fs)*(0:size(y,1)-1);
Compute the times corresponding to the histogram bins.
hisTimes = histTimeStep/2 + histTimeStep*(0:nTBins);
Compute the weighting factors (equation 5.47 in [2]).
W = zeros(size(impTimes,2),numel(FVect)); EdgeF = getBandedgeFrequencies(FVect,fs); BW = diff(EdgeF); for k=1:size(TFHist,1) gk0 = floor((k-1)*fs*histTimeStep)+1; gk1 = floor(k*fs*histTimeStep); yy = y(gk0:gk1,:).^2; val = sqrt(TFHist(k,:)./sum(yy,1)).*sqrt(BW/(fs/2)); for iRay=gk0:gk1 W(iRay,:)= val; end end
Create the impulse response.
y = y.*W; ip = sum(y,2); ip = ip./max(abs(ip));
Auralization
Plot the impulse response.
figure plot((1/fs)*(0:numel(ip)-1),ip) grid on xlabel("Time (s)") ylabel("Impulse Response")
Apply the impulse response to an audio signal.
[audioIn,fs] = audioread("FunkyDrums-44p1-stereo-25secs.mp3");
audioIn = audioIn(:,1);
Simulate the received audio by filtering with the impulse response.
audioOut = filter(ip,1,audioIn); audioOut = audioOut/max(audioOut);
Listen to a few seconds of the original audio.
T = 10; sound(audioIn(1:T*fs),fs) pause(T)
Listen to a few seconds of the received audio.
sound(audioOut(1:T*fs),fs)
References
[1] "Overview of geometrical room acoustic modeling techniques", Lauri Savioja, Journal of the Acoustical Society of America 138, 708 (2015).
[2] "Physically Based Real-Time Auralization of Interactive Virtual Environments", Dirk Schröder, Aachen, Techn. Hochsch., Diss., 2011.
[3] "Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality", Michael Vorlander, Second Edition, Springer.
[4] https://www.acoustic.ua/st/web_absorption_data_eng.pdf
[5] "Scattering in Room Acoustics and Related Activities in ISO and AES", Jens Holger Rindel, 17th ICA Conference, Rome, Italy, September 2001.
Helper Functions
function plotRoom(roomDimensions,rx,tx,figHandle) % PLOTROOM Helper function to plot 3D room with receiver/transmitter points figure(figHandle) X = [0;roomDimensions(1);roomDimensions(1);0;0]; Y = [0;0;roomDimensions(2);roomDimensions(2);0]; Z = [0;0;0;0;0]; figure; hold on; plot3(X,Y,Z,"k",LineWidth=1.5); plot3(X,Y,Z+roomDimensions(3),"k",LineWidth=1.5); set(gca,"View",[-28,35]); for k=1:length(X)-1 plot3([X(k);X(k)],[Y(k);Y(k)],[0;roomDimensions(3)],"k",LineWidth=1.5); end grid on xlabel("X (m)") ylabel("Y (m)") zlabel("Z (m)") plot3(tx(1),tx(2),tx(3),"bx",LineWidth=2) plot3(rx(1),rx(2),rx(3),"ro",LineWidth=2) end function X=RandSampleSphere(N) % RANDSAMPLESPHERE Return random ray directions % Sample the unfolded right cylinder z = 2*rand(N,1)-1; lon = 2*pi*rand(N,1); % Convert z to latitude z(z<-1) = -1; z(z>1) = 1; lat = acos(z); % Convert spherical to rectangular co-ords s = sin(lat); x = cos(lon).*s; y = sin(lon).*s; X = [x y z]; end function [surfaceofimpact,displacement] = getImpactWall(ray_xyz,ray_dxyz,roomDims) % GETIMPACTWALL Determine which wall the ray encounters surfaceofimpact = -1; displacement = 1000; % Compute time to intersection with x-surfaces if (ray_dxyz(1) < 0) displacement = -ray_xyz(1) / ray_dxyz(1); if displacement==0 displacement=1000; end surfaceofimpact = 0; elseif (ray_dxyz(1) > 0) displacement = (roomDims(1) - ray_xyz(1)) / ray_dxyz(1); if displacement==0 displacement=1000; end surfaceofimpact = 1; end % Compute time to intersection with y-surfaces if ray_dxyz(2)<0 t = -ray_xyz(2) / ray_dxyz(2); if (t<displacement) && t>0 surfaceofimpact = 2; displacement = t; end elseif ray_dxyz(2)>0 t = (roomDims(2) - ray_xyz(2)) / ray_dxyz(2); if (t<displacement) && t>0 surfaceofimpact = 3; displacement = t; end end % Compute time to intersection with z-surfaces if ray_dxyz(3)<0 t = -ray_xyz(3) / ray_dxyz(3); if (t<displacement) && t>0 surfaceofimpact = 4; displacement = t; end elseif ray_dxyz(3)>0 t = (roomDims(3) - ray_xyz(3)) / ray_dxyz(3); if (t<displacement) && t>0 surfaceofimpact = 5; displacement = t; end end surfaceofimpact = surfaceofimpact + 1; displacement = displacement * ray_dxyz; end function N = getWallNormalVector(surfaceofimpact) % GETWALLNORMALVECTOR Get the normal vector of a surface switch surfaceofimpact case 1 N = [1 0 0]; case 2 N = [-1 0 0]; case 3 N = [0 1 0]; case 4 N = [0 -1 0]; case 5 N = [0 0 1]; case 6 N = [0 0 -1]; end end function bef = getBandedgeFrequencies(cf,fs) G = 2; BandsPerOctave = 1; fbpo = .5/BandsPerOctave; bef = [ cf(1)*G^-fbpo, ... sqrt(cf(1:end-1).*cf(2:end)), ... min(cf(end)*G^fbpo, fs/2) ]; end function visualizeGeneralRoom(tri,txinates,rxinates) figure; trisurf(tri, ... 'FaceAlpha', 0.3, ... 'FaceColor', [.5 .5 .5], ... 'EdgeColor', 'none'); view(60, 30); hold on; axis equal; grid off; xlabel('x'); ylabel('y'); zlabel('z'); % Plot edges fe = featureEdges(tri,pi/20); numEdges = size(fe, 1); pts = tri.Points; a = pts(fe(:,1),:); b = pts(fe(:,2),:); fePts = cat(1, reshape(a, 1, numEdges, 3), ... reshape(b, 1, numEdges, 3), ... nan(1, numEdges, 3)); fePts = reshape(fePts, [], 3); plot3(fePts(:, 1), fePts(:, 2), fePts(:, 3), 'k', 'LineWidth', .5); hold on tx = txinates; rx = rxinates; scatter3(tx(1), tx(2), tx(3), 'sb', 'filled'); scatter3(rx(1,1), rx(1,2), rx(1,3), 'sr', 'filled'); end