説明

信号分解装置、方法、及びプログラム

【課題】複数の信号源からの信号各々が混合された観測信号を分解する場合に、分解成分と信号源との対応関係も把握する。
【解決手段】周波数領域変換部12で、観測信号を時間周波数成分xijに変換し、前処理部14で、時間周波数成分xijを、センサ間の相対的な関係を保持したまま一部が非負値となる時間周波数成分x’ijに変換する。初期値設定部16で、hik,tik,vkjの初期値を設定する。分解部18で、信号源とセンサとの伝達関数を表すベクトルhikとtikとvkjとの積和と、前処理後の時間周波数成分x’ijとを近似させるためのコスト関数Cを最小化するように、hik,tik,vkjのいずれかを更新する。時間領域変換部20で、コスト関数Cが最小化したときのパラメータを用いて信号分解すると共に、分解成分と信号源との対応関係を表すベクトルhikを出力する。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、信号分解装置、方法、及びプログラムに係り、特に、複数の信号源から受信した時系列信号を、周波数スペクトルパターンに基づいて分解する信号分解装置、方法、及びプログラムに関する。
【背景技術】
【0002】
従来、信号の分解に関しては、信号の時間周波数表現に基づく各種の分解方法が広く知られている。信号内に繰り返して現れるパターン等のような高度な構造を抽出する分解手法として非負値行列分解(NMF:Non‐negative Matrix Factorization)が知られている(例えば、非特許文献1参照)。NMFは、自動採譜や、モノラル混合信号からの音源の分離に適用されている(例えば、非特許文献2参照)。
【0003】
ここで、従来のNMFについて説明する。一個または複数の信号が混合され、一個のセンサで観測されたとする。その観測信号に対して短時間フーリエ変換を適用して、時間周波数表現xijを得る。ここで、i=1,・・・,Iは周波数ビンの番号、j=1,・・・,Jは時間フレームの番号を表す。xijはフーリエ変換の結果であるため複素数となるが、NMFを適用するために、このxijの絶対値|xij|を取り非負値とする。これらを要素とするI×J行列X、すなわち(X)ij=|xij|を考え、この行列Xに対してNMFを適用すると、下記(1)式のように、二つの非負値行列TとVとの積に近似的に分解される。
【0004】
【数1】

【0005】
ここで、T=[t,・・・,t]はI×K行列、V=[v,・・・,vはK ×J行列であり、・は転置を示す。Kは、NMFを適用する際に事前に決めておくパラメータであり、観測信号の分解数(何個に分解したいか)を指定するものである。このような分解は、下記(2)式のコスト関数を最小化することにより求められる。
【0006】
【数2】

【0007】
このNMFの結果において、tはk番目の分解結果の周波数スペクトルパターン、vはそのパターンがどの時間フレームで活性化されているかを表現している。NMFによる分解の後、各kに関して行列tを時間周波数表現と捉え、短時間フーリエ変換の逆変換を適用することで、観測信号をK個の信号に分解することができる。
【先行技術文献】
【非特許文献】
【0008】
【非特許文献1】D. D. Lee, H. S. Seung, “Learning the part of objects by non-negative matrix factorization,” Nature Vol.401, pp. 788-791, 1999.
【非特許文献2】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.
【発明の概要】
【発明が解決しようとする課題】
【0009】
従来のNMFによる信号の分解では、周波数スペクトルパターンに基づいて観測信号をK個に分解することができる。ただし、複数の信号源からの信号各々が混合された観測信号を分解する場合でも、同様にK個に分解されるが、どの分解成分がどの信号源に対応するか、という情報までは得ることができない、という問題がある。
【0010】
本発明は、上記の課題を解決するためになされたもので、複数の信号源からの信号各々が混合された観測信号を分解する場合に、分解成分と信号源との対応関係も取得することができる信号分解装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【0011】
上記目的を達成するために、本発明の信号分解装置は、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分xij(iは周波数ビン、jは時間フレームのインデックスである。)に変換し、該時間周波数成分xijを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換する変換手段と、基底となる周波数スペクトルパターンを表す非負値のパラメータtik、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータvkj、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータhikm(mは受信手段のインデックスである。)を要素とするベクトルhikの各パラメータの初期値を設定する設定手段と、前記ベクトルhikと前記パラメータtikと前記パラメータvkjとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分x’ijとの相違を表すコスト関数が小さくなり、かつ前記ベクトルhikについては、前記時間周波数成分x’ijに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する更新手段と、を含んで構成されている。
【0012】
本発明の信号分解装置によれば、変換手段が、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ受信手段の個数に対応した次元数の時間周波数成分xij(iは周波数ビン、jは時間フレームのインデックスである。)に変換し、時間周波数成分xijを、受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換する。
【0013】
そして、設定手段が、基底となる周波数スペクトルパターンを表す非負値のパラメータtik、周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータvkj、及び信号源各々と受信手段各々との間の伝達関数を表すパラメータhikm(mは受信手段のインデックスである。)を要素とするベクトルhikの各パラメータの初期値を設定する。
【0014】
そして、更新手段が、ベクトルhikとパラメータtikとパラメータvkjとの積和で表されたモデルと、変換手段により変換された時間周波数成分x’ijとの相違を表すコスト関数が小さくなり、かつベクトルhikについては、時間周波数成分x’ijに保持された受信手段間の相対的な関係に拘束されるように、各パラメータを更新する。
【0015】
このように、信号源各々と受信手段各々との間の伝達関数を表すベクトルhikを導入してNMFを適用することにより、得られたベクトルhikに基づいて、分解成分と信号源との対応関係を取得することができる。
【0016】
また、前記変換手段は、前記時間周波数成分xijに、該時間周波数成分xijの共役複素数を複素ベクトルの偏角を保持したまま絶対値を1にした値を乗じて、前記時間周波数成分x’ijに変換し、前記更新手段は、前記ベクトルhikを更新する際に、前記コスト関数が小さくなるベクトルhikを求め、求めたベクトルhikに、該ベクトルhikの共役複素数を複素ベクトルの偏角を保持したまま絶対値を1にした値を乗じたベクトルhikを用いて更新することができる。
【0017】
また、前記変換手段は、前記時間周波数成分xijと該時間周波数成分xijの共役転置成分との外積を、前記時間周波数成分x’ijとし、前記更新手段は、前記コスト関数を、前記ベクトルhikと該ベクトルhikの共役転置成分との外積を含んで表し、前記ベクトルhikを更新する際に、各パラメータを更新するために導入されたパラメータ行列を固有値分解して導出された最大固有値に対応する固有ベクトルを用いて更新することができる。
【0018】
また、本発明の信号分解方法は、変換手段と、設定手段と、更新手段とを含む信号分解装置における信号分解方法であって、前記変換手段は、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分xij(iは周波数ビン、jは時間フレームのインデックスである。)に変換し、該時間周波数成分xijを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換し、前記設定手段は、基底となる周波数スペクトルパターンを表す非負値のパラメータtik、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータvkj、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータhikm(mは受信手段のインデックスである。)を要素とするベクトルhikの各パラメータの初期値を設定し、前記更新手段は、前記ベクトルhikと前記パラメータtikと前記パラメータvkjとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分x’ijとの相違を表すコスト関数が小さくなり、かつ前記ベクトルhikについては、前記時間周波数成分x’ijに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する方法である。
【0019】
また、本発明の信号分解プログラムは、コンピュータを、上記信号分解装置を構成する各手段として機能させるためのプログラムである。
【発明の効果】
【0020】
以上説明したように、本発明の信号分解装置、方法、及びプログラムによれば、信号源各々と受信手段各々との間の伝達関数を表すベクトルhikを導入してNMFを適用することにより、得られたベクトルhikに基づいて、分解成分と信号源との対応関係を取得することができる、という効果が得られる。
【図面の簡単な説明】
【0021】
【図1】本実施の形態に係る信号分解装置の構成を示す概略図である。
【図2】分解部でのコスト関数最小化の処理を説明するための図である。
【図3】本実施の形態に係る信号分解装置における信号分解処理ルーチンの内容を示すフローチャートである。
【図4】観測信号の時間周波数表現を示す図である。
【図5】更新回数に対するコスト関数を示すグラフである。
【図6】非負値行列Tを示すグラフである。
【図7】非負値行列Vを示すグラフである。
【図8】伝達関数を表すベクトルhikについて、複素数の偏角arg(h・h)を示すグラフである。
【発明を実施するための形態】
【0022】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【0023】
まず、本実施の形態の信号分解装置による信号分解処理の概略について説明する。
【0024】
本実施の形態の信号分解装置は、空間的位置の異なる複数の信号源からの信号各々が混合された時系列信号(時間領域の観測信号)を、空間的位置の異なる複数のセンサ(受信手段)で観測し、これらの観測信号を分解する。このような状況では、信号源から各センサへの伝達関数(具体的には、減衰量及び時間遅れ量)に相当するものが信号源毎に異なるため、この違いを利用することで、分解成分と信号源との対応関係を取得する。
【0025】
具体的には、分解結果となる信号k(k=1,・・・,K)毎、さらに周波数ビンi毎に、各センサへの伝達関数を、パラメータhikmを要素とするベクトルhik=[hik1,・・・,hikMで表現する。ここで、mはセンサの番号、Mはセンサの数である。また、複数のセンサで観測された観測信号を用いるため、周波数ビン番号i、時間フレーム番号jの観測信号をM次元ベクトルxij=[xij1,・・・,xijMで表現する。
【0026】
そして、周波数スペクトルパターンを表すパラメータtikと、周波数スペクトルパターンの各時間フレームのゲインを表すパラメータvkiと、上記の伝達関数を表すベクトルhikとの積和を、観測信号を表すxijに近似させることにより、各パラメータを推定する。ここで推定されたベクトルhikに基づいて、分解成分と信号源との対応関係を取得することができる。各パラメータの推定は、後述するコスト関数を最小化する一連の手続きとして提供される。
【0027】
次に、第1の実施の形態の信号分解装置の構成について説明する。
【0028】
第1の実施の形態の信号分解装置10は、CPUと、RAMと、後述する信号分解処理ルーチンを実行するためのプログラムを記憶したROMとを備えたコンピュータで構成されている。このコンピュータは、機能的には、図1に示すように、周波数領域変換部12と、前処理部14と、初期値設定部16と、分解部18と、時間領域変換部20とを含んだ構成で表すことができる。
【0029】
周波数領域変換部12は、時系列信号としての観測信号を入力として、時間周波数成分xij(i(i=1,・・・,I)は周波数ビン、j(j=1,・・・J)は時間フレームのインデックスを示す。)を計算し、前処理部14に出力する。より詳細には、周波数領域変換部12は、観測信号を入力として、短時間フーリエ変換(Short-Time Fourier Transform;STFT)を用いて時間周波数解析を行うことにより時間周波数成分xijを計算する。なお、時間周波数成分xijは、ウェーブレット変換を用いて計算してもよい。
【0030】
前処理部14は、M個のセンサから基準センサを1つ選定する。ここでは、一般性を失うことなく、1番目のセンサとして話を進める。基準センサにおける時間周波数表現が全て非負値となるように、下記(3)式で示す前処理を全ての周波数ビン番号i、及び時間フレーム番号jについて行う。
【0031】
【数3】

【0032】
ここで、・は、共役複素数を示し、sign(・)は、複素数の偏角を保存したまま絶対値を1にする操作である。また、x’ijは、前処理後の時間周波数成分を表す。これにより、時間周波数成分xijは、センサ間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換される。
【0033】
初期値設定部16は、後述する分解部18で用いるhik,tik,vkjの初期値hik(0),tik(0),vkj(0)を設定する。初期値としては、乱数を用いて適当な値を設定することができる。
【0034】
分解部18は、ベクトルhikとパラメータtikとパラメータvkjとの積和で表されるモデルを、x’ijに近似させるためのコスト関数を定義し、このコスト関数を最小化するhik,tik,vkjを算出することで、結果的に観測信号を分解する。
【0035】
より具体的には、下記(4)式のようなコスト関数を定義する。
【0036】
【数4】

【0037】
ここで、‖a‖=Σm=1・aは、ベクトルのノルムの二乗である。また、tik、vkjはそれぞれ非負値であり、全てのi、j、kにより通常のNMFと同様の非負値行列T及びVを構成する。
【0038】
コスト関数を最小化する手続きは様々なものが考えられるが、本実施の形態では、補助関数を用いた方法について説明する。
【0039】
まず、(4)式のコスト関数Cに対して、下記(5)式に示す補助関数C’を定義する。
【0040】
【数5】

【0041】
ここで、rijk、sijkは新たに導入された補助変数であり、下記(6)式及び(7)式を満たすものである。
【0042】
【数6】

【0043】
このように定義された補助関数は、下記(8)式を常に満たす。
【0044】
【数7】

【0045】
さらに、下記(9)式を満たす補助変数により、C=C’となり、元のコスト関数と補助関数とが一致する。
【0046】
【数8】

【0047】
補助関数を用いたコスト関数の最小化方法は、以下の手続きを繰り返すことにより行う。
【0048】
(i) 補助変数を(9)式に従って更新し、補助関数をコスト関数に一致させる。
【0049】
(ii)補助関数を減少させるように、本変数であるhik,tik,vkjのいずれかを更新する。
【0050】
このコスト関数最小化の手続きを、図2を参照して、より詳細に説明する。
[補助変数の更新]
まず、(9)式に従って、補助変数を更新する。なお、rijk については、(6)式を満たす値であればよいが、下記(10)式に従って更新することが効率的である。
【0051】
【数9】

【0052】
[本変数の更新]
次に、本変数の更新について説明する。本実施の形態では、本変数の更新式の簡素化のため、‖hik‖=1の制約を設ける。もともと、hikとtikとの間には大きさの任意性があるため(例えばベクトルhikの長さを短くする必要があれば、tikを大きくすることで対処できる)、このような制約を設けても、コスト関数Cを最小化する自由度は変化しない。
【0053】
補助関数をそれぞれの本変数で偏微分し、これらを0に設定した式を解くことで、本変数の更新式が得られる。本変数であるパラメータtik及びvkjに関する更新式は、下記(11)式及び(12)式となる。
【0054】
【数10】

【0055】
ここで、Re[・]は、複素数の実部のみを取り出すものであり、・は、複素ベクトルの共役転置を示す。また、^x’ijは、(10)式で定義されているものである。
【0056】
本変数であるベクトルhikに関する更新式は、補助関数を最小化するものとして、下記(13)式に従って更新した後、下記(14)式に従って、前処理の条件を満たすための操作を行い、さらに、上記の‖hik‖=1の制約を満たすため、下記(15)式の操作を行う。
【0057】
【数11】

【0058】
[終了条件チェック]
補助変数の更新及び本変数の更新を行った後に、終了条件チェックを行う。終了条件の判定としては、例えば、一定回数(例えば100回)更新処理を行った場合に、終了条件を満たすと判定したり、コスト関数の変化量が所定の閾値以下になった場合に、終了条件を満たすと判定したりすることができる。
【0059】
時間領域変換部20は、分解部18で算出されたhik,tik,vkjに基づいて、分解されたk番目の信号毎に、短時間フーリエ変換の逆変換によって時間周波数表現を時間領域の信号に変換して、出力する。また、時間領域変換部20は、分解成分と信号源との対応関係を表す情報として、分解部18で算出されたベクトルhikの値もあわせて出力する。
【0060】
次に、第1の実施の形態の信号分解装置10の作用について説明する。まず、分解対象の観測信号が信号分解装置10に入力され、図示しないメモリに格納される。そして、信号分解装置10において、図3に示す信号分解処理ルーチンが実行される。
【0061】
ステップ100で、メモリから、あるフレーム内の観測信号を読み込み、観測信号を、短時間フーリエ変換により、時間周波数成分xijに変換する。
【0062】
次に、ステップ102で、上記ステップ100で変換した時間周波数成分xijに対して、(3)式で示す前処理を施して、時間周波数成分x’ijに変換する。これにより、時間周波数成分xijは、センサ間の相対的な関係を保持したまま一部を非負値とした時間周波数成分x’ijに変換される。
【0063】
次に、ステップ104で、乱数を用いて、hik,tik,vkjの初期値hik(0),tik(0),vkj(0)を設定する。
【0064】
次に、ステップ106で、上記ステップ104で設定した初期値hik(0),tik(0),vkj(0)を用いて、(9)式に従って補助変数sijk及びrijkを更新する。これにより、(4)式で示すコスト関数Cと(5)式で示す補助関数C’とが一致する。
【0065】
次に、ステップ108で、本変数であるhik,tik,vkjのいずれかを更新する。パラメータtikを更新する場合には、(11)式に従って更新する。また、パラメータvkjを更新する場合には、(12)式に従って更新する。また、ベクトルhikを更新する場合には、(13)式に従って更新した後、(14)式に従って、前処理の条件を満たすための操作を行い、さらに、‖hik‖=1の制約を満たすため、(15)式の操作を行う。hik,tik,vkjのいずれかの更新により、補助関数C’を減少させることができる。また、ベクトルhikについては、複素ベクトルの偏角の情報を保持した前処理の条件に拘束されて、更新される。
【0066】
次に、ステップ110で、一定回数(例えば100回)更新処理を行った、コスト関数の変化量が所定の閾値以下になった等の終了条件が成立するか否かを判定する。終了条件が成立していない場合には、ステップ106へ戻って処理を繰り返す。ステップ106では、前のステップ108で更新されたhik,tik,vを用いて、(9)式に従って補助変数sijk及びrijkを更新する。これにより、C=C’となるが、前のステップ108の処理で、補助関数C’が減少するようにhik,tik,vのいずれかが更新されているため、コスト関数Cも減少する。すなわち、ステップ110で、終了条件が成立したと判定されるまでステップ106及び108の処理を繰り返すことにより、コスト関数Cを最小化させることができる。
【0067】
上記ステップ110で、終了条件が成立したと判定された場合には、ステップ112へ移行し、上記ステップ108で算出されたhik,tik,vkjに基づいて、分解されたk番目の信号毎に、短時間フーリエ変換の逆変換によって時間周波数表現を時間領域の信号に変換して、出力する。また、分解成分と信号源との対応関係を表す情報として、上記ステップ108で算出されたベクトルhikの値もあわせて出力して、処理を終了する。
【0068】
次に、第2の実施の形態について説明する。第2の実施の形態の信号分解装置では、パラメータの更新式を導出するためのコスト関数が第1の実施の形態とは異なる。なお、第2の実施の形態の信号分解装置について、第1の実施の形態の信号分解装置10と同様の構成については、同一の符号を付して説明を省略する。
【0069】
図1に示すように、第2の実施の形態の信号分解装置210を構成するコンピュータは、機能的には、周波数領域変換部12と、前処理部214と、初期値設定部16と、分解部218と、時間領域変換部20とを含んだ構成で表すことができる。
【0070】
前処理部214は、全ての周波数ビン番号i、及び時間フレーム番号jについて、下記(16)式に示す時間周波数成分xijのベクトルの外積を計算する。
【0071】
【数12】

【0072】
ここで、x’’ijは、前処理後の時間周波数成分を表す。これにより、時間周波数成分xijは、センサ間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換される。
【0073】
分解部218は、ベクトルhikとパラメータtikとパラメータvkjとの積和で表されるモデルを、x’’ijに近似させるためのコスト関数を定義し、このコスト関数を最小化するhik,tik,vkjを算出することで、結果的に観測信号を分解する。
【0074】
より具体的には、下記(17)式のようなコスト関数を定義する。
【0075】
【数13】

【0076】
ここで、ここで、‖A‖=Σm=1Σn=1mm・amnは、行列のフロベニウスノルムの二乗である。また、第1の実施の形態におけるコスト関数Cと同様に、tik、vkjはそれぞれ非負値であり、全てのi、j、kにより通常のNMFと同様の非負値行列T及びVを構成する。
【0077】
第2の実施の形態においても、補助関数を用いてコスト関数を最小化する手続きについて説明する。
【0078】
まず、(17)式のコスト関数Cに対して、下記(18)式に示す補助関数C’を定義する。
【0079】
【数14】

【0080】
ここで、rijk及びSijkは、第1の実施の形態の補助関数C’のrijk及びsijkと同様に、新たに導入された補助変数であり、(6)式及び下記(19)式を満たすものである。
【0081】
【数15】

【0082】
このように定義された補助関数は、(8)式を常に満たし、さらに、下記(20)式を満たす補助変数により、C=C’となり、元のコスト関数と補助関数とが一致する。
【0083】
【数16】

【0084】
このコスト関数最小化の手続きを、図2を参照して、より詳細に説明する。
【0085】
[補助変数の更新]
第1の実施の形態の場合と同様に、(20)式に従って、補助変数を更新する。
【0086】
[本変数の更新]
コスト関数Cにおける本変数であるパラメータtik及びvkjに関する更新式は、下記(21)式及び(22)式となる。
【0087】
【数17】

【0088】
ここで、^x’’ijは、式(19) で定義されているものである。
【0089】
本変数であるベクトルhikに関する更新式は、上記の‖hik‖=1の制約を満たしながら補助関数を最小化するものとして、下記(23)式が導かれる。
【0090】
【数18】

【0091】
この解は、行列の最大固有値に対応する固有ベクトルとして与えられるため、行列j ^x’’ijijkを固有値分解して最大固有値を導出し、その最大固有値に対応する固有ベクトルをhikとする。
【0092】
第2の実施の形態の信号分解装置210の作用については、用いるコスト関数が異なる点を除いて、第1の実施の形態の信号分解装置10における信号分解処理と同様であるため、説明を省略する。
【0093】
次に、本実施の形態の信号解析装置の効果について、下記表1に示す条件の下で行った実験について説明する。
【0094】
【表1】

【0095】
本実験では、2つの楽器音を同時に鳴らした時の混合音を2つのマイクロホンで観測し、その観測信号を時間周波数表現に変換した。図4に、1つ目のマイクロホンでの観測信号の時間周波数表現を示す。これをパラメータK=10として、第1の実施の形態で説明したコスト関数Cを用いて分解した。また、ランダムな初期値で本変数であるhik,tik,vkjを初期化し、100回の繰り返しを通じてコスト関数の最小化を行った。図5に示す通り、良好にコスト関数が最小化された。
【0096】
結果として、図6、7及び8に示す非負値行列tik及びvkjとベクトルhikが得られた。非負値行列tik及びvkjは通常のNMFでも得られるものである。特徴的なパターンを示してはいるが、信号源の物理的な位置に関わる情報は含んでいない。
【0097】
一方、本実施の形態のポイントである信号源からセンサまでの伝達関数を表現したベクトルhikは、信号源の物理的な位置に関わる情報を含んでいる。図7は、ベクトルh=[h,hについて、複素数の偏角arg(h・h)を示したものである。なお、対応するtの値が0.1を下回るものは表示していない。これは、信号源から2つのセンサ(マイクロホン)への到達時間差と周波数とを乗じたものに比例する値を示している。従って、これらから、k=9に関わる分解成分のみがマイクロホン1に先に到達し、その他の分解成分はマイクロホン2に先に到達していることが分かる。
【0098】
また、上記のようにマイクロホンへの到達時間差(順序)を判定し、それに基づいて、分解成分と信号源との対応関係を得る方法の他に、図8の結果をk毎の点の分布と見て、その分布をクラスタリングすることで、10個を1個と9個とに分解することができる。
【0099】
以上説明したように、本実施の形態の信号分解装置によれば、複数のセンサを用いて観測した際の、信号源の空間的位置による各センサへの伝達関数を表すベクトルhikを導入してNMFを適用することにより、得られたベクトルhikのk毎の集合を見ることにより、分解成分と信号源との対応関係を取得することができる。
【0100】
また、上記第1の実施の形態の変形例として、(9)式に(10)式のrijkを代入すると、下記(24)式となる。
【0101】
【数19】

【0102】
さらに、更新式の簡素化のため‖hik‖=1の制約を設けた場合、下記(25)式の関係を導くことができる。
【0103】
【数20】

【0104】
(25)式を、(11)式及び(12)式に適用すると、下記(26)式及び(27)式となり、sijkの計算を行う必要がなく、eijのみで十分となる更新式を導くことができる。
【0105】
【数21】

【0106】
次に、ベクトルhikに関する更新式は、(13)式の代わりに、補助関数を最小化するものとして、下記(28)式に従って更新した後、前処理の条件を満たすための操作((14)式)を行い、さらに上記の‖hik‖=1の制約を満たすために(15)式の操作を行う。
【0107】
【数22】

【0108】
また、同様に、第2の実施の形態の変形例として、(20)式に式(10)のrijkを代入すると、下記(29)式となる。
【0109】
【数23】

【0110】
さらに、更新式の簡素化のため‖hik‖=1の制約を設けた場合、下記(30)式の関係を導くことができる。
【0111】
【数24】

【0112】
これを(21)式及び(22)式に適用すると下記(31)式及び(32)式となり、Sijkの計算を行う必要がなく、Eijのみで十分となる更新式を導くことができる。
【0113】
【数25】

【0114】
次に、ベクトルhikに関する更新式は、上記の‖hik‖=1の制約を満たしながら補助関数を最小化するものとして、下記(33)式が導かれる。
【0115】
【数26】

【0116】
この解は、行列の最大固有値に対応する固有ベクトルとして与えられるため、下記(34)式に示す行列を固有値分解して最大固有値を導出し、その最大固有値に対応する固有ベクトルをhikとする。
【0117】
【数27】

【0118】
上記の簡素化によると、第1の実施の形態の変形例の補助変数の更新式は、下記(35)式、第2の実施の形態の変形例の補助変数の更新式は、下記(36)式となる。
【0119】
【数28】

【0120】
(35)式または(36)式によると、補助関数を用いたコスト関数の最小化方法は、以下の手続きに変形することができる。
【0121】
(i) 補助変数を(35)式(または(36)式)に従って更新し、補助関数をコスト関数に一致させる。
【0122】
(ii)補助関数を減少させるように、本変数であるhik,tik,vkjのいずれかを更新する。
【0123】
なお、本実施の形態の信号分解装置を利用して、例えば、様々な楽器が演奏されている状況を複数のマイクロホン(音に対するセンサ)で観測し、各楽器の音に分解することが可能となる。
【0124】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【0125】
例えば、上述の信号分解装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、WWWシステムを利用している場合であれば、ホームページ提供環境(あるいは表示環境)も含むものとする。
【0126】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【0127】
10、210 信号分解装置
12 周波数領域変換部
14、214 前処理部
16 初期値設定部
18、218 分解部
20 時間領域変換部

【特許請求の範囲】
【請求項1】
複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分xij(iは周波数ビン、jは時間フレームのインデックスである。)に変換し、該時間周波数成分xijを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換する変換手段と、
基底となる周波数スペクトルパターンを表す非負値のパラメータtik、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータvkj、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータhikm(mは受信手段のインデックスである。)を要素とするベクトルhikの各パラメータの初期値を設定する設定手段と、
前記ベクトルhikと前記パラメータtikと前記パラメータvkjとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分x’ijとの相違を表すコスト関数が小さくなり、かつ前記ベクトルhikについては、前記時間周波数成分x’ijに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する更新手段と、
を示す信号分解装置。
【請求項2】
前記変換手段は、前記時間周波数成分xijに、該時間周波数成分xijの共役複素数を複素ベクトルの偏角を保持したまま絶対値を1にした値を乗じて、前記時間周波数成分x’ijに変換し、
前記更新手段は、前記ベクトルhikを更新する際に、前記コスト関数が小さくなるベクトルhikを求め、求めたベクトルhikに、該ベクトルhikの共役複素数を複素ベクトルの偏角を保持したまま絶対値を1にした値を乗じたベクトルhikを用いて更新する
請求項1記載の信号分解装置。
【請求項3】
前記変換手段は、前記時間周波数成分xijと該時間周波数成分xijの共役転置成分との外積を、前記時間周波数成分x’ijとし、
前記更新手段は、前記コスト関数を、前記ベクトルhikと該ベクトルhikの共役転置成分との外積を含んで表し、前記ベクトルhikを更新する際に、各パラメータを更新するために導入されたパラメータ行列を固有値分解して導出された最大固有値に対応する固有ベクトルを用いて更新する
請求項1記載の信号分解装置。
【請求項4】
変換手段と、設定手段と、更新手段とを含む信号分解装置における信号分解方法であって、
前記変換手段は、複数の受信手段により受信した空間的位置が異なる複数の信号源からの時系列信号を、所定の時間フレーム毎の周波数スペクトルを表し、かつ前記受信手段の個数に対応した次元数の時間周波数成分xij(iは周波数ビン、jは時間フレームのインデックスである。)に変換し、該時間周波数成分xijを、前記受信手段間の相対的な関係を保持したまま一部を非負値に変換した時間周波数成分x’ijに変換し、
前記設定手段は、基底となる周波数スペクトルパターンを表す非負値のパラメータtik、前記周波数スペクトルパターンの各時間フレームでのゲインを表す非負値のパラメータvkj、及び前記信号源各々と前記受信手段各々との間の伝達関数を表すパラメータhikm(mは受信手段のインデックスである。)を要素とするベクトルhikの各パラメータの初期値を設定し、
前記更新手段は、前記ベクトルhikと前記パラメータtikと前記パラメータvkjとの積和で表されたモデルと、前記変換手段により変換された時間周波数成分x’ijとの相違を表すコスト関数が小さくなり、かつ前記ベクトルhikについては、前記時間周波数成分x’ijに保持された前記受信手段間の相対的な関係に拘束されるように、前記各パラメータを更新する
信号分解方法。
【請求項5】
コンピュータを、請求項1〜請求項3のいずれか1項記載の信号分解装置を構成する各手段として機能させるための信号分解プログラム。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図5】
image rotate

【図4】
image rotate

【図6】
image rotate

【図7】
image rotate

【図8】
image rotate


【公開番号】特開2012−242493(P2012−242493A)
【公開日】平成24年12月10日(2012.12.10)
【国際特許分類】
【出願番号】特願2011−110463(P2011−110463)
【出願日】平成23年5月17日(2011.5.17)
【出願人】(000004226)日本電信電話株式会社 (13,992)