JP4769238B2 - 信号分離装置、信号分離方法、プログラム及び記録媒体 - Google Patents

信号分離装置、信号分離方法、プログラム及び記録媒体 Download PDF

Info

Publication number
JP4769238B2
JP4769238B2 JP2007218612A JP2007218612A JP4769238B2 JP 4769238 B2 JP4769238 B2 JP 4769238B2 JP 2007218612 A JP2007218612 A JP 2007218612A JP 2007218612 A JP2007218612 A JP 2007218612A JP 4769238 B2 JP4769238 B2 JP 4769238B2
Authority
JP
Japan
Prior art keywords
signal
unit
frequency
class
permutation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2007218612A
Other languages
English (en)
Other versions
JP2009053349A (ja
Inventor
宏 澤田
章子 荒木
昭二 牧野
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nippon Telegraph and Telephone Corp
Original Assignee
Nippon Telegraph and Telephone Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nippon Telegraph and Telephone Corp filed Critical Nippon Telegraph and Telephone Corp
Priority to JP2007218612A priority Critical patent/JP4769238B2/ja
Publication of JP2009053349A publication Critical patent/JP2009053349A/ja
Application granted granted Critical
Publication of JP4769238B2 publication Critical patent/JP4769238B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Obtaining Desirable Characteristics In Audible-Bandwidth Transducers (AREA)
  • Circuit For Audible Band Transducer (AREA)

Description

本発明は、信号処理の技術分野に属し、特に、複数の信号が空間内で混合されたものから、源信号をできるだけ正確に復元する信号分離技術に関する。
［ブラインド信号分離］
まず、ブラインド信号分離の定式化を行う。扱う信号はあるサンプリング周波数fsでサンプリングされ、離散的に表現されるものとする。また、N個の信号が混合されてM個のセンサで観測されたとする。以下では、信号の発生源からセンサまでの距離により信号が減衰・遅延し、また壁や床などによる反射／残響が発生する状況を扱う。このような状況での混合は、源信号sn(t)(n=1,...,N)を発した信号源nからセンサm(m=1,...,M)へのインパルス応答hmn(r)による畳み込み混合
となる。ここでtはサンプリング時間を、rは掃引（時間シフトした信号のサンプル値それぞれに異なる係数を作用させる操作）のための変数を、それぞれ示している。一般的なインパルス応答hmn(r)の形状は、適当な時間経過後にパルス的な強い応答を持ち、時間と共に減衰していくものである。ブラインド信号分離の目的は、源信号s1(t),...,sN(t)やインパルス応答h11(r),...,h1N(r),...,hM1(r),...,hMN(r)を知らずに、観測信号x1(t),...,xM(t)のみから、源信号s1(t),...,sN(t)にそれぞれ対応する分離信号y1(t),...,yN(t)を求めることにある。
［周波数領域における信号分離］
畳み込み混合の問題は扱いが繁雑である。よって、上述の式（１）に短時間離散フーリエ変換(DFT: Discrete Fourier Transform)を施して、信号を周波数領域に変換した上で分離の操作を行うことが有効である。上述の式（１）に短時間離散フーリエ変換を適用して周波数毎の時間系列を求めると以下のようになる。
ここでfは周波数であり、f=0, fs/L ,・・・, fs(L-1)/Lと離散化されている（fsはサンプリング周波数）。また、τは時間インデックスであり、jは虚数単位である。また、g(r)は窓関数である。ハニング窓g(r)=(1+cos(2π・r /L))/2などのg(0)にパワーの中心を持つ窓関数を用いることで、Xm(f,τ)は時間τを中心とする観測信号xm(t)の周波数特性を表現する。なお、Xm(f,τ)はＬサンプルにわたる情報を含んでいるため、すべての時間tを時間インデックスτとしてXm(f,τ)を求める必要はなく、適当な間隔で時間インデックスτを設定してXm(f,τ)を求める。
式(1)で示される時間領域での畳み込み混合を周波数領域での表現に変換すると、
と各周波数での単純混合に近似でき、分離の操作が単純になる。なお、Hmn(f)は信号源nからセンサmまでの周波数応答であり、Sn(f,τ)は式(2)と同様な式に従って源信号sn(t)に短時間離散フーリエ変換を施したものである。式(3)をベクトル表記すると、
となる。ここで、X(f,τ)=[X1(f,τ),...,XM(f,τ)]TはXm(f,τ)を要素とする観測信号ベクトルであり、Hn(f)=[H1n(f),...,HMn(f)]Tは信号源nからセンサmまでの周波数応答Hmn(f)を要素とするベクトルである。なお、[・]Tは[・]の転置を示す。
［スパース性に基づく信号分離］
ブラインド信号分離法の一つにスパース性に基づいて信号分離を行う方法がある（例えば、特許文献１等参照）。このスパース性に基づく信号分離の場合、信号源の数Nとセンサの数Mの関係にかかわらず（M≧2であればN>MでもN≦Mでも良い）、同一の仕組みにより分離の処理が可能である。これは独立成分分析（ICA: Independent Component Analysis）を用いたブラインド信号分離（例えば、特許文献2、非特許文献5等参照）と対比される。ICAを用いる場合には、信号源の数Nがセンサの数Mを超えない（N≦M）ことが強く望まれる。これに対して、スパース性に基づく信号分離では、センサの数に関する要求条件がより緩くなっており、より広い適用範囲が見込まれる。
ただし、スパース性による信号分離が有効に働くためには、対象となる源信号がスパース性を持つことが条件となる。スパース性とは、ほとんどの場合において信号の振幅が零に近く、大きな振幅となるのは稀であるという性質である。例えば、周波数領域での音声信号にはスパース性を十分に確認できる。スパース性を満たす源信号s1(t),...,sN(t)の場合、式(4)の混合過程は、さらに
X(f,τ)=Hp(f)・Sp(f,τ) …(5)
と近似表現できる。ここで添字ｐは、時間周波数(f,τ)に依存したものとなる。ほとんどの場合において信号源の振幅が零に近いため、個々の時間周波数(f,τ)において最も振幅の大きい源信号Sp(f,τ)に関わる項だけで、式(4)が近似されている。
スパース性による信号分離では、各時間周波数(f,τ)において、どの源信号Sp(f,τ)の振幅が最も大きいかを推定する。言い替えると、観測信号ベクトルX(f,τ)をN個のクラスC1,...,CNに分類（クラスタリング）し、クラスCkには、源信号Sk(f,τ)が最も支配的な観測信号ベクトルX(f,τ)が属するようにする。ここで、サンプルXを観測した後における、サンプルXがクラスCkに属する事象の事後確率をP(Ck |X)で表記する。そのような事後確率をなんらかの方法で推定できれば、分離信号Yn(f,τ)は、例えば時間周波数マスキング
により構成できる。ここで、Jは分離信号を構成するために用いる基準センサの添字であり、1からMの範囲から選択される。
［クラス分類／事後確率の計算方法］
クラス分類或いは事後確率計算の方法として、信号源の方向や位置に相当する値を推定し、それに基づいてすべての時間周波数(f,τ)に関する観測信号ベクトルX(f,τ)を一気にクラス分類したり事後確率計算したりする方法が提案されている（例えば、特許文献１、特許文献3、非特許文献1、非特許文献2、非特許文献3等参照）。これらの方法では、1)信号源毎に推定された方向や位置に相当する値に従って、すべての時間周波数(f,τ)に関する観測信号ベクトルX(f,τ)をN個のクラスに分類したり事後確率を計算したりするプロセスと、2)分類された観測信号を元に、方向や位置に相当する値を信号源毎に再推定するプロセスとを行う。これらは、反射や残響の影響が比較的少ない場合には有効に働く。
また、観測信号ベクトルX(f,τ)をN個のクラスヘ分類する操作を、周波数毎に行う方法も提案されている。この場合は、ある周波数でのi番目のクラスと別の周波数でのi番目のクラスが、同じ信号源に対応するものかどうかが不明となる。従って、その後、同一信号源に対応するクラスを全周波数に渡って同定する必要がある。この問題は、ICAを用いたブラインド信号分離におけるパーミュテーション問題とほぼ同じである。これに対し、各周波数での分類結果から各信号源の方向や位置に相当する値を推定し、その推定結果に基づいてパーミュテーション問題を解決する方法（例えば、特許文献２、特許文献３、非特許文献４、非特許文献５等参照）や、周波数毎の分離信号エンベロープの相関係数の類似度に基づいてパーミュテーション問題を解決する方法 (例えば、非特許文献５、非特許文献６、非特許文献７、非特許文献８等）がこれまで用いられてきた。
WO2005/024788 WO2004/079388 WO2006/085537 O. Yilmaz and S. Rickard, "Blind separation of speech mixtures via time-frequency masking," IEEE Trans. Signal Processing, vol. 52, no. 7, pp. 1830-1847, July 2004. M. Mandel, D. Ellis, and T. Jebara, "An EM Algorithm for Localizing Multiple Sound Sources in Reverberant Environments," Advances in Neural Information Processing Systems, vol. 19, http://books.nips.cc/papers/files/nips19/NIPS2006_0202.pdf, 2006. S. Araki, H. Sawada, R. Mukai and S. Makino, "Underdetermined Blind Sparse Source Separation for Arbitrarily Arranged Multiple Sensors," Signal Processing., vol. 87, no.8, 99. 1833-1847, 2007. S. Winter, W. Kellermann, H. Sawada, and S. Makino, "MAP-Based Underdetermined Blind Source Separation of Convolutive Mixtures by Hierarchical Clustering and L1-Norm Minimization," EURASIP Journal on Advances in Signal Processing, 2007, Article ID 24717. H. Sawada, R. Mukai, S. Araki, S. Makino, " A robust and precise method for solving the permutation problem of frequency-domain blind source separation," IEEE Trans. Speech and Audio Processing, vol. 12, no. 5, pp. 530-538, Sep. 2004. R.K. Olsson and L.K. Hansen, "Blind Separation of More Sources than Sensors in Convolutive Mixtures," Proc. ICASSP 2006, May 2006, vol. V, pp. 657-660. J. Anemuller, B. Kollmeier, "Amplitude Modulation Decorrelation for Convolutive Blind Source Separation," in Proc. ICA 2000, June 2000, pp. 215-220. N. Murata, S. Ikeda, and A. Ziehe, "An Approach to Blind Source Separation Based on Temporal Structure of Speech Signals," Neurocomputing, vol. 41, pp. 1-24, Oct. 2001.
しかし、すべての時間周波数(f,τ)に関する観測信号ベクトルX(f,τ)を一気にクラス分類したり事後確率計算したりする方法の場合、反射や残響の影響が強い環境では信号源の方向や位置が正確に推定できず信号分離性能が劣化してしまう。
また、観測信号ベクトルX(f,τ)をN個のクラスヘ分類する操作を周波数毎に行い、各周波数での分類結果から各信号源の方向や位置に相当する値を推定し、その推定結果に基づいてパーミュテーション問題を解決する場合にも、反射や残響の影響が強い環境では信号源の方向や位置が正確に推定できず信号分離性能が劣化してしまう。
一方、観測信号ベクトルX(f,τ)をN個のクラスヘ分類する操作を周波数毎に行い、周波数毎の分離信号エンベロープの相関係数の類似度に基づいてパーミュテーション問題を解決する場合には、反射や残響の影響をそれほど受けることなく信号分離を行うことができる。しかし、この従来方法の場合、源信号が全周波数に渡って同じような振幅のエンベロープを持たない限り、周波数全体に渡って一貫性のあるパーミュテーション問題の解を得ることはできない。以下、このことを詳細に説明する。
この従来方法では、式(6)などに従って分離信号Yn(f,τ)を周波数f毎に計算した後、それらのエンベロープをvi f(τ)=|Yn(f,τ)|として計算する。そして、それらの類似度を相関係数によって表現し、同一の源信号に対応するエンベロープ間の相関係数が最も大きくなると仮定してパーミュテーション問題を解決する。なお、系列長（時間インデックスτの数）がそれぞれTである２つの系列vi f(τ)，vk g(τ)の相関係数は、
として計算される。ここで、
は、それぞれ、相関、平均、標準偏差である。また、相関係数は、−１から１までの値を取り、２つの系列が等しいときには１となる。また、
は、T個のτにそれぞれ対応するα(τ)の和を意味する。
図１８（ａ）に、信号源が３つ存在する場合の２つの周波数f=766Hz,g=906Hzにおける分離信号のエンベロープ系列v1 f(τ),...,v3 f(τ)，v1 g(τ),...,v3 g(τ)を例示する。なお、図１８（ａ）の横軸は時間（時間インデックスτ）を示し、縦軸はエンベロープを示す。また、エンベロープ系列v1 f(τ),...,v3 f(τ)，v1 g(τ),...,v3 g(τ)は、同じ信号源に対応するものが同じ添字となるようにパーミュテーションが揃えられている（例えば、v1 f(τ)とv1 g(τ)とは同じ信号源に対応する系列である）。そして、３つの信号源にそれぞれ対応する系列を濃い実線、薄い実線、破線で区別してある。
ここで、図１８（ａ）のエンベロープ系列間の相関係数を求めると以下のようになる。
一般に、周波数fとgが隣接や倍音の関係にあれば、同じ信号源に対応する分離信号のエンベロープ系列の相関係数の値は、異なる信号源に対応する分離信号のエンベロープ系列の相関係数の値よりも格段に大きくなる。しかし、式(8)の例では、周波数fとgが隣接や倍音の関係になっていないため、同じ信号源に対応する分離信号のエンベロープ系列の相関係数の値が、異なる信号源に対応する分離信号のエンベロープ系列の相関係数の値よりもさほど大きくなっていない。特に、相関係数ρ(v1 f，v1 g)とρ(v1 f，v2 g)との間では大小関係が逆転している。これでは周波数全体に渡って一貫性のあるパーミュテーション問題の解を得ることはできない。
本発明はこのような点に鑑みてなされたものであり、反響や残響の影響が強い環境であっても、また、信号源の全周波数に渡って同じような振幅のエンベロープを持たない場合であっても、高精度にブラインド信号分離を行うことができる技術を提供することを目的とする。
本発明では、まず、周波数領域変換部が、源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時間）を、周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換する。
次に、クラス分類部が、周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tを周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f)(n=1,...,N、N≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出する。この処理は周波数f毎に独立に行われるため、反響や残響の影響が強い場合であっても処理の精度はさほど低下しない。また、当該クラス分類部で算出された事後確率P(Cn(f)|X(f,τ))に対応するクラスCn(f)の番号nは信号源に対応する。しかし、クラスCn(f)の番号nと信号源との対応関係は周波数f毎に相違する可能性が高い。
次に、パーミュテーション問題解決部が、対応する周波数fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成する。ここで、同じ信号源から発せられた信号が支配的な観測信号ベクトルX(f,τ)に対応する事後確率P(Cn(f)|X(f,τ))は、周波数fが相違する場合であっても類似度が大きい。この性質は全周波数中の多くの組合せについて妥当なものである。パーミュテーション問題解決部は、この性質を利用してパーミュテーション問題を解決する。
そして、分離部が、パーミュテーション問題解決部で生成された事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yn(f,τ)として抽出する。
本発明では、反響や残響の影響が強い環境であっても、また、信号源の全周波数に渡って同じような振幅にエンベロープを持たない場合であっても、高精度にブラインド信号分離を行うことができる。
以下、本発明を実施するための最良の形態を図面を参照して説明する。
〔原理〕
まず、本形態の信号分離の原理について説明する。
本形態では、周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]TをN個のクラスへ分類（クラスタリング）する操作を周波数f毎に行う。従ってパーミュテーション問題を解決する必要がある。本形態では、この分類操作において観測信号ベクトルX(f,τ)がクラスCn(f)に属する事後確率P(Cn(f)|X(f,τ))を明示的に計算しておき、この事後確率の系列を用いてパーミュテーション問題を解決する。このように周波数毎の分類操作において事後確率を明示的に計算し、その事後確率の系列を用いてパーミュテーション問題を高精度に解決する部分に本形態の特徴がある。
すなわち、本形態では、事後確率P(Cn(f)|X(f,τ))の系列（「アクティブ系列」と呼ぶ）
vn f(τ)= P(Cn(f)|X(f,τ)) …(9)
の類似度を求め、同一の源信号に対応するアクティブ系列間の類似度が最も大きくなると仮定してパーミュテーション問題を解決する。
図１８（ｂ）は、前述の図１８（ａ）と同じ観測信号に基づき生成されたアクティブ系列vn f(τ)= P(Cn(f)|X(f,τ))を示したグラフである。なお、図１８（ｂ）の横軸は時間（時間インデックスτ）を示し、縦軸は事後確率を示す。また、アクティブ系列v1 f(τ),...,v3 f(τ)，v1 g(τ),...,v3 g(τ)は、同じ信号源に対応するものが同じ添字となるようにパーミュテーションが揃えられている。そして、３つの信号源にそれぞれ対応する系列を濃い実線、薄い実線、破線で区別してある。
ここで、図１８（ｂ）と前述の図１８（ａ）と比較すれば分かるように、同じ信号源に対応するアクティブ系列の相関関係は、同じ信号源に対応するエンベロープ系列の相関関係よりも強いことが分かる。
例えば、図１８（ａ）に例示した周波数fのエンベロープ系列の時間３秒直後では、破線で示されたエンベロープ系列の振幅が格段に大きくなっており、濃い実線や薄い実線で示されたエンベロープ系列の振幅は零に近い。しかし、このようなエンベロープ系列の特徴は、図１８（ａ）の周波数ｇのエンベロープ系列の時間３秒直後には顕著に表れていない。これに対し、図１８（ｂ）に示す周波数fのアクティブ系列の時間３秒直後と、周波数ｇのアクティブ系列の時間３秒直後とは互いに類似した振幅を持つ。
また、図１８（ｂ）に示すアクティブ系列について式(7)に示した相関係数を求めると
となる。このように、同じ信号源に対応するアクティブ系列の相関係数の値は、異なる信号源に対応するアクティブ系列の相関係数の値よりも格段に大きくなる。
なお、上記の例では二つの周波数f=766Hz及びg=906Hzの組み合わせを選択して、それらに関する相関係数を計算した結果を示した。しかし、同じ信号源に対応する系列間の相関係数の値が大きくなるという傾向は、多くの周波数の組み合わせに関し、エンベロープ系列よりもアクティブ系列のほうが顕著である。
以上より、本形態では、事後確率の系列であるアクティブ系列の類似度を用いることで、エンベロープ系列を用いていた従来技術よりも高精度にパーミュテーション問題を解決することができる。
〔第１実施形態〕
＜信号分離装置の構成＞
図１は、本形態の信号分離装置１０の機能構成の全体を例示したブロック図である。また、図２（ａ）は、図１に示したクラス分類部１２０の機能構成の詳細を例示したブロック図である。また、図２（ｂ）は、図１に示したパーミュテーション問題解決部１３０の機能構成の詳細を例示したブロック図である。また、図３は、図２（ｂ）のクラスタリング部１３２の機能構成の詳細を例示したブロック図である。また、図４は、本形態の信号分離装置１０を構成するハードウェアの構成を例示したブロック図である。なお、各図において、実線の矢印はデータの流れを示し、点線の矢印は論理的な情報の流れを示す。しかし、制御部１６０やメモリ１７０等、一部の構成に対するデータの流れの表記は省略する。
以下、これらの図を用い、本形態の信号分離装置の構成を説明する。
［ハードウェア構成］
図４に例示するように、この例の信号分離装置１０は、ＣＰＵ（Central Processing Unit）１０ａ、入力部１０ｂ、出力部１０ｃ、補助記憶装置１０ｆ、ＲＡＭ（Random Access Memory）１０ｄ、ＲＯＭ（Read Only Memory）１０ｅ及びバス１０ｇを有している。
この例のＣＰＵ１０ａは、制御部１０ａａ、演算部１０ａｂ及びレジスタ１０ａｃ有し、レジスタ１０ａｃに読み込まれた各種プログラムに従って様々な演算処理を実行する。また、補助記憶装置１０ｆは、本形態の信号分離処理を実行するための信号分離プログラムを格納した信号分離プログラム領域１０ｆａ及びセンサで観測された時間領域の混合信号等の各種データが格納されるデータ領域１０ｆｂを有している。また、ＲＡＭ１０ｄは、信号分離プログラムが書き込まれる信号分離プログラム領域１０ｄａ及び各種データが書き込まれるデータ領域１０ｄｂを有している。また、この例のバス１０ｇは、ＣＰＵ１０ａ、入力部１０ｂ、出力部１０ｃ、補助記憶装置１０ｆ、ＲＡＭ１０ｄ及びＲＯＭ１０ｅを通信可能に接続している。
［ハードウェアとソフトウェアとの協働］
この例のＣＰＵ１０ａは、読み込まれたＯＳ（Operating System）プログラムに従い、補助記憶装置１０ｆの信号分離プログラム領域１０ｆａに格納されている信号分離プログラムを、ＲＡＭ１０ｄの信号分離プログラム領域１０ｄａに書き込む。同様にＣＰＵ１０ａは、補助記憶装置１０ｆのデータ領域１０ｆｂに格納されている時間領域の混合信号等の各種データをＲＡＭ１０ｄのデータ領域１０ｄｂに書き込む。さらに、ＣＰＵ１０ａは、この信号分離プログラムや各種データが書き込まれたＲＡＭ１０ｄ上のアドレスをレジスタ１０ａｃに格納する。そして、ＣＰＵ１０ａの制御部１０ａａは、レジスタ１０ａｃに格納されたこれらのアドレスを順次読み出し、読み出したアドレスが示すＲＡＭ１０ｄ上の領域からプログラムやデータを読み出し、そのプログラムが示す演算を演算部１０ａｂに順次実行させ、その演算結果をレジスタ１０ａｃに格納していく。
このようにＣＰＵ１０ａに信号分離プログラムが読み込まれることにより図１から図３に例示する機能構成を具備する信号分離装置１０が構築される。
図１に例示するように、本形態の信号分離装置１０は、メモリ１００，１７０と、周波数領域変換部１１０と、クラス分類部１２０と、パーミュテーション問題解決部１３０と、分離部１４０と、時間領域変換部１５０と、信号分離装置１０全体を制御する制御部１６０とを有する。また、図２（ａ）に例示するように、この例のクラス分類部１２０は、ノルム正規化部１２１と、モデル化部１２２とを有し、モデル化部１２２は、初期パラメータ設定部１２２ａと、事後確率計算部１２２ｂと、パラメータ推定部１２２ｃと、演算制御部１２２ｄとを有する。また、図２（ｂ）に例示するように、パーミュテーション問題解決部１３０は、アクティブ系列生成部１３１と、クラスタリング部１３２と、並び替え部１３３とを有する。また、図３に例示するように、この例のクラスタリング部１３２は、大域的最適化部１３２ａと、局所的最適化部１３２ｂとを有する。また、大域的最適化部１３２ａは、初期パラメータ設定部１３２ａａと、順列生成部１３２ａｂと、セントロイド算出部１３２ａｃと、演算制御部１３２ａｄとを有し、局所的最適化部１３２ｂは、順列生成部１３２ｂａと、演算制御部１３２ｂｂとを有する。
ここでメモリ１００，１７０は、レジスタ１０ａｃ、補助記憶装置１０ｆのデータ領域１０ｆｂ或いはＲＡＭ１０ｄのデータ領域１０ｄｂ等に相当する。また、周波数領域変換部１１０、クラス分類部１２０、パーミュテーション問題解決部１３０、分離部１４０、時間領域変換部１５０及び制御部１６０は、ＣＰＵ１０ａにＯＳプログラムや信号分離プログラムが読み込まれることにより構成されるものである。
＜信号分離方法＞
図５は、本形態の信号分離方法の全体を説明するためのフローチャートである。また、図６は、図５のクラス分類過程（ステップＳ２）の詳細を説明するためのフローチャートである。また、図７（ａ）は、図５のパーミュテーション問題解決過程（ステップＳ３）の詳細を説明するためのフローチャートであり、図７（ｂ）は、図７（ａ）のクラスタリング過程（ステップＳ２２）の詳細を説明するためのフローチャートである。また、図８（ａ）は、図７（ｂ）の大域的最適化過程（ステップＳ３１）の詳細を説明するためのフローチャートであり、図８（ｂ）は、図７（ｂ）の局所的最適化過程（ステップＳ３２）の詳細を説明するためのフローチャートである。以下、これらの図を用い、本形態の信号分離方法を説明する。なお、各演算は、制御部１６０の制御のもと実行される。また、明示しないか限り、各演算過程で生成された演算結果は逐一メモリ１７０に格納され、必要に応じて読み出されて他の演算に用いられる。
［信号分離方法の全体］
まず、源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時刻）が信号分離装置１０（図１）のメモリ１００の記憶領域１０１に格納される。なお、観測信号xm(t)は、サンプリング周波数fsでサンプリングされた離散値である。
次に、周波数領域変換部１１０が、記憶領域１０１から観測信号xm(t)を読み込み、それらを周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換し、生成された周波数領域の観測信号Xm(f,τ)をｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０２に格納する（周波数領域変換過程／ステップＳ１）。
次に、クラス分類部１２０が、記憶領域１０２から周波数領域の観測信号Xm(f,τ)を読み込み、これらを要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tをスパース性に基づき周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f) (n=1,...,N、N≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出する。算出された事後確率P(Cn(f)|X(f,τ))は、ｎ，ｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０３に格納される（クラス分類過程／ステップＳ２）。
次に、パーミュテーション問題解決部１３０が、記憶領域１０３から事後確率P(Cn(f)|X(f,τ))を読み込み、同じ信号源に対応するクラスがすべての周波数に渡って同じ添字（クラスCn(f)の番号n）を持つように、クラスCn(f)の番号nを付け替える。すなわち、パーミュテーション問題解決部１３０は、対応する周波数fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成する。生成された事後確率P’(Ck(f)|X(f,τ))は、ｋ，ｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０４に格納される（パーミュテーション問題解決過程／ステップＳ３）。
次に、分離部１４０が、記憶領域１０２から周波数領域の観測信号Xm(f,τ)を読み込み、記憶領域１０４から事後確率P’(Ck(f)|X(f,τ))を読み込み、それを用い、周波数領域の分離信号Yn(f,τ)を抽出する。すなわち、分離部１４０が、事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yk(f,τ)として抽出する。抽出された分離信号Yk(f,τ)は、ｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０５に格納される（分離過程／ステップＳ４）。
最後に、時間領域変換部１５０が、記憶領域１０５から分離信号Yk(f,τ)を読み込み、分離信号Yk(f,τ)を時間領域の分離信号yk(t)に変換してメモリ１００の記憶領域１０６に格納する（時間領域変換過程／ステップＳ５）。
［周波数領域変換過程（ステップＳ１）の詳細］
周波数領域変換過程（ステップＳ１）は、例えば、式(2)に従った短時間フーリエ変換によって行う。
［クラス分類過程（ステップＳ２）の詳細］
本形態のクラス分類過程では、クラスCn(f)の代表ベクトルであるセントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率のモデルP(Cn(f)|X(f,τ), θ(f))（θ(f)はパラメータ集合）を生成し、事後確率P(Cn(f)|X(f,τ), θ(f))の計算とパラメータ集合θ(f)の推定とを、所定の終了条件を満たすまで交互に繰り返し、事後確率P(Cn(f)|X(f,τ))を求める。
以下では、混合ガウス分布で観測信号ベクトルX(f,τ)の集合をモデル化する方法を例示する。まず、ガウス分布の混合数を、仮定した源信号の数Nとし、観測信号ベクトルX(f,τ)がクラスCn(f)に属する事象の確率密度関数を、ガウス分布
でモデル化する。なお、an(f)は、ｎ番目のクラスCn(f)に属する観測信号ベクトルX(f,τ)のセントロイドである。本形態では、an(f)のノルムが所定値（例えば１）に正規化されている。また、σn(f)は、ｎ番目のクラスCn(f)に属する観測信号ベクトルX(f,τ)の標準偏差であり、(σn(f))2は、ｎ番目のクラスCn(f)に属する観測信号ベクトルX(f,τ)の分散である。また、・Hは・の複素共役転置である。また、‖・‖は・のノルムを示す。
ここで、{ (an(f))H・X(f,τ)}・an(f)は、セントロイドan(f)が張る部分空間への観測信号ベクトルX(f,τ)の直交射影である。そのため、式(11)の‖X(f,τ)-{ (an(f))H・X(f,τ)}・an(f)‖は、観測信号ベクトルX(f,τ)とセントロイドan(f)が張る部分空間との最少距離を示し、これが小さいほど観測信号ベクトルX(f,τ)がクラスCn(f)に属する尤度が高くなる。
次に、各ガウス分布p(X(f,τ)|an(f),σn(f))の混合比をαn(f)（0<αn(f)<1とα1+…+αＮ =1とを満たす）とし、パラメータ集合をθ(f)={a1(f),σ1(f),α1(f),...,aN(f),σN(f),αN(f)}とすると、混合ガウス分布による密度関数は、
p(X(f,τ)|θ(f))=Σn=1 Nαn(f)・p(X(f,τ)|an(f),σn(f)) …(12)
と表現される。
ここで、周波数f毎にＴ個の観測信号ベクトルX(f,τ)のサンプルが得られたとする。最尤推定の原理では、対数尤度の和
Στ Tlog p(X(f,τ)|θ(f))=Στ TlogΣn=1 Nαn(f)・p(X(f,τ)|an(f),σn(f))
を最大化するパラメータ集合θ(f)を求めるが、この形では、対数の中に確率密度関数p(X(f,τ)|an(f),σn(f))の和が含まれているため、計算が困難となる。そこで、本形態では、EMアルゴリズム（例えば、「汪金芳,田栗正章,手塚集,樺島祥介,上田修功,「計算統計I確率計算の新しい手法」,統計科学のフロンティア11,ISBN4-00-006851-2」等参照）を用いてパラメータ推定を行う。EMアルゴリズムでは、対数尤度の代わりにいわゆるQ関数
Q(f,θ(f))=Στ TΣn=1 N {P(Cn(f)|X(f,τ),θ(f))・logαn(f)・p(X(f,τ)|an(f),σn(f))} …(13)
を最大化するパラメータ集合θ(f)を求める。ここで、P(Cn(f)|X(f,τ),θ(f))は、観測信号ベクトルX(f,τ)を周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f)となる事象の、観測信号ベクトルX(f,τ)を得た後における事後確率であり、ベイズの定理により、
P(Cn(f)|X(f,τ),θ(f))=αn(f)・p(X(f,τ)|an(f),σn(f))/p(X(f,τ)|θ(f))
…(14)
と書き下せる。
すなわち、本形態のクラス分類過程では、現在のパラメータ集合θ(f)を固定したまま、式 (11)(12)(14)を用いて、すべての観測信号ベクトルX(f,τ)とクラスCn(f)に対し、周波数f毎に事後確率P(Cn(f)|X(f,τ),θ(f))を計算する事後確率計算過程（E-step）と、事後確率P(Cn(f)|X(f,τ),θ(f))を固定したまま、式(13)のQ関数Q(f,θ(f))が最大となるパラメータ集合θ(f)を計算するパラメータ推定過程（M-step）とを、所定の終了条件を満たすまで繰り返し、終了条件を満たした時点の事後確率P(Cn(f)|X(f,τ),θ(f))を事後確率P(Cn(f)|X(f,τ))として出力する。
なお、この例のパラメータ推定過程でのセントロイドan(f)は、クラスCn(f)に関する相関行列
R=Στ T P(Cn(f)|X(f,τ),θ(f))・X(f,τ)・X H(f,τ)
の最大固有値として算出される。また、分散(σn(f))2は、
として計算される。また、混合比αn(f)は、
として計算される。
また、事後確率計算過程とパラメータ推定過程とは、所定の終了条件を満たすまで繰り返されるが、その最終結果は初期値によって異なったものになることがあるため、初期値の設定は重要である。一般的には、狭い範囲に集中するセントロイドを初期値に設定したり、実際のサンプルからあまりにも乖離した初期値を設定したりすることは避けるべきである。適切にクラスタリングが行われない場合があるからである。好ましい初期値の設定方法には、特に制限はないが、例えば、以下のような方法を例示できる。
まず、セントロイドa1(f),...,aN(f)の初期値には、例えば、観測信号ベクトルX(f,τ)のＴ個のサンプルからランダムに選択したN個のサンプルを用いる。また、分散(σn(f))2の初期値には、例えば、観測信号ベクトルX(f,τ)毎にセントロイドan(f)との２乗距離の最小値（観測信号ベクトルX(f,τ)と何れかのセントロイドan(f)との２乗距離）を求め、それらを平均した値
を用いる。また、混合比αn(f) の初期値は、例えば、αn(f)=1/Nと設定する。
また、上述のようにセントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率をモデル化する場合、その最尤推定に用いるサンプルである観測信号ベクトルX(f,τ)のノルム‖X(f,τ)‖は所定値（例えば１）に正規化されていることが望ましい。本来、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいてクラス分類を行う場合、それらのノルムが各n、f、τにおいて一定値に正規化されていないと、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離を、各n、f、τに渡って厳密に比較評価し、クラス分類を行うことはできない。すなわち、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率をモデル化する場合、サンプルである観測信号ベクトルX(f,τ)のノルムが大きいほどセントロイドan(f)との距離が大きくなり、適切な最尤推定を行うことができない場合がある。これは、事後確率P(Cn(f)|X(f,τ))の推定精度を低下させる。よって、本形態では、尤推定に用いるサンプルである観測信号ベクトルX(f,τ)のノルムを、例えば、
により１に正規化する。なお、α←βは、βの値をαの値とすることを意味する。
以下、本形態のクラス分類過程（ステップＳ２）の詳細を図６のフローチャートに沿って説明する。なお、ここでは周波数fに関する処理のみを説明するが、クラス分類過程（ステップＳ２）の処理は周波数f毎に独立に行われる。
まず、クラス分類部１２０のノルム正規化部１２１（図２）にメモリ１００から読み込まれた観測信号ベクトルX(f,τ)が入力され、ノルム正規化部１２１が当該観測信号ベクトルX(f,τ)のノルムを所定値に正規化する。具体的には、ノルム正規化部１２１は、例えば、式(16)に従って観測信号ベクトルX(f,τ)のノルムを１に正規化する。また、正規化された測信号ベクトルX(f,τ)は、メモリ１７０に格納される（ノルム正規化過程／ステップＳ１１）。
次に初期パラメータ設定部１２２ａが、メモリ１７０からノルムが正規化された観測信号ベクトルX(f,τ)を読み込み、初期パラメータ設定部１２２ａは、これらを用いてパラメータ集合θ(f)={a1(f),σ1(f),α1(f),...,aN(f),σN(f),αN(f)}の初期値を設定する。この初期値の設定は、例えば、先に例示した初期値の設定方法に従って行う。また、設定されたパラメータ集合θ(f)の初期値は、メモリ１７０に格納される（初期パラメータ設定過程／ステップＳ１２）。なお、Nの値については、信号源数が既知なのであればその数をNとして用いてもよいし、厳密な信号源数が未知の場合には、経験則等から仮定できる信号源数をNとして用いてもよい。
次に、事後確率計算部１２２ｂは、メモリ１７０からノルムが正規化された観測信号ベクトルX(f,τ)と最新のパラメータ集合θ(f)とを読み込み、パラメータ集合θ(f)を固定値として、式 (11)(12)(14)を用いて、すべての観測信号ベクトルX(f,τ)とクラスCn(f)に対し、周波数f毎に事後確率P(Cn(f)|X(f,τ),θ(f))を算出する。算出された事後確率P(Cn(f)|X(f,τ),θ(f))は、メモリ１７０に格納される（事後確率計算過程／ステップＳ１３）。
次に、パラメータ推定部１２２ｃが、メモリ１７０から観測信号ベクトルX(f,τ)と最新の事後確率P(Cn(f)|X(f,τ),θ(f))を読み込み、読み込んだ事後確率P(Cn(f)|X(f,τ),θ(f))を固定値として、式(13)に示したQ関数Q(f,θ(f))〔「各観測信号ベクトルX(f,τ)に対応するΣn=1 N {P(Cn(f)|X(f,τ),θ(f))・logαn(f)・p(X(f,τ)|an(f),σn(f))}を周波数f毎に独立に加算したQ(f,θ(f))」に相当〕がそれぞれ最大となるパラメータ集合θ(f)を算出する。この算出方法は、前述した通りである。算出されたパラメータ集合θ(f)は、メモリ１７０に格納される（パラメータ推定過程／ステップＳ１４）。
次に、演算制御部１２２ｄが、所定の終了条件を満たした否かを判定する（終了条件判定過程／ステップＳ１５）。なお「所定の終了条件」としては、例えば、以下を例示できる。
・ステップＳ１３で固定値として用いたパラメータ集合θ(f)からなるベクトルと、ステップＳ１４で新たに算出されたパラメータ集合θ(f)からなるベクトルとの距離が所定値以下（又は未満）であること。
・ステップＳ１３で算出された最新の事後確率P(Cn(f)|X(f,τ),θ(f))と、それよりも１つ前に算出された事後確率P(Cn(f)|X(f,τ),θ(f))（ステップＳ１２又は１つ前のループのステップＳ１３で作成された事後確率）との差の合計が所定値以下（又は未満）であること。
・ステップＳ１３とＳ１４の処理を所定回数繰り返したこと。
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ１３に戻される。一方、所定の終了条件を満たしたと判定された場合、すべてのCn(f)及びX(f,τ)にそれぞれ対応する最新の事後確率P(Cn(f)|X(f,τ),θ(f))がメモリ１７０からクラス分類部１２０に読み込まれ、これらが事後確率P(Cn(f)|X(f,τ))として出力される（事後確率出力過程／ステップＳ１６）。
［パーミュテーション問題解決過程（ステップＳ３）の詳細］
次に、図７（ａ）を用い、パーミュテーション問題解決過程（ステップＳ３）の詳細を説明する。
まず、パーミュテーション問題解決部１３０（図２（ｂ））のアクティブ系列生成部１３１にメモリ１００から読み込まれた各事後確率P(Cn(f)|X(f,τ))が入力される。アクティブ系列生成部１３１は、式(9)に従い、各事後確率P(Cn(f)|X(f,τ))に対応するアクティブ系列vn f(τ)を生成して出力する（アクティブ系列生成過程／ステップＳ２１）。
各アクティブ系列vn f(τ)はクラスタリング部１３２に入力され、クラスタリング部１３２はそれらのクラスタリングを行う。ここでのクラスタリングは一般的なものとは少し異なり、クラスタリング結果は周波数毎の順列Πfとして表現される。より具体的には、クラスタリング部１３２は、対応する周波数fが異なる系列vn f(τ)間の類似度を指標として、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成する。生成された順列Πfはメモリ１７０に格納される（クラスタリング過程／ステップＳ２２）。
次に、並び替え部１３３がメモリ１００から事後確率P(Cn(f)|X(f,τ))を読み込み、メモリ１７０から順列Πfを読み込む。そして、並び替え部１３３は、順列Πfに従い、周波数f毎に事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、事後確率P’(Ck(f)|X(f,τ))を生成する。生成された事後確率P’(Ck(f)|X(f,τ))はメモリ１７０に格納される（並び替え過程／ステップＳ２３）。
［クラスタリング過程（ステップＳ２２）の詳細］
次に、図７（ｂ）を用い、クラスタリング過程（ステップＳ２２）の詳細について説明する。
まず、クラスタリング部１３２の大域的最適化部１３２ａにアクティブ系列vn f(τ)が入力される。大域的最適化部１３２ａは、すべての異なる周波数f∈Fの組合せに対応するアクティブ系列vn f(τ)間の類似度を指標とし、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成する。なお、Fは取り扱う全周波数ビンの集合を意味する。生成された順列Πfはメモリ１７０に格納される（大域的最適化過程／ステップＳ３１）。ここで、前述したように、アクティブ系列vn f(τ)間の類似度は、多くの周波数の組み合わせに関し、同じ信号源に対応する系列間の類似度が大きくなるという傾向をもつ。また、大域的最適化過程では、すべての異なる周波数f∈Fの組合せに対応する系列vn f(τ)間の類似度を指標として順列Πfを生成する。よって、この大域的最適化過程により、周波数全体に渡って一貫性のあるパーミュテーション問題の解となる順列Πfを求めることができる。
次に、局所的最適化部１３２ｂにアクティブ系列vn f(τ)と大域的最適化過程（ステップＳ３１）で生成された順列Πfとが入力される。局所的最適化部１３２ｂは、対応する周波数f∈Fが異なる系列vn f(τ)間の類似度のうち、特定の周波数の組合せに対応する系列間の類似度のみを指標として用い、大域的最適化部で生成された順列Πfを更新し、新たな順列Πfを生成する。生成された順列Πfはメモリ１７０に格納される（局所的最適化過程／ステップＳ３２）。ここで、「特定の周波数の組合せ」として、同じ信号源に対応するアクティブ系列vn f(τ)間の類似度が特に大きくなる周波数の組み合わせを選択することにより、大域的最適化過程で生成された順列Πfを、パーミュテーション問題をより高精度に解決できる順列Πfに補正することができる。なお、「特定の周波数の組合せ」としては、周波数差が所定範囲内にある周波数の組合せや、倍音関係にある周波数の組合せを例示できる。
このように、本形態のクラスタリング過程では、大域的最適化過程（ステップＳ３１）と局所的最適化過程（ステップＳ３２）を順に適用して、各周波数fでの順列Πfを算出することとしたため、周波数全体に渡って一貫性のある高精度な順列Πfを得ることができる。
［大域的最適化過程（ステップＳ３１）の詳細］
次に、大域的最適化過程（ステップＳ３１）の詳細を例示する。
この例では、異なる周波数f∈Fの組合せに対応するアクティブ系列vn f(τ)間の類似度を直接的に指標として用いるのではなく、各信号源に対応するセントロイドcn(τ)とアクティブ系列vn f(τ)との類似度を指標とすることで、異なる周波数f∈Fの組合せに対応するアクティブ系列vn f(τ)間の類似度を間接的に指標として用いる。これにより、大域的最適化過程の演算精度と演算効率が向上する。
具体的には、この例の大域的最適化過程では、信号源毎にセントロイドcn(τ)を推定し、目的関数
を最大化する順列Πfを求める。この目的関数は、順列Πfによって並び替えられたアクティブ系列vk f(τ)｜k=Πf(n)とｎ番目の信号源に対応するセントロイドcn(τ)との類似度をすべての信号源及び全ての周波数で足し合わせたものである。
なお、ρ(vk f,cn)は、アクティブ系列vn f(τ)とセントロイドcn(τ)との類似度を示す関数値であり、例えば、アクティブ系列vn f(τ)とセントロイドcn(τ)との相関係数（式(7)参照）である。しかし、アクティブ系列vn f(τ)とセントロイドcn(τ)との類似度を示すのであれば、別の関数値をρ(vk f,cn)として用いてもよい。例えば、
としてもよいし、その他のアクティブ系列vn f(τ)とセントロイドcn(τ)との距離D(vk f(τ),cn(τ))に対して単調減少の関係にある関数をρ(vk f,cn)としてもよい。また、系列長１のアクティブ系列vn f(τ)とセントロイドcn(τ)との類似度を示す関数値をρ(vk f,cn)として用いてもよい。例えば、
ρ(vk f,cn)= vk f(τ)・cn(τ) …(19)
としてもよい。
式(17)の目的関数は、よく知られたk-means法（例えば、「R.O. Duda, P. E. hart, and D. G. Stork, Pattern Classification, Wiley Interscience, 2nd edition, 2000」等参照）と同じように、セントロイドcn(τ)と順列Πfとを交互に最適化することで最大化することができる。以下、図８（ａ）を用い、本形態の大域的最適化過程（ステップＳ３１）の詳細を例示する。
まず、クラスタリング部１３２（図３）の初期パラメータ設定部１３２ａａにメモリ１７０から読み込まれたアクティブ系列vn f(τ)が入力され、初期パラメータ設定部１３２ａａは、これらのアクティブ系列vn f(τ)を用いてセントロイドcn(τ)の初期値を設定する。初期パラメータ設定部１３２ａａは、例えば、読み込まれたアクティブ系列vn f(τ)から時間インデックスτ毎にN個のサンプルを選択し、それらをN個のセントロイドcn(τ)の初期値とする。生成されたセントロイドcn(τ)の初期値はメモリ１７０に格納される（初期パラメータ設定過程／ステップＳ４１）。
次に、順列生成部１３２ａｂが、メモリ１７０からアクティブ系列vn f(τ)とセントロイドcn(τ)とを読み込み、すべての周波数f∈Fに対応するアクティブ系列vn f(τ)とセントロイドcn(τ)との類似度を指標とし、順列Πfを周波数f毎に生成する。本形態の順列生成部１３２ａｂは、周波数f毎に、アクティブ系列vn f(τ)とセントロイドcn(τ)との類似度を最大化させる順列Πfを以下ように決定する。なお、argmaxΠαは、αを最大にする順列Πを意味する。
そして、このように生成された順列Πfはメモリ１７０に格納される（順列生成過程／ステップＳ４２）。
次に、セントロイド算出部１３２ａｃが、メモリ１７０からアクティブ系列vn f(τ)(n=1,...,N)と最新の順列Πfとを読み込み、順列Πfに従って周波数f毎にアクティブ系列vn f(τ)と番号nとの対応関係を並び替えたアクティブ系列vk f(τ)(k=1,...,N)の周波数方向の平均値又は代表値をセントロイドck(τ)として算出する。なお、並び替えたアクティブ系列vk f(τ)の周波数方向の平均値をセントロイドck(τ)とする場合、セントロイド算出部１３２ａｃは、例えば、
によってセントロイドck(τ)を算出する。ここで、|F|は取り扱う全周波数ビンの集合Fの要素数を意味する。また、並び替えたアクティブ系列vk f(τ)の周波数方向の代表値をセントロイドck(τ)とする場合、セントロイド算出部１３２ａｃは、例えば、kとτとが同じアクティブ系列vk f(τ)の集合毎にアクティブ系列vk f(τ)の振幅の頻度を求め、当該集合内で頻度が最大となる振幅を持つ何れかのサンプルを、k，τに対応するセントロイドck(τ)とする。
そして、このように生成されたセントロイドck(τ)はメモリ１７０に格納される（セントロイド算出過程／ステップＳ４３）。
次に、演算制御部１３２ａｄが、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ４４）。なお「所定の終了条件」としては、例えば、以下を例示できる。
・ステップＳ４３で生成されたセントロイドck(τ)と、前回のループで生成されたセントロイドck(τ)（前回のループで生成されたセントロイドck(τ)が存在しない場合はセントロイドck(τ)の初期値） との距離が所定値以下（又は未満）であること。
・ステップＳ４２で生成された順列Πfと前回のループで生成された順列Πfとの相違箇所が所定個以下（又は未満）であること。
・ステップＳ４２とＳ４３の処理を所定回数繰り返したこと。
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ４２に戻される。一方、所定の終了条件を満たしていると判定された場合、大域的最適化部１３２ａは、メモリ１７０に格納されている最新の順列Πfを出力する（順列出力過程／ステップＳ４５）。
［局所的最適化過程（ステップＳ３２）の詳細］
次に、図８（ｂ）を用い、局所的最適化過程（ステップＳ３２）の詳細を例示する。
まず、局所的最適化部１３２ｂ（図３）の順列生成部１３２ｂaに、順列Πfとアクティブ系列vn f(τ)とが入力される。なお、１回目のループの際に入力される順列Πfは大域的最適化過程（ステップＳ３１）で出力された順列であり、２回目以降のループの際に入力される順列Πfは前回の順列生成過程（ステップＳ４６）で生成された順列である。
順列生成部１３２ｂaは、入力された順列Πfの一部である順列Πg(g∈R(f))を固定し、順列Πgに従ってアクティブ系列vn g(τ)と番号nとの対応関係を並び替えたアクティブ系列vk’ g(τ)と、順列Πに従ってアクティブ 系列vn f(τ)と番号nとの対応関係を並び替えたアクティブ系列vk f(τ)との類似度の全信号源についての総和を最大にする、当該順列Πを、新たな順列Πfとして算出する（順列生成過程／ステップＳ４６）。この処理は各周波数fについて行われ、例えば、以下の式に従って行われる。
なお、R(f)は、周波数fと特定の関係にある周波数gの集合であり、周波数fとgとの組み合わせが「特定の周波数の組合せ」に相当する。すなわち、R(f)は、アクティブ系列vk f(τ)とvk’ g(τ)とが同じ信号源に相当するものであったときに、これらの類似度が特に大きくなるような周波数gを要素に持つ集合であることが望ましい。典型的には、周波数fの近傍の周波数の集合A(f)と、周波数fと倍音関係にある周波数の集合H(f)とに対し、R(f)= A(f)∪H(f)であることが望ましい。また、周波数fの近傍の周波数の集合A(f)は、例えば、
Α(f)={f-3Δf, f-2Δf, f-Δf, f+Δf, f+2Δf, f+3Δf}
として定義できる。ここで、Δf＝(1/L)fsは、隣り合う周波数ビン間の周波数の差である。 また、倍音関係にある周波数の集合Ηは、例えば、
Η(f)={round(f/2)-Δf, round(f/2), round(f/2)+Δf, 2f-Δf, 2f, 2f+Δf}
として定義できる。ここで、round(・)は、周波数の集合Fから・に最も近い周波数を意味する。その他、R(f)=A(f)又はR(f)=H(f)とする構成も可能である。
次に、演算制御部１３２ｂｂが、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ４７）。なお「所定の終了条件」としては、例えば、以下を例示できる。
・ステップＳ４６で更新された順列Πfの更新箇所が所定数以下（又は未満）であること。
・ステップＳ４６の処理を所定回数繰り返したこと。
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ４６に戻される。一方、所定の終了条件を満たしていると判定された場合、局所的最適化部１３２ｂは最新の順列Πfを出力する（順列出力過程／ステップＳ４８）。
［並び替え過程（ステップＳ２３）の詳細］
前述のように、並び替え部１３３は、順列Πfに従い、周波数f毎に事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、事後確率P’(Ck(f)|X(f,τ))を生成する。具体的には、例えば、
P’(Ck(f)|X(f,τ))←P(Cn(f)|X(f,τ))|n=Πf(k), ∀k,f,τ …(24)
に従い、事後確率P’(Ck(f)|X(f,τ))を生成する（パーミュテーション問題解決過程（ステップS３）の説明終わり）。
［分離過程（ステップＳ４）の詳細］
前述のように、分離過程では、分離部１４０（図１）が、事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yk(f,τ)として抽出する。例えば、以下に例示する時間周波数マスキングによる方法が、比較的簡単な分離方法となる。なお、m’=1,...,Mであり、観測信号X m’(f,τ)はセンサm’に対応する観測信号である。
その他、L1ノルム最小化規範やその近似方法による分離方法を用いてもよい（例えば、特許文献１や非特許文献４等参照）。
［時間領域変換過程（ステップＳ５）の詳細］
最後に、時間領域変換部１５０が、分離信号Yk(f,τ)を時間領域の分離信号yk(t)に変換する。この処理は、例えば、短時間逆フーリエ変換等によって行う。
〔第２実施形態〕
次に、本発明の第２実施形態について説明する。
本形態は第１実施形態のクラス分類部及びクラス分類過程の変形例である。以下では、クラス分類部及びクラス分類過程の相違点のみを説明する。
＜クラス分類部及びクラス分類過程の変形例１＞
図９は、第１実施形態のクラス分類部１２０の変形例であるクラス分類部２２０の機能構成を示したブロック図である。
クラス分類部１２０とクラス分類部２２０との構成上の相違点は、クラス分類部１２０がパラメータ集合θ(f)の初期値を設定する初期パラメータ設定部１２２ａを具備するモデル化部１２２を有していたのに対し、クラス分類部２２０が事後確率P(Cn(f)|X(f,τ))の初期値を設定する初期パラメータ設定部２２２ａを具備するモデル化部２２２を有する点である。
また、クラス分類部１２０とクラス分類部２２０との処理上の相違点は、クラス分類部１２０では初期パラメータ設定部１２２ａがステップＳ１２（図６）でパラメータ集合θ(f)の初期値を生成していたのに対し、クラス分類部２２０では初期パラメータ設定部２２２ａが事後確率P(Cn(f)|X(f,τ))の初期値を生成する点と、ステップＳ１２とＳ１３との順序が逆になる点である。
なお、事後確率P(Cn(f)|X(f,τ))の初期値の設定方法については特に制限はないが、広範囲に分布する観測信号ベクトルX(f,τ)が同一のクラスCn(f)に属するような初期値設定は好ましくない。適切なクラス分類がなされない場合があるからである。好ましい初期値の設定方法には、特に制限はないが、例えば、以下の(1)〜(4)のような手順を例示できる。
(1)各観測信号ベクトルX(f,τ)のクラス分けパターンをランダムに複数パターン生成する。
(2)各クラス分けパターンに対し、クラスCn(f)毎のセントロイドの組み合わせを算出する。
(3)セントロイドの組み合わせに毎に、異なるクラスのセントロイド間の内積を求め、内積が最小となる組み合わせを選択する。
(4)選択されたセントロイドの組み合わせに対応するクラス分けパターンに従い、事後確率P(Cn(f)|X(f,τ))の初期値を設定する。例えば、そのクラス分けパターンに従うと観測信号ベクトルX(f,τ)がクラスCn(f)に属することになる場合、P(Cn(f)|X(f,τ))=1とし、P(Cn'(f)|X(f,τ))=0 (n'≠n)として事後確率P(Cn(f)|X(f,τ))の初期値を設定する。
＜クラス分類部及びクラス分類過程の変形例２＞
この変形例では、観測信号ベクトルX(f,τ)のノルムを正規化することなくモデル化部の処理を行う。第１実施形態で述べたように、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率をモデル化する場合、サンプルである観測信号ベクトルX(f,τ)のノルムが正規化されていないと事後確率P(Cn(f)|X(f,τ))の推定精度が低下してしまう。しかし、高い推定精度が要求されない用途に用いる場合や観測信号ベクトルX(f,τ)のノルムが安定している場合などには、クラス分類過程での観測信号ベクトルX(f,τ)のノルムの正規化を省略してもよい。
また、セントロイドan(f)と観測信号ベクトルX(f,τ)の方向のみの類似度に基づいて事後確率をモデル化する場合には、観測信号ベクトルX(f,τ)のノルムを正規化する必要はない。例えば、セントロイドan(f)と観測信号ベクトルX(f,τ)とのコサイン距離
｜XH（f,τ）・an(f)｜/（‖X（f,τ）‖・‖an(f)‖） …(26)
を用い、前述の式(11)の替わりに、
でモデル化してもよい。なお、式(26)の｜α｜はαの絶対値である。
図１０（ａ）は、第１実施形態のクラス分類部１２０（図２（ａ））からノルム正規化部１２１を排除したクラス分類部３２０の機能構成を示したブロック図である。また、図１０（ｂ）は、図９に示したクラス分類部２２０からノルム正規化部１２１を排除したクラス分類部４２０の機能構成を示したブロック図である。このようにクラス分類部がノルム正規化部を具備しない構成であってもよい。
＜クラス分類部及びクラス分類過程の変形例３＞
この変形例では、第１実施形態のモデル化部１２２のように最尤推定によって事後確率P(Cn(f)|X(f,τ))を推定するのではなく、よく知られたK-means法によって観測信号ベクトルX(f,τ)をクラスタリングして各観測信号ベクトルX(f,τ)が属するクラスCn(f)を１つずつ推定する。そして、観測信号ベクトルX(f,τ)がクラスCn(f)に属する場合の事後確率P(Cn(f)|X(f,τ))を１とし、観測信号ベクトルX(f,τ)がクラスCn(f)に属しない場合の事後確率P(Cn(f)|X(f,τ))を０とする。すなわち、この場合の事後確率P(Cn(f)|X(f,τ))は０と１のみをとる。
図１１は、このようなクラス分類部５２０の変形例を示したブロック図であり、図１２は、このような場合のクラス分類過程を説明するためのフローチャートである。なお、これまで説明したのと同じ構成部分については同じ符号を付した。
図１１に示すように、この例のクラス分類部５２０は、ノルム正規化部１２１と、初期パラメータ設定部５２２と、事後確率計算部５２３と、セントロイド算出部５２４と、演算制御部５２５とを有する。以下、この例のクラス分類部５２０によって行われるクラス分類過程を説明する。
まず、ノルム正規化部１２１に観測信号ベクトルX(f,τ)が入力され、ノルム正規化部１２１は、各観測信号ベクトルX(f,τ)のノルムを所定値（例えば１）に正規化した各観測信号ベクトルX (f,τ)を生成してメモリ１７０（図１）に格納する（ノルム正規化過程／ステップＳ６１）。
次に、初期パラメータ設定部６２２がメモリ１７０からノルムが正規化された観測信号ベクトルX(f,τ)を読み込み、これらを用い、周波数f毎に各クラスCi(f)(i=1,...,N)のセントロイドci(f)の初期値を設定してメモリ１７０に格納する（初期パラメータ設定過程／ステップＳ６２）。例えば、これらの観測信号ベクトルX(f,τ)からランダムにN個のサンプルを選択し、それらをセントロイドci(f)の初期値とする。
次に、事後確率計算部５２３が、メモリ１７０からセントロイドci(f)とノルムが正規化された観測信号ベクトルX(f,τ)とを読み込み、各セントロイドci(f)を固定値として、観測信号ベクトルX(f,τ)毎に、 観測信号ベクトルX(f,τ)との距離が最も近いセントロイドcB(f)(B⊂{1,...,N})を選択し、
P(Cn(f)|X(f,τ))=1 (if n=B)
P(Cn(f)|X(f,τ))=0 (if n≠B)
として、事後確率を算出してメモリ１７０に格納する（事後確率計算過程／ステップＳ６３）。
次に、セントロイド算出部５２４が、メモリ１７０から事後確率P(Cn(f)|X(f,τ))と正規化された観測信号ベクトルX(f,τ)とを読み込み、読み込んだ事後確率P(Cn(f)|X(f,τ))を固定値として、各クラスCn(f)のセントロイドcn(f)を算出してメモリ１７０に格納する（セントロイド算出過程／ステップＳ６４）。例えば、クラスCn(f)に関する相関行列
R=Στ T P(Cn(f)|X(f,τ))・X(f,τ)・X H(f,τ)
の最大固有値としてセントロイドcn(f)が算出される。
次に、演算制御部５２５が、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ６５）。なお「所定の終了条件」としては、例えば、以下を例示できる。
・ステップＳ６４で生成されたセントロイドcn(f)と、前回のループで生成されたセントロイドcn(f)（前回のループで生成されたセントロイドcn(f)が存在しない場合はセントロイドcn(f)の初期値） との距離が所定値以下（又は未満）であること。
・ステップＳ６３とＳ６４の処理を所定回数繰り返したこと。
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ６３に戻される。一方、所定の終了条件を満たしたと判定された場合、すべてのCn(f)及びX(f,τ)にそれぞれ対応する最新の事後確率P(Cn(f)|X(f,τ))が出力される（事後確率出力過程／ステップＳ６６）。
また、図１３（ａ）のクラス分類部６２０のように、初期パラメータ設定部６２２が、正規化された観測信号ベクトルX(f,τ)を用いて事後確率P(Cn(f)|X(f,τ))の初期値を生成し、セントロイド算出部５２４がステップＳ６４と同様に各クラスCn(f)のセントロイドcn(f)を算出し、事後確率計算部５２３がステップＳ６３と同様に事後確率P(Cn(f)|X(f,τ))を算出し、所定の終了条件を満たすまでセントロイド算出部５２４と事後確率計算部５２３との処理を繰り返すこととしてもよい。
また、図１３（ｂ）のクラス分類部７２０のように、図１１のクラス分類部５２０からノルム正規化部１２１を排除し、観測信号ベクトルX(f,τ)のノルムを正規化せずに事後確率P(Cn(f)|X(f,τ))を算出してもよい。また、図１３（ｃ）のクラス分類部８２０のように、図１３（ａ）のクラス分類部６２０からノルム正規化部１２１を排除し、観測信号ベクトルX(f,τ)のノルムを正規化せずに事後確率P(Cn(f)|X(f,τ))を算出してもよい。
この場合、事後確率計算部５２３は、観測信号ベクトルX(f,τ)との距離が最も近いセントロイドcB(f)(B⊂{1,...,N})を選択するのではなく、例えば、測信号ベクトルX(f,τ)毎に、観測信号ベクトルX(f,τ)とのコサイン距離
cosθ=｜XH（f,τ）・cB(f)｜/（‖X（f,τ）‖・‖cB(f)‖） …(28)
が最も近いセントロイドcB(f)を選択し、
P(Cn(f)|X(f,τ))=1 (if n=B)
P(Cn(f)|X(f,τ))=0 (if n≠B)
として、事後確率を算出してメモリ１７０に格納する。
〔第３実施形態〕
次に、本発明の第３実施形態について説明する。
本形態は、第１実施形態のパーミュテーション問題解決部が具備するクラスタリング部及びそのクラスタリング過程の変形例である。以下では、クラスタリング部及びクラスタリング過程の相違点のみを説明する。
＜クラスタリング部及びクラスタリング過程の変形例１＞
図１４（ａ）は、第１実施形態のクラスタリング部１３２（図３）の変形例であるクラスタリング部２３２の機能構成を示したブロック図である。
クラスタリング部１３２とクラスタリング部２３２との相違点は、クラスタリング部１３２が大域的最適化部１３２ａと局所的最適化部１３２ｂを有していたのに対し、クラスタリング部２３２が大域的最適化部１３２ａを有するが局所的最適化部１３２ｂを有しない点である。すなわち、クラスタリング部２３２は、大域的最適化部１３２ａによって算出された順列Πfをクラスタリング部２３２の出力とする。
＜クラスタリング部及びクラスタリング過程の変形例１＞
図１４（ｂ）は、第１実施形態のクラスタリング部１３２（図３）の変形例であるクラスタリング部３３２の機能構成を示したブロック図である。また、図１５は、クラスタリング部３３２の局所的最適化部３３２ｂが行う局所的最適化過程を説明するためのフローチャートである。
クラスタリング部１３２とクラスタリング部３３２との相違点は、クラスタリング部１３２が大域的最適化部１３２ａと局所的最適化部１３２ｂを有していたのに対し、クラスタリング部３３２が局所的最適化部１３２ｂを有するが大域的最適化部１３２ａを有しない点である。すなわち、クラスタリング部３３２は、局所的最適化部１３２ｂのみによって算出された順列Πfをクラスタリング部３３２の出力とする。その相違から、クラスタリング部３３２は、さらに初期順列設定部３３２ｂａを有する。
以下、図１５を用い、クラスタリング部３３２の局所的最適化部１３２ｂが行う局所的最適化過程を説明する。
まず、初期順列設定部３３２ｂａが初期順列Πfを生成してメモリ１７０（図１）に格納する（初期順列設定過程／ステップＳ７１）。次に、局所的最適化部３３２ｂ（図１４（ｂ））の順列生成部１３２ｂaに、順列Πfとアクティブ系列vn f(τ)とが入力される。なお、１回目のループの際に入力される順列Πfは初期順列設定過程（ステップＳ７１）で出力された初期順列Πfであり、２回目以降のループの際に入力される順列Πfは前回の順列生成過程（ステップＳ７２）で生成された順列である。
順列生成部１３２ｂaは、第１実施形態の順列生成過程（ステップＳ４６）と同様に新たな順列Πfとして算出する（順列生成過程／ステップＳ７２）。
次に、演算制御部１３２ｂｂが、第１実施形態の終了条件判定過程（ステップＳ７３）と同様に、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ７３）。
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ７２に戻される。一方、所定の終了条件を満たしていると判定された場合、局所的最適化部３３２ｂは最新の順列Πfを出力する（順列出力過程／ステップＳ７４）。
＜クラスタリング部及びクラスタリング過程の変形例２＞
第１実施形態のクラスタリング部１３２は、各信号源に対応するセントロイドcn(τ)とアクティブ系列vn f(τ)との類似度を指標とすることで、異なる周波数f∈Fの組合せに対応するアクティブ系列vn f(τ)間の類似度を間接的に指標として用い、順列Πfを生成していた。
しかし、異なる周波数f∈Fの組合せに対応するアクティブ系列vn f(τ)間の類似度を直接的に指標として用いて順列Πfを生成してもよい。この場合には、例えば、すべての周波数f∈Fに対応するアクティブ系列vn f(τ)間の類似度を指標として階層的クラスタリングを行い、 そのクラスタリング結果を用いて順列Πfを生成する。以下、このような手法を例示する。
図１６（ａ）（ｂ）は、それぞれ、すべての周波数f∈Fに対応するアクティブ系列vn f(τ)間の類似度を指標として階層的クラスタリングを行い、そのクラスタリング結果を用いて順列Πfを生成する大域的最適化部４３２ａを具備するクラスタリング部４３２，５３２を示したブロック図である。
図１６（ａ）（ｂ）に示すように、ラスタリング部４３２，５３２の相違点は、ラスタリング部４３２がさらに局所的最適化部１３２ｂを具備し、大域的最適化部４３２ａで生成された順列Πfを局所的最適化部１３２ｂで補正して出力するものであるのに対し、ラスタリング部５３２が大域的最適化部４３２ａで生成された順列Πfをそのまま出力するものである点である。また、大域的最適化部４３２ａは、階層的クラスタリング部４３２ａａと、演算制御部４３２ａｂと、クラス選択部４３２ａｃと、順列生成部４３２ａｄとを有する。なお、局所的最適化部１３２ｂについては第１実施形態で説明済であるため、以下では、局域的最適化部４３２ａの説明のみを行う。
図１７は、大域的最適化部４３２ａが行う大域的最適化過程（ステップＳ３１）を説明するためのフローチャートである。以下、この図に従って大域的最適化部４３２ａが行う大域的最適化過程を説明する。
まず、階層的クラスタリング部４３２ａａに各アクティブ系列vn f(τ)が入力され、階層的クラスタリング部４３２ａａは、すべての周波数f∈Fに対応するアクティブ系列vn f(τ)間の類似度を指標として凝集型の階層的クラスタリングを１階層分行い、各アクティブ系列vn f(τ)が属するクラスの情報CLq{vn f(τ)}を生成してメモリ１７０に格納する（階層的クラスタリング過程／ステップＳ８１） 。なお、CLq{vn f(τ)}は、アクティブ系列vn f(τ)がクラスCLqに属することを意味する。また、アクティブ系列vn f(τ)間の類似度としては、アクティブ系列vn f(τ)間の相関係数（式(7)参照）や、式(18)や式(19)のセントロイドcn(τ)をアクティブ系列に置換した関数等を用いることができる。
次に、演算制御部４３２ａｂがメモリ１７０から最新のクラスの情報CLq{vn f(τ)}を読み込み、クラスの総数が仮定された信号源Ｎと定数βとの和以下であるか否かを判定する（終了条件判定過程／ステップＳ８２）。なお、（クラスの総数）≦Ｎを満たすか否かではなく、（クラスの総数）≦Ｎ＋βを満たすか否かを判定する理由は、生成されたクラスがすべて信号源に対応するとは限らないからである。すなわち、クラスタリング精度によっては、何れの信号源にも対応しないクラスが誤って生成される可能性がある。そのため、凝集型の階層的クラスタリングをクラスの総数がＮ以下になるまで実行すると、信号源に対応する正しいクラスが統合され、正しいクラスの数が信号源の数N未満となる可能性があるからである。よって、ここでは、何れの信号源にも対応しないクラスの想定数以上の定数βを設定しておき、信号源に対応する正しいクラスどうしが統合される事態を防止する。
ここで、（クラスの総数）≦Ｎ＋βを満たさないと判定された場合には、処理がステップＳ７１に戻される。一方、（クラスの総数）≦Ｎ＋βを満たすと判定された場合には、クラス選択部４３２ａｃがメモリ１７０から最新のクラスの情報CLq{vn f(τ)}を読み込み、メンバーの多い方から順にN個のクラスを選択し、選択したクラスの情報CL’u{vn f(τ)}(u=1,...,N)をメモリ１７０に格納する（クラス選択過程／ステップＳ８３） 。この判定は、信号源に対応する正しいクラスのメンバー数は、何れの信号源にも対応しない誤ったクラスのメンバー数よりも大きい、という仮定に基づくものである。
次に、順列生成部４３２ａｄがメモリ１７０から、クラス選択部４３２ａｃで選択されたクラスの情報CL’u{vn f(τ)}を読み込み、これらを用い、vn f(τ)をvu f(τ)の値として並び替える順列Πfを周波数f毎に生成してメモリに格納する（順列出力過程／ステップＳ８４） 。
〔実験結果１〕
図１９及び図２０は、音声信号を対象として第１実施形態の信号分離を行った場合のスペクトログラムとクラス分類結果とを示した図である。なお、各図の横軸は時間であり、縦軸は周波数である。また、図１９(a)に示すN=3個の源信号Sn(f,τ)(n=1,2,3)が混ざり合い、その結果、図１９(b)に示すM=2個のセンサでの観測信号Xm(f,τ)(m=1,2)が得られた場合を例示する。
まず、クラス分類部１２０によって周波数f毎にこのような観測信号ベクトルX(f,τ)をN=3個のクラスに分類することにより、図２０（ａ）に示す事後確率P(Cn(f)|X(f,τ))が得られる。なお、図２０（ａ）では、各(f,τ)に対応する事後確率P(Cn(f)|X(f,τ))の大きさ（0≦P(Cn(f)|X(f,τ))≦1）を色の濃度で表現している。色が薄いほど事後確率P(Cn(f)|X(f,τ))が小さく、色が濃いほど事後確率P(Cn(f)|X(f,τ))が大きい。
次に、パーミュテーション解決部１３０により、事後確率とクラスとのの対応関係を並び替えると、図２０（ｂ）に示す事後確率P’(Ck(f)|X(f,τ)) (k=1,2,3)が得られる。図２０（ｂ）より、事後確率のパーミュテーション問題が解決されていることが分かる。
その後、分離部１４０により、観測信号Xm(f,τ)と事後確率P’(Ck(f)|X(f,τ))とを用いて分離信号Yk(f,τ)を生成すると、図２０（ｃ）の結果が得られる。
〔実験結果２〕
次に、第１実施形態の効果を示すために、図２１（ａ）に示す実験条件と図２１（ｂ）に示す３つのマイクロホンと４つのスピーカの配置を用いて実験を行った。４つの音を同時に鳴らした時の混合音を３つのマイクロホンで観測し、その観測信号のみからそれぞれの音に対応する分離信号を算出するという問題設定である。様々な音声信号の組合せで評価できるように、スピーカからマイクロホンまでのインパルス応答を測定し、音声信号をインパルス応答に畳み込んで混合することで観測信号を生成した。分離性能は、signal-to-interference ratio(SIR)の改善量で評価した。これは、各出力i毎に、出力SIRと入力SIRの差OutputSIRi-InputSIRiとして計算される。入力SIRと出力SIRは、それぞれ以下の式で計算される。
ここで、J∈｛1,...,M｝はある選択された基準センサの番号を示す。また、源信号skのうち出力yi(t)に出てきた成分をyik(t)と表記する。この定義により、yi(t)=Σk=1 Nyik(t)が満たされる。
４つの音声の組合せを８通り用いて実験を行った。また、スピーカからマイクロホンまでの２種類の距離（60, 120cm）と６種類の残響時間（130, 200, 270, 320, 380, 450ms）を試した。
図２１（ｃ）に、４種類のパーミュテーション解決法を用いた結果を、SIR改善量の全出力に関する平均値で示す。“TDOA”と“Envelope”は,それぞれ、従来の技術に相当するものであり、信号の方向や位置に相当する値（センサ間到達時間差、TDOA: Time Difference Of Arreival）を推定することに基づくものと、分離信号エンベロープの相関係数に基づくものである。“Posterior”は、第１実施形態の事後確率の系列を用いるものである。“Optima1”と書かれたものは、信号源に関する情報を用いて最適なパーミュテーションを算出したものである。現実的な状況では、そのような情報は得られないが、性能の上限を示す目的で掲載した。
図２１（ｃ）の結果から以下のことが考察できる。“TDOA”は、スピーカからマイクロホンまでの距離が短い（60cm）場合や残響時間が短い場合（130ms）には適度に良い結果を出しているが、スピーカからマイクロホンまでの距離が長く(120cm)、残響の影響が大きい場合には性能が劣化している。“Envelope”は、多くの場合、それほど良い結果にはなっていない。“Posterior”は、“Optima1”以外の現実的な方法の中で最も良い性能を達成しており、第１実施形態の効果が確認できる。
〔その他の変形例等〕
なお、本発明は上述の各実施形態に限定されるものではない。例えば、上述の各種初期パラメータの設定（ステップＳ４１等）には観測信号ベクトルを用いることとしたが、初期パラメータを固定値とする構成であってもよい。
また、上述の実施形態では、得られた周波数領域の分離信号を時間領域に変換することとしたが、得られた周波数領域の分離信号をそのまま出力する構成であってもよい。
さらに、上述の実施形態では、時間領域と周波数領域との変換に短時間フーリエ変換を用いることとしたが、wavelet変換、DFTフィルタバンク、ポリフェイズフィルタバンクなどを用い、この変換を行うこととしてもよい（例えば、「R. E. Crochiere, L. R. Rabiner, "Multirate Digital Signal Processing." Eaglewood Cliffs, NJ: Prentice-Hall,1983 （ISBN 0-13-605162-6）」参照）。
また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（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）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
本技術により、様々な妨害信号が発生する実環境において、目的の信号を精度良く取り出すことが可能となる。音信号に対する応用例としては、音声認識器のフロントエンドとして働く音源分離システムなどが挙げられる。話者とマイクが離れた位置にあり、マイクが話者の音声以外を収音してしまうような状況でも、そのようなシステムを使うことで、話者の音声のみを取り出して正しく音声を認識することができる。
図１は、第１実施形態の信号分離装置の機能構成の全体を例示したブロック図である。 図２（ａ）は、図１に示したクラス分類部の機能構成の詳細を例示したブロック図である。また、図２（ｂ）は、図１に示したパーミュテーション問題解決部の機能構成の詳細を例示したブロック図である。 図３は、図２（ｂ）のクラスタリング部の機能構成の詳細を例示したブロック図である。 図４は、第１実施形態の信号分離装置を構成するハードウェアの構成を例示したブロック図である。 図５は、第１実施形態の信号分離方法の全体を説明するためのフローチャートである。 図６は、図５のクラス分類過程（ステップＳ２）の詳細を説明するためのフローチャートである。 図７（ａ）は、図５のパーミュテーション問題解決過程（ステップＳ３）の詳細を説明するためのフローチャートである。また、図７（ｂ）は、図７（ａ）のクラスタリング過程（ステップＳ２２）の詳細を説明するためのフローチャートである。 図８（ａ）は、図７（ｂ）の大域的最適化過程（ステップＳ３１）の詳細を説明するためのフローチャートである。また、図８（ｂ）は、図７（ｂ）の局所的最適化過程（ステップＳ３２）の詳細を説明するためのフローチャートである。 図９は、第１実施形態のクラス分類部の変形例であるクラス分類部の機能構成を示したブロック図である。 図１０（ａ）（ｂ）は、第１実施形態のクラス分類部（図２（ａ））からノルム正規化部を排除したクラス分類部の機能構成を示したブロック図である。 図１１は、クラス分類部の変形例を示したブロック図である。 図１２は、クラス分類過程の変形例を説明するためのフローチャートである。 図１３（ａ）（ｂ）（ｃ）は、クラス分類部の変形例を示したブロック図である。 図１４（ａ）（ｂ）は、第１実施形態のクラスタリング部（図３）の変形例を示したブロック図である。 図１５は、局所的最適化過程の変形例を説明するためのフローチャートである。 図１６（ａ）（ｂ）は、それぞれ、すべての周波数f∈Fに対応するアクティブ系列vn f(τ)間の類似度を指標として階層的クラスタリングを行い、そのクラスタリング結果を用いて順列Πfを生成する大域的最適化部を具備するクラスタリング部を示したブロック図である。 図１７は、大域的最適化過程（ステップＳ３１）の変形例を説明するためのフローチャートである。 図１８（ａ）は、信号源が３つ存在する場合の２つの周波数f=766Hz,g=906Hzにおける分離信号のエンベロープ系列v1 f(τ),...,v3 f(τ)，v1 g(τ),...,v3 g(τ)を例示したグラフである。図１８（ｂ）は、図１８（ａ）と同じ観測信号に基づき生成されたアクティブ系列vn f(τ)= P(Cn(f)|X(f,τ))を示したグラフである。 図１９は、音声信号を対象として第１実施形態の信号分離を行った場合のスペクトログラムとクラス分類結果とを示した図である。 図２０は、音声信号を対象として第１実施形態の信号分離を行った場合のスペクトログラムとクラス分類結果とを示した図である。 図２１（ａ）は実験条件を示した表であり、図２１（ｂ）は、実験に用いた３つのマイクロホンと４つのスピーカの配置を示した図である。また、図２１（ｃ）は、実験結果を示したグラフである。
符号の説明
１０ 信号分離装置

Claims (10)

  1. 源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時間）を、周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換する周波数領域変換部と、
    周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tを周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f) (n=1,...,N、N≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出するクラス分類部と、
    対応する周波数fが異なる上記事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成するパーミュテーション問題解決部と、
    上記パーミュテーション問題解決部で生成された上記事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される上記周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yn(f,τ)として抽出する分離部と、
    を有することを特徴とする信号分離装置。
  2. 請求項１に記載の信号分離装置であって、
    上記クラス分類部は、
    クラスCn(f)に属する観測信号ベクトルX(f,τ)の代表ベクトルをセントロイドan(f)とし、クラスCn(f)に属する観測信号ベクトルX(f,τ)の標準偏差をσn(f)とし、観測信号ベクトルX(f,τ)がクラスCn(f)に属する事象の確率密度関数をp(X(f,τ)|an(f),σn(f))とし、p(X(f,τ)|an(f),σn(f))の混合比をαn(f)とし、パラメータ集合θ(f)={a1(f),σ1(f),α1(f),...,aN(f),σN(f),αN(f)}とし、p(X(f,τ)|θ(f))=Σn=1 Nαn(f)・p(X(f,τ)|an(f),σn(f))とした場合における、事後確率P(Cn(f)|X(f,τ),θ(f))=αn(f)・p(X(f,τ)|an(f),σn(f))／p(X(f,τ)|θ(f))を、上記パラメータ集合θ(f)を固定値として各観測信号ベクトルX(f,τ)について算出する事後確率計算部と、
    各観測信号ベクトルX(f,τ)に対応するΣn=1 N {P(Cn(f)|X(f,τ),θ(f))・logαn(f)・p(X(f,τ)|an(f),σn(f))}を周波数f毎に独立に加算したQ(f,θ(f))がそれぞれ最大となるパラメータ集合θ(f)を、事後確率P(Cn(f)|X(f,τ),θ(f))を固定値として算出するパラメータ推定部と、
    所定の終了条件を満たすまで上記事後確率計算部の処理と上記パラメータ推定部の処理とを交互に実行させる第１演算制御部と、を有し、
    上記事後確率P(Cn(f)|X(f,τ))は、
    上記終了条件を満たした際に事後確率計算部で算出されていた最新の事後確率P(Cn(f)|X(f,τ),θ(f))である、
    ことを特徴とする信号分離装置。
  3. 請求項１又は２に記載の信号分離装置であって、
    上記パーミュテーション問題解決部は、
    対応する周波数fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成するクラスタリング部と、
    上記クラスタリング部で生成された順列Πfに従い、周波数f毎に上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替え、上記事後確率P’(Ck(f)|X(f,τ))を生成する並び替え部と、
    を有することを特徴とする信号分離装置。
  4. 請求項３に記載の信号分離装置であって、
    上記クラスタリング部は、
    取り扱う全周波数の集合をFとした場合における、すべての異なる周波数f∈Fの組合せに対応する事後確率P(Cn(f)|X(f,τ))間の類似度を指標とし、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成する大域的最適化部を有する、
    ことを特徴とする信号分離装置。
  5. 請求項４に記載の信号分離装置であって、
    上記大域的最適化部は、
    上記順列Πfに従って周波数f毎に上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替えた事後確率P’(Ck(f)|X(f,τ))の周波数方向の平均値又は代表値をセントロイドck(τ)として算出するセントロイド算出部と、
    上記事後確率P(Cn(f)|X(f,τ))とセントロイドcn(τ)との類似度を指標とし、上記順列Πfを周波数f毎に生成する順列生成部と、
    所定の終了条件を満たすまでセントロイド算出部の処理と上記順列生成部の処理とを交互に実行させる第２演算制御部と、
    を有することを特徴とする信号分離装置。
  6. 請求項４又は５に記載の信号分離装置であって、
    上記クラスタリング部は、
    対応する周波数f∈Fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度のうち、特定の周波数の組合せに対応する事後確率間の類似度のみを指標として用い、上記大域的最適化部で生成された順列Πfを更新し、新たな順列Πfを生成する局所的最適化部を更に有する、
    ことを特徴とする信号分離装置。
  7. 請求項６に記載の信号分離装置であって、
    上記特定の周波数の組合せは、
    周波数差が所定範囲内にある周波数の組合せ、及び／又は、倍音関係にある周波数の組合せである、
    ことを特徴とする信号分離装置。
  8. 周波数領域変換部が、源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時間）を、周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換する周波数領域変換過程と、
    クラス分類部が、周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tを周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f) (n=1,...,N、Ｎ≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出するクラス分類過程と、
    パーミュテーション問題解決部が、対応する周波数fが異なる上記事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成するパーミュテーション問題解決過程と、
    分離部が、上記周波数領域の観測信号Xm(f,τ)と上記パーミュテーション問題解決部で生成された上記事後確率P’(Ck(f)|X(f,τ))とを用い、周波数領域の分離信号Yn(f,τ)を抽出する分離過程と、
    を有することを特徴とする信号分離方法。
  9. 請求項１から７の何れかに記載の信号分離装置としてコンピュータを機能させるためのプログラム。
  10. 請求項９に記載のプログラムを格納したコンピュータ読取り可能な記録媒体。
JP2007218612A 2007-08-24 2007-08-24 信号分離装置、信号分離方法、プログラム及び記録媒体 Active JP4769238B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2007218612A JP4769238B2 (ja) 2007-08-24 2007-08-24 信号分離装置、信号分離方法、プログラム及び記録媒体

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2007218612A JP4769238B2 (ja) 2007-08-24 2007-08-24 信号分離装置、信号分離方法、プログラム及び記録媒体

Publications (2)

Publication Number Publication Date
JP2009053349A JP2009053349A (ja) 2009-03-12
JP4769238B2 true JP4769238B2 (ja) 2011-09-07

Family

ID=40504498

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2007218612A Active JP4769238B2 (ja) 2007-08-24 2007-08-24 信号分離装置、信号分離方法、プログラム及び記録媒体

Country Status (1)

Country Link
JP (1) JP4769238B2 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104429099A (zh) * 2012-04-20 2015-03-18 阿尔卡米斯 用于控制扬声器的工作温度的方法和装置
CN112133321A (zh) * 2020-09-23 2020-12-25 青岛科技大学 一种基于盲源分离的水声信号高斯/非高斯噪声抑制方法

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6059072B2 (ja) * 2013-04-24 2017-01-11 日本電信電話株式会社 モデル推定装置、音源分離装置、モデル推定方法、音源分離方法及びプログラム
EP2830043A3 (en) 2013-07-22 2015-02-18 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. Method for Processing an Audio Signal in accordance with a Room Impulse Response, Signal Processing Unit, Audio Encoder, Audio Decoder, and Binaural Renderer
JP6290803B2 (ja) * 2015-02-24 2018-03-07 日本電信電話株式会社 モデル推定装置、目的音強調装置、モデル推定方法及びモデル推定プログラム
KR101750017B1 (ko) * 2015-11-05 2017-06-23 중앙대학교 산학협력단 무선 이미지 센서 네트워크를 위한 k-평균 클러스터링 기반의 데이터 압축 시스템 및 방법
CN113345456B (zh) * 2021-05-31 2023-06-06 北京小米移动软件有限公司 回声分离方法、装置及存储介质

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104429099A (zh) * 2012-04-20 2015-03-18 阿尔卡米斯 用于控制扬声器的工作温度的方法和装置
CN112133321A (zh) * 2020-09-23 2020-12-25 青岛科技大学 一种基于盲源分离的水声信号高斯/非高斯噪声抑制方法

Also Published As

Publication number Publication date
JP2009053349A (ja) 2009-03-12

Similar Documents

Publication Publication Date Title
US11900947B2 (en) Method and system for automatically diarising a sound recording
EP3479377B1 (en) Speech recognition
JP4769238B2 (ja) 信号分離装置、信号分離方法、プログラム及び記録媒体
JP3949150B2 (ja) 信号分離方法、信号分離装置、信号分離プログラム及び記録媒体
Arberet et al. A robust method to count and locate audio sources in a multichannel underdetermined mixture
US9601119B2 (en) Systems and methods for segmenting and/or classifying an audio signal from transformed audio information
JPWO2006085537A1 (ja) 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体
US9437208B2 (en) General sound decomposition models
JP6195548B2 (ja) 信号解析装置、方法、及びプログラム
Jasa et al. Nested sampling applied in Bayesian room-acoustics decay analysis
JP6538624B2 (ja) 信号処理装置、信号処理方法および信号処理プログラム
JP5994639B2 (ja) 有音区間検出装置、有音区間検出方法、及び有音区間検出プログラム
JP5974901B2 (ja) 有音区間分類装置、有音区間分類方法、及び有音区間分類プログラム
CN110992977B (zh) 一种目标声源的提取方法及装置
JP4630203B2 (ja) 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体、並びに、信号到来方向推定装置、信号到来方向推定方法、信号到来方向推定プログラム及び記録媒体
Jleed et al. Acoustic environment classification using discrete hartley transform features
Shah et al. Inferring room semantics using acoustic monitoring
Wohlmayr et al. Model-based multiple pitch tracking using factorial HMMs: Model adaptation and inference
Cipli et al. Multi-class acoustic event classification of hydrophone data
Kim et al. A subspace approach based on embedded prewhitening for voice activity detection
JP6734237B2 (ja) 目的音源推定装置、目的音源推定方法及び目的音源推定プログラム
JP6114053B2 (ja) 音源分離装置、音源分離方法、およびプログラム
Gburrek et al. On source-microphone distance estimation using convolutional recurrent neural networks
Adiloğlu et al. A general variational Bayesian framework for robust feature extraction in multisource recordings
JP4676920B2 (ja) 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20090729

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20110519

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20110607

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20110617

R150 Certificate of patent or registration of utility model

Ref document number: 4769238

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20140624

Year of fee payment: 3

S531 Written request for registration of change of domicile

Free format text: JAPANESE INTERMEDIATE CODE: R313531

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350