音源分離装置、音源分離方法及びプログラム
【課題】到達時間差δを効率的に推定する方法を与え全探索操作を不要とし、高速な音源分離技術を提供する。
【解決手段】周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と音源スペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定する。観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、二個のマイクの間隔をＤとし、原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、スペクトルｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、到達時間差δｆ，ｋを
として推定する。推定された到達時間差δｆ，ｋが内包する±πの不定性を補正する。
【解決手段】周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と音源スペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定する。観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、二個のマイクの間隔をＤとし、原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、スペクトルｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、到達時間差δｆ，ｋを
として推定する。推定された到達時間差δｆ，ｋが内包する±πの不定性を補正する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は信号処理の技術分野に属する。特に本発明は複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を推定し、分離抽出する音源分離技術に関する。特に、原信号やそれらがどのように混ざったかの情報を用いずに複数の原信号とノイズとが混在している観測信号のみから、それぞれの原信号を推定する、ブラインド音源分離技術に属する。
【背景技術】
【０００２】
非特許文献１が音源分離の従来技術として知られている。非特許文献１では、音源ｋから発せられる原信号の、二個のマイクへの到達時間差δｋを
【０００３】
【数１】
【０００４】
として推定している。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】和泉洋介，小野順貴，嵯峨山茂樹，“スパースな混合モデルに基づく雑音・残響環境下の劣決定ブラインド音源分離”，電子情報通信学会総合大会講演論文集，２００８年３月
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、従来技術では、到達時間差δｋの解析的な更新式は与えられていないため、多くの計算コストを要する全探索操作によって、到達時間差δｋを推定する必要がある。
【０００７】
本発明は、到達時間差δの性質に着目し、到達時間差δを効率的に推定する方法を与え、従来技術において必要であった全探索操作を不要とし、高速な音源分離技術を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の課題を解決するために、本発明の第一の態様によれば、複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する。原信号の音源のインデックスをｋとし、周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、二個のマイクへの原信号の到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と原信号のスペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定する。観測信号ｘｆ，ｔと推定されたパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝とから分離信号ｙｆ，ｔ，ｋを生成する。観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、二個のマイクの間隔をＤとし、原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、スペクトルｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、到達時間差δｆ，ｋを
【０００９】
【数２】
【００１０】
として推定する。推定された到達時間差δｆ，ｋが内包する±πの不定性を補正する。
【発明の効果】
【００１１】
本発明は、到達時間差δを効率的に推定する方法を与え、高速な音源分離ができるという効果を奏する。
【図面の簡単な説明】
【００１２】
【図１】第一実施形態に係る音源分離装置の機能ブロック図。
【図２】第一実施形態に係る音源分離装置の処理フローを示す図。
【図３】Ｅステップ計算部の機能ブロック図。
【図４】Ｍステップ計算部の機能ブロック図。
【図５】パラメタ推定部の処理フローを示す図。
【図６】シミュレーションを行った環境を示す図。
【図７】図７Ａは音源数が２つの場合のシミュレーション結果を、図７Ｂは音源数が３つの場合のシミュレーション結果を示す図。
【図８】音源スペクトル推定部の処理を時間差推定部の処理の前に行う場合のＭステップ計算部の機能ブロック図。
【発明を実施するための形態】
【００１３】
以下、本発明の実施形態について説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
【００１４】
＜第一実施形態に係る音源分離装置１＞
図１は本実施形態に係る音源分離装置１の機能ブロック図を、図２はその処理フローを示す。
【００１５】
音源分離装置１は周波数領域変換部１１とパラメタ推定部１２と分離信号生成部１３と時間領域変換部１４とを含み、さらにパラメタ推定部１２はＥステップ計算部１２１とＭステップ計算部１２２とを含む。
【００１６】
まず、観測信号について説明する。信号原（以下「音源」とも言う）がＫ個あり、ｋ（ｋ＝１，２，…，Ｋ）を音源のインデックスとし、音源ｋから発せられる信号（以下「原信号」という）をｓｋ（ｔ）とする。複数の原信号ｓ１（ｔ），・・・，ｓＫ（ｔ）がノイズとともに二個のマイクＬ，Ｒで観測される状況で、マイクＬで観測される時間領域の観測信号をｘＬ（ｔ）とし、マイクＲで観測される時間領域の観測信号をｘＲ（ｔ）とし、二個のマイクＬ，Ｒで観測される時間領域の観測信号をｘ（ｔ）＝［ｘＬ（ｔ），ｘＲ（ｔ）］Ｔとする。「Ｔ」は転置を表す。ｔはフレーム番号及びそのフレーム番号に対応する時刻を表す。Ｔを時間フレームの総数とすると、ｔ＝０，１，…，Ｔ−１である。ここで周波数領域の観測信号をｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔと表記する。なお、ｆはサンプリング周波数ｆｓをＦ等分した離散点であり、ｆ∈｛０，ｆｓ／Ｆ，…，（Ｆ−１）ｆｓ／Ｆ｝である。以降、断りのない場合、観測信号とは、周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔを指し、時間領域の観測信号の場合はそれを明記する。
ここで、観測信号は、
【００１７】
【数３】
【００１８】
と表されると仮定する。ここで、ｓｆ，ｔ，ｋは原信号ｓｋ（ｔ）のスペクトル（以下「音源スペクトル」ともいう）を、ｎｆ，ｋ＝［ｎｆ，ｋ，Ｌ，ｎｆ，ｋ，Ｒ］はマイクＬ，Ｒにおける加算的雑音を表す。またｂｆ，ｋ＝［ｂｆ，ｋ，Ｌ，ｂｆ，ｋ，Ｒ］Ｔは音源ｋに関するステアリングベクトル（音源ｋの方向を特定するベクトルであり、方向ベクトルともいう）であり、原信号ｓｋ（ｔ）のマイクＬ、Ｒへの到達時間差をδｋとすると、
【００１９】
【数４】
【００２０】
である（非特許文献１参照）。本実施形態では、原信号の観測時間内においては、音源及びマイクは固定されており、またＫ個の音源は全て、異なる位置に配置されているとする。すなわち、ステアリングベクトルｂｆ，ｋは時間ｔに依らず、ｋの値によって異なる値を取るものと仮定する。音源分離の目的は、観測信号ｘｆ，ｔのみを用いて、全ての音源スペクトルｓｆ，ｔ，ｋを推定することである。
【００２１】
本実施形態に係る音源分離装置１は、時間領域の観測信号ｘ（ｔ）＝［ｘＬ（ｔ），ｘＲ（ｔ）］Ｔを入力とし、時間領域の分離信号（推定された各原信号）ｙｋ（ｔ）を出力する。以下、各部の処理内容を説明する。
【００２２】
＜周波数領域変換部１１＞
まず、周波数領域変換部１１は、マイクＬ、Ｒで収音した時間領域の観測信号ｘ（ｔ）＝［ｘＬ（ｔ），ｘＲ（ｔ）］Ｔを入力とし、これを短時間フーリエ変換等により周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔに変換し（ｓ１）、パラメタ推定部１２及び分離信号生成部１３に出力する。
【００２３】
＜パラメタ推定部１２＞
次に、パラメタ推定部１２は、観測信号ｘｆ，ｔから、音源分離のために必要なパラメタθを推定する（ｓ２）。なおθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝であり、δｆ，ｋは原信号の、二個のマイクへの到達時間差を表す。σｆ２は雑音のパワースペクトルを、ｓｆ，ｔ，ｋは音源スペクトルを、ｐ（ｋ｜θ）は音源存在確率（混合信号中の音源ｋの寄与率）を表す。
【００２４】
本実施形態では、上記パラメタを推定するために、以下の２つの仮定を用いる。
【００２５】
（仮定１）は、雑音ｎが、平均０、共分散行列σｆ２Ｖｆの正規分布に従う定常雑音でモデル化できるという仮定である。ここでσｆ２は周波数ｆにおける雑音のパワーであり、Ｖｆは例えば拡散性雑音の場合
【００２６】
【数５】
【００２７】
で与えられる（非特許文献１参照）。ここでｃは原信号の速度（音速等）、Ｄは二個のマイクの間隔である。
【００２８】
（仮定２）は、原信号がスパースな信号（つまり、成分のうち０でないもの（非零の成分）がまばらである信号、言い換えると、ほとんどの（または多くの）成分が０である信号）であるという仮定である。すなわち、ある時間周波数（ｆ，ｔ）において、たかだか１つの原信号のみが支配的であると仮定する。（仮定２）によると、式（２）は
【００２９】
【数６】
【００３０】
と表記できる（非特許文献１参照）。
上記（仮定１）、（仮定２）に基づけば、到達時間差δｆ，ｋとなる方向に存在する音源ｋから発せられる原信号が到来してｘｆ，ｔが観測される尤度は、
【００３１】
【数７】
【００３２】
で与えられる（非特許文献１参照）。「Ｈ」はエルミート転置を表す。ここでｂｆ，ｋは式（３）で表される。これより、到達時間差δｆ，ｋ、音源スペクトルｓｆ，ｔ，ｋ及び雑音のパワースペクトルσｆ２は、対数尤度関数
【００３３】
【数８】
【００３４】
を最大化するパラメタとして最尤推定により求める（非特許文献１参照）。但し、上記において、ｐ（ｋ｜δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ）は音源存在確率を表し、混合信号中の音源ｋの寄与率であり、
【００３５】
【数９】
【００３６】
である（非特許文献１参照）。
具体的には、期待値最大化法（以下「ＥＭアルゴリズム」ともいう）を適用し、Ｑ関数
【００３７】
【数１０】
【００３８】
を最大とするパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝を、以下のＥステップ及びＭステップの繰り返しにより求める（非特許文献１参照）。ここでｍｆ，ｔ，ｋは後述する式（１１）で、ｐ（ｘｆ，ｔ｜ｋ，θ’）は前述の式（６）で与えられ、θ’は、１回前の繰り返しで得られているパラメタを意味する。
【００３９】
この繰り返し計算を、パラメタ推定部１２にて行なう。以下、パラメタ推定部１２の処理の詳細を説明する。
【００４０】
図３はＥステップ計算部１２１の機能ブロック図を、図４はＭステップ計算部１２２の機能ブロック図を、図５はパラメタ推定部１２の処理フローを示す。
【００４１】
まず、パラメタ推定部１２は、各パラメタを初期化する（ｓ２０）。パラメタのうちσｆ２、δｆ，ｋ、及びｐ（ｋ｜θ）を初期化する。例えば、σｆ２＝｜ｘｆ，ｔ＝０，Ｌ｜２，δｆ，ｋ＝（Ｄ／ｃ）ｃｏｓαｋ（ｃは原信号の速度（音速等）、Ｄはマイク間隔，αｋは音源ｋの方向の初期値（−π／２〜π／２の間の適当な値）），ｐ（ｋ｜θ）＝１／Ｋとする。さらに、初期化したパラメタδｆ，ｋとｘｆ，ｔ＝０とを用いて、式（３）及び後述する式（１９）に基づき、ｓｆ，ｔ，ｋを初期化する。更新回数ｎ＝０とする。なお、最大更新回数Ｎ及び収束判定閾値Δは、当該装置の設計者や利用者等により予め設定されているものとする。
【００４２】
Ｅステップ計算部１２１は事後確率推定部１２１１とＱ関数計算部１２１２とを含む（図３参照）。Ｍステップ計算部１２２は時間差推定部１２２１と音源スペクトル推定部１２２２と雑音パワー推定部１２２３と音源存在確率推定部１２２４とを含み、時間差推定部１２２１はさらに逆正接計算部１２２１１と時間補正部１２２１２とを備える（図４参照）。事後確率推定部１２１１とＱ関数計算部１２１２とにおける処理を併せてＥステップと呼び、時間差推定部１２２１と音源スペクトル推定部１２２２と雑音パワー推定部１２２３と音源存在確率推定部１２２４とにおける処理を併せてＭステップと呼ぶ。
【００４３】
（事後確率推定部１２１１）
Ｅステップ計算部１２１の事後確率推定部１２１１は、観測信号ｘｆ，ｔと、一回前の繰り返しで得られているパラメタθ’（但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の事後確率推定においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて事後確率ｍｆ，ｔ，ｋ＝ｐ（ｋ｜ｘｆ，ｔ，θ’）を以下の式（１１）により求め（ｓ２２）、Ｑ関数計算部１２１２とＭステップ計算部１２２とに出力する。
【００４４】
【数１１】
【００４５】
なお、ｐ（ｘｆ，ｔ｜ｋ，θ’）は式（６）により与えられる。なお、事後確率ｍｆ，ｔ，ｋは観測信号ｘｆ，ｔが音源ｋに帰属する事後確率を表す。
【００４６】
（逆正接計算部１２２１１）
Ｍステップ計算部１２２の時間差推定部１２２１の逆正接計算部１２２１１は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの音源スペクトルｓｆ，ｔ，ｋである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の逆正接計算においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、到達時間差δｆ，ｋ（より詳しく言うと到達時間差δｆ，ｋに２πｆを乗じた値２πｆδｆ，ｋ）を以下の式（１３）により推定し（ｓ２５）、時間補正部１２２１２に出力する。
【００４７】
【数１２】
【００４８】
ψｓｋ（但し添え字ｓｋはｓｋを表す）及びψξｋ（但し添え字ξｋはξｋを表す）は、それぞれｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相を表す。なお、式（１３）の導出については後述する。
【００４９】
なお、従来技術では式（１）に基づきδｋを離散全探索によって推定するため多くの計算コストを要していたが、本実施形態では式（１３）に基づきδｆ，ｋを算出するため全探索を要せず計算コストを小さくできる。また、式（１３）に基づき周波数ｆ毎に到達時間差δｆ，ｋを求める点が従来技術とは異なる。
【００５０】
（時間補正部１２２１２）
Ｍステップ計算部１２２の時間差推定部１２２１の時間補正部１２２１２は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、一回前の繰り返しで得られているパラメタθ’（但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、式（１３）にて推定した到達時間差δｆ，ｋが内包する±πの不定性を補正する（ｓ２６）。この補正は、式（１３）の左辺２πｆδｆ，ｋが−πからπの値を取るのに対し、式（１３）の右辺の逆正接が−π／２からπ／２の値しか返すことができないため、２πｆδｆ，ｋが−π〜−π／２及びπ／２〜πの範囲を取る場合の値を正しく求めるために必要である。式（１３）で得られた値を２πｆδ’ｆ，ｋと記載すると、補正は以下のように行なう。
【００５１】
【数１３】
【００５２】
ここでξｆ，ｔ，ｋ、φｆはそれぞれ式（１４）、（１５）で与えられ、ψｓｋ（但し添え字ｓｋはｓｋを表す）、ψξｋ（但し添え字ξｋはξｋを表す）はそれぞれｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相を表す。なお、式（１７）、（１８）の導出については後述する。
【００５３】
さらに、時間補正部１２２１２は、補正した値２πｆδｆ，ｋを２πｆで除算し、到達時間差δｆ，ｋを求め、音源スペクトル推定部１２２２と雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。
【００５４】
（音源スペクトル推定部１２２２）
Ｍステップ計算部１２２の音源スペクトル推定部１２２２は、到達時間差δｆ，ｋと観測信号ｘｆ，ｔを入力とし、これらの値を用いて、音源スペクトルｓｆ，ｔ，ｋを以下の式（１９）により推定し（ｓ２７）、雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。
【００５５】
【数１４】
【００５６】
ここで、ｂｆ，ｋ、Ｖｆはそれぞれ式（３）、（４）で表される。
【００５７】
（雑音パワー推定部１２２３）
Ｍステップ計算部１２２の雑音パワー推定部１２２３は、到達時間差δｆ，ｋと音源スペクトルｓｆ，ｔ，ｋと観測信号ｘｆ，ｔと事後確率ｍｆ，ｔ，ｋとを入力とし、これらの値を用いて、雑音のパワースペクトルσｆ２を以下の式（２０）により推定し（ｓ２８）、Ｑ関数計算部１２１２に出力する。
【００５８】
【数１５】
【００５９】
なお、Ｔは時間フレームの総数である。
【００６０】
（音源存在確率推定部１２２４）
Ｍステップ計算部１２２の音源存在確率推定部１２２４は、事後確率ｍｆ，ｔ，ｋを入力とし、音源存在確率を以下の式（２１）により推定し（ｓ２９）、Ｑ関数計算部１２１２に出力する。
【００６１】
【数１６】
【００６２】
（Ｑ関数計算部１２１２）
Ｅステップ計算部１２１のＱ関数計算部１２１２は、観測信号ｘｆ，ｔと、パラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝と、事後確率ｍｆ，ｔ，ｋと、を入力とし、これらの値を用いてＱ関数を上述の式（１０）により求める（ｓ３０）。
【００６３】
ｓ２０〜ｓ３０までの処理を終えると、パラメタ推定部１２は、（条件１）Ｑ関数の値の変化量（｜Ｑ（θ｜θｎ−１）−Ｑ（θ｜θｎ）｜）が所定の収束判定閾値Δより小さくなるか、または、（条件２）更新回数ｎが所定の最大更新回数Ｎ（例えばＮ＝２０）以上か否かを判定する（ｓ３１）。
【００６４】
パラメタ推定部１２は、（条件１）、（条件２）の何れかを満たしたときは、パラメタ推定部１２はその時点で取得している最新の事後確率ｍｆ，ｔ，ｋと到達時間差δｆ，ｋを分離信号生成部１３に出力する。
【００６５】
パラメタ推定部１２は、（条件１）、（条件２）の何れも満たさないときは、ＥステップとＭステップを繰り返す（ｓ３１、ｓ２１）。なお、図示しない記憶部にパラメタθとＱ関数の値Ｑ（θ｜θｎ）とを記憶しておき、次の繰り返しの際に用いる。
【００６６】
＜分離信号生成部１３＞
分離信号生成部１３は、事後確率ｍｆ，ｔ，ｋと到達時間差δｆ，ｋと観測信号ｘｆ，ｔとを入力とし、以下の式（２２）により、分離信号ｙｆ，ｔ，ｋを生成し（ｓ３）、時間領域変換部１４へ出力する。
【００６７】
【数１７】
【００６８】
なお、音源スペクトルｓｆ，ｔ，ｋは、音源の発する原信号のスペクトルを推定したものであるが、この音源スペクトルを単純に時間領域の信号に変換した場合には、他の音源の発する原信号が残ることがある。上述の式（２２）によって、音源ｋの発する原信号のみを抽出、分離することができる。
【００６９】
＜時間領域変換部１４＞
時間領域変換部１４は、分離信号ｙｆ，ｔ，ｋを入力とし、周波数領域変換部１１において行った周波数領域変換方法に対応する時間領域変換方法（例えば短時間フーリエ逆変換）で、分離信号ｙｆ，ｔ，ｋを時間領域の分離信号ｙｋ（ｔ）に変換し（ｓ４）、音源分離装置１の出力値として出力する。
【００７０】
＜本実施形態のポイント＞
以下、本実施形態のポイントを説明し、式（１３）、（１７）、（１８）の導出方法を説明する。
【００７１】
本実施形態では、
（性質１）到達時間差δｆ，ｋがＲチャネルとＬチャネルの位相差に影響を与える値であること
（性質２）位相差が周期的な値を取る量であること
という２つの性質を利用して、到着時間差δｆ，ｋを推定する。ここで（性質１）は、式（２）にてノイズｎｆ，ｔが十分に小さい場合を考えれば明らかである。（性質２）は、ある２つの位相を表わす量の差Θが、−π≦Θ＜πの範囲の値だけではなく、Θ±２πＭ（Ｍは任意の整数）という周期的な不定性を内包する値を取る性質を持つことを意味する。
【００７２】
（性質２）について、さらに以下にて説明を行なう。式（６）
【００７３】
【数１８】
【００７４】
の右辺における、ｅｘｐのカッコの中のベクトル及び行列ｘ，ｂ，Ｖを、それぞれの成分で表して整理すると、
【００７５】
【数１９】
【００７６】
となる（Ｃはδｆ，ｋに依らない定数）。式（３２）は、位相や角度の分布のように周期的な値を取る変数に対する分布であるＶｏｎ Ｍｉｓｅｓ分布
【００７７】
【数２０】
【００７８】
と同じ形をしていることが分かる（参考文献１参照）。
［参考文献１］Ｃ．Ｍ．ビショップ著，元田ら訳，“パターン認識と機械学習（上）”，シュプリンガー・ジャパン，２００６．
【００７９】
ここで−π＜ｘ≦π、μは分布の平均（−π＜μ≦π）、κ＞０は分布の集中度パラメタ（正規分布での（１／分散）に相当）、Ｉ０（ｘ）は０次の第１種ベッセル関数である。
【００８０】
すなわち、式（３３）の変数ｘが式（３２）のψξｋ−ψｓｋに対応し、式（３３）の平均μが式（３２）の２πｆδｆ，ｋに対応し、式（３３）の集中度κが式（３２）の（｜ξｆ，ｔ，ｋ｜｜ｓｆ，ｔ，ｋ｜）／（σｆ２（１−φｆ２））に対応する。
【００８１】
説明をより直感的にするため、雑音のパラメタがφｆ＝０である場合を考える（これは、雑音が式（４）に示す拡散性雑音ではなく、分散σｆ２のガウス雑音が観測ｘＬ，ｘＲにそれぞれ乗ることを意味する）。このとき、式（１４）によりξｆ，ｔ，ｋ＝ｘｆ，ｔ，Ｒとなるため、式（３２）において、ψξｋはψｘｆ，ｔ，Ｒ（但し、添え字ｘｆ，ｔ，Ｒはｘｆ，ｔ，Ｒを表し、ψｘｆ，ｔ，ＲはマイクＲの観測信号ｘｆ，ｔ，Ｒの位相を表す）となる。また、式（３）、（５）により（但し式（５）においてｎｆ，ｔ＝０）、ｓｆ，ｔ，ｋ＝ｘｆ，ｔ，Ｌとなるため、式（３２）においてψｓｋはψｘｆ，ｔ，Ｌ（但し、添え字ｘｆ，ｔ，Ｌはｘｆ，ｔ，Ｌを表し、ψｘｆ，ｔ，ＬはマイクＬの観測信号ｘｆ，ｔ，Ｌの位相を表す）となる。また、前述の通りξｆ，ｔ，ｋ＝ｘｆ，ｔ，Ｒであり、式（３）、（５）よりｘｆ，ｔ，Ｒ＝ｅｊ２πｆδｋ，ｆ・ｓｆ，ｔ，ｋ（但し、ｅの添え字δｋ，ｆはδｋ，ｆを表す）となるため、式（３２）において、｜ξｆ，ｔ，ｋ｜＝｜ｅｊ２πｆδｋ，ｆ・ｓｆ，ｔ，ｋ｜＝｜ｓｆ，ｔ，ｋ｜となる。よって、式（３２）は
【００８２】
【数２１】
【００８３】
となる。これとＶｏｎ Ｍｉｓｅｓ分布との解釈を合わせると、
（１）ＲチャネルとＬチャネルの位相差ψｘｆ，ｔ，Ｒ−ψｘｆ，ｔ，Ｌという周期的な値を取る変数の分布が、平均μ＝２πｆδｆ，ｋを取る、
（２）Ｖｏｎ Ｍｉｓｅｓ分布の集中度κが、ＳＮ比｜ｓｆ，ｋ，ｔ｜２／σｆ２に対応するようになる。すなわち、ＳＮ比が低い（条件が悪い）と、ＲチャネルとＬチャネルの位相差の値の集中度が下がる（＝分散が大きくなる）。これは、ＳＮ比が低い条件においては位相差の測定値がばらつく現象と対応する、
の２点が言える。
【００８４】
本実施形態では、（性質１）、（性質２）を利用しており、実施の手続きとしては、ＲチャネルとＬチャネルの位相差に関する量が、周期的な値を取る変数に対する分布であるＶｏｎ Ｍｉｓｅｓ分布で表現できることに着目し、Ｖｏｎ Ｍｉｓｅｓ分布のパラメタに対する最尤推定により、分布の平均値μ＝２πｆδｆ，ｋを推定することで、信号到達時間差パラメタδｆ，ｋを推定する。
【００８５】
式（３１）のｐ（ｘｆ，ｔ｜ｋ，θ）をＱ関数の式（１０）に代入し、
【００８６】
【数２２】
【００８７】
を解くことで、式（１３）が得られる。
また時間補正部１２２１２における関数Ｆの式（１７）、（１８）は、Ｑ関数の２階微分
【００８８】
【数２３】
【００８９】
である。時間補正部１２２１２では、式（１３）で得られた解δ’ｆ，ｋがＱ関数の極大値・極小値のどちらを与えるかを、Ｆ（２πｆδ’ｆ，ｋ）の値の正負にて調べ、δ’ｆ，ｋが極小値を与える場合に、上述の＋πまたは−πの補正を行なう。
【００９０】
＜効果＞
従来法においては、時間差推定部において解析的な更新式が与えられていなかったため、多くの計算コストを要する全探索操作が必要であった。よって、時間差推定部の計算コストを削減し、高速な音源分離手段を提供することが課題である。本実施形態では、到着時間差δｆ，ｋが、マイクＲとマイクＬの位相差に影響を与える値であることと、位相差が周期的な値を取る性質を持つことに着目し、到着時間差δｆ，ｋを推定する。これにより、従来必要であった全探索操作が不要となるため、高速な音源分離手段を提供することが可能となる。
【００９１】
＜シミュレーション結果＞
発明の効果を示すため実験を行なった。図６に示す部屋において、音源数は２つ（７０度と１５０度）又は３つ（３０，７０，１５０度）とした。音源は、英語話者音声を用い、音源数２及び３の場合それぞれにおいて１０通りの音源組合せにて実験を行なった。
【００９２】
雑音としては、平均０、共分散行列σｆ２Ｖｆ（Ｖｆは式（４）により与えられる）のガウスノイズをＳＮ比約２５ｄＢにて重畳した。部屋の残響時間は１３０ｍｓ、サンプリング周波数は８ｋＨｚ，短時間フーリエ変換の窓長及びシフト長はそれぞれ６４ｍｓ，１６ｍｓとした。
【００９３】
従来法では到達時間差δｋを推定する際に、音源位置Θ（図６参照）を０度から１８０度まで１度きざみで変化させ、それに対応する１８１種類の到達時間差δｋの値について全探索を行なった。
【００９４】
図７に、ＳＩＲ（信号対雑音比、雑音には他話者音声含む）、ＳＤＲ（信号対歪み比）、及びパソコン（IntelXeon（登録商標）X5650 2.67GHz(6Core)×2CPU）における計算時間を示す。図７Ａは音源が二つの場合（７０度と１５０度）を、図７Ｂは音源が三つの場合（３０度，７０度，１５０度）を示す。ＳＩＲとＳＤＲは大きな値であるほど良い性能であることを示す。数字は、１０通りの音源組合せの平均値である。図７Ａ、Ｂに示す通り、本実施形態は、１／１０程度の計算時間で、従来法とほぼ同程度の分離性能を達成できることが見てとれる。
【００９５】
＜その他の変形例＞
本実施形態のポイントは上述の通り、到達時間差の推定方法である。従って、他の処理やパラメタの推定方法については、上記の実施形態に限定されるものではなく、他の従来技術を用いてもよい。
【００９６】
本実施形態では、各部間で直接データを受け渡しているが、図示しない記憶部を介して、各データを読み書きしてもよい。
【００９７】
また、本実施形態では、雑音パワー推定部１２２３は音源スペクトルｓｆ，ｔ，ｋを音源スペクトル推定部１２２２から取得しているが、到達時間差δｆ，ｋと観測信号ｘｆ，ｔとを用いて式（１９）に基づき雑音パワー推定部１２２３で計算する構成としてもよい。
【００９８】
本実施形態では、拡散性雑音の場合を想定しているが、他の特性を持つ雑音であってもよい。その場合、雑音の特性に応じて、式（４）や式（１５）のｓｉｎｃ（２πｆＤ／ｃ）を適宜変更すればよい。
【００９９】
分離信号生成部１３は、事後確率ｍｆ，ｔ，ｋと音源スペクトルｓｆ，ｔ，ｋとを入力とし、式（２２）に代えて、以下の式により、分離信号ｙｆ，ｔ，ｋを生成してもよい。
【０１００】
【数２４】
【０１０１】
この場合、パラメタ推定部１２は、Ｑ関数の値の変化量（｜Ｑ（θ｜θｎ−１）−Ｑ（θ｜θｎ）｜）が所定の収束判定閾値Δより小さくなったとき、または、更新回数ｎが所定の最大更新回数Ｎ以上になったとき（ｓ３１）に取得している最新の事後確率ｍｆ，ｔ，ｋと音源スペクトルｓｆ，ｔ，ｋを分離信号生成部１３に出力する。このような構成により、式（２２）と同様に分離信号ｙｆ，ｔ，ｋを生成することができる（式（１９）参照）。
【０１０２】
Ｍステップの計算順序は、本実施形態の計算順序に限らない。例えば、時間差推定部１２２１と音源スペクトル推定部１２２２の計算順序はどちらを先に行ってもよい。図８は音源スペクトル推定部１２２２の音源スペクトル推定処理（ｓ２７）を行った後に、時間差推定部１２２１の到達時間差推定処理（ｓ２５、ｓ２６）を行う場合の機能ブロック図を示す。以下、図８を用いて説明する。
【０１０３】
Ｍステップ計算部１２２の音源スペクトル推定部１２２２は、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの到達時間差δｆ，ｋである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の音源スペクトル計算においては、前述の初期化したパラメタ）と観測信号ｘｆ，ｔを入力とし、これらの値を用いて、音源スペクトルｓｆ，ｔ，ｋを式（１９）により推定し、雑音パワー推定部１２２３とＱ関数計算部１２１２と時間差推定部１２２１と時間補正部１２２１２とに出力する。
【０１０４】
Ｍステップ計算部１２２の時間差推定部１２２１の逆正接計算部１２２１１は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、音源スペクトルｓｆ，ｔ，ｋとを入力とし、これらの値を用いて、到達時間差δｆ，ｋを式（１３）により推定し、時間補正部１２２１２に出力する。
【０１０５】
Ｍステップ計算部１２２の時間差推定部１２２１の時間補正部１２２１２は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、音源スペクトルｓｆ，ｔ，ｋと、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの雑音のパワースペクトルσｆ２である。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、式（１３）にて推定した到達時間差δｆ，ｋが内包する±πの不定性を式（１６）に基づき補正する。さらに、時間補正部１２２１２は、補正した値２πｆδｆ，ｋを２πｆで除算し、到達時間差δｆ，ｋを求め、雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。
【０１０６】
上述の構成とすることで、第一実施形態と同様の効果を得ることができる。なお、同様に、一回前の繰り返しで得られているパラメタθ’のうちの到達時間差δｆ，ｋに基づき雑音パワー推定部１２２３の雑音パワー推定処理（ｓ２８）を行った後に、時間差推定部１２２１の到達時間差推定処理（ｓ２５、ｓ２６）を行ってもよい。
【０１０７】
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【０１０８】
＜プログラム及び記録媒体＞
上述した音源分離装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置（各種実施例で図に示した機能構成をもつ装置）として機能させるためのプログラム、またはその処理手順（各実施例で示したもの）の各過程をコンピュータに実行させるためのプログラムを、ＣＤ−ＲＯＭ、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。
【技術分野】
【０００１】
本発明は信号処理の技術分野に属する。特に本発明は複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を推定し、分離抽出する音源分離技術に関する。特に、原信号やそれらがどのように混ざったかの情報を用いずに複数の原信号とノイズとが混在している観測信号のみから、それぞれの原信号を推定する、ブラインド音源分離技術に属する。
【背景技術】
【０００２】
非特許文献１が音源分離の従来技術として知られている。非特許文献１では、音源ｋから発せられる原信号の、二個のマイクへの到達時間差δｋを
【０００３】
【数１】
【０００４】
として推定している。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】和泉洋介，小野順貴，嵯峨山茂樹，“スパースな混合モデルに基づく雑音・残響環境下の劣決定ブラインド音源分離”，電子情報通信学会総合大会講演論文集，２００８年３月
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、従来技術では、到達時間差δｋの解析的な更新式は与えられていないため、多くの計算コストを要する全探索操作によって、到達時間差δｋを推定する必要がある。
【０００７】
本発明は、到達時間差δの性質に着目し、到達時間差δを効率的に推定する方法を与え、従来技術において必要であった全探索操作を不要とし、高速な音源分離技術を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の課題を解決するために、本発明の第一の態様によれば、複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する。原信号の音源のインデックスをｋとし、周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、二個のマイクへの原信号の到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と原信号のスペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定する。観測信号ｘｆ，ｔと推定されたパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝とから分離信号ｙｆ，ｔ，ｋを生成する。観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、二個のマイクの間隔をＤとし、原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、スペクトルｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、到達時間差δｆ，ｋを
【０００９】
【数２】
【００１０】
として推定する。推定された到達時間差δｆ，ｋが内包する±πの不定性を補正する。
【発明の効果】
【００１１】
本発明は、到達時間差δを効率的に推定する方法を与え、高速な音源分離ができるという効果を奏する。
【図面の簡単な説明】
【００１２】
【図１】第一実施形態に係る音源分離装置の機能ブロック図。
【図２】第一実施形態に係る音源分離装置の処理フローを示す図。
【図３】Ｅステップ計算部の機能ブロック図。
【図４】Ｍステップ計算部の機能ブロック図。
【図５】パラメタ推定部の処理フローを示す図。
【図６】シミュレーションを行った環境を示す図。
【図７】図７Ａは音源数が２つの場合のシミュレーション結果を、図７Ｂは音源数が３つの場合のシミュレーション結果を示す図。
【図８】音源スペクトル推定部の処理を時間差推定部の処理の前に行う場合のＭステップ計算部の機能ブロック図。
【発明を実施するための形態】
【００１３】
以下、本発明の実施形態について説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
【００１４】
＜第一実施形態に係る音源分離装置１＞
図１は本実施形態に係る音源分離装置１の機能ブロック図を、図２はその処理フローを示す。
【００１５】
音源分離装置１は周波数領域変換部１１とパラメタ推定部１２と分離信号生成部１３と時間領域変換部１４とを含み、さらにパラメタ推定部１２はＥステップ計算部１２１とＭステップ計算部１２２とを含む。
【００１６】
まず、観測信号について説明する。信号原（以下「音源」とも言う）がＫ個あり、ｋ（ｋ＝１，２，…，Ｋ）を音源のインデックスとし、音源ｋから発せられる信号（以下「原信号」という）をｓｋ（ｔ）とする。複数の原信号ｓ１（ｔ），・・・，ｓＫ（ｔ）がノイズとともに二個のマイクＬ，Ｒで観測される状況で、マイクＬで観測される時間領域の観測信号をｘＬ（ｔ）とし、マイクＲで観測される時間領域の観測信号をｘＲ（ｔ）とし、二個のマイクＬ，Ｒで観測される時間領域の観測信号をｘ（ｔ）＝［ｘＬ（ｔ），ｘＲ（ｔ）］Ｔとする。「Ｔ」は転置を表す。ｔはフレーム番号及びそのフレーム番号に対応する時刻を表す。Ｔを時間フレームの総数とすると、ｔ＝０，１，…，Ｔ−１である。ここで周波数領域の観測信号をｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔと表記する。なお、ｆはサンプリング周波数ｆｓをＦ等分した離散点であり、ｆ∈｛０，ｆｓ／Ｆ，…，（Ｆ−１）ｆｓ／Ｆ｝である。以降、断りのない場合、観測信号とは、周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔを指し、時間領域の観測信号の場合はそれを明記する。
ここで、観測信号は、
【００１７】
【数３】
【００１８】
と表されると仮定する。ここで、ｓｆ，ｔ，ｋは原信号ｓｋ（ｔ）のスペクトル（以下「音源スペクトル」ともいう）を、ｎｆ，ｋ＝［ｎｆ，ｋ，Ｌ，ｎｆ，ｋ，Ｒ］はマイクＬ，Ｒにおける加算的雑音を表す。またｂｆ，ｋ＝［ｂｆ，ｋ，Ｌ，ｂｆ，ｋ，Ｒ］Ｔは音源ｋに関するステアリングベクトル（音源ｋの方向を特定するベクトルであり、方向ベクトルともいう）であり、原信号ｓｋ（ｔ）のマイクＬ、Ｒへの到達時間差をδｋとすると、
【００１９】
【数４】
【００２０】
である（非特許文献１参照）。本実施形態では、原信号の観測時間内においては、音源及びマイクは固定されており、またＫ個の音源は全て、異なる位置に配置されているとする。すなわち、ステアリングベクトルｂｆ，ｋは時間ｔに依らず、ｋの値によって異なる値を取るものと仮定する。音源分離の目的は、観測信号ｘｆ，ｔのみを用いて、全ての音源スペクトルｓｆ，ｔ，ｋを推定することである。
【００２１】
本実施形態に係る音源分離装置１は、時間領域の観測信号ｘ（ｔ）＝［ｘＬ（ｔ），ｘＲ（ｔ）］Ｔを入力とし、時間領域の分離信号（推定された各原信号）ｙｋ（ｔ）を出力する。以下、各部の処理内容を説明する。
【００２２】
＜周波数領域変換部１１＞
まず、周波数領域変換部１１は、マイクＬ、Ｒで収音した時間領域の観測信号ｘ（ｔ）＝［ｘＬ（ｔ），ｘＲ（ｔ）］Ｔを入力とし、これを短時間フーリエ変換等により周波数領域の観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔに変換し（ｓ１）、パラメタ推定部１２及び分離信号生成部１３に出力する。
【００２３】
＜パラメタ推定部１２＞
次に、パラメタ推定部１２は、観測信号ｘｆ，ｔから、音源分離のために必要なパラメタθを推定する（ｓ２）。なおθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝であり、δｆ，ｋは原信号の、二個のマイクへの到達時間差を表す。σｆ２は雑音のパワースペクトルを、ｓｆ，ｔ，ｋは音源スペクトルを、ｐ（ｋ｜θ）は音源存在確率（混合信号中の音源ｋの寄与率）を表す。
【００２４】
本実施形態では、上記パラメタを推定するために、以下の２つの仮定を用いる。
【００２５】
（仮定１）は、雑音ｎが、平均０、共分散行列σｆ２Ｖｆの正規分布に従う定常雑音でモデル化できるという仮定である。ここでσｆ２は周波数ｆにおける雑音のパワーであり、Ｖｆは例えば拡散性雑音の場合
【００２６】
【数５】
【００２７】
で与えられる（非特許文献１参照）。ここでｃは原信号の速度（音速等）、Ｄは二個のマイクの間隔である。
【００２８】
（仮定２）は、原信号がスパースな信号（つまり、成分のうち０でないもの（非零の成分）がまばらである信号、言い換えると、ほとんどの（または多くの）成分が０である信号）であるという仮定である。すなわち、ある時間周波数（ｆ，ｔ）において、たかだか１つの原信号のみが支配的であると仮定する。（仮定２）によると、式（２）は
【００２９】
【数６】
【００３０】
と表記できる（非特許文献１参照）。
上記（仮定１）、（仮定２）に基づけば、到達時間差δｆ，ｋとなる方向に存在する音源ｋから発せられる原信号が到来してｘｆ，ｔが観測される尤度は、
【００３１】
【数７】
【００３２】
で与えられる（非特許文献１参照）。「Ｈ」はエルミート転置を表す。ここでｂｆ，ｋは式（３）で表される。これより、到達時間差δｆ，ｋ、音源スペクトルｓｆ，ｔ，ｋ及び雑音のパワースペクトルσｆ２は、対数尤度関数
【００３３】
【数８】
【００３４】
を最大化するパラメタとして最尤推定により求める（非特許文献１参照）。但し、上記において、ｐ（ｋ｜δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ）は音源存在確率を表し、混合信号中の音源ｋの寄与率であり、
【００３５】
【数９】
【００３６】
である（非特許文献１参照）。
具体的には、期待値最大化法（以下「ＥＭアルゴリズム」ともいう）を適用し、Ｑ関数
【００３７】
【数１０】
【００３８】
を最大とするパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝を、以下のＥステップ及びＭステップの繰り返しにより求める（非特許文献１参照）。ここでｍｆ，ｔ，ｋは後述する式（１１）で、ｐ（ｘｆ，ｔ｜ｋ，θ’）は前述の式（６）で与えられ、θ’は、１回前の繰り返しで得られているパラメタを意味する。
【００３９】
この繰り返し計算を、パラメタ推定部１２にて行なう。以下、パラメタ推定部１２の処理の詳細を説明する。
【００４０】
図３はＥステップ計算部１２１の機能ブロック図を、図４はＭステップ計算部１２２の機能ブロック図を、図５はパラメタ推定部１２の処理フローを示す。
【００４１】
まず、パラメタ推定部１２は、各パラメタを初期化する（ｓ２０）。パラメタのうちσｆ２、δｆ，ｋ、及びｐ（ｋ｜θ）を初期化する。例えば、σｆ２＝｜ｘｆ，ｔ＝０，Ｌ｜２，δｆ，ｋ＝（Ｄ／ｃ）ｃｏｓαｋ（ｃは原信号の速度（音速等）、Ｄはマイク間隔，αｋは音源ｋの方向の初期値（−π／２〜π／２の間の適当な値）），ｐ（ｋ｜θ）＝１／Ｋとする。さらに、初期化したパラメタδｆ，ｋとｘｆ，ｔ＝０とを用いて、式（３）及び後述する式（１９）に基づき、ｓｆ，ｔ，ｋを初期化する。更新回数ｎ＝０とする。なお、最大更新回数Ｎ及び収束判定閾値Δは、当該装置の設計者や利用者等により予め設定されているものとする。
【００４２】
Ｅステップ計算部１２１は事後確率推定部１２１１とＱ関数計算部１２１２とを含む（図３参照）。Ｍステップ計算部１２２は時間差推定部１２２１と音源スペクトル推定部１２２２と雑音パワー推定部１２２３と音源存在確率推定部１２２４とを含み、時間差推定部１２２１はさらに逆正接計算部１２２１１と時間補正部１２２１２とを備える（図４参照）。事後確率推定部１２１１とＱ関数計算部１２１２とにおける処理を併せてＥステップと呼び、時間差推定部１２２１と音源スペクトル推定部１２２２と雑音パワー推定部１２２３と音源存在確率推定部１２２４とにおける処理を併せてＭステップと呼ぶ。
【００４３】
（事後確率推定部１２１１）
Ｅステップ計算部１２１の事後確率推定部１２１１は、観測信号ｘｆ，ｔと、一回前の繰り返しで得られているパラメタθ’（但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の事後確率推定においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて事後確率ｍｆ，ｔ，ｋ＝ｐ（ｋ｜ｘｆ，ｔ，θ’）を以下の式（１１）により求め（ｓ２２）、Ｑ関数計算部１２１２とＭステップ計算部１２２とに出力する。
【００４４】
【数１１】
【００４５】
なお、ｐ（ｘｆ，ｔ｜ｋ，θ’）は式（６）により与えられる。なお、事後確率ｍｆ，ｔ，ｋは観測信号ｘｆ，ｔが音源ｋに帰属する事後確率を表す。
【００４６】
（逆正接計算部１２２１１）
Ｍステップ計算部１２２の時間差推定部１２２１の逆正接計算部１２２１１は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの音源スペクトルｓｆ，ｔ，ｋである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の逆正接計算においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、到達時間差δｆ，ｋ（より詳しく言うと到達時間差δｆ，ｋに２πｆを乗じた値２πｆδｆ，ｋ）を以下の式（１３）により推定し（ｓ２５）、時間補正部１２２１２に出力する。
【００４７】
【数１２】
【００４８】
ψｓｋ（但し添え字ｓｋはｓｋを表す）及びψξｋ（但し添え字ξｋはξｋを表す）は、それぞれｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相を表す。なお、式（１３）の導出については後述する。
【００４９】
なお、従来技術では式（１）に基づきδｋを離散全探索によって推定するため多くの計算コストを要していたが、本実施形態では式（１３）に基づきδｆ，ｋを算出するため全探索を要せず計算コストを小さくできる。また、式（１３）に基づき周波数ｆ毎に到達時間差δｆ，ｋを求める点が従来技術とは異なる。
【００５０】
（時間補正部１２２１２）
Ｍステップ計算部１２２の時間差推定部１２２１の時間補正部１２２１２は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、一回前の繰り返しで得られているパラメタθ’（但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、式（１３）にて推定した到達時間差δｆ，ｋが内包する±πの不定性を補正する（ｓ２６）。この補正は、式（１３）の左辺２πｆδｆ，ｋが−πからπの値を取るのに対し、式（１３）の右辺の逆正接が−π／２からπ／２の値しか返すことができないため、２πｆδｆ，ｋが−π〜−π／２及びπ／２〜πの範囲を取る場合の値を正しく求めるために必要である。式（１３）で得られた値を２πｆδ’ｆ，ｋと記載すると、補正は以下のように行なう。
【００５１】
【数１３】
【００５２】
ここでξｆ，ｔ，ｋ、φｆはそれぞれ式（１４）、（１５）で与えられ、ψｓｋ（但し添え字ｓｋはｓｋを表す）、ψξｋ（但し添え字ξｋはξｋを表す）はそれぞれｓｆ，ｔ，ｋ及びξｆ，ｔ，ｋの位相を表す。なお、式（１７）、（１８）の導出については後述する。
【００５３】
さらに、時間補正部１２２１２は、補正した値２πｆδｆ，ｋを２πｆで除算し、到達時間差δｆ，ｋを求め、音源スペクトル推定部１２２２と雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。
【００５４】
（音源スペクトル推定部１２２２）
Ｍステップ計算部１２２の音源スペクトル推定部１２２２は、到達時間差δｆ，ｋと観測信号ｘｆ，ｔを入力とし、これらの値を用いて、音源スペクトルｓｆ，ｔ，ｋを以下の式（１９）により推定し（ｓ２７）、雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。
【００５５】
【数１４】
【００５６】
ここで、ｂｆ，ｋ、Ｖｆはそれぞれ式（３）、（４）で表される。
【００５７】
（雑音パワー推定部１２２３）
Ｍステップ計算部１２２の雑音パワー推定部１２２３は、到達時間差δｆ，ｋと音源スペクトルｓｆ，ｔ，ｋと観測信号ｘｆ，ｔと事後確率ｍｆ，ｔ，ｋとを入力とし、これらの値を用いて、雑音のパワースペクトルσｆ２を以下の式（２０）により推定し（ｓ２８）、Ｑ関数計算部１２１２に出力する。
【００５８】
【数１５】
【００５９】
なお、Ｔは時間フレームの総数である。
【００６０】
（音源存在確率推定部１２２４）
Ｍステップ計算部１２２の音源存在確率推定部１２２４は、事後確率ｍｆ，ｔ，ｋを入力とし、音源存在確率を以下の式（２１）により推定し（ｓ２９）、Ｑ関数計算部１２１２に出力する。
【００６１】
【数１６】
【００６２】
（Ｑ関数計算部１２１２）
Ｅステップ計算部１２１のＱ関数計算部１２１２は、観測信号ｘｆ，ｔと、パラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝と、事後確率ｍｆ，ｔ，ｋと、を入力とし、これらの値を用いてＱ関数を上述の式（１０）により求める（ｓ３０）。
【００６３】
ｓ２０〜ｓ３０までの処理を終えると、パラメタ推定部１２は、（条件１）Ｑ関数の値の変化量（｜Ｑ（θ｜θｎ−１）−Ｑ（θ｜θｎ）｜）が所定の収束判定閾値Δより小さくなるか、または、（条件２）更新回数ｎが所定の最大更新回数Ｎ（例えばＮ＝２０）以上か否かを判定する（ｓ３１）。
【００６４】
パラメタ推定部１２は、（条件１）、（条件２）の何れかを満たしたときは、パラメタ推定部１２はその時点で取得している最新の事後確率ｍｆ，ｔ，ｋと到達時間差δｆ，ｋを分離信号生成部１３に出力する。
【００６５】
パラメタ推定部１２は、（条件１）、（条件２）の何れも満たさないときは、ＥステップとＭステップを繰り返す（ｓ３１、ｓ２１）。なお、図示しない記憶部にパラメタθとＱ関数の値Ｑ（θ｜θｎ）とを記憶しておき、次の繰り返しの際に用いる。
【００６６】
＜分離信号生成部１３＞
分離信号生成部１３は、事後確率ｍｆ，ｔ，ｋと到達時間差δｆ，ｋと観測信号ｘｆ，ｔとを入力とし、以下の式（２２）により、分離信号ｙｆ，ｔ，ｋを生成し（ｓ３）、時間領域変換部１４へ出力する。
【００６７】
【数１７】
【００６８】
なお、音源スペクトルｓｆ，ｔ，ｋは、音源の発する原信号のスペクトルを推定したものであるが、この音源スペクトルを単純に時間領域の信号に変換した場合には、他の音源の発する原信号が残ることがある。上述の式（２２）によって、音源ｋの発する原信号のみを抽出、分離することができる。
【００６９】
＜時間領域変換部１４＞
時間領域変換部１４は、分離信号ｙｆ，ｔ，ｋを入力とし、周波数領域変換部１１において行った周波数領域変換方法に対応する時間領域変換方法（例えば短時間フーリエ逆変換）で、分離信号ｙｆ，ｔ，ｋを時間領域の分離信号ｙｋ（ｔ）に変換し（ｓ４）、音源分離装置１の出力値として出力する。
【００７０】
＜本実施形態のポイント＞
以下、本実施形態のポイントを説明し、式（１３）、（１７）、（１８）の導出方法を説明する。
【００７１】
本実施形態では、
（性質１）到達時間差δｆ，ｋがＲチャネルとＬチャネルの位相差に影響を与える値であること
（性質２）位相差が周期的な値を取る量であること
という２つの性質を利用して、到着時間差δｆ，ｋを推定する。ここで（性質１）は、式（２）にてノイズｎｆ，ｔが十分に小さい場合を考えれば明らかである。（性質２）は、ある２つの位相を表わす量の差Θが、−π≦Θ＜πの範囲の値だけではなく、Θ±２πＭ（Ｍは任意の整数）という周期的な不定性を内包する値を取る性質を持つことを意味する。
【００７２】
（性質２）について、さらに以下にて説明を行なう。式（６）
【００７３】
【数１８】
【００７４】
の右辺における、ｅｘｐのカッコの中のベクトル及び行列ｘ，ｂ，Ｖを、それぞれの成分で表して整理すると、
【００７５】
【数１９】
【００７６】
となる（Ｃはδｆ，ｋに依らない定数）。式（３２）は、位相や角度の分布のように周期的な値を取る変数に対する分布であるＶｏｎ Ｍｉｓｅｓ分布
【００７７】
【数２０】
【００７８】
と同じ形をしていることが分かる（参考文献１参照）。
［参考文献１］Ｃ．Ｍ．ビショップ著，元田ら訳，“パターン認識と機械学習（上）”，シュプリンガー・ジャパン，２００６．
【００７９】
ここで−π＜ｘ≦π、μは分布の平均（−π＜μ≦π）、κ＞０は分布の集中度パラメタ（正規分布での（１／分散）に相当）、Ｉ０（ｘ）は０次の第１種ベッセル関数である。
【００８０】
すなわち、式（３３）の変数ｘが式（３２）のψξｋ−ψｓｋに対応し、式（３３）の平均μが式（３２）の２πｆδｆ，ｋに対応し、式（３３）の集中度κが式（３２）の（｜ξｆ，ｔ，ｋ｜｜ｓｆ，ｔ，ｋ｜）／（σｆ２（１−φｆ２））に対応する。
【００８１】
説明をより直感的にするため、雑音のパラメタがφｆ＝０である場合を考える（これは、雑音が式（４）に示す拡散性雑音ではなく、分散σｆ２のガウス雑音が観測ｘＬ，ｘＲにそれぞれ乗ることを意味する）。このとき、式（１４）によりξｆ，ｔ，ｋ＝ｘｆ，ｔ，Ｒとなるため、式（３２）において、ψξｋはψｘｆ，ｔ，Ｒ（但し、添え字ｘｆ，ｔ，Ｒはｘｆ，ｔ，Ｒを表し、ψｘｆ，ｔ，ＲはマイクＲの観測信号ｘｆ，ｔ，Ｒの位相を表す）となる。また、式（３）、（５）により（但し式（５）においてｎｆ，ｔ＝０）、ｓｆ，ｔ，ｋ＝ｘｆ，ｔ，Ｌとなるため、式（３２）においてψｓｋはψｘｆ，ｔ，Ｌ（但し、添え字ｘｆ，ｔ，Ｌはｘｆ，ｔ，Ｌを表し、ψｘｆ，ｔ，ＬはマイクＬの観測信号ｘｆ，ｔ，Ｌの位相を表す）となる。また、前述の通りξｆ，ｔ，ｋ＝ｘｆ，ｔ，Ｒであり、式（３）、（５）よりｘｆ，ｔ，Ｒ＝ｅｊ２πｆδｋ，ｆ・ｓｆ，ｔ，ｋ（但し、ｅの添え字δｋ，ｆはδｋ，ｆを表す）となるため、式（３２）において、｜ξｆ，ｔ，ｋ｜＝｜ｅｊ２πｆδｋ，ｆ・ｓｆ，ｔ，ｋ｜＝｜ｓｆ，ｔ，ｋ｜となる。よって、式（３２）は
【００８２】
【数２１】
【００８３】
となる。これとＶｏｎ Ｍｉｓｅｓ分布との解釈を合わせると、
（１）ＲチャネルとＬチャネルの位相差ψｘｆ，ｔ，Ｒ−ψｘｆ，ｔ，Ｌという周期的な値を取る変数の分布が、平均μ＝２πｆδｆ，ｋを取る、
（２）Ｖｏｎ Ｍｉｓｅｓ分布の集中度κが、ＳＮ比｜ｓｆ，ｋ，ｔ｜２／σｆ２に対応するようになる。すなわち、ＳＮ比が低い（条件が悪い）と、ＲチャネルとＬチャネルの位相差の値の集中度が下がる（＝分散が大きくなる）。これは、ＳＮ比が低い条件においては位相差の測定値がばらつく現象と対応する、
の２点が言える。
【００８４】
本実施形態では、（性質１）、（性質２）を利用しており、実施の手続きとしては、ＲチャネルとＬチャネルの位相差に関する量が、周期的な値を取る変数に対する分布であるＶｏｎ Ｍｉｓｅｓ分布で表現できることに着目し、Ｖｏｎ Ｍｉｓｅｓ分布のパラメタに対する最尤推定により、分布の平均値μ＝２πｆδｆ，ｋを推定することで、信号到達時間差パラメタδｆ，ｋを推定する。
【００８５】
式（３１）のｐ（ｘｆ，ｔ｜ｋ，θ）をＱ関数の式（１０）に代入し、
【００８６】
【数２２】
【００８７】
を解くことで、式（１３）が得られる。
また時間補正部１２２１２における関数Ｆの式（１７）、（１８）は、Ｑ関数の２階微分
【００８８】
【数２３】
【００８９】
である。時間補正部１２２１２では、式（１３）で得られた解δ’ｆ，ｋがＱ関数の極大値・極小値のどちらを与えるかを、Ｆ（２πｆδ’ｆ，ｋ）の値の正負にて調べ、δ’ｆ，ｋが極小値を与える場合に、上述の＋πまたは−πの補正を行なう。
【００９０】
＜効果＞
従来法においては、時間差推定部において解析的な更新式が与えられていなかったため、多くの計算コストを要する全探索操作が必要であった。よって、時間差推定部の計算コストを削減し、高速な音源分離手段を提供することが課題である。本実施形態では、到着時間差δｆ，ｋが、マイクＲとマイクＬの位相差に影響を与える値であることと、位相差が周期的な値を取る性質を持つことに着目し、到着時間差δｆ，ｋを推定する。これにより、従来必要であった全探索操作が不要となるため、高速な音源分離手段を提供することが可能となる。
【００９１】
＜シミュレーション結果＞
発明の効果を示すため実験を行なった。図６に示す部屋において、音源数は２つ（７０度と１５０度）又は３つ（３０，７０，１５０度）とした。音源は、英語話者音声を用い、音源数２及び３の場合それぞれにおいて１０通りの音源組合せにて実験を行なった。
【００９２】
雑音としては、平均０、共分散行列σｆ２Ｖｆ（Ｖｆは式（４）により与えられる）のガウスノイズをＳＮ比約２５ｄＢにて重畳した。部屋の残響時間は１３０ｍｓ、サンプリング周波数は８ｋＨｚ，短時間フーリエ変換の窓長及びシフト長はそれぞれ６４ｍｓ，１６ｍｓとした。
【００９３】
従来法では到達時間差δｋを推定する際に、音源位置Θ（図６参照）を０度から１８０度まで１度きざみで変化させ、それに対応する１８１種類の到達時間差δｋの値について全探索を行なった。
【００９４】
図７に、ＳＩＲ（信号対雑音比、雑音には他話者音声含む）、ＳＤＲ（信号対歪み比）、及びパソコン（IntelXeon（登録商標）X5650 2.67GHz(6Core)×2CPU）における計算時間を示す。図７Ａは音源が二つの場合（７０度と１５０度）を、図７Ｂは音源が三つの場合（３０度，７０度，１５０度）を示す。ＳＩＲとＳＤＲは大きな値であるほど良い性能であることを示す。数字は、１０通りの音源組合せの平均値である。図７Ａ、Ｂに示す通り、本実施形態は、１／１０程度の計算時間で、従来法とほぼ同程度の分離性能を達成できることが見てとれる。
【００９５】
＜その他の変形例＞
本実施形態のポイントは上述の通り、到達時間差の推定方法である。従って、他の処理やパラメタの推定方法については、上記の実施形態に限定されるものではなく、他の従来技術を用いてもよい。
【００９６】
本実施形態では、各部間で直接データを受け渡しているが、図示しない記憶部を介して、各データを読み書きしてもよい。
【００９７】
また、本実施形態では、雑音パワー推定部１２２３は音源スペクトルｓｆ，ｔ，ｋを音源スペクトル推定部１２２２から取得しているが、到達時間差δｆ，ｋと観測信号ｘｆ，ｔとを用いて式（１９）に基づき雑音パワー推定部１２２３で計算する構成としてもよい。
【００９８】
本実施形態では、拡散性雑音の場合を想定しているが、他の特性を持つ雑音であってもよい。その場合、雑音の特性に応じて、式（４）や式（１５）のｓｉｎｃ（２πｆＤ／ｃ）を適宜変更すればよい。
【００９９】
分離信号生成部１３は、事後確率ｍｆ，ｔ，ｋと音源スペクトルｓｆ，ｔ，ｋとを入力とし、式（２２）に代えて、以下の式により、分離信号ｙｆ，ｔ，ｋを生成してもよい。
【０１００】
【数２４】
【０１０１】
この場合、パラメタ推定部１２は、Ｑ関数の値の変化量（｜Ｑ（θ｜θｎ−１）−Ｑ（θ｜θｎ）｜）が所定の収束判定閾値Δより小さくなったとき、または、更新回数ｎが所定の最大更新回数Ｎ以上になったとき（ｓ３１）に取得している最新の事後確率ｍｆ，ｔ，ｋと音源スペクトルｓｆ，ｔ，ｋを分離信号生成部１３に出力する。このような構成により、式（２２）と同様に分離信号ｙｆ，ｔ，ｋを生成することができる（式（１９）参照）。
【０１０２】
Ｍステップの計算順序は、本実施形態の計算順序に限らない。例えば、時間差推定部１２２１と音源スペクトル推定部１２２２の計算順序はどちらを先に行ってもよい。図８は音源スペクトル推定部１２２２の音源スペクトル推定処理（ｓ２７）を行った後に、時間差推定部１２２１の到達時間差推定処理（ｓ２５、ｓ２６）を行う場合の機能ブロック図を示す。以下、図８を用いて説明する。
【０１０３】
Ｍステップ計算部１２２の音源スペクトル推定部１２２２は、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの到達時間差δｆ，ｋである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の音源スペクトル計算においては、前述の初期化したパラメタ）と観測信号ｘｆ，ｔを入力とし、これらの値を用いて、音源スペクトルｓｆ，ｔ，ｋを式（１９）により推定し、雑音パワー推定部１２２３とＱ関数計算部１２１２と時間差推定部１２２１と時間補正部１２２１２とに出力する。
【０１０４】
Ｍステップ計算部１２２の時間差推定部１２２１の逆正接計算部１２２１１は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、音源スペクトルｓｆ，ｔ，ｋとを入力とし、これらの値を用いて、到達時間差δｆ，ｋを式（１３）により推定し、時間補正部１２２１２に出力する。
【０１０５】
Ｍステップ計算部１２２の時間差推定部１２２１の時間補正部１２２１２は、観測信号ｘｆ，ｔと、事後確率ｍｆ，ｔ，ｋと、音源スペクトルｓｆ，ｔ，ｋと、一回前の繰り返しで得られているパラメタθ’（より詳しく説明するとθ’のうちの雑音のパワースペクトルσｆ２である。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタ）とを入力とし、これらの値を用いて、式（１３）にて推定した到達時間差δｆ，ｋが内包する±πの不定性を式（１６）に基づき補正する。さらに、時間補正部１２２１２は、補正した値２πｆδｆ，ｋを２πｆで除算し、到達時間差δｆ，ｋを求め、雑音パワー推定部１２２３とＱ関数計算部１２１２とに出力する。
【０１０６】
上述の構成とすることで、第一実施形態と同様の効果を得ることができる。なお、同様に、一回前の繰り返しで得られているパラメタθ’のうちの到達時間差δｆ，ｋに基づき雑音パワー推定部１２２３の雑音パワー推定処理（ｓ２８）を行った後に、時間差推定部１２２１の到達時間差推定処理（ｓ２５、ｓ２６）を行ってもよい。
【０１０７】
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【０１０８】
＜プログラム及び記録媒体＞
上述した音源分離装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置（各種実施例で図に示した機能構成をもつ装置）として機能させるためのプログラム、またはその処理手順（各実施例で示したもの）の各過程をコンピュータに実行させるためのプログラムを、ＣＤ−ＲＯＭ、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。
【特許請求の範囲】
【請求項１】
複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離装置であって、
前記原信号の音源のインデックスをｋとし、周波数領域の前記観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、前記二個のマイクへの前記原信号の到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と原信号のスペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定するパラメタ推定手段と、
前記観測信号ｘｆ，ｔと推定されたパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝とから分離信号ｙｆ，ｔ，ｋを生成する分離信号生成手段とを含み、
前記パラメタ推定手段は、
観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、前記二個のマイクの間隔をＤとし、前記原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、前記スペクトルｓｆ，ｔ，ｋ及び前記ξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、前記到達時間差δｆ，ｋを
【数２５】
として推定する逆正接計算手段と、
推定された前記到達時間差δｆ，ｋが内包する±πの不定性を補正する補正手段とを有する、
音源分離装置。
【請求項２】
複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離方法であって、
前記原信号の音源のインデックスをｋとし、周波数領域の前記観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、前記二個のマイクへの前記原信号の到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と原信号のスペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定するパラメタ推定ステップと、
前記観測信号ｘｆ，ｔと推定されたパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝とから分離信号ｙｆ，ｔ，ｋを生成する分離信号生成ステップとを含み、
前記パラメタ推定ステップは、
観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、前記二個のマイクの間隔をＤとし、前記原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、前記スペクトルｓｆ，ｔ，ｋ及び前記ξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、前記到達時間差δｆ，ｋを
【数２６】
として推定する逆正接計算ステップと、
推定された前記到達時間差δｆ，ｋが内包する±πの不定性を補正する補正ステップとを有する、
音源分離方法。
【請求項３】
コンピュータを請求項１記載の音源分離装置として機能させるためのプログラム。
【請求項１】
複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離装置であって、
前記原信号の音源のインデックスをｋとし、周波数領域の前記観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、前記二個のマイクへの前記原信号の到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と原信号のスペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定するパラメタ推定手段と、
前記観測信号ｘｆ，ｔと推定されたパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝とから分離信号ｙｆ，ｔ，ｋを生成する分離信号生成手段とを含み、
前記パラメタ推定手段は、
観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、前記二個のマイクの間隔をＤとし、前記原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、前記スペクトルｓｆ，ｔ，ｋ及び前記ξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、前記到達時間差δｆ，ｋを
【数２５】
として推定する逆正接計算手段と、
推定された前記到達時間差δｆ，ｋが内包する±πの不定性を補正する補正手段とを有する、
音源分離装置。
【請求項２】
複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離方法であって、
前記原信号の音源のインデックスをｋとし、周波数領域の前記観測信号ｘｆ，ｔ＝［ｘｆ，ｔ，Ｌ，ｘｆ，ｔ，Ｒ］Ｔから、前記二個のマイクへの前記原信号の到達時間差δｆ，ｋと雑音のパワースペクトルσｆ２と原信号のスペクトルｓｆ，ｔ，ｋと音源存在確率ｐ（ｋ｜θ）とを推定するパラメタ推定ステップと、
前記観測信号ｘｆ，ｔと推定されたパラメタθ＝｛δｆ，ｋ，σｆ２，ｓｆ，ｔ，ｋ，ｐ（ｋ｜θ）｝とから分離信号ｙｆ，ｔ，ｋを生成する分離信号生成ステップとを含み、
前記パラメタ推定ステップは、
観測信号ｘｆ，ｔが音源ｋに帰属する期待値を示す事後確率をｍｆ，ｔ，ｋとし、前記二個のマイクの間隔をＤとし、前記原信号の速度をｃとし、φｆ＝ｓｉｎｃ（２πｆＤ／ｃ）とし、ξｆ，ｔ，ｋ＝［ｘｆ，ｔ，Ｒ−φｆ（ｘｔ，ｔ，Ｌ−ｓｆ，ｔ，ｋ）］とし、前記スペクトルｓｆ，ｔ，ｋ及び前記ξｆ，ｔ，ｋの位相をそれぞれψｓｋ及びψξｋとし、前記到達時間差δｆ，ｋを
【数２６】
として推定する逆正接計算ステップと、
推定された前記到達時間差δｆ，ｋが内包する±πの不定性を補正する補正ステップとを有する、
音源分離方法。
【請求項３】
コンピュータを請求項１記載の音源分離装置として機能させるためのプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【公開番号】特開２０１３−９７１７６（Ｐ２０１３−９７１７６Ａ）
【公開日】平成２５年５月２０日（２０１３．５．２０）
【国際特許分類】
【出願番号】特願２０１１−２４００５４（Ｐ２０１１−２４００５４）
【出願日】平成２３年１１月１日（２０１１．１１．１）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
【公開日】平成２５年５月２０日（２０１３．５．２０）
【国際特許分類】
【出願日】平成２３年１１月１日（２０１１．１１．１）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
[ Back to top ]