第23回回路とシステム軽井沢ワークショップ The 23rd Workshop on Circuits and Systems in Karuizawa, April 19-20, 2010 高次の可変ディジタルフィルタを用いた狭帯域雑音除去システムの改良 Improvement of Narrow-Band Noise Suppression Systems Using a High-Order Variable Digital Filter 熊本裕樹越田俊介阿部正英東北大学大学院工学研究科電子工学専攻川又政征 Yuki KUMAMOTO Shunsuke KOSHITA Masahide ABE Masayuki KAWAMATA Department of Electronic Engineering, Graduate School of Engineering, Tohoku University 1 はじめにディジタルフィルタを用いた適応信号処理システムの例として，狭帯域雑音の除去システムがある．狭帯域雑音除去とは，広帯域の信号に対して狭帯域雑音が付加された場合に，狭帯域雑音を除去することである．狭帯域雑音の除去システムを実現するディジタルフィルタとして，適応ノッチフィルタ (Adaptive Notch Filter: ANF)[1] がよく用いられる．しかし，ANF を用いて狭帯域雑音の除去システムを実現する場合に，ある問題点が生じる．それは，単一正弦波以外の狭帯域雑音を除去する場合に、ANF では雑音を十分には除去しきれないことである．本稿では，単一正弦波以外の狭帯域雑音のことを帯域幅を持った雑音と呼ぶ．帯域幅を持った雑音の除去を，帯域阻止型適応フィルタを用いて実現できる．そして，帯域幅を持った雑音の除去性能を向上するために，帯域阻止型適応フィルタの遷移域の形状が急峻であることが望まれる．つまり，狭帯域雑音の除去システムを実現する帯域阻止型適応フィルタの次数を高くすることで帯域幅を持った雑音の除去性能を向上することができる．なぜなら，ディジタルフィルタには次数を高くすることで遷移域の形状が急峻になる性質があるからである．しかし，文献 [1] における ANF の次数は 2 次に固定されており，次数を高くすることができない．したがって，帯域幅を持った雑音を除去する場合には，次数を高くすることができる帯域阻止型適応フィルタを選択することが有効である．そこで，著者らはこれまでの研究において高次の可変帯域阻止フィルタ（Variable Band-Stop Filter: VBSF）を用いた狭帯域雑音の除去システムを提案している．文献 [2] において，6 次の VBSF を用いた狭帯域雑音の除去システムが提案しており，ANF[1] の次数よりも高い次数の帯域阻止型適応フィルタが用いられている．文献 [2] における VBSF の実現法は，周波数変換 [3] であり，プロトタイプフィルタの遅延素子を 2 次の全域通過フィルタで置換することで VBSF を実現している．そして，帯域幅を持った雑音が広帯域の信号に付加されている場合において，提案法の雑音除去性能が ANF よりも向上している結果が示されている．しかし，文献 [2] の手法では，VBSF の次数が高くなるにつれて VBSF のフィルタ係数の記述及び適応アルゴリズムの導出が困難になる問題がある．なぜなら，VBSF のフィルタ係数をプロトタイプフィルタのフィルタ係数と可変パラメータ α で記述する必要があり，次数が高くなるとこの記述が複雑化するからである．本稿では，これを解決する VBSF の実現方法を提案し，提案法に基づいた適応アルゴリズムを簡潔な形で導出する．まず，VBSF の内部状態を差分方程式で記述することで，高次の VBSF を容易に実現できる再帰的な式を導出する．次に，提案法に基づいた適応アルゴリズムを簡潔な形で導出する．最後に，提案法を用いた計算例を示す． 2 VBSF を用いた狭帯域雑音除去システム本章では，まず狭帯域雑音除去システムのモデルを示す．次に，VBSF について述べ，最後に狭帯域雑音除去システムを構築する際の一般的な高次の VBSF の実現法とその問題点について述べる． − 143 − [3]． z −1 ← T (z) 図 1: VBSF を用いた狭帯域雑音除去システムのモデル (2) 式 (2) において，帯域阻止フィルタを実現する T (z) は 2 次の全域通過伝達関数である．プロトタイプ低域通過フィルタの伝達関数 H(z) に対して，式 (2) を適用し，VBSF を実現する．したがって，プロトタイプ低域通過フィルタの伝達関数 H(z) と VBSF の伝達関数 HBS (z) との関係は以下のようになる．本稿におけるプロトタイプフィルタは，低域通過フィルタである． HBS (z) = H(z)|z −1 =T (z) (3) VBSF の阻止域幅は可変であるが，本稿では固定とする．そして，本稿における VBSF の阻止域幅は，プロトタイプフィルタの阻止域幅と等しくする．また，VBSF の中心周波数は可変パラメータにより変化させる．この要求を満たす VBSF を実現する T (z) は，以下の式で与えられる [4]． T (z) = VBSF を用いた狭帯域雑音除去システムのモデル VBSF を用いた狭帯域雑音除去を実現するシステムのモデルを図 1 に示す．図 1 において，時刻 n における入力信号 x(n) を次式で与える． x(n) = xs (n) + xv (n) (1) 式 (1) における xs (n) と xv (n) は，それぞれ広帯域信号と狭帯域雑音を表す．xs (n) は，平均 0，分散 σs2 の白色ガウス過程とする．xv (n) は，帯域幅 B の狭帯域雑音である．このようにして与えた入力信号の周波数スペクトルを図 2 に示す． 2.2 (4) 式 (4) における可変パラメータ α と VBSF の中心周波数 ω0 との関係式は以下のように与えられる [4]．図 2: 入力信号の周波数スペクトル 2.1 −αz −1 + z −2 1 − αz −1 VBSF α = cos ω0 (5) 本稿における VBSF は，式 (3) と式 (4) を用いてプロトタイプフィルタを周波数変換することで，阻止域幅固定，中心周波数可変の VBSF を実現できる． 2.3 高次の VBSF の一般的な実現法とその問題点本節では，VBSF の一般的な実現法について述べ， VBSF の次数を高くした場合の問題点を指摘する．本節では議論を簡単にするため，プロトタイプフィルタの次数は 2 次とする．つまり，プロトタイプフィルタの伝達関数を次式で与える． H(z) = b0 + b1 z −1 + b2 z −2 1 + a1 z −1 + a2 z −2 (6) 式 (6) に対して式 (4) の周波数変換を適用すると， VBSF の伝達関数 HBS (z) は次式で表される． VBSF は，可変フィルタに分類される．可変フィ HBS (z) = ルタは，少ないパラメータの調節で周波数特性を変化させることができるフィルタである．可変フィルタの周波数特性を変化させるパラメータは，可変パ = ラメータと呼ばれる． b0 + b1 T (z) + b2 T 2 (z) 1 + a1 T (z) + a2 T 2 (z) ( −2 ) ( −2 ) −1 −1 2 b0 + b1 z−αz−αz + b2 z−αz−αz −1 +1 −1 +1 ( ) ( ) −2 −1 −2 −1 2 1 + a1 z−αz−αz + a2 z−αz−αz −1 +1 −1 +1 本稿では，周波数変換 [3] を用いて VBSF を実現する．周波数変換は，以下の代入式に基づいている b′0 + b′1 z −1 + b′2 z −2 + b′3 z −3 + b4 z −4 (7) 1 + a′1 z −1 + a′2 z −2 + a′3 z −3 + a′4 z −4 − 144 − = 式 (7) における分母係数と分子係数はそれぞれ次式のようになる． a′1 = −α(a1 + 2) a′2 = a1 + α2 (1 + a1 + a2 ) a′3 = −α(a1 + 2a2 ) a′4 = a2 b′0 = b0 b′1 = −α(2b0 + b1 ) b′2 = b1 + α2 (b0 + b1 + b2 ) b′3 = −α(b1 + 2b2 ) b′4 = b2 (8) 式 (7)–式 (8) からわかるように，VBSF の一般的な実現法では周波数変換後の伝達関数 HBS (z) のフィルタ係数を，プロトタイプフィルタ H(z) のフィルタ係数と可変パラメータ α に基づく式で記述する必要がある．したがって，プロトタイプフィルタの次数が高次になるほど HBS (z) のフィルタ係数の記述が困難になる．これが適応アルゴリズムの導出を困難にする原因となる．なぜなら，適応アルゴリズムを導出する際に HBS (z) のフィルタ係数を α で微分する必要があり [5]，プロトタイプフィルタの次数が高くなるほど，この微分の式が非常に複雑になるからである．本稿では，これらの問題を解決する手法を提案する． 3 図 3(b) のブロック図に基づいて式 (7)–式 (8) の手順を踏まずに VBSF を実現する方法を導出する．まず，VBSF を実現するための変数の定義を以下のように与える．図 3(a) に示されるように，プロトタイプフィルタの伝達関数の分母係数と分子係数をそれぞれ，[1, a1 , a2 , · · · , aN ]，[b0 , b1 , b2 , · · · , bN ] とおく．図 3(b) に示されるように，時刻 n におけるサブフィルタの入出力の値を uk (n) とおく．ここで，k = 0, 1, · · · , N である．そして，周波数変換する際の全域通過フィルタの分母係数と分子係数を，式 (4) より [1, −α(n), 0] と [0, −α(n), 1] で与える．ここで，α(n) は時刻 n における可変パラメータの値を表す．上記の定義に基づいて，VBSF の記述法を導出する．k 番目のサブフィルタの出力信号 uk (n) は，入力信号が uk−1 (n) の 2 次の IIR フィルタの出力信号であることから，差分方程式で記述できる．したがって，uk (n) と y(n) は次式のように記述できる． u0 (n) = − 式 (3) で表される周波数変換は，プロトタイプフィルタの伝達関数 H(z) の z −1 を全域通過フィルタの伝達関数 T (z) で置換することを意味する．また，周波数変換をブロック図で考えると，プロトタイプフィルタの遅延素子が全域通過フィルタ T (z) で置き換えられることを意味する．本稿では，T (z) として式 (4) を用いるので，VBSF はプロトタイプフィルタの遅延素子を式 (4) で表される 2 次の全域通過フィルタに置換することで実現される．図 3(a) に，N 次のプロトタイプフィルタのブロック図を示す．図 3(b) に，図 3(a) の遅延素子を 2 次 am um (n) + x(n) (9) uk (n) = −α(n)uk−1 (n − 1) + uk−1 (n − 2) + α(n)uk (n − 1) (10) y(n) = VBSF の実現方法 N ∑ m=1 提案法本章では，式 (7)–式 (8) のような手順を踏まずに VBSF を実現する手法を提案する．この手法に基づく適応アルゴリズムを簡潔な形で導出する．本稿では，プロトタイプフィルタのフィルタ構造を直接型 II とする． 3.1 の全域通過フィルタ T (z) で置換したブロック図を示す．本稿では，図 3(b) における点線部で囲まれた部分をサブフィルタと呼び，サブフィルタの入出力の関係式を差分方程式で記述し，VBSF を実現する． N ∑ b0 um (n) (11) m=0 式 (9)–式 (11) を用いることで VBSF を実現できる．これより，式 (7)–式 (8) の手順を踏まずに，VBSF を実現できることが示された．式 (9)–式 (11) を用いることにより，高次の VBSF を容易に実現することができる．すなわち，式 (9)–式 (11) における N を変化させるだけで VBSF を実現できる． 3.2 適応アルゴリズム式 (9)–式 (11) で記述された VBSF を用いた狭帯域雑音除去のための適応アルゴリズムを導出する．本稿で導出する適応アルゴリズムは勾配法に基づいて [ ] おり，評価関数 E y 2 (n) を最小化するアルゴリズムである．ここで，E [∗] は期待値演算を示す． [ ] まず，E y 2 (n) における期待値演算を瞬時値 y 2 (n) に置き換える．すると，VBSF を実現するために必要な可変パラメータ α の更新式は次式のように記述 − 145 − (a) プロトタイプフィルタのブロック図 (b) VBSF のブロック図図 3: プロトタイプフィルタと VBSF のブロック図できる．と式 (10) をそれぞれ α(n) で微分する． α(n + 1) = α(n) − µy(n) ∂y(n) ∂α(n) (12) ここで，µ はステップサイズパラメータである． = − ∂uk (n) ∂α(n) = −uk−1 (n − 1) − α(n) ∂y(n) ∂α(n) = = = ∂ ∂α(n) N ∑ N ∑ ∂b0 um (n) + ∂α(n) m=0 N ∑ m=0 b0 ∂um (n) ∂α(n) m=0 b0 ∂um (n) (14) ∂α(n) ∂uk−1 (n − 1) ∂α(n) (18) 式 (18) において，uk−1 (n−1)，uk−1 (n−2) と uk (n−1) を α(n) で偏微分する項があるが，この演算をすることは困難なので，次の近似を適用する． ∂α(n − 1) ≈1 ∂α(n) (15) ここで，式 (14) から式 (15) を得るために次の性質を利用した． ∂b0 =0 ∂α(n) (17) (13) m=0 N ∑ ∂um (n) ∂α(n) + } b0 um (n) m=1 am ∂uk−1 (n − 2) + uk (n − 1) ∂α(n) ∂uk (n − 1) +α(n) ∂α(n) 式 (12) の右辺第二項における ∂y(n)/∂α(n) は，式 (11) より以下のように表される． { N ∑ ∂u1 (n) ∂α(n) 式 (19) を用いることで，次の近似式が成り立つ． ∂uk (n − 1) ∂α(n) (16) 式 (16) は，周波数変換前後ではプロトタイプフィルタのフィルタ係数は変化しないことから導かれる．式 (15) における ∂um (n)/∂α(n) を求めるために，式 (9) (19) = ≈ ∂uk (n − 1) ∂α(n − 1) ∂α(n − 1) ∂α(n) ∂uk (n − 1) (20) ∂α(n − 1) 式 (20) の近似式を利用することで，∂y(n)/∂α(n) を − 146 − 求めるために必要な式は以下のように記述できる． ≈ −uk−1 (n − 1) − α(n) ∂uk−1 (n − 1) ∂α(n − 1) ∂uk−1 (n − 2) + + uk (n − 1) ∂α(n − 2) ∂uk (n − 1) +α(n) ∂α(n − 1) ∂u1 (n) ∂α(n) ∂y(n) ∂α(n) ≈ − N ∑ m=1 ≈ N ∑ m=0 am ∂um (n) ∂α(n) ∂um (n) bm ∂α(n) (21) 30 20 Power[dB] ∂uk (n) ∂α(n) 40 10 0 −10 (22) −20 −30 0 1 40 Power[dB] 式 (21)–式 (23) より，プロトタイプフィルタの次数が高次であっても，式 (21)–式 (23) における N を変化させるだけで容易に適応アルゴリズムを構築できる．つまり，式 (21)–式 (23) を用いることによって，高次の VBSF に基づく適応アルゴリズムを容易に構築できる． 0.4 0.6 0.8 Normalized frequency 図 4: 入力信号の周波数スペクトル式 (21)–式 (23) を用いることにより，∂y(n)/∂α(n) が再帰的に求められる．このようにして得られた ∂y(n)/∂α(n) を式 (12) に代入することにより，可変パラメータが更新される． 4 0.2 (23) 20 0 −20 計算例 0 0.5 1 Power[dB] 本章では，提案法に基づく狭帯域雑音除去システ Normalized frequency ムを用いた狭帯域除去の計算例を示す．比較として， (a)VBSF の出力信号の周波数スペクトル ANF に基づく計算例も示す．ANF の適応アルゴリズムは，文献 [1] における Lattice Gradient Algorithm 40 を用いる．入力信号は，式 (1) で表される信号とする．xs (n) は，平均 0，分散 0.01 の白色ガウス過程 20 である．xv (n) は，平均 0，分散 1 の白色ガウス過程を帯域通過フィルタに通すことで生成した．ただ 0 し，xv (n) を生成するために用いた白色ガウス過程と xs (n) は互いに独立の白色ガウス過程である．xv (n) を生成するために用いた帯域通過フィルタは，中心 −20 周波数 0.3π[rad]，帯域幅 0.1π[rad] の 20 次バタワース型 IIR フィルタとした．つまり，狭帯域雑音の帯域 0 0.5 1 Normalized frequency 幅は，0.1π[rad] である．このようにして与えた入力信号 x(n) の周波数スペクトルを図 4 に示す．VBSF (b) ANF の出力信号の周波数スペクトルの次数は，6 次とした．ANF の次数は 2 次である [1]．図 5: 各フィルタの出力信号の周波数スペクトル VBSF と ANF の阻止域幅は，0.15π[rad] とした．ここで，VBSF と ANF の阻止域幅を狭帯域雑音の帯域幅と比べて 0.05π[rad] 広くとるのは，それぞれのフィ狭帯域雑音除去システムにおけるステップサイズパルタの遷移帯域幅が 0 ではないからである．VBSF ラメータ µ は，ともに 0.005 とした．VBSF と ANF を用いた狭帯域雑音除去システムと ANF を用いたの初期中心周波数は，ともに 0.5π[rad] とした． − 147 − Normalized center frequency 心周波数の平均をとったものを示している．図 6(b) も同様に，実験を 200 回独立に行い，ANF の中心周波数の平均をとったものを示している．図 6 より， VBSF と ANF ともに，中心周波数が狭帯域雑音の中心周波数に収束していることがわかる．収束速度は，本章で示した例では VBSF を用いた場合の方が速い結果が得られた．しかし，ANF を用いた場合の方が VBSF を用いた場合と比べて収束速度が速いことが多い． 0.5 0.45 0.4 0.35 0.3 狭帯域雑音の中心周波数 0.25 0 5000 Iteration 10000 Normalized center frequecy (a)VBSF の中心周波数の時間変化 0.5 0.45 0.4 0.35 0.3 狭帯域雑音の中心周波数 0.25 0 5000 Iteration 10000 5 まとめ狭帯域雑音の除去性能を向上するために，VBSF の次数を高くすることが求められる．しかし，文献 [2] では周波数変換により実現される VBSF の次数が高くなるほど，VBSF のフィルタ係数の記述が複雑になるという問題点があった．この問題点は適応アルゴリズムの構築も困難にするという問題があった．本稿では，VBSF の内部状態を差分方程式で記述することによって，高次の VBSF を容易に実現できる手法を提案した．そして，上記の提案法に基づく適応アルゴリズムを簡潔な形で導出した．最後に，狭帯域雑音除去の計算例を示し，高次の VBSF を用いた狭帯域雑音除去システムの有効性を示した．参考文献 (b) ANF の中心周波数の時間変化 [1] P. A. Regalia, “Adaptive IIR Filtering,” Marcel Dekker, 1995. 図 6: 各フィルタの中心周波数の時間変化図 5 は，それぞれ提案法に基づいた VBSF と ANF の各出力信号の定常状態における周波数スペクトルを表している．図 5 より，6 次の VBSF を適応フィルタとした方が ANF を適応フィルタとしたよりも狭帯域雑音が除去できていることがわかる．また，図 5(b) と図 5(a) を比べて図 5(b) の方が広帯域信号のパワーがより減衰されていることがわかる．広帯域信号は，保存されるべき信号なので減衰されてしまうことは望ましくない．これは，2 次の ANF の遷移帯域幅が 6 次の VBSF の遷移帯域幅と比べて広いことが原因である．これらの結果から，帯域幅を持った狭帯域雑音を除去する場合には，著者らが提案する高次の VBSF を用いた狭帯域雑音除去システムが有効であることがわかる．図 6 は，それぞれ提案法に基づいた VBSF と ANF の各中心周波数の時間変化と狭帯域雑音の中心周波数の関係を表している．図 6(a) における VBSF の中心周波数は，実験を 200 回独立に行い，VBSF の中 [2] 越田俊介, 阿部正英, 川又政征, “高次の可変ディジタルフィルタを用いた狭帯域信号検出システムに関する検討,” 平成 21 年度電気関係学会東北支部連合大会講演論文集, p. 266, Aug. 2009. [3] A. G. Constantinides, “Spectral transformations for digital filters,” Proc. IEE, vol. 117, no. 8, pp. 1585–1590, Aug. 1967. [4] 樋口龍雄, 川又政征, “ディジタル信号処理,” 昭晃堂, 2000. [5] K. Hashimoto and M. Kawamata, “Fast adaptive detection of sinusoidal signals using variable digital filters and all-pass filters,” Proceedings of IEEE International Symposium on Intelligent Signal Processing and Communication Systems, pp. 73–77, Nov. 2001. − 148 −