ディジタルフィルタってなんだろう

by K.I

2008/07/xx

Index

概要
移動平均
ローパスフィルタ
- RCフィルタをディジタルフィルタに変換
- ステップ応答
- 伝達関数
- 周波数応答
自分で計算してみよう
- デジタルフィルタ関連の記事

概要

ディジタルフィルタは、以前から興味があったけど、DSPとか使わないと出来ないかなぁと思っていた。
PastelMagicの掲示板で、ごく低速の信号でディジタルフィルタを使うという話題が出たので、ちょっと試してみたくなった。
Marsさんという方が紹介していた式（4msサンプリングで50Hzならば、ということだが。。。）
どんなフィルタになるのか全然イメージが湧かないけど、低速なら適当に試してみることが出来そうだ。
ちょっと模式的に描いてみる。

まぁ、わかんないなりに考えてみよう。

これは4msサンプリングとすれば、4ms×5＝20ms間隔で相殺してる。
つまり、1/20ms＝50Hzの信号を消してしまうと考えられる。
あとは、前回のデータの1/5の割合で加算してるから、感覚的にローパスになっている感じがする。

でもこれ以上は、わかんないなぁ。

とにかく、いろいろ調べてみよう。こことか、ここなんか参考になりそうだ。

[top]

移動平均

サンプリングしたデータについて処理するとしたら、最初に思いつくのは単純に平均をとることだよね。例えば、

これなら自分にもなんとなく分かる。ローパスになってる気がする。

こんな風に、有限の範囲のデータに係数を掛けて加算して作るフィルタをFIRフィルタというようだ。移動平均も、FIRフィルタの一種ということらしい。

\[ y(n) = ( 5y(n-1) + x(n) - x(n-5) ) / 5 \]

こうしてみると、これは立派なデジタルフィルタみたいだ。意識せずにデジタルフィルタを使っていたんだなぁ。

[top]

ローパスフィルタ

なんとなくローパスになってるのは分かるんだけど、特性は正確には分からない。

やっぱり、理論的に知りたいよね。

ディジタル信号処理の基礎（CQ出版）を買ってきたので、それを元に勉強しよう。

この本の最初にRCフィルタの式が出ている。

\[ \frac{dy(t)}{dt} + \frac{1}{CR}y(t) = \frac{1}{CR}x(t) \]

うーん、あまり見たことない式だなぁ。

RCフィルタをディジタルフィルタに変換

RCフィルタ（いわゆる積分回路）をちょっと復習してみよう。

直列RC回路の電圧E(t)は、以下の式で表される。これなら分かる。

\[ E(t) = RI(t) + \frac{q(t)}{C} \]

ここで、電流I(t)は電荷q(t)の変化量だから、

\[ I(t) = \frac{dq(t)}{dt} \]\[ R\frac{dq(t)}{dt} + \frac{q(t)}{C} = E(t) \]

さらに容量Cの両端の電圧をEc(t)とすれば

\[ CEc(t) = q(t) \]\[ RC\frac{dEc(t)}{dt} + Ec(t) = E(t) \]

つまりx=E、y=Ecとすれば、、なるほど、こうなるわけね。

\[ \frac{dy(t)}{dt} + \frac{1}{CR}y(t) = \frac{1}{CR}x(t) \]

また本の説明を追っかけてみよう。

この微分方程式を差分で置き換えるってことなんだけど、いきなり変換してるから良く分からん。
ちょっと悩んだけど、こんな感じで良いのかな。

\[ \frac{y(nT)-y((n-1)T)}{T} + \frac{1}{CR}y(nT) = \frac{1}{CR}x(nT) \]\[ CRy(nT) - CRy((n-1)T) + Ty(nT) = Tx(nT) \]\[ (CR+T)y(nT) - CRy((n-1)T) = Tx(nT) \]\[ y(nT) = \frac{CR}{CR+T}y((n-1)T) + \frac{T}{CR+T}x(nT) \]\[ = \frac{CR}{CR+T}y((n-1)T) + (1-\frac{CR}{CR+T})x(nT) \]

で、係数をaで置き換える

\[ \frac{CR}{CR+T} = a \]\[ y(nT) = ay((n-1)T) + (1-a)x(nT) \]

これが差分方程式ってことなんだ。。それで、Tは一定なので省略して、

\[ y[n] = ay[n-1] + (1-a)x[n] \]

と表すわけね。そして、ブロック図で描くと

あぁ、これは最初のフィルタの後半部分と一緒だね。なるほど～やっぱりローパスフィルタだったんだ。

うーん、RCフィルタをディジタルフィルタで記述出来るっていうのは、自分的にはちょっとしたカルチャーショックだなぁ。

ステップ応答

さらに、本をなぞってみよう。

単位ステップ信号の応答を考える。これは矩形波の応答と同じだから実際にも良くある信号だ。

これはお馴染みの積分回路の過渡応答ってやつだね。微分方程式を解くとこうなる。

\[ y(t) = 1 - exp \left( -\frac{t}{CR} \right), \qquad t \ge 0 \]

ディジタルフィルタの場合は、やはり差分方程式を解いて、以下のようになるらしい。

もちろん、初期値x=1,y=0として、差分方程式に当てはめて順に数値計算しても同じ。

\[ y[n] = 1 - a^n+1 \] ステップ応答

アナログでもディジタルでも同じ¹ってところが面白いなぁ。

伝達関数

さぁ、いよいよどんな特性をもったフィルタなのかってことを知ることが出来そうだ。

伝達関数の定義は以下のようになっている。

\[ H(z) = \frac{Y(z)}{X(z)} \]

H(z)が伝達関数、X(z)が入力信号x[n]のz変換、Y(z)が出力信号y[n]のz変換

ここで、z変換というのが出てきたが、以下の性質を持っている \[ X(z) = Z(x[n]) \]\[ Z(ax_1[n] + bx_2[n]) = aZ(x_1[n]) + bZ(x_2[n]) \]\[ Z(x[n-k]) = z^{-k}Z(x[n]) \]

これを基に、差分方程式をｚ変換して、伝達関数を求める。

\[ y[n] = ay[n-1] + (1-a)x[n] \]\[ Y(z) = az^{-1}Y(z) + (1-a)X(z) \]\[ H(z) = \frac{Y(z)}{X(z)} = \frac{1-a}{1-az^{-1}} \]

ふーん、なるほどね。って、なんだか全然分からないけど。。。

周波数応答

やっぱり、周波数応答を見ないと分かんないよね。難しい理屈は置いといて、zを以下のように置き換えれば、周波数応答が求められるらしい。

\[ z = exp(j{\omega}T) \]\[ \omega = 2{\pi}f \]\[ H(e^{j{\omega}T}) = \frac{1-a}{1-ae^{-j{\omega}T}} \]

なんか見慣れた式になってきた。これは複素数が入っているので、極形式、つまり絶対値と偏角で表す。

\[ e^{jx} = \cos{x} + j\sin{x} \]\[ e^{-jx} = \cos{x} - j\sin{x} \]\[ \cos{x} = \frac{e^{jx} + e^{-jx}}{2} \]\[ \sin{x} = \frac{e^{jx} - e^{-jx}}{2j} \]

オイラーの公式。これを使って展開する。

\[ H(e^{j{\omega}T}) = \left| H(\omega) \right| e^{-j\theta(\omega)} \]\[ \left| H(\omega) \right| = \left| \frac{1-a}{1-ae^{-j{\omega}T}} \right| = \frac{\left| 1-a \right|}{\sqrt{1+a^2-2a\cos{{\omega}T}}} \]\[ \theta(\omega) = tan^{-1}\frac{Im\left[H(e^{j{\omega}T}\right]}{Re\left[H(e^{j{\omega}T}\right]} = -tan^{-1}\frac{a\sin{{\omega}T}}{1-a\cos{{\omega}T}} \]

へー、こんな式になるんだ。

これで周波数応答を計算出来る。素晴らしい。

ゲイン（絶対値）

位相（偏角）

ここ²を参考に（というより、そのまんまだが。。。）Javaで計算してみました³。（T=0.001,a=0.5）

サンプリング周波数が１kHzだから、500Hzぐらいまでしか使えないというように見れば、ちゃんとローパスになってることが分かる。
ん、でも標本化周波数で同じ特性が繰返されるはずじゃないのかなぁ。何か間違っているのか⁴な？

サーバクラッシュで、Javaが使えなくなったので、グラフはMaximaで改めて計算しなおしました。（130813追記）

Maximaでの計算方法については、デジタルフィルタを計算するを参照

¹ディジタルの場合は、aがマイナスもアリなのでちょっと違うらしい。
²リンク先が無くなってしまったので、やり方も分からなくなってしまいました。。
³計算式を完全に理解してないので、間違ってるかも。
⁴これは、Javaで計算した時はそう思ったんだけど、Maximaで計算したグラフは、ちゃんと繰り返しになってるので、計算誤差があったと思われる。

[top]

自分で計算してみよう

じゃ、最初の式はどうなんだろう。

\[ y[n] = ( 5y[n-1] + x[n] - x[n-5] ) / 5 \]\[ Y(z) = z^{-1}Y(z) + ( X[z] - z^{-5}X(z) ) / 5 \]\[ H(z) = \frac{Y(z)}{X(z)} = \frac{1-z^{-5}}{5(1-z^{-1})} \]\[ H(e^{j{\omega}T}) = \frac{1-e^{-5j{\omega}T}}{5(1-e^{-j{\omega}T})} \]

あぁ、計算出来ない→計算苦手。複素指数関数を三角関数に変換するソフトでもないかなぁ。

Mathematicaとかにはありそうなんだけどね。それだけでも計算出来れば良いんだけど。。。

これは、Maximaなら無理やり数値計算してしまうことができることがわかったけど、それはまた別の機会に。。。（130813追記）

デジタルフィルタを計算するに纏めたので、そちらも参照ください（140804追記）

デジタルフィルタ関連の記事

デジタルフィルタを計算する→Maximaを使って、デジタルフィルタの特性を計算してみた
fftwでいろいろ実験→fftwを使って、いろいろ遊んでみました
プログラムでデジタルフィルタ→プログラムで、いろいろなデジタルフィルタを作ってfftwでFFT掛けてみた

[top] [電子工作関連に戻る]

　⇒　Disqusの広告がうるさすぎるので基本は非表示にしました