JP3452339B2

JP3452339B2 - 適応的制御方法

Info

Publication number: JP3452339B2
Application number: JP20616395A
Authority: JP
Inventors: 雅史田中; 昭二牧野; 順治小島
Original assignee: Nippon Telegraph and Telephone Corp
Current assignee: Nippon Telegraph and Telephone Corp
Priority date: 1995-08-11
Filing date: 1995-08-11
Publication date: 2003-09-29
Anticipated expiration: 2015-08-11
Also published as: JPH0954590A

Description

【発明の詳細な説明】

【０００１】

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

【０００２】

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

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

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

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

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

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

【０００８】ここで、この系の計算量を考えてみる。ま
ず、プレフィルタ係数計算部２１で行う式（８）の演算
は、ｐ行ｐ列の行列Ｘ_L,p(k) ^TＸ_L,p(k) の逆行列の
計算を含んでいるので、ｐの３乗に比例する積和演算量
（以降、Ｏ（p³) と表す）が必要となる。これは、ｐが
大きくなると、たいへん大きな演算量となる。次に、フ
ィルタ係数更新部２２で行う式（９）の演算はＸ
_L,p(k) がＬ行ｐ列の行列、Ｇ(k) がｐ次のベクトルで
あることより、ｐＬ回の積和演算が必要である。更にフ
ィルタリング部１１で行う式（10）の演算はＬ次のベク
トルの内積であるので、Ｌ回の積和演算が必要である。
従って、合計（p+1)Ｌ＋Ｏ（p³) の積和演算量が必要と
なる。例えば、ｐ＝２０、Ｌ＝１０００の場合には、２
１０００回以上もの積和演算が必要となる。

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

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

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

【００１２】Ｚ(k+1) ＝Ｚ(k) ＋Ｘ_L(k-p+1) ｓ_p(k) （12）Ｒ_p-1(k+1)＝Ｒ_p-1(k)＋ｘ(k+1) Ｘ_p-1(k)−ｘ(k-L+1) Ｘ_p-1(k-L) （13）ｙ(k+1) ＝Ｘ_L(k+1) ^TＺ(k+1) ＋Ｒ_p-1(k+1)^TＳ_p-1(k) （14）ただし、Ｓ_p(k) およびＳ_p-1(k-1)は、それぞれｐ次お
よびｐ−１次の平滑化係数ベクトルであって、Ｓ_p(k) ＝［ｓ₁(k), ｓ₂(k)，…，ｓ_p(k) ］^T （15）Ｓ_p-1(k) ＝［ｓ₁(k), ｓ₂(k)，…，ｓ_p-1(k) ］^T （16）と表される。式（11）からわかるように、Ｓ_p(k) はプ
レフィルタ係数ベクトルＧ(k) を用いて再帰的に計算さ
れる。また、Ｚ(k) は、近似的フィルタ係数ベクトルで
あり、Ｘ_L(k-p+1) とｓ_p(k) を用いて再帰的に計算さ
れる。更に、Ｒ_p- ₁(k+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）このとき、数学的な証明は省略するが、図６Ａの構成、
または式（11）−（14）の演算により、入力信号ｘ(k)
および誤差信号ｅ(k) を用いて出力信号ｙ(k+1) を合成
することができる。なお、式（13）（14）の時間パラメ
ータがｋ＋１となっているのに図６Ａではｋとなってい
る。これは時刻ｋにおいて式（13），（14），（11），
（12）の順に計算することを表したものである。

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

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

【００１６】ブロック処理を射影法に適用した従来のブ
ロック射影法では、Ｅ_p(k),Ｇ(k)の計算、そして、Ｈ
(k+N) の更新は各々式(7.1) ，（８），（９）と同様に
行われる。Ｅ_p(k+N-1) ＝［ｄ(k+N-1),ｄ(k+N-2),…，ｄ(k+N-p）］^T −Ｘ_L,p(k+N-1) ^TＨ(k) （18）Ｇ(k+N-1) ＝［Ｘ_L,p(k+N-1) ^TＸ_L,p(k+N-1) ］^-1Ｅ_p(k+N-1) （19）Ｈ(k+N) ＝Ｈ(k) ＋μＸ_L,p(k+N-1) Ｇ(k+N-1) （20）出力ｙの計算は式（10）による１時刻の出力ｙ(k+1) 計
算の代わりに、次式でＮ時刻分ｙ(k+N) ，…，ｙ(k+2N-
1)を一度に計算する。

【００１７】Ｙ_N(k+2N-1)＝Ｘ_L,N(k+2N-1)^TＨ(k+N) （21）式（18）の右辺第２項の一部または全要素は式（21）で
ｋ＋２Ｎをｋ＋Ｎとした式で計算される。式（18），
（19），（20），（21）が従来のブロック射影法の計算
手順であり、図５Ｂにその機能ブロックを表わす。同図
において、プレフィルタ係数計算部２１、フィルタ係数
更新部２２、ブロックフィルタリング部４１、ブロック
誤差計算部４２はそれぞれ式（19），（20），（21），
（18）の演算を実行する。またこの従来のブロック射影
法の計算手順を図８に示す。

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

【００１９】Ｘ_L,N(k+N) ^TＨ(k+N) を効率的に計算す
るには、次式のように、Ｘ_L,NをＮ×Ｎの部分行列、Ｈ
をＮ要素の部分ベクトルに分解する。ただし、Ｌ＝ＮＭ
（Ｍは整数）となるように、Ｈに値が０の要素を加え、
Ｘを長くする。Ｘ_L,N(k+N) ^T=[Ｘ_N,N(k+N) ^T，Ｘ_N,N(k) ，…，Ｘ_N,N(k+N-(M-1)N)^T］^T Ｈ(k+N) ^T＝［Ｈ(k+N) ^T _0,Ｈ(k+N) ^T _1,…，Ｈ(k+N) ^T _M-1］^T (22) このように分解すると、Ｘ_L,N(k+N) ^TＨ(k+N) は、Ｘ
の部分行列、Ｈの部分ベクトル間の積の和として計算で
きる。

【００２０】Ｘ_L,N(k+N) ^TＨ(k+N) ＝ΣＸ_N,N(k+N-iN)^TＨ(k+N) _i（23） Σはｉ＝０からＭ−１まで Σの中の各々の正方行列とベクトル積は、入力Ｘ_N,N(k
+N-iN)とフィルタＨ(k+N) _iの畳み込みとみなすことが
できる。畳み込みを高速に行う方法として、ＦＦＴを用
いる方法、Winogradの高速畳み込み法に基づく方法（詳
しくは文献 Blahut, FAST ALGORITHMS FOR DIGITAL SIG
NAL PROCESSING, ADDISON-WESLEY, 1987，第３，９章）
がある。これらの、高速畳み込み法を式（23）に適用す
ることでＹ_N(k+N) の演算量が低減可能である。

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

【００２２】ｘ＊ｈ＝Ｆ^-1［Ｆ［ｘ＊ｈ］］＝Ｆ^-1［Ｆ［ｘ］★Ｆ［ｈ］］（25）ここで、Ｆ^-1［］はフーリエ逆変換を表す。フーリエ
変換に高速フーリエ変換（ＦＦＴ）を用いることで、高
速に畳み込みを行うことが可能である。式（20），（2
1）に適用すると次式になる。Ｈ(k+N) ＝Ｈ(k) ＋μＦ^-1［Ｆ［Ｘ_L,p(k+N-1) ］★Ｆ［Ｇ(k+N-1) ］］（26）Ｙ_N(k+2N-1)＝Ｆ^-1［Ｆ［Ｘ_L,N(k+2N-1)］★Ｆ［Ｈ(k+N) ］］（27）演算の低減量は、用いる手法、ブロック長Ｎ、フィルタ
長Ｌ、そしてハードウェアの機能に依存する。Ｎ＝数１
０、Ｌ＝数１００で、一般的なＤＳＰ（デジタルシグナ
ルプロセッサ）を用いた場合には、直接式（23）の左辺
を計算するのに比べて、Winogradの高速畳み込みに基づ
く方法で、約５０〜８０％の演算量となる。

【００２３】

【発明が解決しようとする課題】従来のブロック射影法
には、ブロック長Ｎを長くするにつれて、適応フィルタ
の収束が遅くなるという問題点がある。適応フィルタを
毎時刻修正する射影法に比べ、従来のブロック射影法で
は、ブロック長Ｎに一度しか、適応フィルタを修正しな
い。このために、適応フィルタの収束が遅くなる。

【００２４】

【課題を解決するための手段】この発明では、上記の問
題点に対して、高速射影法における、相関ベクトルＲ
_p-1(k)と平滑係数ベクトルＳ_p(k) の計算範囲を拡張す
ることで解決を図り、適応フィルタの収束速度の低下を
伴うことなく、射影法の演算量のブロック処理による低
減を可能にするものである。

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

【００２６】

【発明の実施の形態】以下、この発明を詳細に説明す
る。先に述べた毎時刻適応フィルタを修正する高速射影
法において、式（11）からＺ(k+i),(i>=0)はＺ(k) を用
いて表すことができる。Ｚ(k+1+i) ＝Ｚ(k+i) ＋ｓ_p(k+i) Ｘ_L(k+1+i-p) ＝Ｚ(k+i-1) ＋ｓ_p(k+i-1) Ｘ_L(k+i-p) ＋ｓ_p(k+i) Ｘ_L(k+1+i-p) … ＝Ｚ(k) ＋Σⁱ _j=0ｓ_p(k+i-j) Ｘ_L(k+1+i-p-j) ＝Ｚ(k) ＋Ｘ_L,i+1(k+1+i-p) Ｓ_p+i(k+i) _p,P+i （29）であり、特に、ｉ＝Ｎ−１のときに近似的フィルタ係数
Ｚ(k) からＺ(k+N) への更新式になる。

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

【００２８】また、拡張平滑係数ベクトルＳ_p+i(k+i) は、式（16）
のＳ_p-1(k+i) とＳ_p+i(k+i) _p,p+iを結合することで
次式で計算される。

【００２９】式（31）と平滑係数ベクトルの回帰式（11）を順に式
（32）に代入することで、拡張平滑ベクトルＳ_p+i(k+
i) の回帰式が得られる。

【００３０】ところで、式（14）において、ｋをk+i-1 とすること
で、次式が得られる。

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

【００３２】Ｒ_N+p-2(k+i) ^T _p,N+p-2＝Ｘ_L(k+i) ^TＸ_L,N-1(k+i-p) ＝Ｒ_N+p-2(k+i-1) ^T _p,N+p-2 ＋ｘ(k+i) Ｘ_N-1(k+i-p)−ｘ(k+i-L) Ｘ_N-1(k+i-p-L) （37）また、拡張相関ベクトルは式（17）のＲ_p-1(k+i)とＲ
_p+i-1(k+i) _p,p+i-1を結合することで計算される。

【００３３】拡張相関ベクトルの時間更新は、式（13）と式（37）に
おいてｉ＝Ｎ−１とした式を組み合わせて次のようにか
ける。Ｒ_N+p-2(k+i) ＝ｘ(k+i) Ｘ_N+p-2(k+i-1)−ｘ(k+i-L) Ｘ_N+p-2(k+i-L-1)（39) 拡張平滑係数ベクトル、拡張相関ベクトルを用いると、
式（36）は次のようにかける。

【００３４】ｙ(k+i) ＝ｙ′(k+i) ＋Ｒ_i+p-1(k+i) ^TＳ_i+p-1(k+i-1) （40）ｙ′(k+i) ＝Ｘ_L(k+1) ^TＺ(k) （41）ｙ(k) ，Ｅ_p(k) ，Ｇ(k) ，Ｓ_p(k) は高速射影法の手
順で計算できるので、ｉ＝０として式（40）を用いてｙ
(k) が計算できる。式(7.1),(7.2) を用いてｙ(k) から
ｅ(k+1) ，Ｅ_p(k) が計算でき、Ｅ_p(k) から式（８）
を用いてＧ(k)が計算できる。Ｇ(k) から、式（11），
（31），（32) を用いて、Ｓ_p+1(k)を計算でき、式（4
0）により、ｙ(k+1) が計算できる。このように、順に
ｙ(k+i) ，（ｉ＝０，１，…，Ｎ−１）は計算すること
ができ、時刻ｋ＋Ｎ−１において式（30）からＺ(k+N)
が計算できる。

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

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

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

【００３８】計算部５６よりの拡張相関ベクトルと計算
部５７よりの拡張平滑係数ベクトルとの内積演算、つま
り式（40）の右辺第２項の演算がなされる。平滑係数拡
張部５４よりの拡張部分ベクトルが近似的フィルタ係数
ブロック更新部５１へ入力されて式（30）が演算され
る。この演算結果の近似的フィルタ係数と入力ベクトル
との内積演算がブロックフィルタリング部４１で行わ
れ、その演算結果と、内積演算計算部３５からの演算結
果とが出力合成部３４で加算されて出力信号とされる。

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

【００４０】以上説明したこの発明法の特長は、まず、
近似的フィルタＺ(k) 、フィルタの出力ｙ(k) は、毎時
刻フィルタ係数を修正する高速射影法と同じ値になり、
従来のブロック射影法のようにフィルタの収束速度が劣
化しないことである、そして、演算量は、ブロック処理
であることにより、式（30）の右辺第２項と式（34）の
右辺第１項における、各入力信号からなる行列とベクト
ルの積に対して、従来のブロック射影法の説明で述べた
高速畳み込み算法を適用することで削減されることであ
る。実施例第１の例は、音響エコーキャンセラである。図３Ａはこ
の発明の方法を適用した音響エコーキャンセラの機能構
成を示す。入力信号ｘ(k) として音声信号がフィルタリ
ング部１１、フィルタ係数更新部１２、スピーカ６１へ
供給されマイクロホン６２の出力が所望信号ｄ(k) とし
て差回路１３へ供給される。スピーカ６１から音響信号
通路、いわゆるエコー伝達経路６３を通じてマイクロホ
ン６２へ達する信号が音響エコーであり、適応フィルタ
の出力信号ｙ(k) に推定エコー、推定誤差ｅ(k) に残留
エコーが対応している。

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

【００４２】図３Ｂは図３Ａの構成におけるこの発明の
効果を示すために行った実験の結果であって、縦軸がエ
コーの消去量（ｄＢ）であり、横軸が音響エコーキャン
セラが動作してからの時間である。曲線７１，７２，７
３はそれぞれ、高速射影法、従来のブロック射影法、こ
の発明によるエコー消去量の時間変化を示すものであ
る。ただし、高速射影法の曲線７１とこの発明法の曲線
７２は一致しているので、区別ができない。実験条件
は、入力信号は音声、適応フィルタ長Ｌ＝１０００，射
影の次数ｐ＝１６、ブロック長Ｎ＝１６である。結果の
曲線は５０回の試行の平均である。図３Ｂから、この発
明法は毎時刻適応フィルタを修正する高速射影法と同じ
収束特性を有し、従来のブロック射影法より、収束速度
が速く４秒以降の定常状態ではエコー消去量が大きいこ
とがわかる。これより、この発明法の従来法に対する優
位性がわかる。

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

【００４４】以上、２つの適用例を示した。どちらの場
合にもエコー、騒音を消去するという目的を達成するた
めには、推定エコー、推定騒音が得られれば十分であ
り、エコー、騒音の伝達経路そのものを求める必要はな
い。従って、どちらの適用例でも近似的フィルタ係数を
使用するこの発明法が適用でき、従来装置と比較して演
算量を削減できる。また、この発明は、図５Ａに示した
ような、構成において誤差信号ｅ(k) のパワーを最小と
する基本機能を有している。従って、解決すべき問題が
図５Ａに示した誤差信号のパワーの最小化としてモデル
化できる全ての応用例に対して、この発明を適用するこ
とが可能である。Ｎは２以上であればよく、例えば、３
６，６４などが考えられ、Ｎを大とする程、演算量が減
少するが、動作を開始してから最初の出力が現われる時
間が長くなる、よってこの発明を適用する対象に応じて
好ましいＮを決定することになる。

【００４５】

【発明の効果】以上説明したように、この発明では、適
応アルゴリズムの１つである射影法に拡張相関ベクトル
Ｒ_N+p-1(k)、拡張平滑化係数ベクトルＳ_i+p-1(k+i),(ｉ
＝１，…，Ｎ）を導入することにより、収束速度を低下
させること無く、ブロック処理による演算量低減を可能
とすることができる。その結果、この発明は収束の遅い
従来のブロック射影法等の手法が利用されていた音響エ
コーキャンセラ装置に対して、収束速度を低下させるこ
となく、ブロック処理による演算量の低減が可能とな
る。即ち、同一のハードウェア規模で収束は速いが多少
演算量の必要なより高い次数の射影法の利用、フィルタ
長を長くすることが可能となり、より高速なエコー低減
や定常時のエコーの低減の実現が達成される。

【図面の簡単な説明】

【図１】この発明の方法における処理手順を示す流れ
図。

【図２】この発明の方法を機能的に示すブロック図。

【図３】Ａは音響エコーキャンセラを示す構成図、Ｂは
音響エコーキャンセラに適用した従来、この発明の適応
フィルタの収束特性を示す図である。

【図４】騒音制御装置を示す構成図。

【図５】Ａは適応フィルタの一般的機能構成を示すブロ
ック図、Ｂは適応アルゴリズムに射影法を適用した適応
フィルタの機能構成を示す図である。

【図６】Ａは適応アルゴリズムに高速射影法を適用した
適応フィルタの機能構成を示すブロック図、Ｂは適応ア
ルゴリズムにブロック射影法を適用した適応フィルタの
機能構成を示すブロック図である。

【図７】適応アルゴリズムに高速射影法を適用した処理
手順を示す流れ図。

【図８】適応アルゴリズムにブロック射影法を適用した
処理手順を示す流れ図。

───────────────────────────────────────────────────── フロントページの続き (51)Int.Cl.⁷ 識別記号ＦＩＨ０４Ｂ 3/23 Ｈ０４Ｍ 1/60 ＣＨ０４Ｍ 1/60 9/08 9/08 Ｇ１０Ｋ 11/16 Ｈ (56)参考文献特開平７−92980（ＪＰ，Ａ) 特開平５−244043（ＪＰ，Ａ) 特開平６−132782（ＪＰ，Ａ) 特開平７−212278（ＪＰ，Ａ) (58)調査した分野(Int.Cl.⁷，ＤＢ名) G10K 11/178 G05B 13/02 G10L 21/02 H03H 17/02 601 H03H 21/00 H04B 3/23 H04M 1/60 H04M 9/08

Claims

(57)【特許請求の範囲】

【請求項１】入力信号を可変フィルタ処理し、その出
力信号と所望信号との誤差信号を求め、その誤差信号の
パワーを最小とするように、上記可変フィルタの特性を
適応的に制御する方法において、上記入力信号ベクトルの自己相関行列の逆行列と、上記
誤差信号ベクトルとの積に相当する演算を行ってプレフ
ィルタ係数ベクトルを求める第１ステップと、上記プレフィルタ係数ベクトルを平滑化して平滑係数ベ
クトルを得る第２ステップと、上記平滑係数ベクトルの要素を加重係数として、上記入
力信号ベクトルの方向に近似的フィルタ係数ベクトルを
修正する第３ステップと、上記入力信号ベクトルと上記近似的フィルタ係数ベクト
ルとの内積演算を行う第４ステップと、上記入力信号ベクトルの自己相関ベクトルを求める第５
ステップと、上記平滑係数ベクトルと上記自己相関ベクトルとの内積
演算を行う第６ステップと、上記第４ステップで得られた内積演算結果と、上記第６
ステップで得られた内積演算結果とを加算して上記出力
信号を得る第７ステップとを有するものであって、上記第３ステップの修正及び上記第４ステップの演算は
Ｎ時刻（Ｎは２以上の整数）ごとに行い、上記第３、第４ステップ以外の各ステップは各時刻ごと
に演算を行い、上記第３、第４ステップの演算が行われるごとに、それ
が次に演算されるまでの各時刻において、上記第２ステ
ップで得る上記平滑係数ベクトルの次数を１ずつ増加
し、次の第３、第４ステップの演算時刻での演算でもと
の次数に戻し、上記第５ステップでは上記平滑係数ベクトルの最大次数
と等しい次数の上記自己相関ベクトルを演算し、上記第６ステップでは上記平滑係数ベクトルの次数と等
しい次数だけ上記第５ステップで得られた自己相関ベク
トルから用いることを特徴とする適応的制御方法。
【請求項２】請求項１の適応制御方法において、上記入力信号をｘ(k) 、上記第６ステップでの相関ベク
トルをＲ_p+N-2(k+i)、上記第６ステップでの平滑係数ベ
クトルをＳ_p+i(k+i),（ｉ＝０，…，Ｎ−１）、その平
滑係数ベクトルのｐ番目からｐ＋Ｎ−１番目までの要素
からなるベクトルをＳ_p+N-1(k+N-1)_p,p+N-1、上記近似
的フィルタ係数をＺ(k) 、上記出力信号をｙ(k+i) とす
ると、上記出力信号、上記相関ベクトルの更新、上記平
滑係数ベクトルの更新、上記近似的フィルタ係数の更新
をそれぞれ次式ｙ(k+i) ＝Ｘ_L(k+i) ^TＺ(k) ＋Ｒ_p+i-1(k+i)^TＳ_p+i(k+i) Ｒ_p+N-2(k+i)＝Ｒ_p+N-2(k+i-1)＋ｘ(k+i) Ｘ_L(k+i-1) −ｘ(k+i-L) Ｘ_L(k+i-L-1) Ｚ(k+N) ＝Ｚ(k) ＋Ｘ_L,N(k+N-p) Ｓ_N+p-1(k+N-1)_p,N+p-1 によって実行することを特徴とする適応的制御方法。