Fourier spectral method for acoustic simulation
with domain enclosed by curved boundary
Yu KOHASE(1) , Tsubasa KUSANO(2) , Kohei YATABE(3) , Yasuhiro OIKAWA(4)(1)

Waseda University, Japan, y.kohase@asagi.waseda.jp(2)

Waseda University, Japan, tsubasa.k@suou.waseda.jp(3)

Waseda University, Japan, k.yatabe@asagi.waseda.jp(4)

Waseda University, Japan, yoikawa@waseda.jp

Abstract
Curved geometries often make sound propagation complex. Such complexity may cause acoustical problems
including flutter echo and sound focusing. When designing the geometry, acoustic simulation can be helpful
to prevent such problems. Fourier spectral method (FSM) is a simulation method by approximation using theFourier basis. Although FSM has many advantages such as its high convergence rate, the application range ofthe conventional FSM is limited to a simply-shaped domain with a specific boundary condition. In a realistic
acoustics room setting, there exist a lot of room shapes beyond the scope of FSM. At the same time, the Fourier
extension (FE) has been proposed for approximating a function on a complicated domain by the Fourier basis.
It can be expected that the FE expands the application range of FSM. In this paper, we introduce FE into FSM
for extending it to make a computational domain enclosed by curved boundaries tractable.
Keywords: Fourier spectral method, curved boundary, Fourier extension, function approximation

1INTRODUCTION

Acoustic wave-based simulation has been widely studied for predicting or understanding acoustical phenomena.
In wave-based method, a solution of the wave equation is approximated by arbitrary polynomial, which explains
the wave properties well. Then, they are studied by finite difference time domain method (FDTD) [1–3] andfinite element method (FEM) [4].
In the wave-based simulation for acoustics field, a complex-shaped domain is often focused because it makes
sound field complex, and its prediction is difficult more than simple one. Sometimes curved boundary, associated with a complex-shaped domain, causes acoustical problems such as flutter echo or sound focusing. Some
researches and simulations are conducted for preventing such problems [3, 5–7].
Fourier spectral method (FSM) is a simulation method by approximation using the Fourier basis. FSM is often
used for fluid mechanics and earth science, and it has several advantages including its high convergence rate, but
the application range of the conventional FSM is limited to a simply-shaped domain with a specific boundary
condition because of characteristics of the Fourier basis. Since many real objects have a complex-shaped interior
space, the conventional FSM cannot be applied to such space [8–12].
At the same time, the Fourier extension (FE) has been proposed for approximating a function on a complicateddomain by the Fourier basis. It can approximate a nonperiodic function which is treated as a periodic function
in an extended domain. Then, it can be expected that FE expands the application range of FSM [13, 14].
In this paper, we introduce FE into FSM for handling a domain enclosed by curved boundary. Since this is a
first step, the Dirichlet boundary condition is considered in this paper. As a result, FSM becomes possible tohandle a curved boundary defined on non-grid points. Some influence of parameter selection is evaluated anddiscussed by numerical simulation on a calabash-shaped domain.

4498

2

FOURIER SPECTRAL METHOD

Spectral method is a numerical simulation method with high accuracy, because it globally approximates a function by high order polynomials. In spectral method, an unknown function u(x) is approximated by N basis
functions φ(x),N−1
Xu(x) ≈ uN (x) =ak φk (x).
(1)
k=0

One example of basis functions φ(x) is the Fourier basis used in the inverse discrete Fourier transform (iDFT),u(x) =N−1X2πû(k)eikx N ,(2)k=0where û(k) is Fourier transform of u(x), i means imaginary unit, and k means wave number. In this case, thedifferential operator can be expressed asN −1 X
∂nu
=∂xn

k=0

2πikN

n

2πû(k)eikx N .

(3)This derivative is more accurate than the other methods as the finite difference method. In addition, it can beefficiently computed by the fast Fourier transform (FFT) algorithm. This differentiation via the Fourier bases isused to approximate the spatial derivative of a partial differential equation. However, owing to the periodicity of
the Fourier bases, the standard spectral method can only be applied to simple domain and boundary condition,
which is the main limitation [8–11].
2.1 Example of Fourier Spectral Method
In acoustic simulation, distribution of sound pressure can be expressed by
P n = P (n∆t),(4)where ∆t is discrete time steps, n ∈ N is the time index. In this paper, a grid points set constructing the P ison the Cartesian grid. Let the trapezoidal rule is considered for approximating the time integration. At first, thewave equation is denoted by∂2P= c2 4 P,
(5)∂t2
where c means sound speed, and 4 means the Laplace operator. When v is partial derivative of P with respect
to t (v = ∂P/∂t),P n+1 − P n
1 n+1=v+ vn(6)∆t2
is obtained by discrete integration using the trapezoidal rule. In this paper, initial condition of v (condition ofv at t = 0) is 0. Then, the wave equation can be represented using v:
v n+1 − v nc2=LP n+1 + LP n ,∆t2(7)

where L is the second order spatial differential operator approximating the Laplace operator. Then, a scheme isobtained by Eq. (6) and Eq. (7),((2 )
2 )
∆tc
∆tc
1−L P n+1 = 1 +L P n + ∆tv n .
(8)224499Finally, Eq. (8) is calculated in FSM by,("(#2 )2 )N−1N−1X
X
2π
∆tc
∆tc
n+1ikx 2πnn1−L̂ P̂(k)e N =1+L̂ P̂ (k) + ∆tv̂ (k) eikx N ,22
k=0

(9)

k=0where L̂ is the frequency-domain representation of L [14].

3FOURIER EXTENSION

Because of the Fourier bases’s periodicity, simulations over a complex-shaped domain are difficult for the standard FSM. FE is a technique to approximate a nonperiodic function f on [0, 1]2 by the Fourier bases, wherethe periodicity of the Fourier bases is treated by extending the computational domain to [0, T ]2 , and T (> 1) isan extension parameter. Therefore, FE is expected to widen the application range of the FSM.
In this paper, we consider two-dimensional case. At first, we define some domains shown in Fig. 1. In thisfigure, PR and PR̂ are extended spatial domain and frequency domain, respectively, PΩ and PΛ are set of NΩ
grid points and set of NΛ bases, respectively, and a red line means the boundary PδΩ . Then, the main problemof FE is formulated by a least squares problem,
2a = argmin

Xc∈CNΛx∈PΩXf (x) −ck φk (x)T

2πφk (x) = eik x T,(10)k∈PΛ

where kT is the transpose of k, c is the Fourier coefficients, and a is the solution set of c. The solution a canbe found by collocation, through solving the rectangular system,
Aa = b, A ∈ CNΩ ×NΛ ,b ∈ CNΩ ,(11)where A is the Fourier extension operator, a subblock of multi-dimensional unitary iDFT matrix. Then, AS
means the Fourier extension operator over spatial domain S. In addition, the boundary of the computational
domain is contained inside of the extended domain. In this technique, over-sampling rate q = NΩ /NΛ is one ofimportant parameter, q > 1 is necessary for better accuracy [13, 14].

𝑃#𝑃#$𝑃%𝑃"𝑎 Spatial domain𝑏 Frequency domain

Figure 1. Relation between the extended domain and original domain. The spatial domain Ω encompassing thegrid points set PΩ , and the frequency domain Λ encompassing the frequency points set PΛ .

4500

3.1 Fourier spectral method using Fourier extension
We introduce FE into FSM. According to Eq (9), this matrix system is nAΩ
no−12 "(# AΩ I +2 )L̂AδΩ2 I − ∆tc2∆tcI−L̂ P̂n+1 = ..

2.
no−1 2AδΩk I − ∆tcL̂2o∆tc 2L̂2wh2..
.

P̂n + ∆tvn 
,


(12)

whk

where h means a boundary value, and w means a weight for boundary domain δΩ. In this paper, w = 1 forsimplicity.

4NUMERICAL EXPERIMENT

4.1 Experimental Condition
In this paper, we calculate a sound field enclosed by the curved boundary shown in Fig. 1, which is given byπr(θ) = 30 cos(2θ + ) + 2.5 ,2

θ = [0, 2π),(13)and set the homogeneous Dirichlet boundary condition to it. Other experimental conditions are summarized inTable 1, and ζb is defined as ratio of points on boundary NδΩ to NΩ ,ζb =NδΩ
.
NΩ

(14)

Table 1. Simulation condition.
NR216

NΩ

214

NΛ

26 , 28 , 210ζb1, 5, 10

Spatial discretization interval [m]1Time discretization interval [ms]2Sound speed [m/s]340

4.2 Evaluation Error
In this paper, we evaluated three kinds of errors; Extension error, boundary error and energy error. In FE
process, we solved the matrix equation using iterative method with error tolerance 1.0 × 10−10 . we defined theextension error,kb − AΩ ak2
EΩ =.
(15)
kbk2
where k · kp is the `p -norm. In calculation without error, boundary value on RδΩ is constant in every time steps.
In this calculation, priority of boundary value in optimization is controlled by boundary weight. We defined theboundary error,EδΩ = kh − AδΩ ak2 .
(16)

4501

This system have no theoretical decay factor, energy Q around RΩ must be conserved. The energy aroundspatial domain is defined as(τ + σ)Q=,(17)
Qmax
where2

τ=1 ∂P,2 ∂t 2σ=c2
2k∇Pk2 ,2(18)and Qmax means the max value of Q over all the simulation time steps [15].
4.3 Simulation
The typical results of the wave propagation are shown in Fig. 2 with NΛ = 28 and ζb = 5. Fig. 2(a) illustrates
the initial condition of this simulation, (b)-(f) show the radiation of the wave with time evolution.
At first, we performed simulation, with changing ζb , then, the obtained results are shown in Fig. 3. In thisresults, when ζb is 1, EδΩ is higher than other cases, and Q is decreasing with time evolution. The cases ofζb = 5 and ζb = 10 suggest similar trends, but ζb = 10 is more accurate than ζb = 5.
In addition, the results of changing NΛ with ζb = 5 are in Fig. 4. According to Fig. 4(b), a trends of EδΩare almost same over every parameters. NΛ = 26 is best condition about EΩ , but its energy is considerably
decreasing with time evolution. Then, the case of NΛ = 28 has smaller value in EΩ than the case of NΛ = 210 ,but about energy Q conservation, the case of NΛ = 210 is better than the case of NΛ = 28 . According to Fig. 3and Fig. 4, ζb has influence of the boundary error EδΩ , and NΩ has influence of the extension error EΩ .
250

250 250

250

200

200

200

150

150

150

100

1005050200

1005000500100 150 200 250(a) 𝑡 = 0050100 150 200 250

00.30.20500.10(c) 𝑡 = 0.4(b) 𝑡 = 0.2150100 150 200 250

250250250-0.1200200200-0.2150150

1001005050150100100500050100 150 200 250(d) 𝑡 = 0.6050050

100 150 200 2500(e) 𝑡 = 0.8-0.3050

100 150 200 250(f) 𝑡 = 1.0Figure 2. Snapshots of the P at t seconds, and their experimental conditions are NΛ = 28 and ζb = 5.

05050

100150200250

CONCLUSION

In this paper, we introduced FE into the FSM, and implemented simulations over the domain enclosed by thecurved boundary. At first, the scheme of FSM for the wave equation is formulated by trapezoidal rule, and FE
as a least squares problem is defined. In addition, the matrix system should be calculated directly was obtainedby above all. Finally, we evaluated three kinds of error with changing dominant parameters for accuracy, andwe confirmed the existence of parameter combinations having higher accuracy, NΛ = 28 , ζ = 5 are one of4502

better combination for parameters over this paper. Then, ζb has influence of boundary error EδΩ , and NΩ hasinfluence of extension error EΩ .
10.015

10.02

0.8

0.6

Q0.005

00.25

0.5

0.75

01Time [s]QEE0.60.01

0.005

00.80.015

0.01

0.25

0.50.75b0.2bb0.4
00.4

0100.25

0.75

=5
= 101

Time [s]Time [s](a) Extension error

0.5=1(b) Boundary error(c) Energy

0.2

Figure 3. Experimental results with NΛ = 28 and changing ζb .
00.015

00.2
10.02

0.6

0.8

1QEE0.6 0.80.01

N = 260.40.005
0.005

N = 280.60.2

N = 210

Q01Time [s]

0.8

0.015

0.01

0.4

00.25

0.50.75

1000.25

0.50.75

Time [s]Time [s](a) Extension error(b) Boundary error

100 0.4

0.25

0.50.75

1Time [s](c) Energy

0.2

Figure 4. Experimental results with ζb = 5 and changing NΛ .
000.2

0.4

0.60.8

1Time [s]REFERENCES[1] B. Hamilton, “Finite difference and finite volume methods for wave-based modelling of room acoustics,”Ph.D. dissertation, Acoustics and Audio Group, Univ. Edinburgh (2016).
[2] S. Bilbao, “Modeling of complex geometries and boundary conditions in finite difference/finite volume time
domain room acoustics simulation,“ IEEE Trans. Audio, Speech, Lang. Process. 21, 1524-1533 (2013).
[3] D. Takeuchi, K. Yatabe and Y. Oikawa, “Source directivity approximation for finite-difference time-domain
simulation by estimating initial value,“ J. Acoust. Soc. Am. 145, 2638-2649 (2019).
[4] T. Okuzono, T. Otsuru and K. Sakagami, “Applicability of an explicit time-domain finite-element method
on room acoustics simulation,“ Acoust. Sci. & Tech. 36, 377-380 (2015).
[5] L. Savioja, and N. Xiang, “Introduction to the Special Issue on Room Acoustic Modeling and Auralization,“ J. Acoust. Soc. Am. 145, 2597-2600 (2019).
[6] Y. Lu and H. Song, “Forecasting flutter echo: The role of largest local Lyapunov exponents,“ Noise ControlEngr. J. 67, 56-68 (2019).

4503

[7] S. Yon, M. Tanter, and M. Fink, “Sound focusing in rooms: The time-reversal approach, “J. Acoust. Soc.
Am. 113,1533-1543 (2003).
[8] J. P. Boyd, “Chebyshev and Fourier Spectral Methods,“ Dover Publications Inc. (USA), second edition(2000).
[9] L. N. Trefethen, “Spectral Methods in MATLAB,“ SIAM (USA), (2000).
[10] J. Shen, T. Tang and L.-L. Wang, “Spectral Methods Algorithms, Analysis and Applications, “ SpringerVerlag (DE), (2011).
[11] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, “Spectral Methods in Fluid Dynamics, “ SpringerVerlag (DE), (1988).
[12] L. Q. Chen and J. Shen, “Applications of semi-implicit Fourier-spectral method to phase field equations,“Comput. Phys. Commun., 108, 147-158 (1998).
[13] R. Matthysen and D. Huybrechs, “Function approximation on arbitrary domains using Fourier extension
frames,“ arXiv preprint arXiv:1706.04848 (2017).
[14] R. Matthysen, “Introducing Redundancy into Numerical Computations,“ Ph.D. dissertation, KU Leuven(2018).
[15] S. Bilbao, “Numerical Sound Synthesis,“. John Wiley and Sons(UK), (2009).

4504