TOP > 国内特許検索 > システム推定方法及びプログラム及び記録媒体、システム推定装置 > 明細書

明細書 :システム推定方法及びプログラム及び記録媒体、システム推定装置

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4444919号 (P4444919)
登録日 平成22年1月22日(2010.1.22)
発行日 平成22年3月31日(2010.3.31)
発明の名称または考案の名称 システム推定方法及びプログラム及び記録媒体、システム推定装置
国際特許分類 H03H  21/00        (2006.01)
H04B   3/23        (2006.01)
H04R   3/02        (2006.01)
FI H03H 21/00
H04B 3/23
H04R 3/02
請求項の数または発明の数 12
全頁数 29
出願番号 特願2005-513012 (P2005-513012)
出願日 平成16年8月5日(2004.8.5)
国際出願番号 PCT/JP2004/011568
国際公開番号 WO2005/015737
国際公開日 平成17年2月17日(2005.2.17)
優先権出願番号 2003291614
優先日 平成15年8月11日(2003.8.11)
優先権主張国 日本国(JP)
審査請求日 平成18年10月3日(2006.10.3)
特許権者または実用新案権者 【識別番号】503360115
【氏名又は名称】独立行政法人科学技術振興機構
発明者または考案者 【氏名】西山 清
個別代理人の代理人 【識別番号】100107010、【弁理士】、【氏名又は名称】橋爪 健
審査官 【審査官】東 昌秋
参考文献・文献 特開2002-135171(JP,A)
特開平7-158625(JP,A)
K.NISHIYAMA,Robust estimation of a single complex sinusoid in white noise-H∞ filtering approach,IEEE Transactions on Signal Processing,1999年10月,vol.47, no.10,pp.2853-2856
K.NISHIYAMA,H∞-learning of layered neural networks,IEEE Transactions on Neural Networks,2001年11月,vol.12, no.6,pp.1265-1277
調査した分野 H03H 21/00
H05B 13/02
H04B 3/23
H04B 7/005
H04R 3/00
H04S 7/00
G10K 11/178
特許請求の範囲 【請求項1】
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として、システム雑音w及び観測雑音vを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kを用いて、次式(20)~(22)により求める、又は、次式(20)と、式(21)及び(22)においてJ-1およびJを削除した式により求めるハイパーHフィルタを実行するステップと、
【数1】
JP0004444919B2_000050t.gif
ここで、
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
:観測信号
:システムのダイナミックス
s,k:フィルタゲイン
:観測行列
Σ^k|k:x^k|kの誤差の共分散行列に対応
Θ(k):J-ユニタリ行列
e,k:補助変数
処理部が、ハイパーHフィルタによって求められた状態xの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列H、又は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項2】
前記ハイパーHフィルタを実行するステップは、
処理部が、Σ^k+1|k1/2を前記式(22)を用いて計算するステップと、
処理部が、Σ^k|kの初期条件とCの初期条件のもとで、フィルタゲインKs,kを前記式(21)を用いて計算するステップと、
処理部が、前記式(20)のHフィルタのフィルタ方程式を更新するステップと、
処理部が、前記式(22)を用いて計算するステップと、前記式(21)を用いて計算するステップと、前記更新するステップとを、時刻kを進ませて繰り返し実行するステップと
を含む請求項に記載のシステム推定方法。
【請求項3】
処理部は、前記存在条件を次式に従い計算する請求項に記載のシステム推定方法。
【数2】
JP0004444919B2_000051t.gif

【請求項4】
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として、システム雑音w及び観測雑音vを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kを用いて、次式により求めるハイパーHフィルタを実行するステップと
【数3】
JP0004444919B2_000052t.gif
ここで、
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
:観測信号
s,k:フィルタゲイン
:観測行列
Θ(k):J-ユニタリ行列
e,k:補助変数
処理部が、ハイパーHフィルタによって求められた状態xの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列H、又は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項5】
前記ハイパーHフィルタを実行するステップは、
処理部が、Re,k+1、Rr,k+1及びLk+1の初期条件のもとで、Kを前記式(63)を用いて計算するステップと、
処理部が、フィルタゲインKs,kを前記式(62)を用いて計算するステップと、
処理部が、前記式(61)のHフィルタのフィルタ方程式を更新するステップと、
処理部は、前記式(63)を用いて計算するステップと、前記式(62)を用いて計算するステップと、前記更新するステップを、時刻kを進ませて繰り返し実行するステップと
を含む請求項に記載のシステム推定方法。
【請求項6】
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として、システム雑音w及び観測雑音vを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法であって、
処理部が、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kを用いて、次式により求めるハイパーHフィルタを実行するステップと、
【数4】
JP0004444919B2_000053t.gif
ここで、
:観測信号
:システムのダイナミックス
:観測行列
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
s,k:フィルタゲイン;ゲイン行列Kから得られる。
e,k、L:補助変数
処理部が、ハイパーHフィルタによって求められた状態xの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列H、又は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法。
【請求項7】
処理部は、前記存在条件を次式に従い計算する請求項又は又はに記載のシステム推定方法。
【数5】
JP0004444919B2_000054t.gif
ただし、前記忘却係数ρ及び前記上限値γは、次式の関係である。
0<ρ=1-χ(γ)≦1 (ただし、χ(γ)は、χ(1)=1、χ(∞)=0を満たすγの単調減衰関数)
【請求項8】
さらに、次式により時刻kの状態推定値x^k|kから出力信号の推定値zk|kを求めるようにした請求項又は又はに記載のシステム推定方法。
k|k=Hx^k|k
【請求項9】
前記Hフィルタ方程式を適用し、状態推定値x^k|k=[h^[k],・・・,h^[k]]を求め、
擬似エコーを次式のように推定し、
求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現する請求項又は又はに記載のシステム推定方法。
【数6】
JP0004444919B2_000055t.gif

【請求項10】
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として、システム雑音w及び観測雑音vを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムであって、
処理部が、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kを用いて、次式により求めるハイパーHフィルタを実行するステップと、
【数7】
JP0004444919B2_000056t.gif
ここで、
:観測信号
:システムのダイナミックス
:観測行列
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
s,k:フィルタゲイン;ゲイン行列Kから得られる。
e,k、L:補助変数
処理部が、ハイパーHフィルタによって求められた状態xの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列H、又は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラム。
【請求項11】
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として、システム雑音w及び観測雑音vを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時にコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部が、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力するステップと、
処理部が、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部が、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kを用いて、次式により求めるハイパーHフィルタを実行するステップと、
【数8】
JP0004444919B2_000057t.gif
ここで、
:観測信号
:システムのダイナミックス
:観測行列
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
s,k:フィルタゲイン;ゲイン行列Kから得られる。
e,k、L:補助変数
処理部が、ハイパーHフィルタによって求められた状態xの推定値を記憶部に記憶するステップと、
処理部が、求められた観測行列H、又は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部が、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
をコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項12】
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として、システム雑音w及び観測雑音vを含む外乱に対するフィルタ誤差の割合を示し且つ忘却係数ρで重み付けされたエネルギーゲインの最大値を、予め与えられた上限値γに対応する項より小さく抑えるように定めた推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び/又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
処理部が、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力する手段と、
処理部が、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定する手段と、
処理部が、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、フィルタゲインKs,kを、前記忘却係数ρとゲイン行列Kを用いて、次式により求めるハイパーHフィルタを実行する手段と
【数9】
JP0004444919B2_000058t.gif
ここで、
:観測信号
:システムのダイナミックス
:観測行列
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
s,k:フィルタゲイン;ゲイン行列Kから得られる。
e,k、L:補助変数
処理部が、ハイパーHフィルタによって求められた状態xの推定値を記憶部に記憶する手段と、
処理部が、求められた観測行列H、又は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算する手段と、
処理部が、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶する手段と、
を備えた前記システム推定装置。
発明の詳細な説明 【技術分野】
本発明は、システム推定方法及びプログラム及び記録媒体、システム推定装置に係り、特に、H評価基準に基づいて開発されたハイパーHフィルタの高速Hフィルタリングアルゴリズムを用いて、状態推定のロバスト化と忘却係数の最適化を同時に実現するシステム推定方法及びプログラム及び記録媒体、システム推定装置に関する。
【背景技術】
一般に、システム推定とは、入出力データに基づいてシステムの入出力関係の数理モデル(伝達関数、あるいはインパルス応答など)のパラメータを推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、非特許文献1.1993年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
(基本原理)
図8に、システム推定のための構成図の例を示す(未知システムはIIR(Infinite Impulse Response)フィルタで表現してもよい)。
このシステムは、未知システム1、適応フィルタ2を備える。また、適応フィルタ2は、FIRディジタルフィルタ3、適応アルゴリズム4を有する。
以下に、未知システム1を同定する出力誤差方式の一例を説明する。ここで、uは未知システム1の入力、dは所望信号であるシステムの出力、d^はフィルタの出力である。(なお、「^」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。)
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差e=d-d^を最小にするように適応アルゴリズムによってFIRディジタルフィルタ3の係数を調節する。
また、従来、システムのパラメータ(状態)の推定には、誤差共分散行列の更新式(リカッチ方程式)に基づくカルマンフィルタが広く用いられて来た。詳細は、非特許文献2.S.Haykin:Adaptive filter theory,Prentice-Hall(1996)などに示されている。
以下にカルマンフィルタの基本原理について説明する。
次式のように、状態空間モデルで表された線形システム
JP0004444919B2_000002t.gifの状態xの最小分散推定値x^k|kは、状態の誤差共分散行列Σ^k|k-1を用いて次のように得られる。
JP0004444919B2_000003t.gif ただし、
JP0004444919B2_000004t.gif:状態ベクトルまたは単に状態;未知であり、これが推定の対象となる。
:観測信号;フィルタの入力となり、既知である。
:観測行列;既知である。
:観測雑音;未知である。
ρ:忘却係数;一般に試行錯誤で決定される。
:フィルタゲイン;行列Σ^k|k-1から得られる。
Σ^k|k:x^k|kの誤差の共分散行列に対応;リカッチ方程式により得られる。
Σ^k+1|k:x^k+1|kの誤差の共分散行列に対応;リカッチ方程式により得られる。
Σ^1|0:初期状態の共分散行列に対応;本来未知であるが、便宜上εIが用いられる。
また、本発明者は、既に高速Hフィルタによるシステム同定アルゴリズムを提案した(特許文献1参照)。これは、システム同定のために新たにH評価基準を定め、これに基づくハイパーHフィルタの高速アルゴリズムを開発すると共に、この高速Hフィルタリングアルゴリズムに基づく高速時変システム同定方法である。高速Hフィルタリングアルゴリズムは、単位時間ステップ当たり計算量O(N)で急激に変化する時変システムの追跡が可能である。上限値の極限で高速カルマンフィルタリングアルゴリズムと完全に一致する。このようなシステム同定により時不変および時変システムの高速実時間同定および推定を実現することができる。
なお、システム推定の分野で通常知られる方法として、例えば、非特許文献2、3参照のこと。
(エコーキャンセラへの適用例)
国際電話など長距離電話回線では,信号増幅などの理由から4線式回線が用いられている。一方、加入者回線は比較的短距離なので、2線式回線が使用されている。
図9に通信系とエコーについての説明図を示す。2線式回線と4線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Bからの信号(音声)は話者Aのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は4線式回線に漏れ、増幅された後、再び受信者(話者A)に戻ると云った現象が起こる。これがエコー(echo)である。エコーは、伝送距離が長くなるにつれて(遅延時間が長くなるにつれて)影響が大きくなり、著しく通話の品質を劣化させる(パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する)。
図10に、エコーキャンセラの原理図を示す。
そこで、図示のようにエコーキャンセラ(echo canceller)を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
エコーパスのインパルス応答の推定は、残留エコーeの平均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(N)(あるいはO(N))であるため、状態ベクトルxの次元Nが大きい場合、1ステップ当たりの演算回数が急激に増大し、実時間処理には適さなかった。
本発明は、以上の点に鑑み、忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム推定方法を提供することを目的とする。
本発明は、上記課題を解決するために、新たに考案したH最適化手法を用いて忘却係数が最適決定可能な状態推定アルゴリズムを導出する。さらに、常に正定であるべき誤差共分散行列の代わりに、その因数行列を更新することによって数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。
本発明の第1の解決手段によると、
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定方法及びプログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
処理部は、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力するステップと、
処理部は、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定するステップと、
処理部は、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーHフィルタを実行するステップと、
x^k|k=Fk-1x^k-1|k-1+Ks,k(y-Hk-1x^k-1|k-1
ここで、
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
:システムのダイナミックス
s,k:フィルタゲイン
処理部は、ハイパーHフィルタに関する求められた値を記憶部に記憶するステップと、
処理部は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算するステップと、
処理部は、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、
を含む前記システム推定方法、各ステップをコンピュータに実行させるためのシステム推定プログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
また、本発明の第2の解決手段によると、
次式で表される状態空間モデルに対して、
k+1=F+G
=H+v
=H
ここで、
:状態ベクトルまたは単に状態
:システム雑音
:観測雑音
:観測信号
:出力信号
:システムのダイナミックス
:駆動行列
評価基準として忘却係数ρで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値γに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数ρの最適化を同時に行うためのシステム推定装置であって、
推定アルゴリズムを実行する処理部と、
前記処理部により読み取り及び/又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、
を備え、
前記処理部は、上限値γ、フィルタの入力である観測信号y、観測行列Hを含む値を記憶部又は入力部から入力すること、
前記処理部は、前記上限値γに従い、状態空間モデルに関連する忘却係数ρを決定すること、
前記処理部は、記憶部から初期値又はある時刻の観測行列Hを含む値を読み取り、前記忘却係数ρを用いて次式で表されるハイパーHフィルタを実行すること、
x^k|k=Fk-1x^k-1|k-1+Ks、k(y-Hk-1x^k-1|k-1
ここで、
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値
:システムのダイナミックス
s,k:フィルタゲイン
前記処理部は、ハイパーHフィルタに関する求められた値を記憶部に記憶すること、
前記処理部は、観測行列HとフィルタゲインKs,iにより、前記上限値γ及び前記忘却係数ρに基づく存在条件を計算すること、
前記処理部は、上限値γを小さくしていき前記ハイパーHフィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること
を備えた前記システム推定装置が提供される。
本発明の推定方法は忘却係数を最適に決定することが可能であり、かつアルゴリズムは単精度でも安定に動作可能であるため、低コストで高い性能が実現できる。一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。
【図面の簡単な説明】
図1は、本実施の形態に関するハードウェアの構成図である。
図2は、Hフィルタのロバスト化と忘却係数ρの最適化についてのフローチャートである。
図3は、図2中のHフィルタ(S105)のアルゴリズムについてのフローチャートである。
図4は、定理2の平方根アレイアルゴリズムの説明図である。
図5は、定理3の数値的に安定な高速アルゴリズムのフローチャートである。
図6は、インパルス応答{h23の値を示す図である。
図7は、定理3の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果である。
図8は、システム推定のための構成図である。
図9は、通信系とエコーについての説明図である。
図10は、エコーキャンセラの原理図である。
【発明を実施するための最良の形態】
以下に、本発明の実施の形態について説明する。
1.記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
:状態ベクトルまたは単に状態;未知であり、これが推定の対象となる。
:初期状態;未知である。
:システム雑音;未知である。
:観測雑音;未知である。
:観測信号;フィルタの入力となり、既知である。
:出力信号;未知である。
:システムのダイナミックス;既知である。
:駆動行列;実行時に既知となる。
:観測行列;既知である。
x^k|k:観測信号y~yまでを用いた時刻kの状態xの推定値;フィルタ方程式によって与えられる。
x^k+1|k:観測記号y~yで用いた時刻k+1の状態xk+1の推定値;フィルタ方程式によって与えられる。
x^0|0:状態の初期推定値;本来未知であるが、便宜上0が用いられる。
Σ^k|k:x^k|kの誤差の共分散行列に対応;リカッチ方程式によって与えられる。
Σ^k+1|k:x^k+1|kの誤差の共分散行列に対応;リカッチ方程式によって与えられる。
Σ^1|0:初期状態の共分散行列に対応;本来未知であるが、便宜上εIが用いられる。
s,k:フィルタゲイン;行列Σ^k|k-1から得られる。
ρ:忘却係数;定理1~3の場合、γが決まればρ=1-χ(γ)より自動的に決定される。
f,i:フィルタ誤差
e,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)
次式のような状態空間モデルを考える。
JP0004444919B2_000005t.gifこのような状態空間モデルに対して、次式のようなH評価基準を提案する。
JP0004444919B2_000006t.gif このH評価基準を満たす状態推定値x^k|k(あるいは出力推定値zk|k)は、次のレベルγのハイパーHフィルタによって与えられる。
JP0004444919B2_000007t.gifここで、
JP0004444919B2_000008t.gifである。なお、式(11)はフィルタ方程式、式(12)はフィルタゲイン、式(13)はリカッチ方程式をそれぞれ示す。
また、駆動行列Gは次のように生成される。
JP0004444919B2_000009t.gif また、上述の高速Hフィルタの追従能力を向上するためには、上限値γは次の存在条件を満たすように出来るだけ小さく設定する。
JP0004444919B2_000010t.gifただし、χ(γ)は、χ(1)=1、χ(∞)=0を満たすγの単調減衰関数である。
定理1の特徴は状態推定のロバスト化と忘却係数ρの最適化を同時に行っている点にある。
図2に、Hフィルタのロバスト化と忘却係数ρの最適化についてのフローチャートを示す。ここで、
ブロック「EXC>0」:Hフィルタの存在条件、
Δγ:正の実数、である。
まず、処理部101は、記憶部105又は入力部102から上限値γを読み出し又は入力する(S101)。この例では、γ>>1が与えられる。処理部101は、式(15)によって忘却係数ρを決定する(S103)。その後、処理部101は、忘却係数ρに基づき、式(10)~式(13)のハイパーHフィルタを実行する(S105)。処理部101は、式(17)(あるいは、後述の式(18))の右辺(これをEXCとする)を計算し(S107)、その存在条件がすべての時刻で満たされれば(S109)、γをΔγだけ小さくして同じ処理を繰り返す(S111)。一方、あるγで存在条件を満たさなくなったとき(S109)、そのγにΔγを加えたものをγの最適値γopとして出力部103に出力及び/又は記憶部105に記憶する(S113)。なお、この例では、Δγを加えているが、それ以外の予め設定された値を加えてもよい。この最適化のプロセスをγ-イタレーションと呼ぶ。なお、処理部101は、Hフィルタ計算ステップS105及び存在条件の計算ステップS107等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部105に記憶し、また、記憶部105から読み出すようにしてもよい。
ハイパーHフィルタが存在条件を満たすとき、式(9)の不等式は常に満たされる。よって、式(9)の分母の外乱エネルギーが限定される場合、式(9)の分子の2乗推定誤差の総和は有界となり、ある時刻以降の推定誤差が0になる。これは、γをより小さく出来れば、状態xの変化に推定値x^k|kが速やかに追従できることを意味する。
ここで、定理1のハイパーHフィルタのアルゴリズムは通常のHフィルタのものとは異なることに注意されたい。また、γ→∞のとき、ρ=1、G=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(N)の計算量が必要であった。しかし、次の条件を用いれば計算量O(N)で定理1のHフィルタの存在性、すなわち式(9)を検証することができる。
系1:スカラー存在条件
次の存在条件を用いれば計算量O(N)でハイパーHフィルタの存在性が判定できる。
JP0004444919B2_000011t.gif ここで、
JP0004444919B2_000012t.gifただし、Ks,iは式(12)で求めたフィルタゲインである。
(証明)
以下に系1の証明を説明する。
2×2の行列Re,kの特性方程式
JP0004444919B2_000013t.gifを解けば、Re,kの固有値λが次のように得られる。
JP0004444919B2_000014t.gif もし、
JP0004444919B2_000015t.gifであれば、行列Re,kの2つの固有値の1つは正となり、もう1つは負となり、行列RとRe,kは同じイナーシャをもつ。これより、
JP0004444919B2_000016t.gifを用いれば、式(18)の存在条件が得られる。ここで、Hs,kの計算量はO(N)である。
4.数値的に安定な状態推定アルゴリズム
上述のハイパーHフィルタは、Σ^k|k-1∈RN×Nを更新するため、単位時間ステップ当たりの計算量はO(N)となる、すなわち、Nに比例する算術演算が必要となる。ここで、Nは状態ベクトルxの次元である。よって、xの次元が増加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列Σ^k|k-1は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度(例:32bit)でも動作可能な状態推定アルゴリズムの開発が望まれる。
そこで、次に、
=R1/2T/2
e,k=R1/2e,kT/2e,k
Σ^k|k-1=Σ^1/2k|k-1Σ^T/2k|k-1
に着目して、数値的に安定化した定理1のHフィルタ(平方根アレイアルゴリズム)を定理2に示す。ただし、ここでは簡単のためF=Iとしたが、F≠Iの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパーHフィルタを示す。
(定理2)
JP0004444919B2_000017t.gif ただし、
JP0004444919B2_000018t.gifJP0004444919B2_000019t.gif位行列である。また、K(:,1)は行列Kの1列目の列ベクトルを表す。
なお、式(21)、(22)において、J-1およびJは削除可能である。
図4に、定理2の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図2に示した定理1のフローチャート中のHフィルタの計算(S105)で用いることができる。
本推定アルゴリズムは、Σ^k|k-1をリカッチ型の更新式で求める代わりに、その因数行列Σ^1/2k|k-1∈RN×N(Σ^k|k-1の平方根行列)をJ-ユニタリ変換に基づく更新式で求めている。このとき生じる1-1ブロック行列と2-1ブロック行列からフィルタゲインKs,kを図示のように求めている。このため、Σ^k|k-1=Σ^1/2k|k-1Σ^T/2k|k-1>0となり、Σ^k|k-1の正定性は保証され、数値的に安定化できる。なお、定理2のHフィルタの単位ステップ当たりの計算量はO(N)のままである。
なお、図4において、J-1は削除可能である。
まず、処理部101は、式(22)の左辺の行列式の各要素に含まれる項を記憶部105から読み出し又は内部メモリ等から得て、J-ユニタリ変換を実行する(S301)。処理部101は、求めた式(22)の右辺の行列式の要素からシステムゲインK、Ks,kを式(21)に基づき計算する(S303、S305)。処理部191は、式(20)に基づき状態推定値x^k|kを計算する(S307)。
5.状態推定のための数値的に安定な高速アルゴリズム
上述のように、定理2のHフィルタの単位ステップ当たりの計算量はO(N)のままである。そこで、計算量の対策として、k+1Ψ,=[u(k),・・・,u(0),0,・・・,0]のとき、=[x,0の1ステップ予測誤差の共分散行列Σk+1|k
JP0004444919B2_000020t.gifを満たすことを利用して、Σk+1|k代わりに次元の低い(すなわちL)を更新する
JP0004444919B2_000021t.gif3が得られる。
(定理3)
JP0004444919B2_000022t.gifJP0004444919B2_000023t.gif ただし、
JP0004444919B2_000024t.gifなお、定理3の証明は、後述する。
上式は、K(=ρ-1/2)の代わりにKについても整理することができる。
さらに、次のJ-ユニタリ行列
JP0004444919B2_000025t.gifを用いれば定理4の高速化した状態推定アルゴリズムが得られる。ただし、Ψはシフト行列を表す。
(定理4)
JP0004444919B2_000026t.gif ただし、
JP0004444919B2_000027t.gifであり、diag{・}は対角行列、Re,k+1(1,1)は行列Re,k+1の1-1成分をそれぞれ表す。また、上式はKの代わりにKに関しても整理できる。
本高速アルゴリズムは、次の因数分解
JP0004444919B2_000028t.gifにおけるL∈R(N+1)×2の更新によってフィルタゲインKs,kを求めているので、単位ステップ当たりの計算量はO(N+1)で済む。ここで、次式に注意されたい。
JP0004444919B2_000029t.gif 図5に、定理の数値的に安定な高速アルゴリズムのフローチャートの一例を示す。この高速アルゴリズムは図2のHフィルタの計算ステップ(S105)に組み込まれ、γ-イタレーションによって最適化される。よって、存在条件が満たされる間はγは除々に減少されるが、満たされなくなった時点で、図示のようにγは増加される。
フィルタリングアルゴリズムは以下のように要約することができる。
[ステップ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で存在条件を満たさなければγを増加し、ステップS401に戻る。
[ステップS411]処理部101は、式(25)のHフィルタのフィルタ方程式を更新する。
[ステップS415]処理部101は、Rr.k+1を式(30)を用いて再帰的に計算する。また、処理部101は、Lk+1を式(28)、(31)を用いて再帰的に計算する。
[ステップS419]処理部101は、時刻kを進ませて(k=k+1)、ステップS403に戻り、データがある限り続ける。
なお、処理部101は、Hフィルタ計算ステップS405~S415及び存在条件の計算ステップS409等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部105に記憶し、また、記憶部105から読み出すようにしてもよい。
6.エコーキャンセラ
つぎに、エコーキャンセリング問題の数理モデルを作成する。
まず、受信信号{u}がエコーパスへの入力信号となることを考慮すれば、エコーパスの(時変)インパルス応答{h[k]}により、エコー{d}の観測値{y}は次式で表される。
JP0004444919B2_000030t.gifここで、u,yはそれぞれ時刻t(=kT;Tは標本化周期)における受信信号とエコーを表し、vは時刻tにおける平均値0の回線雑音とし、h[k],i=0,・・・,N-1は、時変インパルス応答であり、そのタップ数Nは既知とする。このとき、インパルス応答の推定値{h^[k]}が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
JP0004444919B2_000031t.gifこれをエコーから差し引けば(y-d^≒0)、エコーをキャンセルすることができる。ただし、k-i<0のときuk-1=0とする。
以上より、問題は直接観測可能な受信信号{u}とエコー{y}からエコーパスのインパルス応答{h[k]}を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにHフィルタを適用するには、まず式(32)を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、問題がインパルス応答{h[k]}を推定することであるから、{h[k]}を状態変数xとし、w程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
JP0004444919B2_000032t.gif ただし、
JP0004444919B2_000033t.gif このような状態空間モデルに対するハイパーおよび高速Hフィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
7.インパルス応答に対する評価
(動作の確認)
エコーパスのインパルス応答が時間的に不変であり(h[k]=h)、かつそのタップ数Nが48である場合について、シミュレーションを用いて、本高速アルゴリズムの動作を確認する。
JP0004444919B2_000034t.gifなお、図6は、ここでのインパルス応答{h}の値を示す図である。
ここで、インパルス応答{hi=023は、図示の値を採用し、その他{h2447は0とする。また、υは平均値0、分散σ=1.0×10-6の定常なガウス白色雑音とし、標本化周期Tを便宜上1.0とする。
また、受信信号{u}は次のように2次のARモデルで近似する。
JP0004444919B2_000035t.gifただし、α=0.7,α=0.1とし、w′は平均値0、分散σw=0.04の定常なガウス白色雑音とする。
(インパルス応答の推定結果)
図7に、定理4の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図7(b)の縦軸は、
√{Σi=047(h-x^(i+1))
を表す。
これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし、ρ=1-χ(γ)、χ(γ)=γ-2、x^0|0=0、Σ^1|0=20Iとし、計算は倍精度で行った。また、存在条件を確認しつつ、γ=5.5と設定とした。
8.定理の証明
8-1.定理2の証明
次の関係式
JP0004444919B2_000036t.gifが成り立つとき、両辺の2×2ブロック行列の各項を比較すれば次式が得られる。
JP0004444919B2_000037t.gifこれは定理1のF=Iのときの式(13)のリカッチ方程式と一致する。ただし、
JP0004444919B2_000038t.gif一方、AJA=BJBが成り立つとき、BはJ-ユニタリ行列Θ(k)を用いてB=AΘ(k)と表すことができる。よって、式(40)より定理1のリカッチ方程式は次式と等価である。
JP0004444919B2_000039t.gif なお、式(40)、(45)において、J-1は削除可能である。
8-2.定理3の証明
次のようにブロック三角化するJ-ユニタリ行列Θ(k)が存在すると仮定する。
JP0004444919B2_000040t.gifのように決定することができる。ただし、Sは対角成分が1または-1をとる対角行列とする。
(1,1)-ブロック行列
JP0004444919B2_000041t.gif(2,1)-ブロック行列
JP0004444919B2_000042t.gif(2,2)-ブロック行列
JP0004444919B2_000043t.gif これより、
JP0004444919B2_000044t.gif8-3.定理4の証明
観測行列Hがシフト特性をもち、かつ
JP0004444919B2_000045t.gifのとき、定理2と同様な方法によって次の関係式が得られる。
JP0004444919B2_000046t.gif ただし、
JP0004444919B2_000047t.gifとなるようにRr,k+1を決定する。次に、式(46)の3行目にRr,k+1の更新式を新たに追加すれば、最終的に次式が得られる。
JP0004444919B2_000048t.gifこの両辺の3×2ブロック行列の各項の対応から次のゲイン行列Kの更新式が得られる。
JP0004444919B2_000049t.gif【産業上の利用可能性】
一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができる。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9