JP6106611B2 - モデル推定装置、雑音抑圧装置、音声強調装置、これらの方法及びプログラム - Google Patents

Description

本発明は、複数のマイクロホンで観測された、拡散性雑音を含む観測信号を用いて、拡散性雑音を含まない信号であるノイズフリー信号を推定する雑音抑圧技術、及び、雑音及び残響を含まないクリーン信号を推定する音声強調技術に関する。さらに、観測信号を、拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとして、雑音抑圧時または音声強調時等に用いるモデルパラメータの推定技術に関する。

＜拡散性雑音下における雑音抑圧及び音声強調＞
拡散性雑音とは、多数の雑音源、または空間的な広がりをもつ雑音源により生じる、様々な方向から到来する雑音を指す。多数の雑音源による拡散性雑音の例としては、食堂での人々の話し声や食器の音があり、空間的な広がりをもつ雑音源による拡散性雑音の例としては、電車内での車体の振動による雑音がある。このような拡散性雑音は、ビームフォーミングのような、従来の雑音抑圧技術による抑圧が困難であり、雑音抑圧技術の適用範囲を大幅に制限する要因となっていた。また、拡散性雑音は、少数の点音源から生じる雑音と比べて、モデル化がより難しく、実環境における音声強調技術の応用範囲を大きく制限する要因となっていた。

観測信号を、拡散性雑音と拡散性雑音を含まない信号との和でモデル化する場合、モデル化したときのモデルパラメータとして、拡散性雑音の共分散行列Φ^v _t及び音源信号の共分散行列Φ^x _tを含み、共分散行列Φ^v _t及びΦ^x _tの推定精度を上げることで、雑音抑圧性能及び音声強調性能を向上させることができる。

信号と拡散性雑音の空間的な性質を利用して、共分散行列Φ^v _t及びΦ^x _tを推定する従来技術として非特許文献１が知られている。

非特許文献１では、次式で定義されるコスト関数を最小化することにより、共分散行列Φ^v _t及びΦ^x _tを推定する。

ここで、D_EU(Φ^y _t,Φ^x _t+Φ^v _t):=||Φ^y _t-(Φ^x _t+Φ^v _t)||_F ²は、観測信号を用いて得られる観測信号の共分散行列Φ^y _tと共分散行列のモデルΦ^x _t+Φ^v _tとのユークリッド距離の平方であり、||・||_Fはフロベニウスノルムである。なお、A:=Bは、「Aの意味するところをBと定義する」ことを意味し、式中「:=」を

とも記載する。

N. Ito, E. Vincent, N. Ono, and S. Sagayama, "General algorithms for estimating spectrogram and transfer functions of target signal for blind suppression of diffuse noise," in Proc. IEEE International Workshop on Machine Learning for Signal Processing (MLSP), 2013.

上記のように、非特許文献１では、簡単のため、式(a1)のように、ユークリッド距離の平方D_EUという単純な尺度により、Φ^y _tとΦ^x _t+Φ^v _tとの間の距離を測っている。しかしながら、より高い音声品質、ならびに、より高精度な音声認識の実現のためには、このような尺度は必ずしも最適ではない。なぜなら、音声品質の向上であれ、音声認識精度の向上であれ、人間の聴覚も、音声認識の特徴量(例えばMFCC: mel - frequency cepstral coefficients)も、ともに対数的であるため、ユークリッド距離の平方ではなく、対数パワー領域での誤差を小さくする必要があるからである。実際、例えば、振幅が大きい信号に多少の歪みが生じても、聴感上さほど気にならないのに対し、振幅が小さい信号に少しでも歪みが生じると、敏感に知覚される。

本発明は、このような背景に鑑みてなされたものであり、音響信号により適した距離尺度に基づく、拡散性雑音の存在下での雑音抑圧技術、音声強調技術、雑音抑圧時及び音声強調時等に用いるモデルパラメータの推定技術を提供することを目的とする。

上記の課題を解決するために、本発明の第一の態様によれば、観測信号を、非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、モデル化するときのモデルパラメータは拡散性雑音の共分散行列を特定するモデルパラメータと、ノイズフリー信号の共分散行列を特定するモデルパラメータとを含み、モデル推定装置は、観測信号から得られる共分散行列と、拡散性雑音の共分散行列とノイズフリー信号の共分散行列との和との板倉-齋藤距離が小さくなるように、モデルパラメータを更新するモデル推定部を含む。

上記の課題を解決するために、本発明の他の態様によれば、観測信号を、非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなるブロックの番号をt∈{1,…,T}とし、ブロック内のフレームの番号をu∈{1,…,U}とし、ブロックtごとの観測信号y_tuの共分散行列をΦ^y _tとし、ブロックtごとのノイズフリー信号x_tuの共分散行列をΦ^x _tとし、ブロックtごとの拡散性雑音v_tuの共分散行列をΦ^v _tとし、モデル推定装置は、観測信号y_tから得られる共分散行列Φ^y _tと、Φ^x _t＋Φ^v _tとの板倉-齋藤距離の、すべてのブロックtについての総和が小さくなるように、共分散行列Φ^x _tと共分散行列Φ^v _tとを更新するモデル推定部を含む。ノイズフリー信号は音源信号である。

上記の課題を解決するために、本発明の他の態様によれば、フレームごとの観測信号を非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなる時間区間をブロックとして、モデル推定装置は、複数ブロック分の観測信号から得られる共分散行列と、複数ブロック分の拡散性雑音の共分散行列と複数ブロック分のノイズフリー信号の共分散行列との和と、の板倉-齋藤距離が小さくなるように、複数ブロック分の拡散性雑音の共分散行列を特定するパラメータと複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータとを更新するモデル推定部を含む。ノイズフリー信号は、音源から発せられる直接音と初期反射成分とからなる信号であるクリーン信号に残響が重畳された信号である。複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータは、フレームごとのノイズフリー信号に含まれる残響をK次の自己回帰過程によりモデル化したときの自己回帰過程の予測係数行列G_k(k=1,2,…,K)と、フレームごとのクリーン信号の共分散行列と、を含む。

上記の課題を解決するために、本発明の他の態様によれば、観測信号を、非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、モデル化するときのモデルパラメータは拡散性雑音の共分散行列を特定するモデルパラメータと、ノイズフリー信号の共分散行列を特定するモデルパラメータとを含み、モデル推定方法は、観測信号から得られる共分散行列と、拡散性雑音の共分散行列とノイズフリー信号の共分散行列との和との板倉-齋藤距離が小さくなるように、モデルパラメータを更新するモデル推定ステップを含む。

上記の課題を解決するために、本発明の他の態様によれば、観測信号を、非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなるブロックの番号をt∈{1,…,T}とし、ブロック内のフレームの番号をu∈{1,…,U}とし、ブロックtごとの観測信号y_tuの共分散行列をΦ^y _tとし、ブロックtごとのノイズフリー信号x_tuの共分散行列をΦ^x _tとし、ブロックtごとの拡散性雑音v_tuの共分散行列をΦ^v _tとし、モデル推定方法は、観測信号y_tから得られる共分散行列Φ^y _tと、Φ^x _t＋Φ^v _tとの板倉-齋藤距離の、すべてのブロックtについての総和が小さくなるように、共分散行列Φ^x _tと共分散行列Φ^v _tとを更新するモデル推定ステップを含む。ノイズフリー信号は音源信号である。

上記の課題を解決するために、本発明の他の態様によれば、フレームごとの観測信号を非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなる時間区間をブロックとして、モデル推定方法は、複数ブロック分の観測信号から得られる共分散行列と、複数ブロック分の拡散性雑音の共分散行列と複数ブロック分のノイズフリー信号の共分散行列との和と、の板倉-齋藤距離が小さくなるように、複数ブロック分の拡散性雑音の共分散行列を特定するパラメータと複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータとを更新するモデル推定ステップを含む。ノイズフリー信号は、音源から発せられる直接音と初期反射成分とからなる信号であるクリーン信号に残響が重畳された信号である。複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータは、フレームごとのノイズフリー信号に含まれる残響をK次の自己回帰過程によりモデル化したときの自己回帰過程の予測係数行列G_k(k=1,2,…,K)と、フレームごとのクリーン信号の共分散行列と、を含む。

本発明によれば、従来技術よりも推定精度よく雑音抑圧時及び音声強調時等に用いるモデルパラメータの推定することができるという効果を奏する。

第一実施形態に係るモデル推定装置の機能ブロック図。 第一実施形態に係るモデル推定装置の処理フローの例を示す図。 第二実施形態に係る雑音抑圧装置の機能ブロック図。 第二実施形態に係る雑音抑圧装置の処理フローの例を示す図。 第二実施形態に係る雑音抑圧装置の効果を説明するための図。 第三実施形態に係るモデル推定装置の機能ブロック図。 第三実施形態に係るモデル推定装置の処理フローの例を示す図。 第四実施形態に係る音声強調装置の機能ブロック図。 第四実施形態に係る音声強調装置の処理フローの例を示す図。 第四実施形態に係る音声強調装置の効果を説明するための図。 評価指標による客観評価の結果を示す図。 第五実施形態に係る雑音抑圧装置の機能ブロック図。 第五実施形態に係る雑音抑圧装置の処理フローの例を示す図である。

以下、本発明の実施形態について説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号「^」、「~」等は、本来直後の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直前に記載する。式中においてはこれらの記号は本来の位置に記述している。

＜第一実施形態＞
本実施形態では、観測信号から拡散性雑音抑圧時に用いるモデルパラメータを推定する。なお、拡散性雑音抑圧とは，上述の拡散性雑音を抑圧する技術である。

＜拡散性雑音抑圧の問題設定＞
本実施形態では、特に断りのない限り、各信号を短時間フーリエ変換（short-time Fouriertransform:STFT）などの時間周波数領域で表現する。フレームの総数をIで表し、フレームの番号をi∈{1,…,I}で表す。また、本実施形態では、簡潔さのため、周波数ビンの番号の表示を省略するが、この省略が混乱をもたらすおそれは小さい。本実施形態では、周波数ビン毎に独立に処理を行うからである。

M個のマイクロホンによる観測信号をy_i∈C^M、音源信号をx_i∈C^M、拡散性雑音をv_i∈C^Mにより表す。ただし、Cは複素数の全体の集合である。この記法を用いて、拡散性雑音抑圧の定義を言い換えると、「Y:={y_i}^I _i=1が与えられたときに、X:={x_i}^I _i=1を推定する問題」ということになる。ここで、記法{y_i}^I _i=1は、{y_i|i=1,…,I}を表す。ここで、x_i、v_i、y_iは、いずれも、M個のマイクロホン位置における各種信号を並べたM次元のベクトルである。 本実施形態では、観測信号を収録する音場は、静止した単一の点音源に起因する音源信号と、拡散性雑音とからなり、観測信号y_iは、次式で表されるとモデル化する。
y_i=x_i+v_i (b1)
つまり、観測信号y_iを、拡散性雑音v_iと拡散性雑音を含まない信号（ノイズフリー信号）である音源信号x_iとの和でモデル化する。ここで、簡単のため、音源信号x_iと拡散性雑音v_iとについて、以下の仮定を置く。
・時間的独立性:{x_i}^I _i=1は独立な系列である。すなわち、任意の相異なるi、jに対し（i∈{1,…,I}、j∈{1,…,I}、i≠j）、音源信号x_iと音源信号x_jとは独立である。また、{v_i}^I _i=1も独立な系列である。
・相互独立性:{x_i}^I _i=1と{v_i}^I _i=1とは、互いに独立である。すなわち、任意のi、jに対し（i∈{1,…,I}、j∈{1,…,I}）、音源信号x_iと拡散性雑音v_jとは独立である。
・ガウス性:音源信号x_iと拡散性雑音v_iとは、平均が0の複素ガウス分布に従う。
ここで、平均μ∈C^M、共分散行列Σ∈C^M×Mの複素ガウス分布の確率密度関数は、次式で与えられる。

ただし、det・は行列・の行列式、・^Hは行列・のエルミート転置、exp[・]は自然対数の底数の・乗である。

さらに、短時間定常性をモデル化するため、I個のフレームを、連続するフレームからなる、T個のブロックに区分し、各ブロックにおいて、音源信号x_i及び拡散性雑音v_iの共分散行列は一定と仮定する。ブロックの番号をt=1,…,T、各ブロックt内のフレームの番号をu=1,…,Uとするとき、iと(t,u)との対応

すなわち、一般に、iをUで割ったときの商および剰余をそれぞれqおよびrとするとき、

である。このように定められたt、uを用いて、x_i→x_tu、v_i→v_tuのように番号を振りなおすとき、短時間定常性は次のように表すことができる。
・短時間定常性:音源信号x_tuおよび拡散性雑音v_tuの共分散行列は、ブロックの番号tには依存するが、各ブロック内でのフレームの番号uには依存せず、Φ^x _t、Φ^v _t∈C^M×Mと置ける。

なお、フレームの総数Iは、Uの整数倍とは限らないため、端数は最後のブロックで調整するものとする。すなわち、最後のブロックは、I-U(T-1)個のフレームからなる。

上述の仮定のもと、観測信号y_tuが与えられたときの、音源信号x_tuの最大事後確率推定量を導くことができる。これは、観測信号y_tuが与えられたときの音源信号x_tuの事後確率p(x_tu|y_tu;Θ)を最大化するx_tuと定義される。そこで、この事後確率を導くため、まず同時確率を導くと、

となる。式(b5)およびベイズの定理より、事後確率は、

ここで、

であり、

は、x_tuに依存しない定数を除いて等しいことを表す。よって、

であるから、x_tuの最大事後確率推定量は、(b8)で与えられる。なお、最大事後確率推定量(b10)は、最小平均二乗誤差推定量であるマルチチャネルウィーナーフィルタ(参考文献１〜４参照)と一致している。
(参考文献１) N.Q.K. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank spatial covariance model”, IEEE Trans. ASLP, Sep. 2010, vol. 18, no. 7, pp. 1830-1840.
(参考文献２)H.L. Van Trees, "Optimum Array Processing", John Wiley & Sons, NewYork, 2002.
(参考文献３)K.U. Simmer, J. Bitzer, and C. Marro, “Post-filtering techniques”, in
Microphone Arrays, M. Brandstein and D.Ward, Eds., pp. 39-60. Springer, Berlin Heidelberg, 2001.
(参考文献４)S. Doclo and M. Moonen, “GSVD-based optimal filtering for single and multimicrophone speech enhancement”, IEEE Trans. SP, Sep. 2002, vol. 50, no. 9,pp. 2230-2244.
これは、ガウス分布の平均と最頻値が一致することに起因する。

マルチチャネルウィーナーフィルタ(b8)の設計は、共分散行列Φ^x _tおよびΦ^v _tの推定に帰着する。この推定問題は、観測信号y_tuを用いて推定できる観測信号y_tuの共分散行列を

と置くと、次のように定式化できる。

ここで、D(・,・)は何らかの距離尺度である。また、Σ_uは、t=1,…,T-1に対しては、Σ^U _u=1を意味し、t=Tに対しては、Σ^I-U(T-1) _u=1を意味する。(b12)において、共分散行列Φ^x _tとΦ^v _tとは共に未知である。したがって、(b12)は劣決定の逆問題であり、制約条件がない場合、(b12)の最適解は無数に存在する。実際、^Φ^x _t、^Φ^v _tが最適解ならば、Φ^x _t+Φ^v _t=^Φ^x _t+^Φ^v _tを満たす任意の共分散行列Φ^x _t、Φ^v _tは最適解である。したがって、この逆問題を解くためには、共分散行列Φ^x _tとΦ^v _tとの性質を利用して、これらを適切にモデル化する必要がある。

＜共分散行列を用いた信号と雑音のモデル化＞
共分散行列Φ^x _tとΦ^v _tとをモデル化するために、本実施形態では、信号と雑音の空間的な性質を利用する。

クリーン信号の信号源が、静止した単一の点音源であるという仮定の下では、音源信号の空間的な特性は、時間によらず一定と考えられる。この空間的な特性の時間的な不変性に着目して、次式のような共分散行列Φ^x _tのモデルを考えることができる。

ここで、φ^x _t∈R(Rは実数の全体の集合)は、スペクトルに対応する、時間依存のパラメータであり、B∈C^M×Mは、空間的な特性を表す、時間非依存のパラメータである。以下、φ^x _t及びBをそれぞれ信号パワー及び信号コヒーレンス行列ともいう。

非特許文献１でも(b13)と類似の共分散行列のモデルが用いられているが、本実施形態における音源信号の共分散行列のモデルでは、(b13)において、信号コヒーレンス行列Bに何の制約も課さないのに対し、非特許文献１に記載の音源信号の共分散行列のモデルでは、(b13)において、信号コヒーレンス行列Bがランク１行列hh^Hであるという制約を課している点が異なる（ｈはステアリングベクトル）。ランク１行列の逆行列は定義できないため、非特許文献１の音源信号の共分散行列のモデルに対しては、音源信号をガウス分布でモデル化することができず、したがって、音響信号に適した尺度である板倉-齋藤距離を用いることができない。これに対し、本実施形態では、信号コヒーレンス行列Bに上記のような制約を課さないため、板倉-齋藤距離を適用することができる。

音源信号とは異なり、拡散性雑音は、不特定多数の音源により生じるため、たとえそれらの音源がすべて静止していたとしても、各フレームで異なる音源がアクティブとなることにより、雑音全体の空間的な特性は時間的に変化する。また、多数の話者が同時に発話する場合のように、拡散性雑音は、スペクトルも時間変化することが多い。このように、拡散性雑音は、一般にスペクトルも空間的な特性も時間的に変化するため（非定常であるため）、時間的な不変性に基づかずにモデル化する必要がある。これに対し、非特許文献１及び参考文献５では、等方性という空間的な不変性に着目して、拡散性雑音をモデル化することが提案された。
（参考文献５）N. Ito, "Robust Microphone Array Signal Processing against Diffuse Noise", Ph.D. thesis, the University of Tokyo, 2012.
拡散性雑音は、不特定多数の音源に起因するため、方向依存性は小さく、等方的とみなせる。そこで、この等方性のモデルとして、非特許文献１及び参考文献５では、「拡散性雑音の二点間のクロススペクトルは、その二点間の方向には依存せず、その距離のみで決まる」と仮定する。この仮定のもとで、「Φ^v _tは、M×Mエルミート行列全体がなす実ベクトル空間Hの、ある低次元部分空間Vに属する」ことが示せる（非特許文献１及び参考文献５参照）。本実施形態においても非特許文献１及び参考文献５と同じ拡散性雑音のモデルを仮定する。すなわち、我々の拡散性雑音モデルは、次式で表される。
Φ^v _t∈V (b14)
Vは、行列がなすベクトル空間Hの部分空間であるから、行列部分空間と呼ぶ。換言すると、Vは、行列のなす線型空間における部分空間である。

行列部分空間Vの具体形は、アレイ配置や雑音場に関する付加的な仮定などによって、複数提案されている。まず、任意のアレイ配置に適用可能で、雑音に関する付加的な仮定も置かない、最も一般的な場合が、実数値雑音共分散（real-valued noise covariance）モデル

である。また、結晶型（参考文献５及び６参照）と呼ばれる、特別なクラスのアレイ配置の場合には、ブラインド雑音無相関化（blind noise decorrelation:BND）モデル

が適用できる。ここで、Uは、既知の定数ユニタリ行列である。
（参考文献６）N. Ito, H. Shimizu, N. Ono, and S. Sagayama, “Diffuse noise suppression using crystal-shaped microphone arrays”, 2011, vol. 19, no. 7, pp. 2101-2110.
また、波長に比べ十分距離の大きい2点間では、拡散性雑音は無相関である、と仮定すると、マイクロホン間の距離を大きくした極限で、空間的無相関雑音（spatially uncorrelated noise）モデル

が成り立つ。最後に、マイクロホン間の相対的な位置関係が既知で、3次元実空間R³の任意の方向から、スペクトルの等しい雑音が、平面波として到来するという仮定の下では、定数雑音コヒーレンス（fixed noise coherence）モデル

が成立する。ここで、Γはコヒーレンス行列と呼ばれ、その(m,n)成分は、

で与えられる。ここで、sinc(・)はsinc関数である。また、L_mnはm番目とn番目のマイクロホンの距離、cは音速であり、fは観測信号y_tuに対応する周波数であり、L_mnとしては実測値を用い、cとしては、例えば、近似値c=340m/sを用いればよい。

これらの雑音モデルの比較については、参考文献５を参照されたい。本実施形態の方法は、これらの全ての雑音モデルに適用可能である。

これらの拡散性雑音の共分散行列のモデルは、拡散性雑音の共分散行列が時間的に変化することを許しているため、非定常の拡散性雑音も扱うことができるという大きな利点がある。従来の多くの雑音抑圧技術では、定常雑音か、少数の点音源から生じる方向性雑音を仮定することが多かった。しかしながら、駅・空港、乗り物、オフィス、工場、カフェ、居酒屋、繁華街といった、さまざまな実環境において、非定常の拡散性雑音が存在する。定常雑音か、少数の点音源から生じる方向性雑音を仮定する従来技術では、このような雑音は扱えないため、このことが、雑音抑圧技術の実環境での応用範囲を大きく制限していた。これに対し、本実施形態における雑音の共分散行列のモデル化では、雑音の空間的な性質に着目することで、このような非定常の拡散性雑音もモデル化できる。

本実施形態における雑音の共分散行列のモデル化に基づくモデル推定の方法として、すでに非特許文献１において、ユークリッド距離を評価尺度とする手法が提案されていた。しかし、既述の通り、非特許文献１では音源信号のモデルに制約を設けているため板倉‐齋藤距離を最小化することはできない。つまり、非定常の拡散性雑音を扱えるようなモデルにおいて、板倉-齋藤距離を最適化する手法の導出が可能であるかどうかは、知られていなかった。

＜拡散性雑音抑圧の従来技術＞
従来技術の課題についてより詳細に説明する。

拡散性雑音抑圧は、音源信号の共分散行列Φ^x _tと拡散性雑音の共分散行列Φ^v _tとを推定する問題に帰着させることができる。すなわち、これらを推定することができれば、マルチチャネルウィーナーフィルタ(b8)により、観測信号y_tuを用いて、音源信号x_tuを推定することができる。この推定値の精度は、共分散行列Φ^x _tとΦ^v _tとの推定精度に依存するため、共分散行列Φ^x _tとΦ^v _tとを正確に推定することが非常に重要である。

非特許文献１では、この問題を解くために、次式で定義されるコスト関数を用いていた。

簡単のため、フレーム数Iが、ブロックあたりのフレーム数Uで割り切れ、端数が生じない場合で説明した。さらに、残響はフレーム長に比べて十分短いなどの仮定の下、前述の雑音モデルの下で、(b20)の最小化により共分散行列Φ^x _tとΦ^v _tとを推定するための反復アルゴリズムを導くことができる。

しかし、前述の通り、ユークリッド距離の平方ではなく、対数パワー領域での誤差を小さくする必要がある。

＜本実施形態のポイント＞
本実施形態では、非特許文献１のように、ユークリッド距離の平方の総和(b20)を最小化するのではなく、最尤法により音源信号の共分散行列Φ^x _tと拡散性雑音の共分散行列Φ^v _tとを推定する。後述のように、この最尤推定は、次式のコスト関数の最小化と等価である。

このコスト関数は、観測信号y_tから得られる共分散行列Φ^y _tと、Φ^x _t＋Φ^v _tとの板倉-齋藤距離の、すべてのブロックtについての総和を表し、本実施形態では、この値が小さくなるように共分散行列Φ^x _tと共分散行列Φ^v _tとを更新する。ここで、D_ISは、次式のマルチチャネル板倉-齋藤ダイバージェンスである（参考文献７〜９参照）。

（参考文献７）H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Multichannel extensions of non-negative matrix factorization with complex-valued data”, IEEE Trans. ASLP, May 2013, vol. 21, no. 5, pp. 971-982.
（参考文献８）K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto, “Infinite positive semidefinite tensor factorization for source separation of mixture signals”, in Proc. International Conference on Machine Learning (ICML), Jun. 2013, pp. 576-584.
（参考文献９）B. Kulis, M. Sustik, and I. Dhillon, “Low-rank kernel learning with Bregman matrix divergences”, Journal of Machine Learning Research, Feb. 2009, vol. 10, pp. 341-376.

つまり、本実施形態では、観測信号y_ｔを拡散性雑音v_ｔと拡散性雑音を含まない信号であるノイズフリー信号x_ｔとの和でモデル化し、観測信号y_ｔから得られる共分散行列Φ^ｙ _tと、拡散性雑音の共分散行列Φ^v _tとノイズフリー信号の共分散行列Φ^x _tとの和（Φ^x _t＋Φ^v _t）と、の板倉-齋藤距離が小さくなるように各モデルパラメータ（ここではΦ^v _tとΦ^x _t）を更新する。本実施形態では、拡散性雑音の共分散行列を特定するモデルパラメータはΦ^v _tそのものであり、ノイズフリー信号の共分散行列を特定するモデルパラメータはΦ^x _tそのものである。拡散性雑音の共分散行列Φ^v _tは、観測信号y_ｔを拡散性雑音v_ｔと音源信号x_ｔとの和でモデル化したときの拡散性雑音v_ｔの統計的性質を表すモデルパラメータであり、音源信号の共分散行列Φ^x _tは上記モデル化における音源信号x_ｔの統計的性質を表すモデルパラメータであるともいえる。

また、簡単のため、フレーム数Iが、ブロックあたりのフレーム数Uで割り切れ、端数が生じない場合で説明した。スカラーk>0に対して、D_IS(kA,kB)=D_IS(A,B)が成り立つことから、D_ISは、上記の音響信号の対数性を反映した、適切な尺度となっていることが分かる。なお、(b22)は、NTT(登録商標)の板倉・齋藤(参考文献１０参照)により提案され、音声処理における有効性が示されてきた、板倉-齋藤ダイバージェンスを、複数チャネルの場合へ拡張した尺度とみなせる。
(参考文献１０)F. Itakura and S. Saito, “Analysis synthesis telephony based on the maximum likelihood method”, in Rep. 6th International Congress on Acoustics, 1968, pp. C-17-C-20.

＜本実施形態における共分散行列推定の原理＞
前述のように、本実施形態では、最尤法によりΘ:={{Φ^x _t}^T _t=1,{Φ^v _t}^T _t=1}={{φ^x _t}^T _t=1,B,{Φ^v _t}^T _t=1}を推定する（(b13）参照)。目的関数である対数尤度は、次式により与えられる。

以下では、簡単のため、フレーム数Iが、ブロックあたりのフレーム数Uで割り切れ、端数が生じない場合で説明するが、端数が生じる場合への拡張は容易である。この場合、上式は、以下のように単純化できる。

したがって、J(Θ)の最大化は、(b21)の最小化と等価である。拘束条件(b14)に注意して、Θを推定する問題は、以下のように定式化される。

J(Θ)は、log det(Φ^x _t+Φ^v _t)および(Φ^x _t+Φ^v _t)^-1なる項を含むため、各パラメータに関するJ(Θ)の偏微分は複雑であり、(b27)の最適解を閉形式で求めることは困難である。一方、音源信号x_tuを隠れ変数とみなし、expectation-maximization(EM)アルゴリズム(参考文献１１参照)を用いれば、効率的な最適化が可能である。
（参考文献１１）A.P. Dempster, N.M. Laird, and D.B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm”, Journal of the Royal Statistical Society: Series B (Methodological), 1977, vol. 39, no. 1, pp. 1-38.
EMアルゴリズムは、次のEステップとMステップを交互に反復するものであり、(b27)の局所解への収束が保証される。
・Eステップ:現在のパラメータの推定値Θ'と(b10)とを用いて、音源信号x_tuの事後確率p(x_tu|y_tu;Θ')を計算する。本実施形態における(b10)はガウス分布であり、平均μ^x|y _tuおよび共分散行列Φ^x|y _tにより完全に決定されるため、Θ'と(b8)(b9)とを用いて、推定値(μ^x|y _tu)'および(Φ^x|y _t)'を計算すれば十分である。ただし、(・)'は、現在のパラメータの推定値Θ'を用いて計算することを示す。
・Mステップ:Eステップで得た事後確率p(x_tu|y_tu;Θ')に関する、同時分布の対数logp(x_tu,y_tu;Θ)の期待値として定義されるQ関数

が最大となるようΘを更新する。但し、<・>_{p(x_tu|y_tu;Θ')}（下付添え字のp(x_tu|y_tu;Θ')は、p(x_tu|y_tu;Θ')を表す）は、事後確率p(x_tu|y_tu;Θ')を用いて計算した、音源信号x_tuに関する期待値を表す。換言すると、モデルパラメータの集合Θが与えられたときの観測信号y_tuの集合Yと音源信号x_tuの集合Xとの同時分布の対数log(p(X,Y;Θ))の集合Xに関する期待値が最大となるように、共分散行列Φ^x _tと共分散行列Φ^v _tとを更新する。

このEMアルゴリズムにより、パラメータの推定値Θ'は、(b27)の局所最適解に収束することが保証されている。

Mステップにおけるパラメータの更新式を導くため、Q関数の具体形を求める。(b5)を(b28)に代入すれば、

ここで、定数項を無視し、なおかつ、<・>_{p(x_tu|y_tu;Θ')}を<・>と略記した。なお、Tr[・]は、行列・の行列の対角成分の和を表す。(b10)より、(b29)における<x_tux_tu ^H>および<(y_tu-x_tu)(y_tu-x_tu)^H>は、Eステップで計算した(μ^x|y _tu)'および(Φ^x|y _t)'を用いて、下記のように計算できる。

Q関数(b29)を、各パラメータに関して最大化することにより、更新式を導ける。

まず、Q関数(b29)の信号パワーφ^x _tに関する偏微分を0と置くことにより((b13)参照)、

これを信号パワーφ^x _tについて解くと、

ただし、

と置いた。次に、信号コヒーレンス行列Bに関する偏微分より、次式を得る((b13)参照)。

これをBについて解くと、次式を得る。

拡散性雑音の共分散行列Φ^v _tの更新式は、行列部分空間Vの具体形に依存する。以下では、まず、

なる性質を持つクラスの行列部分空間Vに対しては、統一的な更新式を導出できることを示す。

このクラスは、実数値雑音共分散モデルV_real、ブラインド雑音無相関化モデルV_BND、空間的無相関雑音モデルV_uncorを含む。次に、上記のクラスに含まれない定数雑音コヒーレンスモデルV_cohに対する更新式を、別途導出する。

行列部分空間Vが(b40)を満たすとき、(Φ^v _t)^-1∈Vに注意すると、Q関数(b29)の共分散行列Φ^v _tに依存する項は、次のように変形できる。

ただし、P[・]は、Hの標準的な内積<A,B>:=Tr[AB]に関する、行列部分空間Vへの正射影を表す。(b41)におけるP[・]内の項は、一般に、行列部分空間Vに属する成分と、これに直交する成分との両方を持つが、(Φ^v _t)^-1∈Vとの内積により後者は消えるため、(b41)のようになる。(b41)を最大化する共分散行列Φ^v _t∈Vを求めるために、拘束条件Φ^v _t∈Vを一旦忘れて、(b41)をΦ^v _tで偏微分すると、次式を得る。

ただし、

である。P[・]の定義より明らかなように、(b42)は、確かにΦ^v _t∈Vを満たす。

上記では、偏微分を用いて更新式(b42)を導いたが、以下のように、より直観的に導くこともできる。(b41)を-Uで割れば、共分散行列Φ^v _tに依存しない定数項を除いて、下記のマルチチャネル板倉-齋藤ダイバージェンス(参考文献７〜９参照)に一致する。

したがって、(b41)の最大化は、(b45)の最小化と等価である。D_IS(・,・)は非負であり、二つの引数が等しいときに限り0となる。P[^Φ^v _t]は共分散行列Φ^v _tの実行可能領域である行列部分空間Vに属するから、(b45)はΦ^v _t=P[^Φ^v _t]のときに最小となる。

次に、行列部分空間Vとして定数雑音コヒーレンスモデルを用いる場合、すなわち、

と表される場合の共分散行列Φ^v _tの更新式を導出する。ここで、φ^v _tは未知の雑音パワーであり、Γは既知と仮定する雑音コヒーレンス行列である。Q関数(b29)に(b46)を代入して、未知の雑音パワーφ^v _tに関して偏微分すると、(b35)の導出と同様に、次式を得る。

以上で導出したEMアルゴリズムの一回の反復を以下にまとめる。
[Eステップ]

[Mステップ]

なお、Mステップにおける、行列部分空間Vへの直交射影は、次式により計算できる。

ここで、Q_dは、行列部分空間Vの正規直交基底であり、DはV行列部分空間の次元である。なお、Q_dの具体形については、非特許文献１、参考文献５を参照されたい。

＜本実施形態に係るモデル推定装置１００＞
図１は第一実施形態に係るモデル推定装置１００の機能ブロック図、図２はその処理フローの例を示す図である。モデル推定装置１００は、周波数領域変換部１１０、事後確率更新部１２０、パラメータ更新部１３０、パラメータ保持部１４０を含む。また、事後確率更新部１２０、パラメータ更新部１３０、パラメータ保持部１４０をまとめて、モデル推定部１５０と称する。

＜周波数領域変換部１１０＞
周波数領域変換部１１０は、M個のマイクロホンで取得した、時間領域の観測信号~y_τ ^(m)(m=1,…,M)を受け取り、これらの値を用いて、短時間フーリエ変換などの時間周波数変換により、周波数領域の観測信号y_tを計算し（ｓ１１０）、事後確率更新部１２０及びパラメータ更新部１３０に出力する。ここで、M>1であり、yの上の〜は時間領域における表現であることを表し、τは時間領域におけるサンプルの番号である。また、前述の通り、フレームは、ブロックの番号tとブロック内でのフレームの番号uとによって指定し、周波数ビンの番号は省略している。

＜事後確率更新部１２０＞
事後確率更新部１２０は、周波数領域の観測信号y_tuを受け取り、パラメータ保持部１４０から現在のパラメータ(要は、EMアルゴリズムの反復処理によって更新する中で、現在のパラメータであることを意味し、「直近に求めたパラメータ」「最新のパラメータ」と言い換えてもよい)の集合の推定値Θ'を取り出す。これらの値を用いて、観測信号y_tuが与えられたときの音源信号x_tuの事後確率p(x_tu|y_tu;Θ')を更新し、パラメータ更新部１３０に出力する。上述の通り、ガウス分布の場合、実際には、事後確率の平均μ^x|y _tuおよび共分散行列Φ^x|y _tを更新して出力すれば十分である。

以下、事後確率更新部１２０における処理を、詳細に説明する。図１のように、事後確率更新部１２０は、事後平均更新手段１２１と、事後共分散行列更新手段１２２とを含む。なお、事後確率更新部１２０における初めての処理に先立ち、パラメータの集合Θの初期値を設定し（ｓ１）、パラメータ保持部１４０に用意しておく。この初期値は、例えば、非特許文献１により推定した、信号のステアリングベクトル^h、音源信号のパワースペクトログラム^φ^x _t、雑音共分散行列^Φ^v _tを用いて、以下のように計算すればよい。

ここで、(b50)の_εI_M×Mの加算は、信号コヒーレンス行列Bの初期値が特異行列となってアルゴリズムが不安定になるのを防ぐ役割がある。_εは正数で、例えば10^-4とすればよく、I_M×MはM×Mの単位行列を表す。

事後平均更新手段１２１は、周波数領域の観測信号y_tuと、現在のパラメータの集合の推定値Θ'に含まれる音源信号の共分散行列Φ^x _tと拡散性雑音の共分散行列Φ^v _tとを用いて、事後確率の平均μ^x|y _tuを次式により更新して出力する（ｓ１２１）。

事後共分散行列更新手段１２２は、現在のパラメータの集合の推定値Θ'に含まれる共分散行列Φ^x _tとΦ^v _tとを用いて、事後確率の共分散行列Φ^x|y _tを、次式により更新して出力する（ｓ１２２）。

＜パラメータ更新部１３０＞
パラメータ更新部１３０は、観測信号y_tuと、事後確率更新部１２０で計算した事後確率の平均μ^x|y _tuと共分散行列Φ^x|y _tとを受け取る。また、パラメータ保持部１４０から現在のパラメータの集合の推定値Θ'を取り出す。パラメータ更新部１３０は、これらの値を用いて、パラメータの集合Θを更新して、パラメータ保持部１４０に出力する。また、所定の更新回数を終了している場合には（ｓ２）、更新したパラメータの集合Θを最終的なパラメータの推定値として出力する。なお、更新の対象をパラメータの集合Θと呼び、その際に用いる、過去に求めた（更新した）パラメータの集合を推定値Θ'と呼ぶ。以下、パラメータ更新部１３０における処理を、詳細に説明する。図１に示すように、パラメータ更新部１３０は、信号パワー更新手段１３１、信号コヒーレンス行列更新手段１３２、信号共分散行列更新手段１３３、雑音共分散行列更新手段１３４を含む。

信号パワー更新手段１３１では、事後確率の平均μ^x|y _tu、共分散行列Φ^x|y _t、信号コヒーレンス行列Bを用いて、信号パワーφ^x _tを更新して出力する（ｓ１３１）。まず、

により、^Φ^x _tを計算する。次に、

により、信号パワーφ^x _tを更新する。

信号コヒーレンス行列更新手段１３２は、信号パワー更新手段１３１で計算した^Φ^x _tと信号パワーφ^x _tとを用いて、信号コヒーレンス行列Bを次式により更新して出力する（ｓ１３２）。

信号共分散行列更新手段１３３は、信号パワー更新手段１３１で計算したφ^x _tと、信号コヒーレンス行列更新手段１３２で計算した信号コヒーレンス行列Bとを用いて、音源信号の共分散行列Φ^x _tを次式により更新して出力する（ｓ１３３）。

雑音共分散行列更新手段１３４は、事後確率の平均μ^x|y _tuと共分散行列Φ^x|y _t、および観測信号y_tuを用いて、拡散性雑音の共分散行列Φ^v _tを更新して出力する（ｓ１３４）。まず、

により、^Φ^v _tを計算する。次に、

により、拡散性雑音の共分散行列Φ^v _tを更新して出力する。ここで、P[^Φ^v _t]は、Q_dをVの正規直交基底として、

により計算され、コヒーレンス行列Γの(m,n)成分γ_mnは、

により計算される。L_mnはm番目とn番目のマイクロホンの距離、cは音速、ｆは観測信号y_tuの周波数である。＜拡散性雑音抑圧の問題設定＞で説明したとおり、本実施形態のモデル推定装置には、周波数ビンごとの観測信号が入力され、周波数ビンごとに処理が行われる。つまり、観測信号y_tuはある周波数ビンの観測信号であり、ｆはその周波数ビンの周波数を表す。

パラメータ更新部１３０は、上述の信号パワー更新手段１３１、信号コヒーレンス行列更新手段１３２、信号共分散行列更新手段１３３、雑音共分散行列更新手段１３４の各部の処理により得られたパラメータの集合Θ:={{Φ^x _t}^T _t=1,{Φ^v _t}^T _t=1}を出力する。

＜パラメータ保持部１４０＞
パラメータ保持部１４０は、パラメータ更新部１３０での更新処理により得られたパラメータの集合Θを受け取り、保持し（ｓ１４０）、事後確率更新部１２０とパラメータ更新部１３０とにおける次回の処理の際に提供する。

＜効果＞
このような構成により、従来技術よりも推定精度よく雑音抑圧時に用いるモデルパラメータの推定することができる。

＜第二実施形態＞
第一実施形態と異なる部分を中心に説明する。

本実施形態では、第一実施形態で推定したパラメータの集合Θ:={{Φ^x _t}^T _t=1,{Φ^v _t}^T _t=1}をを用いて、観測信号y_tuに含まれる雑音を抑圧する雑音抑圧装置２００について説明する。

図３は雑音抑圧装置２００の機能ブロック図、図４はその処理フローの例を示す図である。雑音抑圧装置２００は、周波数領域変換部１１０、モデル推定部１５０、平均計算部２６０、時間領域変換部２７０を含む。

周波数領域変換部１１０、モデル推定部１５０の処理内容は第一実施形態で説明した通りである。

モデル推定部１５０は、反復終了後のパラメータの集合Θ:={{Φ^x _t}^T _t=1,{Φ^v _t}^T _t=1}を出力し、平均計算部２６０での処理に供する。

なお、図４において、事後確率更新部１２０、パラメータ更新部１３０及びパラメータ保持部１４０における処理をまとめて、モデル推定処理（ｓ１５０）という。

＜平均計算部２６０＞
平均計算部２６０は、周波数領域変換部１１０から出力された観測信号y_tuとモデル推定部から出力されたパラメータの集合Θ:={{Φ^x _t}^T _t=1,{Φ^v _t}^T _t=1}を用いて、式(b52)により、
観測信号y_tuが与えられたときの音源信号x_tuの事後確率の平均μ^x|y _tuを得て出力する。

＜時間領域変換部２７０＞
時間領域変換部２７０は、平均計算部２６０の出力である事後確率の平均μ^x|y _tuを受け取り、この値に対し、逆短時間フーリエ変換（inverse STFT）などの時間周波数変換の逆変換を適用し、時間領域における音源信号の推定値である、時間領域の信号^x_τ∈R^Mに変換し（ｓ２７０）、これを雑音抑圧装置２００の出力値として出力する。

＜効果＞
以上のように、本実施形態では、非特許文献１におけるようなユークリッド距離の平方ではなく、より音響信号に適したマルチチャネル板倉-齋藤ダイバージェンスを用いることにより、対数的な音響信号に、より適した処理を行うことが可能である。さらに、本実施形態は、観測信号の確率的生成モデルに基づくため、確率的生成モデルに基づく他の音声強調手法（例:音源分離(参考文献１及び７参照)、残響除去（参考文献１２参照）との統合に向いている（つまり、第一実施形態で求めたモデルパラメータは拡散性雑音状況下の音源分離、残響除去等の技術に用いることができる）。
（参考文献１２）T. Nakatani, T. Yoshioka, K. Kinoshita, M. Miyoshi, and B.-H. Juang, “Speech dereverberation based on variance-normalized delayed linear prediction”, IEEE Trans. ASLP, Sep. 2010, vol. 18, no. 7, pp. 1717-1731.
このような統合により、様々な音響的事象を含む実環境において、広く適用可能な新しい音声強調の枠組みを構築することができると期待されるため、大きな利点である。

＜実験結果＞
本実施形態の効果を確かめるため実験を行った。観測信号としては、REVERB challenge(参考文献１３参照)のデータベース中の、AMI_WSJ20-Array1-*_T10c030x.wav（*=1,…,8）を用いた。
(参考文献１３)K. Kinoshita, M. Delcroix, T. Yoshioka, T. Nakatani, E.A.P. Habets,R. Hab-Umbach, V. Leutnant, A. Sehr, W. Kellermann, R. Maas, S. Gannot, and B. Raj, "The REVERB challenge: A common evaluation framework for dereverberation and recognition of reverberant speech", in Proc. WASPAA, Oct. 2013.

これは、雑音・残響を含む実環境において、一人の話者の音声を、8チャネルのマイクロホンアレイを用いて収録した信号である。データ長は2秒で、サンプリング周波数は16kHzであった。

本実施形態と非特許文献１のそれぞれにより推定した共分散行列を用いて設計した、マルチチャネルウィーナーフィルタによる雑音抑圧性能を比較した。短時間フーリエ変換および逆短時間フーリエ変換のフレーム長、フレームシフト、窓は、それぞれ、1024点（64ms）、32点（2ms）、Hamming窓とした。また、ブロック長は32フレームとした。アルゴリズムの反復回数は20回とした。雑音抑圧性能の評価指標としては、ケプストラム距離（cepstrum distance:CD）を用いた（定義は参考文献１３参照）。その際、参照信号(雑音・残響をほぼ含まない信号)としては、ヘッドセットによる収録音AMI_WSJ20-Headset1_T10c030x.wavを用いた。ケプストラム距離が小さい程、参照信号に近く、雑音抑圧性能が高いことを示す。図５に、観測信号、参照信号、処理音のスペクトログラムとケプストラム歪み（括弧内の数字）を示す。スペクトログラムの横軸は時間（秒）、縦軸は周波数（kHz）を表す。(a)は観測信号、(b)は非特許文献１による処理音、(c)は本実施形態による処理音、(d)は参照信号である。比較のため、手法間の差異が顕著であった0-4kHzの範囲に限って図示する。(b)の非特許文献１に比べて、(c)の本実施形態による処理音は、より拡散性雑音を効果的に抑圧できたことが分かる。これは、本実施形態による処理音のケプストラム距離3.9dBが、非特許文献１のケプストラム距離4.6dBよりも小さいことからもわかる。

＜第三実施形態＞
第一実施形態と異なる部分を中心に説明する。なお、必要に応じて、第一実施形態で定義した記号の一部を再定義する。

本実施形態では、拡散性雑音下での音声強調技術、特に、残響の影響を受けた単一のノイズフリー音と、拡散性雑音とからなるとモデル化される観測信号を用いて、雑音と残響の影響が除去されたクリーン音成分を推定する際に用いるモデルパラメータ（クリーン信号の共分散行列、拡散性雑音の共分散行列、予測係数行列）を推定する。音源の位置は、未知であってもよく、また、それは静止しているとモデル化する。なお、音源がゆるやかに移動する場合には、音源が静止していると近似しうる複数のブロックに観測信号を分割し、ブロックバッチ処理を適用することにより、性能の低下を抑制することができる。

＜拡散性雑音下での音声強調の問題設定＞
本実施形態では、特に断りのない限り、各信号を短時間フーリエ変換（short-time Fourier transform: STFT）などの時間周波数領域で表現する。フレームの総数をTで表し、フレームの番号をt∈{1,…,T}で表す。

以下のように変数を定義する。
・y_t∈C^M:M個のマイクロホンによる観測信号
・r_t∈C^M:残響の影響を含む、ノイズフリー信号
・x_t∈C^M:残響の影響を含まない、クリーン信号、すなわち、音源から発せられる直接音と初期反射成分とからなる信号
・v_t∈C^M:拡散性雑音
・G^k∈C^M×M:残響をモデル化する自己回帰過程の予測係数行列（k∈{1,…,K}はタップの番号、Kは予測次数）

ここで、y_t、r_t、x_t、v_tは、M個のマイクロホン位置における各種信号を並べたM次元のベクトルである。
このとき、観測信号のモデルは、数学的には次式により表現できる。

ここで、遅延Δは、直近のΔブロックを残響除去のための線形予測から除外することにより、残響に起因する時間相関のみを無相関化し、音声固有の相関を白色化するのを防ぐ効果がある（参考文献１２参照）。また、t≦0に対しては、r_t=0と定める。第一実施形態も第三実施形態も、観測信号を、非定常の拡散性雑音と拡散性雑音を含まない信号（ノイズフリー信号）との和でモデル化する点は共通するが、第一実施形態では残響を考慮せずに、拡散性雑音と音源信号との和でモデル化するのに対し、本実施形態では残響を考慮して、拡散性雑音とクリーン信号に残響が重畳された信号r_t（残響の影響を含む、(c2)参照）との和でモデル化する点が異なる。つまり、第一実施形態ではノイズフリー信号を音源信号と仮定しているのに対し、第三実施形態ではノイズフリー信号がクリーン信号に残響が重畳された信号であると仮定してモデル化している点が異なる。

次に、統計的信号処理の枠組みで音声強調を行うために、上記のモデルを、確率的生成モデルに拡張する。本実施形態では、クリーン信号x_t,ノイズフリー信号r_t,拡散性雑音v_t,観測信号y_tを確率変数とみなし、予測係数行列G_kを確定的な未知パラメータとみなす。{x_t}_1≦t≦Tが与えられると(c2)により{r_t}_1≦t≦Tが定まり、この{r_t}_1≦t≦Tと{v_t}_1≦t≦Tが与えられると(c1)により{y_t}_1≦t≦Tが定まる。したがって、(c1)(c2)を生成モデル化することは、{x_t}_1≦t≦Tと{v_t}_1≦t≦Tの確率分布を定義することに帰着する。この確率分布として、本実施形態では、簡単のため、以下の仮定を満たすクラスを考える。
・時間的独立性:{x_t}_1≦t≦Tは独立な系列である。すなわち、任意の相異なるt,t'に対し（t∈{1,…,T}、t'∈{1,…,T}、t≠t'）、x_tとx_t'とは独立である。また、{v_t}_1≦t≦Tも独立な系列である。
・相互独立性:2つの系列{x_t}_1≦t≦T,{v_t}_1≦t≦Tは、互いに独立である。すなわち、任意のt,t'に対し（t∈{1,…,T}、t'∈{1,…,T}）、x_tとv_t'とは独立である。
・ガウス性:x_tおよびv_tは、平均がともに0、共分散行列がΦ^x _tおよびΦ^v _tの複素ガウス分布に従う。

ここで、平均μ∈C^M、共分散行列Σ∈C^M×Mの複素ガウス分布の確率密度関数は、次式で与えられる。

以上の仮定により、観測信号の生成モデルを立てることは、結局、クリーン信号x_tおよび拡散性雑音v_tの共分散行列Φ^x _tおよびΦ^v _tのモデル化に帰着する。この詳細については、後述する。

以上のモデル化では、各フレームで異なる共分散行列Φ^x _t、Φ^v _tを用いていたが、信号の短時間定常性をモデル化するために、モデルを拡張する。各フレームが、フレームシフトをU倍にすることにより得られるU個のサブフレームよりなると仮定する。共分散行列Φ^x _t、Φ^v _tは、フレームの番号tには依存するが、各フレーム内のサブフレームの番号u∈{1,…,U}には依存しないと仮定する。また、各信号は、x_tu,r_tu,v_tu,y_tu∈C^Mのように、tとuを用いて表記し、同一フレームtに属するサブフレームu∈{1,…,U}に対するこれらの信号は、互いに独立で同一の分布に従うもの（independent and identically distributed: i.i.d.)とする。(c1)(c2)は、次式のように書き改められる。

上で定義した記法を用いて、本実施形態で扱う音声強調の問題を改めて述べると、観測信号Y:={y_tu}_{1≦t≦T,1≦u≦U}が与えられたときに、クリーン信号X:={x_tu}_{1≦t≦T,1≦u≦U}を推定する問題である。特に、そこでの主要な問題は、上記のモデルのパラメータ
Θ:={Θ_x,Θ_v,Θ_g} (c6)
Θ_x:={Φ^x _t}_1≦t≦T (c7)
Θ_v:={Φ^v _t}_1≦t≦T (c8)
Θ_g:={G_k}_1≦k≦K (c9)
の推定であり、音声強調技術の性能は、その推定精度に大きく左右される。

＜信号と拡散性雑音の共分散行列のモデル化＞
上述の通り、観測信号の生成モデルを立てることが、共分散行列Φ^x _tとΦ^v _tのモデル化に帰着する。非特許文献１では、信号と拡散性雑音の空間的な性質が利用されている。

クリーン音が、静止した単一の点音源から発せられるという仮定の下では、クリーン音の空間的な特性は、時間によらず一定である。ゆえに、クリーン信号のモデルは、第一実施形態と同じく(b13)を用いる。第一実施形態と同様に、本実施形態においても信号コヒーレンス行列Bは非特許文献１で用いているランク１行列であるという制約を課さない。

拡散性雑音については、第一実施形態と同じく(b14)によりモデル化する。

Vの具体形は、アレイ配置や雑音場に関する付加的な仮定などによって、複数提案されている。第一実施形態で説明したように、以下の四つのモデルが考えられる。
・実数値雑音共分散（real-valued noise covariance）モデル

・ブラインド雑音無相関化（blind noise decorrelation: BND）モデル

・空間的無相関雑音（spatially uncorrelated noise）モデル

・定数雑音コヒーレンス（fixed noise coherence）モデル

ここで、B^vは雑音コヒーレンス行列であり、その(m,n)成分は、

で与えられる。

これらの雑音モデルの比較については、参考文献５を参照されたい。本実施形態の方法は、これらの全ての雑音モデルに適用可能である。

＜本実施形態のポイント＞
上記の目的を実現するために、本実施形態では、最尤法によりモデルパラメータ（クリーン信号の共分散行列、拡散性雑音の共分散行列、予測係数行列）を推定する。この最尤法は、ユークリッド距離に基づく尺度と比べ、音声により適した尺度である、板倉-齋藤行列ダイバージェンス（参考文献７〜９）の最小化と等価である。すでに説明した通り、非特許文献１は残響が存在しない場合、すなわち本実施形態のモデルにおいてG_k=0とした場合において、ユークリッド距離に基づいてモデルパラメータを推定するものであった。比較のため、非特許文献１と同じ条件であるG_k=0の場合について述べると、本実施形態における最尤推定は、次式のコスト関数の最小化と等価である。

ここで、D_ISは、次式で定義される。

スカラーk>0に対して、D_IS(kA,kB)=D_IS(A,B)が成り立つことから、D_ISは、上記の音響信号の対数性を反映した、適切な尺度となっていることが分かる。

G_k≠0の場合にも、同様のことが言える。なお、(c19)は、NTTの板倉・齋藤により提案され（参考文献１０参照）、音声処理における有効性が示されてきた、板倉-齋藤ダイバージェンスを、複数チャネルの場合へ拡張した尺度とみなせるため、本実施形態では、板倉-齋藤行列ダイバージェンスと呼ぶ。

＜原理＞
＜目的関数＞
本実施形態では、パラメータの集合Θ:={{φ^x _t}_1≦t≦T,B^x,{φ^v _t}_1≦t≦T,{G_k}_1≦k≦K}を推定するために、最尤法に基づき、次式の対数尤度を最大化する。

すなわち、拘束条件(b14)に注意して、

によりΘを推定する。

前述の観測信号生成モデルに基づき、対数尤度(c20)の具体的な表式を導出する。{y_tu}_1≦t≦Tは、残響に起因する時間相関を持つため、独立な系列ではない。

そこで、{y_tu}_1≦t≦Tを並べて拡張したベクトル

の確率密度関数を導出することにする。まず、(c1)(c2)のモデルを、拡張されたベクトルを用いて書き直すと、次式のようになる。
~y_u=~r_u+~v_u (c23)
F^H~r_u=~x_u (c24)
ここで、~r_u,~v_u,~x_uは、~y_uと同様に定義する。つまり、
~r_u:=[r_Tu ^T…r_1u ^T]^T∈C^MT×1
~v_u:=[v_Tu ^T…v_1u ^T]^T∈C^MT×1
~x_u:=[x_Tu ^T…x_1u ^T]^T∈C^MT×1
である。また、F∈C^MT×MTは、ブロックテプリッツ行列であって、T²個のM×M行列をブロックに持ち、(i,j)番目のブロックは次式で与えられる。

ただし、Iは単位行列を表す。つまり、以下のように表される。

(c24)より、F^Hは雑音除去オペレータ、F^-H:=(F^H)^-1は残響付加オペレータと解釈できる。つまり、F^-Hはクリーン信号に重畳される残響の特徴を表すオペレータ（パラメータ）ともいえる。ここで、F^Hは対角要素がすべて1の上三角行列であり、したがってdet(F^H)=1であるから、F^Hは正則であることに注意する。次に、拡張されたベクトル~x_u,~v_uの確率密度関数は、前述の観測信号生成モデル、および~x_u,~v_uの定義より、次式で与えられる。
p(~x_u;Θ)=N_C(~x_u;0,~Φ^x) (c26)
p(~v_u;Θ)=N_C(~v_u;0,~Φ^v) (c27)
ここで、~Φ^x∈C^MT×MTは、ブロック対角行列であって、T²個のM×M行列をブロックに持ち、i∈{1,…,T}番目の対角ブロックはΦ^x _T-i+1に等しい。~Φ^vも同様に定義される。つまり、以下のように表される。

よって、観測モデル(c23)(c24)、および~x_u,~v_uの確率密度関数(c26)(c27)より、対数尤度は、

である。ここで、前述の観測信号生成モデルより、{y_tu}_1≦t≦Tは独立な系列ではないが、{~y_u}_1≦u≦Uは独立な系列であることに注意する。

なお、(c28)は、次式のように変形できる。

ここで、D_ISは前述の板倉-齋藤ダイバージェンスであり、

である。したがって、J(Θ)=logp(Y;Θ)の最大化は、D_IS(~Φ^y,F^-H~Φ^xF^-1+~Φ^v)の最小化と等価である。よって、本実施形態では、~Φ^yと、F^-H~Φ^xF^-1+~Φ^vとの板倉-齋藤距離が小さくなるように、共分散行列Φ^x _tと共分散行列Φ^v _tと予測係数行列G_kとを更新する。特に、G_k=0、すなわちF=Iの場合には、前述のように(c18)の最小化に帰着する。

つまり、本実施形態においても、第一実施形態と同様に、観測信号y_ｔを拡散性雑音v_tと拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化する。ただし、本実施形態におけるノイズフリー信号は、音源から発せられる直接音と初期反射成分とからなるクリーン信号x_tに残響が重畳された信号である。そして、本実施形態では、観測信号y_tから得られる共分散行列~Φ^ｙと、拡散性雑音の共分散行列~Φ^vとノイズフリー信号の共分散行列F^-H~Φ^xF^-1との和（F^-H~Φ^xF^-1+~Φ^v）と、の板倉-齋藤距離が小さくなるように各モデルパラメータ（ここではΦ^v _tとΦ^x _tとG_k）を更新する。上述の通り、~Φ^v はΦ^v _tにより特定される共分散行列なので、拡散性雑音の共分散行列を特定するモデルパラメータはΦ^v _tである。また、~Φ^x はΦ^ｘ _t により特定される共分散行列であり、F^-HとF^-1は予測係数行列G_kにより特定されるオペレータであるので、拡散性雑音を含まない信号(ノイズフリー信号)を特定するモデルパラメータはΦ^x _tとG_kである。なお、~Φ^vは観測信号を非定常の拡散性雑音とノイズフリー信号との和でモデル化したときの拡散性雑音の統計的性質を表しており、F^-H~Φ^xF^-1はノイズフリー信号の統計的性質を表しているといえる。

＜Q関数＞
R:={r_tu}_{1≦t≦T,1≦u≦U}を隠れ変数とみなし、expectation-maximization(EM)アルゴリズムにより（参考文献１１）、パラメータの集合Θの更新式を導くことができる。パラメータの集合Θの現在の推定値をΘ'で表すと、Q関数は、

で定義される。これを計算するために、まず、完全データの対数尤度関数を計算すると、

次に、隠れ変数の事後確率を計算する。上式と、ベイズの定理により、

となる。したがって、

である。ただし、

である。次式により、μ^r|y _tuおよびΦ^r|y _tt'を定義する。

(c35)および(c39)より、Q関数は

と書ける。ここで、

と定義した。

＜更新式の導出＞
式(c45)のQ関数を、各パラメータに関して最大化することにより、Mステップにおける更新式が得られる。信号パワーφ^x _tに関する偏微分より、

よって、

次に、信号コヒーレンス行列B^xに関する偏微分より、

よって、

拡散性雑音の共分散行列Φ^v _tの更新式については、雑音抑圧技術（第一実施形態参照）の場合と同様に、

を得る。なお、第一実施形態のように、「V=V_BND,V_real,V_uncor」を「Vが

なる性質を持つ場合」に置き換えてもよい。最後に、^Φ^x _tが~Gを含むことに注意して、~G^*に関する偏微分より（・^*は・の複素共役を表す）、次式を得る。
~G=D^-1N (c58)
ただし、DおよびNは、次式により定義される。

以上で導出したEMアルゴリズムの一回の反復を、以下にまとめる。
[Eステップ]

[Mステップ]

＜本実施形態に係るモデル推定装置３００＞
図６は第三実施形態に係るモデル推定装置３００の機能ブロック図、図７はその処理フローの例を示す図である。

モデル推定装置３００は、モデル推定装置の機能構成例を示す。モデル推定装置３００は、周波数領域変換部１１０、事後確率更新部３２０、パラメータ更新部３３０及びパラメータ保持部１４０を含む。また、事後確率更新部３２０、パラメータ更新部３３０、パラメータ保持部１４０をまとめて、モデル推定部３５０と称する。図７を参照して、モデル推定装置３００の動作例を、手続きの順に従って説明する。周波数領域変換部１１０、及びパラメータ保持部１４０の処理内容は第一実施形態と同様である。

＜事後確率更新部３２０＞
事後確率更新部３２０は、周波数領域変換部１１０で計算した、周波数領域の観測信号y_tuを受け取り、パラメータ保持部１４０に保持されている、現在のパラメータの集合の推定値Θ'を取り出し、これらの値を用いて、観測信号y_tuが与えられたときのノイズフリー信号r_tuの事後確率の平均μ^r|y _tuおよび共分散行列Φ^r|y _tt'を更新し、更新した値をパラメータ更新部３３０に出力する。

以下、事後確率更新部３２０における処理を、詳細に説明する。事後確率更新部３２０は、事後平均更新手段３２１と、事後共分散行列更新手段３２２とを含む。

なお、事後確率更新部３２０における初めての処理に先立ち、パラメータの集合Θの初期値を設定し（ｓ１）、パラメータ保持部１４０に保持しておく。この初期値は、例えば、非特許文献１により推定した、信号のステアリングベクトル^h、信号のパワースペクトログラム^φ^x _t、雑音共分散行列^Φ^v _tを用いて、以下のように計算すればよい。

ここで、_εは正数（例えば10^-4）であり、B^xの初期値が特異行列となってアルゴリズムが不安定になるのを防ぐ役割がある。

事後平均更新手段３２１は、周波数領域の観測信号y_tuと、現在のパラメータの集合の推定値Θ'とを用いて、事後確率の平均μ^r|y _tuを更新して出力する（ｓ３２１）。まず、定義に従って、Fおよび~Φ^xを更新する（(c25)参照、~Φ^x∈C^MT×MTは、ブロック対角行列であって、T²個のM×M行列をブロックに持ち、i∈{1,…,T}番目の対角ブロックはΦ^x _T-i+1に等しい）。次に、

により~Φ^rを更新する。次に、(c57)によりΦ^v _tを求める。上述の通り

なので、この処理は~Φ^vを更新することに相当する。つまり、(c57)によりΦ^v _tを求めることで~Φ^vを更新する。次に、~Φ^rと~Φ^vを用いて

により~μ^r|y _uを更新して出力する。ここで、

である。 (c43)の定義より、~μ^r|y _uはμ^r|y _tuを統合した行列であるので、~μ^r|y _uを更新することはμ^r|y _tuも同時に更新することを意味する。

事後共分散行列更新手段３２２では、~Φ^rと~Φ^vとを用いて、

により~Φ^r|yを更新する。(c44)の定義より、~Φ^r|yはΦ^r|y _tt'を統合した行列であるので、~Φ^r|yを更新することはΦ^r|y _tt'も同時に更新することを意味する（ｓ３２２）。

＜パラメータ更新部３３０＞
パラメータ更新部３３０は、事後確率更新部３２０で計算した事後確率の平均μ^r|y _tuと共分散行列Φ^r|y _tt'、ならびに、観測信号y_tuを受け取り、パラメータ保持部１４０から保持されている現在のパラメータの集合の推定値Θ'を取り出し、これらの値を用いて、パラメータの集合Θを更新し、パラメータの集合Θをパラメータ保持部１４０に格納する。また、所定の更新回数を終了している場合には（ｓ２）、更新したパラメータの集合Θを最終的なパラメータの推定値として出力する。以下、パラメータ更新部３３０における処理を、詳細に説明する。

図６のように、パラメータ更新部３３０は、信号パワー更新手段３３１、信号コヒーレンス行列更新手段３３２、信号共分散行列更新手段３３３、雑音共分散行列更新手段３３４及び予測係数行列更新手段３３５を含む。

信号パワー更新手段３３１では、事後確率更新部３２０から出力された事後確率の平均μ^r|y _tuと共分散行列Φ^r|y _tt'、および、パラメータ集合の現在の推定値Θ'を用いて、信号パワーφ^x _tを計算して出力する（ｓ３３１）。まず、

により、Ψ_tt'を計算する。次に、

により^Φ^x _tを更新して出力する。最後に、

により、信号パワーφ^x _tを更新して出力する。

信号コヒーレンス行列更新手段３３２では、信号パワー更新手段３３１で計算した^Φ^x _tと信号パワーφ^x _tとを用いて、信号コヒーレンス行列B^xを、次式により更新して出力する（ｓ３３２）。

信号共分散行列更新手段３３３では、信号パワー更新手段３３１で計算した信号パワーφ^x _tと、信号コヒーレンス行列更新手段３３２で計算した信号コヒーレンス行列B^xとを用いて、クリーン信号の共分散行列Φ^x _tを次式により更新して出力する（ｓ３３３）。

雑音共分散行列更新手段３３４では、事後確率更新部３２０から出力された事後確率の平均μ^r|y _tuと共分散行列Φ^r|y _tt'、パラメータの集合の現在の推定値Θ'、および観測信号y_tuを用いて、拡散性雑音の共分散行列Φ^v _tを更新して出力する（ｓ３３４）。まず、

により、^Φ^v _tを計算する。次に、

により、拡散性雑音の共分散行列Φ^v _tを更新する。ここで、P[^Φ^v _t]は、{Q_d}_1≦d≦DをVの正規直交基底として、

により計算され、信号コヒーレンス行列B^vの各成分b^v _mnは、

により計算される。

予測係数行列更新手段３３５では、信号パワー更新手段３３１で計算したΨ_tt'および信号パワーφ^x _tを用いて、予測係数行列G_kを更新して出力する（ｓ３３５）。まず、

によりDを更新する。次に、

によりNを更新する。そして、
~G←D^-1N (c80)
により~Gを更新する。(c52)の定義より、~GはG_kを統合した行列なので、~Gを更新することはG_kも同時に更新することを意味する。

＜効果＞
このような構成により、従来技術よりも推定精度よく音声強調時に用いるモデルパラメータの推定することができる。

＜第四実施形態＞
第三実施形態と異なる部分を中心に説明する。本実施形態では、第三実施形態に係るモデル推定装置によって推定されたモデルパラメータを用いて、雑音と残響の影響が除去されたクリーン信号を推定する音声強調装置４００について説明する。

図８は第四実施形態に係る音声強調装置４００の機能ブロック図、図９はその処理フローの例を示す図である。音声強調装置４００は、周波数領域変換部１１０、モデル推定部３５０、平均計算部４５０、強調音声生成部４６０、時間領域変換部４７０を含む。

なお、図９において、事後確率更新部３２０、パラメータ更新部３３０及びパラメータ保持部１４０における処理をまとめて、モデル推定処理（ｓ３５０）という。

周波数領域変換部１１０及びモデル推定部３５０における処理は第三実施形態で説明した通りである。

＜平均計算部４５０＞
平均計算部４５０は、周波数領域変換部１１０から出力された観測信号y_tuとモデル推定部１５０から出力されたパラメータの集合Θを受け取り、式（c66）及び式(c43)により、観測信号y_tuが与えられたときのノイズフリー信号r_tuの事後確率の平均μ^r|y _tuを得て出力する（ｓ４５０）。

＜強調音声生成部４６０＞
強調音声生成部４６０は、平均計算部４５０から出力された事後確率の平均μ^r|y _tuおよびモデル推定部１５０から出力されたパラメータの集合Θを受け取り、これらの値を用いて、次式により、強調音声^x_tuを生成（計算）し（ｓ４６０）、時間領域変換部４７０に出力する。

ここで、F^Hは、

により、パラメータの集合Θに含まれる予測係数行列G_kを用いて構成される行列である。なお、強調音声^x_tuは、観測信号y_tuが与えられたときのクリーン信号の事後確率の平均に相当する。

＜時間領域変換部４７０＞
最後に、時間領域変換部４７０は、強調音声^x_tuに対し、逆短時間フーリエ変換（inverse STFT）などの時間周波数変換の逆変換を適用し、時間領域におけるクリーン信号の推定値である、時間領域の信号^x_τ∈R^Mに変換し（ｓ４７０）、音声強調装置４００の出力値として出力する。

＜効果＞
以上のように、本実施形態では、非特許文献１におけるようなユークリッド距離の平方ではなく、より音響信号に適した板倉-齋藤行列ダイバージェンスを用いることにより、対数的な音響信号に、より適した処理を行うことが可能である。さらに、本実施形態は、観測信号の確率的生成モデルに基づくため、確率的生成モデルに基づく他の音声強調手法（例:音源分離、参考文献１、７参照）との統合に向いている。このような統合により、様々な音響的事象を含む実環境において、広く適用可能な新しい音声強調の枠組みを構築することができると期待されるため、大きな利点である。

＜実験結果＞
本実施形態の効果を確認するために、本実施形態の音声強調方法による、雑音および残響抑圧実験をおこなった。

観測信号としては、参考文献１３記載のREVERB challengeのReal Dataデータベースに含まれる、8チャネルの実環境収録音AMIWSJ20-Array1-*T10c020c.wav（*=1,…,8はマイクロホンの番号）から、長さ2秒の区間を切り出したものを用いた。

この観測信号は、残響時間RT60〜0.7sの残響、および雑音を含み、サンプリング周波数は16kHzであった。

比較手法としては、第四実施形態の音声強調方法において、Eステップの1〜4を、トリビアルな更新則μ^r|y _tu←y_tuおよびΦ^r|y _tt'←0で置き換えて、雑音モデルを無効にしたアルゴリズムを用いた。このアルゴリズムは、信号の空間相関をモデル化し、全てのチャネルを用いて信号分散を計算する点を除けば、参考文献１２の周波数領域NDLP法（variance-normalized delayed linear prediction)と等価である。

評価指標としては、参考文献１３にならい、ケプストラム距離（cepstrum distance:CD）、対数尤度比（log-likelihood ratio: LLR)、周波数重み付き区分信号対雑音比frequency-weighted segmental signal-to-noise ratio: FWSegSNR）、および信号対残響変調エネルギー比（speech-to-reverberation modulation energy ratio: SRMR）を用いた。これらの指標の計算のための参照信号として、ヘッドセットによる録音AMI WSJ20-Headset1 T10c020c.wav を用いた。FWSegSNR とSRMRが大きいほど、また、CDとLLRが小さいほど、音声強調性能が高いとみなす。

分析条件は、フレーム長を1024点、フレームシフトを256点（それぞれ64ms,16msに相当）、窓をハミング窓とした。予測次数はK=3、予測遅延はΔ=3、反復数は20とした。

図１０にシミュレーション結果のスペクトログラムを示す。スペクトログラムの横軸は時間（秒）、縦軸は周波数（kHz）を表す。(a)は観測信号、(b)は比較手法による処理音、(c)は第四実施形態の音声強調法で、残響と雑音の両者を考慮したモデルを学習するが、マルチチャネルウィーナーフィルタによる雑音抑圧のみを行ったもの、(d)は第四実施形態の音声強調法による、雑音と残響の同時抑圧を行ったもの、(e)は参照信号である。差異が最も明白な0-4kHzの範囲を示した。(a)の観測信号に比べ、(b)の比較手法は、残響を抑圧できており、特に1sおよび1.7s付近で顕著であるが、雑音はほとんど抑圧できていない。(c)は、(a)に対して雑音を抑圧できており、特に信号の最初の非音声区間で顕著であるが、残響はほとんど抑圧できていない。(d)は、雑音と残響を効果的に抑圧できている。

図１１に評価指標による客観評価の結果を示す。ラベル(a)-(d)は、図１０で定義した通りである。太字および斜体は、各指標について一番目および二番目に良い数値を表す。第四実施形態の音声強調方法により、CD,FWSegSNR,SRMRでは最も良い結果、LLRでも2番目に良い結果がえられた。これらの結果は、第四実施形態の音声強調方法の有効性を示している。

＜第五実施形態＞
第四実施形態と異なる部分を中心に説明する。

本実施形態では、第三実施形態に係るモデル推定装置によって推定されたモデルパラメータを用いて、雑音を除去したノイズフリー信号r_tuを推定する雑音抑圧装置５００について説明する。

図１２は第五実施形態に係る雑音抑圧装置５００の機能ブロック図、図１３はその処理フローの例を示す図である。雑音抑圧装置５００は、周波数領域変換部１１０、モデル推定部３５０、平均計算部４５０、時間領域変換部５７０を含む。

周波数領域変換部１１０及びモデル推定部３５０及び平均計算部４５０における処理は第四実施形態で説明した通りである。

＜時間領域変換部５７０＞
時間領域変換部５７０は、平均計算部４５０の出力である事後確率の平均μ^r|y _tuを受け取り、この値に対し、逆短時間フーリエ変換（inverse STFT）などの時間周波数変換の逆変換を適用し、時間領域におけるノイズフリー信号の推定値である、時間領域の信号^r_τ∈R^Mに変換し（ｓ５７０）、これを雑音抑圧装置５００の出力値として出力する。

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

＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (17)

1. 観測信号を、非定常の拡散性雑音と前記拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、モデル化するときのモデルパラメータは拡散性雑音の共分散行列を特定するモデルパラメータと、ノイズフリー信号の共分散行列を特定するモデルパラメータとを含み、
   観測信号から得られる共分散行列と、拡散性雑音の共分散行列とノイズフリー信号の共分散行列との和との板倉-齋藤距離が小さくなるように、前記モデルパラメータを更新するモデル推定部を含む、
   モデル推定装置。
2. 観測信号を、非定常の拡散性雑音と前記拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなるブロックの番号をt∈{1,…,T}とし、ブロック内のフレームの番号をu∈{1,…,U}とし、ブロックtごとの観測信号y_tuの共分散行列をΦ^y _tとし、ブロックtごとのノイズフリー信号x_tuの共分散行列をΦ^x _tとし、ブロックtごとの拡散性雑音v_tuの共分散行列をΦ^v _tとし、
   観測信号y_tから得られる共分散行列Φ^y _tと、Φ^x _t＋Φ^v _tとの板倉-齋藤距離の、すべてのブロックtについての総和が小さくなるように、前記共分散行列Φ^x _tと前記共分散行列Φ^v _tとを更新するモデル推定部を含み、
   前記ノイズフリー信号は音源信号である、
   モデル推定装置。
3. フレームごとの観測信号を非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなる時間区間をブロックとして、
   複数ブロック分の観測信号から得られる共分散行列と、複数ブロック分の拡散性雑音の共分散行列と複数ブロック分のノイズフリー信号の共分散行列との和と、の板倉-齋藤距離が小さくなるように、前記複数ブロック分の拡散性雑音の共分散行列を特定するパラメータと前記複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータとを更新するモデル推定部を含み、
   前記ノイズフリー信号は、音源から発せられる直接音と初期反射成分とからなる信号であるクリーン信号に残響が重畳された信号であり、
   前記複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータは、フレームごとの前記ノイズフリー信号に含まれる残響をK次の自己回帰過程によりモデル化したときの自己回帰過程の予測係数行列G_k(k=1,2,…,K)と、フレームごとの前記クリーン信号の共分散行列と、を含む、
   モデル推定装置。
4. 請求項１または２または３記載のモデル推定装置であって、
   前記拡散性雑音の共分散行列は、フレームごとに、行列のなす線型空間における部分空間Vに属するとしてモデル化されたものである、
   モデル推定装置。
5. 請求項２記載のモデル推定装置であって、
   観測信号をモデル化するときのモデルパラメータの集合をΘとし、前記モデル推定部は、モデルパラメータの集合Θが与えられたときの前記観測信号y_tuの集合Yと前記音源信号x_tuの集合Xとの同時分布の対数log(p(X,Y;Θ))の前記集合Xに関する期待値が最大となるように、前記共分散行列Φ^x _tと前記共分散行列Φ^v _tとを更新する、
   モデル推定装置。
6. 請求項５記載のモデル推定装置であって、
   前記モデル推定部は、
   前記観測信号y_tuと現在の前記共分散行列Φ^x _tと前記共分散行列Φ^v _tとを用いて、
   により、観測信号y_tuが与えられたときの音源信号x_tuの事後確率の平均μ^x|y _tuと共分散行列Φ^x|y _tとを更新する事後確率更新部と、
   行列・のエルミート転置を・^Hとし、観測信号を取得する際に用いるマイクロホンの個数をMとし、行列・の対角成分の和をTr[・]とし、空間的な特性を表すパラメータを表す信号コヒーレンス行列をBとし、前記平均μ^x|y _tuと前記共分散行列Φ^x|y _tとを用いて、
   を求め、
   により信号パワーφ^x _tを更新する信号パワー更新手段と、
   前記^Φ^x _tと前記信号パワーφ^x _tとを用いて、
   により、前記信号コヒーレンス行列Bを更新する信号コヒーレンス行列更新手段と、
   前記信号パワーφ^x _tと前記信号コヒーレンス行列Bとを用いて、
   により前記共分散行列Φ^x _tを更新する信号共分散行列更新手段と、
   前記部分空間Vの正規直交基底を{Q_d}_1≦d≦Dとし、m∈{1,…,M}とし、n∈{1,…,M}とし、m番目のマイクロホンとn番目のマイクロホンとの距離をL_mnとし、音速をcとし、前記観測信号y_tuに対応する周波数をfとし、sinc関数をsinc(・)とし、(m,n)成分として
   を持つ雑音コヒーレンス行列をΓとし、前記平均μ^x|y _tuと前記共分散行列Φ^x|y _tと前記観測信号y_tuとを用いて、
   を求め、前記部分空間Vが
   なる性質を持つ場合、
   により、前記共分散行列Φ^v _tを更新し、前記部分空間Vとして定数雑音コヒーレンスモデルを用いる場合、
   により前記共分散行列Φ^ｖ _tを更新する雑音共分散行列更新手段と、を含む、
   モデル推定装置。
7. 請求項３記載のモデル推定装置であって、
   観測信号をモデル化するときのモデルパラメータの集合をΘとし、前記モデル推定部は、モデルパラメータの集合Θが与えられたときの前記観測信号の集合Yと前記ノイズフリー信号の集合Ｒとの同時分布の対数log(p(R,Y;Θ))の前記集合Rに関する期待値が最大となるように、前記モデルパラメータの各々を更新する、
   モデル推定装置。
8. 請求項７記載のモデル推定装置であって、
   前記ブロックの番号をt∈{1,…,T}（Tはブロックの総数）とし、ブロック内のフレームの番号をu∈{1,…,U}とし、観測信号y_tuのブロックtごとの共分散行列をΦ^y _tとし、クリーン信号x_tuのブロックtごとの共分散行列をΦ^x _tとし、拡散性雑音v_tuのブロックtごとの共分散行列をΦ^v _tとし、行列・のエルミート転置を・^Hとし、
   とし、
   k∈{1,…,K}とし、i∈{1,…,T}とし、単位行列をIとし、所定の遅延をΔとし、j∈{1,…,T}とし、行列FはT²個のM×M行列をブロックに持つブロックテプリッツ行列であり、前記Fの(i,j)番目のブロックは
   であるとして、
   前記複数ブロック分のノイズフリー信号の共分散行列はF^-H~Φ^xF^-1である、
   モデル推定装置。
9. 請求項８記載のモデル推定装置であって、
   前記モデル推定部は、
   行列・の転置を・^Tとし、~y_u:=[y_Tu ^T…y_1u ^T]^Tとし、t'∈{1,…,T}とし、現在の前記予測係数行列G_kを用いて前記Fを更新し、現在の前記共分散行列Φ^x _tを用いて前記~Φ^xを更新し、前記Fと前記~Φ^xとを用いて、
   により、~Φ^rを更新し、現在の前記共分散行列Φ^v _tを用いて前記~Φ^vを更新し、前記~Φ^rと前記~Φ^vと前記観測信号y_tuとを用いて、
   を求め、前記~μ_u ^r|yを用いて、
   により、観測信号y_tuが与えられたときのノイズフリー信号r_tuの事後確率の平均μ_tu ^r|yを更新し、前記~Φ^rと前記~Φ^vとを用いて、
   を求め、前記~Φ^r|yを用いて、
   により、観測信号y_tuが与えられたときのブロックtおよびブロックt’におけるノイズフリー信号r_tuおよびr_t'uの事後共分散行列Φ_tt' ^r|yを更新する事後確率更新部と、
   ~G:=[G₁ ^T…G_K ^T]^Tとし、行列・の対角成分の和をTr[・]とし、空間的な特性を表すパラメータを表す信号コヒーレンス行列をB^xとし、前記平均μ_tu ^r|yと前記共分散行列Φ_tt' ^r|yとを用いて、
   を求め、前記~Gと前記Ψ_tt'とを用いて、
   を求め、前記^Φ^x _tと現在の信号コヒーレンス行列B^xとを用いて、
   により信号パワーφ^x _tを更新する信号パワー更新手段と、
   前記^Φ^x _tと前記信号パワーφ^x _tとを用いて、
   により、前記信号コヒーレンス行列B^xを更新する信号コヒーレンス行列更新手段と、
   前記信号パワーφ^x _tと前記信号コヒーレンス行列B^xとを用いて、
   により前記共分散行列Φ^x _tを更新する信号共分散行列更新手段と、
   前記共分散行列Φ^v _tは、行列のなす線型空間における部分空間Vに属するとしてモデル化されたものであり、部分空間Vの正規直交基底を{Q_d}_1≦d≦Dとし、m∈{1,…,M}とし、n∈{1,…,M}とし、m番目のマイクロホンとn番目のマイクロホンとの距離をL_mnとし、音速をcとし、前記観測信号y_tuに対応する周波数をfとし、sinc関数をsinc(・)とし、(m,n)成分として
   を持つ雑音コヒーレンス行列をB^vとし、前記平均μ_tu ^r|yと前記共分散行列Φ_tt' ^r|yとを用いて、
   を求め、前記部分空間Vが
   なる性質を持つ場合、
   により、前記共分散行列Φ^v _tを更新し、前記部分空間Vが定数雑音コヒーレンスモデルである場合、
   により、前記共分散行列Φ^v _tを更新する雑音共分散行列更新手段と、
   前記Ψ_tt'と前記信号パワーφ^x _tとを用いて、
   を求め、前記~Gを用いて、前記予測係数行列G_kを更新する予測係数行列更新手段と、を含む、
   モデル推定装置。
10. 請求項２または３のモデル推定装置により推定されたモデルパラメータと、前記観測信号y_tuとを用いて、前記観測信号y_tuが与えられたときの前記ノイズフリー信号の事後確率の平均を求める平均計算部と、
    前記事後確率の平均を時間領域の信号に変換し、時間領域における前記ノイズフリー信号の推定値を求める時間領域変換部と、を含む
    雑音抑圧装置。
11. 請求項３または８または９のいずれかのモデル推定装置により推定されたモデルパラメータと、前記観測信号y_tuとを用いて、前記観測信号y_tuが与えられたときの前記ノイズフリー信号の事後確率の平均μ^r|y _tuを求める平均計算部と、
    請求項３または８または９のいずれかのモデル推定装置により推定された前記予測係数行列G_kと前記平均計算部で求めた事後確率の平均μ_tu ^r|yとを用いて、前記観測信号y_tuが与えられたときの前記クリーン信号の事後確率の平均^x_tuを計算する強調音声生成部と、
    前記平均^x_tuを時間領域の信号に変換し、時間領域における前記クリーン信号の推定値を求める時間領域変換部と、を含む
    音声強調装置。
12. 観測信号を、非定常の拡散性雑音と前記拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、モデル化するときのモデルパラメータは拡散性雑音の共分散行列を特定するモデルパラメータと、ノイズフリー信号の共分散行列を特定するモデルパラメータとを含み、
    観測信号から得られる共分散行列と、拡散性雑音の共分散行列とノイズフリー信号の共分散行列との和との板倉-齋藤距離が小さくなるように、前記モデルパラメータを更新するモデル推定ステップを含む、
    モデル推定方法。
13. 観測信号を、非定常の拡散性雑音と前記拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなるブロックの番号をt∈{1,…,T}とし、ブロック内のフレームの番号をu∈{1,…,U}とし、ブロックtごとの観測信号y_tuの共分散行列をΦ^y _tとし、ブロックtごとのノイズフリー信号x_tuの共分散行列をΦ^x _tとし、ブロックtごとの拡散性雑音v_tuの共分散行列をΦ^v _tとし、
    観測信号y_tから得られる共分散行列Φ^y _tと、Φ^x _t＋Φ^v _tとの板倉-齋藤距離の、すべてのブロックtについての総和が小さくなるように、前記共分散行列Φ^x _tと前記共分散行列Φ^v _tとを更新するモデル推定ステップを含み、
    前記ノイズフリー信号は音源信号である、
    モデル推定方法。
14. フレームごとの観測信号を非定常の拡散性雑音と拡散性雑音を含まない信号であるノイズフリー信号との和でモデル化するものとし、１つ以上のフレームからなる時間区間をブロックとして、
    複数ブロック分の観測信号から得られる共分散行列と、複数ブロック分の拡散性雑音の共分散行列と複数ブロック分のノイズフリー信号の共分散行列との和と、の板倉-齋藤距離が小さくなるように、前記複数ブロック分の拡散性雑音の共分散行列を特定するパラメータと前記複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータとを更新するモデル推定ステップを含み、
    前記ノイズフリー信号は、音源から発せられる直接音と初期反射成分とからなる信号であるクリーン信号に残響が重畳された信号であり、
    前記複数ブロック分のノイズフリー信号の共分散行列を特定するパラメータは、フレームごとの前記ノイズフリー信号に含まれる残響をK次の自己回帰過程によりモデル化したときの自己回帰過程の予測係数行列G_k(k=1,2,…,K)と、フレームごとの前記クリーン信号の共分散行列と、を含む、
    モデル推定方法。
15. 請求項１３のモデル推定方法により推定されたモデルパラメータと、前記観測信号y_tuとを用いて、前記観測信号y_tuが与えられたときの前記ノイズフリー信号の事後確率の平均を求める平均計算ステップと、
    前記事後確率の平均を時間領域の信号に変換し、時間領域における前記ノイズフリー信号の推定値を求める時間領域変換ステップと、を含む
    雑音抑圧方法。
16. 請求項１４のモデル推定方法により推定されたモデルパラメータと、前記観測信号y_tuとを用いて、前記観測信号y_tuが与えられたときの前記ノイズフリー信号の事後確率の平均μ^r|y _tuを求める平均計算ステップと、
    請求項１４のモデル推定方法により推定された前記予測係数行列G_kと前記平均計算ステップで求めた事後確率の平均μ_tu ^r|yとを用いて、前記クリーン信号の事後確率の平均^x_tuを計算する強調音声生成ステップと、
    前記平均^x_tuを時間領域の信号に変換し、時間領域における前記クリーン信号の推定値を求める時間領域変換ステップと、を含む
    音声強調方法。
17. 請求項１〜９の何れかのモデル推定装置、または、請求項１０の雑音抑圧装置、または、請求項１１の音声強調装置として、コンピュータを機能させるためのプログラム。