一からデジタルフィルタ設計、Pythonで実装してみる（その２）

デジタルフィルタを一から設計してPythonで実装しよう！（タイトルまんま）
…という企画の続き。

残る作業は
・数式内の定数を決める
・差分方程式を元にコーディング

たったこれだけの簡単なお仕事。


まず、数式内の定数を定める。

これは一番初めに決めた「信号から1000Hz以下の周波数成分をブっこ抜く」、という部分とデジタルフィルタの伝達関数から求めていく。

まずデジタルフィルタの伝達関数Hd(z)から周波数特性を得るため、z=e^(j2πfT)とする。
そしてその関数を用いて周波数特性を求めグラフにしたものがこちら。





遮断周波数(振幅が最大値に比べて√2/2になる周波数)が1000Hz付近になるようなRとCになっている。
R=1000、C=0.2uで決定。


これを差分方程式に代入してIIRフィルタとしてコーディングしてみる。

#coding:utf-8

from pylab import *
from scipy import *

def LPF(array):
    output = []
    
    R = 1000.
    C = 0.0000002
    Ts = 1 / 50000.

    for i in range(len(array)):
        if i == 0:
            output.append(1/(R*C)*Ts*array[0])
        else:
            output.append(1/(R*C)*Ts*array[i] + exp(-1/(R*C)*Ts)*output[i-1])

    return output

Fs = 50000.
Ts = 1 / Fs
t = arange(0,1,Ts)
f = 500

s1 = sin(2 * pi * f * t)
s2 = sin(2 * pi * 2 * f * t)
s3 = sin(2 * pi * 5 * f * t)

plot(t[0:200],LPF(s1)[0:200])
plot(t[0:200],LPF(s2)[0:200])
plot(t[0:200],LPF(s3)[0:200])
xlabel("time(s)")
legend(("500Hz","1000Hz","2500Hz"))

show()

結果がコチラ。





周波数特性と比較してみるとその通りにフィルタリングされている。

実際に周波数特性を求めたり、フィルタを作って分かったこと。
それは、アナログのローパス→インパルス不変法というデジタルフィルタの作り方だと、良い特性を持つLPFが作れないということ。
遮断周波数を過ぎてもすぐに減衰しない。

ということで実用的なデジタルフィルタを構築するには他の方式も勉強しておく必要があるようだ。