JPWO2007094463A1 - 信号歪み除去装置、方法、プログラム及びそのプログラムを記録した記録媒体 - Google Patents

Abstract

所定の繰り返し終了条件を満たした場合には、観測信号に逆フィルタを適用した結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、観測信号に逆フィルタを適用した結果をアドホック信号として出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、各フレームのアドホック信号に対してこのフレームに対応する予測誤差フィルタを適用して得る各信号（イノベーション推定値）を結合した全フレームでのイノベーション推定値が、その全サンプル間で独立となる逆フィルタを求めてこれを出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段とを備えた信号歪み除去装置とする。

Description

この発明は、信号の歪み（ひずみ：distortion）を除去する技術に関する。

信号は反射や残響などが存在する環境で観測されると、本来の信号に反射や残響などが重畳された信号として観測される。以下、本来の信号を「原信号」と云い、観測された信号を「観測信号」と云うことにする。また、反射や残響などに代表される、原信号に重畳された歪みを「伝達特性」と云うことにする。このため、観測信号から原信号固有の特徴を抽出することが困難になる。この不都合を解消すべく、かねてより信号歪み除去処理技術が各種試みられてきた。信号歪み除去処理は、原信号に重畳した伝達特性を観測信号から取り除く処理である。

従来の信号歪み除去処理の一例として非特許文献１に開示されている信号歪み除去方法を図１５を用いて説明する。予測誤差フィルタ計算部（９０１）は、観測信号をフレーム化処理して、各フレームに含まれる観測信号に対して線形予測分析を行い、予測誤差フィルタを計算する。この明細書では、フィルタはディジタルフィルタであり、信号のサンプル値に対して作用するいわゆるフィルタ係数を求める意味で単にフィルタを計算するなどということがある。予測誤差フィルタ適用部（９０２）は、各フレームごとに上記計算された予測誤差フィルタを当該フレームの観測信号に適用する。逆フィルタ計算部（９０３）は、予測誤差フィルタ適用後の信号に対して逆フィルタを適用して得られる信号の正規化尖度が最大となるような逆フィルタを計算する。逆フィルタ適用部（９０４）は、上記計算された逆フィルタを観測信号に適用することで信号歪み除去後の信号（復元信号）を得る。
B.W.Gillespie, , H.S.Malvar, ,and D.A.F.Florencio, ,"Speech dereverberation via maximum-kurtosis subband adaptive filtering," IEEE International Conference on Acoustics, Speech, and Signal Processing, pp.3701-3704, 2001.

上記の従来的な信号歪み除去方法は、観測信号の各フレーム内ではサンプル間の相関は原信号固有の特性の寄与が大きく、フレームを跨ぐサンプル間の相関は伝達特性による寄与が大きいことを仮定している。上記従来方法は、この仮定に基づいて、フレーム化処理された観測信号に予測誤差フィルタを適用して観測信号中の原信号固有の特性の寄与を低減している。

しかし、この仮定は粗い近似であるため、逆フィルタの精度は不十分である。つまり、観測信号から求まる予測誤差フィルタは伝達特性の影響を受けているので、原信号固有の特性のみを正しく取り除くことができない。このため、予測誤差フィルタ適用後の信号から求める逆フィルタの精度は劣化する。結果として、観測信号に逆フィルタを適用して得る信号は、本来の原信号との比較で精度の良いものではなかった。
そこで本発明は、精度の良い逆フィルタを求めることで、伝達特性に由来する歪みを観測信号から除去して精度の良い復元信号を得ることを目的とする。

上記課題を解決するため、本発明の信号歪み除去装置は、所定の繰り返し終了条件を満たした場合には、観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、観測信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、観測信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段と、を備える。
この発明では、観測信号に伝達特性を除去するための逆フィルタを適用して得るアドホック信号に対して、このアドホック信号に基づいて求めた予測誤差フィルタを適用して得る信号（イノベーション推定値系列）が、その全サンプル間で独立となるような逆フィルタを求める。そして、所定の繰り返し終了条件を満たしたときの逆フィルタを観測信号に適用することで復元信号を得る。

上記の信号歪み除去装置では、予測誤差フィルタ計算手段は、各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの予測誤差フィルタを出力するものであり、逆フィルタ計算手段は、上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタとして求め、この逆フィルタを出力するものであるとしてもよい。
この構成は、イノベーション系列の独立性の尺度として相互情報量を規定し、これを最小化する予測誤差フィルタと逆フィルタを交代変数法で求めるものである。この詳細は後述する。

あるいは、上記の信号歪み除去装置では、予測誤差フィルタ計算手段は、各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの予測誤差フィルタを出力するものであり、逆フィルタ計算手段は、上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタとして求め、この逆フィルタを出力するものであるとしてもよい。
この構成は、イノベーション系列の独立性の尺度として相互情報量を規定し、これを最小化する予測誤差フィルタと逆フィルタを交代変数法で求めるものであるが、信号の高次統計量を用いることなく予測誤差フィルタと逆フィルタを交代変数法で求めることができる。

上記の信号歪み除去装置では、プリ・ホワイトニング処理を前置させ、プリ・ホワイトニング処理で得られた白色化信号に対して、上記同様の処理を行う装置構成とすることができる。具体的には、観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算手段と、白色化フィルタを観測信号に適用して白色化信号を出力する白色化フィルタ適用手段と、所定の繰り返し終了条件を満たした場合には、白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、白色化信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、白色化信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段と、を備えた信号歪み除去装置とすることができる。

上記課題を解決するため、本発明の方法は、逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、観測信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、観測信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、予測誤差フィルタ計算手段が、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、逆フィルタ計算手段が、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、制御手段が、繰り返し終了条件を満たすまで逆フィルタ適用ステップ、予測誤差フィルタ計算ステップ、逆フィルタ計算ステップを繰り返し実行させる制御ステップと、を有する信号歪み除去方法とする。

また、上記の信号歪み除去方法では、プリ・ホワイトニング処理を前置させ、プリ・ホワイトニング処理で得られた白色化信号に対して、上記同様の処理を行う方法とすることができる。具体的には、白色化フィルタ計算手段が、観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算ステップと、白色化フィルタ適用手段が、白色化フィルタを観測信号に適用して白色化信号を出力する白色化フィルタ適用ステップと、逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、白色化信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、白色化信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、予測誤差フィルタ計算手段が、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、逆フィルタ計算手段が、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値系列という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、制御手段が、繰り返し終了条件を満たすまで逆フィルタ適用ステップ、予測誤差フィルタ計算ステップ、逆フィルタ計算ステップを繰り返し実行させる制御ステップと、を有する信号歪み除去方法とする。

本発明の信号歪み除去装置としてコンピュータを機能させる信号歪み除去プログラムによって、コンピュータを信号歪み除去装置として作動処理させることができる。そして、この信号歪み除去プログラムを記録した、コンピュータに読み取り可能なプログラム記録媒体によって、他のコンピュータを信号歪み除去装置として機能させることや、信号歪み除去プログラムを流通させることなどが可能になる。

本発明では、観測信号から求まる予測誤差フィルタを用いて観測信号中の原信号固有の特性の寄与を低減するのではなく、観測信号に（仮の）逆フィルタを適用して得られるアドホック信号（仮の復元信号）から求まる予測誤差フィルタを用いて観測信号中の原信号固有の特性を低減する。アドホック信号から求まる予測誤差フィルタは、伝達特性の影響を受けにくいので、原信号固有の特性をより正確に取り除くことが可能である。このような予測誤差フィルタをアドホック信号に適用して得る信号（イノベーション推定値系列）が全サンプル間で独立となるように求められた逆フィルタは精度良く伝達特性を除去可能なものであるから、このような逆フィルタを観測信号に適用することで、伝達特性に由来する歪みが除去された精度の良い復元信号を得ることができる。

本発明の原理を説明するためのモデル機構を表したブロック線図。 第１実施形態に係る信号歪み除去装置（１）のハードウェア構成例を示す図。 第１実施形態に係る信号歪み除去装置（１）の機能構成例を示す機能ブロック図。 信号歪み除去装置（１）の逆フィルタ計算部（１３）の機能構成例を示す機能ブロック図。 第１実施形態における信号歪み除去処理の流れを示す処理フロー図。 第２実施形態に係る信号歪み除去装置（１）の機能構成例を示す機能ブロック図。 第２実施形態における信号歪み除去処理の流れを示す処理フロー図。 観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ_１とＤ_５０値の関係を示す図。 Ａは残響を含まない音声のスペクトログラム、Ｂは残響を含む音声のスペクトログラム、Ｃは残響除去後の音声のスペクトログラム。 Ａは残響除去音声のＬＰＣスペクトル歪みの時間変動を説明するためのグラフ、Ｂは対応する区間における原音声信号の抜粋。 第３実施形態に係る信号歪み除去装置（１）の逆フィルタ計算部（１３）の機能構成例を示す機能ブロック図。 第３実施形態における信号歪み除去処理の流れを示す処理フロー図。 観測信号長Ｎを３秒、４秒、５秒、１０秒としたときの、ＲＡＳＴＩの値を表示した図。 残響除去前後におけるエネルギー減衰曲線の例を示した図。 従来技術を説明するための機能ブロック図。

§１ 本発明の理論
以下、実施形態の説明に先立ち、本発明の理論を説明する。
以下の説明では、特に断りのない限り、信号源は１つとする。

１．１ 信号
本発明の対象となる信号は、人の音声、音楽、生体信号、測定対象物の物理量をセンサで観測した電気信号などの信号を広く包含する。より好ましくは、自己回帰（Autoregressive：ＡＲ）過程として表現することができる、あるいは表現することが好ましい信号であればよい。例えば音声信号は、通常、区分定常な自己回帰過程として表現される信号、すなわち独立同一分布（i.i.d. : Independent and Identically Distributed）信号に音韻性を表すＡＲ系を作用させた信号として看做される（参考文献１参照）。
以下、信号の代表例として音声信号を挙げて本発明の理論を説明する。
（参考文献１） L.R.Rabiner, R.W.Schafer, "Digital Processing of Speech Signals", Bell Laboratories, Incorporated, 1978.

１．２ 音声信号のモデル化
まず、原信号である音声信号ｓ（ｔ）を、以下の３つの条件を満足する信号としてモデル化する。

［条件１］音声信号ｓ（ｔ）は、区分定常なＡＲ過程で生成される。
この[条件１]から、ＡＲ過程の次数をＰ、定常とみなせる区間長をＷサンプルとして音声信号ｓ（ｔ）をフレーム化すると、第ｉフレームの音声信号ｓ_ｉ（ｎ）は、式（１）のように表される。式（２）は、第ｉフレームの音声信号ｓ_ｉ（ｎ）のサンプルと、フレーム化前の音声信号ｓ（ｔ）のサンプルとの対応を示している。つまり、第ｉフレームのｎ番目のサンプルは、音声信号ｓ（ｔ）において、（ｉ−１）Ｗ＋ｎ番目のサンプルに相当する。式（１）および式（２）において、ｂ_ｉ（ｋ）は線形予測係数、ｅ_ｉ（ｎ）はイノベーションを表す。但し、１≦ｎ≦Ｗ、１≦ｔ≦Ｎ、Ｎは全サンプル数である。以下、特に断りの無い限り、パラメータｎは１フレームのサンプル番号を表し、パラメータｔは全てのサンプル番号を表す。また、全フレーム数はＦとする。

なお、第ｉフレームにおけるｎ番目のイノベーションｅ_ｉ（ｎ）についても、第ｉフレームのｎ番目のイノベーションｅ_ｉ（ｎ）と、フレーム化前の音声信号ｓ（ｔ）に対するイノベーションｅ（ｔ）との対応を示すことができる。この場合、第ｉフレームのｎ番目のイノベーションｅ_ｉ（ｎ）は、イノベーションｅ（ｔ）において（ｉ−１）Ｗ＋ｎ番目のイノベーションに相当し、ｅ_ｉ（ｎ）＝ｅ（（ｉ−１）Ｗ＋ｎ）が成り立つ。

式（１）をｚ変換する。左辺のｚ変換をＳ_ｉ（Ｚ）とし、右辺第二項のｚ変換をＥ_ｉ（Ｚ）とし、Ｂ_ｉ（ｚ）＝Σ_ｋ＝１ ^Ｐｂ_ｉ（ｋ）ｚ^−ｋとすれば、右辺第一項は、Ｂ_ｉ（ｚ）Ｓ_ｉ（Ｚ）となる。従って、式（１）のｚ変換は、（１−Ｂ_ｉ（ｚ））Ｓ_ｉ（Ｚ）＝Ｅ_ｉ（Ｚ）である。なお、ｚ^−１は時間領域では１タップ遅延素子に相当する。以降、時間領域信号（タップ重み係数）を小文字で、ｚ領域信号（伝達関数）を大文字でそれぞれ表す。１−Ｂ_ｉ（ｚ）は最小位相性を満足しなければならず、『１−Ｂ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ』ことが要求される。

［条件２］第ｉフレームに属するイノベーション系列ｅ_ｉ（１），…，ｅ_ｉ（Ｗ）は独立且つ同一分布に属する。イノベーション系列ｅ_ｉ（１），…，ｅ_ｉ（Ｗ）の確率分布の平均及び歪度（３次キュムラント）は０、尖度（４次キュムラント）は正である。さらに、異なるフレームｉ、ｊ〔ｉ≠ｊ〕に属するイノベーションｅ_ｉ（ｎ）とｅ_ｊ（ｎ′）同士も独立である。ただし、これらは必ずしも同一分布に属するとは限らない。
［条件３］予測誤差フィルタ１−Ｂ_ｉ（ｚ）は、相異なるフレーム間で共通する零点をもたない。

式（１）および式（２）から、音声信号ｓ（ｔ）は、式（３）のように表される。[・]はガウス記号である。

このとき、[条件２]は、『イノベーション過程ｅ（ｔ）は時間的に独立な信号である。また、その統計的性質（あるいは統計量）はフレーム内では定常である。』と表現できる。また、[条件３]は、『線形予測係数｛ｂ_ｉ（ｋ）｝_ｋ＝１ ^Ｐは、時不変な極を持たない』と表現できる。

１．３ 観測信号のモデル化
次に、Ｍ個のマイクロホンで音声信号を観測して観測信号を得たときの観測信号をモデル化する。但しＭは、Ｍ≧１の整数である。
ｍ番目（１≦ｍ≦Ｍ）のマイクロホンで観測される残響を含む観測信号ｘ_ｍ（ｔ）を、音源からｍ番目のマイクロホンに至る経路の伝達関数Ｈ_ｍ（ｚ）のタップ重み係数｛ｈ_ｍ（ｋ）；０≦ｋ≦Ｋ；Ｋはインパルス応答の持続時間とする。｝を用いて式（４）のようにモデル化する。ここでは、音声信号の場合の伝達特性の代表例として残響を挙げて、伝達特性を残響に言い換えて説明する。但し、伝達特性を残響に限定する趣旨ではない。

Ｍ個の観測信号についてまとめて表現すれば、式（５）のように表すことができる。但し、式（５）において、ｘ（ｔ）＝［ｘ_１（ｔ），…，ｘ_Ｍ（ｔ）］^Ｔ、ｈ（ｋ）＝［ｈ_１（ｋ），…，ｈ_Ｍ（ｋ）］^Ｔである。

１．４ 信号歪み除去の原理
信号歪み除去後の復元信号ｙ（ｔ）は、多チャネル逆フィルタ｛Ｇ_ｍ（ｚ）；１≦ｍ≦Ｍ｝のタップ重み係数｛ｇ_ｍ（ｋ）；１≦ｍ≦Ｍ，０≦ｋ≦Ｌ；Ｌは逆フィルタの次数｝を用いて式（６）により計算される。本発明においては、逆フィルタ係数であるｇ_ｍ（ｋ）を観測信号ｘ_１（ｔ），…，ｘ_Ｍ（ｔ）のみから推定する。

１．５ 本発明の基本原理
本発明の基本原理は、伝達関数｛Ｈ_ｍ（ｚ）；１≦ｍ≦Ｍ｝の逆フィルタ｛Ｇ_ｍ（ｚ）；１≦ｍ≦Ｍ｝とＡＲ系のフィルタ｛１／（１−Ｂ_ｉ（ｚ））；１≦ｉ≦Ｆ｝の逆フィルタである予測誤差フィルタ｛１−Ａ_ｉ（ｚ）；１≦ｉ≦Ｆ｝とを並行して推定することを主な特徴とする。

この基本原理を説明するため、上述のモデル機構を組み込んだ系全体の構成線図を図１に示す。上述のモデル化に拠れば、原信号ｓ（ｔ）は、フレームごとのイノベーション系列ｅ_ｉ（１），…，ｅ_ｉ（Ｗ）に対してＡＲ系のフィルタ１／（１−Ｂ_ｉ（ｚ））を適用して得られる信号ｓ_ｉ（ｎ）のフレーム結合と看做すことができ、観測信号ｘ（ｔ）は、原信号ｓ（ｔ）に対して伝達関数Ｈ（ｚ）が作用したものと言える。そして、信号歪み除去処理は、観測信号ｘ（ｔ）に対して逆フィルタＧ（ｚ）を作用させて復元信号ｙ（ｔ）を得る処理となる。このとき、信号歪み除去処理で得られた復元信号ｙ（ｔ）をフレーム分割して、それぞれに対して、それぞれの信号に基づいて求めた予測誤差フィルタ１−Ａ_ｉ（ｚ）を適用して得られるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）はイノベーション系列ｅ_ｉ（１），…，ｅ_ｉ（Ｗ）に一致することが望ましい。もし、予測誤差フィルタ１−Ａ_ｉ（ｚ）の出力信号ｄ_ｉ（ｎ）がｄ_ｉ（ｎ）＝ｅ_ｉ（ｎ）〔１≦ｉ≦Ｆ，１≦ｎ≦Ｗ〕を満たすならば、［条件３］の条件下でΣ_ｍ＝１ ^ＭＨ_ｍ（ｚ）Ｇ_ｍ（ｚ）＝１となることが示せる（数学的証明については、参考文献Ａを参照されたい。）。つまり、ｓ（ｔ）＝ｙ（ｔ）が言える。このとき、１−Ａ_ｉ（ｚ）は１−Ｂ_ｉ（ｚ）に等しくなる。
（参考文献Ａ） Takuya Yoshioka, Takafumi Hikichi, Masato Miyoshi, Hiroshi G. Okuno: Robust Decomposition of Inverse Filter of Channel and Prediction Error Filter of Speech Signal for Dereverberation, Proceedings of the 14th European Signal Processing Conference (EUSIPCO 2006), CD-ROM Proceedings, Florence, 2006.

しかし、実際には、イノベーションｅ_ｉ（ｎ）〔１≦ｉ≦Ｆ，１≦ｎ≦Ｗ〕を信号歪み除去装置への入力信号として利用できない。図１に示す系において、各イノベーション系列ｅ_ｉ（ｎ）から観測信号ｘ（ｔ）を得る一連の過程は、モデル過程であって、実際には各イノベーション系列ｅ_ｉ（ｎ）、フィルタ１／（１−Ｂ_ｉ（ｚ））や伝達関数Ｈ（ｚ）を知ることはできないか知ることが困難であり、利用できる情報は観測信号ｘ（ｔ）のみである。そこで、上記［条件２］に基づいて、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を結合して得る全フレームでのイノベーション推定値系列が、その全サンプル間で独立になるように、つまりイノベーション推定値系列ｄ_１（１），…，ｄ_１（Ｗ），…，ｄ_ｉ（１），…，ｄ_ｉ（Ｗ），…，ｄ_Ｆ（１），…，ｄ_Ｆ（Ｗ）が独立となるように逆フィルタＧ_ｍ（ｚ）と予測誤差フィルタ１−Ａ_ｉ（ｚ）を推定する。

このことを、従来手法との対比で述べる。従来手法は、逆フィルタを、「観測信号に基づいて求めた予測誤差フィルタを観測信号に適用し、予測誤差フィルタ適用後の信号に対して逆フィルタを適用して得られる信号の正規化尖度が最大となる逆フィルタを求めよ」という問題の解として得ていた。これに対して、本発明は、逆フィルタを、「観測信号に逆フィルタを適用して得る信号に対して、当該信号に基づいて求めた予測誤差フィルタを適用して得る信号が、全サンプル間で独立となる逆フィルタを求めよ」という問題の解として得る。この問題で留意しなければならないことは、予測誤差フィルタが、観測信号に逆フィルタを適用して得る信号に基づいて求められるため、逆フィルタだけでなく予測誤差フィルタも一緒に求めることになるということである。
この問題は、ＩＣＡ（Independent Component Analysis）と同様の考え方によって定式化することができる。ここでは相互情報量を最小化する観点から説明を行うが、例えば対数尤度を最大化する方針で定式化することも可能である。いずれにしても問題に対するアプローチ方法の違いに過ぎない。

独立性の尺度として相互情報量（Kullback-Leibler情報量）を用いると、解くべき問題は式（７）のように定式化される。ただし、ｇ＝［ｇ_１ ^Ｔ，…，ｇ_Ｍ ^Ｔ］^Ｔ，ｇ_ｍ＝［ｇ_ｍ（０），…，ｇ_ｍ（Ｌ）］^Ｔ，ａ＝［ａ_１ ^Ｔ，…，ａ_Ｆ ^Ｔ］^Ｔ，ａ_ｉ＝［ａ_ｉ（１），…，ａ_ｉ（Ｐ）］^Ｔとし、ａ_ｉ（ｋ）は予測誤差フィルタ係数を表す。Ｉ（Ｕ_１，…，Ｕ_ｎ）は確率変数Ｕ_ｉ間の相互情報量を表す。またｇおよびａに記号＾を付したものは、得るべき最適解である。Ｔは転置を表す。

拘束条件
［１］ ‖ｇ‖＝１ （但し‖・‖はノルムを表す。）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

相互情報量Ｉは、イノベーション推定値系列ｄ_１（１），…，ｄ_１（Ｗ），…，ｄ_ｉ（１），…，ｄ_ｉ（Ｗ），…，ｄ_Ｆ（１），…，ｄ_Ｆ（Ｗ）の振幅が定数倍されても変化しない。式（７）の拘束条件[１]は、この振幅の不定性を排除するための条件である。式（７）の拘束条件[２]は、上記[条件１]に対応して、予測誤差フィルタを最小位相系に制限するための条件である。以下、Ｉを、イノベーション推定値系列を入力としそれらの間の相互情報量を出力する関数と看做して、損失関数と呼称することにする。

１．６ 損失関数の導出
式（７）の最適化を実行するためには、損失関数Ｉ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））を有限長の信号系列｛ｄ_ｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝から推定しなければならない。（多変量）確率変数Ｕの微分エントロピーをＤ（Ｕ）と表記すると、Ｉ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））は式（８）で定義される。ただし、ｄ＝［ｄ_Ｆ ^Ｔ，…，ｄ_１ ^Ｔ］^Ｔ、ｄ_ｉ＝［ｄ_ｉ（Ｗ），…，ｄ_ｉ（１）］^Ｔである。

ｙ＝［ｙ_Ｆ ^Ｔ，…，ｙ_１ ^Ｔ］^Ｔ，ｙ_ｉ＝［ｙ_ｉ（Ｗ），…，ｙ_ｉ（１）］^Ｔとおくと、ｄはｙを用いて、ｄ＝Ａｙと表される。ただし、行列Ａは、式（９）および式（１０）で表される。

よって、Ｄ（ｄ）は式（１１）のように表される。

多変量確率変数Ｕの共分散行列をΣ（Ｕ）と表記すると、式（１１）右辺第二項について、Σ（ｄ）＝Ｅ｛ｄｄ^Ｔ｝＝ＡＥ｛ｙｙ^Ｔ｝Ａ^Ｔ＝ＡΣ（ｙ）Ａ^Ｔが成立するから、式（１２）が成り立つ。

式（１１）、式（１２）を式（８）に代入すると、式（１３）を得る。ただし、σ（Ｕ）^２は確率変数Ｕの分散を表す。

式（１３）でＪ（Ｕ）は（多変量）確率変数Ｕのネゲントロピー（negentropy）である。ネゲントロピーはＵの非ガウス性の度合いを表す非負の値をとり、Ｕがガウス分布に従う場合に限り０をとる。Ｃ（Ｕ_１，…，Ｕ_ｎ）は式（１４）で定義される。Ｃ（Ｕ_１，…，Ｕ_ｎ）は確率変数Ｕ_ｉ間の相関の度合いを表す非負の値をとり、これらが無相関の場合に限り０をとる。

ところで、ｓ＝［ｓ_Ｆ ^Ｔ，…，ｓ_１ ^Ｔ］^Ｔ，ｓ_ｉ＝［ｓ_ｉ（Ｗ），…，ｓ_ｉ（１）］^ＴとおくとＪ（ｙ）＝Ｊ（ｓ）＝ｃｏｎｓｔａｎｔとなるため（証明略）、式（１３）は更に式（１５）のように簡単化できる。

以上から、式（１６）の最適化問題を解けばよいことになる。

拘束条件
［１］ ‖ｇ‖＝１ （但し‖・‖はノルムを表す。）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

１．７ 交代変数法による最適化
式（１６）について、交代変数の方法により、ｇとａを最適化する。すなわち、ｒ回目の繰り返しにおけるｇ及びａの推定値をそれぞれｇ＾^（ｒ）、ａ＾^（ｒ）と表せば、式（１７）および式（１８）の交互の最適化により更新された推定値ｇ＾^{（ｒ＋１）}，ａ＾^{（ｒ＋１）}を得る。なお、ｇ＾およびａ＾は、記号＾がｇ、ａのそれぞれの上に付されたものを表す。例えば繰り返し回数の上限をＲ_１とすれば、Ｒ_１回目で得られるｇ＾^{（Ｒ１＋１）}、ａ＾^{（Ｒ１＋１）}が式（１６）の最適解である。上付き文字のＲ１は、Ｒ_１である。

拘束条件
［１］ ｇ＝ｇ＾^（ｒ）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

拘束条件
［１］ ａ＝ａ＾^{（ｒ＋１）}
［２］ ‖ｇ‖＝１

式（１７）の意図するところは、伝達特性を打ち消すための逆フィルタの現在の推定値に基づいて原信号に固有の特性を打ち消すための予測誤差フィルタを推定することである。同様に、式（１８）の意図するところは、予測誤差フィルタの現在の推定値に基づいて逆フィルタを推定することである。イノベーション推定値系列ｄ_１（１），…，ｄ_１（Ｗ），…，ｄ_ｉ（１），…，ｄ_ｉ（Ｗ），…，ｄ_Ｆ（１），…，ｄ_Ｆ（Ｗ）が互いにより独立になるようにこれら２種類の最適化を繰り返すことで、逆フィルタと予測誤差フィルタを並行して推定することが可能になっている。したがって、ここでの繰り返しは逆フィルタの高精度な推定のために重要である。但し、図８から明らかなように処理する観測信号長が長くなる程、繰り返し回数は１回でも或る程度の効果が得られることが見て取れる。従って、この発明では、繰り返し回数は１回でもよい。

１．８ ａの最適化
本発明では、式（１７）の最適化を以下のように行う。
まず注意すべきことは、Ｃ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））はｄ_ｉ（ｎ）の二次の統計量に関連するのに対して、Ｊ（ｄ_ｉ（ｎ））はｄ_ｉ（ｎ）の高次の統計量に関連する値である。二次の統計量は信号の振幅情報のみ提供するが、高次の統計量は位相情報も提供する。したがって、一般に、高次統計量を含む最適化は、非最小位相系を導く可能性がある。そこで、１−Ａ_ｉ（ｚ）が最小位相系であるという拘束条件から、ａの最適化においては式（１９）の最適化問題を解く。

拘束条件
［１］ ｇ＝ｇ＾^（ｒ）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

Ｃ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））は式（２０）で与えられる。

ここで、行列Ａは式（９）および式（１０）に示すように上三角行列でその対角成分がすべて１であるから、ｌｏｇ ｄｅｔ Ａ＝０である。これを式（１２）に代入することで式（２１）の関係を得る。

よって、式（１９）は、式（２２）の最適化問題と等価である。なお、式（２２）は、上記[条件２]を反映した表現になっていることに留意しなければならない。式（２２）を説明すれば、式（２２）は、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるａを求めよ、と云っている。

拘束条件
［１］ ｇ＝ｇ＾^（ｒ）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

式（２２）で表される最適化問題を解くことは、観測信号にｇ＾^（ｒ）で与えられる逆フィルタを適用して得られるアドホック信号に対して、各フレームにおいて線形予測分析を行うことと等価であり、必ず最小位相予測誤差フィルタを得ることができる。線形予測分析に関しては、上記参考文献１を参照されたい。

なお、式（２２）では、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大きい。この場合、対数関数は単調増加関数であるから、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めることができる。

１．９ ｇの最適化
本発明では、式（１８）の最適化を以下のように行う。
前述したとおり、Ｃ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））は｛ｄ_ｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝の相関の度合いに関わる指標であるが、（ｒ＋１）回目のａの最適化において最小化されているため、Σ_ｉ＝１ ^ＦΣ_ｎ＝１ ^ＷＪ（ｄ_ｉ（ｎ））に比べて無視できる。そこでｇの最適化においては、式（２３）の最適化問題を解く。

拘束条件
［１］ ａ＝ａ＾^{（ｒ＋１）}
［２］ ‖ｇ‖＝１

Ｊ（ｄ_ｉ（ｎ））は、［条件２］に基づいて、式（２４）によって近似できる。この詳細は参考文献２を参照されたい。ただし、確率変数Ｕについて、κ_４（Ｕ）はＵの尖度（４次キュムラント）を表す。式（２４）の右辺を第ｉフレームにおける正規化尖度という。
（参考文献２） A.Hyvarinen, J.Karhunen, E.Oja, "INDEPENDENT COMPONENT ANALYSIS", John Wiley & Sons, Inc. 2001.

［条件２］から音声信号のイノベーションの尖度は正であるため、κ_４（ｄ_ｉ（ｎ））／σ（ｄ_ｉ（ｎ））^４は正である。従って、式（２３）の最適化問題は、式（２５）の最適化問題に帰着する。σ（ｄ_ｉ（ｎ）），κ_４（ｄ_ｉ（ｎ））は、［条件１］で述べた音声信号の局所的な定常性に基づいて、各フレーム内のサンプルから計算される。式（２６）では、１／Ｗを付加しているが、これは後の計算の便宜に過ぎず、式（２５）でｇの最適解を求めるにあたり影響を及ぼすものではない。式（２５）および式（２６）から、正規化尖度の全フレームでの総和が最大となるときのｇが、ｇ＾^{（ｒ＋１）}となる。なお、式（２５）および式（２６）は、上記[条件２]を反映した表現になっていることに留意しなければならない。式（２５）および式（２６）を説明すれば、これらは、第ｉフレームにおける正規化尖度を全フレームで加算した値が最大となるｇを求めよ、と云っている。

拘束条件
［１］ ａ＝ａ＾^{（ｒ＋１）}
［２］ ‖ｇ‖＝１

式（２５）に従ってｇの最適解を求めるには、Ｑをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（２７）で表される更新則に従って求められる。ｇ′をｇ′のノルムで除しているのは上記拘束条件[２]を課すためである。η（ｕ）は学習率を表す。ｕは、ｇの最適化における更新回数を表す。

式（２７）において、∇Ｑ_ｇは式（２８）および式（２９）で与えられる。

式（２９）において、ｄ_ｉ（ｎ）は式（３０）で、ｖ_ｍｉ（ｎ）は式（３１）および式（３２）で与えられる。ｘ_ｍｉ（ｎ）は、ｍ番目のマイクロホンで観測された観測信号のｉ番目のフレームの信号である。

§２ 二次統計量に基づく信号歪み除去
上述の従来的手法の信号歪み除去方法は、比較的長時間の観測信号（例えば２０秒程度である。）を要する。これは、一般に、正規化尖度のような高次統計量を計算するためには大量の観測信号のサンプルが必要となるからである。しかし、実際にはそうした長時間の観測信号を利用できない場合が多い。このため、従来的手法の信号歪み除去方法の適用分野は極めて限られていた。
また高次統計量の計算は比較的複雑であるため、従来的手法の信号歪み除去方法では装置の構成が複雑になりやすい。
そこで、観測信号がより短時間（例えば３秒から５秒程度である。）の場合にも有効であり、かつ計算が従来に比して容易な信号歪み除去の原理を説明する。この原理は、信号の二次統計量のみを用いるものであり、§１で説明した本発明の基本原理より派生する。

２．１ 二次統計量に基づく信号歪み除去の原理
二次統計量に基づく信号歪み除去では、上述の３つの条件に、次の２つの条件を設定する。
[条件４] Ｍ≧２である。すなわち、複数本のマイクロホンを用いる。
[条件５] Ｈ_ｍ＝｛ｈ_ｍ（ｋ）｝_ｋ＝０ ^Ｋは相異なるマイクロホンｍの間で共通の零点を持たない。

上記の式（１６）の最適化問題では、高次の統計量に関する値であるネゲントロピーＪおよび確率変数間の相関の度合いを示す指標Ｃを含む値を最小化するｇおよびａを求めた。
確率変数間の相関の度合いを示す指標Ｃは、二次の統計量で規定される。そこで、解くべき最適化問題を式（３３）で定式化する。

式（２１）を参酌すれば、式（３３）の最適化問題は、式（３４）の最適化問題に転化される。なお、式（３４）は、上記[条件２]を反映した表現になっていることに留意しなければならない。式（３４）を説明すれば、式（３４）は、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるｇおよびａを求めよ、と云っている。

ところで、上記の[条件４]および[条件５]が成立する場合、多チャンネルの観測信号は、音源からの原信号によって駆動されるＡＲ系として捉えることができる（参考文献３参照）。このことは、逆フィルタＧの先頭タップを式（３５）のように固定できることを意味する。但し、ｍ＝１に相当するマイクロホンは、最も音源に近いマイクロホンである。
（参考文献３） K. Aded-Meraim, E. Moulines, and P. Loubaton. Prediction error method for second-order blind identification. IEEE Trans. Signal Processing, Vol. 45, No.3, pp. 694-705, 1997.

式（３４）および式（３５）で規定されるｇを係数とする逆フィルタＧを、式（６）に従って観測信号ｘ（ｔ）に適用することで伝達特性が除去された復元信号ｙ（ｔ）を得る。

２．２ ａの最適化
式（３４）について、交代変数の方法により、ｇとａを最適化する。
逆フィルタの係数ｇ_ｍ（ｋ）を固定した状態で予測誤差フィルタの係数ａ_ｉ（ｋ）に関して式（３４）の損失関数を最小化する。
このとき、次の２点に注意する。１点目は、ｇ＝［ｇ_１ ^Ｔ，…，ｇ_Ｍ ^Ｔ］^Ｔは固定されているので、逆フィルタＧの出力である復元信号ｙ（ｔ）は予測誤差フィルタの最適化において不変であるということである。２点目は、ｉ番目のフレームにおける予測誤差フィルタの係数ａ_ｉ（１），…，ａ_ｉ（Ｐ）は、ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）にのみ寄与するということである。

従って、各々のフレームにおいて、Σ_ｎ＝１ ^Ｗlogσ（ｄ_ｉ（ｎ））^２を最小化するように、予測誤差フィルタの係数ａ_ｉ（１），…，ａ_ｉ（Ｐ）を推定すればよい。[条件２]からイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散はフレーム内で定常であるから、Σ_ｎ＝１ ^Ｗlogσ（ｄ_ｉ（ｎ））^２の最小化は、Ｗ＊σ（ｄ_ｉ（ｎ））^２の最小化と等価である。記号＊は乗算を表す。分散σ（ｄ_ｉ（ｎ））^２は、＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗとして計算できる。但し、＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗは、１フレーム分のイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を使って計算したｄ_ｉ（ｎ）の２乗平均を表す。結局、係数ａ_ｉ（ｋ）は、Ｗ＊＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗ、つまりｄ_ｉ（ｎ）の２乗総和が最小になるときのａ_ｉ（ｋ）として推定される。このような係数ａ_ｉ（ｋ）は、線形予測分析の手法によって計算される。

なお、ここでは、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めるとして説明したが、これに限定する趣旨ではない。既述のとおり、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めることができる。

２．３ ｇの最適化
予測誤差フィルタの係数ａ_ｉ（ｋ）を固定した状態で逆フィルタの係数ｇ_ｍ（ｋ）に関して式（３４）の損失関数を最小化する。
逆フィルタの係数ｇ_ｍ（ｋ）に関する損失関数の最小化には、勾配法を用いる。[条件２]を用いると、式（３４）の最適化問題は、式（３６）の最適化問題に転化される。

式（３６）に従ってｇの最適解を求めるには、Σ_ｉ＝１ ^Ｆlog＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（３７）で表される更新則に従って求められる。δは学習率を表す。１≦ｍ≦Ｍ、１≦ｋ≦Ｌである。なお、式（３７）では、式（３５）の条件から‖ｇ‖＝１なる拘束条件が課されないことに留意しなければならない。同様に式（３５）の条件からｋが取る値の範囲は１≦ｋ≦Ｌである。

上記式（２９）あるいは上記非特許文献１に記載された式（３）と比べて明らかなように、式（３７）の右辺第二項は、二次の統計量で表されており、この計算に高次の統計量を必要としない。このため、高次統計量を計算するには短い時間の観測信号の場合にも有効であり、計算自体も容易である。

なお、式（３６）では、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのｇをｇ＾として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大である。この場合、対数関数は単調増加関数であるから、第ｉフレームにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の全フレームでの総和が最小となるときのｇをｇ＾として求めることができる。なお、この場合には式（３７）で示した更新則は適用できなくなり、改めてΣ_ｉ＝１ ^Ｆ＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この結果得られる更新則はＩＣＡと同様の考え方で定式化できるから略する。

§３ プリ・ホワイトニング
本発明の信号歪み除去には、プリ・ホワイトニングを適用することができる。観測信号をプリ・ホワイトニングすることで、最適化計算の安定化、とくに更新則の高速な収束が可能となる。
各マイクロホンで得られた観測信号系列全体｛ｘ_ｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆ_ｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。
式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号ｘ_ｍ（ｔ）に適用する。ｗ_ｍ（ｔ）は、ｍ番目のマイクロホンで得た観測信号ｘ_ｍ（ｔ）を白色化した信号を表す。

このとき、式（３１）および式（３８）は式（４０）に、式（３２）は式（４１）に、変更すればよい。

§４ 実施形態
以下、本発明の実施形態を図面を参照して説明する。本発明の実施形態として、後述の各実施形態に限定するものではなく、各セクションで説明した原理を実現する実施形態であればよい。

＜第１実施形態＞
本発明の第１実施形態を実施する場合、以下の手順に従ってセンサで得た観測信号を処理する。ここでは、実施形態を具体的に説明する観点から信号として音声信号を例に挙げて説明する。
なお、第１実施形態の説明に先立ち、観測信号およびフレーム化処理について概説する。

((観測信号))
図示しないセンサ（例えばマイクロホン）によって得られたアナログ信号（このアナログ信号には伝達特性に由来する歪みが重畳されている。）は、例えば８,０００Ｈｚのサンプリングレートでサンプリングされ、適宜量子化された離散信号に変換される。以下、この離散信号を観測信号ということにする。アナログ信号から観測信号へのＡ／Ｄ変換などを実行するために必要となる構成要素（手段）は、いずれも公知技術の常套手段によって達成されるから、説明および図示を略する。

((フレーム化処理))
図示しない信号切出手段が、離散信号から、時間軸方向に一定時間幅で始点を移動させながら、所定時間長の離散信号をフレームとして切り出す。例えば２００サンプル点（８,０００Ｈｚ×２５ｍｓ）長の離散信号を、８０サンプル点（８,０００Ｈｚ×１０ｍｓ）ずつ始点を移動させながら切り出す。各フレームの切り出しは、離散信号に公知の窓関数（例えば、ハミング窓、ガウス窓、方形窓など）を適用すればよい。窓関数の適用によるフレームの切り出しは公知の常套手段によって達成される。

本発明の第１実施形態である信号歪み除去装置（１）をコンピュータ（汎用機）で実現する場合のハードウェア構成例を説明する。
図２に例示するように、信号歪み除去装置（１）は、キーボード、ポインティングデバイスなどが接続可能な入力部（１１）、液晶ディスプレイ、ＣＲＴ（Cathode Ray Tube）ディスプレイなどが接続可能な出力部（１２）、信号歪み除去装置（１）外部に通信可能な通信装置（例えば通信ケーブル、ＬＡＮカード、ルータ、モデムなど）が接続可能な通信部（１３）、ＤＳＰ（Digital Signal Processor）（１４）〔ＣＰＵ（Central Processing Unit）でも良い。またキャッシュメモリやレジスタ（１９）などを備えていてもよい。〕、メモリであるＲＡＭ（１５）、ＲＯＭ（１６）やハードディスク、光ディスク、半導体メモリなどである外部記憶装置（１７）並びにこれらの入力部（１１）、出力部（１２）、通信部（１３）、ＤＳＰ（１４）、ＲＡＭ（１５）、ＲＯＭ（１６）、外部記憶装置（１７）間のデータのやり取りが可能なように接続するバス（１８）を有している。また必要に応じて、信号歪み除去装置（１）に、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＤＶＤ（Digital Versatile Disc）などの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。

信号歪み除去装置（１）の外部記憶装置（１７）には、信号歪み除去のためのプログラムおよびこのプログラムの処理において必要となるデータ（観測信号）などが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶され、他のプログラムの処理に供されるときに、ＲＡＭや外部記憶装置などから読み込まれる。

より具体的には、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭなど〕には、観測信号に逆フィルタを適用する処理のためのプログラム、観測信号に逆フィルタを適用して得られた信号から予測誤差フィルタを求める処理のためのプログラム、予測誤差フィルタから逆フィルタを求める処理のためのプログラム、およびこれらのプログラムの処理において必要となるデータ（フレーム単位の観測信号など）が記憶されている。その他、これらのプログラムに基づく処理を制御するための制御プログラムも適宜に保存しておく。

第１実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、制御部）を実現することで、信号歪み除去が実現される。

そこで次に、図３〜図５を参照して、信号歪み除去装置（１）における信号歪み除去処理の流れを順次説明する。
大まかな処理の手順は、（ａ）観測信号ｘ（ｔ）に対して逆フィルタを適用した信号（以下、アドホック信号という。）を求め、（ｂ）アドホック信号から予測誤差フィルタを求め、（ｃ）この予測誤差フィルタから逆フィルタを求め、（ｄ）前記（ａ）、（ｂ）、（ｃ）の処理を繰り返して最適な逆フィルタを求め、（ｅ）最適化された逆フィルタを観測信号に対して適用した信号を復元信号ｙ（ｔ）として得る。
（ｂ）は上述のａの最適化に相当し、（ｃ）は上述のｇの最適化に相当し、（ｄ）は、式（１７）および式（１８）に相当する。（ｄ）の処理の繰り返し回数は予め定めた回数Ｒ_１とする。つまり、１≦ｒ≦Ｒ_１とする。また、（ｃ）の処理でｇを最適化する更新則の更新回数は予め定めた回数Ｒ_２とする。つまり、１≦ｕ≦Ｒ_２とする。（ｄ）の処理、つまり（ａ）、（ｂ）、（ｃ）の一連の処理を１回行うたびに、更新則によるＲ_２回の更新が行なわれる。実施形態では、回数Ｒ_１は、予め定めた回数とするが、これに限定されず、例えば、ｒ回目のｇを算出したときの式（２６）のＱの値とｒ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。同様に、回数Ｒ_２は、予め定めた回数とするが、これに限定されず、例えば、ｕ回目のｇを算出したときの式（２６）のＱの値とｕ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。

（ステップＳ１００）
逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタを観測信号ｘ（ｔ）＝［ｘ_１（ｔ），…，ｘ_ｍ（ｔ），…，ｘ_Ｍ（ｔ）］^Ｔに適用することで、アドホック信号ｙ（ｔ）を求める。アドホック信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、ここでは後述のＲ_１回の処理を経て求められた復元信号ではないことを明示するためアドホック信号と呼称する。ここでｔは、全てのサンプル番号を示し、１≦ｔ≦Ｎとする。Ｎは全サンプル数である。第１実施形態では、マイクロホンの個数Ｍは１以上であればよい。

逆フィルタの係数列｛ｇ_ｍ（ｋ）；０≦ｋ≦Ｌ｝として、繰り返し回数Ｒ_１の初回には予め定められた初期値を、２回目以降には後述する逆フィルタ計算部（１３）によって求められた逆フィルタｇ＾^{（ｒ＋１）}を用いる。

（ステップＳ１０１）
予測誤差フィルタ計算部（１５）は、フレーム化処理を行うフレーム化処理部（１５１）とフレーム予測誤差フィルタ計算部（１５２）によって構成される。そして、フレーム予測誤差フィルタ計算部（１５２）は、第ｉフレームのアドホック信号から予測誤差フィルタを求める第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）からなる。ただし、ｉは、１≦ｉ≦Ｆを満たす整数である。

フレーム化処理部（１５１）は、逆フィルタ適用部（１４）で求められたアドホック信号｛ｙ（ｔ）；１≦ｔ≦Ｎ｝をフレーム化処理する。フレーム化処理は、例えば式（４３）のように、Ｗ点分を切り出す窓関数をＷ点ずつシフトさせて適用することにより行う。｛ｙ_ｉ（ｎ）；１≦ｎ≦Ｗ｝はｉ番目のフレームに含まれるアドホック信号列を表す。

そして、第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）は、式（２２）に従って、第ｉフレームのアドホック信号列｛ｙ_ｉ（ｎ）；１≦ｎ≦Ｗ｝に対してＰ次の線形予測分析を行い、予測誤差フィルタの係数列｛ａ_ｉ（ｋ）；１≦ｋ≦Ｐ｝を計算する。この算出方法は、上記参考文献１を参照されたい。ここで得られたａ_１（１），…，ａ_１（Ｐ），…，ａ_ｉ（１），…，ａ_ｉ（Ｐ），…，ａ_Ｆ（１），…，ａ_Ｆ（Ｐ）は、式（２２）のａ＾^{（ｒ＋１）}を与える。

（ステップＳ１０２）
逆フィルタ計算部（１３）の機能構成例を、図４を参照して説明する。逆フィルタ計算部（１３）は勾配計算部（１３１）、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。更に、勾配計算部（１３１）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）とを備えて構成される。ここで更新用逆フィルタは、式（２７）のｇ^〈ｕ〉に相当する。

第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘ_ｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘ_ｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖ_ｍｉ（ｎ）を計算する（式（３１）を参照）。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。

第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｙ_ｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用して予測誤差フィルタ適用後のイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。爾後、第２の予測誤差フィルタ適用部（１３１２）は、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）を入力とする。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。

勾配ベクトル計算部（１３１３）は、信号ｖ_ｍｉ（ｎ）とイノベーション推定値ｄ_ｉ（ｎ）とを用いて現在の更新用逆フィルタｇ^〈ｕ〉の勾配ベクトル∇Ｑ_ｇを計算する（式（２８）および式（２９）を参照）。有限個のサンプルｖ_ｍｉ（ｎ）およびｄ_ｉ（ｎ）を用いて式（２９）を演算するときは、期待値Ｅをサンプルから求めればよい。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。

逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇ^〈ｕ〉、学習率η（ｕ）、勾配ベクトル∇Ｑ_ｇを入力として、式（２７）に従って、ｕ＋１回目の更新用逆フィルタｇ^{〈ｕ＋１〉}を求める。式（２７）は、求められたｇ^{〈ｕ＋１〉}を新たなｇ^〈ｕ〉と見立てて更新を行なうことを意味する。

更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇ^{〈ｕ＋１〉}、つまり新たなｇ^〈ｕ〉および観測信号ｘ（ｔ）を入力として、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。具体的には、式（４２）のｇ_ｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。なお、更新用逆フィルタ適用後信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、ここでは後述のＲ_１回の処理を経て求められた復元信号ではなく、更新則を行なうために算出される信号であることを明示するため更新用逆フィルタ適用後信号と呼称する。

制御部（６００）の制御によってＲ_２回の更新が行なわれた結果として得られたｇ^{〈Ｒ２＋１〉}は、式（２５）のｇ＾^{（ｒ＋１）}に相当する。上付き文字のＲ２は、Ｒ_２である。逆フィルタ計算部（１３）は、ｇ＾^{（ｒ＋１）}を出力する。

制御部（５００）の制御によって、上述の一連の処理を１回行うごとにｒに１を加算してｒがＲ_１に等しくなるまで、つまり上述の一連の処理をＲ_１回繰り返すことで（ステップＳ１０３）、ｇ＾^{（Ｒ１＋１）}を得る。上付き文字のＲ１は、Ｒ_１である。このｇ＾^{（Ｒ１＋１）}が、式（１６）の最適解とされる。そこで、ｇ＾^{（Ｒ１＋１）}を得た段階で、逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタｇ＾^{（Ｒ１＋１）}を観測信号ｘ（ｔ）＝［ｘ_１（ｔ），…，ｘ_Ｍ（ｔ）］^Ｔに適用することで、復元信号ｙ（ｔ）を得ることができる（ステップＳ１０４）。

＜第２実施形態＞
第２実施形態は、第１実施形態の変形例に相当する。具体的には、§３で述べたプリ・ホワイトニングを行なう形態である。そこで、第１実施形態と異なる部分について図６および図７を参照して説明を加える。なお、プリ・ホワイトニングは観測信号に対して行なうプリ・プロセスであるから、ここで説明するプリ・ホワイトニングを行なう形態は、後述の第３実施形態にも適用可能である。
第２実施形態では、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭなど〕に、白色化フィルタを求める処理のためのプログラム、白色化フィルタを観測信号に適用する処理のためのプログラムも記憶されている。
第２実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、白色化フィルタ計算部、白色化フィルタ適用部）を実現することで、信号歪み除去が実現される。

（ステップＳ１００ａ）
白色化フィルタ計算部（１１）は、各マイクロホンで得られた観測信号全体｛ｘ_ｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆ_ｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。この計算は線形予測分析と同じであり上記参考文献１を参照されたい。白色化フィルタの係数は、白色化フィルタ適用部（１２）の入力となる。

（ステップＳ１００ｂ）
白色化フィルタ適用部（１２）は、式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号に適用して、白色化信号ｗ_ｍ（ｔ）を得る。既述のとおり、式（３１）は式（４０）に変更すればよいので、第１実施形態において、逆フィルタ計算部（１３）、とくに第１の予測誤差フィルタ適用部（１３１１）による処理を式（３１）ではなく式（４０）による計算処理に改めればよい。また、第１実施形態において、逆フィルタ適用部（１４）による処理を、式（４２）ではなく式（４４）による計算処理に改めればよい。ステップＳ１００ｂの処理の後、第１実施形態のステップＳ１００〜Ｓ１０４の処理を行うが、これらの処理では第１実施形態の各処理における観測信号をステップＳ１００ｂの処理で得られた白色化信号に読み替えて第１実施形態と同様の処理を行う。この意味で、図７では、第１実施形態のステップＳ１００〜Ｓ１０４の各処理に相当する処理を示す符号に記号′を付している。

＜実施例１＞
発明者らは第２実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、マイクロホンの数Ｍ＝４、白色化フィルタの次数Ｘ＝５００、逆フィルタの次数Ｌ＝１０００、窓関数の切出しサンプル数（１フレームのサンプル数）Ｗ＝２００、予測誤差フィルタの次数Ｐ＝１６、繰り返し回数Ｒ_１＝１０、逆フィルタ計算部の更新回数Ｒ_２＝２０とした。学習率η（ｕ）は、初期値を０．０５に設定し、もし式（２７）によって式（２６）の値が減少するならば、η（ｕ）の値を順次半減させることで、式（２６）の値が必ず増大するようにした。図６に示した逆フィルタ適用部（１４）へ入力する初期逆フィルタは、式（４５）のように設定した。

本発明の第２実施形態の効果を、信号歪み除去の指標としてＤ_５０値（インパルス応答の全エネルギーに対する初期の５０ｍｓｅｃまでのエネルギーの比）を用いて評価した。連続発話データベースから男女各一名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
図８は、男声および女声について観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ_１（各信号長Ｎを持つ観測信号を図６に示す逆フィルタ適用部（１４）と、予測誤差フィルタ計算部（１５）と、逆フィルタ計算部（１３）を一巡する処理を実行して逆フィルタを求める回数）とＤ_５０値の関係を示している。いずれの場合においても、繰り返し回数を増加させるとＤ_５０値が向上しており、繰り返し処理の効果が顕著に見て取れる。特に観測信号長が５〜１０秒程度の比較的短い長さであっても、繰り返し処理によってＤ_５０値が大きく向上したことが分かる。

また、本発明の第２実施形態の効果を、音声スペクトグラムの比較から検証した。
図９Ａは観測信号長が１分の場合の残響を含まない音声（原音声）、図９Ｂは観測信号長が１分の場合の残響を含む音声（観測音声）、図９Ｃは観測信号長が１分の場合の残響除去後の音声（復元音声）のスペクトログラムの抜粋を示している。図９Ａと図９Ｃとの対比および図９Ｂと図９Ｃとの対比から、観測信号に含まれる残響が抑制され、原音声固有の特徴である調波構造やフォルマント構造が回復されたことが分かる。

また、本発明の第２実施形態の効果を、ＬＰＣスペクトル歪みを用いて検証した。
図１０Ｂは原音声の波形、図１０Ａは原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）および原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）を示している。図１０Ａおよび図１０Ｂの各横軸は秒単位の時間を表し、両図でタイムスケールを揃えている。図１０Ｂの縦軸は振幅値を表している。但し、原信号の相対的な大きさが分かればよいので、この趣旨から図１０Ｂの縦軸では単位を明示していない。図１０Ａの縦軸は、ＬＰＣスペクトル歪みＳＤ（ｄＢ）を表している。
図１０Ａから、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりも常に小さいことがわかる。なお、観測音声ではＬＰＣスペクトル歪みの平均が５．３９ｄＢ、分散が４．２０ｄＢであったのに対して、復元音声では平均が２．３８ｄＢ、分散が２．００ｄＢであった。
また、図１０Ａと図１０Ｂとの対比から、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）が大きな値を示す区間（例えば約１．０秒〜約１．２秒の区間を参照）は、原音声の波形の振幅値がほぼ０であることがわかる。実際、この区間では発声がなく無音区間である。このため、実際に知覚される歪みはかなり小さくなっていた。つまり、発声区間における原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりもかなり小さく、このため原音声のスペクトルを高い精度で復元できたことが結論付けられる。

＜第３実施形態＞
第３実施形態は、第１実施形態の変形例に相当する。具体的には、§２で述べた二次統計量に基づく信号歪み除去処理を行なう形態である。そこで、第１実施形態と異なる部分について図１１および図１２を参照して説明を加える。但し、第３実施形態では、マイクロホンの個数Ｍは２以上とする。

ステップＳ１００の処理およびステップＳ１０１の処理は、第１実施形態と同じである。

ステップＳ１０１の処理に続いて、ステップＳ１０２ａの処理を行う。
第３実施形態に係る逆フィルタ計算部（１３）の機能構成例を、図１１を参照して説明する。
逆フィルタ計算部（１３）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）と、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。ここで更新用逆フィルタは、式（３７）のｇ_ｍ（ｋ）に相当する。

第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘ_ｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘ_ｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖ_ｍｉ（ｎ）を計算する（式（３８）を参照）。具体的には、フレーム化処理部（４０２Ｂ）が、入力された観測信号ｘ_ｍ（ｔ）に対してフレーム化処理を行い、観測信号ｘ_ｍ（ｔ）のｉ番目のフレームの信号ｘ_ｍｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０４ｉ）が信号ｘ_ｍｉ（ｎ）を入力として、式（３８）に従って信号ｖ_ｍｉ（ｎ）を出力する。但し、１≦ｉ≦Ｆである。

第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、各フレームに対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用して予測誤差フィルタ適用後のイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。具体的には、フレーム化処理部（４０２Ａ）が、初期値の場合を除き、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）に対してフレーム化処理を行い、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、ｉ番目のフレームの信号ｙ_ｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０３ｉ）が信号ｙ_ｉ（ｎ）を入力として、式（３０）に従ってイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を出力する。但し、１≦ｉ≦Ｆである。

勾配ベクトル計算部（１３１３）は、信号ｖ_ｍｉ（ｎ）とイノベーション推定値ｄ_ｉ（ｎ）とを用いて現在の更新用逆フィルタｇ_ｍ（ｋ）の勾配ベクトルを計算する（式（３７）の右辺第二項を参照）。具体的には、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、相互相関計算部（４０５ｉ）は信号ｖ_ｍｉ（ｎ）とイノベーション推定値ｄ（ｎ）との相互相関〈ｄ_ｉ（ｎ）ｖ_ｍｉ（ｎ−ｋ）〉_ｎ＝１ ^Ｗを計算する。また、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、分散計算部（４０６ｉ）は、イノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散〈ｄ_ｉ（ｎ）^２〉_ｎ＝１ ^Ｗを求める。各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、除算部（４０７ｉ）は、〈ｄ_ｉ（ｎ）ｖ_ｍｉ（ｎ−ｋ）〉_ｎ＝１ ^Ｗ／〈ｄ_ｉ（ｎ）^２〉_ｎ＝１ ^Ｗを求める。加算部（４０７）は、除算部（４０７１）〜（４０７Ｆ）の全フレームに亘る総和、つまり式（３７）の右辺第二項を求める。

逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇ_ｍ（ｋ）、学習率δ、勾配ベクトルを入力として、式（３７）に従って、ｕ＋１回目の更新用逆フィルタｇ_ｍ（ｋ）′を求める。式（３７）は、求められたｇ_ｍ（ｋ）′を新たなｇ_ｍ（ｋ）と見立てて更新を行なうことを意味する。

更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇ_ｍ（ｋ）′、つまり新たなｇ_ｍ（ｋ）および観測信号ｘ（ｔ）を入力として、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。具体的には、式（４２）のｇ_ｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。

ステップＳ１０２ａの処理に続いて、ステップＳ１０３およびステップＳ１０４の処理を行うが、第１実施形態と同じであるから説明を略する。

＜実施例２＞
発明者らは第３実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、Ｍ＝４，Ｌ＝１０００，Ｗ＝２００，Ｐ＝１６，Ｒ_１＝６，Ｒ_２＝５０とした。学習率δは、初期値を０．０５に設定し、Σ_ｉ＝１ ^Ｆlog＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗの値が増加するならば、学習率δの値を順次半減させることでΣ_ｉ＝１ ^Ｆlog＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗの値が必ず減少するように設定した。逆フィルタの初期推定値は、ｇ_ｍ（ｋ）＝０，１≦ｍ≦Ｍ，１≦ｋ≦Ｌとして設定した。

本発明の第３実施形態の効果を、音声明瞭度を表すＲＡＳＴＩ（参考文献５を参照）を残響除去の指標として評価した。連続発話データベースから男女各五名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
（参考文献５） H. kuttruff. Room acoustics. Elsevier Applied Science, third edition, P.237 1991.

図１３は、観測信号長Ｎを３秒、４秒、５秒、１０秒としたときの、ＲＡＳＴＩの値を表示したものである。図１３に示すように、観測信号が３〜５秒のように短時間の場合でも、高い残響除去性能を示していることが分かる。
図１４は、残響除去前後におけるエネルギー減衰曲線の例である。直接音が到達してから５０ミリ秒後の反射音のエネルギーが１５ｄＢ低減されていることが分かる。

本発明は、様々な信号処理システムの性能向上に寄与する要素技術であるところ、例えば音声認識システム、テレビ会議システム、補聴器、音楽情報処理システム等に利用することができる。

この発明は、信号の歪み（ひずみ：distortion）を除去する技術に関する。

信号は反射や残響などが存在する環境で観測されると、本来の信号に反射や残響などが重畳された信号として観測される。以下、本来の信号を「原信号」と云い、観測された信号を「観測信号」と云うことにする。また、反射や残響などに代表される、原信号に重畳された歪みを「伝達特性」と云うことにする。このため、観測信号から原信号固有の特徴を抽出することが困難になる。この不都合を解消すべく、かねてより信号歪み除去処理技術が各種試みられてきた。信号歪み除去処理は、原信号に重畳した伝達特性を観測信号から取り除く処理である。

従来の信号歪み除去方法の一例として非特許文献１に開示されている信号歪み除去処理を図１５を用いて説明する。予測誤差フィルタ計算部（９０１）は、観測信号をフレーム化処理して、各フレームに含まれる観測信号に対して線形予測分析を行い、予測誤差フィルタを計算する。この明細書では、フィルタはディジタルフィルタであり、信号のサンプル値に対して作用するいわゆるフィルタ係数を求める意味で単にフィルタを計算するなどということがある。予測誤差フィルタ適用部（９０２）は、各フレームごとに上記計算された予測誤差フィルタを当該フレームの観測信号に適用する。逆フィルタ計算部（９０３）は、予測誤差フィルタ適用後の信号に対して逆フィルタを適用して得られる信号の正規化尖度が最大となるような逆フィルタを計算する。逆フィルタ適用部（９０４）は、上記計算された逆フィルタを観測信号に適用することで信号歪み除去後の信号（復元信号）を得る。
B.W.Gillespie, , H.S.Malvar, ,and D.A.F.Florencio, ,"Speech dereverberation via maximum-kurtosis subband adaptive filtering," IEEE International Conference on Acoustics, Speech, and Signal Processing, pp.3701-3704, 2001.

上記の従来的な信号歪み除去方法は、観測信号の各フレーム内ではショートラグ(short−lag)の自己相関は原信号固有の特性の寄与が大きく、フレームを跨ぐロングラグ（long-lag）の自己相関は伝達特性による寄与が大きいことを仮定している。上記従来方法は、この仮定に基づいて、フレーム化処理されたフレーム単位の観測信号に予測誤差フィルタを適用して観測信号中の原信号固有の特性の寄与を低減している。

しかし、この仮定は粗い近似であるため、推定された逆フィルタの精度は不十分である。つまり、観測信号から求まる予測誤差フィルタは伝達特性の影響を受けているので、原信号固有の特性のみを正しく取り除くことができない。このため、予測誤差フィルタ適用後の信号から求める逆フィルタの精度は劣化する。結果として、観測信号に逆フィルタを適用して得る信号は、本来の原信号の正確な推定値ではない。
そこで本発明は、伝達特性に由来する歪みを観測信号から除去して精度の良い復元信号を得ることを目的とする。

上記課題を解決するため、本発明の信号歪み除去装置は、所定の繰り返し終了条件を満たした場合には、観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、観測信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、観測信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段と、を備える。
この発明では、観測信号に伝達特性を除去するための逆フィルタを適用して得るアドホック信号に対して、このアドホック信号に基づいて求めた予測誤差フィルタを適用して得る信号（イノベーション推定値系列）が、その全サンプル間で独立となるような逆フィルタを求める。そして、所定の繰り返し終了条件を満たしたときの逆フィルタを観測信号に
適用することで復元信号を得る。

上記の信号歪み除去装置では、予測誤差フィルタ計算手段は、各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの予測誤差フィルタを出力するものであり、逆フィルタ計算手段は、上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタとして、各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタを求め、この逆フィルタを出力するものであるとしてもよい。
この構成は、イノベーション系列のサンプル間の独立性の尺度として相互情報量を規定し、これを最小化する予測誤差フィルタと逆フィルタを交代変数法で求めるものである。この詳細は後述する。

あるいは、上記の信号歪み除去装置では、予測誤差フィルタ計算手段は、各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの予測誤差フィルタを出力するものであり、逆フィルタ計算手段は、上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタとして、各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタを求め、この逆フィルタを出力するものであるとしてもよい。
この構成は、イノベーション系列のサンプル間の独立性の尺度として相互情報量を規定し、これを最小化する予測誤差フィルタと逆フィルタを交代変数法で求めるものであるが、信号の高次統計量を用いることなく予測誤差フィルタと逆フィルタを交代変数法で求めることができる。

上記の信号歪み除去装置では、プリ・ホワイトニング処理を前置させ、プリ・ホワイトニング処理で得られた白色化信号に対して、上記同様の処理を行う装置構成とすることができる。具体的には、観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算手段と、白色化フィルタを観測信号に適用して白色化信号を出力する白色化フィルタ適用手段と、所定の繰り返し終了条件を満たした場合には、白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、白色化信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、白色化信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段と、を備えた信号歪み除去装置とすることができる。

上記課題を解決するため、本発明の方法は、逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、観測信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、観測信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、予測誤差フィルタ計算手段が、アドホック信号を
フレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、逆フィルタ計算手段が、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、制御手段が、繰り返し終了条件を満たすまで逆フィルタ適用ステップ、予測誤差フィルタ計算ステップ、逆フィルタ計算ステップを繰り返し実行させる制御ステップと、を有する信号歪み除去方法とする。

また、上記の信号歪み除去方法では、プリ・ホワイトニング処理を前置させ、プリ・ホワイトニング処理で得られた白色化信号に対して、上記同様の処理を行う方法とすることができる。具体的には、白色化フィルタ計算手段が、観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算ステップと、白色化フィルタ適用手段が、白色化フィルタを観測信号に適用して白色化信号を出力する白色化フィルタ適用ステップと、逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、白色化信号に適用して、この結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、白色化信号に逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、予測誤差フィルタ計算手段が、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、逆フィルタ計算手段が、各フレームのアドホック信号に対して当該フレームに対応する予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値系列という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、制御手段が、繰り返し終了条件を満たすまで逆フィルタ適用ステップ、予測誤差フィルタ計算ステップ、逆フィルタ計算ステップを繰り返し実行させる制御ステップと、を有する信号歪み除去方法とする。

本発明の信号歪み除去装置としてコンピュータを機能させる信号歪み除去プログラムによって、コンピュータを信号歪み除去装置として作動処理させることができる。そして、この信号歪み除去プログラムを記録した、コンピュータに読み取り可能なプログラム記録媒体によって、他のコンピュータを信号歪み除去装置として機能させることや、信号歪み除去プログラムを流通させることなどが可能になる。

本発明では、観測信号から求まる予測誤差フィルタを用いて観測信号中の原信号固有の特性の寄与を低減するのではなく、観測信号に（仮の）逆フィルタを適用して得られるアドホック信号（仮の復元信号）から求まる予測誤差フィルタを用いて観測信号中の原信号固有の特性を低減する。アドホック信号から求まる予測誤差フィルタは、伝達特性の影響を受けにくいので、原信号固有の特性をより正確に取り除くことが可能である。このような予測誤差フィルタをアドホック信号に適用して得る信号（イノベーション推定値系列）が全サンプル間で独立となるように求められた逆フィルタは精度良く伝達特性を除去可能なものであるから、このような逆フィルタを観測信号に適用することで、伝達特性に由来する歪みが除去された精度の良い復元信号を得ることができる。

§１ 本発明の理論
以下、実施形態の説明に先立ち、本発明の理論を説明する。
以下の説明では、特に断りのない限り、信号源は１つとする。

１．１ 信号
本発明の対象となる信号は、人の音声、音楽、生体信号、測定対象物の物理量をセンサで観測した電気信号などの信号を広く包含する。より好ましくは、自己回帰（Autoregressive：ＡＲ）過程として表現することができる、あるいは表現することが好ましい信号であればよい。例えば音声信号は、通常、区分定常な自己回帰過程として表現される信号、すなわち独立同一分布（i.i.d. : Independent and Identically Distributed）信号に音韻性を表すＡＲ系を作用させた信号として看做される（参考文献１参照）。
以下、信号の代表例として音声信号を挙げて本発明の理論を説明する。
（参考文献１） L.R.Rabiner, R.W.Schafer, "Digital Processing of Speech Signals", Bell Laboratories, Incorporated, 1978.

１．２ 音声信号のモデル化
まず、原信号である音声信号ｓ（ｔ）を、以下の３つの条件を満足する信号としてモデル化する。

［条件１］音声信号ｓ（ｔ）は、区分定常なＡＲ過程で生成される。
この[条件１]から、ＡＲ過程の次数をＰ、定常とみなせる区間長をＷサンプルとして音声信号ｓ（ｔ）をフレーム化すると、第ｉフレームの音声信号ｓ_ｉ（ｎ）は、式（１）のように表される。式（２）は、第ｉフレームの音声信号ｓ_ｉ（ｎ）のサンプルと、フレーム化前の音声信号ｓ（ｔ）のサンプルとの対応を示している。つまり、第ｉフレームのｎ番目のサンプルは、フレーム化前の音声信号ｓ（ｔ）において、（ｉ−１）Ｗ＋ｎ番目のサンプルに相当する。式（１）および式（２）において、ｂ_ｉ（ｋ）は線形予測係数、ｅ_ｉ（ｎ）はイノベーションを表す。但し、１≦ｎ≦Ｗ、１≦ｔ≦Ｎ、Ｎは全サンプル数である。以下、特に断りの無い限り、パラメータｎは１フレームのサンプル番号を表し、パラメータｔは全てのサンプル番号を表す。また、全フレーム数はＦとする。

なお、第ｉフレームにおけるｎ番目のイノベーションｅ_ｉ（ｎ）についても、第ｉフレームのｎ番目のイノベーションｅ_ｉ（ｎ）と、フレーム化前の音声信号ｓ（ｔ）に対するイノベーションｅ（ｔ）との対応を示すことができる。この場合、第ｉフレームのｎ番目のイノベーションｅ_ｉ（ｎ）は、フレーム化前のイノベーションｅ（ｔ）において（ｉ−１）Ｗ＋ｎ番目のイノベーションに相当し、ｅ_ｉ（ｎ）＝ｅ（（ｉ−１）Ｗ＋ｎ）が成り立つ。

式（１）をｚ変換する。左辺のｚ変換をＳ_ｉ（Ｚ）とし、右辺第二項のｚ変換をＥ_ｉ（Ｚ）とし、Ｂ_ｉ（ｚ）＝Σ_ｋ＝１ ^Ｐｂ_ｉ（ｋ）ｚ^−ｋとすれば、右辺第一項は、Ｂ_ｉ（ｚ）Ｓ_ｉ（Ｚ）となる。従って、式（１）のｚ変換は、（１−Ｂ_ｉ（ｚ））Ｓ_ｉ（Ｚ）＝Ｅ_ｉ（Ｚ）である。なお、ｚ^−１は時間領域では１タップ遅延素子に相当する。以降、時間領域信号（タップ重み係数）を小文字で、ｚ領域信号（伝達関数）を大文字でそれぞれ表す。１−Ｂ_ｉ（ｚ）は最小位相性を満足しなければならず、『１−Ｂ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ』ことが要求される。

［条件２］第ｉフレームに属するイノベーションｅ _ｉ （１），…，ｅ_ｉ（Ｗ）は独立且つ
同一分布に属する。イノベーションｅ _ｉ （１），…，ｅ_ｉ（Ｗ）の確率分布の平均及び歪度（３次キュムラント）は０、尖度（４次キュムラント）は正である。さらに、異なるフレームｉ、ｊ〔ｉ≠ｊ〕に属するイノベーションｅ_ｉ（ｎ）とｅ_ｊ（ｎ′）同士も独立である。ただし、これらは必ずしも同一分布に属するとは限らない。
［条件３］予測誤差フィルタ１−Ｂ_ｉ（ｚ）は、相異なるフレーム間で共通する零点をもたない。

式（１）および式（２）から、音声信号ｓ（ｔ）は、式（３）のように表される。[・]は床関数を表す。

このとき、[条件２]は、『イノベーション過程ｅ（ｔ）は時間的に独立な信号である。また、その統計的性質（あるいは統計量）はフレーム内では定常である。』と表現できる。また、[条件３]は、『線形予測係数｛ｂ_ｉ（ｋ）｝_ｋ＝１ ^Ｐ で表されるＡＲ系は、時不変な極を持たない』と表現できる。

１．３ 観測信号のモデル化
次に、Ｍ個のマイクロホンで音声信号を観測して観測信号を得たときの観測信号をモデル化する。但しＭは、Ｍ≧１の整数である。
ｍ番目（１≦ｍ≦Ｍ）のマイクロホンで観測される残響信号ｘ_ｍ（ｔ）を、音源からｍ番目のマイクロホンに至る経路の伝達関数Ｈ_ｍ（ｚ）のタップ重み係数｛ｈ_ｍ（ｋ）；０≦ｋ≦Ｋ；Ｋはインパルス応答の持続時間とする。｝を用いて式（４）のようにモデル化する。ここでは、音声信号の場合の伝達特性の代表例として残響を挙げて、伝達特性を残響に言い換えて説明する。但し、伝達特性を残響に限定する趣旨ではない。

Ｍ個の観測信号についてまとめて表現すれば、式（５）のように表すことができる。但し、式（５）において、ｘ（ｔ）＝［ｘ_１（ｔ），…，ｘ_Ｍ（ｔ）］^Ｔ、ｈ（ｋ）＝［ｈ_１（ｋ），…，ｈ_Ｍ（ｋ）］^Ｔである。

１．４ 信号歪み除去の原理
信号歪み除去後の復元信号ｙ（ｔ）は、多チャネル逆フィルタ｛Ｇ_ｍ（ｚ）；１≦ｍ≦Ｍ｝のタップ重み係数｛ｇ_ｍ（ｋ）；１≦ｍ≦Ｍ，０≦ｋ≦Ｌ；Ｌは逆フィルタの次数｝を用いて式（６）により計算される。本発明においては、逆フィルタ係数であるｇ_ｍ（ｋ）を観測信号ｘ_１（ｔ），…，ｘ_Ｍ（ｔ）のみから推定する。

１．５ 本発明の基本原理
本発明の基本原理は、伝達関数｛Ｈ_ｍ（ｚ）；１≦ｍ≦Ｍ｝の逆フィルタ｛Ｇ_ｍ（ｚ）；１≦ｍ≦Ｍ｝とＡＲフィルタ｛１／（１−Ｂ_ｉ（ｚ））；１≦ｉ≦Ｆ｝の逆フィルタである予測誤差フィルタ｛１−Ａ_ｉ（ｚ）；１≦ｉ≦Ｆ｝とを並行して推定することを主な特徴とする。

この基本原理を説明するため、上述のモデル機構を組み込んだ系全体の構成線図を図１に示す。上述のモデル化に拠れば、原信号ｓ（ｔ）は、フレームごとのイノベーション系列ｅ_ｉ（１），…，ｅ_ｉ（Ｗ）に対してＡＲフィルタ１／（１−Ｂ_ｉ（ｚ））を適用して得られる信号ｓ _１ （ｎ）,…，ｓ _Ｆ （ｎ）のフレーム結合と看做すことができ、観測信号ｘ（ｔ）は、原信号ｓ（ｔ）に対して伝達関数Ｈ（ｚ）が作用したものと言える。そして、信号歪み除去処理は、観測信号ｘ（ｔ）に対して逆フィルタＧ（ｚ）を作用させて復元信号ｙ（ｔ）を得る処理となる。このとき、信号歪み除去処理で得られた復元信号ｙ（ｔ）をフレーム分割して、それぞれに対して、それぞれの信号に基づいて求めた予測誤差フィルタ１−Ａ_ｉ（ｚ）を適用して得られるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）はイノベーション系列ｅ_ｉ（１），…，ｅ_ｉ（Ｗ）に一致することが望ましい。もし、予測誤差フィルタ１−Ａ_ｉ（ｚ）の出力信号ｄ_ｉ（ｎ）がｄ_ｉ（ｎ）＝ｅ_ｉ（ｎ）〔１≦ｉ≦Ｆ，１≦ｎ≦Ｗ〕を満たすならば、［条件３］の条件下でΣ_ｍ＝１ ^ＭＨ_ｍ（ｚ）Ｇ_ｍ（ｚ）＝１となることが示せる（数学的証明については、参考文献Ａを参照されたい。）。つまり、ｓ（ｔ）＝ｙ（ｔ）が言える。このとき、１−Ａ_ｉ（ｚ）は１−Ｂ_ｉ（ｚ）に等しくなる。
（参考文献Ａ） Takuya Yoshioka, Takafumi Hikichi, Masato Miyoshi, Hiroshi G. Okuno: Robust Decomposition of Inverse Filter of Channel and Prediction Error Filter of Speech Signal for Dereverberation, Proceedings of the 14th European Signal
Processing Conference (EUSIPCO 2006), CD-ROM Proceedings, Florence, 2006.

しかし、実際には、イノベーションｅ_ｉ（ｎ）〔１≦ｉ≦Ｆ，１≦ｎ≦Ｗ〕を信号歪み除去装置への入力信号として利用できない。図１に示す系において、各イノベーション系列ｅ_ｉ（ｎ）から観測信号ｘ（ｔ）を得る一連の過程は、モデル過程であって、実際には各イノベーション系列ｅ_ｉ（ｎ）、ＡＲフィルタ１／（１−Ｂ_ｉ（ｚ））や伝達関数Ｈ _ｍ （ｚ）を知ることはできないか知ることが困難であり、利用できる情報は観測信号ｘ（ｔ）のみである。そこで、上記［条件２］に基づいて、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を結合して得る全フレームでのイノベーション推定値系列が、その全サンプル間で独立になるように、つまりイノベーション推定値系列ｄ_１（１），…，ｄ_１（Ｗ），…，ｄ_ｉ（１），…，ｄ_ｉ（Ｗ），…，ｄ_Ｆ（１），…，ｄ_Ｆ（Ｗ）が独立となるように逆フィルタＧ_ｍ（ｚ）と予測誤差フィルタ１−Ａ_ｉ（ｚ）を推定する。

ここに述べた本願の着想は、下記の点で従来手法と区別される。従来手法は、逆フィルタを、「観測信号に基づいて求めた予測誤差フィルタを観測信号に適用し、予測誤差フィルタ適用後の信号に対して逆フィルタを適用して得られる信号の正規化尖度が最大となる逆フィルタを求めよ」という問題の解として得ていた。これに対して、本発明は、逆フィルタを、「観測信号に逆フィルタを適用して得る信号に対して、当該信号に基づいて求めた予測誤差フィルタを適用して得る信号が、全サンプル間で独立となる逆フィルタを求め
よ」という問題の解として得る。この問題で留意しなければならないことは、予測誤差フィルタが、観測信号に逆フィルタを適用して得る信号に基づいて求められるため、逆フィルタだけでなく予測誤差フィルタも一緒に求めることになるということである。
この問題は、ＩＣＡ（Independent Component Analysis）と同様の考え方によって定式化することができる。ここでは相互情報量を最小化する観点から説明を行うが、例えば最尤推定法に基づいて定式化することも可能である。いずれにしても問題の定式化の違いに過ぎない。

独立性の尺度として相互情報量（Kullback-Leibler情報量）を用いると、解くべき問題は式（７）のように定式化される。ただし、ｇ＝［ｇ_１ ^Ｔ，…，ｇ_Ｍ ^Ｔ］^Ｔ，ｇ_ｍ＝［ｇ_ｍ（０），…，ｇ_ｍ（Ｌ）］^Ｔ，ａ＝［ａ_１ ^Ｔ，…，ａ_Ｆ ^Ｔ］^Ｔ，ａ_ｉ＝［ａ_ｉ（１），…，ａ_ｉ（Ｐ）］^Ｔとし、ａ_ｉ（ｋ）は予測誤差フィルタ係数を表す。Ｉ（Ｕ_１，…，Ｕ_ｎ）は確率変数Ｕ_ｉ間の相互情報量を表す。またｇおよびａに記号＾を付したものは、得るべき最適解である。Ｔは転置を表す。

拘束条件
［１］ ‖ｇ‖＝１ （但し‖・‖はノルムを表す。）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

相互情報量Ｉは、イノベーション推定値系列ｄ_１（１），…，ｄ_１（Ｗ），…，ｄ_ｉ（１），…，ｄ_ｉ（Ｗ），…，ｄ_Ｆ（１），…，ｄ_Ｆ（Ｗ）の振幅が定数倍されても変化しない。式（７）の拘束条件[１]は、この振幅の不定性を排除するための条件である。式（７）の拘束条件[２]は、上記[条件１]に対応して、予測誤差フィルタを最小位相系に制限するための条件である。以下、Ｉを、イノベーション推定値系列を入力としそれらの間の相互情報量を出力する関数と看做して、損失関数と呼称することにする。

１．６ 損失関数の導出
式（７）の最適化を実行するためには、損失関数Ｉ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））を有限長の信号系列｛ｄ_ｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝から推定しなければならない。（多変量）確率変数Ｕの微分エントロピーをＤ（Ｕ）と表記すると、Ｉ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））は式（８）で定義される。ただし、ｄ＝［ｄ_Ｆ ^Ｔ，…，ｄ_１ ^Ｔ］^Ｔ、ｄ_ｉ＝［ｄ_ｉ（Ｗ），…，ｄ_ｉ（１）］^Ｔである。

ｙ＝［ｙ_Ｆ ^Ｔ，…，ｙ_１ ^Ｔ］^Ｔ，ｙ_ｉ＝［ｙ_ｉ（Ｗ），…，ｙ_ｉ（１）］^Ｔとおくと、ｄはｙを用いて、ｄ＝Ａｙと表される。ただし、行列Ａは、式（９）および式（１０）で表される。

よって、Ｄ（ｄ）は式（１１）のように表される。

多変量確率変数Ｕの共分散行列をΣ（Ｕ）と表記すると、式（１１）右辺第二項について、Σ（ｄ）＝Ｅ｛ｄｄ^Ｔ｝＝ＡＥ｛ｙｙ^Ｔ｝Ａ^Ｔ＝ＡΣ（ｙ）Ａ^Ｔが成立するから、式（１２）が成り立つ。

式（１１）、式（１２）を式（８）に代入すると、式（１３）を得る。ただし、σ（Ｕ）^２は確率変数Ｕの分散を表す。

式（１３）でＪ（Ｕ）は（多変量）確率変数Ｕのネゲントロピー（negentropy）である。ネゲントロピーはＵの非ガウス性の度合いを表す非負の値をとり、Ｕがガウス分布に従う場合に限り０をとる。Ｃ（Ｕ_１，…，Ｕ_ｎ）は式（１４）で定義される。Ｃ（Ｕ_１，…，Ｕ_ｎ）は確率変数Ｕ_ｉ間の相関の度合いを表す非負の値をとり、これらが無相関の場合に限り０をとる。

ところで、ｓ＝［ｓ_Ｆ ^Ｔ，…，ｓ_１ ^Ｔ］^Ｔ，ｓ_ｉ＝［ｓ_ｉ（Ｗ），…，ｓ_ｉ（１）］^ＴとおくとＪ（ｙ）＝Ｊ（ｓ）＝ｃｏｎｓｔａｎｔとなるため（証明略）、式（１３）は更に式（１５）のように簡単化できる。

以上から、式（７）の最適化問題は式（１６）の最適化問題に等価である。

拘束条件
［１］ ‖ｇ‖＝１ （但し‖・‖はノルムを表す。）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

１．７ 交代変数法による最適化
式（１６）について、交代変数の方法により、ｇとａを最適化する。すなわち、ｒ回目の繰り返しにおけるｇ及びａの推定値をそれぞれｇ＾^（ｒ）、ａ＾^（ｒ）と表せば、式（１７）および式（１８）の交互の最適化により更新された推定値ｇ＾^{（ｒ＋１）}，ａ＾^{（ｒ＋１）}を得る。なお、ｇ＾およびａ＾は、記号＾がｇ、ａのそれぞれの上に付されたものを表す。例えば繰り返し回数の上限をＲ_１とすれば、Ｒ_１回目で得られるｇ＾^{（Ｒ１＋１）}、ａ＾^{（Ｒ１＋１）}が式（１６）の最適解である。上付き文字のＲ１は、Ｒ_１である。

拘束条件
［１］ ｇ＝ｇ＾^（ｒ）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

拘束条件
［１］ ａ＝ａ＾^{（ｒ＋１）}
［２］ ‖ｇ‖＝１

式（１７）の意図するところは、伝達特性を打ち消すための逆フィルタの現在の推定値に基づいて原信号に固有の特性を打ち消すための予測誤差フィルタを推定することである。他方、式（１８）の意図するところは、予測誤差フィルタの現在の推定値に基づいて逆フィルタを推定することである。イノベーション推定値系列ｄ_１（１），…，ｄ_１（Ｗ），…，ｄ_ｉ（１），…，ｄ_ｉ（Ｗ），…，ｄ_Ｆ（１），…，ｄ_Ｆ（Ｗ）が互いにより独立になるようにこれら２種類の最適化を繰り返すことで、逆フィルタと予測誤差フィルタを並行して推定することが可能になっている。したがって、ここでの繰り返しは逆フィルタの高精度な推定のために重要である。但し、図８から明らかなように処理する観測信号長が長くなる程、繰り返し回数は１回でも或る程度の信号歪み除去が達成されることが見て取れる。従って、この発明では、繰り返し回数は１回でもよい。

１．８ ａの最適化
本発明では、式（１７）の最適化を以下のように行う。
まず注意すべきことは、Ｃ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））はｄ_ｉ（ｎ）の２次の統計量に関連するのに対して、Ｊ（ｄ_ｉ（ｎ））はｄ_ｉ（ｎ）の高次の統計量に関連する値である。２次の統計量は信号の振幅情報のみ提供するが、高次の統計量は位相情報も提供する。したがって、一般に、高次統計量を含む最適化は、非最小位相系を導く可能性がある。そこで、１−Ａ_ｉ（ｚ）が最小位相系であるという拘束条件から、ａの最適化においては式（１９）の最適化問題を解く。

拘束条件
［１］ ｇ＝ｇ＾^（ｒ）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

Ｃ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））は式（２０）で与えられる。

ここで、行列Ａは式（９）および式（１０）に示すように上三角行列でその対角成分がすべて１であるから、ｌｏｇ ｄｅｔ Ａ＝０である。式（２０）を式（１２）に代入することで式（２１）の関係を得る。

よって、式（１９）は、式（２２）の最適化問題と等価である。式（２２）は、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるａを求めよ、と云っている。

拘束条件
［１］ ｇ＝ｇ＾^（ｒ）
［２］ １−Ａ_ｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。

式（２２）で表される最適化問題を解くことは、観測信号にｇ＾^（ｒ）で与えられる逆フィルタを適用して得られるアドホック信号に対して、各フレームにおいて線形予測分析を行うことと等価であり、必ず最小位相予測誤差フィルタを得ることができる。線形予測分析に関しては、上記参考文献１を参照されたい。

なお、式（２２）では、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大きい。この場合、対数関数は単調増加関数であるから、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めることができる。

１．９ ｇの最適化
本発明では、式（１８）の最適化を以下のように行う。
前述したとおり、Ｃ（ｄ_１（１），…，ｄ_Ｆ（Ｗ））は｛ｄ_ｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝の相関の度合いに関わる指標であるが、（ｒ＋１）回目のａの最適化において最小化されているため、Σ_ｉ＝１ ^ＦΣ_ｎ＝１ ^ＷＪ（ｄ_ｉ（ｎ））に比べて無視できる。そこでｇの最適化においては、式（２３）の最適化問題を解く。

拘束条件
［１］ ａ＝ａ＾^{（ｒ＋１）}
［２］ ‖ｇ‖＝１

Ｊ（ｄ_ｉ（ｎ））は、［条件２］に基づいて、式（２４）によって近似できる。この詳細は参考文献２を参照されたい。ただし、確率変数Ｕについて、κ_４（Ｕ）はＵの尖度（４次キュムラント）を表す。式（２４）の右辺を第ｉフレームにおける正規化尖度という。
（参考文献２） A.Hyvarinen, J.Karhunen, E.Oja, "INDEPENDENT COMPONENT ANALYSIS", John Wiley & Sons, Inc. 2001.

［条件２］から音声信号のイノベーションの尖度は正であるため、κ_４（ｄ_ｉ（ｎ））／σ（ｄ_ｉ（ｎ））^４は正である。従って、式（２３）の最適化問題は、式（２５）の最適化問題に帰着する。σ（ｄ_ｉ（ｎ）），κ_４（ｄ_ｉ（ｎ））は、［条件１］で述べた音声信号の局所的な定常性に基づいて、各フレーム内のサンプルから計算される。式（２６）では、１／Ｗを付加しているが、これは後の計算の便宜に過ぎず、式（２５）でｇの最適解を求めるにあたり影響を及ぼすものではない。式（２５）および式（２６）から、正規化尖度の全フレームでの総和が最大となるときのｇが、ｇ＾^{（ｒ＋１）}となる。これらは、第ｉフレームにおける正規化尖度を全フレームで加算した値が最大となるｇを求めよ、と云っている。

拘束条件
［１］ ａ＝ａ＾^{（ｒ＋１）}
［２］ ‖ｇ‖＝１

式（２５）に従ってｇの最適解を求めるには、Ｑをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（２７）で表される更新則に従って求められる。ｇ′をｇ′のノルムで除しているのは上記拘束条件[２]を課すためである。η（ｕ）は学習率を表す。ｕは、ｇの最適化における更新回数を表す。

式（２７）において、∇Ｑ_ｇは式（２８）および式（２９）で与えられる。

式（２９）において、ｄ_ｉ（ｎ）は式（３０）で、ｖ_ｍｉ（ｎ）は式（３１）および式（３２）で与えられる。ｘ_ｍｉ（ｎ）は、ｍ番目のマイクロホンで観測された観測信号のｉ番目のフレームの信号である。

§２ 二次統計量に基づく信号歪み除去
上述の従来的手法の信号歪み除去方法は、比較的長時間の観測信号（例えば２０秒程度である。）を要する。これは、一般に、正規化尖度のような高次統計量を計算するためには大量の観測信号のサンプルが必要となるからである。しかし、実際にはそうした長時間の観測信号を利用できない場合が多い。このため、従来的手法の信号歪み除去方法の適用
分野は極めて限られていた。
また高次統計量の計算は比較的複雑であるため、従来的手法の信号歪み除去方法では装置の構成が複雑になりやすい。
そこで、観測信号がより短時間（例えば３秒から５秒程度である。）の場合にも有効であり、かつ計算が従来に比して容易な信号歪み除去の原理を説明する。この原理は、信号の二次統計量のみを用いるものであり、§１で説明した本発明の基本原理より派生する。

２．１ 二次統計量に基づく信号歪み除去の原理
二次統計量に基づく信号歪み除去では、上述の３つの条件に、次の２つの条件を設定する。
[条件４] Ｍ≧２である。すなわち、複数本のマイクロホンを用いる。
[条件５] Ｈ_ｍ＝｛ｈ_ｍ（ｋ）｝_ｋ＝０ ^Ｋは相異なるマイクロホンの間で共通の零点を持たない。

上記の式（１６）の最適化問題では、高次の統計量に関する値であるネゲントロピーＪおよび確率変数間の相関の度合いを示す指標Ｃを含む値を最小化するｇおよびａを求めた。
確率変数間の相関の度合いを示す指標Ｃは、二次の統計量で規定される。そこで、解くべき最適化問題を式（３３）で定式化する。

式（２１）を参酌すれば、式（３３）の最適化問題は、式（３４）の最適化問題に転化される。式（３４）は、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるｇおよびａを求めよ、と云っている。

ところで、上記の[条件４]および[条件５]が成立する場合、多チャンネルの観測信号は、音源からの原信号によって駆動されるＡＲ系として捉えることができる（参考文献３参照）。このことは、逆フィルタＧの先頭タップを式（３５）のように固定できることを意味する。但し、ｍ＝１に相当するマイクロホンは、最も音源に近いマイクロホンである。（参考文献３） K. Aded-Meraim, E. Moulines, and P. Loubaton. Prediction error method for second-order blind identification. IEEE Trans. Signal Processing, Vol.
45, No.3, pp. 694-705, 1997.

式（３４）および式（３５）で規定されるｇを係数とする逆フィルタＧを、式（６）に従って観測信号ｘ（ｔ）に適用することで伝達特性が除去された復元信号ｙ（ｔ）を得る。

２．２ ａの最適化
式（３４）について、交代変数の方法により、ｇとａを最適化する。
逆フィルタの係数ｇ_ｍ（ｋ）を固定した状態で予測誤差フィルタの係数ａ_ｉ（ｋ）に関して式（３４）の損失関数を最小化する。
このとき、次の２点に注意する。１点目は、ｇ＝［ｇ_１ ^Ｔ，…，ｇ_Ｍ ^Ｔ］^Ｔは固定されているので、逆フィルタＧの出力である復元信号ｙ（ｔ）は予測誤差フィルタの最適化において不変であるということである。２点目は、ｉ番目のフレームにおける予測誤差フィルタの係数ａ_ｉ（１），…，ａ_ｉ（Ｐ）は、ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）にのみ寄与するということである。

従って、各々のフレームにおいて、Σ_ｎ＝１ ^Ｗlogσ（ｄ_ｉ（ｎ））^２を最小化するように、予測誤差フィルタの係数ａ_ｉ（１），…，ａ_ｉ（Ｐ）を推定すればよい。[条件２]から第ｉフレームのイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散はフレーム内で定常であるから、Σ_ｎ＝１ ^Ｗlogσ（ｄ_ｉ（ｎ））^２の最小化は、Ｗ＊σ（ｄ_ｉ（ｎ））^２の最小化と等価である。記号＊は乗算を表す。分散σ（ｄ_ｉ（ｎ））^２は、＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗとして計算できる。但し、＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗは、１フレーム分のイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を使って計算したｄ_ｉ（ｎ）の２乗平均を表す。結局、係数ａ_ｉ（ｋ）は、Ｗ＊＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗ、つまりｄ_ｉ（ｎ）の２乗総和が最小になるときのａ_ｉ（ｋ）として推定される。このような係数ａ_ｉ（ｋ）は、線形予測分析の手法によって計算される。

なお、ここでは、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めるとして説明したが、これに限定する趣旨ではない。既述のとおり、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾^{（ｒ＋１）}として求めることができる。

２．３ ｇの最適化
予測誤差フィルタの係数ａ_ｉ（ｋ）を固定した状態で逆フィルタの係数ｇ_ｍ（ｋ）に関して式（３４）の損失関数を最小化する。
逆フィルタの係数ｇ_ｍ（ｋ）に関する損失関数の最小化には、勾配法を用いる。[条件２]を用いると、式（３４）の最適化問題は、式（３６）の最適化問題に転化される。

式（３６）に従ってｇの最適解を求めるには、Σ_ｉ＝１ ^Ｆlog＜ｄ_ｉ（ｎ）^２＞_ｎ＝１
^Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（３７）で表される更新則に従って求められる。δは学習率を表す。１≦ｍ≦Ｍ、１≦ｋ≦Ｌである。なお、式（３７）では、式（３５）の条件から‖ｇ‖＝１なる拘束条件が課されないことに留意しなければならない。同様に式（３５）の条件からｋが取る値の範囲は１≦ｋ≦Ｌである。

上記式（２９）あるいは上記非特許文献１に記載された式（３）と比べて明らかなように、式（３７）の右辺第二項は、二次の統計量で表されており、この計算に高次の統計量を必要としない。このため、高次統計量を計算するには短い時間の観測信号の場合にも有効であり、計算自体も容易である。

なお、式（３６）では、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのｇをｇ＾として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大である。この場合、対数関数は単調増加関数であるから、第ｉフレームそれぞれにおけるイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散の全フレームでの総和が最小となるときのｇをｇ＾として求めることができる。なお、この場合には式（３７）で示した更新則は適用できなくなり、改めてΣ_ｉ＝１ ^Ｆ＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この結果得られる更新則はＩＣＡと同様の考え方で定式化できるから略する。

§３ プリ・ホワイトニング
本発明の信号歪み除去には、プリ・ホワイトニングを適用することができる。観測信号をプリ・ホワイトニングすることで、最適化計算の安定化、とくにフィルタ係数の推定値の高速な収束が可能となる。
各マイクロホンで得られた観測信号系列全体｛ｘ_ｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆ_ｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。
式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号ｘ_ｍ（ｔ）に適用する。ｗ_ｍ（ｔ）は、ｍ番目のマイクロホンで得た観測信号ｘ_ｍ（ｔ）を白色化した信号を表す。

このとき、式（３１）および式（３８）は式（４０）に、式（３２）は式（４１）に、変更すればよい。

§４ 実施形態
以下、本発明の実施形態を図面を参照して説明する。本発明の実施形態として、後述の各実施形態に限定するものではなく、各セクションで説明した原理を実現する実施形態であればよい。

＜第１実施形態＞
本発明の第１実施形態を実施する場合、以下の手順に従ってセンサで得た観測信号を処理する。ここでは、実施形態を具体的に説明する観点から信号として音声信号を例に挙げて説明する。
なお、第１実施形態の説明に先立ち、観測信号およびフレーム化処理について概説する。

((観測信号))
図示しないセンサ（例えばマイクロホン）によって得られたアナログ信号（このアナログ信号には伝達特性に由来する歪みが重畳されている。）は、例えば８,０００Ｈｚのサンプリングレートでサンプリングされ、適宜量子化された離散信号に変換される。以下、この離散信号を観測信号ということにする。アナログ信号から観測信号へのＡ／Ｄ変換などを実行するために必要となる構成要素（手段）は、いずれも公知技術の常套手段によって達成されるから、説明および図示を略する。

((フレーム化処理))
図示しない信号フレーム化手段が、離散信号から、時間軸方向に一定時間幅でフレームの始点を移動させながら、所定時間長の離散信号を切り出す。例えば２００サンプル点（８,０００Ｈｚ×２５ｍｓ）長の離散信号を、８０サンプル点（８,０００Ｈｚ×１０ｍｓ）ずつ始点を移動させながら切り出す。切り出された信号は、離散信号に公知の窓関数（例えば、ハミング窓、ガウス窓、方形窓など）が適用される。窓関数の適用によるフレーム化は公知の常套手段によって達成される。

本発明の第１実施形態である信号歪み除去装置（１）をコンピュータ（汎用機）で実現する場合のハードウェア構成例を説明する。
図２に例示するように、信号歪み除去装置（１）は、キーボード、ポインティングデバイスなどが接続可能な入力部（１１）、液晶ディスプレイ、ＣＲＴ（Cathode Ray Tube）ディスプレイなどが接続可能な出力部（１２）、信号歪み除去装置（１）外部に通信可能な通信装置（例えば通信ケーブル、ＬＡＮカード、ルータ、モデムなど）が接続可能な通信部（１３）、ＤＳＰ（Digital Signal Processor）（１４）〔ＣＰＵ（Central Processing Unit）でも良い。またキャッシュメモリやレジスタ（１９）などを備えていてもよい。〕、メモリであるＲＡＭ（１５）、ＲＯＭ（１６）やハードディスク、光ディスク、半導体メモリなどである外部記憶装置（１７）並びにこれらの入力部（１１）、出力部（１２）、通信部（１３）、ＤＳＰ（１４）、ＲＡＭ（１５）、ＲＯＭ（１６）、外部記憶装置（１７）間のデータのやり取りが可能なように接続するバス（１８）を有している。
また必要に応じて、信号歪み除去装置（１）に、ＣＤ−ＲＯＭ（Compact Disc Read Only
Memory）、ＤＶＤ（Digital Versatile Disc）などの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。

信号歪み除去装置（１）の外部記憶装置（１７）には、信号歪み除去のためのプログラムおよびこのプログラムの処理において必要となるデータ（観測信号）などが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶され、他のプログラムの処理に供されるときに、ＲＡＭや外部記憶装置などから読み込まれる。

より具体的には、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭなど〕には、観測信号に逆フィルタを適用する処理のためのプログラム、観測信号に逆フィルタを適用して得られた信号から予測誤差フィルタを求める処理のためのプログラム、予測誤差フィルタから逆フィルタを求める処理のためのプログラム、およびこれらのプログラムの処理において必要となるデータ（フレーム単位の観測信号など）が記憶されている。その他、これらのプログラムに基づく処理を制御するための制御プログラムも適宜に保存しておく。

第１実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、制御部）を実現することで、信号歪み除去が実現される。

そこで次に、図３〜図５を参照して、信号歪み除去装置（１）における信号歪み除去処理の流れを順次説明する。
大まかな処理の手順は、（ａ）観測信号ｘ（ｔ）に対して逆フィルタを適用した信号（以下、アドホック信号という。）を求め、（ｂ）アドホック信号から予測誤差フィルタを求め、（ｃ）この予測誤差フィルタから逆フィルタを求め、（ｄ）前記（ａ）、（ｂ）、（ｃ）の処理を繰り返して最適な逆フィルタを求め、（ｅ）最適化された逆フィルタを観測信号に対して適用した信号を復元信号ｙ（ｔ）として得る。
（ｂ）は上述のａの最適化に相当し、（ｃ）は上述のｇの最適化に相当し、（ｄ）は、式（１７）および式（１８）に相当する。（ｄ）の処理の繰り返し回数は予め定めた回数Ｒ_１とする。つまり、１≦ｒ≦Ｒ_１とする。また、（ｃ）の処理でｇを最適化する更新則の更新回数は予め定めた回数Ｒ_２とする。つまり、１≦ｕ≦Ｒ_２とする。（ｄ）の処理、つまり（ａ）、（ｂ）、（ｃ）の一連の処理を１回行うたびに、更新則によるＲ_２回の更新が行なわれる。実施形態では、回数Ｒ_１は、予め定めた回数とするが、これに限定されず、例えば、ｒ回目のｇを算出したときの式（２６）のＱの値とｒ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。同様に、回数Ｒ_２は、予め定めた回数とするが、これに限定されず、例えば、ｕ回目のｇを算出したときの式（２６）のＱの値とｕ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。

（ステップＳ１００）
逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタを観測信号ｘ（ｔ）＝［ｘ_１（ｔ），…，ｘ_ｍ（ｔ），…，ｘ_Ｍ（ｔ）］^Ｔに適用することで、アドホック信号ｙ（ｔ）を求める。アドホック信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、こ
こでは後述のＲ_１回の処理を経て求められた復元信号ではないことを明示するためアドホック信号と呼称する。ここでｔは、全てのサンプル番号を示し、１≦ｔ≦Ｎとする。Ｎは全サンプル数である。第１実施形態では、マイクロホンの個数Ｍは１以上であればよい。

逆フィルタの係数列｛ｇ_ｍ（ｋ）；０≦ｋ≦Ｌ｝として、繰り返し回数Ｒ_１の初回には予め定められた初期値を、２回目以降には後述する逆フィルタ計算部（１３）によって求められた逆フィルタｇ＾^{（ｒ＋１）}を用いる。

（ステップＳ１０１）
予測誤差フィルタ計算部（１５）は、フレーム化処理を行うフレーム化処理部（１５１）とフレーム予測誤差フィルタ計算部（１５２）によって構成される。そして、フレーム予測誤差フィルタ計算部（１５２）は、第ｉフレームのアドホック信号から予測誤差フィルタを求める第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）からなる。ただし、ｉは、１≦ｉ≦Ｆを満たす整数である。

フレーム化処理部（１５１）は、逆フィルタ適用部（１４）で求められたアドホック信号｛ｙ（ｔ）；１≦ｔ≦Ｎ｝をフレーム化処理する。フレーム化処理は、例えば式（４３）のように、Ｗ点分を切り出す窓関数をＷ点ずつシフトさせて適用することにより行う。｛ｙ_ｉ（ｎ）；１≦ｎ≦Ｗ｝はｉ番目のフレームに含まれるアドホック信号列を表す。

そして、第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）は、式（２２）に従って、第ｉフレームのアドホック信号列｛ｙ_ｉ（ｎ）；１≦ｎ≦Ｗ｝に対してＰ次の線形予測分析を行い、予測誤差フィルタの係数列｛ａ_ｉ（ｋ）；１≦ｋ≦Ｐ｝を計算する。この算出方法は、上記参考文献１を参照されたい。ここで得られたａ_１（１），…，ａ_１（Ｐ），…，ａ_ｉ（１），…，ａ_ｉ（Ｐ），…，ａ_Ｆ（１），…，ａ_Ｆ（Ｐ）は、式（２２）のａ＾^{（ｒ＋１）}を与える。

（ステップＳ１０２）
逆フィルタ計算部（１３）の機能構成例を、図４を参照して説明する。逆フィルタ計算部（１３）は勾配計算部（１３１）、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。更に、勾配計算部（１３１）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）とを備えて構成される。ここで更新用逆フィルタは、式（２７）のｇ^〈ｕ〉に相当する。

第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘ_ｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘ_ｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖ_ｍｉ（ｎ）を計算する（式（
３１）を参照）。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。

第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｙ_ｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用してイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。爾後、第２の予測誤差フィルタ適用部（１３１２）は、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）を入力とする。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。

勾配ベクトル計算部（１３１３）は、信号ｖ_ｍｉ（ｎ）とイノベーション推定値ｄ_ｉ（ｎ）とを用いて現在の更新用逆フィルタｇ^〈ｕ〉の勾配ベクトル∇Ｑ_ｇを計算する（式（２８）および式（２９）を参照）。有限個のサンプルｖ_ｍｉ（ｎ）およびｄ_ｉ（ｎ）を用いて式（２９）を演算するときは、期待値Ｅをサンプルから求めればよい。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。

逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇ^〈ｕ〉、学習率η（ｕ）、勾配ベクトル∇Ｑ_ｇを用いて、式（２７）に従って、ｕ＋１回目の更新用逆フィルタｇ^{〈ｕ＋１〉}を求める。式（２７）は、求められたｇ^{〈ｕ＋１〉}を新たなｇ^〈ｕ〉と見立てて更新を行なうことを意味する。

更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇ^{〈ｕ＋１〉}、つまり新たなｇ^〈ｕ〉および観測信号ｘ（ｔ）を用いて、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。つまり、式（４２）のｇ_ｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。なお、更新用逆フィルタ適用後信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、ここでは後述のＲ_１回の処理を経て求められた復元信号ではなく、更新則を行なうために算出される信号であることを明示するため更新用逆フィルタ適用後信号と呼称する。

制御部（６００）の制御によってＲ_２回の更新が行なわれた結果として得られたｇ^{〈Ｒ２＋１〉}は、式（２５）のｇ＾^{（ｒ＋１）}に相当する。上付き文字のＲ２は、Ｒ_２である。逆フィルタ計算部（１３）は、ｇ＾^{（ｒ＋１）}を出力する。

制御部（５００）の制御によって、上述の一連の処理を１回行うごとにｒに１を加算してｒがＲ_１に等しくなるまで、つまり上述の一連の処理をＲ_１回繰り返すことで（ステップＳ１０３）、ｇ＾^{（Ｒ１＋１）}を得る。上付き文字のＲ１は、Ｒ_１である。このｇ＾^{（Ｒ１＋１）}が、式（１６）の最適解とされる。そこで、ｇ＾^{（Ｒ１＋１）}を得た段階で、逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタｇ＾^{（Ｒ１＋１）}を観測信号ｘ（ｔ）＝［ｘ_１（ｔ），…，ｘ_Ｍ（ｔ）］^Ｔに適用することで、復元信号ｙ（ｔ）を得ることができる（ステップＳ１０４）。

＜第２実施形態＞
第２実施形態は、第１実施形態の変形例に相当する。具体的には、§３で述べたプリ・ホワイトニングを行なう形態である。そこで、第１実施形態と異なる部分について図６および図７を参照して説明を加える。なお、プリ・ホワイトニングは観測信号に対して行なうプリ・プロセスであるから、ここで説明するプリ・ホワイトニングを行なう形態は、後述の第３実施形態にも適用可能である。
第２実施形態では、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭ
など〕に、白色化フィルタを求める処理のためのプログラム、白色化フィルタを観測信号に適用する処理のためのプログラムも記憶されている。
第２実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、白色化フィルタ計算部、白色化フィルタ適用部）を実現することで、信号歪み除去が実現される。

（ステップＳ１００ａ）
白色化フィルタ計算部（１１）は、各マイクロホンで得られた観測信号全体｛ｘ_ｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆ_ｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。この計算は線形予測分析と同じであり上記参考文献１を参照されたい。白色化フィルタの係数は、白色化フィルタ適用部（１２）の入力となる。

（ステップＳ１００ｂ）
白色化フィルタ適用部（１２）は、式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号に適用して、白色化信号ｗ_ｍ（ｔ）を得る。既述のとおり、式（３１）は式（４０）に変更すればよいので、第１実施形態において、逆フィルタ計算部（１３）、とくに第１の予測誤差フィルタ適用部（１３１１）による処理を式（３１）ではなく式（４０）による計算処理に改めればよい。また、第１実施形態において、逆フィルタ適用部（１４）による処理を、式（４２）ではなく式（４４）による計算処理に改めればよい。ステップＳ１００ｂの処理の後、第１実施形態のステップＳ１００〜Ｓ１０４の処理を行うが、これらの処理では第１実施形態の各処理における観測信号をステップＳ１００ｂの処理で得られた白色化信号に読み替えて第１実施形態と同様の処理を行う。このことを明示するため、図７では、第１実施形態のステップＳ１００〜Ｓ１０４の各処理に相当する処理を示す符号に記号′を付している。

＜実施例１＞
発明者らは第２実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、マイクロホンの数Ｍ＝４、白色化フィルタの次数Ｘ＝５００、逆フィルタの次数Ｌ＝１０００、窓関数の切出しサンプル数（１フレームのサンプル数）Ｗ＝２００、予測誤差フィルタの次数Ｐ＝１６、繰り返し回数Ｒ_１＝１０、逆フィルタ計算部の更新回数Ｒ_２＝２０とした。学習率η（ｕ）は、初期値を０．０５に設定し、もし式（２７）によって式（２６）の値が減少するならば、式（２６）の値が増大するまでη（ｕ）の値を再帰的に半減した。図６に示した逆フィルタ適用部（１４）へ入力する初期逆フィルタは、式（４５）のように設定した。

本発明の第２実施形態の効果を、信号歪み除去の指標としてＤ_５０値（インパルス応答
の全エネルギーに対する初期の５０ｍｓｅｃまでのエネルギーの比）を用いて評価した。連続発話データベースから男女各一名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
図８は、男声および女声について観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ_１ （図６に示す逆フィルタ適用部（１４）と、予測誤差フィルタ計算部（１５）と、逆フィルタ計算部（１３）を一巡する処理を実行して逆フィルタを求める回数）とＤ_５０値の関係を示している。いずれの場合においても、繰り返し回数を増加させるとＤ_５０値が向上しており、繰り返し処理の効果が顕著に見て取れる。特に観測信号長が５〜１０秒程度の比較的短い長さであっても、繰り返し処理によってＤ_５０値が大きく向上したことが分かる。

また、本発明の第２実施形態の効果を、音声スペクトグラムの比較から検証した。
図９Ａは１分間の観測信号を用いて得られた残響を含まない音声（原音声）のスペクトログラムの抜粋、図９Ｂは１分間の観測信号を用いて得られた残響を含む音声（観測音声）のスペクトログラムの抜粋、図９Ｃは１分間の観測信号を用いて得られた残響除去後の音声（復元音声）のスペクトログラムの抜粋を示している。図９Ａと図９Ｃとの対比および図９Ｂと図９Ｃとの対比から、観測信号に含まれる残響が抑制され、原音声固有の特徴である調波構造やフォルマント構造が回復されたことが分かる。

また、本発明の第２実施形態の効果を、ＬＰＣスペクトル歪みを用いて検証した。
図１０Ｂは原音声の波形、図１０Ａは原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）および原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）を示している。図１０Ａおよび図１０Ｂの各横軸は秒単位の時間を表し、両図でタイムスケールを揃えている。図１０Ｂの縦軸は振幅値を表している。但し、原信号の相対的な大きさが分かればよいので、この趣旨から図１０Ｂの縦軸では単位を明示していない。図１０Ａの縦軸は、ＬＰＣスペクトル歪みＳＤ（ｄＢ）を表している。
図１０Ａから、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりも常に小さいことがわかる。なお、観測音声ではＬＰＣスペクトル歪みの平均が５．３９ｄＢ、分散が４．２０ｄＢであったのに対して、復元音声では平均が２．３８ｄＢ、分散が２．００ｄＢであった。
また、図１０Ａと図１０Ｂとの対比から、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）が大きな値を示す区間（例えば約１．０秒〜約１．２秒の区間を参照）は、原音声の波形の振幅値がほぼ０であることがわかる。実際、この区間では発声がなく無音区間である。このため、実際に知覚される歪みはかなり小さくなっていた。つまり、発声区間における原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりもかなり小さく、このため原音声のスペクトルを高い精度で復元できたことが結論付けられる。

＜第３実施形態＞
第３実施形態は、第１実施形態の変形例に相当する。具体的には、§２で述べた二次統計量に基づく信号歪み除去処理を行なう形態である。そこで、第１実施形態と異なる部分について図１１および図１２を参照して説明を加える。但し、第３実施形態では、マイクロホンの個数Ｍは２以上とする。

ステップＳ１００の処理およびステップＳ１０１の処理は、第１実施形態と同じである。

ステップＳ１０１の処理に続いて、ステップＳ１０２ａの処理を行う。
第３実施形態に係る逆フィルタ計算部（１３）の機能構成例を、図１１を参照して説明
する。
逆フィルタ計算部（１３）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）と、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。ここで更新用逆フィルタは、式（３７）のｇ_ｍ（ｋ）に相当する。

第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘ_ｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘ_ｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖ_ｍｉ（ｎ）を計算する（式（３８）を参照）。具体的には、フレーム化処理部（４０２Ｂ）が、入力された観測信号ｘ_ｍ（ｔ）に対してフレーム化処理を行い、観測信号ｘ_ｍ（ｔ）のｉ番目のフレームの信号ｘ_ｍｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０４ｉ）が信号ｘ_ｍｉ（ｎ）を入力として、式（３８）に従って信号ｖ_ｍｉ（ｎ）を出力する。但し、１≦ｉ≦Ｆである。

第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｙ _ｉ （ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａ_ｉ（ｋ）を適用してイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。具体的には、フレーム化処理部（４０２Ａ）が、初期値の場合を除き、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）に対してフレーム化処理を行い、ｉ番目のフレームの信号ｙ_ｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０３ｉ）が信号ｙ_ｉ（ｎ）を入力として、式（３０）に従ってイノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）を出力する。但し、１≦ｉ≦Ｆである。

勾配ベクトル計算部（１３１３）は、信号ｖ_ｍｉ（ｎ）とイノベーション推定値ｄ_ｉ（ｎ）とを用いて現在の更新用逆フィルタｇ_ｍ（ｋ）の勾配ベクトルを計算する（式（３７）の右辺第二項を参照）。具体的には、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、相互相関計算部（４０５ｉ）は信号ｖ_ｍｉ（ｎ）とイノベーション推定値ｄ _ｉ （ｎ）との相互相関〈ｄ_ｉ（ｎ）ｖ_ｍｉ（ｎ−ｋ）〉_ｎ＝１ ^Ｗを計算する。また、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、分散計算部（４０６ｉ）は、イノベーション推定値ｄ_ｉ（１），…，ｄ_ｉ（Ｗ）の分散〈ｄ_ｉ（ｎ）^２〉_ｎ＝１ ^Ｗを求める。各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、除算部（４０７ｉ）は、〈ｄ_ｉ（ｎ）ｖ_ｍｉ（ｎ−ｋ）〉_ｎ＝１ ^Ｗ／〈ｄ_ｉ（ｎ）^２〉_ｎ＝１ ^Ｗを求める。加算部（４０８）は、除算部（４０７１）〜（４０７Ｆ）の出力の全フレームに亘る総和、つまり式（３７）の右辺第二項を求める。

逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇ_ｍ（ｋ）、学習率δ、勾配ベクトルを用いて、式（３７）に従って、ｕ＋１回目の更新用逆フィルタｇ_ｍ（ｋ）′を求める。式（３７）は、求められたｇ_ｍ（ｋ）′を新たなｇ_ｍ（ｋ）と見立てて更新を行なうことを意味する。

更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇ_ｍ（ｋ）′、つまり新たなｇ_ｍ（ｋ）および観測信号ｘ（ｔ）を用いて、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。具体的には、式（４２）のｇ_ｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新
用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。

ステップＳ１０２ａの処理に続いて、ステップＳ１０３およびステップＳ１０４の処理を行うが、第１実施形態と同じであるから説明を略する。

＜実施例２＞
発明者らは第３実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、Ｍ＝４，Ｌ＝１０００，Ｗ＝２００，Ｐ＝１６，Ｒ_１＝６，Ｒ_２＝５０とした。学習率δは、初期値を０．０５に設定し、Σ_ｉ＝１ ^Ｆlog＜ｄ_ｉ（ｎ）^２＞_ｎ＝１ ^Ｗの値が増加するならば、Σ _ｉ＝１ ^Ｆ log＜ｄ _ｉ （ｎ） ^２ ＞ _ｎ＝１ ^Ｗ が減少するまで、学習率δの値を順次半減させた。逆フィルタの初期推定値は、ｇ_ｍ（ｋ）＝０，１≦ｍ≦Ｍ，１≦ｋ≦Ｌとして設定した。

本発明の第３実施形態の効果を、音声明瞭度を表すＲＡＳＴＩ（参考文献５を参照）を残響除去の指標として評価した。連続発話データベースから男女各五名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
（参考文献５） H. kuttruff. Room acoustics. Elsevier Applied Science, third edition, P.237 1991.

図１３は、Ｎを３秒、４秒、５秒、１０秒とする各観測信号のＲＡＳＴＩ値を表示したものである。図１３に示すように、観測信号が３〜５秒のように短時間の場合でも、高い残響除去性能を示していることが分かる。
図１４は、残響除去前後におけるエネルギー減衰曲線の例である。直接音が到達してから５０ミリ秒後の反射音のエネルギーが１５ｄＢ低減されていることが分かる。

本発明は、様々な信号処理システムの性能向上に寄与する要素技術であるところ、例えば音声認識システム、テレビ会議システム、補聴器、音楽情報処理システム等に利用することができる。

本発明の原理を説明するためのモデル機構を表したブロック線図。 第１実施形態に係る信号歪み除去装置（１）のハードウェア構成例を示す図。 第１実施形態に係る信号歪み除去装置（１）の機能構成例を示す機能ブロック図。 信号歪み除去装置（１）の逆フィルタ計算部（１３）の機能構成例を示す機能ブロック図。 第１実施形態における信号歪み除去処理の流れを示す処理フロー図。 第２実施形態に係る信号歪み除去装置（１）の機能構成例を示す機能ブロック図。 第２実施形態における信号歪み除去処理の流れを示す処理フロー図。 観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ_１とＤ_５０値の関係を示す図。 Ａは残響を含まない音声のスペクトログラム、Ｂは残響を含む音声のスペクトログラム、Ｃは残響除去後の音声のスペクトログラム。 Ａは残響除去音声のＬＰＣスペクトル歪みの時間変動を説明するためのグラフ、Ｂは対応する区間における原音声信号の抜粋。 第３実施形態に係る信号歪み除去装置（１）の逆フィルタ計算部（１３）の機能構成例を示す機能ブロック図。 第３実施形態における信号歪み除去処理の流れを示す処理フロー図。 ３秒、４秒、５秒、１０秒の各観測信号のＲＡＳＴＩ値を表示した図。 残響除去前後におけるエネルギー減衰曲線の例を示した図。 従来技術を説明するための機能ブロック図。

Claims (14)

1. 観測信号の信号歪みを除去して復元信号を得る信号歪み除去装置であって、
   所定の繰り返し終了条件を満たした場合には、上記観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記観測信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記観測信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、
   上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、
   上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、
   上記繰り返し終了条件を満たすまで上記逆フィルタ適用手段、上記予測誤差フィルタ計算手段、上記逆フィルタ計算手段を繰り返し実行させる制御手段と、
   を備えた信号歪み除去装置。
2. 上記予測誤差フィルタ計算手段は、
   上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
   上記逆フィルタ計算手段は、
   上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタとして求め、この逆フィルタを出力するものである
   ことを特徴とする請求項１に記載の信号歪み除去装置。
3. 上記予測誤差フィルタ計算手段は、
   上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
   上記逆フィルタ計算手段は、
   上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタとして求め、この逆フィルタを出力するものである
   ことを特徴とする請求項１に記載の信号歪み除去装置。
4. 観測信号の信号歪みを除去して復元信号を得る信号歪み除去装置であって、
   上記観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算手段と、
   上記白色化フィルタを上記観測信号に適用して白色化信号を出力する白色化フィルタ適用手段と、
   所定の繰り返し終了条件を満たした場合には、上記白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記白色化信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記白色化信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、
   上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、
   上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、
   上記繰り返し終了条件を満たすまで上記逆フィルタ適用手段、上記予測誤差フィルタ計算手段、上記逆フィルタ計算手段を繰り返し実行させる制御手段と、
   を備えた信号歪み除去装置。
5. 上記繰り返し終了条件は、
   繰り返し回数がＲ_１回（但しＲ_１は、Ｒ_１≧１を満たす整数である。）となることである
   ことを特徴とする請求項１から請求項４のいずれかに記載の信号歪み除去装置。
6. 上記観測信号は、信号歪みを含む音声信号である
   ことを特徴とする請求項１から請求項５のいずれかに記載の信号歪み除去装置。
7. 観測信号の信号歪みを除去して復元信号を得る信号歪み除去方法であって、
   逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、上記観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記観測信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記観測信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、
   予測誤差フィルタ計算手段が、上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、
   逆フィルタ計算手段が、上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、
   制御手段が、上記繰り返し終了条件を満たすまで上記逆フィルタ適用ステップ、上記予測誤差フィルタ計算ステップ、上記逆フィルタ計算ステップを繰り返し実行させる制御ステップと、
   を有する信号歪み除去方法。
8. 上記予測誤差フィルタ計算ステップは、
   上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
   上記逆フィルタ計算ステップは、
   上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタとして求め、この逆フィルタを出力するものである
   ことを特徴とする請求項７に記載の信号歪み除去方法。
9. 上記予測誤差フィルタ計算ステップは、
   上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
   上記逆フィルタ計算ステップは、
   上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタとして求め、この逆フィルタを出力するものである
   ことを特徴とする請求項７に記載の信号歪み除去方法。
10. 観測信号の信号歪みを除去して復元信号を得る信号歪み除去方法であって、
    白色化フィルタ計算手段が、上記観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算ステップと、
    白色化フィルタ適用手段が、上記白色化フィルタを上記観測信号に適用して白色化信号を出力する白色化フィルタ適用ステップと、
    逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、上記白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記白色化信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記白色化信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、
    予測誤差フィルタ計算手段が、上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、
    逆フィルタ計算手段が、上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、
    制御手段が、上記繰り返し終了条件を満たすまで上記逆フィルタ適用ステップ、上記予測誤差フィルタ計算ステップ、上記逆フィルタ計算ステップを繰り返し実行させる制御ステップと、
    を有する信号歪み除去方法。
11. 上記繰り返し終了条件は、
    繰り返し回数がＲ_１回（但しＲ_１は、Ｒ_１≧１を満たす整数である。）となることである
    ことを特徴とする請求項７から請求項１０のいずれかに記載の信号歪み除去方法。
12. 上記観測信号は、信号歪みを含む音声信号である
    ことを特徴とする請求項７から請求項１１のいずれかに記載の信号歪み除去方法。
13. 請求項１から請求項６のいずれかに記載された信号歪み除去装置としてコンピュータを機能させるための信号歪み除去プログラム。
14. 請求項１３に記載の信号歪み除去プログラムを記録した、コンピュータに読み取り可能な記録媒体。

Images