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

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

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

JupyterLabは、Jupyter notebookの後継の対話型開発環境(IDE)です。データの可視化がインタラクティブで、プラグイン作成により新しいコンポーネントの追加および既存のコンポーネントも統合可能。サーバに閉じているため、データ分析に向いています。

Python

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

Q&A

解決済

PythonでCSVデータにハイパスフィルタをかけたいが適用されない

NMKN
NMKN

総合スコア8

JupyterLab

JupyterLabは、Jupyter notebookの後継の対話型開発環境(IDE)です。データの可視化がインタラクティブで、プラグイン作成により新しいコンポーネントの追加および既存のコンポーネントも統合可能。サーバに閉じているため、データ分析に向いています。

Python

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

2回答

0グッド

0クリップ

749閲覧

投稿2022/10/12 10:10

前提

PythonでCSVデータを解析しています。データを用いてスペクトログラムを作ったところノイズを除去する必要があると考えました。そこで500Hz以下を通さないハイパスフィルタをかけたいと考えてプログラムしたのですが上手くいきません。ですのでご教授いただければと思います。

実現したいこと

PythonでCSVデータにハイパスフィルタをかけてノイズを除去したい

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

エラーは出ておりませんがそもそものデータにフィルタが適用されずデータが生み出されません。
下記の画像の通り下側の画像にフィルターをかけた生データがプロットされるはずなのですがされておりません。下記のサイトを参考にプログラムしました。
(https://watlab-blog.com/2021/08/05/csv-digital-filter/)
イメージ説明

エラーメッセージ import numpy as np from scipy import signal import pandas as pd import matplotlib.pyplot as plt # ハイパスフィルタ def highpass(x, samplerate, fp, fs, gpass, gstop): fn = samplerate / 2 #ナイキスト周波数 wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化 ws = fs / fn #ナイキスト周波数で阻止域端周波数を正規化 N, Wn = signal.buttord(wp, ws, gpass, gstop) #オーダーとバターワースの正規化周波数を計算 b, a = signal.butter(N, Wn, "high") #フィルタ伝達関数の分子と分母を計算 y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける return y # csvから列方向に順次フィルタ処理を行い保存する関数 def csv_filter(in_file, out_file, type): df = pd.read_csv(in_file, encoding='SHIFT-JIS') # ファイル読み込み dt = df.T.iloc[0,1] # 時間刻み # データフレームを初期化 df_filter = pd.DataFrame() df_filter[df.columns[0]] = df.T.iloc[0] # ハイパスの設定----------------------------------------------------------------------------- fp_hp = 520 # 通過域端周波数[Hz] fs_hp = 500 # 阻止域端周波数[Hz] gpass = 3 # 通過域端最大損失[dB] gstop = 40 # 阻止域端最小損失[dB] # 列方向に順次フィルタ処理をするコード for i in range(len(df.T)-1): data = df.T.iloc[i+1] # フィルタ処理するデータ列を抽出 # フィルタ処理の種類を文字列で読み取って適切な関数を選択する if type == 'hp': # ハイパスフィルタを実行 print('wave=', i, ':Highpass.') data = highpass(x=data, samplerate=1 / dt, fp=fp_hp, fs=fs_hp, gpass=gpass, gstop=gstop) else: # 文字列が当てはまらない時はパス(動作テストでフィルタかけたくない時はNoneとか書いて実行するとよい) pass data = pd.Series(data) # dataをPandasシリーズデータへ変換 df_filter[df.columns[i + 1] + '_filter'] = data # 保存用にデータフレームへdataを追加 df_filter.to_csv(out_file) # フィルタ処理の結果をcsvに保存 return df, df_filter # 関数を実行してcsvファイルをフィルタ処理するだけの関数を実行 # type='lp', 'hp', 'bp', 'bs':LowPass, HighPass, BandPass, BandStop df, df_filter = csv_filter(in_file='csv in/signals.csv', out_file='filter.csv', type='hp') # ここからグラフ描画------------------------------------- # フォントの種類とサイズを設定する。 plt.rcParams['font.size'] = 14 plt.rcParams['font.family'] = 'Times New Roman' # 目盛を内側にする。 plt.rcParams['xtick.direction'] = 'in' plt.rcParams['ytick.direction'] = 'in' # グラフの上下左右に目盛線を付ける。 fig = plt.figure(figsize=(10, 7)) ax1 = fig.add_subplot(211) ax1.yaxis.set_ticks_position('both') ax1.xaxis.set_ticks_position('both') ax2 = fig.add_subplot(212) ax2.yaxis.set_ticks_position('both') ax2.xaxis.set_ticks_position('both') # 軸のラベルを設定する。 ax1.set_xlabel('Time [s]') ax1.set_ylabel('Amplitude_Original') ax2.set_xlabel('Time [s]') ax2.set_ylabel('Amplitude_Filtered') # データプロットの準備とともに、ラベルと線の太さ、凡例の設置を行う。 size = len(df.T)-1 for i in range(size): ax1.plot(df.T.iloc[0], df.T.iloc[i+1], label=df.columns[i+1], lw=1) ax2.plot(df_filter.T.iloc[0], df_filter.T.iloc[i + 1], label=df_filter.columns[i + 1], lw=1) ax1.legend() ax2.legend() # レイアウト設定 fig.tight_layout() # グラフを表示する。 plt.show() plt.close() ### 該当のソースコード ```ここに言語名を入力 Python

試したこと

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

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

以下のような質問にはグッドを送りましょう

  • 質問内容が明確
  • 自分も答えを知りたい
  • 質問者以外のユーザにも役立つ

グッドが多くついた質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

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

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

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

下記のような質問は推奨されていません。

  • 間違っている
  • 質問になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

適切な質問に修正を依頼しましょう。

PondVillege

2022/10/12 10:49

CSVデータの形式は当たっていますか?データのうち,左上5x5の範囲でいいので教えていただけると助かります.
jbpb0

2022/10/12 12:26 編集

> dt = df.T.iloc[0,1] # 時間刻み は、いくつ(何秒)でしょうか? > def highpass(x, samplerate, fp, fs, gpass, gstop): の関数内の、 > return y のすぐ上に print(y) を追加して実行して、結果表示を確認してみてください [nan nan nan ... nan nan nan] となってませんでしょうか?
NMKN

2022/10/12 12:44

ps_aux_grep様 ヘッダーも含めますと time. signal 0 0.00347 0.000005 -0.0118 0.00001. -0.0118 のようになっています。参照サイトにありましたサンプルCSVを参考に同様の形式にしております。
NMKN

2022/10/12 12:48

jbpb0様 おっしゃったことを実行したところ、jbpb0様のおっしゃったようにnanが続く結果となりました。この場合どのようなことが起きているのかわからないのですがご教授していただくことは可能でしょうか?
PondVillege

2022/10/12 14:41 編集

CSVデータに欠損がある場合にFFTを行うと,全てNaNになります.FFTとはいえ,実態は全ての値に対する畳み込みの演算です.1つでもNaNがあるなら,それを埋めなくてはなりません. ひとまず,NaNがあるかの確認をお願いします.具体的には,FFTしたいデータを取得data = df.T.iloc[i+1]した次の行で print(np.sum(np.isnan(data))) を記入して,NaNの個数を確認しましょう.0であれば他の原因を考えることになります.NaNがあれば,埋める方法を解答にて提示します. ちなみに,示していただいたCSVデータの一部ですが,かなり怪しいように見えます. CSVはComma Separated Valuesのはずですが,CommaではなくDotによるseparateや空白のみによるseparateが行われているように見えます.正しく掲示いただいてますでしょうか.
jbpb0

2022/10/12 14:37

元データにNaNが無くても、この条件ではNaNになると思います
NMKN

2022/10/12 14:43

ps_aux_grep様 提示されたコードを試したところ、0と表記されました。 jbpb0様 そもそものコード自体に問題があるという考え方ということでよろしいでしょうか?
NMKN

2022/10/12 14:54

ps_aux_grep様 こちらjupyter lab でprintしたものの一部をコピペしたものになります。 time signal 0 0.000000 0.00347 1 0.000005 -0.01180 2 0.000010 -0.01180 3 0.000015 -0.01180 4 0.000020 -0.01180
jbpb0

2022/10/13 00:36

https://stackoverflow.com/questions/50254625/why-it-returns-nan-while-using-scipy-signal-filtfilt と同じことが起きてるのだと思います > def highpass(x, samplerate, fp, fs, gpass, gstop): の関数内の、 > return y のすぐ上に print(wp) print(ws) print(N) を追加して実行して、結果表示を確認したら、フィルタの次数「N」が非常に大きくなってますよね 時間刻み「dt」で決まる「samplerate」に対する「fp」と「fs」の設定により、「wp」と「ws」の違いがわずかであり、非常に急峻なフィルタを設定してるためだと思います 「samplerate」を小さく(「dt」を大きく)するとか、「samplerate」に対する「fp」と「fs」の設定を見直す(「fp」と「fs」の差を大きくする)とかして、フィルタの次数「N」がもっと小さくなれば、「NaN」が出ずに計算できるようになると思います たとえば、下記の変更のどちらかを行ってみてください fp_hp = 520 ↓ 変更 fp_hp = 1000 または fs_hp = 500 ↓ 変更 fs_hp = 250 数値は適当なので、上記変更でうまくいかなければ、「fp」と「fs」の差がもっと大きくなるように変えてみてください そのあたりを変えたくないなら、下記の変更を試してみてください b, a = signal.butter(N, Wn, "high") y = signal.filtfilt(b, a, x) ↓ 変更 sos = signal.butter(N, Wn, "high", output='sos') y = signal.sosfiltfilt(sos, x) 参考 https://program4ptotat.livedoor.blog/archives/13572699.html
jbpb0

2022/10/13 01:49 編集

> 元データと見比べたところハイパスフィルタが適用されていないようです。 コードの最後に下記を追加したら、「df」と「df_filter」の周波数特性の比較ができます df_f = np.abs(np.fft.fft(df.values[:, 1], norm='ortho')) df_filter_f = np.abs(np.fft.fft(df_filter.values[:, 1], norm='ortho')) freq = np.fft.fftfreq(len(df_f), d=df.T.iloc[0,1]) fig, axs = plt.subplots(2, 1) axs[0].plot(freq, df_f, '.') axs[0].set_xlim(0, 1500) axs[1].plot(freq, df_filter_f, '.') axs[1].set_xlim(0, 1500) plt.show()
NMKN

2022/10/13 06:08

1T2R3M4様 別サイトでの質問、クローズさせていただきました。知らなかったのでありがとうございます。

回答2

0

ベストアンサー

python

1 b, a = signal.butter(N, Wn, "high") 2 y = signal.filtfilt(b, a, x)

↓ 変更

python

1 sos = signal.butter(N, Wn, "high", output='sos') 2 y = signal.sosfiltfilt(sos, x)

参考
A.2P Pythonによるデジタルフィルタ

 
【追記】
コードの最後に下記を追加したら、「df」と「df_filter」の周波数特性の比較ができます

python

1df_f = np.abs(np.fft.fft(df.values[:, 1], norm='ortho')) 2df_filter_f = np.abs(np.fft.fft(df_filter.values[:, 1], norm='ortho')) 3freq = np.fft.fftfreq(len(df_f), d=df.T.iloc[0,1]) 4 5fig, axs = plt.subplots(2, 1) 6axs[0].plot(freq, df_f, '.') 7axs[0].set_xlim(0, 1500) 8axs[1].plot(freq, df_filter_f, '.') 9axs[1].set_xlim(0, 1500) 10plt.show()

 
下記のようにして、適当にデータをでっち上げて確認してみました

python

1 df = pd.read_csv(in_file, encoding='SHIFT-JIS')

↓ 変更

python

1 t = np.arange(0, 5.5, 0.000005) 2 y = np.random.rand(len(t))*2-1 3 df = pd.DataFrame(np.vstack([t, y]).T)

「df」と「df_filter」の周波数特性の比較

投稿2022/10/13 00:47

編集2022/10/13 02:15
jbpb0

総合スコア7521

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

下記のような回答は推奨されていません。

  • 間違っている回答
  • 質問の回答になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

このような回答には修正を依頼しましょう。

回答へのコメント

jbpb0

2022/10/13 02:38 編集

上記の方法の場合は、 > N, Wn = signal.buttord(wp, ws, gpass, gstop) で計算した「N」と「Wn」をそのまま使って大丈夫だと思います そうすれば、質問のコードの「# ハイパスの設定」のところで設定してるフィルタ特性をそのまま実現できます
NMKN

2022/10/13 06:07

詳細なご回答ありがとうございます。jbpb0様のご提示して頂いたフィルタ特性を確認するコードで確認したところjbpb0様と同様のようなフィルタ特性が確認されました。ところでご相談なのですが私の質問文にグラフがあり、そちらを見ていただくと4.8sあたりで急峻な値があると思われます。フィルタ適用をしたところその値が小さくなっていました。このことはフィルタが適用されてると判断しても良いという考え方でよろしいのでしょうか? またご提示頂いた周波数特性の比較グラフについては横軸が周波数、縦軸がデータをFFTしたものと考えているのですがよろしいでしょうか?長くなってしまい申し訳ありません。
jbpb0

2022/10/13 07:38 編集

> 4.8sあたりで急峻な値があると思われます。フィルタ適用をしたところその値が小さくなっていました。このことはフィルタが適用されてると判断しても良いという考え方でよろしいのでしょうか? まずは、私の回答のようにノイズを入力にして、フィルター適用前後の周波数特性の比較をしてください ノイズは、いろんな周波数の情報を持ってるので、フィルターの効果が分かりやすいです ノイズ入力でフィルターの特性が把握できたら、他の入力(この質問のデータとか)を使ってみてください その入力データによっては、周波数特性を見ても効果が分かりにくいかもしれませんが、フィルターの効果はノイズの時と同じです > 周波数特性の比較グラフについては横軸が周波数、縦軸がデータをFFTしたものと考えているのですがよろしいでしょうか? そうです ただし、fft後のデータは(位相情報を含む)複素数で、縦軸はそれの絶対値(位相情報を消した振幅情報のみ)です
NMKN

2022/10/13 12:49

詳細なご返信ありがとうございます。フィルタは適用できたと考えられます。最後に私は初学者で頓珍漢なことをお聞きしたりなどご迷惑をかけたと思いますが丁寧なご回答ありがとうございました。そこで備忘録を学習用で記録しておきたいのですが、「今回のそもそもの問題は通過域周波数と阻止域周波数の違いがわずかであったため、フィルタの次数が大きくなってしまっていたのが原因である。またフィルタの適用に関しては適当なノイズで通してみて周波数特性を確認してフィルタの特性を把握することが望ましい。」以上のような内容で相違ないでしょうか?何度も同じようなことを聞くような形になってしまい申し訳ありません。
jbpb0

2022/10/14 08:25 編集

> 「今回のそもそもの問題は通過域周波数と阻止域周波数の違いがわずかであったため、フィルタの次数が大きくなってしまっていたのが原因である。 > またフィルタの適用に関しては適当なノイズで通してみて周波数特性を確認してフィルタの特性を把握することが望ましい。」 は合ってますが、不足してることがあります 上記の二つの文章の間に、下記を追加してください フィルタの次数が大きくて b, a = signal.butter(N, Wn, "high") y = signal.filtfilt(b, a, x) が正常に計算できない場合でも、 sos = signal.butter(N, Wn, "high", output='sos') y = signal.sosfiltfilt(sos, x) ならば計算できる場合がある 公式ページ https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html にも下記の記載があり、「sos」の使用を推奨してるので、問題無さそうでも「sos」を使う方がいいようだ 「Default is ‘ba’ for backwards compatibility, but ‘sos’ should be used for general-purpose filtering.」 「If the transfer function form [b, a] is requested, numerical problems can occur since the conversion between roots and the polynomial coefficients is a numerically sensitive operation, even for N >= 4. It is recommended to work with the SOS representation.」 「It’s recommended to use second-order sections format when filtering, to avoid numerical error with transfer function (ba) format」 (「second-order sections」は「sos」のことです)
NMKN

2022/10/15 06:42

なるほど、参考サイトも併せて確認しておきます。非常にためになりました。ありがとうございます!

0

N, Wn = signal.buttord(wp, ws, gpass, gstop)にて生成された,NWnに問題があったようです.次のようにして直しましょう.

Python

1def highpass(x, samplerate, fp): 2 b, a = signal.butter(5, fp, btype = 'high', fs = samplerate) 3 y = signal.filtfilt(b, a, x) 4 return y

確認のため,次のコードで波形を得てみました.

Python

1from scipy import signal, fftpack 2import numpy as np 3from matplotlib import pyplot as plt 4 5fs = 200000 6x = np.linspace(0, 4 * np.pi, fs * 3) 7data = np.sin(2 * np.pi * x) 8for i in range(500): 9 data += np.sin(2 * np.pi * x * np.random.randint(1, fs // 200) + np.random.randn() * np.pi) 10 11def highpass(x): 12 b, a = signal.butter(5, 500, btype = 'high', fs = fs) 13 w, h = signal.freqz(b, a, fs = fs) 14 return signal.filtfilt(b, a, x) 15 16def fft(x): 17 spec = fftpack.fft(x) 18 amp = abs(spec) 19 amp = amp / (len(x) / 2) 20 freq = np.linspace(0, fs, len(x)) 21 return amp[:len(x) // 2], freq[:len(x) // 2] 22 23fig, axs = plt.subplots(2, 2) 24axs[0, 0].plot(x, data) 25axs[0, 0].set_title('original wave') 26axs[0, 0].set_ylabel('amp') 27 28amp, freq = fft(data) 29axs[0, 1].plot(freq[:fs // 20], amp[:fs // 20]) 30axs[0, 1].set_title('original spectrum') 31axs[0, 1].set_ylabel('amp') 32 33hp = highpass(data) 34axs[1, 0].plot(x, hp) 35axs[1, 0].set_title('high passed wave') 36axs[1, 0].set_xlabel('time[s]') 37axs[1, 0].set_ylabel('amp') 38 39amp, freq = fft(hp) 40axs[1, 1].plot(freq[:fs // 20], amp[:fs // 20]) 41axs[1, 1].set_title('high passed spectrum') 42axs[1, 1].set_xlabel('frequency[Hz]') 43axs[1, 1].set_ylabel('amp') 44 45plt.show()

イメージ説明
阻止域端周波数の設定はNの大きさに依存しますが,これを大きくしすぎるとまたnanになるので,一応のこと応急処置的な回答になります.

投稿2022/10/12 16:53

編集2022/10/12 16:55
PondVillege

総合スコア1066

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

下記のような回答は推奨されていません。

  • 間違っている回答
  • 質問の回答になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

このような回答には修正を依頼しましょう。

回答へのコメント

NMKN

2022/10/12 21:46

ひとまず提案して頂いた修正コードで試したところ引数のエラーが出たため def highpass()の部分にfs gpass gstopを追加したところコードが作動し、下記のグラフに描画されない問題は解決したのですが、元データと見比べたところハイパスフィルタが適用されていないようです。修正箇所を下記に記しておきました。もしかしたら私の修正が間違っている可能性もあるためよろしければご教授いただければ幸いです。 """def highpass(x, samplerate, fp, fs, gpass, gstop): fn = samplerate / 2 #ナイキスト周波数 wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化 ws = fs / fn #ナイキスト周波数で阻止域端周波数を正規化 N, Wn = signal.buttord(wp, ws, gpass, gstop) #オーダーとバターワースの正規化周波数を計算 b, a = signal.butter(N, Wn, "high") #フィルタ伝達関数の分子と分母を計算 y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける print(y) return y""" def highpass(x, samplerate, fp, fs, gpass, gstop): b, a = signal.butter(5, fp, btype = 'high', fs = samplerate) y = signal.filtfilt(b, a, x) return y
PondVillege

2022/10/13 00:17

> 元データと見比べたところハイパスフィルタが適用されていないようです ちゃんと周波数空間上で比較しましたか? 私の例示した画像でも分かる通り,低周波がカットされたかどうかなんてパッと見でわからないと思います. 私のコードではdata = np.sin(2 * np.pi * x)となっていますが,そちらのdata = df.T.iloc[i+1]が使えるよう df = pd.read_csv(in_file, encoding='SHIFT-JIS') data = df.T.iloc[1] にでもしてデータを取り出して周波数空間上でプロットさせてみてはいかがでしょうか
NMKN

2022/10/13 05:57

ps_aux_grep様 流石に私の横着でした。生データ上で比較してもわかるはずがありませんね^^; 周波数空間上なりスペクトログラムを描くなりで確認したいと思います。

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

ただいまの回答率
86.02%

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

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

質問する

関連した質問

同じタグがついた質問を見る

JupyterLab

JupyterLabは、Jupyter notebookの後継の対話型開発環境(IDE)です。データの可視化がインタラクティブで、プラグイン作成により新しいコンポーネントの追加および既存のコンポーネントも統合可能。サーバに閉じているため、データ分析に向いています。

Python

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