JP6389805B2 - 伝達系推定装置、その方法及びプログラム - Google Patents

Description

本発明は、電気信号や、音、振動が伝播する伝達系の特性を推定する技術に関する。

電気信号や音、振動が伝播する未知の伝達系の特性を推定することは、当該伝達系への任意の入力に対する当該伝達系の出力を予測する上で有益である。伝達系を推定する問題は、有限のパラメータで表現可能なモデルを予め与え、当該パラメータの値を推定する問題として扱うことが多い。予め与えるモデルは、推定対象の伝達系の物理的素性（入力と出力の関係）に基づき選定される。

ここで対象とするのは、未知の伝達系への入力とそれに対する出力とが線形な関係を有し、その伝達系の特性が、FIR(Finite Impulse Response)フィルタによってモデル化できる場合における、当該フィルタ係数の推定の問題とする。

このとき、伝達系への入力信号x(n)に対し、観測される出力信号(以下「観測出力信号」ともいう)y(n)は、

と表すことができ、ここでh₀,h₁,…,h_Lが推定すべき未知の有限長インパルス応答Hの各要素である。ただし、次数Lは0以上の整数であり、nは離散時間(時刻)を表す。また、v(n)は、観測出力信号y(n)に混入し得る外乱信号である。

いま、時刻nにおける各要素h₀,h₁,…,h_Lの各推定値w₀(n),w₁(n),…,w_L(n)を要素としてもつ特性ベクトルを

とする。ただし、^Tは転置を表す。このとき、時刻nからn−K+1までの間に観測されるK個の観測出力信号y(n),y(n-1),…,y(n-K+1)を要素とするK次元の出力信号ベクトル

は、次式により生成される。

なお、

である。有限長インパルス応答Hは便宜上、時不変としているが、実際は、時間変動を伴う場合も対象とされうる。また、Kは、1以上の整数であり、特性ベクトルW(n)を更新する際、各更新ステップにおいて、直近Kサンプルの観測出力信号y(n),y(n-1),…,y(n-k+1)の情報が用いられる。すなわち、特性ベクトルW(n)を用いて生成される推定出力信号の系列X(n)^TW(n)とのK次元の残差信号ベクトル

を用いて、特性ベクトルW(n)を

のように更新できることが知られている（例えば、非特許文献１参照）。ここで、スカラー値μ、行列Λは、それぞれ、収束特性の調整のために用いられ、μはステップサイズ、ΛはK=1のとき正則化係数、K≧2のとき正則化行列と呼ばれる。正則化行列は、Λ=εIのように、K×Kの単位行列Iに正のスカラー値εを乗じることで得られる対角行列であり、非対角要素は零となるように与えられる。また、対角要素を等しくεとせずに、要素毎に異なる値を与えることも、非特許文献１などにおいて提案されている。（式１０）の更新式は、K=1のとき（正則化付きの）学習同定法(NLMS)、K≧2のとき（正則化付きの）アフィン射影アルゴリズム(APA)と呼ばれる。（式１０）における右辺第２項を更新ベクトル

と呼ぶことにする。この更新ベクトルΔW(n)は、K≧2の場合に、K=1の場合とは異なる特徴を有する。すなわち、K=1の場合は、（式１１）の更新ベクトルΔWは、時刻nに対応する入力ベクトル(入力信号の系列)X’(n)と等しい方向が与えられるのに対し、K≧2の場合には、入力ベクトル(入力信号の系列)X’(n)から過去の入力ベクトル(入力信号の系列)X’(n−1),…, X’(n−K+1)と相関のある成分が除去された方向に更新されることとなる。特に、音声信号のように相関の高い入力信号が与えられた場合において、この相関除去の効果により、K=1の場合と比べ、K≧2において、高速な特性推定が可能となることが知られている。

Y.-S. Choi, H.-C. Shin, and W.-J. Song,;"Adaptive Regularization Matrix for Affine Projection Algorithm", IEEE Transactions on Circuits and Systems II: Express Briefs, Vol. 54, No. 12, pp. 1087-1091, 2007.

しかしながら、従来技術では外乱信号の相関については、考慮されておらず、外乱信号が低周波ノイズや音声信号など周期性成分を含み、相関が高い場合には、推定精度が劣化することがある。

本発明は、外乱信号の相関の大きさに応じて、適切な特性ベクトルの推定を可能とする技術を提供することを目的とする。

上記の課題を解決するために、本発明の一態様によれば、伝達系推定装置は、推定対象の未知伝達系の伝達特性の推定値をフィルタ係数として持つ特性ベクトルと、未知伝達系の入力信号の系列とを用いて、未知伝達系の出力信号を推定し、その推定値である推定出力信号を得る出力推定部と、外乱信号が含まれ得る未知伝達系の出力信号である観測出力信号と推定出力信号との差分である残差信号を算出する残差算出部と、特性ベクトルに任意の更新候補ベクトルを加えたベクトルと入力信号の系列とを用いて生成しうる推定出力候補信号の系列を、観測出力信号の系列から差し引いた場合に得られる想定残差信号系列が外乱信号の系列と一致する確率と、任意の更新候補ベクトルの生起確率と、の積を最大とする更新候補ベクトルを更新ベクトルとして選定し、選定された更新ベクトルを用いて、特性ベクトルを更新する更新ベクトル算出部とを含む。想定残差信号系列が外乱信号の系列と一致する確率は、外乱信号の系列の分散共分散行列の非対角項が非零の要素を取りうる多次元正規分布に依存し、任意の更新候補ベクトルの生起確率は、任意の更新候補ベクトルの分散共分散行列の非対角項が全て零となる多次元正規分布に依存する。

上記の課題を解決するために、本発明の他の態様によれば、伝達系推定方法は、出力推定部が、推定対象の未知伝達系の伝達特性の推定値をフィルタ係数として持つ特性ベクトルと、未知伝達系の入力信号の系列とを用いて、未知伝達系の出力信号を推定し、その推定値である推定出力信号を得る出力推定ステップと、残差算出部が、外乱信号が含まれ得る未知伝達系の出力信号である観測出力信号と推定出力信号との差分である残差信号を算出する残差算出ステップと、更新ベクトル算出部が、特性ベクトルに任意の更新候補ベクトルを加えたベクトルと入力信号の系列とを用いて生成しうる推定出力候補信号の系列を、観測出力信号の系列から差し引いた場合に得られる想定残差信号系列が外乱信号の系列と一致する確率と、任意の更新候補ベクトルの生起確率と、の積を最大とする更新候補ベクトルを更新ベクトルとして選定し、選定された更新ベクトルを用いて、特性ベクトルを更新する更新ベクトル算出ステップとを含む。想定残差信号系列が外乱信号の系列と一致する確率は、外乱信号の系列の分散共分散行列の非対角項が非零の要素を取りうる多次元正規分布に依存し、任意の更新候補ベクトルの生起確率は、任意の更新候補ベクトルの分散共分散行列の非対角項が全て零となる多次元正規分布に依存する。

本発明によれば、外乱信号の相関の大きさに応じて、適切な特性ベクトルの推定を可能とするという効果を奏する。

第一実施形態に係る伝達系推定装置の機能ブロック図。 第一実施形態に係る伝達系推定装置の処理フローの例を示す図。

以下、本発明の実施形態について、説明する。なお、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。

＜本実施形態のポイント＞
従来の未知伝達系の推定技術として、入力信号の相関の影響を反映し、伝達系の特性（インパルス応答）の推定精度を改善する手法はあった。しかし、未知伝達系からの出力信号の観測時に混入する外乱信号の相関については考慮されていなかった。本実施形態では、特性ベクトルおよび観測出力信号の振幅値を確率変数とみなし、それらがガウス分布（正規分布）に従うとして、伝達特性の推定問題を分析する。これにより、外乱信号の相関の影響が、従来対角行列により与えられていた正則化行列の非対角要素に現れることを明らかにし、非対角要素に適切な相関係数を持たせた新たな推定技術を考案する。このような構成により、低周波ノイズや音声信号などの周期性や相関の高い性質を有する外乱信号混入時においても、精度高く未知伝達系の特性推定を可能とする。

＜第一実施形態の原理＞
従来技術による特性ベクトルの推定精度が、特に相関の高い外乱信号に対して劣化する原因を明らかにするため、K≧2の場合における（式１０）の更新式の性質について、以下のように、特性ベクトルW(n)および観測出力信号の系列Y(n)の振幅値を確率変数とみなして分析する。

まず、観測出力信号の系列Y(n)が観測された条件のもとでの、更新されるべき特性ベクトルW(n+1)についての事後確率について、ベイズの定理により

が成り立つ。ここで、簡単のため、時刻nについての表記を省略した。

ここで、p(Y|W)は、未知伝達系の特性ベクトルがWであった場合における観測出力信号の系列Yが取り得る値の確率（観測出力信号の系列Yのうち、特性ベクトルWに起因する未知伝達系の出力は与えられたものとして、それ以外の外乱信号の系列Vの取り得る値の確率に相当する、言い換えると、特性ベクトルWに任意の更新候補ベクトルΔW^を加えたベクトル(W+ΔW^)と入力信号の系列Xとを用いて生成しうる推定出力候補信号の系列｛X^T(W+ΔW^)｝を、観測出力信号の系列Yから差し引いた場合に得られる想定残差信号の系列R^が外乱信号の系列Vと一致する確率に相当する）を表し、p(W)は更新前の特性ベクトルWが取りうる値の事前確率（任意の更新候補ベクトルΔW^の生起確率に相当する）である。

いま、外乱信号の系列Vおよび、特性ベクトルWがともにガウス分布に従うとすると、（式１２）は、

と書き表すことができる。ここで、AおよびBは、外乱信号の系列Vおよび特性ベクトルWそれぞれの分散共分散行列の逆行列に相当し、

は、特性ベクトルの各要素の振幅値が取りうる平均値を要素として持つベクトルに相当し、ここでは、１ステップ前に得られる特性ベクトルWを

とみなす。（式１３）について、両辺の対数を取ると、

となる。ここで、Cは特性ベクトルWに依存しない定数項とする。Jの最大値を与える特性ベクトルWは、（式１２）の事後確率の最大値を与える。そこで、

を満足する特性ベクトルWを計算すると

が得られる。（式１６）について、行列A、Bそれぞれを、対角要素の値が全て等しい対角行列として

と近似した上で、（式１０）と（式１６）を対比させると、（式１０）において、μ=1、ε=β/αとすることで、両者が一致することが分かる。つまり、本来、外乱信号の系列V及び特性ベクトルWそれぞれの分散共分散行列の逆行列として与えられる行列A及びBそれぞれを、対角要素が全て等しい値を持つ対角行列によって近似することで、（式１０）が導出されたと解釈することができる。特性ベクトルWの各要素は、未知伝達系の有限インパルス応答の各要素(サンプル)に相当し、各要素(サンプル)間の値は独立性が高いと考えると、特性ベクトルWの各要素間の相関は小さく、このため、行列Bを対角行列とみなすことが妥当な場合も多い。そこで、任意の更新候補ベクトルΔW^の生起確率((式１２)のp(W)に相当する確率)は、任意の更新候補ベクトルΔW^の分散共分散行列の非対角項が全て零となる多次元正規分布に依存するものとする。

一方、行列Aは、外乱信号の系列Vについての分散共分散行列の逆行列に相当する。外乱信号の系列Vの要素間の相関は、例えば、低周波ノイズや音声信号も外乱信号となり得ることから、必ずしも小さいとは限らない。つまり、行列Aを対角行列とみなすことは、外乱信号が各サンプル毎に無相関な信号とみなすことに相当し、このため、（式１０）により実現される従来技術では、外乱信号の相関が高い場合において、未知伝達系の特性ベクトルWの推定精度が低下すると考えられる。そこで、想定残差信号の系列R^(n)が外乱信号の系列V(n)と一致する確率((式１２)のp(Y|W)に相当する確率)は、外乱信号の系列V(n)の分散共分散行列A^-1の非対角項が非零の要素を取りうる多次元正規分布に依存するものとする。

本実施形態では、外乱信号の相関に応じた特性ベクトルの推定を可能とするため、（式１６）において、行列Aについては、対角化近似を用いずに、

を特性ベクトルの更新式として用いる。

すなわち、従来技術を実現する（式１０）における正則化行列としてεIの代わりに、βA^-1を用いることを特徴とする。ここで、β^-1は、特性ベクトルWを確率変数とみなしたときの各要素の分散が互いに等しいとした場合のその分散の値に相当する。A^-1は、前述のとおり、外乱信号の系列の分散共分散行列に相当し、非対角要素には、外乱信号のサンプル間の相関に応じた値が現れる。

以下、（式１９）により、W(n)を更新する構成について説明する。

＜第一実施形態に係る伝達系推定装置＞
図１は本実施形態による伝達系推定装置２００の機能ブロック図を、図２はその処理フローの例を示す。

伝達系推定装置２００は、フィルタ係数保持部２０１と、出力推定部２０２と、残差算出部２０３と、外乱分散共分散行列保持部２０４と、更新ベクトル算出部２０５とを含む。

伝達系推定装置２００は、推定対象の未知伝達系１の入力信号x(n)と、入力信号x(n)に対する未知伝達系１の出力信号y(n)とを受け取り、特性ベクトルW(n)を出力する。なお、未知伝達系１の伝達特性（有限長インパルス応答）をH={h₀,h₁,…,h_L}とすると、特性ベクトルW(n)は、h₀,h₁,…,h_Lのそれぞれの時刻nにおける推定値w₀(n),w₁(n),…,w_L(n)を要素として持つベクトルである。

＜出力推定部２０２＞
出力推定部２０２は、入力信号x(n)を受け取り、後述するフィルタ係数保持部２０１からフィルタ係数を取り出す。なお、フィルタ係数は、未知伝達系１の伝達特性の推定値を要素として持つ特性ベクトルW(n)の各要素である。言い換えると、特性ベクトルは未知伝達系１の伝達特性の推定値をフィルタ係数として持つ。

出力推定部２０２は、フィルタ係数と入力信号の系列X(n)とを用いて、その畳み込み演算により、未知伝達系１の出力信号を推定し、その推定値である推定出力信号（の系列）D(n)=X(n)^TW(n)を得（Ｓ２０２）、出力する。

＜残差算出部２０３＞
残差算出部２０３は、外乱信号v(n)が含まれ得る未知伝達系１の出力信号である観測出力信号y(n)と、推定出力信号（の系列）D(n)とを受け取る。観測出力信号y(n)の系列Y(n)と推定出力信号（の系列）D(n)との差分である残差信号の系列R(n)を算出し（Ｓ２０３、（式９）参照）、出力する。

＜外乱分散共分散行列保持部２０４＞
外乱分散共分散行列保持部２０４は、外乱信号の分散共分散行列の近似行列であって、非対角要素には外乱信号のサンプル間の相関に応じた値が現れる行列（以下、単に「外乱信号の分散共分散行列」または「行列」ともいう）A^-1を保持する（Ｓ２０４）。

(i)例えば、外乱分散共分散行列保持部２０４は、直近のKサンプルの残差信号の系列R(n)=[r(n),r(n-1),…,r(n-K+1)]^Tを受け取り、次式により、行列A^-1(n)を逐次算出し、保持する。

なお、Eは期待値を表し、実際には、過去に計算した結果との平均処理により計算する。

また、（式２０）の近似行列の計算において、行列の全て要素を計算する必要はなく、対角要素である(k,k)成分（ただしk=1,…,K）に加え、(k,k+1)成分（ただしk=1,…,K-1）や、(k-1,k)成分（ただしk=2,…,K）など、外乱信号の相関の影響が現れやすい対角要素近傍の非対角要素のみを計算して与え、それ以外の非対角要素は零としてもよい。このような構成により演算量を低減することができる。

(ii)また、例えば、外乱分散共分散行列保持部２０４は、リアルタイムの伝達系推定に先立ち、未知伝達系１において、入力信号の系列X(n)を与えていない状態において、外乱信号の系列V(p)=[v(p),v(p-1),…,v(p-K+1)]^Tを予め観測しておき、次式により、行列A^-1を予め算出し、保持しておく。

上述の(i)において、行列A^-1(n)は残差信号の系列R(n)から算出される分散共分散行列から逐次見積もられるといえる。また、上述の(ii)において、行列A^-1はあらかじめ見込値として固定的に与えられるといえる。逐次見積もられる行列A^-1(n)及び固定的に与えられる行列A^-1を併せて行列A^-1と表現する。上述の(i)の方法によって、より柔軟に伝達特性の変化に対応することができる。一方、(ii)の方法により、リアルタイムの処理量を減らすことができる。使用状況に合わせて何れかの方法により、行列A^-1を求め、保持すればよい。

＜更新ベクトル算出部２０５＞
更新ベクトル算出部２０５は、入力信号x(n)と残差信号の系列R(n)とを受け取り、行列A^-1を外乱分散共分散行列保持部２０４から取り出す。更新ベクトル算出部２０５は、行列A^-1と、入力信号の系列X(n)と、残差信号の系列R(n)とを用いて、特性ベクトルW(n)に任意の更新候補ベクトルΔW^(n)を加えたベクトル(W(n)+ΔW^(n))と入力信号の系列X(n)とを用いて生成しうる推定出力候補信号の系列{X^T(n)(W(n)+ΔW^(n))}を、観測出力信号の系列Y(n)から差し引いた場合に得られる想定残差信号の系列R^(n)が外乱信号の系列V(n)と一致する確率((式１２)のp(Y|W)に相当する確率)と、任意の更新候補ベクトルΔW^(n)の生起確率((式１２)のp(W)に相当する確率)と、の積((式１２)のp(Y|W)p(W)に相当)が最大となるよう選定された更新候補ベクトルを更新ベクトルΔW(n)とし、選定された更新ベクトルΔW(n)を用いて、特性ベクトルW(n)を更新し（Ｓ２０５）、更新後の特性ベクトルW(n+1)を出力する。

まず、更新ベクトルΔW(n)を選定する二つの方法について説明する。

(i)更新ベクトル算出部２０５は、特性ベクトルW(n)を確率変数とみなしたときの各要素の分散を全要素について平均した値の逆数を特性ベクトルW(n)に対応する値βとして逐次算出する。例えば、直近のJ個の特性ベクトルW(n),W(n-1),…,W(n-J+1)のそれぞれのl番目の要素w_l(n),w_l(n-1),…,w_l(n-J+1)からl番目の要素の分散1/β_lを求め、その平均値の逆数、または、その逆数の平均値をβとする。すなわち、

または、

また、「特性ベクトルW(n)を確率変数とみなしたときの各要素の分散」は、「特性ベクトルの更新のために加算される更新ベクトルの各要素の大きさ」から近似することができる。例えば、1ステップ前の更新ベクトルΔW(n-1)の各要素をΔw_l(n-1)とすると、次式によりl番目の要素1/β_lを算出することもでき、このように得られた1/β_lを用いで、上述にように、βを逐次算出することもできる。

更新ベクトル算出部２０５は、算出した値βと、入力信号の系列X(n)と残差信号の系列R(n)と行列A^-1を用いて、次式により、更新ベクトルΔW(n)を算出する。

(ii)更新ベクトル算出部２０５は、あらかじめ固定値としてβを与えられる。例えば、更新ベクトル算出部２０５は、リアルタイムの伝達系推定に先立ち、予め未知伝達系１の特性ベクトルを算出しておき、第二特性ベクトルとする。第二特性ベクトルは何らかの方法（例えば本実施形態の方法や非特許文献１の方法）で予め求めておけばよい。更新ベクトル算出部２０５は、第二特性ベクトルを確率変数とみなしたときの各要素の分散を全要素について平均した値の逆数を第二特性ベクトルに対応する値βとして予め算出しておく。つまり、逐次得られる特性ベクトルW(n)に代えて、予め求めておいた第二特性ベクトルを用いて、βを予め算出し、固定値として用いる。

更新ベクトル算出部２０５は、あらかじめ固定値として与えられた値βと、入力信号の系列X(n)と残差信号の系列R(n)と行列A^-1を用いて、（式Ａ）により、更新ベクトルΔW(n)を算出する。

上述の(i),(ii)の処理が、想定残差信号の系列R^(n)が外乱信号の系列V(n)と一致する確率と、任意の更新候補ベクトルΔW^(n)の生起確率との積が最大となる更新候補ベクトルを更新ベクトルΔW(n)として選定する処理に相当する。上述の(i)において、特性ベクトルに対応する値βは、特性ベクトルW(n)の各要素の分散から逐次見積られるといえる。別の言い方をすると、更新ベクトルの各要素の大きさから逐次見積られるといえる。また、上述の(ii)において、第二特性ベクトルに対応する値βは、あらかじめ見込値として固定的に与えられるといえる。上述の(i)の方法によって、より柔軟に伝達特性の変化に対応することができる。一方、(ii)の方法により、リアルタイムの処理量を減らすことができる。使用状況に合わせて何れかの方法により、βを求め、保持すればよい。

次に、更新ベクトル算出部２０５は、（式１９）により、選定された更新ベクトルΔW(n)を用いて、特性ベクトルW(n)を更新する。

なお、（式１９）では、正則化行列として、（式１０）のΛ(=εI)に代えて、βA^-1を用いる。このような構成とすることで、特性ベクトルW(n)の各要素の分散の近似値の値が大きいほど、行列A^-1の各要素の大きさが小さくなるようにスケーリングしたものを正則化行列とすることができる。

このような構成により、更新ベクトル算出部２０５は、正則化行列をβA^-1として、（式１１）に代えて、（式A）により特性ベクトルの更新ベクトルΔWを算出し、（式１０）に代えて、（式１９）のように、後述するフィルタ係数保持部２０１の特性ベクトルW(n)を更新する。

＜フィルタ係数保持部２０１＞
フィルタ係数保持部２０１は、特性ベクトルW(n+1)を受け取り、特性ベクトルW(n+1)からなるフィルタ係数を保持する（Ｓ２０１）。

＜効果＞
本実施形態によれば、未知伝達系の出力信号の観測時に重畳される外乱信号の相関の影響も考慮された未知伝達系の特性ベクトルの推定を可能とする。そのため、音響エコーキャンセラにおけるダブルトークのように、音声信号が外乱信号に相当する場合や、低周波騒音等の周期性や相関の強い外乱信号が混入する場合においても、精度高く未知の伝達特性を推定することが可能となる。その結果、エコーやハウリングの発生を抑え、快適なハンズフリー通話の提供などにつながる。

＜変形例＞
入力信号と出力信号とを、複数の周波数成分の信号に分割し、分割された周波数成分毎の信号に対して、本実施形態を適用し、分割された周波数成分毎の伝達特性を推定してもよい。

＜第二実施形態＞
第一実施形態と異なる部分を中心に説明する。

伝達系推定装置３００は、フィルタ係数保持部２０１と、出力推定部２０２と、残差算出部２０３と、外乱分散共分散行列保持部２０４と、更新ベクトル算出部３０５とを含む（図１参照）。

＜更新ベクトル算出部３０５＞
更新ベクトル算出部３０５は、入力信号x(n)と残差信号の系列R(n)とを受け取り、行列A^-1を外乱分散共分散行列保持部２０４から取り出す。更新ベクトル算出部３０５は、第一実施形態の場合と同様に、行列A^-1と、入力信号の系列X(n)と、残差信号の系列R(n)とを用いて、想定残差信号の系列R^(n)が外乱信号の系列V(n)と一致する確率と、任意の更新候補ベクトルΔW^(n)の生起確率との積が最大となるよう選定された更新候補ベクトルを更新ベクトルΔW(n)とし、選定された更新ベクトルΔW(n)を用いて、特性ベクトルW(n)を更新し（図２のＳ３０５）、更新後の特性ベクトルW(n+1)を出力する。本実施形態では、後述のように、特性ベクトルW(n)の各要素の値が取りうる確率が、対角要素が互いに異なる値を取る対角行列に基づく分散共分散行列B^-1により表される多次元正規分布に従うものとして扱うことにより、第一実施形態とは異なる更新ベクトルΔW(n)が選定される。更新ベクトルΔW(n)の二つの選定方法について説明する。

(i)更新ベクトル算出部３０５は、特性ベクトルW(n)に対する更新ベクトルΔW(n-1)の各要素の大きさに対応する値1/β_lを対角要素とする対角行列B^-1を逐次算出する。例えば、

とする。ただし、（式２４）の更新ベクトルΔW(n-1)は、１ステップ前に得られる特性ベクトルW(n)に対する更新ベクトルΔW(n-1)、つまり、特性ベクトルW(n-1)を更新するために加算される更新ベクトルΔW(n-1)である。（式２３）では、行列B^-1の対角要素1/β_lを更新ベクトルΔW(n-1)の各要素Δw_l(n-1)の二乗の平均値として、近似的に計算している。なお、前述の通り、「特性ベクトルW(n)を確率変数とみなしたときの各要素の分散」は、「特性ベクトルの更新のために加算される更新ベクトルの各要素の大きさ」から近似することができるので、別の言い方をすると、更新ベクトル算出部３０５は、特性ベクトルW(n)を確率変数とみなしたときの各要素の分散に対応する値1/β_lを対角要素とする対角行列をB^-1として逐次算出する。また、例えば、直近のJ個の特性ベクトルW(n),W(n-1),…,W(n-J+1)のそれぞれのl番目の要素w_l(n),w_l(n-1),…,w_l(n-J+1)からl番目の要素の分散δ_lを求め、β_l=1/δ_lとしてもよい。

更新ベクトル算出部３０５は、算出した対角行列B^-1と、入力信号の系列X(n)と残差信号の系列R(n)と行列A^-1を用いて、次式により、更新ベクトルΔW(n)を算出する。

(ii)更新ベクトル算出部２０５は、あらかじめ固定値として対角行列B^-1を与えられる。例えば、更新ベクトル算出部３０５は、リアルタイムの伝達系推定に先立ち、予め未知伝達系１の特性ベクトルを算出しておき、第二特性ベクトルとする。更新ベクトル算出部３０５は、第二特性ベクトルの更新のために加算される更新ベクトルを求め、その各要素の大きさ（例えば二乗の平均値）に対応する値を対角要素とする対角行列B^-1を予め算出しておく。つまり、逐次得られる特性ベクトルW(n)に代えて、予め求めておいた第二特性ベクトルを用いて、B^-1を予め算出し、固定値として用いる。別の言い方をすると、更新ベクトル算出部３０５は、「特性ベクトルW(n)を確率変数とみなしたときの各要素の分散」（の近似値である「特性ベクトルの更新のために加算される更新ベクトルの各要素の大きさ」）に対応する値1/β_lを対角要素とする対角行列B^-1を予め算出しておく。

更新ベクトル算出部３０５は、あらかじめ固定値として与えられた対角行列B^-1と、入力信号の系列X(n)と残差信号の系列R(n)と行列A^-1を用いて、（式Ｂ）により、更新ベクトルΔW(n)を算出する。

上述の(i),(ii)の処理が、想定残差信号の系列R^(n)が外乱信号の系列V(n)と一致する確率と、任意の更新候補ベクトルΔW^(n)の生起確率との積が最大となる更新候補ベクトルを更新ベクトルΔW(n)として選定する処理に相当する。

なお、行列B^-1の対角要素1/β_lの逆数の平均、または平均値の逆数をとることで、第一実施形態のβを求めることができる。

次に、更新ベクトル算出部３０５は、（式２１）により、選定された更新ベクトルΔW(n)を用いて、特性ベクトルW(n)を更新する。

このような形式とすることで、外乱信号の相関の影響に加え、特性ベクトルの要素毎の分散（別の言い方をすると、更新ベクトルの要素毎の大きさ）の違いも考慮した更新が可能となる。ここでは、（式１６）に基づく更新式の導出において、行列B^-1を（式１８）で近似する代わりに、非負の要素を対角要素にもつ対角行列とする（（式２２）参照）。

上述の(i)において、対角行列B^-1は、特性ベクトルの更新のために加算された過去の更新ベクトルΔWの大きさ（別の言い方をすると、特性ベクトルの分散）から逐次見積られるといえる。また、上述の(ii)において、対角行列B^-1は、あらかじめ見込値として固定的に与えられるといえる。上述の(i)の方法によって、より柔軟に伝達特性の変化に対応することができる。一方、(ii)の方法により、リアルタイムの処理量を減らすことができる。使用状況に合わせて何れかの方法により、対角行列B^-1を求め、保持すればよい。

＜効果＞
このような構成とすることで、第一実施形態と同様の効果を得ることができる。さらに、本実施形態では、各要素に重み付けすることができ、より柔軟に特性ベクトルを更新することができる。

＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。

＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (6)

1. 推定対象の未知伝達系の伝達特性の推定値をフィルタ係数として持つ特性ベクトルと、前記未知伝達系の入力信号の系列とを用いて、前記未知伝達系の出力信号を推定し、その推定値である推定出力信号を得る出力推定部と、
   外乱信号が含まれ得る前記未知伝達系の出力信号である観測出力信号と前記推定出力信号との差分である残差信号を算出する残差算出部と、
   前記特性ベクトルに任意の更新候補ベクトルを加えたベクトルと前記入力信号の系列とを用いて生成しうる推定出力候補信号の系列を、前記観測出力信号の系列から差し引いた場合に得られる想定残差信号系列が前記外乱信号の系列と一致する確率と、前記任意の更新候補ベクトルの生起確率と、の積を最大とする更新候補ベクトルを更新ベクトルとして選定し、選定された前記更新ベクトルを用いて、前記特性ベクトルを更新する更新ベクトル算出部とを含み、
   前記想定残差信号系列が前記外乱信号の系列と一致する確率は、前記外乱信号の系列の分散共分散行列の非対角項が非零の要素を取りうる多次元正規分布に依存し、前記任意の更新候補ベクトルの生起確率は、前記任意の更新候補ベクトルの分散共分散行列の非対角項が全て零となる多次元正規分布に依存する、
   伝達系推定装置。
2. 請求項１の伝達系推定装置であって、
   n及びpを離散時間のインデックス、Pを1以上の整数の何れか、Kを2以上の整数の何れかとし、^Tは転置、E[]は期待値を表すものとし、前記外乱信号の系列の分散共分散行列A^-1は、
   (i)直近のKサンプルの前記残差信号の系列R(n)=[r(n),r(n-1),…,r(n-K+1)]^Tを用いて、
   

   

   
   により、逐次算出される、
   または、
   (ii)前記未知伝達系に前記入力信号の系列を与えていない状態において、予め観測された外乱信号の系列V(p) =[v(p),v(p-1),…,r(v-K+1)]^Tを用いて、
   

   

   
   により、予め算出される、
   伝達系推定装置。
3. 請求項１または２の伝達系推定装置であって、
   (i) 前記任意の更新候補ベクトルの分散共分散行列の対角項の各要素は、互いに等しく、前記特性ベクトルの各要素の分散を全要素について平均した値として逐次算出され、
   または、
   (ii) 前記任意の更新候補ベクトルの分散共分散行列の対角項の各要素は、互いに等しく、予め固定値として与えられ、
   離散時間のインデックスをn、特性ベクトルをW(n)、ステップサイズをμ、入力信号の系列をX(n)、前記任意の更新候補ベクトルの分散共分散行列の対角項の互いに等しい各要素の値を1/β、残差信号の系列をR(n)とし、前記更新ベクトル算出部は、
   

   

   
   により、前記特性ベクトルW(n)を更新する、
   伝達系推定装置。
4. 請求項１または２の伝達系推定装置であって、
   (i) 前記任意の更新候補ベクトルの分散共分散行列の対角項の各要素は、前記特性ベクトルの各要素の分散として逐次算出され、
   または、
   (ii) 前記任意の更新候補ベクトルの分散共分散行列の対角項の各要素は、予め固定値として与えられ、
   離散時間のインデックスをn、特性ベクトルをW(n)、ステップサイズをμ、入力信号の系列をX(n)、前記任意の更新候補ベクトルの分散共分散行列に相当する対角行列をB^-1、残差信号の系列をR(n)とし、前記更新ベクトル算出部は、
   

   

   
   により、前記特性ベクトルW(n)を更新する、
   伝達系推定装置。
5. 出力推定部が、推定対象の未知伝達系の伝達特性の推定値をフィルタ係数として持つ特性ベクトルと、前記未知伝達系の入力信号の系列とを用いて、前記未知伝達系の出力信号を推定し、その推定値である推定出力信号を得る出力推定ステップと、
   残差算出部が、外乱信号が含まれ得る前記未知伝達系の出力信号である観測出力信号と前記推定出力信号との差分である残差信号を算出する残差算出ステップと、
   更新ベクトル算出部が、前記特性ベクトルに任意の更新候補ベクトルを加えたベクトルと前記入力信号の系列とを用いて生成しうる推定出力候補信号の系列を、前記観測出力信号の系列から差し引いた場合に得られる想定残差信号系列が前記外乱信号の系列と一致する確率と、前記任意の更新候補ベクトルの生起確率と、の積を最大とする更新候補ベクトルを更新ベクトルとして選定し、選定された前記更新ベクトルを用いて、前記特性ベクトルを更新する更新ベクトル算出ステップとを含み、
   前記想定残差信号系列が前記外乱信号の系列と一致する確率は、前記外乱信号の系列の分散共分散行列の非対角項が非零の要素を取りうる多次元正規分布に依存し、前記任意の更新候補ベクトルの生起確率は、前記任意の更新候補ベクトルの分散共分散行列の非対角項が全て零となる多次元正規分布に依存する、
   伝達系推定方法。
6. 請求項１から請求項４の何れかに記載の伝達系推定装置として、コンピュータを機能させるためのプログラム。

Images