ディジタルフィルタ初歩*1 - Kimura

ディジタルフィルタ初歩*1
2008/05/16
鳴海 智博
フィルタ解析
ディジタルフィルタの種類
一般的なディジタルフィルタの入力 x と出力 y の関係は，
y(k) + b1 y(k − 1) + · · · + bM y(k − M ) = a0 x + a1 x(k − 1) + · · · + aN y(k − N )
(1)
と表せる．
b1 · · · bM が存在するときには，過去の出力が回帰的に使用され，回帰フィルタ（recursive filter）と呼ばれる．これら
が全てゼロのときには，過去の入力のみが使用され，非回帰フィルタ（nonrecursive filter）と呼ばれる．なお，フィル
タ係数が有限個の場合には，回帰フィルタでは一般に無限の過去からの情報が生き続けるので，IIR フィルタ（infinite
impulse response filter）と呼ばれ，非回帰フィルタではフィルタ係数の個数より以前のサンプルのデータは計算に用い
られず出力に影響しないので，FIR フィルタ（finite impulse response）と呼ばれる．
回帰フィルタは少ない個数の係数でシャープなカットオフを実現できるなど有利な点も多いが，安定性が問題になっ
たり，波形の観察に重要な属性である線形位相遅れフィルタが簡単に実現しにくいなどの問題点がある．以下ではわか
りやすい非回帰フィルタのみに限定して考える．
フィルタと z 変換
ディジタルフィルタは，その係数の組，
a 0 , a1 , · · · , aN
(2)
で定義される．測定値に対してフィルタを適用することは，
y(k) =
N
ai x(k − 1)
(3)
i=0
という畳み込み演算を行うことである．z 変換を適用すれば，時間領域の畳み込み演算に対応する z 領域での伝達関数
H(z) を考えることが出来る．z 変換では，現在よりも k サンプル前のデータは z −k を掛けることで表されることを用
いれば，H(z) は，
H(z) =
N
i=0
で与えられる．
*1
工学の世界では"デジタル"ではなく"ディジタル"と書くのが普通．
–i–
ai z −1
(4)
フィルタの周波数応答と線形位相特性
上式で表されるフィルタの周波数特性は，
xk = ejkωT
(5)
という周期入力に対するフィルタの応答を考えれば求めることが出来る．ただし，T はサンプリング間隔である．出力
も同じ周波数の周期関数であり，
y(k) = F (ω)ejkωT
(6)
とおけば，
F (ω) =
N
ak e−jkωT
(7)
k=0
となる．言い換えれば，フィルタの周波数応答は，z に対して ejωT
線形位相特性
ところで，フィルタの重み係数が，
ak = aN −k
(8)
という性質を持つ場合，すなわち中心に対して前後対称である場合には，
N = 2M
(9)
とするとき，
F (ω) = aM e−jMωT +
=e
−jMωT
M−1
ak e−jMωT ej(M−k) ωT + e−j(M−k)ωT
k=0
aM + 2
M−1
(10)
ak cos(M − k)ωT
(11)
k=0
である．括弧の中は実数なので振幅を表している．このフィルタによる位相遅れは M ωT であり，位相遅れが周波数
に比例することがわかる．このような特性を線形位相特性と呼ぶ．
フィルタを通すことによって，時間軸で波形が遅れる量は，(周波数 x 位相遅れ) であるから，線形位相特性を持つ
フィルタでは，全ての周波数で波形は同一時間遅れる．いいかえれば，フィルタを通すことによって，通過帯域の周波
数成分によって形成される波形の姿は変わらないことが分かる．これは，波形から現象を解析するに当たって，非常に
有用な性質である．以後は，線形位相特性を与える対称フィルタについてのみ考察する．
– ii –
フィルタの遅れ
ここで考えている長さ 2M + 1 のフィルタは，k = M + 1 に関して前後対称である．このようなフィルタで畳み込
み演算を行うと，出力は M サンプル分入力よりも遅れて現れることに注意しなければならない．実時間の流れに影響
されないオフライン処理であれば，フィルタの係数を −M から M までと定義すれば，位相遅れのないフィルタが得ら
れる．
周波数応答の周期性
フィルタ出力の絶対値は，複素共役な入力に対しても変わらないことから，フィルタの周波数特性は ω = 0 に対して
対称な偶関数であり，またオイラー面における 1 回転すなわち
2π
T
(12)
1
ωs
=
2π
T
(13)
ωs =
周波数で表せば，
fs =
を周期とする周期関数である．したがってまた，周波数応答のゲイン特性は fs /2 で折り返す偶関数であり，意味のあ
る情報は，
0<ω<
ωs
2
(14)
の区間にのみ存在するということができる．この意味で，
m=
2ω
2f
2f
=
=
ωs
fs
1/T
(15)
を無次元周波数と呼ぶ．
周波数応答の解析
式 11 より，ゲイン特性は，
ωT = πm
(16)
より，
M−1
|F (ω)| = aM + 2
ak cos(M − k)πm
(17)
k=0
で与えられる．また，m = 1 における遅れは，
∠F (ω) = πM
である．
– iii –
(18)
単純移動平均フィルタ
2M + 1 の区間での移動平均は最も簡単なローパスフィルタであり，
ak = 1/2M + 1
(19)
k = 0, · · · , 2M
(20)
で定義される．Figure 1 には，M = 5 および M = 10 の移動平均フィルタのゲイン特性を示す．
Figure 1 単純移動平均フィルタ (M = 5 および 10)
通過領域以外のサイドローブが 10% 以上大きく，簡単ではあるが，フィルタとしての特性はよくないことが分かる．
三角フィルタ
もう少しなだらかな重み関数として，三角波形状の重み，
ak = a2(M+1)−k = k/
k = 1, · · · , M + 1
2M+1
ai
(21)
i=1
(22)
を考える．M = 10 の場合，重み関数は Figure 2 のようになる．このフィルタの周波数特性は Figure 3 のようになり，
単純移動平均よりもかなり改善されることが分かる．
– iv –
Figure 2 三角フィルタの係数
Figure 3 三角フィルタの周波数特性
Hamming フィルタ
もう一つの例として，次によって定義される Hamming 窓関数と呼ばれるものを，移動平均の重みとして用いるフィ
ルタを考える．
ak = a2(M+1)−k = [0.54 + 0.46 cos ((M − k + 1)π/M )] /
2M+1
(23)
i=1
M = 10 の場合，次のような値を取る
Figure 4 Hamming フィルタの係数
Hamming 窓関数は，元来そのままフィルタとして用いるものではなく，次節で述べるように理想フィルタを有限項
–v–
で近似する際に，性質を改良するために用いるものではあるが，Figure 3 のように，そのままでフィルタとして用いて
も，単純移動平均より良い性質を示す．単純移動平均に比べて，サイドローブ特性がかなり改善されて 1% 以下になっ
Figure 5 Hamming フィルタのゲイン特性 (M = 5 および 10)
ている．
フィルタの設計
理想フィルタ
理想的なローパスフィルタとは，あるカットオフ周波数 ωc より低周波は完全に通し，それよりも高周波は全く通さ
ないというゲイン特性を持つものである．すなわち，
|H(ω)| = 1
ω ≤ ωc
(24)
|H(ω)| = 0
ω > ωc
(25)
先に述べたように，フィルタの係数が 0 を中心として −M から M までと考えれば位相遅れを考慮しないでよいので，
H(ω) = 1
ω ≤ ωc
(26)
H(ω) = 0
ω > ωc
(27)
となる．
H(ω) は周期 ω の周期関数であるから，フーリエ級数で表現でき，
H(ω) =
∞
h(n)e−jnωT
n=−∞
– vi –
(28)
遅れのある，すなわち係数が 0 から始まっている有限長フィルタの周波数応答の式 7 と比較すれば，上式におけるフー
リエ係数 h(n) が，無限長フィルタの係数を表していることが分かる．逆フーリエ展開より，
h(n) =
1
ωs
=
1
ωs
ωs /2
H(ω)ejnωT dω
−ωs /2
ωc
jnωT
−ωc
e
dω
(29)
(30)
であり，解は，
h(n) =
ωc T sin(nωc T )
π
nωc T
(31)
で与えられる．ただし，h(0) = ωc T /π である．後の便宜のために，m = 2ω/ωc および mc = 2ωc /ωs を用いて書き直
せば，
h(n) = mc
sin(πnm)
πnmc
(32)
である．
理想フィルタの有限近似
上式で表されるフィルタは n = −∞ ∼ ∞ までの無限個の係数を持つが，これを有限個の係数でうち切って近似する場
合について考える．簡単のために，2 以上の整数 Nc を用いて，
mc =
1
Nc
(33)
であるとすれば，
h(n) =
1 sin(πn/Nc )
Nc πn/Nc
(34)
であり，係数は Nc ごとに 0 となる．Figure 6 は，Nc = 5 の場合の重み係数を示す．なお，有限個でうち切った場合
にも係数の総和が 1 となるように，
M
1 sin(πn/Nc ) h(n) =
/
h(k)
Nc πn/Nc
(35)
k=−M
によって求め，また 1 から始まるようにずらしてある．この場合，M = Nc ，M = 2Nc ，M = 3Nc ，M = 4Nc でそ
れぞれうち切った場合の周波数応答は，次のようになる． 係数の個数を増すほど，カットオフはシャープになるが，同
時にいわゆるギブス現象のために，カットオフ直前でリンギングが発生している．
窓関数によるギブス現象の抑制
ギブス現象は，重み関数に窓関数と呼ばれる特殊な重み関数を乗ずることによって抑制できることが知られている．
簡単でよく用いられるのが，Hamming 窓関数であり，遅れなしの形では
ak = a−k = 0.54 + 0.46 cos π
– vii –
k
M
(36)
Figure 6 Nc = 5, M = 20 の場合の重み係数の分布
で与えられる．M = Nc , 2Nc , 3Nc , 4Nc でうち切った有限フィルタに Hamming 窓関数を乗じた結果の周波数応答を，
Figure 8 に示す．なお，式 36 を乗じた結果の係数は，総和が 1 になるように定数倍してある． いずれの場合も，ギブ
ス現象がほとんどなくなっているだけでなく，サイドロープも減少して，大いに改良されたフィルタ特性が得られてい
る．なお，M = 2Nc の場合の重み関数は，次のような姿をしている．単に無限フィルタをうち切っただけの Figure 6
と比較すると，Hamming 窓の効果がよく分かる．
実用的なフィルタ
サンプリング間隔に対するカットオフ周波数の関係，それとフィルタ長さの関係はどのようであれば，実用的な性能
のフィルタとして用いることができるだろうか．サンプリング間隔とフィルタのカットオフ周波数の関係は無次元カッ
– viii –
Figure 7 M = Nc , 2Nc , 3Nc , 4Nc でうち切った理想フィルタ
トオフ周波数
mc =
2ωc
2fc
=
= 2fc T
ωs
fs
(37)
1
mc
(38)
または，
Nc =
で表される．たとえば，1 秒サンプリングの場合，カットオフ周波数が 0.1Hz で無次元カットオフ周波数は 0.2，Nc は
5 になる．Nc が整数の場合，フィルタ係数は Nc ごとに 0 となる性質をもつので，フィルタ長は中心から片側の長さが
Nc の何倍あるかで評価できる．片側フィルタ長が kNc のとき，最後の 0 は不要なので，フィルタの全長は 2kNc − 1
である．以下に，2, 4, 8 の Nc について，k が 2,4 の場合についての周波数応答を示す．
これらから，次のようなことが分かる．
– ix –
Figure 8 Hamming 処理を施した有限フィルタ (M = Nc , 2Nc , 3Nc , 4Nc )
Figure 9 M=2Nc の場合の重み関数
• Nc でカットオフ周波数が変化するが，他の特性には目立った影響を与えない．ただし，が 2 では，最初のサイ
ドローブに達する前に周波数応答の 1 周期に到達している．4 程度以上が望ましいといえるだろう．
• k はカットオフのシャープネスを定める．2 では周波数が 4 倍ほど変化する間にゲインは 1% 程度に落ちて
いる (-20dB/Oct) が，4 では 1.5 倍くらいの間に 1% 程度になっており，非常にするどい．二次のフィルタで
も-12dB/Oct であることを考えれば，-20dB/Oct でも実用的には十分であろうが，いかにもディジタルフィルタ
というような鋭い切れ味を実現したければ， に選べば良いであろう．
–x–
– xi –
• すなわち，サンプリング周期とカットオフ周波数は，
Nc =
1
≥4
2T fc
(39)
となるように選ぶ．1 秒サンプリングならば，カットオフ周波数は 1/8Hz 以下．このとき，フィルタの片幅を少
なくとも 2Nc ，余裕があれば 4Nc に選ぶ．
付録：フィルタ係数の例
Hamming 処理を行った有限フィルタの係数の例と周波数応答を示す．カットオフ周波数は T = 1 秒について示す．
1 --- -0.00078996
2 --- -0.00145500
3 --- -0.00187362
4 --- -0.00158819
5 --- 0.00000000
6 --- 0.00306148
7 --- 0.00673462
8 --- 0.00898611
9 --- 0.00728419
10 --- -0.00000000
11 --- -0.01204878
12 --- -0.02477581
13 --- -0.03143140
14 --- -0.02475258
15 --- 0.00000000
16 --- 0.04235435
17 --- 0.09557911
18 --- 0.14800081
19 --- 0.18643814
20 --- 0.20055302
21 --- 0.18643814
22 --- 0.14800081
23 --- 0.09557911
24 --- 0.04235435
25 --- 0.00000000
26 --- -0.02475258
27 --- -0.03143140
28 --- -0.02477581
29 --- -0.01204878
– xii –
30 --- -0.00000000
31 --- 0.00728419
32 --- 0.00898611
33 --- 0.00673462
34 --- 0.00306148
35 --- 0.00000000
36 --- -0.00158819
37 --- -0.00187362
38 --- -0.00145500
39 --- -0.00078996
Bibliography
1. Bozic, S. M., Digital and Kalman Filtering, Edward Arnold Ltd, 1979.
2. 中村尚五, ビギナーズ・ディジタルフィルタ, 東京電機大学出版局, 1989.
– xiii –