FIRフィルタを計算してみる
by K.I
2015/12/15～
Index
- フィルタで、位相を変えずにゲインを変化させるのは出来ないって思ってたんだけど、
- デジタルフィルタでは、インパルス応答に対して対称の係数のフィルタにすれば良いらしい。
- FIRフィルタって、そういうフィルタだから位相は変わらないのかなぁ
- FIRフィルタを使ったことなかったんだけど、ちょっと気になる
- ついでに、IIRフィルタの計算もしてみたんだけど、
- C言語でFIRフィルタのプログラムを作って、サンプルレート変換を試してみたので、
[top]
- サンプリングしたデータに係数を掛けて、加算するだけ。
\[
y[n] = a_0x[n] + a_1x[n-1] + a_2x[n-2] + a_3x[n-3] + a_4x[n-4] + \ \cdots \ + a_{N-1}x[n-(N-1)]
\]
\[
Y(z) = a_0 X(z) + a_1 z^{-1} X(z) + a_2 z^{-2} X(z) + a_3 z^{-3} X(z) + a_4 z^{-4} X(z) + \ \cdots \ + a_{N-1} z^{-(N-1)} X(z)
\]
\[
H(z) = \frac{Y(z)}{X(z)} = a_0 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4} + \ \cdots \ + a_{N-1} z^{-(N-1)}
\]
\[
H(e^{j{\omega}T}) = a_0 + a_1 e^{-j{\omega}T} + a_2 e^{-2j{\omega}T} + a_3 e^{-3j{\omega}T} + a_4 e^{-4j{\omega}T} + \ \cdots \ + a_{N-1} e^{-(N-1)j{\omega}T}
\]
- FIRフィルタで、LPFを作ってみる。
- 窓関数法ってやつで、フィルタ長は長すぎると面倒なので、N=13に設定。
- 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]);
- 位相特性
- 直線的に変化しているので、多分、遅延はあるけど、歪は少ないってことなのかな
- まぁ、少なくともちゃんと計算出来ているんじゃないかな。
- FIRの係数を配列に入れたら、伝達関数ももっとスマートに計算できるかもしれない。
\[
y[n] = \sum_{k=0}^{N-1} a[k] x[n-k]
\]
\[
H(e^{j{\omega}T}) = \sum_{k=0}^{N-1} a[k] e^{-kj{\omega}T}
\]
sum(a[k]*exp(-k*j*\omega*T),k,0,N-1);
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]);
- 結果は同じになる。このほうが、係数とフィルタ長だけ指定するだけなので簡単。
- 配列じゃなくて、リストで係数を定義するようにした
- 結果は変わらない
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を使う1ようにした。
- 訂正：リストは配列と同様に使えるんだけど、添字が１から始まること、あとNのような変数を使えないのでsumの式を修正
- 配列のプログラムを実行後、そのまま実行したため間違いに気がつかなかった。そのため kill(all)も追加した。（170615）
- 係数をリストにしたので、係数のプロットが簡単2になった。
x: makelist(n,n,0,12);
wxplot2d([discrete,x,a],[style,impulses]);
x: makelist(n,n,0,12);
wxplot2d([discrete,x,a],[style,[points,2,1,1]]);
1Windows用の、wxMaximaじゃないと使えないけれども。
2ということで、以前無理やりプロットしたんだけど、それは消した。
[top]
- このIIRフィルタ（１セクション）の差分方程式は、この構成からはよく分からないけど、
\[
y_i = \frac{1}{a_0} [ - \sum_{j=1}^{N-1} a_j y_{i-j} + \sum_{k=0}^{N-1} b_k x_{i-k} ]
\]
- それに、構成の似ている ノッチフィルタの差分方程式から類推すると、
- たぶん差分方程式は、以下のような式になるんじゃないかと思う3。
\[
y[n] = k ( - (b_1 y[n-1] + b_2 y[n-2]) + (a_0 x[n] + a_1 x[n-1] + a_2 x[n-2]) )
\]
\[
Y(z) = k ( - (b_1 z^{-1} Y(z) + b_2 z^{-2} Y(z) ) + (a_0 X(z) + a_1 z^{-1} X(z) + a_2 z^{-2} X(z) ) )
\]
\[
(1/k + b_1 z^{-1} + b_2 z^{-2}) Y(z) = (a_0 + a_1 z^{-1} + a_2 z^{-2}) X(z)
\]
\[
H(z) = \frac{Y(z)}{X(z)} = \frac{a_0 + a_1 z^{-1} + a_2 z^{-2}}{1/k + b_1 z^{-1} + b_2 z^{-2}}
\]
\[
H(e^{j{\omega}T}) = \frac{a_0 + a_1 e^{-j{\omega}T} + a_2 e^{-2j{\omega}T}}{1/k + b_1 e^{-j{\omega}T} + b_2 e^{-2j{\omega}T}}
\]
- 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特性は単純に落ちるだけで、FIRのLPFの特性とはかなり違う。
- 位相特性
- IIRのLPFも、周波数特性は直線的なので、歪は少ないようだ。
- 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]);
- BPFの場合、ナイキスト周波数だけでなく、DC（f=0）も計算出来ないようなので、f=1～499Hzとした
- ゲイン特性
- 一応、BPFにはなっているが、通過域が平坦になっていない。。なんか間違えてるのか。。。
3自信はない。。
[top]
- MaximaでのFIRフィルタの計算を試してみたが、ちゃんと計算できそう。
- 計算方法は移動平均と同様なので簡単だったが、設計は難しそうだ。
- でもIIRのBPFのグラフが、 DF-Designでの特性グラフとかなり違っている
- IIRの伝達関数が、なにか間違っているみたいだけど、なにを間違えているのか良くわからない。。
- IIRは、BPFがおかしいから、多分LPFも怪しい。
- IIRのブロックダイアグラムからの、差分方程式の導出が間違っているのか
- それとも、セクションを直列接続する場合の計算方法がおかしいのか
- 或いは、なにか他にも間違っているのか。。
- まぁ、当初の目的であるFIRフィルタの計算はちゃんと出来てるっぽいので、今はよしとする。
- いままで位相の変化って、あまり気にしたことなかったので、ちょっと考えてみよう
- フィルタじゃなくても、伝送路を通るだけでも位相は変化する。
- 考えてみると、当然か。。
- 単純な伝送路で、ある周波数で１波長分の変化があるとすると、
- ２倍の周波数では２波長分、３倍の周波数では３波長分位相が変化する。
- 逆に言うと、周波数に比例して位相が変化するフィルタならば、
- 単純な伝送路を通った場合と同じことになり、位相が遅れているだけで、位相の歪は無いと言っても良いと思う。
- ディジタルフィルタ設計プログラム DF-Design →石川高専の山田洋士研究室のページ・FIR/IIRの設計・解析プログラム
[top]
[電子工作関連に戻る]
⇒ Disqusの広告がうるさすぎるので基本は非表示にしました
