パラメタ推定装置、音源分離装置、方向推定装置、それらの方法、プログラム
【課題】空間的エイリアジングの問題が生じても、適切な音源分離を行うための確率分布のパラメタを推定する。
【解決手段】複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成し（３０、Ｓ１０２）、前記周波数観測信号の前記収音手段間の位相差を計算し（３４、Ｓ１０４）、前記位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する（２００、Ｓ１０６）。
【解決手段】複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成し（３０、Ｓ１０２）、前記周波数観測信号の前記収音手段間の位相差を計算し（３４、Ｓ１０４）、前記位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する（２００、Ｓ１０６）。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は複数の音源信号が混合された観測信号を例えば分離するために用いる確率分布モデルのパラメタを求めるパラメタ推定装置、求められたパラメタを用いた音源分離装置、方向推定装置、それらの方法、プログラムに関する。
【背景技術】
【０００２】
従来技術の音源分離装置５００の機能構成例を図１に示す。この従来技術の詳細は、例えば、非特許文献１に記載されている。Ｎ個の音源１０ｎ（ｎ＝１、．．．、Ｎ）それぞれから同時に発せられる音源信号ｓｎ（ｔ）を２個の収音手段である第１収音手段２１、第２収音手段２２で、ある収録時間内（例えば、５秒間）に観測する。この観測状況を状況Ｘとする。この状況Ｘの下、第１収音手段２１で収音された観測信号をｘ１（ｔ）とし、第２収音手段２２で収音された観測信号をｘ２（ｔ）とし、観測された観測信号をＸ（ｔ）＝［ｘ１（ｔ）、ｘ２（ｔ）］Ｔとする。収音手段とは例えばマイクロホンのことであり、「Ｔ」は行列の転置を表し、ｔを離散時刻とし、ｔ＝０、．．．、Ｔ−１とする。周波数領域変換部３０は、観測信号Ｘ（ｔ）を周波数領域に変換することで、周波数毎の時系列信号である観測信号ベクトルＸ（ｆ，ｔ）＝［ｘ１（ｆ，ｔ），ｘ２（ｆ，ｔ）］Ｔに変換する。Ｘ（ｆ，ｔ）を以下では、周波数観測信号ベクトルという。周波数領域への変換は、例えば短時間フーリエ変換を用いれば良い。ｆは周波数を示し、ｆ∈｛０，ｆｓ／Ｆ，．．．，（Ｆ−１）ｆｓ／Ｆ｝であり、ｆｓはサンプリング周波数を示し、Ｆは周波数帯域の数を示す。
ここで、周波数観測信号ベクトルは、以下の式（１）で表されると仮定する。
【０００３】
【数１】
【０００４】
ただし、ｊ（＝１，２）は収音手段のインデックスであり、ｊ＝１が第１収音手段２１について示し、ｊ＝２が第２収音手段２２について示し、ｈｊｎ（ｆ）は音源１０ｎから収音手段（第１収音手段２１、第２収音手段２２までの周波数応答を示し、ｓｎ（ｆ，ｔ）は音源信号ｓｎ（ｔ）を周波数領域で表現した信号であり、以後では、周波数音源信号ｓｎ（ｆ，ｔ）という。
【０００５】
音源分離を行うために、音源がスパースである、すなわち、周波数音源信号ｓｎ（ｆ，ｔ）はまれにしか大きな値をとらず、ｓｎ（ｆ，ｔ）とｓｎ’（ｆ，ｔ）（ｎ≠ｎ’）は各時間周波数スロット（ｆ，ｔ）で互いに重ならない、ということを仮定する。これは、互いに異なる音声信号などで確認される性質である。これを仮定すると前記式（１）は、
ｘｊ（ｆ，ｔ）≒ｈｊｎ（ｆ）ｓｎ（ｆ，ｔ） （２）
となる。ここで、周波数音源信号ｓｎ（ｆ，ｔ）は時間周波数スロット（ｆ，ｔ）にて支配的な音源信号である。
【０００６】
また、上記式（１）における周波数応答ｈｊｎ（ｆ）が、音源信号ｓｎ（ｔ）が収音手段に到達する間に、部屋の影響を受けない、つまり、無響モデルであると仮定すると、以下の式（３）で表される。
【０００７】
【数２】
【０００８】
ここで、μｎ（ｆ，ｔ）は音源１０ｎが第１収音手段２１と第２収音手段２２に到達する時間差である。また、図２に第１収音手段２１、第２収音手段２２から見た音源の方向角度φｎ（ｆ，ｔ）について示す。図２からも明らかなように、
μｎ（ｆ，ｔ）＝Ｄｃｏｓφｎ（ｆ，ｔ）／ｖ （４）
となる。ただし、Ｄは第１収音手段２１、第２収音手段２２との距離であり、ｖは音速である。
【０００９】
次に、位相差計算部３４は、第１収音手段２１で収音された周波数観測信号ｘ１（ｆ，ｔ）、第２収音手段２２で収音された周波数観測信号ｘ２（ｆ，ｔ）の位相差ｄｆｔを以下の式（５）により計算する。
ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］ （５）
ここで、前記式（２）（３）より音源１０ｎが支配的な時間周波数スロット（ｆ，ｔ）では以下の式（６）、（７）が成り立つ。
【００１０】
【数３】
【００１１】
次に、分類部７０は、位相差ｄｆｔを用いてクラスタリング処理を行う。まず、分類部７０は、以下の式（８）により位相差ｄｆｔを周波数ｆで正規化する。
【００１２】
【数４】
【００１３】
分類部７０は、式（８）についてクラスタリングを行う。スパース性の仮定より定められる前記式（２）より、音源１０ｎだけが支配的な時間周波数（ｆ，ｔ）では時間差μｎ（ｆ，ｔ）が求まっており、音源１０ｎ’だけが支配的な時間周波数（ｆ，ｔ）では時間差μｎ’（ｆ，ｔ）が求まっているので、音源１０ｎの方向角度φｎ（ｆ，ｔ）と音源１０ｎ’の方向角度φｎ７（ｆ，ｔ）とが異なっていれば、両者を分類（クラスタリング）できる。そして、それぞれのクラスタが各音源成分に対応する。
【００１４】
音源分離部７２では、音源１０ｎに対応するクラスタを形成している時間周波数（ｆ，ｔ）では１をとり、それ以外の時間周波数では０をとるマスクＬｎ（ｆ，ｔ）を生成する。このマスクＬｎ（ｆ，ｔ）を全ての音源１０ｎについて生成する。
【００１５】
そして、音源分離部７２は、マスクＬｎ（ｆ，ｔ）を観測信号の１つ（ここではｘ１（ｆ，ｔ））に乗算することで、つまり以下の式（１０）を演算することで、分離信号ｙｎ（ｆ，ｔ）を得る。
ｙｎ（ｆ，ｔ）＝ｘ１（ｆ，ｔ）Ｌｎ（ｆ，ｔ） （１０）
そして、時間領域変換部７４は、分離信号ｙｎ（ｆ，ｔ）を時間領域に変換することで、時間領域分離信号ｙｎ（ｔ）を得る。
【先行技術文献】
【非特許文献】
【００１６】
【非特許文献１】S．Araki、H．Sawada、R．Mukai、and S．Makino、”Underdetermined blind sparse source separation for arbitrarily arranged multiple sensors、”Signal Processing、vol．87、pp．1833−1847、Feb．2007．
【発明の概要】
【発明が解決しようとする課題】
【００１７】
従来技術では、第１収音手段２１および第２収音手段２２との間の収音手段間隔Ｄが狭く、音源信号の周波数ｆが低い場合には、非特許文献１に記載されているように良好に動作する。しかし、収音手段間隔Ｄが広い場合や、周波数ｆが高い場合にはクラスタリング処理がうまく動作しないため、従来技術では適切に音源分離を行うことは難しかった。図３Ａ〜Ｆに音源１０ｎが２つの場合、つまり、ｎ＝１，２となる場合に、収音手段間隔Ｄおよびサンプリング周波数ｆｓを変えた場合のクラスタリング処理の模様を示す。図３Ａ〜Ｆにおいて、実線を音源１０１の音源信号を示し、破線を音源１０２の音源信号を示し、横軸はそれぞれの音源信号ｓ１（ｆ，ｔ）、音源信号ｓ２（ｆ，ｔ）の周波数ｆを示す。また、図３Ａ、Ｂの縦軸はそれぞれ、サンプリング周波数ｆｓを８ｋＨｚ、収音手段間隔Ｄを４ｃｍとした（Ｄが狭い）場合の、位相差ｄｆｔ、正規化後の位相差ｄｆｔ／２πｆを示し、図３Ｃ、Ｄの縦軸はそれぞれ、サンプリング周波数ｆｓを８ｋＨｚ、収音手段間隔Ｄを１０ｃｍとした（Ｄが広い）場合の、位相差ｄｆｔ、正規化後の位相差ｄｆｔ／２πｆを示し、図３Ｅ、Ｆの縦軸はそれぞれ、サンプリング周波数ｆｓを１６ｋＨｚ、収音手段間隔Ｄを４ｃｍとした（Ｄが狭い）場合の、位相差ｄｆｔ、正規化後の位相差ｄｆｔ／２πｆを示す。この図３を用いて、従来技術の問題点を説明する。
【００１８】
図３Ｂに示すように、収音手段間隔Ｄが狭く、周波数ｆが低い場合には、全ての周波数ｆで、位相差ｄｆｔ／２πｆが同じ値をとっており、これをクラスタリングすると２個のクラスタが形成され、適切に音源分離を行うことができる。
【００１９】
一方、図３Ｄに示すように収音手段間隔Ｄが狭い場合には、ｆ＞３０００Ｈｚの範囲では、位相差ｄｆｔ／２πｆが一定でなくなる。このようなｄｆｔ／２πｆをクラスタリングしても、音源ごとのクラスタは形成されず、適切に音源分離を行うことができない。この理由は、前記式（７）の右辺が収音手段間隔Ｄが大きいために±πの範囲を超えてしまうにも関わらず、式（６）のａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］の演算が−πからπの範囲の値を算出するためである。これを具体的に説明すると、ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］のａｒｇ演算が−πからπの間の値しか返せないため
２πｆμｎ（ｆ，ｔ）＝ｄｆｔ＋２πｋ （１１）
＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］＋２πｋ
（１２）
の関係を満たす−π＜ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］＜πとして返されてしまうためである。これを「空間的エイリアジングの問題」や、「２πｋの不定性がある」という。ここで、ｋはμｎ（ｆ，ｔ）の値が既知であれば、一意に決まる整数であり、不定性係数ｋという。ただし、一般的には、μｎ（ｆ，ｔ）の値は未知であるため、不定性係数ｋは推定すべき整数である。
【００２０】
この空間的エイリアジングの問題は、図３Ｅ、図３Ｆ（特に、図３Ｆ）に示すように、収音手段間隔Ｄが狭くても（Ｄ＝４ｃｍ）、周波数（図３Ｆの例では、周波数ｆが高くなる（周波数ｆが６０００Ｈｚ以上）場合には、生じる問題である。
【００２１】
このように従来技術では、空間的エイリアジングの問題が生じるような、収音手段間隔Ｄが狭く、周波数ｆが高くなる場合には、適切なクラスタリングができないため、良好な音源分離ができなかった。
【課題を解決するための手段】
【００２２】
上記の課題を解決するために、この出願のパラメタ推定装置は、周波数領域変換部と、位相差計算部と、推定部と、を有する。周波数領域変換部は、複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成する。位相差計算部は、周波数観測信号の収音手段間の位相差を計算する。推定部は、位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する。
【発明の効果】
【００２３】
本発明のパラメタ推定装置により推定されたパラメタθを用いて、音源分離を行うことで、空間的エイリアジングの問題が生じる場合においても、良好に音源分離を行うことができる。
【図面の簡単な説明】
【００２４】
【図１】従来の音源分離装置の機能構成例を示した図。
【図２】収音手段間距離Ｄなどを説明した図。
【図３】従来の音源分離装置のクラスタリングの様子を示した図。
【図４】本実施例のパラメタ推定装置の機能構成例を示した図。
【図５】本実施例のパラメタ推定装置の処理フローを示した図。
【図６】位相差の分布に正規分布モデルを当てはめることを示した図。
【図７】本実施例の推定部の機能構成例を示した図。
【図８】本実施例の推定部の処理フローを示した図。
【図９】本実施例の音源分離装置の機能構成例を示した図。
【図１０】本実施例の方向推定装置の機能構成例を示した図。
【図１１】実験を行った部屋を示した図
【図１２】実験結果を示した図。
【発明を実施するための形態】
【００２５】
以下に、発明を実施するための最良の形態を示す。なお、同じ機能を持つ構成部や同じ処理を行う過程には同じ番号を付し、重複説明を省略する。上述の課題を解決する音源分離装置については実施例２で説明し、実施例１では当該音源分離装置に用いられるパラメタ推定装置について説明する。実施例３では、当該パラメタ推定装置を用いた方向推定装置について説明する。
【実施例１】
【００２６】
図４に実施例１のパラメタ推定装置１００の機能構成例を示し、図５に処理フローを示す。このパラメタ推定装置１００を用いた音源分離装置３００の機能構成例を図９および実施例２に示す。また、前記状況Ｘの下、第１収音手段２１で収音された観測信号ｘ１（ｔ）および第２収音手段２２で収音された観測信号ｘ２（ｔ）は、周波数領域変換部３０に入力される。周波数領域変換部３０は、観測信号Ｘ（ｔ）＝［ｘ１（ｔ），ｘ２（ｔ）］Ｔを周波数領域に変換することで周波数観測信号（ベクトル）Ｘ（ｆ，ｔ）＝［ｘ１（ｆ，ｔ），ｘ２（ｆ，ｔ）］Ｔを生成する（ステップＳ１０２）。
【００２７】
周波数観測信号Ｘ（ｆ，ｔ）は入力作成部３３に入力される。入力作成部３３は、パワー推定部３２および位相差計算部３４とで構成される。パワー推定部３２は以下の式（１３）により時間周波数スロット（ｆ，ｔ）ごとにパワーａｆｔを計算する。
ａｆｔ＝│ｘｊ（ｆ，ｔ）│ ただしｊ＝１，２ （１３）
位相差計算部３４は時間周波数スロット（ｆ，ｔ）ごとに、第１収音手段２１と第２収音手段２２との間の位相差ｄｆｔを計算する（ステップＳ１０４）。具体的には、前記式（５）により計算する。念のため式（５）を以下に示す。
ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］ （５）
【００２８】
次に、推定部２００は、位相差ｄｆｔの分布（ヒストグラム）に当てはまり、かつ、周波数依存性のある確率分布モデルのパラメタθを推定する（ステップＳ１０６）。換言すれば、各音源に対応するクラスタを確率分布モデル（例えば、正規分布）で当てはめる（近似するまたはモデル化する）場合のパラメタ推定を行う。「周波数依存性のある確率分布モデル」については後ほど詳細に述べる。また、推定した確率分布モデルのパラメタθは、例えば、以下で説明する音源分離処理（実施例２）や方向推定処理（実施例３）で用いる。以下の説明では、確率分布モデルとして正規分布を用いた場合の説明を行うが、用いる確率分布モデルは、正規分布モデルに限らない。
【００２９】
ここでは、まず、各音源に対応するクラスタを正規分布モデルで当てはめる例の説明を行う。まず、音源数が１（クラスタが１個）の場合を考える。このとき、位相差ｄｆｔの分布を書くと、図６Ａに示すようになる。図６Ａからも理解できるように、１つの山（＝クラスタ）ができる。このクラスタに以下の式（１４）で示される正規分布モデルを当てはめる。
【００３０】
【数５】
【００３１】
ここで、μは正規分布のパラメタθのうち「平均」を示し、σ２は正規分布のパラメタθのうち「分散」を示す。特に、μは前記式（３）で説明した、音源からの音が第１収音手段２１、第２収音手段２２に到達する時間の時間差であり、後述する正規分布モデルの当てはめを行うことで推定する値である。また上述の通り、ｋは不定性係数であり、後述する正規分布モデルの当てはめを行うことで推定する必要がある。また。Ｋｆは不定性係数ｋの最大値であり、周波数ｆによって異なる値をとりうる。具体的には、例えば、以下の式（１５）で表される。
【００３２】
【数６】
【００３３】
また、全ての周波数ｆで同じ値（ただし、十分大きな値）を用いても良い。
次に、音源数が２以上（すなわちクラスタが２つ以上）の場合を考える。図６Ｂに音源数が２の場合の位相差ｄｆｔの分布を示す。図６Ｂからも理解されるように、２つの山（＝クラスタ）ができることがわかる。このクラスタにそれぞれ１つの正規分布モデルを当てはめることを考える。つまり、分布全体を以下の式（１６）で示される混合正規分布モデルＧ（ＧＭＭ）で当てはめることを考える。
【００３４】
【数７】
【００３５】
ここで、ｍは正規分布のインデックスを示しつまり、（ｍ＝１，．．．，Ｍ）であり、図６Ｂの例では、Ｍ＝２となり、μｍはｍ番目の正規分布の平均を示し、σｍ２はｍ番目の正規分布の分散を示し、特に、μｍは音源からの音が第１収音手段２１、第２収音手段２２に到達する時間の時間差の後述する正規分布モデルの当てはめを行うことで推定する値である。また、αｍは、ｍ番目の正規分布の混合重みであり、Σｍ＝１Ｍαｍ＝１であり、０≦αｍ≦１である。また、混合正規分布を構成するＭ個の正規分布をΨｍ（ｍ＝１，．．．，Ｍ）とし、正規分布の数Ｍ（以下、「混合数Ｍ」という。）は、音源数Ｎが既知の場合には、Ｍ＝Ｎとすることができる。また、音源数が未知の場合には、Ｍを十分大きな数（例えば、Ｍ＝１０）を用いれば良い。本実施例では、複数音源の分離などの処理を行うために、位相差ｄｆｔに当てはまる前記式（１６）に示す混合正規分布を用いてパラメタθを推定する。以下、推定部２００のパラメタθの推定処理について詳細に説明する。
【００３６】
図７に推定部２００の機能構成例を示し、図８に推定部２００の処理フローを示す。パラメタθを混合正規分布の平均μ、分散σ２、混合重みαをまとめたものを示し、つまり、θ＝（μｍ，σｍ，αｍ）＝（μ１，σ１，α１，．．．，μｍ，σｍ，αｍ，．．．，μＭ，σＭ，αＭ）となる。また、ｒを更新回数とし、θに更新回数の概念を付与したもの、つまり、ｒ回更新したθをθｒとすると、θｒ＝（μ１ｒ，σ１ｒ，α１ｒ，．．．，μｍｒ，σｍｒ，αｍｒ，．．．，μＭｒ，σＭｒ，αＭｒ）となる。また、記憶部１６には予め用いる正規分布モデルのモデル数Ｍと混合正規分布モデルの各パラメタの初期値θ０が記憶されている。事前分布情報保持部１１０には、ハイパーパラメタω（後述する）、重みパラメタｃが保持されている。
【００３７】
推定部２００には、パワー推定部３２よりのパワーａｆｔが重み係数ａｆｔとして入力され、位相差計算部３４よりの位相差ｄｆｔも入力される。または全てのａｆｔについてａｆｔ＝１としても良く、この場合は、パワー推定部３２はなくてもよい。また、重み係数ａｆｔを各時間周波数（ｆ，ｔ）における観測信号のパワーや信号の瞬時的ＳＮ比などとすることも出来る。従って、パワー推定部３２を設けない場合であっても、重み係数ａｆｔを入力部３５から入力することができる。
【００３８】
まず、初期設定として、ｒ＝０（つまり更新回数が０）、ｒ＝０のときの混合正規分布のパラメタθ０の値、用いる正規分布のモデル数Ｍ、不定性係数ｋの範囲であるＫｆ、更新回数閾値Ｒまたは差閾値Δ（後述する）を設定する（ステップＳ２）。更新回数閾値Ｒまたは差閾値Δは、後述する収束判定処理の際に用いられる。これらの初期設定は、入力部３５からの入力により行われる。
【００３９】
事後確率計算部１２は、位相差ｄｆｔと、現在の確率分布モデル（混合正規分布モデル）のパラメタθｒ（＝（μｍｒ，σｍｒ，αｍｒ ｍ＝１，．．．，Ｍ））から、Ｍ個の正規分布Ψｍごとに事後確率ｐ（ｍ，ｋ│ｄｆｔ，θｒ）を計算する（ステップＳ６）。またパラメタ保持部１８には、現在の混合正規分布のθｒが保持されている。事後確率計算部１２は具体的には例えば、以下の式（１７）（１８）により計算する。
【００４０】
【数８】
【００４１】
次に、更新部１４は、位相差ｄｆｔと事後確率ｐ（ｍ，ｋ│ｄｎ，θｒ）を用いて、現在の混合正規分布の各パラメタθｒを更新する（ステップＳ８）。以下、更新処理について詳細に説明する。更新部１４は更新処理の際に、ハイパーパラメタω、重みパラメタｃを事前分布情報保持部１１０から取り出す。この実施例では、音源数Ｎが未知、つまり、正規分布の数Ｍが未知の場合であっても、適用可能にするために、正規分布のパラメタθの混合重みαｍに適切な事前分布を与え、例えばＥＭアルゴリズムにてパラメタθの更新処理を行う。この実施例１では、混合重みαｍの事前分布として、ディリクレ分布を考える。ディリクレ分布の詳細は、参考文献１である「Ｃ．Ｍ．ビショップ著（元田、栗田他訳） 「パターン認識と機械学習（上）」、シュプリンガー・ジャパン２００７年 ｐ．７４−ｐ．７７」等に記載されている。ディリクレ分布は例えば以下の式（１９）で表される。
【００４２】
【数９】
【００４３】
ここで、αは混合重み行列であり、α＝｛α１，．．．，αｍ，．．．，αＭ｝で表され、ΣｍＭαｍ＝１、０≦αｍ≦１という条件を満たす。これは混合正規分布のパラメタである混合重みの条件と同じであることに注意されたい。またβ（ω）は正規化項（ベータ分布）であり、ここで、ハイパーパラメタωを１より小さい正の値（例えば、０．９）に設定すると、αｍのごく少数のみが十分に大きな値を持ち、残りは０に近い値をとるようになる。求められたαｍを前記式（１６）で用いられている混合重みαｍの事前分布として用いることで、混合正規分布モデルＧのうちの少数の正規分布のみに十分大きな混合重みがかかり、その他の正規分布モデルの混合重みは０に近くなる。結果として、なるべく少数の正規分布による当てはめが可能である。従って、１つのクラスタに複数の正規分布が当てはまるような現象を防ぐことが出来、音源数未知数の場合でも、それぞれのクラスタに１つずつ正規分布を当てはめることができる。
【００４４】
次に、この事前分布を含みながら、パラメタ更新を行うためのＥＭアルゴリズムを導出する。ここで、正規分布のインデックスｍと不定性係数ｋは位相差ｄｆｔから推定すべき変数であるため、ＥＭアルゴリズムにおける隠れ変数として扱う。これにより不定性係数ｋは隠に自動推定されるため、位相差ｄｆｔにおける２πｋの不定性を自動的に扱うことが可能になる。まず、最尤推定のためのコスト関数Ｌ（θ）は次のように与えられる。
【００４５】
【数１０】
【００４６】
また、重みパラメタｃは、式（２２）の第１項と第２項の重みをコントロールするパラメタであり、上述のように、事前分布情報保持部１１０に保持される。
【００４７】
【数１１】
【００４８】
となる。ここで、式（２４）のＥ［Ｈ］は式Ｈの期待値を示し、式（２５）中のｐ（ｍ，ｋ│ｄｆｔ，θｒ）は式（１８）で表される事後確率分布である。
【００４９】
【数１２】
【００５０】
図７中の更新部１４中の平均更新手段１４２が式（２６）より現在の平均μｍｒを更新することで更新後の平均μｍｒ＋１を出力する。分散更新手段１４４が式（２７）より分散（σｍ２）ｒを更新することで更新後の分散（σｍ２）ｒ＋１を出力する。混合重み更新手段１４６が式（２８）により混合重みαｍを更新することで更新後の混合重みαｍｒ＋１を出力する。パラメタ算出手段が、更新後の平均μｍｒ＋１、分散（σｍ２）ｒ＋１、混合重みαｍｒ＋１についての更新後のパラメタθｒ＋１を算出する（ステップＳ８）。
【００５１】
各パラメタの更新処理が数回行われ（ステップＳ４）、更新部１４内の収束判定手段１５０は、更新されたθｒ＋１に対して、予め定められた規則により、各パラメタ値が収束しているか否かの収束判定を行う（ステップＳ１０）。各パラメタ値が収束していると判断した場合には、更新されたパラメタθｒ＋１を出力する。また、各パラメタ値が収束していないと判断した場合には、更新されたパラメタθｒ＋１を現在の確率分布モデルの平均、分散、混合重みとしてパラメタ保持部１８に保持させる。を繰り返す。そして、収束判定手段１５０が、各パラメタ値が収束していると判断するまで、ステップＳ４〜ステップＳ１０の処理（平均更新手段１４２、分散更新手段１４４、混合重み更新手段１４６の処理）を繰り返す。
【００５２】
ここで収束判定に用いる予め定められた規則の例を説明する。更新回数閾値Ｒを用いる例を説明すると、更新部１４内のカウント手段（図示せず）は更新回数ｒをカウントし、更新回数ｒが更新回数閾値Ｒ（例えば３０）を超えた場合には、十分更新しており、収束していると判断して、パラメタ算出手段１４８は、更新後のパラメタθを出力する。また、差閾値Δを用いる例を説明すると、以下の式（２９）の式を満たす場合には、収束していると判断して、パラメタ算出手段１４８は、更新後のパラメタθを出力する。
│Ｑ（θ│θｒ＋１）−Ｑ（θ│θｒ）│＜Δ （２９）
このようにして、推定部２００は、各音源に対応するクラスタを正規分布モデルで当てはめたときのパラメタθ＝（μｍ，σｍ，αｍ）を出力する。
【００５３】
この実施例１では混合重みαｍのみに事前分布を導入したが、各ガウス分布の平均μｍと分散σｍに対しても事前分布を導入することで、より精度の高い混合正規分布の当てはめを実現できる。また、各ガウス分布の各パラメタである平均μｍ、分散σｍ、混合重みαｍ、に事前分布を導入した場合の当てはめ処理には、ＥＭアルゴリズムの他、不定性係数ｋを隠に自動推定するアルゴリズムであれば、何でも良い。これらの拡張は当業者であれば、上記参考文献１などを参照すれば、容易に実現できるため、ここでは省略する。
【００５４】
また、音源数Ｎが既知であり、混合数Ｍ＝Ｎとできれば、前記式（１９）で示されるディリクレ分布を用いる必要はない。すなわち、この場合には、前記式（２８）でハイパーパラメタω＝１とすればよい。例えば、ユーザが入力部３５から音源数Ｎが既知であるか未知であるかを示す情報である音源数情報を入力し、音源数情報が既知である旨の情報であれば、ディリクレ分布を用いず、音源数情報が未知である旨の情報であれば、ディリクレ分布を用いれば良い。
【００５５】
図３で説明したように、空間的エイリアジングの影響は周波数毎に異なる、つまり、周波数依存性があるといえる。よって、確率分布モデルでモデル化する場合には、この周波数依存性を考えることが必要である。
【００５６】
これを式（１４）（１６）に示す混合正規分布モデルについて検討する。特にｅｘｐ項の分子の式「−（ｄｆｔ＋２πｋ−２πｆμｍ）２」を検討すると、まず、２πｋについては、上述したように、パラメタ推定は例えばＥＭアルゴリズムを用いて行うが、この際、不定性係数ｋはＥＭアルゴリズムにおける隠れ変数として扱うことができる。よって不定性係数ｋはＥＭアルゴリズムでデータから（確率的に）自動推定されるため、位相差ｄｆｔの２πｋの不定性を自動的に扱うことができる。また、ｄｆｔについては式（１１）、（１２）の通り、周波数領域の観測信号ｘ１（ｆ，ｔ）、ｘ２（ｆ，ｔ）の位相差であるため、周波数依存性のある値であるといえる。２πｆμｍについても周波数ｆを含んでいるので、周波数依存性のある値であるといえる。つまり、式（１４）（１６）に示す混合正規分布モデルは、周波数依存性のあるモデルであるといえる。また、推定部２００で用いる確率分布モデルは、周波数依存性があれば、式（１４）（１６）に限らない。
【００５７】
また、式（１１）において、実測値（観測された値）は式（１１）の右辺に示すｄｆｔ＋２πｋであり、当てはめ後（モデル化後）の値は、式（１１）左辺の２πｆμｎ（ｆ，ｔ）である。式（１６）のｅｘｐ項の分子「（ｄｆｔ＋２πｋ−２πｆμｎ）２」は、実数値とモデル化後の値との二乗誤差（モデル化誤差）を意味するとも捉えることができる。観測された位相差ｄｆｔのヒストグラムに当てはまる正規分布は、このモデル化誤差を最も小さくする正規分布ということになり、前記式（２６）〜（２８）を用いて、当該正規分布を求めている。
【００５８】
従来技術では、前記式（８）に示すように、ｄｆｔ／２πｆのように、周波数正規化した量をクラスタリングしていた。そのため、図３Ｄ、Ｆに示すように、正しくクラスタリングできなかった。しかし、本実施例では、周波数依存性のあるｄｆｔの分布をそのまま前記式（１６）で示す混合正規分布に当てはめる。具体的には、本実施例のパラメタ推定装置１００は、前記式（１１）（１２）で与えられる空間的エイジアリングの問題（２πｋの不定性）を陽に定式化し、位相差ｄｆｔのヒストグラムに当てはまり、周波数依存性のある混合正規分布（式（１４）や式（１６））のパラメタを求める。この求められたパラメタを用いて例えば音源分離を行うと、周波数依存性の高い空間的エイリアジングの問題を扱うことが可能となり、適切な音源分離などを行うことができる。
【００５９】
また、このパラメタ推定装置１００で求められたパラメタθは、実施例２で説明する音源分離処理や、実施例３で説明する音源方向推定処理のほか、様々な観測信号処理に用いられる。
【実施例２】
【００６０】
この実施例２では、実施例１で説明したパラメタ推定装置１００により推定されたパラメタθを用いて、音源分離を行う音源分離装置について説明する。図９に実施例２の音源分離装置３００の機能構成例を示す。またパラメタ推定装置で推定されたパラメタθを決定後パラメタという。
【００６１】
パラメタ推定装置１００よりの決定後パラメタθは有効音源推定部４０に入力される。有効音源推定部４０は、音源に該当する確率分布モデルを示す音源該当情報を求める。有効音源推定部４０による音源該当情報の生成手法は以下の３つの手法により求められる。ここで、音源該当情報とは、例えば、音源に該当する確率分布モデルのインデックスｍをいう。
【００６２】
まず、第１の手法として、音源数Ｎが既知であり、パラメタ推定装置１００で用いられる混合正規分布の混合数Ｍ＝Ｎとしている場合には、混合正規分布を構成する全ての正規分布が音源に該当するので、全ての正規分布のインデックスｍ（＝１，．．．，Ｍ）を出力する。また音源数Ｎが未知の場合には、下記の第２手法、第３の手法により求められる。
【００６３】
第２の手法として、有効音源推定部４０は、決定後パラメタθの混合重みαｍ（ｍ＝１，．．．，Ｍ）のうち、混合重みが予め定められた第１閾値ε１（例えば１０−６）よりも大きな値である正規分布を音源に該当する正規分布と判断して、当該正規分布（以下、「音源該当正規分布」という。）のインデックスｍ’を出力する。何故なら、パラメタ推定装置１００の演算が十分収束している場合には、決定後パラメタθの中の混合重みαｍのうち十分大きな値を持つ個数は位相差ｄｆｔのヒストグラム中の分布の山の個数と等しくなるからである。また、音源該当正規分布の数をＭ’（つまり、ｍ’＝１，．．．，Ｍ’）とする。
【００６４】
また、第３の手法として、第２の手法においてパラメタ推定装置１００の演算が十分に収束していない場合は、有効音源推定部４０は、次のような推定処理を行うことが好ましい。有効音源推定部４０は、混合重みαｍが第１閾値ε１よりも大きく、かつ分散σｍ２が予め定められた第２閾値ε２（例えばπ／５）よりも小さい正規分布を音源に該当する正規分布と判断して、音源該当正規分布のインデックスｍ’を音源該当情報として出力する。
【００６５】
また、音源数Ｎが既知であるか未知であるかについての情報である音源数情報（つまり、第１の手法を用いるか、または第２、第３の手法を用いるか）は、ユーザに入力部４７から入力させればよい。
【００６６】
次に、マスク作成部４２は、音源該当情報（ここでは、音源該当正規分布のインデックスｍ’）が示す確率分布モデル（正規分布モデル）を周辺化することでマスクΩｍ’（ｆ，ｔ）を作成する。マスクΩｍ’（ｆ，ｔ）は、各音源該当正規分布Ψｍ’ごとに、かつ、各時間周波数スロット（ｆ，ｔ）ごとに求められる。具体的には、Ｍ’個の音源該当正規分布に関する事後確率ｐ（ｍ’，ｋ│ｄｆｔ，θｒ）を周辺化することで、周辺化事後確率ｐ（ｍ’│ｄｎ，θｒ）（＝マスクΩｍ’（ｆ，ｔ））を求める。
【００６７】
【数１３】
【００６８】
式（３０）中のｐ（ｍ’，ｋ│ｄｆｔ，θ）については、マスク作成部４２がパラメタ推定装置１００中の事後確率計算部１２から前記式（１８）の結果を抽出すればよい。
【００６９】
そして、分離部４４は、周波数観測信号にマスクΩｍ’（ｆ，ｔ）を乗算することで、分離信号ｙｍ’（ｆ，ｔ）を求める。具体的には以下の式（３１）により分離信号ｙｍ’（ｆ，ｔ）を求める。例えば、分離部４４は、マスクΩｍ’（ｆ，ｔ）を観測信号の１つ（ここでは、周波数観測信号ｘ１（ｆ，ｔ））に乗算し、分離信号ｙｍ’（ｆ，ｔ）を得る。つまり、以下の式（３１）により求められる。
ｙｍ’（ｆ，ｔ）＝ｘ１（ｆ，ｔ）Ωｍ’（ｆ，ｔ） （３１）
分離部４４よりの分離信号ｙｍ’（ｆ，ｔ）は、時間領域変換部４６に入力される。そして、時間領域変換部４６は、分離信号ｙｍ’（ｆ，ｔ）を時間領域に変換して時間領域分離信号ｙｍ’（ｔ）を求め、出力する。
【００７０】
ここで、従来の音源分離装置５００（図１参照）と、実施例２の音源分離装置３００（図９参照）の対応関係を以下に示す。
音源分離装置５００の音源分離部７２
→音源分離装置３００のマスク作成部４２と分離部４４を統合したもの
音源分離装置５００の周波数領域変換部３０と位相差計算部３４と分類部７０を統合
したもの→音源分離装置３００のパラメタ推定装置１００
また、従来の音源分離装置５００では、音源数が既知である場合が多いため、音源分離装置５００は、音源分離装置３００の有効音源推定部４０に対応するものを有していなかった。
【００７１】
このように、実施例１で説明したパラメタ推定装置１００により推定されたパラメタは周波数依存性の高い空間的エイリアジングの問題を扱うことができるパラメタである。この実施例２の音源分離装置３００は、パラメタ推定装置１００で推定されたパラメタθを用いて音源分離を行うことから、空間的エイリアジングの問題が生じるような条件下であっても、適切な音源分離を行うことができる。
【実施例３】
【００７２】
実施例１で説明したパラメタ推定装置１００により推定されたパラメタθを用いて、音源の方向を推定することもできる。この実施例３では、実施例１で説明したパラメタ推定装置１００を用いた、音源の方向を推定する方向推定装置４００を説明する。図１０に方向推定装置４００の機能構成例を示す。方向推定部６０は、有効音源推定部５０と方向出力部５２とで構成されている。
【００７３】
パラメタ推定装置１００よりの決定後パラメタθは、有効音源推定部５０および方向出力部５２に入力される。有効音源推定部５０は、実施例２で説明したように、音源に該当する確率分布モデルを示す音源該当情報（例えば、音源該当正規分布のインデックスｍ’）を求める。有効音源推定部５０の処理が終了すると、方向出力部５２は、音源該当方向情報分布モデルのインデックスｍ’｛ｍ’＝１，．．．，Ｍ’｝に対応する平均μｍ’をパラメタ推定装置１００から取り出し、推定すべき音源方向として当該平均μｍ’を抽出する。この平均μｍ’は、前記式（４）の左辺のμｎ（ｆ，ｔ）に相当する。従って、方向出力部５２は、前記式（４）の右辺のφｎを求めるために、以下の式（３２）を行う。
φｍ’＝ａｒｃｃｏｓ（μｍ’・ｖ／Ｄ） （３２）
【００７４】
このように、実施例１で説明したパラメタ推定装置１００により推定されたパラメタは周波数依存性の高い空間的エイリアジングの問題を扱うことができるパラメタである。この実施例３の方向推定装置４００は、パラメタ推定装置１００で推定されたパラメタθを用いて音源の方向推定を行うことから、空間的エイリアジングの問題が生じるような条件下であっても、適切な音源の方向推定を行うことができる。
【００７５】
［実験結果］
次に、実施例２で説明した音源分離装置（以下、「本願法」という。）と従来技術で説明した音源分離装置（以下、「従来法」という。）との効果の違いを説明する。まず、図１１を用いて、実験条件について説明する。長手方向４．４５ｍ（＝Ｌｂ、以下、「長手辺」という。）、短手方向３．５５ｍ（＝Ｌａ、以下、「短手辺」という。）、高さ２．５ｍの室内に、第１収音手段２１、第２収音手段２２が部屋の短手方向に一直線上に配置されている。第１収音手段２１と第２収音手段２２との収音手段間隔Ｄは２０ｃｍである。第１収音手段２１と第２収音手段２２（この実験例では両方ともマイクロホン）とを結ぶ線分の中央の点をＣとする。点Ｃから部屋の短手方向の辺までの距離Ｌｂは２．２５ｍとし、点Ｃから部屋の長手方向の辺までの距離Ｌｃは１．７５ｍとする。点Ｃを中心とし、半径０．５５ｍの円をＲとし、図１１記載の円Ｒの円周上に３つの音源（スピーカ）を配置させる。詳細には、点Ｃを通り長手辺と垂直に交わる直線と、円Ｒとが交わる箇所を角度０度とした場合に、時計と反対周りの円周方向の様々な角度（図１１では、４５度、９０度、１３５度）に配置させる。マイクロホンは高さは１．３９ｍとし、スピーカの高さは１．３５ｍとし、サンプリング周波数を１６ｋＨｚとした。これは、音源信号の周波数が８５０Ｈｚ以上で空間的エイリアジング現象が起こる条件である。
【００７６】
図１２に従来法と本願法の音源分離性能を信号対妨害音比（Ｓｉｇｎａｌ ｔｏ ｉｎｔｅｒｆｅｒｅｎｃｅ ｒａｔｉｏ：ＳＩＲ）の改善量を評価した。この実験では、３つのスピーカの配置角度や音声組み合わせを２０通り変更し、それぞれの場合のＳＩＲを求め、平均した値を評価した。図１２からも明らかなように、従来法ではＳＩＲ改善平均量は５．１ｄＢであるが、本願法では、１０．６ｄＢとなり、本願法の方がＳＩＲ改善平均量が大きく、本願法は従来法よりも制度の高い音源分離を可能とすることが理解されよう。
【００７７】
＜ハードウェア構成＞
本発明は上述の実施の形態に限定されるものではない。また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００７８】
また、上述の構成をコンピュータによって実現する場合、パラメタ推定装置１００、音源分離装置３００、方向推定装置４００、が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、処理機能がコンピュータ上で実現される。
【００７９】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記憶しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記憶装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記憶装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００８０】
また、このプログラムの流通は、例えば、そのプログラムを記憶したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００８１】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記憶されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【００８２】
また、本実施例で説明したパラメタ推定装置１００、音源分離装置３００、方向推定装置４００、は、ＣＰＵ（Central Processing Unit）、入力部、出力部、補助記憶装置、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）及びバスを有している（何れも図示せず）。
【００８３】
ＣＰＵは、読み込まれた各種プログラムに従って様々な演算処理を実行する。補助記憶装置は、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、ＲＡＭは、ＳＲＡＭ(Static Random Access Memory)、ＤＲＡＭ (Dynamic Random Access Memory)等である。また、バスは、ＣＰＵ、入力部、出力部、補助記憶装置、ＲＡＭ及びＲＯＭを通信可能に接続している。
【００８４】
＜ハードウェアとソフトウェアとの協働＞
本実施例の単語追加装置は、上述のようなハードウェアに所定のプログラムが読み込まれ、ＣＰＵがそれを実行することによって構築される。以下、このように構築される各装置の機能構成を説明する。
パラメタ推定装置１００、音源分離装置３００、方向推定装置４００、の入力部、出力部は、所定のプログラムが読み込まれたＣＰＵの制御のもと駆動するＬＡＮカード、モデム等の通信装置である。その他の構成部は、所定のプログラムがＣＰＵに読み込まれ、実行されることによって構築される演算部である。記憶部は前記補助記憶装置として機能する。
【技術分野】
【０００１】
本発明は複数の音源信号が混合された観測信号を例えば分離するために用いる確率分布モデルのパラメタを求めるパラメタ推定装置、求められたパラメタを用いた音源分離装置、方向推定装置、それらの方法、プログラムに関する。
【背景技術】
【０００２】
従来技術の音源分離装置５００の機能構成例を図１に示す。この従来技術の詳細は、例えば、非特許文献１に記載されている。Ｎ個の音源１０ｎ（ｎ＝１、．．．、Ｎ）それぞれから同時に発せられる音源信号ｓｎ（ｔ）を２個の収音手段である第１収音手段２１、第２収音手段２２で、ある収録時間内（例えば、５秒間）に観測する。この観測状況を状況Ｘとする。この状況Ｘの下、第１収音手段２１で収音された観測信号をｘ１（ｔ）とし、第２収音手段２２で収音された観測信号をｘ２（ｔ）とし、観測された観測信号をＸ（ｔ）＝［ｘ１（ｔ）、ｘ２（ｔ）］Ｔとする。収音手段とは例えばマイクロホンのことであり、「Ｔ」は行列の転置を表し、ｔを離散時刻とし、ｔ＝０、．．．、Ｔ−１とする。周波数領域変換部３０は、観測信号Ｘ（ｔ）を周波数領域に変換することで、周波数毎の時系列信号である観測信号ベクトルＸ（ｆ，ｔ）＝［ｘ１（ｆ，ｔ），ｘ２（ｆ，ｔ）］Ｔに変換する。Ｘ（ｆ，ｔ）を以下では、周波数観測信号ベクトルという。周波数領域への変換は、例えば短時間フーリエ変換を用いれば良い。ｆは周波数を示し、ｆ∈｛０，ｆｓ／Ｆ，．．．，（Ｆ−１）ｆｓ／Ｆ｝であり、ｆｓはサンプリング周波数を示し、Ｆは周波数帯域の数を示す。
ここで、周波数観測信号ベクトルは、以下の式（１）で表されると仮定する。
【０００３】
【数１】
【０００４】
ただし、ｊ（＝１，２）は収音手段のインデックスであり、ｊ＝１が第１収音手段２１について示し、ｊ＝２が第２収音手段２２について示し、ｈｊｎ（ｆ）は音源１０ｎから収音手段（第１収音手段２１、第２収音手段２２までの周波数応答を示し、ｓｎ（ｆ，ｔ）は音源信号ｓｎ（ｔ）を周波数領域で表現した信号であり、以後では、周波数音源信号ｓｎ（ｆ，ｔ）という。
【０００５】
音源分離を行うために、音源がスパースである、すなわち、周波数音源信号ｓｎ（ｆ，ｔ）はまれにしか大きな値をとらず、ｓｎ（ｆ，ｔ）とｓｎ’（ｆ，ｔ）（ｎ≠ｎ’）は各時間周波数スロット（ｆ，ｔ）で互いに重ならない、ということを仮定する。これは、互いに異なる音声信号などで確認される性質である。これを仮定すると前記式（１）は、
ｘｊ（ｆ，ｔ）≒ｈｊｎ（ｆ）ｓｎ（ｆ，ｔ） （２）
となる。ここで、周波数音源信号ｓｎ（ｆ，ｔ）は時間周波数スロット（ｆ，ｔ）にて支配的な音源信号である。
【０００６】
また、上記式（１）における周波数応答ｈｊｎ（ｆ）が、音源信号ｓｎ（ｔ）が収音手段に到達する間に、部屋の影響を受けない、つまり、無響モデルであると仮定すると、以下の式（３）で表される。
【０００７】
【数２】
【０００８】
ここで、μｎ（ｆ，ｔ）は音源１０ｎが第１収音手段２１と第２収音手段２２に到達する時間差である。また、図２に第１収音手段２１、第２収音手段２２から見た音源の方向角度φｎ（ｆ，ｔ）について示す。図２からも明らかなように、
μｎ（ｆ，ｔ）＝Ｄｃｏｓφｎ（ｆ，ｔ）／ｖ （４）
となる。ただし、Ｄは第１収音手段２１、第２収音手段２２との距離であり、ｖは音速である。
【０００９】
次に、位相差計算部３４は、第１収音手段２１で収音された周波数観測信号ｘ１（ｆ，ｔ）、第２収音手段２２で収音された周波数観測信号ｘ２（ｆ，ｔ）の位相差ｄｆｔを以下の式（５）により計算する。
ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］ （５）
ここで、前記式（２）（３）より音源１０ｎが支配的な時間周波数スロット（ｆ，ｔ）では以下の式（６）、（７）が成り立つ。
【００１０】
【数３】
【００１１】
次に、分類部７０は、位相差ｄｆｔを用いてクラスタリング処理を行う。まず、分類部７０は、以下の式（８）により位相差ｄｆｔを周波数ｆで正規化する。
【００１２】
【数４】
【００１３】
分類部７０は、式（８）についてクラスタリングを行う。スパース性の仮定より定められる前記式（２）より、音源１０ｎだけが支配的な時間周波数（ｆ，ｔ）では時間差μｎ（ｆ，ｔ）が求まっており、音源１０ｎ’だけが支配的な時間周波数（ｆ，ｔ）では時間差μｎ’（ｆ，ｔ）が求まっているので、音源１０ｎの方向角度φｎ（ｆ，ｔ）と音源１０ｎ’の方向角度φｎ７（ｆ，ｔ）とが異なっていれば、両者を分類（クラスタリング）できる。そして、それぞれのクラスタが各音源成分に対応する。
【００１４】
音源分離部７２では、音源１０ｎに対応するクラスタを形成している時間周波数（ｆ，ｔ）では１をとり、それ以外の時間周波数では０をとるマスクＬｎ（ｆ，ｔ）を生成する。このマスクＬｎ（ｆ，ｔ）を全ての音源１０ｎについて生成する。
【００１５】
そして、音源分離部７２は、マスクＬｎ（ｆ，ｔ）を観測信号の１つ（ここではｘ１（ｆ，ｔ））に乗算することで、つまり以下の式（１０）を演算することで、分離信号ｙｎ（ｆ，ｔ）を得る。
ｙｎ（ｆ，ｔ）＝ｘ１（ｆ，ｔ）Ｌｎ（ｆ，ｔ） （１０）
そして、時間領域変換部７４は、分離信号ｙｎ（ｆ，ｔ）を時間領域に変換することで、時間領域分離信号ｙｎ（ｔ）を得る。
【先行技術文献】
【非特許文献】
【００１６】
【非特許文献１】S．Araki、H．Sawada、R．Mukai、and S．Makino、”Underdetermined blind sparse source separation for arbitrarily arranged multiple sensors、”Signal Processing、vol．87、pp．1833−1847、Feb．2007．
【発明の概要】
【発明が解決しようとする課題】
【００１７】
従来技術では、第１収音手段２１および第２収音手段２２との間の収音手段間隔Ｄが狭く、音源信号の周波数ｆが低い場合には、非特許文献１に記載されているように良好に動作する。しかし、収音手段間隔Ｄが広い場合や、周波数ｆが高い場合にはクラスタリング処理がうまく動作しないため、従来技術では適切に音源分離を行うことは難しかった。図３Ａ〜Ｆに音源１０ｎが２つの場合、つまり、ｎ＝１，２となる場合に、収音手段間隔Ｄおよびサンプリング周波数ｆｓを変えた場合のクラスタリング処理の模様を示す。図３Ａ〜Ｆにおいて、実線を音源１０１の音源信号を示し、破線を音源１０２の音源信号を示し、横軸はそれぞれの音源信号ｓ１（ｆ，ｔ）、音源信号ｓ２（ｆ，ｔ）の周波数ｆを示す。また、図３Ａ、Ｂの縦軸はそれぞれ、サンプリング周波数ｆｓを８ｋＨｚ、収音手段間隔Ｄを４ｃｍとした（Ｄが狭い）場合の、位相差ｄｆｔ、正規化後の位相差ｄｆｔ／２πｆを示し、図３Ｃ、Ｄの縦軸はそれぞれ、サンプリング周波数ｆｓを８ｋＨｚ、収音手段間隔Ｄを１０ｃｍとした（Ｄが広い）場合の、位相差ｄｆｔ、正規化後の位相差ｄｆｔ／２πｆを示し、図３Ｅ、Ｆの縦軸はそれぞれ、サンプリング周波数ｆｓを１６ｋＨｚ、収音手段間隔Ｄを４ｃｍとした（Ｄが狭い）場合の、位相差ｄｆｔ、正規化後の位相差ｄｆｔ／２πｆを示す。この図３を用いて、従来技術の問題点を説明する。
【００１８】
図３Ｂに示すように、収音手段間隔Ｄが狭く、周波数ｆが低い場合には、全ての周波数ｆで、位相差ｄｆｔ／２πｆが同じ値をとっており、これをクラスタリングすると２個のクラスタが形成され、適切に音源分離を行うことができる。
【００１９】
一方、図３Ｄに示すように収音手段間隔Ｄが狭い場合には、ｆ＞３０００Ｈｚの範囲では、位相差ｄｆｔ／２πｆが一定でなくなる。このようなｄｆｔ／２πｆをクラスタリングしても、音源ごとのクラスタは形成されず、適切に音源分離を行うことができない。この理由は、前記式（７）の右辺が収音手段間隔Ｄが大きいために±πの範囲を超えてしまうにも関わらず、式（６）のａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］の演算が−πからπの範囲の値を算出するためである。これを具体的に説明すると、ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］のａｒｇ演算が−πからπの間の値しか返せないため
２πｆμｎ（ｆ，ｔ）＝ｄｆｔ＋２πｋ （１１）
＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］＋２πｋ
（１２）
の関係を満たす−π＜ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］＜πとして返されてしまうためである。これを「空間的エイリアジングの問題」や、「２πｋの不定性がある」という。ここで、ｋはμｎ（ｆ，ｔ）の値が既知であれば、一意に決まる整数であり、不定性係数ｋという。ただし、一般的には、μｎ（ｆ，ｔ）の値は未知であるため、不定性係数ｋは推定すべき整数である。
【００２０】
この空間的エイリアジングの問題は、図３Ｅ、図３Ｆ（特に、図３Ｆ）に示すように、収音手段間隔Ｄが狭くても（Ｄ＝４ｃｍ）、周波数（図３Ｆの例では、周波数ｆが高くなる（周波数ｆが６０００Ｈｚ以上）場合には、生じる問題である。
【００２１】
このように従来技術では、空間的エイリアジングの問題が生じるような、収音手段間隔Ｄが狭く、周波数ｆが高くなる場合には、適切なクラスタリングができないため、良好な音源分離ができなかった。
【課題を解決するための手段】
【００２２】
上記の課題を解決するために、この出願のパラメタ推定装置は、周波数領域変換部と、位相差計算部と、推定部と、を有する。周波数領域変換部は、複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成する。位相差計算部は、周波数観測信号の収音手段間の位相差を計算する。推定部は、位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する。
【発明の効果】
【００２３】
本発明のパラメタ推定装置により推定されたパラメタθを用いて、音源分離を行うことで、空間的エイリアジングの問題が生じる場合においても、良好に音源分離を行うことができる。
【図面の簡単な説明】
【００２４】
【図１】従来の音源分離装置の機能構成例を示した図。
【図２】収音手段間距離Ｄなどを説明した図。
【図３】従来の音源分離装置のクラスタリングの様子を示した図。
【図４】本実施例のパラメタ推定装置の機能構成例を示した図。
【図５】本実施例のパラメタ推定装置の処理フローを示した図。
【図６】位相差の分布に正規分布モデルを当てはめることを示した図。
【図７】本実施例の推定部の機能構成例を示した図。
【図８】本実施例の推定部の処理フローを示した図。
【図９】本実施例の音源分離装置の機能構成例を示した図。
【図１０】本実施例の方向推定装置の機能構成例を示した図。
【図１１】実験を行った部屋を示した図
【図１２】実験結果を示した図。
【発明を実施するための形態】
【００２５】
以下に、発明を実施するための最良の形態を示す。なお、同じ機能を持つ構成部や同じ処理を行う過程には同じ番号を付し、重複説明を省略する。上述の課題を解決する音源分離装置については実施例２で説明し、実施例１では当該音源分離装置に用いられるパラメタ推定装置について説明する。実施例３では、当該パラメタ推定装置を用いた方向推定装置について説明する。
【実施例１】
【００２６】
図４に実施例１のパラメタ推定装置１００の機能構成例を示し、図５に処理フローを示す。このパラメタ推定装置１００を用いた音源分離装置３００の機能構成例を図９および実施例２に示す。また、前記状況Ｘの下、第１収音手段２１で収音された観測信号ｘ１（ｔ）および第２収音手段２２で収音された観測信号ｘ２（ｔ）は、周波数領域変換部３０に入力される。周波数領域変換部３０は、観測信号Ｘ（ｔ）＝［ｘ１（ｔ），ｘ２（ｔ）］Ｔを周波数領域に変換することで周波数観測信号（ベクトル）Ｘ（ｆ，ｔ）＝［ｘ１（ｆ，ｔ），ｘ２（ｆ，ｔ）］Ｔを生成する（ステップＳ１０２）。
【００２７】
周波数観測信号Ｘ（ｆ，ｔ）は入力作成部３３に入力される。入力作成部３３は、パワー推定部３２および位相差計算部３４とで構成される。パワー推定部３２は以下の式（１３）により時間周波数スロット（ｆ，ｔ）ごとにパワーａｆｔを計算する。
ａｆｔ＝│ｘｊ（ｆ，ｔ）│ ただしｊ＝１，２ （１３）
位相差計算部３４は時間周波数スロット（ｆ，ｔ）ごとに、第１収音手段２１と第２収音手段２２との間の位相差ｄｆｔを計算する（ステップＳ１０４）。具体的には、前記式（５）により計算する。念のため式（５）を以下に示す。
ｄｆｔ＝ａｒｇ［ｘ１（ｆ，ｔ）／ｘ２（ｆ，ｔ）］ （５）
【００２８】
次に、推定部２００は、位相差ｄｆｔの分布（ヒストグラム）に当てはまり、かつ、周波数依存性のある確率分布モデルのパラメタθを推定する（ステップＳ１０６）。換言すれば、各音源に対応するクラスタを確率分布モデル（例えば、正規分布）で当てはめる（近似するまたはモデル化する）場合のパラメタ推定を行う。「周波数依存性のある確率分布モデル」については後ほど詳細に述べる。また、推定した確率分布モデルのパラメタθは、例えば、以下で説明する音源分離処理（実施例２）や方向推定処理（実施例３）で用いる。以下の説明では、確率分布モデルとして正規分布を用いた場合の説明を行うが、用いる確率分布モデルは、正規分布モデルに限らない。
【００２９】
ここでは、まず、各音源に対応するクラスタを正規分布モデルで当てはめる例の説明を行う。まず、音源数が１（クラスタが１個）の場合を考える。このとき、位相差ｄｆｔの分布を書くと、図６Ａに示すようになる。図６Ａからも理解できるように、１つの山（＝クラスタ）ができる。このクラスタに以下の式（１４）で示される正規分布モデルを当てはめる。
【００３０】
【数５】
【００３１】
ここで、μは正規分布のパラメタθのうち「平均」を示し、σ２は正規分布のパラメタθのうち「分散」を示す。特に、μは前記式（３）で説明した、音源からの音が第１収音手段２１、第２収音手段２２に到達する時間の時間差であり、後述する正規分布モデルの当てはめを行うことで推定する値である。また上述の通り、ｋは不定性係数であり、後述する正規分布モデルの当てはめを行うことで推定する必要がある。また。Ｋｆは不定性係数ｋの最大値であり、周波数ｆによって異なる値をとりうる。具体的には、例えば、以下の式（１５）で表される。
【００３２】
【数６】
【００３３】
また、全ての周波数ｆで同じ値（ただし、十分大きな値）を用いても良い。
次に、音源数が２以上（すなわちクラスタが２つ以上）の場合を考える。図６Ｂに音源数が２の場合の位相差ｄｆｔの分布を示す。図６Ｂからも理解されるように、２つの山（＝クラスタ）ができることがわかる。このクラスタにそれぞれ１つの正規分布モデルを当てはめることを考える。つまり、分布全体を以下の式（１６）で示される混合正規分布モデルＧ（ＧＭＭ）で当てはめることを考える。
【００３４】
【数７】
【００３５】
ここで、ｍは正規分布のインデックスを示しつまり、（ｍ＝１，．．．，Ｍ）であり、図６Ｂの例では、Ｍ＝２となり、μｍはｍ番目の正規分布の平均を示し、σｍ２はｍ番目の正規分布の分散を示し、特に、μｍは音源からの音が第１収音手段２１、第２収音手段２２に到達する時間の時間差の後述する正規分布モデルの当てはめを行うことで推定する値である。また、αｍは、ｍ番目の正規分布の混合重みであり、Σｍ＝１Ｍαｍ＝１であり、０≦αｍ≦１である。また、混合正規分布を構成するＭ個の正規分布をΨｍ（ｍ＝１，．．．，Ｍ）とし、正規分布の数Ｍ（以下、「混合数Ｍ」という。）は、音源数Ｎが既知の場合には、Ｍ＝Ｎとすることができる。また、音源数が未知の場合には、Ｍを十分大きな数（例えば、Ｍ＝１０）を用いれば良い。本実施例では、複数音源の分離などの処理を行うために、位相差ｄｆｔに当てはまる前記式（１６）に示す混合正規分布を用いてパラメタθを推定する。以下、推定部２００のパラメタθの推定処理について詳細に説明する。
【００３６】
図７に推定部２００の機能構成例を示し、図８に推定部２００の処理フローを示す。パラメタθを混合正規分布の平均μ、分散σ２、混合重みαをまとめたものを示し、つまり、θ＝（μｍ，σｍ，αｍ）＝（μ１，σ１，α１，．．．，μｍ，σｍ，αｍ，．．．，μＭ，σＭ，αＭ）となる。また、ｒを更新回数とし、θに更新回数の概念を付与したもの、つまり、ｒ回更新したθをθｒとすると、θｒ＝（μ１ｒ，σ１ｒ，α１ｒ，．．．，μｍｒ，σｍｒ，αｍｒ，．．．，μＭｒ，σＭｒ，αＭｒ）となる。また、記憶部１６には予め用いる正規分布モデルのモデル数Ｍと混合正規分布モデルの各パラメタの初期値θ０が記憶されている。事前分布情報保持部１１０には、ハイパーパラメタω（後述する）、重みパラメタｃが保持されている。
【００３７】
推定部２００には、パワー推定部３２よりのパワーａｆｔが重み係数ａｆｔとして入力され、位相差計算部３４よりの位相差ｄｆｔも入力される。または全てのａｆｔについてａｆｔ＝１としても良く、この場合は、パワー推定部３２はなくてもよい。また、重み係数ａｆｔを各時間周波数（ｆ，ｔ）における観測信号のパワーや信号の瞬時的ＳＮ比などとすることも出来る。従って、パワー推定部３２を設けない場合であっても、重み係数ａｆｔを入力部３５から入力することができる。
【００３８】
まず、初期設定として、ｒ＝０（つまり更新回数が０）、ｒ＝０のときの混合正規分布のパラメタθ０の値、用いる正規分布のモデル数Ｍ、不定性係数ｋの範囲であるＫｆ、更新回数閾値Ｒまたは差閾値Δ（後述する）を設定する（ステップＳ２）。更新回数閾値Ｒまたは差閾値Δは、後述する収束判定処理の際に用いられる。これらの初期設定は、入力部３５からの入力により行われる。
【００３９】
事後確率計算部１２は、位相差ｄｆｔと、現在の確率分布モデル（混合正規分布モデル）のパラメタθｒ（＝（μｍｒ，σｍｒ，αｍｒ ｍ＝１，．．．，Ｍ））から、Ｍ個の正規分布Ψｍごとに事後確率ｐ（ｍ，ｋ│ｄｆｔ，θｒ）を計算する（ステップＳ６）。またパラメタ保持部１８には、現在の混合正規分布のθｒが保持されている。事後確率計算部１２は具体的には例えば、以下の式（１７）（１８）により計算する。
【００４０】
【数８】
【００４１】
次に、更新部１４は、位相差ｄｆｔと事後確率ｐ（ｍ，ｋ│ｄｎ，θｒ）を用いて、現在の混合正規分布の各パラメタθｒを更新する（ステップＳ８）。以下、更新処理について詳細に説明する。更新部１４は更新処理の際に、ハイパーパラメタω、重みパラメタｃを事前分布情報保持部１１０から取り出す。この実施例では、音源数Ｎが未知、つまり、正規分布の数Ｍが未知の場合であっても、適用可能にするために、正規分布のパラメタθの混合重みαｍに適切な事前分布を与え、例えばＥＭアルゴリズムにてパラメタθの更新処理を行う。この実施例１では、混合重みαｍの事前分布として、ディリクレ分布を考える。ディリクレ分布の詳細は、参考文献１である「Ｃ．Ｍ．ビショップ著（元田、栗田他訳） 「パターン認識と機械学習（上）」、シュプリンガー・ジャパン２００７年 ｐ．７４−ｐ．７７」等に記載されている。ディリクレ分布は例えば以下の式（１９）で表される。
【００４２】
【数９】
【００４３】
ここで、αは混合重み行列であり、α＝｛α１，．．．，αｍ，．．．，αＭ｝で表され、ΣｍＭαｍ＝１、０≦αｍ≦１という条件を満たす。これは混合正規分布のパラメタである混合重みの条件と同じであることに注意されたい。またβ（ω）は正規化項（ベータ分布）であり、ここで、ハイパーパラメタωを１より小さい正の値（例えば、０．９）に設定すると、αｍのごく少数のみが十分に大きな値を持ち、残りは０に近い値をとるようになる。求められたαｍを前記式（１６）で用いられている混合重みαｍの事前分布として用いることで、混合正規分布モデルＧのうちの少数の正規分布のみに十分大きな混合重みがかかり、その他の正規分布モデルの混合重みは０に近くなる。結果として、なるべく少数の正規分布による当てはめが可能である。従って、１つのクラスタに複数の正規分布が当てはまるような現象を防ぐことが出来、音源数未知数の場合でも、それぞれのクラスタに１つずつ正規分布を当てはめることができる。
【００４４】
次に、この事前分布を含みながら、パラメタ更新を行うためのＥＭアルゴリズムを導出する。ここで、正規分布のインデックスｍと不定性係数ｋは位相差ｄｆｔから推定すべき変数であるため、ＥＭアルゴリズムにおける隠れ変数として扱う。これにより不定性係数ｋは隠に自動推定されるため、位相差ｄｆｔにおける２πｋの不定性を自動的に扱うことが可能になる。まず、最尤推定のためのコスト関数Ｌ（θ）は次のように与えられる。
【００４５】
【数１０】
【００４６】
また、重みパラメタｃは、式（２２）の第１項と第２項の重みをコントロールするパラメタであり、上述のように、事前分布情報保持部１１０に保持される。
【００４７】
【数１１】
【００４８】
となる。ここで、式（２４）のＥ［Ｈ］は式Ｈの期待値を示し、式（２５）中のｐ（ｍ，ｋ│ｄｆｔ，θｒ）は式（１８）で表される事後確率分布である。
【００４９】
【数１２】
【００５０】
図７中の更新部１４中の平均更新手段１４２が式（２６）より現在の平均μｍｒを更新することで更新後の平均μｍｒ＋１を出力する。分散更新手段１４４が式（２７）より分散（σｍ２）ｒを更新することで更新後の分散（σｍ２）ｒ＋１を出力する。混合重み更新手段１４６が式（２８）により混合重みαｍを更新することで更新後の混合重みαｍｒ＋１を出力する。パラメタ算出手段が、更新後の平均μｍｒ＋１、分散（σｍ２）ｒ＋１、混合重みαｍｒ＋１についての更新後のパラメタθｒ＋１を算出する（ステップＳ８）。
【００５１】
各パラメタの更新処理が数回行われ（ステップＳ４）、更新部１４内の収束判定手段１５０は、更新されたθｒ＋１に対して、予め定められた規則により、各パラメタ値が収束しているか否かの収束判定を行う（ステップＳ１０）。各パラメタ値が収束していると判断した場合には、更新されたパラメタθｒ＋１を出力する。また、各パラメタ値が収束していないと判断した場合には、更新されたパラメタθｒ＋１を現在の確率分布モデルの平均、分散、混合重みとしてパラメタ保持部１８に保持させる。を繰り返す。そして、収束判定手段１５０が、各パラメタ値が収束していると判断するまで、ステップＳ４〜ステップＳ１０の処理（平均更新手段１４２、分散更新手段１４４、混合重み更新手段１４６の処理）を繰り返す。
【００５２】
ここで収束判定に用いる予め定められた規則の例を説明する。更新回数閾値Ｒを用いる例を説明すると、更新部１４内のカウント手段（図示せず）は更新回数ｒをカウントし、更新回数ｒが更新回数閾値Ｒ（例えば３０）を超えた場合には、十分更新しており、収束していると判断して、パラメタ算出手段１４８は、更新後のパラメタθを出力する。また、差閾値Δを用いる例を説明すると、以下の式（２９）の式を満たす場合には、収束していると判断して、パラメタ算出手段１４８は、更新後のパラメタθを出力する。
│Ｑ（θ│θｒ＋１）−Ｑ（θ│θｒ）│＜Δ （２９）
このようにして、推定部２００は、各音源に対応するクラスタを正規分布モデルで当てはめたときのパラメタθ＝（μｍ，σｍ，αｍ）を出力する。
【００５３】
この実施例１では混合重みαｍのみに事前分布を導入したが、各ガウス分布の平均μｍと分散σｍに対しても事前分布を導入することで、より精度の高い混合正規分布の当てはめを実現できる。また、各ガウス分布の各パラメタである平均μｍ、分散σｍ、混合重みαｍ、に事前分布を導入した場合の当てはめ処理には、ＥＭアルゴリズムの他、不定性係数ｋを隠に自動推定するアルゴリズムであれば、何でも良い。これらの拡張は当業者であれば、上記参考文献１などを参照すれば、容易に実現できるため、ここでは省略する。
【００５４】
また、音源数Ｎが既知であり、混合数Ｍ＝Ｎとできれば、前記式（１９）で示されるディリクレ分布を用いる必要はない。すなわち、この場合には、前記式（２８）でハイパーパラメタω＝１とすればよい。例えば、ユーザが入力部３５から音源数Ｎが既知であるか未知であるかを示す情報である音源数情報を入力し、音源数情報が既知である旨の情報であれば、ディリクレ分布を用いず、音源数情報が未知である旨の情報であれば、ディリクレ分布を用いれば良い。
【００５５】
図３で説明したように、空間的エイリアジングの影響は周波数毎に異なる、つまり、周波数依存性があるといえる。よって、確率分布モデルでモデル化する場合には、この周波数依存性を考えることが必要である。
【００５６】
これを式（１４）（１６）に示す混合正規分布モデルについて検討する。特にｅｘｐ項の分子の式「−（ｄｆｔ＋２πｋ−２πｆμｍ）２」を検討すると、まず、２πｋについては、上述したように、パラメタ推定は例えばＥＭアルゴリズムを用いて行うが、この際、不定性係数ｋはＥＭアルゴリズムにおける隠れ変数として扱うことができる。よって不定性係数ｋはＥＭアルゴリズムでデータから（確率的に）自動推定されるため、位相差ｄｆｔの２πｋの不定性を自動的に扱うことができる。また、ｄｆｔについては式（１１）、（１２）の通り、周波数領域の観測信号ｘ１（ｆ，ｔ）、ｘ２（ｆ，ｔ）の位相差であるため、周波数依存性のある値であるといえる。２πｆμｍについても周波数ｆを含んでいるので、周波数依存性のある値であるといえる。つまり、式（１４）（１６）に示す混合正規分布モデルは、周波数依存性のあるモデルであるといえる。また、推定部２００で用いる確率分布モデルは、周波数依存性があれば、式（１４）（１６）に限らない。
【００５７】
また、式（１１）において、実測値（観測された値）は式（１１）の右辺に示すｄｆｔ＋２πｋであり、当てはめ後（モデル化後）の値は、式（１１）左辺の２πｆμｎ（ｆ，ｔ）である。式（１６）のｅｘｐ項の分子「（ｄｆｔ＋２πｋ−２πｆμｎ）２」は、実数値とモデル化後の値との二乗誤差（モデル化誤差）を意味するとも捉えることができる。観測された位相差ｄｆｔのヒストグラムに当てはまる正規分布は、このモデル化誤差を最も小さくする正規分布ということになり、前記式（２６）〜（２８）を用いて、当該正規分布を求めている。
【００５８】
従来技術では、前記式（８）に示すように、ｄｆｔ／２πｆのように、周波数正規化した量をクラスタリングしていた。そのため、図３Ｄ、Ｆに示すように、正しくクラスタリングできなかった。しかし、本実施例では、周波数依存性のあるｄｆｔの分布をそのまま前記式（１６）で示す混合正規分布に当てはめる。具体的には、本実施例のパラメタ推定装置１００は、前記式（１１）（１２）で与えられる空間的エイジアリングの問題（２πｋの不定性）を陽に定式化し、位相差ｄｆｔのヒストグラムに当てはまり、周波数依存性のある混合正規分布（式（１４）や式（１６））のパラメタを求める。この求められたパラメタを用いて例えば音源分離を行うと、周波数依存性の高い空間的エイリアジングの問題を扱うことが可能となり、適切な音源分離などを行うことができる。
【００５９】
また、このパラメタ推定装置１００で求められたパラメタθは、実施例２で説明する音源分離処理や、実施例３で説明する音源方向推定処理のほか、様々な観測信号処理に用いられる。
【実施例２】
【００６０】
この実施例２では、実施例１で説明したパラメタ推定装置１００により推定されたパラメタθを用いて、音源分離を行う音源分離装置について説明する。図９に実施例２の音源分離装置３００の機能構成例を示す。またパラメタ推定装置で推定されたパラメタθを決定後パラメタという。
【００６１】
パラメタ推定装置１００よりの決定後パラメタθは有効音源推定部４０に入力される。有効音源推定部４０は、音源に該当する確率分布モデルを示す音源該当情報を求める。有効音源推定部４０による音源該当情報の生成手法は以下の３つの手法により求められる。ここで、音源該当情報とは、例えば、音源に該当する確率分布モデルのインデックスｍをいう。
【００６２】
まず、第１の手法として、音源数Ｎが既知であり、パラメタ推定装置１００で用いられる混合正規分布の混合数Ｍ＝Ｎとしている場合には、混合正規分布を構成する全ての正規分布が音源に該当するので、全ての正規分布のインデックスｍ（＝１，．．．，Ｍ）を出力する。また音源数Ｎが未知の場合には、下記の第２手法、第３の手法により求められる。
【００６３】
第２の手法として、有効音源推定部４０は、決定後パラメタθの混合重みαｍ（ｍ＝１，．．．，Ｍ）のうち、混合重みが予め定められた第１閾値ε１（例えば１０−６）よりも大きな値である正規分布を音源に該当する正規分布と判断して、当該正規分布（以下、「音源該当正規分布」という。）のインデックスｍ’を出力する。何故なら、パラメタ推定装置１００の演算が十分収束している場合には、決定後パラメタθの中の混合重みαｍのうち十分大きな値を持つ個数は位相差ｄｆｔのヒストグラム中の分布の山の個数と等しくなるからである。また、音源該当正規分布の数をＭ’（つまり、ｍ’＝１，．．．，Ｍ’）とする。
【００６４】
また、第３の手法として、第２の手法においてパラメタ推定装置１００の演算が十分に収束していない場合は、有効音源推定部４０は、次のような推定処理を行うことが好ましい。有効音源推定部４０は、混合重みαｍが第１閾値ε１よりも大きく、かつ分散σｍ２が予め定められた第２閾値ε２（例えばπ／５）よりも小さい正規分布を音源に該当する正規分布と判断して、音源該当正規分布のインデックスｍ’を音源該当情報として出力する。
【００６５】
また、音源数Ｎが既知であるか未知であるかについての情報である音源数情報（つまり、第１の手法を用いるか、または第２、第３の手法を用いるか）は、ユーザに入力部４７から入力させればよい。
【００６６】
次に、マスク作成部４２は、音源該当情報（ここでは、音源該当正規分布のインデックスｍ’）が示す確率分布モデル（正規分布モデル）を周辺化することでマスクΩｍ’（ｆ，ｔ）を作成する。マスクΩｍ’（ｆ，ｔ）は、各音源該当正規分布Ψｍ’ごとに、かつ、各時間周波数スロット（ｆ，ｔ）ごとに求められる。具体的には、Ｍ’個の音源該当正規分布に関する事後確率ｐ（ｍ’，ｋ│ｄｆｔ，θｒ）を周辺化することで、周辺化事後確率ｐ（ｍ’│ｄｎ，θｒ）（＝マスクΩｍ’（ｆ，ｔ））を求める。
【００６７】
【数１３】
【００６８】
式（３０）中のｐ（ｍ’，ｋ│ｄｆｔ，θ）については、マスク作成部４２がパラメタ推定装置１００中の事後確率計算部１２から前記式（１８）の結果を抽出すればよい。
【００６９】
そして、分離部４４は、周波数観測信号にマスクΩｍ’（ｆ，ｔ）を乗算することで、分離信号ｙｍ’（ｆ，ｔ）を求める。具体的には以下の式（３１）により分離信号ｙｍ’（ｆ，ｔ）を求める。例えば、分離部４４は、マスクΩｍ’（ｆ，ｔ）を観測信号の１つ（ここでは、周波数観測信号ｘ１（ｆ，ｔ））に乗算し、分離信号ｙｍ’（ｆ，ｔ）を得る。つまり、以下の式（３１）により求められる。
ｙｍ’（ｆ，ｔ）＝ｘ１（ｆ，ｔ）Ωｍ’（ｆ，ｔ） （３１）
分離部４４よりの分離信号ｙｍ’（ｆ，ｔ）は、時間領域変換部４６に入力される。そして、時間領域変換部４６は、分離信号ｙｍ’（ｆ，ｔ）を時間領域に変換して時間領域分離信号ｙｍ’（ｔ）を求め、出力する。
【００７０】
ここで、従来の音源分離装置５００（図１参照）と、実施例２の音源分離装置３００（図９参照）の対応関係を以下に示す。
音源分離装置５００の音源分離部７２
→音源分離装置３００のマスク作成部４２と分離部４４を統合したもの
音源分離装置５００の周波数領域変換部３０と位相差計算部３４と分類部７０を統合
したもの→音源分離装置３００のパラメタ推定装置１００
また、従来の音源分離装置５００では、音源数が既知である場合が多いため、音源分離装置５００は、音源分離装置３００の有効音源推定部４０に対応するものを有していなかった。
【００７１】
このように、実施例１で説明したパラメタ推定装置１００により推定されたパラメタは周波数依存性の高い空間的エイリアジングの問題を扱うことができるパラメタである。この実施例２の音源分離装置３００は、パラメタ推定装置１００で推定されたパラメタθを用いて音源分離を行うことから、空間的エイリアジングの問題が生じるような条件下であっても、適切な音源分離を行うことができる。
【実施例３】
【００７２】
実施例１で説明したパラメタ推定装置１００により推定されたパラメタθを用いて、音源の方向を推定することもできる。この実施例３では、実施例１で説明したパラメタ推定装置１００を用いた、音源の方向を推定する方向推定装置４００を説明する。図１０に方向推定装置４００の機能構成例を示す。方向推定部６０は、有効音源推定部５０と方向出力部５２とで構成されている。
【００７３】
パラメタ推定装置１００よりの決定後パラメタθは、有効音源推定部５０および方向出力部５２に入力される。有効音源推定部５０は、実施例２で説明したように、音源に該当する確率分布モデルを示す音源該当情報（例えば、音源該当正規分布のインデックスｍ’）を求める。有効音源推定部５０の処理が終了すると、方向出力部５２は、音源該当方向情報分布モデルのインデックスｍ’｛ｍ’＝１，．．．，Ｍ’｝に対応する平均μｍ’をパラメタ推定装置１００から取り出し、推定すべき音源方向として当該平均μｍ’を抽出する。この平均μｍ’は、前記式（４）の左辺のμｎ（ｆ，ｔ）に相当する。従って、方向出力部５２は、前記式（４）の右辺のφｎを求めるために、以下の式（３２）を行う。
φｍ’＝ａｒｃｃｏｓ（μｍ’・ｖ／Ｄ） （３２）
【００７４】
このように、実施例１で説明したパラメタ推定装置１００により推定されたパラメタは周波数依存性の高い空間的エイリアジングの問題を扱うことができるパラメタである。この実施例３の方向推定装置４００は、パラメタ推定装置１００で推定されたパラメタθを用いて音源の方向推定を行うことから、空間的エイリアジングの問題が生じるような条件下であっても、適切な音源の方向推定を行うことができる。
【００７５】
［実験結果］
次に、実施例２で説明した音源分離装置（以下、「本願法」という。）と従来技術で説明した音源分離装置（以下、「従来法」という。）との効果の違いを説明する。まず、図１１を用いて、実験条件について説明する。長手方向４．４５ｍ（＝Ｌｂ、以下、「長手辺」という。）、短手方向３．５５ｍ（＝Ｌａ、以下、「短手辺」という。）、高さ２．５ｍの室内に、第１収音手段２１、第２収音手段２２が部屋の短手方向に一直線上に配置されている。第１収音手段２１と第２収音手段２２との収音手段間隔Ｄは２０ｃｍである。第１収音手段２１と第２収音手段２２（この実験例では両方ともマイクロホン）とを結ぶ線分の中央の点をＣとする。点Ｃから部屋の短手方向の辺までの距離Ｌｂは２．２５ｍとし、点Ｃから部屋の長手方向の辺までの距離Ｌｃは１．７５ｍとする。点Ｃを中心とし、半径０．５５ｍの円をＲとし、図１１記載の円Ｒの円周上に３つの音源（スピーカ）を配置させる。詳細には、点Ｃを通り長手辺と垂直に交わる直線と、円Ｒとが交わる箇所を角度０度とした場合に、時計と反対周りの円周方向の様々な角度（図１１では、４５度、９０度、１３５度）に配置させる。マイクロホンは高さは１．３９ｍとし、スピーカの高さは１．３５ｍとし、サンプリング周波数を１６ｋＨｚとした。これは、音源信号の周波数が８５０Ｈｚ以上で空間的エイリアジング現象が起こる条件である。
【００７６】
図１２に従来法と本願法の音源分離性能を信号対妨害音比（Ｓｉｇｎａｌ ｔｏ ｉｎｔｅｒｆｅｒｅｎｃｅ ｒａｔｉｏ：ＳＩＲ）の改善量を評価した。この実験では、３つのスピーカの配置角度や音声組み合わせを２０通り変更し、それぞれの場合のＳＩＲを求め、平均した値を評価した。図１２からも明らかなように、従来法ではＳＩＲ改善平均量は５．１ｄＢであるが、本願法では、１０．６ｄＢとなり、本願法の方がＳＩＲ改善平均量が大きく、本願法は従来法よりも制度の高い音源分離を可能とすることが理解されよう。
【００７７】
＜ハードウェア構成＞
本発明は上述の実施の形態に限定されるものではない。また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００７８】
また、上述の構成をコンピュータによって実現する場合、パラメタ推定装置１００、音源分離装置３００、方向推定装置４００、が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、処理機能がコンピュータ上で実現される。
【００７９】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記憶しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記憶装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記憶装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００８０】
また、このプログラムの流通は、例えば、そのプログラムを記憶したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００８１】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記憶されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【００８２】
また、本実施例で説明したパラメタ推定装置１００、音源分離装置３００、方向推定装置４００、は、ＣＰＵ（Central Processing Unit）、入力部、出力部、補助記憶装置、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）及びバスを有している（何れも図示せず）。
【００８３】
ＣＰＵは、読み込まれた各種プログラムに従って様々な演算処理を実行する。補助記憶装置は、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、ＲＡＭは、ＳＲＡＭ(Static Random Access Memory)、ＤＲＡＭ (Dynamic Random Access Memory)等である。また、バスは、ＣＰＵ、入力部、出力部、補助記憶装置、ＲＡＭ及びＲＯＭを通信可能に接続している。
【００８４】
＜ハードウェアとソフトウェアとの協働＞
本実施例の単語追加装置は、上述のようなハードウェアに所定のプログラムが読み込まれ、ＣＰＵがそれを実行することによって構築される。以下、このように構築される各装置の機能構成を説明する。
パラメタ推定装置１００、音源分離装置３００、方向推定装置４００、の入力部、出力部は、所定のプログラムが読み込まれたＣＰＵの制御のもと駆動するＬＡＮカード、モデム等の通信装置である。その他の構成部は、所定のプログラムがＣＰＵに読み込まれ、実行されることによって構築される演算部である。記憶部は前記補助記憶装置として機能する。
【特許請求の範囲】
【請求項１】
複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成する周波数領域変換部と、
前記周波数観測信号の前記収音手段間の位相差を計算する位相差計算部と、
前記位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する推定部と、を有するパラメタ推定装置。
【請求項２】
請求項１記載のパラメタ推定装置であって、
前記推定部は、
現在の前記確率分布モデルの各パラメタを保持するパラメタ保持部と、
前記推定部は、前記位相差と、前記現在の確率分布モデルの各パラメタと、を用いて確率分布モデルごとに事後確率を計算する事後確率計算部と、
確率分布モデルの各パラメタ値を更新する更新部と、を有するパラメタ推定装置。
【請求項３】
請求項１または２記載のパラメタ推定装置と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定部と、
前記音源該当情報が示す確率分布モデルごとのマスクを作成するマスク作成部と、
前記周波数観測信号に前記マスクを乗算することで、分離信号を求める分離部と、
前記分離信号を時間領域に変換する時間領域変換部と、を有することを特徴とする音源分離装置。
【請求項４】
請求項１または２記載のパラメタ推定装置と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定部と、
前記音源該当情報が示す確率分布モデルの平均を出力する方向出力部と、を有することを特徴とする方向推定装置。
【請求項５】
複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成する周波数領域変換過程と、
前記周波数観測信号の前記収音手段間の位相差を計算する位相差計算過程と、
前記位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する推定過程と、を有するパラメタ推定方法。
【請求項６】
請求項５記載のパラメタ推定方法であって、
前記推定過程は、
現在の前記確率分布モデルの各パラメタを保持するパラメタ保持過程と、
前記推定過程は、前記位相差と、前記現在の確率分布モデルの各パラメタと、を用いて確率分布モデルごとに事後確率を計算する事後確率計算過程と、
確率分布モデルの各パラメタ値を更新する更新過程と、を有するパラメタ推定方法。
【請求項７】
請求項５または６記載のパラメタ推定方法の各過程と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定過程と、
前記音源該当情報が示す確率分布モデルごとのマスクを作成するマスク作成過程と、
前記周波数観測信号に前記マスクを乗算することで、分離信号を求める分離過程と、
前記分離信号を時間領域に変換する時間領域変換過程と、を有することを特徴とする音源分離方法。
【請求項８】
請求項５または６記載のパラメタ推定方法の各過程と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定過程と、
前記音源該当情報が示す確率分布モデルの平均を出力する方向出力過程と、を有することを特徴とする方向推定方法。
【請求項９】
請求項５または６記載のパラメタ推定方法、または請求項７記載の音源分離方法、または請求項８記載の方向推定方法の各過程をコンピュータに実行させるためのプログラム。
【請求項１】
複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成する周波数領域変換部と、
前記周波数観測信号の前記収音手段間の位相差を計算する位相差計算部と、
前記位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する推定部と、を有するパラメタ推定装置。
【請求項２】
請求項１記載のパラメタ推定装置であって、
前記推定部は、
現在の前記確率分布モデルの各パラメタを保持するパラメタ保持部と、
前記推定部は、前記位相差と、前記現在の確率分布モデルの各パラメタと、を用いて確率分布モデルごとに事後確率を計算する事後確率計算部と、
確率分布モデルの各パラメタ値を更新する更新部と、を有するパラメタ推定装置。
【請求項３】
請求項１または２記載のパラメタ推定装置と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定部と、
前記音源該当情報が示す確率分布モデルごとのマスクを作成するマスク作成部と、
前記周波数観測信号に前記マスクを乗算することで、分離信号を求める分離部と、
前記分離信号を時間領域に変換する時間領域変換部と、を有することを特徴とする音源分離装置。
【請求項４】
請求項１または２記載のパラメタ推定装置と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定部と、
前記音源該当情報が示す確率分布モデルの平均を出力する方向出力部と、を有することを特徴とする方向推定装置。
【請求項５】
複数の音源それぞれからの音源信号が混合され、２個の収音手段で収音された観測信号を周波数領域に変換することで周波数観測信号を生成する周波数領域変換過程と、
前記周波数観測信号の前記収音手段間の位相差を計算する位相差計算過程と、
前記位相差の分布に当てはまり、周波数依存性のある確率分布モデルのパラメタを推定する推定過程と、を有するパラメタ推定方法。
【請求項６】
請求項５記載のパラメタ推定方法であって、
前記推定過程は、
現在の前記確率分布モデルの各パラメタを保持するパラメタ保持過程と、
前記推定過程は、前記位相差と、前記現在の確率分布モデルの各パラメタと、を用いて確率分布モデルごとに事後確率を計算する事後確率計算過程と、
確率分布モデルの各パラメタ値を更新する更新過程と、を有するパラメタ推定方法。
【請求項７】
請求項５または６記載のパラメタ推定方法の各過程と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定過程と、
前記音源該当情報が示す確率分布モデルごとのマスクを作成するマスク作成過程と、
前記周波数観測信号に前記マスクを乗算することで、分離信号を求める分離過程と、
前記分離信号を時間領域に変換する時間領域変換過程と、を有することを特徴とする音源分離方法。
【請求項８】
請求項５または６記載のパラメタ推定方法の各過程と、
音源に該当する確率分布モデルを示す音源該当情報を求める有効音源推定過程と、
前記音源該当情報が示す確率分布モデルの平均を出力する方向出力過程と、を有することを特徴とする方向推定方法。
【請求項９】
請求項５または６記載のパラメタ推定方法、または請求項７記載の音源分離方法、または請求項８記載の方向推定方法の各過程をコンピュータに実行させるためのプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【公開番号】特開２０１０−１８７０６６（Ｐ２０１０−１８７０６６Ａ）
【公開日】平成２２年８月２６日（２０１０．８．２６）
【国際特許分類】
【出願番号】特願２００９−２８２７０（Ｐ２００９−２８２７０）
【出願日】平成２１年２月１０日（２００９．２．１０）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
【公開日】平成２２年８月２６日（２０１０．８．２６）
【国際特許分類】
【出願日】平成２１年２月１０日（２００９．２．１０）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
[ Back to top ]