(57)【要約】【課題】目的信号のインパルス応答ベクトル等の事前知識がなくても、ビームフォーマを用いて、混合信号から目的信号を抽出する。【解決手段】Ｎ個の信号源から発せられた信号をＭ個のセンサ４ｍで観測し、１個以上の信号を抽出する場合において、ただし、Ｎ、Ｍは２以上の整数であり、センサ４ｍで観測された観測信号ｘｍ（ｔ）を周波数領域の信号ｘｎ（ｆ、τ）に変換し（５）、ｘｎ（ｆ、τ）に対し、正規化を行い（２２）、正規化されたｘｎ（ｆ、τ）をＮ個のクラスタＣｎにクラスタリングし（２４）、Ｃｎから、不要信号のみが含まれる観測信号の相関行列ＲｎＪ（ｆ）を推定し（２５）、クラスタの情報ＣｎとＲｎＪ（ｆ）とからビームフォーマｗｎ（ｆ）を計算し（２８）、ｗｎ（ｆ）を用い、ｘｎ（ｆ、τ）から目的信号ｙｎ（ｆ、τ）を抽出し（３０）、ｙｎ（ｆ、τ）を時間領域の信号に変換する（３２）。【選択図】図３
2008060635.tif 000002
【特許請求の範囲】
【請求項１】
Ｎ個の信号源から発せられた信号をＭ個のセンサで観測し、観測された信号のうち、１個以上の信号を抽出する信号抽出装置において、ただし、Ｎ、Ｍは２以上の整数であり、
上記Ｍ個のセンサで観測された観測信号を周波数領域の信号に変換する周波数領域変換部と、
上記周波数領域の信号に対し、正規化を行い、正規化観測信号ベクトルを算出する正規化部と、
上記正規化観測信号ベクトルをＮ個のクラスタにクラスタリングするクラスタリング部と、
上記クラスタの情報から、不要信号のみが含まれる観測信号の相関行列である不要信号相関行列を推定する不要信号相関行列推定部と、
上記クラスタの情報と、上記不要信号相関行列と、からビームフォーマを計算するビームフォーマ計算部と、
上記ビームフォーマを用い、上記周波数領域の信号から目的信号を抽出する目的信号抽出部と、
上記抽出された上記目的信号を時間領域の信号に変換する時間領域変換部と、を具備することを特徴とするブラインド信号抽出装置。
【請求項２】
請求項１記載のブラインド信号抽出装置であって、
上記ビームフォーマ計算部は、
上記クラスタのセントロイド情報から上記目的信号のインパルス応答を推定するインパルス応答推定部と、
上記インパルス応答と上記不要信号相関行列を用いて適応型ビームフォーマを計算する適応型ビームフォーマ計算部と、により構成され、
上記目的信号抽出部は、上記適応型ビームフォーマを用い、上記周波数領域の信号から上記目的信号を抽出するものであることを特徴とするブラインド信号抽出装置。
【請求項３】
請求項１記載のブラインド信号抽出装置であって、
上記不要信号相関行列推定部は、上記クラスタの情報から選択されたＫ個の不要信号の相関行列である不要信号相関行列を推定するものであり、ただし、ＫはＫ≦Ｍ−１を満たす整数であり、
更に、上記クラスタの情報から目的信号と上記選択されたＫ個の不要信号のみを含むビームフォーマ入力信号を推定する入力信号推定部を有し、
上記ビームフォーマ計算部は、
上記クラスタのセントロイド情報から上記目的信号のインパルス応答を推定するインパルス応答推定部と、
上記インパルス応答と上記不要信号相関行列とを用いて、適応型ビームフォーマを計算する適応型ビームフォーマ計算部と、により構成され、
上記目的信号抽出部は、上記適応型ビームフォーマを用い、上記ビームフォーマ入力信号から目的信号を抽出するものであることを特徴とするブラインド信号抽出装置。
【請求項４】
請求項２又は３記載のブラインド信号抽出装置であって、
更に、上記Ｍ個のセンサの位置を表すセンサ位置情報を記憶しているセンサ位置情報記憶部を備え、
上記インパルス応答推定部は、
上記センサ位置情報と上記クラスタの上記セントロイド情報とを用いて信号の到来方向を推定する到来方向推定部と、
上記推定された信号の到来方向と上記センサ位置情報からインパルス応答を計算するインパルス応答計算部と、で構成されていることを特徴とするブラインド信号抽出装置。
【請求項５】
請求項１記載のブラインド信号抽出装置であって、
更に、観測信号から観測信号の相関行列である観測信号相関行列を推定する観測信号相関行列推定部を備え、
上記ビームフォーマ計算部は、
上記クラスタの情報から上記目的信号が含まれる観測信号の相関行列である目的信号相関行列を推定する目的信号相関行列推定部と、
上記不要信号相関行列と上記目的信号相関行列とから最大利得ビームフォーマを計算する最大利得ビームフォーマ計算部と、
上記最大利得ビームフォーマに対し、上記観測信号相関行列を用いて、周波数歪みを補正し、補正ビームフォーマを計算する補正部と、を具備し、
上記目的信号抽出部は、上記補正ビームフォーマを用い、上記周波数領域の信号から目的信号を抽出するものである、
ことを特徴とするブラインド信号抽出装置。
【請求項６】
請求項１記載のブラインド信号抽出装置であって、
更に、観測信号から観測信号の相関行列である観測信号相関行列を推定する観測信号相関行列推定部を備え、
上記ビームフォーマ計算部は、
上記クラスタの情報から上記目的信号が含まれる観測信号の相関行列である目的信号相関行列を推定する目的信号相関行列推定部と、
上記不要信号相関行列と上記目的信号相関行列とから最大利得ビームフォーマを計算する最大利得ビームフォーマ計算部とを具備し、
上記観測信号相関行列と上記最大利得ビームフォーマとを用いて、補正ベクトルを計算する補正ベクトル計算部を備え、
上記目的信号抽出部は、
上記最大ビームフォーマを用い、上記周波数領域の信号から歪みを含む目的信号を抽出する信号抽出部と、
上記抽出された歪みを含む目的信号に対し、上記補正ベクトルを用いて、歪み補正をして、上記目的信号を出力する歪み補正部とを具備する
ことを特徴とするブラインド信号抽出装置。
【請求項７】
Ｎ個の信号源から発せられた信号をＭ個のセンサで観測し、観測された信号のうち、１個以上の信号を抽出する信号抽出方法において、ただし、Ｎ、Ｍは２以上の整数であり、
周波数領域変換手段が、上記Ｍ個のセンサで観測された観測信号を周波数領域の信号に変換する周波数領域変換過程と、
正規化手段が、上記周波数領域の信号に対し、正規化を行い、正規化観測信号ベクトルを算出する正規化過程と、
クラスタリング手段が、上記正規化観測信号ベクトルをＮ個のクラスタにクラスタリングするクラスタリング過程と、
不要信号相関行列推定過程が、上記クラスタの情報から、不要信号のみが含まれる観測信号の相関行列である不要信号相関行列を推定する不要信号相関行列推定過程と、
ビームフォーマ計算手段が、上記クラスタの情報と、上記不要信号相関行列と、からビームフォーマを計算するビームフォーマ計算過程と、
目的信号抽出手段が、上記ビームフォーマを用い、上記周波数領域の信号から目的信号を抽出する目的信号抽出過程と、
時間領域変換手段が、上記抽出された上記目的信号を時間領域の信号に変換する時間領域変換過程と、を有することを特徴とするブラインド信号抽出方法。
【請求項８】
請求項７記載のブラインド信号抽出方法であって、
上記ビームフォーマ計算過程は、
インパルス応答推定手段が、上記クラスタのセントロイド情報から上記目的信号のインパルス応答を推定するインパルス応答推定過程と、
適応型ビームフォーマ計算手段が、上記インパルス応答と上記不要信号相関行列を用いて適応型ビームフォーマを計算する適応型ビームフォーマ計算過程と、を有し、
上記目的信号抽出過程は、上記適応型ビームフォーマを用い、上記周波数領域の信号から上記目的信号を抽出する過程であることを特徴とするブラインド信号抽出方法。
【請求項９】
請求項７記載のブラインド信号抽出方法であって、
上記不要信号相関行列推定過程は、上記クラスタの情報から選択されたＫ個の不要信号の相関行列である不要信号相関行列を推定する過程であり、ただし、ＫはＫ≦Ｍ−１を満たす整数であり、
更に、入力信号推定手段が、上記クラスタの情報から目的信号と上記選択されたＫ個の不要信号のみを含むビームフォーマ入力信号を推定する入力信号推定過程を有し、
上記ビームフォーマ計算過程は、
インパルス応答手段が、上記クラスタのセントロイド情報から上記目的信号のインパルス応答を推定するインパルス応答推定過程と、
適応型ビームフォーマ計算手段が、上記インパルス応答と上記不要信号相関行列とを用いて、適応型ビームフォーマを計算する適応型ビームフォーマ計算過程とを有し、
上記目的信号抽出過程は、上記適応型ビームフォーマを用い、上記ビームフォーマ入力信号から目的信号を抽出する過程であることを特徴とするブラインド信号抽出方法。
【請求項１０】
請求項８又は９記載のブラインド信号抽出方法であって、
上記インパルス応答推定過程は、
到来方向推定手段が、センサ位置情報記憶部中の上記センサ位置情報と上記クラスタの上記セントロイド情報とを用いて信号の到来方向を推定する到来方向推定過程と、ここで、センサ位置情報記憶部とは、上記Ｍ個のセンサの位置を表すセンサ位置情報を記憶しているものであり、
インパルス応答計算手段が、上記推定された信号の到来方向と上記センサ位置情報からインパルス応答を計算するインパルス応答計算過程と、を有することを特徴とするブラインド信号抽出方法。
【請求項１１】
請求項７記載のブラインド信号抽出方法であって、
更に、観測信号相関行列推定手段が、観測信号から観測信号の相関行列である観測信号相関行列を推定する観測信号相関行列推定過程を備え、
上記ビームフォーマ計算過程は、
目的信号相関行列推定手段が、上記クラスタの情報から上記目的信号が含まれる観測信号の相関行列である目的信号相関行列を推定する目的信号相関行列推定過程と、
最大利得ビームフォーマ計算手段が、上記不要信号相関行列と上記目的信号相関行列とから最大利得ビームフォーマを計算する最大利得ビームフォーマ計算過程と、
補正手段が、上記最大利得ビームフォーマに対し、上記観測信号相関行列を用いて、周波数歪みを補正し、補正ビームフォーマを計算する補正過程と、を有し、
上記目的信号抽出過程は、上記補正ビームフォーマを用い、上記周波数領域の信号から目的信号を抽出する過程である、
ことを特徴とするブラインド信号抽出方法。
【請求項１２】
請求項７記載のブラインド信号抽出方法であって、
更に、観測信号相関行列推定手段が、観測信号から観測信号の相関行列である観測信号相関行列を推定する観測信号相関行列推定過程を備え、
上記ビームフォーマ計算過程は、
目的信号相関行列推定手段が、上記クラスタの情報から上記目的信号が含まれる観測信号の相関行列である目的信号相関行列を推定する目的信号相関行列推定過程と、
最大利得ビームフォーマ計算手段が、上記不要信号相関行列と上記目的信号相関行列とから最大利得ビームフォーマを計算する最大利得ビームフォーマ計算過程とを有し、
更に、補正ベクトル計算手段が、上記観測信号相関行列と上記最大利得ビームフォーマとを用いて、補正ベクトルを計算する補正ベクトル計算過程を有し、
上記目的信号抽出過程は、
信号抽出手段が、上記最大ビームフォーマを用い、上記周波数領域の信号から歪みを含む目的信号を抽出する信号抽出過程と、
歪み補正手段が、上記抽出された歪みを含む目的信号に対し、上記補正ベクトルを用いて、歪み補正をして、上記目的信号を出力する歪み補正過程とを有する
ことを特徴とするブラインド信号抽出方法。
【請求項１３】
請求項７〜１２の何れかに記載のブラインド信号抽出方法の各過程をコンピュータに実行させるためのブラインド信号抽出プログラム。
【請求項１４】
請求項１３記載のブラインド信号抽出プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、必要である源信号（目的信号）のみを直接観測することが出来ず、目的信号に他のノイズ、干渉信号などが重畳されて観測されるという状況において目的信号を推定して抽出するブラインド信号抽出装置、方法、プログラム、および記録媒体に関するものである。
【背景技術】
【０００２】
ここでは、まず観測信号のモデル化と信号の周波数領域の定義を行い、次に従来技術について簡単に述べる。
［観測信号］
全ての信号はあるサンプリング周波数ｆｓでサンプリングされ、離散的に表現されるものとする。Ｎ個（Ｎは２以上の整数）の源信号が混合されて、Ｍ個（Ｍは２以上の整数）のセンサで観測されたとする。この発明は、信号の発生源からセンサまでの距離により信号が減衰・遅延し、また壁などにより、信号が反射するなどして伝送路歪みが発生しうる状況を扱う。このような状況では、複数の信号源からの源信号ｓｎ（ｔ）（ｎ＝１、．．．、Ｎ）が複数のセンサで観測信号ｘｍ（ｔ）（ｍ＝１、．．．、Ｍ）として観測され、各信号源ｎからセンサｍまでのインパルス応答をｈｍｎ（ｕ）（ｕは時間を表す）とする。センサｍでの観測信号ｘｍ（ｔ）は各源信号ｓｎ（ｔ）に対し、対応するインパルス応答ｈｍｎ（ｕ）が畳込み混合され、次式で表される。
ｘｍ(t)＝Σｎ＝１ＮΣｕ＝０∞ｈｍｎ(ｕ)ｓｎ(ｔ−ｕ)        （１）
ここでは、源信号ｓ１（ｔ）、．．．、ｓＮ（ｔ）やインパルス応答ｈ１１（ｕ）、．．．、ｈ１Ｎ（ｕ）、．．．、ｈＭ１（ｕ）、．．．、ｈＭＮ（ｕ）についての情報を事前に得られない状況を考える。この状況で、観測信号ｘ１（ｔ）、．．．、ｘＭ（ｔ）のみを用いて源信号ｓ１（ｔ）、．．．、ｓＮ（ｔ）を分離抽出することがこの発明の広義の目的である。
［周波数領域表現］
この発明では、周波数領域において、各操作を行う。そのため、センサでの観測信号ｘｍ（ｔ）にＬ点（Ｌは任意の整数）に、公知の技術である例えば、短時間フーリエ変換を適用して周波数ごとの時間系列
ｘｍ(ｆ、τ)＝Σu=-L/2L/2-1ｘｍ(τ＋ｕ)ｇ(ｕ)ｅ−ｉ２πｆｕ        （２）
を求める。ここで、ｆは周波数であり、ｆ＝０、ｆｓ／Ｌ、．．．、ｆｓ（Ｌ−１）／Ｌと離散化されており、τは任意の時間であり、上述の通り、ｆｓはサンプリング周波数である。ｇ（ｕ）は例えばハニング窓などの窓関数である。
【０００３】
式（１）で示される時間領域での畳み込み混合は、周波数領域では、
ｘｍ(ｆ、τ)＝Σｎ＝１Ｎｈｍｎ(ｆ)ｓｎ(ｆ、τ)      （３）
と各周波数での単純混合に近似表現される。ここで、ｈｍｎ(ｆ)は信号源ｎからセンサｍまでの周波数成分ｆについての周波数応答（インパルス応答）、ｓｎ(ｆ、τ)は式（２）と同様の式に従って、源信号ｓｎ（ｔ）に短時間フーリエ変換を施したものであり、以下も同様とする。センサ１〜Ｍの観測信号ｘ１（ｆ、τ）、．．．、ｘＭ(ｆ、τ)を式（３）を用いて、ベクトルで表記すると、
ｘ(ｆ、τ)＝Σｎ＝１Ｎｈｎ(ｆ)ｓｎ(ｆ、τ)          （４）
となる。ここで、ｘ(ｆ、τ)は、ｘ(ｆ、τ)＝［ｘ１(ｆ、τ)、．．．、ｘｍ(ｆ、τ)、．．．、ｘＭ(ｆ、τ)］Ｔとなる観測信号ベクトルであり、ｈｎ(ｆ)は、ｈｎ(ｆ)＝［ｈ1ｎ(ｆ)、．．．、ｈｍｎ(ｆ)、．．．、ｈＭｎ(ｆ)］Ｔであり、信号源から各センサへの周波数応答をまとめたベクトルである。また、［Ａ］ＴはベクトルＡの転置ベクトルを示す。以下の説明も同様とする。
［代表的な従来技術］
混合信号から目的とする信号を抽出する代表的な信号抽出手法として、適応型ビームフォーマ（ａｄａｐｔｉｖｅ  ｂｅａｍｆｏｒｍｅｒ：ＡＢＦ）が非特許文献１等に記載され、広く知られている。
【０００４】
従来の適応型ビームフォーマ（以下、従来型ビームフォーマという）の機能構成例を図１に示す。複数のセンサ４ｍで観測された信号ｘｍ（ｔ）が周波数領域変換部５に入力される。周波数領域変換部５で信号ｘｍ（ｔ）が周波数領域信号ｘｍ（ｆ、τ）に変換される。ｘｍ（ｆ、τ）（ｍ＝１、．．．、Ｍ）は全て従来型ビームフォーマ６へ入力される。
【０００５】
従来型ビームフォーマ６は複数のセンサを用いたシステムにおいて、目的信号ｓｎ（ｔ）を強調し、不要信号ｓ１（ｔ）、．．．、ｓｎ−１（ｔ）、．．．、ｓｎ+１（ｔ）、．．．、ｓＮ（ｔ）をできるだけ抑圧するフィルタｗｎ（ｆ）＝［ｗ１ｎ（ｆ）、．．．、ｗｍｎ（ｆ）、．．．、ｗＭｎ（ｆ）］Ｔを推定することで実現される。
【０００６】
従来型ビームフォーマ６を設計する際には、「目的信号発生源から各センサへのインパルス応答ベクトルｈｎ（ｆ）もしくは、その近似であるステアリングベクトル
ａｎ(ｆ)＝［ｅｘｐ（−ｉ２πｆτ１ｎ），…，ｅｘｐ（−ｉ２πｆτＭｎ)］Ｔ    （５）
が既知である」ということを仮定する。ここでτｍｎは信号源ｎがセンサｍに達する時刻と原点０に達する時間差である。従来は、図２に示すように直線状に配置したセンサシステムを用いることが多く、信号源ｎの方向をθｎ、センサ４ｍのセンサ４１を基準とした座標ｄｍとすると、上述のτｍｎは、τｍｎ＝ｄｍｃｏｓθｎ／ｃで与えられる。ここでｃは信号の速度である。
【０００７】
図１に説明を戻すと、不要信号を抑圧する従来型ビームフォーマ６として、以下の式で表される出力パワーＡ’(ｗｎ(ｆ))を最小にするフィルタ群（ベクトル）ｗｎ(ｆ)を推定する。
Ａ’(ｗｎ(ｆ))＝Ｅ｛｜ｙｎ｜２(ｆ、τ)｝
＝Ｅ｛ｙｎ(ｆ、τ)ｙｎ＊(ｆ、τ)｝
＝Ｅ｛ｗｎＨ(ｆ)ｘ(ｆ、τ）ｘＨ（ｆ、τ）ｗｎ(ｆ)｝
＝ｗｎＨ(ｆ)Ｒｘ(ｆ)ｗｎ(ｆ)                    （６）
ここで、Ｅ｛・｝は時間τに関する平均操作、Ａ＊はＡの複素共役、Ｒｘ（ｆ）＝Ｅ｛ｘ(ｆ、τ）ｘＨ（ｆ、τ）｝は観測信号の相関行列、［Ａ］Ｈは行列（ベクトル）Ａの共役転置行列（ベクトル）を示し、ｙｎ（ｆ、τ）は従来型ビームフォーマ６の出力であり、以下の式（７）で表すことが出来る。
ｙｎ(ｆ、τ)＝ｗｎＨ(ｆ)ｘ(ｆ、τ)      （７）
ここで、意味のない解(ｗｎ(ｆ)＝０＝［０、．．．、０］Ｔ)を回避するために、目的信号が無歪みで得られるという以下の式に示す拘束条件を付与する。
ｗｎＨ(ｆ)ｈｎ(ｆ)＝１                                    （８）
これにより、式（８）を満たし、かつ上記式（６）のＡ’(ｗｎ(ｆ))の値が最小となるｗｎ(ｆ)の値を求める問題はＬａｇｕｒａｎｇｅの未定乗数ｐを用いて、以下の式（９）で表すことができる。
Ａ(ｗｎ(ｆ))＝Ａ’(ｗｎ(ｆ))＋ｐ(ｗｎＨ(ｆ)ｈｎ(ｆ)−１)    （９）
式（９）を解くことにより、従来型ビームフォーマ６は、
【数１】
2008060635.tif 000003
で得られる。
【０００８】
従来型ビームフォーマ６（従来の適応型ビームフォーマ）では、式（１０）におけるインパルス応答ベクトルｈｎ（ｆ）は実測してインパルス応答記憶部１０に記憶させておいたインパルス応答ベクトルｈｎ（ｆ）を読み出して用いることが理想である。しかし、代わりに上記式（５）に示すステアリングベクトルａｎ（ｆ）をステアリングベクトル記憶部１２に記憶させ、読み出されたステアリングベクトルａｎ（ｆ）を上記式（１０）のインパルス応答ベクトルｈｎ（ｆ）の代わりに用いることが広く行われている。
【０００９】
しかし、実環境において、インパルス応答ベクトルｈｎ（ｆ）やステアリングベクトルａｎ（ｆ）が正しく与えられることは稀であり、上記式（９）に示すＡ(ｗｎ(ｆ))の最小化が必ずしも不要信号のみの最小化にはならないことが多い。このことから、混合信号（観測信号）の相関行列Ｒｘ（ｆ）の代わりに、不要信号のみの時間区間における信号ξ（ｆ、τ）の相関行列ＲＪ（ｆ）＝Ｅ｛ξ（ｆ、τ）ξＨ（ｆ、τ）｝を用いることが非特許文献２などで広く行われている。これは、Ｒｘ（ｆ）を用いる場合よりも高い性能を実現することが知られている。即ち、従来型ビームフォーマにおいては、不要信号のみの時間区間（目的音不在の時間区間）における相関行列が精度よく推定できることが望ましい。
【００１０】
従来型ビームフォーマ６は上記式（１０）のｗｎ（ｆ）と観測信号ベクトルｘ（ｆ、τ）により、上記式（７）により、出力信号ベクトルｙｎ（ｆ、τ）が出力される。出力信号ｙｎ（ｆ、τ）は時間領域変換部８に入力され、周波数領域から時間領域に変換され、ｙｎ（ｔ）が生成される。
【非特許文献１】Ｈａｙｋｉｎ，Ｓ．  適応フィルタ理論  科学技術出版  ２００１  ６９０頁−６９３頁
【非特許文献２】大賀寿郎  山崎芳男  金田豊  音響システムとディジタル処理  電子情報通信学会編  コロナ社  １９０頁−１９１頁
【発明の開示】
【発明が解決しようとする課題】
【００１１】
上述の通り、従来型ビームフォーマでは目的信号源から各センサへのインパルス応答ベクトルｈｎ（ｆ）もしくはその近似であるステアリングベクトルａｎ（ｆ）が必要である。すなわち目的信号に関する事前知識が必要であるという難点がある。更にそれらは、実環境では正しく得ることが困難であり、事前知識と使用環境でのインパルス応答ベクトルｈｎ（ｆ）がずれてしまった場合の従来型ビームフォーマの性能は著しく低下する。
また、高い性能を得るためには、不要信号のみの時間区間における信号の相関行列ＲＪ（ｆ）を推定する必要があるが、不要信号が非定常な信号である場合には、それは非常に困難である。
【課題を解決するための手段】
【００１２】
Ｎ個の信号源から発せられた信号をＭ個のセンサで観測し、観測された信号のうち、１個以上の信号を抽出する信号抽出装置において、ただし、Ｎ、Ｍは２以上の整数であり、上記Ｍ個のセンサで観測された観測信号を周波数領域の信号に変換し、上記周波数領域の信号に対し、正規化を行い、正規化観測信号ベクトルを算出し、上記正規化観測信号ベクトルを上記Ｎ個のクラスタにクラスタリングし、
上記クラスタの情報から、不要信号のみが含まれる観測信号の相関行列である不要信号相関行列を推定し、上記クラスタの情報と、上記不要信号相関行列と、からビームフォーマを計算し、上記ビームフォーマを用い、上記周波数領域の信号から上記目的信号を抽出し、上記抽出された上記目的信号を時間領域の信号に変換する。
【発明の効果】
【００１３】
上記の構成により、事前にインパルス応答ベクトルあるいは、ステアリングベクトルを測定しておくことなく、また不要信号が非定常な信号であっても、目的信号を精度よく、分離抽出することが可能となる。
【発明を実施するための最良の形態】
【００１４】
以下に発明を実施するための最良の形態を示す。
【実施例１】
【００１５】
この発明の機能構成例を図３に示し、この発明の主要な処理の流れを図４に示す。図１と同一機能構成部分には同一参照符号を付け、重複説明を省略する。以下も同様とする。
【００１６】
また、この発明では、信号のスパース性を仮定する。スパースとは、信号が殆どの時刻τにおいて、０であることを示す。信号のスパース性は、例えば、音声信号で確認される。信号のスパース性を仮定することで、複数の源信号が存在しても、各時間周波数ポイント（ｆ、τ）では互いに重ならず、高々１つしか存在しないと仮定することが出来る。即ち上記式（４）は以下の式で表すことができる。
【００１７】
ｘ(ｆ、τ)＝ｈｎ(ｆ)ｓｎ(ｆ、τ)      （１１）
ここで、ｈｎ(ｆ)はインパルス応答ベクトル、ｓｎ(ｆ、τ)は（ｆ、τ）に存在する源信号を表す。
【００１８】
センサ４ｍで収音されたそれぞれの観測信号ｘｍ（ｔ）（ｍ＝１、．．．Ｍ）は周波数領域変換部５に入力される。周波数領域変換部５でそれぞれの観測信号ｘｍ（ｔ）は、例えば、公知の技術である上記の短時間フーリエ変換などで、時間領域から周波数領域に変換され、ｘｍ（ｆ、τ）に変換される（ステップＳ２）。更に、ｘｍ（ｆ、τ）は、観測信号ベクトルｘ（ｆ、τ）として出力される。観測信号ベクトルｘ（ｆ、τ）は正規化部２２に入力され、正規化観測信号ベクトルが算出される（ステップＳ４）。具体的には、観測信号ベクトルｘ（ｆ、τ）＝［ｘ１（ｆ、τ）、．．．、ｘｍ（ｆ、τ）］Ｔに対し、偏角の正規化を以下の式で行う。
【数２】
2008060635.tif 000004
また、ノルム正規化を以下の式で行う。
【００１９】
ｘ−（ｆ、τ）←ｘ−（ｆ、τ）／‖ｘ−（ｆ、τ）‖      （１３）
ここで、ｘ−（ｆ、τ）は正規化された観測信号ベクトルｘ（ｆ、τ）を表し、ａｒｇ（ｒ）はｒの偏角を表し、ｉは虚数単位を表し、│ｒ│はｒの絶対値を表し、‖ｒ‖はｒのノルムを表し、Ｑは基準とするセンサの番号（Ｑ∈｛1、．．．、Ｍ｝）を表し、ｃは信号の速度を表し、αは任意の正の定数を表す。αについては、α＝４ｄｍａｘが最も好ましい。ただし、ｄｍａｘは、基準として選択された任意のセンサＱと他のセンサとの距離の最大値を表す。また、αは他の数値でもよい。
【００２０】
上記式（１１）〜（１３）より、正規化された観測信号ベクトルｘ−（ｆ、τ）は以下の式で表すことができる。
【数３】
2008060635.tif 000005
ここで、Ａｎ＝（Σｍ＝１Ｍ│ｈｍｎ│２）１／２であり、信号ｓｎ（ｆ、τ）に関するインパルス応答情報にのみ依存することが分かる。
【００２１】
正規化された全ての時間周波数の観測信号ベクトルｘ−（ｆ、τ）はクラスタリング部２４に入力され、Ｎ個のクラスタにクラスタリングされる（ステップＳ６）。このクラスタリングは、例えば、ｋ−ｍｅａｎｓ法を用いて効果的に行うことができる。また、詳細は「Ｒ．Ｏ．Ｄｕｄａ、Ｐ．Ｅ．Ｈａｒｔ、ａｎｄ  Ｄ．Ｇ．Ｓｔｏｒｋ，Ｐａｔｔｅｒｎ  Ｃｌａｓｓｉｆｉｃａｔｉｏｎ，Ｗｉｌｅｙ  Ｉｎｔｅｒｓｃｉｅｎｃｅ,  ２ｎｄ ｅｄｉｔｉｏｎ，２０００．」に記載されている。以下にクラスタリングの方法を具体的に説明する。  
【００２２】
記憶部２６にはあらかじめ、上記式（１４）で示される正規化観測信号ベクトルｘ−（ｆ、τ）とセントロイドの初期値ｃｊｎ（ｊ＝０、ｎ＝１、...、Ｎ）が記憶されている。
クラスタリング部２４は、記憶部２６から正規化観測信号ベクトルｘ−（ｆ、τ）を読み込み、これらをクラスタリングしてＮ個のクラスタＣ１、...、ＣＮを生成する。すなわち、Ｍ次元複素ベクトルである正規化された観測信号ベクトルｘ−（ｆ、τ）をＭ次元複素空間で以下の手順で直接クラスタリングする。
【００２３】
１．クラスタのセントロイドの初期値ｃｊｎを記憶部２６から読み込む。セントロイドの初期値ｃｊｎは、正規化観測信号ベクトルｘ−（ｆ、τ）と同じ次元のベクトル（Ｍ次元複素ベクトル）である。なお、セントロイドの初期値ｃ０ｎの選び方については後述する。
２．ｊ＋１を新たなｊとする。
３．すべての時間周波数（ｆ、τ）における正規化観測信号ベクトルｘ−（ｆ、τ）を、最も近いセントロイドｃｊ−１ｎで代表されるクラスタＣｎに割り当てる。すなわち、各正規化ベクトルｘ−（ｆ、τ）に対して、‖ｘ−（ｆ、τ）−ｃｊ−１ｎ‖が最も小さくなるようにｎを選ぶ。
４．各クラスタＣｎに割りあてられた正規化観測信号ベクトルｘ−（ｆ、τ）の平均値を計算し、そのノルムを１にすることでセントロイドを更新する。すなわち、各クラスタＣｎに割りあてられた正規化観測信号ベクトルｘ−（ｆ、τ）に対して、
ｃｎｊ＝Ｅ｛ｘ（ｆ、τ）｝ｎ／‖Ｅ｛ｘ−（ｆ、τ）｝ｎ‖    （１５）
の演算を行うことにより、セントロイドを更新する。ここで、Ｅ｛・｝ｎは、クラスタＣｎのメンバに対する平均操作を表す。
５．セントロイドｃｊｎが収束するまで、手順２−５を繰り返す。最後に収束したセントロイドを、ｃｎ（ｎ＝１、...、Ｎ）として、記憶部２６に記憶する。以上が、クラスタリング手順である。
【００２４】
次に、セントロイドの初期値の選び方の例を説明する。
《初期値設定方法１》
正規化観測信号ベクトルｘ−（ｆ、τ）の中からランダムにＮ個のベクトルを選び、それをセントロイドの初期値ｃ０ｎ（ｎ＝１、...、Ｎ）とする。
《初期値設定方法２》
セントロイドは、式（１１）〜（１５）において、全てのｍ、ｎについて、│ｈｍｎ（ｆ）│＝１と仮定すると、以下の式（１６）のように書けるので、これを用いる。
｛ｃｎ｝ｑ＝Ｅ［ｘ−（ｆ、τ）］ｎ
＝ｅｘｐ［ｉ２π（ｄｍ−ｄＱ）Ｔｖｎ／α］／Ｍ１／２ 
＝ｅｘｐ［ｉ２π‖ｄｍ−ｄＱ‖ＴｃｏｓΘｎｍＱ／α］／Ｍ１／２ 
（１６）
ここで、ｄｍはセンサ４ｍの位置ベクトルを表し、ｖｎ＝ｃｏｓΘｎｍＱは、センサ４ｍと基準として選択したセンサ４Ｑを結ぶ軸に対する信号ｓｎ（ｔ）の到来方向ベクトルを表し、図５において、太いベクトルで示されているものである。また、ｖｎは単位ベクトルであり、‖ｖｎ‖＝１である。
【００２５】
センサ位置ｄｍ（ｍ＝１、...、Ｍ）は、記憶部２６において保持されている値を、方位θｎと仰角φｎ（ｎ＝１、...、Ｎ）を適当に与える。ここで、センサ位置ｄｍ、方位θｎと仰角φｎは初期値であるので、適当な値で良い。例えば、θｎ＝２πｎ／Ｎ、φｎ＝０とすると、空間的に散らばった初期値が得られる。
【００２６】
図３に説明を戻すと、クラスタリング部２４により求まったクラスタＣｎは各源信号ｓｎ（ｆ、τ）に対応している。また、そのセントロイドｃ−ｎ＝Ｅ｛ｘ−（ｆ、τ）｝ｘ−（ｆ、τ）∈Ｃｎは上記式（１４）から理解されるように、源信号ｓｎ（ｆ、τ）に関するインパルス応答情報を表すことが分かる。
【００２７】
各クラスタＣｎはそれぞれ、不要信号相関行列推定部２５に入力される。不要信号相関行列推定部２５では各クラスタＣｎごとにその情報から、源信号ｓｎ（ｆ、τ）に対する不要信号区間の相関行列、つまり、不要信号のみが含まれる観測信号の相関行列である不要信号相関行列ＲｎＪ（ｆ）を以下の式（１７）（１８）で推定する（ステップＳ８）。
【数４】
2008060635.tif 000006
推定された不要信号相関行列ＲｎＪ（ｆ）とクラスタリング部２４よりのクラスタのセントロイド情報ｃ―ｎは、ビームフォーマ計算部２８に入力される。ビームフォーマ計算部２８ではクラスタの情報Ｃｎと、不要信号相関行列ＲｎＪ（ｆ）とからビームフォーマｗｎ（ｆ）を計算する（ステップＳ１０）。ビームフォーマｗｎ（ｆ）の具体的な計算方法については、実施例２で詳細に説明する。
【００２８】
計算されたビームフォーマｗｎ（ｆ）は目的信号抽出部３０に入力される。目的信号抽出部３０では、ビームフォーマｗｎ（ｆ）を用いて、以下の式（１９）を計算して、周波数領域の観測信号ｘ（ｆ、τ）から目的信号ｙｎ（ｆ、τ）を抽出する（ステップＳ１２）。
ｙｎ(ｆ、τ)＝ｗｎ(ｆ)Ｈｘ(ｆ、τ)              （１９）
式（１７）〜（１９）を全てのｎ（ｎ＝１、．．．、Ｎ）に対して行うことで、Ｎ個全ての信号を抽出する。ｙｎ(ｆ、τ)は全て、時間領域変換部３２に入力される。目的信号抽出部３０で抽出された目的信号ｙｎ（ｆ、τ）は、例えば公知の技術である短時間逆フーリエ変換などで時間領域の信号ｙｎ（ｔ）に変換される（ステップＳ１４）。
【実施例２】
【００２９】
次に、この発明の実施例２について説明する。実施例２は、実施例１で説明したビームフォーマ計算部２８をより詳細に構成した例である。図６に実施例２のビームフォーマ計算部２８とこれに関係する部分の機能構成例を示す。図６に図示されていない部分は実施例１で説明したものと同様の処理を行うものとし、以下の実施例についても同様とする。ビームフォーマ計算部２８はインパルス応答推定部４０と適応型ビームフォーマ計算部４２とで構成されている。
【００３０】
インパルス応答推定部４０では、クラスタＣｎのセントロイド情報ｃ―ｎから目的信号のインパルス応答ベクトルｈｎ（ｆ）を推定する。具体的には、クラスタリング部２４よりのセントロイド情報ｃ―ｎについて逆正規化を行うことで源信号ｓｎ（ｆ、τ）に対するインパルス応答ベクトルの推定を行う。
【００３１】
まず、クラスタＣｎのセントロイドは上記式（１５）であり、ｘ−（ｆ、τ）は式（１４）で表すことができる。ここで、全てのｍ、ｎに対して、│ｈｍｎ│＝１と仮定するとセントロイドｃ―ｎのｍ番目の成分ｃ―ｍｎは、以下の式（２０）が成り立つ。
【数５】
2008060635.tif 000007
式（２０）を、ｈｍｎ（ｆ）について解くと、以下の式（２１）を得ることが出来る。
【００３２】
【数６】
2008060635.tif 000008
式（２０）の右辺はインパルス応答ｈｍｎ（ｆ）が正規化されたものになっているが、式（２１）は式（２０）から逆にインパルス応答ベクトルｈｎ（ｆ）について求め直しているので、式（２１）は逆正規化と呼ぶ。
【００３３】
図６の説明に戻ると、推定されたインパルス応答ベクトルｈｎ（ｆ）と、不要信号相関行列推定部２５よりのＲｎＪ（ｆ）とが、適応型ビームフォーマ計算部４２に入力される。
【００３４】
適応型ビームフォーマ計算部４２は、インパルス応答ベクトルｈｎ（ｆ）と上記不要信号相関行列ＲｎＪ（ｆ）を用いて適応型ビームフォーマｗｎ（ｆ）を計算する。具体的には、以下の式（２２）により適応型ビームフォーマｗｎ（ｆ）を計算することが出来る。
【数７】
2008060635.tif 000009
上記式（２２）は、上記式（１０）のＲｘ（ｆ）をＲｎＪ（ｆ）に置き換えることで得ることができる。
【００３５】
目的信号抽出部３０では、式（２２）の適応型ビームフォーマｗｎ（ｆ）と、周波数領域の信号ｘ（ｆ、τ）とを上記式（１９）に適用して、目的信号ｙｎ（ｆ、τ）を抽出する。
【００３６】
次に実施例２の変形例として実施例２’を示す。インパルス応答推定部４０による、上記式（２２）を用いたインパルス応答ベクトルｈｎ（ｆ）の推定を出力する代わりに、ステアリングベクトルａｎ（ｆ）を推定して、出力させることも考えられる。信号ｓｎ（ｆ、τ）のステアリングベクトルを上記式（５）と同様に
ａｎ（ｆ）＝［ｅｘｐ（−ｉ２πｆτ１ｎ）、・・・、ｅｘｐ（−ｉ２πｆτＭｎ）］Ｔ    (２３)
とすると、ステアリングベクトルａｎ（ｆ）はインパルス応答ベクトルｈｎ（ｆ）の推定であるから、上記式（２１）と上記式（２３）の位相項を比較すると、τｍｎ（ｍ＝１、．．．、Ｍ）は以下の式（２４）で推定できる。
τ＾ｍｎ＝αｃ−１ａｒｇ［ｃ−ｍｎｃ−Ｑｎ］／２π      （２４）
この式（２４）の計算をインパルス応答推定部４０’で行う。
【００３７】
τ＾ｍｎを用いたステアリングベクトルａｎ（ｆ）をインパルス応答ベクトルｈｎ（ｆ）の推定として出力する。すなわち上記式（２３）（２４）からステアリングベクトルａｎ（ｆ）は以下の式（２７）で表すことができ、インパルス応答推定部４０’からインパルス応答ｈ＾ｎ（ｆ）として出力される。
ａｎ（ｆ）＝［ｅｘｐ（−ｉ２πｆτ＾１ｎ）、・・・、ｅｘｐ（−ｉ２πｆτ＾Ｍｎ）］Ｔ≒ｈ＾ｎ（ｆ）    （２５）
式（２５）のｈ＾ｎ（ｆ）とＲｎＪから、上記式（２２）で適応型ビームフォーマを計算する。
【００３８】
また、上記式（２２）を使用して不要信号相関行列推定部２５で適応型ビームフォーマｗｎ（ｆ）を推定する際に、不要信号相関行列ＲｎＪ（ｆ）を使用して、音響伝達特性つまり、インパルス応答ベクトルｈｎ（ｆ）やステアリングベクトルａｎ（ｆ）が既知の場合はそれらを使用しても良い。また、式（２２）中の不要信号相関行列ＲｎＪ（ｆ）の代わりに、観測信号の相関行列Ｒｘ（ｆ）を使用しても良い。これらのことは、以下の実施例においても同様である。
【実施例３】
【００３９】
実施例２で示した適応型ビームフォーマでは、Ｎ≦Ｍの場合には、高い性能を得ることが出来るが、Ｎ＞Ｍの場合には、性能が限られることが問題であった。具体的には、適応型ビームフォーマは不要信号の数がＭ−１個以下であれば、効果的に目的信号ｙｎ（ｆ、τ）を抽出できるが、Ｍ個以上であると、その効果が不十分であることが知られている。そこで実施例３では、Ｎ＞Ｍの場合、即ち、Ｎ−１（＞Ｍ−１）個の不要信号がある場合にも目的信号を抽出できることを示す。実施例３の機能構成例を図７に示す。実施例３は、実施例２と比較して、不要信号選択部４９、入力信号推定部５０が追加され、不要信号相関行列推定部２５、インパルス応答推定部４０、適応型ビームフォーマ計算部４２、の処理が変更されている。
【００４０】
不要信号選択部４９で、Ｋ個の不要信号についての不要信号相関行列ＲｎＪ（ｆ）を推定する。ここで、Ｋは、Ｋ≦Ｍ−１を満たす整数とする。つまり、目的信号ｓｎ（ｆ、τ）に相当するクラスタＣｎ以外の不要信号に相当するクラスタＣＬ（Ｌ＝１、．．．、ｎ−１、ｎ＋１、．．．、Ｎ）からＫ個のクラスタを選び、これらのクラスタをＣＪとする。
Ｋ個のクラスタの選び方として、クラスタＣＬ中で、クラスタメンバが多いものから順にＫ個のクラスタを選ぶ方法や、以下の式（２６）で表されるξＬ（ｆ、τ）のパワーが大きいものから順にＫ個のクラスタを選ぶ方法等が考えられる。
【数８】
2008060635.tif 000010
【００４１】
不要信号選択部４９で選択されたＫ個のクラスタＣＪを用いて、Ｋ個の不要信号についての不要信号相関行列ＲｎＪ（ｆ）を以下の式（２７）（２８）で計算する。
【数９】
2008060635.tif 000011
【００４２】
また、入力信号推定部５０において、不要信号選択部４９で選択したＫ個の不要信号と目的信号Ｃｎとが混合したビームフォーマ入力信号ベクトルｘ（ｆ、τ）を推定する。これは、不要信号クラスタＣＪと目的信号クラスタＣｎを用いて、以下の式（２９）で得ることができる。
【数１０】
2008060635.tif 000012
適応型ビームフォーマ計算部４２では、上記式（２７）（２８）で得られた不要信号相関行列ＲｎＪ（ｆ）とインパルス応答推定部４０または４０’よりのインパルス応答ベクトルｈｎ（ｆ）を用いて、上記式（２２）を用いて、適応型ビームフォーマｗｎ（ｆ）を計算する。
【００４３】
目的信号抽出部３０では、入力信号推定部５０よりのビームフォーマ入力信号ｘ（ｆ、τ）が入力される。目的信号抽出部３０では適応型ビームフォーマ計算部４２よりの適応型ビームフォーマｗｎ（ｆ）と、ビームフォーマ入力信号ベクトルｘ（ｆ、τ）を用いて、上記式（１９）を計算して、目的信号ｙｎ（ｆ、τ）を抽出する。
【００４４】
実施例３では上述の通り、Ｎ＞Ｍの場合を説明した。しかし、Ｎ≦Ｍの場合でも実施できる。この場合は、実施例２と比較して、入力信号推定部５０と不要信号選択部４９の処理の分だけ余計にコストがかかる。
【実施例４】
【００４５】
実施例４は、センサの位置情報が既知の場合の実施例である。実施例４の機能構成例の一部を図８に示す。実施例２もしくは３で説明したインパルス応答推定部４０は到来方向推定部６０とインパルス応答計算部６２により構成されている。また、Ｍ個のセンサの位置を表すセンサ位置情報を記憶しているセンサ位置情報記憶部６４を有する。
【００４６】
到来方向推定部６０で、信号の到来方向を推定する。到来方向推定部６０には、クラスタリング部２４よりのセントロイド情報ｃ−ｎとセンサ位置情報記憶部６４よりの各センサの位置を表す３次元ベクトルｄｍ（ｍ＝１、．．．、Ｍ）が入力される。信号ｓｎの到来方向を表す長さ１の３次元ベクトルをｖｎ（ｎ＝１、．．．、Ｎ）とすると、信号ｓｎの到来方向の推定値はセントロイド情報ｃ−ｎを用いて、以下の式（３０）で計算することができる。
【００４７】
ｖｎ＝αＤ＋ａｒｇ［ｃ−ｎ］／２π      （３０）
ここで、Ｄ＝［ｄ１−ｄＱ、．．．、ｄｍ−ｄＱ、．．．、ｄＭ−ｄＱ］Ｔであり、ｄＱは基準として、任意に選択したセンサ４Ｑの位置を表す３次元ベクトルであり、Ｄ＋は、Ｄの一般化逆行列を表す。
【００４８】
次にインパルス応答計算部６２において、信号ｓｎの到来方向とセンサ位置情報とを用いて、インパルス応答を計算する。インパルス応答計算部６２には、到来方向推定部６０よりの信号ｓｎの到来方向の推定値ｑｎと、センサ位置情報記憶部６４よりのセンサ位置情報ｄＱが入力される。インパルス応答計算部６２は以下の式（３１）で表す信号ｓｎについてのステアリングベクトルａｎ＾（ｆ）の推定値を求める。このステアリングベクトルａｎ＾（ｆ）をインパルス応答ベクトルの推定値ｈｎ（ｆ）として計算する。
【数１１】
2008060635.tif 000013
そして、ステアリングベクトルａｎ＾（ｆ）（インパルス応答ベクトルの推定値ｈｎ（ｆ））がインパルス応答計算部６２から出力され、適応型ビームフォーマ計算部４２に入力される。
【実施例５】
【００４９】
本実施例では、適応型ビームフォーマの代わりに、最大利得ビームフォーマを用いる構成を示す。最大利得ビームフォーマとは、センサアレイ出力における目的信号を最大にしつつセンサアレイ出力における不要信号成分を最小にするようなフィルタｗｎ（ｆ）をビームフォーマとする方法である。（Ｄ．Ｈ．Ｊｏｈｎｓｏｎ  ａｎｄ  Ｄ．Ｅ．Ｄｕｄｇｅｏｎ，“Ａｒｒａｙ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｃｏｎｃｅｐｔｓ  ａｎｄ  Ｔｅｃｈｎｉｑｕｅｓ”，Ｐｒｅｎｔｉｃｅ  Ｈａｌｌ，１９９３．）
最大利得ビームフォーマにおいては、センサアレイ出力中の目的信号成分と不要信号成分を推定することが１つのポイントとなるが、不要信号が非定常信号である場合、不要信号を推定することは非常に困難であるという問題があった。実施例５では、この問題をスパース性の仮定を用いることで解決する。つまり（１）目的信号のみの観測信号の相関行列である目的信号相関行列ＲＴｎ（ｆ）と、不要信号のみの観測信号の相関行列である不要信号相関行列ＲＪｎ（ｆ）とを推定すること、（２）最大利得ビームフォーマ計算部で、目的信号相関行列ＲＴｎ（ｆ）と不要信号相関行列ＲＪｎ（ｆ）とから最大利得ビームフォーマｗｎ（ｆ）を推定することにより解決することが出来る。
【００５０】
また、最大利得ビームフォーマでは「目的信号の歪みを最小にする」という上記式拘束条件（８）が無いため、各周波数ｆにおいて、様々な、ゲイン特性を持つビームフォーマｗｎ（ｆ）が構成される。これは、例えば、音声信号のような広帯域信号に最大利得ビームフォーマを適用した場合、出力がｗｎ（ｆ）の周波数特性により歪んでしまうことを意味する。このため、従来は最大利得ビームフォーマを広帯域信号に用いることは困難であった。実施例５では、観測信号ベクトルｘ（ｆ、τ）と最大利得ビームフォーマｗｎ（ｆ）の出力信号との誤差が最小となるように、最大利得ビームフォーマｗｎ（ｆ）を補正することでこれを解決する。
【００５１】
まず、最大利得ビームフォーマの原理を簡単に説明する。上述したように「センサアレイ出力における目的信号を最大にしつつセンサアレイ出力における不要信号成分を最小にする」との条件より評価関数は以下の式（３２）となる。
【数１２】
2008060635.tif 000014
【００５２】
ここで、分母は不要信号の出力パワー、分子は目的信号の出力パワーであり、ＲＴｎ（ｆ）は目的信号のみの観測信号の相関行列、ＲＪｎ（ｆ）は不要信号のみの観測信号の相関行列である。また、（ＲＪｎ（ｆ））１／２＝ＥＦ１／２ＥＨで表すことが出来、ここで、Ｅ＝［ｅ１、．．．ｅｍ］であり、ｅｉはＲＪｎ（ｆ）の固有ベクトルであり、Ｆ＝ｄｉａｇ（λ１、．．．、λＭ）であり、λｉはｅｉに対応するＲＪｎの固有値とし、ｗ〜＝（ＲＪｎ（ｆ））１／２ｗｎとすると、上記式（３２）は以下の式（３３）に変えることが出来る。
【数１３】
2008060635.tif 000015
【００５３】
ここで「児玉、須田、“システム制御のためのマトリクス理論、コロナ社、１９９５」に記載のレイリー商の定理より、ｇ（ｗ〜）の最大値は、（ＲＪｎ（ｆ））−１／２（ＲＴｎ（ｆ））（ＲＪｎ（ｆ））−１／２の最大固有値λで与えられ、対応する固有ベクトルをｅとすると、最大値はｍａｘｇ（ｗ〜）＝λ＝ｇ（ｅ）となる。すなわち求める最大利得ビームフォーマｗｎは以下の式（３４）（３５）で表すことができる。
ｗ〜＝ｅ      （３４）
ｗｎ＝（ＲＪｎ（ｆ））−１／２ｅ    （３５）
実施例５の機能構成例を図９に示す。実施例１と比較して、観測信号相関行列推定部７２が追加され、ビームフォーマ計算部２８は、目的信号相関行列推定部７０、固有ベクトル計算部７４、最大利得ビームフォーマ計算部７６、補正ベクトル計算部７８、補正部８０、とで構成されている。
【００５４】
目的信号相関行列推定部７０で、クラスタの情報から目的信号ｓｎ（ｆ、τ）のみの時間区間の相関行列を以下の式（３６）（３７）で推定する。
【数１４】
2008060635.tif 000016
ここで、Ｃｎは目的信号に対応するクラスタである。不要信号相関行列推定部２５よりの不要信号相関行列ＲＪｎ（ｆ）と、目的信号相関行列ＲＴｎ（ｆ）とが、固有ベクトル計算部７４に入力される。固有ベクトル計算部７４で（ＲＪｎ（ｆ））−１／２（ＲＴｎ（ｆ））（ＲＪｎ（ｆ））−１／２の最大固有ベクトルｅｎ（ｆ）を上記で説明したレイリー商の定理より、計算する。
【００５５】
不要信号相関行列推定部２５よりのＲＪｎ（ｆ）と固有ベクトル計算部７４よりのｅｎ（ｆ）とが最大利得ビームフォーマ計算部７６に入力される。最大利得ビームフォーマ計算部７６では、以下の式（３８）より最大利得ビームフォーマｗｎ(ｆ)を計算する。
ｗｎ(ｆ)＝（ＲＪｎ(ｆ)）−1/2ｅｎ(ｆ)              （３８）
この式（３８）は、上記式（３５）に基づいている。
【００５６】
一方、観測信号相関行列推定部７２で、観測信号ベクトルｘ（ｆ、τ）の相関行列である観測信号相関行列Ｒｘ（ｆ）を以下の式（３９）を用いて推定する。
Ｒｘ（ｆ）＝Ｅ｛ｘ（ｆ、τ）ｘＨ（ｆ、τ）｝      （３９）
補正ベクトル計算部７８に、最大利得ビームフォーマ計算部７６よりの最大利得ビームフォーマｗｎ（ｆ）と、観測信号相関行列推定部７２よりの観測信号相関行列Ｒｘ（ｆ）が入力される。補正ベクトル計算部７８では、最大利得ビームフォーマｗｎ（ｆ）を補正するための補正ベクトルαｎ（ｆ）を生成する。この補正は、最大利得ビームフォーマｗｎ（ｆ）が出力に与える歪みが最小になるよう最大利得ビームフォーマｗｎ（ｆ）を変換する。例えば、以下の式（４０）であらわされる観測信号ベクトルｘ（ｆ、τ）と出力信号ベクトルｙｎ（ｆ、τ）との誤差Ａを最小にする補正ベクトルαｎ（ｆ）を計算する。
Ａ(αｎ(ｆ))＝Ｅ｛‖ｘ(ｆ、τ)−αｎ(ｆ)ｙｎ(ｆ、τ)‖２｝    （４０）
ここで、ｙｎ(ｆ、τ)は最大利得ビームフォーマｗｎ（ｆ）の出力ｙｎ（ｆ、τ）＝ｗｎ（ｆ）ｘ（ｆ、τ）である。
上記式（４０）の右辺を展開すると、
Ａ(αｎ(ｆ))＝｛Ｅ［‖ｘ(ｆ、τ)‖］｝２−αｎ(ｆ)Ｅ［ｘＨ（ｆ、τ）ｙｎ(ｆ、τ)］−αｎＨ(ｆ)Ｅ［ｙｎ(ｆ、τ)＊ｘ（ｆ、τ）］
＋αｎαｎＨＥ［│ｙｎ(ｆ、τ)│２］    （４１）
式（４１）において、両辺をαｎＨ(ｆ)で偏微分すると、以下の式（４２）になる。
∂Ａ(αｎ(ｆ))／∂  αｎＨ(ｆ)＝
−Ｅ［ｙｎ(ｆ、τ)＊ｘ（ｆ、τ）］＋αｎＥ［│ｙｎ(ｆ、τ)│２］    （４２）
上記式（４２）の左辺を０とおき、αｎについて求めると、以下の式（４３）になる。
αｎ(ｆ)＝Ｅ［ｙｎ(ｆ、τ)＊ｘ（ｆ、τ）］／Ｅ［│ｙｎ(ｆ、τ)│２］    
（４３）
ここで、上記式（１９）と上記式（３９）より上記式（４３）は以下の式（４４）になる。
【００５７】
【数１５】
2008060635.tif 000017
ここで、上述したように、Ｒｘ（ｆ）は観測信号ベクトルｘ（ｆ、τ）の相関行列である。上記式（４４）から理解されるように、最大利得ビームフォーマｗｎ（ｆ）と観測信号ベクトルｘ（ｆ、τ）を用いて、補正ベクトル計算部７８では、補正ベクトルαｎ(ｆ)が計算される。
【００５８】
補正部８０は、最大利得ビームフォーマｗｎ（ｆ）に対し、補正ベクトルαｎ(ｆ)を用いて、周波数歪みを補正し、補正ビームフォーマを計算する。具体的には以下の式（４５）により補正して補正ビームフォーマｗｎ’(ｆ)を求めることが出来る。
ｗｎ’(ｆ)＝［αｎ(ｆ)］Ｂｗｎ(ｆ)    （４５）
ここで、Ｂは任意のセンサの番号であり、Ｂ∈｛１、．．．、Ｍ｝であり、［ｑ］ＢはベクトルｑのＢ番目の要素であることを示している。
【００５９】
目的信号抽出部３０では、補正ビームフォーマｗｎ’(ｆ)を用いて、以下の式（４６）で目的信号ｙｎ（ｆ、τ）を抽出する。
ｙｎ（ｆ、τ）＝ｗｎ’Ｈ(ｆ)ｘ（ｆ、τ）      （４６）
また、実施例５の変形例の機能構成例を図１０に示す。ビームフォーマ計算部２８は目的信号相関行列推定部７０、固有ベクトル計算部７４、最大利得ビームフォーマ計算部７６、とで構成され、目的信号抽出部３０は信号抽出部８１と歪み補正部８２とで構成されている。
【００６０】
最大利得ビームフォーマ計算部７６よりの最大利得ビームフォーマｗｎ（ｆ）と、周波数領域変換部５よりの観測信号ベクトルｘ（ｆ、τ）とは、信号抽出部８１に入力される。信号抽出部８１では、以下の式（４７）を計算して、歪みを含んだ目的信号ｙｎ（ｆ、τ）を抽出する。
ｙｎ（ｆ、τ）＝ｗｎＨ(ｆ)ｘ（ｆ、τ）      （４７）
歪みを含んだ目的信号ｙｎ（ｆ、τ）は歪み補正部８２に入力される。
また、補正ベクトル計算部７８よりの補正ベクトルαｎ(ｆ)も歪み補正部８２に入力される。歪み補正部８２では、以下の式（４８）で出力信号を変換することで、歪みを補正して補正出力信号ｙｎ’（ｆ、τ）を出力する。
ｙｎ’(ｆ、τ)＝［αｎ(ｆ)］Ｂｙｎ(ｆ、τ)    （４８）
なお、以上で説明した実施例１〜５では、全てのｎについて信号を抽出するとしてきたが、単独の信号（１つのｎ）についてのみ、ビームフォーマを構成するだけでもよい。目的信号の選択については、例えば、データベース上の目的信号のインパルス応答ベクトルｈｄと発明法により全ての音源ｎについて推定されたインパルス応答ベクトルｈｎを比較して、最もｈｄに近いｈｎを持つ音源ｎを選ぶことで選択できる。例えば、ｍｉｎｎ（ｈ１・ｈｎ）などのアルゴリズムが考えられる。選ばれたｎについてのみ実施例２〜５で説明したビームフォーマ計算部２８による上記式（２４）などを用いたビームフォーマを構成すれば、目的信号についての適応型ビームフォーマを得ることができる。
【００６１】
［実験結果］
上記実施例の効果を示すために、実験を行った。図１１に示す部屋で測定したインパルス応答を複数の音声に畳み込み混合することで、混合信号を模した。実験条件は図１１に示す通りである。長辺が８８０ｃｍ、短辺が３７５ｃｍ、高さが２４０ｃｍ、残響は１２０ｍｓの室内において、底面の長辺から２００ｃｍ、短辺から２８２ｃｍの位置に３つのセンサ４１、４２、４３を配置した。長辺と平行軸をｘ、短辺と平行軸をｙとし、図１２に示すように、３つのセンサ４１、４２、４３をｙ軸上に２個、ｘ軸上に１個、辺の長さ４ｃｍの正三角形の頂点につまり２次元に配した場合の実験を行う。またセンサとしてはマイクロホンを用いた。
４通りの音声組み合わせについて、信号対不要信号比（ＳＩＲ）と信号対歪み比（ＳＤＲ）を評価した。なお、単位はｄＢである。
【００６２】
４つの音源をセンサ位置におけるｘ軸とｙ軸の交点を中心とし、ｘ軸の＋方向を０度とし、左回りに３０度、３１５度方向とセンサ位置を中心と半径５０ｃｍの円との各交差点上にそれぞれの音源を、２２５度、３１５度の方向と半径８０ｃｍの円との交差点上に、それぞれ音源を配置させる。実施例２の効果を確かめる実験では、１２０度、２２５度、３１５度方向の音源を用い、Ｎ（源信号の数）＝Ｍ（センサの数）＝３とした。また実施例３の効果を確かめる実験ではＮ＝４、Ｍ＝３とした。
【００６３】
図１３にこの実験の結果を示す。実施例３は従来法、実施例２、実施例２’、実施例４、実施例５において、図７記載の入力信号推定部５０を設けた場合を示す。従来法では、図１記載の適応型ビームフォーマ６を表す上記式（１０）において、ｈｎ（ｆ）に既知のステアリングベクトルａｎ（ｆ）を与えたものを用いた。この場合、Ｎ＝Ｍの場合も、Ｎ＞Ｍの場合も、共に高い性能を得られなかった。これは、残響のある環境での実験であるため、与えられたステアリングベクトルａｎ（ｆ）が残響の影響まで、考慮できなかったことが主な原因として考えられる。また、Ｎ＞Ｍの場合に十分なＳＩＲが得られないのは、適応型ビームフォーマの限界つまりＭ−１個の不要信号しか効果的に抑圧できないことを示している。
【００６４】
従来法に対し、Ｎ＝Ｍの場合、Ｎ＞Ｍの場合であっても、上記実施例はＳＩＲ、ＳＤＲの値を比較すると、従来法よりも高い性能を持つことが理解される。
【００６５】
以上の各実施形態の他、本発明であるブラインド信号抽出装置は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、ブラインド信号抽出装置において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【００６６】
また、この発明のブラインド信号抽出装置における処理をコンピュータによって実現する場合、ブラインド信号抽出装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、ブラインド信号抽出装置における処理機能がコンピュータ上で実現される。
【００６７】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｃ）、ＤＶＤ−ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｃ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＣＤ−Ｒ（Ｒｅｃｏｒｄａｂｌｅ）／ＲＷ（ＲｅＷｒｉｔａｂｌｅ）等を、光磁気記録媒体として、ＭＯ（Ｍａｇｎｅｔｏ−Ｏｐｔｉｃａｌ  ｄｉｓｃ）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Ｅｌｅｃｔｒｏｎｉｃａｌｌｙ  Ｅｒａｓａｂｌｅ  ａｎｄ  Ｐｒｏｇｒａｍｍａｂｌｅ−Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）等を用いることができる。
【００６８】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００６９】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Ａｐｐｌｉｃａｔｉｏｎ  Ｓｅｒｖｉｃｅ  Ｐｒｏｖｉｄｅｒ）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００７０】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ブラインド信号抽出装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【００７１】
この発明は、オーディオ分野の応用として、音声認識機の入力マイクロホンと話者が離れた位置にあるためマイクロホンが目的話者音声以外の音まで収音してしまうような状況でも、目的音声を分離抽出することで、認識率の高い音声認識系の構築が可能になる。
【図面の簡単な説明】
【００７２】
【図１】従来技術のシステムの機能構成例を示すブロック図。
【図２】直線状に配置したセンサシステムを用いた場合において、音源ｎが任意のセンサｊに達する時刻と原点０に達する時刻との時間差τを説明するための図。
【図３】この発明の実施例１のシステムの機能構成例を示すブロック図。
【図４】この発明の実施例１の主な処理の流れを示すフローチャート。
【図５】任意の２つのセンサであるセンサｍとセンサＱとにおいて、上記式（１６）で使用するｃｏｓΘｎｍＱを説明するための図。
【図６】この発明の実施例２のシステムの機能構成例の一部を示すブロック図。
【図７】この発明の実施例３のシステムの機能構成例の一部を示すブロック図。
【図８】この発明の実施例４のシステムの機能構成例の一部を示すブロック図。
【図９】この発明の実施例５のシステムの機能構成例の一部を示すブロック図。
【図１０】この発明の実施例５の変形例のシステムの機能構成例の一部を示すブロック図。
【図１１】従来の技術とこの発明の技術との比較実験を真上から見た図。
【図１２】図１１の３つのセンサ４１、４２、４３の位置関係の詳細を示す図。
【図１３】従来の技術とこの発明の技術の効果を比較した実験結果を示す図。
【図１】
2008060635.tif 000018
【図２】
2008060635.tif 000019
【図３】
2008060635.tif 000020
【図４】
2008060635.tif 000021
【図５】
2008060635.tif 000022
【図６】
2008060635.tif 000023
【図７】
2008060635.tif 000024
【図８】
2008060635.tif 000025
【図９】
2008060635.tif 000026
【図１０】
2008060635.tif 000027
【図１１】
2008060635.tif 000028
【図１２】
2008060635.tif 000029
【図１３】
2008060635.tif 000030