(57)【要約】【課題】複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができるようにする。【解決手段】確率更新部２４１、時変ゲイン更新部２４２、及び分離行列更新部２４３により、各周波数に対する各時間フレームにおいて調波構造を有するパワースペクトルテンプレートの各々が音源毎に選択される確率を要素にもつ三次元配列Π＾を用いて、観測時間周波数成分の三次元配列Ｘが与えられたときの三次元配列Π＾、各音源のパワーを要素にもつ二次元配列Σ＾、及び各周波数ωの分離行列Ｗω＾の事後確率を表す目的関数を最大化するように、三次元配列Π＾、二次元配列Σ＾、及び各周波数ωの分離行列Ｗω＾を更新することを繰り返す。【選択図】図１
2014048398.tif 000002
【特許請求の範囲】
【請求項１】
Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘmτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘを出力する時間周波数解析手段と、
Ｍ個の音源ｋの時間周波数成分Ｙkτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ(l)ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π(l)kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ2kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘτω（＝（Ｘ1τω，・・・，ＸMτω））に作用させて音源信号の時間周波数成分Ｙτω（＝（Ｙ1τω，・・・，ＹKτω））を得るための分離行列Ｗω＾の各々の初期値を設定するパラメータ初期値設定手段と、
（k、τ、ω、ｚkτ）の全ての組み合わせにおける、ｚkτとσ2kγが与えられたときの、λ(ｚkτ)ω・σ2kγを分散とするガウス分布で表されるＹkτωの確率密度関数、及びπkτが与えられたときのｚkτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π(l)kτの事前確率、及び各周波数ωに対する前記分離行列Ｗω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾の事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾を更新するパラメータ更新手段と、
各周波数ωの前記分離行列Ｗω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新する音源信号推定値更新手段と、
予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う終了判定手段と、
を含む音響信号解析装置。
【請求項２】
前記目的関数を、（ｋ，τ、ω）の全ての組み合わせについて観測時間周波数成分Ｘmτωが各インデックスｌのパワースペクトルテンプレートλ(l)ωに帰属する確率を示す補助変数γ(l)kτωを用いた補助関数とし、
前記パラメータ更新手段は、
前記三次元配列Π＾、二次元配列Σ＾、及び前記複数のパワースペクトルテンプレートλ(l)ωに基づいて、（ｋ，τ、ω、ｌ）の全ての組み合わせの各々について前記補助変数γ(l)kτωを更新する補助変数更新手段と、
前記補助変数γ(l)kτωに基づいて、前記三次元配列Π＾を更新する確率更新手段と、
前記補助変数γ(l)kτω及び前記複数のパワースペクトルテンプレートλ(l)ωに基づいて、前記二次元配列Σを更新するパワー更新手段と、
前記補助変数γ(l)kτω、前記二次元配列Σ＾、前記三次元配列Ｘ＾、及び前記複数のパワースペクトルテンプレートλ(l)ωに基づいて、各周波数ωの前記分離行列Ｗω＾を更新する分離行列更新手段と、
を含む請求項１記載の音響信号解析装置。
【請求項３】
前記確率π(l)kτの事前分布を、ディレクレ分布とした請求項１又は２記載の音響信号解析装置。
【請求項４】
時間周波数解析手段によって、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘmτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘを出力し、
パラメータ初期値設定手段によって、Ｍ個の音源ｋの時間周波数成分Ｙkτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ(l)ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π(l)kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ2kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘτω（＝（Ｘ1τω，・・・，ＸMτω））に作用させて音源信号の時間周波数成分Ｙτω（＝（Ｙ1τω，・・・，ＹKτω））を得るための分離行列Ｗω＾の各々の初期値を設定し、
パラメータ更新手段によって、（k、τ、ω、ｚkτ）の全ての組み合わせにおける、ｚkτとσ2kγが与えられたときの、λ(ｚkτ)ω・σ2kγを分散とするガウス分布で表されたＹkτωの確率密度関数、及びπkτが与えられたときのｚkτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π(l)kτの事前確率、及び各周波数ωに対する前記分離行列Ｗω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾の事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾を更新し、
音源信号推定値更新手段によって、各周波数ωの前記分離行列Ｗω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新し、
終了判定手段によって、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う
音響信号解析方法。
【請求項５】
請求項１〜請求項３の何れか１項に記載の音響信号解析装置の各手段としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号解析装置、方法、及びプログラムに係り、特に、複数のマイクロホンから出力される音響信号の時系列データから、各音源の信号に分離する音響信号解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
音源の成分と音源からマイクロホンまでの伝達特性がともに未知のもとで、マイクロホン入力信号から個々の音源成分を分離抽出する技術をブラインド音源分離(Blind Source Separation;BSS)という。BSSでは観測信号だけから音源信号とその混合過程を推定する必要があるため、通常は音源に関して何らかの仮定を置き、これにより立てられる規準をもとに両未知変数を推定する最適化問題として定式化される。例えば、観測信号数が音源数以上の場合には、音源信号成分が優ガウス分布に従うという仮定のもとで分離フィルタを最尤推定する独立成分分析(Independent Component Analysis;ICA) と呼ぶ方法が有名である（非特許文献１）。周波数領域におけるICAは、音響信号のような時間遅れを有する混合系に対するBSSの有効なアプローチの1つであるが、周波数帯域毎に信号分離が行われるため、帯域毎の分離結果を同一音源毎にまとめる、いわゆるパーミュテーション問題を解く必要があった。
【０００３】
近年提案されている、ベクトル型変数を用いた独立成分分析(以下、独立ベクトル分析(Independent Vector Analysis;IVA))は、パーミュテーション問題を生じない分離手法として知られている（非特許文献２）。IVAにおいては、各音源のフレーム毎の複素スペクトルをベクトル型変数Yk,τ＾=(Yk,τ,1,...,Yk,τ,ω,...,Yk,τ,N)Tとして扱い(ただし、kは音源、τは時間フレーム、ωは周波数を表すインデックスである。)、そのノルムである
【０００４】
【数１】
2014048398.tif 000003
が優ガウス分布に従うという仮定のもとで最尤となる各帯域の分離フィルタが推定される。また、優ガウス分布の具体的な分布として時変ガウス分布(分散が時刻ごとに変化することを許容したガウス分布)を仮定したIVAを実現する方法も提案されている（非特許文献３）。
【０００５】
IVAによるブラインド音源分離の問題は、観測信号の時間周波数表現XTω＾=(X1Tω,...,XMTω)Tに対し、分離信号
【０００６】
【数２】
2014048398.tif 000004
の各要素が統計的に独立になるように分離行列Wω＾の推定する問題として定式化される。ただし、Ｙτω＾＝(Y1τω,...,Ykτω,...,YKτω)Tであり、各音源の時間周波数成分を表わしている。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】A. Hyv¨arinen, J. Karhumen, and E. Oja, Independent Component Analysis, John Wiley & Sons, 2001.
【非特許文献２】T. Kim, T. Eltoft and T. Lee, “Independent Vector Analysis: An Extension of ICA to Multivariate Components,” Proc. ICA, pp.165−172, 2006.
【非特許文献３】小野拓磨、小野順貴、嵯峨山茂樹、“音源のアクティベーションを事前情報とした独立ベクトル析による音源分離、”日本音響学会秋季研究発表会講演集、pp.613−614,Sep. 2011.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
以上のように、IVAでは、複素スペクトルをベクトル化したもののノルムが優ガウス分布に従うと仮定されるが、優ガウス分布として時変ガウス分布(フレームごとに分散が変化することを許容したガウス分布)を仮定した場合、この仮定は、各音源に対してどの周波数でもパワーの時間包絡が等しいと仮定したことになっている。しかし、一般に対象となる音源は音楽や音声など調波性という特別な性質を持つ場合がある。この場合、調波成分間の谷の部分での振幅包絡は必ずしも他の周波数と一致しない一方、基本周波数や倍音間のパワーの依存関係は音源分離の重要な手掛かりになり得る。また、仮定に適合しない信号（ピッチや周期といった調波性をもつ音源信号）では分離の精度に限界がある。
【０００９】
また、IVAでは、YTω＾の振幅包絡がどの周波数間でも類似する、という仮定ができるだけ満たされるようにWω＾が決定されるため、音源が調波構造を有する場合にはこの仮定が大きく崩れ、調波成分の谷間に該当する箇所に他音源の成分が混入したり、パーミュテーション不整合が起こることがあった。よって、もし各音源が調波構造を有するならば、音源モデルには、調波構造形のパワースペクトル密度を仮定した方が良いはずである。しかし、そのためには音源の基本周波数の情報が必要となるが、音源の基本周波数の情報は通常観測することができない。
【００１０】
本発明は、上記の事情を鑑みてなされたもので、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる音響信号解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記の目的を達成するために本発明に係る音響信号解析装置は、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘmτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘを出力する時間周波数解析手段と、Ｍ個の音源ｋの時間周波数成分Ｙkτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ(l)ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π(l)kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ2kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘτω（＝（Ｘ1τω，・・・，ＸMτω））に作用させて音源信号の時間周波数成分Ｙτω＾（＝（Ｙ1τω，・・・，ＹKτω））を得るための分離行列Ｗω＾の各々の初期値を設定するパラメータ初期値設定手段と、（k、τ、ω、ｚkτ）の全ての組み合わせにおける、ｚkτとσ2kγが与えられたときの、λ(ｚkτ)ω・σ2kγを分散とするガウス分布で表されたＹkτωの確率密度関数、及びπkτが与えられたときのｚkτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π(l)kτの事前確率、及び各周波数ωに対する前記分離行列Ｗω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾の事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾を更新するパラメータ更新手段と、各周波数ωの前記分離行列Ｗω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新する音源信号推定値更新手段と、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う終了判定手段と、を含んで構成されている。
【００１２】
本発明に係る音響信号解析方法は、時間周波数解析手段によって、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘmτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘ＾を出力し、パラメータ初期値設定手段によって、Ｍ個の音源ｋの時間周波数成分Ｙkτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ(l)ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π(l)kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ2kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘτω＾（＝（Ｘ1τω，・・・，ＸMτω））に作用させて音源信号の時間周波数成分Ｙτω＾（＝（Ｙ1τω，・・・，ＹKτω））を得るための分離行列Ｗω＾の各々の初期値を設定し、パラメータ更新手段によって、（k、τ、ω、ｚkτ）の全ての組み合わせにおける、ｚkτとσ2kγが与えられたときの、λ(ｚkτ)ω・σ2kγを分散とするガウス分布で表されたＹkτωの確率密度関数、及びπkτが与えられたときのｚkτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π(l)kτの事前確率、及び各周波数ωに対する前記分離行列Ｗω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗωの事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗω＾を更新し、音源信号推定値更新手段によって、各周波数ωの前記分離行列Ｗω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新し、終了判定手段によって、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う。
【００１３】
本発明に係るプログラムは、上記の音響信号解析装置の各手段としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１４】
以上説明したように、本発明の音響信号解析装置、方法、及びプログラムによれば、各周波数に対する各時間フレームにおいて調波構造を有するパワースペクトルテンプレートの各々が音源毎に選択される確率を要素にもつ三次元配列Π＾を用いて、観測時間周波数成分の三次元配列Ｘが与えられたときの三次元配列Π＾、各音源のパワーを要素にもつ二次元配列Σ＾、及び各周波数ωの分離行列Ｗω＾の事後確率を表す目的関数を最大化するように、三次元配列Π＾、二次元配列Σ＾、及び各周波数ωの分離行列Ｗω＾を更新することを繰り返すことにより、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる、という効果が得られる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施の形態に係る音響信号解析装置の構成を示す概略図である。
【図２】パワースペクトルテンプレートを示す図である。
【図３】本発明の実施の形態に係る音響信号解析装置における音響信号解析処理ルーチンの内容を示すフローチャートである。
【図４】（Ａ）ＳＤＲによる部屋E2Aにおける分離性能の評価結果を示すグラフ、及び（Ｂ）ＳＤＲによる部屋JR2における分離性能の評価結果を示すグラフである。
【発明を実施するための形態】
【００１６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する手法では、スペクトルテンプレートが時変ガウス分布における分散パラメータに組み込まれた音源の生成モデルを立てて、それをもとに、各音源において各フレームでどのスペクトルテンプレートが選ばれるべきかということと、各帯域の分離フィルタとを同時に推定する。
【００１７】
＜発明の原理＞
まず、本発明の原理について説明する。まず、音源の生成モデルについて説明する。
【００１８】
各フレームにおける基本周波数を潜在変数と見なし、潜在変数(基本周波数インデックスzkT)に応じて調波構造形のパワースペクトルテンプレートが一つ選択され、そのパワースペクトルをもとに音源の周波数成分が決定されるという生成プロセスを仮定する。これにより、従来のIVAにおいて仮定される音源モデルに比べ、調波性を有する音源を適切にモデル化できるはずであろうと考えられる。
【００１９】
l=1,...,Lを調波構造テンプレートのインデックスとし、λ(l)ωを特定の基本周波数の調波構造テンプレートとする（図２参照）。以上より、音源のGauss性を仮定すると、zkTが与えられたもとでの音源の生成モデルは、以下の（３）式で表される。
【００２０】
【数３】
2014048398.tif 000005
【００２１】
ただし、NC(・;μ,σ2)は平均μ、分散σ2の複素Gauss分布を表し、σkT2はフレームτにおける音源kのパワーである。一方、zkτ=lが選択される確率をπlkτとすると、潜在変数zkτは、以下の（４）式に示す離散分布で表される。
【００２２】
【数４】
2014048398.tif 000006
【００２３】
更に、πkτをスパースにするように誘導する目的でπkτのハイパー事前分布を以下の（５）式に示すDirichlet分布と仮定する。
【００２４】
【数５】
2014048398.tif 000007
【００２５】
以上の音源の生成モデルを、上記（２）式に組み込んだIVAの問題は、以下の（６）式に示す目的関数を最大化する問題に帰着する。
【００２６】
【数６】
2014048398.tif 000008
【００２７】
ただし、X＾={Xmτω},Π＾={π(l)kτ},Σ＾={σ2kτ},W＾={Wω},Θ＾={Π,Σ,W}である。Tは、時間フレーム数である。なお、記号に付された「＾」は、当該記号が行列または多次元配列またはベクトルであることを表わしている。
【００２８】
上記の（６）式の目的関数の大域最適解は解析的に得ることはできないが、補助関数法により局所解を効率的に探索することができる。ここで、Σlγ(l)kτω=1となる補助変数γkτωを用いれば、以下の（７）式のように補助関数を設計できる。
【００２９】
【数７】
2014048398.tif 000009
【００３０】
＜補助変数の更新＞
補助変数についての更新式はΣlγ(l)kτω=1の制約のもと、∂Q/∂γ(l)kτω=0を解けば以下の（８）式のとおりに得られる。
【００３１】
【数８】
2014048398.tif 000010
【００３２】
ただし、λはテンプレートであり、ｌは、テンプレートのインデックスを表す。
【００３３】
＜パラメータの更新＞
π(l)kτとσ2kτの更新式については、∂Q/∂π(l)kτ=0,∂Q/∂σ2kτ=0を解くことにより以下の（９）式、（１０）式のとおりに得られる。
【００３４】
【数９】
2014048398.tif 000011
【００３５】
Wω＾については、以下の（１１）式〜（１４）式のように1行ずつ更新する。
【００３６】
【数１０】
2014048398.tif 000012
【００３７】
ただし、ek^は、k番目のベクトル要素が1でそれ以外の要素は0の単位ベクトル（＝[0,...,1,...,0]T）である。
【００３８】
＜音源信号推定値の更新＞
以上で更新された分離行列を用いて音源信号の推定を以下の（１５）式のとおりに更新する。
【００３９】
【数１１】
2014048398.tif 000013
【００４０】
ただし、Yτ,ω＾=(Y1,τ,ω,...,YK,τ,ω)Tである。
【００４１】
＜システム構成＞
次に、Ｍ個（Ｍ≧２）のマイクロホンから得られた音響信号を解析して、既知のＫ個（Ｋ＜Ｍ）の音源信号に分離する音響信号解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００４２】
図１に示すように、本発明の実施の形態に係る音響信号解析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００４３】
音響信号解析装置１００は、入力部１０と、演算部２０と、記憶部３０と、出力部４０とを備えている。
【００４４】
入力部１０により、Ｍ個のマイクロホンから出力された音響信号（多チャンネル信号）の時系列データが入力される。記憶部３０は、入力部１０により入力された音響信号の時系列データを記憶する。また、記憶部３０は、後述する各処理での結果を記憶すると共に、各基本周波数の調波構造テンプレートλ(l)ω（l=1,...,L、ω＝１,・・・,Ｎ）、パラメータβを記憶している。
【００４５】
演算部２０は、時間周波数解析部２１と、初期設定部２２と、補助変数更新部２３と、パラメータ更新部２４と、時間周波数成分推定部２５と、終了判定部２６と、信号変換部２７とを備えている。また、パラメータ更新部２４は、確率更新部２４１と、時変ゲイン更新部２４２と、分離行列更新部２４３とを備えている。
【００４６】
時間周波数解析部２１は、各マイクロホンの時系列信号としての観測された音響信号を入力として、時間周波数成分（観測時間周波数成分）Ｘmτω（ｍ＝１，・・・，Ｍ、ω＝１,・・・,Ｎ,τ＝１,・・・,Ｔは、それぞれマイクロホン、周波数、時間フレームに対応するインデックスを示す。）を各（ｍ，ω，τ）の要素にもつ三次元配列Ｘを計算する。また、計算した時間周波数成分Ｘτωを、記憶部３０に記憶しておく。より詳細には、時間周波数解析部２１は、各マイクロホンｍについて、当該マイクロホンの音響信号の時系列データを入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより時間周波数成分Ｘmτωを計算し、時間周波数成分Ｘmτωを格納した行列（振幅スペクトログラム）Ｘ＾＝（Ｘmτω）M×N×Ｔを出力する。なお、時間周波数成分Ｘmτωは、ウェーブレット変換を用いて計算してもよい。
【００４７】
初期設定部２２は、後述する処理で用いる各パラメータπ(l)kτ、σ2kτ、Wω＾の初期値を設定する。初期値として、乱数を用いて適当な値を設定すればよい。また、初期設定部２２は、Wω＾の初期値と時間周波数成分Ｘmτωとに基づいて、上記（２）式に従って、分離信号Ｙτω＾の初期値を計算する。
【００４８】
補助変数更新部２３は、（ｌ、k、τ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているπ(l)kτ、σ2kτ、Ｙkτω、λ(l)ωに基づいて、上記（８）式に従って、観測時間周波数成分Ｘmτωが各インデックスｌのパワースペクトルテンプレートλ(l)ωに帰属する確率を示す補助変数γ(l)kτωを更新し、記憶部３０に格納する。
【００４９】
確率更新部２４１は、（ｌ、k、τ）の全ての組み合わせの各々について、記憶部３０に記憶されているγ(l)kτω、βに基づいて、上記（９）式に従って、パワースペクトルテンプレートλ(l)ωが選択される確率π(l)kτを更新し、記憶部３０に格納する。
【００５０】
時変ゲイン更新部２４２は、（k、τ）の全ての組み合わせの各々について、記憶部３０に記憶されているγ(l)kτω、Ｙkτω、λ(l)ωに基づいて、上記（１０）式に従って、フレームτにおける音源kのパワーσkT2を更新し、記憶部３０に格納する。
【００５１】
分離行列更新部２４３は、（k、τ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているγ(l)kτω、λ(l)ω、σkτ2に基づいて、上記（１１）式に従って、パラメータ~σkτω2を更新し、記憶部３０に格納する。また、分離行列更新部２４３は、（k、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているＸτω、~σkτω2に基づいて、上記（１２）式に従って、中間パラメータＶkωを更新し、記憶部３０に格納する。
【００５２】
そして、分離行列更新部２４３は、（k、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているWω＾、Ｖkωに基づいて、上記（１３）式、（１４）式に従って、ｗkω＾（音源ｋの音響信号の周波数成分ωに対する各マイクロホンｍの重みを要素にもつ重みベクトル）を更新し、記憶部３０に格納する。
【００５３】
時間周波数成分推定部２５は、（τ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているWω＾、Ｘτω＾に基づいて、上記（１５）式に従って、各音源の時間周波数成分を表わすベクトルＹτω＾＝(Y1τω,...,Ykτω,...,YKτω)Tを更新し、記憶部３０に格納する。
【００５４】
終了判定部２６は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、補助変数更新部２３、パラメータ更新部２４、及び時間周波数成分推定部２５の各処理を繰り返す。終了判定部２６は、終了条件を満足したと判定した場合には、信号変換部２７による処理に移行する。信号変換部２７は、記憶部３０に記憶されているＹτω＾を、各音源の音源信号に変換し、出力部４０により、音源信号を出力する。
【００５５】
終了条件としては、繰り返し回数ｓが予め定めた回数Ｓに達したことを用いればよい。なお、s-1回目のパラメータを用いたときの目的関数の値とs回目のパラメータを用いたときの目的関数の値との差が、予め定めた閾値よりも小さくなったことを、終了条件として用いてもよい。
【００５６】
＜音響信号解析装置の作用＞
次に、本実施の形態に係る音響信号解析装置１００の作用について説明する。まず、解析対象の信号として各マイクロホンからの音響信号の時系列データが音響信号解析装置１００に入力され、記憶部３０に格納される。そして、音響信号解析装置１００において、図３に示す音響信号解析処理ルーチンが実行される。
【００５７】
まず、ステップＳ１０１において、記憶部３０から、マイクロホン毎に、各時間フレームτ内の音響信号を読み込み、当該音響信号に対して、短時間フーリエ変換を用いた時間周波数分析を行った結果から、観測時間周波数成分Ｘmτωを各（ｍ，τ，ω）の要素にもつ三次元配列Ｘ＾を生成して、記憶部３０に記憶する。
【００５８】
そして、ステップＳ１０２において、乱数を用いて、パラメータΘ＾＝｛Π＾、Σ＾、Ｗ＾｝の初期値を設定して、記憶部３０に記憶すると共に、各音源信号の時間周波数成分Ｙkτωを各（ｋ，τ，ω）の要素にもつ三次元配列Ｙ＾を生成して、記憶部３０に記憶する。
【００５９】
次にステップＳ１０３では、上記ステップＳ１０２で設定されたパラメータΠ＾、Σ＾、Ｙ＾、又は後述するステップＳ１０４、Ｓ１０５、Ｓ１０７で更新されたパラメータΠ＾、Σ＾、Ｙ＾に基づいて、上記（８）式に従って、補助係数γ(l)kτωを各（ｌ，k，τ，ω）の組み合わせについて算出して、記憶部３０に格納する。
【００６０】
そして、ステップＳ１０４では、上記ステップＳ１０３で更新された補助係数γ(l)kτωに基づいて、上記（９）式に従って、パワースペクトルテンプレートλ(l)ωが選択される確率π(l)kτを各（ｌ，k，τ）の組み合わせについて算出して、記憶部３０に格納する。
【００６１】
ステップＳ１０５では、上記ステップＳ１０２で設定されたパラメータＹ＾、又はステップＳ１０７で更新されたパラメータＹ＾と、上記ステップＳ１０３で更新された補助係数γ(l)kτωに基づいて、上記（１０）式に従って、時間フレームτにおける音源kのパワーσkτ2を各（k，τ）の組み合わせについて算出して、記憶部３０に格納する。
【００６２】
次のステップＳ１０６では、上記ステップＳ１０２で設定されたパラメータＷ＾、又は前回のステップＳ１０６で更新されたパラメータＷ＾と、上記ステップＳ１０１で生成された三次元配列Ｘ＾と、上記ステップＳ１０５で更新されたパラメータΣ＾と、上記ステップＳ１０３で更新された補助係数γ(l)kτωと、基づいて、上記（１１）式〜（１４）式に従って、ｗkωを各（k，ω）の組み合わせについて算出して、記憶部３０に格納する。
【００６３】
そして、ステップＳ１０７では、上記ステップＳ１０６で更新されたパラメータＷと、上記ステップＳ１０１で生成された三次元配列Ｘ＾とに基づいて、上記（１５）式に従って、各音源信号の時間周波数成分Ｙkτωを各（ｋ，τ，ω）の要素にもつ三次元配列Ｙを算出して、記憶部３０に記憶する。
【００６４】
次のステップＳ１０８では、終了条件として、繰り返し回数ｓが、Ｓに到達したか否かを判定し、繰り返し回数ｓがＳに到達していない場合には、終了条件を満足していないと判断して、上記ステップＳ１０３へ戻り、上記ステップＳ１０３〜ステップＳ１０７の処理を繰り返す。一方、繰り返し回数ｓがＳに到達した場合には、終了条件を満足したと判断し、ステップＳ１０９で、上記ステップＳ１０７で最終的に更新された三次元配列Ｙに基づいて、各音源の音源信号を算出し、出力部４０により出力して、音響信号解析処理ルーチンを終了。
【００６５】
＜実験結果＞
次に、本実施の形態に係る手法の有用性を示す目的で、単旋律楽器を音源として用いたシミュレーションによる実験を行った結果について説明する。以下の表１に実験条件を示す。
【００６６】
【表１】
2014048398.tif 000014
【００６７】
予め録音された残響時間の異なる二種類の室内インパルス応答（参考文献（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,” Proc. LREC, pp. 965−968, 2000.）を参照）を単旋律楽器音(トランペット、サックス)とそれぞれ畳み込み加算することで観測信号とした。部屋E2A(残響時間T60=300ms)では、{−20°,0°}の方向から、部屋JR2(T60=470ms)では、{−30°,30°}の方向から音源が到来するとした。
【００６８】
本実施の形態に係る手法（以下、提案法と称する。）に用いた重みテンプレートは220Hz(A3)から半音ごとに24種類の倍音成分が1/√nで減衰する調波構造テンプレートと全ての周波数で同じ重みのテンプレート1種類の合計25 種類を用いた。提案法の分離行列の初期値として、時変Gauss分布に基づく独立ベクトル分析により分離行列を40回更新した値を用い、その後に提案法を40回反復し、復元信号を得た。提案法はπ(l)kTの事前分布としてDirichlet分布を考慮した手法(Prop(dir))と事前分布を一様にした手法(Prop(flat))の2種類を行った。Dirichlet分布を仮定した実験では、フレーム長が大きくなると事前分布の効用が小さくなるため効果を調節する目的でp(πkT)を周波数bin数乗した。また、提案法40回反復の間にp(πkT)のパラメータβを1から0.95まで徐々に小さくするアニーリング処理を行った。このとき更新の間にπ(l)kτが負になることがあったが、その場合十分小さい正数に置き換えた。
【００６９】
比較対象として、音源の生成モデルに時変Gauss分布を仮定した独立ベクトル分析による手法(IVA(TVG))音源の生成モデルにラプラス分布を仮定した独立ベクトル分析による手法(IVA(lap))（上記の非特許文献２参照）、独立成分分析による手法(ICA)(参考文献（H. Sawada, R. Mukai, S. Araki and S. Makino, “A Robust Approach to the Permutation Problem of Frequency-domain Blind source Separation,” Proc. ICASSP, pp. 381−384, 2003.）を参照)の3種類をそれぞれ分離行列の更新に80回行った。
【００７０】
SDR（参考文献（E. Vincent, R. Gribonval and C. F´evotte, “Performance Measurement in Blind Audio Source Separation,” IEEE Trans. ASLP, pp.1462−1469, 2006.）を参照）による評価の結果を、図４に示す。ＳＤＲによる部屋E2A(残響時間T60=300ms)における分離性能の評価結果を図４（Ａ）に示し、ＳＤＲによる部屋JR2(T60=470ms)における分離性能の評価結果を図４（Ｂ）に示す。2種類の環境どちらにおいても従来法に比べ提案法が優位な結果を得られた。特に、部屋E2AにおいてDirichlet事前分布を考慮することで、6dB程度の改善が見られた。
【００７１】
以上説明したように、本発明の実施の形態に係る音響信号解析装置によれば、各周波数及び各時間フレームにおいて調波構造を有するパワースペクトルテンプレートの各々が音源毎に選択される確率を要素にもつ三次元配列Π＾を用いて、観測時間周波数成分の三次元配列Ｘが与えられたときの三次元配列Π＾、各音源のパワーを要素にもつ二次元配列Σ＾、及び各周波数ωの分離行列Ｗω＾の事後確率を表す目的関数を最大化するように、三次元配列Π＾、二次元配列Σ＾、及び各周波数ωの分離行列Ｗω＾を更新することを繰り返すことにより、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる。
【００７２】
また、調波構造をもつパワースペクトルのテンプレートを複数個用意しておき、スペクトルテンプレートが時変ガウス分布における分散パラメータに組み込まれた音源の生成モデルを立てて、それをもとに、各音源において各時間フレームでどのスペクトルテンプレートが選ばれるべきかということと、各帯域の分離フィルタとを同時に推定することにより、音源信号が調波構造をもつ場合に、音源分離の精度を向上させることができる。
【００７３】
また、上記の非特許文献３とのアルゴリズムにおける相違点は、（１）γ(l)kτω、π(l)kTというパラメータとその更新式(上記（８）式、（９）式)が新たに加わったこと、（２）各音源のパワーσ2kτ が、上記（１０）式のようにγ(l)kTωを用いて算出されること、（３）中間変数Vkωが、上記（１１）式に従って算出される各音源の各時刻におけるパワースペクトル推定値~σ2kτω(周波数ωごとに異なる値をとりうる点が従来技術との重要な差異)を用いて上記（１２）式のように更新されることである。
【００７４】
また、複数のマイクロホンで取得した音響信号から混在する個々の音源信号を分離することができるため、ハンズフリーテレビ会議システムや会議録コンテンツの自動作成システムなどの応用が期待される。
【００７５】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００７６】
例えば、上述の音響信号解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【００７７】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【００７８】
１０   入力部
２０   演算部
２１   時間周波数解析部
２２   初期設定部
２３   補助変数更新部
２４   パラメータ更新部
２５   時間周波数成分推定部
２６   終了判定部
２７   信号変換部
３０   記憶部
４０   出力部
１００ 音響信号解析装置
２４１ 確率更新部
２４２ 時変ゲイン更新部
２４３ 分離行列更新部
【図１】
2014048398.tif 000015
【図２】
2014048398.tif 000016
【図３】
2014048398.tif 000017
【図４】
2014048398.tif 000018