Modelling wave propagation characteristics in an environment

p-0002The present invention relates to a method and apparatus for modelling and simulating wave propagation characteristics of an environment. In particular, but not exclusively, the present invention relates to a method for modelling wave propagation characteristics (such as acoustic characteristics) of an environment (for example a real or virtual space) using at least two modelling methods, each of which is utilised over a respective portion of a time and/or frequency domain.

p-0003Architects and composers have long known the dramatic gravitas that can be imposed upon music or spoken word through the creative use of acoustic spaces. This has been exploited for hundreds, if not thousands of years in places of worship and great concert halls by composers such as Mozart, Berlioz and Mahler. In more modern times the psychoacoustic cues that impart this sense of space are usually generated by sound engineers to give an acoustic image that differs from the listening environment into which the sound is actually delivered. Such techniques are often used in theatre, cinema and more recently, interactive multimedia applications and give the composer/audio engineer/sound designer the artistic freedom to place the sound source and listener where required with absolute and complete control.

p-0004Acousticians, engineers and practitioners have sought to better understand the underlying physical principles that govern the behaviour of sound in an enclosed or open space. This knowledge, when combined with modern computing techniques has yielded numerous models that aim to synthesize a Room Impulse Response (RIR) suitable to impart the notion of a particular acoustic environment upon an anechoic (or dry) sound. These techniques can be traced back to simple electro-mechanical devices, the first digital feedback networks, through to more modern techniques described below.

p-0005One option for acquiring RIRs is to make a direct measurement in the desired space although this can be both time consuming and logistically difficult. Another option (which is the only possibility if the RIR of a virtual space is required) is to use a computational model to produce a virtual RIR that is an accurate analogy of the equivalent real-world space.

p-0006Geometrical acoustic techniques are the most popular solution for virtual space modelling with known architectural acoustics modelling techniques making use of ray-tracing [A. Krokstad, S. Str{acute over (ø)}m &amp; S. Sr{acute over (ø)}sdal, “Calculating the Acoustical Room Response by use of a Ray Tracing Technique”, <i>J. of Sound Vibration, vol. </i>8, No. 1, pp. 118-125, 1968] and beam tracing [A. Farina, “Pyramid Tracing Vs Ray Tracing for the Simulation of Sound Propagation in Large Rooms”, <i>Computational Acoustics </i>&amp; <i>its Environmental Applications</i>, pp. 109-116, Computational Mechanics Publications, Southampton, UK, 1995] techniques as well as others to derive RIRs as well as other more general acoustic properties. This is achieved by calculating a sufficient proportion of all the possible propagation paths that exist between a sound source and receiver through the geometric interaction of a sound-ray and the surfaces present in the model. Although these techniques produce RIRs appropriate for a digital implementation using delay lines, ray-tracing techniques have a number of limitations. The discrete, clearly defined reflection patterns that result from these geometrical methods have to be convolved with Head-Related Transfer Functions (HRTF) data in order to make them suitable for auralization purposes and they do not take into account wave interference effects. As such, these geometric methods are valid for high frequencies only and are less appropriate for low frequencies where the wave based properties of sound propagation and the presence of sparsely distributed modal frequencies tend to dominate. They are further limited in their ability to successfully model diffraction effects and hence by extension, sound occlusion due to objects being present in the propagation path, resulting in potential spatialisation errors. The rays used in ray-tracing have no cross-sectional area, whereas the rays used in beam-tracing are often conical or tetrahedral in shape and expand in area as they travel away from the sound source. This allows greater geometrical coverage than ray-tracing for the same number of rays or beams, resulting in quicker detection of valid source-receiver paths.

p-0007Finite Element Models and Boundary Element Models offer iterative methods for calculating the resonant frequencies present within an enclosed space. Although accurate these methods are computationally intensive, depending on dense mesh structures to produce results across the audible spectrum. These techniques have been used to create RIRs of virtual spaces, but alternative modelling techniques that produce equally valid results with less computational overhead and greater flexibility in terms of implementation and realization are currently more common. One such method uses the multi-dimensional digital waveguide mesh (DWM).

p-0008The digital waveguide mesh (DWM) was first proposed by Van Duyne and Smith [S. A. Van Duyne and J. O. Smith, “Physical Modeling with the 2-D Digital Waveguide Mesh”, <i>Proc. Int. Computer Music Conf</i>., Tokyo, Japan, 1993, pp. 40-47] as an extension to 1-D digital waveguide sound synthesis appropriate for modeling plates and membranes, potentially leading to full 3-D object modeling. The 1-D digital waveguide model is commonly used to model string and wind instruments, an approach similar to the Kelly-Lochbaum 1-D transmission line simulation of the vocal tract. Both of these 1-D models are based on a sampled traveling-wave implementation of the d'Alembert solution to the wave equation using bi-directional digital delay lines and scattering junctions. The reader is referred to J. O. Smith, “Physical Modeling using Digital Waveguides”, <i>Computer Music Journal</i>, vol. 16, no. 4, pp. 74-87, 1992 and J. O. Smith, “Principles of digital waveguide models of musical instruments”, <i>Applications of Digital Signal Processing to Audio and Acoustics</i>, M. Kahrs and K. Brandenburg, Eds., pp. 417-466, Kluwer Academic Publishers, Boston, 1998 for a thorough treatment and discussion of this area and a full derivation of the terms involved. Additionally, the interested reader is also directed to S. Bilbao, <i>Wave and Scattering Methods for Numerical Simulation</i>, John Wiley and Sons, Ltd., Chichester, UK, 2004 where the DWM is examined within a rigorous theoretical and comparative framework for more established, yet related wave scattering numerical simulation techniques.

p-0009From the basic 1-D digital waveguide model, higher dimension mesh structures are constructed using bi-directional delay line waveguide elements and scattering junctions which act as a regular grid of spatial and temporal sampling points within the modelled domain.

p-0010In particular since the DWM was first applied to the problem of room acoustics simulation by Savioja et al. in 1994 [L. Savioja, T. J Rinne and T. Takala, “Simulation of Room Acoustics with a 3-D Finite Difference Mesh”, <i>Proc. Int. Computer Music Conf</i>., Denmark, 1994, pp. 463-466], related research has been focused in a number of areas. There has been some considerable exploration of the properties and use of mesh topologies in 2-D or 3-D. Significant effort has also been applied to the accurate simulation of DWM boundaries with a view to how the physical properties of real materials might be modelled.

p-0011Whilst significant effort has been applied to the use of 2-D or 3-D mesh topologies, it is to be understood that such modelling techniques of high quality DWM based RIR synthesis for reverberation, auralization and virtual environment modelling demands a full 3-D render at high sample rate. Unfortunately this brings with it prohibitive execution times and costly hardware to bring the run times anywhere near acceptable levels.

p-0012A further reason for the lack of DWM methods in more general room acoustics modelling is that no intuitive software exists for the non-specialist and specialist alike to experiment with. At present the generation of RIRs using DWM models requires meticulous and time consuming editing of computer code to set up all but the simplest room geometries and mesh topologies. This means that any research scientist wishing to do work in the field has a steep learning curve to overcome as they write and maintain their own code base, making it difficult to quickly try out new ideas and share knowledge.

p-0013It is an aim of the present invention to at least partly mitigate the above-mentioned problems.

p-0014It is an aim of embodiments of the present invention to provide a method for modelling wave propagation characteristics of an environment in a way that a model having an acceptable level of accuracy can be utilised in a reasonable period of time. Embodiments of the present invention provide a model able to run close to or at real time limits.

p-0015It is an aim of embodiments of the present invention to provide a method and apparatus for modelling wave propagation characteristics of an environment using a hybrid simulation technique. Use is made of a combination of high resolution modelling techniques and coarser resolution modelling techniques as appropriate. This strikes a balance between an acceptable level of resolution and thus realism whilst maintaining an acceptable level of computational time and hardware resources required.

p-0016According to a first aspect of the present invention there is provided a method for simulating wave propagation in an environment, comprising the steps of: <ul><li id="ul0001-0001" num="0000"><ul><li id="ul0002-0001" num="0016">via a first model of the environment generated with a first modelling method, simulating a first response of the environment in a time and/or frequency domain over a first portion of the domain; and</li><li id="ul0002-0002" num="0017">via a further model of the environment generated with a further modelling method, simulating a further response of the environment in the time and/or frequency domain over at least one further portion of the domain.</li></ul></li></ul>

p-0017According to a second aspect of the present invention there is provided a computer program comprising computer program code means adapted to perform all the steps of any one of claims <b>1</b> to <b>27</b> when said program is run on a computer.

p-0018Embodiments of the present invention provide a method for modelling wave propagation characteristics of an environment, such as acoustic characteristics of a bounded virtual space, using a hybrid modelling technique. One modelling method is used to model a response of the environment over a first portion of a time and/or frequency domain whilst at least one further modelling method is used to model a response of the same environment over a further portion of the time and/or frequency domain. The first portion of the time and/or frequency domain is selected so as to contain information having a highest level of perceptual relevance. Alternatively, the first modelling method used has a level of accuracy or computational efficiency higher than that of the other modelling method/s used. The further portion of the time/frequency domain can include all or part of the first portion or may be separate and distinct.

p-0019Embodiments of the present invention will now be described hereinafter, by way of example only, with reference to the accompanying drawings in which:

p-0020<figref idrefs="DRAWINGS">FIG. 1</figref> illustrates a scattering junction with connected waveguides:

p-0021<figref idrefs="DRAWINGS">FIG. 2</figref><i>a </i>illustrates two W-DWM nodes connected via a waveguide;

p-0022<figref idrefs="DRAWINGS">FIG. 2</figref><i>b </i>illustrates two K-DWM nodes connected via a waveguide;

p-0023<figref idrefs="DRAWINGS">FIG. 3</figref> illustrates a 2-D rectilinear hybrid DWM structure;

p-0024<figref idrefs="DRAWINGS">FIG. 4</figref> illustrates stages during simulating wave propagation in an environment:

p-0025<figref idrefs="DRAWINGS">FIG. 5</figref> illustrates configuration;

p-0026<figref idrefs="DRAWINGS">FIG. 6</figref> illustrates an object to be placed in an environment;

p-0027<figref idrefs="DRAWINGS">FIG. 7</figref> illustrates simulation;

p-0028<figref idrefs="DRAWINGS">FIG. 8</figref> illustrates a simulation;

p-0029<figref idrefs="DRAWINGS">FIG. 9</figref> illustrates a simulation;

p-0030<figref idrefs="DRAWINGS">FIG. 10</figref> illustrates a computer system;

p-0031<figref idrefs="DRAWINGS">FIG. 11</figref> illustrates functional elements of a computer system;

p-0032<figref idrefs="DRAWINGS">FIG. 12</figref> illustrates an environment;

p-0033In the drawings like reference numerals refer to like parts.

p-0034Throughout this specification reference will be made to the term “environment”. It is to be understood that the term environment is to be broadly construed to include any region in which waves may propagate. The environment may, for example, be a space which is bounded or open on one or more sides. The environment may be empty or may include matter such as air or fluid or solid. The environment may include a combination of these. The environment may also include identifiable physical objects such as a chair or desk located in a particular locale.

p-0035Reference will also be made throughout this specification to the term “modelling”. Modelling involves the process of producing a model of, for example, part or parts or the whole of an environment. The model may be of an environment's physical characteristics alone or may be a model of wave propagation characteristics of the environment. That is to say a model of how waves may propagate through the environment. The model may alternatively be a combination of wave propagation characteristics and physical characteristics.

p-0036Reference will be made throughout this specification to the term “wave propagation”. It is to be understood that this term is to be broadly construed so as to cover, for example, the propagation of sound waves through an environment or other forms of wave flowing or otherwise moving through a medium.

p-0037Reference is made throughout this specification to the term “simulating”. It is to be understood that this term is to be broadly construed and includes steps carried out to estimate how properties or parameters of one feature will react when under the influence of other effects.

p-0038Reference will be made throughout this specification to an “excitation source signal”. An excitation source signal is a signal that can be used to extract acoustic characteristics of an environment. When such a signal is applied to a model of wave propagation characteristics of an environment, a response from the model will contain a time domain spread or frequency response of characteristics of the model. An excitation source signal can be an impulse or swept sign wave or other type of excitation source.

p-0039It is to be understood that a “virtual space” is a model of a space. Such a model may be a representation of a real/physical space or an imaginary space. However, in both circumstances the model is not a real space, it is not tangible and is therefore virtual.

p-0040A “real/physical space” is a space that exists in real life. It is tangible. It has physical, acoustic, and many other characteristics associated with it as does any real space or object.

p-0041<figref idrefs="DRAWINGS">FIG. 1</figref> illustrates a scattering junction J with N connected waveguides which can be used to illustrate how a model of an environment can be constructed. In particular <figref idrefs="DRAWINGS">FIG. 1</figref> illustrates how a digital waveguide mesh (DWM) <b>100</b> can be constructed using bi-directional delay line waveguide elements <b>101</b> and scattering junctions <b>102</b> which can be joined together as a regular grid of spatial and temporal sampling points within a modelled environment. It will be understood that in referring to waveguides there is no need for a specific element to be identifiable. Rather features having the functionality of a waveguide are sufficient.

p-0042Reference will be made throughout this specification to the modelling of characteristics of an environment. Such an environment may be a real environment or virtual environment and may be bounded or open.

p-0043The sound pressure in a waveguide element <b>101</b> is represented by p_i, the volume velocity by v_iand the impedance of the waveguide element by Z_iwhere p_i/v_i=Z_i. The admittance Y_iis the inverse of Z_i, such that Y_i=1/Z_i. The input to a waveguide is termed p_i⁺ and the output p_i⁻. The signal p_J,i⁺ therefore represents the incoming signal to scattering junction J along the waveguide element from the opposite junction i. Similarly, the signal p_J,i⁻ represents the outgoing signal from junction J along the waveguide to the opposite junction i. <figref idrefs="DRAWINGS">FIG. 1</figref> shows the general case of a scattering junction J with N neighbours, i=1, 2, . . . , N, with each connected unit waveguide element having an associated admittance Y_i.

p-0044From the d'Alembert solution to the 1-D wave equation the total sound pressure p_Jin a waveguide element connected to junction J can also be defined as the sum of the travelling waves in this element, or alternatively as the sum of the input and output:

p-0045<br /><i>p</i>_J<i>=p</i>_J,i⁺<i>+p</i>_J,i⁻  (Eq1)

p-0046By determining that for a lossless junction J the sum of the input velocities is equal to the sum of the output velocities, and that the sound pressures in all crossing waveguides are equal, the sound pressure p_Jat junction J for N connected waveguides can be expressed as:

p-0047<maths id="MATH-US-00001" num="00001"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mi>p</mi><mi>J</mi></msub><mo>=</mo><mfrac><mrow><mn>2</mn><mo></mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo></mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo></mo><mrow><msub><mi>Y</mi><mi>i</mi></msub><mo>·</mo><msubsup><mi>p</mi><mrow><mi>J</mi><mo>,</mo><mi>i</mi></mrow><mo>+</mo></msubsup></mrow></mrow></mrow><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo></mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo></mo><msub><mi>Y</mi><mi>i</mi></msub></mrow></mfrac></mrow></mtd><mtd><mrow><mo>(</mo><mrow><mi>Eq</mi><mo></mo><mstyle><mspace width="0.8em" height="0.8ex" /></mstyle><mo></mo><mn>2</mn></mrow><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

p-0048As the waveguide elements in a DWM are equivalent to bi-directional unit-delay lines, the input to scattering junction J at time index n, p_J,i⁺ (n), is equal to the output from neighbouring junction i into the connecting waveguide at the previous time step, p_i,J⁻(n−1). Expressing this relationship in the z-domain gives:

p-0049<br /><i>p</i>_J,i⁺<i>=z</i>⁻¹<i>·p</i>_i,J⁻  (Eq3)

p-0050Hence from (Eq2) junction pressure values are calculated according to input values from immediate neighbours, output values are calculated using (Eq1) and then propagated to neighbours via the bi-directional waveguide elements, becoming inputs at the next iteration according to (Eq3). From (Eq1), (Eq2) and (Eq3), via an appropriate linear transformation, it is possible to derive an equivalent formulation in terms of junction pressure values only:

p-0051<maths id="MATH-US-00002" num="00002"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mi>p</mi><mi>J</mi></msub><mo>=</mo><mrow><mfrac><mrow><mn>2</mn><mo></mo><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo></mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo></mo><mrow><msub><mi>Y</mi><mi>i</mi></msub><mo>·</mo><msub><mi>p</mi><mi>i</mi></msub><mo>·</mo><msup><mi>z</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup></mrow></mrow></mrow><mrow><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo></mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo></mo><msub><mi>Y</mi><mi>i</mi></msub></mrow></mfrac><mo>-</mo><mrow><msub><mi>p</mi><mi>J</mi></msub><mo>·</mo><msup><mi>z</mi><mrow><mo>-</mo><mn>2</mn></mrow></msup></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mrow><mi>Eq</mi><mo></mo><mstyle><mspace width="0.8em" height="0.8ex" /></mstyle><mo></mo><mn>4</mn></mrow><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

p-0052This expression can also be derived directly from a finite difference time domain (FDTD) formulation of the wave equation. A digital waveguide model generally refers to a representation of acoustic signal propagation via two directional wave components and schemes implemented in this way are termed W-models or W-DWMs. A linear transformation of a W-DWM leads to this alternative implementation as a Kirchhoff variable DWM (K-DWM), depending on physical quantities only rather than sampled travelling-wave components. In this form, and under certain conditions, a K-DWM can be computationally equivalent to an FDTD simulation. <figref idrefs="DRAWINGS">FIG. 2</figref><i>a </i>shows a graphical representation of two adjacent W-nodes <b>201</b>, <b>202</b> connected by a waveguide element <b>203</b>, as defined by Eq1-Eq3, while <figref idrefs="DRAWINGS">FIG. 2</figref><i>b </i>shows a graphical representation of two adjacent K-nodes <b>204</b>, <b>205</b> connected by a waveguide element <b>206</b> as defined by Eq4.

p-0053W-DWMs are able to provide a better fit at the boundary of a more complex geometric structure due to the larger number of possible air-node connections and the additional modelling flexibility they offer. However, W-DWMs are computationally more inefficient than the equivalent K-DWM, relying on two independent passes through each element in the mesh data structure and greater memory resources. It is therefore advantageous to combine the fast and efficient K-DWM with the more flexible W-DWM. To enable these two different mesh types to be combined, a KW-pipe transform can be used as an interface between the two types of DWM nodes. A KW-pipe is an all-pass network that is delay-free when moving from K-node to W-node, and that introduces delay in the opposite direction from W-node to K-node, so allowing physical variable K-DWM models to be effectively connected to scattering W-DWM models. <figref idrefs="DRAWINGS">FIG. 3</figref> shows a 2-D rectilinear hybrid DWM structure in which K-DWM and W-DWM structures are incorporated and combined within the same mesh structure via a KW-pipe interface. In this hybrid DWM the W-DWM/KW-pipe <b>301</b> models a boundary <b>302</b>, while the K-DWM nodes <b>303</b> model geometrically simple areas of uniform structure such as air, or the inside of objects.

p-0054Embodiments of the present invention provide a cross-platform (Windows, Mac-Intel, Mac-PPC, Linux or others) DWM room acoustics modelling application, implemented in the style of an Integrated Development Environment (IDE) with a hierarchical data and file structure. It will be appreciated that embodiments of the present invention are not restricted to only acoustic modelling. Rather embodiments can model wave propagation characteristics generally, such as disturbance propagation through a fluid medium. A workspace contains one or more projects; each project contains the specification for the dimensions and contents of a particular real or virtual space or other such environment, defined as a Geometry. For each geometry different surface sets can be specified, various source and receiver configurations can be set up and a number of DWMs can be defined. This data can be browsed using the workspace navigator pane on a GUI or other such user interface.

p-0055<figref idrefs="DRAWINGS">FIG. 4</figref> illustrates stages of simulating wave propagation such as acoustic wave propagation through an environment in accordance with embodiments of the present invention. The modelling process includes a configuration stage <b>400</b> during which dimensions and contents of a particular environment are specified as a geometry. Subsequently, during a simulation stage <b>401</b>, the defined environment is populated with a 2-D or 3-D DWM either across a full region of the time and/or frequency domain or over a portion of the time and/or frequency domain as selected. Once the environment is defined an input signal is input at least one location in the environment and a response is measured at least one receiver location. A further portion of the time and/or frequency domain may be modelled with a 2-D DWM or Ray-trace or other user selected modelling method. A third, post processing, stage <b>402</b> may be utilised to improve the results of the simulation.

p-0056<figref idrefs="DRAWINGS">FIG. 5</figref> illustrates a configuration block <b>500</b> used in the configuration stage <b>400</b>. The configuration block may be configured in hardware or in software.

p-0057The size and shape of a room and all the objects in it are specified in terms of a number of planar surfaces together with their reflective/absorptive and diffusive properties. Standard Cartesian co-ordinates are used and so a surface is defined by specifying the locations of its corners relative to the system origin. Embodiments of the present invention can of course utilise other techniques and co-ordinate systems for defining surfaces and locations. By combining multiple surfaces into groups, complex shapes can be created and these can then be stored as Models and re-used as required. Files are written in a simple scripting language that allows variable manipulation, loops and conditional statements, enabling sophisticated models to be defined. Input parameters can also be used to control their properties. For example, a basic model of an auditorium chair <b>600</b> can be defined as open <b>601</b>, closed <b>602</b> or occupied <b>603</b> as shown in <figref idrefs="DRAWINGS">FIG. 6</figref>. The occupied chair model uses a separately defined model for a sitting person. This enables a basic model to be defined as different objects.

p-0058As shown in <figref idrefs="DRAWINGS">FIG. 5</figref>, in this way the specification for the dimensions and contents of a particular virtual or real space are provided. For each geometry different surface sets can be specified, various source and receiver configurations can be set up and a number of DWMs can be defined. A Geometry can be defined and loaded into a project from either a Script File <b>501</b> or an imported CAD File <b>502</b> exported from another computer aided design (CAD) or 3-D graphics package. In the script file, the size and shape of a room and all the objects in it are specified in terms of a number of planar surfaces together with their reflective/absorptive and diffusive properties. This is input as surface definition <b>503</b>. By combining multiple surfaces into groups, complex shapes can be created and stored as Models <b>504</b>. These, combined with the objects <b>505</b> defined as per above, are used to generate the Script File.

p-0059Materials <b>506</b> are assigned to surfaces in the geometry file, but once loaded the materials on a boundary surface may be easily changed and saved as surface sets, allowing the same geometry to have several different profiles, and the effect of using different building or furnishing materials to be explored. These may be manipulated and edited via a GUI <b>507</b>. A user extensible library of materials that can be applied to a surface (adapted from standard literature) is provided, which stores the absorption coefficients of each material for any number of incident angles at eight octave bands and the diffusion coefficient at the same octave bands for any number of surface sizes. Material properties may also be edited and new materials created from within the main program.

p-0060One or more sources <b>508</b> and one or more receivers <b>509</b> are placed according to their co-ordinates. Using the scripting language, groups of transducers can be created along with an associated geometry to produce, for example, a spaced stereo pair, a binaural RIR incorporating a dummy head, an Ambisonic first-order 4-channel B-format receiver or some higher order B-format receiver. For each source the input signal is specified in a wave or text file.

p-0061Once a virtual space has been defined the structure is populated with a model such as a 2-D or 3-D DWM. Plug-in architecture is used to allow different DWMs to be selected. A Mesh topology plug-in <b>510</b> is used to fit the desired basic N-port scattering junction structure into a cubic array representation, used internally to map out the space. Parameters such as the mesh sample rate <b>511</b> and Air Temperature/Humidity <b>512</b> variables are set. Mesh sample rate affects the valid bandwidth of the output RIRs. Air Temperature and Humidity values affect the Air Absorption Filter settings in the Post-Processing block. Once the topology is decided upon, a compatible mesh-type plug-in <b>513</b> may be selected and this refers to how the mesh is actually implemented, which may be as a K-, W-, or hybrid KW-DWM. The particular implementation of the mesh plug-in algorithm will affect the performance and ultimately the quality of the final RIR over and above theoretical considerations. Embodiments of the present invention utilise a single modelling method to model wave propagation characteristics of an environment. Other embodiments utilise a hybrid modelling method utilising one modelling method such as 3-D DWM and at least one further modelling method such as 2-D DWM or ray-tracing.

p-0062Operation of the simulation stage is illustrated in more detail with respect to <figref idrefs="DRAWINGS">FIG. 7</figref>. <figref idrefs="DRAWINGS">FIG. 7</figref> illustrates a simulation block <b>700</b> used to simulate the environment being tested. The simulation block may be configured in hardware or software.

p-0063The meshing algorithm <b>701</b> fills the room geometry with a uniform spatial grid of DWM nodes. As the sampling grid arrangement of a DWM will vary according to the topology plug-in used and the arbitrary geometry of the defined virtual space, a generalized, flexible approach to filling the space is used: <ul><li id="ul0003-0001" num="0000"><ul><li id="ul0004-0001" num="0061">1. A single air-node is placed at a user-defined seed-point within the geometry model and this node reproduces itself by sending out “creepers” along its ports into the surrounding empty space.</li><li id="ul0004-0002" num="0062">2. An identification is made to see if a boundary or existing air-node has been encountered. If not then a new air-node is created at the end of the creeper.</li><li id="ul0004-0003" num="0063">3. If a creeper discovers a surface then a boundary-node is created at the end of the creeper rather than an air-node. These boundary-nodes are incapable of reproduction and hence do not send out any creepers of their own.</li><li id="ul0004-0004" num="0064">4. If a creeper discovers another air-node then a connection between the two is established.</li><li id="ul0004-0005" num="0065">5. All new air-nodes will then send out their own creepers (and hence the space begins to fill with basic N-port air-nodes).</li></ul></li></ul>

p-0064In this manner the meshing process will continue until the edge of the modelled space, or object geometry boundaries are encountered. When a boundary-node is created it inherits the properties of the particular surface associated with its geometrical location.

p-0065Although it is possible to run a 2-D or 3-D DWM, or a ray-trace simulation directly from this point embodiments of the present invention add additional flexibility and parallel modelling options. The user is provided with information identifying system parameters based on the current system resources, with the maximum available memory input <b>703</b> determining the maximum bandwidth of the final RIRs. An input RIR Bandwidth <b>704</b> allows the user to reduce this value further, so improving execution time. An input Execution time <b>705</b> may be reduced or extended based on the time the user is prepared to wait for the RIRs to be generated. As a minimum execution time is determined by an Early/Late threshold input <b>706</b> the threshold input is set to determine where a cut in the time and/or frequency domain is made. This may be dependent upon dimensions of the environment being modelled as a larger space will take longer for direct sound and early reflections from one or more sources to reach a receiver. For a reduced bandwidth or execution time 3-D DWM simulation, the remainder of the output may be replaced by a 2-D DWM simulation with similar parameters, or a ray-trace, with quality determined by the number of rays input <b>707</b> (provided by the user) actually traced. The Input block <b>708</b> loads the appropriate excitation files into the simulation prior to execution ready for the main simulation <b>709</b> which may be carried out as a 3-D DWM <b>710</b> and/or a 2-D DWM <b>711</b> and/or a ray-trace <b>712</b> or other modelling method or any combination of these according to the selected option. By updating the state of each node according to the DWM algorithm used one sample at a time the excitation signals propagate through the mesh. Alternatively in the ray-trace case virtual sound rays are fired off into the space until they encounter a surface where they are reflected/absorbed/diffused as appropriate. A visualisation block <b>713</b> allows the DWM process to be viewed using an animation where the pressure at each node is represented by either a coloured bar or point. The ray-trace option may also be viewed if required. The simulation is terminated either when the required number of samples have been generated or when the signal level falls below a specified value.

p-0066An optional post processing stage <b>402</b> may then be utilised to process the signals output from the simulation block. An output of the post processing block provides a final complete RIR wave file for each simulation run.

p-0067<figref idrefs="DRAWINGS">FIG. 7</figref> illustrates schematically how wave propagation through an environment can be simulated. Via a graphical user interface <b>507</b>, a user first defines an environment, such as a bounded virtual space identifying surface characteristics including dimension and surface materials. This is the configuration stage <b>400</b> noted above with respect to <figref idrefs="DRAWINGS">FIGS. 4 and 5</figref>. A user then selects one or more modelling methods to produce models indicating wave propagation characteristics of the space. It is possible for a user to select only one model, for example, a high accuracy 3-D digital waveguide mesh. However, it is understood that modelling using such a modelling type is computationally expensive and time-consuming. A user may therefore utilise a hybrid modelling scheme which utilises more than one modelling type to estimate wave propagation through an environment. By identifying parts of a response which are of most significance and which are therefore modelled with a modelling type having a high level of accuracy and modelling a remainder with a modelling type of less accuracy, a compromise is reached which provides a satisfactory level of accuracy but in an acceptable level of time or having an acceptable hardware requirement. For example, a user may select to use a 3-D digital waveguide mesh <b>710</b> and a ray-trace <b>712</b>. Each model will be used to model wave propagation caused by one or more input signals over an associated portion of a time and/or frequency band. For example, it has been appreciated that typically up to the first 100 milliseconds of the impulse response of a space contains the most perceptually relevant information. That is to say, a person can determine more or less accurately information relating to a space by listening only to the first 100 milliseconds. The tail (or latter part) of a signal contains less perceptually relevant information that is less critical in terms of low level content. Therefore a model such as the 3-D digital waveguide mesh <b>710</b> is typically utilised over the first 100 milliseconds in the time domain. This portion of the time domain response contains the most critical relevant information in terms of perceiving the size and shape of the space. A ray-trace model is utilised to model the wave propagation during a remainder portion of the time domain. This portion of the time domain includes less perceptually relevant information and therefore a cruder (but quicker and less processing intensive) modelling technique can be used. It is to be understood that embodiments of the present invention can utilise two or more models to provide a hybrid modelling methodology to produce a result representing how waves will propagate in an environment. It is also to be understood that the portion of the time and/or frequency domain which is selected by a user can be chosen to strike a balance between obtaining acceptably accurate results in a reasonable time or using a reasonable amount of processing power and thus cost.

p-0068As illustrated in <figref idrefs="DRAWINGS">FIG. 8</figref>, an excitation source signal <b>801</b> is used to stimulate each model. A response from each model is then combined. The combined response provides a room impulse response RIR of the environment. The output from the models is then convolved or otherwise combined with an input signal <b>803</b> which may be an anechoic signal. An anechoic signal includes no reverberant information. Effectively the anechoic signal provides a signal identifying a sound which is to be modelled in an environment. The output from the models when convolved with that sound provide a final output <b>804</b> indicating how the anechoic signal will be affected by the environment. In the example of the anechoic signal being an acoustic signal the final output <b>804</b> will provide an indication of how that sound would be heard by a person in the selected environment. Post processing <b>402</b> may be carried out on the final output signal <b>804</b> to improve the quality of the result.

p-0069<figref idrefs="DRAWINGS">FIG. 9</figref> illustrates a further embodiment of the present invention in which a user selects two or more modelling types when implementing a hybrid modelling strategy. In <figref idrefs="DRAWINGS">FIG. 9</figref> the anechoic signal <b>901</b> is input directly into these models. Each model is used to determine a response of the environment to the anechoic signal and the response from each model is combined to provide an output response <b>904</b>. It will be appreciated that such a simulating scheme provides an alternative scheme for providing accurate results but requires a drastically increased run time or hardware which is more expensive.

p-0070<figref idrefs="DRAWINGS">FIG. 10</figref> illustrates a computer system that can be used to operate the present invention as a software computer program or hardware component. <figref idrefs="DRAWINGS">FIG. 11</figref> shows functional elements of the computer system of <figref idrefs="DRAWINGS">FIG. 10</figref> as a block diagram.

p-0071It is possible to run the present invention using a computer system or laptop computer. When the present invention is embodied as a computer program, the program can be stored within a memory device of the computer <b>1001</b>. A user can activate and control the program using a user interface such as the computer's keyboard <b>1002</b> and mouse <b>1003</b> or any other computer user interface. The program can be viewed using a visual display unit (VDU) such as a monitor <b>1004</b>. The central processing unit (CPU) <b>1101</b> of the computer will control and perform the processing required by the program as well as commands and control signals input by the user through user interfaces of the computer.

p-0072The user can define an environment to be simulated by controlling a graphical user interface (GUI) of the program using a computer user interface; the GUI can be viewed on a monitor <b>1004</b>. The CPU <b>1101</b> will process the input commands of the user and store data relating to the defined space within a memory of the computer such as the internal memory <b>1102</b>.

p-0073An input sound source can be input into the computer through the computer's audio input <b>1103</b>. The input sound can then be stored within a memory device such as the internal memory <b>1102</b> of the computer. When the user selects for a simulation of the environment to be performed, the CPU <b>1101</b> will perform the simulation using the parameters of the defined environment stored in memory as well as the input sound. The final audio output can then be stored into internal memory <b>1102</b> or stored on a peripheral memory device <b>1104</b> such as a CD-ROM. The sound before and after processing can be listened to using a sound reproduction system such as speakers <b>1105</b> or headphones. As an alternative it is possible for pre-recorded/pre-generated sounds to be stored, for example on a hard disk of a computer or a CD-ROM. As a still further alternative the input sounds may be generated by the computer in real time using a sound synthesis algorithm.

p-0074Results of a series of tests will now be described identifying advantageous features of embodiments of the present invention.

p-0075The main objective of these tests has been to firstly demonstrate how embodiments of the present invention can be used to successfully synthesize a virtual environment such that the RIRs produced are appropriate for auralization purposes. A further objective is to evaluate the computational savings offered and what affect they may have on the resultant RIRs compared with a full 3-D render. These tests are based on a model of an actual room and data will be compared with RIR measurements obtained from within this space and with a purely geometric acoustic simulation produced using a leading commercial room acoustics software package.

p-0076The room <b>1200</b> that has been selected for measurement and modelling is a rehearsal room in the Music Department at the University of York in the UK. The room is illustrated in <figref idrefs="DRAWINGS">FIG. 12</figref>. It is relatively small, with a volume of approximately 27 m³, implying that a high sample rate 3-D DWM model can be computed reasonably. It has no soft furnishings, and all items of furniture have been removed for the purposes of this experiment. The main surfaces in the room consist of painted plasterboard, cork floor tiles, standard ceiling tiles and glass windows. Hence despite being small, the room is quite bright sounding and the lack of absorbing materials, together with parallel walls, indicate a highly modal response. A plan view of the room with the single source S and receiver positions R<b>1</b>-R<b>4</b> marked as shown in <figref idrefs="DRAWINGS">FIG. 12</figref>. Both source and receivers are set at a height of 1.5 m with the floor to ceiling height of the room being 2.49 m.

p-0077RIR measurements for this test room are obtained for five different cases as follows:

Case 1. Actual Test Space.

p-0078Case 1 is achieved by using a 4-channel Soundfield SPS422B B-Format microphone in each of the four receiver positions, R<b>1</b>-R<b>4</b>, to capture the response of the space to a 15 s 22 Hz-22 kHz logarithmic sine-sweep excitation, with the loudspeaker positioned at S. Deconvolution of these responses with the inverse of the sine-sweep signal yields the actual RIR. For the purposes of these experiments only the first-order W-channel B-Format omnidirectional pressure based RIR is used.

Case 2. Geometric Acoustic Model

p-0079Based on physical measurements of the test space. Ideal impulse is applied at the input and outputs are based on a first-order W-channel omnidirectional RIR for each receiver position R<b>1</b>-R<b>4</b>. 20,000 rays were used for the simulation.

Case 3. 3-D DWM Simulation

p-0080Based on physical measurements of the test space. Input function applied for a 3-D rectilinear mesh and the outputs at each receiver point are taken from a single air-node, equivalent to a W-channel RIR. Mesh sample rate selected such that the final bandwidth is valid to 10 kHz.

Case 4. 2-D DWM Simulation

p-0081Based on physical measurements of the test space. Input function applied for a 2-D triangular mesh and the outputs at each receiver point are taken from a single air-node as before. Mesh sample rate selected such that the final bandwidth is valid to 22050 kHz.

Case 5. 3-D12-D/R-T DWM Hybrid Simulation

p-0082Based on physical measurements of the test space. Input functions applied as appropriate and the outputs at each receiver point are taken from a single air-node as before. The early part is a 3-D simulation valid to 10 kHz and is truncated at 38 ms. A ray-trace model is then used for the 10-18 kHz region for the same duration. 40,000 rays are used. The reverberant tail is based on a 2-D DWM model valid to 18 kHz. Post processing is applied to obtain the final RIR wave file. Case 5 thus illustrates a a-part simulation in accordance with an embodiment of the present invention.

p-0083For both the DWM and the geometric acoustic models building materials were initially identified as closely as possible and appropriate absorption/reflection coefficients applied. Unfortunately the results were very non-ideal and resulted in predictions for the amount of time, in seconds, for reverb to decay 60 dB below the initial value (RT₆₀) that were significantly shorter than the actual measurements. Given that the space has little absorption and demonstrates strong modal effects this is not surprising. Hence a different approach was adopted and instead DWM absorption/reflection coefficients were optimised according to octave band RT₆₀values for the actual test room and then applied to all surfaces. The values applied are shown in Table 1. Note that when the values were applied directly in the geometric acoustics package, reverberation times of 3-4 s were observed and hence values were also optimised in a similar manner. Note also that for this ‘optimised’ and hence highly reflective room, even small changes in the reflection factor values will have a potentially large influence in the final reverberation time values. In all cases, no diffusion modelling has been applied as all surfaces in the real space are hard, flat and smooth, and applied diffusion can also act to reduce RT60 values further.

p-0084<tables id="TABLE-US-00001" num="00001"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="1"><colspec colname="1" colwidth="217pt" align="center" /><thead><row><entry namest="1" nameend="1" rowsep="1">TABLE 1</entry></row></thead><tbody valign="top"><row><entry namest="1" nameend="1" align="center" rowsep="1" /></row><row><entry>Octave band reflection coefficients as applied in DWM </entry></row><row><entry>models and in the geometric acoustics/ray-tracing cases.</entry></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="3"><colspec colname="1" colwidth="77pt" align="center" /><colspec colname="2" colwidth="56pt" align="center" /><colspec colname="3" colwidth="84pt" align="center" /><tbody valign="top"><row><entry /><entry /><entry>Geometric</entry></row><row><entry>Octave</entry><entry>Input</entry><entry>Acoustic</entry></row><row><entry>Frequency</entry><entry>Reflection</entry><entry>Reflection</entry></row><row><entry>Band (Hz)</entry><entry>Coefficients</entry><entry>Coefficients</entry></row><row><entry namest="1" nameend="3" align="center" rowsep="1" /></row><row><entry> 62</entry><entry>0.99</entry><entry>0.95 </entry></row><row><entry>125</entry><entry>0.97</entry><entry>0.945</entry></row><row><entry>250</entry><entry>0.97</entry><entry>0.900</entry></row><row><entry>500</entry><entry>0.97</entry><entry>0.900</entry></row><row><entry>1 k</entry><entry>0.96</entry><entry>0.900</entry></row><row><entry>2 k</entry><entry>0.94</entry><entry>0.900</entry></row><row><entry>4 k</entry><entry>0.89</entry><entry>0.900</entry></row><row><entry>8 k</entry><entry>0.85</entry><entry>0.900</entry></row><row><entry namest="1" nameend="3" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

Results

p-0085RT₆₀values for each RIR are calculated in octave bands from the quantity T30 according to ISO3382 using the Aurora plug-in suite for Adobe Audition. Once generated these values are averaged across the four receiver positions to arrive at a final value to indicate the overall behaviour of the space. The results for each case are summarised in Table 2. Note first of all that there is a missing result for the 16 k octave band in Case 3 as the 3-D DWM example only has a valid bandwidth up to 10 kHz. This also affects the value in the 8 k band as it increases slightly to 0.44 s on an otherwise downwards trend. Also despite a reasonable match to the test room, the Case 2 values should not be relied upon in the lowest octave frequency bands as geometric acoustic modelling techniques are not capable of accurately simulating the acoustics of a room in this region.

p-0086<tables id="TABLE-US-00002" num="00002"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="1"><colspec colname="1" colwidth="217pt" align="center" /><thead><row><entry namest="1" nameend="1" rowsep="1">TABLE 2</entry></row></thead><tbody valign="top"><row><entry namest="1" nameend="1" align="center" rowsep="1" /></row><row><entry>Octave band ISO3382-T30 values for each test </entry></row><row><entry>case scenario as described in Section 6.3.</entry></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="7"><colspec colname="offset" colwidth="14pt" align="left" /><colspec colname="1" colwidth="21pt" align="center" /><colspec colname="2" colwidth="42pt" align="center" /><colspec colname="3" colwidth="28pt" align="center" /><colspec colname="4" colwidth="42pt" align="center" /><colspec colname="5" colwidth="28pt" align="center" /><colspec colname="6" colwidth="42pt" align="center" /><tbody valign="top"><row><entry /><entry>Freq.</entry><entry>Case</entry><entry>Case</entry><entry>Case</entry><entry>Case</entry><entry>Case</entry></row><row><entry /><entry>(Hz)</entry><entry>1</entry><entry>2</entry><entry>3</entry><entry>4</entry><entry>5</entry></row><row><entry /><entry namest="offset" nameend="6" align="center" rowsep="1" /></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="7"><colspec colname="offset" colwidth="14pt" align="left" /><colspec colname="1" colwidth="21pt" align="char" char="." /><colspec colname="2" colwidth="42pt" align="char" char="." /><colspec colname="3" colwidth="28pt" align="char" char="." /><colspec colname="4" colwidth="42pt" align="char" char="." /><colspec colname="5" colwidth="28pt" align="char" char="." /><colspec colname="6" colwidth="42pt" align="char" char="." /><tbody valign="top"><row><entry /><entry>31</entry><entry>1.63</entry><entry>2.55</entry><entry>1.22</entry><entry>1.20</entry><entry>1.37</entry></row><row><entry /><entry>62</entry><entry>2.20</entry><entry>2.28</entry><entry>1.31</entry><entry>1.26</entry><entry>1.32</entry></row><row><entry /><entry>125</entry><entry>1.30</entry><entry>1.77</entry><entry>1.14</entry><entry>1.24</entry><entry>1.28</entry></row><row><entry /><entry>250</entry><entry>0.97</entry><entry>1.13</entry><entry>0.99</entry><entry>1.05</entry><entry>1.09</entry></row><row><entry /><entry>500</entry><entry>0.91</entry><entry>0.92</entry><entry>0.71</entry><entry>1.00</entry><entry>1.00</entry></row><row><entry /><entry> 1 k</entry><entry>0.83</entry><entry>0.81</entry><entry>0.60</entry><entry>0.91</entry><entry>0.92</entry></row><row><entry /><entry> 2 k</entry><entry>0.75</entry><entry>0.76</entry><entry>0.48</entry><entry>0.73</entry><entry>0.74</entry></row><row><entry /><entry> 4 k</entry><entry>0.68</entry><entry>0.65</entry><entry>0.40</entry><entry>0.54</entry><entry>0.53</entry></row><row><entry /><entry> 8 k</entry><entry>0.53</entry><entry>0.52</entry><entry>0.44</entry><entry>0.46</entry><entry>0.43</entry></row><row><entry /><entry>16 k</entry><entry>0.37</entry><entry>0.49</entry><entry>—</entry><entry>0.47</entry><entry>0.41</entry></row><row><entry /><entry namest="offset" nameend="6" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0087Table 3 presents comparative performance results for each of the test cases. Clearly the hybrid solutions bring the total elapsed time required to produce a complete RIR down to reasonable levels compared with a full 3-D render, which is even in this case, only valid to 10 kHz. Note that the 2-D DWM offers a good compromise when computational resources are limited.

p-0088<tables id="TABLE-US-00003" num="00003"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="1"><colspec colname="1" colwidth="217pt" align="center" /><thead><row><entry namest="1" nameend="1" rowsep="1">TABLE 3</entry></row></thead><tbody valign="top"><row><entry namest="1" nameend="1" align="center" rowsep="1" /></row><row><entry>Performance data for each of the 4 simulated cases, </entry></row><row><entry>showing total time elapsed to produce the complete</entry></row><row><entry>required RIR data and the total system memory used.</entry></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="5"><colspec colname="1" colwidth="63pt" align="center" /><colspec colname="2" colwidth="28pt" align="center" /><colspec colname="3" colwidth="49pt" align="center" /><colspec colname="4" colwidth="28pt" align="center" /><colspec colname="5" colwidth="49pt" align="center" /><tbody valign="top"><row><entry /><entry>Case 2</entry><entry>Case 3</entry><entry>Case 4</entry><entry>Case 5</entry></row><row><entry namest="1" nameend="5" align="center" rowsep="1" /></row><row><entry>Total Time</entry><entry>00:54</entry><entry>64:56</entry><entry>00:26</entry><entry>03:24</entry></row><row><entry>(Hrs:Mins)</entry><entry /><entry /><entry /><entry /></row><row><entry>Memory</entry><entry>—</entry><entry>645</entry><entry>8</entry><entry>653</entry></row><row><entry>Used (Mb)</entry></row><row><entry namest="1" nameend="5" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0089Throughout the description and claims of this specification, the words “comprise” and “contain” and variations of the words, for example “comprising” and “comprises”, means “including but not limited to”, and is not intended to (and does not) exclude other moieties, additives, components, integers or steps.

p-0090Throughout the description and claims of this specification, the singular encompasses the plural unless the context otherwise requires. In particular, where the indefinite article is used, the specification is to be understood as contemplating plurality as well as singularity, unless the context requires otherwise.

p-0091Features, integers, characteristics, compounds, chemical moieties or groups described in conjunction with a particular aspect, embodiment or example of the invention are to be understood to be applicable to any other aspect, embodiment or example described herein unless incompatible therewith.