お約束の module import （実際には、~/.pystartup に入っている）

In [35]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

scipy の下の signal(.singaltool) を import する。何故か、scipy.signal とはできない

In [36]:
import scipy.signal as sps

テスト用信号を作成 Xin=0.2cos⁡(2πfnt)+cos⁡(2πfst)

In [37]:
t = np.arange(0.0, 2.0, 0.001)
Xin = 0.2 * np.cos(2*np.pi*50*t) + np.cos(2*np.pi*5*t)
plt.plot(t, Xin)
plt.ylim(-2, 2)
plt.ylabel("Amplitude")
plt.xlabel("Time (sec)")
plt.show()

scipy.signal.iirfilter を使った notch filter 設計のための notch_filter_ba を定義（やっている事は、ナイキスト周波数で正規化しているだけ）

In [38]:
def notch_filter_ba (n, fs, fc, bw, rp, rs, ftype):
    fnyq = fs/2.0
    fl = fc - bw/2
    fh = fc + bw/2
    wn = [fl/fnyq, fh/fnyq]
    b, a = sps.iirfilter(n, wn, rp=rp, rs=40, btype='bandstop', ftype=ftype) # PA2
    return b, a
In [39]:
b, a = notch_filter_ba (2, 1000, 50, 4, 1, 40, 'ellip')
In [40]:
Xout = sps.lfilter(b, a, Xin)
plt.plot(t, Xout)
plt.ylim(-2, 2)
plt.ylabel("Amplitude")
plt.xlabel("Time (sec)")
plt.show()

少しの間、ノイズが残るものの、その後は完全に除去できている事、また、振幅が若干小さくなっている。

伝送特性を算出・プロットしてみる。

In [45]:
w, h = sps.freqz(b, a) 
f = w * 1000 / (2 * np.pi)
logh = 20 * np.log10(np.abs(h))
plt.plot(f, logh)
plt.ylabel("Attenuation (dB)")
plt.xlabel("Frequency (Hz)")
plt.show()
plt.plot(f, logh)
plt.xlim(30, 70)
plt.ylabel("Attenuation (dB)")
plt.xlabel("Frequency (Hz)")
plt.show()

以上の計算を、今度は 2次 x 2 の「type-2 チェビチェフフィルター」に対して繰り返してみる。（全部のセルをコピーして、notc_filter_ba に与えるパラメータだけ変更する。）

In [42]:
b, a = notch_filter_ba (2, 1000, 50, 1, 1, 40, 'cheby2')
In [43]:
Xout = sps.lfilter(b, a, Xin)
plt.plot(t, Xout)
plt.ylim(-2, 2)
plt.ylabel("Amplitude")
plt.xlabel("Time (sec)")
plt.show()
In [44]:
w, h = sps.freqz(b, a) 
f = w * 1000 / (2 * np.pi)
logh = 20 * np.log10(np.abs(h))
plt.plot(f, logh)
plt.ylabel("Attenuation (dB)")
plt.xlabel("Frequency (Hz)")
plt.show()
plt.plot(f, logh)
plt.xlim(30, 70)
plt.ylabel("Attenuation (dB)")
plt.xlabel("Frequency (Hz)")
plt.show()