FIRフィルタを計算してみる

by K.I

2015/12/15〜

Index

概要
FIRフィルタ
- FIRフィルタの式
- LPFの設計値
- Maximaで計算
- ゲインと位相のグラフ
- FIRフィルタの計算（改良版）
- リストでFIR計算（170215追記）
- FIR係数のプロット（170215追記）
IIRフィルタ
- IIRフィルタの式
- IIRのLPF設計値
- IIRのLPF特性
- IIRのLPF特性グラフ
- IIRのBPF特性
- IIRのBPF特性グラフ（ちょっと変かも）
まとめ
- 位相の変化
- 参考

概要

Maximaで、以前、いろいろなデジタルフィルタを計算してみたけど、

今度は、FIRフィルタの計算をしてみようと思う

フィルタで、位相を変えずにゲインを変化させるのは出来ないって思ってたんだけど、

デジタルフィルタでは、インパルス応答に対して対称の係数のフィルタにすれば良いらしい。

FIRフィルタって、そういうフィルタだから位相は変わらないのかなぁ

FIRフィルタを使ったことなかったんだけど、ちょっと気になる

でもFIRフィルタの設計は良くわからないので、

石川高専の DF-Designのお世話になりました。。

構成的には、移動平均フィルタと同じような感じだから、

やっぱり、数値計算ソフトの Maximaを使って、計算してみよう。

ついでに、IIRフィルタの計算もしてみたんだけど、

これは残念ながら、うまく行かなかったっぽい。。

C言語でFIRフィルタのプログラムを作って、サンプルレート変換を試してみたので、

そちらも参照してみてください → FIRフィルタでサンプルレート変換（170131追記）

ここで計算したFIRフィルタをFPGAで実際に動かしてみた →FPGAでFIRフィルタ（170615追記）

[top]

FIRフィルタ

FIRフィルタの構成

サンプリングしたデータに係数を掛けて、加算するだけ。

移動平均も、この一種ということになる。

FIRフィルタの式

FIRフィルタの差分方程式は、

z変換で、伝達関数を求める

LPFの設計値

FIRフィルタで、LPFを作ってみる。

といっても設計はよく分からないので、こちらにお任せ

----- From here -----
13
  4.036409219416838e-03 
 -2.208366021329011e-18 
 -2.346162858786038e-02 
 -3.367761421639100e-02 
  7.203267485172524e-02 
  2.840739092673040e-01 
  4.000000000000000e-01 
  2.840739092673040e-01 
  7.203267485172524e-02 
 -3.367761421639100e-02 
 -2.346162858786038e-02 
 -2.208366021329011e-18 
  4.036409219416838e-03 

Filter type                      : LPF
Window type                 : Hamming
Filter length                   : 13
Normalized cut-off freq  : 0.200000
Date and time                 : 2015/12/15 18:47
----- End -----

窓関数法ってやつで、フィルタ長は長すぎると面倒なので、N=13に設定。

Maximaで計算

Maximaで、FIRフィルタのN=13の場合の式を作る

kill(all)$
a0:   4.036409219416838e-03 ;
a1:  -2.208366021329011e-18 ;
a2:  -2.346162858786038e-02 ;
a3:  -3.367761421639100e-02 ;
a4:   7.203267485172524e-02 ;
a5:   2.840739092673040e-01 ;
a6:   4.000000000000000e-01 ;
a7:   2.840739092673040e-01 ;
a8:   7.203267485172524e-02 ;
a9:  -3.367761421639100e-02 ;
a10: -2.346162858786038e-02 ;
a11: -2.208366021329011e-18 ;
a12:  4.036409219416838e-03 ;
a0+a1*exp(-j*omega*T)+a2*exp(-2*j*omega*T)+a3*exp(-3*j*omega*T)+a4*exp(-4*j*omega*T)+a5*exp(-5*j*omega*T)+a6*exp(-6*j*omega*T)+a7*exp(-7*j*omega*T)+a8*exp(-8*j*omega*T)+a9*exp(-9*j*omega*T)+a10*exp(-10*j*omega*T)+a11*exp(-11*j*omega*T)+a12*exp(-12*j*omega*T);

続けて、サンプリング周波数1kHz（T=0.001）で計算

subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
H:%;
plot2d(cabs(H),[f,0,1000],[logy]);
plot2d(carg(H),[f,0,1000]);

ゲインと位相のグラフ

ゲイン特性

十分、急峻なLPFの特性が得られている

位相特性

直線的に変化しているので、多分、遅延はあるけど、歪は少ないってことなのかな

まぁ、少なくともちゃんと計算出来ているんじゃないかな。

FIRフィルタの計算（改良版）

FIRの係数を配列に入れたら、伝達関数ももっとスマートに計算できるかもしれない。

まず、差分方程式は

それで、伝達関数は、

Maximaで伝達関数を記述すると、

sum(a[k]*exp(-k*j*\omega*T),k,0,N-1);

あらためて、Maximaで計算してみる

kill(all)$
array(a,13);
a[0]:   4.036409219416838e-03 ;
a[1]:  -2.208366021329011e-18 ;
a[2]:  -2.346162858786038e-02 ;
a[3]:  -3.367761421639100e-02 ;
a[4]:   7.203267485172524e-02 ;
a[5]:   2.840739092673040e-01 ;
a[6]:   4.000000000000000e-01 ;
a[7]:   2.840739092673040e-01 ;
a[8]:   7.203267485172524e-02 ;
a[9]:  -3.367761421639100e-02 ;
a[10]: -2.346162858786038e-02 ;
a[11]: -2.208366021329011e-18 ;
a[12]:  4.036409219416838e-03 ;
sum(a[k]*exp(-k*j*\omega*T),k,0,N-1);
subst(13,N,%);
subst(%i,j,%);
subst(2*%pi*f,\omega,%);
subst(0.001,T,%);
H:%;
plot2d(cabs(H),[f,0,1000],[logy]);
plot2d(carg(H),[f,0,1000]);

結果は同じになる。このほうが、係数とフィルタ長だけ指定するだけなので簡単。

リストでFIR計算（170215追記）

配列じゃなくて、リストで係数を定義するようにした

結果は変わらない

kill(all)$
a: [
4.036409219416838e-03,
-2.208366021329011e-18,
-2.346162858786038e-02,
-3.367761421639100e-02,
7.203267485172524e-02,
2.840739092673040e-01,
4.000000000000000e-01,
2.840739092673040e-01,
7.203267485172524e-02,
-3.367761421639100e-02,
-2.346162858786038e-02,
-2.208366021329011e-18,
4.036409219416838e-03];

sum(a[k+1]*exp(-k*j*\omega*T),k,0,12);
subst(%i,j,%);
subst(2*%pi*f,\omega,%);
subst(0.001,T,%);
H:%;
wxplot2d(cabs(H),[f,0,1000],[logy]);
wxplot2d(carg(H),[f,0,1000]);

あとplot2dの代わりに、wxplot2dを使う¹ようにした。

これだと画像の保存が右クリックで出来るので便利

訂正：リストは配列と同様に使えるんだけど、添字が１から始まること、あとNのような変数を使えないのでsumの式を修正

配列のプログラムを実行後、そのまま実行したため間違いに気がつかなかった。そのため kill(all)も追加した。（170615）

FIR係数のプロット（170215追記）

係数をリストにしたので、係数のプロットが簡単²になった。

x: makelist(n,n,0,12);
wxplot2d([discrete,x,a],[style,impulses]);

pointsで表示した方が良いかも

x: makelist(n,n,0,12);
wxplot2d([discrete,x,a],[style,[points,2,1,1]]);

¹Windows用の、wxMaximaじゃないと使えないけれども。
²ということで、以前無理やりプロットしたんだけど、それは消した。

[top]

IIRフィルタ

同様に、IIRフィルタも計算してみた

でも、こちらはちょっと怪しい。。

DF-DesignのIIRフィルタ（１セクション）の構成は、こんなふうになっている。

IIRフィルタの式

このIIRフィルタ（１セクション）の差分方程式は、この構成からはよく分からないけど、

NIの IIRフィルタとFIRフィルタの解説によると、IIRフィルタの差分方程式は、こんな式になるらしい。

それに、構成の似ているノッチフィルタの差分方程式から類推すると、

たぶん差分方程式は、以下のような式になるんじゃないかと思う³。

z変換で書き換えて、伝達関数を求める

IIRのLPF設計値

LPFは、３セクション構成になるらしい。

やはり、 DF-Designで計算

----- From here -----
3 

2.23801748954850749e-01
1.04793004180597144e-01
0.00000000000000000e+00
2.46824032729850495e+00
2.46824032729850495e+00
0.00000000000000000e+00

5.53109465901169481e-02
2.31444265788747300e-01
1.16419411230618330e-01
6.09221031330263063e+00
1.21844206266052613e+01
6.09221031330263063e+00

7.60854166202676457e-02
3.18373386689702720e-01
5.35740052613985052e-01
6.09221031330263063e+00
1.21844206266052613e+01
6.09221031330263063e+00
Filter type                : LPF(IIR Butterworth)
Order of transfer function : 5
Number of sections         : 3
Sampling freq.      [Hz]   : 1000.000000
Passband edge freq. [Hz]   : 250.000000
Stopband edge freq. [Hz]   : 400.000000
Passband attenuation[dB]   : 0.500000
Stopband attenuation[dB]   : 35.000000
----- End -----

IIRのLPF特性

Maximaを使って、数値計算させてみた

サンプリング周期Tを1kHzとして、特性を求める

直列接続なので、単純に伝達関数を掛け算してみたが、これで良いのかな？

kill(all)$
(a0+a1*exp(-j*omega*T)+a2*exp(-2*j*omega*T))/(1/k+b1*exp(-j*omega*T)+b2*exp(-2*j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
define(H(f,k,b1,b2,a0,a1,a2),%);
H1(f):=H(f,
2.23801748954850749e-01,
1.04793004180597144e-01,
0.00000000000000000e+00,
2.46824032729850495e+00,
2.46824032729850495e+00,
0.00000000000000000e+00);
H2(f):=H(f,
5.53109465901169481e-02,
2.31444265788747300e-01,
1.16419411230618330e-01,
6.09221031330263063e+00,
1.21844206266052613e+01,
6.09221031330263063e+00);
H3(f):=H(f,
7.60854166202676457e-02,
3.18373386689702720e-01,
5.35740052613985052e-01,
6.09221031330263063e+00,
1.21844206266052613e+01,
6.09221031330263063e+00);
plot2d(cabs(H1(f)*H2(f)*H3(f)),[f,0,499],[logy]);
plot2d(carg(H1(f)*H2(f)*H3(f)),[f,0,499]);

Maximaの記述方法が、なんか滅茶苦茶だけど。。

最初、プロット範囲を0〜1000Hzにしたら、エラーになった。

Maxima encountered a Lisp error:
 Error in PROGN [or a callee]: Zero is the logarithmic singularity.
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.

どうやらナイキスト周波数付近で計算が発散するっぽいので、0〜499Hzにしたら出来た。

IIRのLPF特性グラフ

ゲイン特性

IIRのLPF特性は単純に落ちるだけで、FIRのLPFの特性とはかなり違う。

位相特性

IIRのLPFも、周波数特性は直線的なので、歪は少ないようだ。

IIRのBPF特性

IIRのBPFも同様に設計値を当てはめて見る

これは４セクション構成になる。（サンプリング周波数は1kHz）

kill(all)$
(a0+a1*exp(-j*omega*T)+a2*exp(-2*j*omega*T))/(1/k+b1*exp(-j*omega*T)+b2*exp(-2*j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
define(H(f,k,b1,b2,a0,a1,a2),%);
H1(f):=H(f,
9.115605050469990e-02,
 -2.503362204683822e-01,
 4.504947873400970e-01,
 1.000000000000000e+00,
 0.000000000000000e+00,
 -1.000000000000000e+00);
H2(f):=H(f,
 1.000000000000000e+00,
 2.503362204683819e-01,
 4.504947873400968e-01,
 1.000000000000000e+00,
 0.000000000000000e+00,
 -1.000000000000000e+00);
H3(f):=H(f,
 1.189175782425170e-01,
 -6.445233834458046e-01,
 7.545210372320261e-01,
 1.000000000000000e+00,
 0.000000000000000e+00,
 -1.000000000000000e+00);
H4(f):=H(f,
 1.000000000000000e+00,
 6.445233834458044e-01,
 7.545210372320261e-01,
 1.000000000000000e+00,
 0.000000000000000e+00,
 -1.000000000000000e+00);
plot2d(cabs(H1(f)*H2(f)*H3(f)*H4(f)),[f,1,499],[logy]);
plot2d(carg(H1(f)*H2(f)*H3(f)*H4(f)),[f,1,499]);