(57)【要約】      （修正有）【課題】時間周波数領域において、最小２乗型同時対角化問題の解法と電力比の相関と多数決を利用したパーミュテーション問題の解法を用いて観測信号から未知信号源信号を推定するブラインド信号分離方法とその装置を提供する。【解決手段】窓関数を乗算した観測信号の短時間離散フーリエ変換後、観測信号の空間相関行列をエポック時刻毎に求め、制約条件付きフォワードモデル型最小２乗型同時対角化問題の対象にする。ラグランジェの未定乗数法と反復法を組み合わせた手順の収束後、累乗法を１回用いて混合行列を求め、混合行列のランクが落ちていた場合には特異値分解を用いて基底と特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求め、同一信号源から発生した信号の周波数ビン間の電力比相関性と多数決を組み合せてパーミュテーション問題を解決する。【選択図】図１
2015210512.tif 000002
【特許請求の範囲】
【請求項１】
未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみからブラインドで推定する方法であって、時間周波数領域においてエポック時刻毎に観測信号ベクトルから空間相関行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる空間相関行列を求め、そのノルムで全ての同エポック時刻の空間相関行列を正規化した後、正規化された空間相関行列から混合行列を推定するために、制約条件付きフォワードモデル型最小２乗型同時対角化問題をラグランジュの未定乗数法と反復法の組み合わせにより解法する手順と、手順の収束後、累乗法を１回用いて最小２乗型同時対角化問題の解となる混合行列を求め、混合行列のランクが落ちていた場合は、混合行列がフルランクになるように特異値分解を用いて基底を補い、次いで、ランク落ちした混合行列の特異値とその総和が変わらないように特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求め、同一信号源から発生した信号の周波数ビンの電力比に相関があることを利用してパーミュテーション問題を、多数決を利用して解法することによって信号分離精度の高い分離行列を求めることを特徴とするブラインド信号分離方法。
【請求項２】
請求項１記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、信号源とブラインド信号分離装置出力間の伝達関数を表した混合行列と分離行列の縦続接続モデルにラグランジュの未定乗数を導入した制約条件付きフォワードモデル型最小２乗型同時対角化問題とバックワード型最小２乗型同時対角化問題を導入し、更に分離行列に遅延を与えることによって信号源とブラインド信号分離装置出力間の因果的なモデルを推定することを特徴とするブラインド信号分離方法。
【請求項３】
請求項１乃至請求項２のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、基準周波数ビンを複数選択し、基準周波数ビン間において最も電力比の相関が大きいパーミュテーション行列を推定し、複数の基準周波数ビンから基準周波数ビンを１つ選択し、全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を推定する手順を、全ての基準周波数ビンが１度選択されるまで繰り返した後、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定し、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用することを特徴とするブラインド信号分離方法。
【請求項４】
請求項１乃至請求項３のいずれか１項に記載のブラインド信号分離方法を用いて信号源分離を行うように構成されていることを特徴とするブラインド信号分離方法を用いたブラインド信号分離装置。
【発明の詳細な説明】
【技術分野】
【０００１】  本発明は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみから推定するブラインド信号分離方法に係わり、特に、最小２乗型同時対角化問題の解法を用いて高い精度で信号を分離することができるブラインド信号分離装置に関する。
【背景技術】
【０００２】  複数の未知信号源信号が未知の畳み込み混合系により混在されて観測されるとき、観測信号を分離して混在前の未知信号源信号を推定する処理をブラインド信号分離という。ブラインド信号分離では、未知信号源信号間の統計的独立性のみを条件として、観測信号から未知信号源信号を推定する方法であり、信号源の位置或いは観測信号の到来方向の推定を必ずしも必要としない方法である。【０００３】  勾配法に基づく周波数領域の信号分離方法が非特許文献１で提案されている。勾配法に基づき周波数ビン毎に分離行列を更新することによって分離信号を推定する。【０００４】  最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法が非特許文献２と非特許文献３で提案されている。これらの方法は、最小２乗型同時対角化問題の解法による混合行列の推定、最小２乗型一般化逆行列を用いた混合行列からの分離行列の推定、スケーリング問題の解法、パーミュテーション問題の解法の４つの手順から成る。【０００５】  周波数ビン間の電力比の相関に基づきパーミュテーション行列を推定する方法が非特許文献４で提案されている。
【先行技術文献】
【０００６】
【非特許文献１】  Ｌ．Ｐａｒｒａ  ａｎｄ  Ｃ．Ｓｐｅｎｃｅ，”Ｃｏｎｖｏｌｕｔｉｖｅ  ｂｌｉｎｄ  ｓｅｐａｒａｔｉｏｎ  ｏｆ  ｎｏｎ−ｓｔａｔｉｏｎａｒｙ  Ｓｏｕｒｃｅｓ”，ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｓｐｅｅｃｈ  ａｎｄ  Ａｕｄｉｏ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．８，ｎｏ．３，ｐｐ．３２０−３２７，Ｍａｙ，２０００
【非特許文献２】  Ａ．Ｙｅｒｅｄｏｒ，”Ｎｏｎ−ｏｒｔｈｏｇｏｎａｌ  ｊｏｉｎｔ  ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  ｉｎ  ｔｈｅ  ｌｅａｓｔ−ｓｑｕａｒｅｓ  ｓｅｎｓｅ  ｗｉｔｈ  ａｐｐｌｉｃａｔｉｏｎ  ｉｎ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ”，ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５０，ｎｏ．７，ｐｐ．１５４５−１５５３，Ｊｕｌｙ  ２００２
【非特許文献３】  Ｋ．Ｒａｈｂａｒ  ａｎｄ  Ｊ．Ｐ．Ｒｅｉｌｌｙ，”Ａ  ｆｒｅｑｕｅｎｃｙ  ｄｏｍａｉｎ  ｍｅｔｈｏｄ  ｆｏｒ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ  ｏｆ  ｃｏｎｖｏｌｕｔｉｖｅ  ａｕｄｉｏ  ｍｉｘｔｕｒｅｓ”，ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｓｐｅｅｃｈ  ａｎｄ  Ａｕｄｉｏ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．１３，ｎｏ．５，ｐｐ．８３２−８４４，Ｓｅｐｔ．２００５
【非特許文献４】  Ｈ．Ｓａｗａｄａ，Ｓ．Ａｒａｋｉ，ａｎｄ  Ｓ．Ｍａｋｉｎｏ，”Ｍｅａｓｕｒｉｎｇ  ｄｅｐｅｎｄｅｎｃｅ  ｏｆ  ｂｉｎ−ｗｉｓｅ  ｓｅｐａｒａｔｅｄ  ｓｉｇｎａｌｓ  ｆｏｒ  ｐｅｒｍｕｔａｔｉｏｎ  ａｌｉｇｎｍｅｎｔ  ｉｎ  ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ  ＢＳＳ”，Ｐｒｏｃ．ＩＥＥＥ  Ｉｎｔ．Ｓｙｍｐ．Ｃｉｒｃｕｉｔｓ  Ｓｙｓｔ．，ｐｐ．３２４７−３２５０，Ｍａｙ  ２００７
【発明の概要】
【発明が解決しようとする課題】
【０００７】  従来のブラインド信号分離方法で使用されている勾配法は演算量が少ないが、収束が遅く、実環境下で十分な信号分離精度を得ることができない。【０００８】  従来のブラインド信号分離方法では、時間周波数領域で周波数ビンとエポック毎、即ち、ポイント毎に観測信号の空間相関行列をそのフロベニウスノルムで正規化しているため、未知信号源信号が音声の場合、音声の特徴を表す周波数特性が失われ、周波数ビンに対して振幅特性が一定になり、信号分離精度が劣化する。【０００９】  従来のブラインド信号分離方法では、時間周波数領域で観測信号の空間相関行列から最小２乗型同時対角化問題の解法を用いて混合行列を推定し、次いで最小２乗型一般化逆行列を用いて混合行列から分離行列を推定しているため、最小２乗型同時対角化問題に分離行列の推定が含まれず、仮に混合行列を精度良く推定できても分離行列が精度良く推定できるとは限らず、信号分離精度が劣化する場合がある。【００１０】  従来のブラインド信号分離方法では、時間周波数領域で観測信号の空間相関行列から制約条件を課すことなく最小２乗型同時対角化問題の解法を用いて混合行列を推定し、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定する前に、塁乗法により制約条件を課しているので、手順が収束に要する反復回数は大幅に増加する。【００１１】  更に、従来のブラインド信号分離方法において混合行列の最小２乗型一般化逆行列を単に分離行列として取り扱うと、加法性雑音が存在する場合、雑音の影響が大きくなって信号分離性能の劣化につながる。【００１２】  同一信号源から発生した信号の隣接周波数ビンに相関があることを利用した従来のパーミュテーション問題の解法では、１度誤りが生じると、これ以降の解法が誤り続ける確率が非常に高くなる。【００１３】  本発明はこのような事情を鑑みてなされたものであり、信号源とブラインド信号分離装置出力間を表した混合行列と分離行列の縦続接続モデルを最小２乗型同時対角化問題に取り込み、更にモデルに因果性を与えることによって信号分離性能の向上と反復手順の収束高速化の両立、更に、同一信号源から発生した信号の周波数ビン間の電力比相関性と多数決を組み合わせによってパーミュテーション問題を正確に解法することを目的とする。
【課題を解決するための手段】
【００１４】  このような目的に応えるために本発明（請求項１記載の発明）に係るブラインド信号分離方法は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみからブラインドで推定する方法であって、時間周波数領域においてエポック時刻毎に観測信号ベクトルから空間相関行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる空間相関行列を求め、そのノルムで全ての同エポック時刻の空間相関行列を正規化した後、正規化された空間相関行列から混合行列を推定するために、制約条件付きフォワードモデル型最小２乗型同時対角化問題をラグランジュの未定乗数法と反復法の組み合わせにより解法する手順と、手順の収束後、累乗法を１回用いて最小２乗型同時対角化問題の解となる混合行列を求め、混合行列のランクが落ちていた場合は、混合行列がフルランクになるように特異値分解を用いて基底を補い、次いで、ランク落ちした混合行列の特異値とその総和が変わらないように特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求め、同一信号源から発生した信号の周波数ビンに相関があることを利用してパーミュテーション問題を、多数決を利用して解法することによって信号分離精度の高い分離行列を求めることを特徴とする。【００１５】  本発明（請求項２記載の発明）に係るブラインド信号分離方法は、信号源とブラインド信号分離装置出力間の伝達関数を表した混合行列と分離行列の縦続接続モデルにラグランジュの未定乗数を導入した制約条件付きフォワードモデル型最小２乗型同時対角化問題とバックワード型最小２乗型同時対角化問題を導入し、更に分離行列に遅延を与えることによって信号源とブラインド信号分離装置出力間の因果的なモデルを推定することを特徴とする。【００１６】  本発明（請求項３記載の発明）に係るブラインド信号分離方法は、基準周波数ビンを複数選択し、基準周波数ビン間において最も電力比の相関が大きいパーミュテーション行列を推定し、複数の基準周波数ビンから基準周波数ビンを１つ選択し、全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を推定する手順を、全ての基準周波数ビンが１度選択されるまで繰り返した後、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定し、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用することを特徴とする。【００１７】  すなわち、本発明によれば、観測信号の空間相関行列を制約条件付き最小２乗型同時対角化問題の対象に、ラグランジュの未定乗数法と反復法を組み合わせることによって最小２乗型同時対角化問題の近似解、即ち、混合行列を求めることができ、推定性能の向上と反復回数の低減を両立させる。【００１８】  また、本発明によれば、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって当該周波数ビンに割り当てられるパーミュテーション行列を決定し、全体の信号分離性能を向上させる。
【発明の効果】
【００１９】  本発明によれば、時間周波数領域においてエポック時刻毎に観測信号ベクトルから空間相関行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる空間相関行列を求め、そのノルムで全ての同エポック時刻の空間相関行列を正規化した後、これに制約条件付き最小２乗型同時対角化問題を適用する。制約条件付き最小２乗型同時対角化問題では、ラグランジュの未定乗数法と反復法を組み合わせることによって近似解、即ち、混合行列を求める。本発明に係るブラインド信号分離方法及びブラインド信号分離装置では、混合行列のランクが落ちていた場合は、混合行列がフルランクになるように特異値分解を用いて基底と特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求めるため、少ない反復回数で近似解に収束するという効果がある。【００２０】  また、本発明に係るブラインド信号分離方法及びブラインド信号分離装置では、信号源とブラインド信号分離装置出力間の伝達関数を表した混合行列と分離行列の縦続接続モデルにラグランジュの未定乗数を導入した制約条件付きフォワードモデル型最小２乗型同時対角化問題とバックワード型最小２乗型同時対角化問題を導入し、更に分離行列に遅延を与えることによって信号源とブラインド信号分離装置出力間の因果的なモデルを推定することによって雑音の影響を最小にする分離行列を推定できるという効果がある。【００２１】  更に、パーミュテーション問題の解法では、基準周波数ビンを複数選択し、基準周波数ビン間において最も電力比の相関が大きいパーミュテーション行列を推定し、複数の基準周波数ビンから基準周波数ビンを１つ選択し、全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を推定する手順を全ての基準周波数ビンが１度選択されるまで繰り返した後、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定し、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用することを特徴とすることによって本発明に係るブラインド信号分離方法及びブラインド信号分離装置には、信号分離性能を高めることができるという効果がある。
【図面の簡単な説明】
【００２２】
【図１】  本発明に係るブラインド信号分離方法の実施の形態を示す図である。
【図２】  本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角化行列のラグランジェ未定乗数法の適用について説明するためのフローチャートである。
【図３】  本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角行列の解法について説明するためのフローチャートである。
【図４】  本発明に係るブラインド信号分離方法のパーミュテーション問題の解法について説明するためのフローチャートである。
【図５】  本発明に係る信号分離方法の実施例１における信号源（スピーカ）とマイクロホンの位置関係を表す平面図である。
【図６】  本発明に係る信号分離方法の実施例１における部屋の残響時間と信号分離性能の関係を示す図である。
【図７】  本発明に係る信号分離方法の実施例１におけるＳＮＲと信号分離性能の関係を示す図である。
【図８】  実施例２における従来のパーミュテーション問題の解法によるパーミュテーション行列の割り当て結果を示す図である。
【図９】  実施例２における本発明に係るブラインド信号分離方法によるパーミュテーション行列の割り当て結果を示す図である。
【図１０】  分離行列の推定に本発明に係るブラインド信号分離方法を使用して、パーミュテーション問題の解法には非特許文献４を使用した場合と、本発明に係るブラインド信号分離方法の実施例２におけるパーミュテーション行列の割り当て正答率と信号分離性能を比較する図である。
【発明を実施するための形態】
【００２３】  本発明に係るブラインド信号分離方法の実施の形態について図面を参照して説明する。
【００２４】１．畳み込み混合モデル
図１に示すように、時刻ｔにおいてＮ個の信号源１１、１２、…、１Ｎから発せられた信号源信号ｓｊ（ｔ）が畳み込み混合されてｘｉ（ｔ）として観測される。ｓｊ（ｔ）は平均０で互いに統計的独立な非定常信号である。ｈｉｊ（ｔ）は信号源１ｊからマイクロホン２ｉまでの経路の時不変なインパルス応答で、因巣的で非最小位相系である。また、ｎｉ（ｔ）はマイクロホン２ｉに加わる平均０、分散σ２のガウス性白色雑音で、ｓｊ（ｔ）と統計的独立である。時刻ｔにおいてＪ個のマイクロホン２１、２２、…、２Ｊで観測される観測信号ｘｉ（ｔ）は式（１）で表される。ここで、Ｊ≧Ｎ≧２とする。
【００２５】
【数１】
2015210512.tif 000003
ここで、＊は畳み込み演算を表す。
【００２６】  観測信号ｘｉ（ｔ）を３１で短時間フーリエ変換すると、フレーム時刻ｍにおける観測信号は式（２）により表される。式（２）において、ｗｉｎ（ｔ）は窓関数、Ｋは短時間フーリエ変換の点数、Ｔｂは２つの重複窓間のシフトサイズ、ωｋ＝２πｋ／Ｋ、ｋ＝０，１，…，Ｋ−１、をそれぞれ表す。離散フーリエ変換の点数Ｋがｈｉｊ（ｔ）のインパルス応答長より十分に大きいとき観測信号は式（３）により近似される。ここで、ｈｉｊ（ｔ）のＫ点離散フーリエ変換をｈｉｊ（ωｋ）、ｓｊ（ｔ）に窓関数を乗算した後、Ｋ点短時間フーリエ変換で時間周波数領域に変換したフレーム時刻ｍの信号源信号をｓｊ（ωｋ，ｍ）、同様に、ｎｉ（ｔ）に窓関数を乗算した後、Ｋ点短時間フーリエ変換で時間周波数領域に変換したフレーム時刻ｍの雑音をｎｉ（ωｋ，ｍ）とそれぞれ表記している。また、式（３）において、ｘ（ωｋ，ｍ）はフレーム時刻ｍに各マイクロホンでの観測信号ベクトル、ｓ（ωｋ，ｍ）はフレーム時刻ｍに各信号源信号ベクトル、混合行列Ｈ（ωｋ）はＮ個の信号源からＪ個のマイクロホンまでの混合行列、ｎ（ωｋ，ｍ）は雑音ベクトルでそれぞれ式（４）、（７）、（５）、（８）により定義される。入手可能なエポック時刻の総数をＭとすると、１≦ｍ≦Ｍとなる。信号源信号の共分散行列はＰｓ（ωｋ，ｍ）＝Ｅ［ｓ（ωｋ，ｍ）ｓ（ωｋ，ｍ）Ｈ］∈ＲＮ×Ｎで、対角行列となる。Ｅ［・］と上付き添字Ｈは期待値と複素共役転置をそれぞれ表す。また、上付き添字ＴとＲＮ×Ｎは転置とＮ×Ｎの実数空間を表す。
【００２７】
【数２】
2015210512.tif 000004
【００２８】  信号を分離するには５１、５２、…、５Ｋで周波数ビン毎に式（９）を満足する分離行列Ｗ（ωｋ）を推定し、６０で信号源の割り当てを定めるパーミュテーション行列Π（ωｋ）∈ＲＮ×Ｎを決定する。周波数ビン毎に独立にΠ（ωｋ）を決定しても信号が完全に分離する保証はなく、同一信号源から発生した信号の隣接または近接周波数ビンに相関があることを利用してパーミュテーション行列Π（ωｋ）を決定する。
【００２９】
【数３】
2015210512.tif 000005
ここで、Ｄ（ωｋ）∈ＣＮ×Ｎは周波数ビン毎に異なる任意の対角行列である。
【００３０】  スケーリング問題とパーミュテーション問題を順に解法した後、７１、７２、…、７Ｋでｘ（ωｋ，ｍ）に左から分離行列Ｗ（ωｋ）を乗算すると、周波数ビンωｋにおける分離信号ｙ（ωｋ，ｍ）は式（１０）で表される。尚、スケーリング問題の解法については後述する。式（１０）を８０で短時間逆フーリエ変換と重複加算によって時間領域に変換すると分離信号ｙｉ（ｔ）が求められる。雑音の分散σ２が十分に小さいとき、ｙｉ（ｔ）≒ｓｉ（ｔ）になる。尚、分離信号ベクトルｙ（ωｋ，ｍ）は式（１１）により表される。
【００３１】
【数４】
2015210512.tif 000006
【００３２】  本発明に関するブラインド信号分離方法について、図１乃至図４を参照して詳細に説明する。図２乃至図４は、図１の３１における短時間フーリエ変換後、ブラインド信号分離システム４０において本発明により周波数ビン毎に推定される分離行列の算出手順を示したものである。図４は、分離行列の算出後、ブラインド信号分離システム４０において本発明によりパーミュテーション行列の算出手順を示したものである。
【００３３】  観測信号ｘ（ωｋ，ｍ）の共分散行列Ｐｘ（ωｋ，ｍ）∈ＣＪ×Ｊは式（１２）で与えられる。式（１３）の制約条件を課して式（１４）を満足する対角化行列Ｂ（ωｋ）と対角行列Λ（ωｋ，ｍ）を求めると、式（９）よりＢ（ωｋ）とＷ（ωｋ）の関係は式（１５）で与えられる。ただし、ＣＪ×ＪはＪ×Ｊの複素空間を表す。
【００３４】
【数５】
2015210512.tif 000007
ただし、Ｉは単位行列である。
【００３５】２．最小２乗型同時対角化問題とその解法
観測信号ｘ（ωｋ，ｍ）の共分散行列Ｐｘ（ωｋ，ｍ）の推定値Ｐｘ（ωｋ，ｍ）を正規化して、式（１７）を最小にする対角化行列Ｂ（ωｋ）と対角行列Λ（ωｋ，ｍ）を求める。式（１７）は最小２乗型同時対角化問題の解法として知られている。
【数６】
2015210512.tif 000008
【００３６】  本発明では、制約条件付き最小２乗型同時対角化問題を解法することによって混合行列Ｂ（ωｋ）を推定した後、分離行列Ｗ（ωｋ）を求める。次いで、対角行列Λ（ωｋ，ｍ）を推定するために、最小２乗法を用いて分離行列Ｗ（ωｋ）を用いた評価量を最小化する。本発明では、混合行列と対角行列の推定を交互に繰り返す。音声は低域周波数帯にフォルマントと呼ばれる振幅スペクトルのピークを有している。この音声波形の特徴を失うことなく、式（１９）によって観測信号ｘ（ωｋ，ｍ）の共分散行列Ｐｘ（ωｋ，ｍ）を正規化することが、本発明の特徴の一つである。
【００３７】２．１  時間周波数領域における観測信号の共分散行列の正規化
時間領域観測信号は式（１８）の短時間フーリエ変換によって時間周波数領域に変換される。
【００３８】
【数７】
2015210512.tif 000009
式（１７）において、ｗ（ｔ）は窓関数、Ｋは短時間フーリエ変換の点数、Ｔｓは２つの重複窓間のシフトサイズ、Ｔｓはエポックサイズ、Ｎｓ＋１は各エポックにおける総重複フレーム数でＫ＋ＮｓＴｓ≦Ｔｂ、ｌ＝０，１，…，Ｎｓをそれぞれ表す。
【００３９】  ステップＳ１０１においてエポック時刻ｍにおける時間周波数領域観測信号の共分散行列を式（１９）によって推定される。
【００４０】
【数８】
2015210512.tif 000010
【００４１】  混合行列Ｂ（ωｋ）と分離行列Ｗ（ωｋ）を縦続に接続したとき、そのインパルス応答は式（１５）を最小にすることによって求められる。式（２１）の評価量ｅ（ωｋ）をＷ（ωｋ）によって微分すると、分離行列Ｗ（ωｋ）は式（２２）によって求められる。Ｂ（ωｋ）ＨＢ（ωｋ）のランクがＮのときのみ、式（２１）の評価量ｅ（ωｋ）は零になる。一方、Ｂ（ωｋ）ＨＢ（ωｋ）のランクがＮ未満のとき、ｅ（ωｋ）は零より大きくなる。そこで、制約条件‖ｂｊ（ωｋ）‖２＝１に制約条件ｒａｎｋ（Ｂ（ωｋ）ＨＢ（ωｋ））＝Ｎを付け加え、混合行列Ｂ（ωｋ）の推定のための最小２乗型同時対角化問題を解法する。ここで、ｂｊ（ωｋ）はＢ（ωｋ）のｊ番目の列ベクトル、‖・‖２はユークリッドノルム、ｒａｎｋ（Ａ）は行列Ａのランクをそれぞれ表す。
【００４２】
【数９】
2015210512.tif 000011
ここで、‖・‖Ｆはフロベニウスノルムを表す。
【００４３】  ステップＳ１０２においてＰｘ（ωｋ，ｍ）を式（２３）によって正規化する。
【００４４】
【数１０】
2015210512.tif 000012
【００４５】２．２  対角化行列の解法
制約条件‖ｂｊ（ωｋ）‖２＝１を課したフォワードモデル型最小２乗型同時対角化問題を周波数ビンωｋ毎に解くことによって、対角化行列Ｂ（ωｋ）、即ち、混合行列を求める。
評価量を式（２４）に示す。式（２４）はフォワードモデル型最小２乗型同時対角化問題として知られている。
【００４６】
【数１１】
2015210512.tif 000013
ここで、γｉはラグランジェの未定乗数を表す。
【００４７】  ベクトル表現を用いると、式（２４）の評価関数は式（２６）のように表現することができる。ここで、ｒｘ（ωｋ，ｍ）、Ｇ（ωｋ）、ｄ（ωｋ，ｍ）、Ｇ（ωｋ）ｄ（ωｋ，ｍ）はそれぞれ式（２８）〜（３１）により表される。ただし、ｖｅｃ｛Ａ｝は行列Ａの列を積み重ね
2015210512.tif 000014
れぞれ表す。λｉは対角行列Ａのｉ番目の要素を表す。
【００４８】
【数１２】
2015210512.tif 000015
【００４９】  反復法によって混合行列Ｂ（ωｋ）を求めるためにｚｉ（ωｋ）とＴ（ωｋ）を式（３２）、（３３）によって定義され、ステップＳ１０３において作成される。
【００５０】
【数１３】
2015210512.tif 000016
【００５１】  ステップＳ１０４においてｇｊ（ωｋ）を求める際、ｇｊ（ωｋ）（ｊ≠ｉ）を式（３４）のように、ｚｊ（ωｋ）を式（３４）のようにそれぞれ定数に設定して、Ｆｊ（ωｋ）を計算する。ここで、ｇｊ（ωｋ）はＧ（ωｋ）のｊ番目の列ベクトルである。
【００５２】
【数１４】
2015210512.tif 000017
【００５３】  式（２６）の制約条件付きフォワードモデル型最小２乗型同時対角化問題を式（３５）のように書き直すことができ、その近似解ｇｊ（ωｋ）はステップＳ１０５において式（３６）のラグランジェの未定乗数法によって求められる。
【００５４】
【数１５】
2015210512.tif 000018
ここで、ｕｎｖｅｃ｛Ａ｝は、Ｊ２×１の列ベクトルＡをＪ×Ｊの行列に変換することを表す。
【００５５】  ステップＳ１０７においてｊ＝１，２，…，Ｎについて誤差の限界がεＧの近似解ｇｊ（ωｋ）を反復法によって推定した後、ステップＳ１０８において累乗法を１回用いて式（３７）を最小にするｂｊ（ωｋ）を算出する。次いで、ステップＳ１０９において式（３８）のようにＢ（ωｋ）を特異値分解する。ここで、ｔｒ［Ａ］は行列Ａのトレースを表す。
【００５６】
【数１６】
2015210512.tif 000019
【００５７】  ステップＳ１１０においてＢ（ωｋ）ＨＢ（ωｋ）のランクがＮ未満のとき、式（２１）の評価量ｅ（ωｋ）を零にするために、Ｂ（ωｋ）を式（４４）の行列によって置き換える。ここで、ステップＳ１１１において正規直交基底ｖ１（ωｋ），ｖ２（ωｋ），…，ｖｒ（ωｋ）によって張られる空間に直交する空間の正規直交基底ｖｒ＋１（ωｋ），ｖｒ＋２（ωｋ），…，ｖＮ（ωｋ）、同様に、ステップＳ１１２において正規直交基底ｕ１（ωｋ），ｕ２（ωｋ），…，ｕｒ（ωｋ）によって張られる空間に直交する空間の正規直交基底ｕｒ＋１（ωｋ），ｕｒ＋２（ωｋ），…，ｕＮ（ωｋ）はそれぞれ求められる。
【００５８】
【数１７】
2015210512.tif 000020
ここで、δ（ωｋ）＞０とする。ステップＳ１１３において追加される特異値δ（ωｋ）は、式（４４）の右辺の√Ｎ／（√Ｎ＋δ（ωｋ）Ｎ）によって条件ｔｒ［Σ（ωｋ）］＝√Ｎを満足するように設定される。
【００５９】    ２．３  対角行列の解法
ステップＳ１１４においてＢ（ωｋ）から分離行列Ｗ（ωｋ）を式（２２）の最小２乗型一般化逆行列によって求める。式（１４）、（１５）よりΛ（ωｋ，ｍ）の左からＷ（ωｋ）Ｂ（ωｋ）、右からＢ（ωｋ）ＨＷ（ωｋ）Ｈをそれぞれ乗算すると、式（４７）を得る。
【数１８】
2015210512.tif 000021
【００６０】  Ｐｘ（ωｋ，ｍ）の推定値を使用して誤差Ψ（ωｋ，ｍ）を式（４８）で定義すると、式（４９）のバックワード型最小２乗型同時対角化問題に最小２乗法を適用すると、対角行列Λ（ωｋ，ｍ）はステップＳ１１５において式（５０）で推定される。式（５０）においてｄｉａｇ［Ａ］は行列Ａの対角行列を表す。
【００６１】
【数１９】
2015210512.tif 000022
【００６２】  誤差の限界がεＣの近似解ｇｊ（ωｋ）と近似解Λ（ωｋ，ｍ）を推定するまで、上記のアルゴリズムはステップＳ１１６において繰り返される。
【００６３】３．パーミュテーション問題の解法
ステップＳ１１７において基準周波数ビンを複数選択し、基準周波数ビン間において電力比の相関に基づきパーミュテーション行列を推定するためにΞ（ωｋ）を式（５３）によりステップＳ１１８で算出する。
【００６４】
【数２０】
2015210512.tif 000023
ここで、Ｔｒ（・）は行列のトレースを表す。また、Ｑは行列の各行に１となる要素が１箇所、その他の要素は０で、１となる要素の位置が他の行と重複しない行列の集合である。
【００６５】  ステップＳ１１９において基準周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を式（５４）によって推定する。
【００６６】  ステップＳ１２０において複数の基準周波数ビンから１つの基準周波数ビンを任意に選択し、ステップＳ１２１において全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を式（５４）によって推定する。
【００６７】  ステップＳ１２２において選択された全ての基準周波数ビンが１度選択されるまで、上記のパーミュテーション行列の推定手順を繰り返す。この結果、基準周波数ビンを除き、各周波数ビンに複数のパーミュテーション行列が割り当てられることになる。
【００６８】  ステップＳ１２３において多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定する。ただし、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用する。式（５５）のように観測信号ｘ（ωｋ，ｍ）にΠ（ωｋ）Ｗ（ωｋ）を左から乗算して分離信号ｙ（ωｋ，ｍ）を得る。
【実施例１】
【００６９】４．１  評価データ
図５のように４．４５×３．５５×２．５メートルの部屋に３個の信号源（スピーカ）１１、１２、１３を半径１．２メートルの円の円周上に、円の中心に位置する一辺が２０センチメートルの正三方形の頂点に３個のマイクロホン２１、２２、２３をそれぞれ配置した。尚、図５は信号源（スピーカ）とマイクロホンの位置関係を示す平面図である。部屋の残響時間は１００ミリ秒から９００ミリ秒に設定し、標本化周波数８ｋＨｚ、量子化ビット数１６ビットで信号源とマイクロホンの間のインパルス応答は人工的に発生させた。実験条件は、１０００秒の音声データ、Ｋ＝８１９２点の短時間フーリエ変換、エポック当たり重複率９９％の２３個のフレームの使用、窓関数にはハニング窓を用いた。ＳＮＲは５ｄＢ間隔で０〜３０ｄＢの範囲で変化させた。マイクロホン２１、２２、２３のＳＮＲの設定方法については４．２で説明する。本発明に係るブラインド信号分離方法では、εＧ＝εＣ＝１０−６、δ（ωｋ）＝σｒ（ωｋ）を用いている。スケーリング問題は周波数ビン毎に分離行列の行ベクトルを正規化することによって解法した。Ｃ言語で作成したプログラムをインテル製コアｉ７−２６００  ３．４ＧＨｚプロセッサを用いて実行した。信号源信号からマイクロホンまでの経路は時不変のインパルス応答で、因果的で非最小位相系であるので、因果的な分離行列を実現するために、Π（ωｋ）−１Ｄ（ωｋ）−１Ｗ（ωｋ）にｅ−ｊπｋを乗算した後、逆離散フーリエ変換をして分離フィルタのインパルス応答を得た。【００７０】４．２  評価指標
ブラインド信号分離方法の信号分離性能を次の方法で評価した．式（５６）によって観測信号における所望信号源信号と干渉信号の電力の比、式（５７）によって出力信号における所望信号源信号と干渉信号の電力の比をそれぞれ計算し、ブラインド信号分離装置の各出力の信号分離性能を求める。各出力の平均を信号分離性能とした。γｉｊ（ｔ）は式（５８）のΓ（ωｋ）のｉ行ｊ列の要素を、ｗｉｊ（ｔ）はＷ（ωｋ）の要素をそれぞれ離散逆フーリエ変換したものである。また、分離行列の推定アルゴリズムにおいて収束に要した反復回数と計算時間も評価指標とする。ＳＮＲは、最適な分離行列ｅ−ｊπｋＤ（ωｋ）−１（Ｈ（ωｋ）ＨＨ（ωｋ））−１Ｈ（ωｋ）Ｈとパーミュテーション行列Πｏｐｔ（ωｋ）を使用して観測信号から信号源信号を分離した後、分離信号ｙｉ（ｔ）に含まれる雑音と干渉信号の電力と所望信号源信号の電力の比によって計算した。最適なパーミュテーション行列Πｏｐｔ（ωｋ）は式（５９）によって求めた。また、非ブラインド法は、受信信号を使用して分離行列を計算した後、混合行列が入手可能であるとして、式（６０）によってパーミュテーション行列を求めた。即ち、推定した分離行列に最適なパーミュテーション行列を求めることになり、ブラインド信号分離装置の性能の上限を与えることになる。【００７１】
【数２２】
2015210512.tif 000024
ここで、Ｃｏｐｔ（ωｋ）＝ｅ−ｊπｋＤ（ωｋ）−１（Ｈ（ωｋ）ＨＨ（ωｋ））−１Ｈ（ωｋ）ＨＨ（ωｋ）、Ｃ（ωｋ）＝Ｗ（ωｋ）Ｈ（ωｋ）である。【００７２】４．３  評価対象
勾配法を用いたバックワードモデル型ブラインド信号分離方法（非特許文献１）、最小２乗型同時対角化問題の解法を用いた２種類のフォワードモデル型ブラインド信号分離方法（非特許文献２、非特許文献３）を比較対象とする。従来のブラインド信号分離方法（非特許文献１、非特許文献２、非特許文献３）と本発明に係るブラインド信号分離方法における分離行列の推定精度を比較するため、パーミュテーション行列の推定法は共通の手法（非特許文献４）を使用した。尚、基準周波数ビンの番号には６１４を用いた。【００７６】４．４  評価結果
部屋の残響時間と信号分離性能の関係を図６に、ＳＮＲと信号分離性能の関係を図７にそれぞれ示す。尚、図６においてＳＮＲは２０ｄＢに設定している。両図において太字の数字が最も優れた性能を表している。両図の信号分離性能から明らかなように、計算時間では非特許文献１より劣るものの、本発明に係るブラインド信号分離方法が従来のブラインド信号分離方法よりも最も高い信号分離性能（高い出力ＳＩＲ）を最も少ない反復回数で得ることができた。この要因はラグランジェの未定乗数法を最小２乗型同時対角化問題に導入したこと、推定した混合行列がランク落ちしていた場合、補空間を補い分離行列を推定したことが高い信号分離性能の実現に貢献したと考えられる。また、非ブラインド法の出力ＳＩＲ、即ち、ブラインド信号分離装置の上限に近い値を、本発明に係るブラインド信号分離装置が実現できることが分かる。
【実施例２】
【００７７】５．１  評価データ
図５のように４．４５×３．５５×２．５メートルの部屋に３個の信号源（スピーカ）１１、１２、１３を半径１．２メートルの円の円周上に、円の中心に位置する一辺が２０センチメートルの正三方形の頂点に３個のマイクロホン２１、２２、２３をそれぞれ配置した。部屋の残響時間は７００ミリ秒に設定し、標本化周波数８ｋＨｚ、量子化ビット数１６ビットで信号源とマイクロホンの間のインパルス応答は人工的に発生させた。実験条件は、１０００秒の音声データ、Ｋ＝８１９２点の短時間フーリエ変換、エポック当たり重複率８０％の２個のフレームの使用、窓関数にはハニング窓を用いた。ＳＮＲは２０ｄＢに設定した。本発明に係るブラインド信号分離方法では、εＧ＝εＣ＝１０−６、δ（ωｋ）＝σｒ（ωｋ）、基準周波数ビンの番号は６１６、６１７、６１８を用いている。スケーリング問題は周波数ビン毎に分離行列の行ベクトルを正規化することによって解法した。信号源信号からマイクロホンまでの経路は時不変のインパルス応答で、因果的で非最小位相系であるので、因果的な分離行列を実現するために、Π（ωｋ）−１Ｄ（ωｋ）−１Ｗ（ωｋ）にｅ−ｊπｋ乗算した後、逆離散フーリエ変換をして分離フィルタのインパルス応答を得た。【００７８】５．２  評価指標
信号源とマイクロホンの個数が共に３である場合、式（６１）に示す６種類のパーミュテーション行列の何れか１つが各周波数ビンに割り当てられる。割り当てられたパーミュテーション行列が、任意のパーミュテーション行列に一致する割合と信号分離性能を計算する。【００７９】
【数２３】
2015210512.tif 000025
【００８０】５．３  評価対象
同一信号源から発生した信号の周波数ビン間の電力比に相関があることを利用したパーミュテーション問題の解法（非特許文献４）を比較対象とする。基準周波数ビンの番号には６１４を用いた。従来のブラインド信号分離方法（非特許文献１、非特許文献２、非特許文献３）と本発明に係るブラインド信号分離方法におけるパーミュテーション行列の推定精度を比較するため、分離行列の推定法は共通の手法（本発明に係る分離行列推定法）を使用した。【００８１】５．４  評価結果
図８と図９に各周波数ビンに割り当てられたパーミュテーション行列の番号を×印で示す。パーミュテーション行列Πｉとパーミュテーション行列の番号ｉの関係を式（６１）に示している。図８と図９では、各周波数ビンでパーミュテーション行列の番号３に割り当てられると未知信号源への割り当てが揃うことになる。したがって、番号３を除く他の番号への割り当ては間違いになる。低周波数帯域（０〜２ｋＨｚ）と全周波数帯域におけるパーミュテーション行列の割り当て結果を図１０にまとめる。本発明に係るパーミュテーション行列の推定法が非特許文献４の方法に比べ正答率が向上していることが分かる。また、信号分離性能においても、本発明に係るパーミュテーション行列の推定法が高い出力ＳＩＲを達成することができた。
【符号の説明】
【００８２】１１〜１Ｎ…信号源、２１〜２Ｊ…マイクロホン、３１…短時間フーリエ変換、４０…ブラインド信号分離システム、５１、５２、…、５Ｋ…最小２乗型同時対角化問題の解法、６０…パーミュテーション問題の解法、７１、７２、…、７Ｋ…畳み込み演算、８０…離散逆フーリエ変換と重複加算
【図１】
2015210512.tif 000026
【図２】
2015210512.tif 000027
【図３】
2015210512.tif 000028
【図４】
2015210512.tif 000029
【図５】
2015210512.tif 000030
【図６】
2015210512.tif 000031
【図７】
2015210512.tif 000032
【図８】
2015210512.tif 000033
【図９】
2015210512.tif 000034
【図１０】
2015210512.tif 000035