JP3475318B2 - 適応的制御方法 - Google Patents

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】この発明は、音響エコーキャ
ンセラ、能動騒音制御等において、適応フィルタの推定
誤差のパワーを最小にする適応制御方法に関する。

【０００２】

【従来の技術】この節では、まず、適応的制御に用いら
れる、適応フィルタと適応アルゴリズムについて説明す
る。次に、その適応フィルタを修正するための適応アル
ゴリズムの中の射影法を説明し、さらに、射影法の演算
量低減法として、演算手順を工夫した高速射影法と適応
フィルタ修正を毎時刻行わないブロック処理に基づく射
影法について説明する。最後に、ブロック処理に基づく
射影法の問題点を述べる。 ３．２．１ 適応フィルタと適応アルゴリズム 図６Ａは適応フィルタを用いた適応制御系のブロック図
を示す。図において、入力信号をｘ（ｋ）、出力信号を
ｙ（ｋ）、所望信号をｄ（ｋ）、誤差信号をｅ（ｋ）で
表わす。ただし、この明細書では、時間の表現は離散時
間として整数ｋで時間を表すものとする。また、適応フ
ィルタはＦＩＲ形フィルタであって、時変係数｛ｈ
₁ （ｋ），ｈ₂ （ｋ），…，ｈ_L （ｋ）｝を持つものと
する。ただし、Ｌはフィルタ係数の数である。また、フ
ィルタ係数はベクトルとして、 Ｈ（ｋ）＝［ｈ₁ （ｋ），ｈ₂ （ｋ），…，ｈ_L （ｋ）］^T (1) と表す。ただし、［ ］^T はベクトルまたは行列の転置
を表す。適応フィルタは、誤差信号ｅ（ｋ）と入力信号
ｘ（ｋ）に基づいて、フィルタ係数Ｈ（ｋ）を調整し
て、誤差信号ｅ（ｋ）のパワーを最小にする。この時の
係数修正の手順は適応アルゴリズムとして知られてい
る。

【０００３】次に、適応アルゴリズムについて説明す
る。適応アルゴリズムとは、例えば図６Ａに示すように
入力信号ｘ（ｋ）がフィルタリング部１１に供給され、
フィルタリング部１１の出力信号ｙ（ｋ）と所望信号ｄ
（ｋ）との誤差信号ｅ（ｋ）が差回路１３でとられ、誤
差信号ｅ（ｋ）と入力信号ｘ（ｋ）とがフィルタ係数更
新部１２に入力され、フィルタリング部１１のフィルタ
係数Ｈ（ｋ＋１）の更新がなされる。適応アルゴリズム
はこのような適応フィルタの内部の演算手順であって、
時刻ｋにおいて、フィルタ係数Ｈ（ｋ）をフィルタ係数
更新部１２で修正して、誤差信号ｅ（ｋ）のパワーがな
るべく小さくなるような係数Ｈ（ｋ＋１）を得る方法で
ある。その修正は、次式で表される。

【０００４】 Ｈ（ｋ＋１）＝Ｈ（ｋ）＋μΔＨ（ｋ） (2.1) ここで、ΔＨ（ｋ）は修正ベクトル、μは修正の大きさ
を制御する量でステップサイズと呼ばれる。現在、最も
多く利用されている適応アルゴリズムとしては、ＬＭＳ
法や学習同定法などが知られている。ＬＭＳ法において
は、修正ベクトルは誤差信号ｅ（ｋ）と入力ベクトルｘ
（ｋ）との積として、 ΔＨ（ｋ）＝ｅ（ｋ）Ｘ_L （ｋ） (2.2) と表され、修正式は、 Ｈ（ｋ＋１）＝Ｈ（ｋ）＋μｅ（ｋ）Ｘ_L （ｋ） (2.3) となる。ただし、 Ｘ_L （ｋ）＝［ｘ（ｋ），ｘ（ｋ−１），…，ｘ（ｋ−Ｌ＋１）］^T (3.1) で、これを入力（信号）ベクトルまたは単に、入力と呼
ぶ。また、次の時刻ｋ＋１のフィルタ出力ｙ（ｋ＋１）
は、図６Ａのフィルタリング部１１で次式による内積演
算により合成される。

【０００５】 ｙ（ｋ＋１）＝Ｘ_L （ｋ＋１）^T Ｈ（ｋ＋１） (3.2) 式（２．３）に基づいてフィルタ係数を修正して、式
（３．２）に基づいて次の時刻の出力を合成するという
適応フィルタの動作に必要な積和演算量は２Ｌ（Ｌはフ
ィルタ係数の数）である。この２Ｌという演算量は、現
在知られている毎時刻フィルタ係数を修正する適応アル
ゴリズムのなかでは最小の演算量となっている。 ３．２．２ 射影法 このように、ＬＭＳ法や学習同定法は演算量が少ないと
いう長所を持っている。しかし、その反面、有色信号や
周期性雑音などに対しては、収束速度が遅い（適切なフ
ィルタ係数を得るまでの時間がかかる）という欠点をも
っている。

【０００６】ＬＭＳ法や学習同定法の欠点を克服する方
法として、１９８４年に射影法が提案されたが、射影法
では演算量が多くなるという欠点をもっていた。しかし
近年、この射影法の欠点を改善して、低演算量で実行す
る高速射影法が提案されている。ここではその概要を以
下に説明する。射影法において、誤差を小さくするため
には、フィルタ係数の修正は、次のｐ個の式 ｄ（ｋ）＝Ｘ_L （ｋ）^T Ｈ（ｋ＋１） ｄ（ｋ＋１）＝Ｘ_L （ｋ＋１）^T Ｈ（ｋ＋１） (4) … ｄ（ｋ−ｐ＋１）＝Ｘ_L （ｋ−ｐ＋１）^T Ｈ（ｋ＋１） を満たすように行われる。ただし、ｐ（ｐ＜Ｌ）は射影
の次数と呼ばれる。この式（４）の意味は、修正を行っ
た後の係数Ｈ（ｋ＋１）を過去の入力Ｘ_L （ｋ−ｉ＋
１）^T と畳み込んだ結果Ｘ_L （ｋ−ｉ＋１）^T Ｈ（ｋ＋
１）が、過去の所望信号ｄ（ｋ−ｉ＋１）と一致する、
という関係が過去ｐ個の入力に対して成立するように、
Ｈ（ｋ＋１）を決定するというものである。このよう
に、過去の入力に対して所望信号と同一の値を出力する
ようにフィルタ係数を定めれば、未来の時刻の入力Ｘ_L
（ｋ＋１）に対する出力ｙ（ｋ＋１）も所望信号ｄ（ｋ
＋１）に近い値が出力できるものと考えられる。その結
果、未来の誤差信号ｅ（ｋ＋１）＝ｄ（ｋ＋１）−ｙ
（ｋ＋１）も小さな値となるものと考えられる。さて、
式（２．１）を式（４）に代入してΔＨ（ｋ）について
解けば、 ΔＨ（ｋ）＝Ｘ_L,p (k) ［Ｘ_L,p (k) ^T Ｘ_L,p (k) ］^-1Ｅ_p (k) (5) と表される。ただし、 Ｘ_L,p (k) ＝［Ｘ_L (k),Ｘ_L (k-1),…, Ｘ_L (k-p+1) ］ (6) Ｅ_p (k) ＝［d(k),d(k-1),…,d(k-p+1) ］^T −Ｘ_L,p (k) ^T Ｈ(k) (7.1) である。Ｘの下付き添え字Ｌ，ｐはＸがＬ行ｐ列である
ことを示す。また、このＥ_p （ｋ）は、次式のように、
再帰的にも表すことができる。

【０００７】 ただし、 Ｅ_p-1 （ｋ−１）＝［ｅ（ｋ−１），（１−μ）ｅ（ｋ−２），…， （１−μ）^p-2 ｅ（ｋ−ｐ＋１）］^T (7.3) である。

【０００８】式（５），（２．１）より、入力Ｘ
_L （ｋ）と誤差信号ｅ（ｋ）を用いてフィルタ係数を修
正して、係数Ｈ（ｋ＋１）を求め、これを用いて時刻ｋ
＋１の出力ｙ（ｋ＋１）を合成する手順は、次式のよう
に表すことができる。 Ｇ（ｋ）＝［Ｘ_L,p (k) ^T Ｘ_L,p (k) ］^-1Ｅ_p (k) (8) Ｈ（ｋ＋１）＝Ｈ（ｋ）＋μＸ_Lp（ｋ）Ｇ（ｋ） (9) ｙ（ｋ＋１）＝Ｘ_L （ｋ＋１）^T Ｈ（ｋ＋１） (10) ただし、Ｇ（ｋ）は式（５）の後半を表し、プレフィル
タ係数と呼ばれる。この式（８），（９），（１０）が
従来の射影法の原理的手順を表す式であり、これを機能
ブロック図で表したものを図６Ｂに示す。図６Ｂにおい
てプレフィルタ係数計算部２１で入力信号ｘ（ｋ）と誤
差信号ｅ（ｋ）とから式（８）の演算を行い、得られた
プレフィルタ係数Ｇ（ｋ）と入力信号ｘ（ｋ）とによ
り、フィルタ係数更新部２２で式（９）の演算を行いフ
ィルタ係数を修正し、その修正されたフィルタ係数と入
力信号とによりフィルタリング部１１で式（１０）の演
算を実行する。

【０００９】ここで、この系の演算量を考えてみる。ほ
とんどのＤＳＰでは、加算と積和演算がそれぞれ１クロ
ックで実行され、また、加算と積和演算が計算の最も大
きな部分を占めるので、ここでは、演算量として加算と
積和演算の数として評価する。まず、プレフィルタ係数
計算部２１で行う式（８）の演算は、ｐ行ｐ列の行列Ｘ
_L,P （ｋ）^T Ｘ_L,P （ｋ）の逆行列の計算を含んでいる
ので、ｐの３乗に比例する積和演算量（以降、Ｏ
（ｐ³ ）と表す）が必要となる。これは、ｐが大きくな
ると、たいへん大きな演算量となる。次に、フィルタ係
数更新部２２で行う式（９）の演算はＸ_L,p （ｋ）がＬ
行ｐ列の行列、Ｇ（ｋ）がｐ次のベクトルであることよ
り、ｐＬ回の積和演算が必要である。更に、フィルタ実
行部１１で行う式（１０）の演算はＬ次のベクトルの内
積であるので、Ｌ回の積和演算が必要である。従って、
合計（ｐ＋１）Ｌ＋Ｏ（ｐ³ ）の積和演算量が必要とな
る。例えば、ｐ＝２０、Ｌ＝１０００の場合には、２
１，０００回以上もの積和演算が必要となる。

【００１０】次に、射影法の演算量低減法として、演算
を工夫した高速射影法と適応フィルタ更新を毎時刻には
行わず、あるブロック間隔で行うブロック処理に基づく
射影法について説明する。 ３．２．３ 高速射影法 これに対して近年提案された高速射影法では、演算の工
夫により積和演算量を２Ｌ＋２０ｐにまで低減すること
が可能である。以下、高速射影法を簡単に説明する。

【００１１】まず、高速射影法では、式（８）の逆行列
演算を直接行わずに、線形予測法を利用してプレフィル
タ係数Ｇ（ｋ）を計算することで、Ｏ（ｐ³ ）の積和演
算量を１６ｐの積和演算量に低減している。この演算量
の低減に関しては本発明とは直接関連しないので、詳細
は省略する。次に、高速射影法では、平滑化係数Ｓ
_i（ｋ）と近似的フィルタ係数Ｚ（ｋ）を利用すること
で、図６Ｂ中のフィルタ係数更新部２２およびフィルタ
リング部１１で行っている（ｐ＋１）Ｌの積和演算を２
Ｌ＋４ｐの積和演算量にまで低減している。この方法に
ついて図７Ａに基づいて説明を行う。

【００１２】図７Ａにおいて、プレフィルタ係数平滑部
３１は、プレフィルタ係数を平滑化して平滑係数ベクト
ルを得、近似的フィルタ係数更新部３２では、平滑係数
ベクトルの要素を荷重係数として入力信号ベクトルの方
向に近似的フィルタ係数ベクトルを修正し、その近似的
フィルタ係数ベクトルと入力信号ベクトルとの内積をフ
ィルタリング部３６で行い、相関計算部３３で入力信号
ベクトルの自己相関ベクトルを求め、これと平滑係数ベ
クトルとの内積演算が内積計算部３５で行われ、この出
力とフィルタリング部１１の出力とが出力信号合成部３
４で加算されて出力信号とされる。この各構成要素で
は、以下の演算が行われる。即ち、プレフィルタ係数平
滑部３１では式（１１）の演算、近似的フィルタ係数更
新部３２では式（１２）の演算、相関計算部３３では式
（１３）の演算、内積演算計算部３５では、式（１４）
の右辺第２項の内積演算、出力合成部３４では式（１
４）の右辺の第１項と第２項との加算が行われる。な
お、図７Ａに対応する高速射影法の計算手順を図８に示
す。

【００１３】 Ｚ（ｋ＋１）＝Ｚ（ｋ）＋Ｘ_L （ｋ−ｐ＋１）ｓ_p （ｋ） (12) Ｒ_p-1 （ｋ＋１）＝Ｒ_p-1 （ｋ）＋ｘ（ｋ＋１）Ｘ_p-1 （ｋ） −ｘ（ｋ−Ｌ＋１）Ｘ_p-1 （ｋ−Ｌ） (13) ｙ（ｋ＋１）＝Ｘ_L （ｋ＋１）^T Ｚ（ｋ＋１） ＋Ｒ_p-1 （ｋ＋１）^T Ｓ_p-1 （ｋ） (14) ただし、Ｓ_P （ｋ）およびＳ_p-1 （ｋ−１）は、それぞ
れｐ次およびｐ−１次の平滑化係数ベクトルであって、 Ｓ_p （ｋ）＝［ｓ₁ （ｋ），ｓ₂ （ｋ），…，ｓ_p （ｋ）］^T (15) Ｓ_p-1 （ｋ）＝［ｓ₁ （ｋ），ｓ₂ （ｋ），…，ｓ_p-1 （ｋ）］^T (16) と表される。式（１１）からわかるように、Ｓ_p （ｋ）
はプレフィルタ係数ベクトルＧ（ｋ）を用いて再帰的に
計算される。また、Ｚ（ｋ）は、近似的フィルタ係数ベ
クトルであり、Ｘ_L （ｋ−ｐ＋１）とｓ_p （ｋ）を用い
て再帰的に計算される。また、Ｒ_p-1 （ｋ＋１）は、次
式で表される相関ベクトルである。

【００１４】 Ｒ_p-1(k+1)＝［Ｘ_L (k+1) ^T Ｘ_L (k),Ｘ_L (k+1) ^T Ｘ_L (k-1),…, Ｘ_L (k+1) ^T Ｘ_L (k-p-2) ］^T (17) このとき、数学的な証明は省略するが、図７Ａの構成、
または、式（１１）−（１４）の演算により、入力信号
ｘ（ｋ）および誤差信号ｅ（ｋ）を用いて出力信号ｙ
（ｋ＋１）を合成することができる。なお、式（１３）
（１４）の時間パラメータがｋ＋１となっているのに図
７Ａではｋとなっている。これは時刻ｋにおいて式（１
３），（１４），（１１），（１２）の順に計算するこ
とを表したものである。

【００１５】さてここで、以下の２点が指摘できる。ま
ず第１は、以上の方法の演算量に関して、式（１１）で
はｐ、式（１２）ではＬ、式（１３）では２ｐ、式（１
４）ではＬ＋ｐ、合計２Ｌ＋４ｐの積和演算量で実行で
きる。このことは、図６Ｂに示した従来の射影法のフィ
ルタ係数更新部２２およびフィルタリング部１１で行う
式（９）（１０）の演算が（ｐ＋１）Ｌの積和演算を必
要とするのに比べて大幅な低減となっている。例えばｐ
＝２０，Ｌ＝１０００の場合、射影法が２１０００であ
るのに対して、高速射影法では２０８０と、約１／１０
の演算量となっている。

【００１６】第２の点は、図７Ａに示した高速射影法で
は、フィルタ係数Ｈ（ｋ＋１）を陽に求めることなく出
力信号ｙ（ｋ）を合成している点が、図６Ｂに示した従
来の射影法との大きな相違点である。より具体的には、
フィルタ係数Ｈ（ｋ＋１）の代わりに、近似的フィルタ
係数Ｚ（ｋ＋１）を利用し、また、Ｈ（ｋ＋１）とＺ
（ｋ＋１）との差を補正するために、相関ベクトルＲ
_p-1 （ｋ＋１）を利用している。そして、このことが演
算量削減に効果をなしている。 ３．２．４ 従来のブロック射影法 これまで説明してきた高速射影法では、フィルタ係数Ｈ
（ｋ）の代わりに近似的フィルタ係数Ｚ（ｋ）を毎時刻
更新することで、演算量の低減が可能になった。これと
は別に、フィルタ係数Ｈ（ｋ）をＮ（＞１）時刻毎に更
新する（ブロック処理）ことで、演算量を削減すること
が可能である。

【００１７】ブロック処理の手法は、出力ｙ（ｋ）が毎
時刻フィルタ係数を修正する射影法と等しいか否かとい
う観点で２つに分類される。出力ｙ（ｋ）が毎時刻フィ
ルタ係数を修正する射影法と異なる場合は、演算量は削
減されるが、収束速度は遅くなる。一方、毎時刻フィル
タ係数を修正する射影法と同じ出力ｙ（ｋ）をもつブロ
ック処理の場合は、演算量の削減の程度は、毎時刻フィ
ルタ係数を修正する射影法と異なる出力ｙ（ｋ）を持つ
ブロック処理法に比べ、小さくなるが、収束速度の劣化
はみられない利点がある。

【００１８】＜Ａ＞ 出力ｙが毎時刻フィルタ係数を修
正する射影法と異なるブロック射影法 ブロック処理を射影法に適用した従来のブロック射影法
では、Ｅｐ（ｋ），Ｇ（ｋ）の計算、そして、Ｈ（ｋ＋
Ｎ）の更新は各々式（７．１），（８），（９）と同様
に行われる。

【００１９】 Ｅ_p （ｋ＋Ｎ−１）＝ ［ｄ（ｋ＋Ｎ−１），ｄ（ｋ＋Ｎ−２），…，ｄ（ｋ＋Ｎ−ｐ）］^T −Ｘ_L,p （ｋ＋Ｎ−１）^T Ｈ（ｋ） (18) Ｇ（ｋ＋Ｎ−１）＝ ［Ｘ_L,p （ｋ＋Ｎ−１）^T Ｘ_L,p （ｋ＋Ｎ−１）］^-1Ｅ_p （ｋ＋Ｎ−１） (19) Ｈ（ｋ＋Ｎ）＝Ｈ（ｋ）＋μＸ_L,p （ｋ＋Ｎ−１）Ｇ（ｋ＋Ｎ−１） (20) 出力ｙの計算は式（１０）による１時刻の出力ｙ（ｋ＋
１）計算の代わりに、次式でＮ時刻分ｙ（ｋ＋Ｎ），
…，ｙ（ｋ＋２Ｎ−１）を一度に計算する。

【００２０】 Ｙ _N （ｋ＋２Ｎ−１）＝Ｘ_L,N （ｋ＋２Ｎ−１）^T Ｈ（ｋ＋Ｎ） (21) 式（１８）の右辺第２項の一部または全要素は式（２
１）でｋ＋２Ｎをｋ＋Ｎとした式で計算される。式（１
８），（１９），（２０），（２１）が従来の出力ｙ
（ｋ）が毎時刻フィルタ係数を修正する射影法と異なる
ブロック射影法の計算手順であり、図７Ｂにその機能ブ
ロック図を表わす。図７Ｂにおいて、プレフィルタ係数
計算部２１、フィルタ係数更新部２２、ブロックフィル
タリング部４１、ブロック誤差計算部４２は、それぞれ
式（１９），（２０），（２１），（１８）の演算を実
行する。なお、図７Ｂに対応する従来のブロック射影法
の計算手順を図９に示す。

【００２１】ブロック処理による演算量が低減するには
２つの側面がある。１つは、フィルタ係数Ｈ（ｋ）の修
正が射影法では毎時刻行われるのに対して、ブロック射
影法ではＮ（＞１）時刻に１回に減るというものであ
る。もう１つは、式（２０）の右辺のＸ_L,P （ｋ＋Ｎ−
１）Ｇ（ｋ＋Ｎ−１）と式（２１）右辺のＸ_L,N （ｋ＋
２Ｎ−１）^T Ｈ（ｋ＋Ｎ）が高速に計算可能なことによ
る演算量の低減である。後者の内のＸ_L,N （ｋ＋２Ｎ−
１）^T Ｈ（ｋ＋Ｎ）の計算について説明する。

【００２２】Ｘ_L,N （ｋ＋２Ｎ−１）^T Ｈ（ｋ＋Ｎ）を
効率的に計算するには、次式のように、Ｘ_L,N （ｋ＋２
Ｎ−１）をＮ行Ｎ列の部分行列、Ｈ（ｋ＋Ｎ）をＮ要素
の部分ベクトルに分解する。ただし、Ｌ＝ＮＭ（Ｍは整
数）となるように、Ｈは値が０の要素をつけ加える。 Ｘ_L,N （ｋ＋２Ｎ−１）^T ＝［Ｘ_N,N （ｋ＋２Ｎ−１）^T ，Ｘ_N,N （ｋ＋Ｎ−１），…， Ｘ_N,N （ｋ＋２Ｎ−１−（Ｍ−１）Ｎ）^T ］^T (22) Ｈ（ｋ＋Ｎ）^T ＝［Ｈ（ｋ＋Ｎ）₀ ^T ，Ｈ（ｋ＋Ｎ）₁ ^T ，…， Ｈ（ｋ＋Ｎ）_M-1 ^T ］^T このように分解すると、Ｘ_L,N （ｋ＋２Ｎ−１）^T Ｈ
（ｋ＋Ｎ）は、Ｘの部分行列、Ｈの部分ベクトル間の積
の和として計算できる。

【００２３】 Ｘ_L,N （ｋ＋２Ｎ−１）^T Ｈ（ｋ＋Ｎ）＝ Σ_i=0 ^M-1 Ｘ_N,N （ｋ＋２Ｎ−１−ｉＮ）^T Ｈ（ｋ＋Ｎ）_i (23) Σの中の各々の正方行列Ｘ_N,N （ｋ＋２Ｎ−１−ｉＮ）
とベクトルＨ（ｋ＋Ｎ）_i の積は、入力信号のベクトル
Ｘ_2N（ｋ＋２Ｎ−１−ｉＮ）とフィルタ係数Ｈ（ｋ＋
Ｎ）_i の畳み込みとみなすことができる。畳み込みを高
速に行う方法として、ＦＦＴを用いる方法、Ｗｉｎｏｇ
ｒａｄの高速畳み込み法に基づく方法（詳しくは文献Ｂ
ｌａｈｕｔ，ＦＡＳＴ ＡＬＧＯＲＩＴＨＭＳ ＦＯＲ
ＤＩＧＩＴＡＬ ＳＩＧＮＡＬ ＰＲＯＣＥＳＳＩＮ
Ｇ，ＡＤＤＩＳＯＮ−ＷＥＳＴＬＥＹ，１９８７，第
３，９章）がある。これらの、高速畳み込み法を式（２
３）に適用することで、Ｙ_N （ｋ＋Ｎ）の演算量が低減
可能である。

【００２４】ＦＦＴを用いる方法の高速の畳み込みの原
理を簡単に説明する。よく知られているように、時間領
域での畳み込みは周波数領域での積になる。 Ｆ［ｘ＊ｈ］＝Ｆ［ｘ］★Ｆ［ｈ］ (24) ここで、Ｆ［ ］はフーリエ変換を、＊は畳み込みを、
★は要素ごとの積を表す。この関係を利用すると、ｘと
ｈの畳み込みは次式で計算することができる。

【００２５】 ｘ＊ｈ＝Ｆ^-1［Ｆ［ｘ＊ｈ］］ ＝Ｆ^-1［Ｆ［ｘ］★Ｆ［ｈ］］ (25) ここで、Ｆ^-1［ ］はフーリエ逆変換を表す。フーリエ
変換に高速フーリエ変換（ＦＦＴ）を用いることで、高
速に畳み込みを行うことが可能である。式（２０）、
（２１）に適用すると次式になる。

【００２６】 Ｈ（ｋ＋Ｎ）_m ＝Ｈ（ｋ）_m ＋μＦ^-1［Ｆ［Ｘ_2p（ｋ＋Ｎ−１−ｍＮ）］ ★Ｆ［Ｇ（ｋ＋Ｎ−１）］］ ｍ＝０，１，…，Ｌ／ｐ−１ (26) Ｙ_N （ｋ＋２Ｎ−１） ＝Σ_m=0 ^M-1 Ｆ^-1［Ｆ［Ｘ_2N（ｋ＋２Ｎ−１−ｍＮ）］ ★Ｆ［Ｈ（ｋ＋Ｎ）_m ］］ ＝Ｆ^-1［Σ_m=0 ^M-1 Ｆ［Ｘ_2N（ｋ＋２Ｎ−１−ｍＮ）］ ★Ｆ［Ｈ（ｋ＋Ｎ）_m ］］ (27) 式（２１）に対応する式（２７）を例にとると、データ
数２ＮのＦＦＴが、Ｘ _2N（ｋ＋２Ｎ−１）に１回（Ｘ_2N
（ｋ＋２Ｎ−１−ｍＮ），ｍ＞１のＦＦＴは時刻ｋ＋２
Ｎ−１−ｍＮで行われるので回数に入れない）、Ｈ（ｋ
＋Ｎ）_m にＭ回、データ数２Ｎの逆ＦＦＴが１回、そし
て、複素数の積和演算が２ＮＭ回必要である。この演算
量を実数の加算と積和演算の合計数に置き換えると、デ
ータ長が２ＮのＦＦＴ、逆ＦＦＴは約８Ｎｌｏｇ₂ （２
Ｎ）、１回の複素数の積和演算は４となる。これらよ
り、式（２７）に必要な演算量は約（Ｍ＋２）８Ｎｌｏ
ｇ₂Ｎ＋８ＮＭとなり、ほぼ、ＭＮｌｏｇ₂ Ｎに比例す
る。一方、式（２７）を式（２１）のように直接計算し
た場合の演算量は、ＮＬ（＝Ｎ² Ｍ）となるので、ＦＦ
Ｔを用いることで演算量が削減されることがわかる。

【００２７】＜Ｂ＞ 毎時刻フィルタ係数を修正する射
影法と同じ出力を持つブロック処理射影法 先に３．２．４＜Ａ＞において説明した、ブロック射影
法では、ブロック長Ｎに一度しか、適応フィルタを修正
しないので、ブロック長Ｎを長くするにつれて、適応フ
ィルタの収束が遅くなるという欠点がある。

【００２８】この欠点に対して、高速射影法における、
相関ベクトルＲ_p-1 （ｋ）と平滑化係数ベクトルＳ
_p （ｋ）の計算範囲を拡張することで解決を図り、適応
フィルタの収束速度の低下を伴うことなく、射影法の演
算量のブロック処理による低減を可能にするのが、ここ
で説明する毎時刻フィルタ係数を修正する射影法と同じ
出力を持つブロック処理である。

【００２９】以下この方法の説明を行う。先に述べた毎
時刻適応フィルタを修正する高速射影法において、式
（１２）からＺ（ｋ＋ｉ＋１），（ｉ＞＝０）はＺ
（ｋ）を用いて表すことができる。 Ｚ（ｋ＋１＋ｉ）＝Ｚ（ｋ＋ｉ）＋ｓ_p （ｋ＋ｉ）Ｘ_L （ｋ＋１＋ｉ−ｐ） ＝Ｚ（ｋ＋ｉ−１）＋ｓ_p （ｋ＋ｉ−１）Ｘ_L （ｋ＋ｉ−ｐ） ＋ｓ_p （ｋ＋ｉ）Ｘ_L （ｋ＋１＋ｉ−ｐ） … ＝Ｚ（ｋ）＋Σ_j=0 ⁱ ｓ_p （ｋ＋ｉ−ｊ）Ｘ_L （ｋ＋１＋ｉ−ｐ−ｊ） ＝Ｚ（ｋ）＋Ｘ_L,i+1 （ｋ＋１＋ｉ−ｐ）Ｓ_p+i （ｋ＋ｉ）_p 〜_p+i (28) であり、特に、ｉ＝Ｎ−１のときに近似的フィルタ係数
Ｚ（ｋ）からＺ（ｋ＋Ｎ）への更新式になる。

【００３０】 Ｚ（ｋ＋Ｎ）＝ Ｚ（ｋ）＋Ｘ_L,N （ｋ＋Ｎ−ｐ）Ｓ_p+N-1 （ｋ＋Ｎ−１）_p 〜_p+N-1 (29) 式（２８）における、Ｓ_p+i （ｋ＋ｉ）は式（１５）の
次数ｐのＳ_p （ｋ）に対し次数ｐ＋ｉとされたものであ
り、これを拡張平滑化係数ベクトルと呼び、Ｓ _p+i （ｋ
＋ｉ）_p 〜_p+i は拡張平滑化ベクトルＳ_p+i （ｋ＋ｉ）
のｐ番目の要素からｐ＋ｉ番目の要素からなるベクトル
を意味する。Ｓ_p+i （ｋ＋ｉ）_p 〜_p+iは次式で再帰的
に計算される。

【００３１】 また、拡張平滑化係数ベクトルＳ_p+i （ｋ＋ｉ）は、式
（１６）のＳ_p-1 （ｋ＋ｉ）とＳ_p+i （ｋ＋ｉ）_p 〜
_p+i を結合することで次式で計算される。

【００３２】 式（３０）と平滑化係数ベクトルの回帰式（１１）を順
に式（３１）に代入することで、拡張平滑ベクトルＳ
_p+i （ｋ＋ｉ）の回帰式が得られる。

【００３３】 ところで、式（１４）において、ｋをｋ＋ｉ−１とする
ことで、次式が得られる。

【００３４】 ｙ（ｋ＋ｉ）＝Ｘ_L （ｋ＋ｉ）^T Ｚ（ｋ＋ｉ） ＋Ｒ_p-1 （ｋ＋ｉ）^T Ｓ_p-1 （ｋ＋ｉ−１） (33) また、式（２８）の右辺から、Ｚ（ｋ＋ｉ）をＺ（ｋ）
で表わす式が得られる。 Ｚ（ｋ＋ｉ）＝ Ｚ（ｋ）＋Ｘ_L,i （ｋ＋ｉ−ｐ）Ｓ_p+i-1 （ｋ＋ｉ−１）_p 〜_p+i-1 (34) 式（３４）を式（３３）に代入すると、次式となる。 ｙ（ｋ＋ｉ）＝ Ｘ_L （ｋ＋ｉ）^T ［Ｚ（ｋ） ＋Ｘ_L,i （ｋ＋ｉ−ｐ）Ｓ_p+i-1 （ｋ＋ｉ−１）_p 〜_p+i-1 ］ ＋Ｒ_p-1 （ｋ＋ｉ）^T Ｓ_p-1 （ｋ＋ｉ−１） ＝Ｘ_L （ｋ＋ｉ）^T Ｚ（ｋ） ＋Ｘ_L （ｋ＋ｉ）^T Ｘ_L,i （ｋ＋ｉ−ｐ） Ｓ_p+i-1 （ｋ＋ｉ−１）_p 〜_p+i-1 ＋Ｒ_p-1 （ｋ＋ｉ）^T Ｓ_p-1 （ｋ＋ｉ−１） ＝Ｘ_L （ｋ＋ｉ）^T Ｚ（ｋ） ＋Ｒ_p+i-1 （ｋ＋ｉ）^T _p 〜_p+i-1 Ｓ_p+i-1 （ｋ＋ｉ−１）_p 〜_p+i-1 ＋Ｒ_p-1 （ｋ＋ｉ）^T Ｓ_p-1 （ｋ＋ｉ−１） (35) ここで、Ｒ_p+i-1 （ｋ＋ｉ）を拡張相関ベクトルと呼
び、Ｒ_p+i-1 （ｋ＋ｉ） _p 〜_p+i-1 は拡張相関ベクトル
のｐ番目の要素からｐ＋ｉ−１番目の要素からなるベク
トルを意味する。Ｒ_p+i-1 （ｋ＋ｉ）_p 〜_p+i-1 はＲ
_N+p-2 （ｋ＋ｉ）_p〜_N+p-2 の一部またはすべての要素
であり、Ｒ_N+p-2 （ｋ＋ｉ）_p 〜_N+p-2 は次式で回帰的
に計算される。

【００３５】 Ｒ_N+p-2 （ｋ＋ｉ）^T _p 〜_N+p-2 ＝Ｘ_L （ｋ＋ｉ）^T Ｘ_L,N-1 （ｋ＋ｉ−ｐ） ＝Ｒ_N+p-2 （ｋ＋ｉ−１）^T _p 〜_N+p-2 ＋ｘ（ｋ＋ｉ）Ｘ_N-1 （ｋ＋ｉ−ｐ） −ｘ（ｋ＋ｉ−Ｌ）Ｘ_N-1 （ｋ＋ｉ−ｐ−Ｌ） (36) また、拡張相関ベクトルは式（１７）のＲ_p-1 （ｋ＋
ｉ）とＲ_p+i-1 （ｋ＋ｉ）_p 〜_p+i-1 を結合することで
計算される。

【００３６】 拡張相関ベクトルの時間更新は、式（１３）と式（３
６）を式（３７）に代入して次のようにかける。

【００３７】 Ｒ_N+p-2 （ｋ＋ｉ）＝ Ｒ_N+p-2 （ｋ＋ｉ−１）＋ｘ（ｋ＋ｉ）Ｘ_N+p-2 （ｋ＋ｉ−１） −ｘ（ｋ＋ｉ−Ｌ）Ｘ_N+p-2 （ｋ＋ｉ−Ｌ−１） (38) 拡張平滑化係数ベクトル、拡張相関ベクトルを用いる
と、式（３５）は次のようにかける。 ｙ（ｋ＋ｉ）＝ｙ′（ｋ＋ｉ）＋Ｒ_i+p-1 （ｋ＋ｉ）^T Ｓ_i+p-1 （ｋ＋ｉ−１） (39) ｙ′（ｋ＋ｉ）＝Ｘ_L （ｋ＋ｉ）^T Ｚ（ｋ） (40) ｙ（ｋ），Ｅ_p （ｋ），Ｇ（ｋ），Ｓ_p （ｋ）は高速射
影法の手順で計算できるので、ｉ＝０として式（３９）
を用いてｙ（ｋ）が計算できる。式（７．２）を用いて
ｙ（ｋ）からｅ（ｋ），Ｅ_p （ｋ）が計算でき、Ｅ
_p （ｋ）からＧ（ｋ）が式（８）により計算できる。Ｇ
（ｋ）から、式（１１），（３０），（３１）を用い
て、Ｓ_p+1 （ｋ）を計算でき、式（３９）により、ｙ
（ｋ＋１）が計算できる。このように、順にｙ（ｋ＋
ｉ），（ｉ＝０，１，…，Ｎ−１）は計算することがで
き、時刻ｋ＋Ｎ−１において式（２９）からＺ（ｋ＋
Ｎ）が計算できる。

【００３８】まとめると、図１１に示すように、まず初
期化として、Ｒ_p+N-1 （０）＝Ｘ_L （０）^T ［Ｘ_L （−
１），Ｘ_L （−２），…，Ｘ_L （−ｐ−Ｎ＋１）］^T と
し、Ｅ_p-1 （０）＝Ｓ_p-1 （０）＝Ｘ_L （０）＝［０，
０，…，０］^T とし（Ｓ₁₀₀ ）、次に、ｋ＝１とし（Ｓ
₁₀₁ ）、次に入力ベクトルと近似的ベクトルとの内積演
算 ｙ′（ｋ＋ｉ）＝Ｘ_L （ｋ＋ｉ）^T Ｚ（ｋ），ｉ＝０，
１，…，Ｎ−１ を行う（Ｓ₁₀₂ ）。その後、ｉ＝０，１，…，Ｎ−１に
対して、式（３８），（３９），（７．２），（８），
（３２）の演算を各時刻ごとに行い、ｉ＝Ｎ−１の演算
が終了すると（Ｓ₁₀₃ ）、式（２９）を演算して、Ｚ
（ｋ）からＺ（ｋ＋Ｎ）に更新する（Ｓ₁₀₄ ）。次に、
ｋをｋ＋ＮとしてステップＳ₁₀₂ に戻る（Ｓ₁₀₅ ）。

【００３９】この手順に対応する機能ブロック図は図１
０に示し、図７Ａ、図７Ｂと対応する部分と同一符号を
付けて示す。入力信号ベクトルから相関計算部３３で式
（１３）が演算され、また入力信号ベクトルから相関計
算拡張部５２で式（３６）が演算され、これら相関計算
部３３からの相関ベクトルと、相関計算拡張部５２から
の相関ベクトルの拡張部のベクトルとが相関ベクトル連
結部３９で式（３７）により連結されて拡張相関ベクト
ルが得られる。要は、これら相関計算部３３、相関計算
拡張部５２、相関ベクトル連結部５３は式（３８）を演
算するものであり、拡張相関ベクトル演算部５６を構成
する。

【００４０】プレフィルタ係数計算部２１で入力信号ベ
クトルと、誤差信号ベクトルとにより式（８）を演算し
てプレフィルタ係数を演算し、このプレフィルタ係数は
プレフィルタ係数平滑部３１で式（３０）の演算により
求められ、この拡張部分のベクトルと平滑部３１よりの
平滑係数ベクトルとが平滑係数ベクトル連結部５５で式
（３１）の演算により連結されて拡張平滑係数ベクトル
が得られる。要は、プレフィルタ係数平滑部３１、平滑
係数拡張部５４、平滑係数ベクトル連結部５５は式（３
２）を演算する拡張平滑係数ベクトル計算部５７を構成
している。

【００４１】拡張相関ベクトル演算部５６よりの拡張相
関ベクトルと拡張平滑係数ベクトル計算部５７よりの拡
張平滑係数ベクトルとの内積演算、つまり、式（３９）
の右辺第２項の演算がなされる。平滑係数拡張部５４よ
りの拡張部分ベクトルが近似的フィルタ係数ブロック更
新部５１へ入力されて式（３４）が演算される。この演
算結果の近似的フィルタ係数と入力ベクトルとの内積演
算がブロックフィルタリング部４１で行われ、その演算
結果と内積演算計算部３５からの演算結果とが出力合成
部３４で加算されて出力信号とされる。ブロックフィル
タリング部４１、近似的フィルタ係数ブロック更新部５
１はＮ時刻ごとに演算され、その他の各部は各時刻ごと
に行われ、拡張平滑係数ベクトル部５７では、ブロック
フィルタリング部４１で演算され、次の演算まで各時刻
ごとに字数がｐから１次づつ拡張され、このことが繰り
返され、拡張相関ベクトル計算部５６では、常時、拡張
平滑係数ベクトルの最大次数ｐ＋Ｎ−１と同一次数の相
関ベクトルが演算され内積演算計算部３５に対しては、
拡張平滑係数ベクトル計算部５７よりのその時刻での拡
張平滑係数ベクトルの字数と同一次数分が入力される。

【００４２】以上説明した、毎時刻フィルタ係数を修正
する射影法と同じ出力を持つブロック射影法の演算量を
求める。まず、同方法の手順を示した図１１のＳ
₁₀₂ は、３．２．４＜Ａ＞で説明したように、ＦＦＴを
用いた高速畳み込みが利用できその演算量は、約８（Ｍ
＋２）Ｎｌｏｇ₂ ２Ｎ＋８ＮＭであり、１時刻当たりで
は、８（Ｍ＋２）ｌｏｇ₂ ２Ｎ＋８Ｍである。Ｓ₁₀₃ で
は、１時刻当たり、式（３８）に２（Ｎ＋ｐ−２）回、
式（３９）にｐ＋Ｎ／２回、式（７．２）にｐ回、式
（８）には１５ｐ回、式（３２）には、ｐ回の演算が必
要であり、合計すると２．５Ｎ＋２０ｐとなる。Ｓ₁₀₄
は３．２．４＜Ａ＞で説明した要領でＦＦＴを用いて高
速に計算できる。つまり、式（２９）は式（２０）と同
様に次のようにかけるので、 Ｚ（ｋ＋Ｎ）_m ＝ Ｚ（ｋ）_m ＋μＦ^-1［Ｆ［Ｘ_2N（ｋ＋Ｎ−ｐ−ｍＮ）］ ★Ｆ［Ｓ（ｋ＋Ｎ−１）_p,N+p-1 ］］ ｍ＝０，１，…，Ｌ／Ｎ−１ (41) データ長が２ＮのＸ_2N（ｋ＋Ｎ−ｐ）をＦＦＴするのに
１回、Ｓ（ｋ＋Ｎ−１）にＮ個０を付加したベクトルを
長さ２ＮのデータのＦＦＴが１回、そして、逆ＦＦＴが
Ｌ／Ｎ（＝Ｍ）回必要である。また、複素数の積和演算
が２ＭＮ回必要である。これら、ＦＦＴ、逆ＦＦＴ、複
素数の積和演算に要する演算量は約８（Ｍ＋２）Ｎｌｏ
ｇ₂ ２Ｎ＋８ＮＭであり、１時刻当たりでは、８（Ｍ＋
２）ｌｏｇ ₂ ２Ｎ＋８Ｍである。以上より、Ｓ₁₀₂ 、Ｓ
₁₀₃ そしてＳ₁₀₄ に要する１時刻当たりの演算量をＴと
かくと、 Ｔ＝１６（Ｍ＋２）ｌｏｇ₂ （２Ｎ）＋２．５Ｎ＋１６Ｍ＋２０ｐ (42) と、見積られる。

【００４３】

【発明が解決しようとする課題】従来の毎時刻フィルタ
係数を修正する射影法（以下、従来法）と同じ出力ｙ
（ｋ）を持つブロック射影法では、演算量Ｔを最小にす
るブロック長が存在し、それ以上に長いブロック長で
は、演算量の低減効果が小さくなるか、または無くなる
という問題点がある。この理由は、式（３８）で計算さ
れる拡張相関ベクトルＲ_N+p-2 の時間更新に要する演算
量２Ｎと式（３９）右辺第２項で計算される拡張相関ベ
クトルと拡張平滑係数ベクトルの積に要する演算量０．
５Ｎがブロック長Ｎに比例しているので、ブロック長Ｎ
を長くするとともに増加することにある。式を用いれ
ば、式（４２）のＴのブロック長Ｎに関して微係数が正
となることで分かる。

【００４４】 ∂Ｔ／∂Ｎ＝１６（（−Ｌ／Ｎ² ）ｌｏｇ₂ （２Ｎ） ＋（（Ｌ／Ｎ）＋２）／（Ｎｌｎ２）） ＋２．５−１６（Ｌ／Ｎ² ） ＝Ｎ^-2（１６（−Ｌｌｏｇ₂ （２Ｎ） ＋（Ｌ＋２Ｎ）／ｌｎ２）＋２．５Ｎ² −１６Ｌ） (43) Ｎ^-2は常に正であり、１６（−Ｌｌｏｇ₂ （２Ｎ）〜＋
（Ｌ＋２Ｎ）／ｌｎ２）＋２．５Ｎ² −１６ＬはＬ＞＞
Ｎでは負の値をとるが、Ｎが大きくなるとＮ、Ｎ ² の項
に比べてｌｏｇ₂ （２Ｎ）の項が無視できて、下に凸な
Ｎの２次式とみなせるので、値は正になる。つまり、演
算量はＮ＜＜Ｌでは減少するが、あるＮ以上では、増加
する。

【００４５】フィルタ長Ｌが与えられたときに、演算量
Ｔが最小となるＮは、式（４３）の右辺が０となるブロ
ック長Ｎである。式（４３）の右辺を０とおいた方程式
は、以下に示すようにＮに関する再帰計算により求める
ことができる。まず、式（４３）の右辺を０とおいた方
程式をＮの２次式とみなしてとく。 Ｎ＝［16/ln2+4√(1/ln2)²+2.5L(log₂(2N)-1/ln2+1) ］/2.5 ＝9+1.6 √2.5Llog₂(2N) ＝9+2.5 √Llog₂(2N) (44) 式（４４）の左辺のブロック長Ｎを右辺に代入仕直すこ
とで、再帰計算を行う。ｉ回目の再帰計算によるブロッ
ク長ＮをＮ（ｉ）と表わすことにすると式（４４）は次
式になる。

【００４６】 Ｎ（ｉ＋１）＝９＋２．５√Ｌｌｏｇ₂ （２Ｎ（ｉ）） (45)

【００４７】

【課題を解決するための手段】この発明では、上記の課
題に対して、従来法における、式（４０）で計算される
出力計算、式（２９）で計算される近似的フィルタ係数
Ｚの更新において異なるブロック長を用いることで解決
を図り、従来法では得られない演算量の低減や遅延量の
短縮を可能にするものである。

【００４８】この発明によれば、入力信号を可変フィル
タ処理し、その出力信号と所望信号との誤差信号を求
め、その誤差信号のパワーを最小とするように、上記可
変フィルタの特性を適応的に制御する方法において、上
記入力信号ベクトルの自己相関行列の逆行列と、上記誤
差信号ベクトルとの積に相当する演算を行ってプレフィ
ルタ係数ベクトルを求める第１ステップと、上記プレフ
ィルタ係数ベクトルを平滑化して平滑係数ベクトルを得
る第２ステップと、上記平滑係数ベクトルの要素を加重
係数として、上記入力信号ベクトルの方向に近似的フィ
ルタ係数ベクトルを修正する第３ステップと、上記入力
信号ベクトルと上記近似的フィルタ係数ベクトルとの内
積演算を行う第４ステップと、上記入力信号ベクトルの
自己相関ベクトルを求める第５ステップと、上記平滑係
数ベクトルと上記自己相関ベクトルとの内積演算を行う
第６ステップと、上記第４ステップで得られた内積演算
結果と、上記第６ステップで得られた内積演算結果とを
加算して上記出力信号を得る第７ステップとを有するも
のであって、上記第３ステップの修正及び上記第４ステ
ップの演算はそれぞれＮ₁ 、Ｎ₂ 時刻（Ｎ₁ 、Ｎ₂ は２
以上の整数、Ｎ₁ ≠Ｎ₂ ）ごとに行い、上記第３、第４
ステップ以外の各ステップは各時刻ごとに演算を行い、
上記第３ステップの演算が行われるごとに、次に上記第
３ステップの演算されるまでの各時刻において、上記第
２ステップで得る上記平滑係数ベクトルの次数を１ずつ
増加し、次の第３ステップの演算時刻での演算でもとの
次式に戻し、上記第５ステップでは上記平滑係数ベクト
ルの最大次数と等しい次数の上記自己相関ベクトルを演
算し、上記第６ステップでは上記平滑係数ベクトルの次
数と等しい次数だけ上記第５ステップで得られた自己相
関ベクトルから用いることを特徴とする。

【００４９】

【発明の実施の形態】以下、この発明の実施の形態を説
明する。一般にブロック処理を行う適応アルゴリズムに
おいては、ブロック処理により演算量を低減可能な計算
ステップが２箇所ある。先に述べた毎時刻適応フィルタ
を修正する射影法と同じ出力を持つブロック射影法にお
いては、図１１中のＳ₁₀₂ で式（４０）により計算され
る出力計算と図１１中Ｓ₁₀₄ で式（２９）で計算される
近似的フィルタ係数Ｚの更新である。ところで、先に述
べた従来法において課題となった、長いブロック長Ｎに
対して演算低減効果が小さくなる理由は式（３８）で計
算される拡張相関ベクトルＲ_N+p-2 の時間更新に要する
演算量２Ｎと式（３９）右辺第２項で計算される拡張相
関ベクトルと拡張平滑係数ベクトルの積に要する演算量
０．５Ｎがブロック長Ｎに比例しているので、ブロック
長Ｎを長くするとともに増加することであった。これら
ブロック長に比例して増加する演算は近似的フィルタの
更新をブロック長の間だけ保留しているために生じてい
る。そこで、上記の２箇所のブロック処理において異な
るブロック長を用いることで上記の課題を解決すること
ができる。

【００５０】二つのブロック処理におけるブロック長を
独立したことにより、この発明法の処理は従来法と処理
順序または処理のタイミングにおいて異なったものとな
る。出力の計算でのブロック長をＮ₁ 、近似的フィルタ
係数更新でのブロック長をＮ₂ と書くと、式（４０），
（３８），（３９），（７．２），（８），（３２）そ
して（２９）で与えられる従来法での処理に対応して、
本発明法の処理は以下のようになる。

【００５１】 ｙ′（ｋ₁ ＋ｉ）＝Ｘ_L （ｋ₂ ＋ｉ）^T Ｚ（ｋ₂ ）， ｉ＝０，…，Ｎ₁ −１ (46) Ｒ_N2+p-2（ｋ₁ ）＝Ｒ_N2+p-2（ｋ−１） ＋ｘ（ｋ₁ ）Ｘ_N2+p-2（ｋ₁ −１） −ｘ（ｋ₁ −Ｌ）Ｘ_N2+p-2（ｋ₁ −Ｌ−１） (47) ｙ（ｋ₁ ）＝ｙ′（ｋ₁ ）＋Ｒ_n2+p-1（ｋ₁ ）^T Ｓ_n2+p-1（ｋ₁ −１） (48) ｅ（ｋ₁ ）＝ｄ（ｋ₁ ）−ｙ（ｋ₁ ） (49) Ｚ（ｋ₂ ＋Ｎ₂）＝Ｚ（ｋ₂ ） ＋Ｘ_L,N2（ｋ₂ ＋Ｎ₂ −ｐ）Ｓ_N2+p-1（ｋ＋Ｎ₂ −１）_p-N2+p-1 (53) ここで、ｋ₂ は近似的フィルタ係数Ｚの修正に関する時
刻の指標であり、ｋ₁ はその他の変数の時刻の指標であ
る。また、ｎ₂ ＝ｋ₁ ｍｏｄ Ｎ₂ −１なる値をもつ変
数である。図１はこの発明法の計算手順を示している。
まず、初期化として、 Ｒ_p+N2-1（０）＝Ｘ_L (0) ^T ［Ｘ_L (-1), Ｘ_L (-2), …, Ｘ_L (-p-N₂+1) ］^T Ｓ_p-1 （０）＝［０，０，…，０］^T Ｚ_L-1 （１）＝［０，０，…，０］^T Ｅ_p-1 （０）＝［ｙ（０），ｙ（−１），…，ｙ（−ｐ＋２）］^T とし（Ｓ₀ ）、次にｋ１＝ｋ２＝１，ｎ ₂ ＝０とし（Ｓ
₁ ）、次に、もし、ｋ₁ｍｏｄ Ｎ₁ ＝１であるかの判
定をし（Ｃ₁ ）、判定Ｃ₁ が真であれば、入力ベクトル
と近似的フィルタ係数との内積演算 ｙ′（ｋ₁ ＋ｉ）＝Ｘ_L （ｋ₂ ＋ｉ）^T Ｚ（ｋ₂ ）， ｉ＝０，…，Ｎ₁ −１ （Ｓ₂ ）を行う。その後、式（４７），（４８），（４
９），（５０），（５１），（５２）と順次計算をおこ
ない（Ｓ₃ ）、次に、ｋ₁ とｎ₂ を１増やし（Ｓ₄ ）、
次に、もし、ｋ₁ ｍｏｄ Ｎ₂ ＝１であるかの判定をし
（Ｃ₂ ）、判定Ｃ₂ が真であれば、近似的フィルタ係数
Ｚの修正（Ｓ₅ ） Ｚ（ｋ₂ ＋Ｎ₂ ）＝Ｚ（ｋ₂ ）＋Ｘ_L,N2（ｋ₂ ＋Ｎ₂ −
ｐ）Ｓ_N2+p-1（ｋ＋Ｎ₂ −１）_p-N2+p-1 を行い、ｋ₂ をＮ₂ だけ増やし、ｎ₂ を０にする
（Ｓ₆ ）。次に判定Ｃ１に戻る。

【００５２】この手順に対応する機能ブロック図は図２
に示し、図７Ａ、図７Ｂと対応する部分と同一符号を付
けて示す。図２のユニットは、従来法の機能ブロック図
である図１１中のユニットと機能は同じだが動作の順序
が異なっている。この発明の演算量を見積る。まず、同
方法の手順を示した図１のＳ₂ は、従来法で説明したよ
うに、ＦＦＴを用いた高速畳み込みが利用でき、その演
算量は、約Ｔ₀ ＝８（Ｍ₁ ＋２）Ｎ₁ ｌｏｇ₂ ２Ｎ₁ ＋
８Ｎ₁ Ｍ₁ であり、１時刻当たりでは、Ｔ₀ ／Ｎ₁ ＝８
（Ｍ₁ ＋２）ｌｏｇ₂ ２Ｎ₁ ＋８Ｍ₁ である。ここで、
Ｌ＝Ｎ₁ Ｍ₁ である。Ｓ₃ では、１時刻当たり、式（４
７）に２（Ｎ₂ ＋ｐ−２）回、式（４８）にｐ＋Ｎ₂ ／
２回、式（５０）にｐ回、式（５１）には１５ｐ回、式
（５２）には、ｐ回の演算が必要であり、合計すると
２．５Ｎ₂ ＋２０ｐとなる（Ｔ_c ＝Ｎ₂ （２．５Ｎ₂ ＋
２０ｐ）と定義しておく）。Ｓ₅ に要する演算量は従来
法で説明したように、Ｌ＝Ｎ₂ Ｍ₂ とすれば、約Ｔ_u ＝
８（Ｍ₂ ＋２）Ｎ₂ ｌｏｇ₂ ２Ｎ₂ ＋８Ｎ₂ Ｍ₂ であ
り、１時刻当たりでは、Ｔ_u ／Ｎ₂ ＝８（Ｍ₂ ＋２）ｌ
ｏｇ₂ ２Ｎ₂ ＋８Ｍ₂ である。以上より、この発明法の
実行に要する演算量をＴとかくと、ＴはＳ₂ 、Ｓ₃ 、そ
してＳ₅ に要する１時刻当たりの演算量の和として、 Ｔ＝Ｔ₀ ／Ｎ₁ ＋Ｔ_C ／Ｎ₂ ＋Ｔ_u ／Ｎ₂ ＝８（Ｍ₁ ＋２）ｌｏｇ₂ （２Ｎ₁ ） ＋８（Ｍ₂ ＋２）ｌｏｇ₂ （２Ｎ₂ ） ＋２．５Ｎ₂ ＋８Ｍ₁ ＋８Ｍ₂ ＋２０ｐ (54) と、見積られる。図３に式（４２）によって計算した演
算量Ｔとブロック長Ｎ₁，Ｎ₂ との関係の例を示した。
条件は、フィルタ長はＬ＝４０９６、ブロック長Ｎ₁ ，
Ｎ₂ は１２８，２５６，５１２，１０２４，２０４８で
ある。同図で、ΔはＮ₁ ＝Ｎ₂ なる従来法により可能な
演算量であり、ブロック長が５１２で演算量が最低で、
ブロック長がそれより短くても長くても演算量が増加す
る。一方、□はＮ₁ ≠Ｎ₂ なるこの発明法による演算量
であり、ブロック長Ｎ₁ を長くするとともに、演算量は
低減し、ブロック長Ｎ₁ ＝１０２４，２０４８では、等
ブロック長の従来法の場合の最小演算量より少なくする
ことが可能である。Ｎ₂ を変化させた場合、Ｎ₂ ＝２５
６で何れのＮ₁ に対しても最小となる。図３より、この
発明が従来法に比べて演算量が少なくなるのは、ブロッ
ク長Ｎ₁ が、式（４５）で計算される従来法の最小演算
量を与えるブロック長Ｎより長い場合である。Ｌ＝４０
９６のこの例では、式（４５）の再帰計算をＮ（０）＝
５１２なる初期値で開始し、Ｎ＝４８７を得る。

【００５３】また、通常は演算量が規定され、その範囲
で可能なだけ長いフィルタ長を実装するので、フィルタ
長を演算量の関数として表わす方が便利なこともある。
式（４２）にＬ＝Ｎ₁ Ｍ₁ 、Ｌ＝Ｎ₂ Ｍ₂ を代入しフィ
ルタ長Ｌについて解くと次式を得る。 Ｌ＝［T-16log₂(2N₁)-16log₂(2N₂)-2.5N₂-20p ] ／［8(log₂(2N₁)+1)/N₁+8(log₂(2N₂)+1)/N₂] (55) 図４に式（５５）によって計算した、フィルタ長Ｌとブ
ロック長Ｎ₁ ，Ｎ₂ との関係の例を示した。条件は演算
量はＴ＝８０００、ブロック長Ｎ₁ ，Ｎ₂ は１２８，２
５６，５１２，１０２４，２０４８，４０９６である。
同図で、ΔはＮ₁＝Ｎ₂ なる従来法により可能なフィル
タ長であり、ブロック長Ｎ₁ が１０２４でフィルタ長が
最長で、ブロック長がそれより短くても長くても演算量
が増加する。特に、ブロック長Ｎ₁ が４０９６以上では
Ｌ＞０となる解がない。一方、□はＮ₁ ≠Ｎ₂ なるこの
発明法によるフィルタ長であり、ブロック長Ｎ₁ を長く
するとともに、フィルタ長が長くなり、ブロック長Ｎ₁
＝２０４８以上では、等ブロック長の従来法の場合の最
長フィルタ長より長いフィルタ長を実装することが可能
である。

【００５４】つぎに、ブロック処理を行うことによるフ
ィルタ出力ｙ（ｋ）が遅延する時間を考える。この遅延
時間をＤと表わすと、図５より、次式で見積られること
がわかる。 Ｄ＝Ｎ₁ ＋Ｔ₀ ／Ｔ ，Ｎ₁ ＞Ｎ₂ ＝（２Ｔ₀ ＋Ｔ_c ＋Ｔ_u ）／Ｔ ，Ｎ ₁ ＜Ｎ₂ (56) 図５Ｃに式（５６）によって計算した出力遅延時間Ｄと
ブロック長Ｎ₂ との関係の例を示した。条件はフィルタ
長Ｌ＝４０９６、ブロック長Ｎ₁ ，Ｎ₂ は１２８，２５
６，５１２，１０２４，２０４８である。同図で、Δは
Ｎ₁ ＝Ｎ₂ なる従来法による出力遅延時間であり、□は
Ｎ₁ ≠Ｎ₂ なるこの発明法による出力遅延時間であり、
ブロック長Ｎ₁ は各々のＮ₂ について、下から１２８，
２５６，５１２，１０２４，２０４８である。Ｎ₂ が２
５６以上では、この発明法は従来法より短い出力遅延時
間を実現できることがわかる。

【００５５】以上説明したこの発明法の特徴は、Ｎ₁ ≠
Ｎ₂ とすることにより、従来の毎時刻フィルタを修正す
る射影法と同じ出力をもつブロック射影法より少ない演
算量で実行できる、いいかえると、同じ演算量でより長
いフィルタ長を使用でき、その結果、誤差信号のパワー
がより減少することである。また、出力遅延時間を従来
法より短くすることもでき、この場合はＮ₂ ＞Ｎ₁ であ
って、フィルタ係数の更新周期が長くなるため、途中で
エコーキャンセラの場合はダブルトークなどの異常を検
出したら、その時の更新を中止する処理も可能となる。

【００５６】実施例 第１の例は、音響エコーキャンセラである。図１２Ａは
この発明の方法を適応した音響エコーキャンセラの機能
構成を示す。入力信号ｘ（ｋ）として音声信号がフィル
タリング部１１、フィルタ係数更新部１２、スピーカ６
１へ供給され、マイクロホン６２の出力が所望信号ｄ
（ｋ）として差回路１３へ供給される。スピーカ６１か
ら音響信号通路、いわゆるエコー伝達経路６３を通じて
マイクロホン６２へ達する信号が音響エコーであり、適
応フィルタ（フィルタリング部）１１の出力信号ｙ
（ｋ）に推定エコー、推定誤差ｅ（ｋ）に残留エコーが
対応している。

【００５７】音響エコーキャンセラは、テレビ会議シス
テム等の拡声通話系において、ハウリングや、エコーを
防止する装置である。音響エコーキャンセラの場合、適
応フィルタ１１は、スピーカ６１からマイクロホン６２
の間の伝達経路６３を模擬する。適応フィルタ１１の出
力は推定されたエコーであり、マイクロホン６２で受け
たエコーから差し引くことでエコーが除去される。スピ
ーカ６１からマイクロホン６２の間の伝達経路６３を適
応フィルタ１１で模擬する場合に必要なフィルタの長さ
は数百から数千にもなり、多くの演算量が必要となる。
また、スピーカ６１からマイクロホン６２の間の伝達経
路６３が人の動き、ドアの開閉等で変化するので、その
変化に追従できる高速な適応アルゴリズムが必要とされ
ている。

【００５８】第２の適用例は騒音制御である。図１２Ｂ
に騒音制御の原理を示す。騒音制御の目的は、騒音源８
１から騒音伝達経路８２を径て観測される騒音を、スピ
ーカ８６から負の音（音圧がｙ（ｋ）と表される時に、
−ｙ（ｋ）と表される音圧を、ｙ（ｋ）に対する負の音
と呼ぶ）を出すことにより消去しようとするものであ
る。この例では、適応フィルタ（フィルタリング部）１
１への入力信号ｘ（ｋ）は騒音源８１の騒音をモニタマ
イクロホン８４で観測した騒音であり、所望信号ｄ
（ｋ）は、観測地点（マイクロホン８３）における騒音
であり、適応フィルタ１１の出力ｙ（ｋ）は推定騒音で
あり、誤差信号ｅ（ｋ）は観測地点で、スピーカ８６か
ら放射された推定騒音と騒音源から到来した騒音の音圧
が合成され、マイクロホン８３の出力に得られる。適応
フィルタ１１が、騒音源８１から観測マイクロホン８３
までの伝達経路８２を良好に模擬していれば、スピーカ
８６から放射される音圧−ｙ（ｋ）によって騒音ｄ
（ｋ）は消去される。

【００５９】以上、２つの適用例を示した。どちらの場
合にもエコー、騒音を消去するという目的を達成するた
めには、推定エコー、推定騒音が得られれば十分であ
り、エコー、騒音の伝達経路そのものを求める必要はな
い。従って、どちらの適用例でも近似的フィルタ係数を
使用するこの発明法が適用でき、従来装置と比較して演
算量を削減できる。更に演算量が多くなるが、出力ｙ
（ｋ）が速く得られるようにすることができる。また、
この発明は、図６Ａに示したような、構成において誤差
信号ｅ（ｋ）のパワーを最小とする基本機能を有してい
る。従って、解決すべき問題が図６Ａに示した誤差信号
のパワーの最小化としてモデル化できる全ての応用例に
対して、この発明を適用することが可能である。

【００６０】

【発明の効果】以上説明したように、この発明では、適
応アルゴリズムの１つである毎時刻フィルタを修正する
射影法と同じ出力をもつブロック射影法において２つの
相違なるブロック長を導入し、動作手順を変更すること
により、収束速度を低下させること無く、演算量低減を
可能とする。その結果、この発明は収束の遅い従来のブ
ロック射影法等の手法が利用されていた音響エコーキャ
ンセラ装置に対して、収束速度を低下させることなく、
演算量の低減が可能となる。即ち、同一のハードウェア
規模で収束は速いが多少演算量の必要なより高い次数の
射影法の利用や、フィルタ長を長くすることが可能とな
り、より高速なエコー低減や定常時のエコーの低減の実
現が達成される。またこの発明によれば、演算量は低減
できないが出力ｙ（ｋ）が得られる遅延時間を短かくす
ることができ、その場合はフィルタの更新時間が長くな
り、この間にダブルトークなどの異常の発生を検出処理
を行い、検出すればフィルタの次の更新を中止するなど
の処理を行うことができる。

【図面の簡単な説明】

【図１】この発明の実施例を示す流れ図。

【図２】この発明の実施例を適用した装置の機能構成を
示す図。

【図３】フィルタ長を一定とした時のブロック長Ｎ₁ ，
Ｎ₂ と演算量Ｔの関係を示す図。

【図４】演算量を一定とした時のブロック長Ｎ₁ ，Ｎ₂
とフィルタ長Ｌとの関係を示す図。

【図５】Ａ及びＢはそれぞれ演算実行時間と出力遅延時
間とを示す図、Ｃはブロック長Ｎ₁ ，Ｎ₂ と出力遅延時
間Ｄとの関係を示す図である。

【図６】Ａは適応フィルタを用いた適応制御系を示す
図、Ｂは従来の射影法を用いた制御の機能構成を示す図
である。

【図７】Ａは従来の高速射影法を用いた制御の機能構成
を示す図、Ｂは従来のブロック射影法を用いた制御の機
能構成を示す図である。

【図８】従来の高速射影法を用いた制御の計算手順を示
す流れ図。

【図９】従来のブロック法を用いた制御の計算手順を示
す流れ図。

【図１０】従来の毎時刻フィルタを修正する射影法と同
じ出力をもつブロック射影法を用いた制御の機能構成を
示す図。

【図１１】従来の毎時刻フィルタを修正する射影法と同
じ出力をもつブロック射影法を用いた制御の計算手順を
示す流れ図。

【図１２】Ａはこの発明を音響エコーキャンセラに適用
した機能構成を示す図、Ｂはこの発明を騒音制御に適用
した機能構成を示す図である。

───────────────────────────────────────────────────── フロントページの続き (56)参考文献 特開 平９−54590（ＪＰ，Ａ) 特開 昭59−139717（ＪＰ，Ａ) (58)調査した分野(Int.Cl.⁷，ＤＢ名) H03H 17/00 - 17/08 H03H 21/00 H04B 3/20 - 3/23

Claims (1)

(57)【特許請求の範囲】
1. 【請求項１】 入力信号を可変フィルタ処理し、その出
   力信号と所望信号との誤差信号を求め、その誤差信号の
   パワーを最小とするように、上記可変フィルタの特性を
   適応的に制御する方法において、 上記入力信号ベクトルの自己相関行列の逆行列と、上記
   誤差信号ベクトルとの積に相当する演算を行ってプレフ
   ィルタ係数ベクトルを求める第１ステップと、 上記プレフィルタ係数ベクトルを平滑化して平滑係数ベ
   クトルを得る第２ステップと、 上記入力信号ベクトルを列要素に持つ行列と上記平滑係
   数ベクトルの部分ベクトルの積を前回の近似的フィルタ
   係数ベクトルに加えることで近似的フィルタ係数ベクト
   ルをブロック処理により更新修正する第３ステップと、 上記入力信号ベクトルと上記近似的フィルタ係数ベクト
   ルとの内積演算をブロック処理により行う第４ステップ
   と、 上記入力信号ベクトルの自己相関ベクトルを求める第５
   ステップと、 上記平滑係数ベクトルと上記自己相関ベクトルとの内積
   演算を行う第６ステップと、 上記第４ステップで得られた内積演算結果と、上記第６
   ステップで得られた内積演算結果とを加算して上記出力
   信号を得る第７ステップとを有するものであって、 上記第３ステップの修正及び上記第４ステップのブロッ
   ク処理による演算はそれぞれＮ₁ 、Ｎ₂ 時刻（Ｎ₁ 、Ｎ
   ₂ は２以上の整数、Ｎ₁ ≠Ｎ₂ ）ごとに行い、 上記第３、第４ステップ以外の各ステップは各時刻ごと
   に演算を行い、 上記第３ステップの演算が行われるごとに、次に上記第
   ３ステップの演算されるまでの各時刻において、上記第
   ２ステップで得る上記平滑係数ベクトルの次数を１ずつ
   増加し、次の第３ステップの演算時刻での演算でもとの
   次数に戻し、 上記第５ステップでは上記平滑係数ベクトルの最大次数
   と等しい次数の上記自己相関ベクトルを演算し、 上記第６ステップでは上記平滑係数ベクトルの次数と等
   しい次数だけ上記第５ステップで得られた自己相関ベク
   トルから用いることを特徴とする適応的制御方法。