Proceedings of the Institute of Acoustics

A TIME DOMAIN VERSION OF THE TRANSFER MATRIX METHOD K Steele University of Salford, Salford, UK JA Hargreaves University of Salford, Salford, UK 1 INTRODUCTION

Porous absorbers are ubiquitous in acoustics applications. They are, for example, widely used to control reverberation in rooms and cavities, or to provide attenuation in ducts 1 . Because of this, they have attracted significant research effort to better understand and model their properties, typically informed by the material’s micro-structure. However, the complexity of these models, and the fact that many of the parameters they require are not readily measurable, means that simpler empirically- fitted models are still widely used in Room Acoustics applications. Probably the best known of these is the 1970 model by Delany and Bazley 2 . This is applicable for fibrous materials and estimates acoustic properties based solely on the material’s flow resistivity, which can be readily measured 1 . Like the vast majority of porous material models, this is a frequency- domain model that assumes all quantities are time-harmonic with some frequency. It returns the characteristic impedance and characteristic wavenumber of the material as an equivalent fluid, both of which are complex-valued functions of frequency. To obtain the surface impedance or reflection coefficient for a particular configuration of material layers, the Transfer Matrix Method (TMM) is used. This is also formulated in the frequency-domain based on a time- harmonic assumption 1 . This paper develops a time-domain version of this model, i.e., one that is free from the time-harmonic assumption and that can model the response of the material to transient acoustic signals. A key motivation for this is the increasing popularity of time-domain numerical models of Room Acoustics, such as Finite Difference Time Domain (FDTD) 3 . These algorithms find application in auralisation of spaces but require that material boundary conditions can also be implemented in the time-domain, usually as digital filters. Typically, this has been achieved by fitting a filter to the surface impedance 3 , or reflection coefficient, of the front face of the configuration of porous material. This is the simplest approach, but simple models can struggle with some common configurations. For example, Mondet et al showed that a simple model can reduce uncertainty in measurement 4 , but saw that it struggled to capture the interference pattern that occurs in fairly open materials on a rigid backing. The approach herein overcomes this issue because it uses filters to represent individual processes – propagation through a material layer and reflections at interface – rather than attempt to the final surface impedance, which includes the cumulative effect of several processes. This produces a network of filters, but it is anticipated that each filter can be simpler and lower order without impacting accuracy. The filters fitted are low order Infinite Impulse Response (IIR) filters, which can produce a long impulse response – as is required for strong frequency dependence – with minimal computation and storage cost. This is also a key advantage of filter fitting over inverse fast Fourier transform of frequency domain results. The latter would produce Finite Impulse Response (FIR) filters, which have much higher computation and storage costs for the same frequency domain behaviour. Other time-domain models of acoustic properties of porous materials exist but have been based on different mathematical models of the material. Many of these are focussed on modifying the differential equations that are then simulated numerically using an FDTD scheme. Zhao et al. adopted a similar approach to ours in 2018, using IIR filters to represent bulk modulus and effective density 5 , but their filters are fitted to measurement rather than being derived from a material model. The 2007 method of Wilson et al begins from a time-domain material model based on relaxation processes 6 ,

Vol. 43. Pt. 1. 2021

Proceedings of the Institute of Acoustics

and then discretised these to yield an FIR filter that can be used within an FDTD model of the material. A similar numerical strategy was used by Umnova and Turo in 2009 to implement their time-domain model 7 . Usually in these papers the aim is to produce an FDTD scheme in 2D or 3D where the user can specify how the material properties vary versus position. The method presented here also uses FDTD to model propagation within the material but does so in a simpler way mimicking the TMM. Herein, the aim is to replace the TMM with a network of filters, so the FDTD blocks are 1D ‘pipes’ that only support unidirectional propagation. The most similar approach to ours is probably the 2005 method by Fellah et al for double-layered porous media 8 , which uses transmission and reflectance filters at interfaces between layers with different properties, as is developed herein. Our method differs to that of Fellah et al , however, because it simulates layers in a more modular fashion and is based on Delany and Bazley model 2 , due to its popularity in Room Acoustics. However, this model must be modified before it can be transferred to the time domain. The first of these modifications was by Miki 9 in 1990, who repeated the fit to experimental data while ensuring that the models are physically realisable. But in 2014, Dragna and Blanc-Benon 10 showed that even Miki’s variant is not passive for a rigidly-backed layer, prompting them to propose further modifications. It is their ‘Modified Miki model’ that is the basis of the method presented herein. Section 2 presents the detail of the model and its implementation, including the filter network (section 2.1), the material model and it’s approximation by a digital filter (section 2.2) and incorporation into a 1D unidirectional FDTD model (section 2.3). Section 3 presents results from the method, including validation against the standard frequency-domain TMM. Results are shown for a single porous layer (section 3.1), two layers of porous material (section 3.2) and a single layer of porous material behind which is an air gap (section 3.3). All of these are simulated with a rigid backing, as is typical in Room Acoustics. Finally, section 0 draws conclusions and identifies avenues for further research. 2 MODEL AND IMPLEMENTATION

2.1 Network Development

The starting point for our method is the same as it is for the TMM; the reflection and transmission of a normally-incident plane wave at a change in media properties. This is considered in many Acoustics textbooks and is depicted in Figure 1 (note that throughout this paper the e !"#$ time convention is used). To the left of the dashed line, the medium has characteristic impedance 𝑧 % and wavenumber 𝑘 % . To the right of it, these become 𝑧 & and 𝑘 & respectively. An incident plane wave 𝑝 ' with amplitude 𝐴 ' arrives from the left. Its interaction with the impedance discontinuity causes a reflected wave 𝑝 ( with amplitude 𝐴 ( as well as a transmitted wave 𝑝 ) with amplitude 𝐴 ) that continue rightwards. The derivation of the TMM goes on to use these amplitudes to relate the surface impedance at this interface to the surface impedance at the interface to the next layer 1 . But in the time domain it is often easier to work with transmittance and reflectance instead of surface impedance, because the former two are guaranteed to the causal 11 . Hence, transmission and reflection coefficients at the interface are instead considered. These are defined as 𝑇= 𝐴 ) 𝐴 ' ⁄ and 𝑅= 𝐴 ( 𝐴 ' ⁄ respectively. By the continuity of pressure and particle velocity at the interface it can be shown that:

Figure 1: The transmission / reflection of a normally-incident plane wave at a

change in media properties.

For porous materials, characteristic acoustic impedance is a function of frequency. Since at least one of these layers will be a porous material, 𝑧 % and/or 𝑧 & will be frequency dependent hence 𝑇 and 𝑅 will also both be frequency dependent. In the time-domain they will therefore manifest as filters that have some transmission, or reflection, impulse response. These will be implemented as IIR filters. In our implementation, IIR filters will be fitted to 𝑧 % and 𝑧 & and then filters for 𝑇 and 𝑅 are found from these.

Vol. 43. Pt. 1. 2021

Proceedings of the Institute of Acoustics

The model in Figure 1 is sufficient when the material on the right is of infinite thickness. In this case there is no wave arriving back from material 2. But in the TMM, waves arriving from material 2 need to be accounted for. These occur due to the presence of other layers or a rigid backing behind. To support this, the structure in Figure 1 is summed with a reflected version of itself giving the structure in Figure 2b. Waves travel inwards and outwards so this filter network has an Inward channel (IN), depicted at the top, and an Outward channel (OUT), depicted at the bottom. Each block in the network has two inlet ports and two outlet ports – one of each for both Inward and Outward propagation. Each block is therefore a two-port network, albeit one based on reflection and transmission processes rather than impedance, as is more common in acoustic analysis. The reflection and transmission processes are also distinguished by the subscripts “IN” and “OUT”. T *+ and R *+ match 𝑇 and 𝑅 in eq. 1. T ,-. and R ,-. are found from the same equation simply by switching the role of 𝑧 % and 𝑧 & . The other blocks are as follows. Block a is where waves enter (IN) and leave (OUT) the network from the external air domain. This is where this network would interface to a room acoustic model and is also where the results shown in section 3 have been computed. Block c models the propagation of the wave through the material layer. This is a simple delay if the layer contains air, but there will also be attenuation and dispersion when a porous material is present. Block d is total reflection at a rigid backing, which can be implemented trivially simply as a direct connection. If multiple layers need to be modelled, blocks b and c are repeated to model these layers and interfaces. The main challenge lies in implementing the filters in the boundary and propagation blocks (b & c). How to achieve this is considered in the following two sections. 2.2 Material Model and its Implementation with Digital Filters

Figure 2: A filter network to model a single layer of porous material on a rigid back. Sections: (a) External – input and

output, (b) The boundary block, (c) Propagation through

the material, (d) Rigid backing – total reflection

The Delany and Bazley model 2 remains popular in Room Acoustics but is unsuitable for implementation in the time-domain due to the reasons given in section 1. Hence the closely related ‘Modified Miki Model’ of Dragna and Blanc-Benon is instead chosen 10 . This can be stated as:

Here 𝑧 / and 𝑘 / are respectively the characteristic impedance and wavenumber of the porous material in the layer being considered. 𝑧 0 = 𝜌 0 𝑐 0 and 𝑘 0 = 𝜔𝑐 0 ⁄ are respectively the characteristic impedance and wavenumber of air, where 𝜌 0 is the density of air in kg/m 3 , 𝑐 0 is the sound speed in air in m/s, and 𝜔 is the angular frequency in rads/s. 𝜎 is the flow resistivity of the material in rayls/m. The form of eq. 2 and 3 is different and rather more compact compared to how the Delaney and Bazley and Miki models are normally written. Notably they have been expressed as functions of the Laplace variable 𝑠= j𝜔 because this will be more amenable to filter fitting. Usually there are separate terms for the real and imaginary parts of 𝑧 / and 𝑘 / , but the matching of the exponent values for each of these allows them to be combined via the identity 4,9 (𝑗𝜔) 4 = cos(𝛾𝜋2 ⁄ )𝜔 4 + j sin(𝛾𝜋2 ⁄ )𝜔 4 . This exponent matching across real and imaginary parts was due to Miki 9 , with the further matching of exponents between 𝑧 / and 𝑘 / seen in eq. 2 and 3 due to Dragna and Blanc-Benon 10 .

Vol. 43. Pt. 1. 2021

Proceedings of the Institute of Acoustics

The raising of 𝑠= j𝜔 to an integer power is associated with integer order differentiation (or integration if the exponent is negative). The raising of 𝑠 to a negative fractional power in eq. 2 and 3 is therefore identified as a fractional integral. Dragna and Blanc-Benon 10 give a closed-form expression for this as a convolution kernel, which could be converted to an FIR filter using the method of Wilson et al 6 . Instead we chose to follow the method of Oustaloup et al 12 , which gives IIR filters of controllable order. This fits an integer order polynomial transfer function, which can be implemented with a digital filter, to a fractional power of 𝑠 , as occurs in eq. 2 and 3. The parameters of the method are the order of the fit 𝑁 , a lower frequency limit ω 5 and a transition frequency 𝜔 6 . These were set equivalent to 20 Hz and 20 kHz respectively in all the results that follow. The Oustaloup algorithm produces a set of 2𝑁 poles and 2𝑁 zeroes in the 𝑠 domain that are purely real (non-oscillatory). These were converted into an order 2𝑁+ 1 digital filter using the Impulse Invariant Method. This was chosen because it provides an accurate fit across all frequencies up to Nyquist.

os Imaginary Part os 1 os Real Part ost foo Bon Bor a9 a90 Foal Part

a)

b)

Figure 3: The impedance filter for 𝜎= 5000 rayls/m. a) The frequency response for two filter fits,

os Imaginary Part os 1 os Real Part ost foo Bon Bor a9 a90 Foal Part

one where N=4 and a lower order one, N=1, compared with the analytical impedance. b) The

pole zero plot for N = 4 with inset magnified section showing the edge of the unit circle.

The performance of this method is illustrated by the fit to 𝑧 / shown in Figure 3a. This has been computed for σ = 5000 Rayls/m with a digital sampling frequency of 44.1 kHz (used throughout). The analytical result is the direct evaluation of eq. 2. The frequency and phase response of the fitted digital filter is shown for 𝑁= 1 and 𝑁= 4 . In the middle of the frequency range the fit to both magnitude and phase is very good. At low frequencies the analytical model tends towards infinite amplitude and non- zero phase. This cannot be realised with a digital filter, hence the fit naturally diverges. It is also likely that the analytical target 𝑧 7 given in this range is not realistic anyway, since this family of models (Delany and Bazley, Miki) are known to be inaccurate at low frequencies. Some additional divergence in phase can be seen at high frequencies. This is likely to be due to the Impulse Invariant Method approaching Nyquist, where some aliasing of the transfer function specification is possible. Relatively little difference can be seen between the 𝑁= 1 and 𝑁= 4 fits, supporting the initial premise that lower order filters are satisfactory if they are modelling individual processes. Figure 3b illustrates the numerical challenges when using higher values of 𝑁 . The Oustaloup algorithm will always produce at least one pole close to the unit circle in the vicinity of (1,0); this is necessary to capture the increase in the magnitude of 𝑧 7 towards 0Hz in Figure 3a. However, as 𝑁 is increased the poles and zeros increase in number and start to cluster, ultimately leading to filter instability due to finite numerical precision, limiting practical values of 𝑁 . 2.3 Finite Difference Time Domain (FDTD) Algorithm

It has been seen that the transmission and reflection filters T and R can be found directly from 𝑧 / . Unfortunately, finding the propagation filter M is less straightforward. When a time-harmonic wave propagates through a porous layer of thickness 𝑑 , it’s amplitude and phase are scaled by e 8'9 ! : . Because 𝑘 / is complex, this incorporates both delay and attenuation, and because 𝑘 / is frequency dependent, the process is dispersive. Ultimately, the fact that 𝑘 / is in the exponent of e 8'9 ! : means that a filter for M cannot be directly fitted. An FDTD model of the material layer is required instead.

Vol. 43. Pt. 1. 2021

Proceedings of the Institute of Acoustics

It is common in the TMM to assume that porous materials are locally reacting and that all propagation is surface-normal. Because of this, only a 1D FDTD scheme was required. Much of the other literature on porous material modelling with FDTD uses a staggered grid scheme 5,6 , but here a simpler second- order central difference scheme was used 3 . Notably, this avoids a dependency of the scheme on 𝑧 / that is not expected given that the frequency domain form of M above only involves 𝑘 / . It states that pressure 𝑝 ;

Here 𝑋 denotes the spatial node spacing and 𝑇= 1 f = ⁄ is the time-step duration, where f = is the digital sampling frequency in Hz. The pressure signal is injected into the FDTD scheme at node 𝑖= 0 and read out at node 𝑖= 𝐼 , where 𝐼× 𝑋= 𝑑 , the material thickness. To avoid reflections from the grid truncation, Mur’s absorbing boundary condition 13 was applied at this location. This was used to update a ‘ghost node’ at node 𝑖= 𝐼+ 1 , which then fed into the standard update equation (eq. 4). The update equation for this is:

The stencil for this is shown in orange in Figure 4. It updates the ghost node (red) for time step 𝑛 , after which eq. 4 (green and blue) is used to update all nodes 0 < 𝑖≤𝐼 for time step 𝑛+ 1 . Finally, 𝑝 0

𝐼

𝐼−1

𝐼−2

Figure 4: Implementation of

Mur’s absorbing boundary condition via a ghost node

Results were generated by feeding a transient pulse into the P *+% port of the network and comparing the resulting output signal at port P ,-.% to it. This output signal is either presented directly, to allow the transient response of the material configuration to be seen, or both are fast Fourier transformed and then divided to compute the frequency domain reflection coefficient, to allow comparison with the standard TMM. The excitation pulse was a Ricker (Mexican hat) wavelet with a pulse width of 0.05ms and a centre time of 0.4ms. This combination gave minimal onset discontinuity at 0s and provided sufficient energy over the frequency range of interest (20 Hz – 10 kHz). It is shown in Figure 5. For most of the results that follow, 𝑁= 4 was used for the filter fit. The exception to this was the double layer cases in section 3.2, where 𝑁= 1 was required to prevent instability arising in the boundary filters.

Figure 5: The excitation pulse fed into the network

Vol. 43. Pt. 1. 2021

& Amplitude & Os 15 Time (s) 25 35 10"

Proceedings of the Institute of Acoustics

3.1 Single Layer Results

a)

The first results were computed for a single-layer configuration, using the network exactly as shown in Figure 2. Three different materials were simulated, with flow resistivity 𝜎 set to 5,000, 15,000 and 25,000 rayls/m. All samples were 10cm thick and on a rigid backing. The reflected pressure signals output by the network (in response to the excitation pulse in Figure 5) are shown in Figure 6. All three absorbers exhibit an immediate reflection from the front face at 0.4ms, due to the impedance mismatch between the porous material ( 𝑧 / ) and air ( 𝑧 0 ). As expected, the reflection’s amplitude is larger for greater values of flow resistivity due to the larger impedance mismatch. This reflection is produced solely by the reflectance filter R *+ in the boundary block (Figure 2b). The pulse then enters the material propagation block (Figure 2c), as is implemented by the FDTD scheme in section 2.3. The delay and attenuation that occur as it propagates through this is dictated by the characteristic wavenumber of the material 𝑘 / and the layer thickness 𝑑 . For case a ( 𝜎= 5000 Rayls/m), a second pulse can be clearly seen at 1ms. This is the reflection from the rigid backing. It appears so clearly because this material is quite ‘open’ and does not attenuate the wave very significantly as it propagates. It arrives with a delay of roughly 6ms, consistent with 2𝑑/𝑐 0 . In case b ( 𝜎= 15000 Rayls/m), a far smaller, and slightly later, 2 nd pulse can be observed, but for case c ( 𝜎= 25000 Rayls/m) none is visible. These denser materials are much more effective at attenuating the wave as it propagates through them, but their effectiveness as porous absorbers is compromised by the increased amplitude of their front reflections, due to the large impedance mismatch they present. To validate the new time domain model against the standard TMM, normal-incidence reflection and absorption coefficients were plotted for all examples. An example result is shown in Figure 7 for the 10cm single layer with 𝜎= 5000 Rayls/m. Agreement is seen to be excellent in the majority of the frequency range. Deviation at low frequencies is mostly in magnitude and is likely due to the lower frequency limit of the material filter fitting, as discussed in section 2.2. At high frequencies, differences are seen in phase and notch frequencies in reflection coefficient. This suggests an inaccuracy in wave speed in the FDTD scheme, which is expected since these are known to exhibit dispersion error 3 .

b)

c)

Figure 6: Reflected pressure signals from

three 10cm single layer absorbers

a)

b)

2 02 2 oa a Bo JX g SR 0 os 1 i 2 28 9 88 Time (6) g — 10cm of 15000 rayis 2 oa = a Bo & 01 0.08 +1 18 2 28 @ a5 Time (6) go 106m of 25000 rayis m"? Zo Bo & 0 os 1 8 2 25 9 65 4 Time () 10%

c)

Figure 7: Frequency domain validation results for a single 10cm layer with 𝜎 = 5000 rayls/m,

showing absorption coefficient (a) and reflection coefficient magnitude (b) and phase (c).

Vol. 43. Pt. 1. 2021

Proceedings of the Institute of Acoustics

3.2 Double Layer Results

a)

The results from three double layer absorbers are shown in Figure 8. For these, the boundary layer and propagation blocks (b&c) in Figure 2 are repeated. All the absorbers consist of a top 5cm layer with 𝜎= 5000 rayls/m. Behind this are layers of differing flow resistivity and/or thickness, followed by a rigid backing. Comparison between the first two designs shows the effect of an impedance discontinuity between the layers. The back layers are both 5cm thick, but the first (a) is very dense with 𝜎= 50000 rayls/m and the second (b) is more open with 𝜎= 10000 rayls/m. The first result (a) has strong similarities with the rigidly back results in Figure 6a. Despite the second layer being a porous material, it is so dense that almost 100% reflection is occurring. Figure 6b shows that using a less dense second layer mitigates this, but now a later reflection around 1ms occurs due to the rigid backing. Finally, the third design in Figure 6c doubles the thickness of the second layer, delaying and attenuating this backing reflection. 3.3 Air Gap Results

Tine) 108

b)

c)

Figure 8: Reflected pressure signals from

three double layer absorbers

a)

Finally, three absorbers with a single porous layer and an air gap were modelled, a common and cost-effective design. This again required duplication of the boundary and propagation blocks in Figure 2 b&c, but the air layer does not require an FDTD model. Instead, a simple delay filter is sufficient to capture its behaviour. In all cases the front porous material layer is 5cm thick. Designs a and c show the effect of its density; 𝜎= 5000 rayls/m for a and 𝜎= 25000 rayls/m for c. Both have a 5cm air gap. Again, it is seen that the impedance mismatch controls the initial reflection. But compared to Figure 6 a&c, the reflections from the rigid backing are stronger since half the material has been removed and replaced with air. The effect of air gap thickness is shown by comparing the responses from designs a and b. These have the same porous material but differing air gap; 5cm and 15cm respectively. This bigger gap spreads out the reflections in time due to the extra propagation delay it allows. What is especially interesting is the third pulse in b that appears around 2.5ms. This involves a second order reflection in the air cavity. This pulse has propagated: i) through the material, ii) through the air gap reflecting off the rigid backing, iii) reflected off the rear of the porous material, iv) travelled through the air gap and reflected off the rigid backing a second time, before v) finally propagating back through the material and emerging into the room from the front of the device. This pulse has low amplitude, so is unlikely to be very important. But nonetheless, these results illustrate how time domain reflection responses can reveal details that are not obvious in results from time harmonic models.

b)

c)

Figure 9: Reflected pressure signals from

three absorbers with air gaps

Amplitude 18 2 25 3 35 Time () 5m of 5000 rayis m and Sem of 10000 rayis Amplitude Os 18 2 25 3 35 4 Time () “a Amplitude u —— Sem of 5000 rayis m'" and 10cm of 10000 rayis m”? is 2 +28 @ 35 4 Time (6) 108

Vol. 43. Pt. 1. 2021

Proceedings of the Institute of Acoustics

4 CONCLUSIONS

A time domain version of the Transfer Matrix Method (TMM) has been presented. This comprises a modular network of filters that intuitively represents transmission and reflection at material interfaces, and propagation through layers. Filters were fitted to the Modified Miki Model of Dragna and Blanc- Benon 10 using the method of Oustaloup et al 12 . The propagation ‘filter’ used a 1D FDTD scheme. Results were validated against the standard frequency domain TMM by fast Fourier transform. Agreement was excellent in the majority of the frequency range, though divergence occurs at low and high frequencies due to filter fitting limitations (and possibly FDTD dispersion error). Issues with filter stability were experienced for higher order filter fits, but these did not appear to be necessary. The validity constraints of the Delany and Bazley model 1 would also apply to this model in real world use. Notably, the results have shown how time domain models can reveal details that are not obvious when time harmonic models are used. Many Room Acoustic designers are concerned about temporal detail, so methods that show this are likely to be of interest, even though they are mathematically equivalent to existing frequency domain approaches. Further work could involve adding support for non-normal incidence, comparing to other time domain material models or attempting to eliminate the need for the FDTD model by direct inverse Laplace transform of the propagation operator. 5 REFERENCES

1. Cox, T. J., & D’Antonio, P. (2009). Acoustic Absorbers and Diffusers: Theory, Design and

Application (2nd ed.). Taylor & Francis. 2. Delany, M. E., & Bazley, E. N. (1970). Acoustical properties of fibrous absorbent materials.

Applied Acoustics , 3 (2), 105–116. 3. Kowalczyk, K., & Walstijn, M. van. (2008). Modeling Frequency-Dependent Boundaries as Digital

Impedance Filters in FDTD and K-DWM Room Acoustics Simulations. Journal of the Audio Engineering Society , 56 (7/8), 569–583. 4. Mondet, B., Brunskog, J., Jeong, C. H., & Rindel, J. H. (2020). From absorption to impedance:

Enhancing boundary conditions in room acoustic simulations. Applied Acoustics , 157 . 5. Zhao, J., Bao, M., Wang, X., Lee, H., & Sakamoto, S. (2018). An equivalent fluid model based finite-difference time-domain algorithm for sound propagation in porous material with rigid frame. The Journal of the Acoustical Society of America , 143 (1), 130. 6. Wilson, D. K., Ostashev, V. E., Collier, S. L., Symons, N. P., Aldridge, D. F., & Marlin, D. H.

(2007). Time-domain calculations of sound interactions with outdoor ground surfaces. Applied Acoustics , 68 (2), 7. Umnova, O., & Turo, D. (2009). Time domain formulation of the equivalent fluid model for rigid

porous media. The Journal of the Acoustical Society of America , 125 (4), 1860–1863. 8. Fellah, Z. E. A., Wirgin, A., Fellah, M., Sebaa, N., Depollier, C., & Lauriks, W. (2005). A time- domain model of transient acoustic wave propagation in double-layered porous media. The Journal of the Acoustical Society of America , 118 (2), 661. 9. Miki, Y. (1990). Acoustical properties of porous materials-Modifications of Delany-Bazley

models. J. Acoust. Soc. Jpn.(E) , 11 (1), 19–24. 10. Dragna, D., & Blanc-Benon, P. (2014). Physically Admissible Impedance Models for Time-

Domain Computations of Outdoor Sound Propagation. Acta Acustica United with Acustica , 100 (3), 401–410. 11. Hargreaves, J. A. (2011). Simulating transient scattering from obstacles with frequency-

dependent surface impedance. Proceedings of Forum Acusticum . 12. Oustaloup, A., Levron, F., Mathieu, B., & Nanot, F. M. (2000). Frequency-band complex

noninteger differentiator: characterization and synthesis. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications , 47 (1), 25–39. 13. Mur, G. (1981). Absorbing Boundary Conditions for the Finite-Difference Approximation of the

Time-Domain Electromagnetic-Field Equations. IEEE Transactions on Electromagnetic Compatibility . https://doi.org/10.1109/TEMC.1981.303970

Vol. 43. Pt. 1. 2021