(57)【特許請求の範囲】
【請求項１】
時不変又は時変システムの高速実時間同定を行うシステム同定方法において、次式で表されるＨ∞フィルタ方程式を用い、
【数１】
0004067269.tif 000046
Ｈ∞評価基準として、｛  （フィルタ誤差に対応する項／評価関数の重み（ρ））  ／  ［（初期状態の誤差に対応する項）＋（システム雑音に対応する項）＋（観測雑音に対応する項／評価関数の重み（ρ））］  ｝の最大値が、予め与えられた上限値（γｆ２）より小さく抑えるように定めることにより、外乱に対して頑健なフィルタリングアルゴリズムとしたシステム同定方法。
ただし、
ｘｋ：  状態ベクトルまたは単に状態  ；  未知、これが推定の対象となる。
ｘ０：  初期状態  ；  未知である。
ｗｋ：  システム雑音  ；  未知である。
ｖｋ：  観測雑音  ；  未知である。
ｙｋ：  観測信号  ；  フィルタの入力となり、既知である。
ｚｋ：  出力信号  ；  未知である。
Ｈｋ：  観測行列  ；  既知である。
ｘ＾ｋ｜ｋ：  観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの状態推定値  ；  フィルタ方程式（１２）によって与えられる。
ｘ＾０｜０：状態の初期推定値  ；  本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ＋１：フィルタゲイン  ；行列Ｐ＾ｋ＋１｜ｋから得られる。
Σｗｋ  ：  システム雑音の共分散行列に対応  ；  理論上は既知であるが、フィルタの実行前には未知である。
Ｐ＾ｋ｜ｋ−１：  ｘ＾ｋ｜ｋ−１の誤差の共分散行列に対応  ；  リカッチ方程式によって与えられる。
Ｐ＾１｜０：  初期状態の誤差の共分散行列に対応  ；  本来未知であるが、便宜上ε０Ｉが用いられる。
さらに、
「フィルタ誤差」：  
（１１）式のｚｖｋ｜ｋと（９）式のｚｋ＝Ｈｋｘｋの差（（１５）式のｅｆ，ｉに対応）の重み付きノルムの２乗の和。
「初期状態の誤差」：
時刻ｋ＝０のときの、以下に示す状態空間モデルの、（７）式の状態ｘkの初期値ｘ０とその推定値の差の重み付きノルム。
「システム雑音に対応する項」：
以下に示す状態空間モデルの、（７）式のシステム雑音ｗｋの重み付きノルムの２乗の和。
「観測雑音に対応する項」：
以下に示す状態空間モデルの、（８）式の観測雑音ｖｋの重み付きノルムの２乗の和。
状態空間モデル
【数２】
0004067269.tif 000047
【請求項２】
さらに、次式により時刻ｋの状態推定値ｘ＾ｋ｜ｋから出力信号ｚｋを求めるようにした請求項１に記載のシステム同定方法。
【数３】
0004067269.tif 000048
【請求項３】
ゲイン行列Ｋｋ、補助変数Ａｋ、Ｓｋ、Ｄｋ、状態推定値ｘ＾ｋ｜ｋの再帰式の初期条件をそれぞれ以下のように定めるステップと、
【数４】
0004067269.tif 000049
時刻ｋにおける補助変数ＡｋとＳｋを以下のように再帰的に決定するステップと、
【数５】
0004067269.tif 000050
ゲイン行列Ｋｋに補助変数を含む行を増やした第２のゲイン行列Ｋ∪ｋを以下のように計算するステップと、
【数６】
0004067269.tif 000051
第２のゲイン行列Ｋ∪ｋを以下のように分割し、第１及び第２の分割ゲイン行列を求めるステップと、
【数７】
0004067269.tif 000052
第１及び第２の分割ゲイン行列を含む次式により、Ｄｋを決定し、時刻ｋ＋１におけるゲイン行列Ｋｋ＋１を求め、時刻ｋ＋１におけるフィルタゲインＫｓ，ｋ＋１を求めるステップと、
【数８】
0004067269.tif 000053
ここで、ηｋ∈Ｒ２×１，Ｄｋ∈ＲＮ×１，
Ｋｋ＋１∈ＲＮ×２，Ｋｓ，ｋ＋１∈ＲＮ×１，
０＜ρ＝１−γｆ−２≦１，γｆ＞１
上記で求められたフィルタゲインＫｓ，ｋ＋１に基づき、前記フィルタ方程式（１２）を以下のように更新するステップと、
【数９】
0004067269.tif 000054
時刻を進ませて、各前記ステップを繰り返すためのステップと
を含む請求項２に記載のシステム同定方法。
【請求項４】
さらに、高速処理に適した存在条件として、次式を用いることにより、計算量Ο（Ｎ）で前記高速Ｈ∞フィルタの存在性を検査することを特徴とする請求項１乃至３のいずれかに記載のシステム同定方法。
【数１０】
0004067269.tif 000055
【請求項５】
前記Ｈ∞フィルタ方程式を適用し、状態推定値ｘ＾ｋ｜ｋを求め、擬似エコーを次式のように推定し、求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現することを特徴とする請求項１又は３又は４に記載のシステム同定方法。
【数１１】
0004067269.tif 000056
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、システム同定方法に係り、特に、新たなＨ∞評価基準に基づいて開発された変形Ｈ∞フィルタの高速Ｈ∞フィルタリングアルゴリズムを用いて時変システムを高速に実時間同定するものである。
【０００２】
【従来の技術】
一般に、システム同定(system identification)とは、入出力データに基づいてシステムの入出力関係の数理モデル(伝達関数、あるいはインパルス応答など)を推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
【０００３】
（基本原理）
図１４に、システム同定のための構成図を示す。このシステムは、未知システム１、適応フィルタ２を備える。また、適応フィルタ２は、ＦＩＲディジタルフィルタ３、適応アルゴリズム４を有する。
以下に、未知システム１を同定する出力誤差方式の一例を説明する。ここで、ｕｋは未知システム１の入力、ｄｋは所望信号であるシステムの出力、ｄ＾ｋはフィルタの出力である。（なお、「＾」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。）
【０００４】
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差ｅｋ＝ｄｋ−ｄ＾ｋを最小にするように適応アルゴリズムによってＦＩＲディジタルフィルタ３の係数を調節する。
【０００５】
図１５に、インパルス応答の調節機構についての構成図を示す。
ここで、適応アルゴリズムの一例として、その簡便さより、次のＬＭＳアルゴリズム(least mean square algorithm)が広く用いられている。
【０００６】
【数６】
0004067269.tif 000002
【０００７】
また、一般に、時変システムの同定には、例えば、収束性が早いカルマンフィルタが適している。
【０００８】
【数７】
0004067269.tif 000003
【０００９】
ここで、求めるべきインパルス応答{ｈｉ}は状態推定値ｘ＾ｋ｜ｋとして得られ、システムへの入力{ｕｋ}は観測行列Ｈｋの要素として用いられている。
また、σ２ｗ＝０のときのカルマンフィルタに対して、観測行列Ｈｋのシフト特性（Ｈｋ＋１（ｉ＋１）＝Ｈｋ（ｉ））を利用して、単位時間ステップ当たりの計算量をＮに比例した演算回数、すなわちＯ(Ｎ)までに軽減した高速カルマンフィルタリングアルゴリズムが知られている。詳細は、１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」など参照。
【００１０】
（エコーキャンセラへの適用例）
国際電話など長距離電話回線では,信号増幅などの理由から４線式回線が用いられている。一方、加入者回線は比較的短距離なので、２線式回線が使用されている。
図１６に通信系とエコーについての説明図を示す。２線式回線と４線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Ｂからの信号(音声)は話者Ａのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は４線式回線に漏れ、増幅された後、再び受信者(話者Ａ)に戻ると云った現象が起こる。これがエコー（echo）である。エコーは、伝送距離が長くなるにつれて(遅延時間が長くなるにつれて)影響が大きくなり、著しく通話の品質を劣化させる(パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する)。
【００１１】
図１７に、エコーキャンセラの原理図を示す。
そこで、図示のようにエコーキャンセラ(echo canceller) を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
【００１２】
エコーパスのインパルス応答の推定は、残留エコーｅｋの平均２乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Ａからの信号(音声)である。一般に、話者２人が同時に話し始めた(ダブルトーク)ときはインパルス応答の推定を中断する。また、ハイブリッドトランスのインパルス応答長は50［ｍｓ］程度なので、サンプリング周期を125［μｓ］とするとエコーパスのインパルス応答の次数は実際は400程度となる。
【００１３】
【発明が解決しようとする課題】
従来技術では、適応アルゴリズムとして、その簡便さより、ＬＭＳアルゴリズム(least mean square algorithm)が広く用いられて来たが、収束性が非常に遅いため急激に変化するような時変システムの同定は不可能であった。
【００１４】
また、追従性の優れたカルマンフィルタは、計算量がＯ（Ｎ２）あるいはＯ（Ｎ３）であり、それが、タップ数Ｎと共に急速に増加してしまい、高いタップ数Ｎが要求される現実の問題の実時間処理は困難であった。この対策として、タップ数Ｎに対して単位時間ステップ当たり計算量Ｏ（Ｎ）で同定可能な高速カルマンフィルタが提案されているが、その定常的な特性(システム雑音が考慮できない点)から時変システムの同定は不可能であった。
【００１５】
本発明は、以上の点に鑑み、新たなＨ∞評価基準に基づいて開発した変形Ｈ∞フィルタの高速Ｈ∞フィルタリングアルゴリズムを用いて、時不変および時変システムの高速実時間同定および推定を実現することを目的とする。また、本発明は、本アルゴリズムの特殊な場合として高速カルマンクイルタリングアルゴリズムを含み、また、時変システムの追従性を支配するシステム雑音の共分散を理論的に決定することを目的とする。また、本発明は、突然回線が切り替わるような激しく変化する時変システムのエコーキャンセラなどのように、入力信号が不連続に変化する場合においても、非常に有効な高速時変システム同定方法を提供することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム同定方法を提供することを目的とする。
【００１６】
【課題を解決するための手段】
本発明では、上述の課題を解決するために、新たにＨ∞評価基準を考案し、これに基づく変形Ｈ∞フィルタの高速アルゴリズムを開発すると共に、この高速Ｈ∞フィルタリングアルゴリズムに基づく高速時変システム同定方法を提案する。本発明による高速アルゴリズムは、単位時間ステップ当たり計算量Ｏ（Ｎ）で急激に変化する時変システムの追跡が可能である。また、γｆ＝∞の極限で高速カルマンフィルタリングアルゴリズムと完全に一致すると云った便利な特性をもっている。
【００１７】
【発明の実施の形態】
以下に、本発明の実施の形態について説明する。なお、詳細は、例えば、K. Nishiyama :"Derivation of A Fast Algorithm of Modified H∞ Filters", IEEE international Conference on Industrial Electronics, Control and Instrumentation, October, 2000  に示されている。
【００１８】
１．記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
ｘｋ：  状態ベクトルまたは単に状態  ；  未知、これが推定の対象となる。
ｘ０：  初期状態  ；  未知である。
ｗｋ：  システム雑音  ；  未知である。
ｖｋ：  観測雑音  ；  未知である。
ｙｋ：  観測信号  ；  フィルタの入力となり、既知である。
ｚｋ：  出力信号  ；  未知である。
Ｈｋ：  観測行列  ；  既知である。
Ｌｋ：  出力行列  ；  既知である。
ｘ＾ｋ｜ｋ：  観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの状態推定値；  フィルタ方程式によって与えられる。
ｘ＾０｜０：状態の初期推定値  ；  本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ＋１：フィルタゲイン  ；行列Ｐ＾ｋ｜ｋ−１から得られる。
Σｗｋ  ：  システム雑音の共分散行列に対応  ；  理論上は既知であるが、実際には未知である。
Ｐ＾ｋ｜ｋ−１：  ｘ＾ｋ｜ｋ−１の誤差の共分散行列に対応  ；  リカッチ方程式によって与えられる。
Ｐ＾１｜０：  初期状態の誤差の共分散行列に対応  ；  本来未知であるが、便宜上ε０Ｉが用いられる。
σ2ｖ  ：  観測雑音の分散  ；  理論上は既知として扱われるが、実際には未知である。
σ2ｗ  ：  システム雑音の分散  ；  理論上は既知として扱われるが、実際には未知である。
【００１９】
なお、記号の上に付される”＾”は、推定値の意味であり、"∪"は、行列を１行増やしたことを表す。”また、”〜”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、Ｌ、Ｈ、Ｐ、Ｋ等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
２．変形Ｈ∞フィルタ
【００２０】
つぎに、次式（７）〜（９）のような状態空間モデルを考える。
【数８】
0004067269.tif 000004
【００２１】
ここで、エコーキャンセラなどを想定し、Ｌｋ＝Ｈｋ（Ｈｋ＝[ｕｋ  ｕｋ−２・・・  ｕｋ−Ｎ＋１]）  とする。このような状態空間モデルに対して、次式（１０）のようなＨ∞評価基準（新たに左辺にγｆが入っている）を提案する。
【００２２】
【数９】
0004067269.tif 000005
【００２３】
この評価基準を満たすレベルγｆの変形Ｈ∞フィルタは、ρやΣｗｋがγｆに依存しないと仮定すれば、システム同定の分野で通常知られる方法を適用することによって、次の式（１１）〜式(１４)で与えられる。なお、その通常知られる方法として、例えば、
B. Hassibi, A. H. Sayed, and T. Kailath:  "Linear Estimation in Krein Spaces - Part I: Theory,"  IEEE Trans. Automatic Control, 41, 1,  pp.18-33,1996.、
B. Hassibi, A. H. Sayed, and T. Kailath:  "Linear Estimation in Krein Spaces - Part II: Applications,"  IEEE Trans. Automatic Control, 41, 1,  pp.34-49,1996.
等を参照のこと。
【００２４】
【数１０】
0004067269.tif 000006
【００２５】
また、評価基準の重みρは予めに決められた上限値γｆに依存するので、上述のアルゴリズムは通常のＨ∞フィルタとは本質的に異なる。本アルゴリズムはρで重みづけされた外乱(初期状態ｘ０,システム雑音｛ｗｉ｝,観測雑音｛υｉ｝からフィルタ誤差｛ｅｆ，ｉ｝への最大エネルギーゲインをγｆ２より小さく押えているので、外乱に対してロバスト(頑健)なフィルタリングアルゴリズムとなる。この性質が時変システムの追従特性に反映される。また、γｆ→∞のとき、ρ＝１、Σｗｋ＝０となり、変形Ｈ∞フィルタは通常のカルマンフィルタと一致する。
【００２６】
変形Ｈ∞フィルタの主な計算上の負担は、Ｎ２またはＮ３に比例する計算量を必要とするＰ＾ｋ＋１｜ｋ∈ＲＮ×Ｎの更新の際に生じる。すなわち、単位時間ステップ当たりＯ（Ｎ２）の算術演算が必要となる。ここで、タップ数Ｎは状態ベクトルｘｋの次元と一致する。ゆえに、ｘｋの次元が増加するにつれて変形Ｈ∞フィルタの実行に要する計算時間は急速に増大する。この計算上の負担を克服するために変形Ｈ∞フィルタの高速アルゴリズムの導出が必要となる。
【００２７】
３．高速Ｈ∞フィルタリングアルゴリズム
変形Ｈ∞フィルタの計算量は、式（１４）のリカッチ方程式（誤差の共分散方程式）の計算に支配される。よって、変形Ｈ∞フィルタを高速に処理するためには、リカッチ方程式を用いずに式（１３）のフィルタゲインを直接決定できれば、大幅に計算量を削減できる。
しかし、フィルタゲインＫｓ，ｋ∈ＲＮ×１を求める高速アルゴリズムの導出は困難なため、次のように定義されるゲイン行列を高速に計算するアルゴリズムを導出することを考える。
【００２８】
【数１１】
0004067269.tif 000007
【００２９】
このとき、次の補題が成り立つ。
補題1
行列Ｐｋは式（１４）のリカッチ方程式を満たす。これより、ゲイン行列Ｋｋが求まれば、次の補題からフィルタゲインＫｓ，ｋが直ちに得られる。
【００３０】
補題2
変形Ｈ∞フィルタのフィルタゲインはＫｓ，ｋは、ゲイン行列Ｋｋを用いて次のように得られる。実際、ゲイン行列Ｋｋは次の再帰的方法によって高速に計算できる。
【００３１】
【数１２】
0004067269.tif 000008
【００３２】
補題3
ゲイン行列Ｋｋは次のように更新される。
【００３３】
【数１３】
0004067269.tif 000009
【００３４】
ここで、ｍｋ∈ＲＮ×２とμｋ∈Ｒ１×２は行列Ｋ∪ｋ＝Ｑ∪ｋ−１Ｃ∪ｋの次の分割によって得られる。
【００３５】
【数１４】
0004067269.tif 000010
【００３６】
また、補助変数Ａｋ∈ＲＮ×１，Ｓｋ∈ＲおよびＢｋＦ−１ｋ∈ＲＮ×１も同様に得られる。
結論として、高速Ｈ∞フィルタリングアルゴリズムは以下のように要約することができる。
図１に、高速アルゴリズムのフローチャートを示す。なお、Lは最大データ数を示す。
[ ステップ０]  再帰式の初期条件を以下のようにする。ここでε０は充分に大きい正の定数である。
【００３７】
【数１５】
0004067269.tif 000011
【００３８】
[ ステップ１]  時刻ｋと最大データ数Lとを比較する。時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。）
[ ステップ２]  ＡｋとＳｋを以下のように再帰的に決定する。
【００３９】
【数１６】
0004067269.tif 000012
【００４０】
[ ステップ３]  Ｋ∪ｋを以下のように計算する。
【００４１】
【数１７】
0004067269.tif 000013
【００４２】
[ ステップ４]  Ｋ∪ｋを以下のように分割する。
【００４３】
【数１８】
0004067269.tif 000014
【００４４】
[ ステップ５]  Ｄｋを決定し、Ｋｋ＋１を通して、ゲイン行列Ｋｓ，ｋ＋１を以下のように得る。
【００４５】
【数１９】
0004067269.tif 000015
【００４６】
ここで、ηｋ∈Ｒ２×１，Ｄｋ∈ＲＮ×１，
Ｋｋ＋１∈ＲＮ×２，Ｋｓ，ｋ＋１∈ＲＮ×１，
０＜ρ＝１−γｆ−２≦１，γｆ＞１  である。
[ ステップ６]  Ｈ∞フィルタのフィルタ方程式を以下のように更新する。
【００４７】
【数２０】
0004067269.tif 000016
【００４８】
[ ステップ７]  時刻ｋを進ませて（ｋ＝ｋ＋１）、 ステップ２に戻り、データがある限り続ける。
（高速処理に適した存在条件：補題６）
また、次の存在条件を用いれば、計算量Ο（Ｎ）で高速Ｈ∞フィルタの存在性が検査できる。
【００４９】
【数２１】
0004067269.tif 000017
【００５０】
４．本高速アルゴリズムの計算量
つぎに、高速Ｈ∞フィルタリングアルゴリズムの計算量が、変形Ｈ∞フィルタリングアルゴリズムの計算量と比べて如何に減少するかを考察する。ここで、式の計算量の評価は乗算回数のみに注目し、以下のような方法で計算した。
（Ｊ×Ｋ行列）×（Ｋ×Ｌ行列）の乗算回数＝Ｊ×Ｋ×Ｌ（回）
ただし、行列やベクトルが３つ以上乗算されるときは特に図に示さない限り左から計算されるものとする。
【００５１】
（変形Ｈ∞フィルタリングアルゴリズムの計算量）
図２及び図３に、変形Ｈ∞フィルタリングアルゴリズムの各部分の計算量の説明図を示す。ただし、Ｎはタップ数（状態ベクトルｘｋのディメンジョン）である。ここで、図３（ａ）においてＲｅ，ｋからＲｅ，ｋ−１を求めるための計算は無視する。同様に、図２（ａ）において（Ｈｋ＋１Ｐ＾ｋ＋１|ｋＨＴｋ＋１＋１）の部分から（Ｈｋ＋１Ｐ＾ｋ＋１|ｋＨＴｋ＋１＋１）−１を求める計算も無視する。
【００５２】
図２（ａ）、図３（ａ）および図３（ｂ）より、Ｋｓ|ｋ＋１、Ｒｅ，ｋおよびＰ＾ｋ＋１|ｋは計算量がタップ数の２乗に比例することがわかる。よって、変形Ｈ∞フィルタリングアルゴリズム全体の計算量は単位時間ステップ当りΟ（Ｎ２）となる。
【００５３】
また、図４に、行列計算の順番を変えた場合の計算量の説明図を示す。すなわち、リカッチ方程式において、右辺第２項の部分の行列計算の順番を変えた場合の計算量を図４に示す。
上述の部分の計算量がインパルス応答の次数の３乗に比例するので、Ｐ＾ｋ＋１|ｋの計算量もインパルス応答の次数の３乗に比例する。これに伴い、Ｈ∞フィルタ全体の計算量もタップ数の２乗から３乗にまで増加する。
【００５４】
しかし、いずれのアルゴリズムにせよタップ数の２乗または３乗に比例する計算量をもつので、タップ数の増加と共にフィルタの実行にかかる計算上の負担が著しく増加する。実際、通信工学の分野に利用する場合などはタップ数が例えば４００程度であるため、このアルゴリズムでの実用はかなり困難なものとなる。
【００５５】
（高速Ｈ∞フィルタリングアルゴリズムの計算量）
次に、図５及び図６に、高速Ｈ∞フィルタリングアルゴリズムの計算量の説明図を示す。ここで、図５（ｂ）のＫｕｋの式においてＳｋ−１はＳｋから求められるが、その計算は無視する。同様に、図６（ａ）のＤｋの式において、[１−μｋＷｋηｋ]から[１−μｋＷｋηｋ]−１を求める計算も無視する。
【００５６】
図５及び図６より、本高速アルゴリズム全体を通して計算量が単位時間ステップ当りΟ（Ｎ）になっている。よって、高速Ｈ∞フィルタリングアルゴリズムは計算量はタップ数に比例することがわかる。また、この場合、高速Ｈ∞フィルタを１回実行するためにかかる計算量(乗算回数)は単位ステップ当り２８Ｎ＋１６であり、高速カルマンフィルタの１２Ｎ＋３と比べて約２倍の計算量(乗算回数)が必要である。
【００５７】
以上のように、変形Ｈ∞フィルタリングアルゴリズムではタップ数の２乗または３乗に比例する計算量が、本高速アルゴリズムではタップ数の１乗に比例するまで減少させることができる。
【００５８】
５．エコ−キャンセラ
つぎに、エコーキャンセラの例を取り上げ、本発明の効果を検証する。
まず、受信信号{ｕｋ}がエコーパスへの入力信号となることを考慮すれば、エコーパスの(時変)インパルス応答{ｈｉ[ｋ]}により、エコー{ｄｋ}の観測値{ｙｋ}は次式で表される。
【００５９】
【数２２】
0004067269.tif 000018
【００６０】
ここで、ｕｋ，ｙｋはそれぞれ時刻ｔｋ（＝ｋＴ；Ｔはサンプリング周期)における受信信号とエコーを表し、ｖｋは時刻ｔｋにおける平均値０の回線雑音とし、ｈｉ[ｋ]，ｉ＝０，・・・，Ｎ−１  は緩やかな変動を想定した時変インパルス応答であり、そのタップ数Ｎは既知とする。このとき、インパルス応答の推定値{ｈ＾ｉ[ｋ]}が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
【００６１】
【数２３】
0004067269.tif 000019
【００６２】
これをエコーから差し引けば（ｙｋ−ｄ＾ｋ≒０）、エコーをキャンセルすることができる。ただし、ｋ−ｉ＜０のときｕｋ−１＝０とする。
以上より、問題は直接観測可能な受信信号{ｕｋ}とエコー{ｙｋ}からエコーパスのインパルス応答{ｈｉ[ｋ]}を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにＨ∞フィルタを適用するには、まず式（２４）を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、問題がインパルス応答{ｈｉ[ｋ]}を推定することであるから、{ｈｉ[ｋ]}を状態変数ｘｋとし、ｗｋ程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
【００６３】
【数２４】
0004067269.tif 000020
【００６４】
このような状態空間モデルに対する変形および高速Ｈ∞フィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
以上より、インパルス応答の推定値{ｈ＾ｉ[ｋ]}が得られれば、これより疑似エコーが次のように逐次求めることができる。
【００６５】
【数２５】
0004067269.tif 000021
【００６６】
よって、これを実際のエコーから差し引き、エコーを打ち消せばエコーキャンセラが実現できる。このとき、推定誤差であるｅｋ＝ｙｋ−ｄ＾ｋは残留エコーと呼ばれる。
６．時不変インパルス応答に対する評価
【００６７】
（推定精度に対する評価）
エコーパスのインパルス応答が時間的に不変であり（ｈｉ[ｋ]＝ｈｉ）、かつそのタップ数Ｎが２４である場合について、シミュレーションを用いて変形Ｈ∞フィルタと高速Ｈ∞フィルタを評価する。
【００６８】
【数２６】
0004067269.tif 000022
【００６９】
なお、図７は、ここでのインパルス応答{ｈｉ}の値を示す図である。また、υｋは平均値０、分散σｖ２＝1.0×10-6の定常なガウス白色雑音とし、サンプリング周期Ｔを便宜上1.0とする。
また、受信信号{ｕｋ}は次のように２次のＡＲモデルで近似する。
ｕｋ＝α１ｕｋ−１＋α２ｕｋ−２＋ｗｋ′    (31)
ただし、α１＝0.7，α２＝0.1とし、ｗｋ′は平均値０、分散σｗ′２=0.04の定常なガウス白色雑音とする。
【００７０】
ここで、変形Ｈ∞フィルタと高速Ｈ∞フィルタとを比較する。
図８に変形Ｈ∞フィルタと高速Ｈ∞フィルタによるインパルス応答の推定結果の説明図を示す（初期値 :ｘ＾０｜０＝０，εｏ=20）。図８（ａ）、（ｂ）はγｆ＝105のときの両者の推定結果であり、図８(c)，(d)はγｆ＝2.0のときの推定結果（ｘ＾100｜100）である。これより、両者の推定精度に対する性能は等しいことがわかる。すなわち、高速化することで推定精度の低下は生じない。ここで、γｆが小さすぎるとフィルタの存在条件を満たさないので注意が必要である。また、γｆ＝1.0×105の場合は高速カルマンフィルタの結果とほぼ一致していた。以上より、高速Ｈ∞フィルタリングアルゴリズムは、高速カルマンフィルタリングアルゴリズムを含んでおり、かつγｆを調節することによって、収束を早くできることがわかる。
【００７１】
（計算時間の評価）
つぎに、エコーパスのインパルス応答が時間的に不変であり、かつタップ数を24,48,96,192,384と増加させたときの変形Ｈ∞フィルタと高速Ｈ∞フィルタの計算時間を評価する。ただし、１回の測定では結果にばらつきがあるので、４回測定した平均を結果として用いた。また、シミュレーションに用いるインパルス応答{ｈｉ}は図７の値とし、それ以降(２４≦ｋ＜Ｎ)のインパルス応答{ｈｉ}は０とする。ただし、フィルタの実行はステップ数１００までとする。計算時間は、ワークステーション(sparc,60MHz, 32MB)上のMATLABのコマンドetimeによって計測された。
【００７２】
図９に、計算時間の測定結果の図を示す。ここで、変形Ｈ∞フィルタ(2)はリカッチ方程式において、計算量がタップ数の２乗に比例するように行列計算を行ったものであり、変形Ｈ∞フィルタ(1)は計算量がタップ数の３乗に比例する行列計算を行ったものである(図３(ｂ)および図４参照)。このように、変形Ｈ∞フィルタは行列計算の順序によって計算量がタップ数の２乗または３乗に比例することになるが、いずれにせよ実用的ではない。
【００７３】
７．時変インパルス応答に対する評価
（追従性能の評価）
システム(インパルス応答)が時間的に変化したときの各アルゴリズムの追従性能についてエコーキャンセラの例を用いて評価する。ただし、インパルス応答のタップ数は48とし、{ｈｉ}は、図７の値を元に図１０(a)のように時間的に変化した場合を想定する。ただし、ｖｋは平均値０、σｖ２＝1.0×10−６の定常なガウス白色雑音とし、便宜上サンプリング周期をＴ＝１とする。また、受信信号{ｕｋ}は次のように２次のARモデルで近似する。
ｕｋ＝α１ｕｋ−１＋α２ｕｋ−２＋ｗｋ′    (32)
【００７４】
ここで、α１=0.7、α２=0.1とし、ｗｋ′は平均値0、分散σ２ｗ´=0.04の定常なガウス白色雑音とする。
図１０及び図１１に、各アルゴリズムのシミュレーション結果の図を示す。これは、高速Ｈ∞フィルタ（高速 HF）、高速カルマンフィルタ（高速 KF）およびＬＭＳアルゴリズム（ＬＭＳ）の時変システムの追従性能を示すものである。図１０(ｂ)はγｆ=2.0の場合の高速Ｈ∞フィルタによる推定値であり、図１１(ａ)は高速カルマンフィルタによる推定値である。ただし、高速Ｈ∞フィルタの初期値はｘ＾０｜０＝０，εｏ＝２０とし、高速カルマンフィルタの初期値も同様に設定した。また、図１１(ｂ)にＬＭＳアルゴリズムによる推定値を示す。ただし、初期値はｈ＾ｏ＝０、ステップサイズは安定かつ早い収束を与えるようにμ＝0.5とした。これより、高速Ｈ∞フィルタの追従性能が飛躍的に優れており、インパルス応答の変化後約３０ステップで推定値が安定していることがわかる。一方、高速カルマンフィルタとＬＭＳアルゴリズムについては全く追従できていない。
【００７５】
一般に、システム雑音を伴わないＨ∞フィルタの追従性が低下する原因は、Ｐ＾ｋ｜ｋ―１の対角成分の減少によりフィルタゲインの値が小さくなり推定値の更新量が減少するからである。つまり、ステップ数が経過するとほとんど推定値を更新しなくなる。よって、カルマンフィルタやＨ∞フィルタの追従性を向上させるためには、行列のＰ＾ｋ｜ｋ―１の対角成分の値に外部から適当な値を加えてやれば良い。しかし、直接導入したのでは観測行列Ｈｋのシフト特性を利用した高速アルゴリズムを導出できない。この問題を重みρ＝１−γｆ−２をＨ∞評価基準に導入することによって理論的に解決したことが本発明の大きな特徴のひとつである。この重みρは高速Ｈ∞フィルタリングアルゴリズムのＳｋの更新式の中に次のように現れる。
【００７６】
（高速Ｈ∞フィルタの補助変数Ｓｋの更新）
高速Ｈ∞フィルタの補助変数Ｓｋは、次式の通りである。
Ｓｋ＝ρＳｋ−１＋ｅＴｋＷｋｅ〜ｋ，０＜ρ＝１−γｆ−２≦１
また、高速Ｈ∞フィルタリングアルゴリズムにおいて、Ｓｋは、Ｋuｋの式でＳｋ−１の形で用いられる。よって、より大きな値の更新を行うためには、よりＳｋ−１が大きくなければならない。つまりＳｋを小さく保ち大きな値の更新を保持する必要がある。ρの存在はＳｋの急激な増大を防ぎ、結果的にシステム雑音を付加することと等価となり、追従性の向上につながっている。また、重みρはρ＝１−γｆ−２で定義されているので、シミュレーションで確認されたとおりγｆを変化させることで追従性を変化させることができる。
【００７７】
図１２に、γｆとρの関係図を示す。これによるとγｆ＝3.0のときρ＝0.8889なのでＳｋ−１の89％がＳｋに伝えられることになる。しかし、あまりγｆを小さくしすぎるとＳｋ−１の影響が著しく低下すると同時にフィルタの存在条件を満たさなくなるため、注意が必要である。また、γｆが大きい場合はρ≒１となりＳｋの増大を全く抑制しないため追従性が低下する。特にγf＝∞（ρ＝１）のとき本高速アルゴリズムは高速カルマンフィルタリングアルゴリズムと完全に一致する。
【００７８】
（計算時間の評価）
図１３に、高速Ｈ∞フィルタ、高速カルマンフィルタおよびＬＭＳアルゴリズムにおけるインパルス応答のタップ数(tap number)と計算時間[s]の関係図を示す。なお、フィルタの実行ステップ数:300, γｆ=3.0とした。そして、高速Ｈ∞フィルタ、高速カルマンフィルタリングアルゴリズムおよびＬＭＳアルゴリズムに対して、図１０及び図１１の例においてタップ数を48,96,192,384と増加させたときの計算時間を計測した。ただし、１回の測定では結果にばらつきがあるので、４回測定しその平均をとった。
【００７９】
いずれのアルゴリズムも計算量がタップ数の１乗に比例することが確認できる。また、タップ数が多い場合、高速Ｈ∞フィルタリングアルゴリズムの計算時間は高速カルマンフィルタリングアルゴリズムと比べて、約２倍弱であり、実用的なＬＭＳアルゴリズムと比較しても４倍程度であることがわかる。追従性を考慮すれば、高速Ｈ∞フィルタリングアルゴリズムの有効性は十分に高いと言えよう。
【００８０】
８．補題の証明
ここで、上述の補題の証明について説明する。
（補題１の証明）
Ｐｋの逆行列をとれば、式（３３）となる。さらに、逆行列の補助定理を用いれば、式（３４）に示すように、行列Ｐｋに関する再帰式が得られる。
【００８１】
【数２７】
0004067269.tif 000023
【００８２】
ここで、ＰｋをＰ＾ｋ＋１｜ｋと見倣せば、上式が式（１４）のリカッチ方程式を満たすことがわかる。
（補題２の証明）
ゲイン行列Ｋｋが次のように整理できる。
【００８３】
【数２８】
0004067269.tif 000024
【００８４】
さらに、Ｇｋ＝（ρ＋ＨｋＰｋ−１ＨＴｋ）/（１＋ＨｋＰｋ−１ＨＴｋ）と、ＨｋＫ〜ｋ＝ＨｋＰｋ−１ＨＴｋ/（１＋ＨｋＰｋ−１ＨＴｋ）を用いれば、ゲイン行列Ｋｋの第１ブロック列から式（１８）のようにフィルタゲインを得ることができる。
【００８５】
（補題３の証明）
ゲイン行列Ｋｉ，ｉ＝０，・・・，ｋ  が与えられたと仮定し、次のＫｋ＋１を求めることにする。
Ｑｋ＋１Ｋｋ＋１＝Ｃｋ＋１Ｔ    (36)
まず、Ｃｋのシフト特性を利用するため、新たに式（３７）と式（３８）を導入する。このとき、ＱＵｋは式（３９）のように再帰的に表され、かつ次の式（４０）のように分割される。
【００８６】
【数２９】
0004067269.tif 000025
【００８７】
この表記を用いれば、時間ステップｋおよびｋ＋１の式（３６）は、次式に含まれる。
【００８８】
【数３０】
0004067269.tif 000026
【００８９】
これらの表記に基づけば、Ｋｋを直接求める変わりに、次式を満たすＫＵｋ∈Ｒ（Ｎ＋１）×２を求める方が便利である。
【００９０】
【数３１】
0004067269.tif 000027
【００９１】
そのため、式（４１）から得られる式（４５）を用いれば、Ｋ∪ｋ∈Ｒ（Ｎ＋１）×２を式（４６）のように整理できる。
【００９２】
【数３２】
0004067269.tif 000028
【００９３】
ここで、Ｋ∪ｋは、ｍｋ∈ＲＮ×２とμｋ∈Ｒ１×２に分割される。また、αＴｋ−ｃＴｋ＝−（ｃｋＴ＋ＡｋＴＣｋＴ）に注意されたい。さらに、Ｑ∪ｋは逆行列が存在すると仮定し、補助変数Ａｋ∈ＲＮ×１とＳｋ∈Ｒは次式を満たす。
【００９４】
【数３３】
0004067269.tif 000029
【００９５】
ここで、上式の下ブロックはＴｋ＋ＱｋＡｋ＝０またはＴｋＴ＝−ＡｋＴＱｋＴを与える。
次に、ＣＴｋの上ブロックに影響することなく、式（４６）のμｋを消去するために、次の式（４８）のような補助変数Ｂｋ∈ＲＮ×１とＦｋ∈Ｒを導入し、さらに式（４６）のＫ∪ｋの分割からＢ∪ｋＦｋ−１μｋを引けば、式（４９）を得る。
【００９６】
【数３４】
0004067269.tif 000030
【００９７】
さらに、左から式（４９）の左辺にＱ∪ｋを掛ければ、次のように整理できる。
【００９８】
【数３５】
0004067269.tif 000031
【００９９】
上式の左辺に式（４９）を代入すれば、式（４３）は次式で表される。
【０１００】
【数３６】
0004067269.tif 000032
【０１０１】
これは式（４２）と同じ形であり、式（５１）の上ブロックから次式（５２）を導くことができる。
Ｑｋ＋１（ｍｋ−ＢｋＦｋ−１μｋ）＝Ｃｋ＋１Ｔ    (52)
ここで、式（３６）と式（５２）を比較すれば、ゲイン行列Ｋｋの更新式を得ることができる。
【０１０２】
（補題４）
補助変数ＡｋとＳｋは次のように得られる。
【０１０３】
【数３７】
0004067269.tif 000033
【０１０４】
ただし、Ａ−１＝０、Ｓ−１＝１/ε０である。
（証明）ＡｋとＳｋの式（５５）と式（３９）を用いると、式（５６）を得る。
【０１０５】
【数３８】
0004067269.tif 000034
【０１０６】
一方、式（４１）の両辺にＷｋ［Ｃｋ＋ＣｋＡｋ−１］を掛ければ、次式を得る。
【０１０７】
【数３９】
0004067269.tif 000035
【０１０８】
式（５６）から式（５７）を引けば、次式（５８）が成り立つ。
【０１０９】
【数４０】
0004067269.tif 000036
【０１１０】
これを式（４７）と比較すれば、αｋＴ＝ＴｋＴ  Ｋｋ＝−ＡｋＴ  ＣｋＴより、式（５３）と式（５４）を得る。
【０１１１】
（補題５）
補助変数Ｄｋ＝ＢｋＦｋ−１は、次式（５９）のように得られる。また、Ｆｋは、次式（６０）で更新される。
【０１１２】
【数４１】
0004067269.tif 000037
【０１１３】
ただし、ηｋ＝Ｃ∪ｋＤ∪ｋ−１＝ｃｋ―Ｎ＋Ｃｋ＋１Ｄｋ−１,  Ｄ−１＝０,Ｆ−１＝０  である。
（証明）ＢｋとＦｋを更新するため、式（６１）を用いれば、式（６２）が成り立つ。
【０１１４】
【数４２】
0004067269.tif 000038
【０１１５】
上式を式（６１）と同じ形に変形するため、式（６２）から  Ｃ∪ｋＴＷｋＣ∪ｋＢｋ−１∪  を引けば、次式を得る。
【０１１６】
【数４３】
0004067269.tif 000039
【０１１７】
この最後の式と式（４８）を比較すれば、Ｂ∪ｋに関する再帰式を得る。
【０１１８】
【数４４】
0004067269.tif 000040
【０１１９】
これより、ＢｋとＦｋが更新される。
しかし、それらはＤ∪ｋ＝Ｂ∪ｋＦｋ−１あるいはＤｋ＝ＢｋＦｋ−１∈ＲＮ×１としてだけ用いられるので、式（４８）と式（６４）を式（６５）を用いて書き換えた方が便利である。また、この行列Ｄｋは式（６６）を満たす。
【０１２０】
【数４５】
0004067269.tif 000041
【０１２１】
次に、式（６３）にＦｋ−１−１を掛ければ、式（６７）となり、さらにＤ∪ｋ−１＝Ｂ∪ｋ−１Ｆｋ−１−１を用いれば次式（６８）のように整理できる。
【０１２２】
【数４６】
0004067269.tif 000042
【０１２３】
これより、式（６８）に［１−μｋＷｋＣ∪ｋＤ∪ｋ−１］−１を掛ければ、次式が得られる。
【０１２４】
【数４７】
0004067269.tif 000043
【０１２５】
これを式（６６）と比較すれば、最終的にＤｋとＦｋの更新式が得られる。
（補題６（高速処理に適した存在条件）の証明）
上述したように、式（２２）、（２３）の存在条件を用いれば、計算量Ο（Ｎ）で高速Ｈ∞フィルタの存在が検査できる。その証明を以下に示す。
次式（６９）で示す２×２の行列Ｒｅ，ｋの特性方程式を解けば、Ｒｅ，ｋの固有値λｉが次式（７０）のように得られる。
【０１２６】
【数４８】
0004067269.tif 000044
【０１２７】
もし、次式（７１）が成り立てば、行列Ｒｅ，ｋの２つの固有値の１つは正となり、もう１つは負となり、行列ＲｋとＲｅ，ｋは同じイナーシャをもつ。これより、次式（７２）を用いれば、式（２２）の存在条件が得られる（即ち、式（７１））。ここで、ＨｋＫ〜ｋの計算がΟ（Ｎ）回の掛け算を必要としている。
【０１２８】
【数４９】
0004067269.tif 000045
【０１２９】
【発明の効果】
本発明によると、以上のように、新たなＨ∞評価基準に基づいて開発した変形Ｈ∞フィルタの高速アルゴリズム（高速Ｈ∞フィルタリングアルゴリズム）を用いて、時不変および時変システムの高速実時間同定および推定を実現することができる。また、本発明によると、本アルゴリズムの特殊な場合として高速カルマンクフィルタリングアルゴリズムを含み、また、時変システムの追従性を支配するシステム雑音の共分散に対応する項を理論的に決定することができる。また、本発明によると、突然回線が切り替わるような激しく変化する時変システムのエコーキャンセラなどのように、システム（インパルス応答）が時間的に不連続に変化する場合において、特に、非常に有効な高速時変システム同定方法を提供することができる。また、本発明によると、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム同定方法を提供することができる。
【図面の簡単な説明】
【図１】  高速アルゴリズムのフローチャート。
【図２】  変形Ｈ∞フィルタリングアルゴリズムの各部分の計算量の説明図（１）。
【図３】  変形Ｈ∞フィルタリングアルゴリズムの各部分の計算量の説明図（２）。
【図４】  リカッチ再帰法の計算量の説明図。
【図５】  高速Ｈ∞フィルタリングアルゴリズムの計算量の説明図（１）。
【図６】  高速Ｈ∞フィルタリングアルゴリズムの計算量の説明図（２）。
【図７】  インパルス応答{ｈｉ}の値を示す図。
【図８】  変形Ｈ∞フィルタと高速Ｈ∞フィルタによるインパルス応答の推定結果の比較説明図。
【図９】  計算時間の測定結果の図。
【図１０】  各アルゴリズムのシミュレーション結果の図（１）。
【図１１】  各アルゴリズムのシミュレーション結果の図（２）。
【図１２】  γｆとρの関係図。
【図１３】  高速Ｈ∞フィルタ、高速カルマンフィルタおよびＬＭＳアルゴリズムにおけるインパルス応答のタップ数(tap number)と計算時間[s]の関係図。
【図１４】  システム同定のための構成図。
【図１５】  インパルス応答の調節機構についての構成図。
【図１６】  通信系とエコーについての説明図。
【図１７】  エコーキャンセラの原理図。
【符号の説明】
１  未知システム
２  適応フィルタ
３  ＦＩＲディジタルフィルタ
４  適応アルゴリズム
【図１】
0004067269.tif 000057
【図２】
0004067269.tif 000058
【図３】
0004067269.tif 000059
【図４】
0004067269.tif 000060
【図５】
0004067269.tif 000061
【図６】
0004067269.tif 000062
【図７】
0004067269.tif 000063
【図８】
0004067269.tif 000064
【図９】
0004067269.tif 000065
【図１０】
0004067269.tif 000066
【図１１】
0004067269.tif 000067
【図１２】
0004067269.tif 000068
【図１３】
0004067269.tif 000069
【図１４】
0004067269.tif 000070
【図１５】
0004067269.tif 000071
【図１６】
0004067269.tif 000072
【図１７】
0004067269.tif 000073