Pythonでバンドパスフィルタ

Question

Pythonで信号をバンドパスフィルタにかけたいと考えています。

fs = 250        # サンプリング周波数
nyq = fs / 2    # ナイキスト周波数
fc1 = 4 / fs # カットオフ周波数
fc2 = 50 / fs   # カットオフ周波数
numtaps = 125   # フィルタ係数の数
dt = 0.04  # サンプリング間隔
t = np.arange(0, N*dt, dt) # 時間軸
freq = np.linspace(0, 1.0/dt, N) # 周波数軸

# フィルタ係数
b = scipy.signal.firwin(numtaps, [fc1, fc2], pass_zero=False)
# scipyによるFIRフィルタ
data_pca = scipy.signal.lfilter(b, 1, data_pca)


F = np.fft.fft(data_pca)
Amp = np.abs(F)

plt.subplot(211)
plt.plot(t, F)
plt.xlim(0, 30)
plt.subplot(212)
plt.plot(freq, Amp)

plt.show()

このようなコードで実行して４Hzから５０Hzの周波数だけ残したと思っています。

実行すると、このようになりました。
下の図を見るとフィルタがかかっていないように思えます。
正しければ、４から５０Hzの周波数帯だけが残ると考えています。

どのようにすれば正しくフィルタをかけることができますか。

フィルタ前の周波数スペクトル

よろしくお願いします。

Comments

[KSwordOfHaste]

波形のみ表示しても周波数成分がどのようになっているか不明瞭です。FFTなどによる周波数スペクトルを比較すべきではないでしょうか？

[ozwk]

フィルタ適用前の周波数スペクトルはありますか?

[gymgym]

周波数スペクトルの画像を追加しました。まだ、足りない点があれば教えてください。よろしくお願いします。

[KSwordOfHaste]

最初の自分のコメントは指摘ミスでした。大変失礼しました。元々フィルター後の周波数スペクトルをプロットされてたんですね。ozwkさんコメントのとおりフィルタ適用前の(data_pca = scipy.signal.lfilter(b, 1, data_pca)をコメント化したときの）実行結果をそのまま追加すればよいと思います。ところで3枚目の図は自分にはよくわかりませんでした・・・

[gymgym]

回答ありがとうございます。コメント化したときの実行結果を更新しておきました。よろしくお願います。

[ozwk]

fs=250Hzならdt=0.004じゃないですかね

[KSwordOfHaste]

dt=1.0/fsと計算で求めた方が間違いがないですね！

[gymgym]

ありがとうございます。とても参考になりました。

Answer 1

えと回答ではないのですが・・・

FFTなどパラメーターが間違ってると混乱するので最初にわかりやすい波形を与えてスケールなどが間違っていないかを確認するとよいと思います。

import numpy as np
from matplotlib import pyplot as plt

N = 1024

fs = 250
dt = 1.0 / fs
t = np.arange(0, N*dt, dt)
freq = np.linspace(0, 1.0/dt, N)
tr = dt * N

y1 =  0.5 * np.sin(t * 2 * np.pi * 50)
y1 += 0.1 * np.sin(t * 2 * np.pi * (50 / 8))
y1 += 0.1 * np.sin(t * 2 * np.pi * (50 / 4))
y1 += 0.1 * np.sin(t * 2 * np.pi * (50 / 2))
y1 += 0.2 * np.sin(t * 2 * np.pi * (50 * 2))

F = np.fft.fft(y1)
Amp = np.abs(F)

plt.subplot(211)
plt.plot(t, y1)

plt.subplot(212)
plt.plot(freq, Amp)

plt.show()

ここから始めると、

• ナイキスト周波数125で左右対称になっている
• 6.25, 12.5, 25, 50, 100の各スペクトルが期待どおりでてる

＝＞パラメーターは間違ってなさそうだというのが確認できると思います。

Comments

[mkgrei]

まさにその基礎から確認すべきでした…

[KSwordOfHaste]

最初コードをよく見ずにおかしなコメントをしてしまい申し訳なかったです。自分も反省しています。

Answer 2

データが無いので、てきとーに生成しています。
実空間・周波数空間の変換がミスっていました。
パラメータはダイレクトに入力して、付随して決まるものは全部計算によって出すと、ミスがすくなるかと。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 200
from scipy import signal

N = 1000
fs = 250
nyq = fs / 2
fc1 = 4.
fc2 = 30.
numtaps = 255
dt = 1. / fs
t = dt*np.arange(N)
freq = np.fft.fftfreq(t.shape[-1])*fs

# フィルタ係数
b = signal.firwin(numtaps, [fc1, fc2], pass_zero=False, nyq=nyq)
#b = signal.firwin(numtaps, [fc1, fc2])

x = np.sum((np.random.random()*np.sin(i*t/(N*dt)) for i in range(N)), axis=1)
x /= np.max(np.abs(x))
# scipyによるFIRフィルタ
kx = signal.lfilter(b, 1, x)

F = np.fft.fft(kx)
print(F.shape)
Amp = np.abs(F)

skip=1
plt.subplot(211)
plt.plot(freq[::skip], F[::skip], marker='.', ls='none', label='t-F')
plt.xlim(0, 50)
plt.legend()
plt.subplot(212)
plt.plot(freq[::skip], Amp[::skip], marker='.', ls='none', label='w-A')
plt.xlim(0, 50)
plt.legend()
plt.show()

Comments

[gymgym]

バンドパスフィルタに対しての質問になってしまうかもしれないのですが、今回の場合１０Hzから５０Hzの周波数帯以外０にすることが目的だと考えています。しかし、下の図をみると３Hzから２３Hzの周波数帯を０にしているように見えます。どのようにしたら、１０Hzから５０Hzの信号のみ通すことができますか。

[mkgrei]

私の見落としですね。
プロットの横軸が正しくないせいですね。
修正します。

[mkgrei]

いろいろと間違っていましたね。
プロットの横軸とか。
恥ずかしい過去は上書きで修正しました。
今度こそ正しいはずです。

ポイントは周波数空間で周波数を指定して、ライブラリにナイキスト周波数周りの処理を任せることですね……
これなら、全体の次元のつじつまを合わせるのにがんばらなくてよくなるはず。

[gymgym]

ありがとうございました。とても参考になりました。