(57)【要約】【課題】音源数が不明であっても、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができるようにする。【解決手段】パラメータ更新部２３により、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源を１個とし、各音源ｋがアクティブになる確率の事前分布を棒折過程でモデル化し、各音源ｋの到来方向のインデックスを示すインジケータの事前分布を離散確率分布でモデル化し、各音源ｋが各インデックスの到来方向となる確率の事前分布をディレクレ分布でモデル化し、音源信号の生成モデルのパラメータの事後分布を近似する変関数のパラメータを変分推論法に基づき繰り返し更新する。時間周波数成分推定部２６は、最終的に更新されたパラメータを用いて、各音源ｋの音源信号の時間周波数成分を推定する。【選択図】図１
2014048399.tif 000002
【特許請求の範囲】
【請求項１】
Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分ｙm,ω,t（ｍはマイクロホン、ｔは時刻、ωは周波数のインデックスである。）を要素にもつ三次元配列ｙ＾を出力する時間周波数解析手段と、
各周波数ωについて、D個（Dは1以上の整数）の音源ｋの各々から各マイクロホンｍへの音源信号の伝達周波数特性ａm,k,ωを要素にもつ三次元配列ａ＾に含まれる各ベクトルａk,ω＾の事後分布のパラメータ（複素正規分布の平均ベクトルｍk,ω＾、共分散行列Γk,ω）、各周波数ωに対する各時刻ｔにおいてアクティブな音源の時間周波数成分ｓω,tを要素にもつ二次元配列ｓ＾に含まれる各要素ｓω,tの事後分布のパラメータ（複素正規分布の平均μω,t、分散σ2ω,t）、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源のインデックスを示すインジケータｚω,tを要素にもつ二次元配列ｚ＾の各要素ｚω,tの事後分布のパラメータ（離散確率分布の各パラメータφk,ω,t）、各音源ｋがアクティブになる確率ｖkを要素にもつベクトルｖ＾の各要素ｖkの事後分布のパラメータ（ベータ分布のパラメータγk,0、γk,1）、各音源ｋの到来方向のインデックスｉを示すインジケータｘkを要素にもつベクトルｘ＾の各要素ｘkの事後分布のパラメータ（離散確率分布の各パラメータψk,i）、各音源ｋが各インデックスｉの到来方向となる確率ρk,iを要素にもつ二次元配列ρ＾の各ベクトルρk＾の事後分布のパラメータ（ディレクレ分布の各パラメータζk,i）の各々の初期値を設定する初期設定手段と、
前記三次元配列ｙ＾が与えられたときの前記三次元配列ａ＾、前記二次元配列ｓ＾、前記二次元配列ｚ＾、前記ベクトルｖ＾、前記ベクトルｘ＾、及び前記二次元配列ρ＾の事後分布ｐ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾｜ｙ＾）と変関数ｑ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、（ｋ、ｍ）の全ての組み合わせの各々における前記平均ベクトルｍk,ω＾及び前記分散Γk,ωと、(ω、ｔ)の全ての組み合わせの各々における前記平均μω,t及び前記分散σ2ω,tと、(ｋ、ω、ｔ)の全ての組み合わせの各々における前記パラメータφk,ω,tと、全てのｋにおける前記パラメータγk,0、γk,1、前記パラメータψk,iと、（ｋ、ｉ）の全ての組み合わせの各々におけるパラメータζk,iとを更新するパラメータ更新手段と、
予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新を繰り返し行う終了判定手段と、
（ｋ、ω、ｔ）の全ての組み合わせの各々について、前記パラメータφk,ω,t及び平均μω,tに基づいて、前記音源ｋの音源信号の時間周波数成分ｓk,ω,tを推定する音源信号推定手段と、
を含む音響信号解析装置。
【請求項２】
前記変関数ｑ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾）を、ｑ（ａ＾）ｑ（ｓ＾）ｑ（）ｑ（ｚ＾）ｑ（ｖ＾）ｑ（ｘ＾）ｑ（ρ＾）で近似し、
前記目的関数を最小化するｑ（ａ＾）を、（ｋ、ω）の全ての組み合わせの各々における前記ベクトルａk,ω＾の複素正規分布を掛け合わせたものとし、
前記目的関数を最小化するｑ（ｓ＾）を、（ω、ｔ）の全ての組み合わせの各々における前記要素ｓω,tの複素正規分布を掛け合わせたものとし、
前記目的関数を最小化するｑ（ｚ＾）を、（ω、ｔ）の全ての組み合わせの各々における前記要素ｚω,tの離散確率分布を掛け合わせたものとし、
前記目的関数を最小化するｑ（ｖ＾）を、全てのｋの各々における前記要素ｖkのベータ分布を掛け合わせたものとし、
前記目的関数を最小化するｑ（ｘ＾）を、全てのｋの各々における前記要素ｘkの離散確率分布を掛け合わせたものとし、
前記目的関数を最小化するｑ（ρ＾）を、全てのｋの各々における前記ベクトルρ＾kのディレクレ分布を掛け合わせたものとした請求項１記載の音響信号解析装置。
【請求項３】
音源数の打ち切り数Dを予め定め、
D＋１以上となるｋ’についてｑ（ｚω,t＝ｋ’）＝０とした請求項２記載の音響信号解析装置。
【請求項４】
前記パラメータ更新手段によって更新された前記（ｋ、ｍ）の全ての組み合わせの各々における前記平均ベクトルｍk,ωのスケールを正規化し、前記音源のインデックスｋについて、音源ｋがアクティブになる周波数ω、時刻ｔにおける前記平均μω,tが降順となるように前記インデックスｋを並び替えるパラメータ調整手段とを更に含み、
前記終了判定手段は、前記予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記パラメータ調整手段による処理を繰り返し行う請求項１又は２記載の音響信号分析装置。
【請求項５】
時間周波数解析手段によって、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分ｙm,ω,t（ｍはマイクロホン、ｔは時刻、ωは周波数のインデックスである。）を要素にもつ三次元配列ｙ＾を出力し、
初期設定手段によって、各周波数ωについて、D個（Dは1以上の整数）の音源ｋの各々から各マイクロホンｍへの音源信号の伝達周波数特性ａm,k,ωを要素にもつ三次元配列ａ＾に含まれる各ベクトルａk,ω＾の事後分布のパラメータ（複素正規分布の平均ベクトルｍk,ω＾、共分散行列Γk,ω）、各周波数ωに対する各時刻ｔにおいてアクティブな音源の時間周波数成分ｓω,tを要素にもつ二次元配列ｓ＾に含まれる各要素ｓω,tの事後分布のパラメータ（複素正規分布の平均μω,t、分散σ2ω,t）、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源のインデックスを示すインジケータｚω,tを要素にもつ二次元配列ｚ＾の各要素ｚω,tの事後分布のパラメータ（離散確率分布の各パラメータφk,ω,t）、各音源ｋがアクティブになる確率ｖkを要素にもつベクトルｖ＾の各要素ｖkの事後分布のパラメータ（ベータ分布のパラメータγk,0、γk,1）、各音源ｋの到来方向のインデックスｉを示すインジケータｘkを要素にもつベクトルｘ＾の各要素ｘkの事後分布のパラメータ（離散確率分布の各パラメータψk,i）、各音源ｋが各インデックスｉの到来方向となる確率ρk,iを要素にもつ二次元配列ρ＾の各ベクトルρk＾の事後分布のパラメータ（ディレクレ分布の各パラメータζk,i）の各々の初期値を設定し、
パラメータ更新手段によって、前記三次元配列ｙ＾が与えられたときの前記三次元配列ａ＾、前記二次元配列ｓ＾、前記二次元配列ｚ＾、前記ベクトルｖ＾、前記ベクトルｘ＾、及び前記二次元配列ρ＾の事後分布ｐ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾｜ｙ＾）と変関数ｑ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、（ｋ、ｍ）の全ての組み合わせの各々における前記平均ベクトルｍk,ω＾及び前記分散Γk,ωと、(ω、ｔ)の全ての組み合わせの各々における前記平均μω,t及び前記分散σ2ω,tと、(ｋ、ω、ｔ)の全ての組み合わせの各々における前記パラメータφk,ω,tと、全てのｋにおける前記パラメータγk,0、γk,1、前記パラメータψkと、（ｋ、ｉ）の全ての組み合わせの各々におけるパラメータζk,iとを更新し、
終了判定手段によって、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新を繰り返し行い、
音源信号推定手段によって、（ｋ、ω、ｔ）の全ての組み合わせの各々について、前記パラメータφk,ω,t及び平均μω,tに基づいて、前記音源ｋの音源信号の時間周波数成分ｓk,ω,tを推定する
音響信号解析方法。
【請求項６】
請求項１〜請求項４の何れか１項に記載の音響信号解析装置の各手段としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号解析装置、方法、及びプログラムに係り、特に、複数のマイクロホンから出力される音響信号の時系列データから、各音源の信号に分離する音響信号解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
音源の成分と音源からマイクロホンまでの伝達特性がともに未知のもとで、マイクロホン入力信号から個々の音源成分を分離抽出する技術をブラインド音源分離(Blind Source Separation;BSS)という。BSSでは観測信号だけから音源信号とその混合過程を推定する必要があるため、通常は音源に関して何らかの仮定を置き、これにより立てられる規準をもとに両未知変数を推定する最適化問題として定式化される。例えば、観測信号数が音源数以上の場合には、音源信号成分が優ガウス分布に従うという仮定のもとで分離フィルタを最尤推定する独立成分分析と呼ぶ方法が有名である。しかし、音源数を仮定しないBSSシステムを実現するためには、マイクロホン数よりも音源数が多い劣決定な問題設定を想定しておく必要があり、独立成分分析をそのまま適用することはできない。劣決定の条件下では、たとえ混合過程が既知であったとしても解が一意に決められないため、音源に関して独立性よりさらに強い仮定が必要となる。
【０００３】
劣決定BSSは難しい不良設定問題であるが、音源数が既知の場合には、時間周波数領域において音声のエネルギーが一部でしか支配的にならないというスパース性と呼ぶ性質を利用した時間周波数領域のBSSアプローチが近年有効なアプローチの一つとして注目されている（非特許文献１）。音声のスパース性とは、音声信号の時間周波数成分が多くの領域でほぼ0となる性質である。このため、複数の音声が同時に発話された状況でも、各時間周波数において音声の時間周波数成分は互いにほとんど重なり合わないと仮定できる場合が多い。この仮定をもとに、目的音声信号の時間周波数成分のみを通過させるような時間周波数マスクをいかにうまく設計するかがこのアプローチにおける問題の焦点となる。また、このアプローチでは、周波数ごとに分離した成分を音源ごとにまとめるためのパーミュテーション整合と呼ばれる問題を解決する必要があり、従来技術では、パーミュテーション整合は周波数ごとの信号分離の後段処理として行われることが多かった。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment,” IEEE Trans. ASLP, vol. 19, no. 3, pp. 516−527, 2010.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
音声信号を対象としたBSSは、ハンズフリーテレビ会議システムや会議録コンテンツの自動作成システムなど、多くの応用が期待されている。例えば会議の場面では参加人数が途中で変化したりドアの開閉音などが突発的に鳴ったりすることがあるように、実環境においては予めあらゆる音源の数を想定しておくことは難しい。従来の多くのBSSアルゴリズムは音源数を仮定して動作するものが多く、仮定した音源数が実際の音源数と異なる場合、高い性能を発揮できない場合がある。
【０００６】
従って、音源数を仮定することなく自律的に音源数を推論しながら動作するBSSシステムの実現が望まれる。音源数が既知の場合には、時間周波数領域において音声のエネルギーが一部でしか支配的にならないというスパース性と呼ぶ性質を利用した時間周波数領域のBSSアプローチが有効であるが、音源数が不明の場合にBSSの問題をいかにして解決するか、が問題となる。
【０００７】
各マイクロホンにおける観測信号は、通常、音源信号の時間遅れを含む畳み込み混合で表されるが、音源からマイクロホンまでのインパルス応答長に対して十分に長い時間窓をもつ時間周波数分解(短時間Fourier変換、ウェーブレット変換など) を用いると、畳み込み混合を近似的に瞬時混合で表すことができる。この観測モデルに基づくBSSは周波数領域BSSと呼ばれ、時間領域の畳み込み混合モデルに基づくBSSに対し、演算量の少ないアルゴリズムを実現できる点や上述の音声のスパース性の仮定を組み込める点などの特長がある一方で、周波数ごとに分離した成分を音源ごとにまとめるためのパーミュテーション整合と呼ばれる問題を扱う必要がある。従来アプローチでは、パーミュテーション整合は周波数ごとの信号分離の後段処理として行われることが多かったが、パーミュテーション整合において手がかりとなる音源の周波数間における関係やスペクトル構造は各周波数における信号分離の精度向上にも寄与しうるはずである。そこで、パーミュテーション整合問題と音源分離をいかに同時解決するか、が問題となる。
【０００８】
本発明は、上記の事情を鑑みてなされたもので、音源数が不明であっても、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる音響信号解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記の目的を達成するために本発明に係る音響信号解析装置は、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分ｙm,ω,t（ｍはマイクロホン、ｔは時刻、ωは周波数のインデックスである。）を要素にもつ三次元配列ｙ＾を出力する時間周波数解析手段と、各周波数ωについて、D個（Dは1以上の整数）の音源ｋの各々から各マイクロホンｍへの音源信号の伝達周波数特性ａm,k,ωを要素にもつ三次元配列ａ＾に含まれる各ベクトルａk,ω＾の事後分布のパラメータ（複素正規分布の平均ベクトルｍk,ω＾、共分散行列Γk,ω）、各周波数ωに対する各時刻ｔにおいてアクティブな音源の時間周波数成分ｓω,tを要素にもつ二次元配列ｓ＾に含まれる各要素ｓω,tの事後分布のパラメータ（複素正規分布の平均μω,t、分散σ2ω,t）、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源のインデックスを示すインジケータｚω,tを要素にもつ二次元配列ｚ＾の各要素ｚω,tの事後分布のパラメータ（離散確率分布の各パラメータφk,ω,t）、各音源ｋがアクティブになる確率ｖkを要素にもつベクトルｖ＾の各要素ｖkの事後分布のパラメータ（ベータ分布のパラメータγk,0、γk,1）、各音源ｋの到来方向のインデックスｉを示すインジケータｘkを要素にもつベクトルｘ＾の各要素ｘkの事後分布のパラメータ（離散確率分布の各パラメータψk,i）、各音源ｋが各インデックスｉの到来方向となる確率ρk,iを要素にもつ二次元配列ρ＾の各ベクトルρk＾の事後分布のパラメータ（ディレクレ分布の各パラメータζk,i）の各々の初期値を設定する初期設定手段と、前記三次元配列ｙ＾が与えられたときの前記三次元配列ａ＾、前記二次元配列ｓ＾、前記二次元配列ｚ＾、前記ベクトルｖ＾、前記ベクトルｘ＾、及び前記二次元配列ρ＾の事後分布ｐ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾｜ｙ＾）と変関数ｑ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、（ｋ、ｍ）の全ての組み合わせの各々における前記平均ベクトルｍk,ω＾及び前記分散Γk,ωと、(ω、ｔ)の全ての組み合わせの各々における前記平均μω,t及び前記分散σ2ω,tと、(ｋ、ω、ｔ)の全ての組み合わせの各々における前記パラメータφk,ω,tと、全てのｋにおける前記パラメータγk,0、γk,1、前記パラメータψk,iと、（ｋ、ｉ）の全ての組み合わせの各々におけるパラメータζk,iとを更新するパラメータ更新手段と、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新を繰り返し行う終了判定手段と、（ｋ、ω、ｔ）の全ての組み合わせの各々について、前記パラメータφk,ω,t及び平均μω,tに基づいて、前記音源ｋの音源信号の時間周波数成分ｓk,ω,tを推定する音源信号推定手段と、を含んで構成されている。
【００１０】
本発明に係る音響信号解析方法は、時間周波数解析手段によって、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分ｙm,ω,t（ｍはマイクロホン、ｔは時刻、ωは周波数のインデックスである。）を要素にもつ三次元配列ｙ＾を出力し、初期設定手段によって、各周波数ωについて、D個（Dは1以上の整数）の音源ｋの各々から各マイクロホンｍへの音源信号の伝達周波数特性ａm,k,ωを要素にもつ三次元配列ａ＾に含まれる各ベクトルａk,ω＾の事後分布のパラメータ（複素正規分布の平均ベクトルｍk,ω＾、共分散行列Γk,ω）、各周波数ωに対する各時刻ｔにおいてアクティブな音源の時間周波数成分ｓω,tを要素にもつ二次元配列ｓ＾に含まれる各要素ｓω,tの事後分布のパラメータ（複素正規分布の平均μω,t、分散σ2ω,t）、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源のインデックスを示すインジケータｚω,tを要素にもつ二次元配列ｚ＾の各要素ｚω,tの事後分布のパラメータ（離散確率分布の各パラメータφk,ω,t）、各音源ｋがアクティブになる確率ｖkを要素にもつベクトルｖ＾の各要素ｖkの事後分布のパラメータ（ベータ分布のパラメータγk,0、γk,1）、各音源ｋの到来方向のインデックスｉを示すインジケータｘkを要素にもつベクトルｘ＾の各要素ｘkの事後分布のパラメータ（離散確率分布の各パラメータψk,i）、各音源ｋが各インデックスｉの到来方向となる確率ρk,iを要素にもつ二次元配列ρ＾の各ベクトルρk＾の事後分布のパラメータ（ディレクレ分布の各パラメータζk,i）の各々の初期値を設定し、パラメータ更新手段によって、前記三次元配列ｙ＾が与えられたときの前記三次元配列ａ＾、前記二次元配列ｓ＾、前記二次元配列ｚ＾、前記ベクトルｖ＾、前記ベクトルｘ＾、及び前記二次元配列ρ＾の事後分布ｐ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾｜ｙ＾）と変関数ｑ（ａ＾、ｓ＾、ｚ＾、ｖ＾、ｘ＾、ρ＾）との間の差異を表すダイバージェンスを表す関数を目的関数として、変分推論法に基づき前記目的関数を最小化するように、（ｋ、ｍ）の全ての組み合わせの各々における前記平均ベクトルｍk,ω＾及び前記分散Γk,ωと、(ω、ｔ)の全ての組み合わせの各々における前記平均μω,t及び前記分散σ2ω,tと、(ｋ、ω、ｔ)の全ての組み合わせの各々における前記パラメータφk,ω,tと、全てのｋにおける前記パラメータγk,0、γk,1、前記パラメータψkと、（ｋ、ｉ）の全ての組み合わせの各々におけるパラメータζk,iとを更新し、終了判定手段によって、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新を繰り返し行い、音源信号推定手段によって、（ｋ、ω、ｔ）の全ての組み合わせの各々について、前記パラメータφk,ω,t及び平均μω,tに基づいて、前記音源ｋの音源信号の時間周波数成分ｓk,ω,tを推定する。
【００１１】
本発明に係るプログラムは、上記の音響信号解析装置の各手段としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１２】
以上説明したように、本発明の音響信号解析装置、方法、及びプログラムによれば、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源を１個とし、各音源ｋがアクティブになる確率をベータ分布で生成し、各音源ｋの到来方向のインデックスを示すインジケータを離散確率分布で生成し、各音源ｋが各インデックスの到来方向となる確率をディレクレ分布で生成する音源信号の生成モデルのパラメータの事後分布に基づく目的関数を最小化するように、変分推論法に基づき各パラメータを更新して、各音源ｋの音源信号の時間周波数成分の事後分布を推定することにより、音源数が不明であっても、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる、という効果が得られる。
【図面の簡単な説明】
【００１３】
【図１】本発明の実施の形態に係る音響信号解析装置の構成を示す概略図である。
【図２】本発明の実施の形態に係る音響信号解析装置における音響信号解析処理ルーチンの内容を示すフローチャートである。
【図３】従来法と提案法について得られたＳＤＲとＫ又はＫ*との関係を示す図である。
【図４】観測信号の合成に利用した室内インパルス応答のチャンネル間位相差を示すグラフである。
【図５】推定した伝達周波数特性より算出されるチャンネル間位相差を示すグラフである。
【発明を実施するための形態】
【００１４】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明の目的は、(1)音源数の推論、(2)音声のスパース性を仮定した劣決定周波数領域BSS、(3)パーミュテーション整合、の問題を一挙に解決することである。本発明で提案する手法は、観測信号が与えられた下で、各時間周波数ビンにおいてどの音源がアクティブであるらしいかの事後確率、各音源がどれだけアクティブになりやすいか(どの程度の音源数が存在しているらしいか)の事後確率、各音源の時間周波数成分の事後分布、各音源のステアリングベクトルの事後分布、各音源がどの方向から到来したらしいかの事後確率を、変分ベイズ法に基づいて推定することによりこれを実現する。
【００１５】
＜発明の原理＞
まず、本発明の原理について説明する。
【００１６】
＜観測モデル＞
まず、観測モデルについて説明する。K個(Kは1以上の整数)の信号源から到来する音源信号をM個(Mは2以上の整数)のマイクロホンで観測する場合を考え、m番目のマイクロホンで観測される信号の時間周波数成分をym(ω,t)、k番目の音源信号の時間周波数成分をsk(ω,t)とし、ｙ＾(ω,t)=(y1(ω,t),...,yM(ω,t))T∈CM，ｓ＾(ω,t)=(s1(ω,t),...,sK(ω,t))T∈CKとする。ただし、1≦ω≦Ω,1≦t≦Tはそれぞれ周波数および時刻に対応するインデックスである。先に述べたとおり、時間周波数領域において観測信号ｙ＾(ω, t)は近似的に以下の（１）式のように、ｓ＾(ω,t)の瞬時混合の形で表すことができる。なお、記号に付された「＾」は、当該記号が行列または多次元配列またはベクトルであることを表わしている。
【００１７】
【数１】
2014048399.tif 000003
【００１８】
信号源kからマイクロホンmまでの伝達周波数特性am,k(ω)を要素にした行列a＾(ω)=(am,k(ω))M×K=(a1(ω),...,aK(ω))∈CM×Kを混合行列と呼び、以下ではこれを時不変と仮定する。n＾(ω,t)は、多数の方向から到来する背景雑音や、フレーム長を超える残響成分など、時不変な伝達特性として表現できない成分を表す。ここで、音声のスパース性が仮定できる場合、各時間周波数ビン(ω,t)においてアクティブ(支配的)となる音源インデックスをzω,t∈{1,...,K}で表すことにすると、上記（１）式は、以下の（２）式のように書き直せる。
【００１９】
【数２】
2014048399.tif 000004
【００２０】
この観測モデルでは、zω,t番目の音源以外の成分はすべて0であると仮定しているので、各時間周波数ビンにおいて音源成分を表す変数は一つだけで十分である。このため上記（２）式ではsk(ω,t)のインデックスkを省いている。すなわち、s(ω,t)は特定の音源の成分ではなく、各時間周波数ビンにおいてアクティブな１個の音源の成分を表す変数である。紙面のスペースの節約のため、以後ωとtを下付き添え字で表記することにする。
【００２１】
＜生成モデル＞
＜観測信号の生成プロセス＞
上述した観測モデルをもとに、観測信号が生成されるプロセスを生成モデルにより記述する。
【００２２】
まず、雑音成分nω,t＾は、平均が0＾、共分散がΣ(n)ωの複素正規分布に従うと仮定すると、もしa1:K,ω＾={a1,ω,...,aK,ω},sω,t，および、各時間周波数ビンでどの音源がアクティブであるか、すなわちzω,tが既知であれば、上記（２）式より、yω,tは、以下の（３）式に従って生成される。
【００２３】
【数３】
2014048399.tif 000005
【００２４】
＜無限音源数混合モデル＞
上記ではzω,tが既知の下での観測信号の生成プロセスを仮定したが、通常は各時間周波数ビンでどの音源がアクティブであるかに関する情報は観測することができない。そこで、本実施の形態では、アクティブな音源インデックスを示すインジケータzω,tを潜在変数と見なし、その生成プロセスをモデル化する。まず、音源数がKの場合、zω,tは、音源インデックスの集合{1,...,K}からいずれかのインデックスが、ある離散分布に従って選ばれる、以下の（４）式に示すプロセスによって生成されると仮定する。
【００２５】
【数４】
2014048399.tif 000006
【００２６】
ただし、π＾=(π,...,πK)は離散分布における各インデックスの選ばれやすさを意味する確率値であり、ΣKk=1πk=1とする。さらに、これらの確率値は、以下の（５）式に示すDirichlet分布よって生成されると仮定する。
【００２７】
【数５】
2014048399.tif 000007
【００２８】
ここまでは有限個の音源数Kを想定していたが、ここで、K→∞の極限をとると、上記（３）式、（４）式、（５）式は、Dirichlet過程混合モデルとなり、以下の（６）式〜（９）式で表すことができる。
【００２９】
【数６】
2014048399.tif 000008
【００３０】
GEM(α0)は、Dirichlet過程（参考文献（T. S. Ferguson, “A Bayesian analysis of some nonparametric problems,” Annals of Statistics, vol.1, no. 2, pp. 209−230, 1973.）を参照）により生成されるπ1,π2,...の一つの具体的な構成方法(棒折過程と呼ぶ（参考文献（J. Sethuraman, “A constructive definition of Dirichlet priors,” Statistica Sinica, vol. 4, pp. 639−650, 1994.）を参照）)であり、以下の（１０）式、（１１）式に従って与えられる。
【００３１】
【数７】
2014048399.tif 000009
【００３２】
以上のプロセスで生成されるπ1,π2,...は、平均的な意味で、kが大きいほどπkが指数的に小さくなるという傾向を持つため、大きいkに対応した音源ほどアクティブになる確率が低くなることを意味する。よって、観測信号からパラメータを推論する際、必要最小限の音源インデックス数の混合モデルで観測信号を説明しようとする効果がもたらされる。以上のyω,tの生成モデルを「無限音源数混合モデル」と呼ぶ。
【００３３】
＜混合DOAモデル＞
上述したモデルにパーミュテーション整合機能を組み込むため、上記（９）式の生成プロセスを以下でモデル化する。ここまで各音源の伝達周波数特性ak,ωを周波数ωごとの独立な変数であるかのように扱っていたが、もし各音源が単一方向から平面波到来すると仮定できるならば、例えばマイクロホン数が2の場合、伝達周波数特性の各ω間の関係は、到来方向(Direction-of-Arrival;DOA)θの関数として以下の（１２）式のように陽に表される。
【００３４】
【数８】
2014048399.tif 000010
【００３５】
ただし、0≦θ＜2π、hωは周波数インデックスωに対応する角周波数(rad/s)，Dをマイクロホンの間隔(m)、Cを音速(m/s)とする。実際には残響や周波数領域の瞬時混合近似などの影響により、ak,ωは上記の理論式からは逸脱することが予想される。そこで、到来方向θkが既知のとき、ak,ωは、
【００３６】
【数９】
2014048399.tif 000011
【００３７】
を中心とした複素正規分布より生成されたものと仮定する。しかし当然ながら到来方向θkは実際には観測することができないため、これを潜在変数と見なすことにすると、ak,ωの生成モデルは、DOAを潜在変数とした混合モデルとなる。これを、上述した無限音源数混合モデルに組み込み、観測信号が与えられた下で全体の生成モデルのパラメータ推論を行うことができれば、音源分離とパーミュテーション整合を同時解決できる可能性がある。
【００３８】
まず、Θ1,...,ΘI(すべて定数)からなるI個のDOA候補値の集合を用意する。例えば、180度をI等分した角度Θi=(i−1)π/I,(i=1,...,I)の集合としよう。各音源のDOAがこのDOA候補値の中から一つ選ばれて決定される、というプロセスを仮定するなら、θkが生成されるプロセスは以下の（１３）式、（１４）式のように記述できる。
【００３９】
【数１０】
2014048399.tif 000012
【００４０】
ただし、ρk＾=(ρk,1,...,ρk,I)である。xk∈{1,...,I}はk番目の音源にどのDOA候補値が割り当てられるかを表すインジケータ変数であり、上記（１３）式はこれが離散分布(各確率値がρk,1,...,ρk,I)から生成されることを意味している。このプロセスにより各音源のDOAが決定され、伝達周波数特性ak(ω)は、以下の（１５）により生成される。
【００４１】
【数１１】
2014048399.tif 000013
【００４２】
また、ここで、ρk＾の事前分布として、以下の（１６）式で表されるDirichlet分布を仮定する。
【００４３】
【数１２】
2014048399.tif 000014
【００４４】
以上の、混合モデルに基づくak(ω)の生成モデルを「混合DOAモデル」と呼ぶ。
【００４５】
＜変分推論アルゴリズム＞
観測信号ｙ＾=ｙ＾1:Ω,1:Tが与えられたもとで、以上の生成モデルのパラメータa＾=a1:∞,1:Ω＾,ｓ＾=s1:Ω,1:T,ｚ＾=z1:Ω,1:T,ｖ＾=v1:∞,ｘ＾=x1:∞,ρ＾=ρ＾1:Kの事後分布p(a＾,ｓ＾,ｚ＾,ｖ＾,ｘ＾,ρ＾｜ｙ＾)を求めたい。この事後分布を解析的に得ることは難しいが、変分推論法に基づき近似分布を反復計算により得ることができる。以下、簡単のため、Σ(n)1:Ω,Σ(a)1:Ω,α0,β0は実験的に定める定数とする。変分推論は、事後分布p(a＾,ｓ＾,ｚ＾,ｖ＾,ｘ＾,ρ＾|ｙ＾)と、以下の（１７）式を満たす非負の変関数q(a＾,ｓ＾,ｚ＾,ｖ＾,ｘ＾,ρ＾)との間の以下の（１８）式に示すKullback-Leiblerダイバージェンスをqに関して最小化することが目的となる。
【００４６】
【数１３】
2014048399.tif 000015
【００４７】
ただし〈f(x)〉q(x)は∫q(x)f(x)dxを表す。そしてqに関して、以下の（１９）式のように近似できると仮定し、q(a＾),q(ｓ＾),q(ｚ＾),q(ｖ＾),q(ｘ＾),q(ρ＾)について反復的にF[q] を最小化することでp(a＾,ｓ＾,ｚ＾,ｖ＾,ｘ＾,ρ＾|ｙ＾)の近似分布を得ようというのが変分推論法の基本的な考え方である。
【００４８】
【数１４】
2014048399.tif 000016
【００４９】
また、q(z)に関して、以下の（２０）式に示すように、打ち切り近似を行う。
【００５０】
【数１５】
2014048399.tif 000017
【００５１】
この近似は、モデルの複雑度(音源数)を固定した、ということではなく、qの関数空間をある領域に限定した、ということを意味している。よって、本来推定したいp(a＾,ｓ＾,ｚ＾,ｖ＾,ｘ＾,ρ＾|ｙ＾)をできるだけ良くqでフィッティングしたければ、Dは大きければ大きいほど良い、ということになる。導出は省略するが、上記（１８）式を上記（１７）式の拘束の下で最小化する各qは、解析的に以下の（２１）式〜（２６）式に示す形として求まる。
【００５２】
【数１６】
2014048399.tif 000018
【００５３】
また、mk,ω＾, Γk,ω, μω,t, σ2ω,t,φω,t＾ = (φ1,ω,t, . . . , φD,ω,t), γk,0, γk,1,ψk＾= (ψk,1, . . . , ψk,I ), ζk,1, . . . , ζk,Iは、以下の（２７）式〜（３７）式で表される。ただし、
【００５４】
【数１７】
2014048399.tif 000019
とする。
【００５５】
【数１８】
2014048399.tif 000020
【００５６】
ただし、Ψはディガンマ関数である。また、記号に付された「*」は複素共役を表す。上記γk,0とγk,1の更新式（（３６）式、（３７）式）では第二項にv0が乗じられているが、これは事前分布をv0乗することに相当し、音声データが長いなどの原因でデータ情報が多くなった場合に事前分布の影響が小さくなりすぎるのを防ぐための措置である。後述する実験ではv0=0.7×Ω×Tとした。
【００５７】
＜システム構成＞
次に、Ｍ個（Ｍは２以上の整数）のマイクロホンから得られた音響信号を解析して、未知のＫ個の音源信号に分離する音響信号解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００５８】
図１に示すように、本発明の実施の形態に係る音響信号解析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００５９】
音響信号解析装置１００は、入力部１０と、演算部２０と、記憶部３０と、出力部４０とを備えている。
【００６０】
入力部１０により、Ｍ個のマイクロホンから出力された音響信号（多チャンネル信号）の時系列データが入力される。記憶部３０は、入力部１０により入力された音響信号の時系列データを記憶する。また、記憶部３０は、後述する各処理での結果を記憶すると共に、定数のパラメータΣ(n)ω＾、Σ(a)ω＾、α0、β0を記憶している。
【００６１】
例えば、マイクロホン数Ｍを2とする。定数のパラメータを、それぞれ以下の（３８）式〜（４１）式に示すように設定する。
【００６２】
【数１９】
2014048399.tif 000021
【００６３】
ただし、I＾はM×Mの単位行列である。また、角度の分割数をI=180とする。v0=0.7×Ω×Tとする。Σ(n)ω＾は
【００６４】
【数２０】
2014048399.tif 000022
(スパース音源モデル)のyω,t(観測信号)からのずれの許容量を表し、Σ(a)ω＾は、ak,ωの
【００６５】
【数２１】
2014048399.tif 000023
(音源k の到来方向が既知のときその方向に対応するステアリングベクトル) からのずれの許容量を表す。また、α0は、大きいインデックスの音源がアクティブになりにくくなる度合いを表し、β0は到来方向xkの分布ρk,1,...,ρk,Iがスパースになる度合いを表す。
【００６６】
演算部２０は、時間周波数解析部２１と、初期設定部２２と、パラメータ更新部２３と、パラメータ調整部２４と、終了判定部２５と、時間周波数成分推定部２６と、信号変換部２７とを備えている。また、パラメータ更新部２３は、アクティブ音源事後分布更新部２３１と、時間周波数成分事後分布更新部２３２と、音源インジケータ事後確率更新部２３３と、音源到来方向事後確率更新部２３４と、ステアリングベクトル事後確率更新部２３５と、音源方向事後確率更新部２３６と、を備えている。
【００６７】
時間周波数解析部２１は、各マイクロホンの時系列信号としての観測された音響信号を入力として、時間周波数成分（観測時間周波数成分）ｙm,ω,t（ｍ＝１，・・・，Ｍ、ω＝１,・・・,Ω,ｔ＝１,・・・,Ｔは、それぞれマイクロホン、周波数、時刻に対応するインデックスを示す。）を各（ｍ，ω，ｔ）の要素にもつ三次元配列ｙ＾を計算する。また、計算した時間周波数成分ｙm,ω,tを、記憶部３０に記憶しておく。より詳細には、時間周波数解析部２１は、各マイクロホンｍについて、当該マイクロホンの音響信号の時系列データを入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより時間周波数成分ｙm,ω,tを計算し、時間周波数成分ｙm,ω,tを格納した行列（振幅スペクトログラム）ｙ＾＝（ｙm,ω,t）M×Ω×Ｔを出力する。なお、時間周波数成分ｙm,ω,tは、ウェーブレット変換を用いて計算してもよい。
【００６８】
初期設定部２２は、Ω、ｌi,ω＾および各変分事後分布qのパラメータ(以後、変分パラメータと称する)の初期値を以下の通りに設定し、記憶部３０に格納する。
【００６９】
まず、DOAが位相に現れる程度に高く、かつエイリアシングが起こらない範囲の周波数帯hΩ0=約300Hz〜約300+120Hzから反復計算を始めるように、Ω＝Ω0を設定する。
【００７０】
また、ｌi,ω＾について、以下の（４２）式に従って、（ｉ，ω）毎に初期値を設定する。
【００７１】
【数２２】
2014048399.tif 000024
【００７２】
また、アクティブな音源信号の平均μω,tについては、各マイクロホンｍで取得した信号の時間周波数成分の平均Σmym,ω,t/Mを初期値として（ｍ，ｔ）毎に設定する。
【００７３】
また、φk,ω,t(音源kが時間周波数(ω,t)においてアクティブである確率を表すパラメータ)を、音源インデックスに反比例するように(大きいインデックスほど低い確率値となるように)、なおかつ、Σkφk,ω,t=1となるように、以下の（４３）式に従って（ｋ、ｍ、ｔ）毎に与える。
【００７４】
【数２３】
2014048399.tif 000025
【００７５】
また、その他の変分パラメータmk,ω＾, Γk,ω, σ2ω,t,γk,0, γk,1,ψk＾= (ψk,1, . . . , ψk,I ), ζk,1, . . . , ζk,Iについては、初期値として、乱数を用いて適当な値を設定すればよい。
【００７６】
時間周波数成分事後分布更新部２３２は、（ω、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているσ2ω,t、φk,ω,t、mk,ω＾, Σ(n)ω＾、ｙm,ω,t、Γk,ωに基づいて、以下の（４４）式、（４５）式に従って、アクティブな音源の時間周波数成分の事後分布のパラメータμω,t、σ2ω,tを更新し、記憶部３０に格納する。
【００７７】
【数２４】
2014048399.tif 000026
【００７８】
ステアリングベクトル事後確率更新部２３５は、（ｋ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているΓk,ω,σ2ω,t、φk,ω,t、mk,ω＾, Σ(n)ω＾、Σ(a)ω＾、μω,t、ｙm,ω,t、ψk、ｌi,ω＾に基づいて、以下の（４６）式、（４７）式に従って、各音源のステアリングベクトルの事後確率のパラメータmk,ω＾、Γk,ω,を更新し、記憶部３０に格納する。
【００７９】
【数２５】
2014048399.tif 000027
【００８０】
音源インジケータ事後確率更新部２３３は、（ｋ、ω、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているΓk,ω,σ2ω,t、φk,ω,t、mk,ω＾, Σ(n)ω＾、μω,t、ｙm,ω,t、γk,0, γk,1に基づいて、以下の（４８）式に従って、アクティブな音源のインジケータ（各時間周波数ビンにおいてどの音源がアクティブであるか）の事後確率のパラメータφk,ω,tを更新し、記憶部３０に格納する。
【００８１】
【数２６】
2014048399.tif 000028
【００８２】
音源到来方向事後確率更新部２３４は、（ｋ、ｉ）の全ての組み合わせの各々について、記憶部３０に記憶されているζk,i、Σ(a)ω＾、ｌi,ω＾、に基づいて、以下の（４９）式に従って、各音源がどの方向から到来したらしいかの事後確率のパラメータψk,iを更新し、記憶部３０に格納する。
【００８３】
【数２７】
2014048399.tif 000029
【００８４】
音源方向事後確率更新部２３６は、（ｋ、ｉ）の全ての組み合わせの各々について、記憶部３０に記憶されているψk,i、β0に基づいて、以下の（５０）式に従って、どの方向に音源があるらしいかの事後確率のパラメータζk,iを更新し、記憶部３０に格納する。
【００８５】
【数２８】
2014048399.tif 000030
【００８６】
アクティブ音源事後分布更新部２３１は、各ｋについて、記憶部３０に記憶されているφk,ω,t、α0、ｖ0に基づいて、以下の（５１）式、（５２）式に従って、各音源がどれだけアクティブになりやすいかの事後分布のパラメータγk,0, γk,1を更新し、記憶部３０に格納する。
【００８７】
【数２９】
2014048399.tif 000031
【００８８】
パラメータ調整部２４は、以下の調整を行う。
【００８９】
まず、以下の（５３）式に従って、（ｋ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているmk,ω＾のスケールを正規化することによって、ステアリングベクトルのスケールを正規化する。
【００９０】
【数３０】
2014048399.tif 000032
【００９１】
ただし、[・]iはベクトルのi番目の要素を表す。
【００９２】
また、反復計算に用いる周波数帯域を、以下の（５４）式に従って、ΔΩ=約120Hzに相当する周波数ビン数分だけ広げる。
【００９３】
【数３１】
2014048399.tif 000033
【００９４】
また、音源インデックスｋについて、音源ｋがアクティブになる周波数ω、時刻ｔにおける時間周波数成分ｓω,tの平均μω,tが降順となるように音源インデックスｋを並び替える。
【００９５】
また、例えば反復回数20回目以降では、更に、各ω毎に、ステアリングベクトルの共分散行列Σ(a)ωを（１／反復回数）倍して小さくする。
【００９６】
終了判定部２５は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、パラメータ更新部２３及びパラメータ調整部２４の各処理を繰り返す。終了判定部２５は、終了条件を満足したと判定した場合には、時間周波数成分推定部２６による処理に移行する。
【００９７】
終了条件としては、反復計算回数が予め定めた回数に達したことを用いればよい。なお、一回のパラメータ更新によるパラメータの変化率がほぼ1になったと見なせたことを、終了条件として用いてもよい。また、一回のパラメータ更新による目的関数F[q]の値の変化率がほぼ1になったと見なせたことを、終了条件として用いてもよい。
【００９８】
時間周波数成分推定部２６は、（ｋ、ω、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているφk,ω,t、μω,tに基づいて、以下の（５５）式に従って、各音源ｋの時間周波数成分の推定値＾sk,ω,tを計算する。
【００９９】
【数３２】
2014048399.tif 000034
【０１００】
信号変換部２７は、時間周波数成分推定部２６によって計算された各音源ｋの時間周波数成分の推定値＾sk,ω,tを、各音源の音源信号に変換し、出力部４０により、音源信号を出力する。
【０１０１】
＜音響信号解析装置の作用＞
次に、本実施の形態に係る音響信号解析装置１００の作用について説明する。まず、解析対象の信号として各マイクロホンからの音響信号の時系列データが音響信号解析装置１００に入力され、記憶部３０に格納される。そして、音響信号解析装置において、図２に示す音響信号解析処理ルーチンが実行される。
【０１０２】
まず、ステップＳ１０１において、記憶部３０から、マイクロホン毎に、各時刻ｔのフレーム内の音響信号を読み込み、当該音響信号に対して、短時間フーリエ変換を用いた時間周波数分析を行った結果から、観測時間周波数成分ｙm,ω,tを各（ｍ，ω，ｔ）の要素にもつ三次元配列ｙ＾を生成して、記憶部３０に記憶する。
【０１０３】
そして、ステップＳ１０２において、Ω、ｌi,ω＾および各変分事後分布qのパラメータμω,t、φk,ω,t、mk,ω＾, Γk,ω, σ2ω,t,γk,0, γk,1, ψk＾= (ψk,1, . . . , ψk,I ), ζk,1, . . . , ζk,Iの初期値を設定し、記憶部３０に格納する。
【０１０４】
次のステップＳ１０３では、上記ステップＳ１０２で設定され、またはステップＳ１０３、Ｓ１０４、Ｓ１０５で更新されたパラメータσ2ω,t、φk,ω,t、mk,ω＾、Γk,ωと、上記ステップＳ１０１で得られたｙm,ω,tとに基づいて、上記（４４）式、（４５）式に従って、アクティブな音源の時間周波数成分の事後分布のパラメータμω,t、σ2ω,tを更新し、記憶部３０に格納する。
【０１０５】
そして、ステップＳ１０４では、上記ステップＳ１０２で設定され、またはステップＳ１０３、Ｓ１０４、Ｓ１０５、Ｓ１０６で更新されたパラメータΓk,ω,σ2ω,t、φk,ω,t、mk,ω＾、μω,t、ψk、ｌi,ω＾と、上記ステップＳ１０１で得られたｙm,ω,tとに基づいて、上記（４６）式、（４７）式に従って、各音源のステアリングベクトルの事後確率のパラメータmk,ω＾、Γk,ω,を更新し、記憶部３０に格納する。
【０１０６】
ステップＳ１０５では、上記ステップＳ１０２で設定され、又はステップＳ１０３、Ｓ１０４、Ｓ１０５、Ｓ１０８で更新されたパラメータΓk,ω,σ2ω,t、φk,ω,t、mk,ω＾、μω,t、γk,0, γk,1と、上記ステップＳ１０１で得られたｙm,ω,tとに基づいて、上記（４８）式に従って、アクティブな音源のインジケータ（各時間周波数ビンにおいてどの音源がアクティブであるか）の事後確率のパラメータφk,ω,tを更新し、記憶部３０に格納する。
【０１０７】
次のステップＳ１０６では、上記ステップＳ１０２で設定され、又はステップＳ１０７で更新されたパラメータζk,i、ｌi,ω＾に基づいて、上記（４９）式に従って、各音源がどの方向から到来したらしいかの事後確率のパラメータψk,iを更新し、記憶部３０に格納する。
【０１０８】
次のステップＳ１０７では、上記ステップＳ１０２で設定され、又はステップＳ１０６で更新されたパラメータψk,iに基づいて、上記（５０）式に従って、どの方向に音源があるらしいかの事後確率のパラメータζk,iを更新し、記憶部３０に格納する。
【０１０９】
ステップＳ１０８では、上記ステップＳ１０２で設定され、又はステップＳ１０５で更新されたパラメータφk,ω,tに基づいて、上記（５１）式、（５２）式に従って、各音源がどれだけアクティブになりやすいかの事後分布のパラメータγk,0, γk,1を更新し、記憶部３０に格納する。
【０１１０】
そして、ステップＳ１０９において、上記ステップ１０４で更新されたパラメータmk,ω＾を正規化し、パラメータΩを更新し、音源インデックスｋの並び替えを行う。また、反復回数20回目以降では、パラメータΣ(a)ωを調整する。
【０１１１】
次のステップＳ１１０では、終了条件として、反復回数が、予め定めた回数に到達したか否かを判定し、反復回数が予め定めた回数に到達していない場合には、終了条件を満足していないと判断して、上記ステップ１０３へ戻り、上記ステップ１０３〜ステップ１０９の処理を繰り返す。一方、反復が予め定めた回数に到達した場合には、終了条件を満足したと判断し、ステップＳ１１１で、上記ステップＳ１０３、Ｓ１０５で最終的に更新されたφk,ω,t、μω,tに基づいて、上記（５５）式に従って、各音源ｋの時間周波数成分の推定値＾sk,ω,tを計算する。
【０１１２】
そして、ステップＳ１１２において、上記ステップＳ１１１で計算された各音源ｋの時間周波数成分の推定値＾sk,ω,tに基づいて、各音源の音源信号を算出し、出力部４０により出力して、音響信号解析処理ルーチンを終了する。
【０１１３】
＜実験結果＞
【０１１４】
次に、本実施の形態に係る手法（以下、提案法と称する）の有用性を示すため、音源分離性能の検証実験を行った結果について説明する。3人の話者(女性2人、男性1 人)の音声信号（参考文献（A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri, H. Kuwabara, and K. Shikano, “ATR Japanese speech database as a tool of speech recognition and synthesis,” Trans. Speech Communication, pp. 357−363, 1990.）を参照）に、室内インパルス応答(残響時間は0ms)（参考文献（S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada, “Acoustical sound database in real environments for sound scene understanding and hands-free speech recognition,” in Proc. LREC ’00, 2000, pp. 965−968.）を参照）を畳み込み加算することで人工的に混合したものを観測信号とした。標本化周波数は16kHzとした。観測信号の時間周波数成分は、短時間Fourier変換(フレーム長は64ms、フレームシフトは16ms)により算出した。Σ(n)ωとΣ(a)ωはそれぞれ＾I,10−1.5×＾Iとした。また、角度の分割数はI=180とした。上記の反復アルゴリズムの実行後、アクティブな音源成分の推定値μω,tに、音源k が各時間周波数ビンでどれだけアクティブらしいかを表す確率値φk,ω,tを乗じたものを、音源kの推定時間周波数成分とした。また、今回の実験では、空間エイリアシングによりζが局所解に陥ってしまう可能性を考慮し、反復計算の初期段階では空間エイリアシングが起こらない低い帯域の観測情報のみを用いてアルゴリズムを実行し、反復回数の増加に従って徐々にその帯域を高帯域に広げていく方法をとった。
【０１１５】
音源分離性能の評価基準として、Signal-to-Distortion Ratio(SDR)を採用した（参考文献（E. Vincent, R. Gribonval, and C. F´evotte, “Performance measurement in blind audio source separation,” IEEE Trans. ASLP, pp. 1462−1469, 2006.）を参照）。まず、提案法の第一のポイントである、必要な音源数を適応的に推論する効果を検証するため、打ち切りレベルDによって分離性能がどう影響するかを確認した。また、従来方法では音源数を仮定する必要があったが、仮定する音源数Kによって従来方法による分離性能がどう変化するかも併せて確認した。その結果を図３に示す。従来方法では、実際の音源数と仮定する音源数が一致する場合は分離性能が高かったが、一致しない場合は著しく低かったのに対して、提案法では、先に述べたとおり、打ち切りレベルDが大きければ大きいほどqが真の事後分布の良い近似が得られるため、Dが大きいほど分離性能が向上することが確認できた。このことは、Dをあらかじめ大きく設定してさえいれば音源数が分かっていなかったとしても高い分離性能が得られる、ということを意味し、提案法の効果を示している。
【０１１６】
次に、混合DOAモデルの効果を見るため、音源の到来方向を正しく推定できているらしいかを確認した。図４は、観測信号の合成に利用した室内インパルス応答のチャンネル間位相差arg([mk,ω＾]2/[mk,ω＾]1)を音源ごとに異なる色でプロットしたものであり、図５は、推定した伝達周波数特性mk,ω＾より算出されるチャンネル間位相差arg([mk,ω]2/[mk,ω]1)を音源ごとに異なる色でプロットしたものである。ただし、[・]iはベクトルのi番目の要素を表す。空間エイリアシングがあっても各音源の到来方向が概ね正しく推定できていることが分かる。
【０１１７】
以上説明したように、本発明の実施の形態に係る音響信号解析装置によれば、各周波数ωに対する各時刻ｔにおいてアクティブとなる音源を１個とし、各音源ｋがアクティブになる確率ｖkの事前分布を棒折過程でモデル化し、各音源ｋの到来方向のインデックスを示すインジケータｘkの事前分布を離散確率分布でモデル化し、各音源ｋが各インデックスの到来方向となる確率ρk＾の事前分布をディレクレ分布でモデル化し、音源信号の生成モデルのパラメータの事後分布を近似する変関数ｑのパラメータを変分推論法に基づき繰り返し更新して、各音源ｋの音源信号の時間周波数成分の事後分布を推定することにより、音源数が不明であっても、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる。
【０１１８】
また、複数のマイクロホンで取得した音響信号から混在する個々の音源信号を分離することができるため、ハンズフリーテレビ会議システムや会議録コンテンツの自動作成システムなどの応用が期待される。
【０１１９】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１２０】
例えば、上述の音響信号解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１２１】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１２２】
１０   入力部
２０   演算部
２１   時間周波数解析部
２２   初期設定部
２３   パラメータ更新部
２４   パラメータ調整部
２５   終了判定部
２６   時間周波数成分推定部
２７   信号変換部
３０   記憶部
４０   出力部
１００ 音響信号解析装置
２３１ アクティブ音源事後分布更新部
２３２ 時間周波数成分事後分布更新部
２３３ 音源インジケータ事後確率更新部
２３４ 音源到来方向事後確率更新部
２３５ ステアリングベクトル事後確率更新部
２３６ 音源方向事後確率更新部
【図１】
2014048399.tif 000035
【図２】
2014048399.tif 000036
【図３】
2014048399.tif 000037
【図４】
2014048399.tif 000038
【図５】
2014048399.tif 000039