(57)【要約】【課題】左右に２チャンネルからなるデジタル音響信号を分析して得られる物理特徴量に基づいて、音像幅を精度よく推定する音像幅推定装置を提供する。【解決手段】  音像幅推定装置１００は、左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割するフィルタバンク７Ｌ、７Ｒと、サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも１つを特徴量とする周波数帯域別特徴量を算出する周波数帯域別特徴量算出手段８ｆと、周波数帯域別特徴量から物理特徴量の代表値を算出する物理特徴量代表値算出手段９と、物理特徴量代表値を推定モデル式に適用して音像幅の推定値を算出する音像幅推定値算出手段１０とを備える。【選択図】図１
2010204325.tif 000002
【特許請求の範囲】
【請求項１】
左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した前記物理特徴量を、前記物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して前記音像幅を推定する音像幅推定装置であって、
前記左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割する周波数帯域分割手段と、
前記周波数帯域分割手段によって分割したサブバンド信号から、前記サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つを前記サブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する周波数帯域別特徴量算出手段と、
前記周波数帯域別特徴量算出手段によって算出した周波数帯域別特徴量に基づいて前記物理特徴量を算出する物理特徴量算出手段と、
前記物理特徴量算出手段によって算出した物理特徴量を、前記推定モデル式に適用して前記音像幅の推定値を算出する推定値算出手段と、
を備えたことを特徴とする音像幅推定装置。
【請求項２】
前記物理特徴量算出手段は、前記周波数帯域別特徴量算出手段によって算出したサブバンド信号ごとの周波数帯域別特徴量の平均、重み付き平均、最大値又は中央値の中の何れか一つを、前記物理特徴量として算出することを特徴とする請求項１に記載の音像幅推定装置。
【請求項３】
左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した前記物理特徴量を、前記物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して前記音像幅を推定する音像幅推定装置であって、
前記左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割する周波数帯域分割手段と、
前記周波数帯域分割手段によって分割したサブバンド信号から、前記サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つを前記サブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する周波数帯域別特徴量算出手段と、
前記周波数帯域別特徴量算出手段によって算出した個々の周波数帯域別特徴量を前記物理特徴量として、前記推定モデル式に適用して前記音像幅の推定値を算出する推定値算出手段と、
を備えたことを特徴とする音像幅推定装置。
【請求項４】
前記周波数帯域分割手段は、周波数帯域幅が１／１２オクターブ以下のサブバンド信号に分割することを特徴とする請求項１乃至請求項３に記載の音像幅推定装置。
【請求項５】
更に、前記推定モデル式における前記重み係数を算出する重み係数算出手段を備え、前記重み係数算出手段は、前記物理特徴量を説明変数とし、前記音像幅を目的変数とする回帰分析によって前記重み係数を算出することを特徴とする請求項１乃至請求項４の何れか一項に記載の音像幅推定装置。
【請求項６】
左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した前記物理特徴量を、前記物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して前記音像幅を推定するために、コンピュータを、
前記左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、予め定められた複数の周波数帯域のサブバンド信号に分割する周波数帯域分割手段、
前記周波数帯域分割手段によって分割したサブバンド信号から、前記サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つを前記サブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する周波数帯域別特徴量算出手段、
前記周波数帯域別特徴量算出手段によって算出した周波数帯域別特徴量に基づいて前記物理特徴量を算出する物理特徴量算出手段、
前記物理特徴量算出手段によって算出した物理特徴量を、前記推定モデル式に適用して前記音像幅の推定値を算出する推定値算出手段、
として機能させることを特徴とする音像幅推定プログラム。
【請求項７】
左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した前記物理特徴量を、前記物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して前記音像幅を推定するために、コンピュータを、
前記左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、予め定められた複数の周波数帯域のサブバンド信号に分割する周波数帯域分割手段、
前記周波数帯域分割手段によって分割したサブバンド信号から、前記サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つを前記サブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する周波数帯域別特徴量算出手段、
前記周波数帯域別特徴量算出手段によって算出した個々の周波数帯域別特徴量を前記物理特徴量として、前記推定モデル式に適用して前記音像幅の推定値を算出する推定値算出手段、
として機能させることを特徴とする音像幅推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、左右２チャンネルで採取したデジタル音響信号を分析して得た物理特徴量に基づいて、聴覚特性である音像幅を推定する音像幅推定装置及び音像幅推定プログラムに関する。
【背景技術】
【０００２】
音響が人間に与える心理的効果の大きさは、主観評価によって定量することができる。この主観評価によって得られる心理的効果の大きさを、音響信号を採取・分析して得られた物理特徴量に基づいて客観評価する手法が多く試みられている。
その中で、人間の聴覚特性の一つである音像幅と物理特徴量との関係についても多くの研究が進められている。音響分析の分野では広く知られ、また受け入れられている物理特徴量として、ＩＡＣＣ（interaural cross-correlation；両耳間相互相関度）がある。一般的には、ＩＡＣＣが小さくなると、音像幅は広がると考えられており、種々の周波数帯域におけるＩＡＣＣの分析についての研究が数多くなされている（例えば、非特許文献１参照）。
また、音像幅と、物理特徴量であるＩＴＤ（interaural time differences；両耳間時間差）及びＩＬＤ（interaural level differences；両耳間レベル差）の変動との関係について報告されている（非特許文献２及び非特許文献３参照）。
更に、例えば、特許文献１においては、音響信号からＩＡＣＦ（interaural cross-correlation function；両耳間相互相関関数）の最大振幅であるＩＡＣＣと、この最大振幅の幅ＷＩＡＣＣとに基づいて、見かけの音源の幅（ＡＳＷ）を評価する手法について記載されている（段落００５０参照）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００３−５７１０８号公報
【非特許文献】
【０００４】
【非特許文献１】Masayuki Morimoto and Kazuhiro Iida, “Appropriate frequency bandwidth in measuring interaural cross-correlation as a physical measure of auditory source width”, Acourstical Science and Technology, 日本, 日本音響学会, 2005, Vol.26, No.2, p.179-184
【非特許文献２】Russell Mason and Francis Rumsey, “A comparison of objective measurements for predicting selected subjective spatial attributes”, Audio Engineering Society 112th Convention Paper 5591, Germany, 2002
【非特許文献３】Jens Blauert and Werner Lindemann, “Auditory spaciousness:Some further psychoacoustic analyses”, Journal of Acoustical Society of America, USA, 1986, Vol.80, No.2, p.533-542
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、左右２チャンネルで採取した音響信号を分析して得られる物理的特徴量に基づく従来の客観評価手法では、用いる物理特徴量と主観評価値との相関が必ずしも高くなく、任意の音源から発せられる音像幅を精度よく評価することができなかった。
【０００６】
そこで、本発明はかかる課題に鑑みてなされたものであり、本発明の目的は、左右に２チャンネルからなるデジタル音響信号を分析して得られる物理特徴量に基づいて、音像幅を精度よく推定する音像幅推定装置を提供することである。
【課題を解決するための手段】
【０００７】
前記した目的を達成するために、請求項１に記載の音像幅推定装置は、左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した物理特徴量を、物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して音像幅を推定する音像幅推定装置であって、周波数帯域分割手段と、周波数帯域別特徴量算出手段と、物理特徴量算出手段と、推定値算出手段と、を備えて構成した。
【０００８】
かかる構成によれば、音像幅推定装置は、周波数帯域分割手段によって、左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割する。次に、音像幅推定装置は、周波数帯域別特徴量算出手段によって、周波数帯域分割手段で分割したサブバンド信号から、サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つをサブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する。続いて、音像幅推定装置は、物理特徴量算出手段によって、周波数帯域別特徴量算出手段で算出した周波数帯域別特徴量に基づいて物理特徴量を算出する。そして、音像幅推定装置は、推定値算出手段によって、物理特徴量算出手段で算出した物理特徴量を、推定モデル式に適用して音像幅の推定値を算出する。
これによって、音像幅推定装置は、物理特徴量を用いた音像幅の客観評価を行う。
【０００９】
請求項２に記載の音像幅推定装置は、請求項１に記載の音像幅推定装置において、物理特徴量算出手段は、周波数帯域別特徴量の平均、重み付き平均、最大値又は中央値の中の何れか一つを、物理特徴量として算出する構成とした。
【００１０】
かかる構成によれば、音像幅推定装置は、物理特徴量算出手段によって、周波数帯域別特徴量算出手段で算出したサブバンド信号ごとの周波数帯域別特徴量の平均、重み付き平均、最大値又は中央値の中の何れか一つを物理特徴量として算出する。そして、音像幅推定装置は、音像幅推定値算出手段によって、物理特徴量算出手段で算出した当該物理特徴量に基づいて音像幅の推定値を算出する。
これによって、音像幅推定装置は、周波数帯域別に算出した特徴量を、当該特徴量の種別ごとに一つの値に集約した物理特徴量に基づいて音像幅の推定値を算出する。
【００１１】
請求項３に記載の音像幅推定装置は、左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した物理特徴量を、物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して音像幅を推定する音像幅推定装置であって、周波数帯域分割手段と、周波数帯域別特徴量算出手段と、推定値算出手段と、を備えて構成した。
【００１２】
かかる構成によれば、音像幅推定装置は、周波数帯域分割手段によって、左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割する。次に、音像幅推定装置は、周波数帯域別特徴量算出手段によって、周波数帯域分割手段で分割したサブバンド信号から、サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つをサブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する。そして、音像幅推定装置は、推定値算出手段によって、周波数帯域別特徴量算出手段で算出した個々の周波数帯域別特徴量を物理特徴量として、推定モデル式に適用して音像幅の推定値を算出する。
これによって、音像幅推定装置は、周波数帯域ごとに算出した物理特徴量を用いた音像幅の客観評価を行う。
【００１３】
請求項４に記載の音像幅推定装置は、請求項１乃至請求項３に記載の音像幅推定装置において、周波数帯域分割手段は、周波数帯域幅が１／１２オクターブ以下のサブバンド信号に分割するよう構成した。
【００１４】
かかる構成によれば、音像幅推定装置は、周波数帯域分割手段によって、左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／１２オクターブ以下のサブバンド信号に分割する。続いて、音像幅推定装置は、周波数帯域別特徴量算出手段によって、１／１２オクターブ以下の周波数帯域幅のサブバンド信号ごとに周波数帯域別特徴量を算出する。音像幅推定装置は物理特徴量算出手段によって、周波数帯域別特徴量算出手段で算出した周波数帯域別特徴量に基づいて物理特徴量を算出する。若しくは音像幅推定装置は、周波数帯域別特徴量算出手段で算出した周波数帯域別特徴量を個々の物理特徴量とする。そして、音像幅推定装置は、推定値算出手段によって、物理特徴量算出手段で算出した物理特徴量若しくは周波数帯域別特徴量算出手段で算出した周波数帯域別特徴量を、推定モデル式に適用して音像幅の推定値を算出する。
これによって、音像幅推定装置は、１／１２オクターブ以下に細かく分割した周波数帯域ごとに算出した特徴量に基づいて音像幅の推定値を算出する。
【００１５】
請求項５に記載の音像幅推定装置は、請求項１乃至請求項４の何れか一項に記載の音像幅推定装置において、更に、重み係数算出手段を備える構成とした。
【００１６】
かかる構成によれば、音像幅推定装置は、重み係数算出手段によって、推定モデル式における重み係数を、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つに基づく物理特徴量若しくは両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つの個々の周波数帯域別特徴量を説明変数とし、音像幅を目的変数とする回帰分析によって予め算出しておく。そして、音像幅推定装置は、推定値算出手段によって、物理特徴量算出手段で算出した物理特徴量若しくは周波数帯域別特徴量算出手段で算出した周波数帯域別特徴量と、重み係数算出手段で予め算出しておいた重み係数とを用いて、推定モデル式によって音像幅の推定値を算出する。
これによって、音像幅推定装置は、回帰分析によって定められた重み係数に従って音像幅の推定値を算出する。
【００１７】
請求項６に記載の音像幅推定プログラムは、左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した物理特徴量を、物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して音像幅を推定するために、コンピュータを、周波数帯域分割手段、周波数帯域別特徴量算出手段、物理特徴量算出手段、推定値算出手段、として機能させることとした。
【００１８】
かかる構成によれば、音像幅推定プログラムは、周波数帯域分割手段によって、左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割する。次に、音像幅推定プログラムは、周波数帯域別特徴量算出手段によって、周波数帯域分割手段で分割したサブバンド信号から、サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つをサブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する。続いて、音像幅推定プログラムは、物理特徴量算出手段によって、周波数帯域別特徴量算出手段で算出した周波数帯域別特徴量に基づいて物理特徴量を算出する。そして、音像幅推定プログラムは、推定値算出手段によって、物理特徴量算出手段で算出した物理特徴量を、推定モデル式に適用して音像幅の推定値を算出する。
これによって、音像幅推定プログラムは、物理特徴量を用いた音像幅の客観評価を行う。
【００１９】
請求項７に記載の音像幅推定プログラムは、左右に２チャンネルからなるデジタル音響信号から物理特徴量を算出し、算出した物理特徴量を、物理特徴量と重み係数とからなる音像幅の推定モデル式に適用して音像幅を推定するために、コンピュータを、周波数帯域分割手段、周波数帯域別特徴量算出手段、推定値算出手段、として機能させることとした。
【００２０】
かかる構成によれば、音像幅推定プログラムは、周波数帯域分割手段によって、左右に２チャンネルからなるデジタル音響信号を、左右のチャンネルの音響信号ごとに、周波数帯域幅が１／６オクターブ以下の複数の周波数帯域のサブバンド信号に分割する。次に、音像幅推定プログラムは、周波数帯域別特徴量算出手段によって、周波数帯域分割手段で分割したサブバンド信号から、サブバンド信号ごとに、両耳間相互相関度、両耳間時間差の時間軸方向における標準偏差又は両耳間レベル差の時間軸方向における標準偏差の中の少なくとも一つをサブバンド信号の左右のチャンネルにおける違いを表す周波数帯域別の特徴量である周波数帯域別特徴量を算出する。そして、音像幅推定プログラムは、推定値算出手段によって、周波数帯域別特徴量算出手段で算出した個々の周波数帯域別特徴量を物理特徴量として、推定モデル式に適用して音像幅の推定値を算出する。
これによって、音像幅推定プログラムは、周波数帯域ごとに算出した物理特徴量を用いた音像幅の客観評価を行う。
【発明の効果】
【００２１】
請求項１又は請求項６に記載の発明によれば、周波数帯域幅が１／６オクターブ以下の周波数帯域ごとに算出した特徴量に基づいて音像幅の推定値を算出するため、安定した精度で音像幅の推定を行うことができる。
請求項２に記載の発明によれば、周波数帯域ごとに算出した特徴量を、当該特徴量の種別ごとに一つの値に集約した物理特徴量に基づいて音像幅の推定値を算出するため、推定モデル式における重み係数の個数を増やすことなく、簡便な計算によって音像幅の推定値を算出することができる。
請求項３又は請求項７に記載の発明によれば、周波数帯域幅が１／６オクターブ以下の周波数帯域ごとに算出した特徴量に基づいて音像幅の推定値を算出するため、精度よく音像幅の推定を行うことができる。
請求項４に記載の発明によれば、周波数帯域幅が１／１２オクターブ以下の周波数帯域ごとに算出した特徴量に基づいて音像幅の推定値を算出するため、より安定した精度で音像幅の推定を行うことができる。
請求項５に記載の発明によれば、推定モデル式における重み係数を、主観評価データと、主観評価データに対応する物理特徴量とを用いた回帰分析によって定めるため、精度よく音像幅の推定を行うことができる。
【図面の簡単な説明】
【００２２】
【図１】本発明に係る第１実施形態の音像幅推定装置の構成を示すブロック図である。
【図２】本発明に係る第１実施形態の音像幅推定装置における演算手段の構成を示すブロック図である。
【図３】本発明に係る第１実施形態の音像幅推定装置の処理の流れを示すフロー図である。
【図４】本発明に係る第１実施形態の音像幅推定装置における推定モデル式の重み係数算出処理の流れを示すフロー図である。
【図５】本発明における音像幅の推定モデル式で用いる物理特徴量と音像幅の主観評価データとの間のピアソン相関分析の結果を示すグラフ図であり、（１）、（２）、（３）及び（４）は、それぞれ音源としてバイオリンのＧ線、Ａ線、Ｄ線及びＥ線の開放弦を用いた場合の結果を示す。
【図６】本発明に係る第２実施形態の音像幅推定装置の構成を示すブロック図である。
【図７】本発明に係る第２実施形態の音像幅推定装置における演算手段の構成を示すブロック図である。
【図８】本発明に係る第２実施形態の音像幅推定装置の処理の流れを示すフロー図である。
【図９】本発明に係る第２実施形態の音像幅推定装置における物理特徴量代表値算出用の重み係数算出処理の流れを示すフロー図である。
【図１０】本発明に係る第３実施形態の音像幅推定装置の構成を示すブロック図である。
【図１１】本発明に係る第３実施形態の音像幅推定装置の処理の流れを示すフロー図である。
【図１２】本発明における音像幅の推定モデル式の重み係数を決定するために用いる、主観評価データを採取するための実験装置の構成例を示す模式図である。
【図１３】本発明に係る音像幅推定装置によって算出した音像幅推定値の例を示すグラフ図である。
【発明を実施するための形態】
【００２３】
以下、本発明の実施形態について、適宜図面を参照して説明する。
［第１実施形態］
まず、図１を参照して、本発明における第１実施形態の音像幅推定装置１００の構成について説明する。図１に示すように、音像幅推定装置１００は、ダミーヘッド１と、マイクロフォン２Ｌ及び２Ｒと、ローパスフィルタ３Ｌ及び３Ｒと、ＡＤ変換器４Ｌ及び４Ｒと、演算手段５と、表示手段１４と、を備えて構成されている。また、演算手段５は、メモリ６Ｌ及び６Ｒと、フィルタバンク（周波数帯域分割手段）７Ｌ及び７Ｒと、周波数帯域別物理特徴量算出手段（周波数帯域別特徴量算出手段）８ｆと、物理特徴量代表値算出手段（物理特徴量算出手段）９と、音像幅推定値算出手段（推定値算出手段）１０と、重み係数記憶手段１１と、推定値重み係数算出手段（重み係数算出手段）１２と、主観評価データ記憶手段１３と、を備えて構成されている。
【００２４】
ダミーヘッド１は、試験対象である音源ＳＳから発生する音響をバイノーラル方式で採取するための模擬頭である。ダミーヘッド１の左右両耳の入り口部には、それぞれマイクロフォン２Ｌ及び２Ｒが取り付けられている。
【００２５】
マイクロフォン２Ｌ及び２Ｒは、ダミーヘッド１のそれぞれ左耳及び右耳の入り口部における音源ＳＳから発生する音響を採取する収音手段である。マイクロフォン２Ｌ及び２Ｒで採取されたアナログ音響信号は、それぞれローパスフィルタ３Ｌ及び３Ｒに入力される。
なお、第１実施形態においては、マイクロフォン２Ｌ及び２Ｒは、ダミーヘッド１の左右両耳の入り口部に配置したが、マイクロフォン２Ｌ及び２Ｒをダミーヘッド１の鼓膜部に配置して収音するようにしてもよい。
また、ダミーヘッド１の替わりに、人間の頭部を模した球体を用いてマイクロフォン２Ｌ及び２Ｒを配置するようにしてもよいし、マイクロフォンスタンドを用いた２点マイクロフォンの形態でマイクロフォン２Ｌ及び２Ｒを配置するようにしてもよい。
【００２６】
ローパスフィルタ３Ｌ及び３Ｒは、それぞれマイクロフォン２Ｌ及び２Ｒによって採取されたアナログ音響信号を入力し、入力したアナログ音響信号からサンプリング周波数ｆｓの１／２を超える高周波数成分をＡＤ変換器４Ｌ及び４Ｒによってデジタル化（サンプリング）する前に帯域制限して、折り返し歪みの発生を防止するためのアンチエイリアシングフィルタである。ローパスフィルタ３Ｌ及び３Ｒは、帯域制限したアナログ音響信号を、それぞれＡＤ変換器４Ｌ及び４Ｒに出力する。
なお、人の可聴周波数の上限は２０ｋＨｚであるから、サンプリング周波数ｆｓは、２０ｋＨｚの２倍である４０ｋＨｚ以上とする必要がある。例えば、サンプリング周波数ｆｓ＝４８ｋＨｚとすると、ローパスフィルタ３Ｌ及び３Ｒによって、ｆｓ／２＝２４ｋＨｚを超える周波数成分を帯域制限するようにすればよい。
【００２７】
ＡＤ変換器４Ｌ及び４Ｒは、それぞれローパスフィルタ３Ｌ及び３Ｒによって帯域制限されたアナログ音響信号を入力し、入力したアナログ音響信号を例えば、サンプリング周波数ｆｓ＝４８ｋＨｚでサンプリングしてデジタル信号に変換する。ＡＤ変換器４Ｌ及び４Ｒは、それぞれデジタル信号に変換した左チャンネルの音響信号ｓｌ（ｎ）及び右チャンネルの音響信号ｓｒ（ｎ）（但し、ｎはサンプリングしたデータの番号を示す）を、それぞれ演算手段５のメモリ６Ｌ及び６Ｒに出力する。
【００２８】
演算手段５は、バイノーラル方式で採取され、ＡＤ変換器４Ｌ及び４Ｒによってデジタル化された音響信号ｓｌ（ｎ）及びｓｒ（ｎ）を入力し、入力した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）を数値演算によって分析することにより音像幅推定値（ハットｙ）を算出する分析手段である。演算手段５は、汎用的なコンピュータを用いて実現することができる。
演算手段５は、算出した音像幅推定値（ハットｙ）を表示手段１４に出力する。
なお、演算手段５の詳細については後記する。
【００２９】
表示手段１４は、演算手段５から入力した音像幅推定値（ハットｙ）を、視認可能に表示する液晶ディスプレイなどの表示装置である。
表示手段１４は、演算手段５から所定の時間間隔ごとに出力される音像幅推定値（ハットｙ）の数値を適宜表示する。なお、表示手段１４は、音像幅推定値（ハットｙ）の経時変化が把握しやすいように、グラフ化して表示するようにしてもよい。
【００３０】
音源ＳＳは、人間に音像幅を誘起させる音響を発生する音響発生手段である。試験対象である音源ＳＳとしては、楽器やスピーカなど任意の音源を用いることができ、音源ＳＳは、１個であっても複数個であってもよい。
【００３１】
次に、演算手段５の各部の構成について説明する。
メモリ６Ｌ及び６Ｒは、それぞれＡＤ変換器４Ｌ及び４Ｒから入力した左チャンネルの音響信号ｓｌ（ｎ）及び右チャンネルの音響信号ｓｒ（ｎ）を記憶する記憶手段である。メモリ６Ｌ及び６Ｒに記憶した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）は、それぞれ適宜にフィルタバンク７Ｌ及び７Ｒによって読み出される。
【００３２】
フィルタバンク（周波数帯域分割手段）７Ｌ及び７Ｒは、それぞれ互いに異なる複数の周波数帯域ｆを通過する特性を有するバンドパスフィルタ群から構成される。ここで、ｆは周波数帯域を識別する番号を示し、ｆ＝１，２，…，Ｆである。また、Ｆは２以上の整数である。
フィルタバンク７Ｌ及び７Ｒは、それぞれメモリ６Ｌ及６Ｒに記憶された左チャンネルの音響信号ｓｌ（ｎ）及び右チャンネルの音響信号ｓｒ（ｎ）を読み出し、読み出した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）の複数の周波数帯域ｆの周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）を、各バンドパスフィルタの出力の組として得るものである。すなわち、フィルタバンク７Ｌ及び７Ｒは、音響信号ｓｌ（ｎ）及びｓｒ（ｎ）を複数の周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）に分割する周波数帯域分割手段である。フィルタバンク７Ｌ及び７Ｒは、音響信号ｓｌ（ｎ）及びｓｒ（ｎ）の各周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）を、それぞれの周波数帯域ｆに対応する周波数帯域別物理特徴量算出手段８ｆに出力する。
【００３３】
フィルタバンク７Ｌ及び７Ｒは、例えば、１／６オクターブバンドフィルタなどの等比帯域フィルタ群で構成することができる。好ましくは、周波数帯域幅が１／６オクターブ以下、更に好ましくは１／１２オクターブ以下の狭帯域の１／ｎオクターブバンドフィルタ（ここで、ｎは１以上の整数）を用いることができる。
なお、フィルタ群を構成する各フィルタは、ＦＩＲ（finite impulse response；有限長インパルス応答）フィルタによって構成することができる。
【００３４】
周波数帯域別物理特徴量算出手段（周波数帯域別特徴量算出手段）８ｆ（ｆ＝１，２，…，Ｆ）は、それぞれフィルタバンク７Ｌ及び７Ｒから音響信号の周波数帯域ｆに対応する左右の周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）を入力し、入力した左右の周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）を分析して、周波数帯域ｆごとの３種類の物理特徴量である周波数帯域別物理特徴量（周波数帯域別特徴量）ｘａ（ｆ）、ｘｔ（ｆ）、ｘｌ（ｆ）を算出して物理特徴量代表値算出手段９に出力する。
【００３５】
物理特徴量代表値算出手段（物理特徴量算出手段）９は、Ｆ個の周波数帯域別物理特徴量算出手段８ｆ（ｆ＝１，２，…，Ｆ）からＦ組の周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）、ｘｌ（ｆ）を入力し、入力した周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）、ｘｌ（ｆ）を物理特徴量の種類ごとに、物理特徴量代表値Ｘａ、Ｘｔ、Ｘｌを算出して音像幅推定値算出手段１０又は推定値重み係数算出手段１２に出力する。
なお、音像幅を推定するための推定モデル式における各物理特徴量代表値Ｘａ、Ｘｔ、Ｘｌに対する重み係数Ｃａ、Ｃｔ及びＣｌを算出する場合は、物理特徴量代表値算出手段９は、物理特徴量代表値Ｘａ、Ｘｔ、Ｘｌを推定値重み係数算出手段１２に出力する。また、推定モデル式と重み係数Ｃａ、Ｃｔ及びＣｌとを用いて音像幅を推定する場合は、物理特徴量代表値算出手段９は、物理特徴量代表値Ｘａ、Ｘｔ、Ｘｌを音像幅推定値算出手段１０に出力する。
【００３６】
ここで、図２を参照（適宜図１参照）して、周波数帯域別物理特徴量算出手段８ｆと物理特徴量代表値算出手段９の詳細な構成について説明する。
図２に示すように、周波数帯域別物理特徴量算出手段８ｆは、窓掛け手段２０Ｌｆ及び２０Ｒｆと、ＣＣＣ（interaural cross-correlation coefficient；両耳間相互相関係数）算出手段２１ｆと、レベル算出手段２２Ｌｆ及び２２Ｒｆと、ＩＡＣＣ算出手段２３ｆと、ＩＴＤ算出手段２４ｆと、ＩＬＤ算出手段２５ｆと、ＩＬＤ標準偏差算出手段２６ｆと、ＩＡＣＣ平均算出手段２７ｆと、ＩＴＤ標準偏差算出手段２８ｆと、を備えて構成されている。
また、物理特徴量代表値算出手段９は、ＩＬＤ標準偏差代表値算出手段３０と、ＩＡＣＣ平均代表値算出手段３１と、ＩＴＤ標準偏差代表値算出手段３２と、を備えて構成されている。
【００３７】
窓掛け手段２０Ｌｆ及び２０Ｒｆは、それぞれフィルタバンク７Ｌ及び７Ｒから対応する周波数帯域ｆの周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）を入力し、入力した周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）に時間窓ｗ（ｎ）を掛けて、順次に所定時間長の信号を切り出す手段である。
窓掛け手段２０Ｌｆ及び２０Ｒｆは、切り出した信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）を、それぞれレベル算出手段２２Ｌｆ及び２２Ｒｆに出力するとともに、左右のチャンネルの信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）を、ＣＣＣ算出手段２１ｆに出力する。
【００３８】
ここで、窓掛け手段２０Ｌｆ及び２０Ｒｆによって周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）から切り出される信号のデータ数をＮ（Ｎは１以上の整数）とすると、時間窓ｗ（ｎ）は、式（１）によって表すことができる。
【００３９】
【数１】
2010204325.tif 000003
【００４０】
なお、時間窓ｗ（ｎ）によって切り出す時間長は、例えば、１０（ｍｓ）〜１００（ｍｓ）とすることができる。
ここで、時間長をｔ（ｍｓ）、ＡＤ変換器４Ｌ及び４Ｒにおけるサンプリング周波数をｆｓ（Ｈｚ）とすると、切り出される信号のデータ数Ｎは、Ｎ＝ｔ１０−３ｆｓとなる。
【００４１】
また、窓掛け手段２０Ｌｆ及び２０Ｒｆは、時間領域において、時間窓ｗ（ｎ）によって、それぞれ周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）に対して移動幅ｄ（ｄは１以上の整数）ずつシフトしながら窓掛けして信号列を切り出す。左チャンネル及び右チャンネルの周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）からｋ番目に切り出される信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）は、それぞれ式（２−１）及び式（２−２）のように表すことができる。
【００４２】
【数２】
2010204325.tif 000004
【００４３】
ここで、両耳間時間差をτ（ｍｓ）、サンプリング周波数をｆｓ（Ｈｚ）とすると、移動幅ｄは、ｄ≧τ１０−３ｆｓとすることができる。すなわち、両耳間時間差τ以上に相当するデータ数ずつ時間窓ｗ（ｎ）によって切り出す位置をシフトするようにすることができる。これによって、後段のＣＣＣ算出手段２１ｆやレベル算出手段２２Ｌｆ及び２２Ｒｆなどの各分析手段によって移動幅ｄに相当する時間幅を時間分解能とした移動分析を行うことができる。
【００４４】
レベル算出手段２２Ｌｆ及び２２Ｒｆは、それぞれ対応する周波数帯域ｆの窓掛け手段２０Ｌｆ及び２０Ｒｆから信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）を入力し、入力したｋ番目の信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）における音響エネルギーレベル（以下、レベルと呼ぶ）ｓｌＥｋ（ｆ）及びｓｒＥｋ（ｆ）を、それぞれ式（３−１）及び式（３−２）によって算出して、対応する周波数帯域ｆのＩＬＤ算出手段２５ｆに出力する。
【００４５】
【数３】
2010204325.tif 000005
【００４６】
ＣＣＣ（interaural cross-correlation coefficient；両耳間相互相関係数）算出手段２１ｆは、それぞれ対応する周波数帯域ｆの窓掛け手段２０Ｌｆ及び２０Ｒｆから信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）を入力し、入力したｋ番目の信号列ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）における両耳間相互相関係数ＣＣＣｋ（τ，ｆ）を、式（４）によって算出して、対応する周波数帯域ｆのＩＡＣＣ算出手段２３ｆ及びＩＴＤ算出手段２４ｆに出力する。
【００４７】
【数４】
2010204325.tif 000006
【００４８】
ＩＡＣＣ（absolute maximum value of the interaural cross-correlation coefficient；両耳間相互相関度）算出手段２３ｆは、対応する周波数帯域ｆのＣＣＣ算出手段２１ｆから両耳間相互相関係数ＣＣＣｋ（τ，ｆ）を入力し、入力した両耳間相互相関係数ＣＣＣｋ（τ，ｆ）における最大振幅である両耳間相互相関度ＩＡＣＣｋ（ｆ）を、式（５−１）によって算出して、対応する周波数帯域ｆのＩＡＣＣ平均算出手段２７ｆに出力する。
なお、ＩＡＣＣ算出手段２３ｆは、両耳間相互相関度ＩＡＣＣｋ（ｆ）を、式（５−１）に替えて、式（５−２）によって算出するようにしてもよい。
【００４９】
【数５】
2010204325.tif 000007
【００５０】
ＩＴＤ（interaural time difference；両耳間時間差）算出手段２４ｆは、対応する周波数帯域ｆのＣＣＣ算出手段２１ｆから両耳間相互相関係数ＣＣＣｋ（τ，ｆ）を入力し、式（６−１）によって、入力した両耳間相互相関係数ＣＣＣｋ（τ，ｆ）において最大振幅を与える時間差τを算出し、算出した時間差τを両耳間時間差ＩＴＤｋ（ｆ）として、対応する周波数帯域ｆのＩＴＤ標準偏差算出手段２８ｆに出力する。
なお、ＩＴＤ算出手段２４ｆは、両耳間時間差ＩＴＤｋ（ｆ）を、式（６−１）に替えて、式（６−２）によって算出するようにしてもよい。
【００５１】
【数６】
2010204325.tif 000008
【００５２】
ＩＬＤ（interaural level difference；両耳間レベル差）算出手段２５ｆは、対応する周波数帯域ｆのレベル算出手段２２Ｌｆ及び２２Ｒｆから、レベルｓｌＥｋ（ｆ）及びｓｒＥｋ（ｆ）を入力し、入力したレベルｓｌＥｋ（ｆ）及びｓｒＥｋ（ｆ）から、式（７）によって両耳間レベル差ＩＬＤｋ（ｆ）を算出して、算出した両耳間レベル差ＩＬＤｋ（ｆ）を、対応する周波数帯域ｆのＩＬＤ標準偏差算出手段２６ｆに出力する。
【００５３】
【数７】
2010204325.tif 000009
【００５４】
ＩＡＣＣ平均算出手段２７ｆは、対応する周波数帯域ｆのＩＡＣＣ算出手段２３ｆから両耳間相互相関度ＩＡＣＣｋ（ｆ）を入力し、窓掛け手段２０Ｌｆ及び２０Ｒｆによって切り出されたすべての区間ｋ（ｋ＝１，２，…，Ｔ）における両耳間相互相関度ＩＡＣＣｋ（ｆ）を入力すると、式（８）によって、時間軸方向における両耳間相互相関度ＩＡＣＣｋ（ｆ）の平均を算出し、算出した平均を周波数帯域別物理特徴量の一つであるｘａ（ｆ）としてＩＡＣＣ平均代表値算出手段３１に出力する。
【００５５】
【数８】
2010204325.tif 000010
【００５６】
なお、周波数帯域別物理特徴量ｘａ（ｆ）は、移動幅ｄごとに算出された両耳間相互相関度ＩＡＣＣｋ（ｆ）の単純平均としたが、これに限定されるものではなく、重み付き平均を用いるようにしてもよいし、最大値又は中央値などを用いるようにしてもよい。
【００５７】
ＩＴＤ標準偏差算出手段２８ｆは、対応する周波数帯域ｆのＩＴＤ算出手段２４ｆから両耳間時間差ＩＴＤｋ（ｆ）を入力し、窓掛け手段２０Ｌｆ及び２０Ｒｆによって切り出されたすべての区間ｋ（ｋ＝１，２，…，Ｔ）における両耳間時間差ＩＴＤｋ（ｆ）を入力すると、式（９）によって、時間軸方向における両耳間時間差ＩＴＤｋ（ｆ）の標準偏差を算出し、算出した標準偏差を周波数帯域別物理特徴量の一つであるｘｔ（ｆ）としてＩＴＤ標準偏差代表値算出手段３２に出力する。
【００５８】
【数９】
2010204325.tif 000011
【００５９】
ＩＬＤ標準偏差算出手段２６ｆは、対応する周波数帯域ｆのＩＬＤ算出手段２５ｆから両耳間レベル差ＩＬＤｋ（ｆ）を入力し、２０Ｌｆ及び２０Ｒｆによって切り出されたすべての区間ｋ（ｋ＝１，２，…，Ｔ）における両耳間レベル差ＩＬＤｋ（ｆ）を入力すると、式（１０）によって、時間軸方向における両耳間レベル差ＩＬＤｋ（ｆ）の標準偏差を算出し、算出した標準偏差を周波数帯域別物理特徴量の一つであるｘｌ（ｆ）としてＩＬＤ標準偏差代表値算出手段３０に出力する。
【００６０】
【数１０】
2010204325.tif 000012
【００６１】
ＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０は、それぞれＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆから周波数帯域ｆごとに算出された周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）を入力し、それぞれ入力した周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の代表値である物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを算出する。ＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０は、それぞれ算出した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを音像幅推定値算出手段１０又は推定値重み係数算出手段１２に出力する。
【００６２】
前記したように、音像幅を推定するための推定モデル式における各物理特徴量代表値Ｘａ、Ｘｔ及びＸｌに対する重み係数Ｃａ、Ｃｔ及びＣｌを算出する場合は、物理特徴量代表値算出手段９の構成要素であるＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０は、それぞれ物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを推定値重み係数算出手段１２に出力する。また、推定モデル式と重み係数Ｃａ、Ｃｔ及びＣｌとを用いて音像幅を推定する場合は、ＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０は、それぞれ物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを音像幅推定値算出手段１０に出力する。
【００６３】
ここで、第１実施形態におけるＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０は、それぞれ周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の代表値として、式（１１−１）、式（１１−２）及び式（１１−３）によって、周波数帯域ｆごとに算出した周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の平均を算出して物理特徴量代表値Ｘａ、Ｘｔ及びＸｌとする。
【００６４】
【数１１】
2010204325.tif 000013
【００６５】
このように、推定モデル式で用いる物理特徴量として、周波数帯域ごとに算出した周波数帯域別物理特徴量を物理特徴量の種別ごとに一つの値に集約した代表値を用いることにより、推定モデル式における重み係数の個数を低減することができ、音像幅推定値（ハットｙ）の算出や重み係数を定めるための主観評価データの採取を簡略化することができる。
【００６６】
図１に戻って（適宜図２参照）、音像幅推定装置１００の構成について説明を続ける。
音像幅推定値算出手段（推定値算出手段）１０は、ＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０から、それぞれ物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを入力するとともに、重み係数記憶手段１１から、予め推定値重み係数算出手段１２によって算出して記憶しておいた重み係数Ｃａ、Ｃｔ及びＣｌを読み出し、式（１２）に示した推定モデル式によって、音像幅の推定値（ハットｙ）を算出して、算出した推定値（ハットｙ）を表示手段１４に出力する。
【００６７】
【数１２】
2010204325.tif 000014
【００６８】
式（１２）に示したように、第１実施形態における音像幅推定値（ハットｙ）は、３つの物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを要素とする３次元ベクトルの絶対値として算出することができる。
【００６９】
重み係数記憶手段１１は、推定値重み係数算出手段１２によって算出した式（１２）に示した推定モデル式の重み係数Ｃａ、Ｃｔ及びＣｌを記憶する記憶手段である。重み係数記憶手段１１に記憶した重み係数Ｃａ、Ｃｔ及びＣｌは、音像幅の推定を行う際に、音像幅推定値算出手段１０によって読み出され、音像幅推定値（ハットｙ）の算出に用いられる。
【００７０】
推定値重み係数算出手段１２は、主観評価データ記憶手段１３に予め記憶しておいた主観評価データｙｉを読み出すとともに、ＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０から、それぞれ当該主観評価データｙｉに対応する物理特徴量代表値Ｘａｉ、Ｘｔｉ及びＸｌｉを入力し、入力した主観評価データｙｉと物理特徴量代表値Ｘａｉ、Ｘｔｉ及びＸｌｉとからなる複数組のデータを用いて、式（１２）に示した推定モデル式の重み係数Ｃａ、Ｃｔ及びＣｌを回帰分析の手法である最小二乗法によって算出し、算出した重み係数Ｃａ、Ｃｔ及びＣｌを重み係数記憶手段１１に記憶する。なお、ｉは、個々の主観評価データを識別する番号である。
【００７１】
ここで、主観評価データｙｉに対応する物理特徴量代表値Ｘａｉ、Ｘｔｉ及びＸｌｉとは、当該主観評価データｙｉを得たときの被験者と同じ音場条件で、ダミーヘッド１に取り付けられたマイクロフォン２Ｌ及び２Ｒを用いて音響信号を採取し、前記した各分析手段を用いて最終的にＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０から出力される物理特徴量代表値Ｘａ、Ｘｔ及びＸｌのことである。
【００７２】
次に、第１実施形態における推定値重み係数算出手段１２による重み係数Ｃａ、Ｃｔ及びＣｌの算出手法について説明する。
第１実施形態では、３つの物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを説明変数とし、音像幅ｙを目的変数とする式（１２）に示した推定モデル式において、回帰分析の手法である最小二乗法によって重み係数Ｃａ、Ｃｔ及びＣｌを算出する。すなわち、音像幅の主観評価データｙｉと推定モデル式によって算出される予測値（ハットｙｉ）との組を予め用意しておき、最小二乗法によって重み係数Ｃａ、Ｃｔ及びＣｌを算出する。
【００７３】
式（１３）に示したように、主観評価データｙｉと推定モデル式によって算出される予測値（ハットｙｉ）との差の二乗和Ｊが最小となる重み係数Ｃａ、Ｃｔ及びＣｌを算出する。なお、式（１３）において、Ｓは主観評価データのデータ数である。
【００７４】
【数１３】
2010204325.tif 000015
【００７５】
ここで、計算の簡略化のため、便宜的に目的変数を音像幅ｙｉの二乗とすると、式（１４）に示した二乗和Ｊが最小となる重み係数Ｃａ、Ｃｔ及びＣｌを算出することになる。
【００７６】
【数１４】
2010204325.tif 000016
【００７７】
式（１４）の推定値（ハットｙｉ）に、式（１２）を代入すると、二乗和Ｊは、式（１５）のように表すことができる。
【００７８】
【数１５】
2010204325.tif 000017
【００７９】
ここで、二乗和Ｊが最小となる条件は、式（１６）に示した二乗和Ｊの各重み係数Ｃａ、Ｃｔ及びＣｌによる偏微分が０になることである。
【００８０】
【数１６】
2010204325.tif 000018
【００８１】
これにより、式（１７）に示した連立方程式が得られる。
【００８２】
【数１７】
2010204325.tif 000019
【００８３】
ここで、式（１７）に示した重み係数Ｃａ、Ｃｔ及びＣｌを変数とする連立方程式において、重み係数Ｃａ、Ｃｔ及びＣｌに対する係数を式（１８）のようにａ１１〜ａ３３及びｂ１〜ｂ３として定義する。
【００８４】
【数１８】
2010204325.tif 000020
【００８５】
式（１８）で定義したａ１１〜ａ３３及びｂ１〜ｂ３を用いると、式（１７）は、式（１９−１）のように表すことができる。そして、式（１９−１）は、式（１９−２）のように変形することができる。
【００８６】
【数１９】
2010204325.tif 000021
【００８７】
ここで、ａ１１〜ａ３３及びｂ１〜ｂ３は、式（１７）に示したように、主観評価データｙｉと、当該主観評価データｙｉを得たときの被験者の位置で採取した音響信号を分析して得られる物理特徴量代表値Ｘａｉ、Ｘｔｉ及びＸｌｉとを用いて算出することができる。
【００８８】
推定値重み係数算出手段１２は、以上説明した手順により、重み係数Ｃａ、Ｃｔ及びＣｌを算出することができる。
【００８９】
主観評価データ記憶手段１３は、重み係数Ｃａ、Ｃｔ及びＣｌを算出するための音像幅の主観評価データｙｉを記憶する記憶手段である。主観評価データ記憶手段１３に記憶した主観評価データｙｉは、推定値重み係数算出手段１２によって読み出され、重み係数Ｃａ、Ｃｔ及びＣｌの算出のために用いられる。
【００９０】
以上、音像幅推定装置１００の構成について説明したが、本発明はこれに限定されるものではない。例えば、音像幅推定装置１００の演算手段５は、一般的なコンピュータにプログラムを実行させ、コンピュータ内の演算装置や記憶装置を動作させることにより実現することができる。このプログラム（音像幅推定プログラム）は、通信回線を介して配布することも可能であるし、ＣＤ−ＲＯＭ等の記録媒体に書き込んで配布することも可能である。
【００９１】
次に、図３を参照（適宜図１及び図２参照）して、音像幅推定装置１００の動作について説明する。
図３に示すように、音像幅推定装置１００は、まず、推定値重み係数算出手段１２によって、式（１２）に示した推定モデル式における重み係数Ｃａ、Ｃｔ及びＣｌを算出して、重み係数記憶手段１１に記憶しておく（ステップＳ１０）。既に重み係数Ｃａ、Ｃｔ及びＣｌが重み係数記憶手段１１に記憶されている場合は、この推定モデル式の重み係数算出処理ステップは省略することができる。なお、推定モデル式の重み係数算出処理ステップの詳細については後記する。
【００９２】
次に、音像幅推定装置１００は、ダミーヘッド１に取り付けられたマイクロフォン２Ｌ及び２Ｒによって、試験対象である音源ＳＳから発生する音響をバイノーラル方式で採取し、採取したアナログ音響信号を、ローパスフィルタ３Ｌ及び３Ｒを介しＡＤ変換器４Ｌ及び４Ｒによって、デジタル信号に変換した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）として、メモリ６Ｌ及び６Ｒに記憶する（ステップＳ１１）。
【００９３】
音像幅推定装置１００は、フィルタバンク７Ｌ及び７Ｒによって、ステップＳ１１でメモリ６Ｌ及び６Ｒに記憶した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）を読み出し、複数の周波数帯域ｆの周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）に分割して、対応する周波数帯域ｆの周波数帯域別物理特徴量算出手段８ｆの窓掛け手段２０Ｌｆ及び２０Ｒｆに出力する（ステップＳ１２）。
なお、ここでは、音像幅推定装置１００は、フィルタバンク７Ｌ及び７Ｒとして、それぞれ１／６オクターブバンドフィルタを用いて演算する。
【００９４】
音像幅推定装置１００は、窓掛け手段２０Ｌｆ及び２０Ｒｆによって、ステップＳ１２でフィルタバンク７Ｌ及び７Ｒから入力した、対応する周波数帯域ｆの周波数帯域成分ｓｌ（ｎ，ｆ）及びｓｒ（ｎ，ｆ）に対して時間窓ｗ（ｎ）を掛けて、所定の移動幅ｄずつシフトした位置の、所定の時間長の音響信号ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）を順次に切り出す。
音像幅推定装置１００は、窓掛け手段２０Ｌｆによって切り出した左チャンネルの音響信号ｙｌｋ（ｎ，ｆ）を対応する周波数帯域ｆのレベル算出手段２２Ｌｆ及びＣＣＣ算出手段２１ｆに順次に出力するとともに、窓掛け手段２０Ｒｆによって切り出した右チャンネルの音響信号ｙｒｋ（ｎ，ｆ）を対応する周波数帯域ｆのレベル算出手段２２Ｒｆ及びＣＣＣ算出手段２１ｆに順次に出力する（ステップＳ１３）。
【００９５】
音像幅推定装置１００は、レベル算出手段２２Ｌｆ及び２２Ｒｆによって、ステップＳ１３でそれぞれ窓掛け手段２０Ｌｆ及び２０Ｒｆから順次に入力した所定の時間長の音響信号ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）から、それぞれレベルｓｌＥｋ（ｆ）及びｓｒｋＥ（ｆ）を算出して、各周波数帯域ｆに対応するＩＬＤ算出手段２５ｆに順次に出力する（ステップＳ１４）。
音像幅推定装置１００は、並行して、ＣＣＣ算出手段によって、ステップＳ１３で窓掛け手段２０Ｌｆ及び２０Ｒｆから順次に入力した音響信号ｙｌｋ（ｎ，ｆ）及びｙｒｋ（ｎ，ｆ）から、両耳間相互相関係数ＣＣＣｋ（ｆ）を算出して、各周波数帯域ｆに対応するＩＡＣＣ算出手段２３ｆ及びＩＴＤ算出手段２４ｆに順次に出力する（ステップＳ１４）。
【００９６】
音像幅推定装置１００は、ＩＡＣＣ算出手段２３ｆによって、ステップＳ１４でＣＣＣ算出手段２１ｆから入力した両耳間相互相関係数ＣＣＣｋ（ｆ）から両耳間相互相関度ＩＡＣＣｋ（ｆ）を算出してＩＡＣＣ平均算出手段２７ｆに順次に出力する（ステップＳ１５）。
音像幅推定装置１００は、並行して、ＩＴＤ算出手段２４ｆによって、ステップＳ１４でＣＣＣ算出手段２１ｆから入力した両耳間相互相関係数ＣＣＣｋ（ｆ）から両耳間時間差ＩＴＤｋ（ｆ）を算出してＩＴＤ標準偏差算出手段２８ｆに順次に出力する（ステップＳ１５）。
音像幅推定装置１００は、更に並行して、ＩＬＤ算出手段２５ｆによって、ステップＳ１４でレベル算出手段２２Ｌｆ及び２２Ｒｆから入力したレベルｓｌＥｋ（ｆ）及びｓｒｋＥ（ｆ）から両耳間レベル差ＩＬＤｋ（ｆ）を算出してＩＬＤ標準偏差算出手段２６ｆに順次に出力する（ステップＳ１５）。
【００９７】
音像幅推定装置１００は、ＩＡＣＣ平均算出手段２７ｆによって、ステップＳ１５でＩＡＣＣ算出手段２３ｆから入力した両耳間相互相関度ＩＡＣＣｋ（ｆ）から時間軸方向における平均を算出し、算出した当該平均を周波数帯域別物理特徴量ｘａ（ｆ）としてＩＡＣＣ平均代表値算出手段３１に順次に出力する（ステップＳ１６）。
音像幅推定装置１００は、並行して、ＩＴＤ標準偏差算出手段２８ｆによって、ステップＳ１５でＩＴＤ算出手段２４ｆから入力した両耳間時間差ＩＴＤｋ（ｆ）から時間軸方向における標準偏差を算出し、算出した当該標準偏差を周波数帯域別物理特徴量ｘｔ（ｆ）としてＩＴＤ標準偏差代表値算出手段３２に順次に出力する（ステップＳ１６）。
音像幅推定装置１００は、更に並行して、ＩＬＤ標準偏差算出手段２６ｆによって、ステップＳ１５でＩＬＤ算出手段２５ｆから入力した両耳間レベル差ＩＬＤｋ（ｆ）から時間軸方向における標準偏差を算出し、算出した当該標準偏差を周波数帯域別物理特徴量ｘｌ（ｆ）としてＩＬＤ標準偏差代表値算出手段３０に順次に出力する（ステップＳ１６）。
【００９８】
音像幅推定装置１００は、ＩＡＣＣ平均代表値算出手段３１によって、ステップＳ１６でＩＡＣＣ平均算出手段２７ｆから入力した周波数帯域別物理特徴量ｘａ（ｆ）の平均を算出し、算出した当該平均を物理特徴量代表値Ｘａとして音像幅推定値算出手段１０に出力する（ステップＳ１７）。
音像幅推定装置１００は、並行して、ＩＴＤ標準偏差代表値算出手段３２によって、ステップＳ１６でＩＴＤ標準偏差算出手段２８ｆから入力した周波数帯域別物理特徴量ｘｔ（ｆ）の平均を算出し、算出した当該平均を物理特徴量代表値Ｘｔとして音像幅推定値算出手段１０に出力する（ステップＳ１７）。
音像幅推定装置１００は、更に並行して、ＩＬＤ標準偏差代表値算出手段３０によって、ステップＳ１６でＩＬＤ標準偏差算出手段２６ｆから入力した周波数帯域別物理特徴量ｘｌ（ｆ）の平均を算出し、算出した当該平均を物理特徴量代表値Ｘｌとして音像幅推定値算出手段１０に出力する（ステップＳ１７）。
【００９９】
音像幅推定装置１００は、音像幅推定値算出手段１０によって、ステップＳ１７でＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０から入力した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌと、ステップＳ１０で推定値重み係数算出手段１２によって重み係数記憶手段１１に記憶しておいた重み係数Ｃａ、Ｃｔ及びＣｌとから、式（１２）によって音像幅推定値（ハットｙ）を算出して表示手段１４に出力する（ステップＳ１８）。
【０１００】
音像幅推定装置１００は、表示手段１４によって、ステップＳ１８で音像幅推定値算出手段１０から入力した音像幅推定値（ハットｙ）を視認可能に表示する（ステップＳ１９）。
以上の処理によって、音像幅推定装置１００は、音像幅を推定することができる。
【０１０１】
次に、図４を参照（適宜図１及び図２参照）して、図３に示した推定モデル式の重み係数算出処理ステップ（ステップＳ１０）における音像幅推定装置１００の動作について説明する。
図４に示すように、音像幅推定装置１００は、まず、予め実施した主観評価によって得られた主観評価データｙｉを不図示の入力手段によって入力し、主観評価データ記憶手段１３に記憶する（ステップＳ３０）。
【０１０２】
次に、音像幅推定装置１００は、マイクロフォン２Ｌ及び２Ｒによって、ステップＳ３０で入力した主観評価データｙｉに対応する音響信号をバイノーラル方式で採取し、採取したアナログ音響信号を、ローパスフィルタ３Ｌ及び３Ｒを介しＡＤ変換器４Ｌ及び４Ｒによって、デジタル信号に変換した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）として、メモリ６Ｌ及び６Ｒに記憶する（ステップＳ３１）。
【０１０３】
ステップＳ３２〜ステップＳ３７までの処理は、それぞれ図３に示した処理におけるステップＳ１２〜ステップＳ１７までの処理と同様であるから説明は省略する。
なお、音像幅推定装置１００は、ステップＳ３０で入力した主観評価データｙｉの個数Ｓに対応して、ステップＳ３１〜ステップＳ３７の処理を繰り返し、推定値重み係数算出手段１２に、Ｓ組の主観評価データｙｉと物理特徴量代表値Ｘａｉ、Ｘｔｉ及びＸｌｉとからなるデータを蓄積する。
【０１０４】
音像幅推定装置１００は、推定値重み係数算出手段１２によって、ステップＳ３７でＩＡＣＣ平均代表値算出手段３１、ＩＴＤ標準偏差代表値算出手段３２及びＩＬＤ標準偏差代表値算出手段３０から入力した物理特徴量代表値Ｘａｉ、Ｘｔｉ及びＸｌｉと、ステップＳ３０で不図示の入力手段から入力して主観評価データ記憶手段１３に記憶しておいた主観評価データｙｉとからなるＳ組のデータを用いて、最小二乗法によって重み係数Ｃａ、Ｃｔ及びＣｌを算出して（ステップＳ３８）、算出した重み係数Ｃａ、Ｃｔ及びＣｌを重み係数記憶手段１１に記憶する（ステップＳ３９）。
以上で、音像幅推定装置１００は、推定モデル式の重み係数算出処理を終了する。
【０１０５】
次に、図５を参照（適宜図１及び図２参照）して、式（１２）に示した本発明における音像幅の推定モデル式で用いる物理特徴量と音像幅の主観評価データとの間のピアソン相関分析の結果について説明する。
【０１０６】
図５の（１）〜（４）は、それぞれバイオリンのＧ線、Ａ線、Ｄ線及びＥ線の開放弦の連続音を音源として用いたピアソン相関分析の結果である。図５の（１）〜（４）において、横軸はフィルタバンク７Ｌ及び７Ｒによって分割する周波数帯域幅を示しており、各図中の左から周波数帯域分割無し（１バンド）、１／１オクターブバンド〜１／９６オクターブバンドとした場合について示している。また、縦軸はピアソン相関係数を示している。「◆」、「□」及び「▲」で示したデータは、音像幅の主観評価データと、それぞれ両耳間相互相関度の周波数帯域についての平均である物理特徴量代表値Ｘａ、両耳間時間差の標準偏差の周波数帯域についての平均である物理特徴量代表値Ｘｔ及び両耳間レベル差の標準偏差の周波数帯域についての平均である物理特徴量代表値Ｘｌとの相関係数を示している。
【０１０７】
何れの結果も、周波数帯域幅を狭くするほど相関が高くなることを示しており、特に１／６オクターブバンド以下で高い相関を示し、１／１２オクターブバンド以下の狭帯域とした場合では、相関係数の値は飽和していることがわかる。
この分析結果より、フィルタバンク７Ｌ及び７Ｒによって分割する周波数帯域幅は、好ましくは１／６オクターブバンド以下、より好ましくは１／１２オクターブバンド以下とすることによって、安定した精度で音像幅を予測できることがわかる。
【０１０８】
［第２実施形態］
次に、図６及び図７を参照して、本発明における第２実施形態の音像幅推定装置１００Ａについて説明する。
図６に示すように、第２実施形態の音像幅推定装置１００Ａは、図１に示した第１実施形態の音像幅推定装置１００とは、演算手段５に替えて演算手段５Ａを備えたことが異なる。詳細には、第２実施形態の音像幅推定装置１００Ａは、図１に示した第１実施形態の音像幅推定装置１００とは、物理特徴量代表値算出手段９及び主観評価データ記憶手段１３に替えて、それぞれ物理特徴量代表値算出手段９Ａ及び主観評価データ記憶手段１３Ａを備えたことと、代表値重み係数記憶手段１５及び代表値重み係数算出手段１６を更に備えたことと、が異なる。
【０１０９】
第１実施形態における物理特徴量代表値算出手段９は、物理特徴量代表値Ｘａ、Ｘｔ及びＸｌとして、それぞれ周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の平均を算出したが、第２実施形態における物理特徴量代表値算出手段９Ａは、物理特徴量代表値Ｘａ、Ｘｔ及びＸｌとして、それぞれ周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の重み付き平均を算出するものである。
【０１１０】
図７に示すように、第２実施形態における演算手段５Ａの物理特徴量代表値算出手段９Ａは、図２に示した第１実施形態における演算手段５の物理特徴量代表値算出手段９とは、ＩＬＤ標準偏差代表値算出手段３０、ＩＡＣＣ平均代表値算出手段３１及びＩＴＤ標準偏差代表値算出手段３２に替えて、ＩＬＤ標準偏差代表値算出手段３０Ａ、ＩＡＣＣ平均代表値算出手段３１Ａ及びＩＴＤ標準偏差代表値算出手段３２Ａを備えたことが異なる。
なお、図１及び図２に示した第１実施形態と同じ構成要素については、同じ符号を付して、説明は適宜省略する。
【０１１１】
第２実施形態におけるＩＬＤ標準偏差代表値算出手段３０Ａは、ＩＬＤ標準偏差算出手段２６ｆから周波数帯域別物理特徴量ｘｌ（ｆ）を入力するとともに、重み係数記憶手段１５から重み係数ｃｌ（ｆ）を読み出し、物理特徴量代表値Ｘｌとして、式（２０−３）によって重み付き平均を算出する。ＩＬＤ標準偏差代表値算出手段３０Ａは、算出した物理特徴量代表値Ｘｌを音像幅推定値算出手段１０に出力する。
第２実施形態におけるＩＡＣＣ平均代表値算出手段３１Ａは、ＩＡＣＣ平均算出手段２７ｆから周波数帯域別物理特徴量ｘａ（ｆ）を入力するとともに、重み係数記憶手段１５から重み係数ｃａ（ｆ）を読み出し、物理特徴量代表値Ｘａとして、式（２０−１）によって重み付き平均を算出する。ＩＡＣＣ平均代表値算出手段３１Ａは、算出した物理特徴量代表値Ｘａを音像幅推定値算出手段１０に出力する。
第２実施形態におけるＩＴＤ標準偏差代表値算出手段３２Ａは、ＩＴＤ標準偏差算出手段２８ｆから周波数帯域別物理特徴量ｘｔ（ｆ）を入力するとともに、重み係数記憶手段１５から重み係数ｃｔ（ｆ）を読み出し、物理特徴量代表値Ｘｔとして、式（２０−２）によって重み付き平均を算出する。ＩＴＤ標準偏差代表値算出手段３２Ａは、算出した物理特徴量代表値Ｘｔを音像幅推定値算出手段１０に出力する。
【０１１２】
【数２０】
2010204325.tif 000022
【０１１３】
主観評価データ記憶手段１３Ａは、図２に示した第１実施形態における主観評価データ記憶手段１３と同様に、代表値重み係数算出手段１２によって音像幅推定値（ハットｙ）を算出するための重み係数Ｃａ、Ｃｔ及びＣｌを算出する際に用いられる主観評価データｙｉを記憶する。加えて、主観評価データ記憶手段１３Ａは、代表値重み係数算出手段１６によって物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを算出するための重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出する際に用いられる主観評価データｙｉを記憶する。重み係数Ｃａ、Ｃｔ及びＣｌを算出する際に用いられる主観評価データｙｉと重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出する際に用いられる主観評価データｙｉとは、同じデータを共用するようにしてもよく、それぞれ異なるデータとしてもよい。
これらの主観評価データｙｉは、不図示の入力手段によって入力され、主観評価データ記憶手段１３Ａに記憶される。
【０１１４】
重み係数算出手段１６は、主観評価データ記憶手段１３Ａから主観評価データｙｉを読み出すとともに、ＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆからそれぞれ当該主観評価データｙｉに対応する３種類の周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）を入力し、物理特徴量の種別ごとに、入力した主観評価データｙｉと周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）とからなる複数組のデータを用いて、式（２１−１）、式（２１−２）及び式（２１−３）に示した音像幅の推定モデル式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）として回帰分析の手法である最小二乗法によって算出する。そして、重み係数算出手段１６は、算出した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を重み係数記憶手段１５に記憶する。なお、ｉは、個々の主観評価データを識別する番号である。
【０１１５】
【数２１】
2010204325.tif 000023
【０１１６】
ここで、主観評価データｙｉに対応する周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）とは、当該主観評価データｙｉを得たときの被験者と同じ音場条件で、ダミーヘッド１に取り付けられたマイクロフォン２Ｌ及び２Ｒを用いて音響信号を採取し、前記した各分析手段を用いて最終的にＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆから出力される周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）のことである。
なお、重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）の算出手法は、前記した第１実施形態における重み係数Ｃａ、Ｃｔ及びＣｌの算出手法と同様であるので、説明は省略する。
【０１１７】
また、第２実施形態においては、周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の代表値として重み付き平均を用いるようにしたが、周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の代表値として、式（２２−１）〜式（２２−３）に示したように、周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の物理特徴量の種別ごとの重み付き最大値や、式（２３−１）〜式（２３−３）に示したように、周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）の物理特徴量の種別ごとの中央値を用いるようにしてもよい。
【０１１８】
【数２２】
2010204325.tif 000024
【０１１９】
【数２３】
2010204325.tif 000025
【０１２０】
但し、式（２３−１）〜式（２３−３）において、ｍｅｄｉａｎ（ａ１，ａ２，…，ａＦ）は、（  ）内の要素ａ１，ａ２，…，ａＦの中の中央値を算出する関数である。
また、これらの式（２２−１）〜式（２２−３）及び式（２３−１）〜式（２３−３）における重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）は、前記した重み付き平均の算出手法と同様の手法により定めることができる。
【０１２１】
重み係数記憶手段１５は、代表値重み係数算出手段１６によって算出した、式（２１−１）、式（２１−２）及び式（２１−３）に示した音像幅の推定モデル式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）、すなわち式（２０−１）〜式（２０−３）に示した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌの算出式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を記憶する記憶手段である。重み係数記憶手段１５に記憶した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）は、ＩＡＣＣ平均代表値算出手段３１Ａ、ＩＴＤ標準偏差代表値算出手段３２Ａ及びＩＬＤ標準偏差代表値算出手段３０Ａによって読み出され、それぞれ物理特徴量代表値Ｘａ、Ｘｔ及びＸｌの算出に用いられる。
【０１２２】
第２実施形態における推定値重み係数算出手段１２は、第１実施形態における推定値重み係数算出手段１２と同様の手法で、式（１２）に示した推定モデル式の重み係数Ｃａ、Ｃｔ及びＣｌを算出し、算出した重み係数Ｃａ、Ｃｔ及びＣｌを重み係数記憶手段１１に記憶する。
【０１２３】
なお、このとき、ＩＡＣＣ平均代表値算出手段３１Ａ、ＩＴＤ標準偏差代表値算出手段３２Ａ及びＩＬＤ標準偏差代表値算出手段３０Ａは、それぞれ代表値重み係数算出手段１６によって予め算出された重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を重み係数記憶手段１５から読み出し、読み出した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）と、それぞれＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆから入力した周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）とから、それぞれ物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを算出し、算出した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを推定値重み係数算出手段１２に出力する。
【０１２４】
次に、図８を参照（適宜図６及び図７参照）して、第２実施形態の音像幅推定装置１００Ａの動作について説明する。
図８に示すように、音像幅推定装置１００Ａは、まず、代表値重み係数算出手段１６によって、式（２０−１）〜式（２０−３）に示した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌの算出式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出して、重み係数記憶手段１５に記憶しておく（ステップＳ５０）。既に重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）が重み係数記憶手段１５に記憶されている場合は、この物理特徴量代表値算出用の重み係数算出処理ステップは省略することができる。なお、物理特徴量代表値算出用の重み係数算出処理ステップの詳細については後記する。
【０１２５】
続いて、音像幅推定装置１００Ａは、推定値重み係数算出手段１２によって、式（１２）に示した推定モデル式における重み係数Ｃａ、Ｃｔ及びＣｌを算出して、重み係数記憶手段１１に記憶しておく（ステップＳ５１）。既に重み係数Ｃａ、Ｃｔ及びＣｌが重み係数記憶手段１１に記憶されている場合は、この推定モデル式の重み係数算出処理ステップは省略することができる。なお、本ステップＳ５１は、図３に示した第１実施形態の音像幅推定装置１００の処理におけるステップＳ１０と同様であるので、詳細な説明は省略する。
【０１２６】
ステップＳ５２〜ステップＳ５７の処理は、それぞれ図３に示した第１実施形態の音像幅推定装置１００の処理におけるステップＳ１１〜ステップＳ１６と同様であるので、説明は省略する。
【０１２７】
音像幅推定装置１００Ａは、ＩＡＣＣ平均代表値算出手段３１Ａによって、ステップＳ５７でＩＡＣＣ平均算出手段２７ｆから入力した周波数帯域別物理特徴量ｘａ（ｆ）と、ステップＳ５０で重み係数記憶手段１５に記憶しておいた重み係数ｃａ（ｆ）とから、式（２０−１）によって重み付き平均を算出し、算出した当該重み付き平均を物理特徴量代表値Ｘａとして音像幅推定値算出手段１０に出力する（ステップＳ５８）。
音像幅推定装置１００Ａは、並行して、ＩＴＤ標準偏差代表値算出手段３２Ａによって、ステップＳ５７でＩＴＤ標準偏差算出手段２８ｆから入力した周波数帯域別物理特徴量ｘｔ（ｆ）と、ステップＳ５０で重み係数記憶手段１５に記憶しておいた重み係数ｃｔ（ｆ）とから、式（２０−２）によって重み付き平均を算出し、算出した当該重み付き平均を物理特徴量代表値Ｘｔとして音像幅推定値算出手段１０に出力する（ステップＳ５８）。
音像幅推定装置１００Ａは、更に並行して、ＩＬＤ標準偏差代表値算出手段３０Ａによって、ステップＳ５０でＩＬＤ標準偏差算出手段２６ｆから入力した周波数帯域別物理特徴量ｘｌ（ｆ）と、ステップＳ５０で重み係数記憶手段１５に記憶しておいた重み係数ｃｌ（ｆ）とから、式（２０−３）によって重み付き平均を算出し、算出した当該重み付き平均を物理特徴量代表値Ｘｌとして音像幅推定値算出手段１０に出力する（ステップＳ５８）。
【０１２８】
音像幅推定装置１００Ａは、音像幅推定値算出手段１０によって、ステップＳ５８でＩＡＣＣ平均代表値算出手段３１Ａ、ＩＴＤ標準偏差代表値算出手段３２Ａ及びＩＬＤ標準偏差代表値算出手段３０Ａから入力した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌと、ステップＳ５１で推定値重み係数算出手段１２によって重み係数記憶手段１１に記憶しておいた重み係数Ｃａ、Ｃｔ及びＣｌとから、式（１２）によって音像幅推定値（ハットｙ）を算出して表示手段１４に出力する（ステップＳ５９）。
【０１２９】
音像幅推定装置１００Ａは、表示手段１４によって、ステップＳ５９で音像幅推定値算出手段１０から入力した音像幅推定値（ハットｙ）を視認可能に表示する（ステップＳ６０）。
以上の処理によって、音像幅推定装置１００Ａは、音像幅を推定することができる。
【０１３０】
次に、図９を参照（適宜図６及び図７参照）して、図８に示した物理特徴量代表値算出用の重み係数算出処理ステップ（ステップＳ５０）における音像幅推定装置１００Ａの動作について説明する。
図９に示すように、音像幅推定装置１００Ａは、まず、予め実施した主観評価によって得られた主観評価データｙｉを不図示の入力手段によって入力し、主観評価データ記憶手段１３Ａに記憶する（ステップＳ７０）。
【０１３１】
次に、音像幅推定装置１００Ａは、マイクロフォン２Ｌ及び２Ｒによって、ステップＳ７０で入力した主観評価データｙｉに対応する音響信号をバイノーラル方式で採取し、採取したアナログ音響信号を、ローパスフィルタ３Ｌ及び３Ｒを介しＡＤ変換器４Ｌ及び４Ｒによって、デジタル信号に変換した音響信号ｓｌ（ｎ）及びｓｒ（ｎ）として、メモリ６Ｌ及び６Ｒに記憶する（ステップＳ７１）。
【０１３２】
ステップＳ７２〜ステップＳ７６までの処理は、それぞれ図８に示した処理におけるステップＳ５３〜ステップＳ５７までの処理と同様であるから説明は省略する。
なお、音像幅推定装置１００Ａは、ステップＳ７０で入力した主観評価データｙｉの個数Ｓに対応して、ステップＳ７１〜ステップＳ７６の処理を繰り返し、代表値重み係数算出手段１６に、Ｓ組の主観評価データｙｉと周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）とからなるデータを蓄積する。
【０１３３】
音像幅推定装置１００Ａは、代表値重み係数算出手段１６によって、ステップＳ７６でＩＡＣＣ平均代表値算出手段３１Ａ、ＩＴＤ標準偏差代表値算出手段３２Ａ及びＩＬＤ標準偏差代表値算出手段３０Ａから入力した周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）と、ステップＳ７０で不図示の入力手段から入力して主観評価データ記憶手段１３Ａに記憶しておいた主観評価データｙｉとからなるＳ組のデータを用いて、最小二乗法によって重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出して（ステップＳ７８）、算出した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を重み係数記憶手段１５に記憶する（ステップＳ７９）。
以上で、音像幅推定装置１００Ａは、物理特徴量代表値算出用の重み係数算出処理を終了する。
【０１３４】
［第３実施形態］
次に、図１０を参照して、本発明における第３実施形態の音像幅推定装置１００Ｂについて説明する。
図１０に示すように、第３実施形態の音像幅推定装置１００Ｂは、図１に示した第１実施形態の音像幅推定装置１００とは、演算手段５に替えて演算手段５Ｂを備えたことが異なる。詳細には、第３実施形態の音像幅推定装置１００Ｂは、図１に示した第１実施形態の音像幅推定装置１００とは、物理特徴量代表値算出手段９を備えていないことと、音像幅推定値算出手段１０、重み係数記憶手段１１及び主観評価データ記憶手段１３に替えて、それぞれ音像幅推定値算出手段１０Ｂ、重み係数記憶手段１１Ｂ及び主観評価データ記憶手段１３Ｂを備えたことと、が異なる。
【０１３５】
第１実施形態における音像幅推定値算出手段１０は、物理特徴量代表値算出手段９によって算出した物理特徴量代表値Ｘａ、Ｘｔ及びＸｌを用いた推定モデル式によって音像幅推定値（ハットｙ）を算出するのに対して、第３実施形態における音像幅推定値算出手段１０Ｂは、周波数帯域別物理特徴量算出手段８ｆによって算出した個々の周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）を物理特徴量として扱い、これらの周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）を用いた推定モデル式によって音像幅推定値（ハットｙ）を算出するものである。
なお、図１に示した第１実施形態と同じ構成要素については、同じ符号を付して、説明は適宜省略する。
【０１３６】
第３実施形態における音像幅推定値算出手段１０Ｂは、周波数帯域別物理特徴量算出手段８ｆのＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆ（図２参照）から、それぞれ周波数帯域別の両耳間相互相関度の時間軸方向の平均である周波数帯域別物理特徴量ｘａ（ｆ）、周波数帯域別の両耳間時間差の時間軸方向の標準偏差である周波数帯域別物理特徴量ｘｔ（ｆ）及び周波数帯域別の両耳間レベル差の時間軸方向の標準偏差である周波数帯域別物理特徴量ｘｌ（ｆ）を入力するとともに、重み係数記憶手段１１Ｂから重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を読み出し、式（２４）に示した推定モデル式によって音像幅推定値（ハットｙ）を算出する。音像幅推定値算出手段１０Ｂは、算出した音像幅推定値（ハットｙ）を表示手段１４に出力する。
【０１３７】
【数２４】
2010204325.tif 000026
【０１３８】
なお、推定モデル式は、式（２４）に限定されるものではなく、例えば、式（２５）のように周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）によって表される他の推定モデル式を用いるようにしてもよい。
【０１３９】
【数２５】
2010204325.tif 000027
【０１４０】
重み係数記憶手段１１Ｂは、推定値重み係数算出手段１２Ｂによって算出した式（２４）に示した推定モデル式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を記憶する記憶手段である。重み係数記憶手段１１Ｂに記憶した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）は、音像幅の推定を行う際に、音像幅推定値算出手段１０Ｂによって読み出され、音像幅推定値（ハットｙ）の算出に用いられる。
【０１４１】
推定値重み係数算出手段１２Ｂは、主観評価データ記憶手段１３Ｂから主観評価データｙｉを読み出すとともに、周波数帯域別物理特徴量算出手段８ｆのＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆ（図２参照）からそれぞれ当該主観評価データｙｉに対応する３種類の周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）を入力し、入力した主観評価データｙｉと周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）とからなる複数組のデータを用いて、式（２４）に示した音像幅の推定モデル式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）として回帰分析の手法である最小二乗法によって算出する。そして、重み係数算出手段１２Ｂは、算出した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を重み係数記憶手段１１Ｂに記憶する。なお、ｉは、個々の主観評価データを識別する番号である。
【０１４２】
ここで、主観評価データｙｉに対応する周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）とは、当該主観評価データｙｉを得たときの被験者と同じ音場条件で、ダミーヘッド１に取り付けられたマイクロフォン２Ｌ及び２Ｒを用いて音響信号を採取し、前記した各分析手段を用いて最終的にＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆ（図２参照）から出力される周波数帯域別物理特徴量ｘａｉ（ｆ）、ｘｔｉ（ｆ）及びｘｌｉ（ｆ）のことである。
なお、重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）の算出手法は、前記した第１実施形態における重み係数Ｃａ、Ｃｔ及びＣｌの算出手法と同様であるので、説明は省略する。
【０１４３】
第３実施形態における推定値重み係数算出手段１２Ｂは、第１実施形態における推定値重み係数算出手段１２と同様の手法で、式（２４）に示した推定モデル式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出し、算出した重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を重み係数記憶手段１１Ｂに記憶する。
【０１４４】
主観評価データ記憶手段１３Ｂは、推定値重み係数算出手段１２Ｂによって音像幅推定値（ハットｙ）を算出するための重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出する際に用いられる主観評価データｙｉを記憶する。主観評価データｙｉは、不図示の入力手段によって入力され、主観評価データ記憶手段１３Ｂに記憶される。
【０１４５】
次に、図１１を参照（適宜図１０参照）して、第３実施形態の音像幅推定装置１００Ｂの動作について説明する。
図１１に示すように、音像幅推定装置１００Ｂは、まず、推定値重み係数算出手段１２Ｂによって、式（２４）に示した音像幅推定値（ハットｙ）の推定モデル式の重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出して、重み係数記憶手段１１Ｂに記憶しておく（ステップＳ９０）。既に重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）が重み係数記憶手段１１Ｂに記憶されている場合は、この推定モデル式の重み係数算出処理ステップは省略することができる。なお、推定モデル式の重み係数算出処理ステップは、図３に示した第１実施形態における推定モデル式の重み係数算出処理ステップとは、説明変数として物理特徴量代表値Ｘａ、Ｘｔ及びＸｌに替えて周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）を用い、重み係数Ｃａ、Ｃｔ及びＣｌに替えて重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）を算出すること以外は、同様であるので、詳細な説明は省略する。
【０１４６】
ステップＳ９１〜ステップＳ９６の処理は、それぞれ図３に示した第１実施形態の音像幅推定装置１００の処理におけるステップＳ１１〜ステップＳ１６と同様であるので、説明は省略する。
【０１４７】
音像幅推定装置１００Ｂは、音像幅推定値算出手段１０Ｂによって、ステップＳ９６で周波数帯域別物理特徴量算出手段８ｆのＩＡＣＣ平均算出手段２７ｆ、ＩＴＤ標準偏差算出手段２８ｆ及びＩＬＤ標準偏差算出手段２６ｆ（図２参照）からそれぞれ入力した周波数帯域別物理特徴量ｘａ（ｆ）、ｘｔ（ｆ）及びｘｌ（ｆ）と、ステップＳ９０で重み係数記憶手段１１Ｂに記憶しておいた重み係数ｃａ（ｆ）、ｃｔ（ｆ）及びｃｌ（ｆ）とから、式（２４）によって音像幅推定値（ハットｙ）を算出し、算出した音像幅推定値（ハットｙ）を表示手段１４に出力する（ステップＳ９７）。
【０１４８】
音像幅推定装置１００Ｂは、表示手段１４によって、ステップＳ９７で音像幅推定値算出手段１０Ｂから入力した音像幅推定値（ハットｙ）を視認可能に表示する（ステップＳ９８）。
以上の処理によって、音像幅推定装置１００Ｂは、音像幅を推定することができる。
【実施例】
【０１４９】
次に、本発明の実施例について説明する。
図１及び図２に示した音像幅推定装置１００において、音源ＳＳとして、バイオリンの各開放弦による連続音を録音して用いた。主観評価は、図１２に示すような実験装置１１０を用いて、被験者ＳＵＢを中心とした被験者ＳＵＢの前方側の半円ＳＣ上に適宜スピーカＳＳ１〜ＳＳ３を配置して、録音しておいたバイオリンの連続音を再生して行った。このとき、音源であるスピーカＳＳ１〜ＳＳ３の配置個数及び被験者ＳＵＢを中心とするスピーカＳＳ１〜ＳＳ３の配置角度θＳＳを調整することにより、被験者ＳＵＢが様々な音像幅を知覚できるように制御した。
音像幅の主観評価値は、被験者ＳＵＢの頭部の中心を視点とする水平方向の角度θに変換した。
【０１５０】
次に、前記した主観評価と同じ音場条件で、被験者ＳＵＢが評価したときと同じ位置にダミーヘッド１を配置し、マイクロフォン２Ｌ及び２Ｒを用いてバイノーラル方式で音響信号を採取した。採取した音響信号に対して、図１及び図２に示した音像幅推定装置１００を用いて、音像幅の推定値（ハットｙ）を算出した。このとき、下限周波数を１５０Ｈｚ、上限周波数を１２ｋＨｚとする１／２４オクターブフィルタによって構成されるフィルタバンク７Ｌ及び７Ｒを用いて周波数帯域の分割を行った。
【０１５１】
本実施例においては、２０種類の音刺激に対する主観評価を行った。その結果を図１３に示す。図１３においては、横軸に音像幅の推定値（ハットｙ）、縦軸に音像幅の主観評価値をとり、結果を示した。図１３に示すように、本発明によって、従来技術による手法に比べて良好に音像幅の推定を行うことができる。
なお、本実施例では、音像幅の推定モデル式として、物理特徴量の線形関数を用いたが、これに限定されるものではなく、物理特徴量の二次関数、べき関数、指数関数などを用いるようにすることもできる。
【符号の説明】
【０１５２】
１      ダミーヘッド
２Ｌ、２Ｒ  マイクロフォン
３Ｌ、３Ｒ  ローパスフィルタ
４Ｌ、４Ｒ  ＡＤ変換器
５、５Ａ、５Ｂ  演算手段
６Ｌ、６Ｒ  メモリ
７Ｌ、７Ｒ  フィルタバンク（周波数帯域分割手段）
８ｆ    周波数帯域別物理特徴量算出手段（周波数帯域別特徴量算出手段）
９、９Ａ  物理特徴量代表値算出手段（物理特徴量算出手段）
１０、１０Ｂ  音像幅推定値算出手段（推定値算出手段）
１１、１１Ｂ  重み係数記憶手段
１２、１２Ｂ  推定値重み係数算出手段（重み係数算出手段）
１３、１３Ａ、１３Ｂ  主観評価データ記憶手段
１４    表示手段
１５    重み係数記憶手段
１６    代表値重み係数算出手段
２０Ｌｆ、２０Ｒｆ  窓掛け手段
２１ｆ  ＣＣＣ算出手段
２２Ｌｆ、２２Ｒｆ  レベル算出手段
２３ｆ  ＩＡＣＣ算出手段
２４ｆ  ＩＴＤ算出手段
２５ｆ  ＩＬＤ算出手段
２６ｆ  ＩＬＤ標準偏差算出手段
２７ｆ  ＩＡＣＣ平均算出手段
２８ｆ  ＩＴＤ標準偏差算出手段
３０、３０Ａ  ＩＬＤ標準偏差代表値算出手段
３１、３１Ａ  ＩＡＣＣ平均代表値算出手段
３２、３２Ａ  ＩＴＤ標準偏差代表値算出手段
１００、１００Ａ、１００Ｂ  音像幅推定装置
ＳＳ、ＳＳ１〜ＳＳ３  音源
ＳＵＢ  被験者
【図１】
2010204325.tif 000028
【図２】
2010204325.tif 000029
【図３】
2010204325.tif 000030
【図４】
2010204325.tif 000031
【図５】
2010204325.tif 000032
【図６】
2010204325.tif 000033
【図７】
2010204325.tif 000034
【図８】
2010204325.tif 000035
【図９】
2010204325.tif 000036
【図１０】
2010204325.tif 000037
【図１１】
2010204325.tif 000038
【図１２】
2010204325.tif 000039
【図１３】
2010204325.tif 000040