JP2017191309A - ブラインド信号分離方法およびその装置 - Google Patents

Abstract

【課題】観測信号から未知信号源信号を推定するブラインド信号分離方法とその装置を提供する。
【解決手段】窓関数を乗算した観測信号の短時間離散フーリエ変換後、観測信号の相互スペクトル密度行列をエポック時刻毎に求め、混合行列とそのエルミート行列の積、並びに観測信号の相互スペクトル密度行列のエルミート性を利用して、それらの行列の対角要素と対角要素の下に位置する要素のみから制約条件付きフォワードモデル型最小２乗型同時対角化問題を導出する。この問題をラグランジュの未定乗数法により能率的に解法し、混合行列を推定する。次いで、混合行列の条件数と分離信号の振幅スペクトルを用いてパーミュテーションの基準周波数ビンを決定め、同一信号源から発生した信号の周波数ビン間の電力比相関性を用いてパーミュテーション問題を解法し、信号分離性能を向上させる。
【選択図】図１

Description

本発明は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみから推定するブラインド信号分離方法に係わり、特に、最小２乗型同時対角化問題の解法を用いて高い精度で信号を分離することができるブラインド信号分離装置に関する。

複数の未知信号源信号が未知の畳み込み混合系により混在されて観測されるとき、観測信号を分離して混在前の未知信号源信号を推定する処理をブラインド信号分離という。ブラインド信号分離では、未知信号源信号間の統計的独立性のみを条件として、観測信号から未知信号源信号を推定する方法であり、信号源の位置或いは観測信号の到来方向の推定を必ずしも必要としない方法である。

最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法が非特許文献２と特許文献１で提案されている。これらの方法は、最小２乗型同時対角化問題の解法による混合行列の推定、最小２乗型一般化逆行列を用いた混合行列からの分離行列の推定、スケーリング問題の解法、パーミュテーション問題の解法の４つの手順から成る。

周波数ビン間の電力比の相関に基づきパーミュテーション行列を推定する方法が非特許文献３で提案されている。

先行技術論文

非特許文献１

Ｖ．Ｍａｕｒａｎｄｉ ａｎｄ Ｅ．Ｍｏｒｅａｕ，“Ａ ｄｅｃｏｕｐｌｅｄ Ｊａｃｏｂｉ−ｌｉｋｅ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｎｏｎ−ｕｎｉｔａｒｙ ｊｏｉｎｔ ｄｉａｇｏｎａｌｉｚａｔｉｏｎ ｏｆ ｃｏｍｐｌｅｘ−ｖａｌｕｅｄ ｍａｔｒｉｃｅｓ”，ＩＥＥＥ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｌｅｔｔｅｒｓ，ｖｏｌ．２１，ｎｏ．１２，ｐｐ．１４５３−１４５６，Ｄｅｃ．２０１４．
Ｍ．Ｒａｊｉｈ，Ｐ．Ｃｏｍｏｎ，ａｎｄ Ｒ．Ａ．Ｈａｒｓｈｍａｎ，"Ｅｎｈａｎｃｅｄ ｌｉｎｅ ｓｅａｒｃｈ：Ａ ｎｏｖｅｌ ｍｅｔｈｏｄ ｔｏ ａｃｃｅｌｅｒａｔｅ ＰＡＲＡＦＡＣ"，ＳＩＡＭ Ｊ．Ｍａｔｒｉｘ Ａｎａｌ．Ａｐｐｌ．，ｖｏｌ．３０，ｎｏ．２，ｐｐ．１１２８−１１４７，Ｓｅｐ．２００８． Ｈ．Ｓａｗａｄａ，Ｓ．Ａｒａｋｉ，ａｎｄ Ｓ．Ｍａｋｉｎｏ，"Ｍｅａｓｕｒｉｎｇ ｄｅｐｅｎｄｅｎｃｅ ｏｆ ｂｉｎ−ｗｉｓｅ ｓｅｐａｒａｔｅｄ ｓｉｇｎａｌｓ ｆｏｒ ｐｅｒｍｕｔａｔｉｏｎ ａｌｉｇｎｍｅｎｔ ｉｎ ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ ＢＳＳ"Ｐｒｏｃ．ＩＥＥＥ Ｉｎｔ．Ｓｙｍｐ．Ｃｉｒｃｕｉｔｓ Ｓｙｓｔ．，ｐｐ．３２４７−３２５０，Ｍａｙ ２００７ 齋藤 晋哉、大石 邦夫、古川 利博，”ブラインド信号分離方法およびその装置”，特開２０１５−２１０５１２号、２０１５−１１−２４

従来のブラインド信号分離方法で使用されている最小２乗型同時対角化問題の解法は、実環境下で十分な信号分離精度を得ることができるが、解法手順が収束に要する演算量が多く、演算コストが増大することが問題となっている。

同一信号源から発生した信号の隣接周波数ビンに相関があることを利用した従来のパーミュテーション問題の解法では、経験に基づき基準周波数ビンを決めており、適切な周波数ビンが選択されたとは言えない。

同一信号源から発生した信号の隣接周波数ビンに相関があることを利用した従来のパーミュテーション問題の解法では、１度誤りが生じると、これ以降の解法が誤り続ける確率が非常に高くなる。

本発明はこのような事情を鑑みてなされたものであり、最小２乗型同時対角化問題の解法手順を用いた従来のブラインド信号分離方法の信号分離精度を損なうことなく、解法手順が収束に要する演算量の削減、更に、パーミュテーション問題の解法において、分離フィルタの出力振幅スペクトルと条件数を用いて基準周波数ビンを適切に決定することを目的とする。

このような目的に応えるために本発明（請求項１記載の発明）に係るブラインド信号分離方法は、互いに統計的に独立な未知信号源信号と未知の畳み込み混合系により、複数の未知信号源信号が混在した観測信号のみからブラインドで空間伝達関数を推定する方法であって、時間周波数領域においてエポック時刻毎に観測信号ベクトルから相互スペクトル密度行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる相互スペクトル密度行列を求め、そのノルムで全ての同エポック時刻の相互スペクトル密度行列を正規化した後、混合行列とそのエルミート行列の積、並びに相互スペクトル密度行列がエルミート行列であることを利用して、それらの行列の対角要素とその対角要素の下に位置する要素のみから混合行列を推定するために、正規化された相互スペクトル密度行列をベクトルに変換し、このベクトルの次元を削減し、混合行列の列ベクトルのノルムが１になることを制約条件としてフォワードモデル型最小２乗型同時対角化問題をラグランジュの未定乗数法により反復することなく解法することにより次元削減ベクトルを推定した後、ベクトルの次元を復元し、次元復元ベクトルを行列に変換することによって混合行列を推定し、一般化逆行列を用いて混合行列から分離行列を算出し、同一信号源から発生した信号の周波数ビンの電力比に相関があることを利用してパーミュテーション問題を解法するために、指定した周波数帯域における分離フィルタの出力信号の振幅スペクトルの平均値を求め、同周波数帯域における混合行列とそのエルミート行列の行列積によって求められる行列の条件数が最小となる周波数ビンにおいて、振幅スペクトルがその平均値以上となっている場合、この周波数ビンをパーミュテーション行列の推定のための基準周波数ビンに選択することを特徴とする。

本発明（請求項２記載の発明）に係るブラインド信号分離方法は、信号源信号と観測信号間の混合行列を推定する際、混合行列の列ベクトルのノルムが１になることをフォワードモデル型最小２乗型同時対角化問題の制約条件に加え、この問題をラグランジュの未定乗数を導入して反復することなく解法することにより、ベクトルの次元削減との相乗効果によって混合行列の推定に要する演算量を削減することを特徴とする。

本発明（請求項３記載の発明）に係るブラインド信号分離方法は、パーミュテーション問題の解法の際、基準周波数ビンの決定に、推定された混合行列とそのエルミート行列の行列積によって求められる行列の条件数、並びに、分離フィルタの出力信号の振幅スペクトルを使用することを特徴とする。

すなわち、本発明によれば、観測信号の空間相関行列を制約条件付き最小２乗型同時対角化問題の対象に、ラグランジュの未定乗数法によって最小２乗型同時対角化問題の近似解、即ち、混合行列を求めることができ、演算量の低減と推定性能の向上を両立させる。

また、本発明によれば、混合行列を用いた条件数に基づき分離性能を評価し、パーミュテーション問題の基準周波数ビンの選定に利用することにより、全体の信号分離性能を向上させる。

本発明によれば、時間周波数領域においてエポック時刻毎に観測信号ベクトルから相互スペクトル密度行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大と相互スペクトル密度行列を求め、そのノルムで全ての同エポック時刻の相互スペクトル密度行列を正規化した後、混合行列とそのエルミート行列の積、並びに相互スペクトル密度行列がエルミート行列であることを利用して、それらの行列の対角要素と対角要素の下に位置する要素のみから混合行列を推定するために、正規化された相互スペクトル密度行列をベクトルに変換し、このベクトルの次元を削減し、混合行列の列ベクトルのノルムが１になることを制約条件としてフォワードモデル型最小２乗型同時対角化問題に加え、ラグランジュの未定乗数法により解法することにより次元削減ベクトルを推定しているので、演算量を削減できるという効果がある。

また、本発明に係るブラインド信号分離方法及びブラインド信号分離装置では、信号源信号と観測信号間の混合行列を推定する際、混合行列の列ベクトルのノルムが１になることをフォワードモデル型最小２乗型同時対角化問題の制約条件に加え、この問題をラグランジュの未定乗数を導入して反復することなく解法することにより、少ない演算量で混合行列を推定できるという効果がある。

更に、パーミュテーション問題の解法では、パーミュテーション問題の解法の際、基準周波数ビンの決定に、推定された混合行列とそのエルミート行列の行列積によって求められる行列の条件数、並びに、分離フィルタの出力信号の振幅スペクトルを使用することで、パーミュテーション行列の割り当て正答率が向上するという効果がある。

本発明に係るブラインド信号分離方法の実施の形態を示す図である。 本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角化行列のラグランジェ未定乗数法の適用について説明するためのフローチャートである。 本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角行列の解法について説明するためのフローチャートである。 本発明に係るブラインド信号分離方法のパーミュテーション問題の解法について説明するためのフローチャートである。 本発明に係る信号分離方法の実施例１における信号源（スピーカ）とマイクロホンの位置関係を表す平面図である。 本発明に係る信号分離方法の実施例１における部屋の残響時間と信号分離性能の関係を示す図である。 実施例２における従来のパーミュテーション問題の解法によるパーミュテーション行列の割り当て結果を示す図である。 実施例２における本発明に係るブラインド信号分離方法によるパーミュテーション行列の割り当て結果を示す図である。 分離行列の推定に本発明に係るブラインド信号分離方法を使用して、パーミュテーション問題の解法には特許文献１を使用した場合と、本発明に係るブラインド信号分離方法の実施例２におけるパーミュテーション行列の割り当て正答率、信号分離性能、反復回数、音質評価を比較する図である。

本発明に係るブラインド信号分離方法の実施の形態について図面を参照して説明する。

１．畳み込み混合モデル
図１に示すように、時刻ｔにおいてＮ個の信号源１１、１２、…、１Ｎから発せられた信号源信号ｓ_ｊ（ｔ）が畳み込み混合されてｘ_ｉ（ｔ）として観測される。ｓ_ｊ（ｔ）は平均０で互いに統計的独立な非定常信号である。ｈ_ｉｊ（ｔ）は信号源１ｊからマイクロホン２ｉまでの経路の時不変なインパルス応答で、因果的で非最小位相系である。また、ｎ_ｉ（ｔ）はマイクロホン２ｉに加わる平均０、分散σ^２のガウス性白色雑音で、ｓ_ｊ（ｔ）と統計的独立である。時刻ｔにおいてＪ個のマイクロホン２１、２２、…、２Ｊで観測される観測信号ｘ_ｉ（ｔ）は式（１）で表される。ここで、Ｊ≧Ｎ≧２とする。

ここで、＊は畳み込み演算を表す。

観測信号ｘ_ｉ（ｔ）を３１で短時間フーリエ変換すると、フレーム時刻ｍにおける観測信号は式（２）により表される。式（２）において、ｗｉｎ（ｔ）は窓関数、Ｋは短時間フーリエ変換の点数、Ｔ_ｂは２つの重複窓間のシフトサイズ、ω_ｋ＝２πｋ／Ｋ、ｋ＝０，１，…，Ｋ−１、をそれぞれ表す。離散フーリエ変換の点数Ｋがｈ_ｉｊ（ｔ）のインパルス応答長より十分に大きいとき観測信号は式（３）により近似される。ここで、ｈ_ｉｊ（ｔ）のＫ点離散フーリエ変換をｈ_ｉｊ（ω_ｋ）、ｓ_ｊ（ｔ）に窓関数を乗算した後、Ｋ点短時間フーリエ変換で時間周波数領域に変換したフレーム時刻ｍの信号源信号をｓ_ｊ（ω_ｋ，ｍ）、同様に、ｎ_ｉ（ｔ）に窓関数を乗算した後、Ｋ点短時間フーリエ変換で時間周波数領域に変換したフレーム時刻ｍの雑音をｎ_ｉ（ω_ｋ，ｍ）とそれぞれ表記している。また、式（３）において、ｘ（ω_ｋ，ｍ）はフレーム時刻ｍに各マイクロホンでの観測信号ベクトル、ｓ（ω_ｋ，ｍ）はフレーム時刻ｍに各信号源信号ベクトル、混合行列Ｈ（ω_ｋ）はＮ個の信号源からＪ個のマイクロホンまでの混合行列、ｎ（ω_ｋ，ｍ）は雑音ベクトルでそれぞれ式（４）、（７）、（５）、（８）により定義される。入手可能なエポック時刻の総数をＭとすると、１≦ｍ≦Ｍとなる。信号源信号の共分散行列はＰ_ｓ（ω_ｋ，ｍ）＝Ｅ［ｓ（ω_ｋ，ｍ）ｓ（ω_ｋ，ｍ）^Ｈ］∈Ｒ^Ｎ×Ｎで、対角行列となる。Ｅ［・］と上付き添字^Ｈは期待値と複素共役転置をそれぞれ表す。また、上付き添字^ＴとＲ^Ｎ×Ｎは転置とＮ×Ｎの実数空間を表す。

信号を分離するには５１、５２、…、５Ｋで周波数ビン毎に式（９）を満足する分離行列Ｗ（ω_ｋ）を推定し、６０で信号源の割り当てを定めるパーミュテーション行列Π（ω_ｋ）∈Ｒ^Ｎ×Ｎを決定する。周波数ビン毎に独立にΠ（ω_ｋ）を決定しても信号が完全に分離する保証はなく、同一信号源から発生した信号の隣接または近接周波数ビンに相関があることを利用してパーミュテーション行列Π（ω_ｋ）を決定する。

ここで、Ｄ（ω_ｋ）∈Ｃ^Ｎ×Ｎは周波数ビン毎に異なる任意の対角行列である。

スケーリング問題とパーミュテーション問題を順に解法した後、７１、７２、…、７Ｋでｘ（ω_ｋ，ｍ）に左から分離行列Ｗ（ω_ｋ）を乗算すると、周波数ビンω_ｋにおける分離信号ｙ（ω_ｋ，ｍ）は式（１０）で表される。尚、スケーリング問題の解法については後述する。式（１０）を８０で短時間逆フーリエ変換と重複加算によって時間領域に変換すると分離信号ｙ_ｉ（ｔ）が求められる。雑音の分散σ^２が十分に小さいとき、ｙ_ｉ（ｔ）≒ｓ_ｉ（ｔ）になる。尚、分離信号ベクトルｙ（ω_ｋ，ｍ）は式（１１）により表される。

本発明に関するブラインド信号分離方法について、図１乃至図４を参照して詳細に説明する。図２乃至図４は、図１の３１における短時間フーリエ変換後、ブラインド信号分離システム４０において本発明により周波数ビン毎に推定される分離行列の算出手順を示したものである。図４は、分離行列の算出後、ブラインド信号分離システム４０において本発明によりパーミュテーション行列の算出手順を示したものである。

観測信号ｘ（ω_ｋ，ｍ）の共分散行列Ｐ_ｘ（ω_ｋ，ｍ）∈Ｃ^Ｊ×Ｊは式（１２）で与えられる。式（１３）の制約条件を課して式（１４）を満足する対角化行列Ｂ（ω_ｋ）と対角行列Λ（ω_ｋ，ｍ）を求めると、式（９）よりＢ（ω_ｋ）とＷ（ω_ｋ）の関係は式（１５）で与えられる。ただし、Ｃ^Ｊ×ＪはＪ×Ｊの複素空間を表す。

ただし、Ｉは単位行列である。

２．最小２乗型同時対角化問題とその解法
観測信号ｘ（ω_ｋ，ｍ）の共分散行列Ｐ_ｘ（ω_ｋ，ｍ）の推定値Ｐ_ｘ（ω_ｋ，ｍ）を正規化して、式（１７）を最小にする対角化行列Ｂ（ω_ｋ）と対角行列Λ（ω_ｋ，ｍ）を求める。式（１７）は最小２乗型同時対角化問題の解法として知られている。

本発明では、制約条件付き最小２乗型同時対角化問題を解法することによって混合行列Ｂ（ω_ｋ）を推定した後、分離行列Ｗ（ω_ｋ）を求める。次いで、対角行列Λ（ω_ｋ，ｍ）を推定するために、最小２乗法を用いて分離行列Ｗ（ω_ｋ）を用いた評価量を最小化する。本発明では、混合行列と対角行列の推定を交互に繰り返す。音声は低域周波数帯にフォルマントと呼ばれる振幅スペクトルのピークを有している。この音声波形の特徴を失うことなく、式（２３）によって観測信号ｘ（ω_ｋ，ｍ）の共分散行列Ｐ_ｘ（ω_ｋ，ｍ）を正規化することが、本発明の特徴の一つである。

２．１ 時間周波数領域における観測信号の共分散行列の正規化
時間領域観測信号は式（１８）の短時間フーリエ変換によって時間周波数領域に変換される。

式（１８）において、ｗ（ｔ）は窓関数、Ｋは短時間フーリエ変換の点数、Ｔ_ｓは２つの重複窓間のシフトサイズ、Ｔ_ｓはエポックサイズ、Ｎ_ｓ＋１は各エポックにおける総重複フレーム数でＫ＋Ｎ_ｓＴ_ｓ≦Ｔ_ｂ、１＝０，１，…，Ｎ_ｓをそれぞれ表す。

ステップＳ１０１においてエポック時刻ｍにおける時間周波数領域観測信号の共分散行列を式（１９）によって推定される。

混合行列Ｂ（ω_ｋ）と分離行列Ｗ（ω_ｋ）を縦続に接続したとき、そのインパルス応答は式（１５）を最小にすることによって求められる。式（２１）の評価量ζ（Ｗ（ω_ｋ））をＷ（ω_ｋ）によって微分すると、分離行列Ｗ（ω_ｋ）は式（２２）によって求められる。Ｂ（ω_ｋ）^ＨＢ（ω_ｋ）のランクがＮのときのみ、式（２１）の評価量ζ（Ｗ（ω_ｋ））は零になる。一方、Ｂ（ω_ｋ）^ＨＢ（ω_ｋ）のランクがＮ未満のとき、ζ（Ｗ（ω_ｋ））は零より大きくなる。そこで、制約条件‖ｂ_ｊ（ω_ｋ）‖_２＝１に制約条件ｒａｎｋ（Ｂ（ω_ｋ）^ＨＢ（ω_ｋ））＝Ｎを付け加え、混合行列Ｂ（ω_ｋ）の推定のための最小２乗型同時対角化問題を解法する。ここで、ｂ_ｊ（ω_ｋ）はＢ（ω_ｋ）のｊ番目の列ベクトル、‖・‖_２はユークリッドノルム、ｒａｎｋ（Ａ）は行列Ａのランクをそれぞれ表す。

ここで、‖・‖_Ｆはフロベニウスノルムを表す。

ステップＳ１０２においてＰ_ｘ（ω_ｋ，ｍ）を式（２３）によって正規化する。

２．２ 対角行列の解法
Ｐ_ｘ（ｍ）の推定値を使用して、式（２４）のバックワード型最小２乗型同時対角化問題に最小２乗法を適用すると、対角行列Λ（ｍ）はステップＳ１０４において式（２５）で推定される。式（２５）においてｄｉａｇ［Ａ］は行列Ａの対角行列を表す。式（２４）からはω_ｋの記述を省略している。

ここで、γ_ｉはラグランジェの未定乗数を表す。

２．３ 混合行列の解法
制約条件‖ｂ_ｊ‖_２＝１を課したフォワードモデル型最小２乗型同時対角化問題を周波数ビンω_ｋ毎に解くことによって、対角化行列Ｂ、即ち、混合行列を求める。評価量を式（２６）に示す。式（２６）はフォワードモデル型最小２乗型同時対角化問題として知られている。

ベクトル表現を用いると、式（２６）の評価関数は式（２７）のように表現することができる。ここで、ｒ_ｘ（ｍ）、Ｇ、ｄ（ｍ）、Ｇｄ（ｍ）はそれぞれ式（２８）〜（３１）により表される。ただし、ｖｅｃ｛Ａ｝は行列Ａの列を積み重ねることによって行列Ａを列ベ

列Λのｉ番目の要素を表す。

ステップＳ１０５においてＤを式（３２）によって計算すると、式（２７）の制約条件付きフォワードモデル型最小２乗型同時対角化問題の近似解ＧはステップＳ１０６において式（３４）のラグランジェの未定乗数法によって求められる。ステップＳ１０５からステップＳ１０８までの実施の形態については第３章で述べる。

ここで、ｖｅｃ^−１｛ａ｝はＪ^２×１の列ベクトルａをＪ×Ｊの行列に変換する演算を表す。また、ステップＳ１０６において計算されたＤの逆行列をＤ^−１と表記すると、ｄ_ｋはＤ^−１のｋ番目の列クトルである。

ステップＳ１０９において累乗法を１回用いて式（４０）を最小にするｂ_ｉを算出する。次いで、ステップＳ１１０において式（４１）のようにＢを特異値分解する。ここで、ｔｒ［Ａ］は行列Ａのトレースを表す。

ステップＳ１１１においてＢ^ＨＢのランクがＮ未満のとき、Ｂを式（４７）の行列によって置き換える。ここで、ステップＳ１１２において正規直交基底ｖ_１，ｖ_２，…，ｖ_ｒによって張られる空間に直交する空間の正規直交基底ｖ_ｒ＋１，ｖ_ｒ＋２，…，ｖ_Ｎ、同様に、ステップＳ１１３において正規直交基底ｕ_１，ｕ_２，…，ｕ_ｒによって張られる空間に直交する空間の正規直交基底ｕ_ｒ＋１，ｕ_ｒ＋２，…，ｕ_Ｎがそれぞれ求められる。

ここで、δ＞０とする。ステップＳ１１４において追加される特異値δは、式（４７）の右辺の√Ｎ／（√Ｎ＋δＮ）によって条件ｔｒ［Σ］＝√Ｎを満足するように設定される。

式（２１）を最小にする分離行列は式（２２）によって推定される。誤差の限界ε_Ｃを下回るまでＧとΛ（ｍ）の推定がステップＳ１１６において繰り返される。

３．最小２乗型同時対角化問題の次元削減法とその解法
最小２乗型同時対角化問題の解法は、対角行列と混合行列の解法から成る。混合行列の解法は対角行列の解法に比べ演算量が多いので、混合行列の推定に要する演算量の低減は、最小２乗型同時対角化問題の解法の演算量の低減に有効である。Ｐ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝はエルミート行列であるので、Ｐ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝の対角要素より上の要素は、Ｐ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝の対角要素より下の要素の複素共役に等しい。この性質を利用すると、混合行列の推定に要する演算量を低減することができる。Ｑ^ｌｏｗ∈Ｒ^{Ｊ２×Ｊ（Ｊ＋１）／２}とＱ^ｕｐ∈Ｒ^{Ｊ２×Ｊ（Ｊ＋１）／２}を式（５０）と式（５１）によってそれぞれ定義する。

３．１ ベクトルの次元削減
Ｑ^ｌｏｗはＰ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝の対角要素を含む下三角の要素をそれぞれ抽出し、次元を削減するために使用される。Ｐ_ｘ（ｍ）の対角要素を含む下三角の要素を取り出し、ベクトルａ_ｘ（ｍ）^ｌｏｗ∈Ｃ^{Ｊ（Ｊ＋１）／２×１}は式（５２）に基づきステップＳ１０３よって生成できる。尚、式（５２）と（５３）は数学的にそれぞれ式（５６）と（５７）に等しいので、行列とベクトルの乗算を実行することなく、Ｐ_ｘ（ｍ）とｂ_ｉの推定値の要素を式（５６）と（５７）に基づき並び替えることによってａ_ｘ（ｍ）^ｌｏｗ、Ｃ^ｌｏｗ∈Ｃ^{Ｊ（Ｊ＋１）／２×１}を算出することができる。

３．１ 対角化行列の解法
Ｑ^ｕｐはＰ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝の上三角の要素をそれぞれ抽出し、次元を削減するために使用される。Ｐ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝の上三角の要素を取り出し、ベクトルａ_ｘ（ｍ）^ｕｐ∈Ｃ^{Ｊ（Ｊ＋１）／２×１}とＣ^ｕｐ∈Ｃ^{Ｊ（Ｊ＋１）／２×１}は式（５８）と（５９）にステップＳ１０３においてよって生成できる。

Ｑ^ｌｏｗとＱ^ｕｐは式（６１）の関係を満たす。式（６４）にＰ_ｘ（ｍ）とｖｅｃ^−１｛ｇ_ｉ｝の対角要素を零にする行列Ａを示す。ａ_ｉは行列Ａのｉ行ｉ列の対角要素を表す。Ａのｉ行ｉ列の対角要素ａ_ｉは式（６４）によって与えられる。行列Ａを用いると、式（５８）と（５９）は式（６５）と（６６）によってそれぞれ表現することもできる。

３．２ 次元削減ベクトルを用いた評価関数
式（５２）、（５３）、（５８）、（５９）を用いると、式（２７）は式（６７）になる。ここで、γ_ｉはラグランジェの未定乗数を表す。また、Ｃ_ＤＬＳ ^ｌｏｗ（Ｃ^ｌｏｗ）とＣ_ＤＬＳ ^ｕｐ（Ｃ^ｕｐ）はそれぞれ式（６９）と（７０）で与えられる。

式（６９）には制約条件‖ｂ_ｊ‖_２＝１が課せられる。一方、この制約条件は式（７０）には課せられない。これは、式（７０）はｖｅｃ^−１｛ｇ_ｉ｝の対角より上の要素に関する最小化問題であるので、対角要素における制約条件は不要となるからである。Ｃ_ＤＬＳ ^ｌｏｗ（Ｃ^ｌｏｗ）とＣ_ＤＬＳ ^ｕｐ（Ｃ^ｕｐ）はそれぞれＣ^ｌｏｗとＣ^ｕｐに関して最小化される。

３．３ 次元削減評価関数の最小化
ステップＳ１０５においてＡ_ｄを式（７１）によって計算すると、Ｃ_ＤＬＳ ^ｌｏｗ（Ｃ^ｌｏｗ）を最小にするＣ^ｌｏｗはステップＳ１０７において式（７２）によって推定される。

式（７５）と（７６）よりＣ^ｌｏｗを用いてＣ_ＤＬＳ ^ｕｐ（Ｃ^ｕｐ）を表現すると、式（７７）を得る。Ｃ_ＤＬＳ ^ｕｐ（Ｃ^ｌｏｗ）をＣ^ｌｏｗで偏微分して式（７２）を代入すると、零ベクトルになる。即ち、Ｃ_ＤＬＳ ^ｌｏｗ（Ｃ^ｌｏｗ）を最小にする解はＣ_ＤＬＳ ^ｕｐ（Ｃ^ｌｏｗ）を最小にする複数の解の一つであるので、Ｃ_ＤＬＳ ^ｕｐ（Ｃ^ｕｐ）を最小化することなく、Ｃ_ＤＬＳ ^ｌｏｗ（Ｃ^ｌｏｗ）を最小化することによってＣ^ｌｏｗの推定値を求めることができる。

３．４ 混合行列の復元
ステップＳ１０８において式（７９）によってｃ_ｉから復元したｇ_ｉのｖｅｃ^−１｛ｇ_ｉ｝を式（８０）に示す。式（８１）と（８２）を用いると、式（８０）は式（８３）によって表現することができる。この式の導出にあたって、式（８０）の右辺第２項が式（８４）によって表現できることを使用している。

ステップＳ１０９において混合行列Ｂのｉ番目の列ベクトルは累乗法を用いて式（８５）を最小化することによって求められる。

式（５６）によりＰ_ｘ（ｍ）をａ_ｘ（ｍ）^ｌｏｗに変換した後、式（７２）によって最小化問題の解Ｃ^ｌｏｗを求め、式（８１）、（８３）、並びに（８５）によって最小解Ｃ^ｌｏｗを混合行列の推定値ｂ_ｉに変換することによって、次元削減ベクトルを利用して演算量を削減することが可能となる。即ち、混合行列の推定に次元削減ベクトルを導入することにより演算量を削減することができる。

４．パーミュテーション問題の解法
ステップＳ１１７において分離フィルタの出力信号の振幅スペクトルを求める。

ステップＳ１１８において指定した周波数帯域における振幅スペクトルの平均値を算出する。

ステップＳ１１９において指定した周波数帯域において、行列Ｂ（ω_ｋ）^ＨＢ（ω_ｋ）の条件数が最小となる周波数ビンで、振幅スペクトルがその平均値以上となっている場合、この周波数ビンをパーミュテーション行列の推定のための基準周波数ビンに選択する。

ステップＳ１２０において式（８９）に基づき、同一信号源から発生した信号の周波数ビンの電力比に相関があることを利用してパーミュテーション行列を推定する。

ここで、Ｔｒ（・）は行列のトレースを表す。また、Ｑは行列の各行に１となる要素が１箇所、その他の要素は０で、１となる要素の位置が他の行と重複しない行列の集合である。

ステップＳ１２１において全周波数ビンが終了するまで、パーミュテーション行列を推定する。式（９０）のように観測信号ｘ（ω_ｋ，ｍ）にΠ（ω_ｋ）Ｗ（ω_ｋ）を左から乗算して分離信号ｙ（ω_ｋ，ｍ）を得る。

４．１ 評価データ
図５のように４．４５×３．５５×２．５メートルの部屋に４個の信号源（スピーカ）１１（３．３５、１．３６、１．２）、１２（２．８３、２．８１、１．２）、１３（１．１４、２．２８、１．２）、１４（１．７２、０．６９、１．２）を半径１．２メートルの円の円周上に、円の中心に位置する一辺が１６．３３センチメートルの正方形の頂点に４個のマイクロホン２１（２．３４、１．７８、１．２）、２２（２．２３、１．８９、１．２）、２３（２．１１、１．７８、１．２）、２４（２．２３、１．６６、１．２）をそれぞれ配置した。符号の後のカッコは３次元の座標を表している。尚、図５は信号源（スピーカ）とマイクロホンの位置関係を示す平面図である。部屋の残響時間は１００ミリ秒から９００ミリ秒に設定し、標本化周波数８ｋＨｚ、量子化ビット数１６ビットで信号源とマイクロホンの間のインパルス応答は人工的に発生させた。実験条件は、１０００秒の音声データ、Ｋ＝８１９２点の短時間フーリエ変換、エポック当たり重複率８０％の２個のフレームの使用、窓関数にはハニング窓を用いた。マイクロホン２１、２２、２３、２４のＳＮＲの設定方法については４．２で説明する。本発明に係るブラインド信号分離方法では、ε_Ｇ＝ε_Ｃ＝１０^−６、δ（ω_ｋ）＝σ_ｒ（ω_ｋ）を用いている。スケーリング問題は周波数ビン毎に分離行列の行ベクトルを正規化することによって解法した。Ｃ言語で作成したプログラムをインテル製コアｉ７−２６００ ３．４ＧＨｚプロセッサを用いて実行した。信号源信号からマイクロホンまでの経路は時不変のインパルス応答で、因果的で非最小位相系であるので、因果的な分離行列を実現するために、Π（ω_ｋ）^−１Ｄ（ω_ｋ）^−１Ｗ（ω_ｋ）にｅ^−ｊπｋを乗算した後、逆離散フーリエ変換をして分離フィルタのインパルス応答を得た。

４．２ 評価指標
ブラインド信号分離方法の信号分離性能を次の方法で評価した．式（９１）によって観測信号における所望信号源信号と干渉信号の電力の比、式（９２）によって出力信号における所望信号源信号と干渉信号の電力の比をそれぞれ計算し、ブラインド信号分離装置の各出力の信号分離性能を求める。各出力の平均を信号分離性能とした。−γ_ｉｊ（ｔ）は式（９３）のΓ（ω_ｋ）のｉ行ｊ列の要素を、ｗ_ｉｊ（ｔ）はＷ（ω_ｋ）の要素をそれぞれ離散逆フーリエ変換したものである。また、分離行列の推定アルゴリズムにおいて収束に要した計算時間も評価指標とする。ＳＮＲは、最適な分離行列ｅ^−ｊπｋＤ（ω_ｋ）^−１（Ｈ（ω_ｋ）^ＨＨ（ω_ｋ））^−１Ｈ（ω_ｋ）^Ｈとパーミュテーション行列Π_ｏｐｔ（ω_ｋ）を使用して観測信号から信号源信号を分離した後、分離信号ｙ_ｉ（ｔ）に含まれる雑音と干渉信号の電力と所望信号源信号の電力の比によって計算した。最適なパーミュテーション行列Π_ｏｐｔ（ω_ｋ）は式（９４）によって求めた。また、非ブラインド法は、受信信号を使用して分離行列を計算した後、混合行列が入手可能であるとして、式（９５）によってパーミュテーション行列を求めた。即ち、推定した分離行列に最適なパーミュテーション行列を求めることになり、ブラインド信号分離装置の性能の上限を与えることになる。

ここで、Ｃ_ｏｐｔ（ω_ｋ）＝ｅ−^ｊπｋＤ（ω_ｋ）^−１（Ｈ（ω_ｋ）^ＨＨ（ω_ｋ））^−１Ｈ（ω_ｋ）^ＨＨ（ω_ｋ）、Ｃ（ω_ｋ）＝Ｗ（ω_ｋ）Ｈ（ω_ｋ）である。

４．３ 評価対象
最小２乗型同時対角化問題の解法を用いた３種類のフォワードモデル型ブラインド信号分離方法（非特許文献１、非特許文献２、特許文献１）を比較対象とする。従来のブラインド信号分離方法（非特許文献１、非特許文献２、特許文献１）と本発明に係るブラインド信号分離方法における分離行列の推定精度を比較するため、パーミュテーション行列の推定法は共通の手法（非特許文献３）を使用した。尚、基準周波数ビンには１１１９Ｈｚを用いた。

４．４ 評価結果
部屋の残響時間と信号分離性能の関係を図６に示す。尚、図６においてＳＮＲは２０ｄＢに設定している。また、太字の数字が最も優れた性能を表している。図６の信号分離性能から明らかなように、特許文献１と等しい信号分離性能（高い出力ＳＩＲ）を最も短時間で得ることができた。この要因はラグランジェの未定乗数法を最小２乗型同時対角化問題に導入したこと、行列のエルミート性を利用してベクトルの次元を削減したことが貢献したと考えられる。また、非ブラインド法の出力ＳＩＲ、即ち、ブラインド信号分離装置の上限に近い値を、本発明に係るブラインド信号分離装置が短い演算時間で実現できることが分かる。

５．１ 評価データ
図５のように４．４５×３．５５×２．５メートルの部屋に４個の信号源（スピーカ）１１（３．３５、１．３６、１．２）、１２（２．８３、２．８１、１．２）、１３（１．１４、２．２８、１．２）、１４（１．７２、０．６９、１．２）を半径１．２メートルの円の円周上に配置した。次いで、円の中心に位置する一辺が１６．３３センチメートルの正方形の頂点に配置した４個のマイクロホン２１（２．３４、１．７８、１．２）、２２（２．２３、１．８９、１．２）、２３（２．１１、１．７８、１．２）、２４（２．２３、１．６６、１．２）に、２つのマイクロホン２５，２６を付け加え、それぞれ座標（２．２３，１．７８，１．３２）、（２．２３，１．７８、１．０９）に配置した。部屋の残響時間は９００ミリ秒に設定し、標本化周波数８ｋＨｚ、量子化ビット数１６ビットで信号源とマイクロホンの間のインパルス応答は人工的に発生させた。実験条件は、１０００秒の音声データ、Ｋ＝８１９２点の短時間フーリエ変換、エポック当たり重複率８０％の２個のフレームの使用、窓関数にはハニング窓を用いた。ＳＮＲは２０ｄＢに設定した。本発明に係るブラインド信号分離方法では、ε_Ｇ＝ε_Ｃ＝１０^−６、δ（ω_ｋ）＝σ_ｒ（ω_ｋ）、パーミュテーション問題の解法において基準周波数ビンに決定するために、１ｋＨｚから１．５ｋＨｚを周波数帯域として用いている。スケーリング問題は周波数ビン毎に分離行列の行ベクトルを正規化することによって解法した。信号源信号からマイクロホンまでの経路は時不変のインパルス応答で、因果的で非最小位相系であるので、因果的な分離行列を実現するために、Π（ω_ｋ）^−１Ｄ（ω_ｋ）^−１Ｗ（ω_ｋ）にｅ^−ｊπｋを乗算した後、逆離散フーリエ変換をして分離フィルタのインパルス応答を得た。

６．２ 評価指標
信号源とマイクロホンの個数が共に４である場合、式（９６）に示す２４種類のパーミュテーション行列の何れか１つが各周波数ビンに割り当てられる。割り当てられたパーミュテーション行列が、任意のパーミュテーション行列に一致する割合と信号分離性能を計算する。

５．３ 評価対象
同一信号源から発生した信号の周波数ビン間の電力比に相関があることを利用したパーミュテーション問題の解法（非特許文献３）を比較対象とする。従来のブラインド信号分離方法（特許文献１）と本発明に係るブラインド信号分離方法におけるパーミュテーション行列の推定精度、反復回数、分離信号の音質を比較する。

５．４ 評価結果
図７と図８に各周波数ビンに割り当てられたパーミュテーション行列の番号を×印で示す。パーミュテーション行列Π_ｉとパーミュテーション行列の番号ｉの関係を式（９６）に示している。尚、０への割り当ては、該当するパーミュテーション行列が無いことを表している。図７では、各周波数ビンでパーミュテーション行列の番号１８に割り当てられると未知信号源への割り当てが揃うことになる。したがって、番号１８を除く他の番号への割り当ては間違いになる。一方、図８では、各周波数ビンでパーミュテーション行列の番号１５に割り当てられると未知信号源への割り当てが揃うことになる。したがって、番号１５を除く他の番号への割り当ては間違いになる。低周波数帯域（０〜２ｋＨｚ）と全周波数帯域におけるパーミュテーション行列の割り当て結果を図９にまとめる。本発明に係るパーミュテーション行列の推定法が特許文献１の方法に比べ正答率が向上していることが分かる。また、信号分離性能においても、本発明に係るパーミュテーション行列の推定法が高い出力ＳＩＲ、分離信号の音質向上（高いＰＥＳＱスコア）を達成することができた。

１１〜１Ｎ…信号源、２１〜２Ｊ…マイクロホン、３１…短時間フーリエ変換、４０…ブラインド信号分離システム、５１、５２、…、５Ｋ…最小２乗型同時対角化問題の解法、６０…パーミュテーション問題の解法、７１、７２、…、７Ｋ…畳み込み演算、８０…離散逆フーリエ変換と重複加算

Claims (4)

1. 互いに統計的に独立な未知信号源信号と未知の畳み込み混合系により、複数の未知信号源信号が混在した観測信号のみからブラインドで空間伝達関数を推定する方法であって、時間周波数領域においてエポック時刻毎に観測信号ベクトルから相互スペクトル密度行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる相互スペクトル密度行列を求め、そのノルムで全ての同エポック時刻の相互スペクトル密度行列を正規化した後、混合行列とそのエルミート行列の積、並びに相互スペクトル密度行列がエルミート行列であることを利用して、それらの行列の対角要素とその対角要素の下に位置する要素のみから混合行列を推定するために、正規化された相互スペクトル密度行列をベクトルに変換し、このベクトルの次元を削減し、混合行列の列ベクトルのノルムが１になることを制約条件としてフォワードモデル型最小２乗型同時対角化問題をラグランジュの未定乗数法により反復することなく解法することにより次元削減ベクトルを推定した後、ベクトルの次元を復元し、次元復元ベクトルを行列に変換することによって混合行列を推定し、一般化逆行列を用いて混合行列から分離行列を算出し、同一信号源から発生した信号の周波数ビンの電力比に相関があることを利用してパーミュテーション問題を解法するために、指定した周波数帯域における分離フィルタの出力信号の振幅スペクトルの平均値を求め、同周波数帯域における混合行列とそのエルミート行列の行列積によって求められる行列の条件数が最小となる周波数ビンにおいて、振幅スペクトルがその平均値以上となっている場合、この周波数ビンをパーミュテーション行列の推定のための基準周波数ビンに選択することを特徴とするブラインド信号分離方法。
2. 請求項１記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、信号源信号と観測信号間の混合行列を推定する際、混合行列の列ベクトルのノルムが１になることをフォワードモデル型最小２乗型同時対角化問題の制約条件に加え、この問題をラグランジュの未定乗数を導入して反復することなく解法することにより、ベクトルの次元削減との相乗効果によって混合行列の推定に要する演算量を削減することを特徴とするブラインド信号分離方法。
3. 請求項１乃至請求項２のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、パーミュテーション問題の解法の際、基準周波数ビンの決定に、推定された混合行列とそのエルミート行列の行列積によって求められる行列の条件数、並びに、分離フィルタの出力信号の振幅スペクトルを使用することを特徴とするブラインド信号分離方法。
4. 請求項１乃至請求項３のいずれか１項に記載のブラインド信号分離方法を用いて信号源分離を行うように構成されていることを特徴とするブラインド信号分離方法を用いたブラインド信号分離装置。

Images