3. セミトランスオーラルと音場模擬技術に基づく
3.5 有限差分時間領域法の高精度化の提案
3.5.2 スペクトル法に基づく数値分散の抑圧
辺は関数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に基づくことから，数値分散の誤差は高周波数帯域ほど大きいこ とが示される．
間領域法の更新式を下記の式(94)として周波数領域の演算で示せる．
u(t+ ∆t, x, y, z) = c2
2[vx(t, x, y, z) +vy(t, x, y, z) +vz(t, x, y, z)]
+2u(t, x, y, z)−u(t−∆t, x, y, z), (94) vx(t, x, y, z) =
∫ ∞
−∞−ωx2Vx(t, ωx, y, z)ejωxxdωx, (95) vy(t, x, y, z) =
∫ ∞
−∞−ωy2Vy(t, x, ωy, z)ejωyydωy, (96) vz(t, x, y, z) =
∫ ∞
−∞−ωz2Vz(t, x, y, ωz)ejωzzdωz, (97) Vx(t, ωx, y, z) =
∫ ∞
−∞
u(t, x, y, z)e−jωxxdx, (98) Vy(t, x, ωy, z) =
∫ ∞
−∞
u(t, x, y, z)e−jωyydy, (99) Vz(t, x, y, ωz) =
∫ ∞
−∞
u(t, x, y, z)e−jωzzdz. (100) ここで空間領域の2次微分を3.5.1項に基づき式(95)∼(100)として導出すること により，従来の差分近似より精度を高めている．また更新方向である時間領域では 従来のコンパクト差分を用いて更新に参照される時刻を最小化することで，従来の 有限差分時間領域法と同等の簡潔な更新式を実現している．
従来の有限差分時間領域法では増幅率の問題から，標本化間隔∆t,∆xと速度cが
c∆t <∆xを満たすようにする条件があるため，位相速度を小さい値にすることが
要求されるが計算コストが発散する傾向があった．数値分散が強く発生する傾向が あった．
図30,31は空間周波数型FDTD法で用いられる，空間音圧分布の2次導関数を示す．
図30(a)は空間の中央付近にエネルギのある音場を示し，その2次導関数(図30(b))
が正確に導出されている．一方図31(a)は演算される空間の端にエネルギのある空 間音圧分布であり，空間周波数型FDTD法で導かれる2次導関数(図31(b))では，負 の空間座標の端に存在したエネルギが正の空間座標の端に回り込んでいる事が確認 できる．従来は音場模擬の対象となる空間に対して，倍の領域を窓関数で切り出し，
対象の空間外のエネルギを0とすることで回り込みを回避していた．しかし，対象 空間の境界でエネルギが急激に変化することで，対象空間の境界近傍における吸音
と反射を適切に扱えない問題があった．そこで次項では，時空間に極座標表現を適 用することで，離散フーリエ変換の要求する周期性が問題とならない計算モデルを 提案する．
G
1,( ω )
G ( ω )
3,L
G ( ω )
2,L
G ( ω )
4,R
G ( ω )
6,R
G ( ω )
5,R
G
2,R( ω ) G ( ω )
5,L
G ( ω )
3,R
G
6,L( ω )
G
1,R(ω ) G (ω
4,L)
(a) Without listener’s head
G (ω)
1,LG ( ω )
3,L
G
2,L( ω )
G (ω)
4,RG
6,R( ω )
G ( ω )
5,R
(b) With listener’s head
図24 頭部近傍スピーカアレイにおける各スピーカと両耳間の伝達関数
(a) Photograph
3 1
2
4
6 5
40˚
Head rest
(b) Overhead view
図25 頭部近傍スピーカアレイ
-1.5 -1 -0.5 0 0.5 1
0 50 100 150 200
Time[msec.]
Amplitude
(a) Left side ear
-1.5 -1 -0.5 0 0.5 1
0 50 100 150 200
Time[msec.]
Amplitude
(b) Right side ear
図26 左側のスピーカから両耳位置までのインパルス応答
x
仮想空間における模擬 実空間における提示
物体近傍において 模擬される波面
頭部近傍において 再現される波面
模擬されるスピーカユニット 実空間のスピーカユニット
無視される音響信号
仮想空間にて 到来する 音響信号
実環境にて 提示される音響信号 到来方向による
音響信号の選択
(a) (b)
図27 頭部に到来する音響信号の予測と提示
(a) Original
(b) Propagated response
図28 有限差分時間領域法によって波束の乱れた応答
0 π2 π
Frequency [rad]
−10
−8
−6
−4
−2 0
Gain
Analytic difference
Numerical difference (2nd-order) Numerical difference (4th-order)
図29 連続の微分と離散の差分の周波数領域において乗算される係数
(a) u(t, x, y)
(b) ∂x∂22u(t, x, y)
図30 中央にエネルギーのある空間音圧分布と離散フーリエ変換を用いて導出され
(a) u(t, x, y)
(b) ∂x∂22u(t, x, y)
図31 境界付近にエネルギーのある空間音圧分布と離散フーリエ変換を用いて導出
表8 数値分散に関する予備実験の計算機シミュレーション条件 Spatial size 5123point
Propagated signal Spherical gaussian wave CFL coefficient 1/2, 1/4, 1/8, 1/16, 1/32, 1/64 Propagated distance 128 points
Interpolation algorithm Cosine interpolation