3. セミトランスオーラルと音場模擬技術に基づく
3.5 有限差分時間領域法の高精度化の提案
3.5.1 数値分散による波束の乱れ
実空間を伝播する音響信号は式(72)の波動方程式によって定式化される．なお，
式(72)の‘u’は空間座標‘x, y, z’および時刻‘t’を引数とした時空間の音圧分布を示 す関数であり，‘c’は音速を示す．音場模擬の一手法である有限差分時間領域法は次 式の波動方程式を解くことで対象音場のインパルス応答を得る．
∂2u
∂t2 = 1 c2
(∂2u
∂x2 + ∂2u
∂y2 +∂2u
∂z2 )
. (72)
このように音圧に対する波動方程式は2次導関数に基づく偏微分方程式で表現さ れ，微分は差分式によって近似が可能である．任意の関数f(x)の引数xについての 1次導関数は, 引数xの有限な変化量∆xに対する関数f(x)の変化量を用いて，式 (73)のように風上差分，風下差分にて近似できる．
∂
∂xf(x)≃ f(x+ ∆x)−f(x)
∆x ≃ f(x)−f(x−∆x)
∆x . (73)
加えて風上差分と風下差分によって得られた式を元に差分式を構成することで，2 次導関数を式(75)のように差分式で近似できる．
∂2
∂x2f(x) = ∂
∂x ( ∂
∂xf(x) )
, (74)
≃ 1
∆x
[f(x+ ∆x)−f(x)
∆x − f(x)−f(x−∆x)
∆x
]
. (75)
有限差分時間領域法では時間を標本化間隔∆t,空間を標本化間隔∆x,∆y,∆zにて 時空間の音圧分布を示す関数u(t, x, y, z)を離散化する．この離散化により式(76)の
ような，時間の進行に対する空間の音圧分布の更新式を構築できる．なお音圧が単位 時間∆tに伝播する距離c∆tが，空間の単位距離∆xを超えると増幅率の問題[116]か ら演算結果が発散する．このことから時空間の標本化では標本化間隔がc∆t/∆x≤1 を満たす条件が課せられる．
u(t+ ∆t, x, y, z) = c2
2[cx(t, x, y, z) +cy(t, x, y, z) +cz(t, x, y, z)]
+2u(t, x, y, z)−u(t−∆t, x, y, z), (76) cx(t, x, y, z) = u(t, x+ ∆x, y, z) +u(t, x−∆x, y, z)
−2u(t, x, y, z), (77)
cy(t, x, y, z) = u(t, x, y+ ∆y, z) +u(t, x, y−∆y, z)
−2u(t, x, y, z), (78)
cz(t, x, y, z) = u(t, x, y, z+ ∆z) +u(t, x, y, z−∆z)
−2u(t, x, y, z). (79)
有限差分時間領域法ではこの更新式によって，音場(音場内の音響信号の伝播)を 計算機で模擬することができる．模擬では初期状態と各時刻の境界における音圧分 布を定義し，離散時間∆t進む毎に変化する空間の音圧分布を式(76)によって逐次 更新することで各時刻の音圧分布を得る．その際，音源座標の音圧を放射される音 響信号によって加圧し，受聴座標の音圧を各時刻に観測することで，これら座標間 のインパルス応答を得ることができる．
しかし有限差分時間領域法では数値分散が生じることが知られている．数値分散 が生じることによって，伝播した音響信号から高周波数帯域の伝播が遅れ，本来同 一波形として伝播するはずの音響信号が波束を乱しつつ伝播する．図28(a)は帯域 制限されたインパルス(sinc (t, x))を示しており，全周波数帯域の位相は0の左右対 称な波形である．反射や拡散などを伴わない指向性のある伝播において，その応答 は原信号の平行移動sinc (t+ ∆t, x−c∆t)であることが波動方程式の一般解として 要求される．しかし有限差分時間領域法を用いてこの音響信号の伝播を模擬すると，
その応答は図28(b)に示されるようにsinc関数の形状をしておらず，波束，すなわ
ち群遅延が歪んでいることが認められる．これは有限差分時間領域法に知られる数 値分散という誤差に起因することが知られており，応答波形の歪みに関する主要な 原因となっている．有限差分時間領域手法は平易な音場模擬アルゴリズムであるが，
数値分散が生じることによって実空間と乖離した音響信号の伝播が模擬され，その 結果として有限差分時間領域法を用いた音場再現の品質低下を引き起こしている．
すなわち音場模擬の精度を改善させることが，高精度な音場再現技術には必要不可 欠である．
連続の式である波動方程式には2次の偏微分方程式である．有限差分時間領域法 では波動方程式の偏微分方程式を2次，4次などの差分式で近似する．この差分近似 が周波数軸で連続の式と異なる特性を示すこと，すなわち数値分散によって帯域毎 の位相速度が歪む問題が生じると考えられる．音圧分布u(t, x, y, z)の空間x軸の微 分について考察する．便宜上，1次元空間の音圧分布u′(x)について数値分散の考察 を行うが，3次元空間の音圧分布u(t, x, y, z)の引数t, x, y, zは直交するため以下の 考察は全ての変数に対して適用可能である．1次元空間の音圧分布を示す関数u′(x) の空間領域の微分を差分近似する場合，その作用は関数u′(x)に2次精度の差分作 用素o2(x)の畳み込みとして考えられる(式(80))．なお，式(80)の記号‘∗’は関数の 畳み込み演算子を意味する．
u′(x)∗o2(x) = u′(x+ ∆x)−2u′(x) +u′(x−∆x), (80) o2(x) = δ(x+ ∆x)−2δ(x) +δ(x−∆x), (81) g(x)∗f(x) =
∫ ∞
−∞
g(x−λ)f(λ)dλ. (82)
空間領域xの畳み込みは，空間周波数ωx領域において乗算として書き下すことが でき，式(80)においてはo2(x)のxに関する空間周波数ωxを引数とした空間周波数 ωx領域の関数O2(ωx)，すなわち2−2 cos (ωx)の乗算と等価である．なお上記の式 はテイラー展開にて式(84)のように書き下せる．
O2(ωx) =
∫ ∞
−∞
o2(x)e−jωxxdx
= 2−2∗cos (ωx) (83)
≃ (
−ωx2+ ω4x 12 − ωx6
72 + ω8x 576· · ·
)
. (84)
つまり，任意の関数u′(x)に対する2次微分を2次の差分近似式で行う演算は，周 波数領域で式(84)の関数を乗算することと等価である．また同様に4次精度の差分 近似を行う場合においても，音圧分布を示す関数u′(x)と差分作用素o4(x)との畳 み込みとして表現できる(式(85)).
u′(x)∗o4(x) = 1
12[−u′(x+ 2∆x) + 16u′(x+ ∆x)−30u′(x)
+16u′(x−∆x)−u′(x−2∆x)], (85) o4(x) = 1
12[−δ(x+ 2∆x) + 16δ(x+ ∆x)−30δ(x)
+16δ(x−∆x)−δ(x−2∆x)]. (86) この畳み込みも2次精度の差分近似と同様に周波数領域において関数
(2 cos (2ωx)−32 cos (ωx) + 30)/12, (87) を乗算することと等価である．上記の式はテイラー展開にて式(90)のように書き下 せる．
O4(ωx) =
∫ ∞
−∞
o4(x)e−jωxxdx (88)
= 2 cos (2ωx)−32 cos (ωx)−30
12 (89)
≃ (
−ωx2− ωx6
90 + ωx8 1008· · ·
)
. (90)
つまり任意の関数u′(x)に対する2次微分を4次の差分近似式で行う演算は，周波 数領域で式(90)の関数を乗算することと等価である．一方で連続の式における理想 的な2次微分が周波数領域にて及ぼす影響を考察する．
u′(x) =
∫ ∞
−∞
U(ωx)ejωxxdωx, (91)
∂
∂xu′(x) =
∫ ∞
−∞
jωxU(ωx)ejωxxdωx, (92)
∂2
∂x2u′(x) =
∫ ∞
−∞−ω2xU(ωx)ejωxxdωx. (93) 関数u(x)はxの空間周波数ωxを引数とした周波数領域の関数U(ωx)の逆フーリエ 変換として扱うことができる(式(91))．そこで式(91)の両辺をxで微分すると，左
辺は関数u′(x)の1次導関数，右辺は逆フーリエ変換の回転子ejωxxに含まれるjωx
が係数としてU(ωx)にかかる．そのうちjは定数項であるが，変数ωxは逆フーリエ 変換の積分変数であるため逆フーリエ変換の積分演算内部に含まれる．これにより 1次微分が周波数領域にて及ぼす影響を式(92)として記述できる．加えて，上記の 式(92)の両辺を再度xで微分することにより，左辺ではu′(x)の2次導関数が得ら れるが，右辺は同様に逆フーリエ変換の回転子からjωxが係数として関数U(ωx)に かかる．これらの式により，関数u′(x)を空間領域xで2次の導関数を導出するこ とは，変数xの周波数ωx領域にて式−ωx2を乗算する演算と等価であることが示さ れる．変数−ω2xは逆フーリエ変換の積分変数であるため逆フーリエ変換の積分演算 内部に含まれる．このとき，2次元微分の理想的な周波数領域の係数−ω2は，2次 の差分近似における係数−ωx2+ωx4/12−ωx6/72 +ωx8/576· · · および4次の差分近似 における係数−ω2x−ωx6/90 +ωx8/1008· · · と異なる．すなわちこの周波数領域にお ける2次微分の係数の差により，有限差分時間領域法では波動方程式の偏微分方程 式を全周波数帯域において適切に解くことができず，数値分散誤差を引き起こして いるものと考えられる．連続の式，2次の差分近似式，4次の差分近似式を用いた2 次微分の周波数領域の作用素(周波数領域の係数)の比較を図29に示す．2次の差分 近似式では連続の式に対してω4x/4−ω6x/12 +ωx8/576· · · の誤差，4次の差分近似式 では連続の式に対してωx6/30 +ωx8/1008· · · の誤差を含んでいることがわかる．この 誤差は角周波数ωxに基づくことから，数値分散の誤差は高周波数帯域ほど大きいこ とが示される．