JP2014048398A - 音響信号解析装置、方法、及びプログラム - Google Patents

Abstract

【課題】複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができるようにする。
【解決手段】確率更新部２４１、時変ゲイン更新部２４２、及び分離行列更新部２４３により、各周波数に対する各時間フレームにおいて調波構造を有するパワースペクトルテンプレートの各々が音源毎に選択される確率を要素にもつ三次元配列Π＾を用いて、観測時間周波数成分の三次元配列Ｘが与えられたときの三次元配列Π＾、各音源のパワーを要素にもつ二次元配列Σ＾、及び各周波数ωの分離行列Ｗ_ω＾の事後確率を表す目的関数を最大化するように、三次元配列Π＾、二次元配列Σ＾、及び各周波数ωの分離行列Ｗ_ω＾を更新することを繰り返す。
【選択図】図１

Description

本発明は、音響信号解析装置、方法、及びプログラムに係り、特に、複数のマイクロホンから出力される音響信号の時系列データから、各音源の信号に分離する音響信号解析装置、方法、及びプログラムに関する。

音源の成分と音源からマイクロホンまでの伝達特性がともに未知のもとで、マイクロホン入力信号から個々の音源成分を分離抽出する技術をブラインド音源分離(Blind Source Separation;BSS)という。BSSでは観測信号だけから音源信号とその混合過程を推定する必要があるため、通常は音源に関して何らかの仮定を置き、これにより立てられる規準をもとに両未知変数を推定する最適化問題として定式化される。例えば、観測信号数が音源数以上の場合には、音源信号成分が優ガウス分布に従うという仮定のもとで分離フィルタを最尤推定する独立成分分析(Independent Component Analysis;ICA) と呼ぶ方法が有名である（非特許文献１）。周波数領域におけるICAは、音響信号のような時間遅れを有する混合系に対するBSSの有効なアプローチの1つであるが、周波数帯域毎に信号分離が行われるため、帯域毎の分離結果を同一音源毎にまとめる、いわゆるパーミュテーション問題を解く必要があった。

近年提案されている、ベクトル型変数を用いた独立成分分析(以下、独立ベクトル分析(Independent Vector Analysis;IVA))は、パーミュテーション問題を生じない分離手法として知られている（非特許文献２）。IVAにおいては、各音源のフレーム毎の複素スペクトルをベクトル型変数Y_k,τ＾=(Y_k,τ,1,...,Y_k,τ,ω,...,Y_k,τ,N)^Tとして扱い(ただし、kは音源、τは時間フレーム、ωは周波数を表すインデックスである。)、そのノルムである

が優ガウス分布に従うという仮定のもとで最尤となる各帯域の分離フィルタが推定される。また、優ガウス分布の具体的な分布として時変ガウス分布(分散が時刻ごとに変化することを許容したガウス分布)を仮定したIVAを実現する方法も提案されている（非特許文献３）。

IVAによるブラインド音源分離の問題は、観測信号の時間周波数表現X_Tω＾=(X_1Tω,...,X_MTω)^Tに対し、分離信号

の各要素が統計的に独立になるように分離行列W_ω＾の推定する問題として定式化される。ただし、Ｙ_τω＾＝(Y_1τω,...,Y_kτω,...,Y_Kτω)^Tであり、各音源の時間周波数成分を表わしている。

A. Hyv¨arinen, J. Karhumen, and E. Oja, Independent Component Analysis, John Wiley & Sons, 2001. T. Kim, T. Eltoft and T. Lee, "Independent Vector Analysis: An Extension of ICA to Multivariate Components," Proc. ICA, pp.165−172, 2006. 小野拓磨、小野順貴、嵯峨山茂樹、"音源のアクティベーションを事前情報とした独立ベクトル析による音源分離、"日本音響学会秋季研究発表会講演集、pp.613−614,Sep. 2011.

以上のように、IVAでは、複素スペクトルをベクトル化したもののノルムが優ガウス分布に従うと仮定されるが、優ガウス分布として時変ガウス分布(フレームごとに分散が変化することを許容したガウス分布)を仮定した場合、この仮定は、各音源に対してどの周波数でもパワーの時間包絡が等しいと仮定したことになっている。しかし、一般に対象となる音源は音楽や音声など調波性という特別な性質を持つ場合がある。この場合、調波成分間の谷の部分での振幅包絡は必ずしも他の周波数と一致しない一方、基本周波数や倍音間のパワーの依存関係は音源分離の重要な手掛かりになり得る。また、仮定に適合しない信号（ピッチや周期といった調波性をもつ音源信号）では分離の精度に限界がある。

また、IVAでは、Y_Tω＾の振幅包絡がどの周波数間でも類似する、という仮定ができるだけ満たされるようにW_ω＾が決定されるため、音源が調波構造を有する場合にはこの仮定が大きく崩れ、調波成分の谷間に該当する箇所に他音源の成分が混入したり、パーミュテーション不整合が起こることがあった。よって、もし各音源が調波構造を有するならば、音源モデルには、調波構造形のパワースペクトル密度を仮定した方が良いはずである。しかし、そのためには音源の基本周波数の情報が必要となるが、音源の基本周波数の情報は通常観測することができない。

本発明は、上記の事情を鑑みてなされたもので、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる音響信号解析装置、方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る音響信号解析装置は、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘ_mτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘを出力する時間周波数解析手段と、Ｍ個の音源ｋの時間周波数成分Ｙ_kτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ^(l) _ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π^(l) _kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ² _kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘ_τω（＝（Ｘ_1τω，・・・，Ｘ_Mτω））に作用させて音源信号の時間周波数成分Ｙ_τω＾（＝（Ｙ_1τω，・・・，Ｙ_Kτω））を得るための分離行列Ｗ_ω＾の各々の初期値を設定するパラメータ初期値設定手段と、（k、τ、ω、ｚ_kτ）の全ての組み合わせにおける、ｚ_kτとσ² _kγが与えられたときの、λ^(ｚkτ) _ω・σ² _kγを分散とするガウス分布で表されたＹ_kτωの確率密度関数、及びπ_kτが与えられたときのｚ_kτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π^(l) _kτの事前確率、及び各周波数ωに対する前記分離行列Ｗ_ω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾の事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾を更新するパラメータ更新手段と、各周波数ωの前記分離行列Ｗ_ω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新する音源信号推定値更新手段と、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う終了判定手段と、を含んで構成されている。

本発明に係る音響信号解析方法は、時間周波数解析手段によって、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘ_mτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘ＾を出力し、パラメータ初期値設定手段によって、Ｍ個の音源ｋの時間周波数成分Ｙ_kτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ^(l) _ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π^(l) _kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ² _kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘ_τω＾（＝（Ｘ_1τω，・・・，Ｘ_Mτω））に作用させて音源信号の時間周波数成分Ｙ_τω＾（＝（Ｙ_1τω，・・・，Ｙ_Kτω））を得るための分離行列Ｗ_ω＾の各々の初期値を設定し、パラメータ更新手段によって、（k、τ、ω、ｚ_kτ）の全ての組み合わせにおける、ｚ_kτとσ² _kγが与えられたときの、λ^(ｚkτ) _ω・σ² _kγを分散とするガウス分布で表されたＹ_kτωの確率密度関数、及びπ_kτが与えられたときのｚ_kτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π^(l) _kτの事前確率、及び各周波数ωに対する前記分離行列Ｗ_ω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ωの事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾を更新し、音源信号推定値更新手段によって、各周波数ωの前記分離行列Ｗ_ω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新し、終了判定手段によって、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う。

本発明に係るプログラムは、上記の音響信号解析装置の各手段としてコンピュータを機能させるためのプログラムである。

以上説明したように、本発明の音響信号解析装置、方法、及びプログラムによれば、各周波数に対する各時間フレームにおいて調波構造を有するパワースペクトルテンプレートの各々が音源毎に選択される確率を要素にもつ三次元配列Π＾を用いて、観測時間周波数成分の三次元配列Ｘが与えられたときの三次元配列Π＾、各音源のパワーを要素にもつ二次元配列Σ＾、及び各周波数ωの分離行列Ｗ_ω＾の事後確率を表す目的関数を最大化するように、三次元配列Π＾、二次元配列Σ＾、及び各周波数ωの分離行列Ｗ_ω＾を更新することを繰り返すことにより、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる、という効果が得られる。

本発明の実施の形態に係る音響信号解析装置の構成を示す概略図である。 パワースペクトルテンプレートを示す図である。 本発明の実施の形態に係る音響信号解析装置における音響信号解析処理ルーチンの内容を示すフローチャートである。 （Ａ）ＳＤＲによる部屋E2Aにおける分離性能の評価結果を示すグラフ、及び（Ｂ）ＳＤＲによる部屋JR2における分離性能の評価結果を示すグラフである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する手法では、スペクトルテンプレートが時変ガウス分布における分散パラメータに組み込まれた音源の生成モデルを立てて、それをもとに、各音源において各フレームでどのスペクトルテンプレートが選ばれるべきかということと、各帯域の分離フィルタとを同時に推定する。

＜発明の原理＞
まず、本発明の原理について説明する。まず、音源の生成モデルについて説明する。

各フレームにおける基本周波数を潜在変数と見なし、潜在変数(基本周波数インデックスz_kT)に応じて調波構造形のパワースペクトルテンプレートが一つ選択され、そのパワースペクトルをもとに音源の周波数成分が決定されるという生成プロセスを仮定する。これにより、従来のIVAにおいて仮定される音源モデルに比べ、調波性を有する音源を適切にモデル化できるはずであろうと考えられる。

l=1,...,Lを調波構造テンプレートのインデックスとし、λ^(l) _ωを特定の基本周波数の調波構造テンプレートとする（図２参照）。以上より、音源のGauss性を仮定すると、z_kTが与えられたもとでの音源の生成モデルは、以下の（３）式で表される。

ただし、N_C(・;μ,σ²)は平均μ、分散σ²の複素Gauss分布を表し、σ_kT ²はフレームτにおける音源kのパワーである。一方、z_kτ=lが選択される確率をπ_lkτとすると、潜在変数z_kτは、以下の（４）式に示す離散分布で表される。

更に、π_kτをスパースにするように誘導する目的でπ_kτのハイパー事前分布を以下の（５）式に示すDirichlet分布と仮定する。

以上の音源の生成モデルを、上記（２）式に組み込んだIVAの問題は、以下の（６）式に示す目的関数を最大化する問題に帰着する。

ただし、X＾={X_mτω},Π＾={π^(l) _kτ},Σ＾={σ² _kτ},W＾={W_ω},Θ＾={Π,Σ,W}である。Tは、時間フレーム数である。なお、記号に付された「＾」は、当該記号が行列または多次元配列またはベクトルであることを表わしている。

上記の（６）式の目的関数の大域最適解は解析的に得ることはできないが、補助関数法により局所解を効率的に探索することができる。ここで、Σ_lγ ^(l) _kτω=1となる補助変数γ_kτωを用いれば、以下の（７）式のように補助関数を設計できる。

＜補助変数の更新＞
補助変数についての更新式はΣ_lγ^(l) _kτω=1の制約のもと、∂Q/∂_γ ^(l) _kτω=0を解けば以下の（８）式のとおりに得られる。

ただし、λはテンプレートであり、ｌは、テンプレートのインデックスを表す。

＜パラメータの更新＞
π^(l) _kτとσ² _kτの更新式については、∂Q/∂π^(l) _kτ=0,∂Q/∂σ² _kτ=0を解くことにより以下の（９）式、（１０）式のとおりに得られる。

W_ω＾については、以下の（１１）式〜（１４）式のように1行ずつ更新する。

ただし、e_k^は、k番目のベクトル要素が1でそれ以外の要素は0の単位ベクトル（＝[0,...,1,...,0]^T）である。

＜音源信号推定値の更新＞
以上で更新された分離行列を用いて音源信号の推定を以下の（１５）式のとおりに更新する。

ただし、Y_τ,ω＾=(Y_1,τ,ω,...,Y_K,τ,ω)^Tである。

＜システム構成＞
次に、Ｍ個（Ｍ≧２）のマイクロホンから得られた音響信号を解析して、既知のＫ個（Ｋ＜Ｍ）の音源信号に分離する音響信号解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。

図１に示すように、本発明の実施の形態に係る音響信号解析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。

音響信号解析装置１００は、入力部１０と、演算部２０と、記憶部３０と、出力部４０とを備えている。

入力部１０により、Ｍ個のマイクロホンから出力された音響信号（多チャンネル信号）の時系列データが入力される。記憶部３０は、入力部１０により入力された音響信号の時系列データを記憶する。また、記憶部３０は、後述する各処理での結果を記憶すると共に、各基本周波数の調波構造テンプレートλ^(l) _ω（l=1,...,L、ω＝１,・・・,Ｎ）、パラメータβを記憶している。

演算部２０は、時間周波数解析部２１と、初期設定部２２と、補助変数更新部２３と、パラメータ更新部２４と、時間周波数成分推定部２５と、終了判定部２６と、信号変換部２７とを備えている。また、パラメータ更新部２４は、確率更新部２４１と、時変ゲイン更新部２４２と、分離行列更新部２４３とを備えている。

時間周波数解析部２１は、各マイクロホンの時系列信号としての観測された音響信号を入力として、時間周波数成分（観測時間周波数成分）Ｘ_mτω（ｍ＝１，・・・，Ｍ、ω＝１,・・・,Ｎ,τ＝１,・・・,Ｔは、それぞれマイクロホン、周波数、時間フレームに対応するインデックスを示す。）を各（ｍ，ω，τ）の要素にもつ三次元配列Ｘを計算する。また、計算した時間周波数成分Ｘ_τωを、記憶部３０に記憶しておく。より詳細には、時間周波数解析部２１は、各マイクロホンｍについて、当該マイクロホンの音響信号の時系列データを入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより時間周波数成分Ｘ_mτωを計算し、時間周波数成分Ｘ_mτωを格納した行列（振幅スペクトログラム）Ｘ＾＝（Ｘ_mτω）_M×N×Ｔを出力する。なお、時間周波数成分Ｘ_mτωは、ウェーブレット変換を用いて計算してもよい。

初期設定部２２は、後述する処理で用いる各パラメータπ^(l) _kτ、σ² _kτ、W_ω＾の初期値を設定する。初期値として、乱数を用いて適当な値を設定すればよい。また、初期設定部２２は、W_ω＾の初期値と時間周波数成分Ｘ_mτωとに基づいて、上記（２）式に従って、分離信号Ｙτω＾の初期値を計算する。

補助変数更新部２３は、（ｌ、k、τ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているπ^(l) _kτ、σ² _kτ、Ｙ_kτω、λ^(l) _ωに基づいて、上記（８）式に従って、観測時間周波数成分Ｘ_mτωが各インデックスｌのパワースペクトルテンプレートλ(l)ωに帰属する確率を示す補助変数γ^(l) _kτωを更新し、記憶部３０に格納する。

確率更新部２４１は、（ｌ、k、τ）の全ての組み合わせの各々について、記憶部３０に記憶されているγ^(l) _kτω、βに基づいて、上記（９）式に従って、パワースペクトルテンプレートλ^(l) _ωが選択される確率π^(l) _kτを更新し、記憶部３０に格納する。

時変ゲイン更新部２４２は、（k、τ）の全ての組み合わせの各々について、記憶部３０に記憶されているγ^(l) _kτω、Ｙ_kτω、λ^(l) _ωに基づいて、上記（１０）式に従って、フレームτにおける音源kのパワーσ_kT ²を更新し、記憶部３０に格納する。

分離行列更新部２４３は、（k、τ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているγ^(l) _kτω、λ^(l) _ω、σ_kτ ²に基づいて、上記（１１）式に従って、パラメータ~σ_kτω ²を更新し、記憶部３０に格納する。また、分離行列更新部２４３は、（k、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているＸ_τω、~σ_kτω ²に基づいて、上記（１２）式に従って、中間パラメータＶ_kωを更新し、記憶部３０に格納する。

そして、分離行列更新部２４３は、（k、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているW_ω＾、Ｖ_kωに基づいて、上記（１３）式、（１４）式に従って、ｗ_kω＾（音源ｋの音響信号の周波数成分ωに対する各マイクロホンｍの重みを要素にもつ重みベクトル）を更新し、記憶部３０に格納する。

時間周波数成分推定部２５は、（τ、ω）の全ての組み合わせの各々について、記憶部３０に記憶されているW_ω＾、Ｘ_τω＾に基づいて、上記（１５）式に従って、各音源の時間周波数成分を表わすベクトルＹ_τω＾＝(Y_1τω,...,Y_kτω,...,Y_Kτω)^Tを更新し、記憶部３０に格納する。

終了判定部２６は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、補助変数更新部２３、パラメータ更新部２４、及び時間周波数成分推定部２５の各処理を繰り返す。終了判定部２６は、終了条件を満足したと判定した場合には、信号変換部２７による処理に移行する。信号変換部２７は、記憶部３０に記憶されているＹ_τω＾を、各音源の音源信号に変換し、出力部４０により、音源信号を出力する。

終了条件としては、繰り返し回数ｓが予め定めた回数Ｓに達したことを用いればよい。なお、s-1回目のパラメータを用いたときの目的関数の値とs回目のパラメータを用いたときの目的関数の値との差が、予め定めた閾値よりも小さくなったことを、終了条件として用いてもよい。

＜音響信号解析装置の作用＞
次に、本実施の形態に係る音響信号解析装置１００の作用について説明する。まず、解析対象の信号として各マイクロホンからの音響信号の時系列データが音響信号解析装置１００に入力され、記憶部３０に格納される。そして、音響信号解析装置１００において、図３に示す音響信号解析処理ルーチンが実行される。

まず、ステップＳ１０１において、記憶部３０から、マイクロホン毎に、各時間フレームτ内の音響信号を読み込み、当該音響信号に対して、短時間フーリエ変換を用いた時間周波数分析を行った結果から、観測時間周波数成分Ｘ_mτωを各（ｍ，τ，ω）の要素にもつ三次元配列Ｘ＾を生成して、記憶部３０に記憶する。

そして、ステップＳ１０２において、乱数を用いて、パラメータΘ＾＝｛Π＾、Σ＾、Ｗ＾｝の初期値を設定して、記憶部３０に記憶すると共に、各音源信号の時間周波数成分Ｙ_kτωを各（ｋ，τ，ω）の要素にもつ三次元配列Ｙ＾を生成して、記憶部３０に記憶する。

次にステップＳ１０３では、上記ステップＳ１０２で設定されたパラメータΠ＾、Σ＾、Ｙ＾、又は後述するステップＳ１０４、Ｓ１０５、Ｓ１０７で更新されたパラメータΠ＾、Σ＾、Ｙ＾に基づいて、上記（８）式に従って、補助係数γ^(l) _kτωを各（ｌ，k，τ，ω）の組み合わせについて算出して、記憶部３０に格納する。

そして、ステップＳ１０４では、上記ステップＳ１０３で更新された補助係数γ^(l) _kτωに基づいて、上記（９）式に従って、パワースペクトルテンプレートλ^(l) _ωが選択される確率π^(l) _kτを各（ｌ，k，τ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１０５では、上記ステップＳ１０２で設定されたパラメータＹ＾、又はステップＳ１０７で更新されたパラメータＹ＾と、上記ステップＳ１０３で更新された補助係数γ^(l) _kτωに基づいて、上記（１０）式に従って、時間フレームτにおける音源kのパワーσ_kτ ²を各（k，τ）の組み合わせについて算出して、記憶部３０に格納する。

次のステップＳ１０６では、上記ステップＳ１０２で設定されたパラメータＷ＾、又は前回のステップＳ１０６で更新されたパラメータＷ＾と、上記ステップＳ１０１で生成された三次元配列Ｘ＾と、上記ステップＳ１０５で更新されたパラメータΣ＾と、上記ステップＳ１０３で更新された補助係数γ^(l) _kτωと、基づいて、上記（１１）式〜（１４）式に従って、ｗ_kωを各（k，ω）の組み合わせについて算出して、記憶部３０に格納する。

そして、ステップＳ１０７では、上記ステップＳ１０６で更新されたパラメータＷと、上記ステップＳ１０１で生成された三次元配列Ｘ＾とに基づいて、上記（１５）式に従って、各音源信号の時間周波数成分Ｙ_kτωを各（ｋ，τ，ω）の要素にもつ三次元配列Ｙを算出して、記憶部３０に記憶する。

次のステップＳ１０８では、終了条件として、繰り返し回数ｓが、Ｓに到達したか否かを判定し、繰り返し回数ｓがＳに到達していない場合には、終了条件を満足していないと判断して、上記ステップＳ１０３へ戻り、上記ステップＳ１０３〜ステップＳ１０７の処理を繰り返す。一方、繰り返し回数ｓがＳに到達した場合には、終了条件を満足したと判断し、ステップＳ１０９で、上記ステップＳ１０７で最終的に更新された三次元配列Ｙに基づいて、各音源の音源信号を算出し、出力部４０により出力して、音響信号解析処理ルーチンを終了。

＜実験結果＞
次に、本実施の形態に係る手法の有用性を示す目的で、単旋律楽器を音源として用いたシミュレーションによる実験を行った結果について説明する。以下の表１に実験条件を示す。

予め録音された残響時間の異なる二種類の室内インパルス応答（参考文献（S. Nakamura, K. Hiyane, F. Asano, T. Nishiura and T. Yamada, “Acoustical Sound Database in Real Environments for Sound Scene Understanding and Hands-Free Speech Recognition,” Proc. LREC, pp. 965−968, 2000.）を参照）を単旋律楽器音(トランペット、サックス)とそれぞれ畳み込み加算することで観測信号とした。部屋E2A(残響時間T60=300ms)では、{−20°,0°}の方向から、部屋JR2(T60=470ms)では、{−30°,30°}の方向から音源が到来するとした。

本実施の形態に係る手法（以下、提案法と称する。）に用いた重みテンプレートは220Hz(A3)から半音ごとに24種類の倍音成分が1/√nで減衰する調波構造テンプレートと全ての周波数で同じ重みのテンプレート1種類の合計25 種類を用いた。提案法の分離行列の初期値として、時変Gauss分布に基づく独立ベクトル分析により分離行列を40回更新した値を用い、その後に提案法を40回反復し、復元信号を得た。提案法はπ^(l) _kTの事前分布としてDirichlet分布を考慮した手法(Prop(dir))と事前分布を一様にした手法(Prop(flat))の2種類を行った。Dirichlet分布を仮定した実験では、フレーム長が大きくなると事前分布の効用が小さくなるため効果を調節する目的でp(π_kT)を周波数bin数乗した。また、提案法40回反復の間にp(π_kT)のパラメータβを1から0.95まで徐々に小さくするアニーリング処理を行った。このとき更新の間にπ^(l) _kτが負になることがあったが、その場合十分小さい正数に置き換えた。

比較対象として、音源の生成モデルに時変Gauss分布を仮定した独立ベクトル分析による手法(IVA(TVG))音源の生成モデルにラプラス分布を仮定した独立ベクトル分析による手法(IVA(lap))（上記の非特許文献２参照）、独立成分分析による手法(ICA)(参考文献（H. Sawada, R. Mukai, S. Araki and S. Makino, “A Robust Approach to the Permutation Problem of Frequency-domain Blind source Separation,” Proc. ICASSP, pp. 381−384, 2003.）を参照)の3種類をそれぞれ分離行列の更新に80回行った。

SDR（参考文献（E. Vincent, R. Gribonval and C. F´evotte, “Performance Measurement in Blind Audio Source Separation,” IEEE Trans. ASLP, pp.1462−1469, 2006.）を参照）による評価の結果を、図４に示す。ＳＤＲによる部屋E2A(残響時間T60=300ms)における分離性能の評価結果を図４（Ａ）に示し、ＳＤＲによる部屋JR2(T60=470ms)における分離性能の評価結果を図４（Ｂ）に示す。2種類の環境どちらにおいても従来法に比べ提案法が優位な結果を得られた。特に、部屋E2AにおいてDirichlet事前分布を考慮することで、6dB程度の改善が見られた。

以上説明したように、本発明の実施の形態に係る音響信号解析装置によれば、各周波数及び各時間フレームにおいて調波構造を有するパワースペクトルテンプレートの各々が音源毎に選択される確率を要素にもつ三次元配列Π＾を用いて、観測時間周波数成分の三次元配列Ｘが与えられたときの三次元配列Π＾、各音源のパワーを要素にもつ二次元配列Σ＾、及び各周波数ωの分離行列Ｗ_ω＾の事後確率を表す目的関数を最大化するように、三次元配列Π＾、二次元配列Σ＾、及び各周波数ωの分離行列Ｗ_ω＾を更新することを繰り返すことにより、複数のマイクロホンから出力された音響信号の時系列データから、音源毎の音源信号に精度よく分離することができる。

また、調波構造をもつパワースペクトルのテンプレートを複数個用意しておき、スペクトルテンプレートが時変ガウス分布における分散パラメータに組み込まれた音源の生成モデルを立てて、それをもとに、各音源において各時間フレームでどのスペクトルテンプレートが選ばれるべきかということと、各帯域の分離フィルタとを同時に推定することにより、音源信号が調波構造をもつ場合に、音源分離の精度を向上させることができる。

また、上記の非特許文献３とのアルゴリズムにおける相違点は、（１）γ^(l) _kτω、π^(l) _kTというパラメータとその更新式(上記（８）式、（９）式)が新たに加わったこと、（２）各音源のパワーσ² _kτ が、上記（１０）式のようにγ^(l) _kTωを用いて算出されること、（３）中間変数V_kωが、上記（１１）式に従って算出される各音源の各時刻におけるパワースペクトル推定値~σ² _kτω(周波数ωごとに異なる値をとりうる点が従来技術との重要な差異)を用いて上記（１２）式のように更新されることである。

また、複数のマイクロホンで取得した音響信号から混在する個々の音源信号を分離することができるため、ハンズフリーテレビ会議システムや会議録コンテンツの自動作成システムなどの応用が期待される。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上述の音響信号解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

１０ 入力部
２０ 演算部
２１ 時間周波数解析部
２２ 初期設定部
２３ 補助変数更新部
２４ パラメータ更新部
２５ 時間周波数成分推定部
２６ 終了判定部
２７ 信号変換部
３０ 記憶部
４０ 出力部
１００ 音響信号解析装置
２４１ 確率更新部
２４２ 時変ゲイン更新部
２４３ 分離行列更新部

Claims (5)

1. Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘ_mτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘを出力する時間周波数解析手段と、
   Ｍ個の音源ｋの時間周波数成分Ｙ_kτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ^(l) _ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π^(l) _kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ² _kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘ_τω（＝（Ｘ_1τω，・・・，Ｘ_Mτω））に作用させて音源信号の時間周波数成分Ｙ_τω（＝（Ｙ_1τω，・・・，Ｙ_Kτω））を得るための分離行列Ｗ_ω＾の各々の初期値を設定するパラメータ初期値設定手段と、
   （k、τ、ω、ｚ_kτ）の全ての組み合わせにおける、ｚ_kτとσ² _kγが与えられたときの、λ^(ｚkτ) _ω・σ² _kγを分散とするガウス分布で表されるＹ_kτωの確率密度関数、及びπ_kτが与えられたときのｚ_kτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π^(l) _kτの事前確率、及び各周波数ωに対する前記分離行列Ｗ_ω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾の事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾を更新するパラメータ更新手段と、
   各周波数ωの前記分離行列Ｗ_ω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新する音源信号推定値更新手段と、
   予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う終了判定手段と、
   を含む音響信号解析装置。
2. 前記目的関数を、（ｋ，τ、ω）の全ての組み合わせについて観測時間周波数成分Ｘ_mτωが各インデックスｌのパワースペクトルテンプレートλ^(l) _ωに帰属する確率を示す補助変数γ^(l) _kτωを用いた補助関数とし、
   前記パラメータ更新手段は、
   前記三次元配列Π＾、二次元配列Σ＾、及び前記複数のパワースペクトルテンプレートλ^(l) _ωに基づいて、（ｋ，τ、ω、ｌ）の全ての組み合わせの各々について前記補助変数γ^(l) _kτωを更新する補助変数更新手段と、
   前記補助変数γ^(l) _kτωに基づいて、前記三次元配列Π＾を更新する確率更新手段と、
   前記補助変数γ^(l) _kτω及び前記複数のパワースペクトルテンプレートλ^(l) _ωに基づいて、前記二次元配列Σを更新するパワー更新手段と、
   前記補助変数γ^(l) _kτω、前記二次元配列Σ＾、前記三次元配列Ｘ＾、及び前記複数のパワースペクトルテンプレートλ^(l) _ωに基づいて、各周波数ωの前記分離行列Ｗ_ω＾を更新する分離行列更新手段と、
   を含む請求項１記載の音響信号解析装置。
3. 前記確率π^(l) _kτの事前分布を、ディレクレ分布とした請求項１又は２記載の音響信号解析装置。
4. 時間周波数解析手段によって、Ｍ個（Ｍは２以上の整数）のマイクロホンｍから出力される音響信号の時系列データを入力として、観測時間周波数成分Ｘ_mτω（ｍはマイクロホン、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｘを出力し、
   パラメータ初期値設定手段によって、Ｍ個の音源ｋの時間周波数成分Ｙ_kτω（ｋは音源、τは時間フレーム、ωは周波数のインデックスである。）を要素にもつ三次元配列Ｙ＾、各音源ｋについて各時間フレームτにおいて調波構造を有するパワースペクトルテンプレートλ^(l) _ω（ｌは、パワースペクトルテンプレートのインデックスである。）が選択される確率π^(l) _kτを要素にもつ三次元配列Π＾、各時間フレームτにおける各音源ｋのパワーσ² _kγを要素にもつ二次元配列Σ＾、及び各周波数ωについて、Ｘ_τω（＝（Ｘ_1τω，・・・，Ｘ_Mτω））に作用させて音源信号の時間周波数成分Ｙ_τω（＝（Ｙ_1τω，・・・，Ｙ_Kτω））を得るための分離行列Ｗ_ω＾の各々の初期値を設定し、
   パラメータ更新手段によって、（k、τ、ω、ｚ_kτ）の全ての組み合わせにおける、ｚ_kτとσ² _kγが与えられたときの、λ^(ｚkτ) _ω・σ² _kγを分散とするガウス分布で表されたＹ_kτωの確率密度関数、及びπ_kτが与えられたときのｚ_kτの確率、（k、τ、ｌ）の全ての組み合わせに対する前記確率π^(l) _kτの事前確率、及び各周波数ωに対する前記分離行列Ｗ_ω＾の行列式を用いて表された、前記三次元配列Ｘ＾が与えられたときの前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾の事後確率を表す目的関数を最大化するように、前記三次元配列Π＾、前記二次元配列Σ＾、及び各周波数ωの前記分離行列Ｗ_ω＾を更新し、
   音源信号推定値更新手段によって、各周波数ωの前記分離行列Ｗ_ω＾及び前記三次元配列Ｘ＾に基づいて、前記三次元配列Ｙ＾を更新し、
   終了判定手段によって、予め定められた終了条件を満たすまで、前記パラメータ更新手段による更新、及び前記音源信号推定値更新手段による更新を繰り返し行う
   音響信号解析方法。
5. 請求項１〜請求項３の何れか１項に記載の音響信号解析装置の各手段としてコンピュータを機能させるためのプログラム。