Thesis Ray Tracing Finland

Room acoustics modeling using the raytracing method: implementation and evaluation

Licentiate Thesis University of Turku Department of Physics 2005 David Oliva Elorza

2

3

A mi abuelo

4

Abstract This thesis discusses room acoustics modeling and the basics and drawbacks of different methods. Then it concentrates on the ray-tracing algorithm, and offers the author's implementation of the method and an evaluation of its accuracy when comparing modeled sound decays and predicted acoustical parameters with the corresponding measured data in three different enclosures. First, a general introduction of the sound propagation principles and a discussion of the state of the art in room acoustics modeling are given. How to implement the ray-tracing algorithm is explained, not only the way adopted by the author but also other possibilities which could improve accuracy. That is followed by a geometrical and an acoustical description of the three studied enclosures, and the comparison of the modeled and measured data. Sound pressure level is predicted with very good accuracy, while reverberation time is more contradictory due to the negation of the wave nature of sound which is inherent in the model. The paper ends with a discussion of the results’ accuracy and the reasons for similarities and/or differences.

5

Preface In August 2004 I started to write this thesis. In principle, it had to include the form of the ray-tracing method that I had implemented during my time in the Finnish Institute of Occupational Health (Työterveyslaitos), together with the evaluation of its accuracy in some indoor enclosures, but it turned out to include more than that. The referred method is good, although it is not perfect and it has several drawbacks. Hence, I finally ended up writing down many things which in principle were outside the scope of the thesis, but that constituted a part of the whole. That is, I discussed diffraction at great length although I had not yet implemented it, but just because I thought I might do so some day. Therefore, it became a very long paper, which could be hard to read for some readers, useless until page 56 or 70 for others, and for the rest of them utterly pointless. I do not know yet into which category my thesis reviewers belong, but maybe to one of them. Their opinion made me work a little harder on the paper, delete about twenty pages, and reorganize most of it. However, that was very positive, not only because this thesis now looks better than it did before, but also because I have improved my knowledge. A lot of time has passed since I typed the first lines of this paper, but on the other hand, I have been doing some other things since then. For instance, the form of the algorithm which is currently on my hard disk differs considerably from the one described in this work, and my expected future work, and even more has already been done. However, that is definitely another story about which I will perhaps still write some day. Just one more thing: if you decide to go on with the reading of this thesis, please, enjoy it!

6

Table of contents ABSTRACT.................................................................................................................................... 4 PREFACE....................................................................................................................................... 5 TABLE OF CONTENTS .............................................................................................................. 6 List of symbols and abbreviations .............................................................................................................8

INTRODUCTION.......................................................................................................................... 9 1. OVERVIEW OF SOUND ....................................................................................................... 13 1.1. Sound: Wave propagation.................................................................................................................13 1.2. Sound: Reflection ...............................................................................................................................16 1.3. Sound: wave effects............................................................................................................................18 1.3.1. Diffraction.................................................................................................................................18 1.3.2. Interference...............................................................................................................................19 1.4. Temporal distributions of reflections; the impulse response..........................................................22 1.5. Diffuse sound field..............................................................................................................................23 1.6. Acoustic levels.....................................................................................................................................24 1.7. Sound absorption ...............................................................................................................................25 1.7.1. Air absorption...........................................................................................................................26 1.7.2. Wall absorption ........................................................................................................................26

2. MODELING METHODS IN ROOM ACOUSTICS........................................................... 28 2.1. Modeling methods..............................................................................................................................28 2.1.1. Empirical and statistical methods...........................................................................................28 2.1.2. Wave-based methods................................................................................................................29 2.1.3. Geometrical acoustics methods ...............................................................................................30 2.1.3.1. Ray-tracing method.....................................................................................................31 2.1.3.2. Image method ..............................................................................................................33 2.1.3.3. Particles........................................................................................................................35 2.1.3.4. Fitting-zone method ....................................................................................................35 2.1.3.5. Cone or pyramids tracing method .............................................................................36 2.1.3.6. Hybrids.........................................................................................................................37 2.2. Diffusion and diffraction in the ray-tracing algorithm...................................................................38 2.2.1. Modeling diffusion in the ray-tracing method .......................................................................38 2.2.1.1. Statistical treatment of diffuse reflections. The scattering coefficient....................38 2.2.1.2. Reflection models with secondary diffuse sources....................................................43 2.2.2. Modeling diffraction in the ray-tracing method....................................................................46 2.3. Considerations about the ray-tracing algorithm.............................................................................49 2.3.1. Frequency range.......................................................................................................................49 2.3.2. Negation of wave nature, interference and diffraction .........................................................49 2.3.3. Angular distribution of reflected sound .................................................................................50 2.3.4. Geometrical description...........................................................................................................51 2.3.5. Frequency dependent factors. The absorption and the scattering coefficients ...................51 2.3.6. Number of rays .........................................................................................................................52 2.3.7. The detection problem .............................................................................................................52 2.3.8. Other problems.........................................................................................................................54 2.4. Round Robin.......................................................................................................................................54

3. IMPLEMENTATION OF THE RAY-TRACING ALGORITHM ................................... 57 3.1. Flowchart of the implemented ray-tracing algorithm ....................................................................58 3.2. Geometrical calculations ...................................................................................................................59 3.2.1. Geometrical description of the enclosure ...............................................................................59 3.2.2. Shooting rays ............................................................................................................................60 3.2.3. Detection of collisions...............................................................................................................61 3.2.3.1. Point-in-polygon test ...................................................................................................61 3.2.3.2. Visibility test ................................................................................................................63

7 3.2.4. Modeling reflections.................................................................................................................64 3.2.4.1. Specular reflection.......................................................................................................65 3.2.4.2. Diffuse reflection .........................................................................................................65 3.2.5. The image method...........................................................................................................................66 3.3. Energy calculations............................................................................................................................66 3.3.1 The instant power of every ray.................................................................................................66 3.3.2. Reception of sound rays ...........................................................................................................67 3.3.3. Calculation of acoustic indices ................................................................................................68 3.3.3.1. Sound pressure level....................................................................................................69 3.3.3.2. Reverberation time......................................................................................................69

4. MATERIALS AND METHODS............................................................................................ 71 4.1. Description of the enclosures.............................................................................................................71 4.1.1. Reverberant room, Room 1......................................................................................................71 4.1.2. Industrial hall, Room 2a and Room 2b....................................................................................71 4.1.3. Auditorium................................................................................................................................72 4.2. Measuring methods............................................................................................................................73 4.2.1. SWL of the sound source .........................................................................................................73 4.2.2. Measurement of SPL and RT..................................................................................................73 4.2.3. Measurement and estimation of the absorption properties of materials.............................74 4.3. Modeling of the enclosures ................................................................................................................74 4.3.1. Reverberant room, Room 1......................................................................................................74 4.3.2. Industrial hall, Room 2a and Room 2b...................................................................................75 4.3.3. Auditorium, Room 3.................................................................................................................77 4.4. Validation methods ............................................................................................................................79

5. RESULTS ................................................................................................................................. 81 5.1. Reverberant room, Room 1 ...............................................................................................................81 5.2. Industrial hall, Room 2a ....................................................................................................................82 5.3. Industrial hall, Room 2b ....................................................................................................................85 5.4. Auditorium, Room 3...........................................................................................................................88

6. DISCUSSION ........................................................................................................................... 92 6.1. Reverberant room, Room 1 ...............................................................................................................92 6.2. Industrial hall, Room 2a ....................................................................................................................93 6.3. Industrial hall, Room 2b ....................................................................................................................95 6.4. Auditorium, Room 3...........................................................................................................................97 6.5. General discussion..............................................................................................................................99

7. CONCLUSIONS AND FUTURE WORK .......................................................................... 101 8. ACKNOWLEDGMENTS..................................................................................................... 103 9. BIBLIOGRAPHY.................................................................................................................. 104 Appendix 1. Measurement of sound absorption in a reverberation room, ISO 354............................................ A1-1 Measurement of the random-incidence scattering coefficient of surfaces, ISO/DIS 17497-1........ A1-3 Appendix 2. Survey of existing commercial software............................................................................................. A2-1

8

List of symbols and abbreviations α δ

∆

ϑi ϑr Θ λ Σ

ϕ ψi ω A Ai c d dB D f F h I IL k l L LI Lp Lw m M N P P Q r r

r R S S S' t T T V W

absorption coefficient scattering coefficient path difference incidence angle reflection angle temperature wavelength wave front azimuth angle phase shift angular frequency area amplitude speed of sound distance decibel diffuse frequency false height intensity insertion loss wavenumber height length acoustic intensity level acoustic pressure level acoustic power level air absorption coefficient Maekawa function Fresnel number, number of rays source strength pressure, point directivity factor distance, radius vector direction receiver specular sound source virtual sound source time period true volume sound power

2D 3D BEM EDT FEM FDTD FIOH GTD ISO RT SPL SWL TO

two dimensional space three dimensional space Boundary Element Method early decay time Finite Element Method Finite Difference Time Domain Method Finnish Institute of Occupational Health geometrical theory of diffraction International Organization for Standardization reverberation time sound pressure level sound power level transition order

9

Introduction There are many reasons for the development and improvement of computer-aided acoustical modeling of rooms. There are also many reasons to think that no existing method is exact or 100 % accurate 1. This thesis is about the second of these topics, and it presents the author's implementation of the ray-tracing algorithm and an evaluation of its accuracy in the prediction of two acoustical parameters when comparing the measured data of three different studied enclosures. Acoustics is important in many different spaces, especially where noise exposure is high and causes hearing loss or disturbance in the workers, and in places where music and speech performance is very important, i.e. concert halls and auditoria. In the first case we could refer to noise and in the second to sound, but the laws that govern both are the same, and we define it as the problem of room acoustics. The adverse effects of sound on human beings are numerous. For instance, high levels of noise and long exposure to it cause hearing loss in industry workers 2. Noise also affects the workers’ concentration, i.e. in hospitals, schools and in open-plan offices 3. On the other hand, good acoustics are very important in concert halls, theatres and auditoria, but also in public buildings. Probably we have all had the experience of being in an airport or bus station where the announcements coming from loudspeakers are difficult, if not impossible, to understand. All these problems have a technical solution. Noise machines can be timbered and the industrial hall can be acoustically treated 4. The ceiling of open-plan offices can be made more absorptive, and the same solution can be applied in schools and hospitals. In concert halls the solutions are normally more complicated, depending on the subjective opinion of musicians and audience, but normally they consist of achieving relatively long reverberation time and improving diffusion 5. Almost the opposite is normally required in public buildings for message announcements. In general, correct acoustic treatment makes a space more comfortable and desirable. However, the problem in all of them is basically economic, although it is also a matter of common sense. If an acoustic treatment costs money, to do it when the construction is already finished will cost much more. And it is senseless to spend one year designing a building and another building it, to finally obtain a bad environment. Not only

10

acoustics but also lighting, humidity, ventilation, temperature, etc. are key elements that affect occupant satisfaction

6

. Acoustic modeling can give information before

construction about the acoustical problems or lucks that a building may have. However, it seems that communication between architects and acousticians has not been very good 7. In this thesis, I hope to shed some light on and bring some sense to the acoustic phenomena involved in room acoustics. Not only for architects or acousticians, but also for the general public who may find it interesting to know a bit more about it and about the methods that are normally applied in room acoustics modeling. Sound is wave motion within matter 8, and the physical phenomena involved in sound wave propagation in enclosures are both numerous and complex, making overall analytical modeling difficult, if not impossible, to perform 9. Because of the large number of parameters to be taken into account for the description of a real situation, only an approximation of it is possible. That is needed because the equation that describes the propagation of sound, the wave equation, is not solvable in real enclosures 10

. So the problem is how to find a good approximation, which simplifies the problem

but not too much, as it seems, from the existing literature, that the greater the simplification is, the less accurate the results become. The use of a digital computer in room acoustics has opened possibilities for the prediction of acoustical behavior of enclosures

11

, so computer modeling of acoustic enclosures has developed in various

forms, although none of them has yet been demonstrated to have 100 % accuracy 1, 12. The thesis is organized in four parts. In the first part, an introduction to the general principles of acoustics is given. It is difficult to set the limits of this part, and some readers may find no need to read it, as it might not tell them anything new. However, it has been considered a fundamental part in this work, because many approaches and their failure or success could not be understood without the knowledge of the wave propagation principles in enclosures. How sound propagates and distributes, what kind of movement it is, and the main effects due to the wave nature of sound, i.e. diffusion, diffraction and interference, are explained. Moreover, the concept of frequency and wavelength, and the description of the impulse response are given space in this chapter. The second part of this work deals with existing modeling methods. Their principles, limitations and accuracy will be the main topics. Existing methods can be separated into three groups or categories: empirical, wave-based, and geometrical acoustics methods.

11

The first are simple but highly inaccurate and are not able to represent many of the effects that are present in wave propagation in complex rooms. The next type, the wavebased methods, are undoubtedly the most exact, but are very slow and impractical for most frequencies. The last group of methods studies the propagation of sound by geometrical laws. Their complexity varies, as does their accuracy, but they have been considered to be preferred by most of the acoustic community, and the author of this thesis has implemented one form of them. For this reason, special attention is given to them and to their limitations. To end the chapter, the results of two international round robins about the accuracy of existing models are reported in detail. The next and third part of this thesis offers a discussion about the algorithmic implementation of the ray-tracing method. If the list of proposed and existing models seems to be never ending, the way in which they are implemented is rarely described. The reasons for such secretiveness remain unknown, but even when the description of the algorithm is not given, their authors are in the habit of claiming their excellent accuracy, their ability to overcome previous limitations, and the fast and interactive possibilities of their models. The writer of this thesis has doubt about some contributions, especially if they offer auralization based on databases or without taking diffusion into account. The author of this work has developed a new form of the raytracing method. So starting from zero, we have had to solve many problems in the implementation of this new geometrical acoustics model. One example is the approximation done to the point-in-polygon test (detailed later), which is repeated in the innermost loop even millions of times, and which has been improved, giving more competitive computation times. All the steps to implement the model are described, with some paid attention to possible limitations or drawbacks. The last and fourth part, which includes chapters 4, 5 and 6, offers a comparison between modeling and measured data in three enclosures. These enclosures have been chosen to cover some of the different types of indoor spaces that need to be modeled in room acoustics. The first space is the reverberant chamber of our laboratory, which is no more than a room where the wall absorption is extremely small, and therefore sound remains for long time. The second space is an industrial hall with different kinds of machines and with a special geometry. The last room is a small auditorium with ten audience rows and some acoustic treatment in walls and ceiling. The geometrical and acoustical descriptions of the enclosures, together with the way they have been

12

modeled, are detailed in chapter 4. The results from the modeling and the comparison with measured data are given in chapter 5. Finally, there is a profound discussion of the accuracy, never perfect but satisfactory, and of the causes and reasons that conditioned it.

13

1. Overview of sound 1.1. Sound: Wave propagation Sound is wave motion within matter, be it gaseous, liquid or solid, and it exhibits, in many respects but not in all, behavior similar to other wave motions that we encounter in nature, i.e. water waves and light waves, whose propagation phenomena are easily observed 8. A fundamental property of all wave motions is that the matter itself does not travel with the wave. Rather, the wave represents a disturbance or modification of the matter, and it is this disturbance that propagates in periodic movements away from the source (called wave front). In sound propagation, the air being the medium of interest in room acoustics, the disturbance is an alteration of the atmospheric pressure over and under its mean value, which produces a periodic movement of the air molecules backand-forth along the same direction in which the wave is propagating. These are called longitudinal waves, while water waves are a case of transversal waves. The propagation of sound through the air is illustrated in Figure 1.1. In the upper part the alteration of the atmospheric pressure is illustrated, while in the lower part it is the motion of the air molecules associated with the propagation of sound. If the intensity of the sound grows, the gradient of the pressure does too, and as a result more air molecules are in movement.

Figure 1.1: Sound motion through the air (longitudinal waves) The velocity of the airborne sound does not depend on the size of the wave, as it does with water surface waves, neither does it have an electro-magnetic nature as is the case

14

with light. If the medium is assumed to be homogeneous and at rest the speed of sound is constant with reference to space and time 13. For air, its magnitude is c = (331.4 + 0.6Θ )

[m/s]

(1.1)

Θ being the temperature in centigrade.

A general principle, first established by Fermat, states that every wave propagates from the source to the receiver by way of the fastest path: not the shortest, but the fastest. If the medium is homogeneous through the space of interest, as we consider it to be for the air, then the speed of sound is uniform through it, and then the fastest path is also the shortest. That is, the propagation occurs along a straight line. The constant distance separating two points of equal phase on adjacent waves is called the wavelength, λ [m], see Figure 1.1. The time elapsed for a point to return to the same vibration state is the period, T [s]. Hence, the speed of propagation of the wave, c, is given by the quotient of λ , the distance traveled by the wave between two instants of equal phase, and T, the time required for this journey. c = λ /T

(1.2)

The reciprocal of the period is called the frequency, f, and it indicates how many periods occur in a unit of time, or in other words, how many times the particles of air have moved back-and-forth in a second. In honor of Heinrich Hertz, the frequency unity is called the Hertz [Hz]. c=λ⋅ f

(1.3)

As the speed of sound is constant, long waves have low vibration movement (low frequencies), and short waves have high vibration movement (high frequencies). The division of the frequency range into high, mid, and low frequencies is quite arbitrary, especially their limits. But it is useful as it allows idealizing the wave nature effects encountered in room acoustics. Octave bands are used in acoustics rather than individual frequencies

14

. When two

sounds have frequencies f1 and f 2 , they are separated by an interval f 2 / f1 and they

15

define a frequency band of height ∆f = f 2 − f1 . When the relation, f 2 / f1 = 2 , one frequency is double than the other, the interval is called octave. The acoustic octave bands are centered in the frequencies that are listed in the first row of Table 1.1. Their corresponding wavelengths are listed in the second row. The filters that are used to analyze signals are able to eliminate the components which are outside the passing band. Table 1.1: Center frequencies of acoustic octave bands and their corresponding wavelength. f [Hz]

63

125

250

500

1000

2000

4000

8000

λ [m]

5.302

2.672

1.336

0.668

0.334

0.167

0.084

0.042

Mathematically the sound propagation in the space is described by the wave equation, Eq. (1.4), also known as the Helmholtz equation, where p [Pa] denotes the sound pressure, and x, y, z are the Cartesian coordinates. However, the wave equation can seldom be solved in an analytic form

10

, ant that is the reason for the existence of all

existing and applied approximations. ⎛ ∂2 p ∂2 p ∂2 p ⎞ ∂2 p c 2 ∆p = c 2 ⎜⎜ 2 + 2 + 2 ⎟⎟ = 2 ∂y ∂z ⎠ ∂t ⎝ ∂x

(1.4)

The propagation of sound through the air is explained by the Huygens propagation principle, which establishes that a wave front ∑ progresses in the space as if every one of its points would emit elemental spherical waves, being the position of the wave front after a certain time ∆t the envelope of such elemental waves, see Figure 1.2.

Figure 1.2: Huygens principle. Propagation of a wave front.

16

In this way, any point source, which spreads its influence equally in all directions without a limit to its range, will obey the inverse square law, also called the divergence law. This comes from strictly geometrical considerations, it applies to diverse phenomena encountered in nature, and it is illustrated in Figure 1.3. The intensity I of the influence at any given distance r is the source strength P divided by the area of the sphere.

Figure 1.3: Divergence law for sound propagation.

1.2. Sound: Reflection The Huygens principle is also applied to sound reflection as is illustrated in Figure 1.4. Imagine an incident plane front ∑ that arrives at a surface with an angle of incidence θ i . As before, every point in the surface after being illuminated by the sound wave front will act like a secondary source radiating semicircular waves. The reflected sound wave can be obtained by connecting the points that are in phase with each other. In Figure 1.4, only six points (P1, P2, …, P6) that act like secondary sources have been represented. When the incident wave front arrives at P6, point P5 is already emitting and its wave has spread the same length as d5-6. This is done in the same way for the other points.

Figure 1.4: Huygens principle. Specular reflection of a plane wave from a smooth flat surface.

17

The reflection illustrated in Figure 1.4 obeys Snell’s law, i.e. the angle of reflection is equal to the angle of incidence, θ i = θ r , like the mirror reflection of light. That is considered as a special case of Fermat’s principle, where the reflection path between the source and receiver via the panel is the one that is traversed in the least time. A specular reflection of sound can be obtained approximately from a plane, rigid surface with dimensions much larger than the wavelength of the incident sound 15. However, Snell’s law fails in the reflections of wave fronts from everyday walls, and this is easy to understand when we apply Huygens’s principle to the reflection of sound over rough surfaces. In the same way as before, we now study the reflection of a plane wave by a surface whose irregularities are comparable to the size of the wavelength. The case is schematically illustrated in Figure 1.5, where we observe how the wave front is ‘broken’ in pieces, which is a consequence of the irregularities of the surface, and a result of the different delay times achieved by the different parts of the wave. This kind of reflection is called a diffuse reflection and it occurs when the size of the irregularities and the wavelength are similar. On the other hand, when the wavelength is very long compared with the irregularities, the distance to the next wave front is so great that the disturbance can be considered negligible. In the opposite case, when the wavelength is very short, the distance between points with the same phase is smaller than the size of the irregularities, and because the reconstruction of the wave front is done over points with the same phase, the effect may not be noticeable. For a correct reconstruction of the wave front of Figure 1.5, edge diffraction and interference effects should also be taken into account, but to the author’s knowledge, there is no existing algorithm able to solve it mathematically in an efficient and simple way.

Figure 1.5: Huygens principle. Reflection of a plane wave front from a surface, whose irregularities are comparable to the wavelength of the sound wave.

18

In a diffuse reflection the sound incident to a surface is reflected into a wider solid angle than that of a specular reflection

12

, and it is said that the sound has been scattered

16

.

Scattering of sound inevitably occurs when sound is reflected from a surface, not only by the fact that real life surfaces are not perfectly smooth nor uniform, but also because of the effects created by their finite size. This last case corresponds to diffraction effects, which are explained more deeply in the next section. The effect of a diffuse reflection is simplified in Figure 1.6, where the arrows illustrate the direction of the front waves. The ideal reflection follows Snell’s law, i.e. θ i ≡ θ rS , and a diffuse reflection is such that θ i ≠ θ rD . A larger problem lies in the difficulty of determining the direction of the sound wave after such a reflection, considering the possibility that there can be more than one.

Figure 1.6: Reflection of sound by a rough surface.

1.3. Sound: wave effects Effects due to the wave nature of sound are described in the following. These effects depend directly on the wavelength of sound, i.e. diffraction and interference being more pronounced at low frequencies, while scattering due to the roughness of surfaces is more pronounced at high frequencies. This last case has just been discussed.

1.3.1. Diffraction Diffraction is the local change in the direction of the propagation of sound waves passing the edge of an obstacle. This is one of the most important acoustic phenomena caused by the wave nature of sound 17, and also one of the most complex to solve. The diffraction phenomenon depends significantly on the ratio of the sound wavelength to the size of the obstacle. The longer the wavelength, the stronger the sound diffracts. Diffraction theory has been extensively studied in optics. But, being an effect that

19

depends basically on the wavelength of the wave, its application and representation in acoustics are quite difficult. When a sound wave encounters an obstacle, which is small in relation to its wavelength, the wave passes round it almost as if it did not exist, forming very little shadow. But, if the frequency of the sound is sufficiently high and the wavelength is therefore sufficiently short, a noticeable shadow is formed

18

. Due to the consequences for

acoustics and room acoustics, diffraction effects can be divided into three groups; by a barrier, by an edge, and by an aperture, and are illustrated in Figure 1.7, where the ideal cases for low and high frequencies are represented. For instance, the most obvious effect of diffraction by a barrier is the presence of sound behind it. On the other hand, diffraction by an edge or a corner causes diffusion or scattering of sound. Diffraction by an aperture is a general cause of energy losses, and it is not treated in greater detail in this study.

Figure 1.7: Diffraction by a barrier, by an aperture and by an edge. Upper cases for ideal low frequencies. Lower cases for ideal high frequencies.

1.3.2. Interference In room acoustics there are many important interference phenomena, and these can be only simulated if reflections of pressure waves with amplitude and phase are considered 19

. Although, as we saw in section 1.1, sound propagation corresponds to the motion of

20

longitudinal waves, we can visualize interference by the principle of superposition applied to sinusoidal waves (however, the atmospheric pressure gradient corresponds to such movement). The principle of superposition may be applied to waves whenever two or more are traveling through the same medium at the same time. The waves pass through each other without being disturbed. The net displacement of the medium at any point in space or time is simply the sum of the individual wave displacements. For instance, if two waves with the same amplitude ym and wavelength λ = 1 / k are traveling in the same direction along the same line x, by applying the principle of superposition the resulting medium displacement may be written as:

φ⎞ ⎛φ ⎞ ⎛ y ( x, t ) = ym sin(kx − ωt ) + ym sin(kx − ωt + φ ) = 2 ym cos⎜ ⎟ sin ⎜ kx − ωt + ⎟ 2⎠ ⎝2⎠ ⎝ (1.11) which is a traveling wave whose amplitude depends of the phase φ . When two waves are in-phase, φ = 0 , they interfere constructively and the result has twice the amplitude of the individual waves. If the two waves have opposite-phase, φ = π , they interfere destructively and cancel each other out. However, if two waves are traveling in opposite directions they create a standing wave, which appears to be still vibrating in place. So if two waves of the same kind are traveling in opposite directions, the principle of superposition will lead to: y ( x, t ) = ym sin( kx − ωt ) + ym sin( kx + ωt ) = 2 ym sin kx cos ωt

(1.12)

This wave is no longer a traveling wave because the position and time dependence have been separated. The amplitude of the wave, 2 y m sin kx , does not change position (absence of propagation), but stands still and oscillates up and down according to cos ωt . Characteristic of standing waves are locations with maximum displacement (antinodes) and locations with zero displacement or nearly complete cancellation (nodes). Nodes and antinodes, if present, have a fixed distribution in the space. Standing waves are commonly set up in rooms by low frequency sound with a wavelength being a multiple of half the length of one side, Eq. (1.13).

21

λn =

2L n

n=1, 2, 3…

(1.13)

To ensure that interference effects are stable in time the room walls need to be highly reflective. Standing waves contribute to higher ambient noise levels, and therefore are to be avoided in acoustic design, particularly in industry. They should also be avoided in concert halls to ensure uniform frequency response and good diffusion. In typical industrial spaces, interference effects are reduced by the fittings in the space and the large number of sound sources which are likely to be non-coherent 20. However, interference is not only observable between parallel walls, as there are also interference effects at reflecting boundaries of reverberant sound fields

19, 22

. Even in a

reverberant or diffuse field, where at all points there is equal mean energy flows in all directions, Waterhouse

21

stated that the sound energy is distributed into interference

patterns at the reflecting boundaries. Thus, the mean energy density is not uniform at all points in the field. The phases of the wave trains near the reflecting surfaces are not entirely random, and this lack of phase-randomness gives rise to the interference patterns. More simply explained, the incident sound interferes with the reflected sound. In their studies, Waterhouse

21

and Waterhouse and Cook

23

found the energy level at

the surface to be 2.2 dB higher than at points further away where the interference patterns are negligible, with the largest departure from uniformity at the corners, with levels 9 dB higher that at remote points. In Figure 1.8a we observe a wave front arriving at a corner and the area where both fronts, incident and reflected, will interact. In Figure 1.8b we observe the interference pattern at a reflecting boundary when plane sound waves are randomly incident over the positive quadrant. The numbers on the contours give the theoretical departure in dB from the steady level, as a function of distance in wavelengths ( z / λ and y / λ ) from two rigid reflecting planes xy and xz, intersecting at right angles, both planes being perpendicular to the plane of the paper.

22

Figure 1.8: a) Incident wave front to a corner. b) Interference pattern at the corner. After Waterhouse and Cook 22.

1.4. Temporal distributions of reflections; the impulse response The impulse response is the time domain response of an enclosure to an idealized infinitely short impulse (Dirac function). In room acoustics it is also understood as the temporal distribution of reflections. The sound propagation from a sound source to a receiver is characterized by the impulse response, and the spatial information on all possible propagation paths between the two positions is reflected in this signal. The first detection always corresponds with the direct sound, and after it, first and multiple reflections are received. Their time delays with respect to the direct sound are according to the lengths of the traveled paths, and their pressure intensities depend on the air absorption of sound, and also on the absorption characteristics of the surfaces involved in every path. Both ways of sound attenuation are explained in section 1.7. For every frequency there is a different impulse response, although it is possible to obtain them in octave bands. Every detection is also characterized by its phase, and if that is neglected, then we are talking about an energy response. In Figure 1.9, an idealized energy response is presented, which has been separated into three parts: direct sound, early or first reflections, and late reverberation. The selection of the limit between early and late reflections is quite arbitrary. In speech it is normally 50 ms, while in music it is 100 ms.

23

Figure 1.9: An idealized energy response. The measured and/or modeled impulse responses give a lot of information about the acoustic characteristics of the enclosure. From these, it is possible to calculate many acoustic parameters, i.e. reverberation time RT and early decay time EDT

24

. The

calculation of these parameters is based on the early-to-late energy ratio of the impulse response. It may be desirable, for instance, and especially in acoustics design, to listen to a virtual environment before it has been constructed. To do that, the auralization process has to be complete, and the impulse response is a crucial part of it. The process used to obtain an auralized signal is not described in this thesis, but the reader can obtain basic information about it from references 45 and 48.

1.5. Diffuse sound field A sound field in an enclosure is composed of many single waves. Each of them has its own particular amplitude, phase, and direction. However, when different conditions are fulfilled it is possible to ignore these, and then we are talking about the existence in the enclosure of a diffuse sound field. Waterhouse stated in 195521 that a diffuse sound field is characterized by the following three properties: 1. Uniform distribution of sound energy, i.e. equal mean energy density at all points. 2. Equal mean energy flow in all directions at all points in the field. 3. Random phase relation between the wave trains converging on any point. Forty years later, Hodgson arrived at the same conclusion

25

, and resumed that in a diffuse sound

field, at any position, sound waves are incident from all directions, with equal intensities

24

and with random phase relations 26. Perfect diffusion is achieved in a box of volume d3, when d << λ . From another point of view, Beranek 27 considered that two requirements need to be met for diffusion of the sound to be good: First, the reverberation time must be fairly long, because the sound will have died out after relatively few reflections if the room is not reverberant. Since sound to be diffuse, must undergo many reflections in a room, a long reverberation time contributes to diffusion. Second, the ceiling and walls of the hall must be irregular so that the sound waves are scattered when they are reflected from them. The role of sound diffusion in the good acoustics of an auditorium has been difficult to prove, as diffusion has been difficult to measure, and the role of diffusing elements or room surfaces difficult to quantify 28. It happens that while some acoustic engineers can blame the disappointing acoustics of certain major halls on a lack of surface diffusion, one eminent concert hall designer regularly claims that too much diffusion is detrimental to the sound quality of the upper strings 28. On the other hand, the existence of diffuse sound field in an enclosure allows the application of simplified prediction models, and also of some important simplifications, e.g. negation of interference.

1.6. Acoustic levels In acoustics, logarithmic scales are used to express sound levels, as the human audible range for the intensity varies between 10-12 and 10 watts per square meter, and as the human hearing has not a linear but a logarithmic response to the change in sonic perturbation 13. The unit used is the decibel, dB. A level is the magnitude of a signal in reference to an arbitrary reference. In other words, a level is the logarithm of the relation between a certain magnitude and a reference one of the same class. It is needed, to express the base of the logarithm, the reference magnitude and the kind of level. One thing worth mentioning about the utility of acoustic levels, is that changes of 5 dB give the same change in the hearing feeling whatever the reference level, but a change of 0.01 Pa is a big change at low levels, while not perceptible at high ones.

25

The acoustic power level Lw of a sound source is: ⎛W ⎞ LW = 10 log⎜⎜ ⎟⎟ ⎝ W0 ⎠

where W0 = 1 pW

(1.14)

where I0 = 1 pW/m2

(1.15)

The acoustic intensity level LI is: ⎛ I LI = 10 log⎜⎜ ⎝ I0

⎞ ⎟⎟ ⎠

The acoustic pressure level LP is: ⎛ p⎞ L p = 20 log⎜⎜ ⎟⎟ ⎝ p0 ⎠

(1.16)

where p0 = 20 µPa , and the factor 2 in the exponent of the logarithm is obtained as the sonic power is proportional to the square of the pressure. The combination of two levels, i.e. LI1 and LI2, can be done graphically using a diagram to add levels in dB 13, or analytically, by solving from Eq. (1.15) the intensities I1 and I2. The total intensity is the sum of both, and from it we obtain the corresponding level. That is: I1 = I 0 ⋅10 0.1L I 1 and I 2 = I 0 ⋅ 10 0.1L I 2

(

I T = I1 + I 2 = I 0 10 0.1L I 1 + 10 0.1L I 2

(

(1.17)

)

(1.18)

)

L IT = 10 log 10 0.1LI 1 + 10 0.1LI 2 = 10 log

IT I0

(1.19)

1.7. Sound absorption There are two mechanisms of sound absorption in room acoustics. The first is a consequence of atmospheric absorption effects absorbent properties of surfaces.

29

. The second is a consequence of the

26

1.7.1. Air absorption The attenuation coefficient of the air, m, depends on the ambient atmospheric temperature and pressure, molar concentration of water vapor, and frequency. While a pure-tone propagates through the air over a distance L, its power W decreases exponentially according to Eq. (1.20), where Wini is the initial power. W ( L) = Wini e − mL

(1.20)

The formula that allows the calculation of m is not considered in this work, but it can be found in reference 100. In Figure 1.10, the values of m for a temperature of 20 0C and a relative humidity of 60 % are given. Although we have illustrated the values of m in decibels per meter, we can see how the absorption of sound increases as the frequency does, and thus the exponent in Eq. (1.20) tends to zero. Atmospheric Attenuation of Sound in dB/m (Temperature 20°C, Humidity 60 %) 0.10 0.09

Attenuation, dB/m

0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 1

10

100

1000

10000 Frequency, Hz

Figure 1.10: Atmospheric attenuation of sound in dB/m

1.7.2. Wall absorption The second mechanism of sound absorption in room acoustics is caused by the absorption characteristics of the surfaces. In this case, absorption is defined as the dissipation of sound energy on striking a physical surface. At every reflection, a portion

α of its energy or power is absorbed as described in Eq. (1.21). This factor α is called the absorption coefficient. The portion of sound energy not absorbed is either bounced back as reflected sound, or transmitted through the material. However, when studying the acoustics of an enclosure, the factor which is transmitted can be considered as if being absorbed.

27

W = Wini (1 − α )

As described by Yuzawa

(1.21) 30

the sound absorption depends upon the angle of

incidence, α = α (θ ) . Figure 1.11 illustrates the results for feather glass. It is important to notice how the values for this material remain more or less stable until angles approach grazing incidence, and the extrapolation which is done between 70 and 90 degrees. At grazing incidence the value is zero.

Figure 1.11: Oblique incidence sound absorption of feather glass. After Yuzawa 30. The standardized method to measure the absorption coefficients of materials is described in Appendix 1.

28

2. Modeling methods in room acoustics Different modeling methods applied in room acoustics are discussed in this chapter. Because the final algorithm used is based on the ray-tracing method, special attention is given to geometrical acoustics methods. In Appendix 2, a review of different existing commercial software for room acoustical modeling is given.

2.1. Modeling methods The physical phenomena involved in sound wave propagation inside enclosures are both numerous and complex, making overall analytical modeling difficult. Because of the large number of parameters to be taken into account in the description of a real situation, only an approximation of it is possible 9. There are several approximations, some better than others, but still no one has been proved to efficiently solve both the high and the low frequency problem. From the computational point of view, there are three different approaches; Empirical methods, wave-based methods, and geometrical methods

10

. All of them are overviewed in this chapter, although only the last one is

treated in detail. In 1968, Krokstad et al

47

carried out the first application of ray-tracing computer

modeling. They studied the space and time distribution of the early reflections over the audience area of a concert hall, and also stated that the use of the digital computer opens up possibilities for the prediction of the acoustical behavior of concert halls. Since then, many different approaches and implementations have been carried out, but although computer modeling of acoustic enclosures has developed into various forms, none has yet demonstrated 100 % accuracy 1, 12. The problem lies in the wave nature of sound and in the complexity of its propagation and reflection in enclosures. Eq. (1.4) has a different solution for each frequency.

2.1.1. Empirical and statistical methods Empirical methods belong to the first group of approximations, the Sabine and Eyring models being the most popular ones. Their methods offer a rough estimation in octave bands of parameters like sound pressure level, SPL [dB], Eq. (2.1), and reverberation time, RT [s], Eq. (2.2). They consider the existence of a diffuse sound field in the space. In that way, only the size of the room and the importance of absorption, but not of the location of each of the surfaces is considered 21. The Sabine and Eyring theories are still

29

used to predict noise levels in factories 32. Diffuse field theory is used by practitioners to predict sound fields in rooms of every type, but what is often forgotten is the fact that the theory is based on assumptions which limit its applicability. If the theoretical assumptions do not hold in the case of a particular room, the predictions may not be accurate

26

. For instance, diffuse field theory can not be applied in highly absorptive

rooms. Other methods and refinements have been proposed, e.g. by Hodgon 33, Kuttruff 34

, Thompson et al.35, Cotana

36

and others

37, 38

. However, those methods have their

own approximations and limitations, and have not proved to be efficient over a variety of enclosures. Nor are empirical methods suitable for auralization. 1−α ⎤ ⎡ 1 SPL (r ) = L w + 10 log ⎢ +4 2 A ⎥⎦ ⎣ 4πr

RT = 0.16

with

(2.1)

V A

(2.2)

ASabine = 4mV + ∑α i Si

(2.3a)

AEyring = 4mV − ∑ Si ln (1 − α i )

(2.3b)

i

and

i

where Lw [dB] is the source power level (re. 1 pW), r [m] the distance from the source, m [m-1] the air absorption coefficient, V [m3] the volume of the enclosure, α i the absorption coefficient of the surface i, and Si [m2] its corresponding area.

2.1.2. Wave-based methods Wave-based methods can not solve the wave equation, but they try to numerically approximate its solution. These are, FEM (finite element method)

39

, BEM (boundary

element method) 40, and FDTD (finite-difference time-domain method). Basically, they divide the space into small elements or nodes. These elements interact with each other according to basics of wave movement phenomena, for instance, as linear combinations of a finite number of parameters. The size of these elements has to be much smaller than the size of the wavelength (at least six nodes) for every particular frequency, and that is a problem at high frequencies. The number of natural modes in a room increases approximately with the third power of the frequency, which means that for practical use,

30

wave models are typically restricted to low frequencies and small spaces 1. Wave models are also characterized by their creation of very accurate results at single frequencies. The most difficult part in their application is the definition of the boundary conditions, as a complex impedance of the surfaces is required, and these are sometimes very difficult data to find from the existing literature 10. FEM is primarily used when only the interior acoustic field of an enclosure is to be computed. In FEM, the complete space is represented by elements, creating matrices that are large but sparsely filled. Each of these elements only interacts with the elements directly adjacent to it. On the other hand, BEM can be used to compute the interior, exterior, or both fields simultaneously. In BEM, only the boundaries of the space are discretized, but all nodes are linked to each other, forming a dense coefficient matrix. The main principle in the FDTD technique is that derivatives in the wave equation are replaced by corresponding finite differences. Those methods are only suitable for small enclosures and for low frequencies due to their extremely long computation times. They are mostly applied in studies of the automotive industry where control of vibrational sound is the main interest, although they have also been applied, e.g. in studying the performance of noise barriers 41.

2.1.3. Geometrical acoustics methods The next group of methods follows high frequency approximations based on geometrical propagation paths 42. In principle, they consider the propagation of sound through the air in straight lines, avoid the wave nature of sound, and model in one or another way reflections from boundaries. There are many of them

43, 16

. They can trace

rays, particles, cones, even pyramids, and model in one or another way diffraction and/or diffusion. However, none has proved to be exact

12

and some of them are not

even good enough. Those methods are also called energy-based methods, because in the end, the only thing they consider is the energy which was assigned to every ray or portion of wave (instead of amplitudes), and after all, they are also a statistical approach. However, geometrical acoustics methods generally have an acceptable degree of accuracy in the prediction of different acoustic parameters. They are also able to achieve the temporal and spatial distribution of reflections, and so they are suitable for auralisation

44, 45, 46

. It has been said that geometrical room acoustics can reflect only a

partial aspect of the acoustical phenomena involved in a room, but that, however, is an

31

aspect of great importance 13. In the following we will look at the principles of every method, also at some of their drawbacks. Special interest will be paid to the ray-tracing method, although many parts in the implementation are also common to the other geometrical methods. As mentioned in Kuttruff’s work

13

, to apply wave theory in the resolution of sound

propagation is more than a difficult task. Instead, we can apply a greatly simplified way of description -just as in geometrical optics- by employing the limiting case of vanishingly small wavelengths, i.e. the limiting case of very high frequencies. Therefore, the wave nature of sound is neglected, and the propagation of sound is studied like the propagation of sound rays (portions of spherical waves with vanishing aperture). This assumption is valid when the wavelength of sound is small compared to the area of the surfaces of the room, and large compared with their roughness. This is frequently met in room acoustics. At a medium frequency of 1000 Hz corresponding to a wavelength of 34 cm, the linear dimensions of the walls and ceiling, and also the distances covered by the sound waves are usually larger that the wavelength by orders of magnitude. All phenomena due to the wave nature, such as diffraction and interference, are ignored

10

, since propagation in straight lines is its main postulate 13,

i.e., when several sound field components are superimposed their mutual phase relations are not taken into account; instead, their energy densities or their intensities are simply added. 2.1.3.1. Ray-tracing method The concept of ray tracing is very simple and is represented in 2D in Figure 2.1. The sound power emitted by a sound source is described by a finite number of rays, which will be considered as carriers of power (or of energy or intensity). These rays travel through space at the speed of sound, and are reflected after every collision with the room boundaries. During that time, their energy decreases as a consequence of the sound absorption of the air and of the walls involved in the propagation path. When the rays cross predefined volumes, the receivers, an energy calculation process is performed, and those data stored. Finally the impulse response at every receiver is obtained, and with it all desirable acoustic parameters.

32

Figure 2.1: 2D representation of the ray-tracing algorithm. If the sound source has omni-directional characteristics, the directions of the virtual sound rays are created from homogeneous random distributions, although special directivity characteristics of sound sources can also be accounted for. The possible path for every new ray emitted from the sound source needs to be found. Considering the ray like a vector that changes its direction upon every reflection, it is possible to determine which and only which surface will be collided with next and where. This is called the point-in-polygon test, see section 3.2.2, and it is a process that is repeated many more times. A wise implementation of this loop is recommended when complex room shapes also involve a large amount of surfaces. After this the receivers are checked to see whether they are crossed by the ray. The energy that the ray ‘impinges’ to a receiver depends directly on the length that it has traveled inside it, and this energy is stored in the time domain with its corresponding time delay. This is how the correspondent energy response is obtained at that point. The attenuation of the energy is calculated according to the absorption of the air and the length of the traveled path, and according to the wall absorption properties. The collision with the surface can be considered specular, following Snell’s law, or diffuse, following a direction other than specular. Different ways to model diffusion are detailed in section 2.2.1. The sound ray now has a new direction and the next surface to be collided with has to be found. Here the loop goes back to the starting point and the process is repeated until the energy of every ray has fallen under a certain limit. The selection of this limit is arbitrary depending on which parameters are being sough. In the case of calculation of RT or detailed echograms, the limit should be one millionth of the initial energy (a decay of 60 dB). If only SPL is desired, the limit can be reduced to 20 or 30 dB. The overall process has to be repeated again for every frequency band, unless the directional distribution of the diffuse reflected sound is considered constant in the frequency range. In any case, the air and wall absorption properties are frequency-dependent.

33

We can see already that the divergence law (1/r2 law) is implicit in the ray-tracing algorithm. It is taken into account by the constant area of the receiver, which comprises a solid angle decreasing proportionally to 1/r2 as the distance to the source increases 47. Although a description of the ray-tracing algorithm has been given, its detailed implementation is described in Chapter 3. There, for instance, it is explained how to detect the collided surfaces (the point-in-polygon test) and the receivers (detection problem), how the sound is absorbed, or how the energy calculation is performed. Ray tracing has been shown to predict noise levels with very good accuracy, although the accuracy is always low at low frequencies

32, 48, 49, 50

. That is logical since the wave

nature of sound is neglected, that nature being the reason for important effects like diffraction or interference. 2.1.3.2. Image method To model an ideal impulse response all the possible sound reflection paths should be discovered

10

and this is guaranteed by the image method 42. It can determine to a high

degree of accuracy in terms of level, arrival time and direction, the early reflections 51, which are the most important in the subjective hearing feeling of listeners 52. However, they could also create high order specular reflections which will be much more accurate than would be possible with a real sound wave 51. A disadvantage of this method is that it requires very long calculation times for higher order reflections in room geometries that are not simple rectangular box shapes. Increasing reflection order leads to an exponential increase in image sources, and therefore, an exponential increase of computation time, which makes the method ineffective 10, 53. In the image method, reflected paths from the real source are replaced by direct paths from reflected mirror images of the source. This process is quite simple and it is illustrated in Figure 2.2 (subscripts used in this figure are f: floor, c: ceiling, w: wall). The sound source, S, is mirrored to the other side of every surface, so virtual sources are created, e.g. through the ceiling S’c1 and S’c2, and through the floor S’f. The reflected path is then replaced by a virtual path that goes from the virtual source to the receiver location R. Then we need to find out if the path really exists, and that is done by the visibility test, which has two parts. The first part checks whether the collision point of

34

the virtual path with the surface is inside the surface or not, in a similar way as is done in the ray-tracing method. For instance, in Figure 2.2, we observe how the virtual source S’c1 can not directly see the receiver while S’c2 can. The second part of the visibility test checks if there are any obstacles that could obstruct the path. This is done by a backwards reconstruction of the propagation path, starting from the receiver location and finishing in the sound source, via all the walls that have contributed in the reflections. In this way, it would be found that the virtual source S’f can not contribute to the impulse response obtained at R. Multiple reflections are found in the same way by mirroring virtual images through other surfaces, i.e. the path ‘source-ceiling2-backwallreceiver’ is obtained when mirroring the virtual source S’c2 through the back wall, from where S’c2-w is obtained, and then tracing all the paths backwards. Descendants of invisible sound sources must still be considered, as high-order reflections may generate visible virtual sources. It is important to notice, finally, that the locations of the image sources are not dependent on the listener's position, and only the visibility of each image source may change when the listener moves. In the same way as with ray-tracing, and due to the computational requirements for finding high order images, the implementation of the image method needs to be as wise as possible. There are simplified algorithms for simple-shaped spaces extension to arbitrary polyhedra was described by Borish 54.

Figure 2.2: 2D representation of the image method.

13

, and the

35

Although it has been stated that the image method is more efficient than ray tracing 55, this is not true and it is negated by several experimental studies 16, 49, 50, 56, 57. 2.1.3.3. Particles The particle-tracing method is a simplification of the ray-tracing method. The procedure is the same, but particles are traced instead of tracing rays. These also travel in straight lines and are reflected by the room boundaries. The difference lies in the way sound is attenuated. Instead of decreasing the energy of the particles at reflections, the particle disappears or continues. To do this, a random number between 0 and 1 is simply created and compared with the absorption coefficient of the corresponding surface. If the number is smaller than the absorption coefficient the particle disappears, otherwise it continues. The point of doing this could be to speed up the algorithm and save computational time. But this method requires a very large number of particles to achieve realistic results, and those come close to the results achieved by ray tracing when the number of particles tends to infinite. However, the simplification is far from representing the reality, and still has all the limitations common to ray-tracing models. 2.1.3.4. Fitting-zone method The fitting-zone method is another feature of the ray-tracing algorithm. Instead of carefully modeling the interior of a space, the space is subdivided into a number of subvolumes, and to each of them is assigned a fitting scattering cross-section density and a fitting absorption coefficient

32

. It is possible to consider the fitting density

homogeneous in the space (isotropic zoning), or to concentrate it close to the floor surface or in certain areas (anisotropic zoning)58. Then, two kinds of approaches can be used when the ray tracing is performed. The first possibility is to follow each of the rays that emanate from the source until they strike a surface or obstacle. The ray is then redirected according to the appropriate reflection law (specular reflection in the case of a surface, diffuse reflection follows random or Poisson distribution in the case of an obstacle). The other possibility is to compute random numbers every time a ray crosses a sub-volume and compare them with the corresponding fitting density. When the random number is smaller that the density the ray is randomly redirected, otherwise it is just left to go through. In both possibilities the power of the ray decreases according to the traveled path (air and fitting absorption).

36

In 1989, Hodgson

32

carried out a study on factory noise prediction with this approach,

comparing predicted and measured sound pressure levels. Although he found differences at individual points in the range from –7 to +6 dB, with low accuracy at lower frequencies, and with an equal tendency to overestimate and underestimate, he concludes that the model is shown to give excellent prediction accuracy. A similar study was published in the same year by Ondet and Barbry 9, who studied an application of this method while studying sound propagation in fitted workshops. Although they demonstrated the influence of fittings in noise levels and affirmed that predicted results agreed closely with measured ones, the results were only obtained at 2 kHz and did not prove the method to be efficient. Dance

58, 59, 60

also studied the influence and

representation of fittings in a computer prediction of sound propagation in enclosed spaces. He considered that computer models for the prediction of sound fields in rooms were complex and required a skilled acoustician to use them effectively, and hence there was a need for simpler models. Considering that, typically, a workroom is rectangular, with uniform wall absorption and all fittings located in the floor, Dance suggested models that require a minimum amount of input data, reducing the scope of the model to only two or three acoustic parameters. From the data shown in his studies it appears that the more detailed the description of the room, the more accurate and consistent results were obtained. This implies that a large simplification of the environment leads to inaccurate results. 2.1.3.5. Cone or pyramids tracing method It is also possible to trace cones, approximate cones, pyramids, etc, in practice portions of solid angles, with a sphere coverage that spreads from the sound source with divergence. In this way, computation time can be reduced, although the implementation can be very difficult. The idea is illustrated in 2D in Figure 2.3. A sound source S and two receivers, R1 and R2, are inside a five-side polygon abcde . The vertices of the cones are equivalent to mirror sources (section 2.1.3.2) and mirrored in the familiar way, i.e. Sa. The cones are hence, mirror image sources with built-in visibility borders 61. In Figure 2.3, we can see how R1 is outside the illuminated region of Sa-d, while R2 is inside it. We also see how the cone of Sa is clipped by side d.

37

Figure 2.3: 2D representation of cone-tracing method in a five-side polygon. S is mirrored through side ‘a’, giving Sa. Sa is mirrored through side ‘d’, giving Sa-d. Sa-d illuminates R2 but not R1. One of the first applications of this method was carried out by the commercial software Ramsete 62, where triangular beams were generated at the sound source. The central axis of each pyramid was traced as in ray tracing, being specular-reflected when a surface was hit. The three corners follow the axis, all being reflected at the same plane. It considers point receivers, and detection occurs when the receiver is inside the pyramid being traced. 2.1.3.6. Hybrids There is the possibility of combining different features of different methods, i.e. hybrids. For instance, the image method could be applied to calculate very accurately the first reflections, and then switch to ray tracing for the later part of the impulse response. Another possibility

63

includes the natural combination of a rigorous

numerical method for low frequencies, i.e. BEM, and a more efficient method for higher frequencies, i.e. the ray-tracing method.

38

2.2. Diffusion and diffraction in the ray-tracing algorithm This thesis centers its attention now on the ray-tracing algorithm. In section 2.1.3.1 its basic concept was introduced, i.e. how sound rays are shot from the sound source and followed through the enclosure. Before we can describe the way to implement the algorithm, two aspects have to be discussed, and they are the different ways to model diffusion and diffraction of sound.

2.2.1. Modeling diffusion in the ray-tracing method It is presently admitted that accounting for the scattering or diffusion of sound by surfaces improves room acoustics predictions Kuttruff

13

, Hodgson

56, 57

, Lam

16

, and others

47

. The investigations carried out by

49, 50

, have shown that room acoustic

calculations which are based only on specular reflections overestimate reverberation times in situations where diffusion is significant. For instance, Howarth and Lam

12

studied, in seven enclosures, the relation of the prediction of acoustic parameters with the way diffusion was modeled. They found that a sufficient portion of reflected energy needed to be modeled as diffuse for predictive accuracy, and that diffusion coefficients had to be defined in the frequency domain. They stated: Accounting for scattering in the modeling of all reflections orders (since the first reflection) would improve accuracy in the time domain and assist in the prediction of acoustic parameters such as reverberation time. Allowing the directivity of scattered energy to be modeled would improve spatial accuracy. 2.2.1.1. Statistical treatment of diffuse reflections. The scattering coefficient δ Diffuse reflections, if treated at all, are usually treated in a crude manner 43. The major problem lies in the difficulty of knowing the angle diffusion properties of a surface (directivity pattern of the scattered sound). Quantification of these would require very precise laboratory measurements, and the fact that every frequency has a different diffuse behavior for every angle of incidence, would make the amount of data very difficult to handle. Therefore, simulating the scattering effects of diffuse surfaces by simplified methods has been accepted. A first approximation in ray-tracing models is to assume four kind of reflections combinations 43. Those are; S-S, S-D, D-S, D-D, where S represents ‘specular’ and D ‘diffuse’, see Figure 2.4. Due to the wavelength of sound, perfect specular reflections might not exist, and in this way, the approximation D-S may fail.

39

Figure 2.4: Four reflection combinations possible when both specular and diffuse reflections are taken into account by the ray-tracing method . S-specular, D-diffuse. After Dalenbäck43. The use of the scattering coefficient is the most popular method to weight the relation between specular and diffuse reflections in geometrical acoustics. This coefficient δ is defined as the ratio between reflected sound power in non-specular directions and the total reflected sound power 16, 51, Eq. (2.4) and Figure 2.5.

δ=

Wref . Non − specular Wref . Non − specular + Wref .specular

(2.4)

Figure 2.5: Energy components of reflected sound. Therefore, to every surface two octave band dependent factors can be assigned, an absorption coefficient α , and a scattering coefficient δ . The relationship of these factors is illustrated graphically in Figure 2.5, and it is described by Eq. (2.5), being on it the terms, absorbed, scattered, and specular.

α + δ (1 − α ) + (1 − δ )(1 − α ) = 1

(2.5)

The selection of δ is made according to the geometrical properties of every surface, and somehow, it is just a statistical way to define diffusion. When a ray encounters a diffusing surface, a random number in the range [0,1] is generated. If the number is smaller than δ the ray direction is randomized to simulate a diffuse reflection,

40

otherwise the reflection is specular

13, 62, 64

. In principle, the stochastic nature of the

algorithm needs one pass per frequency band, since the diffusion factor of a surface is frequency dependent

62

. However, today’s knowledge about which scattering

coefficients are realistic and how to determine them is very limited, and some authors consider it sufficient to characterize each surface by only one scattering coefficient, constant for all frequencies

1, 15, 16, 51, 65

. The author of this thesis thinks that this

assumption can perhaps be accepted if both causes of diffusion are of similar importance, as is illustrated in Figure 2.6. A positive point of this approximation is that the computation time will improve when only one computation pass is needed.

Figure 2.6: Scattering coefficient as a function of two factors. Roughness and finite size of surfaces. Keränen et al. 4, 66 carried out a study of noise control in six industrial spaces. Based on the idea of a single scattering coefficient for all frequency bands, they applied the following approximation for the selection of the scattering coefficient, Table 2.1, which combines both finite size and roughness of surfaces. Other studies have suggested the advantage of weighting diffusion coefficients at low frequencies according to surface dimensions 12, 67. Table 2.1: Suggestion for choosing scattering coefficient accounting for surface size and roughness. After Keränen et al. 4. Scattering coefficient, δ

Description of the surface

0.1,…, 0.19

Large, plain surfaces

0.2,…, 0.39

Large partially fitted surfaces

0.4,…, 0.59

Small or fitted surfaces

0.6,…, 0.89

Large densely fitted surfaces

0.9,…, 1.00

Small densely fitted surfaces

41

Another aspect of the scattering coefficient δ is that it not only depends on the frequency, but also on the angle of sound incidence. However, analogue to the randomincidence absorption coefficient obtained in reverberant rooms (see Appendix 1), and angular average of the scattering coefficient, i.e. the random-incidence scattering coefficient, can be defined as well. Vorländer and Mommertz proposed an experimental method to obtain δ

65, 68

, which has been accepted, adopted and standardized by the

International Organization for Standardization in ISO/DIS 17497-1

15

. The first

published application of this new method has already been carried out by Jeon et al 69. The experimental method to determine δ is also explained later in Appendix 1. Different methods for measuring the scattering coefficient of a surface have been proposed by Farina 70, and by Haan and Kwon 71. Another approach was given in the modified scattered diffuse energy model (SU model)16. Upon each reflection, the diffuse-reflection coefficient δ is used to define the fraction of energy diffusely scattered into non-specular angles, while the remaining energy is reflected specularly in the usual way. The diffuse energy is then assumed to decay exponentially and contribution to the receivers is calculated. This method has an important drawback, which is that the estimation of the decay constant for the diffuse energy is estimated from Sabine or Eyring’s formula since the energy is assumed to be diffuse, while it is possible for a room to have different decays throughout the space. A different approach was suggested by Dalenbäck

72

. This algorithm is based on

approximate cone tracing, where purely specular contributions are accounted for in a primary pass. During the primary pass, the energy to be diffusely reflected is stored at square wall patches and the energy to be specularly reflected continues. During the secondary and subsequent passes, these patches in turn act as sources for diffusely reflected sound, and energy that will be diffusely reflected is, once again, separately stored at the patches. The problem is in principle recursive, but by using a finite number of patches and exploring the properties of diffuse reflection, recursion is abolished. A list is then assigned to every patch. When a primary ray encounters a patch, to that list is added the arrival time and the diffuse energy fraction δ (1 − α ) . When a ray encounters a receiver, the whole list of reflections associated with the patch, now acting as a diffuse source, is added to the echogram. The algorithm contains some approximations that can lead to errors. For instance, since the centre of the patch is used as secondary source instead of the exact point of impact of the ray, the time origin of the diffuse reflection is

42

somehow translated in time. The largest uncertainty will be for grazing incidence, but this can be diminished for the most important surfaces by increasing the patch density. Refinements are made, creating a higher number of rays for patches that are excited earlier, and modeling the patch density dependent on the factor δ (1 − α ) (high density for hard surfaces with strong diffusion). Dalenbäck established the long computation times required by this model. Stephenson

53, 62

radiosity method

tries to avoid the explosion of computation time by applying the 72

, where the energy flows are not only split but also re-unified. In

principle, the radiosity method is not consistent with ray tracing, as the coincidence probability is very small, and the computation time will explode before recombination occurs. The model traces pyramids, which can have many edges and can be clipped by obstacles. A good point of the algorithm is the suggestion to model diffraction, where a pyramid can be clipped by the edge of an obstacle into three new ones. Stephenson admitted the difficulty of implementing the algorithm, recognized problems of dynamical storage and access to pyramid data, and affirmed that the total algorithm was yet to be tested. A user-defined transition order, TO, has also been suggested, i.e. as in the commercial software ODEON 2.5 73, 74. Reflections with lower orders than TO are considered purely specular. For larger orders only the usual ray tracing is done with an evaluation of δ . However, the concept of transition order is not physically satisfactory since diffuse reflections should occur even at the very first reflection rather than suddenly being switched on at the transition order 16, 43. Experimental studies do not support the use of TO12. Lam 16 compared the accuracy of three different computer models and their dependence on δ in the prediction of the acoustical parameters of two rooms very different in shape (cube and disproportionate), but composed of the same wall material. He found that the three models need different values of δ to achieve convergence with the measured values, and what is worst, that the selection of δ in all of them changed with the shape of the room, even though they were made of the same type of material. That made him conclude that there is much work still to be done on solving the relationship between the physical scattering properties of real surfaces and the scattering coefficient used in computer models. Moreover, as stated by Dalenbäck et al 43, a disturbing fact is that in

43

some halls the choice of diffusion factors has a greater impact on the estimation of acoustic parameters than the uncertainty in absorption factors. A final possibility to model diffuse ray reflections would be to model the reflected rays as having a certain dependence on the specular direction, i.e. forcing the direction of the reflected ray to have a direction of ± 20° with the angle obtained by Snell’s law.

More information about different ways to account for diffusion in computer models can be found in Dalenbäck et al 43, where a good overview is given. 2.2.1.2. Reflection models with secondary diffuse sources

A different way to look at a diffuse reflection process is to consider every point where reflection occurs as a new secondary sound source, which emits the diffuse reflected sound in the form of new rays with some directional characteristics

45, 62

. However,

direct implementation of a split up ray algorithm would result in an exponential increase in the calculation time, since even secondary rays, upon encountering a diffusing surface have to be split up. Because the details of the sound become less audible with time, it is perhaps not very useful to expend so much computational effort to follow an exponentially increasing number of particles, even while the total energy remains constant

54, 61

. Instead, a secondary source can be created at each subsequent wall

reflection, which will radiate into a hemisphere inside the room as an elementary area source. In this way, all visible receivers are irradiated with diffuse energy, but then the energy is re-grouped into a ray and which is traced forward along a new direction, again by creating a random number and comparing it with δ , to choose between specular or diffuse reflection. This is called the Naylor method

43, 75

. With this approach a high

density of reflections at a receiver is achieved without the need of tracing a very high number of rays. The rays are only losing energy due to wall and air absorption, and their reflected directions are randomized if the walls encountered are diffuse. Both the high reflection density and the spatial information make the approach feasible for auralization. In an ideal diffuse reflection the sound is reflected in all directions over the hemisphere, and all directions are equally probable, see Eq. (2.6). Integrating this equation into spherical coordinates over the hemisphere we obtain 2π , this being therefore the influence at any point a distance r away from that secondary source 1 / 2πr 2 .

44

Ir I r , max

=1

(2.6)

Many room acoustics prediction models consider the diffusion of sound to follow Lambert’s cosine law

13, 72

. Lambert’s law for the diffuse reflection states that the

directional distribution of the reflected intensity is proportional to cosθ r , Eq. (2.7). Ir I r , max

= cos θ r

(2.7)

Integrating now Eq. (2.7) over the hemisphere we obtain π , being therefore the influence at any point a distance r away form the secondary source cos θ r / πr 2 . Although Lambert’s cosine law has been used in many room acoustics prediction methods, its use has been debated since it originates from the reflection of light. On the other hand, Dalenbäck

72

considers his proposed diffuse model to be more efficient if

exit angles are independent of incidence angles as in the Lambert model. Some authors no longer consider Lambert's law a valid approximation 76, 77. Recently, Rindel

76

has proposed a new approximation to model the directional

characteristics of sound reflections, which takes into account the wavelength and the finite area of a rectangular reflecting area, and is valid for all angles of incidence including the grazing 0° . The model is based on sound transmission theory and the use of Babinet’s principle, which states that the reflection from a plane rigid surface is equivalent to the transmission through an opening of the same geometrical shape surrounded by an infinite rigid baffle. In Figure 2.7, the size of the surfaces is 2a ⋅ 2b , and the direction of the incident sound is defined by the angles α 0 and β 0 relative to the x- and y-axis. In a similar way, the direction towards the receiver point is defined by the angles α and β . The wavenumber of sound is k = 2πf / c [m-1], (f is frequency and c the speed of sound). Using the result for sound transmission through a rectangular

opening 78, the intensity of the reflected sound relative to the maximum value is:

Ir I r , max

⎛ sin X sin Y ⎞ =⎜ ⋅ ⎟ Y ⎠ ⎝ X

2

(2.8)

45

where X = ka (cos α − cos α 0 )

(2.9)

Y = kb(cos β − cos β 0 )

Figure 2.7: Definition of angles of incidence and reflection from a rectangular surface. After Rindel 76.

The radiation function (sin X / X ) is valid for a single frequency, but averaging to the 2

frequency band, it is found that the centre frequency is not a bad approximation. Still it has some drawbacks, e.g. X > 2π and X = 0. A mathematical approximation to the octave band average, which does not have the same drawbacks was also suggested by Rindel.

⎛ sin X ⎞ ⎜ ⎟ ⎝ X ⎠

2

≈ (cosh(0.67 ⋅ X ))

−2

(2.10)

The suggested approximate directivity function for different values of the relation ka is plotted for an angle of incidence of 90° in Figure 2.8a, and for an angle of incidence of 30° in Figure 2.8b. In addition, Lambert’s law of diffuse reflection has been added for comparison. It is observed that the reflections have high directivity in the direction of the specular reflection for the high ka values, whereas the reflections are radiated with a nearly uniform directivity for the low ka values. The Lambert model is not particularly close to any of the examples.

46

Figure 2.8: Directivity functions. a) for the angle of incidence α 0 = 90° . b) for the angle of incidence α 0 = 30° . After Rindel 76.

2.2.2. Modeling diffraction in the ray-tracing method Most of room acoustic modeling algorithms rely on geometrical acoustic techniques which completely neglect the diffraction effects. This is because such methods assume that sound propagates as light, i.e. that the wavelength of sound is much shorter than the dimensions of the reflecting objects 17. In section 1.3.1, we introduced the dependence of the wavelength of sound on diffraction effects, i.e. how low frequencies bend at the edge of the barrier. As barriers are frequently used in industrial environments in order to protect workers from noise exposure 20, diffraction effects should be correctly taken into account in order to avoid an over-emphasizing of the barrier effects. Also diffraction appears in auditoriums, especially in seating areas behind balcony fronts cause of the seat-dip effect

12

, or as the

79, 80, 81, 82

. However, as we discuss below, diffraction effects

are very poorly considered by geometrical acoustics models. According to the Geometrical Theory of Diffraction (GTD) proposed by Keller in 1962 83

, an acoustic field, incident upon an edge between two non-coplanar surfaces, i.e. a

wedge, forms a diffracted wave that propagates from the wedge in a cone-shaped pattern, see Figure 2.9.

47

Figure 2.9: Diffraction by a wedge.

Many methods, which are not discussed in this thesis, have been proposed to account for the effect of barriers in a free field. Such methods are not applicable to room acoustics, as also the sound energy reflected from boundaries has to be considered, i.e. the sound reflected from the floor might also suffer diffraction, see Figure 2.10, and the sound reflected from the ceiling obviously decreases the performance of the barrier 84.

Figure 2.10: Diffraction paths in a barrier geometry. The direct sound path is obstructed, and the sound bends at the edge of the barrier.

To simplify this problem, Kurze

85

suggested that the barrier attenuation due to

diffraction is exclusively effective in the case of direct sound. For reflected sound, the barrier is just another obstacle in the room that causes diffuse reflections. Other methods have been proposed to calculate only the effect of the diffracted beam in the shadow region, with the discontinuities in the sound field considered not perceivable 86. Haaren et al 87 validated a ray-tracing model applied in the modeling of noise barriers. Their study can not be considered valid as the validation was done in open field, and only at six points 25 meters away from the barrier. They did not explain how their

48

algorithm works, nor how diffraction was modeled. On the other hand, Svensson et al 88 proposed an analytical method to obtain the impulse response behind a barrier. This method applies the concept of secondary edge sources (every point at the edge of the barrier is considered a sound emitter) and allows the derivation of directivity functions for such edge sources (re-radiated wavefronts spread in all directions with different amplitudes). The applicability of this method relies on the fact that the impulse response inside an enclosure can be written as a sum of geometrical acoustics, hGA, and diffraction components hdiffr. This is good news, but the theoretical development of the method is not easy to follow and an application of it will require a deeper future study. Still, known diffraction formulae assume that rays hit edges exactly, whereas sound rays never do so 53. This problem could be solved as suggested by Dance 20, who, based on an idea of Benedetto and Spagnolo 84, which at the same time was derived from Fresnel and Kirchoff theory 53, proposed three versions of barrier models to apply in ray-tracing methods, and compared predictions of each with measured values. According to the published results, only two of them are interesting. These are: An imaginary plane follows the perimeter of the barrier at a distance λ , where λ is the wavelength of the mid-frequency of the octave band of interest. If a sound ray strikes that ‘diffraction area’ it is randomly redirected propagating along its new path with no attenuation. The second possibility (ten times slower in computation time) is to create a secondary source at the incident point, which is assumed to be omni-directional and to contain all the energy of the incident ray with no attenuation loss. Stephenson has also reviewed these approaches

53

. He said that staying with one particle (first approach) and taking the

directional distribution as a probability distribution would not save computation time to the same extent as more sound rays had to be emitted from the source, requiring even higher computational efforts. Continuing with the experience of Dance

20

, the results

were better for the first approach, with good accuracy at low frequencies, although at 1 kHz, some deviations were observed and attributed to interference effects. Both models failed to predict the negative IL in front of the barrier, and converged with the measured IL at the furthest points, as the barrier effects becomes negligible in the reverberant sound field. The first approach is considered interesting and should one day be implemented by the author of this thesis, together with others, and compared with field data.

49

2.3. Considerations about the ray-tracing algorithm Every model has a certain number of assumptions which might limit its applicability. If the theoretical assumptions do not hold in the case of a particular room for which predictions are to be made, the predictions might not be accurate assumptions which are fulfilled in some space may fail in others

16, 43

26

. In any case,

, and here is where

validations of acoustical models should be checked carefully. The ray-tracing algorithm has been shown to predict noise levels with good accuracy 32, and it is considered one of the most elegant methods to represent diffuse reflections. However, its limitations, restrictions and inherent problems are not few. In this section we discuss and introduce some topics that should be considered when implementing or evaluating the algorithm.

2.3.1. Frequency range The wavelength or the frequency of the sound is not inherent in the ray-tracing model 51. Normally, the only frequency-dependent factors that can be included are the absorption coefficient and the statistical estimation of diffuse properties. Already in their pioneering work Krokstadt et al 11 thought that the most severe restriction of the method is the limited frequency range over which the results are valid. It is an underlying assumption in all methods using sound rays that the wavelength corresponding to the lowest frequency of the sound is small compared with the dimensions of the room and its surfaces. This assumption, i.e. to neglect the undulatory nature of sound, is the reason why experimental studies based on ray tracing

26, 32, 49, 50

have been shown to

predict noise levels and acoustic parameters with good accuracy at high frequencies, while not so good at low ones.

2.3.2. Negation of wave nature, interference and diffraction Therefore, the wave nature of sound is not really considered in the algorithm, and important effects like interference and diffraction are poorly represented or simply ignored. For instance, as we saw in section 1.3.2, in room acoustics there are two kinds of interference effects that should be considered. These are interferences close to the reflecting boundaries, and interferences between highly reflecting parallel walls. The

50

modeling of those could be done only if the sound is considered as what it is; a pressure wave which changes its amplitude and phase at every reflection. But then we would need to use pressure reflection factors allowing a phase shift at reflection 19, and these can only be obtained for certain materials in very ideal or controlled cases, but not in real-life spaces. In this way when several sound field components are superimposed, their mutual phase relations are not taken into account; instead, their energy densities are simply added 13. However, in many cases, e.g. in industrial spaces, the interference effects are negated by the fittings and bodies in the space and the large number of noncoherent sound sources 20. On the other hand diffraction phenomena are neglected in geometrical room acoustics since the propagation in straight lines is its main postulate. In this way, we ignore the effect of diffusion due to the finite size of surfaces, with the consequence of obtaining stronger reflections by geometrical laws than those that occur in reality. There is also one topic which is rarely found in the literature, namely, the seat-dip effect 81, 82

, i.e. the selective attenuation of low frequencies for sound passing over an audience

area at near grazing incidence. This effect, in which multiple diffractions and interferences play a role, appears in the lower frequencies (~ 118 - 178 Hz), and produces a variable attenuation of even 30 dB in magnitude 79, 80. Due to the wave nature of sound, it would be very important for the quality of the simulation model that diffusion and diffraction effects be included in a proper way, in order to achieve more reliable computer models and, in consequence, better results. This would extend the validity of ray tracing to the lower frequencies, where the usually assumed ‘optical limiting case’ is no longer valid and the wavelengths are no longer small compared to the room dimensions.

2.3.3. Angular distribution of reflected sound Generally speaking, the most problematic matter in the geometrical acoustics methods is how to model the angular distribution of reflected sound. It seems that only few options are available, see section 2.2.1, be they specular, diffuse, or a combination of both. However, we should wonder what 'diffuse reflection' really means or what it looks like. Certainly, a lot of experimental work has already been done 69, but to measure and save

51

the data of all existing construction materials has proved to be also a very laborious task. The existing literature 12, 49, 50 tells that models based only on specular reflections lead to large errors and to too deep energy decays. Also, that they over predict the reverberation time, especially in enclosed spaces where the absorption is concentrated on one surface, or where the shape is highly disproportionate, such as large factories 16. On the contrary, approaches that model diffusion by random distribution improve accuracy, but certainly are far from correctly representing the sound distribution of sound after reflection from uneven or finite surfaces.

2.3.4. Geometrical description When modeling a certain enclosure, it is important to decide how precise the description of the enclosures should be. This certainly affects accuracy but also computation time. It has been said that a lot of acoustic experience is necessary to know how detailed a geometrical model must be to provide a basis sufficient for the calculations study

31

50

. One

suggested that to obtain a reasonable degree of precision, it is necessary to

avoid modeling too small surfaces. Another work

49

that compared the efficiency of

different computer models concluded that better results were obtained with simple geometrical descriptions. It a way it is true that the probability of an image of being visible from the receiver decreases with the size of the surfaces involved in its generation

51

, and anyway, what is the point of representing every small detail of the

room when later the reflections are modeled by creating random numbers. However and the other hand, Lisa in his recent work

67

, which studied the acoustic prediction in

ancient open-air theatres, observed that the more complete the geometrical description, the better the results obtained. In this kind of enclosure, where no roof or ceiling is present, this has been previously observed 89.

2.3.5. Frequency dependent factors. The absorption and the scattering coefficients In section 1.7.2, I have introduced the absorption coefficient α , which describes the absorption properties of materials and mentioned that it is an angle-dependent factor. In section 2.2.1.1, I did the same with the scattering coefficient δ , which is used to describe the diffusing capabilities of surfaces. The experimental methods used to

52

measure random sound incidence in the laboratory as described in ISO are introduced in Appendix 1. At this point, we mention that the main reason for avoiding angledependent measurements is that they would need to be very precise, and thus to handle such a large data base of all typical construction materials becomes very difficult. Therefore, estimating the absorption characteristics of a surface by a form that averages over all angles of incidence has been accepted90. Sometimes it is not possible to measure anyone of those coefficients, and then their estimation becomes a guess

16

, which is a recognized source of inaccuracies

12

. The

quality of a computer simulation is then strongly dependent then on the user’s feeling about the absorption and diffusing properties of the modeled surfaces 50. Still, as appointed out in one study

16

differences up to 5 % or even 10 % can be

expected between the absorption coefficients of a certain material measured in a reverberation chamber and the same material placed in the walls of an auditorium.

2.3.6. Number of rays Uncertainties are also introduced in the results by the use of a limited number of discrete rays to represent the sound field 31. As the angle between adjacent rays remains roughly constant, and since all rays are radiated from a point source, this representation gradually becomes less exact with increasing ray length. To achieve convergence criteria (after all, this is a Monte Carlo method) the number of rays needs to be large (even 100,000 in a large auditorium 47). The larger the number of rays, the smaller the solid angle that an individual ray will cover. Again, this will affect the computation time, but a too low number of rays can cause an insufficient number of detected reflections 91.

2.3.7. The detection problem Three kinds of receivers can be used in room acoustics simulation, point, plane and volume, but one basic condition: The traced rays must meet a detection condition at the receiver's position in order to contribute to the impulse response of the room at that point 91. The probability of a ray being detected is a function of the size of the detector, and point receivers are not adequate in ray tracing. An infinitesimally thin line arbitrarily emitted by a source point almost never hits a receiver point

53

. Plane

reception cells are also inadequate, as rays reaching the reception plane at grazing

53

incidence can not be handled properly 9. Spherical receivers allow an omnidirectional sensitivity pattern, as the cross-section of a sphere is a circle of constant size for all directions of incidence

10, 91

. Spatially anisotropic shapes, like cube receivers, show

more complex directional behavior

91

. Ondet and Barbry

9

considered that spherical

receivers are more difficult to implement and require longer computation time than plane detectors. But this is not true, as a sphere is geometrically characterized only by its center and its size (ratio), while a plane needs at least the coordinates of three points, and the detection problem is surely more complex. In section 3.3.4, how to implement the detection of sound rays by spherical receivers will be covered. Lehnert 91 studied the systematic errors of the ray-tracing algorithm due to the detection problem, concluding that one of them was due to the finite size of the receivers (independently of their shape). In order to achieve a constant density of rays across the solid angle of the receivers, he derived a formula, Eq. (2.11), to obtain the size of the receivers r, for a certain number of rays N and the maximum length of a ray, lmax (which is the largest path a ray can travel within the room). However, large detectors can cause more errors as they may detect some rays too far away from the point of interest.

r = l max (2π / N )

1/ 2

(2.11)

Another possibility that has been suggested is Eq. (2.10), where VR and V are the volume of the receiver and of the enclosure, respectively 92:

V R = 10V / N

Lehnert

91

⎛ 15V ⎞ so r = ⎜ ⎟ ⎝ 2πN ⎠

1/ 3

(2.12)

also called attention to the fact that due to the finite size of the detectors,

sound paths may appear which are invalid in terms of the geometrical model, see Figure 2.11. Some of them (first and second) can be interpreted as rough models of diffraction and diffuse reflections, and can as such be tolerated. Others (third) are definitely intolerable.

54

Figure 2.11: Examples of invalid detections; a) The sound source is shaded. b) The reflection point is not located within the wall. c) The receiver is on the back side of the wall.

2.3.8. Other problems Ray-tracing programs execute in ‘batch’ mode, taking seconds or minutes to update new locations, and therefore are not suitable for real-time auralization. Nor can they scale to support large 3D environments, and only do an inspection of certain predefined propagation paths, as a function of the established sound source and receiver locations 42

.

2.4. Round Robin As we have seen throughout this chapter, many approaches have been suggested to model room acoustics, and therefore many commercial softwares have also appeared in the international market. In Appendix 2 a review of software is given. Those programs are in use today for the purpose of research, training and consulting. Different approaches may lead to different results, and then it is necessary to know, which ones are accurate enough to be accepted as able to represent the sound distribution in enclosed spaces. The PTB institute published the results of the first international round robin on room acoustical computer simulations in 1995

49

. Not only could the accuracy of different

models be checked, but also the ability of researchers, acousticians and/or architects to make use of them. Sixteen participants from 7 countries, and 14 different programs for room acoustical computer simulations were compared concerning the results of 8 different acoustical criteria (RT, EDT, D, C, TS, G, LF and LFC) in a speech

55

auditorium. Two source positions and 5 receiver positions were chosen and the criteria were calculated only for the 1 kHz octave band. The project was run in two steps: first the room geometry and the material data provided on the basis of construction plans and material descriptions in words. No information on the measured results was given. In the second step, an estimation of the absorption and diffuse properties was given to the participants. The results of the first phase showed a surprisingly large scatter with a strong tendency to underestimate the absorption coefficients, and thus to overestimate the reverberation time. In the second phase, the agreement was, of course, better, but even with harmonized input data some results varied considerably from others. Only three of the softwares were considered accurate enough. One interesting point for further investigation was the increasing difference between the results from simulations and measurements at increasing source-receiver distances. The reason could be that in none of the programs is the attenuation of sound waves at grazing incidence over seat rows or the audience (seat-dip effect) taken into account. All programs neglected the wave nature of sound. Programs which modeled only specular reflections were the most outlying. The three ‘winners’ were: The image source method (5th order) plus random tail (35 minutes in processor 486), the hybrid model plus random tail (4 hours in processor 386), and ray tracing (200,000 rays) (140 minutes in processor 486). It is curious to note that these three programs did not make a very detailed description of the room, and nor did they need very long computation times, compared with the other participants. In 1996-98, the PTB organized another international ‘round robin 2’ with 16 participants from 9 different countries and using 13 different computer programs

50

.

This time the object under examination was a concert hall in Jönköping/Sweden named the ELMIA hall, which has an approximate volume of 11,000 m3 and 1,100 seats. The results of this round robin were much better than those of ‘round robin 1’ but the complex design of this hall generated many problems. The ground plan of the hall with the source and receiver locations is given in Figure 2.13, also a photo of it.

56

Figure 2.13: a) Ground plan of the Elmia hall with source and receiver locations. b) Photo of the Elmia hall.

Again, the round robin was done in two phases. In the first one, only a number of photos, drawings, and a description of the surface materials and the climatic conditions were given to the participants. Thus, not only the quality of the simulation software is relevant, also the skill of the user as an acoustician. The amount of results data is extremely high, and there is not room for it in this thesis, although it can be obtained from the net 94. As expected, the weak points of the programs are concentrated on the low frequency calculations as diffraction effects have been neglected in the geometrical calculations. In the second phase, the estimated absorption and diffusion were given by the coordinator of the study, so some real differences in the programs should be revealed. For instance, even with the same kind of computer processors, the computation times varied from 42 seconds to 30 hours. In Ref. 94, all the information about the hall geometry and the measured results can be found. This will allow us in the nearest future to know more about the accuracy of our model.

57

3. Implementation of the ray-tracing algorithm In this chapter we present the actual and complete form of our ray-tracing algorithm. From all the methods presented in the previous chapter, the ray-tracing method has been selected not only because it has been shown to predict acoustics with considerably good accuracy 32, but also because its concept is feasible and it has a simple implementation. Some author may claim that his/her proposed way is the best among others, or that he/she already did that already twenty years ago. That is good news, although it has been noticed that normally more time is spent on describing the virtues and accuracy of an algorithm than on the algorithm itself. It was not the aim of this thesis to present an ultimate implementation of the ray-tracing algorithm, nor to state that it is the best method among all the approaches. The propagation of sound may be well defined, but not its representation in room acoustics. Every method has a certain degree of limitations and accuracy, and we are content to comment on the method we chose. It has been also noticed

43

that it seems tempting to extend a method that one already knows

rather than to incorporate a completely new one (or even to rewrite the complete algorithm). So far that is not our case, and we expect that many changes will be made in the future to this research tool. However, the time spent on the development of the model has been very illustrative, and it is time to write down some impressions. The ray-tracing algorithm is an application of geometrical acoustics. The sound power emitted by the sound source is described by a finite number of rays, which are considered as carriers of energy. These rays travel through the air at the speed of sound and are reflected from room boundaries, until their energy, due to the air and wall absorption, has been reduced below a certain limit. During this time, the rays might cross predefined volumes, the receivers, where an energy calculation process is performed. These data build up the energy response, from which desired acoustic parameters can be obtained

24

. Geometry laws are used to follow the paths of all the

rays, and therefore they must be repeated a considerably large number of times, thus causing great computation efforts. To achieve competitive computation times, all the operations involved in the inner loops should be as efficient as possible. A deeper description of the ray-tracing algorithm has already been given in section 2.1.3.1. First, we will look at the schematic flowchart of the model. Then, the algorithm implementation is described. This will be done in two parts, a geometrical and an energetic part, because they are the two main groups of operations involved. In the first

58

part, a description of the enclosure from the geometrical point of view is given, e.g. the surfaces or locations of the sound source and receivers. Also the way the sound distributes, that is, how the rays interact with the surfaces, i.e. the point-in-polygon test, the visibility test, or the way collisions are modeled, i.e. specularly or diffusely. In the second part, the instant energy of every ray is described as a function of the characteristics of the sound source and of their interaction with the medium, i.e. air and walls. Also how those rays ‘impinge’ energy to the receivers and how the energy response and the desired acoustic parameters are obtained. Comments will also be made on which operations are performed in a logarithmic scale while using sound levels or in linear scales.

3.1. Flowchart of the implemented ray-tracing algorithm First, all the input parameters have to be introduced. These are the geometry and acoustic properties of the enclosure, and the locations and characteristics of the sound source and receivers. Then, one by one, create new ray with a direction and an initial power. We trace the ray from the sound source or from the previous collision point, as long as its power Wray is larger than a certain limit Wlimit. Then, we have to check all surfaces, and by applying the visibility test we find out which of them all is collided.

Then, we check receivers to see if the ray crosses any along its actual path. If that happens, the energy that the ray ‘impinges’ to it is calculated, therefore calculate Erec(t). Then, the sound power of the ray is attenuated according to the air and wall absorption, calculate absorption, to finally model reflection by a diffuse or specular reflection.

When all rays have been shot, the impulse response of the receivers is obtained, and from it the calculation of acoustic parameters, SPLRec & RTRec. It should be noticed that the repetition of the process for every octave band has not been represented in the flowchart.

59

Figure 3.1: Flowchart of the implemented ray-tracing algorithm

3.2. Geometrical calculations All the operations and mathematical steps which are related to the geometrical part of the model are described in the following subsections.

3.2.1. Geometrical description of the enclosure Cartesian coordinates are used. A point in the space is then described by its coordinates x, y, z, a line by its canonical form, Eq. (3.1), and a plane by its general equation, Eq.

(3.2). x − x1 y − y1 z − z1 = = l m n

(3.1)

Ax + By + Cz + D = 0

(3.2)

The canonical form of a line has all the information that we need. It contains a point, (x1, y1, z1) (departure point of the ray) and a vector, (l, m, n) (direction of the ray). At every collision with a surface the direction of the ray has to change, also the departure point, which will be the previously calculated collision point at that surface. How to find the collision point of a ray with a surface is discussed in the next section. But first we

60

have to notice that a plane extends to infinite, while a surface does not. A surface then has to be described by a plane and by its corners, or by the segments that define its r edges. The normal vector N ( A, B, C ) is orthogonal to the plane, and if D = 0 the plane lies over the origin of coordinates. In Figure 3.2a, a plane obtained from three points (P1, P2 and P3) and its normal are represented. All the equations in the following use the same nomenclature, even though the given subscripts.

Figure 3.2: a) A plane and its normal obtained from three points. b) Collision point P between a line r and a plane A.

3.2.2. Shooting rays In principle, the way sound rays are emitted can be either randomized or predefined

10

as a function of the directional characteristics of the sound source. Normally, the sound sources used in acoustic measurements have omnidirectional characteristics, which is our case. Therefore, it must be ensured that the generation is almost uniform on the surface of a spherical source 62. The simple assumption of three random generators for the three components of the vector direction of the ray is not completely correct, as that produces a “cube of rays” instead of a sphere. It is possible to ‘cut away’ the corners of the cube discarding each vector with a modulus greater than one. Instead, we have applied a semi-probabilistic generator, in which the sphere surface is mathematically divided into a large number of equal areas (actually 400=20x20)62. This is done by creating two random numbers, rnd1 and rnd2 in the range [0, 1], and projecting their values over the sphere to obtain the components of the ray, as mathematically expressed in Eq. (3.3). The result is equivalent to cutting the sphere with 20 iso-z planes, equally spaced along the z-axis, and then dividing each circle again into 20 parts, as shown in Figure 3.3. Although the single facets have a very different shape, all of them have the

61

same area. In this way, with just two random numbers, rnd1 and rnd2, the directions of 400 rays can be created, each of them discretely shot through one of these sub-areas.

2

l=2

i + rnd 1 ⎛ i + rnd 1 ⎞ j + rnd 2 ⎞ ⎛ −⎜ ⎟ cos⎜ 2π ⎟ 20 20 ⎠ ⎝ 20 ⎠ ⎝

(3.3)

2

i + rnd 1 ⎛ i + rnd 1 ⎞ j + rnd 2 ⎞ ⎛ −⎜ m=2 ⎟ sin⎜ 2π ⎟ 20 20 ⎠ ⎝ 20 ⎠ ⎝ n =1− 2

i + rnd 1 20

for i = 0…19, and j = 0…19, with rnd1∈ [0,1] and rnd 2 ∈ [0,1]

Figure 3.3: Subdivision of the omnidirectional properties of the sound source in facets of equal area. After Farina 62.

3.2.3. Detection of collisions Once the ray has been shot in a certain direction, we have to find out the one and only surface which will be collided with. After that, the ray is reflected and again the same calculation process is performed. How to detect the true reflection points is described in this section. 3.2.3.1. Point-in-polygon test

Mathematically, the collision point xP, yP, zP, between a line, Eq. (3.1), and a plane, Eq. (3.2), is calculated by Eq. (3.4). ⎡ x P ⎤ ⎡ x1 ⎤ ⎡l ⎤ ⎢ y ⎥ = ⎢ y ⎥ − x1 A + y 1 B + z 1 C + D ⎢ m ⎥ ⎢ P ⎥ ⎢ 1⎥ ⎢ ⎥ lA + mB + nC ⎢⎣ z P ⎥⎦ ⎢⎣ z1 ⎥⎦ ⎢⎣ n ⎥⎦

(3.4)

62

Because every surface is bounded, the next question to answer is whether the collision point is inside or outside the geometrical perimeter of the surface. This is called the point-in-polygon test. In principle, this test has to go over every surface after every ray

reflection, so it should be meticulously implemented as it is repeated many times in the inner loop. One definition of whether a point is inside a polygon is the Jordan curve theorem. Essentially, this says that the point is inside the polygon if, for any line that we trace from it, there is an odd number of crossings of the line within the polygon’s edges, see Figure 3.4a. This definition also solves the problem for areas that are inside the polygon.

Figure 3.4: a) Jordan Curve theorem. b) Resolution of Jordan theorem in the plane oxy.

On the other hand, the polygon and the sound ray are naturally defined in three dimensions, and this can lead to a problem that is now discussed. The problem is related to the floating point precision of the used programming language, i.e. with the determination of the decimal digits. As every rational number is rounded, results can appear which differ from the analytical results. For instance, when applying the Jordan curve theorem to perform the polygon test, a line is traced from the collision point

towards one of the polygon edges. This line and the collision point should lie over the surface, but perhaps no longer do because of the rounding. In this way, the computer program fails when it just lightly considers that two lines are not in the same plane and therefore never intersecting with each other. All this leads to likings, sound rays that never hit a surface and escape from the room. To avoid this problem and also to simplify computation, it is worthwhile to project the polygon and the test point into two dimensions. One way to do this is to simply ignore one component. The best component to ignore is usually that which, when ignored,

63

gives the largest area polygon. This is easily done by taking the absolute value of each component of the polygon plane’s normal and finding the largest. The corresponding coordinates of the polygon and test point are then ignored. If the polygon is not parallel to the coordinate planes, oxy, oyz, oxz, in practice, it does not matter which component is ignored. Precomputing this information once for a polygon uses more memory but increases the speed of the intersection test itself. For instance, when we project the line described in Eq. (3.1) over the plane oxy, the z components disappear, see Figure 3.4b. Then we can express that line in a more convenient way, as in Eq. (3.5). Note: Do not confuse with Eq. (3.2). Ax + By + C = 0

(3.5)

The coordinates (x0, y0) of the intersection point between two lines in 2D are obtained by solving their equations simultaneously. If these lines are the line that we trace from the collision point to perform the Jordan theorem and the line defined in one of the edges, which are given by A1 x + B1 y + C1 = 0 and A2 + B2 + C 2 = 0 , then the intersection point is

x0 =

and if

A1

B1

A2

B2

B1

C1

B2

C2

A1

B1

A2

B2

and

y0 =

C1

A1

C2

A2

A1

B1

A2

B2

(3.6)

= 0 , the lines are parallel.

We should repeat this process with all the edges of the polygon and count how many are crossed within the corner limits to complete the point-in-polygon test. 3.2.3.2. Visibility test

It has been explained how to obtain the collision point between a line and a plane, and then how to know if the point is inside the surface’s polygon or not. That is actually the second of the three Boolean functions that compose the visibility test. However, the remaining two are even simpler. The three functions are represented in Figure 3.5, where T represents true and F false. The first ensures that the collided surface is in

64

accordance with the direction of the ray and not in the opposite direction, see Figure 3.5a. The second has already been explained, and the last selects from all still possible solutions only the closest, see Figure 3.5c.

Figure 3.5: The three steps of the visibility test. a) Right direction of the ray. b) Pointin-polygon test. c) Only the closest of all the possible collided surfaces. T represents True and F represents False.

3.2.4. Modeling reflections As we saw in chapter 2, the reflection of sound from real surfaces does not generally follow Snell’s law. This means that the reflections are never perfectly specular, instead they are partially diffuse. This is a consequence of the natural roughness of construction materials and of the finite size of surfaces. Also we said in section 2.2.1.1 that the most popular method for modeling reflections in the ray-tracing algorithm is the use of the scattering coefficient δ , and this is the method we use. When a ray encounters a diffusing surface, a random number in the range [0, 1] is generated. If the number is smaller than δ the ray direction is randomized to simulate diffusion, otherwise the reflection is specular

13, 72

. The problem is; how well will the real situation be

reproduced by just a normal randomization, Figure 3.6a, or by a weighted randomization closer to a specular reflection, Figure 3.6b. We have simply chosen the first possibility where all directions are equally probable in the case that reflection is diffuse.

Figure 3.6: a) Diffuse reflection by normal randomization. b) Diffuse reflection by weighted randomization.

65

Just to illustrate the consequences of the applied method we have presented in Figure 3.7 the possible path followed by a ray. The energy of the ray has been represented by the thickness of the line, which decreases at every reflection, jn. The first three reflections, j1, j2, j3 are specular, while the others are diffuse.

Figure 3.7: Representation of specular and diffuse reflections in a 2D enclosure. 3.2.4.1. Specular reflection

The three directions of the new vector direction of a ray after a specular reflection from a surface are calculated from: ⎡ l´ ⎤ ⎡ l ⎤ ⎢m´⎥ = ⎢m ⎥ − 2 lA + mB + nC ⎢ ⎥ ⎢ ⎥ A2 + B 2 + C 2 ⎢⎣ n´ ⎥⎦ ⎢⎣ n ⎥⎦

⎡ A⎤ ⎢B⎥ ⎢ ⎥ ⎢⎣C ⎥⎦

(3.7)

3.2.4.2. Diffuse reflection

To model a diffuse reflection, some alternatives are available. One way is to make use of Eq. (3.7), and check what components have changed their sign. The sign of those ones should be kept, but the values could be randomized. This is what we have done, although another way would be to also make use of Eq. (3.7), and then introduce some randomness to every one of the components, i.e. by creating random numbers in the interval [-1, 1], dividing them e.g. by 20, and add these values to the components of the vector (this method is closer to what was illustrated in Figure 3.6b). There are some other possibilities, however, a simple randomization of the vector components can lead to obvious errors, and it is important to ensure which components have changed their sign upon reflection. One last thing to be mentioned is that if the scattering coefficient is constant for all the frequencies, the tracing rays process can be done in a single pass. However, the absorption of sound is not constant along the frequency range.

66

3.2.5. The image method If, for the early reflections the image method is desired, its implementation is still simple. Its results can be superimposed on those from ray tracing in order to see similarities or differences between both in the early reflections. Having a point in the space, P(x, y, z), its image, PIm(x’, y’, z’), is calculated by Eq. (3.8), which translates the coordinates of the point to the other side of the plane. We remember that the concepts of the image method were introduced in section 2.1.3.2, as well as the backwards visibility test, to find out whether an image is real. Second and higher order images, e.g. S’c2-w in Figure 2.2, are obtained by successive applications of Eq. (3.8). ⎡ x´⎤ ⎡ x ⎤ ⎡ A⎤ ⎢ y´⎥ = ⎢ y ⎥ − 2 xA + yB + zC + D ⎢ B ⎥ ⎢ ⎥ ⎢ ⎥ A2 + B 2 + C 2 ⎢ ⎥ ⎣⎢ z´ ⎦⎥ ⎣⎢ z ⎦⎥ ⎣⎢C ⎦⎥

(3.8)

3.3. Energy calculations Once the implementation of the geometrical part has been discussed, it is time to introduce the part that includes the energy calculations. Here we notice already that instead of energy, we could talk about power or intensity, but not about levels as they are a logarithmic and not a linear scale, as is discussed in the following.

3.3.1 The instant power of every ray In section 1.6 we introduced the concept of sound levels. But the problem in the operations that involve logarithmic scales is that they are more time-consuming. Certainly all the calculations could be performed in decibels, but to add two levels requires more computation effort than adding two powers. Therefore, it is wise to just convert the sound power level of the sound source into power, then perform all the operations in a linear scale, and finally change back to decibels if needed. For instance, if we know the sound power level of the sound source for a certain frequency band, we can obtain the corresponding power using Eq. (1.14), and obtaining Eq. (3.9), where W0 = 1 pW. W = W0 10 0.1Lw

(3.9)

67

Once we have the sound source described by its power and not by its sound power level, we can assign to every ray its initial power just by dividing the total of the sound source by the number of rays, N. Avoiding all explanations that have already been given, the i’th ray that arrives at a receiver carries with it a certain power Wi , which will be its initial value minus the part that has been absorbed in all previous paths, Eq. (3.10) 92.

Wi =

W0 Q −mL J e ∏ (1 − α j ) N j =1

(3.10)

where W0 is the power of the sound source, Q is the directivity factor of the sound source (if omnidirectional, Q = 1), N is the total number of rays, and m is the air absorption coefficient, L is the total distance traveled by the ray from the sound source to the receiver, J is the number of reflections so far involved in the path, and α j the absorption coefficient of the wall on which the j’th reflection has occurred.

3.3.2. Reception of sound rays Suppose the energy contributed by the ith ray is Ei , and its sound power is Wi , we can get the following formula: E i = Wi Ti

(3.11)

where Ti is the time that the ray has spent on crossing the receiver, which can be obtained from the velocity and the distance that the ray has crossed within the receiver. Let us illustrate how to calculate this distance by looking at Figure 3.8a, where a spherical receiver is crossed by one ray. Three quantities are illustrated there; dri the distance traveled through the receiver, d the shortest distance from the ray to the center of the receiver, and r the radius. Because the distance d from a line to a point can be calculated from Eq. (3.12) 95, and because the radius of the sphere receiver is known, we can apply the Pythagoras

68

theorem and obtain the desired distance as d ri = 2 r 2 − d 2 . If the distance dri is smaller than the radius of the receiver, the ray has crossed it.

⎡ [(a − x1 )m − (b − y1 )l ]2 + [(b − y1 )n − (c − z1 )m]2 + [(c − z1 )l − (a − x1 )n]2 ⎤ d =⎢ ⎥ l 2 + m2 + n2 ⎣⎢ ⎦⎥

1/ 2

(3.12)

Figure 3.8: a) Sphere receiver. b) Cube receiver.

Finally, it can be derived that the relationship between the sound energy of a ray and the sound intensity ‘impinge’ to the sphere receiver is 92:

Ii =

E i ⋅ c Wi ⋅ Ti ⋅ c Wi ⋅ d ri = = Vr Vr Vr

(3.13)

The corresponding arrival time of every ray at detection will determine where the intensity is allocated in the intensity response.

3.3.3. Calculation of acoustic indices Once all the rays have been traced, we can consider the impulse response obtained at every receiver almost complete. The last thing that we have to do is to change back to decibels by the use of Eq. (1.19). The reason for doing this is because the different acoustic indices that we can calculate are obtained using sound levels 24.

69

The impulse response, see Figure 1.9, could be saved on disk in two different ways. The first includes two columns of data: one for the intensity or energy of every pulse, and the other for their corresponding time delay. Another possibility is to divide the impulse response into certain ∆t steps, e.g. 1/44.1 kHz. Thus, only one column (or a row of cells) is needed and that is how we have done. Two acoustic parameters are introduced, the sound pressure level, SPL, and the reverberation time, RT. These parameters provide information about the acoustic quality of an enclosure, and also allow a good comparison of measured and modeled data. Both can be obtained directly from the impulse response, see chapter 5 and Figure 5.1. Other parameters can be obtained from the impulse response and are defined in Ref. 24. 3.3.3.1. Sound pressure level

Sound pressure level, SPL, is the experimental measurement of the amount of pressure in the air being sensed at a given location. In the ray-tracing algorithm, it can be obtained from Eq. (3.14), although it has the same value as the total intensity or energy expressed in dB. ∞ SPL = 120 + log10 ⎛⎜ ∫ I (t )dt ⎞⎟ ⎝ 0 ⎠

(3.14)

3.3.3.2. Reverberation time

Reverberation time, RT, is the time expressed in seconds, that would be required for the sound pressure level to decrease by 60 dB after the sound source has stopped. The accuracy in reverberation time measurements is limited due to random fluctuations in the decay curves, which are produced as a result of the randomness of the excitation signal. The same effect occurs in modeled impulse responses. A popular method minimizes the effect of the fluctuations in decay curves by repeating the experiment many times and then averaging. This is inefficient, especially when modeling. Instead, we can apply a method proposed by Schroeder 96, the Reverse integration of the squared impulse response, also called the integrated impulse response method, which is more

efficient than spatial averaging and forms the basics of the standard ISO 338224. For each octave band a backward integration of the squared impulse response is performed. In an ideal situation with no background noise the integration should start at the end of

70

the impulse response (t → ∞ ) and proceed to the beginning of the squared impulse response. Thus the decay as a function of time is: ∞

∞

t

t

E (t ) = ∫ p (τ )dτ = ∫ p (τ )d ( −τ ) = ∫ p (τ )dτ − ∫ p 2 (τ )dτ 2

2

t

∞

∞

tn

2

0

and

/

0

or ∞

tn

E (t ) = ∫ p 2 (τ )dτ = ∫ p 2 (τ )dτ + ∫ p 2 (τ )dτ = ∫ p 2 (τ )dτ + C t

t

tn

(3.15)

t

In chapter 5, i.e. Figure 5.1, we will see how the Schroeder method applies in modeled impulse responses. The reverberation time is then obtained by a least-squares fit line to the decay obtained by Schroeder over the range from 5 dB to 35 dB below the initial level, represented as T30, or from 5 dB to 25 dB, represented as T20.

71

4. Materials and methods Three different enclosures have been studied. Measured data and modeling results of two acoustic parameters, sound pressure level SPL and reverberation time RT, have been compared in order to evaluate the accuracy of our algorithm. In this chapter the geometrical and the acoustical descriptions of the enclosures are given. Also the way the measurements were performed and how the modeling was carried out are described.

4.1. Description of the enclosures The enclosures have been chosen to cover three different types of indoor spaces. A reverberant room, an industrial hall, and an auditorium. Those spaces are referred in the following as Room 1, Room 2a, Room 2b and Room 3. Previous studies of Room 2 and Room 3 have already been reported 97, 98.

4.1.1. Reverberant room, Room 1 A reverberant room has highly reflecting surfaces (low absorption characteristics), and therefore its reverberation time is fairly long. The dimensions of the enclosure were 2.9 x 7.6 x 3.6 m. The bounding surfaces were of concrete. There were three wooden doors, and a large opening which was also covered with wood. On the ceiling there was an acoustic resonator.

4.1.2. Industrial hall, Room 2a and Room 2b The laboratory of acoustics of the FIOH studied the hall in 1999, during a noise control project 4, 66. For that project, the acoustic characteristics of the enclosure were modified, so there are data available before and after the acoustic treatment was performed. The hall before the acoustical treatment is named as Room 2a and after it as Room 2b. This space was an industrial textile hall with dimensions 40 x 20 x 8.5 m. It had a saddle roof with 5 baulks, 16 weaving machinery pieces of an approximate height of 2 meters, and some shelves. The floor was made of concrete, also the front wall and the two side walls. Special characteristics were a corridor that connected to another room through the front wall, and opposite it, an opening that led to another large hall, which is not included in the model. The acoustic treatment that was carried out consisted of covering the two side walls and the ceiling baulks with 50 mm mineral wool. A photo of the industrial hall is shown in Figure 4.1.

72

Figure 4.1: View of the industrial hall, Room 2b.

4.1.3. Auditorium Photos of the studied auditorium are shown in Figure 4.2. It had an approximate volume of 1050 m3, and a rectangular ground area of 4.6 x 13 m. The seating area was inclined by approximately 20 degrees with a capacity of 180 seats distributed in 10 rows. The walls were of concrete and the floor was covered by wool carpet. There were absorbent elements in the sidewalls. The same material was suspended from the ceiling by 90-200 cm. Tables were located in the ground area. A strong flutter echo was observed between the side walls using handclaps.

Figure 4.2: Views of the auditorium, Room 3. Omnidirectional loudspeaker at location S1 and receiver at location R1.

73

4.2. Measuring methods The measurements in Room 1, Room 2a and Room 2b were done in 1999 by Hongisto and Keränen. The author of this thesis was directly involved in the measurements of Room 3 in 2004.

4.2.1. SWL of the sound source The sound power level of the sound source, B&K OmniPower 4296, was measured before every study according to ISO 3741 in a reverberant chamber

99

. The source is

proved to be omnidirectional. The obtained values for the six studied octave bands are listed in Table 4.1. Table 4.1: Measured sound power level of the sound source (ref. 1 pW) 125

250

500

1000

2000

4000

Source at Room 1

92.8

93.7

93.7

92.4

92.2

84.4

Source at Room 2

105.1

106.1

106.1

104.2

103.8

96.2

Source at Room 3

103.1

111.1

107.2

102.9

103.0

96.6

4.2.2. Measurement of SPL and RT In Room 1, Room 2a and Room 2b, an omnidirectional loudspeaker (B&K 4296), a power amplifier (QSC 1300 W USA) and a pink noise signal generator (Behringer DSP 800) were used as the sound source. Sound pressure levels, SPL, were measured using a precision sound level meter (B&K 2260A) equipped with a ½” microphone (B&K 4189). For the reverberation time measurements a digital tape recorder (Tascam DA-P1) was used to record pistol shots. The shots were recorded via the sound level meter and the recordings were analyzed with a real time analyzer (B&K 2133) using Schroeder’s backward integration method. The RTs were determined using a decay of 20 dB. The final RT values are arithmetic averages of several RT measurements. In Room 3 the sound pressure level was measured in 1/3 octaves with an Investigator B&K 2260A and a calibrated microphone (B&K 4189). The sound source was the same as in previous rooms. The height of the receivers was the same height as the ears of a possible listener sitting in the audience area, 1.2 meters above the floor. At the same locations, room responses were measured with WinMLS2000 using a sound level meter (B&K 2231 with microphone B&K 4165). From the room responses the RTs were calculated using a decay of 20 dB.

74

4.2.3. Measurement and estimation of the absorption properties of materials As we have mentioned above, it is not always possible to reliably know the absorption properties of materials, and then their estimation becomes a guess. That was the case in Room 1 and Room 2, where a library database was used. On the other hand, we had the

chance to study in the laboratory the most important materials present in Room 3 (wall and ceiling absorbent, and audience chairs). The principles and measuring procedure of ISO 354 are described in Appendix 1, although they are introduced now. The reverberation time of a reverberation chamber is measured. Later, the specimen of the material is placed in the room and the reverberation time is measured again. The difference between the two reverberation times allows the calculation of the absorption coefficients of the surface, and these are called random incidence absorption coefficients.

4.3. Modeling of the enclosures The studies were done only for six octave bands, from 125 Hz to 4000 Hz. We did not model the 63 Hz and 8 kHz octaves because loudspeakers are normally inefficient at those frequencies. As mentioned above, the absorption coefficients were estimated or measured, and the scattering coefficients were always estimated to be constant through the frequency range. The absorption coefficients of the air were taken from ISO 96131100, and the temperature and humidity conditions of the enclosures were modeled as they were when measuring. Because the effect on accuracy of the number of rays and of the size of the receivers was unknown, different numbers of these were used.

4.3.1. Reverberant room, Room 1 In Figure 4.3, some views of the modeled room can be seen. Receivers with a radius of 40 cm were located at the same locations as measured, see Table 4.2. In Figure 4.3, the sound source has been represented by a point at a certain height, and the spherical receivers by a point and a circle. The sound power level of the modeled sound source corresponds with the first row of Table 4.1. Four different construction materials were used in the room. Their absorption characteristics, which are listed in Table 4.3, were estimated from the material database. In Table 4.3, the estimation of the scattering coefficient (Scatt.) for every surface has been included. The last column in Table 4.3

75

illustrates the total area covered by each material in m2, and the last row represents the average modeled absorption and scattering distributed throughout the room in m-2. Table 4.2: Locations of source, S1, and receivers, R1, R2,…,R6, in Room 1 (in meters). S1

R1

R2

R3

R4

R5

R6

X

1.45

1.45

1.45

1.45

1.45

1.45

1.45

Y

1.5

2.03

3.04

4.05

5.07

6.08

7.09

Z

1.50

1.50

1.50

1.50

1.50

1.50

1.50

Table 4.3: Description of modeled surfaces in Room 1. Absorption and scattering (Scatt) coefficients, area in m2, and average in m-2. Walls, Floor, Ceiling Doors (wood) Acoustic resonator Opening (wood)

125 0.01 0.14 0.15 0.15

250 0.01 0.1 0.25 0.11

500 0.01 0.06 0.4 0.1

1000 0.02 0.08 0.55 0.07

2000 0.02 0.1 0.6 0.06

4000 0.02 0.1 0.6 0.07

Scatt. 0.2 0.2 0.2 0.2

Average

0.03

0.03

0.03

0.03

0.03

0.04

0.2

Area 84.34 6.32 0.96 10.16

Figure 4.3: 3D and 2D front and side views of the modeled reverberant room, Room 1. (Dimensions in meters).

4.3.2. Industrial hall, Room 2a and Room 2b The appearance of the room model is shown in Figure 4.4. The locations of the sound source and of the receivers are listed in Table 4.4. The radius of the receivers was 50 cm and their centers were at the same level as the top of all the machinery fittings. For simplicity, all the machines were considered to have similar geometry and absorptive properties. The absorption characteristics of all the modeled surfaces were estimated from the database. The acoustic treatment, which consisted of 50 mm mineral wool, was

76

mounted on the sidewalls and in the ceiling baulks. These data are illustrated in Table 4.5 for Room 2a (before acoustic treatment) and in Table 4.6 for Room 2b (after acoustic treatment). The estimation of the scattering coefficient (Scatt.) was done according to Table 2.1. The scattering coefficient of the floor has been modeled high due to the amount of fittings lying over it. The scattering coefficients of other surfaces have been modeled with a value of 0.5 for simplicity. In Table 4.5 and Table 4.6, the last column represents the total area covered by each material in m2, and the last row represents the average modeled absorption and scattering distributed throughout the room in m-2. The sound power of the sound source corresponds with the second row of Table 4.1. Table 4.4: Locations of source, S1, and receivers, R1, R2,…,R6 in Room 2a and Room 2b (in meters). S1

R1

R2

R3

R4

R5

R6

X

8.5

27.0

9.2

2.5

10.0

27.0

29.0

Y

6.0

6.0

1.4

6.0

18.0

13.5

18.8

Z

1.5

1.5

1.5

1.5

1.5

1.5

1.5

Table 4.5: Description of modeled surfaces in Room 2a (before acoustic treatment). Absorption and scattering (Scatt) coefficients, area in m2, and average in m-2. Floor Walls, Bricks Ceiling, Baulks Corridor, Machinery Open wall

125 0.01 0.02 0.01 0.50 0.99

250 0.01 0.03 0.02 0.50 0.99

500 0.01 0.03 0.02 0.50 0.99

1000 0.02 0.04 0.02 0.50 0.99

2000 0.02 0.05 0.02 0.50 0.99

4000 0.02 0.07 0.05 0.50 0.99

Scatt. 0.7 0.5 0.5 0.5 0.5

Average

0.15

0.15

0.15

0.16

0.16

0.17

0.54

Area 800 860 1101 649 165

Table 4.6: Description of modeled surfaces in Room 2b (after acoustic treatment). Absorption and scattering (Scatt) coefficients, area in m2, and average in m-2. Floor Bricks Ceiling Walls, Baulks (acoustic treatment) Corridor, Machinery Open wall Average

125 0.01 0.02 0.01 0.20

250 0.01 0.03 0.02 0.70

500 0.01 0.03 0.02 0.95

1000 0.02 0.04 0.02 0.95

2000 0.02 0.05 0.02 0.95

4000 0.02 0.07 0.05 0.95

Scatt. 0.70 0.50 0.50 0.50

Area 800 56 801 1104

0.50 0.99

0.50 0.99

0.50 0.99

0.50 0.99

0.50 0.99

0.50 0.99

0.50 0.50

649 165

0.20

0.36

0.44

0.44

0.44

0.45

0.54

77

Figure 4.4: 3D and 2D ground and side view of the modeled industrial hall, Room 2a and Room 2b. Dimensions in meters.

4.3.3. Auditorium, Room 3 The measurements and modeling in Room 3, Figure 4.5, were done over ten points, and repeated for two different locations of the sound source, altogether twenty measuring points, see Table 4.7. The size of the receivers was conditioned by their proximity to the audience chairs, as they should not overlap such planes. Their ratios were 38 cm and their height corresponded with the height of a listener in the audience, 120 cm above the floor. Ten different types of materials were used to model the auditorium. Their estimated absorption and scattering coefficients are listed in Table 4.8. This time we had the chance to study, in our laboratory, the absorptive performance of the most important materials, according to ISO 35490. These data correspond with the first three rows of Table 4.8. The wall absorber was measured with a gap of 7 cm above the floor, and the ceiling absorber with a gap of 90 cm. The rest of the materials were estimated from the material database. Again, the last column in Table 4.8 represents the total area covered by each material in m2, and the last row represents the average modeled absorption and scattering distributed throughout the room in m-2. The sound power of the sound source corresponds with the last row of Table 4.1.

78

Table 4.7: Locations of sources and receivers in Room 3 (in meters). S1

S2

R1

R2

R3

R4

R5

R6

X

2.85

1.01

4.15

4.25

5.20

5.97

8.05

8.12

9.01 11.07 12.01 13.12

Y

11.30 7.19 10.12 4.56

3.06 11.26 9.60

3.58

7.41

3.55 10.80 7.98

Z

1.50

1.65

3.00

3.35

4.35

1.50

1.20

1.20

2.10

3.00

R7

R8

R9

4.80

R10

5.20

Table 4.8: Description of modeled surfaces in Room 3. Absorption and scattering (Scatt) coefficients, area in m2, and average in m-2. Wall Absorber Ceiling Absorber Chairs, soft side Chairs, hard side Concrete Walls Hard Parts in Ceiling Stairs (soft carpet) Floor (soft carpet) Tables Doors Average

125 0.09 0.45 0.13 0.09 0.01 0.10 0.03 0.03 0.09 0.09 0.13

250 0.35 0.45 0.45 0.11 0.01 0.10 0.04 0.04 0.11 0.11 0.17

500 0.85 0.52 0.75 0.09 0.02 0.10 0.09 0.09 0.09 0.10 0.24

1000 0.63 0.60 0.94 0.07 0.02 0.10 0.20 0.20 0.07 0.07 0.29

2000 0.38 0.50 0.94 0.06 0.02 0.10 0.34 0.34 0.06 0.06 0.30

4000 0.47 0.60 0.91 0.07 0.03 0.10 0.45 0.45 0.07 0.07 0.35

Scatt. 0.4 0.6 0.8 0.8 0.1 0.4 0.7 0.2 0.4 0.2 0.48

Area 32 192 68 116 210 12 157 94 28 6

Figure 4.5: 3D and 2D ground and side view of the modeled auditorium, Room 3. Dimensions in meters.

79

4.4. Validation methods The validation of any acoustic modeling method has to be done by a comparison with measured data, and by an evaluation of its accuracy. However, it should be pointed out that there are no perfect measurements, and that three types of factors affect the accuracy of a measurement, i.e. human, device and external factors. The first type occurs when two different persons make a set of measurements or when one repeats his/her own ones. It might be that the measurement locations are not the same, or neither is the orientation of the microphones. The second type is due to the devices. The accuracy of the sound source SWL, and the sound level meter accuracy, are aspects to be considered. For instance, the microphone used in Room 3 was a free field microphone that is optimized to sound arriving from directions straight to it. On the other hand, the modeled spherical receivers are considered to have omnidirectional sensitivity. Finally, there are external factors like changes in temperature and humidity conditions, or changes in the background noise. Normally, it is assumed that the accuracy of acoustic measurements in 1/1 bands is within ± 3 dB in a reverberant room 97

, and ± 5 dB in the open field 100, 102.

Sound pressure levels SPL and reverberation times RT were measured. In modeling an impulse response was obtained for every receiver and at every frequency octave. From these SPL and RT were calculated. The SPL at a point is calculated by adding all the energy pulses of the corresponding point response, Eq. (3.14). Both measured and modeled data can be compared. Their difference at a certain receiver location i and frequency f is obtained from Eq. (4.1). This can be averaged over all the receivers to obtain the mean value of the deviation, Eq. (4.2). Applying Eq. (4.3), we can obtain the standard deviation of the difference between what is predicted by modeling and what is measured at any frequency octave. Ai , f = SPLi , f , Modeled − SPLi , f , Measured

ASPL f =

1 N

N

∑A i =1

i, f

(4.1)

(4.2)

80

⎞ ⎛ N N ∑ A − ⎜ ∑ Ai , f ⎟ i =1 ⎠ ⎝ i =1 STD f = ± N ( N − 1) N

2

2 i, f

(4.3)

The comparison between modeled and measured RT is done visually from the plots.

81

5. Results 5.1. Reverberant room, Room 1 In figure 5.1, four modeled point responses in the reverberant room, Room 1, are plotted. Those plots represent the sound decay of the enclosure at a certain position. They have been obtained at two different receivers, R2 (left) and R6 (right), for two different frequencies, f = 125 (upper) and f = 2000 Hz (lower). SPL and RT can be calculated directly from the data of the plot. SPL is obtained by direct application of Eqs. (1.18) and (1.19) over all the intensity pulses. RT is calculated by applying its Schroeder backwards integration method and fitting a line to it by the least-squares method over the range 5 to 35 dB below the initial level, see section 3.3.3.2.

Figure 5.1: Point responses obtained in Room1 in receivers, R2 (left) and R6 (right), at the frequency octaves f = 125 Hz (upper) and f = 2000 Hz (lower). x-axis represents time with a time interval of ∆t = 0.01 seconds, while the y-axis represents SPL in dB.

The difference between modeled and measured SPL in Room 1 has been illustrated in Figure 5.2, where the middle point represents what was obtained by Eq. (4.2) and the

82

error bars what was calculated by Eq. (4.3). The same calculations were repeated to compare also the prediction by the Sabine theory, Eq. (2.1), with the measured data.

Figure 5.2: Room 1, averaged difference in SPL between modeled and measured, and between Sabine and measured. (mean value by Eq. (4.2) and error bars by Eq. (4.3)). N = 1000 rays and r = 40 cm.

5.2. Industrial hall, Room 2a In Figure 5.3, six modeled sound decays of the industrial hall, Room 2a (before acoustic treatment) are plotted. They correspond with the point responses obtained at receiver R4 with r = 50 cm. The three sound decays on the left in Figure 5.3 correspond with the 125 Hz octave band. The three on the right correspond with the 2000 Hz octave band. For both octave bands, and from top to bottom, they illustrate the cases N = 3000, N = 10000 and N = 30000 rays. SPL and RT can be calculated in the same way as described in section 5.1 for Room 1. Figure 5.4 represents the difference between modeled and measured SPL by the use of (4.2) and (4.3) in Room 2a. The modeled results correspond to the case where N = 30000 rays and r = 50 cm. As in figure 5.2, the difference between predicted SPL by the Sabine theory and the measured SPL has also been illustrated. Figure 5.5 represents the difference between modeled and measured reverberation time, RT, in Room 2a. The modeling parameters are the same as in Figure 5.4, N = 30000 rays and r = 50 cm.

83

Figure 5.3: Six point responses obtained in Room 2a at receiver R4 with r = 50 cm. On the left for frequency f = 125 Hz, and on the right for frequency f = 2000 Hz. Three cases, N = 3000, N = 10000, N = 30000 rays. x-axis represents time with a time interval of ∆t = 0.001 seconds, while the y-axis represents SPL in dB.

84

Figure 5.4: Room 2a, averaged difference in SPL between modeled and measured, and between Sabine and measured. (mean value by Eq. (4.2) and error bars by Eq. (4.3)). N = 30000 rays and r = 50 cm.

Figure 5.5: Modeled, measured, and Sabine reverberation time in Room 2a. N = 30000 rays.

It is possible to compare the predictions made by our model in Room 2a with the predictions made by the commercial software ODEON v3.1. This is done very lightly simply with the idea of comparing our model, which is nicknamed dBWorks, with any other existing model. Such a comparison is only possible for enclosures, Room 2a and

85

Room 2b, and it is based on the work of Keränen 4. The prediction of SPL by both

methods is compared with the measured data, in Table 5.1 for Room 2a. The comparison is illustrated by the results of the arithmetic mean Eq. (4.2), and the standard deviation Eq. (4.3). Table 5.1: Difference in Room 2a between modeled SPL by dBWorks and by ODEON with measured SPL. Arithmetic mean (ASPL) by Eq. (4.2) and Standard deviation by Eq. (4.3). 125

250

500

1000

2000

4000

1.6

1.7

1.5

1.0

0.8

0.7

STD

1.4

1.3

0.9

0.6

1.2

1.0

Mean (ASPL)

1.2

1.7

1.3

0.7

0.2

-0.6

STD

1.1

0.8

0.9

0.7

1.1

1.1

dBWorks Mean(ASPL) ODEON

5.3. Industrial hall, Room 2b In Figure 5.6, six modeled sound decays of the industrial hall, Room 2b (after acoustic treatment) are plotted. They correspond with the point responses obtained at receiver R4 with r = 50 cm. The three sound decays on the left in Figure 5.6 correspond with the 125 Hz octave band. The three on the right correspond with the 2000 Hz octave band. For both frequency octaves, and from top to bottom, they illustrate the cases N = 3000, N = 10000 and N = 30000 rays. Note that the time axis has changed with respect to

Figure 5.3. Figure 5.7 represents the difference between modeled and measured SPL in Room 2b, where the middle point is the difference between both data, Eq. (4.2), and the error bars its standard deviation, Eq. (4.3). The modeled results correspond to the case where N = 30000 rays and r = 50 cm. As in Figure 5.2, the difference between predicted SPL by the Sabine theory and measured SPL has been also illustrated. Figure 5.8 represents the difference between modeled and measured reverberation time, RT, in Room 2b, when N = 30000 rays and r = 50 cm.

86

Figure 5.6: Six point responses obtained in Room 2b at receiver R4 with r = 50 cm. On the left for frequency f = 125 Hz, and on the right for frequency f = 2000 Hz. Three cases, N = 3000, N = 10000, N = 30000 rays. x-axis represents time with a time interval of ∆t = 0.01 seconds, while the y-axis represents SPL in dB.

87

Figure 5.7: Room 2b, averaged difference in SPL between modeled and measured, and between Sabine and measured. (mean value by Eq. (4.2) and error bars by Eq. (4.3)). N = 30000 rays and r = 50 cm.

Figure 5.8: Modeled, measured, and Sabine reverberation time in Room 2b. N = 30000 rays.

In the same way as for Room 2a, the SPL prediction by our model, dBWorks, and by the commercial software ODEON v3.1 have been compared in Table 5.2.

88

Table 5.2: Difference in Room 2b between SPL modeled by dBWorks and by ODEON with measured SPL. Arithmetic mean (ASPL) by Eq. (4.2) and Standard deviation by Eq. (4.3). 125

250

500

1000

2000

4000

-0.6

-1.3

-1.1

-0.4

1.2

-6.7

STD

2.7

2.0

2.1

2.1

2.2

1.9

Mean (ASPL)

1.3

-1.4

-1.5

-1.2

-1.1

-2.5

STD

2.0

1.3

1.6

1.4

1.5

1.5

dBWorks Mean (ASPL) ODEON

5.4. Auditorium, Room 3 Figure 5.9 represents the difference between modeled and measured SPL and between Sabine and measured SPL in Room 3. The upper plot represents the case when the sound source was at location S1, and the lower plot when the sound source was at location S2. Both cases for N = 10000 rays and r = 38 cm. Figure 5.10 represents the difference between modeled, measured and Sabine reverberation time, RT, in Room 3. The upper plot represents the case when the sound source was at location S1, and the lower plot when the sound source was at location S2. Both cases for N = 10000 rays and r = 38 cm.

89

Figure 5.9: Room 3, averaged difference in SPL between modeled and measured, and between Sabine and measured. (mean value by Eq. (4.2) and error bars by Eq. (4.3)). N = 10000 rays and r = 38 cm. Top: source S1. Bottom: source S2.

90

Figure 5.10: Modeled, measured, and Sabine reverberation time in Room 3. N = 10000 rays, r = 38 cm. Top: source S1. Bottom: source S2.

91

Figure 5.11: Sound decays at receiver R4 with source S2 in Room 3. Left: f = 125 Hz. Right: f = 1000 Hz. x-axis represents time with a time interval of ∆t = 0.001 seconds, while the y-axis represents SPL in dB.

92

6. Discussion 6.1. Reverberant room, Room 1 As a consequence of the low absorption properties of this enclosure, the sound decays are very smooth with fairly long reverberation times, see Figure 5.1. The decays are steeper for higher frequencies. The gradient of the sound decays are independent of the location. According to these data and with the concepts that were introduced in section 1.5, this enclosure fulfils the first requirement for the existence of a diffuse sound field, according to Beranek 27. The second requirement is about the irregularities of the walls for a good scattering of sound. This is not met at all because all surfaces were quite smooth and because there was no furniture. That may be the reason, as will be discussed in the following, for the inaccuracy in the estimation of the SPL at low frequencies. According to Figure 5.2, the mean value of the difference between modeled and measured SPL, Eq. (4.2), is smaller than ± 1.6 dB, and the standard deviation, Eq. (4.3) is within ± 1.5 dB. In the same figure we can see how at low frequencies (125 and 250 Hz) the model underestimates the SPL. One of the sources of potential errors is found to be in the estimation of absorption coefficients, and that underestimation could be due to an overestimation of the absorption properties of the walls. However, this may not be the reason in this case, since the absorption characteristics of the most important walls are already very low. Instead, we could consider the reason to be in interference effects as a consequence of the low diffuse properties of Room 1. It was introduced, in section 1.3.2, that standing modes may appear between highly reflective parallel walls. Although the interference effect was lightly simulated in this case by considering the low scattering coefficient, see Table 4.2, its consequences are not correctly represented. In section 1.3.2, we also introduced the existence of interference effects at reflecting boundaries in reverberant sound fields. That is another reason for the discrepancies between modeled and measured SPL at low frequencies. Due to the dimensions of Room 1, the receivers are 1.45 meters away from the side walls, which is not enough to

ensure that the receivers are inside a node or an antinode. For a frequency of 125 Hz with a corresponding wavelength of 2.6 meters, y / λ = 0.56 . For this relation, the increment of SPL due the presence of standing modes can be up to 3 dB, see Figure 1.8.

93

It is observed in Figure 5.2 how both modeled and Sabine estimation of the sound pressure level are within 0.3 dB. This is obvious since the Sabine theory considers the existence of a diffuse sound field, and this room partly fulfils these requirements. Figure 5.2 referred to the case where the number of rays N was 1000 and the receivers had a ratio r of 40 cm. Topics of interest in room acoustics modeling are; what is the optimum number of rays, or how big the receivers should be. For this enclosure the running process was repeated with N = 3000 rays obtaining the same values within 0.1 dB. Using same amount of rays, but reducing the size of the receivers from 40 cm to 20 cm, again produced results within 0.1 dB of the results illustrated in Figure 5.2. This suggests that in similar enclosures (small, without barriers or furniture, with low absorption and diffuse properties, and with non-locally reacting surfaces) the number of rays does not need to be very high to achieve convergence, and nor does the size of the receivers. Still further, a final test was carried out with N = 1000 rays and receivers with r = 40 cm. But this time the scattering coefficients of all surfaces were changed from 0.2 to 0.05. These results are again within 0.1 dB of the results of Figure 5.2, suggesting that in this type of enclosure, where the sound remains for a considerably long time, the way reflections are modeled is not a determinant factor.

6.2. Industrial hall, Room 2a Figure 5.3 presents six modeled sound decays at receiver R4 in the industrial hall, Room 2a. The three sound decays on the left of Figure 5.3 correspond to the 125 Hz octave

band. The three on the right correspond to the 2000 Hz octave band. For both octave bands, and from top to bottom, they correspond with the cases N = 3000, N = 10000 and N = 30000 rays. Some comments can be made, and we start by comparing the slope of

the sound decays. It is observed that the upper slopes are not as smooth as the others, especially after one second. In section 2.3.6 the importance of the number of rays for convergence in prediction was discussed. A too low number of rays may cause an insufficient number of detected reflections, and/or also, the importance of the energy of every detected reflection may be too high, as the energy of every ray depends on the total amount of rays, N, Eq. (3.10). Those two factors produce inaccuracies in the determination of the sound decay, and therefore also in the estimation of the

94

reverberation time. However, we have to wait for results in Room 2b to see how these inaccuracies increase when the total absorption increases. Also, we can compare how the gradient of the sound decays in the 125 Hz frequency band are not as steep as in the 2000 Hz frequency band. This difference is the result of a slightly larger air and wall absorption at the higher frequency. According to Figure 5.4, the mean value of the difference between modeled and measured SPL in Room 2a is smaller than 1.7 dB. The standard deviation is within ± 1.5 dB. SPL values are normally overestimated by the model, and larger departures

from measurements are found in the lower frequencies, probably as a consequence of neglecting the wave nature of sound. We should remember that the height of the receivers is the same as the height of the fittings, see Figure 4.4. Important diffraction effects, which have not been modeled, may exist as a consequence of the multiple obstacles in the line of sight between source and receiver, i.e. R4 and R5. It has been previously stated 51 that reflections from small surfaces (diffraction due to the finite size of surfaces, sections 1.3.1 and 2.2.2) are generally much weaker than calculated by the laws of geometrical acoustics. In Figure 5.4, we can also observe the difference between Sabine and measured SPL. The mean value of the difference is within ± 2 dB, but the standard deviation is almost ± 4 dB. The reason for this large difference between predicted SPL by the Sabine theory

and measured SPL is due to the fact that the Sabine theory, see section 2.1, takes into account only the distance from the sound source to the receiver, but not the exact location of the receiver throughout the room, as the Sabine theory assumes the existence a of diffuse field, see section 2.1.1. Going back to the discussion about the importance and effect of the number of rays in the calculation, we point out that in the case N = 3000 rays the main value of the difference between modeled and measured SPL, calculated by Eq. (4.2), was ± 1.9 dB, and the standard deviation, Eq. (4.3), was ± 1.5 dB. In the case of 10000 rays, the values were 2.2 and ± 1.7 dB. In the three cases, the trends were similar along every frequency. It seems that the number of rays is not a determining factor in the evaluation of SPL. The answer is yes and not. In cases like Room 1 and Room 2a, the sound absorption properties of the enclosures are low, this allowing the sound rays to travel for a long enough time, and therefore, to achieve a sufficient number of detected collisions

95

without a too large number of rays. Room 2b, where absorption properties are higher illustrates the opposite case, see section 6.3. Figure 5.5 presented the difference between modeled and measured reverberation time, RT, in Room 2a. It can be observed that the measured data have greater values than the others, but the trend of all of them is the same. A description of this enclosure was given in section 4.1.2. It was commented there, how this room led to another hall by one wall, and that this other hall was not modeled. Instead, the wall was modeled as almost totally absorbent. It may be, that part of the sound that reaches that second hall comes back attenuated up to some level and with a certain time delay. This effect may make the reverberation time of the hall to longer. In consequence, the comparison between measured and modeled RT does not tell much about the model’s accuracy. Another reason could be suggested for the general underestimation of the RT, being an overestimation of the absorption properties of the hall. However, that might not be the reason as the prediction of SPL was indeed very good. Anyway, it should be commented now that the measurements of the RT in that room (year 1999) were taken only in three locations, and the locations were different from the locations where SPL was measured, and from the locations that we have simulated. From Figure 5.5, it is also observed that the prediction of RT by the Sabine theory and by the model lead to similar values. In this respect, it can be stated that the Sabine theory has proved to be efficient only in enclosures where absorption is low 25, 26.

6.3. Industrial hall, Room 2b It is important to notice that in Figure 5.6, the time axis has changed with respect to Figure 5.3. This had to be done, as in Room 2b the absorption was higher than in Room 2a. Thus, a visual comparison between the gradient of the sound decays obtained in

both rooms as a function of the amount of absorption has been sacrificed. However, an evaluation of the acoustic treatment effect can be carried out by comparing the corresponding SPL and RT values, which are listed in both figures. From Tables 4.4 and 4.5 it can be seen how the absorptive treatment carried out in the balks and side walls of the industrial hall was not very effective in the 125 Hz octave band. This is why a reduction of ‘only’ 2 dB in SPL but of 0.5 seconds in RT was obtained at that frequency. Instead, in the 2000 Hz octave, the corresponding reduction

96

is of about 7 dB in SPL and 1.4 seconds in RT, as a consequence of the high absorptive properties of the used materials at that frequency. In Room 1 and Room 2a, the number of rays was not a strongly determining factor for convergence in the prediction of SPL. In this case, Room 2b, the effect is bigger, and as we discuss, it is directly related to the amount of absorption throughout the hall. For instance, in Figure 5.7, which referred to the case N = 30000 rays, we can see that the model has predicted, with high accuracy, the SPL with -1 dB, with a standard deviation of ± 2.2 dB. A special case was found in the 4000 Hz, with a large underestimation of the SPL. It is obvious that, at this frequency, the estimation of the absorption properties of the material introduced in the acoustic treatment is greater than what it really was. We no longer consider that frequency, and we concentrate the discussion on the others. When N = 3000 rays was used, the prediction of SPL was underestimated by the model by –2.8 dB, with a large standard deviation, i.e. of ± 8 dB in the 2000 Hz octave. When the running process was repeated with N = 10000 rays, the prediction was better, with – 2.5 dB and ± 2.5 dB, respectively. From all these data, it can be concluded that in enclosures where absorption is considerably high, the number of rays for convergence in the prediction of SPL needs to be large. In Room 2b, 30,000 rays are considered to be good enough. In Figure 5.8, the prediction of RT in Room 2b by the model and the Sabine theory and the corresponding measured data, have been illustrated. A large underestimation is observed in the model prediction. Again, we could think the reason to be in the connection of that enclosure through an opening to another hall. We remember that the other hall has not been modeled and instead the opening has been modeled as an almost totally absorbent surface. It may be that the sound going through the opening to the other hall, comes back attenuated after a certain time, affecting the sound decay. In Figure 5.8, a large standard deviation in the prediction of RT is found in the 4000 Hz octave. From the results of Figure 5.7 and the values in Table 4.5, it was concluded that a large overestimation was made of the absorptive properties of the acoustic treatment at that particular frequency. But checking the individual sound decays obtained at every receiver, it is found that in receiver R6 a too low number of detected reflections was found, producing a poor estimation of the rate of the sound decay. It could be that the

97

used number of rays, N = 30000, was not enough to accurately predict the reverberation time in this enclosure with considerably high absorptive properties. However, before increasing the number of rays it should be confirmed that the receiver detects the most important reflections, which are the first ones. Because diffuse reflections have been modeled since the first reflection of every ray, it could be that some important propagation paths have been missed. This problem could have been solved by the image method.

6.4. Auditorium, Room 3 According to Figure 5.9, the mean value of the difference between modeled and measured SPL in Room 3 when source is at location S1, is smaller than ± 1.4 dB, and the standard deviation is within ± 1.8 dB. When the sound source is at location S2, these values are ± 1 dB and ± 2 dB. Those deviations are acceptable and also smaller than expected, seeming to confirm the reliability of the estimated absorption coefficients. In a diffuse room, the sound is homogeneously distributed in the space and the reverberation time is constant in the room volume. As flutter echo was observed in this auditorium, the room could not be considered as a diffuse space, even when the audience chairs, the stairs, and the wall and ceiling absorbent material have diffusing properties. The small standard deviation of measured reverberation time seems to confirm the partially diffuse properties of the room. However, it must be noticed that to have a large standard deviation in the prediction of RT does not necessarily mean that the accuracy of the algorithm is bad. It is normal that the sound decay in an enclosure depends on the place where it is measured. The ray-tracing algorithm is a statistical method and can not represent correctly the wave nature of sound. In Figure 5.11, it was observed for the frequencies 125 and 250 Hz, that the difference between modeled and measured data becomes large, probably because the ray-tracing algorithm is not able to represent standing modes and other effects of an undulate nature. Such effects are easily observed at low frequencies, and can make the reverberation time to be large. The uncertainty in the estimation of the reverberation time is smaller at low than at high frequencies. The key to this problem is found by looking at the individual sound decays, i.e. the sound decay of R4 with the sound source at location S2, Figure 5.11. First it can

98

be seen that the shape of the decay depends directly on the amount of absorption distributed in the room. In the 1000 Hz octave sound decay, with a time delay of 0.26 and 0.40 seconds, two very high peaks have appeared and affected the Schroeder curve. This receiver is in the lower part of the room, between the concrete walls that were modeled as non-diffuse. It can be suggested that some rays were trapped between both walls without being attenuated, until they finally crossed the receiver. In this case, according to sound speed, the high peaks would correspond to 6 and 10 collisions between both walls since they were shot. At 125 Hz, the decay is not so steep as at 1 kHz and the peaks are not so clearly observed, even when they also exist, like one with a time delay of 0.37 seconds. This effect may also help to explain the differences in results between the two locations of the sound source. Although the results in both cases, S1 and S2, seem very similar, it could be that the accuracy in the determination of SPL and RT has been somewhat conditioned by the location of the sound source throughout the room. As we can observe in Table 4.6 and Figure 4.5, S1 is close to the left wall of the auditorium, being at the same time between parallel walls. The other source, S2, is more in the front space and has some fittings around it. In the first case, it seems that a larger amount of rays have been trapped between the parallel walls, which were highly reflecting and non-diffuse, affecting the determination of the acoustic parameters, especially in the closest receivers. As an experiment, the running process was repeated considering all collisions to be purely diffuse (Scatt = 1). In this way, with no specular collisions, the effect of some rays being trapped between sidewalls was avoided. This made the high peaks disappear, which led to a larger underestimation of the reverberation time. In consequence, the standard deviation in the prediction of RT throughout the hall became smaller. Another thing that should not be forgotten and could explain the underestimation of the reverberation time at low frequencies is the size and height of the receivers (0.38 and 1.2 meters). In Refs. 97 and 98, it was stated that the size of the receivers was an important factor. Smaller sizes could produce greater standard deviations if there is not a large enough number of detected reflections, see section 2.3.8. In this study, the receivers were placed at the same height of a possible listener sitting in the audience area, and their size was affected by their proximity to the chairs. The effects at low frequencies that occur at such close distances to obstacles might not be reproducible at all by simple ray tracing.

99

Over a certain limit, the number of rays, N, did not play a significant role in the obtained results. 5000 rays produced the same deviations as 10,000 and 15,000.

6.5. General discussion The repeatability of calculations led to the same results and deviations. This could indicate that the deviation of the created algorithm is small, and that all created random numbers, i.e. for shooting rays or evaluating diffusion, correspond to homogeneous random distributions. The computation time is not the most important factor, although it has to be taken into account. Obviously, it depends on the used computer processor, so it changes a lot between one computer and another. Parameters that affect computation time are number of rays, number of surfaces, and amount of absorption. As we described in Chapter 3, the rays are traced up to a certain limit of energy. The faster the energy decreases due to absorption, the sooner this limit is reached. It has been interesting to see how the computation time for Room 1 with N = 3000, has been larger than for Room 2b with N = 10000,

even though the number of surfaces is larger in the second case. The

computation time in the first case was 1 1/2 hours, and in the second case 25 minutes, with a Pentium4 © processor. Further comments can be made when comparing the predictions done by our model, dBWorks, with the predictions done by the commercial software ODEON v3.1. For Room 2a, Table 5.1, the deviations of both methods with respect to measured SPL are

very similar. In Room 2b, Table 5.2, the predictions by ODEON are slightly better than the predictions done by our model. However, from these results we could think that our model is quite close in terms of accuracy to existing models. ODEON has proved to be efficient 49, 50, 73, and it is the most popular prediction method in the Nordic countries. It has been frequently stated in previous research works

11, 26, 32

, and observed in this

study, that perhaps the most severe restriction of the ray-tracing method is the limited frequency range over which the results are valid. It is an underlying assumption in the ray-tracing algorithm that the wavelength corresponding to the lowest frequency of the sound is small compared to the linear dimensions of the room and its surfaces. Although at low frequencies the estimation of SPL has been satisfactory, the estimation of RT has

100

been more problematic. But here, there is also a conflict with previous researches. For instance, Howarth and Lam

12

studied the prediction by ray tracing of SPL and RT in

seven different enclosures. The volumes of the enclosures varied between 2,000 and 110,000 m3, and their shapes were rectangular, asymmetric and fan-shape. When averaging over all the enclosures, reverberation time was over-predicted at low frequencies. Interference effects due to the wave nature of sound are negated by the raytracing algorithm. Such interference effects, i.e. standing modes, contribute to increasing the reverberation time in any enclosure where they are present. From that point of view, the underestimation of RT by our model seems more logical than their overestimation. Another interesting topic deals with the way diffusion is modeled. It is presently admitted that accounting for the scattering of sound by surfaces improves room acoustics prediction

13, 16, 47, 49, 50, 56, 57

. However, it seems that the way diffusion is

modeled is not very clear, and many authors 12, 43, 51, 53, 76, 77 have suggested that there is more work to do in this field. This thesis has described how to implement the raytracing algorithm, and has discussed its accuracy. We should be critical with ourselves, and admit that the way we have implemented diffusion has been very simple and rough. First, we have considered a constant scattering coefficient for all the frequencies. Second, the directivity of the reflected sound has been modeled from random numbers distributions. As will be mentioned in the next chapter, the future work of the author of this thesis is directly involved with the way diffuse reflected sound should/could be modeled. Nor has diffraction been modeled, although it might be present in enclosures Room 2 and Room 3. In Room 2, sound arriving at the edge of all machinery fittings and

all ceiling baulks may be diffracted. In Room 3, diffraction effects may appear at the edges of the audience chairs and seats. Not to consider diffraction in that room may have affected the accuracy of the determination of RT at low frequencies.

101

7. Conclusions and future work This thesis has offered the mathematical implementation of the ray-tracing method applied in room acoustics modeling. A discussion about its weak points has also been presented. The prediction by the model of two acoustic parameters, SPL and RT, has been compared with the corresponding measured data in three enclosures. Two main conclusions can be outlined. First, the ray-tracing algorithm can represent room acoustics with acceptable accuracy. Second, the inaccuracies might be due to negation of the wave nature of sound at low frequencies. Diffraction, diffusion and interference should be better taken into account. We have not discovered very much new, but now we have our modeling tool ready and we are about to achieve the actual state of the art in room acoustics modeling by ray tracing. It is also a tool to test new inventions, ideas, and approaches. The future work of the author involves three main lines. First, different series of laboratory measurements would be carried out, in order, as introduced in the previous chapter, to obtain more clues about the directional characteristics of reflected and diffracted sound. These results will be applied in our model in order to obtain a better representation of the phenomena involved in room acoustics. Second, the development of the model should continue. The results of the mentioned measurements and other suggested methods, e.g. Svensson 88, Rindel 76 and Metchel 101, should be applied in the model. Other acoustic parameters can also be obtained from the impulse responses and should be implemented, i.e. speech transmission index STI, and early decay time EDT. Another very interesting topic for further development is the concept of auralization, listening in a virtual environment. Although it has been mentioned throughout this

work, the concept of auralization has not been detailed as it was outside the scope of this thesis, but the reader is referred to Refs. 56, 57 and 58 for an overview and main principles. The implementation of auralization, and we are not far from it, can open very big doors for work-place and auditoria design, if the acoustics of an environment can be listened to before construction. Third, a proper validation of the model will be carried out. It would be desirable to compare the prediction of our model with measurement data in a wider group of enclosures. Such indoor spaces could be two auditoria, two speech halls, two industrial halls, and an open-plan office.

102

The future work of the author is directly related to this topic, and some approximations should be studied, e.g. the recently presented findings by Rindel 76 described in section 2.2.1.2, or the corner source model suggested by Mechel 103.

103

8. Acknowledgments The Finnish Work Environmental Fund is endlessly appreciated for the financial support throughout this work, project No. 103464. Virtual Space 4D and Tekes are also appreciated. The author would like also to thank Dr Valtteri Hongisto, Mr Jukka Keränen and all the personnel of the Laboratory of Ventilation and Acoustics of the Finnish Institute of Occupational Health (FIOH) for their encouraging and intellectual support.

104

9. Bibliography 1

J. H. Rindel. The use of computer modelling in room acoustics. Journal of Vibroengineering, No 3(4), 219-224, 2000. 2

P. C. Eleftheriou. Industrial noise and its effects in human hearing. Appl.Acoust. 63, 35-42, 2002. 3

R. Helenius, V. Hongisto. The effect of the acoustical improvement of an open-plan office on workers. In Proc. The 33rd International Congress and Exposition on Noise Control Engineering (InterNoise 2004). Paper No. 674. Prague, Czech Republic, 22-25 August 2004. 4

J. Keränen, E. Airo, P. Olkinuora, V. Hongisto. Validity of the ray-tracing method for the application of noise control in workplaces. Acustica - acta acustica 89, 863-874, 2003. 5

T. J. Cox, P. D'Antonio. Engineering art: the science of concert hall acoustics. Interdisciplinary science reviews, 2003, Vol. 28, No. 2, 119-129. 6

C. Wang, J. S. Bradley. Prediction of the speech intelligibility index behind a single screen in an open-plan office. Appl. Acoust. 63, 867-883, 2002. 7

R. E. Apfel. Deaf architects and blind acousticians?. A guide to the principles of sound design. Apple Enterprises Press. 1998. 8

L. Cremer, H. A. Müller. Principles and applications of room acoustics. Volume 1. Applied Science Publishers. 1982. 9

A. M. Ondet, J. L. Barbry. Modeling of sound propagation in fitted workshops using ray-tracing. J. Acout. Soc. Am. 85(2), 787-796, 1989. 10

L. Savioja. Modeling techniques for virtual acoustics. PhD thesis, Helsinki University of Technology. Espoo 1999. 11

A. Krokstad, S. Strom, S. Sorsdal. Calculation of the acoustical room response by the use of a ray tracing technique. J.Sound Vib. 8(1), 118-125, 1968. 12

M. J. Howarth, Y. W. Lam. An assessment of the accuracy of a hybrid room acoustics model with surface diffusion facility. Appl. Acoust. 60, 237-251, 2000. 13

H. Kuttruff. Room acoustics. Elsevier Applied Science, London, UK, 3rd edition, 1991. 14

M. Recuero. Acústica arquitectónica aplicada. Editorial Paraninfo. Madrid, Spain, 1999. 15

ISO/DIS 17497-1 Acoustics - Measurement of the sound scattering properties of surfaces. International Organization for Standardization, Genève Switzerland, 2002 16

Y. W. Lam. A comparison of three diffuse reflection modeling methods used in room acoustics computer models. J. Acoust. Soc. Am. 100(4) Pt. 1, 1996, 2181-2192.

105

17

T. Lokki, P. Svensson, L. Savioja. An efficient auralization of edge diffraction. Proc. AES 21st Int. Conf., 166-172. 1-3 June 2002, St. Petersburg, Rusia. 18

J. R. Hassall, K. Zaveri. Acoustic noise measurements. K. Larsen & Son A/S. Denmark, 1988. 19 J. H. Rindel. Modelling the angle-dependent pressure reflection factor. Appl. Acoust. 38, 223-34, 1993. 20

S. M. Dance, J. P. Roberts, B. M. Shield. Computer prediction of insertion loss due to a single barrier in a non diffuse empty enclosed space. Build. Acoust. 1(2) 1994, pp 125-136. 21

R. V. Waterhouse. Interference patterns in reverberant sound fields. J. Acoust. Soc. Am. 27(2) 1955, 247-258. 22

R. V. Waterhouse, R. K. Cook. Interference patterns in reverberant sound fields. Part II. J. Acoust. Soc. Am. 37(3) 1965 424-428. 23

C. Wang, J. S. Bradley. A mathematical model for a single screen barrier in open-plan offices. Appl. Acoust. 63, 849-866, 2002. 24

ISO 3382 Acoustics - Measurements of the reverberation time of rooms with reference to other acoustical parameters. International Organization for Standardization, Genève, Switzerland, 1997. 25

M. Hodgson. Experimental evaluation of the accuracy of the Sabine and Eyring theories in the case of non-low surface absorption. J. Acoust. Soc. Am. 94(2) Pt. 1 1993, 835-840. 26

M. Hodgson. When is diffuse-field theory applicable?. Appl. Acoust. 49(3) 197-207, 1996. 27

L. L. Beranek. Music, acoustics and architecture. John Wiley, New York, 1962, p. 445. 28

C. H. Haan, F. R. Fricke. An evaluation of the importance of surface diffusivity in concert halls.Appl. Acoust. 51(1) 53-69, 1997. 29

H. E. Bass, H. J. Bauer, L. B. Evans. Atmospheric absorption of sound: Analytical expressions. J. Acoust. Soc. Am. 52 (3,2) 821-825. 1972. 30

M. Yuzawa. A method of obtaining the oblique incident sound absorption coefficient through an on-the-spot measurement. Appl. Acoust. 8, 1975. 31

A. Krokstad, S. Strom, S. Sorsdal. Calculation of the acoustical room response by the use of a ray tracing technique. J. Sound Vib. (1968) 8 (1), 118-125 32

M. Hodgson. Case history: Factory noise prediction using ray tracing - Experimental validation and the effectiveness of noise control measures. Noise Con. Eng. J. 33(3) 1989 97-104.

106 33

M. R. Hodgson. Sound-propagation curves in industrial workrooms: statistical trends and empirical prediction models. J. Building Acoust. 3(1), 25-32, 1996. 34

H. Kuttruff. Sound propagation in working environments. Proc. 5th FASE symposium, Thessaloniki, 17-32, 1985. 35

J. K. Thompson, L. D. Mitchell, C. J. Hurst. A modified room acoustics approach to determine sound-pressure levels in irregularly-proportioned workroom spaces. Proc. Inter-Noise '76, 465-468, 1976. 36

F. Cotana. An improved room acoustics model. Appl. Acoust. 61, 1-25, 2000.

37

M. Hodgson. On the accuracy of models for predicting sound propagation in fitted rooms. J. Acoust. Soc. Am. 88(2) 1990 871-878. 38

M. Hodgson. Experimental evaluation of simplified models for predicting noise levels in industrial workrooms. J. Acoust. Soc. Am. 103(4) 1998 1933-1939. 39

S. Kirkup. The boundary element method in acoustics. Integrated sound software, Todmorden, UK, 1998. 40

C. A. Brebbia, M. H. Aliabadi. Adaptive finite and boundary element methods. CMP, Southampton, UK, 1993. 41

S. N. Y. Gerges, A. J. Calza. Acoustic barriers: Analytical methods, boundary element method and experimental verification. Build. Acoust. 9(3) 2002 167-190. 42

T. Funkhouser, N. Tsingos, I. Carlbom, G. Elko, M. Sondhi, J. West, G. Pingali. P. Min, A.Ngan. A beam tracing method for interactive architectural acoustics. J. Acoust. Soc. Am. 115(2), 739-756, Frebrury 2004. 43

B. I. Dalenbäck, M. Kleiner, P. Svensson. A macroscopic view of diffuse reflection. J. Audio Eng. Soc. Vol. 42, No. 10, 1994 October 44

K. H. Kuttruff. Auralization of impulse responses modeled on the basis of ray-tracing results. J. Audio Eng. Soc., Vol. 41, No. 11, November 1993. 45

H. Lehnert, J. Blauert. Principles of binaural room simulation. Appl. Acoust. 36, 259291, 1992. 46

M. Kleiner, B-I. Dalenbäck, P. Svensson. Auralization - An overview. J. Audio Eng. Soc., Vol. 41, No. 11, November 1993. 47

J. J. Embrechts. Broad spectrum diffusion model for room acoustics ray-tracing algorithms. J. Acoust. Soc. Am. 107(4) 2000 2068-2081 48

V. Hongisto, J. Keränen, P. Larm. Simple model for the acoustical design of openplan offices. Acustica - acta acustica. Vol. 90 (2004) 481-495. 49

M. Vorländer. International Round Robin on room acoustical computer simulations. Proc.of 15th ICA, Trondheim, Norway 26-30 June 1995 689-692.

107 50

I. Bork. A comparison of room simulation software - The 2nd round robin on room acoustical computer simulation. Acustica - Acta acustica. Vol. 86(6), 943-956, 2000. 51

J. H. Rindel. Computer simulation techniques for acoustical design of rooms - How to treat reflections in sound field simulation. ASVA 97, Tokyo 2-4 April 1997. Proceedings p.201-208 52

J. S. Bradley, G. A. Soulodre. The acoustics of concert halls. Pysics World, May 1997. 53

U.M.Stephenson. Quantized pyramidal beam tracing - a new algorithm for room acoustics and noise inmission prognosis. Acustica - Acta Acustica. Vol 82. 1996. 517525. 54

J. Borish. Extension of the image model to arbitrary polyhedra. J. Acoust. Soc. Am., 75(6):1827-1836, 1984. 55

R. N. S. Hammad. Simulation of noise distribution in rectangular rooms by means of computer modelling techniques. Appl. Acoust. 24, 211-28, 1988. 56

M. Hodgson. On the prediction of sound fields in large empty rooms. J. Acoust. Soc. Am. 84(1) 1988 253-261. 57

M. Hodgson. Evidence of diffuse surface reflections in rooms. J. Acoust. Soc. Am. 89(2) 1991 765-771. 58

S. M. Dance. The representation of fittings in the computer prediction of sound propagation in enclosed spaces. Appl. Acoust. 52(1), 71-83, 1997. 59

S. M. Dance, B. M. Shield. Modelling of sound fields in enclosed spaces with absorbent room surfaces. Part I: performance spaces. Appl. Acoust. 58, 1-18, 1999. 60

S. M. Dance. Minimal input models for sound level prediction in fitted enclosed spaces. Appl. Acoust. 63, 359-372, 2002. 61

U. M. Stephenson. Quantized pyramidal beam tracing or a sound-particle-radiosityalgorithm?. Not published. 62

A. Farina. RAMSETE - a new pyramid tracer for medium and large scale acoustic problems. In Proc. Euro-Noise 95. Lyon, France 21-23 march 1995. 63

U.P. Svensson. Modelling acoustic spaces for audio virtual reality. Proc.1st IEEE Benelux Workshop on Model based Processing and Coding of Audio (MPCA-2002), Leuven, Belgium, November 15, 2002 64

M. Hodgson. On the accuracy of models for predicting sound propagation in fitted rooms. J. Acoust. Soc. Am. 88(2) 1990 871-878. 65

M. Vorländer, E. Mommertz. Definition and measurement of random-incidence scattering coefficients. Appl. Acoust. 60, 187-199, 2000.

108 66

V. Hongisto, J. Keränen, E. Airo, P. Olkinuora. Akustinen mallintaminen meluntorjuntasuunnittelussa. Työ ja ihminen tutkimusraportti 20 (in finnish). Helsinki, Finland 2001. 67

M. Lisa, J. H. Rindel. Predicting the acoustics of ancient open-air theaters. The importance of calculation methods and geometrical details. In Proc. Baltic-Nordic Acoustical Meeting (BNAM 2004). Paper No. 19. Marieham, Finland, 8-10 June, 2004. 68

E. Mommertz. Determination of scattering coefficients from the reflection directivity of architechtural surfaces. Appl. Acoust. 60, 201-203, 2000. 69

J. Y. Jeon, S. C. Lee, M. Vorländer. Development of scattering surfaces for concert halls. Appl.Acoust. 65, 341-355, 2004. 70

A. Farina. A new method for measuring the scattering coefficient and the diffusion coefficient of panels. Acustica - Acta acustica. Vol. 86(6), 928-942, 2000 71

C. H. Haan, K. W. Kwon. A method of evaluating surface diffusivity of rooms - I: computer modeling and field measurement. Appl. Acoust. 62, 1313-1327, 2001. 72

B-I. L.Dalenbäck. Room acoustic prediction based on a unified treatment of diffuse and specular reflection. J. Acoust. Soc. Am. 100 (2), Pt.1, August 1996 73

C. Lynge. ODEON Room acoustics program v. 2.5 user manual, industrial, auditorium and combined editions. Technical university of Denmark, Lyngby, Denmark, 1997. 74

G. Naylor: ODEON - Another hybrid room acoustical model. Appl. Acoust. 38, 131143, 1993. 75

G. Naylor. Treatment of early and late reflections in a hybrid computer model for room acoustics. Presented at the 124th Meeting of the Acoustical Society of America. J. Acoust. Soc. Am. (Abstracts), vol. 92, pt. 2, p. 2345 (1992 Oct.), paper 3aAA2. 76

J. H. Rindel. Modeling the directional characteristics of sound reflections. In Proc. Baltic-Nordic Acoustical Meeting (BNAM 2004). Paper No. 25. Marieham, Finland, 8-10 June, 2004. 77

U. P. Svensson. Modeling room acoustics. In Proc. Baltic-Nordic Acoustical Meeting (BNAM 2004). Paper No. i03. Marieham, Finland, 8-10 June, 2004. 78

J. H. Rindel. Transmission of traffic noise through windows - Influence of incident angle on sound insulation theory and experiment. Report No. 9, Laboratoriet for Akustik, DTH, Lyngby, 1975. 79

J.S.Bradley. Some further investigations of the seat dip effect. J. Acoust. Soc. Am. 90 (1), July 1991 80

W. J. Davies, Y. W. Lam. New attributes of seat dip attenuation. Appl.Acoust. 41, 123, 1994.

109 81

T. J. Schultz, B. G. Watters. Propagation of sound across audience seating. J. Acoust. Soc. Am. 36, 885-896 (1964). 82

G. M. Sessler, J. E. West. Sound transmission over theatre seats. J. Acoust. Soc. Am. 36, 1725-1732 (1964). 83

J. Keller. Geometrical theory of diffraction. J. Acoust. Soc. Am. 52(2):116-130, 1962.

84

G. Benedetto, R. Spagnolo. A study of barriers in enclosures by a ray-tracing computer model. Appl. Acoust. 17, 183-199, 1984. 85

U. J. Kurze. Scattering of sound in industrial spaces. J. Sound. Vib. 98(3) 1985 349364. 86

N. Tsingos, I. Carlbom. Geometrical theory of diffraction for modeling acoustics in virtual environments. Proc. SIGGRAPH-2001. 12-17 August 2001. Los Angeles, CA, USA. 87

E. van Haaren and P. H. van Tol. Validation of ray acoustics applied for the modelling of noise barriers. J. Sound Vib. 231(3) 2000 681-688 88

U. P. Svensson, R. I. Fred, J. Vanderkooy. An analytical secondary source model of edge diffraction impulse responses. J. Acoust. Soc. Am. 106 (5) November 1999 (23312344). 89

S. R. Shankland. Acoustics of Greek theatres. Physics Today (1973) 30-35.

90

ISO 354. Acoustics - Measurement of sound absorption in a reverberation room. International Organization for Standardization, Genève Switzerland, 2003. 91

H. Lehnert. Systematic errors of the ray-tracing algorithm. Appl. Acoust, 38, 207221, 1993. 92

Z. Xiangyang, C. Kean, S. Jincai. On the accuracy of the ray-tracing algorithms based on various sound receiver models. Appl. Acoust. 64, 433-441, 2003. 94

http://www.ptb.de/en/org/1/17/173/roundrob2_1.htm

95

I. Brohnshtein, K. Semendiaev. Manual de Matemáticas. URSS, Moscú. 1988.

96

M. R. Schroeder. New method of measuring reverberation time. J. Acoust. Soc. Am. 37 1965 409-12. 97

D. Oliva, J. Keränen, V. Hongisto. New acoustical model for the workplace design. In Proc. Akustiikkapäivät 2003. p. 103-108. 6-7 October 2003, Turku, Finland. 98

D. Oliva, J. Keräanen. P. Larm, V. Hongisto. Virtual space 4D: The room acoustical tool. In Proc. Baltic-Nordic Acoustical Meeting (BNAM 2004). Paper No. 13. Marieham, Finland, 8-10 June, 2004.

110 99

ISO 3741. Acoustics - Determination of sound power levels of noise sources Prediction methods for broad-band sources in reverberation rooms. International Organization for Standardization, Genève Switzerland, 1988. 100

ISO 9613-1. Acoustics - Attenuation of sound during propagation outdoors. International Organization for Standardization, Genève Switzerland, 1993. 101

ISO 3746. Acoustics - Determination of sound power levels of noise sources using sound pressure - Survey method using an enveloping measurement surface over a reflecting plane. International Organization for Standardization, Genève Switzerland, 1995. 102

ISO 3747. Acoustics - Determination of sound power levels of noise sources Survey method using a reference sound source. International Organization for Standardization, Genève Switzerland, 1987. 103

F. P. Metchel. The corner source model. Acustica - acta acustica 86(6).2000.957-969

A1-1

Appendix 1 The absorption and scattering properties of materials are needed in room acoustics modeling. Standard methods are available and they are briefly introduced in the following. Measurement of sound absorption in a reverberation room, ISO 354A1 When a sound source operates in an enclosed space, the level to which reverberant sound builds up, and the subsequent decay of reverberant sound when the source is stopped, are governed by the sound-absorbing characteristics of the boundary surfaces and objects within the space. In general, the fraction of the incident sound power absorbed at a surface depends upon the angle of incidence. In order to relate the reverberation time of an auditorium, office, workshop, etc. to the noise reduction that would be effected by an absorbing treatment, a knowledge of the sound-absorbing characteristics of the surfaces, usually in the form of a suitable average over all angles of incidence, is required. Since the distribution of sound waves in typical enclosures includes a wide and largely unpredictable range of angles, it is convenient, for the purposes of standardization, to take a uniform distribution as the basic condition. If, furthermore, the sound intensity is independent of location within the room, such a distribution is called a diffuse sound field, and the sounds reaching a room surface are said to be at random incidenceA1. Reverberation time is the time that would be required for the sound pressure level to decrease by 60 dB after the sound source has stopped. The quantity is denoted by T and is expressed in seconds. The calculation of the absorption coefficient of a certain material, A [m2], is based on the difference between the reverberation time measured in a reverberation room with a sample of the material, T2, and without it, T1. The equivalent sound absorption area of the test specimen is calculated by the following formula

A = 55.3

V c

⎛1 1⎞ ⎜⎜ − ⎟⎟ ⎝ T2 T1 ⎠

(A.1)

where V [m3] is the volume of the empty reverberation room, c [m/s] is the speed of sound in air, T1 [s] is the reverberation time of the empty reverberation room.

A1-2

The reverberation time of the room in each frequency band is expressed by the arithmetic mean of the total number of reverberation time measurements made in that frequency band, for example, two each of six sound source/microphone combinations. Finally, the sound absorption coefficient α of the test specimen is calculated using the formula

α=

A S

(A.2)

where A [m2] is the equivalent sound absorption area calculated by Eq. (A.1), and S is the area [m2] of the test specimen. A photo of the test arrangement is shown in Figure A1, when the absorption coefficient of a material present in Room 3 was measured in the laboratory of acoustics of the Finnish Institute of Occupational Health, see sections 4.1.3 and 4.2.3. Further information about the test arrangement, i.e. size and diffusion characteristics of the reverberation room, size of the sample, and receiving equipment, is given in more detail in ISO 354A1.

Measurement of the random-incidence scattering coefficient of surfaces, ISO/DIS 17497-1A2 The principle of the measuring technique takes advantage of the fact that an impulse response, measured over a diffuser in the same place but in different orientations, presents differences in amplitude and phase in its late part. When impulse responses measured for different orientations of the diffuser are averaged, the late non-correlated parts will cancel each other out, resulting in a specular impulse response, whose decay is smaller when compared to a single measurement, see Figure A2A3.

A1-3

Figure A2: Impulses responses measured in a scale model reverberation chamber. After Vorländer et alA3. The scattering coefficient is defined as the ratio of the scattered energy to the total reflected energy. In Figure 3, the energy components of the reflected sound are illustrated. From the total incident sound, one part α , is absorbed by the material. The rest is partly scattered and partly specularly reflected. This definition agrees very well with the model of diffuse reflections used in our raytracing program.

Figure A3: Energy components of the reflected sound. The method goes as follows. Impulse responses are measured without and with the sample following ISO 354 and giving the reverberation times T1 and T2, respectively. The random incidence absorption coefficient α of the material is then calculated by Eq. (A.1). Then, a rotating

A1-4

base plate is placed in the reverberation room, and T3 is measured. Finally, the sample of the material is placed on the rotating base, and n measurements of the reverberation time are performed, with the sample rotated between each measurement by ∆ϕ = 360° / n . A value of n = 72 is preferred, corresponding to angular steps of ∆ϕ = 5° . The mean value of all these measurements is T4. The specular absorption coefficient, α SPEC , is calculated by Eq. (A.3)

V⎛ 1 1 ⎞ ⎟ − S ⎝ c 4T4 c3T3 ⎟⎠

α SPEC = 55.3 ⎜⎜

(A.3)

Under ideal conditions, the reverberation times T1 and T3 should be equal. The random-incidence scattering coefficient s is finally calculated using the formula

s=

α SPEC − α 1−α

(A.4)

More details about the test arrangement are given in ISO/DIS 17497A2.

Bibliography: A1

ISO 354. Acoustics - Measurement of sound absorption in a reverberation room. International

Organization for Standardization, Genève Switzerland, 2003. A2

ISO/DIS 17497-1 Acoustics - Measurement of the sound scattering properties of surfaces.

International Organization for Standardization, Genève Switzerland, 2002 A3

M. Vorländer, J-J. Embrechts, L. Geetere, G. Vermeir, M. Averlar Gomes. Case studies in

measurements of random incidence scattering coefficient. Acustica - acta acustica. Vol. 90 (2004) 858-867

A2-1

Appendix 2 Survey of existing commercial software A survey of existing commercial acoustic software has been carried out, in order to find out some more about the current state of the art in room acoustics modeling. A total of 18 programs have been considered. In Table 2.2, their name, price, internet address and home origin are listed. In Table 2.3, their properties are summarized. It is important to notice the large range of prices. The software referred to in Tables 2.2 and 2.3 by the number 17, dBWorks, is the acoustic model developed by the author of this thesis. Table 2.2: General information on commercial acoustic software. Name

Price

Net address

1

ODEON 6.5

5080-14400 € http://www.dat.dtu.dk/~odeon

Intotel, Finland

2

CATT v8

?

http://www.catt.se

Catt, Sweden

3

RayNoise

24 650 €

http://www.lmsintl.com

LMS, Belgium

4

ULYSSES

995 €

http://www.ifbsoft.de

IFBsoft, Germany

5

SoundPlan

6 000 €

http://www.soundplan.com

SP, Sweden

6

PlantNoise

2100 $

http://www.soeh.ubc.ca/hodgson_research/plantnoise.htm UBC, Canada

7

ClassTalk

652 $

http://www.soeh.ubc.ca/hodgson_research/classtalk.htm

UBC, Canada

8

RAMSETE 2 3 750 €

http://www.ramsete.com

Parma YO, Italy

9

Ease 4.0

?

http://www.olsonsound.com/EASE/index.html

Olson Sound Design

10 VisualEars

49 €

http://www.kbacoustics.com/visualears/

KB Acoustics

11 CARA 2.2

49-75 €

http://www.cara.de/

Scanteknik, Finland

12 NEMPEE

155 €

http://www.nempee.com/AcouSoft.htm

Nempee, India

13 Acoustic-X

286 €

http://www.pilchner-schoustal.com/old/acoustic-x

Pilchner Schoustal

14 Acousalle

245 €

http://lesowww.epfl.ch/anglais/Leso_a_frame_sof.html

LESO-PB, Switzerland

15 Sopran

Not in market …

16 RaySoft

?

http://www.oe.com.tw/01/prode2.htm

17 dBWorks

in progress

…

18 MagiCAD

2600-3500€

http://www.progman.fi

01dB, France

Progman, Finland

A2-2 Table 2.3: General information on commercial acoustic software. Commercial name. Applied method. Application field. Calculated acoustic parameters. Name

Method

Application

Acuoustic parameters

1

ODEON 6.5 Image + Ray-tracing

Auditorium, industrial halls SPL, RT, EDT, STI,…, Auralization

2

CATT v8

Image + Cone-tracing

Auditorium, industrial halls SPL, RT, EDT, STI,…, Auralization

3

RayNoise

Image + Cone-tracing

Auditorium, industrial halls SPL, RT, STI,…, Auralization

4

Ulysses

Image + Ray-tracing

Loudspeaker arrangement

SPL, RT, STI,…, Auralization

5

SoundPlan

Image + Sabine diffuse field

Industrial spaces

SPL, RT

6

PlantNoise

Empirical (Heerema & Hodgson)

Industrial spaces

SPL, RT, Auralization

7

ClassTalk

Empirical (Heerema & Hodgson)

Speech and classrooms

RT, STI, Auralization

8

RAMSETE

Pyramid-tracing

Industrial spaces

SPL, RT, Auralization

9

Ease 4.0

Image + Ray-tracing

Loudspeaker arrangement

SPL, RT, EDT, RASTI, Auralization

10 VisualEars

Mode analysis

Listening rooms

SPL

11 CARA 2.2

Mode analysis + diffuse field

Listening rooms

SPL, RT, EDT, Auralization

12 NEMPEE

Diffuse field theory

General spaces

RT

13 Acoustic-X

Mode analysis + Ray-tracing

Listening rooms

SPL, RT

14 Acousalle

Mode analysis + diffuse field

Auditorium, industrial halls RT

15 Sopran

Particle-tracing

Auditorium, industrial halls SPL, RT, EDT,…

16 RaySoft

Ray-tracing

Industrial spaces

17 dBWorks

Image + Ray-tracing

Auditorium, industrial halls SPL, RT, EDT, STI,…, Auralization

18 MagiCAD

Unknown (Sabine?)

Ventilation noise

SPL, RT

SPL