JP5227393B2 - 残響除去装置、残響除去方法、残響除去プログラム、および記録媒体 - Google Patents

Description

この発明は、観測された観測信号から残響信号を取り除く残響除去装置、残響除去方法、残響除去プログラム、および記録媒体に関する。

以下の説明では、音源から発せられた信号を音響信号とし、残響のある部屋で音響信号が発せられ、複数の収音手段（例えばマイクロホン）で収音して得られる信号を観測信号とする。観測信号は、音響信号に残響信号が重畳された信号である。このため、観測信号から本来の音響信号の性質を抽出するのは、困難になると共に、音の明瞭度も低下する。これに対し、残響除去処理は、観測信号から重畳した残響信号を取り除くことで、音響信号本来の性質を抽出しやすくすると共に、音の明瞭度を回復できる。これは、他の様々な音響信号処理システムの要素技術として用いることで、そのシステム全体の性能向上につながる技術である。残響除去処理が要素技術として、性能向上に寄与できるような音響信号処理システムには、例えば以下のようなものがある。
（１）残響信号除去を前処理として用いる音声認識システム
（２）残響信号除去により音声の明瞭度を向上させるＴＶ会議システムなどの通信システム
（３）講演の録音に含まれる残響信号を除去することで、録音された音声の明瞭度を向上させる再生システム
（４）残響信号を除去することで、聞き取りやすさを向上させる補聴器
（５）人が発した声に反応して機械にコマンドを渡す機械制御インターフェース、および機械と人間の対話装置
（６）音響コンテンツの収音時に残響信号を含んで、収音された音響信号の音質を改善するポストプロダクションシステム
（７）音楽コンテンツの残響信号を除去したり、付加したりすることで、音楽コンテンツの音響制御を行う音響エフェクタ
図１に従来の残響除去装置１００の機能構成例を示す（以下、「従来技術１」という。）。残響除去装置１００は、推定部１０４、除去部１０６、音源モデル記憶部１０８により構成される。音源モデル記憶部１０８には、残響信号を含まない音響信号の短時間区間の波形を有限状態機械でモデル化しておくと共に各状態の波形の特徴を信号の自己相関関数で表現した音源モデルを記憶させておく。また、時間領域で観測信号に残響除去フィルタを適用する演算と上記音源モデルに基づき、観測信号から残響除去された信号（理想的な目的信号）の尤もらしさを表現した最適化関数を定義しておく。この最適化関数は、残響除去フィルタ係数と音源モデルの状態時系列をパラメータとして持ち、より適切なフィルタ係数および状態時系列を与えることで、より大きな値をとる関数として設計されている。
以下の説明では、入力される時間領域の観測信号をｘ_ｔ ^（１），．．．，ｘ_ｔ ^（ｑ），．．．，ｘ_ｔ ^（Ｑ）とする。ただし、ｘの右下の添え字「ｔ」は、離散時刻のインデックスを示し、右上のｑ（ｑ＝１，．．．，Ｑ）は収音手段（例えば、「マイクロホン」）のインデックスを示す。以下ではインデックスがｑのマイクロホンを第ｑチャネルのマイクロホンと示す。このことは以下、同様とする。
観測信号ｘ_ｔ ^（ｑ）が入力されると、推定部１０４は、当該観測信号ｘ_ｔ ^（ｑ）と上記最適化関数を用いて、残響除去フィルタを推定する。具体的には、推定部１０４は、上記最適化関数の値を最大にするパラメータを求めることで、残響除去フィルタを推定する。除去部１０６は、推定した残響除去フィルタに観測信号を畳み込むことで、観測信号から残響信号を除去した信号を出力する。この信号を目的信号と呼ぶ。
図２に従来の残響除去装置２００の機能構成例を示す（以下、「従来技術２」という。）。残響除去装置２００は、観測信号をＵ個の周波数帯域に分割する分割部２０２、当該周波数帯域毎の記憶部２０４_ｕ（ｕ＝０，．．．，Ｕ−１）、当該周波数帯域毎の除去部２０６_ｕ、統合部２０８により構成される。
分割部２０２は、観測信号をサブバンド分割することで、Ｕ個の周波数帯域毎に分割されたサブバンド信号を求める。当該分割されたサブバンド信号は時間領域の信号である。また、サブバンド分割の際に、ダウンサンプリング（サンプルの間引き）を行う場合がある。以下の説明では、サブバンド信号をｘ’_ｎ，ｕ ^（ｑ）とする。ただし、ｎをダウンサンプリング後のサンプルのインデックスとし、ｕを周波数帯域のインデックス（ｕ＝０，．．．，Ｕ−１）とする。以下では、第ｑチャネルのマイクロホンで収音された観測信号ｘ_ｔ ^（ｑ）のｕ番目の周波数帯域のサブバンド信号ｘ’_ｎ，ｕ ^（ｑ）について説明する。
上述の通り、Ｕ個の周波数帯域毎に、除去部２０６_ｕ（ｕ＝０，．．．，Ｕ−１）、記憶部２０４_ｕは設けられている。記憶部２０４_ｕには、残響除去フィルタが記憶されている。残響除去フィルタは、予め測定しておいた音源から各マイクロホンまでの室内伝達関数を利用して、この室内伝達関数、分割部２０２によるサブバンド分割処理、除去部２０６_ｕによる残響除去処理、統合部２０８による統合処理の各処理を順次適用した際に得られるシステム全体の入出力関数が、なるべく単位インパルス関数になるように、二乗誤差最小基準に基づき、予め残響除去フィルタの係数を定めておく。
除去部２０６_ｕはサブバンド信号ｘ’_ｎ，ｕ ^（ｑ）に残響除去フィルタを畳み込むことで、サブバンド信号から残響信号を除去する。サブバンド信号から残響信号が除去された周波数帯域毎のサブバンド信号を周波数別目的信号ｓ^〜 _ｎ，ｕとする。そして、統合部２０８は、周波数別目的信号ｓ_ｎ，ｕ ^〜（ｕ＝０，．．．，Ｕ−１）を統合して、目的信号ｓ_ｔ ^〜を求める。
残響除去装置１００、２００の詳細は、非特許文献１、２、３に記載されている。
Ｔ．Ｎａｋａｔａｎｉ，Ｂ．Ｈ．Ｊｕａｎｇ，Ｔ．Ｙｏｓｈｉｏｋａ，Ｋ．Ｋｉｎｏｓｈｉｔａ，Ｍ．Ｄｅｌｃｒｏｉｘ，ａｎｄ Ｍ．Ｍｉｙｏｓｈｉ，″Ｓｔｕｄｙ ｏｎ ｓｐｅｅｃｈ ｄｅｒｅｖｅｒｂｅｒａｔｉｏｎ ｗｉｔｈ ａｕｔｏｃｏｒｒｅｌａｔｉｏｎ ｃｏｄｅｂｏｏｋ，″Ｐｒｏｃ．ＩＥＥＥ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ａｃｏｕｓｔｉｃｓ，Ｓｐｅｅｃｈ，ａｎｄ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ（ＩＣＡＳＳＰ−２００７），ｖｏｌ．Ｉ，ｐｐ．１９３−１９６，Ａｐｒｉｌ ２００７． Ｔ．Ｎａｋａｔａｎｉ，Ｂ．Ｈ．Ｊｕａｎｇ，Ｔ．Ｙｏｓｈｉｏｋａ，Ｋ．Ｋｉｎｏｓｈｉｔａ，Ｍ．Ｍｉｙｏｓｈｉ，″Ｉｍｐｏｒｔａｎｃｅ ｏｆ ｅｎｅｒｇｙ ａｎｄ ｓｐｅｃｔｒａｌ ｆｅａｔｕｒｅｓ ｉｎ Ｇａｕｓｓｉａｎ ｓｏｕｒｃｅ ｍｏｄｅｌ ｆｏｒ ｓｐｅｅｃｈ ｄｅｒｅｖｅｒｂｅｒａｔｉｏｎ，″ＷＡＳＰＡＡ−２００７，２００７ Ｎ．Ｄ．Ｇａｕｂｉｔｃｈ，Ｍ．Ｒ．Ｐ．Ｔｈｏｍａｓ，Ｐ．Ａ．Ｎａｙｌｏｒ，"Ｓｕｂｂａｎｄ Ｍｅｔｈｏｄ ｆｏｒ Ｍｕｌｔｉｃｈａｎｎｅｌ Ｌｅａｓｔ Ｓｑｕａｒｅｓ Ｅｑｕａｌｉｚａｔｉｏｎ ｏｆ Ｒｏｏｍ Ｔｒａｎｓｆｅｒ Ｆｕｎｃｔｉｏｎｓ，"Ｐｒｏｃ．ＩＥＥＥ Ｗｏｒｋｓｈｏｐ ｏｎ Ａｐｐｌｉｃａｔｉｏｎｓ ｏｆ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ ｔｏ Ａｕｄｉｏ ａｎｄ Ａｃｏｕｓｔｉｃｓ（ＷＡＳＰＡＡ−２００７），ｐｐ．１４−１７，２００７．

上記従来技術１の残響除去装置１００では、音響信号の時変性特徴を最適な方法で利用するためには、上記最適化関数の値の最大化の計算を行うために非常に大きなサイズの共分散行列を計算する必要があった。このため、最適化関数の値の最大化に膨大な計算時間が必要であった。共分散行列のサイズが大きくなる理由を説明する。従来技術１で扱う観測信号に関する共分散行列Ｈ（ｒ）は以下の式（１）で表される。
なお、以下の説明では、共分散行列Ｈ（ｒ）とは、従来技術１で扱う観測信号に関する共分散行列であることを示す。ここで、１つの音響信号を２つのマイクロホンで収音したとすると、Ｘ_ｔ−１＝［ｘ⁻ _ｔ−１ ^（１），．．．，ｘ⁻ _ｔ−Ｋ ^（１），ｘ⁻ _ｔ−１ ^（２），．．．，ｘ⁻ _ｔ−Ｋ ^（２）］であり、ｘ⁻ _ｔ ^（１）はｘ_ｔ ^（１）の長さＮの短時間フレームからなる列ベクトルｘ⁻ _ｔ ^（１）＝［ｘ_ｔ ^（１），ｘ_ｔ＋１ ^（１），．．．，ｘ_{ｔ＋Ｎ−１} ^（１）］^Ｔであり、ｘ_ｔ ^（１）、ｘ_ｔ ^（２）はそれぞれ第１チャネル、第２チャネルのマイクロホンで収音された観測信号とする。Ｔは行列、ベクトルの転置を示す。Ｋは予測フィルタ（推定する残響除去フィルタ）の長さである。また、ｒ_ｔは音響信号の短時間フレームからなる列ベクトルｓ⁻ _ｔ＝［ｓ_ｔ，ｓ_ｔ＋１，．．．，ｓ_{ｔ＋Ｎ−１}］^Ｔに関する共分散行列ｒ_ｔ＝Ｅ｛ｓ⁻ _ｔｓ⁻ _ｔ ^Ｔ｝であり、Ｅ｛・｝は期待値関数を示す。一般にｒ_ｔは既知ではないので、音源モデル記憶部１０８に記憶されている音源モデルに基づき、推定部１０４が求めた推定値で代替される。
一般に、予測フィルタ長Ｋは理論的には少なくとも、室内インパルス応答長と同じ長さでなければならない。従って、共分散行列Ｈ（ｒ）のサイズが非常に大きくなる。一方、音響信号が定常信号と仮定すると、上記の共分散行列を相関行列と近似できるため、高速フーリエ変換などの高速な計算方法を用いることが出来るが、音声信号などの時変信号にこの仮定を用いると、残響除去の計算精度が悪くなるという問題があった。このように、残響除去装置１００で精度良く残響除去を行うためには膨大な計算時間を必要とし、また高速に残響除去を行うためには、音響信号が時変信号である場合の残響除去の精度が悪くなるという問題があった。
また上記従来技術２の残響除去装置２００では、残響除去フィルタ（室内伝達関数の逆フィルタ）を予め推定しておかなければならず、あらかじめ室内伝達関数を求めておく必要があった。しかも、室内伝達関数の逆フィルタを用いて、残響除去を行う処理方法は、室内伝達関数の誤差に極めて敏感であり、ある程度以上の誤差が室内伝達関数に含まれていると、残響除去処理によりかえって、音響信号の歪みが増大してしまうという問題があった。さらに室内伝達関数は音源の位置や室温の変化に敏感であり、音源の位置や室温が事前に正確に特定できない場合には、正確な室内伝達関数を定めることはできなかった。このように、残響除去装置２００では、予め精度の良い室内伝達関数を用意する必要があり、しかも、１つの条件下で求めた室内伝達関数は極めて限られた条件でしか残響除去に利用することができなかった。
そこで本発明は次のように残響除去を行う。記憶部に、音響信号を確率密度関数として表現する音源モデルが記憶されている。音響信号を収音して得られた観測信号は複数の周波数帯域のそれぞれに対応する周波数別観測信号に変換される。そして、各周波数帯域における音響信号と観測信号と残響除去フィルタとの関係を表す残響モデルと音源モデルに基づき、各周波数別観測信号を用いて各周波数帯域に対応する残響除去フィルタが推定される。各周波数別観測信号に各残響除去フィルタを適用して、各周波数帯域に対応する周波数別目的信号が求められ、各周波数別目的信号が統合される。

図１は、従来技術１の残響除去装置の機能構成例を示したブロック図である。
図２は、従来技術２の残響除去装置の機能構成例を示したブロック図である。
図３は、実施例１の残響除去装置の機能構成例を示したブロック図である。
図４は、実施例１の残響除去装置の主な処理を示したフローチャートである。
図５は、実施例２の残響除去装置の機能構成例を示したブロック図である。
図６は、実施例２の残響除去装置の主な処理を示したフローチャートである。
図７は、実施例３の残響除去装置の機能構成例を示したブロック図である。
図８は、実施例４の残響除去装置の機能構成例を示したブロック図である。
図９は、実験結果を示す図である。
図１０Ａは、単一のマイクロホンを用いて実施例４に基づき残響除去の効果を実証した実験における観測信号のスペクトルグラムを示す図である。
図１０Ｂは、単一のマイクロホンを用いて実施例４に基づき残響除去の効果を実証した実験結果のスペクトルグラムを示す図である。

以下に、発明を実施するための最良の形態を示す。なお、同じ機能を持つ構成部や同じ処理を行う過程には同じ番号を付し、重複説明を省略する。

図３に実施例１の残響除去装置３００のブロック図を示し、図４に残響除去装置３００の主な処理の流れを示す。図３に示すように、実施例１の残響除去装置３００は、Ｕ個の周波数帯域毎に観測信号を分割する分割部３０２、音源モデル記憶部３０４、周波数帯域毎の推定部３０６_ｕ（ｕ＝０，．．．，Ｕ−１）、周波数帯域毎の除去部３０８_ｕ、統合部３１０を備える。
分割部３０２は、観測信号を周波数帯域毎に分割しながら、観測信号のサンプル数を減少させることで、周波数別観測信号を出力する。この実施例１の分割部３０２は、観測信号に対して短時間分析窓を時間シフトしながら適用するとともに、周波数領域に変換することで周波数帯域毎に分割する。
音源モデル記憶部３０４には、周波数別観測信号の特徴を周波数帯域毎に表現した音源モデルが記憶されている。
推定部３０６_ｕは周波数帯域毎に設けられており、推定部３０６_ｕは音源モデルに関連して定義される観測信号の最適化関数に基づき、周波数別観測信号から残響除去フィルタを推定する。
除去部３０８_ｕは周波数帯域毎に設けられており、周波数別観測信号と残響除去フィルタとを用いて、周波数帯域毎に周波数別目的信号を求める。この実施例１の除去部３０８_ｕは、残響除去フィルタを前記周波数別観測信号に畳み込むことで、周波数別目的信号を求める。
統合部３１０は、周波数別目的信号を統合することで、後述する目的信号を出力する。この実施例１の統合部３１０は、周波数別目的信号を統合し、全周波数帯域を１つにまとめた時間領域の信号に変換することで、後述する目的信号を出力する。
まず、音響信号ｓ_ｔと観測信号ｘ_ｔ ^（ｑ）との関係を説明する。音源から各マイクロホンへの室内伝達関数が共通ゼロ点を持たないと仮定し、音源に一番近いマイクロホンをｑ＝１（第１チャネルのマイクロホン）とする。音響信号と観測信号の関係は以下の式（１１）のように表すことができる。また、詳細は、「Ｍ．Ｍｉｙｏｓｈｉ，“Ｅｓｔｉｍａｔｉｎｇ ＡＲ ｐａｒａｍｅｔｅｒ−ｓｅｔｓ ｆｏｒ ｌｉｎｅａｒ−ｒｅｃｕｒｒｅｎｔ ｓｉｇｎａｌｓ ｉｎ ｃｏｎｖｏｌｕｔｉｖｅ ｍｉｘｔｕｒｅｓ，”Ｐｒｏｃ．ＩＣＡ−２００３，ｐｐ．５８５−５８９，２００３．」に記載されている。
ｈ_０ ^（１）は音源からｑ＝１のマイクロホンへの室内インパルス応答の１タップ目の値であり、ｃ_ｔ ^（ｑ）は予測係数と呼び、推定部３０６_ｕで推定される残響除去フィルタの係数であり、τは離散時間のインデックスであり、Ｋは上述のように、予測フィルタ長（従来技術１で推定する残響除去フィルタのサイズ）である。
ここで、音響信号のゲインを無視すると、右辺の第２項ｈ_０ ^（１）ｓ_ｔは音響信号ｓ_ｔを定数倍した信号であるので、この信号を推定すべき音響信号ｓ_ｔと見なすことが出来る。これにより式（１１）は以下の式（１２）のように書き換えることが出来る。
式（１２）では、現在の観測信号ｘ_ｔ ^（ｑ）は、過去の観測信号の時系列ｘ_ｔ−τ ^（ｑ）から予測され、音響信号ｓ_ｔは、予測の残差信号とみなされる。なお、式（１２）の前提として、第１チャネルのマイクロホン（ｑ＝１）が音源に一番近いと仮定したが、この条件を満たしていない場合も、同じ式（１２）を用いて、観測信号と音響信号との関係を表すことができる。すなわち、第１チャネルのマイクロホン（ｑ＝１）以外のマイクロホンの観測信号に十分な遅延を導入することで、仮想的に音源からの音が最初に到達するマイクロホンを第１チャネルのマイクロホン（ｑ＝１）にすることができ、第１チャネルのマイクロホンを音源に最も近いマイクロホンとして扱うことが出来る。このため、例えば、マイクロホンｑに導入する遅延時間をｄ^（ｑ）タップとすると、式（１２）の表現のままで、ｑ＝１以外の予測係数｛ｃ_１ ^（ｑ），ｃ_２ ^（ｑ），．．．，ｃ_Ｋ ^（ｑ）｝の先頭からｄ^（ｑ）タップに固定値０が代入されていると考えることで、観測信号と音響信号との関係を上記式（１２）同様に表すことができる。
観測信号ｘ_ｔ ^（ｑ）が分割部３０２に入力されると、観測信号を周波数帯域毎に分割しながら、観測信号のサンプル数を減少させることで、周波数別観測信号を出力する（ステップＳ２）。そして、実施例１の分割部３０２は、観測信号に対して短時間分析窓を時間シフトながら適用するとともに、周波数領域に変換することで周波数帯域毎に分割する。例えば、分割部３０２は、短時間フーリエ変換を行う。以下では、分割部３０２は短時間フーリエ変換を行うとして、具体的に説明する。
次に上記式（１２）を一般化して、以下の式（１２’）を考える。
ここで、ｄは現在の観測信号を予測する過去の観測信号に遅延を導入する定数である。ｄ＝１としたとき、式（１２）と一致する。一方、ｄ＞１としたときは、上記式（１２’）は観測信号と音響信号の関係を厳密には表現できなくなる。何故なら、現在の時刻ｔから過去にさかのぼって、ｄタップの音響信号に由来する信号は上記式（１２’）の右辺の過去の信号系列中に含まれていなくなるため、その時間区間の音響信号に由来して現在の観測信号に含まれている残響信号は、過去の観測信号の線形結合では表現できないことになる。「その時間区間の音響信号に由来して現在の観測信号に含まれている残響信号」は室内インパルス応答の最初のｄタップに相当する初期反射音に対応する。従って、上記式（１２’）においては残差信号中に音響信号以外にこの初期反射音が含まれていることを想定している。これを明らかにするために残差信号をｓ_ｔ ^〜と記述している。なお、この明細書では、記号Ａ_α ^〜は、記号Ａの真上に記号^〜が付された組み合わせ文字を表す。
＜周波数信号に対する畳み込み演算＞
次に上記式（１２’）の右辺第１項に含まれる時間領域における畳込みに相当する演算を、周波数領域信号に対して計算する方法を説明する。まず、ある音響信号ｘ_ｔにフィルタ長Ｋのある残響除去フィルタｃ_ｔを時間領域で畳み込んで得られる信号をｙ_ｔとする。時刻ｔ０で始まるｙ_ｔから、窓関数により短時間フレームを時間窓で切り出した信号はｚ変換領域で、以下の式（１３）のように表現できる。
Ｗ_Ｎ（ｙ（ｚ）ｚ^ｔ０）＝Ｗ_Ｎ（ｃ（ｚ）・ｘ（ｚ）ｚ^ｔ０） （１３）
ただし、ｙ（ｚ）＝ｃ（ｚ）・ｘ（ｚ）とし、「・」は畳込みを示し、Ｗ_Ｎ（）は時間領域における長さＮの窓関数に相当する関数とした。Ｗ_Ｎ（ｃ（ｚ））はｃ（ｚ）中の−Ｎ＋１次から０次の項を取り出し、窓の形に比例して各係数を変更し、窓の外の項を除外する。ｚ^ｔ０は時刻ｔ０で始まる短時間フレームを窓関数の中に移動する時間シフトオペレータである。
更に、時刻ｔにおけるフィルタ係数ｃ_ｔから長さＭのフレームを切り出すことを、ｃ_ｔ，Ｍ（ｚ）＝Ｗ_Ｍ ^Ｒ（ｃ（ｚ）ｚ^ｔ）と表現し、Ｗ_Ｍ ^Ｒ（）を長さＭの短時間分析窓（方形窓）を表すものとする。そうすると明らかにｃ（ｚ）＝Σ_τｃ_τＭ，Ｍ（ｚ）となる。上記式（１３）は以下のように書き換えることが出来る。
ここで、式（１４）中のΣ_τｃ_τＭ，Ｍ（ｚ）ｚ^−τＭは、ｃ（ｚ）（式（１３）参照）に対応するものであり、式（１６）中のｘ_{ｔ０−Ｍ＋１−τＭ，Ｍ＋Ｎ−１}（ｚ）はｘ（ｚ）（式（１３）参照）に対応するものである。
また、Ｋ_Ｒ＝＜Ｋ／Ｍ＞とし、＜Ｋ／Ｍ＞はＫ／Ｍ以上の最小の整数を表す。Ｋ_Ｒは推定部３０６_ｕで推定される残響除去フィルタのフィルタ長（タップ数）となる。式（１５）において、窓関数の引数に含まれる項のうち、窓の外の項を除去することで式（１６）が導かれる。
式（１６）中の「ｃ_τＭ，Ｍ（ｚ）ｘ_{ｔ０−Ｍ＋１−τＭ，Ｍ＋Ｎ−１}（ｚ）」は、時間領域におけるフィルタ係数ｃ_τのτＭタップ目から長さＭのフレームを切り出したものと、時間領域の観測信号ｘ_ｔの時刻ｔ０−Ｍ＋１−τＭ時点から長さＭのフレームを切り出したもの、それぞれをｚ領域で乗算したものである。ｚ領域における乗算は畳込み演算に相当するので、上記観測信号ｘ_ｔとフィルタ係数ｃ_ｔの各フレームの時間領域における畳込み演算を表現していることになる。またｃ_τＭ，Ｍ（ｚ）のフレーム長はＭ、ｘ_{ｔ０−Ｍ＋１−τＭ，Ｍ＋Ｎ−１}（ｚ）のフレーム長はＭ＋Ｎ−１であるので、短時間フーリエ変換のポイント数（周波数帯域の数）ＵがＵ≧２Ｍ＋Ｎ−２の時、時間領域の畳込みは短時間フーリエ変換領域の積で厳密に表現される。ここで、音響信号処理でしばしば用いられる近似を利用する。すなわち、短時間分析窓に含まれる信号とフィルタの畳込みは、当該フィルタの長さＭが短時間分析窓長Ｎに比して十分に短い時は、短時間フーリエ変換領域において、その信号とフィルタの積で近似できる。この近似を利用すると、式（１６）はｚ領域の単位円上（短時間フーリエ変換領域に相当）で以下の式（１７）のように書き換えることが出来る。
短時間フーリエ変換表現を用いると、式（１７）は以下の式（１８）のようになる。
ここで、ｎとτは短時間フレームのインデックスであり、Ｙ_ｎ、Ｃ_ｎ、Ｘ_ｎはそれぞれ、ｙ（ｚ）、ｃ（ｚ）、ｘ（ｚ）に対応する時間領域信号から時間窓で切り出した信号の短時間フーリエ変換後の各周波数帯域の値を要素に持つベクトルであり、ｄｉａｇ（Ｘ）はベクトルＸの成分を対角成分に持つ対角行列である。なお、本明細書では、短時間フーリエ変換を以下で表す。ここで、ｔ_τは、フレームτの最初のサンプルの離散時間インデックスを表す。
式（１８）より、時間領域の畳込み演算は、周波数別観測信号の周波数帯域毎の畳込み演算として計算できる。式（１７）においては、Ｍはフレームシフトに相当する値なので、この近似計算においては、フレームシフトＭは窓関数Ｗ_Ｎ（）の窓長Ｎに比して十分に小さい値であることが必要である。
以上で＜補足説明：周波数信号に対する畳み込み演算＞を終わる。
上記式（１２’）の両辺に式（１６）を用いて例えば短時間フーリエ変換を施すことで、以下の式（２２）を得る。
式（２２）は式（２２ａ）に等価である。
ここで、Ｄは式（２２）の遅延ｄに対応し、周波数信号における過去の観測信号に導入する遅延をフレーム数で表したものである。隣接するフレームの周波数信号は、時間領域において相互に重なりを持つ。このため、フレームｎの観測信号（式（２２）の左辺Ｘ_ｎ ^（１））に含まれる音響信号の一部は、その直前のフレームに対応する観測信号にも含まれていることになる。したがって、式（２２）において、直前のフレームを含む過去の観測信号を用いてＸ_ｎ ^（１）を予測すると、音響信号の一部をも予測できてしまう。観測信号の予測可能な部分は、残差信号には含まれないため、音響信号の一部は残響除去によって除去されることになる。これを防ぐために、周波数信号を用いる本発明では、式（２２）のように、現在の観測信号を予測する際に直前のフレームの観測信号を用いず、ある程度の遅延Ｄ以上はなれた過去の観測信号だけを用いるようにする。なおｄ＝ＤＭを満たす時、上記式（１２’）と式（２２）は一致する。以下の説明では、式（２２）を観測信号と音響信号の関係を表現する式として本実施形態を説明する。式（２２）において、Ｘ_ｎ ^（ｑ）は第ｑチャネルのマイクロホンで収音された時間領域信号に関する短時間フーリエ変換に相当する。短時間フーリエ変換は、式（１９）、（２０）に従う。ここで、ｎはフレーム番号を表す。また、周波数帯域ｕ（ｕ＝０，．．．，Ｕ−１）での周波数別観測信号はＸ_ｎ，ｕ ^（ｑ）と表される。このＸ_ｎ，ｕ ^（ｑ）を求めるために、分割部３０２は、短時間分析窓をＭサンプルずつ時間シフトしながら適用するとともに、周波数領域に変換する。これにより、周波数帯域毎に分割した周波数別観測信号Ｘ_ｎ，ｕ ^（ｑ）が得られる。
また、後ほど詳細を説明する推定部３０６_ｕは、周波数別観測信号Ｘ_ｎ，ｕ ^（ｑ）から残響除去のための残響除去フィルタを推定する。当該残響除去フィルタの係数である予測係数Ｃ_τ ^（ｑ）が得られると、目的信号（初期反射音を含む音響信号）Ｓ^〜 _ｎを以下のように推定できる。
式（２３）は周波数帯域ごとにＳ_ｎ ^〜＝［Ｓ_ｎ，０ ^〜，Ｓ_ｎ，１ ^〜，．．．，Ｓ_{ｎ，Ｕ−１} ^〜］の各要素について以下の式（２４）のようにも表すことができる。
ここで式（２５）−（２８）を用いると、式（２４）は式（２９）のようにも表すことができる。
ただしＴはベクトルおよび行列の転置を示す。この実施例では、Ｃ_ｕをｕ番目の周波数帯域の残響除去フィルタという。なお、式（２９）中のＢ_{ｎ−Ｄ，ｕ}Ｃ_ｕ ^Ｔは、各チャネル毎にＢ_ｎ，ｕ ^（ｑ）とＣ_ｕ ^（ｑ）を畳み込んで得られる信号を、全てのｑに関して加算した信号に相当する。推定部３０６_ｕで残響除去フィルタＣ_ｕが推定され、除去部３０８_ｕが式（２９）に基づき、残響信号を除去する。
また、０_Ｄ−１を全ての要素が０のＤ−１次元の行ベクトルとすると、残響除去フィルタＷ_ｕは以下のようにも定義できる。
Ｗ_ｕ＝［１，０_Ｄ−１，Ｃ_ｕ ^（１），０，０_Ｄ−１，Ｃ_ｕ ^（２），．．．，０，０_Ｄ−１，Ｃ_ｕ ^（Ｑ）］この場合、除去部３０８_ｕにより以下の式に基づき残響信号を除去する。
以上のように、推定部３０６_ｕが残響除去フィルタＣ_ｕまたはＷ_ｕを推定できると、除去部３０８_ｕは式（２９）または式（３０）に基づき、残響信号を除去できる。次に、残響除去フィルタの推定の説明の前に、音源モデルについて説明する。
音源モデル記憶部３０４には、周波数別観測信号の特徴を周波数帯域毎に表現した音源モデルが記憶されている。
この実施例の音源モデルは、音響信号の取りうる値の傾向を確率分布で表現する。そして、この確率分布に基づき最適化関数を定義する。音源モデルは例えば、時変正規分布が有効であり、求める周波数別信号Ｓ_ｎ ^〜の確率密度関数は以下のように定義される。
ｐ（Ｓ_ｎ ^〜）＝Ｎ（Ｓ_ｎ ^〜；０，Ψ_ｎ） （３１）
Ψ_ｎ∈Ω_Ψ （３２）
ここで、Ｎ（Ｓ_ｎ ^〜；０，Ψ_ｎ）は平均０、音源モデルの共分散行列Ψ_ｎ＝Ｅ（Ｓ_ｎ ^〜（Ｓ_ｎ ^〜）^＊Ｔ）の多次元複素正規分布を表し、Ψ_ｎは短時間フレームｎごとに異なる値、もしくは同じ値をとる。以下の説明では、Ψ_ｎをモデル共分散行列といい、モデル共分散行列Ψ_ｎは、短時間フレームｎごとに異なる値をとる対角行列であると仮定する。また、「^＊」は、複素共役を表す。Ω_ΨはΨ_ｎが取りうる値全てを含む集合（すなわち、Ψ_ｎのパラメータ空間）を表す。ψ_ｎ，ｕ ^２＝Ｅ（Ｓ_ｎ，ｕ ^〜Ｓ_ｎ，ｕ ^〜＊Ｔ）をΨ_ｎのｕ番目の対角要素を表すものとすると、Ψ_ｎは対角行列なので、確率密度関数は、各周波数帯域ごとに独立に
ｐ（Ｓ_ｎ，ｕ ^〜）＝Ｎ（Ｓ_ｎ，ｕ ^〜；０，ψ_ｎ，ｕ ^２） （３３）
とできる。
周波数帯域毎の推定部３０６_ｕは、音源モデルに関連して定義される観測信号の最適化関数に基づき、周波数別観測信号から残響除去フィルタを推定する（ステップＳ４）。残響除去フィルタの推定の詳細を具体的に説明する。
残響除去フィルタＣ_ｕは、上記式（２５）に示すように、全てのマイクロホンに関する観測信号の予測係数Ｃ_ｕ ^（ｑ）からなるベクトルで表される。予測係数Ｃ_ｕ ^（ｑ）は周波数領域の予測係数である。ψ_ｕ ^２はモデル共分散行列の第ｕ対角要素の時系列を表し、ψ_ｕ ^２＝｛ψ_ｎ，ｕ ^２｝と示す。また、θ_ｕ＝｛Ｃ_ｕ、ψ_ｕ ^２｝を推定パラメータの集合を表すものとする。更に、全周波数帯域の推定パラメータ全体の集合をθ＝｛θ_０，θ_１，．．．，θ_Ｕ−１｝と表す。そして、各周波数帯域ごとの最適化関数として対数尤度関数Ｌ_ｕ（θ_ｕ）および全周波数帯域にわたる最適化関数として対数尤度関数Ｌ（θ）を以下のように定義する。
式（３４）は、式（２９）（３３）に基づき、以下の式（３６）のように表すことができる。
式（３５）の左辺を最大化するパラメータを推定することで、残響除去フィルタの予測係数Ｃ_ｕ ^（ｑ）を求めることが出来る。式（３５）の最大化は、以下の最適アルゴリズムにより実現できる。
１．全ての周波数帯域ｕに関して初期値を例えば以下の式（３７）のように定める。
２．以下の２つの式を収束するまで繰り返す。
２−１．全ての周波数帯域ｕに関して、Ｃ_ｎ，ｕ ^（ｑ）を固定して、最適化関数Ｌ（θ）を最大化するように、モデル共分散行列Ψ_ｎを更新する。
２−２．Ψ_ｎを固定して、全ての周波数帯域ｕに関して、最適化関数Ｌ_ｕ（θ_ｕ）を最大化するように、残響除去フィルタＣ_ｕを更新する。
ただし、上記アルゴリズムの表記において、パラメータＡの値をＢに更新する操作を「Ａ→Ｂ」と記述した。また、「＋」はムーアペンローズの擬似逆行列を表す。なお、上記アルゴリズム中で計算する必要がある観測信号に関する共分散行列Ｈ’（ψ_ｎ，ｕ ^２）は以下の式（４０）のようになる。
この最適化アルゴリズムに基づき、最終的に得られたＣ_ｕを元に残響除去フィルタを構成する。除去部３０８_ｕは、式（２９）または式（３０）に基づき、当該残響除去フィルタＣ_ｕまたはＷ_ｕを周波数別観測信号Ｘ_ｎ，ｕ ^（ｑ）に畳み込むことで、Ｘ_ｎ，ｕ ^（ｑ）から残響信号を除去して、周波数別目的信号Ｓ_ｎ，ｕ ^〜を求める（ステップＳ１２）。
そして、統合部３１０が周波数帯域毎の周波数別目的信号Ｓ_ｎ，ｕ ^〜を統合すると共に、時間領域に変換することで目的信号ｓ_ｔ ^〜を出力する（ステップＳ１４）。具体的には、短時間フーリエ変換のフレームの時系列を時間領域信号に変換する一般的な方法を用いることが出来る。すなわち、各フレームｎごとにＳ_ｎ ^〜＝［Ｓ_ｎ，０ ^〜，Ｓ_ｎ，１ ^〜，．．．，Ｓ_{ｎ，Ｕ−１} ^〜］に短時間逆フーリエ変換を適用して各フレームの時間信号を得ると共に、各フレームの信号をオーバラップ加算することで目的信号ｓ_ｔ ^〜を得る。フレームτの短時間逆フーリエ変換は式（４０ａ）で表される。オーバラップ加算は、短時間逆フーリエ変換を適用して得られる各フレームの時間信号に何らかの時間窓を適用するとともに、分割部で用いたのと同じフレームシフト幅Ｍで信号を加算することで実現される。具体的な計算式は式（４０ｂ）で表される。ここで、ｗ_ｔ ^Ｉは長さＮの時間窓、ｆｌｏｏｒ（ａ）はａ以下の最大の整数を表す。
この実施例１の残響除去装置３００の効果を説明する。この残響除去装置３００による観測信号ｘ_ｔ ^（ｑ）（ｑ＝１，．．．，Ｑ）から残響除去処理を各周波数帯域ごとの演算として近似計算できる。Ｍサンプルずつ時間シフトさせながら長さＮの短時間分析窓を適用して周波数領域信号への変換を行うことで、各周波数帯域毎の残響除去フィルタの長さを短くすることが出来る。そして、残響除去フィルタの推定に必要な共分散行列のサイズを小さくできる。その理由を説明すると、一般的に、残響除去フィルタのサイズと、当該残響除去フィルタを求めるために用いる共分散行列のサイズは等しい。そして、Ｍサンプルずつ時間シフトさせながらＮサンプル分切り取って（長さＮの短時間分析窓を適用して）、周波数領域変換処理を行っているので、従来技術１と比較して畳み込まれる残響除去フィルタのサイズも小さくなる。従って、共分散行列のサイズも小さくなる。このことは、式（１）、式（４０）からも明らかである。つまり、式（１）に示す共分散行列Ｈ（ｒ）のサイズと、式（４０）に示す共分散行列Ｈ’（ψ_ｎ，ｕ ^２）のサイズを比較すると、従来技術１の共分散行列Ｈ（ｒ）のサイズは予測フィルタ長（室内インパルス応答長）Ｋに依存する。しかし、本実施例１で用いた共分散行列Ｈ’（ψ_ｎ，ｕ ^２）は、Ｋ_Ｒ（つまり、＜Ｋ／Ｍ＞）に依存する。何故なら、式（３５）に示すように、共分散行列Ｈ’（ψ_ｎ，ｕ ^２）を構成するＢ_{ｎ−Ｄ，ｕ} ^（ｑ）の要素の数（タップ数）は、Ｋ_Ｒ−Ｄ個だからである。従って、従来技術１と比べると、本実施例１で用いる共分散行列のサイズが小さくできることが理解できよう。残響除去フィルタの推定では共分散行列の計算に加えて、その逆行列の計算が必要であり、これらにかかる計算コストは、残響除去処理全体の計算コストの大部分を占める。更に、この両方の計算コストは、共分散行列のサイズを小さくすることで縮小できる。以上のようにして本実施例では、残響除去処理全体の計算コストを大幅に削減できる。

実施例１では、各周波数帯域ごとに推定された残響除去フィルタを観測信号に畳み込むことで、残響除去を実現した。一方、残響信号を推定し、観測信号のエネルギーと残響信号のエネルギーとの差である差信号を求める方が、実施例１の残響除去方法より残響除去フィルタの推定誤差の影響を受けにくい残響除去を行えることが知られている。例えば、「Ｋ．Ｋｉｎｏｓｈｉｔａ，Ｔ．Ｎａｋａｔａｎｉ，ａｎｄ Ｍ．Ｍｉｙｏｓｈｉ，“Ｓｐｅｃｔｒａｌ ｓｕｂｔｒａｃｔｉｏｎ ｓｔｅｅｒｅｄ ｂｙ ｍｕｌｔｉ−ｓｔｅｐ ｆｏｒｗａｒｄ ｌｉｎｅａｒ ｐｒｅｄｉｃｔｉｏｎ ｆｏｒ ｓｉｎｇｌｅ ｃｈａｎｎｅｌ ｓｐｅｅｃｈ ｄｅｒｅｖｅｒｂｅｒａｔｉｏｎ，”Ｐｒｏｃ．ＩＣＡＳＳＰ−２００６，ｖｏｌ．１，ｐｐ．８１７−８２０，Ｍａｙ，２００６．」に記載されている。この実施例２では、この考え方を適用する。
実施例２の残響除去装置４００について説明する。図５に残響除去装置４００の機能構成例を示し、図６に主な処理の流れを示す。残響除去装置４００は、残響除去装置３００と比較して、除去部３０８_ｕが除去部４０７_ｕに代替されている点で異なる。除去部４０７_ｕは、周波数帯域毎の残響信号生成手段４０８_ｕ、周波数帯域毎の残響信号周波数別パワー生成手段４１０_ｕ、周波数帯域毎の観測信号周波数別パワー生成手段４１２_ｕ、周波数帯域毎の減算手段４１４_ｕ、により構成される。
分割部３０２により観測信号が周波数帯域毎に分割され（ステップＳ２）、推定部３０６_ｕにより、周波数帯域毎の残響除去フィルタが推定されると（ステップＳ４）、残響信号生成手段４０８_ｕは、残響除去フィルタと周波数別観測信号Ｘ_ｎ，ｕ ^（ｑ）を用いて、周波数別残響信号Ｒ_ｎ，ｕを生成する（ステップＳ２２）。具体的には、例えば以下の式（４１）により周波数別残響信号Ｒ_ｎ，ｕを求める。
残響信号周波数別パワー生成手段４１０_ｕは、周波数別残響信号Ｒ_ｎ，ｕの周波数別パワー｜Ｒ_ｎ，ｕ｜^２を求める（ステップＳ２４）。一方、観測信号周波数別パワー生成手段４１２_ｕが例えば、第１チャネルのマイクロホンで収音された周波数別観測信号の周波数別パワー｜Ｘ^（１） _ｎ，ｕ｜^２を求める（ステップＳ２６）。そして、減算手段４１４_ｕが、周波数別残響信号の周波数別パワーと周波数別観測信号の周波数別パワーの差を計算することで差信号｜Ｘ^（１） _ｎ，ｕ｜^２−｜Ｒ_ｎ，ｕ｜^２を求め、当該差信号の計算に用いた周波数別観測信号Ｘ^（１） _ｎ，ｕと当該差信号に基づき、周波数別目的信号を求める（ステップＳ２８）。例えば以下の式に基づき周波数別目的信号Ｓ_ｎ，ｕ ^〜を求める。
ただし、ｍａｘ｛Ａ，Ｂ｝は、Ａ、Ｂのうち大きいほうを選択する関数とし、Ｇ_０は、Ｇ_０＞０であり、パワー減算で信号のエネルギーを抑圧する下限を定めるフロアリング定数とする。そして、統合部４１６が当該周波数別目的信号を時間領域に変換することで、目的信号ｓ_ｔ ^〜を求める（ステップＳ３０）。
この残響除去装置４００は、実施例１の残響除去装置３００より残響除去フィルタに推定誤差が含まれていても音質の劣化の少ない残響除去を行うことが出来る。
また、従来技術の残響除去処理は、時間領域でしか動作させることが出来なかった。しかし、実施例１、２で説明した残響除去装置３００、４００は、周波数領域で動作させるので、ブラインド音源分離やウィーナフィルタなど、周波数領域で動作する他の多くの有用な音声強調技術と組み合わせることが出来る。

図７に実施例３の残響除去装置５００の機能構成例を示す。実施例１の残響除去装置３００と違う主な点を説明すると、（１）残響除去装置３００の分割部３０２は、時間領域の観測信号を時間シフトしながら周波数領域に変換することで、周波数帯域に分割していたのに対し、残響除去装置５００の分割部５０２は、サブバンド分割して周波数帯域に分割する点、（２）残響除去装置３００の除去部、統合部の処理は周波数領域で行っていたのに対し、本実施例の残響除去装置５００の除去部、統合部の処理は時間領域で行っている点で異なる。
サブバンド分割した信号をサブバンド信号とし、サブバンドの数をＶとし、サブバンドのインデックスをｖ（ｖ＝０，．．．，Ｖ−１）とする。推定部５０６_ｖは各サブバンド信号ごとに残響除去フィルタを推定し、除去部５０８_ｖは各サブバンド信号ごとに残響を除去する。統合部５１０により統合されることで目的信号ｓ_ｔ ^〜を求める。分割部５０２によるサブバンド分割処理、統合部５１０による統合処理は、「Ｍ．Ｒ．Ｐｏｒｔｎｏｆｆ，“Ｉｍｐｌｅｍｅｎｔａｔｉｏｎ ｏｆ ｔｈｅ ｄｉｇｉｔａｌ ｐｈａｓｅ ｖｏｃｏｄｅｒ ｕｓｉｎｇ ｔｈｅ ｆａｓｔ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ，”ＩＥＥＥ Ｔｒａｎｓ．ＡＳＳＰ，ｖｏｌ．２４，Ｎｏ．３，ｐｐ．２４３−２４８，１９７６．（以下、「非特許文献Ａ」という。）」や「Ｊ．Ｐ．Ｒｅｉｌｌｙ，Ｍ．Ｗｉｌｂｕｒ，Ｍ．Ｓｅｉｂｅｒｔ，ａｎｄ Ｎ．Ａｈｍａｄｖａｎｄ，“Ｔｈｅ ｃｏｍｐｌｅｘ ｓｕｂｂａｎｄ ｄｅｃｏｍｐｏｓｉｔｉｏｎ ａｎｄ ｉｔｓ ａｐｐｌｉｃａｔｉｏｎ ｔｏ ｔｈｅ ｄｅｃｉｍａｔｉｏｎ ｏｆ ｌａｒｇｅ ａｄａｐｔｉｖｅ ｆｉｌｔｅｒｉｎｇ ｐｒｏｂｌｅｍｓ，”ＩＥＥＥ Ｔｒａｎｓ．Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５０，ｎｏ．１１，ｐｐ．２７３０−２７４３，Ｎｏｖ．２００２」などに記載されている。以下の説明では、非特許文献Ａの技術を用いて説明する。当該非特許文献Ａには、後述する式（５０）が記載されている。また、主な処理の流れは、図４と同様なので、省略する。
まず、音響信号と観測信号の関係を説明する。分割部５０２は、観測信号にサブバンド分割を行い、Ｖ個の周波数帯域毎（サブバンド）に分割する。この分割を非特許文献Ａの定義に従い、式で表すと以下の式（５０）のようになる。
ここで、各サブバンドにおいて、観測信号の周波数シフトおよび低域通過フィルタを適用して得られる信号のサンプルインデックスをｔ（サブバンド処理される前の観測信号の離散時刻と同じ）とし、第ｑチャネルのマイクロホンで収音された観測信号に関するｖ（ｖ＝０，．．．，Ｖ−１）番目のサブバンドのｔ番目のサンプルをｘ_ｔ，ｖ ^（ｑ）とする。ｅ^{−ｊ２πｖτ／Ｖ}はｖ番目のサブバンドに対応する周波数シフト演算子であり、ｈ_τは長さ２Ｎ_ｈ＋１の低域通過フィルタの係数である。そして、式（５０）を上記式（１２’）の両辺に適応すると以下の式を得る。
ここで式（５１）の右辺のｓ_ｔ，ｖ ^〜は初期反射音を含む音響信号にサブバンド分割処理を適用して得られる信号である。本実施例ではｓ_ｔ，ｖ ^〜を求めるべき目的信号として扱う。そして、分割部５０２は、サブバンド分割を行うと共に各サブバンド信号に対してダウンサンプリングを行う。例えば第１チャネルのマイクロホンで収音された観測信号ｘ_ｔ，ｖ ^（１）および音響信号ｓ_ｔ，ｖの各時系列をγ個のサンプル間隔でダウンサンプリング（サンプルの間引き）を行った信号のサンプルのインデックスをｂとし、ダウンサンプリング後に得られるサブバンド信号をｘ_ｂ，ｖ’^（ｑ）やｓ_ｂ，ｖ ^〜’と示す。ダウンサンプリングされた信号のサンプルインデックスｂに対応する、ダウンサンプリングする前の信号のサンプルインデックスをｔ_ｂとする。そうすると、以下の式（５２）のように表すことができる。
一方、ｈ_τは低域通過フィルタなのでこの低域通過フィルタのカットオフ周波数の２倍以上のサンプリング周波数でダウンサンプリングが行われる場合は、アップサンプリングにより高精度にダウンサンプリングする前の信号に復元できる。このアップサンプリングは、例えば以下の手順で行われる。
手順１．ダウンサンプリングされた信号の各サンプル間に、γ−１個の「０」を挿入する。
手順２．低域通過フィルタを適用する。
手順２．では有限長インパルス応答フィルタを用いることが一般的である。これはアップサンプリングにより復元される信号は、ダウンサンプリングされた信号の線形結合で表現できることを意味する。
この関係を用いると式（５２）の右辺の記載ｘ_{ｔｂ−τ，ｖ} ^（ｑ）は以下の式（５３）のように表現できる。
β_τ，ｋはアップサンプリングにおける低域通過フィルタの係数に対応して決まる係数、ｋ_０はアップサンプリングに用いる低域通過フィルタのフィルタリングの遅延、ｋ_０＋ｋ_１＋１はアップサンプリングに用いる低域通過フィルタのフィルタ長に相当する。式（５３）を式（５２）に代入して整理すると、以下の式（５４）を得る。
ここで、α_ｋ，ｖ ^（ｑ）は、式（５３）を式（５２）に代入して整理した時に、ｘ’_{ｂ−ｋ，ｖ} ^（ｑ）の項の係数となるものを表している。ｄ’はα_ｋ，ｖ ^（ｑ）によるフィルタリングの遅延を示し、Ｋ’はα_ｋ，ｖ ^（ｑ）によるフィルタリングのフィルタ長を示す。式（５２）（５３）および間引き間隔γの関係に基づき、ｄ’≒ｄ／γ―ｋ_０、Ｋ’≒Ｋ／γ＋ｋ_１と定めることが出来る。ｄ’≧１の場合、式（５４）は各サブバンド信号に対して、α_ｋ，ｖ ^（ｑ）を予測係数（推定部５０６_ｖで推定される残響除去フィルタの係数）として、過去の観測信号から現在の観測信号を予測した場合の残差信号が初期反射音を含む音響信号となる関係を表している。以下の説明では、式（５４）を各サブバンド信号における観測信号と音響信号との関係を表す式として扱う。
ここで式（５５）−（５８）を定着する。
この場合、式（５４）は、式（５９）のように表現することができる。
この実施例３では、α_ｖをｖ番目のサブバンド信号に対する残響除去フィルタとし、除去部５０８_ｖは上記式（５９）に基づき残響信号の除去を行う。なお、０_ｄ’−１を全ての要素が０のｄ’−１次元の行ベクトルとすると、残響除去フィルタｗ_ｖは以下の式（６０）のようにも表すことができる。
この場合、除去部５０８_ｖは式（６１）に基づき、残響信号の除去を行う。
次に、推定部５０６_ｖによる残響除去フィルタの推定手法について説明する。この実施例の音源モデル記憶部５０４に記憶されている音源モデルは、実施例１、２同様、音響信号の取りうる傾向を確率分布で表現しており、これに基づき最適化関数を定義する。音源モデルとしては、例えば、時変正規分布が有効である。以下の説明では、最も単純な音源モデルとして、各サブバンド間で信号が独立であるモデルを導入する。また、各サブバンド信号は周波数スペクトルが平坦で、信号のエネルギーのみが時間的に変化する時変白色正規過程であると仮定する。
上記式（３１）（３２）同様、パラメータ空間を定義し、以下のように変更する。このとき、ｓ_ｂ ^〜’＝［ｓ_ｂ，０ ^〜’，ｓ_ｂ，１ ^〜’，．．．，ｓ_{ｂ，Ｖ−１} ^〜’］^Ｔの確率密度関数は以下のように定義できる。
ｐ（ｓ_ｂ ^〜’）＝Ｎ（ｓ_ｂ ^〜’；０，Ψ_ｂ’） （３１’）
Ψ_ｂ’∈Ω_Ψ’ （３２’）
ここで、Ｎ（ｓ_ｂ ^〜’；０，Ψ_ｂ’）は平均０、音源モデルの共分散行列Ψ_ｂ’＝Ｅ（ｓ_ｂ ^〜’（ｓ_ｂ ^〜’）^＊Ｔ）の多次元複素正規分布を表し、Ψ_ｂ’はサンプルｂごとに異なる値、もしくは同じ値をとる。以下の説明では、Ψ_ｂ’をモデル共分散行列と呼び、モデル共分散行列Ψ_ｂ’は、サンプルごとに異なる値をとる対角行列であると仮定する。Ω_Ψ’はΨ_ｂ’が取りうる値全てを含む集合（すなわち、Ψ_ｂ’のパラメータ空間）を表す。ψ_ｂ，ｖ’^２＝Ｅ（ｓ_ｂ，ｖ ^〜’（ｓ_ｂ，ｖ ^〜’）^＊）はΨ_ｂ’のｖ番目の対角要素である。Ψ_ｂ’は対角行列なので、確率密度関数は、各サブバンドごとに独立にｐ（ｓ_ｂ，ｖ ^〜’）＝Ｎ（ｓ_ｂ，ｖ ^〜’；０，ψ_ｂ，ｖ’^２）とできる。ψ_ｖ’^２はモデル共分散行列の第ｖ対角要素の時系列を表し、ψ_ｖ’^２＝｛ψ_ｂ，ｖ’^２｝と示す。また、θ_ｖ＝｛α_ｖ，ψ_ｖ’^２｝をサブバンドｖに関する推定パラメータの集合を表すものとする。更に、全サブバンドの推定パラメータ全体の集合をθ’＝｛θ_０，θ_１，．．．，θ_Ｖ−１｝と表す。そして、各サブバンドごとの最適化関数として対数尤度関数Ｌ_ｖ（θ_ｖ）および全サブバンドにわたる最適化関数として対数尤度関数Ｌ’（θ’）を以下のように定義する。
式（６３）は式（５９）、式（３１’）に基づき、式（６４）のように表すことができる。
式（６４）を最大化するパラメータを推定することで、残響除去フィルタの係数の推定値を得ることができる。式（６４）の最大化は、以下の最適化アルゴリズムにより実現できる。
１．全てのサブバンドｖに関して、初期値を以下の式（６５）のように定める。
２．以下の２つの式を収束するまで繰り返す。
２−１．全てのサブバンドｖに関して、α_ｂ，ｖ ^（ｑ）を固定して、最適化関数Ｌ’（θ’）を最大化するように、モデル共分散行列Ψ_ｂ’を更新する。
２−２．Ψ_ｂ’を固定して、全てのサブバンド_ｖに関して、最適化関数Ｌ_ｖ（θ_ｖ）を最大化するように、残響除去フィルタ係数α_ｖを更新する。
最終的に得られたα_ｖをもとに推定部５０６_ｖは残響除去フィルタを構成し、除去部５０８_ｖは上記式当該残響除去フィルタにより上記式（５９）または（６１）に基づいて残響信号を除去することで、周波数別目的信号ｓ_ｂ，ｖ ^〜’を求める。そして、統合部５１０は、周波数別目的信号ｓ_ｂ，ｖ ^〜’をアップサンプリング処理と共に各サブバンド信号を統合することで、目的信号ｓ_ｔ ^〜を求める。
以上説明したように、サブバンド処理では、観測信号を周波数帯域ごとの時間領域信号に分割後にγ個間隔でダウンサンプリングすることで各周波数帯域の時間領域信号のサンプリング周波数を１／γにすることが出来る。
本実施例では、各周波数帯域毎の時間領域信号に対して個別に残響除去処理を行い、これらを統合することで、全周波数帯域にわたる残響除去を実現する。時間領域の信号に対して、ダウンサンプリングする場合としない場合を比較すると、ダウンサンプリングする場合の方が残響除去フィルタの推定に扱う共分散行列のサイズを小さく出来る。何故なら、共分散行列のサイズは、残響除去フィルタのフィルタ長で決まるものであり、残響除去フィルタのフィルタ長Ｋは部屋のインパルス応答のタップ数に対応して決まるものであり、物理的に同じ時間長のインパルス応答はサンプリング周波数が小さくなると少ないタップ数になるためである。換言すれば、γ個間隔でダウンサンプリングを行うことで、残響除去フィルタのフィルタ長はＫ’（＝Ｋ／γ＋ｋ_１）になり、従来技術の残響除去フィルタのフィルタ長Ｋより小さくなる。
残響除去フィルタのフィルタ長が小さくなると、上述したように、残響除去フィルタ推定の際に用いる共分散行列のサイズを小さく出来るので、残響除去フィルタの推定処理の計算コストを削減できる。
また、当該ダウンサンプリングが、低域通過フィルタのカットオフ周波数の２倍以上のサンプリング周波数で行われる場合は、当該ダウンサンプリング処理と共に行ったサブバンド分割処理により求められたサブバンド信号は、アップサンプリングにより高精度に復元できるという性質を有する。従って、統合部５１０による統合処理の際にアップサンプリングをしても、目的信号が劣化することはない。

図８に実施例４の残響除去装置６００の機能構成例を示す。残響除去装置６００は、残響除去装置５００と比較して、除去部５０８_ｖが除去部６０７_ｖに代替されている点で異なる。当該代替により、残響除去装置５００と比較して、残響除去フィルタの推定誤差の影響を受けにくい残響除去を行うことが出来る。理由は、実施例２で説明した通りである。除去部６０７_ｖは、実施例２で説明した除去部４０７_ｕと対応している。除去部６０７_ｖは、周波数帯域毎の残響信号生成手段６０８_ｖ、周波数帯域毎の残響信号周波数別パワー生成手段６１０_ｖ、周波数帯域毎の観測信号周波数別パワー生成手段６１２_ｖ、周波数帯域毎の減算手段６１４_ｖ、により構成される。
残響信号生成手段６０８_ｖは、残響除去フィルタα_ｖと観測信号ｘ_ｔ，ｖ ^（ｑ）を用いて、周波数別残響信号ｒ_ｂ，ｖを求める。具体的には以下の式（７０）により求められる。
ｒ_ｂ，ｖ＝Ｆ_{ｂ−ｄ’，ｖ}・α_ｖ ^Ｔ （７０）
そして、残響信号周波数別パワー生成手段６１０_ｖが、周波数別残響信号の周波数別パワー｜ｒ_ｂ，ｖ｜^２を求める。また、観測信号周波数別パワー生成手段６１２_ｖが、第１チャネルのマイクロホンにより収音された観測信号ｘ_ｂ，ｖ ^（１）の周波数別パワー｜ｘ_ｂ，ｖ ^（１）｜^２を求める。そして、減算手段６１４_ｖが、周波数別残響信号の周波数別パワーと周波数別観測信号の周波数別パワーの差を計算することで差信号｜ｘ_ｂ，ｖ ^（１）｜^２−｜ｒ_ｂ，ｖ｜^２を求め、当該差信号の計算に用いた周波数別観測信号ｘ_ｂ，ｖ ^（１）と当該差信号に基づき、周波数別目的信号を求める（ステップＳ２８）。例えば以下の式に基づき周波数別目的信号ｓ_ｂ，ｖ ^〜’を求める。例えば、以下の式により周波数別目的信号ｓ_ｂ，ｖ ^〜’は求められる。
ただし、ｍａｘ｛Ａ，Ｂ｝は、Ａ、Ｂのうち大きいほうを選択する関数とし、Ｇ_０は、Ｇ_０＞０であり、パワー減算で信号のエネルギーを抑圧する下限を定めるフロアリング定数とする。
そして、それぞれの周波数別目的信号ｓ_ｂ，ｖ’^〜（ｖ＝０，．．．，Ｖ−１）は統合部５１０により統合され、目的信号ｓ_ｔ ^〜として出力される。
残響除去装置６００のような構成にすることで、残響除去装置５００と比較して、残響除去フィルタの推定誤差の影響をあまり受けることなく残響信号の除去を行うことが出来る。

実施例１〜４で説明した残響除去装置３００〜６００では、事前に全ての信号が得られるバッチ処理を前提に構成している。実施例５として、マイクロホンで収音した観測信号を逐次的に残響信号の除去を行うことも可能である。例えば、推定部により推定される残響除去フィルタは、予め定められた時間間隔で（逐次的に）推定、更新されるようにしておく。当該更新時において、その時刻より、前に得られた観測信号の全てもしくは一部に対して、上記最適化アルゴリズムを適用することで、残響除去フィルタを推定する。この推定と共に、残響除去装置３００の推定部３０６_ｕ（図３参照）、残響除去装置４００の残響信号生成手段４０８_ｕ（図５参照）、残響除去装置５００の推定部５０６_ｖ（図７参照）、残響除去装置６００の残響信号生成手段６０８_ｖ（図８参照）は、各時刻に逐次的に得られる観測信号に対して、それまでに得られた最新の残響除去フィルタをその時刻の観測信号に適用するという構成で対応することができる。この逐次処理により、より正確な残響信号の除去を行うことが出来る。
［音源モデルの具体例］
以下に、実施例１から５に関する音源モデルの具体例について、集合Ω_Ψ、Ω_Ψ’の例を示して説明する。主として、実施例１、２、５について説明する。実施例３、４については以下の説明中の各記号について以下の読み替えを行うことで具体例を構成できることから説明を省略する。
Ω_Ψ→Ω_Ψ’
Ψ_ｕ→Ψ_ｖ’
ψ_ｎ，ｕ→ψ_ｂ，ｖ’
Ｘ_ｎ，ｕ ^（ｑ）→ｘ_ｂ，ｖ ^（ｑ）’
Ｓ_ｎ，ｕ ^〜→ｓ_ｂ，ｖ ^〜’
Ｂ_ｎ，ｕ→Ｆ_ｂ，ｖ
Ｄ→ｄ’
Ｃ_ｕ→α_ｖ
ｉ_ｎ→ｉ_ｂ
式（３８）→式（６６）
式（３９）→式（６７）
３０６_ｕ→５０６_ｖ
（１）１つ目の具体例として、集合Ω_Ψが任意の正定値対角行列からなる集合とする。これは、ψ_ｎ，ｕ ^２が任意の正の値をとることが出来ることを意味する。このとき上記最適化アルゴリズムの中で、式（３８）の更新式は、全ての周波数帯域で個別に計算される以下の更新式（８０）に置き換えることが出来る。なお、式（３９）の更新式については変更はない。
（２）２つ目の具体例を説明する。非特許文献１記載の技術と同様に、音響信号の波形を有限状態機械でモデル化する場合について説明する。このとき、集合Ω_Ψは有限個の正定値対角行列からなる集合となる。各行列は、観測信号の短時間信号に対応する周波数領域信号が取りうる有限個の状態のそれぞれに対応する共分散行列になる。これらの有限個の行列は、事前に残響を含まない環境で収音された音響信号の周波数領域信号やその共分散行列をクラスタリングするなどの手法に基づき構成することが出来る。また、有限個の行列の数をＺとし、そのインデックスをｉ（ｉ＝１，．．．，Ｚ）とし、状態ｉに対応する共分散行列をΨ（ｉ）とする。
そうすると、上記繰り返しアルゴリズムの中で推定すべきパラメータは、共分散行列の代わりにインデックスの値となる。以下、時刻ｎの状態をｉ_ｎとし、状態ｉ_ｎに対応する共分散行列をΨ（ｉ_ｎ）とし、共分散行列Ψ（ｉ_ｎ）の対角要素をψ_ｕ ^２（ｉ_ｎ）とする。各時刻における音源モデルの状態ｉ_ｎは、各周波数帯域毎に決まる値ではなく、全周波数帯域に対して１つ決まる値である。このため、対数尤度関数をもとに定められる最適化関数は、全周波数帯域に対して以下の式（８１）のように定義できる。
ここで、推定パラメータθ＝｛Ｃ，Ｉ｝は、ｉ_ｎの時系列Ｉ＝｛ｉ_１，ｉ_２・・・｝と各周波数帯域ごとの予測係数Ｃ＝｛Ｃ_０，Ｃ_１，．．．，Ｃ_Ｕ−１｝から構成されているものとする。この最適化関数に基づき、前記最適化アルゴリズムのうち、式（３８）の更新式は、全周波数帯域に関する以下の更新式（８２）に置き換えることが出来る。なお、式（３９）の更新式については変更はない。
式（３８）から式（８２）への置き換えにより、推定部３０６_ｕはより正確に、残響除去フィルタの推定を行うことが出来る。
（３）３つ目の具体例を説明する。（２）で説明した状態ｉ_ｎを確率変数と仮定することで、より精密な音源モデルに基づく最適化関数を構成することが出来る。一例として、状態ｉ_ｎが一次のマルコフ過程でモデル化できる場合を説明する。マルコフ過程の仮定によりｐ（Ｉ）＝ｐ（ｉ）Π_ｎｐ（ｉ_ｎ｜ｉ_ｎ−１）と出来る。音源モデルのパラメータは、任意の状態ｉ、ｊに対するｐ（ｉ）、ｐ（ｉ｜ｊ）、および各状態における共分散行列Ψ（ｉ）であり、これらのパラメータは残響を含まない環境で収音された音響信号と共に事前に用意できる。このとき残響信号の除去のための最適化関数は、以下の式のようになる。
式（８３）の最適化関数における推定パラメータθは有限状態機械で定義した推定パラメータと同じである。式（８３）の最適化関数は上記最適化アルゴリズムにおいて、式（３８）の状態の更新式のみを以下の更新式で置き換えることで容易に最大化できる。
なお、上記式（８４）の最大化は、公知の技術であるダイナミックプログラミングを用いることで、効率的に計算できる。
実施例１〜５の説明において、観測信号、音響信号の関係を導いた上記式（１２’）では異なるマイクロホン間で室内伝達関数が共通ゼロ点を持たないこと、また、マイクロホンの本数は２本以上必要であることを仮定した。しかし、本発明で構成した実施例１から５に基づく残響除去法ではこれらの仮定が成立していない場合においても、良好な残響除去が実現できるこが実験的に確認されている。
単一のマイクロホンを用いて実施例４に基づき残響除去装置の効果を実証した実験結果について説明する。対象となる音声は、女性一名が発した５単語の発話列で構成される音声信号である。観測信号は残響のある部屋で測定された１チャンネル室内インパルス応答を畳み込むことで合成した。残響時間（ＲＴ６０）は０．５秒である。図１０に観測信号（図１０Ａ）と本実施例を適用して得られた信号（図１０Ｂ）のスペクトルグラムを示す。図には、最初の２単語のみを表示している。図１０より、残響が効果的に抑制されていることが確認できる。
従って、本発明は、マイクロホンの数がＱ＝１の場合やマイクロホン間で室内伝達関数が共通ゼロ点を持つ場合にも適用できる。また、上記従来技術１の場合、音源に最も近いマイクロホンを第１チャネルのマイクロホンとして既知である仮定したが、本発明の技術の場合は、音源に最も近いマイクロホンが既知であるという仮定は必要としないことが実験的に確認されている。
また、実施例１〜５の分割部の処理は、上述では、短時間フーリエ変換、サブバンド分割を用いた。その他の周波数領域に分割する手法として、観測信号のサンプル数を減少させるようにさえすれば、ウェーブレット変換や離散コサイン変換などを用いても良い。また、それらの変換が周波数帯域の間の信号が無相関にならないような変換であっても相関を近似的に無視することで、同様の効果を得ることができる。
また、残響除去フィルタＣ_ｕ、α_ｖ、の最適化のために、上記式（３９）（Ｃ_ｕの推定の場合）、上記式（６７）（α_ｖの推定の場合）を計算する代わりに、適応フィルタでしばしば用いられる逐次推定アルゴリズムを用いることも出来る。そのような最適化手法としては、公知の技術であるＬＭＳ（Ｌｅａｓｔ Ｍｅａｎ Ｓｑｕａｒｅ）法、ＲＬＳ（Ｒｅｃｕｒｓｉｖｅ Ｌｅａｓｔ Ｓｑｕａｒｅｓ）法、最急降下法、共役勾配法、などが知られている。これにより、１回の繰り返しに必要な計算量を大幅に縮小できる。従って、少ない計算コストで実時間内に少なくとも１回以上の繰り返し推定を行うことが出来る。このため、比較的安価なＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）を用いても、実時間処理を実現できる。１回の繰り返しだけでは必ずしも精度の高い残響除去フィルタは得られないが、時間の経過と共に逐次的に推定精度を改善できる。
＜ハードウェア構成＞
本実施例で説明した、プログラムで機能させる残響除去装置は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、入力部、出力部、補助記憶装置、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）及びバスを有している（何れも図示せず）。
ＣＰＵは、読み込まれた各種プログラムに従って様々な演算処理を実行する。補助記憶装置は、例えば、ハードディスク、ＭＯ（Ｍａｇｎｅｔｏ−Ｏｐｔｉｃａｌ ｄｉｓｃ）、半導体メモリ等であり、ＲＡＭは、ＳＲＡＭ（Ｓｔａｔｉｃ Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＤＲＡＭ（Ｄｙｎａｍｉｃ Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）等である。また、バスは、ＣＰＵ、入力部、出力部、補助記憶装置、ＲＡＭ及びＲＯＭを通信可能に接続している。
＜ハードウェアとソフトウェアとの協働＞
本実施例の残響除去装置は、上述のようなハードウェアに所定のプログラムが読み込まれ、ＣＰＵがそれを実行することによって構築される。以下、このように構築される各装置の機能構成を説明する。
残響除去装置の入力部、出力部は、所定のプログラムが読み込まれたＣＰＵの制御のもと駆動するＬＡＮカード、モデム等の通信装置である。分割部、推定部、処理部は、所定のプログラムがＣＰＵに読み込まれ、実行されることによって構築される演算部である。音源モデル記憶部は上記補助記憶装置として機能する。
［実験結果］
本実施例の残響除去装置の効果を実証した実験結果について説明する。この実験では、実施例１で説明した残響除去装置３００と従来技術で説明した残響除去装置１００を比較した。対象となる音声は、５単語の発話列で構成される音声信号であり、男性と女性、各一名が発した合計２種類の発話列からなる。観測信号は残響のある部屋で測定された２チャネル室内インパルス応答を畳み込むことで合成した、残響時間（ＲＴ６０）は０．５秒である。残響除去は各発話列に対して行い、その性能は残響除去後の信号のケプストラム歪み（ｃｅｐｓｔｒｕｍ ｄｉｓｔｏｒｔｉｏｎ、以下、単に「ＣＤ」と示す。）と残響除去処理の実時間性（ｒｅａｌ ｔｉｍｅ ｆａｃｔｏｒ、以下単に「ＲＴＦ」と示す。）を用いて残響除去性能を評価した。ＣＤは以下で定義される。
ここで、ｃ_ｋ＾とｃ_ｋはおのおの評価する音声信号とクリーン音声信号のケプストラム係数で、Ｄ＝１２とした。この評価尺度で、エネルギー時間パターンとスペクトル包絡の両方に関して、信号に含まれる歪みを評価できる。ＲＴＦは（残響除去処理に要した時間）／（観測信号の時間）とした。実験に用いた残響除去法は何れもリナックスコンピュータ上でプログラミング言語マトラブで実装した。標本化周波数は８ｋＨｚ、短時間分析窓長Ｎは２５６とした。
図９にグラフで示した実験結果を示す。縦軸がＣＤを示し、横軸（対数表示）がＲＴＦを示す。残響除去装置３００（実施例１）については、折れ線で示しており、フレームシフトＭの値を２５６、１２８、６４、３２、１６、８の場合についてのＲＴＦ、ＣＤの関係を示す。残響除去装置１００（従来技術１）については、×印を付す。観測信号は破線で示し、ＣＤの値が約４．１である。
図９から残響除去装置１００では、ＲＴＦ９０に対してＣＤが約２．４である。これに対し、残響除去装置３００では例えばＭ＝６４の場合は、ＣＤが従来技術とほぼ等しい約２．４であるにも関わらず、ＲＴＦが約２．５となっている。この結果より、残響除去装置３００は残響除去装置１００よりも優れていることが理解できよう。また、残響除去装置３００では、ＲＴＦが増加するにつれて、ＣＤが減少していることも理解できよう。
発明の効果
本発明によると、観測信号が複数の周波数帯域のそれぞれに対応する周波数別観測信号に変換され、各周波数別観測信号を用いて各周波数帯域に対応する残響除去フィルタが推定される。各周波数帯域に対応する残響除去フィルタの次数は、観測信号をそのまま用いた場合の残響除去フィルタの次数よりも小さい。これに呼応して、共分散行列のサイズが小さくなるため、残響除去フィルタの推定に係る計算コストを低減することができる。また、各周波数別観測信号を用いて残響除去フィルタを推定するから、予め室内伝達関数が既知である必要が無い。

Claims (8)

1. 音源から発せられた音響信号を収音して得られた観測信号に残響除去フィルタを適用することでこの観測信号から残響信号を除去する残響除去装置であって、
   音響信号を確率密度関数として表現する音源モデルを記憶している音源モデル記憶部と、
   上記観測信号を複数の周波数帯域のそれぞれに対応する周波数別観測信号に変換する分割部と、
   現在の観測信号を、所定の遅延を持つ過去の観測信号に残響除去フィルタを適用して得られる信号に音響信号を加算して得られる信号として表現する自己回帰モデルと、上記音源モデルに基づき、各上記周波数別観測信号を用いて各上記周波数帯域に対応する残響除去フィルタを求める推定部と、
   各上記周波数別観測信号に上記推定部によって得られた各上記残響除去フィルタを適用して、各上記周波数帯域に対応する周波数別目的信号を求める除去部と、
   各上記周波数別目的信号を統合する統合部と
   を含む残響除去装置。
2. 請求項１に記載の残響除去装置であって、
   上記音源モデルは、平均０且つ周波数帯域間で相関を持たない時変複素正規分布モデルである
   ことを特徴とする残響除去装置。
3. 請求項２に記載の残響除去装置であって、
   上記推定部は、上記周波数別目的信号の分散を推定し、この推定された周波数別目的信号の分散で正規化された各上記周波数別観測信号の共分散行列を用いて上記残響除去フィルタを推定する
   ことを特徴とする残響除去装置。
4. 音源から発せられた音響信号を収音して得られた観測信号に残響除去フィルタを適用することでこの観測信号から残響信号を除去する残響除去方法であって、
   音源モデル記憶部に音響信号を確率密度関数として表現する音源モデルが記憶されており、
   上記観測信号を複数の周波数帯域のそれぞれに対応する周波数別観測信号に変換する分割ステップと、
   現在の観測信号を、所定の遅延を持つ過去の観測信号に残響除去フィルタを適用して得られる信号に音響信号を加算して得られる信号として表現する自己回帰モデルと、上記音源モデルに基づき、各上記周波数別観測信号を用いて各上記周波数帯域に対応する残響除去フィルタを求める推定ステップと、
   各上記周波数別観測信号に上記推定ステップで得られた各上記残響除去フィルタを適用して、各上記周波数帯域に対応する周波数別目的信号を求める除去ステップと、
   各上記周波数別目的信号を統合する統合ステップと
   を含む残響除去方法。
5. 請求項４に記載の残響除去方法であって、
   上記音源モデルは、平均０且つ周波数帯域間で相関を持たない時変複素正規分布モデルである
   ことを特徴とする残響除去方法。
6. 請求項５に記載の残響除去方法であって、
   上記推定ステップでは、上記周波数別目的信号の分散を推定し、この推定された周波数別目的信号の分散で正規化された各上記周波数別観測信号の共分散行列を用いて上記残響除去フィルタを推定する
   ことを特徴とする残響除去方法。
7. 請求項１から請求項３のいずれかに記載の残響除去装置としてコンピュータを動作させる残響除去プログラム。
8. 請求項１から請求項３のいずれかに記載の残響除去装置としてコンピュータを動作させるプログラムを記録したコンピュータが読み取り可能な記録媒体。