Impulse Response Plot and Receiver Dataset
The impulse response (IR) in a room acoustics simulation can be analyzed by collecting the ray information using the Receiver dataset and then analyzing the results using the Impulse Response plot. Room acoustic metrics like reverberation times, early energy metrics, and speech intelligibility can also be postprocessed using the Energy Decay subfeature. The IR can be determined for sources and boundary data (absorption coefficients, scattering parameters, source directivity, volume attenuation, and so on) given in octave, 1/3 octave, or 1/6 octave bands.
Preparing a Room Acoustics Simulation
When setting up a Ray Acoustics simulation and preparing it for a room acoustics simulation where the impulse response (IR) is an important output, certain considerations and steps need to be taken. They are as follows:
Create parameters under Global Definitions>Parameters for the band center frequency and the source power, for example, f0 and P0 respectively.
Set up interpolation functions for the properties that depend on the frequency. The data can easily be stored in a Spreadsheet formatted file and interpolation can be set to Nearest neighbor. An example of the format for a .txt file that specifies absorption coefficients for some boundaries in a simulation could be (data is invented):
%f0, wall, entrance, window, floor, diffusers, seats
125 0.1 0.4 0.01 0.02 0.1 0.05
250 0.07 0.2 0.01 0.03 0.1 0.06
500 0.05 0.12 0.02 0.05 0.1 0.07
1000 0.04 0.07 0.02 0.1 0.1 0.15
2000 0.04 0.05 0.02 0.3 0.1 0.35
4000 0.05 0.05 0.03 0.5 0.1 0.4
8000 0.08 0.05 0.04 0.5 0.1 0.6
16000 0.1 0.05 0.05 0.5 0.1 0.6
Create parameters under Global Definitions>Parameters for the receiver radius and number of released rays. These two quantities are linked by the expected error in a given time interval Δt of the impulse response, see Ref. 6 for details. With the room volume V, the receiver radius R, and for an expected error of 1 dB, the number of released rays N needed is
The receiver radius can be taken as R = 0.3 m to match the width of a common seat, and a standard time interval is Δt = 0.01 s. For small rooms or car cabins, it can be useful to set a receiver radius that represents the actual size of a microphone or slightly larger. The value obtained from this calculation should be rounded up to ensure an integer number of released rays. A higher value can also be chosen and will result in increased resolution.
Make sure that either the Compute intensity and power or the Compute power option is selected in the Intensity Computation section. The Compute power option requires fewer degrees of freedom and is more robust. If it is selected, the ray intensity (rac.I and rac.logI) cannot be visualized in ray trajectories plots, but both the Sound Pressure Level Calculation on boundaries and the impulse response can be computed.
Make sure that Count reflections is selected in the Additional Variables section.
Also make sure that the parameter that you created for the band center frequency (say f0) is used as the Ray frequency under the Ray Properties node.
Set up the ray acoustics model with appropriate boundary conditions and sources. Generally, all Wall Conditions should be set to Mixed diffuse and specular reflection, with a default scattering coefficient s = 0.05 for flat surfaces.
Create a Ray Termination node and choose Bounding box, from geometry as Spatial extents of ray propagation. Select Power as Additional termination criteria and enter the threshold P0/N*1e-6. For an omnidirectional source, this expression makes rays disappear once the power they carry has dropped by 60 dB from its initial value.
To be able to postprocess the IR, a Parametric Sweep needs to be used in the study around the Ray Tracing study step. The sweep parameter should be the frequency defined in the parameters. To easily create a Parameter value list representing the center frequencies of the bands use the ISO preferred frequencies entry method.
For the Times, specified under the Ray Tracing study step, only enter 0 and the final simulation time (the final time should be slightly larger than the estimated reverberation time). COMSOL uses small internal steps that accurately account for all reflections. The so-called Extra Time Steps are used when reconstructing the IR and other data (see below).
Receiver Dataset
Use a Receiver 2D () or Receiver 3D dataset (), selected from the More 2D Datasets submenu and the More 3D Datasets submenu, respectively, to collect the data necessary to visualize the impulse response using an Impulse Response Plot. The Impulse Response plot uses the data from a Receiver dataset as input.
Data
Select the appropriate ray dataset in the Dataset list. Then select the frequency parameters that should be used by the receiver; select one frequency for a single band analysis or all for the broad band analysis. This is the typical behavior if a parametric sweep has been set up, as described in Preparing a Room Acoustics Simulation.
Finally, select if the Receiver is Local or pointing to a Receiver feature, for example, Receiver 1 (rac/rec1) referring to the first receiver feature in physics. Pointing to a given Receiver feature defined in the physics will ensure an efficient subsequent evaluation and analysis of the impulse response in the Impulse Repose plot. If the Local option is used, the intersection of all rays and the locally defined receiver (see below) are recomputed. This step may be very time-consuming for models with many rays. It does however offer the flexibility of not having a fixed receiver defined in the geometry.
Receiver
The Receiver section is only visible if the Local option is selected as the Receiver in the Data section. For this option, the receiver behaves as a (transparent) sphere with a given radius and located at a center position (the microphone needs to have a certain finite size such that there is a reasonable probability of rays interacting with it). By changing the center position, the receiver can be moved around without the need to run the simulation again. However, computation of the intersection of rays with this virtual receiver can be time-consuming.
Under Center, specify the x, y, and (3D only) z coordinates for the center of the receiver.
Under Radius, specify the radius of the receiver. From the Radius input list, choose Variable size (large room volume) to determine the radius using an expression (see below), or choose Fixed size to enter a value for the radius in the Radius field (SI unit: m). Different theories exist for the appropriate size of the receiver. For room acoustic applications, it is recommended to use a fixed radius from which the number of rays can be derived, as described in Preparing a Room Acoustics Simulation.
When selected, the following built-in expression determines the radius R of the receiver:
where you enter the values for the Number of rays, N; Room volume, V (SI unit: m2); and Source-receiver distance, dSR (SI unit: m), to determine the radius. See Ref. 19 for details on the expression. This computed receiver size is optimized for modeling of large rooms. However, it will return a different value for each source-receiver pair if there are multiple sources and receivers present in the simulation, and it will not necessarily correspond to the dimensions of a listener.
Directivity
The Directivity section is only visible if the Local option is selected as the Receiver in the Data section. From the Directivity type list, choose Omnidirectional (the default), User defined (dB), or User defined (linear) and enter a directivity expression in the Expression field.
The User defined (dB) option applies the given gain to the received signal (positive or negative gain).
The User defined (linear) option applies a linear gain. If a negative value is used the phase of the arriving ray is inverted.
The expression can depend on the ray direction vector components (rac.nix, rac.niy, rac.niz) in order to set up a complex directivity pattern.
Extra Time Steps
The Extra Time Steps section is only visible if the Local option is selected as the Receiver in the Data section. Use the default All option. The impulse response plot has been tuned to use this option in order to get accurate and fast evaluation of arrival time and ray power.
Advanced
The Advanced section is only visible if the Local option is selected as the Receiver in the Data section. For the Interpolation between time steps list always use the default Linear interpolation option. The Cubic interpolation is not applicable for the impulse response computation and can lead to erroneous results.
Under Normal variables and Other variables, if desired, you can change the default names of the created variables:
The Distance traveled, that is the distance traveled by a ray inside the receiver.
The Directivity, that can be used to visualize the expression for the directivity.
The First ray arrival time.
Evaluation and Export
The main purpose for the Receiver dataset is to be used as input for the Impulse Response Plot. However, the data generated by the dataset can also be exported and used in an external analysis tool or software (ray power and arrival times).
Right-click the dataset and select Add Data to Export, then enter the desired data. For example, export the arrival time of the rays t, the frequency of the rays rac.f, the power of the rays rac.Q, the distance traveled by the rays in receiver re1dist, and so forth. The exported data can be stored and sorted in several formats.
When working with phase alignment of loudspeakers, it can be useful to evaluate the arrival time if the first ray. This is easily done in a Derived Values>Global Evaluation feature that points to the receiver dataset. Simply evaluate the special variable re1first.
See also the Receiver 2D and Receiver 3D section in the COMSOL Multiphysics Reference Manual.
Impulse Response Plot
In room acoustics applications, the impulse response (IR) of a source-receiver configuration represents one of the most important postprocessing results. To create the plot based on the data source defined in the Receiver dataset, add an Impulse Response subnode () to a 1D Plot Group to create the impulse response plot. After defining the characteristics of the impulse response, click Plot () to create the plot. The Energy Decay () subnode can be added to perform an analysis of the impulse response to get room acoustic metrics and the energy/level decay curves. The first rendering of the IR plot can take some time. As the data is cached, any subsequent changes to the plot will be fast.
The Small Concert Hall Acoustics model is a tutorial on how to compute the impulse response using the Ray Acoustics interface. Application Library path Acoustics_Module/Building_and_Room_Acoustics/ small_concert_hall
In the ray tracing methodology it is, as discussed above, typical to characterize sources and boundaries in frequency bands. This also means that some of the frequency content is lost. This frequency content is reconstructed when determining or reconstructing the IR. For each ray that interacts with the receiver, a small piece of signal with appropriate amplitude, arrival time, phase, and frequency content is added to the IR (the filter bank kernel). The sum of all the contributions reconstruct the full IR.
Data
Select if the Source is a Dataset or Function. If Dataset is selected the plot will analyze data from a Ray Acoustics simulation collected by a Receiver dataset. If Function is selected the input data can be impulse response data stored in an interpolation function, for example, from importing a WAV file with measured data.
If the Source is set to Dataset, then from the Dataset list, choose the appropriate Receiver Dataset. Then from the Frequency interpretation list choose an interpretation of the frequency: Octave (the default), 1/3 octave, or 1/6 octave. This selections should coincide with the frequency interpretation used in boundary conditions and sources in the underlying ray acoustics simulation, as described in Preparing a Room Acoustics Simulation. The frequency content and resolution of the impulse response is based on this selection.
Expressions
Set up the variables that are necessary for reconstructing the impulse response (typically use the defaults). These variables are used to define and add the contribution of each ray, that intersects the receiver dataset, to the impulse response.
Define the power variable for the rays (SI unit: kg·m2/s2) in the Power field. The default is rac.Q.
Define the density (SI unit: kg/m3) in the Density field. The default is rac.rho.
Define the speed of sound (SI unit: m/s) in the Speed of Sound field. The default is rac.c.
Define the number of reflections undergone by the rays (SI unit: 1) in the Number of reflections field. The default is rac.Nrefl.
x-Axis Data
From the Transformation list, choose a transformation of the data on the x-axis: None (the default), for no transformation, to show the time domain signal, or Frequency spectrum. The latter will perform an FFT of the IR and depict the power spectra as function of frequency.
Advanced
The following advanced settings are available to control some aspects of the IR reconstruction as well as the underlying filter kernel used.
Define a sampling frequency fs (SI unit: Hz). The default is 44100.
Define the Zero passing length. The default is 11025.
The next two options are legacy options (to achieve behavior for COMSOL Multiphysics releases older than version 5.6).
Select Remove noncausal signal to remove the signal (force it to zero) prior to the arrival of the first ray. Using this option will alter the energy content of the early energy component of the IR.
Select Use fully randomized phase to use the legacy method for reconstructing the IR. The fully randomized phase method (see Ref. 20) can lead to unwanted cancellations of the direct sound and early energy contributions.
The following options control the filter kernel used to reconstruct the impulse response. Advanced options allow to modify the kernel, enter a user-defined analytical kernel, or a user-defined kernel based on an interpolation function. The default is to use a Brick-wall with Kaiser window filter kernel (see Ref. 21).
Select Show the filters to show and inspect the filter bank in the Graphics window, click Plot () to create the plot. To switch between time domain and frequency domain representation change the Transformation option in the x-Axis Data section.
Figure 8-6: Example of octave band filters depicted using the Show the filters option.
Select the Filter kernel to use Brick-wall with Kaiser window (the default) or User defined. When User defined is selected enter the expression for the filter kernel hfc[t]. Use and set up parameters in the Parameters table (as you would use variables). The parameter fc is reserved for the (exact) band center frequency; fr is reserved for the ray frequency (this is typically the nominal band center frequency selected in the study); fs is reserved for the sampling frequency; t is for time; and Np is the padding length. Note that per default the User defined option sets up the brick-wall with Kaiser window filter kernel for octaves.
If the default Brick-wall with Kaiser window is selected, you can control the ripple factor δ (the default is 0.05). The effect on the filter can be seen using the Show the filters option.
Coloring and Style, and Legends
You can also make change in the Coloring and Style section or add and format a legend in the Legends section.
See also the Impulse Response section tn the COMSOL Multiphysics Reference Manual.
Energy Decay and Objective Room Acoustic Metrics
For further analysis of the impulse response, add the Energy Decay () subnode to the Impulse Response plot. Once selections have been made, click Plot () to create the plot.
Display
Select the Band type as Broadband (the default) or Individual bands. For the Individual bands option select the Band frequency, either an individual band or All. For All, all the level or energy decay curves are plotted.
Select the Plot as Energy decay (the default), Level decay, or Modulation transfer function, to define the type of graphical analysis of the impulse response to be plotted in the Graphics window.
Table
Select the room acoustic metrics to be included in the results table, per default all objective quality metrics are selected. See Ref. 8 for the definitions. The options include:
Clarity: C50 and C80
Reverberation times: T20, T30, and T60
The default name of the table is Objective Quality Metrics. The table displays the results in columns with the band center frequency fc in the first column. The results table can also be located under the Results>Tables node. A plot of the results located in the table can be created using the Table Graph plot.
Custom Plots
The energy response depicted as a reflectogram or echogram can be plotted using a Ray (Plot) and depict the expression re1dist*rac.Q/re1vol. This will give the intensity received by the “microphone” for each ray.