JP3475318B2 - 適応的制御方法 - Google Patents
適応的制御方法Info
- Publication number
- JP3475318B2 JP3475318B2 JP5048497A JP5048497A JP3475318B2 JP 3475318 B2 JP3475318 B2 JP 3475318B2 JP 5048497 A JP5048497 A JP 5048497A JP 5048497 A JP5048497 A JP 5048497A JP 3475318 B2 JP3475318 B2 JP 3475318B2
- Authority
- JP
- Japan
- Prior art keywords
- vector
- calculation
- equation
- filter
- block
- 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 - Fee Related
Links
Landscapes
- Filters That Use Time-Delay Elements (AREA)
- Cable Transmission Systems, Equalization Of Radio And Reduction Of Echo (AREA)
Description
【0001】
【発明の属する技術分野】この発明は、音響エコーキャ
ンセラ、能動騒音制御等において、適応フィルタの推定
誤差のパワーを最小にする適応制御方法に関する。
ンセラ、能動騒音制御等において、適応フィルタの推定
誤差のパワーを最小にする適応制御方法に関する。
【0002】
【従来の技術】この節では、まず、適応的制御に用いら
れる、適応フィルタと適応アルゴリズムについて説明す
る。次に、その適応フィルタを修正するための適応アル
ゴリズムの中の射影法を説明し、さらに、射影法の演算
量低減法として、演算手順を工夫した高速射影法と適応
フィルタ修正を毎時刻行わないブロック処理に基づく射
影法について説明する。最後に、ブロック処理に基づく
射影法の問題点を述べる。 3.2.1 適応フィルタと適応アルゴリズム 図6Aは適応フィルタを用いた適応制御系のブロック図
を示す。図において、入力信号をx(k)、出力信号を
y(k)、所望信号をd(k)、誤差信号をe(k)で
表わす。ただし、この明細書では、時間の表現は離散時
間として整数kで時間を表すものとする。また、適応フ
ィルタはFIR形フィルタであって、時変係数{h
1 (k),h2 (k),…,hL (k)}を持つものと
する。ただし、Lはフィルタ係数の数である。また、フ
ィルタ係数はベクトルとして、 H(k)=[h1 (k),h2 (k),…,hL (k)]T (1) と表す。ただし、[ ]T はベクトルまたは行列の転置
を表す。適応フィルタは、誤差信号e(k)と入力信号
x(k)に基づいて、フィルタ係数H(k)を調整し
て、誤差信号e(k)のパワーを最小にする。この時の
係数修正の手順は適応アルゴリズムとして知られてい
る。
れる、適応フィルタと適応アルゴリズムについて説明す
る。次に、その適応フィルタを修正するための適応アル
ゴリズムの中の射影法を説明し、さらに、射影法の演算
量低減法として、演算手順を工夫した高速射影法と適応
フィルタ修正を毎時刻行わないブロック処理に基づく射
影法について説明する。最後に、ブロック処理に基づく
射影法の問題点を述べる。 3.2.1 適応フィルタと適応アルゴリズム 図6Aは適応フィルタを用いた適応制御系のブロック図
を示す。図において、入力信号をx(k)、出力信号を
y(k)、所望信号をd(k)、誤差信号をe(k)で
表わす。ただし、この明細書では、時間の表現は離散時
間として整数kで時間を表すものとする。また、適応フ
ィルタはFIR形フィルタであって、時変係数{h
1 (k),h2 (k),…,hL (k)}を持つものと
する。ただし、Lはフィルタ係数の数である。また、フ
ィルタ係数はベクトルとして、 H(k)=[h1 (k),h2 (k),…,hL (k)]T (1) と表す。ただし、[ ]T はベクトルまたは行列の転置
を表す。適応フィルタは、誤差信号e(k)と入力信号
x(k)に基づいて、フィルタ係数H(k)を調整し
て、誤差信号e(k)のパワーを最小にする。この時の
係数修正の手順は適応アルゴリズムとして知られてい
る。
【0003】次に、適応アルゴリズムについて説明す
る。適応アルゴリズムとは、例えば図6Aに示すように
入力信号x(k)がフィルタリング部11に供給され、
フィルタリング部11の出力信号y(k)と所望信号d
(k)との誤差信号e(k)が差回路13でとられ、誤
差信号e(k)と入力信号x(k)とがフィルタ係数更
新部12に入力され、フィルタリング部11のフィルタ
係数H(k+1)の更新がなされる。適応アルゴリズム
はこのような適応フィルタの内部の演算手順であって、
時刻kにおいて、フィルタ係数H(k)をフィルタ係数
更新部12で修正して、誤差信号e(k)のパワーがな
るべく小さくなるような係数H(k+1)を得る方法で
ある。その修正は、次式で表される。
る。適応アルゴリズムとは、例えば図6Aに示すように
入力信号x(k)がフィルタリング部11に供給され、
フィルタリング部11の出力信号y(k)と所望信号d
(k)との誤差信号e(k)が差回路13でとられ、誤
差信号e(k)と入力信号x(k)とがフィルタ係数更
新部12に入力され、フィルタリング部11のフィルタ
係数H(k+1)の更新がなされる。適応アルゴリズム
はこのような適応フィルタの内部の演算手順であって、
時刻kにおいて、フィルタ係数H(k)をフィルタ係数
更新部12で修正して、誤差信号e(k)のパワーがな
るべく小さくなるような係数H(k+1)を得る方法で
ある。その修正は、次式で表される。
【0004】
H(k+1)=H(k)+μΔH(k) (2.1)
ここで、ΔH(k)は修正ベクトル、μは修正の大きさ
を制御する量でステップサイズと呼ばれる。現在、最も
多く利用されている適応アルゴリズムとしては、LMS
法や学習同定法などが知られている。LMS法において
は、修正ベクトルは誤差信号e(k)と入力ベクトルx
(k)との積として、 ΔH(k)=e(k)XL (k) (2.2) と表され、修正式は、 H(k+1)=H(k)+μe(k)XL (k) (2.3) となる。ただし、 XL (k)=[x(k),x(k−1),…,x(k−L+1)]T (3.1) で、これを入力(信号)ベクトルまたは単に、入力と呼
ぶ。また、次の時刻k+1のフィルタ出力y(k+1)
は、図6Aのフィルタリング部11で次式による内積演
算により合成される。
を制御する量でステップサイズと呼ばれる。現在、最も
多く利用されている適応アルゴリズムとしては、LMS
法や学習同定法などが知られている。LMS法において
は、修正ベクトルは誤差信号e(k)と入力ベクトルx
(k)との積として、 ΔH(k)=e(k)XL (k) (2.2) と表され、修正式は、 H(k+1)=H(k)+μe(k)XL (k) (2.3) となる。ただし、 XL (k)=[x(k),x(k−1),…,x(k−L+1)]T (3.1) で、これを入力(信号)ベクトルまたは単に、入力と呼
ぶ。また、次の時刻k+1のフィルタ出力y(k+1)
は、図6Aのフィルタリング部11で次式による内積演
算により合成される。
【0005】
y(k+1)=XL (k+1)T H(k+1) (3.2)
式(2.3)に基づいてフィルタ係数を修正して、式
(3.2)に基づいて次の時刻の出力を合成するという
適応フィルタの動作に必要な積和演算量は2L(Lはフ
ィルタ係数の数)である。この2Lという演算量は、現
在知られている毎時刻フィルタ係数を修正する適応アル
ゴリズムのなかでは最小の演算量となっている。 3.2.2 射影法 このように、LMS法や学習同定法は演算量が少ないと
いう長所を持っている。しかし、その反面、有色信号や
周期性雑音などに対しては、収束速度が遅い(適切なフ
ィルタ係数を得るまでの時間がかかる)という欠点をも
っている。
(3.2)に基づいて次の時刻の出力を合成するという
適応フィルタの動作に必要な積和演算量は2L(Lはフ
ィルタ係数の数)である。この2Lという演算量は、現
在知られている毎時刻フィルタ係数を修正する適応アル
ゴリズムのなかでは最小の演算量となっている。 3.2.2 射影法 このように、LMS法や学習同定法は演算量が少ないと
いう長所を持っている。しかし、その反面、有色信号や
周期性雑音などに対しては、収束速度が遅い(適切なフ
ィルタ係数を得るまでの時間がかかる)という欠点をも
っている。
【0006】LMS法や学習同定法の欠点を克服する方
法として、1984年に射影法が提案されたが、射影法
では演算量が多くなるという欠点をもっていた。しかし
近年、この射影法の欠点を改善して、低演算量で実行す
る高速射影法が提案されている。ここではその概要を以
下に説明する。射影法において、誤差を小さくするため
には、フィルタ係数の修正は、次のp個の式 d(k)=XL (k)T H(k+1) d(k+1)=XL (k+1)T H(k+1) (4) … d(k−p+1)=XL (k−p+1)T H(k+1) を満たすように行われる。ただし、p(p<L)は射影
の次数と呼ばれる。この式(4)の意味は、修正を行っ
た後の係数H(k+1)を過去の入力XL (k−i+
1)T と畳み込んだ結果XL (k−i+1)T H(k+
1)が、過去の所望信号d(k−i+1)と一致する、
という関係が過去p個の入力に対して成立するように、
H(k+1)を決定するというものである。このよう
に、過去の入力に対して所望信号と同一の値を出力する
ようにフィルタ係数を定めれば、未来の時刻の入力XL
(k+1)に対する出力y(k+1)も所望信号d(k
+1)に近い値が出力できるものと考えられる。その結
果、未来の誤差信号e(k+1)=d(k+1)−y
(k+1)も小さな値となるものと考えられる。さて、
式(2.1)を式(4)に代入してΔH(k)について
解けば、 ΔH(k)=XL,p (k) [XL,p (k) T XL,p (k) ]-1Ep (k) (5) と表される。ただし、 XL,p (k) =[XL (k),XL (k-1),…, XL (k-p+1) ] (6) Ep (k) =[d(k),d(k-1),…,d(k-p+1) ]T −XL,p (k) T H(k) (7.1) である。Xの下付き添え字L,pはXがL行p列である
ことを示す。また、このEp (k)は、次式のように、
再帰的にも表すことができる。
法として、1984年に射影法が提案されたが、射影法
では演算量が多くなるという欠点をもっていた。しかし
近年、この射影法の欠点を改善して、低演算量で実行す
る高速射影法が提案されている。ここではその概要を以
下に説明する。射影法において、誤差を小さくするため
には、フィルタ係数の修正は、次のp個の式 d(k)=XL (k)T H(k+1) d(k+1)=XL (k+1)T H(k+1) (4) … d(k−p+1)=XL (k−p+1)T H(k+1) を満たすように行われる。ただし、p(p<L)は射影
の次数と呼ばれる。この式(4)の意味は、修正を行っ
た後の係数H(k+1)を過去の入力XL (k−i+
1)T と畳み込んだ結果XL (k−i+1)T H(k+
1)が、過去の所望信号d(k−i+1)と一致する、
という関係が過去p個の入力に対して成立するように、
H(k+1)を決定するというものである。このよう
に、過去の入力に対して所望信号と同一の値を出力する
ようにフィルタ係数を定めれば、未来の時刻の入力XL
(k+1)に対する出力y(k+1)も所望信号d(k
+1)に近い値が出力できるものと考えられる。その結
果、未来の誤差信号e(k+1)=d(k+1)−y
(k+1)も小さな値となるものと考えられる。さて、
式(2.1)を式(4)に代入してΔH(k)について
解けば、 ΔH(k)=XL,p (k) [XL,p (k) T XL,p (k) ]-1Ep (k) (5) と表される。ただし、 XL,p (k) =[XL (k),XL (k-1),…, XL (k-p+1) ] (6) Ep (k) =[d(k),d(k-1),…,d(k-p+1) ]T −XL,p (k) T H(k) (7.1) である。Xの下付き添え字L,pはXがL行p列である
ことを示す。また、このEp (k)は、次式のように、
再帰的にも表すことができる。
【0007】
ただし、
Ep-1 (k−1)=[e(k−1),(1−μ)e(k−2),…,
(1−μ)p-2 e(k−p+1)]T (7.3)
である。
【0008】式(5),(2.1)より、入力X
L (k)と誤差信号e(k)を用いてフィルタ係数を修
正して、係数H(k+1)を求め、これを用いて時刻k
+1の出力y(k+1)を合成する手順は、次式のよう
に表すことができる。 G(k)=[XL,p (k) T XL,p (k) ]-1Ep (k) (8) H(k+1)=H(k)+μXLp(k)G(k) (9) y(k+1)=XL (k+1)T H(k+1) (10) ただし、G(k)は式(5)の後半を表し、プレフィル
タ係数と呼ばれる。この式(8),(9),(10)が
従来の射影法の原理的手順を表す式であり、これを機能
ブロック図で表したものを図6Bに示す。図6Bにおい
てプレフィルタ係数計算部21で入力信号x(k)と誤
差信号e(k)とから式(8)の演算を行い、得られた
プレフィルタ係数G(k)と入力信号x(k)とによ
り、フィルタ係数更新部22で式(9)の演算を行いフ
ィルタ係数を修正し、その修正されたフィルタ係数と入
力信号とによりフィルタリング部11で式(10)の演
算を実行する。
L (k)と誤差信号e(k)を用いてフィルタ係数を修
正して、係数H(k+1)を求め、これを用いて時刻k
+1の出力y(k+1)を合成する手順は、次式のよう
に表すことができる。 G(k)=[XL,p (k) T XL,p (k) ]-1Ep (k) (8) H(k+1)=H(k)+μXLp(k)G(k) (9) y(k+1)=XL (k+1)T H(k+1) (10) ただし、G(k)は式(5)の後半を表し、プレフィル
タ係数と呼ばれる。この式(8),(9),(10)が
従来の射影法の原理的手順を表す式であり、これを機能
ブロック図で表したものを図6Bに示す。図6Bにおい
てプレフィルタ係数計算部21で入力信号x(k)と誤
差信号e(k)とから式(8)の演算を行い、得られた
プレフィルタ係数G(k)と入力信号x(k)とによ
り、フィルタ係数更新部22で式(9)の演算を行いフ
ィルタ係数を修正し、その修正されたフィルタ係数と入
力信号とによりフィルタリング部11で式(10)の演
算を実行する。
【0009】ここで、この系の演算量を考えてみる。ほ
とんどのDSPでは、加算と積和演算がそれぞれ1クロ
ックで実行され、また、加算と積和演算が計算の最も大
きな部分を占めるので、ここでは、演算量として加算と
積和演算の数として評価する。まず、プレフィルタ係数
計算部21で行う式(8)の演算は、p行p列の行列X
L,P (k)T XL,P (k)の逆行列の計算を含んでいる
ので、pの3乗に比例する積和演算量(以降、O
(p3 )と表す)が必要となる。これは、pが大きくな
ると、たいへん大きな演算量となる。次に、フィルタ係
数更新部22で行う式(9)の演算はXL,p (k)がL
行p列の行列、G(k)がp次のベクトルであることよ
り、pL回の積和演算が必要である。更に、フィルタ実
行部11で行う式(10)の演算はL次のベクトルの内
積であるので、L回の積和演算が必要である。従って、
合計(p+1)L+O(p3 )の積和演算量が必要とな
る。例えば、p=20、L=1000の場合には、2
1,000回以上もの積和演算が必要となる。
とんどのDSPでは、加算と積和演算がそれぞれ1クロ
ックで実行され、また、加算と積和演算が計算の最も大
きな部分を占めるので、ここでは、演算量として加算と
積和演算の数として評価する。まず、プレフィルタ係数
計算部21で行う式(8)の演算は、p行p列の行列X
L,P (k)T XL,P (k)の逆行列の計算を含んでいる
ので、pの3乗に比例する積和演算量(以降、O
(p3 )と表す)が必要となる。これは、pが大きくな
ると、たいへん大きな演算量となる。次に、フィルタ係
数更新部22で行う式(9)の演算はXL,p (k)がL
行p列の行列、G(k)がp次のベクトルであることよ
り、pL回の積和演算が必要である。更に、フィルタ実
行部11で行う式(10)の演算はL次のベクトルの内
積であるので、L回の積和演算が必要である。従って、
合計(p+1)L+O(p3 )の積和演算量が必要とな
る。例えば、p=20、L=1000の場合には、2
1,000回以上もの積和演算が必要となる。
【0010】次に、射影法の演算量低減法として、演算
を工夫した高速射影法と適応フィルタ更新を毎時刻には
行わず、あるブロック間隔で行うブロック処理に基づく
射影法について説明する。 3.2.3 高速射影法 これに対して近年提案された高速射影法では、演算の工
夫により積和演算量を2L+20pにまで低減すること
が可能である。以下、高速射影法を簡単に説明する。
を工夫した高速射影法と適応フィルタ更新を毎時刻には
行わず、あるブロック間隔で行うブロック処理に基づく
射影法について説明する。 3.2.3 高速射影法 これに対して近年提案された高速射影法では、演算の工
夫により積和演算量を2L+20pにまで低減すること
が可能である。以下、高速射影法を簡単に説明する。
【0011】まず、高速射影法では、式(8)の逆行列
演算を直接行わずに、線形予測法を利用してプレフィル
タ係数G(k)を計算することで、O(p3 )の積和演
算量を16pの積和演算量に低減している。この演算量
の低減に関しては本発明とは直接関連しないので、詳細
は省略する。次に、高速射影法では、平滑化係数S
i(k)と近似的フィルタ係数Z(k)を利用すること
で、図6B中のフィルタ係数更新部22およびフィルタ
リング部11で行っている(p+1)Lの積和演算を2
L+4pの積和演算量にまで低減している。この方法に
ついて図7Aに基づいて説明を行う。
演算を直接行わずに、線形予測法を利用してプレフィル
タ係数G(k)を計算することで、O(p3 )の積和演
算量を16pの積和演算量に低減している。この演算量
の低減に関しては本発明とは直接関連しないので、詳細
は省略する。次に、高速射影法では、平滑化係数S
i(k)と近似的フィルタ係数Z(k)を利用すること
で、図6B中のフィルタ係数更新部22およびフィルタ
リング部11で行っている(p+1)Lの積和演算を2
L+4pの積和演算量にまで低減している。この方法に
ついて図7Aに基づいて説明を行う。
【0012】図7Aにおいて、プレフィルタ係数平滑部
31は、プレフィルタ係数を平滑化して平滑係数ベクト
ルを得、近似的フィルタ係数更新部32では、平滑係数
ベクトルの要素を荷重係数として入力信号ベクトルの方
向に近似的フィルタ係数ベクトルを修正し、その近似的
フィルタ係数ベクトルと入力信号ベクトルとの内積をフ
ィルタリング部36で行い、相関計算部33で入力信号
ベクトルの自己相関ベクトルを求め、これと平滑係数ベ
クトルとの内積演算が内積計算部35で行われ、この出
力とフィルタリング部11の出力とが出力信号合成部3
4で加算されて出力信号とされる。この各構成要素で
は、以下の演算が行われる。即ち、プレフィルタ係数平
滑部31では式(11)の演算、近似的フィルタ係数更
新部32では式(12)の演算、相関計算部33では式
(13)の演算、内積演算計算部35では、式(14)
の右辺第2項の内積演算、出力合成部34では式(1
4)の右辺の第1項と第2項との加算が行われる。な
お、図7Aに対応する高速射影法の計算手順を図8に示
す。
31は、プレフィルタ係数を平滑化して平滑係数ベクト
ルを得、近似的フィルタ係数更新部32では、平滑係数
ベクトルの要素を荷重係数として入力信号ベクトルの方
向に近似的フィルタ係数ベクトルを修正し、その近似的
フィルタ係数ベクトルと入力信号ベクトルとの内積をフ
ィルタリング部36で行い、相関計算部33で入力信号
ベクトルの自己相関ベクトルを求め、これと平滑係数ベ
クトルとの内積演算が内積計算部35で行われ、この出
力とフィルタリング部11の出力とが出力信号合成部3
4で加算されて出力信号とされる。この各構成要素で
は、以下の演算が行われる。即ち、プレフィルタ係数平
滑部31では式(11)の演算、近似的フィルタ係数更
新部32では式(12)の演算、相関計算部33では式
(13)の演算、内積演算計算部35では、式(14)
の右辺第2項の内積演算、出力合成部34では式(1
4)の右辺の第1項と第2項との加算が行われる。な
お、図7Aに対応する高速射影法の計算手順を図8に示
す。
【0013】
Z(k+1)=Z(k)+XL (k−p+1)sp (k) (12)
Rp-1 (k+1)=Rp-1 (k)+x(k+1)Xp-1 (k)
−x(k−L+1)Xp-1 (k−L) (13)
y(k+1)=XL (k+1)T Z(k+1)
+Rp-1 (k+1)T Sp-1 (k) (14)
ただし、SP (k)およびSp-1 (k−1)は、それぞ
れp次およびp−1次の平滑化係数ベクトルであって、 Sp (k)=[s1 (k),s2 (k),…,sp (k)]T (15) Sp-1 (k)=[s1 (k),s2 (k),…,sp-1 (k)]T (16) と表される。式(11)からわかるように、Sp (k)
はプレフィルタ係数ベクトルG(k)を用いて再帰的に
計算される。また、Z(k)は、近似的フィルタ係数ベ
クトルであり、XL (k−p+1)とsp (k)を用い
て再帰的に計算される。また、Rp-1 (k+1)は、次
式で表される相関ベクトルである。
れp次およびp−1次の平滑化係数ベクトルであって、 Sp (k)=[s1 (k),s2 (k),…,sp (k)]T (15) Sp-1 (k)=[s1 (k),s2 (k),…,sp-1 (k)]T (16) と表される。式(11)からわかるように、Sp (k)
はプレフィルタ係数ベクトルG(k)を用いて再帰的に
計算される。また、Z(k)は、近似的フィルタ係数ベ
クトルであり、XL (k−p+1)とsp (k)を用い
て再帰的に計算される。また、Rp-1 (k+1)は、次
式で表される相関ベクトルである。
【0014】
Rp-1(k+1)=[XL (k+1) T XL (k),XL (k+1) T XL (k-1),…,
XL (k+1) T XL (k-p-2) ]T (17)
このとき、数学的な証明は省略するが、図7Aの構成、
または、式(11)−(14)の演算により、入力信号
x(k)および誤差信号e(k)を用いて出力信号y
(k+1)を合成することができる。なお、式(13)
(14)の時間パラメータがk+1となっているのに図
7Aではkとなっている。これは時刻kにおいて式(1
3),(14),(11),(12)の順に計算するこ
とを表したものである。
または、式(11)−(14)の演算により、入力信号
x(k)および誤差信号e(k)を用いて出力信号y
(k+1)を合成することができる。なお、式(13)
(14)の時間パラメータがk+1となっているのに図
7Aではkとなっている。これは時刻kにおいて式(1
3),(14),(11),(12)の順に計算するこ
とを表したものである。
【0015】さてここで、以下の2点が指摘できる。ま
ず第1は、以上の方法の演算量に関して、式(11)で
はp、式(12)ではL、式(13)では2p、式(1
4)ではL+p、合計2L+4pの積和演算量で実行で
きる。このことは、図6Bに示した従来の射影法のフィ
ルタ係数更新部22およびフィルタリング部11で行う
式(9)(10)の演算が(p+1)Lの積和演算を必
要とするのに比べて大幅な低減となっている。例えばp
=20,L=1000の場合、射影法が21000であ
るのに対して、高速射影法では2080と、約1/10
の演算量となっている。
ず第1は、以上の方法の演算量に関して、式(11)で
はp、式(12)ではL、式(13)では2p、式(1
4)ではL+p、合計2L+4pの積和演算量で実行で
きる。このことは、図6Bに示した従来の射影法のフィ
ルタ係数更新部22およびフィルタリング部11で行う
式(9)(10)の演算が(p+1)Lの積和演算を必
要とするのに比べて大幅な低減となっている。例えばp
=20,L=1000の場合、射影法が21000であ
るのに対して、高速射影法では2080と、約1/10
の演算量となっている。
【0016】第2の点は、図7Aに示した高速射影法で
は、フィルタ係数H(k+1)を陽に求めることなく出
力信号y(k)を合成している点が、図6Bに示した従
来の射影法との大きな相違点である。より具体的には、
フィルタ係数H(k+1)の代わりに、近似的フィルタ
係数Z(k+1)を利用し、また、H(k+1)とZ
(k+1)との差を補正するために、相関ベクトルR
p-1 (k+1)を利用している。そして、このことが演
算量削減に効果をなしている。 3.2.4 従来のブロック射影法 これまで説明してきた高速射影法では、フィルタ係数H
(k)の代わりに近似的フィルタ係数Z(k)を毎時刻
更新することで、演算量の低減が可能になった。これと
は別に、フィルタ係数H(k)をN(>1)時刻毎に更
新する(ブロック処理)ことで、演算量を削減すること
が可能である。
は、フィルタ係数H(k+1)を陽に求めることなく出
力信号y(k)を合成している点が、図6Bに示した従
来の射影法との大きな相違点である。より具体的には、
フィルタ係数H(k+1)の代わりに、近似的フィルタ
係数Z(k+1)を利用し、また、H(k+1)とZ
(k+1)との差を補正するために、相関ベクトルR
p-1 (k+1)を利用している。そして、このことが演
算量削減に効果をなしている。 3.2.4 従来のブロック射影法 これまで説明してきた高速射影法では、フィルタ係数H
(k)の代わりに近似的フィルタ係数Z(k)を毎時刻
更新することで、演算量の低減が可能になった。これと
は別に、フィルタ係数H(k)をN(>1)時刻毎に更
新する(ブロック処理)ことで、演算量を削減すること
が可能である。
【0017】ブロック処理の手法は、出力y(k)が毎
時刻フィルタ係数を修正する射影法と等しいか否かとい
う観点で2つに分類される。出力y(k)が毎時刻フィ
ルタ係数を修正する射影法と異なる場合は、演算量は削
減されるが、収束速度は遅くなる。一方、毎時刻フィル
タ係数を修正する射影法と同じ出力y(k)をもつブロ
ック処理の場合は、演算量の削減の程度は、毎時刻フィ
ルタ係数を修正する射影法と異なる出力y(k)を持つ
ブロック処理法に比べ、小さくなるが、収束速度の劣化
はみられない利点がある。
時刻フィルタ係数を修正する射影法と等しいか否かとい
う観点で2つに分類される。出力y(k)が毎時刻フィ
ルタ係数を修正する射影法と異なる場合は、演算量は削
減されるが、収束速度は遅くなる。一方、毎時刻フィル
タ係数を修正する射影法と同じ出力y(k)をもつブロ
ック処理の場合は、演算量の削減の程度は、毎時刻フィ
ルタ係数を修正する射影法と異なる出力y(k)を持つ
ブロック処理法に比べ、小さくなるが、収束速度の劣化
はみられない利点がある。
【0018】<A> 出力yが毎時刻フィルタ係数を修
正する射影法と異なるブロック射影法 ブロック処理を射影法に適用した従来のブロック射影法
では、Ep(k),G(k)の計算、そして、H(k+
N)の更新は各々式(7.1),(8),(9)と同様
に行われる。
正する射影法と異なるブロック射影法 ブロック処理を射影法に適用した従来のブロック射影法
では、Ep(k),G(k)の計算、そして、H(k+
N)の更新は各々式(7.1),(8),(9)と同様
に行われる。
【0019】
Ep (k+N−1)=
[d(k+N−1),d(k+N−2),…,d(k+N−p)]T
−XL,p (k+N−1)T H(k) (18)
G(k+N−1)=
[XL,p (k+N−1)T XL,p (k+N−1)]-1Ep (k+N−1)
(19)
H(k+N)=H(k)+μXL,p (k+N−1)G(k+N−1) (20)
出力yの計算は式(10)による1時刻の出力y(k+
1)計算の代わりに、次式でN時刻分y(k+N),
…,y(k+2N−1)を一度に計算する。
1)計算の代わりに、次式でN時刻分y(k+N),
…,y(k+2N−1)を一度に計算する。
【0020】
Y N (k+2N−1)=XL,N (k+2N−1)T H(k+N) (21)
式(18)の右辺第2項の一部または全要素は式(2
1)でk+2Nをk+Nとした式で計算される。式(1
8),(19),(20),(21)が従来の出力y
(k)が毎時刻フィルタ係数を修正する射影法と異なる
ブロック射影法の計算手順であり、図7Bにその機能ブ
ロック図を表わす。図7Bにおいて、プレフィルタ係数
計算部21、フィルタ係数更新部22、ブロックフィル
タリング部41、ブロック誤差計算部42は、それぞれ
式(19),(20),(21),(18)の演算を実
行する。なお、図7Bに対応する従来のブロック射影法
の計算手順を図9に示す。
1)でk+2Nをk+Nとした式で計算される。式(1
8),(19),(20),(21)が従来の出力y
(k)が毎時刻フィルタ係数を修正する射影法と異なる
ブロック射影法の計算手順であり、図7Bにその機能ブ
ロック図を表わす。図7Bにおいて、プレフィルタ係数
計算部21、フィルタ係数更新部22、ブロックフィル
タリング部41、ブロック誤差計算部42は、それぞれ
式(19),(20),(21),(18)の演算を実
行する。なお、図7Bに対応する従来のブロック射影法
の計算手順を図9に示す。
【0021】ブロック処理による演算量が低減するには
2つの側面がある。1つは、フィルタ係数H(k)の修
正が射影法では毎時刻行われるのに対して、ブロック射
影法ではN(>1)時刻に1回に減るというものであ
る。もう1つは、式(20)の右辺のXL,P (k+N−
1)G(k+N−1)と式(21)右辺のXL,N (k+
2N−1)T H(k+N)が高速に計算可能なことによ
る演算量の低減である。後者の内のXL,N (k+2N−
1)T H(k+N)の計算について説明する。
2つの側面がある。1つは、フィルタ係数H(k)の修
正が射影法では毎時刻行われるのに対して、ブロック射
影法ではN(>1)時刻に1回に減るというものであ
る。もう1つは、式(20)の右辺のXL,P (k+N−
1)G(k+N−1)と式(21)右辺のXL,N (k+
2N−1)T H(k+N)が高速に計算可能なことによ
る演算量の低減である。後者の内のXL,N (k+2N−
1)T H(k+N)の計算について説明する。
【0022】XL,N (k+2N−1)T H(k+N)を
効率的に計算するには、次式のように、XL,N (k+2
N−1)をN行N列の部分行列、H(k+N)をN要素
の部分ベクトルに分解する。ただし、L=NM(Mは整
数)となるように、Hは値が0の要素をつけ加える。 XL,N (k+2N−1)T =[XN,N (k+2N−1)T ,XN,N (k+N−1),…, XN,N (k+2N−1−(M−1)N)T ]T (22) H(k+N)T =[H(k+N)0 T ,H(k+N)1 T ,…, H(k+N)M-1 T ]T このように分解すると、XL,N (k+2N−1)T H
(k+N)は、Xの部分行列、Hの部分ベクトル間の積
の和として計算できる。
効率的に計算するには、次式のように、XL,N (k+2
N−1)をN行N列の部分行列、H(k+N)をN要素
の部分ベクトルに分解する。ただし、L=NM(Mは整
数)となるように、Hは値が0の要素をつけ加える。 XL,N (k+2N−1)T =[XN,N (k+2N−1)T ,XN,N (k+N−1),…, XN,N (k+2N−1−(M−1)N)T ]T (22) H(k+N)T =[H(k+N)0 T ,H(k+N)1 T ,…, H(k+N)M-1 T ]T このように分解すると、XL,N (k+2N−1)T H
(k+N)は、Xの部分行列、Hの部分ベクトル間の積
の和として計算できる。
【0023】
XL,N (k+2N−1)T H(k+N)=
Σi=0 M-1 XN,N (k+2N−1−iN)T H(k+N)i (23)
Σの中の各々の正方行列XN,N (k+2N−1−iN)
とベクトルH(k+N)i の積は、入力信号のベクトル
X2N(k+2N−1−iN)とフィルタ係数H(k+
N)i の畳み込みとみなすことができる。畳み込みを高
速に行う方法として、FFTを用いる方法、Winog
radの高速畳み込み法に基づく方法(詳しくは文献B
lahut,FAST ALGORITHMS FOR
DIGITAL SIGNAL PROCESSIN
G,ADDISON−WESTLEY,1987,第
3,9章)がある。これらの、高速畳み込み法を式(2
3)に適用することで、YN (k+N)の演算量が低減
可能である。
とベクトルH(k+N)i の積は、入力信号のベクトル
X2N(k+2N−1−iN)とフィルタ係数H(k+
N)i の畳み込みとみなすことができる。畳み込みを高
速に行う方法として、FFTを用いる方法、Winog
radの高速畳み込み法に基づく方法(詳しくは文献B
lahut,FAST ALGORITHMS FOR
DIGITAL SIGNAL PROCESSIN
G,ADDISON−WESTLEY,1987,第
3,9章)がある。これらの、高速畳み込み法を式(2
3)に適用することで、YN (k+N)の演算量が低減
可能である。
【0024】FFTを用いる方法の高速の畳み込みの原
理を簡単に説明する。よく知られているように、時間領
域での畳み込みは周波数領域での積になる。 F[x*h]=F[x]★F[h] (24) ここで、F[ ]はフーリエ変換を、*は畳み込みを、
★は要素ごとの積を表す。この関係を利用すると、xと
hの畳み込みは次式で計算することができる。
理を簡単に説明する。よく知られているように、時間領
域での畳み込みは周波数領域での積になる。 F[x*h]=F[x]★F[h] (24) ここで、F[ ]はフーリエ変換を、*は畳み込みを、
★は要素ごとの積を表す。この関係を利用すると、xと
hの畳み込みは次式で計算することができる。
【0025】
x*h=F-1[F[x*h]]
=F-1[F[x]★F[h]] (25)
ここで、F-1[ ]はフーリエ逆変換を表す。フーリエ
変換に高速フーリエ変換(FFT)を用いることで、高
速に畳み込みを行うことが可能である。式(20)、
(21)に適用すると次式になる。
変換に高速フーリエ変換(FFT)を用いることで、高
速に畳み込みを行うことが可能である。式(20)、
(21)に適用すると次式になる。
【0026】
H(k+N)m
=H(k)m +μF-1[F[X2p(k+N−1−mN)]
★F[G(k+N−1)]]
m=0,1,…,L/p−1 (26)
YN (k+2N−1)
=Σm=0 M-1 F-1[F[X2N(k+2N−1−mN)]
★F[H(k+N)m ]]
=F-1[Σm=0 M-1 F[X2N(k+2N−1−mN)]
★F[H(k+N)m ]] (27)
式(21)に対応する式(27)を例にとると、データ
数2NのFFTが、X 2N(k+2N−1)に1回(X2N
(k+2N−1−mN),m>1のFFTは時刻k+2
N−1−mNで行われるので回数に入れない)、H(k
+N)m にM回、データ数2Nの逆FFTが1回、そし
て、複素数の積和演算が2NM回必要である。この演算
量を実数の加算と積和演算の合計数に置き換えると、デ
ータ長が2NのFFT、逆FFTは約8Nlog2 (2
N)、1回の複素数の積和演算は4となる。これらよ
り、式(27)に必要な演算量は約(M+2)8Nlo
g2N+8NMとなり、ほぼ、MNlog2 Nに比例す
る。一方、式(27)を式(21)のように直接計算し
た場合の演算量は、NL(=N2 M)となるので、FF
Tを用いることで演算量が削減されることがわかる。
数2NのFFTが、X 2N(k+2N−1)に1回(X2N
(k+2N−1−mN),m>1のFFTは時刻k+2
N−1−mNで行われるので回数に入れない)、H(k
+N)m にM回、データ数2Nの逆FFTが1回、そし
て、複素数の積和演算が2NM回必要である。この演算
量を実数の加算と積和演算の合計数に置き換えると、デ
ータ長が2NのFFT、逆FFTは約8Nlog2 (2
N)、1回の複素数の積和演算は4となる。これらよ
り、式(27)に必要な演算量は約(M+2)8Nlo
g2N+8NMとなり、ほぼ、MNlog2 Nに比例す
る。一方、式(27)を式(21)のように直接計算し
た場合の演算量は、NL(=N2 M)となるので、FF
Tを用いることで演算量が削減されることがわかる。
【0027】<B> 毎時刻フィルタ係数を修正する射
影法と同じ出力を持つブロック処理射影法 先に3.2.4<A>において説明した、ブロック射影
法では、ブロック長Nに一度しか、適応フィルタを修正
しないので、ブロック長Nを長くするにつれて、適応フ
ィルタの収束が遅くなるという欠点がある。
影法と同じ出力を持つブロック処理射影法 先に3.2.4<A>において説明した、ブロック射影
法では、ブロック長Nに一度しか、適応フィルタを修正
しないので、ブロック長Nを長くするにつれて、適応フ
ィルタの収束が遅くなるという欠点がある。
【0028】この欠点に対して、高速射影法における、
相関ベクトルRp-1 (k)と平滑化係数ベクトルS
p (k)の計算範囲を拡張することで解決を図り、適応
フィルタの収束速度の低下を伴うことなく、射影法の演
算量のブロック処理による低減を可能にするのが、ここ
で説明する毎時刻フィルタ係数を修正する射影法と同じ
出力を持つブロック処理である。
相関ベクトルRp-1 (k)と平滑化係数ベクトルS
p (k)の計算範囲を拡張することで解決を図り、適応
フィルタの収束速度の低下を伴うことなく、射影法の演
算量のブロック処理による低減を可能にするのが、ここ
で説明する毎時刻フィルタ係数を修正する射影法と同じ
出力を持つブロック処理である。
【0029】以下この方法の説明を行う。先に述べた毎
時刻適応フィルタを修正する高速射影法において、式
(12)からZ(k+i+1),(i>=0)はZ
(k)を用いて表すことができる。 Z(k+1+i)=Z(k+i)+sp (k+i)XL (k+1+i−p) =Z(k+i−1)+sp (k+i−1)XL (k+i−p) +sp (k+i)XL (k+1+i−p) … =Z(k)+Σj=0 i sp (k+i−j)XL (k+1+i−p−j) =Z(k)+XL,i+1 (k+1+i−p)Sp+i (k+i)p 〜p+i (28) であり、特に、i=N−1のときに近似的フィルタ係数
Z(k)からZ(k+N)への更新式になる。
時刻適応フィルタを修正する高速射影法において、式
(12)からZ(k+i+1),(i>=0)はZ
(k)を用いて表すことができる。 Z(k+1+i)=Z(k+i)+sp (k+i)XL (k+1+i−p) =Z(k+i−1)+sp (k+i−1)XL (k+i−p) +sp (k+i)XL (k+1+i−p) … =Z(k)+Σj=0 i sp (k+i−j)XL (k+1+i−p−j) =Z(k)+XL,i+1 (k+1+i−p)Sp+i (k+i)p 〜p+i (28) であり、特に、i=N−1のときに近似的フィルタ係数
Z(k)からZ(k+N)への更新式になる。
【0030】
Z(k+N)=
Z(k)+XL,N (k+N−p)Sp+N-1 (k+N−1)p 〜p+N-1
(29)
式(28)における、Sp+i (k+i)は式(15)の
次数pのSp (k)に対し次数p+iとされたものであ
り、これを拡張平滑化係数ベクトルと呼び、S p+i (k
+i)p 〜p+i は拡張平滑化ベクトルSp+i (k+i)
のp番目の要素からp+i番目の要素からなるベクトル
を意味する。Sp+i (k+i)p 〜p+iは次式で再帰的
に計算される。
次数pのSp (k)に対し次数p+iとされたものであ
り、これを拡張平滑化係数ベクトルと呼び、S p+i (k
+i)p 〜p+i は拡張平滑化ベクトルSp+i (k+i)
のp番目の要素からp+i番目の要素からなるベクトル
を意味する。Sp+i (k+i)p 〜p+iは次式で再帰的
に計算される。
【0031】
また、拡張平滑化係数ベクトルSp+i (k+i)は、式
(16)のSp-1 (k+i)とSp+i (k+i)p 〜
p+i を結合することで次式で計算される。
(16)のSp-1 (k+i)とSp+i (k+i)p 〜
p+i を結合することで次式で計算される。
【0032】
式(30)と平滑化係数ベクトルの回帰式(11)を順
に式(31)に代入することで、拡張平滑ベクトルS
p+i (k+i)の回帰式が得られる。
に式(31)に代入することで、拡張平滑ベクトルS
p+i (k+i)の回帰式が得られる。
【0033】
ところで、式(14)において、kをk+i−1とする
ことで、次式が得られる。
ことで、次式が得られる。
【0034】
y(k+i)=XL (k+i)T Z(k+i)
+Rp-1 (k+i)T Sp-1 (k+i−1) (33)
また、式(28)の右辺から、Z(k+i)をZ(k)
で表わす式が得られる。 Z(k+i)= Z(k)+XL,i (k+i−p)Sp+i-1 (k+i−1)p 〜p+i-1 (34) 式(34)を式(33)に代入すると、次式となる。 y(k+i)= XL (k+i)T [Z(k) +XL,i (k+i−p)Sp+i-1 (k+i−1)p 〜p+i-1 ] +Rp-1 (k+i)T Sp-1 (k+i−1) =XL (k+i)T Z(k) +XL (k+i)T XL,i (k+i−p) Sp+i-1 (k+i−1)p 〜p+i-1 +Rp-1 (k+i)T Sp-1 (k+i−1) =XL (k+i)T Z(k) +Rp+i-1 (k+i)T p 〜p+i-1 Sp+i-1 (k+i−1)p 〜p+i-1 +Rp-1 (k+i)T Sp-1 (k+i−1) (35) ここで、Rp+i-1 (k+i)を拡張相関ベクトルと呼
び、Rp+i-1 (k+i) p 〜p+i-1 は拡張相関ベクトル
のp番目の要素からp+i−1番目の要素からなるベク
トルを意味する。Rp+i-1 (k+i)p 〜p+i-1 はR
N+p-2 (k+i)p〜N+p-2 の一部またはすべての要素
であり、RN+p-2 (k+i)p 〜N+p-2 は次式で回帰的
に計算される。
で表わす式が得られる。 Z(k+i)= Z(k)+XL,i (k+i−p)Sp+i-1 (k+i−1)p 〜p+i-1 (34) 式(34)を式(33)に代入すると、次式となる。 y(k+i)= XL (k+i)T [Z(k) +XL,i (k+i−p)Sp+i-1 (k+i−1)p 〜p+i-1 ] +Rp-1 (k+i)T Sp-1 (k+i−1) =XL (k+i)T Z(k) +XL (k+i)T XL,i (k+i−p) Sp+i-1 (k+i−1)p 〜p+i-1 +Rp-1 (k+i)T Sp-1 (k+i−1) =XL (k+i)T Z(k) +Rp+i-1 (k+i)T p 〜p+i-1 Sp+i-1 (k+i−1)p 〜p+i-1 +Rp-1 (k+i)T Sp-1 (k+i−1) (35) ここで、Rp+i-1 (k+i)を拡張相関ベクトルと呼
び、Rp+i-1 (k+i) p 〜p+i-1 は拡張相関ベクトル
のp番目の要素からp+i−1番目の要素からなるベク
トルを意味する。Rp+i-1 (k+i)p 〜p+i-1 はR
N+p-2 (k+i)p〜N+p-2 の一部またはすべての要素
であり、RN+p-2 (k+i)p 〜N+p-2 は次式で回帰的
に計算される。
【0035】
RN+p-2 (k+i)T p 〜N+p-2 =XL (k+i)T XL,N-1 (k+i−p)
=RN+p-2 (k+i−1)T p 〜N+p-2
+x(k+i)XN-1 (k+i−p)
−x(k+i−L)XN-1 (k+i−p−L) (36)
また、拡張相関ベクトルは式(17)のRp-1 (k+
i)とRp+i-1 (k+i)p 〜p+i-1 を結合することで
計算される。
i)とRp+i-1 (k+i)p 〜p+i-1 を結合することで
計算される。
【0036】
拡張相関ベクトルの時間更新は、式(13)と式(3
6)を式(37)に代入して次のようにかける。
6)を式(37)に代入して次のようにかける。
【0037】
RN+p-2 (k+i)=
RN+p-2 (k+i−1)+x(k+i)XN+p-2 (k+i−1)
−x(k+i−L)XN+p-2 (k+i−L−1) (38)
拡張平滑化係数ベクトル、拡張相関ベクトルを用いる
と、式(35)は次のようにかける。 y(k+i)=y′(k+i)+Ri+p-1 (k+i)T Si+p-1 (k+i−1) (39) y′(k+i)=XL (k+i)T Z(k) (40) y(k),Ep (k),G(k),Sp (k)は高速射
影法の手順で計算できるので、i=0として式(39)
を用いてy(k)が計算できる。式(7.2)を用いて
y(k)からe(k),Ep (k)が計算でき、E
p (k)からG(k)が式(8)により計算できる。G
(k)から、式(11),(30),(31)を用い
て、Sp+1 (k)を計算でき、式(39)により、y
(k+1)が計算できる。このように、順にy(k+
i),(i=0,1,…,N−1)は計算することがで
き、時刻k+N−1において式(29)からZ(k+
N)が計算できる。
と、式(35)は次のようにかける。 y(k+i)=y′(k+i)+Ri+p-1 (k+i)T Si+p-1 (k+i−1) (39) y′(k+i)=XL (k+i)T Z(k) (40) y(k),Ep (k),G(k),Sp (k)は高速射
影法の手順で計算できるので、i=0として式(39)
を用いてy(k)が計算できる。式(7.2)を用いて
y(k)からe(k),Ep (k)が計算でき、E
p (k)からG(k)が式(8)により計算できる。G
(k)から、式(11),(30),(31)を用い
て、Sp+1 (k)を計算でき、式(39)により、y
(k+1)が計算できる。このように、順にy(k+
i),(i=0,1,…,N−1)は計算することがで
き、時刻k+N−1において式(29)からZ(k+
N)が計算できる。
【0038】まとめると、図11に示すように、まず初
期化として、Rp+N-1 (0)=XL (0)T [XL (−
1),XL (−2),…,XL (−p−N+1)]T と
し、Ep-1 (0)=Sp-1 (0)=XL (0)=[0,
0,…,0]T とし(S100 )、次に、k=1とし(S
101 )、次に入力ベクトルと近似的ベクトルとの内積演
算 y′(k+i)=XL (k+i)T Z(k),i=0,
1,…,N−1 を行う(S102 )。その後、i=0,1,…,N−1に
対して、式(38),(39),(7.2),(8),
(32)の演算を各時刻ごとに行い、i=N−1の演算
が終了すると(S103 )、式(29)を演算して、Z
(k)からZ(k+N)に更新する(S104 )。次に、
kをk+NとしてステップS102 に戻る(S105 )。
期化として、Rp+N-1 (0)=XL (0)T [XL (−
1),XL (−2),…,XL (−p−N+1)]T と
し、Ep-1 (0)=Sp-1 (0)=XL (0)=[0,
0,…,0]T とし(S100 )、次に、k=1とし(S
101 )、次に入力ベクトルと近似的ベクトルとの内積演
算 y′(k+i)=XL (k+i)T Z(k),i=0,
1,…,N−1 を行う(S102 )。その後、i=0,1,…,N−1に
対して、式(38),(39),(7.2),(8),
(32)の演算を各時刻ごとに行い、i=N−1の演算
が終了すると(S103 )、式(29)を演算して、Z
(k)からZ(k+N)に更新する(S104 )。次に、
kをk+NとしてステップS102 に戻る(S105 )。
【0039】この手順に対応する機能ブロック図は図1
0に示し、図7A、図7Bと対応する部分と同一符号を
付けて示す。入力信号ベクトルから相関計算部33で式
(13)が演算され、また入力信号ベクトルから相関計
算拡張部52で式(36)が演算され、これら相関計算
部33からの相関ベクトルと、相関計算拡張部52から
の相関ベクトルの拡張部のベクトルとが相関ベクトル連
結部39で式(37)により連結されて拡張相関ベクト
ルが得られる。要は、これら相関計算部33、相関計算
拡張部52、相関ベクトル連結部53は式(38)を演
算するものであり、拡張相関ベクトル演算部56を構成
する。
0に示し、図7A、図7Bと対応する部分と同一符号を
付けて示す。入力信号ベクトルから相関計算部33で式
(13)が演算され、また入力信号ベクトルから相関計
算拡張部52で式(36)が演算され、これら相関計算
部33からの相関ベクトルと、相関計算拡張部52から
の相関ベクトルの拡張部のベクトルとが相関ベクトル連
結部39で式(37)により連結されて拡張相関ベクト
ルが得られる。要は、これら相関計算部33、相関計算
拡張部52、相関ベクトル連結部53は式(38)を演
算するものであり、拡張相関ベクトル演算部56を構成
する。
【0040】プレフィルタ係数計算部21で入力信号ベ
クトルと、誤差信号ベクトルとにより式(8)を演算し
てプレフィルタ係数を演算し、このプレフィルタ係数は
プレフィルタ係数平滑部31で式(30)の演算により
求められ、この拡張部分のベクトルと平滑部31よりの
平滑係数ベクトルとが平滑係数ベクトル連結部55で式
(31)の演算により連結されて拡張平滑係数ベクトル
が得られる。要は、プレフィルタ係数平滑部31、平滑
係数拡張部54、平滑係数ベクトル連結部55は式(3
2)を演算する拡張平滑係数ベクトル計算部57を構成
している。
クトルと、誤差信号ベクトルとにより式(8)を演算し
てプレフィルタ係数を演算し、このプレフィルタ係数は
プレフィルタ係数平滑部31で式(30)の演算により
求められ、この拡張部分のベクトルと平滑部31よりの
平滑係数ベクトルとが平滑係数ベクトル連結部55で式
(31)の演算により連結されて拡張平滑係数ベクトル
が得られる。要は、プレフィルタ係数平滑部31、平滑
係数拡張部54、平滑係数ベクトル連結部55は式(3
2)を演算する拡張平滑係数ベクトル計算部57を構成
している。
【0041】拡張相関ベクトル演算部56よりの拡張相
関ベクトルと拡張平滑係数ベクトル計算部57よりの拡
張平滑係数ベクトルとの内積演算、つまり、式(39)
の右辺第2項の演算がなされる。平滑係数拡張部54よ
りの拡張部分ベクトルが近似的フィルタ係数ブロック更
新部51へ入力されて式(34)が演算される。この演
算結果の近似的フィルタ係数と入力ベクトルとの内積演
算がブロックフィルタリング部41で行われ、その演算
結果と内積演算計算部35からの演算結果とが出力合成
部34で加算されて出力信号とされる。ブロックフィル
タリング部41、近似的フィルタ係数ブロック更新部5
1はN時刻ごとに演算され、その他の各部は各時刻ごと
に行われ、拡張平滑係数ベクトル部57では、ブロック
フィルタリング部41で演算され、次の演算まで各時刻
ごとに字数がpから1次づつ拡張され、このことが繰り
返され、拡張相関ベクトル計算部56では、常時、拡張
平滑係数ベクトルの最大次数p+N−1と同一次数の相
関ベクトルが演算され内積演算計算部35に対しては、
拡張平滑係数ベクトル計算部57よりのその時刻での拡
張平滑係数ベクトルの字数と同一次数分が入力される。
関ベクトルと拡張平滑係数ベクトル計算部57よりの拡
張平滑係数ベクトルとの内積演算、つまり、式(39)
の右辺第2項の演算がなされる。平滑係数拡張部54よ
りの拡張部分ベクトルが近似的フィルタ係数ブロック更
新部51へ入力されて式(34)が演算される。この演
算結果の近似的フィルタ係数と入力ベクトルとの内積演
算がブロックフィルタリング部41で行われ、その演算
結果と内積演算計算部35からの演算結果とが出力合成
部34で加算されて出力信号とされる。ブロックフィル
タリング部41、近似的フィルタ係数ブロック更新部5
1はN時刻ごとに演算され、その他の各部は各時刻ごと
に行われ、拡張平滑係数ベクトル部57では、ブロック
フィルタリング部41で演算され、次の演算まで各時刻
ごとに字数がpから1次づつ拡張され、このことが繰り
返され、拡張相関ベクトル計算部56では、常時、拡張
平滑係数ベクトルの最大次数p+N−1と同一次数の相
関ベクトルが演算され内積演算計算部35に対しては、
拡張平滑係数ベクトル計算部57よりのその時刻での拡
張平滑係数ベクトルの字数と同一次数分が入力される。
【0042】以上説明した、毎時刻フィルタ係数を修正
する射影法と同じ出力を持つブロック射影法の演算量を
求める。まず、同方法の手順を示した図11のS
102 は、3.2.4<A>で説明したように、FFTを
用いた高速畳み込みが利用できその演算量は、約8(M
+2)Nlog2 2N+8NMであり、1時刻当たりで
は、8(M+2)log2 2N+8Mである。S103 で
は、1時刻当たり、式(38)に2(N+p−2)回、
式(39)にp+N/2回、式(7.2)にp回、式
(8)には15p回、式(32)には、p回の演算が必
要であり、合計すると2.5N+20pとなる。S104
は3.2.4<A>で説明した要領でFFTを用いて高
速に計算できる。つまり、式(29)は式(20)と同
様に次のようにかけるので、 Z(k+N)m = Z(k)m +μF-1[F[X2N(k+N−p−mN)] ★F[S(k+N−1)p,N+p-1 ]] m=0,1,…,L/N−1 (41) データ長が2NのX2N(k+N−p)をFFTするのに
1回、S(k+N−1)にN個0を付加したベクトルを
長さ2NのデータのFFTが1回、そして、逆FFTが
L/N(=M)回必要である。また、複素数の積和演算
が2MN回必要である。これら、FFT、逆FFT、複
素数の積和演算に要する演算量は約8(M+2)Nlo
g2 2N+8NMであり、1時刻当たりでは、8(M+
2)log 2 2N+8Mである。以上より、S102 、S
103 そしてS104 に要する1時刻当たりの演算量をTと
かくと、 T=16(M+2)log2 (2N)+2.5N+16M+20p (42) と、見積られる。
する射影法と同じ出力を持つブロック射影法の演算量を
求める。まず、同方法の手順を示した図11のS
102 は、3.2.4<A>で説明したように、FFTを
用いた高速畳み込みが利用できその演算量は、約8(M
+2)Nlog2 2N+8NMであり、1時刻当たりで
は、8(M+2)log2 2N+8Mである。S103 で
は、1時刻当たり、式(38)に2(N+p−2)回、
式(39)にp+N/2回、式(7.2)にp回、式
(8)には15p回、式(32)には、p回の演算が必
要であり、合計すると2.5N+20pとなる。S104
は3.2.4<A>で説明した要領でFFTを用いて高
速に計算できる。つまり、式(29)は式(20)と同
様に次のようにかけるので、 Z(k+N)m = Z(k)m +μF-1[F[X2N(k+N−p−mN)] ★F[S(k+N−1)p,N+p-1 ]] m=0,1,…,L/N−1 (41) データ長が2NのX2N(k+N−p)をFFTするのに
1回、S(k+N−1)にN個0を付加したベクトルを
長さ2NのデータのFFTが1回、そして、逆FFTが
L/N(=M)回必要である。また、複素数の積和演算
が2MN回必要である。これら、FFT、逆FFT、複
素数の積和演算に要する演算量は約8(M+2)Nlo
g2 2N+8NMであり、1時刻当たりでは、8(M+
2)log 2 2N+8Mである。以上より、S102 、S
103 そしてS104 に要する1時刻当たりの演算量をTと
かくと、 T=16(M+2)log2 (2N)+2.5N+16M+20p (42) と、見積られる。
【0043】
【発明が解決しようとする課題】従来の毎時刻フィルタ
係数を修正する射影法(以下、従来法)と同じ出力y
(k)を持つブロック射影法では、演算量Tを最小にす
るブロック長が存在し、それ以上に長いブロック長で
は、演算量の低減効果が小さくなるか、または無くなる
という問題点がある。この理由は、式(38)で計算さ
れる拡張相関ベクトルRN+p-2 の時間更新に要する演算
量2Nと式(39)右辺第2項で計算される拡張相関ベ
クトルと拡張平滑係数ベクトルの積に要する演算量0.
5Nがブロック長Nに比例しているので、ブロック長N
を長くするとともに増加することにある。式を用いれ
ば、式(42)のTのブロック長Nに関して微係数が正
となることで分かる。
係数を修正する射影法(以下、従来法)と同じ出力y
(k)を持つブロック射影法では、演算量Tを最小にす
るブロック長が存在し、それ以上に長いブロック長で
は、演算量の低減効果が小さくなるか、または無くなる
という問題点がある。この理由は、式(38)で計算さ
れる拡張相関ベクトルRN+p-2 の時間更新に要する演算
量2Nと式(39)右辺第2項で計算される拡張相関ベ
クトルと拡張平滑係数ベクトルの積に要する演算量0.
5Nがブロック長Nに比例しているので、ブロック長N
を長くするとともに増加することにある。式を用いれ
ば、式(42)のTのブロック長Nに関して微係数が正
となることで分かる。
【0044】
∂T/∂N=16((−L/N2 )log2 (2N)
+((L/N)+2)/(Nln2))
+2.5−16(L/N2 )
=N-2(16(−Llog2 (2N)
+(L+2N)/ln2)+2.5N2 −16L) (43)
N-2は常に正であり、16(−Llog2 (2N)〜+
(L+2N)/ln2)+2.5N2 −16LはL>>
Nでは負の値をとるが、Nが大きくなるとN、N 2 の項
に比べてlog2 (2N)の項が無視できて、下に凸な
Nの2次式とみなせるので、値は正になる。つまり、演
算量はN<<Lでは減少するが、あるN以上では、増加
する。
(L+2N)/ln2)+2.5N2 −16LはL>>
Nでは負の値をとるが、Nが大きくなるとN、N 2 の項
に比べてlog2 (2N)の項が無視できて、下に凸な
Nの2次式とみなせるので、値は正になる。つまり、演
算量はN<<Lでは減少するが、あるN以上では、増加
する。
【0045】フィルタ長Lが与えられたときに、演算量
Tが最小となるNは、式(43)の右辺が0となるブロ
ック長Nである。式(43)の右辺を0とおいた方程式
は、以下に示すようにNに関する再帰計算により求める
ことができる。まず、式(43)の右辺を0とおいた方
程式をNの2次式とみなしてとく。 N=[16/ln2+4√(1/ln2)2+2.5L(log2(2N)-1/ln2+1) ]/2.5 =9+1.6 √2.5Llog2(2N) =9+2.5 √Llog2(2N) (44) 式(44)の左辺のブロック長Nを右辺に代入仕直すこ
とで、再帰計算を行う。i回目の再帰計算によるブロッ
ク長NをN(i)と表わすことにすると式(44)は次
式になる。
Tが最小となるNは、式(43)の右辺が0となるブロ
ック長Nである。式(43)の右辺を0とおいた方程式
は、以下に示すようにNに関する再帰計算により求める
ことができる。まず、式(43)の右辺を0とおいた方
程式をNの2次式とみなしてとく。 N=[16/ln2+4√(1/ln2)2+2.5L(log2(2N)-1/ln2+1) ]/2.5 =9+1.6 √2.5Llog2(2N) =9+2.5 √Llog2(2N) (44) 式(44)の左辺のブロック長Nを右辺に代入仕直すこ
とで、再帰計算を行う。i回目の再帰計算によるブロッ
ク長NをN(i)と表わすことにすると式(44)は次
式になる。
【0046】
N(i+1)=9+2.5√Llog2 (2N(i)) (45)
【0047】
【課題を解決するための手段】この発明では、上記の課
題に対して、従来法における、式(40)で計算される
出力計算、式(29)で計算される近似的フィルタ係数
Zの更新において異なるブロック長を用いることで解決
を図り、従来法では得られない演算量の低減や遅延量の
短縮を可能にするものである。
題に対して、従来法における、式(40)で計算される
出力計算、式(29)で計算される近似的フィルタ係数
Zの更新において異なるブロック長を用いることで解決
を図り、従来法では得られない演算量の低減や遅延量の
短縮を可能にするものである。
【0048】この発明によれば、入力信号を可変フィル
タ処理し、その出力信号と所望信号との誤差信号を求
め、その誤差信号のパワーを最小とするように、上記可
変フィルタの特性を適応的に制御する方法において、上
記入力信号ベクトルの自己相関行列の逆行列と、上記誤
差信号ベクトルとの積に相当する演算を行ってプレフィ
ルタ係数ベクトルを求める第1ステップと、上記プレフ
ィルタ係数ベクトルを平滑化して平滑係数ベクトルを得
る第2ステップと、上記平滑係数ベクトルの要素を加重
係数として、上記入力信号ベクトルの方向に近似的フィ
ルタ係数ベクトルを修正する第3ステップと、上記入力
信号ベクトルと上記近似的フィルタ係数ベクトルとの内
積演算を行う第4ステップと、上記入力信号ベクトルの
自己相関ベクトルを求める第5ステップと、上記平滑係
数ベクトルと上記自己相関ベクトルとの内積演算を行う
第6ステップと、上記第4ステップで得られた内積演算
結果と、上記第6ステップで得られた内積演算結果とを
加算して上記出力信号を得る第7ステップとを有するも
のであって、上記第3ステップの修正及び上記第4ステ
ップの演算はそれぞれN1 、N2 時刻(N1 、N2 は2
以上の整数、N1 ≠N2 )ごとに行い、上記第3、第4
ステップ以外の各ステップは各時刻ごとに演算を行い、
上記第3ステップの演算が行われるごとに、次に上記第
3ステップの演算されるまでの各時刻において、上記第
2ステップで得る上記平滑係数ベクトルの次数を1ずつ
増加し、次の第3ステップの演算時刻での演算でもとの
次式に戻し、上記第5ステップでは上記平滑係数ベクト
ルの最大次数と等しい次数の上記自己相関ベクトルを演
算し、上記第6ステップでは上記平滑係数ベクトルの次
数と等しい次数だけ上記第5ステップで得られた自己相
関ベクトルから用いることを特徴とする。
タ処理し、その出力信号と所望信号との誤差信号を求
め、その誤差信号のパワーを最小とするように、上記可
変フィルタの特性を適応的に制御する方法において、上
記入力信号ベクトルの自己相関行列の逆行列と、上記誤
差信号ベクトルとの積に相当する演算を行ってプレフィ
ルタ係数ベクトルを求める第1ステップと、上記プレフ
ィルタ係数ベクトルを平滑化して平滑係数ベクトルを得
る第2ステップと、上記平滑係数ベクトルの要素を加重
係数として、上記入力信号ベクトルの方向に近似的フィ
ルタ係数ベクトルを修正する第3ステップと、上記入力
信号ベクトルと上記近似的フィルタ係数ベクトルとの内
積演算を行う第4ステップと、上記入力信号ベクトルの
自己相関ベクトルを求める第5ステップと、上記平滑係
数ベクトルと上記自己相関ベクトルとの内積演算を行う
第6ステップと、上記第4ステップで得られた内積演算
結果と、上記第6ステップで得られた内積演算結果とを
加算して上記出力信号を得る第7ステップとを有するも
のであって、上記第3ステップの修正及び上記第4ステ
ップの演算はそれぞれN1 、N2 時刻(N1 、N2 は2
以上の整数、N1 ≠N2 )ごとに行い、上記第3、第4
ステップ以外の各ステップは各時刻ごとに演算を行い、
上記第3ステップの演算が行われるごとに、次に上記第
3ステップの演算されるまでの各時刻において、上記第
2ステップで得る上記平滑係数ベクトルの次数を1ずつ
増加し、次の第3ステップの演算時刻での演算でもとの
次式に戻し、上記第5ステップでは上記平滑係数ベクト
ルの最大次数と等しい次数の上記自己相関ベクトルを演
算し、上記第6ステップでは上記平滑係数ベクトルの次
数と等しい次数だけ上記第5ステップで得られた自己相
関ベクトルから用いることを特徴とする。
【0049】
【発明の実施の形態】以下、この発明の実施の形態を説
明する。一般にブロック処理を行う適応アルゴリズムに
おいては、ブロック処理により演算量を低減可能な計算
ステップが2箇所ある。先に述べた毎時刻適応フィルタ
を修正する射影法と同じ出力を持つブロック射影法にお
いては、図11中のS102 で式(40)により計算され
る出力計算と図11中S104 で式(29)で計算される
近似的フィルタ係数Zの更新である。ところで、先に述
べた従来法において課題となった、長いブロック長Nに
対して演算低減効果が小さくなる理由は式(38)で計
算される拡張相関ベクトルRN+p-2 の時間更新に要する
演算量2Nと式(39)右辺第2項で計算される拡張相
関ベクトルと拡張平滑係数ベクトルの積に要する演算量
0.5Nがブロック長Nに比例しているので、ブロック
長Nを長くするとともに増加することであった。これら
ブロック長に比例して増加する演算は近似的フィルタの
更新をブロック長の間だけ保留しているために生じてい
る。そこで、上記の2箇所のブロック処理において異な
るブロック長を用いることで上記の課題を解決すること
ができる。
明する。一般にブロック処理を行う適応アルゴリズムに
おいては、ブロック処理により演算量を低減可能な計算
ステップが2箇所ある。先に述べた毎時刻適応フィルタ
を修正する射影法と同じ出力を持つブロック射影法にお
いては、図11中のS102 で式(40)により計算され
る出力計算と図11中S104 で式(29)で計算される
近似的フィルタ係数Zの更新である。ところで、先に述
べた従来法において課題となった、長いブロック長Nに
対して演算低減効果が小さくなる理由は式(38)で計
算される拡張相関ベクトルRN+p-2 の時間更新に要する
演算量2Nと式(39)右辺第2項で計算される拡張相
関ベクトルと拡張平滑係数ベクトルの積に要する演算量
0.5Nがブロック長Nに比例しているので、ブロック
長Nを長くするとともに増加することであった。これら
ブロック長に比例して増加する演算は近似的フィルタの
更新をブロック長の間だけ保留しているために生じてい
る。そこで、上記の2箇所のブロック処理において異な
るブロック長を用いることで上記の課題を解決すること
ができる。
【0050】二つのブロック処理におけるブロック長を
独立したことにより、この発明法の処理は従来法と処理
順序または処理のタイミングにおいて異なったものとな
る。出力の計算でのブロック長をN1 、近似的フィルタ
係数更新でのブロック長をN2 と書くと、式(40),
(38),(39),(7.2),(8),(32)そ
して(29)で与えられる従来法での処理に対応して、
本発明法の処理は以下のようになる。
独立したことにより、この発明法の処理は従来法と処理
順序または処理のタイミングにおいて異なったものとな
る。出力の計算でのブロック長をN1 、近似的フィルタ
係数更新でのブロック長をN2 と書くと、式(40),
(38),(39),(7.2),(8),(32)そ
して(29)で与えられる従来法での処理に対応して、
本発明法の処理は以下のようになる。
【0051】
y′(k1 +i)=XL (k2 +i)T Z(k2 ),
i=0,…,N1 −1 (46)
RN2+p-2(k1 )=RN2+p-2(k−1)
+x(k1 )XN2+p-2(k1 −1)
−x(k1 −L)XN2+p-2(k1 −L−1) (47)
y(k1 )=y′(k1 )+Rn2+p-1(k1 )T Sn2+p-1(k1 −1)
(48)
e(k1 )=d(k1 )−y(k1 ) (49)
Z(k2 +N2)=Z(k2 )
+XL,N2(k2 +N2 −p)SN2+p-1(k+N2 −1)p-N2+p-1 (53)
ここで、k2 は近似的フィルタ係数Zの修正に関する時
刻の指標であり、k1 はその他の変数の時刻の指標であ
る。また、n2 =k1 mod N2 −1なる値をもつ変
数である。図1はこの発明法の計算手順を示している。
まず、初期化として、 Rp+N2-1(0)=XL (0) T [XL (-1), XL (-2), …, XL (-p-N2+1) ]T Sp-1 (0)=[0,0,…,0]T ZL-1 (1)=[0,0,…,0]T Ep-1 (0)=[y(0),y(−1),…,y(−p+2)]T とし(S0 )、次にk1=k2=1,n 2 =0とし(S
1 )、次に、もし、k1mod N1 =1であるかの判
定をし(C1 )、判定C1 が真であれば、入力ベクトル
と近似的フィルタ係数との内積演算 y′(k1 +i)=XL (k2 +i)T Z(k2 ), i=0,…,N1 −1 (S2 )を行う。その後、式(47),(48),(4
9),(50),(51),(52)と順次計算をおこ
ない(S3 )、次に、k1 とn2 を1増やし(S4 )、
次に、もし、k1 mod N2 =1であるかの判定をし
(C2 )、判定C2 が真であれば、近似的フィルタ係数
Zの修正(S5 ) Z(k2 +N2 )=Z(k2 )+XL,N2(k2 +N2 −
p)SN2+p-1(k+N2 −1)p-N2+p-1 を行い、k2 をN2 だけ増やし、n2 を0にする
(S6 )。次に判定C1に戻る。
刻の指標であり、k1 はその他の変数の時刻の指標であ
る。また、n2 =k1 mod N2 −1なる値をもつ変
数である。図1はこの発明法の計算手順を示している。
まず、初期化として、 Rp+N2-1(0)=XL (0) T [XL (-1), XL (-2), …, XL (-p-N2+1) ]T Sp-1 (0)=[0,0,…,0]T ZL-1 (1)=[0,0,…,0]T Ep-1 (0)=[y(0),y(−1),…,y(−p+2)]T とし(S0 )、次にk1=k2=1,n 2 =0とし(S
1 )、次に、もし、k1mod N1 =1であるかの判
定をし(C1 )、判定C1 が真であれば、入力ベクトル
と近似的フィルタ係数との内積演算 y′(k1 +i)=XL (k2 +i)T Z(k2 ), i=0,…,N1 −1 (S2 )を行う。その後、式(47),(48),(4
9),(50),(51),(52)と順次計算をおこ
ない(S3 )、次に、k1 とn2 を1増やし(S4 )、
次に、もし、k1 mod N2 =1であるかの判定をし
(C2 )、判定C2 が真であれば、近似的フィルタ係数
Zの修正(S5 ) Z(k2 +N2 )=Z(k2 )+XL,N2(k2 +N2 −
p)SN2+p-1(k+N2 −1)p-N2+p-1 を行い、k2 をN2 だけ増やし、n2 を0にする
(S6 )。次に判定C1に戻る。
【0052】この手順に対応する機能ブロック図は図2
に示し、図7A、図7Bと対応する部分と同一符号を付
けて示す。図2のユニットは、従来法の機能ブロック図
である図11中のユニットと機能は同じだが動作の順序
が異なっている。この発明の演算量を見積る。まず、同
方法の手順を示した図1のS2 は、従来法で説明したよ
うに、FFTを用いた高速畳み込みが利用でき、その演
算量は、約T0 =8(M1 +2)N1 log2 2N1 +
8N1 M1 であり、1時刻当たりでは、T0 /N1 =8
(M1 +2)log2 2N1 +8M1 である。ここで、
L=N1 M1 である。S3 では、1時刻当たり、式(4
7)に2(N2 +p−2)回、式(48)にp+N2 /
2回、式(50)にp回、式(51)には15p回、式
(52)には、p回の演算が必要であり、合計すると
2.5N2 +20pとなる(Tc =N2 (2.5N2 +
20p)と定義しておく)。S5 に要する演算量は従来
法で説明したように、L=N2 M2 とすれば、約Tu =
8(M2 +2)N2 log2 2N2 +8N2 M2 であ
り、1時刻当たりでは、Tu /N2 =8(M2 +2)l
og2 2N2 +8M2 である。以上より、この発明法の
実行に要する演算量をTとかくと、TはS2 、S3 、そ
してS5 に要する1時刻当たりの演算量の和として、 T=T0 /N1 +TC /N2 +Tu /N2 =8(M1 +2)log2 (2N1 ) +8(M2 +2)log2 (2N2 ) +2.5N2 +8M1 +8M2 +20p (54) と、見積られる。図3に式(42)によって計算した演
算量Tとブロック長N1,N2 との関係の例を示した。
条件は、フィルタ長はL=4096、ブロック長N1 ,
N2 は128,256,512,1024,2048で
ある。同図で、ΔはN1 =N2 なる従来法により可能な
演算量であり、ブロック長が512で演算量が最低で、
ブロック長がそれより短くても長くても演算量が増加す
る。一方、□はN1 ≠N2 なるこの発明法による演算量
であり、ブロック長N1 を長くするとともに、演算量は
低減し、ブロック長N1 =1024,2048では、等
ブロック長の従来法の場合の最小演算量より少なくする
ことが可能である。N2 を変化させた場合、N2 =25
6で何れのN1 に対しても最小となる。図3より、この
発明が従来法に比べて演算量が少なくなるのは、ブロッ
ク長N1 が、式(45)で計算される従来法の最小演算
量を与えるブロック長Nより長い場合である。L=40
96のこの例では、式(45)の再帰計算をN(0)=
512なる初期値で開始し、N=487を得る。
に示し、図7A、図7Bと対応する部分と同一符号を付
けて示す。図2のユニットは、従来法の機能ブロック図
である図11中のユニットと機能は同じだが動作の順序
が異なっている。この発明の演算量を見積る。まず、同
方法の手順を示した図1のS2 は、従来法で説明したよ
うに、FFTを用いた高速畳み込みが利用でき、その演
算量は、約T0 =8(M1 +2)N1 log2 2N1 +
8N1 M1 であり、1時刻当たりでは、T0 /N1 =8
(M1 +2)log2 2N1 +8M1 である。ここで、
L=N1 M1 である。S3 では、1時刻当たり、式(4
7)に2(N2 +p−2)回、式(48)にp+N2 /
2回、式(50)にp回、式(51)には15p回、式
(52)には、p回の演算が必要であり、合計すると
2.5N2 +20pとなる(Tc =N2 (2.5N2 +
20p)と定義しておく)。S5 に要する演算量は従来
法で説明したように、L=N2 M2 とすれば、約Tu =
8(M2 +2)N2 log2 2N2 +8N2 M2 であ
り、1時刻当たりでは、Tu /N2 =8(M2 +2)l
og2 2N2 +8M2 である。以上より、この発明法の
実行に要する演算量をTとかくと、TはS2 、S3 、そ
してS5 に要する1時刻当たりの演算量の和として、 T=T0 /N1 +TC /N2 +Tu /N2 =8(M1 +2)log2 (2N1 ) +8(M2 +2)log2 (2N2 ) +2.5N2 +8M1 +8M2 +20p (54) と、見積られる。図3に式(42)によって計算した演
算量Tとブロック長N1,N2 との関係の例を示した。
条件は、フィルタ長はL=4096、ブロック長N1 ,
N2 は128,256,512,1024,2048で
ある。同図で、ΔはN1 =N2 なる従来法により可能な
演算量であり、ブロック長が512で演算量が最低で、
ブロック長がそれより短くても長くても演算量が増加す
る。一方、□はN1 ≠N2 なるこの発明法による演算量
であり、ブロック長N1 を長くするとともに、演算量は
低減し、ブロック長N1 =1024,2048では、等
ブロック長の従来法の場合の最小演算量より少なくする
ことが可能である。N2 を変化させた場合、N2 =25
6で何れのN1 に対しても最小となる。図3より、この
発明が従来法に比べて演算量が少なくなるのは、ブロッ
ク長N1 が、式(45)で計算される従来法の最小演算
量を与えるブロック長Nより長い場合である。L=40
96のこの例では、式(45)の再帰計算をN(0)=
512なる初期値で開始し、N=487を得る。
【0053】また、通常は演算量が規定され、その範囲
で可能なだけ長いフィルタ長を実装するので、フィルタ
長を演算量の関数として表わす方が便利なこともある。
式(42)にL=N1 M1 、L=N2 M2 を代入しフィ
ルタ長Lについて解くと次式を得る。 L=[T-16log2(2N1)-16log2(2N2)-2.5N2-20p ] /[8(log2(2N1)+1)/N1+8(log2(2N2)+1)/N2] (55) 図4に式(55)によって計算した、フィルタ長Lとブ
ロック長N1 ,N2 との関係の例を示した。条件は演算
量はT=8000、ブロック長N1 ,N2 は128,2
56,512,1024,2048,4096である。
同図で、ΔはN1=N2 なる従来法により可能なフィル
タ長であり、ブロック長N1 が1024でフィルタ長が
最長で、ブロック長がそれより短くても長くても演算量
が増加する。特に、ブロック長N1 が4096以上では
L>0となる解がない。一方、□はN1 ≠N2 なるこの
発明法によるフィルタ長であり、ブロック長N1 を長く
するとともに、フィルタ長が長くなり、ブロック長N1
=2048以上では、等ブロック長の従来法の場合の最
長フィルタ長より長いフィルタ長を実装することが可能
である。
で可能なだけ長いフィルタ長を実装するので、フィルタ
長を演算量の関数として表わす方が便利なこともある。
式(42)にL=N1 M1 、L=N2 M2 を代入しフィ
ルタ長Lについて解くと次式を得る。 L=[T-16log2(2N1)-16log2(2N2)-2.5N2-20p ] /[8(log2(2N1)+1)/N1+8(log2(2N2)+1)/N2] (55) 図4に式(55)によって計算した、フィルタ長Lとブ
ロック長N1 ,N2 との関係の例を示した。条件は演算
量はT=8000、ブロック長N1 ,N2 は128,2
56,512,1024,2048,4096である。
同図で、ΔはN1=N2 なる従来法により可能なフィル
タ長であり、ブロック長N1 が1024でフィルタ長が
最長で、ブロック長がそれより短くても長くても演算量
が増加する。特に、ブロック長N1 が4096以上では
L>0となる解がない。一方、□はN1 ≠N2 なるこの
発明法によるフィルタ長であり、ブロック長N1 を長く
するとともに、フィルタ長が長くなり、ブロック長N1
=2048以上では、等ブロック長の従来法の場合の最
長フィルタ長より長いフィルタ長を実装することが可能
である。
【0054】つぎに、ブロック処理を行うことによるフ
ィルタ出力y(k)が遅延する時間を考える。この遅延
時間をDと表わすと、図5より、次式で見積られること
がわかる。 D=N1 +T0 /T ,N1 >N2 =(2T0 +Tc +Tu )/T ,N 1 <N2 (56) 図5Cに式(56)によって計算した出力遅延時間Dと
ブロック長N2 との関係の例を示した。条件はフィルタ
長L=4096、ブロック長N1 ,N2 は128,25
6,512,1024,2048である。同図で、Δは
N1 =N2 なる従来法による出力遅延時間であり、□は
N1 ≠N2 なるこの発明法による出力遅延時間であり、
ブロック長N1 は各々のN2 について、下から128,
256,512,1024,2048である。N2 が2
56以上では、この発明法は従来法より短い出力遅延時
間を実現できることがわかる。
ィルタ出力y(k)が遅延する時間を考える。この遅延
時間をDと表わすと、図5より、次式で見積られること
がわかる。 D=N1 +T0 /T ,N1 >N2 =(2T0 +Tc +Tu )/T ,N 1 <N2 (56) 図5Cに式(56)によって計算した出力遅延時間Dと
ブロック長N2 との関係の例を示した。条件はフィルタ
長L=4096、ブロック長N1 ,N2 は128,25
6,512,1024,2048である。同図で、Δは
N1 =N2 なる従来法による出力遅延時間であり、□は
N1 ≠N2 なるこの発明法による出力遅延時間であり、
ブロック長N1 は各々のN2 について、下から128,
256,512,1024,2048である。N2 が2
56以上では、この発明法は従来法より短い出力遅延時
間を実現できることがわかる。
【0055】以上説明したこの発明法の特徴は、N1 ≠
N2 とすることにより、従来の毎時刻フィルタを修正す
る射影法と同じ出力をもつブロック射影法より少ない演
算量で実行できる、いいかえると、同じ演算量でより長
いフィルタ長を使用でき、その結果、誤差信号のパワー
がより減少することである。また、出力遅延時間を従来
法より短くすることもでき、この場合はN2 >N1 であ
って、フィルタ係数の更新周期が長くなるため、途中で
エコーキャンセラの場合はダブルトークなどの異常を検
出したら、その時の更新を中止する処理も可能となる。
N2 とすることにより、従来の毎時刻フィルタを修正す
る射影法と同じ出力をもつブロック射影法より少ない演
算量で実行できる、いいかえると、同じ演算量でより長
いフィルタ長を使用でき、その結果、誤差信号のパワー
がより減少することである。また、出力遅延時間を従来
法より短くすることもでき、この場合はN2 >N1 であ
って、フィルタ係数の更新周期が長くなるため、途中で
エコーキャンセラの場合はダブルトークなどの異常を検
出したら、その時の更新を中止する処理も可能となる。
【0056】実施例
第1の例は、音響エコーキャンセラである。図12Aは
この発明の方法を適応した音響エコーキャンセラの機能
構成を示す。入力信号x(k)として音声信号がフィル
タリング部11、フィルタ係数更新部12、スピーカ6
1へ供給され、マイクロホン62の出力が所望信号d
(k)として差回路13へ供給される。スピーカ61か
ら音響信号通路、いわゆるエコー伝達経路63を通じて
マイクロホン62へ達する信号が音響エコーであり、適
応フィルタ(フィルタリング部)11の出力信号y
(k)に推定エコー、推定誤差e(k)に残留エコーが
対応している。
この発明の方法を適応した音響エコーキャンセラの機能
構成を示す。入力信号x(k)として音声信号がフィル
タリング部11、フィルタ係数更新部12、スピーカ6
1へ供給され、マイクロホン62の出力が所望信号d
(k)として差回路13へ供給される。スピーカ61か
ら音響信号通路、いわゆるエコー伝達経路63を通じて
マイクロホン62へ達する信号が音響エコーであり、適
応フィルタ(フィルタリング部)11の出力信号y
(k)に推定エコー、推定誤差e(k)に残留エコーが
対応している。
【0057】音響エコーキャンセラは、テレビ会議シス
テム等の拡声通話系において、ハウリングや、エコーを
防止する装置である。音響エコーキャンセラの場合、適
応フィルタ11は、スピーカ61からマイクロホン62
の間の伝達経路63を模擬する。適応フィルタ11の出
力は推定されたエコーであり、マイクロホン62で受け
たエコーから差し引くことでエコーが除去される。スピ
ーカ61からマイクロホン62の間の伝達経路63を適
応フィルタ11で模擬する場合に必要なフィルタの長さ
は数百から数千にもなり、多くの演算量が必要となる。
また、スピーカ61からマイクロホン62の間の伝達経
路63が人の動き、ドアの開閉等で変化するので、その
変化に追従できる高速な適応アルゴリズムが必要とされ
ている。
テム等の拡声通話系において、ハウリングや、エコーを
防止する装置である。音響エコーキャンセラの場合、適
応フィルタ11は、スピーカ61からマイクロホン62
の間の伝達経路63を模擬する。適応フィルタ11の出
力は推定されたエコーであり、マイクロホン62で受け
たエコーから差し引くことでエコーが除去される。スピ
ーカ61からマイクロホン62の間の伝達経路63を適
応フィルタ11で模擬する場合に必要なフィルタの長さ
は数百から数千にもなり、多くの演算量が必要となる。
また、スピーカ61からマイクロホン62の間の伝達経
路63が人の動き、ドアの開閉等で変化するので、その
変化に追従できる高速な適応アルゴリズムが必要とされ
ている。
【0058】第2の適用例は騒音制御である。図12B
に騒音制御の原理を示す。騒音制御の目的は、騒音源8
1から騒音伝達経路82を径て観測される騒音を、スピ
ーカ86から負の音(音圧がy(k)と表される時に、
−y(k)と表される音圧を、y(k)に対する負の音
と呼ぶ)を出すことにより消去しようとするものであ
る。この例では、適応フィルタ(フィルタリング部)1
1への入力信号x(k)は騒音源81の騒音をモニタマ
イクロホン84で観測した騒音であり、所望信号d
(k)は、観測地点(マイクロホン83)における騒音
であり、適応フィルタ11の出力y(k)は推定騒音で
あり、誤差信号e(k)は観測地点で、スピーカ86か
ら放射された推定騒音と騒音源から到来した騒音の音圧
が合成され、マイクロホン83の出力に得られる。適応
フィルタ11が、騒音源81から観測マイクロホン83
までの伝達経路82を良好に模擬していれば、スピーカ
86から放射される音圧−y(k)によって騒音d
(k)は消去される。
に騒音制御の原理を示す。騒音制御の目的は、騒音源8
1から騒音伝達経路82を径て観測される騒音を、スピ
ーカ86から負の音(音圧がy(k)と表される時に、
−y(k)と表される音圧を、y(k)に対する負の音
と呼ぶ)を出すことにより消去しようとするものであ
る。この例では、適応フィルタ(フィルタリング部)1
1への入力信号x(k)は騒音源81の騒音をモニタマ
イクロホン84で観測した騒音であり、所望信号d
(k)は、観測地点(マイクロホン83)における騒音
であり、適応フィルタ11の出力y(k)は推定騒音で
あり、誤差信号e(k)は観測地点で、スピーカ86か
ら放射された推定騒音と騒音源から到来した騒音の音圧
が合成され、マイクロホン83の出力に得られる。適応
フィルタ11が、騒音源81から観測マイクロホン83
までの伝達経路82を良好に模擬していれば、スピーカ
86から放射される音圧−y(k)によって騒音d
(k)は消去される。
【0059】以上、2つの適用例を示した。どちらの場
合にもエコー、騒音を消去するという目的を達成するた
めには、推定エコー、推定騒音が得られれば十分であ
り、エコー、騒音の伝達経路そのものを求める必要はな
い。従って、どちらの適用例でも近似的フィルタ係数を
使用するこの発明法が適用でき、従来装置と比較して演
算量を削減できる。更に演算量が多くなるが、出力y
(k)が速く得られるようにすることができる。また、
この発明は、図6Aに示したような、構成において誤差
信号e(k)のパワーを最小とする基本機能を有してい
る。従って、解決すべき問題が図6Aに示した誤差信号
のパワーの最小化としてモデル化できる全ての応用例に
対して、この発明を適用することが可能である。
合にもエコー、騒音を消去するという目的を達成するた
めには、推定エコー、推定騒音が得られれば十分であ
り、エコー、騒音の伝達経路そのものを求める必要はな
い。従って、どちらの適用例でも近似的フィルタ係数を
使用するこの発明法が適用でき、従来装置と比較して演
算量を削減できる。更に演算量が多くなるが、出力y
(k)が速く得られるようにすることができる。また、
この発明は、図6Aに示したような、構成において誤差
信号e(k)のパワーを最小とする基本機能を有してい
る。従って、解決すべき問題が図6Aに示した誤差信号
のパワーの最小化としてモデル化できる全ての応用例に
対して、この発明を適用することが可能である。
【0060】
【発明の効果】以上説明したように、この発明では、適
応アルゴリズムの1つである毎時刻フィルタを修正する
射影法と同じ出力をもつブロック射影法において2つの
相違なるブロック長を導入し、動作手順を変更すること
により、収束速度を低下させること無く、演算量低減を
可能とする。その結果、この発明は収束の遅い従来のブ
ロック射影法等の手法が利用されていた音響エコーキャ
ンセラ装置に対して、収束速度を低下させることなく、
演算量の低減が可能となる。即ち、同一のハードウェア
規模で収束は速いが多少演算量の必要なより高い次数の
射影法の利用や、フィルタ長を長くすることが可能とな
り、より高速なエコー低減や定常時のエコーの低減の実
現が達成される。またこの発明によれば、演算量は低減
できないが出力y(k)が得られる遅延時間を短かくす
ることができ、その場合はフィルタの更新時間が長くな
り、この間にダブルトークなどの異常の発生を検出処理
を行い、検出すればフィルタの次の更新を中止するなど
の処理を行うことができる。
応アルゴリズムの1つである毎時刻フィルタを修正する
射影法と同じ出力をもつブロック射影法において2つの
相違なるブロック長を導入し、動作手順を変更すること
により、収束速度を低下させること無く、演算量低減を
可能とする。その結果、この発明は収束の遅い従来のブ
ロック射影法等の手法が利用されていた音響エコーキャ
ンセラ装置に対して、収束速度を低下させることなく、
演算量の低減が可能となる。即ち、同一のハードウェア
規模で収束は速いが多少演算量の必要なより高い次数の
射影法の利用や、フィルタ長を長くすることが可能とな
り、より高速なエコー低減や定常時のエコーの低減の実
現が達成される。またこの発明によれば、演算量は低減
できないが出力y(k)が得られる遅延時間を短かくす
ることができ、その場合はフィルタの更新時間が長くな
り、この間にダブルトークなどの異常の発生を検出処理
を行い、検出すればフィルタの次の更新を中止するなど
の処理を行うことができる。
【図1】この発明の実施例を示す流れ図。
【図2】この発明の実施例を適用した装置の機能構成を
示す図。
示す図。
【図3】フィルタ長を一定とした時のブロック長N1 ,
N2 と演算量Tの関係を示す図。
N2 と演算量Tの関係を示す図。
【図4】演算量を一定とした時のブロック長N1 ,N2
とフィルタ長Lとの関係を示す図。
とフィルタ長Lとの関係を示す図。
【図5】A及びBはそれぞれ演算実行時間と出力遅延時
間とを示す図、Cはブロック長N1 ,N2 と出力遅延時
間Dとの関係を示す図である。
間とを示す図、Cはブロック長N1 ,N2 と出力遅延時
間Dとの関係を示す図である。
【図6】Aは適応フィルタを用いた適応制御系を示す
図、Bは従来の射影法を用いた制御の機能構成を示す図
である。
図、Bは従来の射影法を用いた制御の機能構成を示す図
である。
【図7】Aは従来の高速射影法を用いた制御の機能構成
を示す図、Bは従来のブロック射影法を用いた制御の機
能構成を示す図である。
を示す図、Bは従来のブロック射影法を用いた制御の機
能構成を示す図である。
【図8】従来の高速射影法を用いた制御の計算手順を示
す流れ図。
す流れ図。
【図9】従来のブロック法を用いた制御の計算手順を示
す流れ図。
す流れ図。
【図10】従来の毎時刻フィルタを修正する射影法と同
じ出力をもつブロック射影法を用いた制御の機能構成を
示す図。
じ出力をもつブロック射影法を用いた制御の機能構成を
示す図。
【図11】従来の毎時刻フィルタを修正する射影法と同
じ出力をもつブロック射影法を用いた制御の計算手順を
示す流れ図。
じ出力をもつブロック射影法を用いた制御の計算手順を
示す流れ図。
【図12】Aはこの発明を音響エコーキャンセラに適用
した機能構成を示す図、Bはこの発明を騒音制御に適用
した機能構成を示す図である。
した機能構成を示す図、Bはこの発明を騒音制御に適用
した機能構成を示す図である。
─────────────────────────────────────────────────────
フロントページの続き
(56)参考文献 特開 平9−54590(JP,A)
特開 昭59−139717(JP,A)
(58)調査した分野(Int.Cl.7,DB名)
H03H 17/00 - 17/08
H03H 21/00
H04B 3/20 - 3/23
Claims (1)
- 【請求項1】 入力信号を可変フィルタ処理し、その出
力信号と所望信号との誤差信号を求め、その誤差信号の
パワーを最小とするように、上記可変フィルタの特性を
適応的に制御する方法において、 上記入力信号ベクトルの自己相関行列の逆行列と、上記
誤差信号ベクトルとの積に相当する演算を行ってプレフ
ィルタ係数ベクトルを求める第1ステップと、 上記プレフィルタ係数ベクトルを平滑化して平滑係数ベ
クトルを得る第2ステップと、 上記入力信号ベクトルを列要素に持つ行列と上記平滑係
数ベクトルの部分ベクトルの積を前回の近似的フィルタ
係数ベクトルに加えることで近似的フィルタ係数ベクト
ルをブロック処理により更新修正する第3ステップと、 上記入力信号ベクトルと上記近似的フィルタ係数ベクト
ルとの内積演算をブロック処理により行う第4ステップ
と、 上記入力信号ベクトルの自己相関ベクトルを求める第5
ステップと、 上記平滑係数ベクトルと上記自己相関ベクトルとの内積
演算を行う第6ステップと、 上記第4ステップで得られた内積演算結果と、上記第6
ステップで得られた内積演算結果とを加算して上記出力
信号を得る第7ステップとを有するものであって、 上記第3ステップの修正及び上記第4ステップのブロッ
ク処理による演算はそれぞれN1 、N2 時刻(N1 、N
2 は2以上の整数、N1 ≠N2 )ごとに行い、 上記第3、第4ステップ以外の各ステップは各時刻ごと
に演算を行い、 上記第3ステップの演算が行われるごとに、次に上記第
3ステップの演算されるまでの各時刻において、上記第
2ステップで得る上記平滑係数ベクトルの次数を1ずつ
増加し、次の第3ステップの演算時刻での演算でもとの
次数に戻し、 上記第5ステップでは上記平滑係数ベクトルの最大次数
と等しい次数の上記自己相関ベクトルを演算し、 上記第6ステップでは上記平滑係数ベクトルの次数と等
しい次数だけ上記第5ステップで得られた自己相関ベク
トルから用いることを特徴とする適応的制御方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP5048497A JP3475318B2 (ja) | 1997-03-05 | 1997-03-05 | 適応的制御方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP5048497A JP3475318B2 (ja) | 1997-03-05 | 1997-03-05 | 適応的制御方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
JPH10247838A JPH10247838A (ja) | 1998-09-14 |
JP3475318B2 true JP3475318B2 (ja) | 2003-12-08 |
Family
ID=12860199
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP5048497A Expired - Fee Related JP3475318B2 (ja) | 1997-03-05 | 1997-03-05 | 適応的制御方法 |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP3475318B2 (ja) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3536921B2 (ja) | 2001-04-18 | 2004-06-14 | 日本電気株式会社 | 相関行列学習方法、装置及びプログラム |
JP4582680B2 (ja) * | 2001-05-01 | 2010-11-17 | 国立大学法人東京工業大学 | 係数算出支援装置、係数算出装置、これらの係数算出支援装置または係数算出装置を実現するプログラムおよびそのプログラムが格納された記録媒体 |
JP4687948B2 (ja) * | 2004-10-29 | 2011-05-25 | ソニー株式会社 | ディジタル信号処理装置、ディジタル信号処理方法及びプログラム並びに認証装置 |
US11770205B2 (en) * | 2019-09-06 | 2023-09-26 | Nippon Telegraph And Telephone Corporation | Optical signal processing apparatus, optical signal processing method and computer program |
-
1997
- 1997-03-05 JP JP5048497A patent/JP3475318B2/ja not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
JPH10247838A (ja) | 1998-09-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP4398146B2 (ja) | 適応フィルタ | |
EP2327156B1 (en) | Method for determining updated filter coefficients of an adaptive filter adapted by an lms algorithm with pre-whitening | |
CN108172231B (zh) | 一种基于卡尔曼滤波的去混响方法及*** | |
US8385864B2 (en) | Method and device for low delay processing | |
US7774396B2 (en) | Method and device for low delay processing | |
US5774562A (en) | Method and apparatus for dereverberation | |
JP4700871B2 (ja) | 音響エコー及びノイズ除去 | |
US4951269A (en) | Echo canceller with short processing delay and decreased multiplication number | |
JP2000035788A (ja) | 多重チャネル適応フィルタリング | |
JP3475318B2 (ja) | 適応的制御方法 | |
Bourgeois et al. | Time-domain beamforming and blind source separation: speech input in the car environment | |
Jin et al. | A simultaneous equation method-based online secondary path modeling algorithm for active noise control | |
Milani et al. | Analysis and optimal design of delayless subband active noise control systems for broadband noise | |
EP0637803A2 (en) | Method and device for adaptively estimating a transfer function of an unknown system | |
Wu et al. | Meta-learning for adaptive filters with higher-order frequency dependencies | |
Althahab | A new robust adaptive algorithm based adaptive filtering for noise cancellation | |
JP2541044B2 (ja) | 適応フィルタ装置 | |
JP3452339B2 (ja) | 適応的制御方法 | |
JP4041770B2 (ja) | 音響エコー消去方法、その装置、プログラム及びその記録媒体 | |
JPH09261135A (ja) | 音響エコー消去装置 | |
JP3435675B2 (ja) | 適応的制御方法 | |
JP3303898B2 (ja) | 適応的伝達関数推定方法及びそれを使った推定装置 | |
KR20050047374A (ko) | 통신 기기용 하이브리드 소음 제거 시스템 및 방법 | |
Gay | Affine projection algorithms | |
JPH09232913A (ja) | 適応フィルタの係数更新装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080926 Year of fee payment: 5 |
|
FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080926 Year of fee payment: 5 |
|
FPAY | Renewal fee payment (prs date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090926 Year of fee payment: 6 |
|
LAPS | Cancellation because of no payment of annual fees |