3 次元計測と境界要素法による個人の頭部伝達関数の推定—並列計算機の利用—
高根昭一 ∗ ，
松橋太陽 † ，曽根敏夫 ∗∗秋田県立大学 システム科学技術学部 電子情報システム学科†
秋田県立大学大学院 システム科学技術研究科 電子情報システム学専攻
概要：3 次元計測された頭部形状をもとに，境界要素法を用いて個人の頭部伝達関数を推定する手法と，推定された頭部伝達関数の実測結果との比較について報告する．
境界要素法を無限領域の問題に適用するとき，解の一意性について配慮をする必要がある．その方法について述べ，プログラムとして情報シナジーセンターの並列計算機でインプリメントした．頭部伝達関数の推定結果は，実測結果と大局的な傾向に一致をみることができた．

1 はじめに
頭部伝達関数 (Head-Related Transfer Function，以下では HRTF と書く) は，頭の中心に相当する 1 点と，外耳の適切な点 (理想的には鼓膜上) の間の伝達関数であり，通常は自由空間 (音の反射のない空間) において定義される．HRTF は，両耳に到達する音から得る空間的な情報の知覚において重要な役割を果たすことが知られている1) ．また，ある方向の HRTF を音源信号に畳み込んで聴取者の両耳に提示すれば，音源信号がその方向から聞こえるように知覚させることが可能である2) ．現在では，ヘッド ホンを用いた 5.1ch のサラウンド 音場の再生などにこの原理が応用され，民生用の機器にも組込まれている3) ．また，著者らのグループでは，反射をもつ一般的な音場における両耳音圧を，聴取者の移動を考慮して提示できる聴覚ディスプレ イの原理を考案
し 4) ，その原理に基づいて実時間動作するシステムの基礎的なインプ リメンテーションを行っている5) ．この原理では，聴取者の全周方向の HRTF が必要となる．
定義からわかるとおり，HRTF は個人と音源の位置に依存する．すなわち，ある聴取者が音から得る空間的な知覚を，HRTF の合成に基づいて制御しようとするならば，理想的には，何らかの方法で，その HRTF をその個人に対し全方向にわたって得なければならない．特定のダミーヘッド や被験者を対象にして計測された HRTF のデータ
ベース，あるいはそれに準じたものはいくつか存在する6)–8) ものの，あらゆる個人に対して HRTF を計測するのは，被験者への身体的負担や計測を行うための設備・装置の制約などから現実的ではない．このため，個人ごとの HRTF を推定するための手法が必要とされる．有限方向の HRTF から任意方向の HRTF を補間する手法が提案されており9)–11) ，ある方向の HRTF を得るという意味では推定手法ということができるだろう．しかし，補間のもととなる HRTF は何らかの方法で得なければならない．福留12)や戸嶋ら 13) は，被験者の頭部を型取りしたダミーヘッド を作成して，HRTF の測定を行っている．戸嶋らはさらに，型取りしたダミーヘッド を用いて計測した HRTF と被験者本人の HRTF の実測結果を定量的に比較した上で，聴取実験による比較も行っている．この手法では，被験者への身体的負担の面では軽減が期待できるものの，計測を行う空間 (無響室など ) は必要である．
上記の音響測定に基づく頭部伝達関数の推定手法の問題点を解決する可能性のある手法として，頭部形状を 3 次元計測し ，それをデータとして境界要素法 (BoundaryElement Method，以下 BEM と書く) などを用いて HRTF を推定するものがある．頭部形状の 3 次元計測は比較的短時間で終了するため，被験者にかける身体的負担も少なく，数値計算によって HRTF を推定するため，測定で用いる特別な空間も必要としない．このことから，本手法は，他の手法よりも簡便に HRTF を得られるという点で有望と考えられる．この手法に関する研究としては，古くは大野木らによる頭部まわりの音場の数値計算14) ，塩谷15) や吉田ら 16) による耳介の伝達関数に関する検討がある．また，近年では，コンピュータの演算能力の向上に伴い，HRTF をより高い周波数帯域まで計算する研究も行われている．代表的なものとして，Kahana et al. によるもの17) や，Katz による検討18), 19) などがある．Otani et al. は，任意の方向の HRTF を高速に計算する方法を提案している20) ．しかし，これらの検討では，頭部形状の 3 次
元計測の精度との関連は明らかにされておらず，推定結果と実測結果との精度の比較が十分にされているとはいえない．一方，Tao et al. は，頭部形状の簡単化と，計算誤差の影響を考察している22) ．Terai et al. は，計算に用いた頭部形状と同様のモデルをプラスチックで作成し，実測結果と推定結果の比較を行い，よい一致がみられることを報告している21) ．これらの研究では，実際の被験者の HRTF を比較対象としていないため，HRTF としての実用性を明らかにするまでには至っていない．
以上の研究を踏まえ，本報告では，3 次元計測で得られたダミーヘッド と被験者の頭部形状をもとに，BEM による HRTF の推定を行うこととする．推定された HRTFと，実際の被験者を対象とした実測結果との比較を行い，その差について考察するとともに，頭部形状の違いによって生じると考えられる HRTF の差を，数値的に推定できるかど うかについて検討を行う．また，本稿では，計算手法の概略を説明し，BEMを用いた HRTF の推定において注意すべき点を述べる．

2

BEM を用いた HRTF の推定

2.1

境界積分方程式

本稿で解析の対象とする聴取者の頭部周りの音場は，図 1 に示されるような系に抽象化できる．音場を計算する領域を Ω とし，境界を Γ とする．自由空間を対象とするため，Ω は無限の領域をもっている．Γ 上の点 r p における法線ベクトル np の向きは，Ω からみて外向きになるように設定する．Ω 内に音源がないと仮定したとき，次に示
す Kirchhoff-Helmholtz の境界積分方程式が成立する．

ΩΓ

Γ3
Γ2
Γ1nqnpΓq rq
ΓprPq =|rP − rq |開空間音場の概念図(ZZ
Γ(r q )

∂P (r q , ω)∂G(r p , r q , ω)G(r p , r q , ω)
− P (r q , ω)∂nq∂nq(ZZP (r P , ω) =rp

rP
図1

C(r p )P (r p , ω) =ΓMΓ(r q )G(r P , r q , ω)∂P (r q , ω)∂G(r P , r q , ω)
− P (r q , ω)
∂nq
∂nq)

dΓ,(1))dΓ,(2)となる．ただし，r p ∈ Γ，r P ∈ Ω である．すなわち，式 (1) は境界上の点について成立する式で，表面境界積分方程式と呼ばれる．C(r p ) は，r p における境界形状の滑らかさによって決まる定数で，r p から Ω を見込んだ立体角を 4π で割った値となる．r pが滑らかであれば，1/2 となる．式 (2) は，Ω 内の点について成立するもので，内部境界積分方程式と呼ばれる．両式中の G は Helmholtz の Green 関数であり，3 次元問題
では，
e−jk|r p −r q |G(r p , r q , ω) =(3)4π|r p − r q |となる．
境界 Γ を，解析する音の波長に比べて十分に小さい M 個の要素 Γi (i = 1, · · · , M )に離散化して，境界条件を適用した上で式 (1) を要素を構成する節点に対して計算して連立させると，境界 Γ 上の節点の音圧あるいは音圧傾度が求められる．それを式 (2)に代入して，Ω 内の音圧を求める．節点数を N とすると，係数行列の組立てや連立方
程式の計算には，特別な方法を用いなければ O(N 3 ) の計算量が必要となる．すなわ
ち，節点数が 10 倍となれば，計算にかかる時間は数千倍のオーダーとなる．
HRTF を計算するためには，音源が Ω 内になければならないが，相反則を適用することにより，上記のような系を考慮すればよいことがわかる．このことについては，

後の 2.3 節で述べる．BEM を用いた音場解析の詳細については，文献23), 24) を参考にされたい．

2.2

解の一意性の問題

無限領域をもつ問題を BEM で解析するときには，解の一意性の問題 (Uniquenessproblem) に注意しなければならない．具体的には，無限領域問題における境界積分方
程式は，図 1 では塗りつぶされている境界 Γ の内側を領域としてもつ閉領域問題のもつ固有周波数で一意な解をもたない，というものである．これにより，無限領域をもつ問題では，境界上の音圧を求めるときに，式 (1) をそのまま用いることはできない．
これは境界積分方程式がもつ数学的な問題であり，物理的な問題ではない．
この問題の解決法は，大きく分けて 2 つある．1 つは，A. J. Burton と G. F. Miller が提案した方法25) である．これは，式 (1) の両辺を点 r p の法線方向 np で偏微分した式，

1 ∂P (r p , ω)=2∂np(ZZΓ(r q )

∂G(r p , r q , ω) ∂P (r q , ω)
∂ 2 G(r p , r q , ω)
− P (r q , ω)∂np
∂nq
∂np ∂nq)dΓ,(4)と，式 (1) 自身を適切な係数 κ で線形結合した式を用いれば，解の一意性が保証されるというものである．ただし ，式 (4) では，r p 上で境界が十分に滑らかであることを仮定している．この式は，式 (1) を法線方向に微分したものなので，NDF (NormalDerivative Form) と呼ばれる．これに対し，式 (1) は BF (Basic Form) と呼ばれる．もう1 つは，H. A. Schenck が提案したもの26) で，境界積分方程式と，境界内の幾つかの点における内部積分方程式を連立させるという方法であり，CHIEF(Combined Helmholtz
Integral Equation Formulation) と呼ばれている．Burton らの方法は，数学的に洗練されているが，境界積分方程式をさらに微分するため，特異性の高い関数の数値積分が必要となるのが難点である．その一方で，Schenck の方法は，数値的にはインプリメントしやすく，SYSNOISE など ，商用の BEM 解析ソフトでも採用されているが，境界内に置く付加点の位置や数について注意が必要であり，場合によっては，付加点を設けても解が発散してしまう場合がある．これを避けるために，Wu et al. は，境界内部に設けた点ではなく，小さな「領域」において満たされる方程式を境界積分方程式に連立させることを提案している27) ．
本稿では，上記の得失から，Burton らの方法をインプリメントした．特異積分の数値計算手法に関しては，様々な検討がなされているが，インプリメントの容易さから，
Kirkup のもの28) を適用した．

2.3

相反則の適用音場における相反則は，A 点に置いた音源によって B 点に生成される音圧は，B 点に音源を置いたときに A 点に生成される音圧に等しいというものである．これを，HRTFの計算に適用すると便利である．
HRTF を計算する場合，定義にしたがえば，音源を領域 Ω 内の一点に置き，両耳の音圧を計算することになる．音源の位置が変われば，境界上の音圧や音圧傾度は変化するため，音源の位置が変わるたびに境界上の音圧や音圧傾度を解き直す必要がある．
しかし，ここで相反則を適用すれば，外耳道入口に位置に音源を置き，受音点を Ω 内の点としても，計算される音圧は同じはずである．この場合，境界面が振動しているとみなせるので，受音点を変えるときに境界上の音圧や音圧傾度を計算し直す必要はない．Terai et al. は，相反則を適用して HRTF を計算し，実測とよい一致がみられることを確認している21) ．
本稿でも，相反則を適用して HRTF を計算した．つまり，音源を外耳道入口に置いて，様々な受音点において計算された音圧を，頭部中心相当位置に点音源がつくる音圧で割ったものを HRTF の推定値とした．

2.4

呼吸球の放射音場の計算作成したプログラムの動作を確認するため，理論解のわかっている呼吸球の放射音場の計算を行った．一定速度 U で振動する半径 a の呼吸球が，球の中心より距離 r の位置につくる音圧 P (r, ω) は次式のようになる．

ajka −jk(r−a)
P (r, ω) = U ρcer 1 + jka(5)ただし，ρ, c はそれぞれ空気の密度と空気中の音速で，k = cω は波数を表す．これにしたがい，半径 0.1 m の呼吸球の表面上の音圧を計算した結果を図 2 に示す．横軸は周波数で，縦軸は相対レベルである．実線は理論解を示しており，×は式 (1) のみを用いて計算した結果，○は 2.2 節で述べた Kirkup の方法で計算した結果である．直径0.2 m を一波長とする音の周波数である約 1700 Hz 付近で，×は理論値から大きくはずれているのに対し，○は理論値とよく一致している．このことから，作成したプログラムで，解の一意性の問題に対処できることが確認された．

3HRTF の推定結果と実測結果の比較

3.1

頭部形状の 3 次元計測

HRTF 推定のためのデータとなる頭部形状の 3 次元計測について述べる．頭部形状は，ミノルタ製の光学式 3 次元計測装置 VIVID300 を用いて行った．測定系の概略を
図 3 に示す．測定対象は，被験者 2 名とダミーヘッド (高研 SAMRAI) である．光学式の 3 次元計測装置は，レーザー光が当たらない，または反射しない部分の形状は計測できないため，計測の際には頭部を回転し，水平方向 45◦ おきに 8 回，頭頂部を 1 回の計 9 回の測定を行った．一方向の測定にかかる時間は 0.6 s であり，頭部の回転を考Sound Pressure Level at Boundary Surface

Relative Level [dB]

150

145

Laser 3D
scanner SCSI
140

Minolta
VIVID300

PC

Subject

135
AnalyticalBF
BF+NDF
130
0500

1000
Frequency [Hz]

1500

2000

図 3 3 次元スキャナによる頭部形状の計測

図 2 呼吸球の表面上の音圧の計算値と理論値の比較

慮しても全体の測定は 30 分前後で終了する．各方向から測定された 3 次元形状を，あらかじめ測定対象に付けた目印 (シール ) をもとに貼り合わせた．
被験者に関しては，頭髪がレーザー光を反射しないため，ゴム製のスイミングキャップをかぶって計測を行った．それでも計測ができない部分については，貼り合わせ後に適切な大きさの要素で穴埋めした．特に，外耳道内部にはレーザー光が届かないため，外耳道入口はふさいだ．このような手順で作成したデータを，要素モデルと呼ぶこととする．測定結果は，一辺が約 1 mm の三角形の要素で構成され，要素数は約 200000となる．要素数 200000 では計算コストが莫大となるので，これを適宜平滑化して要素数を減らした．例として，要素数を減らした後のダミーヘッド と被験者 1 の要素モデルを図 4 に示す．

(a) ダミーヘッド(b) 被験者 1
図 4 要素モデル (シェーデ ィング表示)

3.2

HRTF 推定の条件

ダミーヘッド については，要素数を 9610 まで減らした要素モデルを用いて計算を行った．被験者 1 については，Terai et al. の検討結果21) を参考に，要素数を 3012，5142，7626 に削減した要素モデルを作成し，計算する周波数によって異なる要素数のモデルを用いた．
相反則より，外耳道入口が振動しているものとし ，水平面上に頭部中心から 1.5 m離れた受音点を 0◦ から 357.5◦ まで 2.5◦ おきに設置した．方位角は，正面前方を 0◦ ，正面後方を 180◦ とし ，被験者の右方向を 90◦ ，左方向を 270◦ とした．本稿の計算結果は，左耳の外耳道入口が振動しているものとするため，90◦ 方向は音源が耳と反対側にあり，270◦ は音源と耳が同じ方向にあることを意味する．
境界条件は，頭部を含め全ての面を音響的に剛とした．これは，音圧傾度を 0 に設定することに対応し，音響インピーダンスを無限大に設定することを意味する．ただし，首の底面は，測定範囲の関係で便宜的に生じる境界面なので，吸音の境界条件として，平面波の音響インピーダンス ρc を設定した．要素は一定要素とした．この場合，要素数と節点数は一致する．

3.3

HRTF 測定の条件

測定は，東北大学電気通信研究所先端音情報研究分野にある球状スピーカアレ イを用いて行った．推定時には，外耳道入口はふさいだものとしているため，測定時もダミーヘッド または被験者の耳はふさいだ状態とした．外耳道入口をふさぐことによって，HRTF の推定/測定結果には，外耳道による共振の影響が表れないと考えられるが，適切な補正を行えば，このような状態で得られる HRTF は有用であることが知られている29) ．また，推定時には，頭部形状のみをデータとして用いているため，胴体などの影響は表れない．実測結果をこれに対応させるため，ダミーヘッドについては，頭部のみを用いて測定を行った．被験者については，吸音材料を胴体にあて，胴体からの反射を極力抑えた条件で測定を行った．測定は，方位角 5◦ ごとに行った．

3.4

推定結果と実測結果の比較

ダミーヘッドについては 9500 Hz まで，被験者 1 については 4000 Hz までの周波数帯域で HRTF の振幅特性を方位角にわたって比較したものを，それぞれ図 5，図 6 に示す．それぞれの図で，(a) は実測結果，(b) は推定結果を示している．
図 5(a) をみると，方位角が 0◦ (360◦ ) から 200◦ 付近まで変化するにしたがって，周波数特性のデ ィップの位置が 2000 Hz 付近から徐々に高くなっている様子がうかがえる．この変化は，推定結果にも表れており，ダミーヘッド の HRTF のもつ特徴を推定できていることがわかる．方位角が 100◦ から 200◦ 付近，つまり音源が後方にある場合の，デ ィップやピークの周波数軸上の位置の方位角による変化にも，実測結果と推定結果の間に傾向の一致がみられる．図 6(a) では，方位角が 150◦ 付近で，低周波数帯域の盛り上がりがみられる．また，方位角 250◦ 付近で，4000 Hz 付近にゆるやかなディップがみられる．これらの傾向は，図 6(b) でも同様にみられる．被験者 1 の HRTFのもつ大まかな特徴については，推定結果と実測結果の間には一致がみられるといえる．これらの特徴は，ダミーヘッド と被験者 1 で異なるものであり，推定結果はその違いを大まかに推定しているといえる．本稿で用いたプログラムを用いることにより，個人の HRTF を推定できる可能性が示唆された．
その一方で，音源の方位角が 90◦ 付近にあるとき，つまり音源と耳が逆方向にあるときに，図 5 でも図 6 でも，HRTF の実測結果と推定結果の間には大きな差がみられる．特に，推定結果には，実測結果にみられない周波数特性のデ ィップが 4 kHz 付近にみられ，そのデ ィップの位置は方位角の減少につれて高い周波数にシフトしているようにみえる．
これらのことから，HRTF の実測結果と推定結果には大まかな傾向の一致がみられるものの，方位角によっては大きな差がみられるところがある．この原因としては，計算時に与える境界条件や 3 次元形状の微妙な計測誤差など ，様々なものが考えられる．しかしながら，人体表面の音響特性を，音響インピーダンスで規定するのは極めて難しく，そこに個人差を含めることは現実的ではないと思われる．現実性を考慮しながら，実測結果と推定結果にある差の原因について検討を深めていきたい．
本稿における計算は，情報シナジーセンターの並列計算機 (par) のバッチジョブとして行った．ジョブのクラスは p16 である．ダミーヘッド のモデルについて，HRTFの 1 Hz あたりの成分を計算する時間は約 100 分であった．

4 まとめと今後の課題
本稿では，3 次元計測と境界要素法による HRTF の推定について述べた．HRTF の推定結果を実測結果との比較を行い，実測結果のもつ大局的な傾向を推定できていることが明らかとなった．また，頭部形状の違いによって HRTF に生じる差も，大局的ながら推定されており，この方法で個人の HRTF を推定できる可能性が示された．
その一方で，実測結果と推定結果が一致しない部分もあった．推定条件は，実測の条件を正確に反映しているとは言えない．この差をなくすことができれば，本稿の手法の有効性はより高くなるといえる．そのためには，境界条件，実測時の胴体部分の影響など ，様々な要因について検討を深める必要がある．
本稿では並列計算機での計算結果を示したが，HRTF をより高い周波数まで高速に計算するには，スーパーコンピュータでのインプ リメンテーションが不可欠である．
現在のところ，ベクトル化において十分な性能を出すには至っていないが，2.2 節で述べた手法をベクトル化に適したアルゴ リズムに変換することが困難である．早急にベクトル化を実現し，境界要素法によって音場を高速に計算できるプログラムとして，ユーザの方々の利用に供することができるものとしたい．

Measured HRTF (SAMRAI, L channel)

Level [dB]200-20-40

cyenquFre-6002000
4000
6000
8000
10000 0z]
[H110 60210 160
260310Azimuth [Degree]

10(a) 実測結果

Estimated HRTF (SAMRAI, L channel)

Level [dB]200-20-4002000
cyen
qu
Fre

40006000z][H図58000
10000 0110 60210 160260310
Azimuth [Degree]

10(b) 推定結果
ダミーヘッド の HRTF の推定結果と実測結果 (振幅特性，〜9500 Hz)

Level [dB]

Measured HRTF (subject1, L channel)

15100-10-200cyen
quFre

10002000z][H30004000 060160 110
210
310 260
Azimuth [Degree]

10(a) 実測結果

Level [dB]

Estimated HRTF (subject1, L channel)
15100-10-20-30-400cyen
quFre

10002000z][H3000図64000 060160 110
210
310 260
Azimuth [Degree]

10(b) 推定結果
被験者 1 の HRTF の推定結果と実測結果 (振幅特性，〜4000 Hz)

謝辞 本稿における頭部伝達関数の推定で用いたプログラムの開発には，情報シナジーセンターの研究開発公募の助成を受けた．また，頭部伝達関数の測定は，東北大学電気通信研究所共同プロジェクト研究 (H14/A06) のもとに，先端音情報研究分野の球状スピーカアレ イを用いて行われた．以上の助成に対し，感謝を申し上げる．

参考文献
1) J. Blauert: Spatial hearing (MIT Press, Cambridge, 1983).
2) D. R. Begault: 3-D SOUND for virtual reality and multimedia (AP Professional, Cambridge, 1994).
3) たとえば，F. Ramsey: Spatial audio (Focal Press, Oxford, 2001).
4) S. Takane, Y. Suzuki, T. Miyajima and T. Sone, “A new theory for high definition virtual acoustic display named ADVISE,” Acoust. Sci. & Tech., 24(5), 276-283(2003).
5) S. Takane, S. Takahashi, Y. Suzuki and T. Miyajima, “Elementary real-time implementation of a virtual acoustic display
based on ADVISE,” Acoust. Sci. & Tech., 24(5), 304-310(2003).
6) B. Gardner and K. Martin, “HRTF measurements of a KEMAR dummy-head microphone,” MIT Media Lab PerceptualComputing - Technical Report #280, 1-7(1994), http://sound.media.mit.edu/KEMAR.html
7) S. Takane, D. Arai, T. Miyajima, K. Watanabe, Y. Suzuki and T. Sone, “A database of Head-Related
Transfer Functions in whole directions on upper hemisphere,” Acoust. Sci. & Tech., 23(3), 160-162(2002),
http://www.ais.riec.tohoku.ac.jp/lab/db-hrtf/8) http://www.itakura.nuee.nagoya-u.ac.jp/HRTF/index-j.html
9) 西野隆典，梶田将司，武田一哉，板倉文忠，“水平面上の頭部伝達関数の補間，” 音響学会誌，55(2)，91-99(1999)．
10) 西野隆典，梶田将司，武田一哉，板倉文忠，“水平方向及び仰角方向に関する頭部伝達関数の補間，” 音響学会誌，
57(11)，685-692(2001)．
11) K. Watanabe, S. Takane and Y. Suzuki, “Interpolation of Head-Related Transfer Functions based on the CommonAcoustical-Pole and Residue model,” Acoust. Sci. & Tech., 24(5), 335-337(2003).
12) 福留公利，“3 次元音場再生のための頭部伝達関数の補間，” 音講論，583-586(1999.9-10)．
13) 戸嶋巌樹，青木茂明，平原達也，“頭部を型取りしたダミーヘッド の頭部伝達関数，” 音講論，563-564(2004.3)．
14) 大野木重夫，加川幸雄，山淵龍夫，財田一也，“境界要素法による頭部まわりの音場シミュレーション，” 音講論，731-732(1989.10)．
15) 塩谷達，“有限要素法を用いた頭部伝達関数の推定に関する基礎的研究，” 東北大学工学部卒業論文 (1989)．
16) 吉田浩子，芹川光彦，佐藤克昌，“有限要素法による実耳モデルを用いた頭部伝達関数の計算，” 音講論集，453-454(1989.10)．
17) Y. Kahana, P. A. Nelson, M. Petyt and S. Choi, “Boundary element simulation of HRTFs and sound fields produced byvirtual acoustic imaging systems,” Proc. 105th AES Convention, No. 4817(1998).
18) B. F. G. Katz, “Boundary element method calculation of individual head-related transfer function. I. Rigid model calculation,” J. Acoust. Soc. Am., 110(5), Pt. 1, 2440-2448(2001).
19) B. F. G. Katz, “Boundary element method calculation of individual head-related transfer function. II. Impedance effects
and comparisons to real measurements,” J. Acoust. Soc. Am., 110(5), Pt. 1, 2449-2455(2001).
20) M. Otani and S. Ise, “A fast calculation method of the head-related transfer functions for multiple source points on theboundary element method,” Acoust. Sci. Tech., 14(5), 259-266(2003).
21) K. Terai and I. Kakuhari, “HRTF calculation with less influence from 3-D modeling error: Making a physical humanhead model from geometric 3-D data,” Acoust. Sci. & Tech., 14(5), 333-334(2003).
22) Y. Tao, A. I. Tew and S. J. Porter, “A study on head-shape simplification using spherical harmonics for HRTF computation at low frequencies,” J. Audio Eng. Soc., 51(9), 799-805(2003).
23) 日本建築学会 編：室内音場予測手法 —理論と応用— (日本建築学会，2001)．
24) Boundary element methods in acoustics, edited by R. D. Ciskowski and C. A. Brebbia (Elsevier Applied Science, London, 1991).
25) A. J. Burton and G. F. Miller, “The application of integral equation methods to the numerical solution of some exterior
boundary-value problems,” Proc. Roy. Soc. Lond., A. 323, 201-210(1971).
26) H. A. Schenck, “Improved integral formulation for acoustic radiation problems,” J. Acoust. Soc. Am., 44(1), 41-58(1968).
27) T. W. Wu and A. F. Seybert, “A weighted residual formulation for the CHIEF method,” J. Acoust. Soc. Am., 90(3),1608-1614(1991).
28) S. Kirkup: The boundary element method in acoustics (Integrated sound software, London, 1998).
29) H. Møller, “Fundamentals of binaural technology,” Applied Acoustics, 36, 171-218(1992).