主成分分析に基づく個人性頭部伝達関数の推定
頭部伝達関数（HRTF：Head-Related Transfer Function）を用いて立体的な音響を提示するバイノーラル再生においては，聴取者本人のHRTFを忠実に再現することが良好な再生をするために重要である。しかし，聴取者本人のHRTFの測定には多大な労力とコストを要するので，主成分分析による基底展開と，部分的な測定からHRTFを合成する簡便な推定手法を検討した。また，HRTFの所要近似精度を主観および客観評価実験を行って求め，提案手法の有効性を検証した。
1. はじめに
頭部伝達関数（HRTF）とは，音波が頭部や身体，耳介などでの反射や回折を経て聴取者の外耳道入口に到来するまでの伝達特性を記述した関数である。ヘッドホンを用いた立体音響再生技術として知られるバイノーラル再生方式*1 は，音響信号にHRTFの時間領域表現である頭部インパルス応答（ HRIR ： Head-Related Impulse Response）と再生系逆フィルター*2 のインパルス応答を畳み込むことによって実現でき，音を使って情報を提示する聴覚ディスプレーでの音像定位などに応用されている。
バイノーラル再生方式において良好な再生をするためには，聴取者本人のHRTF（個人性HRTF）を精度よく近似する必要がある。しかし，HRTFの測定には多大な労力とコストを要するので，線形近似，スプライン近似*3 など，さまざまな補間手法が提案されてきた1)2)3)4)。それらのモデルも時間領域表現，周波数領域表現，極・零モデル*4 など多様にわたり5)6)7)8)，それらの比較検討も行われている3)4)9)。
本稿では，HRTFを簡便に導出する手法の一検討として，少ない測定方向から全方向のHRTFを推定する手法を提案する。提案手法では測定したHRTFを主成分分析してHRTFの振幅周波数特性（以下，HRTF振幅応答）を基底展開*5 するとともに，推定した重み係数を用いて未測定方向のHRTF振幅応答を合成する。また，HRTF推定時の所要近似精度を音源別，方向別に主観評価実験を行って調べ，提案手法の有効性を検証した結果を示す。
2. HRTFの測定
HRTFの測定は当所の音響無響室で行った（1図)。1図に示すように，被験者は頭部をヘッドレストに支持された状態で回転いすに腰掛ける。スピーカーは被験者の頭部を同心円の中心として，円弧状に仰角10°間隔で配置されている。円弧の半径は1.3mである。被験者の腰掛ける回転いすが一定の角度おきに回転し，全方向の測定を行う。
測定信号には，時間引き延ばしパルス*6 の一種であるOATSP（Optimized Aoshima’s Time Stretched Pulse）を用いた。サンプリング周波数は48kHz，信号長は16,384サンプルである。スピーカーごとに再生され，外耳道入口で収音された測定信号は，同期加算による平均化，ローパスフィルターおよびOATSP，測定系それぞれの逆フィルターのインパルス応答との畳み込みを経て，512サンプルのインパルス応答として矩形（くけい）窓を用いて切り出される。512サンプルの応答長（約10ms）はHRIRが収束するのに十分な長さとして一般的に用いられる値である。このインパルス応答を離散フーリエ変換により伝達関数に変換する。このとき，伝達関数の半分は鏡像関係の値となるので，分析には256サンプルのHRTFを用いれば十分である。
被験者は10歳代～50歳代の男性75名，女性88名の合計163名である。測定方向は方位角5°間隔の72方位，仰角－20°～90°まで（一部の被験者では0°～80°まで）の計793（648）方向である。
3. 個人性HRTFの推定
3.1 主成分分析
主成分分析とは，変量間の相関関係を解析して直交基底を求め，変量全体の変動を簡潔に表す合成変量を導出する多変量解析手法である。それぞれの合成変量に集約される変動の程度が寄与率として得られ，情報の次元を圧縮するための手法として用いることができる。
i番目の被験者のj番目方向のHRTFの振幅応答を要素とするベクトルをhijとすると，変量ベクトル
は(1)式および(2)式で表される。なお，HRTF振幅応答ベクトルhijの要素はあらかじめ対数変換されているとする。
ここで，Ndは測定方向の数を，Tは転置行列を表す。また，
は平均HRTF振幅応答ベクトルを表し，分析に用いた変量ベクトルの数Nsを用いて(3)式で計算される。
変量ベクトルの長さとなる変量数NvはHRTF振幅応答のベクトル長Nƒと測定方向の数Ndの積である。変量ベクトルのすべての組み合わせに対して共分散*7 を計算し，分散・共分散行列*8 を作成する。主成分は固有値問題を解くことで得られ，分散・共分散行列の固有ベクトルが主成分に，固有値が寄与率となる。
導出された主成分は被験者の集合のHRTF振幅応答の直交基底であり，これを列に持つ主成分行列Cを用いてi番目の被験者のHRTF振幅応答hiを(4)式で再合成することができる。
ここで，ωiは各主成分に対する重み係数ベクトルである。
3.2 HRTFの推定
導出された平均HRTF振幅応答および主成分をHRTF振幅応答の普遍的な基底であると仮定すると，任意のHRTF振幅応答gkは各主成分の重み係数との積和に平均HRTF振幅応答を加えることで合成できる。
ここで，ωkは各主成分に対応する重み係数ベクトルである。
以下，部分的に測定したHRTF振幅応答から全方向のHRTF振幅応答を推定することを考える。測定したn方向のうちl番目の方向をdlとし，部分変量ベクトル
を(6)式で定義する。
また，主成分行列Cを(7)式で示すようにNdの部分行列に分割し，部分主成分行列Cdを(8)式で定義する。
このとき，重み係数ベクトルωkは近似的に(9)式で計算できる。
ここで，記号 \ は行列の左除算*9 を表す。(9)式は劣決定系*10 であり，(9)式の最小２乗解として求めた重み係数ベクトルωkを再度(5)式に代入して，全方向のHRTF振幅応答を合成する。
次元圧縮は重み係数ベクトルωkを(5)式に代入するときに，主成分数を制限することで行う。すわなち，主成分行列Cの列のうち，寄与率（固有値）の小さいものを削除することによって行う。なお，次元圧縮と近似精度の関係については４章で検討する。
3.3 推定例
水平面内の72方向のHRTF振幅応答を分析対象とした平均HRTF振幅応答を2図に，第１主成分～第３主成分を3図4図5図に示す。ただし，方位は正面を0°として反時計回りに定めた（左側方が90°)。また，SN比を考慮し，各方向の256サンプルのHRTF振幅応答のうち，17kHz以下に相当する182サンプルを分析対象とした。従って，変量数Nvは13,104（＝182×72）である。分析には，163名のうちの80名分のHRTF振幅応答を用いた。右耳の応答は方位角を反転し，被験者の本来の左耳のデータとは別の左耳のデータとして扱った。すなわち，１人の被験者から２つの変量ベクトルを得た。変数の一覧を1表に示す。なお，分析前に，それぞれのHRTF振幅応答をピーク値が－２dBFS（dB Full Scale）となるよう正規化した。6図にHRTF振幅応答の測定値を，7図と8図に72方向と４方向だけの測定値を用いて推定したHRTF振幅応答の例を示す。7図と8図はほぼHRTF振幅応答が推定できていることを示しているが，以下，本手法の推定精度を検討する。
| Ns | 160 |
|---|---|
| Nd | 72 |
| Nf | 182 |
| Nv (＝Nf×Nd) | 13,104 |
| n | 72または4 |
4. 測定方向の検証
HRTFを簡便に推定するという研究の目的からは，測定方向は少ない方が良い。そこで，測定方向の数と測定方向を変えてHRTF振幅応答を推定した。推定したHRTF振幅応答のスペクトル歪（ひず）み（SD：Spectral Distortion)*11 の値SDを(10)式で定義し，この値に基づいて推定に必要な測定方向数を検証した。
ここで，g(i)は測定したHRTF振幅応答の各要素を，g'(i)は推定したHRTF振幅応答の各要素を表す。３章で述べたようにNƒ=182である。評価対象は３章の分析に用いた80名以外の20名の水平面内の72方向のHRTF振幅応答とした。
測定方向数，主成分数ごとの平均SDを9図に示す。測定方向は正面を起点として等間隔である。また，変量ベクトルNsの数は160なので，全主成分の数は１少ない159である。全主成分を用いた場合には，SDは測定方向数36で3.0dB以下，測定方向数９で3.7dB以下の精度でHRTF振幅応答を近似できた。主成分の数を削減すると，測定方向数が多い場合の近似精度が低くなった。しかし，測定方向数が少ない場合の歪みの拡大は緩やかで，測定方向数８～９を境として，全主成分を用いた場合より良くなった。この結果は，寄与率の低い下位の主成分は任意のHRTF振幅応答の基底になるという仮説が成り立たないことを示している。
10図に分析に用いた80名のHRTF振幅応答を推定した結果を示す。9図と比較して，主成分数が多いほどSDが改善されていることがわかる。これは，抽出した主成分が80名の被験者のHRTF振幅応答の直交基底であることによる。従って，9図は提案手法の精度限界を示していると言える。
主成分数を10，測定方向数を２に固定し，測定方向の組み合わせを変えて平均SDを計算した。結果を11図に示す。SDは測定方向（5°，115°）で最小（4.22dB）となり，この測定方向の近傍ではこの方向から離れるに従って大きくなることがわかる。
測定方向数を３とした場合の平均SDの小さい組み合わせの上位1,000組を12図に示す。平均SDの非常に小さい組み合わせ（１位～10位）は２つの領域に，11位～1,000位は４つの領域に集まることがわかる。４つの領域において，SDが極小となる測定方向は（95°，155°，350°)，（5°，65°，130°)，（10°，115°，195°)，（100°，205°，340°）であった。12図には示していないが，1,000位以降の平均SDは４つの領域を拡大する形で広がる。同様に，測定方向数を４とした場合の平均SDの小さい組み合わせの上位10,000組を計算した。ただし，３方向での結果を考慮し，測定方向の１つをほぼ正面方向（340°‐20°）に限定した。平均SDが小さい測定方向の組み合わせは（0°，80°，135°，210°）で極小となる領域と，（20°，120°，210°，310°)，（20°，95°，155°，325°）で極小となる領域とに分布する傾向があった。10,000位以降の平均SDはこれらの領域を拡大する形で広がり，その傾向は３方向の場合と同様である。
以上の結果から，少ない測定方向で良好な近似精度を得るための条件は
- 正面方向の近傍を測定すること
- 測定耳側の後方を測定すること
- 測定方向間の角度を広くすること
であると言える。
(分析に用いた被験者ではない場合)
(分析に用いた被験者の場合)
5. 所要近似精度の評価
提案手法の有効性を検証するために，HRTF推定時の所要近似精度を音源別，水平面内の方向別に主観評価実験を行って調べた。
評価にはMUSHRA（MUlti Stimulus test with Hidden Reference and Anchor）法を用いた。MUSHRA法は音質評価などに用いられる主観評価法の１つで，隠れ基準*12，隠れアンカー*13 を含む複数の評価音を同時に比較できるという特徴がある。
5.1 実験条件
基準音および評価音の合成には，被験者本人のHRTF，主成分分析により得られた方向別の平均HRTF振幅応答および第１主成分～第15主成分を用いた。主成分数は歪みの検知限を調べた予備実験の結果を基に，すべての主成分を用いて合成した評価音の近似精度が十分に良くなるように定めた。被験者を含む108名のHRTF振幅応答を用いて主成分分析を行った。音源信号に被験者本人のHRTFから求めたHRIRと再生ヘッドホンの逆フィルターのインパルス応答とを畳み込んだ信号を基準音とし，左右耳のいずれかを第n主成分までを用いて推定したHRIRに置き換えた信号を第n評価音とした。従って，次数nが大きいほど評価音の近似精度は向上する。また，n=０すなわち平均HRTF振幅応答より推定したHRIRに置き換えた信号を隠れアンカーとした。なお，HRTF位相応答は被験者本人のものを用いた。その他の実験条件を2表に示す。提示音圧レベルは正面方向の基準音が同一レベルとなるよう調整し，HRTFの方向の違いによるレベル差を維持した。被験者は20歳代～30歳代の正常な聴力を有する６名である。実験は左右耳別，方向別に行った。基準音と評価音を自由に聴き比べ，その差異を５から１の間で小数第１位まで評価した（3表)。
| 音源信号 | 白色雑音/楽音/スピーチ |
|---|---|
| 方位 | 水平面内で10°間隔，36方向 |
| 提示音レベル | 63dBA |
| 試行回数 | 左右耳×36方向×延べ18回 |
| 5 | 基準との差がわからない |
|---|---|
| 4 | わかるが気にならない |
| 3 | やや気になる |
| 2 | 気になる |
| 1 | 非常に気になる |
5.2 実験結果
全方向，全被験者の平均評価値を13図に示す。エラーバーは95％信頼区間である。13図から，いずれの音源においても，合成に用いる主成分数nの増加に伴って平均評価値が上昇することがわかる。また，白色雑音，楽音，スピーチの順で評価値が低いことが確認される。許容限を3.5とすると，白色雑音では第10主成分，楽音では第２主成分までを用いれば許容限以上の品質にすることができる。なお，スピーチの場合には平均HRIRを用いた場合においても許容限以上の品質であった。
主成分数を同じにして合成したHRTF振幅応答のSDを比較すると，SDの値は被験者によって異なっている。そこで，方向別のHRTF振幅応答のSDの値と平均評価値の関係をプロットした。結果を14図15図16図に示す。SDの計算は小数第２位までとし，隣接する±10°，±0.5dBの範囲で移動平均をとり平滑化した。右耳の結果は方位角を反転して合算した。14図15図16図においても，SDの増加に伴って評価値が下降する傾向が確認できる。しかし，方向によってその傾きは異なっている。そこで，方向，主成分数を要因とする２元配置分散分析*14 を行った。分散分析表を4表に示す。白色雑音と楽音においては主成分数と方向の要因とそれらの交互作用*15 で，スピーチにおいては主成分数と方向それぞれの要因で有意差があった。評価値の単調増傾向を考慮すると，主成分数の増加に伴い評価値が上昇する傾きが方向によって有意に異なると解釈できる。このことは方向によって所要近似精度が異なることを示唆している。
評価にあたっては，音像の定位感（方向感）のほかに，音色や音の立ち上がり，立ち下りなど，総合的に差異を判断するように指示した。被験者の内省報告では，特に，白色雑音において音色の差異とそれに伴う定位感の変化が多数の被験者において指摘された。また，前後方向と比較して，側方の弁別が困難であることを報告する被験者も複数いた。そこで，音像の定位感と音色についての追加実験を行った。実験条件は方向と試行回数以外は5.1節の実験条件と同じである。17図18図に平均評価値を示す。13図と比較して，主成分数が少ない評価音で平均評価値が大きい。この結果はHRTFの近似精度は音像の定位感だけでなく，音色などの他の要因にも影響があることを示唆している。
19図は14図15図16図の許容限3.5に相当するSDの値である。ただし，スピーチは常に許容限3.5以上の品質であり，3.5に相当するSDの値を求めることができないので，19図には白色雑音と楽音の結果だけを示した。なお，19図に示すSDの値以下では許容限以上の品質である。許容限3.5に相当するSDの値が方位角によって大きく異なる楽音では，許容限3.5に相当するSDの値は30°方向で最小になり，30°～180°方向でやや大きく，音源とは反対の方向（180°～360°）で更に大きくなっている。このことは，正面方向ではHRTFを精度よく推定する必要があること，音源の反対側では高い推定精度が必要ではないことを意味している。許容限3.5に相当するSDの値は音源や方向によって異なるので，HRTFの歪みを評価するための指標としてSDを用いる場合には，音源や方向の違いを考慮する必要があると言える。
| 白色雑音 | ||||
|---|---|---|---|---|
| 平方和 | 自由度 | 平均平方 | F | |
| 主成分数 | 16,799.0 | 16 | 1,049.9 | 1,033.68* |
| 方向 | 707.5 | 35 | 20.2 | 19.9* |
| 交互作用 | 840.2 | 560 | 1.5 | 1.48* |
| 楽音 | ||||
| 平方和 | 自由度 | 平均平方 | F | |
| 主成分数 | 1,827.0 | 9 | 203.0 | 176.18* |
| 方向 | 1,071.9 | 35 | 30.6 | 26.58* |
| 交互作用 | 450.1 | 315 | 1.4 | 1.24* |
| スピーチ | ||||
| 平方和 | 自由度 | 平均平方 | F | |
| 主成分数 | 211.0 | 9 | 23.4 | 38.29* |
| 方向 | 144.1 | 35 | 4.1 | 6.72* |
| 交互作用 | 113.9 | 315 | 0.4 | 0.59 |
*p<.05
6. おわりに
主成分分析に基づくHRTFの推定手法を提案した。また，主観評価実験を行って音源や方向の違いによって所要近似精度が異なることを示した。測定方向数が少ない場合には，正面方向と測定耳側の後方を測定し，測定方向間の角度を広くした組み合わせが適切であることがわかった。
開発した手法は多数の被験者のHRTF振幅応答を基に導出した主成分が任意のHRTF振幅応答に対する基底になると仮定し，少ない測定方向から全方向のHRTF振幅応答を推定しようとするものである。この仮定を満たすためには，主成分分析に基づいて構成した低次元空間が，HRTFの集合の統計的な特徴を十分に包含している必要がある。従って，今後，より多くの被験者のHRTFを測定し，多くのサンプル数を基に解析を行う必要があると言える。
本稿はAcoustical Science and Technology誌に掲載された以下の論文の内容を元に加筆・修正したものである。
K. Matsui and A. Ando：“Estimation of Individualized Head-Related Transfer Function Based on Principal Component Analysis，” Acoust. Sci. & Tech., Vol.30, No.5, pp.338-347（2009）