Small Concert Hall Acoustics

Designing structures and open spaces with respect to sound quality is important for concert halls, outdoor environments, and even the rooms of a house. Simulating acoustics in the high-frequency limit, where the wavelength is smaller than the geometrical features, is best done with ray acoustics.

This tutorial model shows the basic steps and principles used when setting up a model using the Ray Acoustics physics interface. In the model, the acoustics of a small concert hall is analyzed. The model setup includes an omnidirectional sound source, wall boundary conditions for specular and diffuse scattering, surface sound pressure level evaluation, analyzing the impulse response, evaluation of objective room acoustic metrics, and a reflectogram. The metric results are compared to simple analytical estimates.

Model Definition

In this model the acoustics of a generic “small concert hall” are analyzed. The model is that of a listening environment with a volume of 460 m3 and a total surface area of 390 m2. The listening environment is fitted with absorbers and diffusers; their location is not particularly optimized and does probably not follow design rules. The aim of this tutorial is to describe the important modeling steps to perform a room acoustics simulation using ray tracing. The geometry of the concert hall is depicted in Figure 1.

Figure 1: Geometry of the small concert hall.

An omnidirectional source is located at the coordinates (xs, ys, zs) near the stage. It generates a pulse that has an SPL of 100 dB at 1 m from the source. The receiver (microphone) is located at the coordinates (xr, yr, zr). These are parameters found under . The location of the receiver can be changed in postprocessing, while the location of the source needs to be changed before running the model.

The absorption properties of the various surfaces (windows, seats, diffusers, floor, entrance, walls, and absorbers) are generic values taken from Ref. 1 and 2. The data is given in octave band values and imported from the file small_concert_hall_absorption_parameters.txt file into an interpolation function (lookup table). Moreover, the (intensity) attenuation of air is imported from the file small_concert_hall_volume_absorption.txt. Note that when the attenuation coefficient m (loaded from the file) is entered in the section, it is the amplitude attenuation αext = 0.5·m that should be used.

Figure 2: Ray location and SPL after 10 ms.

Results and Discussion

The local wavefront sound pressure level (SPL) is depicted in Figure 2 after 10 ms and in Figure 3 after 20 ms propagation for the 16 kHz frequency band. At this frequency the volume attention is the largest. When the compute intensity option is selected in the Ray Acoustics interface, wavefront curvature, intensity, and SPL is calculated along each ray. They allow visualization of the (spatially) local acoustic properties. However, it is the acoustic power transported by each ray that is important when calculating the impulse response (IR) and when visualizing the sound pressure level at surfaces. This means that the compute power option should always be selected for IR computation, while the compute intensity can be turned off. Only selecting will also make the model run slightly faster and reduce the number of degrees of freedom (DOFs) solved for. The option is also necessary when analyzing the IR.

Figure 3: Ray location and field SPL after 20 ms.

The temporal impulse response (IR) for the default source location configuration used in the model is depicted in Figure 4. The frequency domain (FFT) of the IR is depicted in Figure 5. No smoothing nor windowing is used in the plot. The IR is further analyzed using the subfeature to the plot. The feature allows the computation of the objective room acoustic metrics like clarity C, definition D, early decay time EDT, center time ts, reverberation times T20, T30, T60, and the STI. The level decay curves for the 6 octave bands used in the model (computed by the subfeature) are depicted in Figure 6. The metrics can be found in the table in the model (). The values are depicted graphically in Figure 9 and Figure 10.

When an IR is reconstructed from a ray tracing simulation, information is inferred and put back into the time signal using the temporal filter kernels. The quality of the simulated IR increases with the number of rays (this model uses 20,000) as well as the frequency resolution of the absorption, scattering, and source data (this data can be difficult to get from vendors but can often be simulated). In this model, octave band resolution is used. The Impulse Response plot also allows the use of 1/3 octave and 1/6 octave frequency resolution.

Figure 4: Temporal impulse response computed at the receiver location.

Figure 5: FFT of the impulse response (no smoothing/windowing is used).

Figure 6: Level decay curves for the 6 octave bands used in the model.

The sound pressure level in a cross section located above the seating section is depicted in Figure 7 for the 16 kHz band. It is calculated using the feature, available as a subnode to all boundary conditions. In this case it is added to a transparent surface ( used as ). The feature can be added to all other walls to postprocess the SPL distribution there if necessary.

Figure 7: SPL in a cross section 60 cm above the ground at the location of the audience.

In the ray tracing method, the intensity I and RMS pressure prms of the nth ray detected by the receiver sphere is expressed as

where Vr is the receiver volume, Lr is the distance traveled by the ray inside the receiver, and Qn is the power carried by the ray (see, for example, Ref. 3). The intensity is evaluated using the expression re1dist*rac.Q/re1vol. Plotting this information in a plot as function of the arrival time yields the (discrete time) energy impulse response or reflectogram. It is plotted for the 500 Hz and the 16 kHz octave bands in Figure 8. The slope of the curves (point data) gives a visual indication of the reverberation time of the room. In Figure 8 approximate trend lines have been added manually; their slope (60 dB down which is 6 decades for log10(In)) gives an estimate of the T60 reverberation time. In this case about 0.5 s for the 16 kHz band and 1 s for the 500 Hz band. These values are seen to fit well with the computed values below.

Figure 8: The raw data of the energetic impulse response or reflectogram. The slope represents the reverberation time in the given octave band.

A comparison between the computed value of T60 and simple statistical estimates is shown in Figure 9. To plot the data from the results table the plot is used. The estimated values are calculated using the Sabine and Eyring equations used in statistical room acoustics

where V is the room volume, S is the total surface area, m is the volume attenuation, and

is the average wall absorption (see Ref. 2). The results show quantitative agreement. It is not expected that the Sabine and Eyring predictions match the more detailed ray tracing results exactly.

Selected room acoustic metrics are plotted as function of the octave band center frequency in Figure 10. The definition, the clarity, and the center time metrics are also compared to analytical estimates in the figures. The estimates are based on analytical models of the direct energy, early energy, and late energy (see Ref. 1). They are defined in the .

Definition (D or D50) gives a metric for syllable intelligibility. In this room it is above 50% for all bands which is acceptable. The clarity metric C80 is used to characterize the transparency of music, for concert halls typical values lie between -5 to +3 dB. In this case the design is not optimal for the higher octave bands. The center time is another metric that correlates to speech intelligibility. All metrics should be estimated taking their just noticeable difference (JND) into account.

Figure 9: Reverberation time estimates based on the Sabine and Eyring formula compared to the computed reverberation time.

The speech transmission index STI is a single valued metric for speech intelligibility. It is based on modulation transfer function values (14 frequencies) and 7 octave bands. The can be plotted using the subfeature. The STI is computed as a single value that combines the information in the 7 bands. To get the single metric, change the to . Making these changes this gives a STI value of 0.63 (this indicates good intelligibility). Note that the standard does specify using also the 125 and 250 Hz octave bands (they are not included in this model). The STI values plotted as function of octave band center frequency in Figure 10 are computed based on the apparent signal to noise ratio in each band. When is selected the values in each band contribute with the appropriate weighting.

Figure 10: Objective room acoustic metrics clarity, definition, reverberation times, and STI plotted as function of octave band center frequency. The clarity and definition metrics are compared to analytical estimates. The center time (not shown here) can be seen in the plots in the model.

Notes About the COMSOL Implementation

There are several options that can be selected on the node that are useful when working with ray acoustic models and especially when evaluating the impulse response.

•

When first setting up plots, it is useful to select the option as some plots take a long time to render. Another trick is to use only a few rays initially.

•

Once the plots are set up, then before running the model (with a large number of rays), select the option. Once the model has solved the plots will be rendered. This is very useful when running the model over lunch break or over night since rendering the IR plot often takes longer time than solving the model.

•

Before saving the model, remember to set the list to . Then all plots do not need to be re-rendered once the model is opened again.

References

1. M. Vorländer, Auralization, Fundamentals of Acoustics, Modeling, Simulation, Algorithms and Acoustic Virtual Reality, Springer, 2008.

2. H. Kuttruff, Room Acoustics, CRC Press, 2009.

3. Z. Xiangyang, C. Ke'an, and S. Jincai, “On the accuracy of the ray-tracing algorithms based on various sound receiver models,” Appl. Acoust., vol. 64, pp. 433–441, 2003.

Acoustics_Module/Building_and_Room_Acoustics/small_concert_hall

Modeling Instructions

This section contains the modeling instructions for the Absorptive Muffler model. They are followed by the Geometry Sequence Instructions section.

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

The geometry is set up by importing a geometry sequence. The sequence imports the small concert hall geometry and sets up several selections. The predefined selections simplify the rest of the model setup.

In the toolbar, click

Browse to the model’s Application Libraries folder and double-click the file small_concert_hall_geom_sequence.mph.

In the toolbar, click

Click the

button in the toolbar.

Import the model parameters form a file. The parameters include the band center frequency f0, the location of the source and receiver, as well as the room volume.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file small_concert_hall_parameters.txt.

Proceed and set up interpolation functions for the absorption coefficients of the different surfaces in the concert hall. The data is easily stored in one .txt file. Also define an interpolation function for the intensity attenuation of air (given at 50 % relative humidity and 20 deg. C).

Interpolation 1 (int1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click .

Browse to the model’s Application Libraries folder and double-click the file small_concert_hall_absorption_parameters.txt.

In the text field, type 1.

Click .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
a_walls	1
a_entrance	2
a_windows	3
a_floor	4
a_diffuser	5
a_seats	6
a_absorbers	7

Locate the section. From the list, choose .

Locate the section. In the text field, type Hz.

In the text field, type 1.

Interpolation 2 (int2)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click .

Browse to the model’s Application Libraries folder and double-click the file small_concert_hall_volume_absorption.txt.

Click .

In the text field, type m_air.

Locate the section. From the list, choose .

Locate the section. In the text field, type Hz.

In the text field, type 1/m.

Now import the variables that define the room acoustic quality metric estimates. They include the reverberation time (T60), based on the Sabine and Eyring estimation equations, clarity, definition, and center time. This also requires setting up integration operators for all the surfaces.

Definitions

Variables: Quality Metric Estimates

In the window, expand the node.

Right-click and choose .

In the window for , type Variables: Quality Metric Estimates in the text field.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file small_concert_hall_variables.txt.

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , type intop_windows in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 2 (intop2)

In the toolbar, click

and choose .

In the window for , type intop_seats in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 3 (intop3)

In the toolbar, click

and choose .

In the window for , type intop_diffusers in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 4 (intop4)

In the toolbar, click

and choose .

In the window for , type intop_floor in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 5 (intop5)

In the toolbar, click

and choose .

In the window for , type intop_entrance in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 6 (intop6)

In the toolbar, click

and choose .

In the window for , type intop_walls in the text field.

Locate the section. From the list, choose .

From the list, choose .

Integration 7 (intop7)

In the toolbar, click

and choose .

In the window for , type intop_absorbers in the text field.

Locate the section. From the list, choose .

From the list, choose .

Proceed to set up and define the physics and boundary conditions of the model. To compute the impulse response, it is necessary to model the power along the rays and count the reflections. In the model the intensity is also computed along the rays. The intensity represents spatially local properties of the acoustic field approximated by the rays. The model only uses a surface mesh. Propagation in the unmeshed domains requires the definition of material properties at the interface level (in the section ). Set up boundary conditions for the different walls (boundaries).

Ray Acoustics (rac)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the cext text field, type c0.

In the ρext text field, type rho0.

In the αext text field, type 0.5*m_air(f0).

Multiplication by 0.5 is necessary as the input in COMSOL is defined for the amplitude attenuation and not the intensity attenuation (as given in the interpolation data).

Locate the section. Select the check box.

Ray Properties 1

In the window, under click .

In the window for , locate the section.

In the f text field, type f0.

Wall: Walls

In the toolbar, click

and choose .

In the window for , type Wall: Walls in the text field.

Locate the section. From the list, choose .

Locate the section. In the α text field, type a_walls(f0).

Wall: Entrance

In the toolbar, click

and choose .

In the window for , type Wall: Entrance in the text field.

Locate the section. From the list, choose .

Locate the section. In the α text field, type a_entrance(f0).

Wall: Windows

In the toolbar, click

and choose .

In the window for , type Wall: Windows in the text field.

Locate the section. From the list, choose .

Locate the section. In the α text field, type a_windows(f0).

Wall: Floor

In the toolbar, click

and choose .

In the window for , type Wall: Floor in the text field.

Locate the section. From the list, choose .

Locate the section. In the α text field, type a_floor(f0).

Wall: Diffusers

In the toolbar, click

and choose .

In the window for , type Wall: Diffusers in the text field.

Locate the section. From the list, choose .

In the γs text field, type 1-s_diffuser.

Locate the section. In the αs text field, type a_diffuser(f0).

In the αd text field, type a_diffuser(f0).

In this model the scattering coefficient s_diffuser is constant across the frequency bands. It can of course also be defined as an interpolation function that depends on f0.

Wall: Seats

In the toolbar, click

and choose .

In the window for , type Wall: Seats in the text field.

Locate the section. From the list, choose .

Locate the section. In the α text field, type a_seats(f0).

Wall: Absorbers

In the toolbar, click

and choose .

In the window for , type Wall: Absorbers in the text field.

Locate the section. From the list, choose .

Locate the section. In the α text field, type a_absorbers(f0).

Wall: SPL cross section

In the toolbar, click

and choose .

In the window for , type Wall: SPL cross section in the text field.

Locate the section. From the list, choose .

The option is used here as this surface only is meant for visualizing the SPL in a cross section. Use rendering to select the surface that is located inside the concert hall.

Click the

button in the toolbar.

Select Boundary 40 only.

Sound Pressure Level Calculation 1

In the toolbar, click

and choose .

In the window for , locate the section.

Select the check box.

Release from Grid 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the qx,0 text field, type x_s.

In the qy,0 text field, type y_s.

In the qz,0 text field, type z_s.

Locate the section. From the list, choose .

In the Nw text field, type Nrays.

Locate the section. In the Psrc text field, type P0.

Ray Termination 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the Ith text field, type 1e-14[W/m^2].

Mesh 1

Free Triangular 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Size

In the window, click .

In the window for , locate the section.

From the list, choose .

Size 1

In the window, right-click and choose .

Add a finer mesh on the surface where the SPL is computed.

In the window for , locate the section.

Click

Select Boundary 40 only.

Locate the section. Click the button.

Locate the section. Select the check box.

In the associated text field, type 0.3.

Click

Proceed and solve the model by adding a parametric sweep over the center frequency variable f0. This represents the center frequency of the octave bands analyzed in this model, in order to get a broadband response. The first time you set up and solve the model it can be useful to reduce the number of rays by changing the value of the parameter Nrays to, for example, 1000. This will make solving and postprocessing faster. Remember that the quality of the results in acoustic ray tracing increase for an increasing number of rays and more narrow frequency bands (you need to have wall absorption data with the desired resolution). In the Ray Acoustics interface the impulse response plot can handle octave, 1/3 octave, and 1/6 octave data.

Study 1

Step 1: Ray Tracing

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the text field, type 0 1.6.

For optimal performance only enter 0 and the end time for the simulation. In postprocessing, when reconstructing the impulse response, additional exact time steps for all the wall reflections are used and rendered.

Parametric Sweep

In the toolbar, click

Using the parametric sweep is important as this gives the frequency resolution (here in full octaves). The ray propagation model is solved once per frequency band. The data is collected in postprocessing, by the receiver dataset and the impulse response plot, to generate the broadband impulse response.

In the window for , locate the section.

Click

In the dialog box, choose from the list.

In the text field, type 500.

In the text field, type 20000.

Click .

Solving the model takes a couple of minutes and uses less than 2 GB of RAM (depending on your hardware). This will increase for an increasing number of rays.

In the toolbar, click

Results

Ray Trajectories (rac)

In the window for , locate the section.

From the list, choose .

In the text field, type 10[ms].

In the toolbar, click

Ray Trajectories 1

In the window, expand the node, then click .

In the window for , locate the section.

Find the subsection. From the list, choose .

Color Expression 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type rac.Lp.

In the toolbar, click

This should reproduce the image in Figure 2.

Ray Trajectories (rac)

In the window, click .

In the window for , locate the section.

In the text field, type 20[ms].

In the toolbar, click

This should reproduce the image in Figure 3.

On the node you have several options that are useful when postprocessing ray tracing simulations, where rendering plots can be time consuming. This is especially true for the impulse response plots. While setting up plots, it is useful to select such that the plots are not generated every time you change an option. It is also good practice to save the plots in the model such that they are already rendered when you open your model at a later stage. Finally, once you have set up all the plots and you are ready to run the model again, it can be useful to enable . All plots will then be recomputed after the model is solved, for example, running over lunch.

In the window, click .

In the window for , locate the section.

Select the check box.

Locate the section. From the list, choose .

Now set up the dataset used for creating the impulse response and subsequent analysis.

Receiver 3D - All Bands

In the window, expand the node.

Right-click and choose .

In the window for , type Receiver 3D - All Bands in the text field.

Locate the section. Find the subsection. In the text field, type x_r.

In the text field, type y_r.

In the text field, type z_r.

Find the subsection. In the text field, type Nrays.

In the text field, type Vol.

In the text field, type dsr.

The default option for the receiver radius is good for classical room acoustics applications (like in this example). If ray tracing is used in smaller spaces, like a car cabin, it is recommended to switch to a user defined radius.

Receiver 3D - 500 Hz Band

In the window, right-click and choose .

In the window for , type Receiver 3D - 500 Hz Band in the text field.

Locate the section. From the list, choose .

In the list, select .

Receiver 3D - 16 kHz Band

In the window, right-click and choose .

In the window for , type Receiver 3D - 16 kHz Band in the text field.

Locate the section. From the list, choose .

In the list, select .

Impulse Response

In the toolbar, click

In the window for , type Impulse Response in the text field.

Locate the section. From the list, choose .

Impulse Response 1

In the toolbar, click

and choose .

To get sharper filters you can modify some settings in the section.

In the window for , click to expand the section.

In the Np text field, type 22050.

In the δ text field, type 0.001.

Rendering the impulse response will typically take one or two minutes (depending on hardware, number of rays, and frequency band resolution). Once the plot has been generated (the first time) the necessary ray data is cached and any subsequent changes/updates of the plot is faster.

In the toolbar, click

This should reproduce the impulse response shown in Figure 4. The impulse response is the most important result of this model. The signal can be exported under the node and used for further analysis in an external signal processing tool. The response is reconstructed from the ray data detected by the dataset (arrival time, power, and band center frequency). It has a default sampling frequency of 44100 Hz. This can also be changed under the section in the plot settings window.

Plot 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type small_concert_hall_impulse_response.wav.

Click to expand the section. From the list, choose .

Click to produce a wav file with the impulse response. This file could be used for auralization or analysis with other tools. Remember to disable any subfeature before exporting.

Proceed to analyze the impulse response with the subfeature. This will create a plot of the level/energy decay and a table with the objective quality metrics.

Energy Decay 1

In the window, expand the node.