(57)【特許請求の範囲】
【請求項１】
複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分ｘｉｊ（ｉは周波数ビン、ｊは時間フレームのインデックスである。）に変換し、該時間周波数成分ｘｉｊを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換する変換手段と、
基底となる周波数スペクトルパターンを表す非負値のパラメータｔｉｋ、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータｖｋｊ、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータｈｉｋｍ（ｍは受信手段のインデックスである。）を要素とするベクトルｈｉｋの各パラメータの初期値を設定する設定手段と、
前記ベクトルｈｉｋと前記パラメータｔｉｋと前記パラメータｖｋｊとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分ｘ’ｉｊとの相違を表すコスト関数が小さくなり、かつ前記ベクトルｈｉｋについては、前記時間周波数成分ｘ’ｉｊに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する更新手段と、
を示す信号分解装置。
【請求項２】
前記変換手段は、前記時間周波数成分ｘｉｊに、前記時間周波数成分ｘｉｊのうち、予め選定されたｗ番目の受信手段に対応する次元の要素である時間周波数表現ｘｉｊｗの共役複素数を複素ベクトルの偏角を保持したまま絶対値を１にした値を乗じて、前記時間周波数成分ｘ’ｉｊに変換し、
前記更新手段は、前記ベクトルｈｉｋを更新する際に、前記コスト関数が小さくなるベクトルｈｉｋを求め、求めたベクトルｈｉｋに、前記ベクトルｈｉｋのうち、前記予め選定されたｗ番目の受信手段との間の伝達関数を表すパラメータｈｉｋｗの共役複素数を複素ベクトルの偏角を保持したまま絶対値を１にした値を乗じたベクトルｈｉｋを用いて更新する
請求項１記載の信号分解装置。
【請求項３】
前記変換手段は、前記時間周波数成分ｘｉｊと該時間周波数成分ｘｉｊの共役転置成分との外積を、前記時間周波数成分ｘ’ｉｊとし、
前記更新手段は、前記コスト関数を、前記ベクトルｈｉｋと該ベクトルｈｉｋの共役転置成分との外積を含んで表し、前記ベクトルｈｉｋを更新する際に、各パラメータを更新するために導入されたパラメータ行列を固有値分解して導出された最大固有値に対応する固有ベクトルを用いて更新する
請求項１記載の信号分解装置。
【請求項４】
変換手段と、設定手段と、更新手段とを含む信号分解装置における信号分解方法であって、
前記変換手段は、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分ｘｉｊ（ｉは周波数ビン、ｊは時間フレームのインデックスである。）に変換し、該時間周波数成分ｘｉｊを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換し、
前記設定手段は、基底となる周波数スペクトルパターンを表す非負値のパラメータｔｉｋ、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータｖｋｊ、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータｈｉｋｍ（ｍは受信手段のインデックスである。）を要素とするベクトルｈｉｋの各パラメータの初期値を設定し、
前記更新手段は、前記ベクトルｈｉｋと前記パラメータｔｉｋと前記パラメータｖｋｊとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分ｘ’ｉｊとの相違を表すコスト関数が小さくなり、かつ前記ベクトルｈｉｋについては、前記時間周波数成分ｘ’ｉｊに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する
信号分解方法。
【請求項５】
コンピュータを、請求項１〜請求項３のいずれか１項記載の信号分解装置を構成する各手段として機能させるための信号分解プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号分解装置、方法、及びプログラムに係り、特に、複数の信号源から受信した時系列信号を、周波数スペクトルパターンに基づいて分解する信号分解装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来、信号の分解に関しては、信号の時間周波数表現に基づく各種の分解方法が広く知られている。信号内に繰り返して現れるパターン等のような高度な構造を抽出する分解手法として非負値行列分解（ＮＭＦ：Non‐negative Matrix Factorization）が知られている（例えば、非特許文献１参照）。ＮＭＦは、自動採譜や、モノラル混合信号からの音源の分離に適用されている（例えば、非特許文献２参照）。
【０００３】
ここで、従来のＮＭＦについて説明する。一個または複数の信号が混合され、一個のセンサで観測されたとする。その観測信号に対して短時間フーリエ変換を適用して、時間周波数表現ｘｉｊを得る。ここで、ｉ＝１，・・・，Ｉは周波数ビンの番号、ｊ＝１，・・・，Ｊは時間フレームの番号を表す。ｘｉｊはフーリエ変換の結果であるため複素数となるが、ＮＭＦを適用するために、このｘｉｊの絶対値｜ｘｉｊ｜を取り非負値とする。これらを要素とするＩ×Ｊ行列Ｘ、すなわち（Ｘ）ｉｊ＝｜ｘｉｊ｜を考え、この行列Ｘに対してＮＭＦを適用すると、下記（１）式のように、二つの非負値行列ＴとＶとの積に近似的に分解される。
【０００４】
【数１】
0005520883.tif 000002
【０００５】
ここで、Ｔ＝［ｔ１，・・・，ｔＫ］はＩ×Ｋ行列、Ｖ＝［ｖ１，・・・，ｖＫ］ＴはＫ  ×Ｊ行列であり、・Ｔは転置を示す。Ｋは、ＮＭＦを適用する際に事前に決めておくパラメータであり、観測信号の分解数（何個に分解したいか）を指定するものである。このような分解は、下記（２）式のコスト関数を最小化することにより求められる。
【０００６】
【数２】
0005520883.tif 000003
【０００７】
このＮＭＦの結果において、ｔｋはｋ番目の分解結果の周波数スペクトルパターン、ｖｋＴはそのパターンがどの時間フレームで活性化されているかを表現している。ＮＭＦによる分解の後、各ｋに関して行列ｔｋｖｋＴを時間周波数表現と捉え、短時間フーリエ変換の逆変換を適用することで、観測信号をＫ個の信号に分解することができる。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】D. D. Lee, H. S. Seung, “Learning the part of objects by non-negative matrix factorization,” Nature Vol.401, pp. 788-791, 1999.
【非特許文献２】P. Smaragdis, J. C. Brown, “Non-Negative Matrix Factorization for Music Transcription,”lnProc. 2003 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA2003), pp. 177-180,2003.
【発明の概要】
【発明が解決しようとする課題】
【０００９】
従来のＮＭＦによる信号の分解では、周波数スペクトルパターンに基づいて観測信号をＫ個に分解することができる。ただし、複数の信号源からの信号各々が混合された観測信号を分解する場合でも、同様にＫ個に分解されるが、どの分解成分がどの信号源に対応するか、という情報までは得ることができない、という問題がある。
【００１０】
本発明は、上記の課題を解決するためになされたもので、複数の信号源からの信号各々が混合された観測信号を分解する場合に、分解成分と信号源との対応関係も取得することができる信号分解装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明の信号分解装置は、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分ｘｉｊ（ｉは周波数ビン、ｊは時間フレームのインデックスである。）に変換し、該時間周波数成分ｘｉｊを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換する変換手段と、基底となる周波数スペクトルパターンを表す非負値のパラメータｔｉｋ、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータｖｋｊ、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータｈｉｋｍ（ｍは受信手段のインデックスである。）を要素とするベクトルｈｉｋの各パラメータの初期値を設定する設定手段と、前記ベクトルｈｉｋと前記パラメータｔｉｋと前記パラメータｖｋｊとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分ｘ’ｉｊとの相違を表すコスト関数が小さくなり、かつ前記ベクトルｈｉｋについては、前記時間周波数成分ｘ’ｉｊに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する更新手段と、を含んで構成されている。
【００１２】
本発明の信号分解装置によれば、変換手段が、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ受信手段の個数に対応した次元数の時間周波数成分ｘｉｊ（ｉは周波数ビン、ｊは時間フレームのインデックスである。）に変換し、時間周波数成分ｘｉｊを、受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換する。
【００１３】
そして、設定手段が、基底となる周波数スペクトルパターンを表す非負値のパラメータｔｉｋ、周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータｖｋｊ、及び信号源各々と受信手段各々との間の伝達関数を表すパラメータｈｉｋｍ（ｍは受信手段のインデックスである。）を要素とするベクトルｈｉｋの各パラメータの初期値を設定する。
【００１４】
そして、更新手段が、ベクトルｈｉｋとパラメータｔｉｋとパラメータｖｋｊとの積和で表されたモデルと、変換手段により変換された時間周波数成分ｘ’ｉｊとの相違を表すコスト関数が小さくなり、かつベクトルｈｉｋについては、時間周波数成分ｘ’ｉｊに保持された受信手段間の相対的な関係に拘束されるように、各パラメータを更新する。
【００１５】
このように、信号源各々と受信手段各々との間の伝達関数を表すベクトルｈｉｋを導入してＮＭＦを適用することにより、得られたベクトルｈｉｋに基づいて、分解成分と信号源との対応関係を取得することができる。
【００１６】
また、前記変換手段は、前記時間周波数成分ｘｉｊに、前記時間周波数成分ｘｉｊのうち、予め選定されたｗ番目の受信手段に対応する次元の要素である時間周波数表現ｘｉｊｗの共役複素数を複素ベクトルの偏角を保持したまま絶対値を１にした値を乗じて、前記時間周波数成分ｘ’ｉｊに変換し、前記更新手段は、前記ベクトルｈｉｋを更新する際に、前記コスト関数が小さくなるベクトルｈｉｋを求め、求めたベクトルｈｉｋに、前記ベクトルｈｉｋのうち、前記予め選定されたｗ番目の受信手段との間の伝達関数を表すパラメータｈｉｋｗの共役複素数を複素ベクトルの偏角を保持したまま絶対値を１にした値を乗じたベクトルｈｉｋを用いて更新することができる。
【００１７】
また、前記変換手段は、前記時間周波数成分ｘｉｊと該時間周波数成分ｘｉｊの共役転置成分との外積を、前記時間周波数成分ｘ’ｉｊとし、前記更新手段は、前記コスト関数を、前記ベクトルｈｉｋと該ベクトルｈｉｋの共役転置成分との外積を含んで表し、前記ベクトルｈｉｋを更新する際に、各パラメータを更新するために導入されたパラメータ行列を固有値分解して導出された最大固有値に対応する固有ベクトルを用いて更新することができる。
【００１８】
また、本発明の信号分解方法は、変換手段と、設定手段と、更新手段とを含む信号分解装置における信号分解方法であって、前記変換手段は、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分ｘｉｊ（ｉは周波数ビン、ｊは時間フレームのインデックスである。）に変換し、該時間周波数成分ｘｉｊを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換し、前記設定手段は、基底となる周波数スペクトルパターンを表す非負値のパラメータｔｉｋ、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータｖｋｊ、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータｈｉｋｍ（ｍは受信手段のインデックスである。）を要素とするベクトルｈｉｋの各パラメータの初期値を設定し、前記更新手段は、前記ベクトルｈｉｋと前記パラメータｔｉｋと前記パラメータｖｋｊとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分ｘ’ｉｊとの相違を表すコスト関数が小さくなり、かつ前記ベクトルｈｉｋについては、前記時間周波数成分ｘ’ｉｊに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する方法である。
【００１９】
また、本発明の信号分解プログラムは、コンピュータを、上記信号分解装置を構成する各手段として機能させるためのプログラムである。
【発明の効果】
【００２０】
以上説明したように、本発明の信号分解装置、方法、及びプログラムによれば、信号源各々と受信手段各々との間の伝達関数を表すベクトルｈｉｋを導入してＮＭＦを適用することにより、得られたベクトルｈｉｋに基づいて、分解成分と信号源との対応関係を取得することができる、という効果が得られる。
【図面の簡単な説明】
【００２１】
【図１】本実施の形態に係る信号分解装置の構成を示す概略図である。
【図２】分解部でのコスト関数最小化の処理を説明するための図である。
【図３】本実施の形態に係る信号分解装置における信号分解処理ルーチンの内容を示すフローチャートである。
【図４】観測信号の時間周波数表現を示す図である。
【図５】更新回数に対するコスト関数を示すグラフである。
【図６】非負値行列Ｔを示すグラフである。
【図７】非負値行列Ｖを示すグラフである。
【図８】伝達関数を表すベクトルｈｉｋについて、複素数の偏角ａｒｇ（ｈ１＊・ｈ２）を示すグラフである。
【発明を実施するための形態】
【００２２】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２３】
まず、本実施の形態の信号分解装置による信号分解処理の概略について説明する。
【００２４】
本実施の形態の信号分解装置は、空間的位置の異なる複数の信号源からの信号各々が混合された時系列信号（時間領域の観測信号）を、空間的位置の異なる複数のセンサ（受信手段）で観測し、これらの観測信号を分解する。このような状況では、信号源から各センサへの伝達関数（具体的には、減衰量及び時間遅れ量）に相当するものが信号源毎に異なるため、この違いを利用することで、分解成分と信号源との対応関係を取得する。
【００２５】
具体的には、分解結果となる信号ｋ（ｋ＝１，・・・，Ｋ）毎、さらに周波数ビンｉ毎に、各センサへの伝達関数を、パラメータｈｉｋｍを要素とするベクトルｈｉｋ＝［ｈｉｋ１，・・・，ｈｉｋＭ］Ｔで表現する。ここで、ｍはセンサの番号、Ｍはセンサの数である。また、複数のセンサで観測された観測信号を用いるため、周波数ビン番号ｉ、時間フレーム番号ｊの観測信号をＭ次元ベクトルｘｉｊ＝［ｘｉｊ１，・・・，ｘｉｊＭ］Ｔ  で表現する。
【００２６】
そして、周波数スペクトルパターンを表すパラメータｔｉｋと、周波数スペクトルパターンの各時間フレームのゲインを表すパラメータｖｋｉと、上記の伝達関数を表すベクトルｈｉｋとの積和を、観測信号を表すｘｉｊに近似させることにより、各パラメータを推定する。ここで推定されたベクトルｈｉｋに基づいて、分解成分と信号源との対応関係を取得することができる。各パラメータの推定は、後述するコスト関数を最小化する一連の手続きとして提供される。
【００２７】
次に、第１の実施の形態の信号分解装置の構成について説明する。
【００２８】
第１の実施の形態の信号分解装置１０は、ＣＰＵと、ＲＡＭと、後述する信号分解処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成されている。このコンピュータは、機能的には、図１に示すように、周波数領域変換部１２と、前処理部１４と、初期値設定部１６と、分解部１８と、時間領域変換部２０とを含んだ構成で表すことができる。
【００２９】
周波数領域変換部１２は、時系列信号としての観測信号を入力として、時間周波数成分ｘｉｊ（ｉ（ｉ＝１，・・・，Ｉ）は周波数ビン、ｊ（ｊ＝１，・・・Ｊ）は時間フレームのインデックスを示す。）を計算し、前処理部１４に出力する。より詳細には、周波数領域変換部１２は、観測信号を入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより時間周波数成分ｘｉｊを計算する。なお、時間周波数成分ｘｉｊは、ウェーブレット変換を用いて計算してもよい。
【００３０】
前処理部１４は、Ｍ個のセンサから基準センサを１つ選定する。ここでは、一般性を失うことなく、１番目のセンサとして話を進める。基準センサにおける時間周波数表現が全て非負値となるように、下記（３）式で示す前処理を全ての周波数ビン番号ｉ、及び時間フレーム番号ｊについて行う。
【００３１】
【数３】
0005520883.tif 000004
【００３２】
ここで、・＊は、共役複素数を示し、ｓｉｇｎ（・）は、複素数の偏角を保存したまま絶対値を１にする操作である。また、ｘ’ｉｊは、前処理後の時間周波数成分を表す。これにより、時間周波数成分ｘｉｊは、センサ間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換される。
【００３３】
初期値設定部１６は、後述する分解部１８で用いるｈｉｋ，ｔｉｋ，ｖｋｊの初期値ｈｉｋ（０），ｔｉｋ（０），ｖｋｊ（０）を設定する。初期値としては、乱数を用いて適当な値を設定することができる。
【００３４】
分解部１８は、ベクトルｈｉｋとパラメータｔｉｋとパラメータｖｋｊとの積和で表されるモデルを、ｘ’ｉｊに近似させるためのコスト関数を定義し、このコスト関数を最小化するｈｉｋ，ｔｉｋ，ｖｋｊを算出することで、結果的に観測信号を分解する。
【００３５】
より具体的には、下記（４）式のようなコスト関数を定義する。
【００３６】
【数４】
0005520883.tif 000005
【００３７】
ここで、‖ａ‖２  ＝Σｍ＝１Ｍａｍ・ａｍ＊は、ベクトルのノルムの二乗である。また、ｔｉｋ、ｖｋｊはそれぞれ非負値であり、全てのｉ、ｊ、ｋにより通常のＮＭＦと同様の非負値行列Ｔ及びＶを構成する。
【００３８】
コスト関数を最小化する手続きは様々なものが考えられるが、本実施の形態では、補助関数を用いた方法について説明する。
【００３９】
まず、（４）式のコスト関数Ｃ１に対して、下記（５）式に示す補助関数Ｃ１’を定義する。
【００４０】
【数５】
0005520883.tif 000006
【００４１】
ここで、ｒｉｊｋ、ｓｉｊｋは新たに導入された補助変数であり、下記（６）式及び（７）式を満たすものである。
【００４２】
【数６】
0005520883.tif 000007
【００４３】
このように定義された補助関数は、下記（８）式を常に満たす。
【００４４】
【数７】
0005520883.tif 000008
【００４５】
さらに、下記（９）式を満たす補助変数により、Ｃ１＝Ｃ１’となり、元のコスト関数と補助関数とが一致する。
【００４６】
【数８】
0005520883.tif 000009
【００４７】
補助関数を用いたコスト関数の最小化方法は、以下の手続きを繰り返すことにより行う。
【００４８】
（ｉ）  補助変数を（９）式に従って更新し、補助関数をコスト関数に一致させる。
【００４９】
（ｉｉ）補助関数を減少させるように、本変数であるｈｉｋ，ｔｉｋ，ｖｋｊのいずれかを更新する。
【００５０】
このコスト関数最小化の手続きを、図２を参照して、より詳細に説明する。
［補助変数の更新］
まず、（９）式に従って、補助変数を更新する。なお、ｒｉｊｋ  については、（６）式を満たす値であればよいが、下記（１０）式に従って更新することが効率的である。
【００５１】
【数９】
0005520883.tif 000010
【００５２】
［本変数の更新］
次に、本変数の更新について説明する。本実施の形態では、本変数の更新式の簡素化のため、‖ｈｉｋ‖＝１の制約を設ける。もともと、ｈｉｋとｔｉｋとの間には大きさの任意性があるため（例えばベクトルｈｉｋの長さを短くする必要があれば、ｔｉｋを大きくすることで対処できる）、このような制約を設けても、コスト関数Ｃ１を最小化する自由度は変化しない。
【００５３】
補助関数をそれぞれの本変数で偏微分し、これらを０に設定した式を解くことで、本変数の更新式が得られる。本変数であるパラメータｔｉｋ及びｖｋｊに関する更新式は、下記（１１）式及び（１２）式となる。
【００５４】
【数１０】
0005520883.tif 000011
【００５５】
ここで、Ｒｅ［・］は、複素数の実部のみを取り出すものであり、・Ｈは、複素ベクトルの共役転置を示す。また、＾ｘ’ｉｊは、（１０）式で定義されているものである。
【００５６】
本変数であるベクトルｈｉｋに関する更新式は、補助関数を最小化するものとして、下記（１３）式に従って更新した後、下記（１４）式に従って、前処理の条件を満たすための操作を行い、さらに、上記の‖ｈｉｋ‖＝１の制約を満たすため、下記（１５）式の操作を行う。
【００５７】
【数１１】
0005520883.tif 000012
【００５８】
［終了条件チェック］
補助変数の更新及び本変数の更新を行った後に、終了条件チェックを行う。終了条件の判定としては、例えば、一定回数（例えば１００回）更新処理を行った場合に、終了条件を満たすと判定したり、コスト関数の変化量が所定の閾値以下になった場合に、終了条件を満たすと判定したりすることができる。
【００５９】
時間領域変換部２０は、分解部１８で算出されたｈｉｋ，ｔｉｋ，ｖｋｊに基づいて、分解されたｋ番目の信号毎に、短時間フーリエ変換の逆変換によって時間周波数表現を時間領域の信号に変換して、出力する。また、時間領域変換部２０は、分解成分と信号源との対応関係を表す情報として、分解部１８で算出されたベクトルｈｉｋの値もあわせて出力する。
【００６０】
次に、第１の実施の形態の信号分解装置１０の作用について説明する。まず、分解対象の観測信号が信号分解装置１０に入力され、図示しないメモリに格納される。そして、信号分解装置１０において、図３に示す信号分解処理ルーチンが実行される。
【００６１】
ステップ１００で、メモリから、あるフレーム内の観測信号を読み込み、観測信号を、短時間フーリエ変換により、時間周波数成分ｘｉｊに変換する。
【００６２】
次に、ステップ１０２で、上記ステップ１００で変換した時間周波数成分ｘｉｊに対して、（３）式で示す前処理を施して、時間周波数成分ｘ’ｉｊに変換する。これにより、時間周波数成分ｘｉｊは、センサ間の相対的な関係を保持したまま一部を非負値とした時間周波数成分ｘ’ｉｊに変換される。
【００６３】
次に、ステップ１０４で、乱数を用いて、ｈｉｋ，ｔｉｋ，ｖｋｊの初期値ｈｉｋ（０），ｔｉｋ（０），ｖｋｊ（０）を設定する。
【００６４】
次に、ステップ１０６で、上記ステップ１０４で設定した初期値ｈｉｋ（０），ｔｉｋ（０），ｖｋｊ（０）を用いて、（９）式に従って補助変数ｓｉｊｋ及びｒｉｊｋを更新する。これにより、（４）式で示すコスト関数Ｃ１と（５）式で示す補助関数Ｃ１’とが一致する。
【００６５】
次に、ステップ１０８で、本変数であるｈｉｋ，ｔｉｋ，ｖｋｊのいずれかを更新する。パラメータｔｉｋを更新する場合には、（１１）式に従って更新する。また、パラメータｖｋｊを更新する場合には、（１２）式に従って更新する。また、ベクトルｈｉｋを更新する場合には、（１３）式に従って更新した後、（１４）式に従って、前処理の条件を満たすための操作を行い、さらに、‖ｈｉｋ‖＝１の制約を満たすため、（１５）式の操作を行う。ｈｉｋ，ｔｉｋ，ｖｋｊのいずれかの更新により、補助関数Ｃ１’を減少させることができる。また、ベクトルｈｉｋについては、複素ベクトルの偏角の情報を保持した前処理の条件に拘束されて、更新される。
【００６６】
次に、ステップ１１０で、一定回数（例えば１００回）更新処理を行った、コスト関数の変化量が所定の閾値以下になった等の終了条件が成立するか否かを判定する。終了条件が成立していない場合には、ステップ１０６へ戻って処理を繰り返す。ステップ１０６では、前のステップ１０８で更新されたｈｉｋ，ｔｉｋ，ｖｋを用いて、（９）式に従って補助変数ｓｉｊｋ及びｒｉｊｋを更新する。これにより、Ｃ１＝Ｃ１’となるが、前のステップ１０８の処理で、補助関数Ｃ１’が減少するようにｈｉｋ，ｔｉｋ，ｖｋのいずれかが更新されているため、コスト関数Ｃ１も減少する。すなわち、ステップ１１０で、終了条件が成立したと判定されるまでステップ１０６及び１０８の処理を繰り返すことにより、コスト関数Ｃ１を最小化させることができる。
【００６７】
上記ステップ１１０で、終了条件が成立したと判定された場合には、ステップ１１２へ移行し、上記ステップ１０８で算出されたｈｉｋ，ｔｉｋ，ｖｋｊに基づいて、分解されたｋ番目の信号毎に、短時間フーリエ変換の逆変換によって時間周波数表現を時間領域の信号に変換して、出力する。また、分解成分と信号源との対応関係を表す情報として、上記ステップ１０８で算出されたベクトルｈｉｋの値もあわせて出力して、処理を終了する。
【００６８】
次に、第２の実施の形態について説明する。第２の実施の形態の信号分解装置では、パラメータの更新式を導出するためのコスト関数が第１の実施の形態とは異なる。なお、第２の実施の形態の信号分解装置について、第１の実施の形態の信号分解装置１０と同様の構成については、同一の符号を付して説明を省略する。
【００６９】
図１に示すように、第２の実施の形態の信号分解装置２１０を構成するコンピュータは、機能的には、周波数領域変換部１２と、前処理部２１４と、初期値設定部１６と、分解部２１８と、時間領域変換部２０とを含んだ構成で表すことができる。
【００７０】
前処理部２１４は、全ての周波数ビン番号ｉ、及び時間フレーム番号ｊについて、下記（１６）式に示す時間周波数成分ｘｉｊのベクトルの外積を計算する。
【００７１】
【数１２】
0005520883.tif 000013
【００７２】
ここで、ｘ’’ｉｊは、前処理後の時間周波数成分を表す。これにより、時間周波数成分ｘｉｊは、センサ間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分ｘ’ｉｊに変換される。
【００７３】
分解部２１８は、ベクトルｈｉｋとパラメータｔｉｋとパラメータｖｋｊとの積和で表されるモデルを、ｘ’’ｉｊに近似させるためのコスト関数を定義し、このコスト関数を最小化するｈｉｋ，ｔｉｋ，ｖｋｊを算出することで、結果的に観測信号を分解する。
【００７４】
より具体的には、下記（１７）式のようなコスト関数を定義する。
【００７５】
【数１３】
0005520883.tif 000014
【００７６】
ここで、ここで、‖Ａ‖Ｆ２＝Σｍ＝１ＭΣｎ＝１Ｍａｍｍ・ａｍｎ＊は、行列のフロベニウスノルムの二乗である。また、第１の実施の形態におけるコスト関数Ｃ１と同様に、ｔｉｋ、ｖｋｊはそれぞれ非負値であり、全てのｉ、ｊ、ｋにより通常のＮＭＦと同様の非負値行列Ｔ及びＶを構成する。
【００７７】
第２の実施の形態においても、補助関数を用いてコスト関数を最小化する手続きについて説明する。
【００７８】
まず、（１７）式のコスト関数Ｃ２に対して、下記（１８）式に示す補助関数Ｃ２’を定義する。
【００７９】
【数１４】
0005520883.tif 000015
【００８０】
ここで、ｒｉｊｋ及びＳｉｊｋは、第１の実施の形態の補助関数Ｃ１’のｒｉｊｋ及びｓｉｊｋと同様に、新たに導入された補助変数であり、（６）式及び下記（１９）式を満たすものである。
【００８１】
【数１５】
0005520883.tif 000016
【００８２】
このように定義された補助関数は、（８）式を常に満たし、さらに、下記（２０）式を満たす補助変数により、Ｃ２＝Ｃ２’となり、元のコスト関数と補助関数とが一致する。
【００８３】
【数１６】
0005520883.tif 000017
【００８４】
このコスト関数最小化の手続きを、図２を参照して、より詳細に説明する。
【００８５】
［補助変数の更新］
第１の実施の形態の場合と同様に、（２０）式に従って、補助変数を更新する。
【００８６】
［本変数の更新］
コスト関数Ｃ２における本変数であるパラメータｔｉｋ及びｖｋｊに関する更新式は、下記（２１）式及び（２２）式となる。
【００８７】
【数１７】
0005520883.tif 000018
【００８８】
ここで、＾ｘ’’ｉｊは、式（１９）  で定義されているものである。
【００８９】
本変数であるベクトルｈｉｋに関する更新式は、上記の‖ｈｉｋ‖＝１の制約を満たしながら補助関数を最小化するものとして、下記（２３）式が導かれる。
【００９０】
【数１８】
0005520883.tif 000019
【００９１】
この解は、行列の最大固有値に対応する固有ベクトルとして与えられるため、行列ｊ    ＾ｘ’’ｉｊＳｉｊｋを固有値分解して最大固有値を導出し、その最大固有値に対応する固有ベクトルをｈｉｋとする。
【００９２】
第２の実施の形態の信号分解装置２１０の作用については、用いるコスト関数が異なる点を除いて、第１の実施の形態の信号分解装置１０における信号分解処理と同様であるため、説明を省略する。
【００９３】
次に、本実施の形態の信号解析装置の効果について、下記表１に示す条件の下で行った実験について説明する。
【００９４】
【表１】
0005520883.tif 000020
【００９５】
本実験では、２つの楽器音を同時に鳴らした時の混合音を２つのマイクロホンで観測し、その観測信号を時間周波数表現に変換した。図４に、１つ目のマイクロホンでの観測信号の時間周波数表現を示す。これをパラメータＫ＝１０として、第１の実施の形態で説明したコスト関数Ｃ１を用いて分解した。また、ランダムな初期値で本変数であるｈｉｋ，ｔｉｋ，ｖｋｊを初期化し、１００回の繰り返しを通じてコスト関数の最小化を行った。図５に示す通り、良好にコスト関数が最小化された。
【００９６】
結果として、図６、７及び８に示す非負値行列ｔｉｋ及びｖｋｊとベクトルｈｉｋが得られた。非負値行列ｔｉｋ及びｖｋｊは通常のＮＭＦでも得られるものである。特徴的なパターンを示してはいるが、信号源の物理的な位置に関わる情報は含んでいない。
【００９７】
一方、本実施の形態のポイントである信号源からセンサまでの伝達関数を表現したベクトルｈｉｋは、信号源の物理的な位置に関わる情報を含んでいる。図７は、ベクトルｈ＝［ｈ１，ｈ２］Ｔについて、複素数の偏角ａｒｇ（ｈ１＊・ｈ２）を示したものである。なお、対応するｔの値が０．１を下回るものは表示していない。これは、信号源から２つのセンサ（マイクロホン）への到達時間差と周波数とを乗じたものに比例する値を示している。従って、これらから、ｋ＝９に関わる分解成分のみがマイクロホン１に先に到達し、その他の分解成分はマイクロホン２に先に到達していることが分かる。
【００９８】
また、上記のようにマイクロホンへの到達時間差（順序）を判定し、それに基づいて、分解成分と信号源との対応関係を得る方法の他に、図８の結果をｋ毎の点の分布と見て、その分布をクラスタリングすることで、１０個を１個と９個とに分解することができる。
【００９９】
以上説明したように、本実施の形態の信号分解装置によれば、複数のセンサを用いて観測した際の、信号源の空間的位置による各センサへの伝達関数を表すベクトルｈｉｋを導入してＮＭＦを適用することにより、得られたベクトルｈｉｋのｋ毎の集合を見ることにより、分解成分と信号源との対応関係を取得することができる。
【０１００】
また、上記第１の実施の形態の変形例として、（９）式に（１０）式のｒｉｊｋを代入すると、下記（２４）式となる。
【０１０１】
【数１９】
0005520883.tif 000021
【０１０２】
さらに、更新式の簡素化のため‖ｈｉｋ‖＝１の制約を設けた場合、下記（２５）式の関係を導くことができる。
【０１０３】
【数２０】
0005520883.tif 000022
【０１０４】
（２５）式を、（１１）式及び（１２）式に適用すると、下記（２６）式及び（２７）式となり、ｓｉｊｋの計算を行う必要がなく、ｅｉｊのみで十分となる更新式を導くことができる。
【０１０５】
【数２１】
0005520883.tif 000023
【０１０６】
次に、ベクトルｈｉｋに関する更新式は、（１３）式の代わりに、補助関数を最小化するものとして、下記（２８）式に従って更新した後、前処理の条件を満たすための操作（（１４）式）を行い、さらに上記の‖ｈｉｋ‖＝１の制約を満たすために（１５）式の操作を行う。
【０１０７】
【数２２】
0005520883.tif 000024
【０１０８】
また、同様に、第２の実施の形態の変形例として、（２０）式に式（１０）のｒｉｊｋを代入すると、下記（２９）式となる。
【０１０９】
【数２３】
0005520883.tif 000025
【０１１０】
さらに、更新式の簡素化のため‖ｈｉｋ‖＝１の制約を設けた場合、下記（３０）式の関係を導くことができる。
【０１１１】
【数２４】
0005520883.tif 000026
【０１１２】
これを（２１）式及び（２２）式に適用すると下記（３１）式及び（３２）式となり、Ｓｉｊｋの計算を行う必要がなく、Ｅｉｊのみで十分となる更新式を導くことができる。
【０１１３】
【数２５】
0005520883.tif 000027
【０１１４】
次に、ベクトルｈｉｋに関する更新式は、上記の‖ｈｉｋ‖＝１の制約を満たしながら補助関数を最小化するものとして、下記（３３）式が導かれる。
【０１１５】
【数２６】
0005520883.tif 000028
【０１１６】
この解は、行列の最大固有値に対応する固有ベクトルとして与えられるため、下記（３４）式に示す行列を固有値分解して最大固有値を導出し、その最大固有値に対応する固有ベクトルをｈｉｋとする。
【０１１７】
【数２７】
0005520883.tif 000029
【０１１８】
上記の簡素化によると、第１の実施の形態の変形例の補助変数の更新式は、下記（３５）式、第２の実施の形態の変形例の補助変数の更新式は、下記（３６）式となる。
【０１１９】
【数２８】
0005520883.tif 000030
【０１２０】
（３５）式または（３６）式によると、補助関数を用いたコスト関数の最小化方法は、以下の手続きに変形することができる。
【０１２１】
（ｉ）  補助変数を（３５）式（または（３６）式）に従って更新し、補助関数をコスト関数に一致させる。
【０１２２】
（ｉｉ）補助関数を減少させるように、本変数であるｈｉｋ，ｔｉｋ，ｖｋｊのいずれかを更新する。
【０１２３】
なお、本実施の形態の信号分解装置を利用して、例えば、様々な楽器が演奏されている状況を複数のマイクロホン（音に対するセンサ）で観測し、各楽器の音に分解することが可能となる。
【０１２４】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１２５】
例えば、上述の信号分解装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１２６】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１２７】
１０、２１０  信号分解装置
１２   周波数領域変換部
１４、２１４  前処理部
１６   初期値設定部
１８、２１８  分解部
２０   時間領域変換部
【図１】
0005520883.tif 000031
【図２】
0005520883.tif 000032
【図３】
0005520883.tif 000033
【図５】
0005520883.tif 000034
【図４】
0005520883.tif 000035
【図６】
0005520883.tif 000036
【図７】
0005520883.tif 000037
【図８】
0005520883.tif 000038