JP4067269B2 - システム同定方法 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、システム同定方法に係り、特に、新たなＨ_∞評価基準に基づいて開発された変形Ｈ_∞フィルタの高速Ｈ_∞フィルタリングアルゴリズムを用いて時変システムを高速に実時間同定するものである。
【０００２】
【従来の技術】
一般に、システム同定(system identification)とは、入出力データに基づいてシステムの入出力関係の数理モデル(伝達関数、あるいはインパルス応答など)を推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
【０００３】
（基本原理）
図１４に、システム同定のための構成図を示す。このシステムは、未知システム１、適応フィルタ２を備える。また、適応フィルタ２は、ＦＩＲディジタルフィルタ３、適応アルゴリズム４を有する。
以下に、未知システム１を同定する出力誤差方式の一例を説明する。ここで、ｕ_ｋは未知システム１の入力、ｄ_ｋは所望信号であるシステムの出力、ｄ＾_ｋはフィルタの出力である。（なお、「＾」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。）
【０００４】
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差ｅ_ｋ＝ｄ_ｋ−ｄ＾_ｋを最小にするように適応アルゴリズムによってＦＩＲディジタルフィルタ３の係数を調節する。
【０００５】
図１５に、インパルス応答の調節機構についての構成図を示す。
ここで、適応アルゴリズムの一例として、その簡便さより、次のＬＭＳアルゴリズム(least mean square algorithm)が広く用いられている。
【０００６】
【数６】

【０００７】
また、一般に、時変システムの同定には、例えば、収束性が早いカルマンフィルタが適している。
【０００８】
【数７】

【０００９】
ここで、求めるべきインパルス応答{ｈ_ｉ}は状態推定値ｘ＾_ｋ｜ｋとして得られ、システムへの入力{ｕ_ｋ}は観測行列Ｈ_ｋの要素として用いられている。
また、σ^２ _ｗ＝０のときのカルマンフィルタに対して、観測行列Ｈ_ｋのシフト特性（Ｈ_ｋ＋１（ｉ＋１）＝Ｈ_ｋ（ｉ））を利用して、単位時間ステップ当たりの計算量をＮに比例した演算回数、すなわちＯ(Ｎ)までに軽減した高速カルマンフィルタリングアルゴリズムが知られている。詳細は、１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」など参照。
【００１０】
（エコーキャンセラへの適用例）
国際電話など長距離電話回線では,信号増幅などの理由から４線式回線が用いられている。一方、加入者回線は比較的短距離なので、２線式回線が使用されている。
図１６に通信系とエコーについての説明図を示す。２線式回線と４線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Ｂからの信号(音声)は話者Ａのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は４線式回線に漏れ、増幅された後、再び受信者(話者Ａ)に戻ると云った現象が起こる。これがエコー（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 に示されている。
【００１８】
１．記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
ｘ_ｋ： 状態ベクトルまたは単に状態 ； 未知、これが推定の対象となる。
ｘ_０： 初期状態 ； 未知である。
ｗ_ｋ： システム雑音 ； 未知である。
ｖ_ｋ： 観測雑音 ； 未知である。
ｙ_ｋ： 観測信号 ； フィルタの入力となり、既知である。
ｚ_ｋ： 出力信号 ； 未知である。
Ｈ_ｋ： 観測行列 ； 既知である。
Ｌ_ｋ： 出力行列 ； 既知である。
ｘ＾_ｋ｜ｋ： 観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの状態推定値； フィルタ方程式によって与えられる。
ｘ＾_０｜０：状態の初期推定値 ； 本来未知であるが、便宜上０が用いられる。
Ｋ_{ｓ，ｋ＋１}：フィルタゲイン ；行列Ｐ＾_{ｋ｜ｋ−１}から得られる。
Σ_ｗｋ ： システム雑音の共分散行列に対応 ； 理論上は既知であるが、実際には未知である。
Ｐ＾_{ｋ｜ｋ−１}： ｘ＾_{ｋ｜ｋ−１}の誤差の共分散行列に対応 ； リカッチ方程式によって与えられる。
Ｐ＾_１｜０： 初期状態の誤差の共分散行列に対応 ； 本来未知であるが、便宜上ε_０Ｉが用いられる。
σ² _ｖ ： 観測雑音の分散 ； 理論上は既知として扱われるが、実際には未知である。
σ² _ｗ ： システム雑音の分散 ； 理論上は既知として扱われるが、実際には未知である。
【００１９】
なお、記号の上に付される”＾”は、推定値の意味であり、"∪"は、行列を１行増やしたことを表す。”また、”〜”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、Ｌ、Ｈ、Ｐ、Ｋ等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
２．変形Ｈ_∞フィルタ
【００２０】
つぎに、次式（７）〜（９）のような状態空間モデルを考える。
【数８】

【００２１】
ここで、エコーキャンセラなどを想定し、Ｌ_ｋ＝Ｈ_ｋ（Ｈ_ｋ＝[ｕ_ｋ ｕ_ｋ−２・・・ ｕ_{ｋ−Ｎ＋１}]） とする。このような状態空間モデルに対して、次式（１０）のようなＨ_∞評価基準（新たに左辺にγ_ｆが入っている）を提案する。
【００２２】
【数９】

【００２３】
この評価基準を満たすレベルγ_ｆの変形Ｈ_∞フィルタは、ρやΣ_ｗｋがγ_ｆに依存しないと仮定すれば、システム同定の分野で通常知られる方法を適用することによって、次の式（１１）〜式(１４)で与えられる。なお、その通常知られる方法として、例えば、
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.
等を参照のこと。
【００２４】
【数１０】

【００２５】
また、評価基準の重みρは予めに決められた上限値γ_ｆに依存するので、上述のアルゴリズムは通常のＨ_∞フィルタとは本質的に異なる。本アルゴリズムはρで重みづけされた外乱(初期状態ｘ_０,システム雑音｛ｗ_ｉ｝,観測雑音｛υ_ｉ｝からフィルタ誤差｛ｅ_ｆ，ｉ｝への最大エネルギーゲインをγ_ｆ ^２より小さく押えているので、外乱に対してロバスト(頑健)なフィルタリングアルゴリズムとなる。この性質が時変システムの追従特性に反映される。また、γ_ｆ→∞のとき、ρ＝１、Σ_ｗｋ＝０となり、変形Ｈ_∞フィルタは通常のカルマンフィルタと一致する。
【００２６】
変形Ｈ_∞フィルタの主な計算上の負担は、Ｎ^２またはＮ^３に比例する計算量を必要とするＰ＾_{ｋ＋１｜ｋ}∈Ｒ^Ｎ×Ｎの更新の際に生じる。すなわち、単位時間ステップ当たりＯ（Ｎ^２）の算術演算が必要となる。ここで、タップ数Ｎは状態ベクトルｘ_ｋの次元と一致する。ゆえに、ｘ_ｋの次元が増加するにつれて変形Ｈ_∞フィルタの実行に要する計算時間は急速に増大する。この計算上の負担を克服するために変形Ｈ_∞フィルタの高速アルゴリズムの導出が必要となる。
【００２７】
３．高速Ｈ_∞フィルタリングアルゴリズム
変形Ｈ_∞フィルタの計算量は、式（１４）のリカッチ方程式（誤差の共分散方程式）の計算に支配される。よって、変形Ｈ_∞フィルタを高速に処理するためには、リカッチ方程式を用いずに式（１３）のフィルタゲインを直接決定できれば、大幅に計算量を削減できる。
しかし、フィルタゲインＫ_ｓ，ｋ∈Ｒ^Ｎ×１を求める高速アルゴリズムの導出は困難なため、次のように定義されるゲイン行列を高速に計算するアルゴリズムを導出することを考える。
【００２８】
【数１１】

【００２９】
このとき、次の補題が成り立つ。
補題1
行列Ｐ_ｋは式（１４）のリカッチ方程式を満たす。これより、ゲイン行列Ｋ_ｋが求まれば、次の補題からフィルタゲインＫ_ｓ，ｋが直ちに得られる。
【００３０】
補題2
変形Ｈ_∞フィルタのフィルタゲインはＫ_ｓ，ｋは、ゲイン行列Ｋ_ｋを用いて次のように得られる。実際、ゲイン行列Ｋ_ｋは次の再帰的方法によって高速に計算できる。
【００３１】
【数１２】

【００３２】
補題3
ゲイン行列Ｋ_ｋは次のように更新される。
【００３３】
【数１３】

【００３４】
ここで、ｍ_ｋ∈Ｒ^Ｎ×２とμ_ｋ∈Ｒ^１×２は行列Ｋ^∪ _ｋ＝Ｑ^∪ _ｋ ^−１Ｃ^∪ _ｋの次の分割によって得られる。
【００３５】
【数１４】

【００３６】
また、補助変数Ａ_ｋ∈Ｒ^Ｎ×１，Ｓ_ｋ∈ＲおよびＢ_ｋＦ^−１ _ｋ∈Ｒ^Ｎ×１も同様に得られる。
結論として、高速Ｈ_∞フィルタリングアルゴリズムは以下のように要約することができる。
図１に、高速アルゴリズムのフローチャートを示す。なお、Lは最大データ数を示す。
[ ステップ０] 再帰式の初期条件を以下のようにする。ここでε_０は充分に大きい正の定数である。
【００３７】
【数１５】

【００３８】
[ ステップ１] 時刻ｋと最大データ数Lとを比較する。時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。）
[ ステップ２] Ａ_ｋとＳ_ｋを以下のように再帰的に決定する。
【００３９】
【数１６】

【００４０】
[ ステップ３] Ｋ^∪ _ｋを以下のように計算する。
【００４１】
【数１７】

【００４２】
[ ステップ４] Ｋ^∪ _ｋを以下のように分割する。
【００４３】
【数１８】

【００４４】
[ ステップ５] Ｄ_ｋを決定し、Ｋ_ｋ＋１を通して、ゲイン行列Ｋ_{ｓ，ｋ＋１}を以下のように得る。
【００４５】
【数１９】

【００４６】
ここで、η_ｋ∈Ｒ^２×１，Ｄ_ｋ∈Ｒ^Ｎ×１，
Ｋ_ｋ＋１∈Ｒ^Ｎ×２，Ｋ_{ｓ，ｋ＋１}∈Ｒ^Ｎ×１，
０＜ρ＝１−γ_ｆ ^−２≦１，γ_ｆ＞１ である。
[ ステップ６] Ｈ_∞フィルタのフィルタ方程式を以下のように更新する。
【００４７】
【数２０】

【００４８】
[ ステップ７] 時刻ｋを進ませて（ｋ＝ｋ＋１）、 ステップ２に戻り、データがある限り続ける。
（高速処理に適した存在条件：補題６）
また、次の存在条件を用いれば、計算量Ο（Ｎ）で高速Ｈ_∞フィルタの存在性が検査できる。
【００４９】
【数２１】

【００５０】
４．本高速アルゴリズムの計算量
つぎに、高速Ｈ_∞フィルタリングアルゴリズムの計算量が、変形Ｈ_∞フィルタリングアルゴリズムの計算量と比べて如何に減少するかを考察する。ここで、式の計算量の評価は乗算回数のみに注目し、以下のような方法で計算した。
（Ｊ×Ｋ行列）×（Ｋ×Ｌ行列）の乗算回数＝Ｊ×Ｋ×Ｌ（回）
ただし、行列やベクトルが３つ以上乗算されるときは特に図に示さない限り左から計算されるものとする。
【００５１】
（変形Ｈ_∞フィルタリングアルゴリズムの計算量）
図２及び図３に、変形Ｈ_∞フィルタリングアルゴリズムの各部分の計算量の説明図を示す。ただし、Ｎはタップ数（状態ベクトルｘ _ｋ のディメンジョン）である。ここで、図３（ａ）においてＲ_ｅ，ｋからＲ_ｅ，ｋ ^−１を求めるための計算は無視する。同様に、図２（ａ）において（Ｈ_ｋ＋１Ｐ＾_{ｋ＋１ | ｋ}Ｈ^Ｔ _ｋ＋１＋１）の部分から（Ｈ_ｋ＋１Ｐ＾_{ｋ＋１ | ｋ}Ｈ^Ｔ _ｋ＋１＋１）^−１を求める計算も無視する。
【００５２】
図２（ａ）、図３（ａ）および図３（ｂ）より、Ｋ_{ｓ | ｋ＋１}、Ｒ_ｅ，ｋおよびＰ＾_{ｋ＋１ | ｋ}は計算量がタップ数の２乗に比例することがわかる。よって、変形Ｈ_∞フィルタリングアルゴリズム全体の計算量は単位時間ステップ当りΟ（Ｎ^２）となる。
【００５３】
また、図４に、行列計算の順番を変えた場合の計算量の説明図を示す。すなわち、リカッチ方程式において、右辺第２項の部分の行列計算の順番を変えた場合の計算量を図４に示す。
上述の部分の計算量がインパルス応答の次数の３乗に比例するので、Ｐ＾_{ｋ＋１ | ｋ}の計算量もインパルス応答の次数の３乗に比例する。これに伴い、Ｈ_∞フィルタ全体の計算量もタップ数の２乗から３乗にまで増加する。
【００５４】
しかし、いずれのアルゴリズムにせよタップ数の２乗または３乗に比例する計算量をもつので、タップ数の増加と共にフィルタの実行にかかる計算上の負担が著しく増加する。実際、通信工学の分野に利用する場合などはタップ数が例えば４００程度であるため、このアルゴリズムでの実用はかなり困難なものとなる。
【００５５】
（高速Ｈ_∞フィルタリングアルゴリズムの計算量）
次に、図５及び図６に、高速Ｈ_∞フィルタリングアルゴリズムの計算量の説明図を示す。ここで、図５（ｂ）のＫ^ｕ _ｋの式においてＳ_ｋ ^−１はＳ_ｋから求められるが、その計算は無視する。同様に、図６（ａ）のＤ_ｋの式において、[１−μ_ｋＷ_ｋη_ｋ]から[１−μ_ｋＷ_ｋη_ｋ]^−１を求める計算も無視する。
【００５６】
図５及び図６より、本高速アルゴリズム全体を通して計算量が単位時間ステップ当りΟ（Ｎ）になっている。よって、高速Ｈ_∞フィルタリングアルゴリズムは計算量はタップ数に比例することがわかる。また、この場合、高速Ｈ_∞フィルタを１回実行するためにかかる計算量(乗算回数)は単位ステップ当り２８Ｎ＋１６であり、高速カルマンフィルタの１２Ｎ＋３と比べて約２倍の計算量(乗算回数)が必要である。
【００５７】
以上のように、変形Ｈ_∞フィルタリングアルゴリズムではタップ数の２乗または３乗に比例する計算量が、本高速アルゴリズムではタップ数の１乗に比例するまで減少させることができる。
【００５８】
５．エコ−キャンセラ
つぎに、エコーキャンセラの例を取り上げ、本発明の効果を検証する。
まず、受信信号{ｕ_ｋ}がエコーパスへの入力信号となることを考慮すれば、エコーパスの(時変)インパルス応答{ｈ_ｉ[ｋ]}により、エコー{ｄ_ｋ}の観測値{ｙ_ｋ}は次式で表される。
【００５９】
【数２２】

【００６０】
ここで、ｕ_ｋ，ｙ_ｋはそれぞれ時刻ｔ_ｋ（＝ｋＴ；Ｔはサンプリング周期)における受信信号とエコーを表し、ｖ_ｋは時刻ｔ_ｋにおける平均値０の回線雑音とし、ｈ_ｉ[ｋ]，ｉ＝０，・・・，Ｎ−１ は緩やかな変動を想定した時変インパルス応答であり、そのタップ数Ｎは既知とする。このとき、インパルス応答の推定値{ｈ＾_ｉ[ｋ]}が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
【００６１】
【数２３】

【００６２】
これをエコーから差し引けば（ｙ_ｋ−ｄ＾_ｋ≒０）、エコーをキャンセルすることができる。ただし、ｋ−ｉ＜０のときｕ_ｋ−１＝０とする。
以上より、問題は直接観測可能な受信信号{ｕ_ｋ}とエコー{ｙ_ｋ}からエコーパスのインパルス応答{ｈ_ｉ[ｋ]}を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにＨ_∞フィルタを適用するには、まず式（２４）を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、問題がインパルス応答{ｈ_ｉ[ｋ]}を推定することであるから、{ｈ_ｉ[ｋ]}を状態変数ｘ_ｋとし、ｗ_ｋ程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
【００６３】
【数２４】

【００６４】
このような状態空間モデルに対する変形および高速Ｈ_∞フィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
以上より、インパルス応答の推定値{ｈ＾_ｉ[ｋ]}が得られれば、これより疑似エコーが次のように逐次求めることができる。
【００６５】
【数２５】

【００６６】
よって、これを実際のエコーから差し引き、エコーを打ち消せばエコーキャンセラが実現できる。このとき、推定誤差であるｅ_ｋ＝ｙ_ｋ−ｄ＾_ｋは残留エコーと呼ばれる。
６．時不変インパルス応答に対する評価
【００６７】
（推定精度に対する評価）
エコーパスのインパルス応答が時間的に不変であり（ｈ_ｉ[ｋ]＝ｈ_ｉ）、かつそのタップ数Ｎが２４である場合について、シミュレーションを用いて変形Ｈ_∞フィルタと高速Ｈ_∞フィルタを評価する。
【００６８】
【数２６】

【００６９】
なお、図７は、ここでのインパルス応答{ｈ_ｉ}の値を示す図である。また、υ_ｋは平均値０、分散σ_ｖ ^２＝1.0×10^-6の定常なガウス白色雑音とし、サンプリング周期Ｔを便宜上1.0とする。
また、受信信号{ｕ_ｋ}は次のように２次のＡＲモデルで近似する。
ｕ_ｋ＝α_１ｕ_ｋ−１＋α_２ｕ_ｋ−２＋ｗ_ｋ′ (31)
ただし、α_１＝0.7，α_２＝0.1 とし、ｗ_ｋ′は平均値０、分散σ_ｗ′ ^２=0.04の定常なガウス白色雑音とする。
【００７０】
ここで、変形Ｈ_∞フィルタと高速Ｈ_∞フィルタとを比較する。
図８に変形Ｈ_∞フィルタと高速Ｈ_∞フィルタによるインパルス応答の推定結果の説明図を示す（初期値 :ｘ＾_０｜０＝０，ε _ｏ=20）。図８（ａ）、（ｂ）はγ_ｆ＝10⁵のときの両者の推定結果であり、図８(c)，(d)はγ_ｆ＝2.0のときの推定結果（ｘ＾ _{100 ｜ 100} ）である。これより、両者の推定精度に対する性能は等しいことがわかる。すなわち、高速化することで推定精度の低下は生じない。ここで、γ_ｆが小さすぎるとフィルタの存在条件を満たさないので注意が必要である。また、γ_ｆ＝1.0×10⁵の場合は高速カルマンフィルタの結果とほぼ一致していた。以上より、高速Ｈ_∞フィルタリングアルゴリズムは、高速カルマンフィルタリングアルゴリズムを含んでおり、かつγ_ｆを調節することによって、収束を早くできることがわかる。
【００７１】
（計算時間の評価）
つぎに、エコーパスのインパルス応答が時間的に不変であり、かつタップ数を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と増加させたときの計算時間を計測した。ただし、１回の測定では結果にばらつきがあるので、４回測定しその平均をとった。
【００７９】
いずれのアルゴリズムも計算量がタップ数の１乗に比例することが確認できる。また、タップ数が多い場合、高速Ｈ_∞フィルタリングアルゴリズムの計算時間は高速カルマンフィルタリングアルゴリズムと比べて、約２倍弱であり、実用的なＬＭＳアルゴリズムと比較しても４倍程度であることがわかる。追従性を考慮すれば、高速Ｈ_∞フィルタリングアルゴリズムの有効性は十分に高いと言えよう。
【００８０】
８．補題の証明
ここで、上述の補題の証明について説明する。
（補題１の証明）
Ｐ_ｋの逆行列をとれば、式（３３）となる。さらに、逆行列の補助定理を用いれば、式（３４）に示すように、行列Ｐ_ｋに関する再帰式が得られる。
【００８１】
【数２７】

【００８２】
ここで、Ｐ_ｋをＰ＾_{ｋ＋１｜ｋ}と見倣せば、上式が式（１４）のリカッチ方程式を満たすことがわかる。
（補題２の証明）
ゲイン行列Ｋ_ｋが次のように整理できる。
【００８３】
【数２８】

【００８４】
さらに、Ｇ_ｋ＝（ρ＋Ｈ_ｋＰ_ｋ−１Ｈ^Ｔ _ｋ）/（１＋Ｈ_ｋＰ_ｋ−１Ｈ^Ｔ _ｋ）と、Ｈ_ｋＫ^〜 _ｋ＝Ｈ_ｋＰ_ｋ−１Ｈ^Ｔ _ｋ/（１＋Ｈ_ｋＰ_ｋ−１Ｈ^Ｔ _ｋ）を用いれば、ゲイン行列Ｋ_ｋの第１ブロック列から式（１８）のようにフィルタゲインを得ることができる。
【００８５】
（補題３の証明）
ゲイン行列Ｋ_ｉ，ｉ＝０，・・・，ｋ が与えられたと仮定し、次のＫ_ｋ＋１を求めることにする。
Ｑ_ｋ＋１Ｋ_ｋ＋１＝Ｃ_ｋ＋１ ^Ｔ (36)
まず、Ｃ_ｋのシフト特性を利用するため、新たに式（３７）と式（３８）を導入する。このとき、Ｑ^Ｕ _ｋは式（３９）のように再帰的に表され、かつ次の式（４０）のように分割される。
【００８６】
【数２９】

【００８７】
この表記を用いれば、時間ステップｋおよびｋ＋１の式（３６）は、次式に含まれる。
【００８８】
【数３０】

【００８９】
これらの表記に基づけば、Ｋ_ｋを直接求める変わりに、次式を満たすＫ^Ｕ _ｋ∈Ｒ^{（Ｎ＋１）×２}を求める方が便利である。
【００９０】
【数３１】

【００９１】
そのため、式（４１）から得られる式（４５）を用いれば、Ｋ^∪ _ｋ∈Ｒ^{（Ｎ＋１）×２}を式（４６）のように整理できる。
【００９２】
【数３２】

【００９３】
ここで、Ｋ^∪ _ｋは、ｍ_ｋ∈Ｒ^Ｎ×２とμ_ｋ∈Ｒ^１×２に分割される。また、α^Ｔ _ｋ−ｃ^Ｔ _ｋ＝−（ｃ_ｋ ^Ｔ＋Ａ_ｋ ^ＴＣ_ｋ ^Ｔ）に注意されたい。さらに、Ｑ^∪ _ｋは逆行列が存在すると仮定し、補助変数Ａ_ｋ∈Ｒ^Ｎ×１とＳ_ｋ∈Ｒは次式を満たす。
【００９４】
【数３３】

【００９５】
ここで、上式の下ブロックはＴ_ｋ＋Ｑ_ｋＡ_ｋ＝０またはＴ_ｋ ^Ｔ＝−Ａ_ｋ ^ＴＱ_ｋ ^Ｔを与える。
次に、Ｃ^Ｔ _ｋの上ブロックに影響することなく、式（４６）のμ_ｋを消去するために、次の式（４８）のような補助変数Ｂ_ｋ∈Ｒ^Ｎ×１とＦ_ｋ∈Ｒを導入し、さらに式（４６）のＫ^∪ _ｋの分割からＢ^∪ _ｋＦ_ｋ ^−１μ_ｋを引けば、式（４９）を得る。
【００９６】
【数３４】

【００９７】
さらに、左から式（４９）の左辺にＱ^∪ _ｋを掛ければ、次のように整理できる。
【００９８】
【数３５】

【００９９】
上式の左辺に式（４９）を代入すれば、式（４３）は次式で表される。
【０１００】
【数３６】

【０１０１】
これは式（４２）と同じ形であり、式（５１）の上ブロックから次式（５２）を導くことができる。
Ｑ_ｋ＋１（ｍ_ｋ−Ｂ_ｋＦ_ｋ ^−１μ_ｋ）＝Ｃ_ｋ＋１ ^Ｔ (52)
ここで、式（３６）と式（５２）を比較すれば、ゲイン行列Ｋ_ｋの更新式を得ることができる。
【０１０２】
（補題４）
補助変数Ａ_ｋとＳ_ｋは次のように得られる。
【０１０３】
【数３７】

【０１０４】
ただし、Ａ_−１＝０、Ｓ_−１＝１/ε_０である。
（証明）Ａ_ｋとＳ_ｋ の式（５５）と式（３９）を用いると、式（５６）を得る。
【０１０５】
【数３８】

【０１０６】
一方、式（４１）の両辺にＷ_ｋ［Ｃ_ｋ＋Ｃ_ｋＡ_ｋ−１］を掛ければ、次式を得る。
【０１０７】
【数３９】

【０１０８】
式（５６）から式（５７）を引けば、次式（５８）が成り立つ。
【０１０９】
【数４０】

【０１１０】
これを式（４７）と比較すれば、α_ｋ ^Ｔ＝Ｔ_ｋ ^Ｔ Ｋ_ｋ＝−Ａ_ｋ ^Ｔ Ｃ_ｋ ^Ｔより、式（５３）と式（５４）を得る。
【０１１１】
（補題５）
補助変数Ｄ_ｋ＝Ｂ_ｋＦ_ｋ ^−１は、次式（５９）のように得られる。また、Ｆ_ｋは、次式（６０）で更新される。
【０１１２】
【数４１】

【０１１３】
ただし、η_ｋ＝Ｃ^∪ _ｋＤ^∪ _ｋ−１＝ｃ_ｋ―Ｎ＋Ｃ_ｋ＋１Ｄ_ｋ−１, Ｄ_−１＝０,Ｆ_−１＝０ である。
（証明）Ｂ_ｋとＦ_ｋを更新するため、式（６１）を用いれば、式（６２）が成り立つ。
【０１１４】
【数４２】

【０１１５】
上式を式（６１）と同じ形に変形するため、式（６２）から Ｃ^∪ _ｋ ^ＴＷ_ｋＣ^∪ _ｋＢ_ｋ−１ ^∪ を引けば、次式を得る。
【０１１６】
【数４３】

【０１１７】
この最後の式と式（４８）を比較すれば、Ｂ^∪ _ｋに関する再帰式を得る。
【０１１８】
【数４４】

【０１１９】
これより、Ｂ_ｋとＦ_ｋが更新される。
しかし、それらはＤ ^∪ _ｋ ＝Ｂ ^∪ _ｋ Ｆ _ｋ ^−１ あるいはＤ_ｋ＝Ｂ_ｋＦ_ｋ ^−１∈Ｒ^Ｎ×１としてだけ用いられるので、式（４８）と式（６４）を式（６５）を用いて書き換えた方が便利である。また、この行列Ｄ_ｋ は式（６６）を満たす。
【０１２０】
【数４５】

【０１２１】
次に、式（６３）にＦ_ｋ−１ ^−１を掛ければ、式（６７）となり、さらにＤ^∪ _ｋ−１＝Ｂ^∪ _ｋ−１Ｆ_ｋ−１ ^−１を用いれば次式（６８）のように整理できる。
【０１２２】
【数４６】

【０１２３】
これより、式（６８）に［１−μ_ｋＷ_ｋＣ^∪ _ｋＤ^∪ _ｋ−１］^−１を掛ければ、次式が得られる。
【０１２４】
【数４７】

【０１２５】
これを式（６６）と比較すれば、最終的にＤ_ｋとＦ_ｋの更新式が得られる。
（補題６（高速処理に適した存在条件）の証明）
上述したように、式（２２）、（２３）の存在条件を用いれば、計算量Ο（Ｎ）で高速Ｈ_∞フィルタの存在が検査できる。その証明を以下に示す。
次式（６９）で示す２×２の行列Ｒ_ｅ，ｋの特性方程式を解けば、Ｒ_ｅ，ｋの固有値λ_ｉが次式（７０）のように得られる。
【０１２６】
【数４８】

【０１２７】
もし、次式（７１）が成り立てば、行列Ｒ_ｅ，ｋの２つの固有値の１つは正となり、もう１つは負となり、行列Ｒ_ｋとＲ_ｅ，ｋは同じイナーシャをもつ。これより、次式（７２）を用いれば、式（２２）の存在条件が得られる（即ち、式（７１））。ここで、Ｈ_ｋＫ^〜 _ｋの計算がΟ（Ｎ）回の掛け算を必要としている。
【０１２８】
【数４９】

【０１２９】
【発明の効果】
本発明によると、以上のように、新たなＨ_∞評価基準に基づいて開発した変形Ｈ_∞フィルタの高速アルゴリズム（高速Ｈ_∞フィルタリングアルゴリズム）を用いて、時不変および時変システムの高速実時間同定および推定を実現することができる。また、本発明によると、本アルゴリズムの特殊な場合として高速カルマンクフィルタリングアルゴリズムを含み、また、時変システムの追従性を支配するシステム雑音の共分散に対応する項を理論的に決定することができる。また、本発明によると、突然回線が切り替わるような激しく変化する時変システムのエコーキャンセラなどのように、システム（インパルス応答）が時間的に不連続に変化する場合において、特に、非常に有効な高速時変システム同定方法を提供することができる。また、本発明によると、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム同定方法を提供することができる。
【図面の簡単な説明】
【図１】 高速アルゴリズムのフローチャート。
【図２】 変形Ｈ_∞フィルタリングアルゴリズムの各部分の計算量の説明図（１）。
【図３】 変形Ｈ_∞フィルタリングアルゴリズムの各部分の計算量の説明図（２）。
【図４】 リカッチ再帰法の計算量の説明図。
【図５】 高速Ｈ_∞フィルタリングアルゴリズムの計算量の説明図（１）。
【図６】 高速Ｈ_∞フィルタリングアルゴリズムの計算量の説明図（２）。
【図７】 インパルス応答{ｈ_ｉ}の値を示す図。
【図８】 変形Ｈ_∞フィルタと高速Ｈ_∞フィルタによるインパルス応答の推定結果の比較説明図。
【図９】 計算時間の測定結果の図。
【図１０】 各アルゴリズムのシミュレーション結果の図（１）。
【図１１】 各アルゴリズムのシミュレーション結果の図（２）。
【図１２】 γ_ｆとρの関係図。
【図１３】 高速Ｈ_∞フィルタ、高速カルマンフィルタおよびＬＭＳアルゴリズムにおけるインパルス応答のタップ数(tap number)と計算時間[s]の関係図。
【図１４】 システム同定のための構成図。
【図１５】 インパルス応答の調節機構についての構成図。
【図１６】 通信系とエコーについての説明図。
【図１７】 エコーキャンセラの原理図。
【符号の説明】
１ 未知システム
２ 適応フィルタ
３ ＦＩＲディジタルフィルタ
４ 適応アルゴリズム

Claims (5)

1. 時不変又は時変システムの高速実時間同定を行うシステム同定方法において、次式で表されるＨ_∞フィルタ方程式を用い、
   

   

   Ｈ_∞評価基準として、｛ （フィルタ誤差に対応する項／評価関数の重み（ρ）） ／ ［（初期状態の誤差に対応する項）＋（システム雑音に対応する項）＋（観測雑音に対応する項／評価関数の重み（ρ））］ ｝の最大値が、予め与えられた上限値（γ_ｆ ^２）より小さく抑えるように定めることにより、外乱に対して頑健なフィルタリングアルゴリズムとしたシステム同定方法。
   ただし、
   ｘ _ｋ ： 状態ベクトルまたは単に状態 ； 未知、これが推定の対象となる。
   ｘ _０ ： 初期状態 ； 未知である。
   ｗ _ｋ ： システム雑音 ； 未知である。
   ｖ _ｋ ： 観測雑音 ； 未知である。
   ｙ _ｋ ： 観測信号 ； フィルタの入力となり、既知である。
   ｚ _ｋ ： 出力信号 ； 未知である。
   Ｈ _ｋ ： 観測行列 ； 既知である。
   ｘ＾ _ｋ｜ｋ ： 観測信号ｙ _０ 〜ｙ _ｋ までを用いた時刻ｋの状態ｘ _ｋ の状態推定値 ； フィルタ方程式（１２）によって与えられる。
   ｘ＾ _０｜０ ：状態の初期推定値 ； 本来未知であるが、便宜上０が用いられる。
   Ｋ _{ｓ，ｋ＋１} ：フィルタゲイン ；行列Ｐ＾ _{ｋ＋１｜ｋ} から得られる。
   Σ _ｗｋ ： システム雑音の共分散行列に対応 ； 理論上は既知であるが、フィルタの実行前には未知である。
   Ｐ＾ _{ｋ｜ｋ−１} ： ｘ＾ _{ｋ｜ｋ−１} の誤差の共分散行列に対応 ； リカッチ方程式によって与えられる。
   Ｐ＾ _１｜０ ： 初期状態の誤差の共分散行列に対応 ； 本来未知であるが、便宜上ε _０ Ｉが用いられる。
   さらに、
   「フィルタ誤差」：
   （１１）式のｚ ^ｖ _ｋ｜ｋ と（９）式のｚ _ｋ ＝Ｈ _ｋ ｘ _ｋ の差（（１５）式のｅ _ｆ，ｉ に対応）の重み付きノルムの２乗の和。
   「初期状態の誤差」：
   時刻ｋ＝０のときの、以下に示す状態空間モデルの、（７）式の状態ｘ _k の初期値ｘ _０ とその推定値の差の重み付きノルム。
   「システム雑音に対応する項」：
   以下に示す状態空間モデルの、（７）式のシステム雑音ｗ _ｋ の重み付きノルムの２乗の和。
   「観測雑音に対応する項」：
   以下に示す状態空間モデルの、（８）式の観測雑音ｖ _ｋ の重み付きノルムの２乗の和。
   状態空間モデル
   

   
2. さらに、次式により時刻ｋの状態推定値ｘ＾_ｋ｜ｋから出力信号ｚ _ｋ を求めるようにした請求項１に記載のシステム同定方法。
   

   
3. ゲイン行列Ｋ _ｋ 、補助変数Ａ _ｋ 、Ｓ _ｋ 、Ｄ _ｋ 、状態推定値ｘ＾ _ｋ｜ｋ の再帰式の初期条件をそれぞれ以下のように定めるステップと、
   

   

   時刻ｋにおける補助変数Ａ _ｋ とＳ _ｋ を以下のように再帰的に決定するステップと、
   

   

   ゲイン行列Ｋ _ｋ に補助変数を含む行を増やした第２のゲイン行列Ｋ ^∪ _ｋ を以下のように計算するステップと、
   

   

   第２のゲイン行列Ｋ ^∪ _ｋ を以下のように分割し、第１及び第２の分割ゲイン行列を求めるステップと、
   

   

   第１及び第２の分割ゲイン行列を含む次式により、Ｄ _ｋ を決定し、時刻ｋ＋１におけるゲイン行列Ｋ _ｋ＋１ を求め、時刻ｋ＋１におけるフィルタゲインＫ _{ｓ，ｋ＋１} を求めるステップと、
   

   

   ここで、η _ｋ ∈Ｒ ^２×１ ，Ｄ _ｋ ∈Ｒ ^Ｎ×１ ，
   Ｋ _ｋ＋１ ∈Ｒ ^Ｎ×２ ，Ｋ _{ｓ，ｋ＋１} ∈Ｒ ^Ｎ×１ ，
   ０＜ρ＝１−γ _ｆ ^−２ ≦１，γ _ｆ ＞１
   上記で求められたフィルタゲインＫ _{ｓ，ｋ＋１} に基づき、前記フィルタ方程式（１２）を以下のように更新するステップと、
   

   

   時刻を進ませて、各前記ステップを繰り返すためのステップと
   を含む請求項２に記載のシステム同定方法。
4. さらに、高速処理に適した存在条件として、次式を用いることにより、計算量Ο（Ｎ）で前記高速Ｈ_∞フィルタの存在性を検査することを特徴とする請求項１乃至３のいずれかに記載のシステム同定方法。
   

   
5. 前記Ｈ_∞フィルタ方程式を適用し、状態推定値ｘ＾_ｋ｜ｋを求め、擬似エコーを次式のように推定し、求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現することを特徴とする請求項１又は３又は４に記載のシステム同定方法。
   

   

Images