(57)【特許請求の範囲】
【請求項１】
複数の信号源から発せられた信号が混合された混合信号を当該信号に分離する信号分離装置であって、
複数のセンサで観測された前記混合信号を周波数領域の信号である周波数領域信号に変換する周波数領域変換部と、
前記周波数領域信号に独立成分分析を適用し、周波数毎に分離行列を算出する分離行列算出部と、
前記分離行列の一般化逆行列を算出する逆行列算出部と、
前記一般化逆行列を構成する基底ベクトルを正規化し、正規化基底ベクトルを算出する正規化部と、
前記正規化基底ベクトルをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出するクラスタリング部と、
前記セントロイドを用い、周波数応答モデルのモデルパラメータを抽出するモデルパラメータ抽出部と、
前記モデルパラメータを用い、前記分離行列の行を並び替えるための順列を算出する順列決定部と、
を有し、
前記クラスタリング部は、
空間的エイリアシングが起き得ない周波数に対応する前記正規化基底ベクトルのみをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出する、
ことを特徴とする信号分離装置。
【請求項２】
請求項１に記載の信号分離装置であって、
前記クラスタリング部は、
空間的エイリアシングが起き得ない周波数について、暫定順列に従って選択された前記正規化基底ベクトルを足し合わせ、その加算結果のノルムを所定値に正規化したものを暫定セントロイドとするセントロイド算出部と、
空間的エイリアシングが起き得ない周波数毎に、暫定順列に従って選択された前記正規化基底ベクトルと暫定セントロイドとの距離の総和が最小となるように当該暫定順列を更新する暫定順列決定部と、を有し、
前記セントロイド算出部及び前記暫定順列決定部は、
所定の終了条件が満たされるまで、それぞれの処理を交互に実行し、
前記セントロイドは、
所定の終了条件を満たした場合における前記暫定セントロイドである、
ことを特徴とする信号分離装置。
【請求項３】
請求項１に記載の信号分離装置であって、
空間的エイリアシングが起き得る周波数を含む周波数帯域に対応する前記基底ベクトルから生成された情報と、前記順列とを用い、前記モデルパラメータを最適化するモデルパラメータ最適化部をさらに有し、
前記順列決定部は、
最適化された前記モデルパラメータを用い、前記順列を更新する、
ことを特徴とする信号分離装置。
【請求項４】
請求項３に記載の信号分離装置であって、
前記周波数帯域は、全周波数帯域である、
ことを特徴とする信号分離装置。
【請求項５】
請求項３に記載の信号分離装置であって、
前記基底ベクトルから生成された情報は、
前記基底ベクトルの各要素の偏角を、所定のセンサに対応する要素の偏角を基準として正規化し、ノルムを所定値に正規化することにより生成された位相・ノルム正規化ベクトルであり、
前記モデルパラメータ最適化部は、
前記順列に従って関連付けられる前記位相・ノルム正規化ベクトルの要素と前記周波数応答モデルとの距離の総和から構成されるコスト関数に対し、前記モデルパラメータに関する偏微分を行い、その結果に比例する値を前記モデルパラメータから減算することによって当該モデルパラメータを最適化する、
ことを特徴とする信号分離装置。
【請求項６】
請求項３に記載の信号分離装置であって、
前記モデルパラメータ最適化部及び前記順列決定部は、
所定の終了条件が満たされるまで、それぞれの処理を交互に実行する、
ことを特徴とする信号分離装置。
【請求項７】
複数の信号源から発せられた信号が混合された混合信号を当該信号に分離する信号分離装置であって、
複数のセンサで観測された前記混合信号を周波数領域の信号である周波数領域信号に変換する周波数領域変換部と、
前記周波数領域信号に独立成分分析を適用し、周波数毎に分離行列を算出する分離行列算出部と、
前記分離行列の一般化逆行列を算出する逆行列算出部と、
前記一般化逆行列を構成する基底ベクトルを正規化し、位相・ノルム正規化ベクトルを生成する正規化部と、
前記位相・ノルム正規化ベクトルと順列とを用い、モデルパラメータを最適化するモデルパラメータ最適化部と、
最適化された前記モデルパラメータを用い、前記順列を更新する順列決定部と、
を有することを特徴とする信号分離装置。
【請求項８】
複数の信号源から発せられた信号が混合された混合信号を当該信号に分離する信号分離方法であって、
周波数領域変換部が、前記混合信号を周波数領域の信号である周波数領域信号に変換する過程と、
分離行列算出部が、前記周波数領域信号に独立成分分析を適用し、周波数毎に分離行列を算出する過程と、
逆行列算出部が、前記分離行列の一般化逆行列を算出する過程と、
正規化部が、前記一般化逆行列を構成する基底ベクトルを正規化し、正規化基底ベクトルを算出する過程と、
クラスタリング部が、前記正規化基底ベクトルをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出する過程と、
モデルパラメータ抽出部が、前記セントロイドを用い、周波数応答モデルのモデルパラメータを抽出する過程と、
順列決定部が、前記モデルパラメータを用い、前記分離行列の行を並び替えるための順列を算出する過程と、
を有し、
クラスタリング部が、前記正規化基底ベクトルをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出する過程は、
前記クラスタリング部が、空間的エイリアシングが起き得ない周波数に対応する前記正規化基底ベクトルのみをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出する過程である、
ことを特徴とする信号分離方法。
【請求項９】
請求項１から７の何れかに記載の信号分離装置としてコンピュータを機能させるための信号分離プログラム。
【請求項１０】
請求項９に記載の信号分離プログラムを格納したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理の技術分野に属し、特に複数の源信号が空間内で混合された混合信号から源信号を抽出する技術に関する。
【背景技術】
【０００２】
［ブラインド信号分離］
複数の源信号が混在した混合信号から、源信号を推定して分離する信号分離技術としてブラインド信号分離がある。
まず、このブラインド信号分離の定式化を行う。すべての信号はあるサンプリング周波数ｆｓでサンプリングされ、離散的に表現されるものとする。また、Ｎ個の信号が混合されてＭ個のセンサで観測されたとする。以下では、信号の発生源からセンサまでの距離により信号が減衰・遅延し、また壁などにより信号が反射して伝送路の歪みが発生しうる状況を扱う。このような状況での混合は、信号源ｋからセンサｑへのインパルス応答hqk(r)による畳み込み混合
【数１】
0004456571.tif 000002
となる。ここでｔはサンプリング時刻を示している。また、ｓｋ（ｔ）はサンプリング時刻ｔにおいて信号源から発せられる源信号を示しており、ｘｑ（ｔ）はそのセンサｑで観測された混合信号を示している。また、ｒは掃引のための変数である。また、ｋは各信号源に対応付けられた値（k∈{1,...,N}）であり、ｋに対応する信号源を信号源ｋと表記する。さらに、ｑは各センサに対応付けられた値（q∈{1,...,M}）であり、ｑに対応するセンサをセンサｑと表記する。
【０００３】
一般的なインパルス応答hqk(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)を求めることにある。
【０００４】
［周波数領域］
次に、従来のブラインド信号分離の手順について説明する。
ここでは周波数領域において分離の操作を行う。そのためにセンサｑでの混合信号ｘｑ（ｔ）に、フレームの長さ（サンプル点の数）をＬ点とした短時間離散フーリエ変換を適用し、周波数ごとの時間系列
【数２】
0004456571.tif 000003
を求める。ここでｆは周波数であり、f∈{0,fs/L,…,fs(L-1)/L}と離散化されている（ｆｓはサンプリング周波数）。また、τは離散時刻である。また、ｊは虚数単位である。さらにｇ（ｒ）は窓関数であり、ハニング窓
【数３】
0004456571.tif 000004
などのｇ（０）にパワーの中心を持つ窓関数を用いることで、Ｘｑ（ｆ，τ）は離散時刻τを中心とする混合信号ｘｑ（ｔ）の周波数特性を表現する。なお、Ｘｑ（ｆ，τ）はＬサンプルにわたる情報を含んでいるため、すべての離散時刻τに対してＸｑ（ｆ，τ）を求める必要はない。適当な間隔でフレームをシフトしながら、フレームの中心時刻τごとにＸｑ（ｆ，τ）を求める。
【０００５】
周波数領域で処理を行うと、式（１）で示される時間領域での畳み込み混合が、
【数４】
0004456571.tif 000005
と各周波数での単純混合に近似表現でき、分離の操作が単純になる。ここで、Ｈｑｋ（ｆ）は源信号ｋからセンサｑまでの周波数応答であり、Ｓｋ（ｆ，τ）は式（２）と同様な式に従って源信号ｓｋ（ｔ）に短時間離散フーリエ変換を施したものである。式（３）を、ベクトルを用いて表記すると、
【数５】
0004456571.tif 000006
となる。ここで、X(f,τ)=[X1(f,τ),…,XM(f,τ)]Tは周波数領域信号ベクトル、Hk(f)=[H1k(f),…,HMk(f)]Tは信号源ｋから各センサへの周波数応答をまとめたベクトルである。なお、［＊］Ｔは［＊］の転置ベクトルを示す。
【０００６】
［独立成分分析による信号分離］
次に、独立成分分析(ＩＣＡ: Independent Component Analysis)用いて離散周波数ｆ毎に信号分離を行う。ＩＣＡでは周波数領域信号ベクトルX(f,τ)のみから、Ｎ行Ｍ列の分離行列W(f)及び分離信号ベクトル
Y(f,τ)=W(f)X(f,τ)  …(5)
を算出する。ここで、分離行列Ｗ（ｆ）は、分離信号ベクトルY(f,τ)=[Y1(f,τ),...,YN(f,τ)]Tの各要素（分離信号）Y1(f,τ),...,YN(f,τ)が互いに独立になるように算出される。そのためのアルゴリズムには、非特許文献１に記載されているものなどがある。
【０００７】
また、分離行列Ｗ（ｆ）のムーア・ペンローズ（Moore-Penrose）型一般化逆行列Ｗ＋（ｆ）（Ｍ＝Ｎの場合は逆行列Ｗ−１（ｆ）に一致）
[A1(f),...,AN(f)]=W+(f), Ap(f)=[A1p(f),...,AMp(f)]T  …(6)
を計算して式（５）の両辺に掛け合わせると、観測信号ベクトルＸ（ｆ，τ）の分解
【数６】
0004456571.tif 000007
を得る。Ａｐ（ｆ）は基底ベクトルと呼ばれるものであり、各周波数ｆの分離信号Ｙｐ（ｆ，τ）毎に決まる。
【０００８】
［パーミュテーション問題］
ＩＣＡでは信号の独立性に着目して分離を行うため、得られる分離信号Y1(f,τ),...,YN(f,τ)には、順序の任意性がある。それらの順序が入れ替わっても独立性が保たれるからである。この任意性は、同じ源信号Ｓｋ（ｆ，τ）に対応する分離信号Ｙｐ（ｆ，τ）がすべての周波数ｆで同じ添字ｐになるように解決されなければならない。これは、パーミュテーション（permutation）問題と呼ばれ、周波数領域での信号分離において非常に重要な問題である。このパーミュテーション問題を解決する従来手法として、例えば非特許文献２に示されるものなどがある。このような従来手法は、信号源から複数のセンサへの到達時間差や振幅ゲイン比をＩＣＡの結果から抽出し、信号源の方向や位置に相当する値を算出して、それらの値を分類することに基づいている。
【非特許文献１】A. Hyvarinen and J. Karhunen and E. Oja, "Independent Component Analysis," John Wiley & Sons, 2001, ISBN 0-471-40540-X
【非特許文献２】H. Sawada, S. Araki, R. Mukai, and S. Makino, "Blind extraction of a dominant source from mixtures of many sources using ICA and time-frequency masking," IEEE International Symposium on Circuits and Systems (ISCAS 2005), pp. 5882-5885, May 2005.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
しかしながら、これらの従来技術が有効に働くのは、周波数ｆとセンサの間隔ｄの関係が
f＜c/(2d)  …(8)
を満たす（cは信号の速度）場合に限られる。より高い周波数に対しては、いわゆる空間的エイリアシング（aliasing）が起こる可能性があり、信号源の方向や位置に相当する値が正しく算出されない。例えば、図１０（ａ）に示す状況（実験条件Ａ）で周波数領域のブラインド信号分離を適用した場合、式（８）によると、0＜f＜850Hzを満たす範囲の周波数にしか上記の従来技術は有効ではない。残りの850≦f＜8000Hzに対しては、空間的エイリアシングが起こる可能性がある。図１２は、後述する式（１０）で基準センサＱ＝１として基底ベクトルを位相・ノルム正規化した際の正規化基底ベクトル要素Ａ２１’，Ａ２２’の偏角（argument）を縦軸とし、周波数（Frequency）を横軸としたグラフである。この図を用いて空間的エイリアシングの一例を示す。この図から、原点を通る傾きの異なる２本の直線の近くに点が集まっているが、１本目は２ｋＨｚを越えた辺り、２本目は１ｋＨｚ手前辺りで、それぞれ、巡回的ジャンプを起こしていることが分かる。式（８）を満たす周波数帯域ではジャンプを起こさないため、直線の傾きは簡単に算出でき、信号源の方向は直線の傾きから推定できる。しかしながら、ジャンプが起きた周波数より高い周波数では、直線の傾きを算出するのが容易ではなく、間違った傾きを算出してしまうこともある。その結果、間違った方向が推定されてしまう。これが空間的エイリアシングと呼ばれるものである。
【００１０】
このように、パーミュテーション問題に対する上述の従来技術では、空間的エイリアシングが生じる周波数帯域でのパーミュテーション問題を解決することができなかった。
本発明は、このような点に鑑みてなされたものであり、上述した従来の技術を拡張し、空間的エイリアシングが起こり得る周波数帯域でもパーミュテーション問題を解決できる手法を提供することを目的とする。すなわち、本発明の目的は、空間的エイリアシングが起こらない場合にも適用可能であるが、空間的エイリアシングが発生した場合であっても好適にパーミュテーション問題を解決できる手法を提供することである。
【課題を解決するための手段】
【００１１】
本発明では上記課題を解決するために、周波数領域変換部が、複数のセンサで観測された混合信号を周波数領域の信号である周波数領域信号に変換し、分離行列算出部が、周波数領域信号に独立成分分析を適用し、周波数毎に分離行列を算出し、逆行列算出部が、分離行列の一般化逆行列を算出し、正規化部が、一般化逆行列を構成する基底ベクトルを正規化し、正規化基底ベクトルを算出し、クラスタリング部が、正規化基底ベクトルをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出し、モデルパラメータ抽出部が、セントロイドを用い、周波数応答モデルのモデルパラメータを抽出し、順列決定部が、モデルパラメータを用い、分離行列の行を並び替えるための順列を算出する。なお、「一般化逆行列」とは、通常の逆行列をも含む概念であり、分離行列が正方行列である場合、逆行列算出部は、分離行列の逆行列を算出する。また、「周波数応答モデル」とは、周波数応答をモデル化した関数である。
【００１２】
後述するように、独立成分分析の結果が良好に得られた場合、基底ベクトルは、信号源からセンサへの周波数応答に相当する情報を持っている。そして、基底ベクトルを正規化した正規化基底ベクトルの各要素は、対応する周波数応答のモデルを正規化したものに近づく。よって、クラスタリング部が正規化基底ベクトルをクラスタリングすることにより算出したクラスタのセントロイドも、そのクラスタに対応する正規化された周波数応答モデルに近づき、モデルパラメータ抽出部は、セントロイドを用い、周波数応答モデルのモデルパラメータを抽出できる。そして、モデルパラメータが抽出できれば、周波数応答モデルが全周波数について決定できるため、順列決定部は、全周波数について分離行列の行を並び替えるための順列を算出できる。
【００１３】
また、本発明において好ましくは、セントロイド算出部は、空間的エイリアシングが起き得ない周波数に対応する正規化基底ベクトルのみをクラスタリングすることによりクラスタを形成し、当該クラスタのセントロイドを算出する。これにより、空間的エイリアシングの影響を受けることなくセントロイドを算出でき、モデルパラメータの精度を高めることができる。
【００１４】
また、本発明において好ましくは、クラスタリング部は、空間的エイリアシングが起き得ない周波数について、暫定順列に従って選択された正規化基底ベクトルを足し合わせ、その加算結果のノルムを所定値に正規化したものを暫定セントロイドとする処理と、空間的エイリアシングが起き得ない周波数毎に、暫定順列に従って選択された正規化基底ベクトルと暫定セントロイドとの距離の総和が最小となるように当該暫定順列を更新する処理とを、所定の終了条件が満たされるまで交互の繰り返し、所定の終了条件を満たした際の暫定セントロイドをセントロイドとする。これにより、空間的エイリアシングの影響を受けることなく、容易にセントロイドを算出できる。
【００１５】
また、本発明において好ましくは、モデルパラメータ最適化部が、空間的エイリアシングが起き得る周波数を含む周波数帯域に対応する基底ベクトルから生成された情報と、順列とを用い、モデルパラメータを最適化し、順列決定部が、最適化されたモデルパラメータを用い、順列を更新する。これにより、空間的エイリアシングが起き得ない周波数だけではなく、空間的エイリアシングが起き得る周波数にも対応する正規化基底ベクトルの情報を、上述のモデルパラメータの決定に反映させることができ、モデルパラメータの精度をより高めることができる。特に、全周波数帯域に対応する基底ベクトルから生成された情報と、順列とを用い、モデルパラメータを最適化する構成とすれば、モデルパラメータの精度をより一層高めることができる。
【００１６】
また、本発明において好ましくは、上述の基底ベクトルから生成された情報とは、基底ベクトルの各要素の偏角を、所定のセンサに対応する要素の偏角を基準として正規化し、ノルムを所定値に正規化することにより生成された位相・ノルム正規化ベクトルであり、モデルパラメータ最適化部は、順列に従って関連付けられる位相・ノルム正規化ベクトルの要素と周波数応答モデルとの距離の総和から構成されるコスト関数に対し、モデルパラメータに関する偏微分を行い、その結果に比例する値をモデルパラメータから減算することによって当該モデルパラメータを最適化する。なお、「位相・ノルム正規化ベクトルの要素と周波数応答モデルとの距離」とは、位相・ノルム正規化ベクトルの要素の値と周波数応答モデル値との離れ方を何らかの形で示す関数を意味する。これにより、容易にモデルパラメータを最適化することができる。
【００１７】
また、本発明において好ましくは、上述のモデルパラメータ最適化部及び順列決定部は、所定の終了条件が満たされるまで、それぞれの処理を交互に実行する。これにより、モデルパラメータ及び順列の精度を十分高めることができる。或いは、これらの精度を計算量とのトレードオフで適切なものにすることができる。
【００１８】
また、周波数領域変換部が、複数のセンサで観測された前記混合信号を周波数領域の信号である周波数領域信号に変換し、分離行列算出部が、周波数領域信号に独立成分分析を適用し、周波数毎に分離行列を算出し、逆行列算出部が、分離行列の一般化逆行列を算出し、正規化部が、一般化逆行列を構成する基底ベクトルを正規化し、位相・ノルム正規化ベクトルを生成し、モデルパラメータ最適化部が、位相・ノルム正規化ベクトルと順列とを用い、モデルパラメータを最適化し、順列決定部が、最適化されたモデルパラメータを用い、順列を更新する構成であってもよい。この場合、モデルパラメータの初期値には、例えば、適当な値を用いる。そして、モデルパラメータ最適化部と順列決定部との処理を繰り返すことにより、最適なモデルパラメータと順列とが生成できる。
【発明の効果】
【００１９】
以上のように、本発明では、空間的エイリアシングが起こり得る周波数帯域でもパーミュテーション問題を解決できる。
【発明を実施するための最良の形態】
【００２０】
以下、本発明の実施の形態を図面を参照して説明する。以下では、本発明の原理を説明した後、各実施の形態を説明していく。
〔原理〕
＜概要＞
まず、本発明の概要について説明する。
本発明では、まず、空間的エイリアシングが起き得ない周波数帯域におけるＩＣＡの結果を用い、パーミュテーション問題を解決する。その結果を解析すると、信号が混合された様子を表現するモデルのモデルパラメータとして、信号源ｋ（ｋ∈｛１，...，Ｎ｝）からセンサｑ（ｑ∈｛１，...，Ｍ｝）への正規化された到達時間γｑｋと振幅ゲインλｑｋとを算出できる。そして、これらのモデルパラメータγｑｋ，λｑｋを用い、空間的エイリアシングが起き得る周波数帯域におけるパーミュテーション問題を解決する。
【００２１】
ただし、空間的エイリアシングが起き得ない周波数帯域において算出されたモデルパラメータγｑｋ，λｑｋは、限られた周波数帯域でのＩＣＡの結果のみを用いて算出されたものであるため、あまり正確でない可能性がある。そこで、本発明のより好ましい形態では、空間的エイリアシングが起き得ない周波数以外の周波数をも含む周波数帯域（好ましくは、全ての周波数帯域）でのＩＣＡの結果を用い、モデルパラメータγｑｋ，λｑｋに微小な修正値を加算・減算する処理を繰り返し、モデルパラメータγｑｋ，λｑｋを最適化する。
【００２２】
＜詳細＞
次に、本発明の原理の詳細について説明する。
［パーミュテーション問題の解決手法］
本発明の好ましい形態では、分離行列Ｗ（ｆ）の一般化逆行列Ｗ＋（ｆ）の列を１列ずつ取り出した基底ベクトルA1(f),...,AN(f)を全ての周波数ｆについて用い、分離行列Ｗ（ｆ）の行を並び替えるための順列Πｆを、全ての周波数ｆについて生成する。なお、順列とは、｛１，２，…，Ｎ｝から｛１，２，…，Ｎ｝への全単射な関数を意味する。
【００２３】
前述の式（４）と（７）とを比較すれば分かりやすいが、ＩＣＡの結果が良好に得られた場合、基底ベクトルは信号源からセンサへの周波数応答に相当する情報を持つ。従って、基底ベクトルを解析することで、信号が混合された様子を表現するモデルのモデルパラメータを算出することができる。
本発明では、信号源ｋからセンサｑヘの周波数応答Ｈｑｋ（ｆ）をモデル化するものとして、直接波モデル
【数７】
0004456571.tif 000008
を仮定する。なお、ｅはネピア数であり、ｊは虚数単位である。ここで、γｑｋとλｑｋ＞０とは、それぞれ、信号源ｋからセンサｑヘの到達時間と減衰（振幅ゲイン）とを表現する。ただし、源信号ｓｋ（ｆ，τ）の位相及び振幅と、周波数応答Ｈｑｋ（ｆ）の位相及び振幅とを実際に区別することは困難である。そのため、λｑｋ及びγｑｋを相対的なものとして考え、ある種の正規化を行う。到達時間γｑｋは、ある基準センサＱにおけるものを０、すなわち、到達時間γQｋ＝０とする。また、振幅ゲインλｑｋは、
【数８】
0004456571.tif 000009
と正規化する。なお、αは正の実数であり、例えば、α＝１である。
【００２４】
そして、そのような正規化を基底ベクトルAp(f)=[A11(f),...,AMp(f)]Tにも反映する。これは、基底ベクトルＡｐ（ｆ）の位相とノルムとを正規化することにより行う。すなわち、ある基準センサＱに対応する基底ベクトルＡｐ（ｆ）の要素ＡＱｐ（ｆ）の偏角を０として、各要素A11(f),...,AMp(f)の偏角を正規化し、さらにノルムがα（例えば１）になるように（‖・‖＝α）正規化する。このように正規化されたベクトルを、位相・ノルム正規化ベクトルAp'(f)=[A11'(f),...,AMp'(f)]Tと呼ぶこととする。また、このような正規化は、例えば、
【数９】
0004456571.tif 000010
により達成できる。なお、ａｒｇ［・］は、・の偏角を示す。
【００２５】
以上のように正規化を行うと、ある位相・ノルム正規化ベクトルＡｐ’（ｆ）のｑ番目の要素Ａqｐ’（ｆ）は、ある信号源ｋからセンサｑヘの周波数応答をモデル化した式（９）（以下「周波数応答モデル」という）に近づく。ここで、添字ｐとｋとは、算出すベき順列Πｆにより、p=Πf(k)と関連付けられる。従って、コスト関数
【数１０】
0004456571.tif 000011
を最小にする周波数ｆ∈Ｆ毎の順列Πｆと、式（９）のモデルのモデルパラメータλｑｋ，γｑｋとを算出し、順列Πｆに従って分離行列Ｗ（ｆ）の行を並び替えることにより、各周波数でのパーミュテーション問題を解決することができる。なお、Ｆは考慮すべき周波数を全て集めた集合、すなわち、F={0,fs/L,…,fs(L-1)/L}（ｆｓはサンプリング周波数）である。
【００２６】
［モデルパラメータの決定］
式（１１）の問題を解くにあたり、まず、空間的エイリアシングが起こり得ない周波数帯域ＦＬ⊆Ｆの各周波数ｆに対して順列Πｆ（ｆ∈ＦＬ）を決定するとともに、各信号源ｋから各センサｑヘの周波数応答モデル（式（９））のモデルパラメータλｑｋ及びγｑｋの初期値を算出する。なお、空間的エイリアシングが起こり得ない周波数帯域ＦＬは、
FL={f : 0＜f＜c/(2dmax)}  …(12)
によって算出できる。なお、ｄｍａｘは、基準センサＱと他のセンサｑとの距離の最大値である。
【００２７】
空間的エイリアシングが起こり得ない周波数帯域ＦＬだけを考えた場合、式（１１）のコスト関数Ｇは、以下の手順により簡単に最小化できる。
まず、周波数帯域ＦＬに属するすべての周波数ｆに関して、位相・ノルム正規化ベクトルAp'(f)=[A11'(f),...,AMp'(f)]Tの要素Aqp'(f)を
【数１１】
0004456571.tif 000012
に従い、周波数依存成分を排除する周波数正規化を行う。このように正規化された各要素Aqp''(f)からなるベクトルAp''(f)=[A1p''(f),...,AMp''(f)]Tを正規化基底ベクトルと呼ぶ。
【００２８】
このように正規化された正規化基底ベクトルＡｐ''（ｆ）の各要素を用いると、式（１１）に示したコスト関数Ｇは、
【数１２】
0004456571.tif 000013
と変更される。ここでＨｑｋ’は、周波数応答モデルＨｑｋの周波数依存成分を排除する正規化を、式（１３）と同様に行ったものであり、
【数１３】
0004456571.tif 000014
と表現されるものである。なお、Ｈｑｋ’を正規化周波数応答モデルと呼ぶ。
また、式（１４）をベクトル表記すると、
【数１４】
0004456571.tif 000015
となる。なお、Hk'=[ H1k',..., HMk']Tは、正規化周波数応答モデルＨｑｋ’を要素とするベクトルである。
【００２９】
ここで、全てのｐ及びｆを考慮すると、正規化基底ベクトルＡｐ''（ｆ）は、信号源ｋ毎のクラスタを形成する（例えば、非特許文献２参照）。そして、上述の正規化周波数応答モデルＨｑｋ’を要素とするベクトルHk'=[ H1k',..., HMk']Tは、各クラスタのセントロイド（クラスタの中心ベクトル）とみなすことができる。この場合、以下の式（１６）（１７）の演算を交互に繰返し適用することで、式（１５）のコスト関数Ｇ’を最小にするセントロイドＨｋ’と暫定順列Πｆ’を算出することができる。
【００３０】
【数１５】
0004456571.tif 000016
式（１６）は、その演算時点での暫定順列Πｆ’で並び替えられた正規化基底ベクトルの最適なセントロイドＨｋ’（以下「暫定セントロイド」という）を求める演算である。なお、ここで求められる暫定セントロイドＨｋ’は、正規化基底ベクトルと同様、ベクトルとしてのノルムがα（本形態の例では１）になるように正規化されたものとなる。また、式（１７）は、その演算時点での暫定セントロイドＨｋ’から、周波数ｆ∈ＦＬ毎に最適な暫定順列Πｆ’を算出する演算である。なお、式（１７）におけるargminΠ・は、・を最小にする順列Πを意味する。
【００３１】
式（１６）（１７）の演算は、所定の終了条件を満たすまで繰り返される。そして、所定の条件を満たした際の暫定セントロイドＨｋ’をセントロイドＨｋ’とし、以下の式に従ってモデルパラメータλｑｋ及びγｑｋを算出する。
【数１６】
0004456571.tif 000017
【００３２】
[モデルパラメータの最適化と全周波数での順列の決定]
次に、空間的エイリアシングが起こり得る周波数帯域も含め、全ての周波数でのパーミュテーション問題を解決する。単純な方法としては、これまでに得られたモデルパラメータλｑｋ，γｑｋを用い、式（１１）のコスト関数Ｇに従い、全周波数ｆ∈Ｆについて順列Πｆを
【数１７】
0004456571.tif 000018
と決定すればよい。
【００３３】
しかしながら、周波数帯域ＦＬのみで算出したモデルパラメータλｑｋ，γｑｋは、それほど正確でない可能性がある。より良い方法としては、全周波数ｆ∈Ｆの基底ベクトルの情報を用いて、より正確なモデルパラメータλｑｋ，γｑｋを算出することである。これは、式（１１）のコスト関数Ｇのモデルパラメータλｑｋ及びγｑｋそれぞれに関する偏微分値を用い、以下のような修正を繰り返すことで実現できる。
【数１８】
0004456571.tif 000019
なお、μはステップサイズパラメータを示す実数である。また、
【数１９】
0004456571.tif 000020
は、それぞれ、式（１１）のコスト関数Ｇのモデルパラメータλｑｋ，γｑｋに関する偏微分である。また、ｉｍａｇ［・］とｒｅａｌ［・］は、それぞれ、複素数・の虚部と実部を抽出する関数を意味する。
式（１９）による順列Πｆ算出と、式（２０）によるモデルパラメータλｑｋ，γｑｋの修正とを交互に繰返し適用することで、より良い順列Πｆとモデルパラメータλｑｋ，γｑｋを算出することができる。この繰り返しは、所定の終了条件を満たすまで行われる。
【００３４】
［分離行列の行の並び替え］
その後、以上のように得られたΠｆを用い、分離行列Ｗ（ｆ）の行を並び替える。この並び替えは、
Wkq(f)= Wpq(f) ,p=Πf(k)      …(23)
を全てのｑ∈｛１，...，Ｍ｝及びｋ∈｛１，...，Ｎ｝について実行することで行われる。なお、Ｗｐｑ（ｆ）は、分離行列Ｗ（ｆ）のｐ行ｑ列目の要素を示す。
以上により、空間的エイリアシングが起こり得る周波数帯域でもパーミュテーション問題を解決できる。
【００３５】
〔第１の実施の形態〕
次に、本発明における第１の実施の形態について説明する。本形態は、モデルパラメータλｑｋ，γｑｋの最適化を行うことなく、全周波数についての順列Πｆを決定するものである。
＜ハードウェア構成＞
図１は、第１の実施の形態における信号分離装置１のハードウェア構成を例示したブロック図である。
図１に例示するように、この例の信号分離装置１は、ＣＰＵ（Central Processing Unit）１０、入力部２０、出力部３０、補助記憶装置４０、ＲＡＭ（Random Access Memory）１０ｄ、ＲＯＭ（Read Only Memory）５０及びバス７０を有している。
【００３６】
この例のＣＰＵ１０は、制御部１１、演算部１２及びレジスタ１３有し、レジスタ１３に読み込まれた各種プログラムに従って様々な演算処理を実行する。また、この例の入力部２０は、データが入力される入力ポート、キーボード、マウス等であり、出力部３０は、データを出力する出力ポート、ディスプレイ等である。補助記憶装置４０は、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、本形態の信号分離処理を実行するための信号分離プログラムを格納した信号分離プログラム領域４１及びセンサで観測された時間領域の混合信号等の各種データが格納されるデータ領域４２を有している。また、ＲＡＭ５０は、例えば、ＳＲＡＭ （Static Random Access Memory）、ＤＲＡＭ （Dynamic Random Access Memory）等であり、信号分離プログラムが書き込まれる信号分離プログラム領域５１及び各種データが書き込まれるデータ領域５２を有している。また、この例のバス７０は、ＣＰＵ１０、入力部２０、出力部３０、補助記憶装置４０、ＲＡＭ５０及びＲＯＭ６０を通信可能に接続している。
【００３７】
＜ハードウェアとソフトウェアとの協働＞
この例のＣＰＵ１０は、読み込まれたＯＳ（Operating System）プログラムに従い、補助記憶装置４０の信号分離プログラム領域４１に格納されている信号分離プログラムを、ＲＡＭ５０の信号分離プログラム領域５１に書き込む。同様にＣＰＵ１０は、補助記憶装置４０のデータ領域４２に格納されている時間領域の混合信号等の各種データをＲＡＭ５０のデータ領域５２に書き込む。さらに、ＣＰＵ１０は、この信号分離プログラムや各種データが書き込まれたＲＡＭ５０上のアドレスをレジスタ１３に格納する。そして、ＣＰＵ１０の制御部１１は、レジスタ１３に格納されたこれらのアドレスを順次読み出し、読み出したアドレスが示すＲＡＭ５０上の領域からプログラムやデータを読み出し、そのプログラムが示す演算を演算部１２に順次実行させ、その演算結果をレジスタ１３に格納していく。
【００３８】
図２は、このようにＣＰＵ１０に信号分離プログラムが読み込まれることにより構成される信号分離装置１のブロック図の例示である。また、図３は、図２におけるパーミュテーション問題解決部１５０の詳細を例示したブロック図である。
図２に例示するように、本形態の信号分離装置１は、メモリ１００、周波数領域変換部１２０、分離行列算出部１３０、逆行列算出部１４０、パーミュテーション問題解決部１５０、分離信号生成部１６０、時間領域変換部１７０、制御部１８０及び一時メモリ１９０を有している。
【００３９】
ここで、メモリ１００は、記憶領域１０１〜１１５を有している。また、図３に例示するように、パーミュテーション問題解決部１５０は、周波数帯域決定部１５１、正規化部１５２、クラスタリング部１５３、モデルパラメータ抽出部１５４、順列決定部１５５及び正順列分離行列生成部１５６を有している。また、正規化部１５２は、位相・ノルム正規化部１５２ａ及び周波数正規化部１５２ｂを有している。さらに、クラスタリング部１５３は、順列決定部１５３ａ及びセントロイド算出部１５３ｂを有している。
【００４０】
＜処理＞
次に、本形態における信号分離装置１の処理について説明する。なお、以下では、Ｎ個の源信号が混合され、Ｍ個のセンサで観測された状況を取り扱う。また、前処理において、各センサで観測された時間領域の混合信号ｘｑ（ｔ）（q∈{1,...,M}）がメモリ１００の記憶領域１０１に格納され、基準センサＱを示す値Ｑと、基準センサＱと他のセンサｑとの距離の最大値ｄｍａｘとが記憶領域１０５に格納され、暫定順列Πｆ’の初期値が記憶領域１１０に格納されているものとする。
【００４１】
図４は、第１の実施の形態における信号分離装置１の処理を説明するためのフローチャートである。また、図５（ａ）は、図４のステップＳ４の処理の詳細を説明するためのフローチャートである。また、図５（ｂ）及び図６は、図４のステップＳ５の処理の詳細を説明するためのフローチャートである。以下、これらの図を用いて本形態の処理を説明する。なお、以下の各処理は制御部１８０の制御のもと実行される。また、明記しない限り、処理過程の各データは一時メモリ１９０に逐一格納・抽出される。
【００４２】
［周波数領域変換部１２０の処理］
まず、周波数領域変換部１２０（図２）において、メモリ１００の記憶領域１０１から時間領域の混合信号ｘｑ（ｔ）を読み出し、これらを短時間離散フーリエ変換等によって周波数毎の時系列の信号（以下「周波数領域信号」という）Ｘｑ（ｆ，τ）（q∈{1,...,M}）に変換し、メモリ１００の記憶領域１０２に格納する（ステップＳ１）。
【００４３】
［分離行列算出部１３０の処理］
次に、分離行列算出部１３０が、メモリ１００の記憶領域１０２から周波数領域信号Ｘｑ（ｆ，τ）を読み込む。そして、分離行列算出部１３０は、周波数領域信号Ｘｑ（ｆ，τ）に独立成分分析（ＩＣＡ）を適用し、周波数ｆ∈Ｆ毎に分離行列Ｗ（ｆ）を算出する（ステップＳ２）。算出された分離行列Ｗ（ｆ）は、メモリ１００の記憶領域１０３に格納される。
【００４４】
［逆行列算出部１４０の処理］
次に、逆行列算出部１４０が、メモリ１００の記憶領域１０３から分離行列Ｗ（ｆ）を読み込む。そして、逆行列算出部１４０は、分離行列Ｗ（ｆ）の一般化逆行列Ｗ＋（ｆ）を算出する（ステップＳ３）。なお、一般化逆行列としては、例えば、ムーア・ペンローズ（Moore-Penrose）型擬似逆行列等の誤差最小型の一般化逆行列を用いることが望ましい。また、分離行列Ｗ（ｆ）が正方行列の場合、一般化逆行列Ｗ＋（ｆ）は、分離行列Ｗ（ｆ）の逆行列と一致する。算出された一般化逆行列Ｗ＋（ｆ）は、メモリ１００の記憶領域１０４に格納される。
【００４５】
［パーミュテーション問題解決部１５０の処理］
次に、パーミュテーション問題解決部１５０の正規化部１５２（図３）が、メモリ１００の記憶領域１０４に格納された一般化逆行列Ｗ＋（ｆ）（=[A1(f),..., AN(f)]）を構成する基底ベクトルＡｐ（ｆ）を正規化し、正規化基底ベクトルＡｐ''（ｆ）を算出する（ステップＳ４）。算出された正規化基底ベクトルＡｐ''（ｆ）は、メモリ１００の記憶領域１０８に格納される。以下に、このステップＳ４の詳細を示す。
【００４６】
［ステップＳ４の処理の詳細］
まず、周波数帯域決定部１５１が、メモリ１００の記憶領域１０５から値ｄｍａｘを読み込み、空間的エイリアシングが起き得ない周波数帯域ＦＬを
FL={f: 0＜f＜c/(2dmax)}
によって算出する（ステップＳ２１）。算出された値ＦＬは、メモリ１００の記憶領域１０６に格納される。
【００４７】
次に、正規化部１５２の位相・ノルム正規化部１５２ａが、メモリ１００の記憶領域１０４から基底ベクトルＡｐ（ｆ）を読み込む。そして、位相・ノルム正規化部１５２ａは、基底ベクトルＡｐ（ｆ）の各要素A1p(f),...,AMp(f)の偏角を、所定の基準センサＱに対応する要素AQp(f)の偏角を基準として正規化し、ノルムを所定値（本形態の例では１）に正規化した位相・ノルム正規化ベクトルAp’(f)=[A1p’(f),...,AMp’(f)]T(f∈F)を生成する（ステップＳ２２）。この正規化は、例えば、前述の式（１０）に従って行われるが、特に式（１０）に限定されるものではない。生成された位相・ノルム正規化ベクトルＡｐ’（ｆ）は、メモリ１００の記憶領域１０７に格納される。
【００４８】
次に、正規化部１５２の周波数正規化部１５２ｂが、メモリ１００の記憶領域１０６から値ＦＬを、記憶領域１０７から位相・ノルム正規化ベクトルＡｐ’（ｆ）を、それぞれ読み込む。そして、周波数正規化部１５２ｂは、周波数帯域ＦＬに属する全ての周波数ｆについて、位相・ノルム正規化ベクトルAp’(f)=[A1p’(f),...,AMp’(f)]Tの各要素Ａｑｐ’（ｆ）の偏角を周波数ｆに比例した値で除算する正規化を行い、周波数依存性を排除した正規化基底ベクトルAp''(f)=[A1p''(f),...,AMp''(f)]T(f∈FL)を生成する(ステップＳ２３)。なお、この正規化は、例えば、前述の式（１３）に従って行われるが、特に式（１３）に限定されるものではない。生成された正規化基底ベクトルＡｐ''（ｆ）は、メモリ１００の記憶領域１０８に格納される（［ステップＳ４の処理の詳細］の説明終わり）。
【００４９】
次に、クラスタリング部１５３が、メモリ１００の記憶領域１０８に格納された正規化基底ベクトルＡｐ''（ｆ）をクラスタリングすることによりクラスタＣｋを形成し、当該クラスタＣｋのセントロイドＨｋ’を算出する（ステップＳ５）。算出されたセントロイドＨｋ’は、メモリ１００の記憶領域１０９に格納される。以下に、このステップＳ５の詳細を示す。
【００５０】
［ステップＳ５の処理の詳細（パターンＡ）］
まず、クラスタリング部１５３の順列決定部１５３ａが、メモリ１００の記憶領域１１０から暫定順列Πｆ’の初期値を読み込む（ステップＳ３１）。次に、順列決定部１５３ａは、メモリ１００の記憶領域１０６から値ＦＬを読み込み、記憶領域１０８から正規化基底ベクトルＡｐ''（ｆ）を読み込む。そして、順列決定部１５３ａは、周波数ｆ∈ＦＬについて、暫定順列Πｆ'に従って関連付けられた正規化基底ベクトルＡΠｆ'(k)''（ｆ）を足し合わせ、その加算結果のノルムを所定値（本形態の例では１）に正規化したものを暫定セントロイドＨｋ’とする（ステップＳ３２）。すなわち、順列決定部１５３ａは、例えば、前述の式（１６）に従い、暫定セントロイドＨｋ’を算出する。なお、下付添え字のΠｆ'（ｋ）は、Πｆ'（ｋ）を示す。また、算出された暫定セントロイドＨｋ’は、メモリ１００の記憶領域１０９に格納される。
【００５１】
次に、セントロイド算出部１５３ｂが、メモリ１００の記憶領域１０６から値ＦＬを読み込み、記憶領域１０８から正規化基底ベクトルＡｐ''（ｆ）を読み込み、記憶領域１０９から暫定セントロイドＨｋ’を読み込む。そして、セントロイド算出部１５３ｂは、周波数ｆ∈ＦＬ毎に、暫定順列Πｆ'に従って関連付けられた正規化基底ベクトルＡΠｆ'(k)''（ｆ）と暫定セントロイドＨｋ’との距離の総和（周波数ｆ∈ＦＬ毎の総和）が最小となるように暫定順列Πｆ'を更新する（ステップＳ３３）。すなわち、セントロイド算出部１５３ｂは、例えば、前述の式（１７）に従い、暫定順列Πｆ'を更新する。更新された暫定順列Πｆ'は、メモリ１００の記憶領域１１０に上書きされる。なお、「記憶領域１１０に上書き」とは、記憶領域１１０に格納されていたデータを削除し、その代わりに記憶領域１１０に新たなデータを格納すること、或いは、記憶領域１１０に格納されていたデータを他の記憶領域に移し、その代わりに記憶領域１１０に新たなデータを格納することを意味する。
【００５２】
次に、制御部１８０が、所定の終了条件を満たした否かを判断する（ステップＳ３４）。この終了条件としては、例えば、「ステップＳ３２，Ｓ３３の処理を交互に所定回数繰り返したこと」や、「ステップＳ３３での更新された暫定順列Πｆ'と更新前の暫定順列Πｆ'とが、全ての周波数において同一であること」等を例示できる。ここで、制御部１８０が、所定の終了条件を満たしていないと判断した場合、制御部１８０は、処理をステップＳ３２に戻す。一方、制御部１８０が、所定の終了条件を満たしたと判断した場合には、制御部１８０は、ステップＳ５の処理を終了させる。なお、ステップＳ５の処理が終了した時点でメモリ１００の記憶領域１０９に格納されている暫定セントロイドＨｋ’を、正規のセントロイドＨｋ’として取り扱う（［ステップＳ５の処理の詳細（パターンＡ）］の説明終わり）。
【００５３】
［ステップＳ５の処理の詳細（パターンＢ）］
また、上述のパターンＡの代わりに、以下のようにステップＳ５の処理を実行してもよい。
まず、クラスタリング部１５３の順列決定部１５３ａが、メモリ１００から暫定セントロイドＨｋ’の初期値を読み込む。この暫定セントロイドＨｋ’としては、例えば、ある周波数の正規化基底ベクトルＡｐ''（ｆ）を用いる（ステップＳ３５）。
次に、セントロイド算出部１５３ｂが、メモリ１００の記憶領域１０６から値ＦＬを読み込み、記憶領域１０８から正規化基底ベクトルＡｐ''（ｆ）を読み込み、記憶領域１０９から暫定セントロイドＨｋ’を読み込む。そして、セントロイド算出部１５３ｂは、周波数ｆ∈ＦＬ毎に、暫定順列Πｆ'に従って関連付けられた正規化基底ベクトルＡΠｆ'（ｋ）''（ｆ）と暫定セントロイドＨｋ’との距離の総和（周波数ｆ∈ＦＬ毎の総和）が最小となるように暫定順列Πｆ'を生成する（ステップＳ３６）。すなわち、セントロイド算出部１５３ｂは、例えば、前述の式（１７）に従い、暫定順列Πｆ'を生成する。生成された暫定順列Πｆ'は、メモリ１００の記憶領域１１０に上書きされる。
【００５４】
次に、順列決定部１５３ａは、メモリ１００の記憶領域１０６から値ＦＬを読み込み、記憶領域１０８から正規化基底ベクトルＡｐ''（ｆ）を読み込む。そして、順列決定部１５３ａは、周波数ｆ∈ＦＬについて、暫定順列Πｆ'に従って関連付けられた正規化基底ベクトルＡΠｆ' (k)''（ｆ）を足し合わせ、その加算結果のノルムを所定値（本形態の例では１）に正規化したものを暫定セントロイドＨｋ’とする（ステップＳ３７）。すなわち、順列決定部１５３ａは、例えば、前述の式（１６）に従い、暫定セントロイドＨｋ’を算出する。算出された暫定セントロイドＨｋ’は、メモリ１００の記憶領域１０９に格納される。
【００５５】
次に、制御部１８０が、上述のパターンＡと同様に、所定の終了条件を満たした否かを判断する（ステップＳ３８）。ここで、制御部１８０が、所定の終了条件を満たしていないと判断した場合、制御部１８０は、処理をステップＳ３６に戻す。一方、制御部１８０が、所定の終了条件を満たしたと判断した場合には、制御部１８０は、ステップＳ５の処理を終了させる。なお、ステップＳ５の処理が終了した時点でメモリ１００の記憶領域１０９に格納されている暫定セントロイドＨｋ’を、正規のセントロイドＨｋ’として取り扱う（［ステップＳ５の処理の詳細（パターンＢ）］の説明終わり）。
【００５６】
次に、モデルパラメータ抽出部１５４が、メモリ１００の記憶領域１０９からセントロイドＨｋ’を読み込む。そして、モデルパラメータ抽出部１５４は、セントロイドＨｋ’を用い、周波数応答モデルＨｋのモデルパラメータγｑｋ，λｑｋを抽出する（ステップＳ６）。具体的には、例えば、前述の式（１３）に従って周波数正規化部１５２ｂでの周波数正規化を行っていたのであれば、モデルパラメータ抽出部１５４は、前述した式（１８）に従って周波数応答モデルＨｋのモデルパラメータγｑｋ，λｑｋを算出する。算出されたモデルパラメータγｑｋ，λｑｋは、メモリ１００の記憶領域１１１に格納する。
【００５７】
次に、順列決定部１５５が、全周波数ｆ∈Ｆ毎に、メモリ１００の記憶領域１１１に格納されたモデルパラメータγｑｋ，λｑｋを用い、分離行列Ｗ（ｆ）の行を並び替えるための順列Πｆを算出する（ステップＳ７）。より具体的には、順列決定部１５５は、まず、メモリ１００の記憶領域１０７から位相・ノルム正規化ベクトルＡｐ’（ｆ）を読み込み、記憶領域１１１からモデルパラメータγｑｋ，λｑｋを読み込む。そして、順列決定部１５５は、前述の式（１９）に従って順列Πｆを算出する。算出された順列Πｆは、メモリ１００の記憶領域１１２に格納される。
【００５８】
次に、正順列分離行列生成部１５６が、メモリ１００の記憶領域１０３から分離行列Ｗ（ｆ）を読み込み、記憶領域１１２から順列Πｆを読み込む。そして、正順列分離行列生成部１５６は、順列Πｆを用い、前述の式（２３）を、全てのｑ∈｛１，...，Ｍ｝及びｋ∈｛１，...，Ｎ｝について実行することによって分離行列Ｗ（ｆ）の行を並び替え、正順列分離行列Ｗ’（ｆ）を生成する（ステップＳ８）。生成された正順列分離行列Ｗ’（ｆ）は、メモリ１００の記憶領域１１３に格納される。
【００５９】
［分離信号生成部１６０の処理］
次に、分離信号生成部１６０（図２）が、メモリ１００の記憶領域１０２から周波数領域信号Ｘｑ（ｆ，τ）を読み込み、記憶領域１１３から正順列分離行列Ｗ’（ｆ）を読み込む。そして、分離信号生成部１６０は、正順列分離行列Ｗ’（ｆ）を、周波数領域信号Ｘｑ（ｆ，τ）からなる周波数領域信号ベクトルX(f,τ)=［X1(f,τ),...,XM(f,τ) ］Tに適用し（Y(f,τ)= W’(f)・X(f,τ)）、分離信号ベクトルY(f,τ)=［Y1(f,τ),...,YM(f,τ) ］Tを生成する（ステップＳ９）。生成された分離信号ベクトルＹ（ｆ，τ）は、メモリ１００の記憶領域１１４に格納される。
【００６０】
［時間領域変換部１７０の処理］
次に、時間領域変換部１７０が、メモリ１００の記憶領域１１４から分離信号ベクトルY(f,τ)=[Y1(f,τ),...,YN(f,τ)]Tを読み込む。そして、時間領域変換部１７０は、短時間逆フーリエ変換等によって、分離信号ベクトルY(f,τ)=[Y1(f,τ),...,YN(f,τ)]Tを時間領域の分離信号ベクトルy(t)=[y1(t),...,yN(t)]Tに変換する（ステップＳ１０）。生成された分離信号ベクトルｙ（ｔ）はメモリ１００の記憶領域１１２に格納される。
【００６１】
〔第２の実施の形態〕
次に、本発明における第２の実施の形態について説明する。本形態は、モデルパラメータλｑｋ，γｑｋの最適化を行って、全周波数についての順列Πｆを決定するものである。なお、以下では、第１の実施の形態との相違点を中心に説明し、第１の実施の形態と共通する事項については説明を省略する。
【００６２】
＜構成＞
第１の実施の形態との相違点はパーミュテーション問題解決部のみである。以下では、パーミュテーション問題解決部の構成のみを示す。
図７は、第２の実施の形態におけるパーミュテーション問題解決部２５０の詳細を例示したブロック図である。なお、図７において、第１の実施の形態と共通する部分については、図３と同じ符号を付し、説明を省略する。
図７に例示するように、本形態のパーミュテーション問題解決部２５０は、周波数帯域決定部１５１、正規化部１５２、クラスタリング部１５３、モデルパラメータ抽出部１５４、最適化部２５５及び正順列分離行列生成部１５６を有している。ここで、最適化部２５５は、順列決定部２５５ａ及びモデルパラメータ最適化部２５５ｂを有している。
【００６３】
＜処理＞
次に、本形態における信号分離装置の処理について説明する。なお、以下でも、Ｎ個の源信号が混合され、Ｍ個のセンサで観測された状況を取り扱う。また、第１の実施の形態と同じ前処理が行なわれているものとする。
図８は、第２の実施の形態における信号分離装置の処理を説明するためのフローチャートである。また、図９は、図８に示したステップＳ１０８の処理の詳細を説明するためのフローチャートである。以下、これらの図を用いて本形態の処理を説明する。なお、以下の各処理は制御部１８０の制御のもと実行される。また、明記しない限り、処理過程の各データは一時メモリ１９０に逐一格納・抽出される。
【００６４】
［周波数領域変換部１２０・分離行列算出部１３０・逆行列算出部の処理］
本形態の周波数領域変換部１２０・分離行列算出部１３０・逆行列算出部１４０の処理（ステップＳ１０１〜Ｓ１０３）の処理は、第１の実施の形態のステップＳ１〜Ｓ３と同じである。
【００６５】
［パーミュテーション問題解決部２５０の処理］
周波数帯域決定部１５１、正規化部１５２、クラスタリング部１５３及びモデルパラメータ抽出部１５４の処理（ステップＳ１０５〜Ｓ１０７）は、第１の実施の形態（ステップＳ５〜Ｓ７）と同じである。第１の実施の形態との相違点は、最適化部２５５（図７）が、空間的エイリアシングが起き得る周波数を含む周波数帯域に対応する位相・ノルム正規化ベクトルＡｐ’（ｆ）と、順列Πｆとを用い、モデルパラメータγｑｋ，λｑｋを最適化し、最適化されたモデルパラメータγｑｋ，λｑｋを用い、順列Πｆを更新する処理を行う点である（ステップＳ１０８）。以下に、この処理の詳細を示す。
【００６６】
［ステップＳ１０８の処理の詳細］
まず、最適化部２５５の順列決定部２５５ａが、メモリ１００の記憶領域１０７から位相・ノルム正規化ベクトルAp’(f)を読み込み（f∈F）、記憶領域１１１からモデルパラメータγｑｋ，λｑｋを読み込む。そして、順列決定部２５５ａは、位相・ノルム正規化ベクトルAΠf(k)’(f)の要素AqΠf(k)’(f)と周波数応答モデルλqk・exp(-j2πfγqk)との距離の周波数ｆ毎の総和を、それぞれ最小にする順列Πｆを算出する（ステップＳ１２１）。すなわち、順列決定部２５５ａは、例えば、前述した式（１９）に従って順列Πｆを算出する。なお、なお、AΠｆ（ｋ）’(f)やAｑΠｆ（ｋ）’(f)の下付き添え字のΠｆ（ｋ）は、Πｆ（ｋ）を意味する。また、expはネピア数である。算出された順列Πｆは、メモリ１００の記憶領域１１２に格納される。
【００６７】
次に、最適化部２５５のモデルパラメータ最適化部２５５ｂが、メモリ１００の記憶領域１０７から位相・ノルム正規化ベクトルAp’(f)を読み込み（f∈F）、記憶領域１１１からモデルパラメータγｑｋ，λｑｋを読み込み、記憶領域１１２から順列Πｆを読み込む。そして、モデルパラメータ最適化部２５５ｂは、順列Πｆに従って関連付けられる位相・ノルム正規化ベクトルＡΠｆ（ｋ）’（ｆ）の要素ＡｑΠｆ（ｋ）’（ｆ）と周波数応答モデルλqk・exp(-j2πfγqk)との距離の総和から構成されるコスト関数Ｇに対し、モデルパラメータγｑｋ，λｑｋに関する偏微分を行い、それらの結果に比例する値をモデルパラメータγｑｋ，λｑｋからそれぞれ減算することによって、各モデルパラメータγｑｋ，λｑｋを最適化する（ステップＳ１２２）。すなわち、モデルパラメータ最適化部２５５ｂは、前述の式（２０）に従って、各モデルパラメータγｑｋ，λｑｋを最適化する。更新された各モデルパラメータγｑｋ，λｑｋは、メモリ１００の記憶領域１１１に格納される。
【００６８】
次に、順列決定部２５５ａが、メモリ１００の記憶領域１０７から位相・ノルム正規化ベクトルＡｐ’（ｆ）を読み込み（f∈F）、記憶領域１１１からモデルパラメータγｑｋ，λｑｋを読み込む。そして、順列決定部２５５ａは、位相・ノルム正規化ベクトルＡΠｆ（ｋ）’（ｆ）の要素ＡｑΠｆ（ｋ）’（ｆ）と周波数応答モデルλqk・exp(-j2πfγqk)との距離の周波数ｆ毎の総和を、それぞれ最小にする順列Πｆを算出し、これをメモリ１００の記憶領域１１２に上書きする（ステップＳ１２３）。なお、ここでの順列Πｆの算出も、例えば、前述した式（１９）に従って行われる。
【００６９】
次に、制御部１８０が、所定の終了条件を満たした否かを判断する（ステップＳ１２４）。この終了条件としては、例えば、「ステップＳ１２２，Ｓ１２３の処理を交互に所定回数繰り返したこと」や、「ステップＳ１３３での更新された順列Πｆと更新前の順列Πｆとが、全ての周波数において同一であること」等を例示できる。ここで、制御部１８０が、所定の終了条件を満たしていないと判断した場合、制御部１８０は、処理をステップＳ１２２に戻す。一方、制御部１８０が、所定の終了条件を満たしたと判断した場合には、制御部１８０は、ステップＳ１０８の処理を終了させる（［ステップＳ１０８の処理の詳細］の説明終わり）。
【００７０】
その後、正順列分離行列生成部１５６、分離信号生成部１６０及び時間領域変換部１７０によって、第１の実施の形態のステップＳ８〜Ｓ１０と同じ処理が行われる（ステップＳ１０９〜Ｓ１１１）。
〔実験〕
次に、上述した各形態の効果を示すために、図１０（ａ）に示した実験条件Ａと、さらに図１０（ｂ）に示した実験条件Ｂについて複数の音声を分離する実験を行った結果を示す。なお、比較のため、パーミュテーション問題の解決に関しては、以下の３種類の手順で実験を行った。
【００７１】
手順Ｉ：パーミュテーション問題は、空間的エイリアシングが起きない低域の周波数帯域ＦＬでのみ解決し、その他の高域ではＩＣＡで得られた結果をそのまま用いた。
手順II：全ての周波数ｆ∈Ｆでパーミュテーション問題を解決したが、低域の周波数帯域ＦＬのみを用いた前述の式（１８）に従って算出したモデルパラメータλｑｋ，γｑｋを用いて解決した。
手順II：式（１９）と（２０）の繰返しにより最適化したモデルパラメータλｑｋ，γｑｋを用い、全ての周波数ｆ∈Ｆでパーミュテーション問題を解決した。
【００７２】
図１１（ａ）は、上述の各条件に対応する分離性能を、SIR（Signal-to-Interference Ratio，分離信号における目的信号成分と干渉信号成分のパワー比）で示した表である。この表から、双方の実験条件において、手順II、手順IIIと分離性能が向上していくことが見受けられ、各形態の効果が確認できた。なお、入力SIRとは、センサでの観測信号のSIRである。
次に、式（１９）と（２０）を繰返し適用することで、モデルパラメータλｑｋ，γｑｋがどのように改善されるかの一例を示す。
【００７３】
図１１（ｂ）（ｃ）は、順列Πｆに従って分離行列Ｗ（ｆ）の行を並び替えた正順列分離行列Ｗ’（ｆ）の一般化逆行列Ｗ’+（ｆ）を構成する基底ベクトルＡ２１（ｆ），Ａ２２（ｆ）の偏角をプロットしたグラフである。ここで、図１１（ｂ）は、空間的エイリアシングが起きない低域の周波数帯域ＦＬでモデルパラメータλｑｋ，γｑｋを算出し、式（１９）により順列Πｆを決定した場合の基底ベクトルＡ２１（ｆ），Ａ２２（ｆ）の偏角をプロットしたグラフである。また、図１１（ｃ）は、式（１９）と式（２０）の繰返しにより、全帯域の基底ベクトルの情報を用いてモデルパラメータを改善し、順列Πｆを決定した場合の基底ベクトルＡ２１（ｆ），Ａ２２（ｆ）の偏角をプロットしたグラフである。なお、これらのグラフにおいて、横軸は周波数（Frequency〔Hz〕）を示し、縦軸は基底ベクトルＡ２１（ｆ），Ａ２２（ｆ）の偏角（Argument）を示す。また、白丸マークは基底ベクトルＡ２１（ｆ）に対応するデータを示し、十字マークは基底ベクトルＡ２２（ｆ）に対応するデータを示す。また、一点鎖線はexp(-j2πfγ21)を示し、破線はexp(-j2πfγ22)を示す。また、これらのデータは、図１０（ａ）の実験条件Ａを前提としたものであり、式（１０）において、基準センサＱ＝１として基底ベクトルの位相とノルムとを正規化した。
【００７４】
図１１（ｂ）の場合、遅延時間を示すモデルパラメータλｑｋ，γｑｋがあまり正確でないため、高城で順列決定のエラーが起きている（例えぱ４ｋＨｚを少し越えた辺り）。これに対し、図１１（ｃ）の場合には、モデルパラメータλｑｋ，γｑｋが正確に得られておち、高城でも順列が正しく決定されていることが分かる。
【００７５】
〔各実施の形態の効果〕
以上、上述した各形態の効果をまとめる。各形態の手法により、空間的エイリアシングが起こる高城でも、パーミュテーション問題が正しく解決できる。さらに、第２の実施の形態のように、全ての帯域でのＩＣＡの結果を用いてモデルパラメータλｑｋ，γｑｋを改善することにより、その正確さが高まる。
【００７６】
〔変形例等〕
なお、変形例として、例えば、図１３に示すように、正規化部１５２の位相・ノルム正規化部１５２ａと、最適化部２５５の順列決定部２５５ａ及びモデルパラメータ最適化部２５５ｂと、正順列分離行列生成部１５６とのみによって、パーミュテーション問題解決部３５０を構成してもよい。この場合、正規化部１５２の位相・ノルム正規化部１５２ａは、上述のように一般化逆行列Ａ（ｆ）を構成する基底ベクトルＡｐ（ｐ）を正規化して位相・ノルム正規化ベクトルＡｐ''（ｐ）を生成し、メモリ１００に格納する。そして、適当な値をモデルパラメータλｑｋ，γｑｋの初期値とし、モデルパラメータ最適化部２５５ｂが、メモリ１００から読み込んだ位相・ノルム正規化ベクトルＡｐ''（ｐ）と順列Πｆとを用い、第２の実施の形態と同様にモデルパラメータλｑｋ，γｑｋを最適化してメモリ１００に格納する。また、順列決定部２５５ａは、メモリ１００から読み込んだモデルパラメータλｑｋ，γｑｋを用い、第２の実施の形態と同様に、メモリ１００の順列Πｆを更新する。そして、モデルパラメータ最適化部２５５ｂ及び順列決定部２５５ａの処理を繰り返し、最適な順列Πｆを生成する。そして、正順列分離行列生成部１５６は、生成された順列Πｆを用い、分離行列Ｗ（ｆ）の行を並び替え、正順列分離行列Ｗ’（ｆ）を生成する。この方法でも、空間的エイリアシングが起こり得る周波数帯域でのパーミュテーション問題の解決が可能である。
【００７７】
また、例えば、上記の各実施の形態では、空間的エイリアシングが起こりえない周波数に対応する基底ベクトルのみをクラスタリングし、セントロイドを算出し、それからモデルパラメータを求めた。しかし、空間的エイリアシングが起こりえない周波数に対応する基底ベクトルに、多少、空間的エイリアシングが起こりえる周波数に対応する基底ベクトルが混在したものをクラスタリングし、セントロイドを算出し、それからモデルパラメータを求めてもよい。
また、第２の実施の形態では、全周波数（f∈F）に対応する位相・ノルム正規化ベクトルを用いてモデルパラメータを最適化したが、全周波数に対応する位相・ノルム正規化ベクトルを用いずにモデルパラメータを最適化してもよい。
【００７８】
また、上述の第２の実施の形態では、最適化部２５５の順列決定部２５５ａとモデルパラメータ最適化部２５５ｂの処理を繰り返す構成とした（ステップＳ１２２〜Ｓ１２４）。しかし、これらの処理を繰り返すことなく、１回のステップＳ１２２〜Ｓ１２４の処理のみでモデルパラメータを最適化する構成であってもよい。
また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００７９】
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（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）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００８１】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【００８２】
本発明の産業上の利用分野としては、例えば、数十センチ間隔でマイクロホンを配置したテレビ会議システムや、高いサンプリング周波数でサンプリングした音声信号の信号分離を行うシステム等のように、空間的エイリアシングがおきやすい環境で利用されるシステムを例示できる。本発明により、このように空間的エイリアシングがおきやすい環境でも、全周波数について正しくパーミュテーション問題が解決できる。
【図面の簡単な説明】
【００８３】
【図１】図１は、第１の実施の形態における信号分離装置のハードウェア構成を例示したブロック図である。
【図２】図２は、第１の実施の形態における信号分離装置のブロック図の例示である。
【図３】また、図３は、図２におけるパーミュテーション問題解決部の詳細を例示したブロック図である。
【図４】図４は、第１の実施の形態における信号分離装置の処理を説明するためのフローチャートである。
【図５】図５（ａ）は、図４のステップＳ４の処理の詳細を説明するためのフローチャートであり、図５（ｂ）は、図４のステップＳ５の処理の詳細を説明するためのフローチャートである。
【図６】図６は、図４のステップＳ５の処理の詳細の詳細を説明するためのフローチャートである。
【図７】図７は、第２の実施の形態におけるパーミュテーション問題解決部の詳細を例示したブロック図である。
【図８】図８は、第２の実施の形態における信号分離装置の処理を説明するためのフローチャートである。
【図９】図９は、図８のステップＳ１０８の処理の詳細を説明するためのフローチャートである。
【図１０】図１０（ａ）は実験条件Ａを示した表であり、図４（ｂ）は実験条件Ｂを示した図である。
【図１１】図１１（ａ）は、各条件に対応する分離性能をSIRで示した表である。また、図１１（ｂ）は、空間的エイリアシングが起きない低域の周波数帯域ＦＬでモデルパラメータλｑｋ，γｑｋを算出し、式（１９）により順列Πｆを決定した場合の基底ベクトルＡ２１（ｆ），Ａ２２（ｆ）の偏角をプロットしたグラフである。また、図１１（ｃ）は、式（１９）と式（２０）の繰返しにより、全帯域の基底ベクトルの情報を用いてモデルパラメータを改善し、順列Πｆを決定した場合の基底ベクトルＡ２１（ｆ），Ａ２２（ｆ）の偏角をプロットしたグラフである。
【図１２】図１２は、基準センサＱ＝１として基底ベクトルを位相・ノルム正規化した際の正規化基底ベクトル要素Ａ２１’とＡ２２’の偏角（argument）を縦軸とし、周波数（Frequency）を横軸としたグラフである。
【図１３】図１３は、パーミュテーション問題解決部の変形例を示したブロック図である。
【符号の説明】
【００８４】
１  信号分離装置
【図１】
0004456571.tif 000021
【図２】
0004456571.tif 000022
【図３】
0004456571.tif 000023
【図４】
0004456571.tif 000024
【図５】
0004456571.tif 000025
【図６】
0004456571.tif 000026
【図７】
0004456571.tif 000027
【図８】
0004456571.tif 000028
【図９】
0004456571.tif 000029
【図１０】
0004456571.tif 000030
【図１１】
0004456571.tif 000031
【図１２】
0004456571.tif 000032
【図１３】
0004456571.tif 000033