(57)【特許請求の範囲】
【請求項１】
後部残響を伴う音響信号から後部残響を除去する残響除去装置であって、
Ｍ（Ｍ≧１）個のセンサによってそれぞれ観測されたＭ個のチャネルｍ（ｍ＝１，...，Ｍ）の上記音響信号をそれぞれ複数の時点でサンプリングして得られた離散音響信号値を記憶するメモリと、
長時間区間におけるＭ個のチャネルｍの各離散音響信号値を線形結合した線形予測項と、当該長時間区間より所定時間後の離散時間ｎにおけるチャネルｗ（ｗ＝１，...，Ｍ）の予測誤差項と、の和によって、当該長時間区間より所定時間後の離散時間ｎにおけるチャネルｗ（ｗ＝１，...，Ｍ）の離散音響信号値を表現した線形予測モデルである、チャネルｗのマルチステップ線形予測モデルの各線形予測係数を、複数の上記離散音響信号値を用いて算出するモデル適用部と、
上記チャネルｗのマルチステップ線形予測モデルの各線形予測係数と複数の上記離散音響信号値とを上記チャネルｗのマルチステップ線形予測モデルの上記線形予測項に代入して得た線形予測値を、離散時間ｎにおけるチャネルｗの後部残響予測値として出力する後部残響予測部と、を有する。
【請求項２】
請求項１に記載の残響除去装置であって、
上記モデル適用部は、
各離散時間の上記離散音響信号値から、当該離散時間直前の短時間区間内の各離散音響信号値と自己相関性を持つ自己相関成分を抑制し、擬似白色化した離散音響信号値を生成する擬似白色化部と、
上記擬似白色化した離散音響信号値を用い、上記マルチステップ線形予測モデルの各線形予測係数を算出する第１線形予測係数算出部と、を有し、
上記短時間区間は、上記長時間区間よりも短い。
【請求項３】
請求項２に記載の残響除去装置であって、
上記擬似白色化部は、
上記短時間区間におけるチャネルｍの各離散音響信号値を線形結合した線形予測項と、上記短時間区間直後の離散時間ｎにおける当該チャネルｍの予測誤差項と、の和によって、上記短時間区間直後の離散時間ｎにおける当該チャネルｍの離散音響信号値を表現した線形予測モデルである、チャネルｍの短時間線形予測モデルの各線形予測係数を、上記離散音響信号値を用いて算出する第２線形予測係数算出部と、
上記第２線形予測係数算出部で算出された上記各線形予測係数をチャネルｍの上記短時間線形予測モデルに代入して得られる逆フィルタに当該チャネルｍの上記離散音響信号値を代入し、それによって得られる当該短時間線形予測モデルの上記予測誤差項の値を当該チャネルｍの上記擬似白色化した離散音響信号値として出力する逆フィルタ処理部と、を有する。
【請求項４】
請求項３に記載の残響除去装置であって、
Ｍ≧２であり、
上記第２線形予測係数算出部は、
チャネル毎に上記離散音響信号値の自己相関係数を算出する自己相関係数算出部と、
Ｍ個のチャネルに対して求められたＭ個の上記自己相関係数をチャネル間で平均した平均自己相関係数を算出する自己相関係数平均化部と、
上記平均自己相関係数を用い、上記短時間線形予測モデルの各線形予測係数を算出する方程式演算部と、を有する。
【請求項５】
請求項３に記載の残響除去装置であって、
Ｍ≧２であり、
上記第２線形予測係数算出部は、
上記Ｍ個のセンサのうち、音響信号の音源に最も近い１つのセンサで観測された音響信号を複数の時点でサンプリングして得られた離散音響信号値の自己相関係数を算出する自己相関係数算出部と、
上記自己相関係数を用い、上記短時間線形予測モデルの各線形予測係数を算出する方程式演算部と、を有する。
【請求項６】
請求項１に記載の残響除去装置であって、
各チャネルの上記離散音響信号値を周波数領域の離散音響信号値に変換し、各チャネルの上記後部残響予測値を周波数領域の後部残響予測値に変換する周波数領域変換部と、
上記周波数領域の離散音響信号値の振幅スペクトルと、上記周波数領域の後部残響予測値の振幅スペクトルとの相対値をチャネル毎に求め、当該相対値を各チャネルの後部残響除去信号値の振幅スペクトル予測値として出力する後部残響除去部と、を有する。
【請求項７】
請求項６に記載の残響除去装置であって、
上記後部残響予測部は、
上記モデル適用部で算出された上記各線形予測係数と複数の擬似白色化された上記離散音響信号値とを上記線形予測項に代入して得られた線形予測値を、離散時間ｎにおけるチャネルｗの後部残響予測値として算出し、
上記周波数領域変換部は、
擬似白色化された各チャネルの上記離散音響信号値を周波数領域の離散音響信号値に変換する。
【請求項８】
請求項６に記載の残響除去装置であって、
チャネルｗの上記後部残響除去信号値の振幅スペクトル予測値と、チャネルｗの上記周波数領域の離散音響信号値の位相情報とを用い、チャネルｗの後部残響除去信号値の複素スペクトル予測値を算出する複素スペクトル生成部と、
チャネルｗの上記後部残響除去信号値の複素スペクトル予測値を時間領域に変換したチャネルｗの後部残響除去信号推定値を算出する時間領域変換部とを、さらに有する。
【請求項９】
請求項８に記載の残響除去装置であって、
Ｍ≧２であり、
上記モデル適用部は、
複数のチャネルに対してそれぞれ上記各線形予測係数を算出し、
上記後部残響予測部は、
複数のチャネルに対してそれぞれ上記後部残響予測値を算出し、
上記後部残響除去部は、
複数のチャネルに対してそれぞれ上記後部残響除去信号値の振幅スペクトル予測値を算出し、
上記複素スペクトル生成部は、
複数のチャネルに対してそれぞれ上記後部残響除去信号値の複素スペクトル予測値を算出し、
上記時間領域変換部は、
複数のチャネルに対してそれぞれ上記後部残響除去信号推定値を算出し、
当該残響除去装置は、
各チャネルの上記後部残響除去信号推定値をそれぞれ或る遅延量で遅延させた場合に、遅延後の各チャネルの上記後部残響除去信号推定値のチャネル間相互相関が極大となる、各チャネルの当該遅延量を決定する遅延量算出部を有する。
【請求項１０】
請求項９に記載の残響除去装置であって、
各チャネルの離散時間ｎの上記後部残響除去信号推定値を、それぞれのチャネルに対して算出された上記遅延量だけ遅延させる遅延部と、
上記遅延部で遅延させた離散時間ｎの上記後部残響除去信号推定値の和を、補正残響除去信号値として算出する遅延補正部と、を有する。
【請求項１１】
請求項１に記載の残響除去装置であって、
Ｍ≧２である。
【請求項１２】
請求項１に記載の残響除去装置であって、
上記マルチステップ線形予測モデルは、
ｘｗ（ｎ）をチャネルｗ（ｗ＝１，...，Ｍ）に対応する離散時間ｎの離散音響信号値とし、ｘｍ（ｎ）をチャネルｍ（ｍ＝１，...，Ｍ）に対応する離散時間ｎの離散音響信号値とし、ｅｗ（ｎ）をチャネルｗ及び離散時間ｎに対応する予測誤差とし、Ｎを正の整数とし、[・]をガウス記号とし、αｗ，ｍ（ｐ）を、ｘｗ（ｎ）に対応するｐ番目の線形予測係数とし、Ｄをステップサイズを示す定数とした場合における、
【数２４】
0004774100.tif 000026
である。
【請求項１３】
後部残響を伴う音響信号から後部残響を除去する残響除去方法であって、
Ｍ（Ｍ≧１）個のセンサによってそれぞれ観測されたＭ個のチャネルｍ（ｍ＝１，...，Ｍ）の上記音響信号をそれぞれ複数の時点でサンプリングして得られた離散音響信号値をメモリに記憶する離散音響信号記憶過程と、
長時間区間におけるＭ個のチャネルｍの各離散音響信号値を線形結合した線形予測項と、当該長時間区間より所定時間後の離散時間ｎにおけるチャネルｗの予測誤差項と、の和によって、当該長時間区間より所定時間後の離散時間ｎにおけるチャネルｗの離散音響信号値を表現した線形予測モデルである、チャネルｗのマルチステップ線形予測モデルの各線形予測係数を、複数の上記離散音響信号値を用いて算出するモデル適用過程と、
上記チャネルｗのマルチステップ線形予測モデルの各線形予測係数と複数の上記離散音響信号値とを上記チャネルｗのマルチステップ線形予測モデルの上記線形予測項に代入して得た線形予測値を、離散時間ｎにおけるチャネルｗの後部残響予測値として出力する後部残響予測過程と、を有する。
【請求項１４】
請求項１３に記載の残響除去方法であって、
上記モデル適用過程は、
各離散時間の上記離散音響信号値から、当該離散時間直前の短時間区間内の各離散音響信号値と自己相関性を持つ自己相関成分を抑制し、擬似白色化した離散音響信号値を生成する擬似白色化過程と、
上記擬似白色化した離散音響信号値を用い、上記マルチステップ線形予測モデルの各線形予測係数を算出する第１線形予測係数算出過程と、を有し、
上記短時間区間は、上記長時間区間よりも短い。
【請求項１５】
請求項１４に記載の残響除去方法であって、
上記擬似白色化過程は、
上記短時間区間におけるチャネルｍの各離散音響信号値を線形結合した線形予測項と、上記短時間区間直後の離散時間ｎにおける当該チャネルｍの予測誤差項と、の和によって、上記短時間区間直後の離散時間ｎにおける当該チャネルｍの離散音響信号値を表現した線形予測モデルである、チャネルｍの短時間線形予測モデルの各線形予測係数を、上記離散音響信号値を用いて算出する第２線形予測係数算出過程と、
上記第２線形予測係数算出過程で算出された上記各線形予測係数をチャネルｍの上記短時間線形予測モデルに代入して得られる逆フィルタに当該チャネルｍの上記離散音響信号値を代入し、それによって得られる当該短時間線形予測モデルの上記予測誤差項の値を当該チャネルｍの上記擬似白色化した離散音響信号値として出力する逆フィルタ処理過程と、を有する。
【請求項１６】
請求項１５に記載の残響除去方法であって、
Ｍ≧２であり、
上記第２線形予測係数算出過程は、
チャネル毎に上記離散音響信号値の自己相関係数を算出する自己相関係数算出過程と、
Ｍ個のチャネルに対して求められたＭ個の上記自己相関係数をチャネル間で平均した平均自己相関係数を算出する自己相関係数平均化過程と、
上記平均自己相関係数を用い、上記短時間線形予測モデルの各線形予測係数を算出する方程式演算過程と、を有する。
【請求項１７】
請求項１５に記載の残響除去方法であって、
Ｍ≧２であり、
上記第２線形予測係数算出過程は、
上記Ｍ個のセンサのうち、音響信号の音源に最も近い１つのセンサで観測された音響信号を複数の時点でサンプリングして得られた離散音響信号値の自己相関係数を算出する自己相関係数算出過程と、
上記自己相関係数を用い、上記短時間線形予測モデルの各線形予測係数を算出する方程式演算過程と、を有する。
【請求項１８】
請求項１３に記載の残響除去方法であって、
各チャネルの上記離散音響信号値を周波数領域の離散音響信号値に変換し、各チャネルの上記後部残響予測値を周波数領域の後部残響予測値に変換する周波数領域変換過程と、
上記周波数領域の離散音響信号値の振幅スペクトルと、上記周波数領域の後部残響予測値の振幅スペクトルとの相対値をチャネル毎に求め、当該相対値を各チャネルの後部残響除去信号値の振幅スペクトル予測値として出力する後部残響除去過程と、を有する。
【請求項１９】
請求項１８に記載の残響除去方法であって、
上記後部残響予測過程は、
上記モデル適用部で算出された上記各線形予測係数と複数の擬似白色化された上記離散音響信号値とを上記線形予測項に代入して得られた線形予測値を、離散時間ｎにおけるチャネルｗの後部残響予測値として算出する過程であり、
上記周波数領域変換部では、
擬似白色化された各チャネルの上記離散音響信号値を周波数領域の離散音響信号値に変換する。
【請求項２０】
請求項１８に記載の残響除去方法であって、
チャネルｗの上記後部残響除去信号値の振幅スペクトル予測値と、チャネルｗの上記周波数領域の離散音響信号値の位相情報とを用い、チャネルｗの後部残響除去信号値の複素スペクトル予測値を算出する複素スペクトル生成過程と、
チャネルｗの上記後部残響除去信号値の複素スペクトル予測値を時間領域に変換したチャネルｗの後部残響除去信号推定値を算出する時間領域変換過程とを、さらに有する。
【請求項２１】
請求項２０に記載の残響除去方法であって、
Ｍ≧２であり、
上記モデル適用過程は、
複数のチャネルに対してそれぞれ上記各線形予測係数を算出する過程であり、
上記後部残響予測過程は、
複数のチャネルに対してそれぞれ上記後部残響予測値を算出する過程であり、
上記後部残響除去過程は、
複数のチャネルに対してそれぞれ上記後部残響除去信号値の振幅スペクトル予測値を算出する過程であり、
上記複素スペクトル生成過程は、
複数のチャネルに対してそれぞれ上記後部残響除去信号値の複素スペクトル予測値を算出する過程であり、
上記時間領域変換過程は、
複数のチャネルに対してそれぞれ上記後部残響除去信号推定値を算出する過程であり、
当該残響除去方法は、
各チャネルの上記後部残響除去信号推定値をそれぞれ或る遅延量で遅延させた場合に、遅延後の各チャネルの上記後部残響除去信号推定値のチャネル間相互相関が極大となる、各チャネルの当該遅延量を決定する遅延量算出過程を有する。
【請求項２２】
請求項２１に記載の残響除去方法であって、
各チャネルの離散時間ｎの上記後部残響除去信号推定値を、それぞれのチャネルに対して算出された上記遅延量だけ遅延させる遅延過程と、
上記遅延過程で遅延させた離散時間ｎの上記後部残響除去信号推定値の和を、補正残響除去信号値として算出する遅延補正過程と、を有する。
【請求項２３】
請求項１３に記載の残響除去方法であって、
Ｍ≧２である。
【請求項２４】
請求項１３に記載の残響除去方法であって、
上記マルチステップ線形予測モデルは、
ｘｗ（ｎ）をチャネルｗ（ｗ＝１，...，Ｍ）に対応する離散時間ｎの離散音響信号値とし、ｘｍ（ｎ）をチャネルｍ（ｍ＝１，...，Ｍ）に対応する離散時間ｎの離散音響信号値とし、ｅｗ（ｎ）をチャネルｗ及び離散時間ｎに対応する予測誤差とし、Ｎを正の整数とし、[・]をガウス記号とし、αｗ，ｍ（ｐ）を、ｘｗ（ｎ）に対応するｐ番目の線形予測係数とし、Ｄをステップサイズを示す定数とした場合における、
【数２５】
0004774100.tif 000027
である。
【請求項２５】
請求項１３に記載された残響除去方法の各過程をコンピュータに実行させるための残響除去プログラム。
【請求項２６】
請求項２５に記載の残響除去プログラムを格納したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号処理の技術分野に係わり、特に、残響を含む音響信号から残響を除去する技術に関する。
【背景技術】
【０００２】
残響のある環境で音響信号を収音すると、本来の信号に残響が重畳された信号が観測される。この場合、重畳した残響成分によって音響信号の明瞭性が大きく低下し、音響信号本来の性質を抽出することが困難となる。例えば、残響を含む音声信号を自動音声認識（以下、音声認識）システムによって認識した場合、この残響の影響によって、音声認識システムの認識率は著しく低下してしまう。
残響除去処理は、このような場合に、重畳した残響を取り除き、音響信号を本来の音質に戻すことができる技術である。これにより、例えば、音声信号の明瞭性を回復し、音声認識率等を改善することが可能となる。
【０００３】
長い残響を除去する残響除去処理の従来例として、非特許文献１に示す方法がある。
この従来例では、後部残響が指数関数的に減衰すると仮定し、指数関数を用いて後部残響のエネルギーを推定し、残響除去を行う。すなわち、観測信号を周波数領域信号に変換し、各周波数において周波数領域信号が直接音と後部残響との和であると仮定し、残響エネルギーが指数関数的に減少するモデル（multi-band decay model）を用いて、そのモデルパラメータを推定し、推定された後部残響エネルギーと観測信号エネルギーとにスペクトル減算法（Spectral subtraction）を適用し、後部残響を除去している。なお、モデルパラメータの推定には、後部残響のみが観測される音声の末尾部分を用いている。
【非特許文献１】I. Tashev and D. Allred, "Reverberation Reduction for Improved Speech Recognition" 2005 Joint Workshop on hands-Free Speech Communication and Microphone Arrays.
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかし、上述の従来例では、環境によって、精度の良い残響除去を自動的に行うことができない場合があるという問題点があった。
すなわち、上述の従来例の場合、用いた指数関数が観測音中の後部残響を良くモデリングできていれば、後部残響を精度良く除去できる。しかし、一般的に部屋、話者とマイク間距離により多様に変化する残響を正確に指数関数で近似することは難しい。また、音響信号の末尾部分とそれ以外の部分で残響特性が異なる場合は、この枠組みを用いることはできない。このように、上述の従来例では、後部残響を指数関数でモデル化しているため、そのモデルと観測音中の後部残響がうまく合わない場合や、音響信号の末尾部分とそれ以外の部分で残響特性が異なる場合に精度の良い残響除去を自動的に行うことはできなかった。
【０００５】
本発明はこのような点に鑑みてなされたものであり、どのような環境でも精度の良い残響除去処理を行うことが可能な残響除去技術を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明では、上述の課題を解決するために、まず、Ｍ（Ｍ≧１）個のセンサによってそれぞれ観測されたＭ個のチャネルｍ（ｍ＝１，...，Ｍ）の音響信号をそれぞれ複数の時点でサンプリングして得られた離散音響信号値をメモリに記憶する離散音響信号記憶過程と、チャネルｗ（ｗ＝１，...，Ｍ）のマルチステップ線形予測モデル（長時間区間におけるＭ個のチャネルの各離散音響信号値を線形結合した線形予測項と、予測誤差項と、の和によって、当該長時間区間より所定時間後の離散時間ｎにおけるチャネルｗの離散音響信号値を表現した線形予測モデル）の各線形予測係数を、複数の上記離散音響信号値を用いて算出するモデル適用過程と、上記チャネルｗのマルチステップ線形予測モデルの各線形予測係数と複数の上記離散音響信号値とを、上記チャネルｗのマルチステップ線形予測モデルの上記線形予測項に代入して得られた線形予測値を、離散時間ｎにおけるチャネルｗの後部残響予測値として出力する後部残響予測過程とが実行される。
【０００７】
なお、「モデル」とは、物理的な状態を表現した関係式を意味する。また、「チャネル」とは、センサ毎の処理系列を意味し、同じセンサで観測された音響信号の処理系列は同じチャネルに属する。また、Ｍ≧２の場合、センサのうちの少なくとも一部は、他のセンサと異なる位置に配置される。すなわち、Ｍ≧２の場合、全てのセンサが全く同じ位置に配置されることはない。また、Ｍ≧２の場合における「マルチステップ線形予測モデル」を「マルチチャネルマルチステップ線形予測モデル」と呼ぶ。すなわち、本発明における「マルチステップ線形予測モデル」は、「マルチチャネルマルチステップ線形予測モデル」の上位概念である。
【０００８】
このように本発明では、マルチステップ線形予測モデルの線形予測項から、音響信号の後部残響予測値を算出し（詳細は後述）、当該後部残響予測値を用いることにより音響信号の残響を除去する。
ここで、本発明では、マルチステップ線形予測モデルによって音響信号をモデル化するため、指数関数のみで後部残響をモデル化する場合に比べ、より正確な近似が可能である。すなわち、音響信号の後部残響成分は、過去の各時点の音響信号に起因するものであり、或る時間区間において自己相関性を持つ。そのため、各時点の時系列データを複数時点の時系列データの線形結合で表現するマルチステップ線形予測モデルは、残響信号のモデル化に適しているといえる。
【０００９】
また、本発明のマルチステップ線形予測モデルは、長時間区間におけるＭ（Ｍ≧１）個のチャネルの各離散音響信号値を線形結合した線形予測項と、予測誤差項と、の和によって、当該長時間区間より所定時間後の離散時間ｎにおけるチャネルｗの離散音響信号値を表現した線形予測モデルである。このような「当該長時間区間より所定時間後の離散音響信号値を表現した線形予測モデル」を用いることで、「当該長時間区間直後の離散音響信号値を表現した線形予測モデル」を用いる場合よりも、後部残響信号を精度良く推定できる（詳細は後述）。その結果、適切な後部残響除去を行うことができる。
【００１０】
また、本発明では、マルチステップ線形予測モデルのモデルパラメータの推定に、音響信号の末尾部分のみではなく、離散音響信号の全体を用いるため、音響信号の末尾部分とそれ以外の部分で残響特性が異なる場合でも、適切な後部残響除去を行うことができる。
また、本発明において好ましくはＭ≧２である。マルチチャネルマルチステップ線形予測モデルを用いることにより、室内伝達関数中の最大位相成分が多い環境であっても、精度良く後部残響除去を行うことが可能となる（詳細は後述）。
また、本発明において好ましくは、上記モデル適用過程は、各離散時間の上記離散音響信号値から、当該離散時間直前の短時間区間内の各離散音響信号値と自己相関性を持つ自己相関成分（short-term correlation）を抑制し、擬似白色化（pre-whitening）した離散音響信号値を生成する擬似白色化過程と、上記擬似白色化した離散音響信号値を用い、上記マルチステップ線形予測モデルの各線形予測係数を算出する第１線形予測係数算出過程と、を有する。ここで、上記短時間区間は、上記長時間区間よりも短い。
【００１１】
マルチステップ線形予測モデルは、後部残響成分を線形予測項として推定するモデルである。よって、その線形予測係数の算出に用いる離散音響信号値の自己相関成分は、後部残響成分に起因するもののみであることが理想的である。しかし、音響信号の直接音成分は、後部残響成分に比べ極めて短い時間区間での自己相関性を持つ。擬似白色化は、この短い時間区間での自己相関性を抑制する処理である。擬似白色化した離散音響信号を用いてマルチステップ線形予測モデルの各線形予測係数を算出することは、後部残響除去処理の精度を向上させることになる。
【００１２】
また、本発明において好ましくは、上記擬似白色化過程は、上記短時間区間におけるチャネルｍの各離散音響信号値を線形結合した線形予測項と、予測誤差項と、の和によって、上記短時間区間直後の離散時間ｎにおける当該チャネルｍの離散音響信号値を表現した線形予測モデルである、チャネルｍの短時間線形予測モデルの各線形予測係数を、上記離散音響信号値を用いて算出する第２線形予測係数算出過程と、上記第２線形予測係数算出過程で算出された上記各線形予測係数をチャネルｍの上記短時間線形予測モデルに代入して得られる逆フィルタに当該チャネルｍの上記離散音響信号値を代入し、それによって得られる当該短時間線形予測モデルの上記予測誤差項の値を当該チャネルｍの上記擬似白色化した離散音響信号値として出力する逆フィルタ処理過程と、を有する。これにより、観測された音響信号の直接音成分の自己相関性を抑制できる。
【００１３】
また、この場合に好ましくは、Ｍ≧２であり、上記第２線形予測係数算出過程は、チャネル毎に上記離散音響信号値の自己相関係数を算出する自己相関係数算出過程と、チャネル毎に求められた上記自己相関係数をチャネル間で平均した平均自己相関係数を算出する自己相関係数平均化過程と、上記平均自己相関係数を用い、上記短時間線形予測モデルの各線形予測係数を算出する方程式演算過程と、を有する。
このように、各チャネルで求められた自己相関関数をチャネル間で平均した平均自己相関係数を、短時間線形予測モデルの各線形予測係数の算出に用いることにより、短時間線形予測モデルの各線形予測係数の算出精度が向上し、擬似白色化した離散音響信号の生成精度が向上する。これは、後部残響除去処理の精度向上に貢献する。
【００１４】
上述の平均自己相関係数を用いる代わりに、上記Ｍ個のセンサのうち、音響信号の音源に最も近い１つのセンサで観測された音響信号を複数の時点でサンプリングして得られた離散音響信号値の自己相関係数を算出し、上記自己相関係数を用い、上記短時間線形予測モデルの各線形予測係数を算出してもよい。これにより、精度良く短時間線形予測モデルの各線形予測係数を算出することができ、擬似白色化した離散音響信号の生成精度が向上する。これは、後部残響除去処理の精度向上に貢献する。
また、本発明において好ましくは、各チャネルの上記離散音響信号値を周波数領域の離散音響信号値に変換し、各チャネルの上記後部残響予測値を周波数領域の後部残響予測値に変換する周波数領域変換過程と、上記周波数領域の離散音響信号値の振幅スペクトルと、上記周波数領域の後部残響予測値の振幅スペクトルとの相対値をチャネル毎に求め、当該相対値を各チャネルの後部残響除去信号値の振幅スペクトル予測値として出力する後部残響除去過程と、を有する。これにより、後部残響除去信号値の振幅スペクトル予測値を算出できる。
【００１５】
また、この場合において好ましくは、上記後部残響予測過程では、上記モデル適用部で算出された上記各線形予測係数と複数の擬似白色化された上記離散音響信号値とを上記線形予測項に代入して得られた線形予測値を、離散時間ｎにおけるチャネルｗの後部残響予測値として算出し、上記周波数領域変換部では、擬似白色化された各チャネルの上記離散音響信号値を周波数領域の離散音響信号値に変換する。これにより、擬似白色化された後部残響除去信号値の振幅スペクトル予測値を算出することができる。このような値は、擬似白色化されたデータを必要とするシステムの入力値として好適である。
【００１６】
また、本発明において好ましくは、チャネルｗの上記後部残響除去信号値の振幅スペクトル予測値と、チャネルｗの上記周波数領域の離散音響信号値の位相情報とを用い、チャネルｗの後部残響除去信号値の複素スペクトル予測値を算出する複素スペクトル生成過程と、チャネルｗの上記後部残響除去信号値の複素スペクトル予測値を時間領域に変換したチャネルｗの後部残響除去信号推定値を算出する時間領域変換過程とを、さらに有する。これにより、後部残響が除去された音響信号の推定値（後部残響除去信号推定値）を得ることができる。
【００１７】
また、この場合において好ましくは、Ｍ≧２であり、上記モデル適用過程は、複数のチャネルに対してそれぞれ上記各線形予測係数を算出する過程であり、上記後部残響予測過程は、複数のチャネルに対してそれぞれ上記後部残響予測値を算出する過程であり、上記後部残響除去過程は、複数のチャネルに対してそれぞれ上記後部残響除去信号値の振幅スペクトル予測値を算出する過程であり、上記複素スペクトル生成過程は、複数のチャネルに対してそれぞれ上記後部残響除去信号値の複素スペクトル予測値を算出する過程であり、上記時間領域変換過程は、複数のチャネルに対してそれぞれ上記後部残響除去信号推定値を算出する過程であり、当該残響除去方法は、各チャネルの上記後部残響除去信号推定値をそれぞれ或る遅延量で遅延させた場合に、遅延後の各チャネルの上記後部残響除去信号推定値のチャネル間相互相関が極大となる、各チャネルの当該遅延量を決定する遅延量算出過程を有する。なお、遅延させない後部残響除去信号推定値については遅延量０と決定する。これにより、各チャネルで算出された後部残響除去信号推定値のチャネル相互での遅延量を補正することが可能となる。
【００１８】
そして、この場合には、各チャネルの上記後部残響除去信号推定値を、それぞれのチャネルに対して算出された上記遅延量だけ遅延させ、遅延させた上記後部残響除去信号推定値（遅延量０の後部残響除去信号推定値も含む）の和を、補正残響除去信号値として算出する。これにより、残響除去信号の推定精度が向上する。
また、本発明において好ましくは、上記マルチステップ線形予測モデルは、
ｘｗ（ｎ）をチャネルｗ（ｗ＝１，...，Ｍ）に対応する離散時間ｎの離散音響信号値とし、ｘｍ（ｎ）をチャネルｍ（ｍ＝１，...，Ｍ）に対応する離散時間ｎの離散音響信号値とし、ｅｗ（ｎ）をチャネルｗ及び離散時間ｎに対応する予測誤差とし、Ｎを正の整数とし、[・]をガウス記号とし、αｗ，ｍ（ｐ）を、ｘｗ（ｎ）に対応するｐ番目の線形予測係数とし、Ｄをステップサイズ（遅延）を示す定数とした場合における、
【数１】
0004774100.tif 000002
である。
【発明の効果】
【００１９】
以上のように、本発明では、どのような環境でも精度の良い残響除去を行うことができる。
【発明を実施するための最良の形態】
【００２０】
以下、本発明の実施の形態を図面を参照して説明する。
〔原理１〕
まず、本発明の原理１を説明する。なお、ここでは、原則、ｚ変換表記を用いる（後述の式（５）以外）。
原音響信号s(z)をモデル化すると、以下の式（１）のような、白色信号u(z)（予測誤差項）と短い（次数が小さな）自己回帰（AR: Auto-Regressive）モデルd(z)=1/(1-β(z))とを掛けたものとなる。すなわち、原音響信号s(z)は、d(z)に従った短期的な自己相関性を持つ。
【００２１】
s(z)=u(z)/(1-β(z))
=u(z)・d(z)      …(1)
なお、β(z)は以下のようなＡＲ多項式である。ここでｑは線形予測の次数であり、b(i)はｉ番目の項の線形予測係数である。
【数２】
0004774100.tif 000003
この場合、センサ（例えば、マイクロフォン）で観測される音響信号ｘ1（ｚ）は、
x1(z)=u(z)・[d(z)・ｈ(z)]
=u(z)・g(z)      …(2)
と表される。なお、h(z)は音響信号源からセンサまでの室内伝達関数を示し、g(z)は、
g(z)=d(z)・h(z)      …(3)
を満たすものとする。なお、g(z)を合成伝達関数と呼ぶ。
【００２２】
これより、原音響信号s(z)のd(z)に従った短期的な自己相関性を無視できるのであれば（原音響信号s(z)を白色信号とみなすことができるのであれば）、式（２）は、以下のように近似できる。
x1(z)≒u(z)・h(z)     …(4)
すなわち、d(z)を１とみなすことができるのであれば、センサで観測される音響信号x(z)は、白色信号u(z)に、d(z)よりも長い（次数が大きな）自己回帰過程が掛かったものとしてモデル化できる。
【００２３】
式（４）のようなモデル化は、d(z)を１に近づけられるほど適切なものとなる。原理１では、擬似白色化（Pre-whitening）処理により、センサで観測される音響信号x(z)の短期的な自己相関性を抑制する。これにより、センサで観測される音響信号x(z)のd(z)を１に近づけ、式（４）のモデルへの適用を適切なものとする。しかし、d(z)の自己相関があまり強くない原音響信号s(z)を対象とするのであれば、擬似白色化処理を行わなくても、式（４）のモデルへの適用はある程度適切なものとなる。
この場合、観測される音響信号は、センサ数Ｍ＝１の場合のマルチステップ線形予測モデル（長時間区間における各離散音響信号値を線形結合した線形予測項と、予測誤差項と、の和によって、当該長時間区間より所定時間後の離散音響信号値を表現した線形予測モデル）で近似できる。すなわち、観測される音響信号は、例えば、以下の式（５）のようにモデル化できる。なお、ｎを離散時間とし、x1(n)を（必要に応じて短時間相関を取り除いた）離散時間ｎに対応する音響信号とし、α1,1(p)を線形予測係数とし、Ｎを線形予測係数の数とし、Ｄをステップサイズ（遅延）とし、e1(n)を離散時間ｎに対応する予測誤差とする。また、式（５）の場合、(n-N-D)以上(n-1-D)以下の時間区間がマルチステップ線形予測モデルの定義中の「長時間区間」に相当し、x1(n)が「長時間区間より所定時間後の離散音響信号値」に相当する。
【００２４】
【数３】
0004774100.tif 000004
ここで、式（２）のg(z)の直接音成分をgd(z)とし、後部残饗成分をgr(z)とすると、以下の式（６）が仮定できる。
g(z):=gd(z)+z-D・gr(z)     …(6)
この場合、式（５）の線形予測係数α1,1(p)をｚ変換したα(z)は、以下の式（７）のように表せる。
【００２５】
【数４】
0004774100.tif 000005
と仮定する。なお、gmin(z)とgmax(z)とは、それぞれg(z)の最小位相成分（Ｚ平面上の単位円内のゼロ点に対応する成分）と最大位相成分（Ｚ平面上の単位円外のゼロ点に対応する成分）を意味する。また、前述のようにpre-whitening処理によってg(z)は室内伝達関数h(z)に近似する。また、一般に室内伝達関数h(z)は非最小位相である（最小位相成分だけではなく最大位相成分も有する）。そのため、g(z):=gmin(z)・gmax(z)との仮定は妥当である。また、min[gmax(z)]は、gmax(z)を最小位相化したものを意味し、gmax(z)の全てのゼロ点がそれらの共役逆であるＺ平面上の単位円内にプロットされたものを意味する。
【００２６】
ここで、式（７）の線形予測係数α(z)をpre-whitening処理された観測信号x'(z)=u(z)・g(z)に掛け、式（７）〜（９）を代入すると、以下の式（１０）が得られる。
【００２７】
【数５】
0004774100.tif 000006
ここで、gmax(z)／min[gmax(z)]はオールパスフィルタ（掛け合わされる信号の振幅はそのまま保持され、位相が変更され得る関数）となる。よって、
【数６】
0004774100.tif 000007
は、u(z)と同じ分散を持つ白色信号であることがわかる。また、式（６）に示した通り、z-D・gr(z)は、後部残響成分に起因する成分である。そして、音響信号は、式（２）のように表されるのだから、式（１０）の振幅スペクトル値は、音響信号の後部残響成分の振幅スペクトル値に近似する。すなわち、pre-whitening処理によって短時間相関を取り除いた音響信号を用い、式（５）のマルチステップ線形予測モデルの数値的最適化を行って線形予測係数を求め、これを観測信号に掛け合わせることで、後部残響成分に振幅スペクトル値が近似する信号を求めることができる。そして、この推定した振幅スペクトルを、観測された音響信号x(z)の振幅スペクトルから減算することによって、後部残響が除去された振幅スペクトルを得ることができる。
【００２８】
〔原理２〕
しかし、原理１の方法では、精度の良い後部残響除去を行うことができない場合がある。以下にその理由を説明する。
上述した式（１０）は分析フレーム長を無限長とした場合に成立する関係式である。有限長の分析フレーム単位でみた場合、式（１０）は完全に成立するとは限らない。式（１０）の右辺のオールパスフィルタgmax(z)/min[gmax(z)]の応答は、室内伝達関数h(z)中の最大位相成分が多くなるほど長くなる。そのため、有限長の分析フレーム単位でみた場合、室内伝達関数h(z)中の最大位相成分が多くなりオールパスフィルタgmax(z)/min[gmax(z)]の応答が長くなるほど、式（１０）の左辺と右辺との乖離が大きくなる。
【００２９】
通常、原理１の方法は有限長の分析フレーム単位で実行される。室内伝達関数h(z)中の最大位相成分が少ないのであれば、分析フレーム単位でみた式（１０）も比較的よく近似され、原理１の方法によって精度良く後部残響除去を行うことができる。しかし、一般に、信号源とセンサとの距離が遠くなるほど、室内伝達関数h(z)中の最大位相成分は増加する。この場合には、分析フレーム単位でみた式（１０）の近似が成り立たなくなり、式（１０）が成立することを前提とする原理１の方法による後部残響除去の精度は低下する。
すなわち、信号源からセンサまでの室内伝達関数中に最大位相成分が存在する場合、１つのセンサのみで観測された単一チャネルの音響信号を用いて完全な逆フィルタを構成することはできない。ここで、逆フィルタを構成可能であることと、上記の線形予測係数を算出できることは等価である。このことは、例えば、「M. Miyoshi and Y. Kaneda, "Inverse Filtering of Room Acoustics, " IEEE Trans. on Acoustics," Speech and Signal Processing, 36(2), pp. 145-152, 1988（以下『参考文献１』という）に開示されている。一方、参考文献１には、信号源からセンサまでの室内伝達関数中に最大位相成分が存在する場合に、複数のセンサで観測された複数チャネルの音響信号を用いて逆フィルタを構成できることが開示されている。
【００３０】
そこで原理２では、線形予測モデルとして上述のマルチチャネルマルチステップ線形予測モデル（センサ数Ｍ≧２の場合のマルチステップ線形予測モデル）を採用し、複数のセンサによって観測された複数チャネルの音響信号を用いて線形予測係数を求める。これにより、線形予測モデルを用いた後部残響除去処理が可能となる。なお、原理２の手法において線形予測フィルタを求めることができることは、室内伝達関数中の最大位相成分の大小に依存しない。以下、この原理２を説明する。
センサｍ（ｍ＝１，．．．，Ｍ）で観測される音響信号xm(z)は、以下のようにモデル化できる。なお、gm(z)=d(z)・hm(z)を満たす合成伝達関数とし、hm(z)を音響信号源からセンサｍまでの室内伝達関数とする。
【００３１】
xm(z)=u(z)・(d(z)・hm(z))
=u(z)・gm(z)      …(11)
これより、原音響信号s(z)のd(z)に従った短期的な自己相関性を無視できるのであれば（原音響信号s(z)を白色信号とみなすことができるのであれば）、式（１１）は、以下のように近似できる。
xm(z)≒u(z)・hm(z)     …(12)
すなわち、d(z)を１とみなすことができるのであれば、センサｍで観測される音響信号xm(z)は、白色信号u(z)に、d(z)よりも長い（次数が大きな）自己回帰過程が掛かったものとしてモデル化できる。
【００３２】
式（１２）のようなモデル化は、d(z)を１に近づけられるほど適切なものとなる。原理２では、擬似白色化（Pre-whitening）処理により、センサｍで観測される音響信号xm(z)の短期的な自己相関性を抑制する。これにより、センサｍで観測される音響信号xm(z)のd(z)を１に近づけ、式（１２）のモデルへの適用を適切なものとする。しかし、d(z)の自己相関があまり強くない原音響信号s(z)を対象とするのであれば、擬似白色化処理を行わなくても、式（１２）のモデルへの適用はある程度適切なものとなる。
一方、室内伝達関数hm(z)には、最大位相成分が存在する。前述した『参考文献１』で開示されているように、室内伝達関数hm(z)に最大位相成分が存在する場合、１つのセンサのみで観測された単一チャネルの音響信号のみを用いて逆フィルタを構成することはできず、複数チャネルの音響信号を用いた場合にのみ逆フィルタを構成できる。これを式（１２）のモデルに当てはめて考えると、複数チャネルの音響信号を対象として式（１２）のモデルを構成した場合にのみ、正しい線形予測係数を算出でき、正確な後部残響信号を予測し、除去することができるといえる。以上より、原理２では、マルチチャネルマルチステップ線形予測モデル（長時間区間におけるＭ個（Ｍ≧２）のチャネルの各離散音響信号値を線形結合した線形予測項と、予測誤差項と、の和によって、当該長時間区間より所定時間後の離散時間における離散音響信号値を表現した線形予測モデル）によって、式（１２）のモデルを構成する。マルチチャネルマルチステップ線形予測モデルとしては、以下の式（１３）を例示できる。なお、式（１３）のマルチチャネルマルチステップ線形予測モデルは、時間領域のものである。
【００３３】
【数７】
0004774100.tif 000008
ここで、ｎは離散時間を示す。また、ｘｗ（ｎ）は、センサｗ（ｗ＝１，...，Ｍ）で観測され、離散時間ｎにおいてサンプリングされたチャネルｗの離散音響信号、又は、それらを擬似白色化した離散音響信号に相当する。また、xm(n)は、センサｍで観測され、離散時間ｎにおいてサンプリングされた離散音響信号、又は、それらを擬似白色化した離散音響信号に相当する。さらに、ew(n)は、センサｗ及び離散時間ｎに対応する予測誤差であり、当該予測誤差項以外の式（１３）右辺の項が線形予測項に相当する。また、Ｍはセンサの総数を示す。Ｎ−１はフィルタ長（線形予測次数）である。なお、[・]はガウス記号であり、・を超えない最大の整数を示す。また、αw,m(p)は、xw(n)に対応するｐ番目の線形予測係数である。また、Ｄはステップサイズ（遅延）を示す。また、式（１３）の場合、(n-[N/M]-D)以上(n-1-D)以下の時間区間が、マルチチャネルマルチステップ線形予測モデルの定義に示した「長時間区間」に相当し、xw(n)が「長時間区間より所定時間後の離散音響信号値」に相当する。
【００３４】
原理２では、複数チャネルの離散音響信号又はそれらを擬似白色化した離散音響信号を用い、このようなマルチチャネルマルチステップ線形予測モデルの線形予測係数αw,m(p)を求め、マルチチャネルマルチステップ線形予測モデルの線形予測項の値を求める。この線形予測項の値が、後部残響成分の予測値（後部残響予測値）となる。その後、周波数領域における離散音響信号の振幅スペクトルと後部残響予測値の振幅スペクトルとの相対値を求め、それを後部残響除去信号の振幅スペクトル予測値とする（例えば、「S. F. Boll, "Suppression of acoustic noise in speech using spectral subtraction," IEEE Trans. on Acoustics, Speech and Signal Processing, 27(2), pp. 113-120, 1979」参照）。このような方法により、音声認識処理等にとって重要な、直接音成分の振幅スペクトルを精度良く抽出できる。
【００３５】
〔別の観点からの原理１，２の説明〕
上記の原理１，２を別の観点から説明する。
［問題設定］
まず、時間領域での離散時間ｎに対応する原音響信号ｓ（ｎ）を、Ｖ次のＦＩＲフィルタd(k)と白色信号u(n)を用いて以下のように表現する。
【００３６】
【数８】
0004774100.tif 000009
これを用いると、センサｍで観測された音響信号xm(n)は以下のようにモデル化できる。
【００３７】
【数９】
0004774100.tif 000010
ただし、
【００３８】
【数１０】
0004774100.tif 000011
であり、hm(n)はセンサｍと音源との間のインパルス応答を表す。
式（１５）を行列によって書き換えると以下のようになる。
Xｍ(n)=Gｍ・U(n)       …(17)
なお、
U(n)=[u(n),u(n-1),...,u(n-T-N+1)]T
Xｍ(n)=[xｍ(n),xｍ(n-1),...,xｍ(n-N)]T
gｍ=[gｍ(0),gｍ(1),...,gｍ(T-1)]
【００３９】
【数１１】
0004774100.tif 000012
である、また、[・]Tは行列・の転置を意味する。
［Ｍ＝１（原理１）の場合の後部残響の推定］
前述のように原理１では、観測される音響信号を例えば式（５）のようにモデル化する。式（５）において予測誤差成分e1(n)のエネルギーを最小にする線形予測係数α1,1(p)を求めることは、以下の正規方程式を解くことと等価である。
【００４０】
(E{x1(n-1-D)・x1T(n-1-D)})・Α=E{x1(n-1-D)・x1(n)}     …(18)
なお、E{・}は・の時間平均を示し、Αは式（５）の線形予測係数α1,1(p)の行列Α=[α1,1(1), α1,1(2),...,α1,1(N-1)]Tを示す。
よって、以下のようにΑを求めることができる。
Α=(E{x1(n-1-D)・x1T(n-1-D)})-1・E{x1(n-1-D)・x1(n)}       …(19)
式（１９）の(・)-1内を展開すると以下のようになる。
E{x1(n-1-D)・x1T(n-1-D)}=G1・E{U(n-1-D)・UT(n-1-D)}・G1T
=σu2・G1・G1T       …(20)
ここで、白色信号u(n)の自己相関行列は、E{U(n-1-D)・UT(n-1-D)}=σu2・Iとなると仮定した。なお、σu2はu(n)の分散を示し、Iは単位行列を示す。
【００４１】
また、式（１９）のE{x1(n-1-D)・x1(n)}は、以下のように展開できる。
E{x1(n-1-D)・x1(n)}=G1・E{U(n-1-D)・UT(n)}・g1T=σu2・G1・glate,1T    …(21)
なお、glate,1=[g(D),g(D+1),...,g(T-1),0,...,0]Tである。すなわち、glate,1はg(n)のＤ個目以降の要素を表し、後部残響に相当する。
式（２０）、（２１）を用い、式（１９）のΑは以下のように書き直すことができる。なお、(・)-1は・の逆行列を表す。
Α=(G1・G1 T)-1・G1・glate,1       …(22)
ここで、式（５）の線形予測係数の行列Αを音響信号の行列X1(n)に掛け合わせ、その２乗の平均をとると以下のようになる。
【００４２】
E{(X1T(n)・Α)2}
=‖ΑT・G1・E{U(n)・UT(n)}・G1T・Α‖
=‖σu2・ΑT・G1・G1T・Α‖         …(23)
=‖σu2・glate,1T・G1T・(G1・G1T)-1・G1・glate,1T‖
≦‖σu2・glate,1T‖・‖G1T・(G1・G1T)-1・G1‖・‖glate,1T‖       …(24)
=‖σu2・glate,1T‖2       …(25)
ここで‖・‖は行列・のＬ２ノルムを示す。また、式（２３）の導出のため、白色信号u(n)の自己相関行列は、E{U(n)・UT(n)}=σu2・Iのように展開されると仮定した。十分に長い音響信号に対してこの分解が成立する。また、式（２４）の導出には、式（２２）とコーシー・シュワルツの不等式を用いた。さらに、式（２５）の導出には、‖G1T・(G1・G1T)-1・G1‖が射影行列であり、そのノルムが１になることを用いた。
【００４３】
また、σu2は白色信号u(n)の分散を示し、glate,1は後部残響に対応する成分であるため、式（２５）の‖σu2・glate,1T‖2は後部残響成分のパワーを示す。よって、式（５）の線形予測係数の行列Αを音響信号の行列X1(n)に掛け合わせ、その２乗の平均をとった値は、後部残響成分のパワーを常に正確に推定する値であるとはいえないまでも、後部残響成分のパワーを過大推定する値ではない。
［Ｍ≧２（原理２）の場合の後部残響の推定］
マルチチャネルマルチステップ線形予測モデルは、前述の式（１３）で定式化できる。ここで、式（１３）の［N/M］をＬとすると、式（１３）は以下のようになる。
【００４４】
【数１２】
0004774100.tif 000013
式（２６）において予測誤差成分ew(n)のエネルギーを最小にする線形予測係数αw,m(p)を求めることは、以下の正規方程式を解くことと等価である。
(E{X(n-1-D)・XT(n-1-D)})・Αw=E{X(n-1-D)・X(n)}     …(27)
なお、X(n)=[X1T(n),X2T(n),...,XMT(n)]Tであり、Αwは式（２６）の線形予測係数αw,m(p)の行列でありΑw=[αw,1(1),...,αw,1(L),αw,2(1),...,αw,M(L)]Tである。
【００４５】
よって、Αwは以下のように得られる。なお、(・)+は、行列(・)のムーア・ペンローズ型一般化逆行列を示す。
Αw=(E{X(n-1-D)・XT(n-1-D)})+・E{X(n-1-D)・X(n)}     …(28)
Ｍ＝１の場合と同様に式（２８）を展開すると、Αwは以下のように変形できる。
Αw=(G・GT)+・G・glate,w
=(GT)+・glate,w      …(29)
なお、G=[G1T,G2T,...,GMT]Tであり、glate,w=[gw(D),gw(D+1),...,gw(T-1),0,...,0]Tである。また、Gは列フルランクである。
【００４６】
次に、推定された線形予測係数αw,m(p)の行列Αwを用いて、複数のセンサで観測された多チャンネルの音響信号から後部残響を推定する。そのために、ベクトルX(n)の転置と式（２６）の線形予測係数αw,m(p)の行列Αwとを掛け合わせると以下のようになる。
XT(n)・Αw=UT(n)・GT・Αw
=UT(n)・GT・(GT)+・glate,w      …(30)
=UT(n)・GT・G・(GT・G)-1・glate,w      …(31)
=UT(n)・glate,w      …(32)
なお、式（３０）の導出には式（２９）を用い、式（３１）の導出にはのムーア・ペンローズ型一般化逆行列の定義を用いた。ここで、UT(n)・glate,wは音響信号の後部残響成分を意味する。よって、式（２８）のベクトルの転置と式（２６）の線形予測係数αw,m(p)の行列Αwとを掛け合わせることにより、音響信号の後部残響成分を正確に推定できることがわかる。言い換えると、マルチチャネルマルチステップ線形モデルを用いることにより、常に正確に後部残響成分を推定することができることがわかる。
【００４７】
〔マルチステップ線形予測モデルのＤの値と音響信号の短時間相関〕
次に、マルチステップ線形予測モデルのＤの値と音響信号の短時間相関との関係について説明する。
原理１，２の方法は、式（４）、（１２）の近似が成り立つことを前提にした方法である。つまり、原理１，２の方法では、室内伝達関数hm(n)と式（１５）に示される合成伝達関数gm(n)との差（‖hm(n)‖−‖gm(n)‖）（m≧1）が十分小さい場合に、正確な後部残響除去ができる。
【００４８】
図２３Ａは室内伝達関数値hを縦軸にとり時間（ｍｓ）を横軸にとったグラフである。図２３Ｂは合成伝達関数値gを縦軸にとり時間（ｍｓ）を横軸にとったグラフである。また、図２３Ｃは室内伝達関数hと合成伝達関数gとのエネルギー差を縦軸にとり時間（ｍｓ）を横軸にとったグラフである。
図２３Ａに例示するように、室内伝達関数値hは時間の経過とともに指数減衰する。また、図２３Ｂに例示するように、合成伝達関数値gも時間の経過とともに指数減衰する。また、図２３Ｃに例示するように、室内伝達関数値hや合成伝達関数値gが大きな時間ではそれらのエネルギー差も大きく、室内伝達関数値hや合成伝達関数値gが小さな時間ではそれらのエネルギー差も小さい。すなわち、時間の経過とともに室内伝達関数hと合成伝達関数gとのエネルギー差も小さくなっていく。そして、「或る時間」の経過後には、当該エネルギー差は、音声信号全体のエネルギーに対して無視できるほど小さくなる（所定の閾値以下又は未満となる）。そのため、正確な後部残響除去を行うためには、式（５）、（１３）のマルチステップ線形予測モデルのＤを当該「或る時間」以上に設定することが望ましい。しかし、室内伝達関数hと合成伝達関数gとのエネルギー比やd(z)は未知であり当該「或る時間」も不定である。従って、一般に、経験則から当該「或る時間」を推測し、その推測に基づいてマルチステップ線形予測モデルのＤの値を設定することになる。そして、より望ましくは、この「或る時間」の推測が困難であることを想定し、前述の擬似白色化によってd(z)成分を抑制する。これにより、室内伝達関数hと合成伝達関数gとのエネルギー差を無視でき正確な後部残響除去が可能となり、Ｄの設定可能範囲が広がる。なお、一般にＤの下限値は１であるが、擬似白色化によってd(z)成分を十分抑制できるならＤ＝０であってもかまわない。
また、マルチステップ線形予測モデルのＤの上限値としては、離散時刻ｎでのxm(n)のn+1+D時点での残響成分が所定値（例えば、離散時刻ｎでのxm(n)よりも６０ｄＢ低い値）以上又は超える値を例示できる。
【００４９】
〔第１実施形態〕
次に、本発明の第１実施形態について説明する。第１実施形態はセンサ数ＭがＭ≧２の場合の実施形態である。
【００５０】
＜ハードウェア構成＞
図３は、本実施形態における残響除去装置１０のハードウェア構成を例示したブロック図である。
図３に例示するように、この例の残響除去装置１０は、ＣＰＵ（Central Processing Unit）１１、入力部１２、出力部１３、補助記憶装置１４、ＲＯＭ（Read Only Memory）１５、ＲＡＭ（Random Access Memory）１６及びバス１７を有している。
この例のＣＰＵ１１は、制御部１１ａ、演算部１１ｂ及びレジスタ１１ｃを有し、レジスタ１１ｃに読み込まれた各種プログラムに従って様々な演算処理を実行する。また、入力部１２は、データが入力される入力インターフェース、キーボード、マウス等であり、出力部１３は、データが出力される出力インターフェース等である。補助記憶装置１４は、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、残響除去装置１０としてコンピュータを機能させるためのプログラムが格納されるプログラム領域１４ａ及び各種データが格納されるデータ領域１４ｂを有している。また、ＲＡＭ１６は、ＳＲＡＭ (Static Random Access Memory)、ＤＲＡＭ (Dynamic Random Access Memory)等であり、上記のプログラムが格納されるプログラム領域１６ａ及び各種データが格納されるデータ領域１６ｂを有している。また、バス１７は、ＣＰＵ１１、入力部１２、出力部１３、補助記憶装置１４、ＲＯＭ１５及びＲＡＭ１６を通信可能に接続する。
なお、このようなハードウェアの具体例としては、例えば、パーソナルコンピュータの他、サーバ装置やワークステーション等を例示できる。
【００５１】
＜プログラム構成＞
上述のように、プログラム領域１４ａ，１６ａには、本形態の残響除去装置１０の各処理を実行するための残響除去プログラムが格納される。残響除去プログラムを構成する各プログラムは、単一のプログラム列として記載されていてもよく、また、少なくとも一部のプログラムが別個のモジュールとしてライブラリに格納されていてもよい。また、各プログラムが単体でそれぞれの機能を実現してもよいし、各プログラムがさらに他のライブラリを読み出して各機能を実現するものでもよい。
【００５２】
＜ハードウェアとプログラムとの協働＞
ＣＰＵ１１（図３）は、読み込まれたＯＳ（Operating System）プログラムに従い、補助記憶装置１４のプログラム領域１４ａに格納されている上述のプログラムをＲＡＭ１６のプログラム領域１６ａに書き込む。同様にＣＰＵ１１は、補助記憶装置１４のデータ領域１４ｂに格納されている各種データを、ＲＡＭ１６のデータ領域１６ｂに書き込む。そして、このプログラムやデータが書き込まれたＲＡＭ１６上のアドレスがＣＰＵ１１のレジスタ１１ｃに格納される。ＣＰＵ１１の制御部１１ａは、レジスタ１１ｃに格納されたこれらのアドレスを順次読み出し、読み出したアドレスが示すＲＡＭ１６上の領域からプログラムやデータを読み出し、そのプログラムが示す演算を演算部１１ｂに順次実行させ、その演算結果をレジスタ１１ｃに格納していく。
【００５３】
図１は、このようにＣＰＵ１１に上述のプログラムが読み込まれて実行されることにより構成される残響除去装置１０の機能構成を例示したブロック図である。また、図２Ａは、モデル適用部１０ｂの機能構成の詳細を例示したブロック図であり、図２Ｂは、遅延調節部１０ｉの機能構成の詳細を例示したブロック図である。
図１に例示するように、残響除去装置１０は、メモリ１０ａと、モデル適用部１０ｂと、後部残響予測部１０ｃと、周波数領域変換部１０ｄと、後部残響除去部１０ｅと、複素スペクトル生成部１０ｆと、時間領域変換部１０ｇと、遅延量算出部１０ｈと、遅延調節部１０ｉと、メモリ１０ｊと、制御部１０ｋとを有する。
【００５４】
また、図２Ａに例示するように、モデル適用部１０ｂは、擬似白色化部１００と第１線形予測係数算出部２００とを有しており、擬似白色化部１００は、第２線形予測係数算出部１１０と逆フィルタ処理部１２０とを有している。また、第２線形予測係数算出部１１０は、自己相関係数算出部１１１と、自己相関係数平均化部１１２と、方程式演算部１１３とを有している。また、図２Ｂに例示するように、遅延調節部１０ｉは、遅延部１０ｉａと、遅延補正部１０ｉｂとを有している。
ここで、メモリ１０ａ及びメモリ１０ｊは、補助記憶装置１４、ＲＡＭ１６、レジスタ１１ｃ、その他のバッファメモリやキャッシュメモリ等の何れか、あるいはこれらを併用した記憶領域に相当する。また、モデル適用部１０ｂ、後部残響予測部１０ｃ、周波数領域変換部１０ｄ、後部残響除去部１０ｅ、複素スペクトル生成部１０ｆ、時間領域変換部１０ｇ、遅延量算出部１０ｈ、遅延調節部１０ｉ及び制御部１０ｋは、ＣＰＵ１１に残響除去プログラムを実行させることにより構成されるものである。
【００５５】
また、本形態の残響除去装置１０は、制御部１０ｋの制御のもと各処理を実行する。また、特に示さない限り、演算過程の各データは、逐一、メモリ１０ｊに格納・読み出され、各演算処理が進められる。メモリ１０ａやメモリ１０ｊには、x1(n)、α1,2(p)等の各データが格納されるが、それらはデータ属性、下付き添え字の値〔例えば、データx1(n)の下付添え字「１」〕及び（・）内の各値〔例えば、データx1(n)のｎ〕に対応付けられて格納され、これらを指定することにより、対応するデータを抽出できるものとする。
【００５６】
＜残響除去処理＞
次に、本形態の残響除去処理について説明する。
図４，５は、本形態の残響除去処理の全体を説明するためのフローチャートである。また、図６Ａは、図４のステップＳ１（モデル適用過程）の詳細を説明するためのフローチャートであり、図６Ｂは、図６ＡのステップＳ２１（擬似白色化過程）の詳細を説明するためのフローチャートである。また、図７Ａは、図６ＢのステップＳ３１（第２線形予測係数算出過程）の詳細を説明するためのフローチャートであり、図７Ｂは、図４のステップＳ４の詳細を説明するためのフローチャートである。以下、これらの図を用い、本形態の残響除去処理を説明する。
【００５７】
［前処理］
まず、Ｍ（Ｍ≧２）個のセンサでそれぞれ観測されたＭ個のチャネルｗ（ｗ＝１，...，Ｍ）の音響信号が所定の標本化周波数でサンプリングされ、チャネル毎の離散音響信号値x1(n)…xM(n)が生成される。なお、ｎは離散時間を示す。生成された各チャネルの離散音響信号値x1(n)…xM(n)は、それぞれメモリ１０ａに格納される。なお、本形態では、残響除去を行う全時間区間の離散音響信号値x1(n)…xM(n)を事前に取得し、メモリ１０ａに格納しておき、分析フレーム毎に、以下の各過程を実行する。しかし、離散音響信号値x1(n)…xM(n)の取得をリアルタイムで行いつつ、以下の各過程を実行してもよい。
【００５８】
また、残響除去装置１０が最終的に出力する情報が、後部残響除去信号の振幅スペクトルのみであるのか、それとも位相成分をも有する音響信号なのかを示す情報をメモリ１０ｊに格納しておく。本形態では、フラグ（データ）δをメモリ１０ｊに格納しておく。そして、出力する情報が直接音の振幅スペクトルのみである場合δ＝１とし、位相成分をも有する音響信号である場合δ＝０とする。なお、残響除去装置１０が最終的に出力する情報が後部残響除去信号の振幅スペクトルのみでよい場合とは、例えば、残響除去装置１０が最終的に出力する情報を音声認識システムの入力情報として利用する場合を例示できる。
【００５９】
以下、本形態の後部残響除去処理を説明する。なお、以下では、１つの分析フレームの処理過程のみを説明するが、実際は複数の分析フレームに対して同様な処理が行われる。また、分析フレームとは、複数の離散時間ｎを含む時間区間を意味する。
［モデル適用過程（ステップＳ１）］
モデル適用過程では、モデル適用部１０ｂが、メモリ１０ａから読み込んだ１分析フレーム分のＭチャネルの離散音響信号値x1(n)…xM(n)を用い、式（１３）に示したマルチステップ線形予測モデルの各線形予測係数αw,1(p)…αw,M(p)を算出する（ステップＳ１）。以下、この処理の詳細を階層的に説明する。
【００６０】
［モデル適用過程（ステップＳ１）の詳細（図６Ａ）］
図６Ａに例示するように、モデル適用過程では、まず、擬似白色化（Pre-whitening）部１００（図２Ａ）が、入力された離散音響信号値x1(n)…xM(n)が有する短時間区間での自己相関成分を抑制し、擬似白色化した離散音響信号値x1’(n)…xM’(n)を生成して出力する（擬似白色化過程／ステップＳ２１）。すなわち、各離散時間の上記離散音響信号値x1(n)…xM(n)から、当該離散時間n直前の短時間区間内の各離散音響信号値と自己相関性を持つ自己相関成分を抑制し、擬似白色化した離散音響信号値x1’(n)…xM’(n)を生成する。
【００６１】
前述したように、マルチステップ線形予測モデルは、d(z)に従った短期的な自己相関（短時間区間での自己相関成分）が抑制された離散音響信号によく合致する。よって、このような短期的な自己相関を離散音響信号値x1(n)…xM(n)から抑制することは、後部残響を精度良く推定する上で望ましい。
次に、上述の擬似白色化された離散音響信号値x1’(n)…xM’(n)が第１線形予測係数算出部２００（図２Ａ）に入力され、第１線形予測係数算出部２００は、当該擬似白色化された離散音響信号値x1’(n)…xM’(n)を用い、式（１３）に示したマルチステップ線形予測モデルの各線形予測係数αw,1(p)…αw,M(p)を算出して出力する（第１線形予測係数算出過程／ステップＳ２２）。なお、一例として、式（１３）における遅延Ｄを、例えば３０ｍｓ（標本化周波数１２０００Ｈｚの場合、３６０タップに相当）とし、Ｎを例えば３０００程度とする。また、x1’(n)…xM’(n)を用いてαw,1(p)…αw,M(p)を算出する方法としては、自己相関法（correlation method）や共分散法（covariance method）を例示できる。また、MATLAB（登録商標）等を利用してこの処理を行ってもよい。
【００６２】
［擬似白色化過程（ステップＳ２１）の詳細（図６Ｂ）］
次に、擬似白色化過程（ステップＳ２１）の詳細を説明する。本形態では、一例として、線形予測によって擬似白色化過程を行う。まず、図６Ｂに例示するように、第２線形予測係数算出部１１０が、入力された離散音響信号値x1(n)…xM(n)を用い、短時間線形予測モデルの各線形予測係数b(1)…b(q)を算出して出力する（第２線形予測係数算出過程／ステップS３１）。なお、「短時間線形予測モデル」とは、短時間区間におけるチャネルｗの各離散音響信号値を線形結合した線形予測項と、予測誤差項と、の和によって、短時間区間直後の離散時間ｎにおける当該チャネルｗの離散音響信号値を表現した線形予測モデルを意味する。ここで、「短時間区間」は、マルチステップ線形予測モデルの定義で示した「長時間区間」よりも短い。本形態では、以下の短時間線形予測モデルを用いる。
【００６３】
【数１３】
0004774100.tif 000014
なお、式（３３）におけるxm'(n)の項が予測誤差項に相当し、それ以外の右辺の項が線形予測項に相当する。また、b(i)は、線形予測項のｉ番目の線形予測係数を意味する。また、式（３３）における短時間区間は、離散音響信号値x1(n)…xM(n)の短時間相関成分の系列長やパワーに応じて適宜設定すればよい。一例として、短時間区間を３０ｍｓ（標本化周波数１２０００Ｈｚの場合、ｑ＝３６０）程度とすることができる。この場合、以下のステップＳ３２により、短時間区間３０ｍｓ内で自己相関を持つ初期反射音成分や直接音成分が抑制できる。
【００６４】
次に、逆フィルタ処理部１２０（図２Ａ）に、各線形予測係数b(1)…b(q)と、離散音響信号値x1(n)…xM(n)とが入力される。逆フィルタ処理部１２０は、各線形予測係数b(1)…b(q)を短時間線形予測モデル（式（３３））に代入して得られる逆フィルタ
【００６５】
【数１４】
0004774100.tif 000015
に、離散音響信号値x1(n)…xM(n)を代入し、それによって得られる短時間線形予測モデルの予測誤差項の値を、擬似白色化した離散音響信号値x1’(n)…xM’(n)として算出して出力する（逆フィルタ処理過程／ステップＳ３２）。
［第２線形予測係数算出過程（ステップＳ３１）の詳細（図７Ａ）］
次に、第２線形予測係数算出過程（ステップＳ３１）の詳細を説明する。本形態の例では、自己相関法を用いて第２線形予測係数算出過程を実行する。しかし、共分散法などその他の公知の線形予測係数算出方法を用いて第２線形予測係数算出過程を実行してもよい。
【００６６】
まず、図７Ａに例示するように、まず、自己相関係数算出部１１１（図２Ａ）が、入力された離散音響信号値x1(n)…xM(n)を用い、チャネル毎に離散音響信号値x1(n)…xM(n)の自己相関係数c1(i)…cM(i)を算出して出力する（自己相関係数算出過程／ステップＳ４１）。具体的には、例えば、自己相関係数算出部１１１が、以下の式（３５）に従って、自己相関係数c1(i)…cM(i)を算出して出力する。なお、Ｔは、ｑ（式（３３）、（３４））よりも大きく、１分析フレームが有するサンプル数より小さな自然数である。また、以下の演算は、例えば、上述の離散音響信号値x1(n)…xM(n)に、ｎ＜０，ｎ≧Ｔの範囲で０となるような有限長の窓（ハミング窓など）を乗じた後に行われる。また、i=0,1,...,qである。
【００６７】
【数１５】
0004774100.tif 000016
次に、自己相関係数平均化部１１２（図２Ａ）に各チャネルの自己相関係数c1(i)…cM(i)が入力され、自己相関係数平均化部１１２は、これら自己相関係数c1(i)…cM(i)をチャネル間で平均した平均自己相関係数c(i)を算出して出力する（自己相関係数平均化過程／ステップＳ４２）。この平均自己相関係数c(i)の算出は、例えば、以下の式（３６）に従って行われる。
【００６８】
【数１６】
0004774100.tif 000017
次に、上述のように求められた各平均自己相関係数c(i)が方程式演算部１１３に入力され、方程式演算部１１３は、各平均自己相関係数c(i)を用い、以下のようにYule-Walkerの方程式（正規方程式）の解を求めることにより、短時間線形予測モデルの各線形予測係数b(1)…b(q)を算出して出力する（方程式演算過程／ステップＳ４３）。
【００６９】
【数１７】
0004774100.tif 000018
以上のように、チャネル毎に生成された自己相関係数c1(i)…cM(i)をチャネル間で平均した平均自己相関係数c(i)を用い、各線形予測係数b(1)…b(q)を算出する構成としたため、何れかのチャネルで生成された自己相関係数を用いる場合に比べ、線形予測係数b(1)…b(q)の算出精度が向上し、離散音響信号値x1(n)…xM(n)が具備するd(z)に従った短期的な自己相関をより効果的に抑制することができる。前述のように、これは、後部残響除去の精度向上につながる（モデル適用過程（ステップＳ１）の詳細の説明終わり）。
【００７０】
［後部残響予測過程（ステップＳ２）］
モデル適用過程（ステップＳ１）の後、後部残響予測部１０ｃに、メモリ１０ａから読み込まれた離散音響信号値x1(n)…xM(n)と、モデル適用過程（ステップＳ１）で算出された各線形予測係数αw,1(p)…αw,M(p)とが入力される。そして、後部残響予測部１０ｃは、各線形予測係数αw,1(p)…αw,M(p)と離散音響信号値x1(n)…xM(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値rw(n)(w=1,...,M)として算出して出力する（ステップＳ２）。本形態では、式（１３）のマルチステップ線形予測モデルを用いているため、後部残響予測部１０ｃは、以下の式（３８）に従って後部残響予測値rw(n)を求めて出力する。
【００７１】
【数１８】
0004774100.tif 000019
［周波数領域変換過程（ステップＳ３）］
次に、周波数領域変換部１０ｄに、メモリ１０ａから読み込まれた離散音響信号値x1(n)…xM(n)と、後部残響予測過程（ステップＳ２）で算出された後部残響予測値r1(n)…rM(n)とが入力される。周波数領域変換部１０ｄは、入力された離散音響信号値x1(n)…xM(n)を周波数領域の離散音響信号値X1(f,t)…XM(f,t)に変換し、後部残響予測値r1(n)…rM(n)を周波数領域の後部残響予測値R1(f,t)…RM(f,t)に変換する（ステップＳ３）。本形態では、例えば、窓長３０ｍｓのハニング窓などの有限長の窓関数を用い、短時間フーリエ変換（DFT: Discrete Fourier Transform）等によって、これらの周波数領域への変換を行う。周波数領域変換部１０ｄは、これらの処理により、周波数領域の離散音響信号値X1(f,t)…XM(f,t)の振幅スペクトル|X1(f,t)|…|XM(f,t)|と位相情報arg[X1(f,t)]…arg[XM(f,t)]、及び、周波数領域の後部残響予測値R1(f,t)…RM(f,t)の振幅スペクトル|R1(f,t)|…|RM(f,t)|と位相情報arg[R1(f,t)]…arg[RM(f,t)]とを抽出し、出力する。なお、arg[・]は、・の偏角を意味する。
【００７２】
［後部残響除去過程（ステップＳ４）］
次に、後部残響除去部１０ｅに、周波数領域の離散音響信号値の振幅スペクトル|X1(f,t)|…|XM(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とが入力される。そして、後部残響除去部１０ｅは、周波数領域の離散音響信号値の振幅スペクトル|X1(f,t)|…|XM(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|として出力する（ステップＳ４）。以下に、この処理の詳細を例示する。
【００７３】
［後部残響除去過程（ステップＳ４）の詳細（図７Ｂ）］
図７Ｂに例示するように、まず、後部残響除去部１０ｅが、振幅スペクトル|X1(f,t)|…|XM(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とを用い、各ｍ（ｍ＝１，．．．，Ｍ）に対して
|Xm(f,t)|k-|Rm(f,t)|k・const    …(39）
の演算を行い、各演算結果をメモリ１０ｊに格納する（ステップＳ５１）。なお、constは定数を示し、ｋは自然数を示す。本形態では、const=1.0とし、k=2とする。
【００７４】
次に、制御部１０ｋは、変数ｍに１を代入して、当該ｍをメモリ１０ｊに格納する（ステップＳ５２）。次に、制御部１０ｋは、メモリ１０ｊに格納された式（３９）の演算結果が以下の関係を満たすか否かを判断する（ステップＳ５３）。
|Xm(f,t)|k-|Rm(f,t)|k・const>0    …(40）
ここで、式（４０）の関係を満たすと判断された場合、制御部１０ｋは後部残響除去部１０ｅに命令を与え、後部残響除去部１０ｅは、
|Sm(f,t)|=(|Xm(f,t)|k-|Rm(f,t)|k・const)1/k    …(41）
によって、変数ｍに対応する|Sm(f,t)|を算出して出力する（ステップＳ５４）。一方、式（４０）の関係を満たさないと判断された場合、制御部１０ｋは後部残響除去部１０ｅに命令を与え、後部残響除去部１０ｅは、０又は十分小さい値を変数ｍに対応する|Sm(f,t)|として出力する（ステップＳ５５）。なお、ステップＳ５３〜Ｓ５５の処理は半波整流処理に相当する。ステップＳ５３〜Ｓ５５以外の方法によって半波整流を行ってもよい。
【００７５】
半波整流後、制御部１０ｋはメモリ１０ｊを参照し、変数ｍがＭであるか否かを判断する（ステップＳ５６）。ここで、ｍ＝Ｍでなければ、制御部１０ｋは、ｍ＋１を新たな変数ｍの値としてメモリ１０ｊに格納し（ステップＳ５７）、処理をステップＳ５３に戻す。一方、ｍ＝Ｍであれば、制御部１０ｋは、ステップＳ４の処理を終了させる（後部残響除去過程（ステップＳ４）の詳細の説明終わり）。
［フラグ判定過程（ステップＳ５，Ｓ６）］
ステップＳ４の後、制御部１０ｋが、メモリ１０ｊに格納されているフラグδを読み出し、そのフラグδが、振幅スペクトルのみを出力することを示すフラグであるか否か、すなわち、δ＝１であるか否かを判断する（ステップＳ５）。ここで、δ＝１であれば、制御部１０ｋは、後部残響除去過程（ステップＳ４）で後部残響除去部１０ｅが生成した後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|を、残響除去装置１０の最終的な出力情報として出力し（ステップＳ６）、当該分析フレームの処理を終了させる。このように出力された振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|は、例えば、残響除去装置１０の後段に続く音声認識システム等のアプリケーションに渡され、特徴量に変換される。
【００７６】
一方、δ＝０であれば、制御部１０ｋは、以下のステップＳ７以降の処理を実行させる。
［複素スペクトル生成過程（ステップＳ７）］
複素スペクトル生成過程では、まず、複素スペクトル生成部１０ｆに、後部残響除去部１０ｅから出力（ステップＳ４）された後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|と、周波数領域変換部１０ｄから出力（ステップＳ３）された周波数領域の離散音響信号値の位相情報arg[X1(f,t)]…arg[XM(f,t)]とが入力される。複素スペクトル生成部１０ｆは、これらの情報を用い、以下の式（４２）に従って、後部残響除去信号値の複素スペクトル予測値S1(f,t)…SM(f,t)を算出して出力する（ステップＳ７）。なお、exp(・)は、ネイピア数を底とした指数関数であり、ｊは虚数単位である。
【００７７】
Sm(f,t)=|Sm(f,t)|・exp(j・arg[Xm(f,t)])    …(42）
［時間領域変換過程（ステップＳ８）］
ステップＳ７の後、時間領域変換部１０ｇに、上述の後部残響除去信号値の複素スペクトル予測値S1(f,t)…SM(f,t)が入力される。そして、時間領域変換部１０ｇは、後部残響除去信号値の複素スペクトル予測値S1(f,t)…SM(f,t)を時間領域に変換した後部残響除去信号推定値s1(n)…sM(n)を算出して出力する（ステップＳ８）。なお、時間領域への変換は、例えば、逆フーリエ変換によって行う。
【００７８】
［遅延量算出過程（ステップＳ９）］
ステップＳ８の後、遅延量算出部１０ｈに、後部残響除去信号推定値s1(n)…sM(n)が入力される。そして、遅延量算出部１０ｈは、後部残響除去信号推定値s1(n)…sM(n)のチャネル間相互相関を極大にする後部残響除去信号推定値の遅延量τ1…τMを、各チャネルについて決定する（ステップＳ９）。以下にこの具体例を示す。
［遅延量算出過程（ステップＳ９）の具体例］
まず、遅延量算出部１０ｈは、入力された分析フレーム内の後部残響除去信号推定値s1(n)…sM(n)に対し、以下の式（４３）のようなチャネル間相関関数Ａｍ（τ）の関数値を求める。なお、E{・}は平均演算子である。
【００７９】
Am(τ)=E{s1(n)・sm(n+τ)}     …(43)
次に、遅延量算出部１０ｈは、各ｍについて、チャネル間相関関数Am(τ)を極大（例えば最大）とするτをτｍとして求める。例えば、チャネル間相関関数Am(τ)を最大とするτをτｍとする場合には、遅延量算出部１０ｈは、
τｍ=max{ Am(τ)}    …(44）
を算出して出力する。なお、max{・}は・の最大値を検出する。また、τｍは、チャネルｍの後部残響除去信号推定値の遅延量であり、遅延量にはτｍ＝０も含む（遅延量算出過程（ステップＳ９）の具体例の説明終わり）。
【００８０】
［遅延調節過程（ステップＳ１０，Ｓ１１）］
ステップＳ９の後、各遅延量τ1…τMと、後部残響除去信号推定値s1(n)…sM(n)とが、遅延調節部１０ｉ（図１）に入力される。そして、遅延調節部１０ｉの遅延部１０ｉａ（図２Ｂ）は、各チャネルの後部残響除去信号推定値s1(n)…sM(n)を、それぞれ遅延量τ1…τMだけ遅延させてs1(n+τ1)…sM(n+τM)を算出して出力する（ステップＳ１０）。
次に、s1(n+τ1)…sM(n+τM)が、遅延補正部１０ｉｂ（（図２Ｂ））に入力され、遅延補正部１０ｉｂは、以下の式（４５）に従い、s1(n+τ1)…sM(n+τM)の和を算出し（ステップＳ１１）、この和を補正残響除去信号値s(n)として出力して（ステップＳ１２）、当該分析フレームの処理を終了する。各チャネルの後部残響信号に含まれる誤差成分は統計的に独立であると過程した場合、この操作により誤差を抑圧できることになる。
【００８１】
【数１９】
0004774100.tif 000020
〔第２実施形態〕
次に、本発明の第２実施形態について説明する。第２実施形態は、第１実施形態の変形例である。
第１実施形態の［後部残響予測過程（ステップＳ２）］では、各線形予測係数αw,1(p)…αw,M(p)と擬似白色化していない離散音響信号値x1(n)…xM(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値rw(n)(w=1,...,M)として算出していた。また、［後部残響除去過程（ステップＳ４）］では、擬似白色化していない周波数領域の離散音響信号値の振幅スペクトル|X1(f,t)|…|XM(f,t)|と周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|としていた。
【００８２】
これに対し、第２実施形態では、［後部残響予測過程］において、各線形予測係数αw,1(p)…αw,M(p)と擬似白色化した離散音響信号値x1'(n)…xM'(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値rw(n)(w=1,...,M)として算出する。また、第２実施形態では、［後部残響除去過程］において、擬似白色化後の周波数領域の離散音響信号値の振幅スペクトル|X1'(f,t)|…|XM'(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|する。このようにして得られた後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|は、短時間相関成分が抑制された（擬似白色化された）ものとなる。そのため、このように得られた振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|は、例えば音声認識システムのように、擬似白色化されたデータが必要なシステムへの入力として好適である。このようなシステムにおいて、擬似白色化する前処理が不要となるからである。
【００８３】
これらが第１実施形態と第２実施形態との相違点である。以下では、第１実施形態との相違点を中心に説明し、第１実施形態と共通する事項については説明を省略する。
＜ハードウェア構成＞
第１実施形態で説明したのと同様である。
＜ハードウェアとプログラムとの協働＞
本形態の残響除去装置もコンピュータに所定のプログラムが読み込まれて実行されることにより構成される。図８は、本形態の残響除去装置３１０の機能構成を例示したブロック図である。また、図９は、モデル適用部３１０ｂの機能構成の詳細を例示したブロック図である。なお、図８，９において、第１実施形態と共通する部分については第１実施形態と同じ符号を用いた。
【００８４】
図８に例示するように、残響除去装置３１０は、メモリ１０ａと、モデル適用部３１０ｂと、後部残響予測部３１０ｃと、周波数領域変換部３１０ｄと、後部残響除去部３１０ｅと、複素スペクトル生成部３１０ｆと、時間領域変換部１０ｇと、遅延量算出部１０ｈと、遅延調節部１０ｉと、メモリ１０ｊと、制御部１０ｋとを有する。
また、図９に例示するように、モデル適用部３１０ｂは、擬似白色化部１００と第１線形予測係数算出部２００とを有しており、擬似白色化部１００は、第２線形予測係数算出部１１０と逆フィルタ処理部１２０とを有している。また、第２線形予測係数算出部１１０は、自己相関係数算出部１１１と、自己相関係数平均化部１１２と、方程式演算部１１３とを有している。モデル適用部３１０ｂと第１実施形態のモデル適用部１０ｂとの相違点は、モデル適用部３１０ｂの逆フィルタ処理部１２０が、擬似白色化した離散音響信号値x1’(n)…xM’(n)を後部残響予測部３１０ｃや周波数領域変換部３１０ｄにも転送する点である。
【００８５】
＜残響除去処理＞
次に、本形態の残響除去処理について説明する。
図１０，１１は、本形態の残響除去処理の全体を説明するためのフローチャートである。以下、これらの図を用い、本形態の残響除去処理を説明する。
［前処理］
第１実施形態と同様である。
［モデル適用過程（ステップＳ１０１）］
モデル適用過程では、モデル適用部３１０ｂが、メモリ１０ａから読み込んだ１分析フレーム分のＭチャネルの離散音響信号値x1(n)…xM(n)を用い、式（１３）に示したマルチステップ線形予測モデルの各線形予測係数αw,1(p)…αw,M(p)を算出する（ステップＳ１０１）。この処理は、第１実施形態の［モデル適用過程（ステップＳ１）］と同様であり、離散音響信号値x1(n)…xM(n)を擬似白色化する過程を含む。
【００８６】
［後部残響予測過程（ステップＳ１０２）］
モデル適用過程（ステップＳ１０１）の後、後部残響予測部３１０ｃに、モデル適用過程（ステップＳ１０１）で擬似白色化された離散音響信号値x1’(n)…xM’(n)と、モデル適用過程（ステップＳ１０１）で算出された各線形予測係数αw,1(p)…αw,M(p)とが入力される。
そして、後部残響予測部３１０ｃは、各線形予測係数αw,1(p)…αw,M(p)と擬似白色化された離散音響信号値x1’(n)…xM’(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値rw(n)(w=1,...,M)として算出して出力する（ステップＳ１０２）。式（１３）のマルチステップ線形予測モデルを用いていた場合、後部残響予測部３１０ｃは、以下の式（４６）に従って後部残響予測値rw(n)を求めて出力する。
【００８７】
【数２０】
0004774100.tif 000021
［周波数領域変換過程（ステップＳ１０３）］
次に、周波数領域変換部３１０ｄに、モデル適用過程（ステップＳ１０１）で擬似白色化された離散音響信号値x1’(n)…xM’(n)と、後部残響予測過程（ステップＳ１０２）で算出された後部残響予測値r1(n)…rM(n)とが入力される。周波数領域変換部３１０ｄは、入力された擬似白色化後の離散音響信号値x1’(n)…xM’(n)を周波数領域の離散音響信号値X1’(f,t)…XM’(f,t)に変換し、後部残響予測値r1(n)…rM(n)を周波数領域の後部残響予測値R1(f,t)…RM(f,t)に変換する（ステップＳ１０３）。周波数領域変換部３１０ｄは、これらの処理により、周波数領域の離散音響信号値X1’(f,t)…XM’(f,t)の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と位相情報arg[X1’(f,t)]…arg[XM’(f,t)]、及び、周波数領域の後部残響予測値R1(f,t)…RM(f,t)の振幅スペクトル|R1(f,t)|…|RM(f,t)|と位相情報arg[R1(f,t)]…arg[RM(f,t)]とを抽出し、出力する。
【００８８】
［後部残響除去過程（ステップＳ１０４）］
次に、後部残響除去部３１０ｅに、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とが入力される。そして、後部残響除去部３１０ｅは、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|として出力する（ステップＳ１０４）。
【００８９】
［フラグ判定過程（ステップＳ１０５，Ｓ１０６）］
ステップＳ１０４の後、制御部１０ｋが、メモリ１０ｊに格納されているフラグδを読み出し、そのフラグδが、振幅スペクトルのみを出力することを示すフラグであるか否か、すなわち、δ＝１であるか否かを判断する（ステップＳ１０５）。ここで、δ＝１であれば、制御部１０ｋは、後部残響除去過程（ステップＳ１０４）で後部残響除去部３１０ｅが生成した後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|を、残響除去装置３１０の最終的な出力情報として出力し（ステップＳ１０６）、当該分析フレームの処理を終了させる。一方、δ＝０であれば、制御部１０ｋは、以下のステップＳ１０７以降の処理を実行させる。
【００９０】
［複素スペクトル生成過程（ステップＳ１０７）］
複素スペクトル生成過程では、まず、複素スペクトル生成部３１０ｆに、後部残響除去部３１０ｅから出力（ステップＳ１０４）された後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|と、周波数領域変換部１０ｄから出力（ステップＳ３）された周波数領域の離散音響信号値の位相情報arg[X1’(f,t)]…arg[XM’(f,t)]とが入力される。複素スペクトル生成部３１０ｆは、これらの情報を用い、以下の式（４７）に従って、後部残響除去信号値の複素スペクトル予測値S1(f,t)…SM(f,t)を算出して出力する（ステップＳ１０７）。
【００９１】
Sm(f,t)=|Sm(f,t)|・exp(j・arg[Xm’(f,t)])    …(47)
［時間領域変換過程（ステップＳ１０８）・遅延量算出過程（ステップＳ１０９）・遅延調節過程（ステップＳ１１０，Ｓ１１１）］
時間領域変換過程（ステップＳ１０８）・遅延量算出過程（ステップＳ１０９）・遅延調節過程（ステップＳ１１０，Ｓ１１１）は、第１実施形態の時間領域変換過程（ステップＳ８）・遅延量算出過程（ステップＳ９）・遅延調節過程（ステップＳ１０，Ｓ１１）と同様である。
【００９２】
〔第３実施形態〕
次に、本発明の第３実施形態について説明する。第３実施形態は、第１，２実施形態の変形例である。
第１実施形態で例示した［第２線形予測係数算出過程（ステップS３１）］では、第２線形予測係数算出部１１０が、チャネル毎に生成した自己相関係数c1(i)…cM(i)をチャネル間で平均した平均自己相関係数c(i)を用い、短時間線形予測モデルの各線形予測係数b(1)…b(q)を算出していた。
【００９３】
これに対し、第３実施形態の［第２線形予測係数算出過程（ステップS３１）］では、第２線形予測係数算出部４１０が、Ｍ個のセンサのうち、音響信号の音源に最も近い１つのセンサで観測された音響信号を複数の時点でサンプリングして得られた離散音響信号値の自己相関係数を算出し、当該自己相関係数を用い、短時間線形予測モデルの各線形予測係数を算出する。
この点が第１実施形態との相違点である。そして、この構成は第２実施形態へも適用可能である。以下では、第１，２実施形態との相違点である第２線形予測係数算出部４１０の構成及び［第２線形予測係数算出過程（ステップS３１）］の処理のみを説明し、第１，２実施形態と共通する事項については説明を省略する。
【００９４】
図１２Ａは、本形態の第２線形予測係数算出部４１０の機能構成を示したブロック図である。なお、図１２Ａにおいて、第１実施形態と共通する部分については第１実施形態と同じ符号を用いた。また、図１２Ｂは、本形態の［第２線形予測係数算出過程（ステップＳ３１）］を説明するためのフローチャートである。
図１２Ａに例示するように、本形態の第２線形予測係数算出部４１０は、自己相関係数算出部４１１と方程式演算部１１３とを有する。本形態の第２線形予測係数算出過程では、まず、自己相関係数算出部４１１（図１２Ａ）が、入力された離散音響信号値x1(n)…xM(n)を用い、Ｍ（Ｍ≧２）個のセンサのうち音響信号の音源に最も近い１つのセンサｙ（y=1,...,M）で観測された音響信号を複数の時点でサンプリングして得られた離散音響信号値xy(n)の自己相関係数cy (i)(i=0,1,...,q)を算出する（ステップＳ１４１）。なお、音響信号の音源に最も近い１つのセンサｙの情報は、自己相関係数算出部４１１が具備する固定情報であってもよいし、自己相関係数算出部４１１に与えられる変動情報であってもよい。
【００９５】
次に、上述のように求められた各自己相関係数c(i)が方程式演算部１１３に入力され、方程式演算部１１３は、各自己相関係数c(i)を用い、Yule-Walkerの方程式（正規方程式）の解を求めることにより、短時間線形予測モデルの各線形予測係数b(1)…b(q)を算出して出力する（方程式演算過程／ステップＳ１４２）。
以上のように、本形態では、音響信号の音源に最も近い１つのセンサに対応する音響信号値の自己相関係数を用い、各線形予測係数b(1)…b(q)を算出する構成とした。これにより、他のセンサに対応する音響信号値の自己相関係数を用いる場合に比べて線形予測係数b(1)…b(q)の算出精度が向上し、離散音響信号値x1(n)…xM(n)が具備するd(z)に従った短期的な自己相関をより効果的に抑制することができる。前述のように、これは、後部残響除去の精度向上につながる。
【００９６】
〔第４実施形態〕
次に、本発明の第４実施形態について説明する。第４実施形態は、第１，２実施形態の変形例である。
第１実施形態の［擬似白色化過程（ステップＳ２１）］では、短時間線形予測モデルを用いて離散音響信号値の擬似白色化を行った。
これに対し、第４実施形態の［擬似白色化過程（ステップＳ２１）］では、Cepstral Mean Subtraction(CMS)（例えば、「B. S. Atal, "Effectiveness of linear prediction characteristics of the speech wave for automatic speaker identification and verification," Journal of Acoustical Society of America, 55(6), pp. 1304-1312, 1974.」参照）を用いて離散音響信号値の擬似白色化を行う。
【００９７】
この点が第１実施形態との相違点である。そして、この構成は第２実施形態へも適用可能である。以下では、第１，２実施形態との相違点である擬似白色化部５１０の構成及び［擬似白色化過程（ステップＳ２１）］の処理のみを説明し、第１，２実施形態と共通する事項については説明を省略する。
図１３は、本形態のモデル適用部５００の機能構成を示したブロック図である。なお、図１３において第１実施形態と共通する部分については、第１実施形態と同じ符号を用いた。
【００９８】
図１３に例示するように、本形態のモデル適用部５００は、擬似白色化部５１０と第１線形予測係数算出部２００とを有する。また、擬似白色化部５１０は、周波数領域変換部５１１と、時間平均化部５１２と、減算部５１３と、時間領域変換部５１４とを有する。
図１４は、本形態の［擬似白色化過程（ステップＳ２１）］を説明するためのフローチャートである。以下、この図を用いて、本形態の［擬似白色化過程（ステップＳ２１）］を説明する。
まず、擬似白色化部５１０の周波数領域変換部５１１が、メモリ１０ａから音響信号１分析フレーム分のＭチャネルの離散音響信号値x1(n)…xM(n)を読み込む。そして、周波数領域変換部５１１は、短時間フーリエ変換等によって離散音響信号値x1(n)…xM(n)を周波数領域の離散音響信号値X1(f,t)…XM(f,t)に変換して出力する（ステップＳ２０１）。なお、短時間フーリエ変換によってこの処理を行う場合は、例えば以下の式（４８）を用いる。また、Ｆ［・］は短時間フーリエ変換関数を示し、Ｌｏｇ[・］は対数関数を示す。
【００９９】
Xm(f, t)=Log[F[xm(n)]]        …(48)
次に、時間平均化部５１２に周波数領域の離散音響信号値X1(f,t)…XM(f,t)が読み込まれ、時間平均化部５１２は、以下の式（４９）によって、周波数領域の離散音響信号値X1(f,t)…XM(f,t)の時間平均Xm’(f)を求め、出力する（ステップＳ２０２）。
【０１００】
【数２１】
0004774100.tif 000022
次に、減算部５１３に周波数領域の離散音響信号値X1(f,t)…XM(f,t)とそれらの時間平均E{Xm(f,t)}とが読み込まれ、減算部５１３は、以下の式（５０）によってXm’(f,t)(m=1,...,M)を算出し、出力する（ステップＳ２０３）。
Xm’(f,t)=Xm(f,t) −E{Xm(ｆ,t)}         …(50)
次に、時間領域変換部５１４にX1’(f,t)…XM’(f,t)が読み込まれ、時間領域変換部５１４は、逆フーリエ変換等によってこれらを時間領域に変換し、擬似白色化された離散音響信号値x1’(n)…xM’(n)を算出し、出力する（ステップＳ２０４）。なお、逆フーリエ変換によってこの処理を行う場合は、例えば以下の式（５１）を用いる。また、ｉｎｖＦ［・］は逆フーリエ変換関数を示し、ｅｘｐ[・］はネイピア数を底とした指数関数を表す。
【０１０１】
xm’(n)=invF[exp[Xm’(f,t)]]         …(51)
なお、上述した短時間フーリエ変換関数Ｆ［・］や逆フーリエ変換関数ｉｎｖＦ［・］において窓長２５ｍｓの窓関数を用いた場合、２５ｍｓ以内の初期反射成分及び短時間相関を取り除くことができる。
また、本形態を第２実施形態に適用する場合には、時間領域変換部５１４で生成された擬似白色化後の離散音響信号値x1’(n)…xM’(n)は、後部残響予測部３１０ｃや周波数領域変換部３１０ｄ（図８）にも転送される。
【０１０２】
〔第５実施形態〕
次に、本発明の第５実施形態について説明する。第５実施形態は、第４実施形態の擬似白色化手法を第１実施形態に適用する際の変形例である。
第１実施形態の［周波数領域変換過程（ステップＳ３）］では、離散音響信号値x1(n)…xM(n)を周波数領域の離散音響信号値X1(f,t)…XM(f,t)に変換し、後部残響予測値r1(n)…rM(n)を周波数領域の後部残響予測値R1(f,t)…RM(f,t)に変換していた。しかし、第４実施形態の擬似白色化を行う場合、その過程で（ステップＳ２０１）で周波数領域の離散音響信号値X1(f,t)…XM(f,t)が得られている。
第５実施形態では、第４実施形態の擬似白色化の過程で得られた周波数領域の離散音響信号値X1(f,t)…XM(f,t)を流用し、周波数領域変換過程の処理を簡略化する。
以下では、これまで説明した実施形態との相違点を中心に説明し、それらと共通する部分については説明を省略する。
【０１０３】
＜ハードウェア構成＞
第１実施形態で説明したのと同様である。
＜ハードウェアとプログラムとの協働＞
本形態の残響除去装置もコンピュータに所定のプログラムが読み込まれて実行されることにより構成される。図１５は、本形態の残響除去装置６１０の機能構成を例示したブロック図である。なお、図１５において、これまで説明した実施形態と共通する部分についてはそれらと同じ符号を用いた。
【０１０４】
図１５に例示するように、残響除去装置６１０は、メモリ１０ａと、モデル適用部５００と、後部残響予測部１０ｃと、周波数領域変換部５１０ｄと、後部残響除去部１０ｅと、複素スペクトル生成部１０ｆと、時間領域変換部１０ｇと、遅延量算出部１０ｈと、遅延調節部１０ｉと、メモリ１０ｊと、制御部１０ｋとを有する。
＜残響除去処理＞
次に、本形態の残響除去処理について説明する。
図１６は、本形態の残響除去処理の全体を説明するためのフローチャートである。以下、この図を用い、本形態の残響除去処理を説明する。
【０１０５】
［前処理］
第１実施形態と同様である。
［モデル適用過程（ステップＳ２１１）］
モデル適用過程では、モデル適用部５００が、メモリ１０ａから読み込んだ１分析フレーム分のＭチャネルの離散音響信号値x1(n)…xM(n)を用い、式（１３）に示したマルチステップ線形予測モデルの各線形予測係数αw,1(p)…αw,M(p)を算出する（ステップＳ２１１）。この処理うち、擬似白色化処理は第４実施形態で説明した通りであり、その他の処理は第１実施形態と同様である。
【０１０６】
［後部残響予測過程（ステップＳ２１２）］
モデル適用過程（ステップＳ２１１）の後、後部残響予測部１０ｃに、メモリ１０ａから読み出された離散音響信号値x1(n)…xM(n)と、モデル適用過程（ステップＳ２１１）で算出された各線形予測係数αw,1(p)…αw,M(p)とが入力される。
そして、後部残響予測部１０ｃは、各線形予測係数αw,1(p)…αw,M(p)と離散音響信号値x1(n)…xM(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値rw(n)(w=1,...,M)として算出して出力する（ステップＳ２１２）。
【０１０７】
［周波数領域変換過程（ステップＳ２１３）］
次に、周波数領域変換部５１０ｄに後部残響予測過程（ステップＳ２１２）で算出された後部残響予測値r1(n)…rM(n)が入力される。周波数領域変換部５１０ｄは、入力された後部残響予測値r1(n)…rM(n)を周波数領域の後部残響予測値R1(f,t)…RM(f,t)に変換する（ステップＳ２１３）。周波数領域変換部５１０ｄは、この処理により、周波数領域の後部残響予測値R1(f,t)…RM(f,t)の振幅スペクトル|R1(f,t)|…|RM(f,t)|と位相情報arg[R1(f,t)]…arg[RM(f,t)]とを抽出し、出力する。
【０１０８】
［後部残響除去過程（ステップＳ２１４）］
次に、後部残響除去部１０ｅに、擬似白色化部５１０の周波数領域変換部５１１（図１３）から転送された周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と、周波数領域変換部５１０ｄで生成された周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とが入力される。そして、後部残響除去部１０ｅは、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|として出力する（ステップＳ２１４）。
【０１０９】
［フラグ判定過程（ステップＳ２１５，Ｓ２１６）］
本形態の［フラグ判定過程（ステップＳ２１５，Ｓ２１６）］は、第１実施形態の［フラグ判定過程（ステップＳ５，Ｓ６）］と同様である。
［その他の過程］
その他の過程は、第１実施形態と同様である。ただし、［複素スペクトル生成過程（ステップＳ７）］において、擬似白色化部５１０の周波数領域変換部５１１（図１３）から転送された位相情報arg[X1(f,t)]…arg[XM(f,t)]を用いる点のみが第１実施形態と相違する。
【０１１０】
〔第６実施形態〕
次に、本発明の第６実施形態について説明する。第６実施形態は、第４実施形態の擬似白色化手法を第２実施形態に適用する際の変形例である。
第２実施形態の［周波数領域変換過程（ステップＳ１０３）］では、離散音響信号値x1(n)…xM(n)を周波数領域の離散音響信号値X1(f,t)…XM(f,t)に変換し、後部残響予測値r1(n)…rM(n)を周波数領域の後部残響予測値R1(f,t)…RM(f,t)に変換していた。しかし、第４実施形態の擬似白色化を行う場合、その過程で（ステップＳ２０１）で周波数領域の離散音響信号値X1(f,t)…XM(f,t)が得られている。
【０１１１】
第６実施形態では、第４実施形態の擬似白色化の過程で得られた周波数領域の離散音響信号値X1(f,t)…XM(f,t)を流用し、周波数領域変換過程の処理を簡略化する。
以下では、これまで説明した実施形態との相違点を中心に説明し、それらと共通する部分については説明を省略する。
＜ハードウェア構成＞
第１実施形態で説明したのと同様である。
＜ハードウェアとプログラムとの協働＞
本形態の残響除去装置もコンピュータに所定のプログラムが読み込まれて実行されることにより構成される。
【０１１２】
図１７は、本形態の残響除去装置６２０の機能構成を例示したブロック図である。図１７において、これまで説明した実施形態と共通する部分についてはそれらと同じ符号を用いた。
図１７に例示するように、残響除去装置６２０は、メモリ１０ａと、モデル適用部５００と、後部残響予測部３１０ｃと、周波数領域変換部５１０ｄと、後部残響除去部３１０ｅと、複素スペクトル生成部３１０ｆと、時間領域変換部１０ｇと、遅延量算出部１０ｈと、遅延調節部１０ｉと、メモリ１０ｊと、制御部１０ｋとを有する。
【０１１３】
＜残響除去処理＞
次に、本形態の残響除去処理について説明する。
図１８は、本形態の残響除去処理の全体を説明するためのフローチャートである。以下、この図を用い、本形態の残響除去処理を説明する。
［前処理］
第１実施形態と同様である。
［モデル適用過程（ステップＳ２２１）］
モデル適用過程では、モデル適用部５００が、メモリ１０ａから読み込んだ１分析フレーム分のＭチャネルの離散音響信号値x1(n)…xM(n)を用い、式（１３）に示したマルチステップ線形予測モデルの各線形予測係数αw,1(p)…αw,M(p)を算出する（ステップＳ２２１）。この処理うち、擬似白色化処理は第４実施形態で説明した通りであり、その他の処理は第１実施形態と同様である。
【０１１４】
［後部残響予測過程（ステップＳ２２２）］
モデル適用過程（ステップＳ２２１）の後、後部残響予測部３１０ｃに、モデル適用過程（ステップＳ２２１）で擬似白色化された離散音響信号値x1’(n)…xM’(n)と、モデル適用過程（ステップＳ２２１）で算出された各線形予測係数αw,1(p)…αw,M(p)とが入力される。
そして、後部残響予測部３１０ｃは、各線形予測係数αw,1(p)…αw,M(p)と擬似白色化された離散音響信号値x1’(n)…xM’(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値rw(n)(w=1,...,M)として算出して出力する（ステップＳ２２２）。
【０１１５】
［周波数領域変換過程（ステップＳ２２３）］
次に、周波数領域変換部５１０ｄに後部残響予測過程（ステップＳ２２２）で算出された後部残響予測値r1(n)…rM(n)が入力される。周波数領域変換部５１０ｄは、入力された後部残響予測値r1(n)…rM(n)を周波数領域の後部残響予測値R1(f,t)…RM(f,t)に変換する（ステップＳ２２３）。周波数領域変換部５１０ｄは、この処理により、周波数領域の後部残響予測値R1(f,t)…RM(f,t)の振幅スペクトル|R1(f,t)|…|RM(f,t)|と位相情報arg[R1(f,t)]…arg[RM(f,t)]とを抽出し、出力する。
【０１１６】
［後部残響除去過程（ステップＳ２２４）］
次に、後部残響除去部３１０ｅに、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とが入力される。そして、後部残響除去部３１０ｅは、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|…|RM(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|…|SM(f,t)|として出力する（ステップＳ２２４）。なお、本ステップで使用される周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|…|XM’(f,t)|は、擬似白色化部５１０の減算部５１３（図１３）から転送されたものである。
【０１１７】
［フラグ判定過程（ステップＳ２２５，Ｓ２２６）］
本形態の［フラグ判定過程（ステップＳ２２５，Ｓ２２６）］は、第１実施形態の［フラグ判定過程（ステップＳ５，Ｓ６）］と同様である。
［その他の過程］
その他の過程は、第１実施形態と同様である。ただし、［複素スペクトル生成過程（ステップＳ７）］において、擬似白色化部５１０の減算部５１３（図１３）から転送された位相情報arg[X1’(f,t)]…arg[XM’(f,t)]を用いる点のみが第１実施形態と相違する。
【０１１８】
〔第７実施形態〕
次に、本発明の第７実施形態について説明する。第７実施形態は、Ｍ＝１とし、遅延量算出部１０ｈや遅延調節部１０ｉを不要とした第１〜６実施形態の変形例である。その代表例として、Ｍ＝１とし、第２実施形態に第４実施形態の擬似白色化方法を適用し、遅延量算出部１０ｈや遅延調節部１０ｉが存在しない構成について説明する。しかし、その他第１〜６実施形態又はそれらの組合せにおいてＭ＝１とし、遅延量算出部１０ｈや遅延調節部１０ｉが存在しない構成としてもよい。さらに、遅延量算出部１０ｈや遅延調節部１０ｉは存在するが、Ｍ＝１の場合には、それらを機能させない構成であってもよい。
【０１１９】
また、以下では、これまで説明した実施形態との相違点を中心に説明し、それらと共通する事項については説明を省略する。
＜ハードウェア構成＞
第１実施形態と同様である。
＜ハードウェアとプログラムとの協働＞
本形態の残響除去装置もコンピュータに所定のプログラムが読み込まれて実行されることにより構成される。図１９は、本形態の残響除去装置７１０の機能構成を例示したブロック図である。また、図２０は、図１９のモデル適用部８００の機能構成の詳細を例示したブロック図である。なお、図１９，図２０において、これまで説明した実施形態と共通する部分についてはそれらと同じ符号を用いた。
【０１２０】
図１９に例示するように、残響除去装置７１０は、メモリ１０ａと、モデル適用部８００と、後部残響予測部３１０ｃと、周波数領域変換部３１０ｄと、後部残響除去部３１０ｅと、複素スペクトル生成部３１０ｆと、時間領域変換部１０ｇと、メモリ１０ｊと、制御部１０ｋとを有する。
また、モデル適用部８００は、擬似白色化部８１０と第１線形予測係数算出部２００とを有する。また、擬似白色化部８１０は、周波数領域変換部８１１と、時間平均化部８１２と、減算部８１３と、時間領域変換部８１４とを有する。
【０１２１】
＜残響除去処理＞
次に、本形態の残響除去処理について説明する。
図２１は、本形態の残響除去処理の全体を説明するためのフローチャートである。また、図２２Ａは、図２１のステップＳ３０１（モデル適用過程）の詳細を説明するためのフローチャートであり、図２２Ｂは、図２２ＡのステップＳ３１１（擬似白色化過程）の詳細を説明するためのフローチャートである。
以下、これらの図を用い、本形態の残響除去処理を説明する。
【０１２２】
［前処理］
まず、Ｍ（Ｍ＝１）個のセンサで観測された１チャネルの音響信号が所定の標本化周波数でサンプリングされ、離散音響信号値x1(n)が生成される。生成された各チャネルの離散音響信号値x1(n)は、それぞれメモリ１０ａに格納される。なお、本形態では、残響除去を行う全時間区間の離散音響信号値x1(n)を事前に取得し、メモリ１０ａに格納しておき、分析フレーム毎に、以下の各過程を実行する。しかし、離散音響信号値x1(n)の取得をリアルタイムで行いつつ、以下の各過程を実行してもよい。
【０１２３】
また、残響除去装置７１０が最終的に出力する情報が、後部残響除去信号の振幅スペクトルのみであるのか、それとも位相成分をも有する音響信号なのかを示す情報をメモリ１０ｊに格納しておく。本形態では、フラグ（データ）δをメモリ１０ｊに格納しておく。そして、出力する情報が直接音の振幅スペクトルのみである場合δ＝１とし、位相成分をも有する音響信号である場合δ＝０とする。
以下、本形態の後部残響除去処理を説明する。なお、以下では、１つの分析フレームの処理過程のみを説明するが、実際は複数の分析フレームに対して同様な処理が行われる。
【０１２４】
［モデル適用過程（ステップＳ３０１）］
モデル適用過程では、モデル適用部８００が、メモリ１０ａから読み込んだ１分析フレーム分の離散音響信号値x1(n)を用い、式（５）に示したマルチステップ線形予測モデルの各線形予測係数α1,1(p)を算出する（ステップＳ３０１）。以下、この処理の詳細を階層的に説明する。
［モデル適用過程（ステップＳ３０１）の詳細（図２２Ａ）］
図２２Ａに例示するように、モデル適用過程では、まず、擬似白色化（Pre-whitening）部８１０（図２０）が、入力された離散音響信号値x1(n)が有する短時間区間での自己相関成分を抑制し、擬似白色化した離散音響信号値x1’を生成して出力する（擬似白色化過程／ステップＳ３１１）。すなわち、各離散時間の上記離散音響信号値x1(n)から、当該離散時間n直前の短時間区間内の各離散音響信号値と自己相関性を持つ自己相関成分を抑制し、擬似白色化した離散音響信号値x1’(n)を生成する。
【０１２５】
次に、上述の擬似白色化された離散音響信号値x1’(n)が第１線形予測係数算出部２００（図２０）に入力され、第１線形予測係数算出部２００は、当該擬似白色化された離散音響信号値x1’(n)を用い、式（５）に示したマルチステップ線形予測モデルの各線形予測係数α1,1(p)を算出して出力する（第１線形予測係数算出過程／ステップＳ３１２）。なお、一例として、式（５）における遅延Ｄを、例えば２５ｍｓ（標本化周波数周波数１２０００Ｈｚの場合、３００タップに相当）とし、各線形予測係数α1,1(p)の数Ｎは、例えば、５０００程度とする。また、各線形予測係数α1,1(p)を算出する方法としては、自己相関法（correlation method）や共分散法（covariance method）を例示できる。また、MATLAB（登録商標）等を利用してこの処理を行ってもよい。
【０１２６】
［擬似白色化過程（ステップＳ３１１）の詳細（図２２Ｂ）］
次に、擬似白色化過程（ステップＳ３１１）の詳細を説明する。
本形態では、一例として、Cepstral Mean Subtraction(CMS)を用いて離散音響信号値の擬似白色化を行う。
まず、擬似白色化部８１０の周波数領域変換部８１１が、メモリ１０ａから音響信号１分析フレーム分の１チャネルの離散音響信号値x1(n)を読み込む。そして、周波数領域変換部８１１は、短時間フーリエ変換等によって離散音響信号値x1(n)を周波数領域の離散音響信号値X1(f,t)に変換して出力する（ステップＳ３２１）。なお、短時間フーリエ変換によってこの処理を行う場合は、例えば以下の式（５２）を用いる。また、Ｆ［・］は短時間フーリエ変換関数を示し、Ｌｏｇ[・］は対数関数を示す。
【０１２７】
X1(f, t)=Log[F[x1(n)]]        …(52)
次に、時間平均化部８１２に周波数領域の離散音響信号値X1(f,t)が読み込まれ、時間平均化部８１２は、以下の式（５３）によって、周波数領域の離散音響信号値X1(f,t)の時間平均X1’(f)を求め、出力する（ステップＳ３２２）。
【０１２８】
【数２２】
0004774100.tif 000023
次に、減算部８１３に周波数領域の離散音響信号値X1(f,t)とその時間平均E{X1(f,t)}とが読み込まれ、減算部８１３は、以下の式（５４）によってX1’(f,t)を算出し、出力する（ステップＳ３２３）。
X1’(f,t)=X1(f,t) −E{X1(ｆ,t)}         …(54)
次に、時間領域変換部８１４にX1’(f,t)が読み込まれ、時間領域変換部８１４は、逆フーリエ変換等によってこれらを時間領域に変換し、擬似白色化された離散音響信号値x1’(n)を算出し、出力する（ステップＳ３２４）。なお、逆フーリエ変換によってこの処理を行う場合は、例えば以下の式(５５)を用いる。また、ｉｎｖＦ［・］は逆フーリエ変換関数を示す。
【０１２９】
x1’(n)=invF[exp[X1’(f,t)]]      …(55)
なお、上述した短時間フーリエ変換関数Ｆ［・］や逆フーリエ変換関数ｉｎｖＦ［・］において窓長２５ｍｓの窓関数を用いた場合、２５ｍｓ以内の初期反射成分及び短時間相関を取り除くことができる。
また、本形態の例では、時間領域変換部８１４で生成された擬似白色化後の離散音響信号値x1’(n)は、後部残響予測部３１０ｃや周波数領域変換部３１０ｄ（図１９）にも転送される（［モデル適用過程（ステップＳ３０１）の詳細］の説明終わり）。
【０１３０】
［後部残響予測過程（ステップＳ３０２）］
モデル適用過程（ステップＳ３０１）の後、後部残響予測部３１０ｃに、時間領域変換部８１４で生成された擬似白色化後の離散音響信号値x1’(n)と、モデル適用過程（ステップＳ３０１）で算出された各線形予測係数α1,1(p)とが入力される。
そして、後部残響予測部３１０ｃは、前述の式（１０）のように、各線形予測係数α1,1(p)と擬似白色化後の離散音響信号値x1’(n)とをマルチステップ線形予測モデルの線形予測項に代入して得られた線形予測値を、後部残響予測値r1(n)として算出して出力する（ステップＳ３０２）。本形態では、式（５）のマルチステップ線形予測モデルを用いているため、後部残響予測部３１０ｃは、以下の式（５６）に従って後部残響予測値r1(n)を求めて出力する。
【０１３１】
【数２３】
0004774100.tif 000024
［周波数領域変換過程（ステップＳ３０３）］
次に、周波数領域変換部３１０ｄに、時間領域変換部８１４（図２０）で生成された擬似白色化後の離散音響信号値x1’(n)と、後部残響予測過程（ステップＳ３０２）で算出された後部残響予測値r1(n)とが入力される。周波数領域変換部３１０ｄは、入力された擬似白色化後の離散音響信号値x1’(n)を周波数領域の離散音響信号値X1’(f,t)に変換し、後部残響予測値r1(n)を周波数領域の後部残響予測値R1(f,t)に変換する（ステップＳ３０３）。本形態では、例えば、窓長２５ｍｓのハニング窓などの有限長の窓関数を用い、短時間フーリエ変換（DFT: Discrete Fourier Transform）等によって、これらの周波数領域への変換を行う。周波数領域変換部３１０ｄは、これらの処理により、周波数領域の離散音響信号値X1’(f,t)の振幅スペクトル|X1’(f,t)|と位相情報arg[X1’(f,t)]、及び、周波数領域の後部残響予測値R1(f,t)の振幅スペクトル|R1(f,t)|と位相情報arg[R1(f,t)]とを抽出し、出力する。
【０１３２】
［後部残響除去過程（ステップＳ３０４）］
次に、後部残響除去部３１０ｅに、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|とが入力される。そして、後部残響除去部３１０ｅは、周波数領域の離散音響信号値の振幅スペクトル|X1’(f,t)|と、周波数領域の後部残響予測値の振幅スペクトル|R1(f,t)|とのセンサ毎の相対値を求め、当該相対値を後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|として出力する（ステップＳ３０４）。この処理の詳細は第１実施形態と同様である。
【０１３３】
［フラグ判定過程（ステップＳ３０５，Ｓ３０６）］
ステップＳ３０４の後、制御部１０ｋが、メモリ１０ｊに格納されているフラグδを読み出し、そのフラグδが、振幅スペクトルのみを出力することを示すフラグであるか否か、すなわち、δ＝１であるか否かを判断する（ステップＳ３０５）。ここで、δ＝１であれば、制御部１０ｋは、後部残響除去過程（ステップＳ３０４）で後部残響除去部３１０ｅが生成した後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|を、残響除去装置７１０の最終的な出力情報として出力し（ステップＳ３０６）、当該分析フレームの処理を終了させる。このように出力された振幅スペクトル予測値|S1(f,t)|は、例えば、残響除去装置７１０の後段に続く音声認識システム等のアプリケーションに渡され、特徴量に変換される。
一方、δ＝０であれば、制御部１０ｋは、以下のステップＳ３０７以降の処理を実行させる。
【０１３４】
［複素スペクトル生成過程（ステップＳ３０７）］
複素スペクトル生成過程では、まず、複素スペクトル生成部３１０ｆに、後部残響除去部３１０ｅから出力（ステップＳ３０４）された後部残響除去信号値の振幅スペクトル予測値|S1(f,t)|と、周波数領域変換部３１０ｄから出力（ステップＳ３０３）された周波数領域の離散音響信号値の位相情報arg[X1'(f,t)]とが入力される。複素スペクトル生成部３１０ｆは、これらの情報を用い、以下の式（５７）に従って、後部残響除去信号値の複素スペクトル予測値S1(f,t)を算出して出力する（ステップＳ３０７）。
S1(f,t)=|S1(f,t)|・exp(j・arg[X1(f,t)])    …(57)
【０１３５】
［時間領域変換過程（ステップＳ３０８）］
ステップＳ３０７の後、時間領域変換部１０ｇに、上述の後部残響除去信号値の複素スペクトル予測値S1(f,t)が入力される。そして、時間領域変換部１０ｇは、後部残響除去信号値の複素スペクトル予測値S1(f,t)を時間領域に変換した後部残響除去信号推定値s1(n)を算出して出力する（ステップＳ３０８）。なお、時間領域への変換は、例えば、逆フーリエ変換によって行う。
【０１３６】
〔シミュレーション結果〕
次に、Ｍ＝１の場合における本発明の効果を示すためのシミュレーション結果を示す。ここでは、第２実施形態に第４実施形態の擬似白色化方法を適用した構成でシミュレーションを行った。
このシミュレーションでは、連続発話データセットから女性と男性のそれぞれ５０発話を取り出し、３０００タップのインパルス応答と畳み込み残響環境をシミュレートした。また、式（５）のマルチステップ線形予測モデルのステップサイズ（遅延）Ｄを２５ｍｓとし、線形予測係数α1,1(p)の数Ｎを５０００とした。また、時間領域から周波数領域への変換には、窓長２５ｍｓの短時間フーリエ変換を用いた。
【０１３７】
図２４にこのシミュレーション結果を示す。ここで、図２４Ａ、図２４Ｂは、それぞれ、残響除去前の振幅スペクトラム値及び音声波形を示した図である。また、図２４Ｃ、図２４Ｄは、それぞれ、本発明（Ｍ＝１）による残響除去後の振幅スペクトラム値及び音声波形を示した図である。なお、図２４Ａ、図２４Ｃの縦軸は振幅スペクトラム値を示し、横軸は時間（ｓ）を示す。また、図２４Ｂ、図２４Ｄの縦軸は周波数（Ｈｚ）を示し、横軸は時間（ｓ）を示す。これらの図からも、本発明によって後部残響が精度良く抑圧されることがわかる。
【０１３８】
次に、本発明の効果を音声認識の観点から評価したシミュレーション結果を示す。
このシミュレーションでは、クリーン音声を用いて構築された音響モデルを用いた。表１に、それぞれの認識対象の単語誤り率を示す。残響音声、残響除去音声の単語誤り率は、それぞれ「Rev.」と「Derev.」と表されている。音響モデルがクリーン音声から学習されたにもかかわらず、本発明により認識率が大幅に改善されていることがわかる。
【０１３９】
【表１】
0004774100.tif 000025
〔実験結果〕
次に、本発明の効果を示すための実験結果を示す。この実験は、後部残響除去を行わない場合（処理無）、Ｍ＝１とし、第２実施形態に第４実施形態の擬似白色化方法を適用した方法（第７実施形態）、第１実施形態（Ｍ≧２）で遅延調節を行うことなく１つのチャネル（m=1）で得られた後部残響除去信号推定値を用いた場合（第１実施形態（遅延調節無））、及び、第１実施形態（Ｍ≧２）で遅延調節を行って後部残響除去を行った場合（第１実施形態（遅延調節無））について、各々の音声認識率を測定した。
【０１４０】
図２５Ａは、この実験条件を示す図である。この実験では、縦３．５ｍ、横４．５ｍ、高さ２．５ｍの室内に、４つのマイクロフォン１０１０（Ｍ＝４）を一列に配置し、ｍ＝１のマイクロフォン１０１０（実線）から、０、５ｍ,１．０ｍ,１．５ｍ,２．０ｍの距離に４つのスピーカ１０２０を一直線に配置した場合を想定した。また、連続発話データセットから女性と男性のそれぞれ１００発話を取り出し、これらに、シミュレートした３０００タップのインパルス応答と畳み込んで後部残響音声を作成した。また、音声認識の際には、音響モデル適応処理としてCepstral Mean Subtraction（B.S. Atal, "Effectiveness of linear prediction characteristics of the speech wave for automatic speaker identification and verification," Journal of the Acoustical Society of America, vol. 55(6), pp. 1304-1312, Jun. 1974. )を用いた。なお、Ｍ＝１とし、第２実施形態に第４実施形態の擬似白色化方法を適用した方法（第７実施形態）では、ｍ＝１のマイクロフォン１０１０を使用した。
【０１４１】
図２５Ｂは、上記の４つの場合〔処理無、第７実施形態、第１実施形態（遅延調節無）、第１実施形態（遅延調節有）〕についての音声認識結果（単語誤り率）を示すグラフである。なお、図２５Ｂでは、マイクロフォン１０１０（ｍ＝１）と各スピーカ１０２０との距離（ｍ）を横軸とし、単語誤り率（％）を縦軸としている。
図２５Ｂに例示するように、マイクロフォン１０１０（ｍ＝１）と各スピーカ１０２０との距離が比較的近い場合には、第７実施形態（Ｍ＝１）に対する第１実施形態（Ｍ≧２）（遅延調節無）及び第１実施形態（Ｍ≧２）（遅延調節有）の単語誤り率の改善量はあまり大きくない。しかし、マイクロフォン１０１０（ｍ＝１）と各スピーカ１０２０との距離が離れるにつれ、伝達関数中の最大位相成分（ゼロ点）が増えるため、第７実施形態（Ｍ＝１）に対する第１実施形態（Ｍ≧２）（遅延調節無）の単語誤り率の改善量は顕著になっていく。さらに、第１実施形態（Ｍ≧２）（遅延調節有）場合には、より一層単語誤り率を改善することができる。
【０１４２】
〔変形例等〕
なお、本発明は上述の各実施形態に限定されるものではない。例えば、各実施形態では、後部残響除去部が、短時間フーリエ変換等により各データを周波数領域に変換して各処理を実行した。したし、残響除去装置の出力として要求される信号が直接音の振幅スペクトルのみであるならば、後部残響除去部が、各データをｚ変換し、ｚ領域で各処理を実行してもよい。
また、各実施形態では、擬似白色化部によって離散音響信号値から短時間相関を取り除いた後、各処理を実行した。しかし、短時間相関を取り除いていない離散音響信号値を用いて各処理を実行してもよい。
【０１４３】
また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。また、２以上の実施形態を結合した形態であってもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
【０１４４】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１４５】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、  このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１４６】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【０１４７】
本発明を、各種音響信号処理システムの要素技術として用いることで、そのシステム全体の性能を向上させることができる。本発明が適用可能な音響信号処理システムとしては、例えば、以下のようなものを列挙できる。環境で収録された音声には、常に残響（反射音）が含まれるが、以下にあげるシステムは、そのような状況で用いられることを想定した例である。
・残響環境での音声認識システム
・歌われたり、楽器で演奏されたり、スピーカで演奏された楽曲の残響を除去してメモリ格納しておき、それら楽曲を検索したり、採譜したりする音楽情報処理システム
・人が発した音に反応して機械にコマンドを渡す機械制御インターフェース、及び機械と人間との対話装置
・残響環境下で残響を除去することで聞き取り易さを向上させる補聴器
・残響除去により音声の明瞭度を向上させるＴＶ会議システムなどの通信システム
【図面の簡単な説明】
【０１４８】
【図１】図１は、第１実施形態の残響除去装置の機能構成を例示したブロック図である。
【図２】図２Ａは、モデル適用部の機能構成の詳細を例示したブロック図である。図２Ｂは、遅延調節部の機能構成の詳細を例示したブロック図である。
【図３】図３は、第１実施形態における残響除去装置１０のハードウェア構成を例示したブロック図である。
【図４】図４は、第１実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図５】図５、第１実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図６】図６Ａは、図４のステップＳ１（モデル適用過程）の詳細を説明するためのフローチャートである。図６Ｂは、図６ＡのステップＳ２１（擬似白色化過程）の詳細を説明するためのフローチャートである。
【図７】図７Ａは、図６ＢのステップＳ３１（第２線形予測係数算出過程）の詳細を説明するためのフローチャートである。図７Ｂは、図４のステップＳ４の詳細を説明するためのフローチャートである。
【図８】図８は、第２実施形態の残響除去装置の機能構成を例示したブロック図である。
【図９】図９は、モデル適用部の機能構成の詳細を例示したブロック図である。
【図１０】図１０は、第２実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図１１】図１１は、第２実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図１２】図１２Ａは、第３実施形態の第２線形予測係数算出部の機能構成を示したブロック図である。図１２Ｂは、第３実施形態の［第２線形予測係数算出過程（ステップＳ３１）］を説明するためのフローチャートである。
【図１３】図１３は、第４実施形態のモデル適用部の機能構成を示したブロック図である。
【図１４】図１４は、第４実施形態の［擬似白色化過程（ステップＳ２１）］を説明するためのフローチャートである。
【図１５】図１５は、第５実施形態の残響除去装置の機能構成を例示したブロック図である。
【図１６】図１６は、第５実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図１７】図１７は、第６実施形態の残響除去装置の機能構成を例示したブロック図である。
【図１８】図１８は、第６実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図１９】図１９は、第７実施形態の残響除去装置の機能構成を例示したブロック図である。
【図２０】図２０は、図１９のモデル適用部の機能構成の詳細を例示したブロック図である。
【図２１】図２１は、第７実施形態の残響除去処理の全体を説明するためのフローチャートである。
【図２２】図２２Ａは、図２１のステップＳ３０１（モデル適用過程）の詳細を説明するためのフローチャートである。図２２Ｂは、図２２ＡのステップＳ３１１（擬似白色化過程）の詳細を説明するためのフローチャートである。
【図２３】図２３Ａは室内伝達関数値hを縦軸にとり時間（ｍｓ）を横軸にとったグラフである。図２３Ｂは合成伝達関数値gを縦軸にとり時間（ｍｓ）を横軸にとったグラフである。また、図２３Ｃは室内伝達関数hと合成伝達関数gとのエネルギー差を縦軸にとり時間（ｍｓ）を横軸にとったグラフである。
【図２４】図２４Ａ、図２４Ｂは、それぞれ、残響除去前の振幅スペクトラム値及び音声波形を示した図である。また、図２４Ｃ、図２４Ｄは、それぞれ、本発明（Ｍ＝１の場合）による残響除去後の振幅スペクトラム値及び音声波形を示した図である。
【図２５】図２５Ａは、実験条件を示す図である。図２５Ｂは、音声認識結果（単語誤り率）を示すグラフである。
【符号の説明】
【０１４９】
１０，３１０，６１０，６２０，７１０  残響除去装置
【図１】
0004774100.tif 000028
【図２】
0004774100.tif 000029
【図３】
0004774100.tif 000030
【図４】
0004774100.tif 000031
【図５】
0004774100.tif 000032
【図６】
0004774100.tif 000033
【図７】
0004774100.tif 000034
【図８】
0004774100.tif 000035
【図９】
0004774100.tif 000036
【図１０】
0004774100.tif 000037
【図１１】
0004774100.tif 000038
【図１２】
0004774100.tif 000039
【図１３】
0004774100.tif 000040
【図１４】
0004774100.tif 000041
【図１５】
0004774100.tif 000042
【図１６】
0004774100.tif 000043
【図１７】
0004774100.tif 000044
【図１８】
0004774100.tif 000045
【図１９】
0004774100.tif 000046
【図２０】
0004774100.tif 000047
【図２１】
0004774100.tif 000048
【図２２】
0004774100.tif 000049
【図２３】
0004774100.tif 000050
【図２４】
0004774100.tif 000051
【図２５】
0004774100.tif 000052