(57)【要約】【課題】  大規模な音響系または通信系を高速かつ高精度に同定すること。【解決手段】  フィルタ処理では、処理部１０１は、フィルタの再帰変数を初期化し（Ｓ１０１）、時刻ｋがＬ以下であれば（Ｓ１０３）、拡大観測行列ＣＵｋを入力信号とシフト特性を利用して更新する。処理部１０１は、フィルタを表す式に従いＣｏｒｅ  Ａｌｇｏｒｉｔｈｍを実行し（Ｓ１０７）、計算された状態推定値ｘ＾ｋ｜ｋ（あるいは出力推定値ｚ＾ｋ｜ｋ）を記憶部１０５に記憶、及び／又は出力部１０３若しくはＩ／Ｆ１０６から出力する。処理部１０１は、時刻ｋを進ませて、ステップＳ１０３に戻り、データがある限り又は所定回数繰り返し続ける。【選択図】図２
2019205049.tif 000002
【特許請求の範囲】
【請求項１】
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜（６２）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置。
【数１】
2019205049.tif 000056
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【請求項２】
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置。
【数２】
2019205049.tif 000057
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項３】
請求項２に記載されたシステム同定装置において、
式（６０）を用い、式（５８）、式（５５）、式（４９）及び式（５０）について、式（５８）、式（５５）、式（４９）及び式（５０）を変形した以下の式を用いることを特徴とするシステム同定装置。
【数３】
2019205049.tif 000058
ここで、μはステップサイズ、μ＞０の実数である。
【請求項４】
請求項２に記載されたシステム同定装置において、
前記フィルタは、γｆ＝１．０、ρ＝１．０、Ｓｋ＝Ｓｋ−１＝Ｓ−１＝１／ε０とすることにより、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）ように縮退した式を用いることを特徴とする前記システム同定装置。
【数４】
2019205049.tif 000059
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項５】
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置。
【数５】
2019205049.tif 000060
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項６】
請求項５に記載のシステム同定装置において、
式（５８）、式（５５）、式（４９）及び式（５０）について、式（５８）、式（５５）、式（４９）及び式（５０）を変形した以下の式を用いることを特徴とするシステム同定装置。
【数６】
2019205049.tif 000061
ここで、μはステップサイズ、μ＞０の実数である。
【請求項７】
請求項１乃至３のいずれかに記載のシステム同定装置において、
前記フィルタを実行する処理部
を備え、
前記処理部は、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化する第１ステップと、
前記処理部は、時刻ｋが予め定められたデータ長Ｌを超えているか判断し、超えていればフィルタゲインを初期化してリスタートする、又は、処理を終了する、第２ステップと、
前記処理部は、時刻ｋがデータ長Ｌ以下であれば、入力信号ｕｋを入力部若しくはインタフェース部から入力し又は記憶部から読み出し、拡大観測行列ＣＵｋを入力信号ｕｋとシフト特性ＣＵｋ＝［ｃｋ，Ｃｋ−１］＝［Ｃｋ，ｃｋ−Ｎ］、を利用して更新する第３ステップと、
前記処理部は、フィルタを与える式を式番号が式（６２）から式（４９）の降順に、式（４９）の右辺第２項を緩和係数μ（ここで、μは、μ＞０の実数）で計算する処理を含み、単位時間あたりのフィルタ・アルゴリズムを実行する第４ステップと、
前記処理部は、式（４９）を含むフィルタ方程式を更新し、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力する第５ステップと、
前記処理部は、時刻ｋを進ませて、前記第１ステップ乃至前記第４ステップの処理を繰り返す第６ステップと、
を含むことを特徴とするシステム同定装置。
【請求項８】
請求項４又は５又は６に記載のシステム同定装置において、
前記フィルタを実行する処理部
を備え、
前記処理部は、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化する第１ステップと、
前記処理部は、時刻ｋが予め定められたデータ長Ｌを超えているか判断し、超えていればフィルタゲインを初期化してリスタートする、又は、処理を終了する、第２ステップと、
前記処理部は、時刻ｋがデータ長Ｌ以下であれば、入力信号ｕｋを入力部若しくはインタフェース部から入力し又は記憶部から読み出し、拡大観測行列ＣＵｋを入力信号ｕｋとシフト特性ＣＵｋ＝［ｃｋ，Ｃｋ−１］＝［Ｃｋ，ｃｋ−Ｎ］、を利用して更新する第３ステップと、
前記処理部は、フィルタを与える式を式番号が式（５８）から式（４９）の降順に、式（４９）をステップサイズμ（ここで、μは、μ＞０の実数）で計算する処理を含み、単位ステップサイズあたりのフィルタ・アルゴリズムを実行する第４ステップと、
前記処理部は、式（４９）を含むフィルタ方程式を更新し、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力する第５ステップと、
前記処理部は、時刻ｋを進ませて、前記第１ステップ乃至前記第４ステップの処理を繰り返す第６ステップと、
を含むことを特徴とするシステム同定装置。
【請求項９】
請求項１乃至８のいずれかに記載のシステム同定装置において、
前記評価基準は、前記状態空間モデルに対して、次式（１４）を満たすことを特徴とするシステム同定装置。
【数７】
2019205049.tif 000062
ここで、
ｅｆ，ｉ：フィルタ誤差；ｅｆ，ｉ＝ｚｖｉ｜ｉ−Ｈｉｘｉ，ｚｖｉ｜ｉは出力推定値Ｈｉｘ＾i｜ｉである。
Σ０−１：状態の不確定性を表す重み行列の逆行列；Σ０は既知である。
ρ：忘却係数；γｆが決まればρ＝１−χ（γｆ）より決定される。γｆと切り離して独立に決定することも可能である。
γｆ：アッテネーションレベル；設計の際に与えられる。
【請求項１０】
請求項１乃至９のいずれかに記載のシステム同定装置において、
前記フィルタの処理部は、前記第３ステップにおいて、処理の実行後又は実行前に、次式のスカラー存在条件の判定ステップを実行し、式（３０）を満たさなければ処理を終了する又はリスタートすること
【数８】
2019205049.tif 000063
【請求項１１】
請求項１乃至１０のいずれかに記載のシステム同定装置において、
前記フィルタの処理部は、前記第５ステップにおいて、さらに次式（４８）により出力推定値ｚｖｋ｜ｋ，ｚｖｋ｜ｋ−１を計算し、記憶部に記憶及び／又は出力部若しくはインタフェース部から出力することを特徴とするシステム同定装置。
ｚ＾ｋ｜ｋ＝Ｈｋｘ＾ｋ｜ｋ                                  (48)
ここで、
ｚ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの出力ｚｋの推定値Ｈｋｘ＾ｋ｜ｋ；フィルタ方程式によって与えられる。
ｚ＾ｋ｜ｋ−１：観測信号ｙ０〜ｙｋ−１までを用いた時刻ｋの出力ｚｋの１ステップ予測値Ｈｋｘ＾ｋ−１｜ｋ−１；フィルタ方程式によって与えられる。
【請求項１２】
請求項１乃至１１のいずれかに記載のシステム同定装置において、
前記フィルタを適用し、状態推定値ｘ＾ｋ｜ｋを求め、
擬似エコーｄ＾ｋを次式（６６）のように推定し、これをエコーから差し引くことにより、エコーをキャンセルすることを特徴とするシステム同定装置。
【数９】
2019205049.tif 000064
ここで、
［ｈ＾０［ｋ］、・・・、ｈ＾Ｎ−１［ｋ］］Ｔはｘ＾ｋ−１｜ｋ−１乃至ｘ＾ｋ｜ｋとし、
ｕｋは時刻ｔｋ（＝ｋＴ；Ｔは標本化周期）における受信信号、
Ｎはタップ数（既知）とする。
【請求項１３】
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２)によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とするシステム同定方法。
【数１０】
2019205049.tif 000065
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【請求項１４】
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定方法。
【数１１】
2019205049.tif 000066
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項１５】
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定方法。
【数１２】
2019205049.tif 000067
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項１６】
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムであって、
処理部が、前記フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
前記処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム。
【数１３】
2019205049.tif 000068
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【請求項１７】
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムであって、
前記処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム。
【数１４】
2019205049.tif 000069
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項１８】
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムであって、
処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム。
【数１５】
2019205049.tif 000070
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項１９】
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、前記フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
前記処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【数１６】
2019205049.tif 000071
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【請求項２０】
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【数１７】
2019205049.tif 000072
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【請求項２１】
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【数１８】
2019205049.tif 000073
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、システム同定装置及び方法及びプログラム及び記憶媒体に係り、特に、大規模な音響系または通信系を高速かつ数値的に安定に同定するシステム同定装置及び方法、及び、システム同定を実行させるためのプログラム及び記憶媒体に関する。
【背景技術】
【０００２】
システム同定（ｓｙｓｔｅｍ  ｉｄｅｎｔｉｆｉｃａｔｉｏｎ）とは、入出力データに基づいてシステムの入出力関係の数理モデル（伝達関数やインパルス応答等）を推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生およびアクティブ騒音制御、補聴器などのフィードバックキャンセラがある。従来、システム同定における適応アルゴリズムとして、ＬＭＳ（Ｌｅａｓｔ  Ｍｅａｎ  Ｓｑｕａｒｅ）アルゴリズム、ＲＬＳ（Ｒｅｃｕｒｓｉｖｅ  Ｌｅａｓｔ  Ｓｑｕａｒｅ）アルゴリズムやカルマンフィルタ（Ｋａｌｍａｎ  ｆｉｌｔｅｒ）が広く用いられて来た。一般に、線形システムの入出力関係は次のように表される。
【数１】
2019205049.tif 000003
ここで、ｕｋは入力、ｈｉはシステムのインパルス応答であり、ｙｋは観測時の雑音ｖｋを含んだ出力である。
詳細は、非特許文献１などに示されている。
【０００３】
１．ＬＭＳアルゴリズム
ＬＭＳは、入力ｕｋと出力ｙｋからシステムのインパルス応答ｘｋ＝［ｈ０、・・・、ｈＮ−１］Ｔを次のように推定する。
【数２】
2019205049.tif 000004
ただし、Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔであり、μ＞０はステップサイズと呼ばれる。また、ステップサイズμをμ／（||ＨＴｋ||２＋δ）で置き換えたものをＮＬＭＳアルゴリズムと呼ぶ。ここで、δは無入力でも分母が０とならないようにあらかじめ定める正の定数である。
【０００４】
２．ＲＬＳ
ＲＬＳは、入力ｕｋと出力ｙｋからシステムのインパルス応答ｘｋ＝［ｈ０、・・・、ｈＮ−１］Ｔを次のように推定する。
【数３】
2019205049.tif 000005
（なお、「＾」、「ｖ」は、推定値の意味であり、数式で示すように、本来、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する場合がある。以下同様。）
【０００５】
３．カルマンフィルタ
状態空間モデルで表される線形システム
【数４】
2019205049.tif 000006
の状態ｘｋの最小分散推定値ｘ＾ｋ｜ｋは、次のカルマンフィルタによって得られる。
【数５】
2019205049.tif 000007
ｘｋ：状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。
ｙｋ：観測信号；フィルタの入力となり、既知である。
Ｈｋ：観測行列；既知である。
ｖｋ：観測雑音；未知である。
ρ：忘却係数；一般に試行錯誤で決定される。
Ｋｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾０｜−１：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉが用いられる。
【０００６】
その他、従来、特許文献１〜３、非特許文献２〜６に示す技術がある。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特許第４０６７２６９号公報
【特許文献２】特許第４４４４９１９号公報
【特許文献３】特許第５１９６６５３号公報
【非特許文献】
【０００８】
【非特許文献１】S．Haykin、Adaptive Filter Theory、3rd Edition、Prentice-Hall、1996．
【非特許文献２】K．Nishiyama、Derivation of a Fast Algorithm of Modified H∞ Filters、Proceedings of IEEE International Conference on Industrial Electronics、Control and Instrumentation、RBC-II、pp．462-467、2000．
【非特許文献３】K．Nishiyama、An H∞ Optimization and Its Fast Algorithm for Time-Variant System Identification、IEEE Transactions on Signal Processing、52、5、pp．1335-1342、2004．
【非特許文献４】K．Nishiyama、A Unified View of Adaptive Algorithms for Finite Impulse Response Filters using the H Framework、Signal Processing (Elsevier)、97、pp．55-64、April 2014．
【非特許文献５】B．Hassibi、A．H．Sayed、and T．Kailath、Indefinite-Quadratic Estimation and Control、1st Edition、SIAM、1999．
【非特許文献６】G．Glentis、K．Berberidis、and S．Theodoridis、Efficient least squares adaptive algorithms for FIR transversal filtering、IEEE Signal Processing Magazine、16、4、pp．13-41、1999．
【発明の概要】
【発明が解決しようとする課題】
【０００９】
現在、システム同定において最も広く用いられている適応アルゴリズムはＬＭＳである。ＬＭＳは計算量は少ないが収束速度が非常に遅いと云った課題があった。一方、ＲＬＳやカルマンフィルタのΣ＾ｋ｜ｋ−１やＰｋは本来、正定行列（ｐｏｓｉｔｉｖｅ−ｄｅｆｉｎｉｔｅ  ｓｙｍｍｅｔｒｉｃ  ｍａｔｒｉｘ）であるが単精度（例：３２ｂｉｔ）で計算した場合には、しばしばこの条件を満たさず、これがカルマンフィルタやＲＬＳを数値的に不安定する要因の一つとなった。また、計算量がＯ（Ｎ２）であるため、状態ベクトルｘｋの次元（タップ数）Ｎが大きい場合、１ステップ当たりの演算回数が急激に増大し、実時間処理には必ずしも適するとはいえない場合があった。
本発明は、以上の点に鑑み、大規模な音響系または通信系を高速かつ高精度に同定することを、ひとつの目的とする。また、本発明は、このような同定する方法／装置等をフィードバックキャンセラに適用することを、他の目的とする。
【課題を解決するための手段】
【００１０】
本発明の第１の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜（６２）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置が提供される。
【数６】
2019205049.tif 000008
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【００１１】
本発明の第２の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置が提供される。
【数７】
2019205049.tif 000009
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【００１２】
本発明の第３の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置が提供される。
【数８】
2019205049.tif 000010
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【００１３】
本発明の第４の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２)によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とするシステム同定方法が提供される。
【数９】
2019205049.tif 000011
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【００１４】
本発明の第５の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定方法が提供される。
【数１０】
2019205049.tif 000012
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【００１５】
本発明の第６の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定方法が提供される。
【数１１】
2019205049.tif 000013
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【００１６】
本発明の第７の解決手段によると、
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、前記フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
前記処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
【数１２】
2019205049.tif 000014
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
【００１７】
本発明の第８の解決手段によると、
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
【数１３】
2019205049.tif 000015
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【００１８】
本発明の第９の解決手段によると、
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
【数１４】
2019205049.tif 000016
Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
【発明の効果】
【００１９】
本発明によると、以上のように、大規模な音響系または通信系を高速かつ高精度に同定することができる。また、本発明によると、このような同定する方法／装置等をフィードバックキャンセラに適用することができる。
【図面の簡単な説明】
【００２０】
【図１】γｆ−σｗ２空間の解の分布と性質についての説明図。
【図２】高速Ｈ∞フィルタ族によるｘ＾ｋ｜ｋ、ｚ＾ｋ｜ｋ−１＝Ｈｋｘ＾ｋ−１｜ｋ−１を求めるフローチャート。
【図３】スカラー存在条件のチェックについてのフローチャート。
【図４】通信系とエコーについての説明図。
【図５】エコーキャンセラの原理図。
【図６】本実施の形態（Ｃ−ＦＨＦ）と他の方法のタップ２乗誤差の推移についての説明図。
【図７】音声信号の説明図。
【図８】図７の音声信号を入力としたときの結果についての説明図。
【図９】高速Ｈ∞フィルタの処理と時間計算量の説明図（表１）。
【図１０】Ｃ−高速Ｈ∞フィルタの処理と時間計算量の説明図（表２）。
【図１１】Ｃ−白色化高速Ｈ∞フィルタの処理と時間計算量の説明図（表３）。
【図１２】Ｃ−ＮＬＭＳアルゴリズムのフィルタの処理と時間計算量の説明図（表４）。
【図１３】エコーパスのインパルス応答を示す図（表５）。
【図１４】本実施の形態に関するハードウェアの構成図。
【図１５】Ｃ−白色化高速Ｈ∞フィルタの最適版の一例の処理と時間計算量の説明図。
【図１６】Ｃ−ＮＬＭＳアルゴリズムのフィルタの最適版の一例の処理と時間計算量の説明図。
【発明を実施するための形態】
【００２１】
１．概要
本発明及び／又は本実施の形態の同定法により、入出力データから未知システムを実時間で高速かつ高精度に同定することが可能となる。また、本同定法は入力信号における白色性の仮定が適切に導入できるため、アルゴリズムを効果的に簡略化でき、大幅な計算量の削減が可能である。これらは、例えば補聴器などの閉ループ系のフィードバックキャンセラに有効である。
本発明及び／又は本実施の形態によると、例えば、以下を提供することができる。
１）入出力データから未知システムを実時間で高速かつ高精度に同定する方法・装置・プログラム・記録媒体。
２）入力信号に対して白色性の仮定を導入し、計算量を大幅に削減したシステム同定法。
３）補聴器などの閉ループ系のフィードバックキャンセラ。
また、本発明及び／又は本実施の形態においては、上記課題を解決するために、入力信号の事前線形予測誤差と事後線形予測誤差の変換係数を利用して、先に提案した高速Ｈ∞フィルタ（ＦＨＦ）の計算量Ｏ（１１Ｎ）をＯ（７Ｎ）まで低減するＣ−高速Ｈ∞フィルタ（Ｃｏｎｖｅｒｓｉｏｎ  Ｆａｃｔｏｒ−ｂａｓｅｄ  Ｆａｓｔ  Ｈ∞  Ｆｉｌｔｅｒ：  Ｃ−ＦＨＦ）を導出する。さらに、本発明及び／又は本実施の形態によると入力信号に白色性の仮定を導入して、計算量をＯ（２Ｎ）まで削減するＣ−白色化高速Ｈ∞フィルタ（Ｃｏｎｖｅｒｓｉｏｎ  Ｆａｃｔｏｒ−ｂａｓｅｄ  Ｗｈｉｔｅｎ  Ｆａｓｔ  Ｈ∞  Ｆｉｌｔｅｒ：  Ｃ−ＷＦＨＦ）を提供する。また、その特殊な場合としてＣ−ＮＬＭＳアルゴリズム（Ｃｏｎｖｅｒｓｉｏｎ  Ｆａｃｔｏｒ−ｂａｓｅｄ  Ｎｏｒｍａｌｉｚｅｄ  Ｌｅａｓｔ  Ｍｅａｎ  Ｓｑｕａｒｅ  Ａｌｇｏｒｉｔｈｍ）を与える。
本発明及び／又は本実施の形態による高速同定法は、大規模な音響系や通信系の同定に適用可能であり、スマートスピーカーやテレビ会議システムなどにおける音響エコーキャンセラやアクティブ騒音制御、音場再生などに有効である。また、本同定法は入力信号に対して白色性の仮定が適切に導入できるため、アルゴリズムを効果的に簡略化でき、大幅な計算量の削減が可能である。これらの簡易化されたアルゴリズムは補聴器や高速データ伝送の全二重中継局などの閉ループ系のフィードバックキャンセラに有効である。
【００２２】
２．記号の説明
まず、本実施の形態で用いる主な記号を説明する。
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。（Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ）
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｚ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの出力ｚｋの推定値Ｈｋｘ＾ｋ｜ｋ；フィルタ方程式によって与えられる。
ｚ＾ｋ｜ｋ−１：観測信号ｙ０〜ｙｋ−１までを用いた時刻ｋの出力ｚｋの１ステップ予測値Ｈｋｘ＾ｋ−１｜ｋ−１；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾０｜−１：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉ、ε０＞０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１やゲイン行列Ｋｋから得られる。
ρ：忘却係数；γｆが決まればρ＝１−χ（γｆ）より決定される。γｆと切り離して独立に決定することも可能である。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
Σ０−１：状態の不確定性を表す重み行列の逆行列；便宜上Σ０＝ε０Ｉ，ε０＞０が用いられる。
ｅｆ，ｉ：フィルタ誤差；  ｅｆ，ｉ＝ｚｖｉ｜ｉ−Ｈｉｘｉ，ｚｖｉ｜ｉは出力推定値Ｈｉｘ＾i｜ｉである。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜ｋ：ＫＵｋのＮ＋１番目の１×２行ベクトル；ＫＵｋより得られる。
ｍｋ：ＫＵｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵｋより得られる。
γｆ：アッテネーションレベル；設計の際に与えられる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
Ｏ：計算量（オーダー）；例えば、乗算回数がＮ２＋２Ｎ＋３のときＯ（Ｎ２）と表す。
なお、記号の上に付される”＾”、”ｖ”は、推定値の意味である。また、”〜”、” ”、”Ｕ”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、ｘ、ｗ、Ｈ、Ｇ，Ｋ、Ｒ、Σ、等はベクトルや行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する場合がある。
【００２３】
２．システム同定のハードウェア及びプログラム
本システム同定方法又はシステム同定装置・システムは、その各手順をコンピュータに実行させるためのシステム同定プログラム、システム同定プログラムを記録したコンピュータ読み取り可能な記録媒体、システム同定プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
図１４は、本実施の形態に関するハードウェアの構成図である。
このハードウェアは、中央処理装置（ＣＰＵ）である処理部１０１、入力部１０２、出力部１０３、表示部１０４、記憶部１０５、インタフェース部（Ｉ／Ｆ）１０６を有する。また、処理部１０１、入力部１０２、出力部１０３、表示部１０４、記憶部１０５、インタフェース部（Ｉ／Ｆ）１０６は、スター又はバス等の適宜の接続手段で接続されている。記憶部１０５は、システム同定される「１．記号の説明」で示した既知のデータが必要に応じて記憶され、未知・既知のデータや計算されたＣ−高速Ｈ∞フィルタ等の本実施の形態の各フィルタに関するデータ・その他のデータが処理部１０１により、必要に応じて書き込み及び／又は読み出しされる。処理部１０１は、入出力データに基づいて、本実施の形態の各フィルタを表す式に従い計算処理を実行し、システムの入出力関係の数理モデルを推定（同定）する。処理部１０１は、例えば、状態推定値、出力予測値等を求めることができる。
【００２４】
３．本実施の形態に関連する同定法の手段
通常のＨ∞フィルタと区別するため、モデル集合と重み集合に関してさらに最適化したＨ∞フィルタを、特にハイパーＨ∞フィルタと呼ぶことにする（非特許文献２−４）。このＨ∞フィルタは状態のダイナミックスが未知の場合でも適用できる特長を備えている。
【００２５】
［ハイパーＨ∞フィルタ］
状態空間モデル
【数１５】
2019205049.tif 000017
を満たす状態推定値ｘ＾ｋ｜ｋ（あるいは出力推定値ｚｖｋ｜ｋ）は、次のレベルγｆのハイパーＨ∞フィルタによって与えられる。
【数１６】
2019205049.tif 000018
であり、χ（γｆ）はχ（１）＝１、χ（∞）＝０を満たすγｆの単調減衰関数である。また、駆動行列Ｇｋは次のように生成される。
【数１７】
2019205049.tif 000019
ハイパーＨ∞フィルタの特徴は、エコーキャンセラ（状態推定）に重要な次の五つの性質、
１）  収束性、
２）  エコー経路変動に対する追従能、
３）  近端話者信号に対するロバスト性、
４）  入力信号特性に関する不感度、
５）  数値的安定性
を調整できる点にある。
【００２６】
図１は、γｆ−σｗ２空間の解の分布と性質についての説明図である。
図１の丸印は、ハイパーＨ∞フィルタのパラメータ空間における解を表している。それぞれの解は異なるγｆとσｗ２の値をもつハイパーＨ∞フィルタに対応する。特に、直線γｆ＝∞上の解は忘却係数ＲＬＳアルゴリズムと一致し、点（γｆ，σｗ２）＝（１，０）上ではステップサイズ１のＮＬＭＳアルゴリズムと一致する。よって、ハイパーＨ∞フィルタは従来の適応アルゴリズムであるＮＬＭＳやＲＬＳを含んでいることがわかる。
ハイパーＨ∞フィルタの計算量はＯ（Ｎ２）であり、このままでは実時間処理には適さない。しかし、観測行列Ｈｋがシフト特性Ｈｋ＝［ｕｋ，Ｈｋ−１（１），Ｈｋ−１（２），・・・，Ｈｋ−１（Ｎ−１）］をもつとき、計算量Ｏ（Ｎ）の高速アルゴリズムが開発されている（非特許文献２、非特許文献３）。
【００２７】
［高速Ｈ∞フィルタ］
観測行列Ｈｋがシフト特性を満たすとき、Σ０＝ε０Ｉ＞０をもつハイパーＨ∞フィルタは次式によって実行できる。
【数１８】
2019205049.tif 000020
であり、ｃｋ∈Ｒ２×１はＣｋ＝［ｃｋ，・・・，ｃｋ−Ｎ＋１］の第一列ベクトルであり、ｃｋ−ｉ＝０２×１、ｋ−ｉ＜０と仮定し、初期値はＫ−１＝０Ｎ×２、ａ−１＝０Ｎ×１、Ｓ−１＝１／ε０、ε０＞０、ｂ−１＝０Ｎ×１、ｘ＾−１｜−１＝０Ｎ×１と設定した。ここで、０ｍ×ｎはｍ×ｎ零行列を表す。
このとき、次のスカラー存在条件を満たさなければならない。
【数１９】
2019205049.tif 000021
【００２８】
４．本実施の形態の同定法の手段
４．１  準備
高速Ｈ∞フィルタの新しいアルゴリズムを導出するためにいくつかの命題を示す。その命題の証明には、次のように定義された（Ｎ＋１）×（Ｎ＋１）行列ＱＵｋ  が利用される。
【数２０】
2019205049.tif 000022
この行列ＱＵｋは、観測行列ＣＵｉのシフト特性から次のように二通りに分解できる。
【数２１】
2019205049.tif 000023
その分解されたＱＵｋは次式によって再帰的に求められる。
【数２２】
2019205049.tif 000024
【００２９】
命題１
ゲイン行列Ｋｋ：＝Ｑｋ−１ＣｋＴは疑似ゲイン行列Ｋ〜ｋ：＝Ｑｋ−１−１ＣｋＴを用いて次のように表すことができる。
【数２３】
2019205049.tif 000025
命題２
後方線形予測係数ｂｋ：＝−Ｑｋ−１ｔｋの新たな更新式は次式で与えられる。
【数２４】
2019205049.tif 000026
命題３
事後後方予測誤差ｅｋは事前後方予測誤差ｅ〜ｋを用いて次のように表すことができる。
【数２５】
2019205049.tif 000027
命題４
事後前方予測誤差ｅｋは事前前方予測誤差ｅ〜ｋを用いて次のように表すことができる。
【数２６】
2019205049.tif 000028
【００３０】
命題５
後方線形予測係数ｂｋは事後後方予測誤差ｅｋを用いて次のように再帰的に更新できる。
【数２７】
2019205049.tif 000029
命題６
前方線形予測係数ａ：＝−Ｑｋ−１−１ｔｋは事後前方予測誤差ｅｋを用いて次式によって再帰的に更新できる。
【数２８】
2019205049.tif 000030
命題７
疑似ゲイン行列Ｋ〜ｋを計算するもう一つの形式は次式で与えられる。
【数２９】
2019205049.tif 000031
命題８
事前後方予測誤差ｅ〜ｋはμ〜ｋを用いて次のように表すことができる。
【数３０】
2019205049.tif 000032
命題９
式（３６）の２×２行列Ｒｅ，ｋは次式で再帰的に更新できる。
【数３１】
2019205049.tif 000033
【００３１】
４．２  アルゴリズムの導出
命題１から命題９を用いれば、前方線形予測係数ａｋ  と後方線形予測係数ｂｋ  が事後前方予測誤差ｅｋと事後後方予測誤差ｅｋを用いて更新される新たな高速Ｈ∞フィルタ（Ｃ−高速Ｈ∞フィルタ）が次のように得られる。
【数３２】
2019205049.tif 000034
ただし、各パラメータの設定や初期化は高速Ｈ∞フィルタの場合と同じである。このＣ−高速Ｈ∞フィルタは、変換係数（ｃｏｎｖｅｒｓｉｏｎ  ｆａｃｔｏｒ）χｋを用いた高速Ｈ∞フィルタの事後予測誤差形式に当たる。
【００３２】
４．３  計算量の評価
Ｃ−高速Ｈ∞フィルタは高速Ｈ∞フィルタと理論的に等価であるが、計算量は少なくて済む。Ｃ−白色化高速Ｈ∞フィルタは、入力信号が白色であると仮定して、アルゴリズムを大幅に縮退させたＣ−高速Ｈ∞フィルタである。Ｃ−ＮＬＭＳアルゴリズムはパラメータが特定の値であるときのＣ−白色化高速Ｈ∞フィルタに対応し、理論的にはＮＬＭＳアルゴリズムと等価であるが、計算量はＯ（２Ｎ）で済む。
図９に、高速Ｈ∞フィルタの処理と時間計算量の説明図（表１）を示す。
高速Ｈ∞フィルタの各式の計算に必要な乗算と加算の回数を図９（表１）にまとめた（特許文献１参照、比較のために記載）。高速Ｈ∞フィルタを実行するためには、図２のＣｏｒｅ  Ａｌｇｏｒｉｔｈｍで表の数式を上から順に（式番号が大きい番号から小さい番号へ降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
【００３３】
図１０に、Ｃ−高速Ｈ∞フィルタの処理と時間計算量の説明図（表２）を示す。
図１０（表２）は、Ｃ−高速Ｈ∞フィルタに対する乗算と加算の回数を示している。ここで、演算回数の算出に当たり、Ｋｋ，Ｃｋ，ｅ〜ｋ，ｅｋ，およびＷの構造が利用された。例えば、Ｋｋ−１Ｗｅ〜ｋは、Ｋｋ−１（：，１）＝Ｋｋ−１（：，２）とｅ〜ｋ（１，１）＝ｅ〜ｋ（２，１）を用いて、（１−γｆ−２）Ｋｋ−１（：，１）ｅ〜ｋ（１，１）に置き換えられた。同様に、Ｃｋ−１とａｋ−１の積を求めるに当たり、Ｃｋ−１（２，：）ａｋ−１はＣｋ−１（１，：）ａｋ−１で置換された。同じ組合せの２項の乗算が複数ある場合は最初の計算のみをカウントした。χｋ−１−１は前のステップで計算したχｋ−１を用いた。スカラーによるベクトルの除算は、スカラーの逆数によるベクトルの乗算として処理した。加算回数の場合も同様な方法で算出した。
Ｃ−高速Ｈ∞フィルタは、Ｎが十分に大きい場合、単位サンプル当たり７Ｎ回の乗算と７Ｎ回の加算のみで実行できる。一方、高速Ｈ∞フィルタは、単位サンプル当たり１１Ｎ回の乗算と１０Ｎ回の加算が必要である。両フィルタとも単位サンプル当たり３回の乗算が必要である。また、いずれの高速Ｈ∞フィルタも空間計算量(メモリ)はＯ（６Ｎ）である。
【００３４】
４．４  高速Ｈ∞フィルタ族の処理
図２に、高速Ｈ∞フィルタ族によるｘ＾ｋ｜ｋ、ｚ＾ｋ｜ｋ−１＝Ｈｋｘ＾ｋ−１｜ｋ−１を求めるフローチャートを示す。ここで、高速Ｈ∞フィルタ族とは、例えば、Ｃ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）、Ｃ−白色化高速Ｈ∞フィルタ（Ｃ−ＷＦＨＦ）、Ｃ−ＮＬＭＳアルゴリズム・フィルタ等を示す。高速Ｈ∞フィルタとＣ−高速Ｈ∞フィルタでは、このフローチャートに、後述のように、図３の存在条件のチェック機能が追加される場合がある。また、必要に応じてｚ＾ｋ｜ｋを求める。
【００３５】
以下に、高速Ｈ∞フィルタ族の処理について、フローチャートを参照して説明する。
［ステップＳ１０１、初期化］処理部１０１は、高速Ｈ∞フィルタの再帰変数を初期化する。たとえば、処理部１０１は、記憶部１０５から再帰式の初期状態を読み出し、又は初期状態を入力部１０２又はＩ／Ｆ１０６から入力し、図示のように定める。
［ステップＳ１０３］処理部１０１は、事前に決定したデータ長Ｌを超えているか判断する。たとえば処理部１０１は、時刻ｋとデータ数Ｌと比較する。なお、Ｌはあらかじめ定められた最大データ数を表すことができる。処理部１０１は、時刻ｋがデータ数Ｌ以下であれば次のステップＳ１０５に進む。処理部１０１は、時刻ｋが超えていればフィルタゲインを初期化し、フィルタをリスタートするか、あるいは終了する（Ｓ１１３）。（なお、不要であれば条件文を取り除くことができる。または、必要に応じて再スタートを行ってもよい。）
［ステップＳ１０５、入力］処理部１０１は、時刻ｋがＬ以下であれば、拡大観測行列ＣＵｋを入力信号とシフト特性を利用して更新する。たとえば、処理部１０１は、入力信号ｕｋを入力部１０２又はＩ／Ｆ１０６から入力し、ＣＵｋを図示のように設定する。なお、入力信号ｕｋを記憶部１０５から読み出すようにしてもよい。
ここで、入力信号ｕｋからなる２×（Ｎ＋１）拡大観測行列
【数３３】
2019205049.tif 000035
［ステップＳ１０７、Ｃｏｒｅ  Ａｌｇｏｒｉｔｈｍ］処理部１０１は、Ｃｏｒｅ  Ａｌｇｏｒｉｔｈｍでは図９（表１）から図１２（表４）及び図１５、図１６のいずれかのアルゴリズムを実行する（詳細は後述）。なお、これらのいずれかのアルゴリズムにおいて、状態推定値ｘ＾ｋ｜ｋを緩和係数（ステップサイズ）μ（ここで、μは、μ＞０の実数）で更新するようにしてもよい。この場合、例えば、式（４９）の右辺第２項にμを乗算することで、緩和係数（ステップサイズ）μで更新することができる。また、必要に応じてアルゴリズムの実行後又は実行前等の適宜のタイミングで、図３の存在条件の判定のための処理を挿入するようにしてもよい（詳細は後述）。
［ステップＳ１０９、推定値］処理部１０１は、ステップＳ１０７により計算された状態推定値ｘ＾ｋ｜ｋ（あるいは出力推定値ｚ＾ｋ｜ｋ）を記憶部１０５に記憶、及び／又は出力部１０３若しくはＩ／Ｆ１０６から出力する。
［ステップＳ１１１］処理部１０１は、時刻ｋを進ませて（ｋ＝ｋ＋１）、ステップＳ１０３に戻り、データがある限り又は所定回数繰り返し続ける。
なお、処理部１０１は、ステップＳ１０５〜Ｓ１０９等の各ステップで求めた適宜の中間値及び最終値、存在条件の値等を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。さらに、処理部１０１は、入力部１０２若しくはＩ／Ｆ１０６からの指示により又は予め定められた処理により、各値によるデータ又はグラフ・表・図等を作成して、記憶部１０５に記憶、及び／又は出力部１０３若しくはＩ／Ｆ１０６から出力するようにしてもよい。
【００３６】
図３に、スカラー存在条件のチェックについてのフローチャートを示す。
図２のフローチャートにさらに図３の存在条件のチェック機能が追加してもよい。たとえば、処理部１０１は、ステップＳ１０７の処理の実行後又は実行前等の適宜のタイミングで、図示のスカラー存在条件の判定ステップを実行し、式（３０）を満たさなければ処理を終了する又はリスタートするようにしてもよい。なお、このステップは、他にも、図２のフローチャートの適宜のひとつ又は複数のステップや、時刻ｋを進ませるステップ（Ｓ１１１）の前後等の適宜のタイミングで、処理部１０１が実行することができる。なお、図中、ＥＸＣ（ｋ）は、前述の式（３０）の左辺に相当する。
【００３７】
４．５  白色性仮定の導入
図１１に、緩和係数μを導入したＣ−白色化高速Ｈ∞フィルタの処理と時間計算量の説明図（表３）を示す。
入力信号が白色雑音と仮定すれば、線形予測係数ａｋとｂｋは０となる。このとき、ｅ〜ｋ＝ｃｋ＋Ｃｋ−１ａｋ−１＝ｃｋ、ｅｋ＝ｃｋ＋Ｃｋ−１ａｋ＝ｃｋ、ｅ〜ｋ＝ｃｋ−Ｎ＋Ｃｋｂｋ−１＝ｃｋ−Ｎ、ｅｋ＝ｃｋ−Ｎ＋Ｃｋｂｋ＝ｃｋ−Ｎ、Ｋ〜ｋ＝ｍ〜ｋ、Ｓｋは不要となり、式（４９）、式（５０）、式（５５）、式（５８）、式（６０）、式（６１）、式（６２）は次のように整理できる。ただし、μは緩和係数(ステップサイズ）で、μ＞０の実数であり、事前に与えられる。
なお、ｅ〜ｋの式（５６）、ｅｋの式（５３）は、本来の定義式とは別に、理論的に等価な表現が複数あり、上述のように表すことができる。また、Ｃ−ＷＦＨＦの計算で、後方予測関係の変数（アンダーバーが付いたもの）は、利用しないことができるので、それらの計算は不要となる。
【数３４】
2019205049.tif 000036
これより、Ｃ−高速Ｈ∞フィルタは表３のＣ−白色化高速Ｈ∞フィルタに縮退する。ただし、ｚ＾ｋ｜ｋ＝Ｈｋｘ＾ｋ｜ｋは省略した。
Ｃ−白色化高速Ｈ∞フィルタを実行するためには、図２のＣｏｒｅ  Ａｌｇｏｒｉｔｈｍで図１１（表３）の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
【００３８】
４．６  Ｃ−ＮＬＭＳアルゴリズム
図１２に、Ｃ−ＮＬＭＳアルゴリズムのフィルタの処理と時間計算量の説明図（表４）を示す。
γｆ＝１．０、ρ＝１．０のとき、Ｓｋ＝Ｓｋ−１＝Ｓ−１＝１／ε０となる。これより、Ｃ−白色化高速Ｈ∞フィルタは次のように縮退する。
【数３５】
2019205049.tif 000037
これはステップサイズμ＝１のＮＬＭＳアルゴリズムと理論的に等価である。よって、緩和係数μ（μ＞０の実数）を式（４９）に導入すれば、表４のステップサイズμのＣ−ＮＬＭＳアルゴリズムが得られる。このアルゴリズムの計算量はＯ（２Ｎ）であり、Ｓ−１＝１／ε０はＮＬＭＳアルゴリズムの正則化項δに一致する。
Ｃ−ＮＬＭＳアルゴリズムのフィルタを実行するためには、図２のＣｏｒｅ  Ａｌｇｏｒｉｔｈｍで図１２（表４）の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
【００３９】
４．７  Ｃ−高速Ｈ∞フィルタの白色化の妥当性
高速Ｈ∞フィルタは入力信号の線形予測係数を０に置くことによって、次のように白色化することができる。
【数３６】
2019205049.tif 000038
ここで、Ｋｋ−１（１：Ｎ−１，：）は、Ｋｋ−１の第１行〜第（Ｎ−１）行までからなる行列を表す。
この白色化高速Ｈ∞フィルタ（ＷＦＨＦ）は、γｆ＝∞のとき、次のように縮退する。
【数３７】
2019205049.tif 000039
このとき、フィルタゲインの分母の正規化項  １＋γｆ−２ＨｋＫｋ（：，１）  は、１となり、入力信号ｕｋの振幅が急激に増大するとフィルタゲインが大きくなり過ぎて不安定となる場合が想定される（γｆ＜∞の場合は、分母にＨｋＫｋ（：，１）があるのでゲインレベルが制限される）。これより、従来の高速Ｈ∞フィルタの白色化に課題があることがわかる。また、ＷＦＨＦの計算量はＯ（３Ｎ）であった。
この課題を解決して、白色化を適切に実行することができるようにしたのが、本発明及び／又は本実施の形態のＣ−ＷＦＨＦである。白色化のために、例えば、数多くあるパラメータの中から、線形予測係数ａｋとｂｋを選択して０とし、また、ｅ〜ｋ、ｅｋ、ｅ〜ｋＮ、ｅｋ、の多数の数式から適切なものをあてはめ、さらに、Ｋ〜ｋ＝ｍ〜ｋ、という式を適用し、また、Ｓｋは不要となるようにした。このように、数多く存在するパラメータの中から特定のものを選択して特定の値とし、数多く存在する式の中から特定のパラメータについて特定の式を適用することは、極めて高度かつ困難なことであって、これにより本発明および／本実施の形態は、白色化を適切に行うことを達成したものである。
【００４０】
４．８  Ｃ−白色化高速Ｈ∞フィルタ（Ｃ−ＷＦＨＦ）、Ｃ−ＮＬＭＳアルゴリズム（Ｃ−ＮＬＭＳ）の最適版の検討
Ｃ−ＷＦＨＦ、Ｃ−ＮＬＭＳの実装を考慮した場合における計算式の最適化の一例及び演算量の見積もりを行なった。ここで、最適化とは、特に、「組み込みのための最適化」をいう。ここでは、例えば、組み込まれたプログラムの演算回数が最小になる式（形式）を求めることを、最適化と呼んでいる。なお、「Ｃ−ＦＨＦとの関係性を維持した形式」と「組み込みに適した形式（最適版）」は、理論的に等価である。例えば、「Ｃ−ＦＨＦとの関係性を維持した形式」はそのアルゴリズムの意味がわかり、一方、「組み込みに適した形式」は実際にＣＰＵにプログラムするのに適した形式である。
ＤＳＰ（digital signal processor）等のプロセッサやハードウェアに組み込む際には掛け算（Ｍｕｌ．）、割り算（Ｄｉｖ．）、足し算（Ａｄｄ．）によって必要となるサイクル数が異なるため、それぞれの演算回数からサイクル数の見積もりを行なった。実際にはメモリやレジストリからの読み込みや並列化演算などによって必要なサイクル数は大きく変化するが組み込み方法によるサイクル数の変化は考慮せず算出した。
図１５に、Ｃ−白色化高速Ｈ∞フィルタの最適版の一例の処理と時間計算量の説明図を示す。
このとき、Ｃ−白色化高速Ｈ∞フィルタの式（６０）、式（５８）、式（５５）、式（４９）、式（５０）は、次のように整理できる。ここで、Ｋ~ｋ−１（１：Ｎ−１，１）はＫ~ｋ−１の第１行〜第Ｎ−１行までの第１列の行列を表す。
なお、表３のＣ−白色化高速Ｈ∞フィルタ（Ｃ−ＷＦＨＦ）の式にμ（ステップサイズ、μ＞０の実数）を付けて、Ｃ−ＷＦＨＦがＣ−ＮＬＭＳを完全に含むように表現することができる。よって、Ｃ−ＮＬＭＳは、γｆ＝１，ρ＝１のときのＣ−ＷＦＨＦと完全に一致する。
ここで、最適化の式とは演算回数が少なくなるように式を組み換えものであり、理論的に等価になる。例えば、Ｃ−ＷＦＨＦの式（５８）の第１列の第１成分から第Ｎ成分を抽出して利用したのが、以下のＫ〜ｋ（：，１）の式である。
【数３８】
2019205049.tif 000040
なお、Ｃ−ＷＦＨＦの式（５８）及び（６１）（６２））から、
Ｋ〜ｋ−１（Ｎ，１）＝ｕｋ−Ｎ／Ｓｋ−Ｎ−１
という関係式が導出され、これによりＣ−ＷＦＨＦの式（５５）から上述のＲｅ，ｋ（１，１）の式が導かれる。
よって、実際の計算に必要な演算はＲｅ，ｋの１行１列目のスカラー値であることから、Ｃ−白色化高速Ｈ∞フィルタの式（５５）は上述した式のように省略することができる。
Ｃ−白色化高速Ｈ∞フィルタを実行するためには、図２のＣｏｒｅ  Ａｌｇｏｒｉｔｈｍで図表の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ｄｉｖ．」は割算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
このように、この最適化によると、先に示したＣ−白色化高速Ｈ∞フィルタより計算量を一層削減することができる。
【００４１】
図１６に、Ｃ−ＮＬＭＳアルゴリズムのフィルタの最適版の一例の処理と時間計算量の説明図を示す。
このとき、Ｃ−ＮＬＭＳアルゴリズムのフィルタの式（５８）、式（５５）、式（４９）、式（５０）は、次のように整理できる。なお、例えば、Ｃ−ＮＬＭＳの式（５８）の第１列の第１成分から第Ｎ成分を抽出して利用したのが、以下のＫ〜ｋ（：，１）の式である。
【数３９】
2019205049.tif 000041
（μはステップサイズ、μ＞０の実数）
なお、Ｃ−ＮＬＭＳの式（５８）から、
Ｋ〜ｋ−１（Ｎ，１）＝ε０ｕｋ−Ｎ
という関係式が導出され、これによりＣ−ＮＬＭＳの式（５５）から上述のＲｅ，ｋ（１，１）の式が導かれる。
よって、実際の計算に必要な演算はＲｅ，ｋの１行１列目のスカラー値であることから、Ｃ−ＮＬＭＳアルゴリズムのフィルタの式（５５）は上述した式のように省略することができる。
Ｃ−ＮＬＭＳアルゴリズムのフィルタを実行するためには、図２のＣｏｒｅ  Ａｌｇｏｒｉｔｈｍで図表の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ｄｉｖ．」は割算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
このように、この最適化の一例によると、先に示したＣ−ＮＬＭＳアルゴリズムのフィルタより計算量を一層削減することができる。これはδ＝１／ε０のＮＬＭＳアルゴリズムと一致する。
【００４２】
５．エコーキャンセラ
５．１  準備
以下に、実施の形態として、エコーキャンセラの例を取り上げ、本同定アルゴリズムの動作を確認する。その前に、通信エコーキャンセラについて簡単に説明する。最近注目されているスマートスピーカーなどの音声インターフェイスに不可欠な音響エコーキャンセラも原理は同じである。
図４に、通信系とエコーについての説明図を示す。
国際電話など長距離電話回線では、信号増幅などの理由から４線式回線が用いられている。一方、加入者回線は比較的短距離なので２線式回線が使用されている。２線式回線と４線式回線の接続部には図４のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Ｂからの信号（音声）は話者Ａのみに到達する。しかし、一般には完全に整合をとることはできず、受信信号の一部は４線式回線に漏れ、増幅された後、再び送信者（話者Ｂ）に戻る現象が起こる。これがエコー（ｅｃｈｏ）である。エコーは伝送距離が長くなるにつれて影響が大きくなり、著しく通話の品質を劣化させることが想定される（パルス伝送においては近距離であってもエコーによる通話品質の劣化は問題となることが想定される）。
【００４３】
図５に、エコーキャンセラの原理図を示す。
そこで、図５のようにエコーキャンセラ（ｅｃｈｏ  ｃａｎｃｅｌｌｅｒ）を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、その結果を利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
エコーパスのインパルス応答の推定は、残留エコーｅｋの平均２乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Ａからの信号（音声）である。一般に、話者２人が同時に話し始めた（ダブルトーク）ときはインパルス応答の推定を中断する。
次に、このエコーキャンセリング問題の数理モデルを作成する。まず、受信信号ｕｋがエコーパスへの入力信号となることを考慮すれば、エコーパスのインパルス応答ｈｉ［ｋ］により、エコーｄｋの観測値ｙｋは次式で表される。
【数４０】
2019205049.tif 000042
ここで、ｕｋ、ｙｋはそれぞれ時刻ｋにおける受信信号とエコーの観測値を表し、ｖｋは時刻ｋにおける平均値０の回線雑音とし、タップ数Ｎは既知とする。このとき、インパルス応答の推定値ｈ＾ｉ［ｋ］が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
【数４１】
2019205049.tif 000043
これをエコーから差し引けば（ｙｋ−ｄ＾ｋ≒０）、エコーをキャンセルすることができる。ただし、ｋ−ｉ＜０のときｕｋ−ｉ＝０とする。
以上より、キャンセリング問題は直接観測可能な受信信号ｕｋとエコーの観測値ｙｋからエコーパスのインパルス応答ｈｉ［ｋ］を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにＨ∞フィルタを適用するには、まず式（６５）を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、｛ｈｉ［ｋ］｝Ｎ−１ｉ＝ｏを状態変数ｘｋとし、ｗｋ程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
【数４２】
2019205049.tif 000044
【００４４】
５．２  動作の確認
図１３は、エコーパスのインパルス応答を示す図（表５）である。
次に、エコーパスの同定問題を用いて、Ｃ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）の妥当性を確認すると共に、高速カルマンフィルタ（ＦＫＦ）と高速トランスバーサルフィルタ（ＦＴＦ）に対する優位性を示す。未知システムは長さＮ＝１２８をもつインパルス応答｛ｈｉ（ｋ）｝Ｎ−１ｉ＝０で表され、時刻ｋ＝４０００で一度シフトし、その後はその値を保持した。
このとき、エコーの観測値ｙｋは次式で表される。
【数４３】
2019205049.tif 000045
ただし、｛ｈｉ｝２３ｉ＝０は表５の値を採用し、その他｛ｈｉ｝１２７ｉ＝２４は０とした。また、ｖｋは平均値０、標準偏差σｖ＝１０−４０の定常なガウス白色雑音とした。
また、受信信号ｕｋは次のように２次のＡＲモデルで近似した。
ｕｋ＝α１ｕｋ−１＋α２ｕｋ−２＋ｗ’ｋ                          (71)
ただし、α１＝０．７、α２＝０．１とし、ｗ’ｋは平均値０、標準偏差σｗ’＝０．０２の定常なガウス白色雑音とした。
【００４５】
図６に、本実施の形態（Ｃ−ＦＨＦ）と他の方法のタップ２乗誤差の推移についての説明図を示す。
図６は、Ｃ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）のタップ２乗誤差||ｘ＾ｋ｜ｋ−ｘｋ||２の推移を、高速Ｈ∞フィルタ（ＦＨＦ）、高速カルマンフィルタ（ＦＫＦ）、高速トランスバーサルフィルタ（ＦＴＦ）と比較した結果を示す。１．３＜γｆ＜∞のとき、計算量の違いはあるが、高速Ｈ∞フィルタ（ＦＨＦ）とＣ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）の性能の差はほとんど見られない。興味深いことに、このように観測雑音が非常に小さい場合、高速Ｈ∞フィルタとＣ−高速Ｈ∞フィルタは、高速カルマンフィルタと高速トランスバーサルフィルタと比べて高い収束性能と追従性能を発揮した。
次に、図７に、音声信号の説明図を示す。
また、図８に、図７の音声信号を入力としたときの結果についての説明図を示す。
この場合も、Ｃ−高速Ｈ∞フィルタは高速カルマンフィルタや高速トランスバーサルフィルタより優れた収束性能と追従性能を示した。
【００４６】
６．  命題の証明
以下に、上述の各命題の証明について説明する。
命題１の証明
【数４４】
2019205049.tif 000046
命題２の証明
【数４５】
2019205049.tif 000047
【００４７】
命題３の証明
命題１と命題２を用いれば、事後後方予測誤差ｅｋは次のように事前後方予測誤差ｅ〜ｋを用いて表すことができる。
【数４６】
2019205049.tif 000048
命題４の証明
命題１を用いれば、事後前方予測誤差ｅｋは次のように事前前方予測誤差ｅ〜ｋを用いて表すことができる。
【数４７】
2019205049.tif 000049
命題５の証明
命題１、命題２および命題３を用いれば、後方線形予測係数ｂｋの更新式は次のように書き換えることができる。
【数４８】
2019205049.tif 000050
命題６の証明
命題１と命題４を式（２９）に適用すれば、前方線形予測係数ａｋの新たな更新式が得られる。
【００４８】
命題７の証明
（Ｎ＋１）×（Ｎ＋１）行列ＱＵｋは次のように二通りに分解できる。
【数４８】
2019205049.tif 000051
逆行列の補助定理を上記のように分解されたＱＵｋ−１に適用すれば次式が得られる。
【数５０】
2019205049.tif 000052
式（７９）の両側に右からＣＵｋＴを掛ければ、次のように整理できる。
【数５１】
2019205049.tif 000053
このとき、式（８３）の両辺の各項を比較すれば次式が得られる。
【数５２】
2019205049.tif 000054
【００４９】
命題８の証明
式（８３）の両辺の第２項からμ〜ｋ＝ｅ〜Ｔｋ／Ｓｋ−１が得られ、これより式（４６）が導かれる。
命題９の証明
式（７９）の両辺に両側からＣＵｋとＣＵｋＴを掛ければ次式が得られる。
【数５３】
2019205049.tif 000055
これに命題８を適用すれば証明は完結する。
【符号の説明】
【００５０】
１０１  処理部
１０２  入力部
１０３  出力部
１０４  表示部
１０５  記憶部
１０６  インタフェース部
【図１】
2019205049.tif 000074
【図２】
2019205049.tif 000075
【図３】
2019205049.tif 000076
【図４】
2019205049.tif 000077
【図５】
2019205049.tif 000078
【図６】
2019205049.tif 000079
【図７】
2019205049.tif 000080
【図８】
2019205049.tif 000081
【図９】
2019205049.tif 000082
【図１０】
2019205049.tif 000083
【図１１】
2019205049.tif 000084
【図１２】
2019205049.tif 000085
【図１３】
2019205049.tif 000086
【図１４】
2019205049.tif 000087
【図１５】
2019205049.tif 000088
【図１６】
2019205049.tif 000089