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

Description

本発明は信号処理の技術分野に属する。特に本発明は複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を推定し、分離抽出する音源分離技術に関する。特に、原信号やそれらがどのように混ざったかの情報を用いずに複数の原信号とノイズとが混在している観測信号のみから、それぞれの原信号を推定する、ブラインド音源分離技術に属する。

非特許文献１が音源分離の従来技術として知られている。非特許文献１では、音源ｋから発せられる原信号の、二個のマイクへの到達時間差δ_ｋを

として推定している。

和泉洋介，小野順貴，嵯峨山茂樹，"スパースな混合モデルに基づく雑音・残響環境下の劣決定ブラインド音源分離"，電子情報通信学会総合大会講演論文集，２００８年３月

しかしながら、従来技術では、到達時間差δ_ｋの解析的な更新式は与えられていないため、多くの計算コストを要する全探索操作によって、到達時間差δ_ｋを推定する必要がある。

本発明は、到達時間差δの性質に着目し、到達時間差δを効率的に推定する方法を与え、従来技術において必要であった全探索操作を不要とし、高速な音源分離技術を提供することを目的とする。

上記の課題を解決するために、本発明の第一の態様によれば、複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する。原信号の音源のインデックスをｋとし、周波数領域の観測信号ｘ_ｆ，ｔ＝［ｘ_{ｆ，ｔ，Ｌ}，ｘ_{ｆ，ｔ，Ｒ}］^Ｔから、二個のマイクへの原信号の到達時間差δ_ｆ，ｋと雑音のパワースペクトルσ_ｆ ^２と原信号のスペクトルｓ_{ｆ，ｔ，ｋ}と音源存在確率ｐ（ｋ｜θ）とを推定する。観測信号ｘ_ｆ，ｔと推定されたパラメタθ＝｛δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}，ｐ（ｋ｜θ）｝とから分離信号ｙ_{ｆ，ｔ，ｋ}を生成する。観測信号ｘ_ｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍ_{ｆ，ｔ，ｋ}とし、二個のマイクの間隔をＤとし、原信号の速度をｃとし、φ_ｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξ_{ｆ，ｔ，ｋ}＝［ｘ_{ｆ，ｔ，Ｒ}−φ_ｆ（ｘ_{ｔ，ｔ，Ｌ}−ｓ_{ｆ，ｔ，ｋ}）］とし、スペクトルｓ_{ｆ，ｔ，ｋ}及びξ_{ｆ，ｔ，ｋ}の位相をそれぞれψ_ｓｋ及びψ_ξｋとし、到達時間差δ_ｆ，ｋを

として推定する。推定された到達時間差δ_ｆ，ｋが内包する±πの不定性を補正する。

本発明は、到達時間差δを効率的に推定する方法を与え、高速な音源分離ができるという効果を奏する。

第一実施形態に係る音源分離装置の機能ブロック図。 第一実施形態に係る音源分離装置の処理フローを示す図。 Ｅステップ計算部の機能ブロック図。 Ｍステップ計算部の機能ブロック図。 パラメタ推定部の処理フローを示す図。 シミュレーションを行った環境を示す図。 図７Ａは音源数が２つの場合のシミュレーション結果を、図７Ｂは音源数が３つの場合のシミュレーション結果を示す図。 音源スペクトル推定部の処理を時間差推定部の処理の前に行う場合のＭステップ計算部の機能ブロック図。

以下、本発明の実施形態について説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。

＜第一実施形態に係る音源分離装置１＞
図１は本実施形態に係る音源分離装置１の機能ブロック図を、図２はその処理フローを示す。

音源分離装置１は周波数領域変換部１１とパラメタ推定部１２と分離信号生成部１３と時間領域変換部１４とを含み、さらにパラメタ推定部１２はＥステップ計算部１２１とＭステップ計算部１２２とを含む。

まず、観測信号について説明する。信号原（以下「音源」とも言う）がＫ個あり、ｋ（ｋ＝１，２，…，Ｋ）を音源のインデックスとし、音源ｋから発せられる信号（以下「原信号」という）をｓ_ｋ（ｔ）とする。複数の原信号ｓ_１（ｔ），・・・，ｓ_Ｋ（ｔ）がノイズとともに二個のマイクＬ，Ｒで観測される状況で、マイクＬで観測される時間領域の観測信号をｘ_Ｌ（ｔ）とし、マイクＲで観測される時間領域の観測信号をｘ_Ｒ（ｔ）とし、二個のマイクＬ，Ｒで観測される時間領域の観測信号をｘ（ｔ）＝［ｘ_Ｌ（ｔ），ｘ_Ｒ（ｔ）］^Ｔとする。「^Ｔ」は転置を表す。ｔはフレーム番号及びそのフレーム番号に対応する時刻を表す。Ｔを時間フレームの総数とすると、ｔ＝０，１，…，Ｔ−１である。ここで周波数領域の観測信号をｘ_ｆ，ｔ＝［ｘ_{ｆ，ｔ，Ｌ}，ｘ_{ｆ，ｔ，Ｒ}］^Ｔと表記する。なお、ｆはサンプリング周波数ｆ_ｓをＦ等分した離散点であり、ｆ∈｛０，ｆ_ｓ／Ｆ，…，（Ｆ−１）ｆ_ｓ／Ｆ｝である。以降、断りのない場合、観測信号とは、周波数領域の観測信号ｘ_ｆ，ｔ＝［ｘ_{ｆ，ｔ，Ｌ}，ｘ_{ｆ，ｔ，Ｒ}］^Ｔを指し、時間領域の観測信号の場合はそれを明記する。
ここで、観測信号は、

と表されると仮定する。ここで、ｓ_{ｆ，ｔ，ｋ}は原信号ｓ_ｋ（ｔ）のスペクトル（以下「音源スペクトル」ともいう）を、ｎ_ｆ，ｋ＝［ｎ_{ｆ，ｋ，Ｌ}，ｎ_{ｆ，ｋ，Ｒ}］はマイクＬ，Ｒにおける加算的雑音を表す。またｂ_ｆ，ｋ＝［ｂ_{ｆ，ｋ，Ｌ}，ｂ_{ｆ，ｋ，Ｒ}］^Ｔは音源ｋに関するステアリングベクトル（音源ｋの方向を特定するベクトルであり、方向ベクトルともいう）であり、原信号ｓ_ｋ（ｔ）のマイクＬ、Ｒへの到達時間差をδ_ｋとすると、

である（非特許文献１参照）。本実施形態では、原信号の観測時間内においては、音源及びマイクは固定されており、またＫ個の音源は全て、異なる位置に配置されているとする。すなわち、ステアリングベクトルｂ_ｆ，ｋは時間ｔに依らず、ｋの値によって異なる値を取るものと仮定する。音源分離の目的は、観測信号ｘ_ｆ，ｔのみを用いて、全ての音源スペクトルｓ_{ｆ，ｔ，ｋ}を推定することである。

本実施形態に係る音源分離装置１は、時間領域の観測信号ｘ（ｔ）＝［ｘ_Ｌ（ｔ），ｘ_Ｒ（ｔ）］^Ｔを入力とし、時間領域の分離信号（推定された各原信号）ｙ_ｋ（ｔ）を出力する。以下、各部の処理内容を説明する。

＜周波数領域変換部１１＞
まず、周波数領域変換部１１は、マイクＬ、Ｒで収音した時間領域の観測信号ｘ（ｔ）＝［ｘ_Ｌ（ｔ），ｘ_Ｒ（ｔ）］^Ｔを入力とし、これを短時間フーリエ変換等により周波数領域の観測信号ｘ_ｆ，ｔ＝［ｘ_{ｆ，ｔ，Ｌ}，ｘ_{ｆ，ｔ，Ｒ}］^Ｔに変換し（ｓ１）、パラメタ推定部１２及び分離信号生成部１３に出力する。

＜パラメタ推定部１２＞
次に、パラメタ推定部１２は、観測信号ｘ_ｆ，ｔから、音源分離のために必要なパラメタθを推定する（ｓ２）。なおθ＝｛δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}，ｐ（ｋ｜θ）｝であり、δ_ｆ，ｋは原信号の、二個のマイクへの到達時間差を表す。σ_ｆ ^２は雑音のパワースペクトルを、ｓ_{ｆ，ｔ，ｋ}は音源スペクトルを、ｐ（ｋ｜θ）は音源存在確率（混合信号中の音源ｋの寄与率）を表す。

本実施形態では、上記パラメタを推定するために、以下の２つの仮定を用いる。

（仮定１）は、雑音ｎが、平均０、共分散行列σ_ｆ ^２Ｖ_ｆの正規分布に従う定常雑音でモデル化できるという仮定である。ここでσ_ｆ ^２は周波数ｆにおける雑音のパワーであり、Ｖ_ｆは例えば拡散性雑音の場合

で与えられる（非特許文献１参照）。ここでｃは原信号の速度（音速等）、Ｄは二個のマイクの間隔である。

（仮定２）は、原信号がスパースな信号（つまり、成分のうち０でないもの（非零の成分）がまばらである信号、言い換えると、ほとんどの（または多くの）成分が０である信号）であるという仮定である。すなわち、ある時間周波数（ｆ，ｔ）において、たかだか１つの原信号のみが支配的であると仮定する。（仮定２）によると、式（２）は

と表記できる（非特許文献１参照）。
上記（仮定１）、（仮定２）に基づけば、到達時間差δ_ｆ，ｋとなる方向に存在する音源ｋから発せられる原信号が到来してｘ_ｆ，ｔが観測される尤度は、

で与えられる（非特許文献１参照）。「^Ｈ」はエルミート転置を表す。ここでｂ_ｆ，ｋは式（３）で表される。これより、到達時間差δ_ｆ，ｋ、音源スペクトルｓ_{ｆ，ｔ，ｋ}及び雑音のパワースペクトルσ_ｆ ^２は、対数尤度関数

を最大化するパラメタとして最尤推定により求める（非特許文献１参照）。但し、上記において、ｐ（ｋ｜δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}）は音源存在確率を表し、混合信号中の音源ｋの寄与率であり、

である（非特許文献１参照）。
具体的には、期待値最大化法（以下「ＥＭアルゴリズム」ともいう）を適用し、Ｑ関数

を最大とするパラメタθ＝｛δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}，ｐ（ｋ｜θ）｝を、以下のＥステップ及びＭステップの繰り返しにより求める（非特許文献１参照）。ここでｍ_{ｆ，ｔ，ｋ}は後述する式（１１）で、ｐ（ｘ_ｆ，ｔ｜ｋ，θ’）は前述の式（６）で与えられ、θ’は、１回前の繰り返しで得られているパラメタを意味する。

この繰り返し計算を、パラメタ推定部１２にて行なう。以下、パラメタ推定部１２の処理の詳細を説明する。

図３はＥステップ計算部１２１の機能ブロック図を、図４はＭステップ計算部１２２の機能ブロック図を、図５はパラメタ推定部１２の処理フローを示す。

まず、パラメタ推定部１２は、各パラメタを初期化する（ｓ２０）。パラメタのうちσ_ｆ ^２、δ_ｆ，ｋ、及びｐ（ｋ｜θ）を初期化する。例えば、σ_ｆ ^２＝｜ｘ_{ｆ，ｔ＝０，Ｌ}｜^２，δ_ｆ，ｋ＝（Ｄ／ｃ）ｃｏｓα_ｋ（ｃは原信号の速度（音速等）、Ｄはマイク間隔，α_ｋは音源ｋの方向の初期値（−π／２〜π／２の間の適当な値）），ｐ（ｋ｜θ）＝１／Ｋとする。さらに、初期化したパラメタδ_ｆ，ｋとｘ_{ｆ，ｔ＝０}とを用いて、式（３）及び後述する式（１９）に基づき、ｓ_{ｆ，ｔ，ｋ}を初期化する。更新回数ｎ＝０とする。なお、最大更新回数Ｎ及び収束判定閾値Δは、当該装置の設計者や利用者等により予め設定されているものとする。

Ｅステップ計算部１２１は事後確率推定部１２１１とＱ関数計算部１２１２とを含む（図３参照）。Ｍステップ計算部１２２は時間差推定部１２２１と音源スペクトル推定部１２２２と雑音パワー推定部１２２３と音源存在確率推定部１２２４とを含み、時間差推定部１２２１はさらに逆正接計算部１２２１１と時間補正部１２２１２とを備える（図４参照）。事後確率推定部１２１１とＱ関数計算部１２１２とにおける処理を併せてＥステップと呼び、時間差推定部１２２１と音源スペクトル推定部１２２２と雑音パワー推定部１２２３と音源存在確率推定部１２２４とにおける処理を併せてＭステップと呼ぶ。

（事後確率推定部１２１１）
Ｅステップ計算部１２１の事後確率推定部１２１１は、観測信号ｘ_ｆ，ｔと、一回前の繰り返しで得られているパラメタθ’（但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の事後確率推定においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて事後確率ｍ_{ｆ，ｔ，ｋ}＝ｐ（ｋ｜ｘ_ｆ，ｔ，θ’）を以下の式（１１）により求め（ｓ２２）、Ｑ関数計算部１２１２とＭステップ計算部１２２とに出力する。

なお、ｐ（ｘ_ｆ，ｔ｜ｋ，θ’）は式（６）により与えられる。なお、事後確率ｍ_{ｆ，ｔ，ｋ}は観測信号ｘ_ｆ，ｔが音源ｋに帰属する事後確率を表す。

（逆正接計算部１２２１１）
Ｍステップ計算部１２２の時間差推定部１２２１の逆正接計算部１２２１１は、観測信号ｘ_ｆ，ｔと、事後確率ｍ_{ｆ，ｔ，ｋ}と、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの音源スペクトルｓ_{ｆ，ｔ，ｋ}である。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の逆正接計算においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、到達時間差δ_ｆ，ｋ（より詳しく言うと到達時間差δ_ｆ，ｋに２πｆを乗じた値２πｆδ_ｆ，ｋ）を以下の式（１３）により推定し（ｓ２５）、時間補正部１２２１２に出力する。

ψ_ｓｋ（但し添え字ｓｋはｓ_ｋを表す）及びψ_ξｋ（但し添え字ξｋはξ_ｋを表す）は、それぞれｓ_{ｆ，ｔ，ｋ}及びξ_{ｆ，ｔ，ｋ}の位相を表す。なお、式（１３）の導出については後述する。

なお、従来技術では式（１）に基づきδ_ｋを離散全探索によって推定するため多くの計算コストを要していたが、本実施形態では式（１３）に基づきδ_ｆ，ｋを算出するため全探索を要せず計算コストを小さくできる。また、式（１３）に基づき周波数ｆ毎に到達時間差δ_ｆ，ｋを求める点が従来技術とは異なる。

（時間補正部１２２１２）
Ｍステップ計算部１２２の時間差推定部１２２１の時間補正部１２２１２は、観測信号ｘ_ｆ，ｔと、事後確率ｍ_{ｆ，ｔ，ｋ}と、一回前の繰り返しで得られているパラメタθ’（但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、式（１３）にて推定した到達時間差δ_ｆ，ｋが内包する±πの不定性を補正する（ｓ２６）。この補正は、式（１３）の左辺２πｆδ_ｆ，ｋが−πからπの値を取るのに対し、式（１３）の右辺の逆正接が−π／２からπ／２の値しか返すことができないため、２πｆδ_ｆ，ｋが−π〜−π／２及びπ／２〜πの範囲を取る場合の値を正しく求めるために必要である。式（１３）で得られた値を２πｆδ’_ｆ，ｋと記載すると、補正は以下のように行なう。

ここでξ_{ｆ，ｔ，ｋ}、φ_ｆはそれぞれ式（１４）、（１５）で与えられ、ψ_ｓｋ（但し添え字ｓｋはｓ_ｋを表す）、ψ_ξｋ（但し添え字ξｋはξ_ｋを表す）はそれぞれｓ_{ｆ，ｔ，ｋ}及びξ_{ｆ，ｔ，ｋ}の位相を表す。なお、式（１７）、（１８）の導出については後述する。

さらに、時間補正部１２２１２は、補正した値２πｆδ_ｆ，ｋを２πｆで除算し、到達時間差δ_ｆ，ｋを求め、音源スペクトル推定部１２２２と雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。

（音源スペクトル推定部１２２２）
Ｍステップ計算部１２２の音源スペクトル推定部１２２２は、到達時間差δ_ｆ，ｋと観測信号ｘ_ｆ，ｔを入力とし、これらの値を用いて、音源スペクトルｓ_{ｆ，ｔ，ｋ}を以下の式（１９）により推定し（ｓ２７）、雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。

ここで、ｂ_ｆ，ｋ、Ｖ_ｆはそれぞれ式（３）、（４）で表される。

（雑音パワー推定部１２２３）
Ｍステップ計算部１２２の雑音パワー推定部１２２３は、到達時間差δ_ｆ，ｋと音源スペクトルｓ_{ｆ，ｔ，ｋ}と観測信号ｘ_ｆ，ｔと事後確率ｍ_{ｆ，ｔ，ｋ}とを入力とし、これらの値を用いて、雑音のパワースペクトルσ_ｆ ^２を以下の式（２０）により推定し（ｓ２８）、Ｑ関数計算部１２１２に出力する。

なお、Ｔは時間フレームの総数である。

（音源存在確率推定部１２２４）
Ｍステップ計算部１２２の音源存在確率推定部１２２４は、事後確率ｍ_{ｆ，ｔ，ｋ}を入力とし、音源存在確率を以下の式（２１）により推定し（ｓ２９）、Ｑ関数計算部１２１２に出力する。

（Ｑ関数計算部１２１２）
Ｅステップ計算部１２１のＱ関数計算部１２１２は、観測信号ｘ_ｆ，ｔと、パラメタθ＝｛δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}，ｐ（ｋ｜θ）｝と、事後確率ｍ_{ｆ，ｔ，ｋ}と、を入力とし、これらの値を用いてＱ関数を上述の式（１０）により求める（ｓ３０）。

ｓ２０〜ｓ３０までの処理を終えると、パラメタ推定部１２は、（条件１）Ｑ関数の値の変化量（｜Ｑ（θ｜θ^ｎ−１）−Ｑ（θ｜θ^ｎ）｜）が所定の収束判定閾値Δより小さくなるか、または、（条件２）更新回数ｎが所定の最大更新回数Ｎ（例えばＮ＝２０）以上か否かを判定する（ｓ３１）。

パラメタ推定部１２は、（条件１）、（条件２）の何れかを満たしたときは、パラメタ推定部１２はその時点で取得している最新の事後確率ｍ_{ｆ，ｔ，ｋ}と到達時間差δ_ｆ，ｋを分離信号生成部１３に出力する。

パラメタ推定部１２は、（条件１）、（条件２）の何れも満たさないときは、ＥステップとＭステップを繰り返す（ｓ３１、ｓ２１）。なお、図示しない記憶部にパラメタθとＱ関数の値Ｑ（θ｜θ^ｎ）とを記憶しておき、次の繰り返しの際に用いる。

＜分離信号生成部１３＞
分離信号生成部１３は、事後確率ｍ_{ｆ，ｔ，ｋ}と到達時間差δ_ｆ，ｋと観測信号ｘ_ｆ，ｔとを入力とし、以下の式（２２）により、分離信号ｙ_{ｆ，ｔ，ｋ}を生成し（ｓ３）、時間領域変換部１４へ出力する。

なお、音源スペクトルｓ_{ｆ，ｔ，ｋ}は、音源の発する原信号のスペクトルを推定したものであるが、この音源スペクトルを単純に時間領域の信号に変換した場合には、他の音源の発する原信号が残ることがある。上述の式（２２）によって、音源ｋの発する原信号のみを抽出、分離することができる。

＜時間領域変換部１４＞
時間領域変換部１４は、分離信号ｙ_{ｆ，ｔ，ｋ}を入力とし、周波数領域変換部１１において行った周波数領域変換方法に対応する時間領域変換方法（例えば短時間フーリエ逆変換）で、分離信号ｙ_{ｆ，ｔ，ｋ}を時間領域の分離信号ｙ_ｋ（ｔ）に変換し（ｓ４）、音源分離装置１の出力値として出力する。

＜本実施形態のポイント＞
以下、本実施形態のポイントを説明し、式（１３）、（１７）、（１８）の導出方法を説明する。

本実施形態では、
（性質１）到達時間差δ_ｆ，ｋがＲチャネルとＬチャネルの位相差に影響を与える値であること
（性質２）位相差が周期的な値を取る量であること
という２つの性質を利用して、到着時間差δ_ｆ，ｋを推定する。ここで（性質１）は、式（２）にてノイズｎ_ｆ，ｔが十分に小さい場合を考えれば明らかである。（性質２）は、ある２つの位相を表わす量の差Θが、−π≦Θ＜πの範囲の値だけではなく、Θ±２πＭ（Ｍは任意の整数）という周期的な不定性を内包する値を取る性質を持つことを意味する。

（性質２）について、さらに以下にて説明を行なう。式（６）

の右辺における、ｅｘｐのカッコの中のベクトル及び行列ｘ，ｂ，Ｖを、それぞれの成分で表して整理すると、

となる（Ｃはδ_ｆ，ｋに依らない定数）。式（３２）は、位相や角度の分布のように周期的な値を取る変数に対する分布であるＶｏｎ Ｍｉｓｅｓ分布

と同じ形をしていることが分かる（参考文献１参照）。
［参考文献１］Ｃ．Ｍ．ビショップ著，元田ら訳，“パターン認識と機械学習（上）”，シュプリンガー・ジャパン，２００６．

ここで−π＜ｘ≦π、μは分布の平均（−π＜μ≦π）、κ＞０は分布の集中度パラメタ（正規分布での（１／分散）に相当）、Ｉ_０（ｘ）は０次の第１種ベッセル関数である。

すなわち、式（３３）の変数ｘが式（３２）のψ_ξｋ−ψ_ｓｋに対応し、式（３３）の平均μが式（３２）の２πｆδ_ｆ，ｋに対応し、式（３３）の集中度κが式（３２）の（｜ξ_{ｆ，ｔ，ｋ}｜｜ｓ_{ｆ，ｔ，ｋ}｜）／（σ_ｆ ^２（１−φ_ｆ ^２））に対応する。

説明をより直感的にするため、雑音のパラメタがφ_ｆ＝０である場合を考える（これは、雑音が式（４）に示す拡散性雑音ではなく、分散σ_ｆ ^２のガウス雑音が観測ｘ_Ｌ，ｘ_Ｒにそれぞれ乗ることを意味する）。このとき、式（１４）によりξ_{ｆ，ｔ，ｋ}＝ｘ_{ｆ，ｔ，Ｒ}となるため、式（３２）において、ψ_ξｋはψ_{ｘｆ，ｔ，Ｒ}（但し、添え字ｘｆ，ｔ，Ｒはｘ_{ｆ，ｔ，Ｒ}を表し、ψ_{ｘｆ，ｔ，Ｒ}はマイクＲの観測信号ｘ_{ｆ，ｔ，Ｒ}の位相を表す）となる。また、式（３）、（５）により（但し式（５）においてｎ_ｆ，ｔ＝０）、ｓ_{ｆ，ｔ，ｋ}＝ｘ_{ｆ，ｔ，Ｌ}となるため、式（３２）においてψ_ｓｋはψ_{ｘｆ，ｔ，Ｌ}（但し、添え字ｘｆ，ｔ，Ｌはｘ_{ｆ，ｔ，Ｌ}を表し、ψ_{ｘｆ，ｔ，Ｌ}はマイクＬの観測信号ｘ_{ｆ，ｔ，Ｌ}の位相を表す）となる。また、前述の通りξ_{ｆ，ｔ，ｋ}＝ｘ_{ｆ，ｔ，Ｒ}であり、式（３）、（５）よりｘ_{ｆ，ｔ，Ｒ}＝ｅ^{ｊ２πｆδｋ，ｆ}・ｓ_{ｆ，ｔ，ｋ}（但し、ｅの添え字δｋ，ｆはδ_ｋ，ｆを表す）となるため、式（３２）において、｜ξ_{ｆ，ｔ，ｋ}｜＝｜ｅ^{ｊ２πｆδｋ，ｆ}・ｓ_{ｆ，ｔ，ｋ}｜＝｜ｓ_{ｆ，ｔ，ｋ}｜となる。よって、式（３２）は

となる。これとＶｏｎ Ｍｉｓｅｓ分布との解釈を合わせると、
（１）ＲチャネルとＬチャネルの位相差ψ_{ｘｆ，ｔ，Ｒ}−ψ_{ｘｆ，ｔ，Ｌ}という周期的な値を取る変数の分布が、平均μ＝２πｆδ_ｆ，ｋを取る、
（２）Ｖｏｎ Ｍｉｓｅｓ分布の集中度κが、ＳＮ比｜ｓ_{ｆ，ｋ，ｔ}｜^２／σ_ｆ ^２に対応するようになる。すなわち、ＳＮ比が低い（条件が悪い）と、ＲチャネルとＬチャネルの位相差の値の集中度が下がる（＝分散が大きくなる）。これは、ＳＮ比が低い条件においては位相差の測定値がばらつく現象と対応する、
の２点が言える。

本実施形態では、（性質１）、（性質２）を利用しており、実施の手続きとしては、ＲチャネルとＬチャネルの位相差に関する量が、周期的な値を取る変数に対する分布であるＶｏｎ Ｍｉｓｅｓ分布で表現できることに着目し、Ｖｏｎ Ｍｉｓｅｓ分布のパラメタに対する最尤推定により、分布の平均値μ＝２πｆδ_ｆ，ｋを推定することで、信号到達時間差パラメタδ_ｆ，ｋを推定する。

式（３１）のｐ（ｘ_ｆ，ｔ｜ｋ，θ）をＱ関数の式（１０）に代入し、

を解くことで、式（１３）が得られる。
また時間補正部１２２１２における関数Ｆの式（１７）、（１８）は、Ｑ関数の２階微分

である。時間補正部１２２１２では、式（１３）で得られた解δ’_ｆ，ｋがＱ関数の極大値・極小値のどちらを与えるかを、Ｆ（２πｆδ’_ｆ，ｋ）の値の正負にて調べ、δ’_ｆ，ｋが極小値を与える場合に、上述の＋πまたは−πの補正を行なう。

＜効果＞
従来法においては、時間差推定部において解析的な更新式が与えられていなかったため、多くの計算コストを要する全探索操作が必要であった。よって、時間差推定部の計算コストを削減し、高速な音源分離手段を提供することが課題である。本実施形態では、到着時間差δ_ｆ，ｋが、マイクＲとマイクＬの位相差に影響を与える値であることと、位相差が周期的な値を取る性質を持つことに着目し、到着時間差δ_ｆ，ｋを推定する。これにより、従来必要であった全探索操作が不要となるため、高速な音源分離手段を提供することが可能となる。

＜シミュレーション結果＞
発明の効果を示すため実験を行なった。図６に示す部屋において、音源数は２つ（７０度と１５０度）又は３つ（３０，７０，１５０度）とした。音源は、英語話者音声を用い、音源数２及び３の場合それぞれにおいて１０通りの音源組合せにて実験を行なった。

雑音としては、平均０、共分散行列σ_ｆ ^２Ｖ_ｆ（Ｖ_ｆは式（４）により与えられる）のガウスノイズをＳＮ比約２５ｄＢにて重畳した。部屋の残響時間は１３０ｍｓ、サンプリング周波数は８ｋＨｚ，短時間フーリエ変換の窓長及びシフト長はそれぞれ６４ｍｓ，１６ｍｓとした。

従来法では到達時間差δ_ｋを推定する際に、音源位置Θ（図６参照）を０度から１８０度まで１度きざみで変化させ、それに対応する１８１種類の到達時間差δ_ｋの値について全探索を行なった。

図７に、ＳＩＲ（信号対雑音比、雑音には他話者音声含む）、ＳＤＲ（信号対歪み比）、及びパソコン（IntelXeon（登録商標）X5650 2.67GHz(6Core)×2CPU）における計算時間を示す。図７Ａは音源が二つの場合（７０度と１５０度）を、図７Ｂは音源が三つの場合（３０度，７０度，１５０度）を示す。ＳＩＲとＳＤＲは大きな値であるほど良い性能であることを示す。数字は、１０通りの音源組合せの平均値である。図７Ａ、Ｂに示す通り、本実施形態は、１／１０程度の計算時間で、従来法とほぼ同程度の分離性能を達成できることが見てとれる。

＜その他の変形例＞
本実施形態のポイントは上述の通り、到達時間差の推定方法である。従って、他の処理やパラメタの推定方法については、上記の実施形態に限定されるものではなく、他の従来技術を用いてもよい。

本実施形態では、各部間で直接データを受け渡しているが、図示しない記憶部を介して、各データを読み書きしてもよい。

また、本実施形態では、雑音パワー推定部１２２３は音源スペクトルｓ_{ｆ，ｔ，ｋ}を音源スペクトル推定部１２２２から取得しているが、到達時間差δ_ｆ，ｋと観測信号ｘ_ｆ，ｔとを用いて式（１９）に基づき雑音パワー推定部１２２３で計算する構成としてもよい。

本実施形態では、拡散性雑音の場合を想定しているが、他の特性を持つ雑音であってもよい。その場合、雑音の特性に応じて、式（４）や式（１５）のｓｉｎｃ（２πｆＤ／ｃ）を適宜変更すればよい。

分離信号生成部１３は、事後確率ｍ_{ｆ，ｔ，ｋ}と音源スペクトルｓ_{ｆ，ｔ，ｋ}とを入力とし、式（２２）に代えて、以下の式により、分離信号ｙ_{ｆ，ｔ，ｋ}を生成してもよい。

この場合、パラメタ推定部１２は、Ｑ関数の値の変化量（｜Ｑ（θ｜θ^ｎ−１）−Ｑ（θ｜θ^ｎ）｜）が所定の収束判定閾値Δより小さくなったとき、または、更新回数ｎが所定の最大更新回数Ｎ以上になったとき（ｓ３１）に取得している最新の事後確率ｍ_{ｆ，ｔ，ｋ}と音源スペクトルｓ_{ｆ，ｔ，ｋ}を分離信号生成部１３に出力する。このような構成により、式（２２）と同様に分離信号ｙ_{ｆ，ｔ，ｋ}を生成することができる（式（１９）参照）。

Ｍステップの計算順序は、本実施形態の計算順序に限らない。例えば、時間差推定部１２２１と音源スペクトル推定部１２２２の計算順序はどちらを先に行ってもよい。図８は音源スペクトル推定部１２２２の音源スペクトル推定処理（ｓ２７）を行った後に、時間差推定部１２２１の到達時間差推定処理（ｓ２５、ｓ２６）を行う場合の機能ブロック図を示す。以下、図８を用いて説明する。

Ｍステップ計算部１２２の音源スペクトル推定部１２２２は、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの到達時間差δ_ｆ，ｋである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の音源スペクトル計算においては、前述の初期化したパラメタ）と観測信号ｘ_ｆ，ｔを入力とし、これらの値を用いて、音源スペクトルｓ_{ｆ，ｔ，ｋ}を式（１９）により推定し、雑音パワー推定部１２２３とＱ関数計算部１２１２と時間差推定部１２２１と時間補正部１２２１２とに出力する。

Ｍステップ計算部１２２の時間差推定部１２２１の逆正接計算部１２２１１は、観測信号ｘ_ｆ，ｔと、事後確率ｍ_{ｆ，ｔ，ｋ}と、音源スペクトルｓ_{ｆ，ｔ，ｋ}とを入力とし、これらの値を用いて、到達時間差δ_ｆ，ｋを式（１３）により推定し、時間補正部１２２１２に出力する。

Ｍステップ計算部１２２の時間差推定部１２２１の時間補正部１２２１２は、観測信号ｘ_ｆ，ｔと、事後確率ｍ_{ｆ，ｔ，ｋ}と、音源スペクトルｓ_{ｆ，ｔ，ｋ}と、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの雑音のパワースペクトルσ_ｆ ^２である。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、式（１３）にて推定した到達時間差δ_ｆ，ｋが内包する±πの不定性を式（１６）に基づき補正する。さらに、時間補正部１２２１２は、補正した値２πｆδ_ｆ，ｋを２πｆで除算し、到達時間差δ_ｆ，ｋを求め、雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。

上述の構成とすることで、第一実施形態と同様の効果を得ることができる。なお、同様に、一回前の繰り返しで得られているパラメタθ’のうちの到達時間差δ_ｆ，ｋに基づき雑音パワー推定部１２２３の雑音パワー推定処理（ｓ２８）を行った後に、時間差推定部１２２１の到達時間差推定処理（ｓ２５、ｓ２６）を行ってもよい。

本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。

＜プログラム及び記録媒体＞
上述した音源分離装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置（各種実施例で図に示した機能構成をもつ装置）として機能させるためのプログラム、またはその処理手順（各実施例で示したもの）の各過程をコンピュータに実行させるためのプログラムを、ＣＤ−ＲＯＭ、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。

Claims (3)

1. 複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離装置であって、
   前記原信号の音源のインデックスをｋとし、周波数領域の前記観測信号ｘ_ｆ，ｔ＝［ｘ_{ｆ，ｔ，Ｌ}，ｘ_{ｆ，ｔ，Ｒ}］^Ｔから、前記二個のマイクへの前記原信号の到達時間差δ_ｆ，ｋと雑音のパワースペクトルσ_ｆ ^２と原信号のスペクトルｓ_{ｆ，ｔ，ｋ}と音源存在確率ｐ（ｋ｜θ）とを推定するパラメタ推定手段と、
   前記観測信号ｘ_ｆ，ｔと推定されたパラメタθ＝｛δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}，ｐ（ｋ｜θ）｝とから分離信号ｙ_{ｆ，ｔ，ｋ}を生成する分離信号生成手段とを含み、
   前記パラメタ推定手段は、
   観測信号ｘ_ｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍ_{ｆ，ｔ，ｋ}とし、前記二個のマイクの間隔をＤとし、前記原信号の速度をｃとし、φ_ｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξ_{ｆ，ｔ，ｋ}＝［ｘ_{ｆ，ｔ，Ｒ}−φ_ｆ（ｘ_{ｔ，ｔ，Ｌ}−ｓ_{ｆ，ｔ，ｋ}）］とし、前記スペクトルｓ_{ｆ，ｔ，ｋ}及び前記ξ_{ｆ，ｔ，ｋ}の位相をそれぞれψ_ｓｋ及びψ_ξｋとし、前記到達時間差δ_ｆ，ｋを
   

   

   として推定する逆正接計算手段と、
   推定された前記到達時間差δ_ｆ，ｋが内包する±πの不定性を補正する補正手段と、
   前記到達時間差δ _ｆ，ｋ と前記観測信号ｘ _ｆ，ｔ とを用いて、前記原信号のスペクトルｓ _{ｆ，ｔ，ｋ} を
   

   

   として推定する音源スペクトル推定手段と、
   前記到達時間差δ _ｆ，ｋ と前記音源スペクトルｓ _{ｆ，ｔ，ｋ} と前記観測信号ｘ _ｆ，ｔ と前記事後確率ｍ _{ｆ，ｔ，ｋ} とを用いて、前記雑音のパワースペクトルσ _ｆ ^２ を
   

   

   として推定する雑音パワー推定手段と、
   前記事後確率ｍ _{ｆ，ｔ，ｋ} を用いて、前記音源存在確率ｐ（ｋ｜θ）を
   

   

   として推定する音源存在確率推定手段とを有し、
   前記分離信号生成手段は、前記分離信号ｙ _{ｆ，ｔ，ｋ} を
   

   

   として生成する、
   音源分離装置。
2. 複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離方法であって、
   前記原信号の音源のインデックスをｋとし、周波数領域の前記観測信号ｘ_ｆ，ｔ＝［ｘ_{ｆ，ｔ，Ｌ}，ｘ_{ｆ，ｔ，Ｒ}］^Ｔから、前記二個のマイクへの前記原信号の到達時間差δ_ｆ，ｋと雑音のパワースペクトルσ_ｆ ^２と原信号のスペクトルｓ_{ｆ，ｔ，ｋ}と音源存在確率ｐ（ｋ｜θ）とを推定するパラメタ推定ステップと、
   前記観測信号ｘ_ｆ，ｔと推定されたパラメタθ＝｛δ_ｆ，ｋ，σ_ｆ ^２，ｓ_{ｆ，ｔ，ｋ}，ｐ（ｋ｜θ）｝とから分離信号ｙ_{ｆ，ｔ，ｋ}を生成する分離信号生成ステップとを含み、
   前記パラメタ推定ステップは、
   観測信号ｘ_ｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍ_{ｆ，ｔ，ｋ}とし、前記二個のマイクの間隔をＤとし、前記原信号の速度をｃとし、φ_ｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξ_{ｆ，ｔ，ｋ}＝［ｘ_{ｆ，ｔ，Ｒ}−φ_ｆ（ｘ_{ｔ，ｔ，Ｌ}−ｓ_{ｆ，ｔ，ｋ}）］とし、前記スペクトルｓ_{ｆ，ｔ，ｋ}及び前記ξ_{ｆ，ｔ，ｋ}の位相をそれぞれψ_ｓｋ及びψ_ξｋとし、前記到達時間差δ_ｆ，ｋを
   

   

   として推定する逆正接計算ステップと、
   推定された前記到達時間差δ_ｆ，ｋが内包する±πの不定性を補正する補正ステップと、
   前記到達時間差δ _ｆ，ｋ と前記観測信号ｘ _ｆ，ｔ とを用いて、前記原信号のスペクトルｓ _{ｆ，ｔ，ｋ} を
   

   

   として推定する音源スペクトル推定ステップと、
   前記到達時間差δ _ｆ，ｋ と前記音源スペクトルｓ _{ｆ，ｔ，ｋ} と前記観測信号ｘ _ｆ，ｔ と前記事後確率ｍ _{ｆ，ｔ，ｋ} とを用いて、前記雑音のパワースペクトルσ _ｆ ^２ を
   

   

   として推定する雑音パワー推定ステップと、
   前記事後確率ｍ _{ｆ，ｔ，ｋ} を用いて、前記音源存在確率ｐ（ｋ｜θ）を
   

   

   として推定する音源存在確率推定ステップとを有し、
   前記分離信号生成ステップは、前記分離信号ｙ _{ｆ，ｔ，ｋ} を
   

   

   として生成する、
   音源分離方法。
3. コンピュータを請求項１記載の音源分離装置として機能させるためのプログラム。

Images