(57)【特許請求の範囲】
【請求項１】
未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみからブラインドで推定する方法であって、観測信号に窓関数を乗算した後、これに点数Ｋの短時間離散フーリエ変換を適用して得られる観測信号ベクトルに周波数ビン毎に忘却係数によって指数関数的に重みを付け、重み付けされたベクトルからブロック処理に基づき観測信号の空間相関行列の逆行列をブロック毎に求め、これを正規化した後、１つのブロック又は複数のブロックを当該周波数ビンの最小２乗型同時対角化問題の対象として、最小２乗法、累乗法、反復法の組み合わせによる解法を用いて最小２乗型同時対角化問題の解を求め、同一信号源から発生した信号の隣接周波数ビンに相関があることを利用してパーミュテーション問題を信号分離精度に基づき適応的に解法することによって信号分離精度の高い分離行列を求めることを特徴とするブラインド信号分離方法。
【請求項２】
請求項１記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、累乗法に含まれる正規化によって分離行列の各列ベクトルを正規化するスケーリング問題を代用することができることを特徴とするブラインド信号分離方法。
【請求項３】
請求項１乃至請求項２のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、短いブロック長でセンサの個数が多い応用においてブロック毎の逆行列の演算量が多い場合には、周波数ビン毎に指数関数的に重み付けられた観測信号の空間相関行列の算出とブロック毎の逆行列の演算の代わりに逆行列の補題を用いることを特徴とするブラインド信号分離方法。
【請求項４】
請求項１乃至請求項３のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、ブロック長がＬでＭ個のブロックを対象区間とした最小２乗型同時対角化問題の解法をｎ回適用して分離行列を推定した後、同一周波数ビンにおける各分離信号の相関に基づき、周波数ビンωｋ（ｋ＝０、１、…、Ｋ−１）におけるＮ個の分離信号のペリオドグラムの逆数の２乗和の平方根による正規化を用いて当該周波数ビンにおける信号分離精度μ（ωｋ）（ｎ）を式（１）で求めることを特徴とするブラインド信号分離方法。Ｔｒ（・）は行列のトレースを表す。また、Ｑは行列の各行に１となる要素が１箇所、その他の要素は０で、１となる要素の位置が他の行と重複しない行列の集合である。行列Ｅ（ωｋ）（ｎ）と行列Γ−１（ωｋ，ｎΨＬ＋τＬ）はそれぞれ式（５）、（６）で与えられ、τ＝１、２、…、ＭにおいてΓ−１（ωｋ，ｎΨＬ＋τＬ）は、分離信号のペリオドグラムλｊ（ωｋ，ｎΨＬ＋τＬ）の逆数を対角要素とする対角行列Λ−１（ωｋ，ｎΨＬ＋τＬ）を正規化した行列で、対角要素γｊ−１（ωｋ，ｎΨＬ＋τＬ）を有する。尚、Ｍ個のブロックが最小２乗型同時対角化問題の対象区間であるので、Ψは連続する２組の対象区間の重複シフトサイズとなる。
【数１】
0005509481.tif 000048
【請求項５】
請求項１乃至請求項４のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、信号分離精度が著しく低い周波数ビンを検出して、該当周波数ビンを間引いた後、これらの周波数ビンを詰めてパーミュテーション問題を適応的に解法することを特徴とするブラインド信号分離方法。
【請求項６】
請求項１乃至請求項５のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、信号分離精度が著しく低い周波数ビンを詰めた状態で同一信号源から発生した信号の隣接または近接周波数ビンに相関があることを利用してパーミュテーション問題を適応的に解法することを特徴とするブラインド信号分離方法。
【請求項７】
請求項１乃至請求項６のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、信号分離精度が著しく低い周波数ビンに該当するパーミュテーションに、近接周波数ビンのパーミュテーションをコピーすることを特徴とするブラインド信号分離方法。
【請求項８】
請求項１乃至請求項７のいずれか１項に記載のブラインド信号分離方法を用いて信号源分離を行うように構成されていることを特徴とするブラインド信号分離方法を用いたブラインド信号分離装置。
【発明の詳細な説明】
【技術分野】
【０００１】  本発明は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみから推定するブラインド信号分離方法に係わり、特に、最小２乗型同時対角化問題の解法を用いて高い精度で信号を分離することができるバックワードモデル型ブラインド信号分離方法及びバックワードモデル型ブラインド信号分離装置に関する。
【背景技術】
【０００２】  複数の未知信号源信号が未知の畳み込み混合系により混在されて観測されるとき、観測信号を分離して混在前の未知信号源信号を推定する処理をブラインド信号分離という。ブラインド信号分離では、未知信号源信号間の統計的独立性のみを条件として、観測信号から未知信号源信号を推定する方法であり、信号源の位置或いは観測信号の到来方向の推定を必ずしも必要としない方法である。【０００３】  勾配法に基づく周波数領域の信号分離方法が非特許文献１で提案されている。勾配法に基づき周波数ビン毎に観測信号を分離した後、パーミュテーション問題を解くことによって信号源を割り付ける。【０００４】  非特許文献２により観測信号の相関行列の特異値分解に基づく信号分離と線形予測分析を組み合わせた方法も提案されている。この方法は線形予測フィルタによって信号分離精度の向上を目的としている。【０００５】  最小２乗型同時対角化問題の解法を用いたフォワードモデル型ブラインド信号分離方法が非特許文献３で提案されている。この方法は、最小２乗型同時対角化問題の解法による混合行列の推定、最小２乗型一般化逆行列を用いた混合行列からの分離行列の推定、スケーリング問題の解法、パーミュテーション問題の解法の４つの手順から成る。【０００６】
【非特許文献１】  Ｓ．Ａｒａｋｉ，Ｒ．Ｍｕｋａｉ，Ｓ．Ｍａｋｉｎｏ，Ｔ．Ｎｉｓｈｉｋａｗａ，ａｎｄ  Ｈ．Ｓａｒｕｗａｔａｒｉ，“Ｔｈｅ  ｆｕｎｄａｍｅｎｔａｌ  ｌｉｍｉｔａｔｉｏｎ  ｏｆ  ｆｒｅｑｕｅｎｃｙ  ｄｏｍａｉｎ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ  ｆｏｒ  ｃｏｎｖｏｌｕｔｉｖｅ  ｍｉｘｔｕｒｅｓ  ｏｆ  ｓｐｅｅｃｈ”，ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｓｐｅｅｃｈ  ａｎｄ  Ａｕｄｉｏ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．１１，ｎｏ．２，ｐｐ．１０９−１１６，Ｍａｒ．２００３
【非特許文献２】  中川、高橋、”漏話音声の聞取り防止に重点を置いたＢＳＳ後処理法”，電子情報通信学会論文誌，ｖｏｌ．Ｊ９０−Ａ，ｎｏ．８，ｐｐ．６３３−６４５，２００７
【非特許文献３】  Ｋ．Ｒａｈｂａｒ  ａｎｄ  Ｊ．Ｐ．Ｒｅｉｌｌｙ，”Ａ  ｆｒｅｑｕｅｎｃｙ  ｄｏｍａｉｎ  ｍｅｔｈｏｄ  ｆｏｒ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ  ｏｆ  ｃｏｎｖｏｌｕｔｉｖｅ  ａｕｄｉｏ  ｍｉｘｔｕｒｅｓ”，ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｓｐｅｅｃｈ  ａｎｄ  Ａｕｄｉｏ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．１３，ｎｏ．５，ｐｐ．８３２−８４４，Ｓｅｐｔ，２００５
【発明の概要】
【発明が解決しようとする課題】
【０００７】  従来のブラインド信号分離方法で使用されている勾配法は演算量が少ないが、収束が遅く、実環境下で十分な信号分離性能を得ることができない。【０００８】  また、特異値分解による信号分離と線形予測分析を組み合わせた方法では、未知信号源信号の品質と信号分離性能にトレードオフの関係が存在するので、両立は難しい。【０００９】  更に、最小２乗型同時対角化問題の解法を用いたフォワードモデル型ブラインド信号分離方法では、混合行列の最小２乗型一般化逆行列を用いた分離行列の推定に加法性雑音が考慮されていないため、雑音の影響を無視できない。【００１０】  同一信号源から発生した信号の隣接周波数ビンに相関があることを利用した従来のパーミュテーション問題の解法では、１度誤りが生じると、これ以降の解法が誤り続ける確率が非常に高くなる。【００１１】  本発明はこのような事情を鑑みてなされたものであり、観測信号の空間相関行列の逆行列を最小２乗型同時対角化問題の対象にすることによって信号分離性能の向上と手順の半減の両立、更に、信号分離精度が低い周波数ビンを除き、隣接周波数ビンの相関を利用してパーミュテーション問題を適応的に解法することを目的とする。
【課題を解決するための手段】
【００１２】  このような目的に応えるために本発明（請求項１記載の発明）に係るブラインド信号分離方法は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみからブラインドで推定する方法であって、観測信号に窓関数を乗算した後、これに点数Ｋの短時間離散フーリエ変換を適用して得られる観測信号ベクトルに周波数ビン毎に忘却係数によって指数関数的に重みを付け、重み付けされたベクトルからブロック処理に基づき観測信号の空間相関行列の逆行列をブロック毎に求め、これを正規化した後、１つのブロック又は複数のブロックを当該周波数ビンの最小２乗型同時対角化問題の対象として、最小２乗法、累乗法、反復法の組み合わせによる解法を用いて最小２乗型同時対角化問題の解を求め、同一信号源から発生した信号の隣接周波数ビンに相関があることを利用してパーミュテーション問題を信号分離精度に基づき適応的に解法することによって信号分離精度の高い分離行列を求めることを特徴とするブラインド信号分離方法。【００１３】  本発明（請求項２記載の発明）に係るブラインド信号分離方法は、累乗法に含まれる正規化によって分離行列の各列ベクトルを正規化するスケーリング問題を代用することができることを特徴とする。【００１４】  本発明（請求項３記載の発明）に係るブラインド信号分離方法は、短いブロック長でセンサの個数が多い応用においてブロック毎の逆行列の演算量が多い場合には、周波数ビン毎に指数関数的に重み付けられた観測信号の空間相関行列の算出とブロック毎の逆行列の演算の代わりに逆行列の補題を用いることを特徴とする。【００１５】  本発明（請求項４記載の発明）に係るブラインド信号分離方法離は、ブロック長がＬでＭ個のブロックを対象区間とした最小２乗型同時対角化問題の解法をｎ回適用して分離行列を推定した後、同一周波数ビンにおける各分離信号の相関に基づき、周波数ビンωｋ（ｋ＝０、１、…、Ｋ−１）におけるＮ個の分離信号のペリオドグラムの逆数の２乗和の平方根による正規化を用いて当該周波数ビンにおける信号分離精度μ（ωｋ）（ｎ）を式（４）で求めることを特徴とする。Ｔｒ（・）は行列のトレースを表す。また、Ｑは行列の各行に１となる要素が１箇所、その他の要素は０で、１となる要素の位置が他の行と重複しない行列の集合である。行列Ｅ（ωｋ）（ｎ）と行列Γ−１（ωｋ，ｎΨＬ＋τＬ）はそれぞれ式（５）、（６）で与えられ、τ＝１、２、…、ＭにおいてΓ−１（ωｋ，ｎΨＬ＋τＬ）は、分離信号のペリオドグラムλｊ（ωｋ，ｎΨＬ＋τＬ）の逆数を対角要素とする対角行列Λ−１（ωｋ，ｎΨＬ＋τＬ）を正規化した行列で、対角要素γｊ−１（ωｋ，ｎΨＬ＋τＬ）を有する。尚、Ｍ個のブロックが最小２乗型同時対角化問題の対象区間であるので、Ψは連続する２組の対象区間の重複シフトサイズとなる。
【数２】
0005509481.tif 000002
【００１６】  本発明（請求項５記載の発明）に係るブラインド信号分離方法は、信号分離精度が著しく低い周波数ビンを検出して、該当周波数ビンを間引いた後、これらの周波数ビンを詰めてパーミュテーション問題を適応的に解法することを特徴とする。【００１７】本発明（請求項６記載の発明）に係るブラインド信号分離方法は、信号分離精度が著しく低い周波数ビンを詰めた状態で同一信号源から発生した信号の隣接または近接周波数ビンに相関があることを利用してパーミュテーション問題を適応的に解法することを特徴とする。【００１８】  本発明（請求項７記載の発明）に係るブラインド信号分離方法は、信号分離精度が著しく低い周波数ビンに該当するパーミュテーションに、近接周波数ビンのパーミュテーションをコピーすることを特徴とする。【００１９】  すなわち、本発明によれば観測信号の空間相関行列の逆行列を最小２乗型同時対角化問題の対象にすることにより、最小２乗法、累乗法、反復法を組み合わせることによって最小２乗型同時対角化問題の近似解、即ち、分離行列を適応的に求めることができ、スケーリング問題を解法することなく信号分離性能の向上と演算量の削減を両立させる。【００２０】  また、本発明によれば、低い信号分離精度の周波数ビンを除き、近接周波数ビンの相関を利用してパーミュテーション問題を解法し、信号分離性能を向上させる。
【発明の効果】
【００２１】  本発明によれば、忘却係数を用いて指数関数的に重み付けられた観測信号の空間相関行列の逆行列を求め、これを最小２乗型同時対角化問題に適用する。最小２乗型同時対角化問題の解法では、スケーリング問題を解法することなく最小２乗法、累乗法、反復法を組み合わせることによって最小２乗型同時対角化問題の近似解、即ち、分離行列を適応的に求める。最小２乗型同時対角化問題を適用した従来のフォワードモデル型ブラインド信号分離が混合行列の推定、混合行列の最小２乗型一般化逆行列を用いた分離行列の推定、スケーリング問題の解法、パーミュテーション問題の解法の４つの手順から成っていたものを、本発明に係るブラインド信号分離方法及びブラインド信号分離装置では、バックワードモデル型を用いて分離行列の各列ベクトルを正規化するスケーリング問題の解法を累乗法に含まれる正規化で代用することによって分離行列の推定とパーミュテーション問題の解法の２つの手順に半減できるという効果がある。【００２２】  また、同一信号源から発生した信号の隣接周波数ビンに相関があることを利用してパーミュテーション問題を解法する前に、信号分離精度が著しく低い周波数ビンを検出して当該周波数ビンを間引いた後、これらの周波数ビンを詰めてパーミュテーション問題を適応的に解法する。信号分離精度が著しく低い周波数ビンに該当するパーミュテーションに、近接周波数ビンのパーミュテーションをコピーすることによって、本発明に係るブラインド信号分離方法及びブラインド信号分離装置には、未知信号源信号の品質に影響を与えることなく信号分離性能を高めるという効果がある。【００２３】  更に、本発明に係るブラインド信号分離方法及びブラインド信号分離装置では観測信号の空間相関行列の逆行列を最小２乗型同時対角化問題の対象にすることによって雑音の影響を最小にする分離行列を推定できるという効果がある。
【発明を実施するための最良の形態】
【００２４】  本発明に係るブラインド信号分離方法の実施の形態について図面を参照して説明する。【００２５】１．畳み込み混合モデル
図１に示すように、時刻ｔにおいてＮ個の信号源１１、１２、…、１Ｎから発せられた信号源信号ｓｊ（ｔ）が畳み込み混合されてｘｉ（ｔ）として観測される。ｓｊ（ｔ）は平均０で互いに統計的独立な非定常信号である。また、ｎｉ（ｔ）はセンサ２ｉに加わる平均０、分散σ２のガウス性白色雑音で、ｓｊ（ｔ）と統計的独立である。時刻ｔにおいてＪ個のセンサ２１、２２、…、２Ｊで観測される観測信号ｘｉ（ｔ）は式（７）で表される。ここで、Ｊ≧Ｎ≧２とする。【００２６】
【数３】
0005509481.tif 000003
ここで、＊は畳み込み演算、ｈｉｊ（ｔ）は信号源１ｊからセンサ２ｉまでの信号経路のインパルス応答、Ａはインパルス応答長をそれぞれ表す。観測信号ｘｉ（ｔ）を連続して３１で短時間離散フーリエ変換すると、フレームｍにおける観測信号は式（８）により表される。【００２７】
【数４】
0005509481.tif 000004
上式（８）において、ｗｉｎ（ｔ）は窓関数、Ｋは短時間離散フーリエ変換の点数、Ｔｂ５は２つの重複窓間のシフトサイズ、ωｋ＝２πｋ／Ｋ、ｋ＝０，１，…，Ｋ−１、をそれぞれ表す。ｈｉｊ（ｔ）のＫ点離散フーリエ変換をｈｉｊ（ωｋ）、ｓｊ（ｔ）に窓関数を乗算した後、Ｋ点短時間離散フーリエ変換で周波数領域に変換したフレームｍの信号源信号をｓｊ（ωｋ，ｍ）、同様に、ｎｉ（ｔ）に窓関数を乗算した後、Ｋ点短時間離散フーリエ変換で周波数領域に変換したフレームｍの雑音をｎｉ（ωｋ，ｍ）とそれぞれ表記すると、点数Ｋがｈｉｊ（ｔ）のインパルス応答長Ａより十分に大きいとき観測信号は式（９）により近似される。【００２８】
【数５】
0005509481.tif 000005
上式（９）において、Ｓ（ωｋ，ｍ）はフレームｍに各信号源から出力される信号ベクトルである。信号ベクトルＳ（ωｋ，ｍ）、混合ベクトルｈｉ（ωｋ）はそれぞれ式（１０）、（１１）により表される。【００２９】
【数６】
0005509481.tif 000006
【００３０】  観測信号ベクトルＸ（ωｋ，ｍ）、混合行列Ｈ（ωｋ）、雑音ベクトルＮ（ωｋ，ｍ）をそれぞれ式（１２）、（１３）、（１４）により定義すると、Ｘ（ωｋ，ｍ）は式（１５）によって表される。【００３１】
【数７】
0005509481.tif 000007
【００３２】  信号を分離するには５１、５２、…、５Ｋで周波数ビン毎に式（１６）を満足する分離行列Ｗ（ωｋ）を推定し、６０で信号源の割り当てを定めるパーミュテーション行列Π（ωｋ）を決定する。周波数ビン毎に独立にパーミュテーション行列Π（ωｋ）を決定しても信号が完全に分離する保証はなく、同一信号源から発生した信号の隣接または近接周波数ビンに相関があることを利用してパーミュテーション行列Π（ωｋ）を決定する。【００３３】
【数８】
0005509481.tif 000008
上式（１６）においてＤ（ωｋ）は周波数ビン毎に異なる任意の対角行列である。【００３４】  スケーリング問題とパーミュテーション問題を順に解法した後、７１、７２、…、７ＫでＸ（ωｋ，ｍ）に左から分離行列Ｗ（ωｋ）を乗算すると、周波数ビンωｋにおける分離信号Ｙ（ωｋ，ｍ）は式（１７）で表される。尚、スケーリング問題の解法については後述する。【００３５】
【数９】
0005509481.tif 000009
上式（１７）を８０で逆離散フーリエ変換と重複加算によって時間領域に変換すると分離信号ｙｉ（ｔ）が求められる。雑音の分散σ２が十分に小さいとき、ｙｉ（ｔ）≒ｓｉ（ｔ）になる。尚、分離信号ベクトルＹ（ωｋ，ｍ）は式（１８）により表される。【００３６】
【数１０】
0005509481.tif 000010
【００３７】  本発明に関するブラインド信号分離方法について、図１乃至図４を参照して詳細に説明する。図２乃至図４は、図１の３１における短時間離散フーリエ変換後、ブラインド信号分離システム４０において本発明により周波数ビン毎に推定される分離行列の算出手順を示したものである。図４は、分離行列の算出後、ブラインド信号分離システム４０において本発明によりパーミュテーション行列の算出手順を示したものである。【００３８】２．最小２乗型同時対角化問題
ブロック長Ｌが短く、センサの個数Ｊが多い場合、逆行列の演算と逆行列の算出頻度が増加してフレーム当たりの演算量は増加する。そこで、ステップＳ１０１においてＬが短くＪが大きい場合、ステップＳ１０２で式（１９）の逆行列の補題を用いて逆行列Ｐ−１（ωｋ，ｍ）を再帰的に求め、ブロック処理に基づきＬフレーム毎に記憶したＰ−１（ωｋ，ｎΨＬ＋τＬ）をステップＳ１０５で式（２１）によって正規化する。ただし、Ｉは単位行列、ｃは正の定数、Ｐ−１（ωｋ，０）＝ｃ−１Ｉである。忘却係数βの値は０＜β≦１から選択される。また、Ｍ個のブロックを最小２乗型同時対角化問題の対象区間とすると、Ψは連続する２組の対象区間の重複シフトサイズを表し、［１、Ｍ］から選択される。更に、τ＝１，２，…，Ｍ、ｎは最小２乗型同時対角化問題の適用回数を表し、ｎ＝０，１，２，…とする。【００３９】
【数１１】
0005509481.tif 000011
【００４０】  ステップＳ１０１においてＬが長い、またはＪが少ない場合、忘却係数βを導入してＸ（ωｋ，ｍ）ＸＨ（ωｋ，ｍ）を指数関数的に重み付け、行列Ｐ（ωｋ，ｍ）をステップＳ１０３において式（２０）に基づき求める。【００４１】
【数１２】
0005509481.tif 000012
ここで、上付き添字Ｈは複素共役転置をそれぞれ表す。また、初期値をＰ（ωｋ，０）＝ｃＩのように設定する。ステップＳ１０４においてブロック処理に基づきＬフレーム毎にＰ（ωｋ，ｎΨＬ＋τＬ）の逆行列Ｐ−１（ωｋ，ｎΨＬ＋τＬ）を求め、ステップＳ１０５で式（２１）によって正規化する。【００４２】
【数１３】
0005509481.tif 000013
ここで、‖・‖Ｆはフロベニウスノルムを表す。【００４３】２．１  対角化行列の解法
ブロック処理を導入して、ブロック長がＬ個のフレームで、Ｍ個のブロックを最小２乗型同時対角化問題の対象区間にする。【００４４】  ‖Ｗｉ（ωｋ）（ｎ）‖２＝１を条件に式（２２）の評価関数が最小になるように対角化行列Ｗ（ωｋ）（ｎ）、即ち、分離行列とＭ個の対角行列Λ−１（ωｋ，ｎΨＬ＋τＬ）を最小２乗法、累乗法、反復法を併用して求める。ただし、Ｗｉ（ωｋ）（ｎ）はＷＨ（ωｋ）（ｎ）のｉ番目の列ベクトル、‖・‖２はユークリッドノルムをそれぞれ表す。式（２２）は最小２乗型同時対角化問題として知られている。上付き添字（ｎ）は最小２乗型同時対角化問題の解法の適用回数を表す。【００４５】
【数１４】
0005509481.tif 000014
【００４６】  式（２２）のように‖Ｗｉ（ωｋ）（ｎ）‖２＝１を条件に最小２乗型同時対角化問題を解法することは分離行列ＷＨ（ωｋ）（ｎ）の各列ベクトルＷｉ（ωｋ）（ｎ）を正規化するスケーリング問題の解法を同時に実行していることと等価であり、別途、スケーリング問題を解法する必要はない。条件‖Ｗｉ（ωｋ）（ｎ）‖２＝１は最小２乗型同時対角化問題の解法で使用する累乗法の正規化によって満たされる。【００４７】  ベクトル表現を用いると、式（２２）の評価関数は式（２３）のように表現することができる。ここで、ｒ−１（ωｋ，ｎΨＬ＋τＬ）、Ｇ（ωｋ）（ｎ）、ｄ（ωｋ，ｎΨＬ＋τＬ）、Ｇ（ωｋ）（ｎ）ｄ（ωｋ，ｎΨＬ＋τＬ）はそれぞれ式（２４）〜（２７）により表される。ただし、ｖｅｃ｛Ａ｝は行列Ａの列を積み重ねることによって
0005509481.tif 000015
【００４８】
【数１５】
0005509481.tif 000016
【００４９】  ｇｉ（ωｋ）（ｎ）∈Ωを条件に式（２３）の評価関数を最小にするＧ（ωｋ）（ｎ）とｄ（ωｋ，ｎΨＬ＋τＬ）は、式（２２）の評価関数を最小にする分離行列Ｗ（ωｋ）（ｎ）と対角行列Λ−１（ωｋ，ｎΨＬ＋τＬ）に等しい。ここで、ｇｉ（ωｋ）（ｎ）はＧ（ωｋ）（ｎ）のｉ番目の列ベクトルである。また、空間Ωは式（２８）により表される。ただし、ＣＪ×１はＪ×１の複素空間を表す。【００５０】
【数１６】
0005509481.tif 000017
【００５１】  反復法によって分離行列Ｗ（ωｋ）（ｎ）と対角行列Λ−１（ωｋ，ｎΨＬ＋τＬ）を求めるためにｚｉ（ωｋ）（ｎ）とＴ（ωｋ）（ｎ）を式（２９）、（３０）によって定義する。【００５２】
【数１７】
0005509481.tif 000018
【００５３】  ｚｉ（ωｋ）（ｎ）とＴ（ωｋ）（ｎ）を用いると、式（２３）は式（３１）のように表現できる。【００５４】
【数１８】
0005509481.tif 000019
【００５５】  ｇｉ（ωｋ）（ｎ）とｚｉ（ωｋ）（ｎ）をステップＳ１０６において式（３２）、（３４）により初期化する。【００５６】
【数１９】
0005509481.tif 000020
ただし、ｇｊ（ωｋ）（ｎ）は最小２乗型同時対角化問題の解法をｎ回適用して求められた近似解である。また、初期値ｇｊ（ωｋ）（０）は式（３３）で設定する。【００５７】
【数２０】
0005509481.tif 000021
ここで、ａを任意の実数定数、１Ｊ＝［１，１，…，１］Ｔである。【００５８】
【数２１】
0005509481.tif 000022
ただし、ｚｊ（ωｋ）（ｎ）は最小２乗型同時対角化問題の解法をｎ回適用して求められた近似解である。また、初期値ｚｊ（ωｋ）（０）は式（３５）で設定する。【００５９】
【数２２】
0005509481.tif 000023
ここで、ｂを任意の実数定数、１Ｍ＝［１，１，…，１］Ｔである。【００６０】  ｇｊ（ωｋ）（ｎ）を求める際、ｇｊ（ωｋ）（ｎ）（ｊ≠ｉ）を式（３２）のように、ｚｊ（ωｋ）（ｎ）を式（３４）のようにそれぞれ定数として設定して、Ｆｊ（ωｋ）（ｎ）を計算する。【００６１】
【数２３】
0005509481.tif 000024
【００６２】  式（３７）の評価関数を最小にする近似解ｇｊ（ωｋ）（ｎ，φ）の算出時には、ステップＳ１０７において対象外のｇｊ（ωｋ）（ｎ）とｚｉ（ωｋ）（ｎ）（ｉ≠ｊ）を定数として固定する。近似解ｇｊ（ωｋ）（ｎ，φ）は反復によって求めるため、反復対象ｇｊ（ωｋ）（ｎ）に反復回数φを上付き添字で示している。【００６３】
【数２４】
0005509481.tif 000025
ｇｊ（ωｋ）（ｎ，φ）∈Ωを条件に式（３７）を最小化する。条件ｇｊ（ωｋ）（ｎ，φ）∈Ω
0005509481.tif 000026
はステップＳ１０８において式（３８）により求められる。【００６４】
【数２５】
0005509481.tif 000027
空間Ωで最小２乗法により式（３８）の解を求めると、ｇｊ（ωｋ）（ｎ，φ）∈Ωを条件に式（３７）を最小化する解ｇｊ（ωｋ）（ｎ，φ）と一致する。【００６５】
【数２６】
0005509481.tif 000028
上式（３９）においてＶ（ωｋ）（ｎ，φ）、Ｗｊ（ωｋ）（ｎ，φ）ＷｊＨ（ωｋ）（ｎ，φ）はそれぞれ式（４０）、（４１）により表される。【００６６】
【数２７】
0005509481.tif 000029
ここで、ｍａｔ｛Ａ｝は、Ｊ２×１の列ベクトルＡをＪ×Ｊの行列に変換することを表す。【００６７】
【数２８】
0005509481.tif 000030
Ｗｊ（ωｋ）（ｎ，φ）ＷｊＨ（ωｋ）（ｎ，φ）の列ベクトルによって張られる空間の次元はｄｉｍ（Ｗｊ（ωｋ）（ｎ，φ）ＷｊＨ（ωｋ）（ｎ，φ））＝１であるので、ステップＳ１０９において、累乗法によって求めたＶｊ（ωｋ）（ｎ，φ）の第１固有ベクトルはＷｊ（ωｋ）（ｎ，φ）に一致する。【００６８】  累乗法の初期値をｕ（０）と表記すると、誤差の限界εＰが式（４４）を満足するまで、ｌ＝１から式（４２）〜（４４）の累乗法を繰り返し、近似解Ｗｊ（ωｋ）（ｎ，φ）＝ｕ（１）を得る。【００６９】
【数２９】
0005509481.tif 000031
【００７０】
0005509481.tif 000032
まれる正規化によって分離行列の各列ベクトルを正規化するスケーリング問題を解いていることになる。【００７１】  ステップＳ１１０においてｊ＝１，２，…，Ｎの順にｊを変えて式（３６）、（３８）、累乗法によるＶｊ（ωｋ）（ｎ，φ）の第１固有ベクトルの計算を繰り返す。【００７２】  ステップＳ１１１において誤差の限界がεＧの近似解を求めるために、式（４５）を満足するまでｊ＝１，２，・・・，Ｎの順にｊを変えて式（３６）、（３８）、累乗法によるＶｊ（ωｋ）（ｎ）の第１固有ベクトルの計算を繰り返す。【００７３】
【数３０】
0005509481.tif 000033
【００７４】２．２  対角行列の解法
式（４６）を最小にすることによって対角行列は求められる。【００７５】
【数３１】
0005509481.tif 000034
【００７６】  式（２４）、（２６）、（２７）を用いると、式（４６）はベクトル表現によって式（４７）のように表される。【００７７】
【数３２】
0005509481.tif 000035
【００７８】  式（４７）の評価関数に最小２乗法を適用して、ステップＳ１１２において式（４８）で解を求める。【００７９】
【数３３】
0005509481.tif 000036
ここで、Ｇ＋（ωｋ）（ｎ）は行列Ｇ（ωｋ）（ｎ）の擬似逆行列である。【００８０】  ステップＳ１１３において誤差の限界がεＪの近似解を求めるために、式（４９）を満足するまで式（３６）、（３８）、Ｖｊ（ωｋ）（ｎ）の第１固有ベクトルによるｇｊ（ωｋ）（ｎ）の計算、式（４８）によるｄ（ωｋ，ｎΨＬ＋τＬ）の計算を繰り返す。ただし、｜・｜は絶対値を表す。【００８１】
【数３４】
0005509481.tif 000037
【００８２】３．パーミュテーション問題の解法
周波数ビンωｋにおける信号分離精度を評価するために、分離信号のペリオドグラムλｊ（ωｋ，ｎΨＬ＋τＬ）の逆数λｊ−１（ωｋ，ｎΨＬ＋τＬ）を用いてΛ−１（ωｋ，ｎΨＬ＋τＬ）を式（５０）によってステップＳ１１４で構成する。ただし、ｄｉａｇ（・）は対角行列を表す。【００８３】
【数３５】
0005509481.tif 000038
【００８４】  分離信号の周波数成分は周波数ビン毎に異なるので、この影響を低減するために、ステップＳ１１５でΛ−１（ωｋ，ｎΨＬ＋τＬ）を式（５１）によって正規化する。γｊ−１（ωｋ，ｎΨＬ＋τＬ）はΓ−１（ωｋ，ｎΨＬ＋τＬ）のｊ行ｊ列の対角要素とする。【００８５】
【数３６】
0005509481.tif 000039
【００８６】  周波数ビン間の相関に基づき信号分離精度を評価するためにＥ（ωｋ）（ｎ）を式（５２）によりステップＳ１１６で算出する。【００８７】
【数３７】
0005509481.tif 000040
【００８８】  周波数ビンωｋにおける信号分離精度μ（ωｋ）（ｎ）を式（５３）によりステップＳ１１７で計算する。【００８９】
【数３８】
0005509481.tif 000041
ここで、Ｔｒ（・）は行列のトレースを表す。また、Ｑは行列の各行に１となる要素が１箇所、その他の要素は０で、１となる要素の位置が他の行と重複しない行列の集合である。【００９０】  ステップＳ１１８において全周波数ビンを対象に式（５３）を用いて信号分離精度μ（ωｋ）（ｎ）を求め、信号分離精度の下位ｘ％を検出する。【００９１】  ステップＳ１１９において下位ｘ％に該当した周波数ビンを間引き、該当周波数ビンは間を詰める。【００９２】  間引いた状態で、隣接周波数ビンにおいて既に定められた信号源の割り当て結果Ｅη（ωη±｜φ｜）（ｎ）をステップＳ１２０で式（５４）により求める。尚、基準周波数ビンωηにおけるパーミュテーション行列はΠ（ωη）（ｎ）＝Ｉに設定する。【００９３】
【数３９】
0005509481.tif 000042
【００９４】  隣接する周波数ビンωη±｜φ＋１｜におけるパーミュテーション行列Π（ωη±｜φ＋１｜）（ｎ）を式（５５）によりステップＳ１２１で決定する。【００９５】
【数４０】
0005509481.tif 000043
【００９６】  ステップＳ１２２において下位ｘ％の周波数ビンを除く、他の全ての周波数ビンを対象に式（５４）、（５５）を繰り返し、Π（ωη±｜φ＋１｜）（ｎ）を求める。【００９７】  信号分離精度の下位ｘ％に該当する周波数ビンのパーミュテーション行列は、ステップＳ１２３において最も近い周波数ビンのパーミュテーション行列をコピーする。【００９８】  ステップＳ１２４において式（５６）のように観測信号Ｘ（ωｋ，ｍ）に分離行列Ｗ（ωｋ）（ｎ）を左から乗算した後、更にパーミュテーション行列Π（ωｋ）（ｎ）を左から乗算して分離信号Ｙ（ωｋ，ｍ）を得る。【００９９】
【数４１】
0005509481.tif 000044
【実施例１】
【０１００】４．１  評価データ
図５のように、３個のセンサ（マイクロホン）２１、２２、２３と半径１．５ｍの円の円周上に３個の信号源（スピーカ）１１、１２、１３を残響時間のある会議室に配置した。尚、図５は３個の信号源（スピーカ）１１、１２、１３と３個のセンサ（マイクロホン）２１、２２、２３の位置関係を示す平面図である。信号源１ｊとセンサ２ｉの間の伝達関数ｈｉｊ（ｔ）は、残響時間を１００ｍｓに設定し、信号源（スピーカ）から発したＴＳＰ信号をセンサ（マイクロホン）で受音することによって測定した。音声データは男性３名と女性２名から異なる３名を選択して、３分間の混合データを２組作成した。実験条件は、信号源（スピーカ）とセンサ（マイクロホン）の個数はＪ＝Ｎ＝３、標本化周波数８ｋＨｚ、量子化ビット数１６ビット、短時間離散フーリエ変換の点数Ｋ＝２０４８、シフトサイズＴｂ＝３０、窓関数はハニング窓である。忘却係数β＝０．９９、ブロック長Ｌ＝３２フレーム、Ψ＝１ブロックでＭ＝１５００ブロックを最小２乗型同時対角化問題の対象にした。本発明に係るブラインド信号分離方法では、初期値ａ＝ｂ＝１、εＧ＝１０−５、εＪ＝１０−５、εＰ＝１０−１５を用いている。パーミュテーション問題の解法において基準周波数ビンはｆη＝１．５ｋＨｚに設定している。また、信号分離精度の下位２０％または下位５０％に該当する周波数ビンは間引いている。ＳＮＲは１０ｄＢ間隔で１０〜３０ｄＢの範囲で変化させた。ここで、センサ２１、２２、２３のＳＮＲは式（５７）によってそれぞれ算出した。ここで、Ｅ［・］は期待値を表す。【０１０１】
【数４２】
0005509481.tif 000045
【０１０２】４．２  評価指標
ブラインド信号分離方法の信号分離性能を次の方法で評価した．式（５８）によって観測信号における所望信号源信号と干渉信号の電力の比、式（５９）によって出力信号における所望信号源信号と干渉信号の電力の比をそれぞれ計算した後、式（６０）のようにＳＩＲｙｉとＳＩＲｘｉの差によってブラインド信号分離装置の各出力の信号分離性能を求める。式（６１）による各出力の平均を信号分離性能とした。ｃｉｊ（ｔ）（ｎ）は式（６２）のＣ（ωｋ）（ｎ）の要素を、ｗｉｊ（ｔ）（ｎ）はＷ（ωｋ）（ｎ）の要素をそれぞれ逆離散フーリエ変換したものである。【０１０３】
【数４３】
0005509481.tif 000046
【０１０４】４．３  評価対象
最小２乗型同時対角化問題の解法を用いたフォワードモデル型ブラインド信号分離方法（非特許文献３）を比較対象とする。なお、非特許文献３ではＮＳ＝３２、Ｍ＝１５００エポックを最小２乗型同時対角化問題の対象とした。【０１０５】４．４  評価結果
図６の信号分離性能から明らかなように、本発明に係るブラインド信号分離方法が従来のフォワードモデル型ブラインド信号分離方法よりも高いＳＩＲを得ることができた。この要因は信号の統計的性質を反映する為に忘却係数βを導入し、空間相関行列の逆行列を最小２乗型同時対角化問題を適用することによって雑音の影響を最小にすると同時に、分離行列が精度良く推定でき、パーミュテーション問題の正答率が高いためと考えられる。
【実施例２】
【０１０６】５．１  評価データ
２個の信号源（スピーカ）１２、１４と２個のセンサ（マイクロホン）２１、２２を図５のように残響のある会議室に配置した。尚、図５は２個の信号源（スピーカ）１２、１４と２個のセンサ（マイクロホン）２１、２２の位置関係を示す平面図である。２個の信号源（スピーカ）１２、１４は直線上に配置されている。２名の男性による３分間の音声を音声データとして使用した。実験条件は、信号源（スピーカ）とセンサ（マイクロホン）の個数はＪ＝Ｎ＝２、標本化周波数８ｋＨｚ、量子化ビット数１６ビット、短時間離散フーリエ変換の点数Ｋ＝４０９６、シフトサイズＴｂ＝６０、窓関数はハニング窓である。忘却係数β＝０．９９、ブロック長Ｌ＝４０フレーム、Ψ＝１ブロックでＭ＝５００ブロックを最小２乗型同時対角化問題の対象にした。本発明に係るブラインド信号分離方法では、初期値ａ＝ｂ＝１、εＧ＝１０−４、εＪ＝１０−４、εＰ＝１０−１５を用いている。信号源１ｊとセンサ２ｉの間の伝達関数ｈｉｊ（ｔ）は、残響時間を３００ｍｓに設定し、信号源（スピーカ）から発したＴＳＰ信号をセンサ（マイクロホン）で受音することによって測定した。パーミュテーション問題の解法において基準周波数ビンはｆη＝１．５ｋＨｚに設定している。また、信号分離精度の下位１７．４％に該当する周波数ビンは間引いている。【０１０７】５．２  評価指標
式（６３）で最適なパーミュテーション行列Πｏｐｔ（ωｋ）（ｎ）、本発明に係るブラインド信号分離方法でパーミュテーション行列Π（ωｋ）（ｎ）をそれぞれ求め、誤差‖Πｏｐｔ（ωｋ）（ｎ）−Π（ωｋ）（ｎ）‖Ｆを計算する。誤差が零のときは、本発明に係るブラインド信号分離方法で正しくパーミュテーション問題が解法されたことになり、誤差が零でない場合は正しいパーミュテーション行列が求められなかったことになる。【０１０８】
【数４４】
0005509481.tif 000047
【０１０９】５．３  評価対象
同一信号源から発生した信号の隣接周波数ビンに相関があることを利用したパーミュテーション問題の解法（非特許文献３）を比較対象とする。【０１１０】５．４  評価結果
誤差‖Πｏｐｔ（ωｋ）（ｎ）−Π（ωｋ）（ｎ）‖Ｆが零になる周波数ビンを白、零にならない周波数ビンを黒で描くと図７、８の上段のグラフを得る。灰色は白と黒が交互に現れる領域である。図７、８の上段のグラフの縦軸は周波数、横軸は最小２乗型同時対角化問題の解法の適用回数ｎをそれぞれ表す。図中で白の面積が多くなると、パーミュテーション問題の正答率が高くなる。本発明に係るブラインド信号分離方法で適応的にパーミュテーション問題を解法すると、適用回数の増加に伴い白の面積が増大する。図７の下段のグラフより最小２乗型同時対角化問題の解法の適用回数ｎの増加とともに誤り率が減少することから、パーミュテーション問題の正答率が高くなることがわかる。一方、従来のパーミュテーション問題の解法では、図８の上段のグラフより適用回数を増しても、白の面積は増加も減少もしないばかりでなく灰色の領域が多く本発明に係るブラインド信号分離方法に比べ不安定である。また、図８の下段のグラフより誤り率は約５０％である。これはランダムにパーミュテーション行列を割り当てた結果と変わらないことから、本発明に係るブラインド信号分離方法に比べパーミュテーション問題の正答率は向上しないことがわかる。
【実施例３】
【０１１１】６．１  評価データ
２個の信号源（スピーカ）１２、１４と２個のセンサ（マイクロホン）２１、２２を図５のように残響のある会議室に配置した。尚、図５は２個の信号源（スピーカ）１２、１４と２個のセンサ（マイクロホン）２１、２２の位置関係を示す平面図である。１名の男性による３分間の音声を音声データ、３分間のクラシック音楽を音楽データとしてそれぞれ使用した。実験条件は、信号源（スピーカ）とセンサ（マイクロホン）の個数はＪ＝Ｎ＝２、標本化周波数８ｋＨｚ、量子化ビット数１６ビット、短時間離散フーリエ変換の点数Ｋ＝４０９６、シフトサイズＴｂ＝４０、窓関数はハニング窓である。忘却係数β＝０．９９、ブロック長Ｌ＝４０フレーム、Ψ＝１ブロックでＭ＝９００ブロックを最小２乗型同時対角化問題の対象にした。本発明に係るブラインド信号分離方法では、初期値ａ＝ｂ＝１、εＧ＝１０−４、εＪ＝１０−４、εＰ＝１０−１５を用いている。信号源１ｊとセンサ２ｉの間の伝達関数ｈｉｊ（ｔ）は、残響時間を３００ｍｓに設定し、信号源（スピーカ）から発したＴＳＰ信号をセンサ（マイクロホン）で受音することによって測定した。パーミュテーション問題の解法において基準周波数ビンはｆη＝１．５ｋＨｚに設定している。また、信号分離精度の下位２０％に該当する周波数ビンは間引いている。【０１１２】６．２  評価指標
５〜９秒の範囲の信号源信号ｓ１（ｔ）、ｓ２（ｔ）、観測信号ｘ１（ｔ）、ｘ２（ｔ）、本発明に係るブラインド信号分離方法で観測信号ｘ１（ｔ）、ｘ２（ｔ）を分離した信号ｙ１（ｔ）、ｙ２（ｔ）をそれぞれスペクトログラムで表示する。信号源信号ｓ１（ｔ）、ｓ２（ｔ）と分離信号ｙ１（ｔ）、ｙ２（ｔ）については、発音内容と信号波形も示す。【０１１３】６．３  評価結果
図９の上段は信号源信号、中段は観測信号、下段は分離信号のスペクトログラムをそれぞれ示す。本発明に係るブラインド信号分離方法による分離信号のスペクトログラムｙ１（ωｋ，ｍ）とｙ２（ωｋ，ｍ）はそれぞれ、遅延させた信号源信号のスペクトログラムｓ１（ωｋ，ｍ）とｓ２（ωｋ，ｍ）に類似しており、信号が分離できていることが確認できる。
【図面の簡単な説明】
【０１１４】
【図１】  本発明に係るブラインド信号分離方法の実施の形態を示す図である。
【図２】  本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角化行列の解法について説明するためのフローチャートである。
【図３】  本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角行列の解法について説明するためのフローチャートである。
【図４】  本発明に係るブラインド信号分離方法のパーミュテーション行列の解法について説明するためのフローチャートである。
【図５】  本発明に係る信号分離方法の実施例１〜３における信号源（スピーカ）とセンサ（マイクロホン）の位置関係を表す平面図である。
【図６】  本発明に係る信号分離方法の実施例１における信号分離性能を示す図である。
【図７】  本発明に係る信号分離方法の実施例２におけるパーミュテーション問題の解法結果と誤り率を示す図である。
【図８】  実施例２における従来のパーミュテーション問題の解法結果と誤り率を示す図である。
【図９】  本発明に係る信号分離方法の実施例３における発音内容、信号波形、スペクトログラムを示す図である。
【符号の説明】
【０１１５】１１〜１Ｎ…信号源、２１〜２Ｊ…センサ、３１…短時間離散フーリエ変換、４０…ブラインド信号分離システム、５１、５２、…、５Ｋ…最小２乗型同時対角化問題の解法、６０…パーミュテーション問題の解法、７１、７２、…、７Ｋ…畳み込み演算、８０…逆離散フーリエ変換と重複加算
【図１】
0005509481.tif 000049
【図２】
0005509481.tif 000050
【図３】
0005509481.tif 000051
【図４】
0005509481.tif 000052
【図５】
0005509481.tif 000053
【図６】
0005509481.tif 000054
【図７】
0005509481.tif 000055
【図８】
0005509481.tif 000056
【図９】
0005509481.tif 000057