デジタルフィルタを計算する
by K.I
2014/07/22〜
Index
- 数値計算ソフトの Maximaを使って、デジタルフィルタの特性をいろいろ計算してみた
- 差分方程式さえ分かれば、あとの計算は Maximaにお任せで、特性を求められる
- 計算は苦手なので、ちょっと前からMaximaを時々使うようになった。
- とりあえず、自分で使うための備忘録として、メモしておこうと思う。
- Maximaのメモのつもりで書き始めたが、デジタルフィルタの事しか書いてなかったので、
- 結局、題名はこの名前に変更した。まぁ、Maximaのことしか書いてないとも言えるんだけど。。
- 「プログラムでデジタルフィルタ」で、C言語のプログラムで作ったデジタルフィルタが、
- Maximaの計算結果と合っているか、比較してみたので、そちらも参照ください。
- Maximaじゃないけど、電気回路計算の基礎的なこと
- 電気信号の正弦波を、複素指数関数で表す(オイラーの公式)


- コンデンサやインダクタのインピーダンスを、このように簡単に表せる

- このように複素数で計算することで、計算が簡単になる
- デジタルフィルタの特性を表す、伝達関数の定義は以下のようになっている。
- サンプリングした離散データについて入出力の関係を表している
 = \frac{Y(z)}{X(z)}&chf=bg,s,ffffff00)
- H(z)が伝達関数、X(z)が入力信号x[n]のz変換、Y(z)が出力信号y[n]のz変換
ここで、z変換というのが出てきたが、以下の性質を持っている
 = z^{-k}Z(x[n])&chf=bg,s,ffffff00)
- 伝達関数を解析すれば、フィルタの特性を調べることが出来る。
- まぁ、これだけじゃ何のことだかさっぱり分からないので、実際のフィルタについて考えてみよう
[top]
- じゃ、まずRCフィルタ(RC直列回路)について考えてみる
- 以前、 同様の考察をしてるので、その復習になるけど
- RC直列回路の出力電圧は、抵抗とCのインピーダンスから計算すると、

- RCフィルタは、アナログフィルタだけど、
- サンプリングした離散データだけ考えることで、同様の特性を持つデジタルフィルタになる
- x[t]=Vin, y[t]=Vout とおいて、周期Tでサンプリングしたデータを扱うことにする
- jωは微分として考える(nはnT,n-1は(n-1)Tだけど、数列として考えるのでTは省略)
 %2B (1 - \frac{CR}{CR%2BT}) x(n)&chf=bg,s,ffffff00)
- 係数 CR/(CR+T)をaで置き換えて、差分方程式を求める
 x[n]&chf=bg,s,ffffff00)
 = \frac{Y(z)}{X(z)} = \frac{1-a}{1-az^{-1}}&chf=bg,s,ffffff00)
- 伝達関数の定義から、RCフィルタの差分方程式をz変換して、伝達関数を求める。
 = \frac{1-a}{1-ae^{-j{\omega}T}}&chf=bg,s,ffffff00)
- 手計算だと、オイラーの公式を使って三角関数に展開して、絶対値と偏角を求めるんだけど、
- Maximaを使って、いきなり数値計算させてみよう
(1-a)/(1-a*exp(-j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(0.5,a,%);
H:%;
plot2d(cabs(H),[f,0,10000],[logx]);
plot2d(carg(H),[f,0,10000],[logx]);
- サンプリング周波数1kHz(T=0.001)、a=0.5で計算してみた
- サンプリング周波数が1kHzなので、500Hz以上は無意味だが、
- 折り返して、同じ特性が繰り返されていることがわかる
- 普通、計算プログラムを書いても、計算誤差でこんなにキレイに繰り返しにならなかったりするんだけど、
- Maximaは、出来るだけ誤差が少なくなる様に計算してくれるので有り難い。
[top]
- 移動平均が、LPFになっていることは大体分かるんだけど、どんな特性を持っているんだろう
- そもそも移動平均は、デジタルフィルタなので、差分方程式は簡単に書ける

 = \frac{1 %2B e^{-j{\omega}T} %2B e^{-2j{\omega}T}}{3}&chf=bg,s,ffffff00)
- またMaximaを使って、いきなり数値計算させてみよう
- 今度は、defineで伝達関数を関数として定義した
(1+exp(-j*omega*T)+exp(-2*j*omega*T))/3;
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(0.5,a,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx]);
plot2d(carg(H(f)),[f,0,500],[logx]);
- ゲイン特性
- 横軸(周波数)を対数スケールじゃなくしてみると、
- 位相特性は、実際はリニア(直線的)に変化していることが分かる
- 位相特性は、途中で不連続な変化をしている様に見える
- そんなわけ無いはずなので、複素平面上に周波数の変化に対する伝達関数の軌跡をプロットしてみる
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,0,500]],[x,-1,1],[y,-1,1]);
- 不連続ではないが、周波数が大きくなるにつれて、ゲインが1から下がって行って、
- ゲインが0になった時に、位相が反転している。実際はなめらかに向きを変えているという感じだ
] }{N}&chf=bg,s,ffffff00)
- 前項と同様に計算すると、伝達関数は、こんな感じで良いのかな?
 = \frac{1}{N} \sum_{n=1}^{N}{e}^{-(n-1) j \omega T}&chf=bg,s,ffffff00)
- Maximaで、N=2,4,8,10について、ゲインと位相、そして複素平面上の軌跡を求める
- サンプリング周波数1kHz(T=0.001)、a=0.5で計算
sum(exp(-(n-1)*j*\omega*T),n,1,N)/N;
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(0.5,a,%);
define(H(f,N),%);
plot2d([cabs(H(f,2)),cabs(H(f,4)),cabs(H(f,8)),cabs(H(f,10))],[f,0,500],[y,0.001,1],[logx],[logy]);
plot2d([carg(H(f,2)),carg(H(f,4)),carg(H(f,8)),carg(H(f,10))],[f,0,500]);
plot2d([[parametric,realpart(H(f,2)),imagpart(H(f,2)),[f,0,500]],
[parametric,realpart(H(f,4)),imagpart(H(f,4)),[f,0,500]],
[parametric,realpart(H(f,8)),imagpart(H(f,8)),[f,0,500]],
[parametric,realpart(H(f,10)),imagpart(H(f,10)),[f,0,500]]],[x,-1,1],[y,-1,1],[nticks,100]);
- ゲイン特性(N=2,4,8,10) →周波数、ゲインは対数表示
- 位相特性(N=2,4,8,10) →周波数、位相はリニア表示
- ゲイン特性も位相特性も、それだけ見ると異様な変化をしているように感じるが、
[top]
- じゃ、なんか分からない不明のデジタルフィルタの特性を求めてみよう
- 以前、計算出来ずにギブアップしたものなんだけど、差分方程式から伝達関数までは求めてある
 = \frac{1-e^{-5j{\omega}T}}{5(1-e^{-j{\omega}T})}&chf=bg,s,ffffff00)
- この後、三角関数に展開して数値計算するのが、出来なかったんだよな〜。。
- Maximaなら、伝達関数があれば、ただ計算するだけだ。
(1-exp(-5*j*omega*T))/(5*exp(-j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(0.5,a,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx],[logy]);
plot2d(carg(H(f)),[f,0,500]);
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,0,500]],[x,-1,1],[y,-1,1]);
- これは、ノッチフィルタって奴なのかな?
- とにかく、伝達関数があれば、すぐに特性が確認できるのは有難い。
- このフィルタって、なんか微分回路の特性に似てるし、式はSINCフィルタみたいだな〜と思って、
- 改めて式を見直してみると、Maximaの式の入力、間違えてるじゃん。。。
- やり直し。。(140901追記)
(1-exp(-5*j*omega*T))/(5*(1-exp(-j*omega*T)));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(0.5,a,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx],[logy]);
plot2d(carg(H(f)),[f,0,500]);
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,0,500]],[x,-1,1],[y,-1,1]);
- これはM=5のSINCフィルタを出力がx1になるように5で割ったもの、
- 言い換えれば、5個の移動平均フィルタだったのか。。
- やっと、正体がわかりました。。。
[top]
- 微分を、ある時間のデータの変化量を求めることと考えると、
- 以下の図の様に、ある時間分遅延させたデータを引くことで、微分していることになる

差分方程式をZ変換して、伝達関数を求める
 = 1 - e^{-jM{\omega}T}}&chf=bg,s,ffffff00)
- Maximaで、伝達関数から微分回路の特性を計算してみる。
- T=0.001, M=16 とする
1-exp(-j*M*omega*T);
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(16,M,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx],[logy],[y,0.0001,10]);
plot2d(carg(H(f)),[f,0,500]);
- ゲイン特性
- 微分回路は、移動平均フィルタに似た特性だが、違いもある
- ゲイン特性の高さが、ずっと変わらない(全く同じ繰り返しに見える)
- また、原理上、DC(ω=0)ではゲインが0になる
- 櫛(くし)のような特性なので、くし型フィルタとか、コムフィルタとか呼ばれる。
- 微分回路はフィードフォワード型コムフィルタとも呼ばれる
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,0,500]],[x,-1,1],[y,-1,1]);
- 最初の部分だけ、60Hzぐらいまで、プロットしてみる
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,0,60]],[x,-1,1],[y,-1,1]);
- これが本当だな。この円をぐるぐる回ってるわけか。。
- gnuplotのparametricは、プロット可能な頂点数が決まってるらしい1。これは要注意だな〜
1たぶん設定できるんだろうけど、分からないので。。
[top]
- 積分は、サンプリングしたデータを全て加算すれば良い、
- 1つ前のデータを保持しておいて新しいデータを加算、つまり以下の図の回路で積分となる

差分方程式をZ変換して、伝達関数を求める
 = \frac{1}{1 - e^{-j{\omega}T}}&chf=bg,s,ffffff00)
- Maximaで、伝達関数から積分回路の特性を計算してみる。
- T=0.001とする
1/(1-exp(-j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx],[logy],[y,0.1,1000]);
plot2d(carg(H(f)),[f,0,500]);
- ゲイン特性
- これを見ると、積分回路は周波数が高いとほとんどゲインが無いが、
- 周波数が低くなると、ゲインが大きくなり、DCでゲインは無限大になる
- 周波数特性的にはローパスフィルタとして働く
- 但し、ω=0で無限大のゲインを持つので、DC成分を入力すると簡単に飽和してしまう
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,0,500]],[x,-1,1])
- 積分回路は、DCでは無限大になってしまうので、周波数は1Hzからに変更する
plot2d([parametric,realpart(H(f)),imagpart(H(f)),[f,1,500]],[x,-1,1])
- かなり極端な特性の回路だ。これは、なかなかに使い辛い。。
[top]
- 積分回路を、微分回路と同様に多段に遅延させてみる
- 全部のデータを加算するんじゃなくて、ちょっと間を空けたデータを足す感じ
- 多段という言い方はちょっと変だけど、区別するため以降、便宜上2そう呼ぶことにする
- 積分回路は、遅延1段で使うことが多いけど、多段の場合はどうなるのかな。

 = \frac{1}{1 - e^{-jM{\omega}T}}&chf=bg,s,ffffff00)
- Maximaで、伝達関数から多段積分回路の特性を計算してみる。
- T=0.001、 M=16とする
1/(1-exp(-j*M*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(16,M,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx],[logy],[y,0.1,1000]);
plot2d(carg(H(f)),[f,0,500]);
- ゲイン特性
- 微分回路と、まるっきり逆の特性になっている。
- フィードバック型のコムフィルタとも呼ばれるらしい。
- 複素平面上の軌跡は、1段の積分回路の繰り返しみたいだけど、
- 無限大が周期的に出るみたいでプロットしづらいので、止めた。。
2あくまで区別するために、自分がこの文章の中だけで便宜上そう呼んでいるだけ。一般的な言い方じゃないので注意。
[top]
- ω=0で、振幅特性が無限大になる積分回路と、ゼロになる微分回路を
- この場合、微分回路の方が遅延段数が多いので、
- 基本的には、積分回路よりも、微分回路の規模が大きくなりやすい
- SINCフィルタの伝達関数は、積分回路と微分回路の伝達関数の積で求められる
 = \frac{1 - e^{-jM{\omega}T}}{1 - e^{-j{\omega}T}&chf=bg,s,ffffff00)
- DC(ω=0)の時、積分回路はゲインが無限大、微分回路はゲイン0となるので、
- 矛と盾みたいで、どうなるのかってことなんだけど。。
- SINCフィルタの特性を、Maximaで計算してみよう
- T=0.001、 M=16とする
(1-exp(-j*M*omega*T))/(1-exp(-j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
subst(16,M,%);
define(H(f),%);
plot2d(cabs(H(f)),[f,0,500],[logx],[logy],[y,0.01,100]);
plot2d(carg(H(f)),[f,0,500]);
- ゲイン特性
- ω=0での特性がうまく組合わされて、このような良い感じの振幅特性のフィルタになる。
- あれ?、でもこの特性って、移動平均フィルタに似てないか?
- 移動平均フィルタを個数で割らない場合と、SINCフィルタのゲインが全く同じ様に見える
- 式で証明出来ると良いんだけど、自分には出来ないので、Maximaで比較してみる3
- 3個のデータを足したものと、遅延3の微分回路のSINCフィルタの結果を較べてみる
(1+exp(-j*omega*T)+exp(-2*j*omega*T))$
ratsimp(%);
(1-exp(-j*3*omega*T))/(1-exp(-j*omega*T))$
ratsimp(%);
(1+exp(-j*omega*T)+exp(-2*j*omega*T)+exp(-3*j*omega*T))$
ratsimp(%);
(1-exp(-j*4*omega*T))/(1-exp(-j*omega*T))$
ratsimp(%);
(1+exp(-j*omega*T)+exp(-2*j*omega*T)+exp(-3*j*omega*T)+exp(-4*j*omega*T)+exp(-5*j*omega*T)+exp(-6*j*omega*T)+exp(-7*j*omega*T))$
ratsimp(%);
(1-exp(-j*8*omega*T))/(1-exp(-j*omega*T))$
ratsimp(%);
(1+exp(-j*omega*T)+exp(-2*j*omega*T)+exp(-3*j*omega*T)+exp(-4*j*omega*T)+exp(-5*j*omega*T)+exp(-6*j*omega*T)+exp(-7*j*omega*T)+exp(-8*j*omega*T)+exp(-9*j*omega*T)+exp(-10*j*omega*T)+exp(-11*j*omega*T)+exp(-12*j*omega*T)+exp(-13*j*omega*T)+exp(-14*j*omega*T)+exp(-15*j*omega*T))$
ratsimp(%);
(1-exp(-j*16*omega*T))/(1-exp(-j*omega*T))$
ratsimp(%);
- SINCフィルタの出力を微分回路の遅延Mで割れば、M個の移動平均と全く同じフィルタになる様だ。
- 全く違う考え方で作られたフィルタが、全く同じ特性になるのは、全く不思議だ。。
- ゲインがx1になるように1/Mすると、SINCフィルタの伝達関数はこうなる
 = \frac{1 - z^{-M}}{M(1 - z^{-1})}&chf=bg,s,ffffff00)
- これは、個数Mの移動平均フィルタの伝達関数でもあることになる
- あらためて調べてみると、移動平均とSINCフィルタが同じというのは、常識みたいだ。
- なんか似たような特性だな〜とは思ったことあるけど、ずっと別物だと4思っていた。。
3まぁ、手抜きですね。
4だって、全然違う構成なんだもの。
[top]
- デジタル共振器の構成
- これは、特定の周波数だけ、強調するというものらしい
- この式って、b0入ってないじゃないか。あぁ、b0 = 1 - a1 - a2 なのか。。
- でも、b0掛かってないのにy[n]として使ってるのは変じゃないか?

- 無理やり、こんな風に考えてみたけど、こんなんで良いのか。。。
 = \frac{1 - a_1 - a_2}{1 - a_1 e^{-j{\omega}T} - a_2 e^{-2j{\omega}T}}&chf=bg,s,ffffff00)
- この式で、a2が1よりわずかに小さい場合、
- 共振器の共振周波数F0、共振の帯域幅5B0を使うと、係数a1,a2は以下のように表すことが出来るらしい。

- Maximaを使って、数値計算させてみた
- F0=(サンプリング周波数fs)/10、サンプリング周期Tを1kHzとして、B0=1, 10, 20 の時の特性を求める
(1-a1-a2)/(1-a1*exp(-j*omega*T)-a2*exp(-2*j*omega*T));
subst(2*exp(-%pi*B0*T)*cos(2*%pi*F0*T),a1,%);
subst(-exp(-2*%pi*B0*T),a2,%);
subst(1/(10*T),F0,%);
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
define(H(f,B0),%);
plot2d([cabs(H(f,1)),cabs(H(f,10)),cabs(H(f,20))],[f,0,500],[logx]);
plot2d([carg(H(f,1)),carg(H(f,10)),carg(H(f,20))],[f,0,500],[logx]);
- ゲイン特性
- 非常に鋭い特性のフィルタであることが分かる →フィルタというより、特定の周波数を増幅するという感じだけど
5通常、共振器の帯域幅は、共振周波数における振幅に対して、幅が1/√2(-3dB)になる点での周波数幅で定義される。(ということらしい)
[top]
- ノッチフィルタの構成
- これは、共振器とは全く逆に、特定の周波数のみ取り除くフィルタになる
- いまいち、上のブロックダイアグラムが、
- この式になるのが良く分からないんだけど、まぁ良しとする。
 = \frac{c_0(1 %2B b_1 e^{-j{\omega}T} %2B e^{-2j{\omega}T})}{1 - a_1 e^{-j{\omega}T} - a_2 e^{-2j{\omega}T}}&chf=bg,s,ffffff00)
- a1,a2は、共振器と同じ。b1, c0は、以下のように定義される

- Maximaを使って、数値計算させてみた
- F0=(サンプリング周波数fs)/10、サンプリング周期Tを1kHzとして、B0=1, 10, 20 の時の特性を求める
c0*(1+b1*exp(-j*omega*T)+exp(-2*j*omega*T))/(1-a1*exp(-j*omega*T)-a2*exp(-2*j*omega*T));
subst((1-a1-a2)/(2+b1),c0,%);
subst(2*exp(-%pi*B0*T)*cos(2*%pi*F0*T),a1,%);
subst(-exp(-2*%pi*B0*T),a2,%);
subst(-2*cos(2*%pi*F0*T),b1,%);
subst(1/(10*T),F0,%);
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
define(H(f,B0),%);
plot2d([cabs(H(f,1)),cabs(H(f,10)),cabs(H(f,20))],[f,0,500],[logx]);
plot2d([carg(H(f,1)),carg(H(f,10)),carg(H(f,20))],[f,0,500],[logx]);
- ゲイン特性
- ちゃんと、ノッチになってる。
- 原理的には良くわからないんだけど、これも使い道がありそう。
[top]
- オールパスフィルタの構成
- 振幅を変化させずに、位相のみ変化させるフィルタらしい。
 = \frac{e^{-j{\omega}T} - a}{1 - a e^{-j{\omega}T}}&chf=bg,s,ffffff00)
- Maximaを使って、数値計算させてみた
- サンプリング周期Tを1kHzとして、a=-0.9, -0.5, 0.5, 0.9 の時の特性を求める
(exp(-j*omega*T)-a)/(1-a*exp(-j*omega*T));
subst(%i,j,%);
subst(2*%pi*f,omega,%);
subst(0.001,T,%);
define(H(f,a),%);
plot2d([cabs(H(f,-0.9)),cabs(H(f,-0.5)),cabs(H(f,0.5)),cabs(H(f,0.9))],[f,0,500],[logx],[y,0.99999,1.00001]);
plot2d([carg(H(f,-0.9)),carg(H(f,-0.5)),carg(H(f,0.5)),carg(H(f,0.9))],[f,0,500],[logx]);
- ゲイン特性
- ゲインが全く変化せず、位相だけが変化していることが分かる。
- 位相の変化が、ゲインに全く影響しないのは、何か不思議だなぁ
- 位相を変えずにゲインを変化させるって、出来ないんじゃないだろうか。。逆は出来るのか。。。
- でも、これって、いったい何に使うんだろう?
[top]
- というか、感想だけど、Maximaは計算誤差が出来るだけ少なくなるように計算してくれるので、
- 特に、複素数の計算が簡単に出来るので、電子回路やデジタルフィルタの特性計算も簡単に出来る。
- でも式の意味が分からないまま、計算出来てしまうので、ちょっと怖い気がする
- まぁ、あくまで道具なので、わかる範囲でのみ使うようにするべきかもしれない。
- 便利すぎて、使いすぎると計算能力は落ちるかも。(自分はいまさらだけどね。。)
- まぁ、でもデジタルフィルタって、面白いなぁ。
- 実際には、サンプリングとかいろいろ面倒なことはあるんだけど、
- 計算してる分には、とても楽しい。
- さらに、このMaximaで特性を計算したフィルタについて、 「プログラムでデジタルフィルタ」 という記事で、
- Cでデジタルフィルタを記述、fftwライブラリを使って、FFTを掛けてみた。
- これも計算上ではあるんだけど、実際にデータを処理して、Maximaの計算結果は合ってるのか?ってことで。
[top]
[電子工作関連に戻る]
⇒ Disqusの広告がうるさすぎるので基本は非表示にしました