JP5524316B2 - パラメータ推定装置、エコー消去装置、パラメータ推定方法、及びプログラム - Google Patents

Description

本発明は、電気信号や、音、振動が伝播する伝達系が、非線形な特性を有する系（以下「非線形系」）と線形な特性を有する系（以下「線形系」）との縦続接続とみなせる場合に、その伝達系の未知のパラメータを推定する技術に関する。

電気信号や音、振動が伝播する未知の伝達系の特性を推定することは、その伝達系への任意の入力に対するその伝達系の出力を予測する上で有益である。伝達系を推定する問題は、その伝達系に対して有限のパラメータで表現可能なモデルを予め与え、そのパラメータの値を推定する問題として扱うことが多い。予め与えるモデルは、推定対象の伝達系の物理的素性（入力と出力の関係）に基づき選定される。例えば、入力と出力の関係が線形とみなせる場合、線形のFIR(Finite Impulse Response)フィルタやIIR(Infinite Impulse Response)フィルタが、離散時間領域における伝達系のモデルとして選ばれる。

ここで、図１に示すような伝達系１のパラメータを推定する問題がある。なお、伝達系１は入力と出力との関係が非線形な特性を有する。また、伝達系１において、非線形系２と線形系３とが分離可能で縦続接続型のモデルとして与えることができる。なおかつ、非線形系２の特性は過去の入力及び出力に依存しないハードクリップ特性とみなせ、かつ、線形系３の特性は有限長のインパルス応答とみなせる。

＜非線形系２のモデル＞
図１の非線形系２は、過去の入力及び出力に依存しないハードクリップ特性を有するとき、その入力信号x(n)と出力信号s(n)との関係は、閾値a(>0)を用いて以下のように表現されるハードクリップ関数fを用いて、

と与えられる。ここで、nは離散時間であり、出力信号s(n)は、入力信号x(n)に対しては、同じ離散時間nにおける値にのみ依存するとして関係づけられている。ここで、閾値aが推定すべき未知のパラメータである。以下、閾値aの推定値を非線形パラメータa^ともいう。

＜線形系３のモデル＞
図１の線形系３への入力信号を非線形系２の出力信号s(n)とすると、線形系３の出力信号y(n)を、

と畳込演算の形で表すことができ、ここで有限長のインパルス応答h₀,h₁,…,h_Lが推定すべき未知のパラメータである。以下、有限長のインパルス応答h₀,h₁,…,h_Lの推定値を線形パラメータh^₀,h^₁,…,h^_Lともいう。ただし、次数Lは0以上の整数であり、出力信号s(n)は、実際には直接観測できない内部信号である。

以上より、伝達系１において推定の対象とするパラメータは、式(1)におけるaと式(2)におけるh₀,h₁,…,h_Lである。

＜従来技術に係るパラメータ推定装置６＞
非特許文献１には、このような枠組におけるパラメータ推定方法が示されている。図２は、非特許文献１に基づく従来のパラメータ推定方法を実現する構成例である。以下、パラメータ推定装置６について説明する。

なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号「^」等は、本来直前の文字の真上に記載されるべきものであるが、テキスト記法の制限により、その文字の直後に記載する。式中においてはこれらの記号は本来の位置に記述している。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。

離散時間nにおいて、パラメータ推定装置６の各部において以下の処理を行う。

非線形パラメータ適用部６０１では、時間領域の入力信号x(n)を受け取り、入力信号x(n)と非線形パラメータa^(n)とを用いて、次式により、時間領域の入力信号x(n)に対する非線形系２の出力信号の推定値である非線形系出力推定値s^(n)を求め、出力する。

信号蓄積部６０２は、非線形系出力推定値s^(n)と後述する偏微分値s^'(n)とを受け取り、Lサンプル過去まで遡って、L+1個の非線形系出力推定値s^(n),s^(n-1),…,s^(n-L)と、L+1個の偏微分値s^'(n),s^'(n-1),…,s^'(n-L)とを蓄積する。なお、S^(n)=[s^(n),s^(n-1),…,s^(n-L)]^T、S^'(n)=[s^'(n),s^'(n-1),…,s^'(n-L)]^Tとする。ここで^Tは転置を表す。

線形パラメータ乗算部６０３は、信号蓄積部６０２からL+1個の非線形系出力推定値s^(n),s^(n-1),…,s^(n-L)を取り出し、次式のように、L+1個の非線形系出力推定値s^(n),s^(n-1),…,s^(n-L)とL+1個の線形パラメータh^₀(n),h^₁(n),…,h^_L(n)をそれぞれ乗じた結果を全て加算し、推定対象の伝達系からの出力信号の推定値である伝達系出力推定値y^(n)を求め、出力する。
y^(n)=h^₀(n)s^(n)+h^₁(n)s^(n-1)+…+h^_L(n)s^(n-L) (4)
を出力する。なお、H^(n)=[h^₀(n),h^₁(n),…,h^_L(n)]^Tと標記すると、式(4)は、
y^(n)=H^(n)^TS^(n) (5)
となり、ベクトルの内積演算に相当する。

誤差算出部６０４では、推定対象の伝達系からの出力信号y(n)と伝達系出力推定値y^(n)とを受け取り、次式のように、その差を求め、誤差信号e(n)として出力する。
e(n)=y(n)-y^(n) (6)

微分信号生成部６０５は、次式のように非線形パラメータa^(n)についての非線形系出力推定値s^(n)の偏微分値s^'(n)を求め、出力する。なお、次式において、入力信号x(n)に代えて、非線形系出力推定値s^(n)を用いてもよい。

パラメータ更新部６０６は、誤差信号e(n)を受け取り、信号蓄積部６０２から偏微分値S^'(n)と非線形系出力推定値S^(n)とを取り出す。

ここで、誤差信号e(n)について、

と与えられる評価式について、非線形パラメータa^(n)及び線形パラメータH^(n)についての勾配が、

と得られる。
そのため、パラメータ更新部６０６は、次式により、非線形パラメータa^(n)及び線形パラメータH^(n)を更新し、時刻nにおける新たな更新値を出力する。

なお、ステップサイズであるμ_a、μ_Hは、それぞれ0≦μ_a≦1、0≦μ_H≦1の範囲の中で選定される。

A. Stenger, W. Kellermann, "Nonlinear acoustic echo cancellation with fast converging memoryless preprocessor", IEEE International Conference on Acoustics, Speech, and Signal Processing, 2000. ICASSP '00, 2000, Vol. 2, II805 - II808.

図１に示すような非線形系２と線形系３とが縦続接続された伝達系１を式(1)及び式(2)でモデル化し、それらの非線形パラメータa^(n)及び線形パラメータH^(n)を、それぞれ式(11)及び式(12)に従って逐次更新する従来技術では、伝達系出力推定値y^(n)の計算において、線形パラメータH^(n)を用いた式(4)または式(5)に相当する時間領域での畳込演算の実行に多くの演算量を要する。演算量の大きい時間領域での畳込演算は、周波数領域で効率的に実行可能であることが知られている。一方、非線形パラメータa^(n)を用いた式(3)に相当する非線形系出力推定値s^(n)の計算は、時間領域で計算した方が効率的である。従って、非線形系出力推定値s^(n)の計算を時間領域で実行した後、伝達系出力推定値y^(n)の計算を周波数領域で実行すればよい。しかし、この場合、非線形パラメータa^(n)及び線形パラメータH^(n)の逐次更新において、全てを時間領域で実行することを想定した式(11)、式(12)をそのまま適用することはできない。

本発明は、非線形系２の出力信号s(n)及び線形系３の出力信号y(n)を推定する際の演算量を小さく抑えながらも、非線形パラメータ及び線形パラメータを統一的に更新できるパラメータ推定技術を提供することを目的とする。

上記の課題を解決するために、本発明の第一の態様によれば、パラメータ推定装置は過去の入力及び出力に依存しないハードクリップ特性とみなせる非線形系と有限長のインパルス応答とみなせる線形系との縦続接続とみなせる伝達系の、ハードクリップ特性の閾値と有限長のインパルス応答とを推定する。パラメータ推定装置は、閾値の推定値である非線形パラメータを用いて、時間領域の伝達系の入力信号に対する非線形系の出力信号の推定値である非線形系出力推定値を求める非線形パラメータ適用部と、非線形系出力推定値の非線形パラメータについての偏微分値に相当する値を求める微分信号生成部と、時間領域の非線形系出力推定値及び偏微分値に相当する値をそれぞれ周波数領域に変換し、周波数領域の非線形系出力推定値及び偏微分値に相当する値を求める周波数領域変換部と、周波数領域の有限長のインパルス応答の推定値である線形パラメータと、周波数領域の非線形系出力推定値とを乗じ、伝達系の出力信号の周波数領域の推定値である伝達系出力推定値を求める周波数領域線形パラメータ乗算部と、伝達系の出力信号と、伝達系出力推定値との差を求める誤差算出部と、周波数領域の差と偏微分値に相当する値と非線形系出力推定値とを用いて、非線形パラメータを更新し、線形パラメータを周波数領域で更新する混合型パラメータ更新部と、を含む。

上記の課題を解決するために、本発明の第二の態様によれば、パラメータ推定方法は、過去の入力及び出力に依存しないハードクリップ特性とみなせる非線形系と有限長のインパルス応答とみなせる線形系との縦続接続とみなせる伝達系の、ハードクリップ特性の閾値と有限長のインパルス応答とを推定する。パラメータ推定方法は、閾値の推定値である非線形パラメータを用いて、時間領域の伝達系の入力信号に対する非線形系の出力信号の推定値である非線形系出力推定値を求める非線形パラメータ適用ステップと、非線形系出力推定値の非線形パラメータについての偏微分値に相当する値を求める微分信号生成ステップと、時間領域の非線形系出力推定値及び偏微分値に相当する値をそれぞれ周波数領域に変換し、周波数領域の非線形系出力推定値及び偏微分値に相当する値を求める周波数領域変換ステップと、周波数領域の有限長のインパルス応答の推定値である線形パラメータと、周波数領域の非線形系出力推定値とを乗じ、伝達系の出力信号の周波数領域の推定値である伝達系出力推定値を求める周波数領域線形パラメータ乗算ステップと、伝達系の出力信号と、伝達系出力推定値との差を求める誤差算出ステップと、周波数領域の差と偏微分値に相当する値と非線形系出力推定値とを用いて、非線形パラメータを更新し、線形パラメータを周波数領域で更新する混合型パラメータ更新ステップと、を含む。

本発明によれば、非線形系２の出力信号s(n)を時間領域で推定し、線形系３の出力信号y(n)を周波数領域で推定することで演算効率を高め、なおかつ、非線形パラメータ及び線形パラメータを統一的に更新できるという効果を奏する。

伝達系を説明するための図。 従来技術を用いたパラメータ推定装置の機能ブロック図。 第一実施形態に係るパラメータ推定装置の機能ブロック図。 第一実施形態に係るパラメータ推定装置の処理フローを示す図。 第一実施形態の変形例に係るパラメータ推定装置の機能ブロック図。 第一実施形態の変形例に係るパラメータ推定装置の処理フローを示す図。 第一実施形態のパラメータ推定装置をエコー消去装置に転用した場合の機能ブロック図。 第一実施形態の変形例のパラメータ推定装置をエコー消去装置に転用した場合の機能ブロック図。

以下、本発明の実施形態について説明する。

＜第一実施形態のポイント＞
本実施形態では、時間領域において演算効率の高い非線形系における処理と、周波数領域において演算効率の高い線形系における処理とを効率的に融合させるために、非線形系の推定パラメータを時間領域での形式、線形系の推定パラメータを周波数領域での形式として保持しておきながら、両者を統一的に更新可能な更新式を導出した。その方法について説明する。

時間領域における推定パラメータの更新式である式(11)、式(12)を次式のように一般化する。

ここで、非線形系２の出力信号s(n)の推定計算を時間領域、線形系３の出力信号y(n)の推定計算を周波数領域で行うために適した、非線形パラメータ、線形パラメータの更新式の具体的な形式を考える。まず、式(13)、式(14)について、ステップサイズμ_a＝μ_H＝１として、更新値Δa^(n)、ΔH^(n)について、対応する時刻nについての入出力関係を満足するものとすると、次式の関係が得られる。

これに対し、N個の誤差信号e(n)からなる誤差信号ベクトル
E(n)=[e(n),e(n-1),…,e(n-N+1)] (16)
を、離散フーリエ変換、離散コサイン変換、フィルタバンク等を用いて、周波数領域に変換した場合、式(15)に対応する関係は、次式のようになる（少なくとも、近似的に次式の関係が成り立つとみなせる）。

但し、

ここで、E(ω_k,i)、H^(ω_k,i)、S^(ω_k,i),S^'(ω_k,i)はk番目の離散周波数ω_kに対応するE(n),H^(n),S^(n),S^’(n)の周波数成分を表し（ただし、k=0,1,…,N-1）、iはブロック番号、式(16)のように、信号をN個の要素によるブロック単位としてまとめる時間間隔（例えば、シフト幅であり、1以上N以下の離散時間に対応する時間間隔）毎に割り当てられる番号であり、Δa^(i)はブロックiにおける非線形パラメータa^(i)の時間領域の更新量に相当し、ΔH^(ω_k,i)はブロックiにおける線形パラメータH^(ω_k,i)の周波数領域における更新量に相当する。

さて、式(17)は未知数よりも、式の数の少ない不定方程式であり、解の与え方には自由度が存在する。本実施形態では、式(17)の重み付き最小ノルム解を各パラメータの更新量として採用する。そのために、式(17)に重み行列Wを導入し、

とする。ここで

である。式(18)について、次式により、重み付き最小ノルム解が得られる。但し、^＊は共役転置を表す。

しかしながら、式(19)について、

の逆行列計算をそのまま実行することは多大な演算量を要する。そこで、本実施形態では、
「逆行列の補助定理」
A=B^-1+CD^-1C^*
A^-1=B-BC(D+C^*BC)^-1C^*B (20)
を利用し、上記の逆行列演算を効率化させる。
まず、

と変形した上で、式(20)において、

ここで、S^S^^*は対角行列であるため、対角要素の逆数を取ることで、簡単に逆行列演算が可能である。式(21)を式(19)に代入すると、次式のように、更新量が得られる。

ここで、αはスカラー量である。これより、非線形パラメータの更新式は、

となる。なお、real(A)は、複素数Aの実部のみを返す関数である。非線形パラメータa^(i)は実数値であるため、複素数のαの実部のみを更新に用いるようにしている。また、線形パラメータの更新式は、離散周波数ω_k毎に独立に、

となる。ここで、ステップサイズμ_hは、0≦μ_h≦１であり、離散周波数ω_k毎に異なる値を与えてもよいし、時間変化するように与えてもよい。なお、式(24)に代えて、次式により、線形パラメータを更新してもよい。

なお、δhは、零除算防止のための正の定数である。以下、パラメータ推定装置の構成例を示す。

＜第一実施形態に係るパラメータ推定装置＞
図３は第一実施形態に係るパラメータ推定装置１０の機能ブロック図、図４はその処理フローを示す図である。

パラメータ推定装置１０は、非線形パラメータ適用部１１０、微分信号生成部１２０、信号蓄積部１３０、周波数領域変換部１４０、線形パラメータ乗算部１５０、周波数領域変換部１６０、誤差算出部１７０及び混合型パラメータ更新部１８０を含む。

パラメータ推定装置１０は、伝達系１の入力信号x(n)と出力信号y(n)とを受け取り、非線形系２のハードクリップ特性の閾値aと、線形系３の有限長のインパルス応答h₀,h₁,…,h_Lとを推定し、閾値aの推定値である非線形パラメータa^(i)と、有限長のインパルス応答h₀,h₁,…,h_Lの周波数領域の推定値である線形パラメータH^(ω_k,i)とを出力する。以下、各部の処理内容を説明する。

離散時間nにおいて、パラメータ推定装置１０の非線形パラメータ適用部１１０、微分信号生成部１２０及び信号蓄積部１３０において以下の処理を行う。

＜非線形パラメータ適用部１１０＞
非線形パラメータ適用部１１０では、非線形パラメータa^(i)と時間領域の入力信号x(n)とを受け取り、これらの値を用いて、次式またはそれと等価な式あるいは近似式により、時間領域の入力信号x(n)に対する非線形系２の出力信号の推定値である非線形系出力推定値s^(n)を求め（ｓ１）、微分信号生成部１２０及び信号蓄積部１３０に出力する。

なお、a^(i)は、ブロックi毎（例えばシフト幅毎）に混合型パラメータ更新部１８０から受け取る値である。

＜微分信号生成部１２０＞
微分信号生成部１２０は、例えば次式またはそれと等価な式あるいは近似式により、非線形系出力推定値s^(n)の非線形パラメータa^(i)についての偏微分値に相当する値s^'(n)を求め（ｓ３）、信号蓄積部１３０に出力する。なお、図３に示すように、非線形系出力推定値s^(n)を受け取り、次式において、入力信号x(n)に代えて、非線形系出力推定値s^(n)を用いても、非線形系出力推定値s^(n)を非線形パラメータa^(i)について偏微分し、値s^'(n)を求めることができる。

＜信号蓄積部１３０＞
信号蓄積部１３０は、非線形系出力推定値s^(n)と値s^'(n)とを受け取り、Lサンプル過去まで遡って、L+1個の非線形系出力推定値s^(n),s^(n-1),…,s^(n-L)と、L+1個の値s^'(n),s^'(n-1),…,s^'(n-L)とを蓄積する（ｓ５）。

ブロックi毎（例えばシフト幅毎）に、パラメータ推定装置１０の周波数領域変換部１４０及び１６０、線形パラメータ乗算部１５０、誤差算出部１７０及び混合型パラメータ更新部１８０において以下の処理を行う。

＜周波数領域変換部１４０及び１６０＞
周波数領域変換部１４０は、信号蓄積部１３０から非線形系出力推定値S^(n)=[s^(n),s^(n-1),…,s^(n-L)]と、値S^'(n)=[s^'(n),s^'(n-1),…,s^'(n-L)]とを取り出し、それぞれ周波数領域に変換し、周波数領域の非線形系出力推定値S^(ω_k,i)と、値S^'(ω_k,i)とを求め（ｓ７）、非線形系出力推定値S^(ω_k,i)及び値S^'(ω_k,i)を混合型パラメータ更新部１８０に、非線形系出力推定値S^(ω_k,i)を線形パラメータ乗算部１５０に出力する。なお、周波数領域への変換方法としては、離散フーリエ変換、離散コサイン変換、フィルタバンク等が考えられる。

同様に周波数領域変換部１６０は、伝達系１の出力信号y(n)を受け取り、図示しない信号蓄積部にL+1個の伝達系１の出力信号y(n),y(n-1),…,y(n-L)を蓄積しておき、ブロックi毎に、L+1個の伝達系１の出力信号y(n),y(n-1),…,y(n-L)を周波数領域に変換し、周波数領域の伝達系１の出力信号Y(ω_k,i)を求め（ｓ１１）、誤差算出部１７０に出力する。

＜線形パラメータ乗算部１５０＞
線形パラメータ乗算部１５０は、線形パラメータH^(ω_k,i)と周波数領域の非線形系出力推定値S^(ω_k,i)とを受け取り、次式のように、離散周波数ω_k毎にこれらの値を乗じ、周波数領域の伝達系１の出力信号Y(ω_k,i)の推定値である伝達系出力推定値Y^(ω_k,i)を求め（ｓ９）、誤差算出部１７０に出力する。
Y^(ω_k,i)=H^(ω_k,i)S^(ω_k,i)
なお、この式の計算量が、式(4)、式(5)の計算量に比べて少ないため、パラメータ推定装置１０は従来技術に比べ、計算量を少なくすることができる。このような構成により、非線形系２の出力信号s(n)を時間領域で推定し、線形系３の出力信号y(n)を周波数領域で推定することができる。

＜誤差算出部１７０＞
誤差算出部１７０は、周波数領域の伝達系１の出力信号Y(ω_k,i)と伝達系出力推定値Y^(ω_k,i)とを受け取り、次式のように、離散周波数ω_k毎にその差を誤差信号E(ω_k,i)として求め（ｓ１３）、混合型パラメータ更新部１８０に出力する。
E(ω_k,i)=Y(ω_k,i)-Y^(ω_k,i)

＜混合型パラメータ更新部１８０＞
混合型パラメータ更新部１８０は、周波数領域の誤差信号E(ω_k,i)と値S^'(ω_k,i)と非線形系出力推定値S^(ω_k,i)とを受け取り、これらの値を用いて、非線形パラメータa^(i)を更新し、線形パラメータH^(ω_k,i)を周波数領域で更新し（ｓ１５）、更新後の非線形パラメータa^(i+1)を非線形パラメータ適用部１１０に、更新後の線形パラメータH^(ω_k,i+1)を線形パラメータ乗算部１５０に出力する。

具体的には、式(23)、(24)により、更新する。

と零除算防止のための正の定数δ（大きくてもH^S^'^*[S^S^^*]^-1H^S^'の平均値を超えない程度）を分母に与えてもよい。このような更新式を用いることで、非線形パラメータを時間領域での形式、線形系の推定パラメータを周波数領域での形式として保持しておきながら、両者を統一的に更新できる。また、重み係数wは正数であり、非線形パラメータの推定に重点を置く場合は、１より小さい値を与え、線形パラメータの推定に重点を置く場合は、１より大きい値を与える。利用者は、非線形系及び線形系の推定パラメータの更新重みを自由に制御することができる。

ここで、重み係数wについて、

と時間変化するように与えてもよい。ここで、Nは周波数分割数である。

また、各離散周波数ω_kについて、線形パラメータH^(ω_k,i)のスカラー値の代わりに、長さM(M個の係数を持つ)の(小さな)FIRフィルタで構成してもよい。このとき、重み係数wは、

と与えることもできる。線形パラメータH^(ω_k,i)のノルムが１前後の値となるのに対し、非線形パラメータa^(i)の大きさは振幅に依存し、その取りうる範囲は線形パラメータH^(ω_k,i)に比べ非常に大きいものとなる。非線形パラメータa^(i)の値に応じて、重み係数wを時間変化するように与えることで、重点的に推定するパラメータをより適切に設定することができる。つまり、非線形パラメータa^(i)の値が大きければ重み係数wの値を小さくして非線形パラメータの推定に重点を置き、小さければ重み係数wの値を大きくして線形パラメータの推定に重点を置く。このような構成により、演算的に効率的かつ精度の高い推定処理を可能としている。

なお、パラメータ推定装置１０は、ブロックi毎に、非線形パラメータa^(i+1)及び線形パラメータH^(ω_k,i+1)を推定値として出力してもよいし、非線形パラメータa^(i+1)及び線形パラメータH^(ω_k,i+1)の値が収束した後に、収束したと判断したときの非線形パラメータa^(i+1)及び線形パラメータH^(ω_k,i+1)を推定値として出力してもよい。例えば、以下のような場合に収束したと判断する。

（１）非線形パラメータa^(i)及び線形パラメータH^(ω_k,i)についての各更新前後の差が一定値以下になった場合（二つの差とも一定値以下になるまで二つのパラメータを更新してもよいし、何れか一方の差が一定値以下になった場合には、該当するパラメータのみ更新を停止してもよい）
（２）上記処理ｓ１〜ｓ１５を一定数回（例えば100回）繰り返した場合
（３）（１）または（２）の何れかを満たす場合

また、パラメータ推定装置１０は、線形パラメータH^(ω_k,i+1)を時間領域に変換して出力してもよい。

＜効果＞
本実施形態によるパラメータ推定装置によれば、非線形系出力推定値を時間領域において計算し、線形系出力推定値を周波数領域において計算することで、演算効率を高めることができる。また、式(22)〜(24)により、非線形パラメータ及び線形パラメータを統一的に更新できる。なおかつ、周波数領域における誤差信号、非線形系出力推定値及びその偏微分値に相当する値を統一的に取り扱うことで、非線形系及び線形系の推定パラメータの更新重みを自由に制御可能な更新式を導出することにより、演算的に効率的かつ精度の高い推定処理を可能としている。

＜変形例＞
図５は第一実施形態の変形例の機能ブロック図、図６はその処理フローを示す図である。第一実施形態と異なる部分についてのみ説明する。

パラメータ推定装置２０は、非線形パラメータ適用部１１０、微分信号生成部１２０、信号蓄積部１３０、周波数領域変換部１４０、線形パラメータ乗算部１５０、誤差算出部２７０、時間領域変換部１９１、周波数領域変換部１９３及び混合型パラメータ更新部１８０を含む。つまり、時間領域変換部１９１及び周波数領域変換部１９３を含む点、周波数領域変換部１６０を含まない点、誤差算出部２７０の入出力が第一実施形態と異なる。ブロックi毎に、誤差算出部２７０、時間領域変換部１９１、周波数領域変換部１９３において以下の処理を行う。

＜時間領域変換部１９１＞
時間領域変換部１９１は、線形パラメータ乗算部１５０から伝達系出力推定値Y^(ω_k,i)を受け取り、時間領域に変換し、時間領域の伝達系出力推定値Y^(n)=[y^(n),y^(n-1),…,y^(n-L)]を求め（ｓ２１）、誤差算出部２７０に出力する。なお、時間領域への変換方法としては、周波数領域変換部１４０の変換方法に対応するものを用いればよい。

＜誤差算出部２７０＞
誤差算出部２７０は、伝達系１の出力信号y(n)と伝達系出力推定値Y^(n)=[y^(n),y^(n-1),…,y^(n-L)]とを受け取り、次式のように、離散時間n毎にその差を誤差信号e(n)として求め（ｓ２３）、周波数領域変換部１９３に出力する。
e(n)=y(n)-y^(n)
例えば、誤差算出部２７０は、図示しない信号蓄積部にL+1個の伝達系１の出力信号y(n),y(n-1),…,y(n-L)を蓄積しておき、伝達系出力推定値Y^(n)を受け取る毎に、L+1個の誤差信号e(n),e(n-1),…,e(n-L)からなる誤差信号ベクトルE(n)を求める。また、シフト幅分だけ伝達系１の出力信号y(n)を蓄積しておき、伝達系出力推定値Y^(n)を受け取る毎に、シフト幅分だけ誤差信号e(n)を求め、それ以前に求めた誤差信号と合わせて1ブロック分の誤差信号e(n),e(n-1),…,e(n-L)からなる誤差信号ベクトルE(n)を周波数領域変換部１９３に出力する構成としてもよい。

＜周波数領域変換部１９３＞
周波数領域変換部１９３は、誤差信号ベクトルE(n)を受け取り、周波数領域に変換し、周波数領域の誤差信号E(ω_k,i)を求め（ｓ２５）、混合型パラメータ更新部１８０に出力する。

このような構成により、第一実施形態と同様の効果を得ることができる。

＜第二実施形態＞
第一実施形態及びその変形例のパラメータ推定装置は、例えばエコー消去装置への転用が可能である。図７及び図８はそれぞれ第一実施形態及びその変形例のパラメータ推定装置をエコー消去装置に転用した場合の機能ブロック図である。

エコー消去装置３０及び４０は、スピーカ２０２から再生される音響信号が室内で反響し、同一室内にあるマイクロホン３０２へ回り込むエコーを推定して得られる擬似エコー信号をマイクロホンの収音信号から差し引くことで、エコーの消去された信号を出力するものである。

図７及び図８に示すように、スピーカ２０２で再生するための再生信号x(n)に増幅器２０１でゲインをかけてスピーカ２０２から再生する際、増幅器２０１やスピーカ２０２により信号に歪が生じる場合がある。すなわち、スピーカ２０２で再生するための再生信号x(n)が、歪を発生させる増幅器２０１やスピーカ２０２を経由し、スピーカ２０２から放音される信号が得られるまでの伝達特性は非線形とみなし、非線形系２を構成しているものと想定できる。言い換えると、非線形系２はスピーカ２０２及び増幅器２０１を含む。また、スピーカ２０２から放出された音響信号が、室内の反響路３０１を経由して、マイクロホン３０２に回り込み、収音信号y(n)が得られるまでの伝達特性は線形とみなし、線形系３を構成するものと想定できる。言い換えると、線形系３は反響路３０１及びマイクロホン３０２を含む。

この場合、図３あるいは図５のパラメータ推定装置１０あるいは２０は、それぞれ図７あるいは図８に示すエコー消去装置３０あるいは４０として動作することが可能であり、上述の説明にて入力信号x(n)をスピーカ２０２で再生するための再生信号x(n)に、出力信号y(n)をスピーカ２０２と同一の空間でマイクロホン３０２で収音した収音信号y(n)に読み替え（置き換え）ればよいので、重複説明を省略する。なお、エコー消去装置３０または４０は、非線形パラメータa^(i+1)も線形パラメータH^(ω_k,i)も外部に出力する必要はない代わりに、誤差信号e(n)を送信信号として出力する。このため、図７のエコー消去装置３０の構成においては、E(ω_k,i)をe(n)に変換する時間領域変換部３１０が新たに設けられている。

収音信号y(n)には、伝達系に由来するエコーのほかに、近端の話者の発話音声などが含まれる。近端の話者の音声とエコーが混合している場合には、y(n)から擬似エコー信号である伝達系出力推定値y^(n)を差し引いて得られる誤差信号e(n)は、近端の話者の音声が主となり、通信相手への送信信号となる。しかしながら同時に、誤差信号e(n)の中に存在する近端話者の音声は、非線形パラメータや線形パラメータの推定結果を乱す要因ともなる。そこで、以下のような対処が別途必要である。ステップサイズμ_aまたはμ_hの値を固定値とする場合は、その値を十分に小さく与える。近端話者の音声の存在を検知できる場合は、近端話者の音声が検知された場合にのみ、ステップサイズμ_aまたはμ_hを小さくするように制御する（例えば、μ_a=μ_h=0とする）。

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

＜プログラム及び記録媒体＞
上述したパラメータ推定装置やエコー消去装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置（各種実施形態で図に示した機能構成をもつ装置）として機能させるためのプログラム、またはその処理手順（各実施形態で示したもの）の各過程をコンピュータに実行させるためのプログラムを、ＣＤ−ＲＯＭ、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。

Claims (9)

1. 過去の入力及び出力に依存しないハードクリップ特性とみなせる非線形系と有限長のインパルス応答とみなせる線形系との縦続接続とみなせる伝達系において、前記伝達系の入力信号と出力信号とを用いて、前記ハードクリップ特性の閾値と前記有限長のインパルス応答とを推定するパラメータ推定装置であって、
   前記閾値の推定値である非線形パラメータを用いて、前記伝達系の時間領域の入力信号に対する前記非線形系の出力信号の推定値である時間領域非線形系出力推定値を求める非線形パラメータ適用部と、
   前記時間領域非線形系出力推定値の前記非線形パラメータに関する偏微分値に相当する値を求める微分信号生成部と、
   前記時間領域非線形系出力推定値の系列及び前記偏微分値に相当する値の系列をそれぞれ周波数領域に変換し、周波数領域の周波数領域非線形系出力推定値及び周波数領域の偏微分値に相当する値である周波数領域偏微分値相当値を求める周波数領域変換部と、
   前記有限長のインパルス応答の周波数領域の推定値である線形パラメータと、前記周波数領域非線形系出力推定値とを乗じ、前記伝達系の出力信号の周波数領域の推定値である周波数領域伝達系出力推定値を求める周波数領域線形パラメータ乗算部と、
   前記伝達系の周波数領域の出力信号である周波数領域伝達系出力信号と、前記周波数領域伝達系出力推定値との差である周波数領域差を求める誤差算出部と、
   前記非線形パラメータと前記線形パラメータのそれぞれを、前記周波数領域差と前記周波数領域偏微分値相当値と前記周波数領域非線形系出力推定値とを用いて、周波数領域で更新する混合型パラメータ更新部と、を含む、
   パラメータ推定装置。
2. 請求項１記載のパラメータ推定装置であって、
   iを時間領域の信号をN個の要素によるブロック単位としてまとめる時間間隔毎に割り当てられる番号とし、μ_a及びμ_hをそれぞれ0≦μ_a≦1及び0≦μ_h≦1の範囲の中で選定されるステップサイズとし、real(A)を複素数Aの実部のみを返す関数とし、k=0,1,…,N-1とし、ω_kをk番目の離散周波数とし、a^(i)及びH^(ω_k,i)をそれぞれ前記非線形パラメータ及び前記線形パラメータとし、S^(ω_k,i)、S^'(ω_k,i)及びE(ω_k,i)をそれぞれ前記周波数領域非線形系出力推定値、前記周波数領域偏微分値相当値及び前記周波数領域差とし、^*を共役転置とし、
   重み係数wを正数とし、前記混合型パラメータ更新部において、
   

   

   
   ただし、
   

   

   
   により、前記非線形パラメータa^及び前記線形パラメータH^を更新する、
   パラメータ推定装置。
3. 請求項２記載のパラメータ推定装置であって、
   前記線形パラメータH^(ω_k,i)がスカラー値の場合、前記重み係数wを
   

   

   
   とし、前記線形パラメータH^(ω_k,i)が長さMのＦＩＲフィルタで構成される場合、前記重み係数wを
   

   

   
   とする、
   パラメータ推定装置。
4. 請求項１から請求項３の何れかに記載のパラメータ推定装置であって、
   前記伝達系の時間領域の出力信号の系列を周波数領域に変換し、前記周波数領域伝達系出力信号を求める第二周波数領域変換部をさらに含む、
   パラメータ推定装置。
5. 過去の入力及び出力に依存しないハードクリップ特性とみなせる非線形系と有限長のインパルス応答とみなせる線形系との縦続接続とみなせる伝達系において、前記伝達系の入力信号と出力信号とを用いて、前記ハードクリップ特性の閾値と前記有限長のインパルス応答とを推定するパラメータ推定装置であって、
   前記閾値の推定値である非線形パラメータを用いて、前記伝達系の時間領域の入力信号に対する前記非線形系の出力信号の推定値である時間領域非線形系出力推定値を求める非線形パラメータ適用部と、
   前記時間領域非線形系出力推定値の前記非線形パラメータに関する偏微分値に相当する値を求める微分信号生成部と、
   前記時間領域非線形系出力推定値の系列及び前記偏微分値に相当する値の系列をそれぞれ周波数領域に変換し、周波数領域の周波数領域非線形系出力推定値及び周波数領域の偏微分値に相当する値である周波数領域偏微分値相当値を求める周波数領域変換部と、
   前記有限長のインパルス応答の周波数領域の推定値である線形パラメータと、前記周波数領域非線形系出力推定値とを乗じ、前記伝達系の出力信号の周波数領域の推定値である周波数領域伝達系出力推定値を求める周波数領域線形パラメータ乗算部と、
   前記周波数領域伝達系出力推定値を時間領域に変換し、時間領域の時間領域伝達系出力推定値を求める時間領域変換部と、
   前記伝達系の時間領域の出力信号と、前記時間領域伝達系出力推定値との差である時間領域差を求める誤差算出部と、
   前記時間領域差の系列を周波数領域に変換し、周波数領域の周波数領域差を求める第三周波数領域変換部と、
   前記非線形パラメータと前記線形パラメータのそれぞれを、前記周波数領域差と前記周波数領域偏微分値相当値と前記周波数領域非線形系出力推定値とを用いて、周波数領域で更新する混合型パラメータ更新部と、を含む、
   パラメータ推定装置。
6. 請求項１から請求項５の何れかに記載のパラメータ推定装置を含むエコー消去装置であって、
   前記非線形系はスピーカ及び増幅器を含むものとし、前記線形系は反響路及びマイクロホンを含むものとし、前記伝達系の入力信号は前記スピーカの再生信号であり、前記伝達系の出力信号は前記マイクロホンの収音信号であり、当該エコー消去装置は請求項１〜４の前記周波数領域差を時間領域に変換した値、または、請求項５の前記時間領域差を送信信号として出力する、
   エコー消去装置。
7. 過去の入力及び出力に依存しないハードクリップ特性とみなせる非線形系と有限長のインパルス応答とみなせる線形系との縦続接続とみなせる伝達系において、前記伝達系の入力信号と出力信号とを用いて、前記ハードクリップ特性の閾値と前記有限長のインパルス応答とを推定するパラメータ推定方法であって、
   非線形パラメータ適用部が、前記閾値の推定値である非線形パラメータを用いて、前記伝達系の時間領域の入力信号に対する前記非線形系の出力信号の推定値である時間領域非線形系出力推定値を求める非線形パラメータ適用ステップと、
   微分信号生成部が、前記時間領域非線形系出力推定値の前記非線形パラメータに関する偏微分値に相当する値を求める微分信号生成ステップと、
   周波数領域変換部が、前記時間領域非線形系出力推定値の系列及び前記偏微分値に相当する値の系列をそれぞれ周波数領域に変換し、周波数領域の周波数領域非線形系出力推定値及び周波数領域の偏微分値に相当する値である周波数領域偏微分値相当値を求める周波数領域変換ステップと、
   線形パラメータ乗算部が、前記有限長のインパルス応答の周波数領域の推定値である線形パラメータと、前記周波数領域非線形系出力推定値とを乗じ、前記伝達系の出力信号の周波数領域の推定値である周波数領域伝達系出力推定値を求める周波数領域線形パラメータ乗算ステップと、
   誤差算出部が、前記伝達系の周波数領域の出力信号である周波数領域伝達系出力信号と、前記周波数領域伝達系出力推定値との差である周波数領域差を求める誤差算出ステップと、
   混合型パラメータ更新部が、前記非線形パラメータと前記線形パラメータのそれぞれを、前記周波数領域差と前記周波数領域偏微分値相当値と前記周波数領域非線形系出力推定値とを用いて、周波数領域で更新する混合型パラメータ更新ステップと、を含む、
   パラメータ推定方法。
8. 過去の入力及び出力に依存しないハードクリップ特性とみなせる非線形系と有限長のインパルス応答とみなせる線形系との縦続接続とみなせる伝達系において、前記伝達系の入力信号と出力信号とを用いて、前記ハードクリップ特性の閾値と前記有限長のインパルス応答とを推定するパラメータ推定方法であって、
   非線形パラメータ適用部が、前記閾値の推定値である非線形パラメータを用いて、前記伝達系の時間領域の入力信号に対する前記非線形系の出力信号の推定値である時間領域非線形系出力推定値を求める非線形パラメータ適用ステップと、
   微分信号生成部が、前記時間領域非線形系出力推定値の前記非線形パラメータに関する偏微分値に相当する値を求める微分信号生成ステップと、
   周波数領域変換部が、前記時間領域非線形系出力推定値の系列及び前記偏微分値に相当する値の系列をそれぞれ周波数領域に変換し、周波数領域の周波数領域非線形系出力推定値及び周波数領域の偏微分値に相当する値である周波数領域偏微分値相当値を求める周波数領域変換ステップと、
   線形パラメータ乗算部が、前記有限長のインパルス応答の周波数領域の推定値である線形パラメータと、前記周波数領域非線形系出力推定値とを乗じ、前記伝達系の出力信号の周波数領域の推定値である周波数領域伝達系出力推定値を求める周波数領域線形パラメータ乗算ステップと、
   時間領域変換部が、前記周波数領域伝達系出力推定値を時間領域に変換し、時間領域の時間領域伝達系出力推定値を求める時間領域変換ステップと、
   誤差算出部が、前記伝達系の時間領域の出力信号と、前記時間領域伝達系出力推定値との差である時間領域差を求める誤差算出ステップと、
   第三周波数領域変換部が、前記時間領域差の系列を周波数領域に変換し、周波数領域の周波数領域差を求める第三周波数領域変換ステップと、
   混合型パラメータ更新部が、前記非線形パラメータと前記線形パラメータのそれぞれを、前記周波数領域差と前記周波数領域偏微分値相当値と前記周波数領域非線形系出力推定値とを用いて、周波数領域で更新する混合型パラメータ更新ステップと、を含む、
   パラメータ推定方法。
9. 請求項１から請求項５の何れかに記載のパラメータ推定装置の各部、または、請求項６に記載のエコー消去装置の各部として、コンピュータを機能させるためのプログラム。

Images