音源位置推定装置、音源位置推定方法、及び音源位置推定プログラム
【課題】音声信号の入力と同時に音源位置を実時間で推定する。
【解決手段】信号入力部は複数のチャネルの音声信号を入力し、時間差算出部はチャネル間の音声信号の時間差を算出する時間差算出部と、状態予測部は音源位置を含む音源状態情報である過去の音源状態情報から現在の前記音源状態情報を予測し、状態更新部は前記時間差算出部が算出した時間差と前記状態予測部が予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する。
【解決手段】信号入力部は複数のチャネルの音声信号を入力し、時間差算出部はチャネル間の音声信号の時間差を算出する時間差算出部と、状態予測部は音源位置を含む音源状態情報である過去の音源状態情報から現在の前記音源状態情報を予測し、状態更新部は前記時間差算出部が算出した時間差と前記状態予測部が予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音源位置推定装置、音源位置推定方法、及び音源位置推定プログラムに関する。
【背景技術】
【０００２】
従来から、音源の方向を推定する音源定位技術が提案されている。音源定位技術は、ロボットが周囲の環境を把握し、又は雑音への耐性を強化するために有用である。音源定位技術では、複数のマイクロホンからなるマイクロホンアレイを用い、チャネル間の音波の到来時刻の差を検出し、マイクロホンの配置に基づいて音源の方向を推定する。そのため、各マイクロホンの位置、もしくは音源と各マイクロホン間の伝達関数が既知であること、チャネル間で音声信号を同期収録すること、いずれも必要である。
【０００３】
そこで、非特許文献１に記載の音源定位技術では、空間的に分散配置した複数のマイクロホンを用い、チャネル間で非同期に音源から音声信号を記録する。当該音源定位技術では、記録を終えた音声信号を用いて音源位置及びマイクロホン位置を推定する。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Ｎ．Ｏｎｏ，Ｈ．Ｋｏｈｎｏ，Ｎ．Ｉｔｏ，ａｎｄ Ｓ．Ｓａｇａｙａｍａ、ＢＬＩＮＤ ＡＬＩＧＮＭＥＮＴ ＯＦ ＡＳＹＮＣＨＲＯＮＯＵＳＬＹ ＲＥＣＯＲＤＥＤ ＳＩＧＮＡＬＳ ＦＯＲ ＤＩＳＴＲＩＢＵＴＥＤ ＭＩＣＲＯＰＨＯＮＥ ＡＲＲＡＹ、「２００９ ＩＥＥＥ Ｗｏｒｋｓｈｏｐ ｏｎ Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ ｔｏ Ａｕｄｉｏ ａｎｄ Ａｃｏｕｓｔｉｃｓ」、ＩＥＥＥ、２００９年１０月１８日、ｐｐ．１６１−１６４
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、非特許文献１に記載の音源定位技術では、音声信号の入力と同時に、音源位置を実時間で推定することができない。
【０００６】
本発明は上記の点に鑑みてなされたものであり、音声信号の入力と同時に音源位置を実時間で推定することができる音源位置推定装置、音源位置推定方法、及び音源位置推定プログラムを提供する。
【課題を解決するための手段】
【０００７】
（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、複数のチャネルの音声信号を入力する信号入力部と、チャネル間の音声信号の時間差を算出する時間差算出部と、音源位置を含む音源状態情報である過去の音源状態情報から現在の音源状態情報を予測する状態予測部と前記時間差算出部が算出した時間差と前記状態予測部が予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する状態更新部とを備えることを特徴とする音源位置推定装置である。
【０００８】
（２）本発明のその他の態様は、上述の音源位置推定装置であって、前記状態更新部は、前記誤差に基づいてカルマンゲインを算出し、算出したカルマンゲインに前記誤差を乗ずることを特徴とする。
【０００９】
（３）本発明のその他の態様は、上述の音源位置推定装置であって、前記音源状態情報は、前記音声信号を前記信号入力部に供給する収音部の位置を含むことを特徴とする。
【００１０】
（４）本発明のその他の態様は、上述の音源位置推定装置であって、前記収音部の位置の変化に基づいて、前記音源位置の変化が収束したか否かを判断する収束判定部を備えることを特徴とする。
【００１１】
（５）本発明のその他の態様は、上述の音源位置推定装置であって、前記複数のチャネルの入力信号を、予め定めた音源位置の評価点から前記複数のチャネルの各々に対応する収音部の位置までの位相で補償した信号を加算して得られる評価値を最大にする評価点を定め、定めた評価点と前記音源状態情報が表す音源位置までの距離に基づいて前記音源位置の変化が収束したか否かを判断する収束判定部を備えることを特徴とする。
【００１２】
（６）本発明のその他の態様は、上述の音源位置推定装置であって、前記収束判定部は、前記評価点を、遅延和ビームフォーミング法を用いて定め、定めた評価点と上述の状態更新部が推定した音源状態情報が表す音源位置までの距離に基づいて、前記音源位置の変化が収束したか否かを判断することを特徴とする
【００１３】
（７）本発明のその他の態様は、音源位置推定装置における方法において、前記音源位置推定装置が、複数のチャネルの音声信号を入力する過程と、前記音源位置推定装置が、チャネル間の音声信号の時間差を算出する過程と、前記音源位置推定装置が、音源位置を含む音源状態情報である過去の音源状態情報から現在の前記音源状態情報を予測する過程と前記音源位置推定装置が、前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する過程を有することを特徴とする音源位置推定方法である。
【００１４】
（８）本発明のその他の態様は、音源位置推定装置のコンピュータに、複数のチャネルの音声信号を入力する手順、チャネル間の音声信号の時間差を算出する手順、音源位置を含む音源状態情報である過去の音源状態情報を予測する手順、前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する手順を実行させるための音源位置推定プログラムである。
【発明の効果】
【００１５】
上述の（１）、（７）、（８）の態様によれば、音声信号の入力と同時に音源位置を実時間で推定することができる。
上述の（２）の態様によれば、音源位置の推定誤差が低減されるように音源位置を安定して推定することができる。
上述の（３）の態様によれば、音源位置とマイクロホンの位置を同時に推定することができる。
上述の（４）、（５）、（６）の態様によれば、誤差が収束した音源位置を取得することができる。
【図面の簡単な説明】
【００１６】
【図１】本発明の第１の実施形態に係る音源位置推定装置の構成を示す概略図である。
【図２】本実施形態に係る収音部の配置例を表す平面図である。
【図３】本実施形態に係る収音部における音源の観測時刻を表す図である。
【図４】音源状態情報の予測及び更新の概要を表す概念図である。
【図５】音源及び本実施形態に係る収音部の位置関係の一例を表す概念図である。
【図６】長方形運動モデルの一例を表す概念図である。
【図７】円運動モデルの一例を表す概念図である。
【図８】本実施形態に係る音源位置推定処理を表すフローチャートである。
【図９】本発明の第２の実施形態に係る音源位置推定装置の構成を示す概略図である。
【図１０】本実施形態に係る収束判定部の構成を表す概略図である。
【図１１】本実施形態に係る収束判定処理を表すフローチャートである。
【図１２】推定誤差の時間変化の一例を表す図である。
【図１３】推定誤差の時間変化のその他の例を表す図である。
【図１４】観測時間誤差の一例を表す表である。
【図１５】音源定位状況の一例を表す図である。
【図１６】音源定位状況のその他の例を表す図である。
【図１７】音源定位状況のその他の例を表す図である。
【図１８】収束時間の一例を表す図である。
【図１９】推定された音源位置の誤差の一例を表す図である。
【発明を実施するための形態】
【００１７】
（第１の実施形態）
以下、図面を参照しながら本発明の実施形態について説明する。
図１は、本実施形態に係る音源位置推定装置１の構成を示す概略図である。
音源位置推定装置１は、Ｎ個の（Ｎは、１よりも大きい整数）収音部１０１−１〜１０１−Ｎと、信号入力部１０２、時間差算出部１０３、状態推定部１０４、収束判定部１０５、及び位置出力部１０６を含んで構成される。
状態推定部１０４は、状態更新部１０４１及び状態予測部１０４２を含んで構成される。
【００１８】
収音部１０１−１〜１０１−Ｎは、空気の振動である音波を電気信号であるアナログ音声信号に変換する電気音響変換器を備える。収音部１０１−１〜１０１−Ｎは、変換したアナログ音声信号を信号入力部１０２に出力する。
収音部１０１−１〜１０１−Ｎは、例えば、音源位置推定装置１の筐体の外部に分散配置されていてもよい。この場合、収音部１０１−１〜１０１−Ｎは、各々、生成した１チャネルの音声信号を無線又は有線で信号入力部１０２に出力する。収音部１０１−１〜１０１−Ｎの各々は、例えば、マイクロホンユニットである。
【００１９】
ここで、収音部１０１−１〜１０１−Ｎの配置例について説明する。
図２は、本実施形態に係る収音部１０１−１〜１０１−８の配置例を表す平面図である。
図２において、横方向がｘ軸方向であり、縦方向がｙ軸方向である。
図２に示された縦長の長方形は、高さ方向（ｚ軸方向）の座標が一定である受聴室６０１の水平面を表す。図２において、黒丸は、各々収音部１０１−１〜１０１−８が配置されている位置を表す。
収音部１０１−１は、受聴室６０１の中央に配置されている。収音部１０１−２は、受聴室６０１の中央からｘ軸の正方向にずれた位置に配置されている。収音部１０１−３は、収音部１０１−２よりもｙ軸の正方向にずれた位置に配置されている。収音部１０１−４は、収音部１０１−３よりもｘ軸の負方向、ｙ軸の正方向にずれた位置に配置されている。収音部１０１−５は、収音部１０１−４よりもｘ軸の負方向、ｙ軸の負方向にずれた位置に配置されている。収音部１０１−６は、収音部１０１−５よりもｙ軸の負方向にずれた位置に配置されている。収音部１０１−７は、収音部１０１−６よりもｘ軸の正方向、ｙ軸の負方向にずれた位置に配置されている。収音部１０１−８は、収音部１０１−７よりもｘ軸の正方向、ｙ軸の正方向にずれた位置であって、収音部１０１−２よりもｙ軸の正方向にずれた位置に配置されている。このように、収音部１０１−２〜１０１−８は、収音部１０１−１を中心にｘ、ｙ平面上を反時計回りに順に配置されている。
【００２０】
図１に戻り、信号入力部１０２には、収音部１０１−１〜１０１−Ｎの各々からのアナログ音声信号が入力される。以下の説明では、収音部１０１−１〜１０１−Ｎに各々対応するチャネルを、チャネル１〜Ｎと呼ぶ。信号入力部１０２は、各チャネルのアナログ音声信号をアナログディジタル（Ａ／Ｄ、Ａｎａｌｏｇ−ｔｏ−Ｄｉｇｉｔａｌ）変換して、ディジタル音声信号を生成する。
信号入力部１０２は、変換した各チャネルのディジタル音声信号を時間差算出部１０３に出力する。
【００２１】
時間差算出部１０３は、信号入力部１０２から入力された音声信号についてチャネル間の時間差を算出する。時間差算出部１０３は、例えば、チャネル１の音声信号と、チャネルｎ（ｎは、１よりも大きく、Ｎと等しい又はＮより小さい整数）の音声信号との時間差ｔｎ，ｋ−ｔ１，ｋ（以下、Δｔｎ，ｋと表す）を算出する。ここで、ｋは、離散時刻を表す整数である。時間差算出部１０３は、時間差Δｔｎ，ｋを算出する際、例えば、チャネル１の音声信号とチャネルｎの音声信号の間で時間差を与えて、両者間の相互相関を算出し、算出した相互相関が最大となる時間差を選択する。
【００２２】
ここで、時間差Δｔｎ，ｋについて図３を用いて説明する。
図３は、収音部１０１−１及び１０１−ｎにおける音源の観測時刻ｔ１，ｋ、ｔｎ，ｋをそれぞれ表す図である。
図３において、横軸が時刻ｔ、縦軸が収音部を表す。図３において、Ｔｋは、音源が音波を発生させた時刻（発音時刻）を表す。ｔ１，ｋは、収音部１０１−１に音源から受信した音波が観測される時刻（観測時刻）を表す。ｔｎ，ｋは、収音部１０１−ｎに音源から受信した音波が観測される観測時刻を表す。観測時刻ｔ１，ｋは、発音時刻Ｔｋにチャネル１における観測時刻誤差ｍ１τと音源から収音部１０１−１までの音波の伝搬時間Ｄ１，ｋ／ｃが加わった時刻である。観測時刻誤差ｍ１τとは、チャネル１の音声信号が観測される時刻の絶対時刻に対する差である。観測時刻誤差が生じる原因は、主に収音部１０１−ｎの位置と音源の位置の計測誤差や収音部１０１−ｎに音波が到達した到達時刻の観測誤差である。観測時刻誤差Ｄ１，ｋは、音源から収音部１０１−ｎまでの距離である。ｃは、音速である。観測時刻ｔｎ，ｋは、発音時刻Ｔｋにチャネルｎにおける観測時刻誤差ｍｎτと音源から収音部１０１−ｎまでの音波の伝搬時間Ｄｎ，ｋ／ｃが加わった時刻である。従って、時間差Δｔｎ，ｋ（＝ｔｎ，ｋ−ｔ１，ｋ）は式（１）で表される。
【００２３】
【数１】
【００２４】
音源から収音部１０１−ｎまでの距離Ｄｎ，ｋは、式（２）で表される。
【００２５】
【数２】
【００２６】
式（２）において、（ｘｋ，ｙｋ）は時刻ｋにおける音源の位置を表す。（ｍｎｘ，ｍｎｙ）は、収音部１０１−ｎの位置を表す。
ここで、各チャネルｎの時間差Δｔｎ，ｋを要素とする（Ｎ−１）列のベクトル［Δｔ２，ｋ，…，Δｔｎ，ｋ，…，ΔｔＮ，ｋ］Ｔを観測値ベクトルζｋと呼ぶ。ここで、Ｔは、行列又はベクトルの転置（ｔｒａｎｓｐｏｓｅ）を表す。時間差算出部１０３は、観測値ベクトルζｋを表す時間差情報を状態推定部１０４に出力する。
【００２７】
図１に戻り、状態推定部１０４は、過去（例えば、時刻ｋ−１）の音源状態情報から現在（時刻ｋ）の音源状態情報を予測し、時間差算出部１０３から入力された時間差情報が表す時間差に基づいて音源状態情報を推定する。音源状態情報は、例えば、音源の位置（ｘｋ，ｙｋ）、各収音部１０１−ｎの位置（ｍｎｘ，ｍｎｙ）及び観測時刻誤差ｍｎτを表す情報を含む。状態推定部１０４は、音源状態情報を推定する際、時間差算出部１０３から入力された時間差情報が表す時間差と予測した音源状態情報に基づく時間差との間の誤差を減少させるように音源状態情報を更新する。状態推定部１０４は、音源状態情報の予測及び更新において、例えば、拡張カルマンフィルタ（Ｅｘｔｅｎｄｅｄ Ｋａｌｍａｎ Ｆｉｌｔｅｒ；ＥＫＦ）法を用いる。ＥＫＦ法を用いた、予測及び更新については後述する。なお、状態推定部１０４は、拡張カルマンフィルタ法の代わりに、最小二乗平均誤差（Ｍｉｎｉｍｕｍ Ｍｅａｎ Ｓｑｕｒｅｄ Ｅｒｒｏｒ；ＭＭＳＥ）法、その他の方式を用いてもよい。
状態推定部１０４は、推定した音源状態情報を収束判定部１０５に出力する。
【００２８】
収束判定部１０５は、状態推定部１０４から入力された音源状態情報ηｋ’が表す音源位置の変化が収束したか否か判断する。収束判定部１０５は、音源の推定位置が収束したことを表す音源収束情報を位置出力部１０６に出力する。ここで、記号’は、推定値であることを表す記号である。
収束判定部１０５は、例えば、過去の収音部１０１−ｎの推定位置（ｍｎｘ，ｋ−１’，ｍｎｙ，ｋ−１’）と現在の収音部１０１−ｎの推定位置（ｍｎｘ，ｋ’，ｍｎｙ，ｋ’）の間の平均距離Δηｍ’を算出する。収束判定部１０５は、平均距離Δηｍ’が予め設定された閾値よりも小さくなったとき収束したと判断する。このように、音源の推定位置を直接収束判断に用いないのは、音源位置は未知であり時間経過によって変化するためである。反面、収音部１０１−ｎの推定位置（ｍｎｘ，ｋ’，ｍｎｙ，ｋ’）を収束判断に用いるのは、収音部１０１−ｎの位置は固定であり、音源状態情報が、音源の推定位置の他、収音部１０１−ｎの推定位置にも依存するためである。
【００２９】
位置出力部１０６は、収束判定部１０５から音源収束情報が入力された場合、収束判定部１０５から入力された音源状態情報に含まれる音源位置情報を外部に出力する。
【００３０】
次に、ＥＫＦ法を用いた音源状態情報の予測及び更新の概要を説明する。
図４は、音源状態情報の予測及び更新の概要を表す概念図である。
図４において、黒塗りの星印は、音源位置の真値を表す。白抜きの星印は、音源位置の推定値を表す。黒丸は、それぞれ収音部１０１−１、１０１−ｎの位置の真値を表す。白丸は、それぞれ収音部１０１−１、１０１−ｎの位置の推定値を示す。収音部１０１−ｎの位置を中心とする実線の円４０１は、収音部１０１−ｎの位置の観測誤差の大きさを表す。収音部１０１−ｎの位置を中心とする一点鎖線の円４０２は、後述する更新ステップを経た後の収音部１０１−ｎの位置の観測誤差の大きさを表す。即ち、円４０１及び４０２は、更新ステップでは、観測誤差が低減されるように収音部１０１−ｎの位置を含む音源状態情報が更新されることを表す。観測誤差は、後述する分散共分散行列Ｐｋ’で定量的に表される。音源の位置を中心とする破線の円４０３は、現実の音源の位置と、音源の移動モデルを用いて推定される音源の位置との間のモデル誤差Ｒを表す円である。モデル誤差は、後述する分散共分散行列Ｒで定量的に表される。
【００３１】
ＥＫＦ法は、Ｉ．観測ステップ、ＩＩ．更新ステップ、ＩＩＩ．予測ステップを含む。状態推定部１０４は、これらのステップを繰り返して実行する。
Ｉ．観測ステップでは、状態推定部１０４は、時間差算出部１０３から時間差情報を入力される。状態推定部１０４は、音源からの音声信号に対する収音部１０１−１、１０１−ｎ間の時間差ΔＴｎ，ｋを表す時間差情報ζｋが観測値として入力される。
ＩＩ．更新ステップでは、状態推定部１０４は、観測値ベクトルζｋと音源状態情報ηｋ’に基づく観測値ベクトルζｋ’との観測誤差が低減されるように、音源状態情報の誤差を表す分散共分散行列Ｐｋ’と音源状態情報ηｋ’を更新する。
ＩＩＩ．予測ステップでは、状態予測部１０４２は、真の音源の位置の時間変化を表す運動モデルに基づき、前時刻ｋ−１の音源状態情報ηｋ−１’から現時刻ｋの音源状態情報ηｋ｜Ｋ−１’を予測する。状態予測部１０４２は、前時刻ｋ−１の分散共分散行列Ｐｋ−１’に音源位置の運動モデルと推定位置とのモデル誤差を表す分散共分散行列Ｒに基づいて分散共分散行列Ｐｋ−１’を更新する。
【００３２】
ここで、音源状態情報ηｋ’は、例えば、音源の推定位置（ｘｋ’，ｙｋ’）、収音部１０１−１〜１０１−Ｎの推定位置（ｍ１ｘ，ｋ’，ｍ１ｙ，ｋ’）〜（ｍＮｘ，ｋ’，ｍＮｙ，ｋ’）及び観測時刻誤差の推定値ｍ１τ’〜ｍＮτ’を要素として含む。つまり、音源状態情報ηｋ’は、例えば、ベクトル［ｘｋ’，ｙｋ’，ｍ１ｘ，ｋ’，ｍ１ｙ，ｋ’，ｍ１τ’，…，ｍＮｘ，ｋ’，ｍＮｙ，ｋ’，ｍＮτ’］Ｔで表わされる情報である。このように、ＥＫＦ法を用いることで、予測誤差が徐々に低減されるように、未知である音源位置、収音部１０１−１〜１０１−Ｎの位置及び観測時刻誤差が予測される。
【００３３】
図１に戻り、状態推定部１０４の構成について説明する。
状態推定部１０４は、状態更新部１０４１と状態予測部１０４２とを含んで構成される。
状態更新部１０４１は、時間差算出部１０３から観測値ベクトルζｋを表す時間差情報が入力される（Ｉ．観測ステップ）。状態更新部１０４１は、状態予測部１０４２から入力された音源状態情報ηｋ｜ｋ−１’と共分散行列Ｐｋ｜ｋ−１が入力される。音源状態情報ηｋ｜ｋ−１’は、前時刻ｋ−１の音源状態情報ηｋ−１’から予測された現時刻ｋの音源状態情報を表す。共分散行列Ｐｋ｜ｋ−１の各要素は、音源状態情報ηｋ｜ｋ−１’が表すベクトルにおける各要素間の共分散である。即ち、この共分散行列Ｐｋ｜ｋ−１は、音源状態情報ηｋ｜ｋ−１’の誤差を表す。その後、状態更新部１０４１は、音源状態情報ηｋ｜ｋ−１’を時刻ｋのηｋ｜ｋ−１’に更新し、共分散行列Ｐｋ｜ｋ−１を共分散行列Ｐｋに更新する（ＩＩ．更新ステップ）。状態更新部１０４１は、更新した現時刻ｋの音源状態情報ηｋ’及び共分散行列Ｐｋを状態予測部１０４２に出力する。
【００３４】
次に、更新ステップにおける更新処理について、より詳細に説明する。
状態更新部１０４１は、観測値ベクトルζｋに観測誤差ベクトルδｋを加算し、加算して得られた和に観測値ベクトルζｋを更新する。観測誤差ベクトルδｋは、平均値が０であり予め定めた共分散で分布しているガウス分布に従う乱数ベクトルである。この共分散を各行各列の要素として含む行列を共分散行列Ｑと表す。
【００３５】
状態更新部１０４１は、音源状態情報ηｋ｜ｋ−１’、共分散行列Ｐｋ｜ｋ−１及び共分散行列Ｑに基づいて、例えば、式（３）を用いてカルマンゲインＫｋを算出する。
【００３６】
【数３】
【００３７】
式（３）において、行列Ｈｋは、式（４）で表されるように観測関数ベクトルｈ（ηｋ｜ｋ−１’）の各要素を、音源状態情報ηｋ｜ｋ−１’の各要素で偏微分して得られるヤコビアンである。
【００３８】
【数４】
【００３９】
観測関数ベクトルｈ（ηｋ’）は、式（５）で表される。
【００４０】
【数５】
【００４１】
観測関数ベクトルｈ（ηｋ’）は、音源状態情報ηｋ’に基づく観測値ベクトルζｋ’である。そこで、状態更新部１０４１は、例えば式（５）を用いて、前時刻ｋ−１の音源状態情報ηｋ−１’から予測された現時刻ｋの音源状態情報ηｋ｜ｋ−１’に対する観測値ベクトルζｋ｜ｋ−１’を算出する。
次に、状態更新部１０４１は、現時刻ｋの観測値ベクトルζｋ、算出した観測値ベクトルζｋ｜ｋ−１’及び算出したカルマンゲインＫｋに基づいて、例えば式（６）を用いて現時刻ｋの音源状態情報ηｋ’を算出する。
【００４２】
【数６】
【００４３】
即ち、式（６）は、前時刻ｋ−１の観測値ベクトルζｋ’から推定された現時刻ｋの観測値ベクトルζｋ｜ｋ−１’に、残差値を加算して現時刻ｋの音源状態情報ηｋ’を算出することを表す。加算される残差値は、観測された現時刻ｋの観測値ベクトルζＫから観測値ベクトルζｋ｜ｋ−１’の差にカルマンゲインＫｋを乗じて得られるベクトル値である。
次に、状態更新部１０４１は、カルマンゲインＫｋ、行列Ｈｋ、及び前時刻ｋ−１の共分散行列Ｐｋ−１から予測された現時刻ｋの共分散行列Ｐｋ｜ｋ−１に基づき、例えば式（７）を用いて現時刻ｋの共分散行列Ｐｋを算出する。
【００４４】
【数７】
【００４５】
式（７）において、Ｉは単位行列を表す。即ち、式（７）は、単位行列ＩからカルマンゲインＫｋと行列Ｈｋとの積を減じて得られた行列を乗じて、音源状態情報ηｋ’の誤差の大きさを低減することを表す。
【００４６】
状態予測部１０４２は、状態更新部１０４１から現時刻ｋの音源状態情報ηｋ’及び共分散行列Ｐｋが入力される。状態予測部１０４２は、前時刻ｋ−１の音源状態情報ηｋ−１’から現時刻ｋの音源状態情報ηｋ｜ｋ−１’を予測し、共分散行列Ｐｋ−１から共分散行列Ｐｋ｜ｋ−1を予測する（ＩＩＩ．予測ステップ）。
【００４７】
次に、予測ステップにおける予測処理について、より詳細に説明する。
本実施形態では、例えば、前時刻ｋ−１における音源位置（ｘｋ−１’，ｙｋ−１’）が、現時刻ｋでの間に、移動量（Δｘ，Δｙ）Ｔだけずれるという運動モデルを仮定する。
状態予測部１０４２は、移動量（Δｘ，Δｙ）Ｔに、その誤差を表す誤差ベクトルεｋを加算して、加算して得られた和に移動量（Δｘ，Δｙ）Ｔを更新する。誤差ベクトルεｋは、平均値が０でありガウス分布に従う乱数ベクトルである。このガウス分布の特性を表す共分散を各行各列の要素として含む行列を共分散行列Ｒと表す。
状態予測部１０４２は、前時刻ｋ−１の音源状態情報ηｋ−１’から現時刻ｋの音源状態情報ηｋ｜ｋ−１’を、例えば式（８）を用いて予測する。
【００４８】
【数８】
【００４９】
式（８）において、行列Ｆηは、式（９）で表される２行２＋３Ｎ列の行列である。
【００５０】
【数９】
【００５１】
次に、状態予測部１０４２は、前時刻ｋ−１の共分散行列Ｐｋ−１から現時刻ｋの共分散行列Ｐｋ｜ｋ−１を、例えば式（１０）を用いて予測する。
【００５２】
【数１０】
【００５３】
即ち、式（１０）は、移動量の誤差を表す共分散行列Ｒに、前時刻ｋ−１の共分散行列Ｐｋ−１で表される音源状態情報ηｋ−１’の誤差を加算して現時刻ｋの共分散行列Ｐｋを算出することを表す。
【００５４】
状態予測部１０４２は、算出した時刻ｋの音源状態情報ηｋ｜ｌｋｋ−１’と共分散行列Ｐｋ｜ｋ−１を状態更新部１０４１に出力する。状態予測部１０４２は、算出した時刻ｋの音源状態情報ηｋ｜ｋ−１’を収束判定部１０５に出力する。
【００５５】
なお、上述では、状態推定部１０４は、Ｉ．観測ステップ、ＩＩ．更新ステップ、ＩＩＩ．予測ステップを時刻ｋ毎に実行する旨、説明したが本実施形態では、これには限られない。本実施形態では、状態推定部１０４は、Ｉ．観測ステップ及びＩＩ．更新ステップを時刻ｋ毎に実行し、ＩＩＩ．予測ステップを、時刻ｌ（エル）毎に実行してもよい。時刻ｌは、時刻ｋとは異なる時間間隔毎に計数される離散時刻である。例えば、前時刻ｌ−１から現時刻ｌまでの時間間隔は、前時刻ｋ−１から現時刻ｋまでの時間間隔よりも広くてもよい。これにより、状態推定部１０４の動作と時間差算出部１０３の動作タイミングが異なっても、相互の処理を同期させることができる。
そこで、状態更新部１０４１は、状態予測部１０４２が出力した時刻ｌの音源状態情報ηｌ｜ｌ−１’を対応する時刻ｋの音源状態情報ηｋ｜ｋ−１’として入力されるようにする。状態予測部１０４２が出力した共分散行列Ｐｌ｜ｌ−１を、状態更新部１０４１は、共分散行列Ｐｋ｜ｋ−１として入力されるようにする。また、状態予測部１０４２は、状態更新部１０４１が出力した音源状態情報ηｋ’を対応する前時刻ｌ−１の音源状態情報ηｌ−１’として入力されるようにする。状態更新部１０４１が出力した共分散行列Ｐｋを、状態予測部１０４２は共分散行列Ｐｌ−１として入力されるようにする。
【００５６】
次に、音源及び収音部１０１−ｎの位置関係の一例について説明する。
図５は、音源及び収音部１０１−ｎの位置関係の一例を表す概念図である。
図５において、黒塗りの★印は、前時刻ｋ−１の音源位置（ｘｋ−１，ｙｋ−１）及び現時刻ｋの音源位置（ｘｋ，ｙｋ）を表す。音源位置（ｘｋ−１，ｙｋ−１）を起点とし、音源位置（ｘｋ，ｙｋ）を終点とする一点破線で表される矢印は、移動量（Δｘ，Δｙ）Ｔを表す。
黒塗りの●印は、収音部１０１−ｎの位置（ｍｎｘ，ｍｎｙ）Ｔを表す。音源位置（ｘｋ，ｙｋ）Ｔを起点とし、収音部１０１−ｎの位置（ｍｎｘ，ｍｎｙ）Ｔを終点とする実線の近傍に表わされているＤｎ，ｋは、これらの間の距離を表す。本実施形態では収音部１０１−ｎの真の位置は定数であると仮定されているが、収音部１０１−ｎの予測値には誤差が含まれている。そのため、収音部１０１−ｎの予測値は変数である。また、距離Ｄｎ，ｋの誤差に対する指標が共分散行列Ｐｋである。
【００５７】
次に、音源の運動モデルの一例として長方形運動モデルについて説明する。
図６は、長方形運動モデルの一例を表す概念図である。
長方形運動モデルは、音源が長方形の軌道上を運動することを仮定する運動モデルである。図６において、横軸がｘ座標、縦軸がｙ座標を表す。図６に表される長方形は、音源が運動する軌道を表す。この長方形のｘ座標の最大値がｘｍａｘ、最小値がｘｍｉｎである。ｙ座標の最大値がｙｍａｘ、最小値がｙｍｉｎである。音源は、長方形の一辺の上を直進し、長方形の一頂点に到達したとき、つまり音源のｘ座標がｘｍａｘもしくはｘｍｉｎ、ｙ座標がｙｍａｘもしくはｙｍｉｎに到達したとき運動方向を９０°回転する。
即ち、長方形運動モデルでは、音源の移動方向θｓ，ｌ−１は、ｘの正方向を基準として０°、９０°、１８０°、−９０°の何れかである。音源が辺上を運動する場合、運動方向の変化量ｄθｓ，ｌ−１Δｔは、０°である。ここで、ｄθｓ，ｌ−１は、音源の角速度を表し、Δｔは、前時刻ｌ−１から現時刻ｌまでの時間間隔を表す。音源が頂点に到達した場合、運動方向の変化量ｄθｓ，ｌ−１Δｔは、反時計回りを正値として９０°又は−９０°である。
【００５８】
長方形運動モデルを用いる場合、本実施形態では、音源位置情報を、２次元の直交座標（ｘｌ，ｘｌ）と運動方向θを要素とする３次元のベクトルηｓ，ｌで表してもよい。音源位置情報ηｓ，ｌは、音源状態情報ηｌに含まれる情報である。この場合、状態予測部１０４２は、式（８）の代わりに式（１１）を用いて音源位置情報の予測を行ってもよい。
【００５９】
【数１１】
【００６０】
式（１１）において、δηは、移動量の誤差ベクトルである。誤差ベクトルδηは、平均値が０であり予め定めた共分散で分布するガウス分布に従う乱数ベクトルである。この共分散を、各行各列の要素として含む行列を共分散行列Ｒと表す。
【００６１】
状態予測部１０４２は、その後、現時刻ｌの共分散行列Ｐｌ｜ｌ−１を、例えば式（１０）の代わりに式（１２）を用いて予測する。
【００６２】
【数１２】
【００６３】
式（１）において、行列Ｇｌは、式（１３）で示される行列である。
【００６４】
【数１３】
【００６５】
式（１３）において、行列Ｆは、式（１４）で示される行列である。
【００６６】
【数１４】
【００６７】
式（１４）において、Ｉ３×３は、３行３列の単位行列であり、Ｏ３×３は、３行３Ｎ列の零行列である。
【００６８】
次に、音源の運動モデルの一例として円運動モデルについて説明する。
図７は、円運動モデルの一例を表す概念図である。
円運動モデルは、音源が円軌道上を運動することを仮定する運動モデルである。図７において、横軸がｘ座標、縦軸がｙ座標を表す。図７に表される円は、音源が運動する軌道を表す。円運動モデルでは、運動方向の変化量ｄθｓ，ｌ−１Δｔが、一定値Δθであり、音源方向もこれに応じて変化する。
【００６９】
円運動モデルを用いる場合も、音源位置情報を、２次元の直交座標（ｘｌ，ｘｌ）と運動方向θを要素とする３次元のベクトルηｓ，ｌで表してもよい。この場合、状態予測部１０４２は、式（８）の代わりに式（１５）を用いて音源位置情報の予測を行う。
【００７０】
【数１５】
【００７１】
状態予測部１０４２は、現時刻ｌの共分散行列Ｐｌ｜ｌ−１を、式（１２）を用いて予測する。但し、行列Ｇｌとして、式（１３）に表される行列Ｇｌの代わりに、式（１６）に表される行列Ｇｌを用いる。
【００７２】
【数１６】
【００７３】
次に、本実施形態に係る音源位置推定処理について説明する。
図８は、本実施形態に係る音源位置推定処理を表すフローチャートである。
（ステップＳ１０１）音源位置推定装置１は、取り扱う変数の初期値を設定する。例えば、状態推定部１０４は、観測時刻ｋ、予測時刻ｌを、それぞれ０と設定し、音源状態情報ηｋ｜ｋ−１と共分散行列Ｐｋ｜ｋ−１をそれぞれ予め定めた値に設定する。その後、ステップＳ１０２に進む。
（ステップＳ１０２）信号入力部１０２は、収音部１０１−１〜１０１−Ｎからチャネル毎の音声信号が各々入力される。信号入力部１０２は、音声信号の入力を継続するか否か判断する。入力を継続する場合（ステップＳ１０２ Ｙ）、信号入力部１０２は、入力された音声信号をＡ／Ｄ変換して時間差算出部１０３に出力し、その後、ステップＳ１０３に進む。入力を継続しない場合（ステップＳ１０２ Ｎ）、処理を終了する。
【００７４】
（ステップＳ１０３）時間差算出部１０３は、信号入力部１０２から入力された音声信号についてチャネル間の時間差を算出する。時間差算出部１０３は、算出されたチャネル間の時間差を要素とする観測値ベクトルζｋを表す時間差情報を状態更新部１０４１に出力する。その後、ステップＳ１０４に進む。
（ステップＳ１０４）状態更新部１０４１は、予め定めた時間毎に観測時刻ｋを１増加させて観測時刻ｋを更新する。その後、ステップＳ１０５に進む。
【００７５】
（ステップＳ１０５）状態更新部１０４１は、時間差算出部１０３から入力された時間差情報が表す観測値ベクトルζｋに観測誤差ベクトルδｋを加算して観測値ベクトルζｋを更新する。
状態更新部１０４１は、音源状態情報ηｋ｜ｋ−１’、共分散行列Ｐｋ｜ｋ−１及び共分散行列Ｑに基づいて、例えば、式（３）を用いてカルマンゲインＫｋを算出する。
状態更新部１０４１は、例えば式（５）を用いて、現観測時刻ｋの音源状態情報ηｋ｜ｋ−１’に対する観測値ベクトルζｋ｜ｋ−１’を算出する。
状態更新部１０４１は、現観測時刻ｋの観測値ベクトルζｋ、算出した観測値ベクトルζｋ｜ｋ−１’及び算出したカルマンゲインＫｋに基づいて、例えば式（６）を用いて現観測時刻ｋの音源状態情報ηｋ’を算出する。
状態更新部１０４１は、カルマンゲインＫｋ、行列Ｈｋ、及び共分散行列Ｐｋ｜ｋ−１に基づき、例えば式（７）を用いて現観測時刻ｋの共分散行列Ｐｋを算出する。その後、ステップＳ１０６に進む。
【００７６】
（ステップＳ１０６）状態更新部１０４１は、現観測時刻ｋが、予測処理を行う予測時刻ｌに相当するか否か判断する。例えば、観測及び更新ステップＮ回（Ｎは、１又は１よりも大きい整数、例えば、５）毎に予測ステップを１回行う場合、観測時刻ｋのＮに対する剰余が０であるか判断する。現観測時刻ｋが予測時刻ｌと判断された場合（ステップＳ１０７ Ｙ）、ステップＳ１０７に進む。現観測時刻ｋが予測時刻ｌと判断されない場合（ステップＳ１０７ Ｎ）、ステップＳ１０２に進む。
【００７７】
（ステップＳ１０７）状態予測部１０４２は、状態更新部１０４１が出力した算出した現観測時刻ｋの音源状態情報ηｋ’及び共分散行列Ｐｋを、前予測時刻ｌ−１の音源状態情報ηｌ−１’及び共分散行列Ｐｌ−１として入力される。
状態予測部１０４２は、前予測時刻ｌ−１の音源状態情報ηｌ−１’から現予測時刻ｌの音源状態情報ηｌ｜ｌ−１’を、例えば式（８）、（１１）又は（１５）を用いて算出する。
状態予測部１０４２は、前予測時刻ｌ−１の共分散行列Ｐｌ−１から現予測時刻ｌの共分散行列Ｐｌ｜ｌ−１を、例えば式（１０）又は（１２）を用いて算出する。
状態予測部１０４２は、現予測時刻ｌの音源状態情報ηｌ｜ｌ−１’と共分散行列Ｐｌ｜ｌ−１を状態更新部１０４１に出力する。状態予測部１０４２は、算出した現予測時刻ｌの音源状態情報ηｌ｜ｌ−１’を、収束判定部１０５に出力する。その後、ステップＳ１０８に進む。
【００７８】
（ステップＳ１０８）状態更新部１０４１は、現予測時刻ｌに１を加えて予測時刻を更新する。状態更新部１０４１は、状態予測部１０４２が出力した予測時刻ｌの音源状態情報ηｌ｜ｌ−１’、共分散行列Ｐｌ｜ｌ−１を、観測時刻ｋのηｋ｜ｋ−１’、共分散行列Ｐｋ｜ｋ−１として入力される。その後、ステップＳ１０９に進む。
【００７９】
（ステップＳ１０９）収束判定部１０５は、状態推定部１０４から入力された音源状態情報ηｌ’が表す音源位置の変化が収束したか否か判断する。収束判定部１０５は、例えば、過去の収音部１０１−ｎの推定位置と現在の収音部１０１−ｎの推定位置の間の平均距離Δηｍ’が予め設定された閾値よりも小さくなったとき収束したと判断する。音源位置の変化が収束したと判断された場合（ステップＳ１０９ Ｙ）、収束判定部１０５は、入力された音源状態情報ηｌ’を位置出力部１０６に出力する。その後、ステップＳ１１０に進む。音源位置の変化が収束したと判断されなかった場合（ステップＳ１０９ Ｎ）、ステップＳ１０２に進む。
（ステップＳ１１０）位置出力部１０６は、収束判定部１０５から入力された音源状態情報に含まれる音源位置情報を外部に出力する。その後、ステップＳ１０２に進む。
【００８０】
このように、本実施形態は、複数のチャネルの音声信号を入力し、チャネル間の音声信号の時間差を算出し、過去の音源位置を含む音源状態情報から現在の前記音源状態情報を予測する。また、本実施形態は、算出した時間差と予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を更新する。これにより、音声信号の入力と同時に音源位置を推定することができる。
【００８１】
（第２の実施形態）
以下、図面を参照しながら本発明の実施形態について説明する。第１の実施形態と同一の構成又は同一の処理については、同一の番号を付す。
図９は、本実施形態に係る音源位置推定装置２の構成を示す概略図である。
音源位置推定装置２は、Ｎ個の収音部１０１−１〜１０１−Ｎと、信号入力部１０２、時間差算出部１０３、状態推定部１０４、収束判定部２０５、及び位置出力部１０６を含んで構成される。即ち、音源位置推定装置２は、音源位置推定装置１（図１参照）の収束判定部１０５の代わりに収束判定部２０５を備え、信号入力部１０２が入力された音声信号を収束判定部２０５にも出力する点が異なる。その他の構成については、音源位置推定装置１と同様である。
【００８２】
次に、収束判定部２０５の構成について説明する。
図１０は、本実施形態に係る収束判定部２０５の構成を表す概略図である。
収束判定部２０５は、ステアリングベクトル（ｓｔｅｅｒｉｎｇ ｖｅｃｔｏｒ）算出部２０５１、周波数領域変換部２０５２、出力算出部２０５３、評価点選択部２０５４、及び距離判定部２０５５を含んで構成される。この構成により、収束判定部２０５は、遅延和ビームフォーミング（Ｄｅｌａｙ−ａｎｄ−Ｓｕｍ Ｂｅａｍｆｏｒｍｉｎｇ，ＤＳ−ＢＦ）法によって推定された評価点と状態推定部１０４から入力された音源状態情報に含まれる音源位置を比較する。ここで、収束判定部２０５は、評価点と音源位置に基づいて音源状態情報が収束したか否かを判断する。
【００８３】
ステアリングベクトル算出部２０５１は、状態予測部１０４２から入力された音源状態情報ηｌ｜ｌ−１’が表す収音部１０１−ｎの位置（ｍｎｘ’，ｍｎｙ’）から音源位置の候補（以下、評価点と呼ぶ）ξｓ’’までの距離Ｄｎ，ｌを算出する。ステアリングベクトル算出部２０５１は、距離Ｄｎ，ｌを算出する際、例えば式（２）を用いる。但し、ステアリングベクトル算出部２０５１は、式（２）の（ｘｋ，ｙｋ）の代わりに評価点ξｓ’’の座標（ｘ’’，ｙ’’）を代入する。この評価点ξｓ’’は、例えば、予め定められた格子点であって、音源が配置されうる空間（例えば、図２に示す受聴室６０１）に配置された複数の格子点の１つである。
ステアリングベクトル算出部２０５１は、算出した距離Ｄｎ，ｌに基づく伝搬遅延Ｄｎ，ｌ／ｃと推定された観測時刻誤差ｍｎτ’を加算してチャネル毎の推定観測時刻ｔｎ，ｌ’’を算出する。ステアリングベクトル算出部２０５１は、算出した推定時間差ｔｎ，ｌ’’に基づいて、ステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を、例えば式（１７）を用いて周波数ω毎に算出する。
【００８４】
【数１７】
【００８５】
式（１７）において、ξｍ’は、収音部１０１−１〜１０１−Ｎの位置の集合を表す。従って、ステアリングベクトルＷ（η’，ω）の各要素は、対応するチャネルｎ（ｎは、１と等しい又は１よりも大きく、Ｎと等しい又はＮよりも小さい）における音源から各収音部１０１−ｎまでの伝搬によって生じた位相の遅延を与える伝達関数である。ステアリングベクトル算出部２０５１は、算出したステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を出力算出部２０５３に出力する。
【００８６】
周波数領域変換部２０５２は、信号入力部１０２から入力された各チャネルの音声信号ｓｎに対して時間領域から周波数領域に変換して、チャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を生成する。周波数領域変換部２０５２は、周波数領域に変換する方式として、例えば、離散フーリエ変換（Ｄｉｓｃｒｅｔｅ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ、ＤＦＴ）を用いる。周波数領域変換部２０５２は、生成したチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を出力算出部２０５３に出力する。
【００８７】
出力算出部２０５３は、周波数領域変換部２０５２からチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）が入力され、ステアリングベクトル算出部２０５１からステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）が入力される。出力算出部２０５３は、周波数領域信号Ｓｎ，ｌ（ω）を要素とする入力信号ベクトルＳｌ（ω）とステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）の内積Ｐ（ξｓ’’，ξｍ’，ω）を算出する。入力信号ベクトルＳｌ（ω）は［Ｓ１，ｌ（ω），…，Ｓｎ，ｌ（ω），…，ＳＮ，ｌ（ω）］Ｔと表される。出力算出部２０５３は、例えば式（１８）を用いて内積Ｐ（ξｓ’’，ξｍ’，ω）を算出する。
【００８８】
【数１８】
【００８９】
式（１８）において、＊は、ベクトル又は行列の複素共役転置（ｃｏｍｐｌｅｘ ｃｏｎｊｕｇａｔｅ ｔｒａｎｓｐｏｓｅ）を表す。式（１８）によれば、入力信号ベクトルＳｋ（ω）の各チャネル成分の伝搬遅延による位相が補償され、各チャネル成分がチャネル間で同期する。そして、位相が補償された各チャネル成分がチャネル間で加算される。
出力算出部２０５３は、算出した内積Ｐ（ξｓ’’，ξｍ’，ω）を、例えば式（１９）を用いて予め定めた周波数帯域にわたって累積して、帯域出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出する。
【００９０】
【数１９】
【００９１】
式（５）において、最低周波数ωｌ（例えば、２００Ｈｚ）を表し、最高周波数ωｈ（例えば、７ｋＨｚ）を表す。
出力算出部２０５３は、算出した帯域出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を評価点選択部２０５４に出力する。
評価点選択部２０５４は、出力算出部２０５３から入力された帯域出力信号＜Ｐ（ξｓ’’，ξｍ’）＞の絶対値を評価値として最大となる評価点ξｓ’’を選択する。評価点選択部２０５４は、選択した評価点ξｓ’’を距離判定部２０５５に出力する。
距離判定部２０５５は、評価点選択部２０５４から入力された評価点ξｓ’’と状態予測部１０４２から入力された音源状態情報ηｌ｜ｌ−１’が表す音源位置（ｘｌ｜ｌ−１’，ｙｌ｜ｌ−１’）との間の距離が予め定めた閾値、例えば、上述の格子点の間隔よりも小さい場合に収束したと判定する。距離判定部２０５５は、収束したと判定した場合、音源の推定位置が収束したことを表す音源収束情報を位置出力部１０６に出力する。また、距離判定部２０５５は、入力された音源状態情報を位置出力部１０６に出力する。
【００９２】
次に、収束判定部２０５における収束判定処理について説明する。
図１１は、本実施形態に係る収束判定処理を表すフローチャートである。
（ステップＳ２０１）周波数領域変換部２０５２は、信号入力部１０２から入力された各チャネルの音声信号ｓｎに対して時間領域から周波数領域に変換し、チャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を生成する。周波数領域変換部２０５２は、生成したチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を出力算出部２０５３に出力する。その後、ステップＳ２０２に出力する。
【００９３】
（ステップＳ２０２）ステアリングベクトル算出部２０５１は、状態推定部１０４から入力された音源状態情報が表す収音部１０１−ｎの位置（ｍｎｘ’，ｍｎｙ’）から評価点ξｓ’’までの距離Ｄｎ，ｌを算出する。ステアリングベクトル算出部２０５１は、算出した距離Ｄｎ，ｌに基づく伝搬遅延Ｄｎ，ｌ／ｃに推定された観測時刻誤差ｍｎτ’を加算してチャネル毎の推定観測時刻ｔｎ，ｌ’’を算出する。ステアリングベクトル算出部２０５１は、算出した推定時間差ｔｎ，ｌ’’に基づいて、ステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を周波数ω毎に算出する。ステアリングベクトル算出部２０５１は、算出したステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を出力算出部２０５３に出力する。その後、ステップＳ２０３に出力する。
【００９４】
（ステップＳ２０３）出力算出部２０５３は、周波数領域変換部２０５２からチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）が入力され、ステアリングベクトル算出部２０５１からステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）が入力される。出力算出部２０５３は、周波数領域信号Ｓｎ，ｌ（ω）を要素とする入力信号ベクトルＳｌ（ω）とステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）の内積Ｐ（ξｓ’’，ξｍ’，ω）を、例えば式（１８）を用いて算出する。
出力算出部２０５３は、算出した内積Ｐ（ξｓ’’，ξｍ’，ω）を、例えば式（１９）を用いて予め定めた周波数帯域にわたって累積し、出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出する。出力算出部２０５３は、算出した出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を評価点選択部２０５４に出力する。その後、ステップＳ２０４に進む。
【００９５】
（ステップＳ２０４）出力算出部２０５３は、全ての評価点について出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出したか否かを判断する。全ての評価点について算出したと判断された場合（ステップＳ２０４ Ｙ）、ステップＳ２０６に進む。全ての評価点について算出していないと判断された場合（ステップＳ２０４ Ｎ）、ステップＳ２０５に進む。
【００９６】
（ステップＳ２０５）出力算出部２０５３は、出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出する評価点を、出力信号を算出していない他の評価点に変更する。その後、ステップＳ２０２に進む。
【００９７】
（ステップＳ２０６）評価点選択部２０５４は、出力算出部２０５３から入力された出力信号＜Ｐ（ξｓ’’，ξｍ’）＞の絶対値を評価値として最大となる評価点ξｓ’’を選択する。評価点選択部２０５４は、選択した評価点ξｓ’’を距離判定部２０５５に出力する。その後、ステップＳ２０７に進む。
【００９８】
（ステップＳ２０７）距離判定部２０５５は、評価点選択部２０５４から入力された評価点ξｓ’’と状態推定部１０４から入力された音源状態情報ηｌ｜ｌ−1’が表す音源位置（ｘｌ｜ｌ−１’，ｙｌ｜ｌ−１’）との間の距離が予め定めた閾値、例えば格子点の間隔よりも小さい場合に収束したと判定する。距離判定部２０５５は、収束したと判定した場合、音源の推定位置が収束したことを表す音源収束情報を位置出力部１０６に出力する。また、距離判定部２０５５は、入力された音源状態情報を位置出力部１０６に出力する。その後、処理を終了する。
【００９９】
次に、本実施形態に係る音源位置推定装置２を用いて検証した結果について説明する。
検証において、受聴室として横４ｍ×縦５ｍ×高さ２．４ｍの防音室を用いた。受聴室の内部に、収音部１０１−１〜１０１−Ｎとして８個のマイクロホンをランダムな位置に配置した。受聴室の内部で、実験者は歩行しながら拍手を行う。実験では、この拍手が音源として用いられた。ここで、実験者は５歩進行する毎に１回の拍手を行う。１歩当たりの歩幅は０．３ｍ、時間間隔は０．５秒である。音源の運動モデルとして長方形運動モデル、円運動モデル、各々について想定した。長方形運動モデルを想定した場合、実験者は、横１．２ｍ×縦２．４ｍの長方形の経路上を歩行した。円運動モデルを想定した場合、実験者は、半径１．２ｍの円形の経路上を歩行した。この実験設定のもとで、音源位置推定装置２を、この音源の位置、８個のマイクロホンの位置及び各マイクロホンの観測時刻誤差を推定させた。
【０１００】
音源位置推定装置２の動作条件として、音声信号のサンプリング周波数を１６ｋＨｚとした。処理単位の窓長（ｗｉｎｄｏｗ ｌｅｎｇｔｈ）を５１２サンプル、処理窓のシフト長（ｓｈｉｆｔ ｌｅｎｇｔｈ）を１６０サンプルとした。また、音源から各収音部までの到達時間の観測誤差における標準偏差を０．５×１０−３とし、音源位置の標準偏差を０．１ｍ、音源の観測方向の標準偏差を１度とした。
【０１０１】
図１２は、推定誤差の時間変化の一例を表す図である。
図１２は、運動モデルとして長方形運動モデルを想定した場合における、音源位置の推定誤差、収音部の位置の推定誤差、観測時刻誤差を、各々（ａ）、（ｂ）、（ｃ）に示す。
図１２において、（ａ）の縦軸は音源位置の推定誤差を表し、（ｂ）の縦軸は収音部の位置の推定誤差を表し、（ｃ）の縦軸は観測時刻誤差を表す。但し、（ｂ）に示す推定誤差は、Ｎ個の収音部間における絶対値の平均値である。（ｃ）に示す観測時刻誤差は、Ｎ−１個の収音部間における絶対値の平均値である。（ａ）、（ｂ）、（ｃ）ともに、横軸は時刻を表す。時刻の単位は、拍手の回数である。即ち、横軸の拍手の回数は時刻の目安である。
【０１０２】
図１２によれば、音源位置の推定誤差は、動作開始直後に初期値０．５ｍよりも大きい値２．６ｍとなるが、時間経過に伴いほぼ０に収束する。但し、収束する過程において、時間経過に伴う振動が認められる。この振動は、長方形運動モデルでは音源の移動方向が非線形に変化することが要因であると推定される。音源位置の推定誤差は、拍手回数が１０回以内で、振動による振幅の範囲内に収まる。
収音位置の推定誤差は、初期値０．９ｍから時間経過に伴いほぼ単調に０に収束する。観測時間誤差の推定誤差は、時間経過に伴いほぼ２．４×１０−３ｓと、初期値３．０×１０−３ｓよりも小さい値に収束する。
従って、図１２は、音源位置、収音位置、観測時間誤差ともに、時間経過に伴い高い精度で推定されることを示す。
【０１０３】
図１３は、推定誤差の時間変化のその他の例を表す図である。
図１３は、運動モデルとして円運動モデルを想定した場合における、音源位置の推定誤差、収音部の位置の推定誤差、観測時刻誤差を、各々（ａ）、（ｂ）、（ｃ）に示す。
図１３において、縦軸と横軸の関係は図１２と同様である。
【０１０４】
図１３によれば、音源位置の推定誤差は、初期値３．０ｍから時間経過に伴いほぼ０に収束する。拍手回数が１０回以内で、推定誤差が０に達する。但し、拍手回数が５０回までの間は、長方形運動モデルの場合よりも長い周期で推定誤差が振動する。
収音位置の推定誤差は、時間経過に伴い初期値１．０ｍよりも十分小さい値０．１に収束する。但し、拍手回数１４回付近において音源位置の推定誤差と同時に収音位置の推定誤差も増加する傾向が認められる。
観測時間誤差の推定誤差は、時間経過に伴いほぼ１．１×１０−３ｓと、初期値２．４×１０−３ｓよりも小さい値に収束する。
従って、図１３は、音源位置、収音位置、観測時間誤差ともに、時間経過に伴い高い精度で推定されることを示す。
【０１０５】
図１４は、観測時間誤差の一例を表す表である。
図１４に示す観測時間誤差は、円運動モデルを想定して推定した値であって、時間経過により収束した値である。
図１４は、最左列から右側へ順に、チャネル２〜８の観測時間誤差ｍ２τ〜収音部１０１−８のｍ８τをそれぞれ示す。これらの値の単位は、１０−３秒である。観測時間誤差ｍ２τ〜ｍ８τは、それぞれ、−０．８５、−１．１１、−１．４２、０．８７、−０．９５、−２．８１、−０．１０である。
【０１０６】
図１５は、音源定位状況の一例を表す図である。
図１５において、Ｘ軸は受聴室６０１の横方向の座標軸、Ｙ軸は縦方向の座標軸、Ｚ軸は、帯域出力信号のパワーを表す。原点は、受聴室６０１のＸ−Ｙ平面上の中心を表す。図１５のＸ−Ｙ平面上に、Ｘ＝０又はＹ＝０を表す破線を示す。
図１５に示す、帯域出力信号のパワーは、評価点選択部２０５４が収音部１０１−１〜１０１−Ｎの位置の初期値に基づいて評価点毎に算出した値である。この値は、評価点により大きく異なる。従って、ピーク値をとる評価点が、音源位置として有意でないことを表す。
【０１０７】
図１６は、音源定位状況のその他の例を表す図である。
図１６において、Ｘ軸、Ｙ軸、Ｚ軸の関係は、図１５と同様である。
図１６に示す、帯域出力信号のパワーは、音源が原点に所在している時点であって、収束後の収音部１０１−１〜１０１−Ｎの推定された位置に基づいて評価点毎に算出した値である。この値は、原点においてピーク値をとる。
【０１０８】
図１７は、音源定位状況のその他の例を表す図である。
図１７において、Ｘ軸、Ｙ軸、Ｚ軸の関係は、図１５と同様である。
図１７に示す、帯域出力信号のパワーは、音源が原点に所在しているとき、現実の収音部１０１−１〜１０１−Ｎの位置に基づいて評価点毎に算出した値である。この値は、原点においてピーク値をとる。図１６の結果を考慮すると、収束後の収音部の推定された位置を用いて帯域出力信号のピーク値をとる評価点が、音源位置として正しく推定されることを表す。
【０１０９】
図１８は、収束時間の一例を表す図である。
図１８は、横軸は音源位置が収束するまでの経過時間帯を表し、縦軸は経過時間帯毎の実験回数を示すヒストグラムである。ここで、収束とは、前時刻ｌ−１から現時刻ｌまでの推定された音源位置の変化量が０．０１ｍを下回った時点である。全実験回数は、１００回である。実験毎に、収音部１０１−１〜１０１−８の位置をランダムに変更した。
図１８において、経過時間帯が１０〜１９、２０〜２９、３０〜３９、４０〜４９、５０〜５９，６０〜６９、７０〜７９、８０〜８９、９０〜９９（いずれも拍手回数）の場合、実験回数は、それぞれ２、１６、３１、２４、１２、７、５、２、１である。その他の経過時間帯においては、いずれも実験回数は０回である。
【０１１０】
図１９は、推定された音源位置の誤差の一例を表す図である。
図１９において、横軸は経過時間、縦軸は経過時間毎の音源位置の誤差を表す。図１９は、経過時間毎の平均値同士を結ぶ折れ線グラフと、経過時間毎の最大値及び最小値を結ぶエラーバーを表す。
図１９において、経過時間が０、５０、１００、１５０、２００（いずれも拍手回数）の場合、平均値は、０．９、０．１３、０．１、０．０８、０．０７ｍである。このことも、時間経過とともに誤差が収束することが表される。また、経過時間が０、５０、１００、１５０、２００（いずれも拍手回数）の場合、最大値は、２．２６、０．５、０．４、０．３５、０．３ｍとなり、最小値は、０．４７、０．１０、０．０９、０．０７、０．０６ｍとなる。従って、時間経過とともに最大値と最小値の差が小さくなり、音源位置が安定して推定されることが示される。
【０１１１】
このように、本実施形態によれば、複数のチャネルの入力信号を、予め定めた音源位置の評価点から複数のチャネルの各々に対応するマイクロホンの位置までの位相で補償した信号を加算して得られる評価値を最大にする評価点を定める。また、本実施形態では、定めた評価点と音源状態情報が表す音源位置までの距離に基づいて音源位置の変化が収束したか否かを判断する収束判定部を備える。これにより、音声信号を収録しながら、未知の音源位置を収音部の位置と同時に推定することができる。また、音源位置を安定的に推定でき、推定精度が向上する。
【０１１２】
なお、上述では、音源状態情報が表す音源の位置や収音部１０１−１〜１０１−Ｎの位置が２次元の直交座標系で表される座標値である場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、２次元の直交座標系の代わりに、３次元の直交座標系であってもよいし、極座標系、等、他の変数空間で表される座標系であってもよい。３次元の座標系で表される座標値を扱う場合には、本実施形態ではチャネル数Ｎを少なくとも３よりも大きい整数とする。
【０１１３】
なお、上述では、音源の運動モデルが円運動モデル及び長方形運動モデルである場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、それ以外の運動モデル、例えば、直線運動モデル、正弦波運動モデルであってもよい。
【０１１４】
なお、上述では、位置出力部１０６は、収束判定部１０５から入力された音源状態情報に含まれる音源位置情報を出力する場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、音源状態情報に含まれる音源位置情報、運動方向情報、収音部１０１−１〜１０１−Ｎの位置情報、観測時刻誤差、又はこれらの任意の組み合わせを出力してもよい。
【０１１５】
なお、上述では、収束判定部２０５は、遅延和ビームフォーミング法を用いて推定された評価点と状態推定部１０４から入力された音源状態情報に含まれる音源位置に基づいて音源状態情報が収束したか否かを判断する場合を例にとって説明した。本実施形態では、これには限られない。本実施形態では、遅延和ビームフォーミング法を用いて推定された評価点の代わりに他の方式、例えばＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ Ｓｉｇｎａｌ Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法を用いて推定された音源位置を評価点として用いてもよい。
【０１１６】
なお、上述では、距離判定部２０５５は、入力された音源状態情報を位置出力部１０６に出力する場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、音源状態情報に含まれる音源位置情報の代わりに評価点選択部２０５４から入力された評価点を表す評価点情報を出力してもよい。
【０１１７】
なお、上述した実施形態における音源位置推定装置１、２の一部、例えば、時間差算出部１０３、状態更新部１０４１、状態予測部１０４２、収束判定部１０５、ステアリングベクトル算出部２０５１、周波数領域変換部２０５２、出力算出部２０５３、評価点選択部２０５４、距離判定部２０５５をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、音源位置推定装置１、２に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
また、上述した実施形態における音源位置推定装置１、２の一部、または全部を、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現しても良い音源位置推定装置１、２の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。
【０１１８】
以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
【符号の説明】
【０１１９】
１、２…音源位置推定装置、１０１−１〜１０１−Ｎ…収音部、１０２…信号入力部、
１０３…時間差算出部、１０４…状態推定部、１０４１…状態更新部、
１０４２…状態予測部
１０５、２０５…収束判定部、１０６…位置出力部、
２０５１…ステアリングベクトル算出部、２０５２…周波数領域変換部、
２０５３…出力算出部、２０５４…評価点選択部、２０５５…距離判定部
【技術分野】
【０００１】
本発明は、音源位置推定装置、音源位置推定方法、及び音源位置推定プログラムに関する。
【背景技術】
【０００２】
従来から、音源の方向を推定する音源定位技術が提案されている。音源定位技術は、ロボットが周囲の環境を把握し、又は雑音への耐性を強化するために有用である。音源定位技術では、複数のマイクロホンからなるマイクロホンアレイを用い、チャネル間の音波の到来時刻の差を検出し、マイクロホンの配置に基づいて音源の方向を推定する。そのため、各マイクロホンの位置、もしくは音源と各マイクロホン間の伝達関数が既知であること、チャネル間で音声信号を同期収録すること、いずれも必要である。
【０００３】
そこで、非特許文献１に記載の音源定位技術では、空間的に分散配置した複数のマイクロホンを用い、チャネル間で非同期に音源から音声信号を記録する。当該音源定位技術では、記録を終えた音声信号を用いて音源位置及びマイクロホン位置を推定する。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Ｎ．Ｏｎｏ，Ｈ．Ｋｏｈｎｏ，Ｎ．Ｉｔｏ，ａｎｄ Ｓ．Ｓａｇａｙａｍａ、ＢＬＩＮＤ ＡＬＩＧＮＭＥＮＴ ＯＦ ＡＳＹＮＣＨＲＯＮＯＵＳＬＹ ＲＥＣＯＲＤＥＤ ＳＩＧＮＡＬＳ ＦＯＲ ＤＩＳＴＲＩＢＵＴＥＤ ＭＩＣＲＯＰＨＯＮＥ ＡＲＲＡＹ、「２００９ ＩＥＥＥ Ｗｏｒｋｓｈｏｐ ｏｎ Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ ｔｏ Ａｕｄｉｏ ａｎｄ Ａｃｏｕｓｔｉｃｓ」、ＩＥＥＥ、２００９年１０月１８日、ｐｐ．１６１−１６４
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、非特許文献１に記載の音源定位技術では、音声信号の入力と同時に、音源位置を実時間で推定することができない。
【０００６】
本発明は上記の点に鑑みてなされたものであり、音声信号の入力と同時に音源位置を実時間で推定することができる音源位置推定装置、音源位置推定方法、及び音源位置推定プログラムを提供する。
【課題を解決するための手段】
【０００７】
（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、複数のチャネルの音声信号を入力する信号入力部と、チャネル間の音声信号の時間差を算出する時間差算出部と、音源位置を含む音源状態情報である過去の音源状態情報から現在の音源状態情報を予測する状態予測部と前記時間差算出部が算出した時間差と前記状態予測部が予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する状態更新部とを備えることを特徴とする音源位置推定装置である。
【０００８】
（２）本発明のその他の態様は、上述の音源位置推定装置であって、前記状態更新部は、前記誤差に基づいてカルマンゲインを算出し、算出したカルマンゲインに前記誤差を乗ずることを特徴とする。
【０００９】
（３）本発明のその他の態様は、上述の音源位置推定装置であって、前記音源状態情報は、前記音声信号を前記信号入力部に供給する収音部の位置を含むことを特徴とする。
【００１０】
（４）本発明のその他の態様は、上述の音源位置推定装置であって、前記収音部の位置の変化に基づいて、前記音源位置の変化が収束したか否かを判断する収束判定部を備えることを特徴とする。
【００１１】
（５）本発明のその他の態様は、上述の音源位置推定装置であって、前記複数のチャネルの入力信号を、予め定めた音源位置の評価点から前記複数のチャネルの各々に対応する収音部の位置までの位相で補償した信号を加算して得られる評価値を最大にする評価点を定め、定めた評価点と前記音源状態情報が表す音源位置までの距離に基づいて前記音源位置の変化が収束したか否かを判断する収束判定部を備えることを特徴とする。
【００１２】
（６）本発明のその他の態様は、上述の音源位置推定装置であって、前記収束判定部は、前記評価点を、遅延和ビームフォーミング法を用いて定め、定めた評価点と上述の状態更新部が推定した音源状態情報が表す音源位置までの距離に基づいて、前記音源位置の変化が収束したか否かを判断することを特徴とする
【００１３】
（７）本発明のその他の態様は、音源位置推定装置における方法において、前記音源位置推定装置が、複数のチャネルの音声信号を入力する過程と、前記音源位置推定装置が、チャネル間の音声信号の時間差を算出する過程と、前記音源位置推定装置が、音源位置を含む音源状態情報である過去の音源状態情報から現在の前記音源状態情報を予測する過程と前記音源位置推定装置が、前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する過程を有することを特徴とする音源位置推定方法である。
【００１４】
（８）本発明のその他の態様は、音源位置推定装置のコンピュータに、複数のチャネルの音声信号を入力する手順、チャネル間の音声信号の時間差を算出する手順、音源位置を含む音源状態情報である過去の音源状態情報を予測する手順、前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する手順を実行させるための音源位置推定プログラムである。
【発明の効果】
【００１５】
上述の（１）、（７）、（８）の態様によれば、音声信号の入力と同時に音源位置を実時間で推定することができる。
上述の（２）の態様によれば、音源位置の推定誤差が低減されるように音源位置を安定して推定することができる。
上述の（３）の態様によれば、音源位置とマイクロホンの位置を同時に推定することができる。
上述の（４）、（５）、（６）の態様によれば、誤差が収束した音源位置を取得することができる。
【図面の簡単な説明】
【００１６】
【図１】本発明の第１の実施形態に係る音源位置推定装置の構成を示す概略図である。
【図２】本実施形態に係る収音部の配置例を表す平面図である。
【図３】本実施形態に係る収音部における音源の観測時刻を表す図である。
【図４】音源状態情報の予測及び更新の概要を表す概念図である。
【図５】音源及び本実施形態に係る収音部の位置関係の一例を表す概念図である。
【図６】長方形運動モデルの一例を表す概念図である。
【図７】円運動モデルの一例を表す概念図である。
【図８】本実施形態に係る音源位置推定処理を表すフローチャートである。
【図９】本発明の第２の実施形態に係る音源位置推定装置の構成を示す概略図である。
【図１０】本実施形態に係る収束判定部の構成を表す概略図である。
【図１１】本実施形態に係る収束判定処理を表すフローチャートである。
【図１２】推定誤差の時間変化の一例を表す図である。
【図１３】推定誤差の時間変化のその他の例を表す図である。
【図１４】観測時間誤差の一例を表す表である。
【図１５】音源定位状況の一例を表す図である。
【図１６】音源定位状況のその他の例を表す図である。
【図１７】音源定位状況のその他の例を表す図である。
【図１８】収束時間の一例を表す図である。
【図１９】推定された音源位置の誤差の一例を表す図である。
【発明を実施するための形態】
【００１７】
（第１の実施形態）
以下、図面を参照しながら本発明の実施形態について説明する。
図１は、本実施形態に係る音源位置推定装置１の構成を示す概略図である。
音源位置推定装置１は、Ｎ個の（Ｎは、１よりも大きい整数）収音部１０１−１〜１０１−Ｎと、信号入力部１０２、時間差算出部１０３、状態推定部１０４、収束判定部１０５、及び位置出力部１０６を含んで構成される。
状態推定部１０４は、状態更新部１０４１及び状態予測部１０４２を含んで構成される。
【００１８】
収音部１０１−１〜１０１−Ｎは、空気の振動である音波を電気信号であるアナログ音声信号に変換する電気音響変換器を備える。収音部１０１−１〜１０１−Ｎは、変換したアナログ音声信号を信号入力部１０２に出力する。
収音部１０１−１〜１０１−Ｎは、例えば、音源位置推定装置１の筐体の外部に分散配置されていてもよい。この場合、収音部１０１−１〜１０１−Ｎは、各々、生成した１チャネルの音声信号を無線又は有線で信号入力部１０２に出力する。収音部１０１−１〜１０１−Ｎの各々は、例えば、マイクロホンユニットである。
【００１９】
ここで、収音部１０１−１〜１０１−Ｎの配置例について説明する。
図２は、本実施形態に係る収音部１０１−１〜１０１−８の配置例を表す平面図である。
図２において、横方向がｘ軸方向であり、縦方向がｙ軸方向である。
図２に示された縦長の長方形は、高さ方向（ｚ軸方向）の座標が一定である受聴室６０１の水平面を表す。図２において、黒丸は、各々収音部１０１−１〜１０１−８が配置されている位置を表す。
収音部１０１−１は、受聴室６０１の中央に配置されている。収音部１０１−２は、受聴室６０１の中央からｘ軸の正方向にずれた位置に配置されている。収音部１０１−３は、収音部１０１−２よりもｙ軸の正方向にずれた位置に配置されている。収音部１０１−４は、収音部１０１−３よりもｘ軸の負方向、ｙ軸の正方向にずれた位置に配置されている。収音部１０１−５は、収音部１０１−４よりもｘ軸の負方向、ｙ軸の負方向にずれた位置に配置されている。収音部１０１−６は、収音部１０１−５よりもｙ軸の負方向にずれた位置に配置されている。収音部１０１−７は、収音部１０１−６よりもｘ軸の正方向、ｙ軸の負方向にずれた位置に配置されている。収音部１０１−８は、収音部１０１−７よりもｘ軸の正方向、ｙ軸の正方向にずれた位置であって、収音部１０１−２よりもｙ軸の正方向にずれた位置に配置されている。このように、収音部１０１−２〜１０１−８は、収音部１０１−１を中心にｘ、ｙ平面上を反時計回りに順に配置されている。
【００２０】
図１に戻り、信号入力部１０２には、収音部１０１−１〜１０１−Ｎの各々からのアナログ音声信号が入力される。以下の説明では、収音部１０１−１〜１０１−Ｎに各々対応するチャネルを、チャネル１〜Ｎと呼ぶ。信号入力部１０２は、各チャネルのアナログ音声信号をアナログディジタル（Ａ／Ｄ、Ａｎａｌｏｇ−ｔｏ−Ｄｉｇｉｔａｌ）変換して、ディジタル音声信号を生成する。
信号入力部１０２は、変換した各チャネルのディジタル音声信号を時間差算出部１０３に出力する。
【００２１】
時間差算出部１０３は、信号入力部１０２から入力された音声信号についてチャネル間の時間差を算出する。時間差算出部１０３は、例えば、チャネル１の音声信号と、チャネルｎ（ｎは、１よりも大きく、Ｎと等しい又はＮより小さい整数）の音声信号との時間差ｔｎ，ｋ−ｔ１，ｋ（以下、Δｔｎ，ｋと表す）を算出する。ここで、ｋは、離散時刻を表す整数である。時間差算出部１０３は、時間差Δｔｎ，ｋを算出する際、例えば、チャネル１の音声信号とチャネルｎの音声信号の間で時間差を与えて、両者間の相互相関を算出し、算出した相互相関が最大となる時間差を選択する。
【００２２】
ここで、時間差Δｔｎ，ｋについて図３を用いて説明する。
図３は、収音部１０１−１及び１０１−ｎにおける音源の観測時刻ｔ１，ｋ、ｔｎ，ｋをそれぞれ表す図である。
図３において、横軸が時刻ｔ、縦軸が収音部を表す。図３において、Ｔｋは、音源が音波を発生させた時刻（発音時刻）を表す。ｔ１，ｋは、収音部１０１−１に音源から受信した音波が観測される時刻（観測時刻）を表す。ｔｎ，ｋは、収音部１０１−ｎに音源から受信した音波が観測される観測時刻を表す。観測時刻ｔ１，ｋは、発音時刻Ｔｋにチャネル１における観測時刻誤差ｍ１τと音源から収音部１０１−１までの音波の伝搬時間Ｄ１，ｋ／ｃが加わった時刻である。観測時刻誤差ｍ１τとは、チャネル１の音声信号が観測される時刻の絶対時刻に対する差である。観測時刻誤差が生じる原因は、主に収音部１０１−ｎの位置と音源の位置の計測誤差や収音部１０１−ｎに音波が到達した到達時刻の観測誤差である。観測時刻誤差Ｄ１，ｋは、音源から収音部１０１−ｎまでの距離である。ｃは、音速である。観測時刻ｔｎ，ｋは、発音時刻Ｔｋにチャネルｎにおける観測時刻誤差ｍｎτと音源から収音部１０１−ｎまでの音波の伝搬時間Ｄｎ，ｋ／ｃが加わった時刻である。従って、時間差Δｔｎ，ｋ（＝ｔｎ，ｋ−ｔ１，ｋ）は式（１）で表される。
【００２３】
【数１】
【００２４】
音源から収音部１０１−ｎまでの距離Ｄｎ，ｋは、式（２）で表される。
【００２５】
【数２】
【００２６】
式（２）において、（ｘｋ，ｙｋ）は時刻ｋにおける音源の位置を表す。（ｍｎｘ，ｍｎｙ）は、収音部１０１−ｎの位置を表す。
ここで、各チャネルｎの時間差Δｔｎ，ｋを要素とする（Ｎ−１）列のベクトル［Δｔ２，ｋ，…，Δｔｎ，ｋ，…，ΔｔＮ，ｋ］Ｔを観測値ベクトルζｋと呼ぶ。ここで、Ｔは、行列又はベクトルの転置（ｔｒａｎｓｐｏｓｅ）を表す。時間差算出部１０３は、観測値ベクトルζｋを表す時間差情報を状態推定部１０４に出力する。
【００２７】
図１に戻り、状態推定部１０４は、過去（例えば、時刻ｋ−１）の音源状態情報から現在（時刻ｋ）の音源状態情報を予測し、時間差算出部１０３から入力された時間差情報が表す時間差に基づいて音源状態情報を推定する。音源状態情報は、例えば、音源の位置（ｘｋ，ｙｋ）、各収音部１０１−ｎの位置（ｍｎｘ，ｍｎｙ）及び観測時刻誤差ｍｎτを表す情報を含む。状態推定部１０４は、音源状態情報を推定する際、時間差算出部１０３から入力された時間差情報が表す時間差と予測した音源状態情報に基づく時間差との間の誤差を減少させるように音源状態情報を更新する。状態推定部１０４は、音源状態情報の予測及び更新において、例えば、拡張カルマンフィルタ（Ｅｘｔｅｎｄｅｄ Ｋａｌｍａｎ Ｆｉｌｔｅｒ；ＥＫＦ）法を用いる。ＥＫＦ法を用いた、予測及び更新については後述する。なお、状態推定部１０４は、拡張カルマンフィルタ法の代わりに、最小二乗平均誤差（Ｍｉｎｉｍｕｍ Ｍｅａｎ Ｓｑｕｒｅｄ Ｅｒｒｏｒ；ＭＭＳＥ）法、その他の方式を用いてもよい。
状態推定部１０４は、推定した音源状態情報を収束判定部１０５に出力する。
【００２８】
収束判定部１０５は、状態推定部１０４から入力された音源状態情報ηｋ’が表す音源位置の変化が収束したか否か判断する。収束判定部１０５は、音源の推定位置が収束したことを表す音源収束情報を位置出力部１０６に出力する。ここで、記号’は、推定値であることを表す記号である。
収束判定部１０５は、例えば、過去の収音部１０１−ｎの推定位置（ｍｎｘ，ｋ−１’，ｍｎｙ，ｋ−１’）と現在の収音部１０１−ｎの推定位置（ｍｎｘ，ｋ’，ｍｎｙ，ｋ’）の間の平均距離Δηｍ’を算出する。収束判定部１０５は、平均距離Δηｍ’が予め設定された閾値よりも小さくなったとき収束したと判断する。このように、音源の推定位置を直接収束判断に用いないのは、音源位置は未知であり時間経過によって変化するためである。反面、収音部１０１−ｎの推定位置（ｍｎｘ，ｋ’，ｍｎｙ，ｋ’）を収束判断に用いるのは、収音部１０１−ｎの位置は固定であり、音源状態情報が、音源の推定位置の他、収音部１０１−ｎの推定位置にも依存するためである。
【００２９】
位置出力部１０６は、収束判定部１０５から音源収束情報が入力された場合、収束判定部１０５から入力された音源状態情報に含まれる音源位置情報を外部に出力する。
【００３０】
次に、ＥＫＦ法を用いた音源状態情報の予測及び更新の概要を説明する。
図４は、音源状態情報の予測及び更新の概要を表す概念図である。
図４において、黒塗りの星印は、音源位置の真値を表す。白抜きの星印は、音源位置の推定値を表す。黒丸は、それぞれ収音部１０１−１、１０１−ｎの位置の真値を表す。白丸は、それぞれ収音部１０１−１、１０１−ｎの位置の推定値を示す。収音部１０１−ｎの位置を中心とする実線の円４０１は、収音部１０１−ｎの位置の観測誤差の大きさを表す。収音部１０１−ｎの位置を中心とする一点鎖線の円４０２は、後述する更新ステップを経た後の収音部１０１−ｎの位置の観測誤差の大きさを表す。即ち、円４０１及び４０２は、更新ステップでは、観測誤差が低減されるように収音部１０１−ｎの位置を含む音源状態情報が更新されることを表す。観測誤差は、後述する分散共分散行列Ｐｋ’で定量的に表される。音源の位置を中心とする破線の円４０３は、現実の音源の位置と、音源の移動モデルを用いて推定される音源の位置との間のモデル誤差Ｒを表す円である。モデル誤差は、後述する分散共分散行列Ｒで定量的に表される。
【００３１】
ＥＫＦ法は、Ｉ．観測ステップ、ＩＩ．更新ステップ、ＩＩＩ．予測ステップを含む。状態推定部１０４は、これらのステップを繰り返して実行する。
Ｉ．観測ステップでは、状態推定部１０４は、時間差算出部１０３から時間差情報を入力される。状態推定部１０４は、音源からの音声信号に対する収音部１０１−１、１０１−ｎ間の時間差ΔＴｎ，ｋを表す時間差情報ζｋが観測値として入力される。
ＩＩ．更新ステップでは、状態推定部１０４は、観測値ベクトルζｋと音源状態情報ηｋ’に基づく観測値ベクトルζｋ’との観測誤差が低減されるように、音源状態情報の誤差を表す分散共分散行列Ｐｋ’と音源状態情報ηｋ’を更新する。
ＩＩＩ．予測ステップでは、状態予測部１０４２は、真の音源の位置の時間変化を表す運動モデルに基づき、前時刻ｋ−１の音源状態情報ηｋ−１’から現時刻ｋの音源状態情報ηｋ｜Ｋ−１’を予測する。状態予測部１０４２は、前時刻ｋ−１の分散共分散行列Ｐｋ−１’に音源位置の運動モデルと推定位置とのモデル誤差を表す分散共分散行列Ｒに基づいて分散共分散行列Ｐｋ−１’を更新する。
【００３２】
ここで、音源状態情報ηｋ’は、例えば、音源の推定位置（ｘｋ’，ｙｋ’）、収音部１０１−１〜１０１−Ｎの推定位置（ｍ１ｘ，ｋ’，ｍ１ｙ，ｋ’）〜（ｍＮｘ，ｋ’，ｍＮｙ，ｋ’）及び観測時刻誤差の推定値ｍ１τ’〜ｍＮτ’を要素として含む。つまり、音源状態情報ηｋ’は、例えば、ベクトル［ｘｋ’，ｙｋ’，ｍ１ｘ，ｋ’，ｍ１ｙ，ｋ’，ｍ１τ’，…，ｍＮｘ，ｋ’，ｍＮｙ，ｋ’，ｍＮτ’］Ｔで表わされる情報である。このように、ＥＫＦ法を用いることで、予測誤差が徐々に低減されるように、未知である音源位置、収音部１０１−１〜１０１−Ｎの位置及び観測時刻誤差が予測される。
【００３３】
図１に戻り、状態推定部１０４の構成について説明する。
状態推定部１０４は、状態更新部１０４１と状態予測部１０４２とを含んで構成される。
状態更新部１０４１は、時間差算出部１０３から観測値ベクトルζｋを表す時間差情報が入力される（Ｉ．観測ステップ）。状態更新部１０４１は、状態予測部１０４２から入力された音源状態情報ηｋ｜ｋ−１’と共分散行列Ｐｋ｜ｋ−１が入力される。音源状態情報ηｋ｜ｋ−１’は、前時刻ｋ−１の音源状態情報ηｋ−１’から予測された現時刻ｋの音源状態情報を表す。共分散行列Ｐｋ｜ｋ−１の各要素は、音源状態情報ηｋ｜ｋ−１’が表すベクトルにおける各要素間の共分散である。即ち、この共分散行列Ｐｋ｜ｋ−１は、音源状態情報ηｋ｜ｋ−１’の誤差を表す。その後、状態更新部１０４１は、音源状態情報ηｋ｜ｋ−１’を時刻ｋのηｋ｜ｋ−１’に更新し、共分散行列Ｐｋ｜ｋ−１を共分散行列Ｐｋに更新する（ＩＩ．更新ステップ）。状態更新部１０４１は、更新した現時刻ｋの音源状態情報ηｋ’及び共分散行列Ｐｋを状態予測部１０４２に出力する。
【００３４】
次に、更新ステップにおける更新処理について、より詳細に説明する。
状態更新部１０４１は、観測値ベクトルζｋに観測誤差ベクトルδｋを加算し、加算して得られた和に観測値ベクトルζｋを更新する。観測誤差ベクトルδｋは、平均値が０であり予め定めた共分散で分布しているガウス分布に従う乱数ベクトルである。この共分散を各行各列の要素として含む行列を共分散行列Ｑと表す。
【００３５】
状態更新部１０４１は、音源状態情報ηｋ｜ｋ−１’、共分散行列Ｐｋ｜ｋ−１及び共分散行列Ｑに基づいて、例えば、式（３）を用いてカルマンゲインＫｋを算出する。
【００３６】
【数３】
【００３７】
式（３）において、行列Ｈｋは、式（４）で表されるように観測関数ベクトルｈ（ηｋ｜ｋ−１’）の各要素を、音源状態情報ηｋ｜ｋ−１’の各要素で偏微分して得られるヤコビアンである。
【００３８】
【数４】
【００３９】
観測関数ベクトルｈ（ηｋ’）は、式（５）で表される。
【００４０】
【数５】
【００４１】
観測関数ベクトルｈ（ηｋ’）は、音源状態情報ηｋ’に基づく観測値ベクトルζｋ’である。そこで、状態更新部１０４１は、例えば式（５）を用いて、前時刻ｋ−１の音源状態情報ηｋ−１’から予測された現時刻ｋの音源状態情報ηｋ｜ｋ−１’に対する観測値ベクトルζｋ｜ｋ−１’を算出する。
次に、状態更新部１０４１は、現時刻ｋの観測値ベクトルζｋ、算出した観測値ベクトルζｋ｜ｋ−１’及び算出したカルマンゲインＫｋに基づいて、例えば式（６）を用いて現時刻ｋの音源状態情報ηｋ’を算出する。
【００４２】
【数６】
【００４３】
即ち、式（６）は、前時刻ｋ−１の観測値ベクトルζｋ’から推定された現時刻ｋの観測値ベクトルζｋ｜ｋ−１’に、残差値を加算して現時刻ｋの音源状態情報ηｋ’を算出することを表す。加算される残差値は、観測された現時刻ｋの観測値ベクトルζＫから観測値ベクトルζｋ｜ｋ−１’の差にカルマンゲインＫｋを乗じて得られるベクトル値である。
次に、状態更新部１０４１は、カルマンゲインＫｋ、行列Ｈｋ、及び前時刻ｋ−１の共分散行列Ｐｋ−１から予測された現時刻ｋの共分散行列Ｐｋ｜ｋ−１に基づき、例えば式（７）を用いて現時刻ｋの共分散行列Ｐｋを算出する。
【００４４】
【数７】
【００４５】
式（７）において、Ｉは単位行列を表す。即ち、式（７）は、単位行列ＩからカルマンゲインＫｋと行列Ｈｋとの積を減じて得られた行列を乗じて、音源状態情報ηｋ’の誤差の大きさを低減することを表す。
【００４６】
状態予測部１０４２は、状態更新部１０４１から現時刻ｋの音源状態情報ηｋ’及び共分散行列Ｐｋが入力される。状態予測部１０４２は、前時刻ｋ−１の音源状態情報ηｋ−１’から現時刻ｋの音源状態情報ηｋ｜ｋ−１’を予測し、共分散行列Ｐｋ−１から共分散行列Ｐｋ｜ｋ−1を予測する（ＩＩＩ．予測ステップ）。
【００４７】
次に、予測ステップにおける予測処理について、より詳細に説明する。
本実施形態では、例えば、前時刻ｋ−１における音源位置（ｘｋ−１’，ｙｋ−１’）が、現時刻ｋでの間に、移動量（Δｘ，Δｙ）Ｔだけずれるという運動モデルを仮定する。
状態予測部１０４２は、移動量（Δｘ，Δｙ）Ｔに、その誤差を表す誤差ベクトルεｋを加算して、加算して得られた和に移動量（Δｘ，Δｙ）Ｔを更新する。誤差ベクトルεｋは、平均値が０でありガウス分布に従う乱数ベクトルである。このガウス分布の特性を表す共分散を各行各列の要素として含む行列を共分散行列Ｒと表す。
状態予測部１０４２は、前時刻ｋ−１の音源状態情報ηｋ−１’から現時刻ｋの音源状態情報ηｋ｜ｋ−１’を、例えば式（８）を用いて予測する。
【００４８】
【数８】
【００４９】
式（８）において、行列Ｆηは、式（９）で表される２行２＋３Ｎ列の行列である。
【００５０】
【数９】
【００５１】
次に、状態予測部１０４２は、前時刻ｋ−１の共分散行列Ｐｋ−１から現時刻ｋの共分散行列Ｐｋ｜ｋ−１を、例えば式（１０）を用いて予測する。
【００５２】
【数１０】
【００５３】
即ち、式（１０）は、移動量の誤差を表す共分散行列Ｒに、前時刻ｋ−１の共分散行列Ｐｋ−１で表される音源状態情報ηｋ−１’の誤差を加算して現時刻ｋの共分散行列Ｐｋを算出することを表す。
【００５４】
状態予測部１０４２は、算出した時刻ｋの音源状態情報ηｋ｜ｌｋｋ−１’と共分散行列Ｐｋ｜ｋ−１を状態更新部１０４１に出力する。状態予測部１０４２は、算出した時刻ｋの音源状態情報ηｋ｜ｋ−１’を収束判定部１０５に出力する。
【００５５】
なお、上述では、状態推定部１０４は、Ｉ．観測ステップ、ＩＩ．更新ステップ、ＩＩＩ．予測ステップを時刻ｋ毎に実行する旨、説明したが本実施形態では、これには限られない。本実施形態では、状態推定部１０４は、Ｉ．観測ステップ及びＩＩ．更新ステップを時刻ｋ毎に実行し、ＩＩＩ．予測ステップを、時刻ｌ（エル）毎に実行してもよい。時刻ｌは、時刻ｋとは異なる時間間隔毎に計数される離散時刻である。例えば、前時刻ｌ−１から現時刻ｌまでの時間間隔は、前時刻ｋ−１から現時刻ｋまでの時間間隔よりも広くてもよい。これにより、状態推定部１０４の動作と時間差算出部１０３の動作タイミングが異なっても、相互の処理を同期させることができる。
そこで、状態更新部１０４１は、状態予測部１０４２が出力した時刻ｌの音源状態情報ηｌ｜ｌ−１’を対応する時刻ｋの音源状態情報ηｋ｜ｋ−１’として入力されるようにする。状態予測部１０４２が出力した共分散行列Ｐｌ｜ｌ−１を、状態更新部１０４１は、共分散行列Ｐｋ｜ｋ−１として入力されるようにする。また、状態予測部１０４２は、状態更新部１０４１が出力した音源状態情報ηｋ’を対応する前時刻ｌ−１の音源状態情報ηｌ−１’として入力されるようにする。状態更新部１０４１が出力した共分散行列Ｐｋを、状態予測部１０４２は共分散行列Ｐｌ−１として入力されるようにする。
【００５６】
次に、音源及び収音部１０１−ｎの位置関係の一例について説明する。
図５は、音源及び収音部１０１−ｎの位置関係の一例を表す概念図である。
図５において、黒塗りの★印は、前時刻ｋ−１の音源位置（ｘｋ−１，ｙｋ−１）及び現時刻ｋの音源位置（ｘｋ，ｙｋ）を表す。音源位置（ｘｋ−１，ｙｋ−１）を起点とし、音源位置（ｘｋ，ｙｋ）を終点とする一点破線で表される矢印は、移動量（Δｘ，Δｙ）Ｔを表す。
黒塗りの●印は、収音部１０１−ｎの位置（ｍｎｘ，ｍｎｙ）Ｔを表す。音源位置（ｘｋ，ｙｋ）Ｔを起点とし、収音部１０１−ｎの位置（ｍｎｘ，ｍｎｙ）Ｔを終点とする実線の近傍に表わされているＤｎ，ｋは、これらの間の距離を表す。本実施形態では収音部１０１−ｎの真の位置は定数であると仮定されているが、収音部１０１−ｎの予測値には誤差が含まれている。そのため、収音部１０１−ｎの予測値は変数である。また、距離Ｄｎ，ｋの誤差に対する指標が共分散行列Ｐｋである。
【００５７】
次に、音源の運動モデルの一例として長方形運動モデルについて説明する。
図６は、長方形運動モデルの一例を表す概念図である。
長方形運動モデルは、音源が長方形の軌道上を運動することを仮定する運動モデルである。図６において、横軸がｘ座標、縦軸がｙ座標を表す。図６に表される長方形は、音源が運動する軌道を表す。この長方形のｘ座標の最大値がｘｍａｘ、最小値がｘｍｉｎである。ｙ座標の最大値がｙｍａｘ、最小値がｙｍｉｎである。音源は、長方形の一辺の上を直進し、長方形の一頂点に到達したとき、つまり音源のｘ座標がｘｍａｘもしくはｘｍｉｎ、ｙ座標がｙｍａｘもしくはｙｍｉｎに到達したとき運動方向を９０°回転する。
即ち、長方形運動モデルでは、音源の移動方向θｓ，ｌ−１は、ｘの正方向を基準として０°、９０°、１８０°、−９０°の何れかである。音源が辺上を運動する場合、運動方向の変化量ｄθｓ，ｌ−１Δｔは、０°である。ここで、ｄθｓ，ｌ−１は、音源の角速度を表し、Δｔは、前時刻ｌ−１から現時刻ｌまでの時間間隔を表す。音源が頂点に到達した場合、運動方向の変化量ｄθｓ，ｌ−１Δｔは、反時計回りを正値として９０°又は−９０°である。
【００５８】
長方形運動モデルを用いる場合、本実施形態では、音源位置情報を、２次元の直交座標（ｘｌ，ｘｌ）と運動方向θを要素とする３次元のベクトルηｓ，ｌで表してもよい。音源位置情報ηｓ，ｌは、音源状態情報ηｌに含まれる情報である。この場合、状態予測部１０４２は、式（８）の代わりに式（１１）を用いて音源位置情報の予測を行ってもよい。
【００５９】
【数１１】
【００６０】
式（１１）において、δηは、移動量の誤差ベクトルである。誤差ベクトルδηは、平均値が０であり予め定めた共分散で分布するガウス分布に従う乱数ベクトルである。この共分散を、各行各列の要素として含む行列を共分散行列Ｒと表す。
【００６１】
状態予測部１０４２は、その後、現時刻ｌの共分散行列Ｐｌ｜ｌ−１を、例えば式（１０）の代わりに式（１２）を用いて予測する。
【００６２】
【数１２】
【００６３】
式（１）において、行列Ｇｌは、式（１３）で示される行列である。
【００６４】
【数１３】
【００６５】
式（１３）において、行列Ｆは、式（１４）で示される行列である。
【００６６】
【数１４】
【００６７】
式（１４）において、Ｉ３×３は、３行３列の単位行列であり、Ｏ３×３は、３行３Ｎ列の零行列である。
【００６８】
次に、音源の運動モデルの一例として円運動モデルについて説明する。
図７は、円運動モデルの一例を表す概念図である。
円運動モデルは、音源が円軌道上を運動することを仮定する運動モデルである。図７において、横軸がｘ座標、縦軸がｙ座標を表す。図７に表される円は、音源が運動する軌道を表す。円運動モデルでは、運動方向の変化量ｄθｓ，ｌ−１Δｔが、一定値Δθであり、音源方向もこれに応じて変化する。
【００６９】
円運動モデルを用いる場合も、音源位置情報を、２次元の直交座標（ｘｌ，ｘｌ）と運動方向θを要素とする３次元のベクトルηｓ，ｌで表してもよい。この場合、状態予測部１０４２は、式（８）の代わりに式（１５）を用いて音源位置情報の予測を行う。
【００７０】
【数１５】
【００７１】
状態予測部１０４２は、現時刻ｌの共分散行列Ｐｌ｜ｌ−１を、式（１２）を用いて予測する。但し、行列Ｇｌとして、式（１３）に表される行列Ｇｌの代わりに、式（１６）に表される行列Ｇｌを用いる。
【００７２】
【数１６】
【００７３】
次に、本実施形態に係る音源位置推定処理について説明する。
図８は、本実施形態に係る音源位置推定処理を表すフローチャートである。
（ステップＳ１０１）音源位置推定装置１は、取り扱う変数の初期値を設定する。例えば、状態推定部１０４は、観測時刻ｋ、予測時刻ｌを、それぞれ０と設定し、音源状態情報ηｋ｜ｋ−１と共分散行列Ｐｋ｜ｋ−１をそれぞれ予め定めた値に設定する。その後、ステップＳ１０２に進む。
（ステップＳ１０２）信号入力部１０２は、収音部１０１−１〜１０１−Ｎからチャネル毎の音声信号が各々入力される。信号入力部１０２は、音声信号の入力を継続するか否か判断する。入力を継続する場合（ステップＳ１０２ Ｙ）、信号入力部１０２は、入力された音声信号をＡ／Ｄ変換して時間差算出部１０３に出力し、その後、ステップＳ１０３に進む。入力を継続しない場合（ステップＳ１０２ Ｎ）、処理を終了する。
【００７４】
（ステップＳ１０３）時間差算出部１０３は、信号入力部１０２から入力された音声信号についてチャネル間の時間差を算出する。時間差算出部１０３は、算出されたチャネル間の時間差を要素とする観測値ベクトルζｋを表す時間差情報を状態更新部１０４１に出力する。その後、ステップＳ１０４に進む。
（ステップＳ１０４）状態更新部１０４１は、予め定めた時間毎に観測時刻ｋを１増加させて観測時刻ｋを更新する。その後、ステップＳ１０５に進む。
【００７５】
（ステップＳ１０５）状態更新部１０４１は、時間差算出部１０３から入力された時間差情報が表す観測値ベクトルζｋに観測誤差ベクトルδｋを加算して観測値ベクトルζｋを更新する。
状態更新部１０４１は、音源状態情報ηｋ｜ｋ−１’、共分散行列Ｐｋ｜ｋ−１及び共分散行列Ｑに基づいて、例えば、式（３）を用いてカルマンゲインＫｋを算出する。
状態更新部１０４１は、例えば式（５）を用いて、現観測時刻ｋの音源状態情報ηｋ｜ｋ−１’に対する観測値ベクトルζｋ｜ｋ−１’を算出する。
状態更新部１０４１は、現観測時刻ｋの観測値ベクトルζｋ、算出した観測値ベクトルζｋ｜ｋ−１’及び算出したカルマンゲインＫｋに基づいて、例えば式（６）を用いて現観測時刻ｋの音源状態情報ηｋ’を算出する。
状態更新部１０４１は、カルマンゲインＫｋ、行列Ｈｋ、及び共分散行列Ｐｋ｜ｋ−１に基づき、例えば式（７）を用いて現観測時刻ｋの共分散行列Ｐｋを算出する。その後、ステップＳ１０６に進む。
【００７６】
（ステップＳ１０６）状態更新部１０４１は、現観測時刻ｋが、予測処理を行う予測時刻ｌに相当するか否か判断する。例えば、観測及び更新ステップＮ回（Ｎは、１又は１よりも大きい整数、例えば、５）毎に予測ステップを１回行う場合、観測時刻ｋのＮに対する剰余が０であるか判断する。現観測時刻ｋが予測時刻ｌと判断された場合（ステップＳ１０７ Ｙ）、ステップＳ１０７に進む。現観測時刻ｋが予測時刻ｌと判断されない場合（ステップＳ１０７ Ｎ）、ステップＳ１０２に進む。
【００７７】
（ステップＳ１０７）状態予測部１０４２は、状態更新部１０４１が出力した算出した現観測時刻ｋの音源状態情報ηｋ’及び共分散行列Ｐｋを、前予測時刻ｌ−１の音源状態情報ηｌ−１’及び共分散行列Ｐｌ−１として入力される。
状態予測部１０４２は、前予測時刻ｌ−１の音源状態情報ηｌ−１’から現予測時刻ｌの音源状態情報ηｌ｜ｌ−１’を、例えば式（８）、（１１）又は（１５）を用いて算出する。
状態予測部１０４２は、前予測時刻ｌ−１の共分散行列Ｐｌ−１から現予測時刻ｌの共分散行列Ｐｌ｜ｌ−１を、例えば式（１０）又は（１２）を用いて算出する。
状態予測部１０４２は、現予測時刻ｌの音源状態情報ηｌ｜ｌ−１’と共分散行列Ｐｌ｜ｌ−１を状態更新部１０４１に出力する。状態予測部１０４２は、算出した現予測時刻ｌの音源状態情報ηｌ｜ｌ−１’を、収束判定部１０５に出力する。その後、ステップＳ１０８に進む。
【００７８】
（ステップＳ１０８）状態更新部１０４１は、現予測時刻ｌに１を加えて予測時刻を更新する。状態更新部１０４１は、状態予測部１０４２が出力した予測時刻ｌの音源状態情報ηｌ｜ｌ−１’、共分散行列Ｐｌ｜ｌ−１を、観測時刻ｋのηｋ｜ｋ−１’、共分散行列Ｐｋ｜ｋ−１として入力される。その後、ステップＳ１０９に進む。
【００７９】
（ステップＳ１０９）収束判定部１０５は、状態推定部１０４から入力された音源状態情報ηｌ’が表す音源位置の変化が収束したか否か判断する。収束判定部１０５は、例えば、過去の収音部１０１−ｎの推定位置と現在の収音部１０１−ｎの推定位置の間の平均距離Δηｍ’が予め設定された閾値よりも小さくなったとき収束したと判断する。音源位置の変化が収束したと判断された場合（ステップＳ１０９ Ｙ）、収束判定部１０５は、入力された音源状態情報ηｌ’を位置出力部１０６に出力する。その後、ステップＳ１１０に進む。音源位置の変化が収束したと判断されなかった場合（ステップＳ１０９ Ｎ）、ステップＳ１０２に進む。
（ステップＳ１１０）位置出力部１０６は、収束判定部１０５から入力された音源状態情報に含まれる音源位置情報を外部に出力する。その後、ステップＳ１０２に進む。
【００８０】
このように、本実施形態は、複数のチャネルの音声信号を入力し、チャネル間の音声信号の時間差を算出し、過去の音源位置を含む音源状態情報から現在の前記音源状態情報を予測する。また、本実施形態は、算出した時間差と予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を更新する。これにより、音声信号の入力と同時に音源位置を推定することができる。
【００８１】
（第２の実施形態）
以下、図面を参照しながら本発明の実施形態について説明する。第１の実施形態と同一の構成又は同一の処理については、同一の番号を付す。
図９は、本実施形態に係る音源位置推定装置２の構成を示す概略図である。
音源位置推定装置２は、Ｎ個の収音部１０１−１〜１０１−Ｎと、信号入力部１０２、時間差算出部１０３、状態推定部１０４、収束判定部２０５、及び位置出力部１０６を含んで構成される。即ち、音源位置推定装置２は、音源位置推定装置１（図１参照）の収束判定部１０５の代わりに収束判定部２０５を備え、信号入力部１０２が入力された音声信号を収束判定部２０５にも出力する点が異なる。その他の構成については、音源位置推定装置１と同様である。
【００８２】
次に、収束判定部２０５の構成について説明する。
図１０は、本実施形態に係る収束判定部２０５の構成を表す概略図である。
収束判定部２０５は、ステアリングベクトル（ｓｔｅｅｒｉｎｇ ｖｅｃｔｏｒ）算出部２０５１、周波数領域変換部２０５２、出力算出部２０５３、評価点選択部２０５４、及び距離判定部２０５５を含んで構成される。この構成により、収束判定部２０５は、遅延和ビームフォーミング（Ｄｅｌａｙ−ａｎｄ−Ｓｕｍ Ｂｅａｍｆｏｒｍｉｎｇ，ＤＳ−ＢＦ）法によって推定された評価点と状態推定部１０４から入力された音源状態情報に含まれる音源位置を比較する。ここで、収束判定部２０５は、評価点と音源位置に基づいて音源状態情報が収束したか否かを判断する。
【００８３】
ステアリングベクトル算出部２０５１は、状態予測部１０４２から入力された音源状態情報ηｌ｜ｌ−１’が表す収音部１０１−ｎの位置（ｍｎｘ’，ｍｎｙ’）から音源位置の候補（以下、評価点と呼ぶ）ξｓ’’までの距離Ｄｎ，ｌを算出する。ステアリングベクトル算出部２０５１は、距離Ｄｎ，ｌを算出する際、例えば式（２）を用いる。但し、ステアリングベクトル算出部２０５１は、式（２）の（ｘｋ，ｙｋ）の代わりに評価点ξｓ’’の座標（ｘ’’，ｙ’’）を代入する。この評価点ξｓ’’は、例えば、予め定められた格子点であって、音源が配置されうる空間（例えば、図２に示す受聴室６０１）に配置された複数の格子点の１つである。
ステアリングベクトル算出部２０５１は、算出した距離Ｄｎ，ｌに基づく伝搬遅延Ｄｎ，ｌ／ｃと推定された観測時刻誤差ｍｎτ’を加算してチャネル毎の推定観測時刻ｔｎ，ｌ’’を算出する。ステアリングベクトル算出部２０５１は、算出した推定時間差ｔｎ，ｌ’’に基づいて、ステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を、例えば式（１７）を用いて周波数ω毎に算出する。
【００８４】
【数１７】
【００８５】
式（１７）において、ξｍ’は、収音部１０１−１〜１０１−Ｎの位置の集合を表す。従って、ステアリングベクトルＷ（η’，ω）の各要素は、対応するチャネルｎ（ｎは、１と等しい又は１よりも大きく、Ｎと等しい又はＮよりも小さい）における音源から各収音部１０１−ｎまでの伝搬によって生じた位相の遅延を与える伝達関数である。ステアリングベクトル算出部２０５１は、算出したステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を出力算出部２０５３に出力する。
【００８６】
周波数領域変換部２０５２は、信号入力部１０２から入力された各チャネルの音声信号ｓｎに対して時間領域から周波数領域に変換して、チャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を生成する。周波数領域変換部２０５２は、周波数領域に変換する方式として、例えば、離散フーリエ変換（Ｄｉｓｃｒｅｔｅ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ、ＤＦＴ）を用いる。周波数領域変換部２０５２は、生成したチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を出力算出部２０５３に出力する。
【００８７】
出力算出部２０５３は、周波数領域変換部２０５２からチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）が入力され、ステアリングベクトル算出部２０５１からステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）が入力される。出力算出部２０５３は、周波数領域信号Ｓｎ，ｌ（ω）を要素とする入力信号ベクトルＳｌ（ω）とステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）の内積Ｐ（ξｓ’’，ξｍ’，ω）を算出する。入力信号ベクトルＳｌ（ω）は［Ｓ１，ｌ（ω），…，Ｓｎ，ｌ（ω），…，ＳＮ，ｌ（ω）］Ｔと表される。出力算出部２０５３は、例えば式（１８）を用いて内積Ｐ（ξｓ’’，ξｍ’，ω）を算出する。
【００８８】
【数１８】
【００８９】
式（１８）において、＊は、ベクトル又は行列の複素共役転置（ｃｏｍｐｌｅｘ ｃｏｎｊｕｇａｔｅ ｔｒａｎｓｐｏｓｅ）を表す。式（１８）によれば、入力信号ベクトルＳｋ（ω）の各チャネル成分の伝搬遅延による位相が補償され、各チャネル成分がチャネル間で同期する。そして、位相が補償された各チャネル成分がチャネル間で加算される。
出力算出部２０５３は、算出した内積Ｐ（ξｓ’’，ξｍ’，ω）を、例えば式（１９）を用いて予め定めた周波数帯域にわたって累積して、帯域出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出する。
【００９０】
【数１９】
【００９１】
式（５）において、最低周波数ωｌ（例えば、２００Ｈｚ）を表し、最高周波数ωｈ（例えば、７ｋＨｚ）を表す。
出力算出部２０５３は、算出した帯域出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を評価点選択部２０５４に出力する。
評価点選択部２０５４は、出力算出部２０５３から入力された帯域出力信号＜Ｐ（ξｓ’’，ξｍ’）＞の絶対値を評価値として最大となる評価点ξｓ’’を選択する。評価点選択部２０５４は、選択した評価点ξｓ’’を距離判定部２０５５に出力する。
距離判定部２０５５は、評価点選択部２０５４から入力された評価点ξｓ’’と状態予測部１０４２から入力された音源状態情報ηｌ｜ｌ−１’が表す音源位置（ｘｌ｜ｌ−１’，ｙｌ｜ｌ−１’）との間の距離が予め定めた閾値、例えば、上述の格子点の間隔よりも小さい場合に収束したと判定する。距離判定部２０５５は、収束したと判定した場合、音源の推定位置が収束したことを表す音源収束情報を位置出力部１０６に出力する。また、距離判定部２０５５は、入力された音源状態情報を位置出力部１０６に出力する。
【００９２】
次に、収束判定部２０５における収束判定処理について説明する。
図１１は、本実施形態に係る収束判定処理を表すフローチャートである。
（ステップＳ２０１）周波数領域変換部２０５２は、信号入力部１０２から入力された各チャネルの音声信号ｓｎに対して時間領域から周波数領域に変換し、チャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を生成する。周波数領域変換部２０５２は、生成したチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）を出力算出部２０５３に出力する。その後、ステップＳ２０２に出力する。
【００９３】
（ステップＳ２０２）ステアリングベクトル算出部２０５１は、状態推定部１０４から入力された音源状態情報が表す収音部１０１−ｎの位置（ｍｎｘ’，ｍｎｙ’）から評価点ξｓ’’までの距離Ｄｎ，ｌを算出する。ステアリングベクトル算出部２０５１は、算出した距離Ｄｎ，ｌに基づく伝搬遅延Ｄｎ，ｌ／ｃに推定された観測時刻誤差ｍｎτ’を加算してチャネル毎の推定観測時刻ｔｎ，ｌ’’を算出する。ステアリングベクトル算出部２０５１は、算出した推定時間差ｔｎ，ｌ’’に基づいて、ステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を周波数ω毎に算出する。ステアリングベクトル算出部２０５１は、算出したステアリングベクトルＷ（ξｓ’’，ξｍ’，ω）を出力算出部２０５３に出力する。その後、ステップＳ２０３に出力する。
【００９４】
（ステップＳ２０３）出力算出部２０５３は、周波数領域変換部２０５２からチャネル毎の周波数領域信号Ｓｎ，ｌ（ω）が入力され、ステアリングベクトル算出部２０５１からステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）が入力される。出力算出部２０５３は、周波数領域信号Ｓｎ，ｌ（ω）を要素とする入力信号ベクトルＳｌ（ω）とステアリングベクトルのＷ（ξｓ’’，ξｍ’，ω）の内積Ｐ（ξｓ’’，ξｍ’，ω）を、例えば式（１８）を用いて算出する。
出力算出部２０５３は、算出した内積Ｐ（ξｓ’’，ξｍ’，ω）を、例えば式（１９）を用いて予め定めた周波数帯域にわたって累積し、出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出する。出力算出部２０５３は、算出した出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を評価点選択部２０５４に出力する。その後、ステップＳ２０４に進む。
【００９５】
（ステップＳ２０４）出力算出部２０５３は、全ての評価点について出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出したか否かを判断する。全ての評価点について算出したと判断された場合（ステップＳ２０４ Ｙ）、ステップＳ２０６に進む。全ての評価点について算出していないと判断された場合（ステップＳ２０４ Ｎ）、ステップＳ２０５に進む。
【００９６】
（ステップＳ２０５）出力算出部２０５３は、出力信号＜Ｐ（ξｓ’’，ξｍ’）＞を算出する評価点を、出力信号を算出していない他の評価点に変更する。その後、ステップＳ２０２に進む。
【００９７】
（ステップＳ２０６）評価点選択部２０５４は、出力算出部２０５３から入力された出力信号＜Ｐ（ξｓ’’，ξｍ’）＞の絶対値を評価値として最大となる評価点ξｓ’’を選択する。評価点選択部２０５４は、選択した評価点ξｓ’’を距離判定部２０５５に出力する。その後、ステップＳ２０７に進む。
【００９８】
（ステップＳ２０７）距離判定部２０５５は、評価点選択部２０５４から入力された評価点ξｓ’’と状態推定部１０４から入力された音源状態情報ηｌ｜ｌ−1’が表す音源位置（ｘｌ｜ｌ−１’，ｙｌ｜ｌ−１’）との間の距離が予め定めた閾値、例えば格子点の間隔よりも小さい場合に収束したと判定する。距離判定部２０５５は、収束したと判定した場合、音源の推定位置が収束したことを表す音源収束情報を位置出力部１０６に出力する。また、距離判定部２０５５は、入力された音源状態情報を位置出力部１０６に出力する。その後、処理を終了する。
【００９９】
次に、本実施形態に係る音源位置推定装置２を用いて検証した結果について説明する。
検証において、受聴室として横４ｍ×縦５ｍ×高さ２．４ｍの防音室を用いた。受聴室の内部に、収音部１０１−１〜１０１−Ｎとして８個のマイクロホンをランダムな位置に配置した。受聴室の内部で、実験者は歩行しながら拍手を行う。実験では、この拍手が音源として用いられた。ここで、実験者は５歩進行する毎に１回の拍手を行う。１歩当たりの歩幅は０．３ｍ、時間間隔は０．５秒である。音源の運動モデルとして長方形運動モデル、円運動モデル、各々について想定した。長方形運動モデルを想定した場合、実験者は、横１．２ｍ×縦２．４ｍの長方形の経路上を歩行した。円運動モデルを想定した場合、実験者は、半径１．２ｍの円形の経路上を歩行した。この実験設定のもとで、音源位置推定装置２を、この音源の位置、８個のマイクロホンの位置及び各マイクロホンの観測時刻誤差を推定させた。
【０１００】
音源位置推定装置２の動作条件として、音声信号のサンプリング周波数を１６ｋＨｚとした。処理単位の窓長（ｗｉｎｄｏｗ ｌｅｎｇｔｈ）を５１２サンプル、処理窓のシフト長（ｓｈｉｆｔ ｌｅｎｇｔｈ）を１６０サンプルとした。また、音源から各収音部までの到達時間の観測誤差における標準偏差を０．５×１０−３とし、音源位置の標準偏差を０．１ｍ、音源の観測方向の標準偏差を１度とした。
【０１０１】
図１２は、推定誤差の時間変化の一例を表す図である。
図１２は、運動モデルとして長方形運動モデルを想定した場合における、音源位置の推定誤差、収音部の位置の推定誤差、観測時刻誤差を、各々（ａ）、（ｂ）、（ｃ）に示す。
図１２において、（ａ）の縦軸は音源位置の推定誤差を表し、（ｂ）の縦軸は収音部の位置の推定誤差を表し、（ｃ）の縦軸は観測時刻誤差を表す。但し、（ｂ）に示す推定誤差は、Ｎ個の収音部間における絶対値の平均値である。（ｃ）に示す観測時刻誤差は、Ｎ−１個の収音部間における絶対値の平均値である。（ａ）、（ｂ）、（ｃ）ともに、横軸は時刻を表す。時刻の単位は、拍手の回数である。即ち、横軸の拍手の回数は時刻の目安である。
【０１０２】
図１２によれば、音源位置の推定誤差は、動作開始直後に初期値０．５ｍよりも大きい値２．６ｍとなるが、時間経過に伴いほぼ０に収束する。但し、収束する過程において、時間経過に伴う振動が認められる。この振動は、長方形運動モデルでは音源の移動方向が非線形に変化することが要因であると推定される。音源位置の推定誤差は、拍手回数が１０回以内で、振動による振幅の範囲内に収まる。
収音位置の推定誤差は、初期値０．９ｍから時間経過に伴いほぼ単調に０に収束する。観測時間誤差の推定誤差は、時間経過に伴いほぼ２．４×１０−３ｓと、初期値３．０×１０−３ｓよりも小さい値に収束する。
従って、図１２は、音源位置、収音位置、観測時間誤差ともに、時間経過に伴い高い精度で推定されることを示す。
【０１０３】
図１３は、推定誤差の時間変化のその他の例を表す図である。
図１３は、運動モデルとして円運動モデルを想定した場合における、音源位置の推定誤差、収音部の位置の推定誤差、観測時刻誤差を、各々（ａ）、（ｂ）、（ｃ）に示す。
図１３において、縦軸と横軸の関係は図１２と同様である。
【０１０４】
図１３によれば、音源位置の推定誤差は、初期値３．０ｍから時間経過に伴いほぼ０に収束する。拍手回数が１０回以内で、推定誤差が０に達する。但し、拍手回数が５０回までの間は、長方形運動モデルの場合よりも長い周期で推定誤差が振動する。
収音位置の推定誤差は、時間経過に伴い初期値１．０ｍよりも十分小さい値０．１に収束する。但し、拍手回数１４回付近において音源位置の推定誤差と同時に収音位置の推定誤差も増加する傾向が認められる。
観測時間誤差の推定誤差は、時間経過に伴いほぼ１．１×１０−３ｓと、初期値２．４×１０−３ｓよりも小さい値に収束する。
従って、図１３は、音源位置、収音位置、観測時間誤差ともに、時間経過に伴い高い精度で推定されることを示す。
【０１０５】
図１４は、観測時間誤差の一例を表す表である。
図１４に示す観測時間誤差は、円運動モデルを想定して推定した値であって、時間経過により収束した値である。
図１４は、最左列から右側へ順に、チャネル２〜８の観測時間誤差ｍ２τ〜収音部１０１−８のｍ８τをそれぞれ示す。これらの値の単位は、１０−３秒である。観測時間誤差ｍ２τ〜ｍ８τは、それぞれ、−０．８５、−１．１１、−１．４２、０．８７、−０．９５、−２．８１、−０．１０である。
【０１０６】
図１５は、音源定位状況の一例を表す図である。
図１５において、Ｘ軸は受聴室６０１の横方向の座標軸、Ｙ軸は縦方向の座標軸、Ｚ軸は、帯域出力信号のパワーを表す。原点は、受聴室６０１のＸ−Ｙ平面上の中心を表す。図１５のＸ−Ｙ平面上に、Ｘ＝０又はＹ＝０を表す破線を示す。
図１５に示す、帯域出力信号のパワーは、評価点選択部２０５４が収音部１０１−１〜１０１−Ｎの位置の初期値に基づいて評価点毎に算出した値である。この値は、評価点により大きく異なる。従って、ピーク値をとる評価点が、音源位置として有意でないことを表す。
【０１０７】
図１６は、音源定位状況のその他の例を表す図である。
図１６において、Ｘ軸、Ｙ軸、Ｚ軸の関係は、図１５と同様である。
図１６に示す、帯域出力信号のパワーは、音源が原点に所在している時点であって、収束後の収音部１０１−１〜１０１−Ｎの推定された位置に基づいて評価点毎に算出した値である。この値は、原点においてピーク値をとる。
【０１０８】
図１７は、音源定位状況のその他の例を表す図である。
図１７において、Ｘ軸、Ｙ軸、Ｚ軸の関係は、図１５と同様である。
図１７に示す、帯域出力信号のパワーは、音源が原点に所在しているとき、現実の収音部１０１−１〜１０１−Ｎの位置に基づいて評価点毎に算出した値である。この値は、原点においてピーク値をとる。図１６の結果を考慮すると、収束後の収音部の推定された位置を用いて帯域出力信号のピーク値をとる評価点が、音源位置として正しく推定されることを表す。
【０１０９】
図１８は、収束時間の一例を表す図である。
図１８は、横軸は音源位置が収束するまでの経過時間帯を表し、縦軸は経過時間帯毎の実験回数を示すヒストグラムである。ここで、収束とは、前時刻ｌ−１から現時刻ｌまでの推定された音源位置の変化量が０．０１ｍを下回った時点である。全実験回数は、１００回である。実験毎に、収音部１０１−１〜１０１−８の位置をランダムに変更した。
図１８において、経過時間帯が１０〜１９、２０〜２９、３０〜３９、４０〜４９、５０〜５９，６０〜６９、７０〜７９、８０〜８９、９０〜９９（いずれも拍手回数）の場合、実験回数は、それぞれ２、１６、３１、２４、１２、７、５、２、１である。その他の経過時間帯においては、いずれも実験回数は０回である。
【０１１０】
図１９は、推定された音源位置の誤差の一例を表す図である。
図１９において、横軸は経過時間、縦軸は経過時間毎の音源位置の誤差を表す。図１９は、経過時間毎の平均値同士を結ぶ折れ線グラフと、経過時間毎の最大値及び最小値を結ぶエラーバーを表す。
図１９において、経過時間が０、５０、１００、１５０、２００（いずれも拍手回数）の場合、平均値は、０．９、０．１３、０．１、０．０８、０．０７ｍである。このことも、時間経過とともに誤差が収束することが表される。また、経過時間が０、５０、１００、１５０、２００（いずれも拍手回数）の場合、最大値は、２．２６、０．５、０．４、０．３５、０．３ｍとなり、最小値は、０．４７、０．１０、０．０９、０．０７、０．０６ｍとなる。従って、時間経過とともに最大値と最小値の差が小さくなり、音源位置が安定して推定されることが示される。
【０１１１】
このように、本実施形態によれば、複数のチャネルの入力信号を、予め定めた音源位置の評価点から複数のチャネルの各々に対応するマイクロホンの位置までの位相で補償した信号を加算して得られる評価値を最大にする評価点を定める。また、本実施形態では、定めた評価点と音源状態情報が表す音源位置までの距離に基づいて音源位置の変化が収束したか否かを判断する収束判定部を備える。これにより、音声信号を収録しながら、未知の音源位置を収音部の位置と同時に推定することができる。また、音源位置を安定的に推定でき、推定精度が向上する。
【０１１２】
なお、上述では、音源状態情報が表す音源の位置や収音部１０１−１〜１０１−Ｎの位置が２次元の直交座標系で表される座標値である場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、２次元の直交座標系の代わりに、３次元の直交座標系であってもよいし、極座標系、等、他の変数空間で表される座標系であってもよい。３次元の座標系で表される座標値を扱う場合には、本実施形態ではチャネル数Ｎを少なくとも３よりも大きい整数とする。
【０１１３】
なお、上述では、音源の運動モデルが円運動モデル及び長方形運動モデルである場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、それ以外の運動モデル、例えば、直線運動モデル、正弦波運動モデルであってもよい。
【０１１４】
なお、上述では、位置出力部１０６は、収束判定部１０５から入力された音源状態情報に含まれる音源位置情報を出力する場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、音源状態情報に含まれる音源位置情報、運動方向情報、収音部１０１−１〜１０１−Ｎの位置情報、観測時刻誤差、又はこれらの任意の組み合わせを出力してもよい。
【０１１５】
なお、上述では、収束判定部２０５は、遅延和ビームフォーミング法を用いて推定された評価点と状態推定部１０４から入力された音源状態情報に含まれる音源位置に基づいて音源状態情報が収束したか否かを判断する場合を例にとって説明した。本実施形態では、これには限られない。本実施形態では、遅延和ビームフォーミング法を用いて推定された評価点の代わりに他の方式、例えばＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ Ｓｉｇｎａｌ Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法を用いて推定された音源位置を評価点として用いてもよい。
【０１１６】
なお、上述では、距離判定部２０５５は、入力された音源状態情報を位置出力部１０６に出力する場合を例にとって説明したが、本実施形態ではこれには限られない。本実施形態では、音源状態情報に含まれる音源位置情報の代わりに評価点選択部２０５４から入力された評価点を表す評価点情報を出力してもよい。
【０１１７】
なお、上述した実施形態における音源位置推定装置１、２の一部、例えば、時間差算出部１０３、状態更新部１０４１、状態予測部１０４２、収束判定部１０５、ステアリングベクトル算出部２０５１、周波数領域変換部２０５２、出力算出部２０５３、評価点選択部２０５４、距離判定部２０５５をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、音源位置推定装置１、２に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
また、上述した実施形態における音源位置推定装置１、２の一部、または全部を、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現しても良い音源位置推定装置１、２の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。
【０１１８】
以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
【符号の説明】
【０１１９】
１、２…音源位置推定装置、１０１−１〜１０１−Ｎ…収音部、１０２…信号入力部、
１０３…時間差算出部、１０４…状態推定部、１０４１…状態更新部、
１０４２…状態予測部
１０５、２０５…収束判定部、１０６…位置出力部、
２０５１…ステアリングベクトル算出部、２０５２…周波数領域変換部、
２０５３…出力算出部、２０５４…評価点選択部、２０５５…距離判定部
【特許請求の範囲】
【請求項１】
複数のチャネルの音声信号を入力する信号入力部と、
チャネル間の音声信号の時間差を算出する時間差算出部と、
音源位置を含む音源状態情報である過去の音源状態情報から現在の音源状態情報を予測する状態予測部と
前記時間差算出部が算出した時間差と前記状態予測部が予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する状態更新部とを備える
ことを特徴とする音源位置推定装置。
【請求項２】
前記状態更新部は、
前記誤差に基づいてカルマンゲインを算出し、算出したカルマンゲインに前記誤差を乗ずることを特徴とする請求項１に記載の音源位置推定装置。
【請求項３】
前記音源状態情報は、前記音声信号を前記信号入力部に供給する収音部の位置を含むことを特徴とする請求項１又は２に記載の音源位置推定装置。
【請求項４】
前記収音部の位置の変化に基づいて、前記音源位置の変化が収束したか否かを判断する収束判定部を備えること
を特徴とする請求項３に記載の音源位置推定装置。
【請求項５】
前記複数のチャネルの入力信号を、予め定めた音源位置の評価点から前記複数のチャネルの各々に対応する収音部の位置までの位相で補償した信号を加算して得られる評価値を最大にする評価点を定め、定めた評価点と前記状態更新部が推定した音源状態情報が表す音源位置までの距離に基づいて、前記音源位置の変化が収束したか否かを判断する収束判定部を備えること
を特徴とする請求項３に記載の音源位置推定装置。
【請求項６】
前記収束判定部は、
前記評価点を、遅延和ビームフォーミング法を用いて定め、定めた評価点と請求項２に記載の状態更新部が推定した音源状態情報が表す音源位置までの距離に基づいて、前記音源位置の変化が収束したか否かを判断することを特徴とする請求項５に記載の音源位置推定装置。
【請求項７】
音源位置推定装置における方法において、
前記音源位置推定装置が、複数のチャネルの音声信号を入力する過程と、
前記音源位置推定装置が、チャネル間の音声信号の時間差を算出する過程と、
前記音源位置推定装置が、音源位置を含む音源状態情報である過去の音源状態情報から現在の前記音源状態情報を予測する過程と
前記音源位置推定装置が、前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する過程を有すること
を特徴とする音源位置推定方法。
【請求項８】
音源位置推定装置のコンピュータに、
複数のチャネルの音声信号を入力する手順、
チャネル間の音声信号の時間差を算出する手順、
音源位置を含む音源状態情報である過去の音源状態情報を予測する手順、
前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する手順
を実行させるための音源位置推定プログラム。
【請求項１】
複数のチャネルの音声信号を入力する信号入力部と、
チャネル間の音声信号の時間差を算出する時間差算出部と、
音源位置を含む音源状態情報である過去の音源状態情報から現在の音源状態情報を予測する状態予測部と
前記時間差算出部が算出した時間差と前記状態予測部が予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する状態更新部とを備える
ことを特徴とする音源位置推定装置。
【請求項２】
前記状態更新部は、
前記誤差に基づいてカルマンゲインを算出し、算出したカルマンゲインに前記誤差を乗ずることを特徴とする請求項１に記載の音源位置推定装置。
【請求項３】
前記音源状態情報は、前記音声信号を前記信号入力部に供給する収音部の位置を含むことを特徴とする請求項１又は２に記載の音源位置推定装置。
【請求項４】
前記収音部の位置の変化に基づいて、前記音源位置の変化が収束したか否かを判断する収束判定部を備えること
を特徴とする請求項３に記載の音源位置推定装置。
【請求項５】
前記複数のチャネルの入力信号を、予め定めた音源位置の評価点から前記複数のチャネルの各々に対応する収音部の位置までの位相で補償した信号を加算して得られる評価値を最大にする評価点を定め、定めた評価点と前記状態更新部が推定した音源状態情報が表す音源位置までの距離に基づいて、前記音源位置の変化が収束したか否かを判断する収束判定部を備えること
を特徴とする請求項３に記載の音源位置推定装置。
【請求項６】
前記収束判定部は、
前記評価点を、遅延和ビームフォーミング法を用いて定め、定めた評価点と請求項２に記載の状態更新部が推定した音源状態情報が表す音源位置までの距離に基づいて、前記音源位置の変化が収束したか否かを判断することを特徴とする請求項５に記載の音源位置推定装置。
【請求項７】
音源位置推定装置における方法において、
前記音源位置推定装置が、複数のチャネルの音声信号を入力する過程と、
前記音源位置推定装置が、チャネル間の音声信号の時間差を算出する過程と、
前記音源位置推定装置が、音源位置を含む音源状態情報である過去の音源状態情報から現在の前記音源状態情報を予測する過程と
前記音源位置推定装置が、前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する過程を有すること
を特徴とする音源位置推定方法。
【請求項８】
音源位置推定装置のコンピュータに、
複数のチャネルの音声信号を入力する手順、
チャネル間の音声信号の時間差を算出する手順、
音源位置を含む音源状態情報である過去の音源状態情報を予測する手順、
前記算出した時間差と前記予測した前記音源状態情報に基づく時間差との間の誤差を減少させるように前記音源状態情報を推定する手順
を実行させるための音源位置推定プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１８】
【図１９】
【図１５】
【図１６】
【図１７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１８】
【図１９】
【図１５】
【図１６】
【図１７】
【公開番号】特開２０１２−１６１０７１（Ｐ２０１２−１６１０７１Ａ）
【公開日】平成２４年８月２３日（２０１２．８．２３）
【国際特許分類】
【出願番号】特願２０１１−２７１７３０（Ｐ２０１１−２７１７３０）
【出願日】平成２３年１２月１２日（２０１１．１２．１２）
【出願人】（０００００５３２６）本田技研工業株式会社 (23,863)
【Ｆターム（参考）】
【公開日】平成２４年８月２３日（２０１２．８．２３）
【国際特許分類】
【出願日】平成２３年１２月１２日（２０１１．１２．１２）
【出願人】（０００００５３２６）本田技研工業株式会社 (23,863)
【Ｆターム（参考）】
[ Back to top ]