JP6894402B2 - システム同定装置及び方法及びプログラム及び記憶媒体 - Google Patents

システム同定装置及び方法及びプログラム及び記憶媒体 Download PDF

Info

Publication number
JP6894402B2
JP6894402B2 JP2018098503A JP2018098503A JP6894402B2 JP 6894402 B2 JP6894402 B2 JP 6894402B2 JP 2018098503 A JP2018098503 A JP 2018098503A JP 2018098503 A JP2018098503 A JP 2018098503A JP 6894402 B2 JP6894402 B2 JP 6894402B2
Authority
JP
Japan
Prior art keywords
state
filter
unknown
given
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2018098503A
Other languages
English (en)
Other versions
JP2019205049A (ja
Inventor
清 西山
清 西山
政浩 春原
政浩 春原
信彦 昼間
信彦 昼間
智穂 春田
智穂 春田
誠 舘野
誠 舘野
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Iwate University
Rion Co Ltd
Original Assignee
Iwate University
Rion Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Iwate University, Rion Co Ltd filed Critical Iwate University
Priority to JP2018098503A priority Critical patent/JP6894402B2/ja
Priority to US16/418,579 priority patent/US10863272B2/en
Priority to EP19176004.0A priority patent/EP3573233B1/en
Priority to CN201910428966.8A priority patent/CN110535452A/zh
Publication of JP2019205049A publication Critical patent/JP2019205049A/ja
Application granted granted Critical
Publication of JP6894402B2 publication Critical patent/JP6894402B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; DEAF-AID SETS; PUBLIC ADDRESS SYSTEMS
    • H04R3/00—Circuits for transducers, loudspeakers or microphones
    • H04R3/04—Circuits for transducers, loudspeakers or microphones for correcting frequency response
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00—Adaptive networks
    • H03H21/0012—Digital adaptive filters
    • H03H21/0043—Adaptive algorithms
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/048—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators using a predictor
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G10L21/0216—Noise filtering characterised by the method used for estimating noise
    • G10L21/0232—Processing in the frequency domain
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/02—Frequency selective networks
    • H03H17/0248—Filters characterised by a particular frequency response or filtering method
    • H03H17/0264—Filter sets with mutual related characteristics
    • H03H17/0266—Filter banks
    • H03H17/0269—Filter banks comprising recursive filters
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04B—TRANSMISSION
    • H04B3/00—Line transmission systems
    • H04B3/02—Details
    • H04B3/20—Reducing echo effects or singing; Opening or closing transmitting path; Conditioning for transmission in one direction or the other
    • H04B3/21—Reducing echo effects or singing; Opening or closing transmitting path; Conditioning for transmission in one direction or the other using a set of bandfilters
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04B—TRANSMISSION
    • H04B3/00—Line transmission systems
    • H04B3/02—Details
    • H04B3/20—Reducing echo effects or singing; Opening or closing transmitting path; Conditioning for transmission in one direction or the other
    • H04B3/23—Reducing echo effects or singing; Opening or closing transmitting path; Conditioning for transmission in one direction or the other using a replica of transmitted signal in the time domain, e.g. echo cancellers
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G10L2021/02082—Noise filtering the noise being echo, reverberation of the speech
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00—Adaptive networks
    • H03H21/0012—Digital adaptive filters
    • H03H21/0043—Adaptive algorithms
    • H03H2021/0049—Recursive least squares algorithm
    • H03H2021/005—Recursive least squares algorithm with forgetting factor

Landscapes

  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Multimedia (AREA)
  • Human Computer Interaction (AREA)
  • Quality & Reliability (AREA)
  • Computational Linguistics (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Physics (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Filters That Use Time-Delay Elements (AREA)
  • Cable Transmission Systems, Equalization Of Radio And Reduction Of Echo (AREA)
  • Debugging And Monitoring (AREA)
  • Signal Processing For Digital Recording And Reproducing (AREA)

Description

本発明は、システム同定装置及び方法及びプログラム及び記憶媒体に係り、特に、大規模な音響系または通信系を高速かつ数値的に安定に同定するシステム同定装置及び方法、及び、システム同定を実行させるためのプログラム及び記憶媒体に関する。
システム同定（ｓｙｓｔｅｍ ｉｄｅｎｔｉｆｉｃａｔｉｏｎ）とは、入出力データに基づいてシステムの入出力関係の数理モデル（伝達関数やインパルス応答等）を推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生およびアクティブ騒音制御、補聴器などのフィードバックキャンセラがある。従来、システム同定における適応アルゴリズムとして、ＬＭＳ（Ｌｅａｓｔ Ｍｅａｎ Ｓｑｕａｒｅ）アルゴリズム、ＲＬＳ（Ｒｅｃｕｒｓｉｖｅ Ｌｅａｓｔ Ｓｑｕａｒｅ）アルゴリズムやカルマンフィルタ（Ｋａｌｍａｎ ｆｉｌｔｅｒ）が広く用いられて来た。一般に、線形システムの入出力関係は次のように表される。
Figure 0006894402
ここで、ｕｋは入力、ｈｉはシステムのインパルス応答であり、ｙｋは観測時の雑音ｖｋを含んだ出力である。
詳細は、非特許文献１などに示されている。
１．ＬＭＳアルゴリズム
ＬＭＳは、入力ｕｋと出力ｙｋからシステムのインパルス応答ｘｋ＝［ｈ０、・・・、ｈＮ−１］Ｔを次のように推定する。
Figure 0006894402
 ただし、Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔであり、μ＞０はステップサイズと呼ばれる。また、ステップサイズμをμ／（||ＨＴ ｋ||２＋δ）で置き換えたものをＮＬＭＳアルゴリズムと呼ぶ。ここで、δは無入力でも分母が０とならないようにあらかじめ定める正の定数である。
２．ＲＬＳ
ＲＬＳは、入力ｕｋと出力ｙｋからシステムのインパルス応答ｘｋ＝［ｈ０、・・・、ｈＮ−１］Ｔを次のように推定する。
Figure 0006894402
（なお、「＾」、「ｖ」は、推定値の意味であり、数式で示すように、本来、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する場合がある。以下同様。）
３．カルマンフィルタ
状態空間モデルで表される線形システム
Figure 0006894402
の状態ｘｋの最小分散推定値ｘ＾ｋ｜ｋは、次のカルマンフィルタによって得られる。
Figure 0006894402
ｘｋ：状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。
ｙｋ：観測信号；フィルタの入力となり、既知である。
Ｈｋ：観測行列；既知である。
ｖｋ：観測雑音；未知である。
ρ：忘却係数；一般に試行錯誤で決定される。
Ｋｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式により得られる。
Σ＾０｜−１：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉが用いられる。
その他、従来、特許文献１〜３、非特許文献２〜６に示す技術がある。
特許第４０６７２６９号公報 特許第４４４４９１９号公報 特許第５１９６６５３号公報
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．
現在、システム同定において最も広く用いられている適応アルゴリズムはＬＭＳである。ＬＭＳは計算量は少ないが収束速度が非常に遅いと云った課題があった。一方、ＲＬＳやカルマンフィルタのΣ＾ｋ｜ｋ−１やＰｋは本来、正定行列（ｐｏｓｉｔｉｖｅ−ｄｅｆｉｎｉｔｅ ｓｙｍｍｅｔｒｉｃ ｍａｔｒｉｘ）であるが単精度（例：３２ｂｉｔ）で計算した場合には、しばしばこの条件を満たさず、これがカルマンフィルタやＲＬＳを数値的に不安定する要因の一つとなった。また、計算量がＯ（Ｎ２）であるため、状態ベクトルｘｋの次元（タップ数）Ｎが大きい場合、１ステップ当たりの演算回数が急激に増大し、実時間処理には必ずしも適するとはいえない場合があった。

本発明は、以上の点に鑑み、大規模な音響系または通信系を高速かつ高精度に同定することを、ひとつの目的とする。また、本発明は、このような同定する方法／装置等をフィードバックキャンセラに適用することを、他の目的とする。
本発明の第１の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜（６２）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
本発明の第２の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
本発明の第３の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
前記フィルタは、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えられ、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定装置が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
本発明の第４の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２)によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とするシステム同定方法が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
本発明の第５の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定方法が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
本発明の第６の解決手段によると、
入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
前記フィルタを、
次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与え、
前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
を特徴とする前記システム同定方法が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
本発明の第７の解決手段によると、
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、前記フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
前記処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
本発明の第８の解決手段によると、
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
本発明の第９の解決手段によると、
入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えるステップと、
前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
処理部が、評価基準として、外乱からフィルタ誤差への最大エネルギーゲインを、予め与えられた上限値γｆ ２より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
をコンピュータに実行させるためのシステム同定プログラム、及び、前記プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
Figure 0006894402
 Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

ここで、
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
Ｒｅ，ｋ：補助変数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
Ｗ：重み行列；γｆより決まる。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
γｆ：アッテネーションレベル；設計の際に与えられる。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
本発明によると、以上のように、大規模な音響系または通信系を高速かつ高精度に同定することができる。また、本発明によると、このような同定する方法／装置等をフィードバックキャンセラに適用することができる。
γｆ−σｗ ２空間の解の分布と性質についての説明図。 高速Ｈ∞フィルタ族によるｘ＾ｋ｜ｋ、ｚ＾ｋ｜ｋ−１＝Ｈｋｘ＾ｋ−１｜ｋ−１を求めるフローチャート。 スカラー存在条件のチェックについてのフローチャート。 通信系とエコーについての説明図。 エコーキャンセラの原理図。 本実施の形態（Ｃ−ＦＨＦ）と他の方法のタップ２乗誤差の推移についての説明図。 音声信号の説明図。 図７の音声信号を入力としたときの結果についての説明図。 高速Ｈ∞フィルタの処理と時間計算量の説明図（表１）。 Ｃ−高速Ｈ∞フィルタの処理と時間計算量の説明図（表２）。 Ｃ−白色化高速Ｈ∞フィルタの処理と時間計算量の説明図（表３）。 Ｃ−ＮＬＭＳアルゴリズムのフィルタの処理と時間計算量の説明図（表４）。 エコーパスのインパルス応答を示す図（表５）。 本実施の形態に関するハードウェアの構成図。 Ｃ−白色化高速Ｈ∞フィルタの最適版の一例の処理と時間計算量の説明図。 Ｃ−ＮＬＭＳアルゴリズムのフィルタの最適版の一例の処理と時間計算量の説明図。
１．概要

本発明及び／又は本実施の形態の同定法により、入出力データから未知システムを実時間で高速かつ高精度に同定することが可能となる。また、本同定法は入力信号における白色性の仮定が適切に導入できるため、アルゴリズムを効果的に簡略化でき、大幅な計算量の削減が可能である。これらは、例えば補聴器などの閉ループ系のフィードバックキャンセラに有効である。
本発明及び／又は本実施の形態によると、例えば、以下を提供することができる。
１）入出力データから未知システムを実時間で高速かつ高精度に同定する方法・装置・プログラム・記録媒体。
２）入力信号に対して白色性の仮定を導入し、計算量を大幅に削減したシステム同定法。
３）補聴器などの閉ループ系のフィードバックキャンセラ。

また、本発明及び／又は本実施の形態においては、上記課題を解決するために、入力信号の事前線形予測誤差と事後線形予測誤差の変換係数を利用して、先に提案した高速Ｈ∞フィルタ（ＦＨＦ）の計算量Ｏ（１１Ｎ）をＯ（７Ｎ）まで低減するＣ−高速Ｈ∞フィルタ（Ｃｏｎｖｅｒｓｉｏｎ Ｆａｃｔｏｒ−ｂａｓｅｄ Ｆａｓｔ Ｈ∞ Ｆｉｌｔｅｒ： Ｃ−ＦＨＦ）を導出する。さらに、本発明及び／又は本実施の形態によると入力信号に白色性の仮定を導入して、計算量をＯ（２Ｎ）まで削減するＣ−白色化高速Ｈ∞フィルタ（Ｃｏｎｖｅｒｓｉｏｎ Ｆａｃｔｏｒ−ｂａｓｅｄ Ｗｈｉｔｅｎ Ｆａｓｔ Ｈ∞ Ｆｉｌｔｅｒ： Ｃ−ＷＦＨＦ）を提供する。また、その特殊な場合としてＣ−ＮＬＭＳアルゴリズム（Ｃｏｎｖｅｒｓｉｏｎ Ｆａｃｔｏｒ−ｂａｓｅｄ Ｎｏｒｍａｌｉｚｅｄ Ｌｅａｓｔ Ｍｅａｎ Ｓｑｕａｒｅ Ａｌｇｏｒｉｔｈｍ）を与える。

本発明及び／又は本実施の形態による高速同定法は、大規模な音響系や通信系の同定に適用可能であり、スマートスピーカーやテレビ会議システムなどにおける音響エコーキャンセラやアクティブ騒音制御、音場再生などに有効である。また、本同定法は入力信号に対して白色性の仮定が適切に導入できるため、アルゴリズムを効果的に簡略化でき、大幅な計算量の削減が可能である。これらの簡易化されたアルゴリズムは補聴器や高速データ伝送の全二重中継局などの閉ループ系のフィードバックキャンセラに有効である。
２．記号の説明

まず、本実施の形態で用いる主な記号を説明する。
ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
ｘ０：初期状態；未知である。
ｗｋ：システム雑音；未知である。
ｖｋ：観測雑音；未知である。
ｙｋ：観測信号；フィルタの入力となり、既知である。
ｚｋ：出力信号；未知である。
Ｇｋ：駆動行列；実行時に既知となる。
Ｈｋ：観測行列；既知である。（Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ）
ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
ｚ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの出力ｚｋの推定値Ｈｋｘ＾ｋ｜ｋ；フィルタ方程式によって与えられる。
ｚ＾ｋ｜ｋ−１：観測信号ｙ０〜ｙｋ−１までを用いた時刻ｋの出力ｚｋの１ステップ予測値Ｈｋｘ＾ｋ−１｜ｋ−１；フィルタ方程式によって与えられる。
ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
Σ＾ｋ｜ｋ：ｘ＾ｋ｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾ｋ＋１｜ｋ：ｘ＾ｋ＋１｜ｋの誤差の共分散行列に対応；リカッチ方程式によって与えられる。
Σ＾０｜−１：初期状態の共分散行列に対応；本来未知であるが、便宜上ε０Ｉ、ε０＞０が用いられる。
Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１やゲイン行列Ｋｋから得られる。
ρ：忘却係数；γｆが決まればρ＝１−χ（γｆ）より決定される。γｆと切り離して独立に決定することも可能である。
μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
Σ０ −１：状態の不確定性を表す重み行列の逆行列；便宜上Σ０＝ε０Ｉ，ε０＞０が用いられる。

ｅｆ，ｉ：フィルタ誤差； ｅｆ，ｉ＝ｚｖ ｉ｜ｉ−Ｈｉｘｉ，ｚｖ ｉ｜ｉは出力推定値Ｈｉｘ＾i｜ｉである。
Ｒｅ，ｋ：補助変数；２×２行列である。
χｋ：変換係数；２×２行列である。
Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
γｆ：アッテネーションレベル；設計の際に与えられる。
Ｗ：重み行列；γｆより決まる。
ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
ｅｋ：事後前方予測誤差；２次元列ベクトルである。
ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
Ｏ：計算量（オーダー）；例えば、乗算回数がＮ２＋２Ｎ＋３のときＯ（Ｎ２）と表す。

なお、記号の上に付される”＾”、”ｖ”は、推定値の意味である。また、”〜”、” ”、”Ｕ”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、ｘ、ｗ、Ｈ、Ｇ，Ｋ、Ｒ、Σ、等はベクトルや行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する場合がある。
２．システム同定のハードウェア及びプログラム

本システム同定方法又はシステム同定装置・システムは、その各手順をコンピュータに実行させるためのシステム同定プログラム、システム同定プログラムを記録したコンピュータ読み取り可能な記録媒体、システム同定プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
図１４は、本実施の形態に関するハードウェアの構成図である。
このハードウェアは、中央処理装置（ＣＰＵ）である処理部１０１、入力部１０２、出力部１０３、表示部１０４、記憶部１０５、インタフェース部（Ｉ／Ｆ）１０６を有する。また、処理部１０１、入力部１０２、出力部１０３、表示部１０４、記憶部１０５、インタフェース部（Ｉ／Ｆ）１０６は、スター又はバス等の適宜の接続手段で接続されている。記憶部１０５は、システム同定される「１．記号の説明」で示した既知のデータが必要に応じて記憶され、未知・既知のデータや計算されたＣ−高速Ｈ∞フィルタ等の本実施の形態の各フィルタに関するデータ・その他のデータが処理部１０１により、必要に応じて書き込み及び／又は読み出しされる。処理部１０１は、入出力データに基づいて、本実施の形態の各フィルタを表す式に従い計算処理を実行し、システムの入出力関係の数理モデルを推定（同定）する。処理部１０１は、例えば、状態推定値、出力予測値等を求めることができる。
３．本実施の形態に関連する同定法の手段

通常のＨ∞フィルタと区別するため、モデル集合と重み集合に関してさらに最適化したＨ∞フィルタを、特にハイパーＨ∞フィルタと呼ぶことにする（非特許文献２−４）。このＨ∞フィルタは状態のダイナミックスが未知の場合でも適用できる特長を備えている。
［ハイパーＨ∞フィルタ］

状態空間モデル
Figure 0006894402
 を満たす状態推定値ｘ＾ｋ｜ｋ（あるいは出力推定値ｚｖ ｋ｜ｋ）は、次のレベルγｆのハイパーＨ∞フィルタによって与えられる。
Figure 0006894402
 であり、χ（γｆ）はχ（１）＝１、χ（∞）＝０を満たすγｆの単調減衰関数である。また、駆動行列Ｇｋは次のように生成される。
Figure 0006894402
ハイパーＨ∞フィルタの特徴は、エコーキャンセラ（状態推定）に重要な次の五つの性質、
１） 収束性、
２） エコー経路変動に対する追従能、
３） 近端話者信号に対するロバスト性、
４） 入力信号特性に関する不感度、
５） 数値的安定性
を調整できる点にある。
図１は、γｆ−σｗ ２空間の解の分布と性質についての説明図である。
図１の丸印は、ハイパーＨ∞フィルタのパラメータ空間における解を表している。それぞれの解は異なるγｆとσｗ ２の値をもつハイパーＨ∞フィルタに対応する。特に、直線γｆ＝∞上の解は忘却係数ＲＬＳアルゴリズムと一致し、点（γｆ，σｗ ２）＝（１，０）上ではステップサイズ１のＮＬＭＳアルゴリズムと一致する。よって、ハイパーＨ∞フィルタは従来の適応アルゴリズムであるＮＬＭＳやＲＬＳを含んでいることがわかる。
ハイパーＨ∞フィルタの計算量はＯ（Ｎ２）であり、このままでは実時間処理には適さない。しかし、観測行列Ｈｋがシフト特性Ｈｋ＝［ｕｋ，Ｈｋ−１（１），Ｈｋ−１（２），・・・，Ｈｋ−１（Ｎ−１）］をもつとき、計算量Ｏ（Ｎ）の高速アルゴリズムが開発されている（非特許文献２、非特許文献３）。
［高速Ｈ∞フィルタ］

観測行列Ｈｋがシフト特性を満たすとき、Σ０＝ε０Ｉ＞０をもつハイパーＨ∞フィルタは次式によって実行できる。
Figure 0006894402
 であり、ｃｋ∈Ｒ２×１はＣｋ＝［ｃｋ，・・・，ｃｋ−Ｎ＋１］の第一列ベクトルであり、ｃｋ−ｉ＝０２×１、ｋ−ｉ＜０と仮定し、初期値はＫ−１＝０Ｎ×２、ａ−１＝０Ｎ×１、Ｓ−１＝１／ε０、ε０＞０、ｂ−１＝０Ｎ×１、ｘ＾−１｜−１＝０Ｎ×１と設定した。ここで、０ｍ×ｎはｍ×ｎ零行列を表す。
このとき、次のスカラー存在条件を満たさなければならない。
Figure 0006894402
４．本実施の形態の同定法の手段

４．１ 準備
高速Ｈ∞フィルタの新しいアルゴリズムを導出するためにいくつかの命題を示す。その命題の証明には、次のように定義された（Ｎ＋１）×（Ｎ＋１）行列ＱＵ ｋ が利用される。
Figure 0006894402
 この行列ＱＵ ｋは、観測行列ＣＵ ｉのシフト特性から次のように二通りに分解できる。
Figure 0006894402
 その分解されたＱＵ ｋは次式によって再帰的に求められる。
Figure 0006894402
命題１
ゲイン行列Ｋｋ：＝Ｑｋ −１Ｃｋ Ｔは疑似ゲイン行列Ｋ〜 ｋ：＝Ｑｋ−１ −１Ｃｋ Ｔを用いて次のように表すことができる。
Figure 0006894402
 命題２
後方線形予測係数ｂｋ：＝−Ｑｋ −１ ｔ ｋの新たな更新式は次式で与えられる。
Figure 0006894402
 命題３
事後後方予測誤差ｅ ｋは事前後方予測誤差ｅ 〜 ｋを用いて次のように表すことができる。
Figure 0006894402
 命題４
事後前方予測誤差ｅｋは事前前方予測誤差ｅ〜 ｋを用いて次のように表すことができる。
Figure 0006894402
命題５
後方線形予測係数ｂｋは事後後方予測誤差ｅ ｋを用いて次のように再帰的に更新できる。
Figure 0006894402
 命題６
前方線形予測係数ａ：＝−Ｑｋ−１ −１ ｔ ｋは事後前方予測誤差ｅｋを用いて次式によって再帰的に更新できる。
Figure 0006894402
 命題７
疑似ゲイン行列Ｋ〜 ｋを計算するもう一つの形式は次式で与えられる。
Figure 0006894402
 命題８
事前後方予測誤差ｅ 〜 ｋはμ〜 ｋを用いて次のように表すことができる。
Figure 0006894402
 命題９
式（３６）の２×２行列Ｒｅ，ｋは次式で再帰的に更新できる。
Figure 0006894402
４．２ アルゴリズムの導出
命題１から命題９を用いれば、前方線形予測係数ａｋ と後方線形予測係数ｂｋ が事後前方予測誤差ｅｋと事後後方予測誤差ｅ ｋを用いて更新される新たな高速Ｈ∞フィルタ（Ｃ−高速Ｈ∞フィルタ）が次のように得られる。
Figure 0006894402
ただし、各パラメータの設定や初期化は高速Ｈ∞フィルタの場合と同じである。このＣ−高速Ｈ∞フィルタは、変換係数（ｃｏｎｖｅｒｓｉｏｎ ｆａｃｔｏｒ）χｋを用いた高速Ｈ∞フィルタの事後予測誤差形式に当たる。
４．３ 計算量の評価

Ｃ−高速Ｈ∞フィルタは高速Ｈ∞フィルタと理論的に等価であるが、計算量は少なくて済む。Ｃ−白色化高速Ｈ∞フィルタは、入力信号が白色であると仮定して、アルゴリズムを大幅に縮退させたＣ−高速Ｈ∞フィルタである。Ｃ−ＮＬＭＳアルゴリズムはパラメータが特定の値であるときのＣ−白色化高速Ｈ∞フィルタに対応し、理論的にはＮＬＭＳアルゴリズムと等価であるが、計算量はＯ（２Ｎ）で済む。

図９に、高速Ｈ∞フィルタの処理と時間計算量の説明図（表１）を示す。
高速Ｈ∞フィルタの各式の計算に必要な乗算と加算の回数を図９（表１）にまとめた（特許文献１参照、比較のために記載）。高速Ｈ∞フィルタを実行するためには、図２のＣｏｒｅ Ａｌｇｏｒｉｔｈｍで表の数式を上から順に（式番号が大きい番号から小さい番号へ降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
図１０に、Ｃ−高速Ｈ∞フィルタの処理と時間計算量の説明図（表２）を示す。
図１０（表２）は、Ｃ−高速Ｈ∞フィルタに対する乗算と加算の回数を示している。ここで、演算回数の算出に当たり、Ｋｋ，Ｃｋ，ｅ〜 ｋ，ｅｋ，およびＷの構造が利用された。例えば、Ｋｋ−１Ｗｅ〜 ｋは、Ｋｋ−１（：，１）＝Ｋｋ−１（：，２）とｅ〜 ｋ（１，１）＝ｅ〜 ｋ（２，１）を用いて、（１−γｆ −２）Ｋｋ−１（：，１）ｅ〜 ｋ（１，１）に置き換えられた。同様に、Ｃｋ−１とａｋ−１の積を求めるに当たり、Ｃｋ−１（２，：）ａｋ−１はＣｋ−１（１，：）ａｋ−１で置換された。同じ組合せの２項の乗算が複数ある場合は最初の計算のみをカウントした。χｋ−１ −１は前のステップで計算したχｋ −１を用いた。スカラーによるベクトルの除算は、スカラーの逆数によるベクトルの乗算として処理した。加算回数の場合も同様な方法で算出した。
Ｃ−高速Ｈ∞フィルタは、Ｎが十分に大きい場合、単位サンプル当たり７Ｎ回の乗算と７Ｎ回の加算のみで実行できる。一方、高速Ｈ∞フィルタは、単位サンプル当たり１１Ｎ回の乗算と１０Ｎ回の加算が必要である。両フィルタとも単位サンプル当たり３回の乗算が必要である。また、いずれの高速Ｈ∞フィルタも空間計算量(メモリ)はＯ（６Ｎ）である。
４．４ 高速Ｈ∞フィルタ族の処理

図２に、高速Ｈ∞フィルタ族によるｘ＾ｋ｜ｋ、ｚ＾ｋ｜ｋ−１＝Ｈｋｘ＾ｋ−１｜ｋ−１を求めるフローチャートを示す。ここで、高速Ｈ∞フィルタ族とは、例えば、Ｃ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）、Ｃ−白色化高速Ｈ∞フィルタ（Ｃ−ＷＦＨＦ）、Ｃ−ＮＬＭＳアルゴリズム・フィルタ等を示す。高速Ｈ∞フィルタとＣ−高速Ｈ∞フィルタでは、このフローチャートに、後述のように、図３の存在条件のチェック機能が追加される場合がある。また、必要に応じてｚ＾ｋ｜ｋを求める。
以下に、高速Ｈ∞フィルタ族の処理について、フローチャートを参照して説明する。
［ステップＳ１０１、初期化］処理部１０１は、高速Ｈ∞フィルタの再帰変数を初期化する。たとえば、処理部１０１は、記憶部１０５から再帰式の初期状態を読み出し、又は初期状態を入力部１０２又はＩ／Ｆ１０６から入力し、図示のように定める。
［ステップＳ１０３］処理部１０１は、事前に決定したデータ長Ｌを超えているか判断する。たとえば処理部１０１は、時刻ｋとデータ数Ｌと比較する。なお、Ｌはあらかじめ定められた最大データ数を表すことができる。処理部１０１は、時刻ｋがデータ数Ｌ以下であれば次のステップＳ１０５に進む。処理部１０１は、時刻ｋが超えていればフィルタゲインを初期化し、フィルタをリスタートするか、あるいは終了する（Ｓ１１３）。（なお、不要であれば条件文を取り除くことができる。または、必要に応じて再スタートを行ってもよい。）
［ステップＳ１０５、入力］処理部１０１は、時刻ｋがＬ以下であれば、拡大観測行列ＣＵ ｋを入力信号とシフト特性を利用して更新する。たとえば、処理部１０１は、入力信号ｕｋを入力部１０２又はＩ／Ｆ１０６から入力し、ＣＵ ｋを図示のように設定する。なお、入力信号ｕｋを記憶部１０５から読み出すようにしてもよい。
ここで、入力信号ｕｋからなる２×（Ｎ＋１）拡大観測行列
Figure 0006894402
 ［ステップＳ１０７、Ｃｏｒｅ Ａｌｇｏｒｉｔｈｍ］処理部１０１は、Ｃｏｒｅ Ａｌｇｏｒｉｔｈｍでは図９（表１）から図１２（表４）及び図１５、図１６のいずれかのアルゴリズムを実行する（詳細は後述）。なお、これらのいずれかのアルゴリズムにおいて、状態推定値ｘ＾ｋ｜ｋを緩和係数（ステップサイズ）μ（ここで、μは、μ＞０の実数）で更新するようにしてもよい。この場合、例えば、式（４９）の右辺第２項にμを乗算することで、緩和係数（ステップサイズ）μで更新することができる。また、必要に応じてアルゴリズムの実行後又は実行前等の適宜のタイミングで、図３の存在条件の判定のための処理を挿入するようにしてもよい（詳細は後述）。
［ステップＳ１０９、推定値］処理部１０１は、ステップＳ１０７により計算された状態推定値ｘ＾ｋ｜ｋ（あるいは出力推定値ｚ＾ｋ｜ｋ）を記憶部１０５に記憶、及び／又は出力部１０３若しくはＩ／Ｆ１０６から出力する。
［ステップＳ１１１］処理部１０１は、時刻ｋを進ませて（ｋ＝ｋ＋１）、ステップＳ１０３に戻り、データがある限り又は所定回数繰り返し続ける。
なお、処理部１０１は、ステップＳ１０５〜Ｓ１０９等の各ステップで求めた適宜の中間値及び最終値、存在条件の値等を必要に応じて適宜記憶部１０５に記憶し、また、記憶部１０５から読み出すようにしてもよい。さらに、処理部１０１は、入力部１０２若しくはＩ／Ｆ１０６からの指示により又は予め定められた処理により、各値によるデータ又はグラフ・表・図等を作成して、記憶部１０５に記憶、及び／又は出力部１０３若しくはＩ／Ｆ１０６から出力するようにしてもよい。
図３に、スカラー存在条件のチェックについてのフローチャートを示す。
図２のフローチャートにさらに図３の存在条件のチェック機能が追加してもよい。たとえば、処理部１０１は、ステップＳ１０７の処理の実行後又は実行前等の適宜のタイミングで、図示のスカラー存在条件の判定ステップを実行し、式（３０）を満たさなければ処理を終了する又はリスタートするようにしてもよい。なお、このステップは、他にも、図２のフローチャートの適宜のひとつ又は複数のステップや、時刻ｋを進ませるステップ（Ｓ１１１）の前後等の適宜のタイミングで、処理部１０１が実行することができる。なお、図中、ＥＸＣ（ｋ）は、前述の式（３０）の左辺に相当する。
４．５ 白色性仮定の導入
図１１に、緩和係数μを導入したＣ−白色化高速Ｈ∞フィルタの処理と時間計算量の説明図（表３）を示す。
入力信号が白色雑音と仮定すれば、線形予測係数ａｋとｂｋは０となる。このとき、ｅ〜 ｋ＝ｃｋ＋Ｃｋ−１ａｋ−１＝ｃｋ、ｅｋ＝ｃｋ＋Ｃｋ−１ａｋ＝ｃｋ、ｅ 〜 ｋ＝ｃｋ−Ｎ＋Ｃｋｂｋ−１＝ｃｋ−Ｎ、ｅ ｋ＝ｃｋ−Ｎ＋Ｃｋｂｋ＝ｃｋ−Ｎ、Ｋ〜 ｋ＝ｍ〜 ｋ、Ｓ ｋは不要となり、式（４９）、式（５０）、式（５５）、式（５８）、式（６０）、式（６１）、式（６２）は次のように整理できる。ただし、μは緩和係数(ステップサイズ）で、μ＞０の実数であり、事前に与えられる。
なお、ｅ〜 ｋの式（５６）、ｅ ｋの式（５３）は、本来の定義式とは別に、理論的に等価な表現が複数あり、上述のように表すことができる。また、Ｃ−ＷＦＨＦの計算で、後方予測関係の変数（アンダーバーが付いたもの）は、利用しないことができるので、それらの計算は不要となる。
Figure 0006894402
これより、Ｃ−高速Ｈ∞フィルタは表３のＣ−白色化高速Ｈ∞フィルタに縮退する。ただし、ｚ＾ｋ｜ｋ＝Ｈｋｘ＾ｋ｜ｋは省略した。
Ｃ−白色化高速Ｈ∞フィルタを実行するためには、図２のＣｏｒｅ Ａｌｇｏｒｉｔｈｍで図１１（表３）の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
４．６ Ｃ−ＮＬＭＳアルゴリズム
図１２に、Ｃ−ＮＬＭＳアルゴリズムのフィルタの処理と時間計算量の説明図（表４）を示す。
γｆ＝１．０、ρ＝１．０のとき、Ｓｋ＝Ｓｋ−１＝Ｓ−１＝１／ε０となる。これより、Ｃ−白色化高速Ｈ∞フィルタは次のように縮退する。
Figure 0006894402
これはステップサイズμ＝１のＮＬＭＳアルゴリズムと理論的に等価である。よって、緩和係数μ（μ＞０の実数）を式（４９）に導入すれば、表４のステップサイズμのＣ−ＮＬＭＳアルゴリズムが得られる。このアルゴリズムの計算量はＯ（２Ｎ）であり、Ｓ−１＝１／ε０はＮＬＭＳアルゴリズムの正則化項δに一致する。

Ｃ−ＮＬＭＳアルゴリズムのフィルタを実行するためには、図２のＣｏｒｅ Ａｌｇｏｒｉｔｈｍで図１２（表４）の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
４．７ Ｃ−高速Ｈ∞フィルタの白色化の妥当性

高速Ｈ∞フィルタは入力信号の線形予測係数を０に置くことによって、次のように白色化することができる。
Figure 0006894402
 ここで、Ｋｋ−１（１：Ｎ−１，：）は、Ｋｋ−１の第１行〜第（Ｎ−１）行までからなる行列を表す。

この白色化高速Ｈ∞フィルタ（ＷＦＨＦ）は、γｆ＝∞のとき、次のように縮退する。
Figure 0006894402
 このとき、フィルタゲインの分母の正規化項 １＋γｆ −２ＨｋＫｋ（：，１） は、１となり、入力信号ｕｋの振幅が急激に増大するとフィルタゲインが大きくなり過ぎて不安定となる場合が想定される（γｆ＜∞の場合は、分母にＨｋＫｋ（：，１）があるのでゲインレベルが制限される）。これより、従来の高速Ｈ∞フィルタの白色化に課題があることがわかる。また、ＷＦＨＦの計算量はＯ（３Ｎ）であった。

この課題を解決して、白色化を適切に実行することができるようにしたのが、本発明及び／又は本実施の形態のＣ−ＷＦＨＦである。白色化のために、例えば、数多くあるパラメータの中から、線形予測係数ａｋとｂｋを選択して０とし、また、ｅ〜 ｋ、ｅｋ、ｅ 〜 ｋＮ、ｅ ｋ、の多数の数式から適切なものをあてはめ、さらに、Ｋ〜 ｋ＝ｍ〜 ｋ、という式を適用し、また、Ｓ ｋは不要となるようにした。このように、数多く存在するパラメータの中から特定のものを選択して特定の値とし、数多く存在する式の中から特定のパラメータについて特定の式を適用することは、極めて高度かつ困難なことであって、これにより本発明および／本実施の形態は、白色化を適切に行うことを達成したものである。
４．８ Ｃ−白色化高速Ｈ∞フィルタ（Ｃ−ＷＦＨＦ）、Ｃ−ＮＬＭＳアルゴリズム（Ｃ−ＮＬＭＳ）の最適版の検討

Ｃ−ＷＦＨＦ、Ｃ−ＮＬＭＳの実装を考慮した場合における計算式の最適化の一例及び演算量の見積もりを行なった。ここで、最適化とは、特に、「組み込みのための最適化」をいう。ここでは、例えば、組み込まれたプログラムの演算回数が最小になる式（形式）を求めることを、最適化と呼んでいる。なお、「Ｃ−ＦＨＦとの関係性を維持した形式」と「組み込みに適した形式（最適版）」は、理論的に等価である。例えば、「Ｃ−ＦＨＦとの関係性を維持した形式」はそのアルゴリズムの意味がわかり、一方、「組み込みに適した形式」は実際にＣＰＵにプログラムするのに適した形式である。
ＤＳＰ（digital signal processor）等のプロセッサやハードウェアに組み込む際には掛け算（Ｍｕｌ．）、割り算（Ｄｉｖ．）、足し算（Ａｄｄ．）によって必要となるサイクル数が異なるため、それぞれの演算回数からサイクル数の見積もりを行なった。実際にはメモリやレジストリからの読み込みや並列化演算などによって必要なサイクル数は大きく変化するが組み込み方法によるサイクル数の変化は考慮せず算出した。

図１５に、Ｃ−白色化高速Ｈ∞フィルタの最適版の一例の処理と時間計算量の説明図を示す。
このとき、Ｃ−白色化高速Ｈ∞フィルタの式（６０）、式（５８）、式（５５）、式（４９）、式（５０）は、次のように整理できる。ここで、Ｋ~ ｋ−１（１：Ｎ−１，１）はＫ~ ｋ−１の第１行〜第Ｎ−１行までの第１列の行列を表す。
なお、表３のＣ−白色化高速Ｈ∞フィルタ（Ｃ−ＷＦＨＦ）の式にμ（ステップサイズ、μ＞０の実数）を付けて、Ｃ−ＷＦＨＦがＣ−ＮＬＭＳを完全に含むように表現することができる。よって、Ｃ−ＮＬＭＳは、γｆ＝１，ρ＝１のときのＣ−ＷＦＨＦと完全に一致する。
ここで、最適化の式とは演算回数が少なくなるように式を組み換えものであり、理論的に等価になる。例えば、Ｃ−ＷＦＨＦの式（５８）の第１列の第１成分から第Ｎ成分を抽出して利用したのが、以下のＫ〜 ｋ（：，１）の式である。

Figure 0006894402
なお、Ｃ−ＷＦＨＦの式（５８）及び（６１）（６２））から、
Ｋ〜 ｋ−１（Ｎ，１）＝ｕｋ−Ｎ／Ｓｋ−Ｎ−１
という関係式が導出され、これによりＣ−ＷＦＨＦの式（５５）から上述のＲｅ，ｋ（１，１）の式が導かれる。

よって、実際の計算に必要な演算はＲｅ，ｋの１行１列目のスカラー値であることから、Ｃ−白色化高速Ｈ∞フィルタの式（５５）は上述した式のように省略することができる。

Ｃ−白色化高速Ｈ∞フィルタを実行するためには、図２のＣｏｒｅ Ａｌｇｏｒｉｔｈｍで図表の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ｄｉｖ．」は割算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
このように、この最適化によると、先に示したＣ−白色化高速Ｈ∞フィルタより計算量を一層削減することができる。
図１６に、Ｃ−ＮＬＭＳアルゴリズムのフィルタの最適版の一例の処理と時間計算量の説明図を示す。
このとき、Ｃ−ＮＬＭＳアルゴリズムのフィルタの式（５８）、式（５５）、式（４９）、式（５０）は、次のように整理できる。なお、例えば、Ｃ−ＮＬＭＳの式（５８）の第１列の第１成分から第Ｎ成分を抽出して利用したのが、以下のＫ〜 ｋ（：，１）の式である。

Figure 0006894402
 （μはステップサイズ、μ＞０の実数）

なお、Ｃ−ＮＬＭＳの式（５８）から、
Ｋ〜 ｋ−１（Ｎ，１）＝ε０ｕｋ−Ｎ
という関係式が導出され、これによりＣ−ＮＬＭＳの式（５５）から上述のＲｅ，ｋ（１，１）の式が導かれる。

よって、実際の計算に必要な演算はＲｅ，ｋの１行１列目のスカラー値であることから、Ｃ−ＮＬＭＳアルゴリズムのフィルタの式（５５）は上述した式のように省略することができる。

Ｃ−ＮＬＭＳアルゴリズムのフィルタを実行するためには、図２のＣｏｒｅ Ａｌｇｏｒｉｔｈｍで図表の数式を上から順に（式番号が降順に）計算することにより、一サンプル（単位時刻）当たりのフィルタが実行できる。この式の順序は一般的なフィルタの記述とは逆となっている。なお、各表中、「Ｍｕｌ．」は積算回数、「Ｄｉｖ．」は割算回数、「Ａｄｄ．」は加算回数を、それぞれ示す。
このように、この最適化の一例によると、先に示したＣ−ＮＬＭＳアルゴリズムのフィルタより計算量を一層削減することができる。これはδ＝１／ε０のＮＬＭＳアルゴリズムと一致する。
５．エコーキャンセラ
５．１ 準備
以下に、実施の形態として、エコーキャンセラの例を取り上げ、本同定アルゴリズムの動作を確認する。その前に、通信エコーキャンセラについて簡単に説明する。最近注目されているスマートスピーカーなどの音声インターフェイスに不可欠な音響エコーキャンセラも原理は同じである。

図４に、通信系とエコーについての説明図を示す。
国際電話など長距離電話回線では、信号増幅などの理由から４線式回線が用いられている。一方、加入者回線は比較的短距離なので２線式回線が使用されている。２線式回線と４線式回線の接続部には図４のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Ｂからの信号（音声）は話者Ａのみに到達する。しかし、一般には完全に整合をとることはできず、受信信号の一部は４線式回線に漏れ、増幅された後、再び送信者（話者Ｂ）に戻る現象が起こる。これがエコー（ｅｃｈｏ）である。エコーは伝送距離が長くなるにつれて影響が大きくなり、著しく通話の品質を劣化させることが想定される（パルス伝送においては近距離であってもエコーによる通話品質の劣化は問題となることが想定される）。
図５に、エコーキャンセラの原理図を示す。
そこで、図５のようにエコーキャンセラ（ｅｃｈｏ ｃａｎｃｅｌｌｅｒ）を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、その結果を利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
エコーパスのインパルス応答の推定は、残留エコーｅｋの平均２乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Ａからの信号（音声）である。一般に、話者２人が同時に話し始めた（ダブルトーク）ときはインパルス応答の推定を中断する。
次に、このエコーキャンセリング問題の数理モデルを作成する。まず、受信信号ｕｋがエコーパスへの入力信号となることを考慮すれば、エコーパスのインパルス応答ｈｉ［ｋ］により、エコーｄｋの観測値ｙｋは次式で表される。
Figure 0006894402
 ここで、ｕｋ、ｙｋはそれぞれ時刻ｋにおける受信信号とエコーの観測値を表し、ｖｋは時刻ｋにおける平均値０の回線雑音とし、タップ数Ｎは既知とする。このとき、インパルス応答の推定値ｈ＾ｉ［ｋ］が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
Figure 0006894402
 これをエコーから差し引けば（ｙｋ−ｄ＾ｋ≒０）、エコーをキャンセルすることができる。ただし、ｋ−ｉ＜０のときｕｋ−ｉ＝０とする。
以上より、キャンセリング問題は直接観測可能な受信信号ｕｋとエコーの観測値ｙｋからエコーパスのインパルス応答ｈｉ［ｋ］を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにＨ∞フィルタを適用するには、まず式（６５）を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、｛ｈｉ［ｋ］｝Ｎ−１ ｉ＝ｏを状態変数ｘｋとし、ｗｋ程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
Figure 0006894402
５．２ 動作の確認
図１３は、エコーパスのインパルス応答を示す図（表５）である。
次に、エコーパスの同定問題を用いて、Ｃ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）の妥当性を確認すると共に、高速カルマンフィルタ（ＦＫＦ）と高速トランスバーサルフィルタ（ＦＴＦ）に対する優位性を示す。未知システムは長さＮ＝１２８をもつインパルス応答｛ｈｉ（ｋ）｝Ｎ−１ ｉ＝０で表され、時刻ｋ＝４０００で一度シフトし、その後はその値を保持した。
このとき、エコーの観測値ｙｋは次式で表される。
Figure 0006894402
 ただし、｛ｈｉ｝２３ ｉ＝０は表５の値を採用し、その他｛ｈｉ｝１２７ ｉ＝２４は０とした。また、ｖｋは平均値０、標準偏差σｖ＝１０−４０の定常なガウス白色雑音とした。
また、受信信号ｕｋは次のように２次のＡＲモデルで近似した。
ｕｋ＝α１ｕｋ−１＋α２ｕｋ−２＋ｗ’ｋ (71)
ただし、α１＝０．７、α２＝０．１とし、ｗ’ｋは平均値０、標準偏差σｗ’＝０．０２の定常なガウス白色雑音とした。
図６に、本実施の形態（Ｃ−ＦＨＦ）と他の方法のタップ２乗誤差の推移についての説明図を示す。
図６は、Ｃ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）のタップ２乗誤差||ｘ＾ｋ｜ｋ−ｘｋ||２の推移を、高速Ｈ∞フィルタ（ＦＨＦ）、高速カルマンフィルタ（ＦＫＦ）、高速トランスバーサルフィルタ（ＦＴＦ）と比較した結果を示す。１．３＜γｆ＜∞のとき、計算量の違いはあるが、高速Ｈ∞フィルタ（ＦＨＦ）とＣ−高速Ｈ∞フィルタ（Ｃ−ＦＨＦ）の性能の差はほとんど見られない。興味深いことに、このように観測雑音が非常に小さい場合、高速Ｈ∞フィルタとＣ−高速Ｈ∞フィルタは、高速カルマンフィルタと高速トランスバーサルフィルタと比べて高い収束性能と追従性能を発揮した。

次に、図７に、音声信号の説明図を示す。
また、図８に、図７の音声信号を入力としたときの結果についての説明図を示す。
この場合も、Ｃ−高速Ｈ∞フィルタは高速カルマンフィルタや高速トランスバーサルフィルタより優れた収束性能と追従性能を示した。
６． 命題の証明

以下に、上述の各命題の証明について説明する。

命題１の証明
Figure 0006894402
 命題２の証明
Figure 0006894402
命題３の証明
命題１と命題２を用いれば、事後後方予測誤差ｅ ｋは次のように事前後方予測誤差ｅ 〜 ｋを用いて表すことができる。
Figure 0006894402
 命題４の証明
命題１を用いれば、事後前方予測誤差ｅｋは次のように事前前方予測誤差ｅ 〜 ｋを用いて表すことができる。
Figure 0006894402
 命題５の証明
命題１、命題２および命題３を用いれば、後方線形予測係数ｂｋの更新式は次のように書き換えることができる。
Figure 0006894402
 命題６の証明
命題１と命題４を式（２９）に適用すれば、前方線形予測係数ａｋの新たな更新式が得られる。
命題７の証明
（Ｎ＋１）×（Ｎ＋１）行列ＱＵ ｋは次のように二通りに分解できる。
Figure 0006894402
 逆行列の補助定理を上記のように分解されたＱＵ ｋ−１に適用すれば次式が得られる。
Figure 0006894402
 式（７９）の両側に右からＣＵ ｋ Ｔを掛ければ、次のように整理できる。
Figure 0006894402
 このとき、式（８３）の両辺の各項を比較すれば次式が得られる。
Figure 0006894402
命題８の証明
式（８３）の両辺の第２項からμ〜 ｋ＝ｅ 〜Ｔ ｋ／Ｓ ｋ−１が得られ、これより式（４６）が導かれる。

命題９の証明
式（７９）の両辺に両側からＣＵ ｋとＣＵ ｋ Ｔを掛ければ次式が得られる。
Figure 0006894402
 これに命題８を適用すれば証明は完結する。
１０１ 処理部
１０２ 入力部
１０３ 出力部
１０４ 表示部
１０５ 記憶部
１０６ インタフェース部

Claims (21)

  1. 入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
    評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
    前記フィルタは、
    次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜（６２）によって与えられ、
    前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
    を特徴とする前記システム同定装置。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    χｋ：変換係数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
    ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  2. 入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
    評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
    前記フィルタは、
    次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えられ、
    前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
    を特徴とする前記システム同定装置。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  3. 請求項２に記載されたシステム同定装置において、
    式（６０）を用い、式（５８）、式（５５）、式（４９）及び式（５０）について、式（５８）、式（５５）、式（４９）及び式（５０）を変形した以下の式を用いることを特徴とするシステム同定装置。
    Figure 0006894402
     ここで、μはステップサイズ、μ＞０の実数である。
  4. 請求項２に記載されたシステム同定装置において、
    前記フィルタは、γｆ＝１．０、ρ＝１．０、Ｓｋ＝Ｓｋ−１＝Ｓ−１＝１／ε０とすることにより、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）ように縮退した式を用いることを特徴とする前記システム同定装置。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ
    ここで、
    ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
  5. 入出力データからシステムの高速実時間同定を行うシステム同定装置であって、
    評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを備え、
    前記フィルタは、
    次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えられ、
    前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
    を特徴とする前記システム同定装置。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  6. 請求項５に記載のシステム同定装置において、
    式（５８）、式（５５）、式（４９）及び式（５０）について、式（５８）、式（５５）、式（４９）及び式（５０）を変形した以下の式を用いることを特徴とするシステム同定装置。
    Figure 0006894402
    ここで、μはステップサイズ、μ＞０の実数である。
  7. 請求項１乃至３のいずれかに記載のシステム同定装置において、
    前記フィルタを実行する処理部
    を備え、
    前記処理部は、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化する第１ステップと、
    前記処理部は、時刻ｋが予め定められたデータ長Ｌを超えているか判断し、超えていればフィルタゲインを初期化してリスタートする、又は、処理を終了する、第２ステップと、
    前記処理部は、時刻ｋがデータ長Ｌ以下であれば、入力信号ｕｋを入力部若しくはインタフェース部から入力し又は記憶部から読み出し、拡大観測行列ＣＵ ｋを入力信号ｕｋとシフト特性ＣＵ ｋ＝［ｃｋ，Ｃｋ−１］＝［Ｃｋ，ｃｋ−Ｎ］、を利用して更新する第３ステップと、
    前記処理部は、フィルタを与える式を式番号が式（６２）から式（４９）の降順に、式（４９）の右辺第２項を緩和係数μ（ここで、μは、μ＞０の実数）で計算する処理を含み、単位時間あたりのフィルタ・アルゴリズムを実行する第４ステップと、
    前記処理部は、式（４９）を含むフィルタ方程式を更新し、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力する第５ステップと、
    前記処理部は、時刻ｋを進ませて、前記第１ステップ乃至前記第４ステップの処理を繰り返す第６ステップと、
    を含むことを特徴とするシステム同定装置。
  8. 請求項４又は５又は６に記載のシステム同定装置において、
    前記フィルタを実行する処理部
    を備え、
    前記処理部は、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化する第１ステップと、
    前記処理部は、時刻ｋが予め定められたデータ長Ｌを超えているか判断し、超えていればフィルタゲインを初期化してリスタートする、又は、処理を終了する、第２ステップと、
    前記処理部は、時刻ｋがデータ長Ｌ以下であれば、入力信号ｕｋを入力部若しくはインタフェース部から入力し又は記憶部から読み出し、拡大観測行列ＣＵ ｋを入力信号ｕｋとシフト特性ＣＵ ｋ＝［ｃｋ，Ｃｋ−１］＝［Ｃｋ，ｃｋ−Ｎ］、を利用して更新する第３ステップと、
    前記処理部は、フィルタを与える式を式番号が式（５８）から式（４９）の降順に、式（４９）をステップサイズμ（ここで、μは、μ＞０の実数）で計算する処理を含み、単位ステップサイズあたりのフィルタ・アルゴリズムを実行する第４ステップと、
    前記処理部は、式（４９）を含むフィルタ方程式を更新し、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力する第５ステップと、
    前記処理部は、時刻ｋを進ませて、前記第１ステップ乃至前記第４ステップの処理を繰り返す第６ステップと、
    を含むことを特徴とするシステム同定装置。
  9. 請求項１乃至８のいずれかに記載のシステム同定装置において、
    前記評価基準は、前記状態空間モデルに対して、次式（１４）を満たすことを特徴とするシステム同定装置。
    Figure 0006894402
     ここで、
    ｅｆ，ｉ：フィルタ誤差；ｅｆ，ｉ＝ｚｖ ｉ｜ｉ−Ｈｉｘｉ，ｚｖ ｉ｜ｉは出力推定値Ｈｉｘ＾i｜ｉである。

    Σ０ −１：状態の不確定性を表す重み行列の逆行列；Σ０は既知である。
    ρ：忘却係数；γｆが決まればρ＝１−χ（γｆ）より決定される。γｆと切り離して独立に決定することも可能である。
    γｆ：アッテネーションレベル；設計の際に与えられる。
  10. 請求項１乃至９のいずれかに記載のシステム同定装置において、
    前記フィルタの処理部は、前記第３ステップにおいて、処理の実行後又は実行前に、次式のスカラー存在条件の判定ステップを実行し、式（３０）を満たさなければ処理を終了する又はリスタートすることを特徴とするシステム同定装置。
    Figure 0006894402
  11. 請求項１乃至１０のいずれかに記載のシステム同定装置において、
    前記フィルタの処理部は、前記第５ステップにおいて、さらに次式（４８）により出力推定値ｚｖ ｋ｜ｋ，ｚｖ ｋ｜ｋ−１を計算し、記憶部に記憶及び／又は出力部若しくはインタフェース部から出力することを特徴とするシステム同定装置。
    ｚ＾ｋ｜ｋ＝Ｈｋｘ＾ｋ｜ｋ (48)

    ここで、
    ｚ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの出力ｚｋの推定値Ｈｋｘ＾ｋ｜ｋ；フィルタ方程式によって与えられる。
    ｚ＾ｋ｜ｋ−１：観測信号ｙ０〜ｙｋ−１までを用いた時刻ｋの出力ｚｋの１ステップ予測値Ｈｋｘ＾ｋ−１｜ｋ−１；フィルタ方程式によって与えられる。
  12. 請求項１乃至１１のいずれかに記載のシステム同定装置において、
    前記フィルタを適用し、状態推定値ｘ＾ｋ｜ｋを求め、
    擬似エコーｄ＾ｋを次式（６６）のように推定し、これをエコーから差し引くことにより、エコーをキャンセルすることを特徴とするシステム同定装置。
    Figure 0006894402
     ここで、
    ［ｈ＾０［ｋ］、・・・、ｈ＾Ｎ−１［ｋ］］Ｔはｘ＾ｋ−１｜ｋ−１乃至ｘ＾ｋ｜ｋとし、
    ｕｋは時刻ｔｋ（＝ｋＴ；Ｔは標本化周期）における受信信号、
    Ｎはタップ数（既知）とする。
  13. 入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
    評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
    前記フィルタを、
    次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２)によって与え、
    前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
    を特徴とするシステム同定方法。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    χｋ：変換係数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
    ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  14. 入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
    評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
    前記フィルタを、
    次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与え、
    前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
    を特徴とする前記システム同定方法。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  15. 入出力データからシステムの高速実時間同定を行うシステム同定方法であって、
    評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタを用い、
    前記フィルタを、
    次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与え、
    前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定すること、
    を特徴とする前記システム同定方法。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  16. 入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムであって、
    処理部が、前記フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２）によって与えるステップと、
    前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
    前記処理部が、評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
    前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
    をコンピュータに実行させるためのシステム同定プログラム。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    χｋ：変換係数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
    ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  17. 入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムであって、
    処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えるステップと、
    前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
    前記処理部が、評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
    前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
    をコンピュータに実行させるためのシステム同定プログラム。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  18. 入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムであって、
    処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えるステップと、
    前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
    前記処理部が、評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
    前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
    をコンピュータに実行させるためのシステム同定プログラム。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  19. 入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
    処理部が、前記フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）〜(６２）によって与えるステップと、
    前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
    前記処理部が、評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
    前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
    をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    χｋ：変換係数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    ｍｋ：ＫＵ ｋの第１〜第Ｎ行からなるＮ×２行列；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ａｋ：前方線形予測係数；Ｎ次元列ベクトルである。
    ｂｋ：後方線形予測係数；Ｎ次元列ベクトルである。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  20. 入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
    処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）、（６０）〜（６２）によって与えるステップと、
    前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
    前記処理部が、評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
    前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
    をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
  21. 入出力データからシステムの高速実時間同定をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
    処理部が、フィルタを、次式（１１）〜（１３）で表される状態空間モデルに対して、次式（４９）、（５０）、（５５）、（５８）によって与えるステップと、
    前記処理部が、記憶部から再帰式の初期状態を読み出し、又は、初期状態を入力部若しくはインタフェース部から入力し、前記フィルタの再帰変数を初期化するステップと、
    前記処理部が、評価基準として、システム雑音ｗ ｋ 及び観測雑音ｖ ｋ を含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γ ｆ に対応する項より小さく抑えるように定めることにより、外乱に対して頑健なフィルタとし、前記フィルタを与える前記式を降順に計算することにより、前記システムの状態ｘｋの状態推定値ｘ＾ｋ｜ｋを推定するステップと、
    前記処理部が、状態推定値ｘ＾ｋ｜ｋを記憶部に記憶する及び／又は出力部若しくはインタフェース部から出力するステップ、
    をコンピュータに実行させるためのシステム同定プログラムを記録したコンピュータ読み取り可能な記録媒体。
    Figure 0006894402
     Ｈｋ＝［ｕｋ，・・・，ｕｋ−Ｎ＋１］Ｔ

    ここで、
    ｘｋ：状態ベクトルまたは単に状態；未知、これが推定の対象となる。
    ｘ０：初期状態；未知である。
    ｗｋ：システム雑音；未知である。
    ｖｋ：観測雑音；未知である。
    ｙｋ：観測信号；フィルタの入力となり、既知である。
    ｚｋ：出力信号；未知である。
    Ｇｋ：駆動行列；実行時に既知となる。
    Ｈｋ：観測行列；既知である。
    ｘ＾ｋ｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ（未知システムのインパルス応答）の推定値；フィルタ方程式によって与えられる。
    ｘ＾ｋ＋１｜ｋ：観測信号ｙ０〜ｙｋまでを用いた時刻ｋの状態ｘｋ＋１の１ステップ予測値；フィルタ方程式によって与えられる。
    ｘ＾−１｜−１：状態の初期推定値；本来未知であるが、便宜上０が用いられる。
    Ｋｓ，ｋ：フィルタゲイン；行列Σ＾ｋ｜ｋ−１から得られる。
    Ｒｅ，ｋ：補助変数；２×２行列である。
    Ｎ：状態ベクトルの次元（タップ数）；事前に与えられる。
    μ〜 ｋ：ＫＵ ｋのＮ＋１番目の１×２行ベクトル；ＫＵ ｋより得られる。
    Ｗ：重み行列；γｆより決まる。
    ｅｋ：事後前方予測誤差；２次元列ベクトルである。
    ｅ ｋ：事後後方予測誤差；２次元列ベクトルである。
    ε０：任意の正の定数；ε０＞０であり、事前に与えられる。
    γｆ：アッテネーションレベル；設計の際に与えられる。
    μ：ステップサイズ(緩和係数）；μ＞０の実数；事前に与えられる。
    ｕ ｋ ：未知系の入力；観測行列Ｈ ｋ の成分
    ρ：忘却係数
JP2018098503A 2018-05-23 2018-05-23 システム同定装置及び方法及びプログラム及び記憶媒体 Active JP6894402B2 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2018098503A JP6894402B2 (ja) 2018-05-23 2018-05-23 システム同定装置及び方法及びプログラム及び記憶媒体
US16/418,579 US10863272B2 (en) 2018-05-23 2019-05-21 System identification device, system identification method, system identification program, and recording medium recording system identification program
EP19176004.0A EP3573233B1 (en) 2018-05-23 2019-05-22 System identification device, system identification method, system identification program, and recording medium recording system identification program
CN201910428966.8A CN110535452A (zh) 2018-05-23 2019-05-22 系统辨识装置和方法以及记录有系统辨识程序的存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2018098503A JP6894402B2 (ja) 2018-05-23 2018-05-23 システム同定装置及び方法及びプログラム及び記憶媒体

Publications (2)

Publication Number Publication Date
JP2019205049A JP2019205049A (ja) 2019-11-28
JP6894402B2 true JP6894402B2 (ja) 2021-06-30

Family

ID=66826806

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2018098503A Active JP6894402B2 (ja) 2018-05-23 2018-05-23 システム同定装置及び方法及びプログラム及び記憶媒体

Country Status (4)

Country Link
US (1) US10863272B2 (ja)
EP (1) EP3573233B1 (ja)
JP (1) JP6894402B2 (ja)
CN (1) CN110535452A (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113691237B (zh) * 2021-07-27 2024-01-02 浙江工商大学 一种加权融合鲁棒滤波方法

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SE516835C2 (sv) * 1995-02-15 2002-03-12 Ericsson Telefon Ab L M Ekosläckningsförfarande
JP3654470B2 (ja) * 1996-09-13 2005-06-02 日本電信電話株式会社 サブバンド多チャネル音声通信会議用反響消去方法
JP4067269B2 (ja) 2000-10-24 2008-03-26 独立行政法人科学技術振興機構 システム同定方法
CA2399159A1 (en) * 2002-08-16 2004-02-16 Dspfactory Ltd. Convergence improvement for oversampled subband adaptive filters
US7480595B2 (en) * 2003-08-11 2009-01-20 Japan Science And Technology Agency System estimation method and program, recording medium, and system estimation device
EP2012430B1 (en) 2006-04-14 2014-07-30 Incorporated National University Iwate University System identifying method and program, storage medium, and system identifying device
NL2011975C2 (nl) * 2013-12-17 2015-06-18 Berkin Bv Stromingsmeetapparaat van het thermische type.
US9729968B2 (en) * 2014-03-05 2017-08-08 Texas Instruments Incorporated Method and system for acoustic echo cancellation using cascaded kalman filtering
US9881630B2 (en) * 2015-12-30 2018-01-30 Google Llc Acoustic keystroke transient canceler for speech communication terminals using a semi-blind adaptive filter model
US10229698B1 (en) * 2017-06-21 2019-03-12 Amazon Technologies, Inc. Playback reference signal-assisted multi-microphone interference canceler

Also Published As

Publication number Publication date
EP3573233A3 (en) 2019-12-11
EP3573233A2 (en) 2019-11-27
CN110535452A (zh) 2019-12-03
JP2019205049A (ja) 2019-11-28
US10863272B2 (en) 2020-12-08
US20190364361A1 (en) 2019-11-28
EP3573233B1 (en) 2022-11-16

Similar Documents

Publication Publication Date Title
JP5196653B2 (ja) システム同定方法及びプログラム及び記憶媒体、システム同定装置
CN105391879A (zh) 一种无回声残留双端通话鲁棒的声学回声消除方法
JP4067269B2 (ja) システム同定方法
JP6894402B2 (ja) システム同定装置及び方法及びプログラム及び記憶媒体
CN113114865A (zh) 一种组合函数链接型的核自应非线性回声消除方法
Ahmad et al. A 2-D recursive inverse adaptive algorithm
JP4444919B2 (ja) システム推定方法及びプログラム及び記録媒体、システム推定装置
CN108141202B (zh) 包括自适应模块和校正模块的分区块频域自适应滤波器设备
Paleologu et al. Study of the optimal and simplified Kalman filters for echo cancellation
Dunne et al. QR-based TLS and mixed LS-TLS algorithms with applications to adaptive IIR filtering
Paleologu et al. A Kalman filter with individual control factors for echo cancellation
JP4705554B2 (ja) エコーキャンセル装置、その方法、そのプログラム、およびその記録媒体
JP6343585B2 (ja) 未知伝達系推定装置、未知伝達系推定方法、およびプログラム
Yang et al. Fast implementation of a family of memory proportionate affine projection algorithm
JPS6255739B2 (ja)
Patel et al. Nonlinear acoustic echo cancellation using low-complexity low-rank recursive least-squares algorithms
Nishiyama Computational Improvement of the Fast ${\rm H} _ {\infty} $ Filter Based on Information of Input Predictor
Ganjimala et al. A proportionate type block-oriented functional link adaptive filter for sparse nonlinear systems
JP5016501B2 (ja) Ｆｉｒ型適応フィルタ
JP5281548B2 (ja) 音響エコーキャンセル装置、音響エコーキャンセル方法及びそのプログラム
JPH0799468A (ja) 音響反響除去装置
Maqsood Adaptive noise cancellation using single sensor
JP2012253647A (ja) エコーキャンセラ
JP2004128994A (ja) 適応フィルタシステム
JPH0799469A (ja) 音響反響除去装置

Legal Events

Date Code Title Description
A80 Written request to apply exceptions to lack of novelty of invention

Free format text: JAPANESE INTERMEDIATE CODE: A80

Effective date: 20180605

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20200130

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20200403

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20201111

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20201117

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20201204

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20210601

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20210603

R150 Certificate of patent or registration of utility model

Ref document number: 6894402

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150