Determination of the acoustical diffusion and impedance characteristics of surfaces by Bayesian inversion of a modified Helmholtz equation by Gavin Arthur Mervyn Wyverstone Steininger B.Sc., Simon Fraser University, 2006 A THESIS SUBMITTED IN THE PARTIAL FULFILMENT OF THE REQUIRMENTS OF THE DEGREE OF Master of Applied Science in The Faculty of Graduate Studies (Mechanical Engineering) The University of British Columbia (Vancouver) August 2009 © Gavin Arthur Mervyn Wyverstone Steininger, 2009 Abstract A surface diffusion coefficient is used in architectural-acoustics to evaluate the effectiveness of diffusing surfaces. The inclusion of the diffusion characteristics is also important for the accuracy of room prediction models. Another important parameter is the absorption or impedance of a surface. In settings with significant low-frequency noise, phase effects are important; consequently impedance values of surfaces are necessary for accurate modeling. A review of existing models for specular and diffuse reflection is made. A new diffusion coefficient is defined and included in three new forward models for predicting the steady-state sound-pressure level above a finiteimpedance plane in an otherwise free field. Data are collected for several typical architectural surfaces in an anechoic chamber. Inverse methods are utilized in order to estimate the diffusion coefficient for surfaces given each of the models. This is done without knowledge of the surface impedance, which is simultaneously estimated. The models are compared with each other and with independently measured values of the surface impedance and diffusion. Inversion is found to be a reasonable way of determining the diffusion properties of a surface. ii Contents Abstract.............................................................................................................................. ii Contents ............................................................................................................................ iii List of Tables ..................................................................................................................... v List of Figures................................................................................................................. viii List of Symbols ................................................................................................................. xi Acknowledgments ........................................................................................................... xii 1 Introduction.....................................................................................................................1 1.1 Objective and motivation.......................................................................................... 1 1.2 Outline and methodology.......................................................................................... 2 2 Background concepts and theory ................................................................................. 3 2.1 Fundamental acoustical relations .............................................................................. 3 2.2 Sound intensity and impedance ................................................................................ 4 2.3 Planar waves ............................................................................................................. 5 2.4 Spherical waves ........................................................................................................ 6 2.6 Diffuse reflections................................................................................................... 10 2.7 Inverse theory.......................................................................................................... 16 3 Literature review ......................................................................................................... 24 3.1 Inverse methods ...................................................................................................... 24 3.2 Direct surface-acoustic measures............................................................................ 28 4 Experimental setup, loudspeaker characteristics and data ..................................... 33 4.1 Experimental setup.................................................................................................. 33 4.2 Loudspeaker characteristics .................................................................................... 35 4.3 Data description ...................................................................................................... 40 5 Forward models and inversion procedure................................................................. 46 5.2 The diffuse-reflection energy model....................................................................... 51 5.3 Inversion techniques ............................................................................................... 53 6 Results and Discussion................................................................................................. 59 6.1 Image-surface-model results................................................................................... 59 6.2 Image-surface semi-phase model results ................................................................ 64 iii 6.3 Diffuse-reflecting-energy model results ................................................................. 69 6.4 Discussion ............................................................................................................... 70 7 Conclusions and suggested further work................................................................... 75 7.1 Conclusions............................................................................................................. 75 7.2 Further work............................................................................................................ 76 Bibliography .................................................................................................................... 78 Appendices....................................................................................................................... 80 Appendix A ...................................................................................................................80 Appendix B .................................................................................................................. 84 Appendix C .................................................................................................................. 88 Appendix D .................................................................................................................. 91 Appendix E .................................................................................................................. 92 iv List of Tables 4.1: The loudspeaker amplitude for the six pure-tone frequencies....................................34 4.2: The loudspeaker amplitude for the white noise in octave bands. ...............................35 4.3: Parameter values of the directivity function at each frequency..................................39 6.1: The measured and estimated impedance values (Z = Z re + iZ im ) for the Flat and Soft surfaces. The measured values are taken from Graves and Hodgson [18]. .......................71 A.1: Estimates (Est) and bounds (LB, UB) for DREM with BB data. ..............................80 A.2: Estimates (Est) and bounds (LB, UB) for DREM with flat data. ..............................80 A.3: Estimates (Est) and bounds (LB, UB) for DREM with FW data...............................81 A.4: Estimates (Est) and bounds (LB, UB) for DREM with SB data................................81 A.5: Estimates (Est) and bounds (LB, UB) for DREM with Sine1 data. ..........................81 A.6: Estimates (Est) and bounds (LB, UB) for DREM with Sine2 data. ..........................81 A.7: Estimates (Est) and bounds (LB, UB) for DREM with BB data. ..............................82 A.8: Estimates (Est) and bounds (LB, UB) for DREM with 250 Hz data.........................82 A.9: Estimates (Est) and bounds (LB, UB) for DREM with 500 Hz data.........................82 A.10: Estimates (Est) and bounds (LB, UB) for DREM with 1000 Hz data.....................83 A.11: Estimates (Est) and bounds (LB, UB) for DREM with 2000 Hz data.....................83 A.12: Estimates (Est) and bounds (LB, UB) for DREM with 4000 Hz data.....................83 A.13: Estimates (Est) and bounds (LB, UB) for DREM with 8000 Hz data.....................83 B.1: Estimates (Est) and bounds (LB, UB) for DREM with BB data. ..............................84 B.2: Estimates (Est) and bounds (LB, UB) for DREM with flat data. ..............................84 B.3: Estimates (Est) and bounds (LB, UB) for DREM with FW data...............................85 B.4: Estimates (Est) and bounds (LB, UB) for DREM with SB data................................85 v B.5: Estimates (Est) and bounds (LB, UB) for DREM with Sine1 data............................85 B.6: Estimates (Est) and bounds (LB, UB) for DREM with Sine2 data............................85 B.7: Estimates (Est) and bounds (LB, UB) for DREM with BB data. ..............................86 B.8: Estimates (Est) and bounds (LB, UB) for DREM with 250 Hz data. ........................86 B.9: Estimates (Est) and bounds (LB, UB) for DREM with 500 Hz data. ........................86 B.10: Estimates (Est) and bounds (LB, UB) for DREM with 1000 Hz data. ....................87 B.11: Estimates (Est) and bounds (LB, UB) for DREM with 2000 Hz data. ....................87 B.12: Estimates (Est) and bounds (LB, UB) for DREM with 4000 Hz data. ....................87 B.13: Estimates (Est) and bounds (LB, UB) for DREM with 8000 Hz data. ....................87 C.1: Estimates (Est) and bounds (LB, UB) for DREM with BB data. ..............................88 C.2: Estimates (Est) and bounds (LB, UB) for DREM with flat data. ..............................88 C.3: Estimates (Est) and bounds (LB, UB) for DREM with FW data...............................88 C.4: Estimates (Est) and bounds (LB, UB) for DREM with SB data................................88 C.5: Estimates (EST) and bounds (LB, UB) for DREM with Sine1 data..........................88 C.6: Estimates (Est) and bounds (LB, UB) for DREM with Sine2 data............................88 C.7: Estimates (Est) and bounds (LB, UB) for DREM with Soft data. .............................89 C.8: Estimates (Est) and bounds (LB, UB) for DREM with 250 Hz data. ........................89 C.9: Estimates (Est) and bounds (LB, UB) for DREM with 500 Hz data. ........................89 C.10: Estimates (Est) and bounds (LB, UB) for DREM with 1000 Hz data. ....................89 C.11: Estimates (Est) and bounds (LB, UB) for DREM with 2000 Hz data. ....................89 C.12: Estimates (Est) and bounds (LB, UB) for DREM with 4000 Hz data. ....................90 C.13: Estimates (Est) and bounds (LB, UB) for DREM with 8000 Hz data. ....................90 vi D.1: Displays the estimates (Est) and bounds (LB, UB) of the average absorption of the surfaces for the ISM, ISSPM and DREM models. The values measured by Bibby [23] are also presented.....................................................................................................................91 E.1: Estimates (Est) and bounds (LB, UB) for ISM with BB broad-band data.................92 E.2: Estimates (Est) and bounds (LB, UB) for ISM with flat broad-band data.................92 E.3: Estimates (Est) and bounds (LB, UB) for ISM with FW broad-band data. ...............93 E.4: Estimates (Est) and bounds (LB, UB) for ISM with SB broad-band data. ................93 E.5: Estimates (Est) and bounds (LB, UB) for ISM with Sine1 broad-band data. ............93 E.6: Estimates (Est) and bounds (LB, UB) for ISM with Sine2 broad-band data. ............93 E.7: Estimates (Est) and bounds (LB, UB) for ISM with BB broad-band data.................94 E.8: Estimates (Est) and bounds (LB, UB) for ISSPM with BB broad-band data. ...........94 E.9: Estimates (Est) and bounds (LB, UB) for ISSPM with flat broad-band data. ...........94 E.10: Estimates (Est) and bounds (LB, UB) for ISSPM with FW broad-band data..........95 E.11: Estimates (Est) and bounds (LB, UB) for ISSPM with SB broad-band data...........95 E.12: Estimates (Est) and bounds (LB, UB) for ISSPM with Sine1 broad-band data.......95 E.13: Estimates (Est) and bounds (LB, UB) for ISSPM with Sine2 broad-band data.......95 E.14: Estimates (Est) and bounds (LB, UB) for ISSPM with BB broad-band data. .........96 vii List of Figures 2.1: Illustration of the principle of the image-source model................................................7 2.2: Shows a sound beam emitted at angle θ r from an energized region , of the surface. A semi-spherical shell centered on the energized region is also depicted.........................11 2.3: Illustration of a sound ray diffusely reflecting from a sub-surface ●, with random impedance and roughness. The angles displayed are those in Eq. (2.6.3).........................12 3.1: Diagram of a typical impedance tube with source on the right-hand side and sample on the left-hand side. .........................................................................................29 3.2: Illustrates a typical spherical-decoupling measurement setup. The diagram has been augmented to include the image receivers , which exist only conceptually....................30 4.1: A plan-view representation of the experimental setup. ..............................................34 4.2: The measured ○, and predicted , pure-tone sound-pressure levels in the empty anechoic chamber assuming no loudspeaker directivity, plotted against the direct distance from the source to the measurement location. ...................................................................36 4.3: The measured ○, and predicted , white-noise sound-pressure levels in the empty anechoic chamber assuming no loudspeaker directivity, plotted against the direct distance from the source to the measurement location. ...................................................................37 4.4: The measured ○, and smoothed , loudspeaker directivities (Delta) plotted against the polar angle (theta) in radians........................................................................................38 4.5: The measured ○, and smoothed , loudspeaker directivities (Delta) plotted against the azimuth angle (Phi) in radians. ....................................................................................40 4.6: The experimental setup while the SB surface was being evaluated (the longer obstacles were not present at the time of measurement)....................................................41 4.7: The measured steady-state sound-pressure levels (pure-tone excitation) above the surfaces ○f interest, plotted against the x (Track 1) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 )...........................42 4.8: The measured steady-state sound-pressure levels (pure-tone excitation) above the surfaces ○f interest, plotted against the y (Track 2) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 .)..........................43 viii 4.9: The measured steady-state sound-pressure levels (broad band excitation) above the surfaces ○f interest, plotted against the x (Track 1) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 )...........................44 4.10: The measured steady-state sound-pressure levels (broad band excitation) above the surfaces ○f interest, plotted against the y (Track 2) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 )...........................45 5.1: The path of a sound ray ( ) diffusely reflecting off a conceptual surface (···) which is a rotation of the actual surface ( ). The angles displayed are those in Eq. (5.1.9) through Eq. (5.1.14). ..........................................................................................................49 5.2: An illustration of an image surface ( ) generated by a source and receiver pair above a reflecting surface ( ). The image surface is projected into the xz plane, the yz plane and the xy plane...................................................................................................51 6.1: The parameter estimates and confidence intervals for the ISM for all seven surface samples at the six frequencies (250 Hz, 500 Hz, 1000 Hz, 2000 Hz, 4000 Hz, 8000 Hz) .... ............................................................................................................................................60 6.2: Normalized histograms of the marginal posterior distributions of the parameters for the Flat surface...................................................................................................................61 6.3: Normalized histograms of the marginal posterior distributions of the parameters for the FW surface. ..................................................................................................................61 6.4: Normalized histograms of the marginal posterior distributions of the parameters for the SB surface. ...................................................................................................................62 6.5: Normalized histograms of the marginal posterior distributions of the parameters for the BB surface....................................................................................................................62 6.6: Normalized histograms of the marginal posterior distributions of the parameters for the Sine1 surface. ...............................................................................................................63 6.7: Normalized histograms of the marginal posterior distributions of the parameters for the Sine2 surface. ...............................................................................................................63 6.8: Normalized histograms of the marginal posterior distributions of the parameters for the Soft surface. .................................................................................................................64 6.9: The parameter estimates and confidence intervals for the ISSPM for all seven surface samples at the six frequencies (250 Hz, 500 Hz, 1000 Hz, 2000 Hz, 4000 Hz, 8000 Hz).65 6.10: Normalized histograms of the marginal posterior distributions of the parameters for the Flat surface...................................................................................................................66 ix 6.11: Normalized histograms of the marginal posterior distributions of the parameters for the FW surface. ..................................................................................................................66 6.12: Normalized histograms of the marginal posterior distributions of the parameters for the SB surface. ...................................................................................................................67 6.13: Normalized histograms of the marginal posterior distributions of the parameters for the BB surface....................................................................................................................67 6.14: Normalized histograms of the marginal posterior distributions of the parameters for the Sine1 surface. ...............................................................................................................68 6.15: Normalized histograms of the marginal posterior distributions of the parameters for the Sine2 surface. ...............................................................................................................68 6.16: Normalized histograms of the marginal posterior distributions of the parameters for the Soft surface. .................................................................................................................69 6.17: The parameter estimates and confidence intervals for the DREM for all seven surface samples at the six frequencies (250 Hz, 500 Hz, 1000 Hz, 2000 Hz, 4000 Hz, 8000 Hz).............................................................................................................................70 6.18: The measured values and estimates with confidence intervals for the acousticalabsorption coeffiecent (Alpha) of seven surfaces. Measured value , ISM ,ISSPM and DREM .....................................................................................................................73 x List of Symbols x t r y z ω k c ρ p u Lp I Z1 , Z 2 & Z • Z Z Re Z Im β α α RP RE Q F θ θ φ cartesian coordinates arbitrary point in time cartesian distance angular frequency wave number speed of sound density of a medium sound pressure particle velocity sound pressure level sound intensity characteristic impedance of a medium specific impedance of a medium real coefficients of the specific impedance of a medium imaginary coefficients of the specific impedance of a medium admittance of a medium absorption coefficient random incidence absorption coefficient planar reflection confident planar energy reflection coefficient reflection factor bottom loss factor vector of unknown parameters polar angle azimuth angle xi Acknowledgments To Ms. Main and Ms. Tessman who where there long before the beginning. Dr Hodgson without whom the completion of this document would have been impossible for many, many reasons. Finally, Ms. Ellerton for typing most of it. xii Chapter 1 Chapter 1 Introduction 1.1 Objective and motivation This thesis discusses the determination of the acoustical parameters of room surfaces. These can in some cases be measured directly. In other cases direct measurement is not possible, and indirect methods are of interest. There is widespread use of inverse methods in the discipline of acoustics; however architectural acoustics appears to be underrepresented in their use. This is surprising, as architectural acoustics shares many similar problems with underwater and outdoor acoustics, as they all consider sound in somewhat bounded spaces. This research applied inverse-method techniques from non-architectural fields of acoustics, to facilitate the solution of architectural-acoustic problems. The specific problem of interest is the accurate determination of the acoustical-impedance and diffusion characteristics of typical architectural surfaces. This problem is of great importance in architectural acoustics, as accurate knowledge of the impedance and diffusion characteristics of surfaces is crucial in the architecturalacoustical design process. The inclusion of phase information—that is, complex impedance instead of energy absorption—is important in industrial settings where a significant portion of the background noise is at low frequency. It has been shown that failure to accurately include surface-diffusion characteristics in prediction models leads to inaccurate predictions of room reverberation times. The determination of the acoustical impedance and diffusion characteristics of a surface does not lend itself well to direct measurement. Conventional, direct methods for measuring the acoustical impedance of a surface assume wholly specular reflection. This is a significant problem; as such methods interpret non-specular reflection as increased absorption, producing inaccurate values. Conventional methods for measuring the diffusion characteristics of a surface assume knowledge of the absorption of the surface, and are currently based in diffuse-field theory and do not consider phase. 1 Chapter 1 1.2 Outline and methodology Chapter 2 is a tutorial section on the background concepts of acoustics and inverse methods. This is followed by a review of methods for determining the acoustical properties of surfaces in Chapter 3. Chapter 4 discusses the collection of data and Chapter 5, the implementation of inverse methods. Chapter 6 presents the independent research, which consists of three components: the development of new forward models, the collection of data, and inversion of the models. Three new forward models were developed using the concept of the image-source method, and have been modified to include diffuse reflection. They differ from each other in their inclusion of phase information. The data were collected in a hemi-anechoic chamber with seven different types of architectural surfaces used as the sole reflecting surface. The inversion of the models was done in such in a way as to not only gain point estimates, but also confidence intervals for the estimates. The models are inverted directly and not approximated. Finally, a comparison of the three forward models is made. The criteria for these comparisons were: the ability of a model to predict the observed data, the level of conformity of the estimated parameters to direct measurements, and the level of certainty of the estimates produced. 2 Chapter 2 Chapter 2 Background concepts and theory In order to have a discussion on diffuse sound reflection, it is first necessary to have a strong understanding of the physical principles of acoustic-wave theory. This will be the primary purpose of this chapter; however it will also introduce inverse theory and optimization. The acoustics material unless otherwise cited is from Kuttruff [1]. 2.1 Fundamental acoustical relations The most basic element of acoustics is the acoustical wave. Acoustical waves are longitudinal; that is, they displace particles of the propagation medium in the same direction as the wave propagates. Acoustical waves are also three dimensional. Acoustical waves can be characterized in several ways. First, a wave can be characterized by the displacement ζ of the particles of the medium affected by the wave. The wave can also be characterized by the velocity of the particles u . A third quantity, that characterizes the amplitude of an acoustical wave, is sound pressure p . This value is also commonly reported in decibels dB as sound-pressure level L p . Equation (2.1.1) gives the sound-pressure level in terms of the pressure, where p ∗ is the reference pressure of 20µPa . p L p = 20 log ∗ p (2.1.1) The wave equation is one of the most fundamental acoustical formulas. The wave equation Eq. (2.1.2) can be derived from three equations; these are the continuity equation, the Euler equation and the equation of state [2]. It is assumed that the fluid medium, in which the wave propagates, is homogeneous, isotropic and perfectly elastic. It should also be noted that the linearization approximations used in this derivation are only valid up to 110 dB. 3 Chapter 2 ∂p = c 2 ∇ 2 p or ∂t 2 ∂2 p ∂2 p ∂2 p 1 ∂2 p + + − =0 ∂ 2 x ∂ 2 y ∂ 2 z c 2 ∂ 2t (2.1.2) The Laplacian appears on account of the three-dimensional nature of acoustical waves. Wave fronts may take a variety of geometries; the simplest of these is a plane. The wave front may be non-planar, such as spherical. 2.2 Sound intensity and impedance Another important attribute of a sound field is the sound intensity. Sound intensity is defined as the average rate of flow of acoustical energy through a unit of area normal to the direction of sound propagation [2]. Eq. (2.2.1) gives intensity as the product of particle velocity and sound pressure. I = pu (2.2.1) The characteristic acoustical impedance Z • of a medium is defined as the ratio of sound pressure to the particle velocity, Eq. (2.2.2). For planar waves the impedance is also equal to the product of the density of the material and the speed of sound within it, Eq. (2.2.3). For spherical waves the relation is more complicated as sound pressure and particle velocity are not in phase. As shown in Eq. (2.2.4), if the radius is large then the impedance will be similar to that of a planar wave. This is expected as the curvature of a large sphere is less than that of a small sphere. However for small spheres the second term in the denominator will dominate and the impedance will be much smaller. The models described in Chapter 5 used the specific not characteristic impedance of a surface. The specific impedance Z of a surface is the ratio of the characteristic impedance of the surface to the characteristic impedance of the reference medium; the admittance β of a surface is the inverse of this ratio. 4 Chapter 2 p u (2.2.2) Z • = ρc (2.2.3) Z• = Z• = ρc (2.2.4) 1 1 + ikr 2.3 Planar waves Waves with planar wave fronts are called plane waves. The most important properties of plane waves are that acoustical pressure and particle velocity have the same phase and magnitude at all points on a plane normal to the direction of wave propagation (wavefront). Plane waves usually occur in vents, pipes or other long and narrow spaces. However a plane wave may be used to approximate a non-plane wave if the curvature of the approximated waves is low. This normally occurs if the source of the wave is in the far field, at a large distance from the source relative to the wave length. As noted, the pressure of a plane wave only changes in the direction of its flow. Without loss of generality this can be considered the x axis; consequently Eq. (2.1.1) can be rewritten as Eq. (2.3.1). The function, p( x, t ) , can be solved for the harmonic solution is given in Eq. (2.3.2). This can be interpreted as two waves: the first has amplitude A , and is propagating in the positive x direction, the second has amplitude B and is propagating in the negative x direction. The strategy of describing a sound field in terms of several waves will be repeated throughout this chapter, as this appears to be a conceptually simple description of acoustical phenomena. It should be noted, however, that the combination or separation of waves is arbitrary. ∂ 1 ∂ ∂ 1 ∂ − + p ( x, t ) = 0 ∂x c ∂t ∂x c ∂t { p( x, t ) = Re Ae i (ωt −kx ) + Be i (ωt + kx ) } (2.3.1) (2.3.2) 5 Chapter 2 2.4 Spherical waves A spherical wave is any wave such that the acoustical pressure and particle velocity have the same phase and magnitude at all points on the surface of a sphere centered on the sound source. Not surprisingly, spherical waves are more easily described using a spherical coordinate system. Eq. (2.4.1) gives the wave equation, in which the Laplacian operator is rewritten in spherical form. On account of the spherical-wave property that only radius, not angle, determines pressure, particle velocity, etc. the equation can be simplified to yield Eq. (2.4.2). Eq. (2.4.2) is of the same form as Eq. (2.3.1) and therefore has a solution of the same form, given in Eq. (2.4.4). This can be solved for p yielding Eq. (2.4.5). ∂ 2 (rp ) 2 ∂ (rp ) 1 ∂ (rp ) ∂ (rp ) 1 ∂ 2 (rp ) 1 ∂ 2 (rp ) ( ) + + sin + − 2 = 0 (2.4.1) θ r ∂r ∂θ r 2 sin 2 (θ ) ∂φ 2 ∂r 2 r 2 sin (θ ) ∂θ c ∂t 2 ∂ 2 (rp ) 1 ∂ 2 (rp ) − 2 =0 ∂r 2 c ∂t 2 ∂ 1 ∂ ∂ 1 ∂ ∗ p ( x, t ) = 0 where − + ∂x c ∂t ∂x c ∂t (2.4.2) p ∗ = rp (2.4.3) p ∗ ( x, t ) = Aei (ωt − kx ) + Bei (ωt + kx ) (2.4.4) Aei (ωt −kx ) Bei (ωt + kx ) p( x, t ) = Re + r r (2.4.5) The spherical-wave equation has a similar interpretation as the plane-wave equation. The first term represents a wave of magnitude A propagating in an outward radial direction; the second term represents a wave converging at the origin with magnitude B . The main difference is that the pressure amplitudes of spherical waves do not remain constant as the wave propagates; they attenuate at a rate of 1 r . 6 Chapter 2 2.5 Basic sound reflection Figure 2.1: Illustration of the principle of the image-source model. To conceptualize sound reflection it is convenient to consider an idealized situation, the one displayed in Fig. (2.1). The idealized situation is a single sound source within one, and above another, semi-infinite medium separated by a planar interface. The media have characteristic impedances Z1 and Z 2 , respectively. It is assumed that the sound field at the receiver is of the form of Eq. (2.5.1); that is, it is the combination of the contributions of the direct P1 and the reflected P2 paths. Ptot = P1 + Rp P2 (2.5.1) The quantity Rp in Eq. (2.5.1) is the reflection coefficient. If it is assumed that the acoustic waves are planar, then a solution for Rp can be found using the boundaryelement method. This is done by solving the system of equations Eq. (2.5.2) and Eq. (2.5.3) that are linear in terms of the amplitudes. It should be noted that Eq. (2.5.2) represents the continuity of sound pressure at the interface and that Eq. (2.5.3) requires that the normal particle velocity at the interface is zero. The y terms have been removed, as the equations describe interfaces at y = 0 . The angles θ i , θ r and θ T are the incident, 7 Chapter 2 reflected and transmitted angles. These are defined as the angle between the ray and the normal not the surface. Ai eik1x cos (θi ) + Ar eik1x cos (θ r ) = AT eik2 x cos (θT ) (2.5.2) Z Ai eik1x cos(θi ) − Ar e ik1x cos(θ r ) = 1 AT eik2 x cos(θT ) Z2 (2.5.3) The solution is given in Eq. (2.5.4). It can be helpful to express RP in terms of a product of terms for magnitude and phase ( χ ). If it is assumed that the interface is locally reactive then, by definition, Z 2 (θ i ) will be constant for all values of θ i ; as a consequence the solution simplifies to Eq. (2.5.5). Another important result of this process is that θ r must equal θ i ; that is, the reflection of a planar wave must be specular. Ar Z cos(θ i ) − Z 1 cos(θ T ) = Rp = 2 = R p e iχ Ai Z 2 cos(θ i ) + Z1 cos(θ T ) cos(θ i ) − Rp = cos(θ i ) + Z1 Z2 Z1 Z2 (2.5.4) (2.5.5) If the acoustical waves are assumed to be spherical, then a restricted solution for Rp was found by Sommerfeld [3] applying the boundary-element method to Eq. (2.5.6) and Eq. (2.5.7). The form of Sommerfeld’s solution is given in Eq. (2.5.8). (− ∇ 2 ) − k p = 4πδ (0,0, hs ) for Z ∂p = −ik 1 p ∂z Z2 z=0 for z>0 Ptot = P1 + QP2 = P1 + (Rp + (1 − Rp )F )P2 (2.5.6) (2.5.7) (2.5.8) 8 Chapter 2 For cases where both the source and the receiver are on the interface between the two media, Sommerfeld was able to express the total wave-field solution in the form of Eq. (2.5.9). The form of Eq. (2.5.9) allows for the interpretation of F as a surface-loss factor; that is, the additional attenuation of the sound waves as a result of the interface. An approximate solution for F was found, and demonstrated to be dependent on only two complex values, the first and second numeric distances ( ρ , τ ). This solution is used by Taraldson [3] to find the general solution given in Eq. (2.5.10). It should be noted that the Taraldson solution, while general and exact, is not explicit. P0 = 2 F ∞ F = 1− ∫ 0 e ikr r (2.5.9) e−w u u 2 1 − − ρ 2τ ∂w (2.5.10) 2 Z1 + cos(θ ) Z i ρ = kr2 2 2 Z 1 + 1 cos(θ ) Z2 Z τ = kr2 1 + cos(θ ) 2 Z2 i (2.5.11) (2.5.12) An efficient and reasonably accurate approximation was found by Lawhead and Rubnik [3] and is given in Eq. (2.5.13) and Eq. (2.5.14). This approximation is based on the observation that F has only a weak dependency on τ ; consequently it is removed, simplifying the equation. Although the erfc function looks computationally intensive because of its similarities to the cumulative density of a Gaussian distribution, there are many fast ways of calculating its value. Moreover it should be noted that this approximation is only good if kr2 , the wave number times the distance from the point of incidence to the listener position, is much larger than one. 9 Chapter 2 ( F ≈ 1 + i πρ e − ρ erfc − i ρ erfc( x ) = 2 π ∞ ∫e −x2 ) (2.5.13) ∂x (2.5.14) x So far the sound-pressure level of a point above a surface has only been considered with regards to its location—that is, its distance from the sound source along the direct and specularly-reflected paths. As for all pure tones, the sound-pressure level will have a sinusoidal variation with time. The amplitude of the sinusoid is given by Eq. (2.5.16). e ikr1 e ikr2 +Q = r1 r2 = = [r1 (Q Re cos(kr2 ) − Q Im sin (kr2 )) + r2 cos(kr1 )] [r2 (QIm cos(kr1 ) + QRe sin (kr2 )) + r2 sin (kr1 )] r1 r2 1 r1r2 +i (2.5.15) r1 r2 [r1 (QRe cos(kε r2 ) − QIm sin(kε r2 )) + r2 cos(kε r1 )]2 + [r1 (QIm cos(kε r2 ) + QRe sin (kε r2 )) + r2 sin(kε r1 )]2 (2.5.16) 2.6 Diffuse reflections So far the reflection paradigms that have been considered are derived from physical principles. However it has been shown that these reflection rules are inadequate for modeling real rooms [4]. A possible improvement is to allow for diffuse reflection; the reflected wave is distributed over a range of angles. In order to model diffuse reflection, it appears to be common (since Kuttruff [5]) to use a probability-density function (PDF). These PDFs are used in two similar but distinct ways. The first method is as follows: a sound ray that reflects from a surface is assigned a random reflecting angle θ r with probability of a given θ r dictated by the PDF. The alternative approach is to assume that energy is distributed in all directions, with the proportion of the energy in a given direction determined by the PDF. The former is computationally easier to implement; however the latter seems closer to reality, with the proviso that neither of the two methods is justifiable from physical principles. 10 Chapter 2 A common special case for diffuse reflection is Lambert’s Law reflection. Under Lambert’s Law it is assumed that the reflected sound intensity is the same in all directions. This does not imply that the proportion of energy reflected in all directions is equal or, in other terms, that the PDF for Lambert’s Law reflection is non-uniform. For the case where the reflected azimuth angle is assumed to be the same as that of the incident ray, this is the two-dimensional case. Then the PDF is given in Eq. (2.6.1). It should be noted that PDFs used in this report will be reported up to a multiplicative constant. f (θ R | θ i ) ∝ cos(θ R ) (2.6.1) Equation (2.6.1) can be derived by considering an arbitrary unit of surface that has been excited by a sound beam. If all of the sound energy is reflected into the same angle θ r then the width of the reflected beam will be cos(θ r ) times the width of the incident beam. This is shown in Fig. (2.2). Consequently the intensity of that beam will have increased; however Lambert’s Law reflection assumes that the intensity is the same for all angles of reflection, and Eq. (2.6.1) is used to normalize the intensity. The case where the reflecting sound ray is able to vary with the polar angle but not the azimuth angle is unsatisfactory, as the surface has already demonstrated that it reflects diffusely. In order to expand Eq. (2.6.1) to three-dimensional diffuse reflection, only two observations are necessary. The first is that the ray density does not vary with azimuth Figure 2.2: Shows a sound beam emitted at angle θ r from an energized region , of the surface. A semi-spherical shell centered on the energized region is also depicted. 11 Chapter 2 angle. The second is that the point density does change with polar angle. It is simpler to describe the three-directional reflection in terms of an azimuth angle that varies from − π to π and a normal angle that varies from 0 to π 2 . In this way the PDF for the three- dimensional case is given in Eq. (2.6.2). f (θ R | θ i ) ∝ cos(θ R )sin (θ R ) = sin (2θ R ) (2.6.2) The first point to make about Eq. (2.6.2) is that it has no dependence on the azimuth angle; this is as expected. The second point is that the cos(θ R ) component is present for the same reasons as given in the two-dimensional case. In order to justify the sin (θ R ) component it is necessary to first consider a shell of arbitrary radius around the excited unit of surface. This is shown in Fig. (2.2). The ray density for reflected waves is determined by the radius of a circle parallel to the surface at a height given by θ r . This value is sin (θ R ) . Figure 2.3: Illustration of a sound ray diffusely reflecting from a sub-surface ●, with random impedance and roughness. The angles displayed are those in Eq. (2.6.3). 12 Chapter 2 The pressure caused by diffuse reflection of a planar wave incident on a sub-surface with randomly rough surface and randomly varying impedance is given in Eq. (2.6.3) [6]. The functions β (x, y ) and ε ( x, y ) are the differences between the values of the admittance and height at a point ( x, y ) on the sub-surface, and their average values. These functions are the realization of random variables. PD = A cos(θ i ) cos(θ r ) e ikr4 × πr3 r4 Z 1 Z1 cos(θ i ) + cos(θ r ) + Z 2 Z 2 ∫∫ (ikβ (x, y ) + k 2 (2.6.3) γ 2 ε (x, y ))exp{i(µ x x + µ y y )}∂x∂y µ x = k (sin (θ i ) cos(φi ) − sin (θ r ) cos(φ r )) (2.6.4) µ y = k (sin (θ i )sin (φi ) − sin (θ r )sin (φ r )) (2.6.5) γ 2 = sin 2 (θ i ) + sin 2 (θ r ) − 2 sin (θ i )sin (θ r ) cos(φi − φ r ) (2.6.6) If it is assumed that the β ( x, y ) and ε ( x, y ) are isotropic as well as independent of each other at a given point, the intensity of the reflected wave I D , is given in Eq. (2.6.7). It is also assumed that correlation functions of β ( x, y ) and ε ( x, y ) are given by Eq. (2.6.8) and Eq. (2.6.9), respectively. The σ ∗2 term in both equations is the standard deviation of the β (x, y ) and ε ( x, y ) functions. The w∗ term is the correlation length of the β (x, y ) and ε ( x, y ) functions. This is the distance that can be traveled without the admittance or boundary position substantially changing. 2 ID = cos(θ i ) cos(θ r ) A2 1 × 2 2 2 ρc π r3 r4 Z 1 Z1 cos(θ i ) + cos(θ r ) + Z 2 Z 2 (2.6.7) ik β ( x, y ) exp{i (µ x + µ y )}∂x∂y 2 + k 2γ 2 ε ( x, y ) exp{i (µ x + µ y )}∂x∂y 2 ∫∫ x y x y ∫∫ 13 Chapter 2 1 d 2 Γβ (d ) = σ β exp− 2 2 wβ (2.6.8) 1 d2 Γε (d ) = σ ε2 exp− 2 2 wε (2.6.9) 2 The two integrals in Eq. (2.6.7) are the two-dimensional spatial Fourier transforms of the surface roughness Fε (ω ) and the spatially varying admittance Fβ (ω ) . Additionally, as kγ is the length of the vector µ projected onto the xy plane, the argument of the Fourier transforms is replaced with kγ . The result of these changes is given in Eq. (2.6.10). Under the assumption that the β ( x, y ) and ε ( x, y ) are statistically independent of each other the cross product of the Fourier functions is zero. Using the identity given in Eq. (2.6.11), Eq. (2.6.10) can be rewritten in the form of Eq. (2.6.12). Indeed it is this identity that allows the intensity of the sound field to be calculated explicitly and, consequently, is the motivation for calculating the intensity of the diffuse wave as opposed to directly finding its pressure. 2 ID = ( cos(θ i ) cos(θ r ) A2 1 2 2 2 ρc π r3 r4 Z Z cos(θ i ) + 1 cos(θ r ) + 1 Z 2 Z2 × k 2 Fβ (kγ ) + k 4 γ 4 Fε (kγ ) 2 (2.6.10) ) 2 ∞ 4π 2 2 F (ω ) = ∫ exp{iωτ }Γ(τ )∂τ Area −∞ (2.6.11) 2 ID = cos(θ i ) cos(θ r ) A2 1 2 2 2 Z 1 π r3 r4 Z Z cos(θ i ) + 1 cos(θ r ) + 1 Z 2 Z2 Area × 4π 2 (2.6.12) ∞ 2∞ k ∫ exp{ikγτ }Γβ (τ )∂τ + k 4γ 4 ∫ exp{ikγτ }Γε (τ )∂τ −∞ −∞ 14 Chapter 2 The integrals in Eq. (2.6.12) are explicitly solvable because the limits of integration are infinite. The computation of the first integral is shown in Eq. (2.6.13) through Eq. (2.6.15); the other integral is computed in exactly the same way. Finally the explicit intensity formula is given in Eq. (2.6.16). An interesting note is that this is not identical to the solution presented in Morse and Ingard [6] which appears to be in error; their solution is presented in Eq. (2.6.17). ∞ ∫ exp{ikγτ }Γβ (τ )∂τ = −∞ = −σ β 2 1 τ 2 2 { } exp ik γτ σ exp ∂τ − β 2 ∫ 2 wβ −∞ ∞ kγwβ2 + iτ 1 exp− k 2 γ 2 wβ2 wβ erf i 2 2 wβ 2 π (2.6.13) τ =∞ τ = −∞ 1 = 2π σ β2 exp− k 2 γ 2 wβ2 wβ 2 (2.6.14) (2.6.15) 2 ID = cos (θ i ) cos (θ r ) A 2 2π Area × ρ c 4π 4 r32 r42 Z 1 Z1 cos (θ i ) + cos (θ r ) + Z Z 2 2 (2.6.16) 2 1 1 k w β σ β2 exp − k 2 γ 2 w β2 + k 2 γ 2 wε σ ε2 exp − k 2 γ 2 wε2 2 2 2 ID = cos(θ i ) cos(θ r ) A 2 Area × 2 2 ρc 2πr3 r4 Z 1 Z1 cos(θ i ) + cos(θ r ) + Z 2 Z 2 (2.6.17) 2 2 2 1 1 k wβ σ β exp− k 2γ 2 wβ2 + k 4γ 4 wε2σ ε2 exp− k 2 γ 2 wε2 2 2 15 Chapter 2 2.7 Inverse theory Inverse theory is the estimation of the parameters of a postulated model of a physical system from observed data. In Eq. (2.7.1), f is a model, θ is the set of unknown parameters, X is the observed independent data and Y is the observed dependent data. Consequently, Eq. (2.7.2) expresses the inverse-theory estimate of the set of parameters, θ . It should be noted that, in general, the observed data will consist of a set of independent data points that are matched with a set of dependent data points by the model. Y = f (θ | X ) (2.7.1) θˆ = f −1 (Y | X ) (2.7.2) The solution given in Eq. (2.7.2) will, in general, either not exist (have no closed form) or be under-specified, thus having an infinite number of solutions. Accordingly, inverse theory is solving these two problems. This section will provide some basic solution techniques that will be used later in this thesis. These techniques are numerical estimation, percentile estimation, and statistical techniques (Bayesian inversion and maximum likelihood estimation). Numerical estimation, in general, involves creating an objective function, and maximizing that objective function in an iterative manner. Two common iterative algorithms are the Newton Raphson algorithm (NR) and the Conjugate Gradient method (CG)[7]. The basic equation of the NR is given in Eq. (2.7.3). x n +1 = x n − f (xn ) f ′( x n ) (2.7.3) In Eq. (2.7.3), x n will approach the value needed to make f ( x n ) equal zero as n increases. This is an example of fixed-point iteration. In the context of optimization it is 16 Chapter 2 less common that the objective is to find the root of a function; more often the goal is to find the maximum or minimum of a function. If the function in question is assumed to be smooth then Eq. (2.7.3) can be rewritten as Eq. (2.7.4). This will now locate the extrema of the function f , as such values will only occur if the derivative is equal to zero. x n +1 = x n − f ′( x n ) f ′′( x n ) (2.7.4) The formulas for the CG method are given in Eq. (2.7.5) through Eq. (2.7.7). The scalar multiple β shown is the Flecher-Reeves β ; however other choices are available. The α n is the step size; this term is given the subscript n , as in many applications it is useful to have a varying step size. xn +1 = xn + α n Λxn (2.7.5) Λxn = β n Λxn −1 − ∇f ( xn ) (2.7.6) ∆xnT ∆xn βn = T ∆xn −1∆xn −1 (2.7.7) A significant problem of numerical estimation on its own is that it is difficult to ascertain to what degree of certainty an estimate is valid. Percentile estimation is one solution to this problem; percentile estimation is simply a process to estimate an arbitrary percentile of the distribution of a parameter. Before continuing, it is useful to go through the short proof of this technique. Ax if ( x >= 0) g ( x) = − Bx if ( x < 0) (2.7.8) ∫ g (θ − θˆ )f (θ )∂θ (2.7.9) Error = ∞ −∞ 17 Chapter 2 Eq. (2.7.8) can be conceptualized as an absolute-value function, as it will always return a positive number. However, unlike a conventional absolute-value function, the slopes of the two sides are different. In Eq. (2.7.9), f is the distribution of the unknown parameter θ ; this distribution can also be thought of as the certainty of θ for parameters for which a distribution seems conceptually unsettling. Consequently, Eq. (2.7.9) can be ( ) ( ) thought of as the expectation of g θ − θˆ . ∞ Error = A * ∫ θ − θˆ f (θ )∂θ + B * θˆ ∫ (θˆ − θ )f (θ )∂θ θˆ ∞ ∞ θˆ θˆ θˆ θˆ −∞ −∞ Error = A * ∫ θ f (θ )∂θ − Aθˆ * ∫ f (θ )∂θ + Bθˆ * (2.7.10) −∞ ∫ f (θ )∂θ − B * ∫ θf (θ )∂θ (2.7.11) If the integral in Eq. (2.7.9) is decomposed into two pieces such that the difference of θ and θˆ is always positive, then it yields Eq. (2.7.10). This equation can be further decomposed into four terms, as given in Eq. (2.7.11). The goal of this process so far has been to find the θˆ that minimizes the error value; such a θˆ would be a critical value of () () function Error θˆ . Consequently, if the derivative of Error θˆ is taken with respect to θˆ , and set equal to zero, and this equation is then solved for θˆ , the solution will be the () value of θˆ that minimizes the Error θˆ . ∞ ∞ θˆ θˆ ∂ ˆ ˆ 0= A * ∫ θ f (θ )∂θ − Aθ * ∫ f (θ )∂θ + Bθ * ∫ f (θ )∂θ − B * ∫ θf (θ )∂θ ∂θˆ −∞ −∞ θˆ θˆ [] [] 0 = − A + A * F θˆ + B * F θˆ [] A ⇒ = F θˆ A+ B (2.7.12) (2.7.13) (2.7.14) Eq. (2.7.12) is the derivative of Eq. (2.7.11) with respect to θˆ . F is the cumulative density function (CDF) associated with the probability function, f . The solution is given 18 Chapter 2 in Eq. (2.7.14). This solution can be interpreted as follows: since F is the CDF of θ , θˆ must equal the percentile value associated with the percentage given on the left-hand side of Eq. (2.7.14); that is to say, if A is equal to B , θˆ would be the 50th percentile of θ , the median. Alternatively, if A were equal to five percent of the sum of A and B , θˆ would be the 5th percentile of θ . In the case where a model is over-specified, it is normally helpful to assume that each measurement has a value that is distributed by some assumed distribution function f around the true value. The true value can then be estimated using standard statistical techniques, such as Bayesian inversion, and maximum likelihood estimation. The first of these to be discussed will be maximum likelihood estimation; in order to do this, the likelihood function must first be introduced. A likelihood function is the product of the probability of each element of observed data, given the parameter θ . Consequently, it can be thought of as the probability of obtaining the sampled dataset. For numerical reasons, it is often easier to use the log-likelihood, which is the sum of the logged probability of each element of the dataset. Eq. (2.7.15) and Eq. (2.7.16) give a generic likelihood function and log-likelihood function, respectively. In order to achieve a maximum likelihood estimate, one takes the derivatives of the log-likelihood function, with respect to unknown parameters, and sets them equal to zero. The resultant system of equations can be solved either directly or through numerical methods to find the estimates of the parameters. L = ∏ p( xi | θ ) (2.7.15) log(L ) = ∑ log( p ( xi | θ ) ) (2.7.16) i i The Bayesian approach to statistical estimation is a direct result of Bayes’ rule [8], which is given in Eq. (2.7.17). This states that, for any two events A and B , the probability of 19 Chapter 2 both A and B is equal to the product of the probability of B and the probability of A given B . p( A & B ) ∝ p ( A | B ) ∗ p(B ) (2.7.17) For Bayesian estimation, event A is the sampled data and event B is a particular value of the unknown parameters. In this context Eq. (2.7.17) can be rewritten as Eq. (2.7.18), where π (θ ) is the prior distribution of θ , and L is the likelihood function of data given θ. p( X ,θ ) = L( X | θ ) ∗ π (θ ) (2.7.18) The prior distribution is used to represent prior knowledge of a given parameter. For example, the acoustical absorption α of a surface is required to have a value between 0 and 1. Thus, a reasonable prior distribution would be uniform between 0 and 1, and 0 elsewhere. If specific knowledge of the absorption is known, as would be the case for a sound absorber, then the prior distribution of α could be non-uniform over the range between 0 and 1, assigning a higher probability to values which are closer to the values of the manufacturer’s claim of α . The prior distribution is also often chosen for its algebraic properties and, ideally, is the conjugate distribution for the likelihood function of the data. For example, if the data has a binomial distribution then the prior should have a beta distribution, for most efficient estimation. This is a common criticism of the Bayesian technique, as computational ease is unrelated to prior knowledge. This concern is, in general, ill-founded, as most distributions have a non-informative case that can be used for prior distribution without explicit information, and are sufficiently flexible to represent all possible cases. f (θ ) ∝ L( X | θ ) ∗ π (θ ) (2.7.19) 20 Chapter 2 Because the data have already been observed, it is reasonable to consider it as an arbitrary constant on the left-hand side of Eq. (2.7.18). It is also useful, for computational convenience, to consider Eq. (2.7.18) as being a proportionality relation, as opposed to an equality (the reasons for this will be discussed later in this section.). This results in Eq. (2.7.19), where the left-hand side is the posterior distribution of θ . If the posterior distribution has an explicit form, it can be used for point estimates—normally the median, mean or mode of θ —and for conference intervals on values of θ . In the cases where the posterior distribution does not have an explicit form for these values, there are several techniques to sample from an arbitrary function that can be used to generate an empirical estimation of the posterior distribution. The techniques discussed in this document are based on Markov-Chain Monte-Carlo methods (MCMC). MCMC methods work by generating a series of random values that will approximate an independent, identically-distributed sample for the desired posterior distribution; the approximation will get better as the length of the sequence is increased. Two techniques for implementing MCMC are Gibbs sampling, and the Hastings Metropolis algorithm (HMA). Gibbs sampling is also called alternative conditional sampling. Consider a model with K parameters ( θ1 θ 2 ... θ K ). Each iteration of the Gibbs sampler algorithm loops from 1 to K , sampling the new value of the kth parameter given all of the others. Eq. (2.7.20) gives the distribution of the kth parameter for the tth iteration of the Gibbs sampler algorithm. ( θ t ,k ~ f θ t ,k | θ −t −k1 ) θ −t −k1 = (θ1t ,θ 2t ,...., θ kt −1 ,θ kt −+11...θ Kt −1 ) (2.7.20) (2.7.21) It is helpful to note that that Gibbs sampling may be modified to allow for a kth set of parameters to be updated at each of the iterations. The Gibbs sampler will sample from a distribution that will, in the long term, approximate the desired distribution. The second technique is the Hastings-Metropolis algorithm. The HMA works in a similar manner to Gibbs sampling, but does not require a direct way to sample from the desired 21 Chapter 2 distribution, even one parameter at a time. In order to begin the HMA, a set of starting values for the parameters is required. The only restriction on the starting values is that the posterior density function must be greater then zero for the staring values. The algorithm then iterates T times; each iteration has K sub-steps, one for each of the K parameters of the model. In order to describe the algorithm it is useful to consider each step as composed of three stages. The first stage for the kth step of the tth iteration of the HMA is to select a new candidate value. This can be done in a variety of ways; however the algorithm is fastest if it can be sampled from the desired distribution. As this distribution is unknown, it is normally efficient to select it from a distribution which is either a “likely suspect” or to add random a symmetric random variable to the current value. θ tC,k ~ J x (x;θ t −1,k ) (2.7.22) The second stage is to calculate the probability of accepting the proposed value. The probability is denoted r and is given in Eq. (2.7.23). It can be interpreted as the ratio of the probability of the candidate value to the current value of the kth parameter. The second term is a correction factor to allow the jump distribution to not be symmetric about the current value of the parameter. ( ( ) ( ) ( ) ) p θ tC,k | Y ,θ t ,1 ,θ t , 2 ,..., θ t −1,k −1 ,θ t −1,k +1 ,... J θ t −1,k | θ tC,k r = Min * ,1 C pθ | Y , θ , θ ,..., θ , θ ,... J θ | θ t −1, k t ,1 t,2 t −1, k −1 t −1, k +1 t ,k t −1, k (2.7.23) The final stage is to generate a uniform random value (denoted U ) on the interval [0,1] . The tth value is then selected given the formula in Eq. (2.7.24). θC U ≤ r θ t , k = t ,k θ t −1,k U > r (2.7.24) 22 Chapter 2 A useful attribute of both the Gibbs sampler and the MHA is that they can be used collectively to estimate a set of parameters. Thus, if it is relatively easy to sample from a distribution for only one of the K parameters of a model, and it is computationally expensive to calculate the likelihood function of the data given the parameters, then it is possible to implement Gibbs sampling for the parameter in question and use MHA for the other parameters. 23 Chapter 3 Chapter 3 Literature review As stated in the introduction the objective of the research reported in this thesis was to utilize inverse-method techniques from other disciplines in the investigation of architectural-acoustics problems. The specific architectural-acoustic problem of interest was determining the diffusion of sound as it reflects from a surface of unknown impedance. Thus this section will review published literature in two broad categories: the application of inverse methods, and the measurement of the acoustical properties of surfaces. These categories are not mutually exclusive, as commonly inverse methods are used to estimate the acoustical properties of surfaces. The first of these categories to be reviewed is the application of inverse methods. 3.1 Inverse methods Inverse methods are commonly used in many areas of acoustics; consequently a wide variety of inverse methods are involved. This section first discusses papers that use numerical methods and then papers that use Bayesian inversion. An application of numerical methods is the paper by Taherzadeh and Attenborough [9]. In it the model for the sound field above a semi-infinite impedance surface, Eq. (2.5.8), was inverted to find estimates of the impedance for various grounds. The inversion method used was to find the root of Eq. (3.1.1) with respect to Zˆ , where Qob is the observed reflection factor given source location, frequency and power. The error is optimized using the Newton Raphson algorithm. An interesting aspect of the research is that the impedance estimates are later used as observed data in order to estimate soil characteristics. ( ) Error = Qob − Q Zˆ 2 (3.1.1) 24 Chapter 3 Kanzler and Oelze [10] used the congregate-gradient algorithm for improved scatterersize estimation using backscatter-coefficient measurements with coded excitation and pulse compression. This paper describes a method of estimating the size of acousticallyhard spherical scatterers located in an otherwise homogenous medium. It also should be noted that the inversions done in both Taherzadeh and Attenborough [9] and Kanzler and Oelze [10] papers are exact inversion; this is because no approximations to the forward models were used. An example of an approximate inversion is in the paper by Poole, Frisk, Lynch and Pierce [11]; this paper describes a method for creating estimates of sea-floor acoustic parameters. The estimates are found by solving the system of equations given in Eq. (3.1.3) for θ using the Moore Penrose pseudo-inverse. In Eq. (3.1.3), J is an n by p matrix, where n is the number of data points and p is the number of parameters to be [ ] estimated. The vector θ is composed of the p unknown parameters; Y and X = X1 X2 .. Xp are the measured dependent and independent data, respectively, as defined in section (2.7). ( ∂f X 1 ;θ J = ∂θ1 ) ( ∂f X 2 ;θ ∂θ 2 ) .... ( ) ∂f X p ;θ ∂θ p Y = Jθ (3.1.2) (3.1.3) With respect to inverse methods (that is, ignoring that different forward models are being used to relate different acoustical parameters to different data sets), the difference between Poole, Frisk, Lynch and Pierce [11] and Taherzadeh and Attenborough [9] is that Taherzadeh and Attenborough [9] has a numerical approximation to the minimum error values of the parameters of the forward model, where as Poole, Frisk, Lynch and Pierce [11] has an exact solution for the parameters of a linear approximation of the forward model. 25 Chapter 3 It may seem that having an approximate solution to a problem is better than an exact solution to an approximate problem. There are several advantages to the technique used in Poole, Frisk, Lynch and Pierce [11]. The first of these is computational time; if the forward model is slow it is simply faster to have an explicit solution. Another advantage is that if a Gaussian distribution is assumed for the residuals, the process will also return standard deviations for the parameter estimates. The third advantage of the process used is that the approximation to the model can be polynomial, while still allowing for an exact solution. In order to implement this it must also be possible to take additional partial derivatives of the model. The exact and approximate inversions are not mutually exclusive. Goutsias and Mendel [12] describes the inversion of a model consisting of the linear combination of four nonlinear models. A more extreme example of a hybrid between exact and approximate inversion is Too, Chen and Hwang [13]. In this paper the not-analytically-invertible forward model is approximated by an artificial neural network which is invertible. Equations (3.1.4) through (3.1.6) give the formula for the quantity y k ,l which is the value encoded or stored at the kth neurode in the lth level of the network. The parameters of the artificial neural network are the wm ,k ,l ’s, and the x m ,l ’s are the input values (note that if 1>2 then the input values for a layer will be the output values of the pervious layer and not necessarily the model input values.). In Too, Chen and Hwang [13] the so-called transfer function ϕ (v k ,l ) is the hyperbolic tangent. If the transfer function was the identity function and the network had only two levels (an input level and output level), then the model would be an average of m linear models. As the transfer function is not the identity function and the network has more than two layers it is a model consisting of the linear combination of many non-linear models that are themselves the smaller artificial neural networks (which are linear combinations of smaller neural networks). m u k ,l = ∑ wm ,k ,l x m ,(l −1) (3.1.4) v k ,l = u k ,l + bk ,l (3.1.5) j =1 26 Chapter 3 y k ,l = ϕ (v k ,l ) = tanh (v k ,l ) (3.1.6) Bayesian inversion is a popular tool for creating estimates of unknown parameters, although its use seems to be currently limited to underwater acoustics. A strong example of an underwater-acoustics paper that uses Bayesian inversion is Dosso and Wilmut [14]. This paper describes explicit methods for interpreting the posterior distribution. The parameters estimates are defined as the maximum a posteriori distribution (MAP), as given in Eq. (3.1.7), where f (θ ) is the posterior probability distribution of the parameters θ . The one-dimensional and two-dimensional parameter distributions are given in Eq. (3.1.8) and Eq. (3.1.9), respectively. θˆ = Arg max ( f (θ )) ( ) ( (3.1.7) ) ′ ′ f i θˆi = ∫ δ θ i′ − θˆi f θ ∂θ = ∫ ∫ ...∫ ∫ ...∫ f (θ1′,θ 2′ ...θ i′−1 , θ i′+1 ...θ P′ )∂θ P′ ...∂θ i′+1∂θ i′−1 ...∂θ 2′ ∂θ ′ ( ) ( )( θ i =θˆi ) ′ ′ f i θˆi , θˆ j = ∫ δ θ i′ − θˆi δ θ ′j − θˆ j f θ ∂θ (3.1.8) (3.1.9) Dosso and Wilmut [14] also describes some fast sampling methods. These are based on the observation that the best distribution for sampling a proposed value of an unknown parameter is the posterior distribution. This choice is not possible since if the posterior distribution was known explicitly, then it would not be necessary to estimate it. However if a reasonable approximation of the posterior distribution is used, the sampling speed can be increased. There is also no possibility that the approximate distribution will create error or bias in the estimation process as the MHA will correct any errors in the proposed distribution. The proposed distribution used in Dosso and Wilmut [14] is the posterior distribution that would have resulted if the inverse problem was linearly approximated in similar way to Poole, Frisk, Lynch and Pierce [11]. An additional benefit of this method is that a rotation of the parameters ϑ can be sampled from instead of the parameters themselves. As a consequence, the algorithm is much faster when trying to estimate 27 Chapter 3 highly correlated parameters. The proposed distribution for the rotated parameters is given in Eq. (3.1.10). The rotation matrix U is defined as the column eigenvector matrix of Σ ; this is why Eq. (3.1.12) is valid, where W is the diagonal matrix of eigenvalues corresponding to the eigenvector columns of U . Finally the matrix J is composed of partial derivatives of the forward model in a similar manner to Eq. (3.1.2). ( )= f proposal ϑ 1 C (2π ) P 2 Σ 12 {( ) ( T exp − ϑ − ϑ Σ −1 ϑ − ϑ C C )} ϑ =UTθ ( −1 1 Σ = UWU T = J T Σ data J + Σ −prior (3.1.10) (3.1.11) ) −1 (3.1.12) Of particular interest is the paper de Vries, Joeman and Schreurs [15] as it used inverse methods to estimate the scattering coefficient of an architectural surface. The inversion process used is the boundary imaging method. An impulse response is measured at each point on an array parallel to and offset from the surface of interest. Reflections from the boundary are extrapolated to the position of their origin using their Rayleigh representation integral. 3.2 Direct surface-acoustic measures This section reviews conventional direct methods of measuring surface impedance and diffusion. The first measurement procedures to be considered are those for impedance (or absorption). A well-known way to measure the impedance of a surface is the impedance-tube method, described in ISO 10534-2 [16]. The method utilizes two microphones, a wave guide (or impedance tube), a white-noise source, and a digital frequency-analysis system. The microphones, loudspeaker and sample of interest are arranged as shown in Fig. (3.1). Equation (3.2.1) gives the transfer function H 1, 2 that is used to calculate the impedance values. The quantities S1,1 and S 2, 2 are the auto-spectra of the signals received at the first 28 Chapter 3 Figure 3.1: Diagram of a typical impedance tube with source on the right-hand side and sample on the lefthand side. and second microphones, respectively. The quantities S 2,1 and S1, 2 are the cross-spectra of the signal at the second microphone with respect the first, and of the first microphone with respect to the second, respectively. Equation (3.2.2) gives the reflection coefficient Rp in terms of the transfer function. Finally Eq. (3.2.3) gives the impedance of the surface in terms of reflection coefficient. This is simply solving Eq. (2.5.4) for Rp (θ I ,θ T , Z ) , where both the incident and transmitted angles are zero. H 1, 2 S1, 2 S 2, 2 = S1,1 S 2,1 Rp = 12 H 1, 2 − H Im H Re − H 1, 2 Z= = H Re + iH im (3.2.1) exp{2ikx} (3.2.2) 1 + Rp 1 − Rp (3.2.3) Another common method of measuring the impedance of a surface is the sphericaldecoupling method (SDM). This method is described in De Geetere [17]. The SDM is similar to the impedance-tube method; however it is a free-field method. The transfer function H 1, 2 is defined as for the impedance-tube method and is given in Eq. (3.2.1). The reflection coefficient is given by Eq. (3.2.4). The distances and angles are displayed in Fig. (3.2). Equation (3.2.5) gives the impedance of the surface. It also should be noted that the sample and microphones should be located in the far field of the source. 29 Chapter 3 Figure 3.2: Illustrates a typical sphericaldecoupling measurement setup. The diagram has been augmented to include the image receivers , which exist only conceptually. exp{ikr1 } exp{ikr2 } − H 1, 2 r1 r2′ Rp = exp{ikr1′} exp{ikr2′} − H 1, 2 r1′ r2′ Z= 1 + Rp 1 1 − Rp cos(θ )(1 − ikr0 ) (3.2.4) (3.2.5) A comparison was made by Graves and Hodgson [18] between impedance measurements made using the impedance-tube method and the SDM. While the results of the two methods were found to be different, neither was considered preferred. A possible reason for the discrepancy in the results is that the SDM assumes wholly specular reflection and only approximates the contribution of the spherically-reflecting wave. There seems less reason why the impedance-tube method would be inaccurate; however it is in general impractical to implement, as it requires a disk-shaped sample of the material of interest. Also it only can measure the impedance at normal incidence. If only the absorption averaged over all incident angles is of interest then a third method using diffuse-field theory is possible. This method is given in the ASTM C423-07a standard [19]. The method works by comparing the reverberation times in a parallelepiped reverberation chamber with and without the sample of interest. The total area Α of wholly absorbing material in a parallelepiped is given in Eq. (3.2.6). V is the 30 Chapter 3 volume of the parallelepiped and d is rate of sound decay (dB per second). Let the total area of absorption in the reverberation chamber without the sample of interest be denoted Α1 and the total area of absorption in the reverberation chamber with the sample of interest be denoted Α2 . If the surface area of the sample of interest is S , then the absorption coefficient of the sample is given in Eq. (3.2.7). Α = 0.9210 α= Vd c ( Α2 − Α1 ) + S S Α1 S RC (3.2.6) (3.2.7) Direct methods also exist for measuring the diffusion or scattering characteristics of a surface. The first method to be considered is taken from ISO 17497 1 [20]; this will be referred to as the reverberation-chamber method (RCM). The method described is for measuring the random-incidence scattering coefficient of surfaces caused by surface roughness. The measurement results can be used to describe what proportion of the sound reflecting from a surface deviates from the specular angle. Similar to the previous method for measuring sound absorption, measurements take place in a reverberation chamber. Indeed it can be viewed as an additional step to the previous method. The additional step is that the impulse response of the reverberation chamber is taken a third set of times with the sample of interest rotating. This third measurement gives the apparent randomincidence absorption, denoted α S . The apparent random incidence absorption is larger than the conventional absorption; the increase is due to the apparent absorption caused by scattering and thus measures the scattering coefficient. The formula for α S is given in Eq. (3.2.8). As before, V is the reverberation-chamber volume and S is the surface area. The quantities TS and Te are the reverberation times for the chamber with and without the sample; in both cases the turntable used to rotate the sample should be present and rotating. It is important to note that TS is not calculated by the standard method of Schroeder’s integrated impulse technique [21]; instead 16 impulse responses taken while 31 Chapter 3 the sample is rotating are averaged, then the reverberation time is calculated in the conventional fashion. The scattering coefficient is then calculated using Eq. (3.2.9). V 1 1 + S c s Ts ceTe α S = 55.3 s= 4V − (mS − me ) S α −αs 1−αs (3.2.8) (3.2.9) A similar method to the RCM is described Vorlander and Mommertz [22]. This so-called free-field method (FFM) is similar to the RCM, as they both utilize the assumption that scattered reflections are incoherent for different orientations of the surface. Eq. (3.2.9) is still used to calculate the scattering coefficient s ; however the formulae for α and α s are different and are given in Eq. (3.2.10) and Eq. (3.2.11), respectively. The value Rp,i is the complex reflection coefficient measured at the ith orientation of the surface. α =1− 1 n ∑ Rp , i n i =1 αs = 1− 1 n ∑ Rp,i n i =1 2 (3.2.10) 2 (3.2.11) 32 Chapter 4 Chapter 4 Experimental setup, loudspeaker characteristics and data This chapter consists of three sections. The first describes the experimental setup for making steady-state sound-pressure-level L p measurements. Section 2 gives a detailed description of the characteristics of the sound source used in these tests, as well as the method by which the characteristics were calculated or estimated. Finally, in the third section, the measured L p values are presented. 4.1 Experimental setup This section describes the method by which steady-state sound-level measurements were made. These measurements were used in an inversion process for models described in Chapter 5, in order to create estimates of the model parameters—in particular, surface impedance and diffusion coefficient. A brief description of the experimental setup and procedure is as follows. A 12 ft × 12 ft rigid surface (¾-in-thick painted plywood on studs with 24-in spacing) was constructed in the anechoic chamber. Materials of interest were placed on top of the rigid surface to create other test surfaces. These surfaces were selected to represent typical architectural surfaces. The sound field was irradiated by a single loudspeaker that was placed approximately one meter above the centre of the test surface. A microphone was then suspended above the surface. Thirty-two L p measurements were then taken for each surface at six pure-tone frequencies and once with white noise; each of the 32 measurements was conducted at a different location within the chamber. 33 Chapter 4 Figure 4.1: A plan-view representation of the experimental setup. A plan diagram of the experimental setup is shown in Fig. 4.1. The side length of the sample is 3.65 m. The dimensions of the anechoic chamber were 4.7 m by 4.3 m by 2.3 m. The microphone was placed every 15.02 cm along Track 1 at a height of 16 cm. The microphone was placed every 7.56 cm on Track 2 at a height of 58 cm. If the lower lefthand corner of the sample surface is taken as the reference, with horizontal direction as the x axis and the vertical as the y axis, then the ends of Track 1 were at x = 52 cm, y = 251 cm and x = 228 cm, y = 250 cm. The ends of Track 2 were at x = 73 cm, y = 64 cm and x = 74 cm, y = 179 cm. The source was located at x = 188 cm, y = 170 cm and z = 100 cm. Table 4.1: The loudspeaker amplitude for the six pure-tone frequencies. 34 Chapter 4 4.2 Loudspeaker characteristics Table 4.2: The loudspeaker amplitude for the white noise in octave bands. The amplitude in decibels of the source at the pure-tone test frequencies is given in Table 4.1. The amplitude in decibels for the white noise in octave bands is displayed in Table 4.2. These values were calculated by solving for the A term in the special case with Q equal to zero in Eq. (2.5.8), as shown for one measurement in Eq. (4.2.1). The L p ,Empty term is the observed L p for the empty anechoic chamber. As one A value would have to fit all 32 measurements, the system of equations generated is over-specified; consequently Aˆ , the value of A that minimizes the absolute error, given in Eq. (4.2.2), was used. The same process was used to find the amplitudes of the source for white noise. The confidence intervals are calculated using the percentile-estimation technique described in Chapter 2, Section 7. exp{ikr1 } + A = A − 20 log10 (r1 ) = L p ,Empty 20 log10 r1 Aˆ = ∑ i =32 L p ,Empty ,i + 20 log10 (r1,i ) 32 (4.2.1) (4.2.2) Figures 4.2 and 4.3 display the measured L p ,Empty and the predicted Lˆ p ,Empty ; the prediction is made assuming that there is no loudspeaker directivity. As can be seen, there is poor agreement at the higher frequencies. Although there could be many sources of error, it is assumed to be a result of loudspeaker directivity. 35 Chapter 4 Figure 4.2: The measured ○, and predicted , pure-tone sound-pressure levels in the empty anechoic chamber assuming no loudspeaker directivity, plotted against the direct distance from the source to the measurement location. 36 Chapter 4 Figure 4.3: The measured ○, and predicted , white-noise sound-pressure levels in the empty anechoic chamber assuming no loudspeaker directivity, plotted against the direct distance from the source to the measurement location. The directivity is over both the azimuth and polar angles. For the purpose of the work presented later in this thesis, the directivity variation with azimuth angle will be assumed to be independent of that with polar angle. Additionally they are assumed to be dB additive, and are calculated from Eq. (4.2.3). The functions δ θ (θ ) and δ φ (φ ) are the directivity functions. Figure 4.4 shows the polar-angle directivity of the speaker at the six pure-tone frequencies of interest with a loess smooth. 37 Chapter 4 Figure 4.4: The measured ○, and smoothed , loudspeaker directivities (Delta) plotted against the polar angle (theta) in radians. Aˆ − 20 log(r1 ) + δ θ (θ ) + δ φ (φ ) = LP , Empty (4.2.3) Although it would be possible to use a look-up table to account for the functional value of δ θ (θ ) , this would inhibit any further analytical work from being conducted. Fortunately, knowledge of the directivity is only required in the range (− π 2 , π 2] , and the directivity function seems to have a near quadratic relation with θ in this range. Accordingly, the polar directivity was approximated by a quadratic polynomial for each frequency; an 38 Chapter 4 Table 4.3: Parameter values of the directivity function at each frequency. arbitrary quadratic polynomial is shown in Eq. (4.2.4). The parameters of the polynomials are listed in Table 4.3. δθ (θ ) ≈ Aθ 2 + Bθ + C (4.2.4) Once the polar directivity was accounted for, the azimuth directivity of the loudspeaker could also be considered. Figure 4.5 shows the relationship between the azimuth angle and the loudspeaker directivity. The relationship between the two values is less clear; in general there appears to only be a significant relationship at the higher frequencies (4000 and 8000 Hz); consequently the function δ φ (φ ) was assumed to be uniformly zero for all frequencies below 4000 Hz. For the higher frequencies, the loess curve shown in Fig. 4.5 was used. 39 Chapter 4 Figure 4.5: The measured ○, and smoothed , loudspeaker directivities (Delta) plotted against the azimuth angle (Phi) in radians. 4.3 Data description As noted above, there were seven test surfaces of interest. These were big blocks (BB), small blocks (SB), FORESTwall (FW) [manufactured by Morinwood, Inc., 1511 Fell Street, Victoria, B.C., V8R 4V9; http://www.forestwall.com/], acoustical baffles (Soft), Corelam plywood [available from Corelam Ltd., 715 W 69th Ave, Vancouver BC V6P 2W2; http://www.corelam.com/] aligned in the x direction (Sine1), Corelam plywood aligned in the y direction (Sine2), and the reference surface (Flat). This section describes each of these surfaces, and the sound fields measured above them. 40 Chapter 4 The Flat surface consisted of ¾-in plywood above studs with 24-in spacing. The plywood was painted, to reduce porosity and therefore absorption; however it was composed of five separate panels. The BB surface consisted of 30-cm, cubic, varnished wood blocks randomly scattered over the top of the reference surface. The SB surface consisted of 10cm lengths of unpainted “two by fours” nailed together. These blocks were then randomly distributed over the reference surface. A picture of this configuration is shown in Fig. 4.6. The next test surface was the Soft surface. This consisted of 18 acoustical baffles, approximately 3-cm thick, covering the reference surface. The surfaces Sine1 and Sine2 consisted of 2.8-mm thick plywood pressed into a sinusoidal shape along one axis with a 5-mm amplitude. For each of the surfaces 224 measurements were made. These measurements were made at 32 locations with each of six pure-tone and one white-noise excitations of the sound field. The pure-tone sounds were made at 250, 500, 1000, 2000, 4000 and 8000 Hz. Also, as noted above, the 32 measurement locations were divided along two measurement Figure 4.6: The experimental setup while the SB surface was being evaluated (the longer obstacles were not present at the time of measurement). 41 Chapter 4 tracks. Figures 4.7 and 4.8 display the measured L p values for Tracks 1 and 2, respectively, at the six pure-tone frequencies. Figures 4.9 and 4.10 display the measured L p values for Tracks 1 and 2, respectively, when the sound field was energized by white noise. The frequencies in Figs. 4.9 and 4.10 refer to the ⅓-octave band centered at that frequency. Figure 4.7: The measured steady-state sound-pressure levels (pure-tone excitation) above the surfaces ○f interest, plotted against the x (Track 1) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 ). 42 Chapter 4 Figure 4.8: The measured steady-state sound-pressure levels (pure-tone excitation) above the surfaces ○f interest, plotted against the y (Track 2) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 .) 43 Chapter 4 Figure 4.9: The measured steady-state sound-pressure levels (broad band excitation) above the surfaces ○f interest, plotted against the x (Track 1) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 ). 44 Chapter 4 Figure 4.10: The measured steady-state sound-pressure levels (broad band excitation) above the surfaces ○f interest, plotted against the y (Track 2) coordinate of the measurement location (BB , SB , FW , Flat , Soft , Sine1 , and Sine2 ). 45 Chapter 5 Chapter 5 Forward models and inversion procedure This chapter describes the three heuristic forward models for predicting the RMS soundpressure level above a diffusely-reflecting impedance surface. The models are the imagesurface model (ISM), the image-surface semi-phase model (ISSPM) and the diffusereflection energy model (DREM). Also discussed are the inverse methods that were used to find estimates of parameters for each of the models. Before continuing with an explicit description of the models it seems useful to give a general one. All of the models that were considered predict the RMS sound-pressure level L p at a point above an impedance surface in a half-space irradiated by a single source. The models take as their arguments the three-dimensional coordinates of the receiver, L = [Lx Ly Lz ] . The parameters of the models are the source location T S = [Sx Sy Sz ] , the source amplitude A , the surface impedance Z = Z re + iZ im , and T the surface diffusion coefficient D . The surface diffusion coefficient is a parameter defined for use in these models; it will be described in detail later, however it is similar to the scattering coefficient described in Vorlander and Mommertz [22]. The model will be inverted to solve for the parameters Z and D . It was possible to directly measure the amplitude and source-location parameters. Therefore they are treated as constants, known from their measured values. The estimated parameters were restricted to a physically meaningful range. For the surface impedance this meant that the real component must be positive, while the imaginary component was unrestricted [1]. As the diffusion coefficient was created for the convenience of the models presented in this thesis, its physically meaningful range is less obvious. Thus the most conservative set of bounds was used; these are zero and one. 46 Chapter 5 5.1 Image-surface and image-surface semi-phase models Both ISM and ISSPM assume that each of the observed L p values is a random variable. The distribution of L p ,i , the ith value of L p , is the Gaussian distribution and is given in Eq. (5.1.1). There are two important points to observe about Eq. (5.1.1). First is that the mean of L p ,i , η i , does not in general equal the mean of L p , j , η j . Second, the variance σ 2 of the distribution is constant over all observations. L p ,i ~ (L p ,i − η i )2 exp− 2σ 2 2πσ 1 (5.1.1) The term η i is the systematic component of the model. Equation (5.1.2) is the formula for η ; note that the subscript has been dropped to avoid the notation becoming too cumbersome. η can also be called L p ,predicted . η = 20 log10 ( Ptotal ) + A (5.1.2) The quantity Ptotal is a heuristically modified version of Eq. (2.5.8). The equation has been augmented such that instead of the two waves (direct and reflected) of Eq. (2.5.8), there are now three waves. The new wave is the diffusely-reflected wave. Equation (5.1.3) gives Ptotal for the ISM; this is not the same as for the ISSPM which will be described later. ′ + PDif ′ Ptotal = PD′ + PSpec (5.1.3) As a result of the significant loudspeaker directivity discussed in Chapter 4 each wave must be multiplied by a directivity term ξθ ,φ (θ , φ ) . Equation (5.1.4) gives the formula for 47 Chapter 5 the directivity in terms of δ θ (θ ) and δ φ (φ ) , the directivity terms from Eq. (4.2.3). The values A f , B f and C f are frequency dependent and are listed in Table 4.3. ξ θ ,φ (θ , φ ) = 10 δ θ (θ )+δ φ (φ ) 20 ≈ 10 A f θ 2 + B f θ + C f +δ φ (φ ) 20 ′ Ptotal = ξ θ ,φ (θ D , φ D )PD + ξ θ ,φ (θ Spec , φ Spec )PSpec + PDif (5.1.4) (5.1.5) The inclusion of directivity in the diffusely-reflecting wave is more complicated, as it is a sum of many waves all radiating from the source at different angles. Its directivity term will be described later. The formulae for the direct and reflected pressure contributions are displayed in Eq. (5.1.6) and Eq. (5.1.7), respectively. The F term is still the bottomloss factor, as described in Eq. (2.5.14); RP is the reflection coefficient defined in Eq. (2.5.4). The parameter D is the diffusion coefficient of the surface described in the introduction to this chapter. PD = exp{ikr1 } r1 PSpec = (1 − D )Q exp{ikr2 } r2 Q = (R p + (1 − R p )F ) (5.1.6) (5.1.7) (5.1.8) The integral that defines the diffusely-reflected wave is given in Eq. (5.1.9). It is important to note that although the integral appears to be free of x and y , the location of the reflection of each diffusely-reflecting ray, all of the angles, as well as the reflected path lengths, are functions of x and y . Equation (5.1.10) gives rD , which is the total distance from the source to the point [x y 0] on the surface and on to the receiver T position. The diffuse-reflection coefficient Γ(θ i ,θ R ) is given in Eq. (5.1.11); the angles are shown in Fig. 5.1, Eq. (5.1.12) and Eq. (5.1.13). The function Φ(θ i ,θ R ) is displayed in Eq. (5.1.14). 48 Chapter 5 uby ubx PDif = D ∫ ∫ ξθ φ (θ , Dif exp{ikrD } ∂x∂y rD (5.1.9) (x − Lx )2 + ( y − Ly )2 + (Lz )2 (5.1.10) , φ Dif )Γ(θ i , θ R )Φ (θ i ,θ R ) lb y lbx (x − Sx )2 + ( y − Sy )2 + (Sz )2 Γ(θ i , θ R ) = + cos(θ D (θ i , θ R )) − cos(θ D (θ i , θ R )) + θ D (θ i , θ R ) = cos(θ T (θ i , θ R )) Z cos(θ T (θ i , θ R )) (5.1.11) Z θi + θ R 2 (5.1.12) θ T (θ i ,θ R ) = θ i − θ D (5.1.13) Φ (θ i ,θ R ) = cos(θ i ) sin (2θ R ) (5.1.14) The ISSPM is similar to the ISM—differing in how Ptotal is defined; Eq. (5.1.15) displays the definition of Ptotal for the ISSPM. The distinction is that the diffuse component is always added in phase to the direct and specular components. Figure 5.1: The path of a sound ray ( ) diffusely reflecting off a conceptual surface (···) which is a rotation of the actual surface ( ). The angles displayed are those in Eq. (5.1.9) through Eq. (5.1.14). 49 Chapter 5 Ptotal = PD + PSpec + PDif (5.1.14) Before continuing, it is important to interpret and heuristically justify the ISM and the ISSPM. The inclusion of the parameter D in these models allows the distribution of reflected energy to vary from wholly specular to wholly Lambert’s Law diffuse. A shortcoming of the model is that one of the possible diffusely-reflected angles is the specular angle; thus D is not exactly the proportion of energy that is reflected in a non-specular direction. Both the ISM and the ISSPM are based on the concept of the image source. For a conventional image-source model, as described in Chapter 2, there is only one image source; this source corresponds to the specularly-reflected wave. However if an image source is assigned to each point on the surface then, instead of a singular image source, there is an image surface. This is similar to the ray-tracing model (version DRAYcube) that is used in Hodgson [23]; given enough rays every point of the surface will reflect a diffusely-reflecting ray to the receiver, consequently every point on the surface can be considered to correspond to an image source. Thus a way of conceptualizing or describing the ISM or the ISSPM is as a virtual diffuse-reflecting ray-tracing model. An example of such an image surface is displayed in Fig. 5.2; the black line is the actual 2D surface; the blue grid is the 3D image surface. Each grid point displayed on the image-surface corresponds to an evenly-spaced grid point of the actual surface. As can be seen, the surface elements of the image-surface are not uniformly spaced like the corresponding elements on the actual surface. The function Φ(θ i ,θ R ) is used to attenuate the contributions of the regions of the image sources that have the smallest surface elements. 50 Chapter 5 Figure 5.2: An illustration of an image surface ( ) generated by a source and receiver pair above a reflecting surface ( ). The image surface is projected into the xz plane, the yz plane and the xy plane. The diffuse-reflection coefficient Γ(θ i ,θ R ) can most readily be interpreted by considering Fig. 5.1. The surface is replaced by a virtual surface that has the local slope required to make the incident and reflected angles equal. However the actual surface is assumed to be locally reactive; as a consequence the horizontal speed of sound in the actual surface must be much smaller than the vertical speed of sound. Thus the transmitted ray must be normal to the actual surface in order to accommodate this; the angle of transmission of the virtual surface will not, in general, be zero. The effect of this is that, as the difference between angles θ i and θ R increases, the surface becomes acoustically harder. This is similar to the reflection observed at the plane interface between two fluids [1]. 5.2 The diffuse-reflection energy model The DREM model is similar to the ISM and the ISSPM in that it models the reflected sound as the proportion D of energy reflecting in a non-specular way. It is distinct in that it does not account for the phase of the sound waves. As a consequence the DREM is only used with the white-noise measurements, for which phase effects would be expected to be small. 51 Chapter 5 Similarly to the ISM and the ISSPM, each of the observed L p values is assumed to be a random variable. However the distribution of L p ,i , the i’th value of L p , is not assumed to have an explicit form. As is the case for Eq. (5.1.1), the mean of L p ,i is η i . The variance is σ 2 , which is constant over all observations. ( LP ,i ~ f L p ,i ;η i , σ 2 ) (5.2.1) Equation (5.2.2) is the formula for η ; again the subscript has been dropped to avoid the notation becoming too cumbersome. 2 2 η = 10 log10 (PD2 + PSpec + PDif )+ A (5.2.2) The definitions of PD , PSpec and PDif are given in Eqs. (5.2.3), (5.2.4) and (5.2.5), respectively. The directivity term ξθ ,φ (θ , φ ) is as defined in Eq. (5.1.4). The energy- reflection coefficient RE is given in terms of the normal-incidence energy absorption coefficient of the surface in Eq. (5.2.6). The diffuse-reflection coefficient is displayed in Eq. (2.5.7); the values θ D and θ T are functions of θ i and θ R , defined in Eq. (5.1.12) and Eq. (5.1.13). PD = ξθ ,φ (θ D , φ D ) PSpec = (1 − D )R E uby ubx PDif = D ∫ ∫ lb y lbx (5.2.3) r1 ξθ ,φ (θ Spec , φ Spec ) (5.2.4) r2 ξθ ,φ (θ Dif , φ Dif )ΓE (θ i ,θ R )Φ(θ i ,θ R ) rD α cos(θ ) − 2 −α RE = α cos(θ ) + 2 −α ∂x∂y (5.2.5) (5.2.6) 52 Chapter 5 α cos(θ D ) − cos(θ T ) 2 −α ΓE (θ i , θ R ) = α cos(θ D ) + cos(θ T ) 2 −α (5.2.7) 5.3 Inversion techniques The goal of the inversion performed in this research was not only to provide estimates of the unknown parameters but also to provide confidence intervals for them. Consequently, Bayesian inversion and percentile estimation were used. The ISM and ISSPM have their parameters estimated through Bayesian inversion; the DREM has its parameters estimated through percentile estimation. The Bayesian inversion will be described first. The probability distribution of a single data point under the ISM and the ISSPM is displayed in Eq. (5.1.1). Consequently the likelihood function of the sampled data is displayed in Eq. (5.3.1). 32 L=∏ i =1 (L p ,i − η i )2 exp− 2σ 2 2π σ 1 (5.3.1) Since both models considered here are new, there is no prior information about the parameters used within the models; thus it is reasonable to assume a non-informative prior distribution for each of the parameters. The non-informative prior distribution used was the Jeffrey’s prior, denoted π J (Z , D, σ 2 )[8]. Equation (5.3.2) displays the formula for a Jeffrey’s prior. The function I (θ ) is the Fisher information of the data with respect to the unknown parameters. The function E ( x ) is the expectation of the random variable x . The particular Jeffery’s prior for the ISM and the ISSPM is given in Eq. (5.3.3). It can be interpreted as a non-normalized uniform distribution over the model ( ) parameters Z and D , as well as log σ 2 . 53 Chapter 5 ∂ log(L(θ )) ∂θ π J (θ ) ∝ det[I (θ )] = det E π J (Z , D, σ 2 ) ∝ (5.3.2) 1 (5.3.3) σ2 The posterior distribution of the models is given in Eq. (5.3.4). For computational efficiently it is useful to consider the log-posterior distribution, which is displayed in Eq. (5.3.5). There is no analytical solution for the integral in Eq. (5.1.9), so it is approximated numerically, as shown in Eq. (5.3.6). The approximation can be interpreted as giving a diffusely-reflected image source to every square centimeter of the surface. ( f Z , D, σ ( log f Z , D, σ 2 ) 2 )= σ ∏ 1 2 32 i =1 (L p ,i − η i )2 = −34 log(σ ) − 16 log(2π ) − ∑ 2σ 2 i =1 PDif ≈ D ∑∑ ξθ ,φ (θ Dif , φ Dif 365 365 l =1 j =1 (L p ,i − η i )2 exp− 2σ 2 2π σ 1 32 (5.3.4) )Γ(θ i ,θ R )Φ(θ i ,θ R ) exp{ikrD } 1 rD 100 (5.3.5) 2 (5.3.6) Now that the approximation to the posterior distribution is well defined, it is possible to describe the inversion procedure. It is important to note that this procedure had to be repeated for each surface test sample at each frequency. Additionally it was repeated for both the ISM and the ISSPM. The inversion consisted of three steps. The first step involved finding starting values for the optimization procedure. The second step involved iteratively taking samples of the parameters. The third step involved summarizing the sampled parameters. The starting values were point estimates of the parameters that maximized the posterior likelihood. The estimates were found through repeat application of the conjugate-gradient (CG) algorithm to the function from many random starting locations. When it was 54 Chapter 5 reasonable to conclude that a global minimum had been found, these values, denoted Z 0 , D0 and σ 02 , were used as the starting values of the estimation process. The iterative samples were obtained by two methods. The Z and D values were obtained through the application of the Metropolis-Hastings algorithm, while the σ 2 values were obtained through Gibbs sampling. Gibbs sampling is the faster of the two methods but has more restrictions on its use, that Z and D fail to meet. The sampling process for the arbitrary kth iteration of the inversion process will now be described. The candidate values of Z and D were selected by adding a random value to current iteration k-1th of Z and D . Equations (5.3.7) and (5.3.8) show the candidate values of Z and D , Z C and DC , in terms of their current values. The random variable J Z , which is used to create a candidate value for Z , is defined in Eq. (5.3.9). The distributions of the real and imaginary parts are given in Eq. (5.3.10) and Eq. (5.3.11), respectively. The meta-parameter τ was tuned such that the model accepted approximately 50 percent of the candidate values. The distribution of the random value J D , which is used to create a candidate value for D , is given in Eq. (5.3.12). The meta-parameter W was also tuned such that approximately 50 percent of the candidate values were accepted. Once the kth candidate values for Z and D were found they were accepted or rejected as described in Chapter 2, Section 7. This method of creating a candidate value was favoured over the use of a proposal distribution to reflect the lack of prior knowledge inherent in the problem. ZC = Zk + J Z (5.3.7) DC = Dk + J D (5.3.8) J z = Re( J z ) + Im( J z )i (5.3.9) Re( J Z ) ~ Re( J Z )2 exp− 2τ 2 2π τ 1 (5.3.10) 55 Chapter 5 Im( J Z ) ~ Im( J Z )2 exp− 2τ 2 2π τ 1 1 | −W ≤ J D ≤ W J D ~ 2W 0 | otherwise (5.3.11) (5.3.12) The second part of the kth iteration is to create a new value for the σ 2 . This can be acquired directly since σ 2 has a scaled inverse − χ 2 distribution. A generic example of the distribution is shown in Eq. (5.3.13). The two parameters of this distribution for the kth step of the estimation processes are displayed in Eq. (5.3.14) and Eq. (5.3.15). These formulas can be found by rewriting the posterior distribution in the form given in Eq. (5.1.16). The sampled value is always accepted, as this is Gibbs sampling not MetropolisHastings. ν ( inverse − χ 2 θ ; v, ς 2 ) vς 2 2 ς2 exp− 2 2 2θ ∝ 1 exp− νς = ν ν 1+ 1+ 2θ ν 2 2 Γ θ θ 2 ν = 32 ∑ (L 32 ς2 = 1 σ2 32 ∏ i =1 i =1 (5.3.13) (5.3.14) − ηi k ) 2 p ,i 32 32 (L p,i − η i )2 ∑ 32 i =1 (L p ,i − η i )2 1 1 32 exp− ∝ 2 (1+32 2 ) exp− 2 2σ 2σ 2 2π σ σ (5.3.15) (5.1.16) The third step of the estimation process was to summarize the sampled parameters. Before any summary statistics could be calculated, it was helpful to burn in, or discard, the first half of the sampled parameters. This process protected the estimates from being 56 Chapter 5 affected by a poorly chosen starting value. The parameter point estimate was, in general, taken as the median of the sampled data with respect to the parameter. The confidence intervals were found in a similar manner, simply using different percentiles. For the cases where the parameter distributions exhibited multi-modality, a subjective decision was made on the most representative value and confidence interval. The optimization that was done for the DREM was significantly simpler than the method that has just been described. The simplification is a result of the DREM model only having two unknown parameters: α and D ; the nuisance parameter σ 2 can be ignored. Both of these parameters are bounded in the set [0 1] . As well, the DREM is significantly faster to implement. Thus it is possible, with a small number of points, to make a brute-force calculation for all combinations of parameters. The array of all parameter values is created by two lists of 1000 random, uniform values [ ] on the range 0 1 . The first list gives the candidate α values; the second list gives the candidate value for D . Equation (5.3.17) shows the objective function, or error function, that was evaluated for the kth parameterization of the model. 32 Errork = ∑ L p ,i − η i (α k , Dk ) (5.3.17) i =1 The process was repeated two more times. On each occasion the sample space was quartered. The new region was a square of side length, one half the side length of the previous iteration. The square was centered on the best-performing parameters found so far. The confidence intervals were found using percentile estimation. This process is only defined for one-dimensional problems; consequently the second parameter was held constant at its optimal value while the percentile was estimated. Equation (5.3.18) displays the objective function that was used to find the lower bound of the confidence 57 Chapter 5 interval; Eq. (5.3.19) gives the objective function for the upper bound. These can be interpreted as finding a 95 percent confidence interval. The function I ( x ) is an indicator function; it returns 1 if x is true and 0 otherwise. Errork = ∑ (97.5I (L p ,i ≥ η i (α k , Dk )) − 2.5 I (L p ,i ≤ η i (α k , Dk )))(L p ,i − η i (α k , Dk )) (5.3.18) 32 i =1 Errork = ∑ (2.5 I (L p ,i ≥ η i (α k , Dk )) − 97.5(L p ,i ≤ η i (α k , Dk )))(L p ,i − η i (α k , Dk )) 32 (5.3.19) i =1 The optimization of these objective functions was done though a similar random bruteforce technique to that used to find the point estimates. This time, however, the optimization was one-dimensional and was restricted to parameter values that were greater for the upper bound, or lower for the lower bound, of the point estimate. 58 Chapter 6 Chapter 6 Results and Discussion This chapter describes the estimates of the parameters for the ISM, ISSPM and DREM for the seven test surfaces of interest. In the discussion, Section 4, the validity of the three models is considered; this is done through evaluation of the parameter estimates from each model. The estimated values are compared with existing measurements of the test surfaces. The precision of the estimates is also considered as well as the model’s ability to predict the measured data. Before continuing there are a few points about the inversion process that could be helpful to anyone endeavouring to repeat this work. The runtimes for the estimation process was approximately 20 minutes for each estimate for the ISM and ISSPM; this constitutes 10,000 samples. The DREM estimates took approximately 10 minutes each (this includes the time to find the confidence estimates). All of the calculations were completed on a computer using an Intel® coreTM Duo processor T2050 (1.6 GHz, 533 MHz FSB, 2MB L2 Cache). 6.1 Image-surface-model results The point estimates and confidence intervals for parameters of the ISM for each surface are given in Tables A.1-A.7 and Tables A.8-A.13 located in Appendix A. Tables A.1-A.7 are grouped by sample; Tables A.8-A.13 are grouped by frequency. Figure 6.1 displays the parameter estimates with confidence intervals for Z Re , Z im , D and σ 2 . Figures 6.2 through 6.7 display the marginal empirical posterior distributions for the parameters of the ISM inverted on the seven sets of test data. In order to allow for multiple parameters to be plotted together, the mean value, over all frequencies, of the parameter is subtracted and the values are divided by the data range. 59 Chapter 6 Figure 6.1: The parameter estimates and confidence intervals for the ISM for all seven surface samples at the six frequencies (250 Hz, 500 Hz, 1000 Hz, 2000 Hz, 4000 Hz, 8000 Hz). 60 Chapter 6 Figure 6.2: Normalized histograms marginal posterior distributions parameters for the Flat surface. of of the the Figure 6.3: Normalized histograms marginal posterior distributions parameters for the FW surface. of of the the 61 Chapter 6 Figure 6.4: Normalized histograms marginal posterior distributions parameters for the SB surface. of of the the Figure 6.5: Normalized histograms marginal posterior distributions parameters for the BB surface. of of the the 62 Chapter 6 Figure 6.6: Normalized histograms marginal posterior distributions parameters for the Sine1 surface. of of the the Figure 6.7: Normalized histograms marginal posterior distributions parameters for the Sine2 surface. of of the the 63 Chapter 6 Figure 6.8: Normalized histograms marginal posterior distributions parameters for the Soft surface. of of the the 6.2 Image-surface semi-phase model results The point estimates and confidence intervals for parameters of the ISSPM for each surface are given in Tables B.1-B.7 and Tables B.8-B.14 located in Appendix B. Tables B.1 through B.7 are grouped by sample; Tables B.8 though B.13 are grouped by frequency. Figure 6.9 displays the parameter estimates with confidence intervals for Z Re , Z im , D and σ 2 . Figures 6.10 through 6.16 display the marginal empirical posterior distributions of the parameters of the ISM inverted on the seven sets of data. 64 Chapter 6 Figure 6.9: The parameter estimates and confidence intervals for the ISSPM for all seven surface samples at the six frequencies (250 Hz, 500 Hz, 1000 Hz, 2000 Hz, 4000 Hz, 8000 Hz). 65 Chapter 6 Figure 6.10: Normalized histograms of marginal posterior distributions of parameters for the Flat surface. the the Figure 6.11: Normalized histograms of marginal posterior distributions of parameters for the FW surface. the the 66 Chapter 6 Figure 6.12: Normalized histograms of marginal posterior distributions of parameters for the SB surface. the the Figure 6.13: Normalized histograms of marginal posterior distributions of parameters for the BB surface. the the 67 Chapter 6 Figure 6.14: Normalized histograms of marginal posterior distributions of parameters for the Sine1 surface. the the Figure 6.15: Normalized histograms of marginal posterior distributions of parameters for the Sine2 surface. the the 68 Chapter 6 Figure 6.16: Normalized histograms of marginal posterior distributions of parameters for the Soft surface. the the 6.3 Diffuse-reflecting-energy model results Tables C.1 through C.7 and Tables C.8 through C.13, located in Appendix C, display the point estimates and confidence intervals for the parameters of the DREM. Tables C.1 through C.7 are grouped by sample; Tables C.8 though C.13 are grouped by frequency. Figure 6.17 displays the parameter estimates with confidence intervals for α and D . 69 Chapter 6 Figure 6.17: The parameter estimates and confidence intervals for the DREM for all seven surface samples at the six frequencies (250 Hz, 500 Hz, 1000 Hz, 2000 Hz, 4000 Hz, 8000 Hz). 6.4 Discussion In order to validate the models they need to be evaluated on three criteria. These are the ability of the models to fit the observed data, the size of the confidence intervals around the parameter estimates and the agreement of the estimates with previously measured properties of the surfaces. The first of these comparisons discussed is the agreement of the estimates with independently measured values. As noted in Chapter 5 the diffusion coefficient D used in all three models was defined for them; consequently there is no independent way to validate the estimates of it. However some characteristics of the results in Figures 6.1 and 6.9 suggest that the estimates created by inverting the ISM and ISSPM are reasonable. First the flat surface has the lowest diffusion coefficient while the rough surfaces have higher diffusion coefficients. Also the diffusion tends to increase with frequency as expected; this is especially visible for the estimates made by inverting the ISM. 70 Chapter 6 Table 6.1: The measured and estimated impedance values (Z = Z re + iZ im ) for the Flat and Soft surfaces. The measured values are taken from Graves and Hodgson [18]. The impedances of the Flat and Soft surfaces were measured by Graves and Hodgson [18] using the spherical-decoupling method; the results are displayed in Table 6.1. The estimates from the ISM and ISSPM are also displayed; the DREM does not utilize impedance and so it has been excluded. The measured and predicted values have a reasonable agreement for the Flat surface and a poor agreement for the Soft surface. This result is not surprising as the Soft surface was found to have a large diffusion coefficient, and the spherical-decoupling method assumes wholly specular reflection. The average random-incidence absorption α of the surfaces is not affected by diffusion—hence its inclusion in ISO 2003, 17497-1 [20]. Consequently it is helpful to evaluate the predictions of the models by comparing them with measured α values. These values were measured by Bibby [24]. In order for a comparison to be made the impedance values of the ISM and ISSPM must be converted into an average random-incidence absorption coefficient. This is done using the so called Paris formula displayed in Eq. (6.4.1), which is evaluated in Eq. (6.4.2) and Eq. (6.4.3) [1]. A similar conversion is necessary for the absorption parameter estimated by the DREM; this was done using Eq. (6.4.4). π 2 0 4 Z re cos(θ ) π 2 α = 2 ∫ α (θ )sin (2θ )∂θ = 2 ∫ 0 Z cos 2 (θ ) + 2 Z re cos(θ ) + 1 2 sin (2θ )∂θ (6.4.1) 71 Chapter 6 α = Z sin (µ ) Z + cos(2µ ) arctan cos ( ) µ 2 1 + Z cos(µ ) − cos(µ ) log 1 + 2 Z cos(µ ) + Z ( ) sin µ Z ( 8 Z im Z re µ = arctan 2 −α cos(θ ) α sin (2θ )∂θ α =8∫ 2 0 2 −α cos(θ ) + 1 α 2 ) (6.4.2) (6.4.3) π 2 = (6.4.4) 2 2 − α α 2 − α 2 − α − log1 + 2 + +1 − α α 2 − α α 2 − α α 8 2 Table D.1, in the Appendix D, gives the measured and predicted average absorption coefficient values with confidence intervals for the three models. Figure 6.18 displays these values. The ISSPM appears to give the strongest agreement with the measured α values; it is followed by the ISM. The DREM has little or no agreement with the measured values. The agreement of the ISSPM estimates is highest at mid frequencies (500 Hz, 1000 Hz, and 2000 Hz). The confidence intervals for the parameter estimates are, on an aggregate level, smallest for the inversion of the ISSPM; these are smaller than the equivalent intervals for the ISM for approximately 73% of the estimates. The confidence intervals for the DREM are so large as to make the estimates almost uniformly uninformative. The intervals of the α estimates from the ISSPM are similar to those obtained by direct measurement. This means that if only the absorption coefficient is needed there is no penalty (in certainty) to using inverse methods as opposed to direct measurement. The final criterion was the ability of the models to fit the observed data. This is most easily evaluated by comparing the model variances σ 2 . The estimates of σ 2 are, on an aggregate level, smaller for the ISSPM than for the ISM. The DREM does not have an 72 Chapter 6 explicit variance estimate; however there are modal patterns in the white-noise data that the model is simply incapable of predicting; consequently it can be considered to be the worst at fitting the data. σ 2 increases with frequency; this seems to indicate that none of the models is valid at 8000 Hz. Additionally the variance of the model for the BB data set has a higher variance than the other models; this also seems to suggest that there is an upper limit on the of surface roughness for which the inversion process is valid. Figure 6.18: The measured values and estimates with confidence intervals for the acoustical-absorption coeffiecent (Alpha) of seven surfaces. Measured value , ISM ,ISSPM and DREM . 73 Chapter 6 For all three of the criteria the ISSPM appears to be the best; this result is surprising as the ISM is more physically realistic than the ISSPM. The DREM appears to be wholly inadequate for the purpose of estimating the broad-band absorption and diffusion characteristics of surfaces. 74 Chapter 7 Chapter 7 Conclusions and suggested further work 7.1 Conclusions The objective of this research was to develop an inverse-method approach for determining the acoustical diffusion properties of surfaces with unknown impedance. To this end, three models were created. All three models included a diffusion coefficient which represented the proportion of diffuse reflection. Two of the models, the ISM and the ISSPM, considered phase effects and were used to estimate single-frequency characteristics. The third model, DREM, considered only sound energy and estimate broad-band characteristics. Steady-state sound-pressure-level data were collected above seven test surfaces irradiated by a single source in an otherwise free-field. The source was found to have significant directivity in both the azimuth and polar angles. Consequently, the three forward models were modified to include loudspeaker directivity. The three models were inverted using the collected data; the ISM and ISSPM were inverted using Bayesian inversion, and the DREM was inverted using percentile estimation. The estimates produced by the inversion of the three models were evaluated. The evaluation considered their ability to match independent direct measurements, their ability to produce small confidence intervals, and their ability to model the measured data. The approach of the Bayesian inversion of image-source models was found to be an adequate way of determining the acoustical properties of a surface. In particular, the inversion of the ISSPM was shown to find single-frequency parameter estimates that are in agreement with known values. The inversion of the DREM was found to be wholly inadequate at finding broad-band estimates. 75 Chapter 7 7.2 Further work In order to use the inverse methods described in this research with confidence, it would be necessary to conduct a more systematic set of tests on surfaces with known impedance and diffusion properties. It would be advisable to invest considerable time to either ensure that the loudspeaker used has minimal directivity or to measure the loudspeaker directivity accurately. The measurement locations used in this research were selected out of convenience. However, future work should optimize their location depending on the anticipated ranges of the acoustical parameters of the surface to be tested, as well as the frequency range of interest. An earlier version of the ISM was such that exact knowledge of the measurement and loudspeaker locations was unnecessary. These values were optimized in addition to the impedances and diffusion coefficients. This proved computationally cumbersome; consequently later versions of the model did not allow for unknown measurement and loudspeaker locations. However as computer power increases, it may be helpful to revisit this technique, as it is not trivial to accurately measure the three-dimensional coordinates of a point within an anechoic chamber. The model did not include any method for estimating the contribution of back-scattering from the edges of the sample to the sound field. Although this contribution appears to be significantly smaller than the direct and reflected contributions, future work could consider a method of estimating this contribution. All three models considered in this research only considered Lambert’s Law diffuse reflection. It would be of interest to relax this assumption and create a more flexible or general law of reflection, specifically one that would enable energy to be distributed nonuniformly in the azimuth angle. 76 Chapter 7 The ISM and ISSPM where also inverted (using the same methods that were used on the pure-tone data) on the broad-band-noise data. The parameter estimates that were obtained were inconsistent with the measured values as well as the other estimates; although it should be noted that the variance predicted was lower for the ISM on the broad-band data than it was for the pure-tone data. These values are listed in appendix E. It would be useful to revisit this as the broad-band data is significantly easier to obtain. To improve the results the models could be modified to sum the predictions of several frequencies, although this would significantly decrease the speed at which estimates could be produced. Finally future work could consider imposing a frequency-dependent relationship on the impedance values and diffusion coefficients. In such a case meta-parameters that would determine the impedance given the frequency would be estimated. 77 Bibliography [1] H. Kuttruff, Acoustics: An Introduction, (Taylor & Francis, New York, 2007). [2] E. Kinsler and L. Frey, Fundamentals of Acoustics, (John Wiley & Sons, New York, 1967). [3] G. Taraldson, “A note on reflections of spherical waves,” J. Acoust. Soc. Am. 117, 3389-3392 (2005). [4] M. Hodgson, “Evidence of diffuse surface reflection in rooms,” J. Acoust. Soc. Am. 89, 765-771 (1991). [5] H. Kuttruff, "Simulierte Nachschallkurven in Rechteckräumen mit diffusem Schallfeld," Acustica, 25, 333-342 (1971). [6] P. Morse and K. Igram, Theoretical Acoustics, (Princeton University Press, Princeton, 1986). [7] J. Leader, Numerical Analysis and Scientific Computation, (Addison Wesley, New York, 2004). [8] A. Gelman, et al. Bayesian Data Analysis, 2nd ed. (Chapman and Hall, New York, 2004). [9] S. Taherzadeh and K. Attenborough, “Deduction of ground impedance from measurements of excess attenuation spectra,” J. Acoust. Soc. Am. 105, 2239-2242 (1999). [10] S. Kanzler and M. Oelze, “Improved scatterer size estimation using backscatter coefficient measurements with coded excitation and pulse compression,” J. Acoust. Soc. Am. 123, 4599-4607 (2008). [11] T. Poole, G. Frisk, J. Lynch and A. Pierce, “Geoacoustic inversion by mode amplitude perturbation,” J. Acoust. Soc. Am. 123, 667-678 (2008). [12] J. Goutsias and J. Mendel, “Inverse problems in two-dimensional acoustic media: A linear imaging model,” J. Acoust. Soc. Am. 81, 1471-1885 (1987). [13] G. Too, S. Chen and S. Hwang, “Inversion for the acoustical impedance of a wall by using artificial neural network,” Applied Acoustics, 68, 377-389 (2007). [14] S. Dosso and M. Wilmut, “Uncertainty estimation in simultaneous Bayesian tracking and environmental inversion,” J. Acoust. Soc. Am. 124, 82-97 (2008). 78 [15] D. de Vries, N. Joeman and E. Schreurs, “Measurement-based simulation and analysis of scattering structures,” Proc. ForumAcusticum 2147-2150 (2005). [16] Acoustics – Determination of sound absorption coefficient and impedance in impedance tubes – part 2 transfer function method, ISO, 1998, 10534-2. [17] L. de Geetere, Analysis and improvement of the experimental techniques to assess the acoustical reflection properties of boundary surfaces, PhD Thesis, Katholieke University, Leuven (2004). [18] C. Graves and M. Hodgson, “Measurements of the acoustical characteristics of room surfaces using the spherical-decoupling method,” Internal UBC Acoustics and Noise Control Research Group document (2008). [19] Standard Test Method for Sound Absorption and Sound Absorption Coefficients by the Reverberation Room Method, ASTM, c423-07a. [20] Acoustics - Sound-scattering properties of surfaces, ISO 2003, 17497-1. [21] F. Fahy, Foundations of Engineering Acoustics, (Academic Press, San Diego,2001). [22] M. Vorlander and E. Mommertz, “Definition and measurement of random-incidence scattering coefficients,” Applied Acoustics 60, 187-199 (2000). [23] M. Hodgson, “Prediction of noise in industrial workrooms,” J. Acoustic. Soc. Am. 92, 2452-2452 (1992). [24] C. Bibby, “Acoustical characterization of architectural surfaces,” Internal UBC Acoustics and Noise Control Research Group document (2008). 79 Appendices Appendices A, B, C and D contain the tabulated values plotted in Figures 6.1, 6.9, 6.17 and 6.18, respectively. The values presented in Appendices A, B, and C are displayed grouped both by test sample and grouped by frequency. The values in Appendix D are only grouped by test sample. Appendix A Tabulated parameter estimates and confidence intervals obtained through inversion of the ISM. Table A.1: Estimates (Est) and bounds (LB, UB) for DREM with BB data. Table A.2: Estimates (Est) and bounds (LB, UB) for DREM with flat data. 80 Table A.3: Estimates (Est) and bounds (LB, UB) for DREM with FW data. Table A.4: Estimates (Est) and bounds (LB, UB) for DREM with SB data. Table A.5: Estimates (Est) and bounds (LB, UB) for DREM with Sine1 data. Table A.6: Estimates (Est) and bounds (LB, UB) for DREM with Sine2 data. 81 Table A.7: Estimates (Est) and bounds (LB, UB) for DREM with BB data. Table A.8: Estimates (Est) and bounds (LB, UB) for DREM with 250 Hz data. Table A.9: Estimates (Est) and bounds (LB, UB) for DREM with 500 Hz data. 82 Table A.10: Estimates (Est) and bounds (LB, UB) for DREM with 1000 Hz data. Table A.11: Estimates (Est) and bounds (LB, UB) for DREM with 2000 Hz data. Table A.12: Estimates (Est) and bounds (LB, UB) for DREM with 4000 Hz data. Table A.13: Estimates (Est) and bounds (LB, UB) for DREM with 8000 Hz data. 83 Appendix B Tabulated parameter estimates and confidence intervals obtained through inversion of the ISSPM. Table B.1: Estimates (Est) and bounds (LB, UB) for DREM with BB data. Table B.2: Estimates (Est) and bounds (LB, UB) for DREM with flat data. 84 Table B.3: Estimates (Est) and bounds (LB, UB) for DREM with FW data. Table B.4: Estimates (Est) and bounds (LB, UB) for DREM with SB data. Table B.5: Estimates (Est) and bounds (LB, UB) for DREM with Sine1 data. Table B.6: Estimates (Est) and bounds (LB, UB) for DREM with Sine2 data. 85 Table B.7: Estimates (Est) and bounds (LB, UB) for DREM with BB data. Table B.8: Estimates (Est) and bounds (LB, UB) for DREM with 250 Hz data. Table B.9: Estimates (Est) and bounds (LB, UB) for DREM with 500 Hz data. 86 Table B.10: Estimates (Est) and bounds (LB, UB) for DREM with 1000 Hz data. Table B.11: Estimates (Est) and bounds (LB, UB) for DREM with 2000 Hz data. Table B.12: Estimates (Est) and bounds (LB, UB) for DREM with 4000 Hz data. Table B.13: Estimates (Est) and bounds (LB, UB) for DREM with 8000 Hz data. 87 Appendix C Tabulated parameter estimates and confidence intervals obtained through inversion of the DREM. Table C.1: Estimates (Est) and bounds (LB, UB) for DREM with BB data. Table C.2: Estimates (Est) and bounds (LB, UB) for DREM with flat data. Table C.3: Estimates (Est) and bounds (LB, UB) for DREM with FW data. Table C.4: Estimates (Est) and bounds (LB, UB) for DREM with SB data. Table C.5: Estimates (EST) and bounds (LB, UB) for DREM with Sine1 data. Table C.6: Estimates (Est)and bounds (LB, UB) for DREM with Sine2 data. 88 Table C.7: Estimates (Est) and bounds (LB, UB) for DREM with Soft data. Table C.8: Estimates (Est) and bounds (LB, UB) for DREM with 250 Hz data. Table C.9: Estimates (Est) and bounds (LB, UB) for DREM with 500 Hz data. Table C.10: Estimates (Est) and bounds (LB, UB) for DREM with 1000 Hz data. Table C.11: Estimates (Est) and bounds (LB, UB) for DREM with 2000 Hz data. 89 Table C.12: Estimates (Est) and bounds (LB, UB) for DREM with 4000 Hz data. Table C.13: Estimates (Est) and bounds (LB, UB) for DREM with 8000 Hz data. 90 Appendix D Tabulated average absorption estimates and confidence intervals. Table D.1: Displays the estimates (Est) and bounds (LB, UB) of the average absorption of the surfaces for the ISM, ISSPM and DREM models. The values measured by Bibby [23] are also presented. 91 Appendix E Parameter estimates of the ISM and ISSPM inverted on the broadband data. Table E.1: Estimates (Est) and bounds (LB, UB) for ISM with BB data. Table E.2: Estimates (Est) and bounds (LB, UB) for ISM with flat data. 92 Table E.3: Estimates (Est) and bounds (LB, UB) for ISM with FW data. Table E.4: Estimates (Est) and bounds (LB, UB) for ISM with SB data. Table E.5: Estimates (Est) and bounds (LB, UB) for ISM with Sine1 data. Table E.6: Estimates (Est) and bounds (LB, UB) for ISM with Sine2 data. 93 Table E.7: Estimates (Est) and bounds (LB, UB) for ISM with BB data. Table E.8: Estimates (Est) and bounds (LB, UB) for ISSPM with BB data. Table E.9: Estimates (Est) and bounds (LB, UB) for ISSPM with flat data. 94 Table E.10: Estimates (Est) and bounds (LB, UB) for ISSPM with FW data. Table E.11: Estimates (Est) and bounds (LB, UB) for ISSPM with SB data. Table E.12: Estimates (Est) and bounds (LB, UB) for ISSPM with Sine1 data. Table E.13: Estimates (Est) and bounds (LB, UB) for ISSPM with Sine2 data. 95 Table E.14: Estimates (Est) and bounds (LB, UB) for ISSPM with BB data. 96
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Determination of the acoustical diffusion and impedance...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Determination of the acoustical diffusion and impedance characteristics of surfaces by Bayesian inversion… Steininger, Gavin Arthur Mervyn Wyverstone 2009
pdf
Page Metadata
Item Metadata
Title | Determination of the acoustical diffusion and impedance characteristics of surfaces by Bayesian inversion of a modified Helmholtz equation |
Creator |
Steininger, Gavin Arthur Mervyn Wyverstone |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | A surface diffusion coefficient is used in architectural-acoustics to evaluate the effectiveness of diffusing surfaces. The inclusion of the diffusion characteristics is also important for the accuracy of room prediction models. Another important parameter is the absorption or impedance of a surface. In settings with significant low-frequency noise, phase effects are important; consequently impedance values of surfaces are necessary for accurate modeling. A review of existing models for specular and diffuse reflection is made. A new diffusion coefficient is defined and included in three new forward models for predicting the steady-state sound-pressure level above a finite-impedance plane in an otherwise free field. Data are collected for several typical architectural surfaces in an anechoic chamber. Inverse methods are utilized in order to estimate the diffusion coefficient for surfaces given each of the models. This is done without knowledge of the surface impedance, which is simultaneously estimated. The models are compared with each other and with independently measured values of the surface impedance and diffusion. Inversion is found to be a reasonable way of determining the diffusion properties of a surface. |
Extent | 1330705 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-08-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067455 |
URI | http://hdl.handle.net/2429/11995 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2009_fall_steininger_gavin.pdf [ 1.27MB ]
- Metadata
- JSON: 24-1.0067455.json
- JSON-LD: 24-1.0067455-ld.json
- RDF/XML (Pretty): 24-1.0067455-rdf.xml
- RDF/JSON: 24-1.0067455-rdf.json
- Turtle: 24-1.0067455-turtle.txt
- N-Triples: 24-1.0067455-rdf-ntriples.txt
- Original Record: 24-1.0067455-source.json
- Full Text
- 24-1.0067455-fulltext.txt
- Citation
- 24-1.0067455.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0067455/manifest