質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

2回答

3524閲覧

Pythonでバンドパスフィルタ

gymgym

総合スコア97

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

1クリップ

投稿2017/11/17 06:02

編集2017/11/17 08:50

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

Python

1fs = 250 # サンプリング周波数 2nyq = fs / 2 # ナイキスト周波数 3fc1 = 4 / fs # カットオフ周波数 4fc2 = 50 / fs # カットオフ周波数 5numtaps = 125 # フィルタ係数の数 6dt = 0.04 # サンプリング間隔 7t = np.arange(0, N*dt, dt) # 時間軸 8freq = np.linspace(0, 1.0/dt, N) # 周波数軸 9 10# フィルタ係数 11b = scipy.signal.firwin(numtaps, [fc1, fc2], pass_zero=False) 12# scipyによるFIRフィルタ 13data_pca = scipy.signal.lfilter(b, 1, data_pca) 14 15 16F = np.fft.fft(data_pca) 17Amp = np.abs(F) 18 19plt.subplot(211) 20plt.plot(t, F) 21plt.xlim(0, 30) 22plt.subplot(212) 23plt.plot(freq, Amp) 24 25plt.show()

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

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

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

フィルタ前の周波数スペクトル
イメージ説明

よろしくお願いします。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

KSwordOfHaste

2017/11/17 06:13 編集

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

2017/11/17 07:15

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

2017/11/17 07:36

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

2017/11/17 08:17

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

2017/11/17 08:51

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

2017/11/17 08:59

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

2017/11/17 09:21

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

2017/11/18 11:59

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

回答2

0

ベストアンサー

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

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

Python

1import numpy as np 2from matplotlib import pyplot as plt 3 4N = 1024 5 6fs = 250 7dt = 1.0 / fs 8t = np.arange(0, N*dt, dt) 9freq = np.linspace(0, 1.0/dt, N) 10tr = dt * N 11 12y1 = 0.5 * np.sin(t * 2 * np.pi * 50) 13y1 += 0.1 * np.sin(t * 2 * np.pi * (50 / 8)) 14y1 += 0.1 * np.sin(t * 2 * np.pi * (50 / 4)) 15y1 += 0.1 * np.sin(t * 2 * np.pi * (50 / 2)) 16y1 += 0.2 * np.sin(t * 2 * np.pi * (50 * 2)) 17 18F = np.fft.fft(y1) 19Amp = np.abs(F) 20 21plt.subplot(211) 22plt.plot(t, y1) 23 24plt.subplot(212) 25plt.plot(freq, Amp) 26 27plt.show()

イメージ説明

ここから始めると、

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

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

投稿2017/11/17 09:57

KSwordOfHaste

総合スコア18392

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

mkgrei

2017/11/17 10:10

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

2017/11/17 10:13

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

0

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

python

1import numpy as np 2import matplotlib.pyplot as plt 3import matplotlib as mpl 4mpl.rcParams['figure.dpi'] = 200 5from scipy import signal 6 7N = 1000 8fs = 250 9nyq = fs / 2 10fc1 = 4. 11fc2 = 30. 12numtaps = 255 13dt = 1. / fs 14t = dt*np.arange(N) 15freq = np.fft.fftfreq(t.shape[-1])*fs 16 17# フィルタ係数 18b = signal.firwin(numtaps, [fc1, fc2], pass_zero=False, nyq=nyq) 19#b = signal.firwin(numtaps, [fc1, fc2]) 20 21x = np.sum((np.random.random()*np.sin(i*t/(N*dt)) for i in range(N)), axis=1) 22x /= np.max(np.abs(x)) 23# scipyによるFIRフィルタ 24kx = signal.lfilter(b, 1, x) 25 26F = np.fft.fft(kx) 27print(F.shape) 28Amp = np.abs(F) 29 30skip=1 31plt.subplot(211) 32plt.plot(freq[::skip], F[::skip], marker='.', ls='none', label='t-F') 33plt.xlim(0, 50) 34plt.legend() 35plt.subplot(212) 36plt.plot(freq[::skip], Amp[::skip], marker='.', ls='none', label='w-A') 37plt.xlim(0, 50) 38plt.legend() 39plt.show()

投稿2017/11/17 08:23

編集2017/11/17 10:14
mkgrei

総合スコア8560

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

gymgym

2017/11/17 08:42

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

2017/11/17 08:59

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

2017/11/17 09:27

いろいろと間違っていましたね。 プロットの横軸とか。 恥ずかしい過去は上書きで修正しました。 今度こそ正しいはずです。 ポイントは周波数空間で周波数を指定して、ライブラリにナイキスト周波数周りの処理を任せることですね…… これなら、全体の次元のつじつまを合わせるのにがんばらなくてよくなるはず。
gymgym

2017/11/18 11:59

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問