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