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