(57)【要約】  所定の繰り返し終了条件を満たした場合には、観測信号に逆フィルタを適用した結果を復元信号として出力し、繰り返し終了条件を満たさない場合には、観測信号に逆フィルタを適用した結果をアドホック信号として出力する逆フィルタ適用手段と、アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、各フレームのアドホック信号に対してこのフレームに対応する予測誤差フィルタを適用して得る各信号（イノベーション推定値）を結合した全フレームでのイノベーション推定値が、その全サンプル間で独立となる逆フィルタを求めてこれを出力する逆フィルタ計算手段と、繰り返し終了条件を満たすまで逆フィルタ適用手段、予測誤差フィルタ計算手段、逆フィルタ計算手段を繰り返し実行させる制御手段とを備えた信号歪み除去装置とする。
2007094463.tif 000002
【特許請求の範囲】
【請求項１】
観測信号の信号歪みを除去して復元信号を得る信号歪み除去装置であって、
所定の繰り返し終了条件を満たした場合には、上記観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記観測信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記観測信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、
上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、
上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、
上記繰り返し終了条件を満たすまで上記逆フィルタ適用手段、上記予測誤差フィルタ計算手段、上記逆フィルタ計算手段を繰り返し実行させる制御手段と、
を備えた信号歪み除去装置。
【請求項２】
上記予測誤差フィルタ計算手段は、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算手段は、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタとして求め、この逆フィルタを出力するものである
ことを特徴とする請求項１に記載の信号歪み除去装置。
【請求項３】
上記予測誤差フィルタ計算手段は、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算手段は、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタとして求め、この逆フィルタを出力するものである
ことを特徴とする請求項１に記載の信号歪み除去装置。
【請求項４】
観測信号の信号歪みを除去して復元信号を得る信号歪み除去装置であって、
上記観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算手段と、
上記白色化フィルタを上記観測信号に適用して白色化信号を出力する白色化フィルタ適用手段と、
所定の繰り返し終了条件を満たした場合には、上記白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記白色化信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記白色化信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用手段と、
上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算手段と、
上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算手段と、
上記繰り返し終了条件を満たすまで上記逆フィルタ適用手段、上記予測誤差フィルタ計算手段、上記逆フィルタ計算手段を繰り返し実行させる制御手段と、
を備えた信号歪み除去装置。
【請求項５】
上記繰り返し終了条件は、
繰り返し回数がＲ１回（但しＲ１は、Ｒ１≧１を満たす整数である。）となることである
ことを特徴とする請求項１から請求項４のいずれかに記載の信号歪み除去装置。
【請求項６】
上記観測信号は、信号歪みを含む音声信号である
ことを特徴とする請求項１から請求項５のいずれかに記載の信号歪み除去装置。
【請求項７】
観測信号の信号歪みを除去して復元信号を得る信号歪み除去方法であって、
逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、上記観測信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記観測信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記観測信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、
予測誤差フィルタ計算手段が、上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、
逆フィルタ計算手段が、上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、
制御手段が、上記繰り返し終了条件を満たすまで上記逆フィルタ適用ステップ、上記予測誤差フィルタ計算ステップ、上記逆フィルタ計算ステップを繰り返し実行させる制御ステップと、
を有する信号歪み除去方法。
【請求項８】
上記予測誤差フィルタ計算ステップは、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算ステップは、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタとして求め、この逆フィルタを出力するものである
ことを特徴とする請求項７に記載の信号歪み除去方法。
【請求項９】
上記予測誤差フィルタ計算ステップは、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算ステップは、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタを、上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタとして求め、この逆フィルタを出力するものである
ことを特徴とする請求項７に記載の信号歪み除去方法。
【請求項１０】
観測信号の信号歪みを除去して復元信号を得る信号歪み除去方法であって、
白色化フィルタ計算手段が、上記観測信号を線形予測分析して得た白色化フィルタを出力する白色化フィルタ計算ステップと、
白色化フィルタ適用手段が、上記白色化フィルタを上記観測信号に適用して白色化信号を出力する白色化フィルタ適用ステップと、
逆フィルタ適用手段が、所定の繰り返し終了条件を満たした場合には、上記白色化信号に適用するためのフィルタ（以下、逆フィルタという。）を、上記白色化信号に適用して、この結果を上記復元信号として出力し、上記繰り返し終了条件を満たさない場合には、上記白色化信号に上記逆フィルタを適用して、この結果をアドホック信号として出力する逆フィルタ適用ステップと、
予測誤差フィルタ計算手段が、上記アドホック信号をフレーム化して、各フレームのアドホック信号を線形予測分析して得た各フレームごとの予測誤差フィルタを出力する予測誤差フィルタ計算ステップと、
逆フィルタ計算手段が、上記各フレームのアドホック信号に対して当該フレームに対応する上記予測誤差フィルタを適用して得る各信号（以下、イノベーション推定値という。）を結合した全フレームでのイノベーション推定値（以下、イノベーション推定値系列という。）が、その全サンプル間で独立となる上記逆フィルタを求め、この逆フィルタを出力する逆フィルタ計算ステップと、
制御手段が、上記繰り返し終了条件を満たすまで上記逆フィルタ適用ステップ、上記予測誤差フィルタ計算ステップ、上記逆フィルタ計算ステップを繰り返し実行させる制御ステップと、
を有する信号歪み除去方法。
【請求項１１】
上記繰り返し終了条件は、
繰り返し回数がＲ１回（但しＲ１は、Ｒ１≧１を満たす整数である。）となることである
ことを特徴とする請求項７から請求項１０のいずれかに記載の信号歪み除去方法。
【請求項１２】
上記観測信号は、信号歪みを含む音声信号である
ことを特徴とする請求項７から請求項１１のいずれかに記載の信号歪み除去方法。
【請求項１３】
請求項１から請求項６のいずれかに記載された信号歪み除去装置としてコンピュータを機能させるための信号歪み除去プログラム。
【請求項１４】
請求項１３に記載の信号歪み除去プログラムを記録した、コンピュータに読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、信号の歪み（ひずみ：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.【００１７】
１．２  音声信号のモデル化
まず、原信号である音声信号ｓ（ｔ）を、以下の３つの条件を満足する信号としてモデル化する。【００１８】
［条件１］音声信号ｓ（ｔ）は、区分定常なＡＲ過程で生成される。
この[条件１]から、ＡＲ過程の次数をＰ、定常とみなせる区間長をＷサンプルとして音声信号ｓ（ｔ）をフレーム化すると、第ｉフレームの音声信号ｓｉ（ｎ）は、式（１）のように表される。式（２）は、第ｉフレームの音声信号ｓｉ（ｎ）のサンプルと、フレーム化前の音声信号ｓ（ｔ）のサンプルとの対応を示している。つまり、第ｉフレームのｎ番目のサンプルは、音声信号ｓ（ｔ）において、（ｉ−１）Ｗ＋ｎ番目のサンプルに相当する。式（１）および式（２）において、ｂｉ（ｋ）は線形予測係数、ｅｉ（ｎ）はイノベーションを表す。但し、１≦ｎ≦Ｗ、１≦ｔ≦Ｎ、Ｎは全サンプル数である。以下、特に断りの無い限り、パラメータｎは１フレームのサンプル番号を表し、パラメータｔは全てのサンプル番号を表す。また、全フレーム数はＦとする。
【数１】
2007094463.tif 000003
【００１９】
なお、第ｉフレームにおけるｎ番目のイノベーションｅｉ（ｎ）についても、第ｉフレームのｎ番目のイノベーションｅｉ（ｎ）と、フレーム化前の音声信号ｓ（ｔ）に対するイノベーションｅ（ｔ）との対応を示すことができる。この場合、第ｉフレームのｎ番目のイノベーションｅｉ（ｎ）は、イノベーションｅ（ｔ）において（ｉ−１）Ｗ＋ｎ番目のイノベーションに相当し、ｅｉ（ｎ）＝ｅ（（ｉ−１）Ｗ＋ｎ）が成り立つ。【００２０】
式（１）をｚ変換する。左辺のｚ変換をＳｉ（Ｚ）とし、右辺第二項のｚ変換をＥｉ（Ｚ）とし、Ｂｉ（ｚ）＝Σｋ＝１Ｐｂｉ（ｋ）ｚ−ｋとすれば、右辺第一項は、Ｂｉ（ｚ）Ｓｉ（Ｚ）となる。従って、式（１）のｚ変換は、（１−Ｂｉ（ｚ））Ｓｉ（Ｚ）＝Ｅｉ（Ｚ）である。なお、ｚ−１は時間領域では１タップ遅延素子に相当する。以降、時間領域信号（タップ重み係数）を小文字で、ｚ領域信号（伝達関数）を大文字でそれぞれ表す。１−Ｂｉ（ｚ）は最小位相性を満足しなければならず、『１−Ｂｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ』ことが要求される。【００２１】
［条件２］第ｉフレームに属するイノベーション系列ｅｉ（１），…，ｅｉ（Ｗ）は独立且つ同一分布に属する。イノベーション系列ｅｉ（１），…，ｅｉ（Ｗ）の確率分布の平均及び歪度（３次キュムラント）は０、尖度（４次キュムラント）は正である。さらに、異なるフレームｉ、ｊ〔ｉ≠ｊ〕に属するイノベーションｅｉ（ｎ）とｅｊ（ｎ′）同士も独立である。ただし、これらは必ずしも同一分布に属するとは限らない。
［条件３］予測誤差フィルタ１−Ｂｉ（ｚ）は、相異なるフレーム間で共通する零点をもたない。【００２２】
式（１）および式（２）から、音声信号ｓ（ｔ）は、式（３）のように表される。[・]はガウス記号である。
【数２】
2007094463.tif 000004
【００２３】
このとき、[条件２]は、『イノベーション過程ｅ（ｔ）は時間的に独立な信号である。また、その統計的性質（あるいは統計量）はフレーム内では定常である。』と表現できる。また、[条件３]は、『線形予測係数｛ｂｉ（ｋ）｝ｋ＝１Ｐは、時不変な極を持たない』と表現できる。【００２４】
１．３  観測信号のモデル化
次に、Ｍ個のマイクロホンで音声信号を観測して観測信号を得たときの観測信号をモデル化する。但しＭは、Ｍ≧１の整数である。
ｍ番目（１≦ｍ≦Ｍ）のマイクロホンで観測される残響を含む観測信号ｘｍ（ｔ）を、音源からｍ番目のマイクロホンに至る経路の伝達関数Ｈｍ（ｚ）のタップ重み係数｛ｈｍ（ｋ）；０≦ｋ≦Ｋ；Ｋはインパルス応答の持続時間とする。｝を用いて式（４）のようにモデル化する。ここでは、音声信号の場合の伝達特性の代表例として残響を挙げて、伝達特性を残響に言い換えて説明する。但し、伝達特性を残響に限定する趣旨ではない。
【数３】
2007094463.tif 000005
【００２５】
Ｍ個の観測信号についてまとめて表現すれば、式（５）のように表すことができる。但し、式（５）において、ｘ（ｔ）＝［ｘ１（ｔ），…，ｘＭ（ｔ）］Ｔ、ｈ（ｋ）＝［ｈ１（ｋ），…，ｈＭ（ｋ）］Ｔである。
【数４】
2007094463.tif 000006
【００２６】
１．４  信号歪み除去の原理
信号歪み除去後の復元信号ｙ（ｔ）は、多チャネル逆フィルタ｛Ｇｍ（ｚ）；１≦ｍ≦Ｍ｝のタップ重み係数｛ｇｍ（ｋ）；１≦ｍ≦Ｍ，０≦ｋ≦Ｌ；Ｌは逆フィルタの次数｝を用いて式（６）により計算される。本発明においては、逆フィルタ係数であるｇｍ（ｋ）を観測信号ｘ１（ｔ），…，ｘＭ（ｔ）のみから推定する。
【数５】
2007094463.tif 000007
【００２７】
１．５  本発明の基本原理
本発明の基本原理は、伝達関数｛Ｈｍ（ｚ）；１≦ｍ≦Ｍ｝の逆フィルタ｛Ｇｍ（ｚ）；１≦ｍ≦Ｍ｝とＡＲ系のフィルタ｛１／（１−Ｂｉ（ｚ））；１≦ｉ≦Ｆ｝の逆フィルタである予測誤差フィルタ｛１−Ａｉ（ｚ）；１≦ｉ≦Ｆ｝とを並行して推定することを主な特徴とする。【００２８】
この基本原理を説明するため、上述のモデル機構を組み込んだ系全体の構成線図を図１に示す。上述のモデル化に拠れば、原信号ｓ（ｔ）は、フレームごとのイノベーション系列ｅｉ（１），…，ｅｉ（Ｗ）に対してＡＲ系のフィルタ１／（１−Ｂｉ（ｚ））を適用して得られる信号ｓｉ（ｎ）のフレーム結合と看做すことができ、観測信号ｘ（ｔ）は、原信号ｓ（ｔ）に対して伝達関数Ｈ（ｚ）が作用したものと言える。そして、信号歪み除去処理は、観測信号ｘ（ｔ）に対して逆フィルタＧ（ｚ）を作用させて復元信号ｙ（ｔ）を得る処理となる。このとき、信号歪み除去処理で得られた復元信号ｙ（ｔ）をフレーム分割して、それぞれに対して、それぞれの信号に基づいて求めた予測誤差フィルタ１−Ａｉ（ｚ）を適用して得られるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）はイノベーション系列ｅｉ（１），…，ｅｉ（Ｗ）に一致することが望ましい。もし、予測誤差フィルタ１−Ａｉ（ｚ）の出力信号ｄｉ（ｎ）がｄｉ（ｎ）＝ｅｉ（ｎ）〔１≦ｉ≦Ｆ，１≦ｎ≦Ｗ〕を満たすならば、［条件３］の条件下でΣｍ＝１ＭＨｍ（ｚ）Ｇｍ（ｚ）＝１となることが示せる（数学的証明については、参考文献Ａを参照されたい。）。つまり、ｓ（ｔ）＝ｙ（ｔ）が言える。このとき、１−Ａｉ（ｚ）は１−Ｂｉ（ｚ）に等しくなる。
（参考文献Ａ）  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情報量）を用いると、解くべき問題は式（７）のように定式化される。ただし、ｇ＝［ｇ１Ｔ，…，ｇＭＴ］Ｔ，ｇｍ＝［ｇｍ（０），…，ｇｍ（Ｌ）］Ｔ，ａ＝［ａ１Ｔ，…，ａＦＴ］Ｔ，ａｉ＝［ａｉ（１），…，ａｉ（Ｐ）］Ｔとし、ａｉ（ｋ）は予測誤差フィルタ係数を表す。Ｉ（Ｕ１，…，Ｕｎ）は確率変数Ｕｉ間の相互情報量を表す。またｇおよびａに記号＾を付したものは、得るべき最適解である。Ｔは転置を表す。
【数６】
2007094463.tif 000008
拘束条件
［１］  ‖ｇ‖＝１  （但し‖・‖はノルムを表す。）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。【００３２】
相互情報量Ｉは、イノベーション推定値系列ｄ１（１），…，ｄ１（Ｗ），…，ｄｉ（１），…，ｄｉ（Ｗ），…，ｄＦ（１），…，ｄＦ（Ｗ）の振幅が定数倍されても変化しない。式（７）の拘束条件[１]は、この振幅の不定性を排除するための条件である。式（７）の拘束条件[２]は、上記[条件１]に対応して、予測誤差フィルタを最小位相系に制限するための条件である。以下、Ｉを、イノベーション推定値系列を入力としそれらの間の相互情報量を出力する関数と看做して、損失関数と呼称することにする。【００３３】
１．６  損失関数の導出
式（７）の最適化を実行するためには、損失関数Ｉ（ｄ１（１），…，ｄＦ（Ｗ））を有限長の信号系列｛ｄｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝から推定しなければならない。（多変量）確率変数Ｕの微分エントロピーをＤ（Ｕ）と表記すると、Ｉ（ｄ１（１），…，ｄＦ（Ｗ））は式（８）で定義される。ただし、ｄ＝［ｄＦＴ，…，ｄ１Ｔ］Ｔ、ｄｉ＝［ｄｉ（Ｗ），…，ｄｉ（１）］Ｔである。
【数７】
2007094463.tif 000009
【００３４】
ｙ＝［ｙＦＴ，…，ｙ１Ｔ］Ｔ，ｙｉ＝［ｙｉ（Ｗ），…，ｙｉ（１）］Ｔとおくと、ｄはｙを用いて、ｄ＝Ａｙと表される。ただし、行列Ａは、式（９）および式（１０）で表される。
【数８】
2007094463.tif 000010
【００３５】
よって、Ｄ（ｄ）は式（１１）のように表される。
【数９】
2007094463.tif 000011
【００３６】
多変量確率変数Ｕの共分散行列をΣ（Ｕ）と表記すると、式（１１）右辺第二項について、Σ（ｄ）＝Ｅ｛ｄｄＴ｝＝ＡＥ｛ｙｙＴ｝ＡＴ＝ＡΣ（ｙ）ＡＴが成立するから、式（１２）が成り立つ。
【数１０】
2007094463.tif 000012
【００３７】
式（１１）、式（１２）を式（８）に代入すると、式（１３）を得る。ただし、σ（Ｕ）２は確率変数Ｕの分散を表す。
【数１１】
2007094463.tif 000013
【００３８】
式（１３）でＪ（Ｕ）は（多変量）確率変数Ｕのネゲントロピー（negentropy）である。ネゲントロピーはＵの非ガウス性の度合いを表す非負の値をとり、Ｕがガウス分布に従う場合に限り０をとる。Ｃ（Ｕ１，…，Ｕｎ）は式（１４）で定義される。Ｃ（Ｕ１，…，Ｕｎ）は確率変数Ｕｉ間の相関の度合いを表す非負の値をとり、これらが無相関の場合に限り０をとる。
【数１２】
2007094463.tif 000014
【００３９】
ところで、ｓ＝［ｓＦＴ，…，ｓ１Ｔ］Ｔ，ｓｉ＝［ｓｉ（Ｗ），…，ｓｉ（１）］ＴとおくとＪ（ｙ）＝Ｊ（ｓ）＝ｃｏｎｓｔａｎｔとなるため（証明略）、式（１３）は更に式（１５）のように簡単化できる。
【数１３】
2007094463.tif 000015
【００４０】
以上から、式（１６）の最適化問題を解けばよいことになる。
【数１４】
2007094463.tif 000016
拘束条件
［１］  ‖ｇ‖＝１  （但し‖・‖はノルムを表す。）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。【００４１】
１．７  交代変数法による最適化
式（１６）について、交代変数の方法により、ｇとａを最適化する。すなわち、ｒ回目の繰り返しにおけるｇ及びａの推定値をそれぞれｇ＾（ｒ）、ａ＾（ｒ）と表せば、式（１７）および式（１８）の交互の最適化により更新された推定値ｇ＾（ｒ＋１），ａ＾（ｒ＋１）を得る。なお、ｇ＾およびａ＾は、記号＾がｇ、ａのそれぞれの上に付されたものを表す。例えば繰り返し回数の上限をＲ１とすれば、Ｒ１回目で得られるｇ＾（Ｒ１＋１）、ａ＾（Ｒ１＋１）が式（１６）の最適解である。上付き文字のＲ１は、Ｒ１である。
【数１５】
2007094463.tif 000017
拘束条件
［１］  ｇ＝ｇ＾（ｒ）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。
【数１６】
2007094463.tif 000018
拘束条件
［１］  ａ＝ａ＾（ｒ＋１）
［２］  ‖ｇ‖＝１【００４２】
式（１７）の意図するところは、伝達特性を打ち消すための逆フィルタの現在の推定値に基づいて原信号に固有の特性を打ち消すための予測誤差フィルタを推定することである。同様に、式（１８）の意図するところは、予測誤差フィルタの現在の推定値に基づいて逆フィルタを推定することである。イノベーション推定値系列ｄ１（１），…，ｄ１（Ｗ），…，ｄｉ（１），…，ｄｉ（Ｗ），…，ｄＦ（１），…，ｄＦ（Ｗ）が互いにより独立になるようにこれら２種類の最適化を繰り返すことで、逆フィルタと予測誤差フィルタを並行して推定することが可能になっている。したがって、ここでの繰り返しは逆フィルタの高精度な推定のために重要である。但し、図８から明らかなように処理する観測信号長が長くなる程、繰り返し回数は１回でも或る程度の効果が得られることが見て取れる。従って、この発明では、繰り返し回数は１回でもよい。【００４３】
１．８  ａの最適化
本発明では、式（１７）の最適化を以下のように行う。
まず注意すべきことは、Ｃ（ｄ１（１），…，ｄＦ（Ｗ））はｄｉ（ｎ）の二次の統計量に関連するのに対して、Ｊ（ｄｉ（ｎ））はｄｉ（ｎ）の高次の統計量に関連する値である。二次の統計量は信号の振幅情報のみ提供するが、高次の統計量は位相情報も提供する。したがって、一般に、高次統計量を含む最適化は、非最小位相系を導く可能性がある。そこで、１−Ａｉ（ｚ）が最小位相系であるという拘束条件から、ａの最適化においては式（１９）の最適化問題を解く。
【数１７】
2007094463.tif 000019
拘束条件
［１］  ｇ＝ｇ＾（ｒ）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。【００４４】
Ｃ（ｄ１（１），…，ｄＦ（Ｗ））は式（２０）で与えられる。
【数１８】
2007094463.tif 000020
【００４５】
ここで、行列Ａは式（９）および式（１０）に示すように上三角行列でその対角成分がすべて１であるから、ｌｏｇ ｄｅｔ Ａ＝０である。これを式（１２）に代入することで式（２１）の関係を得る。
【数１９】
2007094463.tif 000021
【００４６】
よって、式（１９）は、式（２２）の最適化問題と等価である。なお、式（２２）は、上記[条件２]を反映した表現になっていることに留意しなければならない。式（２２）を説明すれば、式（２２）は、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるａを求めよ、と云っている。
【数２０】
2007094463.tif 000022
拘束条件
［１］  ｇ＝ｇ＾（ｒ）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。【００４７】
式（２２）で表される最適化問題を解くことは、観測信号にｇ＾（ｒ）で与えられる逆フィルタを適用して得られるアドホック信号に対して、各フレームにおいて線形予測分析を行うことと等価であり、必ず最小位相予測誤差フィルタを得ることができる。線形予測分析に関しては、上記参考文献１を参照されたい。【００４８】
なお、式（２２）では、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大きい。この場合、対数関数は単調増加関数であるから、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めることができる。【００４９】
１．９  ｇの最適化
本発明では、式（１８）の最適化を以下のように行う。
前述したとおり、Ｃ（ｄ１（１），…，ｄＦ（Ｗ））は｛ｄｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝の相関の度合いに関わる指標であるが、（ｒ＋１）回目のａの最適化において最小化されているため、Σｉ＝１ＦΣｎ＝１ＷＪ（ｄｉ（ｎ））に比べて無視できる。そこでｇの最適化においては、式（２３）の最適化問題を解く。
【数２１】
2007094463.tif 000023
拘束条件
［１］  ａ＝ａ＾（ｒ＋１）
［２］  ‖ｇ‖＝１【００５０】
Ｊ（ｄｉ（ｎ））は、［条件２］に基づいて、式（２４）によって近似できる。この詳細は参考文献２を参照されたい。ただし、確率変数Ｕについて、κ４（Ｕ）はＵの尖度（４次キュムラント）を表す。式（２４）の右辺を第ｉフレームにおける正規化尖度という。
（参考文献２）  A.Hyvarinen, J.Karhunen, E.Oja, "INDEPENDENT COMPONENT ANALYSIS", John Wiley & Sons, Inc. 2001.
【数２２】
2007094463.tif 000024
【００５１】
［条件２］から音声信号のイノベーションの尖度は正であるため、κ４（ｄｉ（ｎ））／σ（ｄｉ（ｎ））４は正である。従って、式（２３）の最適化問題は、式（２５）の最適化問題に帰着する。σ（ｄｉ（ｎ）），κ４（ｄｉ（ｎ））は、［条件１］で述べた音声信号の局所的な定常性に基づいて、各フレーム内のサンプルから計算される。式（２６）では、１／Ｗを付加しているが、これは後の計算の便宜に過ぎず、式（２５）でｇの最適解を求めるにあたり影響を及ぼすものではない。式（２５）および式（２６）から、正規化尖度の全フレームでの総和が最大となるときのｇが、ｇ＾（ｒ＋１）となる。なお、式（２５）および式（２６）は、上記[条件２]を反映した表現になっていることに留意しなければならない。式（２５）および式（２６）を説明すれば、これらは、第ｉフレームにおける正規化尖度を全フレームで加算した値が最大となるｇを求めよ、と云っている。
【数２３】
2007094463.tif 000025
拘束条件
［１］  ａ＝ａ＾（ｒ＋１）
［２］  ‖ｇ‖＝１  【００５２】
式（２５）に従ってｇの最適解を求めるには、Ｑをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（２７）で表される更新則に従って求められる。ｇ′をｇ′のノルムで除しているのは上記拘束条件[２]を課すためである。η（ｕ）は学習率を表す。ｕは、ｇの最適化における更新回数を表す。
【数２４】
2007094463.tif 000026
【００５３】
式（２７）において、∇Ｑｇは式（２８）および式（２９）で与えられる。
【数２５】
2007094463.tif 000027
【００５４】
式（２９）において、ｄｉ（ｎ）は式（３０）で、ｖｍｉ（ｎ）は式（３１）および式（３２）で与えられる。ｘｍｉ（ｎ）は、ｍ番目のマイクロホンで観測された観測信号のｉ番目のフレームの信号である。
【数２６】
2007094463.tif 000028
【００５５】
§２  二次統計量に基づく信号歪み除去
上述の従来的手法の信号歪み除去方法は、比較的長時間の観測信号（例えば２０秒程度である。）を要する。これは、一般に、正規化尖度のような高次統計量を計算するためには大量の観測信号のサンプルが必要となるからである。しかし、実際にはそうした長時間の観測信号を利用できない場合が多い。このため、従来的手法の信号歪み除去方法の適用分野は極めて限られていた。
また高次統計量の計算は比較的複雑であるため、従来的手法の信号歪み除去方法では装置の構成が複雑になりやすい。
そこで、観測信号がより短時間（例えば３秒から５秒程度である。）の場合にも有効であり、かつ計算が従来に比して容易な信号歪み除去の原理を説明する。この原理は、信号の二次統計量のみを用いるものであり、§１で説明した本発明の基本原理より派生する。【００５６】
２．１  二次統計量に基づく信号歪み除去の原理
二次統計量に基づく信号歪み除去では、上述の３つの条件に、次の２つの条件を設定する。
[条件４]  Ｍ≧２である。すなわち、複数本のマイクロホンを用いる。
[条件５]  Ｈｍ＝｛ｈｍ（ｋ）｝ｋ＝０Ｋは相異なるマイクロホンｍの間で共通の零点を持たない。【００５７】
上記の式（１６）の最適化問題では、高次の統計量に関する値であるネゲントロピーＪおよび確率変数間の相関の度合いを示す指標Ｃを含む値を最小化するｇおよびａを求めた。
確率変数間の相関の度合いを示す指標Ｃは、二次の統計量で規定される。そこで、解くべき最適化問題を式（３３）で定式化する。
【数２７】
2007094463.tif 000029
【００５８】
式（２１）を参酌すれば、式（３３）の最適化問題は、式（３４）の最適化問題に転化される。なお、式（３４）は、上記[条件２]を反映した表現になっていることに留意しなければならない。式（３４）を説明すれば、式（３４）は、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるｇおよびａを求めよ、と云っている。
【数２８】
2007094463.tif 000030
【００５９】
ところで、上記の[条件４]および[条件５]が成立する場合、多チャンネルの観測信号は、音源からの原信号によって駆動されるＡＲ系として捉えることができる（参考文献３参照）。このことは、逆フィルタＧの先頭タップを式（３５）のように固定できることを意味する。但し、ｍ＝１に相当するマイクロホンは、最も音源に近いマイクロホンである。
（参考文献３）  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.
【数２９】
2007094463.tif 000031
【００６０】
式（３４）および式（３５）で規定されるｇを係数とする逆フィルタＧを、式（６）に従って観測信号ｘ（ｔ）に適用することで伝達特性が除去された復元信号ｙ（ｔ）を得る。【００６１】
２．２  ａの最適化
式（３４）について、交代変数の方法により、ｇとａを最適化する。
逆フィルタの係数ｇｍ（ｋ）を固定した状態で予測誤差フィルタの係数ａｉ（ｋ）に関して式（３４）の損失関数を最小化する。
このとき、次の２点に注意する。１点目は、ｇ＝［ｇ１Ｔ，…，ｇＭＴ］Ｔは固定されているので、逆フィルタＧの出力である復元信号ｙ（ｔ）は予測誤差フィルタの最適化において不変であるということである。２点目は、ｉ番目のフレームにおける予測誤差フィルタの係数ａｉ（１），…，ａｉ（Ｐ）は、ｄｉ（１），…，ｄｉ（Ｗ）にのみ寄与するということである。【００６２】
従って、各々のフレームにおいて、Σｎ＝１Ｗlogσ（ｄｉ（ｎ））２を最小化するように、予測誤差フィルタの係数ａｉ（１），…，ａｉ（Ｐ）を推定すればよい。[条件２]からイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散はフレーム内で定常であるから、Σｎ＝１Ｗlogσ（ｄｉ（ｎ））２の最小化は、Ｗ＊σ（ｄｉ（ｎ））２の最小化と等価である。記号＊は乗算を表す。分散σ（ｄｉ（ｎ））２は、＜ｄｉ（ｎ）２＞ｎ＝１Ｗとして計算できる。但し、＜ｄｉ（ｎ）２＞ｎ＝１Ｗは、１フレーム分のイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を使って計算したｄｉ（ｎ）の２乗平均を表す。結局、係数ａｉ（ｋ）は、Ｗ＊＜ｄｉ（ｎ）２＞ｎ＝１Ｗ、つまりｄｉ（ｎ）の２乗総和が最小になるときのａｉ（ｋ）として推定される。このような係数ａｉ（ｋ）は、線形予測分析の手法によって計算される。【００６３】
なお、ここでは、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めるとして説明したが、これに限定する趣旨ではない。既述のとおり、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めることができる。【００６４】
２．３  ｇの最適化
予測誤差フィルタの係数ａｉ（ｋ）を固定した状態で逆フィルタの係数ｇｍ（ｋ）に関して式（３４）の損失関数を最小化する。
逆フィルタの係数ｇｍ（ｋ）に関する損失関数の最小化には、勾配法を用いる。[条件２]を用いると、式（３４）の最適化問題は、式（３６）の最適化問題に転化される。
【数３０】
2007094463.tif 000032
【００６５】
式（３６）に従ってｇの最適解を求めるには、Σｉ＝１Ｆlog＜ｄｉ（ｎ）２＞ｎ＝１Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（３７）で表される更新則に従って求められる。δは学習率を表す。１≦ｍ≦Ｍ、１≦ｋ≦Ｌである。なお、式（３７）では、式（３５）の条件から‖ｇ‖＝１なる拘束条件が課されないことに留意しなければならない。同様に式（３５）の条件からｋが取る値の範囲は１≦ｋ≦Ｌである。
【数３１】
2007094463.tif 000033
【００６６】
上記式（２９）あるいは上記非特許文献１に記載された式（３）と比べて明らかなように、式（３７）の右辺第二項は、二次の統計量で表されており、この計算に高次の統計量を必要としない。このため、高次統計量を計算するには短い時間の観測信号の場合にも有効であり、計算自体も容易である。【００６７】
なお、式（３６）では、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのｇをｇ＾として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大である。この場合、対数関数は単調増加関数であるから、第ｉフレームにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の全フレームでの総和が最小となるときのｇをｇ＾として求めることができる。なお、この場合には式（３７）で示した更新則は適用できなくなり、改めてΣｉ＝１Ｆ＜ｄｉ（ｎ）２＞ｎ＝１Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この結果得られる更新則はＩＣＡと同様の考え方で定式化できるから略する。【００６８】
§３  プリ・ホワイトニング
本発明の信号歪み除去には、プリ・ホワイトニングを適用することができる。観測信号をプリ・ホワイトニングすることで、最適化計算の安定化、とくに更新則の高速な収束が可能となる。
各マイクロホンで得られた観測信号系列全体｛ｘｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。
式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号ｘｍ（ｔ）に適用する。ｗｍ（ｔ）は、ｍ番目のマイクロホンで得た観測信号ｘｍ（ｔ）を白色化した信号を表す。
【数３２】
2007094463.tif 000034
【００６９】
このとき、式（３１）および式（３８）は式（４０）に、式（３２）は式（４１）に、変更すればよい。
【数３３】
2007094463.tif 000035
【００７０】
§４  実施形態
以下、本発明の実施形態を図面を参照して説明する。本発明の実施形態として、後述の各実施形態に限定するものではなく、各セクションで説明した原理を実現する実施形態であればよい。【００７１】
＜第１実施形態＞
本発明の第１実施形態を実施する場合、以下の手順に従ってセンサで得た観測信号を処理する。ここでは、実施形態を具体的に説明する観点から信号として音声信号を例に挙げて説明する。
なお、第１実施形態の説明に先立ち、観測信号およびフレーム化処理について概説する。【００７２】
((観測信号))
図示しないセンサ（例えばマイクロホン）によって得られたアナログ信号（このアナログ信号には伝達特性に由来する歪みが重畳されている。）は、例えば８,０００Ｈｚのサンプリングレートでサンプリングされ、適宜量子化された離散信号に変換される。以下、この離散信号を観測信号ということにする。アナログ信号から観測信号へのＡ／Ｄ変換などを実行するために必要となる構成要素（手段）は、いずれも公知技術の常套手段によって達成されるから、説明および図示を略する。【００７３】
((フレーム化処理))
図示しない信号切出手段が、離散信号から、時間軸方向に一定時間幅で始点を移動させながら、所定時間長の離散信号をフレームとして切り出す。例えば２００サンプル点（８,０００Ｈｚ×２５ｍｓ）長の離散信号を、８０サンプル点（８,０００Ｈｚ×１０ｍｓ）ずつ始点を移動させながら切り出す。各フレームの切り出しは、離散信号に公知の窓関数（例えば、ハミング窓、ガウス窓、方形窓など）を適用すればよい。窓関数の適用によるフレームの切り出しは公知の常套手段によって達成される。【００７４】
本発明の第１実施形態である信号歪み除去装置（１）をコンピュータ（汎用機）で実現する場合のハードウェア構成例を説明する。
図２に例示するように、信号歪み除去装置（１）は、キーボード、ポインティングデバイスなどが接続可能な入力部（１１）、液晶ディスプレイ、ＣＲＴ（Cathode Ray Tube）ディスプレイなどが接続可能な出力部（１２）、信号歪み除去装置（１）外部に通信可能な通信装置（例えば通信ケーブル、ＬＡＮカード、ルータ、モデムなど）が接続可能な通信部（１３）、ＤＳＰ（Digital Signal Processor）（１４）〔ＣＰＵ（Central Processing Unit）でも良い。またキャッシュメモリやレジスタ（１９）などを備えていてもよい。〕、メモリであるＲＡＭ（１５）、ＲＯＭ（１６）やハードディスク、光ディスク、半導体メモリなどである外部記憶装置（１７）並びにこれらの入力部（１１）、出力部（１２）、通信部（１３）、ＤＳＰ（１４）、ＲＡＭ（１５）、ＲＯＭ（１６）、外部記憶装置（１７）間のデータのやり取りが可能なように接続するバス（１８）を有している。また必要に応じて、信号歪み除去装置（１）に、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＤＶＤ（Digital Versatile Disc）などの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。【００７５】
信号歪み除去装置（１）の外部記憶装置（１７）には、信号歪み除去のためのプログラムおよびこのプログラムの処理において必要となるデータ（観測信号）などが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶され、他のプログラムの処理に供されるときに、ＲＡＭや外部記憶装置などから読み込まれる。【００７６】
より具体的には、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭなど〕には、観測信号に逆フィルタを適用する処理のためのプログラム、観測信号に逆フィルタを適用して得られた信号から予測誤差フィルタを求める処理のためのプログラム、予測誤差フィルタから逆フィルタを求める処理のためのプログラム、およびこれらのプログラムの処理において必要となるデータ（フレーム単位の観測信号など）が記憶されている。その他、これらのプログラムに基づく処理を制御するための制御プログラムも適宜に保存しておく。【００７７】
第１実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、制御部）を実現することで、信号歪み除去が実現される。【００７８】
そこで次に、図３〜図５を参照して、信号歪み除去装置（１）における信号歪み除去処理の流れを順次説明する。
大まかな処理の手順は、（ａ）観測信号ｘ（ｔ）に対して逆フィルタを適用した信号（以下、アドホック信号という。）を求め、（ｂ）アドホック信号から予測誤差フィルタを求め、（ｃ）この予測誤差フィルタから逆フィルタを求め、（ｄ）前記（ａ）、（ｂ）、（ｃ）の処理を繰り返して最適な逆フィルタを求め、（ｅ）最適化された逆フィルタを観測信号に対して適用した信号を復元信号ｙ（ｔ）として得る。
（ｂ）は上述のａの最適化に相当し、（ｃ）は上述のｇの最適化に相当し、（ｄ）は、式（１７）および式（１８）に相当する。（ｄ）の処理の繰り返し回数は予め定めた回数Ｒ１とする。つまり、１≦ｒ≦Ｒ１とする。また、（ｃ）の処理でｇを最適化する更新則の更新回数は予め定めた回数Ｒ２とする。つまり、１≦ｕ≦Ｒ２とする。（ｄ）の処理、つまり（ａ）、（ｂ）、（ｃ）の一連の処理を１回行うたびに、更新則によるＲ２回の更新が行なわれる。実施形態では、回数Ｒ１は、予め定めた回数とするが、これに限定されず、例えば、ｒ回目のｇを算出したときの式（２６）のＱの値とｒ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。同様に、回数Ｒ２は、予め定めた回数とするが、これに限定されず、例えば、ｕ回目のｇを算出したときの式（２６）のＱの値とｕ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。【００７９】
（ステップＳ１００）
逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタを観測信号ｘ（ｔ）＝［ｘ１（ｔ），…，ｘｍ（ｔ），…，ｘＭ（ｔ）］Ｔに適用することで、アドホック信号ｙ（ｔ）を求める。アドホック信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、ここでは後述のＲ１回の処理を経て求められた復元信号ではないことを明示するためアドホック信号と呼称する。ここでｔは、全てのサンプル番号を示し、１≦ｔ≦Ｎとする。Ｎは全サンプル数である。第１実施形態では、マイクロホンの個数Ｍは１以上であればよい。
【数３４】
2007094463.tif 000036
【００８０】
逆フィルタの係数列｛ｇｍ（ｋ）；０≦ｋ≦Ｌ｝として、繰り返し回数Ｒ１の初回には予め定められた初期値を、２回目以降には後述する逆フィルタ計算部（１３）によって求められた逆フィルタｇ＾（ｒ＋１）を用いる。【００８１】
（ステップＳ１０１）
予測誤差フィルタ計算部（１５）は、フレーム化処理を行うフレーム化処理部（１５１）とフレーム予測誤差フィルタ計算部（１５２）によって構成される。そして、フレーム予測誤差フィルタ計算部（１５２）は、第ｉフレームのアドホック信号から予測誤差フィルタを求める第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）からなる。ただし、ｉは、１≦ｉ≦Ｆを満たす整数である。【００８２】
フレーム化処理部（１５１）は、逆フィルタ適用部（１４）で求められたアドホック信号｛ｙ（ｔ）；１≦ｔ≦Ｎ｝をフレーム化処理する。フレーム化処理は、例えば式（４３）のように、Ｗ点分を切り出す窓関数をＷ点ずつシフトさせて適用することにより行う。｛ｙｉ（ｎ）；１≦ｎ≦Ｗ｝はｉ番目のフレームに含まれるアドホック信号列を表す。
【数３５】
2007094463.tif 000037
【００８３】
そして、第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）は、式（２２）に従って、第ｉフレームのアドホック信号列｛ｙｉ（ｎ）；１≦ｎ≦Ｗ｝に対してＰ次の線形予測分析を行い、予測誤差フィルタの係数列｛ａｉ（ｋ）；１≦ｋ≦Ｐ｝を計算する。この算出方法は、上記参考文献１を参照されたい。ここで得られたａ１（１），…，ａ１（Ｐ），…，ａｉ（１），…，ａｉ（Ｐ），…，ａＦ（１），…，ａＦ（Ｐ）は、式（２２）のａ＾（ｒ＋１）を与える。【００８４】
（ステップＳ１０２）
逆フィルタ計算部（１３）の機能構成例を、図４を参照して説明する。逆フィルタ計算部（１３）は勾配計算部（１３１）、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。更に、勾配計算部（１３１）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）とを備えて構成される。ここで更新用逆フィルタは、式（２７）のｇ〈ｕ〉に相当する。【００８５】
第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖｍｉ（ｎ）を計算する（式（３１）を参照）。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。【００８６】
第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｙｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用して予測誤差フィルタ適用後のイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。爾後、第２の予測誤差フィルタ適用部（１３１２）は、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）を入力とする。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。【００８７】
勾配ベクトル計算部（１３１３）は、信号ｖｍｉ（ｎ）とイノベーション推定値ｄｉ（ｎ）とを用いて現在の更新用逆フィルタｇ〈ｕ〉の勾配ベクトル∇Ｑｇを計算する（式（２８）および式（２９）を参照）。有限個のサンプルｖｍｉ（ｎ）およびｄｉ（ｎ）を用いて式（２９）を演算するときは、期待値Ｅをサンプルから求めればよい。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。【００８８】
逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇ〈ｕ〉、学習率η（ｕ）、勾配ベクトル∇Ｑｇを入力として、式（２７）に従って、ｕ＋１回目の更新用逆フィルタｇ〈ｕ＋１〉を求める。式（２７）は、求められたｇ〈ｕ＋１〉を新たなｇ〈ｕ〉と見立てて更新を行なうことを意味する。【００８９】
更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇ〈ｕ＋１〉、つまり新たなｇ〈ｕ〉および観測信号ｘ（ｔ）を入力として、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。具体的には、式（４２）のｇｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。なお、更新用逆フィルタ適用後信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、ここでは後述のＲ１回の処理を経て求められた復元信号ではなく、更新則を行なうために算出される信号であることを明示するため更新用逆フィルタ適用後信号と呼称する。【００９０】
制御部（６００）の制御によってＲ２回の更新が行なわれた結果として得られたｇ〈Ｒ２＋１〉は、式（２５）のｇ＾（ｒ＋１）に相当する。上付き文字のＲ２は、Ｒ２である。逆フィルタ計算部（１３）は、ｇ＾（ｒ＋１）を出力する。【００９１】
制御部（５００）の制御によって、上述の一連の処理を１回行うごとにｒに１を加算してｒがＲ１に等しくなるまで、つまり上述の一連の処理をＲ１回繰り返すことで（ステップＳ１０３）、ｇ＾（Ｒ１＋１）を得る。上付き文字のＲ１は、Ｒ１である。このｇ＾（Ｒ１＋１）が、式（１６）の最適解とされる。そこで、ｇ＾（Ｒ１＋１）を得た段階で、逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタｇ＾（Ｒ１＋１）を観測信号ｘ（ｔ）＝［ｘ１（ｔ），…，ｘＭ（ｔ）］Ｔに適用することで、復元信号ｙ（ｔ）を得ることができる（ステップＳ１０４）。【００９２】
＜第２実施形態＞
第２実施形態は、第１実施形態の変形例に相当する。具体的には、§３で述べたプリ・ホワイトニングを行なう形態である。そこで、第１実施形態と異なる部分について図６および図７を参照して説明を加える。なお、プリ・ホワイトニングは観測信号に対して行なうプリ・プロセスであるから、ここで説明するプリ・ホワイトニングを行なう形態は、後述の第３実施形態にも適用可能である。
第２実施形態では、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭなど〕に、白色化フィルタを求める処理のためのプログラム、白色化フィルタを観測信号に適用する処理のためのプログラムも記憶されている。
第２実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、白色化フィルタ計算部、白色化フィルタ適用部）を実現することで、信号歪み除去が実現される。【００９３】
（ステップＳ１００ａ）
白色化フィルタ計算部（１１）は、各マイクロホンで得られた観測信号全体｛ｘｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。この計算は線形予測分析と同じであり上記参考文献１を参照されたい。白色化フィルタの係数は、白色化フィルタ適用部（１２）の入力となる。【００９４】
（ステップＳ１００ｂ）
白色化フィルタ適用部（１２）は、式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号に適用して、白色化信号ｗｍ（ｔ）を得る。既述のとおり、式（３１）は式（４０）に変更すればよいので、第１実施形態において、逆フィルタ計算部（１３）、とくに第１の予測誤差フィルタ適用部（１３１１）による処理を式（３１）ではなく式（４０）による計算処理に改めればよい。また、第１実施形態において、逆フィルタ適用部（１４）による処理を、式（４２）ではなく式（４４）による計算処理に改めればよい。ステップＳ１００ｂの処理の後、第１実施形態のステップＳ１００〜Ｓ１０４の処理を行うが、これらの処理では第１実施形態の各処理における観測信号をステップＳ１００ｂの処理で得られた白色化信号に読み替えて第１実施形態と同様の処理を行う。この意味で、図７では、第１実施形態のステップＳ１００〜Ｓ１０４の各処理に相当する処理を示す符号に記号′を付している。
【数３６】
2007094463.tif 000038
【００９５】
＜実施例１＞
発明者らは第２実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、マイクロホンの数Ｍ＝４、白色化フィルタの次数Ｘ＝５００、逆フィルタの次数Ｌ＝１０００、窓関数の切出しサンプル数（１フレームのサンプル数）Ｗ＝２００、予測誤差フィルタの次数Ｐ＝１６、繰り返し回数Ｒ１＝１０、逆フィルタ計算部の更新回数Ｒ２＝２０とした。学習率η（ｕ）は、初期値を０．０５に設定し、もし式（２７）によって式（２６）の値が減少するならば、η（ｕ）の値を順次半減させることで、式（２６）の値が必ず増大するようにした。図６に示した逆フィルタ適用部（１４）へ入力する初期逆フィルタは、式（４５）のように設定した。
【数３７】
2007094463.tif 000039
【００９６】
本発明の第２実施形態の効果を、信号歪み除去の指標としてＤ５０値（インパルス応答の全エネルギーに対する初期の５０ｍｓｅｃまでのエネルギーの比）を用いて評価した。連続発話データベースから男女各一名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
図８は、男声および女声について観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ１（各信号長Ｎを持つ観測信号を図６に示す逆フィルタ適用部（１４）と、予測誤差フィルタ計算部（１５）と、逆フィルタ計算部（１３）を一巡する処理を実行して逆フィルタを求める回数）とＤ５０値の関係を示している。いずれの場合においても、繰り返し回数を増加させるとＤ５０値が向上しており、繰り返し処理の効果が顕著に見て取れる。特に観測信号長が５〜１０秒程度の比較的短い長さであっても、繰り返し処理によってＤ５０値が大きく向上したことが分かる。【００９７】
また、本発明の第２実施形態の効果を、音声スペクトグラムの比較から検証した。
図９Ａは観測信号長が１分の場合の残響を含まない音声（原音声）、図９Ｂは観測信号長が１分の場合の残響を含む音声（観測音声）、図９Ｃは観測信号長が１分の場合の残響除去後の音声（復元音声）のスペクトログラムの抜粋を示している。図９Ａと図９Ｃとの対比および図９Ｂと図９Ｃとの対比から、観測信号に含まれる残響が抑制され、原音声固有の特徴である調波構造やフォルマント構造が回復されたことが分かる。【００９８】
また、本発明の第２実施形態の効果を、ＬＰＣスペクトル歪みを用いて検証した。
図１０Ｂは原音声の波形、図１０Ａは原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）および原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）を示している。図１０Ａおよび図１０Ｂの各横軸は秒単位の時間を表し、両図でタイムスケールを揃えている。図１０Ｂの縦軸は振幅値を表している。但し、原信号の相対的な大きさが分かればよいので、この趣旨から図１０Ｂの縦軸では単位を明示していない。図１０Ａの縦軸は、ＬＰＣスペクトル歪みＳＤ（ｄＢ）を表している。
図１０Ａから、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりも常に小さいことがわかる。なお、観測音声ではＬＰＣスペクトル歪みの平均が５．３９ｄＢ、分散が４．２０ｄＢであったのに対して、復元音声では平均が２．３８ｄＢ、分散が２．００ｄＢであった。
また、図１０Ａと図１０Ｂとの対比から、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）が大きな値を示す区間（例えば約１．０秒〜約１．２秒の区間を参照）は、原音声の波形の振幅値がほぼ０であることがわかる。実際、この区間では発声がなく無音区間である。このため、実際に知覚される歪みはかなり小さくなっていた。つまり、発声区間における原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりもかなり小さく、このため原音声のスペクトルを高い精度で復元できたことが結論付けられる。【００９９】
＜第３実施形態＞
第３実施形態は、第１実施形態の変形例に相当する。具体的には、§２で述べた二次統計量に基づく信号歪み除去処理を行なう形態である。そこで、第１実施形態と異なる部分について図１１および図１２を参照して説明を加える。但し、第３実施形態では、マイクロホンの個数Ｍは２以上とする。【０１００】
ステップＳ１００の処理およびステップＳ１０１の処理は、第１実施形態と同じである。【０１０１】
ステップＳ１０１の処理に続いて、ステップＳ１０２ａの処理を行う。
第３実施形態に係る逆フィルタ計算部（１３）の機能構成例を、図１１を参照して説明する。
逆フィルタ計算部（１３）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）と、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。ここで更新用逆フィルタは、式（３７）のｇｍ（ｋ）に相当する。【０１０２】
第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖｍｉ（ｎ）を計算する（式（３８）を参照）。具体的には、フレーム化処理部（４０２Ｂ）が、入力された観測信号ｘｍ（ｔ）に対してフレーム化処理を行い、観測信号ｘｍ（ｔ）のｉ番目のフレームの信号ｘｍｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０４ｉ）が信号ｘｍｉ（ｎ）を入力として、式（３８）に従って信号ｖｍｉ（ｎ）を出力する。但し、１≦ｉ≦Ｆである。【０１０３】
第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、各フレームに対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用して予測誤差フィルタ適用後のイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。具体的には、フレーム化処理部（４０２Ａ）が、初期値の場合を除き、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）に対してフレーム化処理を行い、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、ｉ番目のフレームの信号ｙｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０３ｉ）が信号ｙｉ（ｎ）を入力として、式（３０）に従ってイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を出力する。但し、１≦ｉ≦Ｆである。【０１０４】
勾配ベクトル計算部（１３１３）は、信号ｖｍｉ（ｎ）とイノベーション推定値ｄｉ（ｎ）とを用いて現在の更新用逆フィルタｇｍ（ｋ）の勾配ベクトルを計算する（式（３７）の右辺第二項を参照）。具体的には、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、相互相関計算部（４０５ｉ）は信号ｖｍｉ（ｎ）とイノベーション推定値ｄ（ｎ）との相互相関〈ｄｉ（ｎ）ｖｍｉ（ｎ−ｋ）〉ｎ＝１Ｗを計算する。また、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、分散計算部（４０６ｉ）は、イノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散〈ｄｉ（ｎ）２〉ｎ＝１Ｗを求める。各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、除算部（４０７ｉ）は、〈ｄｉ（ｎ）ｖｍｉ（ｎ−ｋ）〉ｎ＝１Ｗ／〈ｄｉ（ｎ）２〉ｎ＝１Ｗを求める。加算部（４０７）は、除算部（４０７１）〜（４０７Ｆ）の全フレームに亘る総和、つまり式（３７）の右辺第二項を求める。【０１０５】
逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇｍ（ｋ）、学習率δ、勾配ベクトルを入力として、式（３７）に従って、ｕ＋１回目の更新用逆フィルタｇｍ（ｋ）′を求める。式（３７）は、求められたｇｍ（ｋ）′を新たなｇｍ（ｋ）と見立てて更新を行なうことを意味する。【０１０６】
更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇｍ（ｋ）′、つまり新たなｇｍ（ｋ）および観測信号ｘ（ｔ）を入力として、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。具体的には、式（４２）のｇｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。【０１０７】
ステップＳ１０２ａの処理に続いて、ステップＳ１０３およびステップＳ１０４の処理を行うが、第１実施形態と同じであるから説明を略する。【０１０８】
＜実施例２＞
発明者らは第３実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、Ｍ＝４，Ｌ＝１０００，Ｗ＝２００，Ｐ＝１６，Ｒ１＝６，Ｒ２＝５０とした。学習率δは、初期値を０．０５に設定し、Σｉ＝１Ｆlog＜ｄｉ（ｎ）２＞ｎ＝１Ｗの値が増加するならば、学習率δの値を順次半減させることでΣｉ＝１Ｆlog＜ｄｉ（ｎ）２＞ｎ＝１Ｗの値が必ず減少するように設定した。逆フィルタの初期推定値は、ｇｍ（ｋ）＝０，１≦ｍ≦Ｍ，１≦ｋ≦Ｌとして設定した。【０１０９】
本発明の第３実施形態の効果を、音声明瞭度を表すＲＡＳＴＩ（参考文献５を参照）を残響除去の指標として評価した。連続発話データベースから男女各五名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
（参考文献５）  H. kuttruff. Room acoustics. Elsevier Applied Science, third edition, P.237 1991.【０１１０】
図１３は、観測信号長Ｎを３秒、４秒、５秒、１０秒としたときの、ＲＡＳＴＩの値を表示したものである。図１３に示すように、観測信号が３〜５秒のように短時間の場合でも、高い残響除去性能を示していることが分かる。
図１４は、残響除去前後におけるエネルギー減衰曲線の例である。直接音が到達してから５０ミリ秒後の反射音のエネルギーが１５ｄＢ低減されていることが分かる。
【産業上の利用可能性】
【０１１１】
本発明は、様々な信号処理システムの性能向上に寄与する要素技術であるところ、例えば音声認識システム、テレビ会議システム、補聴器、音楽情報処理システム等に利用することができる。
【図１】
2007094463.tif 000040
【図２】
2007094463.tif 000041
【図３】
2007094463.tif 000042
【図４】
2007094463.tif 000043
【図５】
2007094463.tif 000044
【図６】
2007094463.tif 000045
【図７】
2007094463.tif 000046
【図８】
2007094463.tif 000047
【図９】
2007094463.tif 000048
【図１０】
2007094463.tif 000049
【図１１】
2007094463.tif 000050
【図１２】
2007094463.tif 000051
【図１３】
2007094463.tif 000052
【図１４】
2007094463.tif 000053
【図１５】
2007094463.tif 000054
【手続補正書】
【提出日】平成19年11月14日(2007.11.14)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】請求項２
【補正方法】変更
【補正の内容】
【請求項２】
上記予測誤差フィルタ計算手段は、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算手段は、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタとして、上記各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタを求め、この逆フィルタを出力するものである
ことを特徴とする請求項１に記載の信号歪み除去装置。
【手続補正２】
【補正対象書類名】特許請求の範囲
【補正対象項目名】請求項３
【補正方法】変更
【補正の内容】
【請求項３】
上記予測誤差フィルタ計算手段は、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算手段は、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタとして、上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタを求め、この逆フィルタを出力するものである
ことを特徴とする請求項１に記載の信号歪み除去装置。
【手続補正３】
【補正対象書類名】特許請求の範囲
【補正対象項目名】請求項８
【補正方法】変更
【補正の内容】
【請求項８】
上記予測誤差フィルタ計算ステップは、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算ステップは、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタとして、上記各イノベーション推定値の正規化尖度の全フレームでの総和が最大となるときの逆フィルタを求め、この逆フィルタを出力するものである
ことを特徴とする請求項７に記載の信号歪み除去方法。
【手続補正４】
【補正対象書類名】特許請求の範囲
【補正対象項目名】請求項９
【補正方法】変更
【補正の内容】
【請求項９】
上記予測誤差フィルタ計算ステップは、
上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの予測誤差フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの予測誤差フィルタ、を求めるとして、各フレームのアドホック信号の線形予測分析を行ない、各フレームごとの上記予測誤差フィルタを出力するものであり、
上記逆フィルタ計算ステップは、
上記イノベーション推定値系列がその全サンプル間で独立となる上記逆フィルタとして、上記各イノベーション推定値の分散の全フレームでの総和が最小となるときの逆フィルタ、または、上記各イノベーション推定値の分散の対数値の全フレームでの総和が最小となるときの逆フィルタを求め、この逆フィルタを出力するものである
ことを特徴とする請求項７に記載の信号歪み除去方法。
【手続補正５】
【補正対象書類名】明細書
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、信号の歪み（ひずみ：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.
【００１６】
１．２  音声信号のモデル化
まず、原信号である音声信号ｓ（ｔ）を、以下の３つの条件を満足する信号としてモデル化する。
【００１７】
［条件１］音声信号ｓ（ｔ）は、区分定常なＡＲ過程で生成される。
この[条件１]から、ＡＲ過程の次数をＰ、定常とみなせる区間長をＷサンプルとして音声信号ｓ（ｔ）をフレーム化すると、第ｉフレームの音声信号ｓｉ（ｎ）は、式（１）のように表される。式（２）は、第ｉフレームの音声信号ｓｉ（ｎ）のサンプルと、フレーム化前の音声信号ｓ（ｔ）のサンプルとの対応を示している。つまり、第ｉフレームのｎ番目のサンプルは、フレーム化前の音声信号ｓ（ｔ）において、（ｉ−１）Ｗ＋ｎ番目のサンプルに相当する。式（１）および式（２）において、ｂｉ（ｋ）は線形予測係数、ｅｉ（ｎ）はイノベーションを表す。但し、１≦ｎ≦Ｗ、１≦ｔ≦Ｎ、Ｎは全サンプル数である。以下、特に断りの無い限り、パラメータｎは１フレームのサンプル番号を表し、パラメータｔは全てのサンプル番号を表す。また、全フレーム数はＦとする。
【数１】
2007094463.tif 000055
【００１８】
なお、第ｉフレームにおけるｎ番目のイノベーションｅｉ（ｎ）についても、第ｉフレームのｎ番目のイノベーションｅｉ（ｎ）と、フレーム化前の音声信号ｓ（ｔ）に対するイノベーションｅ（ｔ）との対応を示すことができる。この場合、第ｉフレームのｎ番目のイノベーションｅｉ（ｎ）は、フレーム化前のイノベーションｅ（ｔ）において（ｉ−１）Ｗ＋ｎ番目のイノベーションに相当し、ｅｉ（ｎ）＝ｅ（（ｉ−１）Ｗ＋ｎ）が成り立つ。
【００１９】
式（１）をｚ変換する。左辺のｚ変換をＳｉ（Ｚ）とし、右辺第二項のｚ変換をＥｉ（Ｚ）とし、Ｂｉ（ｚ）＝Σｋ＝１Ｐｂｉ（ｋ）ｚ−ｋとすれば、右辺第一項は、Ｂｉ（ｚ）Ｓｉ（Ｚ）となる。従って、式（１）のｚ変換は、（１−Ｂｉ（ｚ））Ｓｉ（Ｚ）＝Ｅｉ（Ｚ）である。なお、ｚ−１は時間領域では１タップ遅延素子に相当する。以降、時間領域信号（タップ重み係数）を小文字で、ｚ領域信号（伝達関数）を大文字でそれぞれ表す。１−Ｂｉ（ｚ）は最小位相性を満足しなければならず、『１−Ｂｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ』ことが要求される。
【００２０】
［条件２］第ｉフレームに属するイノベーションｅｉ（１），…，ｅｉ（Ｗ）は独立且つ
同一分布に属する。イノベーションｅｉ（１），…，ｅｉ（Ｗ）の確率分布の平均及び歪度（３次キュムラント）は０、尖度（４次キュムラント）は正である。さらに、異なるフレームｉ、ｊ〔ｉ≠ｊ〕に属するイノベーションｅｉ（ｎ）とｅｊ（ｎ′）同士も独立である。ただし、これらは必ずしも同一分布に属するとは限らない。
［条件３］予測誤差フィルタ１−Ｂｉ（ｚ）は、相異なるフレーム間で共通する零点をもたない。
【００２１】
式（１）および式（２）から、音声信号ｓ（ｔ）は、式（３）のように表される。[・]は床関数を表す。
【数２】
2007094463.tif 000056
【００２２】
このとき、[条件２]は、『イノベーション過程ｅ（ｔ）は時間的に独立な信号である。また、その統計的性質（あるいは統計量）はフレーム内では定常である。』と表現できる。また、[条件３]は、『線形予測係数｛ｂｉ（ｋ）｝ｋ＝１Ｐで表されるＡＲ系は、時不変な極を持たない』と表現できる。
【００２３】
１．３  観測信号のモデル化
次に、Ｍ個のマイクロホンで音声信号を観測して観測信号を得たときの観測信号をモデル化する。但しＭは、Ｍ≧１の整数である。
ｍ番目（１≦ｍ≦Ｍ）のマイクロホンで観測される残響信号ｘｍ（ｔ）を、音源からｍ番目のマイクロホンに至る経路の伝達関数Ｈｍ（ｚ）のタップ重み係数｛ｈｍ（ｋ）；０≦ｋ≦Ｋ；Ｋはインパルス応答の持続時間とする。｝を用いて式（４）のようにモデル化する。ここでは、音声信号の場合の伝達特性の代表例として残響を挙げて、伝達特性を残響に言い換えて説明する。但し、伝達特性を残響に限定する趣旨ではない。
【数３】
2007094463.tif 000057
【００２４】
Ｍ個の観測信号についてまとめて表現すれば、式（５）のように表すことができる。但し、式（５）において、ｘ（ｔ）＝［ｘ１（ｔ），…，ｘＭ（ｔ）］Ｔ、ｈ（ｋ）＝［ｈ１（ｋ），…，ｈＭ（ｋ）］Ｔである。
【数４】
2007094463.tif 000058
【００２５】
１．４  信号歪み除去の原理
信号歪み除去後の復元信号ｙ（ｔ）は、多チャネル逆フィルタ｛Ｇｍ（ｚ）；１≦ｍ≦Ｍ｝のタップ重み係数｛ｇｍ（ｋ）；１≦ｍ≦Ｍ，０≦ｋ≦Ｌ；Ｌは逆フィルタの次数｝を用いて式（６）により計算される。本発明においては、逆フィルタ係数であるｇｍ（ｋ）を観測信号ｘ１（ｔ），…，ｘＭ（ｔ）のみから推定する。
【数５】
2007094463.tif 000059
【００２６】
１．５  本発明の基本原理
本発明の基本原理は、伝達関数｛Ｈｍ（ｚ）；１≦ｍ≦Ｍ｝の逆フィルタ｛Ｇｍ（ｚ）；１≦ｍ≦Ｍ｝とＡＲフィルタ｛１／（１−Ｂｉ（ｚ））；１≦ｉ≦Ｆ｝の逆フィルタである予測誤差フィルタ｛１−Ａｉ（ｚ）；１≦ｉ≦Ｆ｝とを並行して推定することを主な特徴とする。
【００２７】
この基本原理を説明するため、上述のモデル機構を組み込んだ系全体の構成線図を図１に示す。上述のモデル化に拠れば、原信号ｓ（ｔ）は、フレームごとのイノベーション系列ｅｉ（１），…，ｅｉ（Ｗ）に対してＡＲフィルタ１／（１−Ｂｉ（ｚ））を適用して得られる信号ｓ１（ｎ）,…，ｓＦ（ｎ）のフレーム結合と看做すことができ、観測信号ｘ（ｔ）は、原信号ｓ（ｔ）に対して伝達関数Ｈ（ｚ）が作用したものと言える。そして、信号歪み除去処理は、観測信号ｘ（ｔ）に対して逆フィルタＧ（ｚ）を作用させて復元信号ｙ（ｔ）を得る処理となる。このとき、信号歪み除去処理で得られた復元信号ｙ（ｔ）をフレーム分割して、それぞれに対して、それぞれの信号に基づいて求めた予測誤差フィルタ１−Ａｉ（ｚ）を適用して得られるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）はイノベーション系列ｅｉ（１），…，ｅｉ（Ｗ）に一致することが望ましい。もし、予測誤差フィルタ１−Ａｉ（ｚ）の出力信号ｄｉ（ｎ）がｄｉ（ｎ）＝ｅｉ（ｎ）〔１≦ｉ≦Ｆ，１≦ｎ≦Ｗ〕を満たすならば、［条件３］の条件下でΣｍ＝１ＭＨｍ（ｚ）Ｇｍ（ｚ）＝１となることが示せる（数学的証明については、参考文献Ａを参照されたい。）。つまり、ｓ（ｔ）＝ｙ（ｔ）が言える。このとき、１−Ａｉ（ｚ）は１−Ｂｉ（ｚ）に等しくなる。
（参考文献Ａ）  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情報量）を用いると、解くべき問題は式（７）のように定式化される。ただし、ｇ＝［ｇ１Ｔ，…，ｇＭＴ］Ｔ，ｇｍ＝［ｇｍ（０），…，ｇｍ（Ｌ）］Ｔ，ａ＝［ａ１Ｔ，…，ａＦＴ］Ｔ，ａｉ＝［ａｉ（１），…，ａｉ（Ｐ）］Ｔとし、ａｉ（ｋ）は予測誤差フィルタ係数を表す。Ｉ（Ｕ１，…，Ｕｎ）は確率変数Ｕｉ間の相互情報量を表す。またｇおよびａに記号＾を付したものは、得るべき最適解である。Ｔは転置を表す。
【数６】
2007094463.tif 000060
拘束条件
［１］  ‖ｇ‖＝１  （但し‖・‖はノルムを表す。）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。
【００３１】
相互情報量Ｉは、イノベーション推定値系列ｄ１（１），…，ｄ１（Ｗ），…，ｄｉ（１），…，ｄｉ（Ｗ），…，ｄＦ（１），…，ｄＦ（Ｗ）の振幅が定数倍されても変化しない。式（７）の拘束条件[１]は、この振幅の不定性を排除するための条件である。式（７）の拘束条件[２]は、上記[条件１]に対応して、予測誤差フィルタを最小位相系に制限するための条件である。以下、Ｉを、イノベーション推定値系列を入力としそれらの間の相互情報量を出力する関数と看做して、損失関数と呼称することにする。
【００３２】
１．６  損失関数の導出
式（７）の最適化を実行するためには、損失関数Ｉ（ｄ１（１），…，ｄＦ（Ｗ））を有限長の信号系列｛ｄｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝から推定しなければならない。（多変量）確率変数Ｕの微分エントロピーをＤ（Ｕ）と表記すると、Ｉ（ｄ１（１），…，ｄＦ（Ｗ））は式（８）で定義される。ただし、ｄ＝［ｄＦＴ，…，ｄ１Ｔ］Ｔ、ｄｉ＝［ｄｉ（Ｗ），…，ｄｉ（１）］Ｔである。
【数７】
2007094463.tif 000061
【００３３】
ｙ＝［ｙＦＴ，…，ｙ１Ｔ］Ｔ，ｙｉ＝［ｙｉ（Ｗ），…，ｙｉ（１）］Ｔとおくと、ｄはｙを用いて、ｄ＝Ａｙと表される。ただし、行列Ａは、式（９）および式（１０）で表される。
【数８】
2007094463.tif 000062
【００３４】
よって、Ｄ（ｄ）は式（１１）のように表される。
【数９】
2007094463.tif 000063
【００３５】
多変量確率変数Ｕの共分散行列をΣ（Ｕ）と表記すると、式（１１）右辺第二項について、Σ（ｄ）＝Ｅ｛ｄｄＴ｝＝ＡＥ｛ｙｙＴ｝ＡＴ＝ＡΣ（ｙ）ＡＴが成立するから、式（１２）が成り立つ。
【数１０】
2007094463.tif 000064
【００３６】
式（１１）、式（１２）を式（８）に代入すると、式（１３）を得る。ただし、σ（Ｕ）２は確率変数Ｕの分散を表す。
【数１１】
2007094463.tif 000065
【００３７】
式（１３）でＪ（Ｕ）は（多変量）確率変数Ｕのネゲントロピー（negentropy）である。ネゲントロピーはＵの非ガウス性の度合いを表す非負の値をとり、Ｕがガウス分布に従う場合に限り０をとる。Ｃ（Ｕ１，…，Ｕｎ）は式（１４）で定義される。Ｃ（Ｕ１，…，Ｕｎ）は確率変数Ｕｉ間の相関の度合いを表す非負の値をとり、これらが無相関の場合に限り０をとる。
【数１２】
2007094463.tif 000066
【００３８】
ところで、ｓ＝［ｓＦＴ，…，ｓ１Ｔ］Ｔ，ｓｉ＝［ｓｉ（Ｗ），…，ｓｉ（１）］ＴとおくとＪ（ｙ）＝Ｊ（ｓ）＝ｃｏｎｓｔａｎｔとなるため（証明略）、式（１３）は更に式（１５）のように簡単化できる。
【数１３】
2007094463.tif 000067
【００３９】
以上から、式（７）の最適化問題は式（１６）の最適化問題に等価である。
【数１４】
2007094463.tif 000068
拘束条件
［１］  ‖ｇ‖＝１  （但し‖・‖はノルムを表す。）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。
【００４０】
１．７  交代変数法による最適化
式（１６）について、交代変数の方法により、ｇとａを最適化する。すなわち、ｒ回目の繰り返しにおけるｇ及びａの推定値をそれぞれｇ＾（ｒ）、ａ＾（ｒ）と表せば、式（１７）および式（１８）の交互の最適化により更新された推定値ｇ＾（ｒ＋１），ａ＾（ｒ＋１）を得る。なお、ｇ＾およびａ＾は、記号＾がｇ、ａのそれぞれの上に付されたものを表す。例えば繰り返し回数の上限をＲ１とすれば、Ｒ１回目で得られるｇ＾（Ｒ１＋１）、ａ＾（Ｒ１＋１）が式（１６）の最適解である。上付き文字のＲ１は、Ｒ１である。
【数１５】
2007094463.tif 000069
拘束条件
［１］  ｇ＝ｇ＾（ｒ）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。
【数１６】
2007094463.tif 000070
拘束条件
［１］  ａ＝ａ＾（ｒ＋１）
［２］  ‖ｇ‖＝１
【００４１】
式（１７）の意図するところは、伝達特性を打ち消すための逆フィルタの現在の推定値に基づいて原信号に固有の特性を打ち消すための予測誤差フィルタを推定することである。他方、式（１８）の意図するところは、予測誤差フィルタの現在の推定値に基づいて逆フィルタを推定することである。イノベーション推定値系列ｄ１（１），…，ｄ１（Ｗ），…，ｄｉ（１），…，ｄｉ（Ｗ），…，ｄＦ（１），…，ｄＦ（Ｗ）が互いにより独立になるようにこれら２種類の最適化を繰り返すことで、逆フィルタと予測誤差フィルタを並行して推定することが可能になっている。したがって、ここでの繰り返しは逆フィルタの高精度な推定のために重要である。但し、図８から明らかなように処理する観測信号長が長くなる程、繰り返し回数は１回でも或る程度の信号歪み除去が達成されることが見て取れる。従って、この発明では、繰り返し回数は１回でもよい。
【００４２】
１．８  ａの最適化
本発明では、式（１７）の最適化を以下のように行う。
まず注意すべきことは、Ｃ（ｄ１（１），…，ｄＦ（Ｗ））はｄｉ（ｎ）の２次の統計量に関連するのに対して、Ｊ（ｄｉ（ｎ））はｄｉ（ｎ）の高次の統計量に関連する値である。２次の統計量は信号の振幅情報のみ提供するが、高次の統計量は位相情報も提供する。したがって、一般に、高次統計量を含む最適化は、非最小位相系を導く可能性がある。そこで、１−Ａｉ（ｚ）が最小位相系であるという拘束条件から、ａの最適化においては式（１９）の最適化問題を解く。
【数１７】
2007094463.tif 000071
拘束条件
［１］  ｇ＝ｇ＾（ｒ）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。
【００４３】
Ｃ（ｄ１（１），…，ｄＦ（Ｗ））は式（２０）で与えられる。
【数１８】
2007094463.tif 000072
【００４４】
ここで、行列Ａは式（９）および式（１０）に示すように上三角行列でその対角成分がすべて１であるから、ｌｏｇ ｄｅｔ Ａ＝０である。式（２０）を式（１２）に代入することで式（２１）の関係を得る。
【数１９】
2007094463.tif 000073
【００４５】
よって、式（１９）は、式（２２）の最適化問題と等価である。式（２２）は、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるａを求めよ、と云っている。
【数２０】
2007094463.tif 000074
拘束条件
［１］  ｇ＝ｇ＾（ｒ）
［２］  １−Ａｉ（ｚ）は、複素平面上で単位円の内部に全ての零点をもつ〔１≦ｉ≦Ｆ〕。
【００４６】
式（２２）で表される最適化問題を解くことは、観測信号にｇ＾（ｒ）で与えられる逆フィルタを適用して得られるアドホック信号に対して、各フレームにおいて線形予測分析を行うことと等価であり、必ず最小位相予測誤差フィルタを得ることができる。線形予測分析に関しては、上記参考文献１を参照されたい。
【００４７】
なお、式（２２）では、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大きい。この場合、対数関数は単調増加関数であるから、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めることができる。
【００４８】
１．９  ｇの最適化
本発明では、式（１８）の最適化を以下のように行う。
前述したとおり、Ｃ（ｄ１（１），…，ｄＦ（Ｗ））は｛ｄｉ（ｎ）；１≦ｉ≦Ｆ，１≦ｎ≦Ｗ｝の相関の度合いに関わる指標であるが、（ｒ＋１）回目のａの最適化において最小化されているため、Σｉ＝１ＦΣｎ＝１ＷＪ（ｄｉ（ｎ））に比べて無視できる。そこでｇの最適化においては、式（２３）の最適化問題を解く。
【数２１】
2007094463.tif 000075
拘束条件
［１］  ａ＝ａ＾（ｒ＋１）
［２］  ‖ｇ‖＝１
【００４９】
Ｊ（ｄｉ（ｎ））は、［条件２］に基づいて、式（２４）によって近似できる。この詳細は参考文献２を参照されたい。ただし、確率変数Ｕについて、κ４（Ｕ）はＵの尖度（４次キュムラント）を表す。式（２４）の右辺を第ｉフレームにおける正規化尖度という。
（参考文献２）  A.Hyvarinen, J.Karhunen, E.Oja, "INDEPENDENT COMPONENT ANALYSIS", John Wiley & Sons, Inc. 2001.
【数２２】
2007094463.tif 000076
【００５０】
［条件２］から音声信号のイノベーションの尖度は正であるため、κ４（ｄｉ（ｎ））／σ（ｄｉ（ｎ））４は正である。従って、式（２３）の最適化問題は、式（２５）の最適化問題に帰着する。σ（ｄｉ（ｎ）），κ４（ｄｉ（ｎ））は、［条件１］で述べた音声信号の局所的な定常性に基づいて、各フレーム内のサンプルから計算される。式（２６）では、１／Ｗを付加しているが、これは後の計算の便宜に過ぎず、式（２５）でｇの最適解を求めるにあたり影響を及ぼすものではない。式（２５）および式（２６）から、正規化尖度の全フレームでの総和が最大となるときのｇが、ｇ＾（ｒ＋１）となる。これらは、第ｉフレームにおける正規化尖度を全フレームで加算した値が最大となるｇを求めよ、と云っている。
【数２３】
2007094463.tif 000077
拘束条件
［１］  ａ＝ａ＾（ｒ＋１）
［２］  ‖ｇ‖＝１  
【００５１】
式（２５）に従ってｇの最適解を求めるには、Ｑをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（２７）で表される更新則に従って求められる。ｇ′をｇ′のノルムで除しているのは上記拘束条件[２]を課すためである。η（ｕ）は学習率を表す。ｕは、ｇの最適化における更新回数を表す。
【数２４】
2007094463.tif 000078
【００５２】
式（２７）において、∇Ｑｇは式（２８）および式（２９）で与えられる。
【数２５】
2007094463.tif 000079
【００５３】
式（２９）において、ｄｉ（ｎ）は式（３０）で、ｖｍｉ（ｎ）は式（３１）および式（３２）で与えられる。ｘｍｉ（ｎ）は、ｍ番目のマイクロホンで観測された観測信号のｉ番目のフレームの信号である。
【数２６】
2007094463.tif 000080
【００５４】
§２  二次統計量に基づく信号歪み除去
上述の従来的手法の信号歪み除去方法は、比較的長時間の観測信号（例えば２０秒程度である。）を要する。これは、一般に、正規化尖度のような高次統計量を計算するためには大量の観測信号のサンプルが必要となるからである。しかし、実際にはそうした長時間の観測信号を利用できない場合が多い。このため、従来的手法の信号歪み除去方法の適用
分野は極めて限られていた。
また高次統計量の計算は比較的複雑であるため、従来的手法の信号歪み除去方法では装置の構成が複雑になりやすい。
そこで、観測信号がより短時間（例えば３秒から５秒程度である。）の場合にも有効であり、かつ計算が従来に比して容易な信号歪み除去の原理を説明する。この原理は、信号の二次統計量のみを用いるものであり、§１で説明した本発明の基本原理より派生する。
【００５５】
２．１  二次統計量に基づく信号歪み除去の原理
二次統計量に基づく信号歪み除去では、上述の３つの条件に、次の２つの条件を設定する。
[条件４]  Ｍ≧２である。すなわち、複数本のマイクロホンを用いる。
[条件５]  Ｈｍ＝｛ｈｍ（ｋ）｝ｋ＝０Ｋは相異なるマイクロホンの間で共通の零点を持たない。
【００５６】
上記の式（１６）の最適化問題では、高次の統計量に関する値であるネゲントロピーＪおよび確率変数間の相関の度合いを示す指標Ｃを含む値を最小化するｇおよびａを求めた。
確率変数間の相関の度合いを示す指標Ｃは、二次の統計量で規定される。そこで、解くべき最適化問題を式（３３）で定式化する。
【数２７】
2007094463.tif 000081
【００５７】
式（２１）を参酌すれば、式（３３）の最適化問題は、式（３４）の最適化問題に転化される。式（３４）は、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値を全フレームで加算した値が最小となるｇおよびａを求めよ、と云っている。
【数２８】
2007094463.tif 000082
【００５８】
ところで、上記の[条件４]および[条件５]が成立する場合、多チャンネルの観測信号は、音源からの原信号によって駆動されるＡＲ系として捉えることができる（参考文献３参照）。このことは、逆フィルタＧの先頭タップを式（３５）のように固定できることを意味する。但し、ｍ＝１に相当するマイクロホンは、最も音源に近いマイクロホンである。（参考文献３）  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.
【数２９】
2007094463.tif 000083
【００５９】
式（３４）および式（３５）で規定されるｇを係数とする逆フィルタＧを、式（６）に従って観測信号ｘ（ｔ）に適用することで伝達特性が除去された復元信号ｙ（ｔ）を得る。
【００６０】
２．２  ａの最適化
式（３４）について、交代変数の方法により、ｇとａを最適化する。
逆フィルタの係数ｇｍ（ｋ）を固定した状態で予測誤差フィルタの係数ａｉ（ｋ）に関して式（３４）の損失関数を最小化する。
このとき、次の２点に注意する。１点目は、ｇ＝［ｇ１Ｔ，…，ｇＭＴ］Ｔは固定されているので、逆フィルタＧの出力である復元信号ｙ（ｔ）は予測誤差フィルタの最適化において不変であるということである。２点目は、ｉ番目のフレームにおける予測誤差フィルタの係数ａｉ（１），…，ａｉ（Ｐ）は、ｄｉ（１），…，ｄｉ（Ｗ）にのみ寄与するということである。
【００６１】
従って、各々のフレームにおいて、Σｎ＝１Ｗlogσ（ｄｉ（ｎ））２を最小化するように、予測誤差フィルタの係数ａｉ（１），…，ａｉ（Ｐ）を推定すればよい。[条件２]から第ｉフレームのイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散はフレーム内で定常であるから、Σｎ＝１Ｗlogσ（ｄｉ（ｎ））２の最小化は、Ｗ＊σ（ｄｉ（ｎ））２の最小化と等価である。記号＊は乗算を表す。分散σ（ｄｉ（ｎ））２は、＜ｄｉ（ｎ）２＞ｎ＝１Ｗとして計算できる。但し、＜ｄｉ（ｎ）２＞ｎ＝１Ｗは、１フレーム分のイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を使って計算したｄｉ（ｎ）の２乗平均を表す。結局、係数ａｉ（ｋ）は、Ｗ＊＜ｄｉ（ｎ）２＞ｎ＝１Ｗ、つまりｄｉ（ｎ）の２乗総和が最小になるときのａｉ（ｋ）として推定される。このような係数ａｉ（ｋ）は、線形予測分析の手法によって計算される。
【００６２】
なお、ここでは、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めるとして説明したが、これに限定する趣旨ではない。既述のとおり、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の全フレームでの総和が最小となるときのａをａ＾（ｒ＋１）として求めることができる。
【００６３】
２．３  ｇの最適化
予測誤差フィルタの係数ａｉ（ｋ）を固定した状態で逆フィルタの係数ｇｍ（ｋ）に関して式（３４）の損失関数を最小化する。
逆フィルタの係数ｇｍ（ｋ）に関する損失関数の最小化には、勾配法を用いる。[条件２]を用いると、式（３４）の最適化問題は、式（３６）の最適化問題に転化される。
【数３０】
2007094463.tif 000084
【００６４】
式（３６）に従ってｇの最適解を求めるには、Σｉ＝１Ｆlog＜ｄｉ（ｎ）２＞ｎ＝１
Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この解は、一般的には、式（３７）で表される更新則に従って求められる。δは学習率を表す。１≦ｍ≦Ｍ、１≦ｋ≦Ｌである。なお、式（３７）では、式（３５）の条件から‖ｇ‖＝１なる拘束条件が課されないことに留意しなければならない。同様に式（３５）の条件からｋが取る値の範囲は１≦ｋ≦Ｌである。
【数３１】
2007094463.tif 000085
【００６５】
上記式（２９）あるいは上記非特許文献１に記載された式（３）と比べて明らかなように、式（３７）の右辺第二項は、二次の統計量で表されており、この計算に高次の統計量を必要としない。このため、高次統計量を計算するには短い時間の観測信号の場合にも有効であり、計算自体も容易である。
【００６６】
なお、式（３６）では、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の対数値の全フレームでの総和が最小となるときのｇをｇ＾として求めるが、これに限定する趣旨ではない。上記各式では対数関数の底（base）を明記していないが、一般的には底を１０ないしネイピア数とするのが慣例であり、いずれにしても底は１よりも大である。この場合、対数関数は単調増加関数であるから、第ｉフレームそれぞれにおけるイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散の全フレームでの総和が最小となるときのｇをｇ＾として求めることができる。なお、この場合には式（３７）で示した更新則は適用できなくなり、改めてΣｉ＝１Ｆ＜ｄｉ（ｎ）２＞ｎ＝１Ｗをｇで微分してこれをゼロとしたときの解を求めればよい。この結果得られる更新則はＩＣＡと同様の考え方で定式化できるから略する。
【００６７】
§３  プリ・ホワイトニング
本発明の信号歪み除去には、プリ・ホワイトニングを適用することができる。観測信号をプリ・ホワイトニングすることで、最適化計算の安定化、とくにフィルタ係数の推定値の高速な収束が可能となる。
各マイクロホンで得られた観測信号系列全体｛ｘｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。
式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号ｘｍ（ｔ）に適用する。ｗｍ（ｔ）は、ｍ番目のマイクロホンで得た観測信号ｘｍ（ｔ）を白色化した信号を表す。
【数３２】
2007094463.tif 000086
【００６８】
このとき、式（３１）および式（３８）は式（４０）に、式（３２）は式（４１）に、変更すればよい。
【数３３】
2007094463.tif 000087
【００６９】
§４  実施形態
以下、本発明の実施形態を図面を参照して説明する。本発明の実施形態として、後述の各実施形態に限定するものではなく、各セクションで説明した原理を実現する実施形態であればよい。
【００７０】
＜第１実施形態＞
本発明の第１実施形態を実施する場合、以下の手順に従ってセンサで得た観測信号を処理する。ここでは、実施形態を具体的に説明する観点から信号として音声信号を例に挙げて説明する。
なお、第１実施形態の説明に先立ち、観測信号およびフレーム化処理について概説する。
【００７１】
((観測信号))
図示しないセンサ（例えばマイクロホン）によって得られたアナログ信号（このアナログ信号には伝達特性に由来する歪みが重畳されている。）は、例えば８,０００Ｈｚのサンプリングレートでサンプリングされ、適宜量子化された離散信号に変換される。以下、この離散信号を観測信号ということにする。アナログ信号から観測信号へのＡ／Ｄ変換などを実行するために必要となる構成要素（手段）は、いずれも公知技術の常套手段によって達成されるから、説明および図示を略する。
【００７２】
((フレーム化処理))
図示しない信号フレーム化手段が、離散信号から、時間軸方向に一定時間幅でフレームの始点を移動させながら、所定時間長の離散信号を切り出す。例えば２００サンプル点（８,０００Ｈｚ×２５ｍｓ）長の離散信号を、８０サンプル点（８,０００Ｈｚ×１０ｍｓ）ずつ始点を移動させながら切り出す。切り出された信号は、離散信号に公知の窓関数（例えば、ハミング窓、ガウス窓、方形窓など）が適用される。窓関数の適用によるフレーム化は公知の常套手段によって達成される。
【００７３】
本発明の第１実施形態である信号歪み除去装置（１）をコンピュータ（汎用機）で実現する場合のハードウェア構成例を説明する。
図２に例示するように、信号歪み除去装置（１）は、キーボード、ポインティングデバイスなどが接続可能な入力部（１１）、液晶ディスプレイ、ＣＲＴ（Cathode Ray Tube）ディスプレイなどが接続可能な出力部（１２）、信号歪み除去装置（１）外部に通信可能な通信装置（例えば通信ケーブル、ＬＡＮカード、ルータ、モデムなど）が接続可能な通信部（１３）、ＤＳＰ（Digital Signal Processor）（１４）〔ＣＰＵ（Central Processing Unit）でも良い。またキャッシュメモリやレジスタ（１９）などを備えていてもよい。〕、メモリであるＲＡＭ（１５）、ＲＯＭ（１６）やハードディスク、光ディスク、半導体メモリなどである外部記憶装置（１７）並びにこれらの入力部（１１）、出力部（１２）、通信部（１３）、ＤＳＰ（１４）、ＲＡＭ（１５）、ＲＯＭ（１６）、外部記憶装置（１７）間のデータのやり取りが可能なように接続するバス（１８）を有している。
また必要に応じて、信号歪み除去装置（１）に、ＣＤ−ＲＯＭ（Compact Disc Read Only
Memory）、ＤＶＤ（Digital Versatile Disc）などの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。
【００７４】
信号歪み除去装置（１）の外部記憶装置（１７）には、信号歪み除去のためのプログラムおよびこのプログラムの処理において必要となるデータ（観測信号）などが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶され、他のプログラムの処理に供されるときに、ＲＡＭや外部記憶装置などから読み込まれる。
【００７５】
より具体的には、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭなど〕には、観測信号に逆フィルタを適用する処理のためのプログラム、観測信号に逆フィルタを適用して得られた信号から予測誤差フィルタを求める処理のためのプログラム、予測誤差フィルタから逆フィルタを求める処理のためのプログラム、およびこれらのプログラムの処理において必要となるデータ（フレーム単位の観測信号など）が記憶されている。その他、これらのプログラムに基づく処理を制御するための制御プログラムも適宜に保存しておく。
【００７６】
第１実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、制御部）を実現することで、信号歪み除去が実現される。
【００７７】
そこで次に、図３〜図５を参照して、信号歪み除去装置（１）における信号歪み除去処理の流れを順次説明する。
大まかな処理の手順は、（ａ）観測信号ｘ（ｔ）に対して逆フィルタを適用した信号（以下、アドホック信号という。）を求め、（ｂ）アドホック信号から予測誤差フィルタを求め、（ｃ）この予測誤差フィルタから逆フィルタを求め、（ｄ）前記（ａ）、（ｂ）、（ｃ）の処理を繰り返して最適な逆フィルタを求め、（ｅ）最適化された逆フィルタを観測信号に対して適用した信号を復元信号ｙ（ｔ）として得る。
（ｂ）は上述のａの最適化に相当し、（ｃ）は上述のｇの最適化に相当し、（ｄ）は、式（１７）および式（１８）に相当する。（ｄ）の処理の繰り返し回数は予め定めた回数Ｒ１とする。つまり、１≦ｒ≦Ｒ１とする。また、（ｃ）の処理でｇを最適化する更新則の更新回数は予め定めた回数Ｒ２とする。つまり、１≦ｕ≦Ｒ２とする。（ｄ）の処理、つまり（ａ）、（ｂ）、（ｃ）の一連の処理を１回行うたびに、更新則によるＲ２回の更新が行なわれる。実施形態では、回数Ｒ１は、予め定めた回数とするが、これに限定されず、例えば、ｒ回目のｇを算出したときの式（２６）のＱの値とｒ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。同様に、回数Ｒ２は、予め定めた回数とするが、これに限定されず、例えば、ｕ回目のｇを算出したときの式（２６）のＱの値とｕ＋１回目のｇを算出したときの式（２６）のＱの値との差の絶対値が所定の正の微小値ε以下（あるいはε未満）になったときに繰り返しを中止するようにしてもよい。
【００７８】
（ステップＳ１００）
逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタを観測信号ｘ（ｔ）＝［ｘ１（ｔ），…，ｘｍ（ｔ），…，ｘＭ（ｔ）］Ｔに適用することで、アドホック信号ｙ（ｔ）を求める。アドホック信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、こ
こでは後述のＲ１回の処理を経て求められた復元信号ではないことを明示するためアドホック信号と呼称する。ここでｔは、全てのサンプル番号を示し、１≦ｔ≦Ｎとする。Ｎは全サンプル数である。第１実施形態では、マイクロホンの個数Ｍは１以上であればよい。
【数３４】
2007094463.tif 000088
【００７９】
逆フィルタの係数列｛ｇｍ（ｋ）；０≦ｋ≦Ｌ｝として、繰り返し回数Ｒ１の初回には予め定められた初期値を、２回目以降には後述する逆フィルタ計算部（１３）によって求められた逆フィルタｇ＾（ｒ＋１）を用いる。
【００８０】
（ステップＳ１０１）
予測誤差フィルタ計算部（１５）は、フレーム化処理を行うフレーム化処理部（１５１）とフレーム予測誤差フィルタ計算部（１５２）によって構成される。そして、フレーム予測誤差フィルタ計算部（１５２）は、第ｉフレームのアドホック信号から予測誤差フィルタを求める第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）からなる。ただし、ｉは、１≦ｉ≦Ｆを満たす整数である。
【００８１】
フレーム化処理部（１５１）は、逆フィルタ適用部（１４）で求められたアドホック信号｛ｙ（ｔ）；１≦ｔ≦Ｎ｝をフレーム化処理する。フレーム化処理は、例えば式（４３）のように、Ｗ点分を切り出す窓関数をＷ点ずつシフトさせて適用することにより行う。｛ｙｉ（ｎ）；１≦ｎ≦Ｗ｝はｉ番目のフレームに含まれるアドホック信号列を表す。
【数３５】
2007094463.tif 000089
【００８２】
そして、第ｉフレーム用予測誤差フィルタ計算部（１５２ｉ）は、式（２２）に従って、第ｉフレームのアドホック信号列｛ｙｉ（ｎ）；１≦ｎ≦Ｗ｝に対してＰ次の線形予測分析を行い、予測誤差フィルタの係数列｛ａｉ（ｋ）；１≦ｋ≦Ｐ｝を計算する。この算出方法は、上記参考文献１を参照されたい。ここで得られたａ１（１），…，ａ１（Ｐ），…，ａｉ（１），…，ａｉ（Ｐ），…，ａＦ（１），…，ａＦ（Ｐ）は、式（２２）のａ＾（ｒ＋１）を与える。
【００８３】
（ステップＳ１０２）
逆フィルタ計算部（１３）の機能構成例を、図４を参照して説明する。逆フィルタ計算部（１３）は勾配計算部（１３１）、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。更に、勾配計算部（１３１）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）とを備えて構成される。ここで更新用逆フィルタは、式（２７）のｇ〈ｕ〉に相当する。
【００８４】
第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖｍｉ（ｎ）を計算する（式（
３１）を参照）。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。
【００８５】
第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｙｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用してイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。爾後、第２の予測誤差フィルタ適用部（１３１２）は、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）を入力とする。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。
【００８６】
勾配ベクトル計算部（１３１３）は、信号ｖｍｉ（ｎ）とイノベーション推定値ｄｉ（ｎ）とを用いて現在の更新用逆フィルタｇ〈ｕ〉の勾配ベクトル∇Ｑｇを計算する（式（２８）および式（２９）を参照）。有限個のサンプルｖｍｉ（ｎ）およびｄｉ（ｎ）を用いて式（２９）を演算するときは、期待値Ｅをサンプルから求めればよい。ここで述べた処理の詳細の一例は、後述の第３実施形態の説明に譲る。
【００８７】
逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇ〈ｕ〉、学習率η（ｕ）、勾配ベクトル∇Ｑｇを用いて、式（２７）に従って、ｕ＋１回目の更新用逆フィルタｇ〈ｕ＋１〉を求める。式（２７）は、求められたｇ〈ｕ＋１〉を新たなｇ〈ｕ〉と見立てて更新を行なうことを意味する。
【００８８】
更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇ〈ｕ＋１〉、つまり新たなｇ〈ｕ〉および観測信号ｘ（ｔ）を用いて、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。つまり、式（４２）のｇｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。なお、更新用逆フィルタ適用後信号ｙ（ｔ）は、計算上は復元信号と全く同じであるが、ここでは後述のＲ１回の処理を経て求められた復元信号ではなく、更新則を行なうために算出される信号であることを明示するため更新用逆フィルタ適用後信号と呼称する。
【００８９】
制御部（６００）の制御によってＲ２回の更新が行なわれた結果として得られたｇ〈Ｒ２＋１〉は、式（２５）のｇ＾（ｒ＋１）に相当する。上付き文字のＲ２は、Ｒ２である。逆フィルタ計算部（１３）は、ｇ＾（ｒ＋１）を出力する。
【００９０】
制御部（５００）の制御によって、上述の一連の処理を１回行うごとにｒに１を加算してｒがＲ１に等しくなるまで、つまり上述の一連の処理をＲ１回繰り返すことで（ステップＳ１０３）、ｇ＾（Ｒ１＋１）を得る。上付き文字のＲ１は、Ｒ１である。このｇ＾（Ｒ１＋１）が、式（１６）の最適解とされる。そこで、ｇ＾（Ｒ１＋１）を得た段階で、逆フィルタ適用部（１４）は、式（４２）に従って逆フィルタｇ＾（Ｒ１＋１）を観測信号ｘ（ｔ）＝［ｘ１（ｔ），…，ｘＭ（ｔ）］Ｔに適用することで、復元信号ｙ（ｔ）を得ることができる（ステップＳ１０４）。
【００９１】
＜第２実施形態＞
第２実施形態は、第１実施形態の変形例に相当する。具体的には、§３で述べたプリ・ホワイトニングを行なう形態である。そこで、第１実施形態と異なる部分について図６および図７を参照して説明を加える。なお、プリ・ホワイトニングは観測信号に対して行なうプリ・プロセスであるから、ここで説明するプリ・ホワイトニングを行なう形態は、後述の第３実施形態にも適用可能である。
第２実施形態では、信号歪み除去装置（１）の外部記憶装置（１７）〔あるいはＲＯＭ
など〕に、白色化フィルタを求める処理のためのプログラム、白色化フィルタを観測信号に適用する処理のためのプログラムも記憶されている。
第２実施形態に係る信号歪み除去装置（１）では、外部記憶装置（１７）〔あるいはＲＯＭなど〕に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭ（１５）に読み込まれて、ＤＳＰ（１４）で解釈実行・処理される。その結果、ＤＳＰ（１４）が所定の機能（逆フィルタ適用部、予測誤差フィルタ計算部、逆フィルタ計算部、白色化フィルタ計算部、白色化フィルタ適用部）を実現することで、信号歪み除去が実現される。
【００９２】
（ステップＳ１００ａ）
白色化フィルタ計算部（１１）は、各マイクロホンで得られた観測信号全体｛ｘｍ（ｔ）；１≦ｔ≦Ｎ｝を白色化するフィルタ（白色化フィルタ）の係数｛ｆｍ（ｋ）；０≦ｋ≦Ｘ｝をＸ次の線形予測分析によって計算する。この計算は線形予測分析と同じであり上記参考文献１を参照されたい。白色化フィルタの係数は、白色化フィルタ適用部（１２）の入力となる。
【００９３】
（ステップＳ１００ｂ）
白色化フィルタ適用部（１２）は、式（３９）に従って、上記白色化フィルタを各マイクロホンで得られた観測信号に適用して、白色化信号ｗｍ（ｔ）を得る。既述のとおり、式（３１）は式（４０）に変更すればよいので、第１実施形態において、逆フィルタ計算部（１３）、とくに第１の予測誤差フィルタ適用部（１３１１）による処理を式（３１）ではなく式（４０）による計算処理に改めればよい。また、第１実施形態において、逆フィルタ適用部（１４）による処理を、式（４２）ではなく式（４４）による計算処理に改めればよい。ステップＳ１００ｂの処理の後、第１実施形態のステップＳ１００〜Ｓ１０４の処理を行うが、これらの処理では第１実施形態の各処理における観測信号をステップＳ１００ｂの処理で得られた白色化信号に読み替えて第１実施形態と同様の処理を行う。このことを明示するため、図７では、第１実施形態のステップＳ１００〜Ｓ１０４の各処理に相当する処理を示す符号に記号′を付している。
【数３６】
2007094463.tif 000090
【００９４】
＜実施例１＞
発明者らは第２実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、マイクロホンの数Ｍ＝４、白色化フィルタの次数Ｘ＝５００、逆フィルタの次数Ｌ＝１０００、窓関数の切出しサンプル数（１フレームのサンプル数）Ｗ＝２００、予測誤差フィルタの次数Ｐ＝１６、繰り返し回数Ｒ１＝１０、逆フィルタ計算部の更新回数Ｒ２＝２０とした。学習率η（ｕ）は、初期値を０．０５に設定し、もし式（２７）によって式（２６）の値が減少するならば、式（２６）の値が増大するまでη（ｕ）の値を再帰的に半減した。図６に示した逆フィルタ適用部（１４）へ入力する初期逆フィルタは、式（４５）のように設定した。
【数３７】
2007094463.tif 000091
【００９５】
本発明の第２実施形態の効果を、信号歪み除去の指標としてＤ５０値（インパルス応答
の全エネルギーに対する初期の５０ｍｓｅｃまでのエネルギーの比）を用いて評価した。連続発話データベースから男女各一名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
図８は、男声および女声について観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ１（図６に示す逆フィルタ適用部（１４）と、予測誤差フィルタ計算部（１５）と、逆フィルタ計算部（１３）を一巡する処理を実行して逆フィルタを求める回数）とＤ５０値の関係を示している。いずれの場合においても、繰り返し回数を増加させるとＤ５０値が向上しており、繰り返し処理の効果が顕著に見て取れる。特に観測信号長が５〜１０秒程度の比較的短い長さであっても、繰り返し処理によってＤ５０値が大きく向上したことが分かる。
【００９６】
また、本発明の第２実施形態の効果を、音声スペクトグラムの比較から検証した。
図９Ａは１分間の観測信号を用いて得られた残響を含まない音声（原音声）のスペクトログラムの抜粋、図９Ｂは１分間の観測信号を用いて得られた残響を含む音声（観測音声）のスペクトログラムの抜粋、図９Ｃは１分間の観測信号を用いて得られた残響除去後の音声（復元音声）のスペクトログラムの抜粋を示している。図９Ａと図９Ｃとの対比および図９Ｂと図９Ｃとの対比から、観測信号に含まれる残響が抑制され、原音声固有の特徴である調波構造やフォルマント構造が回復されたことが分かる。
【００９７】
また、本発明の第２実施形態の効果を、ＬＰＣスペクトル歪みを用いて検証した。
図１０Ｂは原音声の波形、図１０Ａは原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）および原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）を示している。図１０Ａおよび図１０Ｂの各横軸は秒単位の時間を表し、両図でタイムスケールを揃えている。図１０Ｂの縦軸は振幅値を表している。但し、原信号の相対的な大きさが分かればよいので、この趣旨から図１０Ｂの縦軸では単位を明示していない。図１０Ａの縦軸は、ＬＰＣスペクトル歪みＳＤ（ｄＢ）を表している。
図１０Ａから、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりも常に小さいことがわかる。なお、観測音声ではＬＰＣスペクトル歪みの平均が５．３９ｄＢ、分散が４．２０ｄＢであったのに対して、復元音声では平均が２．３８ｄＢ、分散が２．００ｄＢであった。
また、図１０Ａと図１０Ｂとの対比から、原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）が大きな値を示す区間（例えば約１．０秒〜約１．２秒の区間を参照）は、原音声の波形の振幅値がほぼ０であることがわかる。実際、この区間では発声がなく無音区間である。このため、実際に知覚される歪みはかなり小さくなっていた。つまり、発声区間における原音声と復元音声とのＬＰＣスペクトル歪みの時系列（図中の実線）は、原音声と観測音声とのＬＰＣスペクトル歪みの時系列（図中の点線）よりもかなり小さく、このため原音声のスペクトルを高い精度で復元できたことが結論付けられる。
【００９８】
＜第３実施形態＞
第３実施形態は、第１実施形態の変形例に相当する。具体的には、§２で述べた二次統計量に基づく信号歪み除去処理を行なう形態である。そこで、第１実施形態と異なる部分について図１１および図１２を参照して説明を加える。但し、第３実施形態では、マイクロホンの個数Ｍは２以上とする。
【００９９】
ステップＳ１００の処理およびステップＳ１０１の処理は、第１実施形態と同じである。
【０１００】
ステップＳ１０１の処理に続いて、ステップＳ１０２ａの処理を行う。
第３実施形態に係る逆フィルタ計算部（１３）の機能構成例を、図１１を参照して説明
する。
逆フィルタ計算部（１３）は、観測信号への予測誤差フィルタ適用部として機能する第１の予測誤差フィルタ適用部（１３１１）と、観測信号に更新用逆フィルタを適用して得られる信号（更新用逆フィルタ適用後信号）への予測誤差フィルタ適用部として機能する第２の予測誤差フィルタ適用部（１３１２）と、勾配ベクトル計算部（１３１３）と、逆フィルタ更新部（１３２）および更新用逆フィルタ適用部（１３３）によって構成される。ここで更新用逆フィルタは、式（３７）のｇｍ（ｋ）に相当する。
【０１０１】
第１の予測誤差フィルタ適用部（１３１１）は、ｍ番目〔１≦ｍ≦Ｍ〕のマイクロホンで観測された観測信号ｘｍ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｘｍｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用して予測誤差フィルタ適用後の信号ｖｍｉ（ｎ）を計算する（式（３８）を参照）。具体的には、フレーム化処理部（４０２Ｂ）が、入力された観測信号ｘｍ（ｔ）に対してフレーム化処理を行い、観測信号ｘｍ（ｔ）のｉ番目のフレームの信号ｘｍｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０４ｉ）が信号ｘｍｉ（ｎ）を入力として、式（３８）に従って信号ｖｍｉ（ｎ）を出力する。但し、１≦ｉ≦Ｆである。
【０１０２】
第２の予測誤差フィルタ適用部（１３１２）は、更新用逆フィルタ適用後信号ｙ（ｔ）をフレーム化して、各フレームにつき、ｉ番目のフレームの信号ｙｉ（ｎ）に対してステップＳ１０１の処理で得られたｉ番目の予測誤差フィルタａｉ（ｋ）を適用してイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を計算する（式（３０）を参照）。なお、更新用逆フィルタ適用後信号ｙ（ｔ）の初期値は、ステップＳ１００の処理で得られた信号とすればよい。具体的には、フレーム化処理部（４０２Ａ）が、初期値の場合を除き、後述する更新用逆フィルタ適用部（１３３）が出力した更新用逆フィルタ適用後信号ｙ（ｔ）に対してフレーム化処理を行い、ｉ番目のフレームの信号ｙｉ（ｎ）を出力する。そして、予測誤差フィルタ適用部（４０３ｉ）が信号ｙｉ（ｎ）を入力として、式（３０）に従ってイノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）を出力する。但し、１≦ｉ≦Ｆである。
【０１０３】
勾配ベクトル計算部（１３１３）は、信号ｖｍｉ（ｎ）とイノベーション推定値ｄｉ（ｎ）とを用いて現在の更新用逆フィルタｇｍ（ｋ）の勾配ベクトルを計算する（式（３７）の右辺第二項を参照）。具体的には、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、相互相関計算部（４０５ｉ）は信号ｖｍｉ（ｎ）とイノベーション推定値ｄｉ（ｎ）との相互相関〈ｄｉ（ｎ）ｖｍｉ（ｎ−ｋ）〉ｎ＝１Ｗを計算する。また、各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、分散計算部（４０６ｉ）は、イノベーション推定値ｄｉ（１），…，ｄｉ（Ｗ）の分散〈ｄｉ（ｎ）２〉ｎ＝１Ｗを求める。各フレーム番号ｉ（１≦ｉ≦Ｆ）に関して、除算部（４０７ｉ）は、〈ｄｉ（ｎ）ｖｍｉ（ｎ−ｋ）〉ｎ＝１Ｗ／〈ｄｉ（ｎ）２〉ｎ＝１Ｗを求める。加算部（４０８）は、除算部（４０７１）〜（４０７Ｆ）の出力の全フレームに亘る総和、つまり式（３７）の右辺第二項を求める。
【０１０４】
逆フィルタ更新部（１３２）は、現在の更新用逆フィルタｇｍ（ｋ）、学習率δ、勾配ベクトルを用いて、式（３７）に従って、ｕ＋１回目の更新用逆フィルタｇｍ（ｋ）′を求める。式（３７）は、求められたｇｍ（ｋ）′を新たなｇｍ（ｋ）と見立てて更新を行なうことを意味する。
【０１０５】
更新用逆フィルタ適用部（１３３）は、逆フィルタ更新部（１３２）によって得られたｇｍ（ｋ）′、つまり新たなｇｍ（ｋ）および観測信号ｘ（ｔ）を用いて、式（４２）に従って、更新用逆フィルタ適用後信号ｙ（ｔ）を求める。具体的には、式（４２）のｇｍ（ｋ）としてｕ＋１回目の更新で得られたｇを用いて計算する。この計算で得られた更新
用逆フィルタ適用後信号ｙ（ｔ）は、第２の予測誤差フィルタ適用部（１３１２）の入力となる。
【０１０６】
ステップＳ１０２ａの処理に続いて、ステップＳ１０３およびステップＳ１０４の処理を行うが、第１実施形態と同じであるから説明を略する。
【０１０７】
＜実施例２＞
発明者らは第３実施形態の実証実験を行ったので、その実験結果を示す。実験条件として、Ｍ＝４，Ｌ＝１０００，Ｗ＝２００，Ｐ＝１６，Ｒ１＝６，Ｒ２＝５０とした。学習率δは、初期値を０．０５に設定し、Σｉ＝１Ｆlog＜ｄｉ（ｎ）２＞ｎ＝１Ｗの値が増加するならば、Σｉ＝１Ｆlog＜ｄｉ（ｎ）２＞ｎ＝１Ｗが減少するまで、学習率δの値を順次半減させた。逆フィルタの初期推定値は、ｇｍ（ｋ）＝０，１≦ｍ≦Ｍ，１≦ｋ≦Ｌとして設定した。
【０１０８】
本発明の第３実施形態の効果を、音声明瞭度を表すＲＡＳＴＩ（参考文献５を参照）を残響除去の指標として評価した。連続発話データベースから男女各五名の発話を取り出し、残響時間０．５秒の残響室で測定したインパルス応答を畳み込むことで観測信号を合成した。
（参考文献５）  H. kuttruff. Room acoustics. Elsevier Applied Science, third edition, P.237 1991.
【０１０９】
図１３は、Ｎを３秒、４秒、５秒、１０秒とする各観測信号のＲＡＳＴＩ値を表示したものである。図１３に示すように、観測信号が３〜５秒のように短時間の場合でも、高い残響除去性能を示していることが分かる。
図１４は、残響除去前後におけるエネルギー減衰曲線の例である。直接音が到達してから５０ミリ秒後の反射音のエネルギーが１５ｄＢ低減されていることが分かる。
【産業上の利用可能性】
【０１１０】
本発明は、様々な信号処理システムの性能向上に寄与する要素技術であるところ、例えば音声認識システム、テレビ会議システム、補聴器、音楽情報処理システム等に利用することができる。
【図面の簡単な説明】
【０１１１】
【図１】本発明の原理を説明するためのモデル機構を表したブロック線図。
【図２】第１実施形態に係る信号歪み除去装置（１）のハードウェア構成例を示す図。
【図３】第１実施形態に係る信号歪み除去装置（１）の機能構成例を示す機能ブロック図。
【図４】信号歪み除去装置（１）の逆フィルタ計算部（１３）の機能構成例を示す機能ブロック図。
【図５】第１実施形態における信号歪み除去処理の流れを示す処理フロー図。
【図６】第２実施形態に係る信号歪み除去装置（１）の機能構成例を示す機能ブロック図。
【図７】第２実施形態における信号歪み除去処理の流れを示す処理フロー図。
【図８】観測信号長Ｎを５秒、１０秒、２０秒、１分、３分に変化させたときの、繰り返し回数Ｒ１とＤ５０値の関係を示す図。
【図９】Ａは残響を含まない音声のスペクトログラム、Ｂは残響を含む音声のスペクトログラム、Ｃは残響除去後の音声のスペクトログラム。
【図１０】Ａは残響除去音声のＬＰＣスペクトル歪みの時間変動を説明するためのグラフ、Ｂは対応する区間における原音声信号の抜粋。
【図１１】第３実施形態に係る信号歪み除去装置（１）の逆フィルタ計算部（１３）の機能構成例を示す機能ブロック図。
【図１２】第３実施形態における信号歪み除去処理の流れを示す処理フロー図。
【図１３】３秒、４秒、５秒、１０秒の各観測信号のＲＡＳＴＩ値を表示した図。
【図１４】残響除去前後におけるエネルギー減衰曲線の例を示した図。
【図１５】従来技術を説明するための機能ブロック図。
【手続補正６】
【補正対象書類名】図面
【補正対象項目名】図１１
【補正方法】変更
【補正の内容】
【図１１】
2007094463.tif 000092
【国際調査報告】
2007094463.tif 000093
2007094463.tif 000094
2007094463.tif 000095
2007094463.tif 000096