In order to incorporate a directive sound source into acoustic simulation using the finite-difference time-domain method (FDTD), this paper proposes an optimization-based method to estimate the initial value which approximates a desired directional pattern after propagation. The proposed method explicitly considers a discretized FDTD scheme and optimizes the initial value directly in the time domain so that every effect of the discretization error of FDTD, including numerical dispersion, is taken into account. It is also able to consider a frequency-wise directivity by integrating the Fourier transform into the optimization procedure, even though the estimated result is defined in the time domain. After the optimization, the obtained result can be utilized in any acoustic simulation based on the same FDTD scheme without modification because the result is represented as the initial value to be propagated and no additional procedure is required.

Acoustic simulation has been widely studied for predicting acoustical phenomena. Most methods are classified into the following two categories: geometrical- and wave-based methods. In geometrical models, including ray tracing1 and the image source method,2 sound is assumed to propagate as a ray, i.e., all of the wave properties of sound are ignored. In contrast, wave-based methods numerically compute an approximate solution to the wave equation which explains the wave properties well. Owing to the computational difficulty of approximating the wave equation, a lot of methodology has been studied, for example, the finite-difference time-domain method (FDTD),3 the digital waveguide mesh (DWM),4 the finite element method (FEM),5 the boundary element method (BEM),6 and the pseudo-spectral time-domain method (PSTD).7 Among these wave-based methods, FDTD is one of the more popular due to its capability of obtaining impulse responses straightforwardly and simplicity of the calculation scheme. Therefore, a number of studies on acoustic simulation using FDTD have been conducted.3,8–22

Among the many factors affecting the simulation accuracy including the boundary condition and numerical dispersion, this paper focuses on the directivity of a sound source in FDTD. As the directivities of sound sources vary depending on the types of the sources,23–26 acoustic simulation working with directive sound sources is important in order to obtain a realistic result. However, most FDTD simulations have been performed with the omnidirectional sound source because of its easiness of implementation in FDTD. The source directivity is decided by how to excite the sound field, and simple excitation methods, including a method using the spatial Gaussian pulse as the initial value, usually result in the omnidirectional pattern. While some researchers have considered directive sound sources for specific simulation methods,27–29 there is little research on the implementation of directive sound sources in FDTD.30–33 

In the existing literature, two types of directivity modeling have been mainly considered: superposition of point sound sources28,30,31 and expansion into spherical harmonics or multipoles.29,32,33 Superposition of point sound sources is a popular methodology in sound field processing34 because any sound field can be approximated by a set of point sources. The same idea has been considered in acoustic simulation to approximate a directive sound source, which utilizes a set of point sound sources distributed around the position of the directive source.28,30 It modeled the source in the frequency domain so that the sound field emitted from a point source admits the closed-form expression. Then, the problem of approximating the directivity is reduced to the problem of finding an appropriate set of coefficients of the predefined point sources. This problem can be solved through the least squares method30 (or its variant)28 easily. However, since these methods consider the approximation problem in the frequency domain, the obtained result must be converted into the time domain by some method whose choice may degrade the final accuracy. For circumventing such situation, some point-source-based methods consider the directivity directly in the time domain.31 Yet, it is not easy to obtain a closed-form expression for a complicated directivity in this way, and thus only a simple source pattern such as cardioid has been considered.

The spherical-harmonics- or multipole-based methods consider the directivity in the time domain, which allows direct implementation into FDTD. The spherical-harmonics-based method realizes a directive sound source by multiplying the spherical harmonics, which approximate the desired directivity, to the spatial Gaussian pulse.32 It can be regarded as a higher-order extension of the usual Gaussian pulse, and therefore the implementation of the approximated directional pattern is straightforward after the directivity is expanded into the coefficients of spherical harmonics. The multipole-based method considers spatial derivatives of a bandlimited Dirac delta function, or bandlimited multipoles, to express a desired directivity by their linear combination.33 It can design directivity having a desired spatial frequency response. However, these time-domain methods did not consider (time-directional) frequency-wise directivity because the spherical harmonics and multipoles serve as a basis of the directional pattern for all frequencies simultaneously. To date, there seems no method which can both control the (time-directional) frequency-wise directivity and be directly implemented in the time domain without the domain conversion.

One subject that the previous research on directive-sound-source simulation has not addressed is the effect of the discretization of FDTD. Although every numerical simulation includes a discretization error, models of the sound source in the conventional methods have not fully taken the effect of discretization into account.30,32 They instead theoretically consider sound propagation and directivity as if there is no numerical error. Then, the obtained directivity should not be optimal in the sense of FDTD. That is, the discretization effect specific to FDTD has partly been considered or totally been ignored.

In this paper, as a first step toward the optimal directive-source representation in FDTD, a method for imposing directivity into FDTD by optimizing the initial value is proposed. The proposed method explicitly considers a discretized FDTD scheme and optimizes the directivity directly in the time domain. Since the sound propagation is represented through the FDTD scheme, both spatial and time discretization effects are included in the optimization. Hence, the obtained result by the proposed method is optimal in the sense of FDTD. Its implementation into the FDTD scheme after the optimization is effortless because the result is an estimated initial value which gives the desired directional pattern after propagation. In addition, the proposed formulation allows us to incorporate frequency-wise directivity, although the estimated result is in the time domain. To the best of our knowledge, except our conference proceeding on the preliminary result,35 a fully optimization-based representation of a directive sound source, which takes every effect of discretization of FDTD into account, is realized for the first time in the literature.

The rest of the paper is organized as follows. First, FDTD is briefly reviewed in Sec. II, and its matrix form is also summarized within Sec. II. Then, based upon the matrix representation of FDTD, the proposed method is introduced as an optimization procedure in Sec. III. Numerical examples are shown in Sec. IV to illustrate the ability of the proposed method, and some discussions on the characteristics of the estimated results are made. Finally, this paper is concluded in Sec. V.

In this paper, only the two-dimensional field is considered for the sake of a simpler explanation. Though extending the proposed method to the three-dimensional field is straightforward, three-dimensional consideration is omitted due to the limitation of space.

In an acoustic field, sound propagation can be represented in the time domain by the wave equation,

(1c22t2)p(r,t)=0,
(1)

where p is sound pressure, c is sound speed, r=(x,y) is a position vector, t is time, and is the Laplace operator. The solution to this equation is uniquely determined by the appropriate initial and boundary values. While the boundary value is related to the reflection at a wall, the initial value (condition of p at t = 0) is related to the sound source. Therefore, we focus on the initial value in this paper for realizing directive sound sources. In order to simplify the problem formulation, a forcing term is not considered as in Eq. (1).

FDTD is a numerical analysis technique for calculating an approximate solution to the differential equations. In FDTD, derivatives are approximated by finite differences to generate a discrete recursive scheme which is often called the update equation. A numerical solution is calculated by iterating this scheme in the time domain, and therefore it is suitable for acoustic simulation because approximate impulse responses can be directly obtained. Various FDTD schemes have been proposed for solving Eq. (1) or the set of equations equivalent to Eq. (1).

The sound pressure can be represented by its values on the Cartesian grid as pi,j[n]p(hXi,hXj,hTn), where hX and hT are the grid sizes of the spatial and time variables, respectively, i,j are the spatial indices, and n is the time index. FDTD approximates the derivatives in Eq. (1) by the finite-difference operators,

Dttpi,j[n]=c2dpi,j[n],
(2)

where Dtt and d are the finite-difference approximations of the second-order time derivative and the Laplace operator, respectively. These operators are arbitrary approximation of the corresponding derivatives, which allows infinite variations of FDTD schemes. While any approximation for Dtt and d can be utilized, a simple example of the time difference is considered here:

Dttpi,j[n]=(1/hT)2(pi,j[n+1]2pi,j[n]+pi,j[n1]).
(3)

Then, an update equation of FDTD can be written as

pi,j[n+1]=c2dpi,j[n]+2pi,j[n]pi,j[n1],
(4)

where all terms except pi,j[n+1] are gathered on the right hand side to obtain a time recursive update rule. For the spatial difference d, one simple example is

dpi,j[n]=(1/hX)2(pi+1,j[n]2pi,j[n]+pi1,j[n]+pi,j+1[n]2pi,j[n]+pi,j1[n]).
(5)

Note that, for calculating the spatial difference dpi,j[n], the boundary condition is necessary because the values of p at i ± 1 and j ± 1 are required. In this paper, the size of the acoustic field is assumed to be large enough so that a pulse emitted from the sound source does not reach to the boundary within the time steps considered in the optimization procedure. That is, any boundary condition is acceptable for the proposed method by setting a large enough boundary.

For simplifying the notation, this paper utilizes a matrix form of FDTD which is explained in this section. Note that, however, while the proposed method will be introduced in the matrix form, it can be calculated without explicitly constructing the huge matrix. Such matrix-free implementation is postponed until Sec. III D.

Let Np be the total number of grid points of a sound field, and p[n]Np be a vector whose elements consist of all sound pressure values pi,j[n] at every grid point. Then, Eq. (4) for all grid points can be summarized as

p[n+1]=c2mtxp[n]+2p[n]p[n1],
(6)

where mtxNp×Np is the matrix representing d, i.e., the elements of mtxp[n] are dpi,j[n] which is calculated as in Eq. (5). By defining a state vector

ζ[n]=[p[n+1]T,p[n]T]T2Np.
(7)

Eq. (6) can be simplified to

ζ[n+1]=Φζ[n]=Φ2ζ[n1]==Φn+1ζ[0],
(8)

where Φ2Np×2Np is a square block matrix,

Φ=[c2mtx+2INpINpINpONp],
(9)

INp and ONp are the identity and zero matrices, respectively, whose sizes are Np × Np. Although this representation is specific to the update equation in Eq. (4), any other FDTD scheme can be represented as Eq. (8) in a similar manner.

This equation indicates that the nth time step of the state vector ζ[n] (sound fields) is calculated as Φnζ[0]. That is, the initial value ζ[0] decides the sound field after propagation by Φn. Therefore, a directive sound source can be implemented by appropriately choosing ζ[0], which requires inverse estimation of the initial value from the desired directivity as illustrated in Fig. 1.

FIG. 1.

(Color online) Schematic of inverse estimation of the initial value from given instances of sound propagation.

FIG. 1.

(Color online) Schematic of inverse estimation of the initial value from given instances of sound propagation.

Close modal

For implementing a directive sound source into FDTD, a method for estimating the initial value, which approximates the desired directivity after propagation, is proposed in this section. The unique feature of the proposed method is that the effect of the discretization of an FDTD scheme is explicitly considered through Φ, while the conventional methods only consider the idealized (theoretical) model of sound propagation. This feature should be interesting because it might be useful to compare FDTD schemes from a directivity point of view, which is left as a future work. Remind that, even though the proposed method is explained through a specific scheme in Eqs. (4) and (6) for the sake of easiness, any FDTD scheme can be expressed by the matrix Φ and vector ζ[n] (see  Appendix A for another example), and therefore the proposed method is not restricted to a specific FDTD scheme.

In FDTD, a spatially concentrated pulse, usually the spatial Gaussian pulse, is utilized as the initial value. For considering a sound source similar to such a condition, the initial value ζ[0] is assumed to be non-zero only inside a small region illustrated as the shaded (blue) disk in Fig. 2(a). That is, the initial value outside the shaded region is set to 0. This assumption is formulated by defining a vector of the initial value inside the region,

ξ=[pin[1]T,pin[0]T]T2Nin,
(10)

whose length is much smaller than ζ[n] in Eq. (7) defined on the whole domain, i.e., NinNp. The initial value ζ[0] is recovered from ξ by padding zero to the outside of the region, which is shortly written as

ζ[0]=Eξ,
(11)

where E{0,1}2Np×2Nin is an embedding matrix whose column has only single 1 and other entries are 0.

FIG. 2.

(Color online) Setting of initial value and observation points. (a) The initial value is considered in the shaded (blue) region, while the observation points for calculating a directional pattern are set on a (green) circle outside the shaded (blue) region. (b) Frequency-wise directivity is illustrated as a two-dimensional diagram, where the vertical axis represents angle from −π to π, and the horizontal axis is frequency from 0 to the sampling frequency fs. The colorbar represents magnitude of the Fourier spectrum at each direction.

FIG. 2.

(Color online) Setting of initial value and observation points. (a) The initial value is considered in the shaded (blue) region, while the observation points for calculating a directional pattern are set on a (green) circle outside the shaded (blue) region. (b) Frequency-wise directivity is illustrated as a two-dimensional diagram, where the vertical axis represents angle from −π to π, and the horizontal axis is frequency from 0 to the sampling frequency fs. The colorbar represents magnitude of the Fourier spectrum at each direction.

Close modal

After propagating the initial value of n steps by FDTD as ΦnEξ(=Φnζ[0]), some directivity is expected to appear if the initial value is corresponding to a directive sound source. Here, the directivity is considered as the sound pressure at the observation points, which are set outside the (blue) region and placed on the (green) dotted circle whose radius is robs, as illustrated in Fig. 2(a). The sound pressure at the observation points in nth time step is represented by an Nobs-dimensional vector pobs[n]Nobs whose length is also much smaller than ζ[n](NobsNp). The relation between these vectors is written similar to Eq. (11) as

pobs[n]=Sζ[n],
(12)

where S{0,1}Nobs×2Np is a selection matrix whose row has only single 1 and other entries are 0. That is, S collects the sound pressure only at the observation points and ignores other entries of ζ[n] (see Fig. 12 in  Appendix B for an illustrative explanation).

Based on the above notations, the relation between the initial value inside the region ξ and the sound pressure at the observation points pobs[n] is expressed as

pobs[n]=SΦnEξ,
(13)

which is an equation for a single time step n. On the other hand, directivity is not decided by a single time instance but by direction-wise energy which is averaged along time. Therefore, the sound pressure pobs[n] must be considered in multiple time steps. To do so, the vectors pobs[n] for multiple n are concatenated to form an nNobs-dimensional vector pobs as

pobs=[pobs[1]T,pobs[2]T,,pobs[n]T]TnNobs,
(14)

which enables us to formulate the simultaneous linear equations involving all time steps at once:

pobs=SdiagΦrepEξ,
(15)

where Sdiag=blkdiag(S,S,,S){0,1}nNobs×2nNp is the block diagonal matrix consisting of the selection matrices, blkdiag(A1,,An) is the operation generating a block diagonal matrix from inputted matrices A1,,An, and Φrep=[ΦT,(Φ2)T,,(Φn)T]T2nNp×2Np is the matrix generating the sound fields of all time steps from the initial value Eξ(=ζ[0]).

This equation indicates that the initial value inside the region ξ can be obtained by, say, the least squares method, if a desired sound pressure at the observation points is given for every time step as pobs. However, directivity is often considered as the direction-wise energy whose information for each time step is not available. That is, only the energy at each observation point pobsengy=m=1n|pobs[m]|2 is available as data. Therefore, the initial value ξ must be estimated from a given direction-wise non-negative data, such as pobsengy, so that the propagated sound SdiagΦrepEξ approximates it.

For summarizing the sound pressure at the observation points pobs in the time domain, the Fourier transform, which sums up the time indices, is utilized. Let fNDFTNobs×nNobs be a matrix, constructed from Nobs discrete Fourier transform matrices (with zero padding) whose size is NDFT × n, such that the product fpobs acts as the Fourier transform at each observation point and results in a vector consisting of Nobs spectra whose length is NDFT (here, NDFT ≥ n is the length of the Fourier transform that can be increased by the zero padding). By applying this matrix to both sides of Eq. (15), the initial value ξ is related to the direction-wise Fourier spectra at the observation points fpobs as

fpobs=fSdiagΦrepEξ,
(16)

where the magnitude of the direction-wise Fourier spectra |fpobs| is schematically illustrated in Fig. 2(b). Since the squared magnitude of these spectra represents the direction- and frequency-wise energy, controlling the energy of fSdiagΦrepEξ results in the control of the frequency-wise directivity.

To control the direction- and frequency-wise energy, a weighting matrix WNDFTNobs×NDFTNobs, a diagonal matrix with non-negative diagonal entries, is introduced. Based on these notations, the estimation problem of the initial value from a given directivity is formulated as the following energy minimization problem:

minξWfSdiagΦrepEξ22s.t.ξ22=1,
(17)

where q22=m|qm|2 is the squared Euclidean norm. The weighting matrix imposes frequency-wise directivity by setting a larger weight to the direction which is not desired, and vice versa. The norm constraint ξ22=1 is required to avoid the trivial solution ξ=0. By solving this problem, the initial value ξ is obtained so that the undesired direction of the propagated wave SdiagΦrepEξ is suppressed frequency-wise. This is schematically illustrated in Fig. 3. Direction-wise signals SdiagΦrepEξ are obtained at the observation points, and its Fourier transform gives the frequency- and direction-wise energy |fSdiagΦrepEξ|2 of the propagated sound. By considering the unwanted direction for each frequency and suppressing it, propagated sound with desired directivity, which is realized by the estimated initial value, is obtained. For suppressing the undesired component, the weighted sum Wq22=m|wmqm|2 is minimized so that the component multiplied with a larger weight is suppressed more than that with a smaller weight. Therefore, a solution to the above optimization problem provides the initial value corresponding to a directive sound source whose directivity is defined by the weights W.

FIG. 3.

(Color online) Two examples of the direction-wise energy (non-directive and directive sound sources) in the frequency domain. Top row: instances of the sound fields p[n]. Middle row: direction-wise signals at the (green) observation points SdiagΦrepEξ. Bottom row: square root of the frequency- and direction-wise energy |fSdiagΦrepEξ| depicted as the right figure in Fig. 2.

FIG. 3.

(Color online) Two examples of the direction-wise energy (non-directive and directive sound sources) in the frequency domain. Top row: instances of the sound fields p[n]. Middle row: direction-wise signals at the (green) observation points SdiagΦrepEξ. Bottom row: square root of the frequency- and direction-wise energy |fSdiagΦrepEξ| depicted as the right figure in Fig. 2.

Close modal

This well-known problem is equivalent to the following minimization problem of the Rayleigh quotient:

minξξTMξξTξ,
(18)

where M is a real positive semidefinite matrix,

M=(WfSdiagΦrepE)*(WfSdiagΦrepE)
(19)
=ETΦrepTSdiagTf*WTWfSdiagΦrepE,
(20)

and * denotes complex conjugate transpose (note that W must be symmetric with respect to the Nyquist frequency so that f*WTWf=f*W2f becomes a real matrix, where f* corresponds to the inverse Fourier transform). A solution to this problem is the eigenvector corresponding to the minimum eigenvalue of M. Since solving the eigenvalue problem of M is effortless if the size of the matrix, decided by the length of the initial value vector ξ, is moderate, the proposed method can easily obtain a directive sound sources after M is calculated.

However, there are two issues which prevent direct implementation of the proposed method mentioned above. First, direct minimization of Eq. (18) might end up with a solution which contains a very high frequency around the spatial Nyquist frequency that may result in some unrealistic sound field after the propagation. Such a situation occurs when the desired directivity is not allowed physically and/or difficult to be approximated by FDTD. To circumvent this issue, use of a spatially smooth basis is additionally proposed in Sec. III C. Another issue is the calculation of the matrix M because the size of each matrix in Eq. (20) is exceedingly huge to store in the memory. This issue is solved by matrix-free implementation, as mentioned in the previous section, which is introduced in Sec. III D.

For controlling the spatial frequency of the estimated initial value, basis expansion is utilized in the proposed method. As the initial value is defined on the circular region, any function defined inside the circle may be used for the expansion. In this paper, an orthonormal basis on the region is chosen for that, in particular, the Dirichlet eigenfunctions of the Laplace operator are considered:

b(r)=αJ0(λ0r),
(21)
bmcos(r,θ)=αmcosJm(λmr)cos(mθ),
(22)
bmsin(r,θ)=αmsinJm(λmr)sin(mθ),
(23)

where Jm is the mth-order Bessel function, α is a normalization factor, m ≥ 1,  ≥ 1, and λm is a coefficient so that λmrinit becomes the th positive root of the Bessel function Jm(λmrinit) = 0 (radius of the region of the initial value is written as rinit). These functions are smooth and zero at the boundary of the region of the initial value. Their spatial frequencies increase as the parameters m and increase, and therefore, by expanding the initial value on them, the maximum spatial frequency of the initial value is decided by the maximum values of m and .

Let the number of basis elements {b,bmcos,bmsin}, decided by the choice of (m, ), be written as NB, and the expansion coefficients corresponding to pin[0] and pin[1] [see Eq. (10)] be denoted by η[0] and η[1], respectively. By considering an Nin × NB matrix BNin×NB consisting of the NB basis elements sampled at the Nin grid points, the initial value ξ can be expanded by the basis as

ξ=Bdiagη(η=[η[0]T,η[1]T]T),
(24)

where Bdiag=blkdiag(B,B)2Nin×2NB. Based on these notations, the proposed method is rewritten as a problem of finding the optimal coefficient η by solving the following optimization problem:

minηWfSdiagΦrepEBdiagη22s.t.η22=1.
(25)

This problem can also be solved through the eigenvalue problem of the corresponding matrix,

(WfSdiagΦrepEBdiag)*(WfSdiagΦrepEBdiag)=BdiagTETΦrepTSdiagTf*WTWfSdiagΦrepEBdiag,
(26)

where the eigenvector corresponding to the minimum eigenvalue of this matrix is a solution to Eq. (25). After obtaining the solution, the initial value to be propagated is recovered via Eq. (24). Note that this basis expansion can reduce the computational cost when the number of the basis elements is sufficiently smaller than that of the grid points of the initial value (NB < Nin).

Since the optimization problems in Eqs. (17) and (25) can be solved through the eigenvalue problems, implementation of the proposed method has no difficulty after the construction of the matrices in Eqs. (20) and (26) because eigenvalue solvers are readily available in most of the programming languages. However, construction of all matrices in Eqs. (20) and (26) is almost impossible owing to the gigantic size of the matrix Φrep which might contain more than a hundred trillion elements for a typical two-dimensional field (and this size explodes in the three-dimensional case). To overcome such troublesome memory requirement, the so-called matrix-free optimization method36–38 is utilized.

The basic idea of the matrix-free method is to implement a linear operator as a function (in terms of programming language) instead of a matrix. For instance, multiplying a diagonal matrix W can be written by the element-wise product which eliminates memory usage for zeros in the off-diagonal elements of W. Likewise, multiplying f can be implemented by the fast Fourier transform algorithm. Sdiag and E can also easily be written as the functions, where Sdiag picks up the nNobs elements from a 2nNp-dimensional vector whose indices coincide with those corresponding to the observation points, and E substitutes a 2Nin vector into the 2Np-dimensional vector consisting of all zeros (see Fig. 12 in  Appendix B). In addition, Φrep can be implemented through the usual FDTD program because it consists of Φ which represents the one step of an FDTD scheme. One may write Bdiag as a function, too. Yet, storing Bdiag should be more beneficial than the matrix-free implementation because B is a dense and reasonably small matrix in typical situations.

In addition to the above matrices, their transposed versions are also required to be implemented as functions. WT and f* are simple because WT= W and f* is the inverse Fourier transform. Transpositions of the other matrices might require time for derivation, and thus they are detailed in  Appendix B.

After all matrices and their transpositions are written as functions, an iterative algorithm can be utilized for solving the eigenvalue problems. An example of iterative solvers is the Lanczos algorithm which only requires calculation of the matrix-vector product. Since the matrix-vector product can be implemented in the matrix-free manner as explained above, the eigenvalue problems can be solved without explicitly constructing the matrices.

The importance of the matrix-free implementation is illustrated by examples of memory usage calculated by the experimental condition similar to that in Sec. IV (summarized in Table I). The memory requirements of the proposed method for two- and three-dimensional sound fields are shown in Table II, where the total memory usage and details of each variable were obtained by actually implementing the proposed method in double precision. From the bottom rows, it is obvious that explicitly storing the matrices is exceedingly demanding for the two-dimensional case and impossible for the three-dimensional case [they are O(hX4) and O(hX6), respectively]. In contrast, the matrix-free method can be implemented with the realistic memory usage, reduced to O(hX2) (2D) and O(hX3) (3D), dominated by the grid of FDTD. Similarly, almost all computations are required by the FDTD calculation as the other operations are relatively cheap. Therefore, the proposed method can be performed when the chosen FDTD scheme is computable.

TABLE I.

Simulation condition.

Size of sound field8 m × 8 m
Shape of region for initial value ξ Circular disk 
Placement of observation points pobs Circle 
Sound speed c 340 m/s 
Width of Cartesian grid hX 0.02 m 
Size of time step hT 1/30 000 s 
Total number of time steps n 360 steps 
Length of Fourier transform NDFT 512 samples 
Size of sound field8 m × 8 m
Shape of region for initial value ξ Circular disk 
Placement of observation points pobs Circle 
Sound speed c 340 m/s 
Width of Cartesian grid hX 0.02 m 
Size of time step hT 1/30 000 s 
Total number of time steps n 360 steps 
Length of Fourier transform NDFT 512 samples 
TABLE II.

Memory requirement for the proposed method.

2D3D
Width of Cartesian grid hX 0.02 0.01 0.02 0.01 
No. elements in ξ 309 1249 4145 33377 
Memory for ξ(KB) 2.47 9.99 33.16 267.02 
No. elements in fPobs 43 008 71 680 1 348 608 3 820 032 
Memory for fPobs(MB) 0.34 0.57 10.79 30.56 
No. grid points of FDTD 4012 8012 4013 8013 
Memory for FDTD (MB) 2.57 10.27 1031.70 8222.76 
Total memory (direct) 68 TB 1186 TB 10 EB 760 EB 
Total memory (matrix-free) 3.02 MB 11.26 MB 1.04 GB 8.26 GB 
2D3D
Width of Cartesian grid hX 0.02 0.01 0.02 0.01 
No. elements in ξ 309 1249 4145 33377 
Memory for ξ(KB) 2.47 9.99 33.16 267.02 
No. elements in fPobs 43 008 71 680 1 348 608 3 820 032 
Memory for fPobs(MB) 0.34 0.57 10.79 30.56 
No. grid points of FDTD 4012 8012 4013 8013 
Memory for FDTD (MB) 2.57 10.27 1031.70 8222.76 
Total memory (direct) 68 TB 1186 TB 10 EB 760 EB 
Total memory (matrix-free) 3.02 MB 11.26 MB 1.04 GB 8.26 GB 

In order to illustrate the properties of the proposed method, numerical simulations were performed. The condition common for all simulations is summarized in Table I, where the second-order scheme explained in Sec. II B was utilized for the FDTD updates. The spatial grid and the time step were chosen so that the FDTD scheme was stable for every simulation performed here.

In the experiments, the diagrams of frequency- and direction-wise magnitude, as in Figs. 2(b) and 3, are utilized for summarizing the results. Therefore, before proceeding to the results, some examples of the diagrams and the corresponding wave propagation are shown in Fig. 4. Since the desired directivity in the leftmost column is represented by 1 (desired) and 0 (undesired), the weighting matrix W was set in the reversed manner, 0 (desired) and 1 (undesired). The radius of the region for the initial value rinit was set to 0.2 m, the radius of the observation points robs was 0.3 m, and the number of the basis elements NB was 50. As in the figures, the desired directivities were approximately realized by the proposed method. The reason that the results were not exact but approximation is because the possible pattern of the directivity is naturally limited by the physical constraint (imposed by the wave equation) and the FDTD scheme. Note that, although the frequency responses of the obtained results seem unbalance as the low-frequency components are smaller than the other band, it is not difficult to compensate that effect by using a post filter because the low-frequency components are not vanished.

FIG. 4.

(Color online) Examples of the directivity diagrams and the corresponding propagated waves. Simulation conditions were rinit = 0.2 m, robs = 0.3 m, and NB = 50. The propagated waves are normalized to ±1 by their maximum absolute value, while the colorbar in the range of ±0.1 is utilized for clearly illustrating the propagated pulses.

FIG. 4.

(Color online) Examples of the directivity diagrams and the corresponding propagated waves. Simulation conditions were rinit = 0.2 m, robs = 0.3 m, and NB = 50. The propagated waves are normalized to ±1 by their maximum absolute value, while the colorbar in the range of ±0.1 is utilized for clearly illustrating the propagated pulses.

Close modal

For investigating the properties of the proposed method, the relation between its parameters (rinit, robs, and NB) and the obtained results are shown in this section. To quantitatively assess the approximation quality, the relative energy ϒ is defined as follows:

ϒ=χdesfpeval22/fpeval22,
(27)

where χdes{0,1}NDFTNobs×NDFTNobs is a diagonal matrix whose diagonal entry is 1 if it corresponds to the desired part (brighter yellow box in the desired directivity diagram) and otherwise 0 (undesired part). For setting the same conditions for all experiments, the directivity diagram was calculated at the evaluation points, as in Fig. 5, whose distance to the origin reval was 1.5 m. This ϒ ∈ [0,1] becomes 1 if all energy is concentrated inside the desired part and 0 if all energy is outside it.

FIG. 5.

(Color online) Setting of evaluation points and the corresponding directivity diagram. (a) The evaluation points on the dashed (red) circle were set constant (1.5 m) for all experiments, while the other geometrical settings (rinit and robs) may vary for each experiment. (b) All directivity diagrams in Sec. IV B, depicted as Fig. 2, were calculated at the evaluation points instead of the observation points.

FIG. 5.

(Color online) Setting of evaluation points and the corresponding directivity diagram. (a) The evaluation points on the dashed (red) circle were set constant (1.5 m) for all experiments, while the other geometrical settings (rinit and robs) may vary for each experiment. (b) All directivity diagrams in Sec. IV B, depicted as Fig. 2, were calculated at the evaluation points instead of the observation points.

Close modal

At first, the relation between the radius of the initial value rinit and the obtained result was investigated. In this experiment, the maximum spatial frequency of the basis elements was set to the (almost) same value for every radius rinit. To do so, the basis elements were chosen automatically for each rinit according to their spatial frequencies measured by the spatial Fourier transform because the spatial frequency of a basis element depends on its radius. The chosen basis elements for each rinit are indicated by the circles in Fig. 6. rinit was changed from 0.05 to 0.5 m, while robs was set to (rinit + 0.1) m. The directivity diagrams of the obtained results are shown in Fig. 7, and their relative energies ϒ are summarized in Fig. 8(i). These figures indicate that the radius of the initial value is important for representing a complex directivity.

FIG. 6.

Illustration of the parameter (m, ) of the chosen basis elements for each rinit corresponding to Figs. 7 and 8.

FIG. 6.

Illustration of the parameter (m, ) of the chosen basis elements for each rinit corresponding to Figs. 7 and 8.

Close modal
FIG. 7.

(Color online) Directivity diagrams obtained by the proposed method (the radius of the initial value rinit varied).

FIG. 7.

(Color online) Directivity diagrams obtained by the proposed method (the radius of the initial value rinit varied).

Close modal
FIG. 8.

(Color online) (i) Relation between approximation accuracy and the radius of the initial value rinit. (ii) Relation between approximation accuracy and the radius of the observation points robs (rinit = 0.2, NB = 50). (iii) Relation between approximation accuracy and the number of the basis elements NB, (rinit = 0.2, robs = 0.3). The alphabets (a)–(c) correspond to each desired directivity shown in Figs. 7 and 9.

FIG. 8.

(Color online) (i) Relation between approximation accuracy and the radius of the initial value rinit. (ii) Relation between approximation accuracy and the radius of the observation points robs (rinit = 0.2, NB = 50). (iii) Relation between approximation accuracy and the number of the basis elements NB, (rinit = 0.2, robs = 0.3). The alphabets (a)–(c) correspond to each desired directivity shown in Figs. 7 and 9.

Close modal

Then, the relation between the radius of the observation points robs and the obtained result was investigated by fixing the other parameters, rinit(= 0.2 m) and NB (= 50). The results are shown in Fig. 8(ii), which indicates that robs has little effect on the approximation accuracy. Note that a larger radius robs may increase the number of observation points Nobs, which results in higher computational cost. Therefore, a smaller radius may be preferable for the calculation.

Finally, the relation between the number of basis elements NB and the obtained result was investigated by fixing rinit(= 0.2 m) and robs (= 0.3 m). The directivity diagrams are shown in Fig. 9, where the utilized basis elements are illustrated in Fig. 10, and its summary is shown in Fig. 8(iii). This result indicates that the basis elements of higher order have less impact on the accuracy comparing to the lower-order elements. As the possible directional pattern is limited by the wave equation, increasing the number of basis elements has little effect after it reaches to a sufficient number. Note that, according to Fig. 8, case (a) resulted in notably worse scores than (b) and (c) for all three experiments (i)–(iii), which indicates that a sharp directivity cannot be achieved by only using low-frequency components.

FIG. 9.

(Color online) Directivity diagrams obtained by the proposed method (the number of the basis elements NB varied).

FIG. 9.

(Color online) Directivity diagrams obtained by the proposed method (the number of the basis elements NB varied).

Close modal
FIG. 10.

Illustration of the parameter (m, ) of the chosen basis elements for each NB corresponding to Figs. 8 and 9.

FIG. 10.

Illustration of the parameter (m, ) of the chosen basis elements for each NB corresponding to Figs. 8 and 9.

Close modal

In summary, the radius of the initial value rinit and the number of the basis elements NB are important parameters for the proposed method, while the radius of the observation points robs has little impact. Therefore, both rinit and NB should be set sufficiently large for a given directivity, while the sufficient number of these values seems to depend on the directivity.

In this section, the proposed method is compared with the conventional method which realizes a directive sound source in FDTD.32 The conventional method approximates the directivity directly using the spherical harmonics (Fourier series) up to 8th order. To set the same number of basis elements for both the conventional and proposed methods, NB = 17 was chosen here, where the basis elements of the proposed method were automatically chosen according to their spatial frequencies as in Sec. IV B. Since the conventional method does not consider the frequency-wise directivity, the weight of the proposed method was set to equally contain all frequencies for comparison. The radius of the initial condition rinit was set to 0.5 m, while the standard deviation of the Gaussian function of the conventional method was also set to 0.5 m.

The obtained results are shown in Fig. 11. The top and middle row indicate the results using the same number of basis elements for both the conventional and proposed method. Therefore, the top two rows give a fair comparison in terms of the number of basis elements. However, in such a case, the basis of the conventional method contains higher angular frequency components than that of the proposed method because the conventional method only considers angular basis while the proposed method considers both radial and angular components [r and θ in Eqs. (21)–(23)]. To compare the results using bases containing the same angular frequency, the bottom row shows the result obtained by the proposed method for NB = 29. That is, comparison between the top and bottom rows is fair in terms of the angular frequency components.

FIG. 11.

(Color online) Comparison of the results obtained by the conventional and proposed methods. The top and middle rows represent the results using the same number of basis elements, while the top and bottom rows are that using the basis containing the same angular frequency. The propagated waves were normalized to ±1 by their maximum absolute value, while the colorbar in the range of ±0.2 is utilized for clearly illustrating the propagated pulses.

FIG. 11.

(Color online) Comparison of the results obtained by the conventional and proposed methods. The top and middle rows represent the results using the same number of basis elements, while the top and bottom rows are that using the basis containing the same angular frequency. The propagated waves were normalized to ±1 by their maximum absolute value, while the colorbar in the range of ±0.2 is utilized for clearly illustrating the propagated pulses.

Close modal

From the figures, difference of the characteristics of each method can be seen. First, from the directivity, it can be seen that the conventional method obtained a good response at the peak directions (0 and 180 deg), while the other direction corresponding to small energy contains larger error. This should be because the conventional method tries to approximate the given directivity directly. In contrast, the proposed method tries to suppress the direction having small energy that resulted in the smaller error for such directions (around 90 and 270 deg) and a larger error for other directions comparing to the conventional method. Second, since the conventional method considers not only the initial value but also the forcing term for t > 0, the propagated wave contains components around the center for all time instances. This property is common for all methods considering a forcing term (and also for all frequency-domain methods) because the source must emit sound at t > 0. In contrast, as the proposed method only considers the initial value at t = 0 and no forcing term for t > 0, its result does not contain the components around the center. In addition, the proposed method obtained the better result when the number of the angular components was set to the same number as the conventional method.

In this paper, a method for imposing directivity into FDTD was proposed. It was formulated as an optimization problem of minimizing the energy of propagated waves at the undesired directions. By reducing the optimization problem into the eigenvalue problem, the optimized result can be obtained by a usual eigenvalue solver. The feature of the proposed method distinct from the other methods is the explicit consideration of the update equation of an FDTD scheme in the optimization, which makes it possible to obtain a result aware of the effect of discretization contained in the FDTD scheme.

While this paper showed the potentiality of the fully optimization-based representation of directive sound sources in FDTD, there is still plenty of room for improving the method toward the optimal directive-source representation. For example, it should be desired to impose a time and spatially varying forcing term because the main limitation of the current method is consideration of only the initial value, which makes treatment of flexible sources, such as moving sources, impossible. Extension of the proposed method by considering multiple time steps in the optimization should be able to allow such forcing term. In addition, a method for choosing and/or optimizing the basis might be important to consider for making the method more efficient. Moreover, as the directivity is controlled only by the weight W in the proposed method, it is desirable to have more flexibility on the control of the directivity. These considerations are remained as the future works.

As mentioned at the beginning of Sec. III, any FDTD scheme other than that in Sec. II B can be utilized in the proposed method. Here, an example of a leapfrog-type FDTD scheme is briefly shown for demonstrating that.

A leapfrog-type FDTD scheme, which aims to approximately solve the set of force equations and continuity equation, may be written as follows:

u[n+1]=u[n]+(hT/ρ)DxTp[n],
(A1)
v[n+1]=v[n]+(hT/ρ)DyTp[n],
(A2)
p[n+1]=p[n]κhT(Dxu[n+1]+Dyv[n+1]),
(A3)

where ρ is density, κ is bulk modulus, u and v are the particle velocities in x and y directions, respectively, and Dx and Dy are matrices corresponding to the finite difference approximations of the first-order spatial derivatives.

This FDTD scheme can also be written in the matrix-vector form in Eq. (8). Let the vector ζ[n] be defined as

ζ[n]=[p[n]T,u[n]T,v[n]T]T.
(A4)

Then, the matrix Φ, defined as

Φ=(IDκhTDxκhTDy(hT/ρ)DxTIO(hT/ρ)DyTOI),
(A5)

provides the update ζ[n+1]=Φζ[n] as in Eq. (8), where

D=(κhT2/ρ)(DxDxT+DyDyT).
(A6)

As mentioned in Sec. III D, all matrices in the proposed method can be implemented as functions so that the memory requirement for storing the matrices is avoided. Here, detailed procedures for writing the matrix-free program are shown in order to aid understanding the concept. All procedures introduced in this appendix can be derived by writing down the matrices and transpose them.

For implementing the proposed method, one has to calculate the multiplication of SdiagΦrepEBdiag, and then apply the combination of the Fourier transform and its inverse with the weighting. After them, multiplication of BdiagTETΦrepTSdiagT must be calculated as in Eq. (26). Fortunately, these matrices do not have to be stored in memory because they can be computed by a set of procedures. Let E and S be the index sets corresponding to the embedding and selection matrices E and S, respectively, as in Fig. 12. Then, the matrix multiplication can be implemented as the functions shown in Algorithms 1 and 2, where the subroutines in Algorithms 3 and 4 are utilized within them. Although these procedures are specific to the FDTD scheme introduced in Sec. II B as an example, any other FDTD scheme as in  Appendix A can also be implemented in a similar manner.

Algorithm 1:

Calculate pobs=SdiagΦrepEBdiagη

Input:η,B,E,S 
Output:pobs 
1: Separate η into η[1] and η[0] 
2: pin[1]=Bη[1], pin[0]=Bη[0] 
3: p[1]=zeropad(pin[1],E), p[0]=zeropad(pin[0],E) 
4: fork = 1, 2,…, n − 1 do 
5:    p[k+1]=FDTD(p[k],p[k1]) 
6:    pobs[k+1]=select(p[k+1],S) 
7: end for 
8: pobs=[pobs[1]T,pobs[2]T,,pobs[n]T]T 
Input:η,B,E,S 
Output:pobs 
1: Separate η into η[1] and η[0] 
2: pin[1]=Bη[1], pin[0]=Bη[0] 
3: p[1]=zeropad(pin[1],E), p[0]=zeropad(pin[0],E) 
4: fork = 1, 2,…, n − 1 do 
5:    p[k+1]=FDTD(p[k],p[k1]) 
6:    pobs[k+1]=select(p[k+1],S) 
7: end for 
8: pobs=[pobs[1]T,pobs[2]T,,pobs[n]T]T 
Algorithm 2:

Calculate η̂=BdiagTETΦrepTSdiagTpobs

Input:pobs,E,S 
Output:η̂ 
1: Separate pobs into pobs[1],pobs[2],,pobs[n] 
2: Set p̂[1]=0Np, p̂[0]=0Np 
3: fork = 1, 2,…, n − 1 do 
4:    p̂[1]=p̂[1]+zeropad(pobs[nk+1],S) 
5:    q=adjointFDTD(p̂[1],p̂[0]) 
6:    p̂[0]=p̂[1], p̂[1]=q 
7: end for 
8: p̂in[1]=select(p̂[1],E), p̂in[0]=select(p̂[0],E) 
9: η̂[1]=BTp̂in[1], η̂[0]=BTp̂in[0] 
10: η̂=[η̂[1]T,η̂[0]T]T 
Input:pobs,E,S 
Output:η̂ 
1: Separate pobs into pobs[1],pobs[2],,pobs[n] 
2: Set p̂[1]=0Np, p̂[0]=0Np 
3: fork = 1, 2,…, n − 1 do 
4:    p̂[1]=p̂[1]+zeropad(pobs[nk+1],S) 
5:    q=adjointFDTD(p̂[1],p̂[0]) 
6:    p̂[0]=p̂[1], p̂[1]=q 
7: end for 
8: p̂in[1]=select(p̂[1],E), p̂in[0]=select(p̂[0],E) 
9: η̂[1]=BTp̂in[1], η̂[0]=BTp̂in[0] 
10: η̂=[η̂[1]T,η̂[0]T]T 
Algorithm 3:

Calculate p[n+1]=FDTD(p[n],p[n1])

Input:p[n],p[n1] 
Output:p[n+1] 
1: for allk, ldo 
2: pk,l[n+1]=c2(hT/hX)2(pk+1,l[n]2pk,l[n]+pk1,l[n]+pk,l+1[n]2pk,l[n]+pk,l1[n])+2pk,l[n]pk,l[n1] 
3: end for 
Input:p[n],p[n1] 
Output:p[n+1] 
1: for allk, ldo 
2: pk,l[n+1]=c2(hT/hX)2(pk+1,l[n]2pk,l[n]+pk1,l[n]+pk,l+1[n]2pk,l[n]+pk,l1[n])+2pk,l[n]pk,l[n1] 
3: end for 
Algorithm 4:

Calculate q=adjointFDTD(p̂[1],p̂[0])

Input:p̂[1],p̂[0] 
Output:q 
1: for allk, ldo 
2: qk,l=c2(hT/hX)2(p̂k+1,l[1]2p̂k,l[1]+p̂k1,l[1]+p̂k,l+1[1]2p̂k,l[1]+p̂k,l1[1])+2p̂k,l[1]+p̂k,l[0] 
3: end for 
Input:p̂[1],p̂[0] 
Output:q 
1: for allk, ldo 
2: qk,l=c2(hT/hX)2(p̂k+1,l[1]2p̂k,l[1]+p̂k1,l[1]+p̂k,l+1[1]2p̂k,l[1]+p̂k,l1[1])+2p̂k,l[1]+p̂k,l[0] 
3: end for 
FIG. 12.

Schematic of the embedding matrix E, the selection matrix S, and their transpositions.

FIG. 12.

Schematic of the embedding matrix E, the selection matrix S, and their transpositions.

Close modal
1.
A.
Krokstad
,
S.
Strom
, and
S.
Sørsdal
, “
Calculating the acoustical room response by the use of a ray tracing technique
,”
J. Sound Vib.
8
(
1
),
118
125
(
1968
).
2.
J. B.
Allen
and
D. A.
Berkley
, “
Image method for efficiently simulating small-room acoustics
,”
J. Acoust. Soc. Am.
65
(
4
),
943
950
(
1979
).
3.
D.
Botteldooren
, “
Finite-difference time-domain simulation of low-frequency room acoustic problems
,”
J. Acoust. Soc. Am.
98
(
6
),
3302
3308
(
1995
).
4.
D.
Murphy
,
A.
Kelloniemi
,
J.
Mullen
, and
S.
Shelley
, “
Acoustic modeling using the digital waveguide mesh
,”
IEEE Signal Process. Mag.
24
(
2
),
55
66
(
2007
).
5.
T.
Shuku
and
K.
Ishihara
, “
The analysis of the acoustic field in irregularly shaped rooms by the finite element method
,”
J. Sound Vib.
29
(
1
),
67
76
(
1973
).
6.
T.
Sakuma
and
Y.
Yasuda
, “
Fast multipole boundary element method for large-scale steady-state sound field analysis. Part I: Setup and validation
,”
Acta Acust. Acust.
88
(
4
),
513
525
(
2002
).
7.
M.
Hornikx
,
R.
Waxler
, and
J.
Forssén
, “
The extended Fourier pseudospectral time-domain method for atmospheric sound propagation
,”
J. Acoust. Soc. Am.
128
(
4
),
1632
1646
(
2010
).
8.
S.
Bilbao
,
Wave and Scattering Methods for Numerical Simulation
(
Wiley
,
London
,
2004
).
9.
K.
Kowalczyk
and
M.
van Walstijn
, “
Wideband and isotropic room acoustics simulation using 2-D interpolated FDTD schemes
,”
IEEE Trans. Audio, Speech, Lang. Process.
18
(
1
),
78
89
(
2010
).
10.
K.
Kowalczyk
and
M.
van Walstijn
, “
A comparison of nonstaggered compact FDTD schemes for the 3D wave equation
,” in
IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)
, test (March
2010
), pp.
197
200
.
11.
K.
Kowalczyk
and
M.
van Walstijn
, “
Room acoustics simulation using 3-D compact explicit FDTD schemes
,”
IEEE Trans. Audio, Speech, Lang. Process.
19
(
1
),
34
46
(
2011
).
12.
S.
Bilbao
, “
Optimized FDTD schemes for 3-D acoustic wave propagation
,”
IEEE Trans. Audio, Speech, Lang. Process.
20
(
5
),
1658
1663
(
2012
).
13.
S.
Bilbao
and
B.
Hamilton
, “
Construction and optimization techniques for high order schemes for the two-dimensional wave equation
,” in
Proc. Int. Congr. Acoust. (ICA)
(
2013
), Vol.
19
.
14.
A.
Southern
,
S.
Siltanen
,
D. T.
Murphy
, and
L.
Savioja
, “
Room impulse response synthesis and validation using a hybrid acoustic model
,”
IEEE Trans. Audio, Speech, Lang. Process.
21
(
9
),
1940
1952
(
2013
).
15.
J. V.
Mourik
and
D.
Murphy
, “
Explicit higher-order FDTD schemes for 3D room acoustic simulation
,”
IEEE/ACM Trans. Audio, Speech, Lang. Process.
22
(
12
),
2003
2011
(
2014
).
16.
J.
Botts
and
L.
Savioja
, “
Spectral and pseudospectral properties of finite difference models used in audio and room acoustics
,”
IEEE Trans. Audio, Speech, Lang. Process.
22
(
9
),
1403
1412
(
2014
).
17.
J.
Botts
and
L.
Savioja
, “
Effects of sources on time-domain finite difference models
,”
J. Acoust. Soc. Am.
136
(
1
),
242
247
(
2014
).
18.
J.
Sheaffer
,
M.
van Walstijn
, and
B.
Fazenda
, “
Physical and numerical constraints in source modeling for finite difference simulation of room acoustics
,”
J. Acoust. Soc. Am.
135
(
1
),
251
261
(
2014
).
19.
C.
Spa
,
A.
Rey
, and
E.
Hernández
, “
A GPU implementation of an explicit compact FDTD algorithm with a digital impedance filter for room acoustics applications
,”
IEEE/ACM Trans. Audio, Speech, Lang. Process.
23
(
8
),
1368
1380
(
2015
).
20.
J.
Saarelma
,
J.
Botts
,
B.
Hamilton
, and
L.
Savioja
, “
Audibility of dispersion error in room acoustic finite-difference time-domain simulation as a function of simulation distance
,”
J. Acoust. Soc. Am.
139
(
4
),
1822
1832
(
2016
).
21.
W.
Huang
and
L.
Yigang
, “
A Higher-order finite-difference time-domain method for sound propagation
,” in
Int. Conf. Audio, Lang. Image Process. (ICALIP)
(July
2016
), pp.
136
140
.
22.
B.
Hamilton
and
S.
Bilbao
, “
FDTD Methods for 3-D Room Acoustics Simulation with High-order accuracy in space and time
,”
IEEE/ACM Trans. Audio, Speech, Lang. Process.
25
(
11
),
2112
2124
(
2017
).
23.
F.
Otondo
and
J. H.
Rindel
, “
The influence of the directivity of musical instruments in a room
,”
Acta Acust. Acust.
90
(
6
),
1178
1184
(
2004
).
24.
P.
Zhu
,
F.
Mo
,
J.
Kang
, and
G.
Zhu
, “
Comparisons between simulated and in-situ measured speech intelligibility based on (binaural) room impulse responses
,”
Appl. Acoust.
97
,
65
77
(
2015
).
25.
A.
Snakowska
,
J.
Jurkiewicz
, and
L.
Gorazd
, “
A hybrid method for determination of the acoustic impedance of an unflanged cylindrical duct for multimode wave
,”
J. Sound Vib.
396
,
325
339
(
2017
).
26.
R. H. C.
Wenmaekers
,
C. C. J. M.
Hak
,
M. C. J.
Hornikx
, and
A. G.
Kohlrausch
, “
Sensitivity of stage acoustic parameters to source and receiver directivity: Measurements on three stages and in two orchestra pits
,”
Appl. Acoust.
123
,
20
28
(
2017
).
27.
H.
Hacihabiboǧlu
,
B.
Günel
, and
A. M.
Kondoz
, “
Time-domain simulation of directive sources in 3-D digital waveguide mesh-based acoustical models
,”
IEEE Trans. Audio, Speech, Lang. Process.
16
(
5
),
934
946
(
2008
).
28.
Y.
Tamura
,
K.
Yatabe
, and
Y.
Oikawa
, “
Least-squares estimation of sound source directivity using convex selector of a better solution
,”
Acoust. Sci. Technol.
38
(
3
),
128
136
(
2017
).
29.
F.
Georgiou
and
M.
Hornikx
, “
Incorporating directivity in the Fourier pseudospectral time-domain method using spherical harmonics
,”
J. Acoust. Soc. Am.
140
(
2
),
855
865
(
2016
).
30.
J.
Escolano
,
J. J.
López
, and
B.
Pueo
, “
Directive sources in acoustic discrete-time domain simulations based on directivity diagrams
,”
J. Acoust. Soc. Am.
121
(
6
),
EL256
EL262
(
2007
).
31.
D. T.
Murphy
,
A.
Southern
, and
L.
Savioja
, “
Source excitation strategies for obtaining impulse responses in finite difference time domain room acoustics simulation
,”
Appl. Acoust.
82
,
6
14
(
2014
).
32.
S.
Sakamoto
and
R.
Takahashi
, “
Directional sound source modeling by using spherical harmonic functions for finite-difference time-domain analysis
,” in
Proceedings of Meetings on Acoustics ICA2013
, ASA (
2013
), Vol.
19
, p.
015129
.
33.
S.
Bilbao
and
B.
Hamilton
, “
Directional source modeling in wave-based room acoustics simulation
,” in
Workshop Appl. Signal Process. Audio Acoust. (WASPAA)
, IEEE (
2017
), pp.
121
125
.
34.
T.
Tachikawa
,
K.
Yatabe
, and
Y.
Oikawa
, “
3D sound source localization based on coherence-adjusted monopole dictionary and modified convex clustering
,”
Appl. Acoust.
139
,
267
281
(
2018
).
35.
D.
Takeuchi
,
K.
Yatabe
, and
Y.
Oikawa
, “
Realizing directional sound source in FDTD method by estimating initial value
,” in
IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP)
(
2018
), pp.
461
465
.
36.
S.
Diamond
and
S.
Boyd
, “
Matrix-free convex optimization modeling
,” in
Optimization and its Applications in Control and Data Sciences
(
Springer
,
New York
,
2016
), pp.
221
264
.
37.
J.
Folberth
and
S.
Becker
, “
Efficient adjoint computation for wavelet and convolution operators [lecture notes]
,”
IEEE Signal Process. Mag.
33
(
6
),
135
147
(
2016
).
38.
N.
Antonello
,
E. D.
Sena
,
M.
Moonen
,
P. A.
Naylor
, and
T.
van Waterschoot
, “
Room impulse response interpolation using a sparse spatio-temporal representation of the sound field
,”
IEEE/ACM Trans. Audio, Speech, Lang. Process.
25
(
10
),
1929
1941
(
2017
).