JP3303898B2 - 適応的伝達関数推定方法及びそれを使った推定装置 - Google Patents
適応的伝達関数推定方法及びそれを使った推定装置Info
- Publication number
- JP3303898B2 JP3303898B2 JP17541794A JP17541794A JP3303898B2 JP 3303898 B2 JP3303898 B2 JP 3303898B2 JP 17541794 A JP17541794 A JP 17541794A JP 17541794 A JP17541794 A JP 17541794A JP 3303898 B2 JP3303898 B2 JP 3303898B2
- Authority
- JP
- Japan
- Prior art keywords
- vector
- transfer function
- estimated
- input signal
- coefficient
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Lifetime
Links
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Complex Calculations (AREA)
- Soundproofing, Sound Blocking, And Sound Damping (AREA)
- Filters That Use Time-Delay Elements (AREA)
- Cable Transmission Systems, Equalization Of Radio And Reduction Of Echo (AREA)
Description
【０００１】
【産業上の利用分野】本発明は、音響エコーキャンセ
ラ、能動騒音制御などにおいて、未知系の伝達関数並び
に未知系の出力を射影アルゴリズムにより適応的に推定
する方法及びそれを使った推定装置に関する。
ラ、能動騒音制御などにおいて、未知系の伝達関数並び
に未知系の出力を射影アルゴリズムにより適応的に推定
する方法及びそれを使った推定装置に関する。
【０００２】
【従来の技術】以下の説明では時間の表現は離散時間k
として説明を行う。例えば、時刻k での信号x の振幅値
はx(k)と表す。図１は未知系の伝達関数の推定を説明す
る図であって、１１は伝達関数推定部を、１２は未知系
を、x(k)は未知系への入力信号を、y(k)は未知系からの
出力信号を表す。未知系の伝達関数h(k)はx(k)とy(k)を
用いて推定される。図２は伝達関数の適応的な推定を説
明する図であり、２１は推定伝達関数修正ベクトル計算
部、２２は推定伝達関数修正部、２３は畳み込み演算部
を表し、これらは推定信号生成部２０を構成している。
未知系の伝達関数H(k)は畳み込み演算部２３を構成する
タップ数ＬのFIR フィルタの伝達関数H ＾(k) として推
定される。より具体的には、FIR フィルタの係数h ＾
1(k) …h ＾L(k)を推定することになる。以降は、“伝
達関数”と“FIR フィルタ係数”は同じものを意味する
とする。簡単にするためにフィルタ係数を次式で定義さ
れる推定伝達関数ベクトルH ＾(k) として表す。
として説明を行う。例えば、時刻k での信号x の振幅値
はx(k)と表す。図１は未知系の伝達関数の推定を説明す
る図であって、１１は伝達関数推定部を、１２は未知系
を、x(k)は未知系への入力信号を、y(k)は未知系からの
出力信号を表す。未知系の伝達関数h(k)はx(k)とy(k)を
用いて推定される。図２は伝達関数の適応的な推定を説
明する図であり、２１は推定伝達関数修正ベクトル計算
部、２２は推定伝達関数修正部、２３は畳み込み演算部
を表し、これらは推定信号生成部２０を構成している。
未知系の伝達関数H(k)は畳み込み演算部２３を構成する
タップ数ＬのFIR フィルタの伝達関数H ＾(k) として推
定される。より具体的には、FIR フィルタの係数h ＾
1(k) …h ＾L(k)を推定することになる。以降は、“伝
達関数”と“FIR フィルタ係数”は同じものを意味する
とする。簡単にするためにフィルタ係数を次式で定義さ
れる推定伝達関数ベクトルH ＾(k) として表す。
【０００３】 H ＾(k) ＝ [h ＾1(k),h＾2(k), …,h＾L(k)]T (1) ここで、Ｔは転置を表す。図２において、未知系１２に
対する入力信号x(k)を畳み込み演算部２３に入力し、式
(2) で与えられる出力y ＾(k) y＾(k) ＝ H＾(k)TX(k) (2) X(k) ＝ [x(k),x(k-1),…,x(k-L+1)]T (3) と未知系１２からの出力y(k)を差分器２４で差を求めて
得られる誤差信号e(k)の２乗の期待値が最小になるよう
に推定伝達関数ベクトルH ＾(k) の推定が行われる。こ
こで、y ＾(k) は未知系の出力の推定値であり、推定伝
達関数H ＾(k) が未知の特性に近ければ、y ＾(k) の値
はy(k)に近くなるものである。
対する入力信号x(k)を畳み込み演算部２３に入力し、式
(2) で与えられる出力y ＾(k) y＾(k) ＝ H＾(k)TX(k) (2) X(k) ＝ [x(k),x(k-1),…,x(k-L+1)]T (3) と未知系１２からの出力y(k)を差分器２４で差を求めて
得られる誤差信号e(k)の２乗の期待値が最小になるよう
に推定伝達関数ベクトルH ＾(k) の推定が行われる。こ
こで、y ＾(k) は未知系の出力の推定値であり、推定伝
達関数H ＾(k) が未知の特性に近ければ、y ＾(k) の値
はy(k)に近くなるものである。
【０００４】例えば講演会場や劇場のような音場におい
て人や物体が移動することにより音響路の伝達関数が変
化するように、実際には未知系の伝達関数は時刻と共に
変化する場合が多い。そこで、未知系の伝達関数の推定
値であるH ＾(k) も未知系音響路の変化を検出し適応的
に修正する必要がある。推定伝達関数修正ベクトル計算
部２１では誤差信号e(k)と未知系１２への入力信号x(k)
から推定伝達関数修正ベクトルδH ＾(k) を計算する。
推定伝達関数修正部２２では、次式のように、推定伝達
関数修正ベクトル計算部２１で得られた推定伝達関数修
正ベクトルδH＾(k) を推定伝達関数ベクトルH ＾(k)
に加えることで推定伝達関数を修正する。
て人や物体が移動することにより音響路の伝達関数が変
化するように、実際には未知系の伝達関数は時刻と共に
変化する場合が多い。そこで、未知系の伝達関数の推定
値であるH ＾(k) も未知系音響路の変化を検出し適応的
に修正する必要がある。推定伝達関数修正ベクトル計算
部２１では誤差信号e(k)と未知系１２への入力信号x(k)
から推定伝達関数修正ベクトルδH ＾(k) を計算する。
推定伝達関数修正部２２では、次式のように、推定伝達
関数修正ベクトル計算部２１で得られた推定伝達関数修
正ベクトルδH＾(k) を推定伝達関数ベクトルH ＾(k)
に加えることで推定伝達関数を修正する。
【０００５】 H ＾(k+1) ＝ H＾(k) ＋ μδH ＾(k) (4) ただし、μはステップサイズと呼ばれ、一回毎に行う修
正の幅を制御する予め選択したスカラ量であり、時間に
よらない定数として扱う。以下の説明では推定伝達関数
修正ベクトルδH ＾(k) を求めるにあたってμ=1とした
場合に付いてのδH ＾(k) を求め、そのδH ＾(k) に対
し必要に応じて所望のスッテプサイズμを乗算して使用
する。応用例によってはμの値を時間と共に変化させる
場合もあるが、その場合にも以下の説明は同様に成立す
る。修正された推定伝達関数ベクトルは畳み込み演算部
２３に転送される。以上が時刻k に置ける伝達関数推定
動作であり時刻k+1 以降も同様の動作を繰り返す。
正の幅を制御する予め選択したスカラ量であり、時間に
よらない定数として扱う。以下の説明では推定伝達関数
修正ベクトルδH ＾(k) を求めるにあたってμ=1とした
場合に付いてのδH ＾(k) を求め、そのδH ＾(k) に対
し必要に応じて所望のスッテプサイズμを乗算して使用
する。応用例によってはμの値を時間と共に変化させる
場合もあるが、その場合にも以下の説明は同様に成立す
る。修正された推定伝達関数ベクトルは畳み込み演算部
２３に転送される。以上が時刻k に置ける伝達関数推定
動作であり時刻k+1 以降も同様の動作を繰り返す。
【０００６】図２で説明した推定伝達関数の修正方法は
適応アルゴリズムとして知られている。適応アルゴリズ
ムとしては、LMS アルゴリズムや学習同定法(NLMS:Norm
alized LMS) などがよく知られているが、ここでは文献
「尾関、梅田：“アフィン部分空間への直行射影を用い
た適応フィルタアルゴリズムとその諸性質”、電子情報
通信学会誌(A),J67-App.126-132,(昭59-2)」で提案され
た射影アルゴリズムを対象とする。
適応アルゴリズムとして知られている。適応アルゴリズ
ムとしては、LMS アルゴリズムや学習同定法(NLMS:Norm
alized LMS) などがよく知られているが、ここでは文献
「尾関、梅田：“アフィン部分空間への直行射影を用い
た適応フィルタアルゴリズムとその諸性質”、電子情報
通信学会誌(A),J67-App.126-132,(昭59-2)」で提案され
た射影アルゴリズムを対象とする。
【０００７】射影アルゴリズムは学習同定法に比べて演
算量は増加するものの、音声信号入力に対して優れた適
応性を持つ。射影アルゴリズムでは上述のようにμ=1の
場合の式(4) から次のｐ個の方程式からなる連立方程式
が満足されるようにベクトルδH ＾(k) を決める。 y(k) ＝ X(k)T( H ＾(k)+ δH ＾(k) ） y(k-1) ＝ X(k-1)T( H ＾(k)+ δH ＾(k) ) ： y(k-p+1) ＝ X(k-p+1)T( H ＾(k)+ δH ＾(k) ) (5) 式(5) は時刻ｋにおいて更新された推定伝達関数ベクト
ルH ＾(k+1) が時刻ｋ以前のｐ個の入力ベクトルX(k),X
(k-1),…,X(k+p-1) に対して、それぞれ未知系の出力と
同じ値y(k),y(k-1),…,y(k-p+1) を出力するようにδH
＾(k) を決めることを示している。このようにすれば、
δH ＾(k) で適応的に更新を繰り返すにつれ推定伝達関
数H ＾(k) の特性が未知系の特性に近づくことが期待さ
れる。
算量は増加するものの、音声信号入力に対して優れた適
応性を持つ。射影アルゴリズムでは上述のようにμ=1の
場合の式(4) から次のｐ個の方程式からなる連立方程式
が満足されるようにベクトルδH ＾(k) を決める。 y(k) ＝ X(k)T( H ＾(k)+ δH ＾(k) ） y(k-1) ＝ X(k-1)T( H ＾(k)+ δH ＾(k) ) ： y(k-p+1) ＝ X(k-p+1)T( H ＾(k)+ δH ＾(k) ) (5) 式(5) は時刻ｋにおいて更新された推定伝達関数ベクト
ルH ＾(k+1) が時刻ｋ以前のｐ個の入力ベクトルX(k),X
(k-1),…,X(k+p-1) に対して、それぞれ未知系の出力と
同じ値y(k),y(k-1),…,y(k-p+1) を出力するようにδH
＾(k) を決めることを示している。このようにすれば、
δH ＾(k) で適応的に更新を繰り返すにつれ推定伝達関
数H ＾(k) の特性が未知系の特性に近づくことが期待さ
れる。
【０００８】ｐは射影の次数と呼ばれる量で、次数ｐが
大きいほど適応性は向上するが、より多くの演算量を必
要とする。従来のNLMS法はｐ＝１の場合に対応してい
る。さて、式(5) 中の第１式を移項すれば、 X(k)TδH ＾(k) ＝ y(k)−X(k)TH ＾(k) ＝ e(k) (6) となる。また、式(5) の第２式を移項し、式(4) 及び
(6) のｋをk-1 とおいて得られる式を利用すれば、 X(k-1)TδH ＾(k) ＝ y(k-1)-X(k-1)TH ＾(k) ＝ y(k-1)-X(k-1)T( H ＾(k-1)+ μδH ＾(k-1) ) ＝ y(k-1)-X(k-1)TH ＾(k-1)- μX(k-1)TδH ＾(k-1) ＝ e(k-1)-μe(k-1) ＝ (1-μ)e(k-1) (7) となる。以下同様に、 X(k-i)TδH ＾(k) ＝ (1-μ)ie(k-i) (8) の関係が成立する。このことより、式(5) はδH ＾(k)
に関する次の連立方程式として書き直すことができる。
大きいほど適応性は向上するが、より多くの演算量を必
要とする。従来のNLMS法はｐ＝１の場合に対応してい
る。さて、式(5) 中の第１式を移項すれば、 X(k)TδH ＾(k) ＝ y(k)−X(k)TH ＾(k) ＝ e(k) (6) となる。また、式(5) の第２式を移項し、式(4) 及び
(6) のｋをk-1 とおいて得られる式を利用すれば、 X(k-1)TδH ＾(k) ＝ y(k-1)-X(k-1)TH ＾(k) ＝ y(k-1)-X(k-1)T( H ＾(k-1)+ μδH ＾(k-1) ) ＝ y(k-1)-X(k-1)TH ＾(k-1)- μX(k-1)TδH ＾(k-1) ＝ e(k-1)-μe(k-1) ＝ (1-μ)e(k-1) (7) となる。以下同様に、 X(k-i)TδH ＾(k) ＝ (1-μ)ie(k-i) (8) の関係が成立する。このことより、式(5) はδH ＾(k)
に関する次の連立方程式として書き直すことができる。
【０００９】 Ｘp(k) δH ＾(k) ＝ E(k) (9) ただし、Ｘp(k) はｐ行Ｌ列の行列、E(k) はｐ行のベク
トルで、それぞれ次式で定義される。 E(k)は誤差ベクトルと呼ぶことにする。さて、通常はｐ
＜Ｌであるので、式(9)はδH ＾(k) に関する不定方程
式となり、δH ＾(k) のノルム最小解は、次式で与え
られる。
トルで、それぞれ次式で定義される。 E(k)は誤差ベクトルと呼ぶことにする。さて、通常はｐ
＜Ｌであるので、式(9)はδH ＾(k) に関する不定方程
式となり、δH ＾(k) のノルム最小解は、次式で与え
られる。
【００１０】 δH ＾(k) ＝ Ｘp(k)T( Ｘp(k)Ｘ p(k)T )-1E(k) ＝ [ X(k)X(k-1), …, X(k-p+1) ]G(k) (12) ただし、 G(k) ＝ Ｃp(k)-1E(k) (13) Ｃp(k)＝ Ｘp(k) Ｘp(k)T (14) である。Ｃp(k)はｐ行ｐ列の行列であり、これをｐ次の
共分散行列または自己相関行列と呼び、G(k)はプレフィ
ルタ係数ベクトルと呼ぶ。また、G(k)の要素をg1(k),g2
(k),…,gp(k)と表せば、式(12)より推定伝達関数修正ベ
クトルδH ＾(k)は、 δH ＾(k) ＝Σg i(k)X(k-i+1) （Σはｉ＝１からｐまで） (15) と表される。
共分散行列または自己相関行列と呼び、G(k)はプレフィ
ルタ係数ベクトルと呼ぶ。また、G(k)の要素をg1(k),g2
(k),…,gp(k)と表せば、式(12)より推定伝達関数修正ベ
クトルδH ＾(k)は、 δH ＾(k) ＝Σg i(k)X(k-i+1) （Σはｉ＝１からｐまで） (15) と表される。
【００１１】以上のことにより、図２における推定伝達
関数修正ベクトル計算部２１は、射影法を用いる場合に
は、図３のように表される。図３において、３１はプレ
フィルタ係数ベクトル計算部を表し、入力信号x(k)及び
誤差信号e(k)を用いて、式(13)に基づいてプレフィルタ
係数G(k)を計算する。３２はプレフィルタリング部であ
り、プレフィルタ係数ベクトル計算部３１より転送され
たプレフィルタ係数G(k)を用いて、式(15)で表されるプ
レフィルタリング操作を行って推定伝達関数修正ベクト
ルδH ＾(k) を合成する。
関数修正ベクトル計算部２１は、射影法を用いる場合に
は、図３のように表される。図３において、３１はプレ
フィルタ係数ベクトル計算部を表し、入力信号x(k)及び
誤差信号e(k)を用いて、式(13)に基づいてプレフィルタ
係数G(k)を計算する。３２はプレフィルタリング部であ
り、プレフィルタ係数ベクトル計算部３１より転送され
たプレフィルタ係数G(k)を用いて、式(15)で表されるプ
レフィルタリング操作を行って推定伝達関数修正ベクト
ルδH ＾(k) を合成する。
【００１２】
【発明が解決しようとする課題】以上説明してきた射影
法の演算量を考える。ただし、演算量とは、１離散単位
時間当りの推定動作に必要な積和（または和）演算量を
表すものとする。まず、図２のタップ数Ｌの畳み込み演
算部２３における式(2) の演算量は、Ｌである。次に、
推定伝達関数修正ベクトル計算部２１内のプレフィルタ
係数ベクトル演算部３１（図３）における式(13)の演算
は代表的な数値計算方であるコレスキー法を用いた場合
に約p3/6である。そして、プレフィルタリング部３２に
おける式(15)の演算量は、(p-1)Lである。最後に、推定
伝達関数修正部２２における式(4)の演算量は、Ｌであ
る。従って、射影法の全演算量NCは次式となる。
法の演算量を考える。ただし、演算量とは、１離散単位
時間当りの推定動作に必要な積和（または和）演算量を
表すものとする。まず、図２のタップ数Ｌの畳み込み演
算部２３における式(2) の演算量は、Ｌである。次に、
推定伝達関数修正ベクトル計算部２１内のプレフィルタ
係数ベクトル演算部３１（図３）における式(13)の演算
は代表的な数値計算方であるコレスキー法を用いた場合
に約p3/6である。そして、プレフィルタリング部３２に
おける式(15)の演算量は、(p-1)Lである。最後に、推定
伝達関数修正部２２における式(4)の演算量は、Ｌであ
る。従って、射影法の全演算量NCは次式となる。
【００１３】 NC ＝ L ＋ p3/6 ＋ (p-1)L ＋ L (16) 一方、学習同定法やLMS アルゴリズムの演算量は約２Ｌ
である。一例として、Ｌ＝500,ｐ＝20（音響エコーキャ
ンセラにおける応用の場合の代表値）の場合の演算量を
計算してみる。学習同定法では、演算量が1000であるの
に対して、射影法では、式(16)より、12000 近くもの演
算量が必要となる。特に式(13)の演算量ｐ3/6 はｐが大
きくなると急激に増加する。このように、射影法は、学
習同定法と比べて良好な収束特性を持つ反面、演算量の
増加が問題となる。
である。一例として、Ｌ＝500,ｐ＝20（音響エコーキャ
ンセラにおける応用の場合の代表値）の場合の演算量を
計算してみる。学習同定法では、演算量が1000であるの
に対して、射影法では、式(16)より、12000 近くもの演
算量が必要となる。特に式(13)の演算量ｐ3/6 はｐが大
きくなると急激に増加する。このように、射影法は、学
習同定法と比べて良好な収束特性を持つ反面、演算量の
増加が問題となる。
【００１４】この発明の目的は演算量を削減した未知系
の伝達関数及び出力の適応的推定方法及びそれを使った
推定装置を提供することである。
の伝達関数及び出力の適応的推定方法及びそれを使った
推定装置を提供することである。
【００１５】
【課題を解決するための手段】上述の課題を解決するた
めに、本発明は第１の観点によりプレフィルタ係数ベク
トル計算部３１の演算量を削減し、第２の観点によりプ
レフィルタリング部３２の演算量を削減する。この発明
の第１の観点によれば、離散的時間ｋにおける未知系の
入力信号x(k)と出力y(k)とからその未知系の伝達関数を
推定し、その推定した伝達関数H ＾(k) を持つ疑似系の
出力信号y ＾(k) の上記未知系の出力信号y(k)に対する
誤差信号e(k)=y(k)-y ＾(k) を算出し、その誤差信号ベ
クトルをE(k)、入力信号x(k)の共分散行列をＣp(k)、プ
レフィルタ係数ベクトルをG(k)とするとき、p 元連立１
次方程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求め、
上記プレフィルタ係数ベクトルG(k)を使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の推定伝達関数修正ベクトルδH ＾
(k) を計算し、その推定伝達関数修正ベクトルδH ＾
(k) と、予め決めた修正ステップサイズμとを使って次
式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正することを繰り返し上記誤差信号e(k)がゼロに近づく
ようにする射影アルゴリズムによる適応的伝達関数推定
方法または装置において、上記プレフィルタ計数ベクト
ルG(k)は、上記入力信号x(k)の前向き線形予測係数ベク
トルF(k)とその２乗予測誤差和f(k)と、後向き線形予測
係数ベクトルB(k)とその２乗予測誤差和b(k)を算出し、
プレフィルタ誘導係数ベクトルをD(k)とするとき、次の
第１及び第２式から成る上記プレフィルタ係数ベクトル
についての回帰式
めに、本発明は第１の観点によりプレフィルタ係数ベク
トル計算部３１の演算量を削減し、第２の観点によりプ
レフィルタリング部３２の演算量を削減する。この発明
の第１の観点によれば、離散的時間ｋにおける未知系の
入力信号x(k)と出力y(k)とからその未知系の伝達関数を
推定し、その推定した伝達関数H ＾(k) を持つ疑似系の
出力信号y ＾(k) の上記未知系の出力信号y(k)に対する
誤差信号e(k)=y(k)-y ＾(k) を算出し、その誤差信号ベ
クトルをE(k)、入力信号x(k)の共分散行列をＣp(k)、プ
レフィルタ係数ベクトルをG(k)とするとき、p 元連立１
次方程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求め、
上記プレフィルタ係数ベクトルG(k)を使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の推定伝達関数修正ベクトルδH ＾
(k) を計算し、その推定伝達関数修正ベクトルδH ＾
(k) と、予め決めた修正ステップサイズμとを使って次
式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正することを繰り返し上記誤差信号e(k)がゼロに近づく
ようにする射影アルゴリズムによる適応的伝達関数推定
方法または装置において、上記プレフィルタ計数ベクト
ルG(k)は、上記入力信号x(k)の前向き線形予測係数ベク
トルF(k)とその２乗予測誤差和f(k)と、後向き線形予測
係数ベクトルB(k)とその２乗予測誤差和b(k)を算出し、
プレフィルタ誘導係数ベクトルをD(k)とするとき、次の
第１及び第２式から成る上記プレフィルタ係数ベクトル
についての回帰式
【００１６】
【数５】
【００１７】により上記係数ベクトルG(k)を求める。こ
の発明の第２の観点によれば、離散的時間ｋにおける未
知系の入力信号x(k)と出力y(k)とからその未知系の伝達
関数を推定し、その推定した伝達関数H ＾(k) を持つ疑
似系の出力信号y ＾(k) の上記未知系の出力信号y(k)に
対する誤差信号e(k)=y(k)-y ＾(k) を算出し、その誤差
信号ベクトルをE(k)、入力信号x(k)の共分散行列をＣ
p(k)、プレフィルタ係数ベクトルをG(k)とするとき、p
元連立１次方程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求め、
上記プレフィルタ係数ベクトルG(k)を使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の推定伝達関数修正ベクトルδH ＾
(k) を計算し、その推定伝達関数修正ベクトルδH ＾
(k) と、予め決めた修正ステップサイズμとを使って次
式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正することを繰り返し上記誤差信号e(k)がゼロに近づく
ようにする射影アルゴリズムによる適応的伝達関数推定
方法または装置において、上記修正ベクトルを計算する
代わりに、上記プレフィルタ係数ベクトルG(k)の要素で
あるプレフィルタ係数g(k)を次式 si(k) ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μg1(k) for i = 1 により平滑化して平滑化係数si(k) を得て、上記修正ベ
クトルにより上記推定伝達関数ベクトルを修正する代わ
りに、上記平滑化係数s(k)を使って次式 Z(k+1)＝Z(k)＋sp(k)X(k-p+1) により近似推定伝達関数Z(k+1)を求め、上記近似推定伝
達関数Z(k+1)と上記入力信号x(k)の畳み込み演算X(k)TZ
(k) を行い、上記平滑化係数s(k)のベクトル、上記入力
信号x(k)のベクトル及び上記入力信号の相関ベクトルを
それぞれ Sp-1(k-1)＝[s1(k-1),s2(k-1),…,sp-1(k-1)]T X p-1(k)＝[x(k-1),x(k-2),…,x(k-p+1)]T Rp-1(k-1)＝[X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T として次式の内積 Sp-1(k-1)TRp-1(k) を計算し、上記畳み込み演算結果X(k)TZ(k) と上記内積
演算結果の和を上記推定信号y ＾(k) として出力する。
の発明の第２の観点によれば、離散的時間ｋにおける未
知系の入力信号x(k)と出力y(k)とからその未知系の伝達
関数を推定し、その推定した伝達関数H ＾(k) を持つ疑
似系の出力信号y ＾(k) の上記未知系の出力信号y(k)に
対する誤差信号e(k)=y(k)-y ＾(k) を算出し、その誤差
信号ベクトルをE(k)、入力信号x(k)の共分散行列をＣ
p(k)、プレフィルタ係数ベクトルをG(k)とするとき、p
元連立１次方程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求め、
上記プレフィルタ係数ベクトルG(k)を使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の推定伝達関数修正ベクトルδH ＾
(k) を計算し、その推定伝達関数修正ベクトルδH ＾
(k) と、予め決めた修正ステップサイズμとを使って次
式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正することを繰り返し上記誤差信号e(k)がゼロに近づく
ようにする射影アルゴリズムによる適応的伝達関数推定
方法または装置において、上記修正ベクトルを計算する
代わりに、上記プレフィルタ係数ベクトルG(k)の要素で
あるプレフィルタ係数g(k)を次式 si(k) ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μg1(k) for i = 1 により平滑化して平滑化係数si(k) を得て、上記修正ベ
クトルにより上記推定伝達関数ベクトルを修正する代わ
りに、上記平滑化係数s(k)を使って次式 Z(k+1)＝Z(k)＋sp(k)X(k-p+1) により近似推定伝達関数Z(k+1)を求め、上記近似推定伝
達関数Z(k+1)と上記入力信号x(k)の畳み込み演算X(k)TZ
(k) を行い、上記平滑化係数s(k)のベクトル、上記入力
信号x(k)のベクトル及び上記入力信号の相関ベクトルを
それぞれ Sp-1(k-1)＝[s1(k-1),s2(k-1),…,sp-1(k-1)]T X p-1(k)＝[x(k-1),x(k-2),…,x(k-p+1)]T Rp-1(k-1)＝[X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T として次式の内積 Sp-1(k-1)TRp-1(k) を計算し、上記畳み込み演算結果X(k)TZ(k) と上記内積
演算結果の和を上記推定信号y ＾(k) として出力する。
【００１８】第１の観点に基づく請求項１の方法及び請
求項６の装置によれば、プレフィルタ係数ベクトルG(k)
を射影次数ｐに比例した演算量１５ｐの積和演算で求め
ることができるので、従来の方法で必要とされる演算量
p3/6に比べて著しく少なく、従って装置の規模は小さく
てすみ、また演算に必要とされる時間も短縮される。第
２の観点に基づく請求項２の方法及び請求項７の装置に
よれば、プレフィルタ係数g(k)の平滑化と推定伝達関数
ベクトルH ＾(k) の近似値Z(k)の蓄積とを行うことによ
り演算量の削減がで、従って装置の規模を小さくするこ
とができ、また演算時間を短縮することができる。
求項６の装置によれば、プレフィルタ係数ベクトルG(k)
を射影次数ｐに比例した演算量１５ｐの積和演算で求め
ることができるので、従来の方法で必要とされる演算量
p3/6に比べて著しく少なく、従って装置の規模は小さく
てすみ、また演算に必要とされる時間も短縮される。第
２の観点に基づく請求項２の方法及び請求項７の装置に
よれば、プレフィルタ係数g(k)の平滑化と推定伝達関数
ベクトルH ＾(k) の近似値Z(k)の蓄積とを行うことによ
り演算量の削減がで、従って装置の規模を小さくするこ
とができ、また演算時間を短縮することができる。
【００１９】請求項３の方法及び請求項９の装置は第１
の観点と第２の観点の技術を組み合わせることにより更
に演算量を削減することができる。請求項４の方法は請
求項１又は３の方法において、また請求項１０の装置は
請求項６又は９の装置において、それぞれ後ろ向き線形
予測係数ベクトルB(k-1)の第ｐ要素が常に１であること
に着目し、回帰式における次項 ［B(k-1)TE(k-1)/b(k-1)］B(k-1) の代わりにg(k-1)B(k-1)を使うことにより演算量を削減
する。
の観点と第２の観点の技術を組み合わせることにより更
に演算量を削減することができる。請求項４の方法は請
求項１又は３の方法において、また請求項１０の装置は
請求項６又は９の装置において、それぞれ後ろ向き線形
予測係数ベクトルB(k-1)の第ｐ要素が常に１であること
に着目し、回帰式における次項 ［B(k-1)TE(k-1)/b(k-1)］B(k-1) の代わりにg(k-1)B(k-1)を使うことにより演算量を削減
する。
【００２０】請求項５の方法は請求項１又は３の方法に
おいて、又請求項１１の装置は請求項６又は９の装置に
おいて、回帰式中の次項 ［F(k)TE(k)/f(k)］F(k) 及び ［B(k-1)TE(k-1)/b(k-1)］B(k-1) のそれぞれの分母に非負数のδF(k)，δB(k-1)をそれぞ
れ加算することによりこれらの項の演算の不安定性を軽
減する。
おいて、又請求項１１の装置は請求項６又は９の装置に
おいて、回帰式中の次項 ［F(k)TE(k)/f(k)］F(k) 及び ［B(k-1)TE(k-1)/b(k-1)］B(k-1) のそれぞれの分母に非負数のδF(k)，δB(k-1)をそれぞ
れ加算することによりこれらの項の演算の不安定性を軽
減する。
【００２１】
【実施例】まず、プレフィルタ係数ベクトル計算部３１
の演算量の削減法を示す。この方法はプレフィルタ係数
ベクトルG(k)を１時刻前のG(k-1)を利用した回帰式によ
り求めることを特徴とする。ここでは、まずG(k)がG(k-
1)から回帰的に求まることを証明する。回帰式は式(11)
の誤差ベクトルE(k)の要素が時間と共に1-μを乗じてず
れていく構造を有していることと、式(14)の入力信号の
共分散行列の逆行列の性質を利用して、導出される。
の演算量の削減法を示す。この方法はプレフィルタ係数
ベクトルG(k)を１時刻前のG(k-1)を利用した回帰式によ
り求めることを特徴とする。ここでは、まずG(k)がG(k-
1)から回帰的に求まることを証明する。回帰式は式(11)
の誤差ベクトルE(k)の要素が時間と共に1-μを乗じてず
れていく構造を有していることと、式(14)の入力信号の
共分散行列の逆行列の性質を利用して、導出される。
【００２２】ここで、式(14)と同様にp-1 次の共分散行
列Ｃp-1(k) を次式で定義する。
列Ｃp-1(k) を次式で定義する。
【００２３】
【数６】
【００２４】共分散行列Ｃp(k)の逆行列とＣp-1(k) の
逆行列の間には次の関係があることが文献「S.Haykin.A
daptive Filter Theory,2nd edition,Prentice-Hall.19
91,pp.577-578」で知られている。
逆行列の間には次の関係があることが文献「S.Haykin.A
daptive Filter Theory,2nd edition,Prentice-Hall.19
91,pp.577-578」で知られている。
【００２５】
【数７】
【００２６】ただし、F(k)は正規方程式 Ｃp(k)F(k)＝
[f(k),0,…,0]Tを満たす前向き線形予測係数ベクトル(p
次) であり、先頭の要素は１である。また、f(k)は前向
きの最小２乗予測誤差和である。B(k-1)は正規方程式Ｃ
p(k-1)B(k-1)=[0, …,0,b(k-1)]Tを満たす後向き線形予
測係数ベクトル(p次) であり、最後の要素は１である。
b(k-1)は後向きの最小２乗予測誤差和である。線形予測
係数F(k)，B(k-1)及び２乗予測誤差和f(k)，b(k-1)は例
えばＦＴＦ(Fast Transversal Filters)アルゴリズムを
利用して少ない演算量で求められることがJ.M.Cioffi
と T.Kailath の文献“Windowed fast transversal ada
ptive filter algorithms with normalization." IEEE
Trans. Acoust,Speech.Signal Processing,vol.ASSP-3
3,no.3,pp.607-625 に示されている。式(13)のプレフィ
ルタ係数ベクトルG(k)を次式に再掲する。
[f(k),0,…,0]Tを満たす前向き線形予測係数ベクトル(p
次) であり、先頭の要素は１である。また、f(k)は前向
きの最小２乗予測誤差和である。B(k-1)は正規方程式Ｃ
p(k-1)B(k-1)=[0, …,0,b(k-1)]Tを満たす後向き線形予
測係数ベクトル(p次) であり、最後の要素は１である。
b(k-1)は後向きの最小２乗予測誤差和である。線形予測
係数F(k)，B(k-1)及び２乗予測誤差和f(k)，b(k-1)は例
えばＦＴＦ(Fast Transversal Filters)アルゴリズムを
利用して少ない演算量で求められることがJ.M.Cioffi
と T.Kailath の文献“Windowed fast transversal ada
ptive filter algorithms with normalization." IEEE
Trans. Acoust,Speech.Signal Processing,vol.ASSP-3
3,no.3,pp.607-625 に示されている。式(13)のプレフィ
ルタ係数ベクトルG(k)を次式に再掲する。
【００２７】 G(k) ＝ Ｃp(k)-1E(k) (20) このG(k)に対応してp-1 次のプレフィルタ誘導係数ベク
トルD(k)を次式で定義する。 D(k) ＝ Ｃp-1(k)-1Ep-1(k) (21) ここで、 である。式(22)においてｋを(k-1) とおいて左右両辺に
(1-μ)を乗算して得られるp-1 次のベクトル(1-μ)Ep-1
(k-1) の全要素は式(11)におけるｐ次のベクトルE(k)の
先頭の要素e(k)以外の要素を構成しているので次式の関
係が成立する。
トルD(k)を次式で定義する。 D(k) ＝ Ｃp-1(k)-1Ep-1(k) (21) ここで、 である。式(22)においてｋを(k-1) とおいて左右両辺に
(1-μ)を乗算して得られるp-1 次のベクトル(1-μ)Ep-1
(k-1) の全要素は式(11)におけるｐ次のベクトルE(k)の
先頭の要素e(k)以外の要素を構成しているので次式の関
係が成立する。
【００２８】 さらに、式(22)のp-1 次のベクトルの全要素は式(11)に
おけるｐ次のベクトルの最後の要素を除く全要素を構成
しているので、E(k)とEp-1(k) の関係は次式のようにも
表せる。
おけるｐ次のベクトルの最後の要素を除く全要素を構成
しているので、E(k)とEp-1(k) の関係は次式のようにも
表せる。
【００２９】 式(18)の両辺のそれぞれの項に右からE(k)を乗ずると、
式(20),(21),(23)より、
式(20),(21),(23)より、
【００３０】
【数８】
【００３１】が得られる。また、式(19)の両辺に右から
E(k-1)を乗ずると、式(20),(21),(24)より、
E(k-1)を乗ずると、式(20),(21),(24)より、
【００３２】
【数９】
【００３３】となる。移項すると、
【００３４】
【数１０】
【００３５】このように、G(k)は、式(25)で表されるよ
うにf(k-1)の値に基づいて計算されるが、式(27)で表さ
れるようにD(k-1)はG(k-1)から計算されている。即ち、
式(25)，(27)はG(k)についての回帰式になっている。さ
て、式(25)の右辺において、第１項の(1-μ)［］を得る
ためにはp-1 回の乗算を必要とし、第２項のF(k)TE(k)
はp 回の積和演算、これをf(k)で除するのに１回の除
算、更に、この結果をF(k) に乗じて右辺第１項に加え
るのにp 回の積和計算が必要である。つまり、演算の種
類を区別しなければ3p回の演算が必要とされる。同様
に、式(27)を得るには約2p回の演算が必要である。従っ
て、式(16)，(17)の演算量は約5pであり、p に比例した
量である。一方、式(25),(27) で必要とされる線形予測
係数ベクトルF(k)，B(k-1)及び、最小２乗予測誤差和f
(k)，b(k-1)は線形予測手法を用いれば、約10p の演算
量で計算できることが前述のJ.M.CioffiとT. Kailathの
文献で知られている。従って、式(25),(27) の回帰式に
基づけば、次数p に比例した演算量でG(k)が計算できる
ことが明らかである。
うにf(k-1)の値に基づいて計算されるが、式(27)で表さ
れるようにD(k-1)はG(k-1)から計算されている。即ち、
式(25)，(27)はG(k)についての回帰式になっている。さ
て、式(25)の右辺において、第１項の(1-μ)［］を得る
ためにはp-1 回の乗算を必要とし、第２項のF(k)TE(k)
はp 回の積和演算、これをf(k)で除するのに１回の除
算、更に、この結果をF(k) に乗じて右辺第１項に加え
るのにp 回の積和計算が必要である。つまり、演算の種
類を区別しなければ3p回の演算が必要とされる。同様
に、式(27)を得るには約2p回の演算が必要である。従っ
て、式(16)，(17)の演算量は約5pであり、p に比例した
量である。一方、式(25),(27) で必要とされる線形予測
係数ベクトルF(k)，B(k-1)及び、最小２乗予測誤差和f
(k)，b(k-1)は線形予測手法を用いれば、約10p の演算
量で計算できることが前述のJ.M.CioffiとT. Kailathの
文献で知られている。従って、式(25),(27) の回帰式に
基づけば、次数p に比例した演算量でG(k)が計算できる
ことが明らかである。
【００３６】以上、説明した論理に基づいたこの発明で
は、以下の事項(A),(B),(C) が演算量の削減、演算の安
定性に有効である。 (A) 式(27)はp 次のベクトルに関する等式であり、ベク
トルの第p 要素（最下の要素）に関しては左辺は０であ
る。また、ベクトルB(k-1)の第ｐ要素は常に１なのでG
(k-1)の第ｐ要素であるgp(k-1) を用いて次式が成立す
る。
は、以下の事項(A),(B),(C) が演算量の削減、演算の安
定性に有効である。 (A) 式(27)はp 次のベクトルに関する等式であり、ベク
トルの第p 要素（最下の要素）に関しては左辺は０であ
る。また、ベクトルB(k-1)の第ｐ要素は常に１なのでG
(k-1)の第ｐ要素であるgp(k-1) を用いて次式が成立す
る。
【００３７】 B(k-1)TE(k-1)／b(k-1)＝ gp(k-1) (28) 従って、式(27)の代わりに次式 を計算してもよい。これによって演算量を更に削減でき
る。ただし、演算誤差を考慮するとgp(k-1) と式(28)の
左辺と２通り求めて平均する方が有利な場合もある。
る。ただし、演算誤差を考慮するとgp(k-1) と式(28)の
左辺と２通り求めて平均する方が有利な場合もある。
【００３８】(B) 式(25),(27) の右辺第２項の分母f
(k),b(k-1) が小さくなると演算が不安定になるので、
次式
(k),b(k-1) が小さくなると演算が不安定になるので、
次式
【００３９】
【数１１】
【００４０】に示すように、分母に非負数のδF(k),δB
(k-1) を加えることで数値演算的な不安定性を軽減でき
る。δF(k),δB(k-1) の実際的な値としては、入力信号
x(k)の平均パワーより４０dB程度小さい（即ち平均パワ
ーの1/10000 程度）所望の一定値としてもよいし、パワ
ーの変動に従って時間k と共に変化させてもよい。
(k-1) を加えることで数値演算的な不安定性を軽減でき
る。δF(k),δB(k-1) の実際的な値としては、入力信号
x(k)の平均パワーより４０dB程度小さい（即ち平均パワ
ーの1/10000 程度）所望の一定値としてもよいし、パワ
ーの変動に従って時間k と共に変化させてもよい。
【００４１】(C) 式(25),(27) で必要とされる線形予測
係数ベクトルF(k)，B(k-1)及び、最小２乗予測誤差和f
(k),b(k-1) を求める線形予測分析では、入力信号x(k)
の分析区間が通常の場合と異なっている。具体的に説明
すると、通常の線形予測分析における分析区間は時刻０
から、現時刻までであり、時刻k-1 がｋに更新された時
には分析区間にx(k)を加えるだけでよい。一方、射影ア
ルゴリズムの場合は分析区間は式(5) のx(k),…,x(k-L-
p+2)であるので、時刻k-1 がｋに更新されたときには、
分析区間にx(k)を加えるだけではなく、x(k-L-p+1)を除
く必要がある。このため、射影アルゴリスムに於ける線
形予測分析には通常の線形予測分析の場合と比べて２倍
の演算量を要する。しかし、入力信号の統計的性質が変
化しないか、または、変化してもその変化が遅いと期待
され得る場合には、線形予測分析の結果が分析区間にあ
まり依存しないので、時刻ｋでは分析区間にx(k)を加え
るだけでの通常の線形予測分析を行えばよく、演算量を
削減することができる。
係数ベクトルF(k)，B(k-1)及び、最小２乗予測誤差和f
(k),b(k-1) を求める線形予測分析では、入力信号x(k)
の分析区間が通常の場合と異なっている。具体的に説明
すると、通常の線形予測分析における分析区間は時刻０
から、現時刻までであり、時刻k-1 がｋに更新された時
には分析区間にx(k)を加えるだけでよい。一方、射影ア
ルゴリズムの場合は分析区間は式(5) のx(k),…,x(k-L-
p+2)であるので、時刻k-1 がｋに更新されたときには、
分析区間にx(k)を加えるだけではなく、x(k-L-p+1)を除
く必要がある。このため、射影アルゴリスムに於ける線
形予測分析には通常の線形予測分析の場合と比べて２倍
の演算量を要する。しかし、入力信号の統計的性質が変
化しないか、または、変化してもその変化が遅いと期待
され得る場合には、線形予測分析の結果が分析区間にあ
まり依存しないので、時刻ｋでは分析区間にx(k)を加え
るだけでの通常の線形予測分析を行えばよく、演算量を
削減することができる。
【００４２】前述の理論的考察に基づいたこの発明の適
応的伝達関数推定方法の実施例を以下で説明するが、そ
の機能ブロックの全体的構成は図２と同じなので、これ
を参照するものとする。この発明が基づく前述の理論は
図２中の推定伝達関数修正ベクトル計算部２１、特にそ
の中の図３に示すプレフィルタ係数ベクトル計算部３１
とプレフィルタリング部３２における計算量の削減に係
わっているので、これらに付いて以下に詳細に説明す
る。
応的伝達関数推定方法の実施例を以下で説明するが、そ
の機能ブロックの全体的構成は図２と同じなので、これ
を参照するものとする。この発明が基づく前述の理論は
図２中の推定伝達関数修正ベクトル計算部２１、特にそ
の中の図３に示すプレフィルタ係数ベクトル計算部３１
とプレフィルタリング部３２における計算量の削減に係
わっているので、これらに付いて以下に詳細に説明す
る。
【００４３】図４は前述の考察に基づくプレフィルタ係
数ベクトル計算部３１の構成である。４１は線形予測部
を、４２はプレフィルタ誘導係数ベクトル修正部を、４
３はプレフィルタ係数ベクトル修正部を、４４は誤差ベ
クトル生成部を表わす。線形予測部４１では、正規方程
式Ｃp(k)F(k)=[f(k),0,…,0]T を満たす前向きの線形予
測係数ベクトルF(k)及びその予測係数F(k)を用いたとき
の２乗予測誤差和f(k)を計算し、また正規方程式Ｃp(k-
1)B(k-1)=[0, …,0,b(k-1)]Tを満たす後向きの線形予測
係数ベクトルB(k),及び、その予測係数B(k) を用いたと
きの２乗予測誤差和b(k) を計算する。これらの計算方
法は前述のJ.M.Cioffi らの文献に示されている方法を
使えばよい。誤差信号ベクトル生成部４４はｐ個の誤差
信号e(k)，e(k-1),…,e(k-p+1)を蓄積し、式(11)の誤差
信号ベクトルE(k)を構成する。プレフィルタ誘導係数ベ
クトル修正部４２ではプレフィルタ係数ベクトルG(k-
1)、後向きの線形予測係数ベクトルB(k-1)、誤差信号ベ
クトルE(k-1)、及び後向きの最小２乗予測誤差和b(k-1)
を用いて、式(27)又は(30)に基づいてプレフィルタ誘導
係数ベクトルD(k-1)を算出する。プレフィルタ係数ベク
トル修正部４３は式(13)を満たすG(k)をプレフィルタ誘
導ベクトルD(k-1)、誤差信号ベクトルE(k)、前向きの線
形予測係数ベクトルF(k)、及び前向きの最小２乗予測誤
差和f(k)を用いて、式(25)又は(29)により算出する。
数ベクトル計算部３１の構成である。４１は線形予測部
を、４２はプレフィルタ誘導係数ベクトル修正部を、４
３はプレフィルタ係数ベクトル修正部を、４４は誤差ベ
クトル生成部を表わす。線形予測部４１では、正規方程
式Ｃp(k)F(k)=[f(k),0,…,0]T を満たす前向きの線形予
測係数ベクトルF(k)及びその予測係数F(k)を用いたとき
の２乗予測誤差和f(k)を計算し、また正規方程式Ｃp(k-
1)B(k-1)=[0, …,0,b(k-1)]Tを満たす後向きの線形予測
係数ベクトルB(k),及び、その予測係数B(k) を用いたと
きの２乗予測誤差和b(k) を計算する。これらの計算方
法は前述のJ.M.Cioffi らの文献に示されている方法を
使えばよい。誤差信号ベクトル生成部４４はｐ個の誤差
信号e(k)，e(k-1),…,e(k-p+1)を蓄積し、式(11)の誤差
信号ベクトルE(k)を構成する。プレフィルタ誘導係数ベ
クトル修正部４２ではプレフィルタ係数ベクトルG(k-
1)、後向きの線形予測係数ベクトルB(k-1)、誤差信号ベ
クトルE(k-1)、及び後向きの最小２乗予測誤差和b(k-1)
を用いて、式(27)又は(30)に基づいてプレフィルタ誘導
係数ベクトルD(k-1)を算出する。プレフィルタ係数ベク
トル修正部４３は式(13)を満たすG(k)をプレフィルタ誘
導ベクトルD(k-1)、誤差信号ベクトルE(k)、前向きの線
形予測係数ベクトルF(k)、及び前向きの最小２乗予測誤
差和f(k)を用いて、式(25)又は(29)により算出する。
【００４４】以上の方法に基づけば、プレフィルタ係数
ベクトル演算部３１で必要とされる演算量はp3/6からは
15p と大幅に低減できる。しかし、プレフィルタリング
部３２における演算量(p-1)Lが残されているために、L
が大きな場合には依然として大きな演算量が必要である
という問題点が残されている。次に、上記の問題に対し
て、推定伝達関数ベクトルH ＾(k+1) の近似値の蓄積、
ならびにプレフィルタリング係数の平滑化によって解決
する方法を説明する。
ベクトル演算部３１で必要とされる演算量はp3/6からは
15p と大幅に低減できる。しかし、プレフィルタリング
部３２における演算量(p-1)Lが残されているために、L
が大きな場合には依然として大きな演算量が必要である
という問題点が残されている。次に、上記の問題に対し
て、推定伝達関数ベクトルH ＾(k+1) の近似値の蓄積、
ならびにプレフィルタリング係数の平滑化によって解決
する方法を説明する。
【００４５】まず、式(4) のk をk-1 とおくと、推定伝
達関数H ＾(k) は、 H ＾(k) ＝ H＾(k-1)+μδH ＾(k-1) (31) と表される。これを式(4) に代入することで、推定伝達
関数H ＾(k+1) は、 H ＾(k+1) ＝ H＾(k-1)+μδH ＾(k)+ μδH ＾(k-1) (32) と表される。更に、式(4) のk をk-2,k-3,…といおて同
様な代入を繰り返せば、推定伝達関数H ＾(k+1) は、 H ＾(k+1) ＝ μδH ＾(k)+ μδH ＾(k-1)+… +μδH ＾(0) (33) と表される。ただし、推定の初期値H ＾(0) は０とし、
この式により、推定伝達関数H ＾(k+1) は、時刻０（伝
達関数推定の開始時刻）から現時刻k までの修正ベクト
ルμδH ＾(k),μδH ＾(k-1),…, μδH ＾(0) を累積
加算したものであることがわかる。
達関数H ＾(k) は、 H ＾(k) ＝ H＾(k-1)+μδH ＾(k-1) (31) と表される。これを式(4) に代入することで、推定伝達
関数H ＾(k+1) は、 H ＾(k+1) ＝ H＾(k-1)+μδH ＾(k)+ μδH ＾(k-1) (32) と表される。更に、式(4) のk をk-2,k-3,…といおて同
様な代入を繰り返せば、推定伝達関数H ＾(k+1) は、 H ＾(k+1) ＝ μδH ＾(k)+ μδH ＾(k-1)+… +μδH ＾(0) (33) と表される。ただし、推定の初期値H ＾(0) は０とし、
この式により、推定伝達関数H ＾(k+1) は、時刻０（伝
達関数推定の開始時刻）から現時刻k までの修正ベクト
ルμδH ＾(k),μδH ＾(k-1),…, μδH ＾(0) を累積
加算したものであることがわかる。
【００４６】さて、修正ベクトルμδH ＾(k) は、式(1
5)で表される。式(15)のk をk-1,k-2,…とおいて式(33)
に代入すれば、 H ＾(k+1) ＝μ[｛g1(k)X(k)+g2(k)X(k-1)+ … +gp(k)X(k-p+1)｝ +｛g1(k-1)X(k-1)+g2(k-1)X(k-2)+ … +gp(k-1)X(k-p)｝ + … +｛g1(0)X(0)+g2(0)X(-1)+ … +gp(0)X(-p+1)｝] ＝μ[g1(k)X(k) +｛g2(k)+g1(k-1)｝X(k-1) +｛g3(k)+g2(k-1)+g1(k-2)｝X(k-2) + … +｛gp-1(k)+gp-2(k-1)+ … +g1(k-p+2)｝X(k-p+2) +｛gp(k)+gp-1(k-1)+ … +g2(k-p+2)+g1(k-p+1)｝X(k-p+1) +｛gp(k-1)+gp-1(k-2)+ … +g2(k-p+1)+g1(k-p)｝X(k-p) +｛gp(k-2)+gp-1(k-3)+ … +g2(k-p)+g1(k-p-1)｝X(k-p-1) + … (34) と表される。
5)で表される。式(15)のk をk-1,k-2,…とおいて式(33)
に代入すれば、 H ＾(k+1) ＝μ[｛g1(k)X(k)+g2(k)X(k-1)+ … +gp(k)X(k-p+1)｝ +｛g1(k-1)X(k-1)+g2(k-1)X(k-2)+ … +gp(k-1)X(k-p)｝ + … +｛g1(0)X(0)+g2(0)X(-1)+ … +gp(0)X(-p+1)｝] ＝μ[g1(k)X(k) +｛g2(k)+g1(k-1)｝X(k-1) +｛g3(k)+g2(k-1)+g1(k-2)｝X(k-2) + … +｛gp-1(k)+gp-2(k-1)+ … +g1(k-p+2)｝X(k-p+2) +｛gp(k)+gp-1(k-1)+ … +g2(k-p+2)+g1(k-p+1)｝X(k-p+1) +｛gp(k-1)+gp-1(k-2)+ … +g2(k-p+1)+g1(k-p)｝X(k-p) +｛gp(k-2)+gp-1(k-3)+ … +g2(k-p)+g1(k-p-1)｝X(k-p-1) + … (34) と表される。
【００４７】この式から以下のことが予想される。まず
第１に、プレフィルタ係数gi(k) は、図３に示すプレフ
ィルタ係数計算部３１において各時刻ｋ毎に計算され、
プレフィルタリング部３３に供給されるが、このプレフ
ィルタ係数gi(k) を平滑化（または平均化）することで
演算量の低減が予想される。第２に、式(34)における+
｛gp(k-1)+… 以降の項（下から２行目以降の項）に
は、時刻ｋで供給されるプレフィルタ係数gi(k) は含ま
れていないので、この部分は時刻k-1 の時の値と同じで
ある。同様に、+｛gp(k)+ … 以降の項（下から３行目
以降の項）は、時刻ｋ以降は変化しない。従って、+｛g
p(k)+ … 以降の項をまとめて、これをH ＾(k+1) の近
似値として蓄積しておけば、この項は時刻ｋ以降におい
ては計算する必要がなく、演算量の低減が予想される。
第１に、プレフィルタ係数gi(k) は、図３に示すプレフ
ィルタ係数計算部３１において各時刻ｋ毎に計算され、
プレフィルタリング部３３に供給されるが、このプレフ
ィルタ係数gi(k) を平滑化（または平均化）することで
演算量の低減が予想される。第２に、式(34)における+
｛gp(k-1)+… 以降の項（下から２行目以降の項）に
は、時刻ｋで供給されるプレフィルタ係数gi(k) は含ま
れていないので、この部分は時刻k-1 の時の値と同じで
ある。同様に、+｛gp(k)+ … 以降の項（下から３行目
以降の項）は、時刻ｋ以降は変化しない。従って、+｛g
p(k)+ … 以降の項をまとめて、これをH ＾(k+1) の近
似値として蓄積しておけば、この項は時刻ｋ以降におい
ては計算する必要がなく、演算量の低減が予想される。
【００４８】次に、以上のことを数式により表現する。
まず、プレフィルタ係数の平滑化は、対応する入力ベク
トルX(k-i)毎に行う。式(34)より、例えば、X(k-1)に対
応する平滑化は、g2(k)+g1(k-1) であり、また、X(k-2)
に対応する平滑化はg3(k)+g2(k-1)+g1(k-2) である。X
(k-i+1)に対応する平滑化の結果（平滑化係数）を定数
μも含めてsi(k) と表せば、 si(k) ＝ μΣ gi-j(k-j) , 1≦i≦p （Σはｊ＝０からｉ−１まで） (35) si(k) ＝ μΣ gi-j(k-j) , p＜i （Σはｊ＝ｉ−ｐからｉ−１まで） (36) と表される。さらに、式(35)は、 si(k) ＝μΣ gi-j(k-j)＋μgi(k) （Σはｊ＝１からｉ−１まで） ＝μΣ gi-j-1(k-j-1)＋μgi(k) （Σはｊ＝０からｉ−２まで） ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μgi(k) for i = 1 (37) と表される。
まず、プレフィルタ係数の平滑化は、対応する入力ベク
トルX(k-i)毎に行う。式(34)より、例えば、X(k-1)に対
応する平滑化は、g2(k)+g1(k-1) であり、また、X(k-2)
に対応する平滑化はg3(k)+g2(k-1)+g1(k-2) である。X
(k-i+1)に対応する平滑化の結果（平滑化係数）を定数
μも含めてsi(k) と表せば、 si(k) ＝ μΣ gi-j(k-j) , 1≦i≦p （Σはｊ＝０からｉ−１まで） (35) si(k) ＝ μΣ gi-j(k-j) , p＜i （Σはｊ＝ｉ−ｐからｉ−１まで） (36) と表される。さらに、式(35)は、 si(k) ＝μΣ gi-j(k-j)＋μgi(k) （Σはｊ＝１からｉ−１まで） ＝μΣ gi-j-1(k-j-1)＋μgi(k) （Σはｊ＝０からｉ−２まで） ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μgi(k) for i = 1 (37) と表される。
【００４９】一方、H ＾(k+1) の近似値をZ(k+1)と表せ
ば、 Z(k+1) ＝ ｛gp(k)+gp-1(k-1)+ … +g2(k-p+2)+g1(k-p+1)｝X(k-p+1) +｛gp(k-1)+gp-1(k-2)+ … +g2(k-p+1)+g1(k-p)｝X(k-p) +｛gp(k-2)+gp-1(k-3)+ … +g2(k-p)+g1(k-p-1)｝X(k-p-1) + … (38) となり、推定伝達関数H ＾(k+1) は式(34),(35),(38)よ
り、 H ＾(k+1) ＝ Z(k+1) ＋ Σ si(k)X(k-i+1) (39) （Σはｉ＝１からｐ−１まで） となる。また、Z(k+1)とZ(k)の間には、 Z(k+1) ＝ ｛gp(k)+gp-1(k-1)+ … +g2(k-p+2)+g1(k-p+1)｝X(k-p+1) + Z(k) ＝ sp(k)X(k-p+1)＋Z(k) (40) の関係が成立する。
ば、 Z(k+1) ＝ ｛gp(k)+gp-1(k-1)+ … +g2(k-p+2)+g1(k-p+1)｝X(k-p+1) +｛gp(k-1)+gp-1(k-2)+ … +g2(k-p+1)+g1(k-p)｝X(k-p) +｛gp(k-2)+gp-1(k-3)+ … +g2(k-p)+g1(k-p-1)｝X(k-p-1) + … (38) となり、推定伝達関数H ＾(k+1) は式(34),(35),(38)よ
り、 H ＾(k+1) ＝ Z(k+1) ＋ Σ si(k)X(k-i+1) (39) （Σはｉ＝１からｐ−１まで） となる。また、Z(k+1)とZ(k)の間には、 Z(k+1) ＝ ｛gp(k)+gp-1(k-1)+ … +g2(k-p+2)+g1(k-p+1)｝X(k-p+1) + Z(k) ＝ sp(k)X(k-p+1)＋Z(k) (40) の関係が成立する。
【００５０】さて、未知系の出力の推定値y ＾(k) は式
(2),(39)より、 y ＾(k) ＝ X(k)TH ＾(k) ＝ X(k)T｛Z(k)＋Σ si(k-1)X(k-i)｝ （Σはｉ＝１からｐ−１まで） ＝ X(k)TZ(k) ＋Σ si(k-1)X(k)TX(k-i) （Σはｉ＝１からｐ−１ まで） ＝ X(k)TZ(k) ＋ Sp-1(k-1)TRp-1(k) (41) と表される。ただし、Sp-1(k-1) は平滑化係数ベクト
ル、Rp-1(k) は相関ベクトルであり、それぞれ式で定義
される。
(2),(39)より、 y ＾(k) ＝ X(k)TH ＾(k) ＝ X(k)T｛Z(k)＋Σ si(k-1)X(k-i)｝ （Σはｉ＝１からｐ−１まで） ＝ X(k)TZ(k) ＋Σ si(k-1)X(k)TX(k-i) （Σはｉ＝１からｐ−１ まで） ＝ X(k)TZ(k) ＋ Sp-1(k-1)TRp-1(k) (41) と表される。ただし、Sp-1(k-1) は平滑化係数ベクト
ル、Rp-1(k) は相関ベクトルであり、それぞれ式で定義
される。
【００５１】 Sp-1(k-1) ＝ [s1(k-1),s2(k-1),…,sp-1(k-1)]T (42) Rp-1(k) ＝ [X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T (43) ベクトルX(k)は式(3) で定義したように、 X(k) ＝ [x(k),x(k-1),…,x(k-L+1)]T (44) であるので、 X(k)TX(k-i) ＝ X(k-1)TX(k-i-1)−x(k-L)x(k-L-i)＋x(k)x(k-i) i = 1,2,…,p-1 (45) の関係が成立し、 Rp-1(k) ＝ Rp-1(k-1)−x(k-L)Xp-1(k-L)＋x(k)Xp-1(k) (46) が成立する。但し、 Xp-1(k) ＝ [x(k-1),x(k-2),…,x(k-p+1)]T (47) である。
【００５２】ここで、推定伝達関数H ＾(k) の近似値Z
(k)の蓄積、及びプレフィルタ係数の平滑化を行った場
合の伝達関数推定手順を、図５に基づいて説明する。ま
ず、時刻ｋにおいて、入力信号x(k)が供給される相関計
算部５２において、このx(k)及び過去の入力値x(k-1),
…,x(k-L)、及び前の時刻の相関ベクトルRp-1(k-1) を
用いて、式(46)に基づいて、相関ベクトルRp-1(k) を計
算する。次に、この相関ベクトルRp-1(k) とプレフィル
タ係数の平滑化係数ベクトルSp-1(k-1) との内積Sp-1(k
-1)TRp-1(k) を内積演算計算部５３において計算する。
また、畳み込み部５４において、蓄積しておいた推定伝
達関数の近似値Z(k)と入力信号との畳み込み演算X(k)TZ
(k) を行う。前記内積演算及び畳み込み演算の結果は、
加算部５７において加算され、未知系出力の推定値y ＾
(k) を合成する。これらの操作は、式(41)の演算に相当
する。
(k)の蓄積、及びプレフィルタ係数の平滑化を行った場
合の伝達関数推定手順を、図５に基づいて説明する。ま
ず、時刻ｋにおいて、入力信号x(k)が供給される相関計
算部５２において、このx(k)及び過去の入力値x(k-1),
…,x(k-L)、及び前の時刻の相関ベクトルRp-1(k-1) を
用いて、式(46)に基づいて、相関ベクトルRp-1(k) を計
算する。次に、この相関ベクトルRp-1(k) とプレフィル
タ係数の平滑化係数ベクトルSp-1(k-1) との内積Sp-1(k
-1)TRp-1(k) を内積演算計算部５３において計算する。
また、畳み込み部５４において、蓄積しておいた推定伝
達関数の近似値Z(k)と入力信号との畳み込み演算X(k)TZ
(k) を行う。前記内積演算及び畳み込み演算の結果は、
加算部５７において加算され、未知系出力の推定値y ＾
(k) を合成する。これらの操作は、式(41)の演算に相当
する。
【００５３】次に図２に示した減算器２４で未知系出力
y(k)と推定値y ＾(k) との誤差e(k)を求め、図４に示し
たプレフィルタ係数ベクトル演算部３１で、プレフィル
タ係数gi(k) を計算する。次に、計算されたプレフィル
タ係数を図５におけるプレフィルタ係数平滑部５１に送
出し、平滑化を行ってｐ個の平滑化係数s1(k),s2(k),
…,sp-1(k),sp(k)を合成する。この平滑化操作は式(37)
に基づいて行う。これらのうちs1(k),s2(k),…,sp-1(k)
は、平滑化係数ベクトルSp-1(k) の要素として内積演算
計算部５３及び推定伝達関数計算部５６に供給する。ま
た、sp(k) は推定伝達関数の近似値蓄積部５５に供給さ
れる。
y(k)と推定値y ＾(k) との誤差e(k)を求め、図４に示し
たプレフィルタ係数ベクトル演算部３１で、プレフィル
タ係数gi(k) を計算する。次に、計算されたプレフィル
タ係数を図５におけるプレフィルタ係数平滑部５１に送
出し、平滑化を行ってｐ個の平滑化係数s1(k),s2(k),
…,sp-1(k),sp(k)を合成する。この平滑化操作は式(37)
に基づいて行う。これらのうちs1(k),s2(k),…,sp-1(k)
は、平滑化係数ベクトルSp-1(k) の要素として内積演算
計算部５３及び推定伝達関数計算部５６に供給する。ま
た、sp(k) は推定伝達関数の近似値蓄積部５５に供給さ
れる。
【００５４】推定伝達関数の近似値蓄積部５５では、sp
(k) 及び入力信号ベクトルX(k)を用いて、近似値の更新
を行う。即ち、これまで蓄積されていた１次近似値Z(k)
にたいしてsp(k)X(k-p+1) を加算し、これをZ(k+1)とし
て蓄積する。この操作は式(40)の演算に対応する。最後
に、推定伝達関数計算部５６において、平滑化係数ベク
トルSp-1(k) 及び入力信号ベクトルX(k)を用いて、式(3
9)の計算を行って推定伝達関数H ＾(k+1)を得る。
(k) 及び入力信号ベクトルX(k)を用いて、近似値の更新
を行う。即ち、これまで蓄積されていた１次近似値Z(k)
にたいしてsp(k)X(k-p+1) を加算し、これをZ(k+1)とし
て蓄積する。この操作は式(40)の演算に対応する。最後
に、推定伝達関数計算部５６において、平滑化係数ベク
トルSp-1(k) 及び入力信号ベクトルX(k)を用いて、式(3
9)の計算を行って推定伝達関数H ＾(k+1)を得る。
【００５５】以上の動作において、各部における演算量
の概算値は、 相関計算部５２ 2p 内積演算部５３ p 畳み込み部５４ L プレフィルタ係数ベクトル演算部３１ 15p プレフィルタ係数平滑部５１ p １次近似値蓄積部５５ L 推定伝達関数計算部５６ Lp となる。ただし、p−1≒p とみなした。これらを足し算
した総演算量は、 (L+19)p＋2L (48) となる。
の概算値は、 相関計算部５２ 2p 内積演算部５３ p 畳み込み部５４ L プレフィルタ係数ベクトル演算部３１ 15p プレフィルタ係数平滑部５１ p １次近似値蓄積部５５ L 推定伝達関数計算部５６ Lp となる。ただし、p−1≒p とみなした。これらを足し算
した総演算量は、 (L+19)p＋2L (48) となる。
【００５６】ここで、注目すべき点は以下の点である。
図２に示した従来の伝達関数推定法では、伝達関数に相
当する推定伝達関数H ＾(k+1) を毎時刻計算し、これに
基づいて未知系の出力の推定値y ＾(k) を合成してき
た。これに対して、推定伝達関数の近似値Z(k)を用いた
本発明では、図５に示したように、推定伝達関数H ＾(k
+1) を計算しなくても、未知系の出力の推定値y ＾(k)
は合成できる。そして、y ＾(k) が合成できれば、上記
一連の動作を行うことができる。従って、推定動作を長
時間行った結果としての伝達関数の推定値H ＾(k+1) の
みが必要な応用例や、また、伝達関数の推定結果そのも
のより未知系の出力の推定値y ＾(k) が必要な応用例
（例えば、時不変系の特性推定や音響エコーキャンセラ
など）においては、毎時刻推定伝達関数を計算する必要
はない。
図２に示した従来の伝達関数推定法では、伝達関数に相
当する推定伝達関数H ＾(k+1) を毎時刻計算し、これに
基づいて未知系の出力の推定値y ＾(k) を合成してき
た。これに対して、推定伝達関数の近似値Z(k)を用いた
本発明では、図５に示したように、推定伝達関数H ＾(k
+1) を計算しなくても、未知系の出力の推定値y ＾(k)
は合成できる。そして、y ＾(k) が合成できれば、上記
一連の動作を行うことができる。従って、推定動作を長
時間行った結果としての伝達関数の推定値H ＾(k+1) の
みが必要な応用例や、また、伝達関数の推定結果そのも
のより未知系の出力の推定値y ＾(k) が必要な応用例
（例えば、時不変系の特性推定や音響エコーキャンセラ
など）においては、毎時刻推定伝達関数を計算する必要
はない。
【００５７】このことより、各時刻で必要とされる総演
算量は、推定伝達関数計算部５６における演算量を除い
て、 19p＋2L (49) となる。従来技術の問題点の記述において述べたよう
に、フィルタのタップ数L=500 、射影の字数p=20とした
場合の従来技術の演算量は約12000 であった。しかし、
本発明の演算量は、式(49)より、1380となり、約1/8 の
演算量低減効果が得られたことがわかる。
算量は、推定伝達関数計算部５６における演算量を除い
て、 19p＋2L (49) となる。従来技術の問題点の記述において述べたよう
に、フィルタのタップ数L=500 、射影の字数p=20とした
場合の従来技術の演算量は約12000 であった。しかし、
本発明の演算量は、式(49)より、1380となり、約1/8 の
演算量低減効果が得られたことがわかる。
【００５８】以上説明したように本発明は、プレフィル
タ係数の回帰的合成、推定伝達関数ベクトルの近似値の
蓄積、ならびにプレフィルタリング係数の平滑化を行う
ことで、従来の射影法の有していた演算量を大幅に低減
するものである。以上の説明においては、ステップサイ
ズμがスカラーとして扱ってきたが、より一般的には、
対角行列Ａを用いてμＡとすることが考えられる。例え
ば、未知系のインパルス応答のエネルギーが指数減衰す
るものである場合には、Ａの要素を次式のように指数的
に減衰する数列にすると伝達関数の推定速度が速まるこ
とがある。
タ係数の回帰的合成、推定伝達関数ベクトルの近似値の
蓄積、ならびにプレフィルタリング係数の平滑化を行う
ことで、従来の射影法の有していた演算量を大幅に低減
するものである。以上の説明においては、ステップサイ
ズμがスカラーとして扱ってきたが、より一般的には、
対角行列Ａを用いてμＡとすることが考えられる。例え
ば、未知系のインパルス応答のエネルギーが指数減衰す
るものである場合には、Ａの要素を次式のように指数的
に減衰する数列にすると伝達関数の推定速度が速まるこ
とがある。
【００５９】 Ａ＝diag(α,αγ,αγ2,…,αγL-1),(0＜γ＜1) (50) この場合には、線形予測部への入力に比率γ1/2 の指数
をかけ、また、式(14)のＣp(k)を次式によって再定義す
る。 このように修正を行えば、以上説明してきた本発明を適
用することができる。
をかけ、また、式(14)のＣp(k)を次式によって再定義す
る。 このように修正を行えば、以上説明してきた本発明を適
用することができる。
【００６０】以下に本発明の伝達関数推定装置の適用例
を説明する。第１の適用例としては、音響機器の伝達関
数測定が上げられる。図６は、スピーカからマイクロホ
ンまでの伝達関数を測定する場合の系を示している。図
において１２１はスピーカ、１２２はマイクロホン、１
１は伝達関数推定装置を表している。マイクロホンの出
力信号y(k)は入力信号x(k)にスピーカの特性が付加され
たものである。スピーカ（音響パスとマイクロホンも含
める）を未知系と考えれば、図６は図１と同一の系とな
り、本発明の伝達関数推定装置１１を図６に示すように
スピーカ１２１の入力とマイクロホン１２２の出力に接
続すれば、スピーカの伝達関数H(k)をＦＩＲフィルタの
フィルタ係数H ＾(k) として推定することができる。測
定は、通常、部屋の音響特性の影響を受けないように無
響室で行われる。
を説明する。第１の適用例としては、音響機器の伝達関
数測定が上げられる。図６は、スピーカからマイクロホ
ンまでの伝達関数を測定する場合の系を示している。図
において１２１はスピーカ、１２２はマイクロホン、１
１は伝達関数推定装置を表している。マイクロホンの出
力信号y(k)は入力信号x(k)にスピーカの特性が付加され
たものである。スピーカ（音響パスとマイクロホンも含
める）を未知系と考えれば、図６は図１と同一の系とな
り、本発明の伝達関数推定装置１１を図６に示すように
スピーカ１２１の入力とマイクロホン１２２の出力に接
続すれば、スピーカの伝達関数H(k)をＦＩＲフィルタの
フィルタ係数H ＾(k) として推定することができる。測
定は、通常、部屋の音響特性の影響を受けないように無
響室で行われる。
【００６１】第２の応用例としては、テレビ会議システ
ムやビジュアルテレホンなどの拡声通話系においてハウ
リングやエコーを防止する音響エコーキャンセラがあげ
られる。図７に音響エコーキャンセラの説明図を示す。
図において１２１は受話スピーカ、１２２は送話マイク
ロホン、１２３は室内音響系を示し、２０は疑似エコー
生成部であり、図２においてこの発明が適用された推定
信号生成部２０と同じ構成である。この実施例の伝達関
数推定装置１１は音響エコーキャンセラとして動作す
る。スピーカとマイクロホンを用いたハンズフリー通話
系では、受話スピーカ１２１から出た相手側の音声が、
伝達関数H(k)の室内音響系１２３を経て送話スピーカ１
２２で受音される。受音された信号y(k)は、相手側に伝
送されて再生される。相手側においては、送話した音声
が戻ってきて再生されるので、この現象は音響エコーと
呼ばれ、快適な通話を妨害する。
ムやビジュアルテレホンなどの拡声通話系においてハウ
リングやエコーを防止する音響エコーキャンセラがあげ
られる。図７に音響エコーキャンセラの説明図を示す。
図において１２１は受話スピーカ、１２２は送話マイク
ロホン、１２３は室内音響系を示し、２０は疑似エコー
生成部であり、図２においてこの発明が適用された推定
信号生成部２０と同じ構成である。この実施例の伝達関
数推定装置１１は音響エコーキャンセラとして動作す
る。スピーカとマイクロホンを用いたハンズフリー通話
系では、受話スピーカ１２１から出た相手側の音声が、
伝達関数H(k)の室内音響系１２３を経て送話スピーカ１
２２で受音される。受音された信号y(k)は、相手側に伝
送されて再生される。相手側においては、送話した音声
が戻ってきて再生されるので、この現象は音響エコーと
呼ばれ、快適な通話を妨害する。
【００６２】音響エコーキャンセラ１１の疑似エコー生
成部２０はスピーカの特性をも含めた室内音響特性H
(k) を推定し、入力信号x(k)と推定された特性H ＾(k)
に基づいて疑似エコー信号y ＾(k) を合成し出力する。
差分器２４は受話信号y(k)から疑似エコー信号y ＾(k)
を差し引く。推定が良好に行われていれば誤差信号e(k)
のパワーが最小となるようにエコーキャンセラ１１が動
作し、y ＾(k) ≒y(k)となって音響エコーは大幅に低減
される。
成部２０はスピーカの特性をも含めた室内音響特性H
(k) を推定し、入力信号x(k)と推定された特性H ＾(k)
に基づいて疑似エコー信号y ＾(k) を合成し出力する。
差分器２４は受話信号y(k)から疑似エコー信号y ＾(k)
を差し引く。推定が良好に行われていれば誤差信号e(k)
のパワーが最小となるようにエコーキャンセラ１１が動
作し、y ＾(k) ≒y(k)となって音響エコーは大幅に低減
される。
【００６３】図７の系を図２と比較すれば、本発明の伝
達関数推定装置が音響エコーキャンセラとして直接適応
できることがわかる。図７における室内音響系１２３が
図２における未知系１２に対応し、図７における疑似エ
コーy ＾(k) および送話信号e(k)が、それぞれ、図２に
おける畳み込み演算部２３の出力y(k)及び誤差信号e(k)
に対応する。
達関数推定装置が音響エコーキャンセラとして直接適応
できることがわかる。図７における室内音響系１２３が
図２における未知系１２に対応し、図７における疑似エ
コーy ＾(k) および送話信号e(k)が、それぞれ、図２に
おける畳み込み演算部２３の出力y(k)及び誤差信号e(k)
に対応する。
【００６４】第３の適用例は騒音制御である。図８に騒
音制御の原理図を示す。図において１３１は騒音源、１
３２は伝達関数H(k)で表された騒音伝達経路、１３３は
観測用マイクロホン、１３４は騒音モニタマイクロホ
ン、２０は疑似騒音生成部、１３６は逆相器、１３７は
スピーカを表す。騒音制御の目的は騒音源１３１から伝
達特性H(k)を経て観測される騒音を、スピーカ１３７か
ら負の疑似騒音（音圧がy(k)と表されるときに、-y(k)
と表される音圧を、y(k)に対する負の音と呼ぶ）を出す
ことにより消去しようとするものである。
音制御の原理図を示す。図において１３１は騒音源、１
３２は伝達関数H(k)で表された騒音伝達経路、１３３は
観測用マイクロホン、１３４は騒音モニタマイクロホ
ン、２０は疑似騒音生成部、１３６は逆相器、１３７は
スピーカを表す。騒音制御の目的は騒音源１３１から伝
達特性H(k)を経て観測される騒音を、スピーカ１３７か
ら負の疑似騒音（音圧がy(k)と表されるときに、-y(k)
と表される音圧を、y(k)に対する負の音と呼ぶ）を出す
ことにより消去しようとするものである。
【００６５】この目的を達成するためには、図２におけ
るこの発明が適用された推定信号生成部２０と同じ構造
の疑似騒音生成部２０において騒音伝達経路１３２の伝
達特性H(k)を推定する。即ち、騒音源１３１の付近でモ
ニタマイクロホン１３４を用いて騒音を検出し、入力信
号x(k)として疑似騒音生成部２０に与える。それによっ
て観測地点（マイクロホン１３３）における騒音信号y
(k)の推定値y ＾(k) を生成する。疑似騒音y ＾(k) を
逆相器１３６で正負を反転して信号-y＾(k) とする。こ
こで、簡単のために、スピーカの特性は無視できるもの
と考えれば、観測点のマイクロホン１３３においてスピ
ーカ１３７からの信号-y＾(k) は騒音信号y(k)と音圧が
合成され、マイクロホン１３３の出力に誤差信号e(k)が
得られる。このとき、騒音伝達経路の推定が良好に行わ
れていれば、信号y ＾(k) は騒音源からの騒音y(k)と類
似したものとなり、スピーカから合成される音圧-y＾
(k) によって騒音y(k)はキャンセルされる。
るこの発明が適用された推定信号生成部２０と同じ構造
の疑似騒音生成部２０において騒音伝達経路１３２の伝
達特性H(k)を推定する。即ち、騒音源１３１の付近でモ
ニタマイクロホン１３４を用いて騒音を検出し、入力信
号x(k)として疑似騒音生成部２０に与える。それによっ
て観測地点（マイクロホン１３３）における騒音信号y
(k)の推定値y ＾(k) を生成する。疑似騒音y ＾(k) を
逆相器１３６で正負を反転して信号-y＾(k) とする。こ
こで、簡単のために、スピーカの特性は無視できるもの
と考えれば、観測点のマイクロホン１３３においてスピ
ーカ１３７からの信号-y＾(k) は騒音信号y(k)と音圧が
合成され、マイクロホン１３３の出力に誤差信号e(k)が
得られる。このとき、騒音伝達経路の推定が良好に行わ
れていれば、信号y ＾(k) は騒音源からの騒音y(k)と類
似したものとなり、スピーカから合成される音圧-y＾
(k) によって騒音y(k)はキャンセルされる。
【００６６】図８の系を図２と比較すれば、マイクロホ
ン１３３が差分器２４と対応し、従って、図８の系と図
２の系との相違点は、推定装置１１の外部で合成した誤
差信号を推定装置に供給するか、または、誤差信号を推
定装置１１の内部で計算するかという点に過ぎず、この
発明の原理は図８に示した騒音制御装置に適用可能であ
る。
ン１３３が差分器２４と対応し、従って、図８の系と図
２の系との相違点は、推定装置１１の外部で合成した誤
差信号を推定装置に供給するか、または、誤差信号を推
定装置１１の内部で計算するかという点に過ぎず、この
発明の原理は図８に示した騒音制御装置に適用可能であ
る。
【００６７】以上、３つの本発明の適用例を示した。図
６に示したスピーカの伝達特性の測定では、スピーカの
特性は時間的に変化しない。従って、測定の途中で測定
結果（伝達特性）は知る必要がなく、一定時間の測定動
作を終了した結果として伝達特性がわかればよい。図７
に示した音響エコーキャンセラの例では、例えば人の移
動やドアの開閉などによって、室内伝達特性は時間的に
変化する。しかし、図７から分かるように、室内伝達特
性（未知系）H(k)の出力y(k)の推定値y ＾(k)が得られ
れば音響エコーキャンセラは目的を達成する。従って、
音響エコーキャンセラにおいては、室内伝達特性そのも
のの推定値H ＾(k) は不用である。図８に示した騒音制
御の場合も、音響エコーキャンセラと同様に、未知系の
出力推定値y ＾(k) が得られれば目的を達成し、室内伝
達特性そのものの推定値は不用である。
６に示したスピーカの伝達特性の測定では、スピーカの
特性は時間的に変化しない。従って、測定の途中で測定
結果（伝達特性）は知る必要がなく、一定時間の測定動
作を終了した結果として伝達特性がわかればよい。図７
に示した音響エコーキャンセラの例では、例えば人の移
動やドアの開閉などによって、室内伝達特性は時間的に
変化する。しかし、図７から分かるように、室内伝達特
性（未知系）H(k)の出力y(k)の推定値y ＾(k)が得られ
れば音響エコーキャンセラは目的を達成する。従って、
音響エコーキャンセラにおいては、室内伝達特性そのも
のの推定値H ＾(k) は不用である。図８に示した騒音制
御の場合も、音響エコーキャンセラと同様に、未知系の
出力推定値y ＾(k) が得られれば目的を達成し、室内伝
達特性そのものの推定値は不用である。
【００６８】以上のことから、これらの適用例において
は、本発明における伝達関数推定装置を適用でき、従来
装置と比較して演算量を大幅に低減できる。また、本発
明は、図２に示したような系構成において誤差信号e(k)
のパワーを最小とするという基本機能を有している。従
って、解決すべき問題が、図２に示した誤差信号のパワ
ーの最小化問題としてモデル化できる全ての応用例に対
して、本発明を適用することが可能である。
は、本発明における伝達関数推定装置を適用でき、従来
装置と比較して演算量を大幅に低減できる。また、本発
明は、図２に示したような系構成において誤差信号e(k)
のパワーを最小とするという基本機能を有している。従
って、解決すべき問題が、図２に示した誤差信号のパワ
ーの最小化問題としてモデル化できる全ての応用例に対
して、本発明を適用することが可能である。
【００６９】上記実施例において、式(50)の対角行列Ａ
の要素は部屋の残響と同じ減衰をするように選ぶと効果
的である。最後に、本発明を用いて構成した音響エコー
キャンセラの実験結果を示す。図９は学習曲線と呼ばれ
るもので、縦軸はエコー消去量、横軸は時間を表してい
る。時間が経過して室内伝達特性の推定が進行すると共
に、エコー消去量は増加する。図において、曲線２１１
は射影の次数p=2 の場合、曲線２１２は射影の次数p=8
の場合、曲線２１３は射影の次数p=32の場合の学習曲線
を表している。図より、ｐが大きいほど、収束速度が早
い（エコー消去量が短時間で増加する）ことが理解でき
る。
の要素は部屋の残響と同じ減衰をするように選ぶと効果
的である。最後に、本発明を用いて構成した音響エコー
キャンセラの実験結果を示す。図９は学習曲線と呼ばれ
るもので、縦軸はエコー消去量、横軸は時間を表してい
る。時間が経過して室内伝達特性の推定が進行すると共
に、エコー消去量は増加する。図において、曲線２１１
は射影の次数p=2 の場合、曲線２１２は射影の次数p=8
の場合、曲線２１３は射影の次数p=32の場合の学習曲線
を表している。図より、ｐが大きいほど、収束速度が早
い（エコー消去量が短時間で増加する）ことが理解でき
る。
【００７０】次に、本発明の効果を表す図を示す。図１
０は、射影の次数ｐと必要な演算量の関係を示す図で、
横軸は射影の次数p を、縦軸は積和演算（和演算も含
む）の回数を表す。図において、曲線２２１は従来法の
演算量を、曲線２２２は本発明を用いた場合の演算量を
示した。但し、伝達関数の推定に用いるＦＩＲフィルタ
のタップ数L=500 とした。この図より、特にｐを大きく
した場合において、本発明は従来法と比べて演算量が大
幅に低減していることが分かる。
０は、射影の次数ｐと必要な演算量の関係を示す図で、
横軸は射影の次数p を、縦軸は積和演算（和演算も含
む）の回数を表す。図において、曲線２２１は従来法の
演算量を、曲線２２２は本発明を用いた場合の演算量を
示した。但し、伝達関数の推定に用いるＦＩＲフィルタ
のタップ数L=500 とした。この図より、特にｐを大きく
した場合において、本発明は従来法と比べて演算量が大
幅に低減していることが分かる。
【００７１】以上説明したように、本発明を利用すれ
ば、射影法による未知系の伝達関数の推定または、未知
系の出力の推定に要する演算量を大幅に低減することが
できる。具体的には、未知系を表現するＦＩＲフィルタ
のタップ長をＬ、射影の次数をｐとすると、従来技術で
は p3/6+(p+1)L、即ちp3に比例する積和演算量が必要で
あったのが、本発明では19p+2Lの演算量にまで削減する
ことができた。
ば、射影法による未知系の伝達関数の推定または、未知
系の出力の推定に要する演算量を大幅に低減することが
できる。具体的には、未知系を表現するＦＩＲフィルタ
のタップ長をＬ、射影の次数をｐとすると、従来技術で
は p3/6+(p+1)L、即ちp3に比例する積和演算量が必要で
あったのが、本発明では19p+2Lの演算量にまで削減する
ことができた。
【００７２】このような演算量の低減は、それに比較し
たハードウェア規模の低減が可能なので、装置の小規模
化、価格の低減化に寄与が大きい。また同一のハードウ
ェア規模とするなら、タップ長Ｌやｐを大きく選ぶこと
ができ、推定動作の高速化や推定精度の向上を得ること
ができる。またこの発明をコンピュータで実施する場
合、演算処理時間を大幅に削減することができる。
たハードウェア規模の低減が可能なので、装置の小規模
化、価格の低減化に寄与が大きい。また同一のハードウ
ェア規模とするなら、タップ長Ｌやｐを大きく選ぶこと
ができ、推定動作の高速化や推定精度の向上を得ること
ができる。またこの発明をコンピュータで実施する場
合、演算処理時間を大幅に削減することができる。
【図１】伝達関数推定を行う一般的構成を示すブロック
図。
図。
【図２】図１における伝達関数推定部１１の処理機能構
成を示すブロック図。
成を示すブロック図。
【図３】図２におけるこの発明を適用する推定伝達関数
修正ベクトル計算部２１の処理機能構成を示すブロック
図。
修正ベクトル計算部２１の処理機能構成を示すブロック
図。
【図４】この発明の第１の観点に基づく図３におけるプ
レフィルタ係数ベクトル計算部３１の処理機能構成を示
すブロック図。
レフィルタ係数ベクトル計算部３１の処理機能構成を示
すブロック図。
【図５】この発明の第２の観点に基づく推定信号生成処
理を示すブロック図。
理を示すブロック図。
【図６】この発明をスピーカの伝達関数測定に適用した
場合の構成ブロック図。
場合の構成ブロック図。
【図７】この発明を音響エコーキャンセラに適用した場
合の構成ブロック図。
合の構成ブロック図。
【図８】この発明を騒音制御に適用した場合の構成ブロ
ック図。
ック図。
【図９】この発明を音響エコーキャンセラに適用した場
合のエコー消去特性を示すグラフ。
合のエコー消去特性を示すグラフ。
【図１０】この発明による伝達関数推定の演算量と射影
次数の関係を従来例と比較して示すグラフ。
次数の関係を従来例と比較して示すグラフ。
───────────────────────────────────────────────────── フロントページの続き (72)発明者 小島 順治 東京都千代田区内幸町１丁目１番６号 日本電信電話株式会社内 (72)発明者 牧野 昭二 東京都千代田区内幸町１丁目１番６号 日本電信電話株式会社内 (58)調査した分野(Int.Cl.7，ＤＢ名) G10K 11/178 G01H 3/00 G05B 13/02 G06F 17/12 H03H 21/00
Claims (12)
- 【請求項１】 離散的時刻ｋにおける未知系の入力信号
x(k)と出力y(k)とからその未知系の伝達関数を推定し、
その推定した伝達関数H ＾(k) を持つ疑似系の出力信号
y ＾(k) の上記未知系の出力信号y(k)に対する誤差信号
e(k)=y(k)-y＾(k) を算出し、その誤差信号ベクトルをE
(k)、入力信号x(k)の共分散行列をＣp(k) 、プレフィル
タ係数ベクトルをG(k)とするとき、p 元連立１次方程式 Ｃp(k)G(k) ＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求め、
上記入力信号x(k)のベクトルをX(k)とし、推定伝達関数
修正ベクトルをδH ＾(k) とし、上記プレフィルタ係数
ベクトルG(k)を使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の上記推定伝達関数に対する上記修正
ベクトルδH ＾(k) を計算し、その推定伝達関数修正ベ
クトルδH ＾(k) と、予め決めた修正ステップサイズμ
とを使って次式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正することを繰り返し上記誤差信号e(k)がゼロに近づく
ようにする射影アルゴリズムによる適応的伝達関数推定
方法において、 上記ｐ元連立１次方程式を解いて上記プレフィルタ係数
ベクトルG(k)を求めるステップは、(a) 上記入力信号x
(k)の前向き線形予測係数ベクトルF(k)とその２乗予測
誤差和f(k)と、後向き線形予測係数ベクトルB(k)とその
２乗予測誤差和b(k)を算出し、(b) プレフィルタ誘導係
数ベクトルをD(k)と表すとき、次の第１及び第２式から
成る上記プレフィルタ係数ベクトルG(k)についての回帰
式 【数１】 により上記プレフィルタ係数ベクトルG(k)を求めるステ
ップである。 - 【請求項２】 離散的時刻ｋにおける未知系の入力信号
x(k)と出力y(k)とからその未知系の伝達関数を推定し、
その推定した伝達関数H ＾(k) を持つ疑似系の出力信号
y ＾(k) の上記未知系の出力信号y(k)に対する誤差信号
e(k)=y(k)-y＾(k) を算出し、その誤差信号ベクトルをE
(k)、入力信号x(k)の共分散行列をＣp(k)、プレフィル
タ係数ベクトルをG(k)とするとき、p 元連立１次方程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求め、
上記プレフィルタ係数ベクトルG(k)を使って次式 δH ＾(k) ＝[x(k),x(k-1),…,x(k-p+1)]G(k) により上記疑似系の推定伝達関数修正ベクトルδH ＾
(k) を計算し、その推定伝達関数修正ベクトルδH ＾
(k) と、予め決めた修正ステップサイズμとを使って次
式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正することを繰り返し上記誤差信号e(k)がゼロに近づく
ようにする射影アルゴリズムによる適応的伝達関数推定
方法において、 上記推定伝達関数修正ベクトルを計算する代わりに、上
記プレフィルタ係数ベクトルG(k)の要素であるプレフィ
ルタ係数g(k)を次式 si(k) ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μg1(k) for i = 1 により平滑化して平滑化係数si(k) を得て、上記推定伝
達関数ベクトルを上記修正ベクトルで修正する代わり
に、上記平滑化係数s(k)を使って次式 Z(k+1)＝Z(k)＋sp(k)X(k-p+1) により近似推定伝達関数Z(k+1)を求め、上記近似推定伝
達関数Z(k+1)と上記入力信号x(k)の畳み込み演算X(k)TZ
(k) を行い、 上記平滑化係数s(k)のベクトル、上記入力信号x(k)のベ
クトル及び上記入力信号の相関ベクトルをそれぞれ Sp-1(k-1)＝[s1(k-1),s2(k-1),…,sp-1(k-1)]T Rp-1(k)＝[X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T として次式の内積 Sp-1(k-1)TRp-1(k) を計算し、 上記畳み込み演算結果X(k)TZ(k) と上記内積演算結果の
和を上記推定信号y ＾(k) として出力する。 - 【請求項３】 請求項１に記載の適応的伝達関数推定方
法において、 上記プレフィルタ係数ベクトルG(k)の要素であるプレフ
ィルタ係数g(k)を次式 si(k) ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μg1(k) for i = 1 により平滑化して平滑化係数si(k) を得て、 上記平滑化係数s(k)を使って次式 Z(k+1) ＝ Z(k)＋sp(k)X(k-p+1) により近似推定伝達関数Z(k+1)を求め、 上記近似推定伝達関数Z(k+1)と上記入力信号x(k)の畳み
込み演算X(k)TZ(k) を行い、 上記平滑化係数s(k)のベクトル、上記入力信号x(k)のベ
クトル及び上記入力信号の相関ベクトルをそれぞれ Sp-1(k-1)＝[s1(k-1),s2(k-1),…,sp-1(k-1)]T Rp-1(k)＝[X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T として次式の内積 Sp-1(k-1)TRp-1(k) を計算し、 上記畳み込み演算結果X(k)TZ(k) と上記内積演算結果の
和を上記推定信号y ＾(k) として出力する。 - 【請求項４】 請求項１または３に記載の適応的伝達関
数推定方法において、上記プレフィルタ係数ベクトルG
(k-1)の最後の要素をgp(k-1) とするとき、上記第２式
の代わりに上記第２式を変形した次式 を用いて上記プレフィルタ係数ベクトルG(k)を計算す
る。 - 【請求項５】 請求項１または３に記載の適応的伝達関
数推定方法において、予め決めた２つの非負数をδ
F(k)，δB(k-1)とし、上記第１及び第２式の代わりにそ
れらを変形した次式 【数２】 により上記プレフィルタ係数ベクトルG(k)を計算する。 - 【請求項６】 離散的時刻ｋにおける未知系の入力信号
x(k)と出力y(k)とから推定した伝達関数H ＾(k) を持つ
畳み込み演算手段により形成された疑似系と、上記入力
信号に対する上記疑似系の出力信号y (k) と上記未知系
の出力との誤差信号e(k)=y(k)-y ＾(k) のベクトルをE
(k)、上記入力信号x(k)の共分散行列をＣp(k)、プレフ
ィルタ係数ベクトルをG(k)と表すとき、p 元連立１次方
程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求める
プレフィルタ係数ベクトル計算手段と、上記入力信号x
(k)のベクトルをX(k)、推定伝達関数修正ベクトルをδH
＾(k) と表し、上記プレフィルタ係数ベクトルG(k)を
使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の上記推定伝達関数に対する上記修正
ベクトルδH ＾(k) を計算するプレフィルタリング手段
と、その推定伝達関数修正ベクトルδH ＾(k) と、予め
決めた修正ステップサイズμとを使って次式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正する推定伝達関数修正手段とを含み、上記推定伝達関
数ベクトルの修正を繰り返し上記誤差信号e(k)がゼロに
近づくようにする射影アルゴリズムによる適応的伝達関
数推定装置において、上記プレフィルタ係数ベクトル計
算手段は： 上記入力信号x(k)の前向き線形予測係数ベクトルF(k)と
その２乗予測誤差f(k)と、後向き線形予測係数ベクトル
B(k)とその２乗予測b(K)を算出する線形予測手段と、 プレフィルタ誘導係数ベクトルをD(k)とするとき、次の
第１及び第２式から成る上記プレフィルタ係数ベクトル
についての回帰式 【数３】 の上記第１式によりプレフィルタ係数ベクトルG(k)を修
正するプレフィルタ係数ベクトル修正手段と、 上記第２式によりプレフィルタ誘導係数ベクトルD(k)を
修正するプレフィルタ誘導係数ベクトル修正手段とを含
む。 - 【請求項７】 離散的時刻ｋにおける未知系の入力信号
x(k)と出力y(k)とから推定した伝達関数H ＾(k) を持つ
畳み込み演算手段により形成された疑似系と、上記入力
信号に対する上記疑似系の出力信号y ＾(k) と上記未知
系の出力との誤差信号e(k)=y(k)-y ＾(k) のベクトルを
G(k)、入力信号x(k)の共分散行列をＣp(k)、プレフィル
タ係数ベクトルをG(k)と表すとき、p 元連立１次方程式 Ｃp(k)G(k)＝ E(k) を解いて、上記プレフィルタ係数ベクトルG(k)を求める
プレフィルタ係数ベクトル計算手段と、上記入力信号x
(k)のベクトルをX(k)、推定伝達関数修正ベクトルをδH
＾(k) と表し、上記プレフィルタ係数ベクトルG(k)を
使って次式 δH ＾(k) ＝[X(k),X(k-1),…,X(k-p+1)]G(k) により上記疑似系の上記推定伝達関数に対する上記修正
ベクトルδH ＾(k) を計算するプレフィルタリング手段
と、その推定伝達関数修正ベクトルδH ＾(k) と、予め
決めた修正ステップサイズμとを使って次式 H ＾(k+1) ＝H ＾(k) ＋μδH ＾(k) により上記疑似系の推定伝達関数ベクトルH ＾(k) を修
正する推定伝達関数修正手段とを含み、上記推定伝達関
数ベクトルの修正を繰り返し上記誤差信号e(k)がゼロに
近づくようにする射影アルゴリズムによる適応的伝達関
数推定装置において、 上記プレフィルタリング手段は上記修正ベクトルを計算
する代わりに、上記プレフィルタ係数ベクトルG(k)の要
素であるプレフィルタ係数g(k)を次式 si(k) ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μg1(k) for i = 1 により平滑化して平滑化係数si(k) として出力するプレ
フィルタ係数平滑手段であり、 上記推定伝達関数修正手段は： 上記修正ベクトルにより上記推定伝達関数ベクトルを修
正する代わりに、上記平滑化係数s(k)を使って次式 Z(k+1) ＝ Z(k)＋sp(k)X(k-p+1) により近似推定伝達関数Z(k+1)を求める近似推定伝達関
数計算手段と、 上記近似推定伝達関数Z(k+1)と上記入力信号x(k)の畳み
込み演算X(k)TZ(k) を行う畳み込み計算手段と、 上記平滑化係数s(k)のベクトル、上記入力信号x(k)のベ
クトル及び上記入力信号の相関ベクトルをそれぞれ Sp-1(k-1)＝[s1(k-1),s2(k-1),…,sp-1(k-1)]T Rp-1(k)＝[X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T として次式の内積 Sp-1(k-1)TRp-1(k) を計算する内積計算手段と、 上記畳み込み演算結果X(k)TZ(k) と上記内積演算結果の
和を上記推定信号y ＾(k) として出力する加算手段とを
含む。 - 【請求項８】 請求項７に記載の適応的伝達関数推定装
置において、上記畳み込み手段のタップ数をＬとする
と、上記内積計算手段は上記入力信号ベクトルからその
上記相関ベクトルを次式 Rp-1(k)＝Rp-1(k-1)−x(k-L)Xp-1(k-L)＋x(k)Xp-1(k) により計算する相関計算手段を含む。 - 【請求項９】 請求項６に記載の適応的伝達関数推定装
置において、上記プレフィルタリング手段は上記プレフ
ィルタ係数ベクトルG(k)の要素であるプレフィルタ係数
g(k)を次式 si(k) ＝ si-1(k-1)＋μgi(k) for 2≦i≦p ＝ μg1(k) for i = 1 により平滑化して平滑化係数si(k) として出力するプレ
フィルタ係数平滑手段であり、 上記推定伝達関数修正手段は： 上記平滑化係数s(k)を使って次式 Z(k+1) ＝ Z(k)＋sp(k)X(k-p+1) により近似推定伝達関数Z(k+1)を求める近似推定伝達関
数計算手段と、 上記近似推定伝達関数Z(k+1)と上記入力信号x(k)の畳み
込み演算X(k)TZ(k) を行う畳み込み計算手段と、 上記平滑化係数s(k)のベクトル、上記入力信号x(k)のベ
クトル及び上記入力信号の相関ベクトルをそれぞれ Sp-1(k-1)＝[s1(k-1),s2(k-1),…,sp-1(k-1)]T Rp-1(k)＝[X(k)TX(k-1),X(k)TX(k-2),…,X(k)TX(k-p+1)]T として次式の内積 Sp-1(k-1)TRp-1(k) を計算する内積計算手段と、 上記畳み込み演算結果X(k)TZ(k) と上記内積演算結果の
和を上記推定信号y ＾(k) として出力する加算手段とを
含む。 - 【請求項１０】 請求項６または９に記載の適応的伝達
関数推定装置において、上記プレフィルタ誘導係数ベク
トル修正手段は上記プレフィルタ係数ベクトルG(k-1)の
最後の要素をgp(k-1) とするとき、上記第２式の代わり
に上記第２式を変形した次式 を用いて上記プレフィルタ係数ベクトルG(k)を計算する
手段である。 - 【請求項１１】 請求項６または９に記載の適応的伝達
関数推定装置において、予め決めた２つの非負数をδ
F(k)，δB(k-1)とし、上記プレフィルタ係数ベクトル修
正手段と上記プレフィルタ誘導係数ベクトル修正手段は
それぞれ上記第１及び第２式の代わりにそれらを変形し
た次式 【数４】 により上記プレフィルタ係数ベクトルG(k)を計算する手
段である。 - 【請求項１２】 請求項６または９に記載の適応的伝達
関数推定装置において、上記未知系の出力信号y(k)と上
記疑似系の出力である上記推定信号y ＾(k)の差y(k)-y
＾(k) を上記誤差信号e(k)として出力する差分手段が設
けられている。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP17541794A JP3303898B2 (ja) | 1993-07-27 | 1994-07-27 | 適応的伝達関数推定方法及びそれを使った推定装置 |
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP18474293 | 1993-07-27 | ||
| JP5-184742 | 1993-07-27 | ||
| JP17541794A JP3303898B2 (ja) | 1993-07-27 | 1994-07-27 | 適応的伝達関数推定方法及びそれを使った推定装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH0792980A JPH0792980A (ja) | 1995-04-07 |
| JP3303898B2 true JP3303898B2 (ja) | 2002-07-22 |
Family
ID=26496705
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP17541794A Expired - Lifetime JP3303898B2 (ja) | 1993-07-27 | 1994-07-27 | 適応的伝達関数推定方法及びそれを使った推定装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3303898B2 (ja) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5860383B2 (ja) * | 2012-11-14 | 2016-02-16 | 日本電信電話株式会社 | 伝達系パラメータ推定装置、伝達系パラメータ推定方法、伝達系パラメータ推定プログラム |
-
1994
- 1994-07-27 JP JP17541794A patent/JP3303898B2/ja not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JPH0792980A (ja) | 1995-04-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Douglas | Introduction to adaptive filters | |
| JP5284475B2 (ja) | 前白色化を伴うｌｍｓアルゴリズムによって適応させられる適応フィルタの更新済みフィルタ係数を決定する方法 | |
| Kuo et al. | Nonlinear adaptive bilinear filters for active noise control systems | |
| US5371789A (en) | Multi-channel echo cancellation with adaptive filters having selectable coefficient vectors | |
| US5774562A (en) | Method and apparatus for dereverberation | |
| JP3008763B2 (ja) | 適応フィルタによるシステム同定の方法および装置 | |
| Albu et al. | Pseudo-affine projection algorithms for multichannel active noise control | |
| US6381272B1 (en) | Multi-channel adaptive filtering | |
| JP2004349806A (ja) | 多チャネル音響エコー消去方法、その装置、そのプログラム及びその記録媒体 | |
| EP0637803B1 (en) | Method and device for adaptively estimating a transfer function of an unknown system | |
| Van Vaerenbergh et al. | A split kernel adaptive filtering architecture for nonlinear acoustic echo cancellation | |
| JP3303898B2 (ja) | 適応的伝達関数推定方法及びそれを使った推定装置 | |
| Ciochina et al. | An optimized proportionate adaptive algorithm for sparse system identification | |
| JP3673727B2 (ja) | 反響消去方法、その装置、そのプログラム及びその記録媒体 | |
| Fan et al. | Effective improvement of under-modeling frequency-domain Kalman filter | |
| CN111883155B (zh) | 回声消除方法、装置及存储介质 | |
| JP3616341B2 (ja) | 多チャネルエコーキャンセル方法、その装置、そのプログラム及び記録媒体 | |
| JPH09261135A (ja) | 音響エコー消去装置 | |
| Hofmann et al. | Recent advances on LIP nonlinear filters and their applications: Efficient solutions and significance-aware filtering | |
| JP3475318B2 (ja) | 適応的制御方法 | |
| JP3147207B2 (ja) | 適応的未知系出力推定方法 | |
| EP1057259A1 (en) | Stable adaptive filter and method | |
| JP3452339B2 (ja) | 適応的制御方法 | |
| Kar et al. | An improved order estimation of MSF for stereophonic acoustic echo cancellation | |
| Hill et al. | An analysis of the exponentiated gradient descent algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090510 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090510 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100510 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100510 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110510 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120510 Year of fee payment: 10 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130510 Year of fee payment: 11 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140510 Year of fee payment: 12 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313531 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| EXPY | Cancellation because of completion of term |