JP4449871B2 - 音声信号分離装置及び方法 - Google Patents

Description

本発明は、複数の信号が混合された音声信号を独立成分分析（Independent Component Analysis；ＩＣＡ）を用いて信号毎に分離する音声信号分離装置及びその方法に関する。

複数の原信号が未知の係数によって線形に混合されているときに、統計的独立性のみを用いて原信号を分離・復元するという独立成分分析（Independent Component Analysis；ＩＣＡ）の手法が信号処理の分野で注目されている。この独立成分分析を応用することで、例えば話者とマイクロホンとが離れた場所にあり、マイクロホンで話者の音声以外の音を拾ってしまうような状況でも、音声信号を分離・復元することが可能となる。

ここで、時間周波数領域の独立成分分析を用いて、複数の信号が混合された音声信号を信号毎に分離する場合について考える。

図１５に示すようにＮ個の音源からそれぞれ異なる音が鳴っており、それらをｎ個のマイクロホンで観測するという状況を想定する。音源が発した音（原信号）がマイクロホンに届くまでには時間遅れや反射などがあるため、ｋ番目（１≦ｋ≦ｎ）のマイクロホンｋで観測される信号（観測信号）ｘ_ｋ（ｔ）は、下記式（１）のように、原信号と伝達関数との畳み込み演算を全音源について総和した式で表される。また、全てのマイクロホンについての観測信号を１つの式で表すと、下記式（２）のようになる。この式（１）、（２）において、ｘ（ｔ）、ｓ（ｔ）はそれぞれ ｘ_ｋ（ｔ）、ｓ_ｋ（ｔ）を要素とする列ベクトルを表し、Ａはａ_ｉｊ（ｔ）を要素とするｎ行Ｎ列の行列を表す。なお、以下ではＮ＝ｎとする。

時間周波数領域の独立成分分析では、上記式（２）のｘ（ｔ）からＡ及びｓ（ｔ）を直接推定するのではなく、ｘ（ｔ）を時間周波数領域の信号に変換し、Ａ及びｓ（ｔ）に対応する信号を時間周波数領域で推定する。以下、その方法について説明する。

信号ベクトルｘ（ｔ）、ｓ（ｔ）を長さＬの窓で短時間フーリエ変換したものをそれぞれＸ（ω，ｔ），Ｓ（ω，ｔ）とし、行列Ａ（ｔ）を同様に短時間フーリエ変換したものをＡ（ω）とすると、時間領域の上記式（２）は時間周波数領域の下記式（３）で表すことができる。但し、ωは周波数binの番号を示し（１≦ω≦Ｍ）、ｔはフレーム番号を示す（１≦ｔ≦Ｔ）。時間周波数領域の独立成分分析では、式（３）のＳ（ω，ｔ）、Ａ（ω）を時間周波数領域で推定することになる。

なお、周波数binの個数は、本来は窓の長さＬと同一であり、各周波数binは、−Ｒ／２からＲ／２まで（Ｒはサンプリング周波数）をＬ等分したそれぞれの周波数成分を表す。但し、負の周波数成分は正の周波数成分の共役複素数であり、Ｘ（−ω）＝conj（Ｘ(ω)）（conj（・）は共役複素数）として求めることができるため、本明細書では０からＲ／２までの非負の周波数成分（周波数binの個数はＬ／２＋１）のみを考え、その周波数成分に１からＭ（Ｍ＝Ｌ／２＋１）までの番号を振っている。

時間周波数領域でＳ（ω，ｔ）、Ａ（ω）を推定するには、先ず、下記式（４）のような式を考える。この式（４）において、Ｙ（ω，ｔ）はｙ_ｋ（ｔ）を長さＬの窓で短時間フーリエ変換したＹ_ｋ（ω，ｔ）を要素とする列ベクトルを表し、Ｗ（ω）はｗ_ｉｊ（ω）を要素とするｎ行ｎ列の行列（分離行列）を表す。

次に、ωを固定してｔを変化させたときにＹ_１（ω，ｔ）〜Ｙ_ｎ（ω，ｔ）が統計的に独立となる（実際には、独立性が最大となる）ようなＷ（ω）を求める。後述のように、時間周波数領域の独立成分分析ではパーミュテーション（permutation）及びスケーリングの不定性があるため、Ｗ（ω）＝Ａ（ω）^−１以外にも解が存在する。統計的に独立となるＹ_１（ω，ｔ）〜Ｙ_ｎ（ω，ｔ）が全てのωについて得られたら、それらを逆フーリエ変換することで、時間領域の分離信号ｙ（ｔ）を得ることができる。

時間周波数領域における従来の独立成分分析の概略を図１６を用いて説明する。ｎ個の音源が発するお互いに独立な原信号をｓ_１〜ｓ_ｎとし、それらを要素とするベクトルをｓとする。マイクロホンで観測される観測信号ｘは、原信号ｓに上記式（２）の畳み込み・混合演算を施したものである。マイクロホンの数ｎが２であるとき、すなわちチャンネル数が２であるときの観測信号ｘの例を図１７（Ａ）に示す。次に、観測信号ｘに対して短時間フーリエ変換を施し、時間周波数領域の信号Ｘを得る。Ｘの要素をＸ_ｋ（ω，ｔ）とすると、Ｘ_ｋ（ω，ｔ）は複素数値をとる。Ｘ_ｋ（ω，ｔ）の絶対値である｜Ｘ_ｋ（ω，ｔ）｜を色の濃淡で表現した図をスペクトログラムという。スペクトログラムの例を図１７（Ｂ）に示す。この図において、横軸は t（フレーム番号）を示し、縦軸はω（周波数bin番号）を示す。なお、以下では時間周波数領域の信号そのもの（絶対値をつける前の信号）も「スペクトログラム」と表現する。続いて、信号Ｘの各周波数binにＷ（ω）を乗算することで、図１７（Ｃ）に示すような分離信号Ｙを得る。そして、分離信号Ｙを逆フーリエ変換することで、図１７（Ｄ）に示すような時間領域の分離信号ｙを得る。

ここで、独立成分分析において、独立性をどのような尺度で表現するか、また、どのようなアルゴリズムで独立性を最大化するかについては、種々のバリエーションが存在する。本明細書では、一例として、独立性をKullback-Leibler情報量（以下、「ＫＬ情報量」という。）で表現し、独立性を最大化するアルゴリズムとして自然勾配法を用いる場合について説明する。

図１８のように、ある周波数binに着目する。Ｙ_ｋ（ω，ｔ）のフレーム番号ｔを１〜Ｔの間で変化させたものをＹ_ｋ（ω）としたとき、分離信号Ｙ_１（ω）〜Ｙ_ｎ（ω）の独立性を表す尺度であるＫＬ情報量Ｉ（Ｙ（ω））を下記式（５）のように定義する。すなわち、各チャンネルについての周波数bin（＝ω）毎のエントロピーＨ（Ｙ_ｋ（ω））の総和から全チャンネルについての周波数bin（＝ω）毎の同時エントロピーＨ（Ｙ（ω））を減算した値をＫＬ情報量Ｉ（Ｙ（ω））と定義する。ｎ＝２のときのＨ（Ｙ_ｋ（ω））とＨ（Ｙ（ω））との関係を図１９に示す。式（５）のうち、Ｈ（Ｙ_ｋ（ω））はエントロピーの定義により下記式（６）の第１項のように書き換えられ、Ｈ（Ｙ（ω））は上記式（４）により式（６）の第２項及び第３項のように展開される。この式（６）において、Ｐ_{Ｙｋ（ω）}（・）はＹ_ｋ（ω，ｔ）の確率密度関数を表し、Ｈ（Ｘ（ω））は観測信号Ｘ（ω）の同時エントロピーを表す。

ＫＬ情報量Ｉ（Ｙ（ω））は、Ｙ_１（ω）〜Ｙ_ｎ（ω）が独立である場合に最小（理想的には０）となる。ここでは、ＫＬ情報量Ｉ（Ｙ（ω））を最小にする分離行列Ｗ（ω）を求めるアルゴリズムとして自然勾配法を用いる。自然勾配法は、Ｉ（Ｙ（ω））を最小化させる方向を下記式（７）で求め、Ｗ（ω）が収束するまで下記式（９）のようにＷ（ω）をその方向に少しずつ変化させるものである。この式（７）において、Ｗ（ω）^ＴはＷ（ω）の転置行列を表す。また、式（９）において、ηは学習係数（正の微小値）を表す。

上記式（７）は上記式（８）のように変形される。この式（８）において、Ｅ_ｔ［・］は時間方向の平均を表す。また、φ（・） は確率密度関数の対数を微分したものであり、スコア関数（又は、「活性化関数」）と称される。スコア関数にはＹ_ｋ（ω）の確率密度関数が含まれているが、ＫＬ情報量の最小値を求めるためには本当の確率密度関数を用いる必要はなく、Ｙ_ｋ（ω）の分布がスーパーガウシアン（super-gaussian）であるかサブガウシアン（sub-gaussian）であるかに応じて、例えば表１に示すような２種類の確率密度関数を切り換えればよいことが知られている。

また、extended infomax法として、表２に示すような２種類の確率密度関数を切り換えるようにしてもよい。

なお、表１、２において、ｈは確率密度関数を−∞〜＋∞の区間で積分した値を１にするための定数である。Ｙ_ｋ（ω）の分布がスーパーガウシアンであるかサブガウシアンであるかは４次のキュムラントκ_４（＝Ｅ_ｔ［Ｙ_ｋ（ω，ｔ）^４］−３Ｅ_ｔ［Ｙ_ｋ（ω，ｔ）^２］^２）の値の正負で決まり、κ_４が正ならスーパーガウシアン、負ならサブガウシアンである。

上記式（８）、（９）を用いた分離処理は、図２０のフローチャートで表される。先ずステップＳ１０１において、周波数bin毎に分離行列Ｗ（ω）を用意し、初期値（例えば単位行列）を代入しておく。次にステップＳ１０２において、全ての周波数binについてのＷ（ω）が収束したか否かを判別し、収束している場合には処理を終了し、収束していない場合にはステップＳ１０３に進む。ステップＳ１０３では、上記式（４）のようなＹ（ω，ｔ）を定義し、ステップＳ１０４では、ＫＬ情報量Ｉ（Ｙ（ω））を最小化させる方向を上記式（８）に従って求める。そしてステップＳ１０５では、上記式（９）に従ってＫＬ情報量Ｉ（Ｙ（ω））を最小化させる方向にＷ（ω）を更新し、ステップＳ１０２に戻る。なお、ステップＳ１０２〜Ｓ１０５の処理は、各周波数binについてＹ（ω）の独立性が十分に高まり、Ｗ（ω）が略々収束するまで繰り返される。

村田昇著，「入門独立成分分析」，東京電気大学出版局 特開２００４−１４５１７２号公報 Mike Davies,"Audio Source Separation", Oxford University press, 2002（http://www.elec.qmul.ac.uk/staffinfo/miked/publications/IMA.ps） Nikolaos Mitianoudis and Mike Davies,"A fixed point solution for convolved audio source separation", IEEE WASPAA01, 2001（http://egnatia.ee.auth.gr/~mitia/pdf/waspaa01.pdf）

ところで、上述した時間周波数領域の独立成分分析では、信号の分離処理を周波数bin毎に行っており、周波数binの間の関係は考慮していない。そのため、分離自体は成功しても、周波数binの間でスケーリング及び分離先の不統一が発生する可能性がある。このうち、スケーリングの不統一については、音源毎に観測信号を推定する方法により解決できる。一方、分離先の不統一とは、例えばω＝１ではＹ_１にＳ_１由来の信号が現れるのに対してω＝２ではＹ_１にＳ_２由来の信号が現れる、というような現象のことであり、パーミュテーション（置換）の問題と呼ばれている。

パーミュテーションが発生している例を図２１に示す。これは、ＷＥＢページ（http://www.ism.ac.jp/~shiro/research/blindsep.html）にある「X_rsm2.wav」というファイルの最初の３２０００サンプルに対してextended infomax法を用いて時間周波数領域で分離を試みた結果である。なお、原信号の一方は“ワン、ツー、スリー”という音声であり、他方は音楽である。上段のスペクトログラムを時間領域の信号へと逆フーリエ変換すると、下段のように、両チャンネルとも両方の信号が混ざった波形となってしまう。このように、周波数bin毎に分離を行うと、観測信号の種類や分離行列Ｗ（ω）の初期値によっては、図２１のような結果になってしまうことが避けられない。

従来、このパーミュテーションの問題を解消するために、後処理により入れ替えを行う方法が知られている。この後処理では、先ず周波数bin毎の分離によって図２１のようなスペクトログラムを得て、その後、何らかの基準に従ってチャンネル間で分離信号の入れ替えを行うことでパーミュテーションの発生していないスペクトログラムを得る。入れ替えの基準としては、（ａ）エンベロープの類似性（非特許文献１を参照）、（ｂ）推定された音源方向（特許文献１の［従来の技術］を参照）、（ｃ）ａとｂとの組合せ（特許文献１を参照）が挙げられる。

しかしながら、上記（ａ）は、周波数binによってはエンベロープの違いが不明瞭なことがあり、そのような場合には入れ替え間違いが発生してしまう。また、入れ替えを１度間違えると、それ以降の周波数binでは全て分離先を間違えてしまうことになる。また、上記（ｂ）は、方向推定の精度に問題があり、さらにマイクロホンの位置情報が必要である。また、両者を組み合わせた上記（ｃ）は、入れ替えの精度は向上しているものの、上記（ｂ）と同様にマイクロホンの位置情報が必要である。また、何れの方法においても、分離と入れ替えという２つのステップを経るため、処理時間が長いという問題がある。処理時間の観点では、分離が完了した時点でパーミュテーションの問題も解消していることが望ましいが、後処理による方法ではそれは難しい。

そこで、非特許文献２，３では、周波数binの間の関係を分離行列Ｗの更新式に反映させる周波数カップリング（frequency coupling）と呼ばれる方法が提案されている。この方法では、下記式（１０）のような確率密度関数と下記式（１１）のような分離行列Ｗの更新式とを用いている（但し、変数の表記法は本明細書と一致させている）。この式（１０）、（１１）において、β_ｋ（ｔ）はＹ_ｋ（ω，ｔ）の各成分の絶対値の間で平均をとった値を表し、β（ｔ）はβ_１（ｔ），・・・，β_ｎ（ｔ）を対角要素とする対角行列を表す。このβ_ｋ（ｔ）の導入により、周波数binの間の関係がΔＷ（ω）に反映される。

しかしながら、上記式（１１）を反復適用して収束させた分離行列Ｗでは、必ずしもパーミュテーションの問題を解消できない。すなわち、パーミュテーション非発生時のＫＬ情報量がパーミュテーション発生時のＫＬ情報量よりも小さくなるという保証がない。実際に上述した「X_rsm2.wav」というファイルの最初の３２０００サンプルに対して上記式（１１）を用いて分離を試みた結果を図２２に示す。図２１と同様に周波数bin毎の分離は成功しており、パーミュテーションも図２１と比較して改善されているものの、依然としてパーミュテーションが発生している。

本発明は、このような従来の実情に鑑みて提案されたものであり、複数の信号が混合された音声信号を独立成分分析を用いて信号毎に分離する際に、分離後の後処理を行うことなくパーミュテーションの問題を解消することが可能な音声信号分離装置及びその方法を提供することを目的とする。

上述した目的を達成するために、本発明に係る音声信号分離装置は、音声信号を含む複数の信号が混合された時間領域の観測信号を独立成分分析を用いて信号毎に分離し、分離信号を生成する音声信号分離装置において、上記時間領域の観測信号を時間周波数領域の観測信号に変換する第１の変換手段と、上記時間周波数領域の観測信号から時間周波数領域の分離信号を生成する分離手段と、上記時間周波数領域の分離信号を時間領域の分離信号に変換する第２の変換手段とを有し、上記分離手段は、上記時間周波数領域の観測信号と初期値が代入された分離行列とから時間周波数領域の分離信号を生成し、この時間周波数領域の分離信号と多次元確率密度関数を用いたスコア関数と上記分離行列とを用いて該分離行列の修正値を計算し、上記修正値を用いて、上記分離行列が略々収束するまで該分離行列を修正し、略々収束した分離行列を用いて上記時間周波数領域の分離信号を生成することを特徴とする。

また、上述した目的を達成するために、本発明に係る音声信号分離方法は、音声信号を含む複数の信号が混合された時間領域の観測信号を独立成分分析を用いて信号毎に分離し、分離信号を生成する音声信号分離方法において、上記時間領域の観測信号を時間周波数領域の観測信号に変換する工程と、上記時間周波数領域の観測信号と初期値が代入された分離行列とから時間周波数領域の分離信号を生成する工程と、この時間周波数領域の分離信号と多次元確率密度関数を用いたスコア関数と上記分離行列とを用いて該分離行列の修正値を計算する工程と、上記修正値を用いて、上記分離行列が略々収束するまで該分離行列を修正する工程と、略々収束した分離行列を用いて生成された時間周波数領域の分離信号を時間領域の分離信号に変換する工程とを有することを特徴とする。

本発明に係る音声信号分離装置及びその方法によれば、音声信号を含む複数の信号が混合された時間領域の観測信号を独立成分分析を用いて信号毎に分離し、分離信号を生成する際に、初期値が代入された分離行列とから時間周波数領域の分離信号を生成し、この時間周波数領域の分離信号と多次元確率密度関数を用いたスコア関数と上記分離行列とを用いて該分離行列の修正値を計算し、上記修正値を用いて、上記分離行列が略々収束するまで該分離行列を修正し、略々収束した分離行列を用いて生成された時間周波数領域の分離信号を時間領域の分離信号に変換することにより、分離後の後処理を行うことなくパーミュテーションの問題を解消することができる。

以下、本発明を適用した具体的な実施の形態について、図面を参照しながら詳細に説明する。この実施の形態は、本発明を、複数の信号が混合された音声信号を独立成分分析を用いて信号毎に分離する音声信号分離装置に適用したものである。特に、本実施の形態における音声信号分離装置は、従来のような１次元確率密度関数を用いて周波数bin毎のエントロピーを計算する代わりに、多次元確率密度関数を用いてスペクトログラム１枚分のエントロピーを計算することにより、分離後の後処理を行うことなくパーミュテーションの問題を解消することができる。以下では先ず、多次元確率密度関数を用いることでパーミュテーションの問題が解消することの理論的根拠及び具体的な計算式について説明し、次いで、本実施の形態における音声信号分離装置の具体的構成について説明する。

先ず、多次元確率密度関数を用いることでパーミュテーションの問題が解消することの理論的根拠について図１を用いて説明する。なお、図１では簡単のため、チャンネル数を２（ｎ＝２）とし、周波数binの総数を３（Ｍ＝３）としているが、任意のｎ，Ｍについて同様の説明が適用可能である。

図１において、周波数bin毎の分離が成功し、しかもパーミュテーションが発生していない場合をケース１とし、周波数bin毎の分離は成功したがω＝２でパーミュテーションが発生している場合をケース２とする。

従来のように、周波数bin毎に計算したＫＬ情報量Ｉ（Ｙ（ω））を最小化する場合、ケース２のω＝２でパーミュテーションが発生しているにも関わらずケース１とケース２とでＩ（Ｙ（２））が同一の値となってしまう。すなわち、ＫＬ情報量Ｉ（Ｙ（ω））と分離行列Ｗ（ω）との関係を模式的に示すと図２（Ａ）のようになり（但し、実際にはＷ（ω）を１本の軸で表すことはできない）、ケース１，２の何れもＫＬ情報量の最小値をとるため、両者を区別することができない。これが、従来手法によりパーミュテーションが発生する本質的な原因である。

これに対して、本実施の形態における音声信号分離装置では、多次元確率密度関数を用いてチャンネル毎にエントロピーを計算し、全チャンネルで１つのＫＬ情報量を計算する（式の詳細な説明は後述する）。このように、本実施の形態では全チャンネルで１つのＫＬ情報量が計算されるため、ケース１とケース２とでＫＬ情報量は異なる値をとる。さらに、適切な多次元確率密度関数を用意すれば、ケース１のＫＬ情報量をケース２のＫＬ情報量よりも小さくすることができる。すなわち、ＫＬ情報量Ｉ（Ｙ）と分離行列Ｗ（ω）との関係を模式的に示すと図２（Ｂ）のようになり、ケース１とケース２とを区別することができる。したがって、従来のような後処理による入れ替えを行わなくても、ＫＬ情報量を最小化するだけで信号を分離すると共にパーミュテーションの発生も防ぐことができる。

なお、本実施の形態においても、全ての周波数binでＹ_１＝Ｓ_２、Ｙ_２＝Ｓ_１のように分離される場合（ケース３とする）とケース１とではＫＬ情報量が同一の値となるため、区別することができない。しかしながら、ケース３にはパーミュテーションが発生していないため、分離結果がケース３であっても問題はない。

ここで、時間周波数領域の独立成分分析に多次元確率密度関数を導入するには、（ａ）分離行列をどのような式で更新するか、（ｂ）複素数への対応、（ｃ）どのような多次元確率密度関数を用いるか、の３点を解決する必要がある。以下では、上記３点について順に説明すると共に、（ｄ）変形例、についても併せて説明する。

（ａ）分離行列Ｗの更新式
上記式（５）〜（９）は１次元確率密度関数を用いた式であるため、そのままでは多次元確率密度関数に適用することはできない。そこで、以下の手順で、多次元確率密度関数を用いた分離行列Ｗの更新式を導出する。

観測信号Ｘと分離信号Ｙとの関係を表した上記式（４）を、全てのω（１≦ω≦Ｍ）について書き下し、それらを１つの式で表現すると、下記式（１２）又は下記式（１５）のようになる（但し、以下では式（１２）を用いる。）。式（１２）のベクトル及び行列をそれぞれ１つの変数で表記すると下記式（１３）のようになる。また、下記式（１２）の同一のチャンネルに由来するベクトルや行列をそれぞれ１つの変数で表記すると下記式（１４）のようになる。この式（１４）において、Ｙ_ｋ（ｔ）はスペクトログラムから１フレーム分を切り出して作った列ベクトルを表し、Ｗ_ｉｊはｗ_ｉｊ（１），・・・，ｗ_ｉｊ（Ｍ）を対角要素とする対角行列を表す。

本実施の形態では、上記式（１２）〜（１４）のＹ_ｋ（ｔ），Ｙ（ｔ）を用いて、ＫＬ情報量Ｉ（Ｙ）を下記式（１６）のように定義する。この式（１６）において、Ｈ（Ｙ_ｋ）は各チャンネルについてのスペクトログラム１枚分のエントロピーを表し、Ｈ（Ｙ）は全チャンネルについてのスペクトログラム１枚分の同時エントロピーを表す。ｎ＝２のときのＨ（Ｙ_ｋ）とＨ（Ｙ）との関係を図３に示す。式（１６）のうち、Ｈ（Ｙ_ｋ）はエントロピーの定義により下記式（１７）の第１項のように書き換えられ、Ｈ（Ｙ）は上記式（１３）により下記式（１７）の第２項及び第３項のように展開される。この式（１７）において、Ｐ_Ｙｋ（・）はＹ_ｋ（１，ｔ），・・・，Ｙ_ｋ（Ｍ，ｔ）のＭ次元確率密度関数を表し、Ｈ（Ｘ）は観測信号Ｘの同時エントロピーを表す。

観測信号Ｘを分離するには、ＫＬ情報量Ｉ（Ｙ）を最小にするような分離行列Ｗを求めればよく、そのような分離行列Ｗは、下記式（１８）、（１９）に従ってＷを少しずつ更新することで求めることができる。

ここで、Ｗの更新は上記式（１２）で非０の要素に対してのみ行えばよい。そこで、ΔＷとＷとから周波数bin＝ωの成分のみを取り出した行列ΔＷ（ω），Ｗ（ω）を下記式（２０）、（２１）のように定義し、下記式（２２）に従ってΔＷ（ω）を計算する。式（２２）を全てのωについて計算すれば、ΔＷの非０の要素は全て計算できたことになる。この式（２２）において、φ_ω（・）は多次元確率密度関数に対応したスコア関数を表し、下記式（２３）を経て下記式（２４）のように計算される。すなわち、多次元確率密度関数の対数をω番目の引数で偏微分することで得られる。

上記式（８）と上記式（２２）との違いは、スコア関数の引数にある。上記式（８）のφ（・）の引数は周波数bin＝ωの成分のみであるため、他の周波数binとの相関を反映させることができない。これに対して上記式（２２）のφ_ω（・）の引数は全ての周波数binの成分であるため、他の周波数binとの相関を反映させることが可能となる。

なお、詳しくは後述するが、Ｙは複素数の信号であるため、実際には上記式（２２）の代わりに複素数に対応させた式を用いる。

ここで、分離行列Ｗの更新を繰り返すと、用いる多次元確率密度関数の種類によっては、要素の値がオーバーフローしてしまうことがある。

そこで、上記式（２２）におけるΔＷの式を以下のように変更し、分離行列Ｗの要素のオーバーフローを防止するようにしても構わない。

上記式（２０）、（２１）における行列ΔＷ（ω），Ｗ（ω）のｋ行目を取り出した行ベクトルΔＷ_ｋ（ω），Ｗ_ｋ（ω）を下記式（２５）、（２６）のように定義する。

Ｗ_ｋ（ω）は、観測信号Ｘのω番目の周波数binからチャンネルｋ、周波数bin＝ωの分離信号Ｙを生成するためのベクトルであるが、信号が分離されたか否かはＷ_ｋ（ω）の要素間の比率（観測信号間の混合比）で決まり、Ｗ_ｋ（ω）の大きさとは関係がない。例えば、観測信号を−１：２で混合するのも−２：４で混合するのも、信号の分離という点では同じことである。図４に示すように、ΔＷ_ｋ（ω）をＷ_ｋ（ω）に直交する成分ΔＷ_ｋ（ω）^［Ｃ］と、Ｗ_ｋ（ω）と平行な成分ΔＷ_ｋ（ω）^［Ｐ］とに分解した場合、ΔＷ_ｋ（ω）^［Ｃ］は信号の分離に寄与するが、ΔＷ_ｋ（ω）^［Ｐ］はＷ_ｋ（ω）を大きくするだけであり、信号の分離には寄与しない。また、Ｗ_ｋ（ω）が大きくなり過ぎると、上述のようにオーバーフローを起こす可能性が高くなる。

したがって、ΔＷ_ｋ（ω）を用いてＷ_ｋ（ω）を更新する代わりに、ΔＷ_ｋ（ω）^［Ｃ］のみを用いてＷ_ｋ（ω）を更新することにより、オーバーフローを防止しつつ、信号を分離することができるようになる。

具体的には、下記式（２７）によってΔＷ_ｋ（ω）^［Ｃ］を計算し、下記式（２８）のようにΔＷ_ｋ（ω）^［Ｃ］からなる行列ΔＷ（ω）^［Ｃ］を用いてＷ（ω）を更新する。

勿論、下記式（２９）のように、Ｗに直交する成分ΔＷ^［Ｃ］を用いてＷを更新するようにしても構わない。また、Ｗと平行な成分ΔＷ^［Ｐ］を全く無視するのではなく、下記式（３０）のように、ΔＷ^［Ｃ］，ΔＷ^［Ｐ］に対してそれぞれ異なる係数η_１，η_２（η_１＞η_２＞０）を乗じて、Ｗを更新するようにしても構わない。

（ｂ）複素数への対応
時間周波数領域の独立成分分析では複素数の信号を扱うため、Ｗの更新式を複素数に対応させる必要がある。ここで、従来の１次元確率密度関数を用いた方法では、前述した式（８）を複素数に対応させた下記式（３１）が提案されている（特開２００３−８４７９３号公報を参照）。この式（３１）において、上付き文字の“Ｈ”は共役転置（ベクトルを転置すると共に要素を共役複素数に置き換える）を表す。

しかしながら、多次元確率密度関数を用いた方法には上記式（３１）を適用することができない。そこで、本実施の形態では、下記式（３２）を新たに考案し、この式（３２）に基づいて分離行列Ｗを更新する。なお、下記式（３３）のφ_kω（・）はＭ個の引数をとる関数として表されているが、これは上記式（２４）のφ_kω（Ｙ_ｋ（ｔ））（Ｍ次元のベクトルを引数とする関数）と等価である。式（３３）のように、関数自体には各引数の絶対値を代入し、関数の返値にω番目の引数の位相成分Ｙ_ｋ（ω，ｔ）／|Ｙ_ｋ（ω，ｔ）|を乗じることで、スコア関数を複素数に対応させることができる。

上記式（３２）においても、上記式（２７）と同様にＷ（ω）に直交する成分ΔＷ（ω）^［Ｃ］を計算ようにしてもよいことは勿論である。

なお、後述の通り、多次元確率密度関数やスコア関数の種類によっては、始めから複素数の入力（引数）に対応しているものもある。そのような関数に対しては上記式（３３）の変形は不要であり、その場合、φハット（＾）はφと同一と見なす。

（ｃ）用いる多次元確率密度関数
多次元確率密度関数として有名なものに下記式（３４）で表される多次元（多変量）正規分布がある。この式（３４）において、ｘはｘ_１，・・・，ｘ_ｄの列ベクトルを表し、μはｘの平均値ベクトルを表し、Σはｘの分散共分散行列を表す。

しかしながら、独立成分分析では正規分布を確率密度関数として用いると信号が分離できないことが知られているため、正規分布以外の多次元確率密度関数を用いる必要がある。そこで、本実施の形態では、以下に説明するように、（ｉ）球状分布、（ii）Ｌ_Ｎノルム、（iii）楕円分布、（iv）Copulaモデル、に基づいて多次元確率密度関数を構築する。

（ｉ）球状分布
球状分布とは、任意の非負関数ｆ（ｘ）（ｘはスカラー）にベクトルのＬ２ノルムを代入して多次元化した確率密度関数のことである。Ｌ２ノルムとは、要素の絶対値の２乗を総和し、その結果の２乗根をとったものである。本実施の形態では、ｆ（ｘ）として主に１次元確率密度関数（指数分布や１／cosh（ｘ）など）を用いる。したがって、球状分布に基づく確率密度関数は下記式（３５）のように表される。この式（３５）において、ｈは全引数について−∞〜＋∞の区間で定積分した結果を１に調整するための定数であるが、スコア関数を求める際に約分されて消えるため、具体的な値を求める必要はない。なお、以下ではｆ（ｘ）の導関数をｆ’（ｘ）と表記する。

上記式（３５）の確率密度関数に対応したスコア関数は、以下の手順で求めることができる。確率密度関数の対数をベクトルｘで偏微分すると、下記式（３６）のような関数ｇ（ｘ）が得られる（但し、ｘはベクトル）。ｇ（ｘ）にＹ_ｋ（ｔ）を代入したｇ（Ｙ_ｋ(ｔ)）は、全ての周波数binのスコア関数を含んでいる。すなわち、ｇ（Ｙ_ｋ(ｔ)）＝［φ_ｋ１（Ｙ_ｋ(ｔ)），・・・，φ_ｋＭ（Ｙ_ｋ(ｔ)）］^Ｔの関係がある。したがって、下記式（３７）のようにｇ（Ｙ_ｋ(ｔ)）からω行目の要素を抽出することで、スコア関数φ_ｋω（Ｙ_ｋ(ｔ)）が得られる。なお、球状分布は要素の絶対値を用いている関係上、始めから複素数の入力にも対応しているため、上記式（３３）の変形は不要である。

ｆ（ｘ）に具体的な数式を代入した例を示す。
ｆ（ｘ）が下記式（３８）のような１次元の指数分布で表されるとする。この式（３８）において、Ｋはスカラー変数ｘの分布の広がりに対応した定数であるが、Ｋ＝１としても構わない。或いは、Ｙ_ｋ（ｔ）のＬ２ノルム||Ｙ_ｋ（ｔ）||_２の分布の広がりに応じてＫの値を変更しても構わない。この式（３８）を球状分布で多次元化すると、下記式（３９）のような確率密度関数が得られ、対応するｇ（Ｙ_ｋ(ｔ)）は下記式（４０）で表される。

また、ｆ（ｘ）が下記式（４１）で表されるとする。この式（４１）において、ｄは正の値である。この式（４１）を球状分布で多次元化すると、下記式（４２）のような確率密度関数が得られ、対応するｇ（Ｙ_ｋ(ｔ)）は下記式（４３）で表される。

（ii）Ｌ_Ｎノルム
上述した任意の非負関数ｆ（ｘ）（ｘはスカラー）にベクトルのＬ_Ｎノルムを代入して多次元化することで、Ｌ_Ｎノルムに基づく多次元確率密度関数を構築することができる。Ｌ_Ｎノルムとは、要素の絶対値のＮ乗を総和し、その結果のＮ乗根をとったものである。Ｙ_ｋ（ｔ）のＬ_Ｎノルム||Ｙ_ｋ（ｔ）||_Ｎを非負関数ｆ（ｘ）に代入して多次元化すると、下記式（４４）のような多次元確率密度関数が得られる。この式（４４）において、ｈは全引数について−∞〜＋∞の区間で定積分した結果を１に調整するための定数であるが、スコア関数を求める際に約分されて消えるため、具体的な値を求める必要はない。上述した球状分布は、このＬ_Ｎノルムに基づく多次元確率密度関数においてＮ＝２とした場合に相当する。

また、上記式（４４）から複素数に対応したスコア関数を導出すると、下記式（４５）が得られる。

上記式（４５）において、ｆ（ｘ）が下記式（４６）のような１次元の指数分布で表されるとすると、下記式（４７）のようなスコア関数が導出される。また、ｆ（ｘ）が下記式（４８）で表されるとすると、下記式（４９）のようなスコア関数が導出される。この式（４６）、（４８）において、Ｋは正の実数であり、ｄ，ｍは自然数である。

上記式（４７）、（４９）においてＮ＝２，ｍ＝１とすると、上述した球状分布の場合と同じスコア関数が得られ、後述のように、パーミュテーションが発生することなく観測信号を分離することができる。しかしながら、上記式（４７）、（４９）においてＮ＝１，ｍ＝１とすると、分離結果にパーミュテーションが発生してしまう。これは、上記式（４７）、（４９）の ||Ｙ_ｋ（ｔ）||_Ｎ ^{（ｍ−Ｎ）}という項がＮ＝ｍの場合は消えてしまい、周波数bin間の相関がΔＷ（ω）にあまり反映されなくなるためと考えられる。また、Ｎ≠ｍであっても、||Ｙ_ｋ（ｔ）||_Ｎ＝０の場合、すなわちｔ番目のフレームに信号が存在しない場合には、演算中に０除算が発生してしまう。

そこで、本実施の形態では、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たすように、スコア関数φ_ｋω（Ｙ_ｋ(ｔ)）の式を変更する。

ここで、スコア関数φ_ｋω（Ｙ_ｋ(ｔ)）の返値が無次元量とは、Ｙ_ｋ（ω，ｔ）の単位を［ｘ］としたとき、スコア関数の分子と分母とで［ｘ］が相殺され、返値には［ｘ］の次元（ｎを非ゼロの値としたときに［ｘ^ｎ］と記述される単位）が含まれないことを表す。

また、スコア関数φ_ｋω（Ｙ_ｋ(ｔ)）の返値の位相がω番目の位相と逆位相であるとは、arg｛φ_ｋω（Ｙ_ｋ(ｔ)）｝＝−arg｛Ｙ_ｋ(ω，ｔ)｝が任意のＹ_ｋ（ω，ｔ）について成立することを表す。但し、arg｛ｚ｝は複素数ｚの位相成分を表す。例えば、大きさｒと位相角 θとを用いてｚ＝ｒ・exp（ｉθ）と表した場合、arg｛ｚ｝＝θである。

なお、本実施の形態では、上記式（２２）、（３２）のように、ΔＷ（ω）＝｛Ｉ_ｎ＋Ｅ_ｔ[...]｝Ｗ（ω）としているため、スコア関数の条件は、返値の位相がω番目の位相と「逆位相」となるが、ΔＷ（ω）＝｛Ｉ_ｎ−Ｅ_ｔ[...]｝Ｗ（ω）とした場合には、スコア関数の符号が反転するため、スコア関数の条件は、返値の位相がω番目の位相と「同位相」となる。何れの場合であっても、スコア関数は、返値の位相がω番目の位相にのみ依存するものであればよい。

上述したスコア関数の返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件は、上記式（３３）を一般化したものであるため、スコア関数がこの２つの条件を満たしている場合には、上記式（３３）の複素数対策は不要である。

以下、具体例を挙げて説明する。
上述の通り、上記式（４７）、（４９）は、Ｌ_Ｎノルムに基づく多次元確率密度関数から導出されたスコア関数である。これらのスコア関数は、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たしているため、Ｎ≠ｍではパーミュテーションが発生することなく観測信号を分離することができる。しかしながら、上述の通り、Ｎ＝ｍでは||Ｙ_ｋ（ｔ）||_Ｎ ^{（ｍ−Ｎ）}という項が消えてしまうため、分離結果にパーミュテーションが発生してしまう。また、Ｎ≠ｍであっても、||Ｙ_ｋ（ｔ）||_Ｎ＝０の場合、すなわちｔ番目のフレームに信号が存在しない場合には、演算中に０除算が発生してしまう。

そこで、Ｎ＝ｍの場合であっても返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たし、さらに演算中に０除算が発生しないように、上記式（４７）、（４９）をそれぞれ下記式（５０）、（５１）のように変更する。この式（５０）、（５１）において、Ｌは正の定数であり、例えばＬ＝１とする。また、ａは０除算を防止するための定数であり、値は非負である。

上記式（５０）、（５１）では、||Ｙ_ｋ（ｔ）||_Ｎの項がＮ＝ｍの場合にも残る。また、||Ｙ_ｋ（ｔ）||_Ｎ＝０の場合にも０除算は発生しない。

上記式（５０）、（５１）において、Ｙ_ｋ（ω，ｔ）の単位を［ｘ］とすると、［ｘ］を持つ量は分子と分母とで同数（何れもＬ＋１回）出現するため、相殺されてスコア関数全体では無次元量となる（tanhは無次元量と見なす）。さらに、これらの式の返値の位相は、−Ｙ_ｋ（ω，ｔ）の位相と等しいため、返値の位相はＹ_ｋ（ω，ｔ）の位相と逆位相となる。したがって、上記式（５０）、（５１）で表されるスコア関数は、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たす。

なお、Ｙ_ｋ（ｔ）のＬ_Ｎノルム||Ｙ_ｋ（ｔ）||_Ｎを計算する際には、複素数の絶対値を求める必要があるが、下記式（５２）、（５３）に示すように、複素数の絶対値を実部又は虚部の絶対値で近似してもよく、下記式（５４）に示すように、両者の和で近似してもよい。

ここで、複素数を実部と虚部とに分解して保持しているシステムにおいて、ｚ＝ｘ＋ｉｙ（ｘ，ｙは実数、ｉは虚数単位）で表される複素数ｚの絶対値は下記式（５５）のように計算される。これに対して実部や虚部の絶対値は、下記式（５６）、（５７）のように計算されるため、計算量が削減される。特に、Ｌ１ノルムの場合には、２乗やルートを用いずに、実数の絶対値と和のみで計算できるため、計算を非常に簡略化することができる。

また、Ｌ_Ｎノルムの値は、Ｙ_ｋ（ｔ）のうちで絶対値の大きな成分によってほぼ決まるため、Ｌ_Ｎノルムの計算の際、Ｙ_ｋ（ｔ）の全ての成分を用いるのではなく、絶対値の大きな成分の上位ｘ％のみを用いるようにしてもよい。この上位ｘ％は、観測信号のスペクトログラムから事前に求めることができる。

（iii）楕円分布
楕円分布とは、下記式（５８）に示すように、列ベクトルｘのマハラノビス距離sqrt（ｘ^ＴΣ^−１ｘ）を任意の非負関数ｆ（ｘ）（ｘはスカラー）に代入することで生成される多次元確率密度関数のことである。Ｙ_ｋ（ｔ）を非負関数ｆ（ｘ）に代入して多次元化すると、下記式（５９）のような多次元確率密度関数が得られる。この式（５９）において、Σ_ｋはＹ_ｋ（ｔ）の分散共分散行列である。

また、上記式（５９）からスコア関数を導出すると、下記式（６０）が得られる。この式（６０）において、（・）_ωは括弧内のベクトルや行列のω行目を抽出することを表す。なお、楕円分布の場合、Ｙ_ｋ（ｔ）の要素が複素数であってもマハラノビス距離は非負の実数しかとらないため、上記式（３３）の複素数対策は不要である。

上記式（６０）において、ｆ（ｘ）が下記式（６１）で表されるとすると、下記式（６２）のようなスコア関数が導出される。この式（６１）において、Ｋは正の実数であり、ｄ，ｍは自然数である。

しかしながら、上記式（６２）を用いて信号を分離しようとすると、分離行列Ｗの更新を繰り返すうちに要素の値がオーバーフローしてしまう。これは、Ｗ←αＷ（α＞１）という更新（新しいＷは前回のＷのスカラー倍）が一度でも発生すると、以降のＷは相似的な拡大しか起こらなくなり、やがて計算機で扱える値の範囲を超えてしまうからである。

そこで、本実施の形態では、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たすように、スコア関数φ_ｋω（Ｙ_ｋ(ｔ)）の式を変更する。

ここで、上記式（６２）で表されるスコア関数は、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たしていない。すなわち、Ｙ_ｋ（ω，ｔ）の単位を［ｘ］とすると分散共分散行列Σ_ｋの単位は［ｘ^２］であるため、スコア関数全体では［１／ｘ］の次元を持つ。また、分子に現れる（Σ_ｋ ^−１Ｙ_ｋ(ｔ)）_ω の演算では、Ｙ_ｋ(ｔ)のうちでＹ_ｋ（ω，ｔ）以外の成分も加算されるため、返値の位相は−Ｙ_ｋ（ω，ｔ）とは異なったものとなる。

そこで、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たし、さらに演算中に０除算が発生しないように、上記式（６２）を下記式（６３）のように変更する。この式（６３）において、Ｌは正の定数であり、例えばＬ＝１とする。また、ａは０除算を防止するための定数であり、値は非負である。

特に、ｆ（ｘ）が上記式（６１）で表され、Ｌ＝１である場合に導出されるスコア関数を下記式（６４）に示す。

なお、Ｙ_ｋ(ｔ)の分布によっては、分散共分散行列Σ_ｋの逆行列が存在しない場合がある。そこで、Σ_ｋの代わりにdiag（Σ_ｋ）（Σ_ｋの対角要素からなる行列）を用いたり、逆行列 Σ_ｋ ^−１の代わりに一般逆行列（例えば、Moore-Penrose 型一般逆行列）を用いても構わない。

（iv）Copulaモデル
Sklarの定理によると、任意の多次元累積分布関数Ｆ（ｘ_１，・・・，ｘ_ｄ）は、ある性質を持ったｄ引数関数Ｃ（ｘ_１，・・・，ｘ_ｄ）と、各引数の周辺分布関数Ｆ_ｋ（ｘ_ｋ）とを用いて、下記式（６５）の右辺のように変形することができる。このＣ（ｘ_１，・・・，ｘ_ｄ）をCopulaという。言い換えると、Copula Ｃ（ｘ_１，・・・，ｘ_ｄ）と任意の周辺分布関数Ｆ_ｋ（ｘ_ｋ）とを組み合わせることで、様々な多次元累積分布関数を構築することができる。なお、Copulaについては、例えば「“COPULAS”（http://gompertz.math.ualberta.ca/copula.pdf）」、「“The Shape of Neural Dependence”（http://wavelet.psych.wisc.edu/Jenison_Reale_Copula.pdf）」、「“Estimation and Model Selection of Semiparametric Copula-Based Multivariate Dynamic Models Under Copula Misspecification”（http://www.nd.edu/~meg/MEG2004/Chen-Xiaohong.pdf）」等のドキュメントに記載されている。

以下、Copulaを用いた多次元確率密度関数の構築法と、分離行列Ｗの更新式とについて説明する。

累積分布関数（Cumulative Distribution Function；ＣＤＦ）の上記式（６５）を全ての引数で偏微分すると、下記式（６６）のような確率密度関数が得られる。この式（６６）において、Ｐ_ｊ（ｘ_ｊ）は引数ｘ_ｊの確率密度関数であり、ｃ'（・）はCopulaを全引数で偏微分したものである。

この確率密度関数の対数をω番目の引数で偏微分すると、下記式（６７）のようなスコア関数が得られる。これが、Copulaを用いた多次元スコア関数の一般式である。この式（６７）において、Ｆ_Ｙｋ(ω)（・）はＹ_ｋ（ω，ｔ）の累積分布関数であり、Ｐ_Ｙｋ(ω)（・）はＹ_ｋ（ω，ｔ）の確率密度関数である。この式（６７）のｃ'（・）、Ｆ_Ｙｋ(ω)（・）、Ｐ_Ｙｋ(ω)（・）に具体的な式を代入することで、様々な多次元スコア関数を構築することができる。

例えば、Copulaの一種に下記式（６８）で表されるClayton's copulaがある。この式（６８）において、αは引数同士の依存度を表すパラメータである。式（６８）を全引数で偏微分すると下記式（６９）が得られ、それを上記式（６７）に代入すると、スコア関数である下記式（７０）が得られる。実際には、さらに上記式（３３）を適用することで、複素数に対応したスコア関数を得る。

Ｆ_Ｙｋ(ω)（・）、Ｐ_Ｙｋ(ω)（・）に具体的な式を代入した例を示す。

各周波数binの分布を指数分布と仮定すると、確率密度関数は下記式（７１）のように表すことができる。この式（７１）において、Ｋは分布の広がりに対応した変数であるが、Ｋ＝１としても構わない。また、指数分布の累積分布関数は下記式（７２）のように表すことができる。なお、上記式（３３）の複素数対策により、式（７２）の引数は非負であるとして構わない。これらを上記式（７０）に代入することにより、スコア関数である下記式（７３）が得られる。

なお、Copulaを用いたスコア関数では、球状分布、Ｌ_Ｎノルム、楕円分布を用いたスコア関数と異なり、周波数bin毎に異なる分布を適用することも可能である。例えば、周波数bin内の信号の分布がスーパーガウシアンであるかサブガウシアンであるかに応じて確率密度関数及び累積分布関数を切り換えることもできる。これは、例えば前述したextended infomax法でスコア関数を−［Ｙ_ｋ（ω，ｔ）＋tanh｛Ｙ_ｋ（ω，ｔ）｝］と−［Ｙ_ｋ（ω，ｔ）−tanh｛Ｙ_ｋ（ω，ｔ）｝］との間で切り換えることに相当する。

具体的には、スーパーガウシアン用の確率密度関数として下記式（７４）に示す指数分布を、累積分布関数として下記式（７５）を用意する。また、サブガウシアン用の確率密度関数として下記式（７６）を、累積分布関数としてWilliams近似と呼ばれる下記式（７７）を用意する。そして、ある周波数binの分布がスーパーガウシアンである場合には式（７４）と式（７６）とを用い、サブガウシアンである場合には式（７５）と式（７７）とを用いる。

（ｄ）変形例
上述した（ｃ）（ii）、（iii）では、Ｌ_Ｎノルム又は楕円分布に基づいてスコア関数を導出した後、返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たすように、スコア関数の式を変更したが、この２つの条件を満たすスコア関数を直接構築しても構わない。

そのようにして構築したスコア関数を下記式（７８）に示す。この式（７８）において、ｇ（ｘ）は以下のｉ）〜iv）の条件を満たす関数である。
ｉ） ｘ≧０においてｇ（ｘ）≧０
ii） ｘ≧０において、ｇ（ｘ）は定数、単調増加関数、又は単調減少関数
iii）ｇ（ｘ）が単調増加又は単調減少である場合、ｘ→∞においてｇ（ｘ）は正の値に収束する
iv） ｇ（ｘ）はｘに対して無次元量

分離に成功するｇ（ｘ）の例を下記式（７９）〜（８３）に示す。この式（７９）〜（８３）において、定数項は上述のｉ）〜iii）の条件を満たすように定める。

さらに一般化したスコア関数を下記式（８４）に示す。このスコア関数は、ベクトルＹ_ｋ(ｔ)を引数とする関数ｆ（Ｙ_ｋ(ｔ)）と、スカラーＹ_ｋ(ω，ｔ)を引数とする関数ｇ（Ｙ_ｋ(ω，ｔ)）と、返値の位相を決定するための項−Ｙ_ｋ(ω，ｔ)との積で表される関数である。但し、ｆ（Ｙ_ｋ(ｔ)）及びｇ（Ｙ_ｋ(ω，ｔ)）は、両者の積が任意のＹ_ｋ(ｔ)、Ｙ_ｋ(ω，ｔ)について以下のｖ）、vi）の条件を満たすように、それぞれ定める。
ｖ） ｆ（Ｙ_ｋ(ｔ)）及びｇ（Ｙ_ｋ(ω，ｔ)）は非負の実数
vi） ｆ（Ｙ_ｋ(ｔ)）及びｇ（Ｙ_ｋ(ω，ｔ)）の次元は［１／ｘ］
（Ｙ_ｋ(ω，ｔ)の単位を［ｘ］とする）

上述のｖ）の条件により、スコア関数の位相は−Ｙ_ｋ（ω，ｔ）と同一となり、スコア関数の返値の位相がω番目の位相と逆位相であるという条件が満たされる。また、上述のvi）の条件により、次元がＹ_ｋ（ω，ｔ）と相殺され、スコア関数が無次元量という条件が満たされる。

以上、多次元確率密度関数及びスコア関数の具体的な計算式について説明したが、以下では本実施の形態における音声信号分離装置の具体的な構成について説明する。

本実施の形態における音声信号分離装置の概略構成を図５に示す。この音声信号分離装置１において、ｎ個のマイクロホン１０_１〜１０_ｎは、ｎ個の音源が発する独立な音を観測し、Ａ／Ｄ（Analogue/Digital）変換部１１は、この信号をＡ／Ｄ変換して観測信号を得る。短時間フーリエ変換部１２は、観測信号を短時間フーリエ変換して観測信号のスペクトログラムを生成する。信号分離部１３は、信号モデル保持部１４に保持された信号モデルを利用して、観測信号のスペクトログラムを独立な信号に基づくスペクトログラムに分離する。信号モデルとは、具体的には上述の多次元確率密度関数のことであり、分離処理において分離信号のエントロピーを計算するために用いられる。但し、実際には、多次元確率密度関数そのものでなく、確率密度関数の対数を各引数で偏微分したスコア関数が信号モデル保持部１４に保持されていればよい。

リスケーリング部１５は、分離信号のスペクトログラムの各周波数binに対してスケールを揃える処理を行う。また、分離処理前に観測信号に対して標準化処理（平均や分散の調整）を施していた場合には元に戻す処理を行う。逆フーリエ変換部１６は、逆フーリエ変換によって分離信号のスペクトログラムを時間領域の分離信号に変換する。Ｄ／Ａ変換部１７は、時間領域の分離信号をＤ／Ａ変換し、ｎ個のスピーカ１８_１〜１８_ｎは、それぞれ独立の音を再生する。

なお、この音声信号分離装置１では、ｎ個のスピーカ１８_１〜１８_ｎを介して音を再生するものとしたが、分離信号を出力し、音声認識等に用いるようにすることも可能である。この場合には、逆フーリエ変換処理を適宜省略しても構わない。

この音声信号分離装置の処理の概略を図６のフローチャートを用いて説明する。先ずステップＳ１において、マイクロホンを介して音声信号を観測し、ステップＳ２において、観測信号を短時間フーリエ変換してスペクトログラムを得る。次にステップＳ３において、観測信号のスペクトログラムに対して各チャンネルの周波数bin毎に標準化を行う。この標準化とは、各周波数binの平均を０に、標準偏差を１に揃える操作である。周波数bin毎に平均値を減算することで平均を０にし、さらに標準偏差で除算することで標準偏差を１にすることができる。なお、多次元確率密度関数として球状分布を用いる場合には、他の方法による標準化も可能である。すなわち、周波数bin毎に平均を０にした後、ベクトルのノルム||Ｙ_ｋ（ｔ）||の１≦ｔ≦Ｔにおける標準偏差を求め、Ｙ_ｋをその値で割ることで、標準化を行うことができる。標準化後の観測信号をＸ’とすると、何れの標準化もＸ’＝Ｐ（Ｘ−μ）と表すことができる。ここで、Ｐは標準偏差の逆数からなる対角行列を表し、μは周波数bin毎の平均値からなるベクトルを表す。

続いてステップＳ４において、標準化された観測信号に対して分離処理を行う。具体的には、分離行列Ｗと分離信号Ｙとを求める。なお、このステップＳ４における処理の詳細は後述する。ステップＳ４で得られた分離信号Ｙは、パーミュテーションは発生していないものの、周波数bin毎にスケールが異なっている。そこでステップＳ５では、リスケーリング処理を行い、周波数binの間のスケールを揃える。ここでは、標準化処理で変更した平均と標準偏差とを元に戻す処理も行う。なお、ステップＳ５における処理の詳細は後述する。続いてステップＳ６において、リスケーリング後の分離信号を逆フーリエ変換によって時間領域の分離信号に変換し、ステップＳ７においてスピーカから再生する。

上述したステップＳ４（図６）における分離処理の詳細を図７及び図８を用いて説明する。図７はバッチ処理、図８はオンライン処理を行う場合におけるフローチャートを示したものである。ここで、バッチ処理とは信号全体を纏めて処理する方法であり、オンライン処理とは１サンプル（時間周波数領域の独立成分分析では１フレーム）入力される毎に逐次的に処理する方法である。なお、図７、図８におけるＸ（ｔ）は標準化された観測信号であり、図６のＸ’（ｔ）に相当する。

始めに、図７のバッチ処理を行う場合における分離処理について説明する。先ずステップＳ１１において分離行列Ｗに初期値を代入しておく。初期値としては、例えば単位行列を代入するようにしてもよく、上記式（２１）の全てのＷ（ω）に共通の行列を代入するようにしてもよい。次にステップＳ１２においてＷが収束したか否かを判別し、収束している場合には処理を終了し、収束していない場合にはステップＳ１３に進む。

続いてステップＳ１３においてその時点での分離信号Ｙを計算し、ステップＳ１４において上記式（３２）に従ってΔＷを計算する。このΔＷは周波数bin毎に計算されるため、ωのループを回し、それぞれのωについて上記式（３２）を適用する。ΔＷを求めたら、ステップＳ１５においてＷを更新し、ステップＳ１２に戻る。

なお、図７ではステップＳ１３，Ｓ１５が周波数binループの外側にある場合について説明したが、これらの処理を周波数binループの内側に移し、前述した図２０のステップＳ１０３，Ｓ１０５のように計算しても構わない。また、図７ではＷが収束するまでＷの更新処理を行うものとして説明したが、十分に大きな所定回数だけ繰り返すようにしても構わない。

次に、図８のオンライン処理を行う場合における分離処理について説明する。バッチ処理との違いは、１サンプル与えられる毎にΔＷを計算すること、及びΔＷの更新式から平均操作Ｅｔ［・］が消えていることである。すなわち、先ずステップＳ２１において、分離行列Ｗに初期値を代入しておく。次にステップＳ２２においてＷが収束したか否かを判別し、収束している場合には処理を終了し、収束していない場合にはステップＳ２３に進む。

続いてステップＳ２３においてその時点での分離信号Ｙを計算し、ステップＳ２４においてΔＷを計算する。このΔＷは周波数bin毎に計算されるため、ωのループを回し、それぞれのωについてΔＷを計算する。上述したように、このΔＷの更新式からは平均操作Ｅｔ［・］が消えている。ΔＷを求めたら、ステップＳ２５においてＷを更新する。ステップＳ２２〜Ｓ２５の処理は、各フレームについてωのループを回しながら全てのフレームについて繰り返される。

なお、ステップＳ２４におけるηは、固定値（例えば０．１）であってもよく、フレーム番号ｔが大きくなるにつれて小さくなるように調整してもよい。後者の場合、初めの方のフレームではηを大きく（例えば１）してＷの収束を速め、終わりの方のフレームではηを小さくして分離信号の急な変動を防ぐようにすることが好ましい。

次に、上述したステップＳ５（図６）におけるリスケーリング処理の詳細を図９を用いて説明する。従来、このリスケーリング処理も周波数bin毎に行っていたが、本実施の形態では、上記式（１３）のＷ，Ｘ，Ｙ等を用いて全ての周波数binに対して同時にリスケーリング処理を行っている。

上述したステップＳ４（図６）の分離処理が終了した時点で分離行列Ｗが求まっている。そこでステップＳ３１では、このＷに観測信号Ｘ'（ｔ）を乗じることで分離信号Ｙ'（ｔ）を得る。ステップＳ３１におけるＰは分散標準化行列である。なお、Ｘ'（ｔ）にＰμを加えているのは、ステップＳ３（図６）で平均を０にした観測信号を元に戻すためである。この段階では、まだスケーリングの問題は解消していない。

次にステップＳ３２において、分離信号から音源毎の観測信号を推定することでスケーリング問題を解決する。以下、その原理を説明する。

前述した図１５のような状況において、仮に音源ｋのみが音（原信号ｋ）を出力しているとする。各マイクロホンで観測される信号（音源毎の観測信号）は、音源ｋの信号に対して各マイクロホンまでの伝達関数を畳み込むことで得られる。ここで、原信号の推定とは異なり、音源毎の観測信号にはスケーリングの不定性がない。これは、原信号の推定では元々小さい原信号が減衰せずにマイクロホンに到達した場合と大きな原信号がマイクロホンに到達するまでに減衰した場合とが区別できないのに対して、音源毎の観測信号では両者を区別する必要がないからである。

推定された原信号でもある分離信号Ｙ’から音源毎の観測信号を推定する手順は以下の通りである。先ず、上記式（１４）の左辺のようにＹ’をチャンネル毎のベクトルＹ_１（ｔ）〜Ｙ_ｎ（ｔ）を使って表現する。次に、Ｙ’の中のＹ_ｋ（ｔ）以外を０ベクトルに置き換えたベクトルを作り、Ｙ_Ｙｋ（ｔ）とする。Ｙ_Ｙｋ（ｔ）は図１５で音源ｋのみが鳴っている状況に相当する。音源毎の観測信号は、Ｘ_Ｙｋ（ｔ）＝（ＷＰ）^−１Ｙ_Ｙｋ（ｔ）を計算することで得られる。この計算は全チャンネルについて繰り返し行われる。なお、Ｘ_Ｙｋ（ｔ）は、上記式（１４）の右辺第２項と同様に、全てのマイクロホンについての観測信号を含んでいる。

後段の処理では、Ｘ_Ｙｋ（ｔ）をそのまま使用してもよく、特定のマイクロホン（例えば１番目のマイクロホン）の観測信号のみを抽出してもよい。また、マイクロホン毎に信号のパワーを計算し、パワーが最大の信号を抽出しても構わない。これは、音源に最も近いマイクロホンで観測された信号を採用することにほぼ相当する。

以上詳細に説明したように、本実施の形態における音声信号分離装置１によれば、従来のような１次元確率密度関数を用いて周波数bin毎のエントロピーを計算する代わりに、多次元確率密度関数を用いてスペクトログラム１枚分のエントロピーを計算することにより、分離後の後処理を行うことなくパーミュテーションの問題を解消することができる。

以下、具体的な分離結果を示す。

球状分布に基づく多次元確率密度関数である上記式（４２）においてＫ＝π／２、ｄ＝１、ｈ＝１として分離した結果を図１０に示す。観測信号は前述した「X_rsm2.wav」というファイルの最初の３２０００サンプルであり、サンプリング周波数は１６ｋＨｚである。また、短時間フーリエ変換では、長さ１０２４のハニング窓をシフト幅１２８で使用している。したがって、周波数binの個数Ｍは１０２４／２＋１＝５１３であり、フレームの総数Ｔは（３２０００−１０２４）／１２８＋１＝２４３である。従来のextended infomax法を用いて分離した結果である図２１ではパーミュテーションが発生しているため後処理が必要であるのに対して、図１０では後処理をしていないにも関わらずパーミュテーションが殆ど発生していない。

また、Ｌ_Ｎノルムに基づくスコア関数である上記式（４９）においてＮ＝Ｋ＝ｄ＝ｍ＝１として分離した結果を図１１（Ａ）に示し、上記式（５１）においてＮ＝Ｋ＝ｄ＝ｍ＝１として分離した結果を図１１（Ｂ）に示す。観測信号は前述した「X_rsm2.wav」というファイルの最初の４００００サンプルであり、サンプリング周波数は１６ｋＨｚである。また、短時間フーリエ変換では、長さ５１２のハニング窓をシフト幅１２８で使用している。返値が無次元量であり、且つ、返値の位相がω番目の位相と逆位相であるという条件を満たしていない上記式（４９）を用いた場合には、図１１（Ａ）中矢印で示すように分離結果にパーミュテーションが発生しているが、この２つの条件を満たした上記式（５１）を用いた場合には、図１１（Ｂ）に示すように、後処理をしていないにも関わらずパーミュテーションは殆ど発生していない。

また、Copulaモデルに基づく多次元確率密度関数である上記式（７３）においてＫ＝１、α＝１として分離した結果を図１２に示す。観測信号及びサンプリング周波数等の条件は図１０と同様である。この場合も、後処理をしていないにも関わらずパーミュテーションは殆ど発生していない。

次に、上述の多次元確率密度関数と観測信号の分離結果とを用いて、図１，２のような状態が実現されているか否かを検証した結果を示す。すわなち、パーミュテーションが発生している状態と発生していない状態とを比較したときに、後者の方がＫＬ情報量が小さくなっているか否かを検証した結果を示す。

手順は以下の通りである。すなわち、先ず図１０に示すスペクトログラムを用意し、この状態のＫＬ情報量を上記式（１７）に従って計算する。なお、この実験においては上記式（１７）の第２，３項は定数と見なすことができ、パーミュテーションの有無には影響されないため、この実験では０としても構わない。次に、周波数binを任意に１本選択し、チャンネル間でその周波数binのデータを交換する。すなわち、人工的にパーミュテーションを発生させる。交換したら、再び上記式（１７）に従ってＫＬ情報量を計算する。周波数binの重複なしでこの操作を周波数binの総数と同じ回数だけ繰り返すと、最終的にはチャンネル間で全ての信号が入れ替わる。その過程を５段階で示したのが図１３（Ａ）〜（Ｅ）である。なお、図１３（Ａ）〜（Ｅ）はそれぞれ周波数binを０％、２５％、５０％、７５％、１００％置換したものである。

この操作の後、縦軸をＫＬ情報量、横軸を操作の回数、すなわち交換した周波数binの本数（パーミュテーションの程度でもある）としてプロットすると、図１４のようなグラフが得られる。但し、周波数binを選択する順番には任意性があるため、図１４では、（ａ）信号成分の大きい順に選択、（ｂ）ω＝１から順に選択、（ｃ）（ｄ）ランダムに選択、という４通りにより実験している。なお、（ａ）の「信号成分の大きい順」とは、下記式（８５）により周波数bin毎（ω毎）に計算される値Ｄ（ω）の大きい順のことであり、図１３もこの尺度に従ったものである。

図１４のグラフでは、４本のプロットは何れも両端が最小値となっている。すなわち、本実施の形態のように多次元確率密度関数を用いて信号を分離することにより、パーミュテーションが発生していない場合（両端）のＫＬ情報量がパーミュテーションが発生しているいかなる場合のＫＬ情報量よりも小さな値をとることが実際のデータからも裏付けられた。

言い換えれば、パーミュテーションの程度と、ある多次元確率密度関数を用いて計算されるＫＬ情報量との関係をプロットしたときに、両端（すなわち、パーミュテーションが発生していない状態）がＫＬ情報量の最小値となるならば、その確率密度関数（或いはその確率密度関数に対応したスコア関数）を用いることで、パーミュテーションを発生させることなく観測信号を分離することができる。

なお、本発明は上述した実施の形態のみに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々の変更が可能であることは勿論である。

例えば、全チャンネルに亘って信号が殆ど存在しない（０に近い成分しか存在しない）周波数binは、分離が成功してもしなくても時間領域の分離信号には殆ど影響しないため、そのような周波数binを省いてスペクトログラムを縮退させることで、計算量を削減し、分離処理を高速化することができる。

スペクトログラムを縮退させる一例としては、観測信号のスペクトログラムを生成した後、周波数bin毎に信号の絶対値が所定の閾値を上回っているか否かの判定を行い、全フレーム且つ全チャンネルにおいて閾値を下回っている周波数binを信号が存在しないと判定してスペクトログラムから除去する方法が挙げられる。但し、後で復元するため、何番目の周波数binを除去したかを記録しておく。信号が存在しない周波数binがｍ本あるとすると、除去後のスペクトログラムはＭ−ｍ本の周波数binを持つ。

また、スペクトログラムを縮退させる他の例としては、周波数bin毎に例えば上記式（５９）に従って信号の強さを計算し、強さの上位Ｍ−ｍ本を採用する（下位ｍ本を除去する）方法が挙げられる。

スペクトログラムを縮退させると、この縮退後のスペクトログラムに対して、標準化、分離処理、リスケーリング処理を行う。さらに、先ほど除去した周波数binを挿入する。なお、除去した信号の代わりに全ての成分が０というベクトルを挿入してもよい。この信号を逆フーリエ変換することで、時間領域の分離信号を得ることができる。

また、上述した実施の形態では、マイクロホンの数と音源数とが一致するものとして説明したが、マイクロホンの数が音源数よりも多い場合にも適用可能である。この場合には、例えば主成分分析（Principal Component Analysis；ＰＣＡ）を用いることで、マイクロホンの数を音源数まで減らすことができる。

また、上述した実施の形態では、分離行列の修正値ΔＷ（ω）を求めるアルゴリズムとして自然勾配法を用いたが、非ホロノームアルゴリズムに基づいてΔＷ（ω）を求めるようにしても構わない。ΔＷ（ω）を計算する式は、Ｂを適切な正方行列として、ΔＷ（ω）＝Ｂ・Ｗ（ω）と書き表すことができる。Ｂの対角成分が常に０となる式を用いている場合、その式を用いた更新式を非ホロノームアルゴリズムを呼ぶ。なお、非ホロノーム自体については、「岩波書店『統計科学のフロンティア5 多変量解析の展開』」等に記載されている。

非ホロノームアルゴリズムに基づくΔＷ（ω）の更新式を下記式（８６）に示す。この非ホロノームアルゴリズムを用いることで、Ｗは直交方向にのみ変化するようになるため、Ｗの演算中のオーバーフローを防止することができる。

多次元確率密度関数を用いることでパーミュテーションの問題が解消することの理論的根拠を説明する図である。 パーミュテーションの発生の有無によるＫＬ情報量の違いを従来と本実施の形態とで比較する図である。 本実施の形態におけるエントロピーと同時エントロピーとを説明する図である。 分離行列Ｗ(ω)の修正値ΔＷ（ω）の行ベクトルΔＷ_ｋ（ω）を、分離行列の行ベクトルＷ_ｋ（ω）に直交する成分ΔＷ_ｋ（ω）^［Ｃ］と平行な成分ΔＷ_ｋ（ω）^［Ｐ］とに分解した様子を示す図である。 本実施の形態における音声信号分離装置の概略構成を示す図である。 同音声信号分離装置の処理の概略を説明するフローチャートである。 バッチ処理を行う場合における分離処理の詳細を説明するフローチャートである。 オンライン処理を行う場合における分離処理の詳細を説明するフローチャートである。 リスケーリング処理の詳細を説明するフローチャートである。 球状分布に基づく多次元確率密度関数を用いて信号を分離した結果を示す図である。 Ｌ_Ｎノルムに基づくスコア関数を用いて信号を分離した結果を示す図である。 Copulaモデルに基づく多次元確率密度関数を用いて信号を分離した結果を示す図である。 得られた分離信号に対して人工的にパーミュテーションを発生させた場合のスペクトログラムの変化を示す図である。 得られた分離信号に対して人工的にパーミュテーションを発生させた場合のＫＬ情報量の変化を示す図である。 Ｎ個の音源から出力された原信号をｎ個のマイクロホンで観測する状況を示す図である。 時間周波数領域における従来の独立成分分析の概略を示す図である。 観測信号及びそのスペクトログラムと分離信号及びそのスペクトログラムとを示す図である。 ある周波数binに着目した場合における観測信号と分離信号とを示す図である。 従来のエントロピーと同時エントロピーとを説明する図である。 従来の分離処理の詳細を説明するフローチャートである。 １次元確率密度関数を用いて信号を分離した結果を示す図である。 周波数カップリングを行い、１次元確率密度関数を用いて信号を分離した結果を示す図である。

符号の説明

１ 音声信号分離装置、１０_１〜１０_ｎ マイクロホン、１１ Ａ／Ｄ変換部、１２ 短時間フーリエ変換部、１３ 信号分離部、１４ 信号モデル保持部、１５ リスケーリング部、１６ 逆フーリエ変換部、１７ Ｄ／Ａ変換部、１８_１〜１８_ｎ スピーカ

Claims (6)

1. 音声信号を含む複数の信号が混合された時間領域の観測信号を独立成分分析を用いて信号毎に分離し、分離信号を生成する音声信号分離装置において、
   上記時間領域の観測信号を時間周波数領域の観測信号に変換する第１の変換手段と、
   上記時間周波数領域の観測信号から時間周波数領域の分離信号を生成する分離手段と、
   上記時間周波数領域の分離信号を時間領域の分離信号に変換する第２の変換手段とを有し、
   上記分離手段は、上記時間周波数領域の観測信号と初期値が代入された分離行列とから時間周波数領域の分離信号を生成し、この時間周波数領域の分離信号と多次元確率密度関数を用いたスコア関数と上記分離行列とを用いて該分離行列の修正値を計算し、上記修正値を用いて、上記分離行列が略々収束するまで該分離行列を修正し、略々収束した分離行列を用いて上記時間周波数領域の分離信号を生成する
   ことを特徴とする音声信号分離装置。
2. 上記時間周波数領域の分離信号は複素信号であり、
   上記スコア関数として、返値の位相成分を１つの引数から計算し、返値の絶対値を１以上の引数から計算するスコア関数を用いることを特徴とする請求項１記載の音声信号分離装置。
3. 上記スコア関数は、返値が無次元量であり、且つ、返値の位相が１つの引数にのみ依存することを特徴とする請求項１記載の音声信号分離装置。
4. 音声信号を含む複数の信号が混合された時間領域の観測信号を独立成分分析を用いて信号毎に分離し、分離信号を生成する音声信号分離方法において、
   上記時間領域の観測信号を時間周波数領域の観測信号に変換する工程と、
   上記時間周波数領域の観測信号と初期値が代入された分離行列とから時間周波数領域の分離信号を生成する工程と、
   この時間周波数領域の分離信号と多次元確率密度関数を用いたスコア関数と上記分離行列とを用いて該分離行列の修正値を計算する工程と、
   上記修正値を用いて、上記分離行列が略々収束するまで該分離行列を修正する工程と、
   略々収束した分離行列を用いて生成された時間周波数領域の分離信号を時間領域の分離信号に変換する工程と
   を有することを特徴とする音声信号分離方法。
5. 上記時間周波数領域の分離信号は複素信号であり、
   上記スコア関数として、返値の位相成分を１つの引数から計算し、返値の絶対値を１以上の引数から計算するスコア関数を用いることを特徴とする請求項４記載の音声信号分離方法。
6. 上記スコア関数は、返値が無次元量であり、且つ、返値の位相が１つの引数にのみ依存することを特徴とする請求項４記載の音声信号分離方法。

Images