特開 2022-120750 到来方向別信号波形分離方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2022120750

(43)【公開日】2022-08-18

(54)【発明の名称】到来方向別信号波形分離方法

(51)【国際特許分類】

   G10K 15/00 20060101AFI20220810BHJP

【ＦＩ】

G10K15/00 L

【審査請求】未請求

【請求項の数】8

【出願形態】書面

(21)【出願番号】P 2021056192

(22)【出願日】2021-02-05

(71)【出願人】

【識別番号】518299688

【氏名又は名称】松倉 敏寛

(72)【発明者】

【氏名】松倉 敏寛

(57)【要約】

【課題】到来方向別に信号波形を分離することができる。
【解決手段】源信号空間事前情報が既知の時に、源信号空間事前情報より推定される混合行列が分離行列の逆行列であるという情報を含む制約条件を有する目的関数を用いて、分離行列を行列ブロックごとに更新した後、得られた分離行列と観測信号から源信号を推定することを特徴とし、源信号空間事前情報が未知の時には、観測信号と源信号の統計的分散に関する事前情報を用いて源信号を推定した後、信号基底ベクトルを特徴量として用い、混合モデルとしてフォンミーゼス・フィッシャー分布を仮定して帰属度及びモデルパラメータを交互に更新することで、推定した源信号の交換不定性問題を解決することを特徴とする。

【特許請求の範囲】

【請求項1】

空気や物体の振動を源信号と定義し、多方向から到来する前記源信号が混ざり合った状態でマイクやセンサで取得される信号を観測信号と定義する場合、
前記観測信号をその到来方向別の前記源信号波形に分離すると共に、前記源信号の到来方向を推定する方向別信号波形分離方法であって、
前記観測信号を到来方向ごとの前記源信号へと分離する行列を分離行列と定義する場合において、前記分離行列を前記観測信号に掛けて前記源信号の推定値である分離信号を算定することを特徴とする分離信号算定部と、前記源信号の到来方向又は波形形状の特徴量をフォンミーゼス・フィッシャー分布に基づくクラスタ分析手法によりカテゴリ別に分類するクラスタリング部及び信号分離結果を基に前記源信号の到来方向を推定する源信号到来方向算定部を演算機能として有することを特徴とし、
複数の前記源信号の内、いくつかの空間事前情報（到来方向事前情報）が得られている場合においては、前記分離信号算定部で、前記観測信号及び前記空間事前情報を利用して前記源信号を到来方向別に分離する空間事前情報に基づく分離信号算定の演算を実行して、分離信号を算定し、
前記源信号の空間事前情報が得られていない場合においては、前記分離信号算定部で前記観測信号と前記源信号の分散に関する事前情報を基にして前記分離信号を算定した後に、前記クラスタリング部及び前記源信号到来方向算定部の演算を実施することを特徴とする、
信号波形分離方法。

【請求項2】

請求項１に記載の前記空間事前情報に基づく分離信号算定のための演算方法であって、
前記分離信号算定方法は、前記源信号到来方向を事前情報として用いて、時間フレームτ、周波数ビンｆごとの分離信号ｙ（ｆ，τ）を算定する方法であって、
センサ配置位置と前記源信号到来方向から、センサから前記源信号までの周波数伝達関数を算定し、
前記周波数伝達関数をセンサ番号が行ベクトル、前記源信号番号が列ベクトルとなるように配置した行列を混合行列Ａ（ｆ）と定義する場合、
前記混合行列とセンサでの前記観測信号ｘ（ｆ，τ）が既知という条件下での前記分離信号の確率密度（事後確率）を最大化することを目標として前記分離信号を算定することを特徴とし、
目的信号以外の雑音と前記分離信号の確率密度がガウス分布に従うと仮定することで前記雑音の分散σ_Ｅ（ｆ）と前記分離信号の分散σ_ｎ（ｆ）との比率により算定される正則化項Σ（ｆ）を前記混合行列の２乗に加えた行列Ａ^Ｈ（ｆ）Ａ（ｆ）＋Σ（ｆ）の逆行列を算定し、
前記逆行列から式（９２－２）により前記分離行列Ｗ（ｆ）を算定し、式（９２－１）により前記分離行列を前記観測信号に掛けることで、前記分離信号を算定することを特徴とする、
前記分離信号算定方法。

【数92】

【請求項3】

請求項１に記載の前記空間事前情報に基づく分離信号算定のための演算方法であって、
前記分離行列は、前記源信号番号を行ベクトル、センサ番号を列ベクトルに有しており、
センサ配置位置と前記源信号到来方向から、センサから前記源信号までの周波数伝達関数を算定し、
前記周波数伝達関数をセンサ番号が行ベクトル、前記源信号番号が列ベクトルとなるように配置した行列を混合行列と定義する場合、
前記混合行列が前記分離行列の逆行列であるという情報を制約条件として設定することを特徴とし、
前記制約条件は、互いに同一の前記源信号番号を有する前記分離行列の行ベクトルと前記混合行列の列ベクトルとの内積が１になり、互いに異なる前記源信号番号を有する前記分離行列の行ベクトルと前記混合行列の列ベクトルの内積が０になるという情報を含むことを特徴とし、
前記制約条件を含む目的関数を最適化することで、前記分離行列を算定し、
前記分離行列を前記観測信号に掛けて前記分離信号を算定することを特徴とする、
前記分離信号を算定するための演算方法。

【請求項4】

請求項１に記載の前記空間事前情報に基づく分離信号算定のための演算方法であって、
正確な前記源信号波形を求めるための目的関数を最適化する前記分離行列を算定する過程において、
前記目的関数を前記分離行列で偏微分して勾配降下方向を算定し、
前記源信号到来方向の事前情報より概略的に算定される前記分離行列として、前記混合行列の逆行列又は前記混合行列をエルミート転置した行列をテンソル成分に持つ行列又は請求項２に記載の式（９２－２）の行列を前記分離行列の初期値として用いて、
前記勾配降下方向に対して前記分離行列を繰り返し更新することで前記目的関数を最適化する前記分離行列を算定し、
前記分離行列を前記観測信号に掛けて前記分離信号を算定することを特徴とする、
前記分離信号を算定するための演算方法。

【請求項5】

請求項３に記載の前記制約条件を含む目的関数を用いて、最適な前記分離行列を探索する方法であって、
正確な前記源信号波形を求めるための目的関数に基づいて前記分離行列を更新する際に、前記目的関数が複雑な解空間を有する場合において、局所解に陥ることなく安定的に前記分離行列の更新を行うための方法であって、
前記分離行列は、前記源信号番号を行ベクトル、センサ番号を列ベクトルに有しており、
前記分離行列を行ベクトルごとに分割した時のテンソルの集合又は任意の範囲で分割した時のテンソルの集合を行列ブロックと定義する場合、前記行列ブロックごとに更新量を調整して前記分離行列を更新することを特徴とし、
前記行列ブロックごとのテンソル成分の２乗和をその行列のブロック単位２乗和と定義し、
前記目的関数を前記分離行列で偏微分して求めた勾配を損失関数と定義し、
前記目的関数の勾配降下方向への更新量を調整するために、前記損失関数に掛ける係数を学習率と定義する場合、
前記損失関数の前記ブロック単位２乗和を求めて、１つ前のタイムステップまでの前記損失関数の前記ブロック単位２乗和の累積値に加算することで、現在のタイムステップにおける前記損失関数の前記ブロック単位２乗和の累積値を算定した後、
現在のタイムステップにおける前記損失関数の前記ブロック単位２乗和の累積値と、１つ前のタイムステップにおける前記分離行列の更新量の前記ブロック単位２乗和の累積値を用いて、前記行列ブロックごとに、前記学習率を算定し、
１つ前のタイムステップにおける前記分離行列の更新量に定数を掛けた値に対して、前記学習率を前記損失関数に掛けた値を加算することで、現在のタイムステップにおける前記分離行列の更新量を前記行列ブロックごとに算定した後、
前記分離行列の更新量の前記ブロック単位２乗和を求めて、１つ前のタイムステップにおける前記分離行列の更新量の前記ブロック単位２乗和の累積値に加算することで、現在のタイムステップにおける前記分離行列の更新量の前記ブロック単位２乗和の累積値を算定するという４段階の変数更新ステップに基づいて、
前記分離行列を更新する演算を何度も繰り返すことで、前記目的関数を最適化する前記分離行列を探索することを特徴とする、
前記分離行列の探索方法。

【請求項6】

請求項１に記載のクラスタリング部における演算方法は、
複数の確率密度関数の線形和で表現される混合モデルにおける未知のパラメータの集まりをクラスタと定義する場合、
周波数ビン・前記源信号ごと又は時間フレーム・周波数ビン・前記源信号ごとに、前記源信号の到来方向又は波形形状の特徴量として信号基底ベクトルを算出し、
前記混合モデルの確率密度関数としてフォンミーゼス・フィッシャー分布を使用し、
前記混合モデルにおける前記信号基底ベクトルの前記クラスタへの帰属確率を帰属度と定義し、
前記クラスタごとの重要度を混合重みと定義し、前記信号基底ベクトルの前記クラスタごとの平均値を平均方向と定義し、前記平均方向周りの前記信号基底ベクトルのばらつき度合を方向集中度と定義する場合、
前記フォンミーゼス・フィッシャー分布の確率密度関数を最大化することを目標として、
前記混合重みと前記平均方向と前記方向集中度を固定して前記帰属度を更新する演算と前記帰属度を固定して前記混合重みと前記平均方向を更新すると共に、近似式又は勾配法又は最尤推定法を用いて、前記方向集中度を更新する演算を交互に繰り返すことで、最適な前記帰属度及び前記混合重みと前記平均方向と前記方向集中度を算定することを特徴とする、
クラスタリング演算方法。

【請求項7】

前記帰属度の値を０又は１の離散値に限定した際に用いることのできる請求項１に記載のクラスタリング部における演算方法であって、
前記フォンミーゼス・フィッシャー分布の確率密度関数を最大化することを目標として、
同一の周波数ビン又は時間フレーム・周波数ビンに属する前記信号基底ベクトルが同一のクラスタには所属しないという制約条件を付与しつつ、前記混合重みと前記平均方向と前記方向集中度を固定して前記帰属度を０又は１のいずれかの値に更新する演算と、
前記帰属度を固定して前記混合重みと前記平均方向を更新すると共に近似式又は勾配法又は最尤推定法を用いて、前記方向集中度を更新する演算とを交互に繰り返すことで、最適な前記帰属度及び前記混合重みと前記平均方向と前記方向集中度を算定することを特徴とする、
クラスタリング演算方法。

【請求項8】

請求項６又は請求項７に記載のクラスタリング演算において、前記方向集中度κ_ｃ又はκを算定するための演算方法であって、
前記方向集中度は、前記フォンミーゼス・フィッシャー分布の式中に含まれる第１種変形ベッセル関数Ｉ_０．５Ｍ（κ_ｃ）の中にも含まれており、
前記信号基底ベクトルｖ_ｎ（ｆ）と前記信号基底ベクトルに対して、その重要度又は信頼度に基づいて設定した基底重みα_ｎ（ｆ）及び前記帰属度ｕ_ｃｎ（ｆ）より、階数の１つ異なる前記第１種変形ベッセル関数の比率の概略値である概略比率を式（９３）又は式（９４）を用いて算定し、

【数93】

【数94】

階数の１つ異なる前記第１種変形ベッセル関数の比率を連分数に展開することで得られる式（９５）又は式（９６）に示す近似式を用いて、前記概略比率より、前記方向集中度の近似値を算定し、

【数95】

【数96】

前記方向集中度の近似値を階数の１つ異なる前記第１種変形ベッセル関数にそれぞれ代入して、階数の１つ異なる前記第１種変形ベッセル関数の比率の詳細値である詳細比率を式（９７）により算定し、

【数97】

前記詳細比率と前記概略比率の誤差を最小化することを目標として、前記詳細比率を前記方向集中度の近似値で微分した値を勾配としたニュートン法の更新式（式（９８））を適用し、前記方向集中度を更新し、

【数98】

前回の更新ステップにおいて更新した前記方向集中度を階数の１つ異なる前記第１種変形ベッセル関数にそれぞれ代入して、前記詳細比率を式（９７）により更新し、更新した前記詳細比率と更新した前記方向集中度を前記ニュートン法の更新式（式（９８））に代入し、前記方向集中度を更新するという式（９７）と式（９８）を交互に用いた演算過程を何度か繰り返すことで、
前記方向集中度を算定することを特徴とする、
クラスタリング演算における前記方向集中度算定方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、到来方向別に信号波形を分離する方法に係り、特に、源信号到来方向事前情報（源信号空間事前情報）に基づく分離信号算定方法及びフォンミーゼス・フィッシャー分布を混合モデルとして用いるクラスタリング方法に関する。

【背景技術】

【発明が解決しようとする課題】

【0002】

一般の道路環境において、車両のエンジン音やタイヤ転動音、道路橋においては、橋梁に車両が進入する衝撃により橋梁構造から発生するジョイント音や車両の走行により生じる桁のたわみ振動に伴う低周波音が問題となっている（非特許文献１）。特に橋梁構造から発生する低周波音は健康被害をもたらすとの報告があり（非特許文献２）、発生源を特定して対策を講じることが重要となる。又、住宅地においては、物の落下やファンの稼働音等による発生源を特定することが困難な不思議音が問題になっている。騒音や不思議音の発生源を調査するためには、数多くの観測データの中から問題となっている音を特定する必要があり、労力がかかる。この他にも工事騒音の予測等に役立てるために、スピーカから建設機械音を発生させ、複数のマイクロフォンアレイを用いて、音の到来方向を特定することで、機械別の騒音予測を行っている事例もあるが、方向特定のみに留まっている。方向別に音の周波数特性を算定して騒音の発生源を調査する時や、複数の環境音の中から問題となっている音を特定する時などにおいて、労力削減のために、人工知能を用いる場合、信号波形の特徴量が必要となる。この時、源信号の到来方向別の時間領域又は周波数領域での波形が得られれば、より高精度に対象とする音の特徴量を抽出することが可能になり、高精度な判定に繋がると考えられる。

【0003】

異常検知の分野では、観測された音の時間波形からパワースペクトルや波形形状の特徴量を抽出し、人工知能に入力することで、異常性を判定する手法が提案されている（非特許文献３）。この異常音検知手法は、機械の故障予測を初めとして、銃声や悲鳴，爆発音の検知による犯罪検出や道路における車の衝突やスリップ音の検知による事故検出など、様々な目的で用いられている（非特許文献４、５、６）。もし事前に方向別に波形を分離することができれば、人工知能に目的とする波形の特徴量とその方向頻度をより高精度に学習させることや学習データに基づいて、異常性を判定させることが可能になると考えられる。

【0004】

従来の騒音や異音の発生源調査には、ミュージック法（非特許文献７）等のサブスペース法がよく用いられてきた。例えばミュージック法は、複数の観測点で得られた信号の相関行列とセンサ配置から定まるステアリングベクトルの情報を基に方向別のＭＵＳＩＣスペクトルを求める手法である。しかし、ＭＵＳＩＣスペクトルは方向別の相対的な音の大きさの比率を表しているため、各方向における到来波の強度比率を評価することはできるが、実際にセンサで観測された音の大きさや時間波形を求めることはできない。ミュージック法を用いて、複数の信号波形が混合した観測信号から分離波形の時系列データを得るには、源信号到来方向の情報から混合行列を推定し、その逆行列を分離行列として用いる必要がある。しかし、この方法では、源信号が互いに近接して存在する場合に混合行列が悪条件になるため、分離信号の算定精度が悪くなる。広帯域信号の分離波形を得る時にも、低周波数側において混合行列は悪条件になる。そのため、ミュージック法により推定した混合行列の逆行列を算定する方法では、全ての源信号が十分離れている状況において、解析対象周波数を狭帯域に限定した状態でしか、分離信号を得ることができない。

【0005】

音楽信号や音声信号の分離においては、複数の混じり合った音の中から、目的音を抽出する手法や混合音信号を因子別に分離する手法が発展してきた。代表的な手法としては、独立成分分析、独立低ランク行列分析がある。これらの源信号分離手法は、源信号の統計的独立性を仮定して、信号分離を行う方法である。独立成分分析や独立低ランク行列分析を周波数領域で用いることにより分離行列を算定した場合でも、広帯域信号における低周波数側での分離信号の推定精度劣化は起きるが、その劣化の程度が、混合行列の逆行列を用いた時よりも少なくなる。そのため、源信号が近接して存在する場合や広帯域信号から分離信号を算定する場合においても適用可能となる。

【0006】

独立成分分析では、源信号の確率分布形状を想定して分離信号を求めている。しかし、確率分布を厳密に定めることはできないため、信号分離結果には常に一定レベルの誤差が含まれる。又、源信号数が観測点数よりも多い場合、正しく信号分離できないという性質があり、目的音以外の複数のノイズが混入した場合、信号分離精度が悪くなるという課題がある。

【0007】

独立成分分析や独立低ランク行列分析には、源信号の推定値が現れる順序を定めることができないというパーティミション問題があり、周波数領域で信号分離する時にパーティミションが起こると、分離信号に大きな誤差が生じるという課題がある。信号分離結果に誤差が含まれる場合、分離結果から得られる方向特徴量も誤差の影響を受けるため、その方向特徴量を用いて、パーティミション問題を解決することが困難になる。

【0008】

源信号到来方向は、分離行列の逆行列である混合行列のテンソル成分を用いて、源信号位置からセンサ位置までの周波数伝達関数を算定し、周波数伝達関数から源信号のセンサ間位相差を推定することで算定する。そのため、分離行列の算定結果に誤差が含まれると、その影響を受けて源信号到来方向推定結果の精度が劣化する

【0009】

このように分離行列の算定誤差は、周波数ビンごとの分離信号の誤差のみにとどまらず、信号基底ベクトルのクラスタ分析の誤差や源信号到来方向の算定誤差につながる。

【0010】

独立成分分析を用いる場合、確率分布形状を厳密に定めることができれば、上記の課題は全て解決するが、未知の源信号の確率分布形状を厳密に定めることは困難であり、そのため独立成分分析による信号分離結果には一定レベルの誤差が含まれる結果となる。

【0011】

パーティミションの解決手法としてこれまで分離行列の指向性パターンを手掛かりにする方法（非特許文献８，９）や分離行列の逆行列の周波数間での連続性を利用する方法（非特許文献１０）が提案されている。しかし、これらの方法では、周波数ごとに信号基底ベクトルが異なるため、広帯域信号における分離が困難になる。この問題を解決するために、分離行列の逆行列を源信号と観測点間の周波数伝達関数で近似することで、全周波数で一貫した源信号ラベルを構成する手法が提案されている（非特許文献１１）。以上の方法では、観測データにノイズが混入するなどの要因で、分離行列（分離信号）に算定誤差が生じた時に、クラスタごとのに信号基底ベクトルもその影響を受けて算定誤差が生じるため、最適解にたどり着けないという問題がある。又、クラスタ算定時に信号基底ベクトルのクラスタごとの平均方向の初期値を定める必要があるが、従来のクラスタ算定方法では、算定結果が初期値に依存し、局所解に陥ることがあるという問題がある。この初期値依存性は、分離行列（分離信号）に算定誤差が生じた時に顕著に現れる。

【0012】

これまで独立成分分析は信号分離に活用されてきたが、信号分離後に高精度に源信号到来方向を推定する手法は確率されていない。その原因として分離行列に基づく信号の到来方向推定方法では、分離行列の算定結果に誤差が含まれる時にその影響を受けて源信号到来方向推定結果の精度が劣化することが考えられる。そのため、これまでは源信号到来方向にはミュージック法等のサブスペース法、信号分離には独立成分分析が用いられることが多かった。

【0013】

ミュージック法等のサブスペース法とカメラ映像を併せて用いた場合、源信号到来方向を高精度に推定できる。このような方法で事前に源信号の空間情報が得られている場合、この空間事前情報を積極的に用いると、センサ－源信号間の周波数伝達関数が求まる。この周波数伝達関数を用いて混合行列を求め、その逆行列により分離行列を算定できる。この分離行列を用いれば、理論上、空間事前情報を信号分離に利用すると、高精度に信号分離できることになる。しかし、源信号の解析対象周波数が広帯域にわたる場合や源信号位置が互いに近接している場合に、混合行列が悪条件となり、分離行列の算定精度が悪くなる。この課題を解決し、信号波形の空間事前情報を用いて、高精度に信号分離する方法は、現段階で確立されていない。

【先行技術文献】

【非特許文献】

【0014】

【非特許文献1】木村真也，小野和行，幸寺駿，金哲佑，川谷充郎：道路橋における交通振動に伴う低周波音伝播特性に関する研究，構造工学論文集，Ｖｏｌ．６４Ａ，ｐｐ３０７－３１４，２０１８．

【非特許文献2】中村俊一，時田保夫，織田厚：低周波音に対する感覚と評価に関する基礎的研究，昭和５５年度文部省科 学研究費「環境科学」特別研究，１９７９．

【非特許文献3】Ａ．Ｉｔｏ，Ａ．Ａｉｂａ，Ｍ．Ｉｔｏ ａｎｄ Ｓ．Ｍａｋｉｎｏ，“Ｄｅｔｅｃｔｉｏｎ ｏｆ ａｂｎｏｒｍａｌ ｓｏｕｎｄ ｕｓｉｎｇ ｍｕｌｔｉ－ｓｔａｇｅ ＧＭＭ ｆｏｒ ｓｕｒｖｅｉｌｌａｎｃｅ ｍｉｃｒｏｐｈｏｎｅ，”Ｐｒｏｃ．Ｉｎｔ．Ｃｏｎｆ．Ｉｎｆｏｒｍａｔｉｏｎ Ａｓｓｕｒａｎｃｅ ａｎｄ Ｓｅｃｕｒｉｔｙ，Ｖｏｌ．１，ｐｐ．７３３－７３６（２００９）．

【非特許文献4】Ｃ．Ｎ．Ｄｏｕｋａｓ ａｎｄ Ｉ．Ｍａｇｌｏｇｉａｎｎｉｓ，“Ｅｍｅｒｇｅｎｃｙ ｆａｌｌ ｉｎｃｉｄｅｎｔｓ ｄｅｔｅｃｔｉｏｎ ｉｎ ａｓｓｉｓｔｅｄ ｌｉｖｉｎｇ ｅｎｖｉｒｏｎｍｅｎｔｓ ｕｔｉｌｉｚｉｎｇ ｍｏｔｉｏｎ，ｓｏｕｎｄ，ａｎｄ ｖｉｓｕａｌ ｐｅｒｃｅｐｔｕａｌ ｃｏｍｐｏｎｅｎｔｓ，”ＩＥＥＥ Ｔｒａｎｓ．Ｉｎｆ．Ｔｅｃｈｎｏｌ．Ｂｉｏｍｅｄ．，１５，２７７－２８９（２０１１）．

【非特許文献5】Ｓ．Ｌｅｃｏｍｔｅ，Ｒ．Ｌｅｎｇｅｌｌ´ｅ，Ｃ．Ｒｉｃｈａｒｄ，Ｆ．Ｃａｐｍａｎ ａｎｄ Ｂ．Ｒａｖｅｒａ，“Ａｂｎｏｒｍａｌ ｅｖｅｎｔｓ ｄｅｔｅｃｔｉｏｎ ｕｓｉｎｇ ｕｎｓｕｐｅｒｖｉｓｅｄ ｏｎｅ－ｃｌａｓｓ ＳＶＭ：Ａｐｐｌｉｃａｔｉｏｎ ｔｏ ａｕｄｉｏ ｓｕｒｖｅｉｌｌａｎｃｅ ａｎｄ ｅｖａｌｕａｔｉｏｎ，”Ｐｒｏｃ．８ｔｈ ＩＥＥＥ Ｉｎｔ．Ｃｏｎｆ．Ａｄｖａｎｃｅｄ Ｖｉｄｅｏ ａｎｄ Ｓｉｇｎａｌ Ｂａｓｅｄ Ｓｕｒｖｅｉｌｌａｎｃｅ，ｐｐ．１２４－１２９（２０１１）．

【非特許文献6】Ｅ．Ｍａｒｃｈｉ，Ｆ．Ｖｅｓｐｅｒｉｎｉ，Ｆ．Ｅｙｂｅｎ，Ｓ．Ｓｑｕａｒｔｉｎｉ ａｎｄ Ｂ．Ｓｃｈｕｌｌｅｒ，“Ａ ｎｏｖｅｌ ａｐｐｒｏａｃｈ ｆｏｒ ａｕｔｏｍａｔｉｃ ａｃｏｕｓｔｉｃ ｎｏｖｅｌｔｙ ｄｅｔｅｃｔｉｏｎ ｕｓｉｎｇ ａ ｄｅｎｏｉｓｉｎｇ ａｕｔｏｅｎｃｏｄｅｒ ｗｉｔｈ ｂｉｄｉｒｅｃｔｉｏｎａｌ ＬＳＴＭ ｎｅｕｒａｌ ｎｅｔｗｏｒｋｓ，”Ｐｒｏｃ．ＩＣＡＳＳＰ ２０１５，ｐｐ．１９９６－２０００（２０１５）．

【非特許文献7】Ｐ．Ｓｔｏｉｃａ ａｎｄ Ａ．Ｎｅｈｏｒａｉ，“ＭＵＳＩＣ，ｍａｘｉｍｕｍ ｌｉｋｅｌｉｈｏｏｄ，ａｎｄ Ｃｒａｍｅｒ－Ｒａｏ ｂｏｕｎｄ，”ＩＥＥＥ Ｔｒａｎｓ．Ａｃｏｕｓｔ．Ｓｐｅｅｃｈ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓ．ｖｏｌ．３７，ｎｏ．５，ｐｐ．７２０－７４１，Ｍａｙ１９８９．

【非特許文献8】Ｈ．Ｓａｒｕｗａｔａｒｉ，Ｓ．Ｋｕｒｉｔａ，Ｋ．Ｔａｋｅｄａ，Ｆ．Ｉｔａｋｕｒａ，Ｔ．Ｎｉｓｈｉｋａｗａ，ＫｉｙｏｈｉｒｏＳｈｉｋａｎｏ：Ｂｌｉｎｄ Ｓｏｕｒｃｅ Ｓｅｐａｒａｔｉｏｎ Ｃｏｍｂｉｎｉｎｇ Ｉｎｄｅｐｅｎｄｅｎｔ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ ａｎｄ Ｂｅａｍｆｏｒｍｉｎｇ，ＥＵＲＡＳＩＰ Ｊｏｕｒｎａｌ ｏｎ Ａｐｐｌｉｅｄ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ，１１３５－１１４６，Ｎｏｖ．Ｊａｎ．２００３．

【非特許文献9】Ｍ．Ｚ．Ｉｋｒａｍ ａｎｄ Ｄ．Ｒ．Ｍｏｒｇａｎ：Ｐｅｒｍｕｔａｔｉｏｎ ｉｎｃｏｎｓｉｓｔｅｎｃｙ ｉｎ ｂｌｉｎｄ ｓｐｅｅｃｈ ｓｅｐａｒａｔｉｏｎ：ｉｎｖｅｓｔｉｇａｔｉｏｎ ａｎｄ ｓｏｌｕｔｉｏｎｓ，ＩＥＥＥ Ｔｒａｎｓ．Ｓｐｅｅｃｈ ａｎｄ Ａｕｄｉｏ Ｐｒｏｃｅｓｓｉｎｇ，Ｖｏｌ．１３，ｎｏ．１，ｐｐ．１～１３，２００５．

【非特許文献10】Ｆ．Ａｓａｎｏ，Ｓ．Ｉｋｅｄａ，Ｍ．Ｏｇａｗａ，Ｈ．Ａｓｏｈ，ａｎｄ Ｎ．Ｋｉｔａｗａｋｉ：Ｃｏｍｂｉｎｅｄ ａｐｐｒｏａｃｈ ｏｆ ａｒｒａｙ ｐｒｏｃｅｓｓｉｎｇ ａｎｄ ｉｎｄｅｐｅｎｄｅｎｔ ｃｏｍｐｏｎｅｎｔ ａｎａｌｙｓｉｓ ｆｏｒ ｂｌｉｎｄ ｓｅｐａｒａｔｉｏｎ ｏｆ ａｃｏｕｓｔｉｃ ｓｉｇｎａｌｓ，ＩＥＥＥ Ｔｒａｎｓ．Ｓｐｅｅｃｈ ａｎｄ Ａｕｄｉｏ Ｐｒｏｃｅｓｓｉｎｇ，Ｖｏｌ．１１，ｎｏ．３，ｐｐ．２０４～２１５，Ｍａｙ．２００３．

【非特許文献11】Ｈ．Ｓａｗａｄａ，Ｓ．Ａｒａｋｉ，Ｒ．Ｍｕｋａｉ，Ｓ．Ｍａｋｉｎｏ：Ｂｌｉｎｄ ｅｘｔｒａｃｔｉｏｎ ｏｆ ａ ｄｏｍｉｎａｎｔ ｓｏｕｒｃｅ ｓｉｇｎａｌ ｆｒｏｍ ｍｉｘｔｕｒｅｓ ｏｆ ｍａｎｙ ｓｏｕｒｃｅｓ，”Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｃｏｎｆｅｒｅｎｃｅ ｏｎ（Ａｃｏｕｓｔ Ｓｐｅｅｃｈ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓ），Ｖｏｌ．３，６１－６４（２００５）．

【非特許文献12】澤田 宏，向井 良，荒木 章子，牧野 昭二：多音源に対する周波数領域ブラインド音源分離，社団法人人工知能学会，ＪＳＡＩ Ｔｅｃｈｎｉｃａｌ Ｒｅｐｏｒｔ，ＳＩＧ－ＣＨａｌｌｅｇｅ－０５２２－３（１０／１４），２００５．

【非特許文献13】Ｈ．Ｓａｗａｄａ，Ｈ．Ｋａｍｅｏｋａ，Ｓ．Ａｒａｋｉ，ａｎｄ Ｎ．Ｕｅｄａ，“Ｍｕｌｔｉｃｈａｎｎｅｌ ｅｘｔｅｎｓｉｏｎｓ ｏｆ ｎｏｎ－ｎｅｇａｔｉｖｅ ｍａｔｒｉｘ ｆａｃｔｏｒｉｚａｔｉｏｎ ｗｉｔｈ ｃｏｍｐｌｅｘ－ｖａｌｕｅｄ ｄａｔａ，”ＩＥＥＥ Ｔｒａｎｓ．ＡＳＬＰ，ｖｏｌ．２１，Ｎｏ．５，ｐｐ．９７１－９８２，２０１３．ＭＮＭＦ１

【非特許文献14】北村 大地，小野 順貴，澤田 宏，亀岡 弘和，猿渡 洋，″独立低ランク行列分析に基づくブラインド音源分離，″ＩＥＩＣＥ Ｔｅｃｈｎｉｃａｌ Ｒｅｐｏｒｔ，ＥＡ２０１７－５６，ｖｏｌ．１１７，ｎｏ．２５５，ｐｐ．７３－８０，Ｏｃｔ．２０１７．

【非特許文献15】Ｎ．Ｑ．Ｋ．Ｄｕｏｎｇ，Ｅ．Ｖｉｎｃｅｎｔ，ａｎｄ Ｒ．Ｇｒｉｂｏｎｖａｌ，“Ｕｎｄｅｒｄｅｔｅｒｍｉｎｅｄ ｒｅｖｅｒｂｅｒａｎｔ ａｕｄｉｏ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ ｕｓｉｎｇ ａ ｆｕｌｌｒａｎｋ ｓｐａｔｉａｌ ｃｏｖａｒｉａｎｃｅ ｍｏｄｅｌ，”ＩＥＥＥ Ｔｒａｎｓ．ＡＳＬＰ，ｖｏｌ．１８，ｎｏ．７，ｐｐ．１８３０－１８４０，２０１０．ＭＮＭＦ２

【非特許文献16】Ｃ．Ｆ▲ｅ▼ｖｏｔｔｅ，Ｎ．Ｂｅｒｔｉｎ，ａｎｄ Ｊ．－Ｌ．Ｄｕｒｒｉｅｕ，“Ｎｏｎｎｅｇａｔｉｖｅ ｍａｔｒｉｘ ｆａｃｔｏｒｉｚａｔｉｏｎ ｗｉｔｈ ｔｈｅ Ｉｔａｋｕｒａ－Ｓａｉｔｏ ｄｉｖｅｒｇｅｎｃｅ．Ｗｉｔｈ ａｐｐｌｉｃａｔｉｏｎ ｔｏ ｍｕｓｉｃ ａｎａｌｙｓｉｓ，”Ｎｅｕｒａｌ Ｃｏｍｐｕｔａｔｉｏｎ，ｖｏｌ．２１，ｎｏ．３，ｐｐ．７９３－８３０，２００９．

【非特許文献17】Ｄ．Ｋｉｔａｍｕｒａ，Ｎ．Ｏｎｏ，Ｈ．Ｓａｗａｄａ，Ｈ．Ｋａｍｅｏｋａ，ａｎｄ Ｈ．Ｓａｒｕｗａｔａｒｉ，“Ｄｅｔｅｒｍｉｎｅｄ ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ ｕｎｉｆｙｉｎｇ ｉｎｄｅｐｅｎｄｅｎｔ ｖｅｃｔｏｒ ａｎａｌｙｓｉｓ ａｎｄ ｎｏｎｎｅｇａｔｉｖｅ ｍａｔｒｉｘ ｆａｃｔｏｒｉｚａｔｉｏｎ，”ＩＥＥＥ／ＡＣＭ Ｔｒａｎｓ．ＡＳＬＰ，ｖｏｌ．２４，ｎｏ．９，ｐｐ．１６２６－１６４１，２０１６．

【非特許文献18】Ｎ．Ｏｎｏ，“Ｓｔａｂｌｅ ａｎｄ ｆａｓｔ ｕｐｄａｔｅ ｒｕｌｅｓ ｆｏｒ ｉｎｄｅｐｅｎｄｅｎｔ ｖｅｃｔｏｒ ａｎａｌｙｓｉｓ ｂａｓｅｄ ｏｎ ａｕｘｉｌｉａｒｙ ｆｕｎｃｔｉｏｎ ｔｅｃｈｎｉｑｕｅ，”Ｐｒｏｃ．ＷＡＳＰＡＡ，ｐｐ．１８９－１９２，２０１１．

【非特許文献19】北村 大地，角野 隼斗，高宗 典玄，高道 慎之介，猿渡 洋，小野 順貴，″独立深層学習行列分析に基づく多チャネル音源分離の実験的評価，″ＩＥＩＣＥ Ｔｅｃｈｎｉｃａｌ Ｒｅｐｏｒ ｔ，ＥＡ２０１７－１０４，ｖｏｌ．１１７，ｎｏ．５１５，ｐｐ．１３－２０，Ｍ ａｒ．２０１８．

【非特許文献20】椋本博学，林和則，金子めぐみ，“Ｋｈａｔｒｉ－Ｒａｏ積拡張アレーを用いた圧縮センシングによる到来角推定法，”信学技報，ＲＣＣ２０１４－１１，ＭＩＣＴ２０１４－１１，ｐｐ．５１－５６，Ｍａｙ ２０１４．

【非特許文献21】Ｊ．Ｃｈｅｎ ａｎｄ Ｘ．Ｈｕｏ，“Ｔｈｅｏｒｅｔｉｃａｌ ｒｅｓｕｌｔｓ ｏｎ ｓｐａｒｓｅ ｒｅｐｒｅｓｅｎｔａｔｉｏｎｓ ｏｆ ｍｕｌｔｉｐｌｅ ｍｅａｓｕｒｅｍｅｎｔ ｖｅｃｔｏｒｓ，”ＩＥＥＥ Ｔｒａｎｓ．Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓ．，ｖｏｌ．５４，ｎｏ．１２，ｐｐ．４６３４－４６４３，２００６．

【非特許文献22】Ｗ．－Ｋ．Ｍａ，Ｔ．－Ｈ．Ｈｓｉｅｈ，ａｎｄ Ｃ．－Ｙ．Ｃｈｉ，“ＤＯＡ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｑｕａｓｉ－ｓｔａｔｉｏｎａｒｙ ｓｉｇｎａｌｓ ｗｉｔｈ ｌｅｓｓ ｓｅｎｓｏｒｓ ｔｈａｎ ｓｏｕｒｃｅｓ ａｎｄ ｕｎｋｎｏｗｎ ｓｐａｔｉａｌ ｎｏｉｓｅ ｃｏｖａｒｉａｎｃｅ：ａ ＫｈａｔｒｉＲａｏ ｓｕｂｓｐａｃｅ ａｐｐｒｏａｃｈ，”ＩＥＥＥ Ｔｒａｎｓ．Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓ．，ｖｏｌ．５８，ｎｏ．４，ｐｐ．２１６８－２２１８，２０１０．

【非特許文献23】Ｊ．Ｍ．Ｋｉｍ，Ｏ．Ｋ．Ｌｅｅ，ａｎｄ Ｊ．Ｃ．Ｙｅ，“Ｃｏｍｐｒｅｓｓｉｖｅ ＭＵＳＩＣ：ｒｅｖｉｓｉｔｉｎｇ ｔｈｅ ｌｉｎｋ ｂｅｔｗｅｅｎ ｃｏｍｐｒｅｓｓｉｖｅ ｓｅｎｓｉｎｇ ａｎｄ ａｒｒａｙ ｓｉｇｎａｌ ｐｒｏｃｅｓｓｉｎｇ，”ＩＥＥＥ Ｔｒａｎｓ．Ｉｎｆ．Ｔｈｅｏｒｙ，ｖｏｌ．５８，ｎｏ．１，ｐｐ．２７８－３０１，Ｊａｎ．２０１２．

【発明の概要】

【発明を実施するための形態】

【0015】

本発明はこのような点に鑑みてなされたものであり、その目的は、分離信号と源信号到来方向を算定する到来方向別信号波形分離方法を提供することにある。

【0016】

複数の源信号が混在する環境において到来方向別に信号波形を分離し、源信号到来方向を推定することを目的とする。

【課題を解決するための手段】

【0017】

上述した課題を解決し、高精度な信号分離と源信号到来方向推定を実現するために、本発明では、源信号空間事前情報（源信号到来方向事前情報）を用いて高精度に分離行列を算定することができる分離信号算定部と分離行列の算定誤差に影響されにくいクラスタリング部及び分離行列算定結果から源信号到来方向を推定する源信号到来方向算定部と算定精度が低い分離信号を補正する分離信号補正部を有することを特徴とする演算方法を提供することを目的とする。

【0018】

本発明の方向別信号波形分離方法は、源信号到来方向が既知の周波数ビンにおいて、到来方向を事前情報として用いて分離行列算定精度を向上させることのできる分離信号算定部と、混合モデルとしてフォンミーゼス・フィッシャー分布を用いて、その確率密度関数を最大化することを目標として、信号基底ベクトルをクラスタリングするクラスタリング部と、周波数ビン・源信号ごとに分離行列（分離信号）の算定精度を評価し、周波数ビンごとの観測信号のパワースペクトル等の各周波数ビンの重要度指標と算定精度評価値を掛けた値を重み値として、全周波数ビンを対象に、周波数ビンごとの源信号到来方向の重み付き平均を取ることで、精度のよい方向推定結果を得ることのできる源信号到来方向算定部を有することを特徴とする。

【0019】

本特許の信号分離演算は、源信号の空間事前情報が得られている場合と得られていない場合とで、分離信号算定部で実施する演算をシステム内で選択的に変更し、空間事前情報が得られていない場合に、クラスタリング部及び源信号到来方向算定部及び分離信号補正部の演算を実施する。

【図面の簡単な説明】

【0020】

【図1】方向別信号波形分離方法の演算を実施する装置の概略図である。

【図2】方向別信号波形分離方法の演算処理の流れを表すフローチャートである。

【図3】源信号到来方向事前情報を用いて概略的に源信号波形を算定する演算処理のフローチャートである。

【図4】源信号到来方向事前情報と分離行列及び源信号統計的分散の推定値を更新して詳細に源信号波形を算定する演算処理のフローチャートである。

【図5】クラスタリング部で行う基底ベクトルのクラスタリング演算の中の初期クラスタリング部で行う演算処理のフローチャートである。

【図6】クラスタリング部で行う基底ベクトルのクラスタリング演算の中のクラスタ決定部で行う演算処理のフローチャートである。

【図7】源信号到来方向算定部で行う演算処理のフローチャートである。

【発明を実施するための形態】

【0021】

以下に図面を用いて、本発明に係る方向別信号分離演算の実施形態を説明する。図１は、到来方向別信号波形分離方法を実行する装置の概略図である。図１に示すように、本実施形態における方向別信号分離方法を実行するこの到来方向別信号波形分離装置では、マイク等を用いて信号を取得し、観測信号分離演算装置に信号データを送信する。観測信号分離演算装置では、観測信号から到来方向別の源信号波形とそれらの到来方向を推定する演算を実施し、その演算過程は、分離信号算定部、基底ベクトル算定部、クラスタリング部、源信号到来方向算定部、分離信号補正部より構成される観測信号分離演算装置で実施される演算により、源信号到来方向とそれに対応する源信号波形を得た後、その演算結果を汎用コンピュータ又は専用のソフトウェアを用いて表示する。

【0022】

図２は、本実施形態における方向別信号分離演算の中の観測信号分離演算装置で実行される演算フローである。図２に示すように、本実施形態における観測信号分離演算は、源信号の到来方向が事前にわかっている時とわかっていない時とで場合分けされる。源信号の到来方向が未知の時、分離信号算定部で観測信号を到来方向別の源信号波形へと分離し、基底ベクトル算定部で得られる波形特徴量をクラスタリング部でカテゴリ別に分類し、分離信号補正部で推定誤差の修正を行って源信号波形を推定すると同時に、源信号到来方向算定部でそれらの到来方向を推定する。源信号の空間情報が事前に得られている時は、分離信号算定部で源信号到来方向の事前情報を用いてより高精度に源信号を到来方向別に算定する。

【0023】

〔信号波形取得部〕
信号波形取得部は２つ以上のマイクロフォンなどで構成される。信号波形を受信し、信号データを観測信号分離演算装置へ送る。

【0024】

〔観測信号分離演算装置における波形分離演算〕
以下に、観測信号分離演算装置で実施される波形分離演算の概要を示す。詳細は後ほど示す。

【0025】

〔分離信号算定部〕
受信センサで取得した波形データを周波数領域に変換し、続いて分離信号算定部にて分離行列を求める。分離行列の生成には、独立成分分析や独立低ランク行列分析などを用いる。これらの手法では、分離演算の出力の同時分布と周辺分布の積とのＫＬ情報量を目的関数に設定し、目的関数が最小となるように、分離行列を更新して算定する。

【0026】

源信号の位置が目視で確認できる場合や事前に特定の周波数ビンに対してミュージック法等を適用して源信号到来方向を求めていた場合など、源信号到来方向が事前にわかっている場合もある。このように源信号到来方向が事前にわかっている時、事前情報を用いると高精度に源信号波形を推定することができる。源信号の空間事前情報が得られている場合の分離行列と分離信号の算定フローを図３に示す。まず源信号と観測センサ間の周波数伝達関数から得られる情報を用いて分離行列と源信号波形を概略的に算定する。図４には分離行列と源信号波形を概略的に推定する方法の一例をまとめていて、本特許の請求項２に記載の事項と対応する。概略的に算定した分離行列と源信号波形を用いて、分離行列と源信号波形を詳細推定する。この分離行列等の詳細推定方法は請求項３及び請求項４に記載の事項と対応する。以下で分離行列等の詳細推定方法を初めに説明する。

【0027】

源信号到来方向が事前にわかっている場合に源信号を詳細推定する方法についてまとめる。混合信号から源信号波形を推定するための目的関数に対して、源信号到来方向の事前情報を制約条件として加えた値を新たに目的関数に設定し、この目的関数が最小となるように、分離行列を更新することで、高精度な信号分離を達成できる。混合信号から源信号波形を推定するための目的関数には、例えば分離演算の出力の同時分布と周辺分布の積とのＫＬ情報量などがある。請求項３では源信号到来方向を事前情報とする制約条件の設定方法についてまとめている。源信号到来方向とセンサ配置位置が既知の時、式（９）を用いて、各センサから源信号までの周波数伝達関数を算定することができる。この周波数伝達関数をセンサ番号が行ベクトル、源信号番号が列ベクトルとなるように配置した行列は、混合行列といい、式（１０）で定義される。混合行列が分離行列の逆行列であるという情報を制約条件として設定することが請求項３に記載の制約条件の特徴である。混合行列が分離行列の逆行列である時、分離行列と混合行列の積は、単位行列になる。このことは、同一の源信号番号を有する分離行列の行ベクトルと混合行列の列ベクトルの内積が１になり、異なる源信号番号を有する分離行列の行ベクトルと混合行列の列ベクトルの内積が０になることを意味する。そのため、式（４１）に示す源信号到来方向に基づく制約条件では、上記の内積の値に対して制約をかけている。この分離行列の行ベクトルと混合行列の列ベクトルの内積に対する制約条件について請求項３では、記載している。グレーティングローブが可視領域に入らないようにするために、センサは、最小波長の半分以下となるように設置されるが、この条件を満たす時、混合行列は、傾向として低周波数になるほど悪条件になる。そのため、異なる源信号番号を有する分離行列の行ベクトルと混合行列の列ベクトルの内積が０になるという情報に基づく制約条件の強さを調整する。すなわち、混合行列が悪条件になるほど、源信号到来方向に基づく制約条件の強さが弱まるように、周波数ビンごとに係数を調整する。この係数の調整により、混合行列の信頼度が高い高周波数側において、分離行列と混合行列の積が、単位行列になるという源信号到来方向に基づく制約条件が重視される。この混合行列の信頼度は、センサ配置と信号の解析対称周波数帯域により異なるため、一概に低周波数で信頼度が落ちると断定することはできない。そのため、状況に応じて、混合行列の信頼度を評価し、式（４２）の係数設定方法を変更することが望ましい。

【0028】

請求項３に記載の制約条件を用いるためには、源信号の分散の大きさを概略的に推定しておく必要がある。分散の推定方法は式（１１）～式（３６）にまとめている。周波数ビン・時間フレームごとの分散を推定するために、概略的に分離信号波形を求めた後、ウィナーマスクを用いて源信号の統計的分散を推定する。この概略的に分離信号波形を算定する方法が請求項２と対応し、具体的な算定手順は次のとおりである。源信号到来方向が既知の時、周波数伝達関数を列ベクトルに並べた混合行列が式（１０）より求まる。この混合行列と観測信号が既知の条件において、分離信号の事後確率の対数尤度は、式（２８）で表される。この事後確率を最大化する分離信号波形を式（２９）に示すように、混合行列に正則化項を加えた行列の逆行列を観測信号に掛けることで算定する。請求項２では、周波数ビンｆごとに分離信号の分散を定義して更新式を求めているが、源信号の分散は時間周波数点ごとに異なるとする表記σ_ｎ（ｆ，τ）や周波数ビン間で一定とする表記であるσ_ｎ（τ）に変更することも可能で、そのような場合も本特許の請求範囲とする。請求項２に記載の分離信号算定方法は、それ自身で分離行列を算定し、算定結果を周波数領域又は時間領域で出力することも可能で、このような場合も本特許の請求範囲とする。

【0029】

源信号到来方向を事前情報として分離行列を算定する場合、その初期値の取り方が重要となる。源信号到来方向が既知の時、センサと源信号の位置関係から式（９）を用いて、周波数伝達関数を算定することができる。式（９）では、基準センサでの位相遅れを０となるようにして、基準センサを視点としているが、視点は任意の位置に設定してよい。この周波数伝達関数をセンサ番号が行ベクトル、源信号番号が列ベクトルとなるように配置することで、式（１０）に示す混合行列が求まる。この混合行列を分離行列の初期値として、解探索を行うと効率的に最適解に至ることができる。又、源信号が互いに近接して存在する場合など、混合行列が悪条件の時は、目的関数に正則化項を加えてから混合行列のムーアペンローズの擬似逆行列を計算することで、安定した分離信号の解が得られる。この方法は、請求項２に記載の内容であり、式（３２）を用いることで実現できる。目的関数を分離行列で偏微分して求まる目的関数の勾配降下方向に分離行列を繰り返し更新して最適解を探索する。このように、混合行列の逆行列又は目的関数に正則化項を加えた混合行列のムーアペンローズの擬似逆行列を分離行列の初期値として、分離行列を勾配法により繰り返し更新する手法を請求項４で記載している。

【0030】

源信号到来方向を事前情報として分離行列を算定する場合、その初期値の取り方が重要となる。源信号到来方向が既知の時、センサと源信号の位置関係から式（９）を用いて、周波数伝達関数を算定することができる。式（９）では、基準センサでの位相遅れを０となるようにして、基準センサを視点としているが、視点は任意の位置に設定してよい。この周波数伝達関数をセンサ番号が行ベクトル、源信号番号が列ベクトルとなるように配置することで、式（１０）に示す混合行列が求まる。この混合行列を分離行列の初期値として、解探索を行うと効率的に最適解に至ることができる。又、源信号が互いに近接して存在する場合など、混合行列が悪条件の時は、混合行列の２乗に正則化項を加えてから逆行列を計算することで、安定した分離信号の解が得られる。この方法は、請求項２に記載の内容であり、式（３０）を用いることで実現できる。解の探索は目的関数を分離行列で偏微分して求まる勾配の降下方向に分離行列を繰り返し更新して行う。このように、混合行列の逆行列又は混合行列の２乗に正則化項を加えた行列の逆行列を分離行列の初期値として、分離行列を勾配法により繰り返し更新する手法を請求項４で記載している。

【0031】

このように、制約条件を加えた目的関数の解空間は複雑化し、分離行列の更新が不安定化する傾向がある。そのため、請求項５に記載の方法で、解空間が複雑化した時に、安定的に分離行列を更新する方法を提案する。まず、分離行列の中を行ベクトルごと又は任意の範囲に区分けして、その１つ１つの範囲を行列ブロックと定義する。通常の独立成分分析では、分離行列の目的関数の勾配降下方向への更新量を学習率という定数パラメータを設定することで、予め決めておく。これに対して、提案方法では、行列ブロックごとに、分離行列の勾配降下方向への更新量を調整するために、学習率を更新する。そのために、式（４７－１）と式（４７－４）を用いて、行列ブロックごとに、現在のタイムステップまでの目的関数の勾配の２乗和の累積値と、分離行列の更新量の２乗和の累積値をそれぞれ算定する。現在のタイムステップにおける目的関数の勾配の２乗和の累積値と、１つ前のタイムステップにおける分離行列の更新量の２乗和の累積値を用いて、行列ブロックごとに、学習率を算定する。１つ前のタイムステップにおける分離行列の更新量に定数を掛けた値に対して、行列ブロックごとに求めた学習率を目的関数の勾配に掛けた値を加算することで、式（４７－３）で、現在のタイムステップにおける分離行列の更新量を算定する。このように、式（４７－１）と式（４７－４）目的関数の勾配の二乗と分離行列更新量の二乗について移動平均をとり、それらの値を基に式（４７－３）で勾配降下方向への更新量を調整するという一連の計算フローを繰り返すことで、最適な分離行列を求める。

【0032】

〔基底ベクトル算定部〕
分離行列を用いて波形データを分離する。又、分離行列又は信号における方向等の特徴量より基底ベクトルを算定する。

【0033】

〔クラスタリング部〕
基底ベクトルを同一の方向又は信号特徴量ごとにクラスタリングし、信号番号とクラスタ番号の対応付けを行う。クラスタごとに時間周波数点ごとの分離信号をまとめることで分離信号の時系列データを得る。

【0034】

まず、クラスタ算定で用いる概念や用語、変数の定義づけを行う。基底ベクトルの分布が独立なＣ個の確率密度関数の線形和で表される混合モデルを仮定する。本実施形態では、混合モデルの確率密度関数として、球面上の正規分布であり、方向データのみを持つ式（５５）に示すフォンミーゼス・フィッシャー分布を使用する。Ｃ個の確率密度関数が有する未知パラメータの集まりをＣ個のモデルパラメータと定義し、ｃ番目の分布に含まれるモデルパラメータの集まりを第ｃクラスタと呼ぶ。Ｃ個の確率密度関数の直列和で表される混合モデルにおいて、ｃ番目の確率密度関数の重要度をクラスタの混合重みπ_ｃ、全基底ベクトルの第ｃクラスタの平均方向μ_ｃ周りのばらつきを方向集中度κ_ｃと定義する。又、周波数ビンｆ・源信号ｎの基底ベクトルが第ｃクラスタに所属する確率を帰属度ｕ_ｃｎ（ｆ）と定義する。フォンミーゼス・フィッシャー分布の混合モデルでは、混合重みπ_ｃ、方向集中度κ_ｃ、平均方向μ_ｃの３種類がモデルパラメータと定義される。又、基底ベクトルのクラスタへの帰属度ｕ_ｃｎ（ｆ）を用いて、各クラスタへの所属確率を表現する。この帰属度が補助関数と定義される。確率密度関数を対数化した値にモデルパラメータと補助関数に関する制約条件にラグランジュ乗数を掛けた値を加算して目的関数を設計する。この目的関数を最大化するモデルパラメータと補助関数の組み合わせが最適解となる。本実施形態では、目的関数を直接最大化する代わりに、補助関数とモデルパラメータの内、一方が既知という前提で固定して、他方に関して目的関数を最大化して、解を更新するステップを交互に繰り返すことで、目的関数が最大となる解を探索する。

【0035】

図５及び図６は、基底ベクトルクラスタリングの計算フローである。本実施形態のクラスタ算定手順は、図５に示す初期クラスタリング部と図６に示すクラスタ決定部から成る。初期クラスタリング部では、補助関数（帰属度）を連続値として扱い、補助関数とモデルパラメータを交互に更新することで、目的関数が最大となる解を算定する。具体的な演算プロセスを説明する。補助関数とモデルパラメータの内、一方が既知という前提で固定して、他方に関して目的関数を偏微分して、解を更新するステップを交互に繰り返すことで、解を更新し、最適解を探索する。まずモデルパラメータが既知という条件の基で固定して、補助関数に関して目的関数を最大化することで、補助関数を更新する演算を実施する。次に補助関数が既知という条件の基で固定して、モデルパラメータを更新する演算を行う。これら２種類の演算を交互に繰り返すことで、最適解を探索する。モデルパラメータを更新する演算は、モデルパラメータの内、混合重みπ_ｃ、平均方向μ_ｃに関してそれぞれ目的関数を偏微分して勾配がゼロとなる解を探索する演算と、方向集中度κ_ｃを近似式又は勾配法又は最尤推定法を用いて算定する演算により構成される。モデルパラメータの内、混合重みπ_ｃと平均方向μ_ｃを更新する演算を式（５９）に示す。又、方向集中度を近似式と勾配法により算定する手順を式（６０）～式（６３）に示す。この方向集中度を求める演算方法は請求項８で記載されているため、後程詳しく述べる。方向集中度は、この他の近似式や最尤推定法を用いて算定することもできるが、請求項７に示すモデルパラメータの更新方法で混合重みと平均方向、帰属度を更新している場合は、請求項６に記載の範囲とする。式（５８）～式（６３）に示すモデルパラメータと補助変数の更新式は、基底ベクトルが周波数ビン・源信号の２次元構造を成す場合について示している。この周波数ビン・源信号の２次元構造のインデックスをもつ変数を更新する演算は、従来にはなかったパターンであり、新規性を主張できるポイントである。又、信号特徴量から求めた基底ベクトルをフォンミーゼス・フィッシャー分布によるクラスタ分析方法で同一源信号ごとにまとめるという試みはこれまでなされていない。更に各々の基底ベクトルに対して重みα_ｎ（ｆ）を掛けて、基底ベクトルの重要度が互いに異なるケースに対応できるように拡張している点も従来にはない新たなクラスタ分析方法である。時間周波数点単位で信号特徴量等から基底ベクトルを求めてクラスタ分析する時は、式（５８）～式（６３）の更新式において、周波数ビンのインデックスを時間フレーム・周波数ビンへと拡張することで対応できる。そのため本特許の請求項６では主に、基底ベクトルが周波数ビン・源信号の２次元構造を成す場合と、時間フレーム・周波数ビン・源信号の３次元構造を成す場合との２パターンを想定している。

【0036】

これらの演算を繰り返すことで得られる初期クラスタリング部における最適解の内、モデルパラメータをクラスタ決定部の初期値に設定する。クラスタ決定部では、モデルパラメータの初期値を取得し、同じ周波数ビンに属する異なる源信号インデックスを有する基底ベクトルは、必ず異なるクラスタに属するという先見情報に基づく制約条件を設けながら、目的関数が最大となる解を探索する。クラスタ決定部では、目的関数が最大となる解を探索する際、補助関数（帰属度）を０又は１の離散値に限定する。すなわち、帰属度の算定において式（５８）を用いず、代わりに、帰属度の値を０又は１の離散値に限定し、先見情報に基づく制約条件を設けながら、目的関数が最大となる帰属度の組み合わせを探索する。帰属度の組み合わせ探索において、総当たりで探索してもよいが、アルゴリズム１に示す集合分割の概念を用いると、解探索を効率化できる。このように、帰属度の値を０又は１に離散化した状態で、同じ周波数ビンに属する異なる源信号インデックスを有する基底ベクトルは、必ず異なるクラスタに属するという先見情報に基づく制約条件を設けながら、モデルパラメータと補助関数を交互更新する演算が本特許の請求項７と対応する。

【0037】

基底ベクトルクラスタリングにおいて方向集中度は、第１種変形ベッセル関数の中にも含まれているため、解析的には求まらない。そこで請求項８では、近似式と勾配法を用い、次の手順に従って、方向集中度を算定する方法を示している。まず式（５８）により、階数の１つ異なる第１種変形ベッセル関数の比率を基底ベクトルと帰属度を用いて、概略的に算定する。この概略値に対して、階数の１つ異なる前記第１種変形ベッセル関数の比率を連分数に展開することで得られる式（６１）又は式（６３）の近似式を適用し、方向集中度の近似値を算定する。その後、方向集中度の近似値を階数の１つ異なる第１種変形ベッセル関数にそれぞれ代入して、階数の１つ異なる第１種変形ベッセル関数の比率の詳細値を算定する。第１種変形ベッセル関数の比率の詳細値を方向集中度の近似値で微分した値を式（６１）又は式（６３）により算定した後、第１種変形ベッセル関数の比率の詳細値と概略値の誤差を最小化することを目標として、ニュートン法を適用した式（６２）により、方向集中度を更新する。更新した方向集中度を再度、階数の１つ異なる第１種変形ベッセル関数にそれぞれ代入して、階数の１つ異なる第１種変形ベッセル関数の比率の詳細値を更新する。更新した階数の１つ異なる第１種変形ベッセル関数の比率の詳細値を更新した方向集中度で微分した値を勾配とし、第１種変形ベッセル関数の比率の詳細値と概略値の誤差を最小化することを目標として、ニュートン法を適用することで求まる式（６２）により、方向集中度を更新する。このニュートン法による方向集中度の更新演算を何度か繰り返すことで、精度良く方向集中度を推定できる。以上で述べた計算手順による方向集中度算定手順が本特許の請求項８に記載の内容と対応する。

【0038】

このように、初期クラスタリング部及びクラスタ決定部における解の探索は、補助関数とモデルパラメータの内、一方が既知という前提で固定して、他方に関して最適な解を探索するステップを交互に繰り返す算定手順により行う。初期クラスタリング部の演算には、源信号数が多くなった等の理由で、補助関数とモデルパラメータの数が増えた時に、クラスタ決定部において、解が局所解に誤って収束するのを防ぐことができるという効果がある。但し、源信号数が少ない場合は乱数により、モデルパラメータの初期値を設定しても支障はなく、その場合は、クラスタ決定部の演算を単独で用いることもできる。これらの最適解探索過程において、方向集中度が第１種変形ベッセル関数の中に入り込んでいて算定に工夫を擁する。初期クラスタリング部では、方向集中度を全てのクラスタで一律とした近似式（６３）を用いるのが望ましいのに対して、クラスタ決定部では、クラスタごとに方向集中度を定義した近似式（６１）を用いるのが望ましい。

【0039】

〔源信号到来方向算定部〕
先ほどクラスタリングにより対応付けを行った信号番号とクラスタ番号の対応関係を用いて、各基底ベクトルから算定される信号の到来方向に対して、クラスタごとに全周波数ビンに対する重み付き平均を取ることで、全周波数点における源信号到来方向を得る。

【0040】

図７は、源信号到来方向算定の計算フローである。分離行列の逆行列を求め、その列ベクトルから式（８）を用いて、周波数ビン・源信号ごとの周波数伝達関数を求める。この源信号と観測点間の情報を含む周波数伝達関数の推定値とセンサの配置位置の情報から、センサ間の位相差に関する情報を基に導かれる式（７４）を用いて、周波数ビン・源信号ごとの源信号到来方向を算定する。又、各周波数ビンにおける信号分離精度が低い信号の信号算定結果や源信号到来方向推定結果に及ぼす影響が少なくなるよう、式（７１）により連続値マスクを算定する。信号番号とクラスタ番号の対応付け結果から式（７５）によりマスク掛けした周波数ビン・源信号ごとの源信号到来方向の重み付き平均を求めることで、全周波数ビンにおける最終的な源信号到来方向を決定する。重み付き平均で用いる基底重みは、式（７６）に示すように周波数ごとの観測値のパワースペクトルを指標とすることもできるが、どの周波数も同じ信頼度なら全て１とすることもできる。他にも基底重みとしては、信号部分空間の固有値の和を指標とすることもでき、この時近似的に反射せずに直接、観測点まで届いた到来波のパワーの総和をとっていることになる。

【0041】

〔分離信号出力部〕
信号波形分離部から源信号ごとに分離された信号情報を受け取り、波形の時系列データまたは周波数領域でのパワースペクトルの形式でモニターに表示する。表示した波形データは、消去若しくは保存をユーザの選択に応じて行う。又、源信号到来方向の推定結果を受け取り、事前にカメラで撮影した画像と重ね合わせることで、源信号到来方向をわかりやすく可視化する。分離波形データの表示や保存を行う装置及び源信号到来方向の推定結果を表示する装置には、汎用コンピュータを用いることもできるし、専用のハードウェアを用いることもできる。前者の場合には、汎用コンピュータのハードウェアとソフトウェアとが協働して構成される。

【0042】

以下に、本実施例における到来方向別源信号分離演算の具体的な演算手順を説明する。演算は、先ほど説明した通り、分離信号算定部、基底ベクトル算定部、クラスタリング部、源信号到来方向算定部に分かれており、算定プロセスごとに説明する。

【0043】

〔分離信号算定部〕
時間領域の源信号ｓ_ｎ（ｔ）、分離信号ｙ_ｎ（ｔ）、観測信号ｘ_ｍ（ｔ）の周波数領域での複素数成分をベクトル表現により次の式で表す。

【0044】

【数1】

ｆ＝１，…，Ｆ、τ＝１，…，Ｔは、それぞれ周波数ビン及び時間フレームのインデックスを表す。又、ｎ＝１，…，Ｎで全周波数帯域における帯域ごとの信号数の最大値、ｍ＝１，…，Ｍで観測センサの番号を表す。

【0045】

時間領域の観測信号ｘ_ｍ（τ）に対して窓の位置を少しずつ移動させながら、窓付きフーリエ変換を施し、時間周波数点ごとの観測データｘ（ｆ，τ）を生成する。時間周波数点におけるｘ（ｆ，τ）を分離行列により分離することで、源信号ｓ（ｆ，τ）の推定値である分離信号ｙ（ｆ，τ）を算定する。

【0046】

分離行列の算定方法としては、独立成分分析や独立低ランク行列分析など様々な方法を用いることができる。波形の事前情報がわからない時、例えば建物内での異音を対象とする場合や街中で騒音波形を方向別に取得したい場合は、独立成分分析が適している。一方、音楽信号や人の声のような低ランク構造をもつ音信号を方向別に分離する場合は、独立低ランク行列分析が適している。源信号到来方向に関する事前情報が得られていない場合においては、独立成分分析による分離行列の更新方法のみを示す。

【0047】

【数2】

周波数ビンごとにＮ×Ｍ次元の行列構造を有する分離行列Ｗ（ｆ）を生成する。

【0048】

独立成分分析による分離行列は、目的信号の統計的な情報を元にした最適化手法を用いて求める。具体的には、式（３）で与えられる分離演算の出力の同時分布と周辺分布の積とのＫＬ情報量を最小化する分離行列を、勾配法を用いて最急降下方向に解探索することで算定する。

【0049】

【数3】

Ｅ［ｘ（ｆ，τ）］は、観測信号の時間フレームτ方向の平均値を算定する処理である。Ｅ［ｘ（ｆ，τ）ｙ^Ｈ（ｆ，τ）］は、観測信号と分離信号の時間フレーム方向の共分散行列を算定する処理を表す。ｐ_ｎ（ｙ_ｎ（ｆ，τ））は、目的信号の統計的な情報を表す確率密度関数である。確率密度関数をｙ_ｎ（ｆ，τ）で微分し、確率密度関数のＮ次元勾配ベクトルであるスコア関数φ（ｆ，τ）を算定する。式（４）に示すＫＬ情報量をＷ^Ｈ（ｆ）で偏微分することにより得られる式の勾配ベクトルを用いて、分離行列を最適化する。

【0050】

【数4】

ηは、分離行列の学習速度を制御する定数であり以降、学習率と称する。Ｉ_ＭはＭ次元単位行列である。

【0051】

【数5】

スコア関数には、様々な種類があり、分離信号の確率密度関数がガウス分布に従うと仮定する場合は、式（５）となる。σ_ｎ（ｆ）は、分離信号の確率分布の統計的分散である。

【0052】

【数6】

分離信号の確率密度関数が双曲線余弦に従うと仮定する場合は、スコア関数は、式（６）となる。この確率密度を用いると、優ガウス分布に従う源信号の高次モーメントを再現できる。

【0053】

【数7】

分離信号の確率密度関数がラプラス分布に従うと仮定する場合は、スコア関数は、式（７）となる。

【0054】

源信号到来方向を得ることに特化した手法として、ミュージック法がある。もし、ミュージック法を適用できる条件が揃っていれば、ある特定の周波数ビンにおける源信号到来方向の組み合わせを知ることができる。又、物理的に源信号位置がわかっている時もある。このように、ある特定の周波数ビンにおいて、源信号到来方向が既知の時にその事前情報を利用することで、信号分離精度向上を図る手法を提案する。この手法を信号分離演算の中にオプションとして組み入れることで、システム内で選択的に用いることができるようにする。

【0055】

基準となるセンサの番号Ｍ′を観測点番号１～Ｍから選択する。以降、このセンサ番号を有するセンサを基準センサと称す。信号伝達モデルとして遠距離場モデルを仮定すると、観測点間での源信号からの距離の差はそれぞれの観測点と源信号との距離に比べて十分に小さいため、観測点間での到来波の振幅減衰の差は微小になり無視できる。このことを考慮して、源信号からセンサに到達するまでの減衰の影響を無視し、源信号と観測点との位置関係及び周波数ビンｆにおける到来波の波長λ（ｆ）をパラメータとして式（８）により周波数伝達関数を近似する。

【0056】

【数8】

ｐ_ｍはｍ番目のセンサの位置ベクトル、ｑ_ｎは基準センサから視た目的信号への方向を表す単位ベクトルである。式（８）において、基準センサの番号をＭ’とし、基準センサとその他センサとの周波数伝達関数の比をとることで、式（９）を得る。このように比をとることで、基準センサにおける波の位相を基準にして、その他のセンサにおける位相を算定していることになる。

【0057】

【数9】

上式においてセンサ間の位置ベクトルの差と源信号到来方向への単位ベクトルとの内積の値を波速で割った値がセンサ間の信号到達時間差（遅延差）となる。

【0058】

Ｎ個の源信号の内、Ｎ’個の源信号到来方向が事前情報としてわかっている場合を考える。源信号到来方向が既知の時、センサ間の遅延差を基に式（９）により周波数伝達関数を算定する。センサＭ’を基準に空間事前情報が得られている源信号に対応する周波数伝達関数を縦方向に並べたＭ次元アレイマニフォールドベクトルを式（１０－２）で定義する。混合系が線形時不変であり、時間周波数領域での複素瞬時混合で表現できると仮定すると、アレイマニフォールドベクトルを用いて、Ｍ×Ｎ’次元複素混合行列Ａ（ｆ）を式（１０－１）で定義できる。

【0059】

【数10】

ここで、到来方向が既知の源信号インデックスをｎ’＝１，…，Ｎ’としている。式（１０－２）で、ベクトル成分の分母に源信号から基準センサまでの周波数伝達関数があるが、この基準センサ番号は、分離信号を出力したいセンサ番号とすればよい。

【0060】

もし源信号から各センサまでの周波数伝達関数を誤差なく算定できるならば、それらを基に作成した混合行列の逆行列を分離行列として分離信号を算定できる。しかし源信号到来方向の観測誤差が生じた時、混合行列は悪条件な方程式で成り立っているので、僅かな誤差が結果に大きな影響を与える。又、遠距離場においてセンサから源信号位置までの距離を正確に測定することは困難なので、信号振幅の減衰の影響を無視して式（９）により周波数伝達関数を基準センサからの相対値で近似している。そのため、混合行列の逆行列を用いると分離波形には必ず誤差が含まれることになる。

【0061】

まず初めに、Ｎ＝Ｎ’すなわち全ての源源信号到来方向が事前にわかっている場合において、源信号の統計的分散を算定し、到来方向の情報に加えることで、分離信号を算定する方法をまとめる。

【0062】

【数11】

これまでと同様に、混合系の近似モデルを仮定すると、分離信号と観測信号の相関行列は、混合行列を用いて式（１１）で関連づけられる。観測信号の相関行列を固有値分解すると次の式が得られる。

【0063】

【数12】

▲Ｑ▼（ｆ）は、Ｍ次元固有ベクトルｑ_ｍ（ｆ）を並べたＭ×Ｍ次元の行列であり、Ｍ個の列ベクトルの内、Ｎ＜＝Ｍを満たす固有ベクトルが信号部分空間を張る。Ａ（ｆ）は、Ｍ個の固有値を値が大きい順に並べた対角行列である。式（１２）の固有値と固有ベクトルにおいて、源信号インデックスｎは周波数ビンごとに異なる源信号と対応していて一貫性がない。この交換不定性問題を解消するために、大きさが上位Ｎ個の固有値に対応する信号部分空間の固有ベクトルをその方向別にクラスタリングする。基底ベクトルを次の式で算定する。

【0064】

【数13】

Ｎ×Ｆ本の基底ベクトルｖ_ｎ（ｆ）をその方向が類視したペアごとにまとめたＮ×Ｆ個の大きさが同じクラスタに分類する。最も簡易なクラスタリング方法は、Ｋ－ｍｅａｎｓ法に基づく方法である。まず同じ源信号インデックスｎを有する基底ベクトルを一つのクラスタとみなして、各クラスタにおけるＭ次元平均方向ベクトルμ_ｃを算定する。基底ベクトルとその平均方向との差が最小すなわち内積が最大となるように、基底ベクトルを並び替え、クラスタを更新する。再度クラスタの平均方向を求める。クラスタに変化がなくなるまでこの交互に演算を繰り返す。演算過程を式（１４）に示す。

【0065】

【数14】

ｃ＝１，…，Ｃはクラスタインデックスである。ｕ_ｃｎ（ｆ）は周波数ビンｆにおける源信号インデックスｎの基底ベクトルが第ｃクラスタに属する時に１、属さない時に０を返す帰属度変数である。基底ベクトルは必ずいずれかのクラスタに所属するため、帰属度変数のクラスタインデックスｃ方向への和は必ず１となる（式（１５－１））。又、帰属度変数の源信号インデックスｎ方向への和は、周波数ビンｆの中に源信号ｎが存在する時に１、存在しない時に０となる。今回のように信号を同数のクラスタに分類する並び替えの問題、すなわちＣ＝Ｎが成り立つ場合は、クラスタ番号と源信号番号が１対１で対応するため、源信号インデックス▲ｎ▼方向への和が１となる（式（１５－２））

【0066】

【数15】

しかし、Ｋ－ｍｅａｎｓ法によるクラスタリングでは、基底ベクトルの本数が多い時、最終的に生成されるクラスタが初期のクラスタに依存してしまい、一意に解が定まらないという問題がある。そこで、後程「クラスタリング部」に示す混合モデルとしてフォンミーゼス・フィッシャー分布を仮定した演算によりクラスタ分析を行い、帰属度変数の集合ｕ_ｃｎ（ｆ）を決定することを推奨する。クラスタリングが完了すると、その結果を用いて固有ベクトルをクラスタごとに算定し、各クラスタに

【0067】

【数16】

式（１６）に示す平均方向ベクトルの成分を用いて基準センサＭ’から視たクラスタごとの源信号到来方向を次式で算定する。

【0068】

【数17】

クラスタごとの源信号到来方向が実際の到来方向ＤＯＡ_ｎと一致するように、全周波数ビンのクラスタインデックスｃを同一規則で並び替える。まず、式（１８）によりクラスタごとの源信号到来方向が実際の源信号到来方向と近い値となるように対応付けを行う。この対応付けによりクラスタ番号ｃを入力すると正しい源信号番号▲ｎ▼を返す順列集合｛Π（ｃ）｝を作成する。

【0069】

【数18】

この順列Π（ｃ）を用いて、固有値と固有ベクトルの並び替え（クラスタリング）を以下により行う。

【0070】

【数19】

並び替えが完了したら、源信号の統計的分散を周波数伝達関数と固有値分解の結果を用いて求める。分離信号が独立であれば、その相関行列は、次の式に示すように源信号の統計的分散を対角に並べた行列と一致する。

【0071】

【数20】

式（１１）と式（１２）の形状を比較すると明らかなように、▲Ｑ▼（ｆ）は無相関化された分離信号の混合行列と対応する。そのため、式（２０）を式（１１）に代入し 式（１２）と係数比較すると、源信号の統計的分散は以下により算定できる。

【0072】

【数21】

源信号到来方向の事前情報を用いて分離信号を概略的に算定する方法についてまとめる。以下に述べる統計的分散推定演算で用いられる分離信号算定方法は本特許の請求項２と対応する。

【0073】

源信号到来方向が既知の時に、到来方向を事前情報として用いて、分離信号を簡易的に算定することを試みる。源信号が混合した状態で、観測信号としてセンサにて受信される場合を考える。この時、周波数領域での混合モデルは以下で表せる。

【0074】

【数22】

源信号の推定値である分離信号は、混合行列の逆行列を用いて算定できる。

【0075】

【数23】

上付文字＋は、Ｍｏｏｒｅ－Ｐｅｎｒｏｓｅ疑似逆行列の表記であり、行列の行数と列数が等しい時に、通常の逆行列となる。解析対称周波数の中の低周波領域や源信号が近接した位置に存在する場合に、混合行列は、悪条件な方程式となる。そのため、式（２３）のように逆行列を用いる方法では、僅かな源信号到来方向の推定誤差が分離信号算定結果に大きな誤差をもたらす。この算定誤差をｅ（ｆ）と定義し、混合行列の悪条件性に対処するために混合モデルを等式制約とした次の最小化問題を解く。

【0076】

【数24】

観測信号がガウス分布に従うと仮定すると、源信号の確率密度関数は、次式で表せる。

【0077】

【数25】

雑音が各周波数ビンで独立に分散σ_Ｅ（ｆ）のガウス分布に従って発生する時、混合行列と分離信号が既知という条件下での観測信号の確率密度関数は以下で与えられる。

【0078】

【数26】

ベイズの定理より、混合行列と観測信号が既知という条件下での分離信号の確率密度関数は、

【0079】

【数27】

となり、Ｐｒ（Ａ（ｆ），ｘ（ｆ，τ））を一定とし、式（２７）に式（２２）、（２５）、（２６）を代入し、対数化すると、式（２８）が得られる。

【0080】

【数28】

上記の対数尤度を最大化する分離信号を求めるために、分離信号ｙ（ｆ，τ）で左辺を微分して極大点すなわち微分値が０となる分離信号を求めると、

【0081】

【数29】

Σ（ｆ）はＮ×Ｎ次元対角行列である。この式を解くには、源信号の分散と雑音の分散の２つを求める必要があるが、これらはいずれも未知である。そこで、雑音の分散を未知数として、式（２９）を微分した値を０と等値すると、

【0082】

【数30】

式（２９）の形状に着目すると、雑音と源信号の分散の比率が正則化項の強さを決定することがわかる。そこで源信号の分散を用いて、雑音の分散の初期値を次式で与える。

【0083】

【数31】

Ｋ_Ｆは観測ノイズの大きさを推定するパラメータであり、０．０１～０．１程度とするのがよい。但し、源信号の位置が互いに離れている場合や解析対称周波数が特定の帯域のみの場合など、混合行列が悪条件でない時は、０．０１より小さな値とすることもできる。ｐｓは、混合行列が低周波数になるほど悪条件となる性質を反映したパラメータで、０～１程度とするのがよい。計算は次の手順で行う。
まず源信号の分散として式（１２）に示す観測信号の相関行列の固有値から求めた値を用いる。源信号数が観測センサ数よりも多い過完備基底の条件では、固有値分解ができないため、源信号の分散構造に対する事前情報が得られない。そのような場合、源信号の分散は全て１とするとよい。雑音分散の初期値は、式（３０）により求める。次に、式（２９）と式（３０）を交互に計算して、分離信号の値と雑音分散の値を交互に更新する。この交互演算を全周波数ビンにおける雑音の分散の変化が小さくなるまで行う。この交互演算が収束した時の分離信号がその事後確率を最大化する最適解となる。

【0084】

【数32】

【0085】

更新が完了したら分離行列の逆行列に相当する次式により振幅補正を行う。これまで説明した演算による分離信号算定方法は、計算量が少なく高速である一方、源信号波形が実際よりも疎構造になる傾向がある。しかし、パワースペクトル形状の特徴をよく捉えることができる点や計算の簡易さを考慮すると、実用的と考えられる。

【0086】

この他にも分離行列の概算的な求め方には様々な方法が考えられる。最も簡易な求め方は、到来方向が既知の源信号のみセンサ間の遅れが補償されるように、基準センサと源信号間での周波数伝達関数を求め、その周波数伝達関数を源信号ごとに並べたベクトルｈ_ｎ（ｆ）の複素共役ベクトルを列ベクトルに並べた行列を分離行列とする方法である（式（３３－１））。この方法は混合行列が悪条件で、逆行列が精度よく算定できない時や到来方向が既知の信号源数Ｎ’が全ての信号源数Ｎよりも少ない時に用いるとよい。又、混合行列のムーアペンローズの擬似逆行列を分離行列とする方法も考えられる（式（３３－２））。この方法は、高周波数帯域などの混合行列の逆行列が精度よく算定できる時に用いるとよい。請求項４ではこれら分離行列の初期値の求め方についてまとめている。

【0087】

【数33】

【0088】

Ｉ_ｐはｐ次元単位行列である。源信号が楽器音や音声のような明確な周波数構造をもつ場合、統計的分散を源信号到来方向の事前情報を正則化項として加えることで更新し、より精度よく分離行列と分離信号を算定することができる。その方法を以降、まとめる。まずこれまで示した概略算定方法により算定される分離行列を用いて、分離信号を算定する。分離信号がＭＮＭＦ（非特許文献１３）の信号源モデルに従うと仮定した場合、分離信号のパワースペクトルが期待値の意味で保存されるため、ウィナーマスクＭ_ｎ ^（Ｇ）（ｆ，τ）は次の式で設計できる。

【0089】

【数34】

基準センサＭ’における観測信号のパワースペクトルをウィナーマスクによりセンサ内の源信号ごとのパワースペクトルに分解し、基準センサ内における源信号の統計的分散を次のいずれかの式で求める。

【0090】

【数35】

一般に上記のような分離信号のパワーの加法性は成り立たないが、観測信号の生成モデルが多変量複素ガウス分布に従う時、パワースペクトルの加法性が期待値の意味で保存される。上記に示すウィナーマスクによる分離信号の分散推定はこの性質を利用したものである。

【0091】

【数36】

上記の算定式は、分離信号のパワーが源信号パワーと等しいという条件式である。源源信号到来方向が全て既知でない場合、ウィナーマスクを設計できないので、上記の算定式を用いるとよい。式（３５）又は式（３６）により算定される源信号の統計的分散と概略算定により求めた分離行列を初期値に設定する。これらを更新して源信号推定精度向上を実現する手法を示す。

【0092】

源信号到来方向の事前情報を目的関数に制約条件として付加することで、信号分離を高精度に行う手法ついてまとめる。ＭＮＭＦでは、源信号が多変量複素ガウス分布に従うと仮定した時の多チャネル観測信号の生成モデル（非特許文献１３、非特許文献１５）に基づく目的関数が用いられる。この目的関数は、ガウス分布に従う源信号分散及びセンサと源信号間の周波数伝達関数が既知の条件における源信号の事後確率と等価である。提案手法では、ＭＮＭＦの目的関数に対して、空間相関行列をランク１行列で近似することで得られる目的関数Ｌ_ＭＮＭＦに、源信号到来方向に関する制約条件ＬＳ（Ｗ（ｆ））を加えた目的関数ＬＬ（Ｗ（ｆ））を定義する。このランク１行列で近似した目的関数Ｌ_ＭＮＭＦは独立低ランク行列分析（非特許文献１４）で用いられている。

【0093】

【数37】

ＭＮＭＦに基づく目的関数Ｌ_ＭＮＭＦに制約条件ＬＳ（Ｗ（ｆ））を加えた式（３７－１）の目的関数が最小となるように、分離行列を更新する。目的関数の勾配降下方向を求めるために、式（３７－１）をＷ^Ｈ（ｆ）で偏微分すると、次式が得られる。

【0094】

【数38】

右辺第１、２項がＭＮＭＦに基づく目的関数の勾配である。上式において右辺第１、２項のみＷ^Ｈ（ｆ）Ｗ（ｆ）を掛けて自然勾配に補正すると、最急降下方向は近似的に以下により与えられる。数学的には、式中の特定項のみを自然勾配に補正する演算は間違いなのだが、正則化項ＬＳ（Ｗ（ｆ））の次元が自然勾配補正後の次元と一致するように設定するため、実用上問題はない。

【0095】

【数39】

Ｉ_Ｍは、Ｍ×Ｍ次元単位行列である。

【0096】

【数40】

ＭＮＭＦに基づく目的関数を用いると、分離信号の確率密度関数の微分値であるスコア関数は式（３９）となる。この分布では、源信号の統計的分散が時間周波数点単位で定義されている。源信号到来方向の事前情報を基に求めた混合行列に、源信号が互いに独立であるという条件を加えて解を探索する。源信号到来方向事前情報の制約条件に基づく正則化項ＬＳ（Ｗ（ｆ））を次の式により表現し、ＭＮＭＦに基づく目的関数に加えることで、分離精度を向上させる。

【0097】

【数41】

ｅ_ｎは▲ｎ▼番目の要素が１、他要素が０のＭ次元単位ベクトルである。正則化係数における上付括弧文字のｃｙｃｌｅ＝１，…，ＣＹＣＬＥは、繰り返し計算の反復回数を表すインデックスである。式（４１－１）の第１項は、同じ源信号番号の分離行列の行ベクトルと混合行列の列ベクトルの内積が１になるという制約条件であり、目的方向から到来する信号を強調する効果がある。第２項は、異なる源信号番号の分離行列の行ベクトルと混合行列の列ベクトルの内積が０になるという制約条件になっていて、目的方向以外から到来する信号を遮断する死角形成効果がある。正則化係数ρ_Ｄ ^{（ｃｙｃｌｅ）}とρ_Ｕ ^{（ｃｙｃｌｅ）}（ｆ）は、各項の制約の強さを支配する定数である。これら２種類の正則化係数は次の関係式で関連付けておくのがよい。

【0098】

【数42】

ｐｓは、０～１程度の定数である。ある特定の２つのセンサを用いて、源信号到来方向を推定したい時に、グレーティングローブが可視領域に入らないようにするためには、センサ間隔は、解析対称の最大周波数で信号伝搬速度（音の場合は音速）を割った値すなわち最小波長の半分以下になるようにすべきである。よって、Ｍ個のセンサで信号を観測し、基準センサＭ’から視た源信号到来方向を推定する場合、基準センサと基準センサから最も距離の近いセンサとの距離が最小波長の半分以下であれば、源信号到来方向が一意に定まることになる。この条件を満たす時、例えば、基準センサと基準センサから最も距離の近いセンサに着目すると、低周波数になるほど可視領域が狭くなる。一方、その他のセンサと基準センサとの間に着目すると、グレーティングローブが可視領域に入らない周波数の中の最大値よりも小さい値の周波数に対して同様の関係が成り立つ。そのため、グレーティングローブが可視領域に入らず源信号到来方向が一意に定まる時、混合行列Ａ（ｆ）は、傾向として低周波数になるほど悪条件になる。この傾向を基に低周波になるほど、目的関数の制約条件に基づく正則化項における死角形成効果を弱める役割を式（４２）の正則化係数ｐｓは担っている。式（４１）の正則化項を分離行列の複素共役で偏微分すると、次の式が得られる。

【0099】

【数43】

Ｏ_ｐｑはｐ×ｑ次元の零行列を表す。制約条件を分離行列で偏微分して得られるＮ×Ｍ次元の行列におけるＮ’×Ｍ次元のＮ’×Ｍ次元の非零要素Ｐ_Ｄ（ｆ）及びＰ_Ｕ（ｆ）は、Ｍ次元アレイマニフォールドベクトルａ_ｎ’（ｆ）を用いて、式（４４）で算定できる。

【0100】

【数44】

Ｒ_ｎ（ｆ）は、空間相関行列をアレイマニフォールドベクトルでランク１近似したＭ×Ｍ次元エルミート行列である。このように到来方向が既知の源信号に対してのみ方向情報に基づく制約条件を加えることで、複数の源信号が混合した波形から到来方向が既知の信号波形を正確に抽出することが可能になる。

【0101】

式（４４）の正則化項を目的関数に加えると、解空間が複雑になり、自然勾配法による最急降下法に基づく分離行列の更新では、局所解に陥ったり、計算が不安定になって発散したりする。そこで、請求項５に記載の新たな分離行列更新手法を提案する。前回の分離行列の更新方向を慣性項として勾配降下方向に加えると共に、勾配降下方向への移動量を調整することで、急激な更新方向の変化を抑制し、分離行列の更新を安定化させることを試みる。分離行列は、行ベクトルなどの行列ブロックごとに更新する。式（３７）に示す目的関数を分離行列の複素共役で偏微分した値の方向から勾配降下方向を求め、分離行列を更新する。行列をいくつかの行列ブロックに分解して、行列ブロックごとに解を更新する。以下、同一の源信号番号ごとに行列ブロックを区切った場合、すなわち行列ブロックとして分離行列の行ベクトルを定義する場合について考える。

【0102】

局所解への誤った収束を防ぐために、行列ブロックごとに、目的関数の勾配に１つ前のタイムステップの分離行列の更新量を定数倍し、慣性項として加える。この慣性項により例えば、解の探索領域に深い谷があるようなケースでも容易に谷から脱出できるため、局所解に陥りづらくなる。さらに、計算の安定化を図るために、目的関数勾配と過去の分離行列更新量をそれぞれ二乗移動平均することで、学習率を調整する。具体的には、目的関数の勾配の二乗和を累積させた値が大きい時に更新量が小さくなり、小さい時に更新量が大きくなるように、さらに１つ前のタイムステップまでの分離行列の更新量の二乗和を累積させた値が大きい時に更新量が大きくなり、小さい時に更新量が小さくなるように調整する。このように分離行列の更新量を調整することは、独立成分分析で用いられる自然勾配法の更新式における学習率を解空間の形状に応じて自動で調整することと等価であり、計算の安定化に繋がる。具体的には、次の式で分離行列更新する。

【0103】

【数45】

分離行列の第ｎ行を式（４５）で定義する。ｅ_ｎはｎ番目の要素が１、他要素が０のＭ次元単位ベクトルである。

【0104】

【数46】

周波数ビンｆにおける式（３９）の目的関数勾配行列の第ｎ行を抜き出したＮ次元ベクトル▽Ｌ_ｎ（ｆ）を式（４６）で定義する。これらの定義式では、行列ブロックを分離行列の行ごとに定義しているが、その他の行列ブロック形状を取ることも可能である。一番単純なのは、分離行列全体を行列ブロックとみなし、分割数を１にする方法である。これらのケースでも請求項５に記載の事項と合致する場合は、本特許の権利範囲とみなす。周波数ビンｆ、源信号ｎごとに次の更新式を用いて分離行列を更新する。

【0105】

【数47】

εは、１０^－６程度の微小な定数、γは、０．９程度の定数である。μ^{（ｃｙｃｌｅ）}は、慣性項の強さを支配する定数で、交互反復計算回数ｃｙｃｌｅに応じて値を調整する。η_Ｍａｘは学習率の上限値であり、０．１程度の値とするのがよい。上付き括弧文字ｌ＝１，…，Ｌは、分離行列更新計算における繰り返しインデックスである。変数Λ_ｎ ^（ｌ）（ｆ）とΓ_ｎ ^（ｌ）（ｆ）は、それぞれタイムステップｌにおける分離行列更新量の二乗和と損失関数（目的関数の勾配）の二乗和の累積値であり、いずれも実変数である。式（４７－１）と式（４７－４）で、目的関数の勾配の二乗と分離行列更新量の二乗について移動平均をとる。式（４７－２）で周波数ビン・源信号ごとに、学習率η_ｎ（ｆ）を目的関数の勾配と１つ前のタイムステップまでの分離行列の更新量に応じて自動調整する。式（４７－２）では学習率の上限値を０．１としているが、適用ケースに応じてこの上限値は任意に定めることができる。式（４７－３）では、右辺第１項の慣性項で局所解からの脱却を図ることで、目的関数が複雑化した時の解探索を可能にしている。さらに右辺第２項で移動平均により学習率を調整することで、目的関数降下方向への分離行列更新量が急激に変化するのを抑えている。分離行列更新後は式（４７－６）で分離信号を算定する。

【0106】

分離行列の更新演算が収束したら、次に源信号の統計的分散を更新する。源信号の統計的分散σ_ｎ（ｆ，τ）の周波数スペクトルをｔ_ｂ（ｆ）により近似し、その時変な励起パターンをγ_ｂ（τ）を用いて再現する。源信号間のスペクトルのランクが大きく異なることを想定し、分割関数β_ｎｂを用いてＢ本の基底をＮ個の源信号に振り分けて、源信号分散を低ランク近似する。分離信号と源信号の統計的分散の誤差評価指標に板倉斎藤擬距離を用いた非負値行列因子分解（ＩＳＮＭＦ）（非特許文献１６）による分散の更新は、補助関数法に基づく次の乗算型更新式（非特許文献１７）により行われる。

【0107】

【数48】

ｂ＝１，…，Ｂは、低ランク近似で用いる基底のインデックスである。この更新式は、源信号の統計的分散が低ランク分解できるという拘束条件を満たしつつ、式（３７）の目的関数を最小化する分散を求める最尤推定と等価である。

【0108】

式（４１）で求まる源信号空間情報に基づく正則化項の勾配と式（３５）又は式（３６）で求めた源信号分散の初期値を用いて、式（３９）により目的関数の勾配を求める。源信号空間情報から式（３３）により求まる分離行列を初期値として、式（４７）により分離行列と分離信号を更新する。更新した分離行列と分離信号を用いて、式（３９）と式（４１）による目的関数勾配の算定と式（４７）による分離行列更新を繰り返す。繰り返しインデックスｌがＬ回に達したら分離行列の更新を終了し、源信号の統計的分散更新に移る。式（４８）を用いて源信号の統計的分散を更新する。統計的分散は、初回更新時（ｃｙｃｌｅ＝１）には３回程度反復更新し、それ以降（ｃｙｃｌｅ≧２）は、１回ずつ更新する。分離行列と統計的分散の更新という２段階の更新を交互に繰り返すことで、信号分離を行う。交互計算を何回繰り返すか事前に決めておいて、繰り返し計算が進むにつれて正則化項の値を小さくしていく。

【0109】

【数49】

例えば式（４９）のように、交互反復回数ｃｙｃｌｅが予め設定した反復回数ＣＹＣＬＥに０．５～０．８程度の定数Ｃｒを掛けた回数を超えると正則化項と慣性項がゼロになるように、各項の係数を単調減少させる。このように正則化項と慣性項を繰り返し計算が進むほど、小さくしていくことで、源信号の空間情報に基づいて、統計的分散を学習させることができる。正則化項の初期値ρ_Ｄ ^（ｌ）は、１０～３０程度の定数であり、源信号が互いに近接して存在する場合などの混合行列が悪条件の時には小さな値、源信号到来方向の信頼度が高い時に大きな値とするとよい。

【0110】

これまで説明したＬ_ＭＮＭＦに源信号空間事前情報に関する制約条件ＬＳ（Ｗ（ｆ））を加えた目的関数ＬＬ（Ｗ（ｆ））を最小化する分離行列の算定フローを以下にまとめる。
ａ．初期値設定
繰り返しインデックス（上付括弧文字）をｃｙｃｌｅ＝１，ｌ＝１と初期化する。各インデックスの更新回数ＣＹＣＬＥ及びＬを予め決めておく。式（３２）～式（３６）により、それぞれ分離行列と源信号分散の初期値を設定する。分離行列更新パラメータを次の式で初期化する。

ｂ．目的関数の勾配算定
分離信号と源信号分散及び式（４１）で求まる源信号空間情報に基づく正則化項の勾配を用いて、式（３９）により目的関数の勾配を求める。
ｃ．分離行列の更新
式（４７）で分離行列と分離信号を更新する。繰り返しインデックスｌがＬ回に達していなければ、ｌを１つ増やしてステップｂに戻る。ｌがＬ回に達したら、次の式で分離行列更新の際に用いる変数の初期値を設定し、繰り返しインデックスｌを１に初期化して源信号分散の更新に移る。

ｄ．源信号分散の更新
式（４８）で源信号分散を更新する。初回更新時（ｃｙｃｌｅ＝１）には３回程度反復更新し、それ以降の更新時（ｃｙｃｌｅ≧２）には、１回ずつ更新する。
ｅ．収束判定
繰り返しインデックスｃｙｃｌｅを１つ増やす。この時、ｃｙｃｌｅ≦ＣＹＣＬＥならばステップｂに戻る。

【0111】

ステップｃの分離行列更新に際して、交互反復回数ｃｙｃｌｅが総反復回数ＣＹＣＬＥの半分を超えて、正則化項が０になった後は、式（３７）に示す目的関数が通常のＭＮＭＦと同じになり、その解空間が単純形状となる。そのため、正則化項が０になった後は、反復射影法（非特許文献１８）を適用して高速かつ安定に分離行列を更新することも可能となる。このような場合でも正則化項が０になるまでの間に請求項３に記載の正則化項を目的関数に加えている場合や、請求項４に記載の分離行列初期値を用いている場合、請求項５に記載の分離行列更新方法を適用している場合は、本特許の請求範囲とみなす。

【0112】

〔基底ベクトル算定部〕
源信号到来方向が既知でない場合、上記分離信号算定部により算定される分離信号の複素ベクトルｙ（ｆ，τ）の振幅は一意に定まらず、不定性が残る。そのため、全ての周波数成分から逆フーリエ変換により時間領域の信号を再構成する際に、周波数ひずみが生じる。そこで分離行列の逆行列をＷ（ｆ）に乗算することで、Ｍ’番目のセンサ対応する分離行列Ｗ^（Ｓ）（ｆ）を求める。

【0113】

【数50】

上記分離信号算定部で求めた分離行列を用いて式（５０）により、分離信号を算定する。

【0114】

前項で算定した分離信号には、ｎ番目の源信号インデックスを有する分離信号が必ずしもｎ番目の源信号に対応しないという交換の不定性という課題が存在する。周波数ビンごとの分離信号における交換の不定性を解決する方法としては、源信号の特徴量を用いて基底ベクトルを構成する手法と、分離行列に含まれる源信号到来方向の情報を用いて基底ベクトルを構成する手法とが考えられる。この内、後者の手法のみを示す。

【0115】

【数51】

周波数伝達関数の基準センサにおける値との比を分離行列の逆行列を用いることで、式（５１）により近似する。

【0116】

【数52】

Ｍ’は基準センサの番号、ｄ_ｍａｘは基準センサと他のセンサの距離の最大値である。式（５２）は、フォンミーゼス・フィッシャー分布を仮定したクラスタリング手法に適用できる基底ベクトルｖ_ｎ（ｆ）の算定式であり、分離行列の逆行列から周波数伝達関数を算定して導くことができる。上記、基底ベクトルは周波数（波長）に依存せず、センサと源信号の相対位置によって値が変化する構成になっているため、周波数ごとの交換不定性解決のためのクラスタリングに用いることができる。又、ノルムが１となるように正規化を行っているため、球面上の分布となり、フォンミーゼス・フィッシャー分布を仮定したクラスタリング手法が適用可能となる。

【0117】

〔クラスタリング部〕
第ｃクラスタにおける未知パラメータの集まりθ_ｃ（以降、モデルパラメータと称する）が既知の基での、Ｍ次元基底ベクトルの条件付確率Ｐｒ（ｖ_ｎ（ｆ）｜θ_ｃ）を考える。この時、ベイズの定理より、第ｃクラスタの重要度を表す混合重みπ_ｃと基底ベクトルの信頼度を表す重みα_ｎ（ｆ）を用いて、基底ベクトルが与えられた条件下でのモデルパラメータの条件付確率の重み付き対数尤度ＬＣは次の式で定義される。

【0118】

【数53】

ｃ＝１，…，Ｃはクラスタインデックスである。対数尤度の下限はＪｅｎｓｅｎの不等式を用いると、次のように与えられる。

【0119】

【数54】

ｕ_ｃｎ（ｆ）は、Ｍ次元基底ベクトルのクラスタへの帰属度合を表す実変数であり、帰属度と称する。帰属度は、周波数ビンｆ、信号ｎごとの和が常に１となる。式（５４）の右辺は、重み付き対数尤度ＬＣの下限を与えるため、右辺を最大化して下限を持ち上げれば、対数尤度ＬＣを最大化する変数の組み合わせが求まる。

【0120】

以降、混合モデルとしてフォンミーゼス・フィッシャー分布を仮定して、式展開する。この時、確率密度関数は次の式で与えられる。

【0121】

【数55】

この時、θ_ｃは、次の式で与えられる。

【0122】

【数56】

ここに、μ_ｃは各クラスタにおける基底ベクトルの平均方向を表すＭ次元ベクトルであり以降、平均方向と称する。κ_ｃは、クラスタごとの基底ベクトルの方向集中度を表す実定数であり以降、方向集中度と称する。変数Ｉ_{ｏ．ｒＭ－１}（κ_ｃ）は、（０．５Ｍ－１）階の第１種変形ベッセル関数である。

【0123】

混合重みと帰属度の変域を考慮しつつ、平均方向が規格化条件を満たすように制約を与えつつ、式（５５）の右辺を最大化するために、次のような最適条件を定義する。

【0124】

【数57】

【0125】

式（５７）の目的関数に３種類の等式制約をラグランジュ変数により課した後、帰属度ｕ_ｃｎ（ｆ）で偏微分して極値を求めることで、式（５８）を導出できる。

【0126】

【数58】

同様にパラメータ集合θ_ｃに含まれる変数の内、π_ｃ及びμ_ｃで偏微分して極値を求めることで、各変数を算定できる。

【0127】

【数59】

【0128】

基底ベクトルの方向集中度は、第１種変形ベッセル関数の中に入り込んでいるため、算定式として値を得ることはできない。そこで、階数の一つ異なる第１種変形ベッセル関数の比率を定義する。

【0129】

【数60】

この比率を連分数に展開した後、一次近似することで、方向集中度κ_ｃを算定する。

【0130】

【数61】

式（６１－１）の近似的に求めた比率と、式（６２－１）の実際に基底ベクトルの方向集中度を第１種変形ベッセル関数に代入して求めた比率との誤差が最小となるように、ニュートン法に基づく解の更新を行う。具体的には、式（６２－１）と式（６２－２）を交互に適用して、最適な基底ベクトルの方向集中度を探索する。更新は２回程度行えばよい。

【0131】

【数62】

基底ベクトルの方向集中度は、次の式変換により、全てのクラスタで同一の値と仮定することもできる。この時、変数に対してクラスタインデックスｃを削除する置換処理を施すことにより、基底ベクトルとその重心ベクトルとの距離が近い、すなわち内積μ_ｃ ^Ｈｖ_ｎ（ｆ）が大きいほど、そのクラスタへの帰属度が大きくなるという単純な関係性が導けるため、基底ベクトルの方向集中度の推定誤差の影響を受けにくくなり、数値計算の安定性が向上する。一方で、源信号数が観測点数を超えている場合や観測ノイズが含まれる場合などの方向集中度がクラスタ間で異なる場合や、基底ベクトルの重みが互いに異なる場合は正確な解が導けない。しかし、近似的にクラスタの位置を求めたい場合には、数値安定性に優れるこの仮定は便利である。よって、源信号数とセンサ数が同数の場合において、初期クラスタ決定時には、一時的に基底ベクトルの重みを全て１とした上で、方向集中度が、全てのクラスタで同一とする仮定を用いることを推奨する。源信号数がセンサ数よりも多い場合は、クラスタの規模が互いに大きく異なる可能性が高いので、初期クラスタ決定時においてもクラスタごとに方向集中度を算定することを推奨する。方向集中度が、全てのクラスタで同一とする仮定を用いた場合、方向集中度は全クラスタの平均値に相当する次の式で近似される。

【0132】

【数63】

【0133】

基底重みα_ｎ（ｆ）は、基底ベクトルの信頼度に基づいて算定される値で様々な求め方が考えられるが、後程、式（７１）で定義する連続値マスクＭ_ｎ（ｆ）をそのまま基底重みとして、用いるのが最も簡便と考えられる。

【0134】

源信号の数が未知の時、最適なクラスタ数を求めることで、源信号の数を推定する。最適なクラスタ数は、次の手順で決定する。まずクラスタ数Ｃを何パターンか設定して式（５８）～式（６３）の繰り返し計算によりクラスタ分析する。その結果を式（６４）により評価し、評価指標ＥＶを最小とするクラスタ数が最適であるとする。

【0135】

【数64】

ρは定数で２又は３程度の値を用いる。

【0136】

式（５８）～式（６３）の更新式で得られる帰属度とモデルパラメータの算定結果を基に、同じ周波数ビンの異なる信号番号を有する基底ベクトル（分離信号）は同じクラスタに属さないという制約を加えながら、クラスタ分割を詳細に行う。従来、このクラスタ分割は評価値を設定し、その値が最良となる基底ベクトルの交換の組み合わせを総当たりで探索していた。しかし、この方法では、源信号の数が増えた時に、探索の組み合わせが多くなり、計算に時間がかかるという課題がある。そこで、最適解探索の範囲を限定することで、高速にクラスタリングする方法を提案する。源信号の数が多くなった時に、クラスタが局所解に収束するのを防ぐために、帰属度とモデルパラメータの初期値として式（５８）～式（６３）の更新式で得られる値を用いる。

【0137】

独立成分分析等による信号分離では、同じ周波数ビンの異なる信号番号を有する分離行列は同じクラスタに属さないという制約を仮定できるため、クラスタ数Ｃは常に源信号数Ｎ以上となる。又、帰属度を０又は１と２値化すると対数尤度下限値は帰属度が０の時に０となり、帰属度が１の時に値を持つ。そこで、帰属度が１の時の対数尤度下限値を式（６５）により、部分対数尤度と定義する。部分対数尤度は、帰属度を補助変数として与える前の対数尤度と同じであり、各周波数ビンにおけるクラスタ・信号ごとの３次元テンソル構造となる。

【0138】

【数65】

確率分布がフォンミーゼス・フィッシャーズ分布に従う時、部分対数尤度は式（６６）となる。

【0139】

【数66】

フォンミーゼス・フィッシャーズ分布において、帰属度が０か１で表され、基底ベクトルの方向集中度がクラスタ間において一定で、基底ベクトルの信頼度が全て同一のため、基底重みα_ｎ（ｆ）が全て１とできるという特定の条件が揃った時、部分対数尤度は以下により、簡易的に算定することもできる。

【0140】

【数67】

周波数ビンごとに、クラスタ番号を行成分、信号番号を列成分にもつ次の行列を考える。周波数ビンを固定して、クラスタ・信号番号の２次元テンソル構造に着目すると、クラスタ・信号番号が重複しないように、選んだ部分対数尤度のＮ個のテンソル和が最大となる組み合わせが、帰属度を０又は１と２値化した時の、対数尤度下限の最大値を与える。そのため、対数尤度を最大化する最適解を探索するには、部分対数尤度において、周波数ビンごとに、クラスタ・信号番号の２次元テンソル構造を考え、行と列が重複しないように選んだテンソルの組み合わせからその和が最大となる組み合わせをみつければよい。この時、和が最大となる部分対数尤度のテンソルの組み合わせの候補は、周波数ビンごとに値が大きい順に並び替えた時、その上位Ｎ－１位のテンソルが含まれている組み合わせから選べばよいことは明らかである。このことを念頭において、最適な組み合わせを部分対数尤度のテンソル及びその行と列番号の集合を用いて探索する手法を以下に示す。

【0141】

【数68】

部分対数尤度、クラスタ番号、信号番号に関する３種類の集合を定義し、それぞれ上付き文字としてＬ，Ｃ，Ｎをつける。
部分対数尤度Ｌ_ｃｎ（ｆ）を周波数ビンごとに大きい順に並べた集合をＢ^Ｌ（ｆ）、対応するクラスタ（式（６８）の行）と源信号（式（６８）の列）の番号の集合をそれぞれＢ^Ｃ（ｆ）及びＢ^Ｎ（ｆ）と定義し、それらの和集合をＢ（ｆ）とおく。又、周波数ビンごとに、部分対数尤度、クラスタ番号、源信号番号の各集合要素の内、ｉ番目の要素を集めた集合をＢ_ｉ（ｆ）と定義する。以降、使用するその他の同じ集合に対しても同様の表記規則を用いる。

【0142】

【数69】

ここで、上付文字ｌ＝１～Ｌは繰り返しインデックスである。又、ｉ＝１，…，Ｎは集合要素の位置を表す。
式（６９）の集合構造を以降、Ｂ（ｆ）｛ｂ（ｆ）｝と短縮表記し、その他の集合にも同じ表記を適用する。最適解の条件を満たす要素すなわち部分対数尤度とそのクラスタ、信号番号を格納する集合をΩ^（ｌ）｛ω^（ｌ）｝、条件を満たさない要素を格納する集合をΨ｛ψ｝を定義する。アルゴリズム内の記号φは空

【0143】

・アルゴリズム１
０．周波数インデックスを初期化する

１．繰り返しインデックスと最適解候補格納用集合を初期化し、探索開始解を登録する。

２．最適解でない要素の格納用集合Ψを初期化する。

３．集合Ｂ（ｆ）の要素に対して、最適解条件を満たすものと満たさないものに分ける。

４．最適解条件を満たさない集合Ψ内の要素が過去の繰り返し計算で最適解集合Ωに含まれているか調べ、含まれている時は集合Ψから該当要素を削除する。

５．集合Ψが空集合ならば、ステップ７へ、そうでなければ、ステップ６へ進む。
６．空集合でない集合Ψの要素の内、要素番号が最小の要素が最適解条件を満たす集合の最初の要素Ω_１に保存して、ステップ２へ戻る。

ステップ２へ戻る。
７．最適解条件を満たす集合Ωの内、評価値の最も高い集合を最適解として帰属度に反映する。

８．ｆ＜Ｆならステップ１へ戻る。ｆ＝Ｆなら繰り返し計算を終了する。

【0144】

帰属度を０又は１と２値化した時に、対数尤度を最大化する最適なクラスタを生成する計算フローを以下にまとめる。初めに式（５８）の帰属度の算定値を用いて、式（５９）～式（６３）を解くと、事後確率最大化条件を満たすモデルパラメータθ_ｃが求まる。モデルパラメータを用いて、部分対数尤度を算定し、対数尤度を最大化する部分対数尤度の組み合わせをアルゴリズム１により、探索する。この最適な組み合わせを帰属度に反映する。この計算を繰り返すことで、対数尤度を最大化する最適なクラスタの組み合わせが特定できる。具体的には、次の繰り返し計算を行う。

【0145】

ｉ） 変数集合θ_Ｃの算定
帰属度が与えられた条件下で、対数尤度を増大させる変数集合θ_Ｃを式（５８）～式（６３）により更新して、算定する。
ｉｉ）部分対数尤度の算定
変数集合θ_Ｃより、各クラスタにおける周波数ビン・信号ごとの部分対数尤度Ｌ_ｃｉ（ｆ）を式（６６）又は式（６７）により算定する。
ｉｉｉ）帰属度の決定
部分対数尤度の算定結果に基づいて、対数尤度を最大化するように基底ベクトルをクラスタへ割り振る組み合わせをアルゴリズム１のステップ２～６により探索し、最適な割り振りの組み合わせから、アルゴリズム１のステップ７により帰属度を０又は１で求める。
ｉｖ）繰り返し計算過程において、０又は１で定義された帰属度又は平均方向の変化がなくなるまでｉ）～ｉｉｉ）の計算を繰り返す。

【0146】

最適化により算定した分離行列には、ノイズ等による源信号波形の推定誤差や学習そのものがうまく進まなかったことによる誤差等の様々な推定誤差が含まれる可能性がある。そこで、分離行列の推定精度を評価し、信頼度の低い分離行列から算定された分離信号の影響を減らすために、連続値マスクＭ_ｎ（ｆ）を定義する。連続値マスクの算定方法は、信号の分離環境やその目的に応じて異なるため、特に指定しないが、方向別に主となる信号データを得る時に有効な算定手順を例として示す。正しく算定された分離行列の逆行列は、源信号ごとの周波数伝達関数を列ベクトルに持つ混合行列と対応する。混合行列の列ベクトルには、白色化空間において互いに直交するという性質がある。この性質を利用して、分離行列が正しく算定されたかを判定する。初めに、分離行列の逆行列を白色化行列Ｖ（ｆ）により、白色化空間に移行したＮ×Ｎ次元の行列Ｇ（ｆ）を各周波数ビンにおいて算定する。白色化行列は、観測信号の相関行列Ｅ［ｘ（ｆ，τ）ｘ^Ｈ（ｆ，τ）］の固有値と固有ベクトルから算定され、一般的に観測信号を白色化すなわち、全周波数帯にわたって一定とし、信号間で無相関化する時に用いられる。

【0147】

【数70】

行列Ｇ（ｆ）の列ベクトルｇ_ｎ（ｆ）が互いになす角度をその内積の大きさにより評価する。評価は周波数ビンごとに行い、同一周波数ビンの他列ベクトルとの内積の２乗和の平均値が一定値を超えた源信号インデックスを異常とみなし、異常度合に応じてマスクを０に近づける。それ以外の正常値については、連続値マスクを１に近づける。この判定をシグモイド関数により連続値として行うと、次の式により連続値マスクが求まる。

【0148】

【数71】

Ｋｉｊｕｎは、列ベクトルの内積の和が異常と判定する基準値であり、０～１の値をとる。Ｋｏｕｂａｉは、正の実数で、シグモイド関数が１から０に切り替わる速さを表し、値が大きいほど０か１のどちらかの値のみを出力する傾向が強まる。

【0149】

【数72】

連続値マスクの値が閾値Ｔｈｒｅｓｈｏｌｄを下回った分離信号を除外するために、バイナリマスクＭ_ｎ ^（Ｂ）（ｆ）を定義する。

【0150】

分離信号に０又は１で定義された帰属度を掛けてクラスタごとに加算し、１つの源信号を１つの行ベクトルにもつ分離信号ｙ_ｃ ^{（ＩＣＡ）}（ｆ，τ）を算定する。

【0151】

【数73】

以上のクラスタリング方法により、効率的に交換の不定性を解決できる。

【0152】

分離信号ｙ_ｃ ^{（ＩＣＡ）}（ｆ，τ）には、算定精度の低い分離信号が含まれており、このまま出力すると源信号推定結果の誤差が大きくなる。そこで源信号到来方向算定方法についてまとめた後、算定精度の低い分離信号の影響を低減させる方法をこの後、分離信号補正部で演算式としてまとめる。

【0153】

〔源信号到来方向算定部〕
まず、源信号到来方向ベクトルｑ＾_ｎ（ｆ）を、観測点座標ベクトルのＭｏｏｒｅ－Ｐｅｎｒｏｓｅ疑似逆行列Ｐ^＋を用いて、推定する（非特許文献１２）。式（７４）において、基準センサの位置が源信号到来方向ベクトルの原点の位置すなわち視点となる。

【0154】

【数74】

行列Ｐ及びベクトルｒ＾_ｎ（ｆ）はそれぞれ、観測点番号として基準センサ以外の（Ｍ－１）通りの位置ベクトルに対して、基準センサを起点とする観測点座標ベクトル及び周波数伝達関数の偏角ａｒｇを羅列したものである。

【0155】

【数75】

源信号及び周波数ビン番号に対して源信号到来方向推定値を加算することでクラスタごとの信号方向ベクトルＱ＾_ｃを式（７５）により算定する。この源信号到来方向推定値に対してカメラにより事前に撮影した画像を重ね合わせることで、信号到来方向の推定精度を向上させることもできる。
Ｍ_ｎ（ｆ）は、連続値マスクであり、周波数ビン・源信号ごとの分離行列の信頼度に基づいて設定する。上記算定部で掛けている連続値マスクは、背景ノイズのほとんどない場合など分離行列の信頼度が高い場合などは、省略することもできる。Ｗｈｅｉｇｈｔ ｆ）は周波数ビン・源信号ごとの源信号到来方向ベクトルの重みであり、どの周波数も同じ信頼度なら全て１となる。重みとして周波数ごとの観測値のパワースペクトルを指標とする時は式（７６）で算定される。他にも重みとしては、信号部分空間の固有値の和を指標とすることもでき、近似的に直接、観測点まで届いた到来波のパワーの総和をとっていることになる。

【0156】

【数76】

ＰＳＤ（ｆ）は、全センサのパワースペクトルを時間フレーム方向に加算した値である。
〔分離信号補正部〕

【0157】

【数77】

クラスタ番号ｋを混合行列の列番号すなわち源信号番号へ変換する順列Π_ｆ（ｋ）を式（７７）で定義する。この順列は、周波数ビンｆにクラスタ番号ｋと対応する源信号が存在しない時に０を返すため、源信号の存在を調べる指標として用いることもできる。周波数ビンｆに存在する源信号の数Ｎ^{（Ａｃｔ）}（ｆ）は、式（７８）で帰属度を源信号番号及びクラスタ番号の方向へ加算することで求まる。

【0158】

【数78】

式（７５）で得られる源信号到来方向推定結果を用いて、Ｍ×Ｎ^{（Ａｃｔ）}（ｆ）次元の混合行列Ａ＾（ｆ）を算定する。

【0159】

【数79】

この混合行列の逆行列を求めることで、簡易的に分離行列と分離信号が算定できる。

【0160】

【数80】

ｅ_ｎ（ｆ）はｎ番目の要素が１、他要素が０のＮ^{（Ａｃｔ）}（ｆ）次元単位ベクトルである。式（８０－３）の上付括弧文字＋は、ムーアペンローズの擬似逆行列を表す。源信号到来方向より得られる分離信号ｙ_ｃ ^{（ＤＯＡ）}（ｆ，τ）をベクトル形式でまとめると、

【0161】

【数81】

【0162】

式（７７）～式（８１）に示す源信号到来方向から分離信号を算定する演算は、クラスタ数▲Ｃ▼と源信号数Ｎが同じ数の時、次の式で簡潔にまとめることができる。

【0163】

【数82】

【0164】

式（７３）に示すＩＣＡによる分離信号算定結果と源信号到来方向に基づく分離信号算定結果を用いて、次式により、算定精度の低い分離信号の影響を軽減した分離信号ｙ^{（ＭＯＤ）}（ｆ，τ）を算定する。

【0165】

【数83】

この式ではバイナリマスクを用いて、ＩＣＡにより算定される分離信号ｙ^{（ＩＣＡ）}（ｆ，τ）の内、算定精度の低い分離信号を取り除き、代わりに源信号到来方向に基づく分離信号ｙ^{（ＤＯＡ）}（ｆ，τ）を加算して、取除分を埋め合わせている。

【0166】

源信号空間事前情報が得られている場合は式（５０）の分離信号、得られていない場合は式（８３）の分離信号を逆フーリエ変換し、窓関数を掛けて重ね合わせることで、分離信号の時間波形が得られる。

【0167】

本発明は、以上の発明の実施の形態に限定されることなく、特許請求の範囲に記載された発明の範囲内で、種々の変更が可能であり、それらも本発明の範囲内に包含される。

【0168】

本特許文面では、分離行列として、独立成分分析及び独立低ランク行列分析を用いた場合を記載しているが、その他の分離行列を用いた場合でも、本発明におけるクラスタリング方法又は源信号到来方向推定方法を用いた場合には、本発明の範囲に属するとみなす。

【0169】

源信号波形の方法として、その源信号の統計的分散情報を人工知能に事前学習させ、高精度に信号分離を行うための手法として、独立深層学習行列分析が存在する。この方法では、分離行列を補助関数法に基づく反復射影法を用いて、更新している。しかし、独立深層学習行列分析には、事前に源信号の情報を一定レベルで得られないと、源信号の統計的分散情報を正しく学習できないという課題がある。源信号到来方向の情報を式（４１）に示す制約条件として加えた式（３７）の事前情報制約付加に基づく分離行列更新方法を、独立深層学習行列分析における反復射影法の代わりに用いることで、分離信号の事前情報が少ない状況でも高精度に源信号波形を推定することができるようになる。そのため独立深層学習行列分析（非特許文献１９）においての適用も考えられる。このような場合でも、目的関数に源信号到来方向を事前情報とした制約条件を付加する場合において、当該制約条件が「同一の源信号番号を有する分離行列の行ベクトルと混合行列の列ベクトルの内積が１になり、異なる源信号番号を有する分離行列の行ベクトルと混合行列の列ベクトルの内積が０になる」という情報を含む場合など、請求項３に記載の事項と合致する場合は、本発明の範囲に属するとみなす。又、１つ前のタイムステップにおける分離行列の更新量に定数を掛けた値を慣性項として用いている場合や、行列ブロックごとに勾配降下方向への更新量すなわち学習率を前回までの分離行列の更新量や目的関数の勾配に応じて調整している場合など、請求項５に記載の事項と合致する場合は、本発明の範囲に属するとみなす。

【0170】

源信号到来方向等の事前情報を制約条件とした分離行列更新方法を独立深層学習行列分析に応用する場合においても、演算の流れは、独立低ランク行列分析と同様であり、信号の統計的分散の更新方法が異なるだけである。基本的な演算の流れを以下に示す。

【0171】

【数84】

分離信号のパワースペクトル構造に関する事前情報が得られている場合、ニューラルネットワークにより信号波形の特徴量から統計的分散を上式により推定することができる。ニューラルネットワークとして独立深層学習行列分析を用いた場合の具体的な分散推定の演算方法は、非特許文献１９を参照されたい。ニューラルネットワークによる分散推定の方法には、様々な手法があるが、いずれの手法を用いた場合でも制約条件の設定方法や分離行列の更新方法は変わらない。まず式（３７）に示す源信号到来方向に関する事前情報制約を付加した目的関数を用い、式（４７）のの分離行列更新ステップを数回、繰り返して分離行列を更新する。更新が完了した分離行列を用いて分離信号を算定し、ニューラルネットワークに入力する。ニューラルネットワークで分離信号波形の入力データと過去に学習した分離信号のパワースペクトル構造を比較し、最も適合する構造を選択してその分散を時間周波数点ごとに求める。得られた分散を用いて、再び上記計算ステップで分離行列を更新し、分離信号を算定する。この演算フローを繰り返すことで、分離信号を正確に推定することができる。この場合、本特許における請求項３に記載の分離行列更新方法と請求項５に制約条件の記載を用いていることになる。これらのいずれかを用いている場合は、本特許の請求範囲とする。

【0172】

図－５及び図－６の計算フローによるクラスタリング計算において、周波数ビンｆ・源信号ｎごとに基底ベクトルを生成しているが、例えば調波信号の分離において、基底ベクトルの生成時に、周波数ビンｆを時間周波数点に置き換えてクラスタ分析を実施することも可能である。このような場合でもクラスタ分析に、フォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５７）に示す最適性条件の基でパラメータを更新している場合など、請求項６、請求項７又は請求項８に記載のクラスタリング方法と合致する場合は、本発明の範囲に属するとみなす。

【0173】

請求項６～請求項８に記載の「クラスタリング部」で提案しているクラスタ分析計算は、ミュージック法（非特許文献７）による周波数ビンごとの源信号到来方向推定結果を解析対称周波数全域にわたって、クラスタリングする際に用いることもできる。その適用手法を以下に示す。

【0174】

【数85】

水平角θと仰角Φを変えながら、式（８５）により源信号到来方向と対応する周波数伝達関数を算定する。式（８６）で算定されるミュージックスペクトルが局所的にピークになる時の水平角θと仰角Φを源信号到来方向として、源信号到来方向に対応する周波数ビンごとの基底ベクトルを式（１３）により、源信号ごとに算定する。

【0175】

【数86】

Ｑ_Ｎ（ｆ）は、観測信号の相関行列の最小の固有値に対応した固有ベクトルである。周波数ビンごとに、源信号到来方向別に算定された基底ベクトルに対して、本特許のクラスタリング部に記載の演算（図－５及び図－６）を適用すると、解析対称周波数全域にわたって、源信号インデックスｎを対応付けることが可能となる。このような場合でも、クラスタ分析に、フォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５８）～式（６３）の最適性条件の基でパラメータを更新している場合など、請求項６、請求項７又は請求項８に記載のクラスタリング方法と合致する場合は、本発明の範囲に属するとみなす。

【0176】

請求項６～請求項８に記載の「クラスタリング部」で提案しているクラスタ分析計算は、ミュージック法を拡張したＫＲ積－圧縮ミュージック法（非特許文献２０）による周波数ビンごとの源信号到来方向推定結果を解析対称周波数全域にわたって、クラスタリングする際に、用いることもできる。その適用手法を以下に示す。

【0177】

観測信号の時間フレームごとのＭ×Ｍ次元の相関行列を用いて、Ｎ^２×Ｔ次元のデータ行列Ｙ^～（ｆ）を次の式で周波数ビンごとに定義する。

【0178】

【数87】

Ｖｅｃ（・）は行列のベクトル化を表す。データ行列の特異値分解を行い、特異値０に対応する左特異ベクトルから成る行列をＱ^～（ｆ）と定義する。この時、一般化ミュージックスペクトルＳｐｅｃ（ｆ，θ，Φ）は次の式で求まる。

【0179】

【数88】

番号の列ベクトルにクロネッカー積を適用する演算であり、一般にＫＲ積と呼ばれる。推定源信号到来方向集合Ｃは、データ行列に対して、事前に直交マッチング追跡（非特許文献２１）などの圧縮センシングのアルゴリズムを適用することで、得られるＮ－Ｔ個の源信号到来方向推定値の集合である

【0180】

水平角θと仰角Φを変えながら、式（８５）により源信号到来方向と対応する周波数伝達関数を算定する。式（８８）で算定されるミュージックスペクトルが局所的にピークになる時の水平角θと仰角Φを源信号到来方向として、源信号到来方向に対応する周波数ビンごとの基底ベクトルを式（１３）により、源信号ごとに算定する。周波数ビンごとに、源信号到来方向別に算定された基底ベクトルに対して、本特許のクラスタリング部に記載の演算（式（５８）～式（６３））を適用すると、解析対称周波数全域にわたって、源信号インデックスｎを対応付けることが可能となる。このような場合でも、クラスタ分析に、フォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５７）に示す最適性条件の基でパラメータを更新している場合など、請求項６、請求項７又は請求項８に記載のクラスタリング方法と合致する場合は、本発明の範囲に属するとみなす。

【0181】

この他にも直交マッチング追跡などの圧縮センシングのアルゴリズムを単体で用いた場合や、ＫＲミュージック法（非特許文献２２）、圧縮ミュージック法（非特許文献２３）を用いた場合にも、周波数ビンごとに、源信号到来方向が得られ、本特許で発明したクラスタリング演算を用いて、解析対称周波数全域にわたって、源信号インデックスｎを対応付けることが可能となる。このような場合でも、クラスタ分析に、フォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５７）に示す最適性条件の基でパラメータを更新している場合など、請求項６、請求項７又は請求項８に記載のクラスタリング方法と合致する場合は、本発明の範囲に属するとみなす。

【0182】

式（５８）～式（６３）のクラスタリング計算は、調波構造を仮定した音声信号の分離に適用することもできる。以下、その適用手法を示す。

【0183】

【数89】

式（８９）の観測信号の時間フレームτと周波数ビンｆを時間周波数点のインデックスｋ＝１，…，Ｋを定義することで、同一次元にしてベクトル化する。この時、時間周波数点の各点ｋにおける観測信号ベクトルは、センサ数をＭとすると、次のＭ次元ベクトルで定義される。

【0184】

【数90】

時間周波数点の各点ｋごとに、ノルムが１のＭ次元埋め込みベクトルｖ^～ _ｎ（ｋ）を考える。

【0185】

【数91】

ｎ＝１，…，Ｎは、源信号のインデックスを表し、各時間周波数点において、多くてもＮ個の源信号がアクティブと仮定している。埋め込みベクトルを用いると、時間周波数点の各点ｋにおいて、支配的な源信号波形を観測信号波形から算定することが可能となる。埋め込みベクトルは時間周波数点ごとに独立して生成されるため、異なる時間周波数点において、源信号インデックスｎの一貫性がなくなる。例えば、源信号数が２つの場合を考えると、時間周波数点１と時間周波数点２で、源信号の組み合わせは２パターン考えられる。そこで、全時間周波数点において、埋め込みベクトルをクラスタリングすることで、対応関係を明確化する。クラスタリングは、全源信号数をクラスタ数Ｃと定義し、ｆ→ｋとインデックスを置き換えて、「クラスタリング部」に示す演算方法すなわち式（５８）～式（６３）及びアルゴリズム１を適用することで、実施することができる。よって調波構造を仮定した音声信号の分離を行う場合においても、クラスタ分析に、混合モデルとしてフォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５７）に示す最適性条件の基でパラメータを更新している場合など、請求項６又は請求項７又は請求項８に記載のクラスタリング演算手法と合致する場合は、本発明の範囲に属するとみなす。

【符号の説明】

【0186】

１ 信号取得部
２ 信号波形分離部
３ 計算結果表示部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】