(57)【要約】忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。まず、処理部は、記憶部又は入力部ら上限値γｆを読み出し又は入力する（Ｓ１０１）。処理部は、式（１５）によって忘却係数ρを決定する（Ｓ１０３）。その後、処理部は、忘却係数ρに基づき、式（１０）〜式（１３）のハイパーＨ∞フィルタを実行する（Ｓ１０５）。処理部１０１は、式（１７）（あるいは、後述の式（１８））の存在条件を計算し（Ｓ１０７）、その存在条件がすべての時刻で満たされれば（Ｓ１０９）、γｆをΔγだけ小さくして同じ処理を繰り返す（Ｓ１１１）。一方、あるγｆで存在条件を満たさなくなったとき（Ｓ１０９）、そのγｆにΔγを加えたものをγｆの最適値γｆｏｐとして出力部に出力及び／又は記憶部に記憶する（Ｓ１１３）。
2005015737.tif 000002
【特許請求の範囲】
【請求項１】
次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部は、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ∞フィルタを実行するステップと、
ｘ＾ｋ｜ｋ＝Ｆｋ−１ｘ＾ｋ−１｜ｋ−１＋Ｋｓ，ｋ（ｙｋ−ＨｋＦｋ−１ｘ＾ｋ−１｜ｋ−１）
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｋｓ，ｋ：フィルタゲイン
処理部は、ハイパーＨ∞フィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項２】
処理部は、前記存在条件を次式に従い計算する請求項１に記載のシステム推定方法。
2005015737.tif 000049
【請求項３】
処理部は、前記存在条件を次式に従い計算する請求項１に記載のシステム推定方法。
2005015737.tif 000050
ここで、
2005015737.tif 000051
【請求項４】
前記忘却係数ρ及び前記上限値γｆは、次式の関係である請求項１に記載のシステム推定方法。
０＜ρ＝１−χ（γｆ）≦１（ただし、χ（γｆ）は、χ（１）＝１、χ（∞）＝０を満たすγｆの単調減衰関数）
【請求項５】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、前記フィルタゲインＫｓ，ｋを、次式により求めることを特徴とする請求項１に記載のシステム推定方法。
2005015737.tif 000052
ここで、
2005015737.tif 000053
なお、式（１６）の右辺はより一般化することもできる。
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｈｋ：観測行列
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応
Ｋｓ，ｋ：フィルタゲイン
ｅｆ，ｉ：フィルタ誤差
Ｒｅ，ｋ：補助変数
【請求項６】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、初期条件に基づき、フィルタゲインＫｓ，ｋを前記式（１２）を用いて計算するステップと、
処理部は、前記式（１１）のＨ∞フィルタのフィルタ方程式を更新するステップと、
処理部は、Σ＾ｋ｜ｋ、Σ＾ｋ＋１｜ｋを前記式（１３）を用いて計算するステップと、
処理部は、前記各ステップを、時刻ｋを進ませて繰り返し実行するステップと
を含む請求項５に記載のシステム推定方法。
【請求項７】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、前記フィルタゲインＫｓ，ｋを、ゲイン行列Ｋｋを用いて、次式により求めることを特徴とする請求項１に記載のシステム推定方法。
2005015737.tif 000054
ただし、
2005015737.tif 000055
2005015737.tif 000056
なお、式（２１）、（２２）において、Ｊ１−１およびＪ１は削除可能である。
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋでを用いた時刻ｋの状態ｘｋの推定値
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
Ｈｋ：観測行列
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応
Θ（ｋ）：Ｊ−ユニタリ行列
Ｒｅ，ｋ：補助変数
【請求項８】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、Ｋｋ、Σ＾ｋ＋１｜ｋ１／２を前記式（２２）を用いて計算するステップと、
処理部は、初期条件に基づき、フィルタゲインＫｓ，ｋを前記式（２１）を用いて計算するステップと、
処理部は、前記式（２０）のＨ∞フィルタのフィルタ方程式を更新するステップと、
処理部は、前記各ステップを、時刻ｋを進ませて繰り返し実行するステップと
を含む請求項７に記載のシステム推定方法。
【請求項９】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、前記フィルタゲインＫｓ，ｋを、ゲイン行列Ｋｋを用いて、次式により求めることを特徴とする請求項１に記載のシステム推定方法。
2005015737.tif 000057
ただし、
2005015737.tif 000058
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
ｙｋ：観測信号
Ｋｓ，ｋ：フィルタゲイン
Ｈｋ：観測行列
Θ（ｋ）：Ｊ−ユニタリ行列
Ｒｅ，ｋ：補助変数
【請求項１０】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、Ｋ−ｋを前記式（６３）を用いて計算するステップと、
処理部は、初期条件に基づき、フィルタゲインＫ−ｓ，ｋを前記式（６２）を用いて計算するステップと、
処理部は、前記式（６１）のＨ∞フィルタのフィルタ方程式を更新するステップと、
処理部は、前記各ステップを、時刻ｋを進ませて繰り返し実行するステップと
を含む請求項９に記載のシステム推定方法。
【請求項１１】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、フィルタゲインＫｓ，ｋを、ゲイン行列Ｋ−ｋを用いて、次式により求める請求項１に記載のシステム推定方法。
2005015737.tif 000059
ただし、
2005015737.tif 000060
なお、上式はＫ−ｋの代わりにＫｋに関しても整理できる。
ここで、
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｈｋ：観測行列
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｋｓ，ｋ：フィルタゲイン；ゲイン行列Ｋ−ｋから得られる。
Ｒｅ，ｋ、Ｌ〜ｋ：補助変数
【請求項１２】
前記ハイパーＨ∞フィルタを実行するステップは、
処理部は、予め定められた初期条件に基づき、Ｋ−ｋ＋１を前記式（２７）を用いて再帰的に計算するステップと、
処理部は、システムゲインＫｓ，ｋを前記式（２６）を用いて計算するステップと、
処理部は、存在条件を計算するステップと、
処理部は、前記存在条件を満たせば、前記式（２５）のＨ∞フィルタのフィルタ方程式を更新し、時刻ｋを進ませて繰り返し各前記ステップを繰り返し実行するステップと、
処理部は、前記存在条件を満たさなければ上限値γｆを増加するステップと
を含む請求項１１に記載のシステム推定方法。
【請求項１３】
さらに、次式により時刻ｋの状態推定値ｘ＾ｋ｜ｋから出力信号の推定値ｚｖｋ｜ｋを求めるようにした請求項１に記載のシステム推定方法。
ｚｖｋ｜ｋ＝Ｈｋｘ＾ｋ｜ｋ
【請求項１４】
前記Ｈ∞フィルタ方程式を適用し、状態推定値ｘ＾ｋ｜ｋを求め、
擬似エコーを次式のように推定し、
求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現する請求項１に記載のシステム推定方法。
2005015737.tif 000061
【請求項１５】
次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムであって、
処理部は、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ∞フィルタを実行するステップと、
ｘ＾ｋ｜ｋ＝Ｆｋ−１ｘ＾ｋ−１｜ｋ−１＋Ｋｓ，ｋ（ｙｋ−ＨｋＦｋ−１ｘ＾ｋ−１｜ｋ−１）
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｆｋ：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
処理部は、ハイパーＨ∞フィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラム。
【請求項１６】
次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部は、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ∞フィルタを実行するステップと、
ｘ＾ｋ｜ｋ＝Ｆｋ−１ｘ＾ｋ−１｜ｋ−１＋Ｋｓ，ｋ（ｙｋ−ＨｋＦｋ−１ｘ＾ｋ−１｜ｋ−１）
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｆｋ：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
処理部は、ハイパーＨ∞フィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１７】
次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び／又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
前記処理部は、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力すること、
前記処理部は、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定すること、
前記処理部は、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ∞フィルタを実行すること、
ｘ＾ｋ｜ｋ＝Ｆｋ−１ｘ＾ｋ−１｜ｋ−１＋Ｋｓ，ｋ（ｙｋ−ＨｋＦｋ−１ｘ＾ｋ−１｜ｋ−１）
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｆｋ−１：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
前記処理部は、ハイパーＨ∞フィルタに関する求められた値を記憶部に記憶すること、
前記処理部は、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算すること、
前記処理部は、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること
を備えた前記システム推定装置。
【発明の詳細な説明】
【技術分野】  本発明は、システム推定方法及びプログラム及び記録媒体、システム推定装置に係り、特に、Ｈ∞評価基準に基づいて開発されたハイパーＨ∞フィルタの高速Ｈ∞フィルタリングアルゴリズムを用いて、状態推定のロバスト化と忘却係数の最適化を同時に実現するシステム推定方法及びプログラム及び記録媒体、システム推定装置に関する。
【背景技術】  一般に、システム推定とは、入出力データに基づいてシステムの入出力関係の数理モデル（伝達関数、あるいはインパルス応答など）のパラメータを推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、非特許文献１．１９９３年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
（基本原理）
図８に、システム推定のための構成図の例を示す（未知システムはＩＩＲ（Ｉｎｆｉｎｉｔｅ  Ｉｍｐｕｌｓｅ  Ｒｅｓｐｏｎｓｅ）フィルタで表現してもよい）。
このシステムは、未知システム１、適応フィルタ２を備える。また、適応フィルタ２は、ＦＩＲディジタルフィルタ３、適応アルゴリズム４を有する。
以下に、未知システム１を同定する出力誤差方式の一例を説明する。ここで、ｕｋは未知システム１の入力、ｄｋは所望信号であるシステムの出力、ｄ＾ｋはフィルタの出力である。（なお、「＾」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。）
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差ｅｋ＝ｄｋ−ｄ＾ｋを最小にするように適応アルゴリズムによってＦＩＲディジタルフィルタ３の係数を調節する。
また、従来、システムのパラメータ（状態）の推定には、誤差共分散行列の更新式（リカッチ方程式）に基づくカルマンフィルタが広く用いられて来た。詳細は、非特許文献２．Ｓ．Ｈａｙｋｉｎ：Ａｄａｐｔｉｖｅ  ｆｉｌｔｅｒ  ｔｈｅｏｒｙ，Ｐｒｅｎｔｉｃｅ−Ｈａｌｌ（１９９６）などに示されている。
以下にカルマンフィルタの基本原理について説明する。
次式のように、状態空間モデルで表された線形システム
2005015737.tif 000003
の状態ｘｋの最小分散推定値ｘ＾ｋ｜ｋは、状態の誤差共分散行列Σ＾ｋ｜ｋ−１を用いて次のように得られる。
2005015737.tif 000004
ただし、
2005015737.tif 000005
ｘｋ：状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。
ｙｋ：観測信号；フィルタの入力となり、既知である。
Ｈｋ：観測行列；既知である。
ｖｋ：観測雑音；未知である。
ρ：忘却係数；一般に試行錯誤で決定される。
Ｋｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾１｜０：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉが用いられる。
また、本発明者は、既に高速Ｈ∞フィルタによるシステム同定アルゴリズムを提案した（特許文献１参照）。これは、システム同定のために新たにＨ∞評価基準を定め、これに基づくハイパーＨ∞フィルタの高速アルゴリズムを開発すると共に、この高速Ｈ∞フィルタリングアルゴリズムに基づく高速時変システム同定方法である。高速Ｈ∞フィルタリングアルゴリズムは、単位時間ステップ当たり計算量Ｏ（Ｎ）で急激に変化する時変システムの追跡が可能である。上限値の極限で高速カルマンフィルタリングアルゴリズムと完全に一致する。このようなシステム同定により時不変および時変システムの高速実時間同定および推定を実現することができる。
なお、システム推定の分野で通常知られる方法として、例えば、非特許文献２、３参照のこと。
（エコーキャンセラへの適用例）
国際電話など長距離電話回線では，信号増幅などの理由から４線式回線が用いられている。一方、加入者回線は比較的短距離なので、２線式回線が使用されている。
図９に通信系とエコーについての説明図を示す。２線式回線と４線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Ｂからの信号（音声）は話者Ａのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は４線式回線に漏れ、増幅された後、再び受信者（話者Ａ）に戻ると云った現象が起こる。これがエコー（ｅｃｈｏ）である。エコーは、伝送距離が長くなるにつれて（遅延時間が長くなるにつれて）影響が大きくなり、著しく通話の品質を劣化させる（パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する）。
図１０に、エコーキャンセラの原理図を示す。
そこで、図示のようにエコーキャンセラ（ｅｃｈｏ  ｃａｎｃｅｌｌｅｒ）を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
エコーパスのインパルス応答の推定は、残留エコーｅｋの平均２乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Ａがらの信号（音声）である。一般に、話者２人が同時に話し始めた（ダブルトーク）ときはインパルス応答の推定を中断する。また、ハイブリッドトランスのインパルス応答長は５０［ｍｓ］程度なので、サンプリング周期を１２５［μｓ］とするとエコーパスのインパルス応答の次数は実際は４００程度となる。
【非特許文献１】  １９９３年電子情報通信学会「ディジタル信号処理ハンドブック」
【非特許文献２】  Ｓ．Ｈａｙｋｉｎ：Ａｄａｐｔｉｖｅ  ｆｉｌｔｅｒ  ｔｈｅｏｒｙ，Ｐｒｅｎｔｉｃｅ−Ｈａｌｌ（１９９６）
【非特許文献３】  Ｂ．Ｈａｓｓｉｂｉ，Ａ．Ｈ．Ｓａｙｅｄ，ａｎｄ  Ｔ．Ｋａｉｌａｔｈ：  ″Ｉｎｄｅｆｉｎｉｔｅ−Ｑｕａｄｒａｔｉｃ  Ｅｓｔｉｍａｔｉｏｎ  ａｎｄ  Ｃｏｎｔｒｏｌ″，ＳＩＡＭ（１９９６）
【特許文献１】  特開２００２−１３５１７１号公報
【発明の開示】
しかしながら、式（１）〜（５）のような従来の忘却係数ρを入れたカルマンフィルタでは、忘却係数ρの値は試行錯誤で決定しなければならず非常に手間が掛かった。さらに、決定された忘却ρ係数の値が果たして最適な値であるかどうか判定する手段も無かった。
また、カルマンフィルタで用いる誤差共分散行列は、本来、零でない任意のベクトルとの２次形式が常に正（以下、「正定」という。）であるがコンピュータにより単精度で計算した場合にはその２次形式が負（以下、「負定」という。）となり、数値的に不安定になることが知られている。また、計算量がＯ（Ｎ２）（あるいはＯ（Ｎ３））であるため、状態ベクトルｘｋの次元Ｎが大きい場合、１ステップ当たりの演算回数が急激に増大し、実時間処理には適さなかった。
本発明は、以上の点に鑑み、忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム推定方法を提供することを目的とする。
本発明は、上記課題を解決するために、新たに考案したＨ∞最適化手法を用いて忘却係数が最適決定可能な状態推定アルゴリズムを導出する。さらに、常に正定であるべき誤差共分散行列の代わりに、その因数行列を更新することによって数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。
本発明の第１の解決手段によると、
次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法及びプログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部は、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ∞フィルタを実行するステップと、
ｘ＾ｋ｜ｋ＝Ｆｋ−１ｘ＾ｋ−１｜ｋ−１＋Ｋｓ，ｋ（ｙｋ−ＨｋＦｋ−１ｘ＾ｋ−１｜ｋ−１）
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｆｋ：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
処理部は、ハイパーＨ∞フィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法、各ステップをコンピュータに実行させるためのシステム推定プログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
また、本発明の第２の解決手段によると、
次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γｆに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び／又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
前記処理部は、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力すること、
前記処理部は、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定すること、
前記処理部は、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーＨ∞フィルタを実行すること、
ｘ＾ｋ｜ｋ＝Ｆｋ−１ｘ＾ｋ−１｜ｋ−１＋Ｋｓ、ｋ（ｙｋ−ＨｋＦｋ−１ｘ＾ｋ−１｜ｋ−１）
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｆｋ：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
前記処理部は、ハイパーＨ∞フィルタに関する求められた値を記憶部に記憶すること、
前記処理部は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算すること、
前記処理部は、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること
を備えた前記システム推定装置が提供される。
本発明の推定方法は忘却係数を最適に決定することが可能であり、かつアルゴリズムは単精度でも安定に動作可能であるため、低コストで高い性能が実現できる。一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。
【図面の簡単な説明】  図１は、本実施の形態に関するハードウェアの構成図である。
図２は、Ｈ∞フィルタのロバスト化と忘却係数ρの最適化についてのフローチャートである。
図３は、図２中のＨ∞フィルタ（Ｓ１０５）のアルゴリズムについてのフローチャートである。
図４は、定理２の平方根アレイアルゴリズムの説明図である。
図５は、定理３の数値的に安定な高速アルゴリズムのフローチャートである。
図６は、インパルス応答｛ｈｉ｝ｉ＝０２３の値を示す図である。
図７は、定理３の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果である。
図８は、システム推定のための構成図である。
図９は、通信系とエコーについての説明図である。
図１０は、エコーキャンセラの原理図である。
【発明を実施するための最良の形態】  以下に、本発明の実施の形態について説明する。
１．記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
ｘｋ：状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｆｋ：システムのダイナミックス；既知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測記号ｙ０〜ｙｋまで用いた時刻ｋ＋１の状態ｘｋ＋１の推定値；フィルタ方程式によって与えられる。
ｘ＾０｜０：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾１｜０：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉが用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
ρ：忘却係数；定理１〜３の場合、γｆが決まればρ＝１−χ（γｆ）より自動的に決定される。
ｅｆ，ｉ：フィルタ誤差
Ｒｅ，ｋ：補助変数
なお、記号の上に付される”＾”、”ｖ”は、推定値の意味である。また、”〜”、”−”、”Ｕ”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、ｘ、ｗ、Ｈ、Ｇ，Ｋ、Ｒ、Σ、等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
２．システム推定のハードウェア及びプログラム
本システム推定方法又はシステム推定装置・システムは、その各手順をコンピュータに実行させるためのシステム推定プログラム、システム推定プログラムを記録したコンピュータ読み取り可能な記録媒体、システム推定プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
図１は、本実施の形態に関するハードウェアの構成図である。
このハードウェアは、中央処理装置（ＣＰＵ）である処理部１０１、入力部１０２、出力部１０３、表示部１０４及び記憶部１０５を有する。また、処理部１０１、入力部１０２、出力部１０３、表示部１０４及び記憶部１０５は、スター又はバス等の適宜の接続手段で接続されている。記憶部１０５は、システム推定される「１．記号の説明」で示した既知のデータが必要に応じて記憶される。また、未知・既知のデータや計算されたハイパーＨ∞フィルタに関するデータ・その他のデータが処理部１０１により、必要に応じて書込み及び／又は読み出しされる。
３．忘却係数が最適決定可能なハイパーＨ∞フィルタ
（定理１）
次式のような状態空間モデルを考える。
2005015737.tif 000006
このような状態空間モデルに対して、次式のようなＨ∞評価基準を提案する。
2005015737.tif 000007
このＨ∞評価基準を満たす状態推定値ｘ＾ｋ｜ｋ（あるいは出力推定値ｚｖｋ｜ｋ）は、次のレベルγｆのハイパーＨ∞フィルタによって与えられる。
2005015737.tif 000008
ここで、
2005015737.tif 000009
である。なお、式（１１）はフィルタ方程式、式（１２）はフィルタゲイン、式（１３）はリカッチ方程式をそれぞれ示す。
また、駆動行列Ｇｋは次のように生成される。
2005015737.tif 000010
また、上述の高速Ｈ∞フィルタの追従能力を向上するためには、上限値γｆは次の存在条件を満たすように出来るだけ小さく設定する。
2005015737.tif 000011
ただし、χ（γｆ）は、χ（１）＝１、χ（∞）＝０を満たすγｆの単調減衰関数である。
定理１の特徴は状態推定のロバスト化と忘却係数ρの最適化を同時に行っている点にある。
図２に、Ｈ∞フィルタのロバスト化と忘却係数ρの最適化についてのフローチャートを示す。ここで、
ブロック「ＥＸＣ＞０」：Ｈ∞フィルタの存在条件、
Δγ：正の実数、である。
まず、処理部１０１は、記憶部１０５又は入力部１０２から上限値γｆを読み出し又は入力する（Ｓ１０１）。この例では、γｆ＞＞１が与えられる。処理部１０１は、式（１５）によって忘却係数ρを決定する（Ｓ１０３）。その後、処理部１０１は、忘却係数ρに基づき、式（１０）〜式（１３）のハイパーＨ∞フィルタを実行する（Ｓ１０５）。処理部１０１は、式（１７）（あるいは、後述の式（１８））の右辺（これをＥＸＣとする）を計算し（Ｓ１０７）、その存在条件がすべての時刻で満たされれば（Ｓ１０９）、γｆをΔγだけ小さくして同じ処理を繰り返す（Ｓ１１１）。一方、あるγｆで存在条件を満たさなくなったとき（Ｓ１０９）、そのγｆにΔγを加えたものをγｆの最適値γｆｏｐとして出力部１０３に出力及び／又は記憶部１０５に記憶する（Ｓ１１３）。なお、この例では、Δγを加えているが、それ以外の予め設定された値を加えてもよい。この最適化のプロセスをγ−イタレーションと呼ぶ。なお、処理部１０１は、Ｈ∞フィルタ計算ステップＳ１０５及び存在条件の計算ステップＳ１０７等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。
ハイパーＨ∞フィルタが存在条件を満たすとき、式（９）の不等式は常に満たされる。よって、式（９）の分母の外乱エネルギーが限定される場合、式（９）の分子の２乗推定誤差の総和は有界となり、ある時刻以降の推定誤差が０になる。これは、γｆをより小さく出来れば、状態ｘｋの変化に推定値ｘ＾ｋ｜ｋが速やかに追従できることを意味する。
ここで、定理１のハイパーＨ∞フィルタのアルゴリズムは通常のＨ∞フィルタのものとは異なることに注意されたい。また、γｆ→∞のとき、ρ＝１、Ｇｋ＝０となり、定理１のＨ∞フィルタのアルゴリズムはカルマンフィルタのアルゴリズムと一致する。
図３に、図２中の（ハイパー）Ｈ∞フィルタ（Ｓ１０５）のアルゴリズムについてのフローチャートを示す。
ハイパーＨ∞フィルタリングアルゴリズムは以下のように要約することができる。
［ステップＳ２０１］処理部１０１は、記憶部１０５から再帰式の初期条件を読み出し、又は、初期条件を入力部１０２から入力し、図示のように定める。なお、Ｌはあらかじめ定められた最大データ数を示す。
［ステップＳ２０３］処理部１０１は、時刻ｋと最大データ数Ｌとを比較する。処理部１０１は、時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、必要に応じて再スタートを行ってもよい。）
［ステップＳ２０５］処理部１０１は、フィルタゲインＫｓ，ｋを式（１２）を用いて計算する。
［ステップＳ２０７］処理部１０１は、式（１１）のハイパーＨ∞フィルタのフィルタ方程式を更新する。
［ステップＳ２０９］処理部１０１は、誤差の共分散行列に対応する項Σ＾ｋ｜ｋ、Σ＾ｋ＋１｜ｋを式（１３）のリカッチ方程式を用いて計算する。
［ステップＳ２１１］時刻ｋを進ませて（ｋ＝ｋ＋１）、ステップＳ２０３に戻り、データがある限り続ける。
なお、処理部１０１は、Ｈ∞フィルタ計算ステップＳ２０５〜Ｓ２０９等の各ステップで求めた適宜の中間値及び最終値、存在条件の値等を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。
（スカラー存在条件）
ところで、式（１７）の存在条件の判定にはＯ（Ｎ２）の計算量が必要であった。しかし、次の条件を用いれば計算量Ｏ（Ｎ）で定理１のＨ∞フィルタの存在性、すなわち式（９）を検証することができる。
系１：スカラー存在条件
次の存在条件を用いれば計算量Ｏ（Ｎ）でハイパーＨ∞フィルタの存在性が判定できる。
2005015737.tif 000012
ここで、
2005015737.tif 000013
ただし、Ｋｓ，ｉは式（１２）で求めたフィルタゲインである。
（証明）
以下に系１の証明を説明する。
２×２の行列Ｒｅ，ｋの特性方程式
2005015737.tif 000014
を解けば、Ｒｅ，ｋの固有値λｉが次のように得られる。
2005015737.tif 000015
もし、
2005015737.tif 000016
であれば、行列Ｒｅ，ｋの２つの固有値の１つは正となり、もう１つは負となり、行列ＲｋとＲｅ，ｋは同じイナーシャをもつ。これより、
2005015737.tif 000017
を用いれば、式（１８）の存在条件が得られる。ここで、ＨｋＫｓ，ｋの計算量はＯ（Ｎ）である。
４．数値的に安定な状態推定アルゴリズム
上述のハイパーＨ∞フィルタは、Σ＾ｋ｜ｋ−１∈ＲＮ×Ｎを更新するため、単位時間ステップ当たりの計算量はＯ（Ｎ２）となる、すなわち、Ｎ２に比例する算術演算が必要となる。ここで、Ｎは状態ベクトルｘｋの次元である。よって、ｘｋの次元が増加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列Σ＾ｋ｜ｋ−１は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度（例：３２ｂｉｔ）でも動作可能な状態推定アルゴリズムの開発が望まれる。
そこで、次に、
Ｒｋ＝Ｒ１／２ｋＪ１Ｒ１／２ｋ、
Ｒｅ，ｋ＝Ｒ１／２ｅ，ｋＪ１Ｒ１／２ｅ，ｋ、
Σ＾ｋ｜ｋ−１＝Σ＾１／２ｋ｜ｋ−１Σ＾１／２ｋ｜ｋ−１
に着目して、数値的に安定化した定理１のＨ∞フィルタ（平方根アレイアルゴリズム）を定理２に示す。ただし、ここでは簡単のためＦｋ＝Ｉとしたが、Ｆｋ≠Ｉの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパーＨ∞フィルタを示す。
（定理２）
2005015737.tif 000018
ただし、
2005015737.tif 000019
なお、式（２１）、（２２）において、Ｊ１−１およびＪ１は削除可能である。
図４に、定理２の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図２に示した定理１のフローチャート中のＨ∞フィルタの計算（Ｓ１０５）で用いることができる。
本推定アルゴリズムは、Σ＾ｋ｜ｋ−１をリカッチ型の更新式で求める代わりに、その因数行列Σ＾１／２ｋ｜ｋ−１∈ＲＮ×Ｎ（Σ＾ｋ｜ｋ−１の平方根行列）をＪ−ユニタリ変換に基づく更新式で求めている。このとき生じる１−１ブロック行列と２−１ブロック行列からフィルタゲインＫｓ，ｋを図示のように求めている。このため、Σ＾ｋ｜ｋ−１＝Σ＾１／２ｋ｜ｋ−１Σ＾１／２ｋ｜ｋ−１＞０となり、Σ＾ｋ｜ｋ−１の正定性は保証され、数値的に安定化できる。なお、定理２のＨ∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ２）のままである。
なお、図４において、Ｊ１−１は削除可能である。
まず、処理部１０１は、式（２２）の左辺の行列式の各要素に含まれる項を記憶部１０５から読み出し又は内部メモリ等から得て、Ｊ−ユニタリ変換を実行する（Ｓ３０１）。処理部１０１は、求めた式（２２）の右辺の行列式の要素からシステムゲインＫｋ、Ｋｓ，ｋを式（２１）に基づき計算する（Ｓ３０３、Ｓ３０５）。処理部１０１は、式（２０）に基づき状態推定値ｘ＾ｋ｜ｋを計算する（Ｓ３０７）。
５．状態推定のための数値的に安定な高速アルゴリズム
上述のように、定理２のＨ∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ２）のままである。そこで、計算量の対策として、Ｈｋ＝Ｈｋ＋１Ψ，Ｈｋ＝［ｕ（ｋ），・・・，ｕ（０），０，・・・，０］のとき、ｘｋ＝［ｘＴｋ，０Ｔ］Ｔの１ステップ予測誤差の共分散行列Σｋ＋１｜ｋが
2005015737.tif 000020
を満たすことを利用して、Σｋ＋１｜ｋの代わりに次元の低いＬｋ（すなわちＬ〜ｋ）を更新する
2005015737.tif 000021
得られる。
（定理３）
2005015737.tif 000022
ただし、
2005015737.tif 000023
なお、定理３の証明は、後述する。
上式は、Ｋ−ｋ（＝Ｐ−１／２Ｋｋ）の代わりにＫｋについても整理することができる。
さらに、次のＪ−ユニタリ行列
2005015737.tif 000024
を用いれば定理４の高速化した状態推定アルゴリズムが得られる。ただし、Ψはシフト行列を表す。
（定理４）
2005015737.tif 000025
ただし、
2005015737.tif 000026
であり、ｄｉａｇ｛・｝は対角行列、Ｒｅ，ｋ＋１（１，１）は行列Ｒｅ，ｋ＋１の１−１成分をそれぞれ表す。また、上式はＫ−ｋの代わりにＫｋに関しても整理できる。
本高速アルゴリズムは、次の因数分解
2005015737.tif 000027
におけるＬ〜ｋ∈Ｒ（Ｎ＋１）×２の更新によってフィルタゲインＫｓ，ｋを求めているので、単位ステップ当たりの計算量はＯ（Ｎ＋１）で済む。ここで、次式に注意されたい。
2005015737.tif 000028
図５に、定理３の数値的に安定な高速アルゴリズムのフローチャートを示す。この高速アルゴリズムは図２のＨ∞フィルタの計算ステップ（Ｓ１０５）に組み込まれ、γ−イタレーションによって最適化される。よって、存在条件が満たされる間はγｆは除々に減少されるが、満たされなくなった時点で、図示のようにγｆは増加される。
Ｈ∞フィルタリングアルゴリズムは以下のように要約することができる。
［ステップＳ４０１］処理部１０１は、再帰式の初期条件を図示のように定める。なお、Ｌは最大データ数を示す。
［ステップＳ４０３］処理部１０１は、時刻ｋと最大データ数Ｌとを比較する。処理部１０１は、時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、再スタートする。）
［ステップＳ４０５］処理部１０１は、フィルタゲインに対応する項Ｋｋ＋１を式（２７）、（３１）を用いて再帰的に計算する。
［ステップＳ４０６］処理部１０１は、Ｒｅ，ｋ＋１を式（２９）を用いて再帰的に計算する。
［ステップＳ４０７］処理部１０１は、さらにＫｓ，ｋを式（２６）、（３１）を用いて計算する。
［ステップＳ４０９］処理部１０１は、ここで、存在条件ＥＸＣ＞０を判定し、存在条件を満たせばステップＳ４１１に進む。
［ステップＳ４１３］一方、処理部１０１は、ステップＳ４０９で存在条件を満たさなければγｆを増加し、ステップＳ４０１に戻る。
［ステップＳ４１１］処理部１０１は、式（２５）のＨ∞フィルタのフィルタ方程式を更新する。
［ステップＳ４１５］処理部１０１は、Ｒｒ，ｋ＋１を式（３０）を用いて再帰的に計算する。また、処理部１０１は、Ｌ〜ｋ＋１を式（２８）、（３１）を用いて再帰的に計算する。
［ステップＳ４１９］処理部１０１は、時刻ｋを進ませて（ｋ＝ｋ＋１）、ステップＳ４０３に戻り、データがある限り続ける。
なお、処理部１０１は、Ｈ∞フィルタ計算ステップＳ４０５〜Ｓ４１５及び存在条件の計算ステップＳ４０９等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。
６．エコーキャンセラ
つぎに、エコーキャンセリング問題の数理モデルを作成する。
まず、受信信号｛ｕｋ｝がエコーパスへの入力信号となることを考慮すれば、エコーパスの（時変）インパルス応答｛ｈｉ［ｋ］｝により、エコー｛ｄｋ｝の観測値｛ｙｋ｝は次式で表される。
2005015737.tif 000029
ここで、ｕｋ，ｙｋはそれぞれ時刻ｔｋ（＝ｋＴ；Ｔは標本化周期）における受信信号とエコーを表し、ｖｋは時刻ｔｋにおける平均値０の回線雑音とし、ｈｉ［ｋ］，ｉ＝０，・・・，Ｎ−１は、時変インパルス応答であり、そのタップ数Ｎは既知とする。このとき、インパルス応答の推定値｛ｈ＾ｉ［ｋ］｝が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
2005015737.tif 000030
これをエコーから差し引けば（ｙｋ−ｄ＾ｋ≒０）、エコーをキャンセルすることができる。ただし、ｋ−ｉ＜０のときｕｋ−１＝０とする。
以上より、問題は直接観測可能な受信信号｛ｕｋ｝とエコー｛ｙｋ｝からエコーパスのインパルス応答｛ｈｉ［ｋ］｝を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにＨ∞フィルタを適用するには、まず式（３２）を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、問題がインパルス応答｛ｈｉ［ｋ］｝を推定することであるから、｛ｈｉ［ｋ］｝を状態変数ｘｋとし、ｗｋ程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
2005015737.tif 000031
ただし、
2005015737.tif 000032
このような状態空間モデルに対するハイパーおよび高速Ｈ∞フィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
７．インパルス応答に対する評価
（動作の確認）
エコーパスのインパルス応答が時間的に不変であり（ｈｉ［ｋ］＝ｈｉ）、かつそのタップ数Ｎが４８である場合について、シミュレーションを用いて、本高速アルゴリズムの動作を確認する。
2005015737.tif 000033
なお、図６は、ここでのインパルス応答｛ｈｉ｝の値を示す図である。
ここで、インパルス応答｛ｈｉ｝ｉ＝０２３は、図示の値を採用し、その他｛ｈｉ｝ｉ＝２４４７は０とする。また、υｋは平均値０、分散σｖ２＝１．０×１０−６の定常なガウス白色雑音とし、標本化周期Ｔを便宜上１．０とする。
また、受信信号｛ｕｋ｝は次のように２次のＡＲモデルで近似する。
2005015737.tif 000034
ただし、α１＝０．７，α２＝０．１とし、ｗｋ′は平均値０、分散σｗ′２＝０．０４の定常なガウス白色雑音とする。
（インパルス応答の推定結果）
図７に、定理３の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図７（ｂ）の縦軸は、
√｛Σｉ＝０４７（ｈｉ−ｘ＾ｋ（ｉ＋１））２｝
を表す。
これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし、ρ＝１−χ（γｆ）、χ（γｆ）＝γｆ−２、ｘ＾０｜０＝０、Σ＾１｜０＝２０Ｉとし、計算は倍精度で行った。また、存在条件を確認しつつ、γｆ＝５．５と設定とした。
８．定理の証明
８−１．定理２の証明
次の関係式
2005015737.tif 000035
が成り立つとき、両辺の２×２ブロック行列の各項を比較すれば次式が得られる。
2005015737.tif 000036
これは定理１のＦｋ＝Ｉのときの式（１３）のリカッチ方程式と一致する。ただし、
2005015737.tif 000037
一方、ＡＪＡＴ＝ＢＪＢＴが成り立つとき、ＢはＪ−ユニタリ行列Θ（ｋ）を用いてＢ＝ＡΘ（ｋ）と表すことができる。よって、式（４０）より定理１のリカッチ方程式は次式と等価である。
2005015737.tif 000038
なお、式（４０）、（４５）において、Ｊ１−１は削除可能である。
８−２．定理３の証明
次のようにブロック三角化するＪ−ユニタリ行列Θ（ｋ）が存在すると仮定する。
2005015737.tif 000039
のように決定することができる。ただし、Ｓは対角成分が１または−１をとる対角行列とする。
（１，１）−ブロック行列
2005015737.tif 000040
（２，１）−ブロック行列
2005015737.tif 000041
（２，２）−ブロック行列
2005015737.tif 000042
これより、
2005015737.tif 000043
８−３．定理４の証明
観測行列Ｈｋがシフト特性をもち、かつ
2005015737.tif 000044
のとき、定理２と同様な方法によって次の関係式が得られる。
2005015737.tif 000045
ただし、
2005015737.tif 000046
となるようにＲｒ，ｋ＋１を決定する。次に、式（４６）の３行目にＲｒ，ｋ＋１の更新式を新たに追加すれば、最終的に次式が得られる。
2005015737.tif 000047
この両辺の３×２ブロック行列の各項の対応から次のゲイン行列Ｋ−ｋの更新式が得られる。
2005015737.tif 000048
【産業上の利用可能性】  一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができる。
【図１】
2005015737.tif 000062
【図２】
2005015737.tif 000063
【図３】
2005015737.tif 000064
【図４】
2005015737.tif 000065
【図５】
2005015737.tif 000066
【図６】
2005015737.tif 000067
【図７】
2005015737.tif 000068
【図８】
2005015737.tif 000069
【図９】
2005015737.tif 000070
【図１０】
2005015737.tif 000071
【手続補正書】
【提出日】平成17年10月27日(2005.10.27)
【手続補正００１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【書類名】特許請求の範囲
【請求項１】削    除
【請求項２】処理部は、前記存在条件を次式に従い計算する請求項７に記載のシステム推定方法。
2005015737.tif 000072
【請求項３】処理部は、前記存在条件を次式に従い計算する請求項７又は９又は１１に記載のシステム推定方法。
2005015737.tif 000073
ここで、
2005015737.tif 000074
ただし、前記忘却係数ρ及び前記上限値γｆは、次式の関係である。
０＜ρ＝１−χ（γｆ）≦１（ただし、χ（γｆ）は、χ（１）＝１、χ（∞）＝０を満たすγｆの単調減衰関数）
【請求項４】（削除）
【請求項５】削    除
【請求項６】削    除
【請求項７】次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として、システム雑音ｗｋ及び観測雑音ｖｋを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γｆに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、前記フィルタゲインＫｓ，ｋを、前記忘却係数ρとゲイン行列Ｋｋ用いて、次式（２０）〜（２２）により求める、又は、次式（２０）と、式（２１）及び（２２）においてＪ１−１およびＪ１を削除した式により求めるハイパーＨ∞フィルタを実行するステップと、
2005015737.tif 000075
ただし、
2005015737.tif 000076
2005015737.tif 000077
位行列である。また、Ｋｋ（：，１）は行列Ｋｋの１列目の列ベクトルを表す。
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋでを用いた時刻ｋの状態ｘｋの推定値
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｋｓ，ｋ：フィルタゲイン
Ｈｋ：観測行列
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応
Θ（ｋ）：Ｊ−ユニタリ行列Ｒｅ，ｋ：補助変数
処理部が、ハイパーＨ∞フィルタによって求められた状態ｘｋの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項８】前記ハイパーＨ∞フィルタを実行するステップは、
処理部が、Σ＾ｋ＋１｜ｋ１／２を前記式（２２）を用いて計算するステップと、
処理部が、Σ＾ｋ｜ｋの初期条件とＣｋの初期条件のもとで、フィルタゲインＫｓ，ｋを前記式（２１）を用いて計算するステップと、
処理部が、前記式（２０）のＨ∞フィルタのフィルタ方程式を更新するステップと、
処理部が、前記式（２２）を用いて計算するステップと、前記式（２１）を用いて計算するステップと、前記更新するステップとを、時刻ｋを進ませて繰り返し実行するステップと
を含む請求項７に記載のシステム推定方法。
【請求項９】次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として、システム雑音ｗｋ及び観測雑音ｖｋ含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γｆに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋ含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Ｈｋ含む値を読み取り、前記フィルタゲインＫｓ，ｋを、前記忘却係数ρとゲイン行列Ｋｋを用いて、次式により求めるハイパーＨ∞フィルタを実行するステップと、
2005015737.tif 000078
2005015737.tif 000079
ただし、
2005015737.tif 000080
ここで、
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
ｙｋ：観測信号
Ｋｓ，ｋ：フィルタゲイン
Ｈｋ：観測行列
Θ（ｋ）：Ｊ−ユニタリ行列Ｒｅ，ｋ：補助変数
処理部が、ハイパーＨ∞フィルタによって求められた状態ｘｋの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項１０】前記ハイパーＨ∞フィルタを実行するステップは、
処理部が、Ｒｅ，ｋ＋１、Ｒｒ，ｋ＋１及びＬ〜ｋ＋１の初期条件のもとで、Ｋ−ｋを前記式（６３）を用いて計算するステップと、
処理部が、フィルタゲインＫｓ，ｋを前記式（６２）を用いて計算するステップと、
処理部が、前記式（６１）のＨ∞フィルタのフィルタ方程式を更新するステップと、
処理部は、前記式（６３）を用いて計算するステップと、前記式（６２）を用いて計算するステップと、前記更新するステップを、時刻ｋを進ませて繰り返し実行するステップと
を含む請求項９に記載のシステム推定方法。
【請求項１１】次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として、システム雑音ｗｋ及び観測雑音ｖｋ含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γｆに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステ厶推定方法であって、
処理部が、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋ含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Ｈｋ含む値を読み取り、フィルタゲインＫｓ，ｋを、前記忘却係数ρとゲイン行列Ｋ−ｋ用いて、次式により求めるハイパーＨ∞フィルタを実行するステップと、
2005015737.tif 000081
ただし、
2005015737.tif 000082
ここで、
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｈｋ：観測行列
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋでを用いた時刻ｋの状態ｘｋの推定値
Ｋｓ，ｋ：フィルタゲイン；ゲイン行列Ｋ−ｋから得られる。
Ｒｅ，ｋ、Ｌ〜ｋ：補助変数
処理部が、ハイパーＨ∞フィルタによって求められた状態ｘｋの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項１２】削    除
【請求項１３】さらに、次式により時刻ｋの状態推定値ｘ＾ｋ｜ｋから出力信号の推定値ｚｖｋ｜ｋを求めるようにした請求項７又は９又は１１に記載のシステム推定方法。
ｚｖｋ｜ｋ＝Ｈｋｘ＾ｋ｜ｋ
【請求項１４】前記Ｈ∞フィルタ方程式を適用し、状態推定値ｘ＾ｋ｜ｋ＝［ｈ＾１［ｋ］，・・・，ｈ＾Ｎ［ｋ］］Ｔを求め、
擬似エコーを次式のように推定し、
求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現する請求項７又は９又は１１に記載のシステム推定方法。
2005015737.tif 000083
【請求項１５】次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋ＧｋＷｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として、システム雑音ｗｋ及び観測雑音ｖｋ含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γｆに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムであって、
処理部が、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、フィルタゲインＫｓ，ｋを、前記忘却係数ρとゲイン行列Ｋ−ｋ諦用いて、次式により求めるハイパーＨ∞フィルタを実行するステップと、
2005015737.tif 000084
ただし、
2005015737.tif 000085
ここで、
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｈｋ：観測行列
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｋｓ，ｋ：フィルタゲイン；ゲイン行列Ｋ−ｋから得られる。
Ｒｅ，ｋ、Ｌ〜ｋ：補助変数
処理部が、ハイパーＨ∞フィルタによって求められた状態ｘｋの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラム。
【請求項１６】次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として、システム雑音ｗｋ及び観測雑音ｖｋを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γｆに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Ｈｋを含む値を読み取り、フィルタゲインＫｓ，ｋを、前記忘却係数ρとゲイン行列Ｋ−ｋ用いて、次式により求めるハイパーＨ∞フィルタを実行するステップと、
2005015737.tif 000086
ただし、
2005015737.tif 000087
ここで、
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｈｋ：観測行列
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｋｓ，ｋ：フィルタゲイン；ゲイン行列Ｋ−ｋから得られる。
Ｒｅ，ｋ、Ｌ〜ｋ：補助変数
処理部が、ハイパーＨ∞フィルタによって求められた状態ｘｋの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γｆを小さくしていき前記ハイパー∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１７】次式で表される状態空間モデルに対して、
ｘｋ＋１＝Ｆｋｘｋ＋Ｇｋｗｋ
ｙｋ＝Ｈｋｘｋ＋ｖｋ
ｚｋ＝Ｈｋｘｋ
ここで、
ｘｋ：状態ベクトルまたは単に状態
ｗｋ：システム雑音
ｖｋ：観測雑音
ｙｋ：観測信号
ｚｋ：出力信号
Ｆｋ：システムのダイナミックス
Ｇｋ：駆動行列
評価基準として、システム雑音ｗｋ及び観測雑音ｖｋを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γｆに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び／又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
処理部が、上限値γｆ、フィルタの入力である観測信号ｙｋ、観測行列Ｈｋ含む値を記憶部又は入力部から入力する手段と、
処理部が、前記上限値γｆに従い、状態空間モデルに関連する忘却係数ρを決定する手段と、
処理部が、記憶部から初期値又はある時刻の観測行列Ｈｋ含む値を読み取り、フィルタゲインＫｓ，ｋを、前記忘却係数ρとゲイン行列Ｋ−ｋを用いて、次式により求めるハイパーＨ∞フィルタを実行する手段と、
2005015737.tif 000088
ただし、
2005015737.tif 000089
ここで、
ｙｋ：観測信号
Ｆｋ：システムのダイナミックス
Ｈｋ：観測行列
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値
Ｋｓ，ｋ：フィルタゲイン；ゲイン行列Ｋ−ｋから得られる。
Ｒｅ，ｋ、Ｌ〜ｋ：補助変数
処理部が、ハイパーＨ∞フィルタによって求められた状態ｘｋの推定値を記憶部に記憶する手段と、
処理部が、求められた観測行列Ｈｉ、又は、観測行列ＨｉとフィルタゲインＫｓ，ｉにより、前記上限値γｆ及び前記忘却係数ρに基づく存在条件を計算する手段と、
処理部が、上限値γｆを小さくしていき前記ハイパーＨ∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶する手段と、
を備えた前記システム推定装置。


【手続補正００２】
【補正対象書類名】明細書
【補正対象項目名】００１０
【補正方法】変更
【補正の内容】
【００１０】
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｆｋ：システムのダイナミックス；既知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋの推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測記号ｙ０〜ｙｋで用いた時刻ｋ＋１の状態ｘｋ＋１の推定値；フィルタ方程式によって与えられる。
ｘ＾０｜０：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾１｜０：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉが用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
ρ：忘却係数；定理１〜３の場合、γｆが決まればρ＝１−χ（γｆ）より自動的に決定される。
ｅｆ，ｉ：フィルタ誤差
Ｒｅ，ｋ：補助変数
なお、記号の上に付される”＾”、”ｖ”は、推定値の意味である。また、”〜”、”−”、”Ｕ”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、ｘ、ｗ等はベクトル、Ｈ、Ｇ，Ｋ、Ｒ、Σ、等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
【手続補正００３】
【補正対象書類名】明細書
【補正対象項目名】００１６
【補正方法】変更
【補正の内容】
【００１６】
る。
４．数値的に安定な状態推定アルゴリズム
上述のハイパーＨ∞フィルタは、Σ＾ｋ｜ｋ−１∈ＲＮ×Ｎを更新するため、単位時間ステップ当たりの計算量はＯ（Ｎ２）となる、すなわち、Ｎ２に比例する算術演算が必要となる。ここで、Ｎは状態ベクトルｘｋの次元である。よって、ｘｋの次元が増加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列Σ＾ｋ｜ｋ−１は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度（例：３２ｂｉｔ）でも動作可能な状態推定アルゴリズムの開発が望まれる。
そこで、次に、
Ｒｋ＝Ｒ１／２ｋＪ１ＲＴ／２ｋ、
Ｒｅ，ｋ＝Ｒ１／２ｅ，ｋＪ１ＲＴ／２ｅ，ｋ、
Σ＾ｋ｜ｋ−１＝Σ＾１／２ｋ｜ｋ−１Σ＾Ｔ／２ｋ｜ｋ−１
に着目して、数値的に安定化した定理１のＨ∞フィルタ（平方根アレイアルゴリズム）を定理２に示す。ただし、ここでは簡単のためＦｋ＝Ｉとしたが、Ｆｋ≠Ｉの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパーＨ∞フィルタを示す。
【手続補正００４】
【補正対象書類名】明細書
【補正対象項目名】００１７
【補正方法】変更
【補正の内容】
【００１７】
（定理２）
2005015737.tif 000090
ただし、
2005015737.tif 000091
2005015737.tif 000092
位行列である。また、Ｋｋ（：，１）は行列Ｋｋの１列目の列ベクトルを表す。
なお、式（２１）、（２２）において、Ｊ１−１およびＪ１は削除可能である。
図４に、定理２の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図２に示した定理１のフローチャート中のＨ∞フィルタの計算（Ｓ１０５）で用いることができる。
本推定アルゴリズムは、Σ＾ｋ｜ｋ−１をリカッチ型の更新式で求める代わりに、その因数行列Σ＾１／２ｋ｜ｋ−１∈ＲＮ×Ｎ（Σ＾ｋ｜ｋ−１の平方根行列）をＪ−ユニタリ変換に基づく更新式で求めている。このとき生じる１−１ブロック行列と２−１ブロック行列からフィルタゲインＫｓ，ｋを図示のように求めている。このため、Σ＾ｋ｜ｋ−１＝Σ＾１／２ｋ｜ｋ−１Σ＾Ｔ／２ｋ｜ｋ−１＞０となり、Σ＾ｋ｜ｋ−１の正定性は保証され、数値的に安定化できる。なお、定理２のＨ∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ２）のままである。
なお、図４において、Ｊ１−１は削除可能である。
まず、処理部１０１は、式（２２）の左辺の行列式の各要素に含まれる項を記憶部１０５から読み出し又は内部メモリ等から得て、Ｊ−ユニタリ変換を実行する（Ｓ３０１）。処理部１０１は、求めた式（２２）の右辺の行列式の要素からシステムゲインＫｋ、Ｋｓ，ｋを
【手続補正００５】
【補正対象書類名】明細書
【補正対象項目名】００１８
【補正方法】変更
【補正の内容】
【００１８】
式（２１）に基づき計算する（Ｓ３０３、Ｓ３０５）。処理部１９１は、式（２０）に基づき状態推定値ｘ＾ｋ｜ｋを計算する（Ｓ３０７）。
５．状態推定のための数値的に安定な高速アルゴリズム
上述のように、定理２のＨ∞フィルタの単位ステップ当たりの計算量はＯ（Ｎ２）のままである。そこで、計算量の対策として、Ｈｋ＝Ｈｋ＋１Ψ，Ｈｋ＝［ｕ（ｋ），・・・，ｕ（０），０，・・・，０］のとき、ｘｋ＝［ｘＴｋ，０Ｔ］Ｔの１ステップ予測誤差の共分散行列Σｋ＋１｜ｋが
2005015737.tif 000093
を満たすことを利用して、Σｋ＋１｜ｋ代わりに次元の低いＬｋ（すなわちＬ〜ｋ）を更新する
2005015737.tif 000094
３が得られる。
（定理３）
2005015737.tif 000095
2005015737.tif 000096
ただし、
2005015737.tif 000097
なお、定理３の証明は、後述する。
上式は、Ｋ−ｋ（＝ρ−１／２Ｋｋ）の代わりにＫｋについても整理することができる。
【手続補正００６】
【補正対象書類名】明細書
【補正対象項目名】００１９
【補正方法】変更
【補正の内容】
【００１９】
さらに、次のＪ−ユニタリ行列
2005015737.tif 000098
を用いれば定理４の高速化した状態推定アルゴリズムが得られる。ただし、Ψはシフト行列を表す。
（定理４）
2005015737.tif 000099
ただし、
2005015737.tif 000100
であり、ｄｉａｇ｛・｝は対角行列、Ｒｅ，ｋ＋１（１，１）は行列Ｒｅ，ｋ＋１の１−１成分をそれぞれ表す。また、上式はＫ−ｋの代わりにＫｋに関しても整理できる。
本高速アルゴリズムは、次の因数分解
2005015737.tif 000101
におけるＬ〜ｋ∈Ｒ（Ｎ＋１）×２の更新によってフィルタゲインＫｓ，ｋを求めているので、単位ス
【手続補正００７】
【補正対象書類名】明細書
【補正対象項目名】００２０
【補正方法】変更
【補正の内容】
【００２０】
テップ当たりの計算量はＯ（Ｎ＋１）で済む。ここで、次式に注意されたい。
2005015737.tif 000102
図５に、定理４の数値的に安定な高速アルゴリズムのフローチャートの一例を示す。この高速アルゴリズムは図２のＨ∞フィルタの計算ステップ（Ｓ１０５）に組み込まれ、γ−イタレーションによって最適化される。よって、存在条件が満たされる間はγｆは除々に減少されるが、満たされなくなった時点で、図示のようにγｆは増加される。
Ｈ∞フィルタリングアルゴリズムは以下のように要約することができる。
［ステップＳ４０１］処理部１０１は、再帰式の初期条件を図示のように定める。なお、Ｌは最大データ数を示す。
［ステップＳ４０３］処理部１０１は、時刻ｋと最大データ数Ｌとを比較する。処理部１０１は、時刻ｋが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、再スタートする。）
［ステップＳ４０５］処理部１０１は、フィルタゲインに対応する項Ｋｋ＋１を式（２７）、（３１）を用いて再帰的に計算する。
［ステップＳ４０６］処理部１０１は、Ｒｅ，ｋ＋１を式（２９）を用いて再帰的に計算する。
［ステップＳ４０７］処理部１０１は、さらにＫｓ，ｋを式（２６）、（３１）を用いて計算する。
［ステップＳ４０９］処理部１０１は、ここで、存在条件ＥＸＣ＞０を判定し、存在条件を満たせばステップＳ４１１に進む。
［ステップＳ４１３］一方、処理部１０１は、ステップＳ４０９で存在条件を満たさなければγｆを増加し、ステップＳ４０１に戻る。
［ステップＳ４１１］処理部１０１は、式（２５）のＨ∞フィルタのフィルタ方程式を更新する。
［ステップＳ４１５］処理部１０１は、Ｒｒ．ｋ＋１を式（３０）を用いて再帰的に計算する。ま
【手続補正００８】
【補正対象書類名】明細書
【補正対象項目名】００２３
【補正方法】変更
【補正の内容】
【００２３】
周期Ｔを便宜上１．０とする。
また、受信信号｛ｕｋ｝は次のように２次のＡＲモデルで近似する。
2005015737.tif 000103
ただし、α１＝０．７，α２＝０．１とし、ｗｋ′は平均値０、分散σｗ′２＝０．０４の定常なガウス白色雑音とする。
（インパルス応答の推定結果）
図７に、定理４の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図７（ｂ）の縦軸は、
√｛Σｉ＝０４７（ｈｉ−ｘ＾ｋ（ｉ＋１））２｝
を表す。
これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし、ρ＝１−χ（γｆ）、χ（γｆ）＝γｆ−２、ｘ＾０｜０＝０、Σ＾１｜０＝２０Ｉとし、計算は倍精度で行った。また、存在条件を確認しつつ、γｆ＝５．５と設定とした。
８．定理の証明
８−１．定理２の証明
次の関係式
2005015737.tif 000104
が成り立つとき、両辺の２×２ブロック行列の各項を比較すれば次式が得られる。
【手続補正００９】
【補正対象書類名】明細書
【補正対象項目名】００２５
【補正方法】変更
【補正の内容】
【００２５】
（１，１）−ブロック行列
2005015737.tif 000105
（２，１）−ブロック行列
2005015737.tif 000106
【手続補正０１０】
【補正対象書類名】明細書
【補正対象項目名】００２６
【補正方法】変更
【補正の内容】
【００２６】
（２，２）−ブロック行列
2005015737.tif 000107
これより、
2005015737.tif 000108
８−３．定理４の証明
観測行列Ｈｋがシフト特性をもち、かつ
2005015737.tif 000109
のとき、定理２と同様な方法によって次の関係式が得られる。
2005015737.tif 000110
ただし、
2005015737.tif 000111
となるようにＲｒ，ｋ＋１を決定する。次に、式（４６）の３行目にＲｒ，ｋ＋１の更新式を新たに追
【手続補正０１１】
【補正対象書類名】図面
【補正対象項目名】図５
【補正方法】変更
【補正の内容】
【図５】
2005015737.tif 000112
【国際調査報告】
2005015737.tif 000113
2005015737.tif 000114
2005015737.tif 000115
2005015737.tif 000116