JP4444919B2 - システム推定方法及びプログラム及び記録媒体、システム推定装置 - Google Patents
システム推定方法及びプログラム及び記録媒体、システム推定装置 Download PDFInfo
- Publication number
- JP4444919B2 JP4444919B2 JP2005513012A JP2005513012A JP4444919B2 JP 4444919 B2 JP4444919 B2 JP 4444919B2 JP 2005513012 A JP2005513012 A JP 2005513012A JP 2005513012 A JP2005513012 A JP 2005513012A JP 4444919 B2 JP4444919 B2 JP 4444919B2
- Authority
- JP
- Japan
- Prior art keywords
- filter
- observation
- processing unit
- matrix
- value
- 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
Links
- 238000000034 method Methods 0.000 title claims description 30
- 239000011159 matrix material Substances 0.000 claims description 100
- 238000012545 processing Methods 0.000 claims description 89
- 238000004422 calculation algorithm Methods 0.000 claims description 53
- 239000013598 vector Substances 0.000 claims description 15
- 238000011156 evaluation Methods 0.000 claims description 14
- 230000014509 gene expression Effects 0.000 claims description 14
- 238000005457 optimization Methods 0.000 claims description 14
- 230000007423 decrease Effects 0.000 claims 6
- 230000004044 response Effects 0.000 description 22
- 238000004364 calculation method Methods 0.000 description 21
- 238000010586 diagram Methods 0.000 description 11
- 238000004891 communication Methods 0.000 description 8
- 230000003044 adaptive effect Effects 0.000 description 7
- 238000001914 filtration Methods 0.000 description 7
- 230000005540 biological transmission Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 230000003321 amplification Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- 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
- 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/0205—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric not using a model or a simulator of the controlled system
- G05B13/024—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric not using a model or a simulator of the controlled system in which a parameter or coefficient is automatically adjusted to optimise the performance
-
- 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
-
- 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
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y10—TECHNICAL SUBJECTS COVERED BY FORMER USPC
- Y10S—TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y10S367/00—Communications, electrical: acoustic wave systems and devices
- Y10S367/901—Noise or unwanted signal reduction in nonseismic receiving system
Landscapes
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Signal Processing (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Physics & Mathematics (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)
- Complex Calculations (AREA)
Description
本発明は、システム推定方法及びプログラム及び記録媒体、システム推定装置に係り、特に、H∞評価基準に基づいて開発されたハイパーH∞フィルタの高速H∞フィルタリングアルゴリズムを用いて、状態推定のロバスト化と忘却係数の最適化を同時に実現するシステム推定方法及びプログラム及び記録媒体、システム推定装置に関する。
【背景技術】
一般に、システム推定とは、入出力データに基づいてシステムの入出力関係の数理モデル(伝達関数、あるいはインパルス応答など)のパラメータを推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、非特許文献1.1993年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
(基本原理)
図8に、システム推定のための構成図の例を示す(未知システムはIIR(Infinite Impulse Response)フィルタで表現してもよい)。
このシステムは、未知システム1、適応フィルタ2を備える。また、適応フィルタ2は、FIRディジタルフィルタ3、適応アルゴリズム4を有する。
以下に、未知システム1を同定する出力誤差方式の一例を説明する。ここで、ukは未知システム1の入力、dkは所望信号であるシステムの出力、d^kはフィルタの出力である。(なお、「^」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。)
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差ek=dk−d^kを最小にするように適応アルゴリズムによってFIRディジタルフィルタ3の係数を調節する。
また、従来、システムのパラメータ(状態)の推定には、誤差共分散行列の更新式(リカッチ方程式)に基づくカルマンフィルタが広く用いられて来た。詳細は、非特許文献2.S.Haykin:Adaptive filter theory,Prentice−Hall(1996)などに示されている。
以下にカルマンフィルタの基本原理について説明する。
次式のように、状態空間モデルで表された線形システム
の状態xkの最小分散推定値x^k|kは、状態の誤差共分散行列Σ^k|k−1を用いて次のように得られる。
ただし、
xk:状態ベクトルまたは単に状態;未知であり、これが推定の対象となる。
yk:観測信号;フィルタの入力となり、既知である。
Hk:観測行列;既知である。
vk:観測雑音;未知である。
ρ:忘却係数;一般に試行錯誤で決定される。
Kk:フィルタゲイン;行列Σ^k|k−1から得られる。
Σ^k|k:x^k|kの誤差の共分散行列に対応;リカッチ方程式により得られる。
Σ^k+1|k:x^k+1|kの誤差の共分散行列に対応;リカッチ方程式により得られる。
Σ^1|0:初期状態の共分散行列に対応;本来未知であるが、便宜上ε0Iが用いられる。
また、本発明者は、既に高速H∞フィルタによるシステム同定アルゴリズムを提案した(特許文献1参照)。これは、システム同定のために新たにH∞評価基準を定め、これに基づくハイパーH∞フィルタの高速アルゴリズムを開発すると共に、この高速H∞フィルタリングアルゴリズムに基づく高速時変システム同定方法である。高速H∞フィルタリングアルゴリズムは、単位時間ステップ当たり計算量O(N)で急激に変化する時変システムの追跡が可能である。上限値の極限で高速カルマンフィルタリングアルゴリズムと完全に一致する。このようなシステム同定により時不変および時変システムの高速実時間同定および推定を実現することができる。
なお、システム推定の分野で通常知られる方法として、例えば、非特許文献2、3参照のこと。
(エコーキャンセラへの適用例)
国際電話など長距離電話回線では,信号増幅などの理由から4線式回線が用いられている。一方、加入者回線は比較的短距離なので、2線式回線が使用されている。
図9に通信系とエコーについての説明図を示す。2線式回線と4線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Bからの信号(音声)は話者Aのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は4線式回線に漏れ、増幅された後、再び受信者(話者A)に戻ると云った現象が起こる。これがエコー(echo)である。エコーは、伝送距離が長くなるにつれて(遅延時間が長くなるにつれて)影響が大きくなり、著しく通話の品質を劣化させる(パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する)。
図10に、エコーキャンセラの原理図を示す。
そこで、図示のようにエコーキャンセラ(echo canceller)を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
エコーパスのインパルス応答の推定は、残留エコーekの平均2乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Aがらの信号(音声)である。一般に、話者2人が同時に話し始めた(ダブルトーク)ときはインパルス応答の推定を中断する。また、ハイブリッドトランスのインパルス応答長は50[ms]程度なので、サンプリング周期を125[μs]とするとエコーパスのインパルス応答の次数は実際は400程度となる。
【非特許文献1】
1993年電子情報通信学会「ディジタル信号処理ハンドブック」
【非特許文献2】
S.Haykin:Adaptive filter theory,Prentice−Hall(1996)
【非特許文献3】
B.Hassibi,A.H.Sayed,and T.Kailath: ″Indefinite−Quadratic Estimation and Control″,SIAM(1996)
【特許文献1】
特開2002−135171号公報
【発明の開示】
しかしながら、式(1)〜(5)のような従来の忘却係数ρを入れたカルマンフィルタでは、忘却係数ρの値は試行錯誤で決定しなければならず非常に手間が掛かった。さらに、決定された忘却ρ係数の値が果たして最適な値であるかどうか判定する手段も無かった。
また、カルマンフィルタで用いる誤差共分散行列は、本来、零でない任意のベクトルとの2次形式が常に正(以下、「正定」という。)であるがコンピュータにより単精度で計算した場合にはその2次形式が負(以下、「負定」という。)となり、数値的に不安定になることが知られている。また、計算量がO(N2)(あるいはO(N3))であるため、状態ベクトルxkの次元Nが大きい場合、1ステップ当たりの演算回数が急激に増大し、実時間処理には適さなかった。
本発明は、以上の点に鑑み、忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム推定方法を提供することを目的とする。
本発明は、上記課題を解決するために、新たに考案したH∞最適化手法を用いて忘却係数が最適決定可能な状態推定アルゴリズムを導出する。さらに、常に正定であるべき誤差共分散行列の代わりに、その因数行列を更新することによって数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。
本発明の第1の解決手段によると、
次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γfに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法及びプログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部は、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーH∞フィルタを実行するステップと、
x^k|k=Fk−1x^k−1|k−1+Ks,k(yk−HkFk−1x^k−1|k−1)
ここで、
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
Fk:システムのダイナミックス
Ks,k:フィルタゲイン
処理部は、ハイパーH∞フィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法、各ステップをコンピュータに実行させるためのシステム推定プログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
また、本発明の第2の解決手段によると、
次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γfに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び/又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
前記処理部は、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力すること、
前記処理部は、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定すること、
前記処理部は、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーH∞フィルタを実行すること、
x^k|k=Fk−1x^k−1|k−1+Ks、k(yk−HkFk−1x^k−1|k−1)
ここで、
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
Fk:システムのダイナミックス
Ks,k:フィルタゲイン
前記処理部は、ハイパーH∞フィルタに関する求められた値を記憶部に記憶すること、
前記処理部は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算すること、
前記処理部は、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること
を備えた前記システム推定装置が提供される。
本発明の推定方法は忘却係数を最適に決定することが可能であり、かつアルゴリズムは単精度でも安定に動作可能であるため、低コストで高い性能が実現できる。一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。
【図面の簡単な説明】
図1は、本実施の形態に関するハードウェアの構成図である。
図2は、H∞フィルタのロバスト化と忘却係数ρの最適化についてのフローチャートである。
図3は、図2中のH∞フィルタ(S105)のアルゴリズムについてのフローチャートである。
図4は、定理2の平方根アレイアルゴリズムの説明図である。
図5は、定理3の数値的に安定な高速アルゴリズムのフローチャートである。
図6は、インパルス応答{hi}i=0 23の値を示す図である。
図7は、定理3の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果である。
図8は、システム推定のための構成図である。
図9は、通信系とエコーについての説明図である。
図10は、エコーキャンセラの原理図である。
【発明を実施するための最良の形態】
以下に、本発明の実施の形態について説明する。
1.記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
xk:状態ベクトルまたは単に状態;未知であり、これが推定の対象となる。
x0:初期状態;未知である。
wk:システム雑音;未知である。
vk:観測雑音;未知である。
yk:観測信号;フィルタの入力となり、既知である。
zk:出力信号;未知である。
Fk:システムのダイナミックス;既知である。
Gk:駆動行列;実行時に既知となる。
Hk:観測行列;既知である。
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値;フィルタ方程式によって与えられる。
x^k+1|k:観測記号y0〜ykで用いた時刻k+1の状態xk+1の推定値;フィルタ方程式によって与えられる。
x^0|0:状態の初期推定値;本来未知であるが、便宜上0が用いられる。
Σ^k|k:x^k|kの誤差の共分散行列に対応;リカッチ方程式によって与えられる。
Σ^k+1|k:x^k+1|kの誤差の共分散行列に対応;リカッチ方程式によって与えられる。
Σ^1|0:初期状態の共分散行列に対応;本来未知であるが、便宜上ε0Iが用いられる。
Ks,k:フィルタゲイン;行列Σ^k|k−1から得られる。
ρ:忘却係数;定理1〜3の場合、γfが決まればρ=1−χ(γf)より自動的に決定される。
ef,i:フィルタ誤差
Re,k:補助変数
なお、記号の上に付される”^”、”v”は、推定値の意味である。また、”〜”、”−”、”U”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、x、w等はベクトル、H、G,K、R、Σ、等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
2.システム推定のハードウェア及びプログラム
本システム推定方法又はシステム推定装置・システムは、その各手順をコンピュータに実行させるためのシステム推定プログラム、システム推定プログラムを記録したコンピュータ読み取り可能な記録媒体、システム推定プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンピュータ、等により提供されることができる。
図1は、本実施の形態に関するハードウェアの構成図である。
このハードウェアは、中央処理装置(CPU)である処理部101、入力部102、出力部103、表示部104及び記憶部105を有する。また、処理部101、入力部102、出力部103、表示部104及び記憶部105は、スター又はバス等の適宜の接続手段で接続されている。記憶部105は、システム推定される「1.記号の説明」で示した既知のデータが必要に応じて記憶される。また、未知・既知のデータや計算されたハイパーH∞フィルタに関するデータ・その他のデータが処理部101により、必要に応じて書込み及び/又は読み出しされる。
3.忘却係数が最適決定可能なハイパーH∞フィルタ
(定理1)
次式のような状態空間モデルを考える。
このような状態空間モデルに対して、次式のようなH∞評価基準を提案する。
このH∞評価基準を満たす状態推定値x^k|k(あるいは出力推定値zv k|k)は、次のレベルγfのハイパーH∞フィルタによって与えられる。
ここで、
である。なお、式(11)はフィルタ方程式、式(12)はフィルタゲイン、式(13)はリカッチ方程式をそれぞれ示す。
また、駆動行列Gkは次のように生成される。
また、上述の高速H∞フィルタの追従能力を向上するためには、上限値γfは次の存在条件を満たすように出来るだけ小さく設定する。
ただし、χ(γf)は、χ(1)=1、χ(∞)=0を満たすγfの単調減衰関数である。
定理1の特徴は状態推定のロバスト化と忘却係数ρの最適化を同時に行っている点にある。
図2に、H∞フィルタのロバスト化と忘却係数ρの最適化についてのフローチャートを示す。ここで、
ブロック「EXC>0」:H∞フィルタの存在条件、
Δγ:正の実数、である。
まず、処理部101は、記憶部105又は入力部102から上限値γfを読み出し又は入力する(S101)。この例では、γf>>1が与えられる。処理部101は、式(15)によって忘却係数ρを決定する(S103)。その後、処理部101は、忘却係数ρに基づき、式(10)〜式(13)のハイパーH∞フィルタを実行する(S105)。処理部101は、式(17)(あるいは、後述の式(18))の右辺(これをEXCとする)を計算し(S107)、その存在条件がすべての時刻で満たされれば(S109)、γfをΔγだけ小さくして同じ処理を繰り返す(S111)。一方、あるγfで存在条件を満たさなくなったとき(S109)、そのγfにΔγを加えたものをγfの最適値γf opとして出力部103に出力及び/又は記憶部105に記憶する(S113)。なお、この例では、Δγを加えているが、それ以外の予め設定された値を加えてもよい。この最適化のプロセスをγ−イタレーションと呼ぶ。なお、処理部101は、H∞フィルタ計算ステップS105及び存在条件の計算ステップS107等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部105に記憶し、また、記憶部105から読み出すようにしてもよい。
ハイパーH∞フィルタが存在条件を満たすとき、式(9)の不等式は常に満たされる。よって、式(9)の分母の外乱エネルギーが限定される場合、式(9)の分子の2乗推定誤差の総和は有界となり、ある時刻以降の推定誤差が0になる。これは、γfをより小さく出来れば、状態xkの変化に推定値x^k|kが速やかに追従できることを意味する。
ここで、定理1のハイパーH∞フィルタのアルゴリズムは通常のH∞フィルタのものとは異なることに注意されたい。また、γf→∞のとき、ρ=1、Gk=0となり、定理1のH∞フィルタのアルゴリズムはカルマンフィルタのアルゴリズムと一致する。
図3に、図2中の(ハイパー)H∞フィルタ(S105)のアルゴリズムについてのフローチャートを示す。
ハイパーH∞フィルタリングアルゴリズムは以下のように要約することができる。
[ステップS201]処理部101は、記憶部105から再帰式の初期条件を読み出し、又は、初期条件を入力部102から入力し、図示のように定める。なお、Lはあらかじめ定められた最大データ数を示す。
[ステップS203]処理部101は、時刻kと最大データ数Lとを比較する。処理部101は、時刻kが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。(不要であれば条件文を取り除くことができる。または、必要に応じて再スタートを行ってもよい。)
[ステップS205]処理部101は、フィルタゲインKs,kを式(12)を用いて計算する。
[ステップS207]処理部101は、式(11)のハイパーH∞フィルタのフィルタ方程式を更新する。
[ステップS209]処理部101は、誤差の共分散行列に対応する項Σ^k|k、Σ^k+1|kを式(13)のリカッチ方程式を用いて計算する。
[ステップS211]時刻kを進ませて(k=k+1)、ステップS203に戻り、データがある限り続ける。
なお、処理部101は、H∞フィルタ計算ステップS205〜S209等の各ステップで求めた適宜の中間値及び最終値、存在条件の値等を必要に応じて適宜記憶部105に記憶し、また、記憶部105から読み出すようにしてもよい。
(スカラー存在条件)
ところで、式(17)の存在条件の判定にはO(N2)の計算量が必要であった。しかし、次の条件を用いれば計算量O(N)で定理1のH∞フィルタの存在性、すなわち式(9)を検証することができる。
系1:スカラー存在条件
次の存在条件を用いれば計算量O(N)でハイパーH∞フィルタの存在性が判定できる。
ここで、
ただし、Ks,iは式(12)で求めたフィルタゲインである。
(証明)
以下に系1の証明を説明する。
2×2の行列Re,kの特性方程式
を解けば、Re,kの固有値λiが次のように得られる。
もし、
であれば、行列Re,kの2つの固有値の1つは正となり、もう1つは負となり、行列RkとRe,kは同じイナーシャをもつ。これより、
を用いれば、式(18)の存在条件が得られる。ここで、HkKs,kの計算量はO(N)である。
4.数値的に安定な状態推定アルゴリズム
上述のハイパーH∞フィルタは、Σ^k|k−1∈RN×Nを更新するため、単位時間ステップ当たりの計算量はO(N2)となる、すなわち、N2に比例する算術演算が必要となる。ここで、Nは状態ベクトルxkの次元である。よって、xkの次元が増加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列Σ^k|k−1は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度(例:32bit)でも動作可能な状態推定アルゴリズムの開発が望まれる。
そこで、次に、
Rk=R1/2 kJ1RT/2 k、
Re,k=R1/2 e,kJ1RT/2 e,k、
Σ^k|k−1=Σ^1/2 k|k−1Σ^T/2 k|k−1
に着目して、数値的に安定化した定理1のH∞フィルタ(平方根アレイアルゴリズム)を定理2に示す。ただし、ここでは簡単のためFk=Iとしたが、Fk≠Iの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパーH∞フィルタを示す。
(定理2)
ただし、
位行列である。また、Kk(:,1)は行列Kkの1列目の列ベクトルを表す。
なお、式(21)、(22)において、J1 −1およびJ1は削除可能である。
図4に、定理2の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図2に示した定理1のフローチャート中のH∞フィルタの計算(S105)で用いることができる。
本推定アルゴリズムは、Σ^k|k−1をリカッチ型の更新式で求める代わりに、その因数行列Σ^1/2 k|k−1∈RN×N(Σ^k|k−1の平方根行列)をJ−ユニタリ変換に基づく更新式で求めている。このとき生じる1−1ブロック行列と2−1ブロック行列からフィルタゲインKs,kを図示のように求めている。このため、Σ^k|k−1=Σ^1/2 k|k−1Σ^T/2 k|k−1>0となり、Σ^k|k−1の正定性は保証され、数値的に安定化できる。なお、定理2のH∞フィルタの単位ステップ当たりの計算量はO(N2)のままである。
なお、図4において、J1 −1は削除可能である。
まず、処理部101は、式(22)の左辺の行列式の各要素に含まれる項を記憶部105から読み出し又は内部メモリ等から得て、J−ユニタリ変換を実行する(S301)。処理部101は、求めた式(22)の右辺の行列式の要素からシステムゲインKk、Ks,kを式(21)に基づき計算する(S303、S305)。処理部191は、式(20)に基づき状態推定値x^k|kを計算する(S307)。
5.状態推定のための数値的に安定な高速アルゴリズム
上述のように、定理2のH∞フィルタの単位ステップ当たりの計算量はO(N2)のままである。そこで、計算量の対策として、H k=H k+1Ψ,H k=[u(k),・・・,u(0),0,・・・,0]のとき、x k=[xT k,0T]Tの1ステップ予測誤差の共分散行列Σ k+1|kが
を満たすことを利用して、Σ k+1|k代わりに次元の低いL k(すなわちL〜 k)を更新する
3が得られる。
(定理3)
ただし、
なお、定理3の証明は、後述する。
上式は、K− k(=ρ−1/2Kk)の代わりにKkについても整理することができる。
さらに、次のJ−ユニタリ行列
を用いれば定理4の高速化した状態推定アルゴリズムが得られる。ただし、Ψはシフト行列を表す。
(定理4)
ただし、
であり、diag{・}は対角行列、Re,k+1(1,1)は行列Re,k+1の1−1成分をそれぞれ表す。また、上式はK− kの代わりにKkに関しても整理できる。
本高速アルゴリズムは、次の因数分解
におけるL〜 k∈R(N+1)×2の更新によってフィルタゲインKs,kを求めているので、単位ステップ当たりの計算量はO(N+1)で済む。ここで、次式に注意されたい。
図5に、定理4の数値的に安定な高速アルゴリズムのフローチャートの一例を示す。この高速アルゴリズムは図2のH∞フィルタの計算ステップ(S105)に組み込まれ、γ−イタレーションによって最適化される。よって、存在条件が満たされる間はγfは除々に減少されるが、満たされなくなった時点で、図示のようにγfは増加される。
H∞フィルタリングアルゴリズムは以下のように要約することができる。
[ステップS401]処理部101は、再帰式の初期条件を図示のように定める。なお、Lは最大データ数を示す。
[ステップS403]処理部101は、時刻kと最大データ数Lとを比較する。処理部101は、時刻kが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。(不要であれば条件文を取り除くことができる。または、再スタートする。)
[ステップS405]処理部101は、フィルタゲインに対応する項Kk+1を式(27)、(31)を用いて再帰的に計算する。
[ステップS406]処理部101は、Re,k+1を式(29)を用いて再帰的に計算する。
[ステップS407]処理部101は、さらにKs,kを式(26)、(31)を用いて計算する。
[ステップS409]処理部101は、ここで、存在条件EXC>0を判定し、存在条件を満たせばステップS411に進む。
[ステップS413]一方、処理部101は、ステップS409で存在条件を満たさなければγfを増加し、ステップS401に戻る。
[ステップS411]処理部101は、式(25)のH∞フィルタのフィルタ方程式を更新する。
[ステップS415]処理部101は、Rr.k+1を式(30)を用いて再帰的に計算する。また、処理部101は、L〜 k+1を式(28)、(31)を用いて再帰的に計算する。
[ステップS419]処理部101は、時刻kを進ませて(k=k+1)、ステップS403に戻り、データがある限り続ける。
なお、処理部101は、H∞フィルタ計算ステップS405〜S415及び存在条件の計算ステップS409等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部105に記憶し、また、記憶部105から読み出すようにしてもよい。
6.エコーキャンセラ
つぎに、エコーキャンセリング問題の数理モデルを作成する。
まず、受信信号{uk}がエコーパスへの入力信号となることを考慮すれば、エコーパスの(時変)インパルス応答{hi[k]}により、エコー{dk}の観測値{yk}は次式で表される。
ここで、uk,ykはそれぞれ時刻tk(=kT;Tは標本化周期)における受信信号とエコーを表し、vkは時刻tkにおける平均値0の回線雑音とし、hi[k],i=0,・・・,N−1は、時変インパルス応答であり、そのタップ数Nは既知とする。このとき、インパルス応答の推定値{h^i[k]}が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
これをエコーから差し引けば(yk−d^k≒0)、エコーをキャンセルすることができる。ただし、k−i<0のときuk−1=0とする。
以上より、問題は直接観測可能な受信信号{uk}とエコー{yk}からエコーパスのインパルス応答{hi[k]}を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにH∞フィルタを適用するには、まず式(32)を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、問題がインパルス応答{hi[k]}を推定することであるから、{hi[k]}を状態変数xkとし、wk程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
ただし、
このような状態空間モデルに対するハイパーおよび高速H∞フィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
7.インパルス応答に対する評価
(動作の確認)
エコーパスのインパルス応答が時間的に不変であり(hi[k]=hi)、かつそのタップ数Nが48である場合について、シミュレーションを用いて、本高速アルゴリズムの動作を確認する。
なお、図6は、ここでのインパルス応答{hi}の値を示す図である。
ここで、インパルス応答{hi}i=0 23は、図示の値を採用し、その他{hi}i=24 47は0とする。また、υkは平均値0、分散σv 2=1.0×10−6の定常なガウス白色雑音とし、標本化周期Tを便宜上1.0とする。
また、受信信号{uk}は次のように2次のARモデルで近似する。
ただし、α1=0.7,α2=0.1とし、wk′は平均値0、分散σw′ 2=0.04の定常なガウス白色雑音とする。
(インパルス応答の推定結果)
図7に、定理4の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図7(b)の縦軸は、
√{Σi=0 47(hi−x^k(i+1))2}
を表す。
これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし、ρ=1−χ(γf)、χ(γf)=γf −2、x^0|0=0、Σ^1|0=20Iとし、計算は倍精度で行った。また、存在条件を確認しつつ、γf=5.5と設定とした。
8.定理の証明
8−1.定理2の証明
次の関係式
が成り立つとき、両辺の2×2ブロック行列の各項を比較すれば次式が得られる。
これは定理1のFk=Iのときの式(13)のリカッチ方程式と一致する。ただし、
一方、AJAT=BJBTが成り立つとき、BはJ−ユニタリ行列Θ(k)を用いてB=AΘ(k)と表すことができる。よって、式(40)より定理1のリカッチ方程式は次式と等価である。
なお、式(40)、(45)において、J1 −1は削除可能である。
8−2.定理3の証明
次のようにブロック三角化するJ−ユニタリ行列Θ(k)が存在すると仮定する。
のように決定することができる。ただし、Sは対角成分が1または−1をとる対角行列とする。
(1,1)−ブロック行列
(2,1)−ブロック行列
(2,2)−ブロック行列
これより、
8−3.定理4の証明
観測行列Hkがシフト特性をもち、かつ
のとき、定理2と同様な方法によって次の関係式が得られる。
ただし、
となるようにRr,k+1を決定する。次に、式(46)の3行目にRr,k+1の更新式を新たに追加すれば、最終的に次式が得られる。
この両辺の3×2ブロック行列の各項の対応から次のゲイン行列K− kの更新式が得られる。
【産業上の利用可能性】
一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができる。
Claims (12)
- 次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として、システム雑音wk及び観測雑音vkを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γfに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kkを用いて、次式(20)〜(22)により求める、又は、次式(20)と、式(21)及び(22)においてJ1 −1およびJ1を削除した式により求めるハイパーH∞フィルタを実行するステップと、
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
yk:観測信号
Fk:システムのダイナミックス
Ks,k:フィルタゲイン
Hk:観測行列
Σ^k|k:x^k|kの誤差の共分散行列に対応
Θ(k):J−ユニタリ行列
Re,k:補助変数
処理部が、ハイパーH∞フィルタによって求められた状態xkの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Hi、又は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。 - 前記ハイパーH∞フィルタを実行するステップは、
処理部が、Σ^k+1|k 1/2を前記式(22)を用いて計算するステップと、
処理部が、Σ^k|kの初期条件とCkの初期条件のもとで、フィルタゲインKs,kを前記式(21)を用いて計算するステップと、
処理部が、前記式(20)のH∞フィルタのフィルタ方程式を更新するステップと、
処理部が、前記式(22)を用いて計算するステップと、前記式(21)を用いて計算するステップと、前記更新するステップとを、時刻kを進ませて繰り返し実行するステップと
を含む請求項1に記載のシステム推定方法。 - 次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として、システム雑音wk及び観測雑音vkを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γfに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kkを用いて、次式により求めるハイパーH∞フィルタを実行するステップと
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
yk:観測信号
Ks,k:フィルタゲイン
Hk:観測行列
Θ(k):J−ユニタリ行列
Re,k:補助変数
処理部が、ハイパーH∞フィルタによって求められた状態xkの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Hi、又は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。 - 前記ハイパーH∞フィルタを実行するステップは、
処理部が、Re,k+1、Rr,k+1及びL〜 k+1の初期条件のもとで、K― kを前記式(63)を用いて計算するステップと、
処理部が、フィルタゲインKs,kを前記式(62)を用いて計算するステップと、
処理部が、前記式(61)のH∞フィルタのフィルタ方程式を更新するステップと、
処理部は、前記式(63)を用いて計算するステップと、前記式(62)を用いて計算するステップと、前記更新するステップを、時刻kを進ませて繰り返し実行するステップと
を含む請求項4に記載のシステム推定方法。 - 次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として、システム雑音wk及び観測雑音vkを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γfに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列K― kを用いて、次式により求めるハイパーH∞フィルタを実行するステップと、
yk:観測信号
Fk:システムのダイナミックス
Hk:観測行列
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
Ks,k:フィルタゲイン;ゲイン行列K― kから得られる。
Re,k、L〜 k:補助変数
処理部が、ハイパーH∞フィルタによって求められた状態xkの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Hi、又は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。 - さらに、次式により時刻kの状態推定値x^k|kから出力信号の推定値zv k|kを求めるようにした請求項1又は4又は6に記載のシステム推定方法。
zv k|k=Hkx^k|k - 次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として、システム雑音wk及び観測雑音vkを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γfに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムであって、
処理部が、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列K― kを用いて、次式により求めるハイパーH∞フィルタを実行するステップと、
yk:観測信号
Fk:システムのダイナミックス
Hk:観測行列
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
Ks,k:フィルタゲイン;ゲイン行列K― kから得られる。
Re,k、L〜 k:補助変数
処理部が、ハイパーH∞フィルタによって求められた状態xkの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Hi、又は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラム。 - 次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として、システム雑音wk及び観測雑音vkを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γfに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列K― kを用いて、次式により求めるハイパーH∞フィルタを実行するステップと、
yk:観測信号
Fk:システムのダイナミックス
Hk:観測行列
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
Ks,k:フィルタゲイン;ゲイン行列K― kから得られる。
Re,k、L〜 k:補助変数
処理部が、ハイパーH∞フィルタによって求められた状態xkの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列Hi、又は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体。 - 次式で表される状態空間モデルに対して、
xk+1=Fkxk+Gkwk
yk=Hkxk+vk
zk=Hkxk
ここで、
xk:状態ベクトルまたは単に状態
wk:システム雑音
vk:観測雑音
yk:観測信号
zk:出力信号
Fk:システムのダイナミックス
Gk:駆動行列
評価基準として、システム雑音wk及び観測雑音vkを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γfに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び/又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
処理部が、上限値γf、フィルタの入力である観測信号yk、観測行列Hkを含む値を記憶部又は入力部から入力する手段と、
処理部が、前記上限値γfに従い、状態空間モデルに関連する忘却係数ρを決定する手段と、
処理部が、記憶部から初期値又はある時刻の観測行列Hkを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列K― kを用いて、次式により求めるハイパーH∞フィルタを実行する手段と
yk:観測信号
Fk:システムのダイナミックス
Hk:観測行列
x^k|k:観測信号y0〜ykまでを用いた時刻kの状態xkの推定値
Ks,k:フィルタゲイン;ゲイン行列K― kから得られる。
Re,k、L〜 k:補助変数
処理部が、ハイパーH∞フィルタによって求められた状態xkの推定値を記憶部に記憶する手段と、
処理部が、求められた観測行列Hi、又は、観測行列HiとフィルタゲインKs,iにより、前記上限値γf及び前記忘却係数ρに基づく存在条件を計算する手段と、
処理部が、上限値γfを小さくしていき前記ハイパーH∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶する手段と、
を備えた前記システム推定装置。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2003291614 | 2003-08-11 | ||
JP2003291614 | 2003-08-11 | ||
PCT/JP2004/011568 WO2005015737A1 (ja) | 2003-08-11 | 2004-08-05 | システム推定方法及びプログラム及び記録媒体、システム推定装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
JPWO2005015737A1 JPWO2005015737A1 (ja) | 2006-10-12 |
JP4444919B2 true JP4444919B2 (ja) | 2010-03-31 |
Family
ID=34131662
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2005513012A Active JP4444919B2 (ja) | 2003-08-11 | 2004-08-05 | システム推定方法及びプログラム及び記録媒体、システム推定装置 |
Country Status (5)
Country | Link |
---|---|
US (1) | US7480595B2 (ja) |
EP (2) | EP1657818B1 (ja) |
JP (1) | JP4444919B2 (ja) |
CN (1) | CN1836372B (ja) |
WO (1) | WO2005015737A1 (ja) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP3573233A2 (en) | 2018-05-23 | 2019-11-27 | National University Corporation Iwate University | System identification device, system identification method, system identification program, and recording medium recording system identification program |
Families Citing this family (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP5196653B2 (ja) | 2006-04-14 | 2013-05-15 | 国立大学法人岩手大学 | システム同定方法及びプログラム及び記憶媒体、システム同定装置 |
JP5115952B2 (ja) * | 2007-03-19 | 2013-01-09 | 学校法人東京理科大学 | 雑音抑圧装置および雑音抑圧方法 |
US8924337B2 (en) * | 2011-05-09 | 2014-12-30 | Nokia Corporation | Recursive Bayesian controllers for non-linear acoustic echo cancellation and suppression systems |
IL218047A (en) | 2012-02-12 | 2017-09-28 | Elta Systems Ltd | Devices and methods are added, for spatial suppression of interference in wireless networks |
KR101716250B1 (ko) * | 2014-03-11 | 2017-03-15 | 메이덴샤 코포레이션 | 구동 트레인의 시험 시스템 |
US10014906B2 (en) | 2015-09-25 | 2018-07-03 | Microsemi Semiconductor (U.S.) Inc. | Acoustic echo path change detection apparatus and method |
US10122863B2 (en) | 2016-09-13 | 2018-11-06 | Microsemi Semiconductor (U.S.) Inc. | Full duplex voice communication system and method |
US10761522B2 (en) * | 2016-09-16 | 2020-09-01 | Honeywell Limited | Closed-loop model parameter identification techniques for industrial model-based process controllers |
CN110069870A (zh) * | 2019-04-28 | 2019-07-30 | 河海大学 | 一种基于自适应无迹h∞滤波的发电机动态状态估计方法 |
CN116679026A (zh) * | 2023-06-27 | 2023-09-01 | 江南大学 | 自适应无偏有限脉冲响应滤波的污水溶解氧浓度估计方法 |
Family Cites Families (13)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPS61200713A (ja) | 1985-03-04 | 1986-09-05 | Oki Electric Ind Co Ltd | デイジタルフイルタ |
US5394322A (en) | 1990-07-16 | 1995-02-28 | The Foxboro Company | Self-tuning controller that extracts process model characteristics |
JP2872547B2 (ja) | 1993-10-13 | 1999-03-17 | シャープ株式会社 | 格子型フィルタを用いた能動制御方法および装置 |
JPH07158625A (ja) * | 1993-12-03 | 1995-06-20 | 英樹 ▲濱▼野 | 吊り下げ部材の取付け用金具 |
JPH07185625A (ja) | 1993-12-27 | 1995-07-25 | Nippon Steel Corp | 帯状鋼板の最低板厚保障のための制御方法 |
SE516835C2 (sv) | 1995-02-15 | 2002-03-12 | Ericsson Telefon Ab L M | Ekosläckningsförfarande |
US5734715A (en) * | 1995-09-13 | 1998-03-31 | France Telecom | Process and device for adaptive identification and adaptive echo canceller relating thereto |
US5987444A (en) | 1997-09-23 | 1999-11-16 | Lo; James Ting-Ho | Robust neutral systems |
US6175602B1 (en) * | 1998-05-27 | 2001-01-16 | Telefonaktiebolaget Lm Ericsson (Publ) | Signal noise reduction by spectral subtraction using linear convolution and casual filtering |
US6487257B1 (en) * | 1999-04-12 | 2002-11-26 | Telefonaktiebolaget L M Ericsson | Signal noise reduction by time-domain spectral subtraction using fixed filters |
US6711598B1 (en) | 1999-11-11 | 2004-03-23 | Tokyo Electron Limited | Method and system for design and implementation of fixed-point filters for control and signal processing |
US6801881B1 (en) | 2000-03-16 | 2004-10-05 | Tokyo Electron Limited | Method for utilizing waveform relaxation in computer-based simulation models |
JP4067269B2 (ja) | 2000-10-24 | 2008-03-26 | 独立行政法人科学技術振興機構 | システム同定方法 |
-
2004
- 2004-05-08 US US10/567,514 patent/US7480595B2/en active Active
- 2004-08-05 JP JP2005513012A patent/JP4444919B2/ja active Active
- 2004-08-05 EP EP04771550.3A patent/EP1657818B1/en active Active
- 2004-08-05 CN CN200480022991.8A patent/CN1836372B/zh active Active
- 2004-08-05 EP EP12007720.1A patent/EP2560281B1/en active Active
- 2004-08-05 WO PCT/JP2004/011568 patent/WO2005015737A1/ja active Application Filing
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP3573233A2 (en) | 2018-05-23 | 2019-11-27 | National University Corporation Iwate University | System identification device, system identification method, system identification program, and recording medium recording system identification program |
US10863272B2 (en) | 2018-05-23 | 2020-12-08 | National University Corporation, Iwate University | System identification device, system identification method, system identification program, and recording medium recording system identification program |
Also Published As
Publication number | Publication date |
---|---|
CN1836372B (zh) | 2010-04-28 |
CN1836372A (zh) | 2006-09-20 |
EP2560281A3 (en) | 2013-03-13 |
EP1657818A1 (en) | 2006-05-17 |
US7480595B2 (en) | 2009-01-20 |
EP2560281B1 (en) | 2015-04-15 |
JPWO2005015737A1 (ja) | 2006-10-12 |
US20070185693A1 (en) | 2007-08-09 |
EP1657818A4 (en) | 2012-09-05 |
WO2005015737A1 (ja) | 2005-02-17 |
EP2560281A2 (en) | 2013-02-20 |
EP1657818B1 (en) | 2015-04-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP4444919B2 (ja) | システム推定方法及びプログラム及び記録媒体、システム推定装置 | |
US8244787B2 (en) | Optimum nonlinear correntropy filter | |
US8380466B2 (en) | System identification method and program, storage medium, and system identification device | |
Nascimento et al. | Adaptive filters | |
JP4067269B2 (ja) | システム同定方法 | |
Scarpiniti et al. | Spline adaptive filters: Theory and applications | |
Nascimento et al. | Adaptive filters | |
US11610598B2 (en) | Voice enhancement in presence of noise | |
Kabzinski et al. | A unified perspective on time-domain and frequency-domain Kalman filters for acoustic system identification | |
Gannot et al. | On the application of the unscented Kalman filter to speech processing | |
JP6894402B2 (ja) | システム同定装置及び方法及びプログラム及び記憶媒体 | |
Caballero-Aguila et al. | A new estimation algorithm from measurements with multiple-step random delays and packet dropouts | |
JP6343585B2 (ja) | 未知伝達系推定装置、未知伝達系推定方法、およびプログラム | |
JPH1168518A (ja) | 係数更新回路 | |
Resende et al. | Analysis of the least mean squares algorithm with reusing coefficient vector | |
JP2002533970A (ja) | 安定適応フィルタおよびその方法 | |
JP2001077730A (ja) | 適応フィルタの係数推定装置 | |
JPH08250981A (ja) | フィルタ係数の推定装置 | |
Timoney et al. | Robustness Analysis of the Adaptive Periodic Noise Canceller Applied to Resonance Cancellation | |
JP4344305B2 (ja) | 未知系推定方法およびこれを実施する装置 | |
JP2006128947A (ja) | 未知系推定方法およびこの方法を実施する装置 | |
JPH08305371A (ja) | 適応的伝達関数推定法及びそれを用いた推定装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20061003 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20091020 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20091201 |
|
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: 20100112 |
|
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20100114 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 4444919 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130122 Year of fee payment: 3 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140122 Year of fee payment: 4 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
S533 | Written request for registration of change of name |
Free format text: JAPANESE INTERMEDIATE CODE: R313533 |
|
R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |