6343771 頭部伝達関数のモデリング装置、その方法及びそのプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6343771

(24)【登録日】2018年6月1日

(45)【発行日】2018年6月20日

(54)【発明の名称】頭部伝達関数のモデリング装置、その方法及びそのプログラム

(51)【国際特許分類】

   G10K 15/00 20060101AFI20180611BHJP

   H04S 7/00 20060101ALI20180611BHJP

【ＦＩ】

   G10K15/00 L

   H04S7/00 310

【請求項の数】6

【全頁数】23

(21)【出願番号】特願2014-161681(P2014-161681)

(22)【出願日】2014年8月7日

(65)【公開番号】特開2016-39493(P2016-39493A)

(43)【公開日】2016年3月22日

【審査請求日】2017年6月30日

(73)【特許権者】

【識別番号】000004352

【氏名又は名称】日本放送協会

(73)【特許権者】

【識別番号】899000079

【氏名又は名称】学校法人慶應義塾

(74)【代理人】

【識別番号】110001807

【氏名又は名称】特許業務法人磯野国際特許商標事務所

(72)【発明者】

【氏名】松井 健太郎

(72)【発明者】

【氏名】足立 修一

(72)【発明者】

【氏名】菅谷 真帆

【審査官】 岩田 淳

(56)【参考文献】

【文献】 特開２００７−１７４２２４（ＪＰ，Ａ）

【文献】 特開２００２−０９５０９６（ＪＰ，Ａ）

【文献】 特開２００８−２１１８３４（ＪＰ，Ａ）

【文献】 米国特許出願公開第２００６／０２７７０３４（ＵＳ，Ａ１）

【文献】 足立 修一，次世代プラント操業に期待される制御技術，計測と制御，日本，社団法人計測自動制御学会，２００８年１１月１０日，第47巻, 第11号，第915-920頁

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ１０Ｋ １５／００−１５／１２

Ｈ０４Ｓ １／００− ７／００

Ｈ０４Ｒ ５／００− ５／０４

(57)【特許請求の範囲】

【請求項1】

スピーカに印加された入力信号と前記スピーカから発せられた音声をマイクロフォンで測定して得られた出力信号とを入出力データとして用いて、漸近推定法により頭部伝達関数をモデリングする頭部伝達関数のモデリング装置であって、
前記入出力データを用いて、予め定められた高次のモデル次数を有した頭部伝達関数及び雑音モデルについての高次モデルを予測誤差法により推定する高次モデル推定手段と、
推定された高次モデルと、周波数領域における評価関数である対数尤度関数とを用いて最尤推定値を導出することで前記高次モデルを低次元化する低次元化手段とを備え、
前記低次元化手段は、
前記高次モデルの周波数伝達関数を求める周波数伝達関数算出手段と、
前記対数尤度関数を最小化することで前記高次のモデル次数よりも低い次数の低次モデルの推定値を求める低次モデル推定手段と、
前記低次の次数を更新して前記対数尤度関数を最小化させることを繰り返すことでそれぞれ推定された各頭部伝達関数の低次モデルと、高次参照モデルと、の間の音像定位知覚に係る特徴量の誤差をそれぞれ求め、前記特徴量の誤差が予め定められた許容条件を満たし且つ最低次数となるときの低次モデルを探索する低次モデル探索手段と、
を備えることを特徴とする頭部伝達関数のモデリング装置。

【請求項2】

前記音像定位知覚に係る特徴量は、前記頭部伝達関数の周波数特性上のピーク及びノッチであるスペクトラルキューの位置であることを特徴とする請求項１に記載の頭部伝達関数のモデリング装置。

【請求項3】

前記低次モデル探索手段は、
前記低次モデル推定手段でそれぞれ推定された各頭部伝達関数の低次モデルと、前記高次参照モデルと、の間のスペクトル歪をそれぞれ求め、前記スペクトル歪が予め定められた第１閾値以下、且つ前記スペクトラルキューの位置についての前記高次参照モデルとの間の誤差が予め定められた第２閾値以下の条件を満たし且つ最低次数となるときの低次モデルを探索することを特徴とする請求項２に記載の頭部伝達関数のモデリング装置。

【請求項4】

前記低次モデル探索手段は、
前記低次の次数を前記高次のモデル次数側から降順に更新して前記スペクトル歪が前記第１閾値より大きくなった場合に前記対数尤度関数の最小化処理を停止し、
前記スペクトル歪が前記第１閾値以下であったときの次数を起点として前記低次の次数を昇順に更新して前記スペクトラルキューの位置についての前記高次参照モデルとの間の誤差が前記第２閾値以下になったときの次数を前記最低次数として決定する、ことを特徴とする請求項３に記載の頭部伝達関数のモデリング装置。

【請求項5】

スピーカに印加された入力信号と前記スピーカから発せられた音声をマイクロフォンで測定して得られた出力信号とを入出力データとして用いて、漸近推定法により頭部伝達関数をモデリングする頭部伝達関数のモデリング方法であって、
前記入出力データを用いて、予め定められた高次のモデル次数を有した頭部伝達関数及び雑音モデルについての高次モデルを予測誤差法により推定する高次モデル推定ステップと、
推定された高次モデルと、周波数領域における評価関数である対数尤度関数とを用いて最尤推定値を導出することで前記高次モデルを低次元化する低次元化ステップと、を有し、
前記低次元化ステップは、
前記高次モデルの周波数伝達関数を求める周波数伝達関数算出ステップと、
前記対数尤度関数を最小化することで前記高次のモデル次数よりも低い次数の低次モデルの推定値を求める低次モデル推定ステップと、
前記低次の次数を更新して前記対数尤度関数を最小化させることを繰り返すことでそれぞれ推定された各頭部伝達関数の低次モデルと、高次参照モデルと、の間の音像定位知覚に係る特徴量の誤差をそれぞれ求め、前記特徴量の誤差が予め定められた許容条件を満たし且つ最低次数となるときの低次モデルを探索する低次モデル探索ステップと、
を有することを特徴とする頭部伝達関数のモデリング方法。

【請求項6】

コンピュータを、請求項１に記載の頭部伝達関数のモデリング装置として機能させるための頭部伝達関数のモデリングプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、頭部伝達関数のモデリングに関し、特に音像定位技術などに用いられる頭部伝達関数のモデリング装置、その方法及びそのプログラムに関する。

【背景技術】

【0002】

頭部伝達関数（ＨＲＴＦ：Head-Related Transfer Function）とは、音響伝達関数であって、具体的には、模擬として頭がない状態での頭部中心に相当する位置から、頭外音源位置を経て、両耳鼓膜位置もしくは外耳道入口までの音響伝達関数のことである。
頭部伝達関数のモデリングとは、頭部伝達関数を有理多項式又はその比として表すことである。頭部伝達関数およびそのモデルは、音響信号の音像を空間内に擬似的に配置させるような音像定位技術など応用は多岐に渡る。

【0003】

一般に、システムの伝達関数Ｇ（ｑ）は、システムのインパルス応答ｇ（ｋ）のｚ変換で定義され、次の式（１）のように表わすことができる。ここで、ｑはシフトオペレータである。

【0004】

【数1】

【0005】

式（１）は無限インパルス応答（ＩＩＲ：Infinite impulse response）モデルを表している。それに対し、一般的に、頭部伝達関数は、そのインパルス応答が十分に収束するような次数Ｍで打ち切ることで、有理多項式を用いて次の式（２）のように、有限インパルス応答（ＦＩＲ：finite impulse response）モデルとしてモデリングされる。

【0006】

【数2】

【0007】

頭部伝達関数のインパルス応答は、スピーカから入力信号を印加し、例えばダミーヘッドの両耳に内蔵したマイクロフォンにより収音を行なうことで測定される。例えば、サンプリング周波数４８ｋＨｚで測定された場合、十分な収束が得られる長さは、およそ５１２サンプルである。そのため、測定信号は、このサンプル数と同じか又はやや長い矩形窓により切り出され、頭部伝達関数は５１２次や１０４８次等の高次有限インパルス応答モデルとなる。

【0008】

頭部伝達関数には、例えば両耳間時間差やレベル差、周波数特性上のスペクトラルキューなど、音像定位知覚に係る多くの特徴量が含まれる。十分に高次の有限インパルス応答モデルによる頭部伝達関数のモデリングにおいて、これらの特徴量が十分に保存されていることは、音像定位実験により確認されている。頭部伝達関数の同定において、スペクトラルキューなどは、音像定位知覚の手がかりとなる。

【0009】

従来、頭部伝達関数のモデリング法が知られている（例えば、非特許文献１，２、特許文献１，２参照）。
非特許文献１には、頭部伝達関数のモデリング法として、極零（pole/zero）モデルにより有理多項式の比で頭部伝達関数をモデリングする方法が開示されている。非特許文献１には、極零モデルが次の式（３）で定義されることが記載されている。

【0010】

【数3】

【0011】

式（３）において、Ｃは定数、ｐ_kは極、ｚ_kは零点であり、このうち極は共振に対応し、零点は時間遅れや反共振に対応する。

【0012】

前記した式（２）で表される有限インパルス応答モデルでは零点のみを用いてインパルス応答を表現している。これに対して、例えば式（３）に示す極零モデルのように零点だけではなく極を用いてモデリングすることで、より少ないパラメータ数で頭部伝達関数を表現することができる。非特許文献１に記載の技術では、極零モデルのパラメータを導出する際に、一次モデルとして有限インパルス応答モデルを求め、それぞれのインパルス応答の誤差、すなわち出力誤差の二乗和を最小化することで、パラメータを導出する。

【0013】

非特許文献２に記載された技術では、さらに、共振を方位によらないものとして、モデルの極を各方位の頭部伝達関数間で共通化した、共通極零（ＣＡＰＺ：Common-acoustical-pole and zero）モデルによるモデリングを行っている。非特許文献２には、音源受音点位置の位置ベクトルをｒ_jとしたとき、共通極零モデルは次の式（４）で定義されることが記載されている。

【0014】

【数4】

【0015】

式（４）において、極ｐ_kが方向に依存せず、零点ｚ_kのみが方向に応じて変化するため、このモデルでは、音源や受音点の変化に応じて変えるパラメータの個数を削減することができる。非特許文献２に記載の技術では、共通極零モデルを導出する際に、複数方向の頭部伝達関数の有限インパルス応答モデルを一次モデルとして求め、その近似として共通極零モデルを導出する。非特許文献２には、このとき、評価関数は、各方位の有限インパルス応答モデルと共通極零モデルの間のインパルス応答の誤差、すなわち出力誤差ε(ｒ_j，ｋ)を用いて次の式（５）のように表わされることが記載されている。

【0016】

【数5】

【0017】

式（５）において、Ｒは方向数、Ｍは有限インパルス応答モデルの次数である。この式（５）におけるＪを最小化することにより、共通極零モデルのパラメータは導出される。

【0018】

また、特許文献１に記載された技術は、複数の頭部伝達関数の有限インパルス応答モデルから主成分分析を用いて抽出した基本ベクトルを、バランスモデル近似技術により極零モデルとして模擬する方法に関するものである。この基本ベクトルは、一つの非方向平均基本ベクトルと複数の方向性基本ベクトルから構成される。ここで、非方向平均基本ベクトルとは、モデリングされた全方向の頭部伝達関数の特徴のうち、音源の方向とは無関係に決定される特徴を代表する基本ベクトルを意味する。一方で、方向性基本ベクトルは、音源の方向により決定される特徴を代表とする基本ベクトルである。特許文献１に記載された技術は、主成分分析とバランスモデル近似技術を用いることで、極零モデルとして少ないパラメータで頭部伝達関数を模擬することを可能としている。

【0019】

また、特許文献２には、頭部伝達関数を表わすパラメータを生成する方法として、頭部インパルス応答信号を周波数領域においてサブバンドに分割し、各サブバンドのパラメータを求める方法が開示されている。この方法では、高速フーリエ変換ビンのグループ化により、周波数領域において少なくとも２つのサブバンドに分割され、サブバンドの信号レベルの二乗平均平方根に基づいて、パラメータが決定される。特許文献２に記載された技術は、サブバンドに分割することにより、頭部伝達関数を用いた演算処理量の低減を可能としている。

【先行技術文献】

【特許文献】

【0020】

【特許文献1】特許第４６８１４６４号公報

【特許文献2】特許第４９２１４７０号公報

【非特許文献】

【0021】

【非特許文献1】F. Asano, Y. Suzuki, and T. Sone, “Role of spectral cues in median plane localization,” J. Acoust. Soc. Am., Vol.88, pp.159-168 (1989)

【非特許文献2】Y. Haneda, S. Makino, and Y. Kaneda, ”Common-Acoustical-Pole and Zero Modeling of Head-Related Transfer Functions,”IEEE Trans. Vol.7, No.2 (1999)

【発明の概要】

【発明が解決しようとする課題】

【0022】

前記した式（２）で表される頭部伝達関数の有限インパルス応答モデルによるモデリングにおいては、５１２次や１０４８次等の高次であるためにパラメータ数が多く、それに伴い、頭部伝達関数を用いた音像定位方式における演算量が多くなる問題がある。加えて、この有限インパルス応答モデルによる頭部伝達関数の模擬では雑音を考慮していないため、雑音の影響を低減するために、同期加算等の前処理が必要となる。しかし、この場合にも雑音の白色性もしくは平均が０となることを仮定しており、低減効果は限定的である。
同様に、非特許文献１，２及び特許文献１，２に記載されたいずれの方法による頭部伝達関数の模擬においても、雑音の影響が考慮されていない。

【0023】

また、従来の技術のうち、十分に高次の有限インパルス応答モデルによる頭部伝達関数のモデリングでは、音像定位知覚の手がかりとなる頭部伝達関数の特徴量が十分に保存されているものの、非特許文献１，２及び特許文献１，２に記載された技術では、少なくともスペクトラルキューが保存されるとは限らない。ここで、スペクトラルキューとは、図１０に示すような頭部伝達関数の周波数特性上の特定のピークＰ１やノッチＮ１，Ｎ２であり、音像を空間内に模擬的に配置させる音像定位技術などにおいては、その正確なモデリングが重要とされている。

【0024】

本発明は、以上のような問題点に鑑みてなされたものであり、音像定位知覚に係る特徴量を保存し且つ雑音の影響を考慮した頭部伝達関数のモデルであって音像定位技術などに利用した際に演算量を低減できるモデルを求める頭部伝達関数のモデリング装置を提供することを課題とする。

【課題を解決するための手段】

【0025】

前記課題を解決するために、本発明に係る頭部伝達関数のモデリング装置は、スピーカに印加された入力信号とスピーカから発せられた音声をマイクロフォンで測定して得られた出力信号とを入出力データとして用いて、漸近推定法により頭部伝達関数をモデリングする頭部伝達関数のモデリング装置であって、高次モデル推定手段と、低次元化手段と、を備え、低次元化手段が、周波数伝達関数算出手段と、低次モデル推定手段と、低次モデル探索手段と、を備えることとした。

【0026】

かかる構成によれば、頭部伝達関数のモデリング装置は、高次モデル推定手段によって、前記入出力データを用いて、予め定められた高次のモデル次数を有した頭部伝達関数及び雑音モデルについての高次モデルのパラメータを予測誤差法により推定する。
この高次モデル推定手段によって、頭部伝達関数の測定における雑音を考慮した、精度の良いモデルを求めることが可能となる。
そして、頭部伝達関数のモデリング装置は、低次元化手段によって、推定された高次モデルと、周波数領域における評価関数である対数尤度関数とを用いて最尤推定値を導出することで前記高次モデルを低次元化する。
ここで、低次元化手段では、周波数伝達関数算出手段によって、前記高次モデルの周波数伝達関数を求める。また、低次元化手段では、低次モデル推定手段によって、前記対数尤度関数を最小化することで前記高次のモデル次数よりも低い次数の低次モデルの推定値を求める。
そして、低次元化手段では、低次モデル探索手段によって、前記低次の次数を更新して前記対数尤度関数を最小化させることを繰り返すことでそれぞれ推定された各頭部伝達関数の低次モデルと、高次参照モデルと、の間の音像定位知覚に係る特徴量の誤差をそれぞれ求め、前記特徴量の誤差が予め定められた許容条件を満たし且つ最低次数となるときの低次モデルを探索する。

【0027】

本発明は、頭部伝達関数のモデリング装置の各手段が処理を実行する頭部伝達関数のモデリング方法で実現することもできる。
本発明は、コンピュータを、頭部伝達関数のモデリング装置の各手段として動作させる頭部伝達関数のモデリングプログラムで実現することもできる。

【発明の効果】

【0028】

本発明によれば、漸近推定法により頭部伝達関数をモデリングすることで、頭部伝達関数として、従来の５１２次や１０２４次といった高次モデルの参照モデルとの誤差が許容される範囲で音像定位知覚に係る特徴量を保存し、従来よりもパラメータ数の少ない、雑音を考慮した低次の頭部伝達関数モデルを求めることができる。そのため、この頭部伝達関数として推定された低次モデルを音像定位技術などの制御対象として用いたときに演算量を低減することができる。

【図面の簡単な説明】

【0029】

【図1】本発明の実施形態に係る頭部伝達関数のモデリング装置の構成を模式的に示すブロック図である。

【図2】図１の入出力データの測定の様子を模式的に示す図面である。

【図3】図１の入出力データの信号波形の一例であって、（ａ）はスピーカへの入力信号の抜粋、（ｂ）はマイクロフォンで測定された出力信号の抜粋を示している。

【図4】頭部伝達関数のモデルの構造を示すブロック線図であり、（ａ）はＦＩＲモデルの構造、（ｂ）はＡＲＸモデルの構造を示している。

【図5】本発明の実施形態に係る頭部伝達関数のモデリング方法を含む音像定位制御の流れを示すフローチャートである。

【図6】図５のＡＲＸモデルの低次元化処理の一例を示すフローチャートである。

【図7】推定された頭部伝達関数の低次モデルと高次参照モデルとの間のスペクトル歪の一例を示すグラフである。

【図8】頭部伝達関数の周波数特性の一例を示すグラフであって、（ｂ）は実施例、（ａ）は実施例と比較例とを重ねて示している。

【図9】本実施形態で推定された頭部伝達関数を適用した音像定位制御の一例を模式的に示すブロック図である。

【図10】スペクトラルキューの一例を示す図である。

【発明を実施するための形態】

【0030】

以下、図面を参照して本発明の頭部伝達関数のモデリング装置を実施するための形態（以下「実施形態」という）について詳細に説明する。

【0031】

図１に示す頭部伝達関数のモデリング装置１は、スピーカに印加された入力信号とスピーカから発せられた音声をマイクロフォンで測定して得られた出力信号とを入出力データとして用いて、漸近推定法により頭部伝達関数をモデリングするものである。

【0032】

＜入出力データ＞
図１に示す入出力データを事前に測定する際には、図２に例示するように音響無響室において、所定方向に設置した例えば１台のスピーカＳＰに対して入力信号ｕ（ｋ）を印加する。ここで、ｋは、音声の連続時間信号のサンプリングを行うときの時間間隔（サンプリング周期）に対応付けられたサンプルを識別する変数である。ｋは離散値であり、その個数がデータ数である。そして、スピーカＳＰから発せられた音声を、ダミーヘッドＤの耳に当たる位置に設置したマイクロフォンで測定する。このときに測定された信号を出力信号ｙ（ｋ）とする。

【0033】

図２では、スピーカＳＰをダミーヘッドＤにとっての正面から右９０°方向（つまり左２７０°方向）に設置してダミーヘッドＤの右耳に向けているが、これは一例である。ダミーヘッドＤを水平面内でスピーカＳＰに対して所定角度だけ相対的に回転させることで、様々な方向から左耳及び右耳に設置したマイクロフォンで音声を測定することができる。なお、スピーカＳＰとダミーヘッドＤとの相対的な距離や高さも可変である。

【0034】

入力信号ｕ（ｋ）としては、擬似白色信号であるＭ系列信号を用いることができる。
一例として、図３（ａ）にシフトレジスタ長が１５（サンプル数＝２¹⁵−１）のＭ系列信号の波形を示す。なお、図３（ａ）には、Ｍ系列信号の途中の一部分（１００サンプル分）の波形を示した。また、図３（ｂ）には、ダミーヘッドＤの正面から左３０°方向に設置したスピーカＳＰに、そのＭ系列を印加したときに測定された出力信号ｙ（ｋ）の一部波形を示す。ただし、図３（ａ）及び図３（ｂ）に示す波形はサンプリング周波数を４８ｋＨｚとした場合に得られたものである。

【0035】

＜漸近推定法＞
漸近推定法（asymptotic method）は、例えばプラント制御のためのモデリング法として公知の手法である。漸近推定法では、まず、システム同定実験によって得られた入出力データに対して雑音モデルを考慮した高次（例えばｎ次とする）モデルのパラメータを推定し、その後に、漸近理論に基づき周波数特性を考慮して低次元化を行なう。

【0036】

［頭部伝達関数のモデリング装置の構成］
頭部伝達関数のモデリング装置１（以下、単にモデリング装置１という）は、一般的なコンピュータと同様に、例えば、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ＨＤＤ（Hard Disk Drive）、入出力インタフェース等を備えている。ここでは、モデリング装置１は、図１に示すように、高次モデル推定手段１０と、高次数設定手段１１と、低次元化手段２０と、を備えている。

【0037】

高次モデル推定手段１０は、入出力データを用いて、予め定められた高次（ｎ次）のモデル次数を有した頭部伝達関数及び雑音モデルについての高次モデルを予測誤差法により推定するものである。

【0038】

ここで、雑音モデルをもつモデルとしては、図４（ｂ）に示すＡＲＸ（AutoRegressive eXogenous）モデル等の公知のモデルを用いることができる。ＡＲＸモデルでは雑音モデル（図中の１／Ａ（ｑ））を考慮しており、本実施形態では、マイクロフォンで測定された出力信号に雑音の影響が含まれていることとして扱う。

【0039】

なお、従来の一般的な頭部伝達関数のモデリングは、図４（ａ）に示すＦＩＲモデルを用いており、音像定位にはＦＩＲフィルタが実装されてきた。ＦＩＲモデルでは雑音モデル（１／Ａ（ｑ））を考慮しておらず、同期加算等の前処理を行うことで、雑音の影響が含まれていないものとした信号を扱っている。

【0040】

モデリング装置１は、後の処理で低次元化を行なうことを前提にして、高次モデル推定手段１０によって、高次モデルのパラメータを推定する。ここで、高次モデルのパラメータを推定するのは、一般的に、モデルが高次であるほど精度が良くなるためである。なお、雑音モデルを考慮しない従来技術では、頭部伝達関数を例えば５１２次や１０４８次等の高次有限インパルス応答モデルとしていた。本実施形態では、高次の次数としては、経験的に三桁の数値を想定している。

【0041】

また、本実施形態において前提とする予測誤差法とは、予測値に基づく予測誤差から構成される評価規範を最小化する推定値を計算する、パラメータ推定法の総称である。この予測誤差法としてはシステム同定に用いられる一般的な手法を用いることができる。予測誤差の大きさの測度として２次関数を用いると、予測誤差法は最小二乗法となる。
ここでは、高次モデル推定手段１０は、高次ＡＲＸモデルを用いてシステム同定理論に基づき最小二乗法でパラメータを推定することとする。

【0042】

高次モデル推定手段１０は、予め測定された入出力データ｛ｕ（ｋ），ｙ（ｋ）；ｋ＝１，２，…，Ｎ｝によって、次の式（６）で定義されるＡＲＸモデルのパラメータを推定する。ここで、ｕ（ｋ）は入力信号、ｙ（ｋ）は出力信号、Ｎはデータ数である。

【0043】

【数6】

【0044】

ただし、式（６）においてＡ^h（ｑ）、Ｂ^h（ｑ）は、システムを構成する要素（以下、単にモデルと呼ぶ）であり、パラメータ｛ａ_i｝，｛ｂ_i｝（ｉ＝１〜ｎ）を用いて、それぞれ以下の式（７）、式（８）で表される。

【0045】

【数7】

【0046】

【数8】

【0047】

また、前記した式（６）において、外乱ｗ（ｋ）としては、平均値が０で分散がσ²_wである正規性白色雑音を仮定している。
式（８）において、Ｌはむだ時間（dead time）を表す。
上記各式において、上添え字ｈはｈｉｇｈの省略形であって、モデルＡ^h（ｑ）、Ｂ^h（ｑ）の次数ｎが十分に高いこと、すなわち高次モデルであることを意味しており、ｈは変数ではない。

【0048】

本実施形態では、高次数設定手段１１が、モデリング装置１の外部から入力されるモデル次数を高次モデル推定手段１０に設定することとした。このモデル次数（高次の次数ｎ）をモデリング装置１の内部に記憶させておき、高次モデル推定手段１０が処理の際に読み出す。ここでは、高次の次数ｎは例えば１００であるものとする。なお、高次の次数ｎは、処理の際にその都度、モデリング装置１の外部から入力するようにしてもよい。

【0049】

高次モデル推定手段１０は、ＡＲＸモデルのパラメータとして、前記式（７）に示すパラメータ｛ａ_i｝（ｉ＝１〜１００）と、前記式（８）に示すパラメータ｛ｂ_i｝（ｉ＝１〜１００）とを最小二乗法により推定する。

【0050】

高次モデル推定手段１０によって前記式（７）から推定された多項式をＡ^^h（ｑ）と表記し、前記式（８）から推定された多項式をＢ^^h（ｑ）と表記する。なお、本明細書において、ある文字の右に配置された記号「^」は、その直前の文字の上に配置されたハット記号を意味することとする。このハット記号は推定値であることを表すものである。

【0051】

このように表記する場合、制御対象とする伝達関数Ｇ（ｑ）についての高次モデルは、Ｇ^^h（ｑ）と表記され、次の式（９ａ）で表される。また、雑音モデルＨ（ｑ）についての高次雑音モデルは、Ｈ^^h（ｑ）と表記され、次の式（９ｂ）で表される。

【0052】

【数9】

【0053】

なお、ＦＩＲモデルの場合、正規性白色雑音もしくは平均０の雑音を仮定しているので、式（９ａ）に対応した関係式がＧ（ｑ）＝Ｂ（ｑ）のように記述され、式（９ｂ）に対応した関係式がＨ（ｑ）＝１のように記述される。比較的高次となるＦＩＲモデルは、後記するように、低次ＡＲＸモデルのリファレンス（以下、高次参照モデルという）として用いることにする。

【0054】

低次元化手段２０は、高次モデル推定手段１０で推定された高次モデル｛Ｇ^^h（ｑ），Ｈ^^h（ｑ）｝と、周波数領域における評価関数である対数尤度関数とを用いて最尤推定値を導出することで高次モデルを低次元化するものである。
この低次元化手段２０は、図１に示すように、周波数伝達関数算出手段２１と、低次モデル推定手段２２と、低次モデル探索手段２３と、を備えている。

【0055】

周波数伝達関数算出手段２１は、高次モデル推定手段１０で推定された既知の高次モデル｛Ｇ^^h（ｑ），Ｈ^^h（ｑ）｝についての周波数伝達関数Ｇ^^h(ｅ^jω)，Ｈ^^h(ｅ^jω)をそれぞれ求めるものである。この周波数伝達関数算出手段２１は、例えば、高次モデルの伝達関数Ｇ^^h（ｑ）から、その周波数応答である高次の周波数伝達関数Ｇ^^h(ｅ^jω)を算出する。

【0056】

低次モデル推定手段２２は、周波数領域における評価関数である対数尤度関数を最小化することで、高次モデル推定手段１０で推定された高次モデルのモデル次数（高次の次数）よりも低い次数の低次モデルの推定値を求めるものである。

【0057】

ここでは、高次モデル推定手段１０に設定されたモデル次数をｎ（＝１００）としている。また、推定しようとする未知の低次モデルのモデル次数をｍ（０＜ｍ＜ｎ）とおく。また、このｍ次の低次モデルのことを、｛Ｇ^^l（ｑ），Ｈ^^l（ｑ）｝と表記する。ここで、上添え字ｌはｌｏｗの省略形であって、Ｇ^^l（ｑ），Ｈ^^l（ｑ）の次数ｍがｎよりも低いこと、すなわち低次モデルであることを意味しており、ｌは変数ではない。
なお、Ｇ^^l（ｑ）は、制御対象とする伝達関数Ｇ（ｑ）についての低次モデルの推定値であり、Ｈ^^l（ｑ）は、雑音モデルＨ（ｑ）についての低次雑音モデルの推定値である。

【0058】

この場合、さらに、推定しようとする未知の低次モデルの伝達関数Ｇ^^l（ｑ）についての周波数伝達関数をＧ^^l(ｅ^jω)とおき、同様に、推定しようとする未知の低次雑音モデルＨ^^l（ｑ）についての周波数伝達関数をＨ^^l(ｅ^jω)とおく。

【0059】

前記した対数尤度関数は、次の式（１０）のＶで表される。なお、漸近理論では、高次モデル（Ｇ^^h(ｅ^jω)）は周波数領域において近似的に正規分布に従うので低次モデル（Ｇ^^l(ｅ^jω)）は漸近的に最尤推定値になる。

【0060】

【数10】

【0061】

式（１０）において、Φ_u（ω）は、入力信号ｕ（ｋ）のパワースペクトル密度関数である。また、Φ_v（ω）は、雑音ｖ（ｋ）のパワースペクトル密度関数である。この雑音ｖ（ｋ）は次の式（１１）で表わされる。

【0062】

【数11】

【0063】

なお、ここで、雑音ｖ（ｋ）は、外乱ｗ（ｋ）として仮定した正規性白色雑音が、雑音成形フィルタ（高次モデル推定手段１０で推定された既知の高次雑音モデルＨ^^h（ｑ））を通過した後の雑音のことである。

【0064】

低次モデル推定手段２２は、推定しようとする未知の低次モデルのモデル次数ｍ（０＜ｍ＜ｎ）が低次モデル探索手段２３によって所定値に設定されたときに、前記式（１０）のＶを最小化させる演算処理を行い、モデル次数ｍが当該設定値のときに、式（１０）のＶを最小化する低次モデルの推定値Ｇ^^l（ｑ）を得る。この低次モデルの推定値Ｇ^^l（ｑ）はモデリング装置１の出力候補である。なお、対数尤度関数Ｖの最小化には、非線形最適化問題を解く必要がある。その求解法は限定しないが、一例について後で説明を行う。

【0065】

＜低次モデル探索手段２３＞
低次モデル探索手段２３は、推定しようとする未知の低次モデルのモデル次数ｍを更新して対数尤度関数Ｖを最小化させることを繰り返すことでそれぞれ推定された各頭部伝達関数の低次モデルＧ^^l（ｑ）と、高次参照モデルと、の間の音像定位知覚に係る特徴量の誤差をそれぞれ求め、この特徴量の誤差が予め定められた許容条件を満たし且つ最低次数となるときの低次モデルを探索するものである。ここで探索された最低次数の低次モデルＧ^^l（ｑ）がモデリング装置１の出力である。

【0066】

本実施形態では、高次参照モデルの一例として、音響分野で一般的に用いられているＦＩＲモデルを用いた。高次参照モデルとしたＦＩＲモデルの次数は、十分な収束を考慮して５１２次とした。

【0067】

以下では、音像定位知覚に係る特徴量としては、一例として、頭部伝達関数の周波数特性上のピーク及びノッチであるスペクトラルキューの位置であることとする。それは、音像を空間内に模擬的に配置させる音像定位技術などにおいては、スペクトラルキューの正確なモデリングが重要とされているからである。

【0068】

本実施形態では、高次ＡＲＸモデルを低次元化して低次ＡＲＸモデルを推定する際に、スペクトラルキューを保存するように構成した。ここで、スペクトラルキューの保存とは、推定しようとする低次ＡＲＸモデルのスペクトラルキューが、高次参照モデル（例えば高次のＦＩＲモデル）のスペクトラルキューを所望の正確さで再現できるように予め定められた許容条件を満たすことをいう。再現性の精度は、音像定位などに応用するときに期待する演算量の所望の低減効果に応じて適宜設定される。両者のスペクトラルキューの中心周波数のずれを所望の許容できる範囲に抑えることにより、頭部伝達関数として推定された低次モデルが所望の確度を有しつつ、その低次モデルを音像定位技術などの制御対象として用いたときの演算量を低減することができる。

【0069】

本実施形態では、低次モデル探索手段２３が、低次モデル推定手段２２でそれぞれ推定された各頭部伝達関数の低次モデルＧ^^l（ｑ）と高次参照モデルとの間のスペクトル歪（ＳＤ:spectral distortion）をそれぞれ求め、スペクトル歪が予め定められた第１閾値以下、且つスペクトラルキューの位置についての高次参照モデルとの間の誤差が予め定められた第２閾値以下の条件を満たし且つ最低次数となるときの低次モデルを探索することとした。ここで、スペクトル歪ＳＤとは、２つの伝達関数の一致度を判定するために、その振幅特性の差をすべての周波数成分で評価した物理指標のことである。このように構成することで、演算量を低減し且つ、推定しようとする低次モデルの精度を保証することができる。

【0070】

そこで、本実施形態では、低次モデル探索手段２３が、図１に示すように、スペクトラル歪算出手段２４と、スペクトラル歪判定手段２５と、音響特徴量算出手段２６と、音響特徴量判定手段２７と、を備えることとした。

【0071】

スペクトラル歪算出手段２４は、２つの頭部伝達関数Ｈ_X，Ｈ_Yが与えられたときに、そのスペクトル歪ＳＤを次の式（１２）によって算出するものである。

【0072】

【数12】

【0073】

式（１２）において、Ｎは頭部伝達関数におけるデータ数を表わす。また、頭部伝達関数Ｈ_Xは例えば推定された頭部伝達関数の低次モデルＧ^^l（ｑ）の周波数伝達関数を表わし、頭部伝達関数Ｈ_Yは例えば高次参照モデルの周波数伝達関数を表わす。

【0074】

スペクトラル歪判定手段２５は、スペクトラル歪算出手段２４によって算出されたスペクトル歪ＳＤが第１閾値以下であるか否かを判定するものである。ここで、第１閾値は特に限定されるものではない。ただし、スペクトル歪ＳＤは、およそ３ｄＢを下回っていれば、２つの頭部伝達関数がよく一致しているとみなすことができるとされている。そこで、本実施形態では、典型的な一例として第１閾値を３ｄＢとした。これにより、スペクトル歪ＳＤが３ｄＢより大きい場合、そのときに設定されているモデル次数ｍにおいて推定された低次モデルＧ^^l（ｑ）は、ＮＧであってモデリング装置１の最終出力とはならない。よって、前記式（１０）や前記式（１２）の演算において、ｍとして取り得るすべての値を選択する前に、モデル次数ｍの更新につれてスペクトル歪ＳＤが３ｄＢよりも大きくなった時点で演算を終了させることができる。

【0075】

音響特徴量算出手段２６は、推定しようとする低次モデルＧ^^l（ｑ）についての周波数伝達関数Ｇ^^l(ｅ^jω)におけるスペクトラルキューの中心周波数と、高次参照モデルについての周波数伝達関数におけるスペクトラルキューの中心周波数との差分を音響特徴量の誤差として算出するものである。

【0076】

音響特徴量判定手段２７は、音響特徴量算出手段２６によって算出された音響特徴量の誤差（スペクトラルキューの中心周波数のずれ）が第２閾値以下であるか否かを判定するものである。ここで、第２閾値は特に限定されるものではない。ただし、本実施形態では、典型的な一例として第２閾値を離散周波数ビン１サンプルとした。つまり、スペクトラルキューの中心周波数において１サンプルの誤差範囲を許容した。具体的には、例えばサンプリング周波数を４８ｋＨｚ、スペクトルをデータ数Ｎ（５１２）の離散フーリエ変換により求めた場合、１サンプルの区間は９３．７５Ｈｚとすることができる。これにより、スペクトラルキューの中心周波数のずれが１サンプルより大きい場合、そのときに設定されているモデル次数ｍにおいて推定された低次モデルＧ^^l（ｑ）は、ＮＧであってモデリング装置１の最終出力とはならない。よって、スペクトラルキューの中心周波数のずれを算出する処理において、ｍとして取り得るすべての値を選択する前に、モデル次数ｍの更新につれてスペクトラルキューの中心周波数のずれが１サンプルよりも大きくなった時点で演算を終了させることができる。

【0077】

図１の低次モデル探索手段２３は、前記したように、推定しようとする未知の低次モデルのモデル次数ｍを更新して前記した式（１０）の対数尤度関数Ｖを最小化させることを繰り返しつつ、スペクトル歪ＳＤとスペクトラルキューの中心周波数のずれについての閾値判定を行う。この際に、モデル次数ｍの更新の仕方は特に限定されるものではない。

【0078】

例えば、スペクトル歪ＳＤの閾値判定と、スペクトラルキューの中心周波数のずれについての閾値判定との一方についてモデル次数ｍを更新しながら処理を行った後に、他方についてモデル次数ｍを更新しながらの閾値判定処理を行ってもよい。または、モデル次数ｍをある値に設定したときに、両方の閾値判定処理を行ってから、モデル次数ｍを更新するようにしてもよい。

【0079】

また、次数ｍを単調に減少させたり、単調に増加させたり、増減させたりしてもよい。
また、次数ｍを１ずつシフトしてもよいし、２以上の所定値ずつシフトしてもよい。必ずしも毎回同じシフト数にする必要もなく、例えば１０ずつシフトした後で、その間を１ずつシフトしてもよい。
また、前記したようにモデリング装置１の最終出力とはならないことが明らかならば、モデル次数ｍの更新の際に、値を取り得る全てのｍ（０＜ｍ＜ｎ）を必ずしも選択しなくてもよい。

【0080】

本実施形態では、一例として、低次モデル探索手段２３が、低次のモデル次数ｍを高次のモデル次数ｎの側から降順に更新してスペクトル歪ＳＤが第１閾値より大きくなった場合（これを以下、反復停止条件という）、対数尤度関数Ｖの最小化処理を停止し、スペクトル歪ＳＤが第１閾値（３ｄＢ）以下であったときの次数を起点として低次のモデル次数ｍを昇順に更新してスペクトラルキューの位置についての高次参照モデルとの間の誤差が第２閾値以下になったときの次数を最低次数として決定することとした。このようにすることで、モデル次数ｍの最低次数を効率よく求めることが可能である。

【0081】

（式（１０）の対数尤度関数Ｖの最小化の求解法）
ここで説明する求解法には、下記のＡ−１，Ａ−２，Ａ−３の３つの手続きがある。

【0082】

Ａ−１．入力のフィルタリング
高次モデル推定手段１０で推定された既知の高次雑音モデルＨ^^h（ｑ）を用いて、次の式（１３）のように入力信号ｕ（ｋ）をフィルタリングする。この高次雑音モデルＨ^^h（ｑ）のフィルタを通過した信号ｕ_f（ｋ）は、頭部伝達関数の高次モデルＧ^^h（ｑ）の入力信号として用いられる。

【0083】

【数13】

【0084】

Ａ−２．出力の計算
前記手続きＡ−１で得られた式（１３）に示すフィルタ通過信号ｕ_f（ｋ）を、頭部伝達関数の高次モデルＧ^^h（ｑ）への新たな入力信号として、高次モデルＧ^^h（ｑ）の出力信号ｙ^_f（ｋ）を次の式（１４）により計算する。

【0085】

【数14】

【0086】

Ａ−３．低次モデルのパラメータ推定
手続きＡ−１，Ａ−２によって、頭部伝達関数の高次モデルＧ^^h（ｑ）の新しい入出力信号｛ｕ_f（ｋ），ｙ^_f（ｋ）；ｋ＝１，２，…，Ｎ｝が得られたならば、この新しい入出力信号を用いて、高次モデルＧ^^h（ｑ）を低次元化した低次モデルのパラメータを出力誤差法により推定する。このとき、出力誤差法の損失関数Ｖ^OEは、次の式（１５）で表される。

【0087】

【数15】

【0088】

ここで、式（１５）におけるＧ^^l（ｑ）が、求めるべき頭部伝達関数の低次モデルであり、次の式（１６）で表わされる。また、低次雑音モデルＨ^^l（ｑ）は次の式（１７）で表される。

【0089】

【数16】

【0090】

ただし、式（１６）において、Ａ^l（ｑ），Ｂ^l（ｑ）は式（１８），式（１９）でそれぞれ表され、式（１７）において、Ｃ^l（ｑ），Ｄ^l（ｑ）は式（２０），式（２１）でそれぞれ表される。

【0091】

【数17】

【0092】

［音像定位制御の流れ］
ここでは、頭部伝達関数およびそのモデルの応用の一例として、音響信号の音像を空間内に擬似的に配置させるような音像定位技術を挙げて説明する。具体的には、モデリング装置１による頭部伝達関数のモデリング方法を含む音像定位制御の全体の流れについて図５を参照（適宜図１〜図３参照）して説明する。
まず、モデリング装置１の処理をする前に、システム同定を行うために必要な入出力データを測定する（ステップＳ１００）。測定方法は、図２を参照して説明した方法を用いることができる。入出力データの具体例は図３（ａ）及び図３（ｂ）に示されている。

【0093】

続いて、モデリング装置１の処理として、入出力データを用いて、漸近推定法により頭部伝達関数をモデリングする（ステップＳ２００）。
漸近推定法による処理（ステップＳ２００）を概説すると、まず、高次ＡＲＸモデルのパラメータを推定する処理を行い（ステップＳ２１０：高次モデル推定ステップ）、その後に、漸近理論に基づきＡＲＸモデルの低次元化処理を行う（ステップＳ２２０：低次元化ステップ）。

【0094】

より詳細には、ステップＳ２１０では、モデリング装置１において、高次モデル推定手段１０が、入出力データを用いて、予め定められた高次（ｎ次）のモデル次数を有した頭部伝達関数及び雑音モデルについての高次モデル（Ｇ^^h（ｑ），Ｈ^^h（ｑ））のパラメータ｛ａ_i｝，｛ｂ_i｝を予測誤差法により推定する。
また、ステップＳ２２０では、モデリング装置１において、低次元化手段２０が、推定された高次モデル（Ｇ^^h（ｑ），Ｈ^^h（ｑ））と、前記式（１０）に示す対数尤度関数Ｖとを用いて最尤推定値を導出することで高次モデルを低次元化する。そして、モデリング装置１は頭部伝達関数の低次モデルを出力する。なお、ＡＲＸモデルの低次元化処理の詳細な流れについては後記する。

【0095】

続いて、音像定位の制御対象に対して、モデリング装置１によって推定された頭部伝達関数の低次ＡＲＸモデルを適用する（ステップＳ３００）。この低次ＡＲＸモデルは、従来よりもパラメータ数の少ないモデルとして求められているので、頭部伝達関数として推定された低次モデルを音像定位技術などの制御対象として用いたときの演算量を従来よりも低減することができる。

【0096】

［ＡＲＸモデルの低次元化処理の詳細な流れ］
次に、モデリング装置１の低次元化手段２０によるＡＲＸモデルの低次元化処理の詳細な流れについて図６を参照（適宜図１〜図３及び図５参照）して説明する。
ＡＲＸモデルの低次元化処理（ステップＳ２２０）では、まず、モデリング装置１の周波数伝達関数算出手段２１が、図５のステップＳ２１０で得られた高次（ｎ次）モデル｛Ｇ^^h（ｑ），Ｈ^^h（ｑ）｝の周波数伝達関数Ｇ^^h(ｅ^jω)，Ｈ^^h(ｅ^jω)を求める（ステップＳ２２１）。

【0097】

そして、低次モデル探索手段２３が、モデル次数ｍ（０＜ｍ＜ｎ）を設定する（ステップＳ２２２）。ここでは、一例としてｎ＝１００としているので、ｍ＝９９を設定することとする。そして、低次モデル推定手段２２は、例えばｎ＝１００、ｍ＝９９の場合において、前記した式（１０）に示す対数尤度関数Ｖを最小化し、例えばｍ＝９９の設定値の場合の低次モデルを推定する（ステップＳ２２３）。これにより、低次モデル探索手段２３が、例えばｍ＝９９の設定値の場合の低次モデルの推定値Ｇ^^l（ｑ）を得る。

【0098】

そして、低次モデル探索手段２３は、モデル次数ｍの更新についての反復停止条件が成立したか否かを判定する（ステップＳ２２４）。具体的には、スペクトラル歪算出手段２４が、９９次のモデル（Ｈ_X）と１００次のモデル（Ｈ_Y）のスペクトル歪ＳＤを前記式（１２）によって算出し、スペクトラル歪判定手段２５が３ｄＢ（第１閾値）より大きいと判定した場合、反復停止条件が成立する。

【0099】

一方、反復停止条件が成立していない場合（ステップＳ２２４：Ｎｏ）、低次のモデル次数ｍを更新し（ステップＳ２２５）、ステップＳ２２３に戻る。具体的には、低次モデル探索手段２３が、モデル次数ｍの値を１だけ減算してｍ＝９８を設定した場合、低次モデル推定手段２２は、例えばｎ＝１００、ｍ＝９８の場合において、前記した式（１０）に示す対数尤度関数Ｖを最小化し、例えばｍ＝９８の設定値の場合の低次モデルを推定する（ステップＳ２２２）。このステップＳ２２４でＮｏの場合の処理は以下同様に減算を行う。

【0100】

ステップＳ２２５によって、モデル次数ｍの値をより低くしてモデルを低次元化し続けると、やがて、低次モデル探索手段２３は、反復停止条件が成立したと判定する（ステップＳ２２４：Ｙｅｓ）。具体的には、スペクトラル歪判定手段２５が３ｄＢ（第１閾値）より大きいと判定する。

【0101】

ステップＳ２２４でＹｅｓの場合、低次モデル探索手段２３は、ステップＳ２２６において、その時点のモデル次数ｍの設定値に「１」を加算する。この加算で得られたモデル次数ｍの値は、低次ＡＲＸモデルにおいてスペクトラル歪ＳＤが第１閾値（３ｄＢ）以下であったときの最低次数である。また、このとき、低次モデル探索手段２３は、音像定位知覚に係る特徴量についての参照（高次参照モデル）との誤差が予め定められた許容条件を満たす最低次数から、低次モデルの次数ｍを決定する（ステップＳ２２６）。
具体的には、音響特徴量算出手段２６が、その時点で設定されているモデル次数ｍの低次モデルＧ^^l（ｑ）についての周波数伝達関数Ｇ^^l(ｅ^jω)におけるスペクトラルキューの中心周波数と、高次参照モデル（５１２次ＦＩＲモデル）についての周波数伝達関数におけるスペクトラルキューの中心周波数との差分を音響特徴量の誤差として算出する。そして、音響特徴量判定手段２７は、スペクトラルキューの中心周波数のずれが１サンプル（第２閾値）以下であるか否かを判定する。

【0102】

スペクトラルキューの中心周波数のずれが１サンプル（第２閾値）より大きい場合、その時点のモデル次数ｍの設定値が低過ぎるので、低次モデル探索手段２３は、その時点のモデル次数ｍの設定値に「１」を加算する。なお、この加算で得られたモデル次数ｍの値の場合、スペクトラル歪ＳＤは当然ながら第１閾値（３ｄＢ）以下である。

【0103】

そして、同様にして、音響特徴量算出手段２６が、スペクトラルキューの中心周波数のずれを算出し、音響特徴量判定手段２７は、スペクトラルキューの中心周波数のずれが１サンプル（第２閾値）以下であるか否かを判定する。スペクトラルキューのずれが大きい場合の処理は以下同様にモデル次数ｍの値の加算を行う。やがて、音響特徴量判定手段２７は、ずれが１サンプル以下になったと判定する。この時点のモデル次数ｍの設定値が、低次モデル探索手段２３で本来探索していた最低次数である。そして、モデリング装置１は、その最低次数のモデル次数を有した頭部伝達関数の低次元ＡＲＸモデルを出力する。

【0104】

［ＡＲＸモデルの低次元化処理の具体例］
本発明の効果を確かめるために、ダミーヘッドＤの頭部中心の位置より１．３ｍ離れた位置で、ダミーヘッドＤの正面から左３０°方向に設置した１つのスピーカＳＰから測定信号を印加する実験を行った。そして、ダミーヘッドＤの左耳に内蔵されたマイクロフォンにより収音を行なって図３（ｂ）に例示した出力信号を得た。このときの入出力データを用いて本実施形態に係るモデリング方法で頭部伝達関数をモデリングした。また、高次参照モデルとして５１２次ＦＩＲモデルを求めた。

【0105】

図６のステップＳ２２３〜Ｓ２２５を繰り返しつつモデル次数ｍの設定値を単調減少させながら、ｎ＝１００の高次モデルとの間で各次数におけるスペクトラル歪ＳＤを求めた。このときのモデル次数ｍに対するＳＤの変化のグラフを図７に示す。モデル次数ｍの設定値を１７次まで下げたときにＳＤが３ｄＢよりも大きくなった。つまり、ＳＤが３ｄＢ以下となる次数は１８次以上であった。なお、モデル次数ｍの更新にあわせて前記式（１０）に示す対数尤度関数Ｖを最小化する処理を８３回行った。

【0106】

続いて、スペクトラルキューの中心周波数を、高次参照モデル（５１２次ＦＩＲモデル）との間で１サンプルの範囲内で捉えることのできるモデルを求めた。高次参照モデル（５１２次ＦＩＲモデル）の周波数特性を図１０に示す。

【0107】

図６のステップＳ２２６に対応させて、既に得られている１８次以上のモデルを探索対象としてスペクトラルキューの中心周波数のずれを求めた。１８次、１９次、２０次の場合、スペクトラルキューのずれが許容範囲を超えていたが、２１次の場合、スペクトラルキューのずれを１サンプルの範囲内で捉えることができた。つまり、モデル次数ｍの更新に係るスペクトラルキューのずれを求める処理は４回だけ行った。

【0108】

上記実験から得られた２１次のＡＲＸモデルの周波数特性を図８（ｂ）に破線で示す。また、図８（ａ）には、図８（ｂ）に示した２１次のＡＲＸモデルの周波数特性（破線）と、図１０に示した高次参照モデル（５１２次ＦＩＲモデル）の周波数特性（実線）とを重ねて表示した。低周波側においてひずみが顕著となるが、これは、モデルの違いと、ゲインの小さい帯域では相対的にＳＤが大きくなることが原因として考えられる。

【0109】

図８（ａ）に示すように、２１次のＡＲＸモデルでは、スペクトラルキューが保存されていることが分かる。加えて、従来用いられてきたＦＩＲモデルが５１２次であることを考慮する高次参照モデルに比べてパラメータ数を大幅に減少させることができることを確かめた。

【0110】

なお、上記実験について、１方向の測定例を挙げたが、実際には、収音し終えたらダミーヘッドを５°時計周りに回転させ、同様にして収音する、という手順を繰り返すことにより、水平面５°間隔７２方向の頭部伝達関数を測定した。そして、水平面７２方向から左耳までの頭部伝達関数をモデリングしたとき、各方位において妥当だと考えられる次数の平均次数を求めると２０次となった。

【0111】

［頭部伝達関数の低次モデルを適用する音像定位制御の具体例］
ここでは、本実施形態で推定された頭部伝達関数の低次モデルを適用する音像定位制御の具体例について説明する。
頭部伝達関数を利用したシステムとして、トランスオーラルシステムと呼ばれる、三次元音響を実現するためのシステムが知られている。図９に、制御点（以下、その識別子をｉ、ただしｉ＝１，…，ｍとする）及び２次音源（以下、その識別子をｊ、ただしｊ＝１，…，ｎとする）を有するトランスオーラル再生システムのブロック図を示す。なお、ｍはｎと等しくてもよい。

【0112】

ここで、制御点は、例えば、図９に示す聴取者（リスナー）９０の右耳位置Ｒや左耳位置Ｌである。一例として、識別子ｉによって、例えばリスナー９０の右耳位置、左耳位置の順に制御点を識別するものとした。シンプルな例では１人のリスナー９０を想定してｍ＝２とすればよいが一般化して説明する。

【0113】

２次音源は、図９に示すスピーカである。ここでは、一例として、リスナー９０の右耳位置Ｌ側から順にスピーカＳＰ_jを識別子ｊによって識別するものとした。音源ＳＳはスピーカＳＰ_j（ｊ＝１，…，ｎ）に信号を出力するものである。なお、音源ＳＳの個数は特に限定されない。スピーカＳＰは、例えばラウドスピーカである。

【0114】

図中、システムの要素Ｇ_ij（ｑ）は、シフトオペレータｑを用いた制御対象の頭部伝達関数を表し、ｊ番目のスピーカ（２次音源）からｉ番目の制御点（耳位置）への音響伝達関数を表す。なお、モデリング装置１で推定する１つの低次元ＡＲＸモデルが１つのＧ_ij（ｑ）に相当する。また、要素Ｘ_i（ｑ）は、各制御点での所望伝達関数を表す。

【0115】

また、要素Ｈ_ji（ｑ）は、クロストーク・キャンセレーションのための制御器として働く。一般に、スピーカを用いてバイノーラル信号を提示する場合、スピーカから同側耳までの信号の伝搬に加え、対側耳への漏洩(クロストーク)も発生する。従って、このクロストークを抑圧し、所望信号のみをそれぞれの耳に伝送する補償処理が必要となってくる。この処理のことをクロストーク・キャンセレーションという。
なお、前記した要素Ｘ_i（ｑ）は、制御器Ｈ_ji（ｑ）によってクロストークを抑圧後に、ｉ番目の各制御点においてリスナー９０の右耳又は左耳にだけ聴かせたい音声信号の伝達関数を表す。

【0116】

従って、システムの入出力信号は、次の式（２２）〜式（２６）のような関係で表される。

【0117】

【数18】

【0118】

ここで、システムの入力信号をスカラーのｕ（ｋ）として、システムの出力信号ｙ（ｋ）を式（２３）に示すように、制御点数に対応してｍ個の要素を有した一次元行列（列ベクトル）で表すこととしている。なお、ｙ（ｋ）の列ベクトルは、式（２３）において行ベクトルの転置Ｔにより表されている。また、ｋ＝１，２，…，Ｎとすると共に、Ｎはデータ数であるものとする。
式（２４）は、制御点数に対応してｍ個の要素（Ｘ_i（ｑ））を有した列ベクトルで表される。なお、この列ベクトルは行ベクトルの転置Ｔにより表されている。
式（２５）は、制御点数と２次音源数とに対応してｍ×ｎ個の要素（Ｇ_ij（ｑ））を有したベクトル（行列）で表される。
式（２６）は、２次音源数と制御点数とに対応してｎ×ｍ個の要素（Ｈ_ji（ｑ））を有したベクトル（行列）で表される。

【0119】

このシステムにおける式（２２）の左辺に示す所望出力信号は、クロストーク・キャンセレーション後、入力信号ｕ（ｋ）に対して前記式（２４）に示す所望伝達関数が作用された信号となるため、次の式（２７）のように記述される。

【0120】

【数19】

【0121】

このようにシフトオペレータｑを用いると、時間領域での畳み込み演算が行列積の形で記述可能となる。そのため、前記式（２６）で定義されたシステムの制御器を求めるには、式（２７）と前記式（２２）とから代数学的な逆行列演算を行えばよい。その結果、次の式（２８）で記述されるようなシステムの制御器を設計することができる。

【0122】

【数20】

【0123】

ただし、式（２８）に示す制御器では不安定となる。これを解決するため、この不安定な制御器を一旦設計した後、その制御器を構成する各伝達関数から不安定極を持つ伝達関数を括り出すことが考えられる。そして、これを遅れ逆システムとして近似することにより、安定な制御器を実現することが可能である。本実施形態では、頭部伝達関数のモデリングは、図４（ｂ）に示すＡＲＸモデルを用いており、音像定位にはＩＩＲフィルタが実装されることになる。

【0124】

以上説明したように、本実施形態に係る頭部伝達関数のモデリング装置によれば、音像定位知覚にかかる特徴量であるスペクトラルキューを保存しつつ、頭部伝達関数を少ないパラメータ数でモデリングすることが可能となる。また、頭部伝達関数の測定における雑音を考慮しているため、より精緻なモデルとなることが期待される。さらに、パラメータ数が少ないため、頭部伝達関数を用いた音像定位方式などにおける演算量を低減することが可能となる。

【0125】

以上、実施形態に基づいて本発明を説明したが、本発明はこれに限定されるものではない。例えば、低次モデル探索手段２３が低次モデルを探索する際の指標に用いる音像定位知覚に係る特徴量としてスペクトラルキューを例示したが、頭部伝達関数に含まれる他の特徴量、例えば両耳間時間差やレベル差を用いてもよい。

【0126】

また、頭部伝達関数のモデリング装置１は、電子回路が各種電子部品や半導体デバイス等によってハードウェア的に構築された回路であってもよいし、当該装置１の各構成の処理を汎用的または特殊なコンピュータ言語によって記述した頭部伝達関数のモデリングプログラムとこれを処理するＣＰＵの協働によって実現するものであってもよい。

【0127】

また、高次参照モデルの一例として５１２次のＦＩＲモデルを挙げて具体的に説明したが、高次参照モデルは、これに限らず、ＦＩＲ以外の例えばＡＲＸモデルであっても、 雑音モデルを考慮した高次モデルであっても構わない。
また、頭部伝達関数およびそのモデルの応用の一例として、音像定位技術を挙げて具体的に説明したが、本発明は、音像定位に限らず、例えば、ラウドネスメーターなど頭部伝達関数を利用した技術全般に適用することができる。

【産業上の利用可能性】

【0128】

本発明に係る頭部伝達関数のモデリング装置は、ヘッドホンを用いた音響再生技術、スピーカによる音響再生技術全般に利用することができる。

【符号の説明】

【0129】

１ 頭部伝達関数のモデリング装置
１０ 高次モデル推定手段
１１ 高次数設定手段
２０ 低次元化手段
２１ 周波数伝達関数算出手段
２２ 低次モデル推定手段
２３ 低次モデル探索手段
２４ スペクトラル歪算出手段
２５ スペクトラル歪判定手段
２６ 音響特徴量算出手段
２７ 音響特徴量判定手段

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】