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.

## I. INTRODUCTION

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 tracing^{1} 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 sources^{28,30,31} and expansion into spherical harmonics or multipoles.^{29,32,33} Superposition of point sound sources is a popular methodology in sound field processing^{34} 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 method^{30} (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.

## II. FINITE-DIFFERENCE TIME-DOMAIN METHOD

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.

### A. Acoustic simulation based on FDTD

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

where *p* is sound pressure, *c* is sound speed, $r=(x,y)$ is a position vector, *t* is time, and $\u25b3$ 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).

### B. Example of FDTD scheme for wave equation

The sound pressure can be represented by its values on the Cartesian grid as $pi,j[n]\u2261p(hX\u2009i,hX\u2009j,hT\u2009n)$, where *h _{X}* and

*h*are the grid sizes of the spatial and time variables, respectively, $i,j\u2208\mathbb{N}$ are the spatial indices, and $n\u2208\mathbb{N}$ is the time index. FDTD approximates the derivatives in Eq. (1) by the finite-difference operators,

_{T}where *D _{tt}* and $\u25b3d$ 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

*D*and $\u25b3d$ can be utilized, a simple example of the time difference is considered here:

_{tt}Then, an update equation of FDTD can be written as

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 $\u25b3d$, one simple example is

Note that, for calculating the spatial difference $\u25b3d\u2009pi,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.

### C. Matrix representation of FDTD scheme

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 *N _{p}* be the total number of grid points of a sound field, and $p[n]\u2208\mathbb{R}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

where $\u25b3mtx\u2208\mathbb{R}Np\xd7Np$ is the matrix representing $\u25b3d$, i.e., the elements of $\u25b3mtx\u2009p[n]$ are $\u25b3d\u2009pi,j[n]$ which is calculated as in Eq. (5). By defining a state vector

Eq. (6) can be simplified to

where $\Phi \u2208\mathbb{R}2Np\xd72Np$ is a square block matrix,

$INp$ and $ONp$ are the identity and zero matrices, respectively, whose sizes are *N _{p}* ×

*N*. 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.

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

## III. PROPOSED METHOD

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 $\zeta [n]$ (see Appendix A for another example), and therefore the proposed method is not restricted to a specific FDTD scheme.

### A. Problem formulation

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 $\zeta [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,

whose length is much smaller than $\zeta [n]$ in Eq. (7) defined on the whole domain, i.e., *N*_{in}$\u226a$*N _{p}*. The initial value $\zeta [0]$ is recovered from $\xi $ by padding zero to the outside of the region, which is shortly written as

where $E\u2208{0,1}2Np\xd72Nin$ is an embedding matrix whose column has only single 1 and other entries are 0.

After propagating the initial value of *n* steps by FDTD as $\Phi nE\xi (=\Phi n\zeta [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 *r*_{obs}, as illustrated in Fig. 2(a). The sound pressure at the observation points in *n*th time step is represented by an *N*_{obs}-dimensional vector $pobs[n]\u2208\mathbb{R}Nobs$ whose length is also much smaller than $\zeta [n]\u2009(Nobs\u226aNp)$. The relation between these vectors is written similar to Eq. (11) as

where $S\u2208{0,1}Nobs\xd72Np$ 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 $\zeta [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 $\xi $ and the sound pressure at the observation points $pobs[n]$ is expressed as

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 *nN*_{obs}-dimensional vector $pobs$ as

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

where $Sdiag=blkdiag(S,S,\u2026,S)\u2208{0,1}nNobs\xd72nNp$ is the block diagonal matrix consisting of the selection matrices, $blkdiag(A1,\u2026,An)$ is the operation generating a block diagonal matrix from inputted matrices $A1,\u2026,An$, and $\Phi rep=[\u2009\Phi T,(\Phi 2)T,\u2026,(\Phi n)T\u2009]T\u2208\mathbb{R}2nNp\xd72Np$ is the matrix generating the sound fields of all time steps from the initial value $E\xi \u2009(=\zeta [0])$.

This equation indicates that the initial value inside the region $\xi $ 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=\u2211m=1n|pobs[m]|2$ is available as data. Therefore, the initial value $\xi $ must be estimated from a given direction-wise non-negative data, such as $pobsengy$, so that the propagated sound $Sdiag\Phi repE\xi $ approximates it.

### B. Proposed method for estimating optimal initial value

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 $f\u2208\u2102NDFTNobs\xd7nNobs$ be a matrix, constructed from *N*_{obs} discrete Fourier transform matrices (with zero padding) whose size is *N*_{DFT} × *n*, such that the product $fpobs$ acts as the Fourier transform at each observation point and results in a vector consisting of *N*_{obs} spectra whose length is *N*_{DFT} (here, *N*_{DFT} ≥ *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 $\xi $ is related to the direction-wise Fourier spectra at the observation points $fpobs$ as

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\Phi repE\xi $ results in the control of the frequency-wise directivity.

To control the direction- and frequency-wise energy, a weighting matrix $W\u2208\mathbb{R}NDFTNobs\xd7NDFTNobs$, 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:

where $\Vert q\Vert 22=\u2211m|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 $\Vert \xi \Vert 22=1$ is required to avoid the trivial solution $\xi =0$. By solving this problem, the initial value $\xi $ is obtained so that the undesired direction of the propagated wave $Sdiag\Phi repE\xi $ is suppressed frequency-wise. This is schematically illustrated in Fig. 3. Direction-wise signals $Sdiag\Phi repE\xi $ are obtained at the observation points, and its Fourier transform gives the frequency- and direction-wise energy $|fSdiag\Phi repE\xi |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 $\Vert Wq\Vert 22=\u2211m|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*.

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

where *M* is a real positive semidefinite matrix,

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 $\xi $, 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.

### C. Representing initial value by smooth orthonormal basis

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:

where *J _{m}* is the

*m*th-order Bessel function,

*α*is a normalization factor,

*m*≥ 1, $\u2113$ ≥ 1, and

*λ*

_{m}_{$\u2113$}is a coefficient so that

*λ*

_{m}_{$\u2113$}

*r*

_{init}becomes the $\u2113$ th positive root of the Bessel function

*J*(

_{m}*λ*

_{m}_{$\u2113$}

*r*

_{init}) = 0 (radius of the region of the initial value is written as

*r*

_{init}). 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 $\u2113$ 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 $\u2113$.

Let the number of basis elements ${b\u2113,bm\u2113cos,bm\u2113sin}$, decided by the choice of (*m*, $\u2113)$, be written as *N _{B}*, and the expansion coefficients corresponding to $pin[0]$ and $pin[1]$ [see Eq. (10)] be denoted by $\eta [0]$ and $\eta [1]$, respectively. By considering an

*N*

_{in}×

*N*matrix $B\u2208\mathbb{R}Nin\xd7NB$ consisting of the

_{B}*N*basis elements sampled at the

_{B}*N*

_{in}grid points, the initial value $\xi $ can be expanded by the basis as

where $Bdiag=blkdiag(B,B)\u2208\mathbb{R}2Nin\xd72NB$. Based on these notations, the proposed method is rewritten as a problem of finding the optimal coefficient $\eta $ by solving the following optimization problem:

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

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 (*N _{B}* <

*N*

_{in}).

### D. Matrix-free implementation for reducing memory usage

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 method^{36–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. *S*_{diag} and *E* can also easily be written as the functions, where *S*_{diag} picks up the *nN*_{obs} elements from a 2*nN _{p}*-dimensional vector whose indices coincide with those corresponding to the observation points, and

*E*substitutes a 2

*N*

_{in}vector into the 2

*N*-dimensional vector consisting of all zeros (see Fig. 12 in Appendix B). In addition, Φ

_{p}_{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

*B*

_{diag}as a function, too. Yet, storing

*B*

_{diag}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. *W ^{T}* and $f*$ are simple because

*W*

^{T}^{ }=

*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(hX\u22124)$ and $O(hX\u22126)$, respectively]. In contrast, the matrix-free method can be implemented with the realistic memory usage, reduced to $O(hX\u22122)$ (2D) and $O(hX\u22123)$ (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.

Size of sound field . | 8 m × 8 m . |
---|---|

Shape of region for initial value $\xi $ | Circular disk |

Placement of observation points $pobs$ | Circle |

Sound speed c | 340 m/s |

Width of Cartesian grid h _{X} | 0.02 m |

Size of time step h _{T} | 1/30 000 s |

Total number of time steps n | 360 steps |

Length of Fourier transform N_{DFT} | 512 samples |

Size of sound field . | 8 m × 8 m . |
---|---|

Shape of region for initial value $\xi $ | Circular disk |

Placement of observation points $pobs$ | Circle |

Sound speed c | 340 m/s |

Width of Cartesian grid h _{X} | 0.02 m |

Size of time step h _{T} | 1/30 000 s |

Total number of time steps n | 360 steps |

Length of Fourier transform N_{DFT} | 512 samples |

. | 2D . | 3D . | ||
---|---|---|---|---|

Width of Cartesian grid h _{X} | 0.02 | 0.01 | 0.02 | 0.01 |

No. elements in $\xi $ | 309 | 1249 | 4145 | 33377 |

Memory for $\xi \u2009(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\u2009(MB)$ | 0.34 | 0.57 | 10.79 | 30.56 |

No. grid points of FDTD | 401^{2} | 801^{2} | 401^{3} | 801^{3} |

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 |

. | 2D . | 3D . | ||
---|---|---|---|---|

Width of Cartesian grid h _{X} | 0.02 | 0.01 | 0.02 | 0.01 |

No. elements in $\xi $ | 309 | 1249 | 4145 | 33377 |

Memory for $\xi \u2009(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\u2009(MB)$ | 0.34 | 0.57 | 10.79 | 30.56 |

No. grid points of FDTD | 401^{2} | 801^{2} | 401^{3} | 801^{3} |

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 |

## IV. NUMERICAL EXPERIMENT

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.

### A. Examples of directivity diagrams and wave propagation

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 *r*_{init} was set to 0.2 m, the radius of the observation points *r*_{obs} was 0.3 m, and the number of the basis elements *N _{B}* 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.

### B. Properties of proposed method

For investigating the properties of the proposed method, the relation between its parameters (*r*_{init}, *r*_{obs}, and *N*_{B}) and the obtained results are shown in this section. To quantitatively assess the approximation quality, the relative energy $\u03d2$ is defined as follows:

where $\chi des\u2208{0,1}NDFTNobs\xd7NDFTNobs$ 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 *r*_{eval} was 1.5 m. This $\u03d2$ ∈ [0,1] becomes 1 if all energy is concentrated inside the desired part and 0 if all energy is outside it.

At first, the relation between the radius of the initial value *r*_{init} 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 *r*_{init}. To do so, the basis elements were chosen automatically for each *r*_{init} 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 *r*_{init} are indicated by the circles in Fig. 6. *r*_{init} was changed from 0.05 to 0.5 m, while *r*_{obs} was set to (*r*_{init} + 0.1) m. The directivity diagrams of the obtained results are shown in Fig. 7, and their relative energies $\u03d2$ are summarized in Fig. 8(i). These figures indicate that the radius of the initial value is important for representing a complex directivity.

Then, the relation between the radius of the observation points *r*_{obs} and the obtained result was investigated by fixing the other parameters, *r*_{init}(= 0.2 m) and *N _{B}* (= 50). The results are shown in Fig. 8(ii), which indicates that

*r*

_{obs}has little effect on the approximation accuracy. Note that a larger radius

*r*

_{obs}may increase the number of observation points

*N*

_{obs}, 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 *N _{B}* and the obtained result was investigated by fixing

*r*

_{init}(= 0.2 m) and

*r*

_{obs}(= 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.

In summary, the radius of the initial value *r*_{init} and the number of the basis elements *N _{B}* are important parameters for the proposed method, while the radius of the observation points

*r*

_{obs}has little impact. Therefore, both

*r*

_{init}and

*N*should be set sufficiently large for a given directivity, while the sufficient number of these values seems to depend on the directivity.

_{B}### C. Comparison with conventional time-domain method

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, *N _{B}* = 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

*r*

_{init}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 *N _{B}* = 29. That is, comparison between the top and bottom rows is fair in terms of the angular frequency components.

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.

## V. CONCLUSIONS

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.

### APPENDIX A: ANOTHER EXAMPLE OF FDTD SCHEME

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:

where *ρ* is density, *κ* is bulk modulus, ** u** and

**are the particle velocities in**

*v**x*and

*y*directions, respectively, and

*D*and

_{x}*D*are matrices corresponding to the finite difference approximations of the first-order spatial derivatives.

_{y}This FDTD scheme can also be written in the matrix-vector form in Eq. (8). Let the vector $\zeta [n]$ be defined as

Then, the matrix Φ, defined as

provides the update $\zeta [n+1]=\Phi \zeta [n]$ as in Eq. (8), where

### APPENDIX B: MATRIX-FREE CALCULATION

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 *S*_{diag}Φ_{rep}*EB*_{diag}, and then apply the combination of the Fourier transform and its inverse with the weighting. After them, multiplication of $BdiagTET\Phi 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.

Input: $\eta ,B,E,S$ |

Output: $pobs$ |

1: Separate $\eta $ into $\eta [1]$ and $\eta [0]$ |

2: $pin[1]=B\eta [1]$, $pin[0]=B\eta [0]$ |

3: $p[1]=zeropad(pin[1],E)$, $p[0]=zeropad(pin[0],E)$ |

4: for k = 1, 2,…, n − 1 do |

5: $p[k+1]=FDTD(p[k],p[k\u22121])$ |

6: $pobs[k+1]=select(p[k+1],S)$ |

7: end for |

8: $pobs=[pobs[1]T,pobs[2]T,\u2026,pobs[n]T]T$ |

Input: $\eta ,B,E,S$ |

Output: $pobs$ |

1: Separate $\eta $ into $\eta [1]$ and $\eta [0]$ |

2: $pin[1]=B\eta [1]$, $pin[0]=B\eta [0]$ |

3: $p[1]=zeropad(pin[1],E)$, $p[0]=zeropad(pin[0],E)$ |

4: for k = 1, 2,…, n − 1 do |

5: $p[k+1]=FDTD(p[k],p[k\u22121])$ |

6: $pobs[k+1]=select(p[k+1],S)$ |

7: end for |

8: $pobs=[pobs[1]T,pobs[2]T,\u2026,pobs[n]T]T$ |

Input: $pobs,E,S$ |

Output: $\eta \u0302$ |

1: Separate $pobs$ into $pobs[1],pobs[2],\u2026,pobs[n]$ |

2: Set $p\u0302[1]=0\u2208\mathbb{R}Np$, $p\u0302[0]=0\u2208\mathbb{R}Np$ |

3: for k = 1, 2,…, n − 1 do |

4: $p\u0302[1]=p\u0302[1]+zeropad(pobs[n\u2212k+1],S)$ |

5: $q=adjointFDTD(p\u0302[1],p\u0302[0])$ |

6: $p\u0302[0]=\u2212p\u0302[1]$, $p\u0302[1]=q$ |

7: end for |

8: $p\u0302in[1]=select(p\u0302[1],E)$, $p\u0302in[0]=select(p\u0302[0],E)$ |

9: $\eta \u0302[1]=BTp\u0302in[1]$, $\eta \u0302[0]=BTp\u0302in[0]$ |

10: $\eta \u0302=[\eta \u0302[1]T,\eta \u0302[0]T]T$ |

Input: $pobs,E,S$ |

Output: $\eta \u0302$ |

1: Separate $pobs$ into $pobs[1],pobs[2],\u2026,pobs[n]$ |

2: Set $p\u0302[1]=0\u2208\mathbb{R}Np$, $p\u0302[0]=0\u2208\mathbb{R}Np$ |

3: for k = 1, 2,…, n − 1 do |

4: $p\u0302[1]=p\u0302[1]+zeropad(pobs[n\u2212k+1],S)$ |

5: $q=adjointFDTD(p\u0302[1],p\u0302[0])$ |

6: $p\u0302[0]=\u2212p\u0302[1]$, $p\u0302[1]=q$ |

7: end for |

8: $p\u0302in[1]=select(p\u0302[1],E)$, $p\u0302in[0]=select(p\u0302[0],E)$ |

9: $\eta \u0302[1]=BTp\u0302in[1]$, $\eta \u0302[0]=BTp\u0302in[0]$ |

10: $\eta \u0302=[\eta \u0302[1]T,\eta \u0302[0]T]T$ |

Input: $p[n],p[n\u22121]$ |

Output: $p[n+1]$ |

1: for all k, l do |

2: $pk,l[n+1]=c2(hT/hX)2(pk+1,l[n]\u22122pk,l[n]+pk\u22121,l[n]+\u2009pk,l+1[n]\u22122pk,l[n]+pk,l\u22121[n])+2pk,l[n]\u2212pk,l[n\u22121]$ |

3: end for |

Input: $p[n],p[n\u22121]$ |

Output: $p[n+1]$ |

1: for all k, l do |

2: $pk,l[n+1]=c2(hT/hX)2(pk+1,l[n]\u22122pk,l[n]+pk\u22121,l[n]+\u2009pk,l+1[n]\u22122pk,l[n]+pk,l\u22121[n])+2pk,l[n]\u2212pk,l[n\u22121]$ |

3: end for |

Input: $p\u0302[1],p\u0302[0]$ |

Output: $q$ |

1: for all k, l do |

2: $qk,l=c2(hT/hX)2(p\u0302k+1,l[1]\u22122p\u0302k,l[1]+p\u0302k\u22121,l[1]+\u2009p\u0302k,l+1[1]\u22122p\u0302k,l[1]+p\u0302k,l\u22121[1])+2p\u0302k,l[1]+p\u0302k,l[0]$ |

3: end for |

Input: $p\u0302[1],p\u0302[0]$ |

Output: $q$ |

1: for all k, l do |

2: $qk,l=c2(hT/hX)2(p\u0302k+1,l[1]\u22122p\u0302k,l[1]+p\u0302k\u22121,l[1]+\u2009p\u0302k,l+1[1]\u22122p\u0302k,l[1]+p\u0302k,l\u22121[1])+2p\u0302k,l[1]+p\u0302k,l[0]$ |

3: end for |