デジタルフィルターの原理と応用

デジタルフィルターの原理と応用 (wrote 2006 12.17 2015.5.6)

このページではデジタルフィルターについて説明します。 原理から実用的なコードまで短く説明していきたいと思います。 これを書こうと思ったきっかけは、 なぜかここに書いてあるような原理と応用が書いてあるページが 見当たらなかったからです。

デジタルフィルターは少し前まではＤＳＰで処理していましたが、 ２０００年頃から徐々にＣＰＵでも処理できるようになってきました。 気軽にプログラムが組める環境になっているのにその原理は、 どうもあまり知られていないようです。

１次フィルターの原理と応用

１次のローパスフィルター（ＬＰＦ）は次のようにプログラムで組むことができます。

void filter_LPF1( f ) float f; { float in,out,a; a=1.0 - exp(-2.0*M_PI*f*dt); out=0; for (;;) { in=get_input_value(); out+=(in - out)*a; put_output_value( out ); } }

f は出力が -3dB になる周波数です。 a はフィルターの時定数をあらわします。a の範囲は０以上１以下です。 仮に a=1 であるとすると

out+=(in - out)*a; (式１－１)

と言う式は

out=in;

となり入力をそのまま出力することになります。 a=0 であれば out=out; となり 0 を出力し続けることになります。 a の値が中間値であるときは、out の値は in の値を追いかけて変化します。 仮に初期値 a=0.5; out=0; in=1; であり in の値は一定だとすると out の値はループを回る毎に 0.5, 0.75, 0.875,... と言うように in の値に 漸近していきます。このような原理でＬＰＦとしての動作をします。

次に、定数 a の算出法について説明します。

まず、電気回路で抵抗とコンデンサーを用いて １次のＬＰＦを構成した場合と比較してみます。回路は

[in]---[抵抗]---[out]---[コンデンサー]---アース

と言うようになります。 この回路の動作を数式で表すと次のようになります。

d(out)/dt=(in - out)*w (式１－２)

フィルターの周波数特性を表す値(角速度w)は w=1/(R*C) で、 その周波数 f は f=2π*w です。この周波数でＬＰＦの出力は -3dB になります。 （この辺りの電気回路の話は多数のページで説明されています。 google で 「RC回路」で検索すると見つけることができます。）

式１－２は式１－１に似ていることが解ります。 そのことから、a と w の関係を考えるといいのではないかと予測がつきます。

まず、in が一定である場合を考えます。式１－２を積分すると

out(t)=in*(1-exp(-w*t)) (式１－３)

となります。

一方、プログラムの一回のループで dt の時間がたつと決めると、 ループのある回の出力とその次の出力は out(t) と out(t+dt) です。 dtは具体的に 0.01 秒とかの数値が入ります。式１－１はプログラムの式ですが、 これを計算式に書き換えると

out(t+dt)=out(t) + (in - out(t))*a (式１－４)

です。式１－４に式１－３を代入すると、

a=1 - exp(-w*dt) (式１－５)

このように a を w と dt で表すことができます。 πを M_PI と書くと w=2.0*M_PI*f なので、

a=1.0 - exp(-2.0*M_PI*f*dt) (式１－５．１)

というように具体的な a の値を得ることができます。 この計算では in を一定としましたが、実際は変化します。 このように条件が異なるので、それが a の有効範囲として影響してきます。 f*dt が１に近付くにつれて誤差が大きくなり、 実際の -3dB点は与えた f よりも低くなってきます。 f*dt<0.05 では誤差は５％以下のようです。 実際に使うときはこの式でおおよその数値を出しておき、 実測して補正するといいと思います。

in の条件を変えて計算すると、さらによくあう解が見つかるかもしれません。

次に、d(out)/dt=lim((out(t+dt) - out(t))/dt) であることを考慮すると、 式１－２と式１－４より

a=2.0*M_PI*f*dt (式１－６）

という関係をえることができます。この式を実際に使うと、 こんどは逆に-3dB点は与えた f よりも高くなります。 そのため、答えは式１－５．１と式１－６の間にあることが予測できます。

a=(2.0*M_PI*f*dt)*0.667 + (1.0 - exp(-2.0*M_PI*f*dt))*0.333 (式１－７）

中間値をとりこのような式が比較的あうようです。

このaの値に相当する値を得ることが出来るソフトウエアがありますがそれらは内部でデジタルフィルターをシミュレートしてフィードバックをかけながら自動で最適値を計算しているのだとおもいます。(wrote2015.5.6)

１次のハイパスフィルター（ＨＰＦ）は ＬＰＦで排除した信号を通過させ、通過させた信号を排除すればいいので、 入力からＬＰＦの出力を引けばいいことが解ります。 そのため、上記のプログラムを多少改造することで得ることができます。

それを含めて実用的な形に書き換えると

typedef struct { enum _filter1_mode { LPF,HPF } mode; float a,val1; } filter1Rec,*filter1; static inline void init_filter1( fi,mode,dt,f ) filter1 fi; enum _filter1_mode mode; float dt,f; { fi->mode=mode; fi->a=1.0 - exp(-2.0*M_PI*f*dt); fi->val1=0; } static inline float calc_filter1( val,fi ) float val; filter1 fi; { float v; fi->val1+=(val - fi->val1) * fi->a; if ( fi->mode==LPF ) v=fi->val1; if ( fi->mode==HPF ) v=val - fi->val1; return( v ); } main() { filter1Rec fi; float v; init_filter1( &fi,LPF,0.001,10.0 ); for (;;) { v=get_input_value(); v=calc_filter1( v,&fi ); put_output_value( v ); } }

このような形になります。プログラムコード上では calc_filter1() の中で 毎回フィルターモードの条件判断をしている形になっていますが、 static inline の関数にしておくことでコンパイラが最適化して 条件判断のコードは消えます。 しかし、コールバックルーチンの中で使う場合は init_filter1() と calc_filter1() が別々の関数から呼ばれることになり、 最適化できなくなる可能性があるので注意が必要です。

２次フィルターの応用

２次フィルターは１次フィルターを２段重ねることで実現できます。 しかしその方法ではＱの値をコントロールできません。 これを、バイカッド型のフィルター(biquad filter)を用いることで実現できます。

aの計算方法は１次も２次も同じですが最適化は目的のプログラムを用いて微調整する必要があります。(wote2015.5.6)

void filter_LPF2( f,q ) float f,q; { float in,val0,val1,val2,a; a=1.0 - exp(-2.0*M_PI*f*dt); val1=0; val2=0; for (;;) { in=get_input_value(); val0=in - val2 - val1/q; val1+=val0 * a; val2+=val1 * a; put_output_value( val2 ); } }

このプログラムでは式１－１に替わって

val1+=val0 * a; (式２－１)

という式を使います。この式は積分器になります。 積分なので正の値を入れ続けると値は増加していってしまいます。 このプログラムでは val0 を通してフィードバックがかかっているので、 そのようなことは起こりにくいのですが、q を大きくすると 中間に位置する val1 の値が大きく振幅することがあります。 ここでは浮動小数を用いているので多少の振動には大丈夫なのですが、 固定小数を用いる場合には注意する必要があります。

a の値は、val1 val2 というように１次のフィルターが２個つらなっている形に なっているので１次のフィルターと同じ方法で計算できます。

この２次のフィルターを１次のフィルターも含めて実用的に書き換えると 次のようになります。

typedef struct { enum _filter_mode { LPF1,HPF1,LPF2,HPF2,BPF2 } mode; float a,rq,val1,val2; } filterRec,*filter; static inline void init_filter( fi,mode,dt,f,q ) filter fi; enum _filter_mode mode; float dt,f,q; { fi->mode=mode; fi->a=1.0 - exp(-2.0*M_PI*f*dt); fi->rq=1.0/q; fi->val1=0; fi->val2=0; } static inline float calc_filter( val,fi ) float val; filter fi; { float v; if ( fi->mode==LPF1 || fi->mode==HPF1 ) { v=val - fi->val1; fi->val1+=v * fi->a; if ( fi->mode==LPF1 ) v=fi->val1; if ( fi->mode==HPF1 ) v=val - fi->val1; } else { v=val - fi->val2 - fi->val1*fi->rq; fi->val1+=v * fi->a; fi->val2+=fi->val1 * fi->a; if ( fi->mode==LPF2 ) v=fi->val2; if ( fi->mode==HPF2 ) v=val - fi->val2 - fi->val1*fi->rq; if ( fi->mode==BPF2 ) v=fi->val1 * fi->rq; } return( v ); }

main() { filterRec fi; float dt; float v; dt=1.0/48000.0; init_filter( &fi,LPF2,dt,1000.0,0.71 ); for (;;) { v=get_input_value(); v=calc_filter( v,&fi ); put_output_value( v ); } }

init_filter() でモードを設定することで動作を切り替えることができます。 mode に与える値は次のものがあります。

LPF1 １次ローパスフィルタ HPF1 １次ハイパスフィルタ LPF2 ２次ローパスフィルタ HPF2 ２次ハイパスフィルタ BPF2 ２次バンドパスフィルタ

calc_filter() の次の２行を見比べてみるとわかるように

if ( fi->mode==HPF1 ) v=val - fi->val1; if ( fi->mode==HPF2 ) v=val - fi->val2 - fi->val1*fi->rq;

２次ハイパスフィルターでは入力信号からローパス信号を引きますが、 さらにバンドパス信号も引く必要があります。

ｎ次フィルターの応用

電気回路では２次をこえるｎ次のフィルターは ２次と１次のフィルターを組み合わせて作ります。 デジタルフィルターでも同様のことができます。 複数組み合わせる場合は、f と q を調整して、 肩特性をコントロールすることができるようになります。

またさらにデジタルフィルターではFIR、IIRフィルターを用いることでプログラム的には簡単に実現できます。しかし、接続やパラメーターの決定方法が難しくなります。z変換の計算上では全く同じ特性の２回路同士であってもz-1の値をより多く足すとぼやけた音になりがちですので、そうならない接続方法を探す必要があります。(wote2015.5.6)

main() { filterRec fi1,fi2,fi3,fi4; float f,dt; float v; dt=1.0/48000.0; f=1000.0; init_filter( &fi1,LPF2,dt,f*0.382,0.593 ); init_filter( &fi2,LPF2,dt,f*0.645,1.183 ); init_filter( &fi3,LPF2,dt,f*0.894,2.453 ); init_filter( &fi4,LPF2,dt,f*1.0342,8.0818 ); for (;;) { v=get_input_value(); v=calc_filter( v,&fi1 ); v=calc_filter( v,&fi2 ); v=calc_filter( v,&fi3 ); v=calc_filter( v,&fi4 ); put_output_value( v ); } }

filterRec と init_filter() と calc_filter() は ２次のフィルターのプログラムと同等です。

このプログラムは「トランジスター技術ＳＰＥＣＩＡＬ No.44」の２６ページの 例題から定数を持ってきた、 リプルが０．１ｄＢの８次チェビシェフ型フィルターです。

肩特性の性質とそれぞれの２次のフィルターに与えるfとqの定数の関係は google で バタワース特性、ベッセル特性、チェビシェフ特性などのキーワードで 検索すると知ることができます。

ΔΣ変換の原理と応用

int の値を出力できるＤＡコンバーターがあった場合、 float で 10.25 の値を連続して出力しようとしたとき、 四捨五入して 10,10,10,10,... と出力することが一般的だと思います。 しかし、10,10,10,11,10,10,10,11,... というように４回に１回の割合で11を出力してやると平均して10.25の数値を 出力したのと同じようなことになります。 ＤＡコンバーターの出力をＬＰＦを通すことで 10.25 に対応する電圧を得ることができます。

出力データーの列のうち１／４の割合で11が混じっていればよいので、 必ずしも、４回に１回11を出力する必要はなく、 10を３０回、11を１０回出力することを繰り返してもよいことになります。 しかし、そうした場合は特性周波数を１／１０にしたＬＰＦを通さないと 前と同じような誤差をもつ値をえることができません。 何らかの方法で分散させる必要がありますが、特別な方法を使わなくとも 次のようなプログラムでうまい具合に分散します。

void filter_delta_sigma( void ) { float in,out,val; val=0; for (;;) { in=get_input_value(); val+=in; out=floor(val - 0.5); val-=out; put_output_value( out ); } }

このプログラムは出力は float ですが floor() で小数以下を切り捨てているので 多ビットの int とみなすことができます。 こうすることで、int 対応のＤＡコンバーターに出力できるようになります。

また、表現を変えると、入力のビット数を減らしてその情報を時間軸に分散させる プログラムであるといえます。float の値の範囲をあわせることで、 例えば２４bitsの信号を１６bitsに変換することができます。

このプログラムでは、10 と 11 は周期的に出力されます。 また、入力に例えば音楽信号を用いた場合では、 音楽のリズムにあわせてノイズが発生します。そこで、

val+=in + noise();

というように元の信号と無関係のノイズを加えると、 周期性は消えて、リズムにあわせたノイズ成分は減少します。

次に、

out=floor(val - 0.5);

の部分を

if ( val>0 ) out=1.0; else out=-1.0;

に置き換えると１ビットの出力にすることができます。これをアレンジして用いて、 オーディオに使われている ６４倍オーバーサンプリングの１ｂｉｔΔΣ変換フィルター を書くと次のようになります。

void filter_1bit_delta_sigma( void ) { float in,out,val,v; filterRec fi_in,fi; int i,n; n=64; dt=1.0/44100.0/(float)n; init_filter( &fi_in,LPF2,dt,44100.0,0.71 ); init_filter( &fi,HPF2,dt,44100.0,0.71 ); val=0; for (;;) { in=get_input_value(); for ( i=0;i<n;i++ ) { val+=calc_filter( &fi_in,in ); if ( calc_filter(&fi,val)>0 ) out=1.0; else out=-1.0; val-=out; put_output_value( out ); } } }

このプログラムのうち、

val+=calc_filter( &fi_in,in );

の calc_filter() の部分がオーバーサンプリング用のフィルターです。

また、１ビット化のための比較部分が

if ( calc_filter(&fi,val)>0 ) out=1.0; else out=-1.0;

このように替わっています。 このcalc_filter()は出力信号に含まれるノイズ成分を高域側に シフトするためのものです。 このフィルターを使わない場合では、ノイズスペクトルは周波数が 上がるにつれて多くなりますが、 これを使うとこのフィルターで制限される領域ではフラットになります。 (しかし、高域のノイズ成分は増加します）

ちなみに、１ビット－マルチビット混合のＤＡコンバーターは

out=floor(calc_filter(&fi,val)/32 + 0.5);

というように元に戻し、桁をあわせることで実現できます。

また、低周波側はマルチビット、 高周波側は１ビットという組み合わせもいいのかもしれません。 クロスオーバー周波数は２０Ｈｚで１次のフィルターを用いると、 特定の周波数にくせが表れることもないと思います。

式からも解るように２次以上のフィルターを用いると合成しても 元の信号に戻らないので必ず１次でないといけません。