JPWO2005015737A1 - システム推定方法及びプログラム及び記録媒体、システム推定装置 - Google Patents

Abstract

忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。まず、処理部は、記憶部又は入力部ら上限値γｆを読み出し又は入力する（Ｓ１０１）。処理部は、式（１５）によって忘却係数ρを決定する（Ｓ１０３）。その後、処理部は、忘却係数ρに基づき、式（１０）〜式（１３）のハイパーＨ∞フィルタを実行する（Ｓ１０５）。処理部１０１は、式（１７）（あるいは、後述の式（１８））の存在条件を計算し（Ｓ１０７）、その存在条件がすべての時刻で満たされれば（Ｓ１０９）、γｆをΔγだけ小さくして同じ処理を繰り返す（Ｓ１１１）。一方、あるγｆで存在条件を満たさなくなったとき（Ｓ１０９）、そのγｆにΔγを加えたものをγｆの最適値γｆｏｐとして出力部に出力及び／又は記憶部に記憶する（Ｓ１１３）。

Description

本発明は、システム推定方法及びプログラム及び記録媒体、システム推定装置に係り、特に、Ｈ_∞評価基準に基づいて開発されたハイパーＨ_∞フィルタの高速Ｈ_∞フィルタリングアルゴリズムを用いて、状態推定のロバスト化と忘却係数の最適化を同時に実現するシステム推定方法及びプログラム及び記録媒体、システム推定装置に関する。

一般に、システム推定とは、入出力データに基づいてシステムの入出力関係の数理モデル（伝達関数、あるいはインパルス応答など）のパラメータを推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、非特許文献１．１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
（基本原理）
図８に、システム推定のための構成図の例を示す（未知システムはＩＩＲ（Ｉｎｆｉｎｉｔｅ Ｉｍｐｕｌｓｅ Ｒｅｓｐｏｎｓｅ）フィルタで表現してもよい）。
このシステムは、未知システム１、適応フィルタ２を備える。また、適応フィルタ２は、ＦＩＲディジタルフィルタ３、適応アルゴリズム４を有する。
以下に、未知システム１を同定する出力誤差方式の一例を説明する。ここで、ｕ_ｋは未知システム１の入力、ｄ_ｋは所望信号であるシステムの出力、ｄ＾_ｋはフィルタの出力である。（なお、「＾」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。）
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差ｅ_ｋ＝ｄ_ｋ−ｄ＾_ｋを最小にするように適応アルゴリズムによってＦＩＲディジタルフィルタ３の係数を調節する。
また、従来、システムのパラメータ（状態）の推定には、誤差共分散行列の更新式（リカッチ方程式）に基づくカルマンフィルタが広く用いられて来た。詳細は、非特許文献２．Ｓ．Ｈａｙｋｉｎ：Ａｄａｐｔｉｖｅ ｆｉｌｔｅｒ ｔｈｅｏｒｙ，Ｐｒｅｎｔｉｃｅ−Ｈａｌｌ（１９９６）などに示されている。
以下にカルマンフィルタの基本原理について説明する。
次式のように、状態空間モデルで表された線形システム

の状態ｘ_ｋの最小分散推定値ｘ＾_ｋ｜ｋは、状態の誤差共分散行列Σ＾_{ｋ｜ｋ−１}を用いて次のように得られる。

ただし、

ｘ_ｋ：状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。
ｙ_ｋ：観測信号；フィルタの入力となり、既知である。
Ｈ_ｋ：観測行列；既知である。
ｖ_ｋ：観測雑音；未知である。
ρ：忘却係数；一般に試行錯誤で決定される。
Ｋ_ｋ：フィルタゲイン；行列Σ＾_{ｋ｜ｋ−１}から得られる。
Σ＾_ｋ｜ｋ：ｘ＾_ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾_{ｋ＋１｜ｋ}：ｘ＾_{ｋ＋１｜ｋ}の誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾_１｜０：初期状態の共分散行列に対応；本来未知であるが、便宜上ε_０Ｉが用いられる。
また、本発明者は、既に高速Ｈ_∞フィルタによるシステム同定アルゴリズムを提案した（特許文献１参照）。これは、システム同定のために新たにＨ_∞評価基準を定め、これに基づくハイパーＨ_∞フィルタの高速アルゴリズムを開発すると共に、この高速Ｈ_∞フィルタリングアルゴリズムに基づく高速時変システム同定方法である。高速Ｈ_∞フィルタリングアルゴリズムは、単位時間ステップ当たり計算量Ｏ（Ｎ）で急激に変化する時変システムの追跡が可能である。上限値の極限で高速カルマンフィルタリングアルゴリズムと完全に一致する。このようなシステム同定により時不変および時変システムの高速実時間同定および推定を実現することができる。
なお、システム推定の分野で通常知られる方法として、例えば、非特許文献２、３参照のこと。
（エコーキャンセラへの適用例）
国際電話など長距離電話回線では，信号増幅などの理由から４線式回線が用いられている。一方、加入者回線は比較的短距離なので、２線式回線が使用されている。
図９に通信系とエコーについての説明図を示す。２線式回線と４線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Ｂからの信号（音声）は話者Ａのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は４線式回線に漏れ、増幅された後、再び受信者（話者Ａ）に戻ると云った現象が起こる。これがエコー（ｅｃｈｏ）である。エコーは、伝送距離が長くなるにつれて（遅延時間が長くなるにつれて）影響が大きくなり、著しく通話の品質を劣化させる（パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する）。
図１０に、エコーキャンセラの原理図を示す。
そこで、図示のようにエコーキャンセラ（ｅｃｈｏ ｃａｎｃｅｌｌｅｒ）を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
エコーパスのインパルス応答の推定は、残留エコーｅ_ｋの平均２乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Ａがらの信号（音声）である。一般に、話者２人が同時に話し始めた（ダブルトーク）ときはインパルス応答の推定を中断する。また、ハイブリッドトランスのインパルス応答長は５０［ｍｓ］程度なので、サンプリング周期を１２５［μｓ］とするとエコーパスのインパルス応答の次数は実際は４００程度となる。
１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」 Ｓ．Ｈａｙｋｉｎ：Ａｄａｐｔｉｖｅ ｆｉｌｔｅｒ ｔｈｅｏｒｙ，Ｐｒｅｎｔｉｃｅ−Ｈａｌｌ（１９９６） Ｂ．Ｈａｓｓｉｂｉ，Ａ．Ｈ．Ｓａｙｅｄ，ａｎｄ Ｔ．Ｋａｉｌａｔｈ： ″Ｉｎｄｅｆｉｎｉｔｅ−Ｑｕａｄｒａｔｉｃ Ｅｓｔｉｍａｔｉｏｎ ａｎｄ Ｃｏｎｔｒｏｌ″，ＳＩＡＭ（１９９６） 特開２００２−１３５１７１号公報

しかしながら、式（１）〜（５）のような従来の忘却係数ρを入れたカルマンフィルタでは、忘却係数ρの値は試行錯誤で決定しなければならず非常に手間が掛かった。さらに、決定された忘却ρ係数の値が果たして最適な値であるかどうか判定する手段も無かった。
また、カルマンフィルタで用いる誤差共分散行列は、本来、零でない任意のベクトルとの２次形式が常に正（以下、「正定」という。）であるがコンピュータにより単精度で計算した場合にはその２次形式が負（以下、「負定」という。）となり、数値的に不安定になることが知られている。また、計算量がＯ（Ｎ^２）（あるいはＯ（Ｎ^３））であるため、状態ベクトルｘ_ｋの次元Ｎが大きい場合、１ステップ当たりの演算回数が急激に増大し、実時間処理には適さなかった。
本発明は、以上の点に鑑み、忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム推定方法を提供することを目的とする。
本発明は、上記課題を解決するために、新たに考案したＨ_∞最適化手法を用いて忘却係数が最適決定可能な状態推定アルゴリズムを導出する。さらに、常に正定であるべき誤差共分散行列の代わりに、その因数行列を更新することによって数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。
本発明の第１の解決手段によると、
次式で表される状態空間モデルに対して、
ｘ_ｋ＋１＝Ｆ_ｋｘ_ｋ＋Ｇ_ｋｗ_ｋ
ｙ_ｋ＝Ｈ_ｋｘ_ｋ＋ｖ_ｋ
ｚ_ｋ＝Ｈ_ｋｘ_ｋ
ここで、
ｘ_ｋ：状態ベクトルまたは単に状態
ｗ_ｋ：システム雑音
ｖ_ｋ：観測雑音
ｙ_ｋ：観測信号
ｚ_ｋ：出力信号
Ｆ_ｋ：システムのダイナミックス
Ｇ_ｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γ_ｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法及びプログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部は、上限値γ_ｆ、フィルタの入力である観測信号ｙ_ｋ、観測行列Ｈ_ｋを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γ_ｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Ｈ_ｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ_∞フィルタを実行するステップと、
ｘ＾_ｋ｜ｋ＝Ｆ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}＋Ｋ_ｓ，ｋ（ｙ_ｋ−Ｈ_ｋＦ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}）
ここで、
ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
Ｆ_ｋ：システムのダイナミックス
Ｋ_ｓ，ｋ：フィルタゲイン
処理部は、ハイパーＨ_∞フィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、観測行列Ｈ_ｉとフィルタゲインＫ_ｓ，ｉにより、前記上限値γ_ｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γ_ｆを小さくしていき前記ハイパーＨ_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法、各ステップをコンピュータに実行させるためのシステム推定プログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
また、本発明の第２の解決手段によると、
次式で表される状態空間モデルに対して、
ｘ_ｋ＋１＝Ｆ_ｋｘ_ｋ＋Ｇ_ｋｗ_ｋ
ｙ_ｋ＝Ｈ_ｋｘ_ｋ＋ｖ_ｋ
ｚ_ｋ＝Ｈ_ｋｘ_ｋ
ここで、
ｘ_ｋ：状態ベクトルまたは単に状態
ｗ_ｋ：システム雑音
ｖ_ｋ：観測雑音
ｙ_ｋ：観測信号
ｚ_ｋ：出力信号
Ｆ_ｋ：システムのダイナミックス
Ｇ_ｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γ_ｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び／又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
前記処理部は、上限値γ_ｆ、フィルタの入力である観測信号ｙ_ｋ、観測行列Ｈ_ｋを含む値を記憶部又は入力部から入力すること、
前記処理部は、前記上限値γ_ｆに従い、状態空間モデルに関連する忘却係数ρを決定すること、
前記処理部は、記憶部から初期値又はある時刻の観測行列Ｈ_ｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ_∞フィルタを実行すること、
ｘ＾_ｋ｜ｋ＝Ｆ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}＋Ｋ_ｓ、ｋ（ｙ_ｋ−Ｈ_ｋＦ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}）
ここで、
ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
Ｆ_ｋ：システムのダイナミックス
Ｋ_ｓ，ｋ：フィルタゲイン
前記処理部は、ハイパーＨ_∞フィルタに関する求められた値を記憶部に記憶すること、
前記処理部は、観測行列Ｈ_ｉとフィルタゲインＫ_ｓ，ｉにより、前記上限値γ_ｆ及び前記忘却係数ρに基づく存在条件を計算すること、
前記処理部は、上限値γ_ｆを小さくしていき前記ハイパーＨ_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること
を備えた前記システム推定装置が提供される。
本発明の推定方法は忘却係数を最適に決定することが可能であり、かつアルゴリズムは単精度でも安定に動作可能であるため、低コストで高い性能が実現できる。一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。

図１は、本実施の形態に関するハードウェアの構成図である。
図２は、Ｈ_∞フィルタのロバスト化と忘却係数ρの最適化についてのフローチャートである。
図３は、図２中のＨ_∞フィルタ（Ｓ１０５）のアルゴリズムについてのフローチャートである。
図４は、定理２の平方根アレイアルゴリズムの説明図である。
図５は、定理３の数値的に安定な高速アルゴリズムのフローチャートである。
図６は、インパルス応答｛ｈ_ｉ｝_ｉ＝_０ ^２３の値を示す図である。
図７は、定理３の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果である。
図８は、システム推定のための構成図である。
図９は、通信系とエコーについての説明図である。
図１０は、エコーキャンセラの原理図である。

以下に、本発明の実施の形態について説明する。
１．記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
ｘ_ｋ：状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。
ｘ_０：初期状態；未知である。
ｗ_ｋ：システム雑音；未知である。
ｖ_ｋ：観測雑音；未知である。
ｙ_ｋ：観測信号；フィルタの入力となり、既知である。
ｚ_ｋ：出力信号；未知である。
Ｆ_ｋ：システムのダイナミックス；既知である。
Ｇ_ｋ：駆動行列；実行時に既知となる。
Ｈ_ｋ：観測行列；既知である。
ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値；フィルタ方程式によって与えられる。
ｘ＾_{ｋ＋１｜ｋ}：観測記号ｙ_０〜ｙ_ｋまで用いた時刻ｋ＋１の状態ｘ_ｋ＋１の推定値；フィルタ方程式によって与えられる。
ｘ＾_０｜０：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Σ＾_ｋ｜ｋ：ｘ＾_ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾_{ｋ＋１｜ｋ}：ｘ＾_{ｋ＋１｜ｋ}の誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾_１｜０：初期状態の共分散行列に対応；本来未知であるが、便宜上ε_０Ｉが用いられる。
Ｋ_ｓ，ｋ：フィルタゲイン；行列Σ＾_{ｋ｜ｋ−１}から得られる。
ρ：忘却係数；定理１〜３の場合、γ_ｆが決まればρ＝１−χ（γ_ｆ）より自動的に決定される。
ｅ_ｆ，ｉ：フィルタ誤差
Ｒ_ｅ，ｋ：補助変数
なお、記号の上に付される”＾”、”ｖ”は、推定値の意味である。また、”〜”、”−”、”Ｕ”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、ｘ、ｗ、Ｈ、Ｇ，Ｋ、Ｒ、Σ、等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
２．システム推定のハードウェア及びプログラム
本システム推定方法又はシステム推定装置・システムは、その各手順をコンピュータに実行させるためのシステム推定プログラム、システム推定プログラムを記録したコンピュータ読み取り可能な記録媒体、システム推定プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
図１は、本実施の形態に関するハードウェアの構成図である。
このハードウェアは、中央処理装置（ＣＰＵ）である処理部１０１、入力部１０２、出力部１０３、表示部１０４及び記憶部１０５を有する。また、処理部１０１、入力部１０２、出力部１０３、表示部１０４及び記憶部１０５は、スター又はバス等の適宜の接続手段で接続されている。記憶部１０５は、システム推定される「１．記号の説明」で示した既知のデータが必要に応じて記憶される。また、未知・既知のデータや計算されたハイパーＨ_∞フィルタに関するデータ・その他のデータが処理部１０１により、必要に応じて書込み及び／又は読み出しされる。
３．忘却係数が最適決定可能なハイパーＨ_∞フィルタ
（定理１）
次式のような状態空間モデルを考える。

このような状態空間モデルに対して、次式のようなＨ_∞評価基準を提案する。

このＨ_∞評価基準を満たす状態推定値ｘ＾_ｋ｜ｋ（あるいは出力推定値ｚ^ｖ _ｋ｜ｋ）は、次のレベルγ_ｆのハイパーＨ_∞フィルタによって与えられる。

ここで、

である。なお、式（１１）はフィルタ方程式、式（１２）はフィルタゲイン、式（１３）はリカッチ方程式をそれぞれ示す。
また、駆動行列Ｇ_ｋは次のように生成される。

また、上述の高速Ｈ_∞フィルタの追従能力を向上するためには、上限値γ_ｆは次の存在条件を満たすように出来るだけ小さく設定する。

ただし、χ（γ_ｆ）は、χ（１）＝１、χ（∞）＝０を満たすγ_ｆの単調減衰関数である。
定理１の特徴は状態推定のロバスト化と忘却係数ρの最適化を同時に行っている点にある。
図２に、Ｈ_∞フィルタのロバスト化と忘却係数ρの最適化についてのフローチャートを示す。ここで、
ブロック「ＥＸＣ＞０」：Ｈ_∞フィルタの存在条件、
Δγ：正の実数、である。
まず、処理部１０１は、記憶部１０５又は入力部１０２から上限値γ_ｆを読み出し又は入力する（Ｓ１０１）。この例では、γ_ｆ＞＞１が与えられる。処理部１０１は、式（１５）によって忘却係数ρを決定する（Ｓ１０３）。その後、処理部１０１は、忘却係数ρに基づき、式（１０）〜式（１３）のハイパーＨ_∞フィルタを実行する（Ｓ１０５）。処理部１０１は、式（１７）（あるいは、後述の式（１８））の右辺（これをＥＸＣとする）を計算し（Ｓ１０７）、その存在条件がすべての時刻で満たされれば（Ｓ１０９）、γ_ｆをΔγだけ小さくして同じ処理を繰り返す（Ｓ１１１）。一方、あるγ_ｆで存在条件を満たさなくなったとき（Ｓ１０９）、そのγ_ｆにΔγを加えたものをγ_ｆの最適値γ_ｆ ^ｏｐとして出力部１０３に出力及び／又は記憶部１０５に記憶する（Ｓ１１３）。なお、この例では、Δγを加えているが、それ以外の予め設定された値を加えてもよい。この最適化のプロセスをγ−イタレーションと呼ぶ。なお、処理部１０１は、Ｈ_∞フィルタ計算ステップＳ１０５及び存在条件の計算ステップＳ１０７等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。
ハイパーＨ_∞フィルタが存在条件を満たすとき、式（９）の不等式は常に満たされる。よって、式（９）の分母の外乱エネルギーが限定される場合、式（９）の分子の２乗推定誤差の総和は有界となり、ある時刻以降の推定誤差が０になる。これは、γ_ｆをより小さく出来れば、状態ｘ_ｋの変化に推定値ｘ＾_ｋ｜ｋが速やかに追従できることを意味する。
ここで、定理１のハイパーＨ_∞フィルタのアルゴリズムは通常のＨ_∞フィルタのものとは異なることに注意されたい。また、γ_ｆ→∞のとき、ρ＝１、Ｇ_ｋ＝０となり、定理１のＨ_∞フィルタのアルゴリズムはカルマンフィルタのアルゴリズムと一致する。
図３に、図２中の（ハイパー）Ｈ_∞フィルタ（Ｓ１０５）のアルゴリズムについてのフローチャートを示す。
ハイパーＨ_∞フィルタリングアルゴリズムは以下のように要約することができる。
［ステップＳ２０１］処理部１０１は、記憶部１０５から再帰式の初期条件を読み出し、又は、初期条件を入力部１０２から入力し、図示のように定める。なお、Ｌはあらかじめ定められた最大データ数を示す。
［ステップＳ２０３］処理部１０１は、時刻ｋと最大データ数Ｌとを比較する。処理部１０１は、時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、必要に応じて再スタートを行ってもよい。）
［ステップＳ２０５］処理部１０１は、フィルタゲインＫ_ｓ，ｋを式（１２）を用いて計算する。
［ステップＳ２０７］処理部１０１は、式（１１）のハイパーＨ_∞フィルタのフィルタ方程式を更新する。
［ステップＳ２０９］処理部１０１は、誤差の共分散行列に対応する項Σ＾_ｋ｜ｋ、Σ＾_{ｋ＋１｜ ｋ}を式（１３）のリカッチ方程式を用いて計算する。
［ステップＳ２１１］時刻ｋを進ませて（ｋ＝ｋ＋１）、ステップＳ２０３に戻り、データがある限り続ける。
なお、処理部１０１は、Ｈ_∞フィルタ計算ステップＳ２０５〜Ｓ２０９等の各ステップで求めた適宜の中間値及び最終値、存在条件の値等を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。
（スカラー存在条件）
ところで、式（１７）の存在条件の判定にはＯ（Ｎ^２）の計算量が必要であった。しかし、次の条件を用いれば計算量Ｏ（Ｎ）で定理１のＨ_∞フィルタの存在性、すなわち式（９）を検証することができる。
系１：スカラー存在条件
次の存在条件を用いれば計算量Ｏ（Ｎ）でハイパーＨ_∞フィルタの存在性が判定できる。

ここで、

ただし、Ｋ_ｓ，ｉは式（１２）で求めたフィルタゲインである。
（証明）
以下に系１の証明を説明する。
２×２の行列Ｒ_ｅ，ｋの特性方程式

を解けば、Ｒ_ｅ，ｋの固有値λ_ｉが次のように得られる。

もし、

であれば、行列Ｒ_ｅ，ｋの２つの固有値の１つは正となり、もう１つは負となり、行列Ｒ_ｋとＲ_ｅ，ｋは同じイナーシャをもつ。これより、

を用いれば、式（１８）の存在条件が得られる。ここで、Ｈ_ｋＫ_ｓ，ｋの計算量はＯ（Ｎ）である。
４．数値的に安定な状態推定アルゴリズム
上述のハイパーＨ_∞フィルタは、Σ＾_{ｋ｜ｋ−１}∈Ｒ^Ｎ×Ｎを更新するため、単位時間ステップ当たりの計算量はＯ（Ｎ^２）となる、すなわち、Ｎ^２に比例する算術演算が必要となる。ここで、Ｎは状態ベクトルｘ_ｋの次元である。よって、ｘ_ｋの次元が増加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列Σ＾_{ｋ｜ｋ−１}は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度（例：３２ｂｉｔ）でも動作可能な状態推定アルゴリズムの開発が望まれる。
そこで、次に、
Ｒ_ｋ＝Ｒ^１／２ _ｋＪ_１Ｒ^１／２ _ｋ、
Ｒ_ｅ，ｋ＝Ｒ^１／２ _ｅ，ｋＪ_１Ｒ^１／２ _ｅ，ｋ、
Σ＾_{ｋ｜ｋ−１}＝Σ＾^１／２ _{ｋ｜ｋ−１}Σ＾^１／２ _{ｋ｜ｋ−１}
に着目して、数値的に安定化した定理１のＨ_∞フィルタ（平方根アレイアルゴリズム）を定理２に示す。ただし、ここでは簡単のためＦ_ｋ＝Ｉとしたが、Ｆ_ｋ≠Ｉの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパーＨ_∞フィルタを示す。
（定理２）

ただし、

なお、式（２１）、（２２）において、Ｊ_１ ^−１およびＪ_１は削除可能である。
図４に、定理２の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図２に示した定理１のフローチャート中のＨ_∞フィルタの計算（Ｓ１０５）で用いることができる。
本推定アルゴリズムは、Σ＾_{ｋ｜ｋ−１}をリカッチ型の更新式で求める代わりに、その因数行列Σ＾^１／２ _{ｋ｜ｋ−１}∈Ｒ^Ｎ×Ｎ（Σ＾_{ｋ｜ｋ−１}の平方根行列）をＪ−ユニタリ変換に基づく更新式で求めている。このとき生じる１−１ブロック行列と２−１ブロック行列からフィルタゲインＫ_ｓ，ｋを図示のように求めている。このため、Σ＾_{ｋ｜ｋ−１}＝Σ＾^１／２ _{ｋ｜ｋ−１}Σ＾^１／２ _{ｋ｜ｋ− １}＞０となり、Σ＾_{ｋ｜ｋ−１}の正定性は保証され、数値的に安定化できる。なお、定理２のＨ_∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ^２）のままである。
なお、図４において、Ｊ_１ ^−１は削除可能である。
まず、処理部１０１は、式（２２）の左辺の行列式の各要素に含まれる項を記憶部１０５から読み出し又は内部メモリ等から得て、Ｊ−ユニタリ変換を実行する（Ｓ３０１）。処理部１０１は、求めた式（２２）の右辺の行列式の要素からシステムゲインＫ_ｋ、Ｋ_ｓ，ｋを式（２１）に基づき計算する（Ｓ３０３、Ｓ３０５）。処理部１０１は、式（２０）に基づき状態推定値ｘ＾_ｋ｜ｋを計算する（Ｓ３０７）。
５．状態推定のための数値的に安定な高速アルゴリズム
上述のように、定理２のＨ_∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ^２）のままである。そこで、計算量の対策として、Ｈ _ｋ＝Ｈ _ｋ＋１Ψ，Ｈ _ｋ＝［ｕ（ｋ），・・・，ｕ（０），０，・・・，０］のとき、ｘ _ｋ＝［ｘ^Ｔ _ｋ，０^Ｔ］^Ｔの１ステップ予測誤差の共分散行列Σ _{ｋ＋１｜ｋ}が

を満たすことを利用して、Σ _{ｋ＋１｜ｋ}の代わりに次元の低いＬ _ｋ（すなわちＬ^〜 _ｋ）を更新する

得られる。
（定理３）

ただし、

なお、定理３の証明は、後述する。
上式は、Ｋ⁻ _ｋ（＝Ｐ^−１／２Ｋ_ｋ）の代わりにＫ_ｋについても整理することができる。
さらに、次のＪ−ユニタリ行列

を用いれば定理４の高速化した状態推定アルゴリズムが得られる。ただし、Ψはシフト行列を表す。
（定理４）

ただし、

であり、ｄｉａｇ｛・｝は対角行列、Ｒ_{ｅ，ｋ＋１}（１，１）は行列Ｒ_{ｅ，ｋ＋１}の１−１成分をそれぞれ表す。また、上式はＫ⁻ _ｋの代わりにＫ_ｋに関しても整理できる。
本高速アルゴリズムは、次の因数分解

におけるＬ^〜 _ｋ∈Ｒ^{（Ｎ＋１）×２}の更新によってフィルタゲインＫ_ｓ，ｋを求めているので、単位ステップ当たりの計算量はＯ（Ｎ＋１）で済む。ここで、次式に注意されたい。

図５に、定理３の数値的に安定な高速アルゴリズムのフローチャートを示す。この高速アルゴリズムは図２のＨ_∞フィルタの計算ステップ（Ｓ１０５）に組み込まれ、γ−イタレーションによって最適化される。よって、存在条件が満たされる間はγ_ｆは除々に減少されるが、満たされなくなった時点で、図示のようにγ_ｆは増加される。
Ｈ_∞フィルタリングアルゴリズムは以下のように要約することができる。
［ステップＳ４０１］処理部１０１は、再帰式の初期条件を図示のように定める。なお、Ｌは最大データ数を示す。
［ステップＳ４０３］処理部１０１は、時刻ｋと最大データ数Ｌとを比較する。処理部１０１は、時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、再スタートする。）
［ステップＳ４０５］処理部１０１は、フィルタゲインに対応する項Ｋ_ｋ＋１を式（２７）、（３１）を用いて再帰的に計算する。
［ステップＳ４０６］処理部１０１は、Ｒ_{ｅ，ｋ＋１}を式（２９）を用いて再帰的に計算する。
［ステップＳ４０７］処理部１０１は、さらにＫ_ｓ，ｋを式（２６）、（３１）を用いて計算する。
［ステップＳ４０９］処理部１０１は、ここで、存在条件ＥＸＣ＞０を判定し、存在条件を満たせばステップＳ４１１に進む。
［ステップＳ４１３］一方、処理部１０１は、ステップＳ４０９で存在条件を満たさなければγ_ｆを増加し、ステップＳ４０１に戻る。
［ステップＳ４１１］処理部１０１は、式（２５）のＨ_∞フィルタのフィルタ方程式を更新する。
［ステップＳ４１５］処理部１０１は、Ｒ_{ｒ，ｋ＋１}を式（３０）を用いて再帰的に計算する。また、処理部１０１は、Ｌ^〜 _ｋ＋１を式（２８）、（３１）を用いて再帰的に計算する。
［ステップＳ４１９］処理部１０１は、時刻ｋを進ませて（ｋ＝ｋ＋１）、ステップＳ４０３に戻り、データがある限り続ける。
なお、処理部１０１は、Ｈ_∞フィルタ計算ステップＳ４０５〜Ｓ４１５及び存在条件の計算ステップＳ４０９等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。
６．エコーキャンセラ
つぎに、エコーキャンセリング問題の数理モデルを作成する。
まず、受信信号｛ｕ_ｋ｝がエコーパスへの入力信号となることを考慮すれば、エコーパスの（時変）インパルス応答｛ｈ_ｉ［ｋ］｝により、エコー｛ｄ_ｋ｝の観測値｛ｙ_ｋ｝は次式で表される。

ここで、ｕ_ｋ，ｙ_ｋはそれぞれ時刻ｔ_ｋ（＝ｋＴ；Ｔは標本化周期）における受信信号とエコーを表し、ｖ_ｋは時刻ｔ_ｋにおける平均値０の回線雑音とし、ｈ_ｉ［ｋ］，ｉ＝０，・・・，Ｎ−１は、時変インパルス応答であり、そのタップ数Ｎは既知とする。このとき、インパルス応答の推定値｛ｈ＾_ｉ［ｋ］｝が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。

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

ただし、

このような状態空間モデルに対するハイパーおよび高速Ｈ_∞フィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
７．インパルス応答に対する評価
（動作の確認）
エコーパスのインパルス応答が時間的に不変であり（ｈ_ｉ［ｋ］＝ｈ_ｉ）、かつそのタップ数Ｎが４８である場合について、シミュレーションを用いて、本高速アルゴリズムの動作を確認する。

なお、図６は、ここでのインパルス応答｛ｈ_ｉ｝の値を示す図である。
ここで、インパルス応答｛ｈ_ｉ｝_ｉ＝０ ^２３は、図示の値を採用し、その他｛ｈ_ｉ｝_ｉ＝_２４ ^４７は０とする。また、υ_ｋは平均値０、分散σ_ｖ ^２＝１．０×１０^−６の定常なガウス白色雑音とし、標本化周期Ｔを便宜上１．０とする。
また、受信信号｛ｕ_ｋ｝は次のように２次のＡＲモデルで近似する。

ただし、α_１＝０．７，α_２＝０．１とし、ｗ_ｋ′は平均値０、分散σ_ｗ′ ^２＝０．０４の定常なガウス白色雑音とする。
（インパルス応答の推定結果）
図７に、定理３の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図７（ｂ）の縦軸は、
√｛Σ_ｉ＝０ ^４７（ｈ_ｉ−ｘ＾_ｋ（ｉ＋１））^２｝
を表す。
これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし、ρ＝１−χ（γ_ｆ）、χ（γ_ｆ）＝γ_ｆ ^−２、ｘ＾_０｜０＝０、Σ＾_１｜０＝２０Ｉとし、計算は倍精度で行った。また、存在条件を確認しつつ、γ_ｆ＝５．５と設定とした。
８．定理の証明
８−１．定理２の証明
次の関係式

が成り立つとき、両辺の２×２ブロック行列の各項を比較すれば次式が得られる。

これは定理１のＦ_ｋ＝Ｉのときの式（１３）のリカッチ方程式と一致する。ただし、

一方、ＡＪＡ^Ｔ＝ＢＪＢ^Ｔが成り立つとき、ＢはＪ−ユニタリ行列Θ（ｋ）を用いてＢ＝ＡΘ（ｋ）と表すことができる。よって、式（４０）より定理１のリカッチ方程式は次式と等価である。

なお、式（４０）、（４５）において、Ｊ_１ ^−１は削除可能である。
８−２．定理３の証明
次のようにブロック三角化するＪ−ユニタリ行列Θ（ｋ）が存在すると仮定する。

のように決定することができる。ただし、Ｓは対角成分が１または−１をとる対角行列とする。
（１，１）−ブロック行列

（２，１）−ブロック行列

（２，２）−ブロック行列

これより、

８−３．定理４の証明
観測行列Ｈ_ｋがシフト特性をもち、かつ

のとき、定理２と同様な方法によって次の関係式が得られる。

ただし、

となるようにＲ_{ｒ，ｋ＋１}を決定する。次に、式（４６）の３行目にＲ_{ｒ，ｋ＋１}の更新式を新たに追加すれば、最終的に次式が得られる。

この両辺の３×２ブロック行列の各項の対応から次のゲイン行列Ｋ⁻ _ｋの更新式が得られる。

一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができる。

【００１０】
ｙ_ｋ：観測信号；フィルタの入力となり、既知である。
ｚ_ｋ：出力信号；未知である。
Ｆ_ｋ：システムのダイナミックス；既知である。
Ｇ_ｋ：駆動行列；実行時に既知となる。
Ｈ_ｋ：観測行列；既知である。
ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値；フィルタ方程式によって与えられる。
ｘ＾_{ｋ＋１｜ｋ}：観測記号ｙ_０〜ｙ_ｋで用いた時刻ｋ＋１の状態ｘ_ｋ＋１の推定値；フィルタ方程式によって与えられる。
ｘ＾_０｜０：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Σ＾_ｋ｜ｋ：ｘ＾_ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾_{ｋ＋１｜ｋ}：ｘ＾_{ｋ＋１｜ｋ}の誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾_１｜０：初期状態の共分散行列に対応；本来未知であるが、便宜上ε_０Ｉが用いられる。
Ｋ_ｓ，ｋ：フィルタゲイン；行列Σ＾_{ｋ｜ｋ−１}から得られる。
ρ：忘却係数；定理１〜３の場合、γ_ｆが決まればρ＝１−χ（γ_ｆ）より自動的に決定される。
ｅ_ｆ，ｉ：フィルタ誤差
Ｒ_ｅ，ｋ：補助変数
なお、記号の上に付される”＾”、”ｖ”は、推定値の意味である。また、”〜”、”−”、”Ｕ”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、ｘ、ｗ等はベクトル、Ｈ、Ｇ，Ｋ、Ｒ、Σ、等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。

【００１６】
る。
４．数値的に安定な状態推定アルゴリズム
上述のハイパーＨ_∞フィルタは、Σ＾_{ｋ｜ｋ−１}∈Ｒ^Ｎ×Ｎを更新するため、単位時間ステップ当たりの計算量はＯ（Ｎ^２）となる、すなわち、Ｎ^２に比例する算術演算が必要となる。ここで、Ｎは状態ベクトルｘ_ｋの次元である。よって、ｘ_ｋの次元が増加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列Σ＾_{ｋ｜ｋ−１}は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度（例：３２ｂｉｔ）でも動作可能な状態推定アルゴリズムの開発が望まれる。
そこで、次に、
Ｒ_ｋ＝Ｒ^１／２ _ｋＪ_１Ｒ^Ｔ／２ _ｋ、
Ｒ_ｅ，ｋ＝Ｒ^１／２ _ｅ，ｋＪ_１Ｒ^Ｔ／２ _ｅ，ｋ、
Σ＾_{ｋ｜ｋ−１}＝Σ＾^１／２ _{ｋ｜ｋ−１}Σ＾^Ｔ／２ _{ｋ｜ｋ−１}
に着目して、数値的に安定化した定理１のＨ_∞フィルタ（平方根アレイアルゴリズム）を定理２に示す。ただし、ここでは簡単のためＦ_ｋ＝Ｉとしたが、Ｆ_ｋ≠Ｉの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパーＨ_∞フィルタを示す。

【００１７】
（定理２）

ただし、

位行列である。また、Ｋ_ｋ（：，１）は行列Ｋ_ｋの１列目の列ベクトルを表す。
なお、式（２１）、（２２）において、Ｊ_１ ^−１およびＪ_１は削除可能である。
図４に、定理２の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図２に示した定理１のフローチャート中のＨ_∞フィルタの計算（Ｓ１０５）で用いることができる。
本推定アルゴリズムは、Σ＾_{ｋ｜ｋ−１}をリカッチ型の更新式で求める代わりに、その因数行列Σ＾^１／２ _{ｋ｜ｋ−１}∈Ｒ^Ｎ×Ｎ（Σ＾_{ｋ｜ｋ−１}の平方根行列）をＪ−ユニタリ変換に基づく更新式で求めている。このとき生じる１−１ブロック行列と２−１ブロック行列からフィルタゲインＫ_ｓ，ｋを図示のように求めている。このため、Σ＾_{ｋ｜ｋ−１}＝Σ＾^１／２ _{ｋ｜ｋ−１}Σ＾^Ｔ／２ _{ｋ｜ｋ−１}＞０となり、Σ＾_{ｋ｜ｋ−１}の正定性は保証され、数値的に安定化できる。なお、定理２のＨ_∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ^２）のままである。
なお、図４において、Ｊ_１ ^−１は削除可能である。
まず、処理部１０１は、式（２２）の左辺の行列式の各要素に含まれる項を記憶部１０５から読み出し又は内部メモリ等から得て、Ｊ−ユニタリ変換を実行する（Ｓ３０１）。処理部１０１は、求めた式（２２）の右辺の行列式の要素からシステムゲインＫ_ｋ、Ｋ_ｓ，ｋを

【００１８】
式（２１）に基づき計算する（Ｓ３０３、Ｓ３０５）。処理部１９１は、式（２０）に基づき状態推定値ｘ＾_ｋ｜ｋを計算する（Ｓ３０７）。
５．状態推定のための数値的に安定な高速アルゴリズム
上述のように、定理２のＨ_∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ^２）のままである。そこで、計算量の対策として、Ｈ _ｋ＝Ｈ _ｋ＋１Ψ，Ｈ _ｋ＝［ｕ（ｋ），・・・，ｕ（０），０，・・・，０］のとき、ｘ _ｋ＝［ｘ^Ｔ _ｋ，０^Ｔ］^Ｔの１ステップ予測誤差の共分散行列Σ _{ｋ＋１｜ｋ}が

を満たすことを利用して、Σ _{ｋ＋１｜ｋ}代わりに次元の低いＬ _ｋ（すなわちＬ^〜 _ｋ）を更新する

３が得られる。
（定理３）

ただし、

なお、定理３の証明は、後述する。
上式は、Ｋ⁻ _ｋ（＝ρ^−１／２Ｋ_ｋ）の代わりにＫ_ｋについても整理することができる。

【００１９】
さらに、次のＪ−ユニタリ行列

を用いれば定理４の高速化した状態推定アルゴリズムが得られる。ただし、Ψはシフト行列を表す。
（定理４）

ただし、

であり、ｄｉａｇ｛・｝は対角行列、Ｒ_{ｅ，ｋ＋１}（１，１）は行列Ｒ_{ｅ，ｋ＋１}の１−１成分をそれぞれ表す。また、上式はＫ⁻ _ｋの代わりにＫ_ｋに関しても整理できる。
本高速アルゴリズムは、次の因数分解

におけるＬ^〜 _ｋ∈Ｒ^{（Ｎ＋１）×２}の更新によってフィルタゲインＫ_ｓ，ｋを求めているので、単位ス

【００２０】
テップ当たりの計算量はＯ（Ｎ＋１）で済む。ここで、次式に注意されたい。

図５に、定理４の数値的に安定な高速アルゴリズムのフローチャートの一例を示す。この高速アルゴリズムは図２のＨ_∞フィルタの計算ステップ（Ｓ１０５）に組み込まれ、γ−イタレーションによって最適化される。よって、存在条件が満たされる間はγ_ｆは除々に減少されるが、満たされなくなった時点で、図示のようにγ_ｆは増加される。
Ｈ_∞フィルタリングアルゴリズムは以下のように要約することができる。
［ステップＳ４０１］処理部１０１は、再帰式の初期条件を図示のように定める。なお、Ｌは最大データ数を示す。
［ステップＳ４０３］処理部１０１は、時刻ｋと最大データ数Ｌとを比較する。処理部１０１は、時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、再スタートする。）
［ステップＳ４０５］処理部１０１は、フィルタゲインに対応する項Ｋ_ｋ＋１を式（２７）、（３１）を用いて再帰的に計算する。
［ステップＳ４０６］処理部１０１は、Ｒ_{ｅ，ｋ＋１}を式（２９）を用いて再帰的に計算する。
［ステップＳ４０７］処理部１０１は、さらにＫ_ｓ，ｋを式（２６）、（３１）を用いて計算する。
［ステップＳ４０９］処理部１０１は、ここで、存在条件ＥＸＣ＞０を判定し、存在条件を満たせばステップＳ４１１に進む。
［ステップＳ４１３］一方、処理部１０１は、ステップＳ４０９で存在条件を満たさなければγ_ｆを増加し、ステップＳ４０１に戻る。
［ステップＳ４１１］処理部１０１は、式（２５）のＨ_∞フィルタのフィルタ方程式を更新する。
［ステップＳ４１５］処理部１０１は、Ｒ_{ｒ．ｋ＋１}を式（３０）を用いて再帰的に計算する。ま

【００２３】
周期Ｔを便宜上１．０とする。
また、受信信号｛ｕ_ｋ｝は次のように２次のＡＲモデルで近似する。

ただし、α_１＝０．７，α_２＝０．１とし、ｗ_ｋ′は平均値０、分散σｗ_′ ^２＝０．０４の定常なガウス白色雑音とする。
（インパルス応答の推定結果）
図７に、定理４の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図７（ｂ）の縦軸は、
√｛Σ_ｉ＝０ ^４７（ｈ_ｉ−ｘ＾_ｋ（ｉ＋１））^２｝
を表す。
これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし、ρ＝１−χ（γ_ｆ）、χ（γ_ｆ）＝γ_ｆ ^−２、ｘ＾_０｜０＝０、Σ＾_１｜０＝２０Ｉとし、計算は倍精度で行った。また、存在条件を確認しつつ、γ_ｆ＝５．５と設定とした。
８．定理の証明
８−１．定理２の証明
次の関係式

が成り立つとき、両辺の２×２ブロック行列の各項を比較すれば次式が得られる。

【００２５】
（１，１）−ブロック行列

（２，１）−ブロック行列

【００２６】
（２，２）−ブロック行列

これより、

８−３．定理４の証明
観測行列Ｈ_ｋがシフト特性をもち、かつ

のとき、定理２と同様な方法によって次の関係式が得られる。

ただし、

となるようにＲ_{ｒ，ｋ＋１}を決定する。次に、式（４６）の３行目にＲ_{ｒ，ｋ＋１}の更新式を新たに追

Claims (17)

1. 次式で表される状態空間モデルに対して、
   ｘ_ｋ＋１＝Ｆ_ｋｘ_ｋ＋Ｇ_ｋｗ_ｋ
   ｙ_ｋ＝Ｈ_ｋｘ_ｋ＋ｖ_ｋ
   ｚ_ｋ＝Ｈ_ｋｘ_ｋ
   ここで、
   ｘ_ｋ：状態ベクトルまたは単に状態
   ｗ_ｋ：システム雑音
   ｖ_ｋ：観測雑音
   ｙ_ｋ：観測信号
   ｚ_ｋ：出力信号
   Ｆ_ｋ：システムのダイナミックス
   Ｇ_ｋ：駆動行列
   評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γ_ｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
   処理部は、上限値γ_ｆ、フィルタの入力である観測信号ｙ_ｋ、観測行列Ｈ_ｋを含む値を記憶部又は入力部から入力するステップと、
   処理部は、前記上限値γ_ｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
   処理部は、記憶部から初期値又はある時刻の観測行列Ｈ_ｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ_∞フィルタを実行するステップと、
   ｘ＾_ｋ｜ｋ＝Ｆ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}＋Ｋ_ｓ，ｋ（ｙ_ｋ−Ｈ_ｋＦ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}）
   ここで、
   ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
   Ｋ_ｓ，ｋ：フィルタゲイン
   処理部は、ハイパーＨ_∞フィルタに関する求められた値を記憶部に記憶するステップと、
   処理部は、求められた観測行列Ｈ_ｉ、又は、観測行列Ｈ_ｉとフィルタゲインＫ_ｓ，ｉにより、前記上限値γ_ｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
   処理部は、上限値γ_ｆを小さくしていき前記ハイパーＨ_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
   を含む前記システム推定方法。
2. 処理部は、前記存在条件を次式に従い計算する請求項１に記載のシステム推定方法。
   

   
3. 処理部は、前記存在条件を次式に従い計算する請求項１に記載のシステム推定方法。
   

   

   ここで、
   

   
4. 前記忘却係数ρ及び前記上限値γ_ｆは、次式の関係である請求項１に記載のシステム推定方法。
   ０＜ρ＝１−χ（γ_ｆ）≦１（ただし、χ（γ_ｆ）は、χ（１）＝１、χ（∞）＝０を満たすγ_ｆの単調減衰関数）
5. 前記ハイパーＨ_∞フィルタを実行するステップは、
   処理部は、前記フィルタゲインＫ_ｓ，ｋを、次式により求めることを特徴とする請求項１に記載のシステム推定方法。
   

   

   ここで、
   

   

   なお、式（１６）の右辺はより一般化することもできる。
   ここで、
   ｘ_ｋ：状態ベクトルまたは単に状態
   ｙ_ｋ：観測信号
   ｚ_ｋ：出力信号
   Ｆ_ｋ：システムのダイナミックス
   Ｈ_ｋ：観測行列
   ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
   Σ＾_ｋ｜ｋ：ｘ＾_ｋ｜ｋの誤差の共分散行列に対応
   Ｋ_ｓ，ｋ：フィルタゲイン
   ｅ_ｆ，ｉ：フィルタ誤差
   Ｒ_ｅ，ｋ：補助変数
6. 前記ハイパーＨ_∞フィルタを実行するステップは、
   処理部は、初期条件に基づき、フィルタゲインＫ_ｓ，ｋを前記式（１２）を用いて計算するステップと、
   処理部は、前記式（１１）のＨ_∞フィルタのフィルタ方程式を更新するステップと、
   処理部は、Σ＾_ｋ｜ｋ、Σ＾_{ｋ＋１｜ｋ}を前記式（１３）を用いて計算するステップと、
   処理部は、前記各ステップを、時刻ｋを進ませて繰り返し実行するステップと
   を含む請求項５に記載のシステム推定方法。
7. 前記ハイパーＨ_∞フィルタを実行するステップは、
   処理部は、前記フィルタゲインＫ_ｓ，ｋを、ゲイン行列Ｋ_ｋを用いて、次式により求めることを特徴とする請求項１に記載のシステム推定方法。
   

   

   ただし、
   

   

   

   なお、式（２１）、（２２）において、Ｊ_１ ^−１およびＪ_１は削除可能である。
   ここで、
   ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋでを用いた時刻ｋの状態ｘ_ｋの推定値
   ｙ_ｋ：観測信号
   Ｆ_ｋ：システムのダイナミックス
   Ｋ_ｓ，ｋ：フィルタゲイン
   Ｈ_ｋ：観測行列
   Σ＾_ｋ｜ｋ：ｘ＾_ｋ｜ｋの誤差の共分散行列に対応
   Θ（ｋ）：Ｊ−ユニタリ行列
   Ｒ_ｅ，ｋ：補助変数
8. 前記ハイパーＨ_∞フィルタを実行するステップは、
   処理部は、Ｋ_ｋ、Σ＾_{ｋ＋１｜ｋ} ^１／２を前記式（２２）を用いて計算するステップと、
   処理部は、初期条件に基づき、フィルタゲインＫ_ｓ，ｋを前記式（２１）を用いて計算するステップと、
   処理部は、前記式（２０）のＨ_∞フィルタのフィルタ方程式を更新するステップと、
   処理部は、前記各ステップを、時刻ｋを進ませて繰り返し実行するステップと
   を含む請求項７に記載のシステム推定方法。
9. 前記ハイパーＨ_∞フィルタを実行するステップは、
   処理部は、前記フィルタゲインＫ_ｓ，ｋを、ゲイン行列Ｋ_ｋを用いて、次式により求めることを特徴とする請求項１に記載のシステム推定方法。
   

   

   ただし、
   

   

   ここで、
   ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
   ｙ_ｋ：観測信号
   Ｋ_ｓ，ｋ：フィルタゲイン
   Ｈ_ｋ：観測行列
   Θ（ｋ）：Ｊ−ユニタリ行列
   Ｒ_ｅ，ｋ：補助変数
10. 前記ハイパーＨ_∞フィルタを実行するステップは、
    処理部は、Ｋ⁻ _ｋを前記式（６３）を用いて計算するステップと、
    処理部は、初期条件に基づき、フィルタゲインＫ⁻ _ｓ，ｋを前記式（６２）を用いて計算するステップと、
    処理部は、前記式（６１）のＨ_∞フィルタのフィルタ方程式を更新するステップと、
    処理部は、前記各ステップを、時刻ｋを進ませて繰り返し実行するステップと
    を含む請求項９に記載のシステム推定方法。
11. 前記ハイパーＨ_∞フィルタを実行するステップは、
    処理部は、フィルタゲインＫ_ｓ，ｋを、ゲイン行列Ｋ⁻ _ｋを用いて、次式により求める請求項１に記載のシステム推定方法。
    

    

    ただし、
    

    

    なお、上式はＫ⁻ _ｋの代わりにＫ_ｋに関しても整理できる。
    ここで、
    ｙ_ｋ：観測信号
    Ｆ_ｋ：システムのダイナミックス
    Ｈ_ｋ：観測行列
    ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
    Ｋ_ｓ，ｋ：フィルタゲイン；ゲイン行列Ｋ⁻ _ｋから得られる。
    Ｒ_ｅ，ｋ、Ｌ^〜 _ｋ：補助変数
12. 前記ハイパーＨ_∞フィルタを実行するステップは、
    処理部は、予め定められた初期条件に基づき、Ｋ⁻ _ｋ＋１を前記式（２７）を用いて再帰的に計算するステップと、
    処理部は、システムゲインＫ_ｓ，ｋを前記式（２６）を用いて計算するステップと、
    処理部は、存在条件を計算するステップと、
    処理部は、前記存在条件を満たせば、前記式（２５）のＨ_∞フィルタのフィルタ方程式を更新し、時刻ｋを進ませて繰り返し各前記ステップを繰り返し実行するステップと、
    処理部は、前記存在条件を満たさなければ上限値γ_ｆを増加するステップと
    を含む請求項１１に記載のシステム推定方法。
13. さらに、次式により時刻ｋの状態推定値ｘ＾_ｋ｜ｋから出力信号の推定値ｚ^ｖ _ｋ｜ｋを求めるようにした請求項１に記載のシステム推定方法。
    ｚ^ｖ _ｋ｜ｋ＝Ｈ_ｋｘ＾_ｋ｜ｋ
14. 前記Ｈ_∞フィルタ方程式を適用し、状態推定値ｘ＾_ｋ｜ｋを求め、
    擬似エコーを次式のように推定し、
    求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現する請求項１に記載のシステム推定方法。
    

    
15. 次式で表される状態空間モデルに対して、
    ｘ_ｋ＋１＝Ｆ_ｋｘ_ｋ＋Ｇ_ｋｗ_ｋ
    ｙ_ｋ＝Ｈ_ｋｘ_ｋ＋ｖ_ｋ
    ｚ_ｋ＝Ｈ_ｋｘ_ｋ
    ここで、
    ｘ_ｋ：状態ベクトルまたは単に状態
    ｗ_ｋ：システム雑音
    ｖ_ｋ：観測雑音
    ｙ_ｋ：観測信号
    ｚ_ｋ：出力信号
    Ｆ_ｋ：システムのダイナミックス
    Ｇ_ｋ：駆動行列
    評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γ_ｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムであって、
    処理部は、上限値γ_ｆ、フィルタの入力である観測信号ｙ_ｋ、観測行列Ｈ_ｋを含む値を記憶部又は入力部から入力するステップと、
    処理部は、前記上限値γ_ｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
    処理部は、記憶部から初期値又はある時刻の観測行列Ｈ_ｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ_∞フィルタを実行するステップと、
    ｘ＾_ｋ｜ｋ＝Ｆ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}＋Ｋ_ｓ，ｋ（ｙ_ｋ−Ｈ_ｋＦ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}）
    ここで、
    ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
    Ｆ_ｋ：システムのダイナミックス
    Ｋ_ｓ，ｋ：フィルタゲイン
    処理部は、ハイパーＨ_∞フィルタに関する求められた値を記憶部に記憶するステップと、
    処理部は、求められた観測行列Ｈ_ｉ、又は、観測行列Ｈ_ｉとフィルタゲインＫ_ｓ，ｉにより、前記上限値γ_ｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
    処理部は、上限値γ_ｆを小さくしていき前記ハイパーＨ_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
    をコンピュータに実行させるためのシステム推定プログラム。
16. 次式で表される状態空間モデルに対して、
    ｘ_ｋ＋１＝Ｆ_ｋｘ_ｋ＋Ｇ_ｋｗ_ｋ
    ｙ_ｋ＝Ｈ_ｋｘ_ｋ＋ｖ_ｋ
    ｚ_ｋ＝Ｈ_ｋｘ_ｋ
    ここで、
    ｘ_ｋ：状態ベクトルまたは単に状態
    ｗ_ｋ：システム雑音
    ｖ_ｋ：観測雑音
    ｙ_ｋ：観測信号
    ｚ_ｋ：出力信号
    Ｆ_ｋ：システムのダイナミックス
    Ｇ_ｋ：駆動行列
    評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γ_ｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
    処理部は、上限値γ_ｆ、フィルタの入力である観測信号ｙ_ｋ、観測行列Ｈ_ｋを含む値を記憶部又は入力部から入力するステップと、
    処理部は、前記上限値γ_ｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
    処理部は、記憶部から初期値又はある時刻の観測行列Ｈ_ｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ_∞フィルタを実行するステップと、
    ｘ＾_ｋ｜ｋ＝Ｆ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}＋Ｋ_ｓ，ｋ（ｙ_ｋ−Ｈ_ｋＦ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}）
    ここで、
    ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
    Ｆ_ｋ：システムのダイナミックス
    Ｋ_ｓ，ｋ：フィルタゲイン
    処理部は、ハイパーＨ_∞フィルタに関する求められた値を記憶部に記憶するステップと、
    処理部は、求められた観測行列Ｈ_ｉ、又は、観測行列Ｈ_ｉとフィルタゲインＫ_ｓ，ｉにより、前記上限値γ_ｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
    処理部は、上限値γ_ｆを小さくしていき前記ハイパーＨ_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
    をコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
17. 次式で表される状態空間モデルに対して、
    ｘ_ｋ＋１＝Ｆ_ｋｘ_ｋ＋Ｇ_ｋｗ_ｋ
    ｙ_ｋ＝Ｈ_ｋｘ_ｋ＋ｖ_ｋ
    ｚ_ｋ＝Ｈ_ｋｘ_ｋ
    ここで、
    ｘ_ｋ：状態ベクトルまたは単に状態
    ｗ_ｋ：システム雑音
    ｖ_ｋ：観測雑音
    ｙ_ｋ：観測信号
    ｚ_ｋ：出力信号
    Ｆ_ｋ：システムのダイナミックス
    Ｇ_ｋ：駆動行列
    評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γ_ｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
    推定アルゴリズムを実行する処理部と、
    前記処理部により読み取り及び／又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
    を備え、
    前記処理部は、上限値γ_ｆ、フィルタの入力である観測信号ｙ_ｋ、観測行列Ｈ_ｋを含む値を記憶部又は入力部から入力すること、
    前記処理部は、前記上限値γ_ｆに従い、状態空間モデルに関連する忘却係数ρを決定すること、
    前記処理部は、記憶部から初期値又はある時刻の観測行列Ｈ_ｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ_∞フィルタを実行すること、
    ｘ＾_ｋ｜ｋ＝Ｆ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}＋Ｋ_ｓ，ｋ（ｙ_ｋ−Ｈ_ｋＦ_ｋ−１ｘ＾_{ｋ−１｜ｋ−１}）
    ここで、
    ｘ＾_ｋ｜ｋ：観測信号ｙ_０〜ｙ_ｋまでを用いた時刻ｋの状態ｘ_ｋの推定値
    Ｆ_ｋ−１：システムのダイナミックス
    Ｋ_ｓ，ｋ：フィルタゲイン
    前記処理部は、ハイパーＨ_∞フィルタに関する求められた値を記憶部に記憶すること、
    前記処理部は、求められた観測行列Ｈ_ｉ、又は、観測行列Ｈ_ｉとフィルタゲインＫ_{ｓ， ｉ}により、前記上限値γ_ｆ及び前記忘却係数ρに基づく存在条件を計算すること、
    前記処理部は、上限値γ_ｆを小さくしていき前記ハイパーＨ_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること
    を備えた前記システム推定装置。

Images