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

Question

前提

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/ツールのバージョンなど）

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

Comments

[ps_aux_grep]

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

[jbpb0]

> 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]

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

[NMKN]

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

[ps_aux_grep]

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]

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

[NMKN]

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

[NMKN]

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

[1T2R3M4]

https://qiita.com/NMKN/questions/7f0850be05df4a9ebca7
以下ご対応ください。
https://teratail.com/help#posted-otherservice

[jbpb0]

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]

> 元データと見比べたところハイパスフィルタが適用されていないようです。

コードの最後に下記を追加したら、「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]

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

Answer 1

    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)

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

　
【追記】
コードの最後に下記を追加したら、「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()

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

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

↓ 変更

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

Comments

[jbpb0]

上記の方法の場合は、

> N, Wn = signal.buttord(wp, ws, gpass, gstop)

で計算した「N」と「Wn」をそのまま使って大丈夫だと思います
そうすれば、質問のコードの「# ハイパスの設定」のところで設定してるフィルタ特性をそのまま実現できます

[NMKN]

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

[jbpb0]

> 4.8sあたりで急峻な値があると思われます。フィルタ適用をしたところその値が小さくなっていました。このことはフィルタが適用されてると判断しても良いという考え方でよろしいのでしょうか？

まずは、私の回答のようにノイズを入力にして、フィルター適用前後の周波数特性の比較をしてください
ノイズは、いろんな周波数の情報を持ってるので、フィルターの効果が分かりやすいです

ノイズ入力でフィルターの特性が把握できたら、他の入力(この質問のデータとか)を使ってみてください
その入力データによっては、周波数特性を見ても効果が分かりにくいかもしれませんが、フィルターの効果はノイズの時と同じです

> 周波数特性の比較グラフについては横軸が周波数、縦軸がデータをFFTしたものと考えているのですがよろしいでしょうか？

そうです
ただし、fft後のデータは(位相情報を含む)複素数で、縦軸はそれの絶対値(位相情報を消した振幅情報のみ)です

[NMKN]

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

[jbpb0]

> 「今回のそもそもの問題は通過域周波数と阻止域周波数の違いがわずかであったため、フィルタの次数が大きくなってしまっていたのが原因である。

> またフィルタの適用に関しては適当なノイズで通してみて周波数特性を確認してフィルタの特性を把握することが望ましい。」

は合ってますが、不足してることがあります
上記の二つの文章の間に、下記を追加してください

フィルタの次数が大きくて
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]

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

Answer 2

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

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

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

from scipy import signal, fftpack
import numpy as np
from matplotlib import pyplot as plt

fs = 200000
x = np.linspace(0, 4 * np.pi, fs * 3)
data = np.sin(2 * np.pi * x)
for i in range(500):
    data += np.sin(2 * np.pi * x * np.random.randint(1, fs // 200) + np.random.randn() * np.pi)

def highpass(x):
    b, a = signal.butter(5, 500, btype = 'high', fs = fs)
    w, h = signal.freqz(b, a, fs = fs)
    return signal.filtfilt(b, a, x)

def fft(x):
    spec = fftpack.fft(x)
    amp = abs(spec)
    amp = amp / (len(x) / 2)
    freq = np.linspace(0, fs, len(x))
    return amp[:len(x) // 2], freq[:len(x) // 2]

fig, axs = plt.subplots(2, 2)
axs[0, 0].plot(x, data)
axs[0, 0].set_title('original wave')
axs[0, 0].set_ylabel('amp')

amp, freq = fft(data)
axs[0, 1].plot(freq[:fs // 20], amp[:fs // 20])
axs[0, 1].set_title('original spectrum')
axs[0, 1].set_ylabel('amp')

hp = highpass(data)
axs[1, 0].plot(x, hp)
axs[1, 0].set_title('high passed wave')
axs[1, 0].set_xlabel('time[s]')
axs[1, 0].set_ylabel('amp')

amp, freq = fft(hp)
axs[1, 1].plot(freq[:fs // 20], amp[:fs // 20])
axs[1, 1].set_title('high passed spectrum')
axs[1, 1].set_xlabel('frequency[Hz]')
axs[1, 1].set_ylabel('amp')

plt.show()

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

Comments

[NMKN]

ひとまず提案して頂いた修正コードで試したところ引数のエラーが出たため 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

[ps_aux_grep]

> 元データと見比べたところハイパスフィルタが適用されていないようです
ちゃんと周波数空間上で比較しましたか？
私の例示した画像でも分かる通り，低周波がカットされたかどうかなんてパッと見でわからないと思います．
私のコードでは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]

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