3. セミトランスオーラルと音場模擬技術に基づく
3.5 有限差分時間領域法の高精度化の提案
3.5.7 計算式の導出
極座標表現された時空間における逐次計算式を導く．式(105)に変形された波動 方程式は因数分解によって次式に変形できる．
(∂2u
∂t2 − 1 c2
∂2u
∂x2 )
= (∂u
∂t + 1 c
∂u
∂x
) (∂u
∂t − 1 c
∂u
∂x )
= 0. (113)
このとき，直交座標系t, xを関係式(106),(107)および極座標系と直交座標系にお けるヤコビアンを適用することで，次式が得られる．
(
sinθ− 1 ccosθ
) ∂
∂r (
sinθ+ 1 ccosθ
) ∂
∂r +
(
sinθ− 1 ccosθ
) ∂
∂r 1 r
(
cosθ− 1 csinθ
) ∂
∂θ + 1
r (
cosθ+ 1 csinθ
) ∂
∂θ (
sinθ− 1 ccosθ
) ∂
∂r + 1
r (
cosθ+ 1 csinθ
) ∂
∂θ 1 r
(
cosθ−1 csinθ
) ∂
∂θ
= 0. (114)
この式を展開することで，次式が得られる．
(
sin2θ− 1 c2cos2θ
)∂2µ
∂r2 + 1 r
(
cos2θ− 1 c2 sin2θ
)∂µ
∂r + 2
r (
cosθsinθ+ 1
c2 cosθsinθ )∂µ
∂r
∂µ
∂θ
− 2 c2
(
cosθsinθ+ 1
c2cosθsinθ ) ∂
∂θ + 1
r2 (
cos2θ− 1 c2sin2θ
) ∂2
∂θ2 = 0. (115)
このとき，2倍角の公式を用いる事で次式を得る．
1
2(cn−cpcos 2θ)∂2µ
∂r2 + 1
2r(cn+cpcos 2θ)∂µ
∂r + 2
rcpsin 2θ∂µ
∂r
∂µ
∂θ − 1
r2cpsin 2θ∂µ
∂θ
+ 1
2r2 (cn+ccos 2θ)∂2µ
∂θ2 = 0, (116)
cp = 1 + 1
c2, (117)
cn = 1− 1
c2. (118)
時間と空間からなる音源からのユークリッドノルムrに沿ったインパルス応答の 更新は，式(116)に於ける時空間の位相θのみで行うことにより達成される．しか
し，式(116)はθに関する1次導関数，2次導関数およびrに関する1次導関数との
合成関数が存在するため，このままでは式(109)形式への変形はできない．そこで本 節では，式(116)を位相θについての周波数領域ωθに写像し，ωθを引数とする関数 の線形和として解く．フーリエ変換を用いてθを周波数領域に写像することで，位 相領域の偏微分は係数の乗算となる．加えて，位相領域は−πからπの定義域で周 期性を伴うため，離散フーリエ変換の周期拡張に基づく巡回生が結果に影響を及ぼ さない．
U(r, ωθ) = Fθ[µ(r, θ)], (119)
−jωθU(r, θ) = Fθ
[ ∂
∂θµ(r, θ) ]
, (120)
ωθ2U(r, θ) = Fθ
[ ∂2
∂θ2µ(r, θ) ]
. (121)
上記の関係式を式(116)に適用することで，極座標表現された時空間におけるユー クリッドノルムrと位相周波数領域θωについての関係式を構築できる．ただし式 (116)において，例えば第4項においてはsin 2θと∂µ/∂θというθに関する関数の積 となっており，位相周波数ωθ領域に写像する場合には周波数領域に写像された両関 数の畳み込みとなる．このとき，位相領域の関数sin 2θとcos 2θは離散フーリエ変 換の利用により下記の関数に写像される．
1
2[δ(ωθ−2δωθ) +δ(ωθ+ 2δωθ)] = Fθ[cos 2θ], (122) j
2[δ(ωθ−2δωθ)−δ(ωθ+ 2δωθ)] = Fθ[sin 2θ]. (123) なお上式のδ(·)はデルタ関数を示す．また時空間における音源からのユークリッ ドノルムrの偏微分は，更新式のためコンパクト差分による近似を行う．
1
2∆r[U(r+ ∆r, ωθ)−U(r−∆r, ωθ)]
= Udr(r, ωθ) (
≃ ∂
∂rU(r, ωθ) )
, (124)
1
∆r2 [U(r+ ∆r, ωθ) +U(r−∆r, ωθ)−2U(r, ωθ)]
= Uar(r, ωθ) (
≃ ∂2
∂r2U(r, ωθ) )
. (125)
得られた式(122)∼(125)を式(116)に適用することで，式(126)が得られる．
C1⊗ 1
∆r2 [U(r+ ∆r, ωθ) +U(r−∆r, ωθ)−2U(r, ωθ)]
+C2⊗ 1
2∆r[U(r+ ∆r, ωθ)−U(r−∆r, ωθ)]
+C3⊗jωθ 1
2∆r[U(r+ ∆r, ωθ)−U(r−∆r, ωθ)]
+C4⊗jωθU(r, ωθ)
+C5⊗ −ω2θU(r, ωθ) = 0. (126)
なお，⊗は時空間における位相周波数領域の畳み込み積分を示し，C1 ∼C5はコン パクト差分で近似されたユークリッドノルム領域の更新式に対する係数として，式
(127)∼(131)で定義される．
C1 = 1
2[cpδ(ωθ−2∆ωθ) + 2cnδ(ωθ)−cpδ(ωθ+ 2∆ωθ)], (127) C2 = 1
2r[cpδ(ωθ−2∆ωθ) + 2δ(ωθ) +cpδ(ωθ+ 2∆ωθ)], (128) C3 = jcp
r [δ(ωθ−2∆ωθ)−δ(ωθ+ 2∆ωθ)], (129) C4 = jcp
r2 [δ(ωθ−2∆ωθ) +δ(ωθ+ 2∆ωθ)], (130) C5 = 1
r2[cp
2δ(ωθ−2∆ωθ) +cnδ(ωθ) + cp
2δ(ωθ+ 2∆ωθ)]. (131)
式(126)では，時空間における原点からのユークリッドノルムrに対する位相周波
数領域の音圧分布U(r, ωθ)の更新を，ユークリッドノルムrに対するコンパクト差 分だけで解いており，位相領域の導関数は位相周波数領域にて計算されるため数値 分散を生じない．また，位相周波数領域への写像には離散フーリエ変換が用いられ るが，位相領域は−πからπまでの値域で周期性を伴うため，離散フーリエ変換の 要求する周期性に基づく問題を回避できる．