TOP > 国内特許検索 > システム同定方法 > 明細書

明細書 :システム同定方法

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4067269号 (P4067269)
公開番号 特開2002-135171 (P2002-135171A)
登録日 平成20年1月18日(2008.1.18)
発行日 平成20年3月26日(2008.3.26)
公開日 平成14年5月10日(2002.5.10)
発明の名称または考案の名称 システム同定方法
国際特許分類 H04B   3/23        (2006.01)
G05B  13/02        (2006.01)
H03H  21/00        (2006.01)
FI H04B 3/23
G05B 13/02 T
H03H 21/00
請求項の数または発明の数 5
全頁数 36
出願番号 特願2000-323958 (P2000-323958)
出願日 平成12年10月24日(2000.10.24)
新規性喪失の例外の表示 特許法第30条第1項適用 IECON-2000(2000年10月22日)第462~467頁に発表
審査請求日 平成16年4月30日(2004.4.30)
特許権者または実用新案権者 【識別番号】503360115
【氏名又は名称】独立行政法人科学技術振興機構
発明者または考案者 【氏名】西山 清
個別代理人の代理人 【識別番号】100107010、【弁理士】、【氏名又は名称】橋爪 健
審査官 【審査官】東 昌秋
参考文献・文献 Nishiyama,Robust estimation of a single complex sinusoid in white noise - H∞ filtering approach,IEEE Transactions on Signal Processing,米国,IEEE,1999年10月,vol.47,No.10,pp.2853-2856
西山 清,高速ロバスト学習アルゴリズム,電子情報通信学会技術研究報告,NC(ニューロコンピューティング),日本,電子情報通信学会,1999年 3月19日,vol.98,no.674,pp.187-192
調査した分野 H04B 3/20-3/23
G05B 13/02
H03H 21/00
特許請求の範囲 【請求項1】
時不変又は時変システムの高速実時間同定を行うシステム同定方法において、次式で表されるHフィルタ方程式を用い、
【数1】
JP0004067269B2_000046t.gif
評価基準として、{ (フィルタ誤差に対応する項/評価関数の重み(ρ)) / [(初期状態の誤差に対応する項)+(システム雑音に対応する項)+(観測雑音に対応する項/評価関数の重み(ρ))] }の最大値が、予め与えられた上限値(γ)より小さく抑えるように定めることにより、外乱に対して頑健なフィルタリングアルゴリズムとしたシステム同定方法。
ただし、
: 状態ベクトルまたは単に状態 ; 未知、これが推定の対象となる。
: 初期状態 ; 未知である。
: システム雑音 ; 未知である。
: 観測雑音 ; 未知である。
: 観測信号 ; フィルタの入力となり、既知である。
: 出力信号 ; 未知である。
: 観測行列 ; 既知である。
x^k|k: 観測信号y~yまでを用いた時刻kの状態xの状態推定値 ; フィルタ方程式(12)によって与えられる。
x^0|0:状態の初期推定値 ; 本来未知であるが、便宜上0が用いられる。
s,k+1:フィルタゲイン ;行列P^k+1|kから得られる。
Σwk : システム雑音の共分散行列に対応 ; 理論上は既知であるが、フィルタの実行前には未知である。
P^k|k-1: x^k|k-1の誤差の共分散行列に対応 ; リカッチ方程式によって与えられる。
P^1|0: 初期状態の誤差の共分散行列に対応 ; 本来未知であるが、便宜上εIが用いられる。
さらに、
「フィルタ誤差」:
(11)式のzk|kと(9)式のz=Hの差((15)式のef,iに対応)の重み付きノルムの2乗の和。
「初期状態の誤差」:
時刻k=0のときの、以下に示す状態空間モデルの、(7)式の状態xkの初期値xとその推定値の差の重み付きノルム。
「システム雑音に対応する項」:
以下に示す状態空間モデルの、(7)式のシステム雑音wの重み付きノルムの2乗の和。
「観測雑音に対応する項」:
以下に示す状態空間モデルの、(8)式の観測雑音vの重み付きノルムの2乗の和。
状態空間モデル
【数2】
JP0004067269B2_000047t.gif

【請求項2】
さらに、次式により時刻kの状態推定値x^k|kから出力信号を求めるようにした請求項に記載のシステム同定方法。
【数3】
JP0004067269B2_000048t.gif

【請求項3】
ゲイン行列K、補助変数A、S、D、状態推定値x^k|kの再帰式の初期条件をそれぞれ以下のように定めるステップと、
【数4】
JP0004067269B2_000049t.gif
時刻kにおける補助変数AとSを以下のように再帰的に決定するステップと、
【数5】
JP0004067269B2_000050t.gif
ゲイン行列Kに補助変数を含む行を増やした第2のゲイン行列Kを以下のように計算するステップと、
【数6】
JP0004067269B2_000051t.gif
第2のゲイン行列Kを以下のように分割し、第1及び第2の分割ゲイン行列を求めるステップと、
【数7】
JP0004067269B2_000052t.gif
第1及び第2の分割ゲイン行列を含む次式により、Dを決定し、時刻k+1におけるゲイン行列Kk+1を求め、時刻k+1におけるフィルタゲインKs,k+1を求めるステップと、
【数8】
JP0004067269B2_000053t.gif
ここで、η∈R2×1,D∈RN×1
k+1∈RN×2,Ks,k+1∈RN×1
0<ρ=1-γ-2≦1,γ>1
上記で求められたフィルタゲインKs,k+1に基づき、前記フィルタ方程式(12)を以下のように更新するステップと、
【数9】
JP0004067269B2_000054t.gif
時刻を進ませて、各前記ステップを繰り返すためのステップと
を含む請求項に記載のシステム同定方法。
【請求項4】
さらに、高速処理に適した存在条件として、次式を用いることにより、計算量Ο(N)で前記高速Hフィルタの存在性を検査することを特徴とする請求項1乃至3のいずれかに記載のシステム同定方法。
【数10】
JP0004067269B2_000055t.gif

【請求項5】
前記Hフィルタ方程式を適用し、状態推定値x^k|kを求め、擬似エコーを次式のように推定し、求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現することを特徴とする請求項1又は3又は4に記載のシステム同定方法。
【数11】
JP0004067269B2_000056t.gif
発明の詳細な説明
【0001】
【発明の属する技術分野】
本発明は、システム同定方法に係り、特に、新たなH評価基準に基づいて開発された変形Hフィルタの高速Hフィルタリングアルゴリズムを用いて時変システムを高速に実時間同定するものである。
【0002】
【従来の技術】
一般に、システム同定(system identification)とは、入出力データに基づいてシステムの入出力関係の数理モデル(伝達関数、あるいはインパルス応答など)を推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、1993年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。
【0003】
(基本原理)
図14に、システム同定のための構成図を示す。このシステムは、未知システム1、適応フィルタ2を備える。また、適応フィルタ2は、FIRディジタルフィルタ3、適応アルゴリズム4を有する。
以下に、未知システム1を同定する出力誤差方式の一例を説明する。ここで、uは未知システム1の入力、dは所望信号であるシステムの出力、d^はフィルタの出力である。(なお、「^」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。)
【0004】
未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差e=d-d^を最小にするように適応アルゴリズムによってFIRディジタルフィルタ3の係数を調節する。
【0005】
図15に、インパルス応答の調節機構についての構成図を示す。
ここで、適応アルゴリズムの一例として、その簡便さより、次のLMSアルゴリズム(least mean square algorithm)が広く用いられている。
【0006】
【数6】
JP0004067269B2_000002t.gif
【0007】
また、一般に、時変システムの同定には、例えば、収束性が早いカルマンフィルタが適している。
【0008】
【数7】
JP0004067269B2_000003t.gif
【0009】
ここで、求めるべきインパルス応答{h}は状態推定値x^k|kとして得られ、システムへの入力{u}は観測行列Hの要素として用いられている。
また、σ=0のときのカルマンフィルタに対して、観測行列Hのシフト特性(Hk+1(i+1)=H(i))を利用して、単位時間ステップ当たりの計算量をNに比例した演算回数、すなわちO(N)までに軽減した高速カルマンフィルタリングアルゴリズムが知られている。詳細は、1993年電子情報通信学会「ディジタル信号処理ハンドブック」など参照。
【0010】
(エコーキャンセラへの適用例)
国際電話など長距離電話回線では,信号増幅などの理由から4線式回線が用いられている。一方、加入者回線は比較的短距離なので、2線式回線が使用されている。
図16に通信系とエコーについての説明図を示す。2線式回線と4線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者Bからの信号(音声)は話者Aのみに到達する。しかし、一般に整合を完全とするのはむずかしく、受信信号の一部は4線式回線に漏れ、増幅された後、再び受信者(話者A)に戻ると云った現象が起こる。これがエコー(echo)である。エコーは、伝送距離が長くなるにつれて(遅延時間が長くなるにつれて)影響が大きくなり、著しく通話の品質を劣化させる(パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する)。
【0011】
図17に、エコーキャンセラの原理図を示す。
そこで、図示のようにエコーキャンセラ(echo canceller) を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。
【0012】
エコーパスのインパルス応答の推定は、残留エコーeの平均2乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者Aからの信号(音声)である。一般に、話者2人が同時に話し始めた(ダブルトーク)ときはインパルス応答の推定を中断する。また、ハイブリッドトランスのインパルス応答長は50[ms]程度なので、サンプリング周期を125[μs]とするとエコーパスのインパルス応答の次数は実際は400程度となる。
【0013】
【発明が解決しようとする課題】
従来技術では、適応アルゴリズムとして、その簡便さより、LMSアルゴリズム(least mean square algorithm)が広く用いられて来たが、収束性が非常に遅いため急激に変化するような時変システムの同定は不可能であった。
【0014】
また、追従性の優れたカルマンフィルタは、計算量がO(N)あるいはO(N)であり、それが、タップ数Nと共に急速に増加してしまい、高いタップ数Nが要求される現実の問題の実時間処理は困難であった。この対策として、タップ数Nに対して単位時間ステップ当たり計算量O(N)で同定可能な高速カルマンフィルタが提案されているが、その定常的な特性(システム雑音が考慮できない点)から時変システムの同定は不可能であった。
【0015】
本発明は、以上の点に鑑み、新たなH評価基準に基づいて開発した変形Hフィルタの高速Hフィルタリングアルゴリズムを用いて、時不変および時変システムの高速実時間同定および推定を実現することを目的とする。また、本発明は、本アルゴリズムの特殊な場合として高速カルマンクイルタリングアルゴリズムを含み、また、時変システムの追従性を支配するシステム雑音の共分散を理論的に決定することを目的とする。また、本発明は、突然回線が切り替わるような激しく変化する時変システムのエコーキャンセラなどのように、入力信号が不連続に変化する場合においても、非常に有効な高速時変システム同定方法を提供することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム同定方法を提供することを目的とする。
【0016】
【課題を解決するための手段】
本発明では、上述の課題を解決するために、新たにH評価基準を考案し、これに基づく変形Hフィルタの高速アルゴリズムを開発すると共に、この高速Hフィルタリングアルゴリズムに基づく高速時変システム同定方法を提案する。本発明による高速アルゴリズムは、単位時間ステップ当たり計算量O(N)で急激に変化する時変システムの追跡が可能である。また、γ=∞の極限で高速カルマンフィルタリングアルゴリズムと完全に一致すると云った便利な特性をもっている。
【0017】
【発明の実施の形態】
以下に、本発明の実施の形態について説明する。なお、詳細は、例えば、K. Nishiyama :"Derivation of A Fast Algorithm of Modified H Filters", IEEE international Conference on Industrial Electronics, Control and Instrumentation, October, 2000 に示されている。
【0018】
1.記号の説明
まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。
: 状態ベクトルまたは単に状態 ; 未知、これが推定の対象となる。
: 初期状態 ; 未知である。
: システム雑音 ; 未知である。
: 観測雑音 ; 未知である。
: 観測信号 ; フィルタの入力となり、既知である。
: 出力信号 ; 未知である。
: 観測行列 ; 既知である。
: 出力行列 ; 既知である。
x^k|k: 観測信号y~yまでを用いた時刻kの状態xの状態推定値; フィルタ方程式によって与えられる。
x^0|0:状態の初期推定値 ; 本来未知であるが、便宜上0が用いられる。
s,k+1:フィルタゲイン ;行列P^k|k-1から得られる。
Σwk : システム雑音の共分散行列に対応 ; 理論上は既知であるが、実際には未知である。
P^k|k-1: x^k|k-1の誤差の共分散行列に対応 ; リカッチ方程式によって与えられる。
P^1|0: 初期状態の誤差の共分散行列に対応 ; 本来未知であるが、便宜上εIが用いられる。
σ2 : 観測雑音の分散 ; 理論上は既知として扱われるが、実際には未知である。
σ2 : システム雑音の分散 ; 理論上は既知として扱われるが、実際には未知である。
【0019】
なお、記号の上に付される”^”は、推定値の意味であり、"∪"は、行列を1行増やしたことを表す。”また、”~”等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載するが、数式で示すように、文字の真上に記載されたものと同一である。また、L、H、P、K等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。
2.変形Hフィルタ
【0020】
つぎに、次式(7)~(9)のような状態空間モデルを考える。
【数8】
JP0004067269B2_000004t.gif
【0021】
ここで、エコーキャンセラなどを想定し、L=H(H=[uk-2・・・ uk-N+1]) とする。このような状態空間モデルに対して、次式(10)のようなH評価基準(新たに左辺にγが入っている)を提案する。
【0022】
【数9】
JP0004067269B2_000005t.gif
【0023】
この評価基準を満たすレベルγの変形Hフィルタは、ρやΣwkがγに依存しないと仮定すれば、システム同定の分野で通常知られる方法を適用することによって、次の式(11)~式(14)で与えられる。なお、その通常知られる方法として、例えば、
B. Hassibi, A. H. Sayed, and T. Kailath: "Linear Estimation in Krein Spaces - Part I: Theory," IEEE Trans. Automatic Control, 41, 1, pp.18-33,1996.、
B. Hassibi, A. H. Sayed, and T. Kailath: "Linear Estimation in Krein Spaces - Part II: Applications," IEEE Trans. Automatic Control, 41, 1, pp.34-49,1996.
等を参照のこと。
【0024】
【数10】
JP0004067269B2_000006t.gif
【0025】
また、評価基準の重みρは予めに決められた上限値γに依存するので、上述のアルゴリズムは通常のHフィルタとは本質的に異なる。本アルゴリズムはρで重みづけされた外乱(初期状態x,システム雑音{w},観測雑音{υ}からフィルタ誤差{ef,i}への最大エネルギーゲインをγより小さく押えているので、外乱に対してロバスト(頑健)なフィルタリングアルゴリズムとなる。この性質が時変システムの追従特性に反映される。また、γ→∞のとき、ρ=1、Σwk=0となり、変形Hフィルタは通常のカルマンフィルタと一致する。
【0026】
変形Hフィルタの主な計算上の負担は、NまたはNに比例する計算量を必要とするP^k+1|k∈RN×Nの更新の際に生じる。すなわち、単位時間ステップ当たりO(N)の算術演算が必要となる。ここで、タップ数Nは状態ベクトルxの次元と一致する。ゆえに、xの次元が増加するにつれて変形Hフィルタの実行に要する計算時間は急速に増大する。この計算上の負担を克服するために変形Hフィルタの高速アルゴリズムの導出が必要となる。
【0027】
3.高速Hフィルタリングアルゴリズム
変形Hフィルタの計算量は、式(14)のリカッチ方程式(誤差の共分散方程式)の計算に支配される。よって、変形Hフィルタを高速に処理するためには、リカッチ方程式を用いずに式(13)のフィルタゲインを直接決定できれば、大幅に計算量を削減できる。
しかし、フィルタゲインKs,k∈RN×1を求める高速アルゴリズムの導出は困難なため、次のように定義されるゲイン行列を高速に計算するアルゴリズムを導出することを考える。
【0028】
【数11】
JP0004067269B2_000007t.gif
【0029】
このとき、次の補題が成り立つ。
補題1
行列Pは式(14)のリカッチ方程式を満たす。これより、ゲイン行列Kが求まれば、次の補題からフィルタゲインKs,kが直ちに得られる。
【0030】
補題2
変形HフィルタのフィルタゲインはKs,kは、ゲイン行列Kを用いて次のように得られる。実際、ゲイン行列Kは次の再帰的方法によって高速に計算できる。
【0031】
【数12】
JP0004067269B2_000008t.gif
【0032】
補題3
ゲイン行列Kは次のように更新される。
【0033】
【数13】
JP0004067269B2_000009t.gif
【0034】
ここで、m∈RN×2とμ∈R1×2は行列K=Q-1の次の分割によって得られる。
【0035】
【数14】
JP0004067269B2_000010t.gif
【0036】
また、補助変数A∈RN×1,S∈RおよびB-1∈RN×1も同様に得られる。
結論として、高速Hフィルタリングアルゴリズムは以下のように要約することができる。
図1に、高速アルゴリズムのフローチャートを示す。なお、Lは最大データ数を示す。
[ ステップ0] 再帰式の初期条件を以下のようにする。ここでεは充分に大きい正の定数である。
【0037】
【数15】
JP0004067269B2_000011t.gif
【0038】
[ ステップ1] 時刻kと最大データ数Lとを比較する。時刻kが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。(不要であれば条件文を取り除くことができる。)
[ ステップ2] AとSを以下のように再帰的に決定する。
【0039】
【数16】
JP0004067269B2_000012t.gif
【0040】
[ ステップ3] Kを以下のように計算する。
【0041】
【数17】
JP0004067269B2_000013t.gif
【0042】
[ ステップ4] Kを以下のように分割する。
【0043】
【数18】
JP0004067269B2_000014t.gif
【0044】
[ ステップ5] Dを決定し、Kk+1を通して、ゲイン行列Ks,k+1を以下のように得る。
【0045】
【数19】
JP0004067269B2_000015t.gif
【0046】
ここで、η∈R2×1,D∈RN×1
k+1∈RN×2,Ks,k+1∈RN×1
0<ρ=1-γ-2≦1,γ>1 である。
[ ステップ6] Hフィルタのフィルタ方程式を以下のように更新する。
【0047】
【数20】
JP0004067269B2_000016t.gif
【0048】
[ ステップ7] 時刻kを進ませて(k=k+1)、 ステップ2に戻り、データがある限り続ける。
(高速処理に適した存在条件:補題6)
また、次の存在条件を用いれば、計算量Ο(N)で高速Hフィルタの存在性が検査できる。
【0049】
【数21】
JP0004067269B2_000017t.gif
【0050】
4.本高速アルゴリズムの計算量
つぎに、高速Hフィルタリングアルゴリズムの計算量が、変形Hフィルタリングアルゴリズムの計算量と比べて如何に減少するかを考察する。ここで、式の計算量の評価は乗算回数のみに注目し、以下のような方法で計算した。
(J×K行列)×(K×L行列)の乗算回数=J×K×L(回)
ただし、行列やベクトルが3つ以上乗算されるときは特に図に示さない限り左から計算されるものとする。
【0051】
(変形Hフィルタリングアルゴリズムの計算量)
図2及び図3に、変形Hフィルタリングアルゴリズムの各部分の計算量の説明図を示す。ただし、Nはタップ数(状態ベクトルxのディメンジョン)である。ここで、図3(a)においてRe,kからRe,k-1を求めるための計算は無視する。同様に、図2(a)において(Hk+1P^k+1|k+1+1)の部分から(Hk+1P^k+1|k+1+1)-1を求める計算も無視する。
【0052】
図2(a)、図3(a)および図3(b)より、K|k+1、Re,kおよびP^k+1|は計算量がタップ数の2乗に比例することがわかる。よって、変形Hフィルタリングアルゴリズム全体の計算量は単位時間ステップ当りΟ(N)となる。
【0053】
また、図4に、行列計算の順番を変えた場合の計算量の説明図を示す。すなわち、リカッチ方程式において、右辺第2項の部分の行列計算の順番を変えた場合の計算量を図4に示す。
上述の部分の計算量がインパルス応答の次数の3乗に比例するので、P^k+1|の計算量もインパルス応答の次数の3乗に比例する。これに伴い、Hフィルタ全体の計算量もタップ数の2乗から3乗にまで増加する。
【0054】
しかし、いずれのアルゴリズムにせよタップ数の2乗または3乗に比例する計算量をもつので、タップ数の増加と共にフィルタの実行にかかる計算上の負担が著しく増加する。実際、通信工学の分野に利用する場合などはタップ数が例えば400程度であるため、このアルゴリズムでの実用はかなり困難なものとなる。
【0055】
(高速Hフィルタリングアルゴリズムの計算量)
次に、図5及び図6に、高速Hフィルタリングアルゴリズムの計算量の説明図を示す。ここで、図5(b)のKの式においてS-1はSから求められるが、その計算は無視する。同様に、図6(a)のDの式において、[1-μη]から[1-μη]-1を求める計算も無視する。
【0056】
図5及び図6より、本高速アルゴリズム全体を通して計算量が単位時間ステップ当りΟ(N)になっている。よって、高速Hフィルタリングアルゴリズムは計算量はタップ数に比例することがわかる。また、この場合、高速Hフィルタを1回実行するためにかかる計算量(乗算回数)は単位ステップ当り28N+16であり、高速カルマンフィルタの12N+3と比べて約2倍の計算量(乗算回数)が必要である。
【0057】
以上のように、変形Hフィルタリングアルゴリズムではタップ数の2乗または3乗に比例する計算量が、本高速アルゴリズムではタップ数の1乗に比例するまで減少させることができる。
【0058】
5.エコ-キャンセラ
つぎに、エコーキャンセラの例を取り上げ、本発明の効果を検証する。
まず、受信信号{u}がエコーパスへの入力信号となることを考慮すれば、エコーパスの(時変)インパルス応答{h[k]}により、エコー{d}の観測値{y}は次式で表される。
【0059】
【数22】
JP0004067269B2_000018t.gif
【0060】
ここで、u,yはそれぞれ時刻t(=kT;Tはサンプリング周期)における受信信号とエコーを表し、vは時刻tにおける平均値0の回線雑音とし、h[k],i=0,・・・,N-1 は緩やかな変動を想定した時変インパルス応答であり、そのタップ数Nは既知とする。このとき、インパルス応答の推定値{h^[k]}が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。
【0061】
【数23】
JP0004067269B2_000019t.gif
【0062】
これをエコーから差し引けば(y-d^≒0)、エコーをキャンセルすることができる。ただし、k-i<0のときuk-1=0とする。
以上より、問題は直接観測可能な受信信号{u}とエコー{y}からエコーパスのインパルス応答{h[k]}を逐次推定する問題に帰着できる。
一般に、エコーキャンセラにHフィルタを適用するには、まず式(24)を状態方程式と観測方程式からなる状態空間モデルで表現しなければならない。そこで、問題がインパルス応答{h[k]}を推定することであるから、{h[k]}を状態変数xとし、w程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。
【0063】
【数24】
JP0004067269B2_000020t.gif
【0064】
このような状態空間モデルに対する変形および高速Hフィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。
以上より、インパルス応答の推定値{h^[k]}が得られれば、これより疑似エコーが次のように逐次求めることができる。
【0065】
【数25】
JP0004067269B2_000021t.gif
【0066】
よって、これを実際のエコーから差し引き、エコーを打ち消せばエコーキャンセラが実現できる。このとき、推定誤差であるe=y-d^は残留エコーと呼ばれる。
6.時不変インパルス応答に対する評価
【0067】
(推定精度に対する評価)
エコーパスのインパルス応答が時間的に不変であり(h[k]=h)、かつそのタップ数Nが24である場合について、シミュレーションを用いて変形Hフィルタと高速Hフィルタを評価する。
【0068】
【数26】
JP0004067269B2_000022t.gif
【0069】
なお、図7は、ここでのインパルス応答{h}の値を示す図である。また、υは平均値0、分散σ=1.0×10-6の定常なガウス白色雑音とし、サンプリング周期Tを便宜上1.0とする。
また、受信信号{u}は次のように2次のARモデルで近似する。
=αk-1+αk-2+w′ (31)
ただし、α=0.7,α=0.1し、w′は平均値0、分散σw′=0.04の定常なガウス白色雑音とする。
【0070】
ここで、変形Hフィルタと高速Hフィルタとを比較する。
図8に変形Hフィルタと高速Hフィルタによるインパルス応答の推定結果の説明図を示す(初期値 :x^0|0=0,ε=20)。図8(a)、(b)はγ=105のときの両者の推定結果であり、図8(c),(d)はγ=2.0のときの推定結果(x^100100である。これより、両者の推定精度に対する性能は等しいことがわかる。すなわち、高速化することで推定精度の低下は生じない。ここで、γが小さすぎるとフィルタの存在条件を満たさないので注意が必要である。また、γ=1.0×105の場合は高速カルマンフィルタの結果とほぼ一致していた。以上より、高速Hフィルタリングアルゴリズムは、高速カルマンフィルタリングアルゴリズムを含んでおり、かつγを調節することによって、収束を早くできることがわかる。
【0071】
(計算時間の評価)
つぎに、エコーパスのインパルス応答が時間的に不変であり、かつタップ数を24,48,96,192,384と増加させたときの変形Hフィルタと高速Hフィルタの計算時間を評価する。ただし、1回の測定では結果にばらつきがあるので、4回測定した平均を結果として用いた。また、シミュレーションに用いるインパルス応答{h}は図7の値とし、それ以降(24≦k<N)のインパルス応答{h}は0とする。ただし、フィルタの実行はステップ数100までとする。計算時間は、ワークステーション(sparc,60MHz, 32MB)上のMATLABのコマンドetimeによって計測された。
【0072】
図9に、計算時間の測定結果の図を示す。ここで、変形Hフィルタ(2)はリカッチ方程式において、計算量がタップ数の2乗に比例するように行列計算を行ったものであり、変形Hフィルタ(1)は計算量がタップ数の3乗に比例する行列計算を行ったものである(図3(b)および図4参照)。このように、変形Hフィルタは行列計算の順序によって計算量がタップ数の2乗または3乗に比例することになるが、いずれにせよ実用的ではない。
【0073】
7.時変インパルス応答に対する評価
(追従性能の評価)
システム(インパルス応答)が時間的に変化したときの各アルゴリズムの追従性能についてエコーキャンセラの例を用いて評価する。ただし、インパルス応答のタップ数は48とし、{h}は、図7の値を元に図10(a)のように時間的に変化した場合を想定する。ただし、vは平均値0、σ=1.0×10-6の定常なガウス白色雑音とし、便宜上サンプリング周期をT=1とする。また、受信信号{u}は次のように2次のARモデルで近似する。
=αk-1+αk-2+w′ (32)
【0074】
ここで、α=0.7、α=0.1とし、w′は平均値0、分散σw´=0.04の定常なガウス白色雑音とする。
図10及び図11に、各アルゴリズムのシミュレーション結果の図を示す。これは、高速Hフィルタ(高速 HF)、高速カルマンフィルタ(高速 KF)およびLMSアルゴリズム(LMS)の時変システムの追従性能を示すものである。図10(b)はγ=2.0の場合の高速Hフィルタによる推定値であり、図11(a)は高速カルマンフィルタによる推定値である。ただし、高速Hフィルタの初期値はx^0|0=0,ε=20とし、高速カルマンフィルタの初期値も同様に設定した。また、図11(b)にLMSアルゴリズムによる推定値を示す。ただし、初期値はh^=0、ステップサイズは安定かつ早い収束を与えるようにμ=0.5とした。これより、高速Hフィルタの追従性能が飛躍的に優れており、インパルス応答の変化後約30ステップで推定値が安定していることがわかる。一方、高速カルマンフィルタとLMSアルゴリズムについては全く追従できていない。
【0075】
一般に、システム雑音を伴わないHフィルタの追従性が低下する原因は、P^k|k—1の対角成分の減少によりフィルタゲインの値が小さくなり推定値の更新量が減少するからである。つまり、ステップ数が経過するとほとんど推定値を更新しなくなる。よって、カルマンフィルタやHフィルタの追従性を向上させるためには、行列のP^k|k—1の対角成分の値に外部から適当な値を加えてやれば良い。しかし、直接導入したのでは観測行列Hのシフト特性を利用した高速アルゴリズムを導出できない。この問題を重みρ=1-γ-2をH評価基準に導入することによって理論的に解決したことが本発明の大きな特徴のひとつである。この重みρは高速HフィルタリングアルゴリズムのSの更新式の中に次のように現れる。
【0076】
(高速Hフィルタの補助変数Sの更新)
高速Hフィルタの補助変数Sは、次式の通りである。
=ρSk-1+e,0<ρ=1-γ-2≦1
また、高速Hフィルタリングアルゴリズムにおいて、Sは、Kuの式でS-1の形で用いられる。よって、より大きな値の更新を行うためには、よりS-1が大きくなければならない。つまりSを小さく保ち大きな値の更新を保持する必要がある。ρの存在はSの急激な増大を防ぎ、結果的にシステム雑音を付加することと等価となり、追従性の向上につながっている。また、重みρはρ=1-γ-2で定義されているので、シミュレーションで確認されたとおりγを変化させることで追従性を変化させることができる。
【0077】
図12に、γとρの関係図を示す。これによるとγ=3.0のときρ=0.8889なのでSk-1の89%がSに伝えられることになる。しかし、あまりγを小さくしすぎるとSk-1の影響が著しく低下すると同時にフィルタの存在条件を満たさなくなるため、注意が必要である。また、γが大きい場合はρ1となりSの増大を全く抑制しないため追従性が低下する。特にγf=∞(ρ=1)のとき本高速アルゴリズムは高速カルマンフィルタリングアルゴリズムと完全に一致する。
【0078】
(計算時間の評価)
図13に、高速Hフィルタ、高速カルマンフィルタおよびLMSアルゴリズムにおけるインパルス応答のタップ数(tap number)と計算時間[s]の関係図を示す。なお、フィルタの実行ステップ数:300, γ=3.0とした。そして、高速Hフィルタ、高速カルマンフィルタリングアルゴリズムおよびLMSアルゴリズムに対して、図10及び図11の例においてタップ数を48,96,192,384と増加させたときの計算時間を計測した。ただし、1回の測定では結果にばらつきがあるので、4回測定しその平均をとった。
【0079】
いずれのアルゴリズムも計算量がタップ数の1乗に比例することが確認できる。また、タップ数が多い場合、高速Hフィルタリングアルゴリズムの計算時間は高速カルマンフィルタリングアルゴリズムと比べて、約2倍弱であり、実用的なLMSアルゴリズムと比較しても4倍程度であることがわかる。追従性を考慮すれば、高速Hフィルタリングアルゴリズムの有効性は十分に高いと言えよう。
【0080】
8.補題の証明
ここで、上述の補題の証明について説明する。
(補題1の証明)
の逆行列をとれば、式(33)となる。さらに、逆行列の補助定理を用いれば、式(34)に示すように、行列Pに関する再帰式が得られる。
【0081】
【数27】
JP0004067269B2_000023t.gif
【0082】
ここで、PをP^k+1|kと見倣せば、上式が式(1)のリカッチ方程式を満たすことがわかる。
(補題2の証明)
ゲイン行列Kが次のように整理できる。
【0083】
【数28】
JP0004067269B2_000024t.gif
【0084】
さらに、G=(ρ+Hk-1)/(1+Hk-1)と、H=Hk-1/(1+Hk-1)を用いれば、ゲイン行列Kの第1ブロック列から式(18)のようにフィルタゲインを得ることができる。
【0085】
(補題3の証明)
ゲイン行列K,i=0・・・,k が与えられたと仮定し、次のKk+1を求めることにする。
k+1k+1=Ck+1 (36)
まず、Cのシフト特性を利用するため、新たに式(37)と式(38)を導入する。このとき、Qは式(39)のように再帰的に表され、かつ次の式(40)のように分割される。
【0086】
【数29】
JP0004067269B2_000025t.gif
【0087】
この表記を用いれば、時間ステップkおよびk+1の式(36)は、次式に含まれる。
【0088】
【数30】
JP0004067269B2_000026t.gif
【0089】
これらの表記に基づけば、Kを直接求める変わりに、次式を満たすK∈R(N+1)×2を求める方が便利である。
【0090】
【数31】
JP0004067269B2_000027t.gif
【0091】
そのため、式(41)から得られる式(45)を用いれば、K∈R(N+1)×2を式(46)のように整理できる。
【0092】
【数32】
JP0004067269B2_000028t.gif
【0093】
ここで、Kは、m∈RN×2とμ∈R1×2に分割される。また、α-c=-(c+A)に注意されたい。さらに、Qは逆行列が存在すると仮定し、補助変数A∈RN×1とS∈Rは次式を満たす。
【0094】
【数33】
JP0004067269B2_000029t.gif
【0095】
ここで、上式の下ブロックはT+Q=0またはT=-Aを与える。
次に、Cの上ブロックに影響することなく、式(46)のμを消去するために、次の式(48)のような補助変数B∈RN×1とF∈Rを導入し、さらに式(46)のKの分割からB-1μを引けば、式(49)を得る。
【0096】
【数34】
JP0004067269B2_000030t.gif
【0097】
さらに、左から式(49)の左辺にQを掛ければ、次のように整理できる。
【0098】
【数35】
JP0004067269B2_000031t.gif
【0099】
上式の左辺に式(49)を代入すれば、式(43)は次式で表される。
【0100】
【数36】
JP0004067269B2_000032t.gif
【0101】
これは式(42)と同じ形であり、式(51)の上ブロックから次式(52)を導くことができる。
k+1(m-B-1μ)=Ck+1 (52)
ここで、式(36)と式(52)を比較すれば、ゲイン行列Kの更新式を得ることができる。
【0102】
(補題4)
補助変数AとSは次のように得られる。
【0103】
【数37】
JP0004067269B2_000033t.gif
【0104】
ただし、A-1=0、S-1=1/εである。
(証明)AとSの式(55)と式(39)を用いると、式(56)を得る。
【0105】
【数38】
JP0004067269B2_000034t.gif
【0106】
一方、式(41)の両辺にW[C+Ck-1]を掛ければ、次式を得る。
【0107】
【数39】
JP0004067269B2_000035t.gif
【0108】
式(56)から式(57)を引けば、次式(58)が成り立つ。
【0109】
【数40】
JP0004067269B2_000036t.gif
【0110】
これを式(47)と比較すれば、α=T=-Aより、式(53)と式(54)を得る。
【0111】
(補題5)
補助変数D=B-1は、次式(59)のように得られる。また、Fは、次式(60)で更新される。
【0112】
【数41】
JP0004067269B2_000037t.gif
【0113】
ただし、η=Ck-1=ck—N+Ck+1k-1, D-1=0,F-1=0 である。
(証明)BとFを更新するため、式(61)を用いれば、式(62)が成り立つ。
【0114】
【数42】
JP0004067269B2_000038t.gif
【0115】
上式を式(61)と同じ形に変形するため、式(62)から Ck-1 を引けば、次式を得る。
【0116】
【数43】
JP0004067269B2_000039t.gif
【0117】
この最後の式と式(48)を比較すれば、Bに関する再帰式を得る。
【0118】
【数44】
JP0004067269B2_000040t.gif
【0119】
これより、BとFが更新される。
しかし、それらはD=B-1あるいは=B-1∈RN×1としてだけ用いられるので、式(48)と式(64)を式(65)を用いて書き換えた方が便利である。また、この行列Dは式(66)を満たす。
【0120】
【数45】
JP0004067269B2_000041t.gif
【0121】
次に、式(63)にFk-1-1を掛ければ、式(67)となり、さらにDk-1=Bk-1k-1-1を用いれば次式(68)のように整理できる。
【0122】
【数46】
JP0004067269B2_000042t.gif
【0123】
これより、式(68)に[1-μk-1-1を掛ければ、次式が得られる。
【0124】
【数47】
JP0004067269B2_000043t.gif
【0125】
これを式(66)と比較すれば、最終的にDとFの更新式が得られる。
(補題6(高速処理に適した存在条件)の証明)
上述したように、式(22)、(23)の存在条件を用いれば、計算量Ο(N)で高速Hフィルタの存在が検査できる。その証明を以下に示す。
次式(69)で示す2×2の行列Re,kの特性方程式を解けば、Re,kの固有値λが次式(70)のように得られる。
【0126】
【数48】
JP0004067269B2_000044t.gif
【0127】
もし、次式(71)が成り立てば、行列Re,kの2つの固有値の1つは正となり、もう1つは負となり、行列RとRe,kは同じイナーシャをもつ。これより、次式(72)を用いれば、式(22)の存在条件が得られる(即ち、式(71))。ここで、Hの計算がΟ(N)回の掛け算を必要としている。
【0128】
【数49】
JP0004067269B2_000045t.gif
【0129】
【発明の効果】
本発明によると、以上のように、新たなH評価基準に基づいて開発した変形Hフィルタの高速アルゴリズム(高速Hフィルタリングアルゴリズム)を用いて、時不変および時変システムの高速実時間同定および推定を実現することができる。また、本発明によると、本アルゴリズムの特殊な場合として高速カルマンクフィルタリングアルゴリズムを含み、また、時変システムの追従性を支配するシステム雑音の共分散に対応する項を理論的に決定することができる。また、本発明によると、突然回線が切り替わるような激しく変化する時変システムのエコーキャンセラなどのように、システム(インパルス応答)が時間的に不連続に変化する場合において、特に、非常に有効な高速時変システム同定方法を提供することができる。また、本発明によると、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができるシステム同定方法を提供することができる。
【図面の簡単な説明】
【図1】 高速アルゴリズムのフローチャート。
【図2】 変形Hフィルタリングアルゴリズムの各部分の計算量の説明図(1)。
【図3】 変形Hフィルタリングアルゴリズムの各部分の計算量の説明図(2)。
【図4】 リカッチ再帰法の計算量の説明図。
【図5】 高速Hフィルタリングアルゴリズムの計算量の説明図(1)。
【図6】 高速Hフィルタリングアルゴリズムの計算量の説明図(2)。
【図7】 インパルス応答{h}の値を示す図。
【図8】 変形Hフィルタと高速Hフィルタによるインパルス応答の推定結果の比較説明図。
【図9】 計算時間の測定結果の図。
【図10】 各アルゴリズムのシミュレーション結果の図(1)。
【図11】 各アルゴリズムのシミュレーション結果の図(2)。
【図12】 γとρの関係図。
【図13】 高速Hフィルタ、高速カルマンフィルタおよびLMSアルゴリズムにおけるインパルス応答のタップ数(tap number)と計算時間[s]の関係図。
【図14】 システム同定のための構成図。
【図15】 インパルス応答の調節機構についての構成図。
【図16】 通信系とエコーについての説明図。
【図17】 エコーキャンセラの原理図。
【符号の説明】
1 未知システム
2 適応フィルタ
3 FIRディジタルフィルタ
4 適応アルゴリズム
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16