(57)【要約】【課題】より高速に伝達系の未知のパラメータを推定する。【解決手段】本発明の伝達系パラメータ推定装置は、出力推定部、誤差算出部、相関計算部、空間射影部、更新部を備える。出力推定部は、入力信号ベクトルＸ（ｎ）と推定パラメータベクトルＨ＾（ｎ）から、推定出力信号ｙ＾（ｎ）を求める。誤差算出部は、出力信号ｙ（ｎ）と推定出力信号ｙ＾（ｎ）の差分である誤差信号ｅ（ｎ）を求める。相関計算部は、１つ前の時間の相関ｇ＾（ｎ−１）と射影係数αｎ−１，βｎ−１から、相関ｇ＾（ｎ）を求める。空間射影部は、相関ｇ＾（ｎ）、推定パラメータベクトルＨ＾（ｎ）を用いて射影係数αｎ，βｎを計算し、射影ベクトルΔＨ＾（ｎ）を求める。更新部は、推定パラメータベクトルＨ＾（ｎ）と射影ベクトルΔＨ＾（ｎ）から推定パラメータベクトルＨ＾（ｎ＋１）を求める。【選択図】図２
2014099738.tif 000002
【特許請求の範囲】
【請求項１】
伝達系のモデルパラメータの推定値の集合を推定パラメータベクトルとして求める伝達系パラメータ推定装置であって、
ｎを離散的な時間を表す整数、Ｋを１以上の整数、ｋを離散的な時間を表す０以上Ｋ以下の整数、ｘ（ｎ）をｎ番目の時間での入力信号、ｙ（ｎ）をｎ番目の時間での出力信号、上付きのＴを転置を示す記号、Ｘ（ｎ）をｘ（ｎ），ｘ（ｎ−１），…，ｘ（ｎ−Ｋ）を要素とする入力信号ベクトル、Ｈ＾（ｎ）をＫ＋１個の要素を持つ推定パラメータベクトル、ｇ＾（ｎ）をｎ番目の時間での相関、αｎ，βｎをｎ番目の時間での射影係数、μを前記推定パラメータベクトルの更新量を調整するステップサイズ、ｇ＾（０）をあらかじめ定めた相関の初期値、α０，β０をあらかじめ定めた射影係数の初期値とし、
入力信号ベクトルＸ（ｎ）と推定パラメータベクトルＨ＾（ｎ）から、推定出力信号ｙ＾（ｎ）を求める出力推定部と、
出力信号ｙ（ｎ）と推定出力信号ｙ＾（ｎ）との差分である誤差信号ｅ（ｎ）を求める誤差算出部と、
ｇ＾（ｎ−１）とαｎ−１，βｎ−１から、相関ｇ＾（ｎ）を
【数１４】
2014099738.tif 000016
のように求める相関計算部と、
【数１５】
2014099738.tif 000017
のように射影ベクトルΔＨ＾（ｎ）を求める空間射影部と、
Ｈ＾（ｎ＋１）＝Ｈ＾（ｎ）＋μΔＨ＾（ｎ）
のように推定パラメータベクトルＨ＾（ｎ＋１）を求める更新部と、
を備える伝達系パラメータ推定装置。
【請求項２】
伝達系のモデルパラメータの推定値の集合を推定パラメータベクトルとして求める伝達系パラメータ推定方法であって、
ｎを離散的な時間を表す整数、Ｋを１以上の整数、ｋを離散的な時間を表す０以上Ｋ以下の整数、ｘ（ｎ）をｎ番目の時間での入力信号、ｙ（ｎ）をｎ番目の時間での出力信号、上付きのＴを転置を示す記号、Ｘ（ｎ）をｘ（ｎ），ｘ（ｎ−１），…，ｘ（ｎ−Ｋ）を要素とする入力信号ベクトル、Ｈ＾（ｎ）をＫ＋１個の要素を持つ推定パラメータベクトル、ｇ＾（ｎ）をｎ番目の時間での相関、αｎ，βｎをｎ番目の時間での射影係数、μを前記推定パラメータベクトルの更新量を調整するステップサイズ、ｇ＾（０）をあらかじめ定めた相関の初期値、α０，β０をあらかじめ定めた射影係数の初期値とし、
出力推定部が、入力信号ベクトルＸ（ｎ）と推定パラメータベクトルＨ＾（ｎ）から、推定出力信号ｙ＾（ｎ）を求める出力推定ステップと、
誤差算出部が、出力信号ｙ（ｎ）と推定出力信号ｙ＾（ｎ）との差分である誤差信号ｅ（ｎ）を求める誤差算出ステップと、
相関計算部が、ｇ＾（ｎ−１）とαｎ−１，βｎ−１から、相関ｇ＾（ｎ）を
【数１６】
2014099738.tif 000018
のように求める相関計算ステップと、
空間射影部が、
【数１７】
2014099738.tif 000019
のように射影ベクトルΔＨ＾（ｎ）を求める空間射影ステップと、
更新部が、
Ｈ＾（ｎ＋１）＝Ｈ＾（ｎ）＋μΔＨ＾（ｎ）
のように推定パラメータベクトルＨ＾（ｎ＋１）を求める更新ステップと、
を有する伝達系パラメータ推定方法。
【請求項３】
請求項１記載の伝達系パラメータ推定装置としてコンピュータを動作させるための伝達系パラメータ推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は電気信号、音、振動が伝搬する伝達系の未知のモデルパラメータを推定する伝達系パラメータ推定装置、伝達系パラメータ推定方法、伝達系パラメータ推定プログラムに関する。
【背景技術】
【０００２】
電気信号や音、振動が伝播する伝達系の未知の特性を推定することは、その伝達系への任意の入力に対する出力を予測する上で有益である。伝達系のパラメータを推定する問題は、有限のパラメータで表現可能なモデルをあらかじめ与え、パラメータの値を推定する問題として扱うことが多い。あらかじめ与えるモデルは、推定対象の伝達系の物理的素性（入力と出力の関係）に基づき選定される。例えば、入力と出力の関係が線形と見なせる場合、線形のＦＩＲ（Finite Impulse Response）フィルタやＩＩＲ（Infinite Impulse Response）フィルタが、離散時間領域における伝達系のモデルとして選ばれる。
【０００３】
本明細書では、未知の伝達系をＦＩＲフィルタでモデル化した場合を考える。このような伝達系の場合、伝達系への入力信号をｘ（ｎ）、出力信号をｙ（ｎ）とすると、
【０００４】
【数１】
2014099738.tif 000003
【０００５】
のように畳込演算の形で入力信号と出力信号との関係を表現でき、ｈ０，ｈ１，ｈ２，…，ｈＫが推定すべき未知のモデルパラメータ（ｋ個前の入力信号ｘ（ｎ−ｋ）に対するインパルス応答）である。また、ＨをＨ＝（ｈ０，ｈ１，…，ｈＫ）Ｔのインパルス応答ベクトルとする。ただし、次数Ｋは０以上の整数であり、ｎは離散時間を示す整数、上付きのＴは転置を示す記号である。
【０００６】
このようなパラメータ推定方法の従来技術として、非特許文献１に示された技術が知られている。図１は、非特許文献１に示された学習同定法またはＮＬＭＳ（Normalized Least Mean Squares）法に基づく、従来の伝達系パラメータ推定装置の機能構成例を示す図である。伝達系パラメータ推定装置９００は、出力推定部９１０、誤差算出部９２０、空間射影部９４０、更新部９５０で構成され、伝達系９９０のパラメータを推定する。
【０００７】
出力推定部９１０は、入力信号ｘ（ｎ）が入力されたときの伝達系９９０の出力信号ｙ（ｎ）の推定値である推定出力信号ｙ＾（ｎ）を
ｙ＾（ｎ）＝Ｘ（ｎ）ＴＨ＾（ｎ）      （２）
のように求める。ただし、Ｘ（ｎ）はＸ（ｎ）＝（ｘ（ｎ），ｘ（ｎ−１），…，ｘ（ｎ−Ｋ））Ｔの入力信号ベクトル、Ｈ＾（ｎ）はＨ＾（ｎ）＝（ｈ０＾（ｎ），ｈ１＾（ｎ），…，ｈＫ＾（ｎ））Ｔの推定パラメータベクトル、ｈｋ＾（ｎ）はｎ−１番目の時間までの情報で推定したｋ個前の入力信号ｘ（ｎ−ｋ）に対するインパルス応答である推定パラメータ、Ｋはインパルス応答長に相当する次数である。
【０００８】
誤差算出部９２０は、出力信号ｙ（ｎ）と推定出力信号ｙ＾（ｎ）から、誤差信号ｅ（ｎ）を、
ｅ（ｎ）＝ｙ（ｎ）−ｙ＾（ｎ）      （３）
のように求める。
【０００９】
空間射影部９４０は、インパルス応答ベクトルＨと推定パラメータベクトルＨ＾（ｎ）との差分ベクトルＨ−Ｈ＾（ｎ）について、入力信号ベクトルへの射影ベクトルΔＨ＾（ｎ）を、
ΔＨ＾（ｎ）＝Ｘ（ｎ）「Ｘ（ｎ）ＴＸ（ｎ）」−１Ｘ（ｎ）Ｔ［Ｈ−Ｈ＾（ｎ）］
＝Ｘ（ｎ）「Ｘ（ｎ）ＴＸ（ｎ）」−１［ｙ（ｎ）−ｙ＾（ｎ）］
＝Ｘ（ｎ）「Ｘ（ｎ）ＴＸ（ｎ）」−１ｅ（ｎ）                （４）
のように求める。ここで、Ｘ（ｎ）「Ｘ（ｎ）ＴＸ（ｎ）」−１Ｘ（ｎ）Ｔは、入力信号ベクトルＸ（ｎ）が構成する部分空間への射影行列を表している。
【００１０】
更新部９５０は、射影ベクトルΔＨ＾（ｎ）を更新ベクトルとして、推定パラメータベクトルＨ＾（ｎ＋１）を、
Ｈ＾（ｎ＋１）＝Ｈ＾（ｎ）＋μΔＨ＾（ｎ）      （５）
のように求める。ただし、μはステップサイズであり、０〜２の範囲で選定して与える。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】J. Nagumo, A. Noda, “A learning method for system identification”, IEEE Transactions on Automatic Control, Volume:12, Issue:3, Page(s):282-287, 1967.
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、従来技術においては、インパルス応答ベクトルＨと推定パラメータベクトルＨ＾（ｎ）との差分ベクトルＨ−Ｈ＾（ｎ）を、Ｘ（ｎ）「Ｘ（ｎ）ＴＸ（ｎ）」−１Ｘ（ｎ）Ｔによって非常に狭い部分空間に射影して、更新ベクトルΔＨ＾（ｎ）を得ていたため、推定の速度が遅くなるという課題がある。
【００１３】
本発明は、より高速な推定が可能な伝達系の未知のパラメータを推定する伝達系パラメータ推定装置を提供することを目的とする。
【課題を解決するための手段】
【００１４】
本発明の伝達系パラメータ推定装置は、伝達系のインパルス応答を推定した推定パラメータベクトルを求める。まず、ｎを離散的な時間を表す整数、Ｋを1以上の整数、ｋを離散的な時間を表す０以上Ｋ以下の整数、ｘ（ｎ）をｎ番目の時間での入力信号、ｙ（ｎ）をｎ番目の時間での出力信号、上付きのＴを転置を示す記号、Ｘ（ｎ）をｘ（ｎ），ｘ（ｎ−１），…，ｘ（ｎ−Ｋ）を要素とする入力信号ベクトル、ｈｋ＾（ｎ）をｎ−１番目の時間までの情報で推定したｋ個前の入力信号ｘ（ｎ−ｋ）に対するインパルス応答に相当する推定パラメータ、Ｈ＾（ｎ）をＫ＋１個の要素を持つ推定パラメータベクトル、ｇ＾（ｎ）をｎ番目の時間での相関、αｎ，βｎをｎ番目の時間での射影係数、μを更新量を調整するステップサイズ、ｇ＾（０）をあらかじめ定めた相関の初期値、α０，β０をあらかじめ定めた射影係数の初期値とする。
【００１５】
本発明の伝達系パラメータ推定装置は、出力推定部、誤差算出部、相関計算部、空間射影部、更新部を備える。出力推定部は、入力信号ベクトルＸ（ｎ）と推定パラメータベクトルＨ＾（ｎ）から、推定出力信号ｙ＾（ｎ）を求める。誤差算出部は、出力信号ｙ（ｎ）と推定出力信号ｙ＾（ｎ）の差分である誤差信号ｅ（ｎ）を求める。相関計算部は、１つ前の時間の相関ｇ＾（ｎ−１）と射影係数αｎ−１，βｎ−１から、相関ｇ＾（ｎ）を
【００１６】
【数２】
2014099738.tif 000004
【００１７】
のように求める。空間射影部は、
【００１８】
【数３】
2014099738.tif 000005
【００１９】
のように射影ベクトルΔＨ＾（ｎ）を求める。更新部は、
Ｈ＾（ｎ＋１）＝Ｈ＾（ｎ）＋μΔＨ＾（ｎ）
のように推定パラメータベクトルＨ＾（ｎ＋１）を求める。
【発明の効果】
【００２０】
本発明の伝達系パラメータ推定装置によれば、未知の真の伝達系パラメータのベクトルＨと推定パラメータベクトルＨ＾（ｎ）との間の内積Ｈ＾（ｎ）ＴＨである相関ｇ＾（ｎ）を再帰的に計算できるようにすることで、入力ベクトルと推定パラメータベクトルの２つのベクトルが構成する部分空間内で、推定パラメータベクトルの更新ベクトルを探索できる。したがって、従来の入力信号ベクトルだけで構成する部分空間よりも広い空間で更新ベクトルを探索できる。よって、本発明の伝達系パラメータ推定装置は、従来よりも高速に、未知の伝達系のパラメータを推定できる。
【図面の簡単な説明】
【００２１】
【図１】従来の伝達系パラメータ推定装置の機能構成例を示す図。
【図２】本発明の伝達系パラメータ推定装置の機能構成例を示す図。
【図３】本発明の伝達系パラメータ推定装置の処理フローの例を示す図。
【発明を実施するための形態】
【００２２】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【実施例１】
【００２３】
＜理論＞
まず、ｎを離散的な時間を表す整数、Ｋを1以上の整数、ｋを離散的な時間を表す０以上Ｋ以下の整数、ｘ（ｎ）をｎ番目の時間での入力信号、ｙ（ｎ）をｎ番目の時間での出力信号、上付きのＴを転置を示す記号、Ｘ（ｎ）をＸ（ｎ）＝（ｘ（ｎ），ｘ（ｎ−１），…，ｘ（ｎ−Ｋ））Ｔの入力信号ベクトル、ｈｋをｋ個前の入力信号ｘ（ｎ−ｋ）に対するインパルス応答、ＨをＨ＝（ｈ０，ｈ１，…，ｈＫ）Ｔのインパルス応答ベクトル、ｈｋ＾（ｎ）をｎ−１番目の時間までの情報で推定したｋ個前の入力信号ｘ（ｎ−ｋ）に対するインパルス応答（推定すべきパラメータ）、Ｈ＾（ｎ）をＨ＾（ｎ）＝（ｈ０＾（ｎ），ｈ１＾（ｎ），…，ｈＫ＾（ｎ））Ｔの推定パラメータベクトルとする。
【００２４】
本発明では、更新ベクトルΔＨ＾（ｎ）を、入力ベクトルと推定パラメータベクトルの２つのベクトルが構成する部分空間内で探索する。そのため、射影行列を
【００２５】
【数４】
2014099738.tif 000006
【００２６】
として、インパルス応答ベクトルＨと推定パラメータベクトルＨ＾（ｎ）との差分ベクトルＨ−Ｈ＾（ｎ）を射影する。すなわち、更新ベクトルΔＨ＾（ｎ）を
【００２７】
【数５】
2014099738.tif 000007
【００２８】
とする。式７の右辺の中の内積Ｈ＾（ｎ）ＴＨには、未知のパラメータであるインパルス応答ベクトルＨが含まれているので、このままではその値を得ることができない。そこで、射影係数ベクトル(αｎ，βｎ)Ｔを以下のようおく。
【００２９】
【数６】
2014099738.tif 000008
【００３０】
そうすると更新式は、式７，８より
【００３１】
【数７】
2014099738.tif 000009
【００３２】
のようになる。なお、μは０〜２の値を与えればよい。そして、式９より、
【００３３】
【数８】
2014099738.tif 000010
【００３４】
なる再帰式が得られる。
ここで、
【００３５】
【数９】
2014099738.tif 000011
【００３６】
を初期値とすれば、式８，１０を再帰的に計算することで、各離散時間ｎにおけるＨ＾（ｎ）ＴＨの値を計算できる。なお、伝達系のパラメータを更新する場合などでは、ｙ（０）／‖Ｘ（０）‖２≠０のときもある。
【００３７】
実際の計算では、インパルス応答ベクトルＨと推定パラメータベクトルＨ＾（ｎ）の内積Ｈ＾（ｎ）ＴＨは、１つのスカラー値として求められるので、内積Ｈ＾（ｎ）ＴＨを相関ｇ＾（ｎ）とする。そして、
【００３８】
【数１０】
2014099738.tif 000012
【００３９】
ｈ０＾（０）＝０，ｈ１＾（０）＝０，…，ｈＫ＾（０）＝０
を初期値として、
【００４０】
【数１１】
2014099738.tif 000013
【００４１】
のように射影ベクトルΔＨ＾（ｎ）を求め、
Ｈ＾（ｎ＋１）＝Ｈ＾（ｎ）＋μΔＨ＾（ｎ）      （１４）
のように推定パラメータベクトルを更新する計算（式１１〜１４）を繰り返せばよい。
【００４２】
＜具体的構成＞
図２に本発明の伝達系パラメータ推定装置の機能構成例を、図３に本発明の伝達系パラメータ推定装置の処理フローの例を示す。伝達系パラメータ推定装置１００は、伝達系９９０のインパルス応答を推定した推定パラメータベクトルＨ＾（ｎ＋１）を求める。伝達系パラメータ推定装置１００は、出力推定部１１０、誤差算出部１２０、相関計算部１３０、空間射影部１４０、更新部１５０を備える。伝達系パラメータ推定装置１００は、入力信号Ｘ（ｎ）が入力されるとステップＳ１１０に進む（Ｓ１０１）。
【００４３】
出力推定部１１０は、入力信号ベクトルＸ（ｎ）と推定パラメータベクトルＨ＾（ｎ）から、推定出力信号ｙ＾（ｎ）を求める（Ｓ１１０）。具体的には、
ｙ＾（ｎ）＝Ｘ（ｎ）ＴＨ＾（ｎ）      （１５）
のように求めればよい。なお、本実施例では、入力信号ベクトルＸ（ｎ）をＸ（ｎ）＝（ｘ（ｎ），ｘ（ｎ−１），…，ｘ（ｎ−Ｋ））Ｔ、推定パラメータベクトルＨ＾（ｎ）をＨ＾（ｎ）＝（ｈ０＾（ｎ），ｈ１＾（ｎ），…，ｈＫ＾（ｎ））Ｔとしているが、同じ結果が得られるようにベクトルの要素の順番を並び替えてもよい。例えば、Ｘ（ｎ）＝（ｘ（ｎ−Ｋ），ｘ（ｎ−Ｋ＋１），…，ｘ（ｎ））Ｔ、Ｈ＾（ｎ）＝（ｈＫ＾（ｎ），ｈＫ−１＾（ｎ），…，ｈ０＾（ｎ））Ｔとしてもよい。誤差算出部１２０は、出力信号ｙ（ｎ）と推定出力信号ｙ＾（ｎ）の差分である誤差信号ｅ（ｎ）を求める（Ｓ１２０）。具体的には、
ｅ（ｎ）＝ｙ（ｎ）−ｙ＾（ｎ）      （１６）
のように求めればよい。
【００４４】
相関計算部１３０は、１つ前の離散時間の相関ｇ＾（ｎ−１）と射影係数αｎ−１，βｎ−１から、相関ｇ＾（ｎ）を
【００４５】
【数１２】
2014099738.tif 000014
【００４６】
のように求める（Ｓ１３０）。空間射影部１４０は、
【００４７】
【数１３】
2014099738.tif 000015
【００４８】
のように射影ベクトルΔＨ＾（ｎ）を求める（Ｓ１４０）。更新部１５０は、
Ｈ＾（ｎ＋１）＝Ｈ＾（ｎ）＋μΔＨ＾（ｎ）      （２０）
のように推定パラメータベクトルＨ＾（ｎ＋１）を求め、出力推定部１１０の推定パラメータベクトルを更新する。そして、ｎに１を加算してステップＳ１０１に戻る。ステップＳ１０１では、入力信号があれば処理を続け、入力信号がなければ処理を終了する。
【００４９】
本発明の伝達系パラメータ推定装置によれば、未知の真の伝達系パラメータのベクトルＨと推定パラメータベクトルＨ＾（ｎ）との間の内積Ｈ＾（ｎ）ＴＨである相関ｇ＾（ｎ）を再帰的に計算できるようにすることで、入力ベクトルと推定パラメータベクトルの２つのベクトルが構成する部分空間内で、推定パラメータベクトルＨ＾（ｎ＋１）の更新ベクトルΔＨ＾（ｎ）を探索できる。したがって、従来の入力信号ベクトルだけで構成する部分空間よりも広い空間で更新ベクトルを探索できる。よって、本発明の伝達系パラメータ推定装置は、従来よりも高速に、未知の伝達系のパラメータを推定できる。
【００５０】
［プログラム、記録媒体］
上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００５１】
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
【００５２】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【００５３】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００５４】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００５５】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【符号の説明】
【００５６】
１００，９００  伝達系パラメータ推定装置
１１０，９１０  出力推定部
１２０，９２０  誤差算出部
１３０  相関計算部
１４０，９４０  空間射影部
１５０，９５０  更新部
９９０  伝達系
【図１】
2014099738.tif 000020
【図２】
2014099738.tif 000021
【図３】
2014099738.tif 000022