(57)【要約】【課題】音情報の各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる。【解決手段】音源からの音情報の分布が複数のピークを持つ場合に、Ｍ個の確率分布モデルを用いて、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせるものであり、現在の確率分布モデルの各パラメタを保持し、音情報と、現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算し、音情報と、Ｍ個の確率分布モデルごとの事後確率を用いて、現在の確率分布モデルの各パラメタを更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタをパラメタ保持部に現在の確率分布モデルの各パラメタとして保持させ、各パラメタのうち、混合重みの事前分布にディリクレ分布を用いる。【選択図】図３
2010145836.tif 000002
【特許請求の範囲】
【請求項１】
音源からの音情報の分布が複数のピークを持つ場合に、Ｍ（Ｍは１以上の整数）個の確率分布モデルを用いて、当該確率分布の各パラメタを更新することで、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる方向情報分布推定装置であって、
現在の確率分布モデルの各パラメタを保持しているパラメタ保持部と、
前記音情報と、前記現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算する事後確率計算部と、
前記音情報と、前記Ｍ個の確率分布モデルごとの事後確率を用いて、前記現在の確率分布モデルの各パラメタを前記更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタを前記パラメタ保持部に前記現在の確率分布モデルの各パラメタとして保持させる更新部と、を備え、
前記更新部において、前記各パラメタのうち、混合重みの事前分布にディリクレ分布を用いる方向情報分布推定装置。
【請求項２】
請求項１記載の方向情報分布推定装置であって、
前記確率分布モデルは、正規分布モデルであり、
前記正規分布モデルの各パラメタは、混合重み、平均、分散、であることを特徴とする方向情報分布推定装置。
【請求項３】
請求項１記載の方向情報分布推定装置であって、
前記確率分布モデルは、フォン・ミーゼス分布モデルであり、
前記フォン・ミーゼス分布モデルの各パラメタは、混合重み、平均、拡散パラメタ、であることを特徴とする方向情報分布推定装置。
【請求項４】
複数の収音手段で入力された音信号を周波数領域に変換することで、周波数領域音信号を求める周波数領域変換部と、
前記周波数領域音信号から音の到来方向情報を求める到来方向推定部と、
前記周波数領域音信号のパワーを求めるパワー推定部と、
前記音の到来方向情報を音情報とし、前記パワーを重み係数として、方向情報分布モデルを求める請求項１〜３何れかに記載の方向情報分布推定装置と、
混合重みが予め定められた第１閾値よりも大きな値である音源該当方向情報分布モデルの個数Ｍ’を測定することで、音源数を求める音源数測定部と、を備える音源数推定装置。
【請求項５】
請求項４記載の音源数推定装置と、
各音源該当方向情報分布モデルの各パラメタのうち、平均を音源方向として出力する音源方向測定部と、を備える音源方向測定装置。
【請求項６】
請求項４記載の音源数推定装置と、
前記Ｍ’個の音源該当方向情報分布モデルごとの周辺化事後確率を求め、当該周辺化事後確率と前記周波数領域音信号とを掛け合わせることで、周波数領域目的信号を求める分離部と、
前記周波数領域目的信号を時間領域に変換することで、目的信号を求める時間領域変換部と、を備える音源分離装置。
【請求項７】
音源からの音情報の分布が複数のピークを持つ場合に、Ｍ（Ｍは１以上の整数）個の確率分布モデルを用いて、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる方向情報分布推定方法であって、
現在の確率分布モデルの各パラメタを保持しているパラメタ保持過程と、
前記音情報と、前記現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算する事後確率計算過程と、
前記音情報と、前記Ｍ個の確率分布モデルごとの事後確率を用いて、前記現在の確率分布モデルの各パラメタを更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタを前記パラメタ保持過程に前記現在の確率分布モデルの各パラメタとして保持させる更新過程と、を有し、
前記更新過程において、前記各パラメタのうち、混合重みの事前分布にディリクレ分布を用いる方向情報分布推定方法。
【請求項８】
請求項７記載の方向情報分布推定方法であって、
前記確率分布モデルは、正規分布モデルであり、
前記正規分布モデルの各パラメタは、混合重み、平均、分散、であることを特徴とする方向情報分布推定方法。
【請求項９】
請求項７記載の方向情報分布推定方法であって、
前記確率分布モデルは、フォン・ミーゼス分布モデルであり、
前記フォン・ミーゼス分布モデルの各パラメタは、混合重み、平均、拡散パラメタ、であることを特徴とする方向情報分布推定方法。
【請求項１０】
複数の収音手段で入力された音信号を周波数領域に変換することで、周波数領域音信号を求める周波数領域変換過程と、
前記周波数領域音信号から音の到来方向情報を求める到来方向推定過程と、
前記周波数領域音信号のパワーを求めるパワー推定過程と、
前記音の到来方向情報を音情報とし、前記パワーを重み係数として、方向情報分布モデルを求める請求項７〜９何れかに記載の方向情報分布推定方法の各過程と、
混合重みが予め定められた第１閾値よりも大きな値である音源該当方向情報分布モデルの個数Ｍ’を測定することで、音源数を求める音源数測定過程と、を有する音源数推定方法。
【請求項１１】
請求項１０記載の音源数推定方法の各過程と、
各音源該当方向情報分布モデルの各パラメタのうち、平均を音源方向として出力する音源方向測定過程と、を有する音源方向測定方法。
【請求項１２】
請求項１０記載の音源数推定方法の各過程と、
前記Ｍ’個の音源該当方向情報分布モデルごとの周辺化事後確率を求め、当該周辺化事後確率と前記周波数領域音信号とを掛け合わせることで、周波数領域目的信号を求める分離過程と、
前記周波数領域目的信号を時間領域に変換することで、目的信号を求める時間領域変換過程と、を有する音源分離方法。
【請求項１３】
請求項７〜９何れかに記載の方向情報分布推定方法、または請求項１０記載の音源数推定方法、または請求項１１記載の音源方向測定方法、または請求項１２記載の音源分離方法、の各過程をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、音響信号処理に用いる方向情報分布推定装置、音源数推定装置、音源方向測定装置、音源分離装置、それらの方法、それらのプログラムに関する。
【背景技術】
【０００２】
従来から、音源の方向情報の分布を推定する技術がある。この技術は、例えば、音響信号処理の分野において、同時に複数の人が発話した音声が混ざった信号が観測されている時に、各人の方向を推定したり、各人の音声を分離抽出する際に重要である。
【０００３】
図１に従来の方向情報分布推定装置１００の機能構成例を示し、図２Ａに音の方向情報について得られたヒストグラムＨの一例を示し、当該方向情報は方向情報分布推定装置１００に入力されるものである。方向情報分布推定装置１００の原理については非特許文献１に記載されている。図２Ａ、後述する図２Ｂ、図２Ｃの横軸は、音源の（音の到来方向）の角度を示す。方向情報分布推定装置１００の目的は、ヒストグラムＨに、正規分布モデルをフィッティングさせる（近似させる）ことである。特にこの技術では、ヒストグラムＨ中にある複数の分布の山それぞれに意味がある場合を考える。例えば、それぞれの分布の山が音声信号源や電波信号源などの推定方向情報を表している場合などである。図２Ａの例では、分布の山は４つ存在し、４つの分布の山をそれぞれａ〜ｄとする。分布の山ａは約−１１５度、分布の山ｂは約−２０度、分布の山ｃは約６０度、分布の山ｄは約１５０度に位置している。
【０００４】
ヒストグラムＨから音源方向を推定する場合には、分布の山の角度の平均値を求めるのであるが、この平均値を求めるためには、各分布の山をそれぞれ１個の確率分布モデルでフィッティングさせることが要求される。
【０００５】
従来の技術では、ヒストグラム全体を例えば混合正規分布モデル（Ｇａｕｓｓｉａｎ  ｍｉｘｔｕｒｅ  ｍｏｄｅｌ：ＧＭＭ）でモデル化していた。混合正規分布Ｇは以下の式（１）により表される。
【０００６】
【数１】
2010145836.tif 000003
【０００７】
図１を用いて、従来の方向情報分布推定装置１００の各構成部の処理を簡単に説明する。まず、方向情報（図２Ａに示すヒストグラムＨ）が事後確率計算部２に入力され、事後確率計算部２でＭ個の正規分布モデルごとに事後確率を求める。そして、Ｍ個の事後確率はパラメタ更新部４に入力される。平均更新手段４２、分散更新手段４４、混合重み更新手段４６がそれぞれ、Ｍ個の事後確率等を用いて、パラメタ保持部８に保持させつつ、μｍ、σｍ、αｍをＥＭアルゴリズムにより更新する。更新過程の詳細は省略する。
【０００８】
そして、例えば、更新回数が閾値Ｔを超えると、収束後パラメタθｖ（θｖ＝（μｍ、σｍ、αｍ）  ｍ＝１、．．．、Ｍ）を出力する。出力された収束後パラメタθｖを用いて、音源方向の測定や音源の分離を行う。
【非特許文献１】M.Mandel,D.Ellis,and T.Jebara,"An EM algorithm for localizing multiple sound sources in reverberant environments,"Proc.Neural Info.Proc.Sys.,2006.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
図２Ａに示すヒストグラムＨに８つの正規分布モデルを用いて、方向情報分布推定装置１００によるフィッティングさせた結果を図２Ｂに示す。図２Ｂの−１１５度付近の箇所（図２Ｂ記載のＰの箇所）に注目されたい。Ｐの箇所では、−１１５度付近の分布も１つの正規分布モデルでフィッティングさせたいにも関わらず、２つの正規分布モデルがフィッティングしてしまう。そして、図２Ｂ記載の正規分布モデルをそれぞれ合計することで、図２Ｃ記載の混合正規分布モデルが求められ、求められた混合正規分布モデルの収束後パラメタθｖが方向情報分布推定装置１００から出力される。この場合であると、Ｐに２つの正規分布モデルがフィッティングしてしまっていることから、正確な収束後パラメタθｖを求めることができないという問題がある。その結果、正確な音源分離や音源方向推定を行うことができない。
【００１０】
本発明では、音源からの音情報として、複数のピークを持つ方向情報が与えられた場合に、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせることのできる方向情報分布推定装置を提供することである。
【課題を解決するための手段】
【００１１】
この発明は音源からの音情報が複数のピークを持つ場合に、Ｍ（Ｍは１以上の整数）個の確率分布モデルを用いて、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる方向情報分布推定装置である。当該方向情報分布推定装置は、パラメタ保持部と、事後確率計算部と、更新部と、を備える。パラメタ保持部は、現在の確率分布モデルの各パラメタを保持している。事後確率計算部は、音情報と、現在の確率分布モデルの各パラメタを用いて、Ｍ個の確率分布モデルごとに事後確率を計算する。更新部は音情報と、Ｍ個の確率分布モデルごとの事後確率を用いて、現在の確率分布モデルの各パラメタを更新し、各パラメタ値が収束していると判断した場合には更新された各パラメタを出力し、各パラメタ値が収束していないと判断した場合には、更新された各パラメタをパラメタ保持部に現在の確率分布モデルの各パラメタとして保持させる。そして、各パラメタのうち、混合重みの事前分布にディリクレ分布を用いる。
【発明の効果】
【００１２】
この発明の方向情報分布推定装置では、確率分布モデルのパラメタである混合重みに事前分布としてディリクレ分布を与えることで、各ピークに対して少数の確率分布モデルのフィッティングが可能であり、結果として各ピークにそれぞれ１つの確率分布モデルをフィッティングさせることができる。
【発明を実施するための最良の形態】
【００１３】
以下に、発明を実施するための最良の形態を示す。なお、同じ機能を持つ構成部や同じ処理を行う過程には同じ番号を付し、重複説明を省略する。
【実施例１】
【００１４】
図３に実施例１の方向情報分布推定装置２００の機能構成例を示し、図４に処理フローを示す。図５Ａに入力されるヒストグラムＨを示し、図５Ｂ、Ｃにそれぞれ、方向情報分布推定装置２００で得られる分布の一例を示す。また、図５ＡのヒストグラムＨは図２Ａと同様であるとする。
【００１５】
この実施例１では、用いるＭ個の確率分布モデルとしてＭ個の正規分布を用いる例を示し、入力される音情報を方向情報Ｄとし、方向情報Ｄの一例としてヒストグラムＨである場合を示す。音情報が方向情報である場合には、横軸は角度、縦軸は度数を表す。そして、方向情報分布推定装置２００はヒストグラムＨが複数のピークを持つ場合に、各ピークにそれぞれ１つの確率分布モデルをフィッティングさせる。正規分布モデルでなくとも、確率分布モデルであれば、どのモデルを用いても良い。ここで、各ピークは、１つの音源方向に対応するものである。
通常、方向情報Ｄには、２πのｋ（ｋは整数）倍の不定性が含まれるため、ここではそれを許すＷｒａｐｐｅｄ  ＧＭＭを用いる。Ｗｒａｐｐｅｄ  ＧＭＭであるＧは以下の式（２）で表すことができる。
【００１６】
【数２】
2010145836.tif 000004
【００１７】
θを混合正規分布の平均μ、分散σ、混合重みαをまとめたものを示し、つまり、θ＝（μｍ、σｍ、αｍ）＝（μ１、σ１、α１、．．．、μｍ、σｍ、αｍ、．．．、μＭ、σＭ、αＭ）となる。また、ｔを更新回数（時刻）とし、θに更新回数の概念を付与したもの、つまり、ｔ回更新したθをθｔとすると、θｔ＝（μ１ｔ、σ１ｔ、α１ｔ、．．．、μｍｔ、σｍｔ、αｍｔ、．．．、μＭｔ、σＭｔ、αＭｔ）となる。また、記憶部１６には予め用いる正規分布モデルのモデル数Ｍと混合正規分布モデルの各パラメタの初期値θ０が記憶されている。事前分布情報保持部１１０には、後述のハイパーパラメタφ、重みパラメタｃが保持されている。
【００１８】
方向情報分布推定装置２００には、Ｎ個の方向情報Ｄ＝｛ｄ１、．．．、ｄｎ、．．．、ｄＮ｝と重み係数Ａ＝｛ａ１、．．．、ａｎ、．．．、ａＮ｝が入力される。重み係数Ａは、方向情報の各要素ｄｎ（ｎ＝１、．．．、Ｎ）に対する重み係数である。この重み係数は例えば、方向情報Ｄが得られる頻度や方向情報Ｄが得られた時の信頼度（取得信号のパワーや信号の瞬時的ＳＮ比など）により与えることができる。または全てのｎについてａｎ＝１としても良い。
【００１９】
まず、ｔ＝０と設定し（つまり更新回数が０）、ｔ＝０のときの混合正規分布のパラメタθ０の値を設定し、用いる正規分布のモデル数Ｍ、ｋの範囲であるＫ、更新回数閾値Ｔまたは差閾値Δ（後述する）を設定する。更新回数閾値Ｔまたは差閾値Δは、後述する収束判定処理の際に用いられる（ステップＳ２）。
【００２０】
事後確率計算部１２は、音情報（この実施例１では方向情報Ｄ）と、現在の混合正規分布のパラメタθｔ（＝（μｍｔ、σｍｔ、αｍｔ  ｍ＝１、．．．、Ｍ））から、Ｍ個の正規分布ごとに事後確率ｐ（ｍ、ｋ│ｄｎ、θｔ）を計算する（ステップＳ６）。またパラメタ保持部１８には、現在の混合正規分布のθｔが保持されている。事後確率計算部１２は具体的には例えば、以下の式（５）により計算する。
【００２１】
【数３】
2010145836.tif 000005
【００２２】
式（５）の右辺の分子「ｐ（ｍ、ｋ、ｄｎ│θｔ）は上記式（３）（４）で表される「ｐ（ｍ、ｋ、ｄｎ│θ）」内のθに更新回数ｔの概念を付与したものである。
【００２３】
次に、更新部１４は、方向情報Ｄと事後確率ｐ（ｍ、ｋ│ｄｎ、θｔ）を用いて、現在の混合正規分布の各パラメタθｔを更新する（ステップＳ８）。以下、更新処理について詳細に説明する。更新部１４は更新処理の際に、ハイパーパラメタφ、重みパラメタｃを事前分布情報保持部１１０から取り出す。この発明では、パラメタθの更新処理は、正規分布のパラメタθの混合重みαｍに適切な事前分布を与え、例えばＥＭアルゴリズムにて行う。この実施例１では、混合重みαｍの事前分布として、ディリクレ分布を考える。ディリクレ分布の詳細は、参考文献１である「Ｃ．Ｍ．ビショップ著（元田、栗田他訳）  「パターン認識と機械学習（上）」、シュプリンガー・ジャパン２００７年  ｐ．７４−ｐ．７７」等に記載されている。ディリクレ分布は例えば以下の式（６）で表される。
【００２４】
【数４】
2010145836.tif 000006
【００２５】
ここで、αは混合重み行列であり、α＝｛α１、．．．、αｍ、．．．、αＭ｝で表され、ΣｍＭαｍ＝１、０≦αｍ≦１という条件を満たす。これは混合正規分布のパラメタである混合重みの条件と同じであることに注意されたい。またβ（φ）は正規化項（ベータ分布）であり、ここで、ハイパーパラメタφを１より小さい正の値（例えば、０．９）に設定すると、αｍのごく少数のみが十分に大きな値を持ち、残りは０に近い値をとるようになる。この性質を式（１）で表される混合正規分布Ｇの混合重みαｍに対して適用することで、混合正規分布Ｇのうちの少数の正規分布のみに十分大きな混合重みがかかり、その他の正規分布の混合重みは０に近くなる。結果として、なるべく少数の正規分布によるフィッティングが可能である。
【００２６】
次に、この事前分布を含みながら、パラメタ更新を行うためのＥＭアルゴリズムを導出する。まず、最尤推定のためのコスト関数Ｌ（θ）は次のように与えられる。
【００２７】
【数５】
2010145836.tif 000007
【００２８】
また、重みパラメタｃは、式（９）の第１項と第２項の重みをコントロールするパラメタである。
【００２９】
【数６】
2010145836.tif 000008
【００３０】
となる。ここで、式（１１）のE［H］は式Hの期待値を示し、式（１２）中のｐ（ｍ、ｋ│ｄｎ、θｔ）は式（５）で表される事後確率分布である。ここで、従来法のＥＭアルゴリズムでは式（１２）中のｌｏｇ（ｐ（α））がないことに注意されたい。
【００３１】
【数７】
2010145836.tif 000009
【００３２】
また上述の通り、この場合には、音情報は音の到来方向を示すＮ個の方向情報ｄｎ（ｎ＝１、．．．、Ｎ）であり、当該方向情報ｄｎには２ｋπ倍（ｋは整数）の不定性が含まれているとし、混合正規分布は、ラップＧＭＭであり、ｃは重みパラメタであり、φはハイパーパラメタであり、Ｋはｋの範囲を示す。
【００３３】
図３中の更新部１４中の平均更新手段１４２が式（１３）より現在の平均μｍｔを更新することで更新後の平均μｍｔ＋１を出力する。分散更新手段１４４が式（１４）より分散σｍを更新することで更新後の分散σｍｔ＋１を出力する。混合重み更新手段１４６が式（１５）により混合重みαｍを更新することで更新後の混合重みαｍｔ＋１を出力する。パラメタ算出手段が、更新後の平均μｍｔ＋１、分散σｍｔ＋１、混合重みαｍｔ＋１についての更新後のパラメタθｔ＋１を算出する。
【００３４】
各パラメタの更新処理が数回行われ（ステップＳ４）、更新部１４内の収束判定手段１５０は、更新されたθｔ＋１に対して、予め定められた規則により、各パラメタ値が収束しているか否かの収束判定を行う（ステップＳ１０）。各パラメタ値が収束していると判断した場合には、更新されたパラメタθｔ＋１を出力する。また、各パラメタ値が収束していないと判断した場合には、更新されたパラメタθｔ＋１を現在の確率分布モデルの平均、分散、混合重みとしてパラメタ保持部１８に保持させる。そして、収束判定手段１５０が、各パラメタ値が収束していると判断するまで、ステップＳ４〜ステップＳ１０の処理を繰り返す。
【００３５】
ここで収束判定に用いる予め定められた規則の例を説明する。更新回数閾値Ｔを用いる例を説明すると、更新部１４内のカウント手段（図示せず）は更新回数ｔをカウントし、更新回数ｔが更新回数閾値Ｔを超えた場合には、十分更新しており、収束していると判断して、更新後のパラメタθを出力する。また、差閾値Δを用いる例を説明すると、以下の式（１６）の式を満たす場合には、収束していると判断して、更新後のパラメタθを出力する。
│Ｑ（θ│θｔ＋１）−Ｑ（θ│θｔ）│＜Δ      （１６）
パラメタ算出手段１４８により算出されるθは図５Ｃの混合正規分布の各パラメタである。
【００３６】
また、この実施例１において、Ｋ＝０とすれば、ラップＧＭＭではなく、通常のＧＭＭによるフィッティングを行うことができる。この場合は、音情報Ｄとして方向情報である必要はない。例えば、音源からの音をＪ個のマイクロホン２０ｊ（ｊ＝１、．．．、Ｊ）で収音した場合であると、マイクロホン２０ｊとマイクロホン２０ｊ’ （ｊ’＝１、．．．、Ｊであり、ｊ≠ｊ’）とのマイクロホン間位相差ｑ’ｊｊ’を音情報Ｄとしても良い。またこの実施例１では混合重みαｍのみに事前分布を導入したが、各ガウス分布の平均μｍと分散σｍに対しても事前分布を導入することで、より精度の高いＧＭＭフィッティング（方向情報分布推定処理）を実現できる。また、各ガウス分布の各パラメタである平均μｍ、分散σｍ、混合重みαｍ、に事前分布を導入した場合にのＧＭＭフィッティングには、ＥＭアルゴリズムのほかにもベイズ推定などの様々なアルゴリズムが知られている。これらの拡張は当業者であれば、上記参考文献１などを参照すれば、容易に実現できるため、ここでは省略する。
【００３７】
この実施例１で説明したように、式（５）のハイパーパラメタφを１より小さな正の値（例えば、０．９）に設定すると、ディリクレ分布の性質からαｍのごく少数のみが十分に大きな値を持ち、残りは０に近い値をとるようになる。上記式（１）に示すＧＭＭの少数の正規分布のみに十分大きな混合重みαｍがかかり、そのほかの正規分布の重みは０に近くなる。この性質を用いることにより、なるべく少数の正規分布によるモデルフィッティングが可能になる。
【実施例２】
【００３８】
実施例１では、Ｍ個の確率分布モデルとして、Ｍ個の正規分布モデルを用いたが、実施例２では、Ｍ個のフォン・ミーゼス（ｖｏｎ  Ｍｉｓｅｓ）分布モデルを用いる。フォン・ミーゼス分布は角度の分布を表す関数であり、フォン・ミーゼス分布モデルの詳細は、参考文献２「Ｋ．Ｖ．Ｍａｒｄｉａ、”Ｓｔａｔｉｓｔｉｃｓ  ｏｆ  Ｄｉｒｅｃｔｉｏｎａｌ  Ｄａｔａ”、Ａｃａｄｅｍｉｃ  Ｐｒｅｓｓ、１９７２、３．４．９節」などに記載されている。フォン・ミーゼス分布を用いる効果は、正規分布モデルを用いた場合と比較すると、ｋおよびＫの値を考慮する必要がないため、演算処理が削減されることである。
【００３９】
この実施例２の方向情報分布推定装置３００の機能構成例、処理フローは図３、図４とほぼ同様であるが、図３中の分散更新手段１４４が拡散パラメタ更新手段１６０に代替されている点が異なる。以下、詳細に説明する。また、フォン・ミーゼス分布モデルのパラメタθをθ＝｛μｍ、кｍ、αｍ｝とし、кｍは拡散パラメタである。
【００４０】
まず、事後確率計算部１２は音情報Ｄ（例えば方向情報Ｄ）とパラメタ保持部１８に保持されている現在のパラメタθｔ＝｛αｍｔ、μｍｔ、кｍｔ｝からＭ個それぞれのフォン・ミーゼス分布モデルに関する事後確率ｐ（ｍ│ｄｎ、θｔ）を求める。
【００４１】
【数８】
2010145836.tif 000010
【００４２】
この式（１７）は式（４）と対応しているものであり、式（１７）中の右辺の分子ｐ（ｍ、ｄｎ│θｔ）は、フォン・ミーゼス分布ｇ（ｄｎ；μｍ、кｍ）である。
【００４３】
【数９】
2010145836.tif 000011
【００４４】
ここで、−π＜ｄｎ≦π、−π＜μｍ≦πとし、また、кｍ＞０である。また、Ｉ０（ｘ）は０次の第１種の変形されたベッセル関数である。
【００４５】
次に、更新部１４は、音情報Ｄと事後確率ｐ（ｍ│ｄｎ、θｔ）を用いて、フォン・ミーゼス分布のパラメタθ、つまり、平均μｍｔ、拡散パラメタкｍｔ、混合重みαｍｔを更新する。以下、詳細に説明する。
【００４６】
平均更新手段１４２は平均μｍｔを例えば以下の式（２０）により更新する。
【００４７】
【数１０】
2010145836.tif 000012
【００４８】
ここで、ａｒｃｔａｎ（ｘ）は−π／２＜μｍ＜π／２の値を返すのが一般的であるから、−π＜μｍ＜πのデータを扱うには、以下の演算も行う。
式（２０）の値が負の場合、μｍｔとμｍｔ＋πの両方について、式（２１）に示すＱ関数の２次導関数を計算し、式（２１）の値が負になるほうをμｍｔ＋１とする。
【００４９】
【数１１】
2010145836.tif 000013
【００５０】
式（２０）の値が正の場合、μｍｔ、μｍｔ−πについて式（２１）を計算し、これが負になる方をμｍｔとして保存する。
拡散パラメタ更新手段１６０は例えば以下の式（２２）により更新する。
【００５１】
【数１２】
2010145836.tif 000014
【００５２】
ここで、Ｉ（кｍｔ＋１）を拡散パラメタ関数とする。кｍｔ＋１は解析的に得られないが、次のように得ることができる。拡散パラメタ関数Ｉ（кｍｔ＋１）は、単調増加関数である。そこで、ある範囲のк（例えば、０≦к≦１００）について、「кｍｔ＋１」と「Ｉ（кｍｔ＋１）」とを対応させたルックアップテーブルを用意しておく。当該ルックアップテーブルは、拡散パラメタ更新手段１６０中の記憶部（図示せず）に記憶させておけばよい。そして、Ｉ（кｍｔ＋１）が求まると、ルックアップテーブルを参照して、Ｉ（кｍｔ＋１）に対応するкｍｔ＋１を出力する。
混合重み更新手段１４６は、例えば、以下の式（２３）により混合重みαｍを更新する。
【００５３】
【数１３】
2010145836.tif 000015
【００５４】
このようにして、更新部１４は分布パラメタθｔ（＝｛αｍ、μｍ、кｍ｝）を更新する。
この実施例２の方向情報分布推定装置３００のように、フォン・ミーゼス分布を用いることで、ｋに関する推定操作が不要であるため、実施例１の方向情報分布推定装置２００と比較して、計算コストやパラメタθの収束時間を削減できる。
【００５５】
［実験結果１］
図５を用いて、実施例１で説明した方向情報分布推定装置２００によるフィッティングの実験結果について説明する。実験条件として８（＝Ｍ）つの正規分布からなる混合正規分布をフィッティングさせ、ハイパーパラメタφを０．９とする。上述のように図５Ａに入力される方向情報ｄｎについての図２Ａと同様のヒストグラムＨを示し、図５Ｂに方向情報分布推定装置２００のフィッティング処理による正規分布の結果を示し、図５Ｃに図５Ｂの正規分布を合計した混合正規分布（ＧＭＭ）を示す。図５Ｂ記載のＰの箇所（−１１５度付近）に注目すると、１つの正規分布でフィッティングできていることが理解されよう。従って、図５Ｃに示す求められる混合正規分布は、正確なものである。従って、実施例３〜５で説明する音源数推定処理、音源方向測定処理、音源分離処理も正確に行うことができる。
【００５６】
一方、上述のように、図１Ｂに示す従来の方向情報分布推定装置１００の実験結果については、Ｐの箇所では、２つの正規分布がフィッティングしてしまい、図１Ｃに示すＧＭＭは、不正確なものとなってしまう。
【実施例３】
【００５７】
この実施例３では、実施例１、２で説明した方向情報分布推定装置２００、３００を用いた音源数測定装置４００について説明する。図６に音源数測定装置４００の機能構成例を示す。この実施例３の音源数測定装置４００は、Ｊ（Ｊは２以上の整数）個の収音手段２０ｊ（例えば、マイクロホン  ｊ＝１、．．．、Ｊ）に接続されている場合を説明する。そして、ある収録時間内（例えば５秒間など）に複数の音源から音が発せられた場合に、当該音をＪ個の収音手段２０ｊで収録したとする（以下、状況Ｘという。）。この実施例３の音源数測定装置４００は、収録音のみを用いて、音を発した音源の数を推定する。
【００５８】
収音手段２０ｊから入力された音信号をｘｉ（ｓ）とし、ｓを離散時刻とする。周波数領域変換部３０は音信号ｘｉ（ｓ）を周波数領域音信号Ｘｊ（ｆ、τ）に変換する。ｆは周波数、τは時間フレーム番号である。また、この実施例３ではｎ＝τＦ＋ｆを考える。ただしＦは周波数領域の数である。
パワー推定部３２は周波数領域音信号から音のパワーを求める。求め方の一例として、パワー推定部３２は、各時間周波数（ｆ、τ）における周波数領域音信号Ｘｊ（ｆ、τ）の信号パワー│Ｘｊ（ｆ、τ）│２を演算し出力する。出力された信号パワー│Ｘｊ（ｆ、τ）│２が、上述した重み係数ａｎとして、以後用いられる。
【００５９】
また、到来方向推定部３４は周波数領域音信号から音の到来方向情報を求める。求め方の例を詳細に説明する。到来方向推定部３４は、収音手段間位相差演算手段３４２、到来方向情報生成手段３４４とで構成されている。まず、収音手段間位相差演算手段３４２が各フレームτ、各周波数ｆにおいて、各収音手段の全ての組み合わせ（マイクロホンペア）について収音手段間位相差ｑ’ｊｊ’（ｆ、τ）を以下の式（２４）により求める。ただし、ｊ＝１、．．．、Ｊであり、ｊ’＝１、．．．、Ｊであり、ｊ≠ｊ’とする。
ｑ’ｊｊ’（ｆ、τ）＝｛ａｒｇ［Ｘｊ（ｆ、τ）Ｘ＊ｊ’（ｆ、τ）］｝／２πｆ
（２４）
【００６０】
ただし、「＊」は複素共役であることを示す。そして、全てのｑ’ｊｊ’（ｆ、τ）を並べたベクトルをＱ’（ｆ、τ）とする。音の到来方向情報Ｑ（ｆ、τ）は音速Ｃと、各収音手段の座標系Ｄを用いて、以下の式（２５）により求められる。
Ｑ（ｆ、τ）＝ＣＤ＋Ｑ’（ｆ、τ）      （２５）
【００６１】
ここでＣは音速であり、「＋」は、Ｍｏｏｒｅ−Ｐｅｎｒｏｓｅの擬似逆行列を表し、Ｄ＝［Ｄ１−ＤＬ、．．．、Ｄｊ−ＤＬ、．．．、ＤＪ−ＤＬ］Ｔであり、Ｄｊは収音手段２０ｊの座標（ｘ、ｙ、ｚ）と並べたベクトルであり、ＬはＪ個の収音手段のうち代表として選ばれた代表収音手段のインデックスである。到来方向情報Ｑ（ｆ、τ）のｘｙｚ座標（ｘＱ、ｙＱ、ｚＱ）は、到来方向水平角（以下、単に「水平角」という。）をΨ（ｆ、τ）とし、到来方向仰角（以下、単に「仰角」という。）をΩ（ｆ、τ）とすると、以下の式（２６）で表すことができる。
Ｑ（ｆ、τ）＝（ｘＱ、ｙＱ、ｚＱ）
＝（ｃｏｓΨ（ｆ、τ）ｃｏｓΩ（ｆ、τ）、
ｓｉｎΨ（ｆ、τ）ｃｏｓΩ（ｆ、τ）、
ｓｉｎΩ（ｆ、τ））                       （２６）
【００６２】
この実施例では、水平角Ψ（ｆ、τ）のみを用いる。求められた到来方向情報Ｑ（ｆ、τ）を方向情報ｄｎとして用いる。また、方向情報ｄｎについてヒストグラムを作成すると、図２Ａに示すヒストグラムＨが得られる。次に、方向情報ｄｎ、方向情報ａｎは方向情報分布推定装置２００（または３００）に入力され、実施例１（または実施例２）で説明した処理により、パラメタθが出力される。以降、出力されたパラメタθを決定後パラメタθとする。
【００６３】
音源数測定部３６は、決定後パラメタθの混合重みαｍ（ｍ＝１、．．．、Ｍ）のうち、混合重みが予め定められた第１閾値ε１（例えば１０−６）よりも大きな値である方向情報分布モデルの個数Ｍ’を測定する。測定された個数Ｍ’を音源数として出力する。何故なら、方向情報分布推定装置２００（３００）の演算が十分収束している場合には、決定後パラメタθの中の混合重みαｍのうち十分大きな値を持つ個数はヒストグラム中の分布の山の数と等しくなるからである。以下の説明では、音源と認められたものについての方向情報分布モデルを音源該当方向情報分布モデルという。
【００６４】
また、方向情報分布推定装置２００（３００）の演算が十分に収束していない場合は、音源数測定部３６は、次のような推定処理を行うことが好ましい。まず、音源数測定装置４００内で、方向情報分布推定装置２００（実施例１で説明）を用いた場合には、音源数測定部３６は、混合重みαｍが第１閾値ε１よりも大きく、かつ分散σｍが予め定められた第２閾値ε２（例えば１５度）よりも小さい方向情報分布モデルを音源該当方向情報分布モデルとして検出し、これら検出された音源該当方向情報分布モデルの個数Ｍ’を測定すればよい。また、音源数測定装置４００内で、方向情報分布推定装置３００（実施例２で説明）を用いた場合には、音源数測定部３６は、混合重みαｍが第１閾値ε１よりも大きく、かつ拡散パラメタкｍが第３閾値（例えば１０）よりも大きい方向情報分布モデルを音源該当方向情報分布モデルとして検出し、これら検出された音源該当方向情報分布モデルの個数Ｍ’を測定すればよい。
【００６５】
従来の方向情報分布推定装置１００は、ヒストグラムの各ピークに対して、正規分布をフィッティングさせると、図２Ｂに示すように、１つのピークに対して、２つの正規分布をフィッティングさせる場合がある。従って、方向情報分布推定装置１００を適用した音源数測定装置であると、誤った音源数測定をしてしまう。しかし、実施例１、２で説明した方向情報分布推定装置２００（または３００）により、図５Ｂに示すように、１つのピークに対して１つの確率分布モデル（例えば、正規分布モデルやフォン・ミーゼス分布）をフィッティングさせることができるので、正確な音源数を測定できる。
【実施例４】
【００６６】
この実施例４では、音源方向測定装置５００について説明する。音源方向測定装置５００は、状況Ｘの場合に、収録音のみを用いて音源の方向を推定する。図７に音源方向測定装置５００の機能構成例を示す。図７の例では、音源方向測定装置５００は音源数測定装置４００（実施例３で説明）と音源方向測定部３８とで構成されている。
【００６７】
音源数測定装置４００の処理が終了すると、音源方向測定部３８は、音源該当方向情報分布モデルのインデックスｍ’｛ｍ’＝１、．．．、Ｍ’｝に対応する平均パラメタμｍ’を方向情報分布推定装置２００から取り出し、推定すべき音源方向として当該平均パラメタμｍ’を出力する。
【００６８】
この実施例４のように、音源方向測定装置５００内に具備する方向情報分布推定装置２００（３００）により、正確な方向情報分布処理がされることから、音源方向測定装置５００は正確な音源方向測定を行うことができる。
【実施例５】
【００６９】
この実施例５では、音源分離装置６００について説明する。音源分離装置６００は、状況Ｘの場合に、収録音のみを用いて音源からの音信号を分離抽出する。図８に音源分離装置６００の機能構成例を示す。図８では、音源分離装置６００は音源数測定装置４００（実施例３で説明）と分離部４０、時間領域変換部４１とで構成されている。
【００７０】
音源数測定装置４００の処理が終了すると、分離部４０は、音源数測定装置４００で定められた音源該当方向情報分布モデルのインデックスｍ’について以下の処理を行う。
【００７１】
音源数測定装置４００が、方向情報分布推定装置２００を具備している場合には、分離部４０は、以下の式（２７）により、Ｍ’個の正規分布(式(５)参照)に関する事後確率ｐ（ｍ’、ｋ│ｄｎ、θｔ）を周辺化することで、周辺化事後確率ｐ（ｍ’│ｄｎ、θｔ）を求める。
ｐ（ｍ’│ｄｎ、θ）＝Σｋ＝−ＫＫｐ（ｍ’、ｋ│ｄｎ、θｔ）      （２７）
また、音源数測定装置４００が、方向情報分布推定装置３００を具備している場合には、周辺化処理を行わず、上記式（１７）の演算結果を用いる。
【００７２】
また、周波数領域変換部３０からの周波数領域音信号Ｘｊ（ｆ、τ）は分離部４０に入力される。分離部４０は周辺化事後確率と周波数領域音信号とを掛け合わせる。つまり、以下の式（２８）を演算することで、ｍ’番目の信号の推定に対応する周波数領域目的信号（分離信号）を出力する。
Ｙｎｍ’＝Ｘｎｐ（ｍ’│ｄｎ、θ）      （２８）
ここでＸｎはＸ１（ｆ、τ）を上述のｎ＝τＦ＋ｆで変形したものである。出力される周波数領域目的信号は時間周波数表現（ｆ、τ）を用いると、ｎ＝τＦ＋ｆより以下の式（２９）で表される事に留意されたい。
Ｙｍ’（ｆ、τ）＝Ｘ１（ｆ、τ）ｐ（ｍ’│Ψ（ｆ、τ）、θ）    （２９）
【００７３】
そして、時間領域変換部４１は周波数領域目的信号Ｙｍ’（ｆ、τ）を時間領域に変換することで、目的信号ｙｍ’（ｔ）を求め、出力する。
この実施例５のように、音源分離装置６００内に具備する方向情報分布推定装置２００（３００）により、正確な方向情報分布処理がされることから、音源分離装置６００は正確な信号分離を行うことができる。
【００７４】
［実験結果２］
次に、方向情報分布推定装置２００を用いた、音源数測定装置４００（実施例３で説明）と音源分離装置６００（実施例５で説明）と（以下、「発明法」という。）、従来の方向情報分布推定装置１００を用いた音源数測定装置、音源分離装置（以下、「従来法」という。）とを比較した実験結果について説明する。まず図９を用いて、実験条件について説明する。長手方向４．４５ｍ（＝Ｌｂ）、短手方向３．５５ｍ（＝Ｌａ）、高さ２．５ｍの室内に、３つのマイクロホンＺ１、Ｚ２、Ｚ３が、正三角形の各頂点に配置される。隣接するマイクロホン同士の間隔は４ｃｍであり、３つのマイクロホンの収音面がそれぞれ外側に向けられる。３つのマイクロホンＺ１、Ｚ２、Ｚ３がなす正三角形の重心は、図９の左下の頂点Ｘから長手方向に２．５６ｍ（＝Ｌｄ）であり、短手方向に１．８ｍ（＝Ｌｃ）の箇所に位置する。また、３つのマイクロホンＺ１、Ｚ２、Ｚ３を囲むように２〜４つのスピーカ（図９の例では４つのスピーカＳ１、Ｓ２、Ｓ３、Ｓ４とする。）が円周Ｒの方向に配置され、音を発しているとする。当該円周Ｒの半径は、５０ｃｍまたは１１０ｃｍであり、音の反響時間は１２８ｍｓであるとする。マイクロホンＺ１、Ｚ２、Ｚ３、スピーカＳ１、Ｓ２、Ｓ３、Ｓ４の高さは全て１．２ｍとする。
【００７５】
実験項目については、（１）音源（スピーカ）が２個、３個、４個の場合について音源の数を測定できるか（音源数測定処理）、（２）音源からの音信号を分離できるか（音源分離処理）、である。これらの項目について、スピーカから発せられる音の音質を変えたり、スピーカがなす円周Ｒの半径を変えるなどして、２０通りの組み合わせについて実験を行った。
【００７６】
このような条件下で、図１０に実験結果を示す。図１０では、音源数測定処理については、２０通りのうち、どの程度の確率で音源数Ｗを判定しているか評価し、音源分離処理については信号対妨害音比（Ｓｉｇｎａｌ  ｔｏ  Ｉｎｔｅｒｆｅｒｅｎｃｅ  Ｒａｔｉｏ：ＳＩＲ）を評価した。図１０からも理解されるように、音源数処理、音源分離処理については従来法では、誤った結果を出力しているが、発明法ではほぼ正確な結果を出していることが理解されよう。
【００７７】
＜ハードウェア構成＞
本発明は上述の実施の形態に限定されるものではない。また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
また、上述の構成をコンピュータによって実現する場合、方向情報分布推定装置２００（３００）、音源数推定装置４００、音源方向測定装置５００、音源分離装置６００が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、処理機能がコンピュータ上で実現される。
【００７８】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（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)等である。また、バスは、ＣＰＵ、入力部、出力部、補助記憶装置、ＲＡＭ及びＲＯＭを通信可能に接続している。
【００８３】
＜ハードウェアとソフトウェアとの協働＞
本発明の方向情報分布推定装置、音源数推定装置、音源方向測定装置、音源分離装置は、コンピュータの記録部に、本発明の各構成部としてを動作させるプログラムを読み込ませ、処理部、入力部、出力部などを動作させることで実現できる。また、コンピュータに読み込ませる方法としては、プログラムをコンピュータ読み取り可能な記録媒体に記録しておき、記録媒体からコンピュータに読み込ませる方法、サーバ等に記録されたプログラムを、電気通信回線等を通じてコンピュータに読み込ませる方法などがある。
【図面の簡単な説明】
【００８４】
【図１】従来の方向情報分布推定装置の機能構成例を示したブロック図。
【図２】Ａは従来の方向情報分布推定装置に入力されるヒストグラムであり、Ｂは従来の方向情報分布推定装置により正規分布をフィッティング処理の結果であり、Ｃはこれら正規分布についての混合正規分布を示す。
【図３】本実施例の方向情報分布推定装置を示した図。
【図４】本実施例の方向情報分布推定装置の処理フローを示した図。
【図５】Ａは本実施例の方向情報分布推定装置に入力されるヒストグラムであり、Ｂは本実施例の方向情報分布推定装置により正規分布をフィッティング処理の結果であり、Ｃはこれら正規分布についての混合正規分布を示す。
【図６】本実施例の音源数測定装置の機能構成例を示したブロック図。
【図７】本実施例の音源方向測定装置の機能構成例を示したブロック図。
【図８】本実施例の音源分離装置の機能構成例を示したブロック図。
【図９】実験２の実験条件を示した図。
【図１０】実験２の結果を示す図。
【図１】
2010145836.tif 000016
【図２】
2010145836.tif 000017
【図３】
2010145836.tif 000018
【図４】
2010145836.tif 000019
【図５】
2010145836.tif 000020
【図６】
2010145836.tif 000021
【図７】
2010145836.tif 000022
【図８】
2010145836.tif 000023
【図９】
2010145836.tif 000024
【図１０】
2010145836.tif 000025