東海大学紀要情報通信学部
vol.3,No2,2010,pp.6-12

論文

緩やかな時変プライマリパスに対応した多入力多出力アクティブノイズキャンセリングシステム
藤野巖*1

A Multi-input Multi-output Active Noise Cancelling Systemfor Slowly Time-varying Primary PathbyIwao FUJINO(Received on Oct. 7, 2010 and accepted on Mar. 1, 2011)

Abstract
In this paper, we present a new structure for multi-input multi-output (MIMO) active noise cancelling (ANC) system.
With this structure, it becomes possible not only to achive good noise cancelling performace, but also to deal with a slowly
time-varying primary acoustic propagation path. In our study, the MIMO ANC system structure is described in state spacemodel. For a feedforward type MIMO ANC system, we divided time signals into block data sequentially, and for eachblock, we estimated system of primary acoustic propagation path by using subspace identification method. Also in order toeliminate the effect of secondary acoustic propagation path, we connected this estimated system to the inverse system of it.
Eventually we can get the control sources for noise cancelling at cancelling points. By computer simulation of a numericalexample, we showed that the proposed MIMO ANC structure works with good performance.

Keywords: ANC, MIMO, Subspace Method, System Identification
キーワード：ANC、MIMO、部分空間法、システム同定
多入力多出力 ANC システムの構成手法は、大きく分けて、適応フィルタを用いて制御信号を作り出す適応フィルタ法と伝達関数の推定よって制御信号を作り出す伝達関数法に大別される。前者では、残留騒音が最小になるように、適応アルゴリズムに基づき、適応システム内部の重みパラメータを繰り返し調整する 1)2)3) 。この構成の最適解および適応アルゴリズムについて、S. J. Elliott が文献 4)において詳細に述べている。後者では、入出力信号のデータに基づき、システムの伝達関数を推定し、それによって、騒音を打ち消すための信号を作り出すことができる。このような構成の報告として、金田らの室内伝達関数を利用したものがある 5)6)。また、著者は以前、部分空間法を用いてシステムを同定する手法によって、高精度に同定作業を行うことができ、より残留騒音の少ないANC システムを実現できることについて確認した 7)。
本論文はシステムの状態空間表現に基づき、 部分空間システム同定法を用いた多入力多出力 ANC の実現法を提案するものである。本研究では、時系列信号を時間ブロックに分割し、それぞれの時間ブロックにおいて、部分空間法を用い、一括処理によってプライマリパスのシステム同定を行う。
このような手法により、従来の適応システム型よりも残留騒音

1．まえがきアクティブノイズキャンセラー(ANC)は電子設備によって能動的に騒音を打ち消す技術である。初期の研究の多くは、音響が一次元的に伝搬し、騒音源及び制御点がともに一個の場合についてモデル化されており、空調ダクトの吹き出し口、ヘッドフォン、携帯電話などへの応用例が数多く報告されている。近年、屋外、室内などの 3 次元空間におけるノイズキャンセリング問題が多く取り上げられ、基本原理と実用開発の両面で盛んに研究が行われている。3 次元空間領域においてノイズキャンセリングを行う場合、複数個の参照マイクと複数個の制御音源が必要となる。複数組の 1 入力 1 出力のシステムを用いることも考えられるけれども、多入力多出力制御システムを用いて、全体的に大局的に多入力多出力のノイズキャンセリングシステムを構成すれば、より良い消音性能が得られるとが考えられる．しかし、その反面、ノイズキャンセリングのアルゴリズムやシステムの実現法は複雑になりがちで、現時点では、十分に満足できる性能まで至ってない。

*1 情報通信学部通信ネットワーク工学科教授―1―−6−

藤野巖

が少なくなるばかりでなく、プライマリパスのシステム特性の緩やかな変動にも追従できる。本文は以下のように構成される。
まず、本研究で対処したい問題を整理し、状態方程式で定式化する。そして、本文で提案する部分空間同定法の利用法を述べ、ANC システムの全体構成とシステム動作の仕組みを示す。その上で、具体的な数値例を用いて、提案したシステムの性能検証を行う。最後に本研究で得られた結論をまとめる。

g(k)として消音ポイントでプライマリ騒音 d(k)を打ち消す。また、本文では、[・]を用いて、ベクトルやマトリクスの次元を示す。例えば、Fig.2 にある d(k):[N]は、d(k)は N 次元ベクトルであることを表す。
Fig.2 のコントローラの実現法は、主に適応フィルタ法と伝達関数法の 2 種類に分別される。以下にそれぞれの概要について述べる。

2．2 コントローラの実現法１：適応フィルタ法
適応フィルタ法の基本構成図を Fig.3 に示す。

2．多入力多出力 ANC システムの構成
2．1 多入力多出力 ANC システムの基本構成
まず１例として、2 参照入力(マイク)信号、2 制御出力(スピーカ)信号からなる ANC システムの配置イメージ図を示す。

Fig.3 Basic structure of feedforward MIMO ANC systemusing adaptive filter as controller
ここでは、W:[M,L]は L 入力 M 出力の適応フィルタの重みパラメータマトリクスを表す。また、HP(z):[N,L]は L 入力 N 出力のプライマリパスの伝達関数マトリクスを表し、HS(z):[N,M]はM 入力 N 出力のセカンダリパスの伝達関数マトリクスを表す 。
適応フィルタ法の場合、コントローラとして適応フィルタを用いる。適応アルゴリズムにしたがって、適応フィルタ重みパラメータを逐次的に調整することによって、残留騒音を最小にすることができる。最も一般的な適応 フィルタの構成は FIRフィルタである。1 チャンネルの騒音制御を行う場合、比較的に高速で、正確な制御を実現することが可能である。しかし、
L 入力 M 出力の重みパラメータマトリクス W:[M,L]の場合、M×L 個の I タップ長の FIR フィルタが必要となり、重みパラメータの数は全部で、M×L×I 個になる。この膨大な数の重み係数を更新するには、より多くの処理時間が必要となる。しかも、収束速度と残留騒音は参照入力間の相関特性に大きく依存する結果が報告されている。また、それを対処するために、白色化と呼ばれる参照入力信号の間の相関を取り除くような前処理が必要となる。適応フィルタ法では、最初からシステムの構成は FIR 形や IIR 形に限定されている。推定対象のシステム構成と一致しない可能性がある。それにもかかわらず、それを前提として、特定の評価関数を最小にするようなアルゴリズムでシステム同定を行っているから、真の最小残留誤差に到達できなく、推定精度が低い。しかも、同定対象が多入力多出力システムの場合、特定の入力から特定の出力までの伝達特性を個別に推定している。

Fig.1 An arrangement image of 2 input 2 output ANC system(dot line: acoustic propagation, solid line: electrical signal)
消音領域にスピーカとマイクを 2 組配置する。外部から来るプライマリ騒音をマイク M1 と M2 で受け取り、スピーカ S1 と S2からそれを打ち消すように、制御音 y1 と y2 を出力する。ここでは、本来の騒音の伝搬経路をプライマリパスと呼ぶ。また、制御スピーカとマイク間の伝搬空間はセカンダリパスと呼ぶ。このような関係をブロック図で表すと以下のようになる。

Fig.2 Basic structure of feedforward MIMO ANC system
Fig.2 において、騒音源からの騒音 u(k)はプライマリパスを経由して、プライマリ騒音 d(k)として消音ポイントにに到達する。
一方、騒音源から参照信号 u(k)を採集し、それを電気信号としてコントローラに入力する。コントローラからスピーカで制御音 y(k)を出力して、その制御音はセカンダリパスを経由して、―2―−7−

2．3 コントローラの実現法２：伝達関数法
伝 達 関 数 法 の 基 本 構 成 図 を Fig.4 に 示 す 。 従 来 で は 、ANC 問題について周波数領域で考えられていた。理想的にノイズキャンセリングが完全に行われたと仮定すれば、 Fig.4
から次式のように、プライマリパス伝達関数 HP(z)はコントロー

緩やかな時変プライマリパスに対応した多入力多出力アクティブノイズキャンセリングシステム

入力多出力システムを一体的に推定でき、非常に高い推定精度に到達できる。前述の適応フィルタ法および伝達関数法
よりも良い性能の ANC が実現できると考えられる。

3．2 提案したシステム構成
本研究は、2.3 で述べた伝達関数法の基本的な考え方に基づきながら、状態方程式表現を用いて、Fig.2 の ANC システムを多入力多出力システムに拡張したものである。具体的には、プライマリパス、セカンダリパスとコントローラの 3 つ構
成要素について、従来の伝達関数マトリクスの代わりに、状態方程式を用いて表現した。こうすることによって、一般的な多入力多出力システムに対応できるようになり、一体的にコントローラを導出して実現することができるようになる。更に、線形システムは一般性を失うことがなく、全域通過システムと最小位相システムの縦続結合として表せる 11)から、セカンダリパ
スを全域通過システム Gall(z)と最小位相システム Gmin(z)に分けることにする。ただし、離散システムの場合、全域通過システムは以下のような固定遅延要素の対角行列となる。また、最小位相システムについて、その逆システムが存在することが保証される。

Fig.4 Basic structure of feedforward MIMO ANC systemusing transfer function description for controller
ラの伝達関数 HC(z)とセカンダリパス伝達関数 HS(z)の積と等しいことが分かる。

H P ( z ) = HC ( z )H S ( z )(1)

そして、セカンダリパスの伝達関数を既知とした場合、形式的にコントローラの伝達関数は次のように表される。

HC ( z ) = H P ( z )H −S1 ( z )⎛ z − Nd1⎜
⎜ 0
G all ( z ) = ⎜⎜ M
⎜ 0
⎝

(2)

この式から 1 入力 1 出力 ANC システムの伝達関数法の基本的な考え方が提案されている 8)。しかし、多入力多出力システムの場合、それぞれの伝達関数はマトリクス関数であり、その逆マトリクスは存在しないケースが多い。 また、伝達関数法
では、多入力多出力システムの場合、個別の伝達関数を要素とする伝達関数マトリクスを用いてシステムを表すので、システム全体を推定するというよりも、伝達関数マトリクスの個別要素を推定している。

zM00 ⎞⎟L0 ⎟⎟OM ⎟ (4)0 z − N dm ⎟⎠L

本論文における問題の設定を以下に示す。
1) プライマリパスは緩やかな時変系とする。
2) セカンダリパスは線形時不変系で、既知とする。
3) 観測データをブロックに分割して処理する。
4) 各ブロック内において、プライマリパスを線形系として扱う。
5) 消音ポイントにおいて、プライマリパス経由の騒音以外にも、観測雑音がある。
上で述べた事項を踏まえ、本研究で用いる ANC システムの構成図を以下のようにまとめる。

3． 提案したシステム構成
3．1 部分空間法によるシステム同定
部分空間法は近年最も注目されているシステム同定法である。部分空間法では、以下の状態方程式によるイノベーションモデルを用いてシステムを表現する。

x ( k + 1) = Ax ( k ) + Bu( k ) + w ( k )y ( k ) = C x ( k ) + Du ( k ) + v ( k )0−Nd 2(3)ただし、u(k)は入力信号、y(k)は出力信号、ｘ(k)は内部状態である。ｗ(k)はプロセス騒音で、v(k)は消音ポイントにおける付加騒音である。A、B、C、D は適切な次元を持つ定数行列である。部分空間法を用いて、有限個入力データとそれに対応する出力データから、一括処理で同定ターゲットシステムの状態方程式の係数行列 A、B、C、D を推定できる。代表的な推定アルゴリズムとして、M. Verhaegen と P. Dewilde による
MOESP 法 9)および P. Van Overschee と B. De Moor による
N4SID10) などがある。
適応フィルタ法や伝達関数法と比べて、部分空間法では一般的な状態方程式表現を用いてシステムを表すので、線形系であれば、どんなシステム構成にも対応できる。その上、多Fig.5 Proposed MIMO ANC system using system
identification―3―−8−藤野巖

DC = D S

部分空間法を用いてプライマリパスの推定ができ、以下のような状態方程式で表されるものとする。

~~x P ( k + 1) = A P x P ( k ) + B P u( k )~~y ( k ) = C P x ( k ) + D P u( k )(16)

本研究では、この状態方程式のシステムを用いて、ブロックごとに制御音 y(k)を生成し、その制御音はセカンダリパスを通過して g(k)=q(k)となり、消音ポイントにおいてプライマリ騒音を打ち消すことができる。
(5)

3．3 データブロックの実現
また、セカンダリパスは既知とし、以下のような状態方程式で表されるものとする。

部分空間法はデータセットを一括処理する方式を取っている。同定対象のシステムの特性が変われば、新たにデータセットを用意してシステム同定を行う必要がある。本研究では 、Fig.6 のように、入力信号と出力信号の時系列データを長さ Lのブロックに分割している。各データブロックについては部分空間法による一括処理を行い、システムの係数行列A、B、C、D を推定する。全体的には、ブロックごとに、繰り返しシステム推定を行うことにしている。このような方法を取ることによって、システム特性の緩やかな変化にも追従できるようになり、より実際問題に対応できると考える。

x S ( k + 1) = A S x S ( k ) + B S u( k )y ( k ) = C S x S ( k ) + D S u( k )(6)
この式から、直接項の係数行列 DS の逆行列 DS-1 が存在すれば、逆システムの状態方程式は次のように導出できる。

x S (k + 1) = ( A S − B S D −S1C S )x S (k ) + B S D −S1y (k )u(k ) = − D −S1C S x S (k ) + D −S1y (k )(7)つまり、逆システムの状態方程式の各係数行列は以下のようになる。

A S = A S − B S D −S1C S

B S = B S D −S1−1S

C S = −D C S

DS = D−1
S

(8)
(9)
(10)
(11)

Fig.5 からコントローラは推定したプライマリパスのシステムとセカンダリパスの最小位相構成の逆システムとの縦続接続からなるので、その状態方程式の係数行列は次のように算出できる。

Fig.6 Data block for using subspace identification method

４．シミュレーション例⎛ x P ( k + 1) ⎞⎛ x (k ) ⎞⎜⎜⎟⎟ = A C ⎜⎜ P ⎟⎟ + BC u( k )⎝ x S ( k + 1) ⎠⎝ x S (k ) ⎠⎛ x (k ) ⎞y ( k ) = CC ⎜⎜ P ⎟⎟ + DC u( k )⎝ x S (k ) ⎠提案したシステム構成の性能を検証するために、具体的なプライマリパスとセカンダリパスの数値例を作成し、コンピュータシミューレーションを行った。その条件設定およびシミューレーション結果を以下に示す。

4．1 条件設定(12)

プライマリパスのシステム構成
この数値例では、以下のような係数行列の状態方程式で表されるプライマリパスを用いた。

上記状態方程式の係数行列は以下のようになる。

~⎛ AP0 ⎞⎟
Α C = ⎜⎜ ~
⎟BC
APS ⎠
⎝ S~
⎛ B ⎞
BC = ⎜⎜ ~P ⎟⎟⎝ BS DP ⎠~⎛ D S CP 0 ⎞⎟CC = ⎜⎜C S ⎟⎠⎝ 0

⎛ 0.5296 − 0.4760 0.1238 ⎞⎜⎟Α P 0 = ⎜ − 0.4760 − 0.0974 0.1354 ⎟⎜ 0.12380.1354 − 0.8233⎟⎠
⎝(13)

(17)

(14)

⎛ − 1.1470 − 0.0376 ⎞⎜⎟B P = ⎜ 1.1910
0.3273 ⎟
⎜
⎟
00⎝
⎠

(15)
―4―−9−(18)緩やかな時変プライマリパスに対応した多入力多出力アクティブノイズキャンセリングシステム⎛ − 0.1867 − 0.5883 − 0.1364 ⎞C P = ⎜⎜⎟
00.1139 ⎟⎠⎝ 0.7258(19)

0 ⎞⎛1.0670
D P = ⎜⎜⎟
0.8000 ⎟⎠⎝ 0(20)

ただし、係数行列 AP 固有値 λi(k)は、この係数行列 AP0 の固有値 λi0 をもとに、総数 Nt 個のデータの間に以下の式にしたがって緩やかに変化するものとする。
Fig.8 Block diagram of secondary acoustic propagation pathN⎛⎞⎜k− t π ⎟2 ⎟λi (k ) = λi 0 cos⎜
N3⎟⎜t⎟⎜⎝ 2⎠

上図から、この部分の入出力関係を以下のように表すことができる。

g1 (k ) = a11 y1 (k − N d ) + a12 y2 (k − N d − 1)
g 2 (k ) = a21 y1 (k − N d − 1) + a22 y1 (k − N d )実際のシミュレーションにおいて、以下のような手順で係数行列 AP0 から時刻 k における係数行列 AP(k)を作っている。

A P 0 = PΛ P 0 P −1N⎞⎛
⎜k− t π ⎟2 ⎟Λ P (k ) = Λ P 0 cos⎜
N3⎟⎜t⎟⎜⎠⎝ 2(24)

そして、セカンダリパスを以下のように、全域通過部分と最小位相部分に分けることができる。
全域通過部分：(21)

f1 (k ) = y1 (k − N d )f 2 (k ) = y2 (k − N d )

(25)

(22)
最小位相部分：A P (k ) = PΛ P P−1(23)g1 (k ) = a11 f1 ( k ) + a12 f 2 (k − 1)
g 2 (k ) = a21 f1 ( k − 1) + a22 f 2 (k )係数行列は初期値の AP0、BP、CP、DP を用いた場合、このプライマリパスの周波数特性を以下のボード図で示す。

(26)
この最小位相部分の入出力関係を、状態方程式を用いて表すと、その係数行列は以下のようになる。

⎛0A S = ⎜⎜⎝0
⎛1B S = ⎜⎜⎝00⎞⎟0 ⎟⎠(27)0⎞⎟1 ⎟⎠(28)⎛ 0 a12 ⎞⎟⎟C S = ⎜⎜⎝ a21 0 ⎠(29)0 ⎞⎛a⎟⎟
D S = ⎜⎜ 11⎝ 0 a22 ⎠Fig.7 Bode graph of primary acoustic propagation path(30)

さらに、この最小位相システムの逆システムを導出できる。その結果、状態方程式の各係数行列は以下のように示される。

セカンダリパスのシステム構成
コントローラから出力される制御音はセカンダリパスを経由して消音ポイントに到達する。本例では、Fig.8 のような構成のセカンダリパスを用いる。

⎛ a a −1
A S = −⎜⎜ 12 11⎝ 0―5―− 10 −⎞⎟a a ⎟⎠0−1
21 22(31)藤野⎛ a −1 0 ⎞⎟B S = ⎜⎜ 11−1 ⎟⎝ 0 a22 ⎠⎛a aC S = −⎜⎜⎝ 0−112 11⎛aD S = ⎜⎜⎝ 0−111巖(32)⎞⎟a a ⎟⎠0−121 220 ⎞
⎟
−1 ⎟
a22
⎠

(33)

(34)

評価パラメータ
ANC の消音精度の数値評価パラメータとして、ブロックごとに以下のような残留騒音のパワーとプライマリ騒音のパワーの比の dB 値を用いる。
LL∑
∑
k =1
J = 10 log10 kL=1
L
2g(k)+1∑∑ g22 ( k )k =1e12 ( k ) +Fig.10 Primary noise (red line) and residual noise (black line)of microphone 2

e22 ( k )

k =1

(35)

4．2 シミュレーションの結果
上記条件のもとで、MATLAB のプログラムを作成し、コンピュータシミュレーションを行った。実際のシミュレーションでは、計算時間の削減とプライマリパスの変化ペースに合わせて ANC システムが作動できれば良いと考えた。具体的に、プログラム上では残留騒音の評価パラメータ J について許容指標を設定した。時間ブロックごとに評価パラメータ J のチェックを行い、許容指標以下になったら、部分空間法によるプライマリパスのシステム同定を行うことにした。以下に示す例では、参照入力は 2 チャンネルの白色ノイズで、観測雑音も 2 チャンネル白色ノイズである。また、SN 比は 36dB、サンプリング周波数は 1000Hz である。データのブロック長は1000 個で、データの総数は全体で 1 分間の信号を構成するように 60000 個とした。さらに、セカンダリパスのパラメータa11=a22=1.0、a12=a22=0.8 とした。以下のグラフは今回のシミュレーション例の結果を示す。Fig.9 はマイク 1 のプライマリ騒音と残留騒音の時間信号である。Fig.10 はマイク 2 のプライマリ騒音と残留騒音の時間信号である。Fig.11 は各時間ブロックにおける評価パラメータ J の値である。

Fig.11 Ratio of primary noise power and residual noise power

5．むすび
本文は、部分空間同定法を用いて、緩やかな時変プライマリパスに対応した ANC システムについて考察した。状態方程式を用いて、ANC システムのプライマリパス、セカンダリパスおよびコントローラを表現することによって、部分空間同定法を導入することができた。また、時間信号を時間のブロックに分割した。各時間ブロック内では、部分空間法を用いて、一括処理によるプライマリパスのシステム同定ができ、非常に良い精度でコントローラを構成することができた。その上、時間ブロックごとに繰り返しプライマリパスのシステム同定を行うことで、プライマリパスの緩やかな変化に追従できるようになった。本研究の性能を検証するために、2 入力 2 出力の ANCシステムの構成例についてコンピュータシミュレーションを行った。その結果、残留騒音対プライマリ騒音の電力比は、全データの時間範囲において幾分の上下があるが、平均で約-26dB、最小で約-21dB、最大で約-36dB の騒音低減を達成することができた。一概に比較できるものではないが、プライマリパスに時不変システムモデルを用いた適応フィルタ法
を 適 用 し た 2 入力 2 出 力 ANC シ ス テ ム構 成 の 例 で は 、-18dB の騒音低減を達成している 12)。これに対して、本研究では、緩やかな時変プライマリパスや観測騒音などマイナス条件にも関わらず、性能的により良い結果が確認された。

Fig.9 Primary noise (red line) and residual noise (black line)of microphone 1―6―− 11 −緩やかな時変プライマリパスに対応した多入力多出力アクティブノイズキャンセリングシステム

4)

S. J. Elliott, “Signal processing for active control”, Chapter 5,
Academic Press, 20015) M. Miyoshi and Y. Kaneda, “Active control of broadband
random noise in a reverberant three-dimensional space,” NoiseControl Eng. J., 36, 85-90, 1991
6) Y. Haneda, S. Makino and Y. Kabeda,”Common Acoustical Poleand zero Modeling of room transfer functions,” IEEE Trans.
Speech and audio processing, vol.2 no.2, pp320-328, 1994
7) I. Fujino, “A multi-input multi-output active noise cancelling
structure using subspace system identification method,”
Proceedings of the 2009 International Symposium on Active
Control of Sound and Vibration, 20098) 梁載冕, 南炫道, 西山栄枝, “ダクト内における能動騒音低減のための適応アルゴリズム ,” 電気設備学会誌,14, 426-433,
19949) M. Verhaegen and P. Dewilde, “Subspace model identification
Part 1. The output-error state-space model identification class ofalgorithm,” Int. J. Control, 56(5) , 1187-1210, 1992
10) P. Van Overschee and B. De Moor, “N4SID: Subspacealgorithms for the identification of combined deterministicstochastic system,” Automatic, 30(1), 75-93, 1994
11) A. V. Oppenheim and R. W. Schafer, “Digital signal processing”,
pp.348-349, Prentice-Hall, New Jersey, 1975
12) S. J. Elliott, “Optimal controllers and adaptive controllers for
multichannel feedforward control of stochastic disturbances,”
IEEE Trans. Signal Processing, vol.48(4), 1053-1060, 2000

謝辞
本研究の取り組みにあたって、英国サウサンプトン大学(University of Southampton) 音 響 振 動 研 究 所 (ISVR:
Institute of Sound and Vibration Research)の S.J. Elliott 教授から大変刺激となるアドバイスを頂き、ここに感謝の意を表します。

参考文献
1)

2)
3)

S. J. Elliott, I. M. Stothers and P. A. Nelson, “A multiple error
LMS algorithm and its application to the active control of soundand vibration,” IEEE Trans. Acoust., Speech Signal Processing,
vol.ASSP-35, 1423-1434, 1987
D. Guicking and M. Bronzel, “Multichannel broadband active
noise control in small enclosures,” Proc. Inter-noise , 12551258, 1990
S. M. Kuo and B. M. Finn, “A general multi-channel filteredLMS algorithm for 3-D active noise control systems,” Proc. 2ndInt. Conf. Recent Developments in Air- And Structure-Borne
Sound Vib., 345-352, 1992

―7―− 12 −