JPWO2009110574A1 - 信号強調装置、その方法、プログラム及び記録媒体 - Google Patents

Abstract

観測信号に含まれる残響の推定値を算出する線形畳み込み演算の回帰係数を含む残響パラメータ推定値と、信号源のパワースペクトルを特定する線形予測係数と予測残差パワーとの推定値を含む信号源パラメータの推定値と、雑音のパワースペクトルの推定値を含む雑音パラメータ推定値とを含むパラメータ推定値の初期値を設定する。その後、最尤推定によって、残響パラメータ推定値及び雑音パラメータ推定値の少なくとも一部を更新する処理と、信号源パラメータ推定値を更新する処理とを、所定の終了条件を満たすまで交互に繰り返す。

Description

本発明は、観測信号中の加法性歪みと乗法性歪みとを抑圧して源信号を強調する技術に関する。

源信号に加法性歪みや乗法性歪みが重畳された観測信号に対し、加法性歪み又は乗法性歪みを抑圧する処理を行い、源信号を強調する信号強調技術がある。まず、信号が音声信号である場合での一般的な音声信号強調技術を説明する。この場合、加法性歪みは室内に存在する雑音に、乗法性歪みは残響に対応する。

図１は、信号強調装置の一般的な構成を示すブロック図である。
まず、マイクロホン等のセンサや音声ファイル等から取得され、標本化・量子化された時間領域の観測音声の波形信号が帯域分割部に入力される。これらの時間領域の観測信号は、帯域分割部において、周波数帯域ごとの狭帯域信号に分割される。すなわち、時間領域の観測信号が時間周波数領域の観測信号に変換される。以下では、周波数帯域ごとに分割された観測信号の集合を観測信号の複素スペクトログラムと呼ぶ。なお、帯域分割部は、短時間フーリエ変換やポリフェーズフィルタバンク等の従来技術によってこの処理を実行する。ただし、この帯域分割を実施せずに、時間領域の観測信号を直接用いて源信号の強調処理を行う方法もある。また、明細書では、信号を表現する領域を明記していない場合、時間周波数領域であると解釈する。

次に、パラメータ推定部において、観測信号の複素スペクトログラムから、観測信号を特徴づける何らかのパラメータが推定される。パラメータの例は、源信号あるいは雑音のパワースペクトルを記述する全極モデルのパラメータや、室内伝達系を記述する自己回帰モデルの回帰係数などである。

そして、源信号推定部において、観測信号の複素スペクトログラムと上記パラメータの推定値とを用い、源信号の複素スペクトログラムの推定値が計算される。最後に、帯域合成部において、源信号の複素スペクトログラムの推定値から時間領域の源信号の推定値が合成される。なお、帯域合成部の処理は帯域分割部の処理に対応する。すなわち、帯域分割部が短時間フーリエ変換を実行するのであれば帯域合成部はオーバーラップ加算合成を行い、帯域分割部がポリフェーズフィルタバンク分析を実行するのであれば帯域合成部はポリフェーズフィルタバンク合成を行う。帯域分割部が省略された場合には、帯域合成部も省略される。

従来の音声信号強調技術は、源信号以外に雑音のみが存在する環境を対象とするものと（例えば、非特許文献１参照）、源信号以外に残響のみが存在する環境を対象とするものに大別される（例えば、非特許文献２参照）。前者は、源信号以外に雑音を含む観測信号から雑音を抑圧する。後者は、源信号以外に残響を含む観測信号から残響を抑圧する。以下に非特許文献１，２でそれぞれ提案されている音声信号強調技術について説明する。なお、以下の説明において、テキスト中で使用する記号「＾」「^〜」等は、文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直後に記載する。

＜非特許文献１の雑音抑圧技術＞
非特許文献１には、源信号に雑音が加算された観測信号から雑音を抑圧する雑音抑圧技術が提案されている。以下に非特許文献１に開示された各処理部の処理を説明する。

非特許文献１の帯域分割部は、観測された観測信号を短時間フーリエ変換によって周波数帯域ごとの狭帯域信号に分割する。また、非特許文献１のパラメータ推定部は、観測信号、すなわち源信号に雑音が重畳された信号を特徴づけるパラメータとして、源信号の全極モデルの信号源パラメータ_sΘ及び雑音モデルの雑音パラメータ_dΘを推定する。

非特許文献１の例では、まず、源信号が存在しない時間区間の観測信号を用い、雑音パラメータの真値_dΘ^〜が計算される（ステップＳ１０１）。次に、信号源パラメータ推定値の初期値_sΘ^⁽⁰⁾が設定される（ステップＳ１０２）。また、繰り返し回数を示すインデックスｉが０に設定される（ステップＳ１０３）。

その後、信号源パラメータの推定値_sΘ^⁽ⁱ⁾と雑音パラメータの真値_dΘ^〜とを用い、信号源パラメータの推定値_sΘ^⁽ⁱ⁾と雑音パラメータの真値_dΘ^〜の組合せと観測信号の複素スペクトログラムＹが与えられた場合における源信号の複素スペクトログラムＳの条件付事後分布p(S|Y,_sΘ^⁽ⁱ⁾,_dΘ^〜)を算出する（ステップＳ１０４）。次に、条件付事後分布p(S|Y,_sΘ^⁽ⁱ⁾,_dΘ^〜)を用い、信号源パラメータの推定値_sΘ^⁽ⁱ⁾を_sΘ^⁽ⁱ⁺¹⁾に更新する（ステップＳ１０５）。そして、終了条件を満たすまで（ステップＳ１０６）、ｉを１ずつ増加させながら（ステップＳ１０７）、ステップＳ１０４とＳ１０５との処理を繰り返し、所定の終了条件が満たされた時点における信号源パラメータの推定値_sΘ^〜(i+1)を信号源パラメータの最終推定値_sΘ^として出力する（ステップＳ１０８）。

その後、源信号推定部が、パラメータ推定部で計算されたパラメータ_dΘ^〜と_sΘ^を用い、Wienerフィルタを用いて、源信号の複素スペクトログラムの推定値を求め、帯域合成部が、オーバーラップ加算合成によって、当該複素スペクトログラムの推定値を時間領域の源信号の推定値に変換する。

＜非特許文献２の残響抑圧技術＞
非特許文献２には、源信号に残響が重畳された観測信号から残響を抑圧する残響抑圧技術が提案されている。以下に非特許文献２に開示された各処理部の処理を説明する。

非特許文献２の残響抑圧技術では、帯域分割処理は実施されない。したがって、非特許文献２のパラメータ推定部及び源信号推定部は、時間領域の観測信号を直接処理する。このパラメータ推定部は、観測信号、すなわち源信号に残響が重畳された信号を特徴づけるパラメータとして、信号源パラメータ_sΘ及び残響パラメータ_gΘを推定する。なお、非特許文献２の残響パラメータは、源信号に残響のみが重畳された時間領域の観測信号に適用され、観測信号に重畳された残響を算出する線形フィルタの回帰係数である。

非特許文献２の例では、まず、残響パラメータの推定値の初期値_gΘ^⁽⁰⁾を設定する（ステップＳ１１１）。また、繰り返し回数を示すインデックスｉを０に設定する（ステップＳ１１２）。

その後、残響パラメータの推定値_gΘ^⁽ⁱ⁾を用い、信号源パラメータの推定値を_sΘ^⁽ⁱ⁺¹⁾に更新する（ステップＳ１１３）。次に、更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾を用い、残響パラメータの推定値を_gΘ^⁽ⁱ⁺¹⁾ に更新する（ステップＳ１１４）。そして、所定の終了条件を満たすまで（ステップＳ１１５）、ｉを１ずつ増加させながら（ステップＳ１１６）、ステップＳ１１３とＳ１１４との処理を繰り返し、所定の終了条件が満たされた時点における信号源パラメータの推定値_sΘ^〜(i+1)を信号源パラメータの最終的な推定値_sΘ^とし、最終的な残響パラメータ推定値_gΘ^⁽ⁱ⁺¹⁾を最終的な推定値_gΘ^として出力する（ステップＳ１１７）。

その後、源信号推定部が、パラメータ推定部で計算された残響パラメータの最終的な推定値_gΘ^を用いて生成した線形フィルタを観測信号に畳み込み、源信号成分を強調した音声信号を算出して出力する。その後、源信号推定部が、パラメータ推定部で計算された残響パラメータの最終的な推定値_gΘ^を用いて生成した線形フィルタを観測信号に畳み込むことで観測信号に含まれる残響を推定し、それを観測信号から減算することで、残響が抑圧された信号を算出して出力する。
Lim, J. S. and Oppenheim, A. V. , "All-pole modeling of degraded speech," IEEE Trans. Acoust. Speech, Signal Process., Vol. 26, No. 3, pp.197-210 (1978). Yoshioka, T., Hikichi, T. and Miyoshi, M., "Dereverberation by Using Time-Variant Nature of Speech Production System, EURASIP J. Advances in Signal Process., Vol. 2007, (2007), Article ID 65698, 15 pages, doi:10.1155/2007/65698.

しかし、雑音と残響がともに存在する環境を対象とした信号強調技術はこれまで存在しなかった。
雑音と残響がともに存在する環境においてＭ（Ｍ≧１）個のセンサ１０００−１〜Ｍで観測された観測信号は、図２に示す系によって生成されたものであるといえる。すなわち、まず、話者などの信号源１０１０から発せられた、雑音や残響を含まない信号（「源信号」と呼ぶ）に対し、残響重畳系（室内伝達系）によって各室内インパルス応答が畳み込まれることで残響が付加される。さらに、残響が付加された信号（「残響重畳信号」と呼ぶ）に対し、雑音重畳系によって雑音が加算される。これにより、雑音と残響を含む信号（「雑音残響重畳信号」と呼ぶ）が生成され、各センサで観測される。

前述の通り、従来の残響抑圧技術は、残響重畳信号が与えられたときに残響パラメータと信号源パラメータを推定した後、推定された残響パラメータに基づいて源信号を回復する。ゆえに図２の系において残響抑圧処理を行うためには、雑音抑圧処理によって雑音残響重畳信号から予め雑音を抑圧して残響重畳信号を求めておかなければならない。一方図２の系において雑音残響重畳信号から効果的に雑音を抑圧するためには、残響重畳信号の特性が既知であることが望ましい。ところが残響重畳信号の特性は、源信号の特性（すなわち、源信号の信号源パラメータ）と室内伝達系（すなわち、残響パラメータ）によって規定されるから、これは残響抑圧処理によって求められるものである。したがって、図２の系において源信号を効果的に強調するためには、雑音抑圧処理と残響抑圧処理を協調して動作させる必要がある。

また、従来の雑音抑圧技術は、源信号に雑音のみ加算された観測信号から雑音を抑圧するものである。そのため、従来の雑音抑圧技術を、雑音と残響を含む雑音残響重畳信号から雑音を抑圧するという上記の雑音抑圧処理にそのまま適用しても、精度よい雑音抑圧は期待できない。また、雑音抑圧処理と残響抑圧処理を単純に結合させるのではなく協調的に動作させることが必要であると述べたが、これをいかにして行うかは自明でない。

このような問題は、音声信号を対象にする場合だけではなく、その他の音響信号、超音波信号その他の信号を対象とする場合にも共通するものである。すなわち、信号源から発せられた加法性歪みや乗法性歪みを含まない信号に、線形畳み込み系によって乗法性歪みが付加され、それによって生成された信号に対し、さらに加法性歪みが加算されて生成された信号から、加法性歪みや乗法性歪みを抑圧し、元の信号を強調する場合一般に共通する問題である。本明細書では、音声信号を対象にする場合との関係を明確にするため、信号源から発せられた加法性歪みや乗法性歪みを含まない信号を「源信号」、源信号に乗法性歪みが付加されて生成された信号を「残響重畳信号」、残響重畳信号に加法性歪みが付加されて生成された信号を「雑音残響重畳信号」、乗法性歪みを付加する線形畳み込み系を「室内伝達系」、加法性歪みを「雑音」、乗法性歪みを「残響」と呼ぶことにする。

本発明のパラメータ推定部では、まず、観測された時間領域信号から変換された時間周波数領域の観測信号を記録部に格納し、初期化部において、観測信号に含まれる残響の推定値を算出する線形畳み込み演算の回帰係数を含む残響パラメータ推定値と、源信号のパワースペクトルを特定する線形予測係数と予測残差パワーの推定値を含む信号源パラメータ推定値と、雑音のパワースペクトルの推定値を含む雑音パラメータ推定値と、を含むパラメータ推定値の初期値を設定する。

次に、観測信号とパラメータ推定値とを第１更新部に入力し、当該第１更新部において、残響パラメータ推定値および雑音パラメータ推定値の少なくとも一部の更新処理、あるいは信号源パラメータ推定値の更新処理、のいずれか一方を行う。更新処理は、パラメータ推定値に関する対数尤度関数の値が増加するように実行される。

また、第１更新部で得られたパラメータ推定値の更新値の少なくとも一部を第２更新部に入力し、第２更新部において、残響パラメータ推定値および雑音パラメータ推定値の少なくとも一部の更新処理、あるいは信号源パラメータ推定値の更新処理のうち、第１更新部で実行されなかったものを実行する。更新処理は、パラメータ推定値の更新値に関する対数尤度関数の値が増加するように実行される。

そして、終了条件判定部において、終了条件が満たされるか否かを判定し、終了条件が満たされない場合、第１更新部と第２更新部の処理が再び実行される。

以上のように、本発明のパラメータ推定部では、第１更新部におけるパラメータの推定値の更新処理と、第２更新部におけるパラメータの推定値の更新処理を、互いに依存させながら繰り返して実行する。これににより、雑音と残響がともに存在する環境における観測信号から、雑音と残響を精度よく抑圧し、源信号を強調することができる。

図１は、音声信号強調装置の一般的な構成を示すブロック図である 図２は、源信号に雑音や残響が付加される系を説明するための図である。 図３は、第１実施形態の信号強調装置の構成を示すブロック図である。 図４は、源信号推定部の詳細構成を示すブロック図である。 図５は、第１実施形態の信号強調方法を説明するためのフローチャートである。 図６は、第２実施形態の信号強調装置の構成を示すブロック図である。 図７は、源信号推定部の詳細構成を示すブロック図である。 図８は、第１実施形態の信号強調方法を説明するためのフローチャートである。 図９は、第３実施形態の信号強調装置の機能構成例を示すブロック図である。 図１０は、第３実施形態の処理を説明するためのフローチャートである 図１１は、第４実施形態のパラメータ推定部の機能構成例を示すブロック図である。 図１２は、第４実施形態のパラメータ推定処理を説明するためのフローチャートである。

以下、図面を参照して本発明の実施の形態を説明する。
まず、本実施形態のパラメータ推定部について述べる。本実施形態のパラメータは、残響パラメータと、信号源パラメータと、雑音パラメータとを含む。残響パラメータは、少なくとも、室内伝達系を多チャンネル自己回帰系としてモデル化したときの回帰行列を含む。なお、この回帰行列からなる多入力多出力インパルス応答を残響重畳信号に畳み込むと、残響重畳信号に含まれる残響が算出される。信号源パラメータは、少なくとも、源信号の短時間パワースペクトル密度を特徴づける線形予測係数と予測残差パワーとを含む。雑音パラメータは、少なくとも、雑音の短時間パワークロススペクトル行列を含む。本実施形態のパラメータ推定部は、残響パラメータと信号源パラメータと雑音パラメータを、ＥＣＭアルゴリズム等のＥＭアルゴリズムの変種を用いて、最尤推定する。

具体的には、本実施形態のパラメータ推定部は、例えば、以下のように表現される。本実施形態のパラメータは、２つの群に分類される。第１パラメータ群は、少なくとも、残響パラメータを含む。第２パラメータ群は、少なくとも、信号源パラメータを含む。雑音パラメータは、第１パラメータ群、第２パラメータ群のいずれに含まれてもよいが、本実施形態では第１パラメータ群に含まれることとする。

まず、観測信号を記憶部に格納する。
初期化部は、第１パラメータ群のパラメータの推定値と、第２パラメータ群のパラメータの推定値とを初期化する。
次に、観測信号と、第１パラメータ群のパラメータの推定値と、第２パラメータ群のパラメータの推定値とが、第１更新部に入力される。第１更新部は、第１パラメータ群と第２パラメータ群のいずれか一方のパラメータ群のパラメータの推定値を固定し、残る一方のパラメータ群のパラメータのうち、少なくとも一部のパラメータの推定値を更新する。第１更新部は、パラメータの推定値に関する対数尤度関数の値が大きくなるように、パラメータの推定値を更新する。

次に、観測信号と、第１パラメータ群のパラメータの推定値と、第２パラメータ群のパラメータの推定値のうちの少なくとも一部が第２更新部に入力される。第２更新部は、第１更新部で更新されたパラメータ群のパラメータの推定値を固定し、第１更新部で固定されたパラメータ群のパラメータのうち、少なくとも一部のパラメータの推定値を更新する。第２更新部は、パラメータの推定値に関する対数尤度関数の値が大きくなるように、パラメータの推定値を更新する。

終了判定条件部は、所定の終了条件が満たされているか否かを判定する。終了条件が満たされていない場合、第１更新部の処理に戻る。終了条件が満たされている場合、その時点のおけるパラメータの推定値を出力する。

〔第１実施形態〕
＜本実施形態のパラメータ推定処理の概要＞
まず、本実施形態のパラメータ推定処理の概要を説明する。
[観測信号記憶処理]
まず、観測信号記憶処理によって、観測信号が記憶部に格納される。
[初期化処理]
次に、初期化処理によって、第１パラメータ群のパラメータの推定値と、第２パラメータ群のパラメータの推定値とが初期化される。

[第１更新処理]
本実施形態の第1更新処理では、第1パラメータ群、すなわち残響パラメータの推定値が固定された状態で、第２パラメータ群、すなわち信号源パラメータの推定値が更新される。本実施形態の第１更新処理は、具体的には、雑音抑圧処理と、信号源パラメータの更新処理とを含む。

《雑音抑圧処理》
雑音抑圧処理では、観測信号とパラメータの推定値を用いて、残響重畳信号の条件付事後分布ｐ（残響重畳信号｜観測信号，パラメータの推定値）を特徴づける複素正規分布の平均と共分散行列が算出される。

この処理は、観測信号から雑音を含まない残響重畳信号の条件付事後分布を求めるという点において、観測信号に含まれる雑音を抑圧していると解釈できる。この雑音抑圧処理は、残響パラメータの推定値と信号源パラメータの推定値を用いて実行されることに注意されたい。このことは、残響の特性が考慮されながら雑音が抑圧されることを意味する。これによって、残響環境において、雑音抑圧を精度よく実施できる。

《信号源パラメータ推定値の更新処理》
信号源パラメータ推定値の更新処理では、残響パラメータの推定値と残響重畳信号の条件付事後分布の平均と共分散行列を用いて、信号源パラメータの推定値が更新される。信号源パラメータの推定値は、パラメータの推定値に関する補助関数の値が最大になるように、更新される。

補助関数は、観測信号と残響重畳信号を所与とした場合のパラメータの推定値に関する対数尤度関数を、残響重畳信号の条件付事後分布ｐ（残響重畳信号｜観測信号，パラメータ推定値）で重み付けした関数を、残響重畳信号について積分して得られる関数である。この重み付け積分により、雑音抑圧処理で算出される残響重畳信号の不確かさを考慮しながら、信号源パラメータの推定値を更新することが可能になっている。

[第２更新処理]
本実施形態の第２更新処理では、第２パラメータ群、すなわち信号源パラメータの推定値が固定された状態で、第１パラメータ群、すなわち残響パラメータの推定値が更新される。残響パラメータの推定値は、パラメータの推定値に関する補助関数の値が最大になるように、更新される。

[終了条件判定処理]
終了条件判定処理では、所定の終了条件が満たされているか否かが判定される。終了条件が満たされていない場合、第１更新処理に戻る。終了条件が満たされている場合、その時点におけるパラメータの推定値を出力する。
以上で述べた処理において、残響重畳信号の条件付事後分布の共分散行列は、雑音の分散に対して単調増加する。すなわち、雑音のレベルが大きいほど、残響重畳信号の条件付事後分布の共分散行列も大きくなる。このことは、本実施形態が、雑音抑圧処理で求められる残響重畳信号の不確かさを妥当な方法で評価していることを示している。

＜本実施形態の原理＞
次に、本実施形態の原理を説明する。
本実施形態は統計的推定の方法論に基づく。まず、信号源パラメータ_sΘ、残響パラメータ_gΘ、及び雑音パラメータ_dΘが規定される必要がある。また、すべてのパラメータの集合がΘ={_sΘ, _gΘ, _dΘ}と表現される。次に、規定したパラメータΘが、観測信号である雑音残響重畳信号の集合Ｙに対応づけられなければならない。なお、雑音残響重畳信号の集合Ｙは、所定の観測区間に属する雑音残響重畳信号の集合である。後述するように、本実施形態の雑音残響重畳信号の集合Ｙは、雑音残響重畳信号の複素スペクトログラムである。

本実施形態では、パラメータΘが与えられた場合における雑音残響重畳信号の集合Ｙの確率密度関数p(Y|Θ)が定式化され、この対応づけが行われる。この定式化により、雑音残響重畳信号の集合Ｙは、未知のパラメータの真値Θ^〜={_sΘ^〜, _gΘ^〜, _dΘ^〜}を前提とした確率密度関数p(Y|Θ^〜)で表される確率分布をとる信号であると捉えることができる。

また、本実施形態では、観測信号である雑音残響重畳信号の集合Ｙからパラメータの真値Θ^〜が最尤推定される。すなわち、雑音残響重畳信号の集合Ｙが観測されたときの尤度関数p(Y|Θ^〜)を最大化するパラメータの値Θ^={_sΘ^, _gΘ^, _dΘ^〜}が求められ、これがパラメータの真値Θ^〜の最終的な推定値とされる。なお、雑音パラメータ_dΘは、源信号が存在しない区間から独立に推定され、その推定値が雑音パラメータの真値_dΘ^〜であると仮定される。したがって、最尤推定法によって推定される値は、信号源パラメータの真値_sΘ^〜、及び残響パラメータの真値_gΘ^〜である。

ところが実際には、確率密度関数p(Y|Θ^〜)を最大化する_sΘ^〜と_gΘ^〜を同時に直接求めることはできない。そこで、本実施形態ではＥＣＭ（Expectation-Conditional Maximization）アルゴリズムが適用される。すなわち、観測信号である雑音残響重畳信号の集合Ｙを用い、雑音残響重畳信号の集合Ｙとパラメータの推定値Θ^との組合せを前提条件とした残響重畳信号の集合Ｘの条件付事後分布p(X|Y,Θ^)の算出処理（Ｅ−ｓｔｅｐ）と、信号源パラメータの推定値_sΘ^の更新処理（ＣＭ−ｓｔｅｐ１）と、残響パラメータの推定値_gΘ^の更新処理（ＣＭ−ｓｔｅｐ２）とが代わる代わる繰り返し実行されて各推定値が更新され、所定の終了条件を充足した時点での各推定値が真値の推定値（最終推定値）とされる。なお、残響重畳信号の集合Ｘは、所定の観測区間に属する残響重畳信号の集合である。後述するように、本実施形態の残響重畳信号の集合Ｘは、残響重畳信号の複素スペクトログラムである。

［観測信号（雑音残響重畳信号）の統計的モデル］
最初になすべきことは、パラメータΘが与えられた場合における雑音残響重畳信号の集合のＹの確率密度関数p(Y|Θ)を定義することである。そのために、観測信号（雑音残響重畳信号）の集合Ｙの統計的モデルが仮定される。本実施形態では、以下に述べる源信号の全極モデル、室内伝達系の自己回帰モデル及び雑音のモデルが仮定される。

なお、以下では、すべての信号が周波数領域で定義される複素スペクトログラムに変換されているものとする。また、複素スペクトログラムのフレーム数をＴ（定数）とし、周波数帯域数をＮ（定数）とする。なお、各説明では短時間フーリエ変換を想定した用語を用いるが、信号の周波数領域への変換には、ポリフェーズフィルタバンク等、帯域幅が一定であるような任意の時間周波数解析方法を用いることができる。

《源信号のモデル》
まず、源信号の全極モデルについて述べる。t(0≦t≦T-1)番目のフレーム、w(0≦w≦N-1)番目の周波数帯域における源信号の離散フーリエ係数（複素数）をS_t,wとおく。なお、t(0≦t≦T-1)は各フレームに対応するインデックスであり、w(0≦w≦N-1)は各周波数帯域に対応するインデックスである。
S_t,wは以下の条件を満たすと仮定される。
１．ω∈{‐π,π}を角周波数として、ｔ番目のフレームにおける源信号のパワースペクトル密度_sλ_t(ω)は、以下のようなＰ次（Ｐ≧１）の全極型スペクトル密度で表される。

なお、{a_t,1,…,a_t,P}と_sσ_t ²とは、それぞれ、源信号を線形予測分析した場合における線形予測係数と予測残差パワーである。また、ｚはｚ変換における複素変数であり、eはネイピア数である。また、jは虚数単位である。よって、信号源パラメータ_sΘは、_sΘ={a_t,1 ,..., a_t,P, _sσ_t ²}_0≦t≦T-1と定義される。ただし、{m_α}_0≦α≦M-1は、m₀, m₁ ,..., m_M-1のＭ個の要素からなる集合を表す。
２．S_t,wは、以下のように、平均０、分散_sλ_t(2πw/N)の複素正規分布にしたがう。

ただし、N_C{x;μ,Σ}は、次式で定義される平均μ、共分散行列Σの複素正規分布にしたがうζ次元確率変数ｘの確率密度関数である。なお、α^Hは、αの複素共役転置（エルミート共役）を意味する。

ただし、|Σ|はΣの行列式を示す。ここで、ζ＝１として式(4)を式(3)に代入するとS_t,wの確率密度関数は次式で表される。

３．(t,w)≠(t',w')ならば、S_t,wとS_t',w'は統計的に独立である。
《室内伝達系のモデル》
次に、室内伝達系のモデルについて述べる。t(0≦t≦T-1)番目のフレーム、w(0≦w≦N-1)番目の周波数帯域における残響重畳信号の離散フーリエ係数をX_t,wとおく。室内伝達系は各周波数帯域において自己回帰系として表現できると仮定される。すなわち、ｗ番目の周波数帯域における自己回帰系の回帰係数をg_1,w, ..., g_Kw,wとおくと、残響重畳信号の離散フーリエ係数X_t,wは次式により生成される。ただし、g_k,w ^*はg_k,wの複素共役値である。

_gΘ={{g_k.w}_1≦k≦Kw}_0≦w≦N-1が残響パラメータ_gΘと定義される。この残響パラメータ_gΘは、次式に示すように、源信号に残響のみが付加された残響重畳信号に適用されて残響重畳信号に含まれる残響を算出する用途に供される。

《雑音のモデル》
次に、雑音のモデルについて述べる。本実施形態では、t(0≦t≦T-1)番目のフレーム、w(0≦w≦N-1)番目の周波数帯域における、雑音と雑音残響重畳信号との離散フーリエ係数がそれぞれD_t,w，Y_t,wとされる。Y_t,wは残響重畳信号X_t,wに雑音D_t,wを加算したものである。
Y_t,w = X_t,w + D_t,w (7)
また、D_t,wが次に述べる条件を満たすと仮定される。
１．雑音は定常であり、そのパワースペクトル密度を_dλ(ω)として（定常であるためフレーム番号ｔには依存しない）、D_t,wは平均０、分散_dλ(2πw/N)の複素正規分布に従う。

ただし、雑音パラメータ_dΘは、_dΘ={_dλ(2πw/N)}_{0≦ｗ≦N-1}と定義される雑音を特徴づけるパラメータである。
２．(t, w)≠(t', w')ならば、D_t,wとD_t',w'とは統計的に独立である。
３．任意の(t, w, t', w')について、S_t,wとD_t',w'とは統計的に独立である。

《雑音残響重畳信号の確率密度関数》
以上の仮定に基づき、雑音残響重畳信号の確率密度関数が定式化される。
本実施形態では、源信号、残響重畳信号及び雑音残響重畳信号の各複素スペクトログラム（源信号、残響重畳信号及び雑音残響重畳信号の各集合に相当）がそれぞれＳ、Ｘ及びＹと表現される。すなわち、
S={S_t,w}_{0≦t≦T-1, 0≦w≦N-1} (9)
X={X_t,w}_{0≦t≦T-1, 0≦w≦N-1} (10)
Y={Y_t,w}_{0≦t≦T-1, 0≦w≦N-1} (11)
と表現される。なお、{m_α,β}_{0≦α≦T-1, 0≦β≦N-1}は、m_0,0 ,..., m_T-1,N-1のT・N個の要素からなる集合を表す。
具体的には、雑音残響重畳信号の複素スペクトログラムＹの確率密度関数（観測信号の集合Ｙが与えられたときのパラメータΘに関する尤度関数に相当）は次のように書ける。

ただし、p(Y,X|Θ)は、以上の仮定に基づいて次式のように書ける。

以上で、パラメータΘ={_sΘ, _gΘ, _dΘ} を用いて雑音残響重畳信号の複素スペクトログラムの確率密度関数p(Y|Θ)が定式化された。

［信号源パラメータ及び残響パラメータの最尤推定］
前述のように、本実施形態では、観測された雑音残響重畳信号の複素スペクトログラムＹから、未知のパラメータの真値Θ^〜が、最尤推定法によって推定される。すなわち、雑音残響重畳信号の集合Ｙが与えられた場合におけるパラメータΘを変数とした尤度関数p(Y|Θ)を最大化するΘが、真値Θ^〜の推定値となる。ただし、本実施形態では、雑音パラメータの真値_dΘ^〜が源信号の存在しない区間から予め独立に推定され、既知となっている為Θ^={_sΘ^, _gΘ^, _dΘ^〜}であり、_sΘ^と_gΘ^が求められることになる。

また、尤度関数p(Y|Θ)を最大化する_sΘ^と_gΘ^を同時に直接求めることはできないから、ＥＣＭアルゴリズムを用いてこれらが計算される。ＥＣＭアルゴリズムの処理の流れを以下に示す。以下の処理では、Ｅ−ｓｔｅｐ、ＣＭ−ｓｔｅｐ１、ＣＭ−ｓｔｅｐ２の３つの処理が代わる代わる繰り返し実行される。そこで、ｉ回目の繰り返しにおけるパラメータの推定値を上付きの添え字(i)を用いて示す。明確さを期するために述べると、Θ^〜，Θ^，Θ^⁽ⁱ⁾はそれぞれ次のように定義される。

《ＥＣＭアルゴリズム》
１．パラメータの推定値の初期値Θ^⁽⁰⁾が決められる。また、繰り返し回数を示すインデックスｉが０にされる。
２．Ｅ−ｓｔｅｐ（雑音抑圧処理）
残響重畳信号の条件付事後分布p(X|Y, Θ^⁽ⁱ⁾)が計算される。
３．ＣＭ−ｓｔｅｐ１（信号源パラメータ推定値の更新処理）
補助関数Q(Θ|Θ^⁽ⁱ⁾)が次式により定義される。

このとき、次の手続きにより、信号源パラメータの推定値が_SΘ^⁽ⁱ⁾から_SΘ^⁽ⁱ⁺¹⁾に更新される。

すなわち、残響パラメータの推定値_gΘ^⁽ⁱ⁾が固定された条件下で補助関数Q(Θ|Θ^⁽ⁱ⁾)を最大化する_SΘ^⁽ⁱ⁺¹⁾が、更新された信号源パラメータの推定値とされる。
４．ＣＭ−ｓｔｅｐ２（残響パラメータ推定値の更新処理）
次の手続きにより、残響パラメータの推定値が更新される。

すなわち、信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾が固定された条件下で補助関数Q(Θ|Θ^⁽ⁱ⁾)を最大化する_gΘ^⁽ⁱ⁺¹⁾が、残響パラメータの更新された推定値とされる。
５．終了条件判定
所定の終了条件を満たしているならば_sΘ^=_sΘ^⁽ⁱ⁺¹⁾，_gΘ^=_gΘ^⁽ⁱ⁺¹⁾として終了。そうでなければ、ｉを１だけ漸増させて「２．Ｅ−ｓｔｅｐ」へ戻る。

《各ｓｔｅｐの計算方法》
以下では、Ｅ−ｓｔｅｐ、ＣＭ−ｓｔｅｐ１及びＣＭ−ｓｔｅｐ２の各計算方法を説明する。
１．Ｅ−ｓｔｅｐの計算方法
源信号、残響重畳信号、雑音残響重畳信号のｗ番目の周波数帯域の離散フーリエ係数系列を、それぞれまとめて次のように表す。

源信号の複素スペクトログラムＳ、残響重畳信号の複素スペクトログラムＸ及び雑音残響重畳信号の複素スペクトログラムＹは、それぞれ、S_w, X_w, Y_wの全周波数帯域（0≦w≦N-1）にわたる集合と等価となる。
式(24)の残響重畳信号の条件付事後分布p(X|Y, Θ^⁽ⁱ⁾)は、次式に示すように周波数帯域wごとに独立な複数の複素正規分布によって表現できる。

なお、平均μ_w(Θ^⁽ⁱ⁾,Y)と共分散行列Σ_w(Θ^⁽ⁱ⁾)は次式で与えられる。

式(29),(30)に現れる各変数はそれぞれ以下のように定義される。なお、式(31)の空欄部分の各要素は０である。

なお、前述のように、雑音が定常であると仮定されているため、
_dλ_T-1 ^〜(2πw/N)= _dλ_T-2 ^〜(2πw/N)=...=_dλ₀ ^〜(2πw/N)=_dλ^〜(2πw/N)
である。また、diag{α_１,...,α_β}は、任意のスカラー値α_１,...,α_βを対角要素とする対角行列である。

式(28)で示されるように、この残響重畳信号の条件付事後分布p(X|Y, Θ^ ⁽ⁱ⁾)は、信号源パラメータ及び残響パラメータ、及び雑音パラメータに基づいて算出される。さらに、式(30),(34)に示すように、この残響重畳信号の集合Xの条件付事後分布p(X|Y, Θ^ ⁽ⁱ⁾)の共分散行列のスケールは、雑音のパワースペクトル（雑音の確率分布を示す複素正規分布の分散）に対して単調増加する値となっている。この場合、雑音のレベルが大きかった場合には残響重畳信号の集合Xの条件付事後分布の共分散行列のスケールも大きくなり、逆に雑音のレベルが小さかった場合には残響重畳信号の集合Xの条件付事後分布の共分散行列のスケールも小さくなる。この振る舞いは極めて自然である。この特徴により、雑音と残響とが存在する環境でのパラメータ推定精度を向上させることができる。

また、後の処理のために、μ_m,w ⁽ⁱ⁾を平均μ_w(Θ^⁽ⁱ⁾,Y)のＴ−ｍ番目の要素とし、μ_m:n,w ⁽ⁱ⁾（m≧n）を平均μ_w(Θ^⁽ⁱ⁾,Y)のＴ−ｍ番目からＴ−ｎ番目の要素で構成される部分ベクトルとし、Σ_(c:m,d:n),w（c≧m, d≧n）を共分散行列Σ_w(Θ^ ⁽ⁱ⁾)の(T-c, T-d)番目の要素から(T-m, T-n)番目の要素（Ｔ−ｄ行目からＴ−ｎ行目かつＴ−ｃ列目からＴ−ｍ列目の各要素）で構成される部分行列とする。
２．ＣＭ−ｓｔｅｐ１の計算方法
ｔ番目のフレームにおける源信号の線形予測係数とその推定値が、それぞれ次のようなベクトルで表現される。

信号源パラメータ_sΘとその推定値_sΘ^は、それぞれ{a_t, _sσ_t ²}及び{a_t^, _sσ^_t ²}の全フレーム（0≦t≦T-1）にわたる集合と等価である。
式(25)による信号源パラメータの更新は、次式に示すa_t及び_sσ_t ²の推定値の更新を全フレーム（0≦t≦T-1）にわたって実行することで実現される。

ただし、_sR_t ⁽ⁱ⁾と_sr_t ⁽ⁱ⁾とV_t,w ⁽ⁱ⁾とは、それぞれ以下のように定義される。

３．ＣＭ−ｓｔｅｐ２の計算方法
ｗ番目の周波数帯域における残響パラメータとその推定値が、それぞれ次のようなベクトルで表現される。

残響パラメータ_gΘとその推定値_gΘ^は、それぞれg_w及びg_w^の全周波数帯域（0≦w≦N-1）にわたる集合と等価となる。
式(26)による残響パラメータの更新は、次式に示すg_wの推定値の更新を全周波数帯域（0≦w≦N-1）にわたって実行することで実現される。

ただし、_xR_w ⁽ⁱ⁾と_xr_w ⁽ⁱ⁾はそれぞれ以下のように定義される。

以上説明したように、本実施形態のパラメータ推定部では、雑音抑圧処理（Ｅ−ｓｔｅｐ）と信号源パラメータ推定値の更新処理（ＣＭ−ｓｔｅｐ１）と残響パラメータ推定値の更新処理（ＣＭ−ｓｔｅｐ２）とが協調的に繰り返して実行され、信号源パラメータ及び残響パラメータの推定値が更新される。Ｅ−ｓｔｅｐとＣＭ−ｓｔｅｐ１とは先に述べた第１更新処理に、ＣＭ−ｓｔｅｐ２は先に述べた第２更新処理に該当する。これにより、雑音と残響がともに存在する環境における観測信号から、雑音と残響とが精度よく抑圧され、源信号が強調される。

＜本実施形態の構成＞
次に、本実施形態の信号強調装置の構成を説明する。
図３は、第１実施形態の信号強調装置１の構成を示すブロック図である。また、図４は、源信号推定部２７の詳細構成を示すブロック図である。
図３に示すように、本実施形態の信号強調装置１は、観測信号記憶部１１、パラメータ記憶部１２、一時記憶部１３、帯域分割部２１、雑音パラメータ推定部２２、初期パラメータ設定部２３、雑音抑圧処理部２４、信号源パラメータ推定値更新部２５、残響パラメータ推定値更新部２６、源信号推定部２７、帯域合成部２８及び制御部２９を有する。また、源信号推定部２７は、残響重畳信号推定部２７ａ及び線形フィルタ適用部２７ｂを有する。なお、雑音パラメータ推定部２２及び初期パラメータ設定部２３は、前述の初期化部に対応する。また、雑音抑圧処理部２４及び信号源パラメータ推定値更新部２５は、前述の第１更新部に対応する。また、残響パラメータ推定値更新部２６は、前述の第２更新部に対応する。

なお、本実施形態の信号強調装置１は、ＣＰＵ（Central Processing Unit）、ＲＡＭ（Random Access Memory）等からなる公知のコンピュータに所定のプログラムが読み込まれることにより構成されるものである。具体的には、観測信号記憶部１１、パラメータ記憶部１２及び一時記憶部１３は、例えば、ＲＡＭ、レジスタ、キャッシュメモリ、若しくは補助記憶装置、又はそれらの少なくとも一部の結合によって構成される記憶部である。また、帯域分割部２１、雑音パラメータ推定部２２、初期パラメータ設定部２３、雑音抑圧処理部２４、信号源パラメータ推定値更新部２５、残響パラメータ推定値更新部２６、源信号推定部２７、帯域合成部２８及び制御部２９は、ＣＰＵに所定のプログラムが読み込まれることにより構成される本装置専用の処理部である。また、制御部２９は、信号強調装置１の各処理を制御する。

＜本実施形態の処理＞
図５は、第１実施形態の信号強調方法を説明するためのフローチャートである。以下、このフローチャートに沿って本実施形態の信号強調方法を説明する。
まず、信号強調装置１の帯域分割部２１に、雑音と残響とが共に存在する環境で観測され、所定の標本化周波数でサンプリングされ量子化された時間領域の観測信号Y_κが入力される。なお、κは離散時刻のインデックスを示す。帯域分割部２１は、短時間フーリエ変換等によって各離散信号Y_κを周波数帯域ごとの狭帯域信号に分割し、周波数領域の観測信号Y_t,wを生成し、観測信号記憶部１１に格納する（ステップＳ１）。なお、式(11)で示したように、Y={Y_t,w}_{0≦t≦T-1, 0≦w≦N-1}を観測信号の複素スペクトログラムと呼ぶ。

次に、雑音パラメータ推定部２２が、観測信号記憶部１１に格納された観測信号Y_t,wのうち、源信号が存在しない区間のものを用い、雑音パラメータの真値_dΘ^〜を推定する。なお、前述のように、本実施形態の雑音パラメータ_dΘは、雑音のパワースペクトル（雑音の確率分布を示す複素正規分布の分散）である。また、本実施形態の仮定では、雑音が定常であり、その振幅の平均が０である。そのため、雑音パラメータの真値_dΘ^〜は、源信号が存在しない区間の観測信号Y_t,wの振幅の２乗平均によって推定することができる。また、源信号が存在しない区間の特定には、例えば、公知の音声区間検出技術を用いる。あるいは、雑音パラメータ推定用に源信号が存在しない観測信号Y_t,wを予め計測しておき、それを用いてもよい。推定された雑音パラメータの最終的な推定値_dΘ^〜は、パラメータ記憶部１２に格納される（ステップＳ２）。

次に、初期パラメータ設定部２３が、信号源パラメータ及び残響パラメータの推定値の初期値_sΘ^⁽⁰⁾,_gΘ^⁽⁰⁾を設定する。例えば、初期パラメータ設定部２３は、観測信号記憶部１１から観測信号Y_t,wを読み込み、それを線形予測して得られた線形予測係数と予測残差パワーとを信号源パラメータの推定値の初期値_sΘ^⁽⁰⁾とし、_gΘ^⁽⁰⁾={{g_k.w^⁽⁰⁾=0}_1≦k≦Kw}_0≦w≦N-1を残響パラメータの推定値の初期値_gΘ^⁽⁰⁾とする。設定された各パラメータの推定値の初期値_sΘ^⁽⁰⁾,_gΘ^⁽⁰⁾は、パラメータ記憶部１２に格納される（ステップＳ３）。

次に、制御部２９が、繰り返し回数を示すインデクスiを0に設定し、一時記憶部１３に格納する（ステップＳ４）。

次に、雑音抑圧処理部２４に、観測信号記憶部１１から読み込まれた観測信号Y_t,wと、信号源パラメータの推定値_sΘ^⁽ⁱ⁾と、パラメータ記憶部１２から読み込まれた雑音パラメータの最終的な推定値_dΘ^〜と、残響パラメータの推定値_gΘ^⁽ⁱ⁾とが入力される。雑音抑圧処理部２４は、これらを用い、観測信号Y_t,wの集合Yとパラメータの推定値Θ^との組合せが与えられた場合における残響重畳信号X_t,wの集合Xの条件付事後分布p(X｜Y,Θ^）を特定する複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,Y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)を算出する（ステップＳ５）。具体的には、前述の式(29)〜(34)を用いて複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,Y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)が算出される。算出された複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,Y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)は、それぞれパラメータ記憶部１２に格納される。

次に、信号源パラメータ推定値更新部２５に、パラメータ記憶部１２から読み込まれた残響パラメータ推定値_gΘ^⁽ⁱ⁾と、複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,Y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)とが入力される。信号源パラメータ推定値更新部２５は、これらを用い、残響パラメータ_gΘを_gΘ^⁽ⁱ⁾として固定した状態で、式(24)に示した補助関数Q(Θ|Θ^⁽ⁱ⁾)の関数値が最大になるように信号源パラメータの推定値_sΘ^⁽ⁱ⁾を更新し、更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾を求める（ステップＳ６）。具体的には、式(36)〜(42)を用い、更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾を算出する。更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾はパラメータ記憶部１２に格納される。

次に、残響パラメータ推定値更新部２６に、パラメータ記憶部１２から読み込まれた信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾と、複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,Y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)とが入力される。残響パラメータ推定値更新部２６は、これらを用い、信号源パラメータ_sΘを_sΘ^⁽ⁱ⁺¹⁾として固定した状態で、式(24)に示した補助関数Q(Θ|Θ^⁽ⁱ⁾)の関数値が最大になるように残響パラメータの更新された推定値_gΘ^⁽ⁱ⁺¹⁾を求める（ステップＳ７）。具体的には、式(44)〜(46)を用い、更新された残響パラメータの推定値_gΘ^⁽ⁱ⁺¹⁾を算出する。更新された残響パラメータの推定値_gΘ^⁽ⁱ⁺¹⁾はパラメータ記憶部１２に格納される。

次に、所定の終了条件を充足するか否かを制御部２９（「終了条件判定部」に対応）が判定する（ステップＳ８）。ここで、所定の終了条件とは、例えば、各パラメータの推定値の更新量〔更新前のパラメータの推定値と更新後のパラメータの推定値との距離（コサイン距離やユークリッド距離等）〕がそれぞれ所定値以下となったことや、繰り返し回数を示すインデックスｉの値が所定値以上になったこと等を例示できる。

ここで、所定の終了条件を充足していなかった場合には、制御部２９は、繰り返し回数を示すインデックスｉの値を１だけ増やし、新たなインデックスｉの値を一時記憶部１３に格納する（ステップＳ９）。そして、ステップＳ１０５に戻る。

一方、所定の終了条件を充足していた場合には、制御部２９は、その時点における信号源パラメータ及び残響パラメータの推定値_sΘ^ ⁽ⁱ⁺¹⁾,_gΘ^⁽ⁱ⁺¹⁾を、信号源パラメータ最終推定値_sΘ^と雑音パラメータ最終推定値_gΘ^とし、それをパラメータ記憶部１２に格納する（ステップＳ１０）。

次に、源信号推定部２７に、観測信号Y_t,wと各パラメータの最終的な推定値_sΘ^,_gΘ^,_dΘ^〜とが入力される。源信号推定部２７は、これらを用い、源信号の推定値S_t,w^を生成する（ステップＳ１１）。そして、S^={S_t,w^}_{0≦t≦T-1, 0≦w≦N-1}が、源信号が強調された信号の複素スペクトログラムとなる。

具体的には、まず、源信号推定部２７の残響重畳信号推定部２７ａ（図４）に、観測信号Y_t,wと各パラメータの最終的な推定値_sΘ^,_gΘ^,_dΘ^〜とが入力される。残響重畳信号推定部２７ａは、これらを用い、観測信号Y_t,wと当該パラメータ推定値Θ^との組合せが与えられた場合における残響重畳信号X_t,wの条件付事後分布p(X｜Y,Θ^）の平均μ_w(Θ^,Y)（0≦w≦N-1）を残響重畳信号の推定値（「残響重畳信号最終推定値」に相当）として算出する。具体的には、前述の式(29)〜(34)でΘ^⁽ⁱ⁾をΘ^に置き換えることで平均μ_w(Θ^,Y)を算出する。算出された残響重畳信号の推定値μ_w(Θ^,Y)は、線形フィルタ適用部２７ｂに送られる。線形フィルタ適用部２７ｂには、算出された残響重畳信号の推定値μ_w(Θ^,Y)と、残響パラメータの最終的な推定値_gΘ^とが入力される。線形フィルタ適用部２７ｂは、入力された残響パラメータの推定値_gΘ^を用いて構成される線形フィルタを残響重畳信号の推定値μ_w(Θ^,Y)に適用し、源信号の推定値S_t,w^（「源信号最終推定値」に相当）を生成する。具体的には、線形フィルタ適用部２７ｂは、以下に従って、源信号の推定値S_t,w^を算出する。ただし、μ_t,wは、残響重畳信号の推定値μ_w(Θ^,Y)のT-t番目の要素である。

算出された源信号の推定値S_t,w^はパラメータ記憶部１２に格納される。
その後、帯域合成部２８に源信号の推定値S_t,w^が入力され、帯域合成部２８は、これを、逆短時間フーリエ変換などによって、時間領域の源信号の推定値S_κ^に変換して出力する（ステップＳ１２）。

＜実験結果＞
次に、本実施形態の処理を行って得られる効果を確認する実験を行った。まず、ASJ-JNASデー夕ベースから１０名（男性５名、女性５名）による発話を抽出した。発話の継続時間はすべて３秒間である。標本化周波数は８ｋＨｚ、量子化ビット数は１６ビットとした。これら源信号に残響時間がおよそ０．５秒の部屋で収録したインパルス応答を畳み込むことで残響重畳信号を合成した。これに、ＳＮＲ（Signal to Noise Ratio）が１０ｄＢとなるように計算機上で合成した定常白色雑音を加算して雑音残響重畳信号とした。

本実施形態の信号強調装置で用いるパラメータは下記の通り設定した。短時間フーリエ変換フレーム長は２５６サンプル、シフト幅は１２８サンプル、窓関数はハニング窓、室内伝達系を表す自己回帰の次数はすべての周波数帯域についてＫ_ｗ＝３０、源信号の線形予測次数はＰ＝１２とした。また、ＥＣＭアルゴリズムの終了条件は、繰り返し回数がｉ＝５回となったこととした。
強調後の源信号の品質は、次式で定義されるSASNR（Segmental Amplitude Signal to Noise Ratio）を用いて評価した。

表１に、話者の性別ごとのSASNRの改善値をまとめる。

表１に示すように、本実施形態の処理により、SASNRを平均で７．７２ｄＢ改善することができた。雑音抑圧処理のみでは、SASNRの平均改善値は４．２６ｄＢに低下した。一方、残響抑圧処理のみでは、SASNRの平均改善値は１．４９ｄＢに低下した。本実験結果から、本実施形態の方法を用いて雑音抑圧処理と残響抑圧処理を協調して動作させることによって、効果的な源信号強調を実現できたことが確認された。

〔第２実施形態〕
次に、本発明の第２実施形態を説明する。第１実施形態では、信号を測定するセンサが１個に限定されていたのに対して、本実施形態では、信号を観測するセンサの個数に制限が設けられない。すなわち、センサの個数ＭはＭ≧１を満たす任意の整数をとる。よって、残響パラメータに含まれる回帰行列は、Ｍ行Ｍ列の正方行列である。それ以外の点については、本実施形態におけるパラメータ推定処理の概要は、第１実施形態におけるパラメータ推定処理の概要と同じである。また、Ｍ＝１であってもよいし、Ｍ≧２であってもよく、Ｍ＝１とした本実施形態は、第１実施形態と等価になる。

＜本形態のパラメータ推定処理の概要＞
本実施形態では、第１更新部は第２パラメータ群のパラメータの推定値を更新し、第２更新部は第１パラメータ群のパラメータの推定値を更新する。
[観測信号記憶処理]
まず、観測信号記憶処理によって、観測信号が記憶部に格納される。
[初期化処理]
次に、初期化処理によって、第１パラメータ群のパラメータの推定値と、第２パラメータ群のパラメータの推定値とが初期化される。

[第１更新処理]
本実施形態の第1更新処理では、第1パラメータ群、すなわち残響パラメータの推定値が固定された状態で、第２パラメータ群、すなわち信号源パラメータの推定値が更新される。本実施形態の第１更新処理は、具体的には、雑音抑圧処理と、信号源パラメータの更新処理とを含む。

《雑音抑圧処理》
雑音抑圧処理では、観測信号とパラメータの推定値を用いて、残響重畳信号の条件付事後分布ｐ（残響重畳信号｜観測信号，パラメータの推定値）を特徴づける複素正規分布の平均と共分散行列が算出される。

この処理は、観測信号から雑音を含まない残響重畳信号の条件付事後分布を求めるという点において、観測信号に含まれる雑音を抑圧していると解釈できる。この雑音抑圧処理は、残響パラメータの推定値と信号源パラメータの推定値を用いて実行されることに注意されたい。このことは、残響の特性が考慮されながら雑音が抑圧されることを意味する。これによって、残響環境において、雑音抑圧を精度よく実施できる。

《信号源パラメータ推定値の更新処理》
信号源パラメータ推定値の更新処理では、残響パラメータの推定値と残響重畳信号の条件付事後分布の平均と共分散行列を用いて、信号源パラメータの推定値が更新される。信号源パラメータの推定値は、パラメータの推定値に関する補助関数の値が最大になるように、更新される。

補助関数は、観測信号と残響重畳信号を所与とした場合のパラメータの推定値に関する対数尤度関数を、残響重畳信号の条件付事後分布ｐ（残響重畳信号｜観測信号，パラメータ推定値）で重み付けした関数を、残響重畳信号について積分して得られる関数である。この重み付け積分により、雑音抑圧処理で算出される残響重畳信号の不確かさを考慮しながら、信号源パラメータの推定値を更新することが可能になっている。

[第２更新処理]
本実施形態の第２更新処理では、第２パラメータ群、すなわち信号源パラメータの推定値が固定された状態で、第１パラメータ群、すなわち残響パラメータの推定値が更新される。残響パラメータの推定値は、パラメータの推定値に関する補助関数の値が最大になるように、更新される。

[終了条件判定処理]
終了条件判定処理では、所定の終了条件が満たされているか否かが判定される。終了条件が満たされていない場合、第１更新処理に戻る。終了条件が満たされている場合、その時点におけるパラメータの推定値を出力する。

以上で述べた処理において、残響重畳信号の条件付事後分布の共分散行列のスケールは、雑音の共分散行列のスケールに対して単調増加する。すなわち、雑音のレベルが大きいほど、残響重畳信号の条件付事後分布の共分散行列のスケールも大きくなる。このことは、本実施形態が、雑音抑圧処理で求められる残響重畳信号の不確かさを妥当な方法で評価していることを示している。

＜本実施形態の原理＞
次に、本実施形態の原理を説明する。以下では、第１実施形態との相違点を中心に説明し、第１実施形態と共通する事項については説明を省略する。なお、本実施形態でも、信号は音声信号などの音響信号に限定されない。

＜本実施形態の原理＞
次に、本実施形態の原理を説明する。本実施形態でもＥＣＭアルゴリズムを適用する。すなわち、観測信号である雑音残響重畳信号の集合yを用い、雑音残響重畳信号の集合yとパラメータの推定値Θ^との組合せを前提条件とした残響重畳信号の集合xの条件付事後分布p(x|y,Θ^)の算出処理（Ｅ−ｓｔｅｐ）と、源信号パラメータの推定値_sΘ^の算出処理（ＣＭ−ｓｔｅｐ１）と、残響パラメータ_gΘの算出処理（ＣＭ−ｓｔｅｐ２）とを代わる代わる繰り返し実行して各推定値を更新し、所定の終了条件を充足した時点での各推定値を真値の推定値（最終推定値）とする。なお、Ｅ−ｓｔｅｐとＣＭ−ｓｔｅｐ１は先に述べた第１更新処理に、ＣＭ−ｓｔｅｐ２は先に述べた第２更新処理に該当する。

なお、本実施形態の残響重畳信号の集合xは、各センサにそれぞれ対応する残響重畳信号の複素スペクトログラムを要素とした集合である。また、本実施形態の雑音残響重畳信号の集合yは、各センサにそれぞれ対応する雑音残響重畳信号の複素スペクトログラムを要素とした集合である。

［観測信号（雑音残響重畳信号）の統計的モデル］
本実施形態でも、まず、パラメータΘが与えられた場合における雑音残響重畳信号集合のyの確率密度関数p(y|Θ)が定義される。そのために、観測信号（雑音残響重畳信号）の集合yの統計的モデルが仮定される。本実施形態では、以下に述べる源信号の全極モデル、室内伝達系の多チャンネル自己回帰モデル及び雑音のモデルが仮定される。

《源信号のモデル》
まず、本実施形態の源信号の全極モデルについて述べる。t(0≦t≦T-1)番目のフレーム、w(0≦w≦N-1)番目の周波数帯域における源信号の離散フーリエ係数（複素数）をS_t,wとおく。また、仮に雑音や残響が存在しない場合に、m(1≦m≦M)番目のセンサで観測されるであろう源信号の離散フーリエ係数をS_t,w ^(m)とおく。また、各S_t,w ^(m)を要素とする次のようなＭ次元の源信号ベクトルが定義される。なお、α^τはαの非共役転置を示す。
s_t,w=[S_t,w ⁽¹⁾,...,S_t,w ^(M)]^τ (49)

ベクトルs_t,wが以下の条件を満たすと仮定される。
１．ω∈{‐π,π}を角周波数として、ｔ番目のフレームにおける源信号のパワースペクトル密度_sλ_t(ω)は、式(1)(2)に示したような全極型スペクトル密度で表される。よって、信号源パラメータ_sΘは、_sΘ={a_t,1 ,..., a_t,P, _sσ_t ²}_0≦t≦T-1と定義される。ただし、{m_α}_0≦α≦M-1は、m₀, m₁ ,..., m_M-1のＭ個の要素からなる集合を表す。
２．s_t,wは、以下のような、平均０_M、共分散行列_sλ_t(2πw/N)I_MのＭ次元複素正規分布にしたがう。

ただし、N_C{x;μ,Σ}は、式(4)で定義される複素正規分布の確率密度関数である。また、０_MとI_Mは、それぞれ、Ｍ次元零ベクトルとＭ次元単位行列を表す。
ここで、ζ＝Ｍとして式(4)を式(50)に代入するとs_t,wの確率密度関数は次式で表される。

ただし、複素ベクトルαに対する||α||²は、次式により定義される。
||α||²=α^H・α (52)
３．(t,w)≠(t',w')ならば、s_t,wとs_t',w'は統計的に独立である。
《室内伝達系のモデル》
次に、本実施形態の室内伝達系のモデルについで述べる。m(1≦m≦M)番目のセンサ、t(0≦t≦T-1)番目のフレーム、w(0≦w≦N-1)番目の周波数帯域における残響重畳信号の離散フーリエ係数をX_t,w ^(m)とおく。また、各X_t,w ^(m)を要素とする次のようなＭ次元の残響重畳信号ベクトルが定義される。
x_t,w=[X_t,w ⁽¹⁾,...,X_t,w ^(M)]^τ (53)
本実施形態では、室内伝達系が各周波数帯域においてＭチャネル自己回帰系として表現できると仮定される。すなわち、ｗ番目の周波数帯域における回帰系の回帰行列を

とおくと、残響重畳信号の残響重畳信号ベクトルx_t,wは次式により生成される。

なお、回帰行列Ｇ_ｋ，ｗは、回帰系の回帰係数g_k,w ^(1,1),..., g_k,w ^(M,M)を要素に持つ以下のようなＭ行Ｍ列の行列である。なお、K_wはＭチャネル自己回帰系の次数を示す。

式(55)を用いると式(54)は以下のように表現される。

本実施形態では、_gΘ={{G_k.w}_1≦k≦Kw}_0≦w≦N-1が残響パラメータ_gΘと定義される。この残響パラメータ_gΘは、次式に示すように、源信号に残響のみが付加された残響重畳信号に適用されて、各センサ位置での源信号を抽出する用途に供される。

《雑音のモデル》
次に、雑音のモデルについて述べる。本実施形態では、m(1≦m≦M)番目のセンサ、t(0≦t≦T-1)番目のフレーム、w(0≦w≦N-1)番目の周波数帯域における、雑音と雑音残響重畳信号との離散フーリエ係数がそれぞれD_t,w ^(m)，Y_t,w ^(m)とされる。また、各D_t,w ^(m)を要素とする次のようなＭ次元の雑音ベクトルが定義される。
d_t,w=[D_t,w ⁽¹⁾,...,D_t,w ^(M)]^τ (58)

同様に、各Y_t,w ^(m)を要素とする次のようなＭ次元の雑音残響重畳信号（観測信号）ベクトルが定義される。
y_t,w=[Y_t,w ⁽¹⁾,...,Y_t,w ^(M)]^τ (59)
雑音残響重畳信号ベクトルy_t,wは、残響重畳信号ベクトルx_t,wに雑音ベクトルd_t,wを加算したものである。
y_t,w = x_t,w + d_t,w (60)

また、d_t,wが次に述べる条件を満たすと仮定される。
１．雑音は定常であり、そのパワークロススペクトル密度を_dΛ(ω)として（定常であるためフレーム番号ｔには依存しない）、d_t,wは平均０_M、共分散行列_dΛ(2πw/N)の複素正規分布に従う。なお、共分散行列_dΛ(2πw/N)のｗ番目の対角要素は、ｗ番目のセンサにおける雑音のパワースペクトル_dλ^(m)(2πw/N)である。

また、本実施形態の雑音パラメータ_dΘは、_dΘ={_dΛ(2πw/N)}_{0≦ｗ≦N-1}と定義される雑音を特徴づけるパラメータである。
２．(t, w)≠(t', w')ならば、d_t,wとd_t',w'とは統計的に独立である。
３．任意の(t, w, t', w')について、s_t,wとd_t,wとは統計的に独立である。

《雑音残響重畳信号の確率密度関数》
以上の仮定に基づき、雑音残響重畳信号の確率密度関数が定式化される。
本実施形態では、各センサにおける源信号の複素スペクトログラムからなる集合（源信号ベクトルの集合に相当）がsと表現される。また、各センサにおける残響重畳信号の複素スペクトログラムからなる集合（残響重畳信号ベクトルの集合に相当）がxと表現される。また、雑音残響重畳信号の複素スペクトログラムからなる集合（雑音残響重畳信号ベクトルの集合に相当）がyと表現される。
すなわち、
s={s_t,w}_{0≦t≦T-1, 0≦w≦N-1} (62)
x={x_t,w}_{0≦t≦T-1, 0≦w≦N-1} (63)
y={y_t,w}_{0≦t≦T-1, 0≦w≦N-1} (64)
と表現される。

具体的には、雑音残響重畳信号ベクトルの集合yの確率密度関数（観測信号ベクトルの集合yが与えられたときのパラメータΘに関する尤度関数に相当）は次のように書ける。

ただし、p(y,x|Θ)は、以上の仮定に基づいて次式のように書ける。

以上で、パラメータΘ={_sΘ, _gΘ, _dΘ} を用いて雑音残響重畳信号の集合の確率密度関数p(y|Θ)が定式化された。

［信号源パラメータ及び残響パラメータの最尤推定］
前述のように、本実施形態では、観測された雑音残響重畳信号の集合のyから、未知のパラメータの真値Θ^〜が、最尤推定法によって推定される。すなわち、雑音残響重畳信号の集合Ｙが与えられた場合におけるパラメータΘを変数とした尤度関数p(Y|Θ)を最大化するΘが、真値Θ^〜の推定値となる。ただし、本実施形態では、雑音パラメータの真値_dΘ^〜が源信号の存在しない区間から予め独立に推定され、既知となっている為Θ^={_sΘ^, _gΘ^, _dΘ^〜}であり、_sΘ^と_gΘ^が求められることになる。

また、尤度関数p(Y|Θ)を最大化する_sΘ^と_gΘ^を同時に直接求めることはできないから、ＥＣＭアルゴリズムを用いてこれらが計算される。ＥＣＭアルゴリズムの処理の流れを以下に示す。以下の処理では、Ｅ−ｓｔｅｐ、ＣＭ−ｓｔｅｐ１、ＣＭ−ｓｔｅｐ２の３つの処理が代わる代わる繰り返し実行される。そこで、ｉ回目の繰り返しにおけるパラメータの推定値を上付きの添え字(i)を用いて示す。明確さを期するために述べると、Θ^〜，Θ^，Θ^⁽ⁱ⁾はそれぞれ次のように定義される。

《ＥＣＭアルゴリズム》
１．パラメータの推定値の初期値Θ^⁽⁰⁾が決められる。また、繰り返し回数を示すインデックスｉが０にされる。
２．Ｅ−ｓｔｅｐ（雑音抑圧処理）
残響重畳信号の条件付事後分布p(x|y, Θ^⁽ⁱ⁾)が計算される。
３．ＣＭ−ｓｔｅｐ１（信号源パラメータ推定値の更新処理）
補助関数Q(Θ|Θ^⁽ⁱ⁾)が次式により定義される。

このとき、次の手続きにより、信号源パラメータの推定値が_SΘ^⁽ⁱ⁾から_SΘ^⁽ⁱ⁺¹⁾に更新される。

すなわち、残響パラメータの推定値_gΘ^⁽ⁱ⁾が固定された条件下で補助関数Q(Θ|Θ^⁽ⁱ⁾)を最大化する_SΘ^⁽ⁱ⁺¹⁾が、更新された信号源パラメータの推定値とされる。
４．ＣＭ−ｓｔｅｐ２（残響パラメータ推定値の更新処理）
次の手続きにより、残響パラメータの推定値が更新される。

すなわち、信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾が固定された条件下で補助関数Q(Θ|Θ^⁽ⁱ⁾)を最大化する_gΘ^⁽ⁱ⁺¹⁾が、更新された残響パラメータの推定値とされる。
５．終了条件判定
所定の終了条件を満たしているならば_sΘ^=_sΘ^⁽ⁱ⁺¹⁾，_gΘ^=_gΘ^⁽ⁱ⁺¹⁾として終了。そうでなければ、ｉを１だけ漸増させて「２．Ｅ−ｓｔｅｐ」へ戻る。

《各ｓｔｅｐの計算方法》
以下では、Ｅ−ｓｔｅｐ、ＣＭ−ｓｔｅｐ１及びＣＭ−ｓｔｅｐ２の各計算方法を説明する。
１．Ｅ−ｓｔｅｐの計算方法
すべてのセンサにおける、源信号、残響重畳信号、雑音残響重畳信号のｗ番目の周波数帯域の離散フーリエ係数系列を、それぞれまとめて次のように表す。

源信号ベクトルの集合s、残響重畳信号ベクトルの集合x及び雑音残響重畳信号ベクトルの集合yは、それぞれ、s_w, x_w, y_wの全周波数帯域（0≦w≦N-1）にわたる集合と等価となる。
式(77)の残響重畳信号の条件付事後分布p(x|y, Θ^⁽ⁱ⁾)は、次式に示すように周波数帯域wごとに独立な複数の複素正規分布によって表現できる。

なお、平均μ_w(Θ^⁽ⁱ⁾,y)と共分散行列Σ_w(Θ^⁽ⁱ⁾)は次式で与えられる。また、平均μ_w(Θ^⁽ⁱ⁾,y)はＭ次元ベクトルである。

式(82),(83)に現れる各変数はそれぞれ以下のように定義される。なお、式(84)の空欄部分の各要素は０である。

なお、bdiag{Ω_１,...,Ω_α}は、任意の正方行列Ω_１,...,Ω_αに対する次のブロック対角行列を示す。

また、前述のように、雑音が定常であると仮定されているため、
_dΛ_T-1 ^〜(2πw/N)= _dΛ_T-2 ^〜(2πw/N)=...=_dΛ₀ ^〜(2πw/N)=_dΛ^〜(2πw/N) (89)
である。
また、後の処理のために、μv_m,w ⁽ⁱ⁾を平均μ_w(Θ^⁽ⁱ⁾,y)のM(T-m-1)+1からM(T-m)番目までの要素で構成される部分ベクトルとし、μv_m:n,w ⁽ⁱ⁾（m≧n）を平均μ_w(Θ^⁽ⁱ⁾,y)のM(T-m-1)+1からM(T-m)番目までの要素で構成される部分ベクトルとする。また、ΣV_{(m1:n1,m2:n2),w} ⁽ⁱ⁾を共分散行列Σ_w(Θ^ ⁽ⁱ⁾)の(M(T-m1-1)+1,M(T-m2-1)+1)番目の要素から(M(T-n1),M(T-n2))番目の要素で構成される部分行列とする。

２．ＣＭ−ｓｔｅｐ１の計算方法
ｔ番目のフレームにおける源信号の線形予測係数とその推定値が、式(35)のようなベクトルで表現される。
信号源パラメータ_sΘとその推定値_sΘ^は、それぞれ{a_t, _sσ_t ²}及び{a_t^, _sσ^_t ²}の全フレーム（0≦t≦T-1）にわたる集合と等価である。
式(78)による信号源パラメータの更新は、式(36)(37)に示したa_t及び_sσ_t ²の推定値の更新を全フレーム（0≦t≦T-1）にわたって実行することで実現される。ただし、本実施形態では、式(41)(42)に代えて

で算出されるV_t,w ⁽ⁱ⁾を用い、式(36)から(40)の計算によって、a_t及び_sσ_t ²の推定値が更新される。なお、正方行列Αに対する式(90)のdavg(Α)は、正方行列Αの対角要素の平均値を表す。
３．ＣＭ−ｓｔｅｐ２の計算方法
ｗ番目の周波数帯域における残響パラメータとその推定値が、それぞれ次のようなベクトルで表現される。

残響パラメータ_gΘとその推定値_gΘ^は、それぞれG_w及びG_w^の全周波数帯域（0≦w≦N-1）にわたる集合と等価となる。
式(78)による残響パラメータの更新は、次式に示すG_wの推定値の更新を全周波数帯域（0≦w≦N-1）にわたって実行することで実現される。

ただし、_xRV_w ⁽ⁱ⁾と_xrv_w ⁽ⁱ⁾はそれぞれ以下のように定義される。

以上説明したように、本実施形態では、雑音抑圧処理（Ｅ−ｓｔｅｐ）と信号源パラメータ推定値の更新処理（ＣＭ−ｓｔｅｐ１）と残響パラメータ推定値の更新処理（ＣＭ−ｓｔｅｐ２）とが協調的に繰り返して実行され、信号源パラメータ及び残響パラメータの推定値が更新される。これにより、雑音と残響がともに存在する環境における観測信号から、雑音と残響とが精度よく抑圧され、源信号が強調される。

＜本実施形態の構成＞
次に、本実施形態の信号強調装置の構成を説明する。
図６は、第２実施形態の信号強調装置１００の構成を示すブロック図である。また、図７は、源信号推定部１２７の詳細構成を示すブロック図である。

図６に示すように、本実施形態の信号強調装置１００は、観測信号記憶部１１１、パラメータ記憶部１１２、一時記憶部１３、帯域分割部１２１、雑音パラメータ推定部１２２、初期パラメータ設定部１２３、雑音抑圧処理部１２４、信号源パラメータ推定値更新部１２５、残響パラメータ推定値更新部１２６、源信号推定部１２７、帯域合成部２８及び制御部２９を有する。また、源信号推定部１２７は、残響重畳信号推定部１２７ａ及び線形フィルタ適用部１２７ｂを有する。なお、雑音パラメータ推定部１２２及び初期パラメータ設定部１２３は、前述の初期化部に対応する。また、雑音抑圧処理部１２４及び信号源パラメータ推定値更新部１２５は、前述の第１更新部に対応する。また、残響パラメータ推定値更新部１２６は、前述の第２更新部に対応する。

なお、本実施形態の信号強調装置１００は、ＣＰＵ、ＲＡＭ等からなる公知のコンピュータに所定のプログラムが読み込まれることにより構成されるものである。具体的には、観測信号記憶部１１１、パラメータ記憶部１１２及び一時記憶部１３は、例えば、ＲＡＭ、レジスタ、キャッシュメモリ、若しくは補助記憶装置、又はそれらの少なくとも一部の結合によって構成される記憶部である。また、帯域分割部１２１、雑音パラメータ推定部１２２、初期パラメータ設定部１２３、雑音抑圧処理部１２４、信号源パラメータ推定値更新部１２５、残響パラメータ推定値更新部１２６、源信号推定部１２７、帯域合成部２８及び制御部２９は、ＣＰＵに所定のプログラムが読み込まれることにより構成される本装置専用の処理部である。また、制御部２９は、信号強調装置１００の各処理を制御する。

＜本実施形態の処理＞
図８は、第２実施形態の信号強調方法を説明するためのフローチャートである。以下、このフローチャートに沿って本実施形態の信号強調方法を説明する。
まず、信号強調装置１００の帯域分割部１２１に、Ｍ個のセンサによってそれぞれ観測され、量子化された時間領域の観測信号Y_κ ^(ｍ)(1≦m≦M)を要素とする観測信号ベクトル[Y_κ ⁽¹⁾,...,Y_κ ^(M)]^τが入力される。帯域分割部１２１は、短時間フーリエ変換等によって観測信号ベクトル[Y_κ ⁽¹⁾,...,Y_κ ^(M)]^τを、時間周波数領域の観測信号ベクトルy_t,w= [Y_t,w ⁽¹⁾,...,Y_t,w ^(M)]^τに変換し、観測信号記憶部１１１に格納する（ステップＳ１０１）。

次に、雑音パラメータ推定部１２２が、観測信号記憶部１１１に格納された観測信号ベクトルy_t,wのうち、源信号が存在しない区間のものを用い、雑音パラメータの真値_dΘ^〜の推定値を計算する。なお、前述のように、本実施形態の雑音パラメータ_dΘは、雑音のパワークロススペクトル（雑音の確率分布を示すＭ次元複素正規分布の共分散行列）である。また、本実施形態では、雑音は定常であり、その振幅の平均は０_Ｍであると仮定している。そのため、雑音パラメータの真値_dΘ^〜は、源信号が存在しない区間の観測信号ベクトルy_t,wを用いて、次式のように推定することができる。

ただし、ηは源信号が存在しない区間のフレーム番号の集合であり、|η|は源信号が存在しない区間のフレーム数である。また、源信号が存在しない区間の特定には、例えば、公知の音声区間検出技術を用いる。あるいは、雑音パラメータ推定用に源信号が存在しない観測信号Y_t,wを予め計測しておき、それを用いてもよい。推定された雑音パラメータの真値_dΘ^〜は、パラメータ記憶部１１２に格納される（ステップＳ１０２）。

次に、初期パラメータ設定部１２３が、信号源パラメータ及び残響パラメータの推定値の初期値_sΘ^⁽⁰⁾,_gΘ^⁽⁰⁾を設定する。例えば、初期パラメータ設定部１２３は、観測信号記憶部１１１から観測信号ベクトルy_t,wを読み込み、その第１要素（すなわち、一番目のセンサで観測された信号）を線形予測分析して得られた線形予測係数と予測残差パワーとを信号源パラメータの推定値の初期値_sΘ^⁽⁰⁾とし、_gΘ^⁽⁰⁾={{G_k.w^⁽⁰⁾=O _M}_1≦k≦Kw}_0≦w≦N-1を残響パラメータの推定値の初期値_gΘ^⁽⁰⁾とする。ただし、O_MはＭ次元零行列である。設定された各パラメータの推定値の初期値_sΘ^⁽⁰⁾,_gΘ^⁽⁰⁾は、パラメータ記憶部１１２に格納される（ステップＳ１０３）。

次に、制御部２９が、繰り返し回数を示すインデクスiを0に設定し、一時記憶部１３に格納する（ステップＳ１０４）。

次に、雑音抑圧処理部１２４に、観測信号記憶部１１１から読み込まれた観測信号ベクトルy_t,wと、信号源パラメータの推定値_sΘ^⁽ⁱ⁾と、パラメータ記憶部１１２から読み込まれた雑音パラメータの真値_dΘ^〜と、残響パラメータの推定値_gΘ^⁽ⁱ⁾とが入力される。雑音抑圧処理部１２４は、これらを用い、観測信号ベクトルy_t,wの集合yとパラメータの推定値Θ^との組合せが与えられた場合における残響重畳信号ベクトルx_t,wの集合xの条件付事後分布p(x｜y,Θ^）を特定する複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)を算出する（ステップＳ１０５）。具体的には、前述の式(82)〜(87)を用いて複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)を算出する。算出された複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)は、それぞれパラメータ記憶部１１２に格納される。

次に、信号源パラメータ推定値更新部１２５に、パラメータ記憶部１１２から読み込まれた残響パラメータ推定値_gΘ^⁽ⁱ⁾と、複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)とが入力される。信号源パラメータ推定値更新部１２５は、これらを用い、残響パラメータ_gΘを_gΘ^⁽ⁱ⁾として固定した状態で、式(77)に示した補助関数Q(Θ|Θ^⁽ⁱ⁾)の関数値が最大になるように信号源パラメータの推定値_sΘ^⁽ⁱ⁾を更新し、更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾を求める（ステップＳ１０６）。具体的には、式(36)〜(40),(90),(91)を用い、更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾を算出する。更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾はパラメータ記憶部１１２に格納される。

次に、残響パラメータ推定値更新部１２６に、パラメータ記憶部１１２から読み込まれた信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾と、複素正規分布の平均μ_w(Θ^⁽ⁱ⁾,y)と、共分散行列Σ_w(Θ^⁽ⁱ⁾)とが入力される。残響パラメータ推定値更新部１２６は、これらを用い、信号源パラメータ_sΘを_sΘ^⁽ⁱ⁺¹⁾として固定した状態で、式(77)に示した補助関数Q(Θ|Θ^⁽ⁱ⁾)の関数値が最大になるように残響パラメータの更新された推定値_gΘ^⁽ⁱ⁺¹⁾を求める（ステップＳ１０７）。具体的には、式(93)〜(95)を用い、残響パラメータの推定値_gΘ^⁽ⁱ⁺¹⁾を算出する。更新された残響パラメータの推定値_gΘ^⁽ⁱ⁺¹⁾はパラメータ記憶部１１２に格納される。

次に、所定の終了条件を充足するか否かを制御部２９（「終了判定部」に対応）が判定する（ステップＳ１０８）。ここで、所定の終了条件とは、例えば、各パラメータの推定値の更新量〔更新前のパラメータの推定値と更新後のパラメータの推定値との距離（コサイン距離やユークリッド距離等）〕がそれぞれ所定値以下となったことや、繰り返し回数を示すインデックスｉの値が所定値以上になったこと等を例示できる。

ここで、所定の終了条件を充足していなかった場合には、制御部２９は、繰り返し回数を示すインデックスｉの値を１だけ増やし、新たなインデックスｉの値を一時記憶部１３に格納する（ステップＳ１０９）。そして、ステップＳ１０５に戻る。

一方、所定の終了条件を充足していた場合には、制御部２９は、その時点における信号源パラメータ及び残響パラメータの推定値_sΘ^ ⁽ⁱ⁺¹⁾,_gΘ^⁽ⁱ⁺¹⁾を信号源パラメータ最終推定値_sΘ^と残響パラメータ最終推定値_gΘ^とし、それをパラメータ記憶部１１２に格納する（ステップＳ１１０）。

次に、源信号推定部１２７に、観測信号Y_t,wと各パラメータの最終的な推定値_sΘ^,_gΘ^,_dΘ^〜とが入力される。源信号推定部１２７は、これらを用い、源信号の推定値S_t,w^を生成する（ステップＳ１１１）。そして、S^={S_t,w^}_{0≦t≦T-1, 0≦w≦N-1}が、源信号が強調された信号の複素スペクトログラムとなる。

具体的には、まず、源信号推定部１２７の残響重畳信号推定部１２７ａ（図７）に、観測信号ベクトルy_t,wと各パラメータの最終的な推定値_sΘ^,_gΘ^,_dΘ^〜とが入力される。残響重畳信号推定部１２７ａは、これらを用い、観測信号ベクトルy_t,wと当該パラメータ推定値Θ^との組合せが与えられた場合における残響重畳信号ベクトルx_t,wの条件付事後分布p(x｜y,Θ^）の平均μ_w(Θ^,y)（0≦w≦N-1）を残響重畳信号ベクトルx_t,wの推定値（「残響重畳信号最終推定値」に相当）として算出する。具体的には、前述の式(82)〜(87)でΘ^⁽ⁱ⁾をΘ^に置き換えることで平均μ_w(Θ^,y)を算出する。算出された残響重畳信号ベクトルx_t,wの推定値μ_w(Θ^,y)は、線形フィルタ適用部１２７ｂに送られる。

線形フィルタ適用部１２７ｂには、算出された残響重畳信号ベクトルx_t,wの推定値μ_w(Θ^,y)と、残響パラメータの最終的な推定値_gΘ^とが入力される。線形フィルタ適用部１２７ｂは、入力された残響パラメータの推定値_gΘ^を用いて構成される線形フィルタを残響重畳信号ベクトルx_t,wの推定値μ_w(Θ^,y)に適用し、源信号ベクトルの推定値s_t,w^を生成する。そして、線形フィルタ適用部１２７ｂは、例えば、源信号ベクトルの推定値s_t,w^の要素を平均し、その平均値を源信号の推定値S_t,w^（「源信号最終推定値」に相当）として出力する。具体的には、線形フィルタ適用部１２７ｂは、例えば、以下に従って、源信号の推定値S_t,w^を算出する。ただし、μv_t,wは、残響重畳信号ベクトルx_t,wの推定値μ_w(Θ^,y)のM(T-t-1)+1からM(T-t)番目までの要素で構成される部分ベクトルである。

ただし、任意のベクトルαに対するavg(α)は、ベクトルαの全要素の平均値を表す。なお、本実施形態では、

の要素の平均値を源信号の推定値S_t,w^としたが、これらの要素の何れかを源信号の推定値S_t,w^としてもよい。
算出された源信号の推定値S_t,w^はパラメータ記憶部１１２に格納される。
その後、帯域合成部２８に源信号の推定値S_t,w^が入力され、帯域合成部２８は、これを、逆短時間フーリエ変換などによって、源信号の推定値S_κ^に変換して出力する（ステップＳ１１２）。

＜実験結果＞
次に、本実施形態の処理を行って得られる効果を確認する実験を行った。男女２話者により発話された音声を用意した。各音声の音響信号に対して、残響時間が約０．５秒の部屋で２個のマイクロホンで収録したインパルス応答を畳み込むことで、残響音声信号を合成した。これに、ＳＮ比が１５ｄＢとなる白色雑音を加算することで、雑音残響音声信号をシミュレートした。

本実施形態を実施するのに必要なパラメータは下記の通り設定した。短時間フーリエ変換のフレーム長は２５６サンプル、シフト幅は１２８サンプル、窓関数はハニング窓、室内伝達系の次数は２５、音声の線形予測次数は１２とした。また、ＥＣＭアルゴリズムの終了条件は，繰り返し回数が３回となった時点とした。強調後の音声信号の品質を評価する尺度として、ケプストラム歪みを用いた。

本実施形態による処理を行う前の信号（雑音残響音声信号）のケプストラム歪みの平均値は，６．９９ｄＢであった．これに対して，本実施形態による処理を行った後の信号のケプストラム歪みの平均値は５．１５ｄＢであり，１．８４ｄＢ改善された。参考までに、マイクロホンを１個だけ用いた場合、ケプストラム歪みの平均値は５．６１ｄＢであった。以上の結果により，本実施形態の効果が確認された。

〔第３実施形態〕
次に、第３実施形態を説明する。
＜本実施形態のパラメータ推定処理の概要＞
まず、本実施形態のパラメータ推定部における処理の概要を説明する。本実施形態では、第２パラメータ群は、信号源パラメータに加えて、少なくとも、ステアリングベクトルを含む。また、本実施形態では、第１更新部は第２パラメータ群の推定値を更新し、第２更新部は第１パラメータ群のパラメータの推定値を更新する。

[観測信号記憶処理]
まず、観測信号記憶処理によって、観測信号が記憶部に格納される。
[初期化処理]
次に、初期化処理によって、第１パラメータ群のパラメータの推定値と、第２パラメータ群のパラメータの推定値とが初期化される。
[第１更新処理]
本実施形態の第１更新処理では、第1パラメータ群、すなわち残響パラメータの推定値が固定された状態で、第２パラメータ群、すなわち信号源パラメータの推定値が更新される。本実施形態の第１更新処理は、具体的には、源信号推定値更新処理、ステアリングベクトル推定値更新処理、信号源パラメータ推定値更新処理を含む。

《源信号推定値更新処理》
源信号推定値更新処理では、まず、観測信号と残響パラメータの推定値を用いて、雑音重畳信号の推定値を算出する。この処理は、雑音残響重畳信号を入力として雑音重畳信号を出力するという点において、残響抑圧処理に相当すると解釈される。

次に、算出された雑音重畳信号の推定値とパラメータの推定値を用いて、源信号の条件付事後分布ｐ（源信号｜雑音重畳信号の推定値，パラメータの推定値）を特徴づける複素正規分布の平均と分散が算出される。この平均と分散は、それぞれ、源信号の推定値と誤差分散に相当する。

《ステアリングベクトル推定値更新処理》
ステアリングベクトル推定値更新処理では、雑音重畳信号推定値と源信号推定値とを用いて、ステアリングベクトルの推定値が更新される。ステアリングベクトルの推定値は、パラメータに関する対数尤度関数が増加するように、更新される。

《信号源パラメータ推定値更新処理》
信号源パラメータ推定値更新処理では、源信号の推定値と誤差分散から、源信号のパワースペクトルの推定値を算出する。このパワースペクトルの推定値に基づいて、信号源パラメータの推定値が更新される。この更新処理は、パラメータに関する対数尤度関数を増加させる。

[第２更新処理]
本実施形態の第２更新処理では、第２パラメータ群、すなわち信号源パラメータ、雑音パラメータ、ステアリングベクトルの各々の推定値が固定された状態で、第１パラメータの群、すなわち残響パラメータの推定値が更新される。本実施形態の第２更新処理は、具体的には、源信号短時間パワースペクトル推定値更新処理、残響パラメータ推定値更新処理、雑音パラメータ推定値更新処理を含む。

《源信号短時間パワースペクトル推定値更新処理》
源信号短時間パワースペクトル推定値更新処では、信号源パラメータ推定値を用いて源信号のパワースペクトルの推定値を更新する。

《雑音パラメータ推定値更新処理》
次に、雑音パラメータ推定値更新処理では、雑音重畳信号の推定値、源信号の推定値、ステアリングベクトルの推定値を用いて、雑音パラメータの推定値を更新する。この更新処理は、パラメータに関する対数尤度関数を増加させる。

《残響パラメータ推定値更新処理》
残響パラメータ推定値更新処理では、観測信号と、更新された源信号のパワースペクトルの推定値と、雑音パラメータの推定値を用いて、残響パラメータの推定値を更新する。残響パラメータの推定値は、信号源パラメータの推定値と雑音パラメータの推定値とステアリングベクトルの推定値とが固定されている条件の下で、パラメータに関する対数尤度関数が最大になるように更新される。

[終了条件判定処理]
終了条件判定処理では、所定の終了条件が満たされているか否かが判定される。終了条件がを満たされていない場合、第１更新処理に戻る。終了条件が満たされている場合、その時点におけるパラメータの推定値を出力する。

〔原理〕
次に、本実施形態の原理を説明する。
本実施形態の信号強調装置の源信号推定部は、観測信号に含まれる残響を線形フィルタ処理で抑圧して雑音重畳信号を推定した後に、Wienerフィルタ等の非線形フィルタ処理により雑音重畳信号から雑音を抑圧する。この手順を実現するために、本実施形態のパラメータ推定部が生成するパラメータが第１，２実施形態のパラメータと異なる。

図２に模式的に示したように、時間領域の観測信号を生成する系は、複数の室内インパルス応答を畳み込む残響重畳系（室内伝達系）と、それぞれの残響重畳系の出力に定常雑音を加算する雑音重畳系とから成る。それらの系によって源信号に残響や雑音が付加され、時間領域の観測信号になる。時間周波数領域の観測信号ベクトルと源信号とを、それぞれｙ_ｔ，ｗ、Ｓ_ｔ，ｗとすると、両者の関係は式（98）で表せる。

ここで、ｄ_ｔ，ｗ＝[D_ｔ，ｗ ^（1），…，D_ｔ，ｗ ^（Ｍ）]^τは雑音ベクトル、ｂ_ｗはＭ次元のステアリングベクトル、Ｇ_ｋ，ｗを室内伝達系に関するｋ次の回帰行列、Ｈは共役転置、τは非共役転置を表す。式（98）は、室内伝達系がｗ番目の周波数帯域において、Ｇ_ｋ，ｗをｋ次の回帰行列にもつＫ_ｗ次のＭチャネル自己回帰系で表せることを意味している。式（99）は式（100）と式（101）に等価変換出来る。

式（101）に示すように、ｖ_ｔ，ｗは、０番目のタップ重み行列が単位行列でｋ番目（ｋ≧１）のタップ重み行列が−Ｇ_ｋ，ｗであるＭ入力Ｍ出力線形フィルタに、雑音ベクトルｄ_ｔ，ｗが入力され得られる出力信号である。すなわち、ｖ_ｔ，ｗは、フィルタ処理された雑音であり、源信号に由来する成分を含まない。本実施形態では、これを単に雑音と呼ぶ。また、式（100）に示すように、φ_ｔ，ｗは、源信号Ｓ_ｔ，ｗとＭ次元のステアリングベクトルｂ_ｗとの積と、雑音ベクトルｖ_ｔ，ｗとの和である。以後φ_ｔ，ｗを雑音重畳信号ベクトルと呼ぶ。また、式（99）に示すように、観測信号ベクトルｙ_ｔ，ｗは、ｋ次の回帰行列がＧ_ｋ，ｗである自己回帰系に雑音重畳信号φ_ｔ，ｗが入力されて得られる残響が重畳された信号である。

本実施形態では、残響パラメータ_gΘは、_gΘ={{G_k,w}_1≦k≦Kw}_0≦w≦N-1と定義される。また、ステアリングベクトルの集合_bΘ={b_w}_0≦w≦N-1も本実施形態におけるパラメータの一部である。さらに、源信号と雑音に関して、第１、２実施形態と同様に、以下の条件を仮定する。

《源信号のモデル》
源信号の短時間パワースペクトル密度はＰ次の全極型の関数で与えられる。すなわち、第ｔフレームにおける源信号のパワースペクトル密度は、式(102)で与えられる。

ω∈[−π，π]は角周波数、ａ_ｔ，ｋは線形予測係数、_ｓσ_ｔ ^２は予測残差パワーである。この信号源パラメータを用い、第ｔフレームの周波数帯域ｗにおける目的音短時間パワースペクトル_ｓλ_ｔ，ｗは、式(104)で表せる。

（ｔ_１，ｗ_１）≠（ｔ_２，ｗ_２）ならばＳ_{ｔ１，ｗ２}とＳ_{ｔ２，ｗ２}は統計的に独立である。源信号Ｓ_ｔ，ｗは、平均０、分散が源信号短時間パワースペクトル_ｓλ_ｔ，ｗに等しい複素正規分布に従う。すなわち、源信号Ｓ_ｔ，ｗの確率密度関数は式(105)で与えられる。

ただし、_ｓΘは、_sΘ={a_t,1,…,a_t,P,_sσ_t ²}_{0≦t≦Ｔ-1}で定義される信号源パラメータである。また、Ｎ{ｘ；μ，Σ}は、式(4)で定義される複素正規分布の確率密度関数である。
《雑音のモデル》
雑音は定常であると仮定すると、雑音の短時間パワースペクトル密度と短時間クロススペクトル密度は時不変である。すなわち、これらはフレーム番号ｔに依存しない。そこで、これらを式(106)のような行列で表現する。

ここで、_ｖλ^{（ｍ，ｍ）}（ω）はｍ番目のマイクロホンに関する雑音の短時間パワースペクトル密度、_ｖλ^{（ｍ1，ｍ2）}（ω）はｍ_１番目のマイクロホンに関する雑音とｍ_２番目のマイクロホンに関する雑音の間のクロススペクトル密度である。ｗ番目の周波数帯域における雑音短時間パワークロススペクトル行列_ｖΛ_ｗは、式(107)により与えられる。

（ｔ_１，ｗ_１）≠（ｔ_２，ｗ_２）ならばｖ_{ｔ１，ｗ1}とｖ_{ｔ２，ｗ２}も統計的に独立である。また、任意の（ｔ_１，ｗ_１，ｔ_２，ｗ_２）について、源信号Ｓ_{ｔ１，ｗ１}と雑音ベクトルｖ_{ｔ２，ｗ２}は統計的に独立である。
雑音ベクトルｖ_ｔ，ｗは、平均Ｏ _Ｍ＝[０，…，０]^τ、共分散行列が雑音短時間パワークロススペクトル行列_ｖΛ_ｗに等しいＭ次元複素正規分布に従う。すなわち、雑音ベクトルｖ_ｔ，ｗの確率密度関数は式（108）で与えられる。

ただし、_ＶΘは、_VΘ=｛_vΛ_w}_{0≦w≦Ｎ−１}で定義される雑音パラメータである。
したがって、本実施形態のパラメータΘは式（109）〜式（113）で定義される。

雑音と残響を含む観測信号が入力された時に、本実施形態のパラメータ推定部は、上記パラメータΘを最尤推定する。さらに、式(102)と式(103)と式(104)に従って、信号源パラメータの推定値から源信号パワースペクトルの推定値を計算する。これらの推定値が源信号推定部に供給される。

また、回帰行列の推定値をＧ_ｋ，ｗ＾、ステアリングベクトルの推定値をｂ_ｗ＾、線形予測係数の推定値をa_t,k＾、予測残差パワーの推定値を_sσ_t＾²、源信号短時間パワースペクトルの推定値を_ｓλ_ｔ，ｗ＾、雑音短時間パワークロススペクトル行列の推定値を_ｖΛ_ｗ＾とおく。
本実施形態の源信号推定部は、まず、式(114)に従って観測信号ベクトルｙ_ｔ，ｗから残響を抑圧して雑音重畳信号ベクトルの推定値残響抑圧信号φ_ｔ，ｗ＾を求める。

次に、源信号推定部は、残響抑圧信号φ_ｔ，ｗ＾に対して多チャネルWienerフィルタを用い、式(115)に示すように源信号Ｓ_ｔ，ｗの最小平均二乗誤差（ＭＭＳＥ）推定値を算出する。

ここでＦ（・）は多チャネルWienerフィルタのゲインベクトルである。
《パラメータの対数尤度関数》
上記した源信号及び雑音と、観測信号ベクトルの生成モデル式（99）と式（100）とに基づき、パラメータΘの対数尤度関数
Ｌ（Θ；ｙ）＝ｌｏｇ ｐ（ｙ｜Θ） (117)
は、式（118）で表せる。

ただし、_φΛ_ｔ，ｗは雑音重畳信号φ_ｔ，ｗの共分散行列を表し、式（119）で与えられる。

式（118）の導出過程を説明する。雑音重畳信号φ_ｔ，ｗの共分散行列が式（119）になることは、例えば参考文献「伊藤信貴他“結晶型マイクロホンアレイを用いたポストフィルタ設計に基づく拡散性雑音抑圧”信学技報ＥＡ2008-13，ｐｐ.43-46,2008」に記載されている。
これと式（99）により、過去の観測信号ベクトルが与えられた下での観測信号ベクトルｙ_ｔ，ｗの条件付確率密度関数が、式（120）で与えられることが分る。

したがって、すべての観測信号ベクトルの集合ｙについての確率密度関数は式（121）で表せる。ただし、ｙ＝{ｙ_ｔ，ｗ}_{０≦ｔ≦Ｔ-1，０≦ｗ≦Ｎ-1}である。

式（121）の両辺の対数を取ることで対数尤度関数、式（118）が導かれる。
＜本実施形態の構成及び処理＞
図９は、第３実施形態の信号強調装置２００の機能構成例を示すブロック図である。図１０は、第３実施形態の処理を説明するためのフローチャートである。

本実施形態の信号強調装置２００は、帯域分割部２２０と、パラメータ推定部３１０と、源信号推定部２３０と、制御部２５０と、帯域合成部２４０と、を有する。源信号推定部２３０は、線形フィルタ処理部２３１と非線形フィルタ処理部２３２とを含む。帯域分割部２２０と帯域合成部２４０とは、第１，２実施形態のものと同じである。信号強調装置２００は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現される専用装置である。

帯域分割部２２０は、時間領域の観測信号を所定数の周波数帯域毎の観測信号ベクトルｙ_ｔ，ｗ（０≦ｔ≦Ｔ−１，０≦ｗ≦Ｎ−１）に分割する（ステップＳ２０１）。パラメータ推定部３１０は、入力された観測信号ベクトルｙ_ｔ，ｗを用いて、残響を推定するための回帰行列Ｇ_ｋ，ｗを含む残響パラメータ_gΘと、源信号を推定するための雑音短時間パワークロススペクトル行列_ｖΛ_ｗを含む雑音パラメータ_vΘと、源信号短時間パワースペクトル_ｓλ_ｔ，ｗを規定する信号源パラメータ_sΘと、ステアリングベクトルｂ_ｗの集合_bΘの各真値をそれぞれ推定する（ステップＳ２０２）。

＜ステップＳ２０２の詳細＞
図１１は、第３実施形態のパラメータ推定部３１０の機能構成例を示すブロック図である。また、図１２は、第３実施形態のパラメータ推定処理を説明するためのフローチャートである。本実施形態のパラメータ推定部３１０は、未知のパラメータΘを最尤推定するために残響パラメータ_ｇΘ、ステアリングベクトル_ｂΘ、信号源パラメータ_ｓΘ、雑音パラメータ_ｖΘのそれぞれの推定値を繰り返し更新する。

パラメータ推定部３１０は、観測信号記録部３１１と、パラメータ推定値初期化部３１２（「初期化部」に相当）と、源信号推定値更新部３１３と、信号源パラメータ推定値更新部３１４と、源信号パワースペクトル推定値更新部３１５と、残響パラメータ推定値更新部３１６と、ステアリングベクトル推定値更新部３１８と、雑音パラメータ推定値更新部３１９と、収束判定部３１７とを有する。

源信号推定値更新部３１３と、ステアリングベクトル推定値更新部３１８と、信号源パラメータ推定値更新部３１４とは、前述した第１更新部に含まれる。また、源信号パワースペクトル推定値更新部３１５と、雑音パラメータ推定値更新部３１９と、残響パラメータ推定値更新部３１６とは、前述した第２更新部に含まれる。

観測信号記録部３１１は、帯域分割部２２０で所定数の周波数帯域に分割された観測信号を記録する。観測信号記録部３１１は、観測区間中のすべての雑音残響重畳信号を記録する。そして、観測信号記録部３１１は、記録した観測信号を源信号推定値更新部３１３と残響パラメータ推定値更新部３１６とパラメータ推定値初期化部３１２とに出力する。

パラメータ推定値初期化部３１２は、入力された観測信号ベクトルｙ_ｔ，ｗを用いて、残響パラメータ_ｇΘ、ステアリングベクトル_ｂΘ、信号源パラメータ_ｓΘ、雑音パラメータ_ｖΘの各初期値を設定する。また、制御部２５０が、繰り返し回数を示すインデックスｉを０にする。

源信号推定値更新部３１３は、入力された観測信号ベクトルｙ_ｔ，ｗと、各パラメータの推定値の初期値_ｇΘ^（０）＾，_ｂΘ^（０）＾，_ｓΘ^（０）＾，_ｖΘ^（０）＾又は更新された各パラメータの推定値_ｇΘ^（ｉ）＾，_ｂΘ^（ｉ）＾，_ｓΘ^（ｉ）＾，_ｖΘ^（ｉ）＾を用いて、源信号の推定値Ｓ_ｔ，ｗ ^（ｉ）＾とその誤差分散と、雑音重畳信号の推定値φ_ｔ，ｗ ^（ｉ）＾を、それぞれＳ_ｔ，ｗ ^{（ｉ＋１）}＾とその誤差分散とφ_ｔ，ｗ ^{（ｉ＋１）}＾に更新する（ステップＳ３０１）。Ｓ_ｔ，ｗ ^{（ｉ＋１）}＾は式(115)を用い、φ_ｔ，ｗ ^{（ｉ＋１）}＾は式(114)を用いて計算される。誤差分散は式（122）を用いて計算される。

ステアリングベクトル推定値更新部３１８には、更新された源信号の推定値Ｓ_ｔ，ｗ ^{（ｉ＋１）}＾と、雑音重畳信号の推定値φ_ｔ，ｗ ^{（ｉ＋１）}＾とが入力される。ステアリングベクトル推定値更新部３１８は、これらを用い、式（123）に従って、更新されたステアリングベクトルの推定値を計算する。式（123）は、雑音ベクトルの平均がＯ_Ｍであるとの仮定に基づいている。

ここで、＊は複素共役を表す。すべての周波数帯域ｗ(０≦ｗ≦Ｎ−１)に渡って式（123）が計算されることで、更新されたステアリングベクトルの推定値_ｂΘ^{（ｉ＋１）}＾が得られる（ステップＳ３０３）。
信号源パラメータ推定値更新部３１４は、源信号の推定値Ｓ_ｔ，ｗ ^{（ｉ＋１）}＾のパワーとその誤差分散ε_ｔ，ｗ ^{（ｉ＋１）}を式（124）に示すように加算してパワースペクトルγ_ｔ，ｗ ^{（ｉ＋１）}を求める。

そして、信号源パラメータ推定値更新部３１４は、求めたパワースペクトルγ_ｔ，ｗ ^{（ｉ＋１）}を用い、Levinson-Durbinアルゴリズムによって、信号源パラメータの推定値を更新する。Levinson-Durbinアルゴリズムは周知の方法であるので詳細な説明は省略するが、式(40)のV_t,w ⁽ⁱ⁾をγ_ｔ，ｗ ^{（ｉ＋１）}に置換し、式(36)から(40)の演算を行うことで、更新された信号源パラメータ（ａ_ｔ，１ ^{（ｉ＋１）}＾，…，ａ_ｔ，Ｐ ^{（ｉ＋１）}＾，_ｓσ_ｔ ^{２（ｉ＋１）}＾）が算出される。そして、すべてのフレーム番号ｔ（０≦ｔ≦Ｔ−１）に渡ってこれらが計算されることで、更新された信号源パラメータ_ｓΘ^{（ｉ＋１）}＾が得られる（ステップＳ３０４）。

源信号パワースペクトル推定値更新部３１５には、更新された信号源パラメータの推定値が入力される。源信号パワースペクトル推定値更新部３１５は、更新された信号源パラメータを用い、源信号の短時間パワースペクトルの推定値を更新する（ステップＳ３０５）。源信号の短時間パワースペクトルの更新された推定値_ｓλ_ｔ，ｗ ^{（ｉ＋１）}＾は、式(102)と式(103)と式(104)を用いて計算される。

雑音パラメータ推定値更新部３１９には、更新された源信号の推定値Ｓ_ｔ，ｗ ^{（ｉ＋１）}＾と雑音重畳信号の推定値φ_ｔ，ｗ ^{（ｉ＋１）}＾とステアリングベクトルの更新値_ｂΘ^{（ｉ＋１）}＾とが入力される。雑音パラメータ推定値更新部３１９は、これらを用い、式（125）に従って、雑音短時間パワークロススペクトル行列の推定値_ｖΛ_ｗ ^{（ｉ＋１）}を、すべての周波数帯域ｗ（０≦ｗ≦Ｎ−１）に渡って計算する。

ここで、Ｔ′は十分小さい値であり、ｔ＝０からｔ＝Ｔ′−１までの区間は、観測信号の冒頭部分である。本実施形態では、冒頭部分のＴ′フレーム（例えば０．３秒間）は雑音のみを含むものと仮定し、その区間に対する計算結果から雑音短時間パワークロススペクトル行列の推定値_ｖΛ_ｗ ^{（ｉ＋１）}＾を更新する（ステップＳ３０６）。

残響パラメータ推定値更新部３１６は、入力された観測信号ベクトルｙ_ｔ，ｗと、更新されたステアリングベクトルの推定値_ｂΘ^{（ｉ＋１）}＾と、源信号短時間パワースペクトルの推定値_ｓλ_ｔ，ｗ ^{（ｉ＋１）}＾と、雑音短時間パワークロススペクトル行列の推定値_ｖΛ_ｗ ^{（ｉ＋１）}＾とを用い、残響パラメータの更新された推定値_ｇΘ^{（ｉ＋１）}＾を求める（ステップＳ３０７）。残響パラメータ推定値更新部３１６は、まず、ｗ番目の周波数帯域における回帰行列の各成分を、式（126）と式（127）に示すように単一のベクトルにまとめる。

式（126）と式（127）の右下の添え字は、それぞれの式が示す行列（あるいはベクトル）の大きさを表す。ここで、ｇ_ｋ，ｗ ^（ｍ）は回帰行列Ｇ_ｋ，ｗのｍ番目の列を表すものとする。以後ｇ_ｗを回帰行列の成分ベクトルと呼ぶ。成分ベクトルｇ_ｗの全周波数帯域に渡る集合{ｇ_ｗ}_{０≦ｗ≦Ｗ−１}は残響パラメータ_ｇΘに一致する。
次に、１フレーム前の観測信号行列ＭＹ_{ｔ−１，ｗ}を式（128）のように定義する。

これらを用い、式（130）に従って、回帰行列の成分ベクトルの更新後の推定値ｇ_ｗ ^{（ｉ＋１）}＾が計算される。

ここで、_φΛ_ｔ，ｗ ^{（ｉ＋１）}＾は式（119）でｂ_ｗ＝ｂ_ｗ ^{（ｉ＋１）}＾，_ｓλ_ｔ，ｗ＝_ｓλ_ｔ，ｗ ^{（ｉ＋１）}＾，_ｖΛ_ｗ＝_ｖΛ_ｗ ^{（ｉ＋１）}＾として得られる値である。すべての周波数帯域ｗ（０≦ｗ≦Ｎ−１）に渡ってこれらが計算されることで残響パラメータの推定値の更新値_ｇΘ^{（ｉ＋１）}＾が得られる。

次に、以上のように更新された残響パラメータの推定値_ｇΘ^{（ｉ＋1）}＾と、ステアリングベクトルの推定値_ｂΘ^{（ｉ＋1）}＾と、信号源パラメータの推定値_ｓΘ^{（ｉ＋１）}＾と、雑音パラメータ_ｖΘ^{（ｉ＋1）}＾とが、収束したか否か（終了条件を充足したか否か）を、収束判定部３１７が判定する（ステップＳ３０８）。例えば、収束判定部３１７は、繰り返し回数ｉが所定数に到達していれば収束していると判定しても良いし、上述の処理が繰り返されるたびに得られる対数尤度関数（式（118））の値の増分が、所定の閾値よりも小さければ収束していると判定しても良い。これらの値が収束するまでステップＳ３０２〜ステップＳ３０７の動作が繰り返され、所定の終了条件が満たされた場合、その時点での残響パラメータの推定値_ｇΘ＾^{（ｉ＋1）}と、ステアリングベクトルの推定値_ｂΘ^{（ｉ＋1）}＾、信号源パラメータの推定値_ｓΘ^{（ｉ＋１）}＾、雑音パラメータ_ｖΘ^{（ｉ＋1）}＾とが、源信号推定部２３０に出力される。この際、パラメータ推定値記録部３２０にこのパラメータの推定値が記録されても良い（ステップＳ２０２の詳細の説明終わり）。

線形フィルタ処理部２３１は、観測信号ベクトルｙ_ｔ，ｗに回帰行列の推定値Ｇ_ｋ，ｗ＾を畳み込み演算して残響を求める。そして、線形フィルタ処理部２３１は、求めた残響を観測信号ベクトルから減算して残響抑圧信号ベクトルφ_ｔ，ｗ＾を生成する（ステップＳ２０３）。非線形フィルタ処理部２３２は、入力された雑音短時間パワークロススペクトル行列の推定値_ｖΛ_ｗ＾と源信号短時間パワースペクトルの推定値_ｓλ_ｔ，ｗ＾とステアリングベクトルの推定値ｂ_ｗ＾と残響抑圧信号φ_ｔ，ｗ＾とを用いて、残響抑圧信号φ_ｔ，ｗ＾から雑音を抑圧した源信号の推定値ｓ_ｔ，ｗ＾を生成する（ステップＳ２０４）。帯域合成部２４０は、源信号の推定値Ｓ_ｔ，ｗ＾を合成して時間領域の源信号の推定値に変換する（ステップＳ２０５）。制御部２５０は、入力される時間領域の観測信号から、残響と雑音が抑圧された時間領域の源信号の推定値が生成されるように、上記各処理部を制御する。

以上のように信号強調装置２００では、線形フィルタ処理部２３１が観測信号ベクトルｙ_ｔ，ｗに含まれる残響を抑圧して残響抑圧信号ベクトルφ_ｔ，ｗ＾を生成し、その後に非線形フィルタ処理部２３２が残響抑圧信号から雑音を抑圧する。この時間領域の源信号の推定値は、観測信号ベクトルを線形フィルタ処理した後に非線形フィルタ処理して得られたものである。そのため、この時間領域の源信号の推定値は、雑音と残響とが十分抑圧された高品質な信号である。

なお、上記では、回帰次数（線形フィルタのフィルタ長）Ｋ_ｗを一つの固定値として説明した。しかし、回帰次数が、周波数帯域の中心周波数に応じて変化しても良い。周波数帯域によって残響時間が異なることは良く知られている。例えば、室内音響の分野においては、５００Ｈｚ以下の周波数帯域の残響時間が長いので、その周波数帯域では回帰次数Ｋ_ｗを大きくし、それ以外の周波数帯域では回帰次数Ｋ_ｗを小さくしてもよい。また、パラメータ推定部３１０内に回帰次数可変部３０１を備え、回帰次数可変部３０１が、周波数帯域に応じて回帰次数、つまり、線形フィルタ処理部２３１のフィルタ長を変化させてもよい。これにより、残響を効率的に抑圧することが可能になる。つまり、線形フィルタ処理部２３１の計算量を削減できる。このような変形は、前述の第１，２実施形態でも可能である。

〔実験結果〕
本実施形態の信号強調方法の効果を確認する目的で実験を行った。実験条件を説明する。源信号には、ＡＳＪ-ＪＮＡＳデータベースから抽出した１０名（男性５名、女性５名）による発話を用いた。これらの音声を残響時間が約０．６秒の部屋でスピーカーから再生し、スピーカーから１．８ｍ離して設置した２個のマイクロホンで録音した。また、同じ部屋、同じマイクロホンで、４箇所に設置したスピーカーから同時に再生したピンクノイズを録音した。その後、録音された残響音声と雑音をＳＮ比が１０ｄＢとなるように加算したものを時間領域の観測信号として用いた。なお、録音時の標本化周波数は８ｋＨｚとした。

本実施形態の 帯域分割部の処理には、ポリフェーズフィルタバンク分析を用いた。帯域分割数は２５６、間引き率は１２８とした。
源信号の線形予測次数はＰ＝１２とした。回帰次数Ｋ_ｗは、観測信号の周波数が１００Ｈｚ未満ならばＫ_ｗ＝５、１００Ｈｚ〜２００ＨｚならばＫ_ｗ＝１０、２００Ｈｚ〜１０００ＨｚならばＫ_ｗ＝３０、１０００Ｈｚ〜１５００ＨｚならばＫ_ｗ＝２０、１５００Ｈｚ〜２０００ＨｚならばＫ_ｗ＝１５、２０００Ｈｚ〜３０００ＨｚならばＫ_ｗ＝１０、３０００Ｈｚ以上ならばＫ_ｗ＝５とした。また、収束判定部は、繰り返し回数が３回で収束したと判定する。

以上の条件で、観測信号そのまま、実施形態１による源信号の推定値、本実施形態による源信号の推定値、のそれぞれの源信号からのＭＦＣＣ距離の平均値を比較した。その結果は、順番に７.３９、５.８１、５.１１であった。このようにこの発明の信号強調方法によるＭＦＣＣ距離が最も近いという結果が得られた。

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

また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。

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

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

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。

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

本発明の利用分野としては、例えば、音声認識システムやテレビ会議システム等での源音声信号の強調処理を例示できる。

その後、残響パラメータの推定値_gΘ^⁽ⁱ⁾を用い、信号源パラメータの推定値を_sΘ^⁽ⁱ⁺¹⁾に更新する（ステップＳ１１３）。次に、更新された信号源パラメータの推定値_sΘ^⁽ⁱ⁺¹⁾を用い、残響パラメータの推定値を_gΘ^⁽ⁱ⁺¹⁾ に更新する（ステップＳ１１４）。そして、所定の終了条件を満たすまで（ステップＳ１１５）、ｉを１ずつ増加させながら（ステップＳ１１６）、ステップＳ１１３とＳ１１４との処理を繰り返し、所定の終了条件が満たされた時点における信号源パラメータの推定値_sΘ^〜(i+1)を信号源パラメータの最終的な推定値_sΘ^とし、残響パラメータの推定値_gΘ^⁽ⁱ⁺¹⁾を残響パラメータの最終的な推定値_gΘ^として出力する（ステップＳ１１７）。

その後、源信号推定部が、パラメータ推定部で計算された残響パラメータの最終的な推定値_gΘ^を用いて生成した線形フィルタを観測信号に畳み込むことで観測信号に含まれる残響を推定し、それを観測信号から減算することで、残響が抑圧された信号を算出して出力する。
Lim, J. S. and Oppenheim, A. V. , "All-pole modeling of degraded speech," IEEE Trans. Acoust. Speech, Signal Process., Vol. 26, No. 3, pp.197-210 (1978). Yoshioka, T., Hikichi, T. and Miyoshi, M., "Dereverberation by Using Time-Variant Nature of Speech Production System, EURASIP J. Advances in Signal Process., Vol. 2007, (2007), Article ID 65698, 15 pages, doi:10.1155/2007/65698.

以上のように、本発明のパラメータ推定部では、第１更新部におけるパラメータの推定値の更新処理と、第２更新部におけるパラメータの推定値の更新処理を、互いに依存させながら繰り返して実行する。これにより、雑音と残響がともに存在する環境における観測信号から、雑音と残響を精度よく抑圧し、源信号を強調することができる。

終了判定条件部は、所定の終了条件が満たされているか否かを判定する。終了条件が満たされていない場合、第１更新部の処理に戻る。終了条件が満たされている場合、その時点におけるパラメータの推定値を出力する。

［信号源パラメータ及び残響パラメータの最尤推定］
前述のように、本実施形態では、観測された雑音残響重畳信号の集合のyから、未知のパラメータの真値Θ^〜が、最尤推定法によって推定される。すなわち、雑音残響重畳信号の集合yが与えられた場合におけるパラメータΘを変数とした尤度関数p(y|Θ)を最大化するΘが、真値Θ^〜の推定値となる。ただし、本実施形態では、雑音パラメータの真値_dΘ^〜が源信号の存在しない区間から予め独立に推定され、既知となっている為Θ^={_sΘ^,_gΘ^, _dΘ^〜}であり、_sΘ^と_gΘ^が求められることになる。

また、尤度関数p(y|Θ)を最大化する_sΘ^と_gΘ^を同時に直接求めることはできないから、ＥＣＭアルゴリズムを用いてこれらが計算される。ＥＣＭアルゴリズムの処理の流れを以下に示す。以下の処理では、Ｅ−ｓｔｅｐ、ＣＭ−ｓｔｅｐ１、ＣＭ−ｓｔｅｐ２の３つの処理が代わる代わる繰り返し実行される。そこで、ｉ回目の繰り返しにおけるパラメータの推定値を上付きの添え字(i)を用いて示す。明確さを期するために述べると、Θ^〜，Θ^，Θ^⁽ⁱ⁾はそれぞれ次のように定義される。

《源信号短時間パワースペクトル推定値更新処理》
源信号短時間パワースペクトル推定値更新処理では、信号源パラメータ推定値を用いて源信号のパワースペクトルの推定値を更新する。

ここで、ｄ_ｔ，ｗ＝[D_ｔ，ｗ ^（1），…，D_ｔ，ｗ ^（Ｍ）]^τは雑音ベクトル、ｂ_ｗはＭ次元のステアリングベクトル、Ｇ_ｋ，ｗを室内伝達系に関するｋ次の回帰行列、Ｈは共役転置、τは非共役転置を表す。式（98）は、室内伝達系がｗ番目の周波数帯域において、Ｇ_ｋ，ｗをｋ次の回帰行列にもつＫ_ｗ次のＭチャネル自己回帰系で表せることを意味している。式（98）は式（99）〜式（101）に等価変換出来る。

ω∈[−π，π]は角周波数、ａ_ｔ，ｋは線形予測係数、_ｓσ_ｔ ^２は予測残差パワーである。この信号源パラメータを用い、第ｔフレームの周波数帯域ｗにおける源信号短時間パワースペクトル_ｓλ_ｔ，ｗは、式(104)で表せる。

また、回帰行列の推定値をＧ_ｋ，ｗ＾、ステアリングベクトルの推定値をｂ_ｗ＾、線形予測係数の推定値をa_t,k＾、予測残差パワーの推定値を_sσ_t＾²、源信号短時間パワースペクトルの推定値を_ｓλ_ｔ，ｗ＾、雑音短時間パワークロススペクトル行列の推定値を_ｖΛ_ｗ＾とおく。
本実施形態の源信号推定部は、まず、式(114)に従って観測信号ベクトルｙ_ｔ，ｗから残響を抑圧して雑音重畳信号ベクトルの推定値（残響抑圧信号）φ_ｔ，ｗ＾を求める。

雑音パラメータ推定値更新部３１９には、更新された源信号の推定値Ｓ_ｔ，ｗ ^{（ｉ＋１）}＾と雑音重畳信号の推定値φ_ｔ，ｗ ^{（ｉ＋１）}＾とステアリングベクトルの更新値_ｂΘ^{（ｉ＋１）}＾とが入力される。雑音パラメータ推定値更新部３１９は、これらを用い、式（125）に従って、雑音短時間パワークロススペクトル行列の推定値 _ｖ Λ _ｗ ^{（ｉ＋１）} ＾を、すべての周波数帯域ｗ（０≦ｗ≦Ｎ−１）に渡って計算する。

式（126）と式（127）の右下の添え字は、それぞれの式が示す行列（あるいはベクトル）の大きさを表す。ここで、ｇ_ｋ，ｗ ^（ｍ）は回帰行列Ｇ_ｋ，ｗのｍ番目の列を表すものとする。以後ｇ_ｗを回帰行列の成分ベクトルと呼ぶ。成分ベクトルｇ_ｗの全周波数帯域に渡る集合{ｇ _ｗ } _{０≦ｗ≦N−１} は残響パラメータ_ｇΘに一致する。
次に、１フレーム前の観測信号行列ＭＹ_{ｔ−１，ｗ}を式（128）のように定義する。

Claims (17)

1. 観測された時間領域信号から変換された時間周波数領域の観測信号を格納する記憶部と、
   前記観測信号に含まれる残響の推定値を算出する線形畳み込み演算の回帰係数を含む残響パラメータ推定値と、源信号のパワースペクトルを特定する線形予測係数と予測残差パワーとの推定値を含む信号源パラメータ推定値と、雑音のパワースペクトルの推定値を含む雑音パラメータ推定値と、を含むパラメータ推定値の初期値を設定する初期化部と、
   前記観測信号と前記パラメータ推定値とが入力され、前記残響パラメータ推定値および雑音パラメータ推定値の少なくとも一部の更新処理、あるいは前記信号源パラメータ推定値の更新処理、のいずれか一方を実行するように構成され、当該更新処理が前記パラメータ推定値に関する対数尤度関数の値が増加するように実行される処理である、第１更新部と、
   前記第１更新部で得られたパラメータ推定値の更新値の少なくとも一部が入力され、前記残響パラメータ推定値および雑音パラメータ推定値の少なくとも一部の更新処理、あるいは前記信号源パラメータ推定値の更新処理のうち、前記第１更新部で実行されなかったものを実行するように構成され、当該更新処理が前記パラメータ推定値の更新値に関する対数尤度関数の値が増加するように実行される処理である、第２更新部と、
   終了条件が満たされるか否かを判定する終了条件判定部と、を有し、
   前記終了条件が満たされない場合、前記第１更新部と前記第２更新部の処理が再び実行される、信号強調装置。
2. 請求項１の信号強調装置であって、
   前記時間領域信号が、Ｍ個のセンサで観測された信号であり、
   前記残響パラメータ推定値が、前記回帰係数を要素にもつＭ行Ｍ列の回帰行列推定値を含み、
   前記雑音パラメータ推定値が、前記雑音のパワースペクトルを対角要素とするＭ行Ｍ列の雑音パワークロススペクトル行列推定値を含み、
   前記パラメータ推定値が、前記残響パラメータ推定値と、前記信号源パラメータ推定値と、前記雑音パラメータ推定値と、Ｍ次元のステアリングベクトル推定値と、を含み、
   前記第１更新部が、
   源信号推定値更新部と、ステアリングベクトル推定値更新部と、信号源パラメータ推定値更新部と、を含み、
   前記源信号推定値更新部は、前記観測信号と前記パラメータ推定値とが入力され、雑音重畳信号推定値と、源信号推定値と、前記源信号推定値の誤差分散とを算出するように構成され、
   前記ステアリングベクトル推定値更新部は、前記雑音重畳信号推定値と前記源信号推定値とが入力され、ステアリングベクトル推定値の更新値を算出するように構成され、
   前記信号源パラメータ推定値更新部は、前記源信号推定値のパワーと前記誤差分散とを加算してパワースペクトルを算出し、前記パワースペクトルを用いて信号源パラメータ推定値の更新値を算出するように構成され、
   前記第２更新部が、源信号パワースペクトル推定値更新部と、雑音パラメータ推定値更新部と、残響パラメータ推定値更新部と、を含み、
   前記源信号パワースペクトル推定値更新部は、前記信号源パラメータ推定値の更新値が入力され、前記信号源パラメータ推定値の更新値に対応する源信号パワースペクトル推定値の更新値を算出するように構成され、
   前記雑音パラメータ推定値更新部は、前記源信号推定値と、前記雑音重畳信号推定値と、前記ステアリングベクトル推定値の更新値とが入力され、前記雑音パラメータ推定値の更新値を生成するように構成され、
   前記残響パラメータ推定値更新部は、前記観測信号と、前記ステアリングベクトル推定値の更新値と、前記源信号パワースペクトル推定値の更新値と、前記雑音パラメータ推定値の更新値とが入力され、前記回帰行列推定値の更新値を算出するように構成される、
   信号強調装置。
3. 請求項２の信号強調装置であって、
   前記雑音パワークロススペクトル行列推定値のｍ行ｍ列（ｍ∈１,...,Ｍ）の要素が、ｍ番目のセンサに対応する前記雑音のパワースペクトルであり、前記雑音パワークロススペクトル行列推定値のｍ１行ｍ２列(ｍ１,ｍ２∈1,...,M)の要素が、ｍ１番目のセンサに対応する前記観測信号の雑音と、ｍ２番目のセンサに対応する前記観測信号の雑音との間のクロススペクトルであり、
   前記雑音重畳信号推定値が、それぞれの要素が各センサに対応する前記観測信号であるＭ次元ベクトルの非共役転置である観測信号ベクトルから、前記回帰行列推定値と前記観測信号ベクトルとの畳み込み演算結果を減じたＭ次元ベクトルであり、
   前記源信号推定値が、前記源信号パワースペクトル推定値と前記雑音パワークロススペクトル行列推定値と前記ステアリングベクトル推定値とに対応するWienerフィルタのゲインベクトルと、前記雑音重畳信号推定値と、の積であり、
   前記源信号推定値の誤差分散が、前記ステアリングベクトル推定値の非共役転置と前記雑音パワークロススペクトル行列推定値の逆行列と前記ステアリングベクトル推定値との積と、前記信号源パラメータ推定値に対応する源信号パワースペクトル推定値の逆数と、の加算値の逆数であり、
   前記ステアリングベクトル推定値の更新値が、前記源信号推定値の複素共役値と前記雑音重畳信号推定値との積和を、前記源信号推定値のパワーの積和で割ったベクトルであり、
   前記雑音パワークロススペクトル行列推定値の更新値が、前記雑音重畳信号推定値から前記源信号推定値と前記ステアリングベクトル推定値の更新値との積を減じた雑音ベクトルと、当該雑音ベクトルの共役転置との積和であり、
   前記回帰行列推定値の更新値の要素からなる成分ベクトルが、前記観測信号を要素とする観測信号行列の共役転置と雑音重畳信号の共分散行列の推定値の逆行列と前記観測信号行列との積和の逆行列と、前記観測信号行列の共役転置と雑音重畳信号の共分散行列の推定値の逆行列と前記観測信号ベクトルとの積和と、の積の共役転置であり、
   前記雑音重畳信号の共分散行列の推定値が、前記源信号パワースペクトル推定値の更新値と前記ステアリングベクトル推定値の更新値と前記ステアリングベクトル推定値の更新値の共役転置との積と、前記雑音パワークロススペクトル行列推定値の更新値との和である、信号強調装置。
4. 請求項２の信号強調装置であって、
   前記残響パラメータ推定値又はその更新値に含まれる回帰行列推定値の回帰次数が、周波数帯域によって異なる、信号強調装置。
5. 請求項２の信号強調装置であって、
   前記観測信号と残響パラメータ最終推定値とが入力され、前記観測信号ベクトルから、前記残響パラメータ最終推定値と前記観測信号との畳み込み演算結果を減じたＭ次元ベクトルである雑音重畳信号最終推定値を生成する線形フィルタ処理部と、
   信号源パラメータ最終推定値によって特定される源信号パワースペクトル最終推定値と、雑音パラメータ最終推定値に含まれる雑音パワークロススペクトル行列最終推定値と、ステアリングベクトル最終推定値と、前記雑音重畳信号最終推定値とが入力され、前記源信号パワースペクトル最終推定値と前記雑音パワークロススペクトル行列最終推定値と前記ステアリングベクトル最終推定値とに対応するWienerフィルタのゲインベクトルと、前記雑音重畳信号最終推定値と、の積を源信号最終推定値とする非線形フィルタ処理部と、を有し、
   前記残響パラメータ最終推定値、前記信号源パラメータ最終推定値、前記雑音パラメータ最終推定値、及び前記ステアリングベクトル最終推定値が、前記終了条件を満たした時点における前記回帰行列推定値の更新値、前記信号源パラメータ推定値の更新値、前記雑音パラメータ最終推定値の更新値、及び前記ステアリングベクトル推定値の更新値を含む、信号強調装置。
6. 請求項１の信号強調装置であって、
   前記観測信号が１個のセンサで観測された信号であり、
   前記残響パラメータ推定値が、前記回帰係数の推定値を含み、
   前記雑音パラメータ推定値が、前記雑音のパワースペクトルの推定値を含み、
   前記パラメータ推定値が、前記信号源パラメータ推定値と、前記残響パラメータ推定値と、前記雑音パラメータ推定値と、を含み、
   前記第１更新部が、
   雑音抑圧処理部と、信号源パラメータ推定値更新部と、を含み、
   前記雑音抑圧処理部は、
   前記観測信号と前記パラメータ推定値とが入力され、所定の観測区間に属する前記観測信号の集合と前記パラメータ推定値との組合せを前提条件とした前記観測区間に属する残響重畳信号の集合の条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）を特定する複素正規分布の平均及び共分散行列を算出するように構成され、
   前記残響重畳信号が、前記観測信号から雑音が取り除かれた信号であり、
   前記信号源パラメータ推定値更新部は、
   前記残響パラメータ推定値と、前記複素正規分布の平均及び共分散行列とが入力され、 信号源パラメータ推定値の更新値を算出するように構成され、
   前記信号源パラメータ推定値の更新値は、残響パラメータが前記残響パラメータ推定値に固定された条件下で、第１補助関数値を最大化する値であり、
   前記第１補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第２パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第２パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値であり、
   前記第２更新部が、
   前記信号源パラメータ推定値の更新値と、前記複素正規分布の平均及び共分散行列とが入力され、残響パラメータ推定値の更新値を算出するように構成された残響パラメータ推定値更新部を含み、
   前記残響パラメータ推定値の更新値は、信号源パラメータが前記信号源パラメータ推定値の更新値に固定された条件下で、第２補助関数値を最大化する値であり、
   前記第２補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値の更新値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第３パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第３パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値である、信号強調装置。
7. 請求項１の信号強調装置であって、
   前記時間領域信号が、Ｍ個のセンサで観測された信号であり、Ｍが２以上であり、
   前記残響パラメータ推定値が、前記回帰係数を要素にもつＭ行Ｍ列の回帰行列推定値を含み、
   前記雑音パラメータ推定値が、前記雑音のパワースペクトルの推定値を対角要素とする、Ｍ行Ｍ列の雑音パワークロススペクトル行列推定値を含み、
   前記パラメータ推定値が、前記信号源パラメータ推定値と、前記残響パラメータ推定値と、前記雑音パラメータ推定値と、を含み、
   前記第１更新部が、雑音抑圧処理部と、信号源パラメータ推定値更新部と、を含み、
   前記雑音抑圧処理部は、
   前記観測信号と前記パラメータ推定値とが入力され、所定の観測区間に属する前記観測信号の集合と前記パラメータ推定値との組合せを前提条件とした前記観測区間に属する前記残響重畳信号の集合の条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）を特定する複素正規分布の平均及び共分散行列を算出するように構成され、
   前記残響重畳信号が、前記観測信号から雑音が取り除かれた信号であり、
   前記信号源パラメータ推定値更新部は、
   前記残響パラメータ推定値と、前記複素正規分布の平均及び共分散行列とが入力され、信号源パラメータ推定値の更新値を算出するように構成され、
   前記信号源パラメータ推定値の更新値は、残響パラメータが前記残響パラメータ推定値に固定された条件下で、第１補助関数値を最大化する値であり、
   前記第１補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第２パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第２パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値であり、
   前記第２更新部が、
   前記信号源パラメータ推定値の更新値と、前記複素正規分布の平均及び共分散行列とが入力され、残響パラメータ推定値の更新値を算出するように構成された残響パラメータ推定値更新部を含み、
   前記残響パラメータ推定値の更新値は、信号源パラメータが信号源パラメータ推定値の更新値に固定された条件下で、第２補助関数値を最大化する値であり、
   前記第２補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値の更新値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第３パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第３パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値である、信号強調装置。
8. 請求項６又は７の信号強調装置であって、
   前記雑音パラメータ推定値は、前記雑音の確率分布を示す複素正規分布の分散である、前記雑音のパワースペクトルの推定値を含み、前記残響重畳信号の集合の条件付事後分布ｐ（残響重畳信号の集合｜観測信号，パラメータ推定値）の共分散行列のスケールは、前記雑音の確率分布を示す複素正規分布の分散に対して単調増加する値である、信号強調装置。
9. 請求項６又は７の信号強調装置であって、
   前記観測信号と、前記終了条件を満たした場合の前記第３パラメータ推定値とが入力され、前記源信号の推定値を生成する源信号推定部を有し、
   前記源信号推定部は、
   前記観測信号と、前記終了条件を満たした場合の前記第３パラメータ推定値とが入力され、前記残響重畳信号の集合の条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の平均を残響重畳信号最終推定値として算出する残響重畳信号推定部と、
   前記残響重畳信号最終推定値と、前記終了条件を満たした場合の前記第３パラメータ推定値が含む前記第２残響パラメータ推定値とが入力され、前記残響重畳信号最終推定値から、前記残響重畳信号最終推定値と当該第２残響パラメータ推定値に含まれる回帰係数又は回帰行列との畳み込み演算結果を減じ、源信号最終推定値を生成する線形フィルタ適用部と、を有する、信号強調装置。
10. 請求項６又は７の信号強調装置であって、
    前記雑音成分のパワースペクトルの推定値が、前記源信号が存在しないと推定される区間の前記観測信号から推定された値である、信号強調装置。
11. 請求項６又は７の信号強調装置であって、
    前記残響パラメータ推定値及び前記残響パラメータ推定値の更新値に含まれる回帰行列推定値の回帰次数が、周波数帯域によって異なる、信号強調装置。
12. (A) 観測された時間領域信号から変換された時間周波数領域の観測信号を記録部に格納するステップと、
    (B) 初期化部において、前記観測信号に含まれる残響の推定値を算出する線形畳み込み演算の回帰係数を含む残響パラメータ推定値と、源信号のパワースペクトルを特定する線形予測係数と予測残差パワーとの推定値を含む信号源パラメータ推定値と、雑音のパワースペクトルの推定値を含む雑音パラメータ推定値と、を含むパラメータ推定値の初期値を設定するステップと、
    (C) 前記観測信号と前記パラメータ推定値とを第１更新部に入力し、当該第１更新部において、前記残響パラメータ推定値および雑音パラメータ推定値の少なくとも一部の更新処理、あるいは前記信号源パラメータ推定値の更新処理、のいずれか一方を、前記パラメータ推定値に関する対数尤度関数の値が増加するように実行するステップと、
    (D) 前記ステップ(C)で得られたパラメータ推定値の更新値の少なくとも一部を第２更新部に入力し、当該第２更新部において、残響パラメータ推定値および雑音パラメータ推定値の少なくとも一部の更新処理、あるいは前記信号源パラメータ推定値の更新処理のうち、前記ステップ(C)で実行されなかったものを、前記パラメータ推定値の更新値に関する対数尤度関数の値が増加するように実行するステップと、
    (E) 終了条件判定部において、終了条件が満たされるか否かを判定するステップと、を有し、
    前記終了条件が満たされない場合、前記第１更新部と前記第２更新部の処理が再び実行される、信号強調方法。
13. 請求項１２の信号強調方法であって、
    前記時間領域信号が、Ｍ個のセンサで観測された信号であり、
    前記残響パラメータ推定値が、前記回帰係数を要素にもつＭ行Ｍ列の回帰行列推定値を含み、
    前記雑音パラメータ推定値が、前記雑音のパワースペクトルを対角要素とするＭ行Ｍ列の雑音パワークロススペクトル行列推定値を含み、
    前記パラメータ推定値が、前記残響パラメータ推定値と、前記信号源パラメータ推定値と、前記雑音パラメータ推定値と、Ｍ次元のステアリングベクトル推定値と、を含み、
    前記第１更新部が、
    源信号推定値更新部と、ステアリングベクトル推定値更新部と、信号源パラメータ推定値更新部と、を含み、
    前記ステップ(C)が、
    (C-1) 前記源信号推定値更新部において、前記観測信号と前記パラメータ推定値とが入力され、雑音重畳信号推定値と、源信号推定値と、前記源信号推定値の誤差分散とを算出するステップと、
    (C-2) 前記ステアリングベクトル推定値更新部において、前記雑音重畳信号推定値と前記源信号推定値とが入力され、ステアリングベクトル推定値の更新値を算出するステップと、
    (C-3) 前記信号源パラメータ推定値更新部において、前記源信号推定値のパワーと前記誤差分散とを加算してパワースペクトルを算出し、前記パワースペクトルを用いて信号源パラメータ推定値の更新値を算出するステップと、を含み、
    前記第２更新部が、源信号パワースペクトル推定値更新部と、雑音パラメータ推定値更新部と、残響パラメータ推定値更新部とを含み、
    前記ステップ(D)が、
    (D-1) 前記信号源パラメータ推定値の更新値を前記源信号パワースペクトル推定値更新部に入力し、前記源信号パワースペクトル推定値更新部において、前記信号源パラメータ推定値の更新値に対応する源信号パワースペクトル推定値の更新値を算出するステップと、
    (D-2) 前記源信号推定値と、前記雑音重畳信号推定値と、前記ステアリングベクトル推定値の更新値を前記雑音パラメータ推定値更新部に入力し、前記雑音パラメータ推定値更新部において、前記雑音パラメータ推定値の更新値を生成するステップと、
    (D-3)前記観測信号と、前記ステアリングベクトル推定値の更新値と、前記源信号パワースペクトル推定値の更新値と、前記雑音パラメータ推定値の更新値とを前記残響パラメータ推定値更新部に入力し、前記残響パラメータ推定値更新部において、前記回帰行列推定値の更新値を算出するステップと、を含む、信号強調方法。
14. 請求項１２の信号強調方法であって、
    前記時間領域信号が１個のセンサで観測された信号であり、
    前記残響パラメータ推定値が、前記回帰係数の推定値を含み、
    前記雑音パラメータ推定値が、前記雑音のパワースペクトルの推定値を含み、
    前記パラメータ推定値が、前記信号源パラメータ推定値と、前記残響パラメータ推定値と、前記雑音パラメータ推定値と、を含み、
    前記第１更新部が、
    雑音抑圧処理部と、信号源パラメータ推定値更新部と、を含み、
    前記ステップ(C)が、
    (C-1) 前記観測信号と前記パラメータ推定値とを前記雑音抑圧処理部に入力し、前記雑音抑圧処理部において、所定の観測区間に属する前記観測信号の集合と前記パラメータ推定値との組合せを前提条件とした前記観測区間に属する残響重畳信号の集合の条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）を特定する複素正規分布の平均及び共分散行列を算出するステップと、
    (C-2) 前記残響パラメータ推定値と、前記複素正規分布の平均及び共分散行列とを前記信号源パラメータ推定値更新部に入力し、前記信号源パラメータ推定値更新部において、信号源パラメータ推定値の更新値を算出するステップと、を含み、
    前記残響重畳信号が、前記観測信号から雑音が取り除かれた信号であり、
    前記信号源パラメータ推定値の更新値は、残響パラメータが前記残響パラメータ推定値に固定された条件下で、第１補助関数値を最大化する値であり、
    前記第１補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第２パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第２パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値であり、
    前記第２更新部が、残響パラメータ推定値更新部を含み、
    前記ステップ(D)が、
    前記信号源パラメータ推定値の更新値と、前記複素正規分布の平均及び共分散行列とを前記残響パラメータ推定値更新部に入力し、前記残響パラメータ推定値更新部において、前記残響パラメータ推定値の更新値を算出するステップを含み、
    前記残響パラメータ推定値の更新値は、信号源パラメータが前記信号源パラメータ推定値の更新値に固定された条件下で、第２補助関数値を最大化する値であり、
    前記第２補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値の更新値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第３パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第３パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値である、信号強調方法。
15. 請求項１２の信号強調方法であって、
    前記時間領域信号が、Ｍ個のセンサで観測された信号であり、Ｍが２以上であり、
    前記残響パラメータ推定値が、前記回帰係数を要素にもつＭ行Ｍ列の回帰行列推定値を含み、
    前記雑音パラメータ推定値が、前記雑音のパワースペクトルの推定値を対角要素とする、Ｍ行Ｍ列の雑音パワークロススペクトル行列推定値を含み、
    前記パラメータ推定値が、前記信号源パラメータ推定値と、前記残響パラメータ推定値と、前記雑音パラメータ推定値と、を含み、
    前記第１更新部が、雑音抑圧処理部と、信号源パラメータ推定値更新部と、を含み、
    前記ステップ(C)が、
    (C-1) 前記観測信号と前記パラメータ推定値とを前記雑音抑圧処理部に入力し、前記雑音抑圧処理部において、所定の観測区間に属する前記観測信号の集合と前記パラメータ推定値との組合せを前提条件とした前記観測区間に属する前記残響重畳信号の集合の条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）を特定する複素正規分布の平均及び共分散行列を算出するステップと、
    (C-2) 前記残響パラメータ推定値と、前記複素正規分布の平均及び共分散行列を前記信号源パラメータ推定値更新部に入力し、前記信号源パラメータ推定値更新部ににおいて、信号源パラメータ推定値の更新値を算出するステップと、を含み、
    前記残響重畳信号が、前記観測信号から雑音が取り除かれた信号であり、
    前記信号源パラメータ推定値の更新値は、残響パラメータが前記残響パラメータ推定値に固定された条件下で、第１補助関数値を最大化する値であり、
    前記第１補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第２パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第２パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値であり、
    前記第２更新部が、残響パラメータ推定値更新部を含み、
    前記ステップ(D)が、
    前記信号源パラメータ推定値の更新値と、前記複素正規分布の平均及び共分散行列とを前記残響パラメータ推定値更新部に入力し、前記残響パラメータ推定値更新部において、前記残響パラメータ推定値の更新値を算出するステップを含み、
    前記残響パラメータ推定値の更新値は、信号源パラメータが前記信号源パラメータ推定値の更新値に固定された条件下で、第２補助関数値を最大化する値であり、
    前記第２補助関数値が、前記観測信号の集合と前記残響重畳信号の集合とが与えられたときの、前記残響パラメータの推定値の更新値と、前記信号源パラメータ推定値の更新値と、前記雑音パラメータ推定値とを含む第３パラメータ推定値に関する尤度関数値ｐ（観測信号の集合，残響重畳信号の集合｜第３パラメータ推定値）の対数関数と、前記条件付事後分布ｐ（残響重畳信号の集合｜観測信号の集合，パラメータ推定値）の積を、当該残響重畳信号の集合について積分した関数の関数値である、信号強調方法。
16. 請求項１２から１５の何れかの信号強調方法の各ステップをコンピュータに実行させるためのプログラム。
17. 請求項１６のプログラムを格納したコンピュータ読み取り可能な記録媒体。

Images