JP6581054B2 - 音源分離装置、音源分離方法及び音源分離プログラム - Google Patents

Description

本発明は、音源分離装置、音源分離方法及び音源分離プログラムに関する。

従来、音声認識などの応用に適用するため、１つ以上の目的音源から出た音響信号と、背景雑音や干渉音が混ざった観測信号と、に基づいて、各目的音源を分離する装置が提案されている（例えば、非特許文献１，２参照）。

N.Ito, S.Araki, T.Yoshioka, and T.Nakatani, "Relaxed Disjointness Based Clustering for Joint Blind Source Separation and Dereverberation", in Proc. Int. Worksh. Acoust. Echo, Noise Contr., pp. 268-272, 2014. P.Smaragdis and J.C.Brown, "Non-Negative Matrix Factorization for Polyphonic Music Transcription", in IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, pp. 177-180, 2003

例えば、非特許文献１に記載の音源分離装置について説明する。図４は、従来の音源分離装置の構成の一例を示すブロック図である。図４には、非特許文献１に開示された音源分離装置の機能構成を示す。図４に示す従来の音源分離装置１Ｐは、時間周波数分析部１１Ｐ、パラメータ推定部１２Ｐ、音源分離部１３Ｐ及び初期値生成部１４Ｐを有する。なお、残響除去部については、記載を省略する。

時間周波数分析部１１Ｐは、Ｎ（Ｎ＞１）個の目的音源に対応する音響信号が混在する状況において、異なる位置で収録されたＭ（Ｍ＞１）個の観測信号を受け取り、各観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出する。その後、時間周波数分析部１１Ｐは、時間周波数点ごとに、全ての観測音響信号をまとめたＭ次元縦ベクトルである観測ベクトルを構成する。言い換えると、時間周波数分析部１１Ｐは、２つ以上のマイクによって録音された音響信号を入力とし、短時間信号分析を適用して観測ベクトルを得る。

パラメータ推定部１２Ｐは、時間周波数分析部１１Ｐから観測ベクトルを受け取る。パラメータ推定部１２Ｐは、初期値生成部１４Ｐからパラメータの初期値を受け取る。そして、パラメータ推定部１２Ｐは、周波数ごとに、観測ベクトルの確率分布を、Ｎ個の音響信号のそれぞれに対応するＮ個の要素分布からなる混合分布でモデル化する。ここで、この分布パラメータのうち、各音源のパワーに相当するパラメータは、各時間周波数点で個別なパラメータとしてモデル化されている。そして、パラメータ推定部１２Ｐでは、分布パラメータ推定部１２１Ｐが、モデル化された混合分布が、時間周波数分析部１１Ｐから受け取った観測ベクトルの分布に近くなるように、混合分布のパラメータを推定する。その後、マスク推定部１２２Ｐは、各要素分布の事後確率を各音響信号に対応するマスクの推定値とする。

音源分離部１３Ｐは、観測ベクトルと目的音源に関するマスクを受け取り、観測ベクトルの１つの要素と、ｎ番目の要素分布に対応するマスクとを、各時間周波数成分において
掛け合わせることによって、ｎ番目の目的音声に対応する分離音を計算する。

初期値生成部１４Ｐは、分布パラメータの初期値を生成し、パラメータ推定部１２Ｐに出力する。

ここで、短時間フーリエ変換などの短時間信号分析を用いて得られる観測ベクトルをｙ_ｆ，ｔと表す。このうち、ｔ及びｆは、それぞれ時間と周波数の番号であり、ｔは、１〜Ｔの整数であり、ｆは、０〜Ｆの整数である。非特許文献１では、目的音源及び雑音は、スパース性を有し、各時間周波数点において高々一つの目的音源だけが存在すると仮定されている。このため、各時間周波数点における観測ベクトルｙ_ｆ，ｔは、以下の（１）式でモデル化できるとされている。

この（１）式のうち、s^（ｎ） _ｆ，ｔは、ｎ番目の音源の時間周波数成分であり、ｄ_ｆ，ｔは、時間周波数点（ｆ，ｔ）において優勢な音源の番号を表す。音源分離部１３Ｐは、この仮定に基づき、各時間周波数点の観測ベクトルｙ_ｆ，ｔを、Ｎ個の目的音源に対応するＮ個のクラスタに分類し、各音源のクラスタに属する成分のみそれぞれ抽出することによって、各目的音源に対応する分離音を得ることができる。以上の流れで音源分離を行うために、パラメータ推定部１２Ｐは、各時間周波数点の観測ベクトルｙ_ｆ，ｔがどれだけ各目的音源のクラスタに属しているらしいかを、以下の方法で推定する。パラメータ推定部１２Ｐは、まず、（２）式を用いて、観測ベクトルｙ_ｆ，ｔの確率分布をＮ個の要素分布からなる混合分布でモデル化する。

ここで、φ^（ｎ） _ｆ，ｔは、ｎ番目の目的音源のパワーに相当するスカラー値のパラメータである。また、Ｒ^（ｎ） _ｆは、ｎ番目の目的音源の時不変な空間相関行列に相当するＭ×Ｍの行列パラメータである。また、α^（ｎ） _ｆは、Σ_ｎα^（ｎ） _ｆ＝１を満たす重みパラメータ（スカラー値）である。上記の確率分布に基づき、各パラメータを推定するための尤度関数は、以下の（３）式で定義できる。

ただし、（３）式において、ｙは、観測ベクトルの集合を表し、Θは、分布パラメータの集合を表し、Ｆ_Ｎｃ（ｙ_ｆ，ｔ；０，φ^（ｎ） _ｆ，ｔＲ^（ｎ） _ｆ）は、観測ベクトルｙ_ｆ，ｔに関する平均０、分散共分散行列φ^（ｎ） _ｆ，ｔＲ^（ｎ） _ｆの複素正規分布の確率密度関数である。パラメータ推定部１２Ｐは、この目的関数が局所最大となる分布パラメータ、すなわち、当該混合分布が受け取った観測ベクトルの分布に近くなるような分布パラメータを求める。その後、パラメータ推定部１２Ｐが、目的音源に対応する要素分布の事後確率を、目的音源に対応するマスクの推定値とするのが、非特許文献１におけるマスク推定方法の基本的考え方である。Expectation-Maximization（ＥＭ）アルゴリズムの原理に基づき、対数尤度関数を局所最大化するためのＱ関数は、以下の（４）式のように定義できる。

ここで、λ^（ｎ） _ｆ，ｔは、ｎ番目の目的音源の事後確率に相当する補助パラメータであり、現在保持している分布パラメータを用いて、Ｅステップにおいて、以下の（５）式のように計算できる。

分布パラメータについては、Ｍステップにおいて上記のＱ関数をそれぞれのパラメータで偏微分したものを０とおくことで、以下の（６）式〜（８）式に示す更新則が導出できる。

パラメータ推定部１２Ｐは、Ｍステップにおける、上記の更新則による分布パラメータ更新と、Ｅステップにおける、更新されたパラメータによるλ^（ｎ） _ｆ，ｔの計算を反復的に行うことで、目的関数を局所最大化するような分布パラメータを推定することができる。そして、非特許文献１に記載の音源分離装置１Ｐでは、推定したλ^（ｎ） _ｆ，ｔを観測信号に掛け合わせることで、ｎ番目の目的音源に対応する分離音を得ている。

続いて、非特許文献２に記載の音源分離装置について説明する。図５は、従来の音源分離装置の構成の他の例を示すブロック図である。図５には、非特許文献２に開示された音源分離装置の機能構成を示す。図５に示す従来の音源分離装置１Ｐ’は、時間周波数分析部１１Ｐ’、パラメータ推定部１２Ｐ’、音源分離部１３Ｐ’及び初期値生成部１４Ｐ’を有する。

時間周波数分析部１１Ｐ’は、Ｎ（Ｎ＞１）個の目的音源に対応する音響信号が混在する状況において、１つのマイクで収録された１個の観測信号を受け取り、観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出する。

パラメータ推定部１２Ｐ’は、時間周波数分析部１１Ｐ’から観測音響信号を受け取る。パラメータ推定部１２Ｐ’は、初期値生成部１４Ｐ’からパラメータの初期値を受け取る。そして、パラメータ推定部１２Ｐ’は、Ｋ個の基底スペクトルと、それぞれの基底スペクトルに対応するアクティベーションの積和（これをモデルスペクトログラムと呼ぶ。）と、によってモデル化する。続いて、パラメータ推定部１２Ｐ’では、パラメータ推定処理部１２１Ｐ’が、モデルスペクトログラムが、時間周波数分析部１１Ｐ’から受け取った観測音響信号に近くなるように、基底スペクトルとアクティベーションとを推定する。

音源分離部１３Ｐ’は、観測音響信号と基底スペクトルとアクティベーションと、を受け取り、各目的音源に対応する分離音を計算する。なお、初期値生成部１４Ｐ’は、初期値生成部１４Ｐと同様の機能を有する。

ここで、短時間フーリエ変換などの短時間信号分析を用いて得られる観測音響信号を、ｘ_ｆ，ｔと表す。ｔ及びｆは、それぞれ時間と周波数の番号であり、ｔは、１〜Ｔの整数であり、ｆは、０〜Ｆの整数である。非特許文献２に記載の音源分離装置１Ｐ’では、観測音響信号を、限られた数のスペクトル（基底スペクトル）と、それぞれの基底スペクトルに対応するアクティベーションの積和（これをモデルスペクトログラムと呼ぶ。）によってモデル化する。このため、各時間周波数点における観測音響信号ｘ_ｆ，ｔは、以下の（９）式でモデル化できるとされている。

この（９）式のうち、ｗ_ｆ，ｋは、ｋ番目の基底スペクトルの周波数ｆでの成分を表し、ｈ_ｋ，ｔは、ｋ番目の基底スペクトルに対応するアクティベーションの時刻ｔでの成分を表す。そして、パラメータ推定部１２Ｐ’は、何らかの距離規準における、観測音響信号とモデルスペクトログラムとの間の距離を最小化することで、観測音響信号を近似するｗ_ｆ，ｋ，ｈ_ｋ，ｔを得る。

なお、Ｄ（ａ｜ｂ）は、ａとｂとの間の何らかの距離規準を表し、具体的には、二乗距離規準、カルバックライブラー距離規準、板倉斉藤距離規準などが用いられる。そして、音源分離装置１Ｐ’では、（１０）式で表される最適化問題を解くことによって得られた基底スペクトルとアクティベーションとを用いて、音源分離部１３Ｐ’が、各目的音源に対応する分離音を得る。具体的には、音源分離部１３Ｐ’として、ウィーナーフィルタなどが用いられる。

非特許文献１に記載の方法では、ｎ番目の音源信号の時間周波数点（ｆ，ｔ）におけるパワーに相当するパラメータφ^（ｎ） _ｆ，ｔは、各時間周波数点で個別に推定されている。しかしながら、非特許文献１に記載の方法では、音源信号の持つスペクトル構造（周波数方向の構造）は、マスクを推定する上でまったく考慮されていない。一方、非特許文献２に記載の方法は、音源信号の持つスペクトル構造に着目した音源分離方法の代表例として挙げられるものの、非特許文献２に記載の方法を非特許文献１に記載の方法に組み込むことは容易ではない。

例えば、非特許文献１に記載の方法を用いて音源を分離した後に、非特許文献２に記載の方法を用いて目的音源成分を抽出する方法は、自然な発想として考えられる。しかしながら、この場合、非特許文献１に記載の方法におけるマスクを推定する部分では、（５）式によって各時間周波数点のマスクが計算されるため、音源信号のもつスペクトル構造は、全く考慮されることがない。同様に、非特許文献２に記載の方法の後に非特許文献１に記載の方法を用いた場合では、前段で非特許文献２に記載の方法を用いるときに、複数マイクによって得られた空間的情報は、全く考慮されずに、音源分離のためのパラメータ（ｗ_ｆ，ｋ，ｈ_ｋ，ｔ）が推定されることになる。

本発明は、上記に鑑みてなされたものであって、音源分離を正確に行うことができる音響分離装置、音響分離方法及び音響分離プログラムを提供することを目的とする。

上述した課題を解決し、目的を達成するために、本発明に係る音響分離装置は、Ｎ（Ｎ＞１）個の目的音源に対応する音響信号が混在する状況において、異なる位置で収録されたＭ（Ｍ＞１）個の観測信号を受け取り、各観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出した後に時間周波数点ごとに全ての観測音響信号をまとめたＭ次元縦ベクトルである観測ベクトルを構成する時間周波数分析部と、時間周波数分析部から観測ベクトルを受け取り、各周波数における観測ベクトルの確率分布を、Ｎ個の音源に対応するＮ個の要素分布からなるとともに分布パラメータのうち各音源のパワーに相当するパラメータが周波数方向に構造を持った基底ベクトルとその時変であるアクティベーションとの積和によってモデル化された混合分布でモデル化し、モデル化した混合分布が、時間周波数分析部から受け取った観測ベクトルの分布に近くなるように混合分布のパラメータを推定したのちに各要素分布の事後確率を各音響信号に対応するマスクの推定値とするパラメータ推定部と、観測ベクトルと、マスクと、基底ベクトル、アクティベーション及び空間相関行列である分布パラメータ、或いは、その一部を受け取り、分離音を計算する音源分離部と、を有する。

本発明によれば、音源分離を正確に行うことができる。

図１は、実施の形態に係る音源分離装置の構成の一例を示すブロック図である。 図２は、図１に示す音源分離装置が実行する音源分離処理の処理手順を示すフローチャートである。 図３は、プログラムが実行されることにより、音源分離装置が実現されるコンピュータの一例を示す図である。 図４は、従来の音源分離装置の構成の一例を示すブロック図である。 図５は、従来の音源分離装置の構成の他の例を示すブロック図である。

以下、図面を参照して、本発明の一実施形態を詳細に説明する。なお、この実施の形態により本発明が限定されるものではない。また、図面の記載において、同一部分には同一の符号を付して示している。

［実施の形態］
実施の形態に係る音源分離装置について、音源分離装置の概略構成、音源分離装置における処理の流れを説明する。まず、図１を参照して、本実施の形態に係る音源分離装置の構成について説明する。

図１は、本実施の形態に係る音源分離装置の構成の一例を示すブロック図である。図１に示すように、実施の形態に係る音源分離装置１は、時間周波数分析部１１、パラメータ推定部１２、音源分離部１３及び初期値生成部１４を有する。

時間周波数分析部１１は、Ｎ（Ｎ＞１）個の目的音源に対応する音響信号が混在する状況において、異なる位置で収録されたＭ（Ｍ＞１）個の観測信号を受け取る。すなわち、時間周波数分析部１１は、多チャンネル観測信号を入力信号として受け取る。そして、時間周波数分析部１１は、各観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出する。この後に、時間周波数分析部１１は、時間周波数点ごとに全ての観測音響信号をまとめたＭ次元縦ベクトルである観測ベクトルを構成する。言い換えると、時間周波数分析部１１は、２つ以上のマイクによって録音された音響信号を入力とし、短時間信号分布を適用して観測ベクトルを得る機能を有する。

パラメータ推定部１２は、初期値生成部１４からパラメータの初期値を受け取る。パラメータ推定部１２は、時間周波数分析部１１から観測ベクトルを受け取る。そして、パラメータ推定部１２は、各周波数における観測ベクトルの確率分布を、Ｎ個の音源に対応するＮ個の要素分布からなる混合分布でモデル化する。ここで、分布パラメータのうち各要素分布における各音源のパワーに相当するパラメータは、周波数方向に構造を持った基底ベクトルとその時変であるアクティベーションとの積和によってモデル化されたものである。

続いて、パラメータ推定部１２は、モデル化した混合分布が、時間周波数分析部１１から受け取った観測ベクトルの分布に近くなるように混合分布のパラメータを推定したのちに、各要素分布の事後確率を各音響信号に対応するマスクの推定値とする。すなわち、パラメータ推定部１２は、観測ベクトルとパラメータの初期値とを受け取り、マスクと分布パラメータ（基底スペクトル、アクティベーション、空間相関行列）を更新する。

パラメータ推定部１２は、分布パラメータ推定部１２１及びマスク推定部１２２を有する。分布パラメータ推定部１２１は、モデル化した混合分布が、時間周波数分析部１１から受け取った観測ベクトルの分布に近くなるように混合分布のパラメータを推定する。そして、マスク推定部１２２は、各要素分布の事後確率を求め、該求めた事後確率を各音響信号に対応するマスクの推定値とする。

音源分離部１３は、観測ベクトルと、マスクと、分布パラメータ、或いは、その一部をパラメータ推定部１２から受け取り、分離音を計算する。分布パラメータは、基底ベクトル、アクティベーション及び空間相関行列である。

初期値生成部１４は、分布パラメータの初期値を生成し、パラメータ推定部１２に出力する。

音源分離装置１が実行する演算処理の概略について説明する。まず、短時間フーリエ変換などの短時間信号分析を用いて得られる観測ベクトルをｙ_ｆ，ｔと表す。ｔ及びｆは、それぞれ時間と周波数の番号であり、ｔは、１〜Ｔの整数であり、ｆは、０〜Ｆの整数である。そして、目的音源及び雑音は、スパース性を有し、各時間周波数点において高々一つの目的音源だけが存在すると仮定されている。本実施の形態では、この仮定を用いて、各時間周波数点における観測ベクトルｙ_ｆ，ｔを、以下の（１１）式でモデル化する。

この（１１）式のうち、s^（ｎ） _ｆ，ｔは、ｎ番目の音源の時間周波数成分であり、ｄ_ｆ，ｔは、時間周波数点（ｆ，ｔ）において優勢な音源の番号を表す。ｒ^（ｎ） _ｆは、周波数ｆにおいて、ｎ番目の音源から各マイクまでの伝達特性が並んだベクトルである。この仮定に基づき、各時間周波数点の観測ベクトルｙ_ｆ，ｔを、それぞれＮ個の目的音源に対応するＮ個のクラスタに分類し、各音源のクラスタに属する成分のみそれぞれ抽出することで、音源分離部１３は、各目的音源に対応する分離音を得ることができる。

以上の流れで音源分離を行うために、パラメータ推定部１２は、各時間周波数点の観測ベクトルがどれだけ各目的音源のクラスタに属しているらしいかを、以下の方法で推定する。まず、パラメータ推定部１２は、観測ベクトルの確率分布をＮ個の要素分布からなる混合分布でモデル化する。このとき、各要素分布は、各音源信号のクラスタに対応しているので、それぞれの音源のパワーに対応するパラメータを持つ。時間周波数点（ｆ，ｔ）におけるｎ番目のパワーに相当するパラメータφ^（ｎ） _ｆ，ｔを、以下の（１２）式のように、基底スペクトルとアクティベーションとの積和で表す。

ここで、ｗ_ｆ，ｋは、ｋ番目の基底スペクトルの周波数ｆでの成分を表す。ｈ_ｋ，ｔは、ｋ番目の基底スペクトルに対応するアクティベーションの時刻ｔでの成分を表す。パラメータ推定部１２は、以上で述べた性質を持つ確率分布が、観測ベクトルの確率分布に近くなるように、分布パラメータを推定することによって、音源信号の周波数方向の構造と空間的情報の双方を考慮しながら分布パラメータとマスクを推定している。この推定したマスクと分布パラメータとは、音源分離部１３で、各音源を分離するために使用される。

［時間周波数分析部の処理］
そこで、時間周波数分析部１１の処理の内容について説明する。例えば、Ｎ個の音源から出た音響信号をＭ個のマイクロホンで収録しているとする。この時、マイクロホンｍで収録された観測信号をｘ^（ｍ）（τ）とすると、ｘ^（ｍ）（τ）は、以下の（１３）式に示すように、各音源信号ｎに由来する音響信号ｓ^（ｍ） _ｎ（τ）の和で構成される。

時間周波数分析部１１は、全てのマイクロホンｍで収録された上記観測信号を受け取り、各観測信号ｘ^（ｍ）（τ）ごとに短時間信号分析を適用して時間周波数ごとの信号特徴量ｘ^（ｍ）（ｔ，ｆ）を求める。短時間信号分析としては、短時間離散フーリエ変換や短時間離散コサイン変換など様々な方法を用いることができる。時間周波数分析部１１は、さらに、各時間周波数で得られた信号特徴量ｘ^（ｍ）（ｔ，ｆ）を全てのマイクロホンｍに関してまとめたベクトルとして、以下の（１４）式に示す観測ベクトルｙ_ｔ，ｆを構成する。時間周波数分析部１１は、パラメータ推定部１２及び音源分離部１３に、観測ベクトルｙ_ｔ，ｆを出力する。

本実施の形態では、非特許文献１に記載の方法と同様に、目的音源及び雑音はスパース性を有し、各時間周波数点において高々一つの目的音源だけが存在すると仮定する。このため、各時間周波数点における観測ベクトルｙ_ｆ，ｔは、以下の（１５）式でモデル化できる。

この（１５）式のうち、s^（ｎ） _ｆ，ｔは、ｎ番目の音源の時間周波数成分を表す。ｄ_ｆ，ｔは、時間周波数点（ｆ，ｔ）において優勢な音源の番号を表す。ｒ^（ｎ） _ｆは、周波数ｆにおいて、ｎ番目の音源から各マイクまでの伝達特性が並んだベクトルである。この仮定に基づき、音源分離部１３は、各時間周波数点の観測ベクトルｙ_ｆ，ｔを、Ｎ個の目的音源に対応するＮ個のクラスタに識別し、各音源のクラスタに属する成分のみそれぞれ抽出することで、各目的音源に対応する分離音を得ることができる。また、以上の流れで音源分離を行うために、パラメータ推定部１２は、各時間周波数点の観測ベクトルがどれだけ各目的音源のクラスタに属しているらしいかを、以下の方法で推定する。

［パラメータ推定部の処理］
そこで、パラメータ推定部１２の処理について説明する。まず、パラメータ推定部１２は、観測ベクトルの確率分布を、以下の（１６）式に示すように、Ｎ個の要素分布からなる混合分布でモデル化する。

ここで、ｗ_ｆ，ｋは、ｋ番目の基底スペクトルの周波数ｆでの成分を表す。ｈ_ｋ，ｔは、ｋ番目の基底スペクトルに対応するアクティベーションの時刻ｔでの成分を表す。Ｒ^（ｎ） _ｆは、ｒ^（ｎ） _ｆｒ^（ｎ） _ｆ ^Ｈに相当する、ｎ番目の目的音源の時不変な空間相関行列であり、Ｍ×Ｍの行列パラメータである。α^（ｎ） _ｆは、Σ_ｎα^（ｎ） _ｆ＝１を満たす重みパラメータ（スカラー値）である。したがって、混合分布の各要素分布は、平均０，分散共分散行列が基底ベクトルとアクティベーションとの積和に空間相関行列を掛けたものであるといえる。そして、上記の確率分布に基づき、各パラメータを推定するための尤度関数は、以下の（１７）式のように定義できる。

ただし、ｙは、観測ベクトルの集合を表す。Θは、分布パラメータの集合を表す。Ｆ_Ｎｃ（ｙ_ｆ，ｔ；０，Σ_ｋｗ^（ｎ） _ｆ，ｋｈ^（ｎ） _ｋ，ｔＲ^（ｎ） _ｆ）は、ｙ_ｆ，ｔに関する平均０、分散共分散行列Σ_ｋｗ^（ｎ） _ｆ，ｋｈ^（ｎ） _ｋ，ｔＲ^（ｎ） _ｆの複素正規分布の確率密度関数である。

パラメータ推定部１２は、この目的関数が局所最大となる分布パラメータ、すなわち、モデル化した混合分布が受け取った観測ベクトルの分布に近くなるような分布パラメータを求める。この後、パラメータ推定部１２は、目的音源に対応する要素分布の事後確率を、目的音源に対応するマスクの推定値とすることで、空間的情報と、音源信号の周波数方向の構造との双方を考慮に入れて、マスクと分布パラメータとを推定する。

まず、非特許文献１に記載の方法と同様に、Expectation-Maximization（ＥＭ）アルゴリズムの原理に基づき、対数尤度関数を局所最大化するためのＱ関数は、以下の（１８）式のように定義できる。

ここで、λ^（ｎ） _ｆ，ｔは、ｎ番目の目的音源の事後確率に相当する補助パラメータであり、現在保持している分布パラメータを用いて、Ｅステップにおいて、以下の（１９）式のように計算できる。なお、以下で示す（１９）式〜（２１）式、（２５）式及び（２６）式で表されるパラメータ更新則を十分反復することによって最終的に推定されたλ^（ｎ） _ｆ，ｔは、マスクとして、音源分離部１３に出力される。

そして、分布パラメータのうち、空間相関行列Ｒ^（ｎ） _ｆと、重みパラメータα^（ｎ） _ｆについては、Ｍステップにおいて、上記のＱ関数をそれぞれのパラメータで偏微分したものを０とおくことで、以下の（２０）式及び（２１）式に示す更新則が導出できる。

次に、本実施の形態では、上記Ｑ関数に対して、対数関数とその接線の間に成り立つ不等式と、逆数関数に対するジェンセン（Ｊｅｎｓｅｎ）の不等式とに基づいて設計された下限関数Ｑ´（Θ，Λ）を以下の（２２）式のように設計する。

パラメータ推定部１２は、（２２）式に示す下限関数Ｑ´（Θ，Λ）を用いてパラメータを最適化する。なお、Λは、補助パラメータγ^（ｎ） _ｆ，ｔとβ^（ｎ） _{ｋ，ｆ，ｔ}との集合を表す。そして、（２２）式の等号成立条件は、以下の（２３）式及び（２４）式で表せる。

（２３）式及び（２４）式の条件下において、下限関数Ｑ′（Θ，Λ）の偏微分を０と置くことで、基底スペクトルとアクティベーションとの更新則は、以下の（２５）式及び（２６）式のように導ける。

パラメータ推定部１２は、上記の（１９）式〜（２１）式、（２５）式及び（２６）式で表されるパラメータ更新則を十分反復することで、観測分布を近似するような分布パラメータ（基底ベクトル、アクティベーション及び空間相関行列）及びマスクを推定することができる。パラメータ推定部１２は、観測分布を近似するような分布パラメータ（基底ベクトル、アクティベーション及び空間相関行列）及びマスクを、音源分離部１３に出力する。

［音源分離部の処理］
次に、音源分離部１３の処理について説明する。音源分離部１３は、以下に示すいずれかの手法に基づいて、パラメータ推定部１２が推定したパラメータを用いて音源を分離する。

まず、第１の手法について説明する。第１の手法として、音源分離部１３は、推定したマスクλ^（ｎ） _ｆ，ｔを、観測信号に掛け合わせることによって、ｎ番目の目的音源に対応する分離音を得る。

そして、第２の手法として、音源分離部１３は、推定した基底スペクトルとアクティベーションとを用いて、単チャンネルウィーナーフィルタを設計し、観測信号に掛け合わせることによって、ｎ番目の目的音源に対応する分離音を得る。

また、第３の手法として、音源分離部１３は、推定した基底スペクトル、アクティベーション、空間相関行列、或いは、空間相関行列のみを用いて、多チャンネルウィーナーフィルタを設計し、観測信号に掛け合わせることによって、ｎ番目の目的音源に対応する分離音を得る。

また、音源分離部１３は、第１の手法を、第２の手法或いは第３の手法と組み合わせた第４の手法を用いて、ｎ番目の目的音源に対応する分離音を得ることもできる。

［音源分離処理］
次に、音源分離装置１が実行する音源分離処理の処理手順について説明する。図２は、図１に示す音源分離装置１が実行する音源分離処理の処理手順を示すフローチャートである。

図２に示すように、まず、音源分離装置１では、時間周波数分析部１１が、Ｎ個の目的音源に対応する音響信号が混在する状況において、異なる位置で収録されたＭ個の観測信号を入力信号として受け取る（ステップＳ１）。言い換えると、時間周波数分析部１１は、多チャンネル観測信号を入力信号として受け取る。そして、時間周波数分析部１１は、各観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出する（ステップＳ２）。続いて、時間周波数分析部１１は、時間周波数点ごとに全ての観測音響信号をまとめた観測ベクトルを構成する（ステップＳ３）。

パラメータ推定部１２は、初期値生成部１４からパラメータの初期値を受け取る（ステップＳ４）。パラメータ推定部１２は、時間周波数分析部１１から観測ベクトルを受け取る。そして、パラメータ推定部１２は、各周波数における観測ベクトルの確率分布を、Ｎ個の音源に対応するＮ個の要素分布からなる混合分布でモデル化する（ステップＳ５）。

続いて、パラメータ推定部１２は、モデル化した混合分布が、時間周波数分析部１１から受け取った観測ベクトルの分布に近くなるように混合分布のパラメータを推定する分布パラメータ推定処理を行う（ステップＳ６）。パラメータ推定部１２は、各要素分布の事後確率を各音響信号に対応するマスクの推定値とするマスク推定処理を行う（ステップＳ７）。なお、実際には、ステップＳ６及びステップＳ７は、所定回反復される。すなわち、ステップＳ６及びステップＳ７では、上述したパラメータ更新則を反復することによって、観測分布を近似するような分布パラメータ（基底ベクトル、アクティベーション及び空間相関行列）及びマスクを推定している。

そして、音源分離部１３は、観測ベクトルと、マスクと、分布パラメータ、或いは、その一部をパラメータ推定部１２から受け取り、分離音を計算し（ステップＳ８）、計算した分離音を出力する。分布パラメータは、基底ベクトル、アクティベーション及び空間相関行列である。

［確認実験］
本実施の形態に対し、確認実験を行った。残響時間０．４４秒の可変残響室において測定した室内インパルス応答に、２つの音声をそれぞれ畳み込み足し合わせることで、擬似的に、Ｍ＝２個のマイクで収録された、混合音を得た。この混合音に対して、以下の（表１）の条件下で、本実施の形態に係る音源分離方法を適用し、分離音を得た。

本実施の形態を適用して音源分離を行った結果、非特許文献１記載の方法を用いて音源分離を行った場合に対して、信号対ひずみ比が０．２６ｄＢ向上し、信号対干渉比が１．９６ｄＢ向上した。なお、本実施の形態における、マスクλ^（ｎ） _ｆ，ｔと、空間相関行列Ｒ^（ｎ） _ｆの初期値は、非特許文献１に記載の方法を用いて求めた。基底スペクトルの初期値は、乱数を用いて設定し、アクティベーションの初期値は、全ての成分を１とした。また、簡単のため、重みパラメータα^（ｎ） _ｆの値は、全てのｎ，ｆについて一様とした。

［実施の形態の効果］
このように、本実施の形態によれば、マスクや基底スペクトル、アクティベーションなどの、音源分離のためのパラメータを正確に推定することによって、音源分離を正確に行うことができる。また、本実施の形態によれば、音源信号の周波数方向の構造を考慮することによって、より正確にマスクを推定することができる。さらに、本実施の形態によれば、空間的情報に基づいて推定されたマスクを考慮することによって、より正確に基底ベクトルやアクティベーションを推定して、音源分離を正確に行うことができる。

したがって、本実施の形態によれば、確認実験の結果でも示したように、従来の方法と比して、信号対ひずみ比及び信号対干渉比を高めた音源分離、すなわち、正確な音源分離を行うことが可能になる。

［他の実施の形態］
［システム構成等］
図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行なわれる各処理機能は、その全部又は任意の一部が、ＣＰＵ（Central Processing Unit）及び当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。

また、本実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的におこなうこともでき、あるいは、手動的におこなわれるものとして説明した処理の全部又は一部を公知の方法で自動的におこなうこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。

［プログラム］
図３は、プログラムが実行されることにより、音源分離装置１が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。

メモリ１０１０は、ＲＯＭ（Read Only Memory）１０１１及びＲＡＭ（Random Access Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic Input Output System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１１００に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１１００に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１１１０、キーボード１１２０に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１１３０に接続される。

ハードディスクドライブ１０９０は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、音源分離装置１の各処理を規定するプログラムは、コンピュータ１０００により実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、音源分離装置１における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤ（Solid State Drive）により代替されてもよい。

また、上述した実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。

なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３及びプログラムデータ１０９４は、ネットワーク（ＬＡＮ、ＷＡＮ等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３及びプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。

以上、本発明者によってなされた発明を適用した実施形態について説明したが、本実施形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。すなわち、本実施形態に基づいて当業者等によりなされる他の実施形態、実施例及び運用技術等は全て本発明の範疇に含まれる。

１，１Ｐ，１Ｐ’ 音源分離装置
１１，１１Ｐ，１１Ｐ’ 時間周波数分析部
１２，１２Ｐ，１２Ｐ’ パラメータ推定部
１３，１３Ｐ，１３Ｐ’ 音源分離部
１４，１４Ｐ，１４Ｐ’ 初期値生成部
１２１，１２１Ｐ 分布パラメータ推定部
１２２，１２２Ｐ マスク推定部
１２１Ｐ’ パラメータ推定処理部

Claims (7)

1. Ｎ（Ｎ＞１）個の目的音源に対応する音響信号が混在する状況において、異なる位置で収録されたＭ（Ｍ＞１）個の観測信号を受け取り、各観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出した後に時間周波数点ごとに全ての観測音響信号をまとめたＭ次元縦ベクトルである観測ベクトルを構成する時間周波数分析部と、
   Ｎ個の音源に対応するＮ個の要素分布の各々は、分布パラメータのうち各音源のパワーに相当するパラメータが周波数方向に構造を持った基底ベクトルとその時変であるアクティベーションとの積和によってモデル化されるものとし、各周波数における観測ベクトルの確率分布は、前記Ｎ個の要素分布の混合分布でモデル化されるものとし、前記時間周波数分析部から観測ベクトルを受け取り、モデル化した前記混合分布が、前記時間周波数分析部から受け取った観測ベクトルの分布に近くなるように前記混合分布のパラメータを推定するパラメータ推定部と、
   前記推定したパラメータを用いて、分離音を計算する音源分離部と、
   を有することを特徴とする音源分離装置。
2. 前記混合分布の各要素分布は、平均を０とし，分散共分散行列を基底ベクトルとアクティベーションとの積和に空間相関行列を掛けたものとするガウス分布であることを特徴とする請求項１に記載の音源分離装置。
3. 前記パラメータ推定部は、対数関数とその接線の間に成り立つ不等式と、逆数関数に対するジェンセン（Ｊｅｎｓｅｎ）の不等式とに基づいて設計された下限関数を用いてパラメータを最適化することを特徴とする請求項１または２に記載の音源分離装置。
4. 目的音源を分離する音源分離装置が行う音源分離方法であって、
   Ｎ（Ｎ＞１）個の目的音源に対応する音響信号が混在する状況において、異なる位置で収録されたＭ（Ｍ＞１）個の観測信号を受け取り、各観測信号に短時間信号分析を適用して時間周波数点ごとの観測音響信号を抽出した後に時間周波数点ごとに全ての観測音響信号をまとめたＭ次元縦ベクトルである観測ベクトルを構成する時間周波数分析工程と、
   Ｎ個の音源に対応するＮ個の要素分布の各々は、分布パラメータのうち各音源のパワーに相当するパラメータが周波数方向に構造を持った基底ベクトルとその時変であるアクティベーションとの積和によってモデル化されるものとし、各周波数における観測ベクトルの確率分布は、前記Ｎ個の要素分布の混合分布でモデル化されるものとし、前記時間周波数分析工程における観測ベクトルを受け取り、モデル化した前記混合分布が、前記観測ベクトルの分布に近くなるように前記混合分布のパラメータを推定するパラメータ推定工程と、
   前記推定されたパラメータを用いて、分離音を計算する音源分離工程と、
   を含んだことを特徴とする音源分離方法。
5. 前記混合分布の各要素分布は、平均を０とし，分散共分散行列を基底ベクトルとアクティベーションとの積和に空間相関行列を掛けたものとするガウス分布であることを特徴とする請求項４に記載の音源分離方法。
6. 前記パラメータ推定工程は、対数関数とその接線の間に成り立つ不等式と、逆数関数に対するジェンセンの不等式とに基づいて設計された下限関数を用いてパラメータを最適化することを特徴とする請求項４または５に記載の音源分離方法。
7. コンピュータを請求項１〜３のいずれか一つに記載の音源分離装置として機能させるための音源分離プログラム。

Images