Thin layers of porous media are widely adopted in sound absorption and noise control applications due to their compact arrangement. Particularly, porous media with low flow resistivity exhibit complex, non-local reaction behavior. Therefore, sound field prediction above these media is computationally challenging. This is due to singularities and branch points in the expression for the reflection coefficient. This paper introduces a framework based on the direct discrete complex image method to analyze the sound field above a rigid-backed, non-locally reacting porous sample. In contrast to traditional complex image applications, the proposed framework avoids the extraction of the quasi-static term and the poles from the reflection coefficient. Instead, the reflection coefficient—including its singularities—is directly approximated in terms of a series of complex exponentials, whose coefficients are determined with the matrix pencil method. This study analyzes the sound field above two test samples made of melamine and rockwool. The sound field computation is efficient and accurate in the near- and in the far-field. Moreover, predicted specific impedances agree well with experimental *in situ* impedance measurements. The proposed framework serves for more applications including object detection in multilayered porous grounds or sound propagation prediction in layered atmosphere.

## I. INTRODUCTION

Porous materials are an effective tool in noise and reverberation control. Sound propagation inside and above such materials generally depends on their complex pore structure. Simplified locally reactive material behavior can be considered for samples with high flow resistivity, i.e., *σ* ≥ 30 000 Ns/m^{4}.^{1} However, many materials in engineering applications, e.g., foams or fibrous materials, are characterized by low flow resistivity and therefore tend to possess non-local reaction behavior. Consequently, the sample's surface acoustic impedance and the sound absorption coefficient depend on the angle of incidence of an impinging sound wave.^{1,2}

There are well-known methods to calculate the sound field caused by a harmonic monopole above a locally reactive surface such as the model of Nobile and Hayek^{3} and the pole subtraction method proposed by Chien and Soroka.^{4} Di and Gilbert^{5} presented an exact formulation for local and extended reaction semi-infinite grounds based on the Laplace transform of the reflection coefficient. Li *et al.*^{6} stated that a non-locally reactive sample of finite thickness with rigid backing can be approximated by a semi-infinite ground with extended reaction behavior, if the sample is thick enough with respect to its intrinsic wave number. However, the need for lightweight and compact design often demands for thin layers. Consequently, non-locally reactive material behavior has to be considered in these cases.

Li *et al.*^{6} and Dragonetti and Romano^{7} tried to overcome the dependence of the surface impedance of a rigid-backed layer with the incidence angle by modifying the well-known models for locally reactive material behavior.^{3–5} The literature shows that this is well accomplished by the modified models, but all of them lack accuracy at low frequencies and when the source is placed near to the sample under investigation, since the spherical character of the radiated wave dominates under these circumstances.^{7,8} Yet, such a source-receiver setup is of particular interest for *in situ* surface impedance estimation.^{1,7,9}

Dragonetti and Romano^{7} showed that the formulation of Allard *et al.*^{10} is the most accurate model for sound field computations above a non-locally reactive sample backed by a rigid wall, since it is based on the velocity potential fields above and inside the porous sample. In this case, the diffractive integral term can be evaluated by using direct numerical integration. However, the integrand has to undergo multiple mathematical transformations prior to the numerical integration,^{1,10} but even after mathematical manipulation, singularities remain in the integrand.^{1} These singularities are notably poles of the reflection coefficient, yielding surface waves in the physical space, which primarily contribute to the far-field.^{11–13} Accuracy and stability of numerical integration in the far-field can be increased by providing information on the integrand's singularities to the integration routine.^{1} Yet, finding these singularities numerically and extracting them is computationally challenging and expensive. This holds in particular for multilayered media, as the singularities cannot be determined analytically.^{14}

Alternatively, the fast field program^{15} can be used to calculate the sound field in the far-field. Li^{16} introduced the integrated exact solution (IES), which combines the direct numerical integration solution for the near-field with the fast field program solution for the far-field. By this, a numerically exact but computationally expensive solution can be obtained.^{17,18} This makes it impractical for optimization tasks, such as the inference of the surface impedance.^{1}

For this reason, recent research concentrates on efficient approximation methods for Sommerfeld integrals. Based on the work of Brekhovskikh and Godin,^{19} Liu and Li^{20,21} presented an improved integration scheme along the steepest descent path to obtain the sound field above a porous sample of finite thickness backed by a rigid wall or a surface impedance. This method is effective for long distances between source and receiver as well as at high frequencies. Further, Li and Tao^{22} proposed an accurate asymptotic solution of the diffraction integral using the modified saddle point method. By this, no singularity extraction from the diffraction integral kernel is required. Yet, this only holds for model setups considering normal incidence. Recently, Wang *et al.*^{23} and Li and Wang^{24} extended the application of the saddle point method in conjunction with the pole subtraction method to estimate the sound field radiated from a point source moving at constant speed above extended-reacting surfaces.

In the field of electromagnetics, the discrete complex image method (DCIM) has attracted interest as an efficient alternative to approximate Sommerfeld integrals^{11,17,25} by computing the integral kernel in terms of a series of damped exponentials.^{26} In the field of acoustics, however, few researchers have addressed the DCIM for sound field analysis with applications limited to either locally reactive behavior^{27–30,32} or extended reaction behavior of semi-infinite grounds.^{27,34,35} The initial work of Li and White^{27} focused on the principle of complex images to evaluate the sound field above an impedance plane and a semi-infinite ground with extended reaction behavior. Fawcett^{34} analyzed the acoustic scattering from a seabed, which was either modeled as a half-space with extended reaction behavior or an acousto-elastic half-space. The DCIM served to approximate the Green's function for these half-space problems and a non-linear least squares optimization was used to find the coefficients of the complex images. Fawcett^{35} extended his work to calculate the Green's function modeling a bounded wave guide and making it applicable to layered half-spaces. Ochmann and Piscoya^{29,30} demonstrated the incorporation of the DCIM for modeling an impedance boundary condition in boundary element simulations. They analyzed the sound field above a locally reactive plane and showed the advantages of the DCIM over conventional steepest descent approaches. Zheng *et al.*^{31} incorporated the principle of complex images into a fast multipole boundary element method to efficiently solve large-scale locally reactive half-space problems. Brandão and Lenzi^{32} reported on the advantages of the DCIM in *in situ* surface impedance inference for a locally reactive material. A general overview of the DCIM and its use for acoustic half-space problems can be found in the recent book of Ochmann and Piscoya.^{33}

All these complex image approaches rely on the traditional DCIM, for which the quasi-static term, i.e., the limiting value of the reflection coefficient for the intrinsic wave number in air tending to zero, $limk0\u21920R\u0303$ (contributing to the near-field) and the poles of the reflection coefficient (contributing to the far-field) have to be identified and extracted from the Sommerfeld kernel.^{14} This method is impractical since these poles, which lead to surface waves in the physical space, are difficult to determine analytically or numerically. Further, errors occur in the near-field, due to the singularity at the origin of the Hankel function.^{11,13,17} Another drawback of the traditional DCIM is the use of polynomial methods, in particular Prony's method,^{36,37} to determine the coefficients for the complex images. Prony's method shows unstable behavior with an increasing number of complex images.^{26} Yet, a multilayered ground demands for an increased number of complex images due to its challenging reflection behavior,^{35} which also holds for a rigid-backed, non-locally reacting sample of finite thickness.

This work introduces a robust and computationally efficient framework to analyze the sound field above a non-locally reacting layer with rigid backing. It is based on the direct DCIM^{14} and uses the matrix pencil method^{26,38} (MPM) to calculate the coefficients for the complex images. To the best of the authors' knowledge, this is the first time the DCIM has successfully been applied to predict the sound field above a rigid-backed, non-locally reacting sample of finite thickness.

This paper is structured as follows: Sec. II presents the general equations for the sound field above a non-locally reactive sample of infinite radial extent backed by a rigid wall. Further, the fundamentals of the direct DCIM are demonstrated. Particular attention is given to the MPM, which is the backbone of the direct DCIM. Results of numerical studies and experimental validation of the proposed method are displayed in Sec. III. Section IV concludes this article and outlines further applications of the proposed framework.

## II. THEORETICAL BACKGROUND

The outlined theory here considers the sound field above a rigid-backed, non-locally reactive sample of finite thickness and infinite radial extent (Sec. II A) and the procedure to estimate the sound field with use of the direct DCIM (Sec. II B). Moreover, the process to calculate the coefficients of the complex images by using the MPM is demonstrated (Sec. II C).

### A. Sound field above a rigid-backed non-locally reactive sample

A radiating monopole above a partly reflective plane of infinite radial extent is a widely used model to estimate the sound absorption properties of porous layers *in situ*.^{1,7,27,39–41} Figure 1 shows the model setup for a porous layer of finite thickness *d* above a rigid backing. The source is located at $rs=(xs,ys,zs)$ in the upper, homogeneous half-space made of air and radiates with angular frequency *ω*. Assuming a harmonic time dependence of $e\u2212i\omega t$ with $i=\u22121$, the sound pressure at a prescribed receiver position $r=(xr,yr,zr)$ can be expressed by^{27}

given a source strength *S*. A unit source strength of *S =* 1 kg/s^{2} is assumed for all further analyses. For the sake of convenience, the problem is reformulated in cylindrical coordinates^{19} and the radial distance between source and receiver $r=(xr\u2212xs)2+(yr\u2212ys)2$ is introduced.

The specific half-space Green's function *G* can be stated as^{27}

with the first term representing the direct contribution from the source to the receiver as a free-field Green's function with the direct distance $D1=r2+(zr\u2212zs)2$ between source and receiver. The intrinsic wave number of air is $k0=\omega /c0$, with *c*_{0} being the speed of sound in air. The zeroth order Bessel function is represented by *J*_{0}. The second term in Eq. (2) accounts for the reflected wave. With the presented method, there is no need to separate this term into a reflective (quasi-static) term and a diffractive integral term. Instead, the integrand remains in its spectral domain representation,^{42} which is expressed in terms of the vertical wave number in air

and the radial wave number *k _{r}*. The vertical and the radial wave number can be interpreted as a decomposition of the wave number vector in vertical and radial coordinate directions. Formally, this manifests in the Sommerfeld identity

^{43}

which represents the decomposition of a spatial spherical wave into a plane wave propagating in vertical direction and a cylidrical wave propagating in radial direction in the wave number domain. The negative (positive) sign in the integrand is selected if the source is located above (below) the receiver.^{27}

The contribution and mathematical complexity of the reflected wave depends on the plane wave reflection coefficient $R\u0303(kr)$. For a single layer of finite thickness *d* backed by a rigid wall, it can be expressed as^{6}

In this equation, *ρ*_{0} and *ρ*_{1} represent the density of air and the complex-valued density of the porous medium, respectively. The vertical wave number in the porous medium is

which depends on the complex-valued intrinsic wave number *k*_{1} of the porous medium. Introducing the reflection coefficient [Eq. (5)] in the integral term in Eq. (2), there is no closed-form solution for the Green's function.^{1} The integrand is mathematically challenging, due to its poles and zeros,^{1} the numerous discontinuities in the tangent function and the highly oscillating behavior of the Bessel function.^{17} This raises the need to manipulate the integrand to make it numerically tractable with conventional methods, e.g., numerical integration.

### B. Direct discrete complex image method

The direct DCIM as introduced by Yuan *et al.*^{14} promises to overcome the need to extract the quasi-static term and the poles from the reflection coefficient $R\u0303(kr)$ [Eq. (5)], which is directly approximated by a series of damped exponentials

with expansion coefficients *a _{n}* and

*b*. Inserting the approximated reflection coefficient [Eq. (7)] into the integrand in Eq. (2) yields

_{n}The integral in the reflected wave term can be integrated exactly using the Sommerfeld identity [Eq. (4)], which leads to a closed-form expression for the specific half-space Green's function

Each summand represents a spherical wave with complex-valued amplitude *a _{n}*, which is radiated from an image source located at the complex-valued distance $Dn=r2+(zr+zs+ibn)2$ from the receiver. As the coordinates of the complex images solely differ in their values for

*b*, all complex images are located at the same real-valued position, but can be interpreted as a line source along the additional imaginary axis. This makes this formulation inherently convenient to include the cylindrical wave nature of surface waves as a contribution from a line source given the total number of image sources

_{n}*N*in Eq. (9) is sufficient.

^{14,44}

Equation (9) serves for an efficient evaluation of the Green's function. Remark that the only approximation error occurs in the approximation of the reflection coefficient [Eq. (7)], while all other operations are analytically exact. Since the reflection coefficient is directly approximated without any terms extracted, Eqs. (7)–(9) differ from their complements in traditional DCIM formulations.^{27–30,32,34,35}

In order to compute the coefficients *a _{n}* and

*b*, it is necessary to sample the reflection coefficient in a real-valued variable.

_{n}^{14}In traditional DCIM applications, a sampling path in the complex $Kz,0$-plane

^{27}

is used. The reflection coefficient has to be evenly sampled in the real-valued variable *t* along this path to allow for the use of Prony's method, while the end of this path *T*_{0} should be chosen larger than the branch points to avoid singularities. This is not necessary for the direct DCIM, as the MPM (Sec. II C) can estimate the behavior beyond *T*_{0} given sufficient information on the reflection coefficient's singularities is included in the sampling path.^{14} For this reason, a modified sampling path

with the scaling parameter *γ* has been proposed by Sarkar *et al.*^{14} for the direct DCIM.

The reflection coefficient should be sampled evenly along the modified path such that all possible singularities (in the multilayered case) are equally accounted for during the sampling and the reconstruction process using the MPM. Figure 2 qualitatively compares two variants of the modified sampling path for the direct DCIM, *C*_{2} for *γ* = 1 and *C*_{12} for $0<\gamma <1$ in Eq. (11), with the path used in traditional DCIM applications $C\u2032$ [Eq. (10)]. The modified path runs parallel to the imaginary axis, which simplifies the relationship between the variable *t* and the vertical wave number in air $Kz,0$ compared to traditional DCIM applications.

Figure 3 shows the reflection coefficient of the melamine foam sample (see Table I) at 1000 Hz sampled along the modified path *C*_{2}, i.e., $\gamma =1.0$ [Fig. 3(a)] and along a shifted modified path *C*_{12} with $\gamma =0.05$ [Fig. 3(b)]. Smaller values of the scaling parameter *γ* move the modified path closer to the imaginary axis, which leads to more variation in the sampled reflection coefficient. This correlates to more information on the reflection coefficient's singularities in the samples used in the MPM, which is pivotal to account for the contribution of surface waves in particular in the far-field. By this, accurate sound field estimations can be achieved for higher distances between source and receiver. In other words, a smaller value of *γ* extends the accuracy of the direct DCIM in the far-field.^{14} However, shifting the modified path closer to the imaginary axis comes at the cost that the sampling path does no longer start at the origin of the complex *k _{r}*-plane [Fig. 2(b)]. Consequently, a compensation sampling path along the real axis—

*C*

_{11}—has to be introduced if $\gamma <1$. Therefore, the procedure to calculate the coefficients for the complex images has to be slightly modified.

^{14}The modified calculation procedure includes the following steps:

. | Flow resistivity . | Porosity . | Tortuosity . | Viscous char. length . | Thermal char. length . |
---|---|---|---|---|---|

. | σ (Ns/m^{4})
. | $\varphi $ (—) . | $\alpha \u221e$ (—) . | Λ (m $10\u22126$) . | $\Lambda \u2032$ (m $10\u22126$) . |

Melamine foam | 9207 | 0.99 | 1.0 | 300 | 600 |

Rockwool | 11 817 | 0.97 | 1.0 | 113 | 122 |

. | Flow resistivity . | Porosity . | Tortuosity . | Viscous char. length . | Thermal char. length . |
---|---|---|---|---|---|

. | σ (Ns/m^{4})
. | $\varphi $ (—) . | $\alpha \u221e$ (—) . | Λ (m $10\u22126$) . | $\Lambda \u2032$ (m $10\u22126$) . |

Melamine foam | 9207 | 0.99 | 1.0 | 300 | 600 |

Rockwool | 11 817 | 0.97 | 1.0 | 113 | 122 |

Select the scaling parameter

*γ*according to the desired maximum radial distance between source and receiver.- Generate
*M*evenly distanced samples of the reflection coefficient along the modified path*C*_{12}and approximate the reflection coefficient with a first series expansion(12)$R\u0303(kr)|C12\u2248\u2211n=1Nane\u2212bnKz,0|C12=\u2211n=1NAneBnt|C12$in terms of the coefficients in the transformed space

*A*and_{n}*B*._{n} Set the primary truncation threshold

*tol*and apply the MPM (Sec. II C) to the samples of Eq. (12) to get the coefficients*a*,_{n}*b*and the number of primary complex image sources_{n}*N*; note that*N*<*M*.- Generate $M\u2032$ samples of the residual along the compensation path
*C*_{11}(Fig. 2) and approximate the residual in terms of a second series expansion(13)$R\u0303(kr)|C11\u2212\u2211n=1Nane\u2212bnKz,0|C11\u2248\u2211n\u2032=1N\u2032a\u2032n\u2032e\u2212b\u2032n\u2032Kz,0|C11.$ Set the compensation truncation threshold $tol\u2032$ and apply the MPM to the samples of Eq. (13) to calculate the coefficients $a\u2032n\u2032,b\u2032n\u2032$ and the number of compensation complex images $N\u2032$; note that $N\u2032<M\u2032$.

- Hence, for $0<\gamma <1$, the expression for the Green's function [Eq. (9)] changes to(14)$G(r,rs)=eik0D1D1+\u2211n=1Naneik0DnDn+\u2211n\u2032=1N\u2032a\u2032n\u2032eik0Dn\u2032Dn\u2032,$
with $Dn\u2032=r2+(z+zS+ib\u2032n\u2032)2$.

In the following, Eq. (14) is denoted as the extended range solution, in which the total number of complex images is the sum of the number of primary image sources *N* (step 3) and compensation image sources $N\u2032$ (step 5).

For applications in the near-field, it suffices to use the modified path *C*_{2} (Fig. 2), i.e., *γ* = 1 in Eq. (11). In this case, no compensation path is required and the two-path procedure reduces to an approximation of the reflection coefficient along a single path. By this means, only steps 1–3 are relevant and the Green's function can then be efficiently evaluated using Eq. (9). We recommend the direct DCIM procedure for near-field analyses in the range of $k0r<10$, as proposed by Li and White.^{27}

### C. Matrix pencil method

The MPM is the backbone of the direct DCIM, since it enables the direct approximation of the reflection coefficient [Eq. (7)]. Compared to Prony's method, the MPM is a single-step process to calculate the expansion coefficients *a _{n}* and

*b*(step 3) and in a second function call $an\u2032$ and $bn\u2032$ (step 5) in case the extended range solution is used.

_{n}^{26}In the following, the procedure for the MPM is outlined for process step 3, but it holds in similar manner for step 5 in which all primary variables have to be replaced by their compensation counterparts, i.e., $M\u2192M\u2032$ and $tol\u2192tol\u2032$.

The MPM takes *M* evenly distanced samples of the reflection coefficient generated along the modified path $R\u0303(kr)|C12=s[j]$ with $j=1,2,\u2026,M$ (or $M\u2032$ evenly distanced samples of the residual along the compensation path *C*_{11}) as an input. In an initial step, these samples are arranged in a rectangular Hankel matrix of dimension $(M\u2212L)\xd7(L+1)$

with *L* denoted as the pencil parameter. We set the pencil parameter to $L=(M/2)\u22121$, according to the recommendation of Sarkar *et al.*^{26,45} Next, a truncated singular value decomposition of the Hankel matrix $Y$ is performed to determine the dominant singular values fulfilling the condition $\sigma j>max(\sigma j)\xd7tol$. Only the columns of the diagonal matrix containing these dominant singular values in descending order and their corresponding right singular vectors are kept in the filtered matrices $\Sigma f$ and $Vf$, while all other singular values are discarded. This truncation is an essential step in the direct DCIM, since the number of dominant singular values defines the number of primary complex images *N* (or compensation images $N\u2032$), i.e., $N=rank(\Sigma f)$. Then, the matrix pencil $Y2\u2212\lambda Y1$ is formed with the submatrices $Y1=U\Sigma fV1fH$ and $Y2=U\Sigma fV2fH$. The matrices $V1fH$ and $V2fH$ result by eliminating the last row and the first row of $Vf$, respectively.^{26} The matrix $U$ contains all left singular vectors of the full singular value decomposition. Finally, the dominant signal poles *z _{n}* are obtained by solving the general eigenvalue problem

where $Y1\u2020$ is the Moore-Penrose pseudoinverse of the first submatrix, which is used to guarantee a stable and efficient eigenvalue analysis, as both submatrices are ill-conditioned.^{38} The damping coefficients in the deformed space

follow directly from the signal poles and the sampling interval $\Delta t$. The unknown amplitudes in the deformed space *A _{n}* follow from the least squares problem

Once the coefficients in the transformed space *A _{n}* and

*B*are found, retransformation yields the original coefficients for the complex images

_{n}Note that the set of coefficients *a _{n}*,

*b*(and $a\u2032n\u2032,b\u2032n\u2032$ as well) is independent of the source and of the receiver position, as stated by Li

_{n}*et al.*

^{27}By this, any source-receiver combination can efficiently be evaluated with the same set of coefficients for the complex images. Moreover, the direct DCIM enables efficient evaluation of the particle velocity with the same set of coefficients by using the linearized Euler equation. Thus, the specific impedance can be evaluated at any receiver position above the porous sample without any additional effort.

Algorithm 1 summarizes the proposed framework to calculate the specific half-space Green's function based on the direct DCIM. It should be mentioned that the presented method is only valid, when source and receiver are located above the porous ground. However, the presented framework can be easily extended to analyze sound propagation within the porous layer. For this purpose, the derivative of the Sommerfeld identity is required.^{14}

SET scaling parameter γ according to desired maximum radial distance |

GENERATE M samples of reflection coefficient along modified path C_{12} [Eq. (11)] |

CALL MPM to compute coefficients a, _{n}b and number of primary complex images _{n}N |

IF $0<\gamma <1$ THEN$\u2192$ extended range solution |

GENERATE $M\u2032$ samples of residual [Eq. (13)] along compensation path C_{11} |

CALL MPM to compute coefficients $a\u2032n\u2032,b\u2032n\u2032$ and number of compensation complex images $N\u2032$ |

COMPUTE Green's function [Eq. (14)] |

ELSE $\u2192$ near-field solution |

COMPUTE Green's function [Eq. (9)] |

END IF |

SET scaling parameter γ according to desired maximum radial distance |

GENERATE M samples of reflection coefficient along modified path C_{12} [Eq. (11)] |

CALL MPM to compute coefficients a, _{n}b and number of primary complex images _{n}N |

IF $0<\gamma <1$ THEN$\u2192$ extended range solution |

GENERATE $M\u2032$ samples of residual [Eq. (13)] along compensation path C_{11} |

CALL MPM to compute coefficients $a\u2032n\u2032,b\u2032n\u2032$ and number of compensation complex images $N\u2032$ |

COMPUTE Green's function [Eq. (14)] |

ELSE $\u2192$ near-field solution |

COMPUTE Green's function [Eq. (9)] |

END IF |

## III. RESULTS

The proposed direct DCIM framework is applied to three test scenarios. First, the accuracy and efficiency of the direct DCIM are numerically validated for normal incidence. Second, numerical prediction of the specific impedance is validated against *in situ* impedance measurements in a hemi-anechoic chamber. Third, the extended range solution of the proposed method is analyzed numerically. Computations are conducted on a machine with eight physical cores operating at 3.30 GHz and with 96 Gb RAM.

Each scenario treats the sound field above a rigid-backed, non-locally reactive sample of finite thickness with source and receiver located at a height of $zs=0.33$ m and $zr=0.015$ m above the surface of the sample. The density and the sound speed of air are $\rho 0=1.21$ kg/m^{3} and $c0=343$ m/s. A melamine foam sample and a rockwool sample are considered. Rigid frame porous material behavior is assumed and the intrinsic wave number of the sample *k*_{1} and its complex density *ρ*_{1} are calculated with the Johnson-Champoux-Allard model.^{46} Since this paper focuses on the capability of the direct DCIM to predict the sound field above the samples, this simpler material model is chosen, even though more elaborate material models have proven a better fit with respect to impedance tube measurements.^{1} Table I lists the parameters characterizing the considered samples, whose measurement is detailed in the work of Brandão *et al.*^{1}

### A. Normal incidence

The first two analyses deal with normal incidence of the radiated sound wave. In these cases, there is no radial offset between source and receiver, hence *r =* 0. Therefore, the near-field solution process for the direct DCIM with *γ* = 1 is used.

#### 1. Numerical validation

A melamine foam sample of thickness *d =* 40 mm is modeled with material properties as listed in Table I. The solution with direct numerical integration (DNI) as presented by Brandão *et al.*^{1} is used as reference. The reflection coefficient is sampled along the modified path *C*_{2} with a sampling rate of $\Delta t=0.1$ in the interval $0<t<T0=7.5$. The truncation tolerance to obtain the dominant singular values is set to $tol=10\u221210$. The sound field is analyzed in a frequency range from 100 to 10 000 Hz in steps of 1 Hz.

Figure 4 shows the dimensionless sound pressure, i.e., the total sound pressure *p* normalized against the free-field sound pressure $pff=S\u2009(eik0D1/D1)$ above the melamine foam sample obtained with the direct DCIM (dotted) and the reference solution (dashed) across the considered frequency range. Its real part is displayed in Fig. 4(a) and the imaginary part in Fig. 4(b). Reflection at the sample's surface and non-local acoustic material behavior lead to interference effects above the sample. The direct DCIM solution yields nearly an identical sound pressure compared to the reference solution across the whole frequency range. The MPM reveals a maximum number of *N =* 11 primary image sources to obtain the presented results. No compensation sources are calculated, since we use the near-field solution (see Algorithm 1) in case of normal incidence (*r =* 0). The proposed direct DCIM framework shows high efficiency in terms of computation time compared to the reference solution (DNI). Given one source-receiver setup, the sound pressure across the whole frequency range as displayed in Fig. 4 is obtained in 25.2 s using the direct DCIM approach compared to 750 s using direct numerical integration.

Since, the direct DCIM solution and the reference solution are visually indistinguishable (Fig. 4), the accuracy of the proposed framework is further assessed with the normalized error in the total sound pressure

individually for real and imaginary part at a discrete frequency *f _{q}*. The absolute error in the sound pressure is normalized with the root mean square across the entire frequency range of either the real part $\u211c{pDNI}RMS=1Nf\u2211q=1Nf(\u211c{pDNI(fq)})2$ or the imaginary part $\u2111{pDNI}RMS=1Nf\u2211q=1Nf(\u2111{pDNI(fq)})2$ of the total sound pressure, with $Nf$ representing the total number of discrete frequencies considered.

The normalized error is displayed in Fig. 5(a) for the real part and in Fig. 5(b) for the imaginary part. Both errors show a slight increase for frequencies above 2500 Hz, but remain below a value of $10\u22125$ for the entire frequency range. Actually, the DCIM solution oscillates around the reference solution, which is a consequence of the approximation in terms of a series of complex exponentials. Yet, the accuracy of the proposed method is validated and its efficiency makes it a viable tool in optimization problems, e.g., in surface impedance estimation or material characterization.

#### 2. Experimental validation

Experimental validation is conducted for a 40 mm thick melamine foam sample with a sample size of 0.6 m × 0.6 m and a 50 mm thick rockwool sample with a sample size of 1.2 m × 0.6 m with properties listed in Table I. Numerical estimations of the specific impedance *Z _{m}* above the samples are compared to

*in situ*impedance measurements from previous studies.

^{1,9}A free-field calibrated pressure-particle velocity sensor was used for measurements in a hemi-anechoic chamber with a volume of 100 m

^{3}and a cut-off frequency of 150 Hz for the absorbing wedges. 30 sequential measurements with high reproducibility have been conducted for each sample under normal incidence given a source height of $zs=0.33$ m and a receiver height of $zs=0.015$ m. The reader is referred to the works of Brandão

*et al.*

^{1,9}for detailed information on the measurement setup. Regarding the direct DCIM framework, all parameters remain the same as in the previous section, i.e., $\gamma =1.0,\u2009\Delta t=0.1,\u2009T0=7.5$, and $tol=10\u221210$.

Figure 6 shows the real (left column) and imaginary parts (right column) of the normalized specific impedance determined above the melamine foam sample [Figs. 6(a) and 6(b)] and the rockwool sample [Figs. 6(c) and 6(d)]. The measured impedances are displayed as the solid black line, the calculated impedances using the direct DCIM framework are displayed in gray, and the reference solution (DNI) is plotted as a black dashed line. The specific impedance shows similar trends for both materials. The numerical models predict the specific impedance well for the melamine foam sample. However, they underestimate the real part for the rockwool sample at low frequencies and show lesser fluctuation than the experimental data at higher frequencies. These deviations can be attributed difficulties in the experimental determination of the macroscopic parameters (Table I) and to the assumption of rigid-frame material behavior assumed with the Johnson-Champoux-Allard model.^{46} Better predictions, in particular in the low frequency regime might be obtained by using more elaborate material models, such as the Johnson-Champoux-Allard-Lafarge model,^{47} which has been identified as the most appropriate model for the materials analyzed,^{1} or the Biot-Allard model.^{48} These models include low-frequency limits of the dynamic material parameters or consider elastic behavior of the porous structure. Moreover, edge effects affect the measured data, in particular in the low frequency regime, since the measured samples are of finite size.^{49,50} In contrast, the numerical setup considers samples of infinite radial extent. Hence, deviations in the low frequency regime also stem from model errors. Again, for the specific impedance, the direct DCIM yields indistinguishable results compared to the numerical integration, which requires integration on both pressure and particle velocity and, therefore, can become computationally expensive. The presented framework shows good performance in estimating the specific impedance and it is inherently designed to serve as the core of a robust impedance inference algorithm.

### B. Extended range analysis

The proposed direct DCIM framework intrinsically also allows for fast and accurate sound field predictions in the far-field, as the coefficients for the complex images are independent of the source and the receiver position. The user only has to adjust the scaling parameter $\gamma <1$ customized to the associated application and the method follows the steps outlined in Algorithm 1 to obtain an accurate approximation of the Green's function.

The sound field above the rigid-backed melamine foam sample of 40 mm thickness is analyzed for varying radial distance between source and receiver. As in the previous analyses, the source and the receiver are located at a height of $zs=0.33$ m and $zr=0.015$ m above the surface of the layer. Studies are conducted at 100, 1000, 5000, and 10 000 Hz. Here, the IES according to Li^{16} serves as a reference to determine the accuracy of the direct DCIM approximation. The end of the modified path for the direct DCIM framework is set to $T0=8.5$ to include more information on the singularities. The sampling rate and the truncation threshold on the modified path *C*_{12} are $\Delta t=0.01$ and $tol=10\u221210$. Note that the results for 100 Hz are obtained with a sampling rate of $\Delta t=0.1$ on *C*_{12}, which allows one to analyze the effect of the number of samples on the calculation time. The residual [Eq. (13)] is sampled with $M\u2032=50$ samples along the compensation path *C*_{11} and a truncation threshold of $tol\u2032=10\u22123$ is prescribed for the MPM on the samples of the residual. Only a few compensation sources are required to approximate the behavior of the residual, since the residual does not contain a lot of variation as the reconstruction with the primary image sources does account for most of the reflection coefficient's variation.^{14} For this reason, the number of samples of the residual $M\u2032$ and the truncation threshold $tol\u2032$ can be set lower than for the primary step (step 3).

Figure 7 shows the estimated excess attenuation $EA=20\u2009log10|p/pff|$ for a variation of the radial distance between source and receiver. The results for 100 Hz are displayed in Fig. 7(a), for 1000 Hz in Fig. 7(b), for 5000 Hz in Fig. 7(c), and for 10 000 Hz in Fig. 7(d). Sound field analysis is conducted for scaling parameters $\gamma =1.0$ (blue), $\gamma =0.4$ (magenta), and $\gamma =0.05$ (orange). Additionally, the IES reference solution is displayed (black dashed). The excess attenuation shows nearly no variation for small radial distances, higher variation with dips at 5000 and 10 000 Hz in an intermediate range and decreases rapidly for larger distances between source and receiver. The direct DCIM framework yields accurate results in the near-field for all considered scaling parameters. The single path solution ($\gamma =1.0$) deteriorates the earliest from the reference solution at each frequency. The solution with $\gamma =0.4$ yields accurate results up to higher radial distances, but only the solution with scaling parameter $\gamma =0.05$ gives accurate results across the entire range. Table II lists the number of primary and compensation complex images *N* and $N\u2032$, the maximum radial distance $max(r)5%$, when the direct DCIM solutions deteriorate 5% relative to the reference solution and the total calculation time *τ* for all direct DCIM solutions displayed in Fig. 7.

f in Hz
. | γ
. | N
. | $N\u2032$ . | $max(r)5%$ . | τ in s
. |
---|---|---|---|---|---|

100 | 1.0 | 10 | n.a. | $11.6\u2009\lambda $ | 0.254 |

0.4 | 12 | 4 | $85.2\u2009\lambda $ | 3.95 | |

0.05 | 20 | 7 | $>408.8\u2009\lambda $ | 5.66 | |

1000 | 1.0 | 11 | n.a. | $7.6\u2009\lambda $ | 2.87 |

0.4 | 12 | 4 | $62.4\u2009\lambda $ | 6.09 | |

0.05 | 15 | 5 | $>2283.6\u2009\lambda $ | 685.0 | |

5000 | 1.0 | 12 | n.a. | $55.6\u2009\lambda $ | 2.92 |

0.4 | 13 | 4 | $111.6\u2009\lambda $ | 6.24 | |

0.05 | 15 | 4 | $>7301.9\u2009\lambda $ | 674.0 | |

10000 | 1.0 | 13 | n.a. | $54.0\u2009\lambda $ | 2.9 |

0.4 | 14 | 4 | $190.8\u2009\lambda $ | 6.15 | |

0.05 | 17 | 4 | $>8057.9\u2009\lambda $ | 710.0 |

f in Hz
. | γ
. | N
. | $N\u2032$ . | $max(r)5%$ . | τ in s
. |
---|---|---|---|---|---|

100 | 1.0 | 10 | n.a. | $11.6\u2009\lambda $ | 0.254 |

0.4 | 12 | 4 | $85.2\u2009\lambda $ | 3.95 | |

0.05 | 20 | 7 | $>408.8\u2009\lambda $ | 5.66 | |

1000 | 1.0 | 11 | n.a. | $7.6\u2009\lambda $ | 2.87 |

0.4 | 12 | 4 | $62.4\u2009\lambda $ | 6.09 | |

0.05 | 15 | 5 | $>2283.6\u2009\lambda $ | 685.0 | |

5000 | 1.0 | 12 | n.a. | $55.6\u2009\lambda $ | 2.92 |

0.4 | 13 | 4 | $111.6\u2009\lambda $ | 6.24 | |

0.05 | 15 | 4 | $>7301.9\u2009\lambda $ | 674.0 | |

10000 | 1.0 | 13 | n.a. | $54.0\u2009\lambda $ | 2.9 |

0.4 | 14 | 4 | $190.8\u2009\lambda $ | 6.15 | |

0.05 | 17 | 4 | $>8057.9\u2009\lambda $ | 710.0 |

The proposed framework yields accurate results for larger distances between source and receiver if a smaller scaling parameter is chosen. Reducing the scaling parameter means that more primary image sources *N* are required to account for the increased variation in the sampled reflection coefficient. Further, smaller values of the scaling parameter lead to larger sampling ranges, $0\u2264t\u2264T0/\gamma $, according to Eq. (11). Since the sampling rate $\Delta t$ is held constant, this results in an increased number of samples *M* for the reflection coefficient along the modified path *C*_{12}. Therefore, the dominant singular values have to be determined from a larger Hankel matrix and a larger least squares problem has to be solved for the amplitudes, which explains increased computation times. Coarser sampling can significantly reduce the calculation time (see results for 100 Hz in Table II), but this comes at the cost of a decreased maximum radial distance. The number of compensation image sources $N\u2032$ is rather unaffected by the choice of the scaling parameter and only few compensation images ($N\u2032\u22645$ in most cases) suffice for accurate predictions. This has similarly been observed by Yuan *et al.*^{14} Note that the frequency hardly affects the calculation time. By adjusting the scaling parameter, the proposed framework can be used for sound propagation predictions into the far-field, e.g., in underwater acoustics above a layered seabed or sound propagation in layered atmosphere. Further, we stress that the computation of *a _{n}*,

*b*and $a\u2032n\u2032,\u2009b\u2032n\u2032$ is independent of the position of sources and receivers and that the evaluation of Eqs. (9) or (14) is straightforward. The operator only has to choose the scaling parameter according to a desired maximum radial distance. The values for $max(r)5%$ in Table II serve as an orientation for such tasks.

_{n}## IV. CONCLUSION

This article introduces the direct discrete complex image method to the field of acoustics. Compared to traditional complex image applications, neither the quasi-static term nor the poles have to be extracted from the reflection coefficient prior to the application of the complex image method. This avoids the cumbersome identification and extraction of singularities of the reflection coefficient. Instead, the reflection coefficient is directly sampled along a modified sampling path and the coefficients of the complex images are calculated with the matrix pencil method, which includes an automated selection of the number of image sources. With these attributes, the direct discrete complex image method allows for a simplified application compared to traditional discrete complex image method applications. Further, the proposed method allows for the analysis of computationally challenging scenarios, such as the sound field above a rigid-backed, non-locally reacting sample of infinite radial extent, which is validated against experimental *in situ* impedance measurements for a melamine foam and a rockwool sample. Accurate sound field predictions are obtained in the near- and in the far-field. The scaling parameter in the modified sampling path, which the user has to prescribe *a priori*, decides up to what distance the approximation yields accurate results. Smaller values for the scaling parameter enable the analysis to larger distances between source and receiver. The proposed method extends the use of the complex image method in the field of acoustics and makes it a viable alternative for computationally efficient predictions of the sound field above layered media, since the reflection coefficient is directly encountered without the need to extract any singularities. Future work includes the incorporation of the direct discrete complex image method in a robust impedance inference framework and half-space sound field analyses using the boundary element method.