以下に、本願に係る信号分析装置、信号分析方法及び信号分析プログラムの実施形態を図面に基づいて詳細に説明する。また、本発明は、以下に説明する実施形態により限定されるものではない。なお、以下では、ベクトル、行列又はスカラーであるAに対し、“^A”と記載する場合は「“A”の直上に“^”が記された記号」と同じであるとする。ベクトル、行列又はスカラーであるAに対し、“~A”と記載する場合は「“A”の直上に“~”が記された記号」と同じであるとする。また、以下では、信号分離は、信号分析の下位概念とする。
[第1の実施形態]
まず、第1の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。なお、第1の実施形態においては、2個の源信号が混在する状況において、それぞれ異なる位置でマイクロホンなどのセンサにより取得されたI個(Iは2以上の整数)の観測信号yi(k)(i=1,・・・,Iは観測信号のインデックス、kはサンプル点のインデックス)と、源信号のサンプル点ごとの分散(パワー)をモデル化するパラメータである分散パラメータvj(k)(j=1,2は源信号のインデックス)と、源信号の空間的特性をモデル化するパラメータである空間共分散行列Rjと、が信号分析装置に入力されるものとする。
なお、本第1の実施形態における「信号」の例としては、音、脳波、無線信号などが挙げられる。本第1の実施形態における「観測信号」とは、複数のセンサ(センサアレイ)により取得された信号を指す。センサアレイの例としては、マイクロホンアレイ、脳波計・脳磁計、アンテナアレイなどが挙げられる。また、本第1の実施形態における「源信号」とは、観測信号の構成要素である信号を指す。観測信号は2個の源信号の和としてモデル化される。
また、本第1の実施形態における「信号分離」の問題設定は、2個の目的信号(たとえば音声)が混在する状況においてこれらに対応する成分を推定する問題設定(たとえば音源分離)に加えて、目的信号(たとえば音声)とノイズ(たとえばテレビから流れる音楽や雑踏における喧噪)が混在する状況においてこれらに対応する成分(または目的信号に対応する成分のみ)を推定する問題設定(ノイズ除去とも呼ばれる)をも含むものとする。前者の場合における源信号は、目的信号(たとえば音声)である。後者の場合における源信号は、目的信号(たとえば音声)とノイズ(たとえばテレビから流れる音楽や雑踏における喧噪)とである。なお、「雑踏における喧噪」は、多数のノイズが混ざったものであるが、この場合まとめて1つの源信号とみなす。
[信号分析装置の構成及び処理]
まず、図1及び図2を用いて、第1の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図1は、第1の実施形態に係る信号分析装置の構成の一例を示す図である。図2は、第1の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。第1の実施形態に係る信号分析装置1は、例えば、ROM(Read Only Memory)、RAM(Random Access Memory)、CPU(Central Processing Unit)等を含むコンピュータ等に所定のプログラムが読み込まれて、CPUが所定のプログラムを実行することで実現される。
図1に示すように、信号分析装置1は、観測信号ベクトル作成部10、信号パラメータ分析部20、同時無相関化部30、シングルチャネルウィーナーフィルタ部40及び同時無相関化逆変換部50を有する。
まず、信号分析装置1の各部の概要について説明する。観測信号ベクトル作成部10は、入力された観測信号yi(k)を取得し(ステップS10)、(4)式に示す、すべての観測信号からなるI次元縦ベクトルである観測信号ベクトルをサンプル点ごとに作成する(ステップS11)。
信号パラメータ分析部20は、入力された空間共分散行列Rjを取得し、一般化固有値問題を解くことにより、一般化固有値行列Λ及び一般化固有ベクトル行列Pを計算する(ステップS12)。
同時無相関化部30は、観測信号ベクトル作成部10からの観測信号ベクトルy(k)と、信号パラメータ分析部20からの一般化固有ベクトル行列Pと、に基づいて、同時無相関化された観測信号ベクトルy´(k)を計算する(ステップS13)。具体的には、同時無相関化部30は、(5)式に示すように、観測信号ベクトルy(k)に、一般化固有ベクトル行列Pのエルミート転置PHを掛けることにより、同時無相関化された観測信号ベクトルy´(k)を計算する。
ここで、「´(上付きのプライム)」は、同時無相関化を施された変数であることを表す。なお、同時無相関化とは、観測信号ベクトルy(k)に含まれる、2個の源信号に対応する成分(源信号成分)を同時に無相関化する線型変換を指す。
シングルチャネルウィーナーフィルタ部40は、入力された分散パラメータvj(k)を取得し、信号パラメータ分析部20からの一般化固有値行列Λと、同時無相関化部30からの同時無相関化された観測信号ベクトルy´(k)と、に基づいて、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を計算する(ステップS14)。具体的には、シングルチャネルウィーナーフィルタ部40は、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を(6),(7)式により計算する。
ただし、ベクトルαに対し、[α]lはαの第l要素を表し、行列Aに対し、[A]lmはAの(l,m)要素を表す。また、変数の上の∧(ハット)は、推定値であることを表す。
同時無相関化逆変換部50は、信号パラメータ分析部20からの一般化固有ベクトル行列Pと、シングルチャネルウィーナーフィルタ部40からの同時無相関化された源信号成分x´j(k)の推定値^x´j(k)と、に基づいて、源信号成分xj(k)の推定値^xj(k)を計算し、出力する(ステップS15)。具体的には、同時無相関化逆変換部50は、(8)式に示すように、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)に対し、同時無相関化部30で掛けた一般化固有ベクトル行列Pのエルミート転置PHの逆行列(PH)−1を掛けることにより、源信号成分xj(k)の推定値^xj(k)を計算し、出力する。
[第1の実施形態の処理の背景]
次に、第1の実施形態の処理の背景について説明する。なお、以下では、観測信号が複素数で表される場合について説明する(観測信号が実数で表される場合の説明は、ほぼ同様だから省略する)。第1の実施形態においては、観測信号ベクトルy(k)は、2個の源信号に対応する成分(源信号成分)x1(k),x2(k)の和として、次に示す(9)式でモデル化される。
源信号成分xj(k)は、平均が0、共分散行列がSj(k)であるとし、2個の源信号成分x1(k),x2(k)は互いに無相関であるとする。2個の源信号成分x1(k),x2(k)は、ともに目的信号(たとえば音声)に対応する源信号成分であっても良いし、一方が目的信号(たとえば音声)に対応する源信号成分で他方がノイズ(たとえばテレビからの雑音や雑踏における喧噪)に対応する源信号成分であっても良い。
観測信号ベクトルy(k)から源信号成分xj(k)を推定する方法としては、マルチチャネルウィーナーフィルタに基づく方法がある。この方法では、源信号成分xj(k)の推定値^xj(k)を(10)式により計算する。
すなわち、この方法では、観測信号ベクトルy(k)に対し、I×Iの行列であるマルチチャネルウィーナーフィルタSj(k)(S1(k)+S2(k))−1を掛けることにより、源信号成分xj(k)の推定値^xj(k)を計算する。
[従来の信号分析方法]
ここで、従来の信号分析方法(例えば、非特許文献1に記載の信号分析方法)では、(10)式に基づいて源信号成分xj(k)の推定値^xj(k)を計算する。その際、源信号成分xj(k)の共分散行列Sj(k)は、次の(11)式によりモデル化される。なお、(11)式におけるEは、期待値演算である。
ここで、vj(k)は、源信号のサンプル点ごと(kごと)の分散をモデル化するパラメータである分散パラメータであり、Rjは、源信号の空間的特性をモデル化するパラメータである空間共分散行列である。従来の信号分析方法では、(10)式に(11)式を代入して得られる(12)式を用いて、源信号成分xj(k)の推定値^xj(k)を計算していた。
上述の従来の信号分析方法の問題点として、膨大な計算量を要するため、大規模なデータセットの処理(サンプル点kの数が膨大になる)や、リアルタイム処理(短時間での処理が要求される)、計算能力の低い機器などへの適用が困難だった。
実際、従来技術では、源信号成分xj(k)の推定値^xj(k)を計算する際に、逆行列計算(具体的には、(v1(k)R1+v2(k)R2)−1の計算)をサンプル点毎に(つまり多数回)行う必要があった。逆行列計算は、I次行列に対しO(I3)という大きな計算量を要する演算である。また、サンプル点の数は、典型的には数百〜数千万程にもなる。そのため、従来の信号分析方法では、源信号成分xj(k)の推定値^xj(k)を計算する際に、膨大な計算量を必要としていた。
[第1の実施の形態に係る信号分析方法]
これに対し、第1の実施形態では、2個の源信号成分x1(k),x2(k)の同時無相関化に基づいて、源信号成分xj(k)の推定値^xj(k)を少ない計算量で計算する方法を提供する。以下、第1の実施形態における考え方について説明する。
各サンプル点kに対し信号α(k)の共分散行列E(α(k)α(k)H)が対角行列であるとき、α(k)は無相関であるという。これは、簡単に言えば、「ベクトルα(k)の異なる要素同士は相関を持たない」ということである。いま仮に、(11)式の源信号成分xj(k)の共分散行列vj(k)Rj(jは1または2)がいずれも対角行列であるとすると、(12)式における逆行列演算(v1(k)R1+v2(k)R2)−1は単なる対角要素ごとの逆数演算により実現できるため、少ない計算量で計算できる。しかしながら、通常、源信号成分xj(k)の異なる要素同士は互いに強い相関を持つため、その共分散行列vj(k)Rjは対角行列ではない。そこで、本第1の実施形態では、2個の源信号成分x1(k),x2(k)を同時無相関化することを考える。
2個の源信号成分x1(k),x2(k)の同時無相関化とは、(13)式及び(14)式に示すように、x1(k),x2(k)に正則行列Vを掛ける線型変換を施して、線型変換後の源信号成分x´1(k),x´2(k)が同時に無相関になるようにする(すなわち、x´1(k),x´2(k)の共分散行列S´1(k),S´2(k)が同時に対角行列になるようにする)ことを指すものとする。この同時無相関化により、源信号成分xj(k)の共分散行列vj(k)Rjはいずれも対角行列に変換される。したがって、同時無相関化により、(12)式における逆行列演算(v1(k)R1+v2(k)R2)−1は、対角行列に対する逆行列演算に変わり、単なる対角要素ごとの逆数演算により実現できるようになるため、少ない計算量で計算できるようになる。
x´1(k),x´2(k)の共分散行列S´1(k),S´2(k)を計算すると、次の(15)〜(19)式のようになる。
同時無相関化の条件は、これらが同時に対角行列になること、すなわち以下の(20)式及び(21)式が同時に成り立つことである。なお、Λ1,Λ2は或る対角行列である。
すなわち、同時無相関化の目標は、空間共分散行列R1,R2に対し、左から正則行列Vを、右からVのエルミート転置VHを、それぞれ掛けたときに、これらが同時に対角行列になるような正則行列Vを求めること(同時対角化)である。
このような行列Vは、一般化固有値問題に基づいて、以下のようにして求めることができる。(22)式に示す、空間共分散行列R1,R2に対する一般化固有値問題を解く。
すると、I個の一般化固有値λ1,・・・,λIと、λ1,・・・,λIにそれぞれ対応するI個の一般化固有ベクトルp1,・・・,pIと、が得られる。ここで、一般化固有ベクトルp1,・・・,pIは、次の(23)式を満たすように取ることができる。
これらの一般化固有ベクトルp1,・・・,pIを用いて、行列Pを次の(24)式により定めると、これは、次の(25)式及び(26)式を満たす(例えば、参考文献1「G. Strang, “Linear Algebra and Its Applications”, Harcourt Brace Jovanovich, 1988.」参照)。
ただし、Iは単位行列、Λは一般化固有値λ1,・・・,λIからなる(27)式の対角行列である。
したがって、行列Vを次の(28)式により定めると、これは同時無相関化の条件である(20)式及び(21)式を満たす(Λ1=Λ、Λ2=Iである)。
以下、(27)式における一般化固有値からなる対角行列Λを一般化固有値行列と呼び、(24)式における一般化固有ベクトルからなる行列Pを一般化固有ベクトル行列と呼ぶことがある。また、上述のように、Pは同時無相関化を実現する行列でもあるから、その意味でPを同時無相関化行列とも呼ぶ。
信号パラメータ分析部20は、空間共分散行列R1,R2に対する一般化固有値λ1,・・・,λIと、λ1,・・・,λIにそれぞれ対応する一般化固有ベクトルp1,・・・,pIであって(23)式を満たすものと、を計算し、一般化固有値λ1,・・・,λIを対角要素とする対角行列である一般化固有値行列Λと、一般化固有ベクトルp1,・・・,pIを列ベクトルとする正則行列である一般化固有ベクトル行列Pと、を、(27)式と(24)式とにより作成する。
上述のように、2個の源信号成分x1(k),x2(k)は、(24)式で定義される行列Pのエルミート転置PHを掛ける線型変換により、同時無相関化される。観測信号ベクトルのモデルを示す(9)式の両辺に対し、この線形変換を施すと、次の(29)式を得る。
ここで、(18)式、(19)式、(25)式、(26)式及び(28)式より、x´1(k),x´2(k)の共分散行列S´1(k),S´2(k)は、次の(30)式及び(31)式のように、ともに対角行列である。
ただし、Λは(27)式で定義される対角行列である。(29)式にマルチチャネルウィーナーフィルタを適用する場合、x´j(k)の推定値^x´j(k)は、次の(32)式で与えられる。
(32)式における共分散行列S´1(k),S´2(k)は、(30)式及び(31)式に示すように対角行列である。したがって、(32)式は、以下の(36)式及び(40)式に示すように非常に簡単になる。このように、同時無相関化により、(12)式における逆行列演算(v1(k)R1+v2(k)R2)−1は、(33)式と(37)式のように対角行列に対する逆行列演算(v1(k)Λ+v2(k)I)−1に変わり、単なる対角要素ごとの逆数演算により実現できるようになるため、少ない計算量で計算できるようになる。
シングルチャネルウィーナーフィルタ部40は、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を(36)式及び(40)式により計算する。
(36)式及び(40)式の処理は、シングルチャネルウィーナーフィルタに基づく処理と解釈することもできる。そのことを説明するために、まず、シングルチャネルウィーナーフィルタとは何であるかを説明する。マルチチャネルウィーナーフィルタが複数(I個)の観測信号を用いて推定を行うのに対し、シングルチャネルウィーナーフィルタは1個の観測信号のみを用いて推定を行う。シングルチャネルウィーナーフィルタでは、1個の観測信号(スカラー)y(k)を、2個の源信号成分(スカラー)x1(k),x2(k)の和として次の(41)式でモデル化する。
ここで、観測信号が1個しかないことに対応して、y(k),x1(k),x2(k)はいずれもスカラーである。源信号成分xj(k)は、平均が0、分散がσ2 j(k)であるとし、2個の源信号成分x1(k),x2(k)は互いに無相関であるとする。シングルチャネルウィーナーフィルタに基づく方法では、源信号成分xj(k)の推定値^xj(k)は、次の(42)式により計算される。
すなわち、この方法では、観測信号y(k)に対し、シングルチャネルウィーナーフィルタσ2 j(k)/(σ2 1(k)+σ2 2(k))を掛けることにより、源信号成分xj(k)の推定値^xj(k)を計算する。
ここで、(29)式を要素ごとに見た場合、第i要素より次の(43)式を得る。
上式に基づくシングルチャネルウィーナーフィルタを用いて、[y´(k)]iから[x´j(k)]iを推定することを考える。一般に、確率ベクトルαの共分散行列がAであるとき、αの第l要素[α]lの分散は、Aの(l,l)要素[A]llであることから、[x´j(k)]iの分散は[S´j(k)]iiである。よって、(30)式及び(31)式より、[x´1(k)]iの分散はv1(k)λiであり、[x´2(k)]iの分散はv2(k)である。したがって、シングルチャネルウィーナーフィルタに基づく[x´j(k)]iの推定値[^x´j(k)]iは、次の(44)式及び(45)式で与えられる。
これらはそれぞれ、(36)式及び(40)式の第i要素に他ならない。
さて、上記のようにして得られた同時無相関化された源信号成分x´j(k)=PHxj(k)の推定値^x´j(k)は、推定が上手く行っていれば、次の(46)式のように真値に近い値を取ると期待される。
上式の両辺に行列(PH)−1を掛けると、次の(47)式を得る。
そこで、次の(48)式により、源信号成分xj(k)の推定値^xj(k)を計算する。
このように、第1の実施形態によれば、源信号成分xj(k)の推定値^xj(k)の計算において、サンプル点ごとの逆行列演算を行う必要がないため、少ない計算量で信号分離を実現することができる。なお、上式は逆行列(PH)−1を含むが、これはkに依存しないため、サンプル点ごとではなく1回計算するだけで良いことに注意する。
従来の信号分析方法における源信号成分xj(k)の推定値を^xj conv(k)で表し、第1の実施形態における源信号成分xj(k)の推定値を^xj prop(k)で表す。^xj conv(k)と^xj prop(k)とは理論上一致することが示される。ここで、^xj conv(k)は、(49)式によって与えられる。
また、(5)式と、(33)式〜(40)式と、(48)式と、により、^xj prop(k)は次の(50)式及び(51)式に等しい。
以下では、j=1の場合について説明する。なお、j=2の場合については、j=1の場合と同様であるから、説明を省略する。(25)式及び(26)式によって、以下の(52)式及び(53)式が導き出される。
これらを、(49)式に代入すると、次の(54)式〜(57)式を得る。
上述のことから、従来の信号分析方法と第1の実施形態とで、源信号成分xj(k)の推定値^xj(k)の計算結果は理論上一致する。ただし、もちろんコンピュータによる計算誤差が生じることはある。
このように、本実施の形態1に係る信号分析装置1では、信号パラメータ分析部20は、2個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列R1,R2に対するpi HR2pm=0(i,mは1以上I以下の互いに異なる任意の整数)を満たす、それぞれ異なるI個(Iは2以上の整数)の一般化固有ベクトルp1,・・・,pIを列ベクトルとする行列である同時無相関化行列Pまたはそのエルミート転置PHを、I個の位置で取得される観測信号による観測信号ベクトルについて2個の源信号に対応する成分を無相関化するためのパラメータとして得る。信号分析装置1では、この無相関化するためのパラメータを用いて同時無相関化を行うことによって、フィルタリングでは、計算量の少ない要素ごとの掛け算を行うだけで足り、源信号成分の推定値の計算量を低減することができる。
同時無相関化部30は、信号パラメータ分析部20が得た同時無相関化行列Pのエルミート転置PHを、入力されたI個の観測信号による観測信号ベクトルに掛けることにより、2個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る。
そして、シングルチャネルウィーナーフィルタ部40は、同時無相関化部30が得た同時無相関化済観測信号ベクトルの各要素に、同時無相関化行列Pに基づく行列PHRjP(jは1以上2以下の整数)の対角成分に基づくシングルチャネルウィーナーフィルタを適用して同時無相関化済源信号成分推定値によるベクトルを得る。
さらに、同時無相関化逆変換部50は、シングルチャネルウィーナーフィルタ部40が得た同時無相関化済源信号成分推定値によるベクトルに、同時無相関化行列Pのエルミート転置PHの逆行列(PH)−1を掛けることにより、源信号成分の推定値を得る。
[第1の実施形態の効果]
このような構成を有する第1の実施形態によれば、上記のように、従来の信号分析方法と比べて、信号分離性能を低下させることなく、より少ない計算量で信号分離を実現することができる。すなわち、第1の実施形態では、観測信号ベクトルに対し、同時無相関化を行ってから、フィルタリングを行い、最後に同時無相関化の逆変換を行うことにより、源信号成分の推定値を計算する。ここで、同時無相関化を実現する行列Pは、空間共分散行列に対する一般化固有値問題を解くことにより、求めることができる。
この結果、第1の実施形態では、従来技術と比べて少ない計算量で源信号成分の推定値を計算できる。これは主に、以下の理由による。
第一に、第1の実施の形態では、同時無相関化により、フィルタリングは要素ごとの掛け算という計算量の小さい処理になり、サンプル点ごとの逆行列演算という計算量の大きい処理が不要になる。第二に、第1の実施の形態では、サンプル点によらない空間共分散行列に対する一般化固有値問題を解いて得られる行列Pを用いて同時無相関化を行うため、一般化固有値問題はサンプル点ごとではなく1回だけ解けば良い。これは、源信号成分の共分散行列が、サンプル点によらない空間共分散行列のスカラー倍であるとモデル化したことにより可能になったことである。第三に、第1の実施の形態では、同時無相関化の逆変換を行う際には、同時無相関化の行列PHの逆行列演算を行う必要があるが、同じ理由により、これはサンプル点ごとではなく1回だけ行えば良い。
さらに、第1の実施の形態では、観測信号ベクトルに対して、同時無相関化を行ってから、フィルタリングを行い、最後に同時無相関化の逆変換を行う、第1の実施形態における処理と、観測信号ベクトルに対して直接フィルタリングを行う従来技術における処理とは、等価である。このため、第1の実施形態によれば、従来技術と比べて信号分離性能を劣化させることなく、従来技術より少ない計算量で源信号成分の推定値を計算することができる。
[第1の実施形態の変形例1]
次に、第1の実施形態の変形例1について説明する。第1の実施形態では、一般化固有ベクトルp1,・・・,pIは、(23)式の制約条件を満たすものとした。本第1の実施形態の変形例1では、一般化固有ベクトルp1,・・・,pIが(23)式以外の制約条件を満たすものとした場合の信号分析装置について説明する。本第1の実施形態の変形例1に係る信号分析装置における信号パラメータ分析部では、空間共分散行列R1,R2に対する(58)式に示す一般化固有値問題を解く。
この結果、I個の一般化固有値λ1,・・・,λIとλ1,・・・,λIにそれぞれ対応するI個の一般化固有ベクトルp1,・・・,pIと、が得られる。ここで、上に述べた実施の形態1の例では、一般化固有ベクトルp1,・・・,pIは、次の(59)式を満たすように取る。
これに対して、第1の実施形態の変形例1では、一般化固有ベクトルp1,・・・,pIを、次の(60)式を満たすように取る(各σlは所定の正数)。
これらの一般化固有ベクトルp1,・・・,pIを用いて、行列Pを次の(61)式により定めると、これは、(62)式及び(63)式を満たす。
ただし、Λ1,Λ2は、次の(64)式及び(65)式に示す対角行列である。
したがって、行列Vを次の(66)式により定めると、これは、同時無相関化の条件(20)式及び(21)式を満たす。ただし、(64)式及び(65)式のΛ1及びΛ2が、(20)式及び(21)式におけるΛ1及びΛ2に相当する。
以下、(27)式における一般化固有値からなる対角行列Λを一般化固有値行列、(61)式における一般化固有ベクトルからなる行列Pを一般化固有ベクトル行列と呼ぶことがある。このようにして求めた行列Pを、同時無相関化行列として用いることができる。本第1の実施形態の変形例1に係る信号分析装置におけるシングルチャネルウィーナーフィルタ部では、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を、次の(67)式及び(68)式により計算する。
他の各部における処理は、上記の第1の実施形態と同様である。
このように、第1の実施形態の変形例1では、信号パラメータ分析部20は、2個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列R1,R2に対するpi HR2pm=0(i,mは1以上I以下の互いに異なる任意の整数)を満たす、それぞれ異なるI個(Iは2以上の整数)の一般化固有ベクトルp1,・・・,pIを列ベクトルとする行列である同時無相関化行列Pまたはそのエルミート転置PHを、I個の位置で取得される観測信号による観測信号ベクトルについて2個の源信号に対応する成分を無相関化するためのパラメータとして得る。
[第1の実施形態の変形例2]
次に、第1の実施形態の変形例2について説明する。第1の実施形態の変形例2では、信号パラメータ分析部における処理のもう一つの変形例について説明する。第1の実施形態では、信号パラメータ分析部20は、固有値分析に基づいて、PHR1PとPHR2Pとが対角行列になる行列Pとそのときの対角行列Λ=PHR1Pを求めていた。これに対し、本第1の実施形態の変形例2では、信号パラメータ分析部20は、最適化に基づいて、Λ1=PHR1Pと、Λ2=PHR2Pとができるだけ、ある対角行列Σ1,Σ2に近づくような行列PとそのときのΛ1=PHR1P及びΛ2=PHR2Pを求める。
すなわち、第1の実施形態の変形例2では、信号パラメータ分析部20は、次のような行列PとそのときのΛ1=PHR1P及びΛ2=PHR2Pを求める((69)式及び(70)式参照)。
具体的には、次の(71)式のコスト関数を行列Pと対角行列Σ1とΣ2に関して最小化することにより、行列Pを求めることができる。
ここで、||A||Fは、行列Aのフロベニウスノルムを表す。例えば、上式のコスト関数の最小化は、参考文献2「A. Yeredor, “Non-Orthogonal Joint Diagonalization in the Least-Squares Sense with Application in Blind Source Separation”, IEEE Transactions on Signal Processing, vol. 50, no. 7, pp. 1545−1553, Jul. 2002.」に記載のアルゴリズムに従って行うことができる。このようにして求めた行列Pを、同時無相関化行列として用いることができる。
また、Λ1=PHR1P及びΛ2=PHR2Pに基づき、第1の実施形態の変形例1と同様にして、シングルチャネルウィーナーフィルタ部40における同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を計算することができる。他の各部における処理は、上に述べた第1の実施形態と同様である。
このように、第1の実施形態の変形例2では、信号パラメータ分析部20は、J個(この変形例ではJは2)の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Rj(jは1以上J以下の整数)について、PHRjPがすべて対角行列に近づくように求めた行列である同時無相関化行列Pまたは/及びそのエルミート転置PHを、それぞれ異なるI個(Iは2以上の整数)の位置で取得される観測信号による観測信号ベクトルについてJ個の源信号に対応する成分を無相関化するためのパラメータとして得る。
[第1の実施形態の変形例3]
次に、第1の実施形態の変形例3について説明する。本第1の実施形態の変形例3では、第1の実施形態の変形例2の方法に基づき、源信号が2個の場合に限らず、源信号が一般にJ個の場合に拡張する場合について説明する。
なお、第1の実施形態の変形例3においては、J個の源信号が混在する状況において、それぞれ異なる位置で取得されたI個の観測信号yi(k)と、源信号のサンプル点ごとの分散(パワー)をモデル化するパラメータである分散パラメータvj(k)(j=1,・・・,Jは源信号のインデックス)と、源信号の空間的特性をモデル化するパラメータである空間共分散行列Rjと、が信号分析装置1に入力されるものとする。
第1の実施の形態の変形例3では、信号パラメータ分析部20は、最適化に基づいて、Λ1=PHR1P,・・・,ΛJ=PHRJPが、できるだけ、ある対角行列Σ1,・・・,ΣJに近づくような行列Pと、そのときのΛ1=PHR1P,・・・,ΛJ=PHRJPを求める。
すなわち、第1の実施形態の変形例3では、信号パラメータ分析部20は、次の(72)式〜(73)式に示すような行列Pと、そのときのΛ1=PHR1P,・・・,ΛJ=PHRJPを求める。
具体的には、次の(74)式のコスト関数を行列Pと各対角行列Σj(jは1以上J以下の整数)に関して最小化することにより、行列Pを求めることができる。
例えば、上式のコスト関数の最小化は、参考文献2に記載のアルゴリズムに従って行うことができる。本第1の実施形態の変形例3では、このようにして求めたPを、同時無相関化行列として用いることができる。また、Λ1=PHR1P,・・・,ΛJ=PHRJPに基づき、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を計算することができる。
次に、源信号がJ個の場合の信号分析装置1の構成と処理について説明する。信号分析装置1は、観測信号ベクトル作成部10、信号パラメータ分析部20、同時無相関化部30、シングルチャネルウィーナーフィルタ部40及び同時無相関化逆変換部50を有する。観測信号ベクトル作成部10は、入力された観測信号yi(k)を取得し、すべての観測信号からなるI次元縦ベクトルである観測信号ベクトルy(k)をサンプル点ごとに作成する。観測信号ベクトルy(k)は、次の(75)式で示される。
信号パラメータ分析部20は、入力された空間共分散行列Rjを取得する。
そして、上述の(76)〜(78)式で定義される各行列Λj(jは1以上J以下の整数)が対角行列になるような同時無相関化行列Pと、そのときのΛ1,・・・,ΛJを計算する。
同時無相関化部30は、観測信号ベクトル作成部10からの観測信号ベクトルy(k)と、信号パラメータ分析部20からの同時無相関化行列Pと、に基づいて、同時無相関化された観測信号ベクトルy´(k)を計算する。具体的には、(79)式に示すように、同時無相関化部30は、観測信号ベクトルy(k)に同時無相関化行列Pのエルミート転置PHを掛けることにより、同時無相関化された観測信号ベクトルy´(k)を計算する。
なお、同時無相関化とは、観測信号ベクトルy(k)に含まれる、J個の源信号に対応する成分(源信号成分)を同時に無相関化する線型変換を指す。
シングルチャネルウィーナーフィルタ部40は、入力された分散パラメータvj(k)を取得し、信号パラメータ分析部20からの行列Λjと、同時無相関化部30からの同時無相関化された観測信号ベクトルy´(k)と、に基づいて、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を計算する。具体的には、シングルチャネルウィーナーフィルタ部40は、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を次の(80)式により計算する。
同時無相関化逆変換部50は、信号パラメータ分析部20からの同時無相関化行列Pと、シングルチャネルウィーナーフィルタ部40からの同時無相関化された源信号成分x´j(k)の推定値^x´j(k)と、に基づいて、源信号成分xj(k)の推定値^xj(k)を計算し、出力する。具体的には、同時無相関化逆変換部50は、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)に対し、(81)式のように、同時無相関化部30で掛けた行列PHの逆行列(PH)−1を掛けることにより、源信号成分xj(k)の推定値^xj(k)を計算し、出力する。
このように、第1の実施形態の変形例3では、信号パラメータ分析部20は、J個(Jは2以上の整数)の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Rj(jは1以上J以下の整数)について、PHRjPがすべて対角行列に近づくように求めた行列である同時無相関化行列Pまたは/及びそのエルミート転置PHを、それぞれ異なるI個(Iは2以上の整数)の位置で取得される観測信号による観測信号ベクトルについてJ個の源信号に対応する成分を無相関化するためのパラメータとして得る。
[第1の実施形態及び第1の実施形態の変形例1〜3の効果]
このように、本第1の実施形態及び本第1の実施形態の変形例1〜3によれば、源信号成分xj(k)の推定値^xj(k)を少ない計算量で計算することができる。特に、第1の実施形態の固有値分析に基づく構成によれば、従来技術と比べて信号分離性能を劣化させることなく、従来技術より少ない計算量で信号分離を行うことができる。
なお、第1の実施形態の変形例2,3では、対角行列Σjに近づくように行列Λjを求める。ただし、行列Λjは実際には対角行列にならないことがある。この場合には、対角行列に近い行列Λj=行列PHRjP」の「対角成分」をλとして使ってシングルチャネルウィーナーフィルタの処理を行う。したがって、第1の実施形態の変形例2,3のシングルチャネルウィーナーフィルタの処理では、行列Λjの非対角成分は用いず、対角成分のみを用いるので、非対角成分は特定しない構成としても良い。
[第2の実施形態]
次に、第2の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。
第1の実施形態で説明したように、分散パラメータvj(k)及び空間共分散行列Rjが既知の場合には、これらに基づいて、源信号成分xj(k)の推定値^xj(k)を計算できる。一方、分散パラメータvj(k)及び空間共分散行列Rjが未知の場合には、これらを推定する必要がある。
そこで、本第2の実施形態においては、2個の源信号が混在する状況において、分散パラメータvj(k)及び空間共分散行列Rjと、源信号成分xj(k)と、を反復法により同時に推定する方法について説明する。第2の実施形態においては、それぞれ異なる位置で取得されたI個の観測信号yi(k)が信号分析装置に入力されるものとする。
[信号分析装置の構成及び処理]
図3及び図4を用いて、第2の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図3は、第2の実施形態に係る信号分析装置の構成の一例を示す図である。図4は、第2の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
図3に示すように、信号分析装置201は、初期化部(図示しない)、観測信号ベクトル作成部210、信号パラメータ分析部220、同時無相関化部230、シングルチャネルウィーナーフィルタ部240、事後共分散行列更新部250、分散パラメータ更新部260、空間共分散行列更新部270、収束判定部(図示しない)及び同時無相関化逆変換部280を有する。
まず、信号分析装置201の各部の概要について説明する。図示しない初期化部は、分散パラメータvj(k)の初期推定値^vj (0)(k)と、空間共分散行列Rjの初期推定値^Rj (0)と、を設定する(ステップS20)。これらは、例えば乱数に基づいて設定すればよい。
観測信号ベクトル作成部210は、入力された観測信号yi(k)を取得し(ステップS21)、すべての観測信号からなるI次元縦ベクトルである観測信号ベクトルy(k)をサンプル点ごとに作成する(ステップS22)。
信号パラメータ分析部220は、空間共分散行列更新部270からの同時無相関化された空間共分散行列R´jの推定値^R´j(ただし、例外として、信号パラメータ分析部220における最初の処理の際には、入力された空間共分散行列Rjの初期推定値^Rj (0))に基づいて、一般化固有値問題を解くことにより、一般化固有値行列Λ及び一般化固有ベクトル行列Pを更新する(ステップS23)。
同時無相関化部230は、観測信号ベクトル作成部210からの観測信号ベクトルy(k)と、信号パラメータ分析部220からの一般化固有ベクトル行列Pと、に基づいて、同時無相関化された観測信号ベクトルy´(k)を更新する(ステップS24)。
シングルチャネルウィーナーフィルタ部240は、信号パラメータ分析部220からの一般化固有値行列Λと、同時無相関化部230からの同時無相関化された観測信号ベクトルy´(k)と、分散パラメータ更新部260からの分散パラメータvj(k)の推定値^vj(k)(ただし、例外として、シングルチャネルウィーナーフィルタ部240における最初の処理の際には、入力された分散パラメータvj(k)の初期推定値^vj (0)(k)と、に基づいて、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を更新する(ステップS25)。
事後共分散行列更新部250は、信号パラメータ分析部220からの一般化固有値行列Λと、分散パラメータ更新部260からの分散パラメータvj(k)の推定値^vj(k)(ただし、例外として、事後共分散行列更新部250における最初の処理の際には、入力された分散パラメータvj(k)の初期推定値^vj (0)(k))と、に基づいて、同時無相関化された事後共分散行列Σ´(k)を更新する(ステップS26)。
分散パラメータ更新部260は、信号パラメータ分析部220からの一般化固有値行列Λと、シングルチャネルウィーナーフィルタ部240からの同時無相関化された源信号成分x´j(k)の推定値^x´j(k)と、事後共分散行列更新部250からの同時無相関化された事後共分散行列Σ´(k)と、に基づいて、分散パラメータvj(k)の推定値^vj(k)を更新する(ステップS27)。
空間共分散行列更新部270は、シングルチャネルウィーナーフィルタ部240からの同時無相関化された源信号成分x´j(k)の推定値^x´j(k)と、事後共分散行列更新部250からの同時無相関化された事後共分散行列Σ´(k)と、分散パラメータ更新部260からの分散パラメータvj(k)の推定値^vj(k)と、に基づいて、同時無相関化された空間共分散行列R´jの推定値^R´jを更新する(ステップS28)。
図示しない収束判定部は、収束したかどうかの判定を行う(ステップS29)。収束判定部が収束していないと判定した場合(ステップS29:No)、ステップS23の信号パラメータ分析部220での処理に戻って、処理が継続される。
一方、収束判定部が収束したと判定した場合(ステップS29:Yes)、同時無相関化逆変換部280での以下の処理が行われる。
具体的には、同時無相関化逆変換部280は、信号パラメータ分析部220からの一般化固有ベクトル行列Pと、シングルチャネルウィーナーフィルタ部240からの同時無相関化された源信号成分x´j(k)の推定値^x´j(k)と、に基づいて、源信号成分xj(k)の推定値^xj(k)を計算し、出力する(ステップS30)。
次に、信号分析装置201の各部の詳細について説明する。なお、初期化部については上述の通りである。
観測信号ベクトル作成部210は、入力された観測信号yi(k)を取得し、すべての観測信号からなるI次元縦ベクトルである観測信号ベクトルy(k)を次の(82)式により作成する。
信号パラメータ分析部220は、信号パラメータ分析部220における最初の処理の際に、初期化部からの空間共分散行列Rjの初期推定値^Rj (0)を受け取る。そして、信号パラメータ分析部220は、^R1 (0),^R2 (0)に対する一般化固有値問題を解くことにより、一般化固有値λ1,・・・,λIと、λ1,・・・,λIにそれぞれ対応する一般化固有ベクトルp1,・・・,pIであって、(83)式を満たすものを計算する。
これによって、信号パラメータ分析部220は、一般化固有値行列Λ及び一般化固有ベクトル行列Pを次の(84)式及び(85)式により更新する。
信号パラメータ分析部220は、信号パラメータ分析部220における2回目以降の処理の際に、空間共分散行列更新部270からの同時無相関化された空間共分散行列R´jの推定値^R´jを受け取って、^R´1,^R´2に対する一般化固有値問題を解くことにより、一般化固有値λ1,・・・,λIと、λ1,・・・,λIにそれぞれ対応する一般化固有ベクトルq1,・・・,qIであって、(86)式を満たすものを計算する。
そして、信号パラメータ分析部220は、一般化固有値行列Λ及び行列Qを(87)式及び(88)式により更新する。
そして、信号パラメータ分析部220は、一般化固有ベクトル行列Pを次の(89)式により更新する。
同時無相関化部230は、(90)式に示すように、観測信号ベクトルy(k)に一般化固有ベクトル行列Pのエルミート転置PHを掛けることにより、同時無相関化された観測信号ベクトルy´(k)を更新する。
シングルチャネルウィーナーフィルタ部240は、同時無相関化された源信号成分x´j(k)の推定値^x´j(k)を(91)式及び(92)式により更新する。
ただし、例外として、シングルチャネルウィーナーフィルタ部240における最初の処理の際に、上の2式において、分散パラメータ更新部260からの分散パラメータvj(k)の推定値^vj(k)の代わりに、初期化部からの分散パラメータvj(k)の初期推定値^vj (0)(k)を用いる。
事後共分散行列更新部250は、同時無相関化された事後共分散行列Σ´(k)を、次の(93)式により更新する。
ただし、例外として、事後共分散行列更新部250における最初の処理の際には、上式において、分散パラメータ更新部260からの分散パラメータvj(k)の推定値^vj(k)の代わりに、初期化部からの分散パラメータvj(k)の初期推定値^vj (0)(k)を用いる。
分散パラメータ更新部260は、分散パラメータvj(k)の推定値^vj(k)を次の(94)式により更新する。
空間共分散行列更新部270は、同時無相関化された空間共分散行列R´jの推定値^R´jを次の(95)式により更新する。(95)式において、Kはサンプル点の数である。
同時無相関化逆変換部280は、源信号成分xj(k)の推定値^xj(k)を次の(96)式により計算し、出力する。
[第2の実施形態に係る信号分析装置の処理の背景]
[従来の信号分析方法]
次に、信号分析装置201の処理の背景について説明する。まず、従来技術について説明する。2個の源信号が混在する状況において、分散パラメータvj(k)及び空間共分散行列Rjと、源信号成分xj(k)と、を反復法により同時に推定する方法が、非特許文献1に開示されている。
ここで、従来技術における観測信号ベクトルy(k)のモデル化について説明する。(9)式における源信号成分xj(k)は、次の(97)式のように平均0、共分散行列vj(k)Rjの複素ガウス分布に従うものとする。
ここで、N(z;μ,Σ)は、平均がμ、共分散行列がΣの複素ガウス分布((98)式)を表す。
また、2個の源信号成分x1(k),x2(k)は互いに無相関であると仮定する。観測信号ベクトルy(k)が与えられた条件下での源信号成分xj(k)の事後確率密度は、次の(99)式で与えられる(例えば、非特許文献1を参照)。
ここで、平均^xj(k)及び共分散行列Σj(k)は、次の(100)式及び(101)式により与えられる。
(100)式における平均^xj(k)は、(10)式におけるマルチチャネルウィーナーフィルタに基づく源信号成分xj(k)の推定値^xj(k)に他ならない(両者に同じ記号を用いているのはそのためである)。したがって、(10)式におけるマルチチャネルウィーナーフィルタに基づく源信号成分xj(k)の推定値^xj(k)は、観測信号ベクトルy(k)が与えられた条件下での源信号成分xj(k)の事後確率密度を最大化する点と解釈できることが分かった。(101)式はjに依らないから、(102)式と書く。
従来技術では、下記第1の処理〜第4の処理により、分散パラメータvj(k)及び空間共分散行列Rjと、源信号成分xj(k)と、を同時に推定する。
第1の処理では、分散パラメータvj(k)の推定値^vj(k)と、空間共分散行列Rjの推定値^Rjと、を初期化する。これらは、例えば、乱数に基づいて初期化すればよい。
第2の処理では、分散パラメータvj(k)の推定値^vj(k)と、空間共分散行列Rjの推定値^Rjと、に基づいて、源信号成分xj(k)を推定する。第1の実施形態では、源信号成分xj(k)の1つの推定値^xj(k)を計算したのに対し、ここでは源信号成分xj(k)の事後確率密度((99)式参照)を更新する。具体的には、(99)式の平均^xj(k)と共分散行列Σ(k)((102)式の通りΣ1(k)及びΣ2(k)に等しい。(101)式も参照。)を次の(103)式及び(104)式により更新する。
第3の処理では、源信号成分xj(k)の事後確率密度の平均^xj(k)と共分散行列Σ(k)に基づいて、分散パラメータvj(k)の推定値^vj(k)と、空間共分散行列Rjの推定値^Rjと、を次の(105)式及び(106)式により更新する。ここで、trはトレースである。
第4の処理では、第2及び第3の処理を、収束するまで交互に反復し、収束した場合に、源信号成分xj(k)の推定値^xj(k)を出力する。
従来技術の問題点として、膨大な計算量を要するため、サンプル点kの数が膨大になる大規模なデータセットの処理、短時間での処理が要求されるリアルタイム処理、或いは、計算能力の低い機器などへの適用が困難だった。
実際、従来技術では、源信号成分xj(k)の推定値^xj(k)と事後共分散行列Σ(k)を計算する際に、逆行列計算(具体的には(v1(k)R1+v2(k)R2)−1の計算)と行列乗算(具体的には(v1(k)R1)(v1(k)R1+v2(k)R2)−1(v2(k)R2)の計算)とを、サンプル点毎に(つまり多数回)行う必要があった。逆行列計算と行列乗算とはともに、I次行列に対し、O(I3)という大きな計算量を要する演算である。また、サンプル点の数は、典型的には数百〜数千万程にもなる。
このため、従来の信号分析方法では、源信号成分xj(k)の推定値^xj(k)と事後共分散行列Σ(k)を計算する際に、膨大な計算量を必要としていた。さらに、従来では、信号分析装置全体として、上記の第4の処理における反復を行う必要がある(反復回数は典型的には数十回にもなる)ため、計算量はさらに膨大となる。
[第2の実施の形態に係る信号分析方法]
これに対し、第2の実施形態では、2個の源信号成分x1(k),x2(k)の同時無相関化に基づいて、源信号成分xj(k)の推定値^xj(k)と事後共分散行列Σ(k)と空間共分散行列Rと分散パラメータvとを少ない計算量で計算する方法を提供する。この結果、信号分析装置201全体としても、少ない計算量で信号分析を実現することができる。以下、第2の実施形態における考え方について説明する。
従来技術(非特許文献1参照)の更新則(103)式〜(106)式において、何回目の反復であるかを明示すると、以下の(107)式〜(110)式のようになる。なお、n=1,・・・,Nは何回目の反復であるかを表すインデックス、Nは反復回数である。
空間共分散行列R1,R2の推定値^R1 (n−1),^R2 (n−1)の一般化固有値問題に対して、第1の実施形態と同様にして定義される一般化固有値行列及び一般化固有ベクトル行列を、Λ(n−1)及びP(n−1)で表す。これらは、(111)式及び(112)式を満たす。
すなわち、Λ(n−1)及びP(n−1)は(113)式及び(114)式を満たす。P(n−1)は同時無相関化(同時対角化)を実現する行列であるから、その意味でP(n−1)を同時無相関化行列とも呼ぶ。
(107)式からサンプル点ごとの逆行列演算を消すために、(107)式に(113)式と(114)式とを代入して整理すると、第1の実施形態と同様に、次の(115)式、(116)式を得る。
同様に、(108)式からサンプル点ごとの逆行列演算を消すために、(108)式に、(113)式と(114)式とを代入すると、次の(118)式〜(121)式を得る。
これで、サンプル点ごとの逆行列演算を消すことができた。しかしながら、(121)式は、まだ、サンプル点ごとの行列の乗算を含んでいる。そこで、Σ(n)(k)の代わりに、(122)式により定義される、Σ´(n)(k)を更新することを考える。
(121)式より、Σ´(n)(k)は、(123)式で表せる。
(123)式には、サンプル点ごとの行列の乗算が含まれない。Σ´(n)(k)は同時無相関化の線型変換を施された事後共分散行列とみなせるから、以下、同時無相関化された事後共分散行列と呼ぶ。
そして、(115)式と(116)式とにおいても、^xj (n)(k)の代わりに、(124)式で定義される^x´j (n)(k)を更新することを考える。
(115)式及び(116)式より、(125)式及び(126)式となる。
^x´
j (n)(k)は同時無相関化された源信号成分である。
次に、(109)式に着目すると、これは、^xj (n)(k)とΣ(n)(k)に依存している。上述のように、第2の実施形態では、^xj (n)(k)とΣ(n)(k)の代わりに^x´j (n)(k)とΣ´(n)(k)とを更新する。そこで、(109)式も、^x´j (n)(k)とΣ´(n)(k)とを用いて計算できる形に書き換える必要がある。(124)式及び(122)式より、^xj (n)(k)とΣ(n)(k)は(127)式及び(128)式と表せる。
これらを(109)式に代入すると(129)式〜(133)式を経て(134)式が得られる。
(110)式についても同様に、^x´j (n)(k)とΣ´(n)(k)とを用いて計算できる形に書き換えるために、(127)式と(128)式とを代入すると、次の(135)式及び(136)式のようになる。
ここでも、^Rj (n)の代わりに、(137)式で定義される^R´j (n)を更新することにする。
これにより、逆行列演算と行列の乗算とをいくつか減らせるため、計算量を削減することができる。また、^R´j (n)は、(136)式より、(138)式のように表せる。
^R´j (n)は同時無相関化の線型変換を施された空間共分散行列とみなせるから、以下、同時無相関化された空間共分散行列と呼ぶ。
最後に、一般化固有値行列Λ及び一般化固有ベクトル行列Pの更新について説明する。最初の反復においては、空間共分散行列R1,R2の初期推定値^R1 (0),^R2 (0)の一般化固有値問題を解くことにより、Λ(0)及びP(0)を求めればよい。
2回目以降の反復においては、同時無相関化された空間共分散行列の推定値^R´1 (n),^R´2 (n)が得られている。これを同時無相関化の逆変換により^R1 (n),^R2 (n)に戻してから、一般化固有値問題を適用することも可能だが、やや無駄な計算が生じる。
そこで、以下のようにすれば、より効率的に計算できる。まず、^R´1 (n),^R´2 (n)の一般化固有値問題を解いて、一般化固有値を対角成分とする対角行列Λ(n)と、一般化固有ベクトルを列ベクトルとする行列Q(n)と、を求める。これらは、(Q(n))H(^R´1 (n))Q(n)=Λ(n)及び(Q(n))H(^R´2 (n))Q(n)=Iを満たすものとする。
次に、P(n)を(139)式により更新する。
このとき、Λ(n)及びP(n)は、^R1 (n),^R2 (n)に対する一般化固有値行列及び一般化固有ベクトル行列になる。これは、次の(140)式〜(147)式のようにして確認できる。
[第2の実施形態の効果]
このように、第2の実施形態によれば、(同時無相関化された)源信号成分の推定値及び(同時無相関化された)事後共分散行列と空間共分散行列と分散パラメータとを従来技術と比べて少ない計算量で計算できるため、従来技術と比べて少ない計算量で信号分離を実現できる。
これは以下の理由による。第一に、第1の実施形態と同様に、同時無相関化に基づくことで、サンプル点ごとの逆行列演算が不要になる。第二に、事後共分散行列Σ(n)(k)そのものではなく、同時無相関化された事後共分散行列Σ´(n)(k)を更新するようにすることで、サンプル点ごとの行列の乗算が不要になる。
また、上記のように、従来の信号分析方法の処理と第2の実施形態の処理とは理論上等価であるため、従来の信号分析方法と第2の実施形態とで信号分析結果は理論上一致する(ただし、もちろんコンピュータによる計算誤差が生じることはある)。したがって、第2の実施形態によれば、従来の信号分析方法と比べて、信号分析性能を低下させることなく、より少ない計算量で信号分析を実現することができる。
また、第1の実施形態とは異なり、第2の実施形態によれば、分散パラメータvj(k)及び空間共分散行列Rjが未知の場合でも、分散パラメータvj(k)及び空間共分散行列Rjと、源信号成分xj(k)と、を同時に推定することができる。
[第3の実施形態]
次に、第3の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。
第3の実施形態においては、2個の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得されたI個(Iは2以上の整数)の観測信号yi(t)が信号分析装置に入力されるものとする。観測信号yi(t)は、周波数分解され、周波数ごとに第2の実施形態に記載の方法で信号分離がなされ、最後に時間領域に変換される。
[信号分析装置の構成及び処理]
図5及び図6を用いて、第3の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図5は、第3の実施形態に係る信号分析装置の構成の一例を示す図である。図6は、第3の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
図5に示すように、信号分析装置301は、周波数領域変換部310、周波数領域信号分離部320及び時間領域変換部330を有する。
まず、信号分析装置301の各部の概要について説明する。周波数領域変換部310は、入力された観測信号yi(t)を受け取り(ステップS40)、短時間フーリエ変換などにより時間周波数領域の観測信号yi(k,f)を計算する(ステップS41)。ここで、kは時間フレームのインデックスであり、fは周波数ビンのインデックスである。以下、yi(k,f)をyi f(k)で表す。
周波数領域信号分離部320は、fごとに第2の実施形態の方法を適用して、周波数領域での信号分離を行う(ステップS42)。すなわち、f=1に対しては、yi 1(k)を入力として、第2の実施形態の方法により源信号成分の推定値^x1 1(k),^x2 1(k)を得る。
そして、f=2に対しては、yi 2(k)を入力として、第2の実施形態の方法により、源信号成分の推定値^x1 2(k),^x2 2(k)を得る。各fについてこれらを繰り返す。
そして、f=Fに対しては、yi F(k)を入力として、第2の実施形態の方法により、源信号成分の推定値^x1 F(k),^x2 F(k)を得る。以下、^xj f(k)を、^xj(k,f)で表す。
時間領域変換部330は、時間周波数領域の源信号成分の推定値^xj(k,f)に基づいて、逆短時間フーリエ変換などにより、時間領域の源信号成分の推定値^xj(t)を計算し、出力する(ステップS43)。
このように、第3の実施形態では、時間周波数領域において、周波数ビンfごとに第2の実施形態に係る方法を適用する。これに対し、非特許文献1には、時間周波数領域において、周波数ビンfごとに第2の実施形態に記載した従来技術を適用する方法が開示されている。
この従来技術の問題点として、膨大な計算量を要するため、大規模なデータセットの処理や、リアルタイム処理(短時間での処理が要求される)、計算能力の低い機器などへの適用が困難だった。実際、従来技術では、逆行列計算と行列乗算とを、周波数ビンごと、フレームごと、反復ごとに(つまり多数回)行う必要があった。このため、従来の信号分析方法では、膨大な計算量を必要としていた。
例えば、後述の実験では、8秒の観測信号に対して、周波数ビンの数は256、フレームの数は249、反復数は10であり、処理に750秒を要した。これが80000秒(およそ22時間)のデータセットになれば、処理に7500000秒(およそ87日)掛かる計算になる。
[第3の実施形態の効果]
本発明の効果を確認するために、従来の方法及び第3の実施形態を用いた音源分離実験について説明する。
図7は、従来の方法及び第3の実施形態を用いた音源分離実験の実験条件を示す図である。図8は、従来の方法及び第3の実施形態を用いた音源分離実験結果(real time factor)を示す図である。図9は、従来の方法及び第3の実施形態を用いた音源分離実験結果(音源分離性能)を示す図である。
8秒間の英語音声に、図7に示す実験室で測定された残響インパルス応答を畳み込んで残響音声を生成し、複数人の残響音声を加算することで観測信号を生成した。このとき、観測信号に対して、各方法を用いて音源分離を行った場合の性能を図8及び図9に示す。図8は、real time factor(観測信号長を1としたときの、処理に要した時間)である。また、図9は、音源分離性能である信号対歪み比(値が大きいほど音源分離性能が高い)である。
図8及び図9に示すように、本第3の実施形態に係る方法によれば、従来の方法と比べて音源分離性能を低下させることなく、250倍以上の計算時間削減が実現できたことが分かる。なお、残響時間が130ミリ秒及び440ミリ秒の場合には、本第3の実施形態に係る方法による信号対歪み比は、従来の方法と比べてわずかに低下しているが、これはコンピュータによる計算の誤差によるものである。このように、本第3の実施形態に係る方法は従来の方法と理論的には等価だが、実際には計算誤差に起因して性能差が生じることもあるものの、それは図8のようにごく僅かな差に過ぎない。
このように、第3の実施形態によれば、従来の方法と比べて音源分離性能を低下させることなく、計算時間の大きな削減が実現できる。
[第4の実施形態]
次に、第4の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。
第4の実施形態は、第2の実施形態において、分散パラメータvj(k)及び空間共分散行列Rjを、補助関数法を用いて推定する構成に変更したものである。
[信号分析装置の構成及び処理]
図10及び図11を用いて、第4の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図10は、第4の実施形態に係る信号分析装置の構成の一例を示す図である。図11は、第4の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
図10に示すように、信号分析装置401は、初期化部(図示しない)、観測信号ベクトル作成部410、信号パラメータ分析部420、同時無相関化部430、分散パラメータ更新部440、空間共分散行列更新部450、収束判定部(図示しない)を有する。
まず、信号分析装置401の各部の概要について説明する。図示しない初期化部は、分散パラメータvj(k)の初期推定値^vj (0)(k)と、空間共分散行列Rjの初期推定値^Rj (0)と、を設定する(ステップS51)。これらは、例えば乱数に基づいて設定すればよい。
観測信号ベクトル作成部410は、入力された観測信号yi(k)を取得し(ステップS52)、すべての観測信号からなるI次元縦ベクトルである観測信号ベクトルy(k)をサンプル点ごとに作成する(ステップS53)。
信号パラメータ分析部420は、空間共分散行列更新部450からの空間共分散行列Rjの推定値^Rjに基づいて、一般化固有値問題を解くことにより、一般化固有値行列Λ及び一般化固有ベクトル行列Pを更新する(ステップS54)。具体的には、信号パラメータ分析部420は、空間共分散行列R1,R2の推定値^R1,^R2の一般化固有値問題に対して、第1の実施形態と同様にして定義される一般化固有値行列Λ及び一般化固有ベクトル行列Pを求める。これらは、(148)式及び(149)式を満たすものとする。
ただし、例外として、信号パラメータ分析部420における最初の処理の際には、信号パラメータ分析部420における上述の処理において、空間共分散行列更新部450からの空間共分散行列Rjの推定値^Rjの代わりに、入力された空間共分散行列Rjの初期推定値^Rj (0)を用いる。
同時無相関化部430は、観測信号ベクトル作成部410からの観測信号ベクトルy(k)と、信号パラメータ分析部420からの一般化固有ベクトル行列Pと、に基づいて、同時無相関化された観測信号ベクトルy´(k)を(90)式により更新する(ステップS55)。
分散パラメータ更新部440は、同時無相関化部430からの同時無相関化された観測信号ベクトルy´(k)と、信号パラメータ分析部420からの一般化固有値行列Λと、に基づいて、分散パラメータvj(k)の推定値^vj(k)を次の(150)式及び(151)式により更新する(ステップS56)。
ここで、右辺の^v1(k)及び^v2(k)は、分散パラメータ更新部440における前回の処理で得られた分散パラメータの推定値(ただし、例外として、分散パラメータ更新部440における最初の処理の際には、初期化部からの分散パラメータの初期推定値)を表す。
空間共分散行列更新部450は、同時無相関化部430からの同時無相関化された観測信号ベクトルy´(k)と、信号パラメータ分析部420からの一般化固有ベクトル行列P及び一般化固有値行列Λと、分散パラメータ更新部440からの分散パラメータvj(k)の推定値^vj(k)と、に基づいて、空間共分散行列Rjの推定値^Rjを更新する(ステップS57)。具体的には、空間共分散行列更新部450は、まず次の(152)式により行列Cjを更新する。
次に、空間共分散行列更新部450は、次の(153)式によりベクトルz(k)を更新する。
次に、空間共分散行列更新部450は、次の(154)式により行列Djを更新する。
次に、空間共分散行列更新部450は、次の(155)式により空間共分散行列Rjの推定値^Rjを更新する。
図示しない収束判定部は、収束したかどうかの判定を行う(ステップS58)。収束判定部が収束していないと判定した場合(ステップS58:No)、信号パラメータ分析部420での処理(ステップS54)に戻って、処理が継続される。
一方、収束判定部が収束したと判定した場合(ステップS58:Yes)、分散パラメータ更新部440が分散パラメータvj(k)の推定値^vj(k)を出力し、空間共分散行列更新部450が空間共分散行列Rjの推定値^Rjを出力する(ステップS59)。
[第4の実施形態の背景]
第4の実施形態の背景について説明する。補助関数法に基づけば、分散パラメータvj(k)及び空間共分散行列Rjは、次のように推定できる。分散パラメータvj(k)及び空間共分散行列Rjの推定値はあらかじめ初期化され、以下に述べる反復処理により更新される。分散パラメータvj(k)の推定値^vj(k)は、次の(156)式により更新される。
空間共分散行列Rjの更新においては、まず(157)式及び(158)式により行列Cj及び行列Djが更新され、次に、(159)式により空間共分散行列Rjの推定値^Rjが更新される。
ここで、半正定値エルミート行列Aの平方根A1/2は、(160)式をAの固有値分解とすると、(161)式により定義される。ここで、UはAの固有ベクトルからなるユニタリ行列、ΛはAの固有値からなる対角行列、Σ1/2はΣの各対角要素[Σ]iiをその平方根([Σ]ii)1/2で置き換えることにより得られる対角行列である。
上記の分散パラメータvj(k)の推定値^vj(k)及び空間共分散行列Rjの推定値^Rjの更新は、収束するまで反復される。
上記の分散パラメータvj(k)の推定値^vj(k)及び空間共分散行列Rjの推定値^Rjの更新則は、参考文献3「Kazuyoshi Yoshii,Ryota Tomioka,Daichi Mochihashi,and Masataka Goto,“Infnite Positive Semidefinite Tensor Factorization for Source Separation of Mixture Signals”,Proceedings of the 30th International Conference on Machine Learning,2013.」、参考文献4「Hiroshi Sawada,Hirokazu Kameoka,Shoko Araki,and Naonori Ueda,“Multichannel Extensions of Non−Negative Matrix Factorization With Complex−Valued Data”,IEEE TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING,VOL.21,No.5,MAY 2013.」に開示されている方法と同様の方法により導出することができる。
上記の方法の問題点として、(156)式〜(158)式の各処理において、サンプル点ごとの逆行列演算を行う必要があるため、膨大な計算量を必要とする。
これに対し、本第4の実施形態では、同時無相関化(同時対角化)に基づいて、サンプル点ごとの逆行列演算を回避することにより、計算量を削減することができる。実際、(148)式及び(149)式を^R1、^R2について解いたものを(156)式〜(158)式に代入することにより、分散パラメータvj(k)及び空間共分散行列Rjの推定値の更新則が以下のように得られる。
分散パラメータvj(k)の推定値^vj(k)の更新則は、次の(162)式及び(163)式のようになる。
ここで、y´(k)は、(164)式の同時無相関化済観測信号ベクトルであり、Λ1は一般化固有値行列Λであり、Λ2は単位行列Iである。
空間共分散行列Rjの推定値^Rjの更新則は次の(165)式及び(166)式のようになる。まず、行列Cjに関しては、次式のようになる。
行列Djに関しては次の(167)式及び(168)式のようになる。
ただし、ベクトルz(k)は、次の(169)式〜(171)式により定義される。
[第4の実施形態の効果]
このように、同時無相関化(同時対角化)により、サンプル点ごとの逆行列演算を不要にし、計算量を削減することができる。具体的には、同時無相関化(同時対角化)に基づき、分散パラメータvj(k)の推定値^vj(k)の更新式を(162)式のように対角行列Λj=PHRjPを用いて表すことができ、その結果、(163)式のように逆行列演算を不要にすることができる。空間共分散行列Rjの推定値^Rjの更新則についても同様である。
なお、(155)式の処理では、(159)式の処理と同様、逆行列演算を行う必要があるが、これはサンプル点ごとではなく各反復にて1回だけ行えば良いため、その計算量は比較的小さい。同様に、(155)式の処理では、(159)式の処理と同様、行列の平方根の計算のために固有値分解を行う必要があるが、これはサンプル点ごとではなく各反復にて2回だけ行えば良いため、その計算量は比較的小さい。また、(166)式、(168)式、及び(155)式の処理では、逆行列演算と同様に計算量の大きい演算である行列の乗算を含んでいるが、これもサンプル点ごとに行う必要はないため、その計算量は比較的小さい。
[第4の実施形態の変形例]
なお、本実施形態においては、源信号が2個の場合について説明したが、第1の実施形態の変形例3と同様の考え方に基づけば、源信号が2個の場合に限らず、源信号が一般にJ個の場合に拡張することができる。
上記のように、第1〜第4の実施形態では、信号パラメータ分析部が、J個(Jは2以上の整数)の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Rj(jは1以上J以下の整数)について、PHRjPがすべて対角行列になるような行列である同時無相関化行列Pまたは/及び、そのエルミート転置PHを、それぞれ異なるI個(Iは2以上の整数)の位置で取得される観測信号による観測信号ベクトルについてJ個の源信号に対応する成分を無相関化するためのパラメータとして得る。信号分析装置1,201,301,401では、この無相関化するためのパラメータを用いて同時無相関化を行うことによって、フィルタリングでは、計算量の少ない要素ごとの掛け算を行うだけで足り、源信号成分の推定値の計算量を低減することができる。
[第5の実施形態]
次に、第5の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。第5の実施形態は、第1の実施形態の変形例3のバリエーションである。第1の実施形態の変形例3では、分散パラメータvj(k)と空間共分散行列Rjが信号分析装置に入力されるとの前提条件がある。これに対し、第5の実施形態では、分散パラメータvj(k)と空間共分散行列Rjとが信号分析装置に入力される必要はなく、分散パラメータvj(k)と、空間共分散行列Rjをモデル化するパラメータである同時無相関化行列P及び同時無相関化された空間共分散行列Λjと、源信号成分xj(k)との同時推定を実現する。
本第5の実施形態では、一般にJ個の源信号が混在し、分散パラメータvj(k)及び空間共分散行列Rjが未知の状況で、分散パラメータvj(k)と、空間共分散行列Rjをモデル化するパラメータである同時無相関化行列P及び同時無相関化された空間共分散行列Λjと、源信号成分xj(k)と、を同時に推定することにより、信号分離を実現する方法について説明する。本第5の実施形態では、それぞれ異なる位置で取得されたI個の観測信号yi(k)が信号分析装置に入力されるものとする。
[信号分析装置の構成及び処理]
図12及び図13を用いて、第5の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図12は、第5の実施形態に係る信号分析装置の構成の一例を示す図である。図13は、第5の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
図12に示すように、信号分析装置501は、初期化部(図示しない)、観測信号ベクトル作成部510、同時無相関化部530、シングルチャネルウィーナーフィルタ部540、事後共分散行列更新部550、分散パラメータ更新部560、空間共分散行列更新部570、信号パラメータ分析部520、収束判定部(図示しない)及び同時無相関化逆変換部580を有する。
まず、信号分析装置501の各部の概要について説明する。図示しない初期化部は、分散パラメータvj(k)の初期推定値^vj (0)(k)と、同時無相関化行列Pの初期推定値^P(0)と、同時無相関化された空間共分散行列Λjの初期推定値^Λj (0)と、を設定する(ステップS61)。これらは、例えば乱数に基づいて設定すればよい。
観測信号ベクトル作成部510は、入力されたそれぞれ異なる位置で取得されたI個の観測信号yi(k)を取得し(ステップS62)、I次元縦ベクトルである観測信号ベクトルy(k)をサンプル点ごとに(82)式により作成する(ステップS63)。
同時無相関化部530は、観測信号ベクトル作成部510からの観測信号ベクトルy(k)と、信号パラメータ分析部520からの同時無相関化行列Pの推定値^P(ただし例外として、同時無相関化部530における最初の処理の際には、初期化部からの同時無相関化行列Pの初期推定値^P(0))と、に基づいて、同時無相関化された観測信号ベクトルy´(k)を次の(172)式により更新する(ステップS64)。
ただし、同時無相関化部530における最初の処理の際には、(172)式の右辺において、^Pを^P(0)で置き換える。
シングルチャネルウィーナーフィルタ部540は、同時無相関化部530からの同時無相関化された観測信号ベクトルy´(k)と、分散パラメータ更新部560からの分散パラメータvj(k)の推定値^vj(k)(ただし、例外として、シングルチャネルウィーナーフィルタ部540における最初の処理の際には、初期化部からの分散パラメータvj(k)の初期推定値^vj (0)(k))と、空間共分散行列更新部570からの同時無相関化された空間共分散行列Λjの推定値^Λj(ただし、例外として、シングルチャネルウィーナーフィルタ部540における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λjの初期推定値^Λj (0))と、に基づいて、同時無相関化された源信号成分xj´(k)の推定値^xj´(k)を次の(173)式により更新する(ステップS65)。
ただし、シングルチャネルウィーナーフィルタ部540における最初の処理の際には、(173)式の右辺において、^Λj(j=1,・・・,J)を^Λj (0)(j=1,・・・,J)に、^vj(k)(j=1,・・・,J;k=1,・・・,K)を^vj (0)(k)(j=1,・・・,J;k=1,・・・,K)に、それぞれ置き換える。
事後共分散行列更新部550は、分散パラメータ更新部560からの分散パラメータvj(k)の推定値^vj(k)(ただし、例外として、事後共分散行列更新部550における最初の処理の際には、初期化部からの分散パラメータvj(k)の初期推定値^vj (0)(k))と、空間共分散行列更新部570からの同時無相関化された空間共分散行列Λjの推定値^Λj(ただし、例外として、事後共分散行列更新部550における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λjの初期推定値^Λj (0))と、に基づいて、同時無相関化された事後共分散行列Σj´(k)を次の(174)式により更新する(ステップS66)。
ただし、事後共分散行列更新部550における最初の処理の際には、(174)式の右辺において、^Λjを^Λj (0)に、^vj(k)を^vj (0)(k)に、それぞれ置き換える。
分散パラメータ更新部560は、シングルチャネルウィーナーフィルタ部540からの同時無相関化された源信号成分xj´(k)の推定値^xj´(k)と、事後共分散行列更新部550からの同時無相関化された事後共分散行列Σj´(k)と、空間共分散行列更新部570からの同時無相関化された空間共分散行列Λjの推定値^Λj(ただし、例外として、分散パラメータ更新部560における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λjの初期推定値^Λj (0))と、に基づいて、分散パラメータvj(k)の推定値^vj(k)を次の(175)式により更新する(ステップS67)。
ただし、分散パラメータ更新部560における最初の処理の際には、(175)式の右辺において、^Λjを^Λj (0)に置き換える。
空間共分散行列更新部570は、シングルチャネルウィーナーフィルタ部540からの同時無相関化された源信号成分xj´(k)の推定値^xj´(k)と、事後共分散行列更新部550からの同時無相関化された事後共分散行列Σj´(k)と、分散パラメータ更新部560からの分散パラメータvj(k)の推定値^vj(k)と、に基づいて、同時無相関化された空間共分散行列Λjの推定値^Λjを次の(176)式により更新する(ステップS68)。
信号パラメータ分析部520は、シングルチャネルウィーナーフィルタ部540からの同時無相関化された源信号成分xj´(k)の推定値^xj´(k)と、事後共分散行列更新部550からの同時無相関化された事後共分散行列Σj´(k)と、分散パラメータ更新部560からの分散パラメータvj(k)の推定値^vj(k)と、空間共分散行列更新部570からの同時無相関化された空間共分散行列Λjの推定値^Λjと、(信号パラメータ分析部520における同時無相関化行列Pの推定値^Pの最初の更新の際には、更に、初期化部からの同時無相関化行列Pの初期推定値^P(0)と、)を受け取って、次の(177)式により同時無相関化行列Pの推定値^Pを更新する(ステップS69)。ただし、eiは単位行列の第i列を表し、行列Aに対して[A]iは、Aの第i列を表す。
ただし、信号パラメータ分析部520における同時無相関化行列Pの推定値^Pの最初の更新の際には、(177)式の右辺において、^Pを^P(0)で置き換える。なお、信号パラメータ分析部520は、(177)式を複数回連続して適用することにより同時無相関化行列Pの推定値^Pを更新してもよい。
図示しない収束判定部は、収束したかどうかの判定を行う(ステップS70)。収束判定部が収束していないと判定した場合(ステップS70:No)、同時無相関化部530によるステップS64の処理に戻って、処理が継続される。
一方、収束判定部が収束したと判定した場合(ステップS70:Yes)、同時無相関化部530における処理(ステップS71)と、シングルチャネルウィーナーフィルタ部540における処理(ステップS72)と、同時無相関化逆変換部580における処理(ステップS73)と、が行われる。同時無相関化部530及びシングルチャネルウィーナーフィルタ部540における処理については上述の通りであるから、同時無相関化逆変換部580での処理について説明する。同時無相関化逆変換部580は、信号パラメータ分析部520からの同時無相関化行列Pの推定値^Pと、シングルチャネルウィーナーフィルタ部540からの同時無相関化された源信号成分xj´(k)の推定値^xj´(k)と、に基づいて、源信号成分xj(k)の推定値^xj(k)を次の(178)式により計算し、出力する(ステップS73)。
[信号分析装置の処理の背景]
次に、本第5の実施形態に係る信号分析装置501の処理の背景について説明する。本第5の実施形態においては、観測信号ベクトルy(k)は、J個の源信号成分x1(k),・・・,xJ(k)の和として、次の(179)式でモデル化される。
源信号成分xj(k)は、(97)式のように平均0、共分散行列Sj(k)=vj(k)Rjの複素ガウス分布に従うとする。ここで、vj(k)はj番目の源信号のサンプル点ごとの分散をモデル化する分散パラメータであり、Rjはj番目の源信号の空間的特性をモデル化するパラメータである空間共分散行列である。分散パラメータvj(k)は正であり、空間共分散行列Rjは正定値エルミート行列であると仮定する。J個の源信号成分x1(k),・・・,xJ(k)は互いに無相関であるとする。
観測信号ベクトルy(k)から源信号成分xj(k)を推定する方法としては、マルチチャネルウィーナーフィルタに基づく方法がある。この方法では、源信号成分xj(k)の推定値^xj(k)を次の(180)式により計算する。
本第5の実施形態のように、分散パラメータvj(k)及び空間共分散行列Rjが未知の状況で(180)式により源信号成分xj(k)の推定値^xj(k)を計算する(すなわち信号分離を実現する)ためには、これらのパラメータを推定する必要がある。
[従来の信号分析方法]
分散パラメータvj(k)及び空間共分散行列Rjが未知の状況で、分散パラメータvj(k)と、空間共分散行列Rjと、源信号成分xj(k)と、を同時に推定することにより、信号分離を実現する方法が、非特許文献1に開示されている。
従来技術では、分散パラメータvj(k)及び空間共分散行列Rjを最尤推定する。この最尤推定は、下記の(181)式の最適化問題として定式化される。
ここで、Ψはパラメータの全体を表し、分散パラメータvj(k)(j=1,・・・,J;k=1,・・・,K)及び空間共分散行列Rj(j=1,・・・,J)からなる。L1(Ψ)は、次の(182)式の対数尤度関数を表す。
下の(183)式は、行列Aが正定値エルミート行列であることを表す。
従来技術では、EMアルゴリズムに基づいて最尤推定を行う。このEMアルゴリズムでは、EステップとMステップとが交互に反復される。
Eステップでは、パラメータΨの現在の推定値^Ψに基づく、源信号成分xj(k)の事後確率p(xj(k)|y(k),^Ψ)を計算する。この事後確率は、次の(184)式のように複素ガウス分布であることが示される。
ただし、平均^xj(k)及び共分散行列Σj(k)は、パラメータの現在の推定値^Ψを用いて次の(185)式及び(186)式のように表される(^vj(k)は分散パラメータvj(k)の現在の推定値であり、^Rjは空間共分散行列Rjの現在の推定値である)。
したがって、Eステップでは、平均^xj(k)及び共分散行列Σj(k)を、(185)式及び(186)式により更新する。なお、(185)式は、(180)式のマルチチャネルウィーナーフィルタによる源信号成分xj(k)の推定値に他ならない。以下、^xj(k)を源信号成分xj(k)の推定値と呼び、Σj(k)を事後共分散行列と呼ぶ。
Mステップでは、Eステップで更新された源信号成分xj(k)の推定値^xj(k)及び事後共分散行列Σj(k)に基づいて、分散パラメータvj(k)の推定値^vj(k)及び空間共分散行列Rjの推定値^Rjを更新する。これらは、下記の更新則((187)式及び(188)式)により更新される。
しかしながら、従来技術では、(185)式の源信号成分xj(k)の推定値^xj(k)の更新と、(186)式の事後共分散行列Σj(k)の更新と、において、逆行列計算と行列乗算とをサンプル点ごとに行う必要があった。具体的には、(189)式に示す逆行列計算(反復ごとサンプル点ごとに1回)と、(190)式に示す行列乗算(反復ごとサンプル点ごとに2J回)とを行う必要があった。
そのため、従来技術では膨大な計算量を要していたため、リアルタイムでの信号分離や、センサネットワーク等における分散処理など、サンプル点の数Kが大きい(観測信号長が長い)場合や計算資源が限られている状況への適用が困難であった。
[第5の実施形態に係る信号分析方法]
これに対し、本第5の実施形態では、J個の源信号成分x1(k),・・・,xJ(k)の同時無相関化に基づいて、Eステップにおける源信号成分xj(k)の推定値^xj(k)及び事後共分散行列Σj(k)の更新を少ない計算量で実現することにより、少ない計算量での信号分離を可能にする。以下、本第5の実施形態における考え方について説明する。
いま仮に、(185)式及び(186)式における源信号成分xj(k)の共分散行列vj(k)Rjがいずれも対角行列であるとすると、(185)式及び(186)式における逆行列演算及び行列乗算は単なる対角要素ごとの逆数演算及び乗算に帰着するため、少ない計算量で計算できる。しかしながら、通常、ベクトルxj(k)の異なる要素(つまり異なる位置で観測されたj番目の源信号)は互いに相関を持つため、その共分散行列vj(k)Rjの非対角要素は0ではない、すなわち共分散行列vj(k)Rjは対角行列ではない。
そこで本第5の実施形態では、J個の源信号成分x1(k),・・・,xJ(k)を、正則行列Pのエルミート転置PHによる次の(191)式の変換により同時無相関化することを考える。
変換後のxj´(k)の共分散行列は、次の(192)式及び(193)式で与えられる。
変換後のxj´(k)(j=1,・・・,J)がすべて無相関になる条件、すなわち、(193)式で与えられるこれらの共分散行列がすべて対角行列になる条件は、次の(194)式〜(196)式である。
既に述べたように、信号数がJ=2の場合には、一般化固有値問題を解くことにより、(194)式〜(196)式を満たす正則行列P及び対角行列Λ1,・・・,ΛJを得ることができる。これに対し、本第5の実施形態では、J=2の場合に限らず一般の信号数Jの場合に適用可能なアプローチとして、未知のJ個の空間共分散行列R1,・・・,RJが、ある正則行列Pにより同時対角化されるという制約条件、すなわち、ある正則行列P及びあるJ個の対角行列Λ1,・・・,ΛJに対し、(194)式〜(196)式が成立するという制約条件の下、正則行列P及び対角行列Λ1,・・・,ΛJを最尤推定する。上述の制約条件は、ある正則行列P及びあるJ個の対角行列Λ1,・・・,ΛJに対し、次の各式((197)式〜(199)式)が成立することと同値である。
本第5の実施形態では、空間共分散行列R1,・・・,RJは、パラメータである正則行列Pとパラメータである対角行列Λ1,・・・,ΛJとを用いて、(197)式〜(199)式のように表される(モデル化される)と仮定する。この場合、空間共分散行列R1,・・・,RJを推定することは、行列P(同時無相関化を実現する行列であるから、以下、同時無相関化行列と呼ぶ)及び対角行列Λ1,・・・,ΛJ(同時無相関化行列Pを用いて(194)式〜(196)式のように空間共分散行列R1,・・・,RJを変換して得られる行列であるから、以下、同時無相関化された空間共分散行列と呼ぶ)を推定することに帰着する。
そこで、本第5の実施形態では、同時無相関化行列Pと、同時無相関化された空間共分散行列Λjと、分散パラメータvj(k)と、を最尤法により同時推定する。これにより、一般の信号数Jに対しても、源信号成分の同時無相関化(源信号成分の共分散行列の同時対角化)が実現し、少ない計算量での信号分離が可能になる。
上述の最尤推定は、次の最適化問題((200)式)として定式化される。
ここで、Θはすべてのパラメータの全体を表し、分散パラメータvj(k)(j=1,・・・,J;k=1,・・・,K)と、同時無相関化行列Pと、同時無相関化された空間共分散行列Λj(j=1,・・・,J)と、からなる。L2(Θ)は、次の(201)式の対数尤度関数であり、(182)式に(197)式〜(199)式を代入することにより得られる。なお、GL(I,C)はすべての正則なI次複素行列からなる集合である。
本第5の実施形態では、EMアルゴリズムにより、上述の最尤推定を実現する。このEMアルゴリズムはEステップとMステップとからなる。
Eステップでは、パラメータΘの現在の推定値^Θに基づいて、源信号成分の事後確率を更新する。従来技術の場合と同様に、xj(k)の事後確率p(xj(k)|y(k),^Θ)は、(184)式の右辺の複素ガウス分布であり、その平均は、(185)式で、その共分散行列は、(186)式で与えられることが示せる。ただし、従来技術の場合とは異なり、(185)式及び(186)式における空間共分散行列Rjの推定値^Rjは、同時無相関化行列Pの推定値^P及び同時無相関化された空間共分散行列Λjの推定値^Λjを用いて、次の(202)式〜(204)式のように表される。
(185)式及び(186)式に(202)式〜(204)式を代入することにより、次の(205)式〜(208)式を得る。
したがって、源信号成分xj(k)の推定値^xj(k)を同時無相関化行列Pの推定値^Pにより基底変換して得られる^xj´(k)(以下、同時無相関化された源信号成分の推定値と呼ぶ)は、(173)式により更新でき、事後共分散行列Σj(k)を同時無相関化行列Pの推定値^Pにより基底変換して得られるΣj´(k)(以下、同時無相関化された事後共分散行列と呼ぶ)は、(174)式により更新できる。その際、同時対角化により、逆行列演算及び行列乗算は、対角行列に対する逆行列演算及び行列乗算になっているため、その計算量は、O(I3)ではなくO(I)であり少ない計算量で計算できる。
Mステップでは、分散パラメータvj(k)の推定値^vj(k)と、同時無相関化された空間共分散行列Λjの推定値^Λjと、同時無相関化行列Pの推定値^Pと、を下記のQ関数の最大化に基づいて更新する((209)式及び(210)式参照)。ただし、(209)式及び(210)式ではΘに依らない定数項は無視してあり、したがって等号はその両辺がΘに依らない定数の差を除いて等しいことを表している。
(210)式のQ関数を分散パラメータvj(k)に関して偏微分して0とおくと、次の(211)式及び(212)式を得る。
これで、(175)式の分散パラメータvj(k)の推定値^vj(k)の更新式が導かれた。
また、(210)式のQ関数を同時無相関化された空間共分散行列Λjに関して偏微分したものが零行列に等しいと置くと、次の(213)式及び(214)式を得る。
これで、(176)式の同時無相関化された空間共分散行列Λjの推定値^Λjの更新則が導かれた。なお、(213)式において、絶対値二乗|・|2はベクトルの要素ごとに計算するものとし、ベクトルαに対してdiag(α)は、αの要素を対角要素とする対角行列を表す。
(209)式のQ関数を同時無相関化行列Pの複素共役P*に関して偏微分すると、次の(215)式を得る。
これが零行列に等しいと置いて両辺のvecを取ると、次の(216)式を得る。
ここで、行列Aに対し、vec(A)は、行列Aの列ベクトルを縦に並べて得られる列ベクトルを表す。これに関し、次の(217)式の公式が成り立つ。
また、「×」を丸で囲んだ記号は、クロネッカ積を表す。ブロック対角行列の性質より、(216)式は次の(218)式〜(220)式と等価である。
そこで、不動点法に基づいて、(177)式により、同時無相関化行列Pの推定値^Pを更新すれば良い。このように不動点法に基づくことで、勾配法や自然勾配法とは異なり学習率の調整が不要であるため、安定な最適化が可能である。
以下、事後確率の詳細な導出を述べる。隠れ変数x1(k),・・・,xJ−1(k)の対数事後確率は、次の(221)式〜(223)式のように書き下せる。
ここで、等号「=」の上に「c」を付した記号は、両辺がx1(k),・・・,xJ−1(k)に依らない定数の差を除いて等しいことを表す。行列Aは、I次の単位行列IをJ−1個横に並べた次の(224)式の行列を表す。また、空間共分散行列Rj(j=1,・・・,J)は(197)式〜(199)式で与えられる。なお、本明細書では、Iを観測信号の個数(スカラー)と単位行列(行列)との2通りの意味で用いているが、Iがスカラーであるか行列であるかは文脈より明らかであるから、これにより混乱が生じるおそれはない。
~x(k)は、x1(k),・・・,xJ−1(k)を縦に並べた次の(225)式のベクトルを表す。
(223)式より、隠れ変数x1(k),・・・,xJ−1(k)の事後確率が次の(226)式の複素ガウス分布であることが導かれる。
ここで、共分散行列~Σ(k)及び平均^~x(k)は次の(227)式〜(230)式及び(231)式〜(233)式で与えられる。但し、(227)式から(228)式への式変形の際に逆行列補題を用いた。
後に分かるように、Eステップでは、行列~Σ全体を計算する必要は実はなく、行列~ΣをI次の正方行列を要素とするJ−1次の正方ブロック行列とみなしたときの対角ブロックのみを計算すれば良い。j番目の対角ブロックをΣj(k)と書くことにすると、これは次の(234)式のように計算できる。明らかに、J=2の場合には、(234)式は(101)式と等価である。
また、^xj(k)を、^~x(k)をI次元縦ベクトルを要素とするサイズ(J−1)×1のブロック行列とみなしたときの(j,1)ブロックと定義すると、これは次の(235)式のように表される。
^xj(k)及びΣj(k)は次の(236)式及び(237)式のようにも表されることに注意する。
ただし、上式においては、j=1,・・・,J−1である。また、E(・|y(k))は、隠れ変数の事後確率p(x1(k),・・・,xJ−1(k)|y(k),Θ)に関する期待値((238)式)である。
これにならい、j=Jに対しても、^xJ(k)及びΣJ(k)を次の(239)式〜(243)式及び(244)式〜(251)式のように定義する。
よって、以下の(252)式及び(253)式は、j=1,・・・,Jに対して成立する。
以下、Q関数の詳細な導出について説明する。観測信号ベクトルy(k)と隠れ変数x1(k),・・・,xJ−1(k)との同時分布の対数は、次の(254)式〜(256)式のように書き下される。ただし、(254)式〜(256)式では定数項は無視してあり、したがって等号はその両辺が定数の差を除いて等しいことを表している。
ここで、次の各式((257)式〜(262)式)が成立することに注意する。
したがって、Q関数は次式((263)式〜(266)式)で与えられる。(264)式は(209)式に他ならない。
ここで、(266)式における絶対値の二乗|・|2は、要素ごとに取るものとする。また、等号「=」の上に「c」を付した記号は、両辺がΘに依らない定数の差を除いて等しいことを表す(先と異なる意味で記号が用いられていることに注意されたい)。
[第5の実施形態の変形例1]
次に、第5の実施形態の変形例1について説明する。本第5の実施形態の変形例1では、信号パラメータ分析部における処理の変形例について説明する。第5の実施形態では、信号パラメータ分析部は、不動点法に基づいて同時無相関化行列Pの推定値^Pを更新していた。これに対し、第5の実施形態の変形例1に掛かる信号パラメータ分析部は、不動点法以外の任意の最適化法(例えば、勾配法、自然勾配法、ニュートン法など)に基づいて、同時無相関化行列Pの推定値^Pを更新する。
例えば、勾配法の場合、本第5の実施形態の変形例1に係る信号パラメータ分析部は、同時無相関化行列Pの推定値^Pを次の(267)式により更新する(これは、(215)式から直ちに導かれる)。
ただし、信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、(267)式の右辺において、^Pを^P(0)で置き換える。μは学習率と呼ばれる正の数であり、例えばμ=0.05などと定めればよい。
また、自然勾配法の場合、本第5の実施形態の変形例1に係る信号パラメータ分析部は、同時無相関化行列Pの推定値を次の(268)式により更新する(これも(215)式から直ちに導かれる)。
ただし、信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、(268)式の右辺において、^Pを^P(0)で置き換える。μは、学習率と呼ばれる正の数であり、例えば、μ=0.05などと定めればよい。自然勾配法によれば、空間の計量が適切に考慮されるため、一般に勾配法よりも安定な更新が可能であることが知られている。
[第5の実施形態の変形例2]
次に、第5の実施形態の変形例2について説明する。本第5の実施形態の変形例2では、信号パラメータ分析部における処理の変形例について説明する。第5の実施形態や第5の実施形態の変形例1では、信号パラメータ分析部は、Q関数に基づいて同時無相関化行列Pの推定値^Pを更新していた。これに対し、本第5の実施形態の変形例2に係る信号パラメータ分析部は、直接、対数尤度関数L2(Θ)に基づいて、同時無相関化行列Pの推定値^Pを更新する。
本第5の実施形態の変形例2に係る信号パラメータ分析部は、同時無相関化部からの同時無相関化された観測信号ベクトルy´(k)と、空間共分散行列更新部からの同時無相関化された空間共分散行列Λjの推定値^Λjと、分散パラメータ更新部からの分散パラメータvj(k)の推定値^vj(k)と、(信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、更に、初期化部からの同時無相関化行列Pの初期推定値^P(0)と、)に基づいて、同時無相関化行列Pの推定値^Pを次式((269)式)により更新する(Eは単位行列を表すものとする)。
ただし、信号パラメータ分析部における^Pの最初の更新の際には、(269)式の右辺において、^Pを^P(0)で置き換える。なお、信号パラメータ分析部は、(269)式を複数回連続して適用することにより、同時無相関化行列Pの推定値^Pを更新してもよい。
(269)式の対数尤度関数L2(Θ)に基づく同時無相関化行列Pの推定値^Pの更新則の導出について説明する。(201)式の対数尤度関数を同時無相関化行列Pの複素共役P*に関して偏微分すると次の(270)式を得る。
(270)式が零行列に等しいと置いて、両辺のvecを取ると、次式((271)式)を得る。
(271)式の右辺の行列は、次の(272)式及び(273)式のように変形できる。
よって次式((274)式及び(275)式)を得る。
(275)式に(276)式を代入することにより、(269)式の不動点法による同時無相関化行列Pの推定値^Pの更新則が得られる。
このように、不動点法に基づくことで、勾配法や自然勾配法とは異なり学習率の調整が不要であるため、安定な最適化が可能である。
また、第5の実施形態の変形例1と同様に、本第5の実施形態の変形例2においても、不動点法に限らず任意の最適化法(勾配法、自然勾配法、ニュートン法など)により同時無相関化行列Pの推定値^Pを更新してよい。以下、そのような例について説明する。
勾配法の場合、本第5の実施形態の変形例2に係る信号パラメータ分析部は、観測信号ベクトル作成部からの観測信号ベクトルy(k)と、空間共分散行列更新部からの同時無相関化された空間共分散行列Λjの推定値^Λjと、分散パラメータ更新部からの分散パラメータvj(k)の推定値^vj(k)と、(信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、更に、初期化部からの同時無相関化行列Pの初期推定値^P(0)と、)に基づいて、同時無相関化行列Pの推定値^Pを次の(277)式により更新する(これは、(270)式より直ちに導かれる)。
ただし、信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、(277)式の右辺において、^Pを^P(0)で置き換える。なお、信号パラメータ分析部は、(277)式を複数回連続して適用することにより同時無相関化行列Pの推定値^Pを更新してもよい。また、信号パラメータ分析部は、(277)式に(276)式を代入することにより得られる、同時無相関化された観測信号ベクトルy´(k)を用いて表した式を、(277)式の代わりに用いても良い。μは、学習率と呼ばれる所定の正の数であり、例えば、μ=0.05などと定めればよい。この勾配法の場合には、(277)式の1回の適用につき逆行列計算が1回のみで行列乗算は不要(対角行列に対する逆行列計算はカウントしていない)であるため、不動点法の場合と比べてより少ない計算量で同時無相関化行列Pの推定値^Pを更新できる。
また、自然勾配法の場合、本第5の実施形態の変形例2に係る信号パラメータ分析部は、観測信号ベクトル作成部からの観測信号ベクトルy(k)と、空間共分散行列更新部からの同時無相関化された空間共分散行列Λjの推定値^Λjと、分散パラメータ更新部からの分散パラメータvj(k)の推定値^vj(k)と、(信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、更に、初期化部からの同時無相関化行列Pの初期推定値^P(0)と、)に基づいて、同時無相関化行列Pの推定値^Pを次の(278)式により更新する(これも(270)式より直ちに導かれる)。
ただし、信号パラメータ分析部における同時無相関化行列Pの推定値^Pの最初の更新の際には、(278)式の右辺において、^Pを^P(0)で置き換える。なお、信号パラメータ分析部は、(278)式を複数回連続して適用することにより同時無相関化行列Pの推定値^Pを更新してもよい。また、信号パラメータ分析部は、(278)式に(276)式を代入することにより得られる、同時無相関化された観測信号ベクトルy´(k)を用いて表した式を、(278)式の代わりに用いても良い。μは、学習率と呼ばれる所定の正の数であり、例えば、μ=0.05などと定めればよい。また、自然勾配法によれば、空間の計量が適切に考慮されるため、一般に勾配法よりも安定な更新が可能であることが知られている。
[第5の実施形態の変形例3]
次に、第5の実施形態の変形例3について説明する。第5の実施形態の変形例3では、第3の実施形態にならい、周波数ごとに第5の実施形態を適用する。
第5の実施形態の変形例3においては、J個の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得されたI個(Iは2以上の整数)の観測信号yi(t)が信号分析装置に入力されるものとする。観測信号yi(t)は、周波数分解され、周波数ごとに第5の実施形態に記載の方法で信号分離がなされ、最後に時間領域に変換される。
図5及び図6を用いて、第5の実施形態の変形例3に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。
周波数領域変換部310(ステップS40、ステップS41)及び時間領域変換部330(ステップS43)での処理は、第3の実施形態と同じである。以下、相違点である周波数領域信号分離部320(ステップS42)での処理について説明する。
周波数領域信号分離部320は、周波数領域変換部310から(279)式に示す時間周波数領域の観測信号を受け取り、fごとに第5の実施形態の方法により周波数領域での信号分離を行うことにより、(280)式に示す時間周波数領域の源信号成分の推定値を計算し、時間領域変換部330に渡す(ステップS42)。
すなわち、f=1に対しては、(281)式に示す観測信号を入力として、第5の実施形態の方法により(282)式の源信号成分の推定値を得る。
そして、f=2に対しては、(283)式に示す観測信号を入力として、第5の実施形形態の方法により(284)式の源信号成分の推定値を得る。
以下、同じ処理を各fに対して繰り返し、最後に、f=Fに対しては、(285)式に示す観測信号を入力として、第5の実施形態の方法により(286)式の源信号成分の推定値を得る。
こうして得られた時間周波数領域の源信号成分の推定値((287)式参照)を時間領域変換部330に渡す。
[第5の実施形態の効果]
本発明の効果を確認するために、従来の方法及び第5の実施形態を用いた音源分離実験について説明する。
図14は、従来の方法及び第5の実施形態を用いた音源分離実験の実験条件を示す図である。図15は、従来の方法及び第5の実施形態を用いた音源分離実験結果(real time factor)を示す図である。図16は、従来の方法及び第5の実施形態を用いた音源分離実験結果(音源分離性能)を示す図である。
8秒の英語音声に、図14に示す実験室で測定された残響インパルス応答を畳み込んで残響音声を生成し、複数人の残響音声を加算することで観測信号を生成した。このとき、観測信号に対して、各方法を用いて信号分離を行った場合の性能を図15及び図16に示す。図15は、real time factor(観測信号長を1としたときの、処理に要した時間)である。また、図16は、音源分離性能である信号対歪み比(値が大きいほど音源分離性能が高い)である。
図15及び図16に示すように、本第5の実施形態に係る方法によれば、音源数が3個以上である場合にも、420倍以上の計算時間削減が実現でき、従来の方法を上回る音源分離性能を得ることができた。
このように、本第5の実施形態によれば、従来の方法と比べて計算時間を大きく削減できるのに加え、従来の方法よりも高い信号分離性能を実現できる。本第5の実施形態は、信号数が2個でない場合にも適用でき、信号数が3個以上の場合にも適用できる。
[第6の実施形態]
次に、第6の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。第6の実施形態では、第4の実施形態の変形例を実施するための具体的な手法について説明する。
本第6の実施形態では、第5の実施形態と同様に、一般にJ個の源信号が混在し、分散パラメータvj(k)及び空間共分散行列Rjが未知の状況で、分散パラメータvj(k)と、空間共分散行列Rjをモデル化するパラメータである同時無相関化された空間共分散行列Λj及び同時無相関化行列Pと、源信号成分xj(k)と、を同時に推定することにより、信号分離を実現する方法について説明する。
なお、第5の実施形態は、EMアルゴリズムに基づいていたのに対し、本第6の実施形態は補助関数法に基づいている。本第6の実施形態では、それぞれ異なる位置で取得されたI個の観測信号yi(k)が信号分析装置に入力されるものとする。
[信号分析装置の構成及び処理]
図17及び図18を用いて、第6の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図17は、第6の実施形態に係る信号分析装置の構成の一例を示す図である。図18は、第6の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
図17に示すように、信号分析装置601は、初期化部(図示しない)、観測信号ベクトル作成部610、同時無相関化部630、シングルチャネルウィーナーフィルタ部640、観測共分散行列更新部650、分散パラメータ更新部660、空間共分散行列更新部670、信号パラメータ分析部620、収束判定部(図示しない)及び同時無相関化逆変換部680を有する。
まず、信号分析装置の各部の概要について説明する。初期化部(図示しない)の処理(ステップS81)、観測信号ベクトル作成部610の処理(ステップS82,S83)、同時無相関化部630の処理(ステップS84)及びシングルチャネルウィーナーフィルタ部640の処理(ステップS85)については、第5の実施形態と同様であるから、説明を省略する。
観測共分散行列更新部650は、分散パラメータ更新部660からの分散パラメータvj(k)の推定値^vj(k)(ただし例外として、観測共分散行列更新部650における最初の処理の際には、初期化部からの分散パラメータvj(k)の初期推定値^vj (0)(k))と、空間共分散行列更新部670からの同時無相関化された空間共分散行列Λjの推定値^Λj(ただし例外として、観測共分散行列更新部650における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λjの初期推定値^Λj (0))と、を受け取って、同時無相関化された観測共分散行列U´(k)を次式(288)式により更新する(ステップS86)。
ただし、観測共分散行列更新部650における最初の処理の際には、上式において、^vj(k)を^vj (0)(k)に、^Λjを^Λj (0)に、それぞれ置き換える。
分散パラメータ更新部660は、シングルチャネルウィーナーフィルタ部640からの同時無相関化された源信号成分xj´(k)の推定値^xj´(k)と、観測共分散行列更新部650からの同時無相関化された観測共分散行列U´(k)と、空間共分散行列更新部670からの同時無相関化された空間共分散行列Λjの推定値^Λj(ただし例外として、分散パラメータ更新部660における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λjの初期推定値^Λj (0))と、を受け取って、分散パラメータvj(k)の推定値^vj(k)を次式((289)式)により更新する(ステップS87)。
ただし、分散パラメータ更新部660における最初の処理の際には、上式において、^Λjを^Λj (0)で置き換える。
空間共分散行列更新部670は、シングルチャネルウィーナーフィルタ部640からの同時無相関化された源信号成分xj´(k)の推定値^xj´(k)と、観測共分散行列更新部650からの同時無相関化された観測共分散行列U´(k)と、分散パラメータ更新部660からの分散パラメータvj(k)の推定値^vj(k)と、を受け取って、同時無相関化された空間共分散行列Λjの推定値^Λjを次の(290)式により更新する(ステップS88)。
信号パラメータ分析部620は、第5の実施形態の変形例2と同様の処理を行って、同時無相関化行列Pの推定値^Pを更新する(ステップS89)。
図示しない収束判定部は、収束したかどうかの判定を行う(ステップS90)。収束判定部が収束していないと判定した場合(ステップS90:No)、同時無相関化部630での処理(ステップS84)に戻って、処理が継続される。
一方、収束判定部が収束したと判定した場合(ステップS90:Yes)、同時無相関化部630における処理(ステップS91)、シングルチャネルウィーナーフィルタ部640における処理(ステップS92)、及び、同時無相関化逆変換部680における処理(ステップS93)が、第5の実施形態と同様にして行われることにより、源信号成分xj(k)の推定値^xj(k)が計算され、出力される。
[信号分析装置の処理の背景]
次に、第6の実施形態に係る信号分析装置の背景について説明する。第5の実施形態で述べたように、分散パラメータvj(k)及び空間共分散行列Rjが未知の状況で、(180)式により、源信号成分xj(k)の推定値^xj(k)を計算する(すなわち信号分離を実現する)ためには、分散パラメータvj(k)及び空間共分散行列Rjを推定する必要がある。
補助関数法に基づいて、(181)式の最適化問題を解くことにより、分散パラメータvj(k)及び空間共分散行列Rjを最尤推定することができる。
この補助関数法では、補助変数Φ及び次式を満たす補助関数L1 −(Ψ,Φ)に基づいて、補助関数L1 −(Ψ,Φ)を補助変数Φについて最大化することによりΦを更新するステップと、補助関数L1 −(Ψ,Φ)が減少しないようにパラメータΨを更新するステップと、が交互に反復される((291)式参照)。
本第6の実施形態では、補助関数として(292)式のL1 −(Ψ,Φ)を用いる。
ただし、右辺における「定数」は、パラメータΨ及び補助変数Φによらない定数である。(292)式が(291)式の補助関数の条件を満たすことは、参考文献3及び参考文献4と同様に示すことができる。なお、補助変数Φは、U(k)及びGj(k)からなり、U(k)は正定値エルミート行列であり、Gj(k)は(293)式を満たすものとする。
補助関数L1 −(Ψ,Φ)を補助変数Φについて最大化することによりΦを更新するステップでは、分散パラメータvj(k)の推定値^vj(k)及び空間共分散行列Rjの推定値^Rjに基づいて、U(k)及びGj(k)を(294)式及び(295)式により更新する。U(k)は観測信号ベクトルy(k)の共分散行列に他ならないから、観測共分散行列と呼ぶ。
補助関数L1 −(Ψ,Φ)が減少しないようにパラメータΨを更新するステップでは、上で更新した補助変数に基づいて、次の(296)式及び(297)式により分散パラメータvj(k)の推定値^vj(k)と、空間共分散行列Rjの推定値^Rjと、を更新する。
ここで、A#Bは、(298)式で定義されるAとBとの幾何平均である。
また、行列Cj及び行列Djは、それぞれ(299)式及び(300)式により計算される。
(294)式〜(300)式の更新式の導出は、参考文献3及び参考文献4と同様にできる。更に、補助変数Gj(k)を代入消去することにより、補助変数Gj(k)が現れない更新則を導くことができる。この場合の更新則(1回の反復)を以下の(301)式〜(305)式に示す。明らかに、J=2の場合には、(301)式〜(305)式による更新は、(156)式〜(159)式による更新と等価である。
しかしながら、上記の方法では、(302)式〜(304)式に含まれる、逆行列計算U(k)−1をサンプル点ごとに行う必要があった。そのため、上記の方法は膨大な計算量を要する。
これに対し、本第6の実施形態では、第5の実施形態と同様に、J個の源信号成分x1(k),・・・,xJ(k)の同時無相関化に基づいて、逆行列計算U(k)−1を少ない計算量で計算できる対角行列に対する逆行列計算に帰着させることにより、少ない計算量での信号分離を実現する。具体的には、本第6の実施形態では、第5の実施形態と同様に、J個の空間共分散行列が(197)式〜(199)式の形に表されるという制約条件の下、同時無相関化行列P及び同時無相関化された空間共分散行列Λ1,・・・,ΛJを最尤推定する。この最尤推定は、(200)式の最適化問題として定式化される。本第6の実施形態では、補助関数法に基づいて、上述の最尤推定を実現する。この方法では、補助関数法の1回の反復を適用することにより分散パラメータvj(k)の推定値^vj(k)及び同時無相関化された空間共分散行列Λjの推定値^Λjを更新するステップと、同時無相関化行列Pの推定値^Pを更新するステップと、が交互に反復される。
補助関数法に基づいて分散パラメータvj(k)の推定値^vj(k)及び同時無相関化された空間共分散行列Λjの推定値^Λjを更新するステップは、(306)式の補助関数L2 −(Θ,Φ)に基づいている。
ただし、右辺の「定数」はパラメータΘ及び補助変数Φによらない定数である。
先と同様に補助変数Gj(k)を代入消去すると、更新則は下記の(307)式〜(311)式のように与えられる。ただし、(310)式の右辺における斜体のDは、行列の非対角成分を0で置き換える演算子である。
これは、更に次の(312)式〜(315)式に書き直せる。
同時無相関化行列Pの推定値^Pの更新則の導出については、第5の実施形態の変形例2と同様であるから説明を省略する。以上で、第6の実施形態における更新則が導かれた。
[第6の実施形態の変形例]
第6の実施形態では対数尤度に基づいて同時無相関化行列Pの推定値^Pを更新するが、補助関数に基づいて同時無相関化行列Pの推定値^Pを更新してもよい。詳細は、第5の実施形態及び第5の実施形態の変形例1と同様だから省略する。
[第6の実施形態の効果]
第6の実施形態によれば、同時無相関化に基づき、計算量を大幅に削減できる。また、EMアルゴリズムは反復回数が少ない場合には必ずしも高い信号分離性能を実現できないが、第6の実施形態によれば、補助関数法に基づき、反復回数が少ない場合にも高い信号分離性能を実現できる。本第6の実施形態は、信号数が3個以上でも適用できる。
[第7の実施形態]
本第7の実施形態では、観測信号に含まれる源信号の総数LがL≧3の場合でも少ない計算量で信号分離を可能にする方法を示す。第5の実施形態や第6の実施形態でもそのような方法を示したが、本実施形態ではもう一つの実現方法を示す。
第5の実施形態や第6の実施形態は、「すべての時間周波数点において観測信号は、L個すべての源信号を含む」というモデル化((179)式参照)の下、L個すべての源信号の同時無相関化(L個すべての源信号の空間共分散行列の同時対角化)に基づいて、L≧3の場合でも少ない計算量で信号分離を実現していた。この同時無相関化(同時対角化)は、「L個すべての源信号の空間共分散行列が1つの同時無相関化行列によって同時対角化できる」という制約条件((197)式〜(199)式参照)に基づいていた。しかしながら、実際にはこの制約条件は現実に即していない場合もあり、そのような場合には、第5の実施形態や第6の実施形態ではモデル化精度が低下し、結果として信号分離性能が低下するおそれがあった。なお、第5の実施形態や第6の実施形態では、L個すべての源信号の同時無相関化(L個すべての源信号の空間共分散行列の同時対角化)に基づいているため、同時無相関化される源信号の最大個数(同時対角化される源信号の空間共分散行列の最大個数)をJとすれば、J=Lである。したがって、第5の実施形態や第6の実施形態では、これらをともにJで表していた。
一方、本第7の実施形態では、源信号の時間周波数領域におけるスパース性に着目する。源信号は、時間周波数領域におけるスパース性、すなわち、「少数の時間周波数成分にエネルギーが集中する」という性質を持つことが多い(例えば、参考文献5「O. Yilmaz and S. Rickard, “Blind separation of speech mixtures via time-frequency masking”, IEEE Transactions on Signal Processing, vol. 52, no. 7, pp. 1830−1847, Jul. 2004.」参照)。この場合、各時間周波数点において観測信号は、L個すべての源信号を含むのではなく、L個よりも少ない源信号しか含まないとみなせることが多い。
そこで、本第7の実施形態では、「各時間周波数点において観測信号は高々2個の源信号しか含まない」というモデル化の下、各時間周波数点において観測信号に含まれる高々2個の源信号の同時無相関化(各時間周波数点において観測信号に含まれる高々2個の源信号の空間共分散行列の同時対角化)に基づいて、L≧3の場合でも少ない計算量で信号分離を実現する。一般化固有値問題に基づいて2個以下の源信号は常に一つの同時無相関化行列によって同時無相関化できる(2個以下の源信号の空間共分散行列は常に一つの同時無相関化行列によって同時対角化できる)から、本第7の実施形態では、第5の実施形態や第6の実施形態とは異なり、源信号の空間共分散行列に特段の制約条件を課す必要がない。
したがって、本第7の実施形態によれば、「L個すべての源信号の空間共分散行列が1つの同時無相関化行列によって同時対角化できる」という制約条件が現実に即していない場合でも高精度なモデル化が可能であり、結果として、そのような場合でも高い信号分離性能を実現可能である。なお、本第7の実施形態では、同時無相関化される源信号の最大個数J(同時対角化される源信号の空間共分散行列の最大個数J)はJ=2である。
本第7の実施形態では、L個(Lは2以上の自然数)の源信号が混在する状況において、それぞれ異なる位置で取得されたI個(Iは2以上の自然数)の観測信号yi(t)が信号分析装置に入力されるものとする。
[信号分析装置の構成及び処理]
図5及び図6を用いて、本第7の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図5は、本第7の実施形態に係る信号分析装置の構成の一例を示す図である。図6は、本第7の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
図5に示すように、信号分析装置301は、周波数領域変換部310、周波数領域信号分離部320、時間領域変換部330からなる。
まず、信号分析装置301の各部の概要について説明する。周波数領域変換部310は、入力された観測信号yi(t)を受け取り(ステップS40)、短時間フーリエ変換などにより時間周波数領域の観測信号yi(k,f)を計算する(ステップS41)。ここで、kは時間フレームのインデックスであり、fは周波数ビンのインデックスである。yi(k,f)をyi f(k)とも表す。
周波数領域信号分離部320は、fごとに信号分離を行って(ステップS42)、時間周波数領域の源信号成分xj(k,f)の推定値^xj(k,f)を得る。^xj(k,f)を^xj f(k)とも表す。
時間領域変換部330は、時間周波数領域の源信号成分xj(k,f)の推定値^xj(k,f)に基づいて、逆短時間フーリエ変換などにより、時間領域の源信号成分xj(t)の推定値^xj(t)を計算し(ステップS43)、出力する。
次に、周波数領域信号分離部320についてより詳細に説明する。図19及び図20を用いて、周波数領域信号分離部320の構成及び周波数領域信号分離処理の処理手順について説明する。図19は、周波数領域信号分離部320の構成の一例を示す図である。図20は、周波数領域信号分離処理の処理手順の一例を示すフローチャートである。
図19に示すように、周波数領域信号分離部320は、初期化部(図示しない)、観測信号ベクトル作成部710、信号パラメータ分析部720、同時無相関化部730、信号存在事後確率更新部740、シングルチャネルウィーナーフィルタ部750、事後共分散行列更新部760、混合重み更新部770、分散パラメータ更新部780、空間共分散行列更新部790、収束判定部(図示しない)及び源信号成分推定部800を有する。これら周波数領域信号分離部320の各部での処理は、周波数ビン毎に行う。
周波数領域信号分離部320は、f←0により周波数ビンのインデックスfを初期化する(ステップS101)。周波数領域信号分離部320は、f←f+1により周波数ビンのインデックスfを1増やす(ステップS102)。
図示しない初期化部は、混合重みαB(f)の初期推定値^αB (0)(f)と、分散パラメータvj(k,f)の初期推定値^vj (0)(k,f)と、空間共分散行列Rj(f)の初期推定値^Rj (0)(f)と、を設定する(ステップS103)。これらは、例えば乱数に基づいて設定すればよい。
観測信号ベクトル作成部710は、周波数領域変換部からの時間周波数領域の観測信号yi(k,f)を受け取り(ステップS104)、観測信号ベクトルy(k,f)を(316)式により作成する(ステップS105)。
ここで、本第7の実施形態における観測信号ベクトルy(k,f)のモデル化について説明する。本実施形態では、各時間周波数点(k,f)において観測信号ベクトルy(k,f)は高々2個の源信号成分しか含まないと仮定する。各時間周波数点(k,f)において観測信号ベクトルy(k,f)に含まれる源信号成分は、そのような源信号成分のインデックスの全体である集合B(k,f)により表すことができる(ただし、集合B(k,f)の元の個数#B(k,f)は#B(k,f)≦2を満たすものとする)。このとき、観測信号ベクトルy(k,f)は、B(k,f)のすべての元jに対する源信号成分xj(k,f)の和として、(317)式によりモデル化される。
例えば、ある時間周波数点(k,f)でB(k,f)={1,3}である場合、(k,f)では、(317)式は、(318)式のようになる。
B(k,f)が取りうる値(集合)の全体を、集合{1,・・・,L}の空でない部分集合からなる集合系(すなわち、冪集合2{1,・・・,L}の部分集合であって空集合を元として含まないもの)Hで表す。例えば、B(k,f)が#B≦2なる任意の空でない集合B⊂{1,・・・,L}になり得る場合、Hは、(319)式及び(320)式で与えられる。
(319)式及び(320)式の集合系HをHmaxで表す。本第7の実施形態では、集合系Hは、あらかじめ与えられており、H⊂Hmaxであるとする。また、#B=2なる集合B∈Hが少なくとも一つはあるものとする。
すなわち、集合系Hの元(集合)Bが存在して、#B=2であるとする。
信号パラメータ分析部720は、空間共分散行列更新部790からの空間共分散行列Rj(f)の推定値^Rj(f)(ただし、例外として、信号パラメータ分析部720での最初の処理の際には、初期化部からの空間共分散行列Rj(f)の初期推定値^Rj (0)(f))を受け取って、L個の源信号成分のうちの任意の2個xj(k,f),xl(k,f)(j,lは1以上L以下の相異なる任意の自然数)に対して、同時無相関化行列P{j,l}(f)と対角行列Λj,{j,l}(f)を計算する(ステップS106)。ここで、同時無相関化行列P{j,l}(f)は2個の源信号成分xj(k,f),xl(k,f)を同時無相関化する行列(すなわち、2個の源信号成分xj(k,f),xl(k,f)の空間共分散行列を同時対角化する行列)であり、添え字である集合{j,l}は、同時無相関化行列P{j,l}(f)によって同時無相関化される源信号成分xj(k,f),xl(k,f)のインデックスを表している。また、対角行列Λj,{j,l}(f)は、源信号成分xj(k,f)の空間共分散行列を同時無相関化行列P{j,l}(f)によって対角化した際の対角行列である。
具体的には、信号パラメータ分析部720は以下の処理を行う。まず、2個の源信号成分xj(k,f),xl(k,f)(j,lは、1≦j<l≦Lなる任意の2つの自然数)の空間共分散行列の推定値に対する一般化固有値問題^Rl(f)p=λ^Rj(f)pを解くことにより、(ただし、例外として、信号パラメータ分析部720での最初の処理の際には、これらの源信号成分の空間共分散行列の初期推定値に対する一般化固有値問題^Rl (0)(f)p=λ^Rj (0)(f)pを解くことにより、)一般化固有値λ1,{j,l}(f),・・・,λI,{j,l}(f)と、一般化固有値λ1,{j,l}(f),・・・,λI,{j,l}(f)にそれぞれ対応する一般化固有ベクトルp1,{j,l}(f),・・・,pI,{j,l}(f)と、を求める。ただし、一般化固有ベクトルp1,{j,l}(f),・・・,pI,{j,l}(f)は、(321)式を満たすように取るものとする。
なお、信号パラメータ分析部720での最初の処理の際には、(321)式において、^Rj(f)を^Rj (0)(f)で置き換える。信号パラメータ分析部720は、次に、一般化固有ベクトルp1,{j,l}(f),・・・,pI,{j,l}(f)を列ベクトルとする行列を同時無相関化行列P{j,l}(f)とし、I次単位行列を対角行列Λj,{j,l}(f)とし、一般化固有値λ1,{j,l}(f),・・・,λI,{j,l}(f)を対角要素とする対角行列を対角行列Λl,{j,l}(f)とする。
同時無相関化部730は、観測信号ベクトル作成部710からの観測信号ベクトルy(k,f)と、信号パラメータ分析部720からの同時無相関化行列P{j,l}(f)と、を受け取って、源信号成分xj(k,f),xl(k,f)(j,lは、1≦j<l≦Lなる任意の2つの自然数)について同時無相関化された観測信号ベクトルy´{j,l}(k,f)を次の(322)式により計算する(ステップS107)。
ここで、添え字である集合{j,l}は、同時無相関化される源信号成分xj(k,f),xl(k,f)のインデックスを表す。
信号存在事後確率更新部740は、信号パラメータ分析部720からの同時無相関化行列P{j,l}(f)及び対角行列Λj,{j,l}(f)と、同時無相関化部730からの同時無相関化された観測信号ベクトルy´{j,l}(k,f)と、混合重み更新部770からの混合重みαB(f)の推定値^αB(f)(ただし、例外として、信号存在事後確率更新部740における最初の処理の際には、初期化部からの混合重み^αB(f)の初期推定値^αB (0)(f))と、分散パラメータ更新部780からの分散パラメータvj(k,f)の推定値^vj(k,f)(ただし、例外として、信号存在事後確率更新部740における最初の処理の際には、初期化部からの分散パラメータvj(k,f)の初期推定値^vj (0)(k,f))と、を受け取って、集合系Hの各元(集合)Bに対して、信号存在事後確率γB(k,f)を(323)式により更新する(ステップS108)。
なお、集合系Hの定義によっては、(323)式において、ただ一つの元からなる集合(一元集合){j}に対する同時無相関化行列P{j}(f)、対角行列Λj,{j}(f)及び同時無相関化された観測信号ベクトルy´{j}(k,f)が現れることがあるが、これらは二元集合に対するP{j,l}(f)、Λj,{j,l}(f)及びy´{j,l}(k,f)を用いて、(324)式、(325)式及び(326)式により定義される。
ここで、ρ:{1,...,L}→{1,...,L}は、「すべてのj∈{1,...,L}に対してρ(j)≠j」なる所定の関数である。このように定義する理由は、次の通りである。
すなわち、1個の源信号成分xj(k,f)の同時無相関化とは、単なるxj(k,f)の無相関化に他ならず、それは2個の源信号成分xj(k,f),xρ(j)(k,f)の同時無相関化により当然実現されるからである。
シングルチャネルウィーナーフィルタ部750は、同時無相関化部730からの同時無相関化された観測信号ベクトルy´{j,l}(k,f)と、信号パラメータ分析部720からの対角行列Λj,{j,l}(f)と、分散パラメータ更新部780からの分散パラメータvj(k,f)の推定値^vj(k,f)(ただし、例外として、シングルチャネルウィーナーフィルタ部750における最初の処理の際には、初期化部からの分散パラメータvj(k,f)の初期推定値^vj (0)(k,f))と、を受け取って、同時無相関化された事後平均m´j,{j,l}(k,f)(j,lは1以上L以下の相異なる任意の自然数)を(327)式により更新する(ステップS109)。
事後共分散行列更新部760は、信号パラメータ分析部720からの対角行列Λj,{j,l}(f)と、分散パラメータ更新部780からの分散パラメータvj(k,f)の推定値^vj(k,f)(ただし、例外として、事後共分散行列更新部760における最初の処理の際には、初期化部からの分散パラメータvj(k,f)の初期推定値^vj (0)(k,f))と、を受け取って、同時無相関化された事後共分散行列Σ´{j,l}(k,f)(j,lは1≦j<l≦Lなる任意の2つの自然数)を次の(328)式により更新する(ステップS110)。
混合重み更新部770は、信号存在事後確率更新部740からの信号存在事後確率γB(k,f)を受け取って、(329)式により混合重みαB(f)の推定値^αB(f)(Bは集合系Hの任意の元(集合))を更新する(ステップS111)。
分散パラメータ更新部780は、信号パラメータ分析部720からの対角行列Λj,{j,l}(f)と、信号存在事後確率更新部740からの信号存在事後確率γB(k,f)と、シングルチャネルウィーナーフィルタ部750からの同時無相関化された事後平均m´j,{j,l}(k,f)と、事後共分散行列更新部760からの同時無相関化された事後共分散行列Σ´{j,l}(k,f)とを受け取って、分散パラメータvj(k,f)の推定値^vj(k,f)を(330)式により更新する(ステップS112)。
なお、集合系Hの定義によっては、(330)式において、一元集合{j}に対する同時無相関化された事後平均m´j,{j}(k,f)及び同時無相関化された事後共分散行列Σ´{j}(k,f)が現れることがあるが、これらは(331)式及び(332)式により定義される。
空間共分散行列更新部790は、信号パラメータ分析部720からの同時無相関化行列P{j,l}(f)と、信号存在事後確率更新部740からの信号存在事後確率γB(k,f)と、シングルチャネルウィーナーフィルタ部750からの同時無相関化された事後平均m´j,{j,l}(k,f)と、事後共分散行列更新部760からの同時無相関化された事後共分散行列Σ´{j,l}(k,f)と、分散パラメータ更新部780からの分散パラメータvj(k,f)の推定値^vj(k,f)と、を受け取って、空間共分散行列Rj(f)の推定値^Rj(f)を(333)式により更新する(ステップS113)。
図示しない収束判定部は、収束したかどうかの判定を行う(ステップS114)。収束判定部が収束していないと判定した場合(ステップS114:No)、信号パラメータ分析部720での処理に戻って、処理が継続される。
一方、収束判定部が収束したと判定した場合(ステップS114:Yes)、源信号成分推定部800での処理が行われる。源信号成分推定部800は、信号パラメータ分析部720からの同時無相関化行列P{j,l}(f)と、信号存在事後確率更新部740からの信号存在事後確率γB(k,f)と、シングルチャネルウィーナーフィルタ部750からの同時無相関化された事後平均m´j,{j,l}(k,f)と、を受け取って、源信号成分xj(k,f)の推定値^xj(k,f)を(334)式により推定し(ステップS115)、出力する。
周波数領域信号分離部320は、f=Fであるかどうか判定する(ステップS116)。f=Fでないと判定された場合(ステップS116:No)、周波数領域信号分離部320でのf←f+1の処理に戻って、処理が継続される。
一方、f=Fであると判定された場合(ステップS116:Yes)、周波数領域信号分離部320での処理は終了する。
[信号分析装置の処理の背景]
次に、本第7の実施形態に係る信号分析装置の処理の背景について説明する。本第7の実施形態では、観測信号ベクトルy(k,f)は、(317)式でモデル化される。
本第7の実施形態では、源信号成分xj(k,f)が平均0、共分散行列Sj(k,f)=vj(k,f)Rj(f)の複素ガウス分布に従うと仮定する。このとき、B(k,f)が値B∈Hを取るという条件下での観測信号ベクトルy(k,f)の確率分布は、次の(335)式で与えられる。
ここで、Φ1はすべてのモデルパラメータをまとめて表したものである。また、集合B(k,f)の確率分布は次の(336)式で与えられる。
ただし、αB(f)は、次の(337)式の制約条件を満たす。
(338)式及び(336)式より、観測信号ベクトルy(k,f)の確率分布は、次の(338)式で与えられる。
以下、αB(f)を混合重みと呼ぶ。パラメータΦ1は、混合重みαB(f)、分散パラメータvj(k,f)及び空間共分散行列Rj(f)からなる。パラメータΦ1は、次の(339)式の尤度関数の最大化により推定できる。
参考文献6「R. Ikeshita, M. Togami, Y. Kawaguchi, Y. Fujita, and K. Nagamatsu, “Local Gaussian model with source-set constraints in audio source separation”, IEEE International Workshop on Machine Learning for Signal Processing, Sept. 2017.」には、EMアルゴリズムに基づいて上式の尤度関数を最大化することにより、パラメータΦ1を推定する方法が開示されている。この従来技術では、EステップとMステップを交互に反復することにより、パラメータΦ1を推定する。
Eステップでは、信号存在事後確率γB(k,f)、事後平均mj,B(k,f)及び事後共分散行列ΣB(k,f)を次の(340)式〜(344)式により計算する。
ここで、信号存在事後確率γB(k,f)は、観測信号ベクトルy(k,f)が与えられたという条件下での集合B(k,f)の条件付き確率分布である事後分布として、次の(345)式で定義される。
ここで、^Φ1はパラメータΦ1の現在の推定値を表す。また、事後平均mj,B(k,f)及び事後共分散行列ΣB(k,f)とは、観測信号ベクトルy(k,f)と集合B(k,f)が与えられたという条件下での源信号成分xj(k,f)の条件付き確率分布である事後分布(複素ガウス分布であることが示される)の平均及び共分散行列として定義される。この事後分布は次の(346)式のように表される。
Mステップでは、混合重みαB(f)の推定値^αB(f)、分散パラメータvj(k,f)の推定値^vj(k,f)及び空間共分散行列Rj(f)の推定値^Rj(f)を、次の(347)式〜(349)式の更新則により更新する。
しかしながら、この従来技術では、(340)式の信号存在事後確率γB(k,f)の更新、(342)式の事後平均mj,B(k,f)の更新、及び(344)式の事後共分散行列Σ{j,l}(k,f)の更新において、行列式計算、逆行列計算及び行列乗算を時間周波数点ごとに行う必要があった。具体的には、(340)式における(350)式の形の行列式計算と、(340)式及び(342)式及び(344)式における(351)式の形の逆行列計算と、(344)式における(352)式の形の行列乗算と、を時間周波数点(k,f)ごとに行う必要があった。
行列式計算、逆行列計算及び行列乗算はいずれもO(I3)のオーダーの計算量を必要とする行列演算であり、かつ時間周波数点の数は通常、非常に大きいため、この従来技術では膨大な計算量を要していた。
いま仮に、(350)式〜(352)式における共分散行列^vj(k,f)^Rj(f)及び^vl(k,f)^Rl(f)がともに対角行列であるとすると、(350)式〜(352)式における行列式計算、逆行列計算及び行列乗算は単なる対角要素に対するスカラー演算に帰着するため、少ない計算量で計算できる。しかしながら、通常、源信号成分(ベクトル)の異なる要素は互いに相関を持つため、共分散行列^vj(k,f)^Rj(f)及び^vl(k,f)^Rl(f)は対角行列ではない。
そこで、本第7の実施形態では、一般化固有値問題に基づいて2つの空間共分散行列^Rj(f),^Rl(f)を同時対角化することを考える。λ1,{j,l}(f),・・・,λI,{j,l}(f)を空間共分散行列^Rl(f),^Rj(f)に対する一般化固有値問題((353)式)の一般化固有値とし、p1,{j,l}(f),・・・,pI,{j,l}(f)を一般化固有値λ1,{j,l}(f),・・・,λI,{j,l}(f)に対応する一般化固有ベクトルとする。
ただし、一般化固有ベクトルp1,{j,l}(f),・・・,pI,{j,l}(f)は、次の(354)式を満たすように取る。
p{j,l}(f)を一般化固有ベクトルp1,{j,l}(f),・・・,pI,{j,l}(f)を列ベクトルとする行列、Λ{j,l}(f)を一般化固有値λ1,{j,l}(f),・・・,λI,{j,l}(f)を対角要素とする対角行列とすると、次の(355)式が成り立つ。
(355)式は、行列P{j,l}(f)が2つの空間共分散行列^Rj(f),^Rl(f)を同時対角化する(2つの源信号成分xj(k,f),xl(k,f)を同時無相関化する)ことを示している。そこで、行列P{j,l}(f)を同時無相関化行列と呼ぶ。
(355)式を空間共分散行列^Rj(f),^Rl(f)について解いて、(340)式、(342)式及び(344)式に代入し、(356)式〜(358)式と定義する(j,lは1以上L以下の互いに異なる任意の整数)と、(362)、(363)式及び(364)式(後述)を得る。
(356)式及び(357)式を事後平均mj,{j,l}(k,f)及び事後共分散行列Σ{j,l}(k,f)について解き、(348)式及び(349)式に代入すると(366)式及び(368)式を得る。
これで、本第7の実施形態に係る信号分析装置における更新則が導出された。
このように、本第7の実施形態に係る信号分析装置では、信号パラメータ分析部720は、2個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列^Rl(f),^Rj(f)に対するpi,{j,l}(f)H^Rj(f)pm,{j,l}(f)=0(i,mは1以上I以下の互いに異なる任意の整数)を満たす、それぞれ異なるI個(Iは2以上の整数)の一般化固有ベクトルp1,{j,l}(f),・・・,pI,{j,l}(f)を列ベクトルとする行列である同時無相関化行列P{j,l}(f)またはそのエルミート転置P{j,l}(f)Hを、I個の位置で取得される観測信号による観測信号ベクトルについて2個の源信号に対応する成分を無相関化するためのパラメータとして得る。
[第7の実施形態の実施例]
本実施例では、集合系Hの具体例について説明する。本実施例では、集合系Hを(359)式とする。
これは、各時間周波数点において観測信号ベクトルは、1番目の源信号成分をつねに含み、2番目〜J番目の源信号成分のうちのちょうど1つを含むことを意味する。例えば、1番目の源信号成分は背景雑音のようなスパースでない源信号成分をモデル化し、2番目〜J番目の源信号成分は音声のようなスパースな源信号成分をモデル化する。
信号パラメータ分析部720、同時無相関化部730、信号存在事後確率更新部740、シングルチャネルウィーナーフィルタ部750、事後共分散行列更新部760、混合重み更新部770、分散パラメータ更新部780、空間共分散行列更新部790及び源信号成分推定部800における処理の更新式を(360)式〜(370)式に示す。本第7の実施形態では、Hの元である集合{1,2},{1,3},・・・,{1,J}は、それぞれインデックス2,3,・・・,Jと1対1に対応しているから、簡単のため以下の記法を用いる。すなわち、y´{1,j}(k,f)をy´j(k,f)と書き、P{1,j}(f)をPj(f)と書き、γ{1,j}(k,f)をγj(k,f)と書き、^α{1,j}(f)を^αj(f)と書き、Λl,{1,j}(f)をΛlj(f)と書き、m´l,{1,j}(k,f)をm´lj(k,f)と書き、Σ´{1,j}(k,f)をΣ´j(k,f)と書く。ただし、j=2,・・・,Jとし、l=1,jとする(下記の更新式でも同様)。なお、本実施例では、信号パラメータ分析部720は、L個の源信号成分のうちの任意の2個xj(k,f),xl(k,f)(j,lは1以上L以下の相異なる任意の自然数)に対して、同時無相関化行列P{j,l}(f)と対角行列Λj,{j,l}(f)を計算する必要はなく、(360)式にあるように、x1(k,f),xj(k,f)(jは2以上L以下の任意の自然数)に対して、同時無相関化行列P{1,j}(f)(すなわちPj(f))と対角行列Λl,{1,j}(f)(すなわちΛlj(f))を計算すれば十分である。
[第8の実施形態]
本第8の実施形態では、第3の実施形態に係る信号分析装置の他の構成例について説明する。図21は、第8の実施形態に係る信号分析装置の構成の一例を示す図である。図22は、第8の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。図21に示すように、第8の実施の形態に係る信号分析装置801は、周波数領域変換部310、周波数領域信号分離部320−1及び時間領域変換部330を有する。
信号分析装置801では、周波数領域信号分離部320−1において、スペクトル基底tn(f)と、アクティベーション基底wn(k)と、分配関数znjと、同時無相関化された空間共分散行列Λj(f)と、同時無相関化行列P(f)との推定値とを計算する。
スペクトル基底tn(f)と、アクティベーション基底wn(k)と、分配関数znjとは、分散パラメータvj(k,f)をモデル化するパラメータである。分散パラメータvj(k,f)は、これらのパラメータを用いて(371)式によりモデル化される。
ここで、nは、スペクトル基底およびアクティベーション基底のインデックスである。Nは、スペクトル基底およびアクティベーション基底の個数である。
周波数領域信号分離部320−1は、初期化部(図示しない)、観測信号ベクトル作成部810、同時無相関化830、シングルチャネルウィーナーフィルタ部840、観測共分散行列更新部850、分散パラメータ更新部860、空間共分散行列更新部870、信号パラメータ分析部820、収束判定部(図示しない)及び同時無相関化逆変換部880を有する。
[信号分析装置の構成及び処理]
図21及び図22を用いて、第8の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。
図示しない初期化部は、分散パラメータvj(k,f)と、スペクトル基底tn(f)と、アクティベーション基底wn(k)と、分配関数znjと、同時無相関化された空間共分散行列Λj(f)と、同時無相関化行列P(f)と、の初期推定値を計算する初期化を行う(ステップS121)。初期推定値は、上付きの(0)を付して、^vj (0)(k,f)などと表す。
周波数領域変換部310は、観測信号を取得し(ステップS122)、観測信号を周波数領域に変換する(ステップS123)ことにより、時間周波数領域の観測信号yi(k,f)を計算する。
観測信号ベクトル作成部810は、観測信号ベクトルy(k,f)を第7の実施形態と同様にして作成する(ステップS124)。
同時無相関化部830は、同時無相関化された観測信号ベクトルy´(k,f)を(372)式により計算し、更新する(ステップS125)。
ここで、^P(f)は信号パラメータ分析部820からの同時無相関化行列の推定値である。
観測共分散行列更新部850は、同時無相関化された観測共分散行列U´(k,f)を(373)式により更新する(ステップS126)。
ここで、^vj(k,f)は、分散パラメータ更新部860からの分散パラメータの推定値である。^Λj(f)は、空間共分散行列更新部870からの同時無相関化された空間共分散行列の推定値である。
分散パラメータ更新部860は、次の手順により分散パラメータの推定値^vj(k,f)を更新する(ステップS127)。
まず、分散パラメータ更新部860は、n番目の基底に対応する同時無相関化された空間共分散行列Σn(f)を(374)式により更新する。(374)式において、^znjは、分配関数の推定値である。
次に、分散パラメータ更新部はスペクトル基底の推定値^tn(f)を(375)式により更新する。(375)式において、^wn(k)はアクティベーション基底の推定値である。
次に、分散パラメータ更新部860は、アクティベーション基底の推定値^wn(k)を(376)式により更新する。
次に、分散パラメータ更新部860は、分配関数の推定値^znjを(377)式により更新する。
最後に、分散パラメータ更新部860は、分散パラメータの推定値^vj(k,f)を(378)式により更新する。
空間共分散行列更新部870は、(379)式により、同時無相関化された空間共分散行列の推定値^Λj(f)を更新する(ステップS128)。
ここで、#は、行列の幾何平均であり、Aj(f)およびBj(f)は、(380)式および(381)式で計算される行列である。
信号パラメータ分析部820は、同時無相関化行列の推定値^P(f)を(382)式により更新する(ステップS129)。
図示しない収束判定部は、収束したかどうかの判定を行う(ステップS130)。収束判定部が収束していないと判定した場合(ステップS130:No)、同時無相関化部830での処理(ステップS125)に戻って、処理が継続される。
一方、収束判定部が収束したと判定した場合(ステップS130:Yes)、シングルチャネルウィーナーフィルタ部840は、同時無相関化された源信号成分の推定値^x´j(k,f)を(383)式により計算する(ステップS131)。
シングルチャネルウィーナーフィルタ部840は、同時無相関化された源信号成分の推定値^x´j(k,f)を更新する(ステップS132)。
同時無相関化逆変換部880は、源信号成分の推定値^xj(k,f)を(384)式により計算する(ステップS133)。
時間領域変換部330は、第7の実施形態と同様の処理を行って、時間領域の源信号成分の推定値^xj(t)を計算し、出力する(ステップS134)。
[第8の実施形態の効果]
第3の実施形態では、音の方向性のみを基に信号処理を行っていた。これに対して、本第8の実施形態では、音の構造(パーツ)も考慮して信号分析を行うため、音源分離の精度を向上させることができる。
そして、第8の実施形態では、パーミテーション問題が原理的に起こらないという効果を奏する。本第8の実施形態では、全ての周波数において、どのようなパーツがあるかを考慮した処理を行うため、周波数ごとに独立に処理しないことから、パーミテーション問題が原理的に起こらない。
[システム構成等]
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行われる各処理機能は、その全部又は任意の一部が、CPU及び当該CPUにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
また、本実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。すなわち、上記学習方法及び音声認識方法において説明した処理は、記載の順にしたがって時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。
[プログラム]
図23は、プログラムが実行されることにより、信号分析装置1,201,301,401,501,601が実現されるコンピュータの一例を示す図である。コンピュータ1000は、例えば、メモリ1010、CPU1020を有する。また、コンピュータ1000は、ハードディスクドライブインタフェース1030、ディスクドライブインタフェース1040、シリアルポートインタフェース1050、ビデオアダプタ1060、ネットワークインタフェース1070を有する。これらの各部は、バス1080によって接続される。
メモリ1010は、ROM1011及びRAM1012を含む。ROM1011は、例えば、BIOS(Basic Input Output System)等のブートプログラムを記憶する。ハードディスクドライブインタフェース1030は、ハードディスクドライブ1090に接続される。ディスクドライブインタフェース1040は、ディスクドライブ1100に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ1100に挿入される。シリアルポートインタフェース1050は、例えばマウス1110、キーボード1120に接続される。ビデオアダプタ1060は、例えばディスプレイ1130に接続される。
ハードディスクドライブ1090は、例えば、OS(Operating System)1091、アプリケーションプログラム1092、プログラムモジュール1093、プログラムデータ1094を記憶する。すなわち、信号分析装置1,201,301,401,501,601の各処理を規定するプログラムは、コンピュータ1000により実行可能なコードが記述されたプログラムモジュール1093として実装される。プログラムモジュール1093は、例えばハードディスクドライブ1090に記憶される。例えば、信号分析装置1,201,301,401,501,601における機能構成と同様の処理を実行するためのプログラムモジュール1093が、ハードディスクドライブ1090に記憶される。なお、ハードディスクドライブ1090は、SSD(Solid State Drive)により代替されてもよい。
また、上述した実施形態の処理で用いられる設定データは、プログラムデータ1094として、例えばメモリ1010やハードディスクドライブ1090に記憶される。そして、CPU1020が、メモリ1010やハードディスクドライブ1090に記憶されたプログラムモジュール1093やプログラムデータ1094を必要に応じてRAM1012に読み出して実行する。
なお、プログラムモジュール1093やプログラムデータ1094は、ハードディスクドライブ1090に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ1100等を介してCPU1020によって読み出されてもよい。あるいは、プログラムモジュール1093及びプログラムデータ1094は、ネットワーク(LAN(Local Area Network)、WAN(Wide Area Network)等)を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール1093及びプログラムデータ1094は、他のコンピュータから、ネットワークインタフェース1070を介してCPU1020によって読み出されてもよい。
以上、本発明者によってなされた発明を適用した実施形態について説明したが、本実施形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。すなわち、本実施形態に基づいて当業者等によりなされる他の実施形態、実施例及び運用技術等はすべて本発明の範疇に含まれる。