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

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

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

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Python

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

Q&A

解決済

2回答

853閲覧

バンドストップフィルタ(ノッチフィルタ)を設計したい。

Riri09020500

総合スコア6

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Python

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

0グッド

0クリップ

投稿2022/11/25 19:03

編集2022/11/28 20:38

前提

以下のコードと同様にして、バンドストップフィルタを設計したいです。

実現したいこと

現状は以下のコードにより、カットオフ周波数間の帯域のみをのこし、それ以外を0にする処理を行なっております。
F2[(freq > fc)] = 0
F2[(freq < fc2)] = 0

逆の処理として、カットオフ周波数間のみを0にする方法を教えていただきたいです。
理想としては F2[(freq > fc) and (freq < fc2)] = 0 のような処理をしたい

発生している問題・エラーメッセージ

エラーメッセージ

該当のソースコード

N = len(df) # サンプル数 dt = 0.004 # サンプリング間隔 # fq1, fq2 = 5, 40 # 周波数 fc = 60 # カットオフ周波数 fc2 = 0.2 t = np.arange(0, N*dt, dt) # 時間軸 freq = np.linspace(0, 1.0/dt, N) # 周波数軸 # F3=np.fft.fft(df[4]) for z in df.columns: F = np.fft.fft(df[z]) # 正規化 + 交流成分2倍 F = F/(N/2) F[0] = F[0]/2 df4[z]=F # 配列Fをコピー F2 = F.copy() # ローパスフィル処理(カットオフ周波数を超える帯域の周波数信号を0にする) F2[(freq > fc)] = 0 F2[(freq < fc2)] = 0 df5[z]=F2 # 高速逆フーリエ変換(時間信号に戻す) f2 = np.fft.ifft(F2) # 振幅を元のスケールに戻す f2 = np.real(f2*N) df3[z]=f2

試したこと

ここに問題に対して試したことを記載してください。

補足情報(FW/ツールのバージョンなど)

ここにより詳細な情報を記載してください。

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

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

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

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

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

melian

2022/11/25 19:21

F2[(fc2 < freq)&(freq < fc)] = 0 でしょうか。
Riri09020500

2022/11/27 06:57

こちらでできました。 ありがとうございます。
guest

回答2

0

質問に記載されてる

現状は以下のコードにより、カットオフ周波数間の帯域のみをのこし、それ以外を0にする処理を行なっております。

のコードが間違ってます

PythonでFFT実装!SciPyのフーリエ変換まとめ
の「補足:ミラーリングと周波数軸について」を見てください
「ミラーリングとナイキスト周波数」の横軸が周波数のグラフの、「ミラーリングされている」と書かれてるナイキスト周波数よりも高い周波数のスペクトルデータは、「負の周波数とは?」のグラフで説明されてるように負の周波数のデータなので、それを考慮した下記のような処理が必要です

python

1N = len(df) # サンプル数 2dt = 0.004 # サンプリング間隔 3# fq1, fq2 = 5, 40 # 周波数 4fc = 60 # カットオフ周波数 5fc2 = 0.2 6t = np.arange(0, N*dt, dt) # 時間軸 7freq2 = np.fft.fftfreq(N, dt) # 周波数軸 8# F3=np.fft.fft(df[4]) 9for z in df.columns: 10 11 F12 = np.fft.fft(df[z]) 12 # 正規化 13 F12 = F12/N 14 df42[z]=F12 15 16 # 配列Fをコピー 17 F22 = F12.copy() 18 19 # フィルタ処理(負の周波数も考慮) 20 F22[(freq2 > fc)] = 0 21 F22[(freq2 < fc2) & (freq2 > -fc2)] = 0 22 F22[(freq2 < -fc)] = 0 23 df52[z]=F22 24 25 # 高速逆フーリエ変換(時間信号に戻す) 26 f22 = np.fft.ifft(F22) 27 28 # 振幅を元のスケールに戻す 29 f22 = np.real(f22*N) 30 31 df32[z]=f22

 
上記で計算した「f22」と、質問に記載のコードで計算した「f2」を、下記のようにしてグラフに重ね書きしたら、結果が異なることが分かると思いますので、ご確認ください

python

1import matplotlib.pyplot as plt 2plt.plot(range(len(f2)), f2, label="f2") 3plt.plot(range(len(f22)), f22, label="f22") 4plt.legend() 5plt.show()

 
この質問の

逆の処理として、カットオフ周波数間のみを0にする方法を教えていただきたいです。

にするには、上記のコードの下記を変更します

python

1 F22[(freq2 > fc)] = 0 2 F22[(freq2 < fc2) & (freq2 > -fc2)] = 0 3 F22[(freq2 < -fc)] = 0

↓ 変更

python

1 F22[(freq2 < fc) & (freq2 > fc2)] = 0 2 F22[(freq2 > -fc) & (freq2 < -fc2)] = 0

 
その場合も、上記のようにして計算した「f22」と、TakaiYさんの回答のコードで計算した「f2」を、グラフに重ね書きしたら、やはり結果が異なることが分かると思いますので、ご確認ください

投稿2022/11/28 11:29

jbpb0

総合スコア7651

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

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

jbpb0

2022/11/28 11:38

上記回答で指摘した負の周波数とは別の観点ですが、 print(1/(dt*N)) を計算した結果と、 print(freq) で表示される「freq」の間隔が一致してないと思いますので、ご確認ください そういう意味でも、 freq = np.linspace(0, 1.0/dt, N) # 周波数軸 の計算は間違ってます
guest

0

ベストアンサー

以下の部分が範囲の両側を0にしているところですよね。

python

1 # ローパスフィル処理(カットオフ周波数を超える帯域の周波数信号を0にする) 2 F2[(freq > fc)] = 0 3 F2[(freq < fc2)] = 0

これを間にしようとして

python

1 # バンドストップフィルタ(周波数範囲の周波数信号を0にする) 2 F2[(freq < fc)] = 0 3 F2[(freq > fc2)] = 0

としてしまうと、全部0になってしまいますよね。
条件として、(freq < fc) かつ (freq > fc2) としなければならないので、こういう場合は「&」演算子を使います。

python

1 # バンドストップフィルタ(周波数範囲の周波数信号を0にする) 2 F2[(freq < fc) & (freq > fc2)] = 0

投稿2022/11/26 06:00

編集2022/11/26 06:02
TakaiY

総合スコア12657

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

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

Riri09020500

2022/11/27 07:00

&演算子で解決だったんですね。ご丁寧にありがとうございます。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問