Butterworthフィルタ
Simple Wiki Based Contents Management System
ソフトウェア関連 >> ツール >> Butterworthフィルタ

Butterworthフィルタ

IMUで姿勢推定やオドメトリの推定をしたいと思ってプログラミングしているとどうしてもセンサデータのばらつきなどを解消しようとデジタルフィルタを使う必要が出てきた。
Butterworthフィルタというのが定番のようなのだが、C/C++での実装があまり見当たらないため自分で作ることにした。
Butterworthフィルタの実装のために必要な式など備忘録として書き留めておく

Butterworth Low passフィルタ

伝達関数 \(H(s)\) は、下記のバターワース多項式の逆数である。(\(m\)は整数)

B ( s ) = k = 1 n / 2 [ s 2 2 s cos ( 2 k + n 1 2 n π ) + 1 ] (for \(n=2m\))

B ( s ) = ( s + 1 ) k = 1 ( n 1 ) / 2 [ s 2 2 s cos ( 2 k + n 1 2 n π ) + 1 ] (for \(n=2m+1\))

また、バターワース多項式の係数配列を下記の式で整理しておく。

B ( s ) = k = 0 n c k s k

次に、\(s\)にプリワーピングを考慮した z 1 の多項式

s = ( 2 T 1 z 1 1 + z 1 ) / ( 2 T tan ( 2 π f c T / 2 ) )
を代入してフィルタ係数\(a, b\)を導出する。

ここで、 y = z 1 w = tan ( 2 π f c T / 2 ) = tan ( π 2 f c F ) とすると

s = 1 y w ( 1 + y )

と書くことができる。

以上のことから、伝達関数\(H(s)\)は、

H ( s ) = 1 B  (  y   ) = 1 k = 0 n c k s k = 1 k = 0 n c k / w k ( 1 y ) k ( 1 + y ) k

となり、この式の部分分数を処理して、

H  (  y   ) = ( 1 + y ) n k = 0 n c k / w k ( 1 y ) k ( 1 + y ) n k

となる。さらに、 a 0 = 1 とするために、 A 0 = k = 0 n c k / w k として正規化すると、

H  (  y   ) = A 0 1 ( 1 + y ) n A 0 1 k = 0 n c k w k ( 1 y ) k ( 1 + y ) n k

となる。 したがって、Butterworth LowPassフィルタのフィルタ係数は、

b ( y ) = ( 1 + y ) n / k = 0 n ( c k / w k )

a ( y ) = k = 0 n c k w k ( 1 y ) k ( 1 + y ) n k / k = 0 n ( c k / w k ) となる。

フィルタ係数 \(a,b\)が決まれば、下式を使ってフィルタリングを行う。

y ( t ) = k = 0 n b k x ( t k ) k = 1 n a k y ( t k )

ここで注意するのは、上式でフィルタリングすると位相のずれが生じる。そのため、オフラインでフィルタリングを行う場合には、セロ位相フィルタリングを行った方が良い。
ゼロ位相フィルタリングは、上のフィルタリングの結果を反転させて、再度フィルタリング処理を行い、その結果の反転をゼロ位相フィルタリングの結果とするというものである。
また、フィルタリングを行う際は、データの始点と終点を補完しておいた方がよい。

Butterworth High passフィルタ

Butterworth High passフィルタのフィルタ係数を求めるには、Lowpassフィルタでの双一次変換のときの\(s\)を
s = w ( 1 + y ) 1 y

として処理すればよい。
したがって、伝達関数\(H(s)\)は、

H ( s ) = 1 B = 1 k = 0 n c k s k = 1 k = 0 n c k w k ( 1 + y ) k ( 1 y ) k

となる。この式もLow Passフィルタの時と同様に部分分数処理と A 0 = k = 0 n c k w k として正規化を行うと

H ( y ) = A 0 1 ( 1 y ) n A 0 1 k = 0 n c k w k ( 1 + y ) k ( 1 y ) n k

となり、フィルタ係数\(a,b\)は、

b ( y ) = ( 1 y ) n / k = 0 n ( c k w k )

a ( y ) = k = 0 n ( c k w k ) ( 1 + y ) k ( 1 y ) n k / k = 0 n ( c k w k )

となる。