計測自動制御学会東北支部

第 295 回研究集会（2015.6.26）

資料番号 295-4

擬似白色化信号を用いたハウリングキャンセラに関する検討
A study on howling canceller using quasi- whitened input signals〇藤村達弘† 工藤憲昌
〇Tatsuhiro Fujimura††田所嘉昭‡

Norimasa Kudoh

八戸高専 機械・電気システム工学専攻

Yoshiaki Tadokoro‡

八戸高専‡豊橋技科大 工学部

Hachinohe National College of Tech. Toyohashi Univ. of Tech.
キーワード：適応フィルタ，周波数推定，補聴器

連絡先：〒039-1192 八戸市田面木上野平 16-1

八戸高専

電気情報工学科

tel:0178-27-7281， e-mail:kudohk-e@hachinohe-ct.ac.jp

1. はじめに補聴器を使用した際に，補聴器内のマイクに入った音は増幅されてからスピーカから出力される．Fig.1 は補聴器の構造を簡便に表したものである．スピーカとマイクの距離が非常に近い為に，スピーカから出た音の一部は外耳道を伝わり，再びマイクの入力となる．
これが続くことにより補聴器内に音響フィードバックループが形成され発振することで，ハウリングが発生する．ハウリングは，補聴器の使用者にとって非常に耳障りな音になる．
本稿では，このハウリングの除去を学習によって自動的に行うことを目的とした適応ハ
ウリングキャンセラを用いた，補聴器の品質改善について検討する．また適応ハウリング
キャンセラでの採用候補として，全域通過フィルタによる新しい周波数推定法を検証し，従来の周波数推定法との比較検討を行う．
以下に本稿の構成を示す．2.ではハウリン
グキャンセラおよびハウリングキャンセラを用いた従来法の補聴器のモデルについて述べる．2-1.では信号の相関とバイアスの関係について説明する．3.では提案法の補聴器のモデルについて概説する．4.では 3.で説明したモデルにおける，周波数推定部での周波数推定法の比較について，4-1.では，NF と BPF を組み合わせた周波数推定法，4-2.では，全域通過フィルタを用いた周波数推定法について説明
1

する．5.は数値例およびシミュレーション条件である．6.はまとめである．

Fig.1 hearing aids
2. ハウリングキャンセラ
本稿では，ハウリングを除去するために，補聴器内に学習機能を持ったハウリングキャンセラを設ける．ストレートフォーワードにハウリングキャンセラを設けた補聴器のブロック図を Fig.2 に示す．

Fig.2 hearing aids using howling canceller入力信号 d(n)は増幅部 G(z)で増幅され出力信号 x(n)となる．ここで出力信号の一部は音響フィードバック経路(以降，未知系 H(z)と呼ぶ)を通り，入力に戻る．y(n)は音響フィードバックを含んだマイクからの信号となる．ハ
ウリングキャンセラ(以降推定系 W(z)と呼ぶ)
は未知系 H(z)を同定し，y(n)に含まれる音響フィードバックを除去する．e(n)はハウリングを除去した後の信号である．
未知系 H(z)は補聴器の使用者の外耳道(耳の穴から鼓膜まで)の構造によって異なるため，推定系 W(z)は LMS(Least mean square)法によって H(z)を推定する必要がある．今回，未知系 H(z)を FIR フィルタの形として設計する．
ブロック図内の各点の式を示す．
x(n) = G(z)e(n)(1)y(n) = d(n) + H(z)G(z)e(n)(2)e(n) = y(n) - W(z)G(z)e(n)= d(n) + H(z)G(z)e(n) - W(z)G(z)e(n)= d(n) + (H(z) - W(z))G(z)e(n)(3)
式(3)よりd(n) = e(n){1 + (H(z) - W(z))G(z)}(4)
よって推定系 W(z)の周波数応答が未知系 H(z)と同じになれば d(n)=e(n)となり，ハウリングが除去されたことになる．これ以降，Fig.2 のような構成でハウリングを除去する方法を「従来法」と呼ぶ．
LMS 法は真値とその推定値との二乗誤差を最小化するように更新式を用いて係数を推定する方法で，この場合，Fig.2 の e(n)の 2 乗を最小にする．d(n)は補聴器の入力であり，y(n)，e(n)はそれぞれ式(2)，(3)で表わされる．LMS法では瞬時勾配を使って更新式を作る．式(2)，(3)を使うと瞬時勾配は式(5)のようになり，W(z)の係数の更新式は式(6)で表される．
𝜕𝑒(𝑛)2𝜕𝑒(𝑛)̂𝑤𝑙 (n) =∇= 2e(n)𝜕𝑤𝑙𝜕𝑤𝑙
𝜕{𝑦(𝑛) − 𝑊(𝑧)𝑥(𝑛)}= 2𝑒(𝑛)𝜕𝑤𝑙= −2𝑒(𝑛)𝑥(𝑛 − 𝑙)(5)̂𝑤𝑙 (n)𝑤̂ 𝑙 (𝑛 + 1) = 𝑤
̂ 𝑙 (𝑛) − 𝜇′∇= 𝑤
̂ 𝑙 (𝑛) + 𝜇𝑒(𝑛)𝑥(𝑛 − 𝑙)(6)
（但し，2 𝜇’ = 𝜇としている）
しかし，従来法では W(z)の係数の真値と推定値との間にバイアス（直流的誤差）が生じ，正しい推定を行うことが出来ない．これは，入力信号 d(n)と鼓膜への出力 x(n)に相関があるためである．

2-1. 信号の相関とバイアスx(n)は補聴器の出力であるが，学習における入力でもある．さらに d(n)は学習における観測雑音となる．例として，Fig.3 のような入力信号 d(n)を入力した場合，補聴器の出力信号x(n)は Fig.4 のような波形となる．
[dB]
f1=1000[Hz]f2=2000[Hz]

Amplitude-10-20-3000.1

0.2
0.30.4normalized frequency

0.5

Fig.3 frequency characteristic of input signal d(n)amplitude[dB]-20-3000.1

0.2
0.3
0.4
normalized frequency

0.5

Fig.4 frequency characteristic of output signal x(n)Fig.3,4 より出力信号は入力信号と周波数特性が似ていて，相関があることがわかる．
推定値と真値の誤差ベクトル C(n)の平均値の式を式(7)に示す．式(7)の E[ ]は平均操作を表し，R は H(z)，W(z)への入力 x(n)の自己相関行列，I は単位行列である．(I-μR)をδ，𝜇𝐸[𝑥𝑛 𝑑𝑛 ]をβとおくと，式(8)のようになる．nを無限としたものを式(9)に示す．第 1 項のδ，つまり(I-μR)は 1 未満であるため，第 1 項はゼロになる．しかし，第 2 項のβ，つまり x(n)d(n)に相関がある場合は，その平均値はゼロにならず式(10)のようになる．ここでδは 1 未満の1 に近い値なので，大きくなった第二項がバイアスとして残る．
E[𝐶𝑛+1 ] = (𝐼 − 𝜇𝑅)𝐸[𝐶𝑛 ] + 𝜇𝐸[𝑥𝑛 𝑑𝑛 ](7)
= δ𝐸[𝐶𝑛 ] + 𝛽(8)𝑛𝑖𝑛lim𝑛→∞ 𝐶n = lim𝑛→∞ (δ 𝐶0 + β ∑𝑖=0 δ )𝑛2𝛽= lim𝑛→∞ (δ 𝐶0 + 1−δ)

(9)

=(10)

𝛽1−δ3. 検討した方法
2-1.で述べたように従来法では d(n)と x(n)との間に相関があるため，推定系にバイアスが生じ，これにより H(z)の同定が難しくなる．
そのため d(n)と x(n)の相関を除去することで，推定値に生じるバイアスを低減させる方法が提案されている[3]．d(n)と x(n)の相関を低減させるため，入力に含まれている複数の周波数を周波数推定により推定し，FIR 形のノッチフィルタで周波数推定により得られた共振周波数を除去し，信号の持つ周波数を除去することで擬似白色化した信号で LMS 法を行う．これにより正確な H(z)のフィルタ係数の推定を行うことができる．文献[3]で提案された補聴器のブロック図を Fig.5 に示す．
Fig.5 の周波数推定部𝐻𝑁 (𝑧)において，入力信号に含まれる周波数を推定する．周波数推定の推定結果をコピーし，e(n)と x(n)に含まれる周波数をノッチフィルタによって除去し，入力信号を疑似白色化する事で相関を低減する．
ノッチフィルタを通った e(n)，x(n)を e’(n)，x’(n)とし，LMS 法の更新式に組み込むことによってできた修正 LMS 法の更新式を式(11)に示す．
𝑤
̂ 𝑙 (𝑛 + 1) = 𝑤
̂ 𝑙 (𝑛) + 𝜇𝑒′(𝑛)𝑥′(𝑛 − 𝑙)(11)amplitude[dB]20100-1000.1
0.2
0.3
0.4
normalized frequency

0.5

Fig.6 frequency characteristic of human voice
4-1. NF および BPF による周波数推定法
ノッチフィルタ(NF）とバンドパスフィルタ(BPF)をトリー型に配置した周波数推定部をFig.7 に示す[2]． NF および BPF の伝達関数
は式(13)，(14)で表される．
𝐻𝑁 (𝑧) =1 − 𝛼I 𝑍 −1 + 𝑍 −2
1 − 𝛾𝛼I 𝑍 −1 + 𝛾 2 𝑍 −2(13)𝐻𝑆 (𝑧) =−(1 − 𝛾)𝑍 −1 + 𝛾(1 − 𝛾)𝑍 −3
1 − 𝛾𝛼I 𝑍 −1 + 𝛾 2 𝑍 −2(14)フィルタのタップ係数𝛼I の真値は2cosωk とし，更新式は式(15)で表される．
(15)𝛼I (𝑛 + 1) = 𝛼I (𝑛) − 𝜇1 𝑒(𝑛)𝑠(𝑛)ここで，𝜇1 はステップサイズ(0<𝜇1 <1)であり，この値により更新量が決まる．e(n)および s(n)はそれぞれ𝐻𝑁 (𝑧)，𝐻𝑆 (𝑧)の出力である．

Fig.5 hearing aids using the method proposed in [3]Fig.7 tree connection

4.周波数推定法の比較
本稿では，3.で述べた周波数推定部𝐻𝑁 (𝑧)において，提案されている二つの周波数推定法について比較，検討を行う[1]，[2]．人間の肉声の周波数特性は Fig.6 のような，ピークを中心に少し広がった共振周波数(フォルマント)をもつため，文献[3]とは異なり，帯域幅を制御できる IIR 形のノッチフィルタを用いる．
本研究における最終目標は，人間の肉声を補聴器の入力信号とするものである．今回は前段階として Fig.3 のような特性をもつ信号でシミュレーションを行う．

4-2. 全域通過フィルタを用いた周波数推定法
全域通過フィルタを用いた 2 次適応ノッチフィルタの縦続接続構成は，狭帯域信号の数にかかわらず，いずれかの狭帯域信号を除去することが証明されている[1]．特に，縦続接続構成の最終段の出力信号で縦続各段のタップ係数を適応制御する方法は，ノッチフィル
タの正弦波信号の周波数に対するバイアスが小さいという利点を有している．
全域通過フィルタを用いる 2 次ノッチフィ
ルタの縦続接続からなる多重適応ノッチフィ
3ガウス性白色雑音∅(𝑛)の分散は𝜎 2 = 0.01とした．
式(19)で与えられる入力信号を用いて，二つの周波数推定法について正確性と収束速度の点で比較を行う．
比較を行う前に，二つの周波数推定法のタップ係数の更新式(15),(18)について，条件を同等としなければならない．それぞれの更新式の第二項は更新量および瞬時勾配の積で表されている．式(15)において更新量は𝜇1 である．
4𝛼式(18)において更新量は 2 (𝑛)+u2(𝑛)となる．

ルタの構成を Fig.8 に示す．この適応ノッチフ
ィルタの伝達関数は式(16)で表される．
𝑁1
𝐻𝑁 (𝑧) = ∏ [ {1 + 𝐻𝐴𝑖 (𝑧)}]
2(16)

i=1

ここで，𝐻𝐴𝑖 (𝑧)は Fig.9 で示すような 2 次全域通過フィルタであり，その伝達関数は式(17)で与えられる．
𝜌 2 + 𝑎𝐼𝐼 𝑧 −1 + 𝑧 −2
𝐻𝐴𝑖 (𝑧) =(17)1 + 𝑎𝐼𝐼 𝑧 −1 + 𝜌 2 𝑧 −2フィルタのタップ係数𝛼II の真値は−2cosωk とし，更新式は式(18)で表される．
𝛼II (𝑛 + 1) = 𝛼II (𝑛) −4𝛼
𝑇𝑟𝑎𝑐𝑒(𝐑)u12

今，入力信号の狭帯域信号数は 2 としたのでFig.10 においてフィルタの縦続接続数は 2 つとなり，各段のフィルタはそれぞれ周波数𝑓1および𝑓2 の信号を除去するノッチフィルタであるため，周波数推定が収束している場合，𝑢1および𝑢2 の振幅はそれぞれ 1.0，0.9 となる．
よって式(18)における更新量は，

𝐸[𝑒𝑁 (𝑛)𝐮 (𝑛 − 1)] (18)

ここで，𝑎はステップサイズ(0< 𝑎 <1)であり，この値により更新量が決まる．また，𝜌は伝達関数の極半径であり，この値により除去帯域幅が決まる．通常，鋭いノッチ特性を得るた
めに𝜌は 1 近傍の値(0< 𝜌<1)に設定する．

4𝛼4𝛼

4𝛼=
=u12 (𝑛) + u22 (𝑛) 1.0 ∗ 1.0 + 0.9 ∗ 0.9 1.81
= 2.21𝛼Fig.8 NFs connected in cascade

同条件の比較を行うためにそれぞれの更新量を等価と考えると，それぞれの更新量の間には，𝜇1 =2.21α という関係が成り立つ．
今回のシミュレーションにおいて縦続接続型の更新式のステップサイズ𝛼 = 0.005 と設定した．よって上式よりトリー型の更新式のステップサイズは𝜇1 = 0.01105 と設定した．
シミュレーション結果をそれぞれ Fig.10,11に示す．縦続接続型は推定し収束するまでに約 800 回の反復を要するのに対し，トリー型は約 300 回の反復で収束するため，収束速度が非常に速いと言える．また，どちらの方法も正確性は高く誤差はほとんど見られなかった．

Fig.9 2nd order APF(i-th stage)
5. 数値例
入力信号を Fig3 のような 2 本の線スペクト
ルと平坦な分布をもつ信号の和として式(19)で与える．ここで，∅(n)はガウス性白色雑音である．
𝑥(𝑛) = 1.0𝑐𝑜𝑠 (+0.9𝑐𝑜𝑠 (2𝜋𝑓1 𝑛2𝜋𝑓1 𝑛) + 1.0𝑠𝑖𝑛 ()𝑓𝑠𝑓𝑠

2𝜋𝑓2 𝑛2𝜋𝑓2 𝑛) + 0.9𝑠𝑖𝑛 () + ∅(𝑛)𝑓𝑠𝑓𝑠(19)

Fig.10 estimation performance oftree connection(𝑓1 = 1000[𝐻𝑧]，𝑓2 = 2000[𝐻𝑧]，𝑓𝑠 = 8000[𝐻𝑧])4Fig.11 estimation performance of
cascade connection
6. まとめ
今回は線スペクトルと白色雑音により音声信号を模した信号を用いたシミュレーションにより，補聴器のハウリングキャンセラに用いる周波数推定部の比較検討を行った．
今回検証した全域通過フィルタを用いた縦続接続型の多重適応ノッチフィルタは周波数の推定値に対するバイアスが小さいという利点を有しているが，収束速度がトリー型に比べると遅く，補聴器のハウリングキャンセラに用いると考えた時に，実用化することは難しいということがわかった．原因として，全域通過フィルタのノッチ周波数部分での振幅が完全に 0 になっていないことが考えられる．
全域通過フィルタの特性上，極半径𝜌を極めて1 に近い値にすることでノッチ周波数部の振幅を 0 に近似することができるが，収束速度が著しく遅くなり，改良するのは難しいと考えられる．
今後はトリー型の周波数推定を軸に，周波数の推定値に対するバイアスの最小化や，収束速度の改善を行っていく．
参考文献[1] 衣笠保智，“全域通過フィルタを用いた多重適応ノッチフィルタの適応アルゴリズムに関する検討” 電子情報通信学会技術研究報告． DSP， ディ
ジタル信号処理 102(334)， 49-53， 2002-09-18[2] J.F. Chicharo, T.S.Ng, “Gradient-based adaptive IIR
notch filtering for frequency estimation” TEEE Trans.
ASSP, PP.769-777, 1990
[3]H.SAKAI, “Analysis of an Adaptive Algorithm for
feedback cancellation in hearing aids for sinusoidalsignals”, Tech. report of IEICE,SIP2007-130

5